基于matlab的IIR数字滤波器的设计

1、基于matlab的IIR数字滤波器的设计 作者： 日期：15 个人收集整理 勿做商业用途 基于matlab的IIR数字滤波器的设计王佳妮(渭南师范学院 物理与电气工程学院 电气工程及其自动化2008年1班)摘 要：IIR数字滤波器在MATLAB环境下的设计方法和实现方法，在无限脉冲响应(IIR)数字滤波器设计中，先进行模拟滤波器的设计,然后进行模拟数字滤波器转换，即采用脉冲响应不变法及双线性Z变化法设计数字滤波器，最后进行滤波器的频带转换.关键词：IIR数字滤波器；matlab；频带转换;引言 数字滤波器是数字信号处理的重要基础,数字信号处理主要是研究数字或符号的序列表示信号波形，并用数字的方

2、式去处理这些序列，把它们改变成在某分量和中意义上更希望的形式，以便估计信号的特征参量，或削弱信号中的多余分量和增强信号中的有用分量。数字滤波器在对信号的过滤、检测与参数估计等处理过程中，是使用最为广泛的一种线性系统。 滤波器的种类很多，从功能上可以分为低通、高通、带通和带阻滤波器，上述每种滤波器又可以分为模拟滤波器和数字滤波器。如果滤波器的输入输出都是数字信号，则这样的滤波器称之为数字滤波器，它通常通过一定的运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分来实现滤波。根据数字滤波器冲激响应的时域特性，可将数字滤波器分为两种，即无限长冲激响应(IIR）滤波器和有限长冲激响应（FIR

3、)滤波器。有数字信号处理的一般理论可知，IIR滤波器的特征是具有无限持续时间的冲激响应，而FIR滤波器使冲激响应只能持续一定的时间.随着信息时代的到来，数字信号处理已经成为当今一门极其重要的学科和技术，并且在通信、语音、图像、自动控制等众多领域得到了广泛的应用。在数字信号处理中,数字滤波器占有极其重要的地位，它具有精度高、可靠性好、灵活性大等特点。现代数字滤波器可以用软件或硬件两种方式来实现.软件方式实现的优点是可以通过滤波器参数的改变去调整滤波器的性能。MATLAB是一种面向科学和工程计算的语言，它集数值分析、矩阵运算、信号处理和图形显示于一体，具有编程效率高、调试手段丰富、扩充能力强等特点

4、.MATLAB的信号处理工具箱具有强大的函数功能，它不仅可以用来设计数字滤波器，还可以使设计达到最忧化，是数字滤波器设计的强有力工具。MATLAB将矩阵运算、数值分析、图形处理、编程技术结合在一起，为用户提供了一个强有力的科学及工程问题的分析计算和程序设计工具,它还提供了专业水平的符号计算、文字处理、可视化建模仿真和实时控制等功能,是具有全部语言功能和特征的新一代软件开发平台。     MATLAB 已发展成为适合众多学科,多种工作平台、功能强大的大型软件。在欧美等国家的高校，MATLAB已成为线性代数、自动控制理论、数理统计、数字信号处理、时间序列分析、动态系统

