EE2703 - Python 9 assignmentAssignment 8 - The Digital Fourier Transform - Solved

You are all studying DSP this semester, and you know the following from this course and the signals course:

•    A time function of the type f(t)u0(t) that grows slower than eαt for some α has a Laplace transform, F(s). The Laplace transform has an inverse transform which involves integrating vertically in the complex plane along a path that is to the right of all poles.

•    A finite energy function of the type f(t) has a Fourier Transform (and its inverse):

F tdt

f 

•    If f(t) is periodic with period 2π, the Fourier Transform collapses to the Fourier Series

f  cnejnt

n=−∞

where

cn   f(t)e−jntdt

for any t0. We can take t0 = 0 for convenience.

•    We can invert this picture and say, suppose f[n] are the samples of some function f(t), then we define the Z transform as

F n

n=−∞

Replacing z with ejθ we get

∞

F(ejθ) = ∑ f[n]e−jnθ n=−∞

So clearly F(z) is like the periodic time function that gives rise to the fourier series whose coefficients are the samples f[n].

F(ejθ) is called the Digital Spectrum of the samples f[n]. It is also called the DTFT of f[n]
•    F(ejθ) is continuous and periodic. f[n] is discrete and aperiodic. Suppose now f[n] is itself periodic with a period N, i.e.,

f[n+N] = f[n] ∀n

Then, it should have samples for its DTFT. This is true, and leads to the Discrete Fourier Transform or the DFT:

Suppose f[n] is a periodic sequence of samples, with a period N. Then the DTFT of the sequence is also a periodic sequence F[k] with the same period N. (You will prove this in the DSP course). So we have

                                                                           −1                                         nk            N−1                  nk

Ff[n]W

                                                                        n=0                                           N            n=0

N−1

1      −nk f[n]   =  ∑ F[k]W

N k=0

Here W   is used simply to make the equations less cluttered.

The values F[k] are what remains of the Digital Spectrum F(ejθ). We can consider them as the values of F(ejθ) for θ = 2πkN, since the first N terms in the expression for F(ej2πk/N) yield

                                                 F(ej2πk/N) =    +...

which is the same as the DFT expression. What about the remaining terms? They are repetitions of the above sum and help build up the delta function that is needed to take us from a continuous transform to discrete impulses.

•    What this means is that the DFT is a sampled version of the DTFT, which is the digital version of the analog Fourier Transform

In this assignment, we want to explore how to obtain the DFT, and how to recover the analog Fourier Tranform for some known functions by the proper sampling of the function.

The DFT in Python
There are two commands in Python, one to compute the forward fourier transform and the other to compute the inverse transform. They are

numpy.fft.fft() numpy.fft.ifft()

If you import pylab, both functions are imported into the local namespace. Let us try this on a random function:

from pylab import * x=rand(100) X=fft(x) y=ifft(X) c_[x,y] print abs(x-y).max()

•    When you run this code, the print statement will give you a value of about 10−15. The starting vector, x, and the vector got by running fft and ifft on x are essentially the same vector.

•    There is another thing to note. The command c_[x,y] tells Python to make a 100×2 matrix by using x and y as the columns. This is a quick and dirty way to print both out side by side. And what this shows is that x is pure real while y is very slightly complex. This is due to the finite accuracy of the CPU so that the ifft could not exactly undo what fft did.

•    Note that I used a 100 point vector, but it is much better to use a number that is 2k.

Let us now take a function we know about: y=sin(x). Let us use the fft and see what spectrum we get. We know that

ejx −e−jx

y = sin(x) =  

2j

So the expected spectrum is

1

                                                             Y(ω) =       [δ (ω−1)−δ (ω +1)]

2j

What do we get? We will use 128 points between 0 and 2π.

from pylab import * x=linspace(0,2*pi,128) y=sin(5*x) Y=fft(y) figure() subplot(2,1,1) plot(abs(Y),lw=2) grid(True) subplot(2,1,2) plot(unwrap(angle(Y)),lw=2) grid(True) show()

3               heg1 3i≡ from pylab import * x=linspace(0,2*pi,128) y=sin(5*x) Y=fft(y) figure() subplot(2,1,1) plot(abs(Y),lw=2) ylabel(r"$|Y|$",size=16) title(r"Spectrum of $\sin(5t)$") grid(True) subplot(2,1,2) plot(unwrap(angle(Y)),lw=2) ylabel(r"Phase of $Y$",size=16) xlabel(r"$k$",size=16) grid(True) savefig("fig9-1.png") show()

 

•    We get spikes as expected. But not where we expected. There is energy at nearby frequencies as well.

•    The spikes have a height of 64, not 0.5. We should divide by N to use it as a spectrum.

•    The phase at the spikes have a phase difference of π, which means they are opposite signs, which is correct

•    The actual phase at the spikes is near but not exactly correct.

•    We haven’t yet got the frequency axis in place

