WebRtc AEC核心算法之一：频域自适应滤波

WebRtc AEC核心算法之一：频域自适应滤波

WebRtc和Speex作为目前开源的语音增强平台，给非科班出身的工程师一探究竟的机会，本文接下来以WebRtc的Aec模块为主线，研究一下比较核心的算法和实现。
一直以来在心目中AEC就是利用自适应滤波器求解梯度的迭代过程，时域的LMS/RLS经典算法在教科书和主流文章中都有介绍，但频域的介绍的就不多了，恰恰两个流行的平台都选用了频域自适应滤波器，首先是因为频域方法避免了时域FIR滤波器抽头太多，计算复杂度较大的缺点，另外就是最近二三十年频域滤波的发展迅速，很多优秀的算法已经产品化，只是我们接触的少不知道而已。在《Acoustic MIMO Signal Processing》第五章 Frequency-Domain Adaptive Filters，有以下介绍：相比于时域方法，频域分块处理的计算复杂度显著的降低了。FFT的运用，令输入信号的时移属性非常适合Toeplitz矩阵结构，意味着设计一个频域自适应滤波算法仅仅是通过被显式的Toeplitz和循环矩阵重写时域的误差尺度而已。（说的真轻松！）据这本书介绍，最早是1978年Dentino等在《Adaptive filtering in the frequency domain》提出了频域滤波器，紧接着Ferrara在1980年的《Fast implementations of LMS adaptive filters》提出了频域最小均方差方法使得频域地维纳解实现了快速收敛。Mansour和Gray在1982年提出了《Unconstrained frequency-domain adaptive filter》，以及后来MDF（Multi delay filter is selected by speex）的提出和应用，这些基本奠定了频域自适应滤波器的基础。《An adaptive filtering algorithm using an orthogonal projection to an affine subspace and its properties》还介绍了 affine project方法，也是热门的话题。

自适应滤波器-从时域到频域

经典的LMS自适应滤波器是采用时域实现的，标准的权重因子$w(n)$迭代表达如下：
$$\begin{gathered} y(n)=w^H(n)u(n) \\ e(n)=d(n)-y(n) \\ w(n+1)=w(n)+\mu u(n)e(n) \end{gathered}$$
$u(n)$是输入矢量，$d(n)$是期望的相应，$w(n)$是FIR滤波器抽头，$\mu$是收敛因子。普通LMS的收敛因子是常数，而NLMS是对$\mu$做了归一化处理，加速收敛速度的变形。但是实际实现过程遇到了两个挑战：

1. 应用于回声消除等延时比较敏感的场景下，FIR的抽头（阶数）数量可观，意味需要较大的存储空间来存放权重因子，更为烦恼的是高阶卷积的运算量巨大，算法复杂度、稳定性等等都会遇到极大的挑战
2. 收敛速度和稳定性一直是自适应滤波器需要追求的一种平衡，能否在频域发现更加优秀的算法取得收敛速度和稳定性的统一考量。
   因此，频域自适应滤波(FDAF)、正交自适应滤波和子带自适应滤波被提出来，有书中归纳为块自适应滤波，因为输入信号序列$u(n)$被分成$L$点的数据块来处理，自适应过程是以块为单位来更新，这区别原始的以采样点为更新基础的算法。
   在进行分块自适应讨论之前，我们先将上面的式子展开，便于理解，假设滤波器抽头为N：
   $$\begin{array}{rl}\left[\begin{array}{c}{w}_0(k+1)\\ w_1(k+1)\\ ....\\ w_{N-1}(k+1)\end{array}\right]& =\left[\begin{array}{c}{w}_0(k)\\ w_1(k)\\ ....\\ w_{N-1}(k)\end{array}\right]+\mu \left[\begin{array}{c}u(k)\\ u(k-1)\\ ....\\ u(k-N+1)\end{array}\right]e(n)\\ & =\left[\begin{array}{c}{w}_0(k)\\ w_1(k)\\ ....\\ w_{N-1}(k)\end{array}\right]+\mu \left[\begin{array}{c}u(k)\\ u(k-1)\\ ....\\ u(k-N+1)\end{array}\right][d(n)-y(n)]\end{array}$$
   从这个公式很容易得出，每计算出一个滤波器输出$y(n)$，都需要更新权重因子，如果处理一帧数据是L个samples，那么随着输出更新L次，权重因子也更新L次，计算量是客观的。

块自适应滤波-block adaptive filtering[5]

