基于LMS自适应滤波器

其实我也没有搞懂为什么在过滤的时候还需要期望信号
可能有参考数据吧
在很多应用中，S或D是可以已知的，例如在信道均衡应用中，先发送一个双方约定好的已知序列，接收方根据这个已知序列，通过自适应算法，求得信道的模型，之后再正常通信。

import numpy as np
import matplotlib.pyplot as plt
import scipy.io as scio


# 信号加噪
def awgn(x, snr):
    snr = 10 ** (snr / 10.0)
    xpower = np.sum(np.abs(x) ** 2) / len(x)
    npower = xpower / snr
    if type(x[0]) != np.complex128:
        return x + np.random.randn(len(x)) * np.sqrt(npower)
    else:
        return x + np.random.randn(len(x)) * np.sqrt(npower / 2) + 1j * np.random.randn(len(x)) * np.sqrt(npower / 2)


# % 输入参数:
# %   xn   输入的信号
# %   dn   所期望的响应
# %   M    滤波器的阶数
# %   mu   收敛因子(步长)
# % 输出参数:
# %   W    滤波器系数矩阵
# %   en   误差序列
# %   yn   滤波器输出
def lmsFunc(xn, dn, M, mu):
    itr = len(xn)
    en = np.zeros((itr, 1))
    W = np.zeros((M, itr))
    yn = np.ones(xn.shape) * np.nan
    for k in range(M, itr):
        if k==20:
            x = xn[k-1::-1]
        else:
            x = xn[k-1:k-M-1:-1]
        try:
            y = np.dot(W[:, k - 2], x)
            print(y)
        except:
            pass
        en[k-1] = dn[k-1] - y
        W[:, k-1] = W[:, k - 2] + 2 * mu * en[k-1] * x
        yn[k]=y;


    # for k in range(M, len(xn) ):
    #     if k == 20:
    #         x = xn[k - 1::-1]
    #     else:
    #         x = xn[k - 1:k - M - 1:-1]
    #     yn[k] = np.dot(W[:, -2], x)

    return yn, W, en


if __name__ == '__main__':
    fs = 1
    f0 = 0.02
    n = 1000
    t = np.arange(n)/fs
    xs = np.cos(2*np.pi*f0*t)
    xy = np.cos(2 * np.pi * f0 * t)
    ws = awgn(xs, 20)
    # data1 = scio.loadmat('xs.mat')
    # data2 = scio.loadmat('ws.mat')
    # xs = data1['xs'].flatten()
    # ws = data2['ws'].flatten()

    M = 20
    xn = ws
    dn = xy
    mu = 0.001
    # % 输入参数:
    # %   xn   输入的信号
    # %   dn   所期望的响应
    # %   M    滤波器的阶数
    # %   mu   收敛因子(步长)# % 输入参数:
    # # %   xn   输入的信号
    # # %   dn   所期望的响应
    # # %   M    滤波器的阶数
    # # %   mu   收敛因子(步长)
    # # % 输出参数:
    # # %   W    滤波器系数矩阵
    # # %   en   误差序列
    # # %   yn   滤波器输出
    # % 输出参数:
    # %   W    滤波器系数矩阵
    # %   en   误差序列
    # %   yn   滤波器输出
    yn, W, en = lmsFunc(xn, dn, M, mu)

    plt.figure()
    plt.subplot(211)
    plt.plot(t, ws)
    plt.subplot(212)
    plt.plot(t, yn)

    plt.figure()
    plt.plot(dn)
    plt.show()