fft.txt

fft 变换C程序 不错的

/*******************************************************************************
                           下面的算法是FFT算法
********************************************************************************/

int **gFFTBitTable = NULL;
const int MaxFastBits = 16;  //2^16=65536=64K

bool IsPowerOfTwo (int x)//判断X是否是2的N次方的函数,绝对高效的代码
{
  if ( x < 2 )
 return false;
   if ( x & (x-1) )  /* 判断X是否是2的N次方,高效! */
 return false;
   return true;
}

int NumberOfBitsNeeded (int PowerOfTwo)
/*   已知PowerOfTwo是2的N次方,该函数就是求n
 *   如果是输入参数是4,则返回2;如果是8,返回3,如果是16,返回4
 */
{
  int i;
  if ( PowerOfTwo < 2 ) {
 fprintf(stderr, "Error: FFT called with size %d\n", PowerOfTwo);
 exit(1);
  }
  for ( i=0; ; i++ )
 if ( PowerOfTwo & (1 << i) )
   return i;
}

int ReverseBits ( int index, int NumBits )
/* 位反转,比如:
 *             输入,index=0x80,NumBits=16,则执行后返回值为0x01
 *             输入,index=0x82,NumBits=16,则执行后返回值为0x41
 */
{
  int i, rev;
  for ( i=rev=0; i < NumBits; i++ ) {
 rev = (rev << 1) | (index & 1);
 index >>= 1;
  }
  return rev;
}

void InitFFT()
{
  gFFTBitTable = new int *[MaxFastBits];
  int len=2;
  for(int b=1; b<=MaxFastBits; b++) {
 gFFTBitTable[b-1] = new int[len];
        for(int i=0; i<len; i++)
   gFFTBitTable[b-1][i] = ReverseBits(i, b);
 len <<= 1;
  }
}

inline int FastReverseBits(int i, int NumBits)
{
  if (NumBits <= MaxFastBits)
   return gFFTBitTable[NumBits-1][i];
  else
 return ReverseBits(i, NumBits);
}

/*
 * Complex Fast Fourier Transform
 * 复数快速傅里叶变换
 */
void FFT2 (
    int       NumSamples,
    int      InverseTransform,    //逆变换
    float    *RealIn,
    float    *ImagIn,
    float    *RealOut,
    float    *ImagOut )
{
  int NumBits;    /* 保存复数需要的位数 */
  int i, j, k, n;
  int BlockSize, BlockEnd;

  double angle_numerator = 2.0 * M_PI;
  float tr, ti;     /* 临时变量实部,虚部 */

  if ( !IsPowerOfTwo(NumSamples) ) {
 fprintf (stderr, "%d 不是2的n次方\n", NumSamples);
 exit(1);
  }

  if (!gFFTBitTable)
 InitFFT();

  if ( InverseTransform )
 angle_numerator = -angle_numerator;

  NumBits = NumberOfBitsNeeded ( NumSamples );

  /*
  **   Do simultaneous data copy and bit-reversal ordering into outputs...
  */
  if (ImagIn==NULL)
     {
          for ( i=0; i < NumSamples; i++ ) {
         j = FastReverseBits ( i, NumBits );
         RealOut[j] = RealIn[i];
         ImagOut[j] = 0;
          }
     }
     else
     {
          for ( i=0; i < NumSamples; i++ ) {
         j = FastReverseBits ( i, NumBits );
         RealOut[j] = RealIn[i];
         ImagOut[j] = ImagIn[i];
          }
     }

  /*
  **   Do the FFT itself...
  */

  BlockEnd = 1;
  for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 ) {

   double delta_angle = angle_numerator / (double)BlockSize;

   float sm2 = sin ( -delta_angle-delta_angle );
   float sm1 = sin ( -delta_angle );
   float cm2 = cos ( -delta_angle-delta_angle );
   float cm1 = cos ( -delta_angle );
   float w = cm1+cm1;
   float ar0, ar1, ar2, ai0, ai1, ai2;

   for ( i=0; i < NumSamples; i += BlockSize ) {
   ar2 = cm2;
   ar1 = cm1;

   ai2 = sm2;
   ai1 = sm1;

   for ( j=i, n=0; n < BlockEnd; j++, n++ ) {
     ar0 = w*ar1 - ar2;
     ar2 = ar1;
     ar1 = ar0;

     ai0 = w*ai1 - ai2;
     ai2 = ai1;
     ai1 = ai0;

     k = j + BlockEnd;
     tr = ar0*RealOut[k] - ai0*ImagOut[k];
     ti = ar0*ImagOut[k] + ai0*RealOut[k];

     RealOut[k] = RealOut[j] - tr;
     ImagOut[k] = ImagOut[j] - ti;

     RealOut[j] += tr;
     ImagOut[j] += ti;
   }
   }

   BlockEnd = BlockSize;
  }

  /*
  **   Need to normalize if inverse transform...
  */

  if ( InverseTransform )
  {
 float denom = (float)NumSamples;
 for ( i=0; i < NumSamples; i++ ) {
   RealOut[i] /= denom;
   ImagOut[i] /= denom;
 }
  }
}

