インパルス応答の測定方法とは？【pythonで音響信号処理】

今回はインパルス応答の意義と取得方法について少しお話します。

応用の１つとしては音源に対してインパルス応答を畳み込むことで任意の響きを付加するなどがあります。

例えてみればリバーブのような形でコンサートホールの響き(インパルス応答)などを付加するような技術です。

本記事もpythonでの音源の読み込み・書き込みが前提となりますので，適宜以下の記事をご参照ください。

【pythonで音響信号処理】wavの仕組みと読み込み・書き出し

では早速本題です。

インパルス応答とは？

まずはインパルス応答について簡単にお話します。
※この辺りの基礎知識をご存知の方は読み飛ばしてもらうと良いかと思います

1つシステム(機構)を考えてみると，”入力$f$に対して何か処理をして出力$g$を返す”
という風に捉えることが出来るかと思います。

例えば自動販売機を考えると，お金を入れてボタンを押すというのが入力で，何か処理(入った金額を確認して飲み物を出す)をして，飲み物が落ちてくる(返ってくる)という感じでしょうか。

これは一般に以下のように図示されます。

ここで$g = Af$と表すことが出来ます。

こういった状況を考えて，入力$f$を時刻0でのみ無限大の値を持つ超関数であるDiracのデルタ関数$\delta(t)$とした出力$h(t) = A\delta (t)$がインパルス応答の定義です。

一般にインパルス応答という概念が用いられるのは(線形)時不変なシステムについてのみであり，これについてはシステムの情報を全て保持します。

線形時不変システム

勿論一般には作用素$A$は線形とは限らない作用素ですが，今回考える系(システム)については線形かつ時不変であるようなシステムです。

これを線形時不変(LTI: Linear Time Invariant)システムと呼びます。
定義を以下に示します。

線形性

$$A(\alpha f_{1}(t) + \beta f_{2}(t)) = \alpha Af_{1}(t) + \beta Af_{2}(t), \; \forall \alpha , \beta \in \mathbb{K}$$

時不変性

$$A f(t – \tau) = g(t – \tau)$$

感覚的には線形性は1次関数(直線)のような性質，時不変性はいつやっても同じ結果を返してくれる性質くらいに考えて貰っていいと思います。

先でインパルス応答の概念は時不変システムについて良く使われると言ったのは，ある時刻での応答を見ても時間が経ってみて測ったら全く違うものになっている(時変な)システムであれば殆どインパルス応答が意味を成さないからですね。

このような性質の良いシステムを考えて意味があるのかと言われるかもしれませんが，事実この世で起こりうる現象の幾つかは線形時不変システムの性質を持っています。

LTIシステムに現れる畳込み

LTIシステムにおいては作用がインパルス応答との畳み込みの形で書けることが知られています。
つまり

\begin{align}
g(t) &= A f(t) \\
&= \int_{-\infty}^{\infty} f(\tau) h(t-\tau) d\tau
\end{align}

これが今回のテーマの後ろ盾になっている訳です。
これを簡単に証明します。

\begin{align}
g(t) &= A f(t) \\
&= A \int_{-\infty}^{\infty} f(\tau) \delta(t-\tau) d\tau \\
&= \int_{-\infty}^{\infty} f(\tau) A \delta(t-\tau) d\tau \\
&= \int \int_{-\infty}^{\infty} f(\tau) h(t-\tau)
\end{align}

証明の1行目から2行目ではデルタ関数$\delta(t)$の性質を使っています。
更に2行目から3行目では線形性(斉次性)を，3行目から4行目では時不変性をそれぞれ使っています。

実装に於いては畳み込みの形で書けるというのは大きなメリットです。
なぜなら畳み込み定理より，フーリエ変換した先の周波数領域では掛け算の形になるからです。

$$\mathscr{F}[g(t)] = \mathscr{F}[(f \ast h) (t)] = \mathscr{F}[f(t)] \mathscr{F}[h(t)]$$

長い時間信号程，時間領域で積分(畳み込み演算)するよりFFT→掛け算の方が高速になります。

インパルス応答の取得方法

ここまでで「インパルス応答」の重要性がわかったかと思いますが，一体どうやって取得すれば良いのでしょうか。

名前の通りインパルス信号の応答ですので，最初に思いつくのはインパルス信号を鳴らしてその応答をマイク等で収録する方法です。

