7
Fast Fourier Transform and Its
Applications
Frequency analysis of digital signals and systems was discussed in Chapter 4. To per-
form frequency analysis on a discrete-time signal, we converted the time-domain
sequence into the frequency-domain representation using the z-transform, the
discrete-time Fourier transform (DTFT), or the discrete Fourier transform (DFT). The
widespread application of the DFT to spectral analysis, fast convolution, and data
transmission is due to the development of the fast Fourier transform (FFT) algorithm
for its computation. The FFT algorithm allows a much more rapid computation of the
DFT, was developed in the mid-1960s by Cooley and Tukey.
It is critical to understand the advantages and the limitations of the DFT and how to
use it properly. We will discuss the important properties of the DFT in Section 7.1. The
development of FFT algorithms will be covered in Section 7.2. In Section 7.3, we will
introduce the applications of FFTs. Implementation considerations such as computa-
tional issues and finite-wordlength effects will be discussed in Section 7.4. Finally,
implementation of the FFT algorithm using the TMS320C55x for experimental purposes
will be given in Section 7.5.
7.1 Discrete Fourier Transform
As discussed in Chapter 4, we perform frequency analysis of a discrete-time signal x(n)
using the DTFT defined in (4.4.4). However, X! is a continuous function of frequency
! and the computation requires an infinite-length sequence x(n). Thus the DTFT cannot
be implemented on digital hardware. We define the DFT in Section 4.4.3 for N samples
of x(n)atN discrete frequencies. The input to the N-point DFT is a digital signal
containing N samples and the output is a discrete-frequency sequence containing N
samples. Therefore the DFT is a numerically computable transform and is suitable for
DSP implementations.
Real-Time Digital Signal Processing. Sen M Kuo, Bob H Lee
Copyright # 2001John Wiley & Sons Ltd
ISBNs: 0-470-84137-0 (Hardback); 0-470-84534-1 (Electronic)
7.1.1 Definitions


NÀ1
n0
x
p
ne
Àj2p=Nkn
, k  0, 1, ..., N À 1,
7:1:1a
where
x
p
n

I
lÀI
xn À lN7:1:1b
is a periodic signal with period N.
In general, the equally spaced frequency samples do not represent the original
spectrum X! uniquely when x(n) has infinite duration. Instead, these frequency
samples correspond to a periodic sequence x
p
n as shown in (7.1.1). When the sequence
x(n) has a finite length N, x
p
n is simply a periodic repetition of x(n). Therefore
the frequency samples X !
k
, k  0, 1, ..., N À 1uniquely represent the
finite-duration sequence x(n). Since xnx

FAST FOURIER TRANSFORM AND ITS APPLICATIONS
and
XN=2

NÀ1
n0
e
Àjpn
xn

NÀ1
n0
À1
n
xn:
The DFT defined in (7.1.2) can also be written as
Xk

NÀ1
n0
xnW
kn
N
, k  0, 1, ...N À 1, 7:1:3
where
W
kn
N
 e
Àj

successive powers W
k
N
, k  0, 1, ..., N À 1are also Nth roots of unity, but in clockwise
direction on the unit circle. It can be shown that W
N=2
N
 e
Àjp
À1, the symmetry
property
W
kN=2
N
ÀW
k
N
,0 k N=2 À 1 7:1:5a
and the periodicity property
W
kN
N
 W
k
N
: 7:1:5b
Figure 7.1illustrates the cyclic property of the twiddle factors for an eight-point DFT.
W
8
0

W
Figure 7.1 Twiddle factors for DFT, N  8 case
DISCRETE FOURIER TRANSFORM
305
The inverse discrete Fourier transform (IDFT) is used to transform the X(k) back into
the original sequence x(n). Given the frequency samples X(k), the IDFT is defined as
xn
1
N

NÀ1
k0
Xke
j2p=Nkn

1
N

