/****************************************************
	FFT()

	參數(shù)：

		TD為時域值
		FD為頻域值
		power為2的冪數(shù)

	返回值：

		無

	說明：

		本函數(shù)實現(xiàn)快速傅立葉變換
****************************************************/
void FFT(COMPLEX * TD, COMPLEX * FD, int power)
{
	int count;
	int i,j,k,bfsize,p;
	double angle;
	COMPLEX *W,*X1,*X2,*X;

	/*計算傅立葉變換點數(shù)*/
	count=1<<power;
	
	/*分配運(yùn)算所需存儲器*/
	W=(COMPLEX *)malloc(sizeof(COMPLEX)*count/2);
	X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	
	/*計算加權(quán)系數(shù)*/
	for(i=0;i<count/2;i++)
	{
		angle=-i*PI*2/count;
		W[i].re=cos(angle);
		W[i].im=sin(angle);
	}
	
	/*將時域點寫入存儲器*/
	memcpy(X1,TD,sizeof(COMPLEX)*count);
	
	/*蝶形運(yùn)算*/
	for(k=0;k<power;k++)
	{
		for(j=0;j<1<<k;j++)
		{
			bfsize=1<<(power-k);
			for(i=0;i<bfsize/2;i++)
			{
				p=j*bfsize;
				X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);
				X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);
			}
		}
		X=X1;
		X1=X2;
		X2=X;
	}
	
	/*重新排序*/
	for(j=0;j<count;j++)
	{
		p=0;
		for(i=0;i<power;i++)
		{
			if (j&(1<<i)) p+=1<<(power-i-1);
		}
		FD[j]=X1[p];
	}
	
	/*釋放存儲器*/
	free(W);
	free(X1);
	free(X2);
}


/****************************************************
	IFFT()

	參數(shù)：

		FD為頻域值
		TD為時域值
		power為2的冪數(shù)

	返回值：

		無

	說明：

		本函數(shù)利用快速傅立葉變換實現(xiàn)傅立葉反變換
****************************************************/
void IFFT(COMPLEX * FD, COMPLEX * TD, int power)
{
	int i, count;
	COMPLEX *x;

	/*計算傅立葉反變換點數(shù)*/
	count=1<<power;

	/*分配運(yùn)算所需存儲器*/
	x=(COMPLEX *)malloc(sizeof(COMPLEX)*count);

	/*將頻域點寫入存儲器*/
	memcpy(x,FD,sizeof(COMPLEX)*count);
	
	/*求頻域點的共軛*/
	for(i=0;i<count;i++)
		x[i].im = -x[i].im;

	/*調(diào)用FFT*/
	FFT(x, TD, power);

	/*求時域點的共軛*/
	for(i=0;i<count;i++)
	{
		TD[i].re /= count;
		TD[i].im = -TD[i].im / count;
	}

	/*釋放存儲器*/
	free(x);
}