インパルス応答の測定
Last update: 2018/8/20
インパルス応答
システムにインパルス(時刻0にのみ値を持ち、あとは0となる信号)を入力した場合の応答を、インパルス応答(Impulse response)と呼びます。離散系では、図1のようなデルタ関数(Delta function) δ(n) を入力 x(n) とした場合の応答 y(n) をインパルス応答と呼びます。
(1)
インパルス応答は、以下のような用途で用いられます。
ソース信号と畳み込むことで、システムのシミュレーションを行う。
フーリエ変換やz変換を行うことにより、伝達関数を求める。
図1:インパルス応答
TSPを用いたインパルス応答の測定
ここでは、Time-stretched pulse(TSP)信号を用いたインパルス応答の計測方法について解説します。上述のようにデルタ関数は、t=0のみに振幅を持つため、信号のエネルギーが小さく、信号対雑音比(SNR)の高い測定が困難です。そこで、デルタ関数と同じ振幅で、SNRを向上させるため、TSP信号が考案されました(参考文献1,2)。TSP信号は、デルタ関数に、周波数の二乗に比例して位相が進む位相回転フィルタをかけたもので、周波数領域で、次式のように定義されます。
(2)
k: 離散周波数
M: TSPの収束を制御する定数(参考文献2)
N: TSP長(点数)
このTSP信号を信号源としたTSP応答を測定し、G(k)=1/H(k) で定義されるTSP逆フィルタをかけることで、インパルス応答を得ます。TSP逆フィルタの役割は、式(2)で時間軸に分散したエネルギーを元に戻すことです。図2はN=16384の場合のTSP信号の時間波形、図3は図2のスペクトログラムです。図3からわかるように、TSP信号は周波数が時間とともに変化するチャープ信号です。図4はTSP逆フィルタの時間波形です。TSP逆フィルタは、TSP信号を時間軸上で反転させたものです。
図2:TSPの時間波形
図3:TSPのスペクトログラム
図3:TSP逆フィルタの時間波形
Kinectへの応用
インパルス応答の測定には、一般的に、同期したA/DとD/Aが必要となります。ここでは、Kinectに応用するため、簡易的な測定法を紹介します。Kinectに同期したD/Aを用意するのは難しいので、ここでは、PCなどからTSP信号を出力します。この場合、このあとの測定手順で述べるように、応答を閾値を使って切り出すことになるので、信号源からの絶対遅延時間は測定できません。
まず、TSP信号を平均回数の分だけ連結した信号を用意します。図4は平均回数4回の連結したTSP信号です(サンプルでは、"16384Repeat4.wav")。これをPCなどからスピーカを介して出力し、Kinectで収録します。収録には、Microsoft社のAudioCaptureRawを用います。図5は、収録例です。このTSP応答から、適当な閾値を用いて、応答長が N(=16384)の応答を切り出します。図5では、切り出し用のフレームを縦の点線で表しています。ここで大切なのは、平均回数+1(=5回)分の応答を切り出し、同期加算することです。これは、最後のTSP信号(4回目)の応答が、次のフレーム(第5フレーム)に入ってしまうためです。TSPは円状畳み込みを用いて設計されているので、第5フレームに入ってしまった応答を、第1フレームの先頭部分に足し合わせることで、つじつまがあうようになっています。図6は、同期加算後のTSP応答です。これにTSP逆フィルタ(サンプルの"16384inv.wav")を円状畳み込みして、図7のようなインパルス応答を得ます。図8はこれを対数表示したものです。この図で、以下を確認します。
インパルス応答がTSP長 N 以内に収束していること
ノイズフロアが低く、SNRが所望の値になっていること
SNRについては、応用によりますが、40-60 dB程度確保されていれば、問題ないでしょう。
図9は、図8を逆2乗積分したものです。これに、単回帰分析を行い、60dB減衰した時間を調べることにより、残響時間RT60を知ることができます。この例では、残響時間は0.24秒です。図10は、インパルス応答の直接音の部分を切り出したもの、図11はこれをフーリエ変換し、位相応答を求めたものです。位相は、チャネル0で正規化してあります。これらから、音源の方向に対応した、時間進み/位相進みが生じているのがわかります。
図4:TSP信号を4回くりかえしたもの
図5:図4に対するTSP応答
図6:TSP応答を同期加算したもの
図7:インパルス応答
図8:図7の対数表示
図9:図8の逆2乗積分
図10: インパルス応答の直接音の部分
図11:チャネル間位相差
サンプルプログラム
このプログラムは、上述のKinectで収録したTSP応答から、インパルス応答を求めるプログラムです。
# モジュールのインポート
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as wav
import scipy.fftpack as ft
import scipy.stats as st
import sys
sys.path.append('../lib')
import kinectSensor as kin
# パラメータ - データにあわせて変更してください
dataRoot = '../Data/TSP/'
dataFile = 'D-45R100'
TSPresponseFile = dataRoot + 'Kinect/TSP_' + dataFile + '.wav' # TSP応答ファイル
dataType = 'int'
InverseTSPFile = dataRoot + "16384inv.wav" # TSP逆フィルタファイル
ImpulseResponseFile = dataRoot + 'Kinect/' + dataFile + '.wav' # インパルス応答ファイル
Threshold = 0.2 # 切り出しの閾値
FlameTopMargin = 4000 # 切り出しのマージン
NumberOfReptition = 4 # 平均回数
SearchOrigin = 10000 # 切り出しサーチの開始点
flagFigSave = False
# TSP応答の読み込み
z, Fs = kin.multiWavRead(TSPresponseFile,dataType)
(M,dataLength) = z.shape
# TSP応答のグラフ
fig = plt.figure()
plt.plot(z[0])
plt.title('TSP Response z[0]')
plt.xlim(0,len(z[0]))
dmax = np.max(np.abs(z[0,:]));
plt.ylim((-dmax*1.2,dmax*1.2))
plt.hlines(Threshold*dmax, 1, dataLength, colors='m', linestyles='dotted')
# Inverse TSP の読み込み
Fs, f = wav.read(InverseTSPFile);
TSPLength = len(f)
# 同期加算
for t in range(SearchOrigin,dataLength):
if np.abs(z[0,t]) > dmax*Threshold:
break;
Origin = t - FlameTopMargin;
y = np.zeros((M,TSPLength));
for k in range(NumberOfReptition+1):
FrameStart = Origin + k*TSPLength
FrameEnd = Origin + k*TSPLength + TSPLength
y = y + z[:,FrameStart:FrameEnd]
(dmin,dmax) = plt.ylim()
for k in range(NumberOfReptition+2):
FrameStart = Origin + k*TSPLength
plt.vlines(FrameStart,dmin,dmax,colors='g',linestyles='dotted')
if flagFigSave:
plt.savefig('Graph/TSPresponse.png',format='png')
# 同期加算のグラフ
fig = plt.figure()
plt.plot(y[0])
plt.xlim(0,len(y[0]))
plt.title('Averaged TSP Response y[0]')
if flagFigSave:
plt.savefig('Graph/AveragedResponse.png',format='png')
# TSP 逆フィルタ
F = ft.fft(f,TSPLength)
h = np.zeros((M,TSPLength))
for m in range(M):
Y = ft.fft(y[m],TSPLength)
H = F*Y;
h[m] = np.real( ft.ifft(H,TSPLength) );
# インパルス応答のグラフ
tm = np.arange(TSPLength)/Fs;
fig = plt.figure()
plt.plot(tm, h[0])
plt.title('Impulse Response h[0]')
plt.xlabel('Time [s]')
plt.xlim(min(tm),max(tm))
if flagFigSave:
plt.savefig('Graph/ImpulseResponse.png',format='png')
# インパルス応答のグラフ - 対数
dmax = np.max(np.abs(h[0]))
plt.figure()
plt.plot(tm,20*np.log10(np.abs(h[0])/dmax));
plt.ylabel('Power [dB]');
plt.xlabel('Time [s]')
plt.title('Log abs of h[0]');
plt.xlim(min(tm),max(tm))
if flagFigSave:
plt.savefig('Graph/LogImpulseResponse.png',format='png')
# インパルス応答のセーブ
np.save(ImpulseResponseFile,h)
# 逆二乗積分と残響時間
c = np.zeros(TSPLength)
for n in range(TSPLength):
c[n] = np.sum(h[0,n:TSPLength]**2)
c = 10*np.log10(c/c[0])
plt.figure()
plt.plot(tm,c)
plt.ylabel('Cumlative sum [dB]')
plt.xlabel('Time [s]')
plt.title('Reverse cumlative sum')
dmax = 5
dmin = -60
plt.ylim(dmin,dmax)
plt.xlim(min(tm),max(tm))
plt.grid()
decayPoint = [-15,-40] # 回帰分析の範囲
index = np.zeros(2,dtype=int)
for k in range(len(decayPoint)):
index[k] = np.argmin(np.abs(c - decayPoint[k]))
plt.plot(tm[index[k]], c[index[k]],'go')
a, b, r_value, p_value, std_err = st.linregress(tm[index[0]:index[1]], c[index[0]:index[1]])
gx = np.array([min(tm),max(tm)])
gy = a*gx + b # 回帰直線
plt.plot(gx,gy,'g--')
reverberationTime = -60/a # 残響時間 RT60
plt.text(max(tm)*0.8, -5, 'RT60=%.2f' % reverberationTime)
if flagFigSave:
plt.savefig('Graph/TSP_RT.png',format='png')
# チャネル間位相差
margin = 20
length = 64
NFFT = length
origin = np.argmax(np.abs(h[0])) - margin
g = np.zeros((M,length),dtype=float)
G = np.zeros((M,length),dtype=complex)
for m in range(M):
g[m] = h[m,origin:origin+length]
G[m] = ft.fft(g[m],NFFT)
Ref = np.copy(G[0]) # 0チャンネルで正規化
for m in range(M):
G[m] /= Ref
plt.figure() # インパルス応答の直接音の部分のグラフ
for m in range(M):
plt.subplot(2,2,m+1), plt.plot(g[m])
plt.title('Channel:'+str(m))
plt.tight_layout()
if flagFigSave:
plt.savefig('Graph/TSP_ImpulseDirect.png',format='png')
plt.figure() # 位相差のグラフ
Nfreq = int(NFFT/2)+1
freq = Fs/NFFT*np.arange(Nfreq)
for m in range(M):
plt.subplot(2,2,m+1), plt.plot(freq,np.angle(G[m][0:Nfreq]))
plt.title('Channel:'+str(m))
plt.ylim(-np.pi,np.pi)
plt.grid()
if m == 2 or m == 3:
plt.xlabel('Frequency [Hz]')
plt.tight_layout()
if flagFigSave:
plt.savefig('Graph/TSP_Phase.png',format='png')
plt.show()
参考文献
N. Aoshima, "Computer-nererated pulse signal applied for sound measurement," J. Acoust. Soc. Am., vol.69, no.5, pp.1484-1488 (1982)
Y. Suzuki, et al., "An optimum computer-generated pulse suitable for the measurement of very long impulse response," J. Acoust. Soc. Am., vol.97, no.2, pp.1119-1123 (1993)