Purdue University: ECE438 - Digital Signal Processing with Applications 1

ECE438 - Laboratory 6:
Discrete Fourier Transform and Fast Fourier
Transform Algorithms (Week 1)
October 6, 2010

1 Introduction
This is the first week of a two week laboratory that covers the Discrete Fourier Transform
(DFT) and Fast Fourier Transform (FFT) methods. The first week will introduce the DFT
and associated sampling and windowing effects, while the second week will continue the
discussion of the DFT and introduce the FFT.
In previous laboratories, we have used the Discrete-Time Fourier Transform (DTFT)
extensively for analyzing signals and linear time-invariant systems.
∞
jω
X
(DTFT) X(e ) = x(n)e−jωn (1)
n=−∞

1 Zπ
(inverse DTFT) x(n) = X(ejω )ejωn dω. (2)
2π −π
While the DTFT is very useful analytically, it usually cannot be exactly evaluated on a com-
puter because equation (1) requires an infinite sum and equation (2) requires the evaluation
of an integral.
The discrete Fourier transform (DFT) is a sampled version of the DTFT, hence it is
better suited for numerical evaluation on computers.
N
X −1
(DFT) XN (k) = x(n)e−j2πkn/N (3)
n=0

1 NX
−1
(inverse DFT) x(n) = XN (k)ej2πkn/N (4)
N k=0

Here XN (k) is an N point DFT of x(n). Note that XN (k) is a function of a discrete integer
k, where k ranges from 0 to N − 1.

Questions or comments concerning this laboratory should be directed to Prof. Charles A. Bouman,
School of Electrical and Computer Engineering, Purdue University, West Lafayette IN 47907; (765) 494-
0340; bouman@ecn.purdue.edu
Purdue University: ECE438 - Digital Signal Processing with Applications 2

In the following sections, we will study the derivation of the DFT from the DTFT, and
several DFT implementations. The fastest and most important implementation is known as
the fast Fourier transform (FFT). The FFT algorithm is one of the cornerstones of signal
processing.

2 Deriving the DFT from the DTFT

