A fast algorithm for the linear canonical transform

Rafael G. Campos

, Jared Figueroa
Facultad de Ciencias F sico-Matematicas, Universidad Michoacana, 58060 Morelia, Michoacan, Me xico
a r t i c l e i n f o
Article history:
Received 19 January 2010
Received in revised form
1 July 2010
Accepted 9 July 2010
Keywords:
Linear canonical transform
Fractional Fourier transform
Quadrature
Hermite polynomials
Fractional discrete Fourier transform
fft
a b s t r a c t
In recent years there has been a renewed interest in nding 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 ONlogN algorithmto 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.
1. Introduction
The linear canonical transform (LCT) of a given function
f(x) is a three-parameter integral transformthat was obtained
in connection with canonical transformations in Quantum
Mechanics [1,2]. It is dened by
L
fa,b,c,dg
f x,y
1

2pib
p
_
1
1
e
i=2bax
2
2xydy
2

f x dx,
for ba0, and by

d
p
e
i=2cdy
2
f dy, if b=0. The four parameters
a, b, c and d appearing in the above expression, are the
elements of a 22 matrix with unit determinant, i.e.,
adbc=1. Therefore, only three parameters are free. Since
this transformis a useful tool for signal processing and optical
analysis, its study and direct computation in digital compu-
ters have become an important issue [310], 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 ONlog N complexity.
In this paper we present an algorithm that takes
ONlogN 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 O1=N
convergent to the continuous Fourier transform for certain
class of functions [13].
2. A discrete fractional Fourier transform
In previous work [1214], we derived a quadrature
formula for the continuous Fourier transform which yields
Contents lists available at ScienceDirect
journal homepage: [Link]/locate/sigpro
Signal Processing
0165-1684/$ - see front matter & 2010 Elsevier B.V. All rights reserved.
doi:10.1016/[Link].2010.07.007

Corresponding author.
E-mail addresses: rcampos@[Link] (R.G. Campos),
jared@[Link] (J. Figueroa).
Signal Processing ] (]]]]) ]]]]]]
Please cite this article as: R.G. Campos, J. Figueroa, A fast algorithm for the linear canonical transform, Signal Process.
(2010), doi:10.1016/[Link].2010.07.007
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
H
n
(x), n=0,1,y, which satises the recurrence equation:
H
n1
x 2nH
n1
x 2xH
n
x, 1
with H
1
x 0. Note that the recurrence equation (1) can
be written as the eigenvalue problem
0 1=2 0
1 0 1=2
0 2 0
^ ^ ^ &
_
_
_
_
_
_
_
_
_
_
_
_
H
0
x
H
1
x
H
2
x
^
_
_
_
_
_
_
_
_
_
_
_
_
x
H
0
x
H
1
x
H
2
x
^
_
_
_
_
_
_
_
_
_
_
_
_
: 2
Let us now consider the eigenproblem associated to the
principal submatrix of dimension N of (2)
H
0 1=2 0 0 0
1 0 1=2 0 0
0 2 0 0 0
^ ^ ^ & ^ ^
0 0 0 0 1=2
0 0 0 N1 0
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
:
It is convenient to symmetrize H by using the similarity
transformation SHS
1
where S is the diagonal matrix
S diag 1,
1

2
p , . . . ,
1

N1!2
N1
_
_

_
_

_
:
Thus, the symmetric matrix H SHS
1
takes the form
0

1
2
_
0 0 0

1
2
_
0

2
2
_
0 0
0

2
2
_
0 0 0
^ ^ ^ & ^ ^
0 0 0 0

N1
2
_
0 0 0

N1
2
_
0
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
:
The recurrence equation (1) and the ChristoffelDarboux
formula [15] can be used to solve the eigenproblem
Hu
k
x
k
u
k
, k 1,2, . . . ,N,
which is a nite-dimensional version of (2). The eigenva-
lues x
k
are the zeros of H
N
(x) and the kth eigenvector u
k
is
given by
c
k
s
1
H
0
x
k
,s
2
H
1
x
k
, . . . ,s
N
H
N1
x
k

