双线性变换法设计数字切比雪夫带通IIR滤波器.docx

1、双线性变换法设计数字切比雪夫带通IIR滤波器摘要在进行DSP系统设计时,往往先采用MATLAB软件等对算法进行仿真,确定最佳算法和 参数。利用MATLAB勺信号处理工具箱可以直接设计数字滤波器，也可以建立模拟原型,离 散化设计数字滤波器。本文介绍了 IIR数字带通滤波器的设计原理、步骤以及在MATLAB中 勺实现方法，并能应用到实际勺信号处理中。关键词：IIR 数字滤波器，数字信号处理(DSP), MATLAB仿真AbstractWhen the DSP system is designed, we often simulate the algorithm and decide the bes

2、t algorithm and parameters on the MATLAB software. We can design a digital filter either using the MATLAB Toolbox directly or creating a analog filter prototype and scatering the digital filter. This paper introduces the design principle of IIR digital bandpass filter, the steps as well as the metho

3、ds in MATLAB ,and it also can be applied to the actual signal processing.Key words: IIR digital filters, Digital signal processing (DSP), MATLAB, Simulation摘 要 IABSTRACT II1绪论 12IIR 数字滤波器设计的原理与方法 12.1 IIR 数字滤波器设计的原理 12.2 IIR 数字滤波器设计的基本方法 33IIR 带通滤波器的 MATLAB设计 33.1IIR 带通滤波器的设计流程 43.2IIR 带通滤波器的设计步骤 54

4、IIR 带通滤波器的仿真结果及波形 85IIR 带通滤波器的仿真结果分析 96总结 11参考文献 13附录 :原程序 14Matlab 课程设计 双线性变换法设计数字切比雪夫带通 IIR 滤波器1绪论在现代通信系统中，由于信号中经常混有各种复杂成分，所以很多信号的处理和分析 都是基于滤波器而进行的。数字滤波器可以理解为是一个计算程序或算法，将代表输入信 号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定 的形式变化。数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤 波器分为两种，即无限长冲激响应 (IIR) 滤波器和有限长冲激响应 (FIR) 滤波

5、器。但是，传统 的数字滤波器的设计使用繁琐的公式计算，改变参数后需要重新计算，从而在设计滤波器 尤其是高阶滤波器时工作量很大。利用 MATLAB信号处理箱(Signal Processing Toolbox可 以快速有效地实现数字滤波器的设计与仿真。本文设计一个 IIR 数字带通滤波器。 IIR 数字滤波器具有无限宽的冲激响应，与模拟滤 波器相匹配，所以 IIR 滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方 法。其设计方法主要有经典设计法、直接设计法和最大平滑滤波器设计法。在对滤波器实 际设计时，整个过程的运算量是很大的。设计完成后对已设计的滤波器的频率响应要进行校核，要得到幅频

6、、相频响应特性， 运算量也是很大的。平时所要设计的数字滤波器，阶数和类型并不一定是完全给定的，很 多时候要根据设计要求和滤波效果不断地调整，以达到设计的最优化。在这种情况下，滤 波器设计就要进行大量复杂的运算，单纯的靠公式计算和编制简单的程序很难在短时间内 完成。利用 MATLAB 强大的计算功能进行计算机辅助设计，可以快速有效地设计数字滤 波器，大大地简化了计算量。2IIR 数字滤波器设计的原理与方法2.1IIR 数字滤波器设计的原理IIR 数字滤波器具有无限持续时间冲激响应 ,需要用递归模型来实现 , 其系统函数为 :(2.1)设计IIR滤波器的任务就是寻求一个因果、物理上可实现的系统函数

7、 H (z)，使其频率响应H(ejw)满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止、通带 衰减和阻带衰减.0利用冲激响应不变法设计数字滤波器时可能会导致频域混叠现象，为了克服这一问 题，需要找到由s平面到z平面的另外的映射关系，这种关系应保证：1)s平面的整个j Q轴仅映射为z平面单位圆上的一周；2)若G(s)是稳定的，由G(s)映射得到的H(z)也应该是稳定的;3)这种映射是可逆的，既能由 G(s)得到H(z)，也能由H(z)得到G(s);4)如果 G(j0)=1，那么 H(ej0) = 1o双线性Z变换满足以上4个条件的映射关系，其变换公式为2 z-1(2.2)n /Ts变

