离散傅里叶变换及其快速算法解说.docx

上机实验一：离散傅里叶变换及其快速算法 一、设计目的 通过编写程序，深入理解快速傅里叶变换算法（FFT）的含义，完成FFT算法的软件实现。 二、设计任务 利用时间抽取算法，编写基2的快速傅立叶变换（FFT）程序，并在FFT程序基础上编写快速傅立叶反变换（IFFT）。 三、设计要求 1、FFT和IFFT子程序相对独立、具有一般性，并加详细注释； 2、验证例5-4，并能得到正确结果； 四、设计条件 C语言 五、编程规则 1）程序输入元素的数目为2的整数次幂，即N为2M幂，整个运算需要M级蝶形运算。 2）输入序列按二进制码位倒置排列，输出序列按自然顺序排列。 3）输出数据占用输入数据的存储单元。 4）每一级含N/2个基本蝶形运算。第L级中有N/2L个群，群与群间隔为2L。同一级中各个群的系数W分布相同。处于第L级的群的系数是2L. 8）对于第L级的蝶形运算，输入数据的间隔为2L-1。 根据上述要求，设计程序源代码如下所示： #includemath.h #includeiostream.h #includeiomanip.h #define swap(a,b) {float T; T=(a);(a)=b;(b)=T;} void fft(float A[],float B[],unsigned M) //蝶形运算程序，A存实部，B存虚部，M是级数 { unsigned long N,I,J,K,L,LE,LE1,P,Q,R; float Wr,Wi,Wlr,Wli,WTr,WTi,theta,Tr,Ti; N=1M;//N=1M表示N=2^M J=0; for (I=0;IN-1;I++)//for循环负责码位倒置存储 { if(JI) { swap(A[I],A[J]);swap(B[I],B[J]); } K=N1;// K=N1表示K=N/2 while (K=2J=K)//while循环表示须向次高位进一位 { J-=K;K=1;//K=1表示K=K/2 } J+=K; } for(L=1;L=M;L++)//for循环为M级FFT运算,外层循环由级数L控制，执行M次 { LE=1L;// LE=1L表示2^L,是群间隔 LE1=LE/2; //每个群的系数W数目 Wr=1.0; Wi=0.0; theta=(-1)*3.1415926536/LE1; Wlr=cos(theta); Wli=sin(theta); for(R=0;RLE1;R++)//中层循环由群系数控制 { for(P=R;PN-1;P+=LE)//R是群系数的编号，P、Q是基本蝶形运算两个输入数据在数组中的编号，循环每次完成同一个系数 的蝶形运算 { Q=P+LE1; Tr=Wr*A[Q]-Wi*B[Q];Ti=Wr*B[Q]+Wi*A[Q]; A[Q]=A[P]-Tr; B[Q]=B[P]-Ti; A[P]+=Tr; B[P]+=Ti;} WTr=Wr; WTi=Wi; Wr=WTr*Wlr-WTi*Wli; Wi=WTr*Wli+WTi*Wlr; } } return; } } void main()//主程序 { int i,M,N,lb; cout请输入转换类别(FFT，请输入1； IFFT，请输入0)endl;//确定转换类别 cinlb; cout请输入序列长度Nendl; cinN; float *A= new float[N]; float *B= new float[N]; M=log(N)/log(2); cout请输入序列的实部endl;//输入序列实部 for(i=0;iN;i++) { cinA[i]; } cout请输入序列的虚部endl;//输入序列虚部 for(i=0;iN;i++) { cinB[i]; } cout setiosflags(ios::fixed);//输出格式控制 cout您输入的序列为endl; cout setiosflags(ios::fixed); for(i=0;iN;i++) coutA[i]+jB[i]endl; coutendl; if(lb==0) { for(i=0;iN;i++) { B[i]=B[i]*(-1); } fft(A,B,M);for(i=0;iN;i++) { B[i]=B[i]*(-1); } for(i=0;iN;i++) { A[i]=A[i]/N;B[i]=B[i]/N; } } if(lb==1)fft(A,B,M); cout转换后的序列为endl;//输出序列