最尤法
Last update: 2018/3/2
最尤法による音源定位
最尤法では、観測データの統計的モデルを作成し、このモデルが観測値と最も適合するよう、モデルのパラメタを最適化します。
観測ベクトル(各マイクロホンの観測信号をフーリエ変換し、ベクトルにしたもの)は次式のようにモデル化されると仮定します(観測値の線形モデル)。ここで、A(θ)は、音源のアレイマニフォールドベクトルであり、これには、音源の方向 θ がパラメタとして含まれます。
(1)
ここで、信号ベクトル sk および雑音ベクトル vk が多次元複素ガウス分布に従う、すなわち、sk 〜N(0,Γ) , vk 〜N(0,Σ)、と仮定すると、観測ベクトル zk も次式の多次元複素ガウス分布に従います(観測値の統計的モデル)。
(2)
R は観測値の共分散行列であり、式(1)を用いて、次式のようにモデル化されます。
(3)
ここで、Γ=E[sk sk H]、Σ=E[vk vk H] は、信号ベクトルおよび雑音ベクトルの共分散行列を表します。これらを用いて、対数尤度は、次式のように表されます。
(4)
Cは、次式で定義される標本共分散行列です。
(5)
最尤法では、(4)で示した尤度が最大となるよう、パラメタ θ を決定します。
サンプルプログラム
MLlocalization.py
このプログラムでは、音源数=1と音源数=2の場合について、対数尤度を求めるものです。図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()
図1:音源数=1の場合の対数尤度
図2(a):音源数=2の場合の対数尤度
図2(b):音源数=2の場合の対数尤度(3次元表示)