算法设计与实现

1、FFT算法研究报告 1、 程序设计背景(FFT算法理解)FFT(fast fourier transformation),快速傅里叶变换。是对DFT算法的改进，其利用了WNnk的周期性、共轭对称性和可约性，使得DFT中有些项可以合并，大大减小了计算量。按输入序列在时间上的次序是属于偶数还是奇数来分解称为“按时间抽取法”(DIT)。另一种是把输出序列X(k)按顺序的奇偶分解为越来越短的序列，称为按频率抽样的FFT算法（DIF）。DIT算法是先作复乘后作加减，而DIF的复乘只出现在减法之后。本次程序采用DIT算法实现FFT。用c语言实现FFT的难点在于数据倒位序的处理，以及各级蝶形运算的实现。倒位

2、序的实现可以使用“反向进位加法”，即倒位序二进制数的下面一个数是上面一个数在最高位加一并由高位向低位进位而得到的。对于点数为N = 2L的FFT运算，可以分解为L阶蝶形图级联，第M阶蝶形图内又分为2(L-M)个蝶形组，每个蝶形组内包含2(M-1)个蝶形。而且旋转因子与蝶形阶数和蝶形分组内的蝶形个数存在关联。因此我们就可以构造循环来实现蝶形运算。2、 FFT算法流程图输入数组指针a，长度n倒位序流程图：n=0? Yi=0Ni<n?N i%2=0?Y NTempi=ai/2temp(n+i)/2=ai Yi=i+1a=tempn=n/2结束Y3、 FFT程序代码#include <st

3、dio.h> #include <stdlib.h> #include <math.h> #include <time.h>#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=a

4、1->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*

5、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<SIZE;i+) datai=rand()%MAX; printf("%lfn",datai); /定义WN的n次方项void WN(double n,double size_n,comp * b) double x=2.0*PI*n/size_n; b->ima

6、g=-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<n;i+) if(i%2=0) tempi/2=ai; else temp(n+i)/2=ai; for(int i=0;i<n;i+) ai=tempi; free(temp); reverse(a, n/2); reverse

7、(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<n;i+) a_compi.real=ai; a_compi.imag=0; /将输入数据赋值给 a_comp for(int i=0;i<k;i+) /依次计算各级蝶

8、形运算 z=0; for(int j=0;j<n;j+) if(j/(2(i-1)%2=1) comp wn; WN(z, n, &wn); Mul_(&a_compj, &wn,&a_compj); z+=2(k-i-2); comp temp; int company=j-(2(i-1); temp.real=a_compcompany.real; temp.imag=a_compcompany.imag; Add_(&temp, &a_compj, &a_compcompany ); Sub_(&temp, &

9、a_compj, &a_compj); else m=0; printf("nnFFT处理结果:n"); for(int i=0;i<n;i+) if(a_compi.imag>=0.0) printf("%lf+%lfjn",a_compi.real,a_compi.imag); else printf("%lf%lfjn",a_compi.real,a_compi.imag); for(int i=0;i<n;i+) bi.imag=a_compi.imag; bi.real=a_compi.real;

10、int main() double arraySIZE; comp bSIZE; Input(array,SIZE); /随机生成16个输入数据 FFT(array, b, SIZE); 4、 与MATLAB自带FFT函数运行结果的比较自编FFT结果：将输入信号在matlab中进行FFT运算，结果如下：由以上截图可知，函数成功地完成了对于序列xn的FFT计算，但是可以看出matlab的准确度更高，应该是matlab计算使用的pi值更为精确。为比较二者的时间，将数据长度重新定义为215Matlab FFT程序： 自编FFT程序： 由于计算matlab用时使用的数据与自编FFT程序所用数据不同，可能会有不可预料的误差