Swept-Sine(TSP)信号を使ったインパルス応答測定(python)

測定用信号を用いたインパルス応答測定

前の記事ではインパルス応答を測定するために拍手などのインパルス信号をつかいましたが、
実際にはある程度の長さがある測定用の信号を使って測定をします。
今回はその中のひとつであるSwept-Sine信号を使った方法に焦点をあてて解説をします。
実装する言語はpythonです。

ちなみに機器同士の同期や遅延の測定や調整などはしていないので、
複数位置で測定して波面を再現とかはできません。

Contents

1. 理論的なこと（っぽいもの）
2. 信号の生成
3. 逆フィルタによる処理
4. 実空間での測定
5. まとめ

1. 理論的なもの

まずSS (Swept-Sine) 信号とは何か、どうやって作るのか、どう使うのか
このパートで紹介したいと思います。

SS信号（もしくはTSP信号とは何か）

インパルス信号はある一点でのみエネルギーを持つ信号でした。（Fig. 1）

これを時間軸に伸長した（伸ばした）信号がTime-Stretched Pulse、
TSP信号と呼ばれるものです。
ちなみにTSPもSwept-Sineも同じものです。

単に時間信号で伸ばすと、Fig. 2のようになる感じがします。
ただ一体これはどんな音でしょう？

単に引き伸ばすだけでは説明不足で、時間とともに周波数があがる音です。
時間波形とスペクトログラムをFig. 3に示します。

Fig. 3に示したような信号をTSP信号だとかSS(Swept-Sine)信号だとか呼びます。
低い方から高い方へ掃引されたsin波という意味でSwept-Sineです。
時間で伸ばされてしまった以上、パルスじゃあないだろうということで、
最近はSwept-Sineの方が主流だそうです。

SS信号の作り方

次にSS信号の作り方です。以下にレシピを紹介します。

1. 周波数領域で徐々に周波数が上がるような数字列を作る
2. 逆FFTして時間領域にする
3. 実部だけ取り出す
4. （程よい位置になるようズラす）
5. 完成！

一番の肝は最初の部分でしょうか。
ただ既に数式として与えられているので、式をプログラムにするだけです。

\[ \rm{upSS}(k)=\left \{
\begin{array}{}
e^{-j2\pi J(k/N)^2} & k = 0, 1, \cdots , N/2 \\
\rm{upSS}(N-k)^* & k = N/2+1, \cdots , N-1 \\
\end{array} \right.
\] \[ \rm{upSS}: 時間とともに上昇するSweptSine信号の周波数成分, j: 虚数単位, N: 信号の長さ, J = 実効長, *: 複素共役
\]

SS信号の使い方

SS信号をインパルス応答測定に使う場合は、
まず再生をして、その音を録音するところからスタートです。

パルス信号を発生させた場合は録音した信号そのまま、
インパルス応答として扱うことができました。
しかしSS信号の場合は録音した音に対して少し処理が必要です。

その処理は、SS信号の逆関数を畳み込むことです。
SS信号の逆関数は（説明を省きますが）時間とともに周波数が降下するSS信号です。
先程はupSS信号と表現していたのでdownSS信号と呼びます。

\[ \rm{downSS}(k)=\left \{
\begin{array}{}
e^{j2\pi J(k/N)^2} & k = 0, 1, \cdots , N/2 \\
\rm{downSS}(N-k)^* & k = N/2+1, \cdots , N-1 \\
\end{array} \right.
\]

式を見ると指数の部分の符号が逆転しています。

これを畳み込むわけですが、
時間領域で畳み込みする方法と周波数領域で乗算する方法があります。
一般的に時間領域の畳み込みは時間がかかるので、後者をとります。

2. 信号の生成

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft 
from IPython.display import Audio 

N = 2**16
J = N/2
k = np.arange(0,N/2,1)
up_ss = np.zeros(N, dtype = np.complex)
up_ss[1:int(N/2+1)] = np.exp(-1j*2*np.pi*J*((k/N)**2))
up_ss_firsthalf = up_ss[1:int(N/2-1)]
up_ss[int(N/2+2):N] = up_ss_firsthalf[::-1].conjugate()

# 逆FFT
up_ss_ifft = ifft(up_ss, n=N)

# 実部取り出し
up_ss_ifft_real = np.real(up_ss_ifft)

# 円状シフト
shift_length = int((N-J)/2)
up_ss_ifft_real_roll = np.roll(up_ss_ifft_real, shift_length)

# 信号の可視化，再生
up_ss_signal = up_ss_ifft_real_roll

plt.plot(up_ss_signal)
plt.xlabel("Sample")
plt.ylabel("Amplitude")
plt.show()

fs = 44100
Audio(up_ss_signal, rate = fs)

こんな感じのコードでSS信号を作れます。
実行すると、Fig. 4と音データ（リンク先で再生できます）を得られます。

3. 逆フィルタによる処理

upSS信号を再生した音を録音したデータをインパルス応答にするために
downSS信号を畳み込みます。

ためしにupSS信号にdownSS信号を畳み込むことでインパルス信号が
現れるか確認したいと思います。

upSS信号の生成の後に以下のコードを実行すると確認ができます。

# down SS信号
down_ss = np.zeros(N, dtype = np.complex)
down_ss[1:int(N/2+1)] = np.exp(1j*2*np.pi*J*((k/N)**2))
down_ss_firsthalf = down_ss[1:int(N/2-1)]
down_ss[int(N/2+2):N] = down_ss_firsthalf[::-1].conjugate()

# 畳み込み
impulse_frequency = up_ss * down_ss
impulse_ifft = ifft(impulse_frequency)
impulse_ifft_real = np.real(impulse_ifft)

# 可視化
plt.plot(impulse_ifft_real)
plt.xlabel("sample")
plt.ylabel("Amplitude")
plt.show()

コードを実行するとFig. 5が得られます。
ある時刻(0 sample)でのみエネルギーが存在していること、
つまりインパルス信号であることが確認できました。

ここまでのコードをGoogle ColabとGithubに公開してあります。

4. 実空間での測定

以前はpythonのコードで再生と録音を同時におこなうこともしていましたが、
pyaudioがpython最新バージョンに対応していないなど面倒です。
したがって今回はwavファイルとしてもっているupSS信号を任意のソフトで
再生することと、録音することを （手動で）ほぼ同時におこないます。

録音した結果をFig. 6に示します。

これにdownSS信号を畳み込んだ結果、
つまりインパルス応答のグラフをFig. 7に示します。

時間軸方向に調整したもののグラフをFig. 8と、音データとして示します。

5. まとめ

今回はSwept-Sine信号を使ったインパルス応答測定をしました。

本文中では触れませんでしたが、測定用信号＋信号処理によって
インパルス応答測定をする手法はいくつか提案されています。
このSwept-Sine信号を使う手法は日本発ということで、日本の研究では
多く使われていて、一方海外だとMLS法が多いような気がします。

またインパルス応答自体、拡散音場で、線形時不変で…など色んな前提条件が
揃わないと厳密に正確な測定ができないことも知っておくといいかと思います。

---