Skip to content Skip to navigation

Connexions

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

Navigation

Lenses

What is a lens?

Definition of a lens

Lenses

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

What is in a lens?

Lens makers point to Connexions 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 Connexions 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.
  • 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.

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 6b - Discrete Fourier Transform and FFT (part 2)

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

User rating (How does the rating system work?)
Ratings

Ratings allow you to judge the quality of modules. If other users have ranked the module then its average rating is displayed below. Ratings are calculated on a scale from one star (Poor) to five stars (Excellent).

How to rate a module

Hover over the star that corresponds to the rating you wish to assign. Click on the star to add your rating. Your rating should be based on the quality of the content. You must have an account and be logged in to rate content.

:
(0 ratings)

Note: Your browser may not currently support MathML. See our browser support page for additional details. You can always view the correct math in the PDF version.

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 second week of a two week laboratory that covers the Discrete Fourier Transform (DFT) and Fast Fourier Transform (FFT). The first week introduced the DFT and associated sampling and windowing effects. This laboratory will continue the discussion of the DFT and will introduce the FFT.

Continuation of DFT Analysis

This section continues the analysis of the DFT started in the previous week's laboratory.

(DFT) X N ( k ) = ∑ n = 0 N - 1 x ( n ) 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 (DFT) X N ( k ) = ∑ n = 0 N - 1 x ( n ) 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 (1)

Shifting the Frequency Range

In this section, we will illustrate a representation for the DFT of Equation 1 that is a bit more intuitive. First create a Hamming window xx of length N=20N=20, using the Matlab command x = hamming(20) . Then use your Matlab function DFTsum to compute the 20 point DFT of xx. Plot the magnitude of the DFT, |X20(k)||X20(k)|, versus the index kk. Remember that the DFT index kk starts at 0 not 1!

INLAB REPORT:

Hand in the plot of the |X20(k)||X20(k)|. Circle the regions of the plot corresponding to low frequency components.

Our plot of the DFT has two disadvantages. First, the DFT values are plotted against kk rather then the frequency ωω. Second, the arrangement of frequency samples in the DFT goes from 0 to 2π2π rather than from -π-π to ππ, as is conventional with the DTFT. In order to plot the DFT values similar to a conventional DTFT plot, we must compute the vector of frequencies in radians per sample, and then “rotate” the plot to produce the more familiar range, -π-π to ππ.

Let's first consider the vector w of frequencies in radians per sample. Each element of w should be the frequency of the corresponding DFT sample X(k)X(k), which can be computed by

ω = 2 π k / N k ∈ [ 0 , ⋯ , N - 1 ] . ω = 2 π k / N k ∈ [ 0 , ⋯ , N - 1 ] . (2)

However, the frequencies should also lie in the range from -π-π to ππ. Therefore, if ω≥πω≥π, then it should be set to ω-2πω-2π. An easy way of making this change in Matlab 5.1 is w(w>=pi) = w(w>=pi)-2*pi .

The resulting vectors X and w are correct, but out of order. To reorder them, we must swap the first and second halves of the vectors. Fortunately, Matlab provides a function specifically for this purpose, called fftshift.

Write a Matlab function to compute samples of the DTFT and their corresponding frequencies in the range -π-π to ππ. Use the syntax

[X,w] = DTFTsamples(x)

where x is an NN point vector, X is the length NN vector of DTFT samples, and w is the length NN vector of corresponding radial frequencies. Your function DTFTsamples should call your function DFTsum and use the Matlab function fftshift.

Use your function DTFTsamples to compute DTFT samples of the Hamming window of length N=20N=20. Plot the magnitude of these DTFT samples versus frequency in rad/sample.

INLAB REPORT:

Hand in the code for your function DTFTsamples. Also hand in the plot of the magnitude of the DTFT samples.

Zero Padding

The spacing between samples of the DTFT is determined by the number of points in the DFT. This can lead to surprising results when the number of samples is too small. In order to illustrate this effect, consider the finite-duration signal

x ( n ) = sin ( 0 . 1 π n ) 0 ≤ n ≤ 49 0 otherwise x ( n ) = sin ( 0 . 1 π n ) 0 ≤ n ≤ 49 0 otherwise (3)

In the following, you will compute the DTFT samples of x(n)x(n) using both N=50N=50 and N=200N=200 point DFT's. Notice that when N=200N=200, most of the samples of x(n)x(n) will be zeros because x(n)=0x(n)=0 for n≥50n≥50. This technique is known as “zero padding”, and may be used to produce a finer sampling of the DTFT.

