高速フーリエ変換の実装を難しそうかなと思っている方が、なんだ簡単じゃないですか！！となるための実装講座です

対象読者さんはどのような方ですか？

FFT（高速フーリエ変換）の定義を知っているものの、その実装が難しそうだと感じて困っている方々です。逆に原理や有用性、理論的な子細にご興味のある方のご期待には応えられないと思います。

目標

FFT に苦手意識のあった方が、最低限動くコードを書くだけなら簡単かも？と感じてくださるまでになれたら、私はとっても嬉しいです。

離散フーリエ変換とは

定義はウィキペディアにあります。（責任放棄）

wikipedia: 離散フーリエ変換

今回採用する定義

最速で実装までたどり着きたいですから、理論的なところはスキップです。

$N = 2 ^ n$ としましょう。$N$ 次多項式を入れると $N$ 次多項式を返してくれる何かがフーリエ変換です。多項式と言いましたが、コンピュータープログラムですから、係数を並べたものだと思ってくださると嬉しいです。

複素係数 $N$ 次多項式 $f$ に対して、そのフーリエ変換 $\mathcal F (f) = \widehat f$ は次のようになります。

\widehat f (q) = \sum _ { i = 0 } ^ { N - 1 } f( \zeta ^ i ) q ^ i

いきなり数式でつらいですね。しかし大丈夫です。 $\widehat f$ は関数なのですが、関数だとではなく係数の並びだと思いましょう。定数項は $f ( 1 )$ です。次は $f ( \zeta )$、その次が $f( \zeta ^ 2 )$ です。

そうです、私たちの知りたいのは多項式ではなく、$f$ に $\zeta$ の冪を代入したときの値の一覧表です。

素直な方法

for-loop を知っていますか？私は知っています。天才ですね。もっと褒めてくださっても良いんですよ？

ところでこの章は、定義を確認するためのものですから、先程の章で十分に理解できた方は読み飛ばしてくださっても構いません。

さて、疑似コードです。

配列 $A$ を入力として受け取ります。
同じ長さの配列 $B$ を作って 0 で初期化です。
$i$ を $0$ から $N$ まで繰り返しです。（ $B _ i$ を計算しましょう。）
- $j$ を $0$ から $N$ まで繰り返しです。（ $j$ 次の項を足したいです。）
  - $B _ i$ に $A _ j \zeta ^ {i j}$ を足します。
$B$ の出来上がりです。

C# で書くと、次のようになります。

static Complex[] FourierTransformation(Complex[] a) {
    int N = a.Length;
    Complex[] b = new Complex[N];
    double circle = 2 * Math.PI;
    for(int i = 0; i &lt; N; i++) {
        for(int j = 0; j &lt; N; j++) {
            b[i] += a[j] * Complex.Exp(new Complex(0, i * j * circle / N));
        }
    }
    return b;
}

フーリエ変換デバックのコツ

フーリエ逆変換を知っていますか？私は知っています（満面の笑み）

$\zeta$ の代わりに逆数 $\zeta ^ {-1}$ を使いましょう。この変換をするとあら不思議、定数倍を除いて元に戻ってしまいます。（ちなみになのですが、標本点の数 $ N $ が多項式の次数よりも少ないと、うまく行きませんから、注意です。）

もう少し正確なお話をすると、逆変換は次のようになっています。

\frac 1 N \sum _ { i = 0 } ^ { N - 1 } f( \zeta ^ { - i } ) q ^ i

普通の変換と見比べてみましょう。$\zeta$ がその逆数にすり替わっているところと、正規化係数 $1 / N$ が掛っているところを除けば、同じです。

逆変換は、先程の実装と殆ど同じですから、同じ関数を使いましょう。第二引数に、逆変換かどうかを入れましょう。次のように変更です。するとコードは次のようになります。

static Complex[] FourierTransformation(Complex[] a, bool reverse) {
    int N = a.Length;
    Complex[] b = new Complex[N];
    double circle = (reverse ? -1 : +1) * 2 * Math.PI;
    for(int i = 0; i &lt; N; i++) {
        for(int j = 0; j &lt; N; j++) {
            b[i] += a[j] * Complex.Exp(new Complex(0, i * j * circle / N));
        }
    }
    return b;
}

これで逆変換の完成です。

これがあれば適当な配列を作って、フーリエ変換をして逆変換をして、戻るかどうかを確認をすると、実装に誤りがるかどうかがすぐにわかります。

実装の前に、数学のお時間です

$f$ の係数列 $( a _ i ) _ { i = 0 } ^ { N - 1 }$ を用いて、$f$ の $\zeta ^ j$ における値を求めましょう。これこそが知りたいことでした。

f( \zeta ^ j ) = \sum _ { i = 0 } ^ { N - 1 } a _ i \zeta ^ { i j }

突然のお知らせなのですが、添え字の偶奇で分解です。するとあら不思議、奇数パートが $\zeta ^ i$ で割れますね。

