MATLAB实现数字巴特沃斯高通IIR滤波器

1、目 录摘要1Abstract21 设计项目要求与说明32 系统设计32.1 设计思路32.2 设计方法对比42.3典型模拟滤波器比较52.4 设计步骤53 仿真程序的设计与调试63.1 数字域指标变换成模拟域指标63.2 数字域频率进行预畸变63.3 模拟滤波器的设计73.4 模拟滤波器变成数字滤波器83.5 理论计算数字滤波器的仿真114.程序调试中出现的问题125. 总结与体会13参考文献14附录一 总程序如下15附录二 设计数字滤波器函数总结17摘要此报告重点介绍了用双线性不变法设计IIR数字滤波器的基本流程，比较了各种设计方法的优缺点，总结了模拟滤波器的性能特征。最后以双线性不变法设计

2、了一个高通巴特沃斯IIR数字滤波器，介绍了设计步骤，然后在Matlab环境下进行了仿真与调试，实现了设计目标。关键词：Matlab 双线性不变法 IIR数字滤波器 巴特沃斯 高通AbstractThis report introduced with emphasis of the basic flow of designing the IIR digit filterby the bilinear political reform, compared with each kind of design methods good and bad points, summarized analog

3、filters performance characteristic. Finally design one by the bilinear political reform to pass Butterworth high IIR digit filter, introduced the design procedure, then has carried on the simulation and the debugging under the Matlab environment, has achieved the project objective. Key word: Matlab

4、bilinearity political reform the IIR digital filter pass high butterworthMatlab课程设计 设计一个数字巴特沃斯高通IIR滤波器1 设计项目要求与说明 课题要求设计一个IIR数字滤波器，高通，采用双线性变换法，用巴特沃斯实现，用matlab软件对其进行仿真与调试。本设计将先说明用双线性法设计IIR数字滤波器的原理，然后写出基于matlab的软件设计流程。在对设计进行调试，分析实验数据。2 系统设计此部分将详细介绍IIR数字滤波器的设计流程，比较各种设计方案的优劣。2.1 设计思路 IIR滤波器设计的主要方法是先设计低通

5、模拟滤波器，然后转换为高通、带通或带阻数字滤波器。对于其他如高通，带通，则通过频率变换转换为设计相应的高通，带通等。在设计的全过程的各个步骤，matlab都提供相应的工具箱函数，使得IIR数字滤波器设计变得非常简单。总的来说，我的设计思路主要有以下两种：思路一：从归一化模拟低通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型的数字滤波器。归一化模拟低通原型数字高,带通或带阻模拟高,带通或带阻 模拟域 冲激响应不变法频率变换 双线性变换法图2-1 先频率变换再离散思路二：先进行双线性变换，将模拟低通原型滤波器变换成数字低通滤波器；然后

6、在Z域内经数字频率变换为所需类型的数字滤波器。 归一化模拟低通原型数字高,带通或带阻数字原型低通 数字域 双线性变换法 频率变换图2-2 先离散再频率变换以上两种思路都可以，我最后选择了第一种思路进行设计，即先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型的数字滤波器。2.2 设计方法对比方案一：冲激响应不变法 冲激相应不变法是从时域出发，要求数字滤波器的激响应h(n)对应于模拟滤波器ha(t)的等间隔抽样，h(n)=ha(nT) ，其中T是抽样周期，因此时域逼近良好。优点：h(n)完全模仿模拟滤波器的单位抽样响应时域逼近良好线性相位模拟滤波

7、器转变为线性相位数字滤波器缺点：对时域的采样会造成频域的“混叠效应”，故有可能使所设计数字滤波器的频率响应与原来模拟滤波器的频率响应相差很大，不能用来设计高通和带阻滤波器。只适用于限带的低通、带通滤波器方案二：双线性变换法双线性变换法是从频域出发，使DF的频率响应与AF的频率响应相似的一种变换法。直接使数字滤波器的频率响应，逼近模拟滤波器的频率响应,进而求得H(z)。优点：避免了频率响应的混迭现象 在特定AF和特定DF处，频率响应是严格相等的，它可以较准确地控制截止频率的位置。 它是一种简单的代数关系，设计十分方便。缺点：除了零频率附近，与之间严重非线性，即线性相位模拟滤波器变为非线性相位数字

