%
% Example of resolution limit of 
% DFT.  Demonstrates use of fft in MATLAB,
% spectral leakage, zero padding, and frequency 
% resolution.  Discussed in class
% 
%  R Kakarala
%  UC Berkeley Extension
%  (c) 2003

freq=linspace(-0.5,0.5,1024);
W=fft(ones(1,32),1024);
subplot(2,2,1),plot(freq,abs(fftshift(W)));
axis([-0.5 0.5 0 35]);
title('Spectrum of 32-point rectangular window');
xlabel('freq (cycles/sample)'),ylabel('|W|');

% compute spectrum of two sinusoids with spacing 100% of 1/32
f1=sqrt(3)/32; f2=f1+1/32;
x=cos(2*pi*f1*(0:31))+cos(2*pi*f2*(0:31));
X=fft(x,1024);
subplot(2,2,2),plot(freq,abs(fftshift(X)));
title(sprintf('spacing 100 percent of 1/32'));
axis([-0.5 0.5 0 35]);
xlabel('freq (cycles/sample)'),ylabel('|X|');

% compute spectrum of two equal-amplitude sinusoids with spacing 75% of 1/32
f1=sqrt(3)/32; f2=f1+0.75/32;
x=cos(2*pi*f1*(0:31))+cos(2*pi*f2*(0:31));
X=fft(x,1024);
subplot(2,2,3),plot(freq,abs(fftshift(X)));
title(sprintf('spacing 75 percent of 1/32'));
axis([-0.5 0.5 0 35]);
xlabel('freq (cycles/sample)'),ylabel('|X|');

% compute spectrum of two equal-amplitude sinusoids with spacing 50% of 1/32
f1=sqrt(3)/32; f2=f1+0.5/32;
x=cos(2*pi*f1*(0:31))+cos(2*pi*f2*(0:31));
X=fft(x,1024);
subplot(2,2,4),plot(freq,abs(fftshift(X)));
title(sprintf('spacing 50 percent of 1/32'));
axis([-0.5 0.5 0 35]);
xlabel('freq (cycles/sample)'),ylabel('|X|');

% compute spectrum of two equal-amplitude sinusoids with spacing 50% of
% 1/32 but with 1024 samples
f1=sqrt(3)/32; f2=f1+0.5/32;
x=cos(2*pi*f1*(0:1023))+cos(2*pi*f2*(0:1023));
X=fft(x,1024);
figure,plot(freq,abs(fftshift(X)));
title(sprintf('spacing 50 percent of 1/32 but with 1024 samples'));
axis([-0.5 0.5 0 600]);
xlabel('freq (cycles/sample)'),ylabel('|X|');