最小阶IIR多陷波滤波器设计

最小阶数全通型无限冲激响应多陷波滤波器的最小二乘设计

1 引言

数字多陷波滤波器在多种应用中都有需求，其中包括导航卫星系统接收机1、心电图2,3和肌电图4信号的预处理、频率估计5,6以及动态称重系统7 。多陷波滤波器的目的是抑制输入信号中某些特定的频谱分量，同时保持其他分量不变。虽然数字多陷波滤波器可以设计为有限冲激响应或无限冲激响应（IIR）滤波器，但与有限冲激响应滤波器相比，IIR多陷波滤波器的阶数显著更低。

众所周知，如果将IIR多陷波滤波器的极点配置在靠近单位圆的位置，则可以获得非常窄的凹口带宽。然而，这会导致瞬态响应持续时间8,9延长以及稳定性裕度10降低。文献中提出了多种IIR多陷波滤波器的设计方法，其中大多数可归类为以下三种方法之一：级联11,12、最优极点配置13–15和全通型16–19。

当考虑级联方法时，IIR多陷波滤波器的传递函数是通过级联多个IIR单陷波滤波器获得的。如果陷波频率之间间隔较远且凹口带宽较窄，则该方法已知能够提供良好的结果。否则，在陷波频率之间会出现不可控增益。8,17

第二种方法，即最优极点配置方法，将IIR多陷波滤波器设计问题表述为以未知极点位置为优化变量的优化问题。王等人14,15提出的方法假设极点半径已知，并利用粒子群优化14和Nelder‐Mead单纯形法15来求解未知的相位角；而在曾祥耀和裴士纲13的研究中，提出了一种迭代二次规划算法来获得未知的极点位置。最优极点配置方法的主要缺点是最大通带增益不相等。

关于全通型方法，首先将多陷波滤波器幅值响应的技术指标转换为相应全通滤波器相位响应的技术指标。然后，如果IIR多陷波滤波器的阶数等于陷波频率数量的两倍（这是IIR多陷波滤波器阶数的最小值），则建立一组线性方程并求解未知滤波器系数。16–18另一方面，如果IIR多陷波滤波器的阶数不是最小阶数，则可通过求解二次规划优化问题来确定未知滤波器系数。19最小阶IIR多陷波滤波器基于全通的方法在实际应用中取得成功的关键在于，这些方法具有简单的系数表达式。此外，基于全通的IIR多陷波滤波器可以采用高效的格型结构实现。

本文提出了一种新的基于最小阶全通的IIR多陷波滤波器的设计方法。该方法的创新之处在于同时利用了关于左截止频率和右截止频率的设计指标，同时严格满足陷波频率位置的技术指标。需要指出的是，佩伊和曾以及王和坎杜尔16,17提出的现有全通型设计方法虽然能够保证陷波频率位置的满足，但并未直接利用右截止频率或左截止频率。另一方面，王等人18提出的全通型设计方法由于假设所有极点具有相同的半径，因此既未直接利用左截止频率也未利用右截止频率。结果表明，与采用佩伊和曾及王和坎杜尔现有设计方法设计的基于最小階全通的IIR多陷波濾波器相比，使用提出的方法设计的滤波器具有更高的稳定性裕度，并且实现的左右截止频率相对于指定值的偏差更小。16,17

本文的其余部分结构如下。在第2节中，提出了全通型IIR多陷波滤波器的设计问题。第3节讨论了一种设计方法，第4节提供了与一些现有设计方法的比较。最后，第5节给出了结论。

2 问题描述

基于最小階全通的IIR多陷波濾波器的传递函数具有以下形式：

$$
HðzÞ¼ \frac{1}{2} \left(1＋AðzÞ\right), \quad ð1Þ
$$

其中A(z)是具有单极点的2K阶稳定全通滤波器

$$
AðzÞ＝z^{-2K} \frac{1＋\sum_{k＝1}^{2K} p_k z^k}{1＋\sum_{k＝1}^{2K} p_k z^{-k}},
\quad ð2Þ
$$

其中K是陷波频率的数量。将$ z= e^{j ω} $代入(1)式，经过一些数学推导后得到

$$
H\left(e^{jω}\right)＝ e^{j θðωÞ/2} \cdot \cos\left(\frac{θðωÞ}{2}\right), \quad ð3Þ
$$

其中 $ θ(ω) $ 是全通滤波器 $ A(z) $ 的相位响应：