8、滤波器 要求模拟滤波器的幅频响应为分段常数型，不然会产生畸变 对于分段常数型AF滤波器，经双线性变换后，仍得到幅频特性为分段常数的DF.但在各个分段边缘的临界频率点产生畸变，这种频率的畸变，可通过频率预畸变加以校正。方案三：频率变换法设计思想：从归一化模拟低通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型的数字滤波器。 先进行双线性变换，将模拟低通原型滤波器变换成数字低通滤波器；然后在Z域内经数字频率变换为所需类型的数字滤波器。 综上所述，频率变换法是基于双线性变换法或冲激响应不变法（用于频带变换），此处选择双线性变换法主要是基于

9、要获得严格的频率响应，以及较准确地控制截止频率的位置。而由于此种方法是一种简单的代数关系，设计也十分方便。2.3典型模拟滤波器比较1、Butterworth巴特沃斯滤波器：它具有单调下降的幅频特性；即最平幅度。2、Chebyshev切比雪夫滤波器：在通带或阻带等波纹，可提高选择性。3.Bessel贝塞尔滤波器：在通带内有较好的线性相位特性。4.Ellipse椭圆滤波器：其选择性相对前三种是最好的。 此处选择巴特沃斯主要是想获得最平稳的幅频响应。而不计较相位特性，而用双线性变换法也会将线性相位变为非线性相位。2.4 设计步骤 如设计一个数字低通滤波器，其技术指标为： 通带临界频率fp ，通带内衰

10、减小于rp；阻带临界频率fs，阻带内衰减大于s；采样频率为FS将指标变为角频率 wp=fp*2*pi;ws= fs*2*pi; （式2-1）将数字滤波器的频率指标Wk由wk=(2/T)tan(Wk/2)转换为模拟滤波器的频率指标wk,由于是用双线性不变法设计，故先采取预畸变。; ; （式2-2）将高通指标转换为低通指标，进而设计高通的s域模型归一化处理 ； （式2-3）； （式2-4）； （式2-5）由式2-3，2-4，2-5计算出N，查表可得模拟低通滤波器的阶数，从而由下式确定模拟高通滤波器的参数。 （式2-6）3 仿真程序的设计与调试3.1 数字域指标变换成模拟域指标 其程序为：fp =

11、400 ; fs= 300; Rp = 1; Rs = 20; wp =fp*2*pi; ws =fs*2*pi;FS=1000;T=1/FS；程序执行结果为：wp=2.5133e+003 ws=1.8850e+003与实际计算结果相符。3.2 数字域频率进行预畸变 其程序为：wp2=2*tan(Wp/2)/T; ws2=2*tan(Ws/2)/T;经过预畸变，可以发现频率变为： wp2= 6.1554e+003ws2= 2.7528e+0033.3 模拟滤波器的设计其程序为%设计模拟滤波器N,Wn = buttord(wp2,ws2,Rp,Rs,s) z,p,k=buttap(N); %创建

12、Buttord低通滤波器原型Bap,Aap=zp2tf(z,p,k); %由零极点转换为传递函数的形式figure(1)freqs(Bap,Aap);%模拟低通滤波器的频率响应title(模拟滤波器（低通原型）的频率响应)Bbs,Abs=lp2hp(Bap,Aap,Wn); %模拟低通变高通figure(2)freqs(Bbs,Abs);title(模拟滤波器的频率响应)程序执行后可以发现其频率响应为： N=4，其波形如下图3-1 图3-1 模拟滤波器的频率响应 由上图分析可得：其符合高通的一般特征，与预期的效果一样。而在此条件下，Butterworth滤波器低通原型的波形如下图3-2。图3-

13、2模拟滤波器（低通原型）的频率响应在设计的过程中，涉及一个频率变换的问题，即将模拟低通原型变为高通，其函数及用法如下：b,a=lp2hp(Bap,Aap,Wn);功能：把模拟滤波器原型转换成截至频率为 Wn 的高通滤波器。其中，Bap，Aap分别为低通传递函数的分子向量和分母向量；b，a分别为高通传递函数的分子向量和分母向量。3.4 模拟滤波器变成数字滤波器其程序为：Bbz,Abz=bilinear(Bbs,Abs,FS); %用双线性变换法设计数字滤波器 freqz(Bbz,Abz,512,FS); 程序运行的结果为：如下图3-2，图3-3图3-2 数字滤波器的频率响应由于使用的是双线性不变

