【OTFFT: High Speed FFT library】

【繰り返しによる Stockham FFT】

　前章、セクション1では Cooley-Tukey のアルゴリズムから単純な変形で Stockham のアルゴリズムを導きました。しかし、この再帰版のアルゴリズムは速度的に Cooley-Tukey に対して有利な点はなく、そのままでは用いる価値はありません。 そこでこの章ではさらなる書き換えを行って、繰り返しによる高速な Stockham のアルゴリズムを導きましょう。

　関数型プログラミング言語の信者の方は、 何でも再帰で書くのが一番みたいに思っているかもしれませんが、 世の中再帰で書いてあると気がつかないこともあります。 再帰バリバリで書かれている リスト３ のアルゴリズムを一部ループに書き換えてみましょう。 一つ一つフーリエ変換を呼び出していた部分をループで繰り返すようにします。

　リスト３ の関数 fft0(n,s,q,x,y), fft1(n,s,q,x,y) は x[q + s*p] をアクセスして処理します。 リスト３ は２分岐していく再帰ですので、 再帰の段数が進むと次のようになります。 

    再帰０段目
    fft0(8,1,0,x,y)
         x[p]

    再帰１段目
    fft1(4,2,0,x,y)                 fft1(4,2,1,x,y)
        x[2p+0]                         x[2p+1]

    再帰２段目
    fft0(2,4,0,x,y) fft0(2,4,2,x,y) fft0(2,4,1,x,y) fft0(2,4,3,x,y)
        x[4p+0]         x[4p+2]         x[4p+1]         x[4p+3]
図１：fft0 と fft1 の再帰の様子

　系列長 8 を例にとり、関数の下にアクセスする要素を表示しています。 再帰２段目 n = 2, s = 4 のときは q = 0,2,1,3 で処理することになります。 処理する順番が 0,1,2,3 ではないですが、 結局 x[4p],x[4p+1],x[4p+2],x[4p+3] 全てを処理することには変わりはありません。n, s は共通ですから q をループ変数として fft0(n,s,q,x,y) を 0 から s-1 まで(つまり、3 まで)回してやれば、 再帰２段目とループが等価になることが分かります。 この事実に基づいて リスト３ の再帰の一部をループに書き換えてやると以下のようになります。

リスト５：再帰のループによる書き換え 
#include <complex>
#include <cmath>

typedef std::complex<double> complex_t;

void fft1(int n, int s, complex_t* x, complex_t* y);

void fft0(int n, int s, complex_t* x, complex_t* y)
// n : 系列長
// s : ストライド
// x : フーリエ変換する系列(入出力)
// y : 作業用配列
{
    const int m = n/2;
    const double theta0 = 2*M_PI/n;

    if (n == 1) {}
    else {
        for (int q = 0; q < s; q++) { // 再帰の代わりに導入したループ
            for (int p = 0; p < m; p++) {
                const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
                const complex_t a = x[q + s*(p + 0)];
                const complex_t b = x[q + s*(p + m)];
                y[q + s*(2*p + 0)] =  a + b;
                y[q + s*(2*p + 1)] = (a - b) * wp;
            }
        }
        fft1(n/2, 2*s, y, x); // 再帰をループで束ねたので呼び出しは1つ
    }
}

void fft1(int n, int s, complex_t* x, complex_t* y)
// n : 系列長
// s : ストライド
// x : フーリエ変換する系列(入力)
// y : 作業用配列兼出力配列(出力)
{
    const int m = n/2;
    const double theta0 = 2*M_PI/n;

    if (n == 1) { for (int q = 0; q < s; q++) y[q] = x[q]; }
    else {
        for (int q = 0; q < s; q++) { // 再帰の代わりに導入したループ
            for (int p = 0; p < m; p++) {
                const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
                const complex_t a = x[q + s*(p + 0)];
                const complex_t b = x[q + s*(p + m)];
                y[q + s*(2*p + 0)] =  a + b;
                y[q + s*(2*p + 1)] = (a - b) * wp;
            }
        }
        fft0(n/2, 2*s, y, x); // 再帰をループで束ねたので呼び出しは1つ
    }
}

void fft(int N, complex_t* x) // フーリエ変換
// N : 系列長
// x : フーリエ変換する系列(入出力)
{
    complex_t* y = new complex_t[N]; // 並べ替え用作業配列の確保
    fft0(N, 1, x, y);
    delete[] y;
    for (int k = 0; k < N; k++) x[k] /= N;
}

void ifft(int N, complex_t* x) // 逆フーリエ変換
// N : 系列長
// x : 逆フーリエ変換する系列(入出力)
{
    for (int p = 0; p < N; p++) x[p] = conj(x[p]);
    complex_t* y = new complex_t[N]; // 並べ替え用作業配列の確保
    fft0(N, 1, x, y);
    delete[] y;
    for (int k = 0; k < N; k++) x[k] = conj(x[k]);
}

　どうでしょう。 ２分岐していく再帰が分岐しない再帰になってかなりスッキリしましたね。 これでもう一段、再帰をループに書き換えられることが分かると思います。しかし、 重要な点はそこではありません。そう、勘の良い人はもう気がついたと思いますが、 変数 p のループと変数 q のループは入れ替え可能です。 実はこのループを入れ替えると決定的な速度の差が現れます。