NÀ1
k0
XkW
Àkn
N
, n  0, 1, ..., N À 1: 7:1:6
This is identical to the DFT with the exception of the normalizing factor 1/N and
the sign of the exponent of the twiddle factors. The IDFT shows that there is no loss
of information by transforming the spectrum X(k) back into the original time sequence
x(n). The DFT given in (7.1.3) is called the analysis transform since it analyzes the signal
x(n)atN frequency components. The IDFT defined in (7.1.6) is called the synthesis
transform because it reconstructs the signal x(n) from the frequency components.

N
1 À ae
Àj2pk=N
, k  0, 1, ..., N À 1:
The DFT and IDFT defined in (7.1.3) and (7.1.6), can be expressed in matrix±vector
form as
X  Wx 7:1:7a
and
x 
1
N
W
Ã
X , 7:1:7b
where x x0 x1 ...xN À 1
T
is the signal vector, the frequency-domain DFT
coefficients are contained in the complex vector X X0 X 1 ...XN À 1
T
, and
the NxN twiddle-factor matrix (or DFT matrix) W is given by
W  W
kn
N
ÂÃ
0 k
,
n NÀ1

11ÁÁÁ 1

T
T
T
R
Q
U
U
U
U
U
U
U
S
,
7:1:8
306
FAST FOURIER TRANSFORM AND ITS APPLICATIONS
and W
Ã
is the complex conjugate of the matrix W. Since W is a symmetric matrix, the
inverse matrix W
À1

1
N
W
Ã
was used to derive (7.1.7b).
Example 7.3: Given xnf1,1,0,0g, the DFT of this four-point sequence can be
computed using the matrix formulation as

T
T
T
T
T
T
R
Q
U
U
U
U
U
U
S
x

1111
1 Àj À1 j
1 À11À1
1 j À1 Àj
P
T
T
T
T
T
T
R
Q

0
1  j
P
T
T
T
T
T
T
R
Q
U
U
U
U
U
U
S
,
where we used symmetry and periodicity properties given in (7.1.5) to obtain
W
0
4
 W
4
4
 1, W
1
4
 W

À4
4
W
À6
4
1 W
À3
4
W
À6
4
W
À9
4
P
T
T
T
R
Q
U
U
U
S
X

1
4
1111
1 j À1 Àj

R
Q
U
U
S
:
As shown in Figure 7.1, the twiddle factors are equally spaced around the unit circle
at frequency intervals of f
s
=N (or 2p=N). Therefore the frequency samples X(k) repre-
sent discrete frequencies
f
k
 k
f
s
N
, k  0, 1, ..., N À 1: 7:1:9
The computational frequency resolution of the DFT is equal to the frequency increment
f
s
=N, and is sometimes referred to as the bin spacing of the DFT outputs. The spacing
DISCRETE FOURIER TRANSFORM
307
of the spectral lines depends on the number of data samples. This issue will be further
discussed in Section 7.3.2.
Since the DFT coefficient X(k) is a complex variable, it can be expressed in polar form
as
XkjXkje
jk

where a and b are arbitrary constants. Linearity is a key property that allows us to
compute the DFTs of several different signals and determine the combined DFT via the
summation of the individual DFTs. For example, the frequency response of a given
system can be easily evaluated at each frequency component. The results can then be
combined to determine the overall frequency response.
308
FAST FOURIER TRANSFORM AND ITS APPLICATIONS
Complex-conjugate property
If the sequence fxn,0 n N À 1g is real, then
XÀkX
Ã
kXN À k,0 k N À 1, 7:1:14
where X
Ã
k is the complex conjugate of X(k). Or equivalently,
XM  kX
Ã
M À k,0 k M, 7:1:15
where M  N=2ifN is even, and M N À 1=2ifN is odd. This property shows that
only the first (M  1) DFT coefficients are independent. Only the frequency compon-
ents from k  0tok  M are needed in order to completely define the output. The rest
can be obtained from the complex conjugate of corresponding coefficients, as illustrated
in Figure 7.2.
The complex-conjugate (or symmetry) property shows that
ReXk  ReXN À k, k  1,2, ..., M À 1 7:1:16
and
ImXk  ÀImXN À k, k  1,2, ..., M À 1: 7:1:17
Thus the DFT of a real sequence produces symmetric real frequency components and
anti-symmetric imaginary frequency components about X(M). The real part of the DFT
output is an even function, and the imaginary part of the DFT output is an odd

