【Six-Step FFT と Eight-Step FFT】

　このページでは、Six-Step FFT アルゴリズムを使った FFT の最適化について説明します。

　まず、離散フーリエ変換について定義しましょう。 入力 \(x_p\) から出力 \(X_k\) を得るサイズ \(N = 2^L,~L=0,1,2,3,\ldots\) の離散フーリエ変換 \(F_N\) は、以下のように定義されます。

　今、\(N = nm\) と因数分解できると仮定して、\(x_p\) を配列 \(x[p]\) で表すと、以下のように \(p_1, p_2\) を用いて２次元的に表現できます。

　ここで、\(p_1 = 0,1,\ldots,n-1\) で \(p_2 = 0,1,\ldots,m-1\) です。

　同様に、\(X_k\) を配列 \(X[k]\) で表すと、以下のように \(k_1, k_2\) を用いて２次元的に表現できます。

　ここで、\(k_1 = 0,1,\ldots,m-1\) で \(k_2 = 0,1,\ldots,n-1\) です。

　今、以下の関係が成り立つことを利用して \(F_N\) を変形すると、

\(F_N\) は以下のような２重の離散フーリエ変換 \(G_N\) に変形できます。

　次に、いくつか用語を定義しましょう。\(x[p_1 + p_2n]\) の \(p_2\) を \(p_2 = q\) に固定した場合、 \(x\) は \(p_1\) のみを添字とする配列と見なせます。 これを \(x\) の \(q\) 行と呼ぶことにします。また、次の関係 \(x[p_1 + p_2n] = y[p_2 + p_1m]\) を満たす \(y\) を \(x\) の転置配列と呼びます。\(y\) を求めることを \(x\) を転置すると言います。 これらの用語を使うと、\(G_N\) は以下の手順で計算されます。

　転置が必要なのは、 \(F_n, F_m\) などをプログラムで実装されたサブルーチンと考えた場合、 配列の添字のストライドが１でないとうまく変換できないと仮定しているからです。 このように、６個のステップで変換できることから、 このアルゴリズムは Six-Step FFT と呼ばれています。 もちろん、\(F_n, F_m\) は FFT で実装します。

　Six-Step FFT にしても、計算量は通常の FFT と変わりません。 では、Six-Step FFT にすると何が嬉しいのでしょうか。 通常の FFT は、 系列のサイズが大きくなるとメモリアクセスの局所性が低下します。しかし、 Six-Step FFT の場合、 大きな系列を複数の小さな系列に分解して FFT するので、 メモリアクセスの局所性が改善されます。 そのため、大きなサイズの系列で Six-Step FFT は性能が向上します。

　Six-Step FFT を実装してみると以下のようになります。

#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 == 2) {
        complex_t* z = eo ? y : x;
        for (int q = 0; q < s; q++) {
            const complex_t a = x[q + 0];
            const complex_t b = x[q + s];
            z[q + 0] = a + b;
            z[q + s] = a - b;
        }
    }
    else if (n >= 4) {
        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 sixstep_fft(int log_N, complex_t* x, complex_t* y)
// log_N : 系列長の対数(２以上の偶数)
// x     : フーリエ変換する系列(入出力)
// y     : 作業用配列
{
    const int N = 1 << log_N;
    const int n = 1 << (log_N/2); // N == n*n とする
    for (int k = 0; k < n; k++) { // x を転置する
        for (int p = k+1; p < n; p++) {
            const complex_t a = x[p+k*n];
            const complex_t b = x[k+p*n];
            x[p+k*n] = b;
            x[k+p*n] = a;
        }
    }
    for (int p = 0; p < n; p++) // x の全ての p 行をフーリエ変換する
        fft0(n, 1, 0, x + p*n, y + p*n);
    for (int p = 0; p < n; p++) { // ひねり係数を乗算し、x を転置する
        const double theta0 = 2*p*M_PI/N;
        for (int k = p; k < n; k++) {
            const double theta = k*theta0;
            const complex_t wkp = complex_t(cos(theta), -sin(theta));
            if (k == p)
                x[p+p*n] *= wkp;
            else {
                const complex_t a = x[k+p*n] * wkp;
                const complex_t b = x[p+k*n] * wkp;
                x[k+p*n] = b;
                x[p+k*n] = a;
            }
        }
    }
    for (int k = 0; k < n; k++) // x の全ての k 行をフーリエ変換する
        fft0(n, 1, 0, x + k*n, y + k*n);
    for (int k = 0; k < n; k++) { // x を転置する
        for (int p = k+1; p < n; p++) {
            const complex_t a = x[p+k*n];
            const complex_t b = x[k+p*n];
            x[p+k*n] = b;
            x[k+p*n] = a;
        }
    }
}

void eightstep_fft(int log_n, complex_t* x, complex_t* y)
// log_n : 系列長の対数(３以上の奇数)
// x     : フーリエ変換する系列(入出力)
// y     : 作業用配列
{
    const int n = 1 << log_n;
    const int m = n/2;
    const double theta0 = 2*M_PI/n;
    for (int p = 0; p < m; p++) {
        const double theta = p*theta0;
        const complex_t wp = complex_t(cos(theta), -sin(theta));
        const complex_t a = x[p+0];
        const complex_t b = x[p+m];
        y[p+0] =  a + b;
        y[p+m] = (a - b) * wp;
    }
    sixstep_fft(log_n-1, y + 0, x + 0);
    sixstep_fft(log_n-1, y + m, x + m);
    for (int p = 0; p < m; p++) {
        x[2*p+0] = y[p+0];
        x[2*p+1] = y[p+m];
    }
}

void fft(int N, complex_t* x) // フーリエ変換
// N : 系列長
// x : フーリエ変換する系列(入出力)
{
    int log_N = 0;
    for (int i = N; i > 1; i >>= 1) log_N++;
    complex_t* y = new complex_t[N];
    if (N <= 1) {}
    else if (N == 2)       fft0(N, 1, 0, x, y);
    else if (!(log_N & 1)) sixstep_fft(log_N, x, y);
    else                   eightstep_fft(log_N, 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]);
    int log_N = 0;
    for (int i = N; i > 1; i >>= 1) log_N++;
    complex_t* y = new complex_t[N];
    if (N <= 1) {}
    else if (N == 2)       fft0(N, 1, 0, x, y);
    else if (!(log_N & 1)) sixstep_fft(log_N, x, y);
    else                   eightstep_fft(log_N, x, y);
    delete[] y;
    for (int k = 0; k < N; k++) x[k] = conj(x[k]);
}

　系列長 \(N\) が \(N = nm\) と因数分解されると仮定しましたが、 この時 \(n = m\) と仮定すると転置の計算が極めて効率的に実行できます。 また、FFT の係数 \(W_N\) の読み込み回数が通常の半分になります。 そのため、このコードでは、 sixstep_fft() サブルーチンにおいて \(n = m\) と仮定しています。

　\(n = m\) にならない時は、\(n = \frac{N}{2}, m = 2\) と置いて、 Six-Step FFT のアルゴリズを適用して分解すると、\(n = uv, u = v\) となって、 \(F_n\) が sixstep_fft() サブルーチンを用いて実行できます。 これを私は Eight-Step FFT と呼んでいます。

　ところで、ここで紹介したアルゴリズムを用いた FFT ライブラリを作ってみました。 興味のある方は、こちらのページ へどうそ。