　再帰版 Stockham のアルゴリズムが遅いのは、再帰の段数が深くなるとストライド s が大きくなって、 配列を大きく飛び飛びにアクセスすることになるからです。 最近のコンピュータはこういうメモリのアクセスをすると速度が低下します。 いわゆるキャッシュに入り切らなくなると言うやつです。それでは、 ループを入れ替えてみましょう。次のようになります。

リスト６：ループを入れ替えた Stockham のアルゴリズム 
#include <complex>
#include <cmath>

typedef std::complex<double> complex_t;

void fft1(int n, int s, complex_t* x, complex_t* y);

void fft0(int n, int s, complex_t* x, complex_t* y)
// n : 系列長
// s : ストライド
// x : フーリエ変換する系列(入出力)
// y : 作業用配列
{
    const int m = n/2;
    const double theta0 = 2*M_PI/n;

    if (n == 1) {}
    else {
        for (int p = 0; p < m; p++) {
            const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
            for (int q = 0; q < s; q++) {
                const complex_t a = x[q + s*(p + 0)];
                const complex_t b = x[q + s*(p + m)];
                y[q + s*(2*p + 0)] =  a + b;
                y[q + s*(2*p + 1)] = (a - b) * wp;
            }
        }
        fft1(n/2, 2*s, y, x);
    }
}

void fft1(int n, int s, complex_t* x, complex_t* y)
// n : 系列長
// s : ストライド
// x : フーリエ変換する系列(入力)
// y : 作業用配列兼出力配列(出力)
{
    const int m = n/2;
    const double theta0 = 2*M_PI/n;

    if (n == 1) { for (int q = 0; q < s; q++) y[q] = x[q]; }
    else {
        for (int p = 0; p < m; p++) {
            const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
            for (int q = 0; q < s; q++) {
                const complex_t a = x[q + s*(p + 0)];
                const complex_t b = x[q + s*(p + m)];
                y[q + s*(2*p + 0)] =  a + b;
                y[q + s*(2*p + 1)] = (a - b) * wp;
            }
        }
        fft0(n/2, 2*s, y, x);
    }
}

void fft(int N, complex_t* x) // フーリエ変換
// N : 系列長
// x : フーリエ変換する系列(入出力)
{
    complex_t* y = new complex_t[N];
    fft0(N, 1, x, y);
    delete[] y;
    for (int k = 0; k < N; k++) x[k] /= N;
}

void ifft(int N, complex_t* x) // 逆フーリエ変換
// N : 系列長
// x : 逆フーリエ変換する系列(入出力)
{
    for (int p = 0; p < N; p++) x[p] = conj(x[p]);
    complex_t* y = new complex_t[N];
    fft0(N, 1, x, y);
    delete[] y;
    for (int k = 0; k < N; k++) x[k] = conj(x[k]);
}

　するとあら不思議、配列へのアクセスがストライド 1 になってしまいました。 ついでに三角関数の取得も一つ外側のループに移動できて、 さらに効率が上がっています。 そう、これこそが一般に言われる Stockham のアルゴリズムなのです。

　さて、このままでも良いのですが、 相互再帰でほとんど同じことを２回書くのは嫌な感じですよね。 そこで今どちらの関数なのかを表すフラグ eo を導入して通常の再帰に書き換えてみましょう。以下のようになります。

リスト７：完成版 Stockham のアルゴリズム 
#include <complex>
#include <cmath>

typedef std::complex<double> complex_t;

void fft0(int n, int s, bool eo, complex_t* x, complex_t* y)
// n  : 系列長
// s  : ストライド
// eo : eo == 0 か false なら x が出力、eo == 1 か true なら y が出力
// x  : フーリエ変換する入力系列(eo == 0 のとき出力)
// y  : 作業用配列(eo == 1 のとき出力)
{
    const int m = n/2;
    const double theta0 = 2*M_PI/n;

    if (n == 1) { if (eo) for (int q = 0; q < s; q++) y[q] = x[q]; }
    else {
        for (int p = 0; p < m; p++) {
            const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
            for (int q = 0; q < s; q++) {
                const complex_t a = x[q + s*(p + 0)];
                const complex_t b = x[q + s*(p + m)];
                y[q + s*(2*p + 0)] =  a + b;
                y[q + s*(2*p + 1)] = (a - b) * wp;
            }
        }
        fft0(n/2, 2*s, !eo, y, x);
    }
}

void fft(int N, complex_t* x) // フーリエ変換
// N : 系列長
// x : フーリエ変換する系列(入出力)
{
    complex_t* y = new complex_t[N];
    fft0(N, 1, 0, x, y);
    delete[] y;
    for (int k = 0; k < N; k++) x[k] /= N;
}

void ifft(int N, complex_t* x) // 逆フーリエ変換
// N : 系列長
// x : 逆フーリエ変換する系列(入出力)
{
    for (int p = 0; p < N; p++) x[p] = conj(x[p]);
    complex_t* y = new complex_t[N];
    fft0(N, 1, 0, x, y);
    delete[] y;
    for (int k = 0; k < N; k++) x[k] = conj(x[k]);
}

　次章、セクション3で、 応用問題として4基底の Stockham のアルゴリズムを Cooley-Tukey のアルゴリズムから導いてみましょう。