第五章 自动回声消除(AEC)原理

回声消除是为了消除机器自身发出的声音，不影响外界传递过去的声音，如在通话场景中，来自远端的经过扬声器放出来的声音被消掉，否则经过麦克风采集和近端说话人信号混在一起被传递到远端，则会出现远端的人听到自己说的话，这种体验不好。在远场音箱交互场景中，这自动回声消除的难度比通话场景要复杂一下，主要的原因是近端说话的语音是淹没在大功率喇叭且存在低音单元之中的。

回声消除原理

回声消除的基本原理是用一个自适应滤波器对未知的回声信道$\omega$进行参数辨识，根据扬声器信号与其产生回声信号相关性为基础，建立远端信号模型，模拟回声路径，通过自适应算法调整，使其冲击响应和真实回声路径相逼近。然后将麦克风接收到的信号减去估计值，即可实现回声消除功能。

$echo = x* \omega \tag {1}$
$d = s + echo \tag {2}$
$\hat y = x * \hat \omega \tag{3}$
$e = d - \hat y \tag {4}$

式中，$\omega$是回声通道的时域冲击响应函数；$x$是远端语音；echo是所得回声；s是近端说话人语音，d为麦克风采集到的信号；$\hat y$为对回声信号的估计值；e为误差。在电话、视频会议中，x指通话另一端的语音信号，而在机器语音识别场景，则指机器自身扬声器发出的声音。

为了消除延迟较大的回声，FIR滤波器的阶数要尽量的大。时域计算诸多不便，使用频域分块自适应滤波算法，提高实时性。涉及到的信号处理术语：

• FFT/IFFT
• 循环卷积和线性卷积的关系;重叠保留法
• 功率谱密度
• 互相关
• LMS/NLMS自适应算法
• 频域分块
• 宽带信号处理

维纳滤波

均方误差(MSE，Mean Square Error)，对于离散时间系统，可定义期待响应$d_k$为一个希望自适应系统的输出$y_k$与之相接近的信号，k为采样时刻。

根据图2， 可以求得误差信号:

$\varepsilon_k = d_k - y_k \tag {5}$

自适应线性组合器输出：

$y_k = W_k^TX_k \tag {6}$

其中：

$X_k = \begin{bmatrix} x_{0k}&x_{1k}&\cdots & x_{Lk} \end{bmatrix}, W_k = \begin{bmatrix} w_{0k}&w_{1k}&\cdots & w_{Lk} \end{bmatrix}$

分别为自适应系统在k时刻的输入信号向量和权向量，系统的均方误差为：

$E(|\varepsilon_k|^2) = E[(d_k - y_k)^{*}(d_k-y_k)] = E(|d_k|^2)+\mathbf W_k^H E[\mathbf X_k^* \mathbf X_k^*]\mathbf W_k - 2Re{\mathbf W_k^T E[d_k^* \mathbf X_k]} \tag {7}$

令：

$\mathbf R = E[\mathbf X \mathbf X^T] =\begin{bmatrix} E[x_0x_0]& E[x_0x_1] & \cdots & E[x_0x_L] \\ E[x_1x_0]& E[x_1x_1] & \cdots & E[x_1x_L]\\ \vdots& \vdots & \ddots & \vdots \\ E[x_Lx_0]& E[x_Lx_1] & \cdots & E[x_Lx_L] \end{bmatrix} \tag {8}$

定义期待响应和输入信号之间的互相关向量为：

$\mathbf P = E[d^*\mathbf X] = \begin{bmatrix} d^*x_0\\ \vdots\\ d^*x_L \end{bmatrix} \tag {9}$

将式9简化成下式：

$\xi(\mathbf w) = E(|d_k|^2)+\mathbf W_k^H \mathbf R \mathbf W_k - 2Re{\mathbf W_k^T \mathbf P} \tag{10}$

理想情况下$E(|\varepsilon|^2)$等于零，是估计值等于观测值，误差越小， 对式10求偏导数，得:

$\nabla = \frac{\partial}{\partial \mathbf W}[\xi(\mathbf W)] = 2 \mathbf R \mathbf W -2 \mathbf P^* \tag {11}$

最佳权向量处的梯度等于0，于是：

$\nabla = 2 \mathbf R \mathbf W_{opt} - 2 \mathbf P^* = 0 \tag {12}$

最小均放误差输出情况下的最佳权向量$\mathbf W_{opt}$满足维纳-霍夫方程：

$\mathbf W_{opt} = \mathbf R^{-1} \mathbf P^* \tag{13}$

值得指出的是MSE准则应用的地方还是很多，在后面阵列波束形成中的一个实例就会用到这一准则。

LMS算法

$\varepsilon_k = d_k - \mathbf X_k^T \mathbf W_k \tag {14}$

式中，$\mathbf X_k$为输入样本向量,使用单次采样数据$|\varepsilon_k|^2$来代替均方误差$\xi_k$，这样梯度估计可表示为如下形式:

$\hat \nabla_k = \frac{\partial}{\partial \mathbf W_k}|\varepsilon_k|^2 = \frac{\partial}{\partial \mathbf w_k}[|d_k|^2 + \mathbf W_k^H \mathbf X_k^* \mathbf X_k^T \mathbf W_k - 2d \mathbf X_k^*] = -2\varepsilon_k \mathbf X_k^* \tag {15}$

基于最速下降法的权向量迭代如下:

$\mathbf W_{k+1} = \mathbf W - \mu \hat \nabla_k = \mathbf W_k + 2 \mu \varepsilon_k \mathbf X_k^* \tag {16}$

其中$\mu$是步长因子，$0 < \mu < \frac {1}{\lambda_{max}}$， $\lambda_{max}$是$\mathbf R_{xx}$的最大特征值。$\mathbf W(k)$收敛于由比值决定，该比值$d=\frac {\lambda_{max}}{\lambda_{min}}$决定,他的比值叫做谱动态范围。大的d值意味着较长的时间才能收敛到最佳权值。
该算法用在语音增强的加性噪声消除功能上时，其工程实践并不完全按照式14意义来实现。

自适应LMS算法在主动降噪中的应用

算法的原理，假设$e(n)$是图一中麦克风采集到的消声后的信号，$d(n)$是真实的环境噪声信号，$x(n)$被称为参考信号，这可以通过麦克（非耳塞上的麦克）来获取， 而$w(n)x(n)$是扬声器发出的信号，用$w(n)x(n)$这一加权值来预测$d(n)$的值，如果预测的准，那么$e(n)$就等于0，即麦克风的震膜没有震动，人也就听不到噪声了。这个算法称为LMS算法。

在语音增强中，其目的是获得纯净的语音信号$s(n)$即上图中的最后输出信号，输入信号有两种，一种是带噪的语音信号$s(n)+x(n)$，另一种是只有噪声的输入$x(n)$，在没有人说话的情况下的输入信号，就仅仅是噪声输入。 这里要使得噪声估计$\hat x(n)$非常接近$x(n)$，对于$\varepsilon = \sum \limits_n s(n)+x(n) -\hat x(n)$，这时如果$\varepsilon^2$最小，则可以估计出的$\varepsilon$最接近$s(n)$。该过程可以概述如下:

1. 首先获取到噪声输入$x(n)$,并存储下来，以64或者128点为总长度不断刷新存储噪声输入$x(n)+s(n)$。
2. 采集带噪声的语音信号。
3. 用采集带噪语音信号减去估计到的噪声信号s(n)+x(n)-\hat x(n)。
4. 用3的输出作为误差，调节噪声权向量W。
   MATLAB实现具体包括如下三个部分：

% Loop over input vector
for ii = 1:length(signal_with_noise)
% Update buffer
%输入噪声估计
noise_buf = obj.update_buf(noise_buf, noise(ii));
% Filter this sample with current coefficient values
%通过权向量估计
filter_output = obj.data_filter(coefs, noise_buf);
% Compute error,相减得到
err = signal_with_noise(ii) - filter_output; 
% Update coefficients
coefs = obj.update_coefs(coefs, noise_buf, obj.filter_params.step_size, 
%用调节权向量
obj.filter_params.leakage, err); 
% Build output vector
%存储输入信号的估计值
dout(ii) = err;

