学习梳理一下用于AEC的LMS频域实现方法以及基础知识

起源

least-mean-square方法是时域的，参考ANC 与 adaptive filter推导，能够得出一个迭代更新公式：$${\mathbf{W}}_{k+1}={\mathbf{W}}_k+2\mu e(k)X(k)$$
回声消除最终关注的是
$$e(n)=d(n)-y(n)=d(n)-{\mathbf{W}}^T\mathbf{X}(n)=d(n)-{\mathbf{X}}^T(n)\mathbf{W}$$
而
$$y(n)={\mathbf{W}}^T\mathbf{X}(n)={\mathbf{X}}^T(n)\mathbf{W}$$是这个滤波器的输出，抛开数学推导，这个实现范式其实简洁明了。但由于时域跟踪，每个采样点更新一次的计算量是非常巨大的，所以在上个世纪七八十年代涌现了频域lms实现的方法，目的是为了降低计算复杂度，从此，旧时王谢堂前燕，飞入寻常百姓家。人们开始慢慢的享受了无回声通话的便宜，这种情况直到AI的出现和普及。

换个维度，到频域

似乎时频变换是现代信号处理绕不过去的一个梗，我能追溯的文章是这篇Adaptive filtering in the frequency domain，这篇文章本身比较精简，唯一引申了The complex LMS algorithm，很显然后者是从数学推导的方式给了一个复数域更新lms参数以及得出误差信号的结果，并没有谈到时频变换的问题，而前者想到了如果把时域信号分块的变换到频域，采用每块更新参数的方法，将会大大减少计算量，按照文中的计算方法，1024点FFT的块更新方法，能节省的计算量非常的可观。但此文仍专注于频域自适应的拓展，即没有考虑线性、圆周卷积这个DFT之后绕不过去的问题，也没有给出滤波器长度和分块的关系，如【论文笔记之 AFiFD】Adaptive Filtering in the Frequency Domain一文的总结，“论文为实现 LMS 自适应滤波器提供了一个新的思路”，应该属于开创性的著作。90年提出的并被speex实现的Multidelay Block Frequency Domain Adaptive Filter引用了这一时期的Block implementation of adaptive digital filters和Fast implementations of LMS adaptive filters两篇论文，这两篇又大概是20世纪80年代初期的论述。两年后，92年一篇On the implementation of a partitioned block frequency domain adaptive filter (PBFDAF) for long acoustic echo cancellation引用了它，而这篇文章为webrtc第一代AEC的PBFDAF实现提供的算法基础。在继续学习之前，有必要恶补一下相关的基础知识。

线性卷积&循环卷积

这可能是要先理清的一个细节，虽然大学学过dft，工作中不断用fft来做事情，但很多特性和数学关系在这里需要重点的关注起来，线性和循环卷积来自matlab网站的介绍，让线性卷积和循环卷积变得等效，而更加简明扼要的论述可以参考Lecture 23: Circular Convolution，当然经典的数字信号处理教科书都会讲述这部分内容，但有时不用到可能就忽略了。

频域FIR的循环卷积等价于时域补0后的线性卷积

线性和循环卷积是为了解决DFT的问题而提出来的，的FIR滤波器搬到频域实现属于基本的数字信号处理问题，这个问题首先要解决的是圆周卷积产生的混叠，最好的办法是找到与线性卷积等价的策略，这个策略就是补0。假设一个长为$L$的有限长序列$x(n)$，经过一个长度$M$冲激响应为$h(n)$的线性系统$(FIR)$，输出$y(n)$的序列，表示如下：
$$y(n)=\sum _{k=0}^{M-1}h(k)x(n-k),x(n)=0,n<0orn\ge L;h(n)=0,n<0orn\ge M$$那么$y(n)$的长度是$L+M-1$。时域卷积等价于频域相乘（这句话背了很多遍，推导嘛。。。）我们可得：
$$Y(\omega )=X(\omega )H(\omega )$$但凡事都要细想一下，三个序列的长度不同，用$DFT$变换得用最大值那个，即$N\ge L+M-1$，此时针对于$L$的有限长序列$x(n)$和长度$M$冲激响应为$h(n)$，后面显然需要补0，否则数据就不对了。数字信号处理――原理、算法与应用（第四版）,JohnG.Proakis，DimitrisG.Manolakis著7章都实例。其中第二个例子讲述了混叠是怎么发生对的，并给出了一个结论，前$M-1$个结果发生混叠污染，而正是利用这一特性，引申出了$\mathbf{长数据分块处理}$的两个重要方法。

