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

TSP信号の生成

インパルス応答

頭部伝達関数を測定したり、スピーカーの応答を測ったり、室内音響指標の計算に使ったり……インパルス応答が測定できると色々と便利です。すでに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)²
• H(k) = exp(j p k²), 0<=k<=N/2
• H(k) = conj(H(N-k)), N/2+1 <= k < N

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

• H(k) = exp(-j p k²), 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)

速度も実用上はまったく問題ないようです。信号長2²¹というとても長い信号（96 kHzで約22秒）でも、Intel i5 (10400)上のSBCLで1秒かからず計算できました。普段は2¹⁸くらい（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')とそっくりで、高い周波数から低い周波数に向かってスイープダウンする信号になっているのが分かります。上記の関数と同様にして、測定関数の逆関数も作ってあげれば良さそうですね。

Common Lispで音の振幅スペクトルの計算

先月末「Common LispでFFTライブラリを使ってみる - 丸井綜研」を書きました。このエントリの最後に「関数化した上でWAVファイルを読み込んで分析してみたい」と書いていたので、今日はそれをやってみました。

WAVファイルを読み込むのには「Common Lispで音ファイルを読み込む - 丸井綜研」でも使ったcl-wavを使います。以下のようにすれば、音を読み込んだあとでWAVに関する情報が得られます。ここではWindowsに最初から入っている44.1 kHz / 16 bit / 2 chの音ファイルを使ってみました。

(ql:quickload :cl-wav)

(defun wavread (filename)
  (wav:read-wav-file
   filename
   :chunk-data-reader
   (wav:wrap-data-chunk-data-samples-reader)))

(setq snd (wavread "c:/windows/media/ding.wav")) ; ファイル名を指定して読み込み
(setq nch (fourth (sixth (second snd))))         ; チャンネル数
(setq fs (sixth (sixth (second snd))))           ; 標本化周波数
(setq nbit (nth 11 (sixth (second snd))))        ; ビット数
(setq data (sixth (last snd)))                   ; 音サンプルのデータ

ただし、音データはWAVファイルに入っていたデータがそのまま入っているようです。今回は2チャンネルのデータなので、LRLR...という順序でデータが並んでいます。なので、たとえば以下のように2つに分割しなければなりません。

(setq dataL (make-array (/ (length data) 2)))
(setq dataR (make-array (/ (length data) 2)))
(loop for ii from 0 to (1- (length dataL))
      do (setf (aref dataL ii) (aref data (* 2 ii)))
         (setf (aref dataR ii) (aref data (1+ (* 2 ii)))))

音をFFTするのですが、長さが2のべき乗になっていないといけないので「Common Lispで畳み込み演算 - 丸井綜研」のときに作った、2のべき乗の長さにゼロ詰めする関数を使います。（以下に再掲）

(defun nextpow2 (n)
  "Find the exponent of the next smallest power of two larger than <n>."
  (ceiling (log n 2)))

(defun zeropad-pow2 (vec)
  "Return the array of length of 2^n with the contents of <vec> in the
   earlier part."
  (let ((xx (make-array (expt 2 (nextpow2 (length vec))))))
    (loop for k from 0 to (1- (length vec))
          do (setf (aref xx k) (aref vec k)))
    xx))

Matlabのlinspaceに相当する関数も用意します。これは指定した2つの数x1とx2を含む等差数列を作る関数で、オプションで個数を指定することもできます（指定しなければデフォルトでは100要素の数列ができます）。このlinspaceはDC～ナイキストの周波数ビンの中心周波数を計算するのに使います。

(defun linspace (x1 x2 &optional (n 100))
  "Generates a list of 100 equally spaced numbers between x1 and x2."
  (let ((x (loop for x from 0 to (1- n)
                 collect (/ x (1- n)))))
    (mapcar #'(lambda (v) (+ v x1))
            (mapcar #'(lambda (v) (* v (- x2 x1)))
                    x))))

(linspace 0 4 5)
;;=> (0 1 2 3 4)

(linspace 0 4 6)
;;=> (0 4/5 8/5 12/5 16/5 4)

(linspace 0.0 4.0 6)
;;=> (0.0 0.8 1.6 2.4 3.2 4.0)

「Common LispでFFTライブラリを使ってみる - 丸井綜研」で書いた手順を一つの関数にまとめたら下のようになりました。音データと標本化周波数を与えると、振幅スペクトルとそれに対応する周波数ビンの中心周波数が多値で戻ってきます。

