We have seen in the preceding sections that the DFT
is a very computationally intensive
operation. In 1965, Cooley and Tukey ([1])
published an algorithm that could
be used to compute the DFT much more efficiently.
Various forms of their algorithm,
which came to be known as the fast Fourier transform (FFT), had actually
been developed much earlier by other mathematicians (even dating back to
Gauss). It was their paper, however, which stimulated a revolution in the field
of signal processing.
It is important to keep in mind at the outset that the FFT is not
a new transform.
It is simply a very efficient way to compute an existing
transform, namely the DFT.
As we saw, a straight forward implementation of the DFT can
be computationally expensive because the number of multiplies
grows as the square of the input length (i.e. N2N2 for an NN point DFT).
The FFT reduces this computation using two simple but important
concepts.
The first concept, known as divideandconquer,
splits the problem into two smaller problems.
The second concept, known as recursion, applies this divideandconquer
method repeatedly until the problem is solved.
Consider the defining equation for the DFT and assume that NN is even,
so that N/2N/2 is an integer:
X
(
k
)
=
∑
n
=
0
N

1
x
(
n
)
e

j
2
π
k
n
/
N
.
X
(
k
)
=
∑
n
=
0
N

1
x
(
n
)
e

j
2
π
k
n
/
N
.
(4)Here we have dropped the subscript of NN in the notation
for X(k)X(k).
We will also use the notation
X
(
k
)
=
DFT
N
x
(
n
)
X
(
k
)
=
DFT
N
x
(
n
)
(5)to denote the NN point DFT of the signal x(n)x(n).
Suppose we break the sum in Equation 4 into two sums, one containing
all the terms for which nn is even, and one containing all the terms for
which nn is odd:
X
(
k
)
=
∑
n
=
0
n
even
N

1
x
(
n
)
e

j
2
π
k
n
/
N
+
∑
n
=
0
n
odd
N

1
x
(
n
)
e

j
2
π
k
n
/
N
.
X
(
k
)
=
∑
n
=
0
n
even
N

1
x
(
n
)
e

j
2
π
k
n
/
N
+
∑
n
=
0
n
odd
N

1
x
(
n
)
e

j
2
π
k
n
/
N
.
(6)We can eliminate the conditions “nn even” and “nn odd” in
Equation 6 by making a
change of variable in each sum. In the first sum, we replace nn by 2m2m.
Then as we sum mm from 0 to N/21N/21, n=2mn=2m will take on all even
integer values between 0 and N2N2. Similarly, in the second sum, we
replace nn by 2m+12m+1. Then as we sum mm from 0 to N/21N/21,
n=2m+1n=2m+1 will take on all odd
integer values between 0 and N1N1.
Thus, we can write
X
(
k
)
=
∑
m
=
0
N
/
2

1
x
(
2
m
)
e

j
2
π
k
2
m
/
N
+
∑
m
=
0
N
/
2

1
x
(
2
m
+
1
)
e

j
2
π
k
(
2
m
+
1
)
/
N
.
X
(
k
)
=
∑
m
=
0
N
/
2

1
x
(
2
m
)
e

j
2
π
k
2
m
/
N
+
∑
m
=
0
N
/
2

1
x
(
2
m
+
1
)
e

j
2
π
k
(
2
m
+
1
)
/
N
.
(7)Next we rearrange the exponent of the complex exponential in the first
sum, and split and rearrange the exponent in the second sum to yield
X
(
k
)
=
∑
m
=
0
N
/
2

1
x
(
2
m
)
e

j
2
π
k
m
/
(
N
/
2
)
+
e

j
2
π
k
/
N
∑
m
=
0
N
/
2

1
x
(
2
m
+
1
)
e

j
2
π
k
m
/
(
N
/
2
)
.
X
(
k
)
=
∑
m
=
0
N
/
2

1
x
(
2
m
)
e

j
2
π
k
m
/
(
N
/
2
)
+
e

j
2
π
k
/
N
∑
m
=
0
N
/
2

1
x
(
2
m
+
1
)
e

