% E06_1.m                28apr21
% IIR & FIR filter design
clear all;
return;

%% Analog approximation functions
for N = 2:10
 [b,a] = butter(N,100,'s');      % order,omegaC
% [b,a] = cheby1(N,1,100,'s');   % order,ripple,omegaP
% [b,a] = cheby2(N,20,100,'s');  % order,ripple,omegaS
% [b,a] = ellip(N,1,20,100,'s'); % order,ripple1,ripple2,omegaC

H = freqs(b,a,0:299);
subplot(211); plot(abs(H)); axis([0 299 0 1.05]); hold on;
subplot(212); zplane(b,a); pause();  % or pzmap(b,a)

end;
subplot(211); hold off;

%% compare Butt. with Cheb.
clf;
N = 4;
[b1,a1] = butter(N,100,'s');      % order,omegaC
[b2,a2] = cheby1(N,3,100,'s');   % order,ripple,omegaC
Ome = 0:1e4;
H1 = 10*log(abs(freqs(b1,a1,Ome)));
H2 = 10*log(abs(freqs(b2,a2,Ome)));
semilogx(Ome,H1,Ome,H2); 
grid on;     % asymptotic slope: 6N dB/oct = 20N dB/dec

%% slide 6.1a.26  -  analog approximation
clear all;
OmeP = 2*pi*1000;
for OmeS = 2*pi*[5000,2500,1100,1010,1001];
   % down to OmeS = 2*pi*2500 order remains 3
   % note: for 1010 Matlab gives N=10, Octave gives N=9 --> 2dB error in response
Rp = 1;
Rs = 40;
[N,OmeN] = ellipord(OmeP,OmeS,Rp,Rs,'s')
[b,a] = ellip(N,Rp,Rs,OmeN,'s');
H = freqs(b,a,[0:2*pi*8000]);
figure(1); plot([0:2*pi*8000]/2/pi,20*log10(abs(H)));
grid on;
figure(2); zplane(b,a);
pause();
end;

%%  slide 6.1.35  -  butterworth notch filter
clear all; clf
om0 = (2*pi*60)/400; % z-plane specifications
omp1 = (2*pi*57)/400;  % (approx)
omp2 = (2*pi*63)/400;

Om0 = tan(om0/2); % warping for s-plane specs. (note T=2)
Omp1 = tan(omp1/2);
Omp2 = tan(omp2/2);
Band = Omp2 - Omp1;

[B,A] = butter(1,1,'s');   % prototype filter: N=1, Om_3dB=1
zplane(B,A); pause;  % beware this is the s-plane
[B1,A1]=lp2bs(B,A,Om0,Band);  % analog notch filter, N=2
zplane(B1,A1); pause;  % beware this is the s-plane
[b,a]=bilinear(B1,A1,.5);   % digital notch filter (Fs=1/T=.5)
zplane(b,a); pause;
freqz(b,a);
b,a
      %%% try also with FDATOOL, same specs.  

%% slide 7.1.21 - FIR filter design, windows
clear all; clf
M = 11;

wind = ones(1,2*M+1);
Wind = abs(fft(wind,1024)); Wind = Wind./max(Wind);
subplot(421); stem(-M:M,wind);
subplot(422); plot((0:512)/512,Wind(1:513)); axis([0 1 0 1]);

wind = hann(2*M+1);
Wind = abs(fft(wind,1024)); Wind = Wind./max(Wind);
subplot(423); stem(-M:M,wind);
subplot(424); plot((0:512)/512,Wind(1:513)); axis([0 1 0 1]);

wind = hamming(2*M+1);
Wind = abs(fft(wind,1024)); Wind = Wind./max(Wind);
subplot(425); stem(-M:M,wind);
subplot(426); plot((0:512)/512,Wind(1:513)); axis([0 1 0 1]);

wind = blackman(2*M+1);
Wind = abs(fft(wind,1024)); Wind = Wind./max(Wind);
subplot(427); stem(-M:M,wind);
subplot(428); plot((0:512)/512,Wind(1:513)); axis([0 1 0 1]);

%% slide 7.1.36  -  FIR filter design, smooth transition approach
clear all; clf
N = 30;
ome_p = .3*pi;
ome_s = .4*pi;
ome_c = (ome_p+ome_s)/2;
n = -N:N;
h = sin(ome_c*n) ./ (pi*n);     % sinc and, below, its spline approx.
h(N+1) = ome_c / pi;
figure(1); stem(n,h); hold on;
H_trunc = abs(fft(h,1024));
figure(2); plot(linspace(0,1,513),H_trunc(1:513)); axis 'tight'; grid on; hold on;

pause();
D = ome_s - ome_p;
P = 2;  % e.g. P = 1,2,3 ... 33
 h = ((2*sin(D*n/2/P)) ./ (D*n/P)).^P .* sin(ome_c*n) ./ (pi*n);
h(N+1) = ome_c / pi;
figure(1); stem(n,h,'d'); hold off;
H_smooth = abs(fft(h,1024));
figure(2); plot(linspace(0,1,513),H_smooth(1:513),'r'); hold off;

return;

