毕业论文：双线性变换法设计数字切比雪夫带通IIR滤波器.doc

《毕业论文：双线性变换法设计数字切比雪夫带通IIR滤波器.doc》由会员分享，可在线阅读，更多相关《毕业论文：双线性变换法设计数字切比雪夫带通IIR滤波器.doc（19页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

2、ameters 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 methods in MATLAB ,and i

3、t also can be applied to the actual signal processing. Key words: IIR digital filters, Digital signal processing (DSP), MATLAB, Simulation目 录摘要IABSTRACTII1绪论12 IIR数字滤波器设计的原理与方法12.1 IIR数字滤波器设计的原理12.2 IIR 数字滤波器设计的根本方法33 IIR带通滤波器的MATLAB 设计33.1 IIR带通滤波器的设计流程43.2 IIR带通滤波器的设计步骤54 IIR带通滤波器的仿真结果及波形85 IIR带通滤

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

5、滤波器尤其是高阶滤波器时工作量很大。利用MATLAB信号处理箱(Signal Processing Toolbox)可以快速有效地实现数字滤波器的设计与仿真。 本文设计一个IIR数字带通滤波器。IIR数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配，所以IIR滤波器的设计可以采取在模拟滤波器设计的根底上进一步变换的方法。其设计方法主要有经典设计法、直接设计法和最大平滑滤波器设计法。在对滤波器实际设计时，整个过程的运算量是很大的。 设计完成后对已设计的滤波器的频率响应要进行校核，要得到幅频、相频响应特性，运算量也是很大的。平时所要设计的数字滤波器，阶数和类型并不一定是完全给定的，很多时候要根据

6、设计要求和滤波效果不断地调整，以到达设计的最优化。在这种情况下，滤波器设计就要进行大量复杂的运算，单纯的靠公式计算和编制简单的程序很难在短时间内完成。利用MATLAB强大的计算功能进行计算机辅助设计，可以快速有效地设计数字滤波器，大大地简化了计算量。2 IIR数字滤波器设计的原理与方法2.1 IIR数字滤波器设计的原理IIR 数字滤波器具有无限持续时间冲激响应,需要用递归模型来实现,其系统函数为: (2.1)设计IIR 滤波器的任务就是寻求一个因果、物理上可实现的系统函数H (z) ，使其频率响应满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止、通带衰减和阻带衰减.。利用冲激响应不

7、变法设计数字滤波器时可能会导致频域混叠现象，为了克服这一问题，需要找到由s平面到z平面的另外的映射关系，这种关系应保证：1) s平面的整个j轴仅映射为z平面单位圆上的一周；2) 假设G(s)是稳定的，由G(s)映射得到的H(z)也应该是稳定的；3) 这种映射是可逆的，既能由G(s)得到H(z)，也能由H(z)得到G(s)；4)如果G(j0)=1，那么1。双线性变换满足以上4个条件的映射关系，其变换公式为 2.2双线性变换的根本思路：首先将整个平面压缩到一条从/Ts变换到2/Ts 的横带里，然后通过标准的变换关系将横带变换到整个平面上去，这样就得到了平面与平面间的一一对应的单值关系。Rezj s

8、平面j1s平面jIm(t)z平面图2.1 双线性变换法S 平面到Z 平面的映射关系2.2 IIR 数字滤波器设计的根本方法IIR 数字滤波器的设计方法有两类,一类是借助于模拟滤波器的设计方法设计出模拟滤波器,利用冲激响应不变法或双线性变换法转换成数字滤波器,再用硬件或软件实现；另一类是直接在频域或时域中进行设计,设计时需要计算机作辅助工具。随着MATLAB 软件尤其是MATLAB 的信号处理工作箱的不断完善,不仅数字滤波器的计算机辅助设计有了可能,而且还可以使设计到达最优化。IIR 数字滤波器设计的根本步骤如下:(1) 根据任务,确定性能指标。在设计一个滤波器之前,首先根据工程实际的需要确定滤