$$
θðωÞ＝－2K ＋2 \arctan \left( \frac{\sum_{k＝1}^{2K} p_k \sin\left(kω\right)}{1＋\sum_{k＝1}^{2K} p_k \cos\left(kω\right)} \right):
\quad ð4Þ
$$

IIR多陷波滤波器的技术指标包括K个陷波频率，表示为$ ω_{n,i} $，$ i＝1,2,…,K $，以及在a分贝衰减处指定的K个凹口带宽 $ Δω_i $，其中$ i＝1,2,…,K $。左右截止频率与指定的凹口带宽相关。

$$
ω_{l,i}＝ ω_{n,i}－ \frac{Δω_i}{2},
\quad
ω_{r,i}＝ ω_{n,i}＋ \frac{Δω_i}{2},
\quad ð5Þ
$$

对于 $ i＝1、2、…、K $，分别地。因此，由于从(3)式导出的IIR多陷波滤波器的幅值响应等于

$$
\left| H\left(e^{jω}\right) \right|＝ \left| \cos\left(\frac{θðωÞ}{2}\right) \right|, \quad ð6Þ
$$

从稳定全通滤波器相位响应的单调递减特性（当 $ θð0Þ＝0 $和$ θðπÞ＝－2Kπ $时），11,16可以得出，如果 $ θ(ω) $满足，则关于陷波、左右截止频率的技术指标均被满足

$$
θ\left(ω_{l,i}\right)＝－2\left(i－1\right)π－ε, \quad ð7Þ
$$

$$
θ\left(ω_{n,i}\right)＝－\left(2i－1\right)π, \quad ð8Þ
$$

$$
θ\left(ω_{r,i}\right)＝－2iπ＋ε, \quad ð9Þ
$$

对于 $ i＝1，2，…，K $，其中

$$
ε＝ 2 \arccos\left(10^{－a/20}\right): \quad ð10Þ
$$

理想情况下，全通滤波器的系数向量，

$$
p ＝ \left[p_1\ p_2\ …\ p_{2K}\right]^T, \quad ð11Þ
$$

应确定使得由(7)–(9)给出的3K个方程得到满足。然而，这通常是不可能的，因为仅有2K个未知系数。

3 设计方法

全通滤波器的相位响应与其系数之间的关系为

$$
\sum_{k＝1}^{2K} p_k \sin\left(\frac{θðωÞ}{2} ＋ \left(K － k\right)ω\right)
＝－\sin\left(\frac{θðωÞ}{2} ＋Kω\right),
\quad ð12Þ
$$

由(4)可得，方程(7)–(9)可以重写为

$$
Ψ_L \cdot p ＝ γ_L , \quad ð13Þ
$$

$$
Ψ_N \cdot p ＝ γ_N , \quad ð14Þ
$$

$$
Ψ_R \cdot p ＝ γ_R, \quad ð15Þ
$$

其中 $ Ψ_L ＝ \left[ψ^{(L)} {ki}\right] $, $ Ψ_N ＝ \left[ψ^{(N)} {ki}\right] $和 $ Ψ_R ＝ \left[ψ^{(R)}_{ki}\right] $是 $ K×2K $ 矩阵，而 $ γ_L ＝ \left[γ^{(L)}_k\right] $, $ γ_N ＝ \left[γ^{(N)}_k\right] $和 $ γ_R ＝ \left[γ^{(R)}_k\right] $是具有元素的 $ K×1 $列矩阵

$$
ψ^{(L)} {ki} ＝ \sin\left(\left(ω {l,k} \left(K － i\right)－ \frac{ε}{2}\right)\right), \quad γ^{(L)} k ＝－\sin\left(Kω {l,k}－ \frac{ε}{2}\right),
$$

$$
ψ^{(N)} {ki} ＝ \cos\left(ω {n,k} \left(K － i\right)\right), \quad γ^{(N)} k ＝－\cos\left(Kω {n,k}\right),
$$

$$
ψ^{(R)} {ki} ＝ \sin\left(\left(ω {r,k} \left(K － i\right)＋ \frac{ε}{2}\right)\right), \quad γ^{(R)} k ＝－\sin\left(Kω {r,k}＋ \frac{ε}{2}\right): \quad ð16Þ
$$

显然，由矩阵表示法(13)–(15)给出的全部3K个线性方程无法同时满足，因为只有2K个未知数。因此，提出以下设计方法：

1. 首先，为了确保精确满足关于陷波频率位置的技术指标，将(14)代入(13)和(15)中。换句话说，进行变量消元。因此，方程(13)和(15)变为

