スペクトログラム

Common Lisp振幅スペクトルの計算ができるようになり、Hann窓もかけられるようになったので、次はスペクトログラムの実装にいきたいと思います。

スペクトログラムは、振幅あるいはパワーの強さを時間と周波数ごとに計算したものです。以下の画像ではAudacityを使ってクラシックギターの撥弦音を上側が時間波形・下側がスペクトログラムとして表示しています。スペクトログラムは横軸が時間、縦軸が周波数、色の濃さが振幅の強さになっています。倍音成分ごとの強さといった、振幅スペクトルの時間変化を直感的に見ることができます。

スペクトログラムを計算するには、時間波形を短時間に区切ってスペクトルを計算します。以下のコードではslicesigというのを1フレームとして、その振幅スペクトルsliceampを計算して、スペクトログラムresに格納してゆきます。時間ビン、周波数ビンの値も返すようにしました。気の利いたプログラムであれば窓関数の指定やフレームの重なり量(あるいはホップサイズ)が指定できるのですが、簡単のために窓関数はHann、重なり量は50パーセントで固定してしまいました。

(defun spectrogram (x &optional (fs 48000) (winsize 1024))
  "Compute amplitude spectrogram from given sound signal <x> (float
vector) and sampling frequency <fs> with <winsize>. Each temporal
slice is processed with Hann window and 50 percent overlap is applied.
The function returns a list containing a spectrogram matrix, time
bins, and frequency bins."
  (let* ((siglen (length x))
         (hopsize (/ winsize 2))
         (sig (zeropad x (* winsize (ceiling (/ siglen winsize)))))
         (nfreq (1+ (/ winsize 2)))
         (ntime (1- (/ (length sig) hopsize)))
         (win (coerce (hann winsize) 'vector))
         (slicesig (make-array winsize))
         (sliceamp (make-array nfreq))
         (res (make-array `(,ntime ,nfreq))))
    (loop for nt from 0 to (1- ntime)
          do (loop for n from 0 to (1- winsize)
                   do (setf (aref slicesig n)
                            (* (aref win n)
                               (aref sig (+ n (* hopsize nt))))))
             (setq sliceamp (first (spectrum slicesig fs)))
             (loop for n from 0 to (1- nfreq)
                   do (setf (aref res nt n) (aref sliceamp n))))
    (list res
          (linspace 0 (/ (length sig) fs) ntime)
          (linspace 0 (/ fs 2) nfreq))))

このままだと結果がわからないので、CSVファイルとして出力してみました。

(let* ((wav (wavread "guitar.wav"))
       (x (first wav))
       (fs (second wav))
       (spec (first (spectrogram x fs))))
  (with-open-file (out #P"out.csv" :direction :output :if-exists :supersede)
    (loop for r from 0 to (1- (array-dimension spec 0))
          do (loop for c from 0 to (1- (array-dimension spec 1))
                   do (when (> c 0) (format out ","))
                      (format out "~F" (aref spec r c)))
             (format out "~%"))))

一応それっぽい値が出ています。Matlabでプロットしてみましょう。

X = csvread('out.csv');
imagesc(log(X'));
axis xy;
colormap gray;

横軸・縦軸はインデックス番号になっていますが、中身のスペクトログラムはちゃんと計算できているようです。

スペクトログラムが計算できたということは、時間フレームごとの周波数分析結果を使った処理ができるということになります。たとえばスペクトル重心の時間変化を計算するといったことができるわけです。いろいろ遊んでいこうと思います。