8、换到2n /Ts的s 二TsZ 1双线性Z变换的基本思路：首先将整个S平面压缩到一条从-横带里，然后通过标准的变换关系 z=esT将横带变换到整个Z平面上去，这样就得到了S 平面与Z平面间的对应的单值关系。图2.1双线性变换法 S平面到Z平面的映射关系2.2IIR 数字滤波器设计的基本方法IIR数字滤波器的设计方法有两类，一类是借助于模拟滤波器的设计方法设计出模拟滤波器,利用冲激响应不变法或双线性变换法转换成数字滤波器 ，再用硬件或软件实现；另一类是直接在频域或时域中进行设计,设计时需要计算机作辅助工具。随着MATLAB软件尤 其是MATLAB的信号处理工作箱的不断完善,不仅数字滤波器的计算机

9、辅助设计有了可能, 而且还可以使设计达到最优化。IIR数字滤波器设计的基本步骤如下：(1)根据任务，确定性能指标。在设计一个滤波器之前，首先根据工程实际的需要确定滤波器 的技术指标如：边界频率：3 p , 3 s , 3 c ;阻带最小衰减As和通带最大衰减Rp；(2)将数字滤波器的技术指标转换成模拟滤波器指标。利用冲激响应不变法与双线性变换 法进行频率间的转换，主要是边界频率 Wp与Ws的转换。(3)用模拟滤波器设计方法得到模拟滤波器的传输函数 Ha ( s);可借助巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer)滤波器、贝塞尔(Bessel

10、) 滤波器等，这些滤波器都有严格的设计公式、现成的曲线和图表供设计人员使用。 映射实现。利用双线性变换法将模拟滤波器 Ha ( s)转换成数字滤波器H( z )。(5)用有限精度算法实现这个系统函数 H ( z )(包括选择运算结构、选择合适的字长、有效 数字处理方法)。(6)用适当的软、硬件技术实现。包括采用通用计算机软件、数字滤波器硬件或者软硬件 结合，确定DF采用的结构将会影响其精度、稳定性、经济性及运算速度等很多重要性质。3IIR带通滤波器的MATLAB设计IIR带通滤波器的设计框图如下：图3.1 IIR 带通滤波器的设计框图3.1 IIR 带通滤波器的设计流程确定带通滤波器的技 术指

11、标设计切比雪夫低通滤波器并求出相关参数双线性变换：模拟低 通滤波器转化为带通 滤波器）求相应的幅频响应与相频响应图3.2 IIR 带通滤波器的设计流程本文设计的IIR带通滤波器是从低通变换过来的，利用的是双线性变换以及切比雪夫II滤波器的原型，其具体的设计流程为上图所示。首先根据题目要求确定带通滤波器的技 术指标，先要进行频率的预畸变，并且归一化频率，再设计出切比雪夫II模拟低通滤波器， 并求出其阶数等相关参数。其次利用双线性变换法设计数字带通滤波器， ，再调用函数进行双线性变换，并求出分子、分母的系数向量。最后通过画图求出其幅频响应、相频响应、 幅度特性曲线与零极点，并画出波形图。最后进行验

12、证，看所设计的滤波器能否达到要求 的指标，若能达到，则说明该滤波器设计符合要求。3.2IIR带通滤波器的设计步骤(1)根据设计流程，首先确定所要设计的数字带通滤波器的相关指标：1通带截止频率 wp1 =0.4 n，wp2 =0.6 n，通带最大衰减 Rp = 2dB;2阻带截止频率 ws1=0.2n 和ws2=0.8n，阻带最小衰减 Rs = 30dB;3取样间隔T = 0.1s。其实现程序如下(程序中pi代表n ):Ts = 0.1; Fs=1/Ts;%取样周期或频率Rp = 2;%通带最大衰减Rs = 30;%阻带最小衰减wp1=0.4*pi;%通带、阻带上、下限截止频率wp2=0.6*p

13、i;ws1=0.2*pi;ws2=0.8*pi;(2)频率的预畸变。双线性变换中无法避免的一个问题即是频率的非线性偏移，因为数字频率的最大值为 n，而模拟频率可以向无穷延伸，两者之间又要保持一一对应的映射关系。双线性变换中 的模拟角频率门与数字角频率之间的关系为：小 2 蛍tan () (3.1)T 2表明S平面与Z平面是单值的一一对应关系，即频率轴是单值变换关系。虽然避免了 脉冲响应不变法的频率响应的混叠现象，但是经过变换后，得到的幅频响应特性各分段边 缘频率不能保持原来的比例关系，必须通过预修正加以校正。做法是将数字频率 按门二2/T*tan(w/2)的关系，变成模拟频率门，利用这组做过修

