#2、IIR滤波器的MATLAB实现

1、2、IIR滤波器的 MATLAB实现2.1 IIR滤波器的设计方法及原理IIR 数字滤波器是一种离散时间系统，其系统函数为:H(Z)k =0N1八akZ丄y(z)XcZ)假设 M N,当M> N寸,系统函数可以看作一个IIR的子系统和一个(M-N)的FIR子系统的级联。IIR数字滤波器的设计实际上是求解滤波器的系数 ak和bk ,它是数 学上的一种逼近问题，即在规定意义上(通常采用最小均方误差准则) 去逼近系 统的特性。如果在S平面上去逼近，就得到模拟滤波器；如果在 z平面上去逼近， 就得到数字滤波器。2.1.1用脉冲相应不变法设计IIR数字滤波器利用模拟滤波器来设计数字滤波器，也就是使

2、数字滤波器能模仿模拟滤波器 的特性，这种模仿可以从不同的角度出发。脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应ha(t)，即将ha(t)进行等间隔采样，使h(n)正好等于ha(t)的采样值，满足h(n)讥(nT)H(z)为h(n)的Z变换,式中,T是采样周期。如果令Ha(s)是ha(t)的拉普拉斯变换, 利用采样序列的Z变换与模拟信号的拉普拉斯变换的关系得X(z)A oO.oO11zstXa(S_jks)'XaT kMTk =aO" 2兀 k s-jk(1-1)图1-1由(1-1 )式，数字滤jImz波器的频率响应间的

