遅延時間推定

Last update: 2018/0825

時間差の推定

図１に示すように、信号源を２つのセンサで観測する場合、信号の伝搬距離に応じた時間差が生じます。これを式で書くと、次式のようになります。

(1)

ここでτ₁、τ₂は伝搬による遅延時間です。相対的な時間差

を用い、x(n) と y(n) に共通な遅延時間を省略すると、 (1)式は、次式のように簡略化して書けます。

(2)

この相対的な遅延時間τ_dは、信号源の位置情報など、有用な情報を含んでいます。このセクションでは、遅延時間推定の方法を解説します。

一般化相互相関関数

遅延時間推定には、相互相関関数がよく用いられます。x(n) と y(n) の相互相関関数は、次式で定義されます。

(3)

(3)式と(2)式を比較するとわかるように、(2)式で生じた時間遅れを、(3)式では、"時間進み"のτを入れて、補正しようとしています。τ=τ_dとなったとき、時間遅れは時間進みで完全に相殺され、(3)式は E[s²(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が高い周波数帯域が選択的に用いられていることがわかります。

参考文献

• 1. 金井浩、"音・振動のスペクトル解析," 第９章、第１０章、コロナ社 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')

図２：観測信号 x(n)

図３：雑音がない場合のGCC

図４：SNR=15dBの場合のGCC

図５：コヒーレンス関数と入力のスペクトル