FFT源程序（仅核心部分）

XJzy

#include<stdio.h>

#include<math.h>

#include <stdlib.h>

#define size_x N

typedef struct

       {

              double real;

              double img;

       }complex;

complex W[size_x/2];

complex x[size_x];

double PI;

void fft();

void initW();

void change();

void add(complexa,complex b,complex *c);

void mul(complexa,complex b,complex *c);

void sub(complexa,complex b,complex *c);

void divi(complexa,complex b,complex *c);

void output();

int main()

{

       system("cls");

       PI=atan(1)*4;

       initW();

       fft();

       output();

       return 0;

}

void fft()

{

       int i=0,j=0,k=0,l=0;

       complex up,down,product;

       change();

       for(i=0;i< (int)( log(size_x)/log(2));i++)

       {

              l=( 1<<i );

              for(j=0;j<size_x;j+=(1<<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 initW()

{

       int i;

       for(i=0;i<size_x/2;i++)

       {

              W.real=cos(2*PI/size_x*i);

              W.img=-1*sin(2*PI/size_x*i);

       }

}

void change()

{

       complex temp;

       int i=0,j=0,k=0,t;

       for(i=0;i<size_x;i++)

       {

              k=i;j=0;

              t=(unsigned) (log(size_x)/log(2));

              while(t--)

              {

                     j=j<<1;

                     j|=(k & 1);

                     k=k>>1;

              }

              if(j>i)

              {

                     temp=x;

                     x=x[j];

                     x[j]=temp;

              }

       }

}

void output()

{

       int i;

       printf("The result are asfollows\n");

       for(i=0;i<size_x;i++)

       {

              printf("%.4f",x.real);

              if(x.img>=0.0001)printf("+%.4fj\n",x.img);

              elseif(fabs(x.img)<0.0001)printf("\n");

              elseprintf("%.4fj\n",x.img);

       }

}

void add(complexa,complex b,complex *c)

{

       c->real=a.real+b.real;

       c->img=a.img+b.img;

}

void mul(complex a,complexb,complex *c)

{

       c->real=a.real*b.real - a.img*b.img;

       c->img=a.real*b.img + a.img*b.real;

}

void sub(complexa,complex b,complex *c)

{

       c->real=a.real-b.real;

       c->img=a.img-b.img;

}