[OTFFT: High Speed FFT library]
command
top page
init page
font spec
font small
font large
toc on
toc off
print
cancel
FontName/Size/Weight

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;
}
```