#include <math.h>
#include <malloc.h>
#include <string.h>

#define pi (double) 3.14159265359

/*復(fù)數(shù)定義*/
typedef struct
{
	double re;
	double im;
}COMPLEX;

/*復(fù)數(shù)加法運算*/
COMPLEX Add(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re+c2.re;
	c.im=c1.im+c2.im;
	return c;
}

/*復(fù)數(shù)減法運算*/
COMPLEX Sub(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re-c2.re;
	c.im=c1.im-c2.im;
	return c;
}

/*復(fù)數(shù)乘運算*/
COMPLEX Mul(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re*c2.re-c1.im*c2.im;
	c.im=c1.re*c2.im+c2.re*c1.im;
	return c;
}

/*快速傅立葉變換
TD為時域值，F(xiàn)D為頻域值，power為2的冪數(shù)*/
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;
	/*分配運算所需存儲器*/
	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);
	/*蝶形運算*/
	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);
}

/*快速傅立葉反變換，利用快速傅立葉變換
FD為頻域值，TD為時域值，power為2的冪數(shù)*/
void IFFT(COMPLEX *FD, COMPLEX *TD, int power)
{
	int i, count;
	COMPLEX *x;

	/*計算傅立葉反變換點數(shù)*/
	count=1<<power;
	/*分配運算所需存儲器*/
	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(x,TD,power);
	/*求時域點的共軛*/
	for(i=0;i<count;i++)
	{
		TD[i].re/=count;
		TD[i].im=-TD[i].im/count;
	}
	/*釋放存儲器*/
	free(x);
}