$$
\left(Υ_L－Φ_L \cdot Φ^{-1}_N \cdot Υ_N\right) \cdot p_x ＝ γ_L－Φ_L \cdot Φ^{-1}_N \cdot γ_N, \quad ð17Þ
$$

$$
\left(Υ_R－Φ_R \cdot Φ^{-1}_N \cdot Υ_N\right) \cdot p_x ＝ γ_R－Φ_R \cdot Φ^{-1}_N \cdot γ_N, \quad ð18Þ
$$

其中 $ Υ_L、Υ_N、Υ_R、Φ_L、Φ_N $和$ Φ_R $是与 $ Ψ_L、Ψ_N $和 $ Ψ_R $相关的 $ K×K $矩阵

$$
Ψ_L ＝ \left[Φ_L\ Υ_L\right], \quad Ψ_N ＝ \left[Φ_N\ Υ_N\right], \quad Ψ_R ＝ \left[Φ_R\ Υ_R\right], \quad ð19Þ
$$

而 $ p $被分块为

$$
p ＝ \begin{bmatrix} p_e \ p_x \end{bmatrix}, \quad ð20Þ
$$

且$ p_e $是被消去的系数的 $ K×1 $向量

$$
p_e ＝Φ^{-1}_N \cdot \left( γ_N －Υ_N \cdot p_x \right), \quad ð21Þ
$$

注意(14)。

2. 其次，求解由(17)和(18)给出的 $ K $个未知数的 $ 2K $个方程组，该方程组在最小二乘意义下求解：

$$
p_x ＝
\begin{bmatrix}
Υ_L－Φ_L \cdot Φ^{-1}_N \cdot Υ_N \
Υ_R－Φ_R \cdot Φ^{-1}_N \cdot Υ_N
\end{bmatrix}^†
\cdot
\begin{bmatrix}
γ_L－Φ_L \cdot Φ^{-1}_N \cdot γ_N \
γ_L－Φ_L \cdot Φ^{-1}_N \cdot γ_N
\end{bmatrix}
, \quad ð22Þ
$$

其中 $ † $表示摩尔‐彭罗斯逆。因此，左右截止频率的技术指标被近似满足。

当 IIR 多陷波滤波器的设计指标满足 $ ω_{n,k} ＝ π－ω_{n,K－k＋1} $ 和 $ Δω_k ＝Δω_{K－k＋1} $，其中 $ k $为 $ ＝ 1,2,…,K $，即设计指标关于 $ ω＝ π/2 $ 具有对称性时，会出现一个有趣的情况。在这种情况下，使用提出的 IIR 多陷波滤波器设计方法可得到 $ p_1 ＝p_2 ＝…＝ p_{2K－1} ＝ 0 $；也就是说，该 IIR 多陷波滤波器仅需要 K 次乘法。

4 与现有设计方法的比较

所提出的设计方法与现有的基于最小阶全通的IIR多陷波滤波器设计方法（佩伊和曾、王和坎杜尔 ^16,17 ）进行了比较，这些方法均能保证满足陷波频率位置的要求。虽然这些方法考虑了3分贝陷波带宽，但可以推广至在任意衰减值下指定凹口带宽的情况。需要注意的是，对于王等人 ¹⁸ 提出的设计方法，这种推广是不可能的，因为其极点半径是从3分贝陷波带宽推导得出的，并且假设相同的陷波带宽会导致所有极点具有相同的半径。以下简要描述用于与提出的方法进行比较的佩伊和曾以及王和坎杜尔 ^16,17 的方法。

1. 使用王和坎杜尔的方法I ¹⁷ 可确保精确满足关于凹口位置和左截止频率的技术指标。该方法与佩伊和曾提出的方法相同 ¹⁶ ；然而，克服了有关正切运算的限制。实际上，对王和坎杜尔中描述方法I的方程进行一些数学推导后 ¹⁷ ，便可得到方程(13)和(14)。因此，利用方法I的推广形式在指定凹口带宽的任意衰减值处获得的传递函数所对应的向量p，由下式确定

$$
p ＝ \left[\Psi_L;\ \Psi_N\right]^{-1} \cdot \left[\gamma_L;\ \gamma_N\right]. \quad ð23Þ
$$

2. 采用王和坎杜尔的方法II ¹⁷ 可确保精确满足关于凹口位置和右截止频率的技术指标。与前述情况类似，这表明：通过将方法II ¹⁷ 推广至指定凹口带宽的任意衰减值所得到的传递函数对应的向量p，可由方程组(14)和(15)的解确定，即

