数字信号处理（曹成茂）第五章 有限长单位脉冲响应（fir）滤波器的设计方法.ppt

第五章 有限长单位脉冲响应（FIR）滤波器的设计方法,序言 §5.1 线性相位FIR数字滤波器的特性,§5.2 窗口设计法（时间窗口法）,§5.3 频率采样法,§5.4 FIR数字滤波器的最优化设计 §5.5 IIR与FIR数字滤器的比较,序言,FIR数字滤波器的差分方程描述 ,对应的系统函数,因为它是一种线性时不变系统，可用卷积和形式表示 ,比较、得：,FIR数字滤波器的特点(与IIR数字滤波器比较)： 优点 ：（1）很容易获得严格的线性相位，避免被处理 的信号 产生相位失真，这一特点在 宽频带信 号处理、阵 列信号处理、数据传输等系统中 非常重要； （2 ）可得到多带幅频特性； （3 ）极点全部在原点（永远稳定），无稳定 性问题； （4 ）任何一个非因果的有限长序列，总可以通过一 定的延时，转变为因果序列， 所以因果性总是 满足； （5）无反馈运算，运算误差小。,缺点：（1）因为无极点，要获得好的过渡带特性，需以较 高的阶数为代价； （2）无法利用模拟滤波器的设计结果，一般无解 析设计 公式，要借助计算机辅助设计程序完成。,§5.1 线性相位FIR数字滤波器的特性,5.1.1 线性相位的条件 线性相位意味着一个系统的相频特性是频率的线性函数，即,式中为常数，此时通过这一系统的各频率分量的时延为一相同的常数，系统的群时延为,线性相位FIR滤波器的DTFT为,式中 H()是正或负的实函数。等式中间和等式右边的实部与虚部应当各自相等，同样实部与虚部的比值应当相等:,将上式两边交叉相乘，再将等式右边各项移到左边，应用三角函数的恒等关系,满足上式的条件是,另外一种情况是，除了上述的线性相位外，还有一附加的相位，即,利用类似的关系，可以得出新的解答为,偶对称,奇对称,图1 线性相位特性,分四种情况,5.1.2 线性相位FIR滤波器的幅度特性,分四种情况 1 h(n) 偶对称，N为奇数 h(n)=h(N-1-n),5.1.2 线性相位FIR滤波器的幅度特性,令 ，则,令,则,由于 偶对称，因此 对这些频率也呈偶对称。,2h(n)偶对称，N为偶数 h(n)=h(N-1-n),令 ，则,或写为：,由于 奇对称，所以 对 也为奇对称，且由于 时， 处必有一零点，因此这种情况不能用于设计 时 的滤波器，如高通、带阻滤波器。,3. h(n)奇对称，N为奇数，h(n)=-h(N-1-n),令 n=m+(N-1)/2，得：,所以,由于 点呈奇对称，所以 对这些点也奇对称。 由于 时， 相当于H（z）在 处有两个零点，不能用于 的滤波器设计，故不能用作低通、高通和带阻滤波器的设计。,4.h(n)奇对称，N为偶数,令,由于 在=0,处为零，所以H()在=0, 2处为零，即H(z)在z=1上有零点，并对=0,2呈奇对称。,四种线性相位FIR滤波器,四种线性相位FIR滤波器,相位特性,辅助序列,适用,幅度特性及特点,N为奇数,N为偶数,N为奇数,N为偶数,低通 高通 带通 带阻,低通 带通,带通,高通 带通,四种线性相位FIR DF特性，参考表5.1 第一种情况 ，偶、奇，四种滤波器都可设计。,第二种情况，偶、偶，可设计低、带通滤波器，不能设计 高通和带阻。 第三种情况，奇、奇，只能设计带通滤波器，其它滤波器 都不能设计。 第四种情况，奇、偶，可设计高通、带通滤波器，不能设 计低通和带阻。,例1 N=5, h (0) = h (1) = h (3) = h (4) = -1/2, h (2) = 2，求幅度函数H ()。 解 为奇数并且h(n)满足偶对称关系 a (0) = h (2) = 2 a (1) = 2 h (3) = -1 a (2) = 2 h (4) = -1 H () = 2 - cos- cos2 = 2- (cos+cos2),小结：,四种FIR数字滤波器的相位特性只取决于h(n)的对称性，而与h(n)的值无关。 幅度特性取决于h(n)。 设计FIR数字滤波器时，在保证h(n)对称的条件下，只要完成幅度特性的逼近即可。 注意：当H()用H()表示时，当H()为奇对称时，其相频特性中还应加一个固定相移。,5.1.3 线性相位FIR滤波器的零点特性,由该式可看出，若z=zi是H（z）的零点，则z=z-1i也一定是H（z）的零点。由于h(n)是实数，H（z）的零点还必须共轭成对，所以z=z*i 及 z=1/z*也必是零点。 所以线性相位滤波器的零点必须是互为倒数的共轭对，即成四出现，这种共轭对共有四种 可能的情况： 既不在单位园上，也不在实轴上，有四个互为倒数的两组共轭 对， zi z*i 1/zi 1/z*i 图4.2(a) 在单位圆上，但不在实轴上，因倒数就是自己的共轭，所以有一对共轭零点， zi,z*i 图4.2（b） 不在单位圆上，但在实轴上，是实数,共轭就是自己，所以有一对互为倒数的零点, zi, 1/zi 图4.2（c） 又在单位圆上，又在实轴上，共轭和倒数都合为一点，所以成单出现，只有两种可能， zi=1或zi=-1 图4.2(d),同样道理，对于第三种FIR滤波器， h(n)奇对称，N为奇数，因 所以z=1，z=-1都是H（z）的单根； 对于第四种滤波器， h(n)奇对称，N为偶数，H（O）=0，所以z=1是H（z）的单根。 所以，h(n)奇对称H(0)=0 N为偶数H()=0 线性相位滤波器是FIR滤波器中最重要的一种，应用最广。实际使用时应根据需用选择其合适 类型，并在设计时遵循其约束条件。,我们从幅度响应的讨论中已经知道，对于第二种FIR滤波器（h(n)偶对称，N为偶数）， ， 即 是 的零点，既在单位圆，又在实轴，所以，必有单根；,§5.2 窗口设计法（时域）,如果希望得到的滤波器的理想频率响应为 ，那么 FIR滤波器的设计就在于寻找一个系统函数，频率响应 去逼近 ，逼近方法有三种： 窗口设计法（时域逼近） 频率采样法（频域逼近） 最优化设计（等波纹逼近） 时间窗口设计法是从单位脉冲响应序列着手，使h(n)逼近理想的单位脉冲响应序列hd(n)。我们知道hd(n)可以从理想频响 通过傅里叶反变换获得,但一般来说，理想频响 是分段恒定，在边界频率处有突变点，所以，这样得到的理想单位脉冲响应hd(n)往往都是无限长序列，而且是非因果的。但FIR的h(n)是有限长的，问题是怎样用一个有限长的序列去近似无限长的hd(n)。最简单的办法是直接截取一段 hd(n) 代替 h(n) 。这种截取可以形象地想象为h(n)是通过一个“窗口”所看到的一段hd(n)，因此 ，h(n)也可表达为h(n)和一个“窗函数”的乘积，即 h(n)=w(n) hd(n) 在这里窗口函数就是矩形脉冲函数RN（n），当然以后我们还可看到，为了改善设计滤波器的特性，窗函数还可以有其它的形式，相当于在矩形窗内对hd(n)作一定的加权处理。,设计步骤：,1）由定义,3）卷积,插值,一.矩形窗口法,则,以一个截止频率为 c的线性相位理想低通滤波器为例,讨论FIR的设计问题。 a. 对于给定的理想低通滤波器 ，计算,：低通滤波器的延时,理想特性的hd(n)和Hd(),这是一个以 为中心的偶对称的无限长非因果序列，如果截取一段n=0N-1的hd(n)作为h(n)，则为保证所得到的是线性相位FIR滤波器，延时 应为h(n)长度N的一半,即,其中,b.计算,c.计算 。 设 为窗口函数的频谱: 用幅度函数和相位函数来表示，则有,其线性相位部分 则是表示延时一半长度,对频响起作用的是它的幅度函数,理想频响也可以写成幅度函数和相位函数的表示形式 Hd(ej)=Hd()e-j 其中幅度函数为 两个信号时域的乘积对应于频域卷积，所以有,如果也以幅度函数和相位函数来表示 H（ej）, 则实际FIR滤波器的幅度函数H（）为 正好是理想滤波器幅度函数与窗函数幅度函数的卷积。,4个特殊频率点看卷积结果：,（1）=0时, H(0)等于,在-c, c内的积分面积,因一般,故H(0)近似为,在-, 内的积分面积,（2）=c时，一半重叠， H(c)=0.5 H(0);,（3） =c 2/N时，第一旁瓣（负数）在通带外，出现正肩峰；,（ 4） =c +2/N 时，第一旁瓣（负数）在通带内，出现负肩峰。,矩形窗的卷积过程,窗口函数对理想特性的影响： 改变了理想频响的边沿特性，形成过渡带，宽为 ， 等于WR()的主瓣宽度。（决定于窗长） 过渡带两旁产生肩峰和余振（带内、带外起伏），取决于 WR()的旁瓣，旁瓣多，余振多；旁瓣相对值大，肩峰强 ，与 N无关。（决定于窗口形状） N增加,过渡带宽减小,肩峰值不变。 因主瓣附近 其中x=N/2,所以N的改变不能改变主瓣与旁瓣的比例关系，只能改变WR（）的绝对值大小和起伏的密度，当N增加时，幅值变大，频率轴变密，而最大肩峰永远为8.95%，这种现象称为吉布斯（Gibbs）效应。,改变窗函数的形状，可改善滤波器的特性，窗函数有 许多种，但要满足以下两点要求：,窗谱主瓣宽度要窄，以获得较陡的过渡带； 相对于主瓣幅度，旁瓣要尽可能小，使能量尽量集中在主瓣中，这样就 可以减小肩峰和余振，以提高阻带衰减和通带平稳性。 但实际上这两点不能兼得，一般总是通过增加主瓣宽度来换取对旁瓣的抑制。,肩峰值的大小决定了滤波器通带内的平稳程度和阻带内的衰减，所以对滤波器的性能有很大的影响。,5.2.2 几种常用的窗函数 1. 矩形窗，上面已讲过，不再细述 2. 汉宁窗（升余弦窗） 利用的移位特性，汉宁窗频谱的幅度函数W（）可用矩形窗的幅度函数表示为：,三部分矩形窗频谱相加，使旁瓣互相抵消，能量集中在主瓣，旁瓣大大减小，主瓣宽度增加1倍，为 。,0,-80,-60,-44,-20,0,矩形窗,Hanning窗,dB,3. 汉明窗（改进的升余弦窗） 它是对汉宁窗的改进，在主瓣宽度（对应第一零点的宽度）相同的情况下，旁瓣进一步减小，可使99.96%的能量集中在窗谱的主瓣内。 4. 布莱克曼窗（三阶升余弦窗） 增加一个二次谐波余弦分量，可进一步降低旁瓣，但主瓣宽度进一步增加，为 。增加N可减少过渡带。 频谱的幅度函数为：,窗口函数的频谱 N=51，A=20lg|W()/W(0)|,四种窗函数的比较,(a)矩形窗；(b) 汉宁窗；(c)汉明窗；(d)布莱克曼窗,几种窗函数的基本参数,5.凯塞窗 以上四种窗函数，都是以增加主瓣宽度为代价来降低旁瓣。凯塞窗则可自由选择主瓣宽度和旁瓣衰减。,I0(x)是零阶修正贝塞尔函数，参数可自由选择，决定主瓣宽度与 旁瓣衰减。越大，w(n)窗越窄，其频谱的主瓣变宽，旁瓣变小。一般取 49。,=5.44 接近汉明 =8.5 接近布莱克曼 =0 为矩形,图2 凯塞窗函数,图1 零阶修正贝塞尔函数,I0(x),x,0,1,不同 值凯塞窗特性,当MN时， hM(n)hd(n) 零阶贝塞尔函数,窗口设计法的主要工作是计算hd(n)和w(n)，当 较为复杂时，hd(n)不容易由反傅里叶变换求得。这时一般可用离散傅里叶变换代替连续傅里叶变换，求得近似值： 令,过渡带宽 At阻带最小衰减,在窗口函数方面，MATLAB可计算以下的窗口函数： W=hanning(N) W=hamming(N) W=blackman(N) W=kaiser(N, beta) 变量是窗函数的长度，beta就是凯塞窗的参数，它可以利用以上讨论的公式计算，是返回的窗口函数，长度为。 对于凯塞窗，MATLAB还提供了一个根据滤波器设计指标计算窗口函数的参数的函数文件，其调用形式如下： M, Wc, beta, ftype=kaiserord(f, a, dev, fs) f为通带和阻带边界频率，最高为fs/2；a为相应频带的幅度值；dev为波动值；fs缺省为2。M是滤波器阶数；Wc边界频率；beta即；ftype是滤波器类型，低通为low。,§5.3 频率采样法,工程上，常给定频域上的技术指标，所以采用频域设计更直接。 一、基本思想 使所设计的FIR数字滤波器的频率特性在某些离散频率点上的值准确地等于所需滤波器在这些频率点处的值，在其它频率处的特性则有较好的逼近。,内插公式,逼近误差,由 得到了H(z) 或 。 要讨论 与 的逼近程度 ，以及 与H（k）的关系？ 由,令,，则,单位圆上的频响为：,这是一个内插公式。,式中 为内插函数 令 则,内插公式表明： 在每个采样点上， 逼近误差为零，频响 严格地与理想频响的采样值 H(k)相等；,在采样点之间，频响由各采样点的内插函数延伸迭加而形成，因而有一定的逼近误差，误差大小与理想频率响应的曲线形状有关，理想特性平滑，则误差小；反之，误差大。在理想频率响应的不连续点附近， 会产生肩峰和波纹。 N增大，则采样点变密，逼近误差减小。,二.设计方法,1）确定,并得到,3）计算,2）计算,三、 约束条件,为了设计线性相位的FIR滤波器，采样值 H（k）要满足一定的约束条件。 前已指出,具有线性相位的FIR滤波器，其单位脉冲响应h(n)是实序列，且满足 ，由此得到的幅频和相频特性，就是对H(k)的约束。（表4.1）。 例如，要设计第一类线性相位FIR滤波器，即N为奇数，h(n)偶对称，则 幅度函数H（）应具有偶对称性：,令 则 必须满足偶对称性： 而 必须取为：,同样，若要设计第二种线性相位FIR滤波器，N为偶数，h(n)偶对称，由于幅度特性是奇对称的，,因此，Hk 也必须满足奇对称性： 相位关系同上， 其它两种线性相位FIR数字滤波器的设计，同样也要满足幅度与相位的约束条件。,例：设计一个FIR数字 LP 滤波器，其理想特性为 采样点数 N=33，要求线性相位。 解：根据表5.1，能设计低通线性相位数字滤波器的只有1、2两种，因N为奇数，所以只能选择第一种。 即 h(n)=h(N-1-n)， 幅频特性关于偶对称，也即 HK 偶对称。 利用 HK 的对称性,求2区间的频响采样值。,根据指标要求，在02内有33个取样点，所以第k点对 应频率为 而截止频率 0.5位于 之间 ，所以，k=08时，取样值为1；根据对称性， 故 k=2532时，取样值也为1，因 k=33 为下一周期，所以0区间有9个值为 1的采样点，2区间有8个值为 1 的采样点，因此：,N1=33; k1=0:(N1-1)/2; Wm1=2*pi*k1./N1; Ad1(1:(N1+1)/2)=1; Ad1(10:17)=0; Hd1=Ad1.*exp(-j*0.5*(N1-1)*Wm1); Hd1=Hd1 conj(fliplr( Hd1(2:(N1+1)/2) ) ); h1=real(ifft(Hd1); w1=linspace(0,pi-0.1,1000); H1=freqz(h1,1,w1); plot(w1/pi,20*log10(abs(H1);grid; axis(0 1 -100 20); xlabel('pi'); ylabel('幅度/db'); title('过渡带不设采样点');,从图上可以看出，其过渡带宽为一个频率采样间隔 2/33，而最小阻带衰减略小于20dB。 对大多数应用场合，阻带衰减如此小的滤波器是不能令人满意的。 增大阻带衰减三种方法： 1）加宽过渡带宽，以牺牲过渡带换取阻带衰减的增加。 例如在本例中可在k=9和k=24处各增加一个过渡带采样点H9=H24=0.5，使过渡带宽增加到二个频率采样间隔4/33，重新计算的H(ej)见图，其阻带衰减增加到约 -40dB。,N1=33; k1=0:(N1-1)/2; Wm1=2*pi*k1./N1; Ad1(1:(N1+1)/2)=1; Ad1(11:17)=0; Ad1(10)=0.5; Hd2=Ad1.*exp(-j*0.5*(N1-1)*Wm1); Hd2=Hd2 conj(fliplr( Hd2(2:(N1+1)/2) ) ); h2=real(ifft(Hd2); w1=linspace(0,pi-0.1,1000); H2=freqz(h2,1,w1); plot(w1/pi,20*log10(abs(H2);grid; axis(0 1 -100 20); xlabel('pi'); title('过渡带设一个采样点');,2）过渡带的优化设计 根据H(ej)的表达式，H(ej)是Hk的线性函数，因此还可以利用线性最优化的方法确定过渡带采样点的值，得到要求的滤波器的最佳逼近（而不是盲目地设定一个过渡带值）。 例如，本例中可以用简单的梯度搜索法来选择H9、H24,使通带或阻带内的最大绝对误差最小化。,要求使阻带内最大绝对误差达到最小（也即最小衰减达到最大），可计算得H9=0.3904。对应的 H(ej)的幅频特性，比H9=0.5时 的阻带衰减大大改善,衰减接近-50dB。如果还要进一步改善阻带衰减，可以进一步加宽过渡区，添上第二个甚至第三个不等于0的频率取样值，当然也可用线性最优化求取这些取样值。,不同过渡带采样值的频率响应,3）增大N 如果要进一步增加阻带衰减，但又不增加过渡带宽，可增加采样点数N。 例如，同样边界频率c=0.5 , 以N=65采样，并在k=17和k=48插入由阻带衰减最优化计算得到的采样值H17=H48=0.5886,在k=18、47处插入经阻带衰减最优化计算获得的采样值H18=H47=0.1065 , 这时得到的 H(ej)， 过渡带为6/65，而阻带衰减增加了20多分贝，达-60dB以上，当然，代价是滤波器阶数增加，运算量增加。,N=65;k=0:(N-1)/2; Wm=2*pi*k./N; Ad(1:(N+1)/2)=1; Ad(18)=0.5886;Ad(19)=0.1065;Ad(20:33)=0; Hd=Ad.*exp(-j*0.5*(N-1)*Wm); Hd=Hd conj(fliplr( Hd(2:(N+1)/2) ) ); h=real(ifft(Hd); w=linspace(0,pi-0.1,1000); H=freqz(h,1,w); plot(w/pi,20*log10(abs(H);grid; axis(0 1 -100 20); xlabel('pi'); title('过渡带设两个采样点，总采样点数提高一倍');,例 用频率采样法设计一个线性相位FIR带通滤波器，设N=32，理想频率特性为,解 N=32为偶数，按第二类线性相位FIR滤波器设计。 频率间隔为,上边界点k在9和10之间；下边界点k在3和4之间。,MATLAB程序如下: N=32; Hk=zeros(1,4) ones(1,6) zeros(1,13) -ones(1,6) zeros(1,3); k=0:N-1; hn=real(ifft(Hk.*exp(-j*pi*(N-1)*k/N); H w=freqz(hn, 1); plot(w/pi, 20*log10(abs(H); axis(0 1 -60 10); grid; xlabel('归一化频率/pi') ylabel('幅度/dB'),小结： 频率采样设计法优点： 直接从频域进行设计，物理概念清楚，直观方便； 适合于窄带滤波器设计，这时频率响应只有少数几个非零值。 典型应用：用一串窄带滤波器组成多卜勒雷达接收机，覆盖不同的频段，多卜勒频偏可反映被测目标的运动速度； 缺点：截止频率难以控制。 因频率取样点都局限在2/N的整数倍点上，所以在指定通带和阻带截止频率时，这种方法受到限制，比较死板。 充分加大N，可以接近任何给定的频率，但计算量和复杂性增加。,§5.4 FIR滤波器的最优化设,前面介绍了FIR数字滤波器的两种逼近设计方法，即窗口法（时域逼近法）和频率采样法（频域逼近法），用这两种方法设计出的滤波器的频率特性都是在不同意义上对给定理想频率特性Hd(ej)的逼近。 说到逼近，就有一个逼近得好坏的问题，对“好”“坏”的恒量标准不同，也会得出不同的结论，我们前面讲过的窗口法和频率采样法都是先给出逼近方法，所需变量，然后再讨论其逼近特性，如果反过来要求在某种准则下设计滤波器各参数，以获取最优的结果，这就引出了最优化设计的概念，最优化设计一般需要大量的计算，所以一般需要依靠计算机进行辅助设计。,最优化设计的前提是最优准则的确定，在FIR滤波器最优化设计中，常用的准则有 最小均方误差准则 最大误差最小化准则。 1) 均方误差最小化准则， 若以E(ej)表示逼近误差，则 那么均方误差为,均方误差最小准则就是选择一组时域采样值，以使均方误差 ，这一方法注重的是在整个-频率区间内总误差的全局最小，但不能保证局部频率点的性能，有些频率点可能会有较大的误差，对于窗口法FIR滤波器设计，因采用有限项的h(n)逼近理想的hd(n)，所以其逼近误差为：,如果采用矩形窗,可以证明，这是一个最小均方误差。 所以，矩形窗窗口设计法是一个最小均方误差FIR设计，根据前面的讨论，我们知道其优点是过渡带较窄，缺点是局部点误差大,或者说误差分布不均匀。 2) 最大误差最小化准则（也叫最佳一致逼近准则）,表示为 其中F是根据要求预先给定的一个频率取值范围，可以是通带，也可以是阻带。最佳一致逼近即选择N个频率采样值 （ 或时域 h(n) 值 ），在给定频带范围内使频响的最大逼近误差达到最小。也叫等波纹逼近。 优点：可保证局部频率点的性能也是最优的，误差分布均匀， 相同指标下，可用最少的阶数达到最佳化。,例如，我们提到的频率采样最优化设计，它是从已知的采样点数N、预定的一组频率取样和已知的一组可变的频率取样（即过渡带取样）出发，利用迭代法（或解析法）得到具有最小的阻带最大逼近误差（即最大的阻带最小衰减）的FIR滤波器。但它只是通过改变过渡带的一个或几个采样值来调整滤波器特性。如果所有频率采样值（或FIR时域序列h(m)）都可调整，显然，滤波器的性能可得到进一步提高。,低通滤波器的误差分配,5.4.1 非线性最优法 用等波纹逼近法设计滤波器需要确定五个参数： M、c、r、1、2 按上图所示的误差容限设计低通滤波器，就是说要在通带 0 p 内以最大误差 1 逼近1，在阻带r 内 以最大误差2逼近零。 要同时确定上述五个参数较困难。常用的两种逼近方法： 1）给定M、1、2，以c和r为变量。 缺点：边界频率不能精确确定。 2）给定M、c和r，以1和2为变量，通过迭代运算 ，使逼近误差1和2 最小，并确定h(n)切比雪 夫最佳一致逼近。 特点：能准确地指定通带和阻带边界频率。,如上图所示的等波动低通滤波器，它对于给定的c和r值有误差值1、2。,在0 内最多有（M+1）个极值。,写出2M个方程，求解出（M+1）个滤波器系数和（M-1）个极值点频率（另外两个极值频率点为=0和= ）,5.4.2 插值解法,插值解法利用迭代算法避免了求一系列的非线性方程。 步骤如下： 1、选定通带的极值点数目和租代的极值点数目，并估计出极值点的频率。 2、利用标准的拉格朗日插值法，计算幅频连续特性的三角多项式，得出这个连续的幅频特性曲线。 3、用新的极值频率点取代原有的极值频率点，重复第二步。,5.4.3 雷米兹交替算法 误差函数 定义逼近误差函数：,为所设计的滤波器与理想滤波器的幅频特性在通带和阻带内的误差值， 是已知的权函数，在不同频带可取不同的值， 所要设计的滤波器的幅频特性 理想滤波器的幅频特性,交替定理（最佳逼近定理） 令F表示闭区间 的任意闭子集，为了使 在 F 上唯一最佳地逼近于 ,其充分必要条件是误差函数 在 F 上至少应有（M+2）次“交替”， 即 其中 ，且 属于F。 1） 至少有 M+2 个极值，且极值正负相间，具有等波纹的性质 ， 2）由于 是常数，所以 的极值也就是 的极值。,逼近方法：固定 k、M、 和 ，以 作为参变量。按照交替定理，如果 F 上的M+2个极值点频率 已知，则由（1）式可得到 M+2 个方程：,为极值点频率对应的误差函数值,注意：极值点频率必须位于 和 区间内。由于 和 固定，因而 和 必为这些极值频率中的一个，设 ，则应有 求解上述方程组可得到全部系数 问题：1）实际情况下，M+2 个极值点频率未知； 2）直接求解上述非线性方程组比较困难。,雷米兹（Remez）算法给出了求解切比雪夫最佳一致逼近问题的方法。,雷米兹（Remez）算法,1）在频率子集 F 上均匀等间隔地选取 M+2 个极值点频率,2）由 求 和 利用重心形式的拉格朗日插值公式，,其中,如在频带 F 上，对所有频率都有 ，则 为所求， 即为极值点频率。,3）对上次确定的极值点频率 中的每一点，在其附近检查是否存在某一频率处有 ， 如有，则以该频率点作为新的局部极值点。对 M+2 个极值点频率依次进行检查，得到一组新的极值点频率。重复步骤1）、2），求出 ，完成一次迭代。 重复上述步骤，直到 的值改变很小，迭代结束，这个 即为所求的 最小值。由最后一组极值点频率求出 ，反变换得到 , 完成设计。 优点： 可准确确定； 逼近误差均匀分布，相同指标下，滤波器所需阶数低。,有一些估算公式可用于决定最佳滤波器长度N：,对于窄带低通滤波器，对滤波器长度N起主要作用：,MATLAB实现,M,fo,ao,w = remezord(f,a,dev,fs) h = remez(M,fo,ao,w) remezord估计滤波器阶数N。remez完成等波纹FIR滤波器设计。,f：如FIR滤波器的B个频带分别为 0fpf1， f2ff3，f2(B-1)-2 f f2(B-1)-1，f2B-2 f f s/2 则f是一个有2B-2个元素的向量(B2)，其值为 f=f1, f2,f3, , f2B-2,a：是一个B个元素的向量,分别表示FIR滤波器在B个频带中的幅度值。一般对通带取值为1，阻带取值为0。,dev：是一个B个元素的向量,分别表示FIR滤波器在B个频带中的波动值。fs: 采样频率,返回参数：fpts，mag为B个频带的2B个边界频率和幅度值。wt为B表示各频带的加权值，缺省时表示各频带的加权值相同。,以上MATLAB程序设计的低通滤波器的频率响应，最小阻带衰减还不到40dB，为满足实际应用的要求可适当增加滤波器的阶数，现取N=33。,例：设计满足下列指标的等波纹线性相位FIR带通滤波器。 fs1=0.2p，fp1=0.3p, fp2=0.6p, fs2=0.7p, d1=0.1，d2=0.01,Fs1=0.2;Fp1=0.3;Fp2=0.6;Fs2=0.7; f=Fs1 Fp1 Fp2 Fs2;a=0 1 0; Rp=0.1;Rs=0.01;dev=Rs Rp Rs; N,fo,ao,w = remezord(f,a,dev); h = remez(N,fo,ao,w);,0,0.2,0.3,0.6,0.7,1,-40,-20,0,Normalized frequency,G,a,i,n,d,b,M=25时（remezord 确定的滤波器阶数）滤波器的增益响应。,0,0.2,0.3,0.6,0.7,1,-40,-20,0,Normalized frequency,G,a,i,n,d,b,M=30时滤波器的增益响应,§5.5 IIR与FIR数字滤器的比较,低,