This activity is meant to be instructive. I won’t ask you to submit your code or plots.

Let be an -element vector. Its DFT is the vector defined by

,

where . (In class yesterday, we had a factor of in front of the above sum. I’m dropping this because MATLAB uses different conventions.) Suppose for some function defined on . Then

.

The sum in this expression is a “left endpoint” Riemann sum approximation to using a subdivision of into equal subintervals of width . Thus, . This means that we can approximate the Fourier coefficients of a function by computing the DFT of the vector with .

The Fast Fourier Transform is an algorithm for computing the DFT very efficiently. In MATLAB, invoking** fft(x)** gives you the DFT of the vector **x**. (Note that MATLAB indexes vectors starting from 1 rather than from 0. From the perspective of the mathematics of the DFT, it’s better to start indexing from 0 as above. Keep this shift in mind.) Let’s repeat some of the calculations of last week’s lab using the **fft** command. The FFT algorithm works best on vector of length . So let **N=2^7**, say, and set **t=0:2*pi/N:(2*pi-2*pi/N);**. We’ll start with the step function **x=t<pi;**. Let **y** be its DFT, i.e., set **y=fft(x);**. Now let’s try to reconstruct x from y. Initialize the vector and, for various small values of **m**, perform:

**s=zeros(1,length(x));**

** for k=1:m**

** s = s + y(k)/N*exp(i*(k-1)*t);**

** end**

** plot(t,real(s))**

Do the graphs look right? Plot **x** and **s** on the same graph: **plot(t,x,t,s)**. Something’s off. What’s going on?

Let’s try a different function. Let **x=[t(1:N/2) t(N/2)-t(1:N/2)];**. (Plot it and figure out why the command you just pasted in gives you the graph you see.) Now run through the above process with this **x**. Still doesn’t look right…

Wait a second. The Fourier series of a function is . We’ve been ignoring the terms with . No wonder our pictures are wrong! Even though for doesn’t make sense, we recall that if is real-valued, then . Let’s take this into account:

**s=zeros(1,length(x))+y(1)/N;**

** for k=2:m**

** s = s + y(k)/N*exp(i*(k-1)*t) + conj(y(k))/N*exp(i*(1-k)*t);**

** end**

** plot(t,real(s),t,x)**

That’s more like it! Try a few more functions **x** like the ones you worked with last week.

Finally, I’d like to acquaint you with a nice command. Check out **stem(y)**, **stem(abs(y))** and compare them with **plot(real(y))** and **plot(abs(y))**. Since DFTs are just **N** element vectors, a “plot” of them should just be **N** points *not* connected by lines. This is what stem gives you. It’s really nice for visualizing discrete signals in the frequency domain, especially when **N** isn’t too big so the circles it plots don’t start overlapping. What interesting symmetries do you observe in these stem plots?