快速傅立叶变换

/*************************************************************************
 *
 * 函数名称：
 *   FFT()
 *
 * 参数:
 *   complex * TD - 指向时域数组的指针
 *   complex * FD - 指向频域数组的指针
 *   r      －2的幂数，即迭代次数
 *
 * 返回值:
 *   无。
 *
 * 说明:
 *   该函数用来实现快速付立叶变换。
 *
 ************************************************************************/
VOID WINAPI FFT(complex * TD, complex * FD, int r)
{
 // 付立叶变换点数
 LONG count;
 
 // 循环变量
 int  i,j,k;
 
 // 中间变量
 int  bfsize,p;
 
 // 角度
 double angle;
 
 complex *W,*X1,*X2,*X;
 
 // 计算付立叶变换点数
 count = 1 << r;
 
 // 分配运算所需存储器
 W  = new complex[count / 2];
 X1 = new complex[count];
 X2 = new complex[count];
 
 // 计算加权系数
 for(i = 0; i < count / 2; i++)
 {
  angle = -i * PI * 2 / count;
  W[i] = complex (cos(angle), sin(angle));
 }
 
 // 将时域点写入X1
 memcpy(X1, TD, sizeof(complex) * 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<   }
  }
  X  = X1;
  X1 = X2;
  X2 = X;
 }
 
 // 重新排序
 for(j = 0; j < count; j++)
 {
  p = 0;
  for(i = 0; i < r; i++)
  {
   if (j&(1<   {
    p+=1<<(r-i-1);
   }
  }
  FD[j]=X1[p];
 }
 
 // 释放内存
 delete W;
 delete X1;
 delete X2;
}