f( \zeta ^ j ) = \sum _ { i = 0 } ^ { N / 2 - 1 } a _ { 2 i } \zeta ^ { 2 i j } + \sum _ { i = 0 } ^ { N / 2 - 1 } a _ { 2 i + 1 } \zeta ^ { ( 2 i + 1) j } \\ = \sum _ { i = 0 } ^ { N / 2 - 1 } a _ { 2 i } \zeta ^ { 2 i j } + \zeta ^ j \sum _ { i = 0 } ^ { N / 2 - 1 } a _ { 2 i + 1 } \zeta ^ { 2 i j }

こちらの数式をよく眺めてみます。前半も後半もフーリエ変換ではありませんか！！嬉しいですね。

というわけで、偶数次の係数ばかりを束ねた多項式を $g$、奇数次の係数ばかりを束ねた多項式を $h$ と書きましょう。ただし単項式パートは、長いですから、畳んでおきましょう。収納上手です。形式的には次です。

g( q ) = \sum _ { i = 0 } ^ { N / 2 - 1 } a _ { 2 i } q ^ i, \quad h( q ) = \sum _ { i = 0 } ^ { N / 2 - 1 } a _ { 2 i + 1 } q ^ i

せっかく定義したのですから、使ってさしあげましょう。先ほど書いた $ f ( \zeta ^ i ) $ の計算式に、代入です。

f ( \zeta ^ i ) = g ( \zeta ^ { 2 i } ) + \zeta ^ i \cdot h( \zeta ^ { 2 i } )

実装のお時間です。

先ほどの式を素直に実装すればよいです。

疑似コードです。

配列 $A$ を入力として受け取ります。
長さが $1$ ならばこのままお返しすればよいですね。
同じ長さの配列 $FFT(A)$ を用意です。
$A$ の偶数番目だけを並べた配列を $B$、奇数番目だけを並べた配列を $C$ としましょう。
再帰呼び出し発動です！ $B, C$ のフーリエ変換 $FFT(B), FFT(C)$ を計算していただきましょう。
マージ開始です。$i$ を $0$ から $N$ まで繰り返しです。（ $FFT(A)_i$ を計算しましょう。）
- 代入です。 $FFT(A) _ i \leftarrow FFT(B) _ i + \zeta ^ i FFT(C) _ i $ です。（あらあら、添字が範囲を超えてしまします。）

さて、この通りに実装をしましょう。疑似コードでは $A$ と $FFT(A)$ のように 2 つの配列を用意してしまいましたが、兼用で良いです。それと注意なのですが、$B$ や $C$ は長さが半分です。そのままイテレートをしたらはみ出してしまいますね。しかし大丈夫です、思い出してください。フーリエ変換というものは、$1$ の冪根の累乗を代入して並べたものです。はみ出してしまっても、そのまま初めに戻ってくれば良いです。

こちらが C# によるソースコードです。たったの 16 行、嬉しいですね。

static void FFT(ref Complex[] a) {
    int N = a.Length;
    if (N == 1) return;
    Complex[] b = new Complex[N / 2];
    Complex[] c = new Complex[N / 2];
    for (int i = 0; i &lt; N; i++) {
        if (i % 2 == 0) b[i / 2] = a[i];
        if (i % 2 == 1) c[i / 2] = a[i];
    }
    FFT(b);
    FFT(c);
    double circle = 2 * Math.PI;
    for (int i = 0; i &lt; N; i++) {
        a[i] = b[i % (N / 2)] + c[i % (N / 2)] * Complex.Exp(new Complex(0, circle * i / N));
    }
}

デバックだいすきクラブのみなさんのために、reverse 付きバージョンも乗せておきますね。

static void FFT(ref Complex[] a, bool reverse){
    int N = a.Length;
    if (N == 1) return;
    Complex[] b = new Complex[N / 2];
    Complex[] c = new Complex[N / 2];
    for (int i = 0; i &lt; N; i++) {
        if (i % 2 == 0) b[i / 2] = a[i];
        if (i % 2 == 1) c[i / 2] = a[i];
    }
    FFT(b, reverse);
    FFT(c, reverse);
    double circle = (reverse ? -1 : +1) * 2 * Math.PI;
    for (int i = 0; i &lt; N; i++) {
        a[i] = b[i % (N / 2)] + c[i % (N / 2)] * Complex.Exp(new Complex(0, circle * i / N));
    }
}

この記事でお話できなかったこと

多項式の積が高速に計算できます。ご興味のある方はぜひとも調べてみてください。

また実は FFT は再帰を使わずにキレイに書くことが出来ます。計算量オーダーこそ同じなのですが、そちらのほうが速いです。FFT で検索して見つかる解説記事では、そういうものを紹介していることが多いのですが、今回は「最速で実装にたどり着く」をテーマに据えたくて、あえて省略させていただきました。もしご興味のお有りの方が多数いらしましたら、続きの記事でご紹介できるかもしれません。