#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define FOREVER while(1){ (void)0; }
#define PI 3.141592654f
//this marco only usable by FFT_Compute() function assuming we already have var realComp, imagineComp defined
//four inputs: real is real component value, imagine is imaginary component value
//kn is muplication of time domain sample location n and frequency domain sample location k
// sum is total number of samples in current stage (stageNumber*2)
#define EXP_NEG_TWOPI_JAY_N_K_OVER_N(real, imagine, kn,sum) { realComp = real*cos(-2.0f*PI*kn/(float)sum)-imagine*sin(-2.0f*PI*kn/(float)sum); \
                                                              imagineComp = imagine*cos(-2.0f*PI*kn/(float)sum)+real*sin(-2.0f*PI*kn/(float)sum);}

typedef unsigned int UINT;

//given a M (M<32) bits number(0 - 2^M) revert bits
void reverseInPlace(UINT * input, UINT m)
{
	UINT out = 0, scrape = 0, k;
	for (k=m; k>=1; k--) {
		//right shift till see the leftmost bit
		scrape =(*input)>>(k-1) & 0x00000001;
		//out is assembled bit by bit from leftmost to rightmost with respect to original number
	        out += ((scrape == 0) ?  0 : (1<<(m-k)));
	}

	*input = out;
}

//FFT_Compute accepts two float array of length 2^m each represent Complex Number array in time domain, real component in realArray
//imaginary component in imagineArray, Output Arrays in frequency domain are stored in place repectively
//computational complexity O(nlogn), n is the length of arrays
void FFT_Compute(float *realArray, float *imagineArray, int m)
{
	//two temp variables used in computation
	float realComp, imagineComp;
	//length of our array is 2^m
	UINT len = 0x00000001 << (m);
	UINT k=0, j=0, temp, magicMask;
	//scrambleIndex array only used in stage 0
	UINT *scrambleIndex;
	//temp variable used for storing intermediate results
	float *realTemp, *imgTemp;

	if (m >= 32) {
		printf("We don't support bit length greater than 31");
		return;
	}

	scrambleIndex = (UINT *)malloc(sizeof(UINT) * len);
	realTemp = (float *)malloc(sizeof(float) * len);
	imgTemp = (float *)malloc(sizeof(float) * len);

	for (k=0; k<len; k++) {
		temp = k;
		reverseInPlace(&temp, m);
		scrambleIndex[k] =  temp;
	}
	//compute output for each stage, intermediate result stored in place
	//j is the stage number, k is the index to sample location

	//stage zero is a bit special cannot be written in general loop
	//the following loop just for stage 0
	for (k=0; k<len; k++) {
		realTemp[k] = realArray[scrambleIndex[k]];
		imgTemp[k] = imagineArray[scrambleIndex[k]];
		//here we only have zero
		/* dont need following since Weight is Identity Matrix
		EXP_NEG_TWOPI_JAY_N_K_OVER_N(realTemp[k], imgTemp[k], 0, len)

		if ((k&0x00000001) != 0) //odd number need to multiply weight
		{
			realTemp[k] = realComp;
			imgTemp[k] = imagineComp;
		} */
	} //copy to output array finish stage 0
	for (k=0; k<len; k++) {
		realArray[k] = realTemp[k];
		imagineArray[k] = imgTemp[k];
	}
	//now start from stage 1 til finish
	for (j=1; j<=m; j++) {
		//first sum then apply weight
		//magicMask to get butterfly location
		magicMask = (0x00000001<<(j-1));
		//compute each sample location in freq domain
		for (k=0; k<len; k++) {
			//adding counter part of sum instead of comparison (<) may also use & (eg. k&magicMask == 0? 1: -1)
			realTemp[k] = ( k < (k^magicMask) ? 1.f :-1.f) * realArray[k] + realArray[k^magicMask];
			imgTemp[k] = ( k < (k^magicMask) ? 1.f :-1.f) * imagineArray[k] + imagineArray[k^magicMask];

			//applying weight only odd subsequence needs weight
			//donot apply weight in last stage
			if (j!=m && (k > (k^(magicMask<<1)))) {
				//apply weight -- by left shift then right shift, we can zero out j leftmost bits
				EXP_NEG_TWOPI_JAY_N_K_OVER_N(realTemp[k], imgTemp[k], ((k<<(32-j))>>(32-j)), (magicMask<<2));
				realTemp[k] = realComp;
			        imgTemp[k] = imagineComp;
			}
		}
		//copy result over finish this stage
		for (k=0; k<len; k++) {
			realArray[k] = realTemp[k];
			imagineArray[k] = imgTemp[k];
	        }
	}

	free(scrambleIndex);
	free(realTemp);
	free(imgTemp);
}

/*
   Direct fourier transform for comparison purpose
   used for verify our Implementation is correct
*/
int DFT(int dir,int m,float *x1,float *y1)
{
   long i,k;
   float arg;
   float cosarg,sinarg;
   float *x2,*y2;

   x2 = (float *)malloc(m*sizeof(float));
   y2 = (float *)malloc(m*sizeof(float));
   if (x2 == NULL || y2 == NULL)
      return(-1);

   for (i=0;i<m;i++) {
      x2[i] = 0;
      y2[i] = 0;
      arg = - dir * 2.0 * 3.141592654 * (float)i / (float)m;
      for (k=0;k<m;k++) {
         cosarg = cos(k * arg);
         sinarg = sin(k * arg);
         x2[i] += (x1[k] * cosarg - y1[k] * sinarg);
         y2[i] += (x1[k] * sinarg + y1[k] * cosarg);
      }
   }

   /* Copy the data back */
   dir = 1;
   if (dir == 1) {
      for (i=0;i<m;i++) {
         x1[i] = x2[i];
         y1[i] = y2[i];
      }
   } else {
      for (i=0;i<m;i++) {
         x1[i] = x2[i];
         y1[i] = y2[i];
      }
   }

   free(x2);
   free(y2);
   return(0);
}

void main()
{
	float testReal[32], testImagine[32];
	int bitCount = 5, i, length = 32;
	for (i=0; i<length; i++) {
		testReal[i] = sin(2.0f*PI*(float)i/(float)length);
		testImagine[i] = 0;
	}

	printf("usage: %d(input number) %d(total bit count)\n");

	FFT_Compute(testReal, testImagine, bitCount);
	//DFT(1, length, testReal, testImagine);
	for (i=0; i<length; i++) {
		printf("Complex Number for sample %d: Real = %f Imaginary = %f\n", i, testReal[i], testImagine[i]);
	}

	FOREVER
}