What could be going wrong? The problem is that the DFT treats the position axis as another frequency axis. So it expects the vector to be on the unit circle starting at 1. Our position vector started at 0 and went to 2π, which is correct. The fft gave an answer in the same value. So we need to shift the π to 2π portion to the left as it represents negative frequency. This can be done with a command called fftshift.

The second thing that is wrong is that both 0 and 2π are the same point. So our 128 points need to stop at the poing just before 2π. The best way to do this is to create a vector of 129 values and drop the last one: x=linspace(0,2*pi,129);x=x[:-1]

The ω axis has a highest frequency corresponding to N. This is because we have N points between 0 and 2π. So ∆x = 2π/N, and the sampling frequency becomes N/2π. Thus ωmax = N. The actual frequencies go upto half that frequency only.

4               heg2 4i≡ from pylab import * x=linspace(0,2*pi,129);x=x[:-1] y=sin(5*x) Y=fftshift(fft(y))/128.0 w=linspace(-64,63,128) figure() subplot(2,1,1) plot(w,abs(Y),lw=2) xlim([-10,10]) ylabel(r"$|Y|$",size=16) title(r"Spectrum of $\sin(5t)$") grid(True) subplot(2,1,2) plot(w,angle(Y),’ro’,lw=2) ii=where(abs(Y)>1e-3) plot(w[ii],angle(Y[ii]),’go’,lw=2) xlim([-10,10]) ylabel(r"Phase of $Y$",size=16) xlabel(r"$k$",size=16) grid(True) savefig("fig9-2.png") show()


The plot is now as follows:

 

•    Things have improved! The peaks are exactly at 0.5, and the remaining values are zero to machine precision.

•    Better, their position along the x axis is correct as our function is sin(5x) and we expect a spike at ω = ±5

•    Only the phase at the peak is meaningful and it is π/2 and −π/2 for the two peaks. The peak for ω = 5 should be 0.5/j, or 0.5exp(−jπ/2).

Suppose we now look at AM modulation. The function we want to analyse is

f(t) = (1+0.1cos(t))cos(10t)

Again, we use 128 points and generate the spectrum using the above code, only changing the definition of f(t).

6 heg3 6i≡ from pylab import * t=linspace(0,2*pi,129);t=t[:-1] y=(1+0.1*cos(t))*cos(10*t) Y=fftshift(fft(y))/128.0 w=linspace(-64,63,128) figure() subplot(2,1,1) plot(w,abs(Y),lw=2) xlim([-15,15]) ylabel(r"$|Y|$",size=16)

title(r"Spectrum of $\left(1+0.1\cos\left(t\right)\right)\cos\left(10t\right)$") grid(True)

subplot(2,1,2) plot(w,angle(Y),’ro’,lw=2) xlim([-15,15]) ylabel(r"Phase of $Y$",size=16) xlabel(r"$\omega$",size=16) grid(True) savefig("fig9-3.png") show()


And we get ...

 

Not what we expected!

•    There are three non-zero amplitudes on both sides, as expected.

•    All have zero phase, since we have cosines.

•    But not three spikes. Just a broader single spike

What went wrong? Well, we did not allow for enough frequencies. Let us use 512 samples instead. But how?

•    If we put more points between 0 and 2π,

t=linspace(0,2*pi,513);t=t[:-1]

it will give the same spacing but a high sampling frequency.

•    What we need to do is to stretch the t axis.

t=linspace(-4*pi,4*pi,513);t=t[:-1]

This will give us tighter spacing between frequency samples, but the same time spacing - i.e., the same sampling frequency.

8 heg4 8i≡ from pylab import * t=linspace(-4*pi,4*pi,513);t=t[:-1] y=(1+0.1*cos(t))*cos(10*t) Y=fftshift(fft(y))/512.0 w=linspace(-64,64,513);w=w[:-1]

figure() subplot(2,1,1) plot(w,abs(Y),lw=2) xlim([-15,15]) ylabel(r"$|Y|$",size=16) title(r"Spectrum of $\left(1+0.1\cos\left(t\right)\right)\cos\left(10t\right)$") grid(True) subplot(2,1,2) plot(w,angle(Y),’ro’,lw=2) xlim([-15,15]) ylabel(r"Phase of $Y$",size=16) xlabel(r"$\omega$",size=16) grid(True) savefig("fig9-4.png") show()


Finally we are there:

 

Perfect! We have three clear spikes. The phases at those spikes are zero to machine precision. The location of the spikes are 9, 10 and 11 radians per sec. The height of the side bands come from

0. 

All the amplitudes are real and positive as seen in the phase plot.

Assignment
1.    Work through the examples above.

2.    Generate the spectrum of sin3t and cos3t. Compare with what is expected.

3.    Generate the spectrum of cos(20t +5cos(t)). Plot phase points only where the magnitude is significant (say greater than 10−3). What do you think is happening?!

4.    The Gaussian exp   is not “bandlimited” in frequency. We want to get its spectrum accurate to 6 digits. Try different time ranges, and see what gets you a frequency domain that is so accurate. (The Fourier Transform of a Gaussian is a

Gaussian in ω. Look up the transform pair to confirm)