一种降低通带群延迟误差的FIR滤波器设计方法

摘要

本发明公开了一种降低通带群延迟误差的FIR滤波器设计方法，包括以下步骤：1)基于最小二乘准则建立待设计滤波器与理想频率响应的通带群延时误差函数，将通带群延时误差函数作为设计模型的目标函数，同时设定约束条件，即滤波器频率响应误差小于某一设定值；2)通过多次迭代最小化滤波器通带最大群延时误差，使FIR滤波器的群延时在通带具有等波纹特性；3)在每次迭代过程中，采用泰勒展开将其简化为非凸的二次优化问题，然后再将其转化为一个等价的凸优化问题，最后利用内点法求解问题，所得结果带入下次迭代，直至迭代终止条件得到满足。本发明将滤波器通带最大群延时误差作为设计问题的优化目标，从而有效提升了滤波器的群延时性能。

说明书

技术领域

本发明涉及一种FIR滤波器的设计方法，特别涉及一种减少通带群延迟误差 的滤波器设计。

背景技术

有限冲激响应（Finite Impulse Response，简称FIR）滤波器是数字信号处理 中常用的滤波器之一，被广泛用于通信、语音/图像/视频、雷达/声纳等信号处理 场合。与无限冲激响应（Infinite Impulse Response，简称IIR）相比，FIR滤波器 具有以下一些优点：

（1）利用对称结构，在保证满足滤波器幅频响应要求的同时，设计者还可获得 具有严格线性相位的FIR滤波器，从而保证通过滤波器的信号不会发生相 位失真；

（2）由于FIR滤波器的极点均位于z平面原点，从而不存在IIR滤波器所面临 的稳定性问题；

（3）由于结构相对简单，许多FIR滤波器设计问题可以表述为一个凸优化问题， 从而设计者可以获得最佳设计结果。

如果待滤波信号的相位并不携带过多的信息，设计过程可以只考虑滤波器的 幅度响应，此时的问题就得到了极大的简化，现有的一些设计方法可以比较可 靠地找到最佳的设计结果。但是，对于另外的一些应用场合，我们希望信号不 同的频率成分经过近似相等的时延到达系统输出端，此时则需要考虑如何在设 计过程中保证群延时满足设计要求。传统的FIR滤波器设计方法通常假设FIR 滤波器具有对称结构，因此设计结果只需要满足幅频要求即可。但是，这类线 性相位滤波器的群延时必须等于滤波器阶数的一半。因此，当滤波器阶数过大 时，相应的FIR滤波器的群延时也会变得很大，而这一特性并不能满足一些实 时性要求较高的应用场合。基于这一现象，另一类FIR滤波器设计方法则放弃 了对称结构。在设计之初，设计者可以给定待设计滤波器的理想频域响应，滤 波器的通带群延时作为理想频率响应的一部分带入了设计过程，从而保证在满 足幅度响应的同时，FIR滤波器的群延时能够得到一定的降低。但是，由于滤波 器的群延时具有非常复杂的数学形式，为了降低设计难度，这类设计方法将频 率响应误差作为整体设计误差。因此，在设计过程中，我们很难单独控制滤波 器的群延时误差。此外，为了简化设计问题，传统的FIR滤波器设计方法通常 假设滤波器的群延时恒为常数，这也不便于我们设计一些具有非线性群延时的 滤波器（例如，数字微分器）。

发明内容

本发明所要解决的问题是针对传统FIR数字滤波器设计方法难于有效控制 通带群延时误差这一缺点，在不影响总体性能的基础上，提出一种有效降低FIR 滤波器群延时误差的方法，使得所设计FIR滤波器的群延时在通带具有近似的 等波纹特性。

为解决上述技术问题，本发明提供一种降低通带群延迟误差的FIR滤波器 设计方法，其特征在于，包括以下步骤：