重叠保留&重叠相加—长数据序列滤波的考量

这里需要强调的时候，对应于时频分析类的算法，此类范式是不能够加窗的。Lecture 16 Linear Filtering with the DFT里有简单的介绍John G. Proakis 数字信号处理的离散傅里叶变换和应用章节有更加详细的论述。这里按照John G. Proakis 数字信号处理讲述的来记录分解一下两种方法的差别。

重叠保留

假设一个长长序列$x(n)$，经过一个长度$M$冲激响应为$h(n)$的线性系统$(FIR)$，输出$y(n)$的序列，为了进行分块分析，把序列分为长度为$L$的数据块，一般的$L>>M$，那么$N=L+M-1$作为FFT长度，对应第m个数据块，冲激响应补$L-1$个0来凑齐，但数据块，则在前面补充上一个数据块的最后$M-1$个点，完成拼接。将两个N点的DFT变换$H(k)$和$X_m(k)$相乘，得出：$${\hat{Y}}_m(k)=H(k)X_m(k);k=0,1,.....,N-1$$于是，我们通过逆变换，得出输出序列：
y ^ m ( n ) = { y ^ m ( 0 )   y ^ m ( 1 )   y ^ m ( 2 )   y ^ m ( 3 )   . . . y ^ m ( M − 1 )   y ^ m ( M )   . . . y ^ m ( N − 1 )   } \hat{y}_m(n) = \{\hat{y}_m(0)\ \hat{y}_m(1)\ \hat{y}_m(2)\ \hat{y}_m(3)\ ... \hat{y}_m(M-1)\ \hat{y}_m(M)\ ... \hat{y}_m(N-1)\ \} y^​m​(n)={y^​m​(0) y^​m​(1) y^​m​(2) y^​m​(3) ...y^​m​(M−1) y^​m​(M) ...y^​m​(N−1) }按照拼接原则，实际有效数据的前$M-1$个数据被混叠污染了，那么剩下的$L$个数据则被完美还原。即$$y_m(n-M)={\hat{y}}_m(n),n=M,M+1,...,N-1$$而为了满足这种拼接方法，每一块处理前的后M-1个点被保留下来，留给下一次用，第一块前面直接补0，得到这样一个送给DFT的数据：
x 1 ( n ) = { 0 , 0 , ⋯   , 0 , ⏟ M − 1 个点 x ( 0 ) , x ( 1 ) , . . . x ( L − 1 ) } ⏟ 第一帧 L 个点 x 2 ( n ) = { x ( L − M + 1 ) , , ⋯   , x ( L − 1 ) , ⏟ 保留 x 1 最后 M − 1 个点 x ( L ) , x ( L + 1 ) , . . . x ( 2 L − 1 ) } ⏟ 第二帧 L 个点 x 3 ( n ) = { x ( 2 L − M + 1 ) , , ⋯   , x ( 2 L − 1 ) , ⏟ 保留 x 2 最后 M − 1 个点 x ( 2 L ) , x ( 2 L + 1 ) , . . . x ( 3 L − 1 ) } ⏟ 第三帧 L 个点 \begin{align} x_1(n) = \begin{matrix} \underbrace{ \{0,0,\cdots,0,} \\ {M-1个点} \end{matrix}\begin{matrix} \underbrace{ x(0),x(1),...x(L-1)\}} \\ {第一帧L个点} \end{matrix} \newline x_2(n) = \begin{matrix} \underbrace{ \{x(L-M+1),,\cdots,x(L-1),} \\ {保留x_1最后M-1个点} \end{matrix}\begin{matrix} \underbrace{ x(L),x(L+1),...x(2L-1)\}} \\ {第二帧L个点} \end{matrix}\newline x_3(n) = \begin{matrix} \underbrace{ \{x(2L-M+1),,\cdots,x(2L-1),} \\ {保留x_2最后M-1个点} \end{matrix}\begin{matrix} \underbrace{ x(2L),x(2L+1),...x(3L-1)\}} \\ {第三帧L个点} \end{matrix} \end{align} x1​(n)= {0,0,⋯,0,​M−1个点​ x(0),x(1),...x(L−1)}​第一帧L个点​x2​(n)= {x(L−M+1),,⋯,x(L−1),​保留x1​最后M−1个点​ x(L),x(L+1),...x(2L−1)}​第二帧L个点​x3​(n)= {x(2L−M+1),,⋯,x(2L−1),​保留x2​最后M−1个点​ x(2L),x(2L+1),...x(3L−1)}​第三帧L个点​​​
实用的情况是一般取每帧N点数据通过N阶滤波器，而FFT长度是2N，上式则变为：
x 1 ( n ) = { 0 , 0 , ⋯   , 0 , ⏟ N 个点 x ( 0 ) , x ( 1 ) , . . . x ( N − 1 ) } ⏟ 第一帧 N 个点 x 2 ( n ) = { x ( 0 ) , , ⋯   , x ( N − 1 ) , ⏟ 保留 x 1 最后 N 个点 x ( N ) , x ( N + 1 ) , . . . x ( 2 N − 1 ) } ⏟ 第二帧 N 个点 x 3 ( n ) = { x ( N ) , , ⋯   , x ( 2 N − 1 ) , ⏟ 保留 x 2 最后 N 个点 x ( 2 N ) , x ( 2 N + 1 ) , . . . x ( 3 N − 1 ) } ⏟ 第三帧 N 个点 \begin{aligned} x_1(n) = \begin{matrix} \underbrace{ \{0,0,\cdots,0,} \\ {N个点} \end{matrix}\begin{matrix} \underbrace{ x(0),x(1),...x(N-1)\}} \\ {第一帧N个点} \end{matrix} \newline x_2(n) = \begin{matrix} \underbrace{ \{x(0),,\cdots,x(N-1),} \\ {保留x_1最后N个点} \end{matrix}\begin{matrix} \underbrace{ x(N),x(N+1),...x(2N-1)\}} \\ {第二帧N个点} \end{matrix}\newline x_3(n) = \begin{matrix} \underbrace{ \{x(N),,\cdots,x(2N-1),} \\ {保留x_2最后N个点} \end{matrix}\begin{matrix} \underbrace{ x(2N),x(2N+1),...x(3N-1)\}} \\ {第三帧N个点} \end{matrix} \end{aligned} x1​(n)= {0,0,⋯,0,​N个点​ x(0),x(1),...x(N−1)}​第一帧N个点​x2​(n)= {x(0),,⋯,x(N−1),​保留x1​最后N个点​ x(N),x(N+1),...x(2N−1)}​第二帧N个点​x3​(n)= {x(N),,⋯,x(2N−1),​保留x2​最后N个点​ x(2N),x(2N+1),...x(3N−1)}​第三帧N个点​​输出则只保留最后的N个点。

