Bad result plotting windowing FFT

Question

im playing with python and scipy to understand windowing, i made a plot to see how windowing behave under FFT, but the result is not what i was specting. the plot is:

the middle plots are pure FFT plot, here is where i get weird things. Then i changed the trig. function to get leak, putting a 1 straight for the 300 first items of the array, the result:

the code:

sign_freq=80
sample_freq=3000
num=np.linspace(0,1,num=sample_freq)
i=0
#wave data:
sin=np.sin(2*pi*num*sign_freq)+np.sin(2*pi*num*sign_freq*2)
while i<1000:
    sin[i]=1
    i=i+1


#wave fft:
fft_sin=np.fft.fft(sin) 
fft_freq_axis=np.fft.fftfreq(len(num),d=1/sample_freq)

#wave Linear Spectrum (Rms)
lin_spec=sqrt(2)*np.abs(np.fft.rfft(sin))/len(num)
lin_spec_freq_axis=np.fft.rfftfreq(len(num),d=1/sample_freq)

#window data:
hann=np.hanning(len(num))

#window fft:
fft_hann=np.fft.fft(hann) 

#window fft Linear Spectrum:
wlin_spec=sqrt(2)*np.abs(np.fft.rfft(hann))/len(num)

#window + sin
wsin=hann*sin

#window + sin fft:
wsin_spec=sqrt(2)*np.abs(np.fft.rfft(wsin))/len(num)
wsin_spec_freq_axis=np.fft.rfftfreq(len(num),d=1/sample_freq)

fig=plt.figure()
ax1 = fig.add_subplot(431)
ax2 = fig.add_subplot(432)
ax3 = fig.add_subplot(433)
ax4 = fig.add_subplot(434)
ax5 = fig.add_subplot(435)
ax6 = fig.add_subplot(436)
ax7 = fig.add_subplot(413)
ax8 = fig.add_subplot(414)

ax1.plot(num,sin,'r')
ax2.plot(fft_freq_axis,abs(fft_sin),'r')
ax3.plot(lin_spec_freq_axis,lin_spec,'r')
ax4.plot(num,hann,'b')
ax5.plot(fft_freq_axis,fft_hann)
ax6.plot(lin_spec_freq_axis,wlin_spec)
ax7.plot(num,wsin,'c')
ax8.plot(wsin_spec_freq_axis,wsin_spec)

plt.show()

EDIT: as asked in the comments, i plotted the functions in dB scale, obtaining much clearer plots. Thanks a lot @SleuthEye !

Answer 1

It appears the plot which is problematic is the one generated by:

ax5.plot(fft_freq_axis,fft_hann)

resulting in the graph:

instead of the expected graph from Wikipedia.

There are a number of issues with the way the plot is constructed. The first is that this command essentially attempts to plot a complex-valued array (fft_hann). You may in fact be getting the warning ComplexWarning: Casting complex values to real discards the imaginary part as a result. To generate a graph which looks like the one from Wikipedia, you would have to take the magnitude (instead of the real part) with:

ax5.plot(fft_freq_axis,abs(fft_hann))

Then we notice that there is still a line striking through our plot. Looking at np.fft.fft's documentation:

The values in the result follow so-called “standard” order: If A = fft(a, n), then A[0] contains the zero-frequency term (the sum of the signal), which is always purely real for real inputs. Then A[1:n/2] contains the positive-frequency terms, and A[n/2+1:] contains the negative-frequency terms, in order of decreasingly negative frequency. [...] The routine np.fft.fftfreq(n) returns an array giving the frequencies of corresponding elements in the output.

Indeed, if we print the fft_freq_axis we can see that the result is:

[ 0.  1.  2. ..., -3. -2. -1.]

To get around this problem we simply need to swap the lower and upper parts of the arrays with np.fft.fftshift:

ax5.plot(np.fft.fftshift(fft_freq_axis),np.fft.fftshift(abs(fft_hann)))

Then you should note that the graph on Wikipedia is actually shown with amplitudes in decibels. You would then need to do the same with:

ax5.plot(np.fft.fftshift(fft_freq_axis),np.fft.fftshift(20*np.log10(abs(fft_hann))))

We should then be getting closer, but the result is not quite the same as can be seen from the following figure:

This is due to the fact that the plot on Wikipedia actually has a higher frequency resolution and captures the value of the frequency spectrum as its oscillates, whereas your plot samples the spectrum at fewer points and a lot of those points have near zero amplitudes. To resolve this problem, we need to get the frequency spectrum of the window at more frequency points. This can be done by zero padding the input to the FFT, or more simply setting the parameter n (desired length of the output) to a value much larger than the input size:

N = 8*len(num)
fft_freq_axis=np.fft.fftfreq(N,d=1/sample_freq)
fft_hann=np.fft.fft(hann, N) 
ax5.plot(np.fft.fftshift(fft_freq_axis),np.fft.fftshift(20*np.log10(abs(fft_hann))))
ax5.set_xlim([-40, 40])
ax5.set_ylim([-50, 80])