#include <stdio.h>
#include <io.h>
#include <math.h>
#include <conio.h>
#include <fcntl.h>
#include <process.h>
#include <string.h>

#define MAX 1024
#define pi 3.14
#define DB double

void FastFourierTransform(DB Re[MAX], DB Im[MAX], int p);
void GenerateData (DB Re[MAX], DB Im[MAX], int p);
void DiscreteFourierTransform(DB DFTR[MAX], DB DFTI[MAX], DB Re[MAX], DB Im[MAX], int p);

void main ()
{
	int i, n, p = 8;
	DB Re[MAX], Im[MAX], s;
	DB DFTR[MAX], DFTI[MAX];

	n = 1 << p;
	system("cls");
	GenerateData (Re, Im, p);

	printf("#### result of DFT ####\n");
	DiscreteFourierTransform(DFTR, DFTI, Re, Im, p);
	FastFourierTransform(Re, Im, p);
	printf("\n\n");
	printf("##### result of FFT #####\n");

	for ( i = 1; i < n; i++){
    	s= sqrt(Re[i]*Re[i] + Im[i]*Im[i]);
	    printf("%3d : %8.3f\n", i, s);
		/*printf("%3d : %8.3f\t %8.3f\n", i, Re[i], Im[i]);*/
	}
	getch();
}

void FastFourierTransform(DB Re[MAX], DB Im[MAX], int p)
{
	int i, ip, j, k, m, me, me1, n, nv2;
	DB uRe, uIm, vRe, vIm, wRe, wIm, tRe, tIm;

	n = 1 << p;		//number of step
	nv2 = n / 2;
	j =0;			//initalize

	//part of bit reverse
	for (i = 0; i < n - 1; i++)
	{
		if (j > i)
		{
			tRe = Re[j];	tIm = Im[j];
			Re[j] = Re[i];	Im[j] = Im[i];
			Re[i] = tRe;	Im[i] = tIm;
		}

		k = nv2;
		while (j >= k) {
			j -= k;
			k /= 2;
		}
		j += k;
	}

	//butterfly loop
	for(m = 1; m <= p; m++)
	{
		me = 1 << m;	me1 = me / 2;
		uRe = 1.0;	uIm = 0.0;
		wRe = cos(pi / me1); wIm =- sin(pi / me1);

		for(j = 0; j < me1; j++)
		{
			for(i = j; i < n; i += me)
			{
				ip = i + me1;
				tRe = Re[ip] * uRe - Im[ip] * uIm;
				tIm = Re[ip] * uIm + Im[ip] * uRe;
				Re[ip] = Re[i] - tRe;
				Im[ip] = Im[i] - tIm;
				Re[i] += tRe;
				Im[i] += tIm;
			}
			vRe = uRe * wRe - uIm * wIm;
			vIm = uRe * wIm + uIm * wRe;
			uRe = vRe;	uIm = vIm;
		}
	}
}

void GenerateData (DB Re[MAX], DB Im[MAX], int p)
{
	int i;
	int n;

	n = 1 << p;

	for (i = 0; i < n; i++)
	{
		Re[i] = i; Im [i] =0;
	}
}

void DiscreteFourierTransform(DB DFTR[MAX], DB DFTI[MAX], DB Re[MAX], DB Im[MAX], int p)
{
	DB tp, s;
	int n;
	int i, j;

	tp = 2 * pi;
	n = 1 << p;

	for (i = 0; i < n; i++)
	{
		DFTR[i] = DFTI[i] = 0.0;
		for (j = 0; j < n; j++)
		{
			DFTR[i] += Re[j] * cos( - (tp / n) * i * j) - Im[j] * sin( - (tp / n) * i * j);
			DFTI[i] += Re[j] * sin( - (tp / n) * i * j) + Im[j] * cos( - (tp / n) * i * j);
		}
		s= sqrt(DFTR[i]*DFTR[i] + DFTI[i]*DFTI[i]);
		printf("%3d : %8.3f\n", i, s);
	}
}