2、IIR滤波器的MATLAB实现.docx

2、IIR滤波器的MAT1.AB实现2.1 HR滤波器的设计方法及原理IIR数字滤波器是一种离散时间系统，其系统函数为:N(Z)X(Z)MH(Z)=T1-zJl=I假设MWN,当M>N时,系统函数可以看作一个HR的子系统和一个(M-N)的FIR子系统的级联。IIR数字滤波器的设计事实上是求解滤波器的系数4和bk,它是数学上的一种靠近问题，即在规定意义上(通常采纳最小均方误差准则)去靠近系统的特性。假如在S平面上去靠近，就得到模拟滤波器；假如在Z平面上去靠近，就得到数字滤波器。2.1.1 用脉冲相应不变法设计HR数字滤波器利用模拟滤波器来设计数字滤波器，也就是使数字滤波器能仿照模拟滤波器的特性，这种仿照可以从不同的角度动身。脉冲响应不变法是从滤波器的脉冲响应动身，使数字滤波器的单位脉冲响应序列h(n)仿照模拟滤波器的冲激响应ha(t),即将ha(t)进行等间隔采样，使h(n)正好等于ha(t)的采样值，满意h(n)=ha(nT)式中,T是采样周期。假如令Ha(三)是ha(t)的拉普拉斯变换，H(Z)为h(n)的Z变换,利用采样序列的Z变换与模拟信号的拉普拉斯变换的关系得X“=隹x<f')=迈4-苧)(1-1)则可看出，脉冲响应不变法将模拟滤波器的S平面变换成数字滤波器的Z平面,这个从S到Z的变换Z=是从S平面变换到Z平面的标准变换关系式。图IT脉冲响应不变法的映射关系由(IT)式，数字滤波器的频率响应和模拟滤波器的频率响应间的关系为Hd)=¥元H(I,="<x>.-211kJT(1-2)这就是说，数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。正如采样定理所探讨的，只有当模拟滤波器的频率响应是限带的，且带限于折叠频率以内时，即“(C)=0-=T2(1-3)才能使数字滤波器的频率响应在折叠频率以内重现模拟滤波器的频率响应，而不产生混叠失真，即2) .“第I6><11(1-4)但是，任何一个实际的模拟滤波器频率响应都不是严格限带的，变换后就会产生周期延拓重量的频谱交叠，即产生频率响应的混叠失真，如图7-4所示。这时数字滤波器的频响就不同于原模拟滤波器的频响，而带有肯定的失真。当模拟滤波器的频率响应在折叠频率以上处衰减越大、越快时，变换后频率响应混叠失真就越小。这时，采纳脉冲响应不变法设计的数字滤波器才能得到良好的效果。对某一模拟滤波器的单位冲激响应h(t)进行采样，采样频率为fs,若使fs增加，即令采样时间间隔(T=lfs)减小，则系统频率响应各周期延拓重量之间相距更远，因而可减小频率响应的混叠效应。脉冲响应不变法优缺点：从以上探讨可以看出，脉冲响应不变法使得数字滤波器的单位脉冲响应完全仿照模拟滤波器的单位冲激响应，也就是时域靠近良好，而且模拟频率。和数字频率3之间呈线性关系3二。因而，一个线性相位的模拟滤波器(例如贝塞尔滤波器)通过脉冲响应不变法得到的仍旧是一个线性相位的数字滤波器。脉冲响应不变法的最大缺点是有频率响应的混叠效应。所以，脉冲响应不变法只适用于限带的模拟滤波器(例如，衰减特性很好的低通或带通滤波器)，而且高频衰减越快，混叠效应越小。至于高通和带阻滤波器，由于它们在高频部分不衰减，因此将完全混淆在低频响应中。假如要对高通和带阻滤波器采纳脉冲响应不变法，就必需先对高通和带阻滤波器加一爱护滤波器，滤掉高于折叠频率以上的频率，然后再运用脉冲响应不变法转换为数字滤波器。当然这样会进一步增加设计困难性和滤波器的阶数。2.1. 2用双线性变换法设计IlR数字滤波器脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到Z平面是多值的映射关系所造成的。为了克服这一缺点，可以采纳非线性频率压缩方法,将整个频率轴上的频率范围压缩到-冗/T11/T之间，再用z=esT转换到Z平面上。也就是说,第一步先将整个S平面压缩映射到Sl平面的-冗/T11/T一条横带里；其次步再通过标准变换关系z=eslT将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消退了多值变换性，也就消退了频谱混叠现象，映射关系如图1-3所示。为了将S平面的整个虚轴j压缩到SI平面jQ1轴上的-11/T到兀/T段上，可以通过以下的正切变换实现C=亲an（空）（1-5）式中,T仍是采样间隔。当Ql由-n/T经过O改变到n/T时，。由-8经过。改变到+8,也即映射了整个JQ轴。将式（1-5）写成2/。储2_/解24=.*。2+"必772将此关系解析延拓到整个S平面和SI平面，令jQ=S,JQl=S1.则得2esjf2-e's',22JslT2l-e-sjT产2+产27%2Jue-sl2I-ZT7+z,Z=再将Sl平面通过以下标准变换关系映射到Z平面Z=ClT从而得到S平面和Z平面的单值映射关系为：(1-6)(1-7)式（1-6）与式（1-7）是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换.式（1-5）与式（1-6）的双线性变换符合映射变换应满意的两点要求。首先,把z=e*,可得(1-8)2-e-j,2（cos=-=J-tan=Tl+/"T，即S平面的虚轴映射到Z平面的单位圆。其次,将s=。+jQ代入式（1-8）,因此由此看出，当。<o时，Z<1;当。>o时，z>io也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也肯定是稳定的。双线性变换法优缺点双线性变换法与脉冲响应不变法相比，其主要的优点是避开了频率响应的混叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个J0轴单C2(=tan-T值地对应于Z平面单位圆一周，即频率轴是单值变换关系。这个关系如式(1-8)所示，重写如下：上式表明，S平面上0与Z平面的3成非线性的正切关系，如图7-7所示。由图7-7看出，在零频率旁边，模拟角频率Q与数字频率3之间的变换关系接近于线性关系；但当Q进一步增加时，3增长得越来越慢，最终当。一8时，3终止在折叠频率G)=Jl处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消退了频率混叠现象。图1-4双线性变换法的频率变换关系但是双线性变换的这个特点是靠频率的严峻非线性关系而得到的，如式(1-8)及图1-4所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必需是分段常数型的，即某一频率段的幅频响应近似等于某一常数(这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图1-5所JO图1-5双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所须要的数字频率上。2.2Butterworth低通滤波器设计推导滤波器最小阶数与设计指标的关系：滤波器阶数就是其系统函数的极点个数。为了避开滤波器的困难程度与我们的要求不匹配，造成不必要的成本奢侈，我们在滤波器设计前先须要确定其合适的阶数，即满意设计要求的最小阶数。以BUtterwOrth滤波器为例，推导其阶数的数学模型。若给出滤波器的设计指标为：通带截止频率阻带截止频率以，通带最大纹波Rp（dB）,阻带最小纹波氏（力5）。因为滤波器幅频特性为：H（.同=-i=1（2-1）/+（一产其中Q.为3dB截止频率，N为滤波器阶数。所以当G=Gp以刚好G=,可得到：(2-2)l”(yj=J1=、HUs=I1=产FF然后由Rp(dB)与RKdB)的关系式可得到：RP=IolgT”(闻w(j)=3TIO10121&=叫两/=回=前(2-3)由上面的(2-2)式和(2-3)式可以联立求得：1+(”产=()需-,"Kn(")N=号l+(k)2=0”110缸1GC(2-4)通过上面的结果我们就可以求得滤波器阶数N为:(2-5)所以滤波器的最小阶数就是大于上式所求得的值的最小整数。IIR避波器的流程框图:起先读入数字滤波器技术指标I将指标转换成归一化模拟低通滤波器的指标设计归一化的模拟低通滤波器阶数N和3db截止频率模拟域频率变换，将G(P)变换成模拟低通滤波器词用双线性变换法将H(三)转换成数字带通滤波器H(Z)输入信号后显示相关结果结束冲激响应不变法是对模拟滤波器的单位冲激响应h等间距采样获得数字滤波器的单位冲激响应，由此得到数字滤波器的系统函数。设定设计指标：通带截止频率Gp=200（M,阻带截止频率以二3000冗，通带最大纹波Rp=3dB,阻带最小纹波Rs=18dB,采样频率Fs=10000Hzo输入信号：Mz）=sin（2rr）+0.5CoS（2%右）Matlab程序如下：wp=2000*pi;ws=3000*pi;fs=10000;R=3;Rs=18;N,wn=buttord(wp,ws,Rp,Rs,s');z,p,k=buttap(N);Bap,Aap=zp2tf(z,p,k);b,a=lp21p(Bap,Aap,wn);bz,az=impinvar(b,a,fs);figure(l);h,w=freqz(bz,az,N,fs);subplot(2,1,1),plot(w,abs（三）);title(,Butterworth1.PF幅频特性)；XlabelC频率(Hz);ylabel('幅值(dB)');gridon；subplot(2,1,2),plot(w,angle（三）);title(,Butterworth1.PF相频特性')；XIabe1('频率(HZ);ylabel('相位(degree);gridon;figure(2);fl=1000;f2=4000;Nl=100;dt=l/fs;n=O:N1-1;t=n*dt;X=Sin(2*pi*fl*t)+0.5*cos(2*pi*f2*t);subplot(2,1,1),plot(t,x);titleC输入信号)gridon;y=filter(bz,az,x);sublot(2,l,2),lot(t,y,r,);title('输出信号'),gridon;Butterworth1.PF幅频特性-4050010001500200025003000350040004500频率(HZ)20-2>aaP)W®通过语句N,wn=buttord(wp,ws,Rp,Rs方)可以得到Matlab估算出的滤波器最小阶数为N=6。用式(25)出来的N=5.1925,取大于它的最小整数得到滤波器的最小阶数为6,两者一样，验证了该数学模型的牢靠性。输入信号和经过滤波后的输出信号图像为:00.0010.0020.0030.0040.0050.0060.0070.0080.0090.01下面是基于Butterworth模拟原型滤波器，运用双线性变换法设计数字低通滤波器的过程。设定与上面不同的设计指标为：通带截止频率纤=0.2冗，阻带截止频率牡二0.3%通带最大纹波RP=IdB,阻带最大纹波Rs=15dB,采样频率Fs=20000Hzo同时，为了检测滤波器的性能，我们设定一样的输入信号x()=Sin(2r)+0.5cos(2M0,其中fl=1000Hz,f2=4000Hz,将该信号与通过滤波器之后产生的输出信号进行比较以测试滤波器的性能。Matlab程序如下：Wp=0.2*pi;WS=O.3*pi;fs=20000;T=lfs;Rp=I;Rs=15;wp=2*tan(W2)T;ws=2*tan(Ws2)T;N,Wn=buttord(wp,ws,Rp,Rs,s,);z,p,k=buttap(N);Bap,Aap=zp2tf(z,p,k);b,a=lp21p(Bap,Aap,Wn);bz,az=bilinear(b,a,fs);figure(1);h,w=freqz(bz,az,N,fs);subplot(2,1,1),plot(w,abs（三）);title(,Butterworth1.PF幅频特性')；XIabelC频率(HZ)NylabelC幅值(dB);gridon;subplot(2,l,2),plot(w,angle（三）);title('Butterworth1.PF相频特性,);XIabeI(濒率(HZ)NylabelC相位(degree);gridon;figure(2);fl=1000;f2=4000;Nl=100;dt=l/fs;n=O:N1-1;t=n*dt;x=sin(2*pi*fl*t)+0.5*cos(2*pi*f2*t);subplot(2,1,1),plot(t,x);title。输入信号'),gridon;y=filter(bz,az,x);subplot(2,l,2),plot(t,y,r-,);title。输出信号),gridon;得到的Butterworth低通滤波器幅频特性和相频特性曲线为:由语句N,Wn=buttord(wp,ws,Rp,Rs,W')我们可以用Matlab得出滤波器的最小阶数N=6,与上面一个滤波器一样。再利用上面我们推导出的式(2-5),将滤波器各项指标代入可以算出N=4.7433,取大于它的最小整数得到滤波器的最小阶数为5。二者并不一样，但差距只有1,说明我们推导出的计算BUtterWOrth滤波器最小阶数的数学模型存在误差，但误差不大。又因为Matlab计算滤波器的最小阶数也是估算出来的，所以这个误差在接受范围之类，我们推导出的数学模型可以认为是正确的。输入信号和经过滤波后的输出信号图像为:输入信号211111O0.511.522.533.544.55通过比较两组幅频特性和相频特性曲线我们可以发觉，其实双线性变换法和冲激响应不变法在设计较为简洁的Butterworth低通滤波器时差别并不大，并且都拥有不错的设计效果。