RでTSP信号を作る

TSP(Time Stretched Pulse)とはインパルス応答測定でよく用いられる信号で、音響機器の周波数特性やホールや部屋の残響測定を調べるときに使われます。matlabやscilabやOctaveで作るスクリプトは検索で出てきたのですが、Rのものは見つからなかったのでtuneRとsingalを使って作ってみました。私の検索力不足で、もしかするとすでに存在しているかもしれませんが・・・

#TSP信号について
TSP信号はサイン波のスイープパルスの一種で、伝達関数等の測定に用いられる信号で、波形は以下の式で表されます。(Optimized ATSP(OATSP)と言います)

\begin{equation}
  H\left( n \right) = \left\{ \begin{array}{l}
    a_0 \exp \left( \frac{j 4 \pi m n ^2}{N ^ 2} \right) & \left( 0 \leq n \leq \frac{N} {2} \right) \\
    a_0 \exp \left( \frac{-j 4 \pi m \left( N - n \right) ^2}{N ^ 2} \right) & \left( \frac{N} {2} - 1 \leq n \leq N - 1 \right)
  \end{array} \right. \\
\end{equation}

ここで$N$は波形全体の長さ(サンプル数)、$a_0$はTSPの振幅、mはTSPのパルス長を決めるパラメータで、$m=N/4$にするとTSPのパルス長は凡そ$N/2$になります。

#TSPによるインパルス応答の測定
具体的にTSPを用いてインパルス応答を測定する方法をざっくり述べると以下の通りです

1. TSP をスピーカーから再生する
2. 適当な場所にマイクを置いて、スピーカーからの再生音を録音する
3. SN比を稼ぐため、1. と 2. を繰り返して録音波形を平均する
4. 録音波形にTSPの逆信号を畳み込んでインパルス応答を求める

TSPの逆信号とは、上記$H(n)$の複素共役を取った波形になりますが、結果としてはTSP信号を時間軸で逆転させるとTSP逆信号になります、

#TSP信号のコーディング
RでTSPのWAVファイルを作成します。上記の式をそのままコーディングしました。

R

library(tuneR)
library(signal)

Fs <- 44100  #サンプリング周波数
bit <- 24    #ビット長
N <- 22100   #作成するwaveオブジェクト全体のサンプル数(0.5秒)
m = N / 4    #TSPパルスの長さ(約0.25秒)

# 0～N/2までの時間における位相(expのカッコ内)w1を作成
k1 <- c(0:(N / 2))
w1 <- 4 * pi * m * (k1 / N) ^ 2

# N/2～N-1までの時間における位相w2を同様に作成
k2 <- (N / 2 + 1) : (N - 1)
w2 <- -4 * pi * m * ((N - k2) / N) ^ 2

# w1とw2を結合させて0～N - 1までの位相wを作成
w <- append(w1, w2)
# 波形を作る
H <- exp(complex(re = 0,im = w))

# TSP信号を作る
# 逆FFTを計算
y <- fft(H, inverse = T)
# 実部のみを取得
y <- Re(y)
# 波形の始まりをずらす
y <- append(y[(m + 1) : length(y)], y[1 : m])
# 逆転させる
tsp <- rev(y)

#TSP逆信号を作る
# Hの複素共役をとる
H = Conj(H)
# 逆FFTを計算
y <- fft(H, inverse = T)
# 実部のみを取得
y <- Re(y)
# 波形の始まりをずらす
y <- append(y[(length(y) - m + 1) : length(y)], y[1 : (length(y) - m)])
# 逆転させる
tsp_inv <- rev(y)

# Waveオブジェクト作成
wav_tsp <- Wave(tsp, samp.rate = Fs, bit = bit)
wav_tsp_inv <- Wave(tsp_inv, samp.rate = Fs, bit = bit)
# 正規化（ノーマライズ)
wav_tsp <- normalize(wav_tsp, as.character(bit))
wav_tsp_inv <- normalize(wav_tsp_inv, as.character(bit))

作成した波形を表示

R

plot(wav_tsp)
plot(wav_tsp_inv)

インパルス応答測定の実際はまたの機会に。