频域LMS自适应滤波

2020-03-02 本文已影响0人 _wjchen_

本文代码位于GitHub。

1. 概述

下图是一个典型的自适应滤波场景，输入信号 $x(n)$ 经过一个位置的系统变换 $h(z)$ 后得到参考信号 $d(n)$ 。自适应滤波器的目标是找出一组滤波器系数 $w(z)$ 来逼近系统 $h(z)$ ，使输入信号经过 $w(z)$ 变换后，与参考信号误差最小。

adaptive_filter.png

若使用FIR滤波器设计 $w(z)$ , 自适应滤波器的输出就是输入信号与滤波器权值的卷积：

$y(n)=\sum_{i=0}^{N-1}w_ix(n-i)$

自适应滤波器的算法就是以误差 $e(n)$ 为最小化目标，迭代求出最优的滤波器系数 $w(z)$ 。

$e(n)=d(n)-y(n)$

2. LMS和NLMS算法

LMS是最广泛应用的自适应滤波算法，以MSE误差为目标函数，以梯度下降为优化算法。并且通常情况下，LMS以最新的输入计算的瞬时梯度替代实际梯度计算，类似于机器学习的随机梯度下降法。

NLMS是使用输入的功率对步长进行归一化的方法，可以取得更好的收敛性能。

时域上实现LMS和NLMS算法的参考资料很多，这里不赘述，下面列出算法迭代步骤。

输入:

输入向量最新样本 $\boldsymbol{x}(n)$

期望输出最新样本 $\boldsymbol{d}(n)$

输出:

滤波器系数 $\boldsymbol{w}$ ,长度为M的FIR滤波器

滤波器输出 $\boldsymbol{y}(n)$

滤波器输出与期望间的误差 $e$

初始化:

滤波器系数 $\boldsymbol{w}(0)=zeros(1, M)$

迭代过程:

for n = 0, 1, 2...

1.读取输入样本 $\boldsymbol{x}(n)$ 和期望输出样本 $\boldsymbol{d}(n)$

2.滤波：

$\boldsymbol{y}(n)=\boldsymbol{w}^T(n)\boldsymbol{x}(n)$

3.计算误差：

$e(n)=\boldsymbol{d}(n) - \boldsymbol{y}(n)$

4.更新系数：

LMS:

$\boldsymbol{w}(n+1)=\boldsymbol{w}(n) + 2\mu e(n)\boldsymbol{x}(n)$

NLMS:

$\boldsymbol{w}(n+1)=\boldsymbol{w}(n) + \frac{\mu_0}{\boldsymbol{x}^T(n)\boldsymbol{x}(n)+\phi} e(n)\boldsymbol{x}(n)$

算法复杂度:

LMS: $2M+1$ 次乘法和 $2M$ 次加法

NLMS: $3M+1$ 次乘法和 $3M$ 次加法

3. Block LMS算法

LMS算法对输入数据是串行处理的，每输入一个样本，都需要进行 $2M+1$ 次乘法和 $2M$ 次加法，对于长度为 $N$ 的信号，计算时间复杂度为 $O(NM)$ 。可以通过将输入数据分段并行处理，并且利用频域FFT来做快速卷积，大大减少计算复杂度。

首先需要将串行的LMS算法转变为分块处理，也就是Block LMS(BLMS)。每次迭代，输入数据被分成长度为 $L$ 的块进行处理。和LMS使用瞬时梯度来进行滤波器参数更新不同，BLMS使用L点的平均梯度来进行参数更新。也就是机器学习里面Stochastic Gradient Descent 和 Mini-Batch Gradient Descent的区别。对第 $k$ 块数据，BLMS算法递推公式为：

$\boldsymbol{w}(k+1)=\boldsymbol{w}(k) + 2\mu_B \frac{\sum_{i=0}^{L-1}e(kL+i)\boldsymbol{x}(kL+i)}{L}$

输入:

输入向量 $\boldsymbol{x}$

期望输出 $\boldsymbol{d}$

输出:

滤波器系数 $\boldsymbol{w}$ ,长度为M的FIR滤波器

滤波器输出 $\boldsymbol{y}$

滤波器输出与期望间的误差 $\boldsymbol{e}$

初始化:

滤波器系数 $\boldsymbol{w}(0)=zeros(1, M)$

迭代过程:

for k = 0, 1, 2 ... N/L，读入第k块数据 $\boldsymbol{x}$ , $\boldsymbol{d}$

1. $\phi = zeros(1,L)$

2. for i = 0, 1, 2 ... L-1

2.1 滤波：

$\boldsymbol{y}(kL+i)=\boldsymbol{w}^T(k)\boldsymbol{x}(kL+i)$

