自适应滤波器分块分段技术详解
前言
先说下目的：对在频域计算FIR自适应滤波器，同时避免使用长滤波器时产生的大延迟的技术进行详细的分解。即要对滤波器分块又分段的过程进行细致的分析。这里假设读者有相应的LMS自适滤波器的基础
一、滤波器分析准备
先从时域LMS滤波器说起，设列向量
\[{\bf{x}}(n) = {\left[ {x(n),x(n - 1),...,x(n - M + 1)} \right]^T}\]
\[{\bf{w}}(n) = {\left[ {{w_0}(n),{w_1}(n),...,{w_{M - 1}}(n)} \right]^T}\]
这里列向量长度为M。考虑通过系数为w(n)的FIR滤波器对序列x(n)滤波，用卷积运算表示为\[y(n) = \sum\limits_{i = 0}^{M - 1} {{w_i}x(n - i)} \]
写成向量形式：
\[y(n) = {{\bf{w}}^T}(n){\bf{x}}(n)\]
简单从公式上，其实并不容易理解FIR滤波器的工作过程，这里换一种表达方式：把向量和矩阵用时间序列索引与符号来表示。应该就会比较好理解一些，对于序列：7，2，-3，-6，12，8，-7，-5，4，6。通过一个长度为4的FIR滤波器，输入向量 随时间序列的迭代变化表示为
\[\begin{array}{*{20}{c}}
{n = 0} \hfill & {{\bf{x}} = {{\left[ {\begin{array}{*{20}{c}}
7 & 0 & 0 & 0 \\
\end{array}} \right]}^T}} \hfill \\
{n = 1} \hfill & {{\bf{x}} = {{\left[ {\begin{array}{*{20}{c}}
2 & 7 & 0 & 0 \\
\end{array}} \right]}^T}} \hfill \\
{n = 2} \hfill & {{\bf{x}} = {{\left[ {\begin{array}{*{20}{c}}
{ - 3} & 2 & 7 & 0 \\
\end{array}} \right]}^T}} \hfill \\
{n = 3} \hfill & {{\bf{x}} = {{\left[ {\begin{array}{*{20}{c}}
{ - 6} & { - 3} & 2 & 7 \\
\end{array}} \right]}^T}} \hfill \\
{n = 4} \hfill & {{\bf{x}} = {{\left[ {\begin{array}{*{20}{c}}
{12} & { - 6} & { - 3} & 2 \\
\end{array}} \right]}^T}} \hfill \\
\end{array}\]
下面准备写成时间序列索引的形式(向量各元素用时间序列对应的索引号代替，如：[0 -1 -2 -3]T)，这样有利于在接下来的分析中看的更清楚一些，对于矩阵，也会沿用这样的方式。
二、滤波器分块处理
先从频域矩阵分块处理说起，分块方法是一种批处理方法，为了说的更清楚一些，这里用一个示例来说明，设滤波器长度M = 6，块长度N = 4，表示每次批处理4个单元
\[\left[ {\begin{array}{*{20}{c}}
\hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill { - 5} \\
\hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} \\
\hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} \\
\hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} \\
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{w_0}} \hfill \\
{{w_1}} \hfill \\
{{w_2}} \hfill \\
{{w_3}} \hfill \\
{{w_4}} \hfill \\
{{w_5}} \hfill \\
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
{{y_0}} \hfill \\
{{y_1}} \hfill \\
{{y_2}} \hfill \\
{{y_3}} \hfill \\
\end{array}} \right]\]
观察输入向量组成的矩阵可以发现，这是一个Toeplitz矩阵。要转到频域进行处理，一个比较可行的方式是把Toeplitz矩阵转为循环矩阵。再利用DFT把循环矩阵对角化。
那这个循环矩阵的大小是多少比较合适呢，由卷积过程可知，卷积后输出长度为L = M + N – 1 = 9，也就是说，用M + N – 1个卷积长度完全可以表达出来，也就是说循环矩阵C的大小为LxL应该是足够的。用循环矩阵表示如下：
\[\left[ {\begin{array}{*{20}{c}}
\hfill { - 5} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} \\
\hfill { - 4} & \hfill { - 5} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} \\
\hfill { - 3} & \hfill { - 4} & \hfill { - 5} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} \\
\hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill { - 5} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} \\
\hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill { - 5} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 \\
\hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill { - 5} & \hfill 3 & \hfill 2 & \hfill 1 \\
\hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill { - 5} & \hfill 3 & \hfill 2 \\
\hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill { - 5} & \hfill 3 \\
\hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill { - 5} \\
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{w_0}} \\
{{w_1}} \\
{{w_2}} \\
{{w_3}} \\
{{w_4}} \\
{{w_5}} \\
0 \\
0 \\
0 \\
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
\times \\
\times \\
\times \\
\times \\
\times \\
{{y_0}} \\
{{y_1}} \\
{{y_2}} \\
{{y_3}} \\
\end{array}} \right]\]
这里输出向量中的X部分是循环卷积的结果，是要舍弃的部分。重写以上过程
\[{\bf{y}}(n) = P_{L \times L}^{01}C{\bf{\hat w}}(n) = P_{L \times L}^{01}{F^{ - 1}}FC{F^{ - 1}}F{\bf{\hat w}}(n)\]
分解为4步：
- 是一个对角矩阵，在编程实现中，可以取循环矩阵C的第一列变换到频域来代替
- 是把后面补0后的滤波器系数变换到频域
- 是一个逆FFT变换
- 是为了只选取最后N个做为滤波器的输出结果，编程中做到这点很方便
再看滤波器的更新过程，这里把步长参数省去
\[w(n + 1) = w(n) + \left[ {\begin{array}{*{20}{c}}
\hfill 0 & \hfill 1 & \hfill 2 & \hfill 3 \\
\hfill { - 1} & \hfill 0 & \hfill 1 & \hfill 2 \\
\hfill { - 2} & \hfill { - 1} & \hfill 0 & \hfill 1 \\
\hfill { - 3} & \hfill { - 2} & \hfill { - 1} & \hfill 0 \\
\hfill { - 4} & \hfill { - 3} & \hfill { - 2} & \hfill { - 1} \\
\hfill { - 5} & \hfill { - 4} & \hfill { - 3} & \hfill { - 2} \\
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{e_0}} \\
{{e_1}} \\
{{e_2}} \\
{{e_3}} \\
\end{array}} \right]\]
用循环矩阵表示梯度向量的计算过程
\[\left[ {\begin{array}{*{20}{c}}
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
\hfill { - 5} & \hfill { - 4} & \hfill { - 3} & \hfill { - 2} & \hfill { - 1} & \hfill 0 & \hfill 1 & \hfill 2 & \hfill 3 \\
\hfill 3 & \hfill { - 5} & \hfill { - 4} & \hfill { - 3} & \hfill { - 2} & \hfill { - 1} & \hfill 0 & \hfill 1 & \hfill 2 \\
\hfill 2 & \hfill 3 & \hfill { - 5} & \hfill { - 4} & \hfill { - 3} & \hfill { - 2} & \hfill { - 1} & \hfill 0 & \hfill 1 \\
\hfill 1 & \hfill 2 & \hfill 3 & \hfill { - 5} & \hfill { - 4} & \hfill { - 3} & \hfill { - 2} & \hfill { - 1} & \hfill 0 \\
\hfill 0 & \hfill 1 & \hfill 2 & \hfill 3 & \hfill { - 5} & \hfill { - 4} & \hfill { - 3} & \hfill { - 2} & \hfill { - 1} \\
\hfill { - 1} & \hfill 0 & \hfill 1 & \hfill 2 & \hfill 3 & \hfill { - 5} & \hfill { - 4} & \hfill { - 3} & \hfill { - 2} \\
\hfill { - 2} & \hfill { - 1} & \hfill 0 & \hfill 1 & \hfill 2 & \hfill 3 & \hfill { - 5} & \hfill { - 4} & \hfill { - 3} \\
\hfill { - 3} & \hfill { - 2} & \hfill { - 1} & \hfill 0 & \hfill 1 & \hfill 2 & \hfill 3 & \hfill { - 5} & \hfill { - 4} \\
\hfill { - 4} & \hfill { - 3} & \hfill { - 2} & \hfill { - 1} & \hfill 0 & \hfill 1 & \hfill 2 & \hfill 3 & \hfill { - 5} \\
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
0 \\
0 \\
0 \\
0 \\
0 \\
{{e_0}} \\
{{e_1}} \\
{{e_2}} \\
{{e_3}} \\
\end{array}} \right]\]
重写滤波器系数更新公式如下：
\[\begin{array}{l}
w(n + 1) = w(n) + P_{L \times L}^{10}{C^T}{{\bf{e}}_L} = w(n) + P_{L \times L}^{10}{F^{ - 1}}F{C^T}{F^{ - 1}}F{{\bf{e}}_L} \\
W(n + 1) = W(n) + FP_{L \times L}^{10}{F^{ - 1}}F{C^T}{F^{ - 1}}E \\
\end{array}\]
这里W向量是滤波器系数的频域表示。同样做4步分解：
- 仍然是一个对解矩阵，可取循环矩阵C的第一列变换到频域后再取共轭来实现
- E是把前面补0后的误差信号变换到频域
- 这步是比较麻烦的地方。如果为了计算量忽略这一步（webrtc的做法），就是不对梯度做约束，而speex是通过变换到时域再置后面部分为0来避免直接矩阵运算的，个人比较喜欢speex的实现方式。当然，由于这是一个定值，不考虑计算量的话，也可以直接在频域进行矩阵的乘法
- 与W(n)在频域相加
至此，已经完成了滤波器分块的频域处理分析，但要注意的是，这里滤波器的长度是L = M + N – 1，分块处理转到频域虽然计算上方便了，但随着时域滤波器系数长度M和分块长度N越来越大，延时也会随之线性增大，接下来，就着手解决这个问题。
三、对滤波器进行分割（分段）
当滤波器长度M很大，且使用一个比M小很多的块长度时，可以把卷积和的运算过程分割为多个块的和，用多个分块滤波器的和来合成最终的结果，这样处理就可以使滤波器的延时大幅度的缩短。以M = 8，P = 2，N = 4为例进行说明这个分割过程。也就是说，如果对一个长度M = 8的滤波器，把滤波器分成2个块，每块4个数据，可以把卷积过程写为
\[\begin{array}{l}
y(n) = \sum\limits_{l = 0}^{P - 1} {{y_l}(n)} \\
{y_l}(n) = \sum\limits_{i = 0}^{N - 1} {{w_{i + lN}}x(n - lN - i)} \\
\end{array}\]
用矩阵的方式表示出来：
\[\left[ {\begin{array}{*{20}{c}}
\hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill | & \hfill { - 4} & \hfill { - 5} & \hfill { - 6} & \hfill { - 7} \\
\hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill | & \hfill { - 3} & \hfill { - 4} & \hfill { - 5} & \hfill { - 6} \\
\hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill | & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill { - 5} \\
\hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill | & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} \\
\end{array}} \right]\left[ {\frac{{\begin{array}{*{20}{c}}
{{w_0}} \\
{{w_1}} \\
{{w_2}} \\
{{w_3}} \\
\end{array}}}{{\begin{array}{*{20}{c}}
{{w_4}} \\
{{w_5}} \\
{{w_6}} \\
{{w_7}} \\
\end{array}}}} \right] \Rightarrow \begin{array}{*{20}{c}}
{\left[ {\begin{array}{*{20}{c}}
\hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} \\
\hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} \\
\hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} \\
\hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 \\
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{w_0}} \\
{{w_1}} \\
{{w_2}} \\
{{w_3}} \\
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
{{y_{\begin{array}{*{20}{c}}
0 & 0 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
0 & 1 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
0 & 2 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
0 & 3 \\
\end{array}}}} \\
\end{array}} \right]} \\
{\left[ {\begin{array}{*{20}{c}}
\hfill { - 4} & \hfill { - 5} & \hfill { - 6} & \hfill { - 7} \\
\hfill { - 3} & \hfill { - 4} & \hfill { - 5} & \hfill { - 6} \\
\hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill { - 5} \\
\hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} \\
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{w_4}} \\
{{w_5}} \\
{{w_6}} \\
{{w_7}} \\
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
{{y_{\begin{array}{*{20}{c}}
1 & 0 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
1 & 1 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
1 & 2 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
1 & 3 \\
\end{array}}}} \\
\end{array}} \right]} \\
\end{array}\]
\[\begin{array}{*{20}{c}}
{{y_0} = {y_{\begin{array}{*{20}{c}}
0 & 0 \\
\end{array}}} + {y_{\begin{array}{*{20}{c}}
1 & 0 \\
\end{array}}}} \\
{{y_1} = {y_{\begin{array}{*{20}{c}}
0 & 1 \\
\end{array}}} + {y_{\begin{array}{*{20}{c}}
1 & 1 \\
\end{array}}}} \\
{{y_2} = {y_{\begin{array}{*{20}{c}}
0 & 2 \\
\end{array}}} + {y_{\begin{array}{*{20}{c}}
1 & 2 \\
\end{array}}}} \\
{{y_3} = {y_{\begin{array}{*{20}{c}}
0 & 3 \\
\end{array}}} + {y_{\begin{array}{*{20}{c}}
1 & 3 \\
\end{array}}}} \\
\end{array}\]
接下来再分别把这两个分割出来的NxN矩阵块分别转换为循环矩阵的形式。（这里只写出第一个，第二个块有兴趣的朋友自己推吧）
\[\left[ {\begin{array}{*{20}{c}}
\hfill { - 3} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} \\
\hfill { - 2} & \hfill { - 3} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} \\
\hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 \\
\hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill 3 & \hfill 2 & \hfill 1 \\
\hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill 3 & \hfill 2 \\
\hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill 3 \\
\hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} \\
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{w_0}} \\
{{w_1}} \\
{{w_2}} \\
{{w_3}} \\
0 \\
0 \\
0 \\
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
\times \\
\times \\
\times \\
{{y_{\begin{array}{*{20}{c}}
0 & 0 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
0 & 1 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
0 & 2 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
0 & 3 \\
\end{array}}}} \\
\end{array}} \right]\]
这样就可以方便的转到频域进行处理，再把每块的处理结果相加，就完成了分段分块的频域卷积运算。滤波器的系数更新也是同理
另外，虽然2N-1个长度的频域复向量足以完成必要的卷积过程，实际算法中FFT长度通常取2N，这样计算更方便，这时分割后的第一个块的循环矩阵如下所示：
\[\left[ {\begin{array}{*{20}{c}}
\hfill { - 4} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} \\
\hfill { - 3} & \hfill { - 4} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} \\
\hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} \\
\hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 \\
\hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill 3 & \hfill 2 & \hfill 1 \\
\hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill 3 & \hfill 2 \\
\hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} & \hfill 3 \\
\hfill 3 & \hfill 2 & \hfill 1 & \hfill 0 & \hfill { - 1} & \hfill { - 2} & \hfill { - 3} & \hfill { - 4} \\
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{w_0}} \\
{{w_1}} \\
{{w_2}} \\
{{w_3}} \\
0 \\
0 \\
0 \\
0 \\
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
\times \\
\times \\
\times \\
\times \\
{{y_{\begin{array}{*{20}{c}}
0 & 0 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
0 & 1 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
0 & 2 \\
\end{array}}}} \\
{{y_{\begin{array}{*{20}{c}}
0 & 3 \\
\end{array}}}} \\
\end{array}} \right]\]
另一个分割块也是同理，这里不再详细列出来，有兴趣的朋友可以自己手画一遍玩玩，上图有一个小细节的哦，循环矩阵的对角元素可以是任意的，不影响最终效果，但通常大家都取前一个块的第一个元素。
剩下的活，就是前面转到频域的处理过程了，不再详述。
最后还有一个问题，是不是分割（段）越多越好，也不再详细分析了，直接给出结果：不是分割数P越大越好。由于分割过程改变了输入向量，也就改变了输入向量相关矩阵特征值的扩散度（条件数）。当P越大时，特征值的扩散度就越大，算法收敛就越慢，也更容易出现发散的可能。
参考文献：
- Advances in Network and Acoustic Echo Cancellation
- Adaptive Filters Theory and Applications Second Edition
- Adaptive Signal Processing Applications to Real-World Problems
