%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% ECE 410
%%% Examples for lecture on 10/31/2001
%%% 
%%% Author:  K.E. Wage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%---------------------------------------------------------------------
%%% Plots of the window functions

N=25;					% window length

%%% Get the window functions
%%%  -- Normalize so that peak of transform is 1

wr=boxcar(N);  wr=wr/sum(wr);		% Rectangular 
wbart=bartlett(N);  wbart=wbart/sum(wbart);	% Bartlett 
whan=hanning(N); whan=whan/sum(whan);	% Hanning 
wham=hamming(N); wham=wham/sum(wham);	% Hamming 
wblack=blackman(N); wblack=wblack/sum(wblack); % Blackman

%%% Get the transforms of the windows
Npts=1024;				% # of points in DFT

[Hwr,w]=fr_resp(wr,Npts);
[Hwbart,w]=fr_resp(wbart,Npts);
[Hwhan,w]=fr_resp(whan,Npts);
[Hwham,w]=fr_resp(wham,Npts);
[Hwblack,w]=fr_resp(wblack,Npts);

%%% Plot the windows 
figure(1)
ni=0:N-1;
plot(ni,wr,'-o',ni,wbart,'-x',ni,whan,'-s',ni,wham,'-d',ni,wblack,'-^')
axis([0 24 0 .1])
fixfont(14,14,12)
xlabel('n');
ylabel('w[n]');
title('Common windows')
legend('Rectangular','Bartlett','Hanning','Hamming','Blackman')
print -deps windows.ps

%%% Plot the transforms
figure(2)
orient tall

subplot(411)
h=plot(w/pi,20*log10(abs(Hwr)),'--m',w/pi,20*log10(abs(Hwbart)),'c-');
set(h(1),'LineWidth',1);
set(h(2),'LineWidth',1);
axis([-1 1 -100 0])
fixfont(14,14,12)
ylabel('20log_{10}|W(e^{j\theta})|')
legend('Rectangular','Bartlett',3)

subplot(412)
h=plot(w/pi,20*log10(abs(Hwr)),'--m',w/pi,20*log10(abs(Hwhan)),'c-');
set(h(1),'LineWidth',1);
set(h(2),'LineWidth',1);
axis([-1 1 -100 0])
fixfont(14,14,12)
ylabel('20log_{10}|W(e^{j\theta})|')
legend('Rectangular','Hanning',3)

subplot(413)
h=plot(w/pi,20*log10(abs(Hwr)),'--m',w/pi,20*log10(abs(Hwham)),'c-');
set(h(1),'LineWidth',1);
set(h(2),'LineWidth',1);
axis([-1 1 -100 0])
fixfont(14,14,12)
ylabel('20log_{10}|W(e^{j\theta})|')
legend('Rectangular','Hamming',3)

subplot(414)
h=plot(w/pi,20*log10(abs(Hwr)),'--m',w/pi,20*log10(abs(Hwblack)),'c-');
set(h(1),'LineWidth',1);
set(h(2),'LineWidth',1);
axis([-1 1 -100 0])
fixfont(14,14,12)
xlabel('Normalized frequency \theta/\pi')
ylabel('20log_{10}|W(e^{j\theta})|')
legend('Rectangular','Blackman',3)

print -deps window_transforms.ps
%---------------------------------------------------------------------
%%% Show the example 

f1=3;					% Sinusoid 1 at 3 Hz
f2=10;					% Sinusoid 2 at 10 Hz
Tw=2;					% 2 second interval
fs=100;					% Sampling frequency

[x1,t]=samp_sine(f1,Tw,fs);
[x2,t]=samp_sine(f2,Tw,fs);
x0=ones(size(x1));			% DC component
x=x0+x1+x2;				% add sinusoids and DC

figure(3)				% Plot time signal
orient tall
subplot(411)
plot(t,x);
xlabel('Time (seconds)')
ylabel('x(t)')

print -deps application_xt.ps

[Hx,f]=fr_resp(x,Npts,fs);		% Get frequency response

figure(4)				% Plot frequency response magnitude
orient tall
subplot(311)
plot(f,abs(Hx));
axis([-50 50 0 210])
xlabel('Frequency (Hz)')
ylabel('Fourier transform magnitude for x(t)')

print -deps application_Xmag.ps


figure(5)				% Zoom in
orient tall
subplot(311)
plot(f,abs(Hx));
axis([-15 15 0 210])
xlabel('Frequency (Hz)')
ylabel('Fourier transform magnitude for x(t)')

print -deps application_Xmag_zoom.ps

%-------------------------------------------------------------
%%% Now look at what happens when we move sinusoid 1 closer to 
%%%  DC value

f1_1hz=1;
f1_halfhz=.5;

[x1,t]=samp_sine(f2,Tw,fs);
[x2,t]=samp_sine(f1_1hz,Tw,fs);
[x3,t]=samp_sine(f1_halfhz,Tw,fs);
x0=ones(size(x1));			% DC component

x=x0+x1+x2;		% contains DC, 1 Hz and 10 Hz sinusoids
y=x0+x1+x3;		% contains DC, 0.5 Hz and 10 Hz sinusoids

[Hx,f]=fr_resp(x,Npts,fs);		% Get frequency responses
[Hy,f]=fr_resp(y,Npts,fs);

figure(6)				% Plot 1 Hz example
orient tall				% resolved!
subplot(311)
plot(f,abs(Hx));
axis([-15 15 0 210])
xlabel('Frequency (Hz)')
ylabel('Fourier transform magnitude')
title('x(t)=1+sin(2\pi(1)t)+sin(2\pi(10)t)')

print -deps example_1hz.ps

figure(7)
orient tall
subplot(311)				% Plot 1/2 Hz example
plot(f,abs(Hy));			% not resolved!
axis([-15 15 0 210])
xlabel('Frequency (Hz)')
ylabel('Fourier transform magnitude')
title('x(t)=1+sin(2\pi(0.5)t)+sin(2\pi(10)t)')

print -deps example_halfhz.ps