TSP信号の生成（Julia）

先週、Common Lisp版のOATSP信号を作りました。

marui.hatenablog.com

Julia版もあると便利かと思ってサクッと実装してみました（OATSPの説明などは上記のページを見てください）。一つの関数で測定信号と逆信号を作れると便利なので、まとめてしまっています。関数本体のコードは後に掲載するとして、まずは使い方です。

音ファイルの保存に必要になるパッケージなどを読み込みます。また、保存するときの標本化周波数は48000 Hzとします。

using SampledSignals, FFTW
using FileIO
import LibSndFile
fs = 48000;

測定信号と逆信号を作り、それぞれWAVファイルとして保存します。ただ、そのままだと振幅が小さいので-3 dBに正規化しています。また、PCM16Sampleのところで16 bitで保存するように指定しています。ここではNに4096、mに1200と、論文に書かれている信号長を指定しましたが、室内の残響時間を測定するときなどにはもっと長い信号を使うことになると思います。

N = 4096;
m = 1200;

# make measurement signal
x = make_oatsp(N, m);
xx = x ./ maximum(abs.(x)) * sqrt(1/2);
save("oatsp.wav", SampleBuf(PCM16Sample.(xx), fs));

# make inverse signal
y = make_oatsp(N, m, inverse=true);
yy = y ./ maximum(abs.(y)) * sqrt(1/2);
save("oatsp_inverse.wav", SampleBuf(PCM16Sample.(yy), fs));

正規化する前の測定信号と逆信号を環状畳み込みして、インパルスに戻るか確かめてみましょう。

using Plots

# check whether inverse works
z = real.(ifft(fft(x) .* fft(y)));
bar(z[1:50], legend=false)

1サンプル目だけに1が入っています。きれいに戻っていますね。

以下、関数本体です。

"""
    x = make_oatsp(N, m; inverse=false)

Makes an Optimum Aoshima's Time Stretched Pulse with specified
duration and signal width.

# References
* Nobuharu Aoshima. ``Computer-generated pulse signal applied for
  sound measurement.'' The Journal of the Acoustical Society of
  America, vol.69, no.5, pp.1484--1488. May 1981.
  https://doi.org/10.1121/1.385782 
* Yoiti Suzuki, Futoshi Asano, Hack-Yoon Kim, and Toshio Sone. ``An
  optimum computer-generated pulse signal suitable for the measurement
  of very long impulse responses.'' The Journal of the Acoustical
  Society of America, vol.97, no.2, pp.1119--1123. February 1995.
  https://doi.org/10.1121/1.412224
"""
function make_oatsp(N, m; inverse=false)
  p = π*m / (N/2)^2
  k = 0 : N÷2
  H = zeros(ComplexF64, N)
  if inverse == false
    H[1 : N÷2+1] =  exp.(1.0im * p * k .* k)
    H[N÷2+2 : N] = conj(reverse(H[2 : N÷2]))
  else
    H[1 : N÷2+1] =  exp.(-1.0im * p * k .* k)
    H[N÷2+2 : N] = conj(reverse(H[2 : N÷2]))  
  end
  Y = real.(ifft(H))

  # circular shift
  n_rot = N÷2 - m
  if inverse == false
    Y = [Y[n_rot+1:end] ; Y[1:n_rot]]
  else
    Y = [Y[end-n_rot+1:end] ; Y[1:end-n_rot]]
  end

  return(Y)
end