NLMS算法

当输入信号幅值较大时，梯度噪声方法会遇到问题，使得能量低的信号算法收敛速度比较慢，将输入信号按照自身的平均能量进行归一化处理，既可以得到归一化的LMS算法(NLMS)，鲁棒性比较好。设输入带噪信号可以表示为：$x(n)$，其迭代算法NLMS公式如下：

$\mathbf W_{n+1} = \mathbf W_n - \mu \hat \nabla_n = \mathbf W_n + \frac{\mu}{N} \frac{e(n)\mathbf x(n)}{\hat \delta_x^2(n)}\tag {17}$

其中，$\hat \delta_x^2(n) = \frac{1}{N} \sum \limits_{n=0}^{N-1}x^2(n-k)$，其中，N是噪声消除器和回波抵消器的长度。(常取512,1024); $\mu$是步长因子。当$\hat \delta_x^2(k)$较小时, $\frac{\mu}{\hat \delta_x^2(n)}$的值可能较大，这是迭代算法变成如下形式：

$\mathbf W_{n+1} = \mathbf W_n - \mu \hat \nabla_n = \mathbf W_n + \frac{\mu}{N} \frac{e(n) \mathbf x(n)}{\sigma + \hat \sigma_x^2(n)} \tag{18}$

其计算过程如下:
参数：M是滤波器抽头系数(阶数)，$\mu$是自适应常数，$0 < \mu < 2 \frac{E[|\mathbf x(n)|^2] E[|\mathbf \varepsilon(n)|^2]}{E[|\mathbf e(n)|^2]}$， 其中$E[|\mathbf \varepsilon(n)|^2] = E[|\mathbf W_{opt} - \hat W(n)|^2]$,是权向量均方偏差，$\mathbf W_{opt}$是最优维纳解，$\hat W(n)$是第n次迭代中得到的估计值。$E[|\mathbf x(n)|^2]$是带噪输入语音信号的功率，$E[|\mathbf e(n)|^2]$是误差信号功率。

计算过程

1.初始化过程
如果知道抽头权向量$\hat W(n)$的先验知识,则用其来初始化$\hat W(0)$,否则令$\hat W(0) = \vec 0$。
2.数据按帧处理
A)给定的$\mathbf x(n)$第n时刻$MX1$抽头输入向量,$d(n)$是第n时刻的期望响应.
B)要计算的：$\hat W(n+1)$是第n+1步抽头权向量估计
3.计算
对$n=0,1,2,...$计算：

$e(n) = d(n) - \mathbf {\hat W^H(n)} \mathbf X(n) \tag {19}$
$\mathbf W_{n+1} = \mathbf W_{n} - \mu \hat \nabla_n = \mathbf W_n + \frac{\mu}{N} \frac{e(n)\mathbf X(n)}{\hat \sigma_x^2(n)} \tag {20}$

此外还有NLMS变种的各种方法,如SE-LMS(signed-error)LMS, SD-LMS(signal-dependent LMS)，LLMS(Leaky LMS), LNLMS(Leaky NLMS).

块自适应滤波

对参考信号x分段并做FFT变换，分别对各段数据做频域滤波，累加后做FFT反变换，并只取后L(L是原始信号的分段后的长度)点为有效的线性卷积结果，得到的是估计信号，将估计信号从回声信号中去除，得残差信号。计算子带步长，调整各段滤波器系数。这一过程表示如下图。

设n时刻输入序列$x(n)$如下:

$$\mathbf X(n) = \begin{bmatrix} x(n)& x(n-1) & \cdots & x(n-M+1) \end{bmatrix}^T \tag{21}$$

对应于长度为M的FIR滤波器在n时刻的抽头权向量为：

