% Impulse Invariance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% specifications wp = 0.2*pi;
ws = 0.3*pi;
d1 = ??????
d2 = ??????
%parameter calculation
N = ???????
N = ceil(N);
%equation 7.17a
Wc = ????????
%poles
k=0:2*N-1;
sk = Wc*exp(????????);
%get only the upper LHP poles: sigma_k <0 & jW_k >0
sk = sk(???????);
%Continuous time filter:
B = Wc^N; % Numerator
A1 = conv([1 -sk(1)],[1 -conj(sk(1))]);
A2 = ???????
A3 = ???????
A = ??????? % Denominator
%Impulse invariance method --> Transformation from analog to digital
[Bz,Az] = impinvar(B,A);
%Frequency response plot
[H,w] = freqz(????,????,512);
% Group delay plot
[Gd,w]= grpdelay(????,????,512);
figure(1);
subplot(311);plot(w/pi,20*log10(abs(H))); % Divide with pi for an easier check
subplot(312);plot(w/pi,abs(H));
subplot(313);plot(w/pi,Gd)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do it using prebuild functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fs = 10000; % sampling frequency in Hz Wp = wp*fs;
Ws = ws*fs;
[Nm,Wcm]=buttord(Wp,Ws,-20*log10(d1),-20*log10(d2),'s'); % Finds optimum N, Wc
[Bm,Am]=butter(Nm,Wcm,'s'); % Constructs the filter
[Bzm,Azm] = impinvar(Bm,Am,fs); % Analog -> Digital
[Hm,w] = freqz(Bzm,Azm,512);
[Gdm,w]= grpdelay(Bzm,Azm,512);
figure(2);
subplot(311);plot(w/pi,20*log10(abs(Hm))); % Divide with pi for easier check subplot(312);plot(w/pi,abs(Hm));
subplot(313);plot(w/pi,Gdm)




% Bilinear Transformation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wp = 0.2*pi; ws = 0.3*pi; d1 = 0.89125; d2 = 0.17783;
Td=1;
Wp = ??????????
Ws = ??????????
Wc = ????????
N = ???????
N = ceil(N);
%equation: 7.34b
Wc = ????????
%poles
k=0:2*N-1;
sk = Wc*exp(???????);
%get only the upper LHP poles: sigma_k <0 & jW_k >0
sk = sk(????????);
%Continuous time filter:
B = Wc^N;
A1 = conv([1 -sk(1)],[1 -conj(sk(1))]);
A2 = ??????????????
A3 = ??????????????
temp = conv(A2, A3);
A = conv(A1,temp));
[Bz,Az] = bilinear(????,????,1/Td);
[H,w] = freqz(????,????,512);
[Gd,w]= grpdelay(????,????,512);
figure(3);
subplot(311);plot(w/pi,20*log10(abs(H))); % Divide with pi for easier check axis([0 1-100 20]);
subplot(312);plot(w/pi,abs(H)); axis([0 1 0 1]); subplot(313);plot(w/pi,Gd);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do it using prebuild functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fs = 1; % sampling frequency in Hz : DEN EINAI SHMANTIKO STOIXEIO
Wp = 2*fs*tan(wp/2);
Ws = 2*fs*tan(ws/2);
[Nm,Wcm]=buttord(Wp,Ws,-20*log10(d1),-20*log10(d2),'s'); % Find optimum N, Wc
[Bm,Am]=butter(Nm,Wcm,'s'); % Construct filter
[Bzm,Azm] = bilinear(Bm,Am,fs); % Analog -> Digital
[Hm,w] = freqz(Bzm,Azm,512);
[Gdm,w]= grpdelay(Bzm,Azm,512);






% Butterworth and elliptic filters examples
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wp = 0.4*pi;
ws = 0.6*pi;
d1 = ??????
d2 = ??????
fs = 10000; % sampling frequency in Hz
Wp = ??????????
Ws = ??????????
% Butterworth
[Nm,Wcm]=buttord(???,???,???,???,'s');
[Bm,Am]=butter(Nm,Wcm,'s');
[Bzm,Azm] = bilinear(Bm,Am,fs);
[Hm,w] = freqz(Bzm,Azm,512);
[Gdm,w]= grpdelay(Bzm,Azm,512); figure(1);
subplot(311);plot(w/pi,20*log10(abs(Hm))); % Divide with pi for easier check axis([0 1-100 20]);
subplot(312);plot(w/pi,abs(Hm)); axis([0 1 0.99 1.010]); subplot(313);plot(w/pi,Gdm);
axis([0 1 0 25]);
figure(2);clf;
zplane(Bzm,Azm);
% Elliptic
[Nm,Wcm]=ellipord(???,???,???,???,'s');
[Bm,Am]=ellip(???,???,???,???,'s');
[Bzm,Azm] = bilinear(Bm,Am,fs);
[Hm,w] = freqz(Bzm,Azm,512);
[Gdm,w]= grpdelay(Bzm,Azm,512);
figure(3);
subplot(311);plot(w/pi,20*log10(abs(Hm))); % Divide with pi for easier check axis([0 1-100 20]);
subplot(312);plot(w/pi,abs(Hm)); axis([0 1 0 1]); subplot(313);plot(w/pi,Gdm); axis([0 1 0 25]);