fft.cpp

图像快速傅立叶变换

/*************************************************************************   
    *   
    *   函数名称：   
    *       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;   
  }