#include<stdio.h>
#include<math.h>
//#include<IFFTC.h>
/**************************************/
/*              定義                  */
/**************************************/
#define N 8                         //fft 點(diǎn)數(shù)
#define pi 3.141592653589793

typedef struct{
    float real;
    float image;
}Complex;

Complex ComMul(Complex p1,Complex p2);
Complex ComAdd(Complex p1,Complex p2);
Complex ComSub(Complex p1,Complex p2);

void FFTcompute(Complex data_in[], Complex data_out[], int len);
void ResetSrc(Complex data_src[], int len, int size);
int BitReverse(int src, int size);
int	ComputeGradeNum(int len);
void CalcW(Complex *W,int len);
void Real2Complex(float data_real[], Complex data_plex[], int len);


//******************************思路說(shuō)明*********************************/
//GradeNum--蝶形運(yùn)算級(jí)數(shù); 
//CellNum--每一級(jí)子fft個(gè)數(shù)；
//Subfftlen--子fft長(zhǎng)度；
//在fft蝶形運(yùn)算中，其運(yùn)算級(jí)數(shù)為log2(N)；
//自左向右每一級(jí)里又有CellNum個(gè)子fft運(yùn)算(CellNum=2^(k-h),基2自左向右減少)；
//每一級(jí)的子fft運(yùn)算長(zhǎng)度為Subfftlen(Subfftlen=2^(h)),基2自左向右增加,根據(jù)fft兩個(gè)source num經(jīng)不同的運(yùn)算得到兩個(gè)目標(biāo)值)；
//CellNum*Subfftlen=N
//----------------------------------------
//1)利用for循環(huán)自左向右計(jì)算每一級(jí)每個(gè)fft
//2)計(jì)算前需對(duì)輸入序列重排
//3)蝶形運(yùn)算為原位運(yùn)算，每次算出來(lái)的結(jié)果直接放入原寄存器
//4)每一級(jí)的比例因子為W(j),       =W(j*CellNum);先把所用到的W值(N/2個(gè))都計(jì)算出來(lái)，用到的時(shí)候直接取 
//                    W(Subfftlen) W(N)
/***********************************************************************/
main()
{

	float data_in[N];
	Complex data_out[N];
	Complex data_in_plex[N];
	int i;
	FILE* fp;


	 if((fp=fopen("C:\MATLAB6p5\work\input","r"))==NULL)
    {
        printf("cannot open file");
  //      exit(0);
    }

    for(i=0;i<N;i++)
        scanf("%f",&data_in[i]);

	Real2Complex(data_in, data_in_plex, N);

	FFTcompute(data_out, data_in_plex, N);
	
//-----------output data---------------------------//
	 for(i=0;i<N;i++)
        printf("%f, %f, \n",data_out[i].real, data_out[i].image);
}