$$\hat {\mathbf W}(n) = \begin{bmatrix} \hat {\mathbf w}_0(n)& \hat {\mathbf w}_1(n-1) & \cdots & \hat {\mathbf w}_{M-1}(n) \end{bmatrix}^T \tag{22}$$

根据FIR滤波器原理:

$y(n) = x(n) \hat w_0(n) + x(n-1) \hat w_1(n) + \cdots _ x(n-M+1)\hat w_{M-1}(n) \tag {23}$

用向量表示如下:

$y(n) = \mathbf X(n)^T \hat {\mathbf W}(n) \tag {24}$

下面对$x(n)$进行分块，设k表示块下标，它与原始样值时间n的关系为：
$n = kL+i, i=0,1,\cdots, L-1; k = 1,2,\cdots$
其中L是块的长度。 第k块的数据为${\mathbf X(kL+i)}_{i=0}^{L-1}$,其矩阵表示形式如下:

$\mathbf A^T(k) = \begin{bmatrix} x(kL)& x(kL+1) & \cdots & x(kL+L-1) \end{bmatrix} \tag {25}$

将滤波器对输入块$A(k)$的响应表示如下:

$y(kL+i) = \hat {\mathbf W}^T(k)\mathbf A(k) = \sum \limits_{j=0}^{M-1}\hat w_j(k)x(kL+i-j), i=0, 1,\cdots, L-1 \tag {26}$

设$d(kL+i)$表示期望信号，误差信号表示如下:

$e(kL+i) = d(kL+i) - y(kL+i) \tag {27}$

考虑滤波器长度M=3，块长度L=3,其三个相邻的数据块是k-1，k， k+1，则第k-1块滤波结果如下：

$$\begin{vmatrix} y(3k-3)\\ y(3k-2)\\ y(3k-1) \end{vmatrix} = \begin{vmatrix} x(3k-3) & x(3k-4) & x(3k-5)\\ x(3k-2) & x(3k-3) & x(3k-4)\\ x(3k-1) & x(3k-2) & x(3k-3) \end{vmatrix} \begin{vmatrix} w_0(k-1)\\ w_1(k-1)\\ w_2(k-1) \end{vmatrix} \tag{28}$$

则第k块滤波结果如下：

$$\begin{vmatrix} y(3k)\\ y(3k+1)\\ y(3k+2) \end{vmatrix} = \begin{vmatrix} x(3k) & x(3k-1) & x(3k-2)\\ x(3k+1) & x(3k) & x(3k-1)\\ x(3k+2) & x(3k+1) & x(3k) \end{vmatrix} \begin{vmatrix} w_0(k)\\ w_1(k)\\ w_2(k) \end{vmatrix} \tag{29}$$

则第k+1块滤波结果如下：

$$\begin{vmatrix} y(3k+3)\\ y(3k+4)\\ y(3k+5) \end{vmatrix} = \begin{vmatrix} x(3k+3) & x(3k+2) & x(3k+1)\\ x(3k+4) & x(3k+3) & x(3k+2)\\ x(3k+5) & x(3k+4) & x(3k+3) \end{vmatrix} \begin{vmatrix} w_0(k+1)\\ w_1(k+1)\\ w_2(k+1) \end{vmatrix} \tag{29}$$

上面的数据矩阵是托伯利兹矩阵，主对角线元素都相同。

权向量调整公式如下

（权向量的调整）=（步长参数）*（抽头输入向量）*（误差信号）
因为在块LMS算法中误差信号随抽样速率而变，对于每一个数据块，我们有不同的用于自适应过程的误差信号值。因此，每一个块的抽头权向量更新公式如下：

$\hat w(k+1) = \hat w(k) + \mu \sum \limits_{i=0}^{L-1}x(kL+i)e(kL+i) \tag {30}$

其梯度向量的估计如下:

$\hat \nabla(k) = - \frac{2}{L}\sum \limits_{i=0}^{L-1}x(kL+i)e(kL+i) \tag {31}$

$\hat \nabla(k)$的无偏估计如下:

