MATLAB讲义第五、六章.docx

《MATLAB讲义第五、六章.docx》由会员分享，可在线阅读，更多相关《MATLAB讲义第五、六章.docx（14页珍藏版）》请在三一办公上搜索。

1、第五、六章IIR滤波器的设计滤波器结构：一、系统函数的表示法及其转换1、表示方法(1) 传递函数法若 h (z)二风)+ 妃)z 一1 + - + b( M) z -M1 + a (2)z t + a (N) z n则a=1 a(2) a(3)a(N)b=b(1) b(2)b(M)(2) 零极点增益法f1 (1 - zz 1)V一1若 H (z) = k +off (1 - zkz 1) k =0则零点向量Z=Z z2z;极点向量P=Z,Z2,，zN k为系统增益。(3) 部分分式法若 h (z) = 1 + ( + k (1) + k (M N +1) z -( m - n )1 一 p (

2、1) z 11 一 p (n) z n则极点向量p=p(1) p p(n)其对应系数向量r=r(1) r(2)r(n)余数多项式系数向量k=k(1) k k(M-N+1)(4) 二阶分式法把H(z)划成二阶因式H (z)=耻(z 顼 & o 心 T + & / 一2 k 1 k k 1 a +a zt +a z-2则其二阶因式为：b01 b11 b21 1 all a21b02 b12 b22 1 a21 a22Sos= b0N b1N b2N 1 a1N a2j2、各表示方法的转换：(1) 由传递函数转换为零极点增益(tf2zp)对应，由零极点增益转换为传递函数(zp2tf)调用方法：z,p

3、,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

4、=1 8 1 3; z,p,k=tf2zp(b,a) z =0-1.5000 + 1.6583i-1.5000 - 1.6583iP =-7.9216-0.0392 + 0.6141i-0.0392 - 0.6141ik =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.00001.0000-0.2500-0.12501

5、.00000.00009.00001.00001.00000.5000g =0.0625则其级联结构方程为:H (z) = 0.0625 1 3 z t + 2 z -21 - 0.25z-i - 0.125z-21 + 9 z -2 1 + z -1 + 0.5 z - 2由此可画出其结构图。例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.82844.00001.0000001.00002.82844.00001.0000001.0000-0.70710.2500

6、1.0000001.00000.70710.25001.000000G =1其传递函数应写成：H(z)= ( 1 + 2.83zT+4z-2 ) (1-2.83zT+4z-2) ( 1 + 0.71zT+0.25z-2 ) ( 1 0.71z-1+0.25z-2)滤波器的设计：一、模拟滤波器的设计1、巴氏模拟原型滤波器的设计：H (JQ)|用buttap函数，用于设计N阶归一化(Q c=1)巴氏模拟低通滤波器。调用格式: z0,p0,k0=buttap(N)其中：N为巴氏滤波器阶数z0,p0,0k为归一化滤波器的零点、极点、增益所设计的N阶巴氏模拟原型滤波器的归一化传递函数为：H (s)=k(

7、s p1)( s p 2)(s pN)若要设计未归一化(Q c尹1)巴氏低通滤波器，则要用Q c乘以p0、z0,而分子k0乘以Q n ，因为极点有n个。 C例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);

8、ylabel(|H(jw)|.八2); title(filer N=,num2str(N);grid;end;2、巴氏滤波器的最小阶数的求法:巴氏模拟原型滤波器的阶数越大，响应曲线在通带内越平缓，在阻带内的衰减越大。如何选择一个合适的阶数，使其刚好能实现系统的最优性能，而实现不复杂。阶数选择函数buttord：调用方式：N,OmegaC=buttord(OmegaP,OmegaS,Rp,Rs， s )N返回的滤波器最小阶数Wc3dB频率OmegaP , OmegaS通带、阻带截止频率，为归一化频率，范围为0,1，对应 n。Rp,Rs通带、阻带内的衰减s表示设计的是模拟滤波器例1：设通带、阻带截

9、止频率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.0000k01 p=p0*Omeg

10、aC;% 去归 化 k=k0*OmegaC八N k =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.00000.00000.00000.00000.00077.7094由结果可见，因为。相当大，0 n达1022数量级，致使向量Ca中的其它系数无法显示。因此最好仍用归一化系数表示滤波器传递函数。即归一化

11、传递函数为： b0=real(poly(z0)b0 =1 a0=real(poly(p0)a0 =1.00003.23615.23615.23613.23611.0000Ha =sssss()5 + 3.2361()4 + 5.2361()3 + 5.2361()2 + 3.2361() +13779237792377923779237792为检验此滤波器在Q p、Q s处的幅频特性，可用freqs函数，调用格式：H= freqs(b,a,w)b,N 一滤波器分子、分母系数向量w自变量的频率向量，不允许取标量，故至少要同时取两个频点若b、a取归一化值b0、a0，w也应归一化为w/OmegaC，

12、即H= freqs(b0,a0,w/ OmegaC)本例中，取OmegaP,OmegaS作为w,即w=OmegaP,OmegaS，并对其归一化为w/OmegaC b0=real(poly(z0) b0 a0=real(poly(p0)a0 =1.00003.23615.23615.23613.23611.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,a 转换的分子分母向量fs采

13、样频率，默认值为1Hzto误差容限，表示转换后的离散系统函数是否有重复的极点。实现思路：先把Ha(s)的分子分母向量ba,aa利用residue函数转换为极点留数Ra,pa及 直接项Ga,再把模拟滤波器的极点映射为数字极点pd=epa*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

14、,pd,Ga)%将极点留数形式转换为传递函数形式bd=real(bd),ad=real(ad)例1:利用脉冲响应不变法，把H(s) = S +1转换为数字滤波器， s 2 + 5s + 6T = 0.1。 b=1,1; a=1 5 6; T=0.1; bz,az=impinvar(b,a,1/T)0.1000-0.0897az =bz =得到H(z)为：、0.1 -0.0897z-1H (z)=1 -1.5595z -1 + 0.6065z - 21.0000-1.55950.6065 subplot(2,1,1),t=0:0.1:3; hat=impulse(b,a,t);%计算模拟系统的脉

15、冲响应ha(t),并乘以T plot(t,hat*T,red),hold on hnT=impz(bz,az,31); %计算数字系统的脉冲响应ha(nT),并以同样时间轴作图 stem(0.1火0:30,hnT,g) subplot(2,1,2),w=0:0.1:10*2火pi; HaS=freqs(b,a,w);plot(w,abs(HaS),red),hold on%模拟系统的频率特性 Ha(s) Hz=freqz(bz,az,w/Fs);plot(w/Fs,abs(Hz),g)% 数字系统的频率特性H(z) plot(w,abs(Hz),blu)%将数字频率放大Fs倍，与模拟频率的幅频

16、特性相比较(=2 n *Fs，范围为02n，所以要与模拟的相比，必须乘以Fs)2、双线性变换法米用 bilinearzd,pd,kd=bilinear(z,p,k,fs)把模拟滤波器的零极点模型转换为数字的零极点模型。bz, az=bilinear(bs,as,fs)把模拟滤波器的传递函数转换为数字的例：H =S +1 ，T = 1和0.1，用双线性变换法设计一个巴氏数字a s 2 + 5s + 6低通滤波器。 bs=1 1; as=1 5 6; T=0.1;Fs=1/T; bz,az=bilinear(bs,as,Fs)bz =0.04150.0040-0.0375T = 0.1 时az =

17、0.0415 + 0.0040z-i - 0.0375z-2az =H (z)=1 -1.5573z -1 + 0.6047 z -21.0000-1.55730.6047T=1;Fs=1/T;bz,az=bilinear(ba,aa,Fs)bz =0.15 + 0.1z-1 - 0.05 z -2 H (z)=1 + 0.2 z -10.15000.1000-0.0500az =1.00000.2000-0.0000 subplot(2,1,1),t=0:0.1:3;at*T, red),hold on:30,hn,*,g) hat=impulse(bs,as,t);plot(t,h hn=impz(bz,az,31);stem(0.1火0 title(ha(t) and h(n); subplot(2,1,2),w=0:0.1:10*2火pi; HaS=freqs(bs,as,w);plot(w,abs(HaS),red),hold on Hz=freqz(bz,az,w/Fs);plot(w/Fs,abs(Hz),g) plot(w,abs(Hz),blu) title(Ha(s) and H(z);