/* Dev-C++ 4.9.9.2 下运行通过 2007.11.2 */
//

/*****************fft programe*********************/
#include<math.h>
#include<stdio.h>
#include<stdlib.h>

struct complex {
    double x;
    double y;
};

#define Num 32
double result[Num+1];
struct complex s[Num+1];

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,nv2,nm1,i,k,j=1,l ;
    /*int f,m,nv2,nm1,i,k,j=N/2,l;*/
    struct complex v,w,t ;
    nv2=N/2 ;
    f=N ;
    for(m=1;(f=f/2)!=1;m++) /* N=2^m */
    {
        ;
    }
    nm1=N-1 ;
   
    /*变址运算*/
    for(i=1;i<=nm1;i++)
    {
        if(i<j)
        {
            t=xin[j]; xin[j]=xin[i]; xin[i]=t ;
        }
        k=nv2 ;
        while(k<j)
        {
            j=j-k ; k=k/2 ;
        }
        j=j+k ;
    }
   
    {
        int le,lei,ip ;
        for(l=1;l<=m;l++)
        {
            le=(int)pow(2,l);
            /* 这里用的是L而不是1  !!!! */
            lei=le/2 ;
            v.x=1.0 ; v.y=0.0 ;
            w.x=cos(M_PI/lei);
            w.y=-sin(M_PI/lei);
            for(j=1;j<=lei;j++)
            {
                /*double p=pow(2,m-l)*j;
                  double ps=2*PI/N*p;
                  w.x=cos(ps);
                  w.y=-sin(ps);*/
                for(i=j;i<=N;i=i+le)
                {
                    /*  w.x=cos(ps); w.y=-sin(ps);*/
                    ip=i+lei ;
                    t=EE(xin[ip],v);
                    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=1;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=1 ;
    /*******赋初值**********/
    for(;i<Num+1;i++)
    {
        s[i].x=i-1;
        s[i].y=0 ;
    }
    /***********************/
   
    FFT(s,Num);
   
    prntft(s,Num);
    getchar();
}

评论