3、关系为百-2 兀 kJz平面十则可看出，脉冲响应不变法将模拟滤波器的 S平面变换成数字滤波器的Z平面, 这个从s到z的变换z=esT是从S平面变换到Z平面的标准变换关系式。3皿T波器的频率响应和模拟-二/ T时瓦HaI k(1-2)这就是说，数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。正如 采样定理所讨论的，只有当模拟滤波器的频率响应是限带的， 且带限于折叠频率 以内时，即(1-3)才能使数字滤波器的频率响应在折叠频率以内重现模拟滤波器的频率响应, 而不产生混叠失真，即(1-4)H(e叫十a(j T I1 I I丿但是，任何一个实际的模拟滤波器频率响应都不是严格限带的，变换后就会产生周

4、期延拓分量的频谱交叠，即产生频率响应的混叠失真，如图7-4所示。这 时数字滤波器的频响就不同于原模拟滤波器的频响，而带有一定的失真。当模拟 滤波器的频率响应在折叠频率以上处衰减越大、越快时，变换后频率响应混叠失 真就越小。这时，采用脉冲响应不变法设计的数字滤波器才能得到良好的效果。H|Ha(jQ)|一图1-2脉冲图1-2脉冲响应不变法中的频响混叠现象对某一模拟滤波器的单,立冲激响应ha(t)进行采样，采样频率为fs ,若使fs脉冲响应不变法优缺点:增加，即令采样时间间隔t=1辱)减小，贝略统频率响应各周期延拓分量 之间相距更远，因而可减小频率H响应的混叠效应从以上讨论可以看法使得数字滤波器的单

5、位脉冲响应完全模仿模拟滤波器的单位冲激响应，也就是时域逼近良好，而且模拟频率Q和数 字频率3之间呈线性关系3 =Q To因而，一个线性相位的模拟滤波器(例如贝塞 尔滤波器)通过脉冲响应不变法得到的仍然是一个线性相位的数字滤波器。脉冲响应不变法的最大缺点是有频率响应的混叠效应。所以，脉冲响应不变法只适用于限带的模拟滤波器(例如，衰减特性很好的低通或带通滤波器)，而 且高频衰减越快，混叠效应越小。至于高通和带阻滤波器，由于它们在高频部分 不衰减，因此将完全混淆在低频响应中。如果要对高通和带阻滤波器采用脉冲响 应不变法，就必须先对高通和带阻滤波器加一保护滤波器，滤掉高于折叠频率以上的频率，然后再使用

6、脉冲响应不变法转换为数字滤波器。当然这样会进一步增加设计复杂性和滤波器的阶数。2.1.2 用双线性变换法设计IIR数字滤波器脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到Z平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-n /Tn /T之间，再用z=esT 转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-n /T n /T 一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平 面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换 性，也就消除了频谱混叠现象，

7、映射关系如图1-3所示。隊繃图1-3双线性变换的映射关系为了将S平面的整个虚轴一j-駐压o上，可以通过以下的正切变o2二 tan TJljlmzn /T到n /T段RezS平面式中,T仍是米样间隔。.q平面(1-5)当Q 1由-n /T经过0变化到n /T时，Q由-经过o变化到+*，也即映射 了整个j Q轴。将式（1-5）写成将此关系解析延拓到整个2 面和_电1平面，令j Q =s，j Q仁s1，则得再将S1平面通过以下标准变换关2系映射到fsT平面22襌戸丁 * 1从而得到S平面禾Tze平面的单值映身T关系为：2丿T 1+eT(1-6)(1-7)2 1 -z s1T 1 z彳T 21 s s

8、 z 2T_式（1-6）与式（1-71 是S平面与T Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换.2g 、小ta n |= j。T i2丿即S平面的虚轴映射到Z平面的单位圆。式（1-5）与式（1-6）的双线性变换符合映射变换应满足的两点要求。 首先,把z=ej °,可得(1-8)2 1-e%. 2s节吋巧其次，将s= c +j Q代入式（1-8 ），得匕由此看出，当(7 <0时，|Z|<1 ；当6 >0时，|Z|>1。也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平

9、面的单位圆上。因此，稳定的模拟滤波器经双线性变换 后所得的数字滤波器也一定是稳定的。双线性变换法优缺点双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混_2 t 皂 1叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个j Q轴单 值地对应于Z平面单位圆一周，即频率轴是单值变换关系。这个关系如式(1-8 )所示，重写如下：上式表明，S平面上Q与Z平面的3成非线性的正切关系，如图7-7所示。由图7-7看出，在零频率附近，模拟角频率Q与数字频率3之间的变换关系 接近于线性关系；但当Q进一步增加时，3增长得越来越慢，最后当时， 3终止在折叠频率3 =冗处，因而双线性变换就不

10、会出现由于高频部分超过折叠 频率而混淆到低频部分去的现象，从而消除了频率混叠现象。护=半tan飙)图1-4双线性变换法的频率变换关系 i靠频率的严重非线性关系而得到的，如式 弊之间的非线性变换关系，就产生了新的问但是双线性变换的这个特点是(1-8 )及图1-4所示。由于这种频题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤 波器，不再保持原频响应必须是分-日_二勻位了;Tt©其次，1非线性关系要求模拟滤波器的幅I|Jr I(这正疋-般典型的低通、高通、带通、带阻 的数字滤波器幅频响应o/示O，即某一频率段的幅频响应近似等于某一常数，不然变换所产生i CO1-5所

11、阻型滤波器的响应特性) 原模拟滤波器的幅频响应会有畸变，女n图argH）图1-5双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变， 这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变， 然后经变换后正 好映射到所需要的数字频率上。2.2 Butterworth低通滤波器设计推导滤波器最小阶数与设计指标的关系：滤波器阶数就是其系统函数的极点个数。为了避免滤波器的复杂程度与我们的要求不匹配，造成不必 要的成本浪费，我们在滤波器设计前先需要确定其合适的阶数，即满足设 计要求的

12、最小阶数。以Butterworth滤波器为例，推导其阶数的数学模型若给出滤波器的设计指标为：通带截止频率 p，阻带截止频率 S，通带最大纹波Rp(dB)，阻带最小纹波Rs(dB)。因为滤波器幅频特性为:H(j(2-1)其中为3dB截止频率，N为滤波器阶数。所以当;：； ：:p以及;I ： s时, 可得到：(2-2)然后由Rp(dB)与Rs(dB)的关系式可得到:1Rp =10lgRs =10igHjp1 H(j s由上面的(2-2)式和(2-3)式可以联立求得:，广1 +(±)2N =1010彳 c 二 閃Rs1+(虫)2N 勺。10%rpH(jdp)2 1Rp1010(2-3)(2

13、-4)通过上面的结果我们就可以求得滤波器阶数 N为:Ig1010 -1p 叱)Rp1010 -1(2-5)所以滤波器的最小阶数就是大于上式所求得的值的最小整数。IIR滤波器的流程框图：开读入数字滤波器技术指标将指标转换成归一化模拟低通滤波器的指标设计归一化的模拟低通滤波器阶数N和3db截止频率模拟域频率变换，将G(P)变换成模拟低通滤波器H(s) 用双线性变换法将H(s)转换成数字带通滤波器H(z) 输入信号后显示相关结果结束2.2.1冲激响应不变法设计Butterworth低通滤波器冲激响应不变法是对模拟滤波器的单位冲激响应h(t)等间距采样获得数字滤波器的单位冲激响应，由此得到数字滤波器的

14、系 统函数。设定设计指标：通带截止频率p=2000n，阻带截止频率s=3000n，通带最大纹波 Rp=3dB，阻带最小纹波 Rs=18dB, 采样频率Fs=10000Hz。输入信号：x(t) = sin(2兀 f|t) + 0.5cos(2 f2t)Matlab程序如下：wp=2000*pi;ws=3000*pi;fs=10000;Rp=3;Rs=18;N,w n=buttord(wp,ws,Rp,Rs,'s'); z,p,k=buttap(N);Bap,Aap=zp2tf(z,p,k);b,a=lp2lp(Bap,Aap,w n);bz,az=imp in var(b,a,f

15、s);figure(1);h,w=freqz(bz,a z, N,fs);subplot(2,1,1),plot(w,abs(h);title('Butterworth LPF 幅频特性'); xlabel('频率(Hz) ');ylabel('幅值(dB) ');grid on; subplot(2,1,2),plot(w,a ngle(h);title('Butterworth LPF 相频特性');xlabel('频率(Hz) ');ylabel('相位(degree) ');grid on;

16、 figure(2);f仁 1000;f2=4000;N1=100;dt=1/fs; n=O:N1-1;t=n*dt;x=si n(2*pi*f1*t)+0.5*cos(2*pi*f2*t);subplot(2,1,1),plot(t,x);title('输入信号'),grid on;y=filter(bz,a z,x); subplot(2,1,2),plot(t,y, 'r-'); title('输出信号'),grid on;得到的Butterworth低通滤波器幅频特性和相频特性曲线为:Butterworth LPF幅频特性5C0 1COO

17、 1500200D25003000350040004500频率陀)Butterworth LPF相濒特性5(里)凹豊500 1GOO 1500 2DOO 25003000350040004500频率(Hk)2 O甘巴0岂日毀通过语句N,wn二buttord(wp,ws,Rp,Rs,'s')可以得到 Matlab 估 算出的滤波器最小阶数为 N=6。用式(2-5)出来的N=5.1925,取 大于它的最小整数得到滤波器的最小阶数为 6,两者一致，验证 了该数学模型的可靠性。输入信号和经过滤波后的输出信号图像为：输入信号0.001 0.002 0.003 0.004 0.005 O.

18、OOG 0.007 0.008 0.0090.010F面是基于Butterworth模拟原型滤波器，使用双线性变换法设计数字低通滤波器的过程。设定与上面不同的 设计指标为：通带截止频率 p=0.2 n , 阻带截止频率s=0.3 n，通带最大纹波Rp=1dB阻带最大纹波Rs=15dB采样 频率Fs=20000Hz同时，为了检测滤波器的性能，我们设定一样的输入信号x(t)二 sin(2 址)0.5cos(2 f2t)，其中 f仁 1000Hz，f2=4000Hz，将该信号与通过滤波器之后产生的输出信号进行比较以测试滤波器的性能。Matlab程序如下：Wp=0.2*pi;Ws=0.3*pi;fs=

19、20000;T=1/fs;Rp=1;Rs=15;wp=2*ta n( Wp/2)/T;ws=2*ta n( Ws/2)/T;N,Wn=buttord(wp,ws,Rp,Rs,'s');z,p,k=buttap(N);Bap,Aap=zp2tf(z,p,k);b,a=lp2lp(Bap,Aap,Wn);bz,az=bili near(b,a,fs);figure(1);h,w=freqz(bz,a z,N,fs);subplot(2,1,1),plot(w,abs(h);title('Butterworth LPF 幅频特性');xlabel('频率(Hz) ');ylabel('幅值(dB) ');grid on;subplot(2,1,2),plot(w,a ngle(h);title('Butterworth LPF 相频特性');xlabel('频率(Hz) ');ylabel('相位(degree) ');grid on;figure(2);f1=1