%用脈沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器
%濾波器的原形用BUTTERWORTH濾波器
%首先定義數(shù)字濾波器的參數(shù)
wp=pi*0.2 ;
ws=pi*0.3 ;
ap=1  ;
as=15  ;
%將數(shù)字濾波器的參數(shù)轉(zhuǎn)為濾波器濾波器的參數(shù)
T=0.001   ;
Fs=1/T;
omegap=wp*Fs;
omegas=ws*Fs;
%計算模擬濾波器的階數(shù)和在ap=3dB時的omgeac
a=10^(abs(as)/10)-1;
b=10^(abs(ap)/10)-1;
c=log10(omegas/omegap);
N=ceil(log10(sqrt(a/b))/c);             %計算階數(shù)N的值（四舍五入）
if(ap==3)
    omegac=omegap;
else
    omegac=omegap/((10^(0.1*abs(ap))-1)^(1/(2*N)));
end                                     %計算3dB是的omegac
%計算極點p、零點z和增益k
[z0,p0,k0]=buttap(N);
z=z0*omegac;
p=p0*omegac;
k=k0*omegac^N;
%計算模擬濾波器傳遞函數(shù)的H(z)中b和a的值
ba=k*real(poly(z));
aa=real(poly(p));
%進(jìn)行脈沖響應(yīng)變換
[ra,pa,ca]=residue(ba,aa);               %用分部分式計算模擬濾波器的C，P和R
                                         %數(shù)字濾波器的脈沖采樣序列為：
                                         %h(n)=h(t)*T=T*∑R*e^(P*T*n)
                                         %H(z)=∑(T*R)/(1-e^(P*T*n)*z^-1)
pd=exp(pa*T);
[bd,ad]=residuez(T*ra,pa,ca);            %得出數(shù)字濾波器傳遞函數(shù)的分子和分母系數(shù)
bd=real(bd);                             %real()消除虛部
ad=real(ad);


%impz(bd,ad);
%[h,w]=freqz(bd,ad);
%subplot(322);
%plot(w,20*log10(abs(h)));
%subplot(323);
%plot(w,angle(h)*180/pi);
%subplot(324);
zplane(bd,ad);
%gd=grpdelay(bd,ad);
%subplot(325);
%plot(gd)

%[H,w]=freqz(bd,ad);
%subplot(221),plot(w,abs(H))
%subplot(222),plot(w,angle(H))
%subplot(223),zplane(bd,ad)