EE2703-Assignment 9 Spectra of Non-Periodic Signals Solved

we want to look at non-periodic signals. Our first signal will be sin . Obtained over 0 to 2π with 64 samples, this function looks as follows:

1            h* 1i≡

heg1 2i

2            heg1 2i≡  (1) from pylab import * t=linspace(-pi,pi,65);t=t[:-1] dt=t[1]-t[0];fmax=1/dt y=sin(sqrt(2)*t) y[0]=0 # the sample corresponding to -tmax should be set zeroo y=fftshift(y) # make y start with y(t=0) Y=fftshift(fft(y))/64.0 w=linspace(-pi*fmax,pi*fmax,65);w=w[:-1] figure() subplot(2,1,1) plot(w,abs(Y),lw=2) xlim([-10,10]) ylabel(r"$|Y|$",size=16) title(r"Spectrum of $\sin\left(\sqrt{2}t\right)$") grid(True) subplot(2,1,2) plot(w,angle(Y),’ro’,lw=2) xlim([-10,10]) ylabel(r"Phase of $Y$",size=16) xlabel(r"$\omega$",size=16) grid(True) savefig("fig10-1.png") show()


 

We expected two spikes, but what we got were two peaks each with two values and a gradually decaying magnitude. The phase is correct though - but take a look at the code. There is one line there:

y[0]=0 # the sample corresponding to -tmax should be set zero

What is this for? An antisymmetric function has a purely imaginary fourier transform. But what about an antisymmetric set of samples? Suppose

                                                         y[0] =           0           sin(0)

N

                                                          y[i] =     −y[N −i]        i = 1,2...     −1

2

 N

y  tN

The DFT of this sequence will give us

                        Y[k]     =                             

                                            n=0                          N

N

N

                                     =              y[n]   exp     kn    −exp     −  kn       +y[   ]exp(πkj)

                                                 =1                                N                             N                     2

N

                                                                                                    k      N

                                    =                                                   +(−1) y[   ]

2

But this is no longer pure imaginary! Since we need that property, we set y[N/2] to zero. That gives us a purely imaginary phase. But what to do about the magnitude? And what went wrong?

To understand what went wrong, let us plot the time function over several time periods.

4            heg2 4i≡ from pylab import * t1=linspace(-pi,pi,65);t1=t1[:-1] t2=linspace(-3*pi,-pi,65);t2=t2[:-1] t3=linspace(pi,3*pi,65);t3=t3[:-1]

# y=sin(sqrt(2)*t) figure(2) plot(t1,sin(sqrt(2)*t1),’b’,lw=2) plot(t2,sin(sqrt(2)*t2),’r’,lw=2) plot(t3,sin(sqrt(2)*t3),’r’,lw=2) ylabel(r"$y$",size=16) xlabel(r"$t$",size=16) title(r"$\sin\left(\sqrt{2}t\right)$") grid(True) savefig("fig10-2.png") show()


 

The blue line connects the points whose DFT we took. The red lines show the continuation of the function. Quite clearly, even though sin  is a periodic function, the portion between −π and π is not the part that can be replicated to generate the function. So which function is the DFT trying to fourier analyse? For that we have to replicate just the blue points. And that is shown below:

5            heg3 5i≡ from pylab import * t1=linspace(-pi,pi,65);t1=t1[:-1] t2=linspace(-3*pi,-pi,65);t2=t2[:-1] t3=linspace(pi,3*pi,65);t3=t3[:-1] y=sin(sqrt(2)*t1) figure(3) plot(t1,y,’bo’,lw=2) plot(t2,y,’ro’,lw=2) plot(t3,y,’ro’,lw=2) ylabel(r"$y$",size=16) xlabel(r"$t$",size=16)

title(r"$\sin\left(\sqrt{2}t\right)$ with $t$ wrapping every $2\pi$ ") grid(True) savefig("fig10-3.png") show()

 

Gibbs Phenomenon
This is something you know very well. The Fourier transform of the box function

f 

is given by

2sin(ωt0)

F(ω) =  ω

The spectrum of the box function decays very slowly, as 2/ω.

Now our function is an odd function with a big jump. So let us consider the periodic ramp:

f(t) = t, −π < t < π

Then the fourier series of this ramp is

f 

Again the coefficients decay very slowly.

The DFT is just like the fourier series, except that both time and frequency are samples. So, if the time samples are like a ramp, the frequency samples will decay as 1/ω. Let us verify this for the ramp itself

7 heg4 7i≡ from pylab import * t=linspace(-pi,pi,65);t=t[:-1] dt=t[1]-t[0];fmax=1/dt y=t y[0]=0 # the sample corresponding to -tmax should be set zeroo y=fftshift(y) # make y start with y(t=0) Y=fftshift(fft(y))/64.0 w=linspace(-pi*fmax,pi*fmax,65);w=w[:-1] figure() semilogx(abs(w),20*log10(abs(Y)),lw=2) xlim([1,10]) ylim([-20,0]) xticks([1,2,5,10],["1","2","5","10"],size=16) ylabel(r"$|Y|$ (dB)",size=16) title(r"Spectrum of a digital ramp") xlabel(r"$\omega$",size=16) grid(True) savefig("fig10-4.png") show()


 

Clearly the spectrum decays as 20 dB per decade, which corresponds to 1ω. The big jumps at nπ force this slowly decaying spectrum, which is why we don’t see the expected spikes for the spectrum of sin .

Windowing
So what do we do? Well the spikes happen at the end of the periodic interval. So we damp the function near there, i.e., we multiply our function sequence f[n] by a “window” sequence w[n]: g(n) = f(n)w(n)