5、仿真等高级课程的基本教学工具。成为攻读学位的本科、硕士、博士生必须掌握的基本技能。在设计研究单位和工业开发部门,MATLAB被广泛的应用于研究和解决各种具体问题。在中国，MATLAB也已日益受到重视,短时间内就将盛行起来，因为无论哪个学科或工程领域都可以从MATLAB中找到合适的功能. 1.1.1 滤波器原理滤波器是对输入信号进行滤波，对于滤波器LSI系统如图11所示，其时域输入输出关系为：y（n）=x（n)h(n） （1-1) 图11 若y(n)，x（n)的傅里叶变换存在，则输入输出的频域关系是： 当输入信号x（n）通过滤波器系统h（n）后，其输出y(n）中不再含有w>wc的频域成分，

6、使|w|wc的信号成分通过,因此，滤波器的形状不同，其滤波后的信号结果也不一样。 若滤波器的输入、输出都是离散时间信号，那么该滤波器的单位冲击响应h(n）也必然是离散的，这种滤波器称为数字滤波器（DF，Digital Filter）。当用硬件实现一个DF时，所需原件是延迟器、乘法器和加法器；而利用MATLAB软件时，它仅需线性卷积程序便可以实现。而模拟滤波器只能用硬件实现。因此比容易实现,且更容易获得理想的滤波性能。数字滤波器还具有以下优点：精度和稳定性高;系统函数容易改变；灵活性高；不存在阻抗匹配问题;便于大规模集成；可实现多维滤波。1。1。2 数字滤波器从功能上分类：可以分为低通滤波器（L

7、P）、高通滤波器(HP）、带通滤波器（BP）、带阻滤波器(BS）。从滤波器的网络结构或者从单位脉冲响应分类:如同模拟滤波器的性能可由g（t)和G（s)来表征一样,数字滤波器的性能完全取决于h(n)和H（z)。因此，无限脉冲响应数字滤波器，简称IIR数字滤波器，它的h(n）序列长度是无限的，即当n趋于无穷大时，h（n）仍有效。数字滤波器的设计与实现,通常按下述步骤进行:1、根据不同用途提出数字滤波器的技术指标、性能要求。2、设计一个稳定的、因果的数字模型H（z)来逼近所要求的技术指标，并用有精度的运算实现所设计的系统。3、设计专用的数字硬件来实现这个数字模型，或者通过电子计算机运行程序软件予以实

8、现。IIR滤波器的设计方法有两类,经常用到的一类设计方法是借助于模拟滤波器的设计方法进行的。其设计思路是：先设计模拟滤波器得到传递函数G(s)，然后将G（s）按某种方法装换为数字滤波器的系统函数H（z）。这一类是基于模拟滤波器的设计方法相对比较成熟,它不仅有完整的设计公式，也有完整的图标供查阅，更可以直接调用中的对应的函数进行设计。另一种是直接在频域或者时域中进行设计,设计时必须使用计算机辅助，直接调用LAB中的程序或函数即可设计。2IIR数字滤波器设计及其MATALB实现IIR数字滤波器技术设计最通用的方法是借助于模拟滤波器的设计方法。模拟滤波器设计已经有了一套相当成熟的方法，它不但有完整的

9、设计公式，而且还有较为完整的图表可供查询,因此充分利用这些已有的资源将会给数字滤波器的设计带来很大方便。滤波器设计最重要的是寻找一个稳定、因果的系统函数去逼近滤波器的技术指标，因此模拟滤波器的设计十分重要.实际中，有三种广泛应用的滤波器，即巴特沃斯低通滤波器、切比雪夫低通滤波器，椭圆低通滤波器，下面，将分别予以介绍。2.1.1 巴特沃斯低通滤波器设计巴特沃斯滤波器拥有最平滑的频率响应，在截止频率以外，频率响应单调下降。在通带中是理想的单位响应，在阻带中响应为0。在截止频率处有3dB的衰减。巴特沃斯滤波器除了具有平滑单调递减的频率响应的优点之外，其过渡带的陡峭程度正比于滤波器的阶数。高阶巴特沃斯

10、滤波器的频率响应近似于理想低通滤波器。巴特沃斯(Butterworth)低通滤波器是将巴特沃斯函数作为滤波器的传递函数，它的平方幅度函数为:式中,为滤波器频率，为3dB截止频率,N表示滤波器的阶次。MATLAB信号处理工具箱为低通模拟巴特沃斯滤波器的产生提供了函数buttap，其调用的格式为：z,p，k=buttap(N），其中,z表示零点,p表示极点，k表示增益，N表示阶次.%巴特沃斯低通模拟滤波器clear all;n=0：0.01:2;for i=1:4switch i case 1 N=2；case 2 N=5; case 3 N=10； case 4 N=20; endz,p，k=b

11、uttap（N); b，a=zp2tf(z,p，k）; H，w=freqs（b，a，n）; magH2=(abs（H）)。2; hold on plot(w，magH2） axis(0 2 0 1）； endxlabel('w/wc);ylabel('|H（jw)|2'）；title(巴特沃斯低通模拟滤波器)；grid on由图21可清楚地分析出,巴特沃斯滤波器拥有平滑的频率响应,在截止频率以外，频率响应单调下降。其过渡带的陡峭程度正比于滤波器的阶数，高阶巴特沃斯滤波器的频率响应近似于理想低通滤波器.由图2-1所示，对N=2阶滤波器和N=20阶滤波器的平方幅度响应进行比

12、较后，证明了高阶巴特沃斯滤波器有着更好的幅度特性,更接近理想低通滤波器。2。1。2 切比雪夫低通滤波器设计巴特沃斯滤波器的频率特性在通带和阻带内部都是随频率单调变化的，如果在通带的边缘能够满足指标，那么在通带的内部会超过设计的指标要求，造成滤波器的阶次N比较高。因此，有效地设计方法是将指标的精确度均匀的分布在整个通带内,或者均匀分布在整个阻带内,或者同时分布在两者之内，这样就可用阶数较低的系统来满足要求。这可通过选择具有等波纹特性的逼近函数来达到。切比雪夫滤波器的振幅特性就具有这种等波纹特性。它有两种形式：振幅特性在通带内是等波纹的、在阻带内是单调的切比雪夫I型滤波器；振幅特性在通带内是单调的

13、、在阻带内是等波纹的切比雪夫II型滤波器。采用何种形式的切比雪夫滤波器取决于实际用途.切比雪夫I型(Chebyshev-I）滤波器平方幅度响应函数为： （22)其中是一个小于1的正数,它与通带波纹有关，越大，波纹也越大;是的切比雪夫多项式，它的定义为： 而切比雪夫II型（ChebyshevII）滤波器平方幅度响应函数为： (24)其中是一个小于1的正数，表示波纹变化情况;为截止频率；N为滤波器的阶次，也是的阶次. 图2-2和图23分别画出了理想时的切比雪夫I型与切比雪夫II型滤波器阶次N为奇数与偶数时的幅频特性。而通过MATLAB信号处理工具箱中的函数cheb1ap及cheb2ap,可以实现切

14、比雪夫滤波器设计，其调用格式为:z，p，k=cheb1ap（N,Rp） z，p,k=cheb2ap(N，Rs）其中，z表示零点,p表示极点，k表示增益，N表示阶次,Rp为通带波纹（dB）， Rs为阻带波纹（dB).切比雪夫I型低通模拟滤波器clear all； n=0：0.01:2;for i=1:2 switch i case 1 N=7； case 2 N=8; endRp=1； z，p,k=cheb1ap（N，Rp); b，a=zp2tf（z,p，k); H,w=freqs（b,a,n)； magH2=(abs（H）).2；%绘制图形posplot=12 num2str(i）; subp

15、lot(posplot） plot（w,magH2)axis（0 2 0 1）；xlabel('w/wc）;ylabel('|H（jw)|2')；title(N=' num2str（N);end切比雪夫II型低通模拟滤波器clear all； n=0：0.01:2;for i=1：2 switch i case 1 N=7； case 2 N=8； endRs=10; z，p,k=cheb2ap(N,Rs）； b，a=zp2tf(z，p，k);H，w=freqs(b，a，n）；magH2=（abs（H）。2；输出图形posplot='12' nu

16、m2str(i);subplot（posplot）plot(w，magH2）axis(0 2 0 1。1);用图24、图25 与图22、2-3分别作比较可清楚地分析得出,切比雪夫滤波器有效地将技术指标的精确度均匀分布，切比雪夫I型滤波器的振幅特性在通带内是等波纹的、在阻带内是单调的;切比雪夫II型滤波器的振幅特性在通带内是单调的、在阻带内是等波纹的。2.1。3 椭圆低通滤波器设计切比雪夫I型滤波器在通带内成等波纹振荡，在阻带内却仍是单调下降的，切比雪夫II型在阻带内是等波纹的，在通带内却是单调下降的。因此过渡带的特性有所提高,但是并不理想。它的主要原因在于两者的系统函数在截止频率附近没有有限个

17、零点，其零点在无限远处。1931年，考尔提出了采样有限零点设计的滤波器，因为这种方法在确定零点位置时与椭圆函数的许多特性有关,所以称之为椭圆低通滤波器。它的平方幅度响应函数为： 式中是雅可比椭圆函数,是与通带衰减有关的函数。滤波器阶次N等于通带和阻带内最大点和最小点的和.MATLAB信号处理工具箱为低通模拟椭圆滤波器的产生提供了函数ellipap，其调用的格式为：z,p,k= ellipap （N，Rp,Rs），其中，z表示零点，p表示极点，k表示增益，N表示阶次,Rp为通带波纹（dB）, Rs为阻带波纹（dB）。%椭圆低通滤波器clear all； n=0：0。01:2；for i=1：4

18、switch i case 1 N=2； case 2 N=3； case 3 N=4； case 4 N=5； endRp=1； Rp-通带纹波系数为1 dBRs=10； %Rs-阻带纹波系数为10dBz,p,k=ellipap（N,Rp，Rs)； 函数ellipap设计椭圆低通滤波器b，a=zp2tf(z，p，k)；H，w=freqs（b,a，n）;magH2=（abs(H）)。2;%输出图形posplot=22num2str（i)；subplot（posplot)plot(w，magH2)axis（0 2 0 1.1);xlabel('w/wc）；ylabel（'H（jw

19、)|2)；title（N=num2str（N);grid onend由图2-6可得出，椭圆滤波器在通带及阻带内均成等波纹振荡,并具有最窄的过渡带,其性质比较优越,但椭圆滤波器的转移函数不但有极点，而且在轴上还有零点,所以椭圆滤波器的设计比较复杂。2.2 模拟数字滤波器变换及其MATLAB实现在设计了模拟低通滤波器后，就可以把它们变成数字滤波器了。这些变换均是复值映射,许多文献对此都有研究，根据数字滤波器所保持的模拟滤波器的不同特性，研究出不同的变换技术。其中，最重要的有两种:脉冲响应不变法（保持脉冲响应不变，又叫冲激响应不变法）和双线性Z变换法（保持系统函数不变).2.2。1 脉冲响应不变法脉

