频域维纳滤波降噪算法——设计思路与工程手记

写在前头

这份笔记记录的是一个面向实时语音通信的频域降噪算法。它的设计目标很明确：在 48kHz 采样率下，以 10ms 帧长实时处理单通道语音，把稳态和准稳态背景噪声压下去，同时尽量不伤语音。算法脱胎于 WebRTC 的 Noise Suppression 模块，但经过了面向嵌入式的重构——所有动态分配被消除，所有 C++ 模板和虚函数被拍平成纯 C 结构体。

我按照信号流经的顺序来谈。先说为什么要这么做，再说具体怎么做。

1. 信号模型

一切降噪算法的出发点都是同一个假设：麦克风采集到的信号 $x(n)$ 是干净语音 $s(n)$ 和加性噪声 $d(n)$ 的叠加：

$$ x(n) = s(n) + d(n) $$

这个模型简单到近乎天真，但它成立的前提条件其实很苛刻：噪声和语音不相关，噪声在短时间尺度上是准平稳的。实际工程中，空调嗡嗡声、风扇声、马路上的低频轰鸣，在几十毫秒内确实变化不大——这就够了。

我们的目标是从 $x(n)$ 中恢复 $\hat{s}(n)$。在频域里，这意味着对每个频率 bin 施加一个增益 $G(k)$：

$$ \hat{S}(k) = G(k) \cdot X(k) $$

问题归结为：怎么求这个 $G(k)$？

2. 为什么选频域

时域降噪（比如 LMS 自适应滤波）需要参考信号，这在单麦克风场景下不现实。频域方法的好处是：

  1. 频率分辨率。噪声和语音在不同频段的能量分布差异很大。白噪声全频段平坦，语音能量集中在 300Hz–4kHz。频域处理可以逐 bin 调节衰减深度。

  2. 短时平稳假设更容易满足。一帧 10ms 内，噪声功率谱几乎不变，这给估计器提供了稳定的观测窗口。

  3. 计算量可控。256 点 FFT 在嵌入式平台上是成熟操作，远比时域长 FIR 滤波器经济。

代价是引入了帧延迟和频谱泄漏，这些后面会谈。

3. 总体处理流程

一帧 480 个采样点（48kHz × 10ms）进来后，经过以下环节：

Fig001 总体处理流程

图 Fig001：算法总体处理流程。左侧为 Band 0 频域主链路，右侧为高频段时域增益处理。

核心思想：只在最低子带做完整的频域处理，高频段只施加一个标量增益。这是工程上的折中——人耳对高频噪声的敏感度低于低频，而高频做逐 bin 处理既费算力又容易引入音乐噪声。

4. 三分频器：多相滤波 + DCT 调制

4.1 为什么要分频

直接对 48kHz 信号做 FFT 当然可以，但有两个问题：

  • 256 点 FFT 在 48kHz 下频率分辨率只有 187.5Hz/bin，低频段太粗。
  • 语音能量集中在 8kHz 以下，高频段几乎是噪声的天下。分频后对低频段做精细处理，高频段粗放处理，算力用在刀刃上。

分频后每个子带等效采样率为 16kHz，256 点 FFT 的频率分辨率变为 62.5Hz/bin——对语音基频（通常 80–400Hz）来说足够细了。

4.2 多相分解

三分频器的原型是一个低通 FIR 滤波器 $H(z)$，通带截止 7kHz，阻带起始 9.2kHz，阻带衰减 40dB。利用 Noble Identity 做多相分解：

$$ H(z) = H_0(z^3) + z^{-1} H_1(z^3) + z^{-2} H_2(z^3) $$

每个子滤波器只有 4 个系数（原型滤波器长度 = 4 × 3 × 子带数 = 48 阶，但多相分解后每相只有 4 tap）。实际实现中，由于 DCT 调制的周期性，12 个多相分支里有 2 个恒为零，只剩 10 个非零分支。

4.3 DCT 调制

三个子带的中心频率分别是 $f_s/12$、$3f_s/12$、$5f_s/12$（即 4kHz、12kHz、20kHz @ 48kHz）。调制矩阵是一个 3 点 DCT：

$$ \begin{bmatrix} 2 & 2 & 2 \\ \sqrt{3} & 0 & -\sqrt{3} \\ 1 & -2 & 1 \end{bmatrix} $$

分析时：先抽取（3:1 降采样），再对每个多相分支做 4-tap 滤波，最后乘以 DCT 系数并累加得到三个子带。综合是逆过程。

4.4 工程代价

