自适应滤波：从 LMS 到 NLMS

LMS 算法

自适应滤波的核心问题是：给定参考信号 x(n) 和期望信号 d(n)，找到一个滤波器系数向量 w，使输出 y(n) 逼近 d(n)。

最小均方（Least Mean Square, LMS）算法是解决这个问题最经典的方法，由 Widrow 和 Hoff 在 1960 年提出。核心思路是每次迭代都沿着误差曲面最陡的方向下降一步——也就是随机梯度下降。

数学模型

滤波器长度为 N，定义：

• 权重向量：w(n) = [w_0(n), w_1(n), ..., w_{N-1}(n)]^T
• 参考信号向量：x(n) = [x(n), x(n-1), ..., x(n-N+1)]^T

滤波器输出是两者的内积：

$$ y(n) = w^T(n) x(n) = \sum_{k=0}^{N-1} w_k(n) x(n-k) $$

误差信号：

$$ e(n) = d(n) - y(n) $$

代价函数为误差的均方值：

$$ J(n) = E[e^2(n)] $$

对 w 求梯度，得到权重更新规则：

$$ \nabla J(n) = -2 e(n) x(n) $$

$$ w(n+1) = w(n) - \mu \nabla J(n) = w(n) + \mu e(n) x(n) $$

其中 \mu 是步长因子，控制收敛速度。

自适应滤波器结构

下图展示了自适应滤波器的闭环工作流程：

flowchart TD
    REF["参考信号 x(n)"] --> W["自适应滤波器 W(z)<br/>权重 w(n)"]
    W -->|"输出 y(n)"| SUM(("⊕"))
    DES["期望信号 d(n)<br/>(误差麦克风)"] --> SUM
    SUM -->|"误差 e(n) = d(n) - y(n)"| UPDATE["权重更新<br/>w(n+1) = w(n) + μ·e(n)·x(n)"]
    REF --> UPDATE
    UPDATE --> W

    classDef signal fill:#2196F3,color:#fff
    classDef proc fill:#9C27B0,color:#fff
    classDef err fill:#f44336,color:#fff
    class REF,DES signal
    class W,UPDATE,SUM proc

参考信号 x(n) 经过自适应滤波器得到输出 y(n)，与期望信号 d(n) 比较产生误差 e(n)，误差再反馈回来指导权重更新。这个闭环持续运行，误差不断减小，直到滤波器收敛到最优解。

公式逐项拆解

LMS 的权重更新公式看起来不过一行，但每一项都有明确的物理含义：

$$ w(n+1) = w(n) + \mu \cdot e(n) \cdot x(n) $$

一个直观的类比：蒙眼下山

想象你被蒙住眼睛站在山坡上，目标是走到山谷最低点：

• e(n) 是你脚下坡度的感受——坡度越陡，误差越大，说明离谷底越远
• μ 是你迈步的幅度——步子太大可能跨过山谷冲到对面坡上（发散），步子太小半天走不到底（收敛慢）
• x(n) 是你脚下的地形方向——告诉你往哪边迈脚
• 每次更新就是：感觉坡度 → 朝正确方向迈一步 → 再感觉 → 再迈一步

数值例子：假设 μ = 0.001，当前 e = 0.5，x = 0.8，则更新量 = 0.001 × 0.5 × 0.8 = 0.0004。如果 e 增大到 5.0，更新量 = 0.001 × 5.0 × 0.8 = 0.004——误差越大，调整越猛。

为什么叫"随机梯度下降"？

严格来说，LMS 的代价函数是均方误差 J(n) = E[e²(n)]，其真实梯度是 ∇J = -2 · E[e(n) · x(n)]。但 LMS 不做期望运算——它直接用当前时刻的瞬时值 -2 · e(n) · x(n) 作为梯度的估计。

这相当于用 一个样本的误差 代替了 全体样本的期望，所以叫"随机"（stochastic）。虽然单次估计有噪声，但多次迭代平均下来，权重的走向与真正的梯度下降一致。这种"以瞬时代期望"的思路，是 LMS 计算量小却有效的核心原因。

收敛条件

LMS 收敛的充分必要条件是步长满足：

$$ 0 < \mu < \frac{2}{N \cdot P_x} $$

其中 P_x = E[x^2(n)] 是输入信号功率。步长越大收敛越快，但超出上限会发散。

收敛条件拆解

0 < μ < 2 / (N · P_x) 这个不等式决定了 LMS 能否正常工作：

• N：滤波器长度。抽头越多，每个抽头分到的"能量"越少，μ 的上限就越低。64 阶滤波器的 μ 上限只有 32 阶的一半。
• P_x：输入信号功率。信号越强，μ 必须越小，否则梯度噪声会放大到失控。
• 取值的直觉：μ 太大→权重一步跨过最优值→来回震荡甚至发散；μ 太小→每次只挪一点点→收敛极慢。

数值例子：假设 N = 64，P_x = 1.0（单位功率信号），则：

1

μ < 2 / (64 × 1.0) = 0.03125

实践中一般取理论上限的 1/10 到 1/2，也就是 0.003 到 0.015，留足安全余量。

收敛过程

LMS 的收敛是一个从"粗调"到"微调"的过程：

• 初始阶段：权重为零（或随机值），误差很大，每次更新的步子也大
• 学习阶段：权重逐渐接近最优值，误差快速下降
• 稳态阶段：权重在最优值附近轻微抖动（称为"失调噪声"），误差不再明显下降

