[comp.lang.c++] mixed radix FFT in C++

christ@tybalt.caltech.edu (Christian L. Keppenne) (08/16/89)

Does someone have a C++ implementation of the general FFT algorithm.
I would prefer the mixed radix case (number of points not necessarily
a power of two) but the standard Cooley Tukey algorithm would do too.
I am sure this must be available somewhere and I will be grateful
to the person who will help me save the time of reprogramming it 
myself.

Thank you in advance.

Christian Keppenne 	christ&tybalt.caltech.edu

jima@hplsla.HP.COM (Jim Adcock) (08/18/89)

#include <complex.h>
#include <stream.h>

// fft a double precision complex array x of size (2 to the m)
// this is like I have posted before, but contains a bug workaround
// necessary for some C++ compilers.

int fft(complex* x, unsigned m)
{
    complex u, w, t;
    unsigned n = 1 << m;
    unsigned nv2 = n >> 1;
    unsigned nm1 = n - 1;
    unsigned j = 0;

    for (unsigned i = 0; i < nm1; ++i)
    {
	if (i < j) {t = x[j]; x[j] = x[i]; x[i] = t;}
	unsigned k = nv2;
	while ((k-1) < j) {j -= k; k >>= 1;}
	j += k;
    }

    for (unsigned l = 1; l <= m; ++l)
    {
        unsigned le = 1 << l;
        unsigned le1 = le >> 1;
        complex u(1.0, 0.0);

        double cs = cos(M_PI/double(le1)); //temporaries introduced to fix bug
	double sn = sin(M_PI/double(le1));

        complex w(cs,sn);

        for (j = 0; j < le1; ++j)
        {
	    for (i = j; i < n; i += le)
	    {
		unsigned ip = i + le1;
		t = x[ip] * u;
		x[ip] = x[i] - t;
		x[i] += t;
            }
	    u *= w;
        }
    }
    return m;
}


#define SZ 16 
#define LGSZ 4

main()
{
    complex x[SZ];
    for (long i = 0; i < SZ; ++i)
    {
	for (long j = 0; j < SZ; ++j) x[j] = 0.0;
	x[i] = 1.0;
	fft(x, LGSZ);
	fft(x, LGSZ);   //two ffts in a row without conjugation causes time
			//reversal and scaling by SZ
	for (j = 0; j < SZ; ++j) 
	{
	  //to make a pretty output force small numbers to zero
	  if (abs(real(x[j])) < 1.0e-10) x[j] = complex(0.0,imag(x[j]));
	  if (abs(imag(x[j])) < 1.0e-10) x[j] = complex(real(x[j]),0.0);

	  cout << "(" << real(x[j]) << "+i" << imag(x[j]) << ")" << " ";
        }
	cout << "\n\n";
    }
}