这个滤波器组不满足完美重建条件。直接分析再综合（中间不做任何处理），SNR 大约只有 9.5dB。这是有意为之——用重建精度换滤波器复杂度和延迟。对降噪应用来说，这个损失可以接受，因为我们的目标本来就不是无损还原。

滤波器组引入 24 个采样点的群延迟（0.5ms @ 48kHz），高频段需要额外延迟补偿来对齐。

5. 帧分割与加窗

5.1 帧参数

分频后的子带信号（160 samples/帧）不能直接送 FFT，需要：

  • 扩展帧：将当前 160 点与上一帧尾部 96 点拼接，形成 256 点的扩展帧。
  • 加窗：对扩展帧施加一个混合窗——前 96 点是升 Hanning，中间 65 点是平坦段（增益为 1），后 95 点是降 Hanning。

这个窗的设计意图：中间平坦段保证帧中心区域的信号不被衰减，两端的 Hanning 渐变抑制频谱泄漏。数学上，窗函数 $w(n)$ 满足：

$$ w(n) = \begin{cases} \sin\left(\frac{\pi}{2} \cdot \frac{n}{96}\right), & 0 \le n < 96 \\ 1, & 96 \le n < 161 \\ \sin\left(\frac{\pi}{2} \cdot \frac{256-n}{96}\right), & 161 \le n < 256 \end{cases} $$

5.2 Overlap-Add

IFFT 之后得到 256 点时域信号，再乘一次同样的窗（分析窗 × 综合窗 = COLA 条件），然后取前 160 点与上一帧残余 96 点叠加：

$$ y(n) = \underbrace{w(n) \cdot \text{IFFT}\{G(k) \cdot \text{FFT}\{w(n) \cdot x_{ext}(n)\}\}}_{\text{当前帧贡献}} + \underbrace{m(n)}_{\text{上帧残余}} $$

其中 $m(n)$ 是上一帧 256 点输出的后 96 点。这保证了帧间连续性，避免块效应。

6. 噪声功率谱估计

这是整个算法里最关键的环节。增益算得准不准，全看噪声估计跟不跟得上。

6.1 分位数噪声跟踪

核心思想：对每个频率 bin，在 log 域跟踪其功率谱的 25% 分位数。为什么是 25%？因为语音出现时功率会跳高，而噪声主导时功率处于低分位。25% 分位数大致对应"没有语音时的噪声底"。

具体做法：维护 $S = 6$ 组并行估计器，每组独立更新但周期错开。每组有一个 log 域分位数估计 $\hat{q}_s(k)$ 和一个概率密度估计 $\hat{f}_s(k)$。每来一帧：

$$ \hat{q}_s(k) \leftarrow \hat{q}_s(k) + \begin{cases} +\alpha \cdot \frac{\delta}{N_s+1}, & \log|X(k)|^2 > \hat{q}_s(k) \\ -(1-\alpha) \cdot \frac{\delta}{N_s+1}, & \log|X(k)|^2 \le \hat{q}_s(k) \end{cases} $$

其中 $\alpha = 0.25$ 是目标分位数，$\delta = 40 / \hat{f}_s(k)$ 是自适应步长（密度越大步长越小，收敛越稳），$N_s$ 是该组已累积的帧数。

6 组估计器以不同周期重置（200 帧为周期，各组偏移 33 帧），保证任何时刻都有一组处于"已收敛"状态。这个设计避免了单组估计器重置时的跳变。

6.2 启动阶段的参数化噪声模型

前 50 帧（500ms），分位数估计器还没收敛。这时用一个参数化模型来填充：

白噪声模型：全频段平坦，电平为所有 bin 的均值乘以过减因子。

粉噪声模型：功率谱随频率下降，建模为 $P(k) = A / k^{\beta}$。在 log 域做线性回归：

$$ \log P(k) = \log A - \beta \cdot \log k $$

用最小二乘拟合 $\log A$（截距）和 $\beta$（斜率），约束 $\beta \in [0, 1]$。

启动阶段，最终噪声估计是分位数估计和参数化模型的加权混合：

$$ \hat{N}(k) = \frac{n}{50} \cdot \hat{N}_{quantile}(k) + \frac{50-n}{50} \cdot \hat{N}_{parametric}(k) $$

其中 $n$ 是当前帧号。随着帧数增加，参数化模型逐渐退出。

6.3 后验更新：语音概率加权的噪声跟踪

分位数估计给出的是"保守下界"。真正的噪声谱还需要一阶 IIR 平滑来跟踪缓慢变化：

$$ \hat{N}(k) \leftarrow \gamma \cdot \hat{N}_{prev}(k) + (1-\gamma) \cdot \big[ P_{ns}(k) \cdot |X(k)| + P_s(k) \cdot \hat{N}_{prev}(k) \big] $$

