数字信号处理教程-MATLAB释义与实现第07章.ppt

1,第七章 FIR滤波器设计,2,FIR滤波器设计,本章的主要内容为：7.1 数字滤波器设计的基本要求 7.2 线性相位滤波器的特性 7.3 设计FIR滤波器的窗函数法 7.4 设计FIR滤波器频率样本法7.5 FIR滤波器的最优设计法 7.6 FIR滤波器设计的一些深入问题,3,滤波器设计概述,从本章起将着重研究信号分析的逆问题，那就是给定了输入输出指标，设计系统。这种可以人为调整参数来满足给定要求的硬件或软件系统，称为滤波器。这种问题没有唯一的解，因为可能有多种类型的滤波器都能满足给定的要求。要人为地进行选择，就要考虑更多的因素,所以设计比之分析更加接近与工程实际。有两种类型的数字滤波系统。本章将讨论FIR数字滤波器的几种基本设计算法；下一章讨论IIR数字滤波器的设计问题。,4,滤波器设计概述,本书中的设计大多数是针对选频型的滤波器。从幅频特性的要求而言，有低通、高通、带通和带阻滤波器；从相频特性来说，往往要求线性相位。在FIR滤波器设计中，还会涉及到像微分器或希尔伯特变换器之类的系统，尽管它们不是选频滤波器，但是也遵循选频滤波器的设计技术。在计算技术高度发展的今天，几乎所有的电子设计软件中都把滤波器设计作为其模块之一。那种用查表、手工计算的方法在工程实践中已经濒于淘汰。,5,滤波器设计概述,数字滤波器的设计要经过三个步骤：确定指标：在设计一个滤波器前，必须有一些指标。这些指标要根据应用确定。模型逼近：一旦确定了技术指标，就可利用前面学习过的基本原理和关系式，提出一个滤波器模型来逼近给定的指标体系。这是滤波器设计所要研究的主要问题。实现：上面的两步的结果得到的滤波器，通常是以差分方程、系统函数或脉冲响应来描述的。根据这个描述，可以如第六章所讨论的那样，用硬件或计算机软件来实现它。,6,滤波器设计概述,数字滤波器的主要指标：在很多实际应用，象语音或音频信号处理中，数字滤波器用来实现选频功能。因此，指标的形式通常为频域中的幅度和相位响应。幅度指标可以两种方式给出。第一种，叫做绝对指标，它提出了对幅度响应函数|H(j)|的要求。这些指标一般可直接用于FIR滤波器。第二种方法叫做相对指标，它以分贝(dB)值的形式提出对幅频特性的要求。其值定义为：,7,FIR滤波器设计概述,相对幅度指标：由于定义中所包含的归一化，滤波器的相对幅频特性最高处的值为0dB，又由于定义中的负号，幅频特性小的地方，其dB值反而是正的。为了说明这些指标，以低通滤波器为例进行讨论。首先看下面两种指标图，它们分别用绝对指标和相对指标表示。,8,FIR滤波器设计概述,9,FIR滤波器设计概述,绝对指标：图a是典型的低通幅特性绝对指标：0,p段叫做通带，1是在理想通带中能接受的振幅波动（或容限）,011。s,段叫做阻带，2是阻带中能接受的振幅波动（或容限），021。p,s叫做过渡带，在此段上对幅度响应通常没有限制，也可以给些弱限制。相对指标（dB）：低通的典型相对技术指标如图b所示，三种频带的定义不变，只是1，2 换成了Rp，As。其中 Rp是通带波动的dB数;As是阻带衰减的dB数；,10,FIR滤波器设计概述,绝对指标和相对指标间的变换关系：可以用MATLAB语句fgp711准确地求出绝对指标和相对指标的对应数值。并画出如图中所示曲线。,11,FIR滤波器设计概述,12,FIR滤波器设计概述,相位指标主要是线性相位条件。理想的滤波器相位应该和频率成正比，也就是()=-。因此满足此式的系统对信号中所有的频率分量都具有相同的时间延迟(无量纲)，因此无波形失真。略低一点的要求是相位和频率成线性关系，即()=0-。IIR滤波器无法做到线性相位，通常用群迟延g()来评价。其定义为相位对频率的导数，它愈接近常数愈好。,13,FIR滤波器设计概述,上述指标是针对低通滤波器的，其他类型的选频滤波器，像高通和带通滤波器，它们的技术指标可作类似的规定，并且可以由低通滤波器演化而得。本章讨论FIR数字滤波器的设计和逼近。FIR滤波器在设计和实现上具有如下的优越性：(a).相位响应可为严格的线性；(b).不存在稳定性问题；(c).只包含实数算法，不需要递推运算。,14,7.2 线性相位滤波器的特性,线性相位条件要求滤波器分子系数bn满足对称性条件。设滤波器的系数长度为N，则这些系数应关于中心点=n0=(N-1)/2对称。偶对称时，h(n)=h(N-n-1)；而奇对称时，h(n)=-h(N-n-1)。再考虑到N可以为奇数或偶数，总共有四种类型的线性相位FIR滤波器。在讨论线性相位滤波器频率响应时需要引进幅特性正负号的概念。以往常设幅特性(Magnitude Response)为正数，因为幅特性的反号可以用相特性加减来补偿。当相位特性要求线性，不得随便增减时，幅特性就必须分出正负,称为符幅特性(Amplitude Response)。,15,线性相位滤波器的特性,为了计算符幅特性，需要直接作解析计算，求出频率特性的解析式，并将它分解为幅度和相位特性，即，如果()满足线性相位条件，这时的A()就是符幅特性。对于长度为N的序列h(n)，先考虑两个对称的系数h(n)和h(N-1-n)，它们的傅立叶变换为（注意运用样本值相等的条件使两项合一）：,16,I型线性相位滤波器的特性,下面分别四种类型进行推导：1类型I，系数对称，N奇数序列h的下标n从0起算，到N-1为止，中点位置是L=(N-1)/2。h(L)是一个孤项，其他可配成L对系数。将它们的符幅特性加起来，得到：,17,I型线性相位滤波器的特性,为了使四种类型滤波器的公式简明统一，引进两个参数。一个是令，它就是(7.1.5)式中的值，用以反映延迟的样本数或群延迟，它可以是分数。另一个是。即把值向下取整，因为要把L用作下标，它必须是整数。当N为奇数时，=L。在推导中，还注意前一个求和号是从0到L-1,后一个求和号是从0到L,那是把孤项也当成成对项求和，然后再减去一个孤项h(L)。这样做的目的也是使四种类型滤波器的公式统一。,18,I型线性相位滤波器的特性,显然，为相角项，相角-,它与成严格的线性关系；A()就是符幅项。可见，其符幅特性由L+1个余弦项叠加组成。n=0时的符幅项为2h(0)cos(L)，在=02频段间波动L个周期；余弦函数在=0，和2处都不等于零，因此类型I线性相位滤波器既可以用作低通滤波器(在=0处，幅特性不为零)；也可以用作高通滤波器(在=处，幅特性不为零)；且可用作带通和带阻滤波器。,19,II型线性相位滤波器的特性,类型II，系数对称，N偶数如果N为偶数，那末全部系数都可以配对，不会出现中心点的单项，一共有N/2组对称出现的系数，频率特性就成为 此时=(N-1)/2将不是整数,也就是说，对称中心将在两个样本点的中间。其相角特性仍为()=-,它仍与成严格的线性关系。,20,II型线性相位滤波器的特性,用整数 来表示求和号的上限。得出其符幅特性的表达式可见，它的符幅特性也由L+1个余弦分量构成。这样，不管N为奇数或偶数时，L都是整数，求和的项数都是L+1。由于=(N-1)/2是分数，所有余弦分量cos(-n)中都含一个0.5。在=处，就会出现符幅特性必定为零的情况。所以类型II不能用作高通和带阻滤波器。,21,III型线性相位滤波器的特性,3类型III，系数反对称，N奇数当系数序列反对称时，在(7.2.1)式的推导过程中，h(N-1-n)项将反号而有h(N-1-n)=-h(n)，两个反对称项合成的结果是：所以只要把类型I 结果中的cos换成sin，得到：,22,III型线性相位滤波器的特性,当N为奇数时，=L=整数。相特性为()=/2-(N-1)/2=/2-,它虽然仍与成线性关系，但多了一个常数项/2。注意反对称序列的对称中心n=L处，序列的值h(L)=0。所以可以放入求和号中，得到它的符幅特性为：正弦函数在=0，和2处都等于零，因此类型III线性相位滤波器的符幅特性在=0，和2处都等于零。它不能用作低通，也不能用作高通滤波器。故不宜用于选频滤波器。,23,IV型线性相位滤波器的特性,4类型IV，系数反对称，N偶数与类型II滤波器相仿，若N为偶数，则共有N/2组数值相同反对称出现的系数，不出现单项，频率特性式就成为此处仍有=(N-1)/2，因此将不是整数,带有小数部分0.5。故对称中心将在两个样本点的中间。它的相角特性为()=/2-,虽然仍与成线性关系，但多了常数项/2，故不通过原点。,24,IV型线性相位滤波器的特性,类型IV的滤波器符幅特性为：因为L是整数，它的符幅特性也由L+1=N/2个正弦分量构成。因为sin(0)=0，在=0处，符幅特性必定为A(0)=0。这样的特性显然不能用作低通和带阻滤波器。实际上这种类型的滤波器也不适用于选频滤波。把这四种类型的线性相位滤波器的相角特性和符幅特性用对比的方法画在一张图上如表。,25,26,27,符幅特性的计算举例,例7.2.1 设系统的系数向量为b=-1,2,4,2,-1，a=1，求其幅频特性及符幅特性，并作比较。解：先用手工计算，因为a=1,说明这是一个FIR滤波器，观察系数向量b，其长度N=5，阶次为M=N-1=4，传递函数为H(z)=-1+2z-1+4z-2+2 z-3 z 4注意h(n)=bn，其频率响应为：调用abs(H)函数只能求出常为正的幅特性。,28,幅特性和符幅特性的算法对比,要求出符幅特性，必须把对称系数配对计算。将n=0,4组成一对，n=1,3组成一对，n=2处是对称中心，h(2)独立计算，于是有：由于对称配置，的高次和低次项搭配成为余弦，合成后的相位项都是，故其符幅特性是这三个幅度项的代数和：把本题的幅特性和符幅特性画成曲线如下图,29,幅特性和符幅特性的结果对比,左上是幅特性，右上是符幅特性；可见符幅特性取绝对值，就得幅特性。下方是相特性，右图是线性相位，左图的第二个突跳与幅特性反号对应。,30,计算符幅特性的MATLAB程序,本书编了一个通用的程序，它能自动判别滤波器的类型并能在(7.2.2a)(7.2.5a)四个公式中，自动选择正确的符幅特性公式进行计算。把这个函数程序命名为amplres.m，其调用方法：A,w,type,tao=amplres(h)其中：h=FIR滤波器的脉冲响应或分子系数向量 A滤波器的符幅特性 w取的频率向量，在0到pi之间取501个点 type线性相位滤波器的类型 tao符幅特性的群迟延,31,计算符幅特性的通用子程序,子程序amplres.m的核心语句为：N=length(h);tao=(N-1)/2;L=floor(n-1)/2);n=1:L+1;w=0:500*pi/500;%滤波器频率向量if all(abs(h(n)-h(N-n+1)1e-10)%若系数对称 A=2*h(n)*cos(N+1)/2-n)*w)-mod(N,2)*h(L+1);%对称条件下计算A（两种类型）elseif all(abs(h(n)+h(N-n+1)1e-10)end,32,线性相位滤波器的特性,在阅读这个程序时，注意以下几点：（1）。首先要注意程序中的n比公式中的n大1，这是因为MATLAB下标从1开始，而公式中的下标从零开始；公式中的求和指标从n=0到L在程序中是n=1到L+1，依此类推。（2）。关于判断序列对称和反对称的条件语句，似乎该用all(h(n)-h(N-n+1)=0，程序中则采用了all(abs(h(n)-h(N-n+1)1e-10)。这是因为生成h(n)和h(N-n+1)中会有计算误差，为了防止造成误判，要允许两者有一点差别；,33,线性相位滤波器的特性,（3）。计算符幅特性时，用行向量2h(n)与列向量cos(-n)相乘，实现求和号下的相乘累加运算。（4）。程序中采用了mod(N,2)函数，来区分奇偶的不同效果，避免用多个条件转移语句。在N为奇数时它等于0，而在N为偶数时它等于1。用mod(N,2)去乘一个变量，相当于在N为偶数时把这项变量取消。所以使适用于四种情况的计算程序程序变得如此简练。,34,符幅特性算例,给出由序列h0=3,-1,2,-3 构成的四种类型的系数向量：(a),h1=3,-1,2,-3,5,-3,2,-1,3;(b),h2=3,-1,2,-3,-3,2,-1,3;(c),h3=3,-1,2,-3,0,3,-2,1,-3;(d),h4=3,-1,2,-3,03,-2,1,-3;分别求出它们的符幅特性曲线，进行比较。解：对(a)的语句如下，其余类推。h1=3,-1,2,-3,5,-3,2,-1,3;A1,w1,typea,tao1=amplres(h1);typea,35,符幅特性算例的结果,运行程序得到typea=1,typeb=2,typec=3 和typed=4，相应的符幅特性曲线见图。可验证表7.2.1:类型I幅特性关于=对称，在=0和=处可以取任何值；类型II幅特性关于=反对称，在=0可以取任何值，在=处必定为零；类型III幅特性关于=反对称，在=0和=处都必定为零；类型IV幅特性关于=对称，在=0处必定为零，而=处可以取任何值；,36,线性相位滤波器的特性,37,线性相位滤波器的零极点分布,零极点应该从FIR滤波器正幂传递函数来判断。所以它具有与零点数目相同的M个极点，它们集中在z平面的原点处，成为M重极点。可以证明，对于具有对称性的系数向量，如果z1是(7.2.13)的分子多项式的根，那末它的倒数q1=1/z1也是这个多项式的根。,38,线性相位滤波器的零极点分布,如果z1是复数零点，根据实系数多项式的性质，它的共軛数z1*也是复数零点；现在又证明了它的倒数1/z1也是零点，因而倒数的共軛数z1*也是复数零点。四个复数零点必成组出现，如图。,39,7.3 设计FIR滤波器的窗函数法,第五章中研究了理想低通滤波器的频率特性应为：振幅(符幅)特性在通带内为1，阻带内为0；在通带内的相位特性与成线性关系。即 对应的理想单位脉冲响应(滤波器系数向量)为理想滤波器之所以不能实现，因为(1)其长度是无限的;(2)它是非因果序列。,40,设计FIR滤波器的窗函数法,为了构造特性接近理想低通滤波器的实际滤波器，应当使它具备如下特性：(1).因果序列，在n0处无样本，起点在n=0处；(2).有限长度，长度为N(阶次为N-1),系数下标从0到N-1;(3).系数对称，以保证线性相位。对称中心应在n0=(N-1)/2处；也即把理想脉冲响应右移n0;用长N对称中心在n0=(N-1)/2的矩形窗RN(n)与hd(n)相乘，截取出的实际系数序列h(n)，如下图,41,矩形窗截取实际滤波器系数,42,矩形窗截取的频率特性计算,上图为理想滤波器的脉冲响应hd(n)和其频谱；中图为矩形窗函数RN(n)和其频谱；下图时域为hd(n)和RN(n)的乘积，其频谱为上两频谱的卷积。出现波动。,43,矩形窗截取的频率特性计算,求理想脉冲函数加窗后的频谱特性(图7.3.1(f)可以利用频域卷积定理。时域相乘h(n)=hd(n)RN(n),对应于频域卷积，做卷积必须考虑函数的正负号，所以其频谱需要用符幅特性来表示：其中Hd(j)为hd(n)的频谱WRN(j)为RN(n)的频谱，用符幅特性表示。,44,矩形窗截取的频率特性计算,分析其幅度，可得截断后滤波器的频率特性这一积分的几何意义见图。实际的H()和理想的Hd()有两点差别：（1）通带内增加了波动，最大的峰值在=c-2/N处。阻带内产生了余振，最大的负峰在=c+2/N处；（2）在理想特性不连续点=c附近形成过渡带。过渡带的宽度近似等于ARN()主瓣宽度。,45,(a).理想滤波器频谱；(b).矩形窗函数符幅特性；此时相当于=0，它与(a)卷积的结果为矩形窗截断区的正面积减负面积，设为1。(c).=c时卷积输出应为阴影面积的积分，此处应为0.5。(d).c时少一块负面积，卷积输出大于1；(e)c时多一块负面积，卷积输出小于0。(f).卷积输出的波动波形。,46,矩形窗截取的实际特性,仔细地分析一下过渡带和波动的特性。我们关心的不是ARN()本身的波动，而是它与理想幅特性卷积后的H()在=c附近的波动。ARN()从-逐步右移，进入-cc积分区间，形成H()的过程实际上就是对ARN()进行积分，积分的结果称为累积符幅函数。因为不能用解析方法计算，这里就用MATLAB进行数值积分。所用的函数为y=cumtrapz(x)，它是按梯形积分方法编程的。程序hc731如下：,47,计算累积符幅特性的程序hc731,N=?;dx=0.01;%设置初值 x=-pi/2:dx:pi/2+1e-10;%自变量数组AR=sin(N*x)./sin(x);%矩形窗符幅特性H=cumtrapz(AR)*dx/pi;%数字积分subplot(2,1,1),plot(x,AR),%画符幅特性曲线subplot(2,1,2),plot(x,H)%画累积符幅特性曲线delta1=max(H)-1，%通带最大波动（绝对值）delta2=abs(min(H),%阻带最大波动（绝对值）%求过渡带宽度xtrans=dx*(find(H=max(H)-find(H=min(H),48,程序hc731中的x换成w,N=?;dw=0.01;%设置初值 w=-pi/2:dw:pi/2+1e-10;%自变量数组AR=sin(N*w)./sin(w);%矩形窗符幅特性H=cumtrapz(AR)*dw/pi;%数字积分subplot(2,1,1),plot(w,AR),%画符幅特性曲线subplot(2,1,2),plot(w,H)%画累积符幅特性曲线delta1=max(H)-1,%通带最大波动（绝对值）delta2=abs(min(H),%阻带最大波动（绝对值）%求过渡带宽度wtrans=dw*(find(H=max(H)-find(H=min(H),49,设计FIR滤波器的窗函数法,比较N=10及100的结果，可以看出N只影响过渡带宽度(分别为0.628和0.063)，不影响吉布斯效应波动的幅度(As分别为22及25dB)。,50,矩形窗截取的实际特性,(1).一般滤波器都要求过渡带愈窄愈好。过渡带与N成反比，减小过渡带的措施是增加N；(2).通带内的波动影响滤波器通带中的信号通过的均匀性；阻带内的波动使某些频率成分在阻带内不能充分抑制，可能使最小衰减不满足技术要求；上述分析说明，增加N对于通带波动和阻带衰减几乎没有影响。矩形窗带来的阻带衰减为21dB左右。要减小波动只有改变窗函数的形状。(3).增大N将增加滤波器的复杂性和计算时间。,51,几种其它窗函数的介绍,为了减少带内波动以及加大阻带的衰减，改变N是无效的，只有改变窗的形状。科学家们力求找到适当的窗函数形状，使其符幅函数的主瓣包含更多的能量，以减小其旁瓣幅度。旁瓣的减小可使通带、阻带波动减小，从而加大阻带衰减。但总是以加宽过渡带为代价的。下面介绍几种常用的窗函数wdx(n)，实际滤波器的系数向量h(n)等于理想滤波器脉冲响应hd(n)与窗函数wdx(n)相乘：h(n)=hd(n)wdx(n),52,几种其它窗函数的介绍,2.三角形窗(Bartlett Window）3.汉宁（Hanning）窗升余弦窗4哈明（Hamming）窗改进的升弦窗,53,MATLAB提供的其它窗函数,N=45;n=1:N;beta=5.658;wdrc=boxcar(N);wdbt=triang(N);wdhn=hanning(N);wdhm=hamming(N);wdbl=blackman(N);wdks=kaiser(N,beta);plot(n,wdrc,wdbt,wdhn,wdhm,wdbl,wdks)legend(rc,bt,hn,hm,bl,ks)%线型标注就可以得到综合表示几种窗函数形状的图,54,几种其它窗函数的介绍,55,几种其它窗函数的介绍,设计窗函数总的原则是使其符幅函数的主瓣包含更多的能量，以减小其旁瓣幅度。所以除了矩形窗具有正负变化的符幅特性以外，其它各种窗函数都没有负的幅特性，其幅特性就等于符幅特性，因此计算时可以调用信号处理工具箱中的freqz函数,进而通过求对数得到以分贝为单位的幅特性。例7.3.3 求出长度为51的布莱克曼窗的符幅特性，累积符幅特性，以及以分贝表示的幅特性和累积符幅特性曲线，用四个子图画在同一张图上。,56,三角窗函数的过渡带和波动,解；解此题的MATLAB程序核心语句如下：N=51;n=1:N;wd=triang(N);给出三角窗A,omega=amplres(wd);求符幅函数Ad=fliplr(A),A(2:end);扩展至负频率区omegad=-fliplr(omega),omega(2:end);Gd=cumtrapz(Ad)*(omega(2)-omega(1)/pi;AdB=20*log10(abs(Ad);将Ad转换为分贝GdB=20*log10(Gd);将Gd转换为分贝分别画出Ad,Gd,AdB,GdB四个变量随omegad的变化曲线如图,57,三角窗函数的过渡带和波动,三角窗的符幅特性A已经没有负值的部分了，所以其累积符幅特性呈单调上升的形状，看不到波动。,58,各种窗函数的过渡带和波动,59,设计FIR滤波器的窗函数法,用窗函数设计FIR滤波器的步骤如下：(1).根据对过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N。(查表7.3.1)(2).根据待求滤波器的理想频率响应求出单位脉冲响应hd(n)。（调用ideallp函数）function hd=ideallp(wc,N);%hd=点0 到 N-1之间的理想脉冲响应%wc=截止频率(弧度);N=理想滤波器的长度tao=(N-1)/2;n=0:(N-1);m=n-tao+eps;%加一个小数偏移hd=sin(wc*m)./(pi*m);,60,理想单位脉冲响应hd(n)求法,(3).计算滤波器的单位脉冲响应h(n)。它是理想脉冲响应和窗函数的乘积。h(n)=hd(n).*Wd(n)(4).验算技术指标是否合格。(调用myfreqz函数）function db,mag,pha,grd,w=myfreqz(b,a);H,w=freqz(b,a,1000,whole);H=(H(1:1:501);w=(w(1:1:501);以下为求db,mag,pha,grd的语句(5).如果不满足要求，可调整窗函数类型或长度，重复计算，直到满足要求为止。,61,设计FIR滤波器的窗函数法,例7.3.4 用矩形窗、汉宁窗和布莱克曼窗设计FIR低通滤波器。设N=11，c=0.2。解：因为题中已给出了窗函数类型及其长度，解此题时可以省去步骤(1)。直接从步骤(2)开始。由此编出的MATLAB程序hc734如下：wc=0.2*pi;N=11;hd=ideallp(wc,N)%用wc=0.2*pi的理想低通wd1=boxcar(N);b1=hd.*wd1;%用矩形窗设计：H1,w=freqz(b1,1);%检验矩形窗设计结果其他两种窗函数程序相仿。,62,设计FIR滤波器的窗函数法,三种窗函数设计的结果画在同一图上。,63,设计FIR滤波器的窗函数法,例7.3.5 根据下列技术指标,p=0.2,Rp=0.25dB;s=0.3,As=50 dB,设计一个数字FIR低通滤波器。并给出所设计的滤波器的频率响应图。解：从表中选出一个适当的窗函数，哈明窗和布莱克曼窗均可提供大于50dB的衰减。选择哈明窗，因为它有较小的过渡带，阶数较低。尽管在设计中，用不到通带波动值Rp=0.25dB，但必须检查设计的实际波动，验证它是否确实在给定容限内。程序hc735给出了设计的步骤。,64,窗函数法设计算例,wp=0.2*pi;ws=0.3*pi;deltaw=ws wp;%过渡带宽 N0=ceil(6.6*pi/deltaw);%计算滤波器长度N=N0+mod(N0+1,2);%确保其长度N为奇数wdham=(hamming(N);%求窗函数wc=(ws+wp)/2;%取均值为截止频率hd=ideallp(wc,N);%求理想脉冲响应h=hd.*wdham;设计的脉冲响应db,mag,pha,grd,w=myfreqz(h,1);以下为对设计结果进行检验和绘图语句,65,窗函数法设计算例,对程序hc735需作的说明：(1).N=N0+mod(N0+1,2);的作用是保证N为奇数，使滤波器始终保持类型I的性质。(2).freqz函数中把0到2分为1000份，频率分辨率是dw=2/1000，Rp=-(min(db(1:wp/dw+1)在下标1到通带下标wp/dw+1的范围内找分贝幅特性的最小值，即通带波动值；(3).As=-round(max(db(ws/dw+1:501)是在阻带下标ws/dw+1到最高数字频率pi下标(501)的范围内找阻带最大衰减值。计算结果：N67，阻带衰减As=52dB，达到要求。,66,窗函数法设计算例,67,窗函数法设计算例,例7.3.6 利用例给出的设计技术指标，选择凯泽窗，设计出所需的低通滤波器。解：设计步骤由程序hc736给出。它与hc735的差别仅在于选N0,beta和凯泽窗函数wdkai三条语句。计算结果为：凯泽窗的参数为N61，4.5513，实际的阻带衰减为52dB。时域和频域的曲线在图中给出。可见它比哈明窗具有较低的阶次，即可以用较少的运算时间，而得到同样的指标。,68,窗函数法设计算例,hc736 中的核心语句N0=ceil(As-7.95)/(14.36*deltaf)+1;%N=N0+mod(N0+1,2);%确保其长度N为奇数beta=0.1102*(As-8.7)%按(7.3.17)式算值wdkai=(kaiser(N,beta);%求凯泽窗wc=(ws+wp)/2;hd=ideallp(wc,N);%求理想脉冲响应h=hd.*wdkai;设计的脉冲响应检验语句同hc735最后带下划线的四条语句，绘出的曲线见下图。,69,窗函数法设计算例,70,窗函数法设计算例,例7.3.7 设计下面的数字带通滤波器。低端阻带边缘:1s=0.2,As=60 dB低端通带边缘:1p=0.35,Rp=1 dB高端通带边缘:2p=0.65,Rp=1 dB高端阻带边缘:2s=0.8,As=60 dB这些参数如图所示。,71,窗函数法设计算例,解：本题设计的是带通滤波器，有两点要注意:(1).理想脉冲响应如何确定。它可以由两个理想低通脉冲响应相减得到，如图所示。(2)存在两个过渡带，在窗设计中必须取得相同。程序hc737给出了设计步骤。,72,窗函数法设计算例,deltaw=min(wp1-ws1),(ws2-wp2);%求过渡带的小者N0=ceil(11*pi/deltaw);%按表求滤波器长度N0N=N0+mod(N0+1,2);%使长度N为奇数wdbla=(blackman(N);%求窗函数wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;%截止频率取通带阻带边缘频率的平均值hd=ideallp(wc2,N)-ideallp(wc1,N);%按图求带通滤波器理想脉冲响应h=hd.*wdbla;%求实际滤波器脉冲响应计算结果为：N61,As=75dB,Rp=0.003。,73,窗函数法设计算例,74,设计FIR滤波器的窗函数法,在信号处理工具箱中提供了窗函数方法的设计函数fir1。它的标准调用格式为：b=fir1(M,wn,type,window)其中 b 为待设计的滤波器系数向量，其长度为N=M+1M为所选的滤波器阶数 wn为滤波器给定的边缘频率，是标量或数组；type为滤波器的类型，如高通、带通、带阻等，window为选定的窗函数类型，缺省时为hamming窗。例如 b=fir1(32,0.7,high,kaiser(33,4)表示用=4的凯泽窗，设计阶数为32（长度为33），边界数字频率为0.7的高通滤波器。键入这个语句，MATLAB将给出设计得的滤波器的系数向量b。,75,设计FIR滤波器的频率样本法,窗函数法是从时域出发的一种设计方法。它从理想脉冲响应出发，进行截断或加权修改。最后使频率响应偏离了理想值。自然会产生一个问题，是不是可以先在频域修改滤波器的预期频率特性，使它接近于实际的需要和可能，然后直接用它的傅立叶反变换作为滤波器系数？答案是肯定的。预期频率特性不是以连续函数的形式，而是在一组离散的频点的样本值上给出，所以称为频率样本法。,76,设计FIR滤波器的频率样本法,在设定预期频率特性时应根据下列基本原则：(1).预期频率特性的样本点数应等于滤波器的长度N，并在单位圆上等间隔分布；(2).预期频率特性应具有共軛对称性，以保证其傅立叶反变换所得系数为实序列；故其幅特性应为偶函数，相特性应为奇函数。(3).预期频率特性的相位特性应该与频率成线性关系；这意味着其符幅特性及其反变换所得序列应该具有对称或反对称的特点。(4).必须恰好取够整个单位圆上的样本，才能正确地求得其h(n)=IDFTH(k).,77,设计FIR滤波器的频率样本法,以上的原则可用数学公式描述如下(见表7.2.1)。为了使h(n)为实序列，且满足h(n)=h(N-n-1)的第I类或第II类线性相位滤波器，H(k)的相特性(k)和符幅特性A(k)应满足如下条件：线性相位条件：符幅特性的对称条件：N为奇数时A(k)=A(N-k-1)N为偶数时A(k)=-A(N-k-1)注意不管符幅特性是对称或反对称，其幅特性|H(k)|都是对称的。,78,设计FIR滤波器的频率样本法,频率样本法的设计过程如下：设定了z域单位圆上的频率样本值如左图，用IDFT求其脉冲响应h(n)。而h(n)的DTFT就是整个单位圆上连续的频率特性，从左图求右图的过程可以保证样本点上的频谱合乎要求，其余频率上则有误差，偏差一般在频率变化最陡处最大。,79,频率样本法的设计算例,例 利用频率样本法设计线性相位低通滤波器，要求截止频率c=/3，采样点数取N=13和N=65作比较。解：用理想低通作为逼近滤波器。它在主值频域的幅特性及相特性如图实线所示，对它的13点采样情况如图中圆点所示。采样得到的数据为(用MATLAB语言表达)：A=ones(1,3),zeros(1,8),ones(1,2);theta=-pi*0:12*12/13;,80,频率样本法的设计算例,81,设计算例的核心语句,N=input(N=);wc=pi/3;N1=fix(wc/(2*pi/N);N2=N-2*N1-1;%样本间距为2*pi/N,通带长度N1,阻带长度N2A=ones(1,N1+1),zeros(1,N2),ones(1,N1);%符幅theta=-pi*0:N-1*(N-1)/N;%相位样本序列H=A.*exp(j*theta);%频率特性样本序列h=real(ifft(H);%反变换求出系数序列db,mag,pha,grd,w=myfreqz(h,1);%检验下为检验和绘图语句,82,频率样本法的设计算例,令N=13，执行hc741所得的曲线。分别给出了(a)幅特性样本、(b)设计出的滤波器系数、(c)滤波器符幅特性，和(d)滤波器对数特性。可见样本点处的幅特性都等于给定值1或0。,83,频率样本法的设计算例,N=13和N=65算出的两种幅特性同画于下图，可知,(1).N愈大，特性愈接近理想情况;(2).N愈大，过渡带宽愈小。但Rp和As没有改善。(3).要想改善Rp和As,最好在过渡带中加一个点，使过渡平缓。,84,设计FIR滤波器的频率样本法,优化设计法 为了得到更大的衰减，必须增大N，并且使过渡带中的样本成为自由样本，也就是说，可以改变它们的样本值，以便在给定的N和过渡带宽下得到最大的衰减。这个问题称为优化问题，利用线性规划技术可以解决它。若不想涉及优化的数学问题，也可以用试凑方法来解决，下面的例子可以说明改变过渡带样本值对设计的影响。,85,频率样本法的设计算例,用加宽过渡带的方法，改进例的设计。解：先把过渡带增加一个样本，把阻带边界上的脉冲由零改为T1。符幅特性赋值语句就成为：A=ones(1,N1+1),T1,zeros(1,N2-2),T1,ones(1,N1);同时，把ws1向右移位：ws1=wp1+2*2*pi/N;加一条对T1的赋值语句，写成键盘输入形式：T1=input(过渡带中单个样本值T1);其它语句不变。运行此程序，并依次输入T1=0.5,0.4,0.3可得到表前三列的结果：,86,频率样本法的设计算例,利用hc742a，输入两位小数点的T1,经过试凑，可知T1取0.38时，得到最大的阻带衰减44分贝。若在过渡带内放两个样本T1,T2，相应地改变符幅序列的赋值语句，同时改变通带和阻带的频率范围语句。还要加两条分别对T1，T2的赋值语句。构成的程序为hc742b。,87,频率样本法的设计算例,从此表可以看到T1=0.59，T2=0.11是一个局部最优的组合，此时可以得到最大的阻带衰减65分贝。这个条件下画出的与图相应的曲线见图。,运行程序hc742b并输入两位数的T1，T2进行试凑，结果为：,88,频率样本法的设计算例,过渡带中放两个样本的设计结果。可以看出，适当放宽了过渡带，带来的好处是大大降低了阻带衰减，滤波器更接近于理想。,89,设计FIR滤波器的频率样本法,MATLAB提供的频率样本法的设计函数fir2。它的典型调用方法为：h=fir2(M,f,A)其中M为FIR数字滤波器的阶数（滤波器的长度为NM1），数组f和A给出它的预期频率响应。算出的滤波器系数为h，其 长度为N。数组f包含各边缘频率，单位为，即0.0f1.0。数组A为各指定频率上预期的幅度响应。fir2函数缺省地使用hamming窗函数。,90,设计FIR滤波器的频率样本法,fir2函数的一个很有用的功能是可以按指定频谱形状，生成信号序列，见下例。例7.4.3 要求生成一个限带序列h(n)，其长度为51，它的频谱应该如图上图。解：先按预期的频率特性写出f,A向量，然后调用fir2函数，用下列程序hc743来完成。N=input(N=);M=N-1;%给定序列长度f=0,0.45,0.5,1;A=0,1,0,0;%预期频谱h=fir2(M,f,A);%生成脉冲响应序列H,w=freqz(h,1);%计算脉冲序列的频谱,91,设计FIR滤波器的频率样本法,1。预期的频谱形状2。用fir2函数设计出的脉冲序列；3。此脉冲序列的频谱,92,7.5 FIR滤波器的最优设计法,窗函数和频率样本两种方法存在一些不足：第一，在设计过程中，不能精确指定通带和阻带频率s和p，只能接受大体合用的值；第二，不管是在窗设计方法中使12，还是在频率样本法中只优化2，均不能同时控制波动系数1和2。第三，理想响应和实际响应之间的逼近误差，在全频带区间上不是均匀分布的。如果使误差均匀分布，就可以得到一个较低阶滤波器。最优化是能克服上述问题的数学方法。它用的数学较深，不过现在可藉计算机之助简化。,93,FIR滤波器的最优设计法,最优化的思想实际上在窗函数法中就有反映。用E()表示理想滤波器Hd()和所设计滤波器H()之间的幅频特性误差，若以这个误差平方的连续积分作为最优化的目标函数，则可以证明，采用矩形窗时，均方误差是最小的，也就是说，矩形窗是连续最小均方误差的最优设计法。不过，这种方法并未得到应用。因吉布斯效应，平方积分误差最小不能保证误差幅度最小。,94,FIR滤波器的最优设计法,本节将介绍两种最优化方法，第一种是离散最小平方积分法，或称离散最小二乘法。思路是使在给定的一些离散频点上，使实际幅频特性与理想幅频特性之间的误差的平方和为最小；第二种是最小最大波动法，也称为等波动法。具有这种性质的滤波器幅特性在通带和阻带上的误差峰值是均匀分布的。其误差具有等波动特性。因而可把波动的幅度控制到最小。,95,FIR滤波器的最优设计法,离散最小二乘法：从频率样本出发，分析实际符幅特性与给定的（理想的）幅特性之间的误差，求出它们在若干样本点的误差平方和（与上述e2的不同就在于用有限个样本点误差平方之和而不是对