=2, and the entire range from Àf
s
=2tof
s
=2 was repeated infinitely in both
directions in the frequency domain. The DFT outputs represent a single period (from
0tof
s
) of the spectrum.
Circular shifts
Let fXkg be the DFT of a given N-periodic sequence fxng, and let y(n) be a circular
shifted sequence defined by
ynxn À m
mod N
, 7:1:21
where m is the number of samples by which x(n) is shifted to the right (or delayed) and
the modulo operation
 j
mod N
 j 
À
iN 7:1:22a
for some integer i such that
0 j
mod N
< N: 7:1:22b
For example, if m  1, xN À 1 replaces x0, x0 replaces x(1), x(1) replaces x(2), etc.
Thus a circular shift of an N-point sequence is equivalent to a linear shift of its periodic
extension.
For a given y(n) in (7.1.21), we have






ze
j2pk=N


I
nÀI
xne
Àj2pk=Nn
: 7:1:24
This is identical to evaluating the discrete-time Fourier transform X! at the N equally
spaced frequencies !
k
 2pk=N, k  0, 1, ..., N À 1. If the sequence x(n) has a finite
duration of length N, the DFT of a sequence yields its z-transform on the unit circle at a
set of points that are 2p=N radians apart, i.e.,
XkXzj
z  e
j
2p
N

k
, k  0, 1, ..., N À 1: 7:1:25
Therefore the DFT is equal to the z-transform of a sequence x(n) of length N, evaluated
at N equally spaced points on the unit circle in the z-plane.

T 1for k T iN, we have Xk0
for k  1,2, ..., N À 1. For k  0,

NÀ1
n0
W
kn
N
 N. Therefore we obtain
XkcNdk; k  0, 1, ..., N À 1:
7.1.3 Circular Convolution
The Fourier transform, the Laplace transform, and the z-transform of the linear con-
volution of two time functions are simply the products of the transforms of the
individual functions. A similar result holds for the DFT, but instead of a linear
convolution of two sequences, we have a circular convolution. If x(n) and h(n) are
real-valued N-periodic sequences, y(n) is the circular convolution of x(n) and h(n)
defined as
ynxnhn, n  0, 1, ..., N À 1, 7:1:26
DISCRETE FOURIER TRANSFORM
311
x(n)
x(n−1)
x(n−2)
x(n−N+1)
h(n)
h(n−1)
h(n−2)
h(n−N+1)
Figure 7.3 Circular convolution of two sequences using the concentric circle approach
where  denotes circular convolution. Then

Example 7.5: Given two 8-point sequences xnf1,1,1,1,1,0,0,0g and
hnf0,0,0,1,1,1,1,1g. Using the circular convolution method illustrated in
Figure 7.3, we can obtain
312
FAST FOURIER TRANSFORM AND ITS APPLICATIONS
n  0, y01  1  1  1  4
1
0
1
1
1
1
110
0
0
0
11
10
n  1, y11  1  1  3
1
0
1
1
1
0
111
0
0
0
11

Example 7.6: Consider the previous example. If these 8-point sequences h(n)and
x(n) are zero-padded to 16 points, the resulting circular convolution is
ynxnhnf0,0,0,1,2,3,4,5,4,3,2,1,0,0,0,0g:
This result is identical to the linear convolution of the two sequences. Thus the
linear convolution discussed in Chapter 5 can be realized by the circular convolu-
tion with proper zero padding.
In MATLAB, zero padding can be implemented using the function zeros. For
example, the 8-point sequence x(n) given in example 7.5 can be zero-padded to 16
points with the following command:
x  [1, 1, 1, 1, zeros(1, 11)];
where the MATLAB function zeros(1, N)generates a row vector of N zeros.
7.2 Fast Fourier Transforms
The DFT is a very effective method for determining the frequency spectrum of a time-
domain signal. The only drawback with this technique is the amount of computation
necessary to calculate the DFT coefficients X(k). To compute each X(k), we need
approximately N complex multiplications and N complex additions based on the DFT
defined in (7.1.3). Since we need to compute N samples of X(k) for k  0, 1, ..., N À 1,
a total of approximately N
2
complex multiplications and N
2
À N complex additions are
required. When a complex multiplication is carried out using digital hardware, it
requires four real multiplications and two real additions. Therefore the number of
arithmetic operations required to compute the DFT is proportional to 4N
2
, which
becomes very large for a large number of N. In addition, computing and storing the
twiddle factors W
kn

