自适应滤波器分块分段技术详解

Posted 2021-03-15 icoolmedia

前言

　　先说下目的：对在频域计算FIR自适应滤波器，同时避免使用长滤波器时产生的大延迟的技术进行详细的分解。即要对滤波器分块又分段的过程进行细致的分析。这里假设读者有相应的LMS自适滤波器的基础

一、滤波器分析准备

　　先从时域LMS滤波器说起，设列向量

[{f{x}}(n) = {left[ {x(n),x(n - 1),...,x(n - M + 1)} ight]^T}]
[{f{w}}(n) = {left[ {{w_0}(n),{w_1}(n),...,{w_{M - 1}}(n)} ight]^T}]

　　这里列向量长度为M。考虑通过系数为w(n)的FIR滤波器对序列x(n)滤波，用卷积运算表示为[y(n) = sumlimits_{i = 0}^{M - 1} {{w_i}x(n - i)} ]

　　写成向量形式：

[y(n) = {{f{w}}^T}(n){f{x}}(n)]

　　简单从公式上，初学者其实不容易理解FIR滤波器的工作过程，这里换一种表达方式：把向量和矩阵用时间序列索引与符号来表示。应该就会比较好理解一些，对于序列：7，2，-3，-6，12，8，-7，-5，4，6。通过一个长度为4的FIR滤波器，输入向量 随时间序列的迭代变化表示为

[egin{array}{*{20}{c}}
{n = 0} hfill & {{f{x}} = {{left[ {egin{array}{*{20}{c}}
7 & 0 & 0 & 0
end{array}} ight]}^T}} hfill
{n = 1} hfill & {{f{x}} = {{left[ {egin{array}{*{20}{c}}
2 & 7 & 0 & 0
end{array}} ight]}^T}} hfill
{n = 2} hfill & {{f{x}} = {{left[ {egin{array}{*{20}{c}}
{ - 3} & 2 & 7 & 0
end{array}} ight]}^T}} hfill
{n = 3} hfill & {{f{x}} = {{left[ {egin{array}{*{20}{c}}
{ - 6} & { - 3} & 2 & 7
end{array}} ight]}^T}} hfill
{n = 4} hfill & {{f{x}} = {{left[ {egin{array}{*{20}{c}}
{12} & { - 6} & { - 3} & 2
end{array}} ight]}^T}} hfill
end{array}]

二、滤波器分块处理

　　下面准备写成时间序列索引的形式(向量各元素用时间序列对应的索引号代替，如：[0 -1 -2 -3]^T)，这样有利于在接下来的分析中看的更清楚一些，对于矩阵，也会沿用这样的方式。

　　先从频域矩阵分块处理说起，分块方法是一种批处理方法，为了说的更清楚一些，这里用一个示例来说明，设滤波器长度M = 6，块长度N = 4，表示每次批处理4个单元

[left[ {egin{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}} ight]left[ {egin{array}{*{20}{c}}
{{w_0}} hfill
{{w_1}} hfill
{{w_2}} hfill
{{w_3}} hfill
{{w_4}} hfill
{{w_5}} hfill
end{array}} ight] = left[ {egin{array}{*{20}{c}}
{{y_0}} hfill
{{y_1}} hfill
{{y_2}} hfill
{{y_3}} hfill
end{array}} ight]]

　　观察输入向量组成的矩阵可以发现，这是一个Toeplitz矩阵。要转到频域进行处理，一个比较可行的方式是把Toeplitz矩阵转为循环矩阵。再利用DFT把循环矩阵对角化。

　　那这个循环矩阵的大小是多少比较合适呢，由卷积过程可知，卷积后输出长度为L = M + N – 1 = 9，也就是说，用M + N – 1个卷积长度完全可以表达出来，也就是说循环矩阵C的大小为LxL应该是足够的。用循环矩阵表示如下：

[left[ {egin{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}} ight]left[ {egin{array}{*{20}{c}}
{{w_0}}
{{w_1}}
{{w_2}}
{{w_3}}
{{w_4}}
{{w_5}}
0
0
0
end{array}} ight] = left[ {egin{array}{*{20}{c}}
imes
imes
imes
imes
imes
{{y_0}}
{{y_1}}
{{y_2}}
{{y_3}}
end{array}} ight]]