9、波器的技术指标如：边界频率：p ，s ，c ；阻带最小衰减As 和通带最大衰减Rp；(2) 将数字滤波器的技术指标转换成模拟滤波器指标。利用冲激响应不变法与双线性变换法进行频率间的转换，主要是边界频率Wp与Ws 的转换。(3) 用模拟滤波器设计方法得到模拟滤波器的传输函数Ha ( s) ；可借助巴特沃斯(Butterworth) 滤波器、切比雪夫(Chebyshev) 滤波器、椭圆(Cauer) 滤波器、贝塞尔(Bessel) 滤波器等，这些滤波器都有严格的设计公式、现成的曲线和图表供设计人员使用。(4) 映射实现。利用双线性变换法将模拟滤波器Ha ( s) 转换成数字滤波器H( z ) 。(

10、5) 用有限精度算法实现这个系统函数H ( z ) (包括选择运算结构、选择适宜的字长、有效数字处理方法 )。(6) 用适当的软、硬件技术实现。包括采用通用计算机软件、数字滤波器硬件或者软硬件结合，确定DF采用的结构将会影响其精度、稳定性、经济性及运算速度等很多重要性质。3 IIR带通滤波器的MATLAB 设计模拟低通滤波器的设计模拟频带变换：低通变高通、带通等滤波器变换：模拟变为数字IIR滤波器IIR带通滤波器的设计框图如下：图 3.1 IIR带通滤波器的设计框图3.1 IIR带通滤波器的设计流程开始确定带通滤波器的技术指标设计切比雪夫低通滤波器并求出相关参数求相应的幅频响应与相频响应双线性

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

12、3.2 IIR带通滤波器的设计步骤1根据设计流程，首先确定所要设计的数字带通滤波器的相关指标：通带截止频率wp1 =0.4，wp2 =0.6，通带最大衰减Rp = 2dB；阻带截止频率ws1=0.2 和ws2=0.8，阻带最小衰减Rs = 30dB；取样间隔T = 0.1s。其实现程序如下程序中pi代表：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;2频率的预畸变。双线性变换中无法防止的一个问题即是频率的非线性偏移，因为数字

13、频率的最大值为，而模拟频率可以向无穷延伸，两者之间又要保持一一对应的映射关系。双线性变换中的模拟角频率与数字角频率之间的关系为： (3.1)说明S平面与Z平面是单值的一一对应关系，即频率轴是单值变换关系。虽然防止了脉冲响应不变法的频率响应的混叠现象，但是经过变换后，得到的幅频响应特性各分段边缘频率不能保持原来的比例关系，必须通过预修正加以校正。做法是将数字频率按2/T*tan(w/2)的关系，变成模拟频率，利用这组做过修正的模拟频率来设计模拟带通滤波器作为模拟原型。Wp1=(2/T)*tan(wp1/2);Wp2=(2/T)*tan(wp2/2);Wp=Wp1,Wp2; 模拟滤波器的通带截止频

