AMAT 415 — Matlab activity 3: Fourier series via FFT

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

Let x=(x_0,\ldots,x_{N-1}) be an n-element vector. Its DFT is the vector y=(y_0,\ldots,y_{N-1}) defined by

\displaystyle{y_n=\sum_{k=0}^{N-1}x_ke(-nk/N)},

where e(x)=e^{2\pi i x}. (In class yesterday, we had a factor of 1/N in front of the above sum. I’m dropping this because MATLAB uses different conventions.) Suppose x_k=f(2\pi k/N) for some function f defined on [0,2\pi]. Then

\displaystyle{y_n=\frac{N}{2\pi}\sum_{k=0}^{N-1}f(2\pi k/N)e^{-in(2\pi k/N)}\frac{2\pi}{N}}.

The sum in this expression is a “left endpoint” Riemann sum approximation to \displaystyle{\int_0^{2\pi}f(t)e^{-int}dt} using a subdivision of [0,2\pi] into N equal subintervals of width 2\pi/N. Thus, y_n\approx N\hat{f}(n). This means that we can approximate the Fourier coefficients of a function f(t) by computing the DFT of the vector x with x_k=f(2\pi k/N).

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 2^n. 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 f is \displaystyle{\sum_{n=-\infty}^\infty \hat{f}(n)e^{int}}. We’ve been ignoring the terms with n<0. No wonder our pictures are wrong! Even though y_{n} for n<0 doesn’t make sense, we recall that if f(t) is real-valued, then \hat{f}(-n)=\overline{\hat{f}(n)}.  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?

This entry was posted in Uncategorized. Bookmark the permalink.

Leave a Reply

Your email address will not be published. Required fields are marked *