(defun spectrum (x fs)
  "Compute amplitude spectrum from given sound signal <x> (float
   array) and sampling frequency <fs>. The funtion returns amplitude
   spectrum and corresponding frequency bin values."
  (let* ((fftres (bordeaux-fft:sfft (zeropad-pow2 x)))
         (ampspec (make-array (1+ (/ (array-dimension fftres 0) 2))))
         (freq (coerce (linspace 0.0 (/ fs 2.0) (length ampspec)) 'vector)))
    (setf (aref ampspec 0) (aref fftres 0))
    (setf (aref ampspec (1- (array-dimension ampspec 0)))
          (aref fftres (1- (array-dimension ampspec 0))))
    (loop for ii from 1 to (- (array-dimension ampspec 0) 2)
          do (setf (aref ampspec ii)
                   (+ (aref fftres ii)
                      (conjugate (aref fftres (- (length fftres) ii))))))
    (loop for ii from 0 to (1- (array-dimension ampspec 0))
          do (setf (aref ampspec ii)
                   (abs (aref ampspec ii))))
    (values ampspec freq)))

出力結果をグラフにすれば上図（これはJuliaで作りました）のように振幅スペクトルの可視化ができますが、まだCommon Lispでグラフ描画をする方法を考えていないので、ここでは使用例として音響特徴量の計算をしてみましょう。

たとえば、スペクトル重心は振幅スペクトルの平均周波数で、「音の明るさ」の感覚との対応が良いことが知られています。これを計算する関数を、定義どおりに以下のように作ってみました。スペクトルデータ（ampspec）の先頭はDC成分、最後はナイキスト成分が入っているので、それらを取り除くためにcdrとbutlastを使っています（べつに取り除かなくても問題ありません）。また、スペクトルデータも周波数もsimple-vector型で得られるので、リストにcoerceしています。

(defun spectral-centroid (data fs)
  (multiple-value-bind (ampspec freq) (spectrum data fs)
    (let ((as (cdr (butlast (coerce ampspec 'list))))
          (ff (cdr (butlast (coerce freq 'list)))))
      (/ (apply #'+ (mapcar #'* ff as))
         (apply #'+ as)))))

(spectral-centroid dataL fs)
;;=> 10158.373793577552d0

(spectral-centroid dataR fs)
;;=> 10241.781910677042d0

Matlab、R、Juliaなどでも同じ計算をするプログラムを作ってありますが、ほぼ同じ値になりました。値の違いは0.05パーセント程度でしたので、演算のどこかで誤差が入ったものと思われます。

少し使ってみて思いましたが、spectrumは多値で返してmultiple-value-bindするより、2つのベクトルが入ったリストで返すほうが使い勝手がいいかもしれませんね。（spectrumの最後のvaluesをlistにするだけの変更です）

Common Lispで畳み込み演算

先日Common LispでFFTが使えるようになったので、それを使って畳み込みをします。畳み込みには環状畳み込み（円状畳み込みや循環畳み込みとも呼びます）と線形畳み込みの2種類がありますので、両方とも作ってみようと思います。

ちょっとした道具

実際の畳み込みのプログラムに入る前に、Matlabを参考にしてnextpow2を作ります。これは与えられた値より大きい最小の2のべき乗の指数を計算するものです。たとえば5の次の2のべき乗は8=2³なので3、12の次の2のべき乗は16=2⁴なので4が返ります。（ceilingは多値を返すので実際には余りも出てきますが、ここでは1つ目の返り値だけを使います。）

(defun nextpow2 (n)
  "Find the exponent of the next smallest power of two larger than n."
  (ceiling (log n 2)))

(nextpow2 5)                            ;=> 3
(nextpow2 12)                           ;=> 4

これを使うと、たとえば任意のベクトルを渡すと、2のべき乗の長さの先頭部分にそのベクトル内容を入れて残りはゼロ詰めしたものを返すような関数が作れます。

(defun zeropad-pow2 (vec)
  "Return the array of length of 2^n with the contents of <vec> in the
   earlier part."
  (let ((xx (make-array (expt 2 (nextpow2 (length vec))))))
    (loop for k from 0 to (1- (length vec))
          do (setf (aref xx k) (aref vec k)))
    xx))

(zeropad-pow2 #(1 2 3))                 ;=> #(1 2 3 0)
(zeropad-pow2 #(1 2 3 4 5))             ;=> #(1 2 3 4 5 0 0 0)

環状畳み込み（circular convolution）

さて、ここからが本題です。ある程度の長さのあるベクトルどうしを畳み込みする際には、高速フーリエ変換を経由すると高速に実行できます。すなわち、vec1とvec2の畳み込みをするときに、それぞれをフーリエ変換してy1とy2を計算してから、y1とy2の要素ごとの積をとり、その結果を逆フーリエ変換するという手順です。なおFFTの制約もあったりして、vec1とvec2は同じ長さで、かつ2のべき乗の長さになっている必要があります。

(defun fftconv-circular (vec1 vec2)
  "Circular convolution using FFT. The two vectors must be of equal
   length as well as 2^n in length."
  (coerce
   (mapcar #'realpart
           (coerce (bordeaux-fft:sifft
                    (let ((y1 (bordeaux-fft:sfft vec1))
                          (y2 (bordeaux-fft:sfft vec2)))
                      (loop for e1 across y1 and e2 across y2
                            collect (* e1 e2))))
                   'list)) 'simple-vector))

(fftconv-circular #(1 2 3 4) #(1 1 0 0)) ;=> #(5.0d0 3.0d0 5.0d0 7.0d0)

上記のコードでは、逆フーリエ変換後に実部だけを取ってくるために一旦ベクトルをリストにしてmapcarをかけています。なんとなく効率よくない気がするので、ベクトルのままループで実部を取るものも作ってみました。

(defun fftconv-circular (vec1 vec2)
  "Circular convolution using FFT. The two vectors must be of equal
   length as well as 2^n in length."
  (let ((z (bordeaux-fft:sifft
            (let ((y1 (bordeaux-fft:sfft vec1))
                  (y2 (bordeaux-fft:sfft vec2)))
              (loop for e1 across y1 and e2 across y2
                    collect (* e1 e2)))))
        (res (make-array  (length vec1))))
    (loop for k from 0 to (1- (length vec1))
          do (setf (aref res k) (realpart (aref z k))))
    res))

どちらのほうが実行速度が速いのか、以下のように乱数列どうしの畳み込みで調べてみました。Intel i5 10400で実行した結果、どちらも14～15秒くらいでした。リストに変換してmapcarするよりも、ループを使って実部をとるほうが3～4%ほど速いことが分かりましたので、今後はそのやりかたで行こうと思います。

（比較のためにMatlab 2022aのSignal Processing Toolboxに入っているcconvを使ってtic; for n=1:100; cconv(rand(2^18,1), rand(2^18,1)); end; tocとしたところ10秒くらいでした。Julia 1.8.1でもDSP.jl 0.7.7を使って@time for n=1:100; conv(rand(2^18), rand(2^18)); endをしたら2.5秒でした。やはりJulia速い……）

(defun make-random-array (N)
  (let ((res (make-array N)))
    (loop for k from 0 to (1- (array-total-size res))
          do (setf (row-major-aref res k) (random 1.0d0)))
    res))

(time
 (dotimes (k 100)
   (let ((x (fftconv-circular
             (make-random-array (expt 2 18))
             (make-random-array (expt 2 18))))))))

線形畳み込み（linear convolution）

線形畳み込みでは、vec1の長さ + vec2の長さ - 1の長さの結果が得られます。たとえば[1 2 3 4]と[1 1]の畳み込み結果は[1 3 5 7 4]になります。bordeaux-fftでは2のべき乗の長さのベクトルしかFFTがかけられないのと、フーリエ領域で要素ごとの積を計算するときには2つのベクトルの長さが一致していなければいけないことから、vec1の長さ + vec2の長さ - 1を超える最小の2のべき乗のベクトルを作り、ゼロ詰めを行います。[1 2 3 4]と[1 1]の例では畳み込み後の長さが4+2-1=5になるので、xx1とxx2はそれぞれ[1 2 3 4 0 0 0 0]と[1 1 0 0 0 0 0 0]というベクトルになります。これらを環状畳み込みに渡してあげると[1 3 5 7 4 0 0 0]という結果がもどってくるので、先頭の5要素を抜き出してあげるというわけです。

(defun fftconv-linear (vec1 vec2)
  "Linear convolution using FFT."
  (let ((xx1 (make-array (expt 2 (nextpow2 (+ (length vec1) (length vec2) -1)))))
        (xx2 (make-array (expt 2 (nextpow2 (+ (length vec1) (length vec2) -1)))))
        (res (make-array (+ (length vec1) (length vec2) -1))))
    (loop for k from 0 to (1- (length vec1))
          do (setf (aref xx1 k) (aref vec1 k)))
    (loop for k from 0 to (1- (length vec2))
          do (setf (aref xx2 k) (aref vec2 k)))
    (setq foo (fftconv-circular xx1 xx2))
    (loop for k from 0 to (1- (length res))
          do (setf (aref res k) (aref foo k)))
    res))

(fftconv-linear #(1 2 3 4) #(1 1)) ;=> #(1.0d0 3.0d0 5.0d0 7.0d0 4.0d0)

Common Lispで音ファイルを読み込む

（2022-12-04追記：cl-wavを便利に使える関数を作りました→Common Lisp用のwavread - 丸井綜研）

Lispで音をやろうと思ってもなかなか情報が出てこないのですが、OpenMusicのような大規模なシステムもCommon Lispで書かれているので無理なことはないはず。Google検索してみると、masatoiさんが音声合成関係の記事をいくつか書かれていたり、t-sinさんが発表スライドを公開されていたりします。

Sounds Like Common Lisp - ゼロからはじめるサウンドプログラミング from t-sin

www.slideshare.net

このスライドの中で紹介されていたpukunuiというシンセサイザーのソースコードを読んでみたところ、cl-wavとcl-portaudioというライブラリがあるのを知ることができました。音ファイルの読み書きと録音再生ができれば、あとはやりたいこと次第でなんとでもなりそうです。

今日は音ファイルの読み込みをやってみます。

cl-wavを試す

WAVファイルの読み込みをするために、cl-wavを試してみます。エレキベースのライン出力を録音した30秒くらいのファイルがあるので、それを読み込みます。

(ql:quickload '(cl-wav)

(setq snd (wav:read-wav-file "EBass_E1.wav"))

このsndというリストには((riff) (fmt) (data))という3つのリストが入っていて、RIFFファイルのチャンク構造がそのままになっています。そのため、first、second、thirdを使って表示してあげると以下のようになります。特に2番目のfmtチャンクに音のメタデータが入っていて、今回のファイルはチャンネル数（NUMBER-OF-CHANNELS）が1、標本化周波数（SAMPLE-RATE）が48000 Hz、ビット数（SIGNIFICANT-BITS-PER-SAMPLE）が16だと読み取れます。

(setf *print-lines* 10) ; limit the number of lines to be displayed in REPL

(first snd)                             ; RIFF chunk
;;=> (:CHUNK-ID "RIFF" :CHUNK-DATA-SIZE 2880036 :FILE-TYPE "WAVE")

(second snd)                            ; fnt chunk (where all information are)
;;=> (:CHUNK-ID "fmt " :CHUNK-DATA-SIZE 16 :CHUNK-DATA
;;    (:COMPRESSION-CODE 1 :NUMBER-OF-CHANNELS 1 :SAMPLE-RATE 48000
;;     :AVERAGE-BYTES-PER-SECOND 96000 :BLOCK-ALIGN 2 :SIGNIFICANT-BITS-PER-SAMPLE
;;     16))

(third snd)                             ; data chunk
;;=> (:CHUNK-ID "data" :CHUNK-DATA-SIZE 2880000 :CHUNK-DATA
;;    #(0 0 0 0 0 0 1 0 255 255 1 0 0 0 255 255 3 0 253 255 5 0 251 255 5 0 253 255
;;      3 0 0 0 0 0 1 0 2 0 254 255 5 0 253 255 3 0 0 0 1 0 2 0 1 0 0 0 1 0 2 0 0 0
;;      3 0 0 0 1 0 3 0 255 255 3 0 1 0 1 0 2 0 2 0 1 0 2 0 0 0 4 0 0 0 3 0 1 0 1 0 ...
;; ずらずらとバイト列が出てきます

dataチャンクはバイト単位の列のようです。fmtチャンクの情報を使って、自力で16ビット単位や24ビット単位にまとめて行くしかないのでしょうか。WAVはリトルエンディアンだという知識も必要です。さいわい、cl-wavのREADME.mdには例として以下のような読み込み方法も書かれていました。chunk-data-readerに読み込み方を指定すれば、-1～+1の範囲に正規化してくれるようです。

CL-USER> (wav:read-wav-file "c:/windows/media/ding.wav" :chunk-data-reader (wav:wrap-data-chunk-data-samples-reader))
((:CHUNK-ID "RIFF" :CHUNK-DATA-SIZE 70052 :FILE-TYPE "WAVE")
 (:CHUNK-ID "fmt " :CHUNK-DATA-SIZE 16 :CHUNK-DATA
  (:COMPRESSION-CODE 1 :NUMBER-OF-CHANNELS 2 :SAMPLE-RATE 44100
   :AVERAGE-BYTES-PER-SECOND 176400 :BLOCK-ALIGN 4 :SIGNIFICANT-BITS-PER-SAMPLE
   16))
 (:CHUNK-ID "data" :CHUNK-DATA-SIZE 70016 :CHUNK-DATA
  #(-3.0517578e-5 -3.0517578e-5 -3.0517578e-5 -6.1035156e-5 -3.0517578e-5
    -6.1035156e-5 -6.1035156e-5 -9.1552734e-5 -6.1035156e-5 -9.1552734e-5
    -6.1035156e-5 -9.1552734e-5 -6.1035156e-5 -9.1552734e-5 -6.1035156e-5
    -6.1035156e-5 -6.1035156e-5 -6.1035156e-5 -3.0517578e-5 -3.0517578e-5 ...

さて、マニュアルを読むと、Limitationのところに「WAVを読み込むことはできるけど、書き出しは未対応です」なんて書いてありました。Pure Common Lispっぽいので可搬性が高くて良いのですが、FLACやOGGを使いたいときもあるのでどうしたものか……。

bodge-sndfileを試す

LibSndFileは音ファイルを扱うための定番のCライブラリで、WAVだけでなく様々なファイルタイプに対応しています。それのCommon Lispラッパーがbodge-sndfileです。（他にもcl-libsndfileというのもあるようです）

(ql:quickload '(sndfile-blob bodge-sndfile))

(setq snd
      (sf:with-open-sound-file (sf-file "EBass_E1.wav")
        (sf:read-short-samples-into-array sf-file)))
;;=> #(0 0 0 1 -1 1 0 -1 3 -3 5 -5 5 -3 3 0 0 1 2 -2 5 -3 3 0 1 2 1 0 1 2 0 3 0 1 3
;;     -1 3 1 1 2 2 1 2 0 4 0 3 1 1 3 1 2 2 1 4 -1 5 -1 5 0 3 1 3 1 3 1 4 0 5 -1 5 1
;;     2 4 0 4 2 1 5 -1 6 1 2 4 0 5 1 3 3 2 2 3 3 1 6 -2 7 -1 5 2 2 4 1 5 2 2 5 1 3
;;     4 1 5 1 4 3 2 5 1 5 2 3 4 1 7 -1 6 3 1 6 2 3 4 4 1 6 2 5 2 5 3 2 6 0 7 2 6 1 ...

こんなかんじで読み込むことができます。リストではなくsimple-vectorにデータだけが入ってくるみたいですね。16ビットのファイルなら-32768（2^{16-1}）～+32767（2^{16-1}-1）の範囲のデータになっているみたいです。

;; find min and max
(loop for x across snd3 minimize x)     ;=> -21134
(loop for x across snd3 maximize x)     ;=>  23196

-1～+1の範囲への正規化はこんな感じでやってあげれば良さそう。

;; normalize according to the number of bits
(coerce
 (let ((mx (expt 2 (1- 16))))       ; change here for other frame size
   (loop for x across snd3
         collect (float (/ x mx))))
 'simple-vector)
;;=> #(0.0 0.0 0.0 3.0517578e-5 -3.0517578e-5 3.0517578e-5 0.0 -3.0517578e-5
;;     9.1552734e-5 -9.1552734e-5 1.5258789e-4 -1.5258789e-4 1.5258789e-4
;;     -9.1552734e-5 9.1552734e-5 0.0 0.0 3.0517578e-5 6.1035156e-5 -6.1035156e-5
;;     1.5258789e-4 -9.1552734e-5 9.1552734e-5 0.0 3.0517578e-5 6.1035156e-5 ...

ファイルの情報を得るには以下のようにすればいいのかな？（ファイルのフォーマット情報、チャンネル数、標本化周波数、サンプル数、です。）

(sf:with-open-sound-file (sf-file "EBass_E1.wav")
  (format t "         Format: ~A~%" (sf:sound-format sf-file))
  (format t "Num of Channels: ~A~%" (sf:sound-channels sf-file))
  (format t "    Sample Rate: ~A~%" (sf:sound-sample-rate sf-file))
  (format t "   Sound Frames: ~A~%" (sf:sound-frames sf-file)))
;;=>          Format: (WAV PCM-16 FILE)
     Num of Channels: 1
         Sample Rate: 48000
        Sound Frames: 1440000

bodge-sndfileはLibSndFileに依存していますが、quickloadの時点で必要なバイナリを一緒に含めてくれる仕組みになっているような雰囲気がします。Windows上のSBCLでは問題なくquickloadできましたが、M1 Macではquickloadの時点で「サポートされていません」と出力されました。今のところx86_64のみに対応しているのが原因のようです。

cl-wavとbodge-sndfileの2つを使ってみましたが、どちらも一長一短です。bodge-sndfileのほうが便利な機能が色々ありそうですが、M1 Macならcl-wav一択かも……。

WAVファイルの書き出しまでやっているPure Common Lispコードとしては以下のページを見つけました。WAVは複雑なファイルフォーマットではないので、これを参考にして自作するのもいいですね。

note.com

エンジンの形式

ボクサーとかVツインとかいろいろあるので、エンジンの形式についてWikipediaへのリンクをまとめてみました。エンジン形式の説明についてはJIS B 0108-1「往復動内燃機関−用語−第1部：機関設計及び運転用語」に掲載されている説明を引用しました。

直列エンジン

一列のシリンダバンクをもつ機関。

U型エンジン

2列の平行なシリンダバンクと2本のクランク軸をもつ機関。（ツインバンク機関）

V型エンジン

2列のシリンダバンクが互いにある角度をもち、1本のクランク軸をもつ機関。

水平対向エンジン（ボクサーエンジン）

2列のシリンダバンクを、クランク軸の両側に対向して配置した機関。

W型エンジン

3列以上のシリンダバンクが互いにある角度をなして1本のクランク軸をもち、両端のバンク間の角度が180度以下の機関。 （扇型機関）

X型エンジン

二つの平面上に四つのシリンダバンクが互いにある角度をなし、各平面上の二つのバンクはクランク軸の反対に対向して配置した1本のクランク軸をもつ機関。

H型エンジン

二つの平面上に四つのシリンダバンクと二本のクランクシャプトをもつ機関。各平面上の二つのバンクは、クランク軸の反対側に対向して配置してある。

星型エンジン

各シリンダローに二個以上のシリンダをもち、それがクランク軸の周囲に等間隔に配置された機関。クランク軸に二つだけのシリンダローがある場合はシリンダを互い違いに配置でき、この場合は星形機関（star engine）と呼ぶ。

ロータリーエンジン

ロータ、ハウジングなどの構成部品が各々軸心の周りに円運動を行う間欠燃焼の容積式機関。

僕が乗っているCRF1100LアフリカツインのエンジンSD08Eは「水冷4ストロークOHC（ユニカム）4バルブ直列2気筒」（総排気量 1082 cm³、内径×行程 92.0×81.4 mm、圧縮比 10.1、最高出力 75 kW、最大トルク 105 N・m）です。

• 水冷エンジン
• 4ストローク
• OHC（ユニカム）
• 4バルブ
• 直列2気筒（270度位相クランク）
• 排気量と行程容積
• 圧縮比
• 出力（仕事率と馬力）
• トルク

一方で、もう一台のSuper Cub 110ProのエンジンJA07Eは「空冷4ストロークOHC単気筒」（総排気量 109 cm³、内径×行程 50.0×55.6 mm、圧縮比 9.0、最高出力 6.0 kW、最大トルク 8.4 N・m）です。

エンジンだけなのに、こんなにも色々な用語が出てくるので、自動車好きじゃないとチンプンカンプンでしょうねぇ。

www.honda.co.jp

本日の給油（スーパーカブ110プロ／JA07）

給油日	オドメーター (km)	給油量 (L)	単価 (円/L)	燃費 (km/L)	距離単価 (円/km)
2021-09-25	12563.7	2.81	148.04	52.10	2.84
2021-11-04	12711.0	2.44	156.97	60.37	2.60
2022-01-03	12814.3	2.18	153.21	47.39	3.23
2022-02-23	12910.1	2.55	158.04	37.57	4.21
2022-03-06	13058.0	2.83	160.07	52.26	3.06
2022-05-26	13182.1	2.46	158.10	50.45	3.09
2022-06-22	13278.4	1.99	165.83	48.39	3.43
2022-07-27	13405.0	2.21	157.92	57.29	2.76
2022-08-31	13546.7	2.92	156.16	48.53	3.22

月に一度くらいのペースで給油をしてますね。今回は燃費が悪かったですが、通勤や買い物などの近距離しか走っていないせいでしょう。まぁ、そんな感じです。