1. はじめに

FFT（高速フーリエ変換）の話題です。
まずは、ただコピーして使えるやつを示します。
そのあと長々と説明します。

(defun fftcore (v n sign)
  ;;  bit reverse
  (let (a b y)
    (dotimes (x n)
      (setq a (ash n -1) y 0 b x)
      (do () ((zerop a))
        (setq y (logior (ash y 1) (logand b 1)))
        (setq b (ash b -1))
        (setq a (ash a -1)))
      (when (< x y)
        (rotatef (elt v x) (elt v y)))))
  ;;  fft
  (do* ((p0 (float (* sign 2.0L0 pi) sign))
        (n2 (ash n -1))
        (a n2 (ash a -1)))
    ((< a 1) v)
    (do* ((n3 (/ n2 a))
          (n4 (ash n3 1))
          (p1 (/ p0 n4))
          (b 0 (1+ b)))
      ((<= a b))
      (do ((n5 (* n4 b))
           (c 0 (1+ c))
           x y z)
        ((<= n3 c))
        (setq x (+ n5 c))
        (setq y (+ x n3))
        (setq z (* (elt v y) (cis (* p1 c))))
        (setf (elt v y) (- (elt v x) z))
        (incf (elt v x) z)))))

2. 使い方

fftcoreは、下記の演算を一括して求めます。

$\displaystyle{V_k = \sum_{n=0}^{N-1} v_n \exp \left( i s \frac{2 \pi n k}{N} \right) }$

$N$ は配列の個数であり、2のべき乗である必要があります。
$v_n$ は変換前の配列の値であり、 $V_k$ は変換後の値です。
$s$ は符号であり、通常 $1$ か $-1$ です。

関数fftcoreの説明をします。

(defun fftcore (v n sign) ...) -> v

vは一次元の配列です。
nは値の個数であり、2のべき乗である必要があります。
たとえば256や1024などになります。
signは符号であり、通常は1.0か-1.0です。
signは符号の他にも意味があり、値の型で演算の精度を決定します。
例えばdouble-floatで演算したい場合は、1.0d0か-1.0d0を指定してください。
配列vは必ずn以上の個数を保有する必要があります。
演算結果は配列vに上書きされます。
関数fftcoreの返却値は、引数vそのものです。

3. 使ってみる

fftcoreの目的は、コピーで手軽に使えることと、必要最小限の機能を提供することです。

ところでみなさんはFFTを使ったことがありますか？
人によっては見たことすらないと思います。
自分は一時期、朝から晩までずっとFFT祭りをしてたことありました。
せっかくなのでその時のやつをCommon Lispで作ってみようと思い公開しました。

ここではフーリエ変換や、離散フーリエ変換の説明はしません。
高速フーリエ変換のアルゴリズムも説明しません。
ただし、使うための入り口だけは説明しなければなりません。

離散フーリエ変換と逆変換の一例を次に示します。

離散フーリエ変換

$\displaystyle{V_k = \frac{1}{N} \sum_{n=0}^{N-1} v_n \exp \left( - i \frac{2 \pi n k}{N} \right) }$

逆変換

$\displaystyle{v_n = \sum_{k=0}^{N-1} V_k \exp \left( i \frac{2 \pi n k}{N} \right) }$

「一例」です。
FFTに詳しい人が言ってたのをなんとなく覚えているのですが、正変換の係数 $1/N$ が逆変換についてるパターンがあったし、正変換と逆変換の $\exp$ の符号が逆になっているのも存在したとのこと。
さらに言うと、係数 $1/N$ は、 $1/\sqrt{N}$ になってたのもあったんだそうです。
たぶん対称性とかそういうのを考えた結果でしょうね。

人によって違うので、全通りあるらしいです。
だから、fftcoreでは係数の乗算をしないし、符号も決め打ちしなかったのです。
わかんないんで自分でやって下さい。
でも、せっかくなんで一例を作ってみます。

