- Email: [email protected]
A fast algorithm for the linear canonical transform
Signal Processing 91 (2011) 1444–1447
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
A fast algorithm for the linear canonical transform Rafael G. Campos , Jared Figueroa ´ticas, Universidad Michoacana, 58060 Morelia, Michoaca ´n, Me´xico Facultad de Ciencias Fı´sico-Matema
a r t i c l e i n f o
abstract
Article history: Received 19 January 2010 Received in revised form 1 July 2010 Accepted 9 July 2010 Available online 15 July 2010
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ðNlogNÞ 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. This chirp-FFT-chirp transform approximates the ordinary Fourier transform more precisely than just the FFT, since it comes from a convergent procedure for non-periodic functions. & 2010 Elsevier B.V. All rights reserved.
Keywords: Linear canonical transform Fractional Fourier transform Quadrature Hermite polynomials Fractional discrete Fourier transform fft
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 1 2 1 2 Lfa,b,c,dg ½f ðxÞ,y ¼ pffiffiffiffiffiffiffiffiffiffi ei=2bðax 2xy þ dy Þ f ðxÞ dx, 2pib 1 pffiffiffi 2 for ba0, and by deði=2Þcdy f ðdyÞ, if b=0. The four parameters a, b, c and d appearing in the above expression, 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 have 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ðNlog NÞ complexity. In this paper we present an algorithm that takes OðNlog 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
Corresponding author.
E-mail addresses: [email protected] (R.G. Campos), jared@fismat.umich.mx (J. Figueroa). 0165-1684/$ - see front matter & 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2010.07.007
In previous work [12–14], we derived a quadrature formula for the continuous Fourier transform which yields
R.G. Campos, J. Figueroa / Signal Processing 91 (2011) 1444–1447
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,y, which satisfies the recurrence equation: Hn þ 1 ðxÞ þ 2nHn1 ðxÞ ¼ 2xHn ðxÞ,
ð1Þ
with H1 ðxÞ 0. Note that the recurrence equation (1) can be written as the eigenvalue problem 1 0 1 0 10 H0 ðxÞ 0 1=2 0 H0 ðxÞ C B C B CB 0 1=2 CB H1 ðxÞ C B H ðxÞ C B1 C ¼ x B 1 C: B CB ð2Þ B B H2 ðxÞ C C B0 2 0 C @ A @ A@ H2 ðxÞ A ^ ^ ^ ^ ^ & Let us now consider the eigenproblem associated to the principal submatrix of dimension N of (2) 0 1 0 1=2 0 0 0 B1 0 1=2 0 0 C B C B C B0 2 0 0 0 C C: H¼B B^ ^ ^ & ^ ^ C B C B C @0 0 0 0 1=2 A 0
0
0
N1
0
It is convenient to symmetrize H by using the similarity transformation SHS1 where S is the diagonal matrix 9 8 > > = < 1 1 S ¼ diag 1, pffiffiffi, . . . , qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : > > N1 2 ; : ðN1Þ!2 Thus, the symmetric matrix H ¼ SHS1 takes the form 0 1 rffiffiffi 1 0 0 0 C B 0 2 B C rffiffiffi B rffiffiffi C B 1 C 2 B 0 0 0 C B 2 C 2 B C rffiffiffi B C B C 2 B 0 C 0 0 0 B C: 2 B C B ^ C ^ ^ & ^ ^ B C r ffiffiffiffiffiffiffiffiffiffi ffi B C B N1 C B 0 C 0 0 0 B 2 C B C ffiffiffiffiffiffiffiffiffiffi ffi r B C N1 @ A 0 0 0 0 2 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
1445
Therefore, sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2N1 ðN1Þ! ð1ÞN þ k ck ¼ : N HN1 ðxk Þ Thus, the components of the orthonormal vectors uk, k= 1,2,y, N, are sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2Nn ðN1Þ! Hn1 ðxk Þ N þk , ð3Þ ðuk Þn ¼ ð1Þ Nðn1Þ! HN1 ðxk Þ n =1,y,N. Let U be the orthogonal matrix whose kth column is uk and let us define the matrix pffiffiffiffiffiffi F z ¼ 2pU 1 DðzÞU, where D(z) is the diagonal matrix D(z)= diag{1,z,z2,y,zN 1} and z is a complex number. Therefore, the components of F z are given by ðF z Þjk ¼
1 pffiffiffiffiffiffi ð1Þj þ k 2N1 ðN1Þ! NX zn Hn ðxj ÞHn ðxk Þ: 2p NHN1 ðxj ÞHN1 ð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] HN ðxÞ C
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 GðN þ 1Þex =2 Np cos 2N þ 1 x : GðN=2 þ 1Þ 2
Thus, the asymptotic form of the zeros of HN(x) are 2kN1 p pffiffiffiffiffiffiffi , xk ¼ 2 2N
ð5Þ
ð6Þ
k= 1,2,y,N. The use of (5) and (6) yields
GðNÞ
HN1 ðxk Þ C ð1ÞN þ k
N þ1 G 2
2
exk =2 ,
N-1,
and the substitution of this asymptotic expression in (4) yields Nþ1 2 N1 2 G pffiffiffiffiffiffi 2 ðx2 þ x2 Þ=2 e j k ðF z Þjk C 2p GðN þ 1Þ 1 X zn Hn ðxj ÞHn ðxk Þ: n 2 n! n¼0 Finally, Stirling’s formula and Mehler’s formula [16] produce " # rffiffiffiffiffiffiffiffiffiffiffi ð1þ z2 Þðx2j þx2k Þ4xj xk z 2 Dxk , ð7Þ exp ðF z Þjk C 1z2 2ð1z2 Þ
ck ðs1 H0 ðxk Þ,s2 H1 ðxk Þ, . . . ,sN HN1 ðxk ÞÞT ,
where Dxk is the difference between two consecutive asymptotic Hermite zeros, i.e.,
where s1,y,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
Dxk ¼ xk þ 1 xk ¼ pffiffiffiffiffiffiffi :
N 1 X
Hn ðxk ÞHn ðxk Þ ¼ 1: ck2 2n n! n¼0
p
2N
ð8Þ
Let us consider now the vector of samples of a given function f(x) f ¼ ðf ðx1 Þ,f ðx2 Þ, . . . ,f ðxN ÞÞT :
1446
R.G. Campos, J. Figueroa / Signal Processing 91 (2011) 1444–1447
The multiplication of the matrix F z by the vector f gives the vector g with entries gj ¼
N X
If we choose k ¼ 4=p, (12) takes the form 2
Fjk ¼
peiðp=2ÞððN1Þ
ðF z Þjk f ðxk Þ
PN
where j = 1,2y,N. This equation is a Riemann sum for the integral rffiffiffiffiffiffiffiffiffiffiffi Z 1 2 ð1þ z2 Þðy2 þ x2 Þ4xyz F z ½f ðxÞ,y ¼ f ðxÞ dx, exp 1z2 1 2ð1z2 Þ where jzj o1. Therefore, if we make yj = xj, F z ½f ðxÞ,yj C
ðF z Þjk f ðxk Þ,
N-1:
ð9Þ
Note that F z ½gðxÞ,y is the continuous fractional Fourier transform [17] of g(x) except for a constant and therefore, F z is a discrete fractional Fourier transform.
!
Lfa,b,c,dg ½f ðxÞ,y ¼
Z
! ixy iax2 exp exp f ðxÞ dx: b 2b 1 1
Thus, for ba0, the LCT of the function f(x) can be represented by the 1/b-scaled Fourier transform of the multiplied by function gðxÞ ¼ expðiax2 =2bÞf ðxÞ, pffiffiffiffiffiffiffiffiffiffi 2 expðidy =2bÞ= 2pib. On the other hand, note that for the case z ¼ 7i, (7) yields a discrete Fourier transform ðF 7 i Þjk C e 7 itj tk Dtk , that can be related to the standard DFT as follows. The use of (6) yields
k¼1
is an approximation of the product of functions pffiffiffiffiffiffiffiffiffiffi 2 evaluated at ðexpðidy =2bÞ= 2pibÞ1 Lfa,b,c,dg ½f ðxÞ,y yj ¼ 4bxj =p. 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
One
gðyj Þ ¼
Two iyj x
e
f ðxÞ dx,
Three
1
1
G ¼ SDF u,
eikyj x f ðxÞ dt ¼ gðkyj Þ,
ð11Þ
1
has the quadrature
PN
k¼1
transform GðyÞ ¼ Lfa,b,c,dg ½f ðxÞ,y evaluated at yj ¼ 4bxj =p. Set up the vectors u, y and s for k=1,2,y,N ! 2 ðk1ÞðN1Þ iaxk 2kN1 f p pffiffiffiffiffiffiffi uk ’exp ip exp N 2b 2 2N pffiffiffiffiffiffiffiffiffi yk ’4bxk =p ¼ 2=N bð2kN1Þ " # pffiffiffiffi p ðN1Þ2 idy2k N1 pexp i ðk1Þ þ ip N 2 N 2b pffiffiffiffiffiffiffiffi sk ’ 2 ibN end for Set up the diagonal matrix S S= diag(s1,s2,y,sN) Let DF be the discrete Fourier transform, ðDF Þjk ¼ expði 2Np jkÞ, j,k= 0,1,2,y,N 1. Obtain the approximation Gj to Gð4b p xj Þ by computing the matrixvector product
a scaled Fourier transform Z
ðS1 FS2 Þjk f ðxk Þ,
Algorithm XFT pffiffiffiffiffiffiffi Input: A list of values f ðpð2jN1Þ=2 2N Þ of a function f(x) at the pffiffiffiffiffiffiffi points xj ¼ ½jðN þ 1Þ=2p= 2N . Output: A list of values Gj approximating the linear canonical
P where we have used (6) and (8). Since N k ¼ 1 ðF i Þjk f ðxk Þ is a quadrature and therefore, an approximation of 1
N X
where S1 and S2 are diagonal matrices whose diagonal pffiffiffiffiffiffiffiffiffiffi 2 2 elements are eidyj =2b = 2pib, and eiaxj =2b , respectively. As it can be seen, the matrix L= S1FS2, 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 XFT algorithm for the linear canonical transform can be given straightforwardly.
2 p p N1 N1 k , ðF i Þjk ¼ e 7 itj tk Dtk ¼ pffiffiffiffiffiffiffi exp i j 2 2 2N 2N ð10Þ
Z
2
Fjk eiaxk =2b f ðxk Þ
k¼1
Firstly, note that if ba0, the LCT can be written as a chirp-FT-chirp transform 2
1
If now we replace f(x) by expðiax2 =2bÞf ðxÞ and take into account (10), we have that
Gðyj Þ ¼ Gð4bxj =pÞ ¼
3. A fast linear canonical transform
idy 2b pffiffiffiffiffiffiffiffiffiffi 2pib
for j,k=0,1,2,y,N 1, and k ¼ 1 Fjk f ðxk Þ is an approximation of gð4yj =pÞ. If now we choose k ¼ 4b=p, but we keep PN the same matrix (13), then k ¼ 1 Fjk f ðxk Þ is an approximation of Z 1 eiðyj =bÞx f ðxÞ dt:
N X
k¼1
exp
½eipððN1Þ=NÞj ½eið2p=NÞjk ½eipððN1Þ=NÞk , ð13Þ
k¼1 " # rffiffiffiffiffiffiffiffiffiffiffi N ð1 þ z2 Þðx2j þ x2k Þ4xj xk z 2 X f ðxk ÞDxk , C exp 1z2 k ¼ 1 2ð1z2 Þ
N X
=NÞ
pffiffiffiffiffiffiffi 2N
F~ jk f ðxk Þ, where
p p2 N1 N1 k : j F~ jk ¼ pffiffiffiffiffiffiffi exp ik 2 2 2N 2N
with a standard FFT algorithm.
4. An application to edge detection ð12Þ
As an untested and unexplored application of the algorithm given before, we use the XFT to find the edges
R.G. Campos, J. Figueroa / Signal Processing 91 (2011) 1444–1447
1447
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 comes from a quadrature formula, it yields accurate results. Fig. 1. (a) Original image. (b) Contoured image obtained by thresholding the two-dimensional XFT of the original image with parameters a1 = 1, d1 = 1, b1 =20, a2 = 1, d2 =1, b2 = 20, and threshold value Gt = 20.
Fig. 2. (c) Contoured image obtained by thresholding the two-dimensional XFT of (a) with parameters a1 = 0.2, d1 =10, b1 = 2.5, a2 = 0.2, d2 =10, b2 = 2, and threshold value Gt = 60. (d) Same parameters as in (c) but threshold value Gt = 100.
in an image. The problem of edge detection has been studied extensively since many years ago and many algorithms have been devised [18–21]. Our purpose in this section is only to point out a possible use of the XFT as an edge detector. Since an edge in an image is essentially a boundary between regions reflecting different amounts of energy, a cut-off to zero in the power spectrum of the image for values greater than a threshold value, let us say Gt, will define the boundary in the spatial domain. This simple idea is applied to the twodimensional XFT of a grayscale image with parameters a1, b1 and d1 for one dimension and a2, b2 and d2 for the other. The result of this procedure is shown in Figs. 1 and 2. In order to compare and contrast the processed images, the intensity of Figs. 1(b) and 2(c) and (d) has been doubled. As expected, different values of the parameters yield different contoured images and high values of the threshold parameter yield imprecise edges.
References [1] M. Moshinsky, C. Quesne, Linear canonical transformations and their unitary representations, J. Math. Phys. 12 (1971) 1772–1783. [2] K.B. Wolf, Integral Transforms in Science and Engineering, Plenum Press, New York, 1979 (Chapter 9–10). [3] J.J. Healy, J.T. Sheridan, Sampling and discretization of the linear canonical transform, Signal Process. 89 (2009) 641–648. [4] A. Koc- , H.M. Ozaktas, C. Candan, M.A. Kutay, Digital computation of linear canonical transforms, IEEE Trans. Signal Process. 56 (2008) 2383–2394. [5] K.K. Sharma, S.D. Joshi, Uncertainty principle for real signals in the linear canonical transform domains, IEEE Trans. Signal Process. 56 (2008) 2677–2683. [6] A. Stern, Sampling of linear canonical transformed signals, Signal Process. 86 (2006) 1421–1425. [7] B.M. Hennelly, J.T. Sheridan, Fast numerical algorithm for the linear canonical transform, J. Opt. Soc. Am. A 22 (2005) 928–937. [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. A 22 (2005) 917–927. [9] B.Z. Li, R. Tao, Y. Wang, New sampling formulae related to linear canonical transform, Signal Process. 87 (2007) 983–990. [10] S.C. Pei, J.J. Ding, Eigenfunctions of linear canonical transform, IEEE Trans. Signal Process. 50 (2002) 11–26. [11] R.G. Campos, J. Rico-Melgoza, E. Cha´vez, XFT: extending the digital application of the Fourier transform, arXiv:0911.0952v1, /http:// arxiv.org/abs/0911.0952v1S; 2009. [12] R.G. Campos, L.Z. Jua´rez, A discretization of the continuous Fourier transform, Il Nuovo Cimento 107B (1992) 703–711. [13] R.G. Campos, A quadrature formula for the Hankel transform, Numer. Algorithms 9 (1995) 343–354. [14] R.G. Campos, F. Domı´nguez Mota, E. Coronado, Quadrature formulas for integrals transforms generated by orthogonal polynomials, IMA J. Numer. Anal., to appear (cf. arXiv.0805.2111v1, 2008, /http://arxiv.org/abs/0805.2111S). [15] G. Szego, Orthogonal Polynomials, Colloquium Publications, American Mathematical Society, Providence, Rhode Island, 1975. [16] A. Erde´lyi, Higher Transcendental Functions, vols. I and II, McGraw-Hill, New York, 1953. [17] V. Namias, The fractional order Fourier transform and its application to quantum mechanics, J. Inst. Math. Appl. 25 (1980) 241–265. [18] D. Ziou, S. Tabbone, Edge detection techniques: an overview, Int. J. Pattern Recognit. Image Anal. 8 (1998) 537–559. [19] D. Marr, E. Hildreth, Theory of edge detection, Proc. R. Soc. London 207 (1980) 187–217. [20] W.Y. Ma, B.S. Manjunath, Edge flow: a framework of boundary detection and image segmentation, in: IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’97), 1997, p. 744. [21] Y. Deng, B.S. Manjunath, H. Shin, Color image segmentation, in: IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’99), vol. 2, 1999, p. 2446.
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
A fast algorithm for the linear canonical transform Rafael G. Campos , Jared Figueroa ´ticas, Universidad Michoacana, 58060 Morelia, Michoaca ´n, Me´xico Facultad de Ciencias Fı´sico-Matema
a r t i c l e i n f o
abstract
Article history: Received 19 January 2010 Received in revised form 1 July 2010 Accepted 9 July 2010 Available online 15 July 2010
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ðNlogNÞ 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. This chirp-FFT-chirp transform approximates the ordinary Fourier transform more precisely than just the FFT, since it comes from a convergent procedure for non-periodic functions. & 2010 Elsevier B.V. All rights reserved.
Keywords: Linear canonical transform Fractional Fourier transform Quadrature Hermite polynomials Fractional discrete Fourier transform fft
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 1 2 1 2 Lfa,b,c,dg ½f ðxÞ,y ¼ pffiffiffiffiffiffiffiffiffiffi ei=2bðax 2xy þ dy Þ f ðxÞ dx, 2pib 1 pffiffiffi 2 for ba0, and by deði=2Þcdy f ðdyÞ, if b=0. The four parameters a, b, c and d appearing in the above expression, 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 have 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ðNlog NÞ complexity. In this paper we present an algorithm that takes OðNlog 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
Corresponding author.
E-mail addresses: [email protected] (R.G. Campos), jared@fismat.umich.mx (J. Figueroa). 0165-1684/$ - see front matter & 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2010.07.007
In previous work [12–14], we derived a quadrature formula for the continuous Fourier transform which yields
R.G. Campos, J. Figueroa / Signal Processing 91 (2011) 1444–1447
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,y, which satisfies the recurrence equation: Hn þ 1 ðxÞ þ 2nHn1 ðxÞ ¼ 2xHn ðxÞ,
ð1Þ
with H1 ðxÞ 0. Note that the recurrence equation (1) can be written as the eigenvalue problem 1 0 1 0 10 H0 ðxÞ 0 1=2 0 H0 ðxÞ C B C B CB 0 1=2 CB H1 ðxÞ C B H ðxÞ C B1 C ¼ x B 1 C: B CB ð2Þ B B H2 ðxÞ C C B0 2 0 C @ A @ A@ H2 ðxÞ A ^ ^ ^ ^ ^ & Let us now consider the eigenproblem associated to the principal submatrix of dimension N of (2) 0 1 0 1=2 0 0 0 B1 0 1=2 0 0 C B C B C B0 2 0 0 0 C C: H¼B B^ ^ ^ & ^ ^ C B C B C @0 0 0 0 1=2 A 0
0
0
N1
0
It is convenient to symmetrize H by using the similarity transformation SHS1 where S is the diagonal matrix 9 8 > > = < 1 1 S ¼ diag 1, pffiffiffi, . . . , qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : > > N1 2 ; : ðN1Þ!2 Thus, the symmetric matrix H ¼ SHS1 takes the form 0 1 rffiffiffi 1 0 0 0 C B 0 2 B C rffiffiffi B rffiffiffi C B 1 C 2 B 0 0 0 C B 2 C 2 B C rffiffiffi B C B C 2 B 0 C 0 0 0 B C: 2 B C B ^ C ^ ^ & ^ ^ B C r ffiffiffiffiffiffiffiffiffiffi ffi B C B N1 C B 0 C 0 0 0 B 2 C B C ffiffiffiffiffiffiffiffiffiffi ffi r B C N1 @ A 0 0 0 0 2 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
1445
Therefore, sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2N1 ðN1Þ! ð1ÞN þ k ck ¼ : N HN1 ðxk Þ Thus, the components of the orthonormal vectors uk, k= 1,2,y, N, are sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2Nn ðN1Þ! Hn1 ðxk Þ N þk , ð3Þ ðuk Þn ¼ ð1Þ Nðn1Þ! HN1 ðxk Þ n =1,y,N. Let U be the orthogonal matrix whose kth column is uk and let us define the matrix pffiffiffiffiffiffi F z ¼ 2pU 1 DðzÞU, where D(z) is the diagonal matrix D(z)= diag{1,z,z2,y,zN 1} and z is a complex number. Therefore, the components of F z are given by ðF z Þjk ¼
1 pffiffiffiffiffiffi ð1Þj þ k 2N1 ðN1Þ! NX zn Hn ðxj ÞHn ðxk Þ: 2p NHN1 ðxj ÞHN1 ð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] HN ðxÞ C
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 GðN þ 1Þex =2 Np cos 2N þ 1 x : GðN=2 þ 1Þ 2
Thus, the asymptotic form of the zeros of HN(x) are 2kN1 p pffiffiffiffiffiffiffi , xk ¼ 2 2N
ð5Þ
ð6Þ
k= 1,2,y,N. The use of (5) and (6) yields
GðNÞ
HN1 ðxk Þ C ð1ÞN þ k
N þ1 G 2
2
exk =2 ,
N-1,
and the substitution of this asymptotic expression in (4) yields Nþ1 2 N1 2 G pffiffiffiffiffiffi 2 ðx2 þ x2 Þ=2 e j k ðF z Þjk C 2p GðN þ 1Þ 1 X zn Hn ðxj ÞHn ðxk Þ: n 2 n! n¼0 Finally, Stirling’s formula and Mehler’s formula [16] produce " # rffiffiffiffiffiffiffiffiffiffiffi ð1þ z2 Þðx2j þx2k Þ4xj xk z 2 Dxk , ð7Þ exp ðF z Þjk C 1z2 2ð1z2 Þ
ck ðs1 H0 ðxk Þ,s2 H1 ðxk Þ, . . . ,sN HN1 ðxk ÞÞT ,
where Dxk is the difference between two consecutive asymptotic Hermite zeros, i.e.,
where s1,y,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
Dxk ¼ xk þ 1 xk ¼ pffiffiffiffiffiffiffi :
N 1 X
Hn ðxk ÞHn ðxk Þ ¼ 1: ck2 2n n! n¼0
p
2N
ð8Þ
Let us consider now the vector of samples of a given function f(x) f ¼ ðf ðx1 Þ,f ðx2 Þ, . . . ,f ðxN ÞÞT :
1446
R.G. Campos, J. Figueroa / Signal Processing 91 (2011) 1444–1447
The multiplication of the matrix F z by the vector f gives the vector g with entries gj ¼
N X
If we choose k ¼ 4=p, (12) takes the form 2
Fjk ¼
peiðp=2ÞððN1Þ
ðF z Þjk f ðxk Þ
PN
where j = 1,2y,N. This equation is a Riemann sum for the integral rffiffiffiffiffiffiffiffiffiffiffi Z 1 2 ð1þ z2 Þðy2 þ x2 Þ4xyz F z ½f ðxÞ,y ¼ f ðxÞ dx, exp 1z2 1 2ð1z2 Þ where jzj o1. Therefore, if we make yj = xj, F z ½f ðxÞ,yj C
ðF z Þjk f ðxk Þ,
N-1:
ð9Þ
Note that F z ½gðxÞ,y is the continuous fractional Fourier transform [17] of g(x) except for a constant and therefore, F z is a discrete fractional Fourier transform.
!
Lfa,b,c,dg ½f ðxÞ,y ¼
Z
! ixy iax2 exp exp f ðxÞ dx: b 2b 1 1
Thus, for ba0, the LCT of the function f(x) can be represented by the 1/b-scaled Fourier transform of the multiplied by function gðxÞ ¼ expðiax2 =2bÞf ðxÞ, pffiffiffiffiffiffiffiffiffiffi 2 expðidy =2bÞ= 2pib. On the other hand, note that for the case z ¼ 7i, (7) yields a discrete Fourier transform ðF 7 i Þjk C e 7 itj tk Dtk , that can be related to the standard DFT as follows. The use of (6) yields
k¼1
is an approximation of the product of functions pffiffiffiffiffiffiffiffiffiffi 2 evaluated at ðexpðidy =2bÞ= 2pibÞ1 Lfa,b,c,dg ½f ðxÞ,y yj ¼ 4bxj =p. 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
One
gðyj Þ ¼
Two iyj x
e
f ðxÞ dx,
Three
1
1
G ¼ SDF u,
eikyj x f ðxÞ dt ¼ gðkyj Þ,
ð11Þ
1
has the quadrature
PN
k¼1
transform GðyÞ ¼ Lfa,b,c,dg ½f ðxÞ,y evaluated at yj ¼ 4bxj =p. Set up the vectors u, y and s for k=1,2,y,N ! 2 ðk1ÞðN1Þ iaxk 2kN1 f p pffiffiffiffiffiffiffi uk ’exp ip exp N 2b 2 2N pffiffiffiffiffiffiffiffiffi yk ’4bxk =p ¼ 2=N bð2kN1Þ " # pffiffiffiffi p ðN1Þ2 idy2k N1 pexp i ðk1Þ þ ip N 2 N 2b pffiffiffiffiffiffiffiffi sk ’ 2 ibN end for Set up the diagonal matrix S S= diag(s1,s2,y,sN) Let DF be the discrete Fourier transform, ðDF Þjk ¼ expði 2Np jkÞ, j,k= 0,1,2,y,N 1. Obtain the approximation Gj to Gð4b p xj Þ by computing the matrixvector product
a scaled Fourier transform Z
ðS1 FS2 Þjk f ðxk Þ,
Algorithm XFT pffiffiffiffiffiffiffi Input: A list of values f ðpð2jN1Þ=2 2N Þ of a function f(x) at the pffiffiffiffiffiffiffi points xj ¼ ½jðN þ 1Þ=2p= 2N . Output: A list of values Gj approximating the linear canonical
P where we have used (6) and (8). Since N k ¼ 1 ðF i Þjk f ðxk Þ is a quadrature and therefore, an approximation of 1
N X
where S1 and S2 are diagonal matrices whose diagonal pffiffiffiffiffiffiffiffiffiffi 2 2 elements are eidyj =2b = 2pib, and eiaxj =2b , respectively. As it can be seen, the matrix L= S1FS2, 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 XFT algorithm for the linear canonical transform can be given straightforwardly.
2 p p N1 N1 k , ðF i Þjk ¼ e 7 itj tk Dtk ¼ pffiffiffiffiffiffiffi exp i j 2 2 2N 2N ð10Þ
Z
2
Fjk eiaxk =2b f ðxk Þ
k¼1
Firstly, note that if ba0, the LCT can be written as a chirp-FT-chirp transform 2
1
If now we replace f(x) by expðiax2 =2bÞf ðxÞ and take into account (10), we have that
Gðyj Þ ¼ Gð4bxj =pÞ ¼
3. A fast linear canonical transform
idy 2b pffiffiffiffiffiffiffiffiffiffi 2pib
for j,k=0,1,2,y,N 1, and k ¼ 1 Fjk f ðxk Þ is an approximation of gð4yj =pÞ. If now we choose k ¼ 4b=p, but we keep PN the same matrix (13), then k ¼ 1 Fjk f ðxk Þ is an approximation of Z 1 eiðyj =bÞx f ðxÞ dt:
N X
k¼1
exp
½eipððN1Þ=NÞj ½eið2p=NÞjk ½eipððN1Þ=NÞk , ð13Þ
k¼1 " # rffiffiffiffiffiffiffiffiffiffiffi N ð1 þ z2 Þðx2j þ x2k Þ4xj xk z 2 X f ðxk ÞDxk , C exp 1z2 k ¼ 1 2ð1z2 Þ
N X
=NÞ
pffiffiffiffiffiffiffi 2N
F~ jk f ðxk Þ, where
p p2 N1 N1 k : j F~ jk ¼ pffiffiffiffiffiffiffi exp ik 2 2 2N 2N
with a standard FFT algorithm.
4. An application to edge detection ð12Þ
As an untested and unexplored application of the algorithm given before, we use the XFT to find the edges
R.G. Campos, J. Figueroa / Signal Processing 91 (2011) 1444–1447
1447
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 comes from a quadrature formula, it yields accurate results. Fig. 1. (a) Original image. (b) Contoured image obtained by thresholding the two-dimensional XFT of the original image with parameters a1 = 1, d1 = 1, b1 =20, a2 = 1, d2 =1, b2 = 20, and threshold value Gt = 20.
Fig. 2. (c) Contoured image obtained by thresholding the two-dimensional XFT of (a) with parameters a1 = 0.2, d1 =10, b1 = 2.5, a2 = 0.2, d2 =10, b2 = 2, and threshold value Gt = 60. (d) Same parameters as in (c) but threshold value Gt = 100.
in an image. The problem of edge detection has been studied extensively since many years ago and many algorithms have been devised [18–21]. Our purpose in this section is only to point out a possible use of the XFT as an edge detector. Since an edge in an image is essentially a boundary between regions reflecting different amounts of energy, a cut-off to zero in the power spectrum of the image for values greater than a threshold value, let us say Gt, will define the boundary in the spatial domain. This simple idea is applied to the twodimensional XFT of a grayscale image with parameters a1, b1 and d1 for one dimension and a2, b2 and d2 for the other. The result of this procedure is shown in Figs. 1 and 2. In order to compare and contrast the processed images, the intensity of Figs. 1(b) and 2(c) and (d) has been doubled. As expected, different values of the parameters yield different contoured images and high values of the threshold parameter yield imprecise edges.
References [1] M. Moshinsky, C. Quesne, Linear canonical transformations and their unitary representations, J. Math. Phys. 12 (1971) 1772–1783. [2] K.B. Wolf, Integral Transforms in Science and Engineering, Plenum Press, New York, 1979 (Chapter 9–10). [3] J.J. Healy, J.T. Sheridan, Sampling and discretization of the linear canonical transform, Signal Process. 89 (2009) 641–648. [4] A. Koc- , H.M. Ozaktas, C. Candan, M.A. Kutay, Digital computation of linear canonical transforms, IEEE Trans. Signal Process. 56 (2008) 2383–2394. [5] K.K. Sharma, S.D. Joshi, Uncertainty principle for real signals in the linear canonical transform domains, IEEE Trans. Signal Process. 56 (2008) 2677–2683. [6] A. Stern, Sampling of linear canonical transformed signals, Signal Process. 86 (2006) 1421–1425. [7] B.M. Hennelly, J.T. Sheridan, Fast numerical algorithm for the linear canonical transform, J. Opt. Soc. Am. A 22 (2005) 928–937. [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. A 22 (2005) 917–927. [9] B.Z. Li, R. Tao, Y. Wang, New sampling formulae related to linear canonical transform, Signal Process. 87 (2007) 983–990. [10] S.C. Pei, J.J. Ding, Eigenfunctions of linear canonical transform, IEEE Trans. Signal Process. 50 (2002) 11–26. [11] R.G. Campos, J. Rico-Melgoza, E. Cha´vez, XFT: extending the digital application of the Fourier transform, arXiv:0911.0952v1, /http:// arxiv.org/abs/0911.0952v1S; 2009. [12] R.G. Campos, L.Z. Jua´rez, A discretization of the continuous Fourier transform, Il Nuovo Cimento 107B (1992) 703–711. [13] R.G. Campos, A quadrature formula for the Hankel transform, Numer. Algorithms 9 (1995) 343–354. [14] R.G. Campos, F. Domı´nguez Mota, E. Coronado, Quadrature formulas for integrals transforms generated by orthogonal polynomials, IMA J. Numer. Anal., to appear (cf. arXiv.0805.2111v1, 2008, /http://arxiv.org/abs/0805.2111S). [15] G. Szego, Orthogonal Polynomials, Colloquium Publications, American Mathematical Society, Providence, Rhode Island, 1975. [16] A. Erde´lyi, Higher Transcendental Functions, vols. I and II, McGraw-Hill, New York, 1953. [17] V. Namias, The fractional order Fourier transform and its application to quantum mechanics, J. Inst. Math. Appl. 25 (1980) 241–265. [18] D. Ziou, S. Tabbone, Edge detection techniques: an overview, Int. J. Pattern Recognit. Image Anal. 8 (1998) 537–559. [19] D. Marr, E. Hildreth, Theory of edge detection, Proc. R. Soc. London 207 (1980) 187–217. [20] W.Y. Ma, B.S. Manjunath, Edge flow: a framework of boundary detection and image segmentation, in: IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’97), 1997, p. 744. [21] Y. Deng, B.S. Manjunath, H. Shin, Color image segmentation, in: IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’99), vol. 2, 1999, p. 2446.