/*************************************************************************   
    *   
    *   函數名稱：   
    *       FFT()   
    *   
    *   參數:   
    *       complex<double>   *   TD -   指向時域數組的指針   
    *       complex<double>   *   FD -   指向頻域數組的指針   
    *       r －2的冪數，即迭代次數   
    *   
    *   返回值:   
    *       無。   
    *   
    *   說明:   
    *       該函數用來實現快速付立葉變換。   
    *   
    ************************************************************************/   
  VOID   WINAPI   FFT(complex<double>   *   TD,   complex<double>   *   FD,   int   r)   
  {   
  //   付立葉變換點數   
  LONG count;   
    
  //   循環變量   
  int i,j,k;   
    
  //   中間變量   
  int bfsize,p;   
    
  //   角度   
  double angle;   
    
  complex<double>   *W,*X1,*X2,*X;   
    
  //   計算付立葉變換點數   
  count   =   1   <<   r;   
    
  //   分配運算所需存儲器   
  W     =   new   complex<double>[count   /   2];   
  X1   =   new   complex<double>[count];   
  X2   =   new   complex<double>[count];   
    
  //   計算加權系數   
  for(i   =   0;   i   <   count   /   2;   i++)   
  {   
  angle   =   -i   *   PI   *   2   /   count;   
  W[i]   =   complex<double>   (cos(angle),   sin(angle));   
  }   
    
  //   將時域點寫入X1   
  memcpy(X1,   TD,   sizeof(complex<double>)   *   count);   
    
  //   采用蝶形算法進行快速付立葉變換   
  for(k   =   0;   k   <   r;   k++)   
  {   
  for(j   =   0;   j   <   1   <<   k;   j++)   
  {   
  bfsize   =   1   <<   (r-k);   
  for(i   =   0;   i   <   bfsize   /   2;   i++)   
  {   
  p   =   j   *   bfsize;   
  X2[i   +   p]   =   X1[i   +   p]   +   X1[i   +   p   +   bfsize   /   2];   
  X2[i   +   p   +   bfsize   /   2]   =   (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   <   r;   i++)   
  {   
  if   (j&(1<<i))   
  {   
  p+=1<<(r-i-1);   
  }   
  }   
  FD[j]=X1[p];   
  }   
    
  //   釋放內存   
  delete   W;   
  delete   X1;   
  delete   X2;   
  }   
  /*************************************************************************   
    *   
    *   函數名稱：   
    *       Fourier()   
    *   
    *   參數:   
    *       LPSTR   lpDIBBits         -   指向源DIB圖像指針   
    *       LONG     lWidth               -   源圖像寬度（象素數）   
    *       LONG     lHeight             -   源圖像高度（象素數）   
    *   
    *   返回值:   
    *       BOOL                               -   成功返回TRUE，否則返回FALSE。   
    *   
    *   說明:   
    *       該函數用來對圖像進行付立葉變換。   
    *   
    ************************************************************************/   
  BOOL   WINAPI   Fourier(LPSTR   lpDIBBits,   LONG   lWidth,   LONG   lHeight)   
  {   
    
  //   指向源圖像的指針   
  unsigned   char* lpSrc;   
    
  //   中間變量   
  double dTemp;   
    
  //   循環變量   
  LONG i;   
  LONG j;   
    
  //   進行付立葉變換的寬度和高度（2的整數次方）   
  LONG w;   
  LONG h;   
    
  int wp;   
  int hp;   
    
  //   圖像每行的字節數   
  LONG lLineBytes;   
    
  //   計算圖像每行的字節數   
  lLineBytes   =   WIDTHBYTES(lWidth   *   8);   
    
  //   賦初值   
  w   =   1;   
  h   =   1;   
  wp   =   0;   
  hp   =   0;   
    
  //   計算進行付立葉變換的寬度和高度（2的整數次方）   
  while(w   *   2   <=   lWidth)   
  {   
  w   *=   2;   
  wp++;   
  }   
    
  while(h   *   2   <=   lHeight)   
  {   
  h   *=   2;   
  hp++;   
  }   
    
  //   分配內存   
  complex<double>   *TD   =   new   complex<double>[w   *   h];   
  complex<double>   *FD   =   new   complex<double>[w   *   h];   
    
  //   行   
  for(i   =   0;   i   <   h;   i++)   
  {   
  //   列   
  for(j   =   0;   j   <   w;   j++)   
  {   
  //   指向DIB第i行，第j個象素的指針   
  lpSrc   =   (unsigned   char*)lpDIBBits   +   lLineBytes   *   (lHeight   -   1   -   i)   +   j;   
    
  //   給時域賦值   
  TD[j   +   w   *   i]   =   complex<double>(*(lpSrc),   0);   
  }   
  }   
    
  for(i   =   0;   i   <   h;   i++)   
  {   
  //   對y方向進行快速付立葉變換   
  FFT(&TD[w   *   i],   &FD[w   *   i],   wp);   
  }   
    
  //   保存變換結果   
  for(i   =   0;   i   <   h;   i++)   
  {   
  for(j   =   0;   j   <   w;   j++)   
  {   
  TD[i   +   h   *   j]   =   FD[j   +   w   *   i];   
  }   
  }   
    
  for(i   =   0;   i   <   w;   i++)   
  {   
  //   對x方向進行快速付立葉變換   
  FFT(&TD[i   *   h],   &FD[i   *   h],   hp);   
  }   
    
  //   行   
  for(i   =   0;   i   <   h;   i++)   
  {   
  //   列   
  for(j   =   0;   j   <   w;   j++)   
  {   
  //   計算頻譜   
  dTemp   =   sqrt(FD[j   *   h   +   i].real()   *   FD[j   *   h   +   i].real()   +     
                    FD[j   *   h   +   i].imag()   *   FD[j   *   h   +   i].imag())   /   100;   
    
  //   判斷是否超過255   
  if   (dTemp   >   255)   
  {   
  //   對于超過的，直接設置為255   
  dTemp   =   255;   
  }   
    
  //   指向DIB第(i<h/2   ?   i+h/2   :   i-h/2)行，第(j<w/2   ?   j+w/2   :   j-w/2)個象素的指針   
  //   此處不直接取i和j，是為了將變換后的原點移到中心   
  //lpSrc   =   (unsigned   char*)lpDIBBits   +   lLineBytes   *   (lHeight   -   1   -   i)   +   j;   
  lpSrc   =   (unsigned   char*)lpDIBBits   +   lLineBytes   *     
  (lHeight   -   1   -   (i<h/2   ?   i+h/2   :   i-h/2))   +   (j<w/2   ?   j+w/2   :   j-w/2);   
    
  //   更新源圖像   
  *   (lpSrc)   =   (BYTE)(dTemp);   
  }   
  }   
    
  //   刪除臨時變量   
  delete   TD;   
  delete   FD;   
    
  //   返回   
  return   TRUE;   
  }