example, if N is a power of 2, then the FFT makes it possible to calculate the DFT with
N log
2
N operations instead of N
2
operations. For N  1024, FFT requires about 10
4
operations instead 10
6
of operations for DFT.
The generic term FFT covers many different algorithms with different features,
advantages, and disadvantages. Each FFT has different strengths and makes different
tradeoffs between code complexity, memory usage, and computation requirements. The
FFT algorithm introduced by Cooley and Tukey requires approximately N log
2
N
multiplications, where N is a power of 2. The FFT can also be applied to cases where
N is a power of an integer other than 2. In this chapter, we introduce FFT algorithms
for the case where N is a power of 2, the radix-2 FFT algorithm.
There are two classes of FFT algorithms: decimation-in-time and decimation-in-
frequency. In the decimation-in-time algorithm, the input time sequence is successively
divided up into smaller sequences, and the DFTs of these subsequences are combined in
a certain pattern to yield the required DFT of the entire sequence with fewer operations.
Since this algorithm was derived by separating the time-domain sequence into succes-
sively smaller sets, the resulting algorithm is referred to as a decimation-in-time algo-
rithm. In the decimation-in-frequency algorithm, the frequency samples of the DFT are
decomposed into smaller and smaller subsequences in a similar manner.
7.2.1 Decimation-in-Time
In the decimation-in-time algorithm, the N-sample sequence fxn, n  0, 1, ..., N À 1g
is first divided into two shorter interwoven sequences: the even numbered sequence

7:2:3
Since
W
2mk
N
 e
Àj
2p
N
2mk
 e
Àj
2p
N=2
mk
 W
mk
N=2
, 7:2:4
Equation (7.2.3) can be written as
Xk

N=2À1
m0
x
1
mW
mk
N=2
 W

X
2
k, k  0, 1, ..., N=2À1
X
1
kÀW
k
N
X
2
k, k  N=2, ..., N À 1,
@
7:2:6
where X
1
kDFTx
1
m and X
2
kDFTx
2
m using the N=2-point DFT.
The important point about this result is that the DFT of N samples becomes a linear
combination of two smaller DFTs, each of N=2 samples. This procedure is illustrated in
Figure 7.4 for the case N  8. The computation of X
1
k and X
2
k requires 2N=2
2

x(2)
x(4)
x(6)
x(1)
x(3)
x(5)
x(7)
N/2-point
DFT
N/2-point
DFT
X
1
(0)
X
1
(1)
X
1
(2)
X
1
(3)
X
2
(0)
X
2
(1)
X

(m−1)th
stage
mth
stage
Figure 7.5 Flow graph for a butterfly computation
316
FAST FOURIER TRANSFORM AND ITS APPLICATIONS
lower half. Each butterfly involves just a single complex multiplication by a twiddle
factor W
k
N
, one addition, and one subtraction. For this first decomposition, the twiddle
factors are indexed consecutively, and the butterfly values are separated by N/2 samples.
The order of the input samples has also been rearranged (split between even and odd
numbers), which will be discussed in detail later.
Since N is a power of 2, N=2 is even. Each of these N=2-point DFTs in (7.2.6) can be
computed via two smaller N=4-point DFTs, and so on. The second step process is
illustrated in Figure 7.6. Note that the order of the input samples has been rearranged
as x(0), x(4), x(2), and x(6) because x(0), x(2), x(4), and x(6) are considered to be the 0th,
1st, 2nd, and 3rd inputs in a 4-point DFT. Similarly, the order of x(1), x(5), x(3), and x(7)
is used in the second 4-point DFT.
By repeating the process associated with (7.2.6), we will finally end up with a set of 2-
point DFTs since N is a power of 2. For example, in Figure 7.6, the N/4-point DFT
became a 2-point DFT since N  8. Since the twiddle factor for the first stage, W
0
N
 1,