これは実際，以前は良く用いられていて「インパルス発生装置」なるものからインパルス信号(に近いもの)を発生させてそれの応答を収録してインパルス応答としていました。

しかしこの方法は実際に理想的なインパルス信号を出力出来るわけではなく，
一瞬で非常に大振幅の信号を発生する事自体が困難で，SN比(信号とノイズの比)の問題が大きいです。

更に収録側のゲインの調整や，聴取者の負担など考えると課題が残る部分があります。

そこで近年は信号処理的にインパルス応答を得る手法が幾つも提案されております。

代表的なものとしてはM系列(Maximum Length Sequence)信号を使用したもの，TSP(Time Stretched Pulse)信号を用いたもの等が挙げられます。
※日本ではTSPが主流(日本人が提案した手法であることも要因かもしれません)ですが，海外では殆どM系列信号によってインパルス応答を推定するようです

本記事では簡単にTSPを用いたインパルス応答測定方法を概説し，実際にコーディングしたいと思います。

TSP信号とは

TSP(Time Stretched Pulse)信号とはその名の通り時間を引き伸ばしたパルス信号です。

実際に生成して聴いてみて頂ければわかるのですが，低音から高音に徐々に周波数が上がっていくようなスイープ音になっています。

ここでTSPのイメージを掴むために少し幾何的な話をしましょう。

インパルス信号をフーリエ変換してみると$\mathscr{F}[\delta(t)] = \int_{-\infty}^{\infty} \delta(t) e^{-i \omega t} dt = 1$

となり全ての周波数成分を位相0で含んでいる事がわかります。つまり時刻0で全周波数の振幅最大が重なり合っている訳です。

これはインパルス信号がある意味で$\cos$の足し合わせで表現できることを表しています。
つまり下のようなイメージです。

$\cos$は周波数に依らずに$x=0$でピークを取るので，全ての周波数の$\cos$の足し合わせは$x=0$で無限に発散します。
その他では周波数によって振動するので無限個の足し合わせで消えてしまいます。

それに対してTSP信号は周波数領域で以下のように定義されます。

$$F(\omega) = e^{-i \omega^{2} t}$$

こちらは周波数が高くなるに連れて位相がずれていっていますね。これによってTSP信号はインパルス信号を時間方向に引き伸ばしたような信号になります。

下図はTSP信号をSTFTした結果になります。時間が経つにつれて周波数が上がっている様子がわかりますね。

この辺りをきちんと考えようと思うと群遅延という概念を扱う必要がありますが，今回は深入りしません。

TSP信号によってインパルス応答を推定する際には，TSP信号を再生しながらその応答(TSP応答)を収録し，TSP応答に対して逆TSP信号を畳み込むことで実現できます。

離散周波数領域におけるTSP信号

離散周波数領域でのTSP信号は以下のようになります。

$$F(k) =\begin{cases}
{e^{-i 2 \pi J(k / N)^{2}}} & {,k=0,1, \cdots, N / 2} \\
{F(N-k)^{*}} & {,k=N / 2+1, \cdots, N}
\end{cases}$$

ここで$N$は信号長，$J$は実効長(実際に信号の存在する長さ)になります。

このままで実行すると信号の開始が時刻0に設定されるため，時間シフトする場合はシフト量を$\tau$とする時

$$\hat{F}(k) = F(k) e^{-i 2 \pi (\tau / N)} ,k=0,1, \cdots, N / 2$$

のようにすることで時間領域で円状シフトできます。後半は同様に複素共役を取ればおっけーです。
ちなみに勿論時間領域信号にした後にシフトしても構いません。

TSP信号生成のためのpythonコード

以下がTSP信号を生成する関数です。

def get_tsp(order=18):
    N = 2**order + 2**(order-1) # length of signal
    J = 2**order # effective length
    shift = int((N-J)/2)
    TSP = np.exp(-1j * (2*np.pi) * J * (np.arange(int(N/2)) / N)**2)
    TSP = TSP * np.exp(-1j * (2*np.pi) * (shift/N) * np.arange(int(N/2)))
    iTSP = 1 / TSP
    tsp = np.fft.irfft(TSP)
    itsp = np.fft.irfft(iTSP)
    return tsp / np.max(np.abs(tsp)), itsp / np.max(np.abs(itsp))

