%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Code for DSP Lecture of March 19, 2005
%
%  Different parts are separated by "pause"
%   You need to press a key to proceed
%   to following part
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%% Study of Windowing
%% Available Windows in Matlab
figure(1);
N= 32;
w0= window(@rectwin,N);
w1= window(@bartlett,N);
w2= window(@hamming,N);
w3= window(@hann,N);
w4= window(@blackman,N);
w5= window(@gausswin,N);

n= 1:N;
plot(n,w0,'k', n,w1,'b', n,w2,'g', n, w3, 'r',n,w4,'m',n,w5,'c');
axis([1 N 0 1.1]);
legend('Boxcar','Bartlett','Hamming', 'Hanning','Blackman','Gaussian');
xlabel('Sample (n)');
ylabel('Amplitude');
title('Comparison of Windowing Functions in the Time Domain');
pause


%%%%%% Now let's look at their Fourier transforms
M= 1024;
fw0= fftshift(fft(w0,M));
fw1= fftshift(fft(w1,M));
fw2= fftshift(fft(w2,M));
fw3= fftshift(fft(w3,M));
fw4= fftshift(fft(w4,M));
fw5= fftshift(fft(w5,M));
figure(2);
n= 1:M;
plot(n,abs(fw0),'k', n,abs(fw1),'b', n,abs(fw2),'g', n,abs(fw3),'r',n,abs(fw4),'m',n,abs(fw5),'c');
legend('Boxcar','Bartlett','Hamming', 'Hanning','Blackman','Gaussian');
xlabel('Frequency (k)');
ylabel('Magnitude (abs scale)');
title('Comparison of Windowing Functions in the Frequency Domain');
figure(3);
semilogy(n,abs(fw0),'k', n,abs(fw1),'b', n,abs(fw2),'g', n,abs(fw3),'r',n,abs(fw4),'m',n,abs(fw5),'c');
legend('Boxcar','Bartlett','Hamming', 'Hanning','Blackman','Gaussian');
xlabel('Frequency (k)');
ylabel('Magnitude (log scale)');
title('Comparison of Windowing Functions in the Frequency Domain');
pause

%%%%%% Effect of Windowing on Spectral Analysis of Sinusoidal Signals
T= 1e-4;            % Sampling rate of 10k
n= (1:64)*T;        % Sampling points
A0= 1;
A1= 0.75;
s1= A0*cos(2*pi/6.*1e4*n)+A1*cos(2*pi/3.*1e4*n);
s2= A0*cos(2*pi/14.*1e4*n)+A1*cos(4*pi/15.*1e4*n);
s3= A0*cos(2*pi/14.*1e4*n)+A1*cos(2*pi/12.*1e4*n);
s4= A0*cos(2*pi/12.*1e4*n)+A1*cos(4*pi/25.*1e4*n);
fs1= fftshift(fft(s1,M));
fs2= fftshift(fft(s2,M));
fs3= fftshift(fft(s3,M));
fs4= fftshift(fft(s4,M));
k=((1:M)*1./M-0.5)*2*pi;
figure(4);
subplot(2,2,1); plot(k,abs(fs1)); title('Fig. 13(b)'); xlabel('Frequency'); ylabel('Magnitude');
subplot(2,2,2); plot(k,abs(fs2)); title('Fig. 13(c)'); xlabel('Frequency'); ylabel('Magnitude');
subplot(2,2,3); plot(k,abs(fs3)); title('Fig. 13(d)'); xlabel('Frequency'); ylabel('Magnitude');
subplot(2,2,4); plot(k,abs(fs4)); title('Fig. 13(e)'); xlabel('Frequency'); ylabel('Magnitude');
pause

%%%%%%% Time dependent Fourier Transform - Spectrogram Example
T= 5e-5;                % Sampling Step
N= 4096;
freq= 0:(N-1);           %Instantaneous frequency 
n= (0:(N-1))*T;          % Sampling points
sig= cos(2*pi*freq.*n);  % signal with time dependent frequency content
%sig(1000:1500)=1;
figure(5);
subplot(3,1,1); plot(sig);
fsig= fftshift(fft(sig));
subplot(3,1,2); plot(abs(fsig));
L= 128;
spectrogram= zeros(L*4,N-L);
for k1=1:N-L-1 
  s1= sig(k1:k1+L-1);    
  fs1= fftshift(fft(s1,L*4));
  spectrogram(:,k1)= abs(fs1);
end
subplot(3,1,3); imagesc(spectrogram); colormap(gray);xlabel('Time'); ylabel('Frequency');title('Spectrogram');
pause


%%%%%%% Hilbert Transform  Sec. 11.4
n= (0:127)*5e-3;
asig= 10.*exp(j*2*pi*10*n);
rsig= real(asig);
frsig= fftshift(fft(rsig));
figure(6);
subplot(2,2,1); plot(n, rsig);
subplot(2,2,2); plot(n, abs(frsig));
hsig= hilbert(rsig);
fhsig= fftshift(fft(hsig));
subplot(2,2,3); plot(n, real(hsig),'b', n, imag(hsig),'r');  
subplot(2,2,4); plot(n, abs(fhsig));












