用脉冲响应不变法设计数字滤波器

皖西学院 《数字信号处理》课程报告     目   用脉冲响应不变法设计数字滤波器 学    院       信息学院      专    业       通信工程专业      班    级         （*** ）班      学生姓名       陈*    孙**      指导教师           吴**          二 0一二年十二月 前言     《数字信号处理》课程设计是在学生完成数字信号处理和MATLAB结合后的基本实验后开设的，本课程设计的目的是为了让学生综合数字信号处理和MATLAB并实现一个较为完整的小型滤波系统。这一点与验证性的基本实验有本质性的区别。开设课程设计环节的主要目的是通过系统设计、软件仿真、程序安排与调试、写实习报告等步骤，使学生初步掌握工程设计的具体步骤和，提高分析问题和解决问题的能力，提高实际应用水平。     IIR数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配，所以IIR滤波器的设计可以采用在模拟滤波器设计的基础上进一步变换的方法。其设计方法主要有间接设计法、直接设计法和最大平滑滤波器设计方法。间接法是借助于模拟滤波器的设计方法进行的。其设计步骤是：先设计过度模拟滤波器得到系统函数，然后将其按某种方法转换成数字滤波器的系统函数。这是因为模拟滤波器的设计方法已经成熟，不仅有完整的设计公式，还有完善的图表和曲线供查阅；另外还有一些优良的滤波器可供我们使用。直接法直接在频域或者时域中设计数字滤波器，由于要解联立方程，设计时需要计算机辅助设计。FIR数字滤波器的单位脉冲响应应是有限长序列。它的设计问题实质上是确定能满足要求的转移序列或脉冲响应的常数问题，它不能采用间接法，设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。 目 录 第1章  绪论    2 1.1课程设计的目的及意义    2 1.2课程设计题目描述及要求    3 1.3数字滤波器的概述    3 1.4数字滤波器的分类    3 1.5数字滤波器的技术指标    4 1.6数字滤波器的设计原理    5 第2章  MATLAB介绍    5 2.1 MATLAB的简介    5 2.2 MATLAB的优势和特点    6 第3章 IIR数字滤波器的设计    6 3.1 IIR数字滤波器的设计概述    6 3.2 IIR数字滤波器的设计思想：    7 3.3脉冲响应不变法设计数字滤波器    7 3.4 巴特沃斯滤波器的设计原理    11 第4章  利用脉冲响应不变法设计数字滤波器的过程    16 4.1课程设计的解题思路及过程    16 4.2    MATLAB程序及仿真    17 第5章 总结    20 参考文献    21 第1章  绪论 1.1课程设计的目的及意义 电子信息工程专业的培养目标是具备电子技术的基本理论和应用技术，能从事电子、信息、通信、电信等领域的工作，具有高素质、宽口径、创新晋升的专业人才。对本专业学生的培养要进行工程素质培养、拓宽专业口径、注重基础和发展潜力。特别是培养学生的创新能力，以实现技术为主线多进行实验技能的培养。学生通过数字信号处理课程设计这一重要环节，可以将本专业的主干课程-数字信号处理从理论学习到实践应用，对数字信号处理技术有较深的了解，进一步增强学生动手能力和适应实际工作的能力。 1.2课程设计题目描述及要求 设计题目：脉冲响应不变法设计数字滤波器。 设计技术指标：p=0.2，p≤2dB；阻带截止频率s=0.6 ，阻带最小衰减s≥15dB。 设计要求：设计模拟巴特沃斯滤波器并脉冲响应不变法转化成数字滤波器。 1.3数字滤波器的概述     数字滤波器一词出现在60年代中期。由于电子计算机技术和大规模集成电路的发展，数字滤波器已可用计算机软件实现，也可用大规模集成数字硬件实时实现。 所谓数字滤波器，是指输入、输出均为数字信号，通过数值运算处理改变输入信号所含频率成分的相对比例，或者滤除某些频率成分的数字器件或程序。因此，数字滤波器的概念和模拟滤波器相同，只是信号的形式和实现滤波方式不同。应用数字滤波器处理模拟信号时，首先须对输入模拟信号进行限带、抽样和模数转换。数字滤波器输入信号的抽样率应大于被处理信号带宽的两倍，其频率响应具有以抽样频率为间隔的周期重复特性，且以折叠频率即1/2抽样频率点呈镜像对称。为得到模拟信号，数字滤波器处理的输出数字信号须经数模转换、平滑。数字滤波器具有高精度、高可靠性、体积小、重量轻、不存在阻抗匹配问题、可程控改变特性或复用、便于集成等优点，可以实现模拟滤波器无法实现的特殊滤波功能。数字滤波器在语言信号处理、图像信号处理、医学生物信号处理以及其他应用领域都得到了广泛应用。 1.4数字滤波器的分类     按照不同的分类方法，数字滤波器有许多种类，按照不同的特性具有不同的分类，具体分类如下： 按功能分：低通、高通、带通、带阻、全通滤波器。如果我们设计这些理想滤波器是不可能实现的，因为它们的单位脉冲响应是非因果且无限长的，我们只能按照某些准则设计滤波器，使之在误差容限内逼近理想滤波器，理想滤波器可作为逼近的。 按实现的网络结构或单位抽样响应长度分类：无限脉冲响应滤波器（IIR滤波器）、有限脉冲响应滤波器（FIR滤波器）。它们的系统函数为：                   按滤波器对信号的处理作用可分为：选频滤波器和其他滤波器。上述低通、高通、带通、带阻滤波器均属于选频滤波器，其他滤波器有微分器、希尔伯特变换器、频谱校正等滤波器。   另外，它还可以被分为线性与非线性、因果与非因果等。   其中，线性时不变的数字滤波器是最基本的类型；而由于数字系统可以对延时器加以利用，因此可以引入一定程度的非因果性，获得比传统的因果滤波器更灵活强大的特性；相对于IIR滤波器，FIR滤波器有着易于实现和系统绝对稳定的优势，因此得到广泛的应用；对于时变系统滤波器的研究则导致了以卡尔曼滤波为代表的自适应滤波理论。 1.5数字滤波器的技术指标 常用的数字滤波器一般属于选频滤波器。假设数字滤波器的频率响应函数H(ejω)用下式表示： 式中，|H(e)|称为幅频特性函数;  θ(ω)称为相频特性函数。     　图1.2表示低通滤波器的幅频特性，ωp和ωs分别称为通带边界频率和阻带截止频率。通带频率范围为0≤|ω|≤ωp，在通带中要求(1－δ1)<|H(ejω)|≤1，阻带频率范围为ωs≤|ω|≤π，在阻带中要求|H(ejω)|≤δ2。从ωp到ωs称为过渡带，过渡带上的频响一般是单调下降的。通常，通带内和阻带内允许的衰减一般用分贝数表示，通带内允许的最大衰减用p表示，阻带内允许的最小衰减用s表示。p 越小,通带波纹越小，通带逼近误差就越小；s越大,阻带波纹越小，阻带逼近误差就越小；ωp与ωs间距越小, 过渡带就越窄。所以低通滤波器的设计指标完全由通带边界频率Ωp、通带最大衰减p 阻带边界频率Ωs和阻带最小衰减s确定。 1.6数字滤波器的设计原理 数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应(IIR)滤波器和有限长冲激响应(FIR)滤波器。IIR滤波器的特征是，具有无限持续时间冲激响应。这种滤波器一般需要用递归模型来实现，因而有时也称之为递归滤波器。FIR滤波器的冲激响应只能延续一定时间，在工程实际中可以采用递归的方式实现，也可以采用非递归的方式实现。数字滤波器的设计方法有多种，如双线性变换法、窗函数设计法、插值逼近法等等。随着MATLAB软件尤其是MATLAB的信号处理工作箱的不断完善，不仅数字滤波器的计算机辅助设计有了可能，而且还可以使设计达到最优化。 第2章  MATLAB介绍 2.1 MATLAB的简介     MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，并且吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。 2.2 MATLAB的优势和特点 （1）友好的工作平台和编程环境MATLAB由一系列工具组成。 （2）简单易用的程序语言。 （3）出色的图形处理功能 （4）应用广泛的模块集合工具箱。 （5）应用软件开发（包括用户界面）在开发环境中，使用户更方便地控制多个文件和图形窗口。 第3章 IIR数字滤波器的设计 3.1 IIR数字滤波器的设计概述 IIR滤波器系统函数的极点可以在单位圆内的任何位置，实现IIR滤波器的阶次较低，所用的存储单元较少，效率高，又由于IIR数字滤波器能够保留一些模拟滤波器的优良特性，因此应用很广。设计IIR数字滤波器的方法主要有基于脉冲响应不变法的IIR数字滤波器设计，基于双线性Z变换法的IIR数字滤波器设计，数字高通、带通及带阻IIR滤波器设计，基于MATLAB函数直接设计IIR数字滤波器。 3.2 IIR数字滤波器的设计思想： 目前，IIR数字滤波器设计最通用的方法是借助于模拟滤波器的设计方法。模拟滤波器设计已经有了一套相当成熟的方法，它不但有完整的设计公式，而且还有较为完整的图表供查询，因此，充分利用这些已有的资源将会给数字滤波器的设计带来很大方便，IIR数字滤波器的设计步骤是： (1)按一定规则将给出的数字滤波器的技术指标转换为模拟滤波器的技术指标。 (2)根据转换后的技术指标设计模拟低通滤波器H(s)。 (3)在按一定规则将H(s)转换为H(z)。 若所设计的数字滤波器是低通的，那么上述设计工作可以结束。若所设计的是高通、带通或者带阻滤波器，那么还有步骤(4)将高通、带通或者带阻数字滤波器的技术指标先转化为低通滤波器的技术指标，然后按上述步骤(2)设计出模拟低通滤波器H(s)，再由冲击响应不变法或双线性变换将H(s)转换为所需的H(z) 设计的关键问题：找出从到的转换方法，主要有冲激响应不变法、阶跃响应不变法、双线性变换法等。下面讨论脉冲响应不变法。 3.3脉冲响应不变法设计数字滤波器 设模拟滤波器的系统函数为，相应的单位冲激响应是，=LT［］。 LT［·］代表拉氏变换，对进行等间隔采样，采样间隔为T，得到，将作为数字滤波器的单位脉冲响应，那么数字滤波器的系统函数H(z)便是h(n)的Z变换。因此脉冲响应不变法是一种时域逼近方法，它使h(n)在采样点上等于h(t)。  但是，模拟滤波的设计结果是H(s)，所以下面基于脉冲响应不变法的思想，导出直接从H(s)到的转换公式。 设模拟滤波器H(s)只有单阶极点，且分母多项式的阶次高于分子多项式的阶次，将H(s)用部分分式表示：                                                     (3.1)   式中si为H(s)的单阶极点。将Ha(s)进行逆拉氏变换, 得到:                                                                                 (3.2)   式中, u(t)是单位阶跃函数。对ha(t)进行等间隔采样，采样间隔为T，得到：                                                       (3.3)   对上式进行Z变换，得到数字滤波器的系统函数H(z)，即                                                                                 (3.4)     对比(3.1)和(3.4)式，H(s)的极点映射到z平面的极点为，系数A不变。下面我们分析从模拟滤波器转换到数字滤波器，s平面和z平面之间的映射关系，从而找到这种转换方法的优缺点。这里以理想采样信号作为桥梁，推导其映射关系。 设的理想采样信号用表示，即                     = 对进行拉氏变换，得到：                                                                                 （3.5）  上式表示采样信号的拉氏变换与相应的序列的Z变换之间的映射关系可用下式表示：                                                           （3.6） 我们知道模拟信号h(t)的傅里叶变换Ha(jΩ)和其采样信号的傅里叶变换              之间的关系满足： 将s=jΩ代入上式，得 设s=,z=re                                                                                                                     （3.7）  由此可得 ①  r与的关系 =0 (s平面虚轴)                r=1 (z平面单位圆) <0 (s左半平面)                r<1 (z平面单位圆内部) >0 (s右半平面)                r>1 (z平面单位圆外部) ②  与的关系： = T =0 (s平面实轴)                =0 (z平面正实轴) = (s平面平行于实轴的直线)  =T(z平面始于原点，辐角为T的辐射线) ：从-/T～/T                ：从-～   (s平面为2/T的一个水平带)      (z平面辐角转了一周，覆盖整个z平面) 图3.1 脉冲响应不变法s平面和z平面之间的映射关系 存在的问题：混叠失真                                                             （3.8） 数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。 图3.2 脉冲响应不变法的频谱混叠现象示意图 假设H(ejΩT)没有频谱混叠现象，即满足                                   由(3.8)式得到：                           |<          (3.9) 上式说明，如果不考虑频谱混叠现象，用脉冲响应不变法设计的数字滤波器可以很好的重现原模拟滤波器的频响特性。但是，H()的幅度与采样间隔成反比，当T很小时， | H()|就会有太高的增益。为避免这一现象，令                         h 那么                         H                  (3.10) (3.10)式称为实用公式，此时                                     < 一般H(s)的极点是一个复数，且以共轭成对的形式出现，将一对复数共轭极点放在一起，可形成一个二阶基本节。如果模拟滤波器的二阶基本节的形式为：                     ，          极点为- 可以推导出相应的数字滤波器二阶基本节(只有实数乘法)的形式为： 如果模拟滤波器二阶基本节的形式为：                         ，          极点为-  相应的数字滤波器二阶基本节的形式为： 利用以上这些变换关系，可以简化设计，使实现结构中无复数乘法器。 3.4 巴特沃斯滤波器的设计原理 巴特沃斯低通滤波器的幅度平方函数|H(jΩ)|2用下式表示：                                                       （3.11） 式中，N称为滤波器的阶数。当Ω=0时，|Ha(jΩ)|=1；Ω=Ωc时，，Ωc是3 dB截止频率。在Ω=Ωc附近，随Ω加大，幅度迅速下降。 幅度特性与Ω和N的关系如图3.3所示。幅度下降的速度与阶数N有关，N愈大，通带愈平坦，过渡带愈窄，过渡带与阻带幅度下降的速度愈快,  总的频响特性与理想低通滤波器的误差愈小。 以s替换jΩ，将幅度平方函数|Ha(jΩ)|2写成s的函数：                                                             （3.12）                    复变量s=+jΩ，此式表明幅度平方函数有2N个极点，极点用下式表示：                                 （3.13） 式中，k=0，1，2，…，2N-1。2N个极点等间隔分布在半径为Ωc的圆上(该圆成为巴特沃斯圆)，间隔是π/N rad。例如N=3，极点间隔为π/3 rad，如图3.4所示。 图3.3 巴特沃斯低通滤波器幅度特性与N的关系 图3.4 三阶巴特沃斯滤波器极点分布图( 为形成因果稳定的滤波器，2N个极点中只取s平面左半平面的N个极点构成Ha(s)，而右半平面的N个极点构成Ha(－s)。Ha(s)的表达式为                                                           （3.14） 设N=3，极点有6个，它们分别为 取s平面左半平面的极点s0、s1、s2组成系统函数Ha(s)，即             由于不同的技术指标对应的边界频率和滤波器幅频特性不同，为使设计公式和图表统一，将频率归一化。巴特沃斯滤波器采用对3 dB截止频率Ωc归一化，归一化后的系统函数为                                                             （3.15）    令p=η+jλ=s/Ωc，λ=Ω/Ωc，λ称为归一化频率，p称为归一化复变量，这样巴特沃斯滤波器的归一化低通原型系统函数为                                                                                                       （3.16） 式中，pk=sk/Ωc,  为归一化极点，用下式表示：                                                             （3.17） 显然,                                                                (3.18)    这样，只要根据技术指标求出阶数N，按照(3.17)式求出N个极点，再按照(3.16)式得到归一化低通原型系统函数Ga(p)，如果给定Ωc，再去归一化，即将p=s/Ωc代入Ga(p)中(或由(3.18)式求出s=Ωp)，便得到期望设计的系统函数Ha(s)。 将极点表示式(3.17)代入(3.6)式，得到Ga(p)的分母是p的N阶多项式，用下式表示：                    (3.19) 归一化原型系统函数Ga(p)的系数bk，k=0，1，…，N－1，以及极点pk，可以由表3.5得到。另外，表中还给出了Ga(p)的因式分解形式中的各系数，这样只要求出阶数N，查表可得到Ga(p)及各极点, 而且可以选择级联型和直接型结构的系统函数表示形式，避免了因式分解运算工作。 表3.5 巴特沃斯归一化低通滤波器参数 由(3.13)式和(3.14)式可知,  只要求出巴特沃斯滤波器的阶数N和3 dB截止频率Ωc，就可以求出滤波器的系统函数Ha(s)。所以，巴特沃斯滤波器的设计实质上就是根据设计指标求阶数N和3 dB截止频率Ωc的过程。下面先介绍阶数N的确定方法。 阶数N的大小主要影响通带幅频特性的平坦程度和过渡带、阻带的幅度下降速度，它由技术指标Ωp、p、Ωs和s确定。将Ω=Ωp代入幅度平方函数(3.11)式中，再将幅度平方函数|Ha(jΩ)|代入，得到：                                                             (3.20) 将Ω=Ωs代入(3.4.11)式中，再将|Ha(jΩ)|2代入，得到：                                                             (3.21)    由(3.20)和(3.21)式得到：                                                                             令                                                               (3.22a)                                                               (3.22b)  则N由下式表示：                                                               (3.22c) 用上式求出的N可能有小数部分，应取大于或等于N的最小整数。关于3 dB截止频率Ωc，如果技术指标中没有给出，可以按照(3.20)式或(3.21)式求出。由(3.20)式得到：                                                               (3.23) 由(3.21)式得到：                                                               (3.24) 请注意，如果采用(3.23)式确定Ω，则通带指标刚好满足要求，阻带指标有富余； 如果采用(3.24)式确定Ωc，则阻带指标刚好满足要求，通带指标有富余。 总结以上，低通巴特沃斯滤波器的设计步骤如下：   (1)  根据技术指标Ωp、 p、Ωs和s，用(3.22)式求出滤波器的阶数N。   (2)  按照(3.17)式，求出归一化极点，将代入(3.16)式，得到归一化低通原型系统函数Ga(p)。也可以根据阶数N直接查表3.5得到和G(p)。 (3)  将Ga(p)去归一化。将p=s/Ωc代入Ga(p)，得到实际的滤波器系统函数                       这里Ωc为3 dB截止频率，如果技术指标没有给出Ωc，可以按照(3.23)式或(3.24)式求出。 第4章  利用脉冲响应不变法设计数字滤波器的过程 4.1课程设计的解题思路及过程 利用模拟巴特沃斯滤波器及脉冲响应不变法设计一数字滤波器，技术指标如下：p=0.2，p≤2dB；阻带截止频率s=0.6 ，阻带最小衰减s≥15dB。 （1）将数字滤波器的指标转换为相应的模拟滤波器指标。设采样周期为T，，，由(3.7)式得到：                           ， dB             ， dB (2）设计相应的模拟滤波器，得到模拟系统函数。 ①确定阶数N。           41.3223           3             3.39，取N=4 ②按照（3.17）式，其极点为                         ，，                         ， 按照（3.16）式，归一化低通原型系统函数为                               上式分母可展开成四阶多项式，或者将共轭极点放在一起，形成因式分解式。这里不如直接查表3.5简单，由N=4直接查表得到： 极点：             -0.3827j0.9239,  -0.9239j0.3827 归一化低通原型系统函数为                     式中，      1.0000，2.6131，3.4142，2.613 分母因式分解形式为               以上公式中的数据均取小数点后四位。 ③为将去归一化，先求3dB截止频率。按照（3.23）式，得到：               0.2139 将代入（3.24）式，得到：                 0.5072 此时算出的比题目中给的小，因此，过渡带小于指标要求。或者说，在0.5072时衰减大于30dB，所以说阻带指标有富余量。     将代入中，得到：                 （3）按照（3.1）和（3.10）式，将模拟滤波器系统函数转换成数字滤波器系统函数H(z)：                 ， 4.2 MATLAB程序及仿真 调用MATLAB信号处理工具箱函数进行设计。设计程序如下。 用脉冲响应不变法设计数字滤波器程序： T=1;      %T=1s fs=1/T; wp=0.2*pi/T;ws=0.6*pi/T;Rp=2;As=30;    %T=1s的模拟滤波器指标 [N,wc]=buttord(wp,ws,Rp,As,'s');    %计算相应的模拟滤波器阶数N和3dB截止频率wc [B,A]=butter(N,wc,'s');        %计算相应的模拟滤波器系统函数 [Bz,Az]=impinvar(B,A,fs);    %用脉冲响应不变法将模拟滤波器转换成数字滤波器 W=linspace(0,pi,100*pi);    %指定一段频率值 hs=freqs(B,A,W);          %计算模拟滤波器的幅频响应 hs0=abs(hs)/abs(hs(1)); hs1=20*log10(hs0); subplot(2,2,1); plot(W/2/pi,hs1);grid on        %绘出巴特沃斯模拟滤波器的幅频特性曲线 xlabel('f/Hz');ylabel('幅度/dB'); title('(a)巴特沃斯模拟滤波器'); axis([0 0.4 -50 0]) Hz=freqz(Bz,Az,W);        %返回频率响应 Hz0=abs(Hz)/abs(Hz(1)); Hz1=20*log10(Hz0); subplot(2,2,2); plot(W/2/pi,Hz1);grid on      %绘出巴特沃斯数字低通滤波器的幅频特性曲线 xlabel('f/Hz');ylabel('幅度/dB)') title('(b)数字滤波器') axis([0 0.4 -50 0]) 程序中，impinvar是脉冲响应不变法的转换函数，[Bz,Az]=impinvar(B，A，fs）实现用脉冲响应不变法将分子和分母多现式系数向量为B和A的模拟滤波器系统函数转换成数字滤波器的系统函数H(z),H(z)的分子和分母多项式系数向量为Bz和Az，fs为采样频率。 取T=1s时的运行结果： N=4 模拟滤波器系统函数分子和分母多项式系数向量B和A： B=[ 0        0        0        0    0.3994] A=[ 1.0000    2.0774    2.1578    1.3129    0.3994] 数字滤波器的系统函数H(z)分子和分母多项式系数向量Bz和Az： Bz=[ -0.0000    0.0385    0.0893    0.0137        0] Az=[  1.0000  -2.0105    1.7768  -0.7502    0.1253] 由Bz和Az写出数字滤波器系统函数： T=1s时巴特沃斯模拟滤波器和数字滤波器的损耗函数分别如下图(a)和(b)所示。 图4.1 巴特沃斯模拟滤波器和数字滤波器的损耗函数 第5章 总结 通过一周数字信号处理的课程设计，学到了很多知识，同时也遇到了一些问题，如：文档中公式如何编辑、怎样自动生成目录、MATLAB中一些库函数的运用等，我们通过查找Word中一些功能的使用，查阅了MATLAB基础教程中函数的使用方法使一些相关问题得以解决。在设计的过程中，我们熟悉了MATLAB的运行环境，了解了更多有关于Matlab软件的知识，初步掌握了MATLAB语言在数字信号处理中一些基本库函数的调用和编写基本程序等应用；通过对数字滤波器的设计，熟悉了滤波器设计的一般原理，对滤波器有了一个感性的认识；学会了数字低通滤波器设计的一般步骤；加深了对脉冲响应不变法的理解和认识。 同时，我们也认识到了自身知识的不足，我们上课学习的是基础理论知识，在实际的操作运用中，我们需要学习的还有很多，通过这次课程设计，使我们对所学知识进一步巩固和加深理解了，我们必须完善自身的知识体系，不断学习，所谓“学无止境”，我们要不断努力， 尽自己的最大努力获取成功。 参考文献 [1]  高西全，丁玉美.数字信号处理.3版.西安：西安电子科技大学出版社，2008 [2]  胡广书.数字信号处理—理论、算法、与实现.北京：清华大学出版社，1998 [3]  薛山.MATLAB基础教程.1版.北京：清华大学出版社，2011 [4]  陈怀琛，吴大正，高西全.MATLAB及在电子信息课程中的应用.北京：电子工业出版社，2006 [5]  王宏.MATLAB 6.5及其在信号处理中的应用.北京：清华大学出版社，2004