python实现LMS、NLMS、RLS、KALMAN等自适应滤波器

最新推荐文章于 2024-11-16 22:54:36 发布

ewan_xu

最新推荐文章于 2024-11-16 22:54:36 发布

阅读量9.8k

点赞数 23

分类专栏：回声消除文章标签：算法 webrtc python

本文链接：https://blog.csdn.net/myangel13141/article/details/108203186

版权

回声消除专栏收录该内容

4 篇文章

订阅专栏

由于使用C实现的自适应滤波器往往较为复杂，代码量较大，而python更加适合科学计算，使用python实现仅需几行代码就可以搞定，更加容易理解，本文使用简单的几行python代码实现8个自适应滤波器，包括lms,nlms,rls,kalman,
frequency domain adaptive filter,partition-block-based frequency adaptive filter,frequency domain kalman filter,partition-block-based frequency domain kalman filter），由简到难，方便学习自适应滤波器原理。

时域自适应滤波器

LMS

lms滤波器算法步骤如下：
假设真实的回声路径 $\boldsymbol h$ 长度为 $L$
$\boldsymbol h = [h_0 h_1 h_2 h_3...h_{L-1}]^T$
在n时刻输入的远端参考信号为：
$\boldsymbol x(n) = [x(n) x(n-1) x(n-2)... x(n-L+1)]^T$
估计的回声为：
$\hat y(n) = \boldsymbol x^T(n).\hat\boldsymbol h$
估计的回声和近端信号间的误差为：
$\hat y(n)$
系数更新公式：
$\hat\boldsymbol h_{n+1} = \hat\boldsymbol h_n + 2 \mu e^* (n) \boldsymbol x(n)$

python代码如下:

import numpy as np

def lms(x, d, N = 4, mu = 0.05):
    L = min(len(x),len(d))
    h = np.zeros(N)
    e = np.zeros(L-N)
    for n in range(L-N):
        x_n = x[n:n+N][::-1]
        d_n = d[n] 
        y_n = np.dot(h, x_n.T)
        e_n = d_n - y_n
        h = h + mu * e_n * x_n
        e[n] = e_n
    return e

NLMS

nlms相对于lms对输入信号做了归一化，系数更新公式变为：
$\hat\boldsymbol h_{n+1} = \hat\boldsymbol h_n + \frac{\mu e^* (n) \boldsymbol x(n)}{ \boldsymbol x^H(n) \boldsymbol x(n)}$

python代码如下:

import numpy as np

def nlms(x, d, N = 4, mu = 0.05):
    L = min(len(x),len(d))
    h = np.zeros(N)
    e = np.zeros(L-N)
    for n in range(L-N):
        x_n = x[n:n+N][::-1]
        d_n = d[n] 
        y_n = np.dot(h, x_n.T)
        e_n = d_n - y_n
        h = h + mu * e_n * x_n / (np.dot(x_n,x_n)+1e-8)
        e[n] = e_n
    return e

RLS

借用wikipedia上的图，rls滤波器算法步骤如下：
在这里插入图片描述

python代码如下：

import numpy as np

def rls(x, d, N = 4, lmbd = 0.999, delta = 0.0002):
    L = min(len(x),len(d))
    lmbd_inv = 1/lmbd
    h = np.zeros((N, 1))
    P = np.eye(N)/delta
    e = np.zeros(L-N)
    for n in range(L-N):
        x_n = np.array(x[n:n+N][::-1]).reshape(N, 1)
        d_n = d[n] 
        y_n = np.dot(x_n.T, h)
        e_n = d_n - y_n
        g = np.dot(P, x_n)
        g = g / (lmbd + np.dot(x_n.T, g))
        h = h + e_n * g
        P = lmbd_inv*(P - np.dot(g, np.dot(x_n.T, P)))
        e[n] = e_n
    return e

Kalman Filter

卡尔曼滤波器回声消除公式如下：:
$E_n = D_n - X_n^Th_{n-1}$
$P_n^- = P_{n-1} + Q$
$K_n = P_n^-X_n(X_nP_n^-X_n^T+R)^{-1}$
$\hat x_n = \hat x_{n-1} + K_nE_n$
$P_n = (I-K_nX_n^T)P_n^-$

对应的代码：

import numpy as np

def kalman(x, d, N = 64, beta=0.9,sgm2u=1e-2,sgm2v=1e-6):
    L = min(len(x),len(d))
    Q = np.eye(N)*sgm2v
    R = np.array([sgm2u]).reshape(1, 1)
    H = np.zeros((N, 1))
    P = np.eye(N)*sgm2v
    I = np.eye(N)

    e = np.zeros(L-N)
    for n in range(L-N):
        x_n = np.array(x[n:n+N][::-1]).reshape(1, N)
        d_n = d[n] 
        y_n = np.dot(x_n, H)
        e_n = d_n - y_n
        R = beta*R + (1-beta)*(e_n**2)
        Pn = P + Q                              
        K = np.dot(Pn, x_n.T) / (np.dot(x_n, np.dot(Pn, x_n.T)) + R)
        H = H + np.dot(K, e_n)
        P = np.dot(I - np.dot(K, x_n), Pn)
        e[n] = e_n
    return e