the 2-point DFT requires only one addition and one subtraction. The 2-point DFT
illustrated in Figure 7.7 is identical to the butterfly network.
Example 7.7: Consider the two-point DFT algorithm which has two input time-

N/4-point
DFT
N/4-point
DFT
N/4-point
DFT
N/4-point
DFT
X(0)
X(1)
X(2)
X(3)
X(4)
X(5)
X(6)
X(7)
W
8
0
W
8
1
W
8
2
W
8
3
W
8

sequence has the unusual order. Actually the order of the input sequence is arranged as
if each index was written in binary form and then the order of binary digits was reversed.
The bit-reversal process is illustrated in Table 7.1for the case N  8. Each of the time
sample indices in decimal is converted to its binary representation. The binary bit streams
are then reversed. Converting the reversed binary numbers to decimal values gives the
reordered time indices. If the input is in natural order the output will be in bit-reversed
order. We can either shuffle the input sequence with a bit-reversal algorithm to get the
output sequence in natural order, or let the input sequence be in natural order and shuffle
the bit-reversed results to obtain the output in natural order. Note that most modern
DSP chips such as the TMS320C55x provide the bit-reversal addressing mode to support
this bit-reversal process. Therefore the input sequence can be stored in memory with the
bit-reversed addresses computed by the hardware.
For the FFT algorithm shown in Figure 7.6, once all the values for a particular stage
are computed, the old values that were used to obtain them are never required again.
Thus the FFT needs to store only the N complex values. The memory locations used for
the FFT outputs are the same as the memory locations used for storing the input data.
This observation is used to produce in-place FFT algorithms that use the same memory
locations for input samples, all intermediate calculations, and final output numbers.
Table 7.1 Example of bit-reversal process, N  8 (3-bit)
Input sample index
&
Bit-reversed sample index
Decimal Binary Binary Decimal
0 000 000 0
1001100 4
20100102
3 011 110 6
4 100 001 1
5 101 101 5
6 110 011 3

N=2k
N

N=2À1
n0
xn
N
2

W
nk
N
:
7:2:7
Using the fact that
W
N=2
N
 e
Àj
2p
N
N=2
 e
Àjp
À1, 7:2:8
Equation (7.2.7) can be simplified to
Xk

N=2À1

7:2:10a
and
X2k  1

N=2À1
n0
xnÀxn
N
2
 !
W
k
N
W
kn
N=2
7:2:10b
for 0 k N=2À1. Let x
1
nxnxn
N
2
ÀÁ
and x
2
nxnÀxn
N
2
ÀÁ
for

x
1
(1)
x
1
(2)
x
1
(3)
x
2
(0)
x
2
(1)
x
2
(2)
x
2
(3)
X(0)
X(2)
X(4)
X(6)
X(1)
X(3)
X(5)
X(7)
W

1
N

NÀ1
k0
X
Ã
kW
kn
N
, n  0, 1, ..., N À 1: 7:2:11
Therefore an FFT algorithm can be used to compute the inverse DFT by first con-
jugating the DFT coefficients X(k) to obtain X
Ã
k, computing the DFT of X
Ã
k use an
FFT algorithm, scaling the results by 1/N to obtain x
Ã
n, and then complex conjugating
x
Ã
n to obtain the output sequence x(n). If the signal is real-valued, the final conjuga-
tion operation is not required.
All the FFT algorithms introduced in this chapter are based on two-input, two-
output butterfly computations, and are classified as radix-2 complex FFT algorithms.
It is possible to use other radix values to develop FFT algorithms. However, these
algorithms do not work well when the length is a number with few factors. In addition,
these algorithms are more complicated than the radix-2 FFT algorithms, and the
320