前言
先说下目的：对在频域计算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步：
- 是一个对角矩阵，在编程实现中，可以取循环矩阵C的第一列变换到频域来代替
- 是把后面补0后的滤波器系数变换到频域
- 是一个逆FFT变换
- 是为了只选取最后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步分解：
- 仍然是一个对解矩阵，可取循环矩阵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个数据，可以把卷积过程写为
[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越大时，特征值的扩散度就越大，算法收敛就越慢，也更容易出现发散的可能。
参考文献：
- Advances in Network and Acoustic Echo Cancellation
- Adaptive Filters Theory and Applications Second Edition
- Adaptive Signal Processing Applications to Real-World Problems