频域自适应滤波器

以上都是时域的自适应滤波器，我们可以通过OLS的方法进行分段然后使用FFT快速计算卷积，这样我们就可以在频域实现自适应滤波器

Frequency Domain Block NLMS Adaptive Filter

我们将NLMS算法基于这种思想进行实现，就得到了频域实现的NLMS自适应滤波器
公式如下：
在这里插入图片描述
对应代码：

import numpy as np
from numpy.fft import rfft as fft
from numpy.fft import irfft as ifft

def fdaf(x, d, M, mu=0.05, beta=0.9):
    H = np.zeros(M+1,dtype=np.complex)
    norm = np.full(M+1,1e-8)

    window =  np.hanning(M)
    x_old = np.zeros(M)

    num_block = len(x) // M
    e = np.zeros(num_block*M)

    for n in range(num_block):
        x_n = np.concatenate([x_old,x[n*M:(n+1)*M]])
        d_n = d[n*M:(n+1)*M]
        x_old = x[n*M:(n+1)*M]
 
        X_n = np.fft.rfft(x_n)
        y_n = ifft(H*X_n)[M:]
        e_n = d_n-y_n

        e_fft = np.concatenate([np.zeros(M),e_n*window])
        E_n = fft(e_fft)

        norm = beta*norm + (1-beta)*np.abs(X_n)**2
        G = mu*E_n/norm
        H = H + X_n.conj()*G

        h = ifft(H)
        h[M:] = 0
        H = fft(h)

        e[n*M:(n+1)*M] = e_n

    return e

Frequency Domain Kalman Filter

类似的，我们也可以得到频域卡尔曼自适应滤波器

import numpy as np
from numpy.fft import rfft as fft
from numpy.fft import irfft as ifft

def fdkf(x, d, M, beta=0.95, sgm2u=1e-2, sgm2v=1e-6):
    Q = sgm2u
    R = np.full(M+1,sgm2v)
    H = np.zeros(M+1,dtype=np.complex)
    P = np.full(M+1,sgm2u)

    window =  np.hanning(M)
    x_old = np.zeros(M)

    num_block = len(x) // M
    e = np.zeros(num_block*M)

    for n in range(num_block):
        x_n = np.concatenate([x_old,x[n*M:(n+1)*M]])
        d_n = d[n*M:(n+1)*M]
        x_old = x[n*M:(n+1)*M]

        X_n = np.fft.rfft(x_n)

        y_n = ifft(H*X_n)[M:]
        e_n = d_n-y_n

        e_fft = np.concatenate([np.zeros(M),e_n*window])
        E_n = fft(e_fft)

        R = beta*R + (1.0 - beta)*(np.abs(E_n)**2)
        P_n = P + Q*(np.abs(H))
        K = P_n*X_n.conj()/(X_n*P_n*X_n.conj()+R)
        P = (1.0 - K*X_n)*P_n 

        H = H + K*E_n
        h = ifft(H)
        h[M:] = 0
        H = fft(h)

        e[n*M:(n+1)*M] = e_n
  
    return e

分区块频域自适应滤波器

频域自适应滤波器每M点计算一次输出和更新滤波器，当抽头数量很大时，会造成比较大的延时，因此可以将滤波器进行分割，将长度为M的自适应滤波器分割为几个小块，这样每次计算滤波器的输入和输出只需要一个小块的长度，大大降低了延时。

Partitioned-Block-Based Frequency Domain Adaptive Filter (PFDAF)

通过将频域块NLMS自适应滤波器进行分块，就得到了也就是众所周知的webrtc AEC和speex中所使用的滤波器，代码相比较上面的滤波器比较复杂，这里就不贴了，感兴趣的同学可以去我的repository中查看：https://github.com/ewan-xu/pyaec

Partitioned-Block-Based Frequency Domain Kalman Filter

同样的我们也可以得到分块的频域卡尔曼滤波器，代码也比较复杂，感兴趣的同学可以去我的repository中查看： https://github.com/ewan-xu/pyaec

另外贴上分区块卡尔曼滤波器实现参考论文，感兴趣的同学可以参考参考：
F. Kuech, E. Mabande, and G. Enzner, “State-space architecture of the partitioned-block-based acoustic echo controller,”
in 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2014, pp. 1295-1299: IEEE

有任何疑问，欢迎加微信交流：xu19315519614