$\hat w(k+1) = \hat w(k) - \frac{1}{2} \mu_B \hat \nabla(k) \tag {32}$

块LMS算法的收敛性

由于时间平均的缘故，它具有估计精度随快长度增加而大幅提高的特性。然而，长度的增加会导致其收敛速度进一步减慢。后文的快速LMS算法加速了这一过程。

平均时长数

$\tau_{mse,av} = \frac{1}{2\mu_B \lambda_{av}} \tag{33}$

其中,$\tau_{av}$是输入自相关矩阵$\mathbf R = E[\mathbf X(n) \mathbf X(n)^T]$
上式中为了使零阶公式成立,$\mu_B$必须小鱼$\frac{1}{\lambda_{max}}$,其中$\lambda_{max}$是相关矩阵的最大特征值.
失调

$\upsilon = \frac{\mu_B}{2L}tr[\mathbf R] \tag {34}$

其中$tr[\mathbf R]$是相关矩阵的迹.

块长的选择

滤波器长度M和块长度L的关系有三种可能： 1.$L = M$，从计算的复杂性上看，最佳。
2.$L < M$，有降低延迟的好处。 3.$L > M$，将产生自适应过程冗余运算。

FLMS(频域LMS)

不论是webRTC还是speex开源的AEC算法都是基于频域来做的。之所以放在频域而非时域实现的主要原因实时性，在16kHz采样率的情况下，屋子里的回声持续时间长达0.1~0.5秒（多次反射），这就要求自适应滤波器的抽头数达到8000之多，工程上在考虑到计算量和延迟因素时基本都选择在频域实现。 FLMS（Fast LMS）的基本思想是将时域块LMS放到频域来计算。利用FFT算法在频域上完成滤波器系数的自适应。快速卷积算法用重叠相加法和重叠存储法。重叠相加法是将长序列分成大小相等的短片段，分别对各个端片段做FFT变换，再将变换重叠的部分相加构成最终FFT结果，重叠存储法在分段时，各个短的段之间存在重叠，对各个段进行FFT变换，最后将FFT变换得结果直接相加即得最终变换结果。当块的大小和权值个数相等时，运算效率达到最高。
根据重叠存储方法，将滤波器M个抽头权值用等个数的零来填补，并采用N点FFT进行计算，其中$N = 2M$，因此，N*1的向量：

$\hat {\mathbf W}(k) = FFT \begin{bmatrix} \hat {\mathbf w}(k)\ \mathbf 0 \end{bmatrix} \tag {35}$

表示FFT补零后的系数，抽头权向量为 $\hat {\mathbf w}(k)$。值得注意的是频域权向量$\hat {\mathbf W}(k)$的长度是时域权向量$\hat {\mathbf w}(k)$长度的两倍。相应的令：

$X(k) = diag\left \{ FFT \begin{bmatrix} x(kM-M),\cdots,x(kM-1)& x(kM),\cdots, x(kM+M-1)\\ K-1, block & K,block \end{bmatrix}\right \} \tag{36}$

表示对输入数据的两个相继字块进行傅里叶变换得到一个$N * N$对角阵.
将重叠存储法应用于3.26得:

$y^T(k) = \begin{bmatrix} y(kM) & y(kM+1) & \cdots & y(kM+M-1) \end{bmatrix} = IFFT[\mathbf X(k) \hat {\mathbf W}(k)], last \mathbf M \tag {37}$

每处理一帧,式37只有最后的M个元素被保留,因为前面的N个元素是循环卷积的结果.
设第$K$块的$M*1$期望响应和误差信号分别如下:

$\mathbf d(k) = \begin{bmatrix} d(kM) & d(kM+1) & \cdots & d(kM+M-1) \end{bmatrix}^T$
$\mathbf e(k) = \begin{bmatrix} e(kM) & e(kM+1) & \cdots & e(kM+M-1) \end{bmatrix}^T = \mathbf d(k) - \mathbf y(k)$

根据式3.37,可将$\mathbf e(k)$变换到频域,即