其中 $P_s(k)$ 是语音存在概率，$P_{ns}(k) = 1 - P_s(k)$。

关键设计：当 $P_s(k) > 0.2$ 时，$\gamma$ 从 0.9 跳到 0.99——几乎冻结更新。这防止了语音段把噪声估计拉高。只有当语音概率很低时，才允许噪声估计快速跟踪。

另外维护一个"保守噪声谱"，只在 $P_s < 0.2$ 时以 0.05 的步长更新。这个保守谱后面会用于语音概率计算中的模板匹配。

7. 语音存在概率

7.1 为什么需要它

维纳滤波的增益公式里，先验 SNR 的估计依赖于"当前帧到底有没有语音"。如果把噪声段误判为语音，噪声估计会被污染；如果把语音段误判为噪声，语音会被削掉。所以需要一个逐帧、逐 bin 的语音存在概率。

7.2 三个特征

算法提取三个全局特征（每帧一个标量，不是逐 bin 的）：

（a）对数似然比 (LRT)

对每个 bin 计算：

$$ \ell(k) = (\gamma_{post}(k) + 1) \cdot \frac{2\gamma_{prior}(k)}{1 + 2\gamma_{prior}(k)} - \log(1 + 2\gamma_{prior}(k)) $$

然后做时间平滑：$\bar{\ell}(k) \leftarrow 0.5 \cdot \bar{\ell}(k) + 0.5 \cdot \ell(k)$

全频段均值 $L = \frac{1}{129}\sum_k \bar{\ell}(k)$ 作为 LRT 特征。

直觉：如果信号里只有噪声，后验 SNR 在 0 附近波动，$\ell(k)$ 趋近于零。有语音时，$\ell(k)$ 显著为正。

（b）频谱平坦度 (Spectral Flatness)

$$ F = \frac{\exp\left(\frac{1}{129}\sum_k \log |X(k)|\right)}{\frac{1}{129}\sum_k |X(k)|} $$

就是几何均值除以算术均值。白噪声的频谱平坦，$F \to 1$；语音有谐波结构，$F$ 明显小于 1。

时间平滑：$F \leftarrow 0.7 F + 0.3 F_{new}$

（c）频谱差异 (Spectral Difference)

当前频谱与保守噪声模板之间的归一化方差差：

$$ D = \frac{\text{Var}(X) - \text{Cov}(X, N)^2 / \text{Var}(N)}{\text{Energy}_{avg}} $$

这本质上是一个偏相关系数——去掉噪声模板能解释的部分后，剩余方差有多大。剩余越大，说明当前信号越不像纯噪声。

7.3 直方图建模与先验模型

三个特征不是直接拿来判决的。算法维护三个 100-bin 直方图，在 500 帧的窗口内累积特征的分布。每 500 帧提取一次模型参数：

  • LRT 直方图 → 确定 LRT 阈值 $\theta_{lrt}$（低波动时设为 1.0，否则取均值的 1.2 倍）
  • 平坦度直方图 → 找主峰位置，乘以 0.9 作为平坦度阈值 $\theta_{flat}$
  • 频谱差异直方图 → 找主峰位置，乘以 1.2 作为差异阈值 $\theta_{diff}$

同时确定三个特征的权重（被拒绝的特征权重置零，剩余特征均分权重）。

7.4 先验概率的 Sigmoid 融合

三个特征各自通过一个 tanh 映射变成 [0, 1] 的指示函数：

$$ I_0 = \frac{1}{2}\left[\tanh\big(w_0 \cdot (L - \theta_{lrt})\big) + 1\right] $$
$$ I_1 = \frac{1}{2}\left[\tanh\big(w_1 \cdot (\theta_{flat} - F)\big) + 1\right] $$
$$ I_2 = \frac{1}{2}\left[\tanh\big(w_2 \cdot (D - \theta_{diff})\big) + 1\right] $$

注意 $I_1$ 里是 $\theta_{flat} - F$（反向），因为平坦度越低越像语音。

加权求和后做一阶平滑得到全局先验语音概率：

$$ P_{prior} \leftarrow 0.9 \cdot P_{prior} + 0.1 \cdot (w_0 I_0 + w_1 I_1 + w_2 I_2) $$

7.5 逐 bin 后验概率

最终，将全局先验和逐 bin 的似然比结合，得到每个 bin 的语音存在概率：

$$ P_s(k) = \frac{1}{1 + \frac{1 - P_{prior}}{P_{prior}} \cdot e^{-\bar{\ell}(k)}} $$

