遅延時間推定
Last update: 2018/0825
時間差の推定
図１に示すように、信号源を２つのセンサで観測する場合、信号の伝搬距離に応じた時間差が生じます。これを式で書くと、次式のようになります。
(1)
ここでτ1、τ2は伝搬による遅延時間です。相対的な時間差
を用い、x(n) と y(n) に共通な遅延時間を省略すると、 (1)式は、次式のように簡略化して書けます。
(2)
この相対的な遅延時間τdは、信号源の位置情報など、有用な情報を含んでいます。このセクションでは、遅延時間推定の方法を解説します。
図１：遅延時間差の推定
一般化相互相関関数
遅延時間推定には、相互相関関数がよく用いられます。x(n) と y(n) の相互相関関数は、次式で定義されます。
(3)
(3)式と(2)式を比較するとわかるように、(2)式で生じた時間遅れを、(3)式では、"時間進み"のτを入れて、補正しようとしています。τ=τdとなったとき、時間遅れは時間進みで完全に相殺され、(3)式は E[s2(n)] となって、最大値をとります。
相互相関関数は、x(n) と y(n) のクロススペクトル
(4)
を用いて次式のように表されます。
(5)
ここで、X(k), Y(k)は、それぞれx(n) と y(n ) の離散フーリエ変換(DFT)を表します。"*"は複素共役を表します。期待値E[ ] は、時間(移動)平均で計算します。(5)式を自己相関関数について表したものは、ウィナー・ヒンチンの定理(Wiener–Khinchin theorem)として、知られています。
(5)式に周波数重みを掛け、一般化した次式は、一般化相互相関関数(generalized corss-correlation function, GCC)と呼ばれます。
(6)
重みの導入により、
雑音が重畳する場合、信号対雑音比(SNR)の高い帯域のみが相関関数の計算に寄与するようにすることができる
信号源に自己相関がある場合、信号源のスペクトルを白色化することにより、自己相関の影響を低減できる
などの利点があります。具体的な重み関数としては、次式のものが提案されています。
(7)
このうち、SCOT法では、(1)式のフーリエ逆変換の中身が、次式で定義されるコヒーレンス関数となります。
(8)
コヒーレンス関数は、x(n) と y(n ) の間の関連性を示す指標として用いられ、SNRを反映しています。このことからも、SNRが高い周波数帯域が選択的に用いられていることがわかります。
参考文献
金井浩、"音・振動のスペクトル解析," 第９章、第１０章、コロナ社 1999
クロススペクトル法のモジュール
crossSpectrum.py
ここでは、クロススペクトル法を用いて、GCCを計算するモジュールについて述べます。コンストラクタの引数を、表１に示します。クロススペクトルの時間平均の回数は、データ長と、フレーム長・フレームシフト幅から自動的に決定されます。コヒーレンス関数や、xとyの間の伝達関数も計算しているので、あとで必要に応じて取り出して使えるようになっています。
表１：コンストラクタの引数
関数GCC()では、式(7)に示した重みを選択して、GCCを計算します。Gxyを逆フーリエ変換すると、負の時間に対応する部分が、データ配列の後半に配置されるので、前半と後半を入れ替える操作を行っています。
# -*- coding: utf-8 -*-
"""
Created on Thu Aug 13 09:17:25 2015
@author: asano
Ver 1.2, 2015/08/26
"""
# Library
import numpy as np
import scipy.fftpack as ft
#############################
# Class: CSM
# x,y: time series data
#############################
class CSM:
def __init__(self,x,y,frameLength,frameShift,window):
# フレーム数の計算
Nsample = len(x)
NFFT = frameLength
Nframe = int((Nsample-(frameLength-frameShift))/frameShift)
print( '[CSM] Number of frames: ', Nframe )
# 窓関数
if window == 'Rect':
w = 1.0
elif window == 'Hamming':
w = np.hamming(NFFT)
elif window == 'Hanning':
w = np.hanning(NFFT)
else:
print( 'window function %s does not exist' % window )
return
# クロススペクトルと自己スペクトル
Gxx = np.zeros(NFFT,dtype=complex)
Gyy = np.zeros(NFFT,dtype=complex)
Gxy = np.zeros(NFFT,dtype=complex)
for n in range(Nframe):
t1 = frameShift*n
t2 = frameShift*n + frameLength
X = ft.fft(w*x[t1:t2],NFFT)
Y = ft.fft(w*y[t1:t2],NFFT)
Gxx += np.conj(X)*X
Gyy += np.conj(Y)*Y
Gxy += np.conj(X)*Y
Gxx = np.real(Gxx)
Gyy = np.real(Gyy)
# コヒーレンス関数
coherence = Gxy/np.sqrt(Gxx*Gyy)
# 伝達関数
transferFunction = Gxy/Gxx
self.Gxx = Gxx
self.Gyy = Gyy
self.Gxy = Gxy
self.NFFT = NFFT
self.coherence = coherence
self.transferFunction = transferFunction
# 一般化相互相関関数 GCC
def GCC(self,weight):
Gxx = self.Gxx
Gyy = self.Gyy
Gxy = self.Gxy
NFFT = self.NFFT
# 重み
if weight == 'FLAT':
Psi = 1.0
elif weight == 'ROTH':
Psi = 1.0/Gxx
elif weight == 'SCOT':
Psi = 1.0/np.sqrt(Gxx*Gyy)
elif weight == 'PHAT':
Psi = 1.0/np.abs(Gxy)
else:
print( 'weight %s does not exist' % weight )
return
# 逆フーリエ変換
rr = np.real(ft.ifft(Psi*Gxy,NFFT))
# 負の時間の部分を入れ替え
r = np.zeros(NFFT)
r[0:int(NFFT/2)-1] = rr[int(NFFT/2)+1:NFFT]
r[int(NFFT/2)-1:NFFT] = rr[0:int(NFFT/2)+1]
return(r)
サンプルプログラム
このサンプルでは、音声信号かランダム雑音を信号源 s(n) とし、これに20ポイントの時間差を与えて、GCCを計算します。観測系では、入出力に雑音が乗ることがよくあるので、この影響も見られるようになっています。図2は、雑音がない場合(SNRを100dB以上に設定）及びSNR=15dBに設定した場合の 観測波形 x(n), y(n) を示しています。
図３は、信号源 s(n) に音声波形を選び、雑音がない場合のGCCを示しています。(a)の重み=FLATは、通常の周波数重みなしの相互相関関数( 式(5) )を、(b)は重み関数にSCOT法を用いた場合のGCC( 式(6) )を示しています。周波数重みがない場合は、ラグτが20ポイント以外のところにもピークが観測されます。これは、音声波形が、図２に示すように、周期波形であるため、音声の周期の整数倍のところにも、ピークが生じるためです。一方、SCOT法を用いた場合は、Gxy が平坦化（白色化）されるため、音声の周期性の影響が取り除かれ、τ= 20 のところにだけピークが観測されています（音声の周期性が完全に除去されていないので、ちいさなピークは見えます）。このことから、周期性のある信号の時間差を測定するには、GCCが有効であることがわかります。
図４は、雑音がある場合のGCCを示しています。図５(b)は、雑音のみのスペクトル（緑）と、観測値 x(n) のスペクトル（青）を表しています。また、図５(a)は x(n) と y(n) のコヒーレンス関数を表しています。図５（b）でSNRの高い周波数においてで、コヒーレンス関数の値が大きくなっています。GCCのスコット法では、コヒーレンス関数を用いているため、SNRの高い周波数帯域が、選択的に、GCCの計算に寄与するようになっています。
# -*- coding: utf-8 -*-
"""
Created on Thu Jun 4 09:50:27 2015
@author: asano
Ver 1.1, 2015/08/13
"""
# Library
import numpy as np
import scipy.io.wavfile as wav
import matplotlib.pyplot as plt
import sys
sys.path.append('../lib')
import crossSpectrum as cs
# Parameter
NumberOfFrame = 4
frameLength = 512 # フレーム長（FFT長に等しい）
frameShift = frameLength
NFFT = frameLength
N = frameLength*NumberOfFrame # データのサンプル数
Tau = 20 # 時間差
SNR = 15.0 # 雑音を入れないときは100に設定
GCCweight = 'SCOT' # GCCに置ける重み関数
windowFunction = 'Hanning' # クロススペクトルに用いられる窓関数
flagFigSave = False
# 音源信号 (音声かランダム雑音かを選択)
#source = 'random'
source = 'speech'
if source == 'speech':
fileName = '../Data/Speech/speech11.wav'
Fs,data= wav.read(fileName)
Origin = 8500 # 音声の母音区間
elif source == 'random':
Origin = Tau
data = np.random.standard_normal(Origin+N)
Fs = 1
# 音源波形の図
if source == 'speech':
Ngraph = Fs
fig=plt.figure()
plt.plot(np.float64(data[1:Ngraph]))
(ymin,ymax) = plt.ylim()
plt.vlines(np.float(Origin),ymin,ymax)
plt.vlines(np.float(Origin+N),ymin,ymax)
# 時間差の挿入
x = np.float64(data[Origin:Origin+N])
y = np.float64(data[Origin-Tau:Origin-Tau+N])
amplitude = np.std(x) # SNR設定のため、振幅の正規化
x /= amplitude
y /= amplitude
# 入出力端に雑音を加算
if SNR < 100.0:
v = np.random.standard_normal((2,N))
gain = 10.0**(-SNR/20.0)
v *= gain
x = x + v[0,:]
y = y + v[1,:]
# x, yの時間波形
fig=plt.figure()
plt.subplot(2,1,1), plt.plot(x), plt.grid(), plt.xlim(0,N)
plt.subplot(2,1,2), plt.plot(y), plt.grid(), plt.xlim(0,N)
plt.xlabel('time [point]')
if flagFigSave:
fileName = 'Graph/observation_' + str(SNR) + '.png'
plt.savefig(fileName,format='png')
fig=plt.figure()
if source == 'speech':
Ngraph = 200
elif source == 'random':
Ngraph = 100
plt.subplot(2,1,1), plt.plot(x[0:Ngraph]), plt.grid()
plt.subplot(2,1,2), plt.plot(y[0:Ngraph]), plt.grid()
plt.xlabel('time [point]')
if flagFigSave:
fileName = 'Graph/observationMagnified_' + str(SNR) + '.png'
plt.savefig(fileName,format='png')
# クロススペクトル法によるGCC
cross = cs.CSM(x,y,frameLength,frameShift,windowFunction)
r = cross.GCC(GCCweight)
fig=plt.figure()
tm = np.arange(-int(NFFT/2)+1,int(NFFT/2)+1)
plt.plot(tm,r)
plt.xlim(-NFFT/2,NFFT/2)
plt.xlabel('time [point]')
plt.grid()
if flagFigSave:
fileName = 'Graph/GCC_' + str(SNR) + '_' + GCCweight + '.png'
plt.savefig(fileName,format='png')
# コヒーレンス関数
coherence = np.abs(cross.coherence[0:int(NFFT/2)+1])
f = np.arange(int(NFFT/2)+1)/float(NFFT)*float(Fs)
fig=plt.figure()
plt.plot(f,coherence)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Coherence')
if flagFigSave:
fileName = 'Graph/Coherence_' + str(SNR) + '.png'
plt.savefig(fileName,format='png')
# 周波数スペクトル
spectrum = 10.0*np.log10(np.abs(cross.Gxx[0:int(NFFT/2)+1]))
fig=plt.figure()
plt.plot(f,spectrum,label='Input')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Spectrum')
if SNR<100.0:
noiseCross = cs.CSM(v[0,:],v[0,:],frameLength,frameShift,windowFunction)
spectrum = 10.0*np.log10(np.abs(noiseCross.Gxx[0:int(NFFT/2)+1]))
plt.plot(f,spectrum,label='Noise')
if flagFigSave:
fileName = 'Graph/Spectrum_' + str(SNR) + '.png'
plt.savefig(fileName,format='png')
(a) 雑音なし
(b) SNR=15dB
図２：観測信号 x(n)
(a) 重み=FLAT
(b) 重み=SCOT
図３：雑音がない場合のGCC
(a) 重み=FLAT
(b) 重み=ROTH
(c) 重み=SCOT
(d) 重み=PHAT
図４：SNR=15dBの場合のGCC
(a) コヒーレンス関数
(b) 入力(青)と雑音(緑）のスペクトル
図５：コヒーレンス関数と入力のスペクトル