EEPW论坛

/**
函数名：void r1fft(x,y,n)
Funtion：实序列快速傅里叶变换（二）
IntPut： x     : 双精度实型一组数据，长度为n，开始时存放要变换数据的实数据，最后存放变换结果的n/2+1的值
                 其存储顺序为[Re(0),Re(1).....Re(n/2),Im(n/2-1)....Im(1)]。其中Re(0)=X(0),Re(n/2)=X(n/2)
                 根据X(k)的共轭对称性，很容易写出后部分的值。
         n     ：整型变量，数据长度，必须是2的证书次幂，n=2m?
**/

void r1fft(double x[],int n)
{
    int i,n1;
    double a,c,e,s,fr,fi,gr,gi,*f,*g;
    f=malloc(n/2*sizeof(double));
    g=malloc(n/2*sizeof(double));
    n1=n/2;
    e=3.141592653589793/n1;
    for(i=0; i<n1; i++)
      {
      f[i]=x[2*i];
      g[i]=x[2*i+1];
      }
    fft(f,g,n1,1);
    x[0]=f[0]+g[0];
    x[n1]=f[0]-g[0];
    for(i=1; i<n1; i++)
      {
      fr=(f[i]+f[n1-i])/2;
      fi=(g[i]-g[n1-i])/2;
      gr=(g[i]+g[n1-i])/2;
      gi=(f[n1-i]-f[i])/2;
      a=i*e;
      c=cos(e);
      s=sin(a);
      x[i]=fr+c*gr+s*gi;
      x[n-i]=fi+c*gi-s*gr;
      }
    free(f);
    free(g);

}