The new spectrum is got by convolving the two fourier transforms:

N−1
Gk = ∑ FnWk−n n=0

Suppose fn is a sinusoid. Then Fk has two spikes. But the two spikes are now smeared out by Wk. So we expect to get broader peaks. But what this also does is to suppress the jump at the edge of the window. The window we will use is called the Hamming window:

w 

Let us look at our time sequence for sin  now ...

8                 heg5 8i≡ from pylab import * t1=linspace(-pi,pi,65);t1=t1[:-1] t2=linspace(-3*pi,-pi,65);t2=t2[:-1] t3=linspace(pi,3*pi,65);t3=t3[:-1] n=arange(64) wnd=fftshift(0.54+0.46*cos(2*pi*n/63)) y=sin(sqrt(2)*t1)*wnd figure(3) plot(t1,y,’bo’,lw=2) plot(t2,y,’ro’,lw=2) plot(t3,y,’ro’,lw=2) ylabel(r"$y$",size=16) xlabel(r"$t$",size=16) title(r"$\sin\left(\sqrt{2}t\right)\times w(t)$ with $t$ wrapping every $2\pi$ ") grid(True) savefig("fig10-5.png") show()


 

The jump is still there, but it is much reduced. There is a little bit of magic in keeping some of the jump - it gives us an extra 10 db of suppression.

Now let us take the DFT of this sequence and see what we get:

10 heg6 10i≡ from pylab import * t=linspace(-pi,pi,65);t=t[:-1] dt=t[1]-t[0];fmax=1/dt n=arange(64) wnd=fftshift(0.54+0.46*cos(2*pi*n/63)) y=sin(sqrt(2)*t)*wnd y[0]=0 # the sample corresponding to -tmax should be set zeroo y=fftshift(y) # make y start with y(t=0) Y=fftshift(fft(y))/64.0

w=linspace(-pi*fmax,pi*fmax,65);w=w[:-1] figure() subplot(2,1,1) plot(w,abs(Y),lw=2) xlim([-8,8]) ylabel(r"$|Y|$",size=16) title(r"Spectrum of $\sin\left(\sqrt{2}t\right)\times w(t)$") grid(True) subplot(2,1,2) plot(w,angle(Y),’ro’,lw=2) xlim([-8,8])

ylabel(r"Phase of $Y$",size=16) xlabel(r"$\omega$",size=16) grid(True) savefig("fig10-6.png")

show()


 

Compare to our first plot and you can see that the magnitude is greatly improved. We√ still have a peak that is two samples wide. But that is because 2 lies between 1 and 2, which are the two fourier components available. If we use four times the number of points we should get better results.

12 heg7 12i≡ from pylab import * t=linspace(-4*pi,4*pi,257);t=t[:-1] dt=t[1]-t[0];fmax=1/dt n=arange(256) wnd=fftshift(0.54+0.46*cos(2*pi*n/256)) y=sin(sqrt(2)*t) # y=sin(1.25*t) y=y*wnd

y[0]=0 # the sample corresponding to -tmax should be set zeroo y=fftshift(y) # make y start with y(t=0) Y=fftshift(fft(y))/256.0 w=linspace(-pi*fmax,pi*fmax,257);w=w[:-1] figure() subplot(2,1,1) plot(w,abs(Y),’b’,w,abs(Y),’bo’,lw=2) xlim([-4,4]) ylabel(r"$|Y|$",size=16) title(r"Spectrum of $\sin\left(\sqrt{2}t\right)\times w(t)$") grid(True) subplot(2,1,2) plot(w,angle(Y),’ro’,lw=2) xlim([-4,4]) ylabel(r"Phase of $Y$",size=16)

xlabel(r"$\omega$",size=16) grid(True) savefig("fig10-7.png") show()


 

Is it better? Well it is quite a bit better since we are now zoomed in and see a lot more detail. But why is it not just a single peak? The reason for that is w(t). Multiplication in time is convolution in frequency and vice versa. So by multiplying with w(t), we got rid of the 1/f decay. But the delta function is now replaced by the shape of the DFT of w[n]. That gives us a factor of two broadening over the peak when there is no window, which is why we still see a peak whose width is two samples.√

Note that it is not because 2 is between 1.25 and 1.5. To verify, there is an alternate function in the above code, namely sin(1.25t). But this gives a broad peak as well. That is because of w[n].

The Assignment
1.   Work through the example codes and understand them.

2.   Consider the function cos3(ω0t). Obtain its spectrum for ω0 = 0.86 with and without a Hamming window.

3.   Write a program that will take a 128 element vector known to contain cos(ω0t +δ) for arbitrary δ and 0.5 < ω0 < 1.5. The values of t go from −π to π. You have to extract the digital spectrum of the signal, find the two peaks at ±ω0, and estimate ω0 and δ.

4.   Suppose the data in Q3 has added “white gaussian noise”. This can be generated byrandn() in python. The extent of this noise is 0.1 in amplitude (i.e., 0.1*randn(N), where N is the number of samples). Repeat the problem and find the ω0 and δ

5.   Plot the DFT of the function

 

for t going from −π to π in 1024 steps. This is known as a “chirped” signal, and its frequency continuously changes from 16 to 32 radians per second. This also means that the period is 64 samples near −π and is 32 samples near +π.

6.   For the same chirped signal, break the 1024 vector into pieces that are 64 sampleswide. Extract the DFT of each and store as a column in a 2D array. Then plot the array as a surface plot to show how the frequency of the signal varies with time.

This is new. So far we worked either in time or in frequency. But this is a “timefrequency” plot, where we get localized DFTs and show how the spectrum evolves in time.