　　这里输出向量中的X部分是循环卷积的结果，是要舍弃的部分。重写以上过程

 [{f{y}}(n) = P_{L imes L}^{01}C{f{hat w}}(n) = P_{L imes L}^{01}{F^{ - 1}}FC{F^{ - 1}}F{f{hat w}}(n)]

　　分解为4步：

1. 是一个对角矩阵，在编程实现中，可以取循环矩阵C的第一列变换到频域来代替
2. 是把后面补0后的滤波器系数变换到频域
3. 是一个逆FFT变换
4. 是为了只选取最后N个做为滤波器的输出结果，编程中做到这点很方便

　　再看滤波器的更新过程，这里把步长参数省去

[w(n + 1) = w(n) + left[ {egin{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}} ight]left[ {egin{array}{*{20}{c}}
{{e_0}}
{{e_1}}
{{e_2}}
{{e_3}}
end{array}} ight]]

　　用循环矩阵表示梯度向量的计算过程

[left[ {egin{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}} ight]left[ {egin{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}} ight]left[ {egin{array}{*{20}{c}}
0
0
0
0
0
{{e_0}}
{{e_1}}
{{e_2}}
{{e_3}}
end{array}} ight]]

　　重写滤波器系数更新公式如下： 

[egin{array}{l}
w(n + 1) = w(n) + P_{L imes L}^{10}{C^T}{{f{e}}_L} = w(n) + P_{L imes L}^{10}{F^{ - 1}}F{C^T}{F^{ - 1}}F{{f{e}}_L}
W(n + 1) = W(n) + FP_{L imes L}^{10}{F^{ - 1}}F{C^T}{F^{ - 1}}E
end{array}]

这里W向量是滤波器系数的频域表示。同样做4步分解：

 

1.  仍然是一个对解矩阵，可取循环矩阵C的第一列变换到频域后再取共轭来实现
2. E是把前面补0后的误差信号变换到频域
3. 这步是比较麻烦的地方。如果为了计算量忽略这一步（webrtc的做法），就是不对梯度做约束，而speex是通过变换到时域再置后面部分为0来避免直接矩阵运算的，个人比较喜欢speex的实现方式。当然，由于这是一个定值，不考虑计算量的话，也可以直接在频域进行矩阵的乘法
4. 与W(n)在频域相加

　　至此，已经完成了滤波器分块的频域处理分析，但要注意的是，这里滤波器的长度是L = M + N – 1，分块处理转到频域虽然计算上方便了，但随着时域滤波器系数长度M和分块长度N越来越大，延时也会随之线性增大，接下来，就着手解决这个问题。

三、对滤波器进行分割（分段）

　　当滤波器长度M很大，且使用一个比M小很多的块长度时，可以把卷积和的运算过程分割为多个块的和，用多个分块滤波器的和来合成最终的结果，这样处理就可以使滤波器的延时大幅度的缩短。以M = 8，P = 2，N = 4为例进行说明这个分割过程。也就是说，如果对一个长度M = 8的滤波器，把滤波器分成2个块，每块4个数据，可以把卷积过程写为

[egin{array}{l}
y(n) = sumlimits_{l = 0}^{P - 1} {{y_l}(n)}
{y_l}(n) = sumlimits_{i = 0}^{N - 1} {{w_{i + lN}}x(n - lN - i)}
end{array}]

　　用矩阵的方式表示出来：

[left[ {egin{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}} ight]left[ {frac{{egin{array}{*{20}{c}}
{{w_0}}
{{w_1}}
{{w_2}}
{{w_3}}
end{array}}}{{egin{array}{*{20}{c}}
{{w_4}}
{{w_5}}
{{w_6}}
{{w_7}}
end{array}}}} ight] Rightarrow egin{array}{*{20}{c}}
{left[ {egin{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}} ight]left[ {egin{array}{*{20}{c}}
{{w_0}}
{{w_1}}
{{w_2}}
{{w_3}}
end{array}} ight] = left[ {egin{array}{*{20}{c}}
{{y_{egin{array}{*{20}{c}}
0 & 0
end{array}}}}
{{y_{egin{array}{*{20}{c}}
0 & 1
end{array}}}}
{{y_{egin{array}{*{20}{c}}
0 & 2
end{array}}}}
{{y_{egin{array}{*{20}{c}}
0 & 3
end{array}}}}
end{array}} ight]}
{left[ {egin{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}} ight]left[ {egin{array}{*{20}{c}}
{{w_4}}
{{w_5}}
{{w_6}}
{{w_7}}
end{array}} ight] = left[ {egin{array}{*{20}{c}}
{{y_{egin{array}{*{20}{c}}
1 & 0
end{array}}}}
{{y_{egin{array}{*{20}{c}}
1 & 1
end{array}}}}
{{y_{egin{array}{*{20}{c}}
1 & 2
end{array}}}}
{{y_{egin{array}{*{20}{c}}
1 & 3
end{array}}}}
end{array}} ight]}
end{array}]

