Fast algorithm for chirp transforms with
zooming-in ability and its applications
Xuegong Deng, Bipin Bihari, Jianhua Gan, Feng Zhao, and Ray T. Chen
Microelectronics Research Center, Department of Electrical and Computer Engineering, The University of Texas,
Austin, Texas 78758
Received July 16, 1999; revised manuscript received December 7, 1999; accepted December 22, 1999
A general fast numerical algorithm for chirp transforms is developed by using two fast Fourier transforms and
employing an analytical kernel. This new algorithm unifies the calculations of arbitrary real-order fractional
Fourier transforms and Fresnel diffraction. Its computational complexity is better than a fast convolution
method using Fourier transforms. Furthermore, one can freely choose the sampling resolutions in both x and
u space and zoom in on any portion of the data of interest. Computational results are compared with ana-
lytical ones. The errors are essentially limited by the accuracy of the fast Fourier transforms and are higher
than the order 10
⫺12
for most cases. As an example of its application to scalar diffraction, this algorithm can
be used to calculate near-field patterns directly behind the aperture, 0 ⭐ z
⬍
d
2
/. It compensates another
algorithm for Fresnel diffraction that is limited to z
⬎
d
2
/N [J. Opt. Soc. Am. A 15, 2111 (1998)]. Experi-
mental results from waveguide-output microcoupler diffraction are in good agreement with the calculations.
© 2000 Optical Society of America [S0740-3232(00)01704-X]
OCIS codes: 350.6980, 070.2590, 050.1590, 050.1940
1. INTRODUCTION
Responses of many physical systems can be described
with chirp transforms (ChT’s). For example, chirps on
laser pulses,
1
scalar diffraction through a first-order opti-
cal system,
2–6
holographic lenses,
7
and Fresnel
transforms
8
(FnT’s) are the most frequently reported
techniques used in scalar diffraction calculations. Re-
cent developments in optical interconnects spurred by
high-speed and huge-capacity optical communications in-
clude problems similar to those just mentioned.
9–12
One
example of such a system is illustrated in Fig. 1, where
the optical signal is transmitted from board to board and
is detected by detectors located several tens of wave-
lengths away. The near-field diffraction pattern at the
detector ends directly affects the performance of the inter-
connection. It is preferable to obtain complete informa-
tion on the evolution of the optical signals along the in-
terconnecting path. In addition, it is easy to show that
fractional-order Fourier transforms
13
(FrFT’s) can also be
treated as chirp transforms. Therefore, an accurate,
simple, and efficient numerical method will be beneficial
in the use of these extensively employed formulas.
Some studies have been done on the numerical evalua-
tion of Fresnel diffraction; two examples are the use of
fast Fourier transforms (FFT’s) and convolution
techniques
14
and a discrete-Fourier-transform- (DFT-)
like matrix method.
15,16
A number of algorithms are also
specifically devised for FrFT’s
17–21
whose properties have
been intensively investigated both mathematically
3–5,22
and physically
23–25
in terms of their applications to opti-
cal beam propagation,
6
imaging,
24–26
diffraction,
27
and
signal and image processing.
28–36
The relations among
FnT’s, DFT’s, and FrFT’s are also well established.
32,37,38
Our aim in this paper is to develop an efficient algorithm
that unifies the evaluations of these formulas.
In this paper we describe a fast numerical algorithm
that is based on the chirp-z transform
39,40
to calculate
chirp transforms. It employs two FFT’s with an analyti-
cal kernel, and its computational complexity is better
than a fast convolution. In addition, one can freely
choose the sampling resolutions in both x space (signal
domain) and u space (transformed, or response domain)
and zoom in on any portion of the data of interest. Zoom-
ing in may be very useful in studying fine structures of
some chirp systems, for example, near-field diffraction.
41
The sampling condition will also be addressed under the
restriction of the Nyquist theorem.
In Section 2, first we compare several physical systems
that can be classified as chirp transforms and establish a
general mathematical description. Then in Subsection
3.A we present a discrete form of the transform. The
sampling condition and discussions on peeping any por-
tion of the data in u space are mentioned at the end of the
Subsection 3.A. Following the discretization of the inte-
gral form of the transform, we use similar techniques in
the chirp-z transform and develop the concrete fast algo-
rithm in Subsections 3.B and 3.C. In Section 4 we give
some numerical examples to demonstrate the effective-
ness of this algorithm. Some closed-form transforms
such as a Gaussian function and rect(x/a) are tested in
Subsection 4.A. In Subsection 4.B the algorithm is ap-
plied to zooming in on the Fresnel diffraction of a rectan-
gular window around the focal plane. Experimental re-
sults are provided for comparison. Near-field diffraction
patterns from a 1-to-48-waveguide fan-out interconnec-
tion layer were measured, and they corresponded well to
762 J. Opt. Soc. Am. A/ Vol. 17, No. 4 / April 2000 Deng et al.
0740-3232/2000/040762-10$15.00 © 2000 Optical Society of America
the results simulated with this algorithm in Subsection
4.C. In Subsection 4.D the versatility of this algorithm is
demonstrated in calculating arbitrary real-order FrFT’s.
2. COMPARISON OF SEVERAL GENERAL
CHIRP TRANSFORM SYSTEMS
A chirp transform of an arbitrary signal f(x)inan
N-dimensional system is expressed by
g
ch
p
共
u
兲
⬅ C
p
兵
f
共
x
兲
其
⫽
冕
⫺⬁
⫹⬁
f
共
x
兲
B
ch
p
共
u,x
兲
dx, (1)
and the integration kernel is
B
ch
p
共
u, x
兲
⬅ A
共
p
兲
exp
冋
⫺2
i
兺
l⫽1
N
共
␣
l
x
l
2
⫹

