Skip to content Skip to navigation

Connexions

You are here: Home » Content » Lab 6a - Discrete Fourier Transform and FFT (part 1)

Navigation

Lenses

What is a lens?

Definition of a lens

Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

Who can create a lens?

Any individual member, a community, or a respected organization.

What are tags? tag icon

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

This content is ...

Affiliated with (What does "Affiliated with" mean?)

This content is either by members of the organizations listed or about topics related to the organizations listed. Click each link to see a list of all content affiliated with the organization.
  • NSF Partnership display tagshide tags

    This module is included inLens: NSF Partnership in Signal Processing
    By: Sidney BurrusAs a part of collection: "Purdue Digital Signal Processing Labs (ECE 438)"

    Click the "NSF Partnership" link to see all content affiliated with them.

    Click the tag icon tag icon to display tags associated with this content.

  • Featured Content display tagshide tags

    This module is included inLens: Connexions Featured Content
    By: ConnexionsAs a part of collection: "Purdue Digital Signal Processing Labs (ECE 438)"

    Click the "Featured Content" link to see all content affiliated with them.

    Click the tag icon tag icon to display tags associated with this content.

Also in these lenses

  • UniqU content

    This module is included inLens: UniqU's lens
    By: UniqU, LLCAs a part of collection: "Purdue Digital Signal Processing Labs (ECE 438)"

    Click the "UniqU content" link to see all content selected in this lens.

  • Lens for Engineering

    This module is included inLens: Lens for Engineering
    By: Sidney Burrus

    Click the "Lens for Engineering" link to see all content selected in this lens.

Recently Viewed

This feature requires Javascript to be enabled.

Tags

(What is a tag?)

These tags come from the endorsement, affiliation, and other lenses that include this content.
 

Lab 6a - Discrete Fourier Transform and FFT (part 1)

Module by: Charles A. Bouman. E-mail the author

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

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.

(DTFT) X ( e j ω ) = ∑ n = - ∞ ∞ x ( n ) e - j ω n (DTFT) X ( e j ω ) = ∑ n = - ∞ ∞ x ( n ) e - j ω n
(1)
(inverse DTFT) x ( n ) = 1 2 π ∫ - π π X ( e j ω ) e j ω n d ω . (inverse DTFT) x ( n ) = 1 2 π ∫ - π π X ( e j ω ) e j ω n d ω .
(2)

While the DTFT is very useful analytically, it usually cannot be exactly evaluated on a computer 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.

(DFT) X N ( k ) = ∑ n = 0 N - 1 x ( n ) e - j 2 π k n / N (DFT) X N ( k ) = ∑ n = 0 N - 1 x ( n ) e - j 2 π k n / N
(3)
(inverse DFT) x ( n ) = 1 N ∑ k = 0 N - 1 X N ( k ) e j 2 π k n / N (inverse DFT) x ( n ) = 1 N ∑ k = 0 N - 1 X N ( k ) e j 2 π k n / N
(4)

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

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.

Deriving the DFT from the DTFT

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)w(n) be a rectangular window of length NN:

w ( n ) = 1 0 ≤ n ≤ N - 1 0 else . w ( n ) = 1 0 ≤ n ≤ N - 1 0 else .
(5)

Then we may define a truncated signal to be

x tr ( n ) = w ( n ) x ( n ) . x tr ( n ) = w ( n ) x ( n ) .
(6)

The DTFT of x tr (n)x tr (n) is given by:

X tr ( e j ω ) = ∑ n = - ∞ ∞ x tr ( n ) e - j ω n = ∑ n = 0 N - 1 x ( n ) e - j ω n . X tr ( e j ω ) = ∑ n = - ∞ ∞ x tr ( n ) e - j ω n = ∑ n = 0 N - 1 x ( n ) e - j ω n .
(7)

We would like to compute X(ejω)X(ejω), but as with filter design, the truncation window distorts the desired frequency characteristics; X(ejω)X(ejω) and X tr (ejω)X tr (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):

X tr ( e j ω ) = 1 2 π ∫ - π π X ( e j σ ) W ( e j ( ω - σ ) ) d σ X tr ( e j ω ) = 1 2 π ∫ - π π X ( e j σ ) W ( e j ( ω - σ ) ) d σ
(8)

where W(ejω)W(ejω) is the DTFT of w(n)w(n). Equation 8 is the periodic convolution of X(ejω)X(ejω) and W(ejω)W(ejω). Hence the true DTFT, X(ejω)X(ejω), is smoothed via convolution with W(ejω)W(ejω) to produce the truncated DTFT, X tr (ejω)X tr (ejω).

We can calculate W(ejω)W(ejω):

W ( e j ω ) = ∑ n = - ∞ ∞ w ( n ) e - j ω n = ∑ n = 0 N - 1 e - j ω n = 1 - e - j ω N 1 - e - j ω , for ω ≠ 0 , ± 2 π , ... N , for ω = 0 , ± 2 π , ... W ( e j ω ) = ∑ n = - ∞ ∞ w ( n ) e - j ω n = ∑ n = 0 N - 1 e - j ω n = 1 - e - j ω N 1 - e - j ω , for ω ≠ 0 , ± 2 π , ... N , for ω = 0 , ± 2 π , ...
(9)

For ω≠0,±2π,...ω≠0,±2π,..., we have:

W ( e j ω ) = e - j ω N / 2 e - j ω / 2 e j ω N / 2 - e - j ω N / 2 e j ω / 2 - e - j ω / 2 = e - j ω ( N - 1 ) / 2 sin ( ω N / 2 ) sin ( ω / 2 ) . W ( e j ω ) = e - j ω N / 2 e - j ω / 2 e j ω N / 2 - e - j ω N / 2 e j ω / 2 - e - j ω / 2 = e - j ω ( N - 1 ) / 2 sin ( ω N / 2 ) sin ( ω / 2 ) .
(10)

