用c实现的fft和ifft

对百度文库中coffealove的代码进行了如下验证
代码如下

#include   <stdio.h>
#include   <math.h>
#include   <stdlib.h>


#define   N   1000
typedef   struct
{
    double   real;
    double   img;
}complex;
void   fft(); /*快速傅里叶变换*/
void   ifft(); /*快速傅里叶逆变换*/
void   initW();
void   change();
void   add(complex   ,complex   ,complex   *);   /*复数加法*/
void   mul(complex   ,complex   ,complex   *);   /*复数乘法*/
void   sub(complex   ,complex   ,complex   *);   /*复数减法*/
void   divi(complex   ,complex   ,complex   *);/*复数除法*/
void   output(); /*输出结果*/

complex   x[N],   *W;/*输出序列的值*/
int   size_x=0;/*输入序列的长度，只限2的N次方*/
double   PI;

int   main()
{
    int   i,method;
    PI=atan(1.0)*4;
    printf("Please input the size of x:\n");
    /*输入序列的长度*/
    scanf("%d",&size_x);
    printf("Please   input   the   data   in   x[N]:(such as:5 6)\n");
    /*输入序列对应的值*/
    for(i=0;i<size_x;i++)
        scanf("%lf %lf",&x[i].real,&x[i].img);
    initW();
    /*选择FFT或逆FFT运算*/
    printf("Use   FFT(0)   or   IFFT(1)?\n");
    scanf("%d",&method);
    if(method==0)
        fft();
    else
        ifft();
    output();
    return   0;
}

/*进行基-2 FFT运算*/
void   fft()
{
    int   i=0,j=0,k=0,l=0;
    complex   up,down,product;
    change();
    for(i=0;i<   log(double(size_x))/log(double(2));i++)
    {
        l=1<<i;
        for(j=0;j<size_x;j+=   2*l   )
        {
            for(k=0;k<l;k++)
            {
                mul(x[j+k+l],W[size_x*k/2/l],&product);
                add(x[j+k],product,&up);
                sub(x[j+k],product,&down);
                x[j+k]=up;
                x[j+k+l]=down;
            }
        }
    }
}


void   ifft()
{
    int   i=0,j=0,k=0,l=size_x;
    complex   up,down;
    for(i=0;i<     log(double(size_x))/log(double(2));i++) /*蝶形运算*/
    {
        l/=2;
        for(j=0;j<size_x;j+=   2*l   )
        {
            for(k=0;k<l;k++)
            {
                add(x[j+k],x[j+k+l],&up);
                up.real/=2;up.img/=2;
                sub(x[j+k],x[j+k+l],&down);
                down.real/=2;down.img/=2;
                divi(down,W[size_x*k/2/l],&down);
                x[j+k]=up;
                x[j+k+l]=down;
            }
        }
    }
    change();
}


void   initW()
{
    int   i;
    W=(complex   *)malloc(sizeof(complex)   *   size_x);
    for(i=0;i<size_x;i++)
    {
        W[i].real=cos(2*PI/size_x*i);
        W[i].img=-1*sin(2*PI/size_x*i);
    }
}


void   change()
{
    complex   temp;
    unsigned   short   i=0,j=0,k=0;
    double   t;
    for(i=0;i<size_x;i++)
    {
        k=i;j=0;
        t=(log(double(size_x))/log(double(2)));
        while(   (t--)>0   )
        {
            j=j<<1;
            j|=(k & 1);
            k=k>>1;
        }
        if(j>i)
        {
            temp=x[i];
            x[i]=x[j];
            x[j]=temp;
        }
    }
}


void   output()   /*输出结果*/
{
    int   i;
    printf("The   result   are   as   follows\n");
    for(i=0;i<size_x;i++)
    {
        printf("%.4f",x[i].real);
        if(x[i].img>=0.0001)
            printf("+%.4fj\n",x[i].img);
        else   if(fabs(x[i].img)<0.0001)
            printf("\n");
        else
            printf("%.4fj\n",x[i].img);
    }
}
void   add(complex   a,complex   b,complex   *c)
{
    c->real=a.real+b.real;
    c->img=a.img+b.img;
}

void   mul(complex   a,complex   b,complex   *c)
{
    c->real=a.real*b.real   -   a.img*b.img;
    c->img=a.real*b.img   +   a.img*b.real;
}
void   sub(complex   a,complex   b,complex   *c)
{
    c->real=a.real-b.real;
    c->img=a.img-b.img;
}
void   divi(complex   a,complex   b,complex   *c)
{
    c->real=(   a.real*b.real+a.img*b.img   )/(
        b.real*b.real+b.img*b.img);
    c->img=(   a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);
}

验证如下