频域LMS自适应滤波

作者: _wjchen_ | 来源:发表于2020-03-02 21:59 被阅读0次

    本文代码位于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-1y。 时域卷积,频域是相乘,因此有

    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算法做快速计算。 将卷积的输入xw都补零成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自适应滤波

        本文链接:https://www.haomeiwen.com/subject/llwikhtx.html