$$
p ＝ \left[\Psi_N;\ \Psi_R\right]^{-1} \cdot \left[\gamma_N;\ \gamma_R\right]. \quad ð24Þ
$$

虽然王和坎杜尔的方法I ¹⁷ 没有直接利用右截止频率，即(15)，但王和坎杜尔的方法II ¹⁷ 没有直接利用左截止频率，即(13)。从这个意义上说，可以将所提出的设计方法视为王和坎杜尔 ¹⁷ 中“缺失”的方法，因为在设计过程中同时使用了(13)和(15)。

使用提出的方法以及佩伊和曾、王和坎杜尔的现有全通型设计方法设计的滤波器，在最大极点半径 $ \rho_{\text{max}} $、实现的左截止频率和右截止频率相对于指定值 $ \delta\omega_{l,i}＝ (\omega^0_{l,i}/\omega_{l,i}－1) \times 100\% $ 和 $ \delta\omega_{r,i}＝ (\omega^0_{r,i}/\omega_{r,i}－1) \times 100\% $ 的相对偏差 $ \omega^0_{l,i} $ 与 $ \omega^0_{r,i} $，以及实现的凹口带宽 $ \Delta\omega^0_i＝ \omega^0_{r,i}－\omega^0_{l,i} $ 等方面进行了比较。注意，若 $ \delta\omega_{l,i}≥0 $ 或 $ \delta\omega_{r,i}≤0 $，则第 $ i $ 个左或右截止频率处的技术指标被精确满足或超额满足；否则，技术指标未充分满足。另一方面，较低的最大极点半径可带来更短的瞬态响应持续时间 ^8,9 和更高的稳定性裕度 ¹⁰ 。显然，若采用方法I设计滤波器，则对于 $ i＝1,2,…,K $，有 $ \delta\omega_{l,i}＝0 $；而若采用方法II ¹⁷ ，则对于 $ i＝1,2,…,K $，有 $ \delta\omega_{r,i}＝0 $。

以第一个例子为例，考虑以下IIR多陷波滤波器的技术指标：$ \omega_{n,1} ＝ 0.3\pi, \omega_{n,2} ＝ 0.7\pi, \Delta\omega_1 ＝ \Delta\omega_2 ＝ 0.1\pi $，$ a＝2 $ dB。使用所提出的方法和王和坎杜尔的现有设计方法设计的IIR多凹口滤波器的实现的左右截止频率相对于指定值的相对偏差（百分比），$ \delta\omega_{l,i} $和 $ \delta\omega_{r,i} $，以及实现的陷波带宽 $ \Delta\omega^0_i $ （$ i＝1，2 $）和最大极点半径如表1所示。可以看出，当使用提出的方法时，第一左截止频率和第二右截止频率的技术指标分别欠满足0.9%和0.3%；而当考虑方法I和方法II ¹⁷ 时，第二右截止频率和第一左截止频率分别欠满足0.58%和1.75%。因此，对于第一个例子，使用方法I ¹⁷ 设计的IIR多陷波滤波器优于其他两个滤波器。另一方面，使用所提出的设计方法设计的滤波器具有最高的稳定性裕度，即最低的最大极点半径，并且优于使用方法II设计的滤波器 ¹⁷ 。此外，使用所提出的设计方法设计的滤波器的实现的陷波带宽窄于指定值，而在使用王和坎杜尔的设计方法时则并非如此。

作为第二个例子，考虑以下IIR多陷波滤波器的技术指标：$ \omega_{n,1} ＝ 0.2\pi, \omega_{n,2} ＝ 0.4\pi, \omega_{n,3} ＝ 0.7\pi, \Delta\omega_1 ＝\Delta\omega_2 ＝ \Delta\omega_3 ＝ 0.1\pi $，$ a＝2.2 $ dB。$ \delta\omega_{l,i},\delta\omega_{r,i} $以及$ {i=1,2,3} $时的 $ \Delta\omega^0_i $，对应于使用提出的方法和王和坎杜尔的现有设计方法设计的IIR多凹口滤波器 ¹⁷ ，如表2所示。对于该组技术指标，使用所提出的设计方法设计的滤波器优于使用王和坎杜尔现有方法设计的两种滤波器 ¹⁷ ，因为在截止频率处未满足规格要求的相对偏差小于0.85%，而分别采用方法I和方法II ¹⁷ 时该偏差为2.58%和6.09%。此外，使用所提出的设计方法获得的滤波器还具有最高的稳定性裕度，其实现的陷波带宽比规定的带宽更窄，而王和坎杜尔的设计方法 ¹⁷ 则不具备这一特性。

