DFT Windowing

This is code I used in a youtube video on DFT windowing – see http://youtu.be/PYd-pzaZpYE

% Tutorial on why DFT windowing works.
% See http://youtu.be/V_dxWuWw8yM for intro video
% Code available at https://dadorran.wordpress.com. 

fs = 1000;
T = 1/fs;
N = 1000; % Create a signal with N samples = 1 second
n = 0:N-1; % sample numbers
t = n*T; % sampling times over a one second duration.
x = cos(2*pi*6.5*t);

x_zpad = [x zeros(1, 99000)]; %zero pad
X_ZPAD_mags= abs(fft(x_zpad));

subplot(2,1,1)
plot(n,x);
xlabel('Samples');
ylabel('Amplitude')
title('Time-Domain Signal to be analysed');

subplot(2,1,2)
bins = 0:length(x_zpad)-1;
plot(bins(1:1500), X_ZPAD_mags(1:1500))
xlabel('Bins');
ylabel('Magnitude')
title('First 1500 Bins of Magnitude Spectrum of Zero-Padded Signal');

% Plot magnitude spectrum against frequency in Hertz
fax = [0:100000-1]*fs/100000;
plot(fax(1:1500), X_ZPAD_mags(1:1500));
title('Magnitude Spectrum of Zero-Padded Signal');
xlabel('Frequency (Hertz)');
ylabel('Magnitude')

%window using a hann window function
han_win = hanning(N)';
x_zpad_win = [x.*han_win zeros(1, 99000)]; %zero pad
X_ZPAD_WIN_mags= abs(fft(x_zpad_win));
hold on
plot(fax(1:1500), X_ZPAD_WIN_mags(1:1500),'r');


%The zero padded signal can be considered as being a sinusoid
% which is multiplied by a rectangular winddow. The folloiwing
% code illustrates this.
Ndisp = 5000; %number of time-domain samples to display
subplot(3,1,1)
plot(x_zpad(1:Ndisp))
title(['Zero Padded Signal (first ' num2str(Ndisp) ' samples)'])
ylabel('Amplitude')
xlabel('Samples');
t = [0:100000-1]*T; % sampling times for 100000 samples.
sinusoid = cos(2*pi*6.5*t);
subplot(3,1,2)
plot(sinusoid(1:Ndisp))
title(['Sinusoidal Signal (first ' num2str(Ndisp) ' samples)'])
ylabel('Amplitude')
xlabel('Samples');
subplot(3,1,3)
win = [ones(1, N) zeros(1, 99000)];
%win = [hanning(N)' zeros(1, 99000)];
plot(win(1:Ndisp))
title(['Rectangular Window (first ' num2str(Ndisp) ' samples)'])
ylabel('Amplitude')
xlabel('Samples');
ylim([-0.5 1.5])

%When you multiply signals in the time-domain it is equivalent to 
% convolution in the frequency domain. Let's take a look at each of 
% these three waveforms in the frequency domain
%update fax to  show negative frequencies
fax = [0:100000-1]*fs/100000;
fax(end-49998:end) = fax(end-49998:end)-1000;
fax(50001) = NaN;
subplot(3,1,1)
plot(fax,abs(fft(x_zpad)))
title(['Double Sided Magnitude Spectrum of  Zero Padded Signal'])
ylabel('Magnitude')
xlabel('Frequency (Hz)');

subplot(3,1,2)
plot(fax,abs(fft(sinusoid)))
title(['Double Sided Magnitude Spectrum of  Sinusoidal Signal'])
ylabel('Magnitude')
xlabel('Frequency (Hz)');

subplot(3,1,3)
plot(fax, abs(fft(win)))
title(['Double Sided Magnitude Spectrum of' ...
    ' Rectangular Window function '])
ylabel('Magnitude')
xlabel('Frequency (Hz)');

% Determine the magnitude spectrums with less zero padding
% Plot magnitude spectrum against frequency in Hz rather than bins
% Observe that the lower res spectrums are a sampled version of the
% higher res ones
X_mags4000 = abs(fft([x zeros(1, 3000)]));
X_mags2000 = abs(fft([x zeros(1, 1000)]));
X_mags1000 = abs(fft(x));

figure
fax = [0:4000-1]*fs/4000;
plot(fax(1:60), X_mags4000(1:60));
hold on
plot(fax(1:60), X_mags4000(1:60),'.');
set(gca,'XTick', fax(1:4:60))
xlabel('Frequency in Hertz')
ylabel('Magnitude')

fax = [0:2000-1]*fs/2000;
plot(fax(1:30), X_mags2000(1:30),'rx');

fax = [0:1000-1]*fs/1000;
plot(fax(1:15), X_mags1000(1:15),'go');

fax = [0:100000-1]*fs/100000;
plot(fax(1:1500), X_ZPAD_mags(1:1500),'k');

%Each of the other magnitude spectrums are a subset or sampled
% version of the high resolution one.
%It should be appreciated from these examples that all DFT's can
% be thought of as being a sampled version of the DTFT (which 
% produces a continuous spectrum that could be sampled at any 
% frequency)
%Understanding the effects of windowing therefore requires an
% undertstanding of what happening for the high res case.



% Comparison of some common window functions
figure
N = 1000;
zpad_len = 99000;
L = N+zpad_len;
zpad = zeros(1,zpad_len);
hamming_freq = fftshift(abs(fft([hamming(N)', zpad])))/(N);
plot([-5 :.01: 5],hamming_freq(L/2-500:L/2+500),'r')
hold on
rect_freq = fftshift(abs(fft([rectwin(N)', zpad])))/(N);
plot([-5 :.01: 5],rect_freq(L/2-500:L/2+500),'g')
hann_freq = fftshift(abs(fft([hanning(N)', zpad])))/(N);
plot([-5 :.01: 5],hann_freq(L/2-500:L/2+500),'c')
blackman_freq = fftshift(abs(fft([blackman(N)', zpad])))/(N);
plot([-5 :.01: 5],blackman_freq(L/2-500:L/2+500),'k')
legend('Hamming', 'Rectangular','Hann','Blackman')
ylabel('Magnitude');
xlabel('Bin width (non zero padded)');