IIR数字滤波器的MATLAB实现.doc

淮北师范大学2011届学士毕业论文 IIR数字滤波器的MATLAB实现IIR数字滤波器的MATLAB实现摘要 文中首先对数字滤波器作了简要的叙述。介绍了数字滤波器的分类情况，并对其工作原理进行了简要描述，指出了数字滤波器设计中应该满足的一定的技术指标。然后文中对IIR数字滤波器的设计基础作了一定的介绍。讨论了用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及优缺点。然后又介绍了用双线性变换法设计IIR数字滤波器的具体步骤。接着简要介绍了MATLAB语言、MATLAB工具箱函数以及用MATLAB语言设计巴特沃斯模拟低通滤波器的设计方法。最后重点介绍用MATLAB工具箱函数设计IIR数字滤波器。文中采用的设计方法是双线性变换法结合MATLAB工具箱函数设计出相应IIR数字滤波器，运行出幅频和相频特性图。关键词 MATLAB；双线性变换法；幅频相频特性图IIR digital filters MATLAB design Zhang CaiSchool of Physics and Electronic Information, Huaibei Normal University , 235000Abstract In this paper it first give a brief description of a digital filter. It describes the classification of digital filters, and briefly describes its working principle, pointing out that the digital filter design should meet a certain degree of technical indicators. Then in the text it describes the basis of the IIR digital filter design . Discussed the principle of designs the IIR digit filter, and the good and bad points with the pulse response not political reform and the bilinearity method of transformation. Then also introduced designs the IIR digit filters concrete step with the bilinearity method of transformation. Then briefly introduced the MATLAB language, the MATLAB toolbox function as well as uses MATLAB the language design Bahtwass simulation low pass filters design method. Finally introduced with emphasis designs the IIR digit filter with the MATLAB toolbox function . In the article uses the design method is the bilinearity transformation method combined with the MATLAB toolbox function to design the corresponding IIR digit filter, running out of amplitude frequency and phase frequency characteristics graph.Keywords MATLAB； Bilinearity method of transformation; amplitude frequency and phase frequency characteristics graph.目 次引言11 数字滤波器21.1 数字滤波器的分类21.2 数字滤波器的工作原理31.3 数字滤波器的设计指标32 IIR数字滤波器的设计方法52.1用脉冲相应不变法设计IIR数字滤波器52.2脉冲响应不变法优缺点：72.3用双线性变换法设计IIR数字滤波器82.4双线性变换法优缺点103 MATLAB语言简介123.1 MATLAB工具箱函数123.2 MATLAB功能和特点134 IIR数字滤波器144.1双线性数字滤波器设计步骤144.2用MATLAB设计模拟低通滤波器154.3 用MATLAB设计IIR滤波器实例17结 论22参考文献23致 谢2425淮北师范大学2011届学士毕业论文 IIR数字滤波器的MATLAB实现引 言随着信息时代与数字技术的发展，数字信号处理已逐渐发展成为当今极其重要的学科与技术领域之一。数字信号处理在通信、语音、图像、自动控制雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。在数字信号处理的基本方法中，通常会涉及到变换、滤波、频谱分析、调制解调和编码解码等处理。其中，滤波是应用非常广泛的一个环节，数字滤波器的相关理论也一直都是人们研究的重点之一。数字滤波器是数字信号处理的重要基础，在对信号的滤波、检测及参数的估计等信号应用中，数字滤波器是使用最为广泛的一种线性系统。数字滤波器根据其单位冲击响应函数的时域特性可分为两类：无限冲击响应(IIR)数字滤波器和有限冲击响应(FIR)数字滤波器。1 数字滤波器 所谓数字滤波器，是指输入输出均为数字信号，通过一定的运算关系，改变输入信号中所含频率成分的相对比例，或滤除某些频率成分的器件或程序1。数字滤波通过数值运算实现滤波，处理精度高，稳定，体积小，质量轻，灵活等优点。1.1 数字滤波器的分类 数字滤波器可以有很多种分类方法，但总体上可分为两大类：经典滤波器和现代滤波器。 经典滤波器的特点是其输入信号中有用的频率成分和希望滤除的频率成分各有不同的频带，通过一个合适的选频滤波器滤除干扰，得到纯净信号，达到滤波的目的。现代滤波器是根据随机信号的一些统计特性，在某种最佳准则下，最大限度地抑制干扰，同时最大限度地回复信号，从而达到最佳滤波的目的。 经典数字滤波器从滤波特性上分类，可以分成低通，高通，带通，和带阻等滤波器。 数字滤波器从实现的网络结构或者从单位脉冲回应长度分类，可以分成无限位脉冲回应（IIR, infinite impulse response）滤波器和有限长单位脉冲回应（FIR，finite impulse response）滤波器。它们的系统函数分别为: (1-1) (1-2) 1.2 数字滤波器的工作原理 数字滤波器是对数字信号实现滤波的线性时不变系统。数字滤波实质上是一种运算过程，实现对信号的运算处理。输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器本质上是一个完成特定运算的数字计算过程，也可以理解为是一台计算机。描述离散系统输出与输入关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。时域离散系统的频域特性:。其中、分别是数字滤波器的输出序列和输入序列的频域特性（或称为频谱特性）,是数字滤波器的单位取样响应的频谱，又称为数字滤波器的频域响应。输入序列的频谱经过滤波后,因此，只要按照输入信号频谱的特点和处理信号的目的， 适当选择，使得滤波后的满足设计的要求，这就是数字滤波器的滤波原理。数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应(IIR)数字滤波器和有限长冲激响应(FIR)数字滤波器。IIR 数字滤波器的特征是，具有无限持续时间冲激响应，需要用递归模型来实现，其差分方程为： 系统函数为：设计IIR滤波器的任务就是寻求一个物理上可实现的系统函数H(z)，使其频率响应H(z)满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止频率、通带衰减系数和阻带衰减系数。1.3 数字滤波器的设计指标 设数字滤波器的传输函数用式（1-3）表示： (1-3) 式中，为幅频特性，为相频特性2。幅频特性表示信号通过滤波器后各频率成分的衰减情况，相频特性则反映各频率成分通过滤波器后在时间上的延时情况。通常，选频滤波器的指标要求都以幅频特性给出，对相频特性不作要求，如果需要对输出波形有严格要求，如语音合成、波形传输等，则要求设计线性相位数字滤波器。数字滤波器的参数指标是、和。和分别称为通带截止频率和阻带截止频率。通带和阻带内允许的衰减一般用分贝数表示，通带内允许的最大衰减用表示，阻带内允许的最小衰减用表示，和分别定义为（1-4），（1-5）： (1-4) (1-5) 式中均假定已被归一化为1.则上式可转化为： 当幅度下降到时，标记，此时，称为3dB通带截止频率。 、和 统称为边界频率，它们是滤波器设计中所涉及到的重要参数。对其他类型的滤波器，（1-4）式和（1-5）中的 应改成 , 为滤波器通带中心频率1。2 IIR数字滤波器的设计方法IIR数字滤波器是一种离散时间系统，其系统函数为： 假设M N，当M N 时,系统函数可以看作一个IIR的子系统和一个(M-N)的FIR子系统的级联。IIR数字滤波器的设计实际上是求解滤波器的系数 和 ，它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。如果在S平面上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。2.1用脉冲相应不变法设计IIR数字滤波器利用模拟滤波器来设计数字滤波器，也就是使数字滤波器能模仿模拟滤波器的特性，这种模仿可以从不同的角度出发。脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应ha(t)，即将ha(t)进行等间隔采样，使h(n)正好等于ha(t)的采样值，满足h(n)=ha(nT)，式中,T是采样周期。如果令Ha(s)是ha(t)的拉普拉斯变换，H(z)为h(n)的Z变换，利用采样序列的Z变换与模拟信号的拉普拉斯变换的关系得（2-1）式： (2-1) 由图1可看出，脉冲响应不变法将模拟滤波器的S平面变换成数字滤波器的 Z 面，这个从s到z的变换是从S平面变换到Z平面的标准变换关系式。 0 S平面 z平面 图1脉冲响应不变法的映射关系由（2-1）式，数字滤波器的频率响应和模拟滤波器的频率响应间的关系为： (2-2) 这就是说，数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。正如采样定理所讨论的，只有当模拟滤波器的频率响应是限带的，且带限于折叠频率以内时，即： (2-3) 才能使数字滤波器的频率响应在折叠频率以内重现模拟滤波器的频率响应，而不产生混叠失真，即： (2-4)但是，任何一个实际的模拟滤波器频率响应都不是严格限带的，变换后就会产生周期延拓分量的频谱交叠，即产生频率响应的混叠失真，如图2所示。这时数字滤波器的频响就不同于原模拟滤波器的频响，而带有一定的失真。当模拟滤波器的频率响应在折叠频率以上处衰减越大、越快时，变换后频率响应混叠失真就越小。这时，采用脉冲响应不变法设计的数字滤波器才能得到良好的效果。 图2脉冲响应不变法中的频响混叠现象对某一模拟滤波器的单位冲激响应ha(t)进行采样，采样频率为fs，若使fs增加，即令采样时间间隔（T=1/fs）减小，则系统频率响应各周期延拓分量之间相距更远，因而可减小频率响应的混叠效应.2.2脉冲响应不变法优缺点从以上讨论可以看出，脉冲响应不变法使得数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，也就是时域逼近良好，而且模拟频率和数字频率之间呈线性关系=T。因而，一个线性相位的模拟滤波器（例如贝塞尔滤波器）通过脉冲响应不变法得到的仍然是一个线性相位的数字滤波器。脉冲响应不变法的最大缺点是有频率响应的混叠效应。所以，脉冲响应不变法只适用于限带的模拟滤波器(例如，衰减特性很好的低通或带通滤波器)，而且高频衰减越快，混叠效应越小。至于高通和带阻滤波器，由于它们在高频部分不衰减，因此将完全混淆在低频响应中。如果要对高通和带阻滤波器采用脉冲响应不变法，就必须先对高通和带阻滤波器加一保护滤波器，滤掉高于折叠频率以上的频率，然后再使用脉冲响应不变法转换为数字滤波器。当然这样会进一步增加设计复杂性和滤波器的阶数。2.3用双线性变换法设计IIR数字滤波器脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-/T/T之间，再用z=esT转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-/T/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图3所示。 0 0 -1 0 1 S平面 S1平面 Z平面 图3双线性变换的映射关系为了将S平面的整个虚轴j压缩到S1平面j1轴上的-/T到/T段上，可以通过以下的正切变换实现 （2-5） 式中,T仍是采样间隔。当1由-/T经过0变化到/T时，由-经过0变化到+，也即映射了整个j轴。将式（2-5）写成 将此关系解析延拓到整个S平面和S1平面，令，则得 再将S1平面通过以下标准变换关系映射到Z平面从而得到S平面和Z平面的单值映射关系为： （2-6） （2-7） 式（2-6）与式（2-7）是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换。式（2-5）与式（2-6）的双线性变换符合映射变换应满足的两点要求。首先,把，可得式（2-8） （2-8） 即S平面的虚轴映射到Z平面的单位圆。其次，将 代入式（2-8），得 因此 由此看出，当0时，|z|0时，|z|1。也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。2.4双线性变换法优缺点双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个j轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。这个关系如式（2-8）所示，重写如下： 上式表明，S平面上与Z平面的成非线性的正切关系，如图4所示。由图4看出，在零频率附近，模拟角频率与数字频率之间的变换关系接近于线性关系；但当进一步增加时，增长得越来越慢，最后当时，终止在折叠频率=处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。 图4双线性变换法的频率变换关系但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（2-8）及图4所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图5所示： 0 0 0 0 0 0 图5双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸变来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。3 MATLAB语言简介MATLAB是Mathworks 公司开发的，目前国际上最流行、应用最广泛的科学与工程软件。MATLAB即Matrix+Laboratory，又称为” 矩阵实验室”，其强项就是强大的矩阵计算。MATLAB是MATLAB产品家族的基础，它提供的数学算法，例如矩阵运算、数值分析算法。MATLAB集成了2D和3D图形功能，以完成相应数值可视化的工作，并且提供了一种交互式的高级编程语言M语言，利用M语言可以通过编写脚本或者函数文件实现用户自己的算法。MATLABCompiler是一种编译工具，它能够将那些利用MATLAB提供的编程语言M语言编写的函数文件编译生成函数库、可执行文件COM组件等，从而可以扩展MATLAB功能，使MATLAB能够同其他高级编程，例如C/C+语言进行混合应用，取长补短，以提高程序的运行效率，丰富程序开发的手段。利用M语言还开发了相应的MATLAB专业工具箱函数供用户直接使用。这些工具箱应用的算法是开放的、可扩展的，用户不仅可以查看其中的算法，还可以针对一些算法进行修改，甚至允许开发自己的算法扩充工具箱的功能。目前MATLAB产品的工具箱函数有40多个，分别涵盖了数据获取、科学计算、控制系统与分析、数字信号处理、数字图像处理、金融财务分析以及生物遗传工程等专业领域3。3.1 MATLAB工具箱函数MATLAB工具箱实际上是用MATLAB的基本语句编程的各种子程序集，用于解决某一方面的专门问题或实现某一类的新算法。MATLAB的工具箱可以增减，不同的工具箱给不同的领域的用户提供了丰富、强大的功能。MATLAB主要有以下工具箱。 控制系统工具箱，系统识别工具箱，信号辨别工具箱，神经网络工具箱，模糊逻辑控制工具箱，小波工具箱，模型预测控制工具箱，通信工具箱，图形处理工具箱，频域系统识别工具箱，优化工具箱，偏微分方程工具箱，财政金融工具箱，系统工具箱等。3.2 MATLAB功能和特点MATLAB集科学与工程计算、图形可视化、图像处理、多媒体处理于一体，并提供了Windows图形界面设计方法。MATLAB语言有以下特点。1. 功能强大MATLAB语言的功能强大体现在以下几个方面。（1） 运算功能强大。MATLAB是以复数矩阵为基本编程单元的程序设计语言，其强大的运算功能使其成为世界顶尖的数学应用软件之一。 MATLAB的数值运算要素不是单个数据，而是矩阵，每个变量都代表一个矩阵，矩阵有mn个元素，每个元素都可视为复数，所有的运算包括加、减、乘、除、函数运算等都对矩阵和复数有效；另外，通过MATLAB的符号工具箱，可以解决在数学、应用科学和工程计算领域中常常遇到的符号计算问题。（2） 功能丰富的工具箱。大量针对各专业的功能工具箱的提供，是MATLAB适用于不同领域。（3） 文字处理功能强大。MATLAB的Notebook为用户提供了强大的文字处理功能，允许用户从Word访问MATLAB的数值计算和可视化结果。通过使用MATLAB的Notebook，用户可以创建MATLAB的程序文档、技术报告、注释文档、手册或教科书。2. 人机界面友好，编程效率高MATLAB的语言规划与笔算式相似，矩阵的行列数无须定义，MATLAB的命令表达方式与标准的数学表达式非常相近，易写、易读并易于在科技人员之间交流。MATLAB是以解释方式工作的，即它对每条语句解释后立即执行，输入算式无须编译立即得到结果，若有错误也立即做出反应，便于编程者立即改正。这些都大大减轻了编程和调试的工作量，提高了编程效率。3. 强大而智能化的作图功能 MATLAB可以方便地将工程计算的结果可视化，使原始数据的关系更加清晰明了，并揭示了数据间的内在联系。MATLAB能根据输入数据自动确定最佳坐标，可规定多种坐标，可设置不同颜色、线性、视角等，并能绘制三维坐标中的曲线和曲面3。4 IIR数字滤波器IIR滤波器设计方法有间接法和直接发，间接法是借助于模拟滤波器的设计方法进行的。其设计步骤是：先设计过渡模拟滤波器得到函数，然后将按某种方法转换成数字滤波器的系统函数。根据以上IIR数字滤波器设计方法，下面运用MATLAB语言设计IIR低通，高通，带通和带阻滤波器。首先我们先介绍模拟低通滤波器的设计，这是因为低通滤波器的设计师设计其他滤波器的基础。4.1双线性数字滤波器设计步骤对于数字高通、带通和带阻滤波器的设计，通用的方法为双线性变换。可以借助于模拟滤波器的频率变换设计一个所需类型的过度模拟滤波器，再通过双线性变换将其转换成所需类型的数字滤波器。流程图如图6所示4；开 始读入所需类型的数字滤波器技术指标将所需类型数字滤波器的边界频率转换成相应类型模拟滤波器的边界频率 将相应类型模拟滤波器技术指标转换成模拟低通滤波器技术指标设计模拟低通滤波器通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器采用双变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器结 束 图6 双线性数字滤波器设计步骤流程图4.2用MATLAB设计模拟低通滤波器这里我们用MATLAB工具箱函数设计低通巴特沃斯滤波器。MATLAB信号处理工具箱函数buttap,buttord,butter是巴特沃斯滤波器设计函数。其5种调用格式如下： 1） Z, P, k=buttap( N ) 该格式用于计算N阶巴特沃斯归一化（3 dB截止频率c=1）模拟低通原型滤波器系数函数的零，极点和增益2。返回长度为N的列向量和，分别给出个零点和极点的位置，表示滤波器增益。2） N, Wc= buttord（wp， ws， Rp， As）该格式用于计算巴特沃斯数字滤波器的阶数N和3 dB截止频率wc。调用参数wp和ws分别为数字滤波器的通带边界和阻带边界频率的归一化值。3） N, wc=buttord(wp, ws, Rp, As, s)该格式用于计算巴特沃斯模拟滤波器的阶数N和3 dB截止频率wc。Wp、ws和wc是实际模拟角频率（rad/s）.4） B, A=buttord(N, wc, ftype)计算N阶巴特沃斯数字滤波器系统函数分子和分母多项式系数向量B和A。5） B, A=butter(N, wc, ftype, s)计算巴特沃斯模拟滤波器系统函数的分子和分母多项式的系数向量B和A。 Ftype=high时，设计3dB截止频率为wc的高通滤波器。缺省ftype是默认设计低通滤波器。Ftype=stop是，设计3dB截止频率为wc的高通滤波器。此时wc为二元向量wcl, wcu,wcl和wcu分别为带阻滤波器的通带3dB下截止频率和上截止频率。缺省ftype时设计带通滤波器，通带为频率区间wclwwcu。例：已知通带截止频率f =5 kHz，通带最大衰减a =2 dB，阻带截止频率f =12 kHz，阻带最小衰减a =30 dB，按照以上技术指标设计巴特沃斯低通滤波器。程序：wp=2*pi*5000;ws=2*pi*12000;Rp=2;As=30; %设置滤波器参数N, wc=buttord(wp,ws,Rp,As,s); %计算滤波器阶数N和3dB截止频率wcB,A=butter(N, wc,s); %计算滤波器系统函数分子分母多项式系数K=0:511;fk=0:14000/512:14000;wk=2*pi*fk; Hk=freqs(B, A, wk);subplot(1, 1, 1);plot(fk/1000, 20*log10(abs(Hk);grid onXlabel(频率(kHz);ylabel(幅度（dB）);axis(0, 14 -40 ,5)运行结果： N=5，wc=3.7792e+004，B=7.7094e+022 A=1 1.2230e+005 7.4785e+009 2.8263e+014 6.6014e+018 7.7094e+022其幅频响应曲线如图7所示：图7 巴特沃斯模拟低通滤波器的幅频响应曲线设计模拟高通、带通和带阻滤波器的一般过程是： (1) 通过频率变化公式，先将希望设计的滤波器指标转换为相应的低通滤波器指标(2) 设计相应的低通滤波器系统函数Q（p）.(3) 对Q（p）进行频率变换，得到希望设计的滤波器系统函数Hd（s）。 4.3 用MATLAB设计IIR滤波器实例1. 设计低通数字滤波器，频率低于0.2rad时，容许幅度误差在1dB以内；在频率0.3到之间的阻带衰减大于15dB。指定模拟滤波器采用巴特沃斯低通滤波器。试双线性变换法设计数字滤波器。程序：T=1;Fs=1/T; wpz=0.2;wsz=0.3; rp=1;rs=15; %设置设计要求参数wp=2*tan(wpz*pi/2);ws=2*tan(wsz*pi/2); %预畸变校正转换指标N,wc=buttord(wp,ws,rp,rs,s); %设计过渡模拟滤波器B,A=butter(N,wc,s); %计算过渡模拟滤波器系统函数的分子分母多项式系数Bz,Az=bilinear(B,A,Fs); %用双线性变换法转换成数字滤波器Freqz(Bz,Az); %作出H(z)的幅频相频图其幅频相频特性图如图8所示：图8 数字低通滤波器幅频相频特性曲线图2. 设计一个数字高通滤波器，通带截止频率fp=400hz，阻带截止频率fs=300hz；通带衰减不大于1dB，阻带衰减不小于15dB。采样频率Fs=1000hz 。程序如下：fp = 400 ; fs = 300; rp = 1; %通带最大衰减Rp=1dBrs = 20; %阻带最小衰减Rs=20dBwp =fp*2*pi; %把数字域滤波器特征换成模拟滤波器特征ws =fs*2*pi;Fs=1000;T=1/Fs %归一化数字频率Wp=wp/(Fs);Ws=ws/(Fs); %频率预畸变：数字域频率wp2=2*tan(Wp/2)/T; %预畸变求滤波器通带临界频率ws2=2*tan(Ws/2)/T; %预畸变求滤波器阻带临界频率%设计模拟滤波器N,Wn = buttord(wp2,ws2,rp,rs,s) z,p,k=buttap(N); %创建Buttord低通滤波器原型Bap,Aap=zp2tf(z,p,k); %由零极点转换为传递函数的形式 Bbs,Abs=lp2hp(Bap,Aap,Wn); %模拟低通变高通%用双线性不变法变换成数字滤波器Bbz,Abz=bilinear(Bbs,Abs,Fs); %双线性变换freqz(Bbz,Abz,512,Fs);其幅频相频特性图如图9所示：图9 数字高通滤波器幅频相频特性曲线图3. 设计一个带通滤波器，采样频率Fs=2000Hz，fpl=300 Hz,fpu=400 Hz,fsl=200 Hz, fsu= 500Hz , rp=3 dB , rs =18 dB。程序如下： fpl=300; fpu=400;fp=fpl fpu; fsl=200 ;fsu=500; fs=fsl fsu; rp=3; rs=18; Fs=2000; wp=fp*2*pi/Fs ; ws=fs*2*pi/Fs; wap=2*Fs*tan(wp./2); was=2*Fs*tan(ws./2); n,wn=buttord(wap,was,rp,rs,s); z,p,k=buttap(n); bp,ap=zp2tf(z,p,k); bw=wap(2)-wap(1) w0=sqrt(wap(1)*wap(2) bs,as=lp2bp(bp,ap,w0,bw) h1,w1=freqs(bp,ap); w2=0:Fs/2-1*2*pi; h2