fft-64bit

``c
/\*
A set of utility programs to compute the Fast Fourier Transform (FFT):

N-1
X[k] = SUM { x[n]exp(-j2 pi nk/N) }
n=0

and inverse Fast Fourier Transform (iFFT):

N-1
x[n] = 1/N SUM { X[k]exp(+j2 pi nk/N) }
k=0

To speed things up, multiplication by 1 and j are avoided. The input, x[],

is an array of complex numbers (pairs of doubles) of length N = 2^n. The
calling program supplies n = log2(N) not the array length, N. The output is
placed in BIT REVERSED order in x[]. Call bit_reverse(x, n) to swap back to
normal order, if needed. However, iFFT(X, n, stbl) requires its input, X[],
to be in bit reversed order making bit reversing unnecessary in some cases,
such as convolution. stbl is an array of doubles of length N/4 containing
the sine function from 0 to pi/2 used to compute the FFT. Call
sin_table(stbl, n) ONLY ONCE before either FFT(x, n, stbl) or iFFT(X, n,
stbl) to set up a sin table for FFT computation.

Written ca. 1985 in THINK C by Robert Bristow-Johnson.

\*/

#define HALFPI 1.570796326794897
#define PI 3.141592653589793
#define TWOPI 6.283185307179586

// #include "complex.h"

typedef struct {
double real;
double imag;
} complex;

#define Re(z) (z).real
#define Im(z) (z).imag

void FFT(complex x, int n, double stbl);
void iFFT(complex x, int n, double stbl);
void bit_reverse(register complex *x, int n);
void sin_table(double *stbl, int n);
void FFT(complex x, int n, double stbl) {
long size;
register long length, step, stepsize, end;
register complex top, bottom; / top & bottom of FFT butterfly /
complex temp;

size = 1L << (n - 2);
end = (long)x + 4 sizeof(complex) size;

length = size;
stepsize = 1;
while (length >= 1) {
top = x;
while ((long)top < end) {
bottom = top + 2 * length;

Re(temp) = Re(top) - Re(bottom); / butterfly: twiddle = 1 /
Im(temp) = Im(top) - Im(bottom);
Re(top) += Re(bottom);
Im(top) += Im(bottom);
*bottom = temp;
top++;
bottom++;

for (step = stepsize; step < size; step += stepsize) {
Re(temp) =
Re(top) - Re(bottom); / butterfly: twiddle = exp(-j theta) /
Im(temp) = Im(top) - Im(bottom);
Re(top) += Re(bottom);
Im(top) += Im(bottom);
Re(bottom) = Re(temp) stbl[size - step] + Im(temp) * stbl[step];
Im(bottom) = Im(temp) stbl[size - step] - Re(temp) * stbl[step];
top++;
bottom++;
}

Re(temp) = Im(top) - Im(bottom); / butterfly: twiddle = -j /
Im(temp) = Re(bottom) - Re(top);
Re(top) += Re(bottom);
Im(top) += Im(bottom);
*bottom = temp;
top++;
bottom++;

for (step = stepsize; step < size; step += stepsize) {
Re(temp) =
Im(top) - Im(bottom); / butterfly: twiddle = -jexp(-j theta) */
Im(temp) = Re(bottom) - Re(top);
Re(top) += Re(bottom);
Im(top) += Im(bottom);
Re(bottom) = Re(temp) stbl[size - step] + Im(temp) * stbl[step];
Im(bottom) = Im(temp) stbl[size - step] - Re(temp) * stbl[step];
top++;
bottom++;
}
top = bottom;
}
length >>= 1;
stepsize <<= 1;
}

top = x;
bottom = x + 1;
while ((long)top < end) {
Re(temp) = Re(top) - Re(bottom); / butterfly: twiddle = 1 /
Im(temp) = Im(top) - Im(bottom);
Re(top) += Re(bottom);
Im(top) += Im(bottom);
*bottom = temp;
top += 2;
bottom += 2;
}
}

void iFFT(complex X, int n, double stbl) {
long size;
register long length, step, stepsize, end;
double scale;
register complex top, bottom; / top & bottom of FFT butterfly /
complex temp;

size = 1L << (n - 2);
end = (long)X + 4 sizeof(temp) size;

scale = 0.25 / size;
top = X;
bottom = X + 1;
while ((long)top < end) {
Re(temp) = (Re(top) - Re(bottom)) scale; / butterfly: twiddle = 1/N */
Im(temp) = (Im(top) - Im(bottom)) * scale;
Re(top) = (Re(top) + Re(bottom)) scale;
Im(top) = (Im(top) + Im(bottom)) scale;
*bottom = temp;
top += 2;
bottom += 2;
}

length = 1;
stepsize = size;
while (stepsize >= 1) {
top = X;
while ((long)top < end) {
bottom = top + 2 * length;

temp = bottom; / butterfly: twiddle = 1 */
Re(bottom) = Re(top) - Re(temp);
Im(bottom) = Im(top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;

for (step = stepsize; step < size;
step += stepsize) { / butterfly: twiddle = exp(+j theta) /
Re(temp) = Re(bottom) stbl[size - step] - Im(bottom) stbl[step];
Im(temp) = Im(bottom) stbl[size - step] + Re(bottom) stbl[step];
Re(bottom) = Re(top) - Re(temp);
Im(bottom) = Im(top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;
}

Re(temp) = -Im(bottom); / butterfly: twiddle = +j */
Im(temp) = Re(*bottom);
Re(bottom) = Re(top) - Re(temp);
Im(bottom) = Im(top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;

for (step = stepsize; step < size;
step += stepsize) { / butterfly: twiddle = +jexp(+j theta) */
Re(temp) = -Im(bottom) stbl[size - step] - Re(bottom) stbl[step];
Im(temp) = Re(bottom) stbl[size - step] - Im(bottom) stbl[step];
Re(bottom) = Re(top) - Re(temp);
Im(bottom) = Im(top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;
}
top = bottom;
}
length <<= 1;
stepsize >>= 1;
}
}

void bit_reverse(register complex *x, int n) {
complex temp;
register long k, i, r, size, count;

size = (1L << n) - 1L;
for (k = 1L; k < size; k++) {
i = k;
r = 0;
for (count = n; count > 0; count--) {
r <<= 1;
r += (i & 0x00000001L);
i >>= 1;
}
if (r > k) {
temp = x[r];
x[r] = x[k];
x[k] = temp;
}
}
}
``