Chris@5: 
Chris@5: /* Public domain FFT implementation from Don Cross. */
Chris@5: 
Chris@5: #include "Cross.h"
Chris@5: 
Chris@5: #include <math.h>
Chris@5: #include <stdio.h>
Chris@5: 
Chris@5: void
Chris@5: fftCross(unsigned int n, int inverse,
Chris@5: 	 const double *ri, const double *ii,
Chris@5: 	 double *ro, double *io)
Chris@5: {
Chris@5:     if (!ri || !ro || !io) return;
Chris@5: 
Chris@5:     unsigned int bits;
Chris@5:     unsigned int i, j, k, m;
Chris@5:     unsigned int blockSize, blockEnd;
Chris@5: 
Chris@5:     double tr, ti;
Chris@5: 
Chris@5:     if (n < 2) return;
Chris@5:     if (n & (n-1)) return;
Chris@5: 
Chris@5:     double angle = 2.0 * M_PI;
Chris@5:     if (inverse) angle = -angle;
Chris@5: 
Chris@5:     for (i = 0; ; ++i) {
Chris@5: 	if (n & (1 << i)) {
Chris@5: 	    bits = i;
Chris@5: 	    break;
Chris@5: 	}
Chris@5:     }
Chris@5: 
Chris@5: #ifdef _MSC_VER
Chris@5:     int *table = (int *)_malloca(n * sizeof(int));
Chris@5: #else
Chris@5:     int table[n];
Chris@5: #endif
Chris@5: 
Chris@5:     for (i = 0; i < n; ++i) {
Chris@5:         m = i;
Chris@5:         for (j = k = 0; j < bits; ++j) {
Chris@5:             k = (k << 1) | (m & 1);
Chris@5:             m >>= 1;
Chris@5:         }
Chris@5:         table[i] = k;
Chris@5:     }
Chris@5: 
Chris@5:     if (ii) {
Chris@5: 	for (i = 0; i < n; ++i) {
Chris@5: 	    ro[table[i]] = ri[i];
Chris@5: 	    io[table[i]] = ii[i];
Chris@5: 	}
Chris@5:     } else {
Chris@5: 	for (i = 0; i < n; ++i) {
Chris@5: 	    ro[table[i]] = ri[i];
Chris@5: 	    io[table[i]] = 0.0;
Chris@5: 	}
Chris@5:     }
Chris@5: 
Chris@5:     blockEnd = 1;
Chris@5: 
Chris@5:     for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
Chris@5: 
Chris@5: 	double delta = angle / (double)blockSize;
Chris@5: 	double sm2 = -sin(-2 * delta);
Chris@5: 	double sm1 = -sin(-delta);
Chris@5: 	double cm2 = cos(-2 * delta);
Chris@5: 	double cm1 = cos(-delta);
Chris@5: 	double w = 2 * cm1;
Chris@5: 	double ar[3], ai[3];
Chris@5: 
Chris@5: 	for (i = 0; i < n; i += blockSize) {
Chris@5: 
Chris@5: 	    ar[2] = cm2;
Chris@5: 	    ar[1] = cm1;
Chris@5: 
Chris@5: 	    ai[2] = sm2;
Chris@5: 	    ai[1] = sm1;
Chris@5: 
Chris@5: 	    for (j = i, m = 0; m < blockEnd; j++, m++) {
Chris@5: 
Chris@5: 		ar[0] = w * ar[1] - ar[2];
Chris@5: 		ar[2] = ar[1];
Chris@5: 		ar[1] = ar[0];
Chris@5: 
Chris@5: 		ai[0] = w * ai[1] - ai[2];
Chris@5: 		ai[2] = ai[1];
Chris@5: 		ai[1] = ai[0];
Chris@5: 
Chris@5: 		k = j + blockEnd;
Chris@5: 		tr = ar[0] * ro[k] - ai[0] * io[k];
Chris@5: 		ti = ar[0] * io[k] + ai[0] * ro[k];
Chris@5: 
Chris@5: 		ro[k] = ro[j] - tr;
Chris@5: 		io[k] = io[j] - ti;
Chris@5: 
Chris@5: 		ro[j] += tr;
Chris@5: 		io[j] += ti;
Chris@5: 	    }
Chris@5: 	}
Chris@5: 
Chris@5: 	blockEnd = blockSize;
Chris@5:     }
Chris@5: 
Chris@5:     if (inverse) {
Chris@5: 
Chris@5: 	double denom = (double)n;
Chris@5: 
Chris@5: 	for (i = 0; i < n; i++) {
Chris@5: 	    ro[i] /= denom;
Chris@5: 	    io[i] /= denom;
Chris@5: 	}
Chris@5:     }
Chris@5: 
Chris@5: #ifdef _MSC_VER
Chris@5:     _freea(table);
Chris@5: #endif
Chris@5: }
Chris@5: