FFT function

/*
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 执行傅立叶变化，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 执行傅立叶变化，-1表示反变换
如果输入序列为实数序列，那么参数double xi[]改为NULL
*/
{
 int i,j;

 if(xi==NULL)//输入为实数序列
 {

  for(i=0;i<N;i++)
  {
   Pr[i] = xr[ num(i,N) ];
   Pi[i] = 0;
  }
 }
 else//输入序列为复数序列
 {

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

 }
 //予处理已经完毕

 //下面进行FFT计算

 int m = int ( log10(N) / log10(2) +0.05 );//蝶形运算的层数
 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 ) );//结合的步长 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;
 }

}