(defun fft (v n &optional (type 'double-float))
  (fftcore v n (coerce -1.0 type))
  (dotimes (x n v)
    (setf (elt v x) (/ (elt v x) n))))

(defun ifft (v n &optional (type 'double-float))
  (fftcore v n (coerce 1.0 type)))

周波数解析もやってみましょうか。
次のテスト用の信号を用意します。

$\displaystyle{A(t) = 1.23 + 2 \sin \omega t + 3 \sin 5 \omega t + 4 \sin \left( 7 \omega t + \frac{\pi}{3} \right) }$

こんな感じ

(defun a-sin-f-phi (a f x n phi)
  (* a (sin (+ (/ (* 2.0d0 pi f x) n) phi))))

(defun make-test-signal (v n)
  (let ((bias 1.23d0) a b c)
    (dotimes (x n v)
      (setq a (a-sin-f-phi 2.0d0 1.0d0 x n 0.0d0))
      (setq b (a-sin-f-phi 3.0d0 5.0d0 x n 0.0d0))
      (setq c (a-sin-f-phi 4.0d0 7.0d0 x n (/ pi 3.0d0)))
      (setf (elt v x) (complex (float (+ bias a b c) 1.0d0) 0.0d0)))))

それでは、フーリエ変換してみましょう。

(defvar *size* 1024)
(let ((v (make-array *size* :element-type '(complex double-float))))
  (make-test-signal v *size*)
  (fft v *size*)
  (dotimes (i 10)
    (let* ((z (elt v i))
           (re (realpart z))
           (im (imagpart z))
           (a (abs z)))
      (format t "~3A: #c(~6F, ~6F) -> ~6F~%" i re im a))))

配列の個数は1024個とりました。
ここは必ず2の何乗の値にしてくださいね。
2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, ...みたいな数のことです。
2の何乗じゃなくてもいいFFTも存在するらしいですが、やりたいなら自分で作って下さい。

実行結果は下記の通り。

0  : #c(  1.23,    0.0) ->   1.23
1  : #c(  -0.0,   -1.0) ->    1.0
2  : #c(  -0.0,   -0.0) ->    0.0
3  : #c(  -0.0,    0.0) ->    0.0
4  : #c(   0.0,    0.0) ->    0.0
5  : #c(  -0.0,   -1.5) ->    1.5
6  : #c(  -0.0,    0.0) ->    0.0
7  : #c(1.7321,   -1.0) ->    2.0
8  : #c(   0.0,   -0.0) ->    0.0
9  : #c(   0.0,   -0.0) ->    0.0

それぞれの要素の実部・虚部の他に、絶対値（振幅）を出力しています。
0番目は直流成分なので、1.23がそのまま入っています。
周波数成分は、512以降にも分散されているので、振幅の値は半分になっています。
1～511番目の振幅を、単純に2倍しましょう。
つまり、

1番目の振幅は2。
5番目の振幅は3。
7番目の振幅は4。

テスト用の信号の振幅があらわされているのが分かります。

4. もともとのソース

もとは、大昔に自分が作ったC言語のソースでした。
とてもひどいものだったので書き直しましたが、こんな感じです。

#define PI 3.141592653589793238462643383276

void fft_double_complex(double complex v[], size_t n, double sign)
{
    size_t x, y, a, b, c, n2, n3, n4, n5;
    double p0;
    double complex p1, z;

    /* bit reverse sort */
    for (x = 0; x < n; x++) {
        y = 0;
        b = x;
        for (a = n >> 1; a != 0; a >>= 1) {
            y = (y << 1) | (b & 1);
            b >>= 1;
        }
        if (x < y) {
            z = v[x]; v[x] = v[y]; v[y] = z;
        }
    }

    /* fft */
    p0 = sign * 2.0 * PI;
    n2 = n >> 1;
    for (a = n2; a >= 1; a >>= 1) {
        n3 = n2 / a;
        n4 = n3 << 1;
        p1 = I * p0 / n4;
        for (b = 0; b < a; b++) {
            n5 = n4 * b;
            for (c = 0; c < n3; c++) {
                x = n5 + c;
                y = x + n3;
                z = v[y] * cexp(p1 * c);
                v[y] = v[x] - z;
                v[x] += z;
            }
        }
    }
}

余談ですが、FreeBSDなら、円周率を求めるときは次のコマンドが便利です。

$ echo 'scale=30; a(1)*4' | bc -l
3.141592653589793238462643383276

Linuxならどうするんだろう。

まあそれはいいとして、さらにもとのソースが存在します。
C99なんて存在すら知らなかったときだったので、複素数の型を使っていません。

こんな感じでした。

#define PI 3.141592653589793238462643383276

void fft_double(double v[], size_t n, double sign)
{
    size_t x, y, x2, y2, x3, y3, a, b, c, n2, n3, n4, n5;
    double p0, p1, theta, w1, w2, z1, z2;

    /* bit reverse sort */
    for (x = 0; x < n; x++) {
        y = 0;
        b = x;
        for (a = n >> 1; a != 0; a >>= 1) {
            y = (y << 1) | (b & 1);
            b >>= 1;
        }
        if (x < y) {
            x2 = x << 1; x3 = x2 + 1;
            y2 = y << 1; y3 = y2 + 1;
            z1 = v[x2]; v[x2] = v[y2]; v[y2] = z1;
            z2 = v[x3]; v[x3] = v[y3]; v[y3] = z2;
        }
    }

    /* fft */
    p0 = sign * 2.0 * PI;
    n2 = n >> 1;
    for (a = n2; a >= 1; a >>= 1) {
        n3 = n2 / a;
        n4 = n3 << 1;
        p1 = p0 / n4;
        for (b = 0; b < a; b++) {
            n5 = n4 * b;
            for (c = 0; c < n3; c++) {
                x = n5 + c;
                y = x + n3;
                x2 = x << 1; x3 = x2 + 1;
                y2 = y << 1; y3 = y2 + 1;
                theta = p1 * c;
                w1 = cos(theta); w2 = sin(theta);
                z1 = v[y2]*w1 - v[y3]*w2;
                z2 = v[y2]*w2 + v[y3]*w1;
                v[y2] = v[x2] - z1;
                v[y3] = v[x3] - z2;
                v[x2] += z1;
                v[x3] += z2;
            }
        }
    }
}

こちらは、doubleの配列を2倍の長さで用意して、実部、虚部、実部、虚部、・・・みたいに交互に値を格納していく方式です。

他にもいろいろ探してたら、多次元高速フーリエ変換なんてのもありました。
ただし、この辺りになると速度的に並列実行しないと厳しいと思います。
今回のようなコピーで作れる高速フーリエ変換なんてやめて、素直に頭いい人たちが作ったライブラリを使うべきなんだと思います。

一応言っておきますが、私は数値解析の専門家じゃないんで、ここのソースは全部怪しいです。
現時点でnptでは動きませんでした。
どうもnptの複素数のexpあたりに問題がありそう（【追記】そうではなくfloat関数にバグがありました。修正済み）。
バグの発見ができてうれしいので、そのうち直します。

もし怪しくてもいいから使いたいという人がいるなら、ご自由にどうぞ。
私は使います。

nptclのブログ

Common Lisp処理系nptの開発メモです。https://github.com/nptcl/npt

Common LispでFFTを使う

1. はじめに

2. 使い方

3. 使ってみる

4. もともとのソース