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

#define pi  3.14159265359

/*復數定義*/
typedef struct
{
	double re;
	double im;
}COMPLEX;

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

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

/*復數乘運算*/
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;
}

//DIF_FFT
//TD為時域值，FD為頻域值，power為2的冪數
void FFT(COMPLEX *TD, COMPLEX *FD, int power)
{
	int count;
	int i,j,k,bfsize,p;
	double angle;
	COMPLEX *W, *X1, *X2, *X;

	/*計算傅立葉變換點數*/
	count=1<<power;
	/*分配運算所需存儲器*/
	W =(COMPLEX *)malloc(sizeof(COMPLEX) * count/2);
	X1=(COMPLEX *)malloc(sizeof(COMPLEX) * count);
	X2=(COMPLEX *)malloc(sizeof(COMPLEX) * count);
	/*計算加權系數*/
	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的冪數*/
void IFFT(COMPLEX *FD, COMPLEX *TD, int power)
{
	int i, count;
	COMPLEX *x;

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

//分裂基FFT的遞歸算法
void L_Atom(double* real,double* imag,int size)
{
    int i,j,half,half2;
    double xa,xb,xc;

    half  = size/2;
    half2 = half/2;
    for(i = 0;i < half;i++)
    {
        xa = real[i];
        xb = imag[i];
        real[i] += real[half + i];
        imag[i] += imag[half + i];
        real[half + i] = xa - real[half + i];
        imag[half + i] = xb - imag[half + i];
    }
    if(size > 2)
    {
        for(i = 0,j = half;i < half2;i++,j++)
        {
            xa = real[j];
            xb = imag[j];
            
            real[j] += imag[j + half2];
            imag[j] -= real[j + half2];

            xc = real[j];
            real[j] = xc * cos(i * pi/half) + imag[j] * sin(i * pi/half);
            imag[j] = imag[j] * cos(i * pi/half) - xc * sin(i * pi/half);
            
            xc = real[j + half2];
            real[j + half2] = xa - imag[j + half2];
            imag[j + half2] = xb + xc;
            
            xc = real[j + half2];
            real[j + half2] = xc * cos(3 * i * pi/half) + imag[j + half2] * sin(3 * i * pi/half);
            imag[j + half2] = imag[j + half2] * cos(3 * i * pi/half) - xc * sin(3 * i * pi/half);
        }
        if(size > 4)
        {
            L_Atom(real + half,imag + half,size/4);
            L_Atom(real + half + half2,imag + half + half2,size/4);
        }
        L_Atom(real,imag,size/2);
    }
}