数字滤波器设计

1.1 第一章 绪论引言 随着信息时代的到来，数字信号处理的理论和技术更加成熟，它应用广泛，渗透到各个重要的科学领域。作为数字信号处理的一项重大突破就是数字滤波技术。数字滤波在通信，语音编码，雷达等许多领域中有着十分广泛的应用，目前，数字滤波器的设计图像处理，数据压缩等方面的应用取得了令人瞩目的进展和成就，近年来迅速发展起来的小波理论，由于其局部分析性能的友谊，在图像处理中的应用研究，尤其是在图像压缩，图像去噪等方面的应用研究，收到了越来越多的关注。由于数字滤波器的概念比较抽象，加上其数值计算又比较繁琐，所以借助好的计算机软件来进行辅助设计，是数字滤波器研究领域的一个发展趋势。这样的软件有很多种，其中最具代表性的就是MATLAB。MATLAB 是美国 MathWorks 公司自 20 世纪 80 年代中期推出的数学软件，优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出。MATLAB 是“ 矩阵实验室”（MATrix LABoratoy） 的缩写，它是由美国 Mathworks 公司于 1984年正式推出的，是一种以矩阵运算为基础的交互式程序语言，专门针对科学、工程计算及绘图的需求。MATLAB 是功能强大的科学及工程计算软件，它不但具有以矩阵计算为基础的强大数学计算和分析功能，而且还具有丰富的可视化图形表现功能和方便的程序计算能力。MATLAB 的应用领域极为广泛，除数学计算和分析外，还被广泛地应用于自动控制、系统仿真、数字信号处理、图形图像分析、数理统计、人工智能、虚拟现实技术、通信工程、金融系统等领域，因此，MATLAB 是面向 21 世纪的计算机程序设计及科学计算语言。 Simulink 是 MATLAB 最重要的组件之一，它提供一个动态系统建模、仿真和综合分析的集成环境。在该环境中，无需大量书写程序，而只需要通过简单直观的鼠标操作，就可构造出复杂的系统。Simulink 具有适应面广、结构和流程清晰及仿真精细、贴近实际、效率高、灵活等优点，并基于以上优点 Simulink 已被广泛应用于控制理论和数字信号处理的复杂仿真和设计。同时有大量的第三方软件和硬件可应用于或被要求应用于 Simulink。MATLAB/Simulink 实现数字滤波器的好处是效率高，直观而且容易修改。有高性能数值计算的高级算法，特别适合矩阵代数领域；有大量事先定义的数学函数，并且有很强的用户自定义函数的能力；有强大的绘图功能以及具有教育、科学和艺术学的图解和可视化的二维、三维图；基于 HTML 的完整的帮助功能；适合个人应用的强有力的面向矩阵(向量)的高级程序设计语言；与其它语言编写的程序结合和输入输出格式化数据的能力； 数字滤波器是指完成信号滤波处理功能的，用有限精度算法实现的离散线性非时变系统，其输入是一组（由模拟信号取样和量化的) 数字量，其输出是经过变换或说处理的另一组数字量。因此，它本身既可以是用数字硬件装配成德一台完成给定运算的专用数字计算机，也可将所需的运算编成程序，让通用计算机来执行。第二章 数字滤波器数字滤波器是指输入输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。它工作在数字信号域，它处理的对象是经由采样器件将模拟信号转换而得到的数字信号。数字滤波器一般由寄存器、延时器、加法器和乘法器等基本数字电路实现。数字滤波器一词出现在60年代中期。由于电子计算机技术和大规模集成电路的发展，数字滤波器已可用计算机软件实现，也可用大规模集成数字硬件实时实现。数字滤波器具有稳定性高、精度高、灵活性大等突出优点。随着数字技术的发展，用数字技术实现滤波器的功能愈来愈受到人们的注意和广泛使用。这里所说的数字滤波器是指理想带通，低通等的频率选择数字滤波器。数字滤波器设计的一个重要步骤是确定一个可实现的传输函数H(z)，这个确定传输函数 H(z)的过程称为数字滤波器设计。数字滤波器的一般设计过程为：（1）按照实际需要，确定滤波器的性能要求（通常在频域内给定数字滤波的性能要求） 。 （2）寻找一满足预定性能要求的离散时间线性系统。 （3）用有限精度的运算实现所设计的系统。 （4）通过模拟，验证所设计的系统是否符合给定性能要求。数字滤波器是对数字信号实现滤波的线性时不变系统。数字滤波实质上是一种运算过程，实现对信号的运算处理。输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器本质上是一个完成特定运算的数字计算过程，也可以理解为是一台计算机。描述离散系统输出与输入关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。时域离散系统的频域特性: ,数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化的过程中，使信号按预定的形式变化。数字滤波器有多种分类，从数字滤波器功能上分可分为低通、高通、带阻、带通滤波器，数字滤波器有低通、高通、带通、带阻和全通等类型。它可以是时不变的或时变的、因果的或非因果的、线性的或非线性的。根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应滤波器（IIR）和有限长冲激响应滤波器（FIR） 。数字滤波器指标：一般来说，滤波器的幅频特性是分段常数的，以低通为例，在通带内逼近于 1，阻带内逼近与 0，实际设计的滤波器并非是锐截止的通带和阻带两个范围，两者之间总有一个过渡带。在设计滤波器时事先给定幅频特性允许误差，在通带范围内幅度响应以误差 逼近于 1，在阻带内幅度响应以误差 逼近于 20。（1）weHrjwcj,2|)(|1式中 wc 和 wr 分别为通带边界频率和阻带边界频率，wr-wc 为过渡带。在具体的技术指标中往往用通带波动 来表示 ，用最小1阻带衰减 At 来表示 ，其具体的对应公式这里就不详述了。2IIR 数字滤波器设计IIR DF 的冲激响应 h(n)是无限长的，其输入输出的关系为：(2)()(inxhnyi系统函数为= (3)nnzhzH)()(kkmrrzab10设计无限长单位脉冲响应（IIR）数字滤波器一般可有三种方法。第一种方法，先设计一个合适的模拟滤波器，然后将其数字话，即将 S 平面映射到 Z 平面得到所需的数字滤波器。模拟滤波器的设计技巧非常成熟，不仅得到的是闭合形式的公式，而且设计系数已经表格化了。因此，由模拟滤波器设计数字滤波器的方法准确，简便，得到普遍采用。对于这种方法，工程上有两种常见得变换法 脉冲响应不变法及双线性变换法。第二种方法，在 Z 平面直接设计 IIR 数字滤波器，给出闭合形式的公式，或者以所希望的滤波器响应作为依据，直接在 Z 平面上通过多次选定极点和零点的位置，以逼近该响应。第三种方法，利用最优化技术设计参数，选定极点和零点在 Z 平面上的合适位置，在某种最优化准则意义上逼近所希望的响应。但一般不能得到滤波器的系数（即零，极点的位置）作为给定响应的闭合形式函数表达式。优化设计需要完成大量的迭代运算，这种设计法实际上也是 IIR 滤波器的直接设计。本文着重介绍由模拟滤波器设计相应的 IIR 数字滤波器的方法。（1）脉冲响应不变法脉冲响应不变法是使数字滤波器的单位脉冲响应序列 h(n)逼近模拟滤波器的冲激响应 ，让 h(n)正好等于 的采样值。()aht ()aht设已有满足要求的模拟滤波器, 则可 ()aHs()aht()()anhT()Hz因为: ,thn的 图 形 的 图 形公式导出:具体转换如下:设(以一阶极点为例)(4)1NkaAHss作拉氏反变换,得 11()()()kNstaakhtFTsAeu采样得 1()()kNsnTakhnTAeu作 Z 变换 ,得(5)1011()()k kNNsnTnsTnk kAHzAeuzez 与 极点关系为: ()z()as ksTkze(6)一般对应关系 sTjTjjzree(7),T图 1 S 平面到 Z 平面变换示意图所以, 模拟系统稳定因果数字系统稳定因果。按照脉冲响应不变法，从 S 平面到 Z 平面的映射不是单值关系，而是先将 在 S 平面沿虚轴作周期严拓，再按照映射关系将()aHs映射到 Z 平面，得到 ，因此，脉冲响应不变法只适用于带()as ()Hz限的滤波器（如低通、带通）。在 Matlab 中利用 M 文件 impinvar 可以对模拟传输函数实行脉冲响应不变法。（2）双线性变换法脉冲响应不变法不适带阻和高通滤波器的设计，因为高频带为通带，前述方法易引起混频。故希望: 平面虚轴z 平面单位圆sjO1rOs平 面 0z平 面 Treks kz一周, 且应有,:0:0,因为 , 所以选变换tan(/2)(8)tK其中 可取任意正常数, 后面将导出 .2/KT设计思路: .设计出模拟滤波tan2, ,psspspsK器 转化成数字滤波器.图 2 数字域频率与模拟域频率的对应关系转化公式推导如下: sin(/2)tan2cojjjK/2/2/1jj jjj jeeKK因只关心频率转换, 故可设 , , 则有sjjze, (称为双线性变换) (9)1zsK所以模拟滤波器转换成数字滤波器的公式为-10 -5 0 5 10-202w2 atan(w) ppO(10)1()()zasKHz由双线性变换公式, 可得 , 视为两复平面变换, 再由Ksz.,jsjzre可得(11)22(1)()Kr从上式可得:时, , 上虚轴 上单位圆周。01rsz时, , 上左半平面 上单位圆内。时, , 上右半平面 上单位圆外。01rsz故若模拟滤波器稳定，则双线性变换后数字滤波器也稳定。由于双线性变换法是一种单值映射，因此消除了频率混叠的现象。双线性变换法的缺点是模拟频率与数字频率间的非线性，这种非线性关系要求被变换的连续系统的幅度响应是分段常数型的（某一频率范围内幅度响应近似于某一常数），不然所映射出的数字频率响应相对于原来的模拟频率响应会产生变形。为解决双线性变换中的频率非线性关系，我们采用预畸的方法，即，其中 K=2/T。tan2K在 Matlab 中利用 M 文件 bilinear 可以对模拟传输函数实行双线性变换法。MATLAB 中 IIR 数字滤波器的设计过程包括两步。第一步，根据给定指标，确定滤波器的阶数 N 和频率缩放因子 Wn。第二步，利用这些参数和给定的波纹参数，确定传输函数的关系。阶数估计：利用双线性变换法设计数字滤波器时，首先要对 IIR 数字滤波器的阶数进行估计，相应的 M 文件为：buttord 用于巴特沃斯滤波器，cheb1ord 用于切比雪夫 1 型滤波器，cheb2ord 用于切比雪夫 2 型滤波器，ellipord 用于椭圆滤波器。滤波器的设计：对于基于双线性变换法的 IIR 滤波器设计，对应于四种逼近技术（即巴特沃斯、切比雪夫 1 型和 2 型及椭圆逼近） ，MATLAB 工具箱中有相应的函数。特别地可以用到下面的 M 文件：butter 用于巴特沃斯滤波器的设计，cheby1 用于切比雪夫 1 型滤波器的设计，cheby2 用于切比雪夫 2 型滤波器的设计，ellip 用于椭圆滤波器的设计。这些函数的输出可以是滤波器传输函数分子和分母的系数向量，也可以是滤波器的零极点向量和标量增益因子。同时，利用 zp2tf 可以由滤波器的零极点向量和标量增益因子得到传输函数分子和分母的系数向量。相应地，利用函数 zp2sos 可以得到传输函数分子和分母系数向量的二次项因子。在计算出传输函数的系数之后，可以利用 M 文件 freqz 来计算频率响应。2.2 FIR数字滤波器设计FIR DF的冲激响应h(n)是有限长的，M阶FIR DF可以表示为：(12)10)()(iinxhny其系统函数为:(13)10)()(MnnzhzH与IIR 数字滤波器的设计不同，FIR滤波器的设计与模拟滤波器的设计没有任何联系。因此，FIR滤波器的设计基于对指定幅度响应的直接逼近，并通常要求其具有线性相位响应。为了保证滤波器具有线性相位特性，滤波器系数必须满足条件：h(n)= h(M-1-n)。目前关于FIR 滤波器的设计方法主要有三种，即窗函数法，频率取样法和切比雪夫等波纹逼近的最优化设计方法。一般应用较多的是第一种和第三种方法。这是因为窗函数法比较简单，可应用现成的窗函数公式，在技术指标要求不严格的情况下市比较灵活的。最优化设计法必须借助计算机计算，但是它能得到最佳的等波纹的线性相位FIR滤波器。目前切比雪夫等波纹的线性相位FIR 滤波器的计算机机助设计程序已经比较完善，由于采用了REMEZ迭代算法，所以设计效率也很高，在应用中越来越占优势。（1）窗函数法 一般设计过程总是先给定一理想的滤波器频率响应 ，然)(jwdeH后设计一个FIR 滤波器，用它的频率响应 来逼近理jnMnjwhe10)(想的 。这种逼近中最直接的方法，是在时域中用FIR滤波器的)(jwdeH单位脉冲响应h(n)去逼近理想的单位脉冲响应 。因而，先由d的IDTFT 导出j )(nhd(14)eHjwnd)(21由于 是矩形频率特性，故 一定是无限长的序列，且是非)(jwdeh因果的。然而FIR 滤波器是有限长的，所以用有限长的h(n) 来逼近无限长的 ，最简单的方法是截取 中最重要的一段，将无限长nh)(nd的 截取成长度为M的有限长序列，等效于再 上施加了一个)(d )(nhd长度为M 的矩形窗口，更为一般的，可以用一个长度为M的窗口函数w(n) 来截取 ，即)(d(15)()(nwhd这一方法通常称为窗函数法，窗口函数的形状及长度M的选择是窗函数法的关键。下面我们一低通为例，了解一下窗函数法的运用:提出希望频率响应函数（低通）图 3 理想低通滤波器的频响线性相位, 具有片断特点, 即O0.25|()|jdHe.1|()0|jj cdceH算出 -1()()d2jjndhnecjjn(无限长)sin(c图 4 理想低通的单位脉冲响应（无限长的一部分）加窗 ,长 , 得 ()wnN(*)()()dhnwn要线性相位, 就要 关于 偶对称,而 关于 偶对称, 故1/2Ndhn要求 ()/所以要求 关于 偶对称.()wn(1)/2N0 10 20 3000.510 10 20 30-0.100.10.20.3图 5 窗函数 图 6 加窗后的单位脉冲响应再回过来检验 是否满足精度要求.()jHe0 1 2 300.51O.2|()|jdHe.()je0 1 2 300.51 0 10 20 3000.510 10 20 30-0.100.10.20.3图7 图 4的脉冲响应的频响 图8 理想频响与实际频响的对比若基本满足, 则依截取的 , 制硬件, 编软件.()hn为便于选择使用, 将 5 种常见的窗函数基本参数如表 1 所示。表 1 5 种常见的窗函数基本参数类型窗函数的旁瓣峰 n过渡带宽度 B加窗后滤波器的阻带最小衰减 srectwin 13 4/N 21bartlet 三角 25 8/N 25hanning 31 8/N 44hamming 41 8/N 53blackman 57 12/N 74（2）频率取样法窗口设计法事从时域出发，把理想的 用一定形状的窗口函)(nhd数截取成有限长的h(n)，以此h(n) 来近似理想的 ，从而频率响应也近似于理想的频率响应 。我们知道一个有限长序列)(jweH)(jwdeH可以通过其频谱的相同长度的等间隔采样值准确地恢复原有的序列。频率采样法便是从频域出发，对理想的频率响应 加以等间隔)(jwdeH采样)()|(2kedkMwjd(16)然后以此 作为实际FIR滤波器的频率特性的离散样本H(k)，即)(kHd(17)1,.20,|)()()(2MkeHkMwjdd由H(k) 通过IDFT可求出有限长序列h(n)为(18)1,.0,)(1)(02nekMnhkj利用M个频率的离散样本H(k)同样可求出FIR滤波器的系统函数H(z)及频率响应 。)(jweH(19)MjMkknn WzzhzH210110)()( ， 其 中令 可得到滤波器的频率响应 。如果设计的是线性相位的jwe )(jweFIR数字滤波器，其采样值H(k)的相位的幅度一定要满足特定的约束条件，这个设计时一定要注意。（3）最优化设计法最优化设计法事以最佳一致逼近（最大误差最小化）理论为基础，利用雷米兹算法设计的具有等波纹特性的设计方法。具体设计步骤如下：对设计指标进行归一化处理。确定remezord 函数所需要的参数。包括归一化边界频率、各频带的幅度要求和波纹要求等。归一化边界频率总是从0开始到1结束，故只需递增列出中间的边界频率；频带幅度要求不含过渡区，个数是边界频率个数的一半加1；波纹要求是频带内幅度允许的波动要求，与分贝间的关系是： (20)(log20)1(log20),1(log20 21spRR利用remezord 函数确定 remez所需参数。调用remez函数进行设计。利用freqz函数验算技术指标是否满足要求。3 数字滤波器的MATLAB设计3.1 直接程序设计法（1）IIR 的直接程序设计法例如欲设计一数字（IIR）带阻滤波器，其数字域指标为：数字阻带边缘频率分别为0.4 和0.7 ，数字通带边缘频率为0.25 和 0.8 ，通带波动为1db最小阻带衰减为40db。此题的MATLAB程序为：ws=0.4*pi 0.7*pi; %数字阻带边缘频率wp=0.25*pi 0.8*pi; %数字通带边缘频率rp=1 ; %通带波动（db）as=40; %阻带衰减n,wn=cheb2ord(wp/pi,ws/pi,rp,as);%根据给定指标，确定滤波器的阶数N和频率缩放因子Wnb,a=cheby2(n,as,ws/pi,stop);%返回的b,a分别为 H（z)的分子、分母。h,w=freqz(b,a,512);%返回的h,w分别为滤波器的频率响应及其频率plot(w/pi,abs(h);%画出频率响应（以w/pi为横轴）grid;xlabel(w/pi);ylabel(幅值);title(频率响应);程序运行结果为：图9 所设计的带阻滤波器的频率响应在设计中如果该滤波器的特性不满足要求，原有的参数必须做相应的调整，在程序中只需对参数做新的设定就可以得到所需的滤波器。接下来我们来看看此题所设计的滤波器的滤波效果：S为含有3个频率成分的信号（归一化频率（w/2 ）分别为0.1、0.3、0.45） ，用所设计的滤波器滤除归一化频率为0.3的成分。n=0:100;s1=sin(pi*0.2*n);s2=sin(pi*0.6*n);s3=sin(pi*0.9*n);s4=s1+s3;s=s1+s2+s3;sf=filter(b,a,s);subplot(311)stem(n,s);title(滤波前的信号);subplot(312);stem(n,sf);title(滤波后的信号 );subplot(313);stem(n,s4);title(想要保留的信号);程序运行的结果为：图10 采用filter 函数进行数字滤波前后信号比较示意图由图可以看出，滤波后的信号与想要保留的信号基本一致（相位有些许偏差，但基本一致） ，所以我们可以说该滤波器基本满足了以上所提出的滤波要求。（2）FIR的 直接程序设计法例如欲设计一个线性相位数字（FIR）带通滤波器，其数字域指标为：数字通带边界频率为0.35 和0.65 ，数字阻带边界频率为0.2和0.8 ，通带波动为1db，最小阻带衰减为 60db。 FIR数字滤波器的窗函数法此题的MATLAB程序为：ws1=0.2*pi;wp1=0.35*pi;wp2=0.65*pi;ws2=0.8*pi;as=60;tr=min(wp1-ws1),(ws2-wp2);M=ceil(11*pi/tr)+1; %滤波器的阶数，程序运行后M=75n=0:1:M-1;r=(M-1)/2;%r为群时延wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;hd=sin(wc2*(n-r)+eps)./(pi*(n-r)+eps)-sin(wc1*(n-r)+eps)./(pi*(n-r)+eps);%hd为理想滤波器的脉冲响应w_bla=(blackman(M);%长度为M的blackman窗h=hd.*w_bla;%h为滤波器的实际脉冲响应stem(n,h);title(滤波器的实际单位脉冲响应);freqz(h,1,512);title(幅度响应和相位响应);图11 所设计的滤波器的实际单位脉冲响应由上图可知滤波器的实际脉冲响应h是偶对称的，即h(n)=h(M-1-n),故该滤波器满足FIR线性相位的条件，该滤波器是线性相位的FIR滤波器。图12 所设计的带通滤波器的幅度和相位响应由滤波器的相位特性也可以看出该滤波器是线性相位的FIR滤波器。接下来我们来看看此题所设计的滤波器的滤波效果：S为含有3个频率成分的信号（归一化频率（w/2 ）分别为0.05、0.2、0.45） ，用所设计的滤波器滤除归一化频率为0.05和0.45的成分。l=0:100;s1=sin(0.1*pi*l);s2=sin(0.4*pi*l);s3=sin(pi*0.9*l);s=s1+s2+s3;sf=filter(h,1,s);subplot(311)stem(l,s);title(滤波前的信号);subplot(312);stem(l,sf);title(滤波后的信号);subplot(313);stem(l,s2);title(想要保留的信号);图13 采用filter 函数进行数字滤波前后信号比较示意图由上图可知滤波后的信号和想要保留的信号的幅度和频率基本一致（滤波后的信号相对于想要保的信号有一个相位延迟，这是线性相位FIR滤波器的群延迟引起的，此滤波器留的群延迟r=(M-1)/2=37） ，所以我们可以说该滤波器基本满足了以上所提出的滤波要求。 FIR数字滤波器的频率采样法此题的MATLAB程序为：M=40;%取滤波器的阶数为40al=(M-1)/2;%群时延n=0:M-1;T2=0.59417456;T1=0.109021;Hrs=zeros(1,5),T1,T2,ones(1,7),T2,T1,zeros(1,9),T1,T2,ones(1,7),T2,T1,zeros(1,4);%采样值的幅值k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1;angH=-al*(2*pi)/M*k1,al*(2*pi)/M*(M-k2);%采样值的相位H=Hrs.*exp(j*angH);h=real(ifft(H,M);%长度为M的单位脉冲响应stem(n,h);title(滤波器的实际单位脉冲响应);freqz(h,1,512);title(幅度响应和相位响应);图14 所设计的滤波器的实际单位脉冲响应由图14可知滤波器的实际脉冲响应h是偶对称的，即h(n)=h(M-1-n),故该滤波器满足FIR线性相位的条件，该滤波器是线性相位的FIR滤波器。图15 所设计的带通滤波器的幅度和相位响应由滤波器的相位特性也可以看出该滤波器是线性相位的FIR滤波器。此滤波器的群延时为al=(M-1)/2=19.5。接下来我们来看看此题所设计的滤波器的滤波效果：S为含有3个频率成分的信号（归一化频率（w/2 ）分别为 0.05、0.2、0.45） ，用所设计的滤波器滤除归一化频率为0.05和0.45的成分。l=0:100;s1=sin(0.1*pi*l);s2=sin(0.4*pi*l);s3=sin(pi*0.9*l);s=s1+s2+s3;sf=filter(h,1,s);subplot(311);stem(l,s);title(滤波前的信号);subplot(312);stem(l,sf);title(滤波后的信号);subplot(313);stem(l,s2);title(想要保留的信号);图16 采用filter 函数进行数字滤波前后信号比较示意图同上面分析相似，滤波后的信号和想要保留的信号的幅度和频率基本一致（滤波后的信号相对于想要保的信号有一个相位延迟，这是线性相位FIR 滤波器的群延迟引起的，此滤波器留的群延迟（r=(M-1)/2=19.5） ，所以我们可以说该滤波器基本满足了以上所提出的滤波要求。 FIR数字滤波器的最优设计法此题的MATLAB程序为：%设计指标ws1=0.2*pi;wp1=0.35*pi;wp2=0.65*pi;ws2=0.8*pi;rp=1;as=60;%设置边界频率和幅度要求F=ws1/pi,wp1/pi,wp1/pi,ws2/pi;A=0,1,0;%设置各频带的波纹要求devp=(10(rp/20)-1)/(10(rp/20)+1);devs=10(-as/20);dev=devs,devp,devs;%确定remez 参数，其中滤波器的阶数为（ N+1） ，程序运行后得到N=26N,Fo,Ao,W=remezord(F,A,dev);%调用remez 函数进行设计h=remez(N,Fo,Ao,W);n=0:N;stem(n,h);title(滤波器的单位冲激响应);freqz(h,1,512);title(幅度响应和相位响应);图17 所设计的滤波器的实际单位脉冲响应由图17可知滤波器的实际脉冲响应h是偶对称的，即h(n)=h(N-n),故该滤波器满足FIR 线性相位的条件，该滤波器是线性相位的FIR 滤波器。图18 所设计的带通滤波器的幅度和相位响应由滤波器的相位特性也可以看出该滤波器是线性相位的FIR滤波器。此滤波器的群延时为al=(N)/2=13。接下来我们来看看此题所设计的滤波器的滤波效果：S为含有3个频率成分的信号（归一化频率（w/2 ）分别为0.05、 0.2、0.45） ，用所设计的滤波器滤除归一化频率为0.05和0.45的成分。l=0:100;s1=sin(0.1*pi*l);s2=sin(0.4*pi*l);s3=sin(pi*0.9*l);s=s1+s2+s3;sf=filter(h,1,s);subplot(311)stem(l,s);title(滤波前的信号);subplot(312);stem(l,sf);title(滤波后的信号);subplot(313);stem(l,s2);title(想要保留的信号);图19 采用filter 函数进行数字滤波前后信号比较示意图同上面分析相似，滤波后的信号和想要保留的信号的幅度和频率基本一致（滤波后的信号相对于想要保的信号有一个相位延迟，这是线性相位FIR 滤波器的群延迟引起的，此滤波器留的群延迟r=(N)/2=13） ，所以我们可以说该滤波器基本满足了以上所提出的滤波要求。3.2 利用信号处理工具箱SPTOOL设计法SPTOOL是信号处理工具箱中一个具有交互式图形用户界面的信号处理工具，专门用于完成常用的数字信号处理任务。通过这个工具，只需要鼠标简单的操纵鼠标（点击或拖动） ，就可以完成载入、观察、分析、实现和设计数字滤波器并进行谱分析等甚至十分复杂的数字信号处理任务，而不需要用户对数字滤波器的设计原理非常熟悉。下面是用SPTOOL 工具设计的IIR 滤波器对上述信号进行滤波。SPTOOL使用步骤：导入信号、滤波器、频谱，设计滤波器，对信号滤波，分析输入输出信号的谱。1 IIR滤波器的设计创建并导入信号源，在MATLAB命令窗口输入命令：fs=4000;（这里取采样频率为4000Hz）;n=0:100;s1=sin(pi*0.2*n);s2=sin(pi*0.6*n);s3=sin(pi*0.9*n);s4=s1+s3;s=s1+s2+s3;此时，变量fs、n、s、s1、s2、s3、s4 将显示在MATLAB 的workspace列表中。在MTLAB命令窗口中输入sptool，即可弹出SPTOOL的主界面 ,如图所示： 图20 SPTOOL 的主界面这里我们首先来设计前面的IIR带阻滤波器，前面的数字指标转化为模拟指标为：通带频率分别为500Hz和1600Hz，阻带频率为800Hz和1400Hz，通带波动为 1db，最小阻带衰减为40db。以上三个信号转换为模拟域的频率分别为400Hz，1200Hz，1800Hz 。1）导入信号。使用【file/import】可以导入信号、滤波器和谱。信号的来源可以是MATLAB工作空间变量或 MAT数据文件，也可以是在DATA 前面的文本框中直接输入信号的数据。这里信号源仍然利用前面程序设计法中的混频信号，将信号s和采样频率fs=4kHz导入并命名为sig1，将信号s2导入命名为sig2,将信号s4导入命名为sig4。2）滤波器的设计。单击Filter栏中的New 按钮打开滤波器设计工具。在界面的最上面制定滤波