For N=50N=50 and N=200N=200, do the following:

  1. Compute the vector x containing the values x(0),⋯,x(N-1)x(0),⋯,x(N-1).
  2. Compute the samples of X(k)X(k) using your function DTFTsamples.
  3. Plot the magnitude of the DTFT samples versus frequency in rad/sample.

INLAB REPORT:

Submit your two plots of the DTFT samples for N=50N=50 and N=200N=200. Which plot looks more like the true DTFT? Explain why the plots look so different.

The Fast Fourier Transform Algorithm

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 divide-and-conquer, splits the problem into two smaller problems. The second concept, known as recursion, applies this divide-and-conquer 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/2-1N/2-1, n=2mn=2m will take on all even integer values between 0 and N-2N-2. Similarly, in the second sum, we replace nn by 2m+12m+1. Then as we sum mm from 0 to N/2-1N/2-1, n=2m+1n=2m+1 will take on all odd integer values between 0 and N-1N-1. 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 even-numbered data points in the original sequence. Similarly, the second sum in Equation 8 is an N/2N/2 point DFT of the odd-numbered data points in the original sequence. To obtain the NN point DFT of the complete sequence, we multiply the DFT of the odd-numbered data points by the complex exponential factor e-j2πk/Ne-j2π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 odd-numbered 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/2-1n=0,...,N/2-1. 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 e-j2πk/Ne-j2π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=e-j2πk/NWNk=e-j2πk/N are commonly known as the twiddle factors.

Figure 1: Divide and conquer DFT of equation (13). The NN-point DFT is computed using the two N/2N/2-point DFT's X0(N/2)(k)X0(N/2)(k) and X1(N/2)(k)X1(N/2)(k) .
Figure 1 (dit1.png)

Figure 1 shows a graphical interpretation of Equation 13 which we will refer to as the “divide-and-conquer 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.

Implementation of Divide-and-Conquer DFT

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:

  1. Separate the samples of x into even and odd points.

    Hint:

    The Matlab command x0 = x(1:2:N) can be used to obtain the “even” points.
  2. Use your function DFTsum to compute the two N/2N/2 point DFT's.
  3. Multiply by the twiddle factors WNk=e-j2πk/NWNk=e-j2πk/N.
  4. Combine the two DFT's to form XX.

Test your function dcDFT by using it to compute the DFT's of the following signals.

  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. Submit the code for your function dcDFT.
  2. 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.)
    HINT:
    Refer to the diagram of Figure 1, and remember to consider the N/2N/2 point DFTs.

Recursive Divide and Conquer

Figure 2: Recursion of the decimation-in-time process. Now each N/2N/2-point is calculated by combining two N/4N/4-point DFT's.
Figure 2 (dit2.png)

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 decimation-in-time 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: 8-Point FFT.
Figure 3 (8pt_fft.png)

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 decimation-in-time 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 2-point 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.

  1. x(n)=δ(n)x(n)=δ(n) for N=8N=8
  2. x(n)=1forN=8x(n)=1forN=8
  3. x(n)=ej2πn/8x(n)=ej2πn/8 for N=8N=8

INLAB REPORT

  1. Submit the code for your functions FFT2, FFT4 and FFT8.
  2. List the output of FFT8 for the case x(n)=1x(n)=1 for N=8N=8.
  3. 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.)
  4. 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 power-of-2 length signal. An outline of the function is as follows:

  1. Determine the length of the input signal.
  2. If N=2, then the function should just compute the 2-pt DFT as in Equation 14, and then return.
  3. If N>2, then the function should perform the FFT steps described previously (i.e. decimate, compute (N/2)-pt DFTs, re-combine), 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 8-pt signals given above, and verify that it returns the same results as FFT8.

INLAB REPORT:

Submit the code for your fft_stage function.

References

  1. J. W. Cooley and J. W. Tukey. (1965, April). An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation, 19(90), 297-301.

Content actions

Give Feedback:

E-mail the module author | Rate module ( How does the rating system work?)

Rating system

Ratings

Ratings allow you to judge the quality of modules. If other users have ranked the module then its average rating is displayed below. Ratings are calculated on a scale from one star (Poor) to five stars (Excellent).

How to rate a module

Hover over the star that corresponds to the rating you wish to assign. Click on the star to add your rating. Your rating should be based on the quality of the content. You must have an account and be logged in to rate content.

(0 ratings)

Download:

Module PDF

Add module to:

My Favorites (?)

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

| A lens (?)

Definition of a lens

Lenses

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

What is in a lens?

Lens makers point to Connexions 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 Connexions 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