FFT程序的设计.doc

1、FFT程序的设计一、设计内容1.1 设计要求 用C程序平台设计快速傅里叶变换的程序，既要有设计的理论内容，也要有相应的流程图，及源代码。1.2 设计目的 1.掌握了解快速傅里叶变换程序的设计。 2.学习C语言平台的使用方法及基本功能。 3.加深对快速傅里叶变换的理解，熟悉快速傅里叶变换子程序调用。 4.熟悉应用快速傅里叶变换实现两个序列线性卷积的方法二、设计思路2.1 设计原理2.1.1 FFT算法原理FFT并不是另一种变换傅里叶变换，而是DFT的一种快速算法。它通过对DFT变换式的一次次分解，将长序列的DFT分解为一系列短序列DFT的组合，从而减少运算量。常用的FFT是以2为基数的，其长度N

2、=2L。当要变换的序列长度不等于2的整数次方时，是为了使用以2为基数的FFT，可以用末位补零的方法，使其长度延长至2的整数次方。2.1.2 FFT的运算量 对于一个N长的序列，直接计算DFT需要N2次复数乘及N(N-1)=N2复数加。 按时间抽选法FFT，当N=2L时，共有L级蝶形，每级都由N/2个蝶形运算组成，每个蝶形有一次复数乘、二次复数加，总的运算量为（N/2）L次复数乘，NL次复数加。2.1.3 用FFT计算线性卷积 用FFT可以实现两个序列的圆周卷积。在一定的条件下，可以使圆周卷积等于线性卷积。一般情况，设两个序列的长度分别为N1和N2，要使圆周卷积等于线性卷积的充分必要条件是FFT

3、的长度N N1+N2-1对于长度不足N的两个序列，分别将他们补零延长到N。2.2 C程序设计法及程序流程图x3(0)m=2x2(0)m=1x1(0)m=0x0(0) FFT算法是离散傅立叶变换的快速算法，其结果应该与离散傅立叶变换所算出的结果一致。图1是采样点数为8的FFT算法的一种形式，本程序也是以这种形式设计的。x3(1)x2(1)x1(1)x0(1)x(0) X(0)x3(2)x2(2)x1(2)x0(2)-1x(4) X(1)x3(3)x2(3)-1x1(3)x0(3)x(2) X(2)x3(4)x2(4)-1x1(4)-1x0(4)x(6) X(3)x3(5)-1x2(5)x1(5)

4、x0(5)x(1) X(4)x3(6)-1x2(6)x1(6)-1x0(6)x(5) X(5)x3(7)-1x2(7)-1x1(7)x0(7)x(3) X(6)-1-1-1x(7) X(7) 图1 FFT算法形式 程序设计的基本思路是在程序开始时刻要求输入采样点数，如果采样点数不是2的整数次方（不包括0次方），则要求输入采样点的数值，并根据采样点数分配响应的数组大小，计算迭代次数。然后对采样点进行逆二进制排序，再按上图所示的算法进行计算，程序流程图如图2：输入采样点数程序开始采样点数是否为2的不为0的整数次方？NY根据输入点数分配数组大小按迭代，分组，碟形运算的顺序进行运算 输出计算结果 程序

5、结束 图2 流程图 本程序在VC+6.0的开发环境下创建Win32 Console Application工程对程序进行设计，下面分别对程序设计中复数类的应用，判断和求迭代次数，逆二进制排序，蝶形运算进行具体说明。2.3 快速傅里叶变换程序设计各类的应用2.3.1 复数类的应用 C+语言本身并不包含复数数据类型，但在经过的不断扩展后，如今的STL（Standard Template Library）中已经包含了相应的复数类，但在尝试后发现其功能仍不能满足本程序设计的要求，所以仍需要在其基础上编写函数。在FFT程序设计中，复数类主要被用来计算两复数的加法和乘法以及旋转因子，其中，式中N=2的i次