这就是贝叶斯后验：先验 odds 乘以似然比。当 $\bar{\ell}(k)$ 很大时（强语音），$P_s(k) \to 1$；当 $\bar{\ell}(k) \approx 0$ 时（纯噪声），$P_s(k) \to P_{prior}$。

8. 维纳滤波器

8.1 增益公式

有了噪声谱 $\hat{N}(k)$ 和先验 SNR $\gamma_{prior}(k)$，维纳增益为：

$$ G(k) = \frac{\gamma_{prior}(k)}{\alpha + \gamma_{prior}(k)} $$

其中 $\alpha$ 是过减因子（over-subtraction factor）。当 $\alpha = 1$ 时是标准维纳滤波；$\alpha > 1$ 时更激进，噪声减得更深但语音失真风险增大。

增益有下限 $G_{min}$（称为 minimum attenuating gain），防止把信号完全静音：

$$ G(k) = \max\left(\frac{\gamma_{prior}(k)}{\alpha + \gamma_{prior}(k)},\; G_{min}\right) $$

8.2 先验 SNR 的定向判决估计

先验 SNR 不是直接算的（因为你不知道干净语音是什么），而是用 Ephraim & Malah (1985) 的 Decision-Directed 方法：

$$ \gamma_{prior}(k) = 0.98 \cdot \underbrace{\frac{|X_{prev}(k)|^2}{\hat{N}_{prev}(k)} \cdot G_{prev}(k)}_{\text{上一帧的"干净"估计}} + 0.02 \cdot \underbrace{\max\left(\frac{|X(k)|^2}{\hat{N}(k)} - 1,\; 0\right)}_{\text{当前帧后验 SNR}} $$

第一项是"上一帧用增益滤波后的信号功率 / 上一帧噪声"——代表对上一帧语音功率的估计。第二项是当前后验 SNR。

为什么 0.98/0.02 这个比例？因为后验 SNR 波动极大（噪声的随机性），如果权重给多了，增益会逐帧剧烈抖动，听起来就是"音乐噪声"。0.98 的时间常数约 50 帧（500ms），足够平滑。代价是对突发语音的响应慢了约 1-2 帧——可以接受。

8.3 启动阶段的滤波器混合

前 50 帧，维纳滤波器也在收敛。这时额外计算一个基于累积频谱的初始滤波器：

$$ G_{init}(k) = \frac{\sum_{n=0}^{N} |X_n(k)| - \alpha \cdot \hat{N}_{parametric}(k)}{\sum_{n=0}^{N} |X_n(k)|} $$

最终滤波器是两者的线性混合：

$$ G(k) = \frac{n}{50} \cdot G_{DD}(k) + \frac{50-n}{50} \cdot G_{init}(k) $$

9. 增益修正：让降噪"听起来自然"

9.1 问题

纯维纳滤波有一个听感问题：在语音段，增益接近 1，没问题；但在语音和噪声的过渡区域，增益变化太快，会产生"呼吸效应"——噪声突然消失又突然出现。

9.2 解决方案

算法计算滤波前后的能量比：

$$ g = \sqrt{\frac{E_{after}}{E_{before}}} $$

然后根据 $g$ 和全局语音概率 $P_{prior}$ 做一个双向缩放：

  • 当 $g > 0.5$（滤波前后能量差别不大，说明是语音段）：$S_1 = 1 + 1.3(g - 0.5)$，轻微放大，补偿维纳滤波对语音的少量衰减。

  • 当 $g < 0.5$（能量被大幅削减，说明是噪声段）：$S_2 = 1 - 0.3(0.5 - \max(g, G_{min}))$，轻微缩小，但不过度——避免静音段出现"抽吸感"。

最终缩放因子：

$$ S = P_{prior} \cdot S_1 + (1 - P_{prior}) \cdot S_2 $$

用语音概率做平滑过渡，而不是硬切换。

10. 高频段处理

Band 1 和 Band 2 不做 FFT，只施加一个标量增益。这个增益从 Band 0 的信息推导：

  1. 取 Band 0 最高 32 个 bin 的平均语音概率 $\bar{P}_s$ 和平均滤波器增益 $\bar{G}$。

  2. 如果中间有 AEC 等模块削了信号，用分析/处理频谱比修正 $\bar{P}_s$。

  3. 通过 tanh 映射得到基础增益：$g_{base} = \frac{1}{2}\big[1 + \tanh(2\bar{P}_s - 1)\big]$

  4. 与低频增益混合：

    • 若 $\bar{P}_s \ge 0.5$（语音为主）：$g = 0.25 g_{base} + 0.75 \bar{G}$
    • 若 $\bar{P}_s < 0.5$（噪声为主）：$g = 0.5 g_{base} + 0.5 \bar{G}$

