Sam's Blog

English

简体中文

English

简体中文

Appearance

第十二章：最小平方滤波

About 1441 wordsAbout 5 min

2026-01-29

Note

最小平方滤波是一种最佳滤波方法：即按最小平方准则（Least-Squares Criterion）设计滤波器，使滤波输出与某期望信号在误差能量意义下最接近。方法具有简单、灵活、有效的特点，是信号处理中的基本滤波方法之一。

§1 最小平方滤波

1.1 最小平方滤波的数学模型

最小平方滤波建模如下：

  • 输入信号：xtx_t
  • 滤波因子：ht=(h−m0,h−m0+1,…,h−m0+m)h_t = (h_{-m_0}, h_{-m_0+1}, \dots, h_{-m_0+m})，共 m+1m+1 项
  • 实际输出：

    yt=ht∗xt=∑s=−m0−m0+mhsxt−sy_t = h_t * x_t = \sum_{s=-m_0}^{-m_0+m} h_s x_{t-s}

  • 期望输出：ztz_t
  • 输出误差：Et=yt−ztE_t = y_t - z_t
  • 误差能量：

    Q=∑t=−∞∞(yt−zt)2Q = \sum_{t=-\infty}^{\infty} (y_t - z_t)^2

  • 相对误差能量：

    E=Q∑t=−∞∞zt2E = \frac{Q}{\sum_{t=-\infty}^{\infty} z_t^2}

最小平方滤波准则

选择滤波因子 hth_t，使误差能量 QQ 最小化。满足该准则的 hth_t 称为最小平方滤波因子。

Note

以上模型中 xt,ht,ztx_t, h_t, z_t 均为实信号。

1.2 最小平方滤波方程

令误差能量 QQ 对各 hnh_n 偏导为零，得到极值条件：

∑s=−m0−m0+mhsrxx(n−s)=rzx(n),n=−m0,−m0+1,…,−m0+m\sum_{s=-m_0}^{-m_0+m} h_s r_{xx}(n-s) = r_{zx}(n), \quad n = -m_0, -m_0+1, \dots, -m_0+m

其中：

  • rxx(t)=∑k=−∞∞xkxk−tr_{xx}(t) = \sum_{k=-\infty}^{\infty} x_k x_{k-t} 是 xtx_t 的自相关函数；
  • rzx(t)=∑k=−∞∞zkxk−tr_{zx}(t) = \sum_{k=-\infty}^{\infty} z_k x_{k-t} 是 ztz_t 与 xtx_t 的互相关函数。

矩阵形式（陶布利兹方程）：

[rxx(0)rxx(1)⋯rxx(m)rxx(1)rxx(0)⋯rxx(m−1)⋮⋮⋱⋮rxx(m)rxx(m−1)⋯rxx(0)][h−m0h−m0+1⋮h−m0+m]=[rzx(−m0)rzx(−m0+1)⋮rzx(−m0+m)]\begin{bmatrix} r_{xx}(0) & r_{xx}(1) & \cdots & r_{xx}(m) \\ r_{xx}(1) & r_{xx}(0) & \cdots & r_{xx}(m-1) \\ \vdots & \vdots & \ddots & \vdots \\ r_{xx}(m) & r_{xx}(m-1) & \cdots & r_{xx}(0) \end{bmatrix} \begin{bmatrix} h_{-m_0} \\ h_{-m_0+1} \\ \vdots \\ h_{-m_0+m} \end{bmatrix} = \begin{bmatrix} r_{zx}(-m_0) \\ r_{zx}(-m_0+1) \\ \vdots \\ r_{zx}(-m_0+m) \end{bmatrix}

该方程的解称为最小平方滤波因子。

频域形式（无限长情形）

若允许 hth_t 为无限长序列，设其频谱为 H(f)H(f)，rxx(t),rzx(t)r_{xx}(t), r_{zx}(t) 的频谱分别为 Rxx(f),Rzx(f)R_{xx}(f), R_{zx}(f)，则频域关系为：

