第章IIR数字滤波器的设计.doc

个人收集整理 勿做商业用途第6章 IIR数字滤波器的设计数字滤波器是数字信号处理技术的重要内容.和模拟滤波器一样，数字滤波器的主要功能是对数字信号进行处理，最常见的处理是保留数字信号中的有用频率成分，去除信号中的无用频率成分。本章首先介绍数字滤波器的原理和分类，使大家对数字滤波器有个整体的认识，然后介绍由模拟滤波器传递函数转换为数字滤波器传递函数的方法-脉冲响应不变法和双线性变换法。在进行MATLAB设计IIR数字滤波器的讲解之前，第三节介绍滤波器特性分析及使用方法.最后大篇幅举例介绍IIR数字滤波器的经典设计、完全设计函数和直接设计方法.6.1 概 述6。1。1 数字滤波器的工作原理数字滤波器是一个离散时间系统,输入x(n）是一个时间序列，输出y(n）也是一个时间序列。如数字滤波器的系统函数为H（z),其脉冲响应为h(n），则在时间域内存在下列的关系。 (6-1)在z域内，输入和输出存在下列关系： （6-2)式中,X（z）、Y(z)分别为输入x(n)和输出y（n）的z变换。同样在频率域内，输入和输出存在下列关系： (6-3）式中，为数字滤波器的频率特性； 和分别为x(n)和y（n)的频谱.为数字角频率，单位rad。通常设计在某些频段的响应值为1,在某些频段的响应为0.和的乘积在频率响应为1的那些频段的值仍为,即在这些频段的振动可以无阻碍地通过滤波器，这些频带为通带。和的乘积在频率响应为0的那些频段的值不管大小如何均为零，即在这些频段里的振动不能通过滤波器，这些频带称为阻带。为讨论方便起见,我们在本章中规定,为数字角频率，单位为弧度(rad),表示模拟角频率，单位弧度/秒(rad/s）.数字角频率在0范围内。一个合适的数字滤波器系统函数H(z）可以根据需要改变输入x(n）的频率特性.经数字滤波器处理后的信号y（n)保留信号x（n)中的有用频率成分，去除无用频率成分。6.1。2 滤波器的分类：按时间域特性，数字滤波器可以分为无限冲激（脉冲）响应数字滤波器（Infinite impulse response digital filter，简称IIR滤波器）和有限冲激（脉冲)响应数字滤波器（Finite impulse response digital filter，简称FIR滤波器）两类。IIR滤波器的传递函数为： (6-4)h(n)为滤波器的脉冲响应，n=0均有值。M和N为分解的分子和分母多项式的系数个数。FIR滤波器的传递函数为： (6-5）该滤波器的脉冲响应h(n）在n=0,1，,N-1的有限个点(N个点）上有值。注意（64)式中分母全为零时，H（z）具有全零点形式，IIR滤波器退化为FIR滤波器。按频率特性来讲，数字滤波器和模拟滤波器一样可分为低通、高通、带通和带阻等。数字滤波器是一个离散时间系统，在频率特性中具有周期性，因此我们讨论的频率范围仅在的范围内，相应的归一化频率在01之间，和1对应于Nyquist频率。和模拟滤波器一样，理想数字滤波器的频率特性在通带内必须满足： (66)式中,K,均为常数。和模拟滤波器一样,数字滤波器的设计目的是使滤波器的频率特性达到所给定的性能指标。其性能指标也包括通带波纹Rp(dB)、阻带衰减Rs(dB）、通带边界频率(Hz）、阻带边界频率（Hz）等。6。1。3 IIR数字滤波器的设计方法IIR滤波器的设计就是根据滤波器的性能指标要求，设计滤波器的分子和分母多项式系数。它和FIR滤波器相比优点是在满足相同性能指标要求的前提下，IIR滤波器的阶数要明显低于FIR滤波器;但IIR滤波器的相位是非线性的。在第5章我们学习了模拟滤波器的设计。在数字滤波器设计中，首先让我们想到的是能否利用模拟滤波器的设计成果进行数字滤波器设计，答案是肯定的。现在我们所讲的IIR经典设计就是将已设计好的模拟滤波器按一定变换原理转换为数字滤波器。该方法先根据滤波器的技术指标设计出相应的模拟滤波器，然后再将设计好的模拟滤波器变换成数字滤波器.在MATLAB中，经典法设计IIR数字滤波器主要采用以下步骤： 模拟离散化/bilinear,impinvar模拟低通滤波器原型/buttap,cheb1ap，cheb2ap,ellipap，besselap频率变换/lp2lp，lp2hp，lp2bp,lp2bsIIR数字滤波器 图61 IIR滤波器的设计步骤由图6-1可见，IIR滤波器设计利用了模拟滤波器的设计成果。第一步和第二步在前面已详细讨论过。第二步完成后，一个达到期望性能指标的模拟滤波器(低通、高通、带通、带阻等)已经设计出来.第三步模拟离散化的主要任务就是把模拟滤波器变换成数字滤波器，即把模拟滤波器的系统函数H（s）映射为数字滤波器的系统函数H（z）.这样，数字滤波器的设计工作全部完成。6.2 模拟滤波器到数字滤波器的转换由前面一节我们清楚了设计IIR滤波器的非常重要的一步是将设计的模拟滤波器转换为数字滤波器。实现模拟滤波器系统传递函数s域至数字滤波器传递函数z域映射的方法有冲激（脉冲）响应不变法和双线性变换法两种方法。下面不加证明地介绍这两种方法设计要点及MATLAB实现。6。2.1 脉冲响应不变法所谓脉冲响应不变法就是使数字滤波器的脉冲响应序列h（n）等于模拟滤波器的脉冲响应ha（t)的采样值，即 （67）式中,T为采样周期。因此数字滤波器的系统函数H(z）可由下式求得: （6-8)Z-表示对的内容进行Z变换，Z变换的内容请参考相应的数字信号处理教材.如果已经获得了满足性能指标的模拟滤波器的传递函数，求与之对应的数字滤波器的传递函数H(z)的方法是：（1） 求模拟滤波器的单位脉冲响应。 (6-9）式中，表示对的Laplace逆变换.Laplace变换内容请参考高等数学的积分变换或信号处理教材。（2） 求模拟滤波器单位冲激响应的采样值，即数字滤波器冲激响应序列h(n)。（3） 对数字滤波器的冲激响应h（n）进行z变换,得到传递函数H(z)。由上述方法推论出更直接地由模拟滤波器系统函数求出数字滤波器系统函数H（z)的步骤是：（1） 利用部分分式展开将模拟滤波器的传递函数H(s）展开成： （6-10）在MATLAB中这步可通过residue函数实现.若调用residue函数的形式为R,P,K=residue(a，b） 形式,则将下式（传递函数形式）： （611）变换为： （6-12）这种形式为极点留数商向量形式，对于本节所讲的特定情况，K为空矩阵。若为b，a=residue（R,P，K)则为上面调用形式的反过程。（2） 将模拟极点pk变换为数字极点即得到数字系统的传递函数 （6-13)其中T为采样间隔.(3)将(6-13）转换为传递函数形式，在该步骤中，可采用R,P,K=residue(b，a）。对于上面的步骤，MATLAB中已经提供了冲激响应不变法设计数字滤波器的函数，调用格式为:bz，az=impinvar（b，a,Fs，Fp）式中，b,a为模拟滤波器分子和分母多项式系数向量;Fs为采样频率（所滤波数据），单位Hz，缺省时为1Hz。Fp为预畸变频率（Prewarped frequency),是一个“匹配频率，在该频率上，频率响应在变换前后和模拟频率可精确匹配。一般设计中可以不考虑。bz，az分别为数字滤波器分子和分母多项式系数向量。前面已提到过,函数输入变量中的表示可添加也可略去的内容。下面我们用例子说明如何使用这个函数。【例61】脉冲响应不变法将模拟滤波器变换为数字滤波器H（z),采样周期为T=0。1s.Samp6_1b=3 2;a=2 3 1；T=0.1; 模拟滤波器分子和分母多项式系数及采样间隔bz1，az1=impinvar（b，a，1/T)程序输出为：bz1 = 0.3000 -0。2807az1 = 2。0000 3。7121 1。7214在应用冲激响应不变法设计数字滤波器时要注意它的特点。脉冲响应不变法由这一基本关系得到数字角频率和模拟角频率满足线性变换关系，T为采样间隔。这使得轴上每隔便映射到z域中的单位圆一周。如果模拟滤波器频率响应是有限带宽的话，通过变换得到的数字滤波器的频率响应非常接近于模拟滤波器的频率响应。由于数字滤波器的频率响应是模拟滤波器频率响应的周期延拓，因此对于高通和带阻滤波器存在混叠效应，会造成频率响应失真,因此这种方法原则上只适用于有限带宽滤波器。对于高通、带阻等滤波器，由于它们高频成分不衰减，势必产生严重的混迭失真。下面所讲的双线性变换法可以弥补这方面的不足。6.2.2 双线性变换法双线性变换法将s平面的整个频率轴映射到z域的一个频率周期中。因此s平面到z平面的映射是非线性的，其单值双线性映射关系为： (614）式中，T为采样周期.因此若已知模拟滤波器的传递函数,将（614）式的第一式代入即可得到数字滤波器的传递函数： (615）在双线性变换中,模拟角频率和数字角频率存在下面关系： (616）可见,模拟角频率和数字角频率之间的关系是非线性的.在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为数字滤波器的零极点和增益。num，den分别为模拟滤波器传递函数分子和分母多项式系数向量，模拟滤波器传递函数具有下面的形式： (6-17）numd和dend分别为数字滤波器传递函数分子和分母多项式系数向量。【例62】用双线性变换法将模拟滤波器变换为数字滤波器H（z），采样周期(间隔）T=0。1s.Samp6_2b=3 2；a=2 3 1;T=0.1; %模拟滤波器分子和分母多项式的系数,采样间隔bz1,az1=bilinear(b,a,1/T） 将模拟滤波器传递函数转换为数字滤波器传递函数程序输出为：bz1 = 0.0720 0.0046 -0.0674az1 = 1.0000 -1.8560 0.8606双线性变换法克服了脉冲响应不变法的频谱混迭问题,其幅值逼近程度好，可适用于高通、带阻等各种类型滤波器的设计。s域和z域对应关系也简单。缺点是频率变换的非线性导致数字滤波器与模拟滤波器在幅度和频率的对应关系上发生畸变。但一般滤波器的幅频响应具有分段常数的特点,即滤波器允许某一频段信号通过，而不允许另外频段的信号通过的特点,故变换后这一特点仍保留，影响不大。由数字边界频率计算模拟边界频率时，不是按线性关系进行的，这就是所谓的“预畸变”。但如果给定预畸变频率为边界频率，经预畸变频率校正则可以保证所要设计的模拟边界频率精确映射在所要求的数字边界频率上。6.3 滤波器特性及使用函数为了形象地看到滤波器的性能，需要绘出滤波器的幅频和相频特性曲线.本节给出MATLAB信号处理工具箱中提供分析滤波器特性的工具函数，包括时域分析（对任意输入的响应、脉冲响应、阶跃响应等)、频域分析（幅值响应、相位响应、零极点位置)和其他特性（如群延迟等)。滤波器时域和频域分析是设计各类滤波器、评价滤波器性能以及滤波器应用的基础，本节将具体介绍这些函数。6.3.1 滤波器的特性函数(1) freqz在上一章我们讲过，对于模拟滤波器，可以用freqs求解滤波器的频率响应。与之对应的函数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)。函数返回值缺省时,绘制幅频响应和相频响应图。该函数适用于下面形式的数字滤波器： (618)函数freqz输出的频率向量在0。为了获得一个滤波器真正的相频特性图,要对相位角进行解缠绕。为此MATLAB提供了一个函数unwrap来解决这个问题，P=unwrap（angle（H)）.(2) impzimpz用于产生数字滤波器的脉冲响应。调用格式为：h,t=impz（b,a,n,Fs）式中，b，a分别为滤波器分子和分母多项式系数向量；n为采样点数;Fs为采样频率，缺省值为1；h为滤波器单位脉冲响应向量；t为和h对应的时间向量。当函数输出缺省时，绘制滤波器脉冲响应图；当n缺省时，函数自动选择n值。（3）零极点图滤波器的零极点位置决定了滤波器稳定性和性能，因此考察滤波器的零极点的位置是分析滤波器特性的重要方面之一。MTALAB信号处理工具箱提供绘制数字滤波器零极点位置图的工具zplane，调用格式为：zplane(z,p）或zplane（b,a)式中,z，p为零极点向量（为复数)，b，a为滤波器分子和分母多项式的系数(为实数)。函数在z平面绘出零点和极点.极点用表示,零点用o表示。【例63】设计一个5阶Butterworth低通数字滤波器，截止频率为0。2,绘制其零极点图.Samp6_3z,p，k=butter（5，0。2）； 设计一个Butterworth数字低通滤波器zplane(z,p）; %绘制其零极点图xlabel(实部);ylabel(虚部)；程序的输出结果为图62。可以明显看出该滤波器零点和极点的分布，从而可以进行进一步的分析。由于该滤波器为一个5阶的Butterworth滤波器，故其极点的个数为5,并且极点为两两共轭，分布在单位圆内。由滤波器的分析理论知,只有极点分布在单位圆内才能保证滤波器的稳定。大家可以换用其他滤波器观看其极点在z平面上的分布。图6-2 例63设计的Butterworth滤波器的零极点图（4）群延迟 信号传输的不失真条件之一为：滤波器相频特性是一条经过原点的直线,即，为常数.但一般滤波器不满足这个条件，衡量实际滤波器相位平均延迟的物理量是群延迟。群延迟定义为信号通过滤波器的延迟随频率变化的函数,即滤波器相频特性图上切线的负斜率： （6-19）MATLAB信号处理工具箱提供计算群延迟函数grpdelay,调用格式为：gd，w=grpdelay（b，a,n,whole)gd，f=grpdelay（b,a,n,whole，Fs）gd=grpdelay（b,a，w)gd=grpdelay(b,a,f，Fs)grpdelay其中，gd为群延迟；其他各项意义同函数freqz，函数输出项缺省时，绘制群延迟图。whole参数表示绘制包括大于Nyquist频率的一个周期的群延迟。【例64】绘制截止频率为0.2的5阶Butterworth低通滤波器的群延迟。Samp6_4b,a=butter(5,0。2）； %设计Butterworth低通滤波器grpdelay(b,a,128) 绘制群延迟曲线xlabel（归一化频率）; ylabel(群延迟/采样数）图6-3 例6-4设计的Butterworth的群延迟图程序运行结果为图63.该图说明，某一频率的波具有最大的延迟。6.3.2 滤波器的使用函数（1） filter函数用来实现数字滤波器对数据的滤波，函数的调用格式为:y=filter（b，a，x）其中，b，a分别为滤波器传递函数H（z)的分子和分母多项式系数。x为滤波器的输入。y为滤波器的输出。y为与x具有相同大小的向量。（2） filtfilt函数实现零相位前向与后向结合的滤波。调用格式为：y=filtfilt（b，a,x）式中，b，a分别为滤波器传递函数H(z)的分子和分母多项式系数。x为滤波器的输入，为值向量。y为滤波器的输出。该函数对序列x进行正常的正向滤波后,将滤波后的输出翻转重新用该滤波器进行滤波，第二次滤波后的输出序列的翻转即得到零相位的滤波输出.这样就可以把延迟后的相位校正至零。但该函数只能用于数字滤波器，FIR滤波器或IIR滤波器均能使用。6.4 经典设计法6。4。1 IIR数字滤波器经典设计法 IIR数字滤波器经典设计法的一般步骤为：（1） 根据给定的性能指标和方法不同,首先对设计性能指标中的频率指标，如数字边界频率进行变换，转换后的模拟频率指标作为模拟滤波器原型设计的性能指标。（2） 估计模拟滤波器最小阶数和截止频率,利用MATLAB工具函数buttord、cheb1ord、cheb2ord、ellipord等。（3） 设计模拟低通滤波器原型。利用MATLAB工具函数buttap、cheb1ap、cheb2ap、ellipap等。（4） 由模拟原型低通滤波器经频率变换获得模拟滤波器（低通、高通、带通、带阻等)，利用MATLAB工具函数lp2lp、lp2hp、lp2bp、lp2bs。（5） 将模拟滤波器离散化获得IIR数字滤波器，利用MATLAB工具函数bilinear或impinvar。后面的几个步骤在前面的章节中已经论述，这里主要介绍第一步，关于设计性能指标的转换。设计IIR滤波器时，给出的性能指标通常分数字指标和模拟指标两种。数字性能指标给出通带截止频率，阻带起始频率，通带波纹Rp，阻带衰减Rs等。数字频率和的取值范围为0，单位弧度.而MATLAB工具函数常采用归一化频率,和的取值范围为01，对应于0，此时需进行转换。模拟性能指标给出通带截止频率,阻带起始频率,通带波纹Rp，阻带衰减Rs等。模拟频率和单位为弧度/秒（rad/s）。MATLAB信号处理工具箱中，设计性能指标的转换应根据不同设计方法进行不同处理.下面就举例介绍这些方法。6.4。2 冲激(脉冲）响应不变法采用来将和变换为和【例65】用脉冲响应不变法设计一个Butterworth低通数字滤波器，使其特征逼近一个低通Butterworth模拟滤波器的下列性能指标：通带截止频率，通带波纹Rp小于3dB，阻带边界频率为,阻带衰减大于15dB，采样频率Fs=10000Hz。假设一个信号，其中f1=1000Hz,f2=4000Hz.试将原信号与通过该滤波器的输出信号进行比较。Samp6_5Wp=2000*2pi；Ws=3000*2*pi； %滤波器截止频率Rp=3;Rs=15; %通带波纹和阻带衰减Fs=10000; %采样频率Nn=128; %调用freqz所用的频率点数N,Wn=buttord（Wp，Ws,Rp,Rs,s）; %模拟滤波器的最小阶数z，p，k=buttap（N）； 设计模拟低通原型Butterworth滤波器Bap，Aap=zp2tf(z，p，k）； 将零点极点增益形式转换为传递函数形式b,a=lp2lp(Bap，Aap，Wn); 进行频率转换bz,az=impinvar（b，a,Fs); %运用脉冲响应不变法得到数字滤波器的传递函数figure(1）H，f=freqz(bz，az,Nn，Fs)； 求解数字滤波器的幅频特性和相频特性subplot(2,1,1）,plot（f，20*log10（abs(H)xlabel（频率/Hz）;ylabel（振幅/dB);grid on；subplot(2，1,2),plot(f，180/piunwrap(angle(H）)xlabel（频率/Hz）；ylabel（相位/o)；grid on;figure(2）f1=1000;f2=4000； 输入信号的频率N=100； 数据长度dt=1/Fs；n=0:N-1；t=n*dt; %采样间隔和时间序列x=sin（2pi*f1t)+0。5*cos（2*pif2*t)； %滤波器输入信号subplot（2，1，1),plot(t，x），title(输入信号) %绘制输入信号y=filtfilt（bz，az，x); 用函数filtfilt对输入信号进行滤波y1=filter（bz，az，x)； %用filter函数对输入信号滤波subplot（2,1,2),plot(t，y，t,y1，：），title（输出信号)，xlabel(时间/s）legend( filtfilt ， filter) %加图例程序的运行结果为图6-4和图65。由图6-4上图可知,在小于2000Hz处的衰减小于3dB，而在大于3000Hz处衰减大于15dB，满足滤波器的设计指标。由图65可见滤波器对含有1000Hz和4000Hz频率成分的信号进行了滤波，滤除了4000Hz的信号。由程序的输出还可以看出，采用filtfilt函数，输出的1000Hz信号(实线）与输入1000Hz的信号相位一致，即经过滤波后并没有改变信号波形形状。而运用filter函数滤波后（虚线)有一些延迟，改变了信号的形状。图 6-4 例6-5所设计Butterworth滤波器的频率响应。上图:幅频特性；下图:相频特性图65例6-5设计滤波器的输入和输出信号输出信号中实线表示运用filtfilt函数的输出，虚线表示运用filter函数的输出【例6-6】用脉冲响应不变法设计Butterworth低通数字滤波器，要求通带频率为，通带波纹小于1dB，阻带在内，幅度衰减大于15dB，采样周期T=0.01s。假设一个信号，其中f1=5Hz，f2=30Hz。试将原信号与经过该滤波器的输出信号进行比较.注意该题给出的通带边界频率和阻带边界频率均为数字频率，因此设计时首先要将其转换为模拟频率。%Samp6_6wp=0.2*pi;ws=0.3pi；Rp=1;Rs=15; %数字滤波器截止频率、通带波纹和阻带衰减T=0.01;Nn=128; 采样间隔Wp=wp/T;Ws=ws/T； %得到模拟滤波器的频率采用脉冲响应不变法的频率转换形式N，Wn=buttord（Wp，Ws，Rp，Rs，s)； 计算模拟滤波器的最小阶数z，p,k=buttap（N）; 设计低通原型数字滤波器Bap,Aap=zp2tf（z，p，k)； 零点极点增益形式转换为传递函数形式b,a=lp2lp(Bap,Aap，Wn)； 低通滤波器频率转换bz,az=impinvar（b，a,1/T）; 脉冲响应不变法设计数字滤波器传递函数figure(1）H,f=freqz（bz，az,Nn,1/T)； 输出幅频响应和相频响应subplot(2,1，1）,plot（f,20log10(abs(H）)；xlabel(频率/Hz)；ylabel(振幅/dB);grid on;subplot（2,1,2),plot(f，180/pi*unwrap(angle(H）)xlabel（频率/Hz);ylabel（相位/o);grid on;figure（2)f1=5;f2=30; %输入信号含有的频率N=100； 数据点数n=0：N-1；t=n*T; %时间序列x=sin（2*pi*f1*t)+0。5cos（2pi*f2t）； %输入信号subplot（2,1，1)，plot（t，x)，title(输入信号）y=filtfilt(bz,az，x）； 对信号进行滤波subplot（2，1,2),plot(t，y）,title（输出信号)，xlabel（时间/s) 程序的运行结果为图66和图6-7。该例所要求的通带频率为,而该滤波器的采样间隔为0。01s,采样频率为100Hz，即2对应于100Hz，则0.2对应于10Hz，即该滤波器的通带范围为010Hz.观看图66上图，其通带范围最大衰减小于1dB，与该分析一致.题意要求阻带在内，幅度衰减大于15dB。其中0.3对应于15Hz,幅度衰减大于15dB.完全符合滤波器设计的要求.测试信号中含有5Hz和30Hz的信号，显然，5Hz可以通过该滤波器，而30Hz的信号不能通过该滤波器。从输出信号可以看出，滤波器滤去了30Hz的高频信号,达到了滤波要求.图6-6 例6-6设计的滤波器的频率特性上图：幅频特性；下图：相频特性图6-7例66设计滤波器的输入和输出信号6.4.3 双线性变换法采用将和变换为和【例6-7】用双线性变换法设计一个椭圆低通滤波器，其性能指标同例66。Samp6_7wp=0.2*pi;ws=0。3pi;Rp=1;Rs=15; %数字滤波器截止频率通带波纹和阻带衰减Fs=100;Ts=1/Fs；Nn=128； %采样频率Wp=2/Tstan(wp/2.);Ws=2/Tstan（ws/2。); %按频率转换公式进行转换N，Wn=ellipord（Wp,Ws,Rp，Rs，s)； %计算模拟滤波器的最小阶数z,p,k=ellipap（N,Rp,Rs)； %设计模拟原型滤波器Bap，Aap=zp2tf（z,p,k)； %零点极点增益形式转换为传递函数形式b，a=lp2lp(Bap，Aap，Wn）； 低通转换为低通滤波器的频率转换bz,az=bilinear（b，a，Fs）; 运用双线性变换法得到数字滤波器传递函数H，f=freqz(bz,az，Nn,Fs); %求出频率特性subplot(2，1,1)，plot(f,20*log10(abs(H）)；xlabel（频率/Hz);ylabel（振幅/dB)；grid on；subplot(2,1，2）,plot(f，180/piunwrap（angle（H)）)xlabel（频率/Hz)；ylabel（相位/o);grid on;图68 例67设计椭圆低通滤波器的频率特性上图:幅频特性;下图：相频特性程序运行结果见图6-8。在10Hz以前，衰减小于1dB，在15Hz以后衰减均大于15dB，即性能指标完全满足滤波器的设计要求。6。5 IIR滤波器的完全设计函数以上介绍了IIR滤波器设计原理和基本方法步骤，并给出了一些例子说明如何用MATLAB编程分步实现这些步骤。由这些步骤可知我们必须多次调用MATLAB信号处理工具箱中的基本工具函数.实际上，MATLAB信号处理工具箱还提供了IIR滤波器设计的完全工具函数,用户只要调用这些工具函数即可一次性完成设计，而不需要调用那些基本工具函数分步实现。IIR滤波器设计的完全工具函数有butter，cheby1,cheby2，ellip。这些工具函数既可用于设计模拟滤波器,也可用于设计数字滤波器。关于这些函数设计模拟滤波器已在前面章节中述及，这里介绍这些函数在IIR数字滤波器中的应用。在这两类滤波器设计中，这些工具函数调用格式基本相同，只是在频率处理上有所不同。在MATLAB滤波器设计工具箱中，数字滤波器采用归一化频率，取值为01之间，归一化频率1对应的数字角频率为，对应的真实频率为采样频率的一半。在应用MATLAB工具函数设计数字滤波器时应注意这一点。数字IIR滤波器的完全设计函数有： b，a=butter(n，wn,ftype) z,p,k=butter(n，wn, ftype) b,a=cheby1（n,Rp,wn，ftype） z,p,k=cheby1（n, Rp，wn,ftype) b，a=cheby2（n，Rs,wn,ftype) z，p,k=cheby2（n， Rs,wn，ftype) b，a=ellip(n,Rp，Rs，wn，ftype） z，p，k=ellip（n, Rp，Rs，wn，ftype)在上面的调用方式中,n为滤波器的阶数，wn为滤波器的截止频率，取值为01。需根据采样频率Fs来定，如滤波器的截止频率为Fc（Hz)，则wn的计算公式为： (6-20）这样就转换为01的归一化频率.其中wp,ws等边界频率都要根据此公式进行转换。 ftype滤波器的类型为： high为高通滤波器，截止频率为wn。 stop为带阻滤波器，截止频率为wn=w1,w2 （ww2). ftype缺省时为低通或带通滤波器.a，b分别为滤波器传递函数分子和分母多项式系数向量;z,p,k分别为滤波器的零极点和增益。Rp,Rs分别为所设计滤波器的通带波纹和阻带衰减，单位为dB。设计好的数字滤波器传递函数具有下面形式： （621）上述函数采用双线性变换法和频率的预畸变处理将模拟滤波器离散化为数字滤波器，同时保证模拟滤波器和数字滤波器在wn(或w1，w2)处具有相同的幅频响应。设计时应注意真实频率和MATLAB归一化数字频率之间的转换,即620式的应用。在进行IIR数字滤波器设计之前，请大家注意在模拟滤波器中我们用求取最小阶数和截止频率的函数如buttord，cheb1ord,cheb2ord,ellipord，这些函数完全可以用于IIR数字滤波器设计中，只不过在模拟滤波器设计中需加可选项s，在数字滤波器中则不加该项。另外，输出的截止频率也是归一化频率（归一化为01).【例6-8】设计一个Butterworth高通数字滤波器，通带边界频率为300Hz,阻带边界频率为200Hz,通带波纹小于1dB，阻带衰减大于20dB，采样频率为1000Hz.假设一个信号，其中f1=10Hz，f2=400Hz。试将原信号与通过该滤波器的输出信号进行比较。Samp6_8Fs=1000； 采样频率wp=3002/Fs;ws=200*2/Fs； 根据采样频率将滤波器边界频率进行转换(619式)Rp=1;Rs=20; 通带波纹和阻带衰减Nn=128； 显示滤波器频率特性的数据长度N,Wn=buttord(wp,ws,Rp，Rs)； %求得数字滤波器的最小阶数和截止频率（归一化频率)b，a=butter(N,Wn,high）; 设计Butterworth高通数字滤波器figure（1）H,f=freqz(b,a,Nn,Fs）; 用Nn点绘出频率特性subplot（2,1，1)，plot(f，20log10（abs（H）)；xlabel(频率/Hz)；ylabel（振幅/dB）；grid on;subplot（2,1,2），plot(f，180/piunwrap（angle(H）xlabel（频率/Hz）；ylabel(相位/o）;grid on;n=0:127；dt=1/Fs;t=ndt; %时间序列f1=10；f2=400； 输入信号频率x=sin（2*pi*f1*t）+0.5*cos（2pi*f2*t）； 输入信号figure（2）subplot（2，1，1）,plot（t,x); title(输入信号） %绘制输入信号y=filter(b，a，x)； %对输入信号进行滤波subplot(2，1，2),plot(t,y）,title（输出信号) 绘制输出信号xlabel(时间/s） 图69例68所设计的高通Butterworth滤波器的频率特性上图:幅频特性；下图：相频特性图6-10 例68所设计滤波器的输入和输出信号程序输出结果为图69和图6-10。由图69上图可以看到所设计滤波器在大于300Hz为通带，其衰减均小于1dB；小于200Hz为阻带，其衰减大于20dB，完全符合设计要求,但相频特性是非线性的（图6-9下图）。由图610可以看出，当滤波器输入10Hz和400Hz两种信号后,滤波器滤除了10Hz的信号，使得400Hz的信号通过了滤波器，起到了滤波的效果。【例6-9】设计一个带通Chebyshev I型数字滤波器,通带为100Hz200Hz，过渡带宽均为50Hz，通带波纹小于1dB，阻带衰减大于30dB，采样频率Fs=1000Hz。假设一个信号，其中f1=30Hz,f2=100Hz，f3=270Hz。试将原信号与经过该