%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% ECE 410
%%% Example for lecture on 10/24/2001
%%% 
%%% Author:  K.E. Wage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%---------------------------------------------------------------------
%% First design the lowpass filter using the window method
%% Specifications:  wp=0.2pi  
%%                  ws=0.35pi
%%                  Hamming window:  41st order (42 points)

wp=.2;                  % Filter specifications
ws=.35;
wc=(wp+ws)/2;

M=41;                     % 41st order filter
h=fir1(M,wc,hamming(M+1),'noscale'); % to meet specifications

%---------------------------------------------------------------------
%% Plots of the lowpass filter impulse and frequency responses

%%%%%%%% Impulse response
figure(1)
orient tall
subplot(311)
n=0:length(h)-1;		% vector of indices to plot against
stem(n,h)
axis([0 41 -.1 .3])
xlabel('n')
ylabel('h[n]')
title('Impulse response of filter h[n]')
fixfont(12,12,12)
print -deps imp_resp.ps

%%%%%%%% Frequency response magnitude -- raw output of FFT
Hfft=(fft(h));
figure(2)
orient tall
subplot(311)
plot(abs(Hfft),'o-')
axis([0 43 0 1.1])
xlabel('index')
ylabel('abs(fft(h))')
title('Magnitude of fft(h)')
fixfont(12,12,12)
print -deps fft_resp.ps

%%%%%%%% Frequency response via freqz
figure(3)
orient tall
subplot(311)
[Hfreqz,wfreqz]=freqz(h,1,2048,'whole');
plot(wfreqz,abs(Hfreqz))
axis([0 2*pi 0 1.1])
xlabel('Frequency \omega (radians)')
ylabel('abs(Hfreqz)')
title('Frequency response computed with [Hfreqz,wfreqz]=freqz(h,1,2048,''whole'')')
fixfont(12,12,12)
print -deps freqz_resp1.ps

figure(4)
orient tall
subplot(311)
[Hfreqz,wfreqz]=freqz(h,1,2048,'whole');
plot(wfreqz/pi,abs(Hfreqz))
axis([0 2 0 1.1])
xlabel('Frequency \omega/\pi')
ylabel('abs(Hfreqz)')
title('Frequency response with frequency axis normalized by \pi')
fixfont(12,12,12)
print -deps freqz_resp2.ps

%%%%%%%% Frequency response via fr_resp
figure(5)
orient tall
subplot(311)
[Hw,w]=fr_resp(h,2048);
plot(w/pi,abs(Hw))
axis([-1 1 0 1.1])
xlabel('Frequency \omega/\pi')
ylabel('abs(Hw)')
title('Frequency response computed with [Hw,w]=fr\_resp(h,2048)');
fixfont(12,12,12)
print -deps fr_resp1.ps

%%%%%%%% Frequency response via fr_resp with optional sampling frequency
figure(6)
orient tall
subplot(311)
fs=1000;				% Sampling rate = 1000 Hz
[Hf,f]=fr_resp(h,2048,fs);
plot(f,abs(Hf))
axis([-fs/2 fs/2 0 1.1])
xlabel('Frequency (Hz)')
ylabel('abs(Hf)')
title('Frequency response computed with [Hf,f]=fr\_resp(h,2048,1000)');
fixfont(12,12,12)
print -deps fr_resp2.ps

%%%%%%%% Frequency response via fr_resp:  dB plot
figure(7)
orient tall
subplot(311)
plot(f,10*log10(abs(Hf)))
axis([-fs/2 fs/2 -50 10])
xlabel('Frequency (Hz)')
ylabel('10log_{10}|Hf|')
title('Frequency response computed with [Hf,f]=fr\_resp(h,2048,1000)');
fixfont(12,12,12)
print -deps fr_resp3.ps

%---------------------------------------------------------------------
%% Set up the signal to be filtered:  x is a sinusoid buried in noise

%% Sinusoid parameters
f0=68;					% frequency
fs=1000;				% sampling rate
Tw=1-1/fs;				% length of signal = 1 second

[x,t]=samp_sine(f0,Tw,fs);		% Make sinusoid
xsig=.1*x;
Nsig=size(x);

%% Make the noise:  take white noise and highpass filter
%%  Highpass filter parameters:  stop=0.22pi   pass=0.35pi  

w=randn(Nsig);
[n,wn]=ellipord(.35,.22,.5,60);
[b,a]=ellip(n,.5,60,wn,'high');
wfilt=filter(b,a,w);

x=xsig+wfilt;

%% Plot the input signal

figure(8)
orient tall
subplot(411)
plot(t,xsig)
axis([0 1 -3 3])
%%xlabel('time (seconds)')
ylabel('x(t)')
title('Sinusoidal signal')
fixfont(12,12,12)

subplot(412)
plot(t,xsig)
axis([0 .1 -.2 .2])
%%xlabel('time (seconds)')
ylabel('x(t)')
title('CLOSEUP:  Sinusoidal signal')
fixfont(12,12,12)

subplot(413)
plot(t,wfilt)
axis([0 1 -3 3])
%%xlabel('time (seconds)')
ylabel('x(t)')
title('Highpass-filtered noise')
fixfont(12,12,12)

subplot(414)
plot(t,x)
axis([0 1 -3 3])
xlabel('time (seconds)')
ylabel('x(t)')
title('Sinusoid plus noise')
fixfont(12,12,12)

print -deps sinusoid_noise.ps

%---------------------------------------------------------------------
%% Filter the signal

y=filter(h,1,x);
figure(9)
orient tall
subplot(411)
plot(t,y)
xlabel('time (seconds)')
ylabel('y(t)')
title('Filtered signal')
fixfont(12,12,12)

print -deps filtered_signal.ps
%---------------------------------------------------------------------