2.1 Truncating the Time-domain Signal
The DTFT usually cannot be computed exactly because the sum in equation (1) is infinite.
However, the DTFT may be approximately computed by truncating the sum to a finite
window. Let w(n) be a rectangular window of length N :
(
1 : 0≤n≤N −1
w(n) = . (5)
0 : else
Then we may define a truncated signal to be
xtr (n) = w(n)x(n) . (6)
The DTFT of xtr (n) is given by:
∞
Xtr (ejω ) =
X
xtr (n)e−jωn (7)
n=−∞

N
X −1
= x(n)e−jωn . (8)
n=0

We would like to compute X(ejω ), but as with filter design, the truncation window
distorts the desired frequency characteristics; X(ejω ) and Xtr (ejω ) are generally not equal.
To understand the relation between these two DTFT’s, we need to convolve in the frequency
domain (as we did in designing filters with the truncation technique):
1 Zπ
Xtr (ejω ) = X(ejσ )W (ej(ω−σ) )dσ (9)
2π −π
where W (ejω ) is the DTFT of w(n). Equation (9) is the periodic convolution of X(ejω ) and
W (ejω ). Hence the true DTFT, X(ejω ), is smoothed via convolution with W (ejω ) to produce
the truncated DTFT, Xtr (ejω ).
We can calculate W (ejω ):
∞
W (ejω ) =
X
w(n)e−jωn
n=−∞

N
X −1
= e−jωn
n=0
(
1−e−jωN
, for ω 6= 0, ±2π, . . .
= 1−e−jω (10)
N, for ω = 0, ±2π, . . .
Purdue University: ECE438 - Digital Signal Processing with Applications 3

For ω 6= 0, ±2π, . . ., we have:

e−jωN/2 ejωN/2 − e−jωN/2

W (ejω ) =
e−jω/2 ejω/2 − e−jω/2
sin(ωN/2)
= e−jω(N −1)/2 . (11)
sin(ω/2)
Notice that the magnitude of this function is similar to sinc(ωN/2) except that it is periodic
in ω with period 2π.

2.2 Frequency Sampling

Equation (8) contains a summation over a finite number of terms. However, we can only
evaluate (8) for a finite set of frequencies, ω. We must sample in the frequency domain to
compute the DTFT on a computer. We can pick any set of frequency points at which to
evaluate (8), but it is particularly useful to uniformly sample ω at N points, in the range
[0, 2π). If we substitute
ω = 2πk/N (12)
for k = 0, 1, . . . (N − 1) in (8), we find that
N

 −1 
jω
X
Xtr (e ) = x(n)e−jωn 
ω= 2πk

N

n=0 ω= 2πk
N

N
X −1
= x(n)e−j2πkn/N
n=0

= XN (k) .

In short, the DFT values result from sampling the DTFT of the truncated signal.

XN (k) = Xtr (ej2πk/N ) (13)

2.3 Windowing Effects

Down load DTFT.m

We willnext investigate the effect of windowing when computing the DFT of the signal
x(n) = cos π4 n truncated with a window of size N = 20.

1. In the same figure, plot the phase and magnitude of W (ejω ), using equations (10) and
(11).

2. Determine an expression for X(ejω ) (the DTFT of the non-truncated signal).

Purdue University: ECE438 - Digital Signal Processing with Applications 4

3. Truncate the signal x(n) using a window of size N = 20 and then use DTFT.m to
compute Xtr (ejω ). Make sure that the plot contains a least 512 points. Hint: Use the
command [X,w] = DTFT(x,512).

4. Plot the magnitude of Xtr (ejω ).

1 12

0.9

10
0.8

0.7
8

0.6

0.5 6

0.4

4
0.3

0.2
2

0.1

0 0
0 5 10 15 20 −4 −2 0 2 4
n ω

Figure 1: The plot of a Hamming window (left) and its DTFT (right).
Purdue University: ECE438 - Digital Signal Processing with Applications 5

INLAB REPORT:

1. Submit the plot of the phase and magnitude of W (ejω ).

2. Submit the analytical expression for X(ejω ).

3. Submit the magnitude plot of Xtr (ejω ).

4. Describe the difference between |Xtr (ejω )| and |X(ejω )|. What is the reason for this
difference?

5. Comment on the effects of using a different window for w(n). For example, what would
you expect your plots to look like if you had used a Hamming window in place of the
truncation (rectangular) window? (See Fig. 1 for a plot of a Hamming window of
length 20 and its DTFT.)

3 The Discrete Fourier Transform

3.1 Computing the DFT
We will now develop our own DFT functions to help our understanding of how the DFT
comes from the DTFT. Write your own Matlab function to implement the DFT of equation
(3). Use the syntax
X = DFTsum(x)
where x is an N point vector containing the values x(0), · · · , x(N − 1) and X is the corre-
sponding DFT. Your routine should implement the DFT exactly as specified by (3) using
for-loops for n and k, and computing the exponentials as they appear.√ Note: In Matlab, j
may be computed with the command j=sqrt(-1). If you use j = −1, remember not to
use j as an index in your for-loop.
Test your routine DFTsum by computing XN (k) for each of the following cases:

1. x(n) = δ(n) for N = 10.

2. x(n) = 1 for N = 10.

3. x(n) = ej2πn/10 for N = 10.

4. x(n) = cos(2πn/10) for N = 10.

Plot the magnitude of each of the DFT’s. In addition, derive simple closed-form analytical
expressions for the DFT (not the DTFT!) of each signal.
Purdue University: ECE438 - Digital Signal Processing with Applications 6

INLAB REPORT:

1. Submit a listing of your code for DFTsum.

2. Submit the magnitude plots.

3. Submit the corresponding analytical expressions.

Write a second Matlab function for computing the inverse DFT of (4). Use the syntax
x = IDFTsum(X)
where X is the N point vector containing the DFT and x is the corresponding time-domain
signal. Use IDFTsum to invert each of the DFT’s computed in the previous problem. Plot
the magnitudes of the inverted DFT’s, and verify that those time-domain signals match
the original ones. Use abs(x) to eliminate any imaginary parts which roundoff error may
produce.

INLAB REPORT:

1. Submit the listing of your code for IDFTsum.

2. Submit the four time-domain IDFT plots.

3.2 Matrix Representation of the DFT

The DFT of (3) can be implemented as a matrix-vector product. To see this, consider the
equation
X = Ax (14)
where A is an N × N matrix, and both X and x are N × 1 column vectors. This matrix
product is equivalent to the summation
N
X
Xk = Akn xn . (15)
n=1

where Akn is the matrix element in the k th row and nth column of A. By comparing equations
(3) and (15) we see that for the DFT,

Akn = e−j2π(k−1)(n−1)/N (16)

The −1’s are in the exponent because Matlab indices start at 1, not 0. For this section,
we need to:

1. Write a Matlab function A = DFTmatrix(N) that returns the N × N DFT matrix A.

Purdue University: ECE438 - Digital Signal Processing with Applications 7

2. Use the matrix A to compute the DFT of the following signals1 . Confirm that the
results are the same as in section 3.1.

(a) x(n) = δ(n) for N = 10

(b) x(n) = 1 for N = 10
(c) x(n) = ej2πn/N for N = 10

INLAB REPORT:

1. Print out the matrix A for N = 5.

2. Hand in the three magnitude plots of the DFT’s.

3. How many multiplies are required to compute an N point DFT using the matrix
method? (Consider a multiply as the multiplication of either complex or real numbers.)

As with the DFT, the inverse DFT may also be represented as a matrix-vector product.

x = BX

For this section,

1. Write an analytical expression for the elements of the inverse DFT matrix B, using
the form of equation (16).

2. Write a Matlab function B = IDFTmatrix(N) that returns the N × N inverse DFT

matrix B.

3. Compute the matrices A and B for N = 5. Then compute the matrix product C =
BA.

INLAB REPORT:

1. Hand in your analytical expression for the elements of B.

2. Print out the matrix B for N = 5.

3. Print out the elements of C = BA. What form does C have? Why does it have this
form?

1
Remember that the symbol ∗ is used for matrix multiplication in Matlab, and that .′ performs a simple
transpose on a vector or matrix. An apostrophe without the period is a conjugate transpose.
Purdue University: ECE438 - Digital Signal Processing with Applications 8

3.3 Computation Time Comparison

Help on cputime function

Although the operations performed by DFTsum are mathematically identical to a matrix

product, the computation times for these two DFT’s in Matlab are quite different. (This is
despite the fact that the computational complexity of two procedures is of the same order!)
This exercise will underscore why you should try to avoid using for loops in Matlab, and
wherever possible, try to formulate your computations using matrix/vector products.
To see this, do the following:

1. Compute the signal x(n) = cos(2πn/10) for N = 512.

2. Compute the matrix A for N = 512.

3. Compare the computation time of X = DFTsum(x) with a matrix implementation X = A*x

by using the cputime function before and after the program execution. Do not include
the computation of A in your timing calculations.

INLAB REPORT:
Report the cputime required for each of the two implementations. Which method is faster?
Which method requires less storage?