MATLAB实现数字巴特沃斯高通IIR滤波器(双线性变换法)

摘要 2Abstract 31设计项目要求与说明 42系统设计 2.1设计思路 42.2设计方法对比 52.3典型模拟滤波器比较 62.4设计步骤 73仿真程序的设计与调试 83.1数字域指标变换成模拟域指标 83.2数字域频率进行预畸变 83.3模拟滤波器的设计 3.4模拟滤波器变成数字滤波器 3.5理论计算数字滤波器的仿真 4.程序调试中出现的问题 5.总结与体会 参考文献 附录一总程序如下 附录二设计数字滤波器函数总结 2此报告重点介绍了用双线性不变法设计IIR数字滤波器的基本流程，比较了各种设计方法的优缺点，总结了模拟滤波器的性能特征。最后以双线性不变法设计了一个高通巴特沃斯IIR数字滤波器，介绍了设计步骤，然后在Matlab环境下进行了仿真与调试，实现了设计目标。数字滤波器巴特沃斯高通大学《数字信号处理》报告3ThisreportintroducedwithemphasisofthebasicflowofdesigningtheIRdigitbythebilinearpoliticalreform,comparedwitheachkindofdesignmethod'sgoodandbadpoints,summarizedanalogfilter'sperfobythebilinearpoliticalreformtopassButterworthhighIRdigitfilter,introducedtheMatlabenvironment,hasachievedtheprojectobjective.Keyword:MatlabbilinearitypoliticalreformtheIIRdigitalfilterpasshighbutterworth大学《数字信号处理》报告4——设计一个数字巴特沃斯高通1设计项目要求与说明课题要求设计一个IR数字滤波器，高通，采用双线性变换法，用巴特沃斯实现，用matlab软件对其进行仿真与调试。本设计将先说明用双线性法设计IR数字滤波器的原理，然后写出基于matlab的软件设计流程。在对设计进行调试，2系统设计此部分将详细介绍IR数字滤波器的设计流程，比较各种设计方案的优劣。IIR滤波器设计的主要方法是先设计低通模拟滤波器，然后转换为高通、带通或带阻数字滤波器。对于其他如高通，带通，则通过频率变换转换为设计相应使得IIR数字滤波器设计变得非常简单。总的来说，我的设计思路主要有以下两思路一：从归一化模拟低通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型归一化模拟低通原型冲激响应不变法双线性变换法大学《数字信号处理》报告5图2-1先频率变换再离散思路二：先进行双线性变换，将模拟低通原型滤波器变换成数字低通滤波器；然后在Z域内经数字频率变换为所需类型的数字滤波器。数字原型低数字原型低通归一化模拟低通原型或带阻数字域图2-2先离散再频率变换以上两种思路都可以，我最后选择了第一种思路进行设计，即先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型的数字滤波器。方案一：冲激响应不变法冲激相应不变法是从时域出发，要求数字滤波器的激响应h(n)对应于模拟滤波器ha(t)的等间隔抽样，h(n)=ha(nT),其中T是抽样周期，因此时域逼近良好。优点：①h(n)完全模仿模拟滤波器的单位抽样响应时域逼近良好②线性相位模拟滤波器转变为线性相位数字滤波器缺点：①对时域的采样会造成频域的“混叠效应”,故有可能使所设计数字滤波器的频率响应与原来模拟滤波器的频率响应相差很大，②不能用来设计高通和带阻滤波器。只适用于限带的低通、带通滤波器方案二：双线性变换法双线性变换法是从频域出发，使DF的频率响应与AF的频率响应相似的一种变换法。直接使数字滤波器的频率响应，逼近模拟滤波器的频率响应，进而求得H(z)。优点：①避免了频率响应的混迭现象②在特定AF和特定DF处，频率响应是严格相等的，它可以较准确地控制截止频率的位置。大学《数字信号处理》报告6③它是一种简单的代数关系，设计十分方便。缺点：①除了零频率附近，o与Ω之间严重非线性，即线性相位模拟滤波器变为非线性相位数字滤波器②要求模拟滤波器的幅频响应为分段常数型，不然会产生畸变③对于分段常数型AF滤波器，经双线性变换后，仍得到幅频特性为分段常数的DF.但在各个分段边缘的临界频率点产生畸变，这种频率的畸变，可通过频率预畸变加以校正。方案三：频率变换法设计思想：①从归一化模拟低通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型的数字滤波器。②先进行双线性变换，将模拟低通原型滤波器变换成数字低通滤波器；然后在Z域内经数字频率变换为所需类型的数字滤波器。综上所述，频率变换法是基于双线性变换法或冲激响应不变法(用于频带变换),此处选择双线性变换法主要是基于要获得严格的频率响应，以及较准确地控制截止频率的位置。而由于此种方法是一种简单的代数关系，设计也十分方1、Butterworth巴特沃斯滤波器：它具有单调下降的幅频特性；即最平幅2、Chebyshev切比雪夫滤波器：在通带或阻带等波纹，可提高选择性。3.Bessel贝塞尔滤波器：在通带内有较好的线性相位特性。4.Ellipse椭圆滤波器：其选择性相对前三种是最好的。此处选择巴特沃斯主要是想获得最平稳的幅频响应。而不计较相位特性，而用双线性变换法也会将线性相位变为非线性相位。大学《数字信号处理》报告如设计一个数字低通滤波器，其技术指标为：通带临界频率fp,通带内衰减小于rp;①将指标变为角频率wp=fp*2*pi;ws=fs*2*pi②将数字滤波器的频率指标{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)大学《数字信号处理》报告83仿真程序的设计与调试3.1数字域指标变换成模拟域指标fp=400;fs=300;Rp=1;ws=fs*2*pi;程序执行结果为：wp=2.5133e+003与实际计算结果相符。3.2数字域频率进行预畸变经过预畸变，可以发现频率变为：3.3模拟滤波器的设计其程序为%设计模拟滤波器[N,Wn]=buttord(wp2,ws2,Rp,Rs,'s')9…………之[z,p,k]=buttap(N);figure(1)freqs(Bap,Aap);%创建Buttord低通滤波器原型%由零极点转换为传递函数的形式%模拟低通滤波器的频率响应title(模拟滤波器(低通原型)的频率响应)[Bbs,Abs]=lp2hp(Bap,Aap,Wn);%模拟低通变高通figure(2)freqs(Bbs,Abs);title(模拟滤波器的频率响应)程序执行后可以发现其频率响应为：N=4,N=4,其波形如下图3-1Bi1sRditYiewInrurtToolsrTr图3-1模拟滤波器的频率响应由上图分析可得：其符合高通的一般特征，与预期的效果一样。而在此条件下，Butterworth滤波器低通原型的波形如下图3-2。EindowHalp模担滤波器的频率响应Frequency(rad/s)Usaktop四燃……………rT大学《数字信号处理》报告图3-2模拟滤波器(低通原型)的频率响应在设计的过程中，涉及一个频率变换的问题，即将模拟低通原型变为高通，其函数及用法如下：[b,a]=1p2hp(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大学《数字信号处理》报告FileEditYiewInsertToolsDe数字滤波器的频率响应0000y图3-2数字滤波器的频率响应由于使用的是双线性不变法设计的，其相位为非线性。此处主要是基于要获得严格的频率响应，以及较准确地控制截止频率的位置，故画出了详细的幅频响应。(如下图3-2)图3-3详细的幅频响应分析该图可知其在0.6(即300Hz)处的衰减为40dB,而在0.8(即400Hz)当阻带衰减变为40dB(之前为20dB),通带不变时，其波形如下图3-4。对当通带变为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.8890e+003az=[00001];bz=[12.6133.414,2.613,1];[Bbs,Abs]=lp2hp(az,bz,Wn)%用双线性不变法处理[Bbz,Abz]=bilinear(Bbs,Abs,1000);其运行结果为：N=3.6947;图形如图3-6图3-6理论计算的滤波器的幅频响应综上所述，本滤波器以四阶即实现了预期的设计目标：采样频率为1000H·,大学《数字信号处理》报告通带临界频率fp=400Hz,通带内衰减小于1dB(αp=1);阻带临界频率fs=300Hz,阻带内衰减大于20dB(αs=25),其在通带内的性能更好。4.程序调试中出现的问题1.在使用巴特沃斯函数获取其阶数时发现，调用函数时有点思路要理清，对于其函数[N,Wn]=buttord(wp2,ws2,Rp,Rs,'s’),若有后面的‘s’,N=4;而缺省时，N=6;经查资料发现，当有‘s’时，其表示的是模拟Butterworth滤波器，而缺省时是数字Butterworth滤波器。2.一个线性是不变系统的典型表达式有状态空间型、传递函数型、零极增益型、极点留数型共四种，在开始的设计过程中我发现在把零极点增益型向传递函数型的转换过程中多用了一步，即先把零极点增益型转换为状态空间型，进行相应的频率变换和双线性变换后，再由状态空间型转换为传递函数型，进而绘制出其频率响应图，但是经仿真分析后发现，性能并不理想，所以并没有采用。3.freqs(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,若没有给定，则按角频率(Angularfrequency)给定f的频率矢量；’whole’表示返回的频率f或w值包含z平面整个单位圆频率矢量，即0~2π;缺省时，频率f或w值包含z平面上半单位圆(0~π)之间点频率向量(单位rad);f为n点频率向量(Hz)。函数返回值缺省时，绘制幅频响应和相频响应图5.总结与体会软件，特别是滤波器设计中的函数基本上处于一大学《数字信号处理》报告种模糊状态。通过在学校的电子资源的期刊网上找了些论文资料，借阅图书，一点一滴的自学，以及和同学不断的交流，最后完成了这次课设，对滤波器的设计有了比较清楚的了解。在课程设计的过程中，我学到了很多东西，比如设计滤波器的一些基本函数的用法，各种模拟滤波器的特性，设计滤波器的一些基本方法。但更为重要的是，我对于解决一个问题的思路更加清晰，找到了属于自己的方法。当然，在设计的过程中，不可能避免的遇到了很多问题，如刚开始思路比较混乱，没有明确的方向。主要是如何将理论计算的模型转换为仿真模型。因为在理论上，将低通转换成高通，一般是变换将高通频率特征转换成低通原型频率特征；而在软件设计中，是对其传递函数进行修改，即变换其z域的表达式，设计初期一直不知道如何将其联系起来。后来发现，其实变换传递函数，也就是变换频率特征，是将变换后的频率代入原低通模型，而后得到高通模型的。总的来说，这次课程设计让我对MATLAB有了更深刻的了解，对数字滤波器的设计流程有了大致的了解，掌握了一些设计滤波器的基本方法，提高了理论用于实践的能力，掌握了更多专业相关的使用知识与技能。同时，也暴露了我很多的不足，在以后的学习中，将进一步发扬有点，克服缺点。大学《数字信号处理》报告参考文献在电子信息课程中的应用(第二版)·电子工业出版社·2006[2]郭仕剑·MATLAB7.X数字信号处理·人民邮电出版社·2006年[3]陈贵明·用MATLAB语言处理数字信号与数字图像·科学出版社·2000年[4]王家文·MATLAB7.0图形图像处理·国防工业出版社·2006年[5]苏金明·MATLAB图形图像·电子工业出版社·2005[6]朱铭锆·DSP应用系统设计[M]·电子工业出版社·2002.附录一总程序如下wp=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')[Bap,Aap]=zp2tf(z,p,kfreqs(Bap,Aap);%通带最大衰减Rp=1dB%阻带最小衰减Rs=20dB%把数字域滤波器特征换成模拟滤波器%预畸变求滤波器通带临界频率%预畸变求滤波器阻带临界频率%创建Buttord低通滤波器原型%由零极点转换为传递函数的形式%模拟低通滤波器的频率响应title(模拟滤波器(低通原型)的频率响应)[Bbs,Abs]=lp2hp(Bap,Aap,Wn);%模拟低通变高通freqs(Bbs,Abs);title('模拟滤波器的频率响应’)%用双线性不变法变换成数字滤波器[Bbz,Abz]=bilinear(Bbs,Abs,FS):%双线性变换%求其频率响应freqz(Bbz,Abz,512,FS);title('数字滤波器的频率响应)%详细显示数字滤波器的幅频响应[hw,w]=freqz(Bbz,Abz,512);plot(w/pi,20*log10(abs(hw)));axis([0,1,-200,10])title('ButterworthTypeHighpassDigitalFilter')%下面将笔算的结果仿真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.8890e+003az=[00001];[Bbs,Abs]=lp2hp(az,bz,Wn)%用双线性不变法处理i[Bbz,Abz]=bilinear(Bbs,Abs,1000);%画图[hw,w]=freqz(Bbz,Abz,512);plot(w/pi,20*log10(abs(hw)));axis([0,1,-200,10])set(gca,'XTickMode','manual','XTick',[0,Ws/(pi),Wp/(pi),1]);gridtitle('理论计算的滤波器的幅频响应xlabel('w/pi');ylabel('幅度(dB));附录二设计数字滤波器函数总结其中Wp和Ws分别是通带和阻带的拐角频率(截止频率),其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹cheby1函数的用法为：在使用cheby1函数设计IR滤波器之前，可使用cheblord函数求出滤波器阶数n和截止频率Wn。cheblord函数可在给定滤波器性能的情况下，选择契比雪夫I型滤波器的最小阶和截止频率Wn。其中Wp和Ws分别是通带和阻带的拐角频率(截止频率),其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹其中n代表滤波器阶数，Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。buttord函数可在给定滤波器性能的情况下，求出巴特沃斯滤波器的最小阶数n,同时给出对应的截止频率Wn。6.[bz,az]=impinvar(b,a[[,Fs[,Fp])式中，b,a为模拟滤波器分子和分母多项式系数向量；Fs为采样频率(所滤波数据),单位Hz,缺省时为1Hz。Fp为预畸变频率(Prewarpedfrequency),是一个“匹配”频率，在该频率上，频率响应在变换前后和模拟频率可精确匹配在MATLAB中，函数bilinear采用双线性变换法实现模拟s域至数字z域的映射，直接用于模拟滤波器变换为数字滤波器。其调用方式为：[zd,pd,kd]=bilinear(z,p,k,Fs)[numd,dend]=bilinear(num,den,Fs)采样频率，单位Hz。zd,pd,kd为数字滤波器的零极点和增益。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,若没有给定，则按角频率(Angularfrequency)给定f的频率矢量；’whole’表示返回的频率f或w值包含z平面整个单位圆频率矢量，即0~2π;缺省时，频率f或w值包含z平面上半单位圆(0~π)之间等间距n个点频率矢量。h为复频率响应；w为n点频率向量(单位rad);f为n点频率向量(Hz)。函数返回值缺省时，绘制幅频响应和相频响应图8.impz用于产生数字滤波器的脉冲响应。调用格式为：[[h,t]=]impz(b,a[,n,Fs])大学《数字信号处理》报告式中，b,a分别为滤波器分子和分母多项式系数向量；n为采样点数；Fs为采样频率，缺省值为1;h为滤波器单位脉冲响应向量；t为和h对应的时间向量。当函数输出缺省时，绘制滤波器脉