OTFFT: High Speed FFT library

[Radix-4 Stockham FFT]

Do you know the Radix-4 Stockham Algorithm? Radix-4 Stockham Algorithm is faster than Radix-2 Stockham Algorithm. In this page, I show the program of Radix-4 Stockham FFT. Please refer to "Introduction to the Stockham FFT" for details of the Stockham algorithm.

#include <complex>
#include <cmath>

typedef std::complex<double> complex_t;

void fft0(int n, int s, bool eo, complex_t* x, complex_t* y)
{
    static const complex_t j = complex_t(0, 1);
    const int n0 = 0;
    const int n1 = n/4;
    const int n2 = n/2;
    const int n3 = n1 + n2;
    const double theta0 = 2*M_PI/n;

    if (n == 1) { if (eo) for (int q = 0; q < s; q++) y[q] = x[q]; }
    else if (n == 2) {
        for (int q = 0; q < s; q++) {
            const complex_t a = x[q + 0];
            const complex_t b = x[q + s];
            y[q + 0] = a + b;
            y[q + s] = a - b;
        }
        fft0(1, 2*s, !eo, y, x);
    }
    else if (n > 2) {
        for (int p = 0; p < n1; p++) {
            const complex_t w1p = complex_t(cos(p*theta0), -sin(p*theta0));
            const complex_t w2p = w1p*w1p;
            const complex_t w3p = w1p*w2p;
            for (int q = 0; q < s; q++) {
                const complex_t a = x[q + s*(p + n0)];
                const complex_t b = x[q + s*(p + n1)];
                const complex_t c = x[q + s*(p + n2)];
                const complex_t d = x[q + s*(p + n3)];
                const complex_t  apc =    a + c;
                const complex_t  amc =    a - c;
                const complex_t  bpd =    b + d;
                const complex_t jbmd = j*(b - d);
                y[q + s*(4*p + 0)] =      apc +  bpd;
                y[q + s*(4*p + 1)] = w1p*(amc - jbmd);
                y[q + s*(4*p + 2)] = w2p*(apc -  bpd);
                y[q + s*(4*p + 3)] = w3p*(amc + jbmd);
            }
        }
        fft0(n/4, 4*s, !eo, y, x);
    }
}

void ifft0(int n, int s, bool eo, complex_t* x, complex_t* y)
{
    static const complex_t j = complex_t(0, 1);
    const int n0 = 0;
    const int n1 = n/4;
    const int n2 = n/2;
    const int n3 = n1 + n2;
    const double theta0 = 2*M_PI/n;

    if (n == 1) { if (eo) for (int q = 0; q < s; q++) y[q] = x[q]; }
    else if (n == 2) {
        for (int q = 0; q < s; q++) {
            const complex_t a = x[q + 0];
            const complex_t b = x[q + s];
            y[q + 0] = a + b;
            y[q + s] = a - b;
        }
        ifft0(1, 2*s, !eo, y, x);
    }
    else if (n > 2) {
        for (int p = 0; p < n1; p++) {
            const complex_t w1p = complex_t(cos(p*theta0), sin(p*theta0));
            const complex_t w2p = w1p*w1p;
            const complex_t w3p = w1p*w2p;
            for (int q = 0; q < s; q++) {
                const complex_t a = x[q + s*(p + n0)];
                const complex_t b = x[q + s*(p + n1)];
                const complex_t c = x[q + s*(p + n2)];
                const complex_t d = x[q + s*(p + n3)];
                const complex_t  apc =    a + c;
                const complex_t  amc =    a - c;
                const complex_t  bpd =    b + d;
                const complex_t jbmd = j*(b - d);
                y[q + s*(4*p + 0)] =      apc +  bpd;
                y[q + s*(4*p + 1)] = w1p*(amc + jbmd);
                y[q + s*(4*p + 2)] = w2p*(apc -  bpd);
                y[q + s*(4*p + 3)] = w3p*(amc - jbmd);
            }
        }
        ifft0(n/4, 4*s, !eo, y, x);
    }
}

void fft(int n, complex_t* x) // Fourier transform
// n : sequence length
// x : input/output sequence
{
    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) // Inverse Fourier transform
// n : sequence length
// x : input/output sequence
{
    complex_t* y = new complex_t[n];
    ifft0(n, 1, 0, x, y);
    delete[] y;
}