T
,
where s
1
,y,s
N
are the diagonal elements of S and c
k
is a
normalization constant that can be determined from the
condition u
k
T
u
k
=1, i.e., from
c
2
k

N1
n 0
H
n
x
k
H
n
x
k

2
n
n!
1:
Therefore,
c
k

2
N1
N1!
N

1
Nk
H
N1
x
k

:
Thus, the components of the orthonormal vectors u
k
,
k=1,2,y, N, are
u
k

n
1
Nk

2
Nn
N1!
Nn1!

H
n1
x
k

H
N1
x
k

, 3
n=1,y,N. Let U be the orthogonal matrix whose kth
column is u
k
and let us dene the matrix
F
z

2p
p
U
1
DzU,
where D(z) is the diagonal matrix
D(z)=diag{1,z,z
2
,y,z
N1
} and z is a complex number.
Therefore, the components of F
z
are given by
F
z

jk

2p
p
1
j k
2
N1
N1!
NH
N1
x
j
H
N1
x
k

N1
n 0
z
n
2
n
n!
H
n
x
j
H
n
x
k
: 4
Next, we want to prove that if N is large enough, (4)
approaches the kernel of the fractional Fourier transform
evaluated at x=x
j
, y=x
k
. To this, we use the asymptotic
expression for H
N
(x) [15]
H
N
xC
GN1e
x
2
=2
GN=21
cos

2N1
_
x
Np
2
_ _
: 5
Thus, the asymptotic form of the zeros of H
N
(x) are
x
k

2kN1

2N
p
_ _
p
2
, 6
k=1,2,y,N. The use of (5) and (6) yields
H
N1
x
k
C1
Nk
GN
G
N1
2
_ _e
x
2
k
=2
, N-1,
and the substitution of this asymptotic expression in (4)
yields
F
z

jk
C

2p
p
2
N1
G
N1
2
_ _ _ _
2
GN1
e
x
2
j
x
2
k
=2

1
n 0
z
n
2
n
n!
H
n
x
j
H
n
x
k
:
Finally, Stirlings formula and Mehlers formula [16]
produce
F
z

jk
C

2
1z
2
_
exp
1z
2
x
2
j
x
2
k
4x
j
x
k
z
21z
2

_ _
Dx
k
, 7
where Dx
k
is the difference between two consecutive
asymptotic Hermite zeros, i.e.,
Dx
k
x
k1
x
k

p

2N
p : 8
Let us consider now the vector of samples of a given
function f(x)
f f x
1
,f x
2
, . . . ,f x
N

T
:
Please cite this article as: R.G. Campos, J. Figueroa, A fast algorithm for the linear canonical transform, Signal Process.
(2010), doi:10.1016/[Link].2010.07.007
R.G. Campos, J. Figueroa / Signal Processing ] (]]]]) ]]]]]] 2
The multiplication of the matrix F
z
by the vector f gives
the vector g with entries
g
j

N
k 1
F
z

jk
f x
k

2
1z
2
_

N
k 1
exp
1z
2
x
2
j
x
2
k
4x
j
x
k
z
21z
2

_ _
f x
k
Dx
k
,
where j =1,2y,N. This equation is a Riemann sum for the
integral
F
z
f x,y

2
1z
2
_
_
1
1
exp
1z
2
y
2
x
2
4xyz
21z
2

_ _
f x dx,
where jzj o1. Therefore, if we make y
j
=x
j
,
F
z
f x,y
j
C

N
k 1
F
z

jk
f x
k
, N-1: 9
Note that F
z
gx,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.
3. A fast linear canonical transform
Firstly, note that if ba0, the LCT can be written as a
chirp-FT-chirp transform
L
fa,b,c,dg
f x,y
exp
idy
2
2b
_ _

2pib
p
_
1
1
exp
ixy
b
_ _
exp
iax
2
2b
_ _
f x dx:
Thus, for ba0, the LCT of the function f(x) can be
represented by the 1/b-scaled Fourier transform of the
function gx expiax
2
=2bf x, multiplied by
expidy
2
=2b=

