FFT算法设计与实现

FFT算法研究报告

1、程序设计背景(FFT算法理解)

FFT(fast fourier transformation),快速傅里叶变换。是对DFT算法的改进，其利用了WNnk的周期性、共轭对称性和可约性，使得DFT中有些项可以合并，大大减小了计算量。

按输入序列在时间上的次序是属于偶数还是奇数来分解称为“按时间抽取法”(DIT)。另一种是把输出序列X(k)按顺序的奇偶分解为越来越短的序列，称为按频率抽样的FFT算法（DIF）。DIT算法是先作复乘后作加减，而DIF的复乘只出现在减法之后。本次程序采用DIT算法实现FFT。

用c语言实现FFT的难点在于数据倒位序的处理，以及各级蝶形运算的实现。倒位序的实现可以使用“反向进位加法”，即倒位序二进制数的下面一个数是上面一个数在最高位加一并由高位向低位进位而得到的。

对于点数为N = 2^L的FFT运算，可以分解为L阶蝶形图级联，第M阶蝶形图内又分为2^(L-M)个蝶形组，每个蝶形组内包含

2^(M-1)个蝶形。而且旋转因子与蝶形阶数和蝶形分组内的蝶形个数存在关联。因此我们就可以构造循环来实现蝶形运算。

2、FFT算法流程图

3、FFT程序代码

#include

#include

#include

#include

#define PI 3.141592653

#define SIZE 16

#define MAX 10

//定义复数结构体

typedef struct Complex{

double real;

double imag;

}comp;

//定义复数乘法，加法，减法

void Add_(comp * a1,comp *a2,comp *b){ b->imag=a1->imag+a2->imag;

b->real=a1->real+a2->real;

}

void Sub_(comp * a1,comp *a2,comp *b){ b->imag=a1->imag-a2->imag;

b->real=a1->real-a2->real;

}

void Mul_(comp * a1,comp *a2,comp *b){ double r1=0.0,r2=0.0;

double i1=0.0,i2=0.0;

r1=a1->real;

r2=a2->real;

i1=a1->imag;

i2=a2->imag;

b->imag=r1*i2+r2*i1;

b->real=r1*r2-i1*i2;

}

//利用srand()、rand()随机生成一个输入并显示数据

void Input(double *data,int n){

printf("输入信号为:\n");

srand((int)time(0));

for(int i=0;i

data[i]=rand()%MAX;

printf("%lf\n",data[i]);

}

}

//定义WN的n次方项

void WN(double n,double size_n,comp * b){

double x=2.0*PI*n/size_n;

b->imag=-sin(x);

b->real=cos(x);

}

//处理FFT的数据位置，输入倒位序，输出自然顺序（正序），用递归结构实现int reverse(double * a,int n){

if(n==1)

return 0;

double * temp=(double *)malloc(sizeof(double)*n);

for(int i=0;i

if(i%2==0)

temp[i/2]=a[i];

else

temp[(n+i)/2]=a[i];

for(int i=0;i

a[i]=temp[i];

free(temp);

reverse(a, n/2);

reverse(a+n/2, n/2);

return 1;

}

//定义FFT函数

void FFT(double * a,comp * b,int n){

reverse(a, n); //对输入信号进行倒位序处理

int k=n;

int m=0;

while (k/=2) { //记录需要蝶形运算的级数

m++; }

k=m;

comp * a_comp=(comp*)malloc(sizeof(comp)*n);

for(int i=0;i

a_comp[i].real=a[i];

a_comp[i].imag=0; //将输入数据赋值给 a_comp }

for(int i=0;i

z=0;

for(int j=0;j

if((j/(2^(i-1)))%2==1){

comp wn;

WN(z, n, &wn);

Mul_(&a_comp[j], &wn,&a_comp[j]);

z+=2^(k-i-2);

comp temp;

int company=j-(2^(i-1));

temp.real=a_comp[company].real;

temp.imag=a_comp[company].imag;

Add_(&temp, &a_comp[j], &a_comp[company ]); Sub_(&temp, &a_comp[j], &a_comp[j]);

}

else

m=0; }

}

printf("\n\nFFT处理结果:\n");

for(int i=0;i

if(a_comp[i].imag>=0.0){

printf("%lf+%lfj\n",a_comp[i].real,a_comp[i].imag);

}

else

printf("%lf%lfj\n",a_comp[i].real,a_comp[i].imag);

for(int i=0;i

b[i].imag=a_comp[i].imag;

b[i].real=a_comp[i].real;

}

}

int main() {

double array[SIZE];