arXiv:0912.1379v1 [math.NA] 8 Dec 2009 A fast algorithm for the linear canonical transform Rafael G. Campos and Jared Figueroa Facultad de Ciencias Fı́sico-Matemáticas, Universidad Michoacana, 58060, Morelia, Mich., México. rcampos@umich.mx, jared@fismat.umich.mx MSC: 65T50; 44A15; 65D32 Keywords: Linear Canonical Transform, Fractional Fourier Transform, Quadrature, Hermite polynomials, Fractional Discrete Fourier Transform, fft Abstract In recent years there has been a renewed interest in finding fast algorithms to compute accurately the linear canonical transform (LCT) of a given function. This is driven by the large number of applications of the LCT in optics and signal processing. The well-known integral transforms: Fourier, fractional Fourier, bilateral Laplace and Fresnel transforms are special cases of the LCT. In this paper we obtain an O(N log N) algorithm to compute the LCT by using a chirp-FFT-chirp transformation yielded by a convergent quadrature formula for the fractional Fourier transform. This formula gives a unitary discrete LCT in closed form. In the case of the fractional Fourier transform the algorithm computes this transform for arbitrary complex values inside the unitary circle and not only at the boundary. In the case of the ordinary Fourier transform the algorithm improves the output of the FFT. 1 1 Introduction The Linear Canonical Transform (LCT) of a given function f (x) is a three-parameter integral transform that was obtained in connection with canonical transformations in Quantum Mechanics [1, 2]. It is defined by Z ∞ i 1 2 2 {a,b,c,d} e 2b (ax −2xy+dy ) f (x)dx, L [f (x), y] = √ 2πib −∞ √ i 2 for b 6= 0, and by de 2 cdy f (dy), if b = 0. The four parameters a, b, c and d appearing in (1), are the elements of a 2 × 2 matrix with unit determinant, i.e., ad − bc = 1. Therefore, only three parameters are free. Since this transform is a useful tool for signal processing and optical analysis, its study and direct computation in digital computers has become an important issue [3]-[10], particularly, fast algorithms to compute the linear canonical transform have been devised [4, 7]. These algorithms use the following related ideas: (a) use of the periodicity and shifting properties of the discrete LCT to break down the original matrix into smaller matrices as the FFT does with the DFT, (b) decomposition of the LCT into a chirp-FFT-scaling transformation and (c) decomposition of the LCT into a fractional Fourier transform followed by a scaling-chirp multiplication. All of these are algorithms of O(N log N) complexity. In this paper we present an algorithm that takes O(N log N) time based in the decomposition of the LCT into a scaling-chirp-DFT-chirp-scaling transformation, obtained by using a quadrature formula of the continuous Fourier transform [11, 12]. Here, DFT stands for the standard discrete Fourier transform. To distinguish this discretization from other implementations, we call it the extended Fourier Transform (XFT). Thus, the quadrature from which the XFT is obtained, uses some asymptotic properties of the Hermite polynomials and yields a fast algorithm to compute the Fourier transform, the fractional Fourier transform and therefore, the LCT. The quadrature formula is O(1/N)-convergent to the continuous Fourier transform for certain class of functions [13]. 2 A discrete fractional Fourier transform In previous work [12], [13], [14], we derived a quadrature formula for the continuous Fourier transform which yields an accurate discrete Fourier transform. For the sake of completeness we give in this section a brief review of the main steps to obtain this formula. Let us consider the family of Hermite polynomials Hn (x), n = 0, 1, . . ., which satisfies the recurrence equation Hn+1 (x) + 2nHn−1 (x) = 2xHn (x), (1) 2 with H−1 (x) ≡ 0. Note that the recurrence equation (1) can be written as the eigenvalue problem      0 1/2 0 · · · H0 (x) H0 (x) 1 0 1/2 · · · H1 (x) H1 (x)      (2) 0 2  H2 (x) = x H2 (x) . 0 · · ·      .. .. .. . . .. .. . . . . . . Let us now consider the eigenproblem associated to the principal submatrix of dimension N of (2)   0 1/2 0 · · · 0 0 1 0 1/2 · · · 0 0    0 2  0 · · · 0 0   H =  .. .. .. .. ..  . .. . . . . . .    0 0 0 ··· 0 1/2 0 0 0 ··· N −1 0 It is convenient to symmetrize H by using the similarity transformation SHS −1 where S is the diagonal matrix ) ( 1 1 . S = diag 1, √ , . . . , p 2 (N − 1)! 2N −1 Thus, the symmetric matrix H = SHS −1  q 1 0 0 2 q q  1 2  2 0 2  q  2  0 0 2   .. .. ..  . . .   0 0  0  0 0 0 takes the form ··· 0 ··· 0 ··· .. . 0 .. . ··· ··· 0 q N −1 2 0       0  . ..   q.  N −1  2   0 0 The recurrence equation (1) and the Christoffel-Darboux formula [15] can be used to solve the eigenproblem Huk = xk uk , k = 1, 2, . . . , N, which is a finite-dimensional version of (2). The eigenvalues xk are the zeros of HN (x) and the kth eigenvector uk is given by ck (s1 H0 (xk ), s2 H1 (xk ), · · · , sN HN −1 (xk ))T , where s1 , . . . , sN are the diagonal elements of S and ck is a normalization constant that can be determined from the condition uTk uk = 1, i.e., from c2k N −1 X n=0 Hn (xk )Hn (xk ) = 1. 2n n! 3 Therefore, r 2N −1 (N − 1)! (−1)N +k . N HN −1 (xk ) Thus, the components of the orthonormal vectors uk , k = 1, 2, . . . , N, are s 2N −n (N − 1)! Hn−1 (xk ) , (uk )n = (−1)N +k N (n − 1)! HN −1 (xk ) ck = (3) n = 1, . . . , N. Let U be the orthogonal matrix whose kth column is uk and let us define the matrix √ Fz = 2πU −1 D(z)U, where D(z) is the diagonal matrix D(z) = diag{1, z, z 2 , . . . , z N −1 } and z is an complex number. Therefore, the components of Fz are given by (Fz )jk = √ N −1 (−1)j+k 2N −1 (N − 1)! X z n Hn (xj )Hn (xk ). 2π N HN −1 (xj )HN −1 (xk ) n=0 2n n! (4) Next, we want to prove that if N is large enough, (4) approaches the kernel of the fractional Fourier transform evaluated at x = xj , y = xk . To this, we use the asymptotic expression for HN (x) [15]) 2 √ Γ(N + 1)ex /2 Nπ HN (x) ≃ cos( 2N + 1 x − ). (5) Γ(N/2 + 1) 2 Thus, the asymptotic form of the zeros of HN (x) are   2k − N − 1 π √ , xk = 2 2N (6) k = 1, 2, . . . , N. The use of (5) and (6) yields HN −1 (xk ) ≃ (−1)N +k Γ(N) x2k /2 , e Γ( N2+1 ) N → ∞, and the substitution of this asymptotic expression in (4) yields (Fz )jk ≃ √ 2π ∞ 2N −1 [Γ( N 2+1 )]2 −(x2 +x2 )/2 X z n e j k H (x )Hn (xk ). n n! n j Γ(N + 1) 2 n=0 Finally, Stirling’s formula and Mehler’s formula [16] produce r 2 (1+z 2 )(x2 j +xk )−4xj xk z 2 − 2(1−z 2 ) ∆xk , (Fz )jk ≃ e 1 − z2 (7) where ∆xk is the difference between two consecutive asymptotic Hermite zeros, i.e., π ∆xk = xk+1 − xk = √ . (8) 2N 4 Let us consider now the vector of samples of a given function f (x) f = (f (x1 ), f (x2 ), . . . , f (xN ))T . The multiplication of the matrix Fz by the vector f gives the vector g with entries r 2 2 +x2 )−4x x z N N X j k j k 2 X − (1+z )(x2(1−z 2) f (xk )∆xk , gj = (Fz )jk f (xk ) ≃ e 1 − z 2 k=1 k=1 where j = 1, 2 . . . , N. This equation is a Riemann sum for the integral r Z ∞ (1+z2 )(y2 +x2 )−4xyz 2 − 2(1−z 2 ) f (x)dx, Fz [f (x), y] = e 2 1 − z −∞ where |z| < 1. Therefore, if we make yj = xj , Fz [f (x), yj ] ≃ N X (Fz )jk f (xk ), k=1 N → ∞. (9) Note that Fz [g(x), y] is the continuous fractional Fourier transform [17] of g(x) except for a constant and therefore, Fz is a discrete fractional Fourier transform. 3 A fast linear canonical transform Firstly, note that if b 6= 0, the LCT can be written as a chirp-FT-chirp transform e L{a,b,c,d}[f (x), y] = √ idy 2 2b 2πib Z ∞ e− ixy b e iax2 2b f (x)dx. −∞ Thus, for b 6= 0, the LCT of the function f (x) can be represented by the 1/b-scaled Fourier idy 2 iax2 e 2b transform of the function g(x) = e 2b f (x), multiplied by √ . 2πib On the other hand, note that for the case z = ±i, (7) yields a discrete Fourier transform (F±i )jk ≃ e±itj tk ∆tk , that can be related to the standard DFT as follows. The use of (6) yields N−1 N−1 π2 π (10) (Fi )jk = e±itj tk ∆tk = √ ei 2N (j− 2 )(k− 2 ) 2N P where we have used (6) and (8). Since N k=1 (Fi )jk f (xk ) is a quadrature and therefore, an approximation of Z g(yj ) = a scaled Fourier transform ∞ eiyj x f (x)dx, −∞ Z ∞ eiκyj x f (x)dt = g(κyj ) −∞ 5 (11) has the quadrature PN k=1 F̃jk f (xk ), where N−1 N−1 π2 π F̃jk = √ eiκ 2N (j− 2 )(k− 2 ) . 2N (12) If we choose κ = 4/π, (12) takes the form Fjk 2 π (N−1) N πei 2 √ = i h 2π i h i h N−1 N−1 e−iπ N j ei N jk e−iπ N k , 2N PN of g(4yj /π). If for j, k = 0, 1, 2, . . . , N − 1, and k=1 Fjk f (xk ) is an approximation PN now we choose κ = 4b/π, but we keep the same matrix (13), then k=1 Fjk f (xk ) is an approximation of Z ∞ yj ei b x f (x)dt. −∞ If now we replace f (x) by e iax2 2b f (x) and take into account (10), we have that N X 2 Fjk eiaxk /2b f (xk ) k=1 is an approximation of the product of functions  idy 2 e 2b √ 2πib −1 L{a,b,c,d} [f (x), y] evaluated at yj = 4bxj /π. Therefore, a discrete (scaled) linear canonical transform L can be given in closed form. If we denote by G(y) the LCT of f (x), then G(yj ) = G(4bxj /π) = N X (S1 F S2 )jk f (xk ), k=1 idyj2 √ where S1 and S2 are diagonal matrices whose diagonal elements are e 2b / 2πib, and 2 eiaxj /2b , respectively. As it can be seen, the matrix L = S1 F S2 , which gives the discrete LTC, i.e., the XFT, consists in a chirp-DFT-chirp transformation, where DFT stands for the standard discrete Fourier transform. Therefore, we can use any FFT to give a fast computation of the linear canonical transform G = Lf . Now, the fast algorithm for the linear canonical transform can be given straightforwardly. 6 Algorithm To compute an approximation G = (G1 , G2 , . . . , GN )T of the linear canonical  trans2j−N −1 π {a,b,c,d} √ form G(y) = L [f (x), y] evaluated at yj = 4bxj /π, where xj = . 2 2N 1. For given N set up the vector u of components   2k − N − 1 2 /2b iax −iπ (k−1)(N−1) N √ e k f π , uk = e 2 2N k = 1, 2, . . . , N. 2. Set yj = 4bxj /π and compute the diagonal matrix S according to Sjk 2 π (N−1) N πei 2 = √ e idyj2 2b 2N N−1 e−iπ N √ 2πib (j−1) δjk , j, k = 1, 2, . . . , N. 2π 3. Let DF be the discrete Fourier transform, i.e., (DF )jk = ei N jk , j, k = 0, 1, 2, . . . , N − 1. Obtain the approximation Gj to G( 4b x ) by computing π j the matrix-vector product G = SDF u, (13) with a standard FFT algorithm. 4 Example For this example we take an integral formula found in [18] that gives √ π G(y) = √ 2 2πib(α + i × e acy 2 2(4b2 α2 +a2 ) i e a2 1/4 ) 4b2 e α(β 2 −αγ) 2 α2 + a 2 4b 2b(α2 dy 2 +2βαy+β 2 a) 4b2 α2 +a2 2 i a − αy e 2 arctan( 2αb ) e 2 +2βay+a2 γ 4b2 α2 +a2 (14) , if f (x) = e−(αy +2βy+γ) , α > 0. Figure 1 shows the exact LTC with α = 1, β = 2, γ = 3, a = 1, b = 2, c = 1/2, and d = 2, compared with the approximation given by the XFT. 2 Figure 2 shows the exact Fresnel transform of f (x) = e−(αy +2βy+γ) with α = 2, β = 1, γ = 3, a = 1, b = 100, c = 0, and d = 1, compared with the approximation given by the XFT. 7 B 1.0 1.0 0.5 0.5 GH yL GH yL A 0.0 -0.5 0.0 -0.5 -1.0 -1.0 -6 -4 -2 0 y 2 4 6 -6 -4 -2 0 y 2 4 6 Figure 1: Real part (A) and imaginary part (B) of the exact linear canonical transform (solid line) compared with the output of the XFT (dashed line) computed with N = 512. The function G(y) is that given in (14) for α = 1, β = 2, γ = 3, a = 1, b = 2, c = 1/2, and d = 2. B 0.004 0.002 0.002 GH yL GH yL A 0.004 0.000 0.000 -0.002 -0.002 -0.004 -4 -2 0 y 2 -0.004 4 -4 -2 0 y 2 4 Figure 2: Real part (A) and imaginary part (B) of the exact Fresnel transform (solid line) compared with the output of the XFT (dashed line) computed with N = 1024. The function G(y) is that given in (14) for α = 2, β = 1, γ = 3, a = 1, b = 100, c = 0, and d = 1. 8 5 Conclusion We have obtained a discrete linear canonical transform and a fast algorithm to compute this transform by projecting the space of functions onto a vector space spanned by a finite number of Hermite functions. The XFT is a discrete LCT given by a unitary matrix in a closed form in which the DFT can be found at the core, surrounded by diagonal transformations, which makes easy to implement it in a fast algorithm. Since this discrete LCT is related to a quadrature formula of the fractional Fourier transform, it yields accurate results. References [1] M. Moshinsky and C. Quesne, “Linear Canonical Transformations and their unitary representations,” J. Math. Phys., vol. 12, pp. 1772-1783, 1971. [2] K.B Wolf, Integral Transforms in Science and Engineering, Ch. 9-10 Plenum Press, 1979. New York: [3] J.J.Healy, J.T.Sheridan, “Sampling and discretization of the linear canonical transform,” Signal Processing, vol. 89, pp. 641-648, 2009. [4] A. Koç, H.M. Ozaktas, C. Candan, M.A. Kutay, “Digital Computation of Linear Canonical Transforms,” IEEE Trans. Sig. Proc., vol. 56, pp. 2383-2394, 2008. [5] K.K.Sharma, S.D. Joshi, “Uncertainty Principle for Real Signals in the Linear Canonical Transform Domains,” IEEE Trans. Sig. Proc., vol. 56, pp. 2677-2683, 2008. [6] A. Stern, “Sampling of linear canonical transformed signals,” Signal Processing, vol. 86, pp. 1421-1425, 2006. [7] B.M. Hennelly, J.T. Sheridan, “Fast numerical algorithm for the linear canonical transform,” J. Opt. Soc. Am., vol. A22, pp. 928-937, 2005. [8] B.M. Hennelly, J.T. Sheridan, “Generalizing, optimizing, and inventing numerical algorithms for the fractional Fourier, Fresnel, and linear canonical transforms,” J. Opt. Soc. Am., vol. A22, pp. 917-927, 2005. [9] B.Z. Li, R. Tao, Y. Wang, “New sampling formulae related to linear canonical transform,” Signal Processing, vol. 87, pp. 983-990, 2007. [10] S.C. Pei, J.J. Ding, “Eigenfunctions of Linear Canonical Transform,” IEEE Trans. Sig. Proc., vol. 50, pp. 11-26, 2002. [11] R.G. Campos, J. Rico-Melgoza, E. Chávez, “RFT: A fast discrete fractional Fourier transform,” unpublished. 9 [12] R. G. Campos, and L.Z. Juárez, “A discretization of the Continuous Fourier Transform”, Il Nuovo Cimento, vol. 107B, pp. 703-711, 1992. [13] R.G. Campos, “A Quadrature Formula for the Hankel Transform,” Numerical Algorithms, vol. 9, pp. 343-354, 1995. [14] R. G. Campos, F. Domı́nguez Mota and E. Coronado, “Quadrature formulas for integrals transforms generated by orthogonal polynomials,” arXiv.0805.2111v1 [math.NA]. [15] G. Szego, Orthogonal Polynomials. Providence, Rhode Island: Colloquium Publications, American Mathematical Society, 1975. [16] A. Erdélyi, Higher Transcendental Functions. Vols. I and II. Hill, 1953. New York: McGraw [17] V. Namias, “The Fractional Order Fourier Transform and its Application to Quantum Mechanics,” J. Inst. Maths. Applics., vol. 25, pp. 241–265, 1980. [18] I.S. Gradshteyn, I.M. Ryzhik, Table of Integrals, Series and Products, Fifth edition New York: Academic Press, 1994, p. 520. 10