TSP 法によるインパルス応答の測定

インパルス応答の測定において、広く使われているのがTSP（Time Stretched Pulse、時間引き伸ばしパルス）法と呼ばれる手法です。本記事では TSP 法の原理と、Python で信号を生成して測定する流れを整理します。

インパルス応答の測定について

インパルス応答の測定方法として最も素朴なのは、単発のインパルス（理想的にはデルタ関数、実際にはパンッという短いクリック音やパルス）を 1 回再生し、マイクで録音した波形をそのままインパルス応答とみなす方法です。

一見よさそうですが、この方法には次のような問題があります。

• インパルスは 1 サンプルにエネルギーが集中するため、十分な SN 比を得るには非常に大きな瞬間音圧が必要になる
• スピーカーやアンプの許容入力を超えると、クリップや非線形歪みが発生し、測定結果がシステムの正しいインパルス応答ではなくなる
• 逆に音圧を抑えると、暗騒音（部屋の環境ノイズ）に埋もれて SN 比が悪化する

つまり「大きな音を出せば歪む、小さな音を出せばノイズに埋もれる」というジレンマがあり、単発インパルスによる測定は原理的に SN 比の上限が低いという課題があります。

TSP 法とは

TSP 法 [1] は、このジレンマを「エネルギーを 1 サンプルに集中させず、長い時間に薄く引き伸ばす」ことで解決する手法です。

TSP 信号は、周波数領域で見るとすべての周波数成分の振幅が等しい（オールパス）一方、位相だけが周波数の 2 乗に比例して変化する信号で、時間領域では低い周波数から高い周波数へ滑らかに掃引していく、振幅がほぼ一定のスイープ信号(下図)として現れます。

この信号を実際にスピーカーから再生すると、単発インパルスのような鋭いピークがないため、同じ最大音圧（スピーカーやアンプが安全に出せる上限）でもエネルギー全体を大きく増やすことができます。系全体のエネルギーが \( N \) 倍に増えれば SN 比はおおよそ \( 10 \times \log_{10}(N) \) [dB] 改善するため、\( N \) を数千〜数万点程度に取れば、単発インパルスでは実現できない SN 比での測定が可能になります。

さらに重要な性質として、TSP 信号は自分自身に対応する「逆フィルタ」を解析的に（数式だけで）作れるため、測定後の波形からインパルス応答を正確に取り出せるという利点があります。

TSP 信号の生成

TSP信号の生成手順は以下のようになります。

Aoshima [1] の基本形の TSP 信号は、周波数領域で次のように定義されます。信号長を \( N \)（偶数）、正の整数パラメータを \( m \) とすると、周波数インデックス \( k = 0, 1, \ldots, N/2 \) に対して

\[ H[k] = \exp\left(-j \, \dfrac{4\pi m k^2}{N^2}\right) \]

とし、残りの \( k = N/2+1, \ldots, N-1 \) は共役対称（\( H[N-k] = \overline{H[k]} \)）になるように埋めます。振幅 \( |H[k]| \) はすべての \( k \) で 1 なので、これは純粋な全域通過（オールパス）フィルタの周波数特性になっています。

\( H[k] \) の位相スペクトルのグラフは下図になります。

この \( H[k] \) を逆 FFT すれば、時間領域の TSP 信号 \( h_{tsp}[n] \) が得られます。

パラメータ \( m \) は掃引の速さを決めるもので、大きくするほど信号が周波数方向に「引き伸ばされ」ます。目安として \( m \) は \( N/4 \) 前後に取ることが多く、これによって信号のエネルギーがブロック長 \( N \) のほぼ全域に均等に分布し、振幅の変動が小さいスイープ信号になります。

TSP 法の測定手順

TSP法の測定手順は次のようになります。