下图展示了这个迭代循环：

flowchart TD
    A["初始化<br/>w(0) = 0"] --> B["处理样本<br/>y = w^T · x"]
    B --> C["计算误差<br/>e = d - y"]
    C --> D["更新权重<br/>w += μ·e·x"]
    D --> E{"收敛?"}
    E -->|"否<br/>误差仍大"| B
    E -->|"是<br/>误差足够小"| F["稳态运行<br/>持续微调"]

    classDef init fill:#FF9800,color:#fff
    classDef loop fill:#2196F3,color:#fff
    classDef done fill:#4CAF50,color:#fff
    class A init
    class B,C,D loop
    class F done

如果画出误差随时间的变化曲线，会看到一条从高处快速下降、然后逐渐平缓的曲线——像滑梯一样，先陡后缓。

下面的动画画出了这条“滑梯”曲线，可以清楚看到三个阶段的过渡：

C 实现

环形缓冲区管理参考信号延迟线，避免每次迭代都移动数据：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

#include <stdint.h>
#include <string.h>
#include <math.h>

#define FILTER_ORDER 64
#define MU 0.001f

typedef struct {
    float w[FILTER_ORDER];
    float x[FILTER_ORDER];
    uint32_t ptr;
} LMS_Filter;

void lms_init(LMS_Filter *f) {
    memset(f->w, 0, sizeof(f->w));
    memset(f->x, 0, sizeof(f->x));
    f->ptr = 0;
}

static float dot_product(const float *a, const float *b, uint32_t n) {
    float sum = 0.0f;
    for (uint32_t i = 0; i < n; i++) sum += a[i] * b[i];
    return sum;
}

float lms_process(LMS_Filter *f, float ref, float err) {
    f->x[f->ptr] = ref;
    float y = 0.0f;
    for (uint32_t i = 0; i < FILTER_ORDER; i++) {
        uint32_t idx = (f->ptr + FILTER_ORDER - i) % FILTER_ORDER;
        y += f->w[i] * f->x[idx];
    }
    for (uint32_t i = 0; i < FILTER_ORDER; i++) {
        uint32_t idx = (f->ptr + FILTER_ORDER - i) % FILTER_ORDER;
        f->w[i] += MU * err * f->x[idx];
    }
    f->ptr = (f->ptr + 1) % FILTER_ORDER;
    return y;
}

lms_process 先写参考信号到环形缓冲区，计算滤波器输出，然后对每个权重执行 w += mu * e * x 更新。ptr 循环移动，等效于一个滑动窗。

NLMS 算法

LMS 有一个实际问题：步长 \mu 对输入信号的幅度非常敏感。信号幅度大时，梯度噪声被放大，可能导致发散；幅度小时收敛又太慢。每次换场景都要重新调 \mu。

归一化 LMS（Normalized LMS, NLMS）解决了这个问题——把步长除以输入信号的瞬时功率，让更新量自动适应信号幅度。

权重更新

NLMS 的权重更新公式：

$$ w(n+1) = w(n) + \frac{\mu}{x^T(n) x(n) + \varepsilon} e(n) x(n) $$

分母 x^T(n) x(n) 是参考信号向量的能量（即各采样点平方和），\varepsilon 是一个很小的正常数，防止分母为零。

NLMS 更新公式拆解

NLMS 的更新公式可以理解为 LMS 的"自适应步长"版本：

$$ w(n+1) = w(n) + \frac{\mu}{x^T(n) x(n) + \varepsilon} \cdot e(n) \cdot x(n) $$

数值例子：假设 μ = 0.5，输入信号 x = [0.1, -0.05, 0.2, ...]，能量 x^T x = 0.1² + (-0.05)² + 0.2² + ... = 0.0525，则等效步长 = 0.5 / (0.0525 + 1e-6) ≈ 9.52。如果信号增强到 x = [0.5, -0.3, 0.8, ...]，能量变为约 0.98，等效步长 ≈ 0.5 / 0.98 ≈ 0.51——信号强时步长自动缩小。

对比 LMS 和 NLMS 的更新量：

NLMS 相当于每一步都根据信号能量调整了等效步长。信号强时步长自动缩小，信号弱时步长放大，梯度噪声水平保持稳定。

步长范围

NLMS 的步长 \mu 在 (0, 2) 范围内保证收敛。相比 LMS 需要根据信号功率反复调试，NLMS 在实际系统中好用得多。

LMS vs NLMS 怎么选？

一句话经验：先试 NLMS，参数调不通再切 LMS。 因为 NLMS 参数少、鲁棒性好，大多数场景一次调好。

NLMS 定点实现

嵌入式系统没有 FPU 时，用 Q31 定点数实现 NLMS。权重和信号都用 int32_t 表示，乘法用 q31_mul 宏做饱和移位：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

#include <stdint.h>
#include <string.h>

#define FILTER_ORDER 32
#define MU_Q31  0x147AE148
#define EPS_Q31 0x0001869F

static inline int32_t q31_mul(int32_t a, int32_t b) {
    return (int32_t)(((int64_t)a * (int64_t)b) >> 31);
}

typedef struct {
    int32_t h[FILTER_ORDER];
    int32_t x[FILTER_ORDER];
    uint16_t ptr;
} NLMS_Filter_Fixed;

定点版本计算能量时用 q31_mul 累加，避免 32 位溢出。MU_Q31 和 EPS_Q31 采用 Q31 格式，等效浮点值约 0.02 和 0.00001。