重叠相加

所有假设与重叠保留一致，只是拼帧的时候，数据块补充M-1个0,并计算DFT,这样得到的结果是没有任何混叠的L+M-1个点。我们通过逆变换，得出输出序列：
y ^ m ( n ) = { y ^ m ( 0 )   y ^ m ( 1 )   y ^ m ( 2 )   y ^ m ( 3 )   . . . y ^ m ( L − 1 )   y ^ m ( L )   . . . y ^ m ( N − 1 )   } \hat{y}_m(n) = \{\hat{y}_m(0)\ \hat{y}_m(1)\ \hat{y}_m(2)\ \hat{y}_m(3)\ ... \hat{y}_m(L-1)\ \hat{y}_m(L)\ ... \hat{y}_m(N-1)\ \} y^​m​(n)={y^​m​(0) y^​m​(1) y^​m​(2) y^​m​(3) ...y^​m​(L−1) y^​m​(L) ...y^​m​(N−1) }
虽然没有混叠，但数据输出长度多出了M-1个，这M-1个信息需要与后面的数据对应相加，完成卷积的交叠拼接。最终得到： y ( n ) = { y 1 ( 0 ) ,   y 1 ( 1 ) ,   . . .   y 1 ( L − 1 ) ,   y 1 ( L ) + y 2 ( 0 ) ,   y 1 ( L + 1 ) + y 2 ( 1 ) ,   . . . y 1 ( N − 1 ) + y 2 ( M − 1 ) , y 2 ( M )   . . . } {y}(n) = \{{y}_1(0),\ {y}_1(1),\ ...\ {y}_1(L-1),\ {y}_1(L)+y_2(0), \ {y}_1(L+1)+y_2(1), \ ...{y}_1(N-1)+y_2(M-1), y_2(M)\ ...\} y(n)={y1​(0), y1​(1), ... y1​(L−1), y1​(L)+y2​(0), y1​(L+1)+y2​(1), ...y1​(N−1)+y2​(M−1),y2​(M) ...}这种情况的拼接关系比较简单： x 1 ( n ) = { x ( 0 ) , x ( 1 ) , . . . x ( L − 1 ) ⏟ 第一帧 L 个点 0 , 0 , ⋯   , 0 , } ⏟ M − 1 个点 x 2 ( n ) = { x ( L ) , x ( L + 1 ) , . . . x ( 2 L − 1 ) ⏟ 第二帧 L 个点 0 , 0 , ⋯   , 0 , } ⏟ M − 1 个点 x 3 ( n ) = { x ( 2 L ) , x ( 2 L + 1 ) , . . . x ( 3 L − 1 ) ⏟ 第三帧 L 个点 0 , 0 , ⋯   , 0 , } ⏟ M − 1 个点 \begin{align} x_1(n) = \begin{matrix} \underbrace{ \{x(0),x(1),...x(L-1)} \\ {第一帧L个点} \end{matrix}\begin{matrix} \underbrace{ 0,0,\cdots,0,\}} \\ {M-1个点} \end{matrix} \newline x_2(n) = \begin{matrix} \underbrace{\{ x(L),x(L+1),...x(2L-1)} \\ {第二帧L个点} \end{matrix}\begin{matrix} \underbrace{ 0,0,\cdots,0,\}} \\ {M-1个点} \end{matrix}\newline x_3(n) = \begin{matrix} \underbrace{\{ x(2L),x(2L+1),...x(3L-1)} \\ {第三帧L个点} \end{matrix}\begin{matrix} \underbrace{ 0,0,\cdots,0,\}} \\ {M-1个点} \end{matrix} \end{align} x1​(n)= {x(0),x(1),...x(L−1)​第一帧L个点​ 0,0,⋯,0,}​M−1个点​x2​(n)= {x(L),x(L+1),...x(2L−1)​第二帧L个点​ 0,0,⋯,0,}​M−1个点​x3​(n)= {x(2L),x(2L+1),...x(3L−1)​第三帧L个点​ 0,0,⋯,0,}​M−1个点​​​

