speex_mdf
2026-04-03 | 分类:未分类 | 评论:0人 | 浏览:10次
- 正文内容
- 我来说两句:(已有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 |
## 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 变学习率机制
### 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;
“`
### 6.4.2 自适应学习率
**学习率调整公式**：
对于每个频率点 $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$：残余误差比
**代码实现**：
“`matlab
for i=1:st.frame_size+1
r = st.leak_estimate*st.Yf(i);
e = st.Rf(i)+1;
if (r>.5*e)
r = .5*e;
end
r = 0.7*r + 0.3*(RER*e);
st.power_1(i) = (r/(e*st.power(i)+10));
end
“`
### 6.4.3 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;
“`
### 6.4.4 双讲处理
**无显式双讲检测**，而是通过连续调整学习率实现鲁棒性。
**双讲场景分析**：
| 场景 | $S_{dd}$ | $S_{ee}$ | 学习率 |
|——|———-|———-|——–|
| 仅远端 | 正常 | 小 | 正常 |
| 双讲 | 大 | 大 | 自动减小 |
| 仅近端 | 大 | 大 | 自动减小 |
**机制**：
$$\alpha = \frac{\beta_0 \cdot S_{yy}}{S_{ee}}$$
双讲时 $S_{ee}$ 增大，$\alpha$ 减小，学习率降低。
## 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
“`
**切换逻辑分析**：
**条件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` 的尾部，用于下一帧的误差）。
如果不更新前景，则检查是否需要将背景重置为前景（即前景优于背景）。若满足条件（负向能量差异大），则重置背景为前景，并更新误差和统计量。
### 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 代码的数学等价性。
