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