14、法设计的，其相位为非线性。此处主要是基于要获得严格的频率响应，以及较准确地控制截止频率的位置，故画出了详细的幅频响应。（如下图3-2）图3-3 详细的幅频响应分析该图可知其在0.6（即300Hz）处的衰减为40dB，而在0.8（即400Hz）处的衰减极小，应小于1dB。由此可见，此设计符合要求设计的参数。而在调试的过程中发现：通带衰减越小，滤波器的性能越好 阻带衰减越大，滤波器的性能越好其曲线也越陡峭，选择性越好，当然所用的滤波器阶数也越高。当阻带衰减变为40dB（之前为20dB），通带不变时，其波形如下图3-4。对比上图3-3可知，其在阻带临界频率处衰减变为了40dB，曲线变陡峭了。图3-4

15、 详细的幅频响应（阻带衰减为40dB）当通带变为5dB时，阻带不变时，其波形如下图3-5。对比图3-3可知，其在通带处的衰减变为了5dB，曲线平滑了一些。图3-5 详细的幅频响应（通带衰减为5dB）3.5 理论计算数字滤波器的仿真wp=0.8*pi;ws=0.6*pi;OmegaP=2*1000*tan(wp/2);OmegaS=2*1000*tan(ws/2);lamdas=OmegaP/OmegaS;N=0.5*log10(10.(20/10)-1)/(10.(1/10)-1)/log10(lamdas);%笔算的结果为N=3.6947；故取N=4%此处为计算高通的传递函数Wn= 4.88

16、90e+003az=0 0 0 0 1;bz=1 2.613 3.414,2.613,1;Bbs,Abs=lp2hp(az,bz,Wn)%用双线性不变法处理Bbz,Abz=bilinear(Bbs,Abs,1000);其运行结果为：N=3.6947；图形如图3-6图3-6 理论计算的滤波器的幅频响应 综上所述，本滤波器以四阶即实现了预期的设计目标：采样频率为1000Hz，通带临界频率fp =400Hz，通带内衰减小于1dB（p=1）；阻带临界频率fs=300Hz，阻带内衰减大于20dB（s=25），其在通带内的性能更好。4.程序调试中出现的问题1.在使用巴特沃斯函数获取其阶数时发现，调用函数时

17、有点思路要理清，对于其函数N,Wn = buttord(wp2,ws2,Rp,Rs,s)，若有后面的s，N=4；而缺省时，N=6；经查资料发现，当有s时，其表示的是模拟Butterworth滤波器，而缺省时是数字Butterworth滤波器。2.一个线性是不变系统的典型表达式有状态空间型、传递函数型、零极增益型、极点留数型共四种，在开始的设计过程中我发现在把零极点增益型向传递函数型的转换过程中多用了一步，即先把零极点增益型转换为状态空间型，进行相应的频率变换和双线性变换后，再由状态空间型转换为传递函数型，进而绘制出其频率响应图，但是经仿真分析后发现，性能并不理想，所以并没有采用。3. freq

18、s(Bbs,Abs)与freqz(Bbz,Abz,512,FS)刚开始时分得不是很清楚，后来发现前者是求模拟域频率响应的，而后者是求数字域频率响应的。后者的用法如下：h,w=freqz(b,a,n,whole);或h,f=freqz(b,a,n,whole,Fs);式中，b,a为数字滤波器分子和分母多项式的系数，n为复数频率的响应点数，为整数，最好为2的幂，缺省时为512；Fs为采样频率，单位Hz。如果给定该值，则f位置输出为频率Hz，若没有给定，则按角频率（Angular frequency）给定f的频率矢量；whole表示返回的频率f或w值包含z平面整个单位圆频率矢量，即02；缺省时，频率

