Zhao et al.
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
2631
Two-dimensional nonseparable linear canonical transform: sampling theorem and unitary discretization Liang Zhao,1 John J. Healy,2 and John T. Sheridan1,* 1
School of Electrical, Electronic and Communications Engineering, Communications and Optoelectronic Research Centre, SFI–Strategic Research Cluster in Solar Energy Conversion, College of Engineering and Architecture, University College Dublin, Belfield, Dublin 4, Ireland 2 Department of Electronic Engineering, Maynooth University, Co. Kildare, Ireland *Corresponding author: [email protected] Received August 4, 2014; revised October 5, 2014; accepted October 6, 2014; posted October 6, 2014 (Doc. ID 214789); published November 11, 2014 The two-dimensional (2D) nonseparable linear canonical transform (NS-LCT) is a unitary, linear integral transform that relates the input and output monochromatic, paraxial scalar wave fields of optical systems characterized by a 4 × 4 ray tracing matrix. In addition to the obvious generalizations of the 1D LCT (which are referred to as separable), the 2D-NS-LCT can represent a variety of nonaxially symmetric optical systems including the gyrator transform and image rotation. Unlike the 1D LCT, the numerical approximation of the 2D-NS-LCT has not yet received extensive attention in the literature. In this paper, (1) we develop a sampling theorem for the general 2D-NS-LCT which generalizes previously published sampling theorems for the 1D case and (2) we determine which sampling rates may be chosen to ensure that the obvious discrete transform is unitary. © 2014 Optical Society of America OCIS codes: (000.4430) Numerical approximation and analysis; (070.2025) Discrete optical signal processing; (070.2590) ABCD transforms; (090.1995) Digital holography. http://dx.doi.org/10.1364/JOSAA.31.002631
1. INTRODUCTION The one-dimensional (1D) linear canonical transform (LCT) [1–5] is a three-parameter group of linear integral transforms, which can be used to model the propagation of the coherent wave fields through first-order optical systems. Among its special cases [6,7] are the Fourier transform (FT), the fractional Fourier transform (FRT), the Fresnel transform (FST), the scaling transform, and the chirp transform. The two-dimensional separable linear canonical transform (2D-S-LCT) can model a wide variety of isotropic optical systems [8], which can be treated as being made up of two independent 1D-LCTs. In the case when the system is fully symmetric and orthogonal, and the parameters for both of the two dimensions are identical, such systems can be represented using only three parameters as in the 1D-LCT case [9]. In the case when the system is still orthogonal but the parameters for the orthogonal dimensions differ, such systems can be represented by six free parameters [9,10]. The two-dimensional nonseparable linear canonical transform (2D-NS-LCT) is a more general integral transform, in which the two dimensions are coupled to each other by four additional cross-parameters, increasing the total number of free parameters to ten [8]. This general nonseparable transform is nonorthogonal, nonaxially symmetric, and anamorphic [8]. The 2D-NS-LCT is able to represent systems not only involving coupling and shearing between components along the different dimensions [9,11,12] but also involving rotations between any arbitrary planes in phase space, such as gyrations [13–17]. The 2D-NS-LCT is important in digital signal processing and optical system modeling [18–26]. 1084-7529/14/122631-11$15.00/0
Discrete transforms that approximate the continuous LCTs allow us to simulate a number of optical propagation problems. Digital simulations to solve these problems require access to a numerical algorithm to approximate the continuous transform efficiently, accurately, and unitarily [27,28]. In the area of optics not only is such an algorithm important when employing iterative techniques to design systems, but it is also significant when solving inverse scatter problems to extract information about an object given a known illumination [27–32], for instance in phase retrieval [33,34]. Traditionally, such problems have been addressed in situations where the FT and FST were sufficient to describe the optical system examined. Unitary algorithms for both the FT [35–40] and FST [41] are widely available. In order to further avoid aliasing during sampling processes, Ding [42] developed a sampling theorem for the general 1D-LCT in 2001 (also see [7]). In 2005, Hennelly and Sheridan [6] discussed how the space–bandwidth product of a signal (and thus the sampling requirements) can be tracked through an LCT system. However the resulting discrete transforms are not in general unitary. In [43–52], these results were shown to lead to a fast and accurate numerical implementation of the 1D-LCT. Recently a sufficient condition on the sampling rates chosen during discretization to ensure that the discrete 1D-LCT (1D-DLCT) is unitary have been presented [52,53]. Proof that these conditions guarantee unitarity is provided in [53]. It should be noted that the sampling theorem developed in [53] is consistent with the Shannon sampling theorem for the FT, the Gori sampling theorem [41] for the FST, and the Ding sampling theorem [42] for the LCT. © 2014 Optical Society of America
2632
J. Opt. Soc. Am. A / Vol. 31, No. 12 / December 2014
Zhao et al.
Numerical calculation of the 2D-NS-LCT has recently received attention [8,17,22–24,54–57]. This paper will focus on developing a sampling theorem for this case and demonstrating how a unitary sampling rate can be chosen. This paper is organized as follows: In Section 2 we introduce the transforms examined for 1D and 2D signals and provide some other important definitions. One powerful technique used to analyze LCTs is to decompose them into, and then apply sequentially, a set of simpler transforms. In Section 3, we introduce this kind of decomposition, discuss some important special cases, and look at the insight one such decomposition gives us about 2D-NS-LCT’s application to discrete signals. In Section 4, the general sampling theorem for the 2DNS-LCT is developed. In Section 5, we consider one DLCT, a discrete transformation that approximates the continuous integral transform. In particular the parameters for which the DLCT is unitary are determined. Finally, a brief conclusion and future work are given in Section 6.
2. DEFINITIONS OF THE CONTINUOUS LCTs A. Definitions of the 1D-LCT and Its Unitary Property The continuous 1D-LCT of a signal gx can be defined by [5,6,53] Gx0 LM fgxgx0 h i ( R 2 0 Dx02 p1 ∞ exp jπAx −2xx gxdx; −∞ B pjB 02 0 D expjπCDx gDx ;
where k1 d11 b22 − d12 b21 ;
k2 2−d11 b12 d12 b11 ;
k3 −d21 b12 d22 b11 ; p2 2a12 b22 − a22 b12 ;
p1 a11 b22 − a21 b12 ; p3 −a12 b21 a22 b11 :
In this case the transformation matrix M of the system is defined as 0
a11 B a21 B M @ c11 c21
a12 a22 c12 c22
b11 b21 d11 d21
1 b12 b22 C C A d12 A C d22
B ; D
B0
; (1)
where given M f A B ; C D g, the parameters of the transform matrix M are restricted to be real numbers with detM 1, and detM represents the determinant of matrix M. For the case of B 0 the LCT is defined as the limit of the B ≠ 0 case with B → 0 [5]. The LCT with parameters fA B; C Dg f0 1; −1 0g corresponds to the FT, while fcos θ sin θ; −sin θ cos θg corresponds to the FRT. The FST results for f 1 z ; 0 1 g. In addition, the transform matrix of the inverse 1D-LCT is given by M −1 f D −B ; −C A g.
(3-3)
where A, B, C, and D are now 2 × 2 submatrices, and detB ≠ 0. We note that the 2D-NS-LCT has 16 parameters, see Eq. (3-3), which must satisfy the constraints ABT BAT , CDT DC T , and ADT − BC T I, or AT C C T A, BT D DT B, and AT D − C T B I, where I is a 2 × 2 identity matrix. Since the above equations are equivalent to six constraints, the number of independent parameters in the matrix M is in fact reduced to ten [8,54–57]. In addition, the inverse 2D-NS-LCT of Gx0 ; y0 is given by gx; y LM −1 fGx0 ; y0 gx; y;
B≠0
(3-2)
where M −1 f DT
−BT ; −C T
(4)
AT g.
2.2 Other Definitions It should be noted that ① When B 0; ⇒ detB 0, the definition of the 2D-NSLCT is [9,24,55,56] Gx0 ; y0 LM fgx; ygx0 ; y0 p detD expjπc11 d11 c12 d12 x02 × expj2πc11 d21 c12 d22 x0 y0
2.1 Unitary Property One of the most important properties of the continuous LCT (both 1D and 2D cases) is the unitary property [3],
× expjπc21 d21 c22 d22 y02 gd11 x0 d21 y0 ; d12 x0 d22 y0 ;
LM −1 fLM fggg g:
where M f DT 0 ; C D g. ② When detB 0 but B ≠ 0, the 2D-NS-LCT is more complicated; see [24].
(2)
B. Definitions of the 2D-NS-LCT The continuous 2D-NS-LCT of a signal gx; y can be defined by [8,54–57] Gx0 ; y0 LM fgx; ygx0 ; y0 ∞ ZZ 1 jπk1 x02 k2 x0 y0 k3 y02 p exp detB j detB −∞ j2π−b22 x0 b12 y0 x b21 x0 − b11 y0 y × exp detB jπp1 x2 p2 xy p3 y2 gx; ydxdy; (3-1) × exp detB
(5) −1
C. Special Cases of the 2D-LCT 1. Special Cases of the 2D-S-LCT In the case when a12 a21 b12 b21 c12 c21 d12 d21 0;
(6)
the 2D-NS-LCT, see Eqs. (3-1) and (5), reduces to the separable transform, the well-known special cases of which are the 2D FT, the 2D FST [21], the 2D FRT [21,22], the 2D chirp transform [18,19,24], and the 2D scaling transform [24]. The transform matrix M of these special cases is shown in Table 1.
Zhao et al.
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
2633
Table 1. Transform Matrices of Some Well-known Cases of the 2D-S-LCT and 2D-NS-LCT 2D-Separable
0
0 B 0 B @ −1 0 0 1 B0 B @0 0
FT
FST
0
cos θx B 0 B @ − sin θx 0 0
FRT
Chirp
Scaling-x
Scaling-y
0 0 0 −1 0 1 0 0
1 0 0 0 z 0 1 0
1 0 1C C 0A 0 1
0 zC C 0A 1
0 0 1 0
Gyrator
Coupling
0 sin θx cos θy 0 0 cos θx − sin θy 0
1 0 B 0 1 B @ −1∕f x 0 0 −1∕f y 0 1∕sx 0 0 B 0 1 0 B @ 0 0 sx 0 0 0 0 1 0 0 B 0 1∕sy 0 B @0 0 1 0 0 0
1 0 sin θy C C A 0 cos θy 1
Coordinate (general)
0 0C C 0A 1 1
Rotation
0 0C C 0A 1 1 0 0C C 0A sy
Shearing-x
Shearing-y
2. Special Cases of the 2D-NS-LCT Three well-known 2D nonseparable transforms are the gyrator transform [13–17,21], the coupling (i.e., multiplication) transform [24], and the homogenous coordinate/affine transform [24–26], which can be defined by LM gyrator fgx;yg 1 jsin θj
Z∞ Z exp −∞
2D-Nonseparable 1 0 cos θ 0 0 sin θ B 0 cos θ sin θ 0 C C B @ 0 − sin θ cos θ 0 A − sin θ 0 0 cos θ 1 0 1 0 0 0 B0 1 0 0C C B @0 τ 1 0A τ 0 0 1 1 0 0 0 a11 a12 B a21 a22 0 0 C B a22 −a21 C A @ 0 0 detA detA −a12 a11 0 0 detA detA 1 0 cos θ sin θ 0 0 B − sin θ cos θ 0 0 C C B @ 0 0 cos θ sin θ A 0 0 − sin θ cos θ 1 0 0 0 1 shx B0 1 0 0C C B @0 0 1 0A 0 0 −shx 1 1 0 1 0 0 0 B shy 1 0 0 C C B @ 0 0 1 −shy A 0 0 0 1
j2πxy x0 y0 cos θ − xy0 yx0 gx;ydxdy; sin θ (7-1)
3. MATRIX DECOMPOSITION AND CONSEQUENCES FOR THE LOCATIONS OF THE SAMPLES IN THE OUTPUT PLANE The ABCD matrices of concatenated systems may be obtained by multiplication; for example, if a wave field propagates through a system with ABCD matrix M 1 and then a second with matrix M 2 , the overall system’s ABCD matrix is given by M 3 M 1 M 1 . We can use this result in two ways: First, we can find the output of a series of systems in a single calculation. Second, we can decompose a calculation into a sequence of simpler operations. A simple example is the venerable spectral method for calculating the FST,
LM Coupling xy fgx; yg expjπτxygx; y;
LM Coupling x0 y0 fGx0 ; y0 g expjπτx0 y0 Gx0 ; y0 ;
(7-2-1)
(7-2-2)
a x0 − a12 y0 −a21 x0 a11 y0 ; : (7-3) LM coordinate fgx; yg g 22 detA detA The transform matrices of these three cases are listed in Table 1. For brevity, the homogenous coordinate transform is referred to as the “coordinate transform” below. The coordinate transform includes some well-known transforms, like rotation [8,20,21,24] and shearing (interferometer) [24]; see Table 1.
1 λz 0 1
0 −1 1 0
1 0 −λz 1
0 −1
1 : 0
(8)
Equation (8) shows that the FST may be calculated by taking a FT, a chirp multiplication, and an inverse FT, all of which may be calculated using simple or mature software packages. Two decompositions have been proposed in the literature to numerically calculate the 2D-NS-LCT, the Iwasawa decomposition [8] and a second decomposition which may be thought of as a generalization of the direct method for calculating the 1D case [55]. This second decomposition is given by M M 4M 3M 2M 1 B 0 0 I 0 −1 T −1 0 B −I DB I where
I 0
I
−1
B A
0 ; I
(9-1)
2634
J. Opt. Soc. Am. A / Vol. 31, No. 12 / December 2014
Zhao et al.
0
1 k1 k2 B detB 2 detB C C; DB−1 B @ A k2 k3 2 detB detB 0 1 p1 p2 B detB 2 detB C C; B−1 A B @ A p2 p3 2 detB detB 0 b −b21 1 22 B detB detB C −1 C; BT B @ −b A b 12
detB
and 1(e). The relevant sampling relationships are discussed in Section 4. The decomposition in Eq. (9) is valid only if detB 0. No general result exists at this time for the detB 0 case. However, if we further specify that B 0, we can decompose the matrix as follows: M M 2M 1 (9-2)
11
detB
and where k1 ∼ k3 and p1 ∼ p3 are as defined in Eq. (3-2). M 1 and M 4 in Eq. (9-1) produce phase terms in the space and LCT domains, respectively; M 3 generates a coordinate transform; and M 2 is an FT. The matrix decompositions are discussed in detail in Appendix A. This decomposition would appear to be an excellent candidate for numerical approximation of the 2D-NS-LCT. It consists of two generalized chirps (M 1 and M 4 ), a 2D FT (M 2 ) and an image coordinate transformation (M 3 ). This last operation indicates an interesting and vital feature of such numerical approximations. Given an input comprised of a rectangular grid of regular samples, we normally expect a discrete transform to relate this input with a rectangular grid of output samples. However, in the discrete 2D-NS-LCT case, it will in general result in a skewed output grid. This is illustrated in Fig. 1. This has significant implications for an appropriate sampling theorem. Consider the 1D case. Given a figure’s width in the input and output domains, it has been shown that the appropriate sampling rates are related to these widths. In the separable 2D case, a similar result is obtained given the input and output widths in both directions; see Figs. 1(b) and 1(d). However, for the nonseparable case, the output samples obtained form a parallelogram of points; see Figs. 1(c)
I CDT
0 I
DT 0
−1
0 ; D
(10-1)
d22 −d21 −d12 d11 where DT−1 f detD detD ; detD detD g. Here M 2 produces a phase term in the LCT domain, while M 1 generates the following coordinate transform,
LM 1 fgx; ygx0 ; y0 gd11 x0 d21 y0 ; d12 x0 d22 y0 : (10-2) Consequently, beginning with a rectangular grid of samples, once again one ends up with a parallelogram arrangement of output samples.
4. SAMPLING THEOREM FOR THE 2D-NS-LCT A. Sampling Effects We consider the impulse sampled version of gx; y, gˆ x; y gx; y
∞ X
∞ ∞ X X
δx − nT x ; y − mT y
n−∞ m−∞ ∞ X
gnT x ; mT y δx − nT x ; y − mT y
n−∞ m−∞
∞ ∞ X X 1 j2πnx j2πmy exp ; gx; y exp T xT y Tx Ty n−∞ m−∞ (15) where T x and T y are the sampling intervals in the x and y directions, respectively. 1. The Effect of Sampling in the detB ≠ 0 Case From Eq. (3), the 2D-NS-LCT of gˆ x; y can be defined by ˆ 0 ; y0 LM fˆgx;ygx0 ; y0 Gx ∞ ∞ X X 1 p H ×G T x T y j detB n−∞ m−∞ nb11 mb12 0 nb21 mb22 0 × x− − ;y − − Tx Ty Tx Ty ∞ ∞ X X 1 p H × Gx0 − xˆ 0n;m ; y0 − yˆ 0n;m ; T x T y j detB n−∞ m−∞ (12-1)
Fig. 1. Input image shape, output image shapes, and their sampling grids: (a) rectangular shape of an input image, with length Lx and height Ly . (b) and (c) are the output image shapes of the 2D S-LCT and NS-LCT, respectively; the output image shape in the separable case is still a rectangle, while the output image shape in the nonseparable case is skewed, that is, becomes a parallelogram (the side lengths of which are different from the values in the separable case). (d) and (e) are the sampling grids of (b) and (c), respectively.
where xˆ 0n;m and
nb11 mb12 ; Tx Ty
and yˆ 0n;m
nb21 mb22 ; Tx Ty
(12-2)
Zhao et al.
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
−jπk1 x0 − xˆ 0 2 k2 x0 − xˆ 0 y0 − yˆ 0 k3 y0 − yˆ 0 2 detB 02 0 0 jπk1 x k2 x y k3 y02 × exp : (12-3) detB
H exp
Equation (12-1) indicates that the 2D-NS-LCT of gˆ x; y, ˆ 0 ; y0 , is made up of modulated shifted replicas namely Gx of the original signal Gx0 ; y0 [7,53,57]. The 2D-NS-LCT of a sampled signal is therefore a periodically replicated version of the 2D-NS-LCT of the original continuous signal before sampling. From Eq. (12-2), the center point (coordinates) of the n; mth replica are designated by xˆ 0 ; yˆ 0 n;m , the location of which is dependent on T x ; T y and on the values of all four parameters of the 2 × 2 matrix B in the 2D-NS-LCT. It should be noted that, in the separable case, the xˆ 0 ; yˆ 0 n;m points are located on a Cartesian coordinate system; see the dashed lines in Fig. 2(a). In the nonseparable case, the xˆ 0 ; yˆ 0 n;m points are located along a series of parallel lines; see Fig. 2(b). 2. Effect of Sampling in the B 0 Case From Eq. (8) [24], as in the former detB ≠ 0 case, the 2D-NS-LCT of gˆ x; y can be defined by
B. Sampling Theorem The sampling rate should be such that the replicas indicated by Eq. (12-1) do not overlap; then the continuous signal may be recovered from the samples taken at this rate. 1. Sampling Theorem for the detB ≠ 0 Case Using Eq. (12), and in order to avoid aliasing, the output LCT image extents should satisfy Lx0 ≤
(13) It should be noted that ① In the case when detB ≠ 0, the 2D-NS-LCT of the ˆ 0 ; y0 , is still a continusampled input signal gˆ x; y, namely Gx ous signal; see Fig. 2b. ② In the case when B 0, the 2D-NS-LCT of the sampled input signal gˆ x; y is a discrete signal. The output S-LCT samples and NS-LCT samples are arranged as shown in Figs. 1(d) and 1(e).
1 Tx
q b211 b221 ;
and Ly0 ≤
1 Ty
q b212 b222 ;
(14-1)
q b211 b221 is the distance between the n-direction q replicas, while T1y b212 b222 is the distance between the where
1 Tx
m-direction replicas [see Fig. 2(b)], in which θx0 arctanbb21 11 and θy0 arctanbb22 . The derivations of these two distances 12 are discussed in Appendix B. 2. Sampling Theorem for the B 0 Case When detB 0, the effect of the LCT on the shape change is discussed in Appendix A. The sampling theorem for this case is
ˆ 0 ; y0 Gx
∞ ∞ X X nd md12 nd md22 δ y0 − 21 − : Gx0 ; y0 δ x0 − 11 − Tx Ty Tx Ty n−∞ m−∞
2635
Lx0 ≤
q Lx d212 d222 detD
;
and Ly0 ≤
q Ly d211 d221 detD
;
(14-2)
where θx0 − arctandd12 and θy0 − arctandd11 ; see Fig. 2(b). 22 21 In conclusion, as long as Eq. (14) is satisfied, Gx0 ; y0 can be reconstructed from the periodically replicated version of ˆ 0 ; y0 . Gx0 ; y0 , namely Gx C. Sampling Theorem for Special Cases The sampling theorems for some well-known special cases are shown in Table 2. The well-known Nyqist sampling theorem for the FT, the Gori sampling theorem [41] for the FST, and the Ding sampling theorem [42] for the 1D-LCT are special cases of the theorem developed here.
5. UNITARY DISCRETE TRANSFORM A. Unitary Discretization During discretization, the continuous coordinates x; y are replaced by the n; mth sample locations x; yn;m , that is, the coordinates nT x ; mT y , while x0 ; y0 are replaced by the n0 ; m0 th sample x0 ; y0 n0 ;m0 , where the integers n, m, n0 , m0 ∈ −N∕2 1; N∕2. N is the number of sampling points, which is assumed to here to be even valued. When N is odd, the analysis is similar and does not warrant separate attention. Fig. 2. (a) and (b) are the outputs of a 2D-S-LCT and a 2D-NS-LCT, respectively, given a discrete input. Were the inputs continuous, we would obtain outputs at locations indicated by the shaded (blue where color is available) tetrahedrons, while the other tetrahedrons indicate the locations of replicas that appear because of the discrete input. Were we to sample the outputs, we should do so on grids corresponding to Figs. 1(d) and 1(e) located in the shaded regions.
1. Unitary Discretization for the detB ≠ 0 Case The corresponding discrete transform is unitary, if the n0 ; m0 th sample x0 ; y0 n0 ;m0 is defined by x0n0 ;m0
n0 b11 m0 b12 ; Lx Ly
and y0n0 ;m0
n0 b21 m0 b22 ; (15) Lx Ly
2636
J. Opt. Soc. Am. A / Vol. 31, No. 12 / December 2014
Zhao et al.
Table 2. Sampling Theorem and Unitary Sampling Rates for Some Well-known Cases of the 2D-S-LCT and 2D-NS-LCT
S-LCT when B ≠ 0
S-LCT when B 0
LCTs
Sampling Theorems
Unitary Rate Chosen
2D-FT
Lx0 ≤ 1∕T x , and Ly0 ≤ 1∕T y
T x0 1∕Lx , and T y0 1∕Ly
2D-FRT
Lx0 ≤ sin θx ∕T x , and Ly0 ≤ sin θy ∕T y
T x0 sin θx ∕Lx , and T y0 sin θy ∕Ly
2D-FST
Lx0 ≤ λz∕T x , and Ly0 ≤ λz∕T y
T x0 λz∕Lx , and T y0 λz∕Ly
2D-S-LCT
Lx0 ≤ b11 ∕T x , and Ly0 ≤ b22 ∕T y
T x0 b11 ∕Lx , and T y0 b22 ∕Ly
2D Chirp
Lx0 ≤ Lx , and Ly0 ≤ Ly
T x0 T x , and T y0 T y
Scaling-x
L ≤ Lx ∕S x , and L ≤ Ly
T x0 T x ∕S x , and T y0 T y
x0
y0
Scaling-y
Lx0 ≤ Lx , and Ly0 ≤ Ly ∕S y
T x0 T x , and T y0 T y ∕S y
NS-LCT when B ≠ 0
Gyrator
NS-LCT when B 0
Coordinate
Lx0 ≤ sin θ∕T y , and Ly0 ≤ sin θ∕T x q q Lx0 ≤ Lx a211 a221 , and Ly0 ≤ Ly a212 a222
T x0 sin θ∕Ly , and T y0 sin θ∕Lx q q T x0 T x a211 a221 , and T y0 T y a212 a222
Lx0 ≤ Lx , and Ly0 ≤ Ly p Lx0 ≤ Lx , and Ly0 ≤ Ly 1 sh2x q Lx0 ≤ Lx 1 sh2y , and Ly0 ≤ Ly
T x0 T x , and T y0 T y p T x0 T x , and T y0 T y 1 sh2x q T x0 T x 1 sh2y , and T y0 T y
Coupling Shearing-x Shearing-y
2. Unitary Discretization for the B 0 Case The corresponding discrete transform is unitary if
where Lx and Ly are the extents (i.e., length and height) of the input image; see Fig. 1(a). The input and output samples each form N × N vectors. The sampling intervals in the space domain are T x Lx ∕N and T y Ly ∕N. The kernel of the discrete transform is a 4D matrix [58–60], namely an N × N × N × N matrix, which is denoted as W M . The final discrete transform is
y0n;m
Gx0n0 ;m0 ; y0n0 ;m0
N∕2 X
N∕2 X
W M n0 ; m0 ; n; mgnT x ; mT y :
nd22 T x − md21 T y ; and detD −nd12 T x md11 T y : detD
x0n;m
(17)
B. Unitary Properties The continuous 2D-NS-LCT is unitary, but discretization can destroy this property. Mathematically, a square matrix W M is unitary if
(16-1)
n−N∕21 m−N∕21
WH M W M I; The operation in Eq. (16-1) corresponds to a matrix multiplication. From Eq. (15), the elements of the kernel W M are given by
(18-1)
where W H M represents the conjugate transpose of W M , and I is an identity matrix, that is,
8 h
2 0
0
0
2 i9 > > 0b 0b 0b 0b 0b n m n b m n b m n b m > > 11 12 11 12 21 22 21 22 = j detB n−N∕21 m−N∕21 > :
× exp
x
y
x
y
x
y
x
y
> > ;
detB
8 nh h 0 0
0
i
0
i o9 0b 0b 0b 0b > > n b m n b m n b m n b m > 11 12 21 22 11 12 21 22 = x
y
x
y
> > :
x
y
x
jπp1 nLx 2 p2 nLx mLy p3 mLy 2 : × exp N 2 detB
It should be noted that in the separable case, the sampling points in the LCT domain are located in a rectangular grid; see Fig. 1(d). In the nonseparable case, the sampling points are located in an irregular grid, namely a skewed grid; see Fig. 1(e).
y
> > ;
detB
(16-2)
WH M W M −1 :
(18-2)
① For the 1D-DLCT, this identity matrix is 2D. ② In the general 2D-NS-DLCT, the matrix is 4D [58–60]. It should be noted that since the two dimensions in the
Zhao et al.
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
2637
2D-S-DLCT can be treated independently, by means of the row–column algorithm [53], the identity matrix in separable cases can be reduced to a 2D matrix [53]. Significantly, unlike the corresponding continuous LCTs in Eq. (2), the DLCTs are not necessarily unitary; that is, without due care discretization can destroy this property. In fact Eq. (16-2) must be multiplied by a scaling factor in order to make it unitary; see Appendix D. First, in order to show how the sampling grid can be skewed by the 2D-NS-LCT, one example is given in Fig. 3. In this example, the sample locations are calculated using Eq. (15), which guarantees that the corresponding discrete transform is unitary. As we discussed in Section 3, the transform can be decomposed in a number of ways. Of particular interest for this discussion is Eq. (9), which shows that we can decompose any 2D LCT where B−1 exists into several elements which do not alter the sampling grid and one coordinate transform. In the example given in Fig. 3, we have chosen LCT parameters such that the coordinate transform is simply a rotation. Consequently, the sampling grid in the LCT domain is rotated, specifically by 60° clockwise; see Fig. 3(b) (where θx0 −60° and θy0 30°). In Section 4 of this paper, we have shown that the sampling rates that result in a unitary discrete transform are integer multiples of the rate given by Ding’s sampling theorem for the LCT. One might reasonably wonder if the gains made by using this transform are simply the result of using an appropriate sampling rate, if indeed unitarity is of real importance. This is not the case. We merely constrain the relationship between the input and output sampling rates; this is signal-independent. In order to demonstrate this, we have generated Fig. 4. This figure shows the error that arises when using the same input signal and LCT used to generate Fig. 3. This signal is sampled and then forward transformed. Finally the resulting output is inverse transformed and compared to the original input. Were the improvement merely a matter of a higher sampling rate reducing aliasing, we would expect to see a monotonic decrease in error with decreasing sampling period. Instead, we observe a “sweet spot” where the transform is unitary. We vary the sampling period of x0 from
Fig. 4. If we propagate to a given LCT domain and then back to the input domain, how does the error vary with sampling period? In this plot, we show the MSE is minimized, and in fact is 0, only at the sampling rate we propose, see Eq. (15) , i.e., the rate at which the DLCT is unitary.
95% to 105% of the rate we propose, which gives unitarity. The mean square error (MSE) between the transformed signal, g¯ x; y, and the original signal, gx; y, is calculated by MSE
(19)
where g¯ is a N × N array of complex values. It can be seen that (1) the minimum MSE value is obtained by choosing the unitary sampling rate defined in Eq. (15); (2) as the difference between the unitary and nonunitary sampling rates increases, the corresponding MSE increases rapidly. C. Special Cases of the Unitary 2D-NS-DLCT 1. Some Special Cases of the Unitary S-LCT For the separable cases when detB ≠ 0, Eq. (15) can be simplified as x0n0 ;m0
n0 b11 ; Lx
and y0n0 ;m0
m0 b22 : Ly
(20)
After dividing each side by N, Eq. (20) becomes T x0
Fig. 3. One example of unitary calculation using the DLCT: (a) the input signal is a 2D Gaussian, 20 × 20 pixels, variances 0.004 and 0.002. The extents (Lx and Ly ) are both 0.05. (b) The DLCT for parameters [0.08 0.09 0.5 0.866; 0.1 −0.0858 −0.8660.5 0.5; −0.495 −0.858 0.05 0.06; 0.856 −0.494 0.07 −0.0558] was calculated using the formula in Eq. (16) and plotted using the contour function in MATLAB. If we decompose the matrix of parameters as per Eq. (9), we find the coordinate transformation matrix is a rotation matrix, which is accounted for in the plotting of the results. According to Subsection 4.B, here arctan−0.866 and θy0 arctanbb22 θx0 arctanbb21 0.5 −60°, 11 12 0.5 30°, that is, the rectangular sampling grid in the LCT arctan0.866 domain is rotated 60° clockwise.
PP jgx; y − g¯ x; yj2 PP ; jgx; yj2
b11 ; Lx
and T y0
b22 : Ly
(21)
Equation (21) is in fact Ding sampling [42]. For the 2D-S-LCT case, we have offered a sufficient condition to ensure it is unitary [53]. Unitary sampling rates for some special cases of the 2D-S-LCT are shown in Table 2. In order to obtain the kernal in the 2D-S-DLCT case when detB ≠ 0, we only need to set b12 b21 0 in Eq. (16-2). However, this approach does not provide a fast way to evaluate the continuous transform. The row–column algorithm provides an efficient way to numerically approximate the 2D-S-LCT [53], in which the 2D transform is reduced to two 1D-LCTs. The corresponding sampling points in the LCT domain are located in a rectangular grid; see Fig. 1(d). 2. Some Special Cases of the Unitary NS-LCT The unitary sampling rates for some special cases are shown in Table 2.
2638
J. Opt. Soc. Am. A / Vol. 31, No. 12 / December 2014
Zhao et al.
0
6. CONCLUSION AND FUTURE WORK The 2D case of the LCT has not received as much attention as the 1D case, and in most 2D cases examined it has been assumed that the optical system is separable, thus permitting 1D analysis. In this paper the general 2D-NS-LCT is explored, namely those cases of the 2D-LCT explicitly not reducible to 1D-LCT analysis. A sampling theorem for this transform has been derived, and it has been shown that the 1D sampling theorems of Nyquist, Gori, and Ding are special cases of this new general theorem. The discrete transform obtained by sampling the LCT kernel directly has been discussed. It has been shown that using the sampling rates derived in this paper, this discrete transform is unitary, which is as indicated important in multistage and iterative calculations. Future work will continue to focus on issues relevant to the numerical approximation of the 2D-NS-LCT. One topic of interest arises because we have found that in some cases, the FFT can be employed to calculate the 2D-NS-LCT efficiently. An example of one such case is shown in Fig. 5. In this case, it is possible to use the samples inside the new rectangular grid (red-dashed line square) instead of those inside the original irregular sampling grid (blue solid line parallelogram) with no loss of information. It seems that under certain circumstances the grids align and that (as is made clear by examining the color coding of the samples in Fig. 5) all of the sampled values appears appropriately in the corresponding replicas—a most fortuitous situation! Currently attempts are being made to identify the general sampling conditions under which the 2D-NS-LCT can be calculated using the fast FT.
APPENDIX A In this appendix, the detailed matrix decompositions of the phase term and the coordinate term in Eq. (9) are presented. First, the coordinate transform can be further decomposed using the following multiplication of the shearing transforms and the scaling transforms:
B
0
0
BT−1
1 shx
B B0 B B B0 @ 0 0
0
0
1
0
0
1
0
−shx
1∕sx B B 0 B ×B B 0 @ 0
0
0
1
0
0 sx 0
0
10
1
0 0
0
1
CB C B 0 C 0C CB shy 1 0 C CB C B 0 0 1 −shy C 0C A@ A 1 0 0 0 1 1 10 1 0 0 0 0 C CB B C 0C CB 0 1∕sy 0 0 C C; CB B 0 1 0C 0C A A@ 0 0 0 0 sy 1 (A1) 1 b22 .
and sy The first where shx , shy two matrices in Eq. (A1) represent shearing operations along x and y axes, respectively [24], while the last two matrices are scaling operations along the x and y axes [24]. Second, the phase term in the space domain can be decomposed into a coupling operation (between x and y) followed by a chirp transform. Similarly, the phase term in the LCT domain can be decomposed into another coupling operation (between x0 and y0 ) and a 2D chirp; see Eq. (A2): I 0 b12 b22
B−1 A 0 B B B B B B B B @
b21 b22 detB , sx
I 1
0
0
1 p2 2 detB
1
0
1
0
0 0
0
1 k2 2 detB
p2 2 detB I 0
B B B B B B B B @
10
CB B 0 0C CB 0 CB p1 B 1 0C CB detB CB CB A@ 0 1 0
0
DB−1 0
0 0
b22 detB,
I
0 k2 2 detB
0
10
1
CB B 0 0C CB 0 CB k1 B 1 0C CB detB CB CB A@ 0 1 0
0
0
1
0
0
1
p3 0 detB
0
0
1
0
0
1
k3 0 detB
0
1
C 0C C C 0C C C C A 1
0
1
C 0C C C 0C C: C C A 1 (A2)
APPENDIX B In this appendix, the mathematical derivation of the sampling theorem for the detB ≠ 0 case developed in this paper is discussed. From Eq. (12), the center point of the zero-order term ˆ 0 ; y0 is defined by xˆ 0 ; yˆ 0 0;0 0; 0. From (n m 0) of Gx Fig. 2(b) we know that Lx0 is the distance between (0,0) and xˆ 0 ; yˆ 0 1;0 bT11x ; bT21x ; Ly0 is the distance between (0,0) and Fig. 5. One example of the 2D-NS-LCT of a sampled image: The original output LCT image parallelogram shape is indicated by a solid bold (blue) line, while its replicas are shown bounded by thin (black) lines. The parallelogram shapes outline the skewed sample grid. In this example all the sampling points are located in a regular Cartesian coordinate system. As a result all the color coded points (yellow, purple, light blue, and red) in the central parallelogram are found suitably arranged inside the new central square shape, bound by the thickdashed (red) line.
xˆ 0 ; yˆ 0 0;1 bT12y ; bT22y , which indicates that s 2 2 q b11 b21 1 b211 b221 ; −0 −0 Lx0 Tx Tx Tx s 2 2 q b12 b22 1 Ly0 b212 b222 −0 −0 Ty Ty Ty
(B1)
(B2)
Zhao et al.
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
2639
APPENDIX C In this appendix, the mathematical derivation of the sampling theorem for the B 0 case developed in this paper is discussed. The shape or boundary of the new output image can be defined by the corner co-ordinate,
x01
x02
x03
x04
y01
y02
y03
y04
1 2 detD
1 2 detD
d22
−d21
L
x
Lx
−Lx
−Lx
Ly −Ly −Ly Ly −d12 d11 d L −d L d22 Lx d21 Ly −d22 Lx d21 Ly 22 x 21 y −d12 Lx d11 Ly
−d12 Lx − d11 Ly
d12 Lx − d11 Ly
−d22 Lx − d21 Ly : d12 Lx d11 Ly
(C1)
Then, q d212 d222 ;
Lx0
q x02 − x03 2 y02 − y03 2
Lx detD
Ly 0
q x01 − x02 2 y01 − y02 2
Ly q d211 d221 : detD
(C2a)
(C2b)
APPENDIX D In this appendix, the mathematical proof regarding how the DLCT is made unitary is presented. The multiplication of W M and W M −1 is shown as W M −1 W M
N X N X
W M −1 n0 m0 kl W M klnm
k1 l1
jπp1 n2 − n02 L2x p2 nm − n0 m0 Lx Ly p3 m2 − m02 L2y 1 exp ; detB N 2 detB
8 n o n o9
> > kb lb kb lb > > 0 0 0 0 11 12 21 22 N X N > N detB > > k1 l1 : ;
jπp1 n2 − n02 L2x p2 nm − n0 m0 Lx Ly p3 m2 − m02 L2y 1 exp ; detB N 2 detB N N X −j2πkn − n0 −j2πkm − m0 X −j2πln − n0 −j2πlm−0 exp exp : exp exp N N N N k1 l1
It can be seen that when n n0 , m m0 , Eq. (D1) can be rearranged as W M −1 W M
N2 I; detB
(D1)
support of the Science Foundation Ireland, the Irish Research Council, and Enterprise Ireland under the National Development Plan.
(D2)
where I is a 4D identity matrix. Thus we modify Eq. (D1) by N dividing every element of the kernel by p , and in this padetB
per we treat V M as this normalized version of the kernel map W detB . This modification makes our trix W M , namely V M M N definition of the 2D-NS-DLCT unitary if we use the sampling specified by Eq. (15). Thus, after modification, V M −1 V M I, and hence V M −1 V −1 M.
ACKNOWLEDGMENTS L. Zhao is supported by a UCD—China Scholarship Council (CSC) joint scholarship and a SPIE Optics and Photonics Education Scholarship. J. J. Healy acknowledges the support of the National University of Ireland (NUI) Postdoctoral Fellowship in the Sciences. The authors also acknowledge the
REFERENCES 1. M. J. Bastiaans, “Wigner distribution function and its application to first-order optics,” J. Opt. Soc. Am. 69, 1710–1716 (1979). 2. S.-C. Pei, “Discrete fractional Fourier transform based on orthogonal projections,” IEEE Trans. Signal Process. 47, 1335–1348 (1999). 3. A. Stern, “Why is the linear canonical transform so little known?” AIP Conf. Proc. 860, 225–234 (2006). 4. B. Deng, R. Tao, and Y. Wang, “Convolution theorems for the linear canonical transform and their application,” Info. Sci. 49, 592–603 (2006). 5. A. Stern, “Uncertainty principles in linear canonical transform domains and some of their implications in optics,” J. Opt. Soc. Am. A 25, 647–652 (2008). 6. B. M. Hennelly and J. T. Sheridan, “Generalizing, optimizing, and inventing numerical algorithms for the fractional Fourier, Fresnel, and linear canonical transforms,” J. Opt. Soc. Am. A 22, 917–927 (2005).
2640 7. 8. 9. 10. 11. 12.
13. 14. 15. 16. 17.
18. 19. 20.
21. 22.
23. 24. 25.
26. 27. 28. 29. 30.
31.
J. Opt. Soc. Am. A / Vol. 31, No. 12 / December 2014 A. Stern, “Sampling of linear canonical transformed signals,” Signal Processing 86, 1421–1425 (2006). A. Koç, H. M. Ozaktas, and L. Hesselink, “Fast and accurate computation of two-dimensional non-separable quadratic-phase integrals,” J. Opt. Soc. Am. A 27, 1288–1302 (2010). T. Alieva and M. J. Bastiaans, “Properties of the linear canonical integral transformation,” J. Opt. Soc. Am. A 24, 3658–3665 (2007). R. Simon and K. B. Wolf, “Structure of the set of paraxial optical systems,” J. Opt. Soc. Am. A 17, 342–355 (2000). E. G. Abramochkin and V. G. Volostnikov, “Generalized Gaussian beams,” J. Opt. A 6, S157–S161 (2004). L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, “Orbital angular momentum of light and the transformation of Laguerre Gaussian laser modes,” Phys. Rev. A 45, 8185–8189 (1992). J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Experimental implementation of the gyrator transform,” J. Opt. Soc. Am. A 24, 3135–3139 (2007). J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Gyrator transform: properties and applications,” Opt. Express 15, 2190–2203 (2007). J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Applications of gyrator transform for image processing,” Opt. Commun. 278, 279–284 (2007). K. Wolf and T. Alieva, “Rotation and gyration of finite two dimensional modes,” J. Opt. Soc. Am. A 25, 365–370 (2008). S.-C. Pei and J.-J. Ding, “Properties, digital implementation, applications, and self image phenomena of the gyrator transform,” in 17th European Signal Processing Conference (EUSIPCO) (Curran Associates, Inc., 2009), pp. 441–444. J. Shamir, “Cylindrical lens systems described by operator algebra,” Appl. Opt. 18, 4195–4202 (1979). H. H. Arsenault, “A matrix representation for non-symmetrical optical systems,” J. Opt. 11, 87–91 (1980). G. Nemes and A. E. Siegman, “Measurement of all ten secondorder moments of an astigmatic beam by the use of rotating simple astigmatic (anamorphic) optics,” J. Opt. Soc. Am. A 11, 2257–2264 (1994). J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Optical system design for orthosymplectic transformations in phase space,” J. Opt. Soc. Am. A 23, 2494–2500 (2006). A. Sahin, H. M. Ozaktas, and D. Mendlovic, “Optical implementations of two-dimensional fractional Fourier transforms and linear canonical transforms with arbitrary parameters,” Appl. Opt. 37, 2130–2141 (1998). A. Sahin, M. A. Kutay, and H. M. Ozaktas, “Nonseparable two-dimensional fractional Fourier transform,” Appl. Opt. 37, 5444–5453 (1998). S.-C. Pei, “Two-dimensional affine generalized fractional Fourier transform,” IEEE Trans. Signal Process. 49, 878–897 (2001). P. Dong and N. P. Galatsanos, “Affine transformation resistant watermarking based on image normalization,” in IEEE International Conference on Image Processing (IEEE, 2002), Vol. 3, pp. 489–492. Z.-J. Liu, H. Chen, T. Liu, P.-F. Li, J.-M. Dai, X.-G. Sun, and S.-T. Liu, “Double-image encryption based on the affine transform and the gyrator transform,” J. Opt. 12, 035407 (2010). Y. Zhang, B.-Z. Dong, B.-Y. Gu, and G.-Z. Yang, “Beam shaping in the fractional Fourier transform domain,” J. Opt. Soc. Am. A 15, 1114–1120 (1998). B. M. Hennelly, D. Kelly, A. Corballis, and J. T. Sheridan, “Phase retrieval using theoretically unitary discrete fractional Fourier transform,” Proc. SPIE 5908, 59080D1 (2005). F.-C. Zhang, I. Yamaguchi, and L. P. Yaroslavsky, “Algorithm for reconstruction of digital holograms with adjustable magnification,” Opt. Lett. 29, 1668–1670 (2004). D.-Y. Wang, J. Zhao, F.-C. Zhang, G. Pedrini, and W. Osten, “High-fidelity numerical realization of multiple-step Fresnel propagation for the reconstruction of digital holograms,” Appl. Opt. 47, D12–D20 (2008). P. F. Almoro and S. G. Hanson, “Object wave reconstruction by speckle illumination and phase retrieval,” J. Eur. Opt. Soc. Rapid 4, 09002 (2009).
Zhao et al. 32. G.-H. Situ, J. P. Ryle, U. Gopinathan, and J. T. Sheridan, “Generalized in-line digital holographic technique based on intensity measurements at two different planes,” Appl. Opt. 47, 711–717 (2008). 33. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35, 237–246 (1972). 34. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). 35. K. B. Wolf, “Concepts from complex vector analysis and the Fourier transform,” in Integral Transforms in Science and Engineering (Plenum, 1979), Chap. 1, p. 16. 36. R. E. Ziemer, W. H. Tranter, and D. R. Fannin, Signals and Systems: Continuous and Discrete (Prentice-Hall, 1983). 37. R. Bracewell, “The discrete Fourier transform and the FFT,” in The Fourier Transform and its Applications (McGraw-Hill, 1999), Chap. 11, p. 261. 38. A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 2nd ed. (Prentice-Hall, 1999), Chaps. 8 and 9, pp. 541–692. 39. J. W. Goodman, Introduction to Fourier Optics (Roberts & Company, 2005). 40. U. Schnars and W. Jüptner, “Digital holography,” in Digital Holography: Digital Hologram Recording, Numerical Reconstruction, and Related Techniques (Springer, 2010), Chap. 3, pp. 46–48. 41. F. Gori, “Fresnel transform and sampling theorem,” Opt. Commun. 39, 293–297 (1981). 42. J.-J. Ding, “Research of fractional Fourier transform and linear canonical transform,” Ph.D. dissertation (National Taiwan University, 2001). 43. B. M. Hennelly and J. T. Sheridan, “Fast numerical algorithm for the linear canonical transform,” J. Opt. Soc. Am. A 22, 928–937 (2005). 44. D. P. Kelly, B. M. Hennelly, and C. McElhinney, “A practical guide to digital holography and generalized sampling,” Proc. SPIE 7072, 707215 (2008). 45. J. J. Healy, B. M. Hennelly, and J. T. Sheridan, “Additional sampling criterion for the linear canonical transform,” Opt. Lett. 33, 2599–2601 (2008). 46. J. J. Healy and J. T. Sheridan, “Sampling and discretization of the linear canonical transform,” Signal Processing 89, 641–648 (2009). 47. D. P. Kelly, B. M. Hennelly, N. Pandey, T. J. Naughton, and W. T. Rhodes, “Resolution limits in practical digital holographic systems,” Opt. Eng. 48, 095801 (2009). 48. J. J. Healy and J. T. Sheridan, “Fast linear canonical transforms,” J. Opt. Soc. Am. A 27, 21–29 (2010). 49. J. J. Healy and J. T. Sheridan, “Reevaluation of the direct method of calculating Fresnel and other linear canonical transforms,” Opt. Lett. 35, 947–949 (2010). 50. D. P. Kelly, J. J. Healy, B. M. Hennelly, and J. T. Sheridan, “Quantifying the 2.5 D imaging performance of digital holographic systems,” J. Eur. Opt. Soc. Rapid 6, 11034 (2011). 51. J. J. Healy and J. T. Sheridan, “Space-bandwidth ratio as a means of choosing between Fresnel and other linear canonical transform algorithms,” J. Opt. Soc. Am. A 28, 786–790 (2011). 52. J. T. Sheridan, L. Zhao, and J. J. Healy, “The condition number and first order optical systems,” in Imaging and Applied Optics Conference, OSA Technical Digest (OSA, 2012), paper ITu3C.4. 53. L. Zhao, J. J. Healy, and J. T. Sheridan, “Unitary discrete linear canonical transform: analysis and application,” Appl. Opt. 52, C30–C36 (2013). 54. J.-J. Ding and S.-C. Pei, “Eigenfunctions and self-imaging phenomena of the two-dimensional non-separable linear canonical transform,” J. Opt. Soc. Am. A 28, 82–95 (2011). 55. J.-J. Ding, S.-C. Pei, and C.-L. Liu, “Improved implementation algorithms of the two-dimensional non-separable linear canonical transform,” J. Opt. Soc. Am. A 29, 1615–1624 (2012). 56. T. Alieva and M. J. Bastiaans, “Alternative representation of the linear canonical integral transform,” Opt. Lett. 30, 3302–3304 (2005).
Zhao et al. 57. L. Zhao, J. J. Healy, and J. T. Sheridan, “Sampling of the two dimensional non-separable linear canonical transform,” Proc. SPIE 9131, 913112 (2014). 58. F. R. Gantmacher, “Matrices and operations on matrices,” in The Theory of Matrices (American Mathematical Society, 1959), Chap. 1, Vol. 1, pp. 1–7.
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
2641
59. F. R. Gantmacher, “Linear operators in an n-dimensional vector space,” in The Theory of Matrices (American Mathematical Society, 1959), Chap. 3, Vol. 1, pp. 55–56. 60. H.-J. Bao, A.-J. Sang, and H.-X. Chen, “Inverse operation of fourdimensional vector matrix,” Int. J. Intell. Syst. Appl. 3, 19–27 (2011).
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
2631
Two-dimensional nonseparable linear canonical transform: sampling theorem and unitary discretization Liang Zhao,1 John J. Healy,2 and John T. Sheridan1,* 1
School of Electrical, Electronic and Communications Engineering, Communications and Optoelectronic Research Centre, SFI–Strategic Research Cluster in Solar Energy Conversion, College of Engineering and Architecture, University College Dublin, Belfield, Dublin 4, Ireland 2 Department of Electronic Engineering, Maynooth University, Co. Kildare, Ireland *Corresponding author: [email protected] Received August 4, 2014; revised October 5, 2014; accepted October 6, 2014; posted October 6, 2014 (Doc. ID 214789); published November 11, 2014 The two-dimensional (2D) nonseparable linear canonical transform (NS-LCT) is a unitary, linear integral transform that relates the input and output monochromatic, paraxial scalar wave fields of optical systems characterized by a 4 × 4 ray tracing matrix. In addition to the obvious generalizations of the 1D LCT (which are referred to as separable), the 2D-NS-LCT can represent a variety of nonaxially symmetric optical systems including the gyrator transform and image rotation. Unlike the 1D LCT, the numerical approximation of the 2D-NS-LCT has not yet received extensive attention in the literature. In this paper, (1) we develop a sampling theorem for the general 2D-NS-LCT which generalizes previously published sampling theorems for the 1D case and (2) we determine which sampling rates may be chosen to ensure that the obvious discrete transform is unitary. © 2014 Optical Society of America OCIS codes: (000.4430) Numerical approximation and analysis; (070.2025) Discrete optical signal processing; (070.2590) ABCD transforms; (090.1995) Digital holography. http://dx.doi.org/10.1364/JOSAA.31.002631
1. INTRODUCTION The one-dimensional (1D) linear canonical transform (LCT) [1–5] is a three-parameter group of linear integral transforms, which can be used to model the propagation of the coherent wave fields through first-order optical systems. Among its special cases [6,7] are the Fourier transform (FT), the fractional Fourier transform (FRT), the Fresnel transform (FST), the scaling transform, and the chirp transform. The two-dimensional separable linear canonical transform (2D-S-LCT) can model a wide variety of isotropic optical systems [8], which can be treated as being made up of two independent 1D-LCTs. In the case when the system is fully symmetric and orthogonal, and the parameters for both of the two dimensions are identical, such systems can be represented using only three parameters as in the 1D-LCT case [9]. In the case when the system is still orthogonal but the parameters for the orthogonal dimensions differ, such systems can be represented by six free parameters [9,10]. The two-dimensional nonseparable linear canonical transform (2D-NS-LCT) is a more general integral transform, in which the two dimensions are coupled to each other by four additional cross-parameters, increasing the total number of free parameters to ten [8]. This general nonseparable transform is nonorthogonal, nonaxially symmetric, and anamorphic [8]. The 2D-NS-LCT is able to represent systems not only involving coupling and shearing between components along the different dimensions [9,11,12] but also involving rotations between any arbitrary planes in phase space, such as gyrations [13–17]. The 2D-NS-LCT is important in digital signal processing and optical system modeling [18–26]. 1084-7529/14/122631-11$15.00/0
Discrete transforms that approximate the continuous LCTs allow us to simulate a number of optical propagation problems. Digital simulations to solve these problems require access to a numerical algorithm to approximate the continuous transform efficiently, accurately, and unitarily [27,28]. In the area of optics not only is such an algorithm important when employing iterative techniques to design systems, but it is also significant when solving inverse scatter problems to extract information about an object given a known illumination [27–32], for instance in phase retrieval [33,34]. Traditionally, such problems have been addressed in situations where the FT and FST were sufficient to describe the optical system examined. Unitary algorithms for both the FT [35–40] and FST [41] are widely available. In order to further avoid aliasing during sampling processes, Ding [42] developed a sampling theorem for the general 1D-LCT in 2001 (also see [7]). In 2005, Hennelly and Sheridan [6] discussed how the space–bandwidth product of a signal (and thus the sampling requirements) can be tracked through an LCT system. However the resulting discrete transforms are not in general unitary. In [43–52], these results were shown to lead to a fast and accurate numerical implementation of the 1D-LCT. Recently a sufficient condition on the sampling rates chosen during discretization to ensure that the discrete 1D-LCT (1D-DLCT) is unitary have been presented [52,53]. Proof that these conditions guarantee unitarity is provided in [53]. It should be noted that the sampling theorem developed in [53] is consistent with the Shannon sampling theorem for the FT, the Gori sampling theorem [41] for the FST, and the Ding sampling theorem [42] for the LCT. © 2014 Optical Society of America
2632
J. Opt. Soc. Am. A / Vol. 31, No. 12 / December 2014
Zhao et al.
Numerical calculation of the 2D-NS-LCT has recently received attention [8,17,22–24,54–57]. This paper will focus on developing a sampling theorem for this case and demonstrating how a unitary sampling rate can be chosen. This paper is organized as follows: In Section 2 we introduce the transforms examined for 1D and 2D signals and provide some other important definitions. One powerful technique used to analyze LCTs is to decompose them into, and then apply sequentially, a set of simpler transforms. In Section 3, we introduce this kind of decomposition, discuss some important special cases, and look at the insight one such decomposition gives us about 2D-NS-LCT’s application to discrete signals. In Section 4, the general sampling theorem for the 2DNS-LCT is developed. In Section 5, we consider one DLCT, a discrete transformation that approximates the continuous integral transform. In particular the parameters for which the DLCT is unitary are determined. Finally, a brief conclusion and future work are given in Section 6.
2. DEFINITIONS OF THE CONTINUOUS LCTs A. Definitions of the 1D-LCT and Its Unitary Property The continuous 1D-LCT of a signal gx can be defined by [5,6,53] Gx0 LM fgxgx0 h i ( R 2 0 Dx02 p1 ∞ exp jπAx −2xx gxdx; −∞ B pjB 02 0 D expjπCDx gDx ;
where k1 d11 b22 − d12 b21 ;
k2 2−d11 b12 d12 b11 ;
k3 −d21 b12 d22 b11 ; p2 2a12 b22 − a22 b12 ;
p1 a11 b22 − a21 b12 ; p3 −a12 b21 a22 b11 :
In this case the transformation matrix M of the system is defined as 0
a11 B a21 B M @ c11 c21
a12 a22 c12 c22
b11 b21 d11 d21
1 b12 b22 C C A d12 A C d22
B ; D
B0
; (1)
where given M f A B ; C D g, the parameters of the transform matrix M are restricted to be real numbers with detM 1, and detM represents the determinant of matrix M. For the case of B 0 the LCT is defined as the limit of the B ≠ 0 case with B → 0 [5]. The LCT with parameters fA B; C Dg f0 1; −1 0g corresponds to the FT, while fcos θ sin θ; −sin θ cos θg corresponds to the FRT. The FST results for f 1 z ; 0 1 g. In addition, the transform matrix of the inverse 1D-LCT is given by M −1 f D −B ; −C A g.
(3-3)
where A, B, C, and D are now 2 × 2 submatrices, and detB ≠ 0. We note that the 2D-NS-LCT has 16 parameters, see Eq. (3-3), which must satisfy the constraints ABT BAT , CDT DC T , and ADT − BC T I, or AT C C T A, BT D DT B, and AT D − C T B I, where I is a 2 × 2 identity matrix. Since the above equations are equivalent to six constraints, the number of independent parameters in the matrix M is in fact reduced to ten [8,54–57]. In addition, the inverse 2D-NS-LCT of Gx0 ; y0 is given by gx; y LM −1 fGx0 ; y0 gx; y;
B≠0
(3-2)
where M −1 f DT
−BT ; −C T
(4)
AT g.
2.2 Other Definitions It should be noted that ① When B 0; ⇒ detB 0, the definition of the 2D-NSLCT is [9,24,55,56] Gx0 ; y0 LM fgx; ygx0 ; y0 p detD expjπc11 d11 c12 d12 x02 × expj2πc11 d21 c12 d22 x0 y0
2.1 Unitary Property One of the most important properties of the continuous LCT (both 1D and 2D cases) is the unitary property [3],
× expjπc21 d21 c22 d22 y02 gd11 x0 d21 y0 ; d12 x0 d22 y0 ;
LM −1 fLM fggg g:
where M f DT 0 ; C D g. ② When detB 0 but B ≠ 0, the 2D-NS-LCT is more complicated; see [24].
(2)
B. Definitions of the 2D-NS-LCT The continuous 2D-NS-LCT of a signal gx; y can be defined by [8,54–57] Gx0 ; y0 LM fgx; ygx0 ; y0 ∞ ZZ 1 jπk1 x02 k2 x0 y0 k3 y02 p exp detB j detB −∞ j2π−b22 x0 b12 y0 x b21 x0 − b11 y0 y × exp detB jπp1 x2 p2 xy p3 y2 gx; ydxdy; (3-1) × exp detB
(5) −1
C. Special Cases of the 2D-LCT 1. Special Cases of the 2D-S-LCT In the case when a12 a21 b12 b21 c12 c21 d12 d21 0;
(6)
the 2D-NS-LCT, see Eqs. (3-1) and (5), reduces to the separable transform, the well-known special cases of which are the 2D FT, the 2D FST [21], the 2D FRT [21,22], the 2D chirp transform [18,19,24], and the 2D scaling transform [24]. The transform matrix M of these special cases is shown in Table 1.
Zhao et al.
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
2633
Table 1. Transform Matrices of Some Well-known Cases of the 2D-S-LCT and 2D-NS-LCT 2D-Separable
0
0 B 0 B @ −1 0 0 1 B0 B @0 0
FT
FST
0
cos θx B 0 B @ − sin θx 0 0
FRT
Chirp
Scaling-x
Scaling-y
0 0 0 −1 0 1 0 0
1 0 0 0 z 0 1 0
1 0 1C C 0A 0 1
0 zC C 0A 1
0 0 1 0
Gyrator
Coupling
0 sin θx cos θy 0 0 cos θx − sin θy 0
1 0 B 0 1 B @ −1∕f x 0 0 −1∕f y 0 1∕sx 0 0 B 0 1 0 B @ 0 0 sx 0 0 0 0 1 0 0 B 0 1∕sy 0 B @0 0 1 0 0 0
1 0 sin θy C C A 0 cos θy 1
Coordinate (general)
0 0C C 0A 1 1
Rotation
0 0C C 0A 1 1 0 0C C 0A sy
Shearing-x
Shearing-y
2. Special Cases of the 2D-NS-LCT Three well-known 2D nonseparable transforms are the gyrator transform [13–17,21], the coupling (i.e., multiplication) transform [24], and the homogenous coordinate/affine transform [24–26], which can be defined by LM gyrator fgx;yg 1 jsin θj
Z∞ Z exp −∞
2D-Nonseparable 1 0 cos θ 0 0 sin θ B 0 cos θ sin θ 0 C C B @ 0 − sin θ cos θ 0 A − sin θ 0 0 cos θ 1 0 1 0 0 0 B0 1 0 0C C B @0 τ 1 0A τ 0 0 1 1 0 0 0 a11 a12 B a21 a22 0 0 C B a22 −a21 C A @ 0 0 detA detA −a12 a11 0 0 detA detA 1 0 cos θ sin θ 0 0 B − sin θ cos θ 0 0 C C B @ 0 0 cos θ sin θ A 0 0 − sin θ cos θ 1 0 0 0 1 shx B0 1 0 0C C B @0 0 1 0A 0 0 −shx 1 1 0 1 0 0 0 B shy 1 0 0 C C B @ 0 0 1 −shy A 0 0 0 1
j2πxy x0 y0 cos θ − xy0 yx0 gx;ydxdy; sin θ (7-1)
3. MATRIX DECOMPOSITION AND CONSEQUENCES FOR THE LOCATIONS OF THE SAMPLES IN THE OUTPUT PLANE The ABCD matrices of concatenated systems may be obtained by multiplication; for example, if a wave field propagates through a system with ABCD matrix M 1 and then a second with matrix M 2 , the overall system’s ABCD matrix is given by M 3 M 1 M 1 . We can use this result in two ways: First, we can find the output of a series of systems in a single calculation. Second, we can decompose a calculation into a sequence of simpler operations. A simple example is the venerable spectral method for calculating the FST,
LM Coupling xy fgx; yg expjπτxygx; y;
LM Coupling x0 y0 fGx0 ; y0 g expjπτx0 y0 Gx0 ; y0 ;
(7-2-1)
(7-2-2)
a x0 − a12 y0 −a21 x0 a11 y0 ; : (7-3) LM coordinate fgx; yg g 22 detA detA The transform matrices of these three cases are listed in Table 1. For brevity, the homogenous coordinate transform is referred to as the “coordinate transform” below. The coordinate transform includes some well-known transforms, like rotation [8,20,21,24] and shearing (interferometer) [24]; see Table 1.
1 λz 0 1
0 −1 1 0
1 0 −λz 1
0 −1
1 : 0
(8)
Equation (8) shows that the FST may be calculated by taking a FT, a chirp multiplication, and an inverse FT, all of which may be calculated using simple or mature software packages. Two decompositions have been proposed in the literature to numerically calculate the 2D-NS-LCT, the Iwasawa decomposition [8] and a second decomposition which may be thought of as a generalization of the direct method for calculating the 1D case [55]. This second decomposition is given by M M 4M 3M 2M 1 B 0 0 I 0 −1 T −1 0 B −I DB I where
I 0
I
−1
B A
0 ; I
(9-1)
2634
J. Opt. Soc. Am. A / Vol. 31, No. 12 / December 2014
Zhao et al.
0
1 k1 k2 B detB 2 detB C C; DB−1 B @ A k2 k3 2 detB detB 0 1 p1 p2 B detB 2 detB C C; B−1 A B @ A p2 p3 2 detB detB 0 b −b21 1 22 B detB detB C −1 C; BT B @ −b A b 12
detB
and 1(e). The relevant sampling relationships are discussed in Section 4. The decomposition in Eq. (9) is valid only if detB 0. No general result exists at this time for the detB 0 case. However, if we further specify that B 0, we can decompose the matrix as follows: M M 2M 1 (9-2)
11
detB
and where k1 ∼ k3 and p1 ∼ p3 are as defined in Eq. (3-2). M 1 and M 4 in Eq. (9-1) produce phase terms in the space and LCT domains, respectively; M 3 generates a coordinate transform; and M 2 is an FT. The matrix decompositions are discussed in detail in Appendix A. This decomposition would appear to be an excellent candidate for numerical approximation of the 2D-NS-LCT. It consists of two generalized chirps (M 1 and M 4 ), a 2D FT (M 2 ) and an image coordinate transformation (M 3 ). This last operation indicates an interesting and vital feature of such numerical approximations. Given an input comprised of a rectangular grid of regular samples, we normally expect a discrete transform to relate this input with a rectangular grid of output samples. However, in the discrete 2D-NS-LCT case, it will in general result in a skewed output grid. This is illustrated in Fig. 1. This has significant implications for an appropriate sampling theorem. Consider the 1D case. Given a figure’s width in the input and output domains, it has been shown that the appropriate sampling rates are related to these widths. In the separable 2D case, a similar result is obtained given the input and output widths in both directions; see Figs. 1(b) and 1(d). However, for the nonseparable case, the output samples obtained form a parallelogram of points; see Figs. 1(c)
I CDT
0 I
DT 0
−1
0 ; D
(10-1)
d22 −d21 −d12 d11 where DT−1 f detD detD ; detD detD g. Here M 2 produces a phase term in the LCT domain, while M 1 generates the following coordinate transform,
LM 1 fgx; ygx0 ; y0 gd11 x0 d21 y0 ; d12 x0 d22 y0 : (10-2) Consequently, beginning with a rectangular grid of samples, once again one ends up with a parallelogram arrangement of output samples.
4. SAMPLING THEOREM FOR THE 2D-NS-LCT A. Sampling Effects We consider the impulse sampled version of gx; y, gˆ x; y gx; y
∞ X
∞ ∞ X X
δx − nT x ; y − mT y
n−∞ m−∞ ∞ X
gnT x ; mT y δx − nT x ; y − mT y
n−∞ m−∞
∞ ∞ X X 1 j2πnx j2πmy exp ; gx; y exp T xT y Tx Ty n−∞ m−∞ (15) where T x and T y are the sampling intervals in the x and y directions, respectively. 1. The Effect of Sampling in the detB ≠ 0 Case From Eq. (3), the 2D-NS-LCT of gˆ x; y can be defined by ˆ 0 ; y0 LM fˆgx;ygx0 ; y0 Gx ∞ ∞ X X 1 p H ×G T x T y j detB n−∞ m−∞ nb11 mb12 0 nb21 mb22 0 × x− − ;y − − Tx Ty Tx Ty ∞ ∞ X X 1 p H × Gx0 − xˆ 0n;m ; y0 − yˆ 0n;m ; T x T y j detB n−∞ m−∞ (12-1)
Fig. 1. Input image shape, output image shapes, and their sampling grids: (a) rectangular shape of an input image, with length Lx and height Ly . (b) and (c) are the output image shapes of the 2D S-LCT and NS-LCT, respectively; the output image shape in the separable case is still a rectangle, while the output image shape in the nonseparable case is skewed, that is, becomes a parallelogram (the side lengths of which are different from the values in the separable case). (d) and (e) are the sampling grids of (b) and (c), respectively.
where xˆ 0n;m and
nb11 mb12 ; Tx Ty
and yˆ 0n;m
nb21 mb22 ; Tx Ty
(12-2)
Zhao et al.
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
−jπk1 x0 − xˆ 0 2 k2 x0 − xˆ 0 y0 − yˆ 0 k3 y0 − yˆ 0 2 detB 02 0 0 jπk1 x k2 x y k3 y02 × exp : (12-3) detB
H exp
Equation (12-1) indicates that the 2D-NS-LCT of gˆ x; y, ˆ 0 ; y0 , is made up of modulated shifted replicas namely Gx of the original signal Gx0 ; y0 [7,53,57]. The 2D-NS-LCT of a sampled signal is therefore a periodically replicated version of the 2D-NS-LCT of the original continuous signal before sampling. From Eq. (12-2), the center point (coordinates) of the n; mth replica are designated by xˆ 0 ; yˆ 0 n;m , the location of which is dependent on T x ; T y and on the values of all four parameters of the 2 × 2 matrix B in the 2D-NS-LCT. It should be noted that, in the separable case, the xˆ 0 ; yˆ 0 n;m points are located on a Cartesian coordinate system; see the dashed lines in Fig. 2(a). In the nonseparable case, the xˆ 0 ; yˆ 0 n;m points are located along a series of parallel lines; see Fig. 2(b). 2. Effect of Sampling in the B 0 Case From Eq. (8) [24], as in the former detB ≠ 0 case, the 2D-NS-LCT of gˆ x; y can be defined by
B. Sampling Theorem The sampling rate should be such that the replicas indicated by Eq. (12-1) do not overlap; then the continuous signal may be recovered from the samples taken at this rate. 1. Sampling Theorem for the detB ≠ 0 Case Using Eq. (12), and in order to avoid aliasing, the output LCT image extents should satisfy Lx0 ≤
(13) It should be noted that ① In the case when detB ≠ 0, the 2D-NS-LCT of the ˆ 0 ; y0 , is still a continusampled input signal gˆ x; y, namely Gx ous signal; see Fig. 2b. ② In the case when B 0, the 2D-NS-LCT of the sampled input signal gˆ x; y is a discrete signal. The output S-LCT samples and NS-LCT samples are arranged as shown in Figs. 1(d) and 1(e).
1 Tx
q b211 b221 ;
and Ly0 ≤
1 Ty
q b212 b222 ;
(14-1)
q b211 b221 is the distance between the n-direction q replicas, while T1y b212 b222 is the distance between the where
1 Tx
m-direction replicas [see Fig. 2(b)], in which θx0 arctanbb21 11 and θy0 arctanbb22 . The derivations of these two distances 12 are discussed in Appendix B. 2. Sampling Theorem for the B 0 Case When detB 0, the effect of the LCT on the shape change is discussed in Appendix A. The sampling theorem for this case is
ˆ 0 ; y0 Gx
∞ ∞ X X nd md12 nd md22 δ y0 − 21 − : Gx0 ; y0 δ x0 − 11 − Tx Ty Tx Ty n−∞ m−∞
2635
Lx0 ≤
q Lx d212 d222 detD
;
and Ly0 ≤
q Ly d211 d221 detD
;
(14-2)
where θx0 − arctandd12 and θy0 − arctandd11 ; see Fig. 2(b). 22 21 In conclusion, as long as Eq. (14) is satisfied, Gx0 ; y0 can be reconstructed from the periodically replicated version of ˆ 0 ; y0 . Gx0 ; y0 , namely Gx C. Sampling Theorem for Special Cases The sampling theorems for some well-known special cases are shown in Table 2. The well-known Nyqist sampling theorem for the FT, the Gori sampling theorem [41] for the FST, and the Ding sampling theorem [42] for the 1D-LCT are special cases of the theorem developed here.
5. UNITARY DISCRETE TRANSFORM A. Unitary Discretization During discretization, the continuous coordinates x; y are replaced by the n; mth sample locations x; yn;m , that is, the coordinates nT x ; mT y , while x0 ; y0 are replaced by the n0 ; m0 th sample x0 ; y0 n0 ;m0 , where the integers n, m, n0 , m0 ∈ −N∕2 1; N∕2. N is the number of sampling points, which is assumed to here to be even valued. When N is odd, the analysis is similar and does not warrant separate attention. Fig. 2. (a) and (b) are the outputs of a 2D-S-LCT and a 2D-NS-LCT, respectively, given a discrete input. Were the inputs continuous, we would obtain outputs at locations indicated by the shaded (blue where color is available) tetrahedrons, while the other tetrahedrons indicate the locations of replicas that appear because of the discrete input. Were we to sample the outputs, we should do so on grids corresponding to Figs. 1(d) and 1(e) located in the shaded regions.
1. Unitary Discretization for the detB ≠ 0 Case The corresponding discrete transform is unitary, if the n0 ; m0 th sample x0 ; y0 n0 ;m0 is defined by x0n0 ;m0
n0 b11 m0 b12 ; Lx Ly
and y0n0 ;m0
n0 b21 m0 b22 ; (15) Lx Ly
2636
J. Opt. Soc. Am. A / Vol. 31, No. 12 / December 2014
Zhao et al.
Table 2. Sampling Theorem and Unitary Sampling Rates for Some Well-known Cases of the 2D-S-LCT and 2D-NS-LCT
S-LCT when B ≠ 0
S-LCT when B 0
LCTs
Sampling Theorems
Unitary Rate Chosen
2D-FT
Lx0 ≤ 1∕T x , and Ly0 ≤ 1∕T y
T x0 1∕Lx , and T y0 1∕Ly
2D-FRT
Lx0 ≤ sin θx ∕T x , and Ly0 ≤ sin θy ∕T y
T x0 sin θx ∕Lx , and T y0 sin θy ∕Ly
2D-FST
Lx0 ≤ λz∕T x , and Ly0 ≤ λz∕T y
T x0 λz∕Lx , and T y0 λz∕Ly
2D-S-LCT
Lx0 ≤ b11 ∕T x , and Ly0 ≤ b22 ∕T y
T x0 b11 ∕Lx , and T y0 b22 ∕Ly
2D Chirp
Lx0 ≤ Lx , and Ly0 ≤ Ly
T x0 T x , and T y0 T y
Scaling-x
L ≤ Lx ∕S x , and L ≤ Ly
T x0 T x ∕S x , and T y0 T y
x0
y0
Scaling-y
Lx0 ≤ Lx , and Ly0 ≤ Ly ∕S y
T x0 T x , and T y0 T y ∕S y
NS-LCT when B ≠ 0
Gyrator
NS-LCT when B 0
Coordinate
Lx0 ≤ sin θ∕T y , and Ly0 ≤ sin θ∕T x q q Lx0 ≤ Lx a211 a221 , and Ly0 ≤ Ly a212 a222
T x0 sin θ∕Ly , and T y0 sin θ∕Lx q q T x0 T x a211 a221 , and T y0 T y a212 a222
Lx0 ≤ Lx , and Ly0 ≤ Ly p Lx0 ≤ Lx , and Ly0 ≤ Ly 1 sh2x q Lx0 ≤ Lx 1 sh2y , and Ly0 ≤ Ly
T x0 T x , and T y0 T y p T x0 T x , and T y0 T y 1 sh2x q T x0 T x 1 sh2y , and T y0 T y
Coupling Shearing-x Shearing-y
2. Unitary Discretization for the B 0 Case The corresponding discrete transform is unitary if
where Lx and Ly are the extents (i.e., length and height) of the input image; see Fig. 1(a). The input and output samples each form N × N vectors. The sampling intervals in the space domain are T x Lx ∕N and T y Ly ∕N. The kernel of the discrete transform is a 4D matrix [58–60], namely an N × N × N × N matrix, which is denoted as W M . The final discrete transform is
y0n;m
Gx0n0 ;m0 ; y0n0 ;m0
N∕2 X
N∕2 X
W M n0 ; m0 ; n; mgnT x ; mT y :
nd22 T x − md21 T y ; and detD −nd12 T x md11 T y : detD
x0n;m
(17)
B. Unitary Properties The continuous 2D-NS-LCT is unitary, but discretization can destroy this property. Mathematically, a square matrix W M is unitary if
(16-1)
n−N∕21 m−N∕21
WH M W M I; The operation in Eq. (16-1) corresponds to a matrix multiplication. From Eq. (15), the elements of the kernel W M are given by
(18-1)
where W H M represents the conjugate transpose of W M , and I is an identity matrix, that is,
8 h
2 0
0
0
2 i9 > > 0b 0b 0b 0b 0b n m n b m n b m n b m > > 11 12 11 12 21 22 21 22 = j detB n−N∕21 m−N∕21 > :
× exp
x
y
x
y
x
y
x
y
> > ;
detB
8 nh h 0 0
0
i
0
i o9 0b 0b 0b 0b > > n b m n b m n b m n b m > 11 12 21 22 11 12 21 22 = x
y
x
y
> > :
x
y
x
jπp1 nLx 2 p2 nLx mLy p3 mLy 2 : × exp N 2 detB
It should be noted that in the separable case, the sampling points in the LCT domain are located in a rectangular grid; see Fig. 1(d). In the nonseparable case, the sampling points are located in an irregular grid, namely a skewed grid; see Fig. 1(e).
y
> > ;
detB
(16-2)
WH M W M −1 :
(18-2)
① For the 1D-DLCT, this identity matrix is 2D. ② In the general 2D-NS-DLCT, the matrix is 4D [58–60]. It should be noted that since the two dimensions in the
Zhao et al.
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
2637
2D-S-DLCT can be treated independently, by means of the row–column algorithm [53], the identity matrix in separable cases can be reduced to a 2D matrix [53]. Significantly, unlike the corresponding continuous LCTs in Eq. (2), the DLCTs are not necessarily unitary; that is, without due care discretization can destroy this property. In fact Eq. (16-2) must be multiplied by a scaling factor in order to make it unitary; see Appendix D. First, in order to show how the sampling grid can be skewed by the 2D-NS-LCT, one example is given in Fig. 3. In this example, the sample locations are calculated using Eq. (15), which guarantees that the corresponding discrete transform is unitary. As we discussed in Section 3, the transform can be decomposed in a number of ways. Of particular interest for this discussion is Eq. (9), which shows that we can decompose any 2D LCT where B−1 exists into several elements which do not alter the sampling grid and one coordinate transform. In the example given in Fig. 3, we have chosen LCT parameters such that the coordinate transform is simply a rotation. Consequently, the sampling grid in the LCT domain is rotated, specifically by 60° clockwise; see Fig. 3(b) (where θx0 −60° and θy0 30°). In Section 4 of this paper, we have shown that the sampling rates that result in a unitary discrete transform are integer multiples of the rate given by Ding’s sampling theorem for the LCT. One might reasonably wonder if the gains made by using this transform are simply the result of using an appropriate sampling rate, if indeed unitarity is of real importance. This is not the case. We merely constrain the relationship between the input and output sampling rates; this is signal-independent. In order to demonstrate this, we have generated Fig. 4. This figure shows the error that arises when using the same input signal and LCT used to generate Fig. 3. This signal is sampled and then forward transformed. Finally the resulting output is inverse transformed and compared to the original input. Were the improvement merely a matter of a higher sampling rate reducing aliasing, we would expect to see a monotonic decrease in error with decreasing sampling period. Instead, we observe a “sweet spot” where the transform is unitary. We vary the sampling period of x0 from
Fig. 4. If we propagate to a given LCT domain and then back to the input domain, how does the error vary with sampling period? In this plot, we show the MSE is minimized, and in fact is 0, only at the sampling rate we propose, see Eq. (15) , i.e., the rate at which the DLCT is unitary.
95% to 105% of the rate we propose, which gives unitarity. The mean square error (MSE) between the transformed signal, g¯ x; y, and the original signal, gx; y, is calculated by MSE
(19)
where g¯ is a N × N array of complex values. It can be seen that (1) the minimum MSE value is obtained by choosing the unitary sampling rate defined in Eq. (15); (2) as the difference between the unitary and nonunitary sampling rates increases, the corresponding MSE increases rapidly. C. Special Cases of the Unitary 2D-NS-DLCT 1. Some Special Cases of the Unitary S-LCT For the separable cases when detB ≠ 0, Eq. (15) can be simplified as x0n0 ;m0
n0 b11 ; Lx
and y0n0 ;m0
m0 b22 : Ly
(20)
After dividing each side by N, Eq. (20) becomes T x0
Fig. 3. One example of unitary calculation using the DLCT: (a) the input signal is a 2D Gaussian, 20 × 20 pixels, variances 0.004 and 0.002. The extents (Lx and Ly ) are both 0.05. (b) The DLCT for parameters [0.08 0.09 0.5 0.866; 0.1 −0.0858 −0.8660.5 0.5; −0.495 −0.858 0.05 0.06; 0.856 −0.494 0.07 −0.0558] was calculated using the formula in Eq. (16) and plotted using the contour function in MATLAB. If we decompose the matrix of parameters as per Eq. (9), we find the coordinate transformation matrix is a rotation matrix, which is accounted for in the plotting of the results. According to Subsection 4.B, here arctan−0.866 and θy0 arctanbb22 θx0 arctanbb21 0.5 −60°, 11 12 0.5 30°, that is, the rectangular sampling grid in the LCT arctan0.866 domain is rotated 60° clockwise.
PP jgx; y − g¯ x; yj2 PP ; jgx; yj2
b11 ; Lx
and T y0
b22 : Ly
(21)
Equation (21) is in fact Ding sampling [42]. For the 2D-S-LCT case, we have offered a sufficient condition to ensure it is unitary [53]. Unitary sampling rates for some special cases of the 2D-S-LCT are shown in Table 2. In order to obtain the kernal in the 2D-S-DLCT case when detB ≠ 0, we only need to set b12 b21 0 in Eq. (16-2). However, this approach does not provide a fast way to evaluate the continuous transform. The row–column algorithm provides an efficient way to numerically approximate the 2D-S-LCT [53], in which the 2D transform is reduced to two 1D-LCTs. The corresponding sampling points in the LCT domain are located in a rectangular grid; see Fig. 1(d). 2. Some Special Cases of the Unitary NS-LCT The unitary sampling rates for some special cases are shown in Table 2.
2638
J. Opt. Soc. Am. A / Vol. 31, No. 12 / December 2014
Zhao et al.
0
6. CONCLUSION AND FUTURE WORK The 2D case of the LCT has not received as much attention as the 1D case, and in most 2D cases examined it has been assumed that the optical system is separable, thus permitting 1D analysis. In this paper the general 2D-NS-LCT is explored, namely those cases of the 2D-LCT explicitly not reducible to 1D-LCT analysis. A sampling theorem for this transform has been derived, and it has been shown that the 1D sampling theorems of Nyquist, Gori, and Ding are special cases of this new general theorem. The discrete transform obtained by sampling the LCT kernel directly has been discussed. It has been shown that using the sampling rates derived in this paper, this discrete transform is unitary, which is as indicated important in multistage and iterative calculations. Future work will continue to focus on issues relevant to the numerical approximation of the 2D-NS-LCT. One topic of interest arises because we have found that in some cases, the FFT can be employed to calculate the 2D-NS-LCT efficiently. An example of one such case is shown in Fig. 5. In this case, it is possible to use the samples inside the new rectangular grid (red-dashed line square) instead of those inside the original irregular sampling grid (blue solid line parallelogram) with no loss of information. It seems that under certain circumstances the grids align and that (as is made clear by examining the color coding of the samples in Fig. 5) all of the sampled values appears appropriately in the corresponding replicas—a most fortuitous situation! Currently attempts are being made to identify the general sampling conditions under which the 2D-NS-LCT can be calculated using the fast FT.
APPENDIX A In this appendix, the detailed matrix decompositions of the phase term and the coordinate term in Eq. (9) are presented. First, the coordinate transform can be further decomposed using the following multiplication of the shearing transforms and the scaling transforms:
B
0
0
BT−1
1 shx
B B0 B B B0 @ 0 0
0
0
1
0
0
1
0
−shx
1∕sx B B 0 B ×B B 0 @ 0
0
0
1
0
0 sx 0
0
10
1
0 0
0
1
CB C B 0 C 0C CB shy 1 0 C CB C B 0 0 1 −shy C 0C A@ A 1 0 0 0 1 1 10 1 0 0 0 0 C CB B C 0C CB 0 1∕sy 0 0 C C; CB B 0 1 0C 0C A A@ 0 0 0 0 sy 1 (A1) 1 b22 .
and sy The first where shx , shy two matrices in Eq. (A1) represent shearing operations along x and y axes, respectively [24], while the last two matrices are scaling operations along the x and y axes [24]. Second, the phase term in the space domain can be decomposed into a coupling operation (between x and y) followed by a chirp transform. Similarly, the phase term in the LCT domain can be decomposed into another coupling operation (between x0 and y0 ) and a 2D chirp; see Eq. (A2): I 0 b12 b22
B−1 A 0 B B B B B B B B @
b21 b22 detB , sx
I 1
0
0
1 p2 2 detB
1
0
1
0
0 0
0
1 k2 2 detB
p2 2 detB I 0
B B B B B B B B @
10
CB B 0 0C CB 0 CB p1 B 1 0C CB detB CB CB A@ 0 1 0
0
DB−1 0
0 0
b22 detB,
I
0 k2 2 detB
0
10
1
CB B 0 0C CB 0 CB k1 B 1 0C CB detB CB CB A@ 0 1 0
0
0
1
0
0
1
p3 0 detB
0
0
1
0
0
1
k3 0 detB
0
1
C 0C C C 0C C C C A 1
0
1
C 0C C C 0C C: C C A 1 (A2)
APPENDIX B In this appendix, the mathematical derivation of the sampling theorem for the detB ≠ 0 case developed in this paper is discussed. From Eq. (12), the center point of the zero-order term ˆ 0 ; y0 is defined by xˆ 0 ; yˆ 0 0;0 0; 0. From (n m 0) of Gx Fig. 2(b) we know that Lx0 is the distance between (0,0) and xˆ 0 ; yˆ 0 1;0 bT11x ; bT21x ; Ly0 is the distance between (0,0) and Fig. 5. One example of the 2D-NS-LCT of a sampled image: The original output LCT image parallelogram shape is indicated by a solid bold (blue) line, while its replicas are shown bounded by thin (black) lines. The parallelogram shapes outline the skewed sample grid. In this example all the sampling points are located in a regular Cartesian coordinate system. As a result all the color coded points (yellow, purple, light blue, and red) in the central parallelogram are found suitably arranged inside the new central square shape, bound by the thickdashed (red) line.
xˆ 0 ; yˆ 0 0;1 bT12y ; bT22y , which indicates that s 2 2 q b11 b21 1 b211 b221 ; −0 −0 Lx0 Tx Tx Tx s 2 2 q b12 b22 1 Ly0 b212 b222 −0 −0 Ty Ty Ty
(B1)
(B2)
Zhao et al.
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
2639
APPENDIX C In this appendix, the mathematical derivation of the sampling theorem for the B 0 case developed in this paper is discussed. The shape or boundary of the new output image can be defined by the corner co-ordinate,
x01
x02
x03
x04
y01
y02
y03
y04
1 2 detD
1 2 detD
d22
−d21
L
x
Lx
−Lx
−Lx
Ly −Ly −Ly Ly −d12 d11 d L −d L d22 Lx d21 Ly −d22 Lx d21 Ly 22 x 21 y −d12 Lx d11 Ly
−d12 Lx − d11 Ly
d12 Lx − d11 Ly
−d22 Lx − d21 Ly : d12 Lx d11 Ly
(C1)
Then, q d212 d222 ;
Lx0
q x02 − x03 2 y02 − y03 2
Lx detD
Ly 0
q x01 − x02 2 y01 − y02 2
Ly q d211 d221 : detD
(C2a)
(C2b)
APPENDIX D In this appendix, the mathematical proof regarding how the DLCT is made unitary is presented. The multiplication of W M and W M −1 is shown as W M −1 W M
N X N X
W M −1 n0 m0 kl W M klnm
k1 l1
jπp1 n2 − n02 L2x p2 nm − n0 m0 Lx Ly p3 m2 − m02 L2y 1 exp ; detB N 2 detB
8 n o n o9
> > kb lb kb lb > > 0 0 0 0 11 12 21 22 N X N > N detB > > k1 l1 : ;
jπp1 n2 − n02 L2x p2 nm − n0 m0 Lx Ly p3 m2 − m02 L2y 1 exp ; detB N 2 detB N N X −j2πkn − n0 −j2πkm − m0 X −j2πln − n0 −j2πlm−0 exp exp : exp exp N N N N k1 l1
It can be seen that when n n0 , m m0 , Eq. (D1) can be rearranged as W M −1 W M
N2 I; detB
(D1)
support of the Science Foundation Ireland, the Irish Research Council, and Enterprise Ireland under the National Development Plan.
(D2)
where I is a 4D identity matrix. Thus we modify Eq. (D1) by N dividing every element of the kernel by p , and in this padetB
per we treat V M as this normalized version of the kernel map W detB . This modification makes our trix W M , namely V M M N definition of the 2D-NS-DLCT unitary if we use the sampling specified by Eq. (15). Thus, after modification, V M −1 V M I, and hence V M −1 V −1 M.
ACKNOWLEDGMENTS L. Zhao is supported by a UCD—China Scholarship Council (CSC) joint scholarship and a SPIE Optics and Photonics Education Scholarship. J. J. Healy acknowledges the support of the National University of Ireland (NUI) Postdoctoral Fellowship in the Sciences. The authors also acknowledge the
REFERENCES 1. M. J. Bastiaans, “Wigner distribution function and its application to first-order optics,” J. Opt. Soc. Am. 69, 1710–1716 (1979). 2. S.-C. Pei, “Discrete fractional Fourier transform based on orthogonal projections,” IEEE Trans. Signal Process. 47, 1335–1348 (1999). 3. A. Stern, “Why is the linear canonical transform so little known?” AIP Conf. Proc. 860, 225–234 (2006). 4. B. Deng, R. Tao, and Y. Wang, “Convolution theorems for the linear canonical transform and their application,” Info. Sci. 49, 592–603 (2006). 5. A. Stern, “Uncertainty principles in linear canonical transform domains and some of their implications in optics,” J. Opt. Soc. Am. A 25, 647–652 (2008). 6. B. M. Hennelly and J. T. Sheridan, “Generalizing, optimizing, and inventing numerical algorithms for the fractional Fourier, Fresnel, and linear canonical transforms,” J. Opt. Soc. Am. A 22, 917–927 (2005).
2640 7. 8. 9. 10. 11. 12.
13. 14. 15. 16. 17.
18. 19. 20.
21. 22.
23. 24. 25.
26. 27. 28. 29. 30.
31.
J. Opt. Soc. Am. A / Vol. 31, No. 12 / December 2014 A. Stern, “Sampling of linear canonical transformed signals,” Signal Processing 86, 1421–1425 (2006). A. Koç, H. M. Ozaktas, and L. Hesselink, “Fast and accurate computation of two-dimensional non-separable quadratic-phase integrals,” J. Opt. Soc. Am. A 27, 1288–1302 (2010). T. Alieva and M. J. Bastiaans, “Properties of the linear canonical integral transformation,” J. Opt. Soc. Am. A 24, 3658–3665 (2007). R. Simon and K. B. Wolf, “Structure of the set of paraxial optical systems,” J. Opt. Soc. Am. A 17, 342–355 (2000). E. G. Abramochkin and V. G. Volostnikov, “Generalized Gaussian beams,” J. Opt. A 6, S157–S161 (2004). L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, “Orbital angular momentum of light and the transformation of Laguerre Gaussian laser modes,” Phys. Rev. A 45, 8185–8189 (1992). J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Experimental implementation of the gyrator transform,” J. Opt. Soc. Am. A 24, 3135–3139 (2007). J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Gyrator transform: properties and applications,” Opt. Express 15, 2190–2203 (2007). J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Applications of gyrator transform for image processing,” Opt. Commun. 278, 279–284 (2007). K. Wolf and T. Alieva, “Rotation and gyration of finite two dimensional modes,” J. Opt. Soc. Am. A 25, 365–370 (2008). S.-C. Pei and J.-J. Ding, “Properties, digital implementation, applications, and self image phenomena of the gyrator transform,” in 17th European Signal Processing Conference (EUSIPCO) (Curran Associates, Inc., 2009), pp. 441–444. J. Shamir, “Cylindrical lens systems described by operator algebra,” Appl. Opt. 18, 4195–4202 (1979). H. H. Arsenault, “A matrix representation for non-symmetrical optical systems,” J. Opt. 11, 87–91 (1980). G. Nemes and A. E. Siegman, “Measurement of all ten secondorder moments of an astigmatic beam by the use of rotating simple astigmatic (anamorphic) optics,” J. Opt. Soc. Am. A 11, 2257–2264 (1994). J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Optical system design for orthosymplectic transformations in phase space,” J. Opt. Soc. Am. A 23, 2494–2500 (2006). A. Sahin, H. M. Ozaktas, and D. Mendlovic, “Optical implementations of two-dimensional fractional Fourier transforms and linear canonical transforms with arbitrary parameters,” Appl. Opt. 37, 2130–2141 (1998). A. Sahin, M. A. Kutay, and H. M. Ozaktas, “Nonseparable two-dimensional fractional Fourier transform,” Appl. Opt. 37, 5444–5453 (1998). S.-C. Pei, “Two-dimensional affine generalized fractional Fourier transform,” IEEE Trans. Signal Process. 49, 878–897 (2001). P. Dong and N. P. Galatsanos, “Affine transformation resistant watermarking based on image normalization,” in IEEE International Conference on Image Processing (IEEE, 2002), Vol. 3, pp. 489–492. Z.-J. Liu, H. Chen, T. Liu, P.-F. Li, J.-M. Dai, X.-G. Sun, and S.-T. Liu, “Double-image encryption based on the affine transform and the gyrator transform,” J. Opt. 12, 035407 (2010). Y. Zhang, B.-Z. Dong, B.-Y. Gu, and G.-Z. Yang, “Beam shaping in the fractional Fourier transform domain,” J. Opt. Soc. Am. A 15, 1114–1120 (1998). B. M. Hennelly, D. Kelly, A. Corballis, and J. T. Sheridan, “Phase retrieval using theoretically unitary discrete fractional Fourier transform,” Proc. SPIE 5908, 59080D1 (2005). F.-C. Zhang, I. Yamaguchi, and L. P. Yaroslavsky, “Algorithm for reconstruction of digital holograms with adjustable magnification,” Opt. Lett. 29, 1668–1670 (2004). D.-Y. Wang, J. Zhao, F.-C. Zhang, G. Pedrini, and W. Osten, “High-fidelity numerical realization of multiple-step Fresnel propagation for the reconstruction of digital holograms,” Appl. Opt. 47, D12–D20 (2008). P. F. Almoro and S. G. Hanson, “Object wave reconstruction by speckle illumination and phase retrieval,” J. Eur. Opt. Soc. Rapid 4, 09002 (2009).
Zhao et al. 32. G.-H. Situ, J. P. Ryle, U. Gopinathan, and J. T. Sheridan, “Generalized in-line digital holographic technique based on intensity measurements at two different planes,” Appl. Opt. 47, 711–717 (2008). 33. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35, 237–246 (1972). 34. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). 35. K. B. Wolf, “Concepts from complex vector analysis and the Fourier transform,” in Integral Transforms in Science and Engineering (Plenum, 1979), Chap. 1, p. 16. 36. R. E. Ziemer, W. H. Tranter, and D. R. Fannin, Signals and Systems: Continuous and Discrete (Prentice-Hall, 1983). 37. R. Bracewell, “The discrete Fourier transform and the FFT,” in The Fourier Transform and its Applications (McGraw-Hill, 1999), Chap. 11, p. 261. 38. A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 2nd ed. (Prentice-Hall, 1999), Chaps. 8 and 9, pp. 541–692. 39. J. W. Goodman, Introduction to Fourier Optics (Roberts & Company, 2005). 40. U. Schnars and W. Jüptner, “Digital holography,” in Digital Holography: Digital Hologram Recording, Numerical Reconstruction, and Related Techniques (Springer, 2010), Chap. 3, pp. 46–48. 41. F. Gori, “Fresnel transform and sampling theorem,” Opt. Commun. 39, 293–297 (1981). 42. J.-J. Ding, “Research of fractional Fourier transform and linear canonical transform,” Ph.D. dissertation (National Taiwan University, 2001). 43. B. M. Hennelly and J. T. Sheridan, “Fast numerical algorithm for the linear canonical transform,” J. Opt. Soc. Am. A 22, 928–937 (2005). 44. D. P. Kelly, B. M. Hennelly, and C. McElhinney, “A practical guide to digital holography and generalized sampling,” Proc. SPIE 7072, 707215 (2008). 45. J. J. Healy, B. M. Hennelly, and J. T. Sheridan, “Additional sampling criterion for the linear canonical transform,” Opt. Lett. 33, 2599–2601 (2008). 46. J. J. Healy and J. T. Sheridan, “Sampling and discretization of the linear canonical transform,” Signal Processing 89, 641–648 (2009). 47. D. P. Kelly, B. M. Hennelly, N. Pandey, T. J. Naughton, and W. T. Rhodes, “Resolution limits in practical digital holographic systems,” Opt. Eng. 48, 095801 (2009). 48. J. J. Healy and J. T. Sheridan, “Fast linear canonical transforms,” J. Opt. Soc. Am. A 27, 21–29 (2010). 49. J. J. Healy and J. T. Sheridan, “Reevaluation of the direct method of calculating Fresnel and other linear canonical transforms,” Opt. Lett. 35, 947–949 (2010). 50. D. P. Kelly, J. J. Healy, B. M. Hennelly, and J. T. Sheridan, “Quantifying the 2.5 D imaging performance of digital holographic systems,” J. Eur. Opt. Soc. Rapid 6, 11034 (2011). 51. J. J. Healy and J. T. Sheridan, “Space-bandwidth ratio as a means of choosing between Fresnel and other linear canonical transform algorithms,” J. Opt. Soc. Am. A 28, 786–790 (2011). 52. J. T. Sheridan, L. Zhao, and J. J. Healy, “The condition number and first order optical systems,” in Imaging and Applied Optics Conference, OSA Technical Digest (OSA, 2012), paper ITu3C.4. 53. L. Zhao, J. J. Healy, and J. T. Sheridan, “Unitary discrete linear canonical transform: analysis and application,” Appl. Opt. 52, C30–C36 (2013). 54. J.-J. Ding and S.-C. Pei, “Eigenfunctions and self-imaging phenomena of the two-dimensional non-separable linear canonical transform,” J. Opt. Soc. Am. A 28, 82–95 (2011). 55. J.-J. Ding, S.-C. Pei, and C.-L. Liu, “Improved implementation algorithms of the two-dimensional non-separable linear canonical transform,” J. Opt. Soc. Am. A 29, 1615–1624 (2012). 56. T. Alieva and M. J. Bastiaans, “Alternative representation of the linear canonical integral transform,” Opt. Lett. 30, 3302–3304 (2005).
Zhao et al. 57. L. Zhao, J. J. Healy, and J. T. Sheridan, “Sampling of the two dimensional non-separable linear canonical transform,” Proc. SPIE 9131, 913112 (2014). 58. F. R. Gantmacher, “Matrices and operations on matrices,” in The Theory of Matrices (American Mathematical Society, 1959), Chap. 1, Vol. 1, pp. 1–7.
Vol. 31, No. 12 / December 2014 / J. Opt. Soc. Am. A
2641
59. F. R. Gantmacher, “Linear operators in an n-dimensional vector space,” in The Theory of Matrices (American Mathematical Society, 1959), Chap. 3, Vol. 1, pp. 55–56. 60. H.-J. Bao, A.-J. Sang, and H.-X. Chen, “Inverse operation of fourdimensional vector matrix,” Int. J. Intell. Syst. Appl. 3, 19–27 (2011).
