/* Dev-C++ 4.9.9.2 下运行通过 2007.11.2 */// 改进版，(把 1~ Num)  -->  (0~ Num-1)/*****************fft programe*********************/#include<math.h>#include<stdio.h>#include<stdlib.h>struct complex {    double x;    double y;};#define STAGE  4#define Num   (1<<STAGE) double result[Num];struct complex s[Num];const double PI = 3.14159265358979323846; /* b1 ,b2 相乘的结果 b3 */struct complex EE(struct complex b1,struct complex b2){    struct complex b3 ;    b3.x=b1.x*b2.x-b1.y*b2.y ;    b3.y=b1.x*b2.y+b1.y*b2.x ;    return(b3);}/* M_PI=3.14159265358979323846 为math.h 中定义的 pi 值 */     void FFT(struct complex *xin,int N){    int f,m,i,k,j,l ;        struct complex v,w,t ;        /*变址运算*/    for(j=0,i=0;i<N-1;i++)    {        if(i<j)        {            t=xin[j]; xin[j]=xin[i]; xin[i]=t ;            // printf("%d<-->%d\n",i,j);        }        k=N/2 ;        while(k<(j+1)&&(j>0))        {            j=j-k ; k=k/2 ;        }        j=j+k ;    }        {        int le,lei,ip ;        for(l=0;l<STAGE;l++)        {            lei = 1<<l; /* 这里用的是L而不是1 */            le = lei<<1 ;            v.x=1.0 ; v.y=0.0 ;            w.x=cos(PI/lei);            w.y=-sin(PI/lei);            for(j=0;j<lei;j++)            {                for(i=j;i<N;i=i+le)                {                    ip=i+lei ;                    t=EE(xin[ip],v); // DIT                    xin[ip].x=xin[i].x-t.x ;                    xin[ip].y=xin[i].y-t.y ;                    xin[i].x=xin[i].x+t.x ;                    xin[i].y=xin[i].y+t.y ;                }                v=EE(v,w);            }        }    }    return ;}/***************** prntft programe *****************/void prntft(struct complex s[],const int n){    int i;    double result[n+1];    printf("n     real(n)      imag(n)      fabs(n)\n") ;    for(i=0;i<n;i++)     /* 存放的位置为 1~~N  ?? */    {        result[i]=sqrt( s[i].x*s[i].x+s[i].y*s[i].y );        printf("%d     %f     %f     %f\n",i,s[i].x,s[i].y,result[i]) ;    }}/*****************main programe********************/main(){    int i;    /*******赋初值**********/    for(i=0;i<Num;i++)    {        s[i].x=i;        s[i].y=0 ;    }    /***********************/        FFT(s,Num);        prntft(s,Num);    getchar();}

评论