【Stockham FFT の種類】

　ここでは、私が Web を検索して見つけた Stockham のアルゴリズムの種類について整理しておきます。まず、 このサイトで説明されている Cooley-Tukey のアルゴリズムからの単純な変形で導ける Stockham のアルゴリズムをタイプ１。 私がなぜこれでフーリエ変換できるのか説明できない Stockham のアルゴリズムをタイプ２とします。

　タイプ１、 タイプ２それぞれに周波数間引き(DIF)と時間間引き(DIT)のアルゴリズムがあります。 すなわち、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)
{
    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)
{
    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)
{
    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]);
}

#include <complex>
#include <cmath>

typedef std::complex<double> complex_t;

void fft0(int n, int s, bool eo, complex_t* x, complex_t* y)
{
    const int m = n/2;
    const double theta0 = 2*M_PI/n;

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

void fft(int N, complex_t* 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)
{
    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]);
}

#include <complex>
#include <cmath>

typedef std::complex<double> complex_t;

void fft0(int n, int s, bool eo, complex_t* x, complex_t* y)
{
    const int m = n/2;
    const double theta0 = M_PI/s;

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

void fft(int N, complex_t* 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)
{
    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]);
}

#include <complex>
#include <cmath>

typedef std::complex<double> complex_t;

void fft0(int n, int s, bool eo, complex_t* x, complex_t* y)
{
    const int m = n/2;
    const double theta0 = M_PI/s;

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

void fft(int N, complex_t* 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)
{
    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]);
}