#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include <complex>
using namespace std;
//*************************************************************************
// *
// * 函數名稱：
// *   FFT()
// *
// * 參數:
// *   complex<double> * TD- 指向時域數組的指針
// *   complex<double> * FD- 指向頻域數組的指針
// *   r－2的冪數，即迭代次數
// *
// * 返回值:
// *   無。
// *
// * 說明:
// *   該函數用來實現快速付立葉變換。
// *
// ************************************************************************/
void FFT(	complex<double> *TD,  complex<double> *FD, 
			complex<double> *X1,  complex<double> *X2,
			complex<double> *Omega,	  int r)
{
	// 付立葉變換點數
	long count;

	// 循環變量
	int i,j,k;

	// 中間變量
	int bfsize,p;

	// 角度
	double angle;

	complex<double> *X;

	// 計算付立葉變換點數
	count = 1 << r;

	// 分配運算所需存儲器
	//Omega  = new complex<double>[count / 2];
	//X1 = new complex<double>[count];
	//X2 = new complex<double>[count];

	// 計算加權系數
	for(i = 0; i < count / 2; i++)
	{
		angle = -i * 3.1415926 * 2 / count;
		Omega[i] = complex<double>(cos(angle), sin(angle));
	}

	// 將時域點寫入X1
	memcpy(X1, TD, sizeof(complex<double>) * count);

	// 采用蝶形算法進行快速付立葉變換
	for(k = 0; k < r; k++)
	{
		for(j = 0; j < 1 << k; j++)
		{
			bfsize = 1 << (r-k);
			for(i = 0; i < bfsize / 2; i++)
			{
				p = j * bfsize;
				X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
				X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * Omega[i * (1<<k)];
			}
		}
		X  = X1;
		X1 = X2;
		X2 = X;
	}

	// 重新排序
	for(j = 0; j < count; j++)
	{
		p = 0;
		for(i = 0; i < r; i++)
		{
			if (j&(1<<i))
			{
				p+=1<<(r-i-1);
			}
		}
		FD[j]=X1[p];
	}

	// 釋放內存
	//delete Omega;
	//delete X1;
	//delete X2;
}