用加权最小二乘法设计指定幅相 FIR 滤波器

用加权最小二乘法设计具有指定幅值和相位的 FIR 滤波器

线性相位滤波器的最小二乘设计是一个特例，其求解在 Matlab/Octave 函数 $\texttt{firls.m}$ 中实现。这里我要描述的是逼近复期望频率响应的一般情形，即所规定的响应不仅包括滤波器的幅值，还包括其相位响应。如果规定的相位响应是线性的，且相应的延迟与所选滤波器长度相匹配，那么所得到的滤波器将具有线性相位。否则，所设计的滤波器将具有非线性相位响应，以逼近规定的相位响应。

所述方法不仅能设计实系数滤波器，还能设计复系数 FIR 滤波器。这只需规定一个非共轭对称的复期望频率响应即可实现。令 $D(e^{j\omega})$ 表示复期望频率响应。若 $$D(e^{j\omega})\neq D^(e^{-j\omega})$$ 则得到复 FIR 滤波器，其中 $^$ 表示复共轭。若 $D(e^{j\omega})$ 共轭对称，即 $D(e^{j\omega}) = D^*(e^{-j\omega})$，则所得滤波器为实值。

总之，下文所述方法可用于设计具有实值或复值系数的线性与非线性相位滤波器。

在区间 $(-\pi,\pi]$ 内的密集频率网格 $\omega_i$，$i=1\ldots M$ 上定义复期望频率 $D(e^{j\omega})$。令 $H(e^{j\omega})$ 为系数为 $h[n]$ 的 FIR 滤波器的频率响应： $$H(e^{j\omega})=\sum_{n=0}^{N-1}h[n]e^{-jn\omega}$$ 其中 $N$ 为滤波器长度（即滤波器阶数等于 $N-1$）。滤波器设计问题现在可以简单地表述为 $$H(e^{j\omega_i})\overset{!}{=}D(e^{j\omega_i}),\quad i=1\ldots M\tag{1}$$ 其中 $\overset{!}{=}$ 表示在某种尚待指定的意义上“应尽可能接近”。我们假设 $(1)$ 是超定系统，即方程数多于未知数，换言之，频率点 $M$ 多于未知滤波器系数 $N$。

通过引入严格正权函数 $W(\omega)$ 来修改 $(1)$ 可能是合乎需要的。为在建立最小二乘问题时书写方便，我们在超定系统中使用 $W(\omega)$ 的平方根： $$\sqrt{W}(\omega)H(e^{j\omega_i})\overset{!}{=}\sqrt{W}(\omega)D(e^{j\omega_i}),\quad i=1\ldots M\tag{2}$$ 在最小化超定系统 $(2)$ 中的误差时，权重较大的频率区域所得的逼近误差会较小。

为建立加权最小二乘逼近问题，将 $(2)$ 写成矩阵形式是有利的。令 $\mathbf C$ 为复 $M\times N$ 矩阵，其行 ${\mathbf c}_i$ 由下式给出 $${\mathbf c}_i=\big[1,e^{-j\omega_i},e^{-2j\omega_i},\ldots,e^{-(N-1)j\omega_i}\big]$$ 将滤波器系数 $h[n]$ 排列成列向量 $\mathbf h$，则在频率 $\omega_i$ 处计算的滤波器频率响应可写为 $H(e^{j\omega_i})={\mathbf c}_i{\mathbf h}$。因此，元素为 $H(e^{j\omega_i})$，$i=1,\ldots,M$ 的 $M\times 1$ 向量由 ${\mathbf C}{\mathbf h}$ 给出。将频率 $\omega_i$ 处的复期望频率响应值排列成向量 $\mathbf d$，我们可以将 $(2)$ 写为 $${\mathbf W}^{1/2}{\mathbf C \mathbf h} \overset{!}{=} {\mathbf W}^{1/2}{\mathbf d}\tag{3}$$ 其中 $\mathbf W$ 是对角元素为 $W(\omega_i)$ 的对角矩阵。

在最小二乘意义上求解 $(3)$ 等价于最小化以下误差度量： $$\begin{align}\epsilon&=\sum_{i=1}^M W(\omega_i)\big|H(e^{j\omega_i})-D(e^{j\omega_i})\big|^2\&=({\mathbf C \mathbf h}-{\mathbf d})^H{\mathbf W}({\mathbf C \mathbf h}-{\mathbf d})\&={\mathbf h}^H{\mathbf C}^H{\mathbf W}{\mathbf C}{\mathbf h}-{\mathbf h}^H{\mathbf C}^H{\mathbf W}{\mathbf d}-{\mathbf d}^H{\mathbf W}{\mathbf C}{\mathbf h}+{\mathbf d}^H{\mathbf d}\tag{4}\end{align}$$ 其中我们用 $^H$ 表示复共轭转置。

最小化 $(4)$ 得到关于未知滤波器系数的如下 $N\times N$ 线性方程组： $${\mathbf C}^H{\mathbf W}{\mathbf C}{\mathbf h} = {\mathbf C}^H{\mathbf W}{\mathbf d} \tag{5}$$ 现在我们可以用某种线性方程求解器求解 $(5)$，也可以直接将超定系统 $(3)$ 作为 QR 求解器的输入。Matlab 和 Octave 通过反斜杠 \ 运算符同时提供这两种可能性。直接用 QR 求解器求解超定系统具有数值上的优势，而求解方系统 $(5)$ 的优势在于不需要存储矩阵，因为系统矩阵 ${\mathbf C}^H{\mathbf W}{\mathbf C}$ 具有 Toeplitz 结构，因而完全由其第一行确定。这样的系统可以用 Levinson 算法高效求解。

一般而言，$(3)$ 和 $(5)$ 的解是复值的。如果我们想设计实值滤波器，那么对每个频率点 $\omega_i>0$，必须存在点 $-\omega_i$ 满足 $D(e^{-j\omega_i})=D^*(e^{j\omega_i})$ 和 $W(-\omega_i)=W(\omega_i)$。另一种更高效的保证实值解的方法是稍微修改系统 $(3)$ 和 $(5)$。通过这些修改，我们只需在区间 $[0,\pi]$ 内指定频率点。超定系统 $(3)$ 变为 $$\left[\begin{array}{c}{\mathbf W}^{1/2}{\mathbf C}_R \{\mathbf W}^{1/2}{\mathbf C}_I \end{array} \right]{\mathbf h}\overset{!}{=} \left[\begin{array}{c}{\mathbf W}^{1/2}{\mathbf d}_R\{\mathbf W}^{1/2}{\mathbf d}_I\end{array}\right]\tag{6}$$ 其中下标 $R$ 和 $I$ 分别表示实部和虚部。$N\times N$ 系统 $(5)$ 变为 $$\textrm{Re}\left{{\mathbf C}^H{\mathbf W}{\mathbf C}\right}{\mathbf h} = \textrm{Re}\left{{\mathbf C}^H{\mathbf W}{\mathbf d}\right} \tag{7}$$

（原文较长，此处译出前 6000 字，完整内容请点「打开原文」。）