H(f)Rxx(f)=Rzx(f)H(f) R_{xx}(f) = R_{zx}(f)

1.3 最小平方滤波的误差能量与相对误差能量

记 Qmin⁡Q_{\min} 为最小误差能量，则有三种等价表达式：

  1. 用 实际输出与期望输出 表示：

    Qmin⁡=∑t=−∞∞zt2−∑t=−∞∞yt2Q_{\min} = \sum_{t=-\infty}^{\infty} z_t^2 - \sum_{t=-\infty}^{\infty} y_t^2

  2. 用 期望输出与滤波因子 表示：

    Qmin⁡=∑t=−∞∞zt2−2∑s=−m0−m0+mhsrzx(s)Q_{\min} = \sum_{t=-\infty}^{\infty} z_t^2 - 2 \sum_{s=-m_0}^{-m_0+m} h_s r_{zx}(s)

  3. 用 输入输出关系 表示：

    Qmin⁡=∑t=−∞∞(yt−zt)ytQ_{\min} = \sum_{t=-\infty}^{\infty} (y_t - z_t) y_t

相对误差能量：

Emin⁡=Qmin⁡∑zt2=1−∑yt2∑zt2E_{\min} = \frac{Q_{\min}}{\sum z_t^2} = 1 - \frac{\sum y_t^2}{\sum z_t^2}

性质：

  • 0≤Emin⁡≤10 \le E_{\min} \le 1
  • Emin⁡=0  ⟺  yt=ztE_{\min}=0 \iff y_t = z_t（完美匹配）
  • Emin⁡=1  ⟺  yt≡0E_{\min}=1 \iff y_t \equiv 0（完全失配，常对应 ht=0h_t=0）

1.4 解最小平方滤波方程中的白噪化方法

直接求解最小平方方程在某些情况下效果不佳，解可能剧烈振荡，滤波结果远离期望输出。

将自相关函数 rbb(n)r_{bb}(n) 修改为：