block分块自适应探讨

频域处理和分块一体不可区分的，分块处理大大减轻了自适应的算法负担，块内只算一次参数更新（输出一块的数据也只在块边界更新一次（相对于sample）），Block implementation of adaptive digital filters被引用的是比较多的一篇，A Unified Approach to Time- and Frequency-Domain Realization of FIR Adaptive Digital Filters也阐述了block lms的一些特性，但这篇文章有点长，而且时间久远，晦涩难懂。由于这个领域的诱人前景，美国OSTI在80年左右召开了一个Block adaptive filtering的会议。讲到这里，我们回顾一下，上一章节的重叠保留&重叠相加方法为块处理数字滤波器奠定了基础，但要实现块处理自适应滤波器，不仅仅要完成时-频-时变换，还有一个自适应参数更新的过程，这个如何更新？在时域还是频域做？不同策略会带来什么样的计算效率？在想实战之前，应该先撸一撸基础的理论和方法，但论文眼花缭乱，从何撸起呢？综合几篇大作，并拜读浩哥依然选三篇，或者说是2+1，1是Block implementation of adaptive digital filters，而2是这篇引用的Adaptive filtering in the frequency domain和Fast implementations of LMS adaptive filters，不妨对比着看看后两篇讲述了啥，尝试解决什么问题。从时间线上来说Adaptive filtering in the frequency domain发表于1978年，文章简短的描述了从利用The complex LMS algorithm算法从频域上计算的可行性和计算效率，但没有追究时频算法的等价性和相关性，以及误差收敛等关键问题，文中的权重值应该都是在频域存在的，计算量非常经济实惠，看着也似乎比较容易实现。1980年发表的Fast implementations of LMS adaptive filters引用了前者，这也说明了一种传承关系，后者推导的时候严格遵守了“重叠保留”方法，有了前文的铺垫，这里不难理解，文中的创新点应该是参数更新部分的时频等价变换这里了。Block implementation of adaptive digital filters侧重点在于论述了分块自适应的数学推导，给出了理论上分块方法的基础，并没有特意强调在频域还是时域。

