基-2FFT算法的软件实现

实验二基-2FFT算法的软件实现一、实验目的1、加深对DFT算法原理和基本性质的理解；2、熟悉FFT算法的流程；3、了解FFT算法的应用。二、基本原理1、DFT算法原理（见教材第三章）2、按时间抽取（DIT）的-2FFT算法（1）算法原理序列x(n)的N(N=2-M)点DFT为knNNnNWnxnxDFTkX10)()]([)(点，k=0,1,…,N-1（2.1）将式（2.1）按n的奇偶性分解为)12(12/212/)12()2()()()(lkNNnlkNNnknNnknNnWlxWlxWnxWnxkX奇数偶数奇数偶数（2.2）令)2()(1lxlx,)12()(2lxlx，因为klNlkNWW2/2,所以式（2.2）可写成)12(2/122/12/012)()()(lkNMnkNklNNlWlxWWlxkX奇数（2.3）式（2.3）说明，按n的奇偶性将x(n)分解为两个N/2长的序列x1(l)和x2(l)，则N点DFT可分解为两个N/2点DFT来计算。用X1(k)和X2(k)分别表示12,...,1,0)()]([)(12/02/12/11NkWlxlxDFTkXNlklNN点（2.4）12,...,1,0)()]([)(12/02/22/22NkWlxlxDFTkXNlklNN点（2.5）将(2.4)式和(2.5)式代入(2.31)式，并利用kNNkNWW2和X1(k)、X2(k)的隐含周期性可得到：12,...1,0)()()2()()()(2121NkkXWkXNkXkXWkXkXkNkN（2.6）这样，将N点DFT的计算分解为计算两个N/2点离散傅立叶变换X1(k)和X2(k)，再计算(2.6)式。为了将如上分解过程用运算流图表示，以便估计其运算量，观察运算规律，总结编程方法，先介绍一种表示(2.6)式的蝶形图。图2.1蝶形运算图图2.28点DFT一次时域抽取分解运算流图根据图2.2可以求得第一次分解后的运算量。图2.2包括两个N/2点DFT和N/2个蝶形，每个N/2点DFT需要2)2/(N次复数乘法和2/)12/(NN次复数加法运算，每个蝶形只有一次复数乘法运算和两次复数加法运算。所以，总的复数乘法次数为总的复数加法次数为2|)1(222)2(212NNNNNN（2.7）22222)12(2NNNN（2.8）N=8点DIT-FFT的运算流图如图2.3(a)所示。根据WkN/m=WkmN，将图2.3(a)转换成如图2.3(b)所示的标准形式的运算流图图2.3N=8点DIT-FFT的运算流图（2）算法效率由图2.3可见，N=2M时，DIT-FFT运算流图由M级蝶形构成，每级有N/2个蝶形。因此，每级需要N/2次复数乘法运算和N次复数加法运算，M级形共需复数乘法次数CM(2)和复数加法次数CA(2)分别为lbNNMNCM22)2(（2.9）CA(2)=N·M=NlbN（2.10）式中，lbN=log2N。直接计算N点DFT的复数乘法次数为N2,复数加法次数为(N-1)N。当N1时,N2/CM(2)1，所以N越大，DIT-FFT运算效率越高。DIT-FFT算法与DFT所需乘法次数与N的关系曲线如图2.4所示。例如，N=210=1024时，DIT-FFT的运算效率为8.20410210241024)2(22MCNFFTDITDFT的乘法次数的乘法次数（2.11）而当N=211=2048时，37.372112048222)2(22MNMNNCNM（2.12）图2.4DIT-FFT与DFT所需乘法次数比较曲线（3）运算规律rNW的确定第m级运算，一个DIT蝶型运算的的两接点“距离”为12m，所以rNmmmmmrNmmmmWkXkXkXWkXkXkX)2()()2()2()()(1111111（2.13）r的求解：（1）将式（2.13）中第一个节点的标号k表示成L（LN2）位二进制数；（2）把此二进制数乘上mL2，即将L位二进制数左移mL位（注意m是第m级运算），把右边空出的位置补0，此数即为所求r的二进制数。序列的倒序DIT-FFT算法的输入序列的排序看起来似乎很乱，但仔细分析就会发现这种倒序是很有规律的。由于LN2，因此顺序数可用L位二进制数(0121nnnnLL)表示。M次偶奇时域抽选过程如图2.5所示。第一次按最低位0n的0和1将x(n)分解为偶奇两组，第二次又按次低位1n的0、1值分别对偶奇组分解；依次类推，第L次按1Ln位分解，最后所得二进制倒序数如图2.5所示。图2.5形成例序的树状图(N=23)形成倒序J后，将原存储器中存放的输入序列重新按倒序排列。设原输入序列x(n)先按自然顺序存入数组A中。例如，对N=8,A(0)，A(1)，A(2)，…，A(7)中依次存放着x(0),x(1),…,x(7)。对x(n)的重新排序(倒序)规律如图2.6所示。图2.6倒序规律三、实验仪器计算机四、实验要求及内容用所学过的编程语言，自行设计出一个按时间抽取的、输入倒位序、输出顺序的基-2FFT算法程序。五、实验报告（1）简述实验目的及原理；（2）画出程序流程框图；（3）主要给出实验内容的程序（要求程序模块化并加注释）。程序流程框图实验的程序#includeiostream.h#includemath.h#includestring.h#definePI3.1415926classPlural//复数类{floatreal;//实部floatimg;//虚部public:Plural(floatr,floati){real=r;img=i;}Plural(floatr)//通过构造函数重载使用数与复数乘法{real=r;img=0;}Plural(){real=0;img=0;}friendPlural*fft(floatX[],intn);//fft（）；friendPluraloperator*(Pluralp1,Pluralp2);//重载乘法运算符friendPluraloperator-(Pluralp1,Pluralp2);//重载减法运算符friendPluraloperator+(Pluralp1,Pluralp2);//重载加法运算符friendPlural*daoxu(Plural*x[],intn);//倒序Pluraloperator=(Pluralp);//重载赋值符号friendPluralW(intN,inti);//生成旋转因子voidshow();//输出real+imgi};Pluraloperator*(Pluralp1,Pluralp2)//复数乘法{returnPlural(p1.real*p2.real-p1.img*p2.img,p1.real*p2.img+p1.img*p2.real);}Pluraloperator+(Pluralp1,Pluralp2)//复数加法{returnPlural(p1.real+p2.real,p1.img+p2.img);}Pluraloperator-(Pluralp1,Pluralp2)//复数加法{returnPlural(p1.real-p2.real,p1.img-p2.img);}PluralPlural::operator=(Pluralp)//重载赋值符号{real=p.real;img=p.img;return*this;}//********************生成旋转因子***************************PluralW(intN,inti){floatr;floatim;r=(float)cos(2*PI*float(i)/float(N));im=(float)sin((-1)*2*PI*float(i)/float(N));returnPlural(r,im);}//*********************输出函数show（）*************************voidPlural::show()//输出real+imgi{if(real==0){if(img==0)cout0;elsecoutimgi;}else{if(img0)coutrealimgi;elseif(img==0)coutreal;elsecoutreal+imgi;}}/************2的n次方*******************/intmi2(intn){intm=1;for(inti=0;in;i++){m*=2;//m=2^n}returnm;}/*************由二进制数转化为十进制数******************/inttwoten(intxu[],inti){intm=0;for(intj=0;ji;j++){m+=xu[j]*mi2(i-j-1);//m=xu[0]*1+xu[1]*2+xu[2]*4+xu[3]*8+xu[4]*16+...}returnm;}/*************对二进制数倒序加一*****************/voidadd(intxu[],inti){xu[0]++;for(intj=0;ji;j++){if(xu[j]==2){xu[j]=0;xu[j+1]++;//例如xu[]=0000,1000,0100,1100，0010......}}}/*************倒序**********************************/voiddaoxu(Pluralx[],intn){floatm=float(n);inti=0;for(i=0;m1;i++)//得到n是2的多少次方{m/=2;//m=log(n)/log(2)}intM=mi2(i);int*xu=newint[i];//定义一个长度为i的数组，当做2进制数for(intj=0;ji;j++){xu[j]=0;}Plural*jia=newPlural[M];//定义一个长度为M的数组jia，以保存倒序for(j=0;jM;j++){intmm;mm=twoten(xu,i);//得到数组xu[]所对应的十进制数jia[j]=x[mm];add(xu,i);//二进制数组最左端加1}for(j=0;jn;j++)//撤销动态生成的数组{x[j]=jia[j];}if(jia)delete[]jia;}//***************快速傅里叶变换FFT()*************************Plural*fft(PluralX[],intm){intn=1;for(;mmi2(n);n++);//得到m等于n次方Plural*A=newPlural[m];for(inti=0;in;i++)//共有n级{for(intj=0;jmi2(n-1-i);j++)//第i级有2^(n-i)族{for(intk=0;kmi2(i);k++)//每族有2^(m-1)个蝶形运算{intm,n;m=k+j*mi2(i+1);n=k+mi2(i)+j*mi2(i+1);Plurala1,a2;if(i==0){a1=X[m]+W(mi2(i+1),k)*X[n];//蝶形运算a2=X[m]-W(mi2(i+1),k)*X[n];//蝶形运算A[m]=a1;A[n]=a2;}else{a1=A[m]+W(mi2(i+1),k)*A[n];//蝶形运算a2=A[m]-W(mi2(i+1),k)*A[n];//蝶形运算A[m]=a1;A[n]=a2;}}}}for(i=0;im;i++)//撤销动态生成的数组{X[i]=A[i