20、冲响应不变法的设计原理是使数字滤波器的单位抽样响应序列h(n)，模仿模拟滤波器的脉冲响应g（t）。设系统传递函数为G(s）的模拟滤波器的单位脉冲响应为g（t)，并将脉冲响应g(t）进行等间隔采样，使得数字滤波器的单位抽样响应h（n)刚好等于g(t)的采样值，即： 其中的为采样周期。G(s)是模拟滤波器的系统传递函数，又令H(z）是数字滤波器的系统传递函数。采样信号的拉式变换与相应的采样序列Z变换的映射关系为：所以系统函数G（s)和H（z)的关系为: （2-8)式(28)的物理意义为首先将模拟滤波器的系统函数G(s)作周期的延拓，在经过式（2-7)的映射变换，映射到Z平面上,从而得到数字滤波器的

21、系统函数H(z)。且模拟和数字频率满足下列关系：=T。经过式(2-7)的映射，s平面的左半平面映射为Z平面的单位圆内，因此，一个因果的和稳定的模拟滤波器映射成因果的和稳定的数字滤波器.经过以上分析，按照脉冲响应不变法，通过模拟滤波器的系统传递函数G（s),可直接求得数字滤波器系统函数H(Z)，其设计具体步骤归纳如下：(1）利用=T(可由关系式推出），将数字滤波器指标，转换为模拟滤波器指标,（2）根据指标，来设计模拟滤波器G（s)（3)利用部分分式展开法,把G（s）展成 (29)（4）最后把模拟极点转换为数字极点，得到数字滤波器： （210)根据上述理论，将举例在MATLAB环境下用函数实现脉冲

