FFT程序的设计

-1-FFT程序的设计一、设计内容1.1设计要求用C程序平台设计快速傅里叶变换的程序，既要有设计的理论内容，也要有相应的流程图，及源代码。1.2设计目的1.掌握了解快速傅里叶变换程序的设计。2.学习C语言平台的使用方法及基本功能。3.加深对快速傅里叶变换的理解，熟悉快速傅里叶变换子程序调用。4.熟悉应用快速傅里叶变换实现两个序列线性卷积的方法二、设计思路2.1设计原理2.1.1FFT算法原理FFT并不是另一种变换傅里叶变换，而是DFT的一种快速算法。它通过对DFT变换式的一次次分解，将长序列的DFT分解为一系列短序列DFT的组合，从而减少运算量。常用的FFT是以2为基数的，其长度N=2L。当要变换的序列长度不等于2的整数次方时，是为了使用以2为基数的FFT，可以用末位补零的方法，使其长度延长至2的整数次方。-2-2.1.2FFT的运算量对于一个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的长度N≥N1+N2-1对于长度不足N的两个序列，分别将他们补零延长到N。2.2C程序设计法及程序流程图FFT算法是离散傅立叶变换的快速算法，其结果应该与离散傅立叶变换所算出的结果一致。图1是采样点数为8的FFT算法的一种形式，本程序也是以这种形式设计的。x(0)X(0)x(4)X(1)x(2)X(2)x(6)X(3)x(1)X(4)x(5)X(5)x(3)X(6)x(7)X(7)图1FFT算法形式程序设计的基本思路是在程序开始时刻要求输入采样点数，如果采样点数不是2的整数次方（不包括0次方），则要求输入采样点的数值，并根据采样点数分配响应的数组大小，计算迭代次数。然后对采样点进行逆二进制排序，再按上图所示的算法进行计算，程x0(1)-1x0(2)x0(3)x0(4)x0(5)x0(6)x0(7)-1-1-1x1(1)x1(2)x1(3)x1(4)x1(5)x1(6)x1(7)08W28W08W28W-1-1-1-1x2(1)x2(2)x2(3)x2(4)x2(5)x2(6)x2(7)08W18W28W38W-1-1-1-1x3(1)x3(2)x3(3)x3(4)x3(5)x3(6)x3(7)x0(0)m=0x1(0)m=1x2(0)m=2x3(0)-3-序流程图如图2：图2流程图本程序在VC++6.0的开发环境下创建Win32ConsoleApplication工程对程序进行设计，下面分别对程序设计中复数类的应用，判断和求迭代次数，逆二进制排序，蝶形运算进行具体说明。2.3快速傅里叶变换程序设计各类的应用2.3.1复数类的应用C++语言本身并不包含复数数据类型，但在经过的不断扩展后，如今的STL（StandardTemplateLibrary）中已经包含了相应的复数类，但在尝试后发现其功能仍不能满足本程序设计的要求，所以仍需要在其基础上编写函数。程序开始输入采样点数采样点数是否为2的不为0的整数次方？根据输入点数分配数组大小按迭代，分组，碟形运算的顺序进行运算输出计算结果程序结束NY-4-在FFT程序设计中，复数类主要被用来计算两复数的加法和乘法以及旋转因子KNW，其中lNjNW2，式中N=2的i次方，i代表第i次迭代，而k代表第k次蝶形运算，由于STL中的Complex类不能进行带参数的复数运算，所以本程序编写了带参数的复数运算cW（intcm,intsign2m），用于计算ikj22，设计的基本思路是化复数乘法和除法运算为几个复数基本单元的加法运算和除法运算，其中加法运算的次数由函数输入参量intcm决定，复数除法运算次数由intsign2m决定。基本加法单元为复数类对象wmt(0,-1*2*pi)，基本除法类对象为pm(2,0)。下面是cW（intcm,intsign2m）的代码：complexdoublecW(intcm,intsign2m){ints;complexdoubleWm(0,0),pm(2,0),Wmt(0,-1*2*pi);for(s=0;scm;s++){Wm=Wm+Wmt;}for(s=0;ssign2m;s++){Wm=Wm/pm;}returnWm;}2.3.1.1判断和求迭代次数本程序编写iternumb(intnumb)函数对采样点数进行判断，如果采样点数不符合2的整数次方或采样点数为1或0，则跳出程序，程序设计基本思路是对输入采样点数的十进制形式进行模2运算和除法运算，在除法运算结果大于1之前，一旦模2运算-5-的结果等于1，则说明输入采样点数不符合要求，而如果符合要求，则把出发结果存入数组当中，函数代码如下：intiternumb(intnumb){intiternumb=0;if((numb==0)||(numb==1)){coutnumerror!\n;exit(0);}while((numb!=0)&&(numb!=1)){if(numb%2){coutnumerror!\n;exit(0);}numb=numb/2;iternumb=iternumb+1;}returniternumb;}2.3.1.2逆二进制排序在逆二进制排序程序中，设置for循环分别将输入数据数组input[i]的索引号i进行模2运算，所得的结果按逆序存入inverse[]数组（存入inverse[]数组的顺序是从数组尾部开始）。-6-然后再将inverse[]中的二进制数转换为十进制数，并以此数为A[j]的索引号，令A[j]=input[i]，从而实现了逆二进制排序，代码如下：for(a=0;anum;a++){coutxa=\n;cintemp;input[a]=temp;}for(a=0;anum;a++){temp2=a;for(b=0;biternum;b++){temp4=temp2%2;inverse[iternum-1-b]=temp4;temp2=temp2/2;}temp1=0;sign=1;for(c=0;citernum;c++){temp1=sign*inverse[c]+temp1;sign=sign*2;}A[temp1]=input[a];}-7-2.3.2蝶形运算蝶形运算是FFT中最基本的一个运算单元。在FFT程序设计中要找到蝶形运算地址与第几次迭代，第几组之间的关系。根据FFT算法的特点，本程序设置了3个for循环的嵌套循环分别表示迭代，分组和蝶形运算，经过总结的得出蝶形运算地址与迭代序号，分组序号间的关系如下：kijjjjjkiiiWjkAjkAjkAWjkAjkAjkA2111112111*]2)1(2[]2)1(2[]2)1(2[*]2)1(2[)]1(2[)]1(2[(1)上式中1A表示前一次迭代运算的结果，i表示迭代序号，j表示分组序号，k表示蝶形运算序号，在程序中分别用a，b，c表示。根据上式(1)的关系，设置sign2=i2，sign1=12i每组表示蝶形运算的次数，sign表示该次迭代的分组数，因为这些变量都与第几次迭代有关，所以最外层的for循环每执行依次，这些变量也相应更新一次。这一部分的代码如下：sign=num/2;sign1=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=A[c+sign2*(b-1)]+A[c+sign2*(b-1)+sign2/2]*exp(W);tempA2=A[c+sign2*(b-1)]-A[c+sign2*(b-1)+sign2/2]*exp(W);A[c+sign2*(b-1)]=tempA1;A[c+sign2*(b-1)+sign2/2]=tempA2;-8-}}sign1=sign1*2;sign2=sign2*2;sign=sign/2;}2.4实验源程序程序一：#includeiostream#includecomplex#includemath.husingnamespacestd;#definepi4*atan(1)intiternumb(intnumb){intiternumb=0;if((numb==0)||(numb==1)){coutnumerror!\n;exit(0);}while((numb!=0)&&(numb!=1)){-9-if(numb%2){coutnumerror!\n;exit(0);}numb=numb/2;iternumb=iternumb+1;}returniternumb;}complexdoublecW(intcm,intsign2m){ints;complexdoubleWm(0,0),pm(2,0),Wmt(0,-1*2*pi);for(s=0;scm;s++){Wm=Wm+Wmt;}for(s=0;ssign2m;s++){Wm=Wm/pm;}returnWm;}intmain(){inta,b,c=0,sign,sign1,sign2,num,iternum,temp1,temp2,temp4;-10-intiternumb(intnumb);doubletemp;complexdoubletempA1,tempA2,W;coutPleaseenterinputpointnumber\n;cinnum;iternum=iternumb(num);double*input=(double*)malloc(sizeof(double)*num);complexdouble*A=(complexdouble*)malloc(sizeof(complexdouble)*num);int*inverse=(int*)malloc(sizeof(int)*iternum);for(a=0;anum;a++){coutxa=\n;cintemp;input[a]=temp;}for(a=0;anum;a++){temp2=a;for(b=0;biternum;b++){temp4=temp2%2;inverse[iternum-1-b]=temp4;temp2=temp2/2;}temp1=0;-11-sign=1;for(c=0;citernum;c++){temp1=sign*inverse[c]+temp1;sign=sign*2;}A[temp1]=input[a];}sign=num/2;sign1=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=A[c+sign2*(b-1)]+A[c+sign2*(b-1)+sign2/2]*exp(W);tempA2=A[c+sign2*(b-1)]-A[c+sign2*(b-1)+sign2/2]*exp(W);A[c+sign2*(b-1)]=tempA1;A[c+sign2*(b-1)+sign2/2]=tempA2;}}sign1=sign1*2;-12-sign2=sign2*2;sign=sign/2;}for(a=0;anum;a++){coutXa=A[a]\n;}return0;}程序二：基-2FFT算法的C++程序#include