1)基于最小二乘准则建立待设计滤波器与理想频率响应的通带群延时误差函 数，将通带群延时误差函数作为设计模型的目标函数，同时设定约束条件，即 滤波器频率响应误差小于某一设定值；

2)通过多次迭代最小化滤波器通带最大群延时误差，使FIR滤波器的群延时 在通带具有近似的等波纹特性；如果是严格的等波纹设计的话，就表示所用设 计方法可以得到模型式(1)的最优设计，但是由于式(1)较为复杂，无法得到最优 设计，所以方法中用了近似的处理手段，获得较为简单的模型（例如(2)和(3)）, 但是此时无法保证由简化模型所得设计为最优设计，所以称之为“近似”等波 纹；

3)在每次迭代过程中，采用泰勒展开将其简化为非凸的二次优化问题，然后 再将其转化为一个等价的凸优化问题，最后利用内点法求解问题，所得结果带 入下次迭代，直至迭代终止条件得到满足。

前述的降低通带群延迟误差的FIR滤波器设计方法，其特征在于：在所述步 骤1)中,所述设计模型为：

$\underset{h}{\min }{E}_{L\infty }^{\mathrm{GD}}\left(h\right)=\underset{\omega \in B_P}{\max }\left|e(h,\omega )\right|---\left(1\right)$

s.t.E_L2(h)＝h^TQh-2q^Th+c≤δ

以上模型所表示的设计任务为：在保证滤波器最小二乘频率响应误差E_L2(h) 小于上限值δ的前提下，通过调节待设计的N阶FIR滤波器系数 $h={\left(\begin{array}{ccccc}{h}_0& h_1& ...& h_{N-1}& h_N\end{array}\right)}^T,$最小化通带内最大群延时误差其中，h_n为第n个滤 波器抽头系数，B_p为通带集合，e(h,ω)为FIR滤波器在频率ω处的群延时误差，

$Q={\int }_{{\Omega }_I}\Psi \left(\omega \right)\mathrm{d\omega }$

$c={\int }_{{\Omega }_I}{\left|H_d\left(\omega \right)\right|}^2\mathrm{d\omega }$

上式中，上标T、H和*分别表示矩阵或向量的转置操作、复数矩阵或向量的 共轭转置操作、及复数的取共轭操作，Ω_I表示设计者感兴趣的频带组合（例如， 通带和阻带），Re{·}表示复数的实部，j表示虚数单位，GD是群延时缩写，c 为常数，H_d(ω)为滤波器理想频域响应。

前述的降低通带群延迟误差的FIR滤波器设计方法，其特征在于：在所述步 骤3）中，具体包括以下步骤：

31)在每次迭代中，先建立设计模型的以下近似设计模型：

$\underset{h}{\min }{E}_{L2}^{\mathrm{GD}}(h,h^{(k-1)})={\int }_{B_p}{W}^{\left(k\right)}\left(\omega \right){\left|e(h,\omega )\right|}^2\mathrm{d\omega }---\left(2\right)$

s.t.h^TQh-2q^Th+c≤δ

其中，群延时误差函数的加权系数，W^(k)(ω)由前一次迭代中所用加权系数 W^(k-1)(ω)及前一次迭代所得设计结果的群延时误差函数的包络|e(h^(k-1),ω)|^p乘积所 得，上标(k)表示迭代次数，h^(k-1)为上一次迭代所得的设计结果；

32)在获得设计模型式（2）后，对|e(h,ω)|²在h^(k-1)处进行二阶泰勒展开，从 而进一步获得以下的简化设计模型

$\underset{h}{\min }{h}^T{G}^{\left(k\right)}h+2g^{\left(k\right)T}h---\left(3\right)$

s.t.h^TQh-2q^Th+c≤δ

其中，

$G^{\left(k\right)}=\frac{1}{2}{\int }_{B_p}{W}^{\left(k\right)}\left(\omega \right)G(h^{(k-1)},\omega )\mathrm{d\omega }---\left(4\right)$

