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


%%%%%%%%%%%%%%%%%%%%%%% IIR Filter Design Examples
%% Butterworth Filter
N= 5;
wn= 0.4;
wn0= 0.2;
[b1,a1] = butter(N, wn);
[h1,w1]= freqz(b1,a1,1024,'whole');
figure(1);
subplot(2,2,1); plot(w1-pi,fftshift(abs(h1))); title('High-Pass Filter');
[b2,a2] = butter(N, wn,'high');
[h2,w2]= freqz(b2,a2,1024,'whole');
subplot(2,2,2); plot(w2-pi,fftshift(abs(h2))); title('High-Pass Filter');
[b3,a3] = butter(N,[wn0 wn]);
[h3,w3]= freqz(b3,a3,1024,'whole');
subplot(2,2,3); plot(w3-pi,fftshift(abs(h3))); title('Band-Pass Filter');
[b4,a4] = butter(N,[wn0 wn],'stop');
[h4,w4]= freqz(b4,a4,1024,'whole');
subplot(2,2,4); plot(w4-pi,fftshift(abs(h4))); title('Band-Stop Filter');
pause

%% Chebyshev Type 1
Rp= 0.5;
N= 5;
wn= 0.4;
wn0= 0.2;
[b1,a1] = cheby1(N, Rp, wn);
[h1,w1]= freqz(b1,a1,1024,'whole');
figure(2);
subplot(2,2,1); plot(w1-pi,fftshift(abs(h1))); title('High-Pass Filter');
[b2,a2] = cheby1(N, Rp, wn,'high');
[h2,w2]= freqz(b2,a2,1024,'whole');
subplot(2,2,2); plot(w2-pi,fftshift(abs(h2))); title('High-Pass Filter');
[b3,a3] = cheby1(N, Rp,[wn0 wn]);
[h3,w3]= freqz(b3,a3,1024,'whole');
subplot(2,2,3); plot(w3-pi,fftshift(abs(h3))); title('Band-Pass Filter');
[b4,a4] = cheby1(N, Rp,[wn0 wn],'stop');
[h4,w4]= freqz(b4,a4,1024,'whole');
subplot(2,2,4); plot(w4-pi,fftshift(abs(h4))); title('Band-Stop Filter');
pause

%% Chebyshev Type 2
Rs= 20;
N= 5;
wn= 0.4;
wn0= 0.2;
[b1,a1] = cheby2(N, Rs, wn);
[h1,w1]= freqz(b1,a1,1024,'whole');
figure(3);
subplot(2,2,1); plot(w1-pi,fftshift(abs(h1))); title('High-Pass Filter');
[b2,a2] = cheby2(N, Rs, wn,'high');
[h2,w2]= freqz(b2,a2,1024,'whole');
subplot(2,2,2); plot(w2-pi,fftshift(abs(h2))); title('High-Pass Filter');
[b3,a3] = cheby2(N, Rs,[wn0 wn]);
[h3,w3]= freqz(b3,a3,1024,'whole');
subplot(2,2,3); plot(w3-pi,fftshift(abs(h3))); title('Band-Pass Filter');
[b4,a4] = cheby2(N, Rs,[wn0 wn],'stop');
[h4,w4]= freqz(b4,a4,1024,'whole');
subplot(2,2,4); plot(w4-pi,fftshift(abs(h4))); title('Band-Stop Filter');
pause

%% Elliptic Function
Rp= 0.5;
Rs= 20;
N= 5;
wn= 0.4;
wn0= 0.2;
[b1,a1] = ellip(N, Rp,Rs, wn);
[h1,w1]= freqz(b1,a1,1024,'whole');
figure(4);
subplot(2,2,1); plot(w1-pi,fftshift(abs(h1))); title('High-Pass Filter');
[b2,a2] = ellip(N, Rp,Rs, wn,'high');
[h2,w2]= freqz(b2,a2,1024,'whole');
subplot(2,2,2); plot(w2-pi,fftshift(abs(h2))); title('High-Pass Filter');
[b3,a3] = ellip(N, Rp,Rs, [wn0 wn]);
[h3,w3]= freqz(b3,a3,1024,'whole');
subplot(2,2,3); plot(w3-pi,fftshift(abs(h3))); title('Band-Pass Filter');
[b4,a4] = ellip(N, Rp,Rs, [wn0 wn],'stop');
[h4,w4]= freqz(b4,a4,1024,'whole');
subplot(2,2,4); plot(w4-pi,fftshift(abs(h4))); title('Band-Stop Filter');
pause


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIR Filter Design
%%%% Standard Filter Type - Windowing Method
b = fir1(48,0.35 );
figure(5)
freqz(b,1,512); title('FIR1 - Lowpass');
b = fir1(48,[0.35 0.65]);
figure(6)
freqz(b,1,512); title('FIR1 - Bandpass');
pause

%%%% Arbitrary Filter Type - Windowing Method
f = [0 0.2 0.4 0.6 0.8 1]; m = [1 0.25 1 0.25 1 0];
b = fir2(24,f,m);               % change order to see different approximation of response
[h,w] = freqz(b,1,128);
figure(7);
plot(f,m,w/pi,abs(h))
legend('Ideal','fir2 Designed')
title('Comparison of Frequency Response Magnitudes')
pause

%%%%% FIR Using Least Squares Solution
F = [0 0.3  0.4 0.6  0.7 0.9]; 
A = [0  1   0  0  0.5 0.5];
b = firls(40,F,A,'hilbert');
figure(8)
for i=1:2:6, 
   plot([F(i) F(i+1)],[A(i) A(i+1)],'b'), hold on
end
[H,f] = freqz(b,1,512,2);
plot(f,abs(H),'r'), grid on, hold off
legend('Ideal','firls Design')