作为第三个例子，考虑以下IIR多陷波滤波器的技术指标：$ \omega_{n,1} ＝ 0.1, \omega_{n,2} ＝ 0.2\pi, \omega_{n,3} ＝ 0.4\pi, \omega_{n,4} ＝ 0.8\pi, \Delta\omega_1 ＝\Delta\omega_2 ＝ 0.06\pi, \Delta\omega_3 ＝ 0.08\pi, \Delta\omega_4 ＝ 0.1\pi $，$ a＝3 $分贝。$ \delta\omega_{l,i}, \delta\omega_{r,i} $以及 $ \Delta\omega^0_i $ （其中$ i=1,2,3,4 $），这些值对应于使用王和坎杜尔提出的和现有的设计方法设计的IIR多凹口滤波器，其数据列于表3中。对于所考虑的IIR多陷波滤波器技术指标，所提出的设计方法远优于其他两种方法，因为仅在第四个左截止频率处有0.02%的技术指标未满足要求。此外，使用所提出的设计方法设计的滤波器的实现的陷波带宽比指定的凹口带宽更窄或相等，而使用王和坎杜尔的设计方法时则并非如此。另一方面，其他两个滤波器在多个地方均未能满足技术指标要求。

表1 示例1：使用提出的方法和王和坎杜尔的现有设计方法设计的IIR多凹口滤波器，其实现的左右截止频率相对于指定值的相对偏差（百分比）、实现的陷波带宽以及最大极点半径 ¹⁷

	提出的方法, $ \rho_{\text{max}} ＝ 0.8814 $	方法I 17 , $ \rho_{\text{max}} ＝ 0.8875 $	方法II 17 , $ \rho_{\text{max}} ＝0.8875 $
i	$ \delta\omega_{l,i} $	$ \delta\omega_{r,i} $	$ \Delta\omega^0_i/\pi $
1	-0.90	-0.78	0.0995
2	0.42	0.30	0.0995

表2 示例2：使用提出的方法和王和坎杜尔的现有设计方法设计的IIR多凹口滤波器，其实现的左右截止频率相对于指定值的相对偏差（百分比）、实现的凹口带宽以及最大极点半径 ¹⁷

	提出的方法, $ \rho_{\text{max}} ＝ 0.8811 $	方法I 17 , $ \rho_{\text{max}} ＝ 0.8855 $	方法II 17 , $ \rho_{\text{max}} ＝0.8929 $
i	$ \delta\omega_{l,i} $	$ \delta\omega_{r,i} $	$ \Delta\omega^0_i/\pi $
1	-0.13	-4.28	0.0895
2	3.54	-0.37	0.0859
3	1.25	0.85	0.0982

表3 示例3：使用提出的方法和王和坎杜尔的现有设计方法设计的IIR多凹口滤波器，其实现的左右截止频率相对于指定值的相对偏差（百分比）、实现的凹口带宽以及最大极点半径 ¹⁷

	提出的方法, $ \rho_{\text{max}} ＝0.9396 $	方法I 17 , $ \rho_{\text{max}} ＝ 0.9088 $	方法II 17 , $ \rho_{\text{max}} ＝0.9287 $
i	$ \delta\omega_{l,i} $	$ \delta\omega_{r,i} $	$ \Delta\omega^0_i/\pi $
1	13.92	-11.11	0.0358
2	6.94	-3.86	0.0393
3	1.34	-0.42	0.0770
4	-0.02	-0.02	0.1000

三个截止频率，见表3。然而，在此示例中，使用所提出的设计方法得到的滤波器具有最低稳定性裕度。所得滤波器的幅频响应（以分贝为单位）如图1所示。

5 结论

本文提出了基于最小阶全通的IIR多陷波滤波器的最小二乘设计方法。首先将多陷波滤波器幅值响应的技术指标表示为线性等式约束，并通过变量消元确保陷波频率位置的精确满足。然后在最小二乘意义下求解剩余的超定线性方程组。与现有的全通型设计方法（佩伊和曾、王和坎杜尔） ^16,17 不同，所提出的设计方法同时利用了左截止频率和右截止频率的技术指标。与现有最小阶全通IIR多陷波滤波器设计方法相比，结果表明，采用提出的方法可以获得具有更高稳定性裕度（即更低的极点最大半径）、实现的左右截止频率相对于指定值的相对偏差更小，以及实现的陷波带宽比指定值更窄的传递函数。