6、方，i代表第i次迭代，而k代表第k次蝶形运算，由于STL中的Complex类不能进行带参数的复数运算，所以本程序编写了带参数的复数运算cW（int cm,int sign2m），用于计算，设计的基本思路是化复数乘法和除法运算为几个复数基本单元的加法运算和除法运算，其中加法运算的次数由函数输入参量int cm决定，复数除法运算次数由int sign2m决定。基本加法单元为复数类对象wmt(0,-1*2*pi)，基本除法类对象为pm(2,0)。下面是cW（int cm,int sign2m）的代码： complex cW(int cm,int sign2m)int s;complex Wm(0,0

7、),pm(2,0),Wmt(0,-1*2*pi);for(s=0;scm;s+)Wm=Wm+Wmt;for(s=0;ssign2m;s+)Wm=Wm/pm;return Wm;2.3.1.1 判断和求迭代次数 本程序编写iternumb(int numb)函数对采样点数进行判断，如果采样点数不符合2的整数次方或采样点数为1或0，则跳出程序，程序设计基本思路是对输入采样点数的十进制形式进行模2运算和除法运算，在除法运算结果大于1之前，一旦模2运算的结果等于1，则说明输入采样点数不符合要求，而如果符合要求，则把出发结果存入数组当中，函数代码如下：int iternumb(int numb)int

8、iternumb=0;if(numb=0)|(numb=1)coutnum error!n;exit(0);while(numb!=0)&(numb!=1) if(numb%2)coutnum error!n;exit(0);numb=numb/2;iternumb=iternumb+1;return iternumb;2.3.1.2 逆二进制排序在逆二进制排序程序中，设置for循环分别将输入数据数组inputi的索引号i进行模2运算，所得的结果按逆序存入inverse数组（存入inverse数组的顺序是从数组尾部开始）。然后再将inverse中的二进制数转换为十进制数，并以此数为Aj的索引号

9、，令Aj= inputi，从而实现了逆二进制排序，代码如下：for(a=0;anum;a+)coutxatemp;inputa=temp;for(a=0;anum;a+)temp2=a;for(b=0;biternum;b+)temp4=temp2%2;inverseiternum-1-b=temp4;temp2=temp2/2;temp1=0;sign=1;for(c=0;citernum;c+)temp1=sign*inversec+temp1;sign=sign*2;Atemp1=inputa;2.3.2 蝶形运算 蝶形运算是FFT中最基本的一个运算单元。在FFT程序设计中要找到蝶形运算

10、地址与第几次迭代，第几组之间的关系。根据FFT算法的特点，本程序设置了3个for循环的嵌套循环分别表示迭代，分组和蝶形运算，经过总结的得出蝶形运算地址与迭代序号，分组序号间的关系如下： (1) 上式中表示前一次迭代运算的结果，i表示迭代序号，j表示分组序号，k表示蝶形运算序号，在程序中分别用a，b，c表示。 根据上式(1)的关系，设置sign2=，sign1=每组表示蝶形运算的次数，sign表示该次迭代的分组数，因为这些变量都与第几次迭代有关，所以最外层的for循环每执行依次，这些变量也相应更新一次。这一部分的代码如下：sign=num/2;sign1=1;sign2=2;for(a=1;ai

11、ternum+1;a+)for(b=1;bsign+1;b+)for(c=0;csign1;c+)W=cW(c,a);tempA1=Ac+sign2*(b-1)+Ac+sign2*(b-1)+sign2/2*exp(W);tempA2=Ac+sign2*(b-1)-Ac+sign2*(b-1)+sign2/2*exp(W);Ac+sign2*(b-1)=tempA1;Ac+sign2*(b-1)+sign2/2=tempA2;sign1=sign1*2;sign2=sign2*2;sign=sign/2;2.4 实验源程序程序一：#include #include #include using

12、 namespace std;#define pi 4*atan(1)int iternumb(int numb)int iternumb=0;if(numb=0)|(numb=1)coutnum error!n;exit(0);while(numb!=0)&(numb!=1) if(numb%2)coutnum error!n;exit(0);numb=numb/2;iternumb=iternumb+1;return iternumb;complex cW(int cm,int sign2m)int s;complex Wm(0,0),pm(2,0),Wmt(0,-1*2*pi);for(

13、s=0;scm;s+)Wm=Wm+Wmt;for(s=0;ssign2m;s+)Wm=Wm/pm;return Wm;int main()int a,b,c=0,sign,sign1,sign2,num,iternum,temp1,temp2,temp4;intiternumb(int numb);double temp;complex tempA1,tempA2,W;coutnum; iternum=iternumb(num);double* input=(double*)malloc(sizeof(double)*num);complex* A=(complex*)malloc(sizeo

14、f(complex)*num);int* inverse=(int*)malloc(sizeof(int)*iternum);for(a=0;anum;a+)coutxatemp;inputa=temp;for(a=0;anum;a+)temp2=a;for(b=0;biternum;b+)temp4=temp2%2;inverseiternum-1-b=temp4;temp2=temp2/2;temp1=0;sign=1;for(c=0;citernum;c+)temp1=sign*inversec+temp1;sign=sign*2;Atemp1=inputa;sign=num/2;sig