$E(k) = FFT\begin{bmatrix} \mathbf 0\\ \mathbf e(k) \end{bmatrix} \tag {38}$

则在更新权值的相关矩阵如下:

$\Phi(k) = \sum \limits_{i=0}^{L-1}x(kL+i)e(kL_i) = IFFT[\mathbf X^T(k) \mathbf E(k)] \tag{39}$,

的最前面M个元素,则抽头更新过程在频域中的表现如下:

$\hat {\mathbf W}(k+1) = \hat {\mathbf W}(k) + \mu FFT\begin{bmatrix} \Phi(k)\\ \mathbf 0 \end{bmatrix} \tag {40}$

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) \tag{41}$

  权值更新如下:

$\hat w_k(n+1) = \hat w_k(n) + \mu \frac{e(n)x^*(n-k)}{\sum_{i=0}^{N-1}|x(n-i)|^2} = \hat w_k(n) + \mu \frac{d(n) - \sum_i \hat w_i(n)x(n-i)x^*(n-k)}{\sum_{i=0}^{N-1}|x(n-i)|^2} \tag {42}$

其中$x(n)$是参考信号，$\hat w_k(n)$是$n$时刻和步长$\mu$的权值更新。 假设滤波后的误差为

$\delta_k(n) = \hat w_k(n) - w_k(n), d(n) = v(n) + \sum_k \hat w_k(n) x(n-k) \tag{43}$

则误差的迭代关系如下：

$\delta_k(n+1) = \delta_k + \mu \frac{v(n) - \sum_i \delta_i(n) x(n-i) x^*(n-k)}{\sum_{i=0}^{N-1}|x(n-1)|^2} \tag {44}$

在每一次调节中，滤波器的误差估计为$\Lambda = \sum_k \delta_k^*(n)\delta_k(n)$,展开后得如下形式：

$\Lambda(n+1) = \sum \limits_{k=0}^{N-1}\begin{vmatrix} \delta_k(n) + \mu \frac{v(n) - \sum_i \delta_i(n) x(n-i) x^*(n-k)}{\sum_{i=0}^{N-1}|x(n-1)|^2} \end{vmatrix} \tag {45}$

如果$x(n)$和$v(n)$是不相关的白噪声信号，则如下式:

$E\{\Lambda(n+1)| \Lambda(n), x(n)\} = \Lambda (n) \begin{bmatrix} 1 - \frac{2\mu}{N} + \frac{\mu^2}{N} + \frac{\mu^2\delta_y^2}{\Lambda(n)\sum_{i=0}^{N-1}|x(n-1)|^2} \end{bmatrix} \tag{46}$

可以通过求解

