#include <math.h>
#define uint unsigned int
#define uchar unsigned char
#define PI 3.1415926

uint x0,x1,x2,x3,x4,x5,x6;
int L,i,j,k,b,p;
uchar w[32];
uint TR,TI,temp;

int dataR[]={5,5,5,5,0,0,0,0,5,5,5,5,0,0,0,0};
int dataI[16];

void CHANGE()
 {                              //i為原始存放位置，最后得invert_pos為倒位序存放位置
  int b0 , b1 , b2 , b3 , b4 , b5 , b6 , all;
      b0 = b1 = b2 = b3 = b4 = b5 = b6 = all = 0;
  b0  = i&0x01; 
  b1  = (i>>1)&0x01; 
  b2  = (i>>2)&0x01;
  b3  = (i>>3)&0x01;
  b4  = (i>>4)&0x01;
  b5  = (i>>5)&0x01;
  b6  = (i>>6)&0x01;
  all = b0*64 + b1*32 + b2*16 + b3*8 + b4*4 + b5*2 + b6;
  dataI[all] = dataR[i];
 }

 
void FFTT()
 {  TR=dataR[k]; 
    TI=dataI[k];
    temp=dataR[k+b]; 
    dataR[k]=dataR[k]+dataR[k+b]*cos(2*PI*p/16)+dataI[k+b]*sin(2*PI*p/16);
    dataI[k]=dataI[k]-dataR[k+b]*sin(2*PI*p/16)+dataI[k+b]*cos(2*PI*p/16);
    dataR[k+b]=TR-dataR[k+b]*cos(2*PI*p/16)-dataI[k+b]*sin(2*PI*p/16);
    dataI[k+b]=TI+temp*sin(2*PI*p/16)-dataI[k+b]*cos(2*PI*p/16);
	
	/*dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];
    dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];
    dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];
    dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p];*/

 }
void FFT(int dataR[],int dataI[])
{
// 以下是數據掉頭功能
  for(i = 0 ; i < 16 ; i++)
   { 
    CHANGE ( ) ;
   }
  for(i=0;i<16;i++)
   {
    dataR[i]=dataI[i]; dataI[i]=0; 
   }
// FFT算法
  for(L=1;L<=4;L++) 						// 第一層循環
  {
     b=1; i=L-1;
     while(i>0) 
    {
      b=b*2; i--;
    } 
       for(j=0;j<=b-1;j++)						 // 第二層循環
         { 
           p=1; i=4-L;
           while(i>0)					
               {
			   p=p*2; i--;
			   }
           p=p*j;
            for(k=j;k<16;k=k+2*b) 				// 第三層循環
              { 
                FFTT();
              }
          }
  } 
  for(i=0;i<32;i++) 							// 32次以下的諧波分析
   {
    w[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);
    w[i]=w[i]/64;
   }
  w[0]=w[0]/2;
}  

void main ()
  {
    FFT(dataR,dataI);
	while(1);
  }