最尤法

Last update: 2018/3/2

最尤法による音源定位

最尤法では、観測データの統計的モデルを作成し、このモデルが観測値と最も適合するよう、モデルのパラメタを最適化します。

観測ベクトル(各マイクロホンの観測信号をフーリエ変換し、ベクトルにしたもの)は次式のようにモデル化されると仮定します（観測値の線形モデル）。ここで、A(θ)は、音源のアレイマニフォールドベクトルであり、これには、音源の方向 θ がパラメタとして含まれます。

(1)

ここで、信号ベクトル s_k および雑音ベクトル v_k が多次元複素ガウス分布に従う、すなわち、s_k 〜N(0,Γ) , v_k 〜N(0,Σ)、と仮定すると、観測ベクトル z_k も次式の多次元複素ガウス分布に従います（観測値の統計的モデル）。

(4)

Cは、次式で定義される標本共分散行列です。

(5)

最尤法では、(4)で示した尤度が最大となるよう、パラメタ θ を決定します。

サンプルプログラム

MLlocalization.py

このプログラムでは、音源数=1と音源数=2の場合について、対数尤度を求めるものです。図１と図２は、それぞれの場合の対数尤度です。

# 最尤推定法

# File name: MLlocalization.py

# Last update: 2018/3/2

# ライブラリ

import numpy as np

import scipy.signal as sig

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

import sys

sys.path.append('../lib')

import kinectSensor as kin

import spatialSpec as ss

# データファイル

observationFile = '../Data/Kinect/Speech0_Noise-45.wav' # 音源数=2のデータ

# パラメータ

dataType = 'int'

NFFT = 512 # FFTの点数

frameLength = NFFT

frameShift = int(frameLength/2)

BlockLength = 4000 # ブロック長

direction = np.arange(-90.0, 91.0, 2.0, dtype='float') # AMVの方向

Ndirection = len(direction)

SoundVelocity = 340.0 # 音速

prm = ss.Param() # 空間相関行列推定用のパラメタ

prm.NFFT = NFFT

prm.NumberOfSensor = 4

prm.FrameLength = frameLength

prm.FrameShift = frameShift

prm.BlockLength = BlockLength

prm.NumberOfSource = 1 # 音源数 1 or 2 .. 音源数を変えてみてください

print('Number of Frames: %d' % prm.getNumberOfFrame() )

if prm.NumberOfSource == 1:

blockNumber = 3 # 使用するブロック番号

freq = 1000 # Hz # 使用する周波数 [Hz]

elif prm.NumberOfSource == 2:

blockNumber = 28

freq = 1800 # Hz

BlockStart = BlockLength*blockNumber

BlockEnd = BlockLength*(blockNumber+1)

# Kinect観測信号のロード

z, Fs = kin.multiWavRead(observationFile,dataType)

(M,dataLength) = z.shape

print('M=%d, dataLength=%d' % (M,dataLength))

# 周波数のインデクス

freqIndex = int( (freq/Fs)*NFFT )

print('Frequency=%.1f' % (Fs/prm.NFFT*freqIndex) )

# 観測信号の表示

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

plt.plot(z[0,:])

(ymin,ymax) = plt.ylim()

plt.vlines(BlockStart,ymin,ymax,colors='g',linestyles='dashed')

plt.vlines(BlockEnd ,ymin,ymax,colors='g',linestyles='dashed')

# 空間相関行列

w = sig.hamming(prm.FrameLength) # 窓関数

R = ss.SpatialCorrelation(z[:,BlockStart:BlockEnd],w,prm)

# アレイ・マニフォールド

AA = kin.getAMatrix(direction,kin.sensorLocation,Fs,NFFT,SoundVelocity)

# 最尤推定法

if prm.NumberOfSource == 1:

Pml = np.zeros(Ndirection) # 最尤推定法の尤度

Pds = np.zeros(Ndirection) # 遅延和法による空間スペクトル

sigma_t = 1.0 # 信号の分散

sigma_n = 1e-4 # 雑音の分散

for n in range(Ndirection):

a = AA[freqIndex,n,:]

C = np.outer(a,np.conj(a))*sigma_t + sigma_n*np.eye(prm.NumberOfSensor)

Pml[n] = - np.trace( np.real(np.dot(np.linalg.inv(C),R[freqIndex,:,:])) )

Pds[n] = np.real( np.dot( np.dot(np.conj(a),R[freqIndex,:,:]), a ) )

plt.figure() # 最尤推定法のグラフ

plt.plot(direction,Pml)

plt.xlabel('Direction [$^\circ$]')

plt.ylabel('Likelihood')

plt.title('ML' )

plt.grid()

plt.figure() # 遅延和法のグラフ

plt.plot(direction,Pds)

plt.xlabel('Direction [$^\circ$]')

plt.ylabel('Spectrum')

plt.title('Delay-and-Sum' )

plt.grid()

elif prm.NumberOfSource == 2:

Pml = np.zeros((Ndirection,Ndirection)) # 最尤推定法の尤度

sigma_t = 1.0 # 信号の分散

sigma_n = 1e-4 # 雑音の分散

Gamma = sigma_t*np.eye(prm.NumberOfSource) # 信号の共分散行列

Sigma = sigma_n*np.eye(prm.NumberOfSensor) # 雑音の許分散行列

dmax = -1e50

for m in range(Ndirection):

for n in range(Ndirection):

A = np.zeros((prm.NumberOfSensor,prm.NumberOfSource),dtype=complex)

a1 = AA[freqIndex,m,:]

a2 = AA[freqIndex,n,:]

A = np.array([a1,a2]).T

C = np.dot( np.dot(A,Gamma), np.conj(A.T) ) + Sigma

Pml[m,n] = - np.trace( np.real(np.dot(np.linalg.inv(C),R[freqIndex,:,:])) )

if Pml[m,n] > dmax:

dmax = Pml[m,n]

idx1 = m

idx2 = n

print('Max at (%d %d) %e' % (direction[idx1], direction[idx2], Pml[idx1,idx2]) )

fig = plt.figure(figsize=(6,6)) # 2次元表示

GX,GY = np.meshgrid(direction,direction)

plt.pcolor(GX,GY,Pml, cmap=plt.cm.coolwarm)

plt.xlabel(r'$\theta_1$ [$^\circ$]')

plt.ylabel(r'$\theta_2$ [$^\circ$]')

plt.vlines(direction[idx1],min(direction),max(direction),colors='g',linestyles='dashed')

plt.hlines(direction[idx2],min(direction),max(direction),colors='g',linestyles='dashed')

fig = plt.figure(figsize=(6,6)) # 3次元表示

ax = fig.gca(projection='3d')

ax.plot_surface(GX, GY, Pml, rstride=1, cstride=1, cmap=plt.cm.coolwarm,

linewidth=0, antialiased=False)

ax.set_xlabel(r'$\theta_1$ [$^\circ$]')

ax.set_ylabel(r'$\theta_1$ [$^\circ$]')

ax.set_zlabel('Log likelihood')

plt.show()