#include <stdio.h>
#include <math.h>
/*routines from Baase for Fast Fourier Transform*/
typedef char BOOLEAN;
#define TRUE 1
#define FALSE 0

typedef struct {
	float real;
	float imag;
	} complex;
typedef complex ComplexArray[512];
typedef float RealArray[512];
ComplexArray omega;  /*nth roots of unity*/
ComplexArray transform;  /*the transform*/
short pi[512];  /*bit reversal permutation*/
short k;   /*2**k elements used in array*/
short i, n;
ComplexArray test,tform;

main()

{
	test1();
}

short twoPower(k)
short k;
{
short twoP;
short i;

twoP = 1;
for (i=1;i<=k;i++)
	twoP *=  2;
return twoP;
}

initPi()
/* stores reversed bits for FFT*/
{
BOOLEAN	bit[12];
short val1, val2, i, j;
short posValue[12];

	posValue[1] = 1;
	for (i=2;i<=(k+1);i++)
		posValue[i] = 2 * posValue[i - 1];
	for (i=1;i<=(k+1);i++)
		bit[i] = FALSE;
	val1 = 0;
	val2 = 0;
	while (!bit[k + 1])
	{
		pi[val1] = val2;
		i = 1;
		while (bit[i])
		{
			bit[i] = FALSE;
			val1 -= posValue[i];
			val2 -= posValue[k + 1 - i];
			i++;
		}
		bit[i] = TRUE;
		val1 += posValue[i];
		val2 += posValue[k + 1 - i];
	}
}


complex addComplex (x,y)
complex x,y;
{
	complex aC;

	aC.real = x.real + y.real;
	aC.imag = x.imag + y.imag;
	return aC;
}

complex subtractComplex (x,y)
complex x,y;
{
	complex sC;

	sC.real = x.real - y.real;
	sC.imag = x.imag - y.imag;
	return sC;
}

complex multiplyComplex (x,y)
complex x,y;
{
	complex mC;

	mC.real = x.real * y.real - x.imag * y.imag;
	mC.imag = x.real * y.imag + x.imag * y.real;
	return mC;
}

initOmega()
{
short i;

omega[0].real = cos(0.0);
omega[0].imag = sin(0.0);
omega[1].real = cos(6.2831854/n);
omega[1].imag = sin(6.2831854/n);
for (i=2;i<n;i++)
	omega[i]=multiplyComplex(omega[i-1],omega[1]);
}

FFT(P,n,transform)
ComplexArray P,transform;
short n;
{
	short l;    /*level number*/
	short num; /*the number of values to be computed at each node of level l*/
	short t; /*the index in transform for the first of these values for a particular node*/
	short j;  /*counts off the pairs of values to be computed for that node*/
	short m;  /*used to pick out correct entry from omega*/
	short p1, power;
	complex xPOdd;

 /*to intialize transform, evaluate polynomials of degree l=2**k-1 at square roots of 1.*/
	for (t=0;t <= (n - 2);t+=2) {
		transform[t]=addComplex(P[pi[t]],P[pi[t + 1]]);
		transform[t + 1]=subtractComplex(P[pi[t]], P[pi[t + 1]]);
	}

    /*The main computation*/
	m = n / 2;
	num = 2;

    /*Begin triply nested loop*/
	for (l=k-2;l>=0;l--) {
		m /= 2;
		num *= 2;
		t = 0;

		power = 1;
		for (p1=1;p1<=l;p1++)
			power *= 2;
		power = (power - 1) * num;

		for (t=0;t <= power;t += num) {
			for (j=0;j<=((num/2)-1);j++) {
				xPOdd=multiplyComplex(omega[m*j],transform[t+(num/2)+j]);
				transform[t+(num/2)+j]=subtractComplex(transform[t+j],xPOdd);
				transform[t+j]=addComplex(transform[t+j],xPOdd);
			}
		}
	}
}

InverseFT (A,n,B)
ComplexArray A,B;
short n;
{
	short i;
	ComplexArray transform;

	FFT(A,n,transform);
	B[0].real = transform[0].real / n;
	B[0].imag = transform[0].imag / n;
	for (i=1;i<=(n - 1);i++) {
		B[i].real = transform[n - i].real / n;
		B[i].imag = transform[n - i].imag / n;
	}
}

test1()
{
	short i;

	k = 9;
	n = 512;

	initPi();

	initOmega();

	for (i=0;i<=511;i++) {
		test[i].real =  2*i;
		test[i].imag = -i;
	}
	
	printf("FFT input\n");
	for (i=0;i<=511;i++)
		printf("%e %e \n",test[i].real,test[i].imag);

	FFT(test,512,tform);

	printf("FFT output\n");
	for (i=0;i<=511;i++)
		printf("%e %e \n",tform[i].real,tform[i].imag);

	InverseFT(tform,512,test);

	printf("inverse FT output\n");
	for (i=0;i<=511;i++)
		printf("%e %e \n",test[i].real,test[i].imag);
}