2pib
p
.
On the other hand, note that for the case z 7i, (7)
yields a discrete Fourier transform F
7i

jk
Ce
7it
j
t
k
Dt
k
,
that can be related to the standard DFT as follows. The use
of (6) yields
F
i

jk
e
7it
j
t
k
Dt
k

p

2N
p exp i
p
2
2N
j
N1
2
_ _
k
N1
2
_ _ _ _
,
10
where we have used (6) and (8). Since

N
k 1
F
i

jk
f x
k
is a
quadrature and therefore, an approximation of
gy
j

_
1
1
e
iy
j
x
f x dx,
a scaled Fourier transform
_
1
1
e
iky
j
x
f x dt gky
j
, 11
has the quadrature

N
k 1
~
F
jk
f x
k
, where
~
F
jk

p

2N
p exp ik
p
2
2N
j
N1
2
_ _
k
N1
2
_ _ _ _
: 12
If we choose k 4=p, (12) takes the form
F
jk

pe
ip=2N1
2
=N

2N
p e
ipN1=Nj
e
i2p=Njk
e
ipN1=Nk
,
13
for j,k=0,1,2,y,N1, and

N
k 1
F
jk
f x
k
is an approxima-
tion of g4y
j
=p. If now we choose k 4b=p, but we keep
the same matrix (13), then

N
k 1
F
jk
f x
k
is an approx-
imation of
_
1
1
e
iy
j
=bx
f x dt:
If now we replace f(x) by expiax
2
=2bf x and take into
account (10), we have that

N
k 1
F
jk
e
iax
2
k
=2b
f x
k

is an approximation of the product of functions

2pib
p

1
L
fa,b,c,dg
f x,y evaluated at
y
j
4bx
j
=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
Gy
j
G4bx
j
=p

N
k 1
S
1
FS
2

jk
f x
k
,
where S
1
and S
2
are diagonal matrices whose diagonal
elements are e
idy
2
j
=2b
=

2pib
p
, and e
iax
2
j
=2b
, respectively. As it
can be seen, the matrix L=S
1
FS
2
, which gives the discrete
LTC, i.e., the XFT, consists in a chirp-DFT-chirp transfor-
mation, 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.
Algorithm XFT
Input:
A list of values f p2jN1=2

2N
p
of a function f(x) at the
points x
j
jN1=2p=

2N
p
.
Output: A list of values G
j
approximating the linear canonical
transform Gy L
fa,b,c,dg
f x,y evaluated at y
j
4bx
j
=p.
One Set up the vectors u, y and s
for k=1,2,y,N
u
k
exp ip
k1N1
N
_ _
exp
iax
2
k
2b
_ _
f p
2kN1
2

2N
p
_ _
y
k
4bx
k
=p

2=N
_
b2kN1
s
k

p
p
exp i
p
2
N1
2
N

idy
2
k
2b
ip
N1
N
k1
_ _
2

ibN
p
end for
Two Set up the diagonal matrix S
S=diag(s
1
,s
2
,y,s
N
)
Three Let D
F
be the discrete Fourier transform,
D
F