2.2 计算误差：

$\boldsymbol{e}(kL+i)=\boldsymbol{d}(kL+i) - \boldsymbol{y}(kL+i)$

2.3 累计梯度：

$\phi=\phi+ \mu \boldsymbol{e}(kL+i)\boldsymbol{x}(kL+i)$

3 更新系数：

$\boldsymbol{w}(k+1)=\boldsymbol{w}(k) + \phi$

算法复杂度:

$2M+1$ 次乘法和 $2M+M/L$ 次加法

4. 快速卷积

前面一节描述的BLMS算法跟LMS算法相比，除了用块平均梯度代替瞬时梯度外，并没有不同。为了提升卷积计算的复杂度，我们需要引入快速卷积。也就是用FFT计算循环卷积，实现线性卷积的快速分块计算。

长度为L的 $x$ 和长度为M的 $w$ ，线性卷积的结果是长度为 $L+M-1$ 的 $y$ 。时域卷积，频域是相乘，因此有

$y(n)=x(n)*w(n)$

$\mathcal{FT}\left[y(n)\right]= \mathcal{FT}\left[x(n)\right] \mathcal{FT}\left[w(n)\right]$

傅立叶变换在频域上是连续的，计算机做傅立叶运算会将频域离散化采样，转化成 $N$ 点离散傅立叶变换(DFT)，并且用FFT算法做快速计算。将卷积的输入 $x$ 和 $w$ 都补零成 $N$ 点，做一个 $N$ 点的DFT，其逆变换就是循环卷积的结果，时域上是周期为 $N$ 的重复信号。

如上图所示，只有取DFT点数 $N>L+M-1$ ，才能防止卷积结果在时域上面混叠。IDFT结果的前 $M+L-1$ 个值就是所需的卷积结果。

于是两组有限长信号的卷积，则转换成3次DFT计算。

$y(n)= \mathcal{IDFT}\left\{\mathcal{DFT}\left[x(n)\right] \mathcal{DFT}\left[w(n)\right]\right\}$

直接计算卷积，大约需要 $N^2$ 次乘法和 $N^2$ 次加法。

采用FFT算法，3次DFT计算需要 $N+\frac{3N}{2}\log_2N$ 次复数乘法和 $3N\log_2N$ 次复数加法，相对直接计算算法复杂度从 $O(N^2)$ 降到了 $O(N\log_2N)$ 。

前面说的是两段有限长信号作一次卷积的流程。回到BLMS，对信号 $x$ 进行分段处理。每段输入长度是 $B$ ，滤波器长度是 $M$ 。

如果在 $w$ 后面补 $K-B$ 个零，而在 $x$ 当前块前部保留前一块的最后 $K-B$ 个点，做 $K$ 点快速卷积，结果中最后 $B$ 点是有效输出，而前 $K-B$ 点可丢弃。这种分块处理长卷积的方法，叫Overlap-and-Save方法。用FFT快速卷积+Overlap-and-Save方法，就可以高效处理BLMS滤波。

5. LMS算法的快速计算：频域自适应滤波

使用快速卷积的方法实现BLMS，就是频域实现的LMS自适应滤波器。称作Fast Block LMS (FBLMS) 或者 Frequency Domain Adaptive Filter (FDAF)。

对于长度为 $M$ 的滤波器，FBLMS一般采用 $2M$ 点FFT，且使用Overlap-and-Savee的快速卷积方法。也就是说，滤波器向量 $w$ 做FFT前须补 $M$ 个零值。

$W$ $= FFT\left[ \begin{matrix} w\\ 0 \end{matrix} \right]$

输入向量的分块长度也设为 $M$ ，则FFT的输入前 $M$ 点是保存前一块的数据，后 $M$ 点是当前块数据。

$FFT\left[ \begin{matrix} x_{k-1}\\ x_k \end{matrix} \right]$

使用overlap-save快速卷积方法实现BLMS中的卷积部分，当前块输出向量 $y$ 为IFFT的后 $M$ 点。

$\left[ \begin{matrix} C\\ y_k \end{matrix} \right]=IFFT\left[ FFT\left[ \begin{matrix} w\\ 0 \end{matrix} \right] FFT\left[ \begin{matrix} x_{k-1}\\ x_k \end{matrix} \right]\right]$

梯度的计算，可以将误差与输入信号放到频域来做，具体推导参考Block Adaptive Filters and Frequency Domain Adaptive Filters

