my_fft.cpp

C++实现的FFT

#include <math.h>
#include "My_ComplexDefines.h"
#include "My_Constants.h"
//   MyFFT(频率域指针引用，信号数组指针引用，总个数，fft or ifft指示)
//                                                  0表示fft,1表示ifft
//time rearrangement
void Rearrange(dComplex*& dcInput, long lSize);
//无返回值，输入时间域数组指针与频率域数组指针与存放中间变量的数组指针
//此外还有数组大小N=2^r
void MyFFT(dComplex*& dcResult, dComplex*& dcInput, long lSize, int iIfftFlag)
{
 long lGroupSize,          //时间抽取法的分组大小的一半
   lBase;               //W的基的次数
 dComplex dcW;             //傅里叶变换中的w
 dComplex dcX;             //中间变量
 if (iIfftFlag == 0) {      //FFT
  //dcW.real(cos(-2.0 * PI / static_cast<double>(lSize)));
  //dcW.imag(sin(-2.0 * PI / static_cast<double>(lSize)));
  dcW = dComplex(cos(-2.0 * PI / static_cast<double>(lSize)), sin(-2.0 * PI / static_cast<double>(lSize)));
 } else {               //IFFT
  //dcW.real(cos(2.0 * PI / static_cast<double>(lSize)));
  //dcW.imag(sin(2.0 * PI / static_cast<double>(lSize)));
  dcW = dComplex(cos(2.0 * PI / static_cast<double>(lSize)), sin(2.0 * PI / static_cast<double>(lSize)));
 }
 long i = lSize, j, lR = 0;
 while (i!=1) {lR++; i/=2;}
 Rearrange(dcInput, lSize);
 for (i=0; i<lSize; i++) {
  dcResult[i] = dcInput[i];
 }
 dComplex *dcTemp = new dComplex [lSize];
 dComplex *dcWarray = new dComplex [lSize];
 for (i=0; i<lSize; i++)
  dcTemp[i] = dcResult[i];
 for (i=0; i<lSize/2; i++) {             //频率数组赋值
  dcWarray[i] = pow(dcW, i);
  dcWarray[i+lSize/2] = -dcWarray[i];
 }
 for (i=1; i<=lR; i++) {
  lGroupSize = 1;
  lBase = 1;
  for (j=0; j<i-1; j++)
   lGroupSize *= 2;
  for (j=0; j<lR-i; j++)
   lBase *= 2;
  for (j=0; j<lSize; j++) {
   if (j%(lGroupSize*2)<lGroupSize) {
    dcX = dcWarray[(lBase*j)%lSize] * dcTemp[j+lGroupSize];
    dcResult[j] = dcTemp[j] + dcX;
   } else {
    dcX = dcWarray[(lBase*j)%lSize] * dcTemp[j];
    dcResult[j] = dcTemp[j-lGroupSize] + dcX;
   }
  }
  for (j=0; j<lSize; j++)
   dcTemp[j] = dcResult[j];
 
 }
 if (iIfftFlag == 1) {           //IFFT要求的系数1/N
  for (i=0; i<lSize; i++) {
   dcResult[i] /= lSize;
  }
 }
 delete []dcTemp;
 delete []dcWarray;
}
void Rearrange(dComplex*& dcInput, long lSize)
{
 long lR = 0;
 long i, j, lTemp;
    long lBase ;  //二进制基
 lTemp = lSize;
 while(lTemp != 1) {
  lR++;
  lTemp /= 2;
 }
 long *lNum = new long [lR];
 long **lIndex;
 lIndex = new long *[lSize];
 for (i=0; i<lSize; i++)
  lIndex[i] = new long [2];
 dComplex *dcTempArray = new dComplex [lSize];
 for (i=0; i<lSize; i++) {
  lIndex[i][0] = i;
  lIndex[i][1] = 0;
  lTemp = i;
  for (j=0; j<lR; j++) {
   lNum[j] = lTemp % 2;
   lTemp /= 2 ;
  }
  lBase = 1;
  for (j=lR-1; j>=0; j--) {
   lIndex[i][1] += lBase * lNum[j];
   lBase *= 2;
  }
 }
 for (i=0; i<lSize-1; i++) {
  for (j=i+1; j<lSize; j++) {
   if (lIndex[i][1] > lIndex[j][1]) {
    lTemp = lIndex[i][1];
    lIndex[i][1] = lIndex[j][1];
    lIndex[j][1] = lTemp;
    lTemp = lIndex[i][0];
    lIndex[i][0] = lIndex[j][0];
    lIndex[j][0] = lTemp;
   }
  }
 }
 for (i=0; i<lSize; i++) {
  dcTempArray[i] = dcInput[lIndex[i][0]];    //时间抽取过程
 }
 for (i=0; i<lSize; i++) {
  dcInput[i] = dcTempArray[i];
 }
 delete []lNum;
 for (i=0; i<lSize; i++)
  delete []lIndex[i];
 delete []lIndex;
 delete []dcTempArray;
}