上一個NB人寫的 傅立葉變換和反變換

 
  //////////////////////////////////////////////////////////  
  //   internal   definitions  
   
  #define   PI   (double)3.14159265359  
   
  /*complex   number*/  
  typedef   struct  
  {  
  double   re;  
  double   im;  
  }COMPLEX;  
   
  /*complex   add*/  
  COMPLEX   Add(COMPLEX   c1,   COMPLEX   c2)  
  {  
  COMPLEX   c;  
  c.re=c1.re+c2.re;  
  c.im=c1.im+c2.im;  
  return   c;  
  }  
   
  /*complex   substract*/  
  COMPLEX   Sub(COMPLEX   c1,   COMPLEX   c2)  
  {  
  COMPLEX   c;  
  c.re=c1.re-c2.re;  
  c.im=c1.im-c2.im;  
  return   c;  
  }  
   
  /*complex   multiple*/  
  COMPLEX   Mul(COMPLEX   c1,   COMPLEX   c2)  
  {  
  COMPLEX   c;  
  c.re=c1.re*c2.re-c1.im*c2.im;  
  c.im=c1.re*c2.im+c2.re*c1.im;  
  return   c;  
  }  
  //////////////////////////////////////////////////////////  
   
   
   
  這是FFT/****************************************************  
  FFT()  
   
  參數：  
   
  TD為時域值  
  FD為頻域值  
  power為2的冪數  
   
  傳回值：  
   
  無  
   
  說明：  
   
  本函數實現快速傅立葉變換  
  ****************************************************/  
  void   FFT(COMPLEX   *   TD,   COMPLEX   *   FD,   int   power)  
  {  
  int   count;  
  int   i,j,k,bfsize,p;  
  double   angle;  
  COMPLEX   *W,*X1,*X2,*X;  
   
  /*計算傅立葉變換點數*/  
  count=1<<power;  
   
  /*分配運算所需儲存空間*/  
  W=(COMPLEX   *)malloc(sizeof(COMPLEX)*count/2);  
  X1=(COMPLEX   *)malloc(sizeof(COMPLEX)*count);  
  X2=(COMPLEX   *)malloc(sizeof(COMPLEX)*count);  
   
  /*計算加權係數*/  
  for(i=0;i<count/2;i++)  
  {  
  angle=-i*PI*2/count;  
  W[i].re=cos(angle);  
  W[i].im=sin(angle);  
  }  
   
  /*將時域點寫入儲存空間*/  
  memcpy(X1,TD,sizeof(COMPLEX)*count);  
   
  /*蝶形運算*/  
  for(k=0;k<power;k++)  
  {  
  for(j=0;j<1<<k;j++)  
  {  
  bfsize=1<<(power-k);  
  for(i=0;i<bfsize/2;i++)  
  {  
  p=j*bfsize;  
  X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);  
  X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);  
  }  
  }  
  X=X1;  
  X1=X2;  
  X2=X;  
  }  
   
  /*重新排序*/  
  for(j=0;j<count;j++)  
  {  
  p=0;  
  for(i=0;i<power;i++)  
  {  
  if   (j&(1<<i))   p+=1<<(power-i-1);  
  }  
  FD[j]=X1[p];  
  }  
   
  /*釋放儲存空間*/  
  free(W);  
  free(X1);  
  free(X2);  
  }  
   
   
  /****************************************************  
  IFFT()  
   
  參數：  
   
  FD為頻域值  
  TD為時域值  
  power為2的冪數  
   
  傳回值：  
   
  無  
   
  說明：  
   
  本函數利用快速傅立葉變換實現傅立葉反變換  
  ****************************************************/  
  void   IFFT(COMPLEX   *   FD,   COMPLEX   *   TD,   int   power)  
  {  
  int   i,   count;  
  COMPLEX   *x;  
   
  /*計算傅立葉反變換點數*/  
  count=1<<power;  
   
  /*分配運算所需儲存空間*/  
  x=(COMPLEX   *)malloc(sizeof(COMPLEX)*count);  
   
  /*將頻域點寫入儲存空間*/  
  memcpy(x,FD,sizeof(COMPLEX)*count);  
   
  /*求頻域點的共軛*/  
  for(i=0;i<count;i++)  
  x[i].im   =   -x[i].im;  
   
  /*調用FFT*/  
  FFT(x,   TD,   power);  
   
  /*求時域點的共軛*/  
  for(i=0;i<count;i++)  
  {  
  TD[i].re   /=   count;  
  TD[i].im   =   -TD[i].im   /   count;  
  }  
   
  /*釋放儲存空間*/  
  free(x);  
  }