Speex回声消除原理深度解析
这里假设读者具有自适应滤波器的基础知识。Speex的AEC是以NLMS为基础，用MDF频域实现，最终推导出最优步长估计：残余回声与误差之比。最优步长等于残余回声方差与误差信号方差之比，这个结论可以记下，下面会用到的。
对于长度为N的NLMS滤波器，误差信号定义为期望信号与估计信号之差，表示如下：
\[e(n) = d(n) - \hat y(n) = d(n) - \sum\limits_{k = 0}^{N - 1} {{{\hat w}_k}(n)x(n - k)} \]
则，滤波器的系数更新方程为：
\[{\hat w_k}(n + 1) = {\hat w_k}(n) + \mu \frac{{e(n){x^*}(n - k)}}{{\sum\nolimits_{i = 0}^{N - 1} {|x(n - i){|^2}} }} = {\hat w_k}(n) + \mu \frac{{(d(n) - \sum\nolimits_i {{{\hat w}_i}(n)x(n - i)} ){x^*}(n - k)}}{{\sum\nolimits_{i = 0}^{N - 1} {|x(n - i){|^2}} }}\]
设滤波器的系数误差为：
\[{\delta _k}(n) = {\hat w_k}(n) - {w_k}(n)\]
且期望信号为本地（近端）语音+残余回声
\[d(n) = v(n) + \sum\nolimits_k {{w_k}(n)x(n - k)} \]
则滤波器的系数更新方程可以重写为
\[{\delta _k}(n + 1) = {\delta _k}(n) + \mu \frac{{(v(n) - \sum\nolimits_i {{\delta _i}(n)x(n - i)} ){x^*}(n - k)}}{{\sum\nolimits_{i = 0}^{N - 1} {|x(n - i){|^2}} }}\]
如果每个时刻的失调定义为：
\[\Lambda (n) = \sum\nolimits_k {\delta _k^*(n){\delta _k}(n)} \]
那么，在每一步的迭代中，滤波器的失调可表示如下：
\[\Lambda (n + 1) = \sum\limits_{k = 0}^{N - 1} {|{\delta _k}(n) + \mu \frac{{(v(n) - \sum\nolimits_i {{\delta _i}(n)x(n - i)} ){x^*}(n - k)}}{{\sum\nolimits_{i = 0}^{N - 1} {|x(n - i){|^2}} }}{|^2}} \]
假设远端信号与近端信号为白噪声，且不相关。
\[\sigma _v^2 = E\{ |v(n){|^2}\} \]
为近端语音信号的方差，则失调的更新方程为
\[E\{ \Lambda (n + 1)|\Lambda (n),x(n)\} = \Lambda (n)\left[ {1 - \frac{{2\mu }}{N} + \frac{{{\mu ^2}}}{N} + \frac{{2{\mu ^2}\sigma _v^2}}{{\Lambda (n)\sum\nolimits_{i = 0}^{N - 1} {|x(n - i){|^2}} }}} \right]\]
这里失调函数
\[E\{ \Lambda (n + 1)|\Lambda (n),x(n)\} \]
为凸函数，对它关于步长求导，并置导数为0，可得：
\[\frac{{\partial E\{ \Lambda (n + 1)\} }}{{\partial \mu }} = \frac{{ - 2}}{N} + \frac{{2\mu }}{N} + \frac{{2\mu \sigma _v^2}}{{\Lambda (n)\sum\nolimits_{i = 0}^{N - 1} {|x(n - i){|^2}} }} = 0\]
最终推出最优步长为：
\[{\mu _{opt}}(n) = \frac{1}{{1 + \frac{{\sigma _v^2}}{{\Lambda (n)/N\sum\nolimits_{i = 0}^{N - 1} {|x(n - i){|^2}} }}}}\]
大家别看最下面的那个分母
\[\Lambda (n)/N\sum\nolimits_{i = 0}^{N - 1} {|x(n - i){|^2}} \]
式子挺长，其实意义很明确，可以近似理解为残余回声的方差，于是输出信号的方差为：近端语音的方差+残余回声的方差，用式子表示如下
\[\sigma _e^2(n) = \sigma _v^2(n) + \sigma _r^2(n)\]
最终，导出最优步长：
\[{\mu _{opt}}(n) = \frac{1}{{1 + \frac{{\sigma _v^2}}{{\sigma _r^2(n)}}}} = \frac{1}{{\frac{{\sigma _r^2(n) + \sigma _v^2}}{{\sigma _r^2(n)}}}} \approx \frac{{\sigma _r^2(n)}}{{\sigma _e^2(n)}}\]
\[{\mu _{opt}}(n) = \min \left( {\frac{{\hat \sigma _r^2(n)}}{{\hat \sigma _e^2(n)}},1} \right)\]
上面的分析是在时域，基于NLMS，可以看到：最优步长等于残余回声方差与误差信号方差之比。其中误差的方差比较好求，残余回声的方差比较难求。下面我们看下上面的结论在频域中如何解决，Speex中在频域的自适应算法为：MDF（multidelay block frequency domain）自适应滤波。
在频域中，设k为频率索引，字母（ell）为帧索引，上面的结论转换到频域，结果如下：
\[{\mu _{opt}}(k,\ell ) \approx \frac{{\sigma _r^2(k,\ell )}}{{\sigma _e^2(k,\ell )}}\]
那么，在频域如何求残余回声的方差呢，我们可以定义一个泄露系数，表示回声相对于估计回声信号的泄露程度，这时残余回声表示为
\[\sigma _r^2(k,\ell ){\rm{ = }}\hat \eta (\ell )\hat \sigma _{\hat Y}^2(k,\ell )\]
根据泄露系数求出残余回声，就可以得到最优步长
\[{\mu _{opt}}(n) = \min \left( {\hat \eta (\ell )\frac{{|\hat Y(k,\ell ){|^2}}}{{|E(k,\ell ){|^2}}},{\mu _{\max }}} \right)\]
也就是说，根据泄露系数，可以估计出远端信号的残余回声，进而可以得到最优步长，那么，带来另一个问题，这里的泄露系数如何估计呢？确定泄露系数的过程，其实就是一元线性回归分析中确定回归系数的过程，具体可以看下回归分析的内容。
\[\hat \eta (\ell ) = \frac{{\sum\nolimits_k {{R_{EY}}(k,\ell )} }}{{\sum\nolimits_k {{R_{YY}}(k,\ell )} }}\]
\[{R_{EY}}(k,\ell ) = (1 - \beta (\ell )){R_{EY}}(k,\ell ) + \beta (\ell ){P_Y}(k){P_E}(k)\]
\[{R_{YY}}(k,\ell ) = (1 - \beta (\ell )){R_{YY}}(k,\ell ) + \beta (\ell ){P_Y}(k){P_Y}(k)\]
\[\beta (\ell ) = {\beta _0}\min (\frac{{\hat \sigma _Y^2(\ell )}}{{\hat \sigma _e^2(\ell )}},1)\]
这里， 是通过递归平均处理方法得到每个频点的自相关、输入信号与误差信号的互相关。最终得到泄露系数，具体实现可以参考speex 的代码实现，相关参数可以参考后面给出来参考论文。
Speex的回声消除原理已经分析完了，最终得出结论是：只有改与泄露系数相关部分的代码，才是对效果影响最大的地方，因为根据泄露系数，最终会估计出滤波器的最优步长。
参考论文：On Adjusting the Learning Rate in Frequency Domain Echo Cancellation With Double-Talk