j
2
π
k
m
/
(
N
/
2
)
.
(8)Now compare the first sum in Equation 8 with the definition
for the DFT given by Equation 4. They have exactly the same form if
we replace NN everywhere in Equation 4 by N/2N/2. Thus the first sum
in Equation 8 is an N/2N/2 point DFT of the evennumbered
data points in the original sequence. Similarly, the second sum
in Equation 8 is an N/2N/2 point DFT of the oddnumbered
data points in the original sequence. To obtain the NN point DFT of the
complete sequence, we multiply the DFT of the oddnumbered data points
by the complex exponential factor ej2πk/Nej2πk/N, and then simply
sum the two N/2N/2 point DFTs.
To summarize, we will rewrite Equation 8 according to this
interpretation.
First, we define two new N/2N/2 point data sequences x0(n)x0(n) and
x1(n)x1(n), which contain the even and oddnumbered data
points from the original NN point sequence, respectively:
x
0
(
n
)
=
x
(
2
n
)
x
1
(
n
)
=
x
(
2
n
+
1
)
,
x
0
(
n
)
=
x
(
2
n
)
x
1
(
n
)
=
x
(
2
n
+
1
)
,
(9)where n=0,...,N/21n=0,...,N/21.
This separation of even and odd points is called decimation in time.
The NN point DFT of x(n)x(n) is then given by
X
(
k
)
=
X
0
(
k
)
+
e

j
2
π
k
/
N
X
1
(
k
)
for
k
=
0
,
.
.
.
,
N

1
.
X
(
k
)
=
X
0
(
k
)
+
e

j
2
π
k
/
N
X
1
(
k
)
for
k
=
0
,
.
.
.
,
N

1
.
(10)where X0(k)X0(k) and X1(k)X1(k) are the N/2N/2 point DFT's of the even and odd points.
X
0
(
k
)
=
DFT
N
/
2
[
x
0
(
n
)
]
X
1
(
k
)
=
DFT
N
/
2
[
x
1
(
n
)
]
X
0
(
k
)
=
DFT
N
/
2
[
x
0
(
n
)
]
X
1
(
k
)
=
DFT
N
/
2
[
x
1
(
n
)
]
(11)While Equation 10 requires less computation
than the original NN point DFT,
it can still be further simplified.
First, note that each N/2N/2 point DFT is periodic
with period N/2N/2.
This means that we need to only compute X0(k)X0(k) and X1(k)X1(k)
for N/2N/2 values of kk rather than the NN values
shown in Equation 10.
Furthermore, the complex exponential factor ej2πk/Nej2πk/N
has the property that

e

j
2
π
k
N
=
e

j
2
π
k
+
N
/
2
N
.

e

j
2
π
k
N
=
e

j
2
π
k
+
N
/
2
N
.
(12)These two facts may be combined to yield a simpler
expression for the NN point DFT:
X
(
k
)
=
X
0
(
k
)
+
W
N
k
X
1
(
k
)
X
(
k
+
N
/
2
)
=
X
0
(
k
)

W
N
k
X
1
(
k
)
for
k
=
0
,
.
.
.
,
N
/
2

1
X
(
k
)
=
X
0
(
k
)
+
W
N
k
X
1
(
k
)
X
(
k
+
N
/
2
)
=
X
0
(
k
)

W
N
k
X
1
(
k
)
for
k
=
0
,
.
.
.
,
N
/
2

1
(13)where the complex constants defined by WNk=ej2πk/NWNk=ej2πk/N
are commonly known as the twiddle factors.
Figure 1
shows a graphical interpretation of Equation 13
which we will refer to as the “divideandconquer DFT”.
We start on the left side with the data separated into
even and odd subsets.
We perform an N/2N/2 point DFT on each subset,
and then multiply the output of the odd DFT by the required twiddle
factors.
The first half of the output is computed by adding
the two branches,
while the second half is formed by subtraction.
This type of flow diagram is conventionally used to describe a
fast Fourier transform algorithm.
In this section, you will implement the DFT transformation
using Equation 13 and the illustration
in Figure 1.
Write a Matlab function with the syntax
X = dcDFT(x)
where x
is a vector of even length NN,
and X
is its DFT.
Your function dcDFT
should do the following:
 Separate the samples of x
into even and odd points.
The Matlab command
x0 = x(1:2:N)
can be used to obtain the “even” points.
 Use your function