22、响应不变法设计一数字低通滤波器。其函数为b,a=impinvar(c，d，T），其中,b表示数字滤波器自变量为的分子多项式,a表示数字滤波器自变量为的分母多项式,c表示模拟滤波器自变量为s的分子多项式，d表示模拟滤波器自变量为s的分母多项式，T表示采样变换参数。设模拟低通巴特沃斯滤波器，通带纹波为Rp=1dB,通带上限角频率p=0。2,阻带下限角频率s=0.3，带阻最小衰减s=15dB，根据该低通模拟滤波器，利用冲激响应不变法设计相应的数字低通滤波器wp=0.2pi；ws=0.3pi;Rp=1；As=15；T=1； 采样周期%性能指标Rip=10（Rp/20）;Atn=10(As/20)；Om

23、gP=wp*T；OmgS=wsT;N，OmgC=buttord（OmgP，OmgS,Rp,As,s）; %选取模拟滤波器的阶数cs，ds=butter（N,OmgC,s'); %设计出所需的模拟低通滤波器b,a=impinvar(cs,ds，T）； 应用脉冲响应不变法进行转换 db,mag,pha,grd,w=freqz_m(b,a); %求得相对，绝对频响及相位、群延迟响应绘制各条曲线subplot(2，2，1）;plot（w/pi,mag）;title(幅频特性）；xlabel('w(/pi）)；ylabel('|H（jw)；axis（0，1,0，1。1）；set(

24、gca，'XTickMode','manual,XTick，0 0。2 0。3 0。5 1）；set（gca，'YTickMode'，manual，'YTick,0 Atn Rip 1）;gridsubplot（2,2,2）；plot(w/pi，db）；title(幅频特性(dB);xlabel(w(/pi)）;ylabel(dB);axis（0，1,-40，5);set（gca，XTickMode,manual'，XTick,0 0。2 0。3 0.5 1）;set(gca,'YTickMode，'manual，'

25、;YTick，40 -As -Rp 0）;gridsubplot(2，2,3);plot(w/pi,pha/pi);title（'相频特性')；xlabel(w（/pi)'）;ylabel('pha(/pi)'）；axis(0，1，-1，1）;set（gca,XTickMode，manual，'XTick,0 0.2 0。3 0。5 1）；gridsubplot(2，2，4);plot（w/pi,grd）；title（'群延迟）；xlabel（'w（/pi);ylabel(Sample)；axis（0,1,0,12）；set(gc

26、a，XTickMode',manual'，'XTick,0 0.2 0.3 0.5 1）；grid%本设计用到的M文件函数functiondb，mag,pha，grd,w=freqz_m（b，a）滤波器幅值响应（绝对、相对）、相位响应及群延迟Usage：db,mag，pha,grd,w=freqz_m(b,a） %500点对应0,pidb 相对幅值响应； mag 绝对幅值响应； pha 相位响应； grd 群延迟响应;%w 采样频率; b 系统函数H（z）的分子项a 系统函数H（z）的分母项H，w=freqz(b，a，500）; 500点的幅频响应mag=abs(H)；

27、db=20*log10（mag/max(mag）；pha=angle（H）;grd=grpdelay(b，a,w）;总结以上，脉冲响应不变法的优点是频率坐标变换是线性的，即=T，如不考虑频率混叠现象，用这种方法设计数字滤波器会很好的重现原模拟滤波器的频率响应。另外一个优点是数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，时域逼近好。但其也具有很大的缺点，若抽样频率不高或其它原因将产生混叠失真，不能重现原模拟滤波器频率响应.所以，脉冲响应不变法适合低通、带通滤波器设计，不适合高通、带阻滤波器的设计。2。2.2 双线性Z变换法利用脉冲响应不变法设计数字滤波器时,由于=T的频率关系是根据推

28、导的,所以是轴每隔2/T便映射到单位圆上一周,引起了频域混叠。为克服这一现象,人们找到了另一种映射关系： 此关系称为双线性Z变换法。双线性Z变换法的基本思路是:首先将整个s平面压缩到s1平面的一条带宽为2/T（丛-/T到/T)的横带里，然后通过标准的变化关系将横带变换成整个Z平面上去，这样就得到s平面与Z平面间的一一对应的单值关系,整个过程如图28所示：图28 双线性Z变换法的映射关系由式（211）得 (2-12)及 （2-13） （2-14)式(2-11)及式(2-12）给出了s和z之间的映射关系，而式（2-13） 和式（214）给出了和之间的映射关系,但这是一种非映射关系，双线性Z变换法正

29、是利用了正切函数的非线性特点,把整个j轴压缩到了单位圆的一周上。在MATLAB中，双线性Z变换可通过bilinear函数实现,其调用格式为：Bz，Az=bilinear（B，A，Fs),其中B，A为模拟滤波器的传递函数G（s）的分子分母多项式系数分量,而Bz，Az为数字滤波器的传递函数H（Z）的分子分母多项式的系数分量。%使用双线性Z变换设计一低通数字滤波器，fp=100HZ，fs=300HZ,p=3dB,s=20dB,抽样频率为Fs=1000HZfp=100；fs=300;Fs=1000;rp=3；rs=20; wp=2pifp/Fs；ws=2pi*fs/Fs;Fs=Fs/Fs %使Fs为一wap=tan（wp/2)；was=tan(ws/2)；n，wn=buttord(wap，was