l
x
l
u
l
⫹
␥
l
u
l
2
兲
册
,
u, x 苸 R
N
, (2)
where A(p) is the amplitude of the kernel and generally a
complex value, i ⬅
冑
⫺1, and p ⬅ (
␣
,

,
␥
) are the pa-
rameters of the transform, each of them being an
N-dimensional vector. This transform maps the signal in
x space (signal domain) to u space (chirp space, or re-
sponse domain). All variables are dimensionless. To
simplify the discussion, we will use the one-dimensional
form of the above equations from now on. However, all
the conclusions and results are applicable to the higher-
dimensional cases as well.
In comparison, the Fresnel formula of scalar diffraction
through a general Gaussian (first-order) optical system
can be expressed with the Collins formula,
2,6,42
in which
the transforming kernel is a special case of chirp trans-
forms [Eq. (1)], namely,
B
FnT
M
共
u, x
兲
⬅
⫺i
B
exp
冋
ik
2B
共
Ax
2
⫺ 2xu ⫹ Du
2
兲
册
, (3)
which is associated with a ray-transfer matrix,
2,42,43
M ⬅
冋
AB
CD
册
, det M ⫽⫾1. (4)
In Eq. (3), ⬅(2
)/k is the wavelength. This formula
connects the paraxial parameters and the Fresnel diffrac-
tion of the system. It can be rewritten as a special case of
chirp transforms when the corresponding parameters are
used:
p ⬅
共
␣
,

,
␥
兲
⫽
⫺k
4
B
共
A,⫺2, D
兲
,
A
共
p
兲
⫽
⫺i
B
. (5)
Similarly, the kernel of a real-value
th-order FrFT,
B
FrFT
, is another special case of the chirp transform and
can be written as
44
B
FrFT
共
u, x
兲
⫽
exp
关
⫺i
共
s⫺
兲
/2
兴
共
2
兩
sin
共
兲
兩
兲
1/2
⫻ exp
再
i
冋
共
u
2
⫹ x
2
兲
2
⫻ cot
⫺ ux csc
册
冎
,
⫽ 2n, (6)
in which
⬅ (
/2)
,
s
⫽
关
(
/2)sign(sin
)
兴
and n is an
integer.
45
Apparently, the corresponding ChT’s param-
eters are
p ⬅
共
␣
,