r~bb(n)={(1+λ)rbb(0),n=0rbb(n),n≠0\tilde{r}_{bb}(n) = \begin{cases} (1+\lambda) r_{bb}(0), & n = 0 \\ r_{bb}(n), & n \ne 0 \end{cases}

其中 λ>0\lambda > 0 为白噪因子（常用 0<λ<0.20 < \lambda < 0.2）。

  • 物理意义：向信号添加一“白噪声背景”，即添加一个常量频谱项 λrbb(0)\lambda r_{bb}(0)，其频谱为常数（平谱），类比白光/白噪声。
  • 效果：改善方程系数矩阵的条件数，使解更稳定、平滑，实际滤波效果显著提升。

§2 最小平方滤波方程的解法及解的性质

最小平方滤波方程是陶布利兹方程（Toeplitz Equation）——一种系数矩阵为陶布利兹矩阵（沿对角线平移不变）的线性方程组。因其特殊结构，有高效递推解法。

2.1 陶布利兹方程及其特点

一般形式（mm 阶）：

[r1r2⋯rmr2r1⋯rm−1⋮⋮⋱⋮rmrm−1⋯r1][x1x2⋮xm]=[g1g2⋮gm]\begin{bmatrix} r_1 & r_2 & \cdots & r_m \\ r_2 & r_1 & \cdots & r_{m-1} \\ \vdots & \vdots & \ddots & \vdots \\ r_m & r_{m-1} & \cdots & r_1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_m \end{bmatrix} = \begin{bmatrix} g_1 \\ g_2 \\ \vdots \\ g_m \end{bmatrix}

特点：已知 (x1,…,xm)(x_1, \dots, x_m) 是上面方程的一组解，则 (xm,…,x1)(x_m, \dots, x_1) 是下面方程的解。

[r1r2⋯rmr2r1⋯rm−1⋮⋮⋱⋮rmrm−1⋯r1][xmxm−1⋮x1]=[gmgm−1⋮g1]\begin{bmatrix} r_1 & r_2 & \cdots & r_m \\ r_2 & r_1 & \cdots & r_{m-1} \\ \vdots & \vdots & \ddots & \vdots \\ r_m & r_{m-1} & \cdots & r_1 \end{bmatrix} \begin{bmatrix} x_m \\ x_{m-1} \\ \vdots \\ x_1 \end{bmatrix} = \begin{bmatrix} g_m \\ g_{m-1} \\ \vdots \\ g_1 \end{bmatrix}

2.2 特殊的陶布利兹方程的递推解法

考虑右侧为单位向量 [1,0,…,0]T[1, 0, \dots, 0]^T 的情形：

∑k=1nr∣j−k∣ak(n)=δj1,j=1,2,…,n\sum_{k=1}^{n} r_{|j-k|} a_k^{(n)} = \delta_{j1}, \quad j=1,2,\dots,n

[r1r2⋯rnr2r1⋯rn−1⋮⋮⋱⋮rnrn−1⋯r1][a1(n)a2(n)⋮an(n)]=[10⋮0]\begin{bmatrix} r_1 & r_2 & \cdots & r_n \\ r_2 & r_1 & \cdots & r_{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ r_n & r_{n-1} & \cdots & r_1 \end{bmatrix} \begin{bmatrix} a_1^{(n)} \\ a_2^{(n)} \\ \vdots \\ a_n^{(n)} \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix}

递推公式（Levinson 递推）：

  • 初值：a1(1)=1/r1a_1^{(1)} = 1 / r_1
  • 对 n=1,2,…,m−1n = 1, 2, \dots, m-1：
    • 计算反射系数：

      pn=∑k=1nrn+1−kak(n)p_n = \sum_{k=1}^{n} r_{n+1-k} a_k^{(n)}

    • 更新解：

      ak(n+1)=11−pn2{ak(n)−pnan+1−k(n),1≤k≤n−pn,k=n+1a_k^{(n+1)} = \frac{1}{1-p_n^2} \begin{cases} a_k^{(n)} - p_n a_{n+1-k}^{(n)}, & 1 \le k \le n \\ -p_n, & k = n+1 \end{cases}

该递推法复杂度为 O(N2)O(N^2)，远优于直接求逆的 O(N3)O(N^3)。

2.3 一般的陶布利兹方程的递推解法

对任意右端 [g1,…,gm]T[g_1, \dots, g_m]^T：

  • 仍采用 Levinson 递推框架；
  • 除计算 ak(n)a_k^{(n)} 外，同步计算主解 xk(n)x_k^{(n)}；
  • 初值：x1(1)=g1/r1x_1^{(1)} = g_1 / r_1
  • 递推：

    xk(n+1)=xk(n)+(gn+1−βn)an+2−k(n+1),k=1,…,n+1x_k^{(n+1)} = x_k^{(n)} + (g_{n+1} - \beta_n) a_{n+2-k}^{(n+1)}, \quad k = 1,\dots,n+1

    其中 βn=∑k=1nrn+1−kxk(n)\beta_n = \sum_{k=1}^{n} r_{n+1-k} x_k^{(n)}

2.4 最小平方滤波方程解的性质

对以下形式的特殊方程

∑s=1nr∣j−s∣as(n)=δj1\sum_{s=1}^{n} r_{|j-s|} a_s^{(n)} = \delta_{j1}

其解 a(n)=(a1(n),…,an(n))a^{(n)} = (a_1^{(n)}, \dots, a_n^{(n)}) 为最小相位信号（即其 Z 变换多项式在单位圆内及圆上无零点）。

Note

一般方程的解不一定是最小相位的。

反例：若右端为 [rm,rm−1,…,r1]T[r_m, r_{m-1}, \dots, r_1]^T，则解为 (0,…,0,1)(0,\dots,0,1)，显然非最小相位。

Copyright Ownership:Staaaaaaaaar

License under:Attribution-NonCommercial-ShareAlike 4.0 International (CC-BY-NC-SA-4.0)