$\left[ \begin{matrix} \phi\\ D \end{matrix} \right]=IFFT\left[ FFT\left[ \begin{matrix} 0\\ e \end{matrix} \right] FFT\left[ \begin{matrix} x_{k-1}\\ x_k \end{matrix} \right]'\right]$

文献[1]中提供的FBLMS框图

输入:

分块输入信号 $\boldsymbol{x_k}$ , 块长度为M

分块参考信号 $\boldsymbol{d_k}$ , 块长度为M

输出:

更新的滤波器系数 $\boldsymbol{w_k}$ ,长度为M的FIR滤波器

分块输出信号 $\boldsymbol{y_k}$ , 块长度为M

滤波器输出与参考信号间的误差 $\boldsymbol{e}$

初始化:

滤波器系数 $\boldsymbol{w}_0=zeros(M,1)$

初始数据块 $\boldsymbol{x}_0=zeros(M,1)$

迭代过程:

for $k = 1, 2 ... N/L$ ，读入第 $k$ 块数据 $\boldsymbol{w_k}$ , $\boldsymbol{d_k}$

1. $\boldsymbol{\phi} = zeros(M,1)$

2. 计算块输出

$\left[ \begin{matrix} C\\ y_k \end{matrix} \right]=IFFT\left[ FFT\left[ \begin{matrix} w_k\\ 0 \end{matrix} \right] FFT\left[ \begin{matrix} x_{k-1}\\ x_k \end{matrix} \right]\right]$

3. 计算误差： $\boldsymbol{e} = \boldsymbol{y}_k - \boldsymbol{d}_k$

4. 计算梯度：

$\left[ \begin{matrix} \phi\\ D \end{matrix} \right] =IFFT\left[ FFT\left[ \begin{matrix} 0\\ e \end{matrix} \right] \overline{FFT\left[ \begin{matrix} x_{k-1}\\ x_k \end{matrix} \right]}\right]$

5. 更新滤波器：

$FFT\left[ \begin{matrix} w_{k+1}\\ 0 \end{matrix} \right] = FFT\left[ \begin{matrix} w_k\\ 0 \end{matrix} \right] + \mu FFT\left[ \begin{matrix} \phi\\ 0 \end{matrix} \right]$

算法复杂度:

$10M\log M+26M$

def process(self, x_b, d_b, update=True): 
    '''
    Process a block data, and updates the adaptive filter (optional)
    
    Parameters
    ----------
    x_b: float
    the new input block signal
    d_b: float
    the new reference block signal
    update: bool, optional
    whether or not to update the filter coefficients
    '''  
    self.x = np.concatenate([self.x[self.B:], x_b])
    
    # block-update parameters
    X = fft(self.x)
    y_2B = ifft( X * self.W)
    y = np.real(y_2B[self.B:])
    
    e = d_b - y
    
    # Update the parameters of the filter
    if self.update:
    E = fft(np.concatenate([np.zeros(self.B), e])) # (2B)
    
    if self.nlms:
    norm = np.abs(X)**2 + 1e-10
    E = E/norm
    # Set the upper bound of E, to prevent divergence
    m_errThreshold = 0.2
    Enorm = np.abs(E) # (2B)
    # print(E)
    for eidx in range(2*self.B):
    if Enorm[eidx]>m_errThreshold:
    E[eidx] = m_errThreshold*E[eidx]/(Enorm[eidx]+1e-10) 
    
    # Compute the correlation vector (gradient constraint)
    phi = np.einsum('i,i->i',X.conj(),E) # (2B)
    phi = ifft(phi)
    phi[self.B:] = 0
    phi = fft(phi)
    
    self.W = self.W + self.mu*phi
    self.w = np.real(ifft(self.W)[:self.B]) 
    
return y, e

6. 减少时延：滤波器分割的频域自适应滤波

前面是将输入信号分块处理，提高算法效率。当FIR滤波器抽头数量很大时，FBLMS每M点计算一次输出和更新滤波器，造成比较大的延时。

一种想法是将滤波器也进行分割，这种改进延时的滤波器有几种名字：
Partitioned Fast Block LMS (PFBLMS)，Frequency Domain Adaptive Filter (PBFDAF)， Multi-DelayBlock Filter(MDF)，都是一回事。

conv.png

如果将长度为 $M$ 的滤波器 $w$ 等分为长度为 $B$ 的小段， $M=P*B$ 。则卷积的结果可以分解为 $P$ 个卷积之叠加。

$y(n)=\sum_{l=0}^{P} y_l(n)$

$y_l(n)=\sum_{i=0}^{B-1} w_{i+lB}x(n-lB-i)$

于是一段线性卷积被分解成 $P$ 个线性卷积，并且可以用FFT+OLS分别计算这 $P$ 个卷积。这样做好处是，每次迭代只需要输入长度为 $B$ 的信号块，保留原buffer中的后 $P-1$ 段，与最新的一段作为新的输入，就可以重复以上的 $P$ 段卷积叠加。每次迭代，可以更新 $B$ 点输出信号。

PFBLMS算法的时延就只有FBLMS的 $1/P$ ，极大地改善了滤波器的可用性。Speex和WebRTC的回声消除代码都使用了这种结构的滤波器。

文献[1]中提供的PFBLMS框图，其中M是我文中的B

输入:

分块输入信号 $\boldsymbol{x_k}$ , 块长度为B

分块参考信号 $\boldsymbol{d_k}$ , 块长度为B

输出:

滤波器系数 $\boldsymbol{w}$ ,长度为M=PB的FIR滤波器

分块输出信号 $\boldsymbol{y_k}$ , 块长度为B

滤波器输出与期望间的误差 $\boldsymbol{e}$

初始化:

滤波器系数 $\boldsymbol{w}_0=zeros(B,P)$

初始数据块 $\boldsymbol{x}_0=zeros(B,P)$

迭代过程:

for k = 0, 1, 2 ... N/B，读入第k块数据 $x_k$ , $d_k$

1. $\boldsymbol{\phi} = zeros(B,1)$

2. 滑动重组输入信号，并且计算FFT。

（实际只需计算最后一列，前P-1列可以使用上一次迭代保留结果）

$Xf_k = FFT\left[ \begin{matrix} x_{k-P}&...& x_{k-2} & x_{k-1}\\x_{k-P+1}&...& x_{k-1} & x_k \end{matrix} \right]$

2. 计算块输出

$\left[ \begin{matrix} C\\ y_k \end{matrix} \right]=IFFT\left[ FFT\left[ \begin{matrix} w_k\\ 0 \end{matrix} \right] Xf_k\right]$

3. 计算误差： $e = y_k - d_k$

4. 计算梯度：

$\left[ \begin{matrix} \phi \\ D \end{matrix} \right] =IFFT\left[ FFT\left[ \begin{matrix} 0\\ e \end{matrix} \right] \overline{Xf_k}\right]$

5. 更新滤波器：

$FFT\left[ \begin{matrix} w_{k+1}\\ 0 \end{matrix} \right] = FFT\left[ \begin{matrix} w_k\\ 0 \end{matrix} \right] + \mu FFT\left[ \begin{matrix} \phi\\ 0 \end{matrix} \right]$

算法复杂度:

与FBLMS相仿

def process(self, x_b, d_b, update=True): 
    '''
    Process a block data, and updates the adaptive filter (optional)
    
    Parameters
    ----------
    x_b: float
    the new input block signal
    d_b: float
    the new reference block signal
    update: bool, optional
    whether or not to update the filter coefficients
    '''
    self.x = np.concatenate([self.x[self.B:], x_b])
    Xf_b = fft(self.x)
    
    # block-update parameters
    self.Xf[1:] = self.Xf[:self.M-1] # Xf : Mx2B  sliding window
    self.Xf[0] = Xf_b # Xf : Mx2B  sliding window
    y_2B = ifft(np.sum(self.Xf * self.Wf, axis=0)) # [Px2B] element multiply [Px2B] , then ifft
    y = np.real(y_2B[self.B:])
    
    e = d_b - y
    
    if update:
    E = fft(np.concatenate([np.zeros(self.B), e])) # (2B)
    
    if self.nlms:
    norm = np.abs(Xf_b)**2 + 1e-6
    E = E/norm
    
    # Set the upper bound of E, to prevent divergence
    m_errThreshold = 0.2
    Enorm = np.abs(E) # (2B)
    # print(E)
    for eidx in range(2*self.B):
    if Enorm[eidx]>m_errThreshold:
    E[eidx] = m_errThreshold*E[eidx]/(Enorm[eidx]+1e-10)    
    
    
    # Update the parameters of the filter
    self.Wf = self.Wf + self.mu*E*self.Xf.conj()
    
    # Compute the correlation vector (gradient constraint)
    if self.constrained:
    waux = ifft(self.Wf, axis=1)
    waux[:, self.B:] = 0
    self.Wf = fft(waux, axis=1)
    
    self.w = np.real(ifft(self.Wf, axis=1)[:, :self.B].flatten())
return y, e

参考文献

[1] B. Farhang-Boroujeny, Adaptive Filters: theory and applications. John Wiley & Sons, 2013.

[2] Block Adaptive Filters and Frequency Domain Adaptive Filters

频域LMS自适应滤波

1. 概述

2. LMS和NLMS算法

3. Block LMS算法

4. 快速卷积

5. LMS算法的快速计算：频域自适应滤波

6. 减少时延：滤波器分割的频域自适应滤波

猜你喜欢

热点阅读