权重向量的时频约束

权重$w(n)$在FDAF的更新是值得推敲的，在推导之前先记住这个英文constraint，中文的意思“限制;约束;限定”，constraint是为了用圆周卷积完美实现线性卷积的时频运算而出现的，参考“重叠保留”法中对数据块和冲击响应的后半部分清零的“constraint”，IFFT之后只保留后半部分数据块儿的“constraint”，权重$w(n)$在时频变换过程中处理的更加烧脑一些，1992年初发表的这篇Frequency-domain and multirate adaptive filtering给出了很细致的推导，结合上文中Fast implementations of LMS adaptive filters两篇文章，相辅相成的理解可能更好，推导之前需要回顾一下卷积、相关的关系，卷积的典型公式为$$y(n)=\sum _{k=0}^{M-1}h(k)x(n-k),x(n)=0,n<0orn\ge L;h(n)=0,n<0orn\ge M$$而相关的公式为：$$y(n)=\sum _{k=0}^{M-1}h(k)x(k),x(n)=0,n<0orn\ge L;h(n)=0,n<0orn\ge M$$在不考量下标长度的情况下，两者就差一个反折（卷积）过程，脑补一下经过FFT-IFFT的过程，混叠发生的位置应该正好是相反的（推导没找到，也懒得自己敲），那么一切参考着卷积的重叠保留，相关的重叠保留都是镜像的。有了这剂预防针，下面的推导可能稍微好理解一下，敲一下Fast implementations of LMS adaptive filters的公式，首先定义抽头延迟线某一时刻$j$的滤波器权重向量$W_j^T$和输入向量$X_j^T$为：$$\begin{gathered} W_j^T=[w_{1,j},w_{2,j},...,w_{n,j}](2) \\ {X}_j^T=[x_j,x_{j-1},...,x_{j-n+1}](3) \end{gathered}$$滤波器输出$y_j$的表示为$$y_j=X_j^T{W}_j(5)$$定义误差信号$error$为${\epsilon }_j$，表示时刻$j$的“期望响应$d_j$”与滤波器输出$y_j$的差：$${\epsilon }_j=d_j-y_j(4)$$权重向量的更新公式为：$$W_{j+1}=W_j+2\mu {\epsilon }_j{X}_j(1)$$进而推广到块方法，在第$k$th个block，滤波器输出为：$$y_{kn+i}=X_{kn+i}^T{W}_{k}0⩽i⩽n-1(6)$$利用公式$(1)$迭代计算块权重更新公式为：$$\begin{gathered} W_{k+1}=W_k+2\mu \sum _{i=0}^{n-1}{\epsilon }_{kn+i}{X}_{kn+i} \\ {W}_{k+1}=W_k+2\mu {\nabla }_k(7) \end{gathered}$$以上为款更新的时域推导，再次基础上做一个时频等效实现，尤其在权重值计算上能够节省很大计算工作量的话，那么$Fast$目的就实现了，为此重新展开公式(6)：$$y_{kn+j}=\sum _{i-0}^{n-1}{w}_{i,k}{x}_{kn+j-i}0⩽j⩽n-1(8)$$要利用重叠保留实现FFT快速运算，首先选取$2n$长度的FFT，对权重值的后半部分补“0”，我们得到频域权重${\omega }_k$：$${\omega }_k^T=FFT[W_k^T\underbrace{\mathrm{0....0}}](9)$$对应我们需要准备2n个点的输入数据： X k T = F F T [ { x ( k − 1 ) n , x ( ( k − 1 ) n + 1 ) , . . . x k n − 1 ⏟ 第 K − 1 帧 n 个点 x k n , ⋯   , x k n + n − 1 } ⏟ 第 K 帧 n 个点 ] ( 10 ) \mathcal{X}_k^T=FFT[\begin{matrix} \underbrace{ \{x_{(k-1)n},x_{((k-1)n+1)},...x_{kn-1}} \\ {第K-1帧n个点} \end{matrix}\begin{matrix} \underbrace{ x_{kn},\cdots,x_{kn+n-1}\}} \\ {第K帧n个点} \end{matrix} ]\quad\quad(10) XkT​=FFT[ {x(k−1)n​,x((k−1)n+1)​,...xkn−1​​第K−1帧n个点​ xkn​,⋯,xkn+n−1​}​第K帧n个点​](10)那么我们可以得到公式（8）的另外一种表达： [ y k n . . . y k n + n − 1 ] T = l a s t n t e r m s o f F F T − 1 { ω k ⊗ X } [y_{kn}...y_{kn+n-1}]^T=last\quad n\quad terms\quad of\quad FFT^{-1}\{\omega_k\otimes\mathcal{X}\} [ykn​...ykn+n−1​]T=lastntermsofFFT−1{ωk​⊗X}要实现权重向量的时频更新，即将公式（7）在频域实现，这时候我们要注意到我们已经将滤波器信号返回到时域，得到了$[y_{kn}...y_{kn+n-1}{]}^T$，那么此时我们能得到一个误差信号$${\epsilon }_{kn+j}=d_{kn+j}-y_{kn+j}0⩽j⩽n-1$$公式（7）的关键是梯度$\nabla$的转换，重写第k帧，第j个向量的梯度为：$${\nabla }_{j,k}=\sum _{i=0}^{n-1}{\epsilon }_{kn+i}{x}_{kn+i-(j-1)}0⩽j⩽n-1$$一个点的梯度需要这么多乘累加，不如再次回到频域吧，这次我们要的相关的结果，所以反向填充0得到误差信号的FFT： E k = F F T [ { 0 , 0 , . . . 0 ⏟ n 个点 d k n − y k n , ⋯   , d k n + n − 1 − y k n + n − 1 } ⏟ 第 K 帧 n 个点误差 ] T ( 13 ) E_k=FFT[\begin{matrix} \underbrace{ \{0,0,...0} \\ {n个点} \end{matrix}\begin{matrix} \underbrace{ d_{kn}-y_{kn},\cdots,d_{kn+n-1}-y_{kn+n-1}\}} \\ {第K帧n个点误差} \end{matrix} ]^T\quad (13) Ek​=FFT[ {0,0,...0​n个点​ dkn​−ykn​,⋯,dkn+n−1​−ykn+n−1​}​第K帧n个点误差​]T(13)时域相关（卷积）等于频域相乘，与$\mathcal{X}$相乘之后，在变换到时域，一样镜像的保留前半部分的数据，得到第k个block的时域梯度 ∇ k = f i r s t n t e r m s   o f F F T − 1 { E k ⊗ X ‾ } ( 14 ) \nabla_{k}= first\quad n\quad terms\quad\ of\quad FFT^{-1}\{E_k\otimes\mathcal{\overline{X}}\}\quad(14) ∇k​=firstnterms ofFFT−1{Ek​⊗X}(14)这是后用复共轭（迷惑），此时又回到了时域的梯度，然后我们时域已经没有了权重值，所以还要奔赴到频域去，为了配合权重值后半部分填0的操作，这个梯度的后半部分也填0，就得到如下的频域权重更新公式 ω i + 1 = ω k + 2 μ F F T [ ∇ k T { 0 , 0 , . . . 0 } ⏟ n 个点 ] T ( 15 ) \omega_{i+1}=\omega_{k}+2\mu FFT[\nabla_{k}^T\begin{matrix} \underbrace{ \{0,0,...0\}} \\ {n个点} \end{matrix}]^T\quad(15) ωi+1​=ωk​+2μFFT[∇kT​ {0,0,...0}​n个点​]T(15)

