/*
FFT function
Write By WangGuanjin at 2005.12.05 16:47
first time Modified By me at 2005.12.06
*/
#include <math.h>
#include <stdio.h>
#define pai 3.1415926

int num(int n,int N)//輸出第 n 個 位置的序號
{
 int i,s;
 int p = int( log10(N) / log10(2) + 0.01 ) ;
 int *a = new int[p];
 for(i=0;i<p;i++) a[i] = 0;
 for(i=0;i<p;i++)
 {
  a[i] = n % 2;
  if( n < 2 ) break;
  n = n / 2;
 }
 a[i] = n;
 s = 0;
 for(i=0;i<p;i++)
  s += a[i]*int( pow(2,p-i-1) );
 return s;

}

void W(double *a,double *b,int N,int k,int flag)
//a,b為返回值實部和虛部
//flag==1 執(zhí)行傅立葉變化，flag==-1,反變換
{

 if(1==flag)
 {
  *a =  cos( -2 * pai * k / N ); //返回實部
  *b =  sin( -2 * pai * k / N ); //返回虛部
 }
 else
 {
  *a =  cos(  2 * pai * k / N );
  *b =  sin(  2 * pai * k / N );
 }

}

void FFT(double *Pr,double *Pi,double xr[],double xi[],int N,int flag)
/*
Pr 返回傅立葉變換的實部
Pi 返回傅立葉變換的虛部
xr 輸入序列的實部
xi 輸入序列的虛部
N  序列的長度
flag 1 執(zhí)行傅立葉變化，-1表示反變換
如果輸入序列為實數(shù)序列，那么參數(shù)double xi[]改為NULL
*/
{
 int i,j;

 if(xi==NULL)//輸入為實數(shù)序列
 {

  for(i=0;i<N;i++)
  {
   Pr[i] = xr[ num(i,N) ];
   Pi[i] = 0;
  }
 }
 else//輸入序列為復(fù)數(shù)序列
 {

  int q;
  for(i=0;i<N;i++)
  {
   q = num(i,N);
   Pr[i] = xr[ q ];
   Pi[i] = xi[ q ];
  }

 }
 //予處理已經(jīng)完畢

 //下面進(jìn)行FFT計算

 int m = int ( log10(N) / log10(2) +0.05 );//蝶形運算的層數(shù)
 int t,k;
 double a = 0,b = 0;
 int *id = new int[N];

 double tempr,tempi,temp1,temp2;

 for( i = 0 ; i < m ; i++ )
 {
  for(j = 0 ; j < N ; j++) id[j] = 0;
  t = int( pow( 2,i ) );//結(jié)合的步長 note!
  k = 0;
  for( j = 0 ; j < N ; j++ )
  {
   if( id[j] ) {k = 0;continue;}
   W(&a,&b,2*t,k,flag);
   tempr   = a*Pr[j+t] - b*Pi[j+t];
   tempi   = a*Pi[j+t] + b*Pr[j+t];
   temp1 = Pr[j] ; temp2 = Pi[j];
   Pr[j+t] = temp1 - tempr;  Pi[j+t]  = temp2 - tempi;
   Pr[j]   = temp1 + tempr;    Pi[j]  = temp2 + tempi;
   id[j] = id[j+t] = 1;
   k = k + 1;

  }

 }
 int ttt = 1;
 if( flag!=1 ) ttt = N;
 for(i=0;i<N;i++)
 {
  Pr[i] = Pr[i] / N ;
  Pi[i] = -Pi[i] / N;
 }

}