FFT(快速傅里叶) c语言版

#include <math.h>  
#include <stdlib.h>  
#define N 1000  
/*定义复数类型*/  
typedef struct{  
double real;  
double img;  
}complex;  
  
  
complex x[N], *W; /*输入序列,变换核*/  
int size_x=0;      /*输入序列的大小，在本程序中仅限2的次幂*/  
double PI;         /*圆周率*/  
void fft();     /*快速傅里叶变换*/  
void initW();   /*初始化变换核*/  
void change(); /*变址*/  
void add(complex ,complex ,complex *); /*复数加法*/  
void mul(complex ,complex ,complex *); /*复数乘法*/  
void sub(complex ,complex ,complex *); /*复数减法*/  
void output();/*输出快速傅里叶变换的结果*/  
  
  
  
  
int main()  
{  
    int i;                             /*输出结果*/  
    system("cls");  
    PI=atan(1)*4;  
    printf("                                        输出DIT方法实现的FFT结果\n");  
    printf("Please input the size of x:\n");//输入序列的大小  
    scanf("%d",&size_x);  
    printf("Please input the data in x[N]:\n");//输入序列的实部和虚部  
    for(i=0;i<size_x;i++)  
    {  
        printf("请输入第%d个序列：",i);  
        scanf("%lf%lf",&x[i].real,&x[i].img);  
    }  
    printf("输出倒序后的序列\n");  
    initW();//调用变换核  
    fft();//调用快速傅里叶变换  
    printf("输出FFT后的结果\n");  
    output();//调用输出傅里叶变换结果函数  
    return 0;  
}  
  
  
  
/*快速傅里叶变换*/  
void fft()  
{  
    int i=0,j=0,k=0,l=0;  
    complex up,down,product;  
    change();  //调用变址函数  
    for(i=0;i< log(size_x)/log(2) ;i++)        /*一级蝶形运算 stage */  
    {     
        l=1<<i;  
        for(j=0;j<size_x;j+= 2*l )     /*一组蝶形运算 group,每个group的蝶形因子乘数不同*/  
        {              
            for(k=0;k<l;k++)        /*一个蝶形运算 每个group内的蝶形运算*/  
            {         
                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;  
            }  
        }  
    }  
}  
  
  
  
/*初始化变换核，定义一个变换核，相当于旋转因子WAP*/  
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);  
    }  
}  
  
  
  
/*变址计算，将x(n)码位倒置*/  
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(size_x)/log(2));  
  while( (t--)>0 )    //利用按位与以及循环实现码位颠倒  
  {  
    j=j<<1;  
    j|=(k & 1);  
    k=k>>1;  
  }  
  if(j>i)    //将x(n)的码位互换  
  {  
    temp=x[i];  
    x[i]=x[j];  
    x[j]=temp;  
  }  
  }  
  output();  
}  
  
  
  
/*输出傅里叶变换的结果*/  
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;  
}
`