14、率 Ws1=(2/T)*tan(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型模拟低通滤波器的零、极点与增益z0,p0,k0=cheb2ap(N,Rs); 利用函数N,OmegaC=cheb2ord(Wp,Ws,Rp,Rs,s)，通过给定滤波器的技术指标Wp、Ws、Rp、Rs，求得滤波器

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

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

17、通变为带通，并且得到模拟带通滤波器的分子、分母系数向量，Omegaw0取为中心频率，BW为带宽。利用函数DigitalB,DigitalA = bilinear(BandB,BandA,Fs)，双线性变换为数字带通滤波器的指标，如分子、分母的系数向量。函数sos,G = tf2sos(DigitalB,DigitalA)，即把z变换传递函数的直接形式转换成级联形式。需要注意的是，这个函数是针对以z的负幂排列的多项式开发的，虽然可以推广到s域，但连续系统传递函数是按s的正幂排列的，要使两者一致，关键是使分子、分母系数向量同长，两序列中各元素的幂次排列一致4求数字带通滤波器的幅频、相频特性、及其群

18、延迟求数字带通滤波器的幅频特性Hz,Wz=freqz(DigitalB,DigitalA,1024,whole);将数字带通滤波器的幅频特性转化为分贝表示dbHz=20*log10(abs(Hz)+eps)/max(abs(Hz); 求数字带通滤波器的相频特性angle(Hz)求数字带通滤波器的群延迟特性grd = grpdelay(DigitalB,DigitalA,Wz); 函数Hz,Wz=freqz(DigitalB,DigitalA,1024,whole)可以求数字带通滤波器的幅频特性，而其幅度即模值的最大值可以归一化为1，那么其模值单位为dB即可以用公式表示为dbHz=20*log1

19、0(abs(Hz)+eps)/max(abs(Hz)。函数angle(Hz)可求得其相频特性，而对于一个滤波器来说，要满足其线性相位，那么其群延迟要为一常数，那么其相位特性必须为一直线，即满足。因此，利用群延迟函数grd = grpdelay(DigitalB,DigitalA,Wz)可以判断所设计的滤波器是否是线性相位，如果不符合，可以更改参数加以较正或者用其他方法重新设计，从而方便了设计。5作图程序见附录。以上是根据切比雪夫II型滤波器设计的带通滤波器，如果以切比雪夫I型滤波器设计，那么只需改变几个参数那么可，其程序如下所示只附改变的参数与程序。模拟滤波器的带宽BW=Wp2-Wp1; 模拟

20、滤波器的中心频率 Omegaw0=sqrt(Wp1*Wp2);求切比雪夫型模拟低通滤波器的阶数与边缘频率N,OmegaC=cheb1ord(Wp,Ws,Rp,Rs,s); 求切比雪夫I型模拟低通滤波器的零、极点与增益 z0,p0,k0=cheb1ap(N,Rp); 4 IIR带通滤波器的仿真结果及波形根据上述对IIR切比雪夫I、II带通滤波器的设计过程，在Matlab软件中得到其幅频响应、相频响应、零极点图以及群延迟特性曲线如下列图所示。图4.1 Chebyshev I 型IIR带通滤波器特性图4.2 Chebyshev II 型IIR数字带通滤波器特性5 IIR带通滤波器的仿真结果分析根据上

21、述的仿真波形，可以看出：1由ChebyshevI 型、ChebyshevII型低通原型变换成带通模型的幅频特性可以用分贝形式表示。对于ChebyshevI 型的幅频响应在通带内为等纹波衰减、阻带内为单调减小的，且通带归一化截止频率在波形上显示为：0.4与0.6，阻带归一化截止频率在波形上显示为：0.2与0.8，与设计要求根本一致。对于ChebyshevII型的幅频响应来说，在通带内单调减小、阻带内为等纹波衰减的，而通带、阻带归一化截止频率为：0.4与0.6，0.2与0.8，均满足设计要求。2此次设计的IIR带通滤波器的阶数可由函数N,OmegaC=cheb1ord(Wp,Ws,Rp,Rs,s)

22、或N,OmegaC=cheb2ord(Wp,Ws,Rp,Rs,s)求得，N=3。那么从ChebyshevI 型带通滤波器的幅频响应曲线上看，当N3时，在归一化后为0-1范围内通带波动只有两个波谷，阻带那么单调减小，这与理论内容一致。对于ChebyshevII 型带通滤波器来说，其幅频响应曲线的通带波动在归一化后为0-1范围内阻带波动只有两个波峰，通带那么单调减小，这与理论内容一致。因此，对于滤波器阶数这一指标来说，满足要求。3对于分贝化后的幅度特性，ChebyshevI 型带通滤波器的通带最大衰减大约为2dB，阻带衰减最小为30dB，而ChebyshevII 型带通滤波器的阻带最小衰减大约为3

23、0dB，通带最小衰减大约为2dB，这与带通滤波器的技术指标相同。因此，在衰减性能上，所设计的滤波器到达了要求。4对于相频响应来说，ChebyshevI 型带通滤波器在通带内的相频响应曲线接近为一条平滑曲线，在阻带内那么存在衰减、畸变。ChebyshevII 型带通滤波器在通带截止频率以内的相频响应曲线几乎为一条直线，但在通带截止频率处有很大角度的转折，使相频特性在通带内产生了的畸变，而在阻带内存在很大的畸变。因此，对于双线性变换法，其相位特性得不到满足，必须用其它方法加以校正或采用其它方法重新设计滤波器，使之满足线性相位。5群延迟是衡量一个滤波器或是整个系统性能指标的一个重要参数。对于利用这两

24、种低通滤波器原型模型设计的这两种带通滤波器，带通滤波器的延迟是同等带宽的低通滤波器延迟的两倍。这个结果是由低通特性向带通特性变换引起的，n阶低通滤波器传递函数总是变换为2n阶带通滤波器传递函数，其滤波器的群延迟特性的中心频率与理论上的很接近。6零、极点图那么反映Z域的传递函数，因为从图上得知，ChebyshevI 型带通滤波器在域单位圆内的传递函数存在3对极点，其每一对那么关于实轴对称，而且存在1对零点，且关于虚轴对称。ChebyshevII 型带通滤波器在域单位圆内的传递函数存在3对极点，其每一对那么关于实轴对称，而且存在3对零点，且关于虚轴对称。从理论上分析，可以求出其传递函数。综上所述，

25、由于软件的精度要求以及所调用函数的特性，在误差允许范围内，上述的滤波器特性的仿真结果根本符合题目的要求，并且在一定的程度上设计的指标远远好于题目的要求。因此，所设计的滤波器到达了题目的要求。6 总结此次Matlab课程设计是在一定的理论根底之上进行的，在先修课程?信号与系统?与?数字信号处理?中，大量有Matlab设计方面的知识，所以做起来还比拟容易，而且经过了自己的亲身实践，学到了许多实践方面的知识。首先，在信号滤波系统中，有时因为模拟滤波器阶数太高，硬件占用空间太大为某些仪器的实现设置了障碍，而对于一些窄带情况下的低通滤波器用模拟手段往往很难实现。在这些情况下，数字滤波器将会是一个很好的解

26、决方法。MATLAB信号处理工具箱提供了丰富而简便的设计、实现FIR 和IIR 的方法，使原来繁琐的程序设计简化成函数的调用，特别是滤波器的表达方式和滤波器之间的相互转换显得十分简便。其次，IIR数字滤波器的设计和模拟滤波器的设计有着紧密的关系。通常要先设计出适当的模拟滤波器，再通过一定的频带变换把它转换成为所需的数字IIR滤波器。此外，任何数字信号处理系统中也还不可防止地用到模拟滤波器，比方A/D变换器前的抗混叠滤波器和D/A变换器后的平滑滤波器，因此模拟滤波器设计也是很重要的。最后，在比拟设计滤波器的方法上应该明确其技术指标以及某些参数的实际意义。比方本文用双线性变换法设计数字带通滤波器时

27、，必须先将频率归一化，并且进行频率预畸变，然后设计模拟滤波器，再利用频率变换法将模拟低通变为模拟带通，最后经过双线性变换法将模拟带通变换为数字带通。如果不进行频率预畸变，那么设计出来的带通滤波器的幅频特性与相频特性将会产生很严重的畸变如图5、图6所示，为设计时没有进行频率预畸变的特性，使设计的结果不满足给定的要求，在实际中会造成很严重的危害。在做本次课程设计的过程中，我深深地感受到了自己所学到知识的有限，明白了只学好课本上的知识是不够的，要通过图书馆和互联网等各种渠道来扩充自己的知识。在实验过程中我们曾经遇到过问题。但是从中我们学习到了如何对待遇到的困难，进一步培养了我们一丝不苟的科学态度和不

28、厌其烦的耐心。所有的这些心得会对我以后的学习和工作有帮助作用，忠心感谢学校给我们提供这次实验时机。图 6.1Chebyshev I 型IIR数字带通滤波器特性无频率预畸变图 6.2 Chebyshev II 型 IIR 数字带通滤波器特性无频率预畸变参考文献1 数字信号处理教程MATLAB释义与实现第二版.陈怀琛编著.北京.电子工业出版社,2021.102 MATLAB及在理工课程中的应用指南.陈怀琛编著.西安.电子科技大学出版社,20003 数字信号处理原理与实现.刘泉,阙大顺编著.北京.电子工业出版社,2005.84 数字信号处理教程(第二版) .北京.清华大学出版社,19975 数字信号

29、处理及其MATLAB实现. Vinay K.Ingle编著.电子工业出版社,19986 Digital Signal Processing Laboratory Using Matlab. Sanjit 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; 模拟滤波器的通带截止频

30、率 Ws1=(2/T)*tan(ws1/2);Ws2=(2/T)*tan(ws2/2);Ws=Ws1,Ws2; 模拟滤波器的阻带截止频率 BW=Ws2-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=cheb

31、1ap(N,Rp); %设计归一化的模拟原型带通滤波器%求原型滤波器的分子系数AnalogB=k0*real(poly(z0); %求原型滤波器的分母系数AnalogA=real(poly(p0); %模拟低通到模拟带通的分子、分母系数的变换BandB,BandA=lp2bp(AnalogB,AnalogA,Omegaw0,BW); %双线性变换：模拟带通与数字带通的分子分母系数的变换 DigitalB,DigitalA = bilinear(BandB,BandA,Fs); %变为二阶节级联结构sos,G = tf2sos(DigitalB,DigitalA); %求数字带通滤波器的幅频、相

32、频特性、及其群延迟求数字带通滤波器的幅频特性Hz,Wz=freqz(DigitalB,DigitalA,1024,whole);将数字带通滤波器的幅频特性转化为分贝表示dbHz=20*log10(abs(Hz)+eps)/max(abs(Hz); 求数字带通滤波器的相频特性angle(Hz)求数字带通滤波器的群延迟特性grd = grpdelay(DigitalB,DigitalA,Wz); 作图subplot(2,3,1);plot(Wz/pi,abs(Hz);title(幅频响应);xlabel(),ylabel(幅度:|Hz|);axis(0,1,0,1.1);set(gca,XTick

33、Mode,manual,XTick,0,ws1/pi,wp1/pi,wp2/pi,ws2/pi,1);grid;subplot(2,3,4);plot(Wz/pi,dbHz);title(模值(dB);xlabel(频率单位：pi);ylabel(分贝dB);axis(0,1,-40,5);set(gca,XTickMode,manual,XTick,0,ws1/pi,wp1/pi,wp2/pi,ws2/pi,1);set(gca,YTickmode,manual,YTick,-50,-30,-2,0);set(gca,YTickLabelMode,manual,YTickLabels,50;

34、30; 2; 0);grid;subplot(2,3,2);plot(Wz/pi,angle(Hz)/pi);title(相频响应);xlabel();ylabel(单位：pi );axis(0,1,-1,1);set(gca,XTickMode,manual,XTick,0,ws1/pi,wp1/pi,wp2/pi,ws2/pi,1);grid;subplot(2,3,5);title(零极点图);ylabel(单位：dB);xlabel(单位：pi );zplane(DigitalB,DigitalA);axis(-1.1,1.1,-1.1,1.1);subplot(2,3,3); plot(Wz/pi,grd); title(群延迟)xlabel(频率单位：pi); ylabel(样本); axis(0,1,0,8)set(gca,XTickMode,manual,XTick,0,ws1/pi,wp1/pi,wp2/pi,ws2/pi,1); %画高刻度线 set(gca,YTickmode,manual,YTick,0:0.5:10); gridset(gcf,color,w) % 置图形背景色为白色