先月末「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にするだけの変更です）