数据都分块了，冲激权重为啥不分个块

这个问题的提出是因为语音回声的冲激响应一般远大于FFT的长度，如果冲激响应也能分块，块块分割，使得实用性得到很大的提高，于是大约1991年底一篇文章横空出世，这就是On the implementation of a partitioned block frequency domain adaptive filter (PBFDAF) for long acoustic echo cancellation，该论文直接被webrtc拿来引用，n年前学过，回头看看大差不差WebRtc AEC核心算法之二：Partitioned block频域自适应滤波，另外发现了一个博主引用了外网的几张图来说明这个分块策略，暂时没看懂是，收藏慢慢理解吧webrtc AEC 线性滤波 PBFDAF（均匀分块频域自适应滤波）介绍，而他引用的是一篇外文Partitioned convolution algorithms for real-time auralization，这篇福就长了，同样收藏慢慢看吧。speex实现的Multidelay Block Frequency Domain Adaptive Filter还要早两年，两篇文章要搞清楚权重如何分块，以及更新策略。

横看成岭(PBFDAF)侧成峰(MDF)，对自适应参数分块的时频推导

借用之前的图，太长的时域激励（滤波器权重）可以分成相等的很多块去考虑。

将$h(n)$等分L个相邻的、N等长的而且不交叠的段：
$$h(n)=\sum _{l=0}^{L-1}{h}_l(n)$$
对于每个段，都有$h_l(n)=h(n)当n=lN,...,lN+N-1$。结合卷积和分块处理的理论，我们可以没N个点分成L个处理单元，各自卷积之后再相加。