このコードではnumpyのirfftという関数を使用することで$N/2$までの前半だけを作成して後半は勝手に複素共役のものを補完してもらっています。

実効長は$2^{order}$の長さに設定するようにしています。

また，実行してwavファイルで保存する際は以下のようなメインを書くと良いかと思います。

if __name__ == "__main__":
    order = 18
    tsp,itsp = get_tsp(order)
    dir_name = "."
    sampwidth = 3
    fs = 48000
    wf.writeWave("{}/tsp_{}.wav".format(dir_name, order), tsp, sampwidth=sampwidth, fs=fs)
    wf.writeWave("{}/itsp_{}.wav".format(dir_name, order), itsp, sampwidth=sampwidth, fs=fs)

「Audacity」によるTSP応答収録

「Audacity」はフリーのオーディオ編集ソフトです。

もし他に使い慣れたDAWなんかをお持ちの場合はそちらを使うことをオススメします。

ここではAudacityによるTSP応答測定の方法を簡単にご紹介しますが，オーディオインターフェイス等を使う場合や，細かい設定変更を行う場合は各自調べてみて下さい。

ではまずはダウンロードしましょう。macの場合はここからダウンロードできます。

ダウンロードしたらdmgファイルを実行してアプリケーションに追加します。

開いたら「ファイル→取り込み→オーディオの取り込み」から生成したTSP信号を取り込みます。

次にTSP応答のためのトラックを作ります。

上と同様にして「トラック→新しく追加→モノラルトラック」で追加します。
以下のような状態になったら大丈夫です。

後は録音の赤丸ボタンを押すだけです。

TSP信号の信号長を過ぎたところで停止ボタンを押しましょう。

上のように収録した波形が見れると思います。

僕のはお風呂場で録ったものなのですが，音が大きすぎてサチってしまってますね笑
今回は手順を説明するだけなので気にせず進めていきます。

最後に収録したトラックを選択して「ファイル→Export→選択したオーディオの書き出し」と進んでwavファイルを出力して下さい。

これでTSP応答を取得できました！

TSP応答をIRに変換しよう

最後に逆TSP信号を用いてTSP応答をインパルス応答に変換しましょう。

信号についてはどちらも今までの過程で持っているので後は畳み込みだけです。

ちょっと思い出すと時間領域の畳み込みは周波数領域で掛け算になっているのでした。

周波数領域の畳み込みが円状畳込みに注意すると以下のようなプログラムコードで実行できることが分かります。

def tspres2ir(tsp_res, itsp):
    N = itsp.shape[0]
    residual = tsp_res.shape[0] - N
    if residual &amp;amp;gt;= N:
        tsp_res[:N] = tsp_res[:N] + tsp_res[N:2*N]
    else:
        tsp_res[:residual] = tsp_res[:residual] + tsp_res[N:N+residual]
    TSP_RES = np.fft.rfft(tsp_res[:N])
    iTSP = np.fft.rfft(itsp)
    IR = TSP_RES * iTSP # convolution
    ir = np.fft.irfft(IR)
    return ir

実はインパルス応答測定では何回か連続でTSP応答を録って加算平均することでS/Nを大きくするような方法もあるのですが，今回は一発で録る方法を採用しました。

複数回録る場合については少しプログラムコードが変わる場合がありますのでご注意ください。
※実際には音場は完全な時不変システムでないため一発で録る場合も多いです

最後にメインのサンプルコードです。

if __name__ == "__main__":
    tsp_res_path = "./tspres_18.wav"
    ir_path = "./ir.wav"
    sampwidth = 3
    fs = 48000
    order = 18
    tsp_res = wf.readWave(tsp_res_path)
    tsp,itsp = get_tsp(order)
    ir = tspres2ir(tsp_res, itsp)
    ir = ir / np.max(ir)
    wf.writeWave(ir_path, ir, sampwidth, fs)

プロットしてみるとインパルス応答になっているのが分かるかと思います。

また，今回は信号全体をインパルス応答として出力していますが，実際には普通の室であればインパルス応答は長くてもせいぜい2秒程度かと思います。

信号を見てもっと短い所までを出力することでデータ量を減らすこともできますので必要に応じてやってみて下さい。