,
␥
兲
⫽
⫺1
4
关
cot
共
兲
, ⫺2 csc
共
兲
, cot
共
兲
兴
,
A
共
p
兲
⫽
exp
关
⫺i
共
s⫺
兲
/2
兴
共
2
兩
sin
共
兲
兩
兲
1/2
. (7)
Therefore a general algorithm for chirp transform applies
to computations of all these cases of Fresnel diffraction,
the Collins formula, and FrFT’s. To preset a suitable
form for numerical calculation of any ChT, in Subsection
3.A we will derive a discrete form of Eq. (1), and in Sub-
section 3.B we will concentrate on developing the fast al-
gorithm.
3. DEVELOPMENT OF A FAST NUMERICAL
ALGORITHM FOR THE GENERAL
CHIRP TRANSFORM
A. Discrete Form of the Chirp Transform
For numerical calculation of the ChT of an arbitrary func-
tion f(x) except those having analytical forms, first Eq. (1)
will be digitized. For simplicity, we will use the simplest
equidistant sampling. The sample steps are denoted
␦
x
for the signal space and
␦
u for the chirp space. The
sample numbers N
x
and N
u
for x space and u space must
be limited; they are not necessarily equal. Therefore the
discretization of Eq. (1) may be given as
x → n
⬘
␦
x ⬅
冉
n ⫺
N
x
2
冊
␦
x, n ⫽ 0, 1, 2,... N
x
⫺ 1,
u → k
⬘
␦
u ⬅
冉
k ⫺
N
u
2
冊
␦
u, k ⫽ 0, 1, 2,... N
u
⫺ 1,
Fig. 1. Schematic top view of the H-tree waveguide used in op-
tical interconnections. The side view of one of its 48 microcou-
plers is illustrated at the right. Near-field diffraction patterns
of the outcoupling are critical to the coupling of the optical signal
to detectors in successive layers and hence to the performance of
the optical interconnections.
Deng et al. Vol. 17, No. 4/ April 2000 / J. Opt. Soc. Am. A 763
and
g
k
p
⫽ g
ch
p
共
k
␦
u
兲
⬇
兺
n⫽0
N
x
⫺1
f
n
B
D
p
共
k, n
兲
␦
x,
B
D
p
共
k, n
兲
⬅ A
共
p
兲
exp
兵
⫺2
i
关
␣
共
n
⬘
␦
x
兲
2
⫹

n
⬘
k
⬘
␦
x
␦
u
⫹
␥
共
k
⬘
␦
u
兲
2
兴
其
, (8)
where f
n
⫽ f(n
⬘
␦
x) and B
D
p
(k, n) ⫽ B
ch
p
(k
⬘
␦
u, n
⬘
␦
x). It
can be seen that the parameterized ChT of f(x) is equiva-
lent to the Fourier transform of the modified function
f
m
(x) ⬅ f(x)exp(⫺2
i
␣
x
2
), with a simultaneous variable
scaling in u space except for the difference of a complex
phasor. This could be inferred from Eq. (2) and written
as
g
ch
p
共
u
兲
⬀ g
共
v
兲
⬅ F
兵
f
m
共
x
兲
其
⫽ F
兵
f
共
x
兲
exp
共
⫺2
i
␣
x
2
兲
其
,
(9)
where v ⫽

u and F denotes a Fourier transform.
46
This is a well-known technique and was used in deriving
several properties of FrFT (for example, see Refs. 13, 37,
and 38). Therefore, according to the Nyquist sampling
theorem, if f
m
(x) is band limited (i.e., there exists a mini-
mum value v
m
, g(v) ⫽ 0,
兩
v
兩
⭓ v
m
⬎
0), the sampling
step in x space must satisfy
␦
x ⭐
1
2v
m
⬅
1
2

u
m
. (10)
This is one implicit form of the Nyquist sampling theorem
in ChT spaces.
34,47
Another explicit form of the sampling
condition will be derived in Subsection 3.B. Under this
restriction, one may correctly calculate the numerical
transform of a given function.
If one is interested in the data peeping in u space and
wants to scrutinize any region of g
ch
p
(u), it may be in-
structive to look at the shifting rule of Fourier trans-
forms. For any nontrivial ChT,

⫽ 0, and any inter-
ested data window center c in u space,
g
ch
p
共
u ⫺ c
兲
⬀ g
共
v ⫺

c
兲
⫽ F
兵
f
m
共
x
兲
exp
关
2
i

