/*
 * A stand-alone implementation of the Fast Fourier Transformation (FFT).
 * Written 09/02/1999 by Mark-André Hopf <mhopf@mark13.de>.
 *
 * Note: I thought this might be useful to someone since all other
 *       FFT implementations need additional libraries.
 */

void fft(double *rbuffer, double *ibuffer, unsigned n)
{
	// presort
	unsigned j=0;
	unsigned b1 = n>>1;
	for(unsigned i=0; i<n; i++) {
		if (i<j) {
			double a=rbuffer[i]; rbuffer[i]=rbuffer[j]; rbuffer[j]=a;
		}
		unsigned m = b1;
		while(j&m) {
			j^=m;
			m>>=1;
		}
		j^=m;
	}

	// fft
	unsigned long mmax, m, i;
	double wr,wi;
	double tempr,tempi;
	
	mmax=1;
	while(mmax<n) {
		double t=0.0;
		double s=M_PI/mmax;
		m=0;
		while(m<mmax) {
			wr=cos(t);
			wi=sin(t);
			t+=s;
			i=m;
			while(i<n) {
				unsigned j=i; 
				i+=mmax;
				tempr=wr*rbuffer[i]-wi*ibuffer[i];
				tempi=wr*ibuffer[i]+wi*rbuffer[i];
				rbuffer[i]=rbuffer[j]-tempr;
				ibuffer[i]=ibuffer[j]-tempi;
				rbuffer[j]+=tempr;
				ibuffer[j]+=tempi;
				i+=mmax;
			}
			m++;
		}
		mmax<<=1;
	}
	return;
}