而下N个点的处理，相当于整体X序列延时N点，或者经过一个FIFO，队列新进入N点，其他过程不变。滤波器的输出$y(n)=x(n)\ast h(n)$可以写成如下形式
$$\begin{array}{rl}y(n)& =x(n)\ast \sum _{l=0}^{L-1}{h}_l(n)\\ & =\sum _{l=0}^{L-1}x(n)\ast h_l(n)\\ & =\sum _{l=0}^{L-1}x(n-lN)\ast h_l(n+lN)\\ & =\sum _{l=0}^{L-1}{y}_l(n)\end{array}$$
这样patitioned block概念就出来了，如果要结合FDAF，就是PBFDAF过程，每新进入一块x(n)，做一次FFT，之前的块只需要延时一个FIFO节拍，最后各自block，多个不同延迟卷积的结果相加，就得到了输出y(n)，这种理解就叫做Multi_delay，MDF和PBFDAF在分块和时频变换这里实际是一回事儿。后面跟踪MDF的思路来推导，我们把L换成M，总长度为$M\ast N$保持不变，FFT的长度$N^{\prime }$要选择大于等于2MN/M的长度，索性我们做一些假设：16k采样率，N=256，M=8，$N^{\prime }$=512，滤波器抽头为M*N=2048，这些点除以采样率=0.128s，也就是这个配置能hold出128ms内的回声信号。采用重叠保留法，当前块的数据变换到频域： X ( M , j ) = d i a g { F F T [ x 0 ( j − 1 ) , x 1 ( j − 1 ) , . . . , x N − 1 ( j − 1 ) , x 0 ( j ) , x 1 ( j ) , . . . , x N − 1 ( j ) ] T } X(M,j)=diag\{FFT[x_0(j-1),x_1(j-1),...,x_{N-1}(j-1),\newline x_0(j),x_1(j),...,x_{N-1}(j)]^T\} X(M,j)=diag{FFT[x0​(j−1),x1​(j−1),...,xN−1​(j−1),x0​(j),x1​(j),...,xN−1​(j)]T}此时的下标$j$表示是一个时刻收到的N个点的数据，和之前$j-1$时刻收到N个点的数据的拼接。更早的延迟块输入向量可以通过FIFO的顺移，从历史得到：$$X(m,j)=X(m+1,j-1)m=1,2,...,M-1$$这里只有M-1个块，和新进入的第M块构成了完整时域序列的频域变换延迟块，最新的块下标是M。权重变换遵守重叠保留法：$$W(m,j)=FFT[h(m,j)\underbrace{\mathrm{0....0}}](9)$$那么最终得到卷积的时域表达为： y ( j ) = l a s t N t e r m s o f F F T − 1 { ∑ m = 1 M X ( m , j ) W ( m , j ) } y(j)=last\quad N\quad terms\quad of\quad FFT^{-1}\{\sum_{m=1}^MX(m,j)W(m,j)\} y(j)=lastNtermsofFFT−1{m=1∑M​X(m,j)W(m,j)}经典的误差信号约束 E ( j ) = F F T [ { 0 , 0 , 0 , . . . 0 , ⏟ 第 M − 1 帧 N 个点 [ d ( j ) − y ( j ) ] T } ⏟ 第 M 帧 N 个点 ] T E(j)=FFT[\begin{matrix} \underbrace{ \{0,0,0,...0,} \\ {第M-1帧N个点} \end{matrix}\begin{matrix} \underbrace{ [d(j)-y(j)]^T\}} \\ {第M帧N个点} \end{matrix}]^T E(j)=FFT[ {0,0,0,...0,​第M−1帧N个点​ [d(j)−y(j)]T}​第M帧N个点​]T权重更新利用这个误差信号和每个块的频域做相乘，再求逆运算得到M（8）组时域的误差：$$\varphi (m,j)=firsthalfofFF{T}^{-1}[X^{\ast }(m,j)E(j)]$$最后再频域： Φ ( m , j ) = F F T [ ϕ ( m , j ) , { 0 , 0 , 0 , . . . 0 } ⏟ N 个点 ] \Phi(m,j) = FFT[\phi(m,j), \begin{matrix} \underbrace{ \{0,0,0,...0\}} \\ {N个点} \end{matrix}] Φ(m,j)=FFT[ϕ(m,j), {0,0,0,...0}​N个点​]获得每个权重组的更新：$$W(m,j+1)=W(m,j)+M{\mu }_B\Phi (m,j)m=1,2,...,M$$不去探究$\mu$值的计算上述公式已经差不读讲明白这种推导了。 PBFDAF如果也按照这个走，那就没有新意了，所以到权重更新这块，他必须有所取舍。所以此文大胆采用了另外一种自适应算法：$$H_i^l(j+1)=H_i^l(j)+Inh(l){\mu }_i^l(j)E_i(j)X_i^{\ast }(j-l+1)fori=1,...,N+1,l=0,...,L-1$$上式中$Inh(l)$是一个键控值，根据一些策略可以决定是否要更新这个块的参数，而这里的E应该和MDF的一样的计算法，可以看出省略了后面两次FFT的过程，文中为此给出了详尽的解释（没太看懂），这样就避免了雷同，这个详细的推导来自于Unconstrained frequency-domain adaptive filter。当然后面还有DTDS等实用内容，就不展开分析了。

非约束的方法能用吗？要注意点什么

为了节省运算量，显然能省一个FFT就省掉一个FFT，这种不按时域-频域等价的自适应方法我们都理解成非约束的办法，在Multidelay block frequency domain adaptive filter的论文里，画了很多图，其中第一张摘录一下：

这张图讲述了误差信号是来自时域还是来自频域的一种拓扑比较，很显然之前的分析都是(a)中的形式，即误差从时域计算，而e(n)其实是回声消除后的信号，也是语音处理的输出。那么从这张图上我们可以看出(b)似乎更加有吸引力，他就频域完成了所有的关键运算，这是不是也能满足目标需求呢？

参考

The Comparison of the Constrained and Unconstrained Frequency-Domain Block-LMS Adaptive Algorithms
线性和循环卷积
Lecture 23: Circular Convolution