cx
兴
其
. (11)
The implementation for the data peeping is straightfor-
ward and needs no further discussion. In Subsection 3.B
we will present the concrete fast procedure for numerical
evaluation of the ChT.
B. Fast Numerical Algorithm for the Chirp Transform
The techniques used in the chirp-z transform
39,40
are also
useful for efficient calculation of the ChT of a given func-
tion f(x). Substituting the following expression into Eq.
(8),
n
⬘
k
⬘
⫽⫺
1
2
关共
n
⬘
⫺ k
⬘
兲
2
⫺ n
⬘
2
⫺ k
⬘
2
兴
, (12)
we get
B
D
p
共
k, n
兲
⫽ A
共
p
兲
P
n
p
Q
k
p
B
convol
p
共
n
⬘
⫺ k
⬘
兲
,
where the phasorlike P
n
p
and Q
k
p
, as well as the modified
kernel B
convol
p
(n
⬘
⫺ k
⬘
) are
P
n
p
⫽ exp
冋
⫺2
i
冉
␣
␦
x
2
⫺

2
␦
x
␦
u
冊
共
n
⬘
兲
2
册
,
Q
k
p
⫽ exp
冋
⫺2
i
冉
␥
␦
u
2
⫺

2
␦
x
␦
u
冊
共
k
⬘
兲
2
册
,
(13)
B
convol
p
共
n
⬘
⫺ k
⬘
兲
⫽ exp
关
i

共
n
⬘
⫺ k
⬘
兲
2
␦
x
␦
u
兴
.
Substituting Eq. (13) into Eq. (8), one can obtain
g
k
p
⫽ A
共
p
兲
␦
x
冋
兺
n⫽0
N
x
⫺1
f
n
共
m
兲
B
convol
p
共
n
⬘
⫺ k
⬘
兲
册
Q
k
p
, (14)
where f
n
(m)
⫽ f
n
P
n
p
is the modified discrete signal.
Therefore the discrete ChT of f
n
can be efficiently calcu-
lated through a fast convolution algorithm by use of FFT’s
such as the techniques in Refs. 17 and 39 which are
straightforward and simple to implement. We can write
the procedure symbolically as
g
k
p
⫽ A
共
p
兲
兵
F
⫺1
兵
F
兵
f
n
共
m
兲
其
F
兵
␦
xB
convol
p
共
n
⬘
兲
其其
Q
k
p
. (15)
F denotes a FFT and F
⫺1
an inverse FFT. However, this
is not the most efficient method. The transforming ker-
nel is actually a generalized Gaussian function whose
Fourier transform has a closed form,
48
i.e.,
F
兵
exp
共
⫺ax
2
兲
其
⫽
冑
a
exp
冋
⫺
共
u
兲
2
a
册
, R
兵
a
其
⭓ 0,
(16)
which, when applied to Eq. (15), becomes
B
ˆ
convol
p
共
k
兲
⫽ F
兵
␦
xB
convol
p
共
n
⬘
兲
其
⫽
冑
i
␦
x

␦
u
exp
冉
⫺
i
u
2

␦
x
␦
u
冊
, R
再
i

冎
⭓ 0.
(17)
For the discrete transform, suppose that the sampling
number in x space is L
x
⭓ N
x
and the Fourier transform
of the kernel on the u grids is
B
ˆ
convol
p
共
k
兲
⫽
冑
i
␦
x

␦
u
exp
冋
⫺
i
1

␦
x
␦
u
冉
k
⬘
L
x
冊
2
册
,
k
⬘
⬅ k ⫺
L
x
2
, k ⫽ 0, 1, 2,..., L
x
⫺ 1. (18)
The discrete ChT of Eq. (15) can be rewritten as
g
k
p
⫽ A
共
p
兲
兵
F
⫺1
兵
F
兵
f
n
共
m
兲
其
B
ˆ
convol
p
共
k
兲
其
Q
k
p
其
. (19)
For implementation of this algorithm, the sampling of the
B
ˆ
convol
p
(k)inu space must be so dense that the kernel is
well approximated. In practice, this condition may be
satisfied by
1