14、正的模拟频率来设计模拟带通滤波 器作为模拟原型。Wp1= (2/T)*ta n(wp1/2);Wp2=(2/T)*ta n(wp2/2);Wp=Wp1,Wp2; %模拟滤波器的通带截止频率Ws仁(2/T)*ta n(ws1/2);Ws2=(2/T)*tan(ws2/2);模拟滤波器的阻带截止频率模拟滤波器的带宽模拟滤波器的中心频率Ws=Ws1,Ws2;BW=Ws2-Ws1;Omegaw0=sqrt(Ws1*Ws2);(2)设计切比雪夫模拟低通滤波器。求模拟低通滤波器的阶数与边缘频率N,OmegaC=cheb2ord(Wp,Ws,Rp,Rs,s);求切比雪夫 II 型模拟低通滤波器的零、极点与增

15、益z0,p0,k0=cheb2ap(N,Rs);利用函数 N,OmegaC=cheb2ord(Wp,Ws,Rp,Rs, s) ，通过给定滤波器的技术指标Wp Ws Rp、Rs，求得滤波器的阶数N与边缘频率OmegaC。Wp Ws、与OmegaC均 在0,1 区间归一化，以n弧度为单位。利用函数z,p,k=cheb2ap(N,Rs) ，来设计一个阶数为N，阻带波动为Rs的归一化切比雪夫II型原型滤波器，得到左半平面零极点。数组 Z中返回零点，数组P中返回极点， 并且返回增益 K。(3)设计归一化的模拟原型带通滤波器：%求原型滤波器的分子系数AnalogB=k0*real(poly(z0);%求原

16、型滤波器的分母系数AnalogA=real(poly(p0);%模拟低通到模拟带通的分子、分母系数的变换BandB,BandA=lp2bp(AnalogB,AnalogA,Omegaw0,BW);%双线性变换：模拟带通与数字带通的分子分母系数的变换DigitalB,DigitalA = bilinear(BandB,BandA,Fs);%变为二阶节级联结构sos,G = tf2sos(DigitalB,DigitalA);利用函数p = poly(A)来计算模拟滤波器的分子、分母系数向量，因其为实数，因此用 real()函数取其实部，即可得模拟滤波器的分子、分母系数向量。这两个函数实现的功能可

17、 以用函数Bs,As=zp2tf(z,p,k)直接求得传递函数的分子、分母系数向量。利用函数BandB,BandA=lp2bp(AnalogB,AnalogA,OmegawO,BW)，将模拟域的低通变为带通，并且 得到模拟带通滤波器的分子、分母系数向量， OmegawO取为中心频率，BW为带宽。利用函数DigitalBQigitalA = bilinear(BandB,BandA,Fs)，双线性变换为数字带通滤波器的指标， 如分子、分母的系数向量。函数sos,G = tf2sos(DigitalB,DigitalA)，即把z变换传递函数的 直接形式转换成级联形式。需要注意的是，这个函数是针对以

18、z的负幕排列的多项式开发的， 虽然可以推广到s域，但连续系统传递函数是按s的正幕排列的，要使两者一致，关键是使 分子、分母系数向量同长，两序列中各元素的幕次排列一致(4)求数字带通滤波器的幅频、相频特性、及其群延迟%求数字带通滤波器的幅频特性H z, Wz=freqz(DigitalB,DigitalA,1024,whole);%将数字带通滤波器的幅频特性转化为分贝表示dbHz=20*log10(abs(Hz)+eps)/max(abs(Hz);%求数字带通滤波器的相频特性 = an gle(Hz)%求数字带通滤波器的群延迟特性grd = grpdelay(DigitalB,DigitalA,

19、Wz);函数Hz,Wz=freqz(DigitalB,DigitalA,1024,whole)可以求数字带通滤波器的幅频特性，而其幅度(即模值)的最大值可以归一化为1,则其模值(单位为dB)即可以用公式表示 为dbHz=20*log10(abs(Hz)+eps)/max(abs(Hz)。函数 = angle(Hz)可求得其相频特性，而 对于一个滤波器来说，要满足其线性相位，则其群延迟 要为一常数，则其相位特性必须为一直线，即满足，-t。因此，利用群延迟函数 grd = grpdelay(DigitalB,DigitalA,Wz)可 以判断所设计的滤波器是否是线性相位，如果不符合，可以更改参数加

20、以较正或者用其他 方法重新设计，从而方便了设计。(5)作图程序见附录。以上是根据切比雪夫II型滤波器设计的带通滤波器，如果以切比雪夫I型滤波器设计， 则只需改变几个参数则可，其程序如下所示(只附改变的参数与程序) 。%模拟滤波器的带宽BW=Wp2-Wp1;%模拟滤波器的中心频率OmegawO=sqrt(Wp1*Wp2);%求切比雪夫I型模拟低通滤波器的阶数与边缘频率N,OmegaC=cheb1ord(Wp,Ws,Rp,Rs,s);%求切比雪夫I型模拟低通滤波器的零、极点与增益zO,pO,kO=cheb1ap(N,Rp);4IIR带通滤波器的仿真结果及波形5根据上述对IIR切比雪夫I、II带通滤

21、波器的设计过程，在 Matlab软件中得到其幅频 响应、相频响应、零极点图以及群延迟特性曲线如下图所示。图4.1 Chebyshev I 型IIR带通滤波器特性图4.2 Chebyshev II 型IIR数字带通滤波器特性6IIR带通滤波器的仿真结果分析根据上述的仿真波形，可以看出：（1） 由ChebyshevI型、ChebyshevII型低通原型变换成带通模型的幅频特性可以用分贝形式表示。对于ChebyshevI型的幅频响应在通带内为等纹波衰减、阻带内为单调减小的， 且通带（归一化）截止频率在波形上显示为： 0.4与0.6,阻带（归一化）截止频率在波形上显示为：0.2与0.8,与设计要求基本

22、一致。对于 ChebyshevII型的幅频响应来说，在通 带内单调减小、阻带内为等纹波衰减的，而通带、阻带（归一化）截止频率为： 0.4与0.6,0.2与0.8，均满足设计要求。（2） 此次设计的IIR带通滤波器的阶数可由函数N,OmegaC=cheb1ord（Wp,Ws,Rp,Rs,s）或N,OmegaC=cheb2ord（Wp,Ws,Rp,Rs,s求得，N=3。则从 ChebyshevI 型带通滤波器的幅 频响应曲线上看，当N = 3时，在0 ：二（归一化后为0-1）范围内通带波动只有两个波谷， 阻带则单调减小，这与理论内容一致。对于 ChebyshevII型带通滤波器来说，其幅频响应 曲

23、线的通带波动在0 ：二（归一化后为0-1 ）范围内阻带波动只有两个波峰，通带则单调减 小，这与理论内容一致。因此，对于滤波器阶数这一指标来说，满足要求。（3） 对于分贝化后的幅度特性，Chebyshevl型带通滤波器的通带最大衰减大约为2dB,阻 带衰减最小为30dB，而Chebyshevll型带通滤波器的阻带最小衰减大约为 30dB，通带最小 衰减大约为2dB，这与带通滤波器的技术指标相同。因此，在衰减性能上，所设计的滤波器达到了要求。（4） 对于相频响应来说，Chebyshevl型带通滤波器在通带内的相频响应曲线接近为一条平滑曲线，在阻带内则存在衰减、畸变。 Chebyshevll型带通滤

24、波器在通带截止频率以内的相频响应曲线几乎为一条直线，但在通带截止频率处有很大角度的转折，使相频特性在 通带内产生了的畸变，而在阻带内存在很大的畸变。因此，对于双线性变换法，其相位特 性得不到满足，必须用其它方法加以校正或采用其它方法重新设计滤波器，使之满足线性 相位。（5） 群延迟是衡量一个滤波器或是整个系统性能指标的一个重要参数。对于利用这两种低通滤波器原型模型设计的这两种带通滤波器，带通滤波器的延迟是同等带宽的低通滤波 器延迟的两倍。这个结果是由低通特性向带通特性变换引起的， n阶低通滤波器传递函数总是变换为2n阶带通滤波器传递函数，其滤波器的群延迟特性的中心频率与理论上的很接 近。（6）

25、 零、极点图则反映Z域的传递函数，因为从图上得知， Chebyshevl型带通滤波器在 Z域单位圆内的传递函数存在3对极点，其每一对则关于实轴对称，而且存在 1对零点，且 关于虚轴对称。Chebyshevll型带通滤波器在Z域单位圆内的传递函数存在 3对极点，其每 一对则关于实轴对称，而且存在 3对零点，且关于虚轴对称。从理论上分析，可以求出其 传递函数。综上所述，由于软件的精度要求以及所调用函数的特性，在误差允许范围内，上述的 滤波器特性的仿真结果基本符合题目的要求，并且在一定的程度上设计的指标远远好于题 目的要求。因此，所设计的滤波器达到了题目的要求。7总结此次 Matlab 课程设计是在

26、一定的理论基础之上进行的，在先修课程信号与系统 与 数字信号处理中，大量有 Matlab 设计方面的知识，所以做起来还比较容易，而且经过 了自己的亲身实践，学到了许多实践方面的知识。首先，在信号滤波系统中，有时因为模拟滤波器阶数太高，硬件占用空间太大为某些 仪器的实现设置了障碍，而对于一些窄带情况下的低通滤波器用模拟手段往往很难实现。 在这些情况下， 数字滤波器将会是一个很好的解决办法。 MATLAB 信号处理工具箱提供了 丰富而简便的设计、实现 FIR 和 IIR 的方法，使原来繁琐的程序设计简化成函数的调用， 特别是滤波器的表达方式和滤波器之间的相互转换显得十分简便。其次， IIR 数字滤

27、波器的设计和模拟滤波器的设计有着紧密的关系。通常要先设计出 适当的模拟滤波器，再通过一定的频带变换把它转换成为所需的数字 IIR 滤波器。此外， 任何数字信号处理系统中也还不可避免地用到模拟滤波器， 比如 A/D 变换器前的抗混叠滤 波器和 D/A 变换器后的平滑滤波器，因此模拟滤波器设计也是很重要的。最后，在比较设计滤波器的方法上应该明确其技术指标以及某些参数的实际意义。比 如本文用双线性变换法设计数字带通滤波器时，必须先将频率归一化，并且进行频率预畸 变，然后设计模拟滤波器，再利用频率变换法将模拟低通变为模拟带通，最后经过双线性 变换法将模拟带通变换为数字带通。如果不进行频率预畸变，那么设

28、计出来的带通滤波器 的幅频特性与相频特性将会产生很严重的畸变（如图 5、图 6 所示，为设计时没有进行频 率预畸变的特性），使设计的结果不满足给定的要求，在实际中会造成很严重的危害。在做本次课程设计的过程中，我深深地感受到了自己所学到知识的有限，明白了只学 好课本上的知识是不够的，要通过图书馆和互联网等各种渠道来扩充自己的知识。在实验 过程中我们曾经遇到过问题。但是从中我们学习到了如何对待遇到的困难，进一步培养了 我们一丝不苟的科学态度和不厌其烦的耐心。所有的这些心得会对我以后的学习和工作有 帮助作用，忠心感谢学校给我们提供这次实验机会。图6.2 Chebyshev II 型IIR 数字带通滤

29、波器特性（无频率预畸变）参考文献1数字信号处理教程 MATLAB 释义与实现（第二版） .陈怀琛编著 .北京 .电子工业 出版社 ,2008.102MATLAB 及在理工课程中的应用指南 .陈怀琛编著 .西安 .电子科技大学出版社 ,20003数字信号处理原理与实现 .刘泉 ,阙大顺编著 .北京 .电子工业出版社 ,2005.84数字信号处理教程 （第二版） .北京.清华大学出版社 ,19975数字信号处理及其 MATLAB实现.Vi nay K.l ngle编著电子工业出版社，19986 Digital Signal Processing Laboratory Using Matlab. S

30、anjit K. Miltra 编著. McGraw-Hill 出版社，2000附录：原程序所设计的数字滤波器的指标Ts = 0.1，Fs=1/Ts， Rp = 2， Rs = 30 wp1=0.4*pi ，wp2=0.6*pi; ws1=0.2*pi ， ws2=0.8*pi;频率的预畸变Wp1=(2/T)*tan(wp1/2);Wp2=(2/T)*tan(wp2/2);模拟滤波器的通带截止频率模拟滤波器的阻带截止频率模拟滤波器的带宽模拟滤波器的中心频率Wp=Wp1,Wp2;Ws1=(2/T)*tan(ws1/2);Ws2=(2/T)*tan(ws2/2);Ws=Ws1,Ws2;BW=Ws2

31、-Ws1;% BW=Wp2-Wp1;Omegaw0=sqrt(Ws1*Ws2);% Omegaw0=sqrt(Wp1*Wp2);求模拟低通滤波器的阶数与边缘频率N,OmegaC=cheb2ord(Wp,Ws,Rp,Rs,s); %N,OmegaC=cheb1ord(Wp,Ws,Rp,Rs,s) 求切比雪夫模拟低通滤波器的零、极点与增益z0,p0,k0=cheb2ap(N,Rs);%z0,p0,k0=cheb1ap(N,Rp); %设计归一化的模拟原型带通滤波器 %求原型滤波器的分子系数 AnalogB=k0*real(poly(z0);%求原型滤波器的分母系数AnalogA=real(poly(p0);%模拟低通到模拟带通的分子、分母系数的变换 BandB,BandA=lp2bp(AnalogB,AnalogA,Omegaw0,BW); %双线性变换：模拟带通与数字带通的分子