15、n1=1;sign2=2;for(a=1;aiternum+1;a+)for(b=1;bsign+1;b+)for(c=0;csign1;c+)W=cW(c,a);tempA1=Ac+sign2*(b-1)+Ac+sign2*(b-1)+sign2/2*exp(W);tempA2=Ac+sign2*(b-1)-Ac+sign2*(b-1)+sign2/2*exp(W);Ac+sign2*(b-1)=tempA1;Ac+sign2*(b-1)+sign2/2=tempA2;sign1=sign1*2;sign2=sign2*2;sign=sign/2;for(a=0;anum;a+)coutXa

16、=Aan;return 0;程序二：基-2FFT算法的C+程序#includeint fft(complex * a,int l)const double pai=3.141592653589793complex u,w,t;unsigned n=1,nv2,nm1,k,le,lei,ip;unsigned i,j,m;double tmp;n1;nm1=n-1;i=0;for(i=0;i,nm1,i+) if(i,j) t=aj aj=ai; ai=t; k=nv2; while(k=1 j+=kle=1for(m=1:m=1:m+) lei=le le=1; u=complex(1,0);

17、 tmp=pai/lei; w=complex(cos(tmp),-sin(tmp); for(j=0;jlei;j+) for(i=j;in;i+=le) ip=i+lei; t=aip*u; aip=ai-t; ai+=t; u*=w; return 0;2.5 程序运行结果输入采样点数为4，数据分别为1,3,5,6.图3 C程序结果 由图3知所得结果为X0=15，X1=-4+3i，X2=-3，X3=-4-3i 在Matlab中进行检验：图 4 Matlab验证结果由图4结果与所设计程序一致。输入采样点数为8分别为2,3,4,5,7,9,10,11. 详见图5：图5 C程序结果在Matla

18、b中检验，如图6：图 6 Matlab验证结果与所设计程序仍然一致。结论：通过对程序的设计以及对运行结果的验证可以说明所设计的程序运行结果正确，完成了对采样点进行FFT变换的要求。（注明：由于程序二中complex头文件不是微软公司的，因此无法用VC运行。）三、心得体会学习数字信号处理已经有一个学期了，现在和大家聊一下我学习过程中的一些经历和想法吧。刚接触时，和很多同学都一样，觉得很难，很深奥，不易懂。可是其中却蕴含着许多技巧，渐渐地，对其产生了兴趣从此不能自拔。很多个日夜就这样陪伴着它度过了。期间也遇到过非常多的问题，也一度被这些问题所困惑等到回过头来，看到自己曾经走过的路，唏嘘不已。经常混

19、迹于论坛里，也看到了很多学成者发的帖子，对自己的帮助很大，也交到了一些关于这方面有所长的朋友，经常上网交流。路学习过来的过程中，帮助最大之一无疑来自于网络了。很多时候，通过网络，我们都可以获取到所需要的学习资料。但是，随着我们学习的深入，我们会慢慢发现，网络提供的东西是有限度的，好像大部分的资料都差不多，或者说是适合大部分的初学者所需，而当我们想更进一步提高时，却发现能够获取到的资料越来越少，这也需要我们的帮助之二老师，通过老师的帮助，我们的知识面得到了扩展与提升。我们这里所做的是有关快速傅里叶变换的，它只是数字信号处理中的一小部分，还有很多诸如此类的知识都是令我们受益匪浅。之前没想到C语言居然可以跟数字信号处理联系在一起，现在都得以发挥了它们的功效，完美的诠释了我们所学的内容，这也给我们启示，在以后的学习道路中，很多东西都可以联系起来，很多看似平凡的东西都有它独到的用途，使它发光发亮。四、参考文献【1】丁玉美. 数字信号处理. 电子科技大学出版社【2】周辉. 数字信号处理及其MATLAB实现. 中国林业出版社【3】程佩青. 数字信号处理教程. 清华大学出版社 【4】余成波. 数字信号处理及MATLAB实现. 清华大学出版社【5】张明照. 应用MATLAB实现信号分析和处理. 科技出版社- 16 -