jk
expi
2p
N
jk, j,k=0,1,2,y,N1. Obtain the
approximation G
j
to G
4b
p
x
j
by computing the matrix-
vector product
G SD
F
u,
with a standard FFT algorithm.
4. An application to edge detection
As an untested and unexplored application of the
algorithm given before, we use the XFT to nd the edges
Please cite this article as: R.G. Campos, J. Figueroa, A fast algorithm for the linear canonical transform, Signal Process.
(2010), doi:10.1016/[Link].2010.07.007
R.G. Campos, J. Figueroa / Signal Processing ] (]]]]) ]]]]]] 3
in an image. The problem of edge detection has been
studied extensively since many years ago and many
algorithms have been devised [1821].
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
reecting 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 G
t
, will dene the boundary in
the spatial domain. This simple idea is applied to the two-
dimensional XFT of a grayscale image with parameters a
1
,
b
1
and d
1
for one dimension and a
2
, b
2
and d
2
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.
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 nite 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 imple-
ment it in a fast algorithm. Since this discrete LCT comes
from a quadrature formula, it yields accurate results.
References
[1] M. Moshinsky, C. Quesne, Linear canonical transformations and
their unitary representations, J. Math. Phys. 12 (1971) 17721783.
[2] K.B. Wolf, Integral Transforms in Science and Engineering, Plenum
Press, New York, 1979 (Chapter 910).
[3] J.J. Healy, J.T. Sheridan, Sampling and discretization of the linear
canonical transform, Signal Process. 89 (2009) 641648.
[4] A. Koc-, H.M. Ozaktas, C. Candan, M.A. Kutay, Digital computation of
linear canonical transforms, IEEE Trans. Signal Process. 56 (2008)
23832394.
[5] K.K. Sharma, S.D. Joshi, Uncertainty principle for real signals in the
linear canonical transform domains, IEEE Trans. Signal Process. 56
(2008) 26772683.
[6] A. Stern, Sampling of linear canonical transformed signals, Signal
Process. 86 (2006) 14211425.
[7] B.M. Hennelly, J.T. Sheridan, Fast numerical algorithm for the linear
canonical transform, J. Opt. Soc. Am. A 22 (2005) 928937.
[8] B.M. Hennelly, J.T. Sheridan, Generalizing, optimizing, and invent-
ing numerical algorithms for the fractional Fourier, Fresnel, and
linear canonical transforms, J. Opt. Soc. Am. A 22 (2005) 917927.
[9] B.Z. Li, R. Tao, Y. Wang, New sampling formulae related to linear
canonical transform, Signal Process. 87 (2007) 983990.
[10] S.C. Pei, J.J. Ding, Eigenfunctions of linear canonical transform, IEEE
Trans. Signal Process. 50 (2002) 1126.
[11] R.G. Campos, J. Rico-Melgoza, E. Cha vez, XFT: extending the digital
application of the Fourier transform, arXiv:0911.0952v1, /http://
[Link]/abs/0911.0952v1S; 2009.
[12] R.G. Campos, L.Z. Jua rez, A discretization of the continuous Fourier
transform, Il Nuovo Cimento 107B (1992) 703711.
[13] R.G. Campos, A quadrature formula for the Hankel transform,
Numer. Algorithms 9 (1995) 343354.
[14] R.G. Campos, F. Domnguez Mota, E. Coronado, Quadrature
formulas for integrals transforms generated by orthogonal poly-
nomials, IMA J. Numer. Anal., to appear (cf. arXiv.0805.2111v1,
2008, /[Link]
[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 applica-
tion to quantum mechanics, J. Inst. Math. Appl. 25 (1980) 241265.
[18] D. Ziou, S. Tabbone, Edge detection techniques: an overview, Int. J.
Pattern Recognit. Image Anal. 8 (1998) 537559.
[19] D. Marr, E. Hildreth, Theory of edge detection, Proc. R. Soc. London
207 (1980) 187217.
[20] W.Y. Ma, B.S. Manjunath, Edge ow: a framework of boundary
detection and image segmentation, in: IEEE Computer Society
Conference on Computer Vision and Pattern Recognition (CVPR97),
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 (CVPR99), vol. 2, 1999, p. 2446.
Fig. 1. (a) Original image. (b) Contoured image obtained by thresholding
the two-dimensional XFT of the original image with parameters a
1
=1,
d
1
=1, b
1
=20, a
2
=1, d
2
=1, b
2
=20, and threshold value G
t
=20.
Fig. 2. (c) Contoured image obtained by thresholding the two-dimen-
sional XFT of (a) with parameters a
1
=0.2, d
1
=10, b
1
=2.5, a
2
=0.2, d
2
=10,
b
2
=2, and threshold value G
t
=60. (d) Same parameters as in (c) but
threshold value G
t
=100.
Please cite this article as: R.G. Campos, J. Figueroa, A fast algorithm for the linear canonical transform, Signal Process.
(2010), doi:10.1016/[Link].2010.07.007
R.G. Campos, J. Figueroa / Signal Processing ] (]]]]) ]]]]]] 4