语音段更信任低频滤波器的逐 bin 结果（权重 0.75），噪声段更信任概率判决（各 0.5）。

高频段还施加了一个延迟补偿（24 samples），对齐三分频器的群延迟。

11. 启动策略

算法有两个启动阶段：

阶段帧数时长行为
短启动0–500–500ms参数化噪声模型渐退，维纳滤波器混合
长启动0–2000–2s分位数估计器收敛，信号能量归一化建立
稳态>200>2s所有模块独立运行，增益修正启用

前 200 帧内，增益修正（第 9 节）被禁用——因为能量统计还不稳定，贸然缩放会引入伪影。

零帧检测：如果输入全零（比如静音段），直接跳过分析。这防止了统计量被零信号拉偏——否则阈值会适应到"零就是噪声"，等真正的信号来了反而全被当语音放过。

12. 降噪档位

算法提供 4 个预设档位，通过两个参数控制降噪深度：

档位过减因子 $\alpha$增益下限 $G_{min}$最大衰减增益修正
01.00.56 dB关
11.00.2512 dB开
21.10.12518 dB开
31.250.0921 dB开

最大衰减 = $-20\log_{10}(G_{min})$ dB。

工程选择建议：嵌入式通话场景用档位 1 或 2。档位 3 在极低 SNR 下会引入明显的语音失真（“水下音”），除非噪声真的很大否则不推荐。

13. 关键设计决策回顾

为什么用 256 点 FFT 而不是 512 点？

160 点帧长（10ms @ 16kHz 等效）加 96 点 overlap 刚好凑 256。512 点会提高频率分辨率但增加算法延迟到 20ms+，对实时通话不可接受。256 点在 62.5Hz 分辨率下已经能区分基频谐波。

为什么噪声估计用 6 组并行？

单组分位数估计器在重置时会跳变。6 组错开重置周期，任何时刻取其中一组的结果，保证输出连续。这是用 6 倍存储换平滑性——在 RAM 紧张的嵌入式上，这是 12KB 的代价（每组 129 bin × 2 个 float 数组 × 6 组）。

为什么高频段不做频域处理？

三个原因：(1) 算力——每多一个频段做 FFT/IFFT，MIPS 翻倍；(2) 高频段语音能量低，逐 bin 增益的估计误差大，容易产生音乐噪声；(3) 人耳对高频噪声的掩蔽效应强，粗放处理听感差异不大。

为什么 Decision-Directed 的平滑因子是 0.98？

这是 Ephraim-Malah 论文里的经典值。0.98 对应约 500ms 时间常数。太小（如 0.9）→ 增益波动大，音乐噪声明显；太大（如 0.999）→ 对语音起始响应太慢，会吃掉辅音。0.98 是大量主观听音实验的折中。

14. 数值实现备忘

  • 幅度谱计算加了 +1 的偏置（$|X(k)| + 1$），防止后续除法出现零除和 log(0)。
  • 所有 log/exp 用快速近似（查表 + 线性插值），不用标准库。在 ARM Cortex-M4 上，logf 要 80+ cycles，查表只要 5 cycles。
  • SNR 计算中的分母加了 0.0001 的 epsilon，这是 Q15 定点时代的遗留习惯，浮点下其实可以更小。
  • 增益下限 $G_{min}$ 同时充当"噪声门"的角色——即使算出来的增益更低，也不会让信号低于这个电平。这避免了完全静音段出现"数字死寂"的不自然感。

15. 算法局限性

坦率地说，这个算法有几个已知的短板：

  1. 非稳态噪声。突然的关门声、键盘敲击，分位数跟踪器来不及反应，会有 200-500ms 的泄漏窗口。
  2. 多说话人。算法假设只有一个目标说话人，第二个人的声音会被部分当作噪声处理。
  3. 音乐噪声。尽管 DD 方法和增益下限都在抑制，极低 SNR（< 0dB）下仍然能听到零星的"鸟鸣声"。
  4. 滤波器组不完美重建。即使不降噪，分析-综合链路也有约 9.5dB 的 SNR 损失。对高保真应用不可接受，对通话场景勉强够用。
  5. 启动时间。完整收敛需要 2 秒。如果设备频繁开关，前 2 秒的降噪效果是不稳定的。

这些是经典频域方法的固有局限。要根本性解决，得上深度学习方案（RNNoise、DTLN 之类），但那就是另一个量级的算力和内存需求了。

笔记完。整理于嵌入式降噪移植项目。