假设输入数据块为向量表达：
$$u(n)=[u(n),u(n-1),...,u(n-N+1){]}^T$$
此时的权重因子向量为：
$$\hat{w}(n)=[\widehat{w_0}(n),\widehat{w_1}(n),...,{\hat{w}}_{M-1}(n){]}^T$$
对于$n$进一步的用如下定义来分解：
$$\begin{array}{rl}n=kL+i,& i=0,1,...L-1,\\ & k=1,2,...,\end{array}$$
$L$表示块长度，这样对于输入数据第$k$块，可以定义为：
$$\begin{array}{rl}{\mathbf{A}}^T(k)& ={u(kL+i)}_{i=0}^{L-1}\\ & =[u(kL),u(kL+1),...,u(kL+L-1)]\end{array}$$
这个时候容易瞢，回头看看$u(n)$可是$M$维数组，所以$A(k)$是一个矩阵，假设$M=6，L=4$，矩阵的格式如下
再举一例，由此我们得到$y(n)$的更新过程如下图$(M=3，L=3)$：
这里很明显的看出$A(k)$是Toeplitz矩阵，主对角线以及平行主对角线的元素均相等。
block adaptive filtering权重因子更新算法是显著区别于经典算法的，下面只列出结果，有兴趣的同学可以参考[4][5]文献的介绍。
$$\begin{array}{rl}w(k+1)& =w(k)+\mu \sum _{i=0}^{L-1}\mathbf{u}(kL+i)e(kL+i)\\ & =w(k)+{\mathbf{A}}^T(k)e(k)\end{array}$$
这里的k是表示第k*L个模块，和之前时域的k不一样了。也就是每个L个samples算一次，我们可以把上式得矩阵表示进一步展开：
$$\begin{array}{rl}\left[\begin{array}{c}{w}_0(k+1)\\ w_1(k+1)\\ ....\\ w_{N-1}(k+1)\end{array}\right]& =\left[\begin{array}{c}{w}_0(k)\\ w_1(k)\\ ....\\ w_{N-1}(k)\end{array}\right]+\mu \left[\begin{array}{c}u(kL)\\ x(kL-1)\\ ....\\ x(kL-N+1)\end{array}\right]e(kL)+\mu \left[\begin{array}{c}u(kL+1)\\ x(kL)\\ ....\\ x(kL-N+2)\\ +\end{array}\right]e(kL+1)+\\ & ...+\mu \left[\begin{array}{c}u(kL+L-1)\\ x(kL+L-2)\\ ....\\ x(kL-N+L)\end{array}\right]e(kL+L-1)\end{array}$$
这样只需每个块L算一次，并且对这个块内的e(n)做了平滑处理。最后提一下，实际应用中往往采用Ｍ（N）＝Ｌ，这是综合考量计算复杂度等因素，对理解算法也有帮助，上文的第二个例子就说明了当M（N）=L之后，很容易设计出$A(k)$的Toeplitz矩阵形式，尤其下面引入频域和傅立叶变换之后。

Fast Block LMS Algorithm & FDAF

回顾上一节的滤波器输出，用公式来表示：
$$\begin{array}{rl}y(kL+i)& ={\hat{w}}^T(k)\mathbf{u}(kL+i)\\ & =\sum _{j=0}^{M-1}{\hat{w}}_j(k)u(kL+i-j),i=0,1,..,L-1.\end{array}$$
以及$\varphi (k)$的定义：
$$\begin{array}{rl}\varphi (k)& =\sum _{i=0}^{L-1}\mathbf{u}(kL+i)e(kL+i)\\ & ={\mathbf{A}}^T(k)e(k)\end{array}$$
这两个公式可以用线性卷积结构实现，而FFT经过特殊的处理就可以满足线性卷积的处理需求，即50%重叠存储法，如果Ｍ个抽头的滤波器，补充Ｍ个0，构建一个N=2M长度的FFT：
$$\hat{W}(k)=FFT\left[\begin{array}{c}\hat{w}(k)\\ 0\end{array}\right]$$
构建频域对角阵
U ( k ) = d i a g { F F T [ u ( k M − M ) , . . . , u ( k M − 1 ) , ⎵ ( k − 1 ) t h b l o c k u ( k M ) , . . . , u ( k M + M − 1 ) , ⎵ ( k ) t h b l o c k ] } \bold{U}(k)=diag\lbrace FFT[\underbrace{u(kM-M),...,u(kM-1),}_{(k-1)th\quad block}\underbrace{u(kM),...,u(kM+M-1),}_{(k)th\quad block}]\rbrace U(k)=diag{FFT[(k−1)thblock u(kM−M),...,u(kM−1),​​(k)thblock u(kM),...,u(kM+M−1),​​]}
最后可以得到：
$$\begin{array}{rl}{y}^T& =[y(kM),y(kM+1),...,y(kM+M-1)]\\ & =lastMelementsofIFFT[\mathbf{U}(k)\hat{\mathbf{W}}(k)]\end{array}$$
而前Ｍ个元素被扔掉了。
定义一维期望信号向量和误差向量为：
$$\begin{array}{rl}d(k)& =[d(kM),d(kM+1),...,d(kM+M-1){]}^T\\ e(k)& =[e(kM),e(kM+1),...,e(kM+M-1){]}^T\\ & =\mathbf{d}(k)-\mathbf{y}(k)\end{array}$$
与输出向量的定义过程类似，定义一个2M维向量，但前Ｍ个用0填充：
$$E(k)=FFT\left[\begin{array}{c}0\\ \mathbf{e}(k)\end{array}\right]$$
将重叠保存的方法应用于线性相关，求得$\varphi (k)$：
$$\varphi (k)=firstMelementsofIFFT[U^H(k)E(k)]$$
而权重因子的频域更新公式如下：
$$\hat{W}(k+1)=\hat{W}(k)+\mu FFT\left[\begin{array}{c}\varphi (k)\\ 0\end{array}\right]$$
频域实现的基本数学过程大致如此。相比较Ns算法中利用谱减规则对频域功率谱分析的侧重点，此处FFT更多的是为了算法效率上的考虑。

参考文献

[1]：《Acoustic MIMO Signal Processing》Yiteng Huang etc. Springer
[2]：《频域自适应滤波算法及应用》华南理工大学-本科毕业设计
[3]：《频域块 LMS 自适应滤波算法的研究》田超 西安理工大学 - 硕士毕业论文
[4]：《Frequency Domain and MutiRate Adaptive Filtering》 John J. shynk. IEEE Digital Signal Process Magzine.
[5]：《Adaptive Filter Theory》Fifth Edition Simon Haykin.