#include <iom128v.h>
#include <math.h>
#define PI 3.1415926

unsigned char i,w[32];
float dataR[]={5,5,5,5,0,0,0,0};
float dataI[]={0,0,0,0,0,0,0,0};

void FFT(float dataR[],float dataI[])
{
   int x0,x1,x2,xx;
   
   int L,j,k,b,p;
   float TR,TI,temp;

/********** following code invert sequence ************/
 for(i=0;i<8;i++)
  { 
   x0=x1=x2=0;
   x0=i&0x01; 
   x1=(i/2)&0x01; 
   x2=(i/4)&0x01;
   xx=x0*4+x1*2+x2;
   dataI[xx]=dataR[i];   //實部的數放在虛部數組中
  }
 for(i=0;i<8;i++)
  { 
   dataR[i]=dataI[i];    //將虛部保存的數傳回實部
   dataI[i]=0; 
  }
/************** following code FFT *******************/
 for(L = 1 ; L <= 3 ; L ++ )            // for (1)  L 代表第幾級
  { 
    b = 1 ; 
    i = L - 1 ;          
    while ( i > 0 )           //計算第L層旋轉因子的個數，即輸入口間距 b= 2^(L-1)  
     {
       b=b*2; 
       i--;
     }                        
   
   for ( j = 0 ; j <= b - 1 ; j ++ )    // for (2)  j 用來控制計算第幾個旋轉因子
    {
       p = 1 ;
       i = 3 - L ;     
       
      while ( i > 0 )                           /* p=pow(2,3-L)*j; */
         {
          p = p * 2 ; 
          i -- ;
         }
      p=p*j;
       
        for ( k = j ; k < 8 ; k = k + 2 * b )  /* for (3) */
           {  p=2*PI*p/8;
              TR = dataR [ k ] ;
              TI = dataI [ k ] ; 
              temp = dataR [ k + b ] ;
              dataR[k]=dataR[k]+dataR[k+b]*cos(p)+dataI[k+b]*sin(p);
              dataI[k]=dataI[k]-dataR[k+b]*sin(p)+dataI[k+b]*cos(p);
              dataR[k+b]=TR-dataR[k+b]*cos(p)-dataI[k+b]*sin(p);
              dataI[k+b]=TI+temp*sin(p)-dataI[k+b]*cos(p);
			  
			  
			  /*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 ] ;*/
           } /* END for (3) */
    
     } /* END for (2) */
  
  } /* END for (1) */

  for ( i = 0 ; i < 32 ; i ++ )
    {                        /* 只需要32次以下的諧波進行分析 */
      w [ i ] = sqrt ( dataR [ i ] * dataR [ i ] + dataI [ i ] * dataI [ i ] ) ;
        if ( i )
           {
             w [ i ] = w [ i ] / 64 ;
           }
     }
   w [ 0 ] = w [ 0 ] / 2 ;

} 

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