$\partial E\{ \Lambda(n+1)\}/ \partial \mu = 0, \Lambda \ne 0$:
$\frac{-2}{N}+ \frac{2\mu}{N}+ \frac{2\mu \sigma_y^2}{\Lambda(n) \sum_{i=0}^{N-1}|x(n-1|^2} = 0 \tag{47}$

求解后得到最优步长：

$\mu_{opt}(n) = \frac{1}{1+ \frac{\sigma_y^2}{\Lambda(n)(1/N)\sum_{i=0}^{N-1}|x(n-i)|^2}} \tag{48}$

期望$\Lambda(n)(1/N) \sum_{i=0}^{N-1}|x(n-i)|^2$等于剩余回声的方差$\sigma_r^2(n)$，如果剩余回声的方差值等于0，则步长因子等于1，$r(n) = y(n) - \hat y(n)$，则有输出信号的方差是：

$\sigma_e^2(n) = \sigma_v^2(n) + \sigma_r^2(n) \tag{49}$

这样可以求得这种情况下的最优步长因子为：

$\mu_{opt} \approx \frac{\sigma_r^2(n)}{\sigma_e^2(n)} \tag {50}$

则最优步长因子如下：

$\hat \mu_{opt}(n) = min(\frac{\hat \sigma_r^2}{\hat \sigma_e^2(n)}, 1) \tag {51}$

当$\Lambda(n) \approx \frac{\sigma_v^2}{\sigma_x^2(\frac{2}{\mu})-1}$时，式45停止迭代。将50带入51得到在滤波器系数不更新情况的剩余回声：

$\delta_r^2(n) \approx min(\frac{1}{2}\delta_r^2(n), \sigma_v^2(n)) \tag {52}$

• 频域解

和时域相比，频域可以使步长因子$\mu(k, l)$按频域划分，$Y(k,l)$和$E(k,l)$分别是频域中的记号，其和时域中的$\hat y(n)$和$e(n)$是对等的关系，k是频域索引，l是帧索引，和3.49类似，可得频域步长因子：

$\mu_{opt}(k,l) \approx \frac{\sigma_r^2(k,l)}{\sigma_e^2(k,l)} \tag {53}$

假设滤波器有一个和频谱无关的泄露（滤波器的误差）系数$\eta(l)$，这将得到：

$\hat \sigma_r^2(k,l) = \hat{\eta}(l) \hat{\sigma}_{\hat Y}^2(k, l) \tag {54}$

$\eta(l)$实际上是滤波器回声返回损失增强ERLE。
为了让步长因子调节的更快，使用瞬时估计，

$\hat \sigma_y(k,l) = |Y(k,l)|^2$和$\hat \sigma_e(k,l) = |E(k,l)|^2 \tag{55}$

这将使得3.50步长因子变为：

$\hat \mu_{opt}(k,l) = min(\hat \eta(l) \frac{\hat Y(k,l)|^2}{|E(k,l)|^2}, \mu_{max}) \tag {56}$

$\mu_{max}$是小于等于1的数，以确保滤波器稳定。

工程中的思考

关于步长因子

对于视频通话这类场景，两个通信终端的时钟偏斜和漂移是不定的，而音箱场景这个是可以在硬件上加以解决的，但是音箱场景的非线性失真却比通信场景严重的，功率放大模块非线性器件带来的谐波失真，在室内四个方向都发声，是得卷积失真，多次反射回声，声音突变等会加剧问题处理的复杂性；
当前绝大部分的AEC算法基本都基于频域分块处理方法，基于LMS/NLMS、RLS(recursive least square), APA(Affine Projection Algorithm)自适应处理方法。该自适应方的基本公式是：

$\mathbf h(n+1) = \mathbf h(n) +\mu e(n) \mathbf x(n) \tag{57}$

这个公式中误差信号e是维纳滤波（相减）后可以计算得到的，而步长因子却无法直接求得，在有些场景中，根据线性代数推导，在稳定性和收敛速度的双重约束下可以得到步长因子$\mu$要满足步长因子小于输入信号协方差矩阵的迹的导数：

$\mu \lt \frac{1}{trace(R_{xx})} \le \frac{1}{\lambda_{max}} \tag {58}$

针对LMS的一个改进是NLMS算法，这个方法根据输入输入信号的功率对步长进行归一化：

$\mu = \frac{\beta}{|x(n)|^2} \tag{59}$

其中$\beta$是归一化步长因子，$0\lt \beta \lt 2$.为了防止分母为零，加小数a，得如下计算公式：

$h(n+1)= h(n) + \frac{\beta}{|x(n)|^2+a} x(n)e(n) \tag{60}$

这样收敛速度和输入信号的功率绝对值无关。但是这两种算法在输入信号相关性在很高时收敛速度会比较慢。更进一步的可以采用自适应NLMS算法：

关于误差信号

LMS算法消除了线性部分，得到的是残余线性部分和非线性部分之和，非线性部分源于外界的噪声，参考源的卷积响应以及喇叭或者传输路径的影响，常用的方法可以

RERL:ERL+ERLE
RERL:residual_echo_return_loss
ERL:echo_return_loss
ERLE:echo_return_loss_enhancement
psd:power spectral density 功率谱密度
x: far end
d: near end
e: error
s: psd
nlp:non-linear processing

ERL vs ERLE

ERL :echo return loss = (mic in power) / (ref power)

ERLE :echo return loss enhancement = (power of residual echo signal) / (mic in power)

ERL和麦克风采集到的带回声信号和参考信号的比值，单位常用dB表示，比值越高，反映的是回声信号越小，大多数的VoIP设备的ERL值在15~20dB之间。

REC（residual echo control）

使用滤波器来近似非线性响应函数，这样也会得到非线性估计，实际上最早开始部分提到的误差信号e是减去了线性和非线性部分得到的误差信号，滤波器系数可以通过最小能量均分来做为准则进行平滑。大部分AEC算法的较大差异在非线性部分/后处理部分，通常希望滤波器阶数较高，能够处理较长的时间长度（即信号经过多次反射到达的场景）。比较有名的处理的滤波器是volterra 滤波器，使用一阶，二阶和三阶滤波器来去除非线性部分，但随着阶数的增加，计算量也呈现指数方式增加，其一阶的表示式如下：

$\mathbf X_1[k] = (x[k], x[k-1], \cdots, x[k-M+1]) \tag {61}$ $\mathbf {\hat h_1}= (\hat h_1[0], \hat h_1[1], \cdots, \hat h_1[M-1]) \tag{62}$

二阶的表示如下：

$\mathbf x_2[k] = (x^2[k], x[k]x[k-1], \cdots, x[k]x[k-M+1], x[k-1]x[k-1], \cdots, x[k-M+1]x[k-M+1]) \tag {63}$

$\mathbf {\hat h}_1 = (\hat h_2[0,0], \hat h_2[0,1], \cdots, \hat h_2[0, M-1], \hat h_2[1,1], \cdots, \hat h_2[M-1, M-1]) \tag{63}$

则误差信号可以表示为：

$e[k] = y[k] -\hat h_0[k] -\mathbf {\hat h}_1[k] \mathbf {x}_1^T[k] -\mathbf {h}_2[k] \mathbf {x}_2^T[k] \tag {64}$

则滤波器的系数更新方程如下：

$\hat h_0[k+1] = {\hat h}_0[k] + \mu_0e[k] \tag{65}$ $\mathbf {\hat h}_1[k+1] = \mathbf {\hat h}_1[k] + \mu_0e[k]\mathbf x_1[k] \tag{66}$ $\mathbf {\hat h}_2[k+1] = \mathbf {\hat h}_2[k] + \mu_0e[k]\mathbf x_2[k] \tag{67}$

这里涉及三个步长更新速率，$\mu_0$, $\mu_1 \lt \frac{2}{||\mathbf x_1||_2^2}$, $\mu_2 \lt \frac{\alpha}{||\mathbf x_1||_2^2 ||\mathbf x_2||_2^2}$,当$0 \lt \alpha \lt 2$时，这一过程是收敛的。

非线性问题

由于扬声器的失真以及传播的非理想性，维纳滤波并不能完美的解决AEC的问题，一般在经过维纳滤波之后会再经过RES和/或NLP,这个模块的核心思想是提取信号的有用信息(如，误差信号和远端信号的互相关性以及误差信号和近端信号的相关性，以及远端和近端信号的相关性)计算相应频点的缩放因子，以实现在输出的误差信号中进一步压制误差信号。

结束语

本章主要阐述了AEC要解决的问题，以及常用的时域和频域解决方法，值得一提的是，产品结构布局影响还是很大的，市场上绝大多数智能音响产品的喇叭和麦克风一般在物理尺寸的两端(圆柱形的外形,一般麦克在顶部,喇叭在中下部，且喇叭开口朝下，通过反射锥以弥补听感上损失)，可以收音和发音方向相反且隔开的。
有了以上理论，可以从开源的实例(WebRTC MATLAB and c code github Address)，下一章以webRTC为实例，剖析MATLAB和C代码，在我的github上有相关开源代码的使用实例。需要注意的是我代码里给定的默认设置值，不一定针对你的场景是最优的，但是效果一定是可以听出来的。