適応フィルタ

(2)

適応フィルタ・アルゴリズム

適応フィルタは、式(1)の解を逐次的に与えます。適応フィルタのアルゴリズムには、様々なものがありますが、ここでは、代表的な３つのアルゴリズムを簡単に紹介します。

LMS法

Least mean square (LMS) 法を正規化したN-LMS法は、最急降下法の近似であり、フィルタ更新式は、次式で与えられます。

(3)ここで、u(n) = [u(n), ..., u(n-M+1]^T はフィルタの入力ベクトル、Mは、フィルタ係数の数を表します。ξ(n) = d(n) - w_n-1^Tu(n) は事前誤差、μは学習の速度を制御するステップサイズパラメタ、α は0で割るのを防ぐための正則化定数です。

APA法

Affine Projection algorithm (APA)法は、NLMS法をブロック入力に拡張した方法であり、フィルタ更新式は、次式で与えられます。

(4)

ここで、ブロック入力は次式で定義されます。

(5)

事前誤差もベクトル化されています。

(6)

RLS法

Recursive least square(RLS)法は、式(1)の正規方程式における自己相関行列 R の逆行列を、次式を用いて逐次的に与えます。

(7)

ここで、

(8)

最適解は、次式で得られます。

(9)アルゴリズムは、これを基に、カルマンフィルタに似た形式に書き替えられています。

参考文献

[1] 浅野太、「音のアレイ信号処理」、コロナ社、2011年

適応フィルタモジュール

adaptiveFilter.py

import numpy as np

import numpy.linalg as lin

########################

# Class LMS

########################

class LMS:

def __init__(self,M,mu,alpha):

self.M = M

self.mu = mu

self.alpha = alpha

def fit(self,x,d):

M = self.M

mu = self.mu

alpha = self.alpha

K = x.size

w = np.zeros(M)

e = np.zeros(K)

for k in np.arange(M,K):

u = x[k:k-M:-1]

ep = d[k] - np.dot(w,u)

w = w + mu/(alpha+lin.norm(u)**2)*ep*u

e[k] = d[k] - np.dot(w,u)

self.w = w

self.e = e

def predict(self,x):

K = x.size

y = np.convolve(self.w,x)[0:K]

return y

########################

# Class APA

########################

class APA:

def __init__(self,M,L,mu,alpha):

self.M = M

self.L = L

self.mu = mu

self.alpha = alpha

def fit(self,x,d):

M = self.M

L = self.L

mu = self.mu

alpha = self.alpha

K = x.size

w = np.zeros(M)

e = np.zeros(K)

for k in np.arange(M+L,K):

U = np.zeros((M,L))

# Block input

for l in range(L):

u = x[k-l:k-l-M:-1]

U[:,L-l-1] = u

Q = lin.inv(alpha*np.eye(L) + np.dot(U.T,U))

# Prior Error

dd = d[k-L+1:k+1]

ep = dd - np.dot(w,U)

# Delta

delta = np.dot(np.dot(U,Q),ep)

w = w + mu*delta

u = x[k:k-M:-1]

e[k] = d[k] - np.dot(w,u)

self.w = w

self.e = e

def predict(self,x):

K = x.size

y = np.convolve(self.w,x)[0:K]

return y

########################

# Class RLS

########################

class RLS:

def __init__(self,M,c):

self.M = M

self.c = c

def fit(self,x,d):

M = self.M

K = x.size

P = (1/self.c) * np.eye(M)

w = np.zeros(M)

e = np.zeros(K)

for k in np.arange(M,K):

# Input vector

u = x[k:k-M:-1]

# Gain vector

Pu = np.dot(P,u)

g = Pu/(1 + np.dot(u,Pu))

# Prior error

ep = d[k] - np.dot(w,u)

# Update

w = w + g*ep

P -= np.outer(g,Pu.T)

# Error

e[k] = d[k] - np.dot(w,u)

self.w = w

self.e = e

def predict(self,x):

K = x.size

y = np.convolve(self.w,x)[0:K]

return y

サンプルプログラム

ここでは、適応フィルタの応用として、図２のような雑音除去問題を考えます。望みの応答には、信号s(n)に雑音v(n)が加わったものが観測されます。例えば、工場などの騒音のある場所で、マイクに向かってしゃべっている場合が、これに当たります。工場の中にある機械からの雑音は、インパルス応答 h(n) により変形されています。もし、工場の中の機械のそばに行って、雑音だけを収録できる場合は、これを適応フィルタの参照入力 u(n) として学習を行い、雑音成分を望みの応答から差し引くことで、信号を回復することができます。この場合、適応フィルタは、インパルス応答 h(n) を学習します。

adaptiveFilterExample.py

import numpy as np

import matplotlib.pyplot as plt

import numpy.linalg as lin

import scipy.stats as st

import adaptiveFilter as af

# Parameters

K = 30000 # データのサンプル数

T = 1000 # サイン波の周期

snr = 1.0 # SNR (信号対雑音比)

mu = 0.001 # ステップサイズパラメタ

alpha = 1e-3 # 正則化定数

M = 8 # フィルタの次数

L = 4 # [APA] アフィン変換の次数

c = 1 # [RLS] 初期化定数

figSize = (8,10)

# 観測値の生成

t = np.arange(K)

signal = np.sin( np.pi/T*t )

noise = st.norm.rvs(size=K)

h = np.array([0.0,1.0,-0.8,0.6,0,0,0,0]) # インパルス応答

d = np.convolve(h,noise)[0:K] + snr*signal

z = noise

plt.figure(figsize=(6,4))

plt.stem(h)

plt.figure(figsize=figSize)

plt.subplot(3,1,1), plt.plot(t,signal), plt.title('Signal')

plt.subplot(3,1,2), plt.plot(t,noise), plt.title('Noise')

plt.subplot(3,1,3), plt.plot(t,d), plt.title('Signal + Noise')

# 正規方程式

R = np.zeros((M,M))

r = np.zeros((M))

for k in range(M,K):

u = z[k:k-M:-1]

R += np.outer(u,u)

r += d[k]*u

w = np.dot(lin.inv(R),r)

y = np.convolve(w,z)[0:K]

plt.figure(figsize=figSize)

plt.subplot(3,1,1), plt.plot(t,d), plt.title('Signal + Noise')

plt.subplot(3,1,2), plt.plot(t,d-y), plt.title('Estimated signal')

# LMS

lms = af.LMS(M,mu,alpha)

lms.fit(z,d)

yLMS = lms.predict(z)

# APA

apa = af.APA(M,L,mu,alpha)

apa.fit(z,d)

yAPA = apa.predict(z)

# RLS

rls = af.RLS(M,c)

rls.fit(z,d)

yRLS = apa.predict(z)

# Graph

titleString = ['LMS','APA','RLS']

plt.figure(figsize=figSize)

gmax = 4.0

for n in range(3):

if n==0:

g = lms.e

elif n==1:

g = apa.e

elif n==2:

g = rls.e

plt.subplot(3,1,n+1), plt.plot(t,g)

plt.title(titleString[n])

plt.xlim(0,K)

plt.ylim(-gmax,gmax)

plt.figure(figsize=figSize)

gmax = 4.0

for n in range(3):

if n==0:

g = lms.e-signal

elif n==1:

g = apa.e-signal

elif n==2:

g = rls.e-signal

plt.subplot(3,1,n+1), plt.plot(t,g)

plt.title(titleString[n])

plt.xlim(0,K)

plt.ylim(-gmax,gmax)

plt.figure(figsize=figSize)

gmax = 1.2

for n in range(3):

if n==0:

g = d-yLMS

elif n==1:

g = d-yAPA

elif n==2:

g = d-yRLS

plt.subplot(3,1,n+1), plt.plot(t,g)

plt.title(titleString[n])

plt.xlim(0,K)

plt.ylim(-gmax,gmax)

plt.show()

上記のサンプルプログラム例では、信号 s(n)をサイン波、雑音をガウス雑音、インパルス応答をh(n)=[0,0, 1.0, -0.8, 0.6] (図３参照）として、学習を行います。図４は、s(n), v(n), d(n)をそれぞれ表します。図５は、望みの応答から、フィルタ出力を引いたe(n)を表します。学習が進むにつれ、雑音が除去され、信号のサイン波が現れてくるのがわかります。図６は、雑音の推定誤差を表します。LMS < APA < RLSの順で、収束が早くなっているのがわかります。これは、時刻nのフィルタ係数 w_n を更新するために用いられるデータのサンプル数に概ね比例する結果となっています。収束が早くなる代償として、計算量もこの順で大きくなります。なお、プログラムでは、結果がわかりやすくなるよう、様々な定数を調整してあります。定数の値によって、結果は変化します。