IIR和FIR数字滤波器的设计及其结构研究2

IIR和FIR数字滤波器的设计及其结构研究【摘要】数字滤波器是指输入、输出均为数字信号，通过数值运算处理改变输入信号所含频率成分的相对比例，或者滤除某些频率成分的数字器件和程序。经典数字滤波器从滤波特性上分类，可以分成低通、高通、带通和带阻等滤波器。根据数字滤波器冲激响应的时域特性，可以分成无限脉冲响应数字滤波器(简称IIR)和有限脉冲响应数字滤波器(简称FIR)，IIR和FIR数字滤波器的设计方法及其结构各不相同。本次课程设计先是对数字滤波器有关理论知识作介绍，在性能指标分析基础上分别对IIR带通数字滤波器和FIR低通数字滤波器运用MATLAB相关函数设计程序，得到幅频特性曲线图像，并对结果进行分析，最后总结课程设计。1、引言随着信息技术的迅猛发展，数字信号处理已成为一个极其重要的学科和技术领域。在通信、语音、图像、自动控制和家用电器等众多领域得到了广泛的应用。数字滤波(DigitalFilter，DF)是数字信号处理的重要环节，它在数字信号处理中占有着重要的地位，它具有可靠性好、精度高、灵活性大、体积小、重量轻等优点。随着数字技术的发展，数字滤波器越来越受到人们的重视，广泛地应用于各个领域。数字滤波器的输入输出信号都是数字信号，它是通过一定的运算过程改变输入信号所含频率成分的相对比例或者滤除某些频率成分来实现滤波的，这种运算过程是由乘法器、加法器和单位延迟器组成的。根据数字滤波器冲激响应的时域特性，可将数字滤波器分为两种即无限长冲激响应(IIR)滤波器和有限长冲激响应(FIR)滤波器。由数字信号处理的一般理论可知IIR滤波器的特征是具有无限持续时间的冲激响应，而FIR滤波器使冲激响应只能持续一定的时间。滤波器的功能是用来移除信号中不需要的部分，比如随机噪声；或取出信号中的有用部分，如位于某段频率范围内的成分。目前随着计算机技术和数字信号处理器芯片的发展，使我们更为便利地识别和提取各种各样的信号。因此研究不同数字滤波器的设计原理和稳定性分析对于满足军事、航空、民营等等各个领域的信号处理要求具有十分重要的意义。第二部分：IIR数字滤波器设计的基本过程（模拟滤波器的设计及方法）、设计方法（脉冲响应不变法和双线性变换法）及两种设计方法的优缺点比较。2、IIR数字滤波器设计方法概述2.1IIR数字滤波器设计 IIR滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的设计方法进行的。模拟滤波器设计是基础，模拟滤波器到数字滤波器的转换是核心。而模拟滤波器的设计都是通过低通滤波器来实现，比较常用的模拟低通滤波器有Butterworth(巴特沃斯)和Chebyshev(切比雪夫)等。将模拟滤波器变换为数字滤波器的主要方法有脉冲响应不变法和双线性变换法。直接法直接在频域或者时域中设计数字滤波器，由于要解联立方程，设计时需要计算机辅助设计。借助于模拟滤波器的数字滤波器设计流程如图2.1所示。图2.1IIR数字滤波器设计过程2.2模拟滤波器设计模拟低通滤波器的设计是设计其他滤波器的基础。模拟高通、带通和带阻滤波器的设计过程是：先将希望设计的各种滤波器的技术指标转换为低通滤波器技术指标，然后设计响应的低通滤波器，最后采用频率转换法将低通滤波器转换成所希望的各种滤波器。模拟滤波器设计流程如图2.2所示。图2.2模拟滤波器设计过程2.2.1Butterworth模拟低通滤波器常用于待设计的原型低通滤波器，设计步骤为：由滤波器的设计指标、、、和式(2-1)确定滤波器的阶数N。 （2-1） 由式(2-2)确定。(2-2)由式(2-3)计算s左半平面的N个极点。k=1,2,...,N(2-3)由式(2-4)确定滤波器的系统函数H(s)。(2-4)2.2.2模拟低通滤波器转变为模拟带通滤波器的设计步骤：(1)确定模拟带通滤波器的技术指标，即：带通上限频率，带通下限频率；通带中心频率，通带宽度；通带最大衰减为，阻带最小衰减为原型低通到带通的变换为式(2-5) (2-5)(3)直接将低通转换成带通。2.2.3脉冲响应不变法的原理及特点假设模拟滤波器的系统函数为H(s)，模拟频率为，频率响应为，单位脉冲响应为h(t)；数字滤波器的系统函数为H(z)，数字频率为，频率响应为，单位取样响应为h(n)。设计步骤如下：将数字滤波器的技术指标转换为模拟滤波器的技术指标。利用模拟频率和数字频率的关系如式(2-6) (2-6)将数字滤波器的频率指标{}转换为模拟滤波器的频率指标{}。设计通带截频{}、通带衰减、阻带截频{}、阻带衰减的模拟滤波器。利用脉冲响应不变法将模拟滤波器的H(s)转换为数字滤波器的H(z)。脉冲响应不变法设计流程如图2.3所示。图2.3脉冲响应不变法设计过程脉冲响应不变法的优、缺点：脉冲响应不变法使得数字滤波器的单位冲激响应能完全模仿模拟滤波器的单位冲激响应，时域逼近良好，而且模拟角频率和数字角频率之间呈线性关系。该方法最大的缺点是有频率响应的混叠效应，所以只适用于限带的模拟滤波器(例如，衰减特性很好的低通或带通滤波器)，而且阻带衰减越快，混叠效应越小。2.2.4双线性变换法的原理及特点双线性变换法的基本思想是，将模拟滤波器的H(s)转换为数字滤波器的H(z)时，不是直接从s域到z域，而是先将非带限的H(s)映射为带限的H(s')，再通过脉冲响应不变法将s'域映射到z域，即H(s)—>H(s')—>H(z)。从频域来看模拟角频率与数字角频率的关系需通过'建立，即—>'—>。设计步骤如下：由式(2-7)将数字滤波器的频率指标{}转换为模拟滤波器的频率指标{}。 (2-7)设计通带截频{}、通带衰减、阻带截频{}、阻带衰减的模拟滤波器。利用双线性变换法将模拟滤波器的H(s)转换为数字滤波器的H(z)。遵循公式如式(2-8)。 (2-8)双线性变换法的优、缺点：双线性变换最突出的优点是避免了频率响应的混叠失真，缺点是频率响应的非线性失真，模拟角频率和数字角频率之间的关系如式(2-7)在零频率附近与之间的关系近似于线性，随着的增加，与之间的关系出现严重非线性，使数字滤波器频率响应不能保真地模仿模拟滤波器频率响应。双线性变换法的非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，否则变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有较大畸变。2.3设计例题要求通带范围为0.25πrad≤ω≤0.45πrad，通带最大衰减为3dB，阻带范围为0≤ω≤0.15πrad和0.55πrad≤ω≤πrad，阻带最小衰减为40dB。1数字带通滤波器技术指标通带上截止频率ωu=0.45πrad，通带下截止频率ωl=0.25πrad，阻带上截止频率ωs2=0.55πrad，阻带下截止频率ωs1=0.15πrad通带内最大衰减αp=3dB，阻带内最小衰减αs=40dB。2模拟带通滤波器技术指标为了计算简单，可设T=1，则有，，（通带中心频率）（带宽）将以上边界频率对带宽B归一化，得到，，，，3模拟归一化低通滤波器技术指标归一化阻带截止频率：归一化通带截止频率：λp=1，αp=3dB,αs=40dB4设计模拟低通滤波器，，取N=7查表得到归一化低通传输函数G(p),5模拟低通转换成模拟带通将归一化模拟低通转换成模拟带通6数字带通滤波器这里仅通过双线性变换法将Ha(s)转换成数字带通滤波器H(z)。下面将(5)、(6)两步合成一步计算：将上式代入上面的转换公式，得将上面的p等式代入G(p)中，得7Matlab仿真利用MATLAB设计IIR数字滤波器有多种方法，典型设计法原理清晰，完全设计法程序简单，最优设计法可满足特殊指标要求，工具设计法形象直观。灵活运用这些方法，在工程上有利于提高设计效率。此次仿真严格按照IIR滤波器设计方法进行的依次为巴特沃斯原型模拟滤波器—>转换的带通模拟滤波器—>脉冲响应不变法设计—>双线性变换法设计。程序见附录，效果如图2.4所示。图2.4IIR的Matlab仿真图由以上运行结果可看出，此数字滤带通滤波器各项技术指标均实现。第三部分：线性相位FIR数字滤波器的基本特性、设计方法（窗函数法和频率取样法）及两种设计方法的优缺点比较。3、FIR数字滤波器设计方法概述3.1FIR数字滤波器设计FIR滤波器的设计是建立在对期望滤波器频率特性的某种近似基础之上的目前有许多方法可以设计FIR滤波器，比如窗函数设计法、频率取样法等。其中窗函数设计法是滤波器设计的主要方法之一，由于运算简便，物理意义直观，已成为工程实际中应用最广泛的方法，常见的窗函数有：矩形窗、三角形窗、布莱克曼窗、切比雪夫窗等。3.2窗函数法设计数字滤波器 窗函数法的基本思想是用一具有有限长度样值响应、并具有线性相位的系统函数逼近理想滤波器的系统函数。就是根据给定的滤波器技术指标，选择滤波器的阶数N和合适的窗函数。用有限长度的窗口函数序列来截取一个无限长的序列获得一个有限长序列，即，并且要满足以下两个条件：(1)窗谱主瓣尽可能地窄，以获得较陡的过渡带；(2)尽量减小窗谱的最大旁瓣的相对幅度，也就是能量尽量集中于主瓣，使峰尖和纹波减小，就可增多阻带的衰减。这就给窗函数序列的形状和长度选择提出了严格的要求。常见的窗函数有：矩形窗、汉宁(Hanning)窗、布莱克曼窗、海明(Hamming)窗等。3.2.1用窗函数法设计FIR滤波器的步骤：根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N(或阶数M=N-1)，窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为Δw，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用Matlab中的窗函数求出窗函数wd(n)。根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)，如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：(3-1)采用离散傅里叶反变换（IDFT）即可求出。(3)用窗函数wd(n)将截断，并进行加权处理，得到(3-2)如果要求线性相位特性，则h(n)还必须满足：(3-3)根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。例如：要设计线性相位低通特性可选择h(n)=h(N-1-n)一类，而不能选h(n)=-h(N-1-n)一类。验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。3.2.2窗函数法的优、缺点：窗函数法的优点是简单，有闭合形式的公式可循，因而很实用。窗函数法是从时域出发，通过一定的窗函数截取有限长的单位脉冲响应来逼近理想单位脉冲响应；窗函数法的缺点是：(1)加窗后，会使频响产生一过渡带，其宽度正好等于窗的频响的主瓣宽度。(2)在处会出现肩峰，肩峰两侧形成起伏振荡，其振荡幅度取决于旁瓣的相对幅度，而振荡的多少则取决于旁瓣的多少。(3)会出现吉布斯(Gibbs)效应。(4)较为复杂时，不容易由反傅里叶变换求得。边界频率因为加窗的影响而不易控制。3.3频率取样法设计数字滤波器频率取样法是从频域出发，对理想的频响进行等间隔取样，以有限个频响采样去近似理想频响应。在实际使用时，为了设计线性相位的FIR滤波器，采样值H(k)要满足一定的约束条件，具有线性相位FIR滤波器，其单位采样响应函数h(n)是实序列，且满足h(n)=±(h-1-n)，由此得到幅频和EQ相频特性就是对H(k)的约束。3.3.1用频率取样法设计FIR滤波器的步骤：(1)根据所要求的滤波器类型，根据Ｎ是偶数还是奇数，指定，在阻带内，。(2)根据构成滤波器的和，并考察的指标是否满足要求。3.3.2频率抽样法的优、缺点：频率取样法设计滤波器的最大优点是直接从频率域进行设计，比较直观，也适合于设计具有任意幅度特性是滤波器，它十分适用于窄带滤波器的设计。频率取样法设计的缺点是由于频率抽样点的分布必须符合一定规律，在规定通、阻带截止频率方面不够灵活。比如当截止频率不是整数倍数时会产生较大逼近误差。3.4设计例题1用Hanning窗函数法和频率取样法设计一个线性相位FIR低通滤波器，并满足性能指标：通带边界频率p=0.5π，阻带边界频率s=0.66π，阻带衰减不小于40dB，通带波纹不大于3dB。2基本解题思路根据频域的采样定理，对给定理想滤波器的频率响应进行等间隔的采样：(k=0,1,……N-1)，把当作待设计的FIR数定滤波器的频率特性的采样值H(k)，即：(k=0,1,……N-l)，由H(k)通过IDFT可得有限长序列h(n)：(n=0,1,……N-l)，利用频域的N个采样点的值，通过下式求出滤波器的系统函数和频率响应：其中是一个内插函数，即：所以，在每个采样点上，频率响应严格与理想特性H(k)一致，在采样点之间，频率响应由各采样点的内插函数延伸叠加而形成，因而有一定的逼近误码差，误差大小与理想频率响应的曲线形状有关。理想特性平滑，则误差小，反之，误差大，并且在理想频率响应的不连续点处会产生肩峰和波纹。3Matlab仿真根据以上分析及所写设计程序，打开Matlab应用软件，输入设计程序(见附录)，最终程序运行结果为图3.1所示。此次仿真严格按照FIR滤波器设计方法进行的依次为单位脉冲响应—>窗函数法设计—>频率取样法设计—>优化法设计。程序，效果如图2.4所示。图3.1FIR的Matlab仿真图由以上运行结果可看出，此数字线性相位FIR低通滤波器各项技术指标均实现。第四部分：IIR和FIR数字滤波器的基本结构研究。4、IIR和FIR数字滤波器的基本结构概述4.1IIR数字滤波器基本结构IIR滤波器的单位脉冲响应h(n)为无限长序列，系统函数H(z)在有限z平面上存在极点，其运算结构的特点是含有反馈环路，即在结构上是递归型的。在给定滤波特性的情况下，IIR滤波器所用的项更少，此外，IIR滤波器还能实现窄带频响。IIR滤波器的实现结构并不唯一，同一系统函数(或差分方程)可以有各种不同的结构形式。基本结构主要有三种，即直接型、级联型和并联型。4.1.1IIR系统的直接实现形式对于IIR系统(4-1)相应的z变换可写成(4-2)式中W(z)、Y(z)对应的差分方程分别是 (4-3)信号流图如图4.1所示。图4.1IIR系统的直接实现由于数字系统的字长总是有限的，因此其系统精度总是有限的。每一个系统的量化误差及乘法器的舍入误差对输出都将有积累效应，以致输出误差偏大，这是直接实现形式的一个缺点。4.1.2IIR系统的级联实现形式将H(z)的分子分母多项式分成一阶或二阶多项式的连乘。考虑到H(z)若有复数极零点，则必为共轭成对出现，作物理实现时，其系数应为实数。因此将它们分解为二阶形式更为合理。若N≥M，N为偶数，则可将H(z)分成N/2个二阶多项式的连乘，若N为奇数，则子系统的数目应为(N+1)/2，其中包含一个一阶子系统。级联型结构中每一个一阶网络决定一个零点、一个极点，每一个二阶网络决定一对零点、一对极点。二阶子系统信号流图如图4.2所示。图4.2二阶子系统信号流图4.1.3IIR系统的并联实现形式将H(z)分解为各因式之和，则每个子系统有着共同的输入x(n)，而其输出之和便是系统的总输出y(n)。并联型结构中，每一个一阶网络决定一个实数极点，每一个二阶网络决定一对共轭极点。由于并联结构的每一个子系统都是独立的，不受其它子系统系数量化误差及乘法舍入误差的影响，因此是所述三种结构中误差最不敏感的结构形式。还有其它的一些数字滤波器，如全通滤波器，梳状滤波器，理想格型滤波器、简单整系数数字滤波器、采样率转换滤波器、平均滤波器、平滑滤波器和低阶低通差分滤波器等，由于不是本论文研究的主要内容，在此就不作介绍，只在后续章节中遇到时再作叙述。4.2FIR数字滤波器基本结构有限脉冲响应（简称FIR）系统的单位脉冲响应h(n)为有限长序列，系统函数H(z)在有限z平面上不存在极点，其运算结构中没有反馈支路，即没有环路。所以，有限脉冲响应滤波器可以设计成在整个频率范围内均可提供精确的线性相位，而且总是可以独立于滤波器系数保持有限输入有限输出（BIBO）稳定，因此在很多领域，这样的滤波器是首选的。FIR滤波器的实现形式通常有以下几种：直接型，级联型，多相实现和线性相位FIR滤波器结构。4.2.1FIR系统的直接实现形式直接型结构的输入输出关系如下：(4-4)通常在这种结构中需要N+1个乘法器和N个两输入加法器来实现。其结构图可以表示为图4.3直接型结构4.2.2FIR系统的级联实现形式级联型结构的输入输出关系如下：(4-5)高阶FIR传输函数可以由一阶或二阶传输函数级联实现，它是通过对式(4-5)进行因式分解得到的。(4-6)其中，当N为偶数时，K=N/2；当N是奇数时，K=(N+1)/2且β=0。由于级联形式是规范型结构，所以需要用N个两输入的加法器和N+1个乘法器来实现N阶有限脉冲响应传递函数。图4.4级联型结构有限脉冲响应滤波器的另一种实现是基于传输函数的多相位分解所得到的并联结构。一般情况下，L支N阶多相分解的传输函数具有以下形式：(4-7)式中(4-8)注意，在有限脉冲响应传输函数的多相实现中，子滤波器也是有限脉冲响应滤波器，并且子滤波器可以用直接或级联型方法实现。第五部分：最后比较一下IIR和FIR数字滤波器的主要优缺点，并对整篇课程论文内容进行概括性总结。5、数字滤波器设计方法总结5.1IIR数字滤波器与FIR数字滤波器比较A.IIR数字滤波器的主要优点是：(1)设计方法简单。通常只要将技术指标代入设计方程组就可以设计出原型滤波器，然后再利用相应的变换公式求得所需要的滤波器系统函数的系数。(2)在满足一定技术要求和幅频响应的情况下，IIR数字滤波器设计成为具有递归运算的环节。所以它的阶次一般比FIR数字滤波器低，所用的存储单元少，滤波器体积也小。B.IIR数字滤波器的主要缺点是：(1)只能设计出有限频段的低、高、带通和带阻等选频滤波器。除幅频特性可以满足技术要求外，它们的相频特性往往是非线性的，这就会使信号产生失真。(2)由于IIR数字滤波器采用了递归型结构，系统存在极点，因此设计系统函数时，必须把所有的极点放在单位圆内，否则系统不稳定。而且有限字长效应所带来的运算误差，可能会使得系统产生寄生振荡。C.FIR数字滤波器的主要优点是：(1)可以设计出具有线性相位的FIR数字滤波器，从而保证信号在传输过程中没有失真。(2)由于FIR数字滤波器没有递归运算，因此不论在理论还是实际应用中，都不会因为有限字长效应所带来的运算误差使得系统不稳定。(3)FIR数字滤波器可以采用快速傅里叶变换实现快速卷积运算，在相同阶数的条件下运算速度快。D.FIR数字滤波器的主要缺点是：虽然可以采用加窗方法或频率采样等简单方法设计FIR数字滤波器，但往往在过渡带上和阻带衰减上难以满足要求，因此不得不多次迭代或者计算机辅助设计，从而使得设计过程变得复杂。(2)在相同频率特性情况下，FIR数字滤波器阶次比较高，因而所需要的存储单元多，从而提高了硬件设计成本。5.2数字滤波器比较概括性总结本次课程设计首先给出了滤波的概念、分类及模拟滤波器设计，接着讨论了无限冲激响应和有限冲激响应数字滤波器的各种设计方法，重点是按照频域技术指标为依据的滤波器设计。对于无限冲激响应，介绍了冲激响应不变法、双线性映射法、IIR滤波器的频率变换设计法、IIR数字滤波器的直接设计法。对于有限冲激响应，介绍了FIR滤波器窗函数设计法、FIR滤波器频率采样设计法。另外通过课程设计我进一步见识到了Matlab功能的强大，内含各种丰富的函数可以让程序设计的简单，但是由于对Matlab不太熟悉，设计程序时还是走了些弯路，但后来还是通过查找资料找到了一个还算简单的设计，也算是得到了收获，我也会在以后的学习中注重对Matlab的学习与使用。6、参考文献：[1]陈后金．数字信号处理．2版[M]．北京：高等教育出版社。2008．11[2]孙强．运用MATLAB实现数字滤波器的设计[J]．电脑学习，2005(2)：32-33．[3]蔡建平．黄晓红，孙丽英，等．基于频率采样法的线性相位滤波器设计及MATI．AB仿真[J]．电气自动化设备，2006，26(7)：59—61[4]李勇，徐震．MATLAB辅助现代工程数字信号处理[M]．西安：西安电子科技大学出版社，2002．附录IIR仿真程序%BW型带通滤波器的指标Wp1=0.25*pi;Wp2=0.45*pi;Ws1=0.15*pi;Ws2=0.55*pi;B=Wp2-Wp1;W02=Wp2*Wp1;W0=sqrt(W02);%确定原型低通滤波器的指标Wp、Ws Wp11=(Wp1*Wp1-W02)/B/Wp1;Wp22=(Wp2*Wp2-W02)/B/Wp2;Ws11=(Ws1*Ws1-W02)/B/Ws1;Ws22=(Ws2*Ws2-W02)/B/Ws2;Wp=max(abs(Wp11),abs(Wp22));Ws=min(abs(Ws11),abs(Ws22));%设置通带最大和最小衰减Ap=1; As=32;%调用butter函数确定巴特沃斯滤波器阶数N,Wc[N,Wc]=buttord(Wp,Ws,Ap,As,'s'); %调用butter函数设计巴特沃斯滤波器[num,den]=butter(N,Wc,'s'); %绘出巴特沃斯模拟滤波器的幅频特性曲线W1=linspace(0,pi,400*pi); %指定一段频率值hf1=freqs(num,den,W1); %计算模拟滤波器的幅频响应subplot(2,2,1);plot(W1/pi,abs(hf1)/abs(hf1(1)));gridon;title('巴特沃斯原型模拟滤波器');xlabel('频率rad/s');ylabel('幅度');%将原型低通转为带通滤波器[numt,dent]=lp2bp(num,den,W0,B);%绘出带通滤波器的幅频特性曲线W1=linspace(0,pi,400*pi); %指定一段频率值hf2=freqs(numt,dent,W1);%计算模拟滤波器的幅频响应subplot(2,2,2);plot(W1/pi,abs(hf2)/abs(hf2(1))); gridon;title('转换的带通模拟滤波器');xlabel('频率rad/s');ylabel('幅度');%利用脉冲响应不变法设计数字带通滤波器T=1; %设置采样周期为1fs=4000; %采样频率为周期倒数wpz=[0.25,0.45];wsz=[0.15,0.55];wp1=wpz/T;ws1=wsz/T;[n1,wc1]=buttord(wp1,ws1,Ap,As,'s');%计算带通滤波器阶数N和3dB截止频率wc1,与前面BW不同 [b1,a1]=butter(n1,wc1,'s');%计算带通滤波器系统函数分子分母多项式系数向量b1,a1 [bz1,az1]=impinvar(b1,a1,fs); [nd1,wdc1]=buttord(wpz,wsz,Ap,As); [bd1,adz1]=butter(nd1,wdc1); hf3=freqz(bd1,adz1,W1)subplot(2,2,3);plot(W1/pi,abs(hf3)/abs(hf3(1))); %绘出巴特沃斯数字低通滤波器的幅频特性曲线gridon;title('脉冲响应不变法实现图');xlabel('频率rad/s');ylabel('幅度');%利用双线性不变法设计数字带通滤波器t=1;fs=4000; wpz=[0.25,0.45]; wsz=[0.15,0.55]; wp=2/t*tan(wpz/2);ws=2/t*tan(wsz/2); [n,wc]=buttord(wp,ws,Ap,As,'s');%计算带通滤波器阶数N和3dB截止频率wc [b,a]=butter(n,wc,'s');%计算带通滤波器系统函数分子分母多项式系数向量b,a [bz,az]=bilinear(b,a,fs); [nd,wdc]=buttord(wpz,wsz,Ap,As); [bd,adz]=butter(nd,wdc); hf4=freqz(bd,adz,W1) subplot(2,2,4) plot(W1/pi,abs(hf