IEEE TRANSACTIONS ON SIGNAL PROCESSING,
VOL.
44,
NO.
9,
SEPTEMBER
1996
2141
Digital Computation
of
the
Fractional Fourier Transform
Haldun
M.
Ozaktas, Orhan Ankan,
Member,
IEEE,
M.
Alper
Kutay,
Student Member,
IEEE,
and Gozde Bozdaki,
Member,
IEEE
Abstract-An algorithm for efficient and accurate computation
of the fractional Fourier transform is given. For signals with
time-bandwidth product
N,
the presented algorithm computes
the fractional transform in
O(
N
log
N)
time. A definition for
the discrete fractional Fourier transform that emerges from our
analysis is also discussed.
I. INTRODUCTION
HE
fractional Fourier transform
[
11-[4] has found many
T
applications in the solution of differential equations [2],
[31,
quantum mechanics and quantum optics
[5]-[
111, optical
diffraction theory and optical beam propagation (including
lasers), and optical systems and optical signal processing
[
11,
[
131-[25], swept-frequency filters [4], time-variant filtering
and multiplexing
[I],
pattern recognition
[26],
and study of
time-frequency distributions [27]. The recently studied Radon
transformation of the Wigner spectrum [281-[30] is also known
to be the magnitude square of the fractional Fourier transform
[ll, [31]. The fractional Fourier transform has been related to
wavelet transforms [l], [32], neural networks [32], and is also
related to various chirp-related operations [l], [33]-[35]. It can
be optically realized much like the usual Fourier transform
[l], [131-[17], [20] and, as we will show in this paper, can
be simulated with a fast digital algorithm. Other applications
that are currently under study or have been suggested include
phase retrieval, signal detection, radar, tomography, and data
compression.
In this paper, we will be concerned with the digital com-
putation of the fractional Fourier transform. We are not only
interested in a numerical method to compute the continuous
transform but also in defining the discrete fractional Fourier
transform and show how it can be used to approximate the
continuous transform. More precisely, we will show that the
samples of the continuous time fractional Fourier transform
of a function can be approximately evaluated in terms of the
samples of the original function in
O(N
log
N)
time, where
N
is the time-bandwidth product of the signal.
In many of the above-mentioned applications, it is possible
to improve performance by use of the fractional Fourier
transform instead of the ordinary Fourier transform. Since
in this paper we show that the fractional transform can be
Manuscript received February 3, 1995, revised January 9, 1996. The
associate editor coordinating the review
of
this paper and approving it for
publication
was
Dr. Bruce Suter.
The authors are with Electrical Engineering, Bilkent University, 06533
Bilkent, Ankara, Turkey.
Publisher Item Identifier
S
1053-587X(96)06680-9.
computed in about the same time as the ordinary transform,
these performance improvements come at no additional cost.
To give one concrete example, in some cases, filtering in a
fractional Fourier domain, rather than the ordinary Fourier
domain, allows one to decrease the mean square error in
estimating a distorted and noisy signal [36].
In Section 11, preliminaries about the fractional Fourier
transform are given, including its relation to the Wigner
distribution. Section
I11
reviews some straightforward yet
inefficient methods
of
computing the fractional Fourier trans-
form. Fast computational algorithms are presented in Section
IV, and simulation examples are given in Section V. Some
alternate methods are discussed in Section VI to better situate
the suggested algorithm among other possibilities. Section
VI1 deals with the issue of defining the discrete fractional
Fourier transform in some detail. The remainder of the paper
constitutes concluding sections.
11. PRELIMINARIES
A.
The
Fractional Fourier Transform
Let
{
3Tf)
(x)
denote the Fourier transform of
f
(s)
.
Integral
powers
3-7
of the operator
3
e
3'
may be defined as its
successive applications. Then, we have
{3'f}(s)
=
f(-x)
and
{34f}(x)
=
f(s).
The ath-order fractional Fourier
transform
{F'"f}(x)
of the function
f(x)
may be defined for
O<(a(<2
as
00
F'"[f(s)]
?E
{F"f}(s)
=
Ba(x,z')f(x')
ds',
1,
Ba(z,
2')
=
A+
exp
[ZT(S'
cot
4
-
2x2'
csc
q5
+
5''
cot
$)I,
exp
(-i~
sgn (sin
4)/4
+
@/a)
1
sin
$(1/2
(1)
A+
where
and
i
is the imaginary unit. The kernel approaches
Bg(z,
s')
E
S(z
-
d)
and
B%2(x,x')
S(z
+
2')
for
a
=
0
and
a
=
f.2,
respectively. The definition is easily extended outside
the interval
[-2,2]
by remembering that
343
is the identity
operator for any integer
j
and that the fractional Fourier
transform operator is additive in index, that is,
3'"13a2
=
.FTa1+'"2.
A
complete set of eigenfunctions of the fractional
1053-587X/96$05.00
0
1996
IEEE
2142
IEEE
TRANSACTIONS
ON
SIGNAL PROCESSING,
VOL
44, NO
9,
SEPTEMBER
1996
Fourier transform are the Hermite-Gaussian functions
where
H,(z)
is the nth-order Hermite polynomial. The spec-
tral expansion of the linear transform kernel is
CO
B,(x,x’)
=
e-2ann’2$,,(Z)t/jn(5’).
(5)
n=O
Two and higher dimensional transforms
[
141-[
181 have sepa-
rable kernels
so
that most results easily generalize to higher
dimensions. Proofs of these and other properties may be found
in [1]-[4],
[14]-[18],
and [31].
The ath fractional Fourier transform
{Faf}(x)
of
the
function
f(x)
will be abbreviatedly denoted by
fa
(x).
As a word on terminology, we believe that ultimately, the
term “Fourier transform” should mean, in general, “fractional
Fourier transform” and that the presently standard Fourier
transform be referred to as the “first-order Fourier trans-
form.” Likewise, DFT should stand for the discrete (fractional)
Fourier transform, etc., and the invention of new acronyms and
abbreviations should be discouraged.
To
avoid confusion, we note that for
a
=
1,
the frac-
tional transform reduces to the ordinary transform defined as
1
f(z’)
exp
(-227rzz’)
dz’.
B.
Relation to the Wigner Distribution and the
Concept
of
Fractional Fourier Domains
The Wigner distribution
Wf(z:,v)
of a signal
f
can be
defined in terms of the time-domain representation
f(x)
of
that signal as
Wf(x,
V)
=
1,
f(z
+
z’/a)f*(z
-
x’/2)e-22Tv2’
dd.
(6)
Roughly speaking,
W(z,v)
is a function that gives the dis-
tribution of signal energy over time and frequency. Properties
of the Wigner distribution may be found in [37] and [38]. We
note the following:
00
00
i,
W(x:,
U)
dv
=
lfb)I2,
(7)
(8)
(9)
The Wigner distribution can also be defined in terms of any
of the fractional transforms of
f(x)
and can be written as a
function of other coordinate variables in the
x-U
plane. Thus,
it should be considered
to
be a geometric entity associated
with the signal
f
in the abstract and not tied to a particular
representation of
f
in a given domain.
It is possible to show that the Wigner distribution of
fa(z)
is merely a rotated version of that of
f(x)
[I], [4],
[lS],
[31]
Wfa
(IC,
v)
=
Wf(z
cos
4
-
v
sin
4,
z
sin
4
+
v
cos
4).
(IO)
CO
W(z,
U)
dz
=
lfl(U)I2,
COCO
W(z,
U)
dx
dv
=
Signal energy.
.i,
.1,
x,
=
v
Fig.
1. Fractional Fourier
domains.
The same property can be stated in the alternative form [l],
[41,
~311
PdW?
V)I>(xa)
=
lfa(xa)l2
(1
1)
where
R,
is the Radon transform operator.
724
takes the
integral projection of the
2-D
function
Wf(z,
U)
onto an axis
making angle
=
with the
z
axis. We will refer to
this axis as the
z,
axis or the ath fractional Fourier domain
(Fig.
1).
The
xo
axis is the usual time domain
z,
and the
z1
axis is the usual frequency domain
U.
Notice that
(7)
and
(8)
are special cases of this equation. In general, the projection of
the Wigner distribution on the ath fractional Fourier domain
gives the magnitude squared of the ath fractional Fourier
transform of the original function.
There is actually nothing special about any of the continuum
of domains; the privileged status we assign to the time and
frequency domains can be interpreted as an arbitrary choice of
the origin of the parameter
a.
All of the fractional transforms,
including the 0th transform (the function itself), are different
functional representations of an abstract signal in different
domains. The unitary transformation between these different
representations is the fractional Fourier transform
[
11.
C.
Compactness in the Time Domain, the
Frequency Domain, and Wigner Space
A function will be referred to as compact if its support is
so.
The support of a function is the subset of the real axis in which
the function is not equal to zero. In other words,
a
function
is compact if and only if its nonzero values are confined to a
finite interval. It is well known that
a
function and its Fourier
transform cannot be both compact (unless they are identically
zero). In practice, however, it seems that we are always
working with a finite time interval and a finite bandwidth.
This discrepancy between our mathematical idealizations and
OZAKTAS
et
al.:
DIGITAL COMPUTATION OF THE FRACTIONAL FOURIER TRANSFORM
2143
the real world is usually not a problem when we work with
signals of large time-bandwidth product. The time-bandwidth
product can be crudely defined as the product of the temporal
extent of the signal and its (double-sided) bandwidth. It is
equal to the number of degrees of freedom and the number of
complex numbers required to uniquely characterize the signal
among others of the same time-bandwidth product.
We will assume that the time-domain representation of our
signal is approximately confined to the interval [-Atla, At121
and that its frequency-domain representation is confined to the
interval [-Af/2,
Af/Z].
With this statement, we mean that a
sufficiently large percentage of the signal energy is confined
to these intervals. For a given class of functions, this can
be ensured by choosing
At
and
Af
sufficiently large. We
then define the time-bandwidth product
N
z
AtA
f,
which is
always greater than unity because of the uncertainty relation.
Let us now introduce the scaling parameter
s
with the
dimension of time and introduce scaled coordinates
x
=
t/s
and
U
=
fs.
With these new coordinates, the time and
frequency domain representations will be confined to intervals
of length Atls and
A
f
s.
Let us choose
s
=
dm
so
that
the lengths of both intervals are now equal to the dimensionless
quantity
dm,
which we will denote by
Ax.
In the newly
defined coordinates, our signal can be represented in both
domains with
N
=
Ax2
samples spaced
Axp1
=
1/fi
apart.
From now
on,
we will assume that this dimensional normal-
ization has been performed and that the coordinates appearing
in the definition of the fractional Fourier transform, Wigner
distribution, etc., are all dimensionless quantities.
If the representation of the signal in the ath domain is
confined to a certain interval around the origin, the Wigner
distribution will be confined to an infinite strip perpendicular
to the
xa
axis defined by that interval. Thus, assuming that
the representation
of
the signal in
all
domains is confined
to an interval of length
Ax
around the origin, is equivalent
to assuming that the Wigner distribution is confined within a
circle of diameter
Ax,
With this, we mean that a sufficiently
large percentage of the energy
of
the signal is contained in
that circle, For any signal, this assumption can be justified by
choosing
Ax
sufficiently large. (Of course, it is in our interest
to
choose it as small as possible to reduce computational
complexity.) For convenience, we will require
Ax
to be an
integer.
Signals whose energy are not concentrated around the origin
of Wigner space might be treated more efficiently than by
simply choosing
Ax
large enough to include them, but this
extension to the “bandpass” case is not treated in this paper.
D.
The Discrete Fourier Transform
The discrete Fourier transform
(DFT)
is a mapping
RN
t
RN.
The matrix elements of this transformation are defined as
(12)
The
DFT
is related to the continuous Fourier transform as
follows [39]:
Assume that a function
f
(x)
and its Fourier transform
fl
(U)
are both confined to the interval
[-Ax/2,
Ax/2]. Then, the
N
=
Ax2
samples of the Fourier transform may be found
by taking the
DFT
of the
N
samples of the original function,
where the sample spacing in both domains is
l/fi
=
l/Ax.
A more precise statement of the above belongs to a class
of
relations known as Poisson formulas [39]:
111.
METHODS
OF
COMPUTING
THE
CONTINUOUS
FRACTIONAL FOURIER TRANSFORM
The defining integral (1) can be rarely evaluated analyti-
cally; therefore, numerical integration is called for. Numerical
integration of quadratic exponentials, which often appear in
diffraction theory, require a very large number of samples
if conventional methods are to be employed, due to the
rapid oscillations of the kernel. The problem is particularly
pronounced when
a
is close to
0
or
52.
If
we assume both
the function and its Fourier transform to be confined to a
finite interval, we can walk around this difficulty as follows:
If a
E
[0.5,1.5]
or a
E
[2.5,3.5], we evaluate the integral
directly. If a
E
(-0.5,0.5)
or
a
E
(1.5,2.5),
we use the
property
Fa
=
F1FT”-l,
noting that in this case, the
(a
-
1)th
transform can be evaluated directly. (Essentially similar issues
are discussed in
[28]-[30].)
Another method
of
evaluating (1) would be to use the
spectral decomposition
of
the kernel (see
(5))
[l], [14]-[161.
This is equivalent to first expanding the function
f(x)
as
C~?o
c,$,
(x),
multiplying the expansion coefficients
c,,
re-
spectively, with
e-1nnr/2,
and summing the components.
Although both ways of evaluating the fractional Fourier
transform may be expected to give accurate results, we do
not consider them further since they take
O(N2)
time.
Iv.
FAST
COMPUTATION
OF
THE
FRACTIONAL
FOURIER
TRANSFORM
The fractional Fourier transform is a member of a more
general class of transformations that are sometimes called
linear canonical transformations or quadratic-phase transforms
[20].
Members of this class of transformations can be broken
down into a succession of simpler operations, such as chirp
multiplication, chirp convolution, scaling, and ordinary Fourier
transformation. Here, we will concentrate on two particular
decompositions that lead to two distinct algorithms.
A.
Method
Z
First, we choose to break down the fractional transform
into a chirp multiplication followed by a chirp convolution
followed by another chirp multiplication [171,
[391.
In this approach, we assume
a
E
[-
1,
11. Manipulating
(l),
we can write
(14)
.fa(x)
=
exp
[-inz2
tan
(dP)]g’(x),
2144
IEEE
TRANSACTIONS
ON
SIGNAL
PROCESSING,
VOL.
44,
NO.
9,
SEPTEMBER
1996
00
g’(x)
=
A4
1
exp
[z~/?(x
-
d)’]g(d)
dd,
(15)
(16)
--oo
g(x)
=
exp
[-ii7x2
tan
($/2)]f(x)
where
g(x)
and
g/(z)
represent intermediate results, and
/3
=
In the first step (see (16)), we multiply the function
f(x)
by
a chirp function.
As
we discuss in the Appendix, the bandwidth
and time-bandwidth product of
g(x)
can be as large as twice
that of
f(x).
Thus, we require the samples
of
g(x)
at intervals
of
1/2Ax.
If the samples of
f(x)
spaced at
l/Ax
are given
to begin with, we can interpolate these and then multiply by
the samples of the chirp function to obtain the desired samples
of
g(z).
There are efficient ways of performing the required
interpolation
[40].
The next step is to convolve
g(x)
with a chirp function, as
given in (15). To perform this convolution, we note that since
g(x)
is
bandlimited, the chirp function can also be replaced
with its bandlimited version without any effect, that is
csc
4.
00
g’(x)
=
A4
[
exp
[i;.P(x
-
z’)’]g(x’)
dx’
J
-00
J
-00
where
AX
H(v)
exp
(i27rvz)
du
(18)
and where
1
H(v)
-e2w/4
exp
(
-i;.v’//P)
(19)
Jp
is the Fourier transform of exp
(irr/3x2).
It is possible
to
express
h(x)
explicitly in terms of the Fresnel integral defined
as
(20)
Now,
(15)
can be sampled, giving
This
convolution can be evaluated using a fast Fourier trans-
form.
Then, after performing the last step (see
(14)),
we obtain the
samples of
fa(z)
spaced at
1/2Ax.
Since we have assumed
that all transforms of
f(x)
are bandlimited to the interval
[-Ax/2,
Ax/2],
we finally decimate these samples by a factor
of 2 to obtain samples of
fa(z)
spaced at
l/Ax.
Overall, the procedure starts with
N
samples spaced at
l/Ax,
which uniquely characterize the function
f(z),
and
returns the same for
fa
(x)
.
If we let
f
and
fa
denote column
vectors with
N
elements containing the samples
of
f(x)
and
fa(z),
the overall procedure can be represented as
fa
=Gf,
(22)
F4
=
DAHl,AJ.
(23)
Here,
D
and
J
are matrices representing the decimation and
interpolation operations
[40].
A
is a diagonal matrix that
corresponds to chirp multiplication, and
Hz,
corresponds to
the convolution operation. We notice that
3’:
allows us to
obtain the samples of the ath transform in terms of the samples
of the original function, which is the basic requirement for a
definition of the discrete fractional Fourier transform matrix.
If we are merely interested in computing and plotting the
fractional Fourier transform of a given continuous
f
(z),
then
the decimation and interpolation steps can be eliminated.
Note that the described algorithm works for
-1
5
a
5
1.
If
a
lies outside this interval, we can use the properties
{F4f}(x)
=
f(x)
and
{F2f}(z)
=
f(-z)
to easily obtain
the desired transform.
B.
Method
I1
We now turn our attention to an alternative method that
does not require Fresnel integrals. The defining equation for
the fractional Fourier transform can be put in the form
00
e-z2Tpzx’
f(.’)]
dx’
JL
{Fat}(.)
=
A+ezTcuz2
(24)
where
a
=
cot
q!J
and
/3
=
csc
4.
We are again assuming
that the Wigner distribution of
f(.)
is zero outside a circle of
diameter
Ax
centered around the origin. (This was_ discussed
in detail in Section 11-C.) Under this assumption, and by
limiting the order
a
to the interval
0.5
5
la1
5
1.5,
the
amount of vertical shear in Wigner space resulting from the
chirp modulation is bounded by
Ax/2.
Then, the modulated
function
eznaxt2
f
(x‘)
is bandlimited to
Ax
in the frequency
domain. Thus,
ezTcux‘2
f(d)
can be represented by Shannon’s
interpolation formula
sinc
2Ax
x’
-
-
))
(25)
( (
2Ax
where
N
=
(Ax)2.
The summation goes from
-N
to
N
since
f(d)
is assumed to be zero outside
[-Ax/2,Ax/2].
By
using (25) and (24) and changing the order of integration
and summation, we obtain
A1
The integral is equal to
e-z2wpz(n/2AX)(
@Ax)
rect
(/3x/2Az).
For the range of
0.5
5
la1
5
1.5,
rect
(/3x/2Ax)
will always
be equal to unity on the support
1x1
5
Ax/2
of
the transformed
function. Hence, we can write
.N
OZAKTAS
et
al.:
DIGITAL COMPUTATION OF
THE
FRACTIONAL FOURIER TRANSFORM
2145
Then, the samples of the transformed function are obtained as
N
-
-
-
’4
,i~(~~(m/2A~)~-2p[mn/(2Az)~]+c~(n/2Az)~)
n=-N
2A2
which is a finite summation, allowing us to obtain the samples
of the fractional transform in terms
of
the samples of the
original function. Direct computation of this form would
require
O(N2)
multiplications. An
O(N
log
N)
algorithm can
be obtained as follows. We put
(28)
into the following form
after some algebraic manipulations:
{3Uf
1
(&)
N
-
-
&ei~(~-~)(m/2Az)’ ,i~@((m-n)/2Az)’
n=-N
2Ax
It can be reco nized that the summation is the convolution
of
einp(n/2Ax)
and the chirp-modulated function
f(
.).
The
convolution can be computed in
O(N
log
N)
time by using
the
FFT.
The output samples are then obtained by a final chirp
modulation. Hence, the overall complexity is
O(
N
log
N).
As
in method
I,
by assuming appropriate
x2
interpolation
and decimation, the procedure starts with
N
samples spaced
at
l/Az,
which uniquely characterizes the function
f
(x)
and
returns the same for
fa(%).
Again, letting
f
and
f,,
denoting
column vectors with
N
elements containing the samples of
f
(2)
and
fa(%),
the overall procedure can be represented as
F
where
K,(m,n)
=
-e
‘4
in(~~(m/2Az)~-2p[mn/(2Az)~]+a(n/2Az)~)
2Ax
for
In1
and
Iml
5
N.
Just as
FI,
we notice that
FYI
also
allows us to obtain the samples of the ath transform in terms
of the samples of the original function.
We has assumed
0.5
5
la1
5
1.5
in deriving this algorithm.
Using the index additivity property of the fractional Fourier
transform, we can extend this range to all values of
a
easily.
For instance, for the range
0
5
a
5
0.5,
we observe that
F-lP.
(33)
3”
=
~u-1+1
-
-
Since
0.5
5
la
-
11
5
1,
the algorithm derived earlier, in
conjunction with an ordinary Fourier transform, can be used
in this case as well. More concretely, since both the signal and
its Fourier transform
is
assumed to be limited
to
the interval
[-Ax/2, Az/2], samples
of
the fractional Fourier transform
-2
-1
0
1
2
(f)
Fig.
2.
(a)
Rectangle function rect(2). The magnitude of its fractional
Fourier transform of orders
(b)a
=
0.25,
(c)
a
=
0.50,
(d)
a
=
0.75,
and
(e)
a
=
1.
(f)
The phase
of
the transform of order
a
=
0.5
is also
shown.
can be related to the samples of the Fourier transform as
N
-
A@
,i.rr(a’(m/ZAz)’-Z~’[mn/(2A~)’]+a’(n/2Am)~)
n=-N
2Ax
.
w
1
(2)
(34)
where
@
=
~(a
-
1)/2,
a’
=
cot
$’,
and
,b”
=
CSC~’.
In
this case
Fa
11
-Fa-lF
-
11
(35)
where
F
is the ordinary
DFT
matrix.
V.
EXAMPLES
Both of the above presented fast methods give results that
are in perfect agreement for all of the examples we tried. We
prefer Method 11, which does not involve Fresnel integrals,
since it
is
faster than the first.
We first tested our algorithm by calculating the fractional
Fourier transform of the Hermite-Gauss functions. We verified
(4)
for the first eight orders with excellent precision.
We then evaluated the Fourier transform of the common
rectangle function (Fig.
2).
It is interesting to see the evolution
of the rectangle function continuously to the sinc function as