␦
x
␦
u
冉
1
L
x
冊
2
Ⰶ 1. (20)
It has been assumed that the FFT’s of the functions in-
volved do exist. However, these functions can be a dis-
tribution (for example, the Dirac delta function) if we take
the Fourier transforms in a general sense. Therefore one
need not resort to mathematically rigorous conditions of
the fast algorithm in practice by abiding by the sampling
condition mentioned above. An alternative is to explic-
itly impose the restrictions of Eq. (17) on the parameter

.
764 J. Opt. Soc. Am. A/ Vol. 17, No. 4 / April 2000 Deng et al.
One can easily verify that all the Fourier transforms are
valid if the chirp transform in Eq. (1) exists.
The complexity of numerically evaluating Eq. (19) is
about the same as that of two FFT’s and no more than six
sequential complex multiplications (including the possible
data-peeping operation), namely, 2L
x
关
log(L
x
) ⫹ O(1)
兴
(L
x
is the number of FFT samples). Compared with the pre-
vious chirp-z techniques,
39,40
such as those used in Ref.
17, whose time complexity is exactly as that of Eq. (15),
namely, 3L
x
关
log(L
x
) ⫹ O(1)
兴
, the improvement in compu-
tation time is about one third. Besides, with an analyti-
cal form, the accuracy will be much better. Some ex-
amples will be presented in Section 4. By using this
algorithm, one can freely set the sampling steps in both x
and u space, which is very useful under some circum-
stances such as examining near-field diffraction patterns
and interpolation of sparse data in u space.
15,39
One can
also choose different numbers of samples in x and u space
by means of the implementation techniques in Subsection
3.C.
C. Implementation of the Fast Algorithm
To implement the algorithm, the most convenient and ef-
ficient way is probably to use the standard FFT’s. Differ-
ent numbers of samples in x and u space will be used in
the following procedure.
Suppose that the original number of samples in x space
is N
x
and the required sampling number in u space is N
u
,
which can be different from N
x
. Thus the length of a
noncircular convolution is L
convol
⫽ N
x
⫹ N
u
⫺ 1, which
is the minimum length that should be used for a FFT in
Eq. (19). The sequence of calculations is summarized be-
low.
1. Find a positive number L
x
⭓ L
convol
that makes
(L
x
⫺ N
x
) an even number. L
x
usually can take the low-
est value that satisfies this requirement.
2. Then compute f
n
and zero-pad f
n
(m)
:
Z
兵
f
n
共
m
兲
其
共
n
兲
⫽
冦
0
0 ⭐ n
⬍
共
L
x
⫺ N
x
兲
2
f
n⫺
共
L
x
⫺N
x
兲
/2
共
m
兲
共
L
x
⫺ N
x
兲
2
⭐ n
⬍
共
L
x
⫹ N
x
兲
2
0
共
L
x
⫹ N
x
兲
2
⭐ n ⭐ L
x
⫺ 1
.
(21)
3. Compute the analytical kernel B
ˆ
convol
p
(k) by means
of Eq. (18) for k
⬘
⫽ k ⫺ L
x
/2, k ⫽ 0, 1, 2,..., L
x
⫺ 1.
4. Perform the fast convolution operation in Eq. (19):
C
共
l
兲
⫽ F
⫺
兵
Z
兵
f
n
共
m
兲
其其
B
ˆ
convol
p
共
k
兲
其
,
l ⫽ 0, 1, 2,..., L
x
⫺ 1. (22)
5. Discard the first and last (L
x
⫺ N
u
)/2 data of C(l),
and for 0 ⭐ k ⭐ N
u
⫺ 1 let
g
k
p
⫽ A
共
p
兲
Q
k
p
C
冉
k ⫹
L
x
⫺ N
u
2
冊
, (23)
in which (L
x
⫺ N
u
) should be an even number. This
constraint is not necessary, though it could be met easily
in practice.
Similar procedures with minor modifications can be de-
veloped for arbitrary N
u
. To simplify our discussion be-
low, the equivalent requirement that N
x
and N
u
be of the
same parity will be used. The performance of the algo-
rithm will be evaluated with some application examples
in Section 4.
4. APPLICATION EXAMPLES
To demonstrate the performance of this fast algorithm,
we examine several chirp systems. Different samplings
and resolutions are used for these systems. In addition,
to unify our discussions, we define a scaling factor,
zoom
⫽
␦
u
␦
x
,
zoom
⬎
0, (24)
which defines the finesse of the zooming in. Its recipro-
cal, 1/
zoom
, is the magnification of the data window. For
applications of data peeping or interpolation, 0
⬍
zoom
⭐ 1, whereas zoom-out effects occur when
zoom
⬎
1 un-
der proper sampling; this topic will be addressed below.
Fig. 2. (a) Errors e
p
⫽ max
关
e(x)
兴
of the ChT’s of a Gaussian
function compared with the analytical transforms for different
zoom factors
zoom
. (b) Typical error curve e(x)(
zoom
⫽ 0.5380).
Deng et al. Vol. 17, No. 4/ April 2000 / J. Opt. Soc. Am. A 765
A. Analytical Chirp Transforms
There are many functions whose ChT’s are analytically
available, such as a generalized Gaussian in Eq. (16).
For simplicity,
␣
⫽
␥
⫽ 0 and R
x
⫽
冑
N
x
were used. The
ChT becomes a Fourier transform, and its space–
bandwidth product
47,49
is N
x
. We tested the fast algo-
rithm with a Gaussian function for which a
2
⫽
冑
N
x
/(2R
x
), where x 苸
关
⫺R
x
/2,R
x
/2
兴
, N
x
⫽ N
u
⫽ 512, and L
x
⫽ N
x
⫹ N
u
⫽ 1024. A normalized error
measurement
e ⫽
兩
g
k
p
⫺ F
兵
exp
共
⫺ax
2
兲
其
兩
max
共
兩
g
k
p
兩
兲
(25)
between the numerical result and the analytical one is
used to check the accuracy. For different values of
zoom
,
the ChT’s and the corresponding e
p
⬅ max(e) are plotted
in Fig. 2(a), and a typical error curve for the ChT is plot-
ted in Fig. 2(b). For 0.01 ⭐
zoom
⬍
2.0, corresponding to
a magnification of the data window of 0.5 ⬃ 100,
e
p
⬇ 10
⫺12
.
A square function, rect(x/a), was used to test the algo-
rithm. Its Fourier transform also has the closed form,
F
兵
rect(x/a)
其
⫽ sin(2
au)/
u. For the same scaling fac-
tors, a ⫽ R
x
/4, the ChT’s as well as the accuracy mea-
surement e
p
are presented in Fig. 3(a). A typical error
curve for the ChT is plotted in Fig. 3(b). The accuracy of
the transform is much lower than that for a Gaussian,
e ⬇ 10
⫺3
. However, this does not decrease the validity
and performance of the algorithm. Comparisons of Fig.
2(b), Fig. 3(b), and Fig. 4 indicate that most of the errors
are accumulated in the FFT’s used in the calculations.
Therefore one can optimize the FFT parameters to
achieve higher accuracy.
50
B. Fresnel Diffraction of a Rectangular Window
The diffraction of a first-order optical system is perhaps
one of the simplest chirp transforms but is not a trivial
chirp one. In this section we calculate the well-known
near-field diffraction of a rectangular window,
rect(x/a, y/b), behind an ideal focusing lens at different
distances. The system is illustrated in Fig. 5(a), whose
space–bandwidth product is discussed in Ref. 49. The el-
ements of its ray-transfer matrix M are A ⫽ 1 ⫺ l
2
/f, B
⫽ l
1
A ⫹ l
2
, C ⫽⫺1/f, and D ⫽ 1 ⫺ l
1
/f. In the simu-
lations N
x
⫽ 2N
u
⫽ 1024, ⫽ 0.85
m, a ⫽ 400
m, b
⫽ 200
m, l
1
⫽ 1.0
⫺8
m, f ⫽ 5 ⫻ 10
3
m, and l
2
were
varied from l
1
to l
1
⫹ f while
zoom
changed correspond-
ingly.
Fig. 3. (a) Errors max(e) of the ChT’s of a rect(x/a) function
compared with the analytical transforms for different zoom fac-
tors
zoom
. (b) Typical error curve e(x)(
zoom
⫽ 0.5380).
Fig. 4. Typical error curve of the FFT for (a) a Gaussian func-
tion, (b) rect(x/a).
766 J. Opt. Soc. Am. A/ Vol. 17, No. 4 / April 2000 Deng et al.
