%Design a IIR filter using bilinear transformation
%wp=0.2*pi,ws=0.3*pi,rp=0.25dB,rs=50dB
clear all;
wp=0.2*pi;ws=0.3*pi;
rp=0.25;rs=50;
Fs=200;   %sampling frequency
Ts=1/Fs;
Wp=2/Ts*tan(wp/2.);Ws=2/Ts*tan(ws/2.);%frequency conversion
[N,Wn]=buttord(Wp,Ws,rp,rs,'s');  %computing the least Degree of analogue filter
[z,p,k]=buttap(N);
[Bap,Aap]=zp2tf(z,p,k);
[b,a]=lp2lp(Bap,Aap,Wn);
[bz,az]=bilinear(b,a,Fs);
Nn=256;%the number of sample point 
figure(1)
[H,f]=freqz(bz,az,Nn,Fs);
subplot(2,1,1),plot(f,20*log10(abs(H)));
xlabel('Frequency/Hz');ylabel('Amplitude/dB');grid on;
subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)));
xlabel('Frequency/Hz');ylabel('Phrase/^o');grid on;
figure(2)
f1=15;f2=40;
dt=1/Fs;t=0:dt:3;
x=sin(2*pi*f1*t)+cos(2*pi*f2*t);
y=filtfilt(bz,az,x);
subplot(2,1,1),plot(t,x),title('Input signal');
subplot(2,1,2),plot(t,y),title('Output signal'),xlabel('Time/s');