/* 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();
}

评论