cfft.txt

在c语言环境下实现快速傅立叶变换的功能

[代码]FFT 快速傅立叶变换 (c语言)

[代码]FFT 快速傅立叶变换 (c语言)

算法参见徐世良《计算机常用算法》第二版 第294页

void KFFT(
BNU::vector<std::complex<double> > &x,
BNU::vector<std::complex<double> > &y,
int inverse)
{
  int n = x.size();
  int k = log((double)n)/log(2.);
  int t ,m ,s , i, j, NV, r;
  double angle;
  std::complex<double> V,PODD;
  for (t = 0 ; t <= n-1 ; t++) {
    m = t;
    s = 0;
    for (i = 0 ; i <=k-1 ; i++) {
      j = int(m/2);
      s = 2*s + m - 2*j;
      m = j;
    }                                         
  if (&x!=&y) {
      y[t] = x[s];
    }else {
       if (s>t) {
           V=x[t];
           y[t]=x[s];
           y[s]=V;
       }
    }
  }
  for (t = 0 ; t <= n-2 ; t+=2) {
    V = y[t];
    y[t] = V + y[t+1];
    y[t+1] = V - y[t+1];
  }
  m = n/2;
  NV = 2;
  for (r = k-2 ; r >=0 ; r-- ) {
    m = m/2;
    NV = 2*NV;
    for (t = 0 ; t <= (m-1)*NV ; t=t+NV ){
      for (j = 0 ; j <= (NV/2)-1 ; j++) {             
  angle = -(2.0*M_PI*inverse*j)/NV;
        PODD = complex<double>(cos(angle),sin(angle))*y[t+j+(NV/2)];
        y[t+j+(NV/2)] = y[t+j] - PODD;
        y[t+j] = y[t+j] + PODD;
      }
    }
  };
for (i = 0 ; i < n; i++) {
   if (inverse==-1)
        y[i] = y[i]/double(n);
}
}