トップ > LISP > TSP信号の生成

TSP信号の生成

LISP 音 研究

インパルス応答

頭部伝達関数を測定したり、スピーカーの応答を測ったり、室内音響指標の計算に使ったり……インパルス応答が測定できると色々と便利です。すでにMatlab、R、Juliaなどでインパルス応答測定のプログラムは書いてあるのですが、勉強中のCommon Lispで行うのが目下の目標です。

まずはインパルス応答を測定するための信号の生成を行います。今日は生成だけで、ファイルへの保存や実際にオーディオデバイスから出力するところなどは後々やります。

代表的な測定信号にはいくつかの種類がありますが（たとえば「PinkTSP信号をJuliaで作る - 丸井綜研」）、今日はシンプルにOATSP (Optimum Aoshima's Time-Stretched Pulse) を作ってみます。参考文献は次の通り。今回はSuzukiらの論文をもとに実装していきます。

  • 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

OATSPの数式

AoshimaもSuzukiも、周波数領域で測定信号Hを設計し、逆フーリエ変換をして時間領域にもってくるという方法です。時間領域で測定信号を測定対象に通して得られた結果と、測定信号の逆関数を畳み込むことで、測定対象のインパルス応答が得られます。Aoshimaの方法だとサブバンド信号のエイリアシングが生じるということで、それの対策をしたのがSuzukiの方法ですが、小難しいところはさておき、数式は以下のようになっています。

  • Nは全体の信号長（サンプル数）
  • mはそのうち信号の占める幅に関係するパラメータ（整数で通常はN/2未満）
  • p = (m π) / (N/2)2
  • H(k) = exp(j p k2), 0<=k<=N/2
  • H(k) = conj(H(N-k)), N/2+1 <= k < N

測定値は逆関数も必要で、次のような数式です。

  • H(k) = exp(-j p k2), 0<=k<=N/2
  • H(k) = conj(H(N-k)), N/2+1 <= k < N

Common Lispで試す

これをそのままCommon LispのREPLで（関数化はせずに）書いてみました。精度は高めでいきたいのでdouble-floatのcomplexを使います。ここではmの値は適当に決めていますが、関数化するときには引数にしたいと思います。

(defconstant j #C(0.0d0 1.0d0))         ; imaginary unit

(defparameter *N* (expt 2 4))           ; signal length in samples

(defparameter m (1- (/ *N* 2)))    ; tsp width (must be less than N/2)
(defparameter p (/ (* m pi) (expt (/ *N* 2) 2)))

(let ((H (make-array *N* :element-type 'complex
                         :initial-element #C(0.0d0 0.0d0))))
  (loop for k from 0 to (/ *N* 2)
        do (setf (aref H k) (exp (* j p k k))))
  (loop for k from (1+ (/ *N* 2)) to (1- *N*)
        do (setf (aref H k) (conjugate (aref H (- *N* k)))))
  H)

試験的にやっているので、まずは16サンプルだけ作ります。実行すると以下のような値が得られます。

#(#C(1.0d0 0.0d0)
  #C(0.9415440651830208d0 0.33688985339222005d0)
  #C(0.19509032201612833d0 0.9807852804032304d0)
  #C(-0.9987954562051724d0 0.049067674327417966d0)
  #C(0.7071067811865474d0 -0.7071067811865477d0)
  #C(-0.6715589548470176d0 0.7409511253549598d0)
  #C(0.9807852804032305d0 -0.19509032201612808d0)
  #C(-0.4275550934302845d0 -0.9039892931234422d0)
  #C(-1.0d0 8.572527594031472d-16)
  #C(-0.4275550934302845d0 0.9039892931234422d0)
  #C(0.9807852804032305d0 0.19509032201612808d0)
  #C(-0.6715589548470176d0 -0.7409511253549598d0)
  #C(0.7071067811865474d0 0.7071067811865477d0)
  #C(-0.9987954562051724d0 -0.049067674327417966d0)
  #C(0.19509032201612833d0 -0.9807852804032304d0)
  #C(0.9415440651830208d0 -0.33688985339222005d0))

TSP信号であれば、フラットな周波数特性を持っているはずです。以下のように調べるとちゃんと振幅1になっていることがわかりました。

(mapcar #'abs (coerce #( 上記の値 ) 'list)
;;=> (1.0d0 1.0d0 1.0d0 1.0d0 1.0d0 0.9999999999999999d0 1.0d0 1.0d0 1.0d0 1.0d0
;;    1.0d0 0.9999999999999999d0 1.0d0 1.0d0 1.0d0 1.0d0)

時間領域にもっていく

時間領域にもっていかないといけないのでIFFTをかけますが、結果の複素数から実数部だけを抜き出す必要があります。実行すると時間領域の信号が得られました。

(ql:quickload 'bordeaux-fft)

(coerce (mapcar #'realpart
                (coerce (bordeaux-fft:sifft H) 'list))
        'simple-vector)
;;=> #(0.09082711803830656d0 0.15284479737458728d0 -0.09082711803830684d0
;;     0.24317632014210552d0 -0.30019142734700194d0 0.3522809059032187d0
;;     -0.18299922203879876d0 -0.16995301543874267d0 0.37991847786317d0
;;     -0.0038532185877515057d0 -0.3799184778631696d0 -0.1699530154387425d0
;;     0.1829992220387991d0 0.35228090590321937d0 0.3001914273470015d0
;;     0.2431763201421058d0)

関数にした

関数にまとめたのが以下のものです。(make-oatsp 4096 1200)のようにすると、時間領域の測定信号が得られます。上ではやりませんでしたが、時間領域にしたときに環状シフトしてあげる必要もあるので、その処理を追加しています。（上記に続いて実行するとdefconstantで作ったjがかち合ってしまうので、Lispを再起動するか(mkunbound 'j)で変数を消しておくといいです。）

(ql:quickload 'bordeaux-fft)

(defun make-oatsp (N m)
  "Make time-domain signal of optimum aoshima's time-stretched pulse."

  ;; design the signal in frequency domain
  (setq H (let ((j #C(0.0d0 1.0d0))
                (p (/ (* m pi) (expt (/ N 2) 2)))
                (H (make-array N :element-type 'complex
                                 :initial-element #C(0.0d0 0.0d0))))
            (loop for k from 0 to (/ N 2)
                  do (setf (aref H k) (exp (* j p k k))))
            (loop for k from (1+ (/ N 2)) to (1- N)
                  do (setf (aref H k) (conjugate (aref H (- N k)))))
            H))

  ;; apply inverse fourier transform
  (setq H2 (coerce (mapcar #'realpart
                           (coerce (bordeaux-fft:sifft H) 'list))
                   'simple-vector))

  ;; circular rotation of time signal
  (let ((H3 (make-array (array-dimensions H2)
                        :element-type 'double-float))
        (n-rot (- (/ N 2) m)))
    (loop for k from 0 to (1- n-rot)
          do (setf (aref H3 (+ (- N n-rot) k)) (aref H2 k)))
    (loop for k from n-rot to (1- N)
          do (setf (aref H3 (- k n-rot)) (aref H2 k)))
    H3))

(make-oatsp 4096 1200)

速度も実用上はまったく問題ないようです。信号長221というとても長い信号（96 kHzで約22秒）でも、Intel i5 (10400)上のSBCLで1秒かからず計算できました。普段は218くらい（48 kHzで約5秒）を使いますし、一度だけ計算すればいいものなので、計算時間は気になりません。

まだCommon Lispでグラフを書いたりWAVに保存したりできていないので、この結果をREPLからテキストでコピペして（1行に1値になるようにするのと1.23d-4を1.23e-4と書き換えるのだけやりました）、その後でJuliaでプロットしました。Julia側では以下の様なことをやっています（Juliaは楽チンで高速で最高です）。

X = parse.(Float64, readlines("hoge.txt"))
snd = SampleBuf(map(PCM16Sample, X), 48000)
save("oatsp.wav", snd)
plot(X, legend=false)
savefig("oatsp.png")

その結果出てきたのが以下のグラフです。Suzukiらの論文のFig. 2(a')とそっくりで、高い周波数から低い周波数に向かってスイープダウンする信号になっているのが分かります。上記の関数と同様にして、測定関数の逆関数も作ってあげれば良さそうですね。