Notice that the magnitude of this function is similar to sinc(ωN/2)sinc(ωN/2) except that it is periodic in ωω with period 2π2π.

Frequency Sampling

Equation 7 contains a summation over a finite number of terms. However, we can only evaluate Equation 7 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 Equation 7, but it is particularly useful to uniformly sample ωω at NN points, in the range [0,2π)[0,2π). If we substitute

ω = 2 π k / N ω = 2 π k / N
(11)

for k=0,1,...(N-1)k=0,1,...(N-1) in Equation 7, we find that

X tr ( e j ω ) ω = 2 π k N = ∑ n = 0 N - 1 x ( n ) e - j ω n ω = 2 π k N = ∑ n = 0 N - 1 x ( n ) e - j 2 π k n / N = X N ( k ) . X tr ( e j ω ) ω = 2 π k N = ∑ n = 0 N - 1 x ( n ) e - j ω n ω = 2 π k N = ∑ n = 0 N - 1 x ( n ) e - j 2 π k n / N = X N ( k ) .
(12)

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

X N ( k ) = X tr ( e j 2 π k / N ) X N ( k ) = X tr ( e j 2 π k / N )
(13)

Windowing Effects

Download DTFT.m for the following section.

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

  1. In the same figure, plot the phase and magnitude of W(ejω)W(ejω), using equations Equation 9 and Equation 10.
  2. Determine an expression for X(ejω)X(ejω) (the DTFT of the non-truncated signal).
  3. Truncate the signal x(n)x(n) using a window of size N=20N=20 and then use DTFT.m to compute X tr (ejω)X tr (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 X tr (ejω)X tr (ejω).
Figure 1: The plot of a Hamming window (left) and its DTFT (right).
Figure 1 (hamming.png)

INLAB REPORT

  1. Submit the plot of the phase and magnitude of W(ejω)W(ejω).
  2. Submit the analytical expression for X(ejω)X(ejω).
  3. Submit the magnitude plot of X tr (ejω)X tr (ejω).
  4. Describe the difference between |X tr (ejω)||X tr (ejω)| and |X(ejω)||X(ejω)|. What is the reason for this difference?
  5. Comment on the effects of using a different window for w(n)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 Figure 1 for a plot of a Hamming window of length 20 and its DTFT.)

The Discrete Fourier Transform

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 Equation 3. Use the syntax

X = DFTsum(x)

where x is an NN point vector containing the values x(0),⋯,x(N-1)x(0),⋯,x(N-1) and X is the corresponding DFT. Your routine should implement the DFT exactly as specified by Equation 3 using for-loops for nn and kk, and compute the exponentials as they appear. Note: In Matlab, "j" may be computed with the command j=sqrt(-1) . If you use j=-1j=-1, remember not to use j as an index in your for-loop.

Test your routine DFTsum by computing XN(k)XN(k) for each of the following cases:

  1. x(n)=δ(n)x(n)=δ(n) for N=10N=10.
  2. x(n)=1x(n)=1 for N=10N=10.
  3. x(n)=ej2πn/10x(n)=ej2πn/10 for N=10N=10.
  4. x(n)=cos(2πn/10)x(n)=cos(2πn/10) for N=10N=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.

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 Equation 4. Use the syntax

x = IDFTsum(X)

where X is the NN 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.

Matrix Representation of the DFT

The DFT of Equation 3 can be implemented as a matrix-vector product. To see this, consider the equation

X = A x X = A x
(14)

where AA is an N×NN×N matrix, and both XX and xx are N×1N×1 column vectors. This matrix product is equivalent to the summation

X k = ∑ n = 1 N A k n x n . X k = ∑ n = 1 N A k n x n .
(15)

where AknAkn is the matrix element in the kthkth row and nthnth column of A . By comparing Equation 3 and Equation 15 we see that for the DFT,

A k n = e - j 2 π ( k - 1 ) ( n - 1 ) / N A k n = e - j 2 π ( k - 1 ) ( n - 1 ) / N
(16)

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

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

    Note:

    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.
  • Use the matrix A to compute the DFT of the following signals. Confirm that the results are the same as in the previous section.
    1. x(n)=δ(n)x(n)=δ(n) for N=10N=10
    2. x(n)=1forN=10x(n)=1forN=10
    3. x(n)=ej2πn/Nx(n)=ej2πn/N for N=10N=10

INLAB REPORT

  1. Print out the matrix AA for N=5N=5.
  2. Hand in the three magnitude plots of the DFT's.
  3. How many multiplies are required to compute an NN 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 = B X x = B X
(17)

For this section,

  1. Write an analytical expression for the elements of the inverse DFT matrix BB, using the form of Equation 16.
  2. Write a Matlab function B = IDFTmatrix(N) that returns the NxN inverse DFT matrix B.
  3. Compute the matrices AA and BB for N=5N=5. Then compute the matrix product C=BAC=BA.

INLAB REPORT

  1. Hand in your analytical expression for the elements of BB.
  2. Print out the matrix BB for N=5N=5.
  3. Print out the elements of C=BAC=BA. What form does CC have? Why does it have this form?

Computation Time Comparison

Click here for help on the 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)x(n)=cos(2πn/10) for N=512N=512.
  2. Compute the matrix AA for N=512N=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?

Content actions

Download module as:

PDF | EPUB (?)

What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

Downloading to a reading device

For detailed instructions on how to download this content's EPUB to your specific device, click the "(?)" link.

| More downloads ...

Add module to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

Definition of a lens

Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

Who can create a lens?

Any individual member, a community, or a respected organization.

What are tags? tag icon

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks