speex
2026-04-17 | 分类:未分类 | 评论:0人 | 浏览:11次
- 正文内容
- 我来说两句:(已有0人参与）
# 1 算法介绍
该代码基于经典的 **MDF（多延迟块频域自适应滤波）** 算法，并采用了 **AUMDF（交替更新MDF）** 变体，以及 **无显式双讲检测** 的变步长机制。代码的核心思想是：**在频域进行快速卷积和滤波器更新，同时通过“前景/背景”双滤波器机制和动态功率控制来保证鲁棒性。**
此实现参考了 Speex 和 Opus 中的 AEC 代码，是实际工程中常用的频域回声消除方案。该算法在VoIP、会议系统等领域得到广泛应用，是频域自适应滤波的经典实现。
本代码实现的是Speex音频编解码器中的回声消除算法，基于以下论文：
1. **MDF算法**：J. S. Soo, K. K. Pang, “Multidelay block frequency adaptive filter,” IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-38, No. 2, February 1990.
2. **变学习率**：J.-M. Valin, “On Adjusting the Learning Rate in Frequency Domain Echo Cancellation With Double-Talk,” IEEE TASLP, Vol. 15, No. 3, pp. 1030-1034, 2007.
**代码亮点解读**
1. **定点模拟**：`float_to_short` 并不是为了节省内存，而是为了复现原版 C 语言（Speex）在 DSP 上运行时因为截断误差带来的特性。
2. **半步更新**：文件头注释中提到 “we only move half the way toward the goal”，体现在 `ss = .35/M`，这是一种动量或软约束技巧，能极大提升频域更新的稳定性。
3. **彻底抛弃了 VAD/GDD**：传统的 AEC 需要一个严格的“双讲检测器”来决定是否冻结滤波器，一旦误判就会要么有回声，要么把近端人声抹掉。本代码完全通过能量的连续统计（`Davg`, `Dvar`）和前景背景切换机制，实现了“润物细无声”的非线性自适应，这是 Valin 2007 论文的最大贡献。
Speex MDF算法代码实现了一个基于 **多延迟块频域（MDF）** 的回声消除器，并采用了 **交替更新 MDF（AUMDF）** 变体和 **可变学习率** 以提高双讲鲁棒性。
是一个高度优化的频域自适应回声消除算法，其核心创新包括：
**算法特点**
1. **MDF/AUMDF结构**：高效处理长回声路径
2. **双滤波器机制**：保证输出稳定性
3. **连续变学习率**：无需显式双讲检测，实现鲁棒的双讲处理
4. **泄漏估计**：精确估计残余回声，指导学习率调整
| 特性 | 描述 |
|——|——|
| **核心算法** | 多延迟块频域自适应滤波器 (MDF) |
| **变体** | 交替更新MDF (AUMDF) |
| **双讲处理** | 无显式检测，连续变学习率 |
| **结构** | 前景/背景双滤波器 |
| **应用** | 声学回声消除 (AEC) |
| 特性 | 优势 |
|——|——|
| **频域处理** | 计算效率高，$O(N \log N)$ vs $O(N^2)$ |
| **MDF结构** | 支持长回声路径，降低FFT大小 |
| **双滤波器** | 提供稳定输出，防止发散 |
| **变学习率** | 无需显式双讲检测，鲁棒性强 |
| **比例自适应** | 加速收敛，适应回声路径特性 |
– **多延迟块（MDF）**：将长滤波器分成多个短块，降低计算复杂度。
– **交替更新（AUMDF）**：每帧只更新部分块的时域约束，减少计算量。
– **前景/背景滤波器**：通过双讲检测，在双讲时暂停背景更新，防止发散；在安静时快速收敛。
– **连续可变学习率**：根据残余回声、双讲和背景噪声调整步长，无需显式双讲检测。
– **鲁棒性处理**：包含 DC 陷波、预加重、饱和检测、异常重置等机制。
## 1.1 关键公式汇总
**频域滤波**
$$\mathbf{Y}(k) = \sum_{j=1}^{M} \mathbf{X}_j(k) \odot \mathbf{W}_j(k)$$
**频域权值更新**
$$\mathbf{W}_j(k+1) = \mathbf{W}_j(k) + p_j \cdot \frac{\mathbf{E}(k) \odot \mathbf{X}_j^*(k)}{\mathbf{P}(k)}$$
**泄漏估计**
$$\hat{\rho} = \frac{P_{ey}}{P_{yy}} = \frac{\sum_i E_h(i) \cdot Y_h^*(i)}{\sum_i |Y_h(i)|^2}$$
**变学习率**
$$\mu_i = \frac{0.7 \cdot \hat{\rho} \cdot Y_i + 0.3 \cdot RER \cdot e_i}{(R_i + 1) \cdot P_i}$$
**滤波器切换条件**
$$\Delta \cdot |\Delta| > S_{ff} \cdot D_{bf}$$
其中 $\Delta = S_{ff} – S_{ee}$。
## 1.2 以下是完整的代码处理流程拆解：
### 第一阶段：初始化阶段__init__
当实例化 `MDF(fs, frame_size, filter_length)` 时，算法进行环境搭建：
1. **参数计算**：
– `M`：将长滤波器（如250ms）拆分成 `M` 个频域块（如按8ms一块）。
– `N`：FFT窗口大小（通常是 `frame_size` 的2倍）。
2. **内存分配**：为频域信号（`X`, `Y`, `E`）、滤波器权重（`W`, `foreground`）、功率谱（`power`）等分配复数/实数数组。
3. **窗函数与比例因子**：
– 生成汉宁窗（`window`），用于重叠保存法中平滑拼接。
– 初始化块比例因子（`prop`），按指数衰减分布，离当前帧越近的块权重越大。
4. **预处理参数**：
– 预加重系数（`preemph = 0.9`），用于高频提升。
– 根据采样率计算直流陷波滤波器（`notch_radius`）的参数。
### 第二阶段：外部帧处理循环 (`main_loop`)
这是与外部交互的接口，处理一整段语音：
1. **量化模拟**：将输入的浮点信号 `u`（远端/扬声器）和 `d`（近端/麦克风）转换为 `[-32768, 32767]` 的16位整数。这模拟了实际DSP芯片的定点运算行为。
2. **分帧送入**：将长信号按 `frame_size` 切片，逐帧送入核心处理函数 `speex_echo_cancellation_mdf`。
3. **结果重组**：将处理后的短整型结果转回浮点数，输出 `e`（残余误差/消除后的信号）和 `y`（估计的回声）。
### 第三阶段：核心帧处理 (`speex_echo_cancellation_mdf`)
这是算法的灵魂，每处理一帧（如8ms）执行一次。流程非常严密，可分为以下几个步骤：
#### 步骤 3.1：时域预处理
– **直流消除**：对麦克风信号 `mic` 通过二阶IIR陷波滤波器（`filter_dc_notch16`），滤除直流干扰。
– **预加重**：对麦克风和远端信号分别进行一阶高通滤波（减去上一帧乘以0.9的值），增强高频，使频谱更平坦。
#### 步骤 3.2：频域变换与排队 (FFT & Queue)
– 将远端信号的新一帧放入缓冲区 `self.x`，对整个 `N` 点窗口做 FFT，得到频域信号。
– **块移位**：将历史频域块 `self.X` 往后移动一位（最老的块被丢弃），将最新的FFT结果放入第0块。形成 $M+1$ 个频域块的队列。
#### 步骤 3.3：前景滤波与初始误差计算
*注：代码维护了两套滤波器权重——`foreground`（前景，用于输出）和 `W`（背景，用于摸索更新）。*
– 使用**前景滤波器** `foreground` 与所有的频域块 `X` 相乘并求和，得到估计回声 `Y`。
– 做 IFFT 还原到时域，利用重叠保存法，取出当前帧的有效部分。
– 计算初始误差：`e_foreground = 麦克风输入 – 前景估计回声`，并计算其能量 `Sff`。
#### 步骤 3.4：背景滤波器更新 (Weight Update – 核心)
如果麦克风没有发生饱和（削波）：
– **频域梯度下降**：结合当前误差 `E`、共轭远端信号、以及最重要的**频率相关功率归一化步长** `power_1`，更新背景权重 `W`。
*(这对应论文中不依赖双讲检测，而是用残余误差和背景噪声连续调节步长的机制)*
– **时域约束 (AUMDF的体现)**：为了防止频域卷积带来的循环卷积效应（块边界处的突变），将背景权重 `W` 做 IFFT 变回时域，**强行将后半部分置零**，然后再 FFT 变回频域。这确保了它等效于一个真实的线性卷积滤波器。且为了节省算力，每次只更新部分块（`2+self.cancel_count) % (M-1)`）。
#### 步骤 3.5：计算背景误差
– 使用刚刚更新过的**背景滤波器** `W` 再次计算估计回声 `y`。
– 计算背景误差能量 `See`（即使用新权重后的残余能量）。
#### 步骤 3.6：前景/背景切换逻辑 (Robustness 核心)
算法通过比较 `Sff`(旧前景误差) 和 `See`(新背景误差) 来决定下一步怎么走：
– **情况A：背景变好 (See < Sff)**：
认为“摸索”有了进步。将背景权重 `W` 提拔为前景权重 `foreground`。并对误差信号应用汉宁窗进行平滑重叠，避免切换带来的听觉突变。清空方差统计。
– **情况B：背景变差 (发散/双讲发生)**：
如果 `See` 远大于 `Sff`（通过方差判定 `Davg1/Davg2`），说明可能遇到了近端说话（双讲）或者突发噪声导致权重更新跑偏了。
**处理**：果断放弃背景权重，将背景 `W` 强行拉回当前的前景 `foreground`。保证输出质量绝不恶化。清空方差统计。
#### 步骤 3.7：后处理与输出
– **去预加重**：将误差信号加上高通滤掉的部分，恢复原始音色。
– **饱和检测**：如果发现原始麦克风信号达到满幅（如 > 32000），标记 `saturated`，在接下来的帧中暂停背景更新，防止误差过大污染权重。
– **异常重置**：如果出现负数能量等违背物理常理的情况（`screwed_up` 累加到50），直接调用 `__init__` 进行硬重启。
#### 步骤 3.8：功率谱与自适应步长更新 (无双讲检测的秘密)
这一步为下一帧的“步骤3.4”准备变量 `power_1`（学习率）：
1. 计算回声返回损耗增强（RER）。
2. **计算泄漏估计**：`leak_estimate`。通过比较残余误差功率和估计回声功率，动态评估当前还有多少比例的回声没消除干净。
3. **更新自适应步长 `power_1`**：
– 基础步长由远端信号能量 `Sxx` 和背景误差 `See` 决定。
– 如果远端声音大、误差小，步长就大（快速收敛）。
– 如果远端声音小或误差突然变大（暗示双讲或噪声），步长自动缩小（停止更新，保护权重）。
*(代码注释掉了一段直接用 `leak_estimate` 调节步长的高级逻辑，启用了其下方一段相对基础的逻辑，但核心思想一致)*
## 1.3 流程图总结
“`text
[远端信号 u] + [近端信号 d]
| |
V V
(预加重) (陷波滤波 + 预加重)
| |
V |
[FFT入队列 X] |
| |
+——+——+
|
V
+—> [前景滤波器 Foreground]
| |
| V
| 初始误差 E_foreground
| |
| V
| +—–+—–+
| | |
| V V
| [更新背景W] [计算背景误差 E_background]
| | |
| | V
| | [对比 E_foreground 和 E_background]
| | |
| | +—–+—–+
| | | |
| | [背景更好] [背景发散(双讲/噪声)]
| | | |
| | [W->前景] [前景->W (放弃更新)]
| | | |
+—–+—–+———–+
|
V
[去预加重 + 饱和保护] —> 输出: 消除回声后的信号 e
|
V
[更新功率谱 & 计算下一帧步长 power_1] —> (反馈给下一帧的 [更新背景W])
“`
# 2 预加重
在数字信号处理（尤其是语音和音频处理）中，**预加重** 是一个极其基础且重要的步骤。你之前阅读的 Speex AEC 代码中，那句 `tmp32 = self.input[i, chan] – (np.dot(self.preemph, self.memD[chan]))`（其中 `preemph=0.9`）就是在做这件事。
## 2.1 预加重
预加重本质上是一个**一阶高通滤波器**。
* **时域公式**：$y[n] = x[n] – \alpha \cdot x[n-1]$
* $x[n]$ 是当前的输入样本
* $x[n-1]$ 是上一个输入样本
* $\alpha$ 是预加重系数（语音处理中通常取 **0.9 到 0.97** 之间，代码里是 0.9）
* **频域特性（Z变换）**：$H(z) = 1 – \alpha z^{-1}$
* 它的频率响应大致以 **+6 dB/倍频程** 的斜率上升。
* **效果**：它会**放大高频信号，衰减低频信号**，使信号的频谱变得“平坦”一些。
人类的语音信号有一个天然的物理特性：**频谱倾斜**。由于声道的形状，语音信号在低频（如元音 a, e）的能量极高，而在高频（如摩擦音 s, sh, f）的能量衰减非常快（大约以 -12 dB/倍频程下降）。
进行预加重有三大核心目的：
1. 为了让 LPC（线性预测编码）更好地工作（最核心原因）
Speex 是基于 CELP（码激励线性预测）算法的。LPC 的核心思想是：**用过去的 $p$ 个样本预测当前的样本**。
* 如果不做预加重，信号的能量全集中在低频，LPC 滤波器的极点会全部挤在低频区域去拼命拟合那巨大的低频能量，导致**对高频的预测精度极差**。
* 经过预加重后，频谱变平坦了，LPC 滤波器就能“平均用力”，在整个频带上都能建立准确的预测模型。这对语音编码的音质、AEC 回声路径的估计准确性至关重要。
2. 均衡信噪比
在实际环境中（比如打电话），背景噪声（风扇声、电路底噪）的频谱通常是相对平坦的，而语音是低频强、高频弱。
如果不做预加重，在 4kHz 以上的高频段，**语音能量可能会比背景噪声还要低**，导致算法无法分辨高频语音和噪声。预加重把高频语音抬上去了，使得整个频带的信噪比更加均匀，有利于后续的噪声抑制（ANS）和回声消除。
3. 避免有限精度带来的数值问题
在定点运算（如 Speex 的 C 语言原生代码）中，低频的大能量数值容易导致整数溢出，而高频的小能量数值又会被截断为 0。预加重拉近了高低频的幅度量级，缓解了定点计算的动态范围压力。
## 2.2 逆预加重
既然把高频抬高了，算法处理完之后，必须把它还原回去，否则声音就不可用了。
* **时域公式**：$x[n] = y[n] + \alpha \cdot x[n-1]$
*(注意：这里用的是还原后的输出 $x[n-1]$ 进行反馈，而不是输入 $y[n-1]$)*
* **频域特性**：$H^{-1}(z) = \frac{1}{1 – \alpha z^{-1}}$
* 这是一个一阶低通滤波器，以 **-6 dB/倍频程** 衰减，刚好抵消预加重的影响。
为什么处理完还要做逆预加重？
在你之前看的 AEC 代码末尾，有这样一句：
`tmp_out = tmp_out + self.preemph * self.memE[chan]`
这就是在做逆预加重。原因很简单：
1. **恢复原始音色**：如果 AEC 输出的是预加重后的信号，听起来会非常尖锐、刺耳（全是被放大的高频气流声），必须通过逆预加重把高频压回去，还原成正常人的声音。
2. **对接下游模块**：AEC 只是前端模块，它的输出通常要送给 AGC（自动增益控制）或交给听筒播放。下游模块假设输入的是原始的、频谱倾斜的语音信号。如果不还原，下游模块（尤其是 AGC）会因为高频能量异常而做出错误的增益控制。
## 2.3 整个流水线的视角
在你的 AEC 代码中，流程是这样的：
1. **输入端**：原始麦克风信号 $\xrightarrow{\text{预加重}}$ 变平坦的信号。
2. **算法核心**：所有复杂的 FFT、滤波器更新（MDF）、回声估计，都在**平坦的频谱**上进行，获得了极高的估计精度。
3. **输出端**：估计出的干净误差信号 $\xrightarrow{\text{逆预加重}}$ 恢复原始频谱倾斜的信号 $\rightarrow$ 交给扬声器播放或后续处理。
这种**“前端抬高 -> 核心处理 -> 后端还原”**的套路，在音频 DSP 中无处不在（比如 MFCC 特征提取中的预加重也是这个道理）。
# 3 DC陷波滤波器
这个DC陷波滤波器的设计精巧地平衡了陷波深度、过渡带宽度和计算效率，是实时音频处理中的经典实现。
这是一个**二阶IIR直流陷波滤波器**，用于消除音频信号中的直流偏置(DC offset)分量。在自适应回声消除等应用中，直流分量会严重影响算法性能，因此需要预处理消除。
DC陷波滤波器 (DC Notch Filter) 原理与公式推导
| 特性 | 说明 |
|——|——|
| **滤波器类型** | 二阶IIR陷波滤波器 |
| **陷波位置** | DC (ω=0, z=1) |
| **零点** | z=1 双重零点 |
| **极点** | 共轭极点对，由r和den2决定 |
| **DC增益** | 0 (完美消除) |
| **高频增益** | ≈1 (透明通过) |
| **关键参数** | notch_radius (控制陷波宽度) |
| **计算复杂度** | ~9次浮点运算/样本 |
| **主要应用** | AEC预处理、音频预处理 |
## 3.1 数学模型推导
设：
– $x[n]$ = `vin` (输入信号)
– $y[n]$ = `out[i]` (输出信号)
– $v[n]$ = `vout` (中间变量)
– $m_0[n]$ = `mem[0]` (状态变量1)
– $m_1[n]$ = `mem[1]` (状态变量2)
– $r$ = `notch_radius` (陷波半径，通常接近1)
– $d$ = `den2` (分母系数)
## 3.2 差分方程提取
从代码逻辑可得状态更新方程：
$$
\begin{aligned}
v[n] &= m_0[n-1] + x[n] \\
m_0[n] &= m_1[n-1] + 2(-x[n] + r \cdot v[n]) \\
m_1[n] &= x[n] – d \cdot v[n] \\
y[n] &= r \cdot v[n]
\end{aligned}
$$
其中：
$$d = r^2 + 0.7(1-r)(1-r)$$
## 3.3 传递函数推导
将状态方程代入，消去中间变量：
**第一步**：将 $v[n] = m_0[n-1] + x[n]$ 代入 $m_0[n]$ 方程：
$$m_0[n] = m_1[n-1] + 2(-x[n] + r \cdot m_0[n-1] + r \cdot x[n])$$
简化：
$$m_0[n] = m_1[n-1] + 2(r-1)x[n] + 2r \cdot m_0[n-1]$$
**第二步**：对 $m_1[n]$ 方程进行代入：
$$m_1[n] = x[n] – d(m_0[n-1] + x[n]) = (1-d)x[n] – d \cdot m_0[n-1]$$
**第三步**：进行Z变换并求解传递函数。
设 $M_0(z)$, $M_1(z)$, $X(z)$, $V(z)$, $Y(z)$ 分别为对应信号的Z变换。
从中间变量 $v[n]$ 的方程：
$$V(z) = z^{-1}M_0(z) + X(z)$$
从 $m_0[n]$ 的Z变换：
$$M_0(z) = z^{-1}M_1(z) + 2(r-1)X(z) + 2r \cdot z^{-1}M_0(z)$$
整理：
$$M_0(z)(1 – 2r \cdot z^{-1}) = z^{-1}M_1(z) + 2(r-1)X(z)$$
从 $m_1[n]$ 的Z变换：
$$M_1(z) = (1-d)X(z) – d \cdot z^{-1}M_0(z)$$
**第四步**：联立求解传递函数 $H(z) = Y(z)/X(z)$。
经过代数运算（过程较繁琐），最终得到：
$$H(z) = \frac{r(1-z^{-1})^2}{1 – 2rz^{-1} + d \cdot z^{-2}}$$
### 3.3.1 传递函数分析
将传递函数写成标准形式：
$$H(z) = \frac{r(1-z^{-1})^2}{1 – 2rz^{-1} + d \cdot z^{-2}}$$
**零点分析**：
分子 $r(1-z^{-1})^2 = r \cdot z^{-2}(z-1)^2$
零点位于 **$z = 1$ (双重零点)**，即在DC频率 ($\omega = 0$) 处形成陷波。
**极点分析**：
分母 $1 – 2rz^{-1} + d \cdot z^{-2} = 0$
极点位置由 $r$ 和 $d$ 决定。由于 $r \approx 1$（典型值0.98-0.99），极点靠近单位圆，但不会导致不稳定。
## 3.4 频率响应特性
### 3.4.1 DC处增益
将 $z = e^{j\omega}$，当 $\omega = 0$ 时 $z = 1$：
$$H(1) = \frac{r(1-1)^2}{1 – 2r + d} = 0$$
**DC增益为零**，完美消除直流分量。
### 3.4.2 高频特性
当 $\omega \to \pi$ (高频)，$z \to -1$：
$$H(-1) = \frac{r(1+1)^2}{1 + 2r + d} = \frac{4r}{1 + 2r + d}$$
对于典型参数（$r = 0.98$, $d \approx 0.96$）：
$$H(-1) \approx \frac{3.92}{1 + 1.96 + 0.96} \approx 1$$
**高频增益接近1**，对音频信号影响极小。
### 3.4.3 过渡带宽度
陷波的宽度由参数 $r$ 控制：
| $r$ 值 | 陷波宽度 | 响应速度 |
|——–|———-|———-|
| 0.99 | 很窄 | 慢 |
| 0.98 | 窄 | 中等 |
| 0.95 | 较宽 | 快 |
$r$ 越接近1，陷波越窄，但收敛越慢（更多的历史样本参与）。
## 3.5 等效方框图
“`
┌─────────────────────────────────────┐
│ 状态空间实现 │
│ │
x[n] ──┬──►│ ┌───────┐ ┌───────┐ │
│ │ │ m₀[n] │ │ m₁[n] │ │
│ │ │ z⁻¹ │ │ z⁻¹ │ │
│ │ └───┬───┘ └───┬───┘ │
│ │ │ │ │
│ │ ┌───▼───┐ ┌───▼───┐ │
│ │ │ + │◄─────┤ + │ │
│ │ └───┬───┘ └───────┘ │
│ │ │ │
│ │ ┌───▼───┐ │
│ │ │ ×r │──────────────────────►│─── y[n]
│ │ └───────┘ │
│ │ │
│ │ 反馈路径：(复杂IIR结构) │
└──►│ · 2(r-1)x[n] 前馈 │
│ · 2r·m₀[n-1] 反馈 │
│ · (1-d)x[n] – d·m₀[n-1] │
└─────────────────────────────────────┘
“`
## 3.6 参数设计依据
den2系数的含义
“`python
den2 = r² + 0.7*(1-r)*(1-r)
“`
这个系数的设计使得极点位置满足：
– 极点实部 $\approx r$（接近单位圆）
– 极点虚部 $\approx \sqrt{0.3}(1-r)$（轻微离轴）
系数 **0.7** 是经验值，提供最佳的陷波宽度与稳定性的折中。
典型参数选择
“`python
# Speex/WebRTC 典型配置
notch_radius = 0.98 # 或 0.99
# 计算得到的 den2
# r = 0.98: den2 = 0.98² + 0.7×0.02² = 0.9604 + 0.00028 = 0.96068
# r = 0.99: den2 = 0.99² + 0.7×0.01² = 0.9801 + 0.00007 = 0.98017
“`
## 3.7 与标准DC Blocker的关系
### 3.7.1 一阶DC Blocker (简单形式)
标准一阶DC blocker：
$$H(z) = \frac{1-z^{-1}}{1-\alpha z^{-1}}$$
其中 $\alpha \approx 0.995$。
**特点**：
– 单零点在 z=1
– 单极点在 z=α
– 计算简单，但陷波不够陡峭
### 3.7.2 二阶DC Notch (本代码)
$$H(z) = \frac{r(1-z^{-1})^2}{1-2rz^{-1}+d \cdot z^{-2}}$$
**特点**：
– 双重零点在 z=1（陷波更深）
– 共轭极点对（过渡更平滑）
– 更好的低频响应
### 3.7.3 频率响应对比
| 特性 | 一阶DC Blocker | 二阶DC Notch |
|——|—————|————–|
| DC衰减 | ∞ (零点) | ∞ (双重零点) |
| 低频衰减斜率 | 20dB/dec | 40dB/dec |
| 计算复杂度 | 2次乘法/样本 | ~6次乘法/样本 |
| 过渡带 | 较宽 | 较窄 |
## 3.8 代码实现细节分析
“`python
def filter_dc_notch16(self, mic, mem):
out = np.zeros_like(mic)
# 预计算分母系数（只需计算一次）
den2 = self.notch_radius**2 + 0.7 * \
(1-self.notch_radius)*(1 – self.notch_radius)
# 逐样本处理（IIR滤波器需要顺序计算）
for i in range(self.frame_size):
vin = mic[i]
# 步骤1：计算中间输出 vout
# vout = mem[0] + vin（状态反馈 + 当前输入）
vout = mem[0] + vin
# 步骤2：更新状态 mem[0]
# mem[0] = mem[1] + 2*(-vin + r*vout)
# 展开：= mem[1] – 2*vin + 2*r*(mem[0] + vin)
# = mem[1] + 2*(r-1)*vin + 2*r*mem[0]
mem[0] = mem[1] + 2*(-vin + self.notch_radius*vout)
# 步骤3：更新状态 mem[1]
# mem[1] = vin – den2*vout
# 部分前馈 + 极点反馈
mem[1] = vin – (den2*vout)
# 步骤4：计算最终输出
# 输出 = r * vout（增益补偿）
out[i] = self.notch_radius * vout
return out, mem
“`
计算复杂度
每样本操作：
– 加法：5次
– 乘法：4次
– 总计：~9次浮点运算/样本
对于帧长N：
– 总运算量：O(9N) ≈ O(N)
状态变量初始化
“`python
mem = [0.0, 0.0] # 通常初始化为零
# 如果已知信号的初始DC偏置，可以预设状态以减少收敛时间
# mem[0] = -initial_dc_offset / notch_radius
“`
## 3.9 应用场景
### 3.9.1 自适应回声消除(AEC)
在AEC系统中，DC分量会导致：
1. 自适应滤波器权值的DC偏置累积
2. 双端检测(DTD)误判
3. 残留回声估计偏差
**解决方案**：在预加重之前先用DC notch消除直流。
### 3.9.2 音频编解码
DC偏置会：
– 降低编码效率（浪费比特）
– 产生低频噪声
– 影响静音检测
### 3.9.3 语音识别
DC偏置会改变能量分布，影响特征提取。
# 4 MDF自适应滤波器中各块的更新比例因子
MDF_adjust_prop 函数，这个函数用于**计算MDF（多延迟块频域）自适应滤波器中各块的更新比例因子**。在AUMDF（交替更新MDF）算法中，不同块的滤波器权值需要根据其”活跃程度”分配不同的学习率，这个比例因子决定了每个块应该获得多大的更新强度。
这个函数体现了**比例自适应**的核心思想：根据滤波器权值的当前状态动态调整学习率，使得算法能够快速收敛活跃区域，同时保持对变化区域的响应能力。
| 特性 | 说明 |
|——|——|
| **核心功能** | 计算MDF各块的更新比例因子 |
| **设计理念** | 能量大的块获得更大的学习率 |
| **关键公式** | $\text{prop}_i = \frac{0.99 \sqrt{1+E_i}}{1+\sum_j p_j}$ |
| **偏置作用** | 保证所有块都有最小更新量 |
| **归一化** | 比例因子之和 < 0.99 |
| **应用场景** | AUMDF频域权值更新 |
| **复杂度** | $O(MCKN)$ |
| 参数 | 含义 | 典型值 |
|——|——|——–|
| `N` | `window_size` | FFT窗口大小（如512, 1024） |
| `M` | 块数量 | 滤波器分块数（如4, 8, 16） |
| `C` | 输出通道数 | 麦克风数量 |
| `K` | 输入通道数 | 扬声器数量 |
| `prop` | 输出比例因子向量 | 维度 M×1 |
在 Speex 的 MDF（多延迟块频域，Multidelay Block Frequency Domain）自适应滤波器中，`MDF_adjust_prop` 函数（或其内部的步长调整逻辑）是整个回声消除算法的**“油门踏板”**。
它的核心设计原理可以用一句话概括：**根据当前的信号状态（远端能量、近端能量、残余误差），动态且平滑地调整滤波器权重更新的步长（比例因子），以在“收敛速度”和“稳态失调”之间找到最优平衡。**
以下是该函数设计原理的深度拆解：
1. 为什么不能使用固定步长？
在自适应滤波器（如 NLMS 或频域 LMS）中，更新公式本质上是：
$$ W_{new} = W_{old} + \mu \cdot (\text{误差信号}) \cdot (\text{远端信号特征}) $$
这里的 $\mu$ 就是步长（对应代码里的 prop 比例因子）。
* **如果 $\mu$ 太大（油门踩到底）：** 滤波器收敛极快，能迅速跟上回声路径的变化。但在稳态时，会引入极大的噪声，导致 ERLE 下降，甚至导致算法发散（声音炸裂）。
* **如果 $\mu$ 太小（油门松开）：** 算法非常稳定，残余噪声小。但一旦人动了一下（回声路径改变），或者刚开始打电话，算法需要极长的时间才能重新收敛，导致用户听到明显的回声。
2. `MDF_adjust_prop` 的四大设计原理
该函数通过复杂的逻辑来动态计算这个 $\mu$，主要基于以下四个原理：
原理一：基于远端信号功率的归一化
这是最基础的保护机制。远端声音大时，更新步长必须小；远端声音小时，步长可以适当放大。
* **逻辑：** 函数会持续估计远端信号 $X(k)$ 在各个频点的功率。步长会除以这个功率（加上一个极小的常数 $\epsilon$ 防止除零）。
* **目的：** 确保无论对方大声喊叫还是小声说话，注入到滤波器权重的“更新能量”是相对恒定的，防止大信号导致滤波器系数溢出。
原理二：基于残余误差与理论回声的比值（最核心的反馈机制）
这是 Speex MDF 算法的精妙之处。它会对比“实际误差”和“理论上应该剩下的误差”。
* **逻辑：** 算法内部有一个“理论残余回声估计”。如果发现**实际算出来的误差功率 > 理论残余回声功率**，说明什么？
* 可能性 A：近端开始说话了（双讲）。
* 可能性 B：回声路径突然发生了剧烈变化（比如拿开了手机）。
* **应对：** 在这种情况下，函数会**增大**比例因子 `prop`，试图快速追赶变化的回声路径（如果是双讲，后续的非线性处理 NLP 会负责切断，这里只管线性收敛）。
* 反之，如果实际误差 $\approx$ 理论误差，说明已经收敛得很好了，函数会**减小** `prop`，进入“精修”模式，降低稳态噪声。
原理三：回声路径变化检测的平滑过渡
当检测到回声路径可能发生变化时，不能瞬间把步长拉到最大，那会导致声音出现严重的突变（咔哒声）。
* **逻辑：** `MDF_adjust_prop` 内部使用了一阶 IIR 低通滤波器（指数平滑）来更新比例因子。
* **公式抽象：** $\mu_{当前} = \alpha \cdot \mu_{目标} + (1 – \alpha) \cdot \mu_{上一次}$
* **目的：** 让步长的增加是一个平滑的爬坡过程，保护听感。
原理四：频域差异化处理
因为 MDF 是在频域工作的，不同频率的特性完全不同。
* **逻辑：** `prop` 往往不是一个全局标量，而是一个频域向量（或者按频带划分）。低频能量大但人耳敏感，高频能量小但容易发散。
* **应对：** 函数在计算比例因子时，会结合各个频段的信噪比特性，对高频施加更严格的步长限制，对低频允许稍微大一点的步长以加快收敛。
## 4.1 代码逐行解析
“`python
def mdf_adjust_prop(self,):
N = self.window_size # FFT大小
M = self.M # 块数量
C = self.C # 输出通道数（麦克风）
K = self.K # 输入通道数（扬声器）
prop = np.zeros((M, 1),) # 初始化比例因子向量
# 遍历每个块
for i in range(M):
tmp = 1 # 初始化能量累加器（加1避免零值）
# 遍历所有输出通道（麦克风）
for chan in range(C):
# 遍历所有输入通道（扬声器）
for speak in range(K):
# 累加当前块所有频点的权值模平方
# W 的维度: [频点 × 扬声器 × 块索引 × 麦克风]
tmp = tmp + np.sum(np.abs(self.W[:N//2+1, speak, i, chan])**2)
# 取平方根得到”权值幅度”
prop[i] = np.sqrt(tmp)
# 找到最大值（但至少为1，避免除零）
max_sum = np.maximum(prop, 1)
# 加入偏置项，避免某些块比例过小
prop = prop + .1 * max_sum
# 计算总和（加1避免零值问题）
prop_sum = 1 + np.sum(prop)
# 归一化并缩放到0.99
prop = 0.99 * prop / prop_sum
return prop
“`
## 4.2 数学公式推导
### 4.2.1 权值能量计算
对于第 $i$ 个块，计算其权值总能量：
$$E_i = \sum_{chan=0}^{C-1} \sum_{speak=0}^{K-1} \sum_{k=0}^{N/2} |W_k^{(speak, i, chan)}|^2$$
其中：
– $k$ 是频点索引（只取正频率部分：$0$ 到 $N/2$）
– $W_k^{(speak, i, chan)}$ 是频域权值
### 4.2.2 比例因子计算
**第一步**：计算能量平方根
$$p_i = \sqrt{1 + E_i}$$
**第二步**：加入偏置
$$p_i \leftarrow p_i + 0.1 \cdot \max_j(p_j, 1)$$
**第三步**：归一化
$$\text{prop}_i = \frac{0.99 \cdot p_i}{1 + \sum_{j=0}^{M-1} p_j}$$
### 4.2.3 比例因子性质
归一化后满足：
$$\sum_{i=0}^{M-1} \text{prop}_i = \frac{0.99 \sum_i p_i}{1 + \sum_i p_i} < 0.99$$
这保证了所有块的比例因子之和略小于1，留有一定的”余地”（稳定性保证）。
## 4.3 设计原理
为什么需要比例因子？
在MDF算法中，长滤波器被分解为多个短块：
“`
时域滤波器: h[0], h[1], …, h[L-1]
↓ 分块
块0: h[0], …, h[N-1] → 频域权值 W_0
块1: h[N], …, h[2N-1] → 频域权值 W_1
…
块M-1: h[(M-1)N], …, h[MN-1] → 频域权值 W_{M-1}
“`
**核心问题**：不同块的权值对滤波器输出的贡献不同。
在回声消除场景中：
– **近端块（索引小）**：对应较小的延迟，通常包含更多的回声能量
– **远端块（索引大）**：对应较大的延迟，可能只包含噪声或很弱的回声
因此，应该根据每个块的”重要性”分配不同的更新强度。
### 4.3.1 能量与重要性的关系
权值能量大的块意味着：
1. 该延迟区间有较强的回声路径响应
2. 该块的调整对输出影响较大
3. 应该获得更大的学习率
### 4.3.2 偏置项的作用
“`python
prop = prop + .1 * max_sum
“`
加入偏置项的目的：
| 无偏置 | 有偏置 |
|——–|——–|
| 能量为0的块比例≈0 | 最小比例≈0.1/M |
| 可能导致某些块”死亡” | 保证所有块都有最小更新 |
| 对新出现的回声响应慢 | 能快速适应变化 |
## 4.4 在AUMDF更新中的应用
权值更新公式，在MDF的频域权值更新中：
$$W_i^{(new)} = W_i^{(old)} + \mu \cdot \text{prop}_i \cdot \frac{X_i^* \odot E}{\|X_i\|^2 + \varepsilon}$$
其中：
– $\mu$ 是全局学习率
– $\text{prop}_i$ 是该块的比例因子
– $X_i$ 是该块对应的频域输入
– $E$ 是频域误差
– $\odot$ 表示元素乘法
比例因子的影响
| 块特性 | 能量 | prop_i | 更新强度 |
|——–|——|——–|———-|
| 活跃块 | 大 | 大 | 强 |
| 不活跃块 | 小 | 小 | 弱 |
| 静默块 | 0 | 最小值 | 最小 |
## 4.5 示例计算
假设 $M=4$ 块，计算得到各块能量：
“`
块0: E=1000 → p=√1001 ≈ 31.6
块1: E=500 → p=√501 ≈ 22.4
块2: E=100 → p=√101 ≈ 10.0
块3: E=10 → p=√11 ≈ 3.3
“`
**加入偏置** ($\max = 31.6$)：
“`
p’ = [31.6 + 3.16, 22.4 + 3.16, 10.0 + 3.16, 3.3 + 3.16]
= [34.76, 25.56, 13.16, 6.46]
“`
**归一化**：
“`
sum = 1 + 34.76 + 25.56 + 13.16 + 6.46 = 80.94
prop = 0.99 × [34.76, 25.56, 13.16, 6.46] / 80.94
= [0.425, 0.312, 0.161, 0.079]
“`
**结果解读**：
| 块 | 能量占比 | prop占比 | 解读 |
|—-|———-|———-|——|
| 0 | 61.7% | 42.5% | 最活跃，获得最大更新 |
| 1 | 30.9% | 31.2% | 较活跃 |
| 2 | 6.2% | 16.1% | 较弱，偏置提升明显 |
| 3 | 0.6% | 7.9% | 最弱，偏置保证最小更新 |
## 4.6 频点范围说明
“`python
self.W[:N//2+1, speak, i, chan]
“`
只取前 `N//2+1` 个频点的原因：
| 频点范围 | 数量 | 说明 |
|———-|——|——|
| 0 | 1 | DC分量 |
| 1 到 N/2-1 | N/2-1 | 正频率 |
| N/2 | 1 | Nyquist频率 |
| N/2+1 到 N-1 | N/2-1 | 负频率（共轭对称） |
对于实数信号，频域权值具有共轭对称性：
$$W_k = W_{N-k}^*$$
因此只需计算正频率部分的能量，避免重复计算。
## 4.7 与其他技术的关联
### 4.7.1 与分块NLMS的对比
| 方法 | 比例分配策略 |
|——|————-|
| 标准NLMS | 所有系数相同学习率 |
| PNLMS（比例NLMS） | 根据每个系数的幅度分配 |
| **MDF adjust_prop** | 根据**每个块**的总能量分配 |
### 4.7.2 与PNLMS的区别
PNLMS (Proportionate NLMS) 为**每个系数**设置独立的学习率：
$$\mu_k \propto |w_k| + \delta$$
而 `mdf_adjust_prop` 为**每个块**设置统一的学习率：
$$\mu_i \propto \sqrt{\sum_k |W_k^{(i)}|^2}$$
**MDF采用块级比例的原因**：
1. 减少计算量（M个块 vs MN个系数）
2. 频域中系数之间有耦合
3. 块级比例更稳定
## 4.8 性能优化考虑
### 4.8.1 当前实现的复杂度
“`python
for i in range(M): # M次迭代
for chan in range(C): # C次迭代
for speak in range(K): # K次迭代
np.sum(…) # O(N) 操作
“`
总复杂度：**$O(M \cdot C \cdot K \cdot N)$**
### 4.8.2 可能的优化
“`python
# 向量化实现
def mdf_adjust_prop_vectorized(self):
# 一次性计算所有块的能量
# W形状: [N//2+1, K, M, C]
W_energy = np.sum(np.abs(self.W[:self.N//2+1, :, :, :])**2, axis=(0, 3))
# 现在形状: [K, M]
prop = np.sqrt(1 + np.sum(W_energy, axis=0)) # 形状: [M]
max_sum = np.maximum(prop, 1)
prop = prop + 0.1 * max_sum
prop_sum = 1 + np.sum(prop)
prop = 0.99 * prop / prop_sum
return prop.reshape(-1, 1)
“`
优化后复杂度：**$O(M \cdot C \cdot K \cdot N)$ 但常数更小**
# 5 MDF算法原理
## 5.1 时域自适应滤波回顾
**标准LMS算法**（时域）：
$$\mathbf{w}(n+1) = \mathbf{w}(n) + \mu \cdot e(n) \cdot \mathbf{x}(n)$$
**问题**：
– 长回声路径需要长滤波器（如250ms @ 16kHz = 4000抽头）
– 时域卷积计算复杂度 $O(L^2)$
– 每个样本都要更新，计算量大
## 5.2 频域自适应滤波 (FDAF)
**核心思想**：利用FFT将时域卷积转换为频域点乘。
**重叠保存法 (Overlap-Save)**：
设滤波器长度为 $L$，帧长为 $N = 2L$：
**时域卷积**：
$$y(n) = \sum_{k=0}^{L-1} w(k) \cdot x(n-k)$$
**频域表示**：
$$\mathbf{Y} = \mathbf{W} \odot \mathbf{X}$$
其中 $\odot$ 表示逐元素乘法，$\mathbf{W}$ 和 $\mathbf{X}$ 为FFT变换后的频域表示。
**复杂度对比**：
| 方法 | 每样本计算量 |
|——|————-|
| 时域LMS | $O(L)$ |
| 频域FDAF | $O(\log L)$ |
## 5.3 MDF算法结构
**问题**：FDAF需要 $N = 2L$ 点FFT，对于长滤波器（如4096抽头），FFT点数过大。
**MDF解决方案**：将长滤波器分解为 $M$ 个较短子滤波器。
**滤波器分解**：
$$\mathbf{w} = [\mathbf{w}_1^T, \mathbf{w}_2^T, \ldots, \mathbf{w}_M^T]^T$$
每个子滤波器长度为帧长 $N/2$。
**频域表示**：
$$\mathbf{W} = [\mathbf{W}_1, \mathbf{W}_2, \ldots, \mathbf{W}_M]$$
**滤波输出**：
$$\mathbf{Y} = \sum_{j=1}^{M} \mathbf{X}_j \odot \mathbf{W}_j$$
**代码对应**：
“`matlab
st.Y(:, chan) = 0;
for j=1:M
st.Y(:, chan) = st.Y(:, chan) + st.X(:, speak, j) .* st.W(:, speak, j, chan);
end
“`
## 5.4 MDF算法公式推导
### 5.4.1 频域自适应滤波推导
**第一步：定义频域信号**
设输入信号帧为 $\mathbf{x}(k) = [x(k), x(k-1), \ldots, x(k-N+1)]^T$，经FFT变换：
$$\mathbf{X}(k) = \text{FFT}\{\mathbf{x}(k)\}$$
**第二步：频域滤波**
频域滤波输出：
$$\mathbf{Y}(k) = \mathbf{W}(k) \odot \mathbf{X}(k)$$
时域输出（IFFT后取有效部分）：
$$\mathbf{y}(k) = \text{IFFT}\{\mathbf{Y}(k)\}[N/2:N-1]$$
**第三步：频域误差**
期望信号帧 $\mathbf{d}(k)$ 与输出误差：
$$\mathbf{e}(k) = \mathbf{d}(k) – \mathbf{y}(k)$$
频域误差（零填充后FFT）：
$$\mathbf{E}(k) = \text{FFT}\{[\mathbf{0}, \mathbf{e}(k)]^T\}$$
**第四步：频域梯度**
时域LMS梯度：
$$\nabla_{\mathbf{w}} J = -e(k) \cdot \mathbf{x}(k)$$
频域梯度（利用循环卷积性质）：
$$\nabla_{\mathbf{W}} J = -\mathbf{E}(k) \odot \mathbf{X}^*(k)$$
**第五步：频域权值更新**
$$\mathbf{W}(k+1) = \mathbf{W}(k) + \mu \cdot \mathbf{E}(k) \odot \mathbf{X}^*(k)$$
### 5.4.2 归一化频域LMS (NFDAF)
**功率归一化**：
$$\mathbf{W}(k+1) = \mathbf{W}(k) + \frac{\mu}{\mathbf{P}(k)} \odot \mathbf{E}(k) \odot \mathbf{X}^*(k)$$
其中功率估计：
$$\mathbf{P}(k) = \lambda \mathbf{P}(k-1) + (1-\lambda) |\mathbf{X}(k)|^2$$
**代码对应**：
“`matlab
st.PHI = [st.power_1; st.power_1(end-1:-1:2)] .* st.prop(j) .* conj(st.X(:, speak, (j+1))) .* st.E(:, chan);
st.W(:, j) = st.W(:, j) + st.PHI;
“`
### 5.4.3 AUMDF变体
**循环卷积问题**：频域自适应滤波会引入循环卷积效应，导致失真。
**AUMDF解决方案**：交替更新各子滤波器，并强制约束。
**约束条件**：权值在时域的后半部分为零（重叠保存法要求）。
**更新规则**（代码中）：
“`matlab
if (j==1 || mod(2+st.cancel_count,(M-1)) == j)
st.wtmp = ifft(st.W(:, speak, j, chan));
st.wtmp(st.frame_size+1:N) = 0; % 强制约束
st.W(:, speak, j, chan) = fft(st.wtmp);
end
“`
**数学表示**：
$$\mathbf{w}_j \leftarrow \text{IFFT}\{\mathbf{W}_j\}$$
$$\mathbf{w}_j[N/2:N-1] = 0$$
$$\mathbf{W}_j \leftarrow \text{FFT}\{\mathbf{w}_j\}$$
# 6 双滤波器机制
## 6.1 前景/背景滤波器结构
代码采用双滤波器结构：
– **前景滤波器** `st.foreground`：用于输出
– **背景滤波器** `st.W`：用于自适应学习
“`matlab
st.W = zeros(N, K, M, C); % 背景滤波器
st.foreground = zeros(N, K, M, C); % 前景滤波器
“`
## 6.2 双滤波器原理
**背景滤波器**：
– 持续进行自适应更新
– 可能受到干扰（双讲、噪声）影响而发散
**前景滤波器**：
– 仅在背景滤波器性能更好时更新
– 提供稳定的输出
**性能评估**：
– $S_{ff}$：前景滤波器的误差能量
– $S_{ee}$：背景滤波器的误差能量
“`matlab
Sff = Sff + sum(abs(st.e(1:st.frame_size, chan)).^2); % 前景误差
See = See + sum(abs(st.e(1:st.frame_size, chan)).^2); % 背景误差
“`
## 6.3 滤波器切换逻辑
**更新条件**（背景→前景）：
$$\Delta^2 > S_{ff} \cdot D_{bf}$$
其中 $\Delta = S_{ff} – S_{ee}$，$D_{bf}$ 为两滤波器响应差异。
**代码实现**：
“`matlab
if (Sff-See)*abs(Sff-See) > (Sff*Dbf)
update_foreground = 1;
elseif (st.Davg1*abs(st.Davg1) > (VAR1_UPDATE*st.Dvar1))
update_foreground = 1;
…
end
“`
**统计验证**：
使用两个时间尺度的统计检验：
– $D_{avg1}$, $D_{var1}$：短时间尺度
– $D_{avg2}$, $D_{var2}$：长时间尺度
$$D_{avg} = \alpha \cdot D_{avg} + (1-\alpha) \cdot (S_{ff} – S_{ee})$$
$$D_{var} = \alpha^2 \cdot D_{var} + (1-\alpha)^2 \cdot S_{ff} \cdot D_{bf}$$
## 6.4 变学习率机制
**无显式双讲检测**，而是通过连续调整学习率实现鲁棒性。
**双讲场景分析**：
| 场景 | $S_{dd}$ | $S_{ee}$ | 学习率 |
|——|———-|———-|——–|
| 仅远端 | 正常 | 小 | 正常 |
| 双讲 | 大 | 大 | 自动减小 |
| 仅近端 | 大 | 大 | 自动减小 |
**机制**：
$$\alpha = \frac{\beta_0 \cdot S_{yy}}{S_{ee}}$$
双讲时 $S_{ee}$ 增大，$\alpha$ 减小，学习率降低。
### 6.4.1 泄漏估计
**核心概念**：泄漏估计 $\hat{\rho}$ 表示回声消除后残余回声与麦克风信号的比例。
**估计方法**：
基于误差信号与滤波输出的相关性：
$$P_{ey} = E[E \cdot Y^*]$$
$$P_{yy} = E[|Y|^2]$$
**泄漏估计**：
$$\hat{\rho} = \frac{P_{ey}}{P_{yy}}$$
**代码实现**：
“`matlab
Pey_cur = Pey_cur + sum(Eh_cur.*Yh_cur);
Pyy_cur = Pyy_cur + sum(Yh_cur.^2);
…
st.leak_estimate = st.Pey/st.Pyy;
“`
这段代码是 Speex MDF 算法中**最核心、最精妙的信号统计模块之一**。它的终极目标是：**实时、精准地估算当前的“回声泄漏比例”（即 `self.Pey / self.Pyy`）**。
这个泄漏比例不仅直接用于你之前问的 RER（理论回声返回损耗增强）计算，还决定了非线性处理（NLP）的抑制深度。
这段代码的设计思想极其高级，它没有直接用原始信号去算相关性，而是使用了**“去除直流分量的动态互相关”**加上**“双讲状态感知的自适应平滑”**。
阶段一：功率谱底噪保护与更新
“`python
See = max(See, N * 100)
for speak in range(K):
Sxx = Sxx + np.sum(self.x[self.frame_size:, speak]**2)
self.Xf = np.abs(self.X[:self.frame_size+1, speak, 0])**2
self.power = ss_1*self.power + 1 + ss*self.Xf[:,None]
“`
* **`See` 下界保护：** `See` 是误差信号的总功率。`N * 100` 是一个极低的能量门限。确保在完全静音时，`See` 不为 0，防止后续作为分母时引发除零崩溃（同样带有定点数防溢出的影子）。
* **远端时频域功率：** `Sxx` 累加了远端信号的时域能量；`self.Xf` 提取了当前帧远端的频域能量谱 $|X(k)|^2$。
* **长期功率平滑（EWMA）：** `self.power` 是远端频域功率的长期估计。
* **公式：** $P_{x,avg}(k) = \beta \cdot P_{x,avg}(k-1) + (1-\beta) \cdot |X(k)|^2 + 1$
* （加 `1` 是为了在定点运算中提供最小偏移量，防止静音时功率归一化失效）。
阶段二：核心技巧——提取“动态变化量”（去除均值）
“`python
Eh_cur = self.Rf – self.Eh
Yh_cur = self.Yf – self.Yh
# … (中间的计算) …
self.Eh = (1-self.spec_average)*self.Eh + self.spec_average*self.Rf
self.Yh = (1-self.spec_average)*self.Yh + self.spec_average*self.Yf
“`
* **原理（极其关键）：** 这里的 `self.Eh` 和 `self.Yh` 分别是误差信号和麦克风信号的**长期平滑频谱（相当于低频分量/直流偏置）**。
* `Eh_cur` 和 `Yh_cur` 是**当前帧减去长期平均值**，得到的是信号的**“高频突变分量”**。
* **为什么要这么做？** 在实际环境中，麦克风里包含恒定的风扇底噪、空调声。如果直接用原始信号算相关性，这些稳态噪声会严重干扰回声路径的估计。通过减去均值，算法**只对“正在发生变化”的声音（即远端刚说出来的话、回声路径的改变）敏感**，完美滤除了稳态背景噪声的干扰。
阶段三：计算瞬时归一化互相关
“`python
Pey_cur = Pey_cur + np.sum(Eh_cur*Yh_cur)
Pyy_cur = Pyy_cur + np.sum(Yh_cur**2)
Pyy = np.sqrt(Pyy_cur)
Pey = Pey_cur/Pyy
“`
* **公式推导：**
* $\text{Pey\_cur} = \sum (E_{dynamic} \cdot Y_{dynamic})$：这是误差变化量与麦克风变化量的**互相关**。只有当误差中的回声与麦克风中的回声同频同相变化时，这个值才大。
* $\text{Pyy\_cur} = \sum (Y_{dynamic}^2)$：这是麦克风变化量的**自能量**。
* **归一化：** $Pey = \frac{\text{Pey\_cur}}{\sqrt{\text{Pyy\_cur}}}$
* 这相当于做了类似协方差归一化的操作，消除了信号本身音量大小对相关性判断的影响，得到一个纯粹的“相似度”指标。
阶段四：双讲感知的自适应平滑因子 $\alpha$ 计算（神来之笔）
“`python
tmp32 = self.beta0*Syy
if tmp32 > self.beta_max*See:
tmp32 = self.beta_max*See
alpha = tmp32 / See
“`
算出了瞬间的 `Pey` 后，不能直接用，必须平滑（IIR滤波）到长期估计 `self.Pey` 中。但平滑速度 $\alpha$ 该怎么定？
* **常规逻辑：** 声音大（`Syy`大），更新快一点，所以 `tmp32 = beta0 * Syy`。
* **防双讲硬钳位（精髓）：** `if tmp32 > beta_max * See`
* 当发生**双讲**（近端也开始说话）时，误差 `See` 会瞬间激增（因为包含了消不掉的近端语音）。
* 此时如果平滑因子 $\alpha$ 还是很大，刚才算出的被近端语音污染的 `Pey` 就会迅速污染长期的 `self.Pey`，导致系统误以为回声路径剧变。
* **解决：** 利用激增的 `See` 强行把 $\alpha$ 压低（`alpha = tmp32 / See`，分母变大，$\alpha$ 变小）。**一旦检测到双讲（误差突增），立刻“冻结”泄漏估计的更新**，用过去的旧值扛过双讲期。这是极其优秀的鲁棒性设计。
阶段五：最终泄漏比例的更新与安全边界
“`python
alpha_1 = 1 – alpha
self.Pey = alpha_1*self.Pey + alpha*Pey
self.Pyy = alpha_1*self.Pyy + alpha*Pyy
if self.Pyy < 1: self.Pyy = 1
if self.Pey < MIN_LEAK * self.Pyy: self.Pey = MIN_LEAK * self.Pyy
if self.Pey > self.Pyy: self.Pey = self.Pyy
“`
* **IIR 平滑更新：** 标准的指数平滑公式，将瞬态估计融合到稳态估计中。
* **三层物理边界钳位（防算法跑飞）：**
1. `Pyy >= 1`：分母防零。
2. `Pey >= MIN_LEAK * Pyy`：**下界保护**。即使当前听起来完全没有回声，也强制假设存在极小比例（如 `MIN_LEAK=0.001`）的泄漏。如果不设下界，在远端停顿的间隙，`Pey` 会趋近于 0；等远端再次开口时，算法会“反应不过来”，导致开头几个字产生严重回声。
3. `Pey <= Pyy`：**上界保护（物理定律）**。误差和麦克风信号的互相关，在数学上绝对不可能超过麦克风信号自身的能量。强行钳位防止异常数值越界。
总结
这段代码不是在算波形，而是在算**“统计概率”**。它通过“去直流提取动态”、“误差感知冻结更新”、“强制最小泄漏假设”三大手段，在极其嘈杂和动态变化的声学环境中，死死咬住真实的回声泄漏比例，为后续的控制逻辑提供了最可靠的“方向盘”。
### 6.4.2 RER（残余误差比）
**定义**：
$$RER = \frac{S_{xx} + \hat{\rho} \cdot S_{yy}}{S_{ee}}$$
**物理意义**：
– 分子：估计的残余回声功率
– 分母：实际误差功率
– 理想情况下 $RER \approx 1$
**代码实现**：
“`matlab
RER = (.0001*Sxx + 3.*st.leak_estimate*Syy) / See;
“`
虽然变量名叫做 `RER`，但在这里它并不是指最终输出的“回声返回损耗增强”指标，而是指 **“当前状态下，理论上能够（且应该）达到的 RER 目标值”**。Speex 拿这个理论目标值作为杠杆，去反推和控制当前的步长（学习率 `prop`）。
回声路径泄漏估计
“`python
self.leak_estimate = self.Pey / self.Pyy
“`
* **公式翻译：** $\hat{g} = \frac{P_{ey}}{P_{yy}}$
* **原理：**
* $P_{yy}$ 是麦克风接收到的总信号（回声 + 近端语音 + 噪声）的功率。
* $P_{ey}$ 是误差信号 $e$ 与麦克风信号 $y$ 的互功率。
* 从线性代数角度看，这其实就是做了一个最简单的“线性回归”，估算出当前回声路径的**平均增益（泄漏比例）** $\hat{g}$。如果远端声音传过来在麦克风端剩下 20%，这个值就大约是 0.2。它是后续预测回声能量的基础。
构建理论 RER 目标（最核心公式）
“`python
RER = (.0001 * Sxx + 3. * self.leak_estimate * Syy) / See
“`
* **公式翻译：** $RER_{target} = \frac{0.0001 \cdot S_{xx} + 3 \cdot \hat{g} \cdot S_{yy}}{S_{ee}}$
* **原理拆解：**
* **分子（预测的回声功率）：** 算法试图猜测“现在麦克风里到底应该有多少纯粹的回声能量”。
* $3 \cdot \hat{g} \cdot S_{yy}$：利用刚才算出的泄漏比例 $\hat{g}$，乘以总功率 $S_{yy}$，再乘以一个经验常数 `3.`（这是一个为了补偿频域分块处理带来的方差而引入的过估计系数），得出预测回声能量。
* $0.0001 \cdot S_{xx}$：这是一个极小的底噪值。当远端完全静音（$S_{xx} \approx 0$）时，防止分子为 0 导致后续除零错误或步长计算异常。
* **分母（实际残余误差功率）：** $S_{ee}$ 是滤波器相减后实际输出的误差功率。
* **物理意义（反馈机制）：**
* **当 RER > 1：** 说明预测的回声能量 > 实际误差能量。意味着滤波器不仅把回声消除了，连近端语音或噪声都消掉了一部分（过度消除）。此时应该**减小学习率**。
* **当 RER < 1：** 说明实际误差比预测的回声大。意味着滤波器还没收敛好，或者回声路径变了。此时应该**增大学习率**。
引入维纳滤波器理论下界
“`python
if RER < Sey * Sey / (1 + See * Syy):
RER = Sey * Sey / (1 + See * Syy)
“`
* **公式翻译：** 如果 $RER < \frac{S_{ey}^2}{1 + S_{ee} \cdot S_{yy}}$，则 $RER = \frac{S_{ey}^2}{1 + S_{ee} \cdot S_{yy}}$
* **原理（非常硬核的数学限制）：**
* $\frac{S_{ey}^2}{S_{ee} \cdot S_{yy}}$ 在统计学上是什么？它是误差信号 $e$ 和麦克风信号 $y$ 之间的**相关系数的平方（$|\rho_{ey}|^2$）**。
* 根据维纳滤波器和线性预测理论，**线性滤波器最多只能消除与输入信号具有线性相关性的部分**。这个相关系数的平方，就是当前状态下线性 AEC 能够达到的**物理极限 ERLE**。
* 如果前面算出的 `RER` 比这个理论极限还要小，说明系统状态极度异常（比如发生了严重的双讲）。此时，强制将 `RER` 钳位到这个理论下界。
* *(注：分母里的 `1+` 同样是定点数防除零的保护措施)*。
强制稳定性上限
“`python
if RER > .5:
RER = .5
“`
* **原理：** 再次出现定点 DSP 的影子。`0.5` 在 Q15 格式中对应 `16384`。
* 无论如何计算，强制将理论 RER 上限限制在 0.5。这是为了从根本上限制后续推导出的学习率 $\mu$ 的最大值，**绝对防止自适应滤波器发散（发散会导致刺耳的尖叫声）**。
### 6.4.3 自适应学习率
**代码实现**：
“`
# if (not self.adapted and self.sum_adapt > M and self.leak_estimate*Syy > .03*Syy):
# self.adapted = 1
#
# if self.adapted:
# for i in range(self.frame_size+1):
# r = self.leak_estimate*self.Yf[i]
# e = self.Rf[i]+1
# if r > .5*e:
# r = .5*e
# r = 0.7*r + 0.3*(RER*e)
# self.power_1[i] = (r/(e*self.power[i]+10))
# else:
adapt_rate = 0
if Sxx > N * 1000:
tmp32 = 0.25 * Sxx
if tmp32 > .25*See:
tmp32 = .25*See
adapt_rate = tmp32 / See
self.power_1 = adapt_rate/(self.power+10)
“`
这段代码堪称 Speex MDF 学习率控制的**“全集”**。它展示了一个非常经典的自适应滤波器工程设计：**双模态变步长机制**。
你提供的代码中，上半部分被注释掉了（`#`），下半部分是实际执行的。这并非代码写错了，而是原作者（或移植者）**故意关闭了稳态的精细调节，强制让算法始终运行在“初始收敛”模式**。
为了让你彻底搞懂，我将同时揭秘被注释掉的“高级模式”和正在运行的“基础模式”，你会发现这背后藏着极高的工程智慧。
#### 6.4.3.1 稳态精细调节
被注释掉的代码：稳态精细调节模式（高级模式）
这段代码在 `self.adapted == True`（算法已收敛）时生效。它的最大特点是：**逐频点计算步长，并且引入了平滑滤波。**
1. 收敛状态切换开关
“`python
# if (not self.adapted and self.sum_adapt > M and self.leak_estimate*Syy > .03*Syy):
# self.adapted = 1
“`
* **原理：** 算法不能刚启动几帧就认为自己收敛了。必须满足三个条件才切换到稳态：
1. `sum_adapt > M`：累计更新次数超过阈值 M（保证经过了足够长的时间）。
2. `leak_estimate > 0.03`：回声泄漏估计大于 3%。如果小于 3%，说明压根没有回声（比如戴着耳机），此时不需要复杂的步长控制。
* **意义：** 防止算法在初始震荡期误入稳态模式。
2. 逐频点步长计算核心循环
“`python
# for i in range(self.frame_size+1):
# r = self.leak_estimate*self.Yf[i] # 预测的该频点回声功率
# e = self.Rf[i]+1 # 实际的该频点误差功率 (+1防除零，定点遗迹)
“`
* **原理：** 进入稳态后，不再用整个频段的总功率（`Sxx`, `See`）算一个全局步长，而是**针对每一个频率点 $i$** 单独算一个步长 `power_1[i]`。因为低频和高频的能量差异极大，必须区别对待。
3. 硬钳位与 IIR 平滑（最精华的部分）
“`python
# if r > .5*e:
# r = .5*e # 防发散硬上限
# r = 0.7*r + 0.3*(RER*e) # IIR 低通平滑滤波
“`
* **防发散：** 预测回声不能超过实际误差的一半，防止双讲时权重跑飞。
* **平滑滤波（极其重要）：** 上一问我们算出了理论目标 `RER`。如果直接用 `RER` 去算步长，由于语音信号的帧间波动，步长会剧烈抖动，导致滤波后的声音有“金属颤音”。
* **这里使用了一阶 IIR 滤波器：** `70% 的历史步长预测 + 30% 的当前 RER 目标`。这让步长的变化变得极其平滑（相当于给油门加了一个阻尼器）。
4. 稳态最终归一化公式
“`python
# self.power_1[i] = (r / (e * self.power[i] + 10))
“`
* **公式抽象：** $\mu(i) = \frac{\text{平滑后的目标更新量}}{\text{实际误差功率} \times \text{远端平均功率} + \text{底限}}$
* **原理：** 分母包含了误差 `e` 和远端功率 `power[i]` 的乘积，这是极其严格的**归一化**。它保证了无论频点能量多高，更新步伐都被死死按住。
**学习率调整公式**：
对于每个频率点 $i$：
$$\mu_i = \frac{r_i}{e_i \cdot P_i}$$
其中：
– $r_i = \hat{\rho} \cdot Y_i + 0.3 \cdot RER \cdot e_i$
– $e_i = R_i + 1$（误差功率）
– $P_i$：输入功率估计
– $RER$：残余误差比
#### 6.4.3.2 初始粗调
未注释的代码：初始粗调模式（当前实际运行）
由于高级模式被注释，当前的代码其实是一个**降级版的、更简单粗暴的初始收敛器**。它不管频点差异，只算一个全局标量步长。
这段代码是 Speex MDF 算法中，将上一阶段算出的理论状态（你之前问的 RER 相关逻辑），**最终转化为频域滤波器实际更新步长（`self.power_1`，即公式里的 $\mu$）**的核心执行代码。
如果说上一段 RER 代码是“参谋部”制定战略目标，这段代码就是“作战部”根据弹药库（信号能量）实际情况，扣动扳机前的最后校准。
“`python
adapt_rate = 0
if Sxx > N * 1000: # 1. 远端能量门槛
tmp32 = 0.25 * Sxx # 2. 基础驱动 (与远端能量成正比)
if tmp32 > .25*See: # 3. 防双讲硬钳位
tmp32 = .25*See
adapt_rate = tmp32 / See # 4. 算出比率
self.power_1 = adapt_rate / (self.power + 10) # 5. 归一化得出最终步长
“`
它完美展现了**“理论步长 + 功率归一化 + 防发散硬保护”**的工程哲学。我们逐段拆解：
第一阶段：能量门槛——没声音不更新
“`python
adapt_rate = 0
if Sxx > N * 1000:
“`
* **变量释义：** `Sxx` 是当前帧远端信号（扬声器出来的声音）在某个频段的功率；`N` 是频域变换的帧长（比如 160 或 256）。
* **原理：** 自适应滤波器不能在远端静音时更新。如果远端没声音，麦克风里的声音全是近端语音或环境噪声，此时更新会把噪声学进回声路径里。
* **`N * 1000` 的奥秘：** 这个阈值不是随便写的。它是一个与帧长成正比的能量底线。它确保远端信号不仅仅是微弱的底噪，而是真正的有效语音/音频。不满足条件，`adapt_rate` 直接锁死为 0，冻结滤波器更新，节省 CPU。
第二阶段：基于误差反馈的初始步长
“`python
tmp32 = 0.25 * Sxx
“`
* **原理：** 只要远端能量够大，就给一个基础步长，大小与远端能量成正比（这里取了 $0.25$ 倍，这是一个经过大量实验调参得出的经验常数）。
第三阶段：防双讲发散的“硬钳位”（极其关键）
“`python
if tmp32 > .25*See:
tmp32 = .25*See
“`
* **变量释义：** `See` 是当前误差信号（麦克风信号减去估计回声后的残留）的功率。
* **原理（安全气囊）：** 这是整个算法防发散的最重要一道物理防线。
* 正常情况下，收敛良好时，误差 `See` 很小，`tmp32` 不会被钳位。
* **异常情况（双讲）：** 当近端突然开始说话，误差 `See` 会瞬间飙升（因为包含了近端语音能量）。如果不加限制，巨大的误差乘以远端信号会把滤波器权重彻底带偏（算法发散，产生刺耳尖叫）。
* **钳位逻辑：** 强行规定，**更新步长的原始驱动力，绝对不能超过误差能量的 25%**。一旦双讲发生，`tmp32` 被强制拉低到与 `See` 同等量级，瞬间踩下刹车，牺牲收敛速度来保命（保系统稳定）。
第四阶段：计算归一化前的原始比率
“`python
adapt_rate = tmp32 / See
“`
* **原理：** 结合上一行代码，这里的数学结果只有两种可能：
1. **正常单讲：** `adapt_rate = (0.25 * Sxx) / See`。这是一个类似信噪比倒数的概念。远端越大、误差越小，这个值越大，更新越积极。
2. **异常双讲：** 因为上一行的钳位，`adapt_rate = (0.25 * See) / See = 0.25`。步长被死死卡在 0.25 这个安全上限。
第五阶段：频域功率归一化（得出最终 $\mu$）
“`python
self.power_1 = adapt_rate / (self.power + 10)
“`
* **变量释义：** `self.power` 是远端信号功率的长期平滑平均值（通常是一阶 IIR 递归求得）；`self.power_1` 就是最终送到滤波器更新公式里的**频域变步长 $\mu(k)$**。
* **原理（标准 NLMS 的频域体现）：** 频域 LMS 的更新公式是 $W_{new} = W_{old} + \mu \cdot X^* \cdot E$。如果不做归一化，能量大的频点更新极快，能量小的频点更新极慢。
* **除以 `self.power` 的作用：** 抹平不同频段的能量差异。使得无论低频（能量大）还是高频（能量小），滤波器系数的收敛轨迹是一致且平滑的。
* **为什么加 `10`？（又是 DSP 的遗迹）：**
* 在纯粹的数学里，这里应该加一个极小的浮点数（如 `1e-10`）防止除零。
* 但在 Speex 的 C 语言原生代码里，这是**定点数（Q15格式）**运算。分母如果太小，除法结果会直接溢出 16 位整数的最大值（32767），导致步长变成一个极其离谱的负数，算法瞬间崩溃。
* 加上整数 `10`（在 Q15 中代表一个微小的非零值），是一道绝对安全的底限，确保 `self.power_1` 永远在一个安全的数值范围内。
总结：这套公式的宏观性格
如果把这套公式拟人化，它是一个**极度保守、谨小慎微的老司机**：
1. **贼喊捉贼不追（`Sxx < N*1000`）：** 远端没动静，绝对不出手。
2. **顺风使劲踩（`tmp32 = 0.25*Sxx`）：** 远端声音大且误差小，加速更新。
3. **遇险急刹车（`tmp32 > 0.25*See` 时钳位）：** 一旦发现误差不对劲（疑似双讲），不管三七二十一，先把步长压到极低的安全值。
4. **弯道慢速过（`/ (self.power + 10)`）：** 不管油门给多大，最终出力都要经过平均功率的衰减和底线的保护，保证平稳不侧翻。
这种基于严密逻辑嵌套的硬保护机制，正是 Speex 能够在早年那些算力孱弱、容易溢出的廉价 DSP 芯片上稳定运行的根本原因。虽然到了 Python 浮点数时代，有些保护看起来有些“笨拙”，但其背后的信号处理鲁棒性设计思想依然是顶级的。
#### 6.4.3.3 总结对比
为什么要注释掉高级模式？（核心洞察）
你一定会问：既然上面那段带平滑、逐频点的代码更高级，为什么要注释掉？
这通常是因为在将 Speex 的 C 代码移植到 Python（或特定平台）时，遇到了**“过度优化”导致的稳定性问题**：
1. **初始态太脆弱：** 如果环境突变（比如突然拿起手机），`self.adapted` 可能没来得及切回 `False`。此时算法还在用稳态那种“慢吞吞、加平滑”的步长去更新，导致回声消除极慢，用户会听到好几秒的回声。
2. **强制激进策略：** 注释掉稳态模式，**强制算法永远保持在“初始粗调模式”**。这种模式对误差 `See` 极其敏感（一有误差就拼命更新），虽然会在稳态时引入一点微小的多余噪声（牺牲了一点极限 ERLE），但换来的是**极快的跟速能力**，在大多数实际通话场景中，用户体验反而更好。
3. **省去状态机维护：** 不用再去判断 `sum_adapt` 和 `leak_estimate` 的阈值，代码鲁棒性大幅提升。
总结对比表
| 维度 | 被注释的高级模式 (`self.adapted=1`) | 未注释的基础模式 (当前运行) |
| :— | :— | :— |
| **适用阶段** | 算法已收敛后的稳态微调 | 算法启动或强制全时运行 |
| **步长维度** | **向量** (逐频点计算 `power_1[i]`) | **标量** (全频段一个值 `power_1`) |
| **平滑处理** | 有 (70%历史 + 30%当前 RER) | 无 (直接根据当前帧计算，响应极快) |
| **分母归一化** | 包含误差功率 $e$ (极严格) | 不包含误差功率 $e$ (较宽松) |
| **优缺点** | 噪声极小，极限 ERLE 高；但反应迟钝 | 噪声略大，但跟迹速度极快，不易卡死 |
这段代码就像一辆同时装了“经济模式”和“赛车模式”的汽车。原作者把经济模式（高级模式）的线剪了，**让你永远踩在赛车模式的油门上**，用油耗（一点点稳态噪声）换取了随时随地地板油的推背感（极速抗回声路径变化能力）。
## 6.5 预处理与后处理
### 6.5.1 预加重
**目的**：降低信号动态范围，改善数值特性。
**预加重滤波器**：
$$x_{pre}(n) = x(n) – \alpha \cdot x(n-1)$$
**代码实现**：
“`matlab
st.preemph = .98;
…
tmp32 = far_end(i, speak) – st.preemph * st.memX(speak);
“`
### 6.5.2 直流陷波滤波器
**目的**：消除直流偏置。
**传递函数**：
$$H(z) = \frac{1 – z^{-1}}{1 – 2r \cos(\omega_0) z^{-1} + r^2 z^{-2}}$$
**代码实现**：
“`matlab
function [out,mem] = filter_dc_notch16(in, radius, len, mem)
for i=1:len
vin = in(i);
vout = mem(1) + vin;
mem(1) = mem(2) + 2*(-vin + radius*vout);
mem(2) = vin – (den2*vout);
out(i) = radius*vout;
end
end
“`
### 6.5.3 窗函数
**Hann窗**：
$$w(n) = 0.5 – 0.5 \cos\left(\frac{2\pi n}{N}\right)$$
**代码实现**：
“`matlab
st.window = .5-.5*cos(2*pi*((1:N)’-1)/N);
“`
**用途**：前景/背景滤波器切换时的平滑过渡。
## 6.6 比例自适应
### 6.6.1 比例因子
**目的**：对不同子滤波器分配不同的更新强度。
**初始化**：
“`matlab
decay = exp(-1/M);
st.prop(1, 1) = .7;
for i=2:M
st.prop(i, 1) = st.prop(i-1, 1) * decay;
end
st.prop = (.8 * st.prop)./sum(st.prop);
“`
**数学表示**：
$$p_i = \frac{0.8 \cdot p_i}{\sum_{j=1}^{M} p_j}$$
**物理意义**：较早的子滤波器（对应较近的回声）获得更大的更新强度。
### 6.6.2 自适应比例调整
**基于权值能量调整**：
“`matlab
function prop = mdf_adjust_prop(W, N, M, C, K)
for i=1:M
tmp = 1;
for chan=1:C
for speak=1:K
tmp = tmp + sum(abs(W(1:N/2+1, K, i, C)).^2);
end
end
prop(i) = sqrt(tmp);
end
max_sum = max(prop, 1);
prop = prop + .1*max_sum;
prop_sum = 1+sum(prop);
prop = .99*prop / prop_sum;
end
“`
**数学表示**：
$$p_i = \frac{0.99 \cdot (\sqrt{E_i} + 0.1 \cdot \max_j \sqrt{E_j})}{1 + \sum_j (\sqrt{E_j} + 0.1 \cdot \max_k \sqrt{E_k})}$$
其中 $E_i = \|\mathbf{W}_i\|^2$ 为第 $i$ 个子滤波器的能量。
## 6.7 完整算法流程
### 6.7.1 算法伪代码
“`
初始化:
设置帧长 N/2, 滤波器块数 M
初始化权值 W = 0, foreground = 0
初始化功率估计 P, 相关性 Pey, Pyy
For each frame k:
1. 预处理:
– DC陷波滤波
– 预加重
2. 频域变换:
X_j = FFT(x_j), j = 1, …, M
3. 前景滤波:
Y_fg = Σ_j X_j ⊙ W_fg_j
y_fg = IFFT(Y_fg)
e_fg = d – y_fg
S_ff = ||e_fg||^2
4. 背景滤波与更新:
Y_bg = Σ_j X_j ⊙ W_bg_j
y_bg = IFFT(Y_bg)
e_bg = d – y_bg
S_ee = ||e_bg||^2
For j = M to 1:
PHI = (1/P) ⊙ prop(j) ⊙ X_j^* ⊙ E
W_j = W_j + PHI
AUMDF约束更新
5. 滤波器切换:
If S_ff > S_ee (统计显著):
foreground = W (背景→前景)
Else if S_ff << S_ee:
W = foreground (前景→背景)
6. 变学习率更新:
计算泄漏估计 ρ
计算RER
更新功率归一化因子
7. 输出:
out = e_fg (去加重后)
End For
“`
“`matlab
% 主要数据结构
st.X % 输入频域信号 [N, K, M+1]
st.W % 背景滤波器权值 [N, K, M, C]
st.foreground % 前景滤波器权值 [N, K, M, C]
st.E % 误差频域信号 [N, C]
st.Y % 输出频域信号 [N, C]
st.power % 功率估计 [frame_size+1, 1]
st.power_1 % 归一化因子 [frame_size+1, 1]
st.prop % 比例因子 [M, 1]
“`
# 7 Speex MDF回声消除算法代码详细分析
– **算法依据**：Soo & Pang (1990) 的 MDF 算法，以及 Valin (2007) 关于频域回声消除中自适应学习率调整的论文。
– **核心思想**：将长滤波器分成多个短块，在频域进行自适应滤波，通过前景/背景滤波器切换和连续调整学习率来应对双讲。
– **输入输出**：
– 参考信号 `u`（扬声器信号）
– 麦克风信号 `d`（近端信号）
– 输出 `e`（消除回声后的信号，即误差）和 `y`（估计回声）
– **关键参数**：
– `fs`：采样率
– `frame_size`：帧长（通常 8ms，须为 2 的幂）
– `filter_length`：滤波器总长度（通常 250ms，须为 2 的幂）
代码整体结构
“`
┌─────────────────────────────────────────────────────────────┐
│ MDF 类结构 │
├─────────────────────────────────────────────────────────────┤
│ __init__() – 初始化所有参数和缓冲区 │
│ main_loop() – 主处理循环 │
│ speex_echo_cancellation_mdf() – 核心AEC处理函数 │
│ filter_dc_notch16() – DC陷波滤波器 │
│ mdf_adjust_prop() – 比例因子调整 │
└─────────────────────────────────────────────────────────────┘
“`
## 7.1 初始化函数分析
“`python
import numpy as np
“`
导入 NumPy 库用于数值计算。
“`python
def float_to_short(x):
x = x*32768.0
x[x < -32767.5] = -32768
x[x > 32766.5] = 32767
x = np.floor(0.5+x)
return x
“`
将浮点数（范围 [-1,1]）转换为 16-bit 整型（范围 -32768~32767），用于内部定点化处理。
### 7.1.1 基本参数设置
“`python
def __init__(self, fs: int, frame_size: int, filter_length: int) -> None:
nb_mic = 1
nb_speakers = 1
self.K = nb_speakers # 扬声器数量
K = self.K
self.C = nb_mic # 麦克风数量
C = self.C
“`
**参数说明**：
| 参数 | 含义 | 典型值 |
|——|——|——–|
| `fs` | 采样率 | 16000 Hz |
| `frame_size` | 帧大小 | 128 (8ms @ 16kHz) |
| `filter_length` | 滤波器长度 | 4096 (256ms @ 16kHz) |
构造函数，初始化所有状态变量。
“`python
nb_mic = 1
nb_speakers = 1
self.K = nb_speakers
K = self.K
self.C = nb_mic
C = self.C
“`
仅支持单通道（单麦克风、单扬声器），`K` 和 `C` 分别表示扬声器通道数和麦克风通道数。
### 7.1.2 MDF核心参数
“`python
self.frame_size = frame_size
self.filter_length = filter_length
self.window_size = frame_size*2 # N = 2 * frame_size
N = self.window_size
self.M = int(np.fix((filter_length+frame_size-1)/frame_size)) # 块数
“`
**数学关系**：
– FFT长度：$N = 2 \times \text{frame\_size}$
– 块数：$M = \lceil \frac{L + \text{frame\_size} – 1}{\text{frame\_size}} \rceil$
– `window_size`：FFT 点数，为帧长的两倍（重叠保留法）。
– `M`：滤波器块数，使总长度 `M * frame_size` 覆盖 `filter_length`。
**示例**：
“`
frame_size = 256, filter_length = 3072
N = 512 (FFT长度)
M = (3072 + 256 – 1) / 256 = 13 块
总滤波器长度 = 13 × 256 = 3328
“`
“`python
self.cancel_count = 0
self.sum_adapt = 0
self.saturated = 0
self.screwed_up = 0
“`
统计计数器：`cancel_count` 处理帧数；`sum_adapt` 累计自适应步长；`saturated` 指示麦克风饱和；`screwed_up` 指示异常状态。
“`python
self.sampling_rate = fs
self.spec_average = (self.frame_size)/(self.sampling_rate)
self.beta0 = (2.0*self.frame_size)/self.sampling_rate
self.beta_max = (.5*self.frame_size)/self.sampling_rate
self.leak_estimate = 0
“`
– `spec_average`：频谱平滑系数（约等于帧长/采样率）。
– `beta0` 和 `beta_max`：用于学习率调整的上下界。
– `leak_estimate`：泄漏估计（回声路径残余）。
### 7.1.3 数据结构初始化
“`python
self.e = np.zeros((N, C),)
self.x = np.zeros((N, K),)
self.input = np.zeros((self.frame_size, C),)
self.y = np.zeros((N, C),)
self.last_y = np.zeros((N, C),)
self.Yf = np.zeros((self.frame_size+1, 1),)
self.Rf = np.zeros((self.frame_size+1, 1),)
self.Xf = np.zeros((self.frame_size+1, 1),)
self.Yh = np.zeros((self.frame_size+1, 1),)
self.Eh = np.zeros((self.frame_size+1, 1),)
“`
分配时域和频域缓冲区：
– `e`、`x`、`y`、`last_y`：时域信号（长度 N，用于重叠保留）。
– `input`：当前帧的麦克风信号（经过预加重后）。
– `Yf`、`Rf`、`Xf`：频谱幅度平方（用于功率估计）。
– `Yh`、`Eh`：平滑后的频谱（用于学习率）。
“`python
# 频域信号存储
self.X = np.zeros((N, K, M+1), dtype=np.complex) # 输入频域信号
self.Y = np.zeros((N, C), dtype=np.complex) # 输出频域信号
self.E = np.zeros((N, C), dtype=np.complex) # 误差频域信号
self.W = np.zeros((N, K, M, C), dtype=np.complex) # 背景滤波器权值
self.foreground = np.zeros((N, K, M, C), dtype=np.complex) # 前景滤波器权值
“`
频域数据结构：
– `X`：参考信号的频域表示，大小为 `(FFT点数, 通道数, 块数+1)`，其中 `X[:, :, 0]` 为当前块的 FFT，`X[:, :, 1:M+1]` 存储历史块（用于 MDF 的延迟线）。
– `Y`：估计回声的频域表示。
– `E`：误差信号的频域表示。
– `W`：自适应滤波器系数（频域）。
– `foreground`：前景滤波器系数（用于双讲时切换）。
**维度解释**：
“`
X: [N, K, M+1]
N: FFT长度
K: 扬声器数
M+1: 延迟块数（+1用于存储历史）
W: [N, K, M, C]
N: FFT长度
K: 扬声器数
M: 滤波器块数
C: 麦克风数
“`
### 7.1.4 比例因子初始化
“`python
self.PHI = np.zeros((frame_size+1, 1),)
self.power = np.zeros((frame_size+1, 1),)
self.power_1 = np.ones((frame_size+1, 1),)
self.window = np.zeros((N, 1),)
self.prop = np.zeros((M, 1),)
self.wtmp = np.zeros((N, 1),)
“`
辅助数组：
– `PHI`：更新量。
– `power`：参考信号功率平滑值。
– `power_1`：归一化步长因子。
– `window`：汉宁窗（用于重叠保留法的时域加窗）。
– `prop`：每个块的自适应比例（用于分配更新步长）。
– `wtmp`：临时时域系数。
“`python
decay = np.exp(-2.4/M)
self.prop[0, 0] = .7
for i in range(1, M):
self.prop[i, 0] = self.prop[i-1, 0] * decay
self.prop = (.8 * self.prop) / np.sum(self.prop)
“`
**数学公式**：
$$\text{prop}_i = \frac{0.8 \cdot \text{prop}_i}{\sum_{j=0}^{M-1} \text{prop}_j}$$
其中初始值：
$$\text{prop}_0 = 0.7, \quad \text{prop}_i = \text{prop}_{i-1} \cdot e^{-2.4/M}$$
**物理意义**：较早的块（对应较近的回声）获得更大的更新强度。
初始化块更新比例 `prop`：指数衰减，总和归一化，用于控制不同块的自适应步长（靠前的块获得更大权重）。
### 7.1.5 窗函数
“`python
self.window = .5 – .5 * np.cos(2*np.pi*(np.arange(1, N+1).reshape(-1, 1)-1)/N)
“`
生成汉宁窗，用于重叠保留法的合成。
**Hann窗公式**：
$$w[n] = 0.5 – 0.5 \cos\left(\frac{2\pi n}{N}\right), \quad n = 0, 1, \ldots, N-1$$
**用途**：前景/背景滤波器切换时的平滑过渡。
### 7.1.6 其它初始化变量
“`python
self.memX = np.zeros((K, 1),)
self.memD = np.zeros((C, 1),)
self.memE = np.zeros((C, 1),)
self.preemph = .9
“`
预加重滤波器状态：`memX`、`memD`、`memE` 用于一阶高通滤波（去直流）和去加重。
“`python
if self.sampling_rate < 12000:
self.notch_radius = .9
elif self.sampling_rate < 24000:
self.notch_radius = .982
else:
self.notch_radius = .992
self.notch_mem = np.zeros((2*C, 1),)
“`
DC notch 滤波器参数，根据采样率选择极点半径。
“`python
self.adapted = 0
self.Pey = 1
self.Pyy = 1
self.Davg1 = 0
self.Davg2 = 0
self.Dvar1 = 0
self.Dvar2 = 0
“`
其他状态变量：
– `adapted`：是否已进入自适应状态（初始为 0，后续根据能量条件开启）。
– `Pey`、`Pyy`：用于泄漏估计的互相关和自相关。
– `Davg1/2`、`Dvar1/2`：用于双讲检测的统计量（均值和方差）。
## 7.2 主循环分析
### 7.2.1 main_loop函数
“`python
def main_loop(self, u, d):
“””MDF core function
Args:
u (array): reference signal (远端/扬声器信号)
d (array): microphone signal (近端/麦克风信号)
“””
assert u.shape == d.shape
u = float_to_short(u) # 转换为16位定点
d = float_to_short(d)
e = np.zeros_like(u) # 误差信号
y = np.zeros_like(u) # 滤波器输出
end_point = len(u)
for n in range(0, end_point, self.frame_size):
# 获取当前帧
u_frame = u[n:n+self.frame_size]
d_frame = d[n:n+self.frame_size]
# 核心处理
out = self.speex_echo_cancellation_mdf(d_frame[:, None], u_frame[:, None])[:,0]
e[n:n+self.frame_size] = out
y[n:n+self.frame_size] = d_frame – out
“`
“`python
e = e/32768.0
y = y/32768.0
return e, y
“`
输出转换为浮点数（范围 [-1,1]）后返回。
逐帧处理：
– 每帧调用 `speex_echo_cancellation_mdf`，输入为列向量（`[:, None]` 增加维度以保持通道维度）。
– 输出为当前帧的误差信号 `out`，保存到 `e`。
– 估计回声 `y` 为 `d_frame – out`。
– 打印处理进度（每帧）。
**处理流程图**：
“`
输入信号 u, d
│
↓
┌─────────────┐
│ 定点转换 │ float_to_short()
└─────────────┘
│
↓
┌─────────────┐
│ 分帧处理 │ for each frame
└─────────────┘
│
↓
┌─────────────────────────────┐
│ speex_echo_cancellation_mdf │
└─────────────────────────────┘
│
↓
输出 e (误差), y (估计回声)
“`
—
## 7.2 核心处理函数分析
核心回声消除函数 `speex_echo_cancellation_mdf`，此函数处理一帧数据，实现 MDF 算法的核心步骤。
### 7.2.1 函数签名与初始化
“`python
def speex_echo_cancellation_mdf(self, mic, far_end):
N = self.window_size # FFT长度
M = self.M # 块数
C = self.C # 麦克风数
K = self.K # 扬声器数
Pey_cur = 1 # 误差-输出相关性
Pyy_cur = 1 # 输出功率
out = np.zeros((self.frame_size, C))
self.cancel_count += 1
ss = .35/M # 平滑参数
ss_1 = 1 – ss
“`
初始化局部变量，递增帧计数器。
### 7.2.2 预处理：DC陷波与预加重
#### 7.2.2.1 DC陷波滤波器
实现一个二阶 IIR 陷波滤波器，用于去除直流分量。使用双二阶结构，滤波器系数根据采样率选择极点半径。
“`python
for chan in range(C):
self.input[:, chan], self.notch_mem[:, chan] = self.filter_dc_notch16(
mic[:, chan], self.notch_mem[:, chan])
“`
**filter_dc_notch16函数**：
“`python
def filter_dc_notch16(self, mic, mem):
out = np.zeros_like(mic)
den2 = self.notch_radius**2 + 0.7 * (1-self.notch_radius)**2
for i in range(self.frame_size):
vin = mic[i]
vout = mem[0] + vin
mem[0] = mem[1] + 2*(-vin + self.notch_radius*vout)
mem[1] = vin – (den2*vout)
out[i] = self.notch_radius * vout
return out, mem
“`
**传递函数分析**：
这是一个二阶IIR陷波滤波器，用于消除直流分量：
$$H(z) = \frac{1 – z^{-1}}{1 – 2r\cos(\omega_0)z^{-1} + r^2 z^{-2}}$$
在DC处 ($\omega_0 = 0$)：
$$H(z) = \frac{1 – z^{-1}}{1 – 2rz^{-1} + r^2 z^{-2}} = \frac{1 – z^{-1}}{(1 – rz^{-1})^2}$$
#### 7.2.2.2 预加重
对麦克风信号应用 DC notch 滤波器（去除直流分量）。
“`python
for i in range(self.frame_size):
tmp32 = self.input[i, chan] – (np.dot(self.preemph, self.memD[chan]))
self.memD[chan] = self.input[i, chan]
self.input[i, chan] = tmp32
“`
一阶预加重（高通）：`input = input – 0.9 * memD`，提升高频分量。
**预加重滤波器**：
$$x_{pre}[n] = x[n] – \alpha \cdot x[n-1]$$
其中 $\alpha = 0.9$（预加重系数）。
**作用**：降低信号动态范围，改善数值稳定性。
### 7.2.3 输入信号缓冲与频域变换
“`python
for speak in range(K):
for i in range(self.frame_size):
# 移位缓冲区
self.x[i, speak] = self.x[i+self.frame_size, speak]
# 预加重
tmp32 = far_end[i, speak] – np.dot(self.preemph, self.memX[speak])
self.x[i+self.frame_size, speak] = tmp32
self.memX[speak] = far_end[i, speak]
# 循环移位历史频域信号 将历史块右移，为当前块腾出位置（`axis=2` 对应块索引）。
self.X = np.roll(self.X, 1, axis=2)
# FFT变换
for speak in range(K):
#对当前参考信号块做 FFT，存入 `X[:,:,0]`（除以 N 是为了与后续逆变换匹配）。
self.X[:, speak, 0] = np.fft.fft(self.x[:, speak]) / N
“`
更新参考信号缓冲区 `x`（长度为 N）：
– 将旧数据前移（丢弃最旧的一帧）。
– 对新的参考信号进行预加重后放入尾部。
– 更新预加重状态。
**数据流图**：
“`
远端信号 far_end
│
↓
┌────────────────┐
│ 预加重 │
└────────────────┘
│
↓
┌────────────────┐
│ 移位缓冲区 x │ [x_old | x_new]
└────────────────┘
│
↓
┌────────────────┐
│ FFT变换 │ X[:, :, 0]
└────────────────┘
│
↓
┌────────────────┐
│ 滚动历史X │ X = roll(X, axis=2)
└────────────────┘
“`
### 7.2.4 前景滤波器计算
“`python
Sxx = 0
for speak in range(K):
Sxx = Sxx + np.sum(self.x[self.frame_size:, speak]**2)
self.Xf = np.abs(self.X[:self.frame_size+1, speak, 0])**2
“`
计算参考信号能量（只计算新数据部分，长度为 `frame_size`）并得到功率谱 `Xf`。
“`python
Sff = 0 # 前景滤波器误差能量
for chan in range(C):
self.Y[:, chan] = 0
for speak in range(K):
for j in range(M):
# 频域卷积累加
self.Y[:, chan] = self.Y[:, chan] + self.X[:, speak, j] * self.foreground[:, speak, j, chan]
# IFFT得到时域输出
self.e[:, chan] = np.fft.ifft(self.Y[:, chan]).real * N
# 重叠保存：取后frame_size个有效点
self.e[:self.frame_size, chan] = self.input[:, chan] – self.e[self.frame_size:, chan]
# 计算误差能量
Sff = Sff + np.sum(np.abs(self.e[:self.frame_size, chan])**2)
“`
用**前景滤波器** `foreground` 计算估计回声：
– 对每个麦克风通道，累加所有扬声器通道和所有块的贡献。
– 逆变换得到时域回声（长度为 N）。
– 截取后半部分（`frame_size` 个点）作为有效回声，与 `input` 相减得到误差 `e`（实际误差仅前 `frame_size` 个点）。
– 计算误差能量 `Sff`。
**频域滤波公式**：
$$\mathbf{Y}^{(fg)} = \sum_{p=0}^{M-1} \mathbf{X}_p \odot \mathbf{W}_p^{(fg)}$$
**时域输出（重叠保存）**：
$$\mathbf{y}^{(fg)} = \text{IFFT}\{\mathbf{Y}^{(fg)}\}[\text{frame\_size}:N]$$
“`python
if self.adapted:
self.prop = self.mdf_adjust_prop()
“`
如果已进入自适应状态，根据当前滤波器系数能量重新计算块更新比例（使能量大的块获得更大步长）。
### 7.2.5 权值更新（背景滤波器）
“`python
if self.saturated == 0:
for chan in range(C):
for speak in range(K):
for j in list(range(M)[::-1]): # 从最后一个块开始
# 计算梯度更新量
self.PHI = np.concatenate([self.power_1, self.power_1[-2:0:-1]], axis=0) * self.prop[j] * np.conj(self.X[:, speak, j+1])[:,None] * self.E[:, chan][:,None]
# 更新权值
self.W[:, speak, j, chan] = self.W[:, speak, j, chan] + self.PHI[:, 0]
else:
self.saturated -= 1
“`
如果麦克风未饱和，则更新**背景滤波器** `W`：
– 遍历所有通道、扬声器、块（反向遍历，避免使用已更新的块？此处是反向但实际影响不大）。
– 计算更新量 `PHI`：
– `self.power_1` 是归一化步长因子（与参考信号功率相关）。
– `self.prop[j]` 是块比例。
– `np.conj(self.X[:, speak, j+1])` 是参考信号频谱的共轭（因为 `X` 中索引 `j+1` 存储的是历史块，`j` 从 0 开始，所以用 `j+1` 对应正确的延迟）。
– `self.E[:, chan]` 是误差频谱。
– 更新 `W`。
– 若饱和，跳过更新并递减计数。
**更新公式**：
$$\mathbf{W}_p \leftarrow \mathbf{W}_p + \text{prop}_p \cdot \frac{1}{\mathbf{P}} \odot \mathbf{X}_p^* \odot \mathbf{E}$$
其中：
– $\mathbf{P}$：功率估计
– $\text{prop}_p$：比例因子
– $\mathbf{X}_p^*$：频域输入的复共轭
– $\mathbf{E}$：频域误差
### 7.2.6 AUMDF时域约束
“`python
for chan in range(C):
for speak in range(K):
for j in range(M):
# AUMDF: 交替约束
if j == 0 or (2 + self.cancel_count) % (M-1) == j:
self.wtmp = np.fft.ifft(self.W[:, speak, j, chan]).real
self.wtmp[self.frame_size:N] = 0 # 强制后N/2个点为零
self.W[:, speak, j, chan] = np.fft.fft(self.wtmp)
“`
**AUMDF约束规则**：
| 块索引 | 约束条件 |
|——–|———-|
| $j = 0$ | 总是约束 |
| $j \geq 1$ | 当 $(2 + k) \mod (M-1) = j$ 时约束 |
**目的**：减少计算量，同时保持稳定性。
对部分块（每帧仅更新一个块，交替更新）施加时域约束：
– 将频域系数逆变换回时域，将后半部分（`frame_size` 点）置零（强制滤波器长度不超过 `frame_size`，对应块长度）。
– 再变换回频域。这种交替更新（AUMDF）减少了计算量并改善了收敛性。
### 7.2.7 背景滤波器输出
“`python
self.Yf = np.zeros((self.frame_size+1, 1),)
self.Rf = np.zeros((self.frame_size+1, 1),)
self.Xf = np.zeros((self.frame_size+1, 1),)
Dbf = 0 # 前景/背景差异
See = 0 # 背景滤波器误差能量
for chan in range(C):
self.Y[:, chan] = 0
for speak in range(K):
for j in range(M):
self.Y[:, chan] = self.Y[:, chan] + self.X[:, speak, j] * self.W[:, speak, j, chan]
self.y[:, chan] = np.fft.ifft(self.Y[:, chan]).real * N
# 计算前景/背景差异
self.e[:self.frame_size, chan] = self.e[self.frame_size:N, chan] – self.y[self.frame_size:N, chan]
Dbf = Dbf + 10 + np.sum(np.abs(self.e[:self.frame_size, chan])**2)
# 计算背景滤波器误差
self.e[:self.frame_size, chan] = self.input[:, chan] – self.y[self.frame_size:N, chan]
See = See + np.sum(np.abs(self.e[:self.frame_size, chan])**2)
“`
用背景滤波器 `W` 重新计算回声（前景滤波器在后续可能被替换）。
更新误差：将上一帧保存的误差（`e[self.frame_size:N]`）减去当前背景滤波器回声的后半部分，得到新的误差（但随后又被覆盖？这里代码似乎有误）。实际上后面又重新计算了误差（用 `input – y`），所以 `Dbf` 计算的是旧误差能量。
### 7.2.8 双滤波器切换逻辑
“`python
# 统计检验更新
VAR1_UPDATE = .5
VAR2_UPDATE = .25
VAR_BACKTRACK = 4
MIN_LEAK = .005
#双讲检测参数。更新用于双讲检测的均值和方差（基于误差能量差异和背景噪声）。
self.Davg1 = .6*self.Davg1 + .4*(Sff – See)
self.Dvar1 = .36*self.Dvar1 + .16*Sff*Dbf
self.Davg2 = .85*self.Davg2 + .15*(Sff – See)
self.Dvar2 = .7225*self.Dvar2 + .0225*Sff*Dbf
update_foreground = 0
# 判断是否更新前景滤波器，判断是否需要将背景滤波器复制到前景（即认为当前背景滤波器优于前景）。条件基于能量差异和统计量。
if (Sff – See) * abs(Sff – See) > (Sff * Dbf):
update_foreground = 1
elif (self.Davg1 * abs(self.Davg1) > (VAR1_UPDATE * self.Dvar1)):
update_foreground = 1
elif (self.Davg2 * abs(self.Davg2) > (VAR2_UPDATE * self.Dvar2)):
update_foreground = 1
“`
**`Dvar1` 和 `Dvar2` 计算过程中的系数（0.36, 0.16, 0.7225, 0.0225），是分别对应均值计算（`Davg1` 和 `Davg2`）中指数加权移动平均（EWMA）系数的“平方”。**
这不是巧合，而是有着严格的数学推导过程。下面为您详细拆解这些系数是怎么来的，以及为什么公式最后是 `Sff * Dbf`。
1. 系数的数学推导：为什么是平方？
在统计学中，如果我们要计算一个随机变量 $X$ 的滑动方差，最标准的公式是：
$$ Variance(X) = E[X^2] – (E[X])^2 $$
在代码中，`Davg1` 实际上是在计算变量 $X = (Sff – See)$ 的期望值 $E[X]$，使用的是一阶低通滤波器（EWMA）：
$$ Davg1_{n} = 0.6 \cdot Davg1_{n-1} + 0.4 \cdot X_n $$
根据信号与系统或随机过程的原理，**如果一个信号经过增益为 $\alpha$ 的线性系统，其方差会变为原来的 $\alpha^2$ 倍**。
因此，`Davg1` 这个序列本身的方差（即围绕其均值的波动程度），可以通过将其自身的递推公式代入方差公式推导出来：
$$ Var(Davg1_n) = (0.6)^2 \cdot Var(Davg1_{n-1}) + (0.4)^2 \cdot Var(X_n) $$
将数字代入：
* $0.6^2 = 0.36$
* $0.4^2 = 0.16$
这就是为什么 `Dvar1` 的系数是 `0.36` 和 `0.16` 的根本原因！**`Dvar1` 实际上是在递归地估计 `Davg1` 的方差。**
同理，对于慢速通道：
* `Davg2` 的系数是 $0.85$ 和 $0.15$
* `Davg2` 的方差系数就是 $0.85^2 = 0.7225$ 和 $0.15^2 = 0.0225$
2. 为什么输入项是 `Sff * Dbf`？
在上面的推导中，我们得出公式最后一项应该是 $Var(X_n)$。但在代码中，这一项变成了 `Sff * Dbf`。这是 Speex 作者基于 AEC 物理场景做的一个绝妙近似。
首先明确变量含义（根据 Speex 源码上下文）：
* **`Sff`**: 远端信号（扬声器播放的声音）的能量。
* **`See`**: 误差信号（麦克风收到的声音 – 估计的回声）的能量。
* **`Dbf`**: 实际上是 `Sff` 的长期平滑值（即远端能量的缓慢变化估计）。在 Speex 源码中，它通常通过类似 `Dbf = 0.99*Dbf + 0.01*Sff` 的方式更新。
那么，$X = Sff – See$ 的方差 $Var(X)$ 怎么会变成 $Sff \times Dbf$ 呢？
1. **单讲时的理想状态**：在没有近端人说话（只有远端声音）时，AEC 的目标是把回声完全消除掉，此时 $See \approx 0$。
2. **变量退化**：此时 $X = Sff – 0 = Sff$。所以 $X$ 的方差就近似等于 $Sff$ 的方差，即 $Var(Sff)$。
3. **方差的近似**：对于能量信号，其方差近似等于其均值的平方，即 $Var(Sff) \approx (Sff)^2$。
4. **避免噪声放大的平滑处理**：直接用瞬时值的平方 $Sff^2$ 会导致方差估计极度不稳定（对突发噪声极其敏感）。因此，Speex 使用了一个平滑后的远端能量 `Dbf` 来代替其中一个 `Sff`。
5. **最终结论**：用 **`Sff * Dbf`** 来作为 $Var(X_n)$ 的低成本、高稳定性的替代品。它既能反映当前帧的能量突变，又包含了长期能量基准。
3. 为什么要有两组参数（1 和 2）？
在双讲检测（DTD）中，最难处理的是**“漏报”**（没检测到双讲导致近端声音被当成回声消掉）和**“误报”**（把回声没消干净的残差当成了近端声音）。Speex 用了快慢两个时间常数来平衡：
* **快速通道（Davg1 / Dvar1）：系数 0.6 / 0.4**
* 反应极快。当近端突然开口说话时，`See` 瞬间变大，`Sff – See` 瞬间变小。
* 它能敏锐地捕捉到这个突变，**用于“快速触发”双讲状态**。
* 缺点是容易受滤波器未收敛或突发噪声干扰而产生误报。
* **慢速通道（Davg2 / Dvar2）：系数 0.85 / 0.15**
* 记忆性强，反应迟钝。它代表了长期的稳态分布。
* 它的作用是**“防误报确认”**。只有当慢速通道的均值也明显偏离其方差边界时，才真正认定是双讲。
* 缺点是有延迟，但不影响整体体验，因为快速通道已经保护了初始的语音切片。
总结
这段代码背后的逻辑是：
1. 提取特征差值 $X = Sff – See$。
2. 用 EWMA 算 $X$ 的均值（快慢两套）。
3. **利用线性系统方差传播特性（系数平方），递归求解均值的方差。**
4. 在求方差时，利用 AEC 单讲时 $See \approx 0$ 的先验知识，用远端能量的乘积 `Sff * Dbf` 巧妙替代了复杂的真实方差计算，避免了乘方运算。
5. 最后通过比较 $Davg$ 和 $\sqrt{Dvar}$（代码未贴出，但在下方逻辑中）的比值，来判定是否发生双讲。
**切换逻辑分析**：
**条件1**：瞬时差异显著
$$(S_{ff} – S_{ee})^2 > S_{ff} \cdot D_{bf}$$
**条件2**：短时间尺度统计显著
$$D_{avg1}^2 > 0.5 \cdot D_{var1}$$
**条件3**：长时间尺度统计显著
$$D_{avg2}^2 > 0.25 \cdot D_{var2}$$
“`python
if update_foreground:
# 背景滤波器更好，更新前景
self.foreground = self.W
# 平滑过渡
for chan in range(C):
self.e[self.frame_size:N, chan] = (self.window[self.frame_size:N][:,0] * self.e[self.frame_size:N, chan]) + (
self.window[:self.frame_size][:,0] * self.y[self.frame_size:N, chan])
else:
reset_background = 0
if (-(Sff-See)*np.abs(Sff-See) > VAR_BACKTRACK*(Sff*Dbf)):
reset_background = 1
if ((-self.Davg1 * np.abs(self.Davg1)) > (VAR_BACKTRACK*self.Dvar1)):
reset_background = 1
if ((-self.Davg2 * np.abs(self.Davg2)) > (VAR_BACKTRACK*self.Dvar2)):
reset_background = 1
if reset_background:
self.W = self.foreground
for chan in range(C):
self.y[self.frame_size:N, chan] = self.e[self.frame_size:N, chan]
self.e[:self.frame_size, chan] = self.input[:, chan] – self.y[self.frame_size:N, chan]
See = Sff
self.Davg1 = 0
self.Davg2 = 0
self.Dvar1 = 0
self.Dvar2 = 0
“`
如果更新前景，则重置统计量，将背景滤波器赋值给前景，并合成新的误差信号（用于下一步输出？这段代码实际上是更新 `e` 的尾部，用于下一帧的误差）。
如果不更新前景，则检查是否需要将背景重置为前景（即前景优于背景）。若满足条件（负向能量差异大），则重置背景为前景，并更新误差和统计量。
这段代码是 Speex AEC 中**双滤波器结构**的核心灵魂所在。
要理解这段代码，首先要明白“为什么要用双滤波器”，然后才能看懂“什么时候切换”。
一、 核心概念：为什么要“双滤波器”？
在声学回声消除（AEC）中，存在一个致命矛盾：
* **为了跟踪回声路径的变化**（比如人拿着手机走动，或者改变握姿），自适应滤波器（比如 NLMS）需要**不断地、快速地更新**。
* **但是，在双讲期间**（近端人和远端人同时说话），误差信号 `See` 里面混入了近端人的语音。如果这时候滤波器还在盲目更新，它会把近端语音当成“没消干净的回声”去学习，导致**滤波器发散**，输出声音破裂。
**双滤波器的解决思路是“主从分离”：**
1. **背景滤波器**：在后台默默运行，不管三七二十一，一直根据算法死命更新。它可能很准，也可能在双讲时发散。
2. **前景滤波器**：真正负责干活，计算最终输出的滤波器。它**绝对不能发散**。
3. **切换逻辑**：只有当我们**确认当前没有双讲（安全期）**，且后台的背景滤波器表现比前台好时，才把后台的权重复制给前台。
二、 代码逻辑解析：判断“是否安全”
这段代码的目的就是判断**“现在是不是安全的无双讲状态，可以更新前景滤波器了？”**
代码中的 `update_foreground = 1` 意思是：“安全！把背景滤波器的结果复制给前景滤波器！”
我们结合上一题的背景，来看这三个条件。注意一个数学等式：`X * abs(X)` 其实就是 $X^2$。
条件 1：瞬时/短时能量判断（最灵敏）
“`python
if (Sff – See) * abs(Sff – See) > (Sff * Dbf):
“`
* **数学转换**：$(Sff – See)^2 > Sff \times Dbf$
* **物理意义**：上一题我们推导过，在单讲（只有回声）且滤波器收敛良好时，误差 $See \approx 0$，此时 $(Sff – 0)^2 = Sff^2$。由于 $Dbf$ 是 $Sff$ 的平滑值，此时左边约等于右边。
* **触发逻辑**：如果发生双讲，近端声音加入导致 $See$ 瞬间变大，$(Sff – See)^2$ 会急剧**缩小**，条件就不成立。
* **结论**：当这个条件成立时，说明当前的误差能量 $See$ 非常小，远端能量完全主导了误差，**证明当前绝对没有双讲**。此时可以放心大胆地切换。
条件 2：快速统计量判断（抗突发干扰）
“`python
elif (self.Davg1 * abs(self.Davg1) > (VAR1_UPDATE * self.Dvar1)):
“`
* **数学转换**：$Davg1^2 > VAR1\_UPDATE \times Dvar1$
* **物理意义**：这就是条件 1 的“平滑升级版”。`Davg1` 是 $(Sff – See)$ 的快速平滑值，`Dvar1` 是它的方差（上一题证明过它近似等于 $Sff \times Dbf$）。
* **触发逻辑**：`VAR1_UPDATE` 是一个小于 1 的常数（类似阈值系数）。如果不使用瞬时值，而是用统计量，即使偶尔有一两帧的噪声导致 $See$ 突然变大，也不会立刻打断切换过程。只要快速统计量证明“近期整体误差很小”，就触发切换。
条件 3：慢速统计量判断（终极把关）
“`python
elif (self.Davg2 * abs(self.Davg2) > (VAR2_UPDATE * self.Dvar2)):
“`
* **数学转换**：$Davg2^2 > VAR2\_UPDATE \times Dvar2$
* **物理意义**：使用时间常数极大的慢速平滑（0.85/0.15）。
* **触发逻辑**：这是最严格的条件。它不看重眼前的几帧，它看重的是过去很长一段时间的稳态表现。只有当长期来看，误差都维持在一个很低的水平，它才会放行。
三、 整体逻辑串联（如果用文字描述这段代码）
这段代码在每一帧运行时，都会像保安一样问三个问题：
1. **保安（看当前帧）**：“这一瞬间，误差是不是几乎为零？” $\rightarrow$ 是的话，放行（条件1）。
2. **保安（看最近几帧）**：“这一瞬间可能有点风噪，但看最近几十毫秒的统计，误差是不是一直压制得很好？” $\rightarrow$ 是的话，放行（条件2）。
3. **保安（看过去几秒）**：“最近几秒的长效统计数据怎么样？系统是不是处于长期稳定的无近端语音状态？” $\rightarrow$ 是的话，放行（条件3）。
只要这三个问题**有一个**回答了“是”，系统就认为：**现在很安全，背景滤波器肯定没有因为双讲而发散，它现在学到的回声路径是最新的、最准的。立刻把它复制给前景滤波器！**
四、 补充细节：为什么要写 `X * abs(X)` 而不写 `X**2`？
这是一个底层 C 语言代码翻译到 Python 时留下的痕迹。
在早期 DSP 芯片或没有硬件浮点乘法器的环境下，计算浮点数的乘方（尤其是平方）是非常耗时的。而计算 `X * abs(X)` 可以利用硬件的符号判断指令快速完成，这在定点数运算中也是一种常见的求平方的替代写法，Speex 作为极度追求性能的 C 库，经常使用这类技巧。
#### 7.2.8.1 使用窗函数计算误差信号
在 Speex 的 MDF（多延迟块频域）算法中，计算误差信号 $e$ 时使用窗函数（如你提供的代码所示），其核心目的只有四个字：**消除块效应**。
具体来说，这涉及到了信号处理中经典的 **重叠相加法**。
下面为你深度解析这段代码背后的原理、数学逻辑以及为什么要这么做。
一、 问题的根源：为什么频域处理会有“咔哒声”？
MDF 是一种分块算法。它不是每进来一个采样点就处理一次，而是攒够一帧（比如 $N=256$ 个点），做一次 FFT，在频域算完结果后，再 IFFT 变回时域。
在标准的“重叠保存法”中，每一帧输出结果的前半段（由于圆周卷积的混叠效应）是要被丢弃的，只保留后半段（`self.frame_size:N`）作为有效输出。
**痛点在于：**
因为自适应滤波器的权重 $W$ 是每一帧都在不断更新的（为了追踪回声路径的变化），这就导致：
* 第 $k$ 帧输出的尾巴，是用权重 $W_k$ 算出来的。
* 第 $k+1$ 帧输出的开头，是用权重 $W_{k+1}$ 算出来的。
如果 $W_k$ 和 $W_{k+1}$ 有差异，那么这两帧在拼接的边界处，信号在时域上**是不连续的（发生突变）**。
人耳对这种微小的时域相位/幅度突变极其敏感，听起来就像是一种连续的“嘎啦嘎啦”的杂音，或者像是声音被蒙在玻璃后面（梳状滤波效应）。
二、 解决方案：加窗重叠相加（OLA）
为了掩盖这种块边界的不连续，Speex 引入了窗函数平滑过渡。
来看你提供代码中最核心的一行：
“`python
self.e[self.frame_size:N, chan] = (
self.window[self.frame_size:N][:,0] * self.e[self.frame_size:N, chan]) + (
self.window[:self.frame_size][:,0] * self.y[self.frame_size:N, chan])
“`
**公式抽象：**
对于输出的第 $n$ 个点（$n$ 位于当前帧的后半段），其最终的误差值 $e_{out}[n]$ 为：
$$ e_{out}[n] = w_{decay}[n] \cdot e_{old}[n] + w_{rise}[n] \cdot y_{new}[n] $$
**变量解析：**
* `self.y[self.frame_size:N]`：这是用**最新**的滤波器权重算出来的、刚刚 IFFT 出来的**新误差信号**。
* `self.e[self.frame_size:N]`：这是**上一帧**（旧帧）遗留在这个缓冲区里的**旧误差信号**。
* `self.window[:self.frame_size]`：窗函数的前半段（通常是**升窗/渐强**，如从 0 渐变到 1）。
* `self.window[self.frame_size:N]`：窗函数的后半段（通常是**降窗/渐弱**，如从 1 渐变到 0）。
**物理过程：**
这就好比音乐播放器里的**“交叉淡入淡出”**。在两帧交接的边界处，我不直接硬切，而是让旧声音慢慢变小（乘降窗），同时让新声音慢慢变大（乘升窗）。这样拼接出来的波形就永远是平滑的，彻底消灭了“咔哒声”。
—
三、 窗函数的数学约束：能量守恒
在这个代码中使用的窗函数（通常是正弦窗或汉宁窗的变体），必须满足一个严格的数学条件：
$$ w_{rise}^2[n] + w_{decay}^2[n] = 1 $$
**为什么要满足这个条件？**
因为声音的能量是与幅度的平方成正比的。
在交接处，旧信号衰减掉的能量，必须完全由新信号增加的能量来补齐。如果不满足这个条件，交叉淡入淡出会导致交接处的总能量发生凹陷（听起来像缺了一块）或凸起（听起来有爆音）。满足平方和为 1 的 OLA，保证了时域能量的完美无缝衔接。
—
四、 为什么是在 `update_foreground` 里做？
看这段代码的上下文，它包裹在 `if update_foreground:` 之中。
Speex MDF 使用了著名的**双滤波器结构**：
1. **后台滤波器：** 疯狂地、不顾吃相地更新权重，尝试追赶回声路径。
2. **前台滤波器：** 权重平时不动。只有当判断后台滤波器比前台好时，才把后台的权重拷贝给前台（`self.foreground = self.W`）。
**前台滤波器输出的信号是真正送给喇叭或网络的，因此对听感要求极高。**
当发生权重切换（`update_foreground = 1`）时，新旧权重的差异是最大的，此时如果不做 OLA 加窗平滑，必然会产生巨大的咔哒声。所以，Speex 强制在前台更新输出的这一刻，使用窗函数对即将输出的误差/残留信号进行一次完美的“缝合”。
总结
这段代码表面上是在算误差 $e$，实际上是在执行一个**高质量的重采样/波形拼接器**的功能。它通过加窗重叠相加（OLA），牺牲了一点点算法的复杂度（多几次乘法），换取了频域自适应滤波器在时域输出上的绝对平滑，是 Speex 能够达到“商用级听感”的关键细节之一。
### 7.2.9 输出处理与去加重
“`python
Sey = 0
Syy = 0
Sdd = 0
for chan in range(C):
for i in range(self.frame_size):
# 计算输出（去加重）
tmp_out = self.input[i, chan] – self.e[i + self.frame_size, chan]
tmp_out = tmp_out + self.preemph * self.memE[chan]
# 检测饱和
if mic[i, chan] <= -32000 or mic[i, chan] >= 32000:
if self.saturated == 0:
self.saturated = 1
out[i, chan] = tmp_out[0]
self.memE[chan] = tmp_out
“`
计算最终输出：
– 误差信号 = 预加重后的麦克风信号 – 背景滤波器估计的回声（`e[i+self.frame_size]` 是之前保存的误差尾部）。
– 对误差进行去加重（`+ 0.9 * memE`）恢复原始动态。
– 检测麦克风是否饱和（幅度超过 ±32000）。
“`python
self.e[self.frame_size:N, chan] = self.e[:self.frame_size, chan]
self.e[:self.frame_size, chan] = 0
“`
将当前误差保存到缓冲区尾部（为下一帧做准备），清空前部。
“`python
Sey = Sey + np.sum(self.e[self.frame_size:N, chan] * self.y[self.frame_size:N, chan])
Syy = Syy + np.sum(self.y[self.frame_size:N, chan]**2)
Sdd = Sdd + np.sum(self.input**2)
“`
计算互相关和自相关（用于后续学习率）。
“`python
self.E = np.fft.fft(self.e,axis=0) / N
self.y[:self.frame_size, chan] = 0
self.Y = np.fft.fft(self.y,axis=0) / N
self.Rf = np.abs(self.E[:self.frame_size+1, chan])**2
self.Yf = np.abs(self.Y[:self.frame_size+1, chan])**2
“`
频域变换，获取功率谱用于平滑。
“`python
if not (Syy >= 0 and Sxx >= 0 and See >= 0):
self.screwed_up = self.screwed_up + 50
out = np.zeros_like(out)
elif Sff > Sdd + N * 10000:
self.screwed_up = self.screwed_up + 1
else:
self.screwed_up = 0
if self.screwed_up >= 50:
print(“Screwed up, full reset”)
self.__init__(self.sampling_rate, self.frame_size, self.filter_length)
“`
检测异常状态（负能量或误差能量异常大），若连续异常则完全重置滤波器（重新初始化）。
### 7.2.10 变学习率计算
“`python
See = max(See, N * 100)
for speak in range(K):
Sxx = Sxx + np.sum(self.x[self.frame_size:, speak]**2)
self.Xf = np.abs(self.X[:self.frame_size+1, speak, 0])**2
self.power = ss_1*self.power + 1 + ss*self.Xf[:,None]
“`
确保 `See` 有下限，计算参考信号功率平滑值 `power`（用于归一化步长）。
“`python
# 计算相关性
Eh_cur = self.Rf – self.Eh
Yh_cur = self.Yf – self.Yh
Pey_cur = Pey_cur + np.sum(Eh_cur * Yh_cur)
Pyy_cur = Pyy_cur + np.sum(Yh_cur**2)
# 更新平滑估计
self.Eh = (1 – self.spec_average) * self.Eh + self.spec_average * self.Rf
self.Yh = (1 – self.spec_average) * self.Yh + self.spec_average * self.Yf
# 计算泄漏估计
Pyy = np.sqrt(Pyy_cur)
Pey = Pey_cur / Pyy
# 自适应平滑因子
alpha = tmp32 / See
alpha_1 = 1 – alpha
self.Pey = alpha_1 * self.Pey + alpha * Pey
self.Pyy = alpha_1 * self.Pyy + alpha * Pyy
# 泄漏估计
self.leak_estimate = self.Pey / self.Pyy
“`
更新平滑频谱，计算当前帧的互相关和自相关，用于后续泄漏估计。
计算泄漏估计因子 `leak_estimate = Pey/Pyy`，通过平滑更新。
**泄漏估计原理**：
$$\rho = \frac{P_{ey}}{P_{yy}} = \frac{E[E \cdot Y^*]}{E[|Y|^2]}$$
表示残余回声与滤波器输出的相关性。
### 7.2.11 学习率计算
“`python
if self.Pyy < 1:
self.Pyy = 1
if self.Pey < MIN_LEAK * self.Pyy:
self.Pey = MIN_LEAK * self.Pyy
if self.Pey > self.Pyy:
self.Pey = self.Pyy
self.leak_estimate = self.Pey/self.Pyy
“`
限制泄漏估计范围。
“`python
RER = (.0001*Sxx + 3.*self.leak_estimate*Syy) / See
if RER < Sey*Sey/(1+See*Syy):
RER = Sey*Sey/(1+See*Syy)
if RER > .5:
RER = .5
“`
计算残余回声比（RER），用于后续可能的学习率调整（但后续代码被注释了，实际使用另一种方式）。
“`python
# 学习率设置
adapt_rate = 0
if Sxx > N * 1000:
tmp32 = 0.25 * Sxx
if tmp32 > .25 * See:
tmp32 = .25 * See
adapt_rate = tmp32 / See
self.power_1 = adapt_rate / (self.power + 10)
“`
计算自适应步长因子 `adapt_rate`，并归一化到 `power_1`（用于更新 `W` 的归一化步长）。
**RER (Residual-to-Error Ratio)**：
$$RER = \frac{0.0001 \cdot S_{xx} + 3 \cdot \rho \cdot S_{yy}}{S_{ee}}$$
“`python
self.sum_adapt = self.sum_adapt+adapt_rate
self.last_y[:self.frame_size] = self.last_y[self.frame_size:N]
if self.adapted:
self.last_y[self.frame_size:N] = mic-out
“`
更新自适应总和，记录上一帧的回声（用于后续？）。
“`python
return out
“`
返回当前帧的输出误差。
## 7.3 比例因子调整函数
根据当前滤波器系数的能量分布更新块比例 `prop`，使能量大的块获得更大的步长，以提高收敛速度。
“`python
def mdf_adjust_prop(self):
N = self.window_size
M = self.M
C = self.C
K = self.K
prop = np.zeros((M, 1))
for i in range(M):
tmp = 1
for chan in range(C):
for speak in range(K):
tmp = tmp + np.sum(np.abs(self.W[:N//2+1, speak, i, chan])**2)
prop[i] = np.sqrt(tmp)
max_sum = np.maximum(prop, 1)
prop = prop + .1 * max_sum
prop_sum = 1 + np.sum(prop)
prop = 0.99 * prop / prop_sum
return prop
“`
**公式**：
$$\text{prop}_i = \frac{0.99 \cdot (\sqrt{E_i} + 0.1 \cdot \max(1, \sqrt{E_{max}}))}{1 + \sum_j (\sqrt{E_j} + 0.1 \cdot \max(1, \sqrt{E_{max}}))}$$
其中 $E_i = \|\mathbf{W}_i\|^2$ 为第 $i$ 个块的能量。
## 7.4 使用示例
“`python
if __name__ == “__main__”:
import soundfile as sf
# 读取音频文件
mic, sr = sf.read(“./glass_wav/wav1_mic.wav”) # 麦克风信号
ref, sr = sf.read(“./glass_wav/wav1_ref.wav”) # 远端参考信号
min_len = min(len(mic), len(ref))
mic = mic[:min_len]
ref = ref[:min_len]
# 创建MDF处理器
# sr=16000, frame_size=256 (16ms), filter_length=3072 (192ms)
processor = MDF(sr, 256, 3072)
# 执行回声消除
e, y = processor.main_loop(ref, mic)
# 保存结果
sf.write(‘./output_1/result.wav’, e, sr)
“`
**参数选择**：
– 采样率 16kHz
– 帧大小 256 样本 (16ms)
– 滤波器长度 3072 样本 (192ms)
– 块数 M = 13
– FFT长度 N = 512
## 7.5 关键算法特性总结
代码特点
| 特点 | 实现方式 |
|——|———-|
| **定点兼容** | float_to_short转换 |
| **多通道支持** | C个麦克风，K个扬声器 |
| **预处理** | DC陷波 + 预加重 |
| **平滑过渡** | Hann窗加权 |
| **异常处理** | 饱和检测、自动重置 |
### 7.5.1 双滤波器机制
“`
┌─────────────────────────────────────────────────────────────┐
│ 双滤波器结构 │
├─────────────────────────────────────────────────────────────┤
│ │
│ 背景 │ 前景 │
│ ↓ ↓ │
│ 持续自适应 稳定输出 │
│ ↓ ↓ │
│ 可能受干扰 性能可靠 │
│ ↓ ↓ │
│ └────────── 性能比较 ──────────┘ │
│ ↓ │
│ 背景更好 → 更新前景 │
│ 前景更好 → 重置背景 │
│ │
└─────────────────────────────────────────────────────────────┘
“`
### 7.5.2 AUMDF约束策略
| 特性 | 说明 |
|——|——|
| **交替更新** | 每帧只约束部分块 |
| **减少计算** | 避免所有块同时约束 |
| **保持稳定** | 块0总是约束 |
### 7.5.3 变学习率机制
**无需显式双讲检测**，通过以下机制实现鲁棒性：
1. **泄漏估计**：$\rho = P_{ey}/P_{yy}$
2. **RER计算**：估计残余回声比例
3. **自适应调整**：根据信号特性调整学习率
# 8 FFT缩放
在 Speex（以及许多其他音频/信号处理代码）中，看到正向 FFT 除以 $N$，逆向 IFFT 乘以 $N$，**这不是代码写错了，而是一种刻意的“缩放约定”**。
这种写法的核心目的有两个：**1. 统一频域幅度的物理意义；2. 历史遗留的定点运算防溢出机制。**
其中的数学逻辑和工程原因：
## 8.1 数学上的完美抵消
首先，我们要知道 NumPy（以及大部分标准数学库）对 FFT 的定义：
* **FFT:** $X[k] = \sum_{n=0}^{N-1} x[n] e^{-j 2\pi kn/N}$ （不缩放）
* **IFFT:** $x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j 2\pi kn/N}$ （默认带有 $\frac{1}{N}$ 的缩放）
现在来看这段代码的执行流程：
1. **频域变换：** `X = np.fft.fft(x) / N`
2. **时域还原：** `y = np.fft.ifft(X) * N`
把 1 代入 2：
$$y = \text{IFFT}\left( \frac{\text{FFT}(x)}{N} \right) \times N$$
因为 NumPy 的 IFFT 内部自带一个 $\frac{1}{N}$，所以：
$$y = \left( \frac{\text{FFT}(x)}{N} \times \frac{1}{N} \right) \times N = \frac{\text{FFT}(x)}{N}$$
等等，这样好像没有还原？
**注意：** $\text{IFFT}(\text{FFT}(x)) = x$。所以：
$$y = \left( \frac{x}{N} \right) \times N = x$$
**结论：** 除以 $N$ 和乘以 $N$ 在数学上**完美抵消**了，信号可以无失真地在时域和频域之间来回转换。
## 8.2 为什么要这么做？（核心原因）
既然不乘除 $N$ 也能正常工作，为什么 Speex 要多此一举？
**原因 A：让频域幅度代表“平均值”而非“总和”（与帧长解耦）**
标准的 FFT 结果是**求和**。如果一段音频能量不变，但你把分析帧长 $N$ 从 160（10ms@16kHz）增加到 320（20ms@16kHz），标准 FFT 出来的频域幅度会直接翻倍（因为多加了 160 个点的值）。
* **除以 $N$ 之后**，FFT 的结果变成了**求平均**。此时频域的幅度不再受帧长 $N$ 的影响，直接反映了该频率在当前帧内的“平均能量”。
* 在 Speex 后续的处理中（比如计算噪声估计、计算感知加权、LPC 分析等），使用“平均幅度”可以让算法的阈值、权重等超参数与窗口大小无关，调参更加直观和鲁棒。
**原因 B：Speex 原本是为“定点 DSP”设计的（历史遗留）**
这是最硬核的工程原因。Speex 最早是用 C 语言写的，专门跑在计算能力极弱、没有浮点运算单元（FPU）的低端嵌入式芯片上，使用的是 **定点数（Q15 或 Q31 格式）**。
* 在定点运算中，数字不能超过最大值（比如 16位有符号整数的最大值是 32767）。
* 标准 FFT 是做加法，$N$ 个点相加，数值会被放大 $N$ 倍。在定点 DSP 中，做标准 FFT **极其容易溢出**。
* 为了防止溢出，定点 DSP 工程师通常会在 FFT 的每一级蝶形运算中右移（相当于除以 2），或者在 FFT 结束后直接除以 $N$。
* 你看到的这段 Python 代码，极有可能是**从 Speex 的 C 语言定点代码直接 1:1 翻译过来的**，保留了其底层防溢出的数学特征。
**原因 C：为什么 IFFT 要乘以 N？**
因为你在 FFT 阶段主动除以了 $N$（把信号缩小了 $N$ 倍），而 NumPy 的 `ifft` 函数内部又很蠢地帮你除以了 $N$（再缩小 $N$ 倍）。
如果我不乘以 $N$，还原出来的时域信号就会变成原来的 $\frac{1}{N^2}$，声音就听不见了。**乘以 $N$ 是为了把 NumPy IFFT 自带的那个 $\frac{1}{N}$ 给抵消掉，并把之前 FFT 人为除掉的 $N$ 给补回来。**
—
## 8.3 总结
这段代码的本质是**改变了 FFT/IFFT 的缩放重心**：
* 标准库做法：能量缩放全压在 IFFT 上（FFT不缩放，IFFT除以N）。
* Speex做法：能量缩放全压在 FFT 上（FFT除以N，IFFT乘N抵消标准库的缩放并还原）。
在 Python 这种浮点数精度极高的环境里，这么做显得有些多余甚至令人困惑，但这正是**移植底层 C 语言音频处理算法时常见的“代码化石”**。它不仅保证了频域特征的一致性，也最大程度还原了原版 C 代码的数学等价性。
# 9 python代码仿真说明
## 9.1 main函数
“`
mic, sr = sf.read(f”./wav/rec.wav”)
ref, sr = sf.read(f”./wav/ref.wav”)
“`
调用后，sr=16000,mic和ref都是(604320，) float64的ndarray数据
“`
processor = MDF(sr, 256, 3072)
e, y = processor.main_loop(ref, mic)
“`
MDF函数中，frame_size=256，filter_length=3072
processor.main_loop返回的第一个值e是回声
## 9.2 MDF __init__函数
K=1是speaker数量，
C=1是mic数量
N=512是windows数量
M=3072/256=12是MDF算法中块的个数
spec_average=256/16000=0.016，即一帧数据的时长
beta0=2*256/16000=0.032，即两帧数据的时长
beta_max=0.5*256/16000，即0.5帧时长
e=np.zeros((N, C),) (512,1)
x=np.zeros((N, K),) (512,1)
input=np.zeros((self.frame_size, C),) (256,1)
y=np.zeros((N, C),) (512,1)
last_y=np.zeros((N, C),) (512,1)
Yf=np.zeros((self.frame_size+1, 1),) (257,1)
Rf=np.zeros((self.frame_size+1, 1),) (257,1)
Xf=np.zeros((self.frame_size+1, 1),) (257,1)
Yh=np.zeros((self.frame_size+1, 1),) (257,1)
Eh=np.zeros((self.frame_size+1, 1),) (257,1)
X=np.zeros((N, K, M+1), dtype=complex) (512,1,13)
Y=np.zeros((N, C), dtype=complex) (512,1)
E=np.zeros((N, C), dtype=complex) (512,1)
W=np.zeros((N, K, M, C), dtype=complex) (512,1,12,1)
foreground=np.zeros((N, K, M, C), dtype=complex) (512,1,12,1)
PHI=np.zeros((frame_size+1, 1),) (257,1)
power= np.zeros((frame_size+1, 1),) (257,1)
power_1=np.ones((frame_size+1, 1),) (257,1)
window=np.zeros((N, 1),) (512,1)
prop=np.zeros((M, 1),) (12,1)
wtmp=np.zeros((N, 1),) (512,1)
memX=np.zeros((K, 1),) (1,1)
memD=np.zeros((C, 1),) (1,1)
memE=np.zeros((C, 1),) (1,1)
notch_mem=np.zeros((2*C, 1),) (2,1)
## 9.3 main_loop函数
u_frame是(256，)，一帧ref数据，float64数据
d_frame是(256，)，一帧mic数据，float64数据
out是(256，1)， 一帧回声消除后的数据，float64数据
## 9.4 speex_echo_cancellation_mdf函数
N=512是windows数量
K=1是speaker数量，
C=1是mic数量
M=3072/256=12是MDF算法中块的个数
out是(256，1)， 一帧回声消除后的数据，float64数据
ss=0.02917
input是(256,1)
notch_mem是(2,1)
### 9.4.1 filter_dc_notch16函数**：
mic和out是（256,）
mem是(2,1)
out是用np.zeros_like(mic)重新申请的内存
### 9.4.2 前景滤波
x是(512，1)
X是(512，1, 13),复数
来源：image processing（微信号/QQ号：1439279），转载请注明出处，谢谢！