$g^{\left(k\right)}=\frac{1}{2}{\int }_{B_p}{W}^{\left(k\right)}\left(\omega \right)g(h^{(k-1)},\omega )\mathrm{d\omega }-G^{\left(k\right)}{h}^{(k-1)}---\left(5\right)$

式（4）和式（5）中，g(h,ω)和G(h,ω)分别为|e(h,ω)|²在h处的梯度和Hessian 矩阵；

33)将简化设计模型式（3）等价转化为以下的凸优化问题：

$\underset{z}{\min }\underset{i=0}{\overset{N}{\Sigma }}({\lambda }_i{z}_i-2|v_i^{\left(k\right)}\left|\sqrt{z_i}\right)$

$s.t.\underset{i=0}{\overset{N}{\Sigma }}{z}_i\le \delta -c+q^T{Q}^{-1}q---\left(6\right)$

z_i≥0,i＝0,1,...N.

其中，λ_i为矩阵的特征值，$v^{\left(k\right)}={\left(\begin{array}{ccc}{v}_0^{\left(k\right)}& ...& v_N^{\left(k\right)}\end{array}\right)}^T$可由以下公式计算所得

$v^{\left(k\right)}=U^{\left(k\right)T}{\tilde{g}}^{\left(k\right)}$

${\tilde{g}}^{\left(k\right)}=Q^{-\frac{1}{2}}{G}^{\left(k\right)}{Q}^{-1}q+Q^{-\frac{1}{2}}{g}^{\left(k\right)}$

上式中，U^(k)由矩阵的特征值向量构成；$\left(\begin{array}{ccc}{v}_0^{\left(k\right)}& ...& v_N^{\left(k\right)}\end{array}\right)$是向量v^(k)的组成元 素，z_i是(6)中待求解变量；

34)将模型式（6）为凸优化问题后，利用高效的内点算法进行求解。

前述的降低通带群延迟误差的FIR滤波器设计方法，其特征在于：式（3）的 最佳设计表示为

$h^{\left(k\right)}=Q^{-\frac{1}{2}}{G}^{\left(k\right)}{y}^{\left(k\right)}+Q^{-1}q---\left(7\right)$

其中，y^(k)中每个元素由下式计算所得

$y_i^{\left(k\right)}=-\mathrm{sign}\left(v_i^{\left(k\right)}\right)\sqrt{z_i^{\left(k\right)}}---\left(8\right)$

上式中，为模型（6）的最优解；

将式（7）中所得最佳设计h^(k)带入下一次迭代计算，直至迭代终止条件得到 满足。

本发明所达到的有益效果：本发明将滤波器通带最大群延时误差作为设计问 题的优化目标，从而有效提升了滤波器的群延时性能。设定滤波器最小二乘标 准下的频率响应误差小于某一设定值。在满足这一约束条件下，设计方法通过 多次迭代最小化滤波器通带最大群延时误差，从而保证FIR滤波器的群延时在 通带具有近似的等波纹特性；在迭代过程中，利用加权最小二乘准则重构设计 问题的目标函数，通过调整目标函数中不同频率点上群延时误差的权重，即在 每次迭代中，W^(k)(ω)由前一次迭代中所用加权系数W^(k-1)(ω)及前一次迭代所得 设计结果的群延时误差函数的包络|e(h^(k-1),ω)|^p计算所得，从而调整加权系数。减 小通带内最大群延时误差；重构后的设计问题目标函数仍然具有复杂的数学形 式，此时采用泰勒展开将其简化为非凸的二次优化问题，然后再将其转化为一 个等价的凸优化问题，利用高效的凸优化求解器最终获得问题的最优解。

附图说明

图1为本发明设计方法的基本流程图；

图2为本发明应用示例1中的低通滤波器幅度响应；

图3为本发明应用示例1中的低通滤波器的通带群延时；

