基于LMS的正弦信号去噪

# -*- coding: utf-8 -*- 
"""
Created on Thu May 30 21:17:42 2013
@author: Timchen525
"""
import numpy as np
import matplotlib.pyplot as plt
from numpy import shape
#定义向量的内积
def multiVector(A,B):
    C=np.zeros(len(A))
    for i in range(len(A)):
        C[i]=A[i]*B[i]
    return sum(C)

#取定给定的反向的个数
def inVector(A,b,a):
    D=np.zeros(b-a+1)
    for i in range(b-a+1):
        D[i]=A[i+a]
    return D[::-1]

#lMS算法的函数
def LMS(xn,dn,M,mu,itr):
    en=np.zeros(itr)
    W=[[0]*M for i in range(itr)]
    for k in range(itr)[M-1:itr]:
        x=inVector(xn,k,k-M+1)
        y=multiVector(W[k-1],x)
        en[k]=dn[k]-y
        W[k]=np.add(W[k-1],2*mu*en[k]*x) #跟新权重
    #求最优时滤波器的输出序列
    yn=np.inf*np.ones(len(xn))
    for k in range(len(xn))[M-1:len(xn)]:
        x=inVector(xn,k,k-M+1)
        yn[k]=multiVector(W[len(W)-1],x)
    return (yn,en)
if __name__=="__main__":
    #参数初始化
    itr=1000 #采样的点数
    M=5 #滤波器的阶数
    mu=3.0342e-005 #步长因子
    t=np.linspace(0,99,itr)
    xs=10*np.sin(0.5*t)
    xn1=np.random.randn(itr) #参考输入端的输入信号
    xn=np.add(xn1,xs) #原始输入端的信号为被噪声污染的正弦信号
    dn=xs #对于自适应对消器，用dn作为期望
    #调用LMS算法
    (yn,en)=LMS(xn,dn,M,mu,itr)
    #画出图形
    plt.figure(1)
    plt.plot(t,xn,label="$xn$")
    plt.plot(t,dn,label="$dn$")
    plt.xlabel("Time(s)")
    plt.ylabel("Volt")
    plt.title("original signal xn and desired signal dn")
    plt.legend()
    plt.figure(2)
    plt.plot(t,dn,label="$dn$")
    plt.plot(t,yn,label="$yn$")
    plt.xlabel("Time(s)")
    plt.ylabel("Volt")
    plt.title("original signal xn and processing signal yn")
    plt.legend()
    plt.figure(3)
    plt.plot(t,en,label="$en$")
    plt.xlabel("Time(s)")
    plt.ylabel("Volt")
    plt.title("error between processing signal yn and desired voltage dn")
    plt.legend()
    plt.show()