第五、六章 IIR 滤波器的设计滤波器结构：一、系统函数的表示法及其转换 1、表示方法 (1)传递函数法若NMz N a z a z M b z b b z H ----++++++=)()2(1)()2()1()(11则a=[1 a(2) a(3)… a(N)]b=[b(1) b(2)… b(M)] (2)零极点增益法若∏∏-=--=---=101101)1()1()(N k kM i iz zzz kz H则 零点向量 Z=[z 1 z 2 z M -1];极点向量 P=[z 1,z 2,…,z N -1] k 为系统增益。

(3)部分分式法 若)(1)1()1()(1)()1(1)1()(N M nz N M k k zn p n r z p r z H ----+-+++-++-=则极点向量 p=[p(1) p(2) … p(n)]其对应系数向量 r=[r(1) r(2) … r(n)]余数多项式系数向量 k=[k(1) k(2) … k(M-N+1)] (4)二阶分式法把H(z)划成二阶因式∏∏==----++++==Nk Nk k k k k k k k z zz z z H z H 112211022110)()(αααβββ 则其二阶因式为：b01 b11 b21 1 a11 a21 b02 b12 b22 1 a21 a22 Sos= …b0N b1N b2N 1 a1N a2N 2、各表示方法的转换：(1)由传递函数转换为零极点增益（tf2zp ） 对应，由零极点增益转换为传递函数（zp2tf ） 调用方法：[z,p,k]=tf2zp(b,a) [b,a]=zp2tf(z,p,k)――a,b 的长度要相等，不等的话要补零(2)由零极点增益转换为二次分式（zp2sos ） 对应，由二次分式转换为零极点增益（sos2zp ）调用方法：[sos,g]=zp2sos(z,p,k),g 为整个系统的增益，即H(z)=g*H1(z)*H2(z)…*HN(z)调用方法：sos2zp(z,p,k) =[sos,g],g 为整个系统的增益,默认为1。

(3)二次分式转换为传递函数(sos2tf) 调用方法：[b,a]=sos2tf(sos)对应，由传递函数转换为二次分式 (tf2sos)，调用方式： [sos,g]=tf2sos(b,a)例1：>> b=[1 3 5 0]; >> a=[1 8 1 3]; >> [z,p,k]=tf2zp(b,a) z =0 -1.5000 + 1.6583i -1.5000 - 1.6583i p =-7.9216 -0.0392 + 0.6141i -0.0392 - 0.6141i k = 1例2：差分方程：16y(n)+12y(n-1)+2y(n-2)-4y(n-3)-y(n-4)=x(n)-3x(n-1)+11x(n-2)-27x(n-3)+18y(n-4)程序：>> b=[1 -3 11 -27 18]; >> a=[16 12 2 -4 -1]; >> [sos,g]=tf2sos(b,a) sos =1.0000 -3.00002.0000 1.0000 -0.2500 -0.1250 1.0000 0.0000 9.0000 1.0000 1.0000 0.5000 g =0.0625则其级联结构方程为：21221215.0191125.025.012310625.0)(-------+++∙--+-∙=zz z z z z z z H 由此可画出其结构图。

例3：H(z)=1+16.0625z-4+z -8>> b=[1 0 0 0 16.0625 0 0 0 1]; >> [sos,G]=tf2sos(b,1) sos =1.0000 -2.8284 4.0000 1.0000 0 0 1.0000 2.8284 4.0000 1.0000 0 0 1.0000 -0.7071 0.2500 1.0000 0 0 1.0000 0.7071 0.2500 1.0000 0 0 G = 1其传递函数应写成：H(z)=（1+2.83z -1+4z -2）(1-2.83z -1+4z -2)（1+0.71z -1+0.25z -2）（1－0.71z -1+0.25z -2）滤波器的设计： 一、模拟滤波器的设计1、巴氏模拟原型滤波器的设计：NCj H 211)(⎪⎪⎭⎫ ⎝⎛ΩΩ+=Ω用buttap 函数，用于设计N 阶归一化（Ωc=1）巴氏模拟低通滤波器。

调用格式： [z0,p0,k0]=buttap(N)其中：N ――为巴氏滤波器阶数z0,p0,0k ――为归一化滤波器的零点、极点、增益 所设计的N 阶巴氏模拟原型滤波器的归一化传递函数为： )()2)(1()()()(pN s p s p s ks P s Z s H ---==若要设计未归一化（Ωc ≠1）巴氏低通滤波器，则要用Ωc 乘以p0、z0，而分子k0乘以ΩN c,因为极点有N 个。

例1：设计阶数为3，5，10，15的巴氏模拟原型滤波器。

并画出幅频响应曲线。

程序for i=1:4;switch icase 1N=3;case 2N=5;case 3N=10;case 4;N=15;end;[z,p,k]=buttap(N);[b,a]=zp2tf(z,p,k);[h,w]=freqs(b,a,n);Ah=abs(h);subplot(2,2,i),plot(w,Ah);axis([0 2 0 1]);xlabel('w/wc');ylabel('|H(jw)|.^2'); title(['filer N=',num2str(N)]);grid;end;2、巴氏滤波器的最小阶数的求法：巴氏模拟原型滤波器的阶数越大，响应曲线在通带内越平缓，在阻带内的衰减越大。

如何选择一个合适的阶数，使其刚好能实现系统的最优性能，而实现不复杂。

阶数选择函数buttord：调用方式：[N,OmegaC]=buttord(OmegaP,OmegaS,Rp,Rs，’s’)N――返回的滤波器最小阶数Wc――3dB频率OmegaP , OmegaS――通带、阻带截止频率，为归一化频率，范围为[0,1]，对应π。

Rp,Rs――通带、阻带内的衰减s――表示设计的是模拟滤波器例1：设通带、阻带截止频率fp=5kHz、fs=12kHz，通带、阻带最大衰减Rp=1dB,Rs=30dB，要求设计巴氏低通滤波器。

>> OmegaP=2*pi*5000;OmegaS=2*pi*12000;>> Rp=1;Rs=30;>> [N,OmegaC]=buttord(OmegaP,OmegaS,Rp,Rs,'s')%确定阶数NN =5OmegaC =3.7792e+004>> [z0,p0,k0]=buttap(N)%确定传递函数z0 =[]p0 =-0.3090 + 0.9511i-0.3090 - 0.9511i-0.8090 + 0.5878i-0.8090 - 0.5878i-1.00001>> p=p0*OmegaC;%去归一化>> k=k0*OmegaC^Nk =7.7094e+022>> p=p0*OmegaCp =1.0e+004 *-1.1678 + 3.5943i-1.1678 - 3.5943i-3.0575 + 2.2214i-3.0575 - 2.2214i-3.7792>> b=k*real(poly(z))b =7.7094e+022>> a=real(poly(p))a =1.0e+022 *0.0000 0.0000 0.0000 0.0000 0.0007 7.7094由结果可见，因为Ωc相当大，ΩN达1022数量级，致使向量ca中的其它系数无法显示。

因此最好仍用归一化系数表示滤波器传递函数。

即归一化传递函数为：>> b0=real(poly(z0))1>> a0=real(poly(p0)) a0 =1.0000 3.2361 5.2361 5.2361 3.2361 1.00001)37792(2361.3)37792(2361.5)37792(2361.5)37792(2361.3)37792(1)(2345+++++=ss s s s s H a 为检验此滤波器在Ωp 、Ωs 处的幅频特性，可用freqs 函数，调用格式：H ＝freqs(b,a,w)b,a ――滤波器分子、分母系数向量w ――自变量的频率向量，不允许取标量，故至少要同时取两个频点 若b 、a 取归一化值b0、a0，w 也应归一化为w/OmegaC ，即 H= freqs(b0,a0,w/ OmegaC )本例中，取OmegaP,OmegaS 作为w,即w=[OmegaP,OmegaS]，并对其归一化为w/OmegaC>> b0=real(poly(z0)) b0 = 1>> a0=real(poly(p0)) a0 =1.0000 3.2361 5.2361 5.2361 3.2361 1.0000>> w=[5000,12000]*2*pi/37792; >> H=freqs(b0,a0,w);键入freqs(b,a)，就可得出系统的频率特性二、IIR 数字滤波器的设计 1、冲激响应不变法采用impinvar 函数[bz,az]=impinvar(b,a,fs,tol)bz,az ――转换的分子分母向量 fs ――采样频率，默认值为1Hztol ――误差容限，表示转换后的离散系统函数是否有重复的极点。

实现思路：先把Ha(s)的分子分母向量ba,aa 利用residue 函数转换为极点留数Ra,pa 及直接项Ga,再把模拟滤波器的极点映射为数字极点pd=e pa*T ，得到H(z)的极点留数形式，再用residuez 把H （z ）转换为传递函数形式。

function [bd,ad]=impinvar_my(ba,aa,Fs)[Ra,pa,Ga]=residue(ba,aa);%把模拟滤波器系数向量变为极点留数形式 T=1/Fs;pd=exp(pa*T);Rd=Ra*T%将模拟极点留数转换为数字极点留数 [bd,ad]=residuez(Rd,pd,Ga)%将极点留数形式转换为传递函数形式 bd=real(bd),ad=real(ad)例1：利用脉冲响应不变法，把651)(2+++=s s s s H 转换为数字滤波器，T ＝0.1。