• TSP 信号を対象のシステム（スピーカー→部屋→マイク、など）に入力し、出力 \( y[n] \) を録音する
• 録音した \( y[n] \) を FFT して \( Y[k] \) を求め、TSP 逆フィルタの伝達関数 \( H_{inv}[k] \) を掛ける
• \( Y[k] H_{inv}[k] \) を逆 FFT した結果がシステムのインパルス応答 \( h[n] \) になる

TSP 信号はオールパスなので、その逆フィルタの伝達関数は位相を反転させるだけで作れます。

\[ H_{inv}[k] = \overline{H[k]} = \exp\left(+j \, \dfrac{4\pi m k^2}{N^2}\right) \]

\( H[k] \times H_{inv}[k] = |H[k]|^2 = 1 \) となるため、\( Y[k] \) に \( H_{inv}[k] \) を掛けると TSP 信号自身の伝達特性がキャンセルされ、システムの伝達特性だけが残ります。これを逆 FFT すれば、理想的にはシステムのインパルス応答 \( h[n] \) がそのまま得られます。

なお同じ TSP 信号を複数回連続再生し、区間ごとに同期加算（コヒーレント平均）を取ると、システムに対して無相関なノイズの影響は平均回数の平方根に反比例して小さくなります。そのため実測ではノイズ環境に応じて数回〜数十回の平均を取るのが一般的みたいです。

プログラム

実際に Python で TSP 信号の生成、逆フィルタの伝達関数の作成、FFT と伝達関数の乗算によるインパルス応答抽出までを実装しました。

import numpy as np