/*****************************************************************/
/*                           子函數(shù)                              */
/*****************************************************************/
/*--------FFT core 算法--------*/
/*In:data_in[]--數(shù)據(jù)列
             N--數(shù)據(jù)列(FFT)長(zhǎng)度              
/*Out:data_out[]--數(shù)據(jù)列after fft
/*-----------------------------*/
void FFTcompute(Complex data_out[], Complex data_in[], int len)
{
	
	Complex W[N/2];                   //比例因子W
	Complex Wmul;
	Complex temp;

	int tmp;
	int	index1, index2;
	int i, j, h;
	int CellNum, Subfftlen, GradeNum;

//-----------Preparation-------------------------//
	CalcW(W, len);								    //計(jì)算所用到的W值(N/2個(gè))
	GradeNum = ComputeGradeNum(len);              //Compute 級(jí)數(shù)
//-----------輸入序列重排(index比特反轉(zhuǎn))--------//
	ResetSrc(data_in, len, GradeNum);
	
//--------- fft computation----------------------//
	for(h=1; h<=GradeNum; h++)
	{
		CellNum = pow(2, (GradeNum-h));                                  
		Subfftlen = pow(2, h);

		for(i=0; i<CellNum; i++)
			for(j=0; j<Subfftlen/2; j++)									//做一次運(yùn)算占用兩個(gè)源數(shù) 
			{

				index1 = j+i*Subfftlen;										// fft數(shù)1index；參見(jiàn)蝶形運(yùn)算圖理解index
				index2 = j+(Subfftlen/2)+i*Subfftlen;						// fft數(shù)2index
				
				Wmul = ComMul(W[CellNum*j], data_in[index2]);                // W()
//				data_in[index1] = ComAdd(data_in[index1], Wmul);
				temp = ComAdd(data_in[index1], Wmul);
				data_in[index2] = ComSub(data_in[index1], Wmul);
				data_in[index1] = temp;
			}
	}


	for(tmp=0; tmp<N; tmp++)
	{
		data_out[tmp] = data_in[tmp];
	}


}
/*--------FFT級(jí)數(shù)--------*/
/*In:len--FFT長(zhǎng)度              
/*Out:M--級(jí)數(shù)
/*-----------------------*/
int	ComputeGradeNum(int len)
{
	int M=0;
	M = int((log10(len))/(log10(2)));
	return M;
}
/*--------W比例因子計(jì)算--------*/
/*In:len--FFT長(zhǎng)度              
/*Out:W--比例因子矩陣
/*-----------------------------*/
void CalcW(Complex *W,int len)
{
	int size = len/2;
	for(int tmp=0; tmp<size; tmp++)
	{
		W[tmp].real = cos(pi*tmp/size);
		W[tmp].image = -1*sin(pi*tmp/size);
	}
	
}
/*--------W比例因子計(jì)算--------*/
/*In:data_src--數(shù)據(jù)列;
          len--數(shù)據(jù)列長(zhǎng)度
		 size--index比特寬度	
/*Out:data_src--重排后數(shù)據(jù)列
/*-----------------------------*/
void ResetSrc(Complex data_src[], int len, int size)
{
	for(int tmp=0; tmp<N; tmp++)
	{
		Complex data_temp;
		int index_reverse = 0;

		index_reverse = BitReverse(tmp, size);
		if(index_reverse>tmp)
		{
			data_temp = data_src[index_reverse];
			data_src[index_reverse] = data_src[tmp];
			data_src[tmp] = data_temp;
		}

	}

}
/*-------------(Num比特反轉(zhuǎn))-----------------*/
/*In:src--source num
/*   size--數(shù)的bit寬度               
/*Out:des--destination num
/*-------------------------------------------*/
int BitReverse(int src, int size)
{
	int temp = src;
	int des = 0;
	int tmp = size-1;
	for(tmp;tmp>=0;tmp--)                  //利用比特移位
	{
		des = ((temp&0x1)<<tmp)|des;
		temp = temp>>1;
	}
	return des;
}
/*--------復(fù)數(shù)計(jì)算(乘，加，減)---------------*/
/*In:Complex p1,p2--source num
/*Out:Res--Computation Result(Complex)
/*-------------------------------------------*/
Complex ComMul(Complex p1,Complex p2)
{
	Complex	Res;
	Res.real = p1.real*p2.real-p1.image*p2.image;
	Res.image = p1.real*p2.image+p1.image*p2.real;
	return Res;
}

Complex ComAdd(Complex p1,Complex p2)
{
    Complex Res;
    Res.real = p1.real + p2.real;
    Res.image = p1.image + p2.image;
    
    return Res;
}

Complex ComSub(Complex p1,Complex p2)
{
    Complex Res;
    Res.real = p1.real - p2.real;
    Res.image = p1.image - p2.image;
    
    return Res;
}
/*--------------實(shí)數(shù)轉(zhuǎn)復(fù)數(shù)---------------*/
/*In:data_real[]--source real num
             len--length of data stream
/*Out:data_plex--destination complex num
/*-------------------------------------------*/
void Real2Complex(float data_real[], Complex data_plex[], int len)
{
	for(int tmp=0; tmp<len; tmp++)
	{
		data_plex[tmp].image = 0;
		data_plex[tmp].real = data_real[tmp];
	
	}
	
}