简介一套将Speex AEC自适应回声消除MDF算法从C语言移植到Matlab的实现源码面向音频处理、语音通信方向的开发者与算法研究人员尤其适合需要快速理解回声消除原理并开展Matlab仿真验证的人群。压缩包共2个文件包含1个Matlab脚本.m和1个说明文档README整体体积仅7KB结构十分精简便于直接阅读、调试与运行验证。该移植版本基于多延迟滤波MDF自适应滤波器思想通过不断更新滤波器系数来追踪回声路径保留原C实现核心流程的同时充分利用Matlab的矩阵运算简化循环和数组操作降低了算法复现门槛README文件给出基本使用说明可帮助使用者快速验证回声消除效果。已有413人浏览学习对正在学习跨语言算法移植或Speex音频技术的开发者具有直接的参考价值。整体简练且聚焦是声学回声消除算法入门与二次开发的实用素材。1. 把 Speex AEC 的 MDF 搬到 Matlab到底在搬什么做回声消除产品迭代时手里往往只有一个跑得很稳的 C 版 Speex AEC效果能接受但你一旦想改双讲检测阈值、想看收敛过程中的频域系数、想对比不同分块长度下的 ERLE就得改 C、重新编译、再跑一遍嵌入式测试一次循环半小时起步。把 Speex AEC 的 MDFMultidelay Frequency Domain算法从 C 移植到 Matlab价值不是换个语言重复一遍而是让每个中间量都变成能画出来的变量参考谱、误差谱、每个分块的滤波器系数、功率估计全部摊开在 workspace 里。真正花时间的也不是自适应滤波的数学而是实谱 FFT 的打包方式、2N 点窗与 50% 重叠、参考信号延迟线以及 ERLE 驱动的步长控制。现在用 Codex 这类工具可以几分钟生成初版代码但 FFT 打包和时域约束这两处工具不会替你判断对错恰恰是整个移植里出错率最高的地方。这篇文章写给做算法验证、参数扫描和 AEC 训练数据生成的工程师按 C 版 Speex 的内部结构一步步拆。2. Speex AEC 的 MDF 算法骨架分块、频域更新与时域约束2.1 为什么是 MDF而不是时域 NLMS回声消除的建模很简单麦克风信号 d(n) 由近端语音 s(n)、远端参考 x(n) 经过回声路径 h 卷积后的回声、以及噪声 v(n) 组成目标是估计出 h 后把回声减掉。最朴素的时域 NLMS 每处理一个样本要 L 次乘加L 是滤波器长度在 8kHz 采样、L2048 时就是每秒 1600 万次乘加DSP 上能跑但代价不低而且时域 NLMS 对步长和参考信号的相关性非常敏感。MDF 的做法是把长度 L 的滤波器切成 P 段每段 N 点LP×N。每帧只处理 N 个新样本做一次 2N 点 FFT频域里做分段卷积再通过时域约束把循环卷积拉回线性卷积。于是自适应滤波的延迟从 L 降到了 N计算量也降了一个量级。单帧的复杂度从 NLMS 的 O(P·N²) 降到 O(N·log N P·N)分块数越多省得越明显。Speex 默认把 frame_size 取为 2568kHz或 51216kHzfilter_length 常见取 frame_size 的 8 到 16 倍也就是 P 在 8 到 16 之间。整个算法每帧的数据流可以用下面这段伪代码概括它同时也是后续 Matlab 移植的主线每帧输入 N 点参考 ref 和 N 点麦克风 mic 1. 把 ref 推进 2N 点缓冲加窗后做 FFT得到参考谱 X 2. 把 X 压入延迟线 Xbuf保留最近 P 帧的参考谱 3. 频域滤波Y sum_p W(p,:) .* Xbuf(p,:) 4. 麦克风信号同样加窗 FFT 得到 D误差谱 E D - Y 5. 自适应更新W mu * conj(Xbuf) .* E ./ SS 是参考功率估计 6. 时域约束对更新量和 W 都做 ifft把第 N 点之后的抽头清零再 fft第 5 步和第 6 步是 MDF 区别于普通频域 LMS 的关键。第 5 步里的 S 不是瞬时功率而是经过递归平滑的功率估计否则参考信号一有静音段步长就抖动第 6 步如果不做频域自适应滤波器会绕进循环卷积的误差里去ERLE 爬到十几 dB 就再也上不去。2.2 mdf 算法的三个核心公式把上一节的伪代码落到数学上一个 bin 一帧的更新可以写清楚。设 k 为频点序号p 为分块序号Xbuf_p(k) 表示延迟了 p 帧的参考谱% 滤波回声估计谱 Y(k) sum_p W_p(k) .* Xbuf_p(k) % 误差谱麦克风谱减去估计谱 E(k) D(k) - Y(k) % 更新参考谱的共轭乘误差谱除以平滑功率 W_p(k) mu .* conj(Xbuf_p(k)) .* E(k) ./ S(k) % 功率估计一阶递归平滑 S(k) beta .* S(k) (1 - beta) .* abs(X(k)).^2公式里最容易写错的是共轭的位置。误差信号 e(n) 和滤波器系数 w 的梯度关系是 X 的共轭乘 E一旦写成 conj(X).E 的反向形式即 X .conj(E)滤波器不但不收敛还会在几个帧内发散。除的是平滑功率 S(k) 而不是瞬时 |X(k)|²这决定了收敛速度和稳态失调之间的平衡beta 取 0.9 左右时二者都能兼顾。2.3 Speex 在 MDF 之上额外加的三样东西Speex 的 mdf.c 不是教科书 MDF 的直译它在核心更新之外加了三个工程上必需的机制。第一是时域约束每帧都做并且对滤波器系数和梯度都做第二是 ERLE 驱动的双讲检测用误差功率和麦克风功率的比值判断当前是否双讲双讲时冻结或降低步长防止近端语音把滤波器拉偏第三是泄漏系数和功率下限滤波器每帧乘一个略小于 1 的因子功率估计设置下限保证数值稳定性。移植时这三样一个都不能省但实现顺序可以调。下面这张表把 Speex 源码里的主要状态变量列出来作为后续映射的参照C 端字段speex_echo_cancel.c / mdf.c含义典型值/量级frame_size每帧样本数 N256 8kHz512 16kHzM / filter_length滤波器总长度 L8×N 到 16×NP分块数等于 L/N8 或 16N状态内FFT 点数等于 2×frame_size512 或 1024beta功率平滑系数0.9 左右step / mu自适应步长自适应数量级 0.10.5power_1功率倒数预计算随输入幅度标定erle / adapt_count双讲检测状态阈值随版本变化具体数值在不同 speexdsp 版本里有差异移植时以手头源码为准。3. 移植第一步把 C 语言结构体映射成 Matlab state跑通顶层循环3.1 speex_echo_cancel.c 的变量结构C 版的核心调用关系很简单speex_echo_state_init(frame_size, filter_length)分配一个 SpeexEchoState 结构体之后每帧调用一次speex_echo_cancellation(st, echo, ref, out)内部再调静态函数 mdf()。状态结构体里的字段表面上很多但按功能分组后只有四类时域缓冲frame、window、last_y、频谱缓冲X、D、Y、E、滤波器与延迟线W、Xbuf、自适应控制power、step、erle、adapt_count。Matlab 端对应关系用下面这张表记录后续所有代码都按这个映射写C 端作用Matlab 端st-frame_sizest-Mst-P帧长、滤波长度、分块数st.frame_sizest.filter_lengthst.num_partst-window2N 点分析/综合窗st.win用 sqrt(hann(2N)) 替代st-W[j][i]第 j 个分块第 i 个频点的滤波系数st.W(j,i)P×2N 复数矩阵st-X[k][i]参考频谱延迟线k0..Pst.Xbuf每帧把新谱压到第 1 行st-Dst-Yst-E麦克风谱、回声估计谱、误差谱局部变量 D、Y、Est-powerst-power_1参考功率估计及其倒数st.power直接用除法不用倒数st-stepst-erlest-adapt_count步长和双讲状态st.must.Sf/st.Sest.adapt_cntC 里那些spx_word32_t交错存储的复数数组到 Matlab 直接变成复数矩阵这一层的简化能让代码量少三分之一。移植时不要试图保留 C 的内存布局保留数据流的顺序就够了。3.2 顶层循环缓冲移位、加窗、FFT、输出重叠相加Speex 的时域处理是标准的 50% 重叠 WOLA加权重叠相加。每帧把之前 N 点与当前 N 点拼成 2N 缓冲加窗后 FFT输出端把误差谱 IFFT 回 2N 点时域前半段与上一帧遗留的后半段叠加得到当前帧输出。C 版窗型用的是自己定义的 spx_window它是 Hann 的一个变体移植时先用标准 sqrt-Hann 也能正常工作窗型只影响过渡带的泄漏不改变分块自适应滤波器的收敛性质。真正要严格对齐的是重叠比例和输出取哪半段这两个错了输出会有咔哒声。先写初始化函数对应speex_echo_state_initfunction st aec_mdf_init(frame_size, filter_length) % 对应 speex_echo_state_init(frame_size, filter_length) % filter_length 必须是 frame_size 的整数倍即分块数 P 为整数 if mod(filter_length, frame_size) ~ 0 error(filter_length 必须是 frame_size 的整数倍); end st.frame_size frame_size; st.filter_length filter_length; st.num_part filter_length / frame_size; % 分块数 P st.nfft 2 * frame_size; % FFT 点数对应 C 里 st-N st.mu 0.3; % 基础步长后面会被 ERLE 逻辑调制 st.beta 0.9; % 功率平滑系数 st.power_floor 1e-4; % 功率下限的相对系数 st.leak 1e-4; % 滤波器每帧泄漏量 st.freeze_len 20; % 连续多少帧 ERLE 差则冻结自适应 st.W zeros(st.num_part, st.nfft); % 频域滤波器复数 st.Xbuf zeros(st.num_part, st.nfft); % 参考谱延迟线 st.power ones(1, st.nfft) * 1e-2; % 功率估计初值 st.Sf 1; % 麦克风功率累积ERLE 用 st.Se 1; % 误差功率累积 st.adapt_cnt 0; % 双讲连续计数 st.ref_tail zeros(frame_size, 1); % 参考信号上一帧 st.mic_tail zeros(frame_size, 1); % 麦克风信号上一帧 st.out_tail zeros(frame_size, 1); % 输出重叠尾 st.win sqrt(hann(st.nfft, periodic)); % 50% 重叠 WOLA 用窗 endinit 里先把所有可变部分清零或置初值这是从 C 移植过来最容易忽略的一步。C 里 calloc 分配的内存天然是零Matlab 如果不显式初始化变量会在第一次赋值时动态扩容调试时很难区分是逻辑错误还是初始化遗漏。power 初值取 1e-2 而不是零是为了让第一帧更新时不出现除零。然后是单帧处理函数function [e, st] aec_mdf_block(ref, mic, st) % 处理一帧对应 C 里一次 speex_echo_cancellation() 调用 % ref / mic 都是 frame_size 的列向量返回误差输出 e N st.frame_size; nfft st.nfft; % 参考信号拼成 2N 缓冲加窗 FFT xbuf [st.ref_tail; ref(:)]; st.ref_tail ref(:); X fft(xbuf .* st.win, nfft); % 麦克风信号同样处理得到 D dbuf [st.mic_tail; mic(:)]; st.mic_tail mic(:); D fft(dbuf .* st.win, nfft); % 新参考谱压入延迟线第一行旧谱向后推 st.Xbuf [X; st.Xbuf(1:end-1, :)]; % 频域滤波P 个分块的加权和 Y sum(st.W .* st.Xbuf, 1); % 误差谱IFFT 回时域做 WOLA 输出 E D - Y; e_full real(ifft(E, nfft)); e e_full(1:N) st.out_tail; st.out_tail e_full(N1:nfft); % 自适应更新见第 4 章 st mdf_update(st, X, D, E); end这一段对应 C 版里从缓冲拼接到调用 mdf() 之前的所有代码。注意延迟线的推进方式[X; st.Xbuf(1:end-1,:)]把最新谱放到第 1 行第 P 行的谱被丢掉正好对应 P 个分块的延迟范围。麦克风和参考共用同一个窗保证频域里 D 和 Y 的加权方式一致。3.3 FFT 打包差异C 的 N1 点实谱 vs Matlab 全谱移植时一个绕不开的坑是 FFT 的数据排布。Speex 用的是 kiss_fft 的实输入接口2N 点实数 FFT 输出被压缩成 N1 个复数第 0 个位置放直流实数第 1 个位置放奈奎斯特频率实数从第 2 个位置开始依次是 Re(bin1)、Im(bin1)、Re(bin2)、Im(bin2) 的交错排列。整个 mdf.c 里的数组长度都是 N1而不是 2N。Matlab 里有两种对应方案。方案 A 是直接用全谱fft 得到 2N 点复数负频率部分不做任何特殊处理。因为输入是实数X、D、E、W 始终满足共轭对称时域约束又进一步强制了这种对称性所以全谱更新不会引入额外的自由度结果和 C 版等价。方案 B 是精确复刻 N1 点打包仅当你需要和 C 版做逐 bit 对比时才值得写。做算法验证用方案 A 就够第 5 章的对齐测试也用方案 A 的中间变量。另外 C 版默认是 float 单精度Matlab 是 double数值对比时两边误差在 1e-5 量级属于正常不用怀疑移植错了。4. Matlab 端 mdf 更新核心功率归一化、梯度约束和双讲步长4.1 mdf_update 主更新循环第 3 章的 block 函数把每帧的频谱算好真正的 MDF 自适应更新集中在 mdf_update 里。它对应 C 版 mdf() 中从功率估计到 W 更新的整段逻辑四件事按顺序做更新功率估计、判双讲、算梯度并约束、更新并约束滤波器系数。function st mdf_update(st, X, D, E) % MDF 自适应更新核心对应 C 版 mdf() 的更新段 N st.frame_size; nfft st.nfft; P st.num_part; % 1) 参考功率递归估计。下限取相对值跟随输入幅度标定 % 这样 16bit 定点标定和浮点 -1~1 标定都能用同一份代码。 Xpow abs(X).^2; st.power st.beta * st.power (1 - st.beta) * Xpow; st.power max(st.power, st.power_floor * mean(st.power)); % 2) ERLE 驱动的双讲检测误差功率相对于麦克风功率的下降量 st.Sf 0.95 * st.Sf 0.05 * sum(abs(D).^2); st.Se 0.95 * st.Se 0.05 * sum(abs(E).^2); erle_db 10 * log10((st.Sf eps) / (st.Se eps)); if erle_db 12 st.adapt_cnt 0; % 对消效果好认为没有双讲 else st.adapt_cnt st.adapt_cnt 1; end % 参考能量过低时同样冻结避免静音段噪声被学进去 ref_energy mean(Xpow); if ref_energy st.power_floor * mean(st.power) eff_mu 0; elseif st.adapt_cnt st.freeze_len eff_mu 0; % 连续多帧对消不掉判为双讲 else eff_mu st.mu; end % 3) 梯度参考谱的共轭乘误差谱。方向反了会直接发散。 grad conj(st.Xbuf) .* E; % P x 2N % 4) 梯度时域约束只保留前 N 个抽头 grad_t real(ifft(grad, nfft, 2)); grad_t(:, N1:nfft) 0; grad fft(grad_t, nfft, 2); % 5) 功率归一化更新st.power 是 1x2N自动广播到 P 行 st.W st.W eff_mu * grad ./ st.power; % 6) 滤波器系数同样做时域约束并加泄漏 Wt real(ifft(st.W, nfft, 2)); Wt(:, N1:nfft) 0; st.W (1 - st.leak) * fft(Wt, nfft, 2); end代码里有三个细节要说明。功率下限用st.power_floor * mean(st.power)这种相对值而不是绝对常数是因为 AEC 的输入标定在不同平台上差异很大C 版如果跑 16bit 定点信号幅度量级是几千到几万浮点版则归一化到 ±1绝对下限两边没法通用。gradient 的约束和 W 的约束分开做二者缺一不可只约束 W 不约束梯度更新量里会带进循环卷积误差稳态 ERLE 会低几个 dB。real(ifft())把时域滤波器的虚部直接丢掉等价于把频域系数投影到实滤波器空间这是频域自适应滤波的标准操作。4.2 双讲检测阈值怎么调双讲检测是 AEC 里最影响主观听感的环节。近端有人说话时误差谱里同时有回声残差和近端语音此时继续更新滤波器近端语音会被当成回声路径学进去导致近端语音被吃掉。上面代码是简化版逻辑erle_db 12认为当前没有双讲连续 20 帧不满足就冻结。12dB 和 20 帧两个阈值对应不同场景会议场景近端说话密集阈值要抬到 15dB 以上、冻结帧数增加到 50纯音乐播放场景近端很少说话阈值可以降到 8dB让滤波器在双讲间隙更快恢复跟踪。Speex 原版不是这种硬切换它用连续的步长调整ERLE 高时步长缓慢增大ERLE 低时步长快速减小并且参考静音时直接置零。硬切换在边界处会让滤波器步长抖动表现为回声忽大忽小但作为移植的第一版硬切换更容易调试曲线平稳后再改成连续调整。还有一点上面的 eff_mu 冻结的是整个滤波器矩阵更精细的做法是按频点或按分块冻结低频段双讲能量集中高频段可以继续更新。4.3 与 C 版的一处重要差异分块能量 spread 归一化speexdsp 1.2 的 mdf.c 里有一块教科书 MDF 没有的逻辑它除了维护当前帧参考谱的功率估计还统计 P 个分块各自的能量分布用 spread 参数对每个分块的步长做加权。原因很直观回声路径的脉冲响应能量通常集中在前几个分块如果所有分块用同一个步长尾部那些能量很小的分块会被噪声持续扰动稳态失调变大。简化版用同一份 st.power 除所有分块白噪声参考下收敛曲线已经足够平滑等到回声路径突变或双讲后的恢复阶段出现抖动再回来补 spread 归一化。4.4 起始参数表下面的参数组合适合 8kHz、frame_size256、filter_length2048 的配置作为第一轮验证的起点参数建议初值调整方向st.mu0.3收敛慢则加大到 0.5稳态抖动大则减到 0.1st.beta0.9回声路径快速变化时降到 0.8st.power_floor1e-4ERLE 天花板偏低时试着降到 1e-5st.leak1e-4长时间运行 W 范数增长时加大到 1e-3freeze_len20近端语音密集场景加到 505. 用合成回声路径量化 ERLE把移植偏差压到最低5.1 测试管线和 ERLE 曲线移植完成后第一件事不是接真实音频而是用合成回声路径做定量验证。合成路径的好处是回声路径 h 已知ERLE 的理论上限可预估任何偏差都能定位。测试流程生成指数衰减的随机脉冲响应把白噪声参考和它卷积得到纯回声加一段系统延迟模拟设备延迟然后跑 AEC 看 ERLE。% test_aec_mdf.m fs 8000; frame_size 256; % 32ms 帧 filter_length frame_size * 8; % 合成回声路径指数衰减随机脉冲长度远小于 filter_length t (0:1023) / fs; h exp(-t / 40e-3) .* randn(1024, 1); h h / norm(h) * 0.5; rng(1); ref randn(fs * 5, 1); % 白噪声参考 mic [zeros(64,1); filter(h,1,ref)]; % 回声加 64 样本延迟 mic mic(1:length(ref)); % 逐帧跑 AEC st0 aec_mdf_init(frame_size, filter_length); e zeros(size(mic)); st st0; nblk floor(length(mic) / frame_size); for n 1:nblk idx (n-1)*frame_size1 : n*frame_size; [e(idx), st] aec_mdf_block(ref(idx), mic(idx), st); end % 逐帧 ERLE50 帧滑动平均 M 50; nblk floor(length(mic) / frame_size); erle zeros(nblk, 1); for n M1:nblk idx (n-1)*frame_size1 : n*frame_size; p_mic mean(mic(idx).^2) eps; p_err mean(e(idx).^2) eps; erle(n) 10*log10(p_mic / p_err); end plot(erle); ylabel(ERLE (dB)); xlabel(帧序号);白噪声参考下ERLE 应在 200 帧内从 0 爬到 25dB 以上稳态在 3040dB 之间。如果参考换成语音稳态 ERLE 会低 610dB这是 LMS 类算法在白谱输入下收敛最优的正常表现不要拿语音的 ERLE 和 C 版白噪声的结果直接比。如果把回声路径延长到接近或超过 filter_lengthERLE 会出现明显天花板这是滤波器长度不够的物理限制不是移植错误。5.2 常见移植偏差诊断跑完测试后对照下面这张表排查症状常见原因检查点ERLE 始终在 03dB曲线不动更新方向错了conj 位置写反打印更新前后 W 的范数应缓慢增长ERLE 爬到 15dB 后缓慢下降缺时域约束或泄漏过大检查 grad 和 W 是否都做了 ifft 清尾ERLE 有天花板比如卡在 20dB回声路径超过 filter_length或功率下限太高缩短 h 到 filter_length 的一半再试输出有周期性咔哒声WOLA 窗不满足 COLA或输出取错半段确认 sqrt(hann) 加 50% 重叠确认 e 取前 N 点与 C 版中间量对不上FFT 缩放或 DC/Nyquist 打包不一致用相同输入 dump 第一帧 X、Y差 2N 倍就是缩放问题5.3 用脉冲响应可视化做最终确认ERLE 曲线只能说明对消效果不能确认滤波器学到的是正确的回声路径。把 st.W 每帧拉回时域能直观确认三个环节。做法是把 W 按行 IFFT每行取前 N 点按分块顺序拼起来Wt zeros(filter_length, 1); for p 1:st.num_part wt real(ifft(st.W(p,:))); Wt((p-1)*frame_size1 : p*frame_size) wt(1:frame_size); end stem(Wt); hold on; stem([zeros(64,1); h], r);学到的脉冲响应应该在延迟 64 样本处与真实 h 对齐分块边界处没有明显的阶跃跳变。如果峰值位置正确但增益偏小检查第一步的 power_floor 相对值是否过高如果峰值位置偏移检查系统延迟和 Xbuf 的索引顺序这是参考信号与麦克风对齐的问题。做这一步时同时在 VSCode 里给 C 版加几行 fprintf把同样输入下的 W dump 出来两边画在同一张图里是定位移植偏差最快的方式。数值差在 1e-5 量级是单双精度差异差到 1e-2 就要回到 FFT 打包和窗函数上找原因。本文还有配套的精品资源点击获取