19、f或w值包含z平面上半单位圆（0）之间等间距n个点频率矢量。h为复频率响应；w为n点频率向量（单位rad）；f为n点频率向量(Hz)。函数返回值缺省时，绘制幅频响应和相频响应图5. 总结与体会在课设之前，我对MATLAB软件，特别是滤波器设计中的函数基本上处于一种模糊状态。通过在学校的电子资源的期刊网上找了些论文资料，借阅图书，一点一滴的自学，以及和同学不断的交流，最后完成了这次课设，对滤波器的设计有了比较清楚的了解。在课程设计的过程中，我学到了很多东西，比如设计滤波器的一些基本函数的用法，各种模拟滤波器的特性，设计滤波器的一些基本方法。但更为重要的是，我对于解决一个问题的思路更加清晰，找到了

20、属于自己的方法。当然，在设计的过程中，不可能避免的遇到了很多问题，如刚开始思路比较混乱，没有明确的方向。主要是如何将理论计算的模型转换为仿真模型。因为在理论上，将低通转换成高通，一般是变换将高通频率特征转换成低通原型频率特征；而在软件设计中，是对其传递函数进行修改，即变换其z域的表达式，设计初期一直不知道如何将其联系起来。后来发现，其实变换传递函数，也就是变换频率特征，是将变换后的频率代入原低通模型，而后得到高通模型的。总的来说，这次课程设计让我对MATLAB有了更深刻的了解，对数字滤波器的设计流程有了大致的了解，掌握了一些设计滤波器的基本方法，提高了理论用于实践的能力，掌握了更多专业相关的使

21、用知识与技能。同时，也暴露了我很多的不足，在以后的学习中，将进一步发扬有点，克服缺点。参考文献1陈怀琛MATLAB在电子信息课程中的应用（第二版）电子工业出版社20062 郭仕剑MATLAB7.X数字信号处理人民邮电出版社2006年3 陈贵明用MATLAB语言处理数字信号与数字图像科学出版社2000年4 王家文MATLAB7.0图形图像处理国防工业出版社2006年5 苏金明MATLAB图形图像电子工业出版社20056朱铭锆DSP 应用系统设计M 电子工业出版社2002.7韩得竹MATLAB 电子仿真与应用M 国防工业出版社2001.8邹鲲MATLAB6. x 信号处理M 清华大学出版社2002

22、.附录一 总程序如下fp = 400 ; fs = 300; Rp = 1; %通带最大衰减Rp=1dBRs = 20; %阻带最小衰减Rs=20dBwp =fp*2*pi; %把数字域滤波器特征换成模拟滤波器特征ws =fs*2*pi;FS=1000;T=1/FS%归一化数字频率Wp=wp/(FS);Ws=ws/(FS);%频率预畸变：数字域频率wp2=2*tan(Wp/2)/T; %预畸变求滤波器通带临界频率ws2=2*tan(Ws/2)/T; %预畸变求滤波器阻带临界频率%设计模拟滤波器N,Wn = buttord(wp2,ws2,Rp,Rs,s) z,p,k=buttap(N); %创

23、建Buttord低通滤波器原型Bap,Aap=zp2tf(z,p,k); %由零极点转换为传递函数的形式figure(1)freqs(Bap,Aap);%模拟低通滤波器的频率响应title(模拟滤波器（低通原型）的频率响应)Bbs,Abs=lp2hp(Bap,Aap,Wn); %模拟低通变高通figure(2)freqs(Bbs,Abs);title(模拟滤波器的频率响应)%用双线性不变法变换成数字滤波器Bbz,Abz=bilinear(Bbs,Abs,FS); %双线性变换%求其频率响应figure(3)freqz(Bbz,Abz,512,FS);title(数字滤波器的频率响应)%详细显示

