前言
先说下目的:对在频域计算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
原文地址:https://www.cnblogs.com/icoolmedia/p/mdf_detailed.html