/*
 * Real Fast Fourier Transform
 *
 * This function was based on the code in Numerical Recipes in C.
 * In Num. Rec., the inner loop is based on a single 1-based array
 * of interleaved real and imaginary numbers.  Because we have two
 * separate zero-based arrays, our indices are quite different.
 * Here is the correspondence between Num. Rec. indices and our indices:
 *
 * i1  <->  real[i]
 * i2  <->  imag[i]
 * i3  <->  real[n/2-i]
 * i4  <->  imag[n/2-i]
 */

void RealFFT(int NumSamples,float *RealIn,float *RealOut,float *ImagOut)
{
  int Half = NumSamples>>1;
  int i;

  float theta = M_PI / float(Half);

  float *tmpReal = new float[Half];
  float *tmpImag = new float[Half];

  for(i=0; i<Half; i++) {
 tmpReal[i] = RealIn[ i<<1 ];
 tmpImag[i] = RealIn[ (i<<1) +1 ];
  }
  FFT2(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);
  float wtemp = float(sin(0.5*theta));
  float wpr = -2.0 * wtemp * wtemp;
  float wpi = float(sin(theta));
  float wr = 1.0+wpr;
  float wi = wpi;
  int i3;
  float h1r, h1i, h2r, h2i;
  for(i=1; i<Half/2; i++) {
 i3 = Half-i;
 h1r = 0.5*(RealOut[i]+RealOut[i3]);
 h1i = 0.5*(ImagOut[i]-ImagOut[i3]);
 h2r = 0.5*(ImagOut[i]+ImagOut[i3]);
 h2i = -0.5*(RealOut[i]-RealOut[i3]);
 RealOut[i]=h1r+wr*h2r-wi*h2i;
 ImagOut[i]=h1i+wr*h2i+wi*h2r;
 RealOut[i3]=h1r-wr*h2r+wi*h2i;
 ImagOut[i3]=-h1i+wr*h2i+wi*h2r;
 wr=(wtemp=wr)*wpr-wi*wpi+wr;
 wi=wi*wpr+wtemp*wpi+wi;
  }
  RealOut[0] = (h1r=RealOut[0])+ImagOut[0];
  ImagOut[0] = h1r-ImagOut[0];
  delete []tmpReal;
  delete []tmpImag;
}

/*
 * PowerSpectrum
 * 功率谱
 *
 * This function computes the same as RealFFT, above, but
 * adds the squares of the real and imaginary part of each
 * coefficient, extracting the power and throwing away the
 * phase.
 *
 * For speed, it does not call RealFFT, but duplicates some
 * of its code.
 * 为了速度,这里没有调用RealFFT,但是复制了他的一些代码
 */

void PowerSpectrum(
    int       NumSamples,
    float    *In,
    float    *Out )
{
  int Half = NumSamples/2;
  int i;
  float theta = M_PI / Half;
  float *tmpReal = new float[Half];
  float *tmpImag = new float[Half];
  float *RealOut = new float[Half];
  float *ImagOut = new float[Half];
  for(i=0; i<Half; i++) {
 tmpReal[i] = In[i<<1];
 tmpImag[i] = In[(i<<1) +1];
  }
  FFT2(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);
  float wtemp = float(sin(0.5*theta));
  float wpr = -2.0 * wtemp * wtemp;
  float wpi = float(sin(theta));
  float wr = 1.0+wpr;
  float wi = wpi;
  int i3;
  float h1r, h1i, h2r, h2i, rt, it;
  for(i=1; i<Half>>1; i++) {
 i3 = Half-i;
 h1r = 0.5*(RealOut[i]+RealOut[i3]);
 h1i = 0.5*(ImagOut[i]-ImagOut[i3]);
 h2r = 0.5*(ImagOut[i]+ImagOut[i3]);
 h2i = -0.5*(RealOut[i]-RealOut[i3]);
 rt = h1r+wr*h2r-wi*h2i;
 it = h1i+wr*h2i+wi*h2r;
 Out[i] = rt*rt + it*it;
 rt = h1r-wr*h2r+wi*h2i;
 it = -h1i+wr*h2i+wi*h2r;
 Out[i3] = rt*rt + it*it;
 wr=(wtemp=wr)*wpr-wi*wpi+wr;
 wi=wi*wpr+wtemp*wpi+wi;
  }

  rt = (h1r=RealOut[0])+ImagOut[0];
  it = h1r-ImagOut[0];
  Out[0] = rt*rt + it*it;

  rt = RealOut[Half>>1];
  it = ImagOut[Half>>1];
  Out[Half>>1] = rt*rt + it*it;

  delete[] tmpReal;
  delete[] tmpImag;
  delete[] RealOut;
  delete[] ImagOut;
}