图4为本发明应用示例2中的带通滤波器幅度响应；

图5为本发明应用示例2中的带通滤波器的通带群延时。

具体实施方式

下面结合附图，对本发明的设计方法做进一步阐述。

一、模型建立

最小二乘准则下，FIR滤波器的误差函数表示为

$E_{L2}\left(h\right)={\left|\right|H\left(e^{\mathrm{j\omega }}\right)-H_d\left(\omega \right)\left|\right|}_2$

(1-1)

$={\int }_{{\Omega }_I}{|H\left(e^{\mathrm{j\omega }}\right)-H_d\left(\omega \right)|}^2\mathrm{d\omega }$

其中，H(e^jω)为所设计的滤波器的频率响应，H_d(ω)为滤波器的理想频率响应，Ω_I表示所需考虑的频率段集合。式(1-1)可以改写为一个关于滤波器系数的二次函 数：

E_L2(h)＝h^TQh-2q^Th+c   (1-2)

其中，

$h={\left(\begin{array}{ccccc}{h}_0& h_1& ...& h_{N-1}& h_N\end{array}\right)}^T$

$Q={\int }_{{\Omega }_I}\Psi \left(\omega \right)\mathrm{d\omega }$

$c={\int }_{{\Omega }_I}{\left|H_d\left(\omega \right)\right|}^2\mathrm{d\omega }$

上式中，上标H和*分别表示矩阵或向量的共轭转置操作与复数的取共轭操作。

对于特定频率ω，滤波器的群延时τ(ω)可以由下式计算所得：

$\tau \left(\omega \right)=\frac{h^T\tilde{\Psi }\left(\omega \right)h}{h^T\Psi \left(\omega \right)h}---(1-3)$

其中，矩阵定义为：

$\tilde{\Psi }\left(\omega \right)=\frac{\mathrm{\Lambda \Psi }\left(\omega \right)+\Psi \left(\omega \right)\Lambda }{2}---(1-4)$

上式中，Λ表示一个对角元素为0、1、…、N的对角矩阵。

利用（1-2）和（1-3），设计问题可以表述为

$\underset{h}{\min }{E}_{L\infty }^{\mathrm{GD}}\left(h\right)=\underset{\omega \in B_P}{\max }\left|e(h,\omega )\right|---\left(1\right)$

s.t.E_L2(h)≤δ

其中，B_p表示通带集合，δ为预先设定的最大误差值，e(h,ω)定义为滤波器在频 率ω处的群延时误差

$e(h,\omega )=\tau \left(\omega \right)-{\tau }_d\left(\omega \right)=\frac{h^T[\tilde{\Psi }\left(\omega \right)-{\tau }_d\left(\omega \right)\Psi \left(\omega \right)]h}{h^T\Psi \left(\omega \right)h}---(1-5)$

上式中，τ_d(ω)为滤波器的理想群延时。

二、简化模型

实际中，很难对上述设计问题进行求解。因此，需要采用迭代算法进行求解。 在当前迭代，首先采用加权最小二乘准则改写模型式（1）中的目标函数

$E_{L2}^{\mathrm{GD}}(h,h^{(k-1)})={\int }_{B_p}{W}^{\left(k\right)}\left(\omega \right){\left|e(h,\omega )\right|}^2\mathrm{d\omega }---(2-1)$

其中，h^(k-1)表示上一次迭代结果，加权系数W^(k)(ω)可由W^(k-1)(ω)和|e(h^(k-1),ω)|^p的 包络乘积所得，p为一常数，用于加快迭代收敛速度。

为了进一步简化模型，利用泰勒展开，将（2-1）中的|e(h,ω)|²近似表示为一 个二次函数

${\left|e(h,\omega )\right|}^2\approx {\left|e(h^{(k-1)},\omega )\right|}^2+g^T(h^{(k-1)},\omega )(h-h^{(k-1)})$

(2-2)