[egin{array}{*{20}{c}}
{{y_0} = {y_{egin{array}{*{20}{c}}
0 & 0
end{array}}} + {y_{egin{array}{*{20}{c}}
1 & 0
end{array}}}}
{{y_1} = {y_{egin{array}{*{20}{c}}
0 & 1
end{array}}} + {y_{egin{array}{*{20}{c}}
1 & 1
end{array}}}}
{{y_2} = {y_{egin{array}{*{20}{c}}
0 & 2
end{array}}} + {y_{egin{array}{*{20}{c}}
1 & 2
end{array}}}}
{{y_3} = {y_{egin{array}{*{20}{c}}
0 & 3
end{array}}} + {y_{egin{array}{*{20}{c}}
1 & 3
end{array}}}}
end{array}]

　　接下来再分别把这两个分割出来的NxN矩阵块分别转换为循环矩阵的形式。（这里只写出第一个，第二个块有兴趣的朋友自己推吧）

[left[ {egin{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}} ight]left[ {egin{array}{*{20}{c}}
{{w_0}}
{{w_1}}
{{w_2}}
{{w_3}}
0
0
0
end{array}} ight] = left[ {egin{array}{*{20}{c}}
imes
imes
imes
{{y_{egin{array}{*{20}{c}}
0 & 0
end{array}}}}
{{y_{egin{array}{*{20}{c}}
0 & 1
end{array}}}}
{{y_{egin{array}{*{20}{c}}
0 & 2
end{array}}}}
{{y_{egin{array}{*{20}{c}}
0 & 3
end{array}}}}
end{array}} ight]]

　　这样就可以方便的转到频域进行处理，再把每块的处理结果相加，就完成了分段分块的频域卷积运算。滤波器的系数更新也是同理

　　另外，虽然2N-1个长度的频域复向量足以完成必要的卷积过程，实际算法中FFT长度通常取2N，这样计算更方便，这时分割后的第一个块的循环矩阵如下所示：

[left[ {egin{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}} ight]left[ {egin{array}{*{20}{c}}
{{w_0}}
{{w_1}}
{{w_2}}
{{w_3}}
0
0
0
0
end{array}} ight] = left[ {egin{array}{*{20}{c}}
imes
imes
imes
imes
{{y_{egin{array}{*{20}{c}}
0 & 0
end{array}}}}
{{y_{egin{array}{*{20}{c}}
0 & 1
end{array}}}}
{{y_{egin{array}{*{20}{c}}
0 & 2
end{array}}}}
{{y_{egin{array}{*{20}{c}}
0 & 3
end{array}}}}
end{array}} ight]]

　　另一个分割块也是同理，这里不再详细列出来，有兴趣的朋友可以自己手画一遍玩玩，上图有一个小细节的哦，循环矩阵的对角元素可以是任意的，不影响最终效果，但通常大家都取前一个块的第一个元素。

　　剩下的活，就是前面转到频域的处理过程了，不再详述。

　　最后还有一个问题，是不是分割（段）越多越好，也不再详细分析了，直接给出结果：不是分割数P越大越好。由于分割过程改变了输入向量，也就改变了输入向量相关矩阵特征值的扩散度（条件数）。当P越大时，特征值的扩散度就越大，算法收敛就越慢，也更容易出现发散的可能。

参考文献：

1. Advances in Network and Acoustic Echo Cancellation
2. Adaptive Filters Theory and Applications Second Edition
3. Adaptive Signal Processing Applications to Real-World Problems