24、数字滤波器的幅频响应hw,w=freqz(Bbz,Abz,512); figure(4)plot(w/pi,20*log10(abs(hw);gridaxis(0,1,-200,10)title( Butterworth Type Highpass Digital Filter)xlabel(w/pi);ylabel(幅度(dB);%下面将笔算的结果仿真wp=0.8*pi;ws=0.6*pi;OmegaP=2*1000*tan(wp/2);OmegaS=2*1000*tan(ws/2);lamdas=OmegaP/OmegaS;N=0.5*log10(10.(20/10)-1)/(10.(1/

25、10)-1)/log10(lamdas);%笔算的结果为N=3.6947；故取N=4%此处为计算高通的传递函数Wn= 4.8890e+003az=0 0 0 0 1;bz=1 2.613 3.414,2.613,1;Bbs,Abs=lp2hp(az,bz,Wn)%用双线性不变法处理Bbz,Abz=bilinear(Bbs,Abs,1000);%画图hw,w=freqz(Bbz,Abz,512); figure(5)plot(w/pi,20*log10(abs(hw);gridaxis(0,1,-200,10)set(gca,XTickMode,manual,XTick,0,Ws/(pi),Wp

26、/(pi),1);gridset(gca,YTickMode,manual,YTick,-200,-Rs,-Rp,10);gridtitle( 理论计算的滤波器的幅频响应)xlabel(w/pi);ylabel(幅度(dB);附录二 设计数字滤波器函数总结1.n,Wn= buttord(Wp,Ws,Rp,Rs) 9 N6 K5 o: N6 Q 其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。cheby1函数的用法为： / , k. P7 P6 P) Y2.b,a=cheby1(n,Rp,Wn,

27、/ftype/) & j1 A8 4 g, y! K, i+ D4 V. A7 0 m* X5 u( T- B5 B在使用cheby1函数设计IIR滤波器之前，可使用cheblord函数求出滤波器阶数n和截止频率Wn。cheblord函数可在给定滤波器性能的情况下，选择契比雪夫I型滤波器的最小阶和截止频率Wn。3.cheblord函数的用法为： * H6 F3 U 0 p2 H3 yj& c- N; YCn,Wn=cheblord(Wp,Ws,Rp,Rs) : 2 q2 j0 n) H5 8 F$ i) & _/ f, A% r( |, a) V, _1 其中Wp和Ws分别是通带和阻带的拐角频

28、率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。4.butter函数的用法为： b,a=butter(n,Wn,/ftype/) 其中n代表滤波器阶数，Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。buttord函数可在给定滤波器性能的情况下，求出巴特沃斯滤波器的最小阶数n，同时给出对应的截止频率Wn。5.Butterworth模拟低通滤波器函数语句：z,p,k=buttap(n). % 其中n为滤波器的阶数，z、 p、k分别为零点、极点和增益6.bz,az=impinvar(b,a,Fs，Fp)式中，b,a

29、为模拟滤波器分子和分母多项式系数向量；Fs为采样频率（所滤波数据），单位Hz，缺省时为1Hz。Fp为预畸变频率（Prewarped frequency）,是一个“匹配”频率,在该频率上，频率响应在变换前后和模拟频率可精确匹配在MATLAB中，函数bilinear采用双线性变换法实现模拟s域至数字z域的映射，直接用于模拟滤波器变换为数字滤波器。其调用方式为： zd,pd,kd=bilinear(z,p,k,Fs) numd,dend=bilinear(num,den,Fs)式中，z,p分别为模拟滤波器零点、极点列向量；k为模拟滤波器的增益；Fs为采样频率，单位Hz。zd,pd,kd为数字滤波器的

30、零极点和增益。num,den分别为模拟滤波器传递函数分子和分母多项式系数向量，模拟滤波器传递函数具有下面的形式： numd和dend分别为数字滤波器传递函数分子和分母多项式系数向量。7.函数freqz用于求数字滤波器的频率响应，其调用格式为：h,w=freqz(b,a,n,whole);或h,f=freqz(b,a,n,whole,Fs);式中，b,a为数字滤波器分子和分母多项式的系数，n为复数频率的响应点数，为整数，最好为2的幂，缺省时为512；Fs为采样频率，单位Hz。如果给定该值，则f位置输出为频率Hz，若没有给定，则按角频率（Angular frequency）给定f的频率矢量；whole表示返回的频率f或w值包含z平面整个单位圆频率矢量，即02；缺省时，频率f或w值包含z平面上半单位圆（0）之间等间距n个点频率矢量。h为复频率响应；w为n点频率向量（单位rad）；f为n点频率向量(Hz)。函数返回值缺省时，绘制幅频响应和相频响应图8.impz用于产生数字滤波器的脉冲响应。调用格式为：h,t=impz(b,a,n,Fs)式中，b,a分别为滤波器分子和分