$+\frac{1}{2}{(h-h^{(k-1)})}^{T}G(h^{(k-1)},\omega )(h-h^{(k-1)})$

其中，g(h,ω)和G(h,ω)分别表示|e(h,ω)|²在h处的梯度与Hessian矩阵。利用式 （2-2），此时的设计问题可以表述为

$\underset{h}{\min }{h}^T{G}^{\left(k\right)}h+2g^{\left(k\right)T}h---\left(3\right)$

s.t.h^TQh-2q^Th+c≤δ

其中，

$G^{\left(k\right)}=\frac{1}{2}{\int }_{B_p}{W}^{\left(k\right)}\left(\omega \right)G(h^{(k-1)},\omega )\mathrm{d\omega }---\left(4\right)$

$g^{\left(k\right)}=\frac{1}{2}{\int }_{B_p}{W}^{\left(k\right)}\left(\omega \right)g(h^{(k-1)},\omega )\mathrm{d\omega }-G^{\left(k\right)}{h}^{(k-1)}---\left(5\right)$

三、问题求解

需要指出的是，由于G^(k)为不定矩阵，因此模型式（3）并非一个凸优化问题。 对于非凸问题，通常存在多个局部最优的设计结果。但是，通过矩阵分解与变 量替代，模型式（3）可等价转化为以下凸优化问题

$\underset{z}{\min }\underset{i=0}{\overset{N}{\Sigma }}({\lambda }_i{z}_i-2|v_i^{\left(k\right)}\left|\sqrt{z_i}\right)$

$s.t.\underset{i=0}{\overset{N}{\Sigma }}{z}_i\le \delta -c+q^T{Q}^{-1}q---\left(6\right)$

z_i≥0,i＝0,1,...N.

其中，λ_i为矩阵的特征值，$v^{\left(k\right)}={\left(\begin{array}{ccc}{v}_0^{\left(k\right)}& ...& v_N^{\left(k\right)}\end{array}\right)}^T$可由以下公式计算所得

$v^{\left(k\right)}=U^{\left(k\right)T}{\tilde{g}}^{\left(k\right)}$

${\tilde{g}}^{\left(k\right)}=Q^{-\frac{1}{2}}{G}^{\left(k\right)}{Q}^{-1}q+Q^{-\frac{1}{2}}{g}^{\left(k\right)}$

上式中，U^(k)由矩阵的特征值向量构成。此时，式（6）为凸优化问题， 可以利用高效的内点算法进行求解。若表示式（6）的最优解，则式（3）的 最佳设计可以表示为

$h^{\left(k\right)}=Q^{-\frac{1}{2}}{G}^{\left(k\right)}{y}^{\left(k\right)}+Q^{-1}q---\left(7\right)$

其中，y^(k)中每个元素可以下式计算所得

$y_i^{\left(k\right)}=-\mathrm{sign}\left(v_i^{\left(k\right)}\right)\sqrt{z_i^{\left(k\right)}}---\left(8\right).$

设计示例：

设计方法的可行性可用以下两个示例进行验证。所有设计均采用归一化频 率，即二分之一采样频率对应于ω＝π。

示例1：示例1为一低通滤波器设计，表1为示例1所用的各项参数。图2 所示为低通滤波器的幅度响应，图3所示为通带内低通滤波器的群延时特性。 可以看到，滤波器的群延时在理想值附近呈现近似等波纹特性。

滤波器阶数N 25 通带 [0,0.3π] 阻带 [0.4π,π] 通带理想群延时 11

表1

滤波器阶数N 40 通带 [0.4π,0.5π] 阻带 [0,0.3π]∪[0.6π,π] 通带理想群延时 15

表2

示例2:示例2为一带通滤波器设计，表2为示例2所用的各项参数。图4所 示为带通滤波器的幅度响应，图5所示为通带内滤波器的群延时特性。可以看 到，通带内滤波器的群延时在理想值附近呈现近似等波纹特性。