% Name: butt_design_bi
% This program is used to design the IIR digital lowpass filter
% using the bilinear transformation
% The analog lowpass filter is the Butterworth
% The specifications of the digital filter are:
% alphap: The maximal passband attenuation in dB
% alphas: The minimum stopband attenuation in dB
% wp: The passband cutoff frequency 
% ws: The stopband cutoff frequency
% N:  The order of the lowpass filter
clear
format short
alphap=3;%input('Type in the data ap=:');
alphas=15;%input('Type in the data as=:');
wp=0.1*pi;%input('Type in the passband cutoff frequency wp=:');
ws=0.3*pi;%input('Type in the stopband cutoff frequency ws=:');
T=1;
omegap=2*tan(wp/2)/T;
omegas=2*tan(ws/2)/T;
%==========================================================================
ksp=sqrt((10^(0.1*alphap)-1)/(10^(0.1*alphas)-1));
lamdsp=omegas/omegap;
N=ceil(-log(ksp)/log(lamdsp));
omegac=omegas*(10^(0.1*alphas)-1)^(-1/(2*N));
fprintf('\n***The order N of Butterworth is %2.0f\n',N);
%===========================================================================
wn=wp/pi;
[z,p,k]=buttap(N);
[b,a]=zp2tf(z,p,k);
fprintf('\n***The coefficient vectors of Butterworth are %2.0f\n');
b
a
%==========================================================================
ba=b;aa=a;
for m=1:N;
    aa(m+1)=a(m+1)*omegac^(m);
    ba(m+1)=b(m+1)*omegac^(m);
end
fprintf('\n***The non-normlized coefficient vectors of Butterworth are %2.0f\n');
ba,aa
%==========================================================================
fprintf('\n***The coefficient vector of IIR digital filter in direct form are %2.0f\n');
[Bz,Az]=bilinear(ba,aa,1/T)
fprintf('\n***The coefficient vector of IIR digital filter in parallel form are %2.0f\n');
[cz,bz,az] = dir2par(Bz,Az)
%==========================================================================

[H,w]=freqz(Bz,Az);H0=abs(H(1));
%==========================================================================
subplot(211)
plot(w/pi,abs(H));
title('Amplitude response of the IIR digital filter'),ylabel('|H|'),axis([0,1,0,1.1]),
grid on
subplot(212)
plot(w/pi,20*log(abs(H)/H0));
ylabel('dB'),axis([0,1,-80,0]),
xlabel('Frequency (pi)'),grid on
%==========================================================================