一种用Remez函数设计FIR数字滤波器的设计方法

摘 要滤波是信号处理中最基本又极为重要的技术，利用滤波器技术可以从复杂的信号中提取出所需要的信号，抑制不需要的信号。绝大多数传感器输出的信号，在使用过程中，都必须进行滤波。所以滤波器是具有一定传输选择特性的、对信号进行加工处理的装置、它允许输入信号中的一些成分通过，抑制或衰减另一些成分。其功能是将输入信号变换为人们所需要的输入信号。在这个极具挑战的时代中，把计算机充分运用到教学及工程计算中，显然具有重要的意义，而随着计算机技术的发展,计算机软件在工程设计领域的应用越来越广。本文介绍了一种用Remez函数设计FIR数字滤波器的设计方法，并利用MATLAB强大的功能进行仿真。仿真结果表明，设计结果的各项性能指标均达到较高要求，并且该方法是一种最优化的设计方法。在同样的技术指标下，用这种方法得到的滤波器要比窗函数法、频率采样法得到滤波器的长度均要小，并且克服了窗函数设计法和频率采样法的缺点，使最大误差最小化，并在整个逼近频段上均匀分布。关键词：FIR数字滤波器；窗函数设计法；频率采样法；等波纹最优化设计法；MATLAB；Remez函数AbstractFiltering is the most basic and extremely important technology in signal processing . Use Filter technology,and you can extract the signal from the complex signals and suppress unwanted signals. The vast majority of the sensor output signals, in the course, must be filtered.Thus,Filter is a signal processing device with a certain transmission select features,which allows some components of the input signals through,inhibition or attenuation of some other components. Its function is to transform the input signals to the required input signals.In this very challenging times, the computer fully applied to the teaching and engineering calculations,is clearly of great significance.With the development of computer technology, computer software is more widely applied in the field of engineering design.This article describes a Remez function of FIR digital filter design methods, and simulation using MATLABs powerful features. The simulation results show that the design results of the performance indicators to meet the higher requirements, and the method is an optimized design method. Under the same technical specifications of the filter obtained by this method than the window function method, frequency sampling method to get the length of the filter are smaller, and to overcome the shortcomings of the window function design method and frequency sampling method, the maximum error minimize and evenly distributed throughout the approximation band.Key words：FIR digital filter；window function design method；frequency sampling method；MATLAB；Remez function目 录序 言1第1章FIR数字滤波器的基本原理及设计方法21.1窗函数设计法21.2频率采样设计法法21.3等波纹最优设计法2第2章 设计步骤4第3章 性能分析比较16第4章 MATLAB的实现194.1 MATLAB函数介绍24.2窗函数设计法24.3等波纹法设计举例2第5章 设计总结19参考文献20致 谢21III序 言数字滤波器的输入输出均为数字信号，信号通过数字滤波器后，可以改变频率成分的相对比例或滤除某些频率成分。数字滤波器可以分为IIR数字滤波器和FIR数字滤波器。文中只介绍FIR数字滤波器的设计，可以根据所给定的频率特性直接设计FIR数字滤波器。FIR数字滤波器在保证幅度特性满足要求的同时，能够做到严格的线性特性。本文分别介绍了窗函数法、频率采样法以及基于remez函数的最优化方法设计FIR滤波器，其中详细介绍了等波纹最优化设计法。在文中对FIR滤波器进行了详细的理论分析，并且对应于每种方法都给出了设计实例。通过编写MATLAB语言程序，运行程序，得到幅频和相频特性图。窗函数法和频率采样法都是容易理解和实现的。然而他们都存在某些缺陷。而等波纹最优化设计法却能够在最大程度上弥补这些缺陷，并且能够利用计算机和MATLAB工具方便的进行模拟和仿真，大大提高了工作效率，本文将对其原理及实现方法进行详细的介绍。- 0 -第一章 FIR数字滤波器的基本原理及设计方法有限长单位脉冲响应数字滤波器（FIRDF，Finite Impulse Response Digital Filter）的最大优点是可以实现线性相位滤波。而IIRDF主要对幅频特性进行逼近，相频特性会存在不同程度非线性。我们知道，无失真传输与滤波处理的条件是，在信号的有效频谱范围内系统幅频响应为常数，相频响应具有线性相位。在数字通信和图像处理与传输等应用场合都要求滤波器具有线性相位特性。另外FIRDF是全零点滤波器，硬件和软件实现结构简单，不用考虑稳定性问题。所以，FIRDF是一种很重要的滤波器，在数字信号处理领域得到广泛应用。当幅频特性指标相同时，FIRDF的阶数比IIRDF高的多，但是同时考虑幅频特性指标和线性相位要求时，IIRDF要附加复杂的相位校正网络，而且难以实现严格线性相位特性。所以，在要求线性相位滤波的应用场合，一般都用FIRDF。FIRDF的设计方法主要有两类：第一类是基于逼近理想滤波器特性的方法，包括窗函数法，频率采样法和等波纹最佳逼近法。第二类是最优设计法，我们主要讨论第一类设计法，侧重与滤波器的设计方法和相应的MATLAB工具箱函数的介绍。 FIR数字滤波器的设计方法有窗函数法、频率采样法和基于remez函数的最优化方法。MATLAB语言中的数字信号处理工具箱，提供了一些滤波器的函数，使FIR滤波器的运算更加方便和快捷。在MATLAB中提供的滤波函数有fir1,此函数以经典的方法实现加窗线性相位FIR数字滤波器设计，可以设计出低通、高通、带通和带阻滤波器；fir2函数设计的FIR滤波器，其滤波的频率特性由矢量f和m决定，f和m分别为滤波器的期望幅频响应的频率相量和幅值相量。firl和remez的基本格式用于设计I型和II型线性相位FIR滤波器，I型和II型的区别是偶函数还是奇函数。freqz()用于求数字滤波器的频率响应。并且提供了各种窗函数的函数，比如，hamming()是海明窗函数，hanning()是汉宁窗函数，kaiser()是凯泽窗函数，使在设计的过程中，不用自己重新设计窗函数。1.1 窗函数法设计FIR数字滤波器设我们所要设计的FIR滤波器的传输函数是(e),(n)是与其对应的单位脉冲响应，因此 （1-1） （1-2)如果我们能够在已知的情况下，求出，经过Z变换可得到滤波器的系统函数。通常情况下理想数字滤波器的单位脉冲相应是无限长的，且是非因果序列。获得有限脉冲响应滤波器的一种可能方法是对截取一段来近似代替，可是这样会改变原来的滤波器指标，出现吉布斯效应误差。窗函数法就是用被称为窗函数的有限加权序列w(n)来修正式(1)的傅里叶基数以求得要求的有限脉冲响应序列，即 （1-3）w(n)是有限长序列，当nN-1时，w(n)=0。这种方法的重点在于选择某种合适的窗函数。要求窗函数主瓣宽度尽可能窄，以获得最小的过渡带；旁瓣相对值尽可能小，以使得通带波纹小，并且阻带衰减大。下面介绍几种常用的窗函数：1.矩形窗(Rectangle Window) （1-4）2.三角形窗(Bartlett Window) （1-5）3.汉宁(Hanning)窗，又称升余弦窗 （1-6）4.汉明(Hamming)窗，又称改进的升余弦窗 （1-7）5.布莱克曼(Blankman)窗，又称二阶升余弦窗 （1-8)6.凯泽(Kaiser)窗 （1-9）其中：是一个可自由选择的参数，I0( x)是第一类修正零阶贝塞尔函数.贝塞尔函数的相关公式均可在教材中查找。上述窗函数的基本参数如下表窗函数旁瓣峰值幅度/db过渡带宽阻带最小衰减/db矩形窗-134/N-21三角形窗-268/N-25汉宁窗-318/N-44汉明窗-408/N-53布莱克曼窗-5712/N-74凯泽窗-5710/N-80窗函数法设计滤波器的步骤：1）根据技术要求确定待求滤波器的单位取样响应。2）根据对过渡带和阻带衰减的要求，选择窗函数的形式，并估计窗口长度N。3）计算滤波器的单位取样响应h(n)： （3-16)式中,是前面所选择好的窗函数。4）检验技术指标是否满足要求。根据下式计算： (3-17)如果不满足要求，根据具体情况重复步骤(2)(3)(4)步，直到满足要求为止。本文以一个FIR滤波器的设计为例说明如何使用MATLAB设计数字滤波器设计实例：用窗函数法设计线性相位FIR低通数字滤波器，要求通带截止频率Wp=0.4,阻带截止频率Ws=0.5, 通带衰减不大于3db，阻带衰减不小于40db。程序如下：Wp=0.4*pi;Ws=0.5*pi;Wdel=Ws-Wp;N=ceil(8*pi/Wdel);Wn=(0.4+0.5)*pi/2;window=hanning(N+1);b=fir1(N,Wn/pi,window);freqz(b,1,512)程序执行后得幅频和相频如下图所示： 图1.11.2 频率采样法设计FIR数字滤波器这种设计方法依据的是系统函数能够从频率响应的样本中求得这样一个事实。再者，这种设计方法非常适合频率采样型结构。对理想滤波器的系统函数Hd(z)进行频率采样以得到系统的理想频响Hd(ejw)的等间隔采样值H(k)。H(k)实际上是所要求的滤波器的单位采样响应(h(n)的离散傅里叶变换（DFT），如下式： (3-18) (3-19)为了减小H（k）的通带边缘由于抽样点的变化而引起的起伏振荡，可以增加过渡点，加宽过渡带以减小通带的起伏。每一个抽样值产生一个与sin()/sin()成正比，并位移（2k）/N的频率响应，而H（k）与内插函数的线性组合就是FIR滤波器的频率响应，增加一点过渡可以使阻带衰减提高到-4454dB,二点过渡衰减-6575,三点过渡衰减-8595dB.如果不能使过渡带太宽，同时要求增大阻带衰减，可以增加取样点数N，但这样会增加计算量、延时和误差。频率取样型FIR滤波器设计步骤：(1)给定理想滤波器频率响应。(2)根据过渡带宽和阻带衰减确定过渡点数和h(n)的长度N。 (3-20)(3)由IFFT计算IDFT得到： (3-21)设计实例：用频率采样法设计一个带通滤波器，满足：低阻带边缘：w1s=0.2*;低通带边缘：w1p=0.35*;高通带边缘：w2p=0.65*;高阻带边缘：w2s=0.8*。设计过渡带中的频率样本值为t1和t2，取t1=0.109021,t2=0.59417456。设计程序如下：M=40;al=(M-1)/2;l=0:M-1;t1=0.109021;t2=0.59417456;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);freqz(h,1,512,1000)实验得幅频相频特性如下图所示： 图1.21.3等波纹最优化方法设计FIR数字滤波器在数字信号处理中, 利用数字滤波器可改变信号中所含频率分量的相对比例或滤除某些频率分量, 使其达到所需要的效果. 其中数字FIR 滤波器由于具有精确的线性相位, 且系统稳定, 所以广泛应用于通信、数字图象处理、语音信号处理、自适应处理、雷达/ 声纳系统等方面. 目前, FIR 滤波器的设计方法主要有窗函数设计法和频率采样设计法. 但是, 这2种方法都不易精确控制通带边界频率p与阻带边界频率s ,所以, 在实际应用中具有一定的局限性.与窗函数和频率采样法比较，由于这种设计法使最大误差均匀化，所以设计的滤波器性能价格比最高。阶数相同时，这种设计法使滤波器的最大逼近误差最小，即通带最大衰减最小，阻带最小衰减最大；指标相同时，这种设计法使阶数最低。第二章 设计步骤等波纹滤波器的最优化设计方法主要有2 种,第 1 种是离散最小二乘法. 它的思路是在给定的一些离散点上,使实际的幅频特性和理想幅频特性之间的误差的平方和为最小. 第 2 种是最优化等波纹设计法,也称为雷米兹法或切比雪夫逼近法. 该类型滤波器幅频特性在通带和阻带上的误差峰值是均匀分布的,其误差具有等波纹特性,因而把波纹的幅度控制到最小, 或在同等指标下减小它的阶次. 第 1 种方法是连续最小的平方法的推广, 容易理解,但它的指标与滤波器没有直接关联,误差平方小的滤波器不能保证没有窄而大的波纹出现,像吉布斯效应 那样. 第 2 种方法直接控制通带波动和阻带衰减, 最具针对性,是滤波器的最优化设计方法. 因此,采用 MATLAB信号处理工具箱提供的函数,运用最优化等波纹设计法实现数字FIR滤波器的设计和仿真.。下面首先构造一种最大值最小的最优FIR滤波器设计问题，并讨论在线性相位FIR滤波器的振幅响应中能够得到的最大值和最小值的总数。利用这个然后讨论一般等波纹FIR滤波器设计算法，这种算法对其解利用多项式内插，该算法即是Parks-McClellan算法。最大值最小问题的确立 对于线性相位FIR数字滤波器的4种情况，其加权误差表示为E()=W()Hd()-P()，在S内。可以证明，切比雪夫近似问题能定义为：确定一组系数a(n),使在通带和阻带内E（）的最大绝对值最小，即minmax丨E()丨现在已经成功的给定了准确的p, s, 1, 2;另外，误差在通带和阻带内都能均匀分布了。极值数目的限制对于一个给定的M点滤波器来说，在误差函数E()内存在多少局部最大值和最小值？这一信息被Parks-McClellan算法用来得到多项式内插。由P()的多项式特性，我们得到E()在S内至多有L+3个极值。另，由交错点定理，最优等波纹滤波器在S内的误差函数不是有(L+2)个就是(L+3)个交错点。其中，大多数等波纹滤波器有(L+2)个交错点。然而对于某些p和s的组合，能够得到(L+3)个交错点。Parks-McClellan算法交错点定理确保最大值最小近似问题的解存在且是唯一的，但是它并没有告诉我们如何求这个解。我们既不知道长度M，也不知道极值频率i，也不知道参数(n)，也不知道最大误差。该算法利用Remez交换算法提供了一种迭代解。它假定滤波器长度M和比值2/1已知，如已选取加权函数，并正确的选取了阶数M，那么得到解时=2。显然，和M是相互关联的；M越大，越小。在滤波器设计数据中，p，s，1，2是已知的，因此，M需先假定知道。所幸，由凯斯公式Parks 和McClellan 算法首先猜设( L + 2)个极值频率i，估计极值频率上的最大误差及其i，由这些新的频率点拟合出一个新的L 阶多项式，重复以上过程，直到找出最优集i和全局最大波纹系数为止. 迭代结果得到多项式系数(n) ,最后算出滤波器脉冲响应 h(n) .由于M是近似的, 可能不等于2 . 若2 ,则增加M，若2 ,则减小M,然后再次用Remez 算法确定一个新的.重复此过程, 直至2,这样就可得到等波纹滤波器.第三章性能分析比较分别采用窗口法、频率取样法和等波纹最优化设计法设计一个数字FIR低通滤波器，其指标要求为： p=0.2，Rp=0.25dB；s=0.3，As=50dB 选取Hamming窗函数，它给出比较小的过渡带，因此阶数较低。MATLAB脚本省略。其幅频响应如下： 图3.1频率采样法所得幅频响应图像如下： 图3.2 采用等波纹最优化设计法设计该低通滤波器，其幅频响应如下： 图1.5比较以上三种滤波器的设计方法，在同样的阶数下，优化设计可以获得最佳的频率特性和衰耗特性，具有通带和阻带平坦，过渡带窄等优点。频率取样法设计过程简单，但阻带衰减明显，若适当选取过渡带样点值，会取得较窗口法略好的衰耗特性。窗口设计法在阶数较低时，阻带特性不满足设计要求，只有当滤波器阶数较高时，使用哈明窗和凯塞窗基本可以达到阻带衰耗要求。FIR 滤波器以其稳定和容易实现严格的线性相位，使信号处理后不产生相位畸变，而在实际中获得广泛的应用。第四章MATLAB的实现4.1MATLAB函数介绍1 Remez函数采用remez算法实现线性相位FIR数字滤波器的等波纹最优化设计。其调用格式为：hn=remez(M，f，m，w)返回单位脉冲响应响亮hn。调用参数含义如下：M：FIRDF阶数，hn长度为N=M+1。f和m：给出希望逼近的幅度特性。f为边界频率向量，0f1，要求为f为单调增向量（即f(k)f（k+1），而且从0开始，以1结束，1对应数字频率=（或模拟频率Fs2，Fs表示时域采样频率）。m是与f对应的幅度向量，m与f长度相等，m（k）表示频点f(k)的频响应值。如果用命令Plot（f,m）画出幅频响应曲线，则k为奇数时，频段f(k),f(k+1)上的幅频响应就是期望逼近的幅频响应值，频段f（k+1），f（k+2）为无关区。简言之，Plot（f，m）命令画出的幅频响应就是期望逼近的幅频响应值，频段f（k+1），f（k+2）为无关区。简言之，Plot(f,m)命令画出的幅频响应曲线中，起始频段为第一段，奇数频段为逼近区，偶数频段为无关区。除了设计选频FIRDF，remez函数还可以设计两种特殊滤波器，调用格式如下。设计希尔伯特变换器的调用格式为：hn=remez（M，f，m，w，hilbert）hn具有奇对称特性：hn(n)= -hn（M+1-n）。在边界频率向量f给定的通带上希望幅度逼近常数1.一般f=a,b，m=1,1,且f中第一个边界频率a不能为0，必须满足0ab1.这是因为希尔伯特变换器属于第二类线性相位滤波器，且滤波器长度N为奇数，所以只能实现带通滤波器。设计数字微分器调用格式为：hn=remez（M,f,m,w,defferentiator）Remez函数的调用参数（M，f，m，w），一般都是调用remezord函数来计算。2.remezord功能：根据逼近指标估算等波纹最佳逼近FIRDF的最低阶数m、误差加权向量M和归一化边界频率向量f，从而是滤波器在满足指标的前提下造价最低。其返回参数作为remez函数的调用参数。其调用格式为： M，fo，mo，w=remezord（f，m，rip，Fs） 参数说明如下。 返回参数作为remez函数的调用参数，设计的滤波器可以满足由参数f，m，rip和Fs满足的指标。f与remez中类似，这里f可以是模拟频率（Hz）或归一化数字频率，但必须以0开始，以Fs2（用归一化频率时对应1）结束，而且其中省略了0和Fs2两个频点。Fs为采样频率，缺省时默认Fs=2Hz。但是这里f的长度（包括省略的0和Fs2两个频点）是m的两倍，即m中的每个元素表示f给定的一个逼近频段上希望逼近的幅度。 注意：省略Fs时，f中必须为归一化频率。有时估算的阶数M略小，使设计结果达不到指标要求，这时要取M+1或M+2（必须注意对滤波器长度N=M+1的奇偶性要求）。所以必须检验设计结果。如果无关区（过渡带）太窄，或截止频率接近零频率和FS2时，设计结果不正确。 rip表示f和m描述的各逼近频段允许的波纹振幅（幅频响应最大偏差），f的长度是rip的两倍。 一般以N,fo,mo,w=remezord（f,m,rip,Fs）返回的参数作为remez的调用参数，计算单位脉冲响应：hn=remez（N，fo,mo,w），对比前面介绍的remez调用参数，可清楚地看出remezord返回参数N，fo,mo和w的含义。 综上所述，调用remez和remezord函数设计线性相位FIRDF，关键是根据设计指标求出remezord函数的调用参数f,m,rip和Fs，其中Fs一般是题目给定的，或根据实际信号处理要求（按照采样定理）确定。 下面给出由给定的各种滤波器设计指标确定remezord调用参数f,m,和rip的公式，编程时直接套用即可。 低通滤波器设计指标逼近通带：0，p，通带最大衰减：pdB；逼近阻带：s，阻带最小衰减：sdB。remezord调用参数：f=p，s，m= 1，0，rip= 1，2其中，f向量省去了起点频率0和终点频率1，1和2分别为通带和阻带波纹幅度，由通带波纹及阻带衰减的相应公式均可换算得到，下面相同。 高通滤波器设计指标逼近通带：p，通带最大衰减：pdB；逼近阻带：0，s，阻带最小衰减：sdB。remezord调用参数： f= s，p，m= 0，1，rip= 2，1 带通滤波器设计指标逼近通带：p1，pu，通带最大衰减：pdB；逼近阻带：0，s1，su，阻带最小衰减：sdB。remezord调用参数： f= s1，p1，pu，us，m= 0，1，0，rip= 2，1，2。 带阻滤波器设计指标逼近阻带：s1，su，阻带最大衰减：sdB；逼近阻带：0，p1，pu,，通带最小衰减：pdB。remezord调用参数： f= p1，s1，su，pu，m= 1，0，1，rip= 1，2，1。4.2设计举例利用等波纹最优化设计法设计FIR带通滤波器，其指标如下：阻带下截止频率s=0.2，通带下截止频率p=0.35，通带上截止频率p=0.65, 阻带上截止频率s=0.8,通带最大衰减p=1dB,阻带最小衰减s=60dB。解：调用remezord和remez函数求解。由调用格式知道，首先要根据设计指标确定remezord函数的调用参数，在直接编写程序调用remezord和remez函数设计得到h(n)。将设计指标带入响应公式得到remezord函数的调用参数。其MATLAB程序较为繁琐，省略不写，时域及频域图如下：图4.1程序运行结果：M=28.即h(n)的长度N=29。h(n)及其耗损函数曲线如图7.4.5所示。而用窗函数设计的FIR滤波器的阶数为80，由此例可见，等波纹最佳逼近法可以使滤波器的阶数大大降低。FIR希尔伯特变换器在频带0.050.95内设计一个希尔伯特变换器。解：因为这是一个宽带希尔伯特变换器，要为该滤波器选取一个奇数长度，现选取M=51。其MATLAB程序如下：f=0.05,0.95;m=1,1;h=remez(50,f,mhilbert);db,mag,pha,grd,w=freqz_m(h,1);subplot(2,1,1);stem(0:50,h);title(Impulse Response);xlabel(n);ylabel(h(n);axis(0,50,-0.8,0.8);set(gca,XTickMode,manual,XTick,0,50)set(gca,YTickMode,manual,YTick,-0.8:0.2:0.8);subplot(2,1,2);plot(w/pi,mag);title(Magnitude Respo