DFTsum
to compute the two
N/2N/2 point DFT's.
 Multiply by the twiddle factors WNk=ej2πk/NWNk=ej2πk/N.
 Combine the two DFT's to form XX.
Test your function dcDFT
by using it to compute the DFT's
of the following signals.
 x(n)=δ(n)x(n)=δ(n) for N=10N=10
 x(n)=1forN=10x(n)=1forN=10
 x(n)=ej2πn/Nx(n)=ej2πn/N for N=10N=10

Submit the code for your function
dcDFT
.
 Determine the number of multiplies that are required in this approach
to computing an NN point DFT.
(Consider a multiply to be one multiplication of
real or complex numbers.)
Refer to the diagram of
Figure 1, and remember
to consider the
N/2N/2 point DFTs.
The second basic concept underlying the FFT algorithm is that of recursion.
Suppose N/2N/2 is also even.
Then we may apply the same decimationintime
idea to the computation of each of the N/2N/2 point DFT's
in Figure 1.
This yields the process depicted in Figure 2.
We now have two stages of twiddle factors
instead of one.
How many times can we repeat the process of decimating the input sequence?
Suppose NN is a power of 2, i.e. N=2pN=2p for some integer pp.
We can then
repeatedly decimate the sequence until each subsequence contains
only two points.
It is easily seen from Equation 4
that the 2 point DFT is a simple sum and difference of values.
X
(
0
)
=
x
(
0
)
+
x
(
1
)
X
(
1
)
=
x
(
0
)

x
(
1
)
X
(
0
)
=
x
(
0
)
+
x
(
1
)
X
(
1
)
=
x
(
0
)

x
(
1
)
(14)Figure 3 shows the flow diagram that results
for an 8 point DFT when we decimate 3 times.
Note that there are 3 stages of twiddle factors (in the first stage, the
twiddle factors simplify to “1”).
This is the flow diagram for the complete decimationintime
8 point FFT algorithm.
How many multiplies are required to compute it?
Write three Matlab functions to compute the
2, 4, and 8 point FFT's using the syntax
X = FFT2(x)
X = FFT4(x)
X = FFT8(x)
The function FFT2
should directly compute
the 2point DFT using Equation 14,
but the functions FFT4
and FFT8
should compute their respective FFT's using the divide and
conquer strategy.
This means that FFT8
should call FFT4
,
and FFT4
should call FFT2
.
Test your function FFT8
by using it to compute the DFT's
of the following signals. Compare these results to the previous ones.
 x(n)=δ(n)x(n)=δ(n) for N=8N=8
 x(n)=1forN=8x(n)=1forN=8
 x(n)=ej2πn/8x(n)=ej2πn/8 for N=8N=8

Submit the code for your functions
FFT2
, FFT4
and FFT8
.
 List the output of
FFT8
for the case x(n)=1x(n)=1 for N=8N=8.
 Calculate the total number of multiplies by twiddle factors
required for your 8 point FFT.
(A multiply is a multiplication by a real or complex
number.)
 Determine a formula for the number of multiplies required for
an N=2pN=2p point FFT. Leave the expression in terms of NN and pp.
How does this compare to the number of multiplies required
for direct implementation when p=10?
If you wrote the FFT4 and FFT8 functions properly, they should have almost the exact same form.
The only difference between them is the length of the input signal,
and the function called to compute the (N/2)pt DFTs.
Obviously, it's redundant to write a separate function for each specific
length DFT when they each have the same form.
The preferred method is to write a recursive function, which means that the function calls itself within the body.
It is imperative that a recursive function has a condition for exiting without
calling itself, otherwise it would never terminate.
Write a recursive function X=fft_stage(x)
that performs one
stage of the FFT algorithm for a powerof2 length signal.
An outline of the function is as follows:
 Determine the length of the input signal.
 If N=2, then the function should just compute the 2pt DFT as in
Equation 14, and then return.
 If N>2, then the function should perform the FFT steps described
previously (i.e. decimate, compute (N/2)pt DFTs, recombine),
calling fft_stage to compute the (N/2)pt DFTs.
Note that the body of this function should look very similar to previous
functions written in this lab.
Test fft_stage on the three 8pt signals given above,
and verify that it returns the same results as FFT8.
Submit the code for your fft_stage function.