上一个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);  
  }