トップ > LISP > Common LispでFFTライブラリを使ってみる

Common LispでFFTライブラリを使ってみる

LISP

Common Lispを使っていろいろと車輪の再実装を続けています。今回はフーリエ変換ライブラリを試してみます。Common Lispにはフーリエ変換のライブラリは入っていないので、Quicklispでfftをキーワードに探していきます。すると以下の4つがヒットしました。

CL-USER> (ql:system-apropos "fft")
#<SYSTEM bordeaux-fft / bordeaux-fft-20150608-http / quicklisp 2019-12-27>
#<SYSTEM fft / fft-20180711-git / quicklisp 2019-12-27>
#<SYSTEM napa-fft3 / napa-fft3-20151218-git / quicklisp 2019-12-27>
#<SYSTEM pfft / fft-20180711-git / quicklisp 2019-12-27>

これらに関係しそうなのは以下のページのようです。fftとpfftは同じライブラリに含まれるもののようです。他にもf2clというプロジェクトにもFFTは含まれているようでした。

どれが良いのかわからないので、なんとなく最初に出てきたBordeaux-FFTを使ってみました。（Napa-FFT3も高機能そうで良いかもです）

インストールと読み込みはQuicklispを使えば簡単でした。（Windows 10とSBCLを使っています）

(ql:quickload 'bordeaux-fft)

長さ16の正弦波（周波数は3 Hzとします）を作って、それにFFTをかけました。sfftというのはsimple fftの略あるいはstupid fftの略らしく、なにも考えずに配列のFFTを計算してくれるものです。fftは配列要素の型を指定したりする必要があるのですが、今回はシンプルなほうを使います。

;; 長さを16とする（2のべき乗にする）
(defparameter *N* 16)

;; 正弦波を入れる配列を作る
(setq sinusoid (make-array *N*))

;; 3 Hzの正弦波を要素ごとに計算
(loop for ii from 0 to (1- *N*)
      do (setf (aref sinusoid ii)
               (sin (* 3 2 pi (/ ii *N*)))))

;; FFTをかけてfftresにバインド
(setq fftres (bordeaux-fft:sfft sinusoid))

この結果、fftresには以下のような16要素の複素数配列が入りました。

#(#C(-2.4322862379963016d-16 0.0d0)
  #C(6.107235716236265d-16 1.4695673325073857d-15)
  #C(2.4963710534292284d-16 4.373545779112952d-16)
  #C(-4.5862415333198d-15 -8.0d0)
  #C(1.588639366831878d-15 -7.216449660063518d-16)
  #C(7.669517010585668d-17 0.0d0)
  #C(4.851509741454889d-16 -6.728684467138615d-16)
  #C(5.6559388869686625d-16 1.8962923073275474d-15)
  #C(-1.4644739508873025d-15 0.0d0)
  #C(-1.3244021561850481d-17 -2.357745752207511d-15)
  #C(4.851509741454891d-16 6.728684467138614d-16)
  #C(7.428289848809507d-16 -1.3322676295501878d-15)
  #C(1.588639366831878d-15 7.216449660063518d-16)
  #C(-1.255572459444331d-15 8.0d0)
  #C(2.496371053429228d-16 -4.373545779112951d-16)
  #C(9.20064081065035d-16 -1.1993546792729692d-16))

最初の#C(-2.43d-16 0.0d0)はDCオフセットの値で、2番目の要素以降が各周波数ビンの値です。4番目の#C(-4.59d-15 -8.0d0)と14番目の#C(-1.26d-15 8.0d0)が3 Hzの振幅値に対応しています。絶対値が8になっているのは正の周波数と負の周波数に分割されているからで、合計すると16になります。これを信号の長さで正規化すると（各要素を信号長16で除算すると）1.0になります。

CL-USER> (abs #C(-4.5862415333198d-15 -8.0d0))
8.0d0
CL-USER> (abs #C(-1.255572459444331d-15 8.0d0))
8.0d0

配列の最初から、DC成分、正の周波数成分、ナイキスト成分、負の周波数成分という順に値が入っているので、正負の周波数成分をまとめます。このあたりは以前Julia版を書いたときに説明したので、そちらを参照してください。

marui.hatenablog.com

ここでは、まずampspecという振幅スペクトルを入れる配列を作ります。この配列の長さはN/2+1となります（今回は16/2+1で長さ9です）。次に、FFT結果（fftres）の配列の0番目（DC成分）とN/2-1番目（ナイキスト成分）をampspecにコピーします。そして、それらの間の成分について、負の周波数成分の複素共役を正の周波数成分に加算します。

(setq ampspec (make-array (1+ (/ (array-dimension fftres 0) 2))))

(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 (- *N* ii))))))

結果、ampspecに得られるのが複素数の周波数スペクトルです。

#(#C(-2.4322862379963016d-16 0.0d0)
  #C(1.5307876526886614d-15 1.5895028004346827d-15)
  #C(4.992742106858457d-16 8.747091558225903d-16)
  #C(-5.841813992764131d-15 -16.0d0)
  #C(3.177278733663756d-15 -1.4432899320127035d-15)
  #C(8.195241549868074d-16 1.3322676295501878d-15)
  #C(9.70301948290978d-16 -1.3457368934277227d-15)
  #C(5.523498671350158d-16 4.254038059535058d-15)
  #C(-1.4644739508873025d-15 0.0d0))

振幅スペクトルにするには、配列の各要素にabsを適用して計算するのですが、配列にはmapcarが使えないのでリストに変換してからmapcarをかけてあげます。（配列にmapをかけるarray-mapという関数がStack Overflowにあったので、それを使ってもいいかもです）

CL-USER> (mapcar #'abs (coerce ampspec 'list))
(2.4322862379963016d-16 2.2067691293413003d-15 1.0071697199260125d-15 16.0d0
 3.489725774217968d-15 1.5641473323680597d-15 1.6590640907420562d-15
 4.289747100668858d-15 1.4644739508873025d-15)

さらに信号長で正規化すれば、所望の振幅スペクトルが得られます。ちゃんと3 Hzに対応した場所に振幅1.0があり、それ以外の要素はゼロ近傍になっていますね。

CL-USER> (mapcar #'(lambda (x) (/ x *N*))
                 (mapcar #'abs (coerce ampspec 'list)))
(1.5201788987476885d-17 1.3792307058383127d-16 6.294810749537578d-17 1.0d0
 2.18107860888623d-16 9.775920827300373d-17 1.0369150567137852d-16
 2.6810919379180364d-16 9.15296219304564d-17)

次は、これを関数化するのと、WAVファイルから読み込んだ音に適用するのをやっていきたいです。