附录D 最优等波纹线性相位FIR滤波器的设计

附录D最优等波纹线性相位FIR滤波器的设计对于线性相位FIR滤波器的设计方法，窗函数和频率采样法是相对简单的方法，然而，它们都有存在不能精确地控制和这类关键频率的问题。本节描述的滤波器设计方法采用切比雪夫等波纹逼近思想，为了将理想幅度特性和实际幅度特性之间的加权逼近误差均匀地分散到滤波器的整个通带和阻带，并且最小化最大误差，则采用切比雪夫逼近方法被视为最优设计准则。所得到的滤波器结构在通带和阻带都有等波纹。下面以低通滤波器的设计为例来说明设计过程，考虑通带截止频率为和阻带频率为的低通滤波器的设计。如图D-1所示，图中给出了一般技术指标，在通带内滤波器幅度特性应满足的条件为图D-1低通滤波器的最佳逼近， (D-1)类似地，在阻带内规定滤波器幅度特性落在范围之间，即， (D-2)式中，表示通带波纹的峰值，表示阻带波纹的峰值。 现在集中考虑四种产生线性相位FIR滤波器的情况，这些在前面已经讨论过，总结如下。（1）情况1：当，且奇数时 ， (D-3)式中 (D-4)（2）情况2：当，且偶数时 ， (D-5)式中， (D-6)进一步对式(D-5)进行整理和重新排列，得到， (D-7)其中，系数和系数线性相关，可以证明两者之间存在如下关系，， (D-8)（3）情况3：当，且奇数时 ， (D-9)式中 (D-10)进一步对式(D-9)进行整理和重新排列，得到， (D-11)其中，系数和系数线性相关，从式(7-2-9)和式(7-2-11)可以推导出两者之间存在如下关系， (D-12)（4）情况4：当，偶数时 ， (D-13)式中 (D-14)与前面情况一样，可以对式(D-13)进行整理和重新排列，得到， (D-15)其中，系数和系数线性相关，可以证明两者之间存在如下关系， (D-16)归纳这四种情况中的表达式，并列于表D-1中。注意到，在情形2、情形3和情形4中所做的整理，允许我们将统一表示为 (D-17)其中 (D-18)且有普遍形式 (D-19)其中，表示滤波器的系数，它与FIR滤波器的单位脉冲响应线性相关。对于情况1来说，求和的上限；对于情况3，求和的上限；对于情况2和情况4，求和上限。 除了在上面给出的关于表达的普通框架外，这里还定义实值理想幅度特性函数和误差加权函数。其中实值理想幅度特性可简单定义为通带内等于1，而阻带内等于0。误差加权函数允许我们在不同的频带内（即通带和阻带）选取相对大小的误差，具体地说，既可以通过在阻带内将进行归一化，也可以在通带内令，即 (D-20)然后，就可以仅仅在通带内选择合适的加权系数，可以约束阻带中的波纹与通带中的波纹之间相对大小。表D-1线性相位FIR滤波器的四种情况根据和的技术指标，可以定义加权逼近误差为 (D-21)为了便于描述，我们定义一个修正的误差加权函数和一个修正的理想幅度特性，分别为 (D-22) (D-23)则对所有四种不同类型的线性相位FIR滤波器来说，加权逼近误差可以表示为 (D-24) 切比雪夫逼近问题基本上就是由加权误差函数确定滤波器参数，使得在所逼近的频带上的最大绝对值可以最小化。实际上就是求如下方程的解 (D-25)式中，表示在其中实施最优化的频带集合（不相交的并集）。基本上，集合由理想滤波器的通带和阻带组成。 Parks和McClellan在1972年解答了这一问题，他们在切比雪夫逼近理论中应用了一个称之为交错定理的定理。下面，在不加以证明的情况下对该定理进行阐述。 交错定理：令为区间的一个闭子集，为了使 是在上为的唯一最佳加权切比雪夫逼近，其充分必要条件是误差函数在上至少呈现个极值曲线频率点，即在上必须至少存在个频率，使得当<<<时，，且， (D-26) 由此可以看到，两个相邻的极值频率点之间，加权误差函数在符号上交替变化，因此，该定理称为交错定理。 下面以低通滤波器为例，对交错定理进行详细说明。设低通滤波器通带为,阻带为的低通滤波器。因为理想幅度特性和误差加权函数是逐段恒定的，所以有因此，加权误差函数各峰值点的频率，对应于满足误差容限时的峰值点。因为是次三角多项式，所以，对于情况1来说，利用三角恒等式则可表示为 (D-27)其中系数与系数和有关。由此得出，在开区间上，至多能有个局部极大值和极小值。另外，和通常是的极值点，并且它们同时也是的极值点。所以，至多能有个极值频率点。此外，和也是的极值点，应为在和点处，具有最大值。故对于理想低通滤波器的唯一最佳近似，在上至多存在个极值频率点。另一方面，交错定理说明在上至少存在个极值频率点，所以低通滤波器加权误差函数不是有就是有个极值频率点。一般来说，如果滤波器结构包含的交错或波纹数目多于，则这种滤波器称为“超波纹”滤波器。当滤波器结构包含有最大数目的交错时，这类滤波器称为“最大波纹滤波器”。 交错定理保证了在式(D-25)中的切比雪夫最优化问题有唯一解，在理想极值频率处，存在方程组， (D-28)式中，表示加权误差函数的最大值。事实上，如果按照式(D-20)所指明的那样选取，则由此得出。将式(D-28)进行整理重排，可以得到，或者，等效为， (D-29)如果把式(D-29)中参数和当做被确定的参数来处理，那么式(D-29)可以按照矩阵形式表示为 (D-30)不过，上述提供的方法是在这些交错点频率给定下得到的。实际上是不知道的，所以直接求解比较困难，只能用逐次迭代的方法来解决，可以利用Remez交换算法的迭代算法来求解。在迭代算法中，我们首先从猜测一组（）个频率开始，接着确定和，然后计算加权误差函数，又可确定另一组（）个极值频率点，重复迭代过程，直到所得到的频率组收敛于最佳的极值频率组为止。虽然，在式(D-30)中的矩阵方程可以在迭代方法中使用，但是矩阵求逆的运算量非常大，所以效率较低。一种更为有效的方法是由Rabiner等（1975）在论文中提出。该方法按照公式 (D-31)解析计算。其中， (D-32)式(D-31)是直接从式(D-30)中的矩阵方程得出的，于是，利用个极值频率点的初始猜测值，就可以计算。现在，由于是形如 的三角多项式，并且已知在点，处，多项式相应的值为， (D-33)所以，利用关于的拉格朗日插值公式，可表示为 (D-34)式中，可由式(D-33)计算出，而为 (D-35)求出后，就可以在频率点密集上计算加权误差函数 (D-36)通常，频率点密集中点的个数等于16就足够了，其中是滤波器的长度。如果对于密集的某些频率来说，，则选取一组新的对应于的前个大峰值频率点，并重复式(D-31)的计算过程。因为选取的这组新的个峰值频率点对应于的峰值，所以，算法在每次迭代中都使增加，直到收敛于上界为止。因此也收敛于切比雪夫逼近问题的最佳解。换言之，当密集中的所有频率点都满足时，通过多项式就已找到最优解。该算法的流程由Remez（1957）提出。一旦通过得到最优解，在不必计算参数的情况下，就可以直接计算单位脉冲响应，事实上，我们已经可以确定具体来说，对于为奇数的情况，在，，计算。然后，依据所设计的滤波器类型，从对应的奇、偶对称情况中选取相应的公式来确定。对基于切比雪夫逼近准则和用Remez交换算法实现的线性相位FIR滤波器来说，由Parks和McClellan（1972）编写的计算机程序是现成的。该程序可以用来设计低通、高通或带通滤波器。利用MATLAB实现最佳逼近设计。交替定理能确保最大-最小逼近问题的解存在且唯一，但并未说明如何去求这个解。因为既不知道滤波器长度，也不知道极值频率，不知道参数，也不知道最大误差。Parks_McClellan利用Remez交换算法给出一种迭代解。该解假定滤波器长度和比值已知，并按式(D-20)选取误差加权函数。如果阶数选取合适，那么得到解时就有。显然，和是有关联的。越大，就越小。在滤波器设计中是已知的，所以必须假定先知道。Kaiser提出了一个简单公式来逼近。， (D-37)最优等波动滤波器的误差在频带上极值数目有个或个，大部分情况有个。Parks_McClellan算法首先设定有个极值频率点，估计这些频率上的最大误差，按照给定的各点，拟合一个阶多项式(D-34)形式。然后在频率点密集上确定个极值频率，这些新频率点又可拟合出一个新的阶多项式。重复以上过程直到找到最佳和全局最大值为止。最后算出滤波器单位脉冲响应。由于是近似的，最大误差可能不等于。如果是这样，需要增加()或减小()，再次用remez算法确定一个新的。重复此过程直至。这样就得到了等波动滤波器。MATLAB提供了一个函数firpm实现Parks_McClellan算法，调用格式通常有[h]=firpm(N,f,A):设计一个阶FIR滤波器，给定频率点数组{}和希望幅度特性。的单位为，即其取值，通带和阻带误差加权函数为1。返回单位脉冲响应。[h]=firpm(M,f,A,W):与上述情况类似。W为误差加权函数，长度为的一半。其他调用可参考helpfirpm命令。求出后必须检查最小阻带衰减，并与给定的比较，然后增加(或减小)滤波器阶数，重复此过程直到得到期望的。【例D-1】设计一个最优化等波动数字FIR低通滤波器，指标为，，，。解MATLAB参考程序如下。运行结果见图D-2。wp=0.3*pi;ws=0.4*pi;Rp=0.4;As=40;%给定指标delta1=(10^(Rp/20)-1)/(10^(Rp/20)+1);delta2=(1+delta1)*(10^(-As/20));%求波动指标weights=[delta2/delta11];deltaf=(ws-wp)/(2*pi);%给定权函数和△f=wp-wsN=ceil((-20*log10(sqrt(delta1*delta2))-13)/(14.6*deltaf)+1);N=N+mod(N,2); %估算阶数NN=N+4f=[0wp/piws/pi1];A=[1100];%给定频率点和希望幅度值h=firpm(N,f,A,weights);%求单位脉冲响应[db,mag,pha,grd,W]=freqz_m(h,[1]);%验证求取幅度特性函数delta_w=