インパルス応答の測定

Last update: 2018/8/20

インパルス応答

システムにインパルス（時刻0にのみ値を持ち、あとは0となる信号）を入力した場合の応答を、インパルス応答(Impulse response)と呼びます。離散系では、図１のようなデルタ関数(Delta function) δ(n) を入力 x(n) とした場合の応答 y(n) をインパルス応答と呼びます。

(1)

インパルス応答は、以下のような用途で用いられます。

• • ソース信号と畳み込むことで、システムのシミュレーションを行う。

  • フーリエ変換やz変換を行うことにより、伝達関数を求める。

図１：インパルス応答

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)で時間軸に分散したエネルギーを元に戻すことです。図２はN=16384の場合のTSP信号の時間波形、図３は図２のスペクトログラムです。図３からわかるように、TSP信号は周波数が時間とともに変化するチャープ信号です。図４はTSP逆フィルタの時間波形です。TSP逆フィルタは、TSP信号を時間軸上で反転させたものです。

Kinectへの応用

インパルス応答の測定には、一般的に、同期したA/DとD/Aが必要となります。ここでは、Kinectに応用するため、簡易的な測定法を紹介します。Kinectに同期したD/Aを用意するのは難しいので、ここでは、PCなどからTSP信号を出力します。この場合、このあとの測定手順で述べるように、応答を閾値を使って切り出すことになるので、信号源からの絶対遅延時間は測定できません。

まず、TSP信号を平均回数の分だけ連結した信号を用意します。図４は平均回数４回の連結したTSP信号です（サンプルでは、"16384Repeat4.wav")。これをPCなどからスピーカを介して出力し、Kinectで収録します。収録には、Microsoft社のAudioCaptureRawを用います。図５は、収録例です。このTSP応答から、適当な閾値を用いて、応答長が N(=16384)の応答を切り出します。図５では、切り出し用のフレームを縦の点線で表しています。ここで大切なのは、平均回数＋１（＝５回）分の応答を切り出し、同期加算することです。これは、最後のTSP信号（４回目）の応答が、次のフレーム（第５フレーム）に入ってしまうためです。TSPは円状畳み込みを用いて設計されているので、第５フレームに入ってしまった応答を、第１フレームの先頭部分に足し合わせることで、つじつまがあうようになっています。図６は、同期加算後のTSP応答です。これにTSP逆フィルタ（サンプルの"16384inv.wav"）を円状畳み込みして、図７のようなインパルス応答を得ます。図８はこれを対数表示したものです。この図で、以下を確認します。

• インパルス応答がTSP長 N 以内に収束していること

• ノイズフロアが低く、SNRが所望の値になっていること

SNRについては、応用によりますが、40-60 dB程度確保されていれば、問題ないでしょう。

図9は、図8を逆２乗積分したものです。これに、単回帰分析を行い、60dB減衰した時間を調べることにより、残響時間RT60を知ることができます。この例では、残響時間は0.24秒です。図10は、インパルス応答の直接音の部分を切り出したもの、図11はこれをフーリエ変換し、位相応答を求めたものです。位相は、チャネル0で正規化してあります。これらから、音源の方向に対応した、時間進み/位相進みが生じているのがわかります。

サンプルプログラム

このプログラムは、上述の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()

参考文献

• 1. N. Aoshima, "Computer-nererated pulse signal applied for sound measurement," J. Acoust. Soc. Am., vol.69, no.5, pp.1484-1488 (1982)

  2. 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)