def generate_tsp(n, m=None):
    """TSP信号 tsp と逆TSP信号 inv_tsp を生成"""
    if m is None:
        m = n // 4

    # 周波数インデックス k = 0, 1, ..., N/2 に対して位相スペクトルを計算
    freq = np.zeros(n, dtype=complex)
    k = np.arange(0, n // 2 + 1)
    phase = -4j * np.pi * m * (k ** 2) / (n ** 2)
    freq[:n // 2 + 1] = np.exp(phase)

    # 残りの周波数成分は共役対称になるように埋める
    freq[n // 2 + 1:] = np.conj(freq[1:n // 2][::-1])

    # 逆FFTでTSP信号を求める
    tsp = np.fft.ifft(freq).real

    # 位相を反転させた逆フィルタの周波数特性から逆TSP信号を求める
    inv_freq = np.conj(freq)
    inv_tsp = np.fft.ifft(inv_freq).real
    return tsp, inv_tsp

def extract_impulse_response(recorded, inv_tsp):
    """録音波形をFFTして逆フィルタの伝達関数を掛け、インパルス応答を抽出"""
    n = len(inv_tsp)

    # 録音波形と逆フィルタをそれぞれFFT
    y = np.fft.fft(recorded, n)
    hinv = np.fft.fft(inv_tsp, n)

    # 周波数領域で乗算してから逆FFT
    return np.fft.ifft(y * hinv).real

generate_tsp は前述の式通り周波数領域で TSP の位相スペクトルを作り、逆 FFT で TSP 信号 tsp と逆 TSP 信号 inv_tsp の両方を time-domain で返します。

測定時は tsp をスピーカーなどに再生し、マイクで録音した波形 recorded を extract_impulse_response に渡すことで、システムのインパルス応答が得られます。

さらに、WAVファイルとして実際に扱えるように soundfile での読み書きを追加した、全体のソースコードは以下です。

import numpy as np
import soundfile as sf

def generate_tsp(n, m=None):
    """TSP信号 tsp と逆TSP信号 inv_tsp を生成"""
    if m is None:
        m = n // 4

    # 周波数インデックス k = 0, 1, ..., N/2 に対して位相スペクトルを計算
    freq = np.zeros(n, dtype=complex)
    k = np.arange(0, n // 2 + 1)
    phase = -4j * np.pi * m * (k ** 2) / (n ** 2)
    freq[:n // 2 + 1] = np.exp(phase)

    # 残りの周波数成分は共役対称になるように埋める
    freq[n // 2 + 1:] = np.conj(freq[1:n // 2][::-1])

    # 逆FFTでTSP信号を求める
    tsp = np.fft.ifft(freq).real

    # 位相を反転させた逆フィルタの周波数特性から逆TSP信号を求める
    inv_freq = np.conj(freq)
    inv_tsp = np.fft.ifft(inv_freq).real
    return tsp, inv_tsp

def extract_impulse_response(recorded, inv_tsp):
    """録音波形をFFTして逆フィルタの伝達関数を掛け、インパルス応答を抽出"""
    n = len(inv_tsp)

    # 録音波形と逆フィルタをそれぞれFFT
    y = np.fft.fft(recorded, n)
    hinv = np.fft.fft(inv_tsp, n)

    # 周波数領域で乗算してから逆FFT
    return np.fft.ifft(y * hinv).real

def save_wav(path, data, fs):
    """信号をピーク正規化してWAVに保存"""
    peak = np.max(np.abs(data))
    if peak > 0:
        data = data / peak * 0.98
    sf.write(path, data, fs)

# パラメータ
N = 65536
M = N // 4    # 掃引パラメータ（目安として N/4）
fs = 48000    # サンプリング周波数

# TSP信号・逆TSP信号を生成してWAVに書き出す
tsp, inv_tsp = generate_tsp(N, M)
save_wav("tsp.wav", tsp, fs)
save_wav("tsp_inv.wav", inv_tsp, fs)

# tsp.wav を再生してマイクで録音した波形（recorded.wav）を読み込む
recorded, fs = sf.read("recorded.wav")

# インパルス応答を抽出してWAVに書き出す
h = extract_impulse_response(recorded, inv_tsp)
save_wav("impulse_response.wav", h, fs)

動作確認

pyroomacoustics（image-source 法による室内音響シミュレータ）で「正解」となる部屋のインパルス応答 \( h_{true}[n] \) を生成しました。5m × 4m × 3m の直方体の部屋（吸音率 0.25、反射次数 20）に音源とマイクを配置し、サンプリング周波数 48kHz でシミュレーションしました。

import pyroomacoustics as pra

N = 65536
M = N // 4    # 掃引パラメータ（目安として N/4）
fs = 48000    # サンプリング周波数
room = pra.ShoeBox([5.0, 4.0, 3.0], fs=fs, materials=pra.Material(0.25), max_order=20)
room.add_source([1.0, 1.0, 1.5])
room.add_microphone([3.5, 2.5, 1.5])
room.compute_rir()
h_true = room.rir[0][0]

save_wav("h_true.wav", h_true, fs)

tsp, inv_tsp = generate_tsp(N)
recorded = np.convolve(tsp, h_true)[:N]  # TSP を再生してマイクで録音した状況を模擬
h_recovered = extract_impulse_response(recorded, inv_tsp)

save_wav("recorded.wav", recorded, fs)
save_wav("impulse_result.wav", h_recovered, fs)

この \( h_{true}[n] \) と TSP 信号を畳み込んで「録音波形」を模擬し、extract_impulse_response で復元したインパルス応答と正解を比較したのが以下です。

シミュレーションの場合、スペクトログラムと波形についてはほとんど同じになることが確認できます。

おわりに

本記事では、TSP法を使用したインパルス応答の測定についてまとめました。インパルスを位相の操作で引き延ばして、SN 比の問題を解決するというのが面白い手法ですね。単純な方法ですが、自分だったら考えつかなさそう。

次は TSP法を使用して実際にインパルス応答を測定してみます。それまでに実環境の測定方法について調べておきます。

■参考文献
[1] N. Aoshima, “Computer-generated pulse signal applied for sound measurement,” J. Acoust. Soc. Am., vol. 69, no. 5, pp. 1484-1488, 1981.
[2] Y. Suzuki, F. Asano, H.-Y. Kim, T. Sone, “An optimum computer-generated pulse signal suitable for the measurement of very long impulse responses,” J. Acoust. Soc. Am., vol. 97, no. 2, pp. 1119-1123, 1995.