C语言实现FFT(快速傅里叶变换）

#include <iom128.h>

#include <intrinsics.h>

/*********************************************************************

快速福利叶变换 C 函数

函数简介：此函数是通用的快速傅里叶变换 C 语言函数，移植性强，以下部分不依

赖硬件。此函数采用联合体的形式表示一个复数，输入为自然顺序的复

数（输入实数是可令复数虚部为 0），输出为经过 FFT 变换的自然顺序的

复数

**********************************************************************/

#include<math.h>

#define PI 3.1415926535897932384626433832795028841971 //定义圆周率值

#define FFT_N 128 //定义福利叶变换的点数

struct compx {float real,imag;}; //定义一个复数结构

struct compx s[FFT_N]; //FFT 输入和输出：从 S[1]开始存放，根

据大小自己定义

/*******************************************************************

函数原型：struct compx EE(struct compx b1,struct compx b2)

struct compx c;

c.real=a.real*b.real-a.imag*b.imag;

c.imag=a.real*b.imag+a.imag*b.real;

return(c);

}

/*****************************************************************

函数原型：void FFT(struct compx *xin,int N)

函数功能：对输入的复数组进行快速傅里叶变换（FFT）

输入参数：*xin 复数结构体组的首地址指针，struct 型

*****************************************************************/

void FFT(struct compx *xin)

{

nv2=FFT_N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法

{

t=xin[j];

xin[j]=xin[i];

xin[i]=t;

}

k=nv2; //求 j 的下一个倒位序

while(k<=j) //如果 k<=j,表示 j 的最高位为 1

{

j=j-k; //把最高位变成 0

k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为 0

}

j=j+k; //把 0 改为 1

}

lei=le/2; //同一蝶形结中参加运算的两点的距离

u.real=1.0; //u 为蝶形结运算系数，初始值为 1

u.imag=0.0;

w.real=cos(PI/lei); //w 为系数商，即当前系数与前一个系数的商

w.imag=-sin(PI/lei);

for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结

{

for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算，即计算系数相同蝶形结

{

ip=i+lei; //i，ip 分别表示参加蝶形运算的两个节点

t=EE(xin[ip],u); //蝶形运算，详见公式

xin[ip].real=xin[i].real-t.real;

xin[ip].imag=xin[i].imag-t.imag;

xin[i].real=xin[i].real+t.real;

u=EE(u,w); //改变系数，进行下一个蝶形运算

}

}

/************************************************************

函数功能：测试 FFT 变换，演示函数使用方法

输出参数：无

void main()

for(i=0;i<FFT_N;i++) //给结构体赋值

s[i].real=sin(2*3.141592653589793*i/FFT_N); //实部为正弦波 FFT_N 点采样，赋值为 1

FFT(s); //进行快速福利叶变换

s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);