Python實作自适應LMS濾波算法

1、生成資料

X = np.linspace(0,4*np.pi,itr,endpoint=True)
    Y = np.sin(X)
    signal_array = Y#[0.0]*noise_size
    noise_array =  np.random.normal(0, 0.3, noise_size)
    """noise_array = []
    for x in range(itr):
        noise_array.append(random.gauss(mu,sigma))"""
    signal_noise_array = signal_array+noise_array  
           

2、LMS算法

主要内容：

根據期望與實際輸出的誤差調節濾波權重值

期望：取滑動視窗内的均值作為期望。

def LMS(xn,dn,M,mu,itr):
    en=sc.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)
        d= x.mean()
        y=multiVector(W[k-1],x)
        en[k]=d -y
        W[k]=np.add(W[k-1],2*mu*en[k]*x) #跟新權重
    
    #求最優時濾波器的輸出序列
    
    yn=sc.inf*sc.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)

           

3、完整源碼

# -*- coding: utf-8 -*-
"""
Created on Wed Mar 27 12:57:42 2019

@author: Administrator
"""


import numpy as np

import matplotlib.pyplot as plt
import random
import scipy as sc


#定義向量的内積

def multiVector(A,B):
    
    C=sc.zeros(len(A))
    
    for i in range(len(A)):
        
        C[i]=A[i]*B[i]
        
    return sum(C)


#取定給定的反向的個數

def inVector(A,b,a):
    
    D=sc.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=sc.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)
        d= x.mean()
        
        y=multiVector(W[k-1],x)
        
        en[k]=d -y
        
        W[k]=np.add(W[k-1],2*mu*en[k]*x) #跟新權重
    
    #求最優時濾波器的輸出序列
    
    yn=sc.inf*sc.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=10000 #采樣的點數
    mu =0
    sigma =0.12
    noise_size  = itr
    
    X = np.linspace(0,4*np.pi,itr,endpoint=True)
    Y = np.sin(X)
    signal_array = Y#[0.0]*noise_size
    noise_array =  np.random.normal(0, 0.3, noise_size)
    """noise_array = []
    for x in range(itr):
        noise_array.append(random.gauss(mu,sigma))"""
    signal_noise_array = signal_array+noise_array  
    
    M=64 #濾波器的階數
    
    mu=0.0001 #步長因子
    
   
    
    xs=signal_noise_array
    
    
    xn=xs #原始輸入端的信号為被噪聲污染的正弦信号
    
    dn=signal_array #對于自适應對消器，用dn作為期望
    
    #調用LMS算法
    
    (yn,en)=LMS(xn,dn,M,mu,itr)
    
    #畫出圖形
    
    plt.figure(1)
    
    plt.plot(xn,label="$xn$")
    
    plt.plot(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(xn,label="$xn$")
    #plt.plot(xn,label="$xn$")
    plt.plot(dn,label="$dn$")
    plt.plot(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(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()
    

           

4、運作結果