使用C语言实现FFT算法

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

#define N 1024

#define PI 3.1415



/*复数结构体类型*/

typedef struct{

double realPart;

double imaginaryPart;

}complexNumber;



complexNumber x[N], W[N/2]; //序列数组以及旋转因子数组

 int size_x=0;               //输入序列的长度，只能为2的次幂 

void fastFourierOperation();//快速傅里叶变换算法 

void initRotationFactor();  //生成旋转因子数组 

void changePosition();      //偶奇变换算法 

void outputArray();         //遍历输出数组

void add(complexNumber ,complexNumber ,complexNumber*); //复数加法 

void mul(complexNumber ,complexNumber ,complexNumber*); //复数乘法 

void sub(complexNumber ,complexNumber ,complexNumber*); //复数减法 





int main()

{

        printf("请输入序列的长度：");              //输入序列的大小

        scanf("%d",&size_x);

        printf("\n请分别输入序列的实部和虚部\n");  //输入序列的实部和虚部

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

        {

                printf("请输入第%d个序列的实部：",i);

                scanf("%lf",&x[i].realPart);

                printf("请输入第%d个序列的虚部：",i);

                scanf("%lf",&x[i].imaginaryPart);

        }

        printf("\n输出原序列：\n");

        outputArray();

        initRotationFactor();  //生成旋转因子数组 

        fastFourierOperation();//调用快速傅里叶变换

        printf("\n输出FFT后的结果：\n");

        outputArray();         //遍历输出数组

        return 0;

}





/*快速傅里叶变换*/

void fastFourierOperation()

{

        int l=0;

        complexNumber up,down,product;

        changePosition();                  //偶奇变换算法

        for(int i=0;i<size_x/2;i++)        //一级蝶形运算

        {   

                l=1<<i;

                for(int j=0;j<size_x;j+= 2*l ) //一组蝶形运算

                {            

                        for(int 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 initRotationFactor() 

{

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

        {

                W[i].realPart=cos(2*PI/size_x*i);//用欧拉公式计算旋转因子

                W[i].imaginaryPart=-1*sin(2*PI/size_x*i);

        }

}





/*偶奇变换算法*/

void changePosition()      

{

        int j=0,k;//待会儿进行运算时需要用到的变量 

        complexNumber temp;



        for (int i=0;i<size_x-1;i++) 

        { 

                if(i<j) 

                {

                        //若倒序数大于顺序数，进行变址（换位）；否则不变，防止重复交换，变回原数

                        temp=x[i];

                        x[i]=x[j];

                        x[j]=temp;

                }

                k=size_x/2;    //判断j的最高位是否为0

                while(j>=k) 

                {              //最高位为1

                        j=j-k;

                        k=k/2;

                }

                j=j+k;        //最高位为0

        }

          printf("\n输出倒序后的结果：\n");

        outputArray();

}





/*遍历输出数组*/

void outputArray()

{

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

        {

                if(x[i].imaginaryPart<0){

                        printf("%.4f-j%.4f\n",x[i].realPart,abs(x[i].imaginaryPart));

                }

                else{

                        printf("%.4f+j%.4f\n",x[i].realPart,x[i].imaginaryPart);

                }

        }

}



/*复数加法的定义*/ 

void add(complexNumber a,complexNumber b,complexNumber *c)

{

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

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

}



/*复数乘法的定义*/ 

void mul(complexNumber a,complexNumber b,complexNumber *c)

{

        c->realPart=a.realPart*b.realPart - a.imaginaryPart*b.imaginaryPart;

        c->imaginaryPart=a.realPart*b.imaginaryPart + a.imaginaryPart*b.realPart;

}



/*复数减法的定义*/ 

void sub(complexNumber a,complexNumber b,complexNumber *c)

{

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

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

}