用matlab实现DFTFFT

1、用matlab实现DFTFFT目录实验目的2实验内容21 .用MATLAB实现DFT22 .用MATLAB实现FFT,分析有限离散序列的FFT33 .通过分别计算时间，得出DFT与FFT的算法差异7实验原理71 .离散傅里叶变换的快速算法FFT72 .FFT提高运算速度的原理83 .理论分析DFT与FFT算法差异10实验步骤11实弱结果11实验分析21实验结论25实验体会26实验目的1 .通过研究DFT,FFT性质，用语言实现DFT.FFT。不使用MATLAB现有的FFT函数,自己编写具体算法。2 .掌握FFT基2时间抽选法，理解其提高减少乘法运算次数提高运算速度的原理。3 .设计实验,得出D

2、FT和FFT算法差异的证明,如复杂度等(精度、不同长度的序列等)。实验内容1.用MATLAB实现DFTN点序列x(n)的DFT为:N-1X(k)=£x(n)W世n=00k<N-1(n-D(nt)VVNW1WV¥NW；CN-nDFT的矩阵为:X(0)X(l)X(2)X(N-1)根据DFT公式与矩阵展开，通过MATLAB实现DFT：2 -N=16;%lengthofsequence3 n=0:N-l;%tiaesample4 -xn=cos(pi*n/6);%generatesequence5 k=0:N-1;%frequencysample6 -WN=exp(-j*2*

3、pi/N)：7 -nk二n'*k:8 -WNnk=IN.*nk;9 -Xk=xn*WNnk;'computeDFT10 %plotsequenceanditsDFI11 -figure(l)12 -stem(xn)13 -figure(2)14 -stem(abs(Xk)2.用Matlab实现FFT编程思想及程序框图：原位计算因为DIT-FFT与DIF-FFT的算法类似，这里我们以DIT-FFT为例。N=2”点的FFT共进行M级运算，且每一级都由N/2个蝶形运算组成，后一级的节点数据由前一级同处一条水平线位置的节点数据产生，所以我们同样可以将后一级的节点数据储存到前一级的节点中

4、，这样的方法叫做原位计算，它大大节省了内存资源，降低了成本，简化了运算。序列的倒序无论是进行DIT-FFT还是DIF-FFT都需要进彳序，包括输入倒序与输出倒序，以一定的方式将数组进行重新排列。倒序的方法：首先由于N=2：我们就可以用M位二进制数来表示节点的顺序，并且按照奇偶时域抽取。然后，如图7所示，第一次按最低位nO的0、1值分解为奇偶组，第二次按次低位n1的0、1值分解为奇偶组，以此类推。最后，所得二进制数所对应的十进制数即为序列倒序后产生的序列。n20b2blbO000(b2blbO)图1序列倒序过程倒序的MATLAB方法：用雷德算法可以对输入信号序列进行倒序重排，流程图如下所示:(开

5、始)I1、LH<N/2J<LHNl<N.2I蝴蝶因子的变化规律在DIT-FFT中，每一级都由N/2个蝶形运算构成,每个蝶形运算包含一个蝴蝶因子,每一级的蝶形因子又有一定的变化规律：设L表示自左而右的运算级次(L=1,2,3,M),序数R,次数K。每个蝶形运算的两个输入量相距B=2XL-1)个点。假设N=8,则M=3,L=1时,B=1f这时有:S二N/2,rr,K=1:N/2则有P=(K-1)*1f所以Wj5二W<，J=0.1,2,3L二2时,B二2,S=N/4,R=0:N/2:N-1,K=1:N/4则有P=(KT)*2,所以Wbw/2,J=o.2L=3时,B=4,S=N

6、/8,R=0:N/4:N-1,K=1:N/8则有P=(K-1)*4,所以二w0/4，J=0所以对于一般情况N=21第L级的旋转因子为P=(KT)*B。从而得到DIT-FFT的Matlab流程图：23456一/891011121314151617181920212223242S26272829303132333435363738394041424344DIT-FFT的MATLAB实现为(以对x(n)=cos(-)进行16点的变化为例)：0clear;clc;- N=16：%离散数组长度- n=0:N-l：%时域采样- xn=cos(pi*n/6):%产生离散扳组- M=nextpow2(N)嘱计

7、算常里M；NEXIP02找到最近以2为底N的对数- A=xn,zeros(l,N-length(xn)(n)输入到A，如果输入的数组长度不是】的整刻次方，填零补位- dispC输入到各存储单元的数据：')，disp(A);- 1=0:%给倒序数喊初值- QforI=O:N-1%这个for循环的结果使原来数组次序改变，如数组123141516经过此循环以后变为%以下序列19513311715210614412816,咨询FF：制序法则- ifI<J：- I=A(I4-1)；A(I4-1)=A(J1):A(J-*-li=T；- K=N/2;- whileJ>=K:- J=J-K

8、;K=K/2：- end- J=J+K:- end- dispC倒序后各存储单元的数裾：disp(A)：%将变换次序以后的数组域值给A%以下是蝴蝶算法.- WN=exp(-j*2*pi/N)：%计算常里- forL=1:M%最外层循环是蝴蝶算法的总行数- dispC运算级次：')，disp(L)：- B=2XL-1)：%乂为第r行的组数(比如第一行有2,(L-r)组，第二行有二(L-2)组最后一行只有一组)- forR=O:B-1:- P=2*(M-L)*R：%以下两个for循环处于潜逃循环的最内层，在同一组内，前一半的计算公式和后一半的计算公式不一样- 巾forK=R:2aL:N-2

9、：- T=A(K+D+A(K+B+1)*WNAP;- A(K+B+1)=A(K+1)-A(K+B+1)*WN*P：- A(K+1)=T;- end- end- disp('本级运算后存储单元的数据：),dispi：A):- end- dispC输出各存储单元的数据：')，Xk=A,%输出运算结果%telapsed=toe(tstart);3.通过分别计算时间，得出DFT与FFT的算法差异:法一：对N二512,1024,2048和4096点的离散时间信号x(n),用Matlab语言编程分别以DFT和FFT计算N个频率样值X(k),比较两者所用时间的大小。设x(n)=cos(?)。

10、0利用TIC,TOC函数，分别跟踪变换时间。法二：利用clock返回当前日期向量的时间，通过etime0函数返回走过的日期向量的时间。用MATLAB实现如下：Nmax=2048;2-fft-tiMe=zeros(l,Nrax)：4-5-6一7-8一forn=l:1:N>axx=rand(l,n);%Uniformlydistributedrandonnumbers.t=clock;fft(x)：fft-tiMe(n)=etime(clock,t);号tineend9-n=l:l:Njaax;top=max工me)：10-plot(n,.3;axis(0?2500,0,top);11-xl

11、abelCN");ylabelClimeinSec.');titleCFFIexecutiontimes")12-text(2100»top/。(N*N)13-text(2100,top*3/4,*o(N*N*3/4?)14-text(2100,top/2,/o(N*N/2),)15-text(2100,top/3,*o(N*N/3),)16-text(2100,top/4,o,N*N/4),)17-text(2100,ime),?o,N*logN)/)实验原理1.离散傅里叶变换的快速算法FFTN点序列x(n)的DFT为：X(k)=x(n)Wj!)k0&l

12、t;k<N-1.九二0由于系数w目=e附位是一个周期函数:W；(Nk)=wk(N-n)=w-nk(1)(2)且是对称的:W+N/2=_wnk(3)快速傅里叶变换算法正是基于这样的基本思想而发展起来的，它的算法基本可以分称两大类：时间抽取法(DIT-FFT)和频率抽取(DIF-FFT)o由于DIF-FFT算法思想基本一致，只是划分方式略有差异，所以这里以DIT-FFT算法为例进行说明。当N是2的整数次方时，称为基2的FFT算法。首先将序列x(n)分解为两组，偶数项为一组，奇数项为一组：r=0.L.,y-1(4)4:x(2r)=xi(r)x(2r+1)=x2(r)将X(r)和X2(r)分别进

13、行N/2点的DFT得&®和X2电,且:=0,1,-1X(k)=Xi(k)+WX2(k)X(+k)=X1(k)-W1!jX2(k)重复这一过程，可得到X(n)的FFTo2 .FFT提高运算速度的原理FFT算法将长序列的DFT分解为短序列的DFToN点的DFT先分解为2个N/2点的DFT,每个N/2点的DFT又分解为N/4点的DFT,以此类推。最小变换的点数即所谓的“基数”，这个基数是2点的DFT(又叫蝴蝶因子)。如公式(5)所示，对于X(k)的计算可以分为两部分，因为在上下两个式子中都出现了Xi(k)与X2(k),因此只需要一次复数乘法与两次复数加法即可分别得到N/2点的DFT

14、,从而得到总的N点的DFT。如图2所示。下图为DIT-FFT与DIF-FFT蝴蝶因子的三种形式:DITDIF图3DIT-FFT和DIF-FFT的蝴蝶因子原始形式；（2）简化形式；（3）优化形式下面以8点的DIT-FFT为例:图48点的DIT-FFT对于8点的FFT,我们需要M=log28二3阶运算，每一阶有四个蝴蝶因子，在每个蝴蝶因子中需要1次复数乘法与两次复数加法，因为每个蝴蝶因子都一样，所以FFT通过重复运算大大简化了算法与复杂度。3 .理论分析DFT与FFT算法差异DFTN-1X(k)=2x(n)Wj5k0<k<N-1n=0DFT算法复杂度的计算：复数乘法复数加法一次x(k)

15、NM-1X(k)(N点DFT)N2N(N-1)对于N=512、1024和8192的DFT,分别需要复数乘法（次）N2=5122=2,8=262144N2=10242=220=1048576N2=81922=226=67108864这是一个巨大的数字，很难在实际中应用。X(k)=Xi(k)+WX2(k)FFT=Xi(k)怵Xz(k)FFT算法复杂度的计算：>对于N=211点的FFT,需要M=log?N阶的重复运算;> 每一阶有N/2个蝴蝶因子；> 每一个蝴蝶因子需要1此复数乘与2次复数加。复数乘法复数加法N点FFTN510g2NNlog2N对于N二512、1024和8192的D

16、FT,分别需要复数乘法（次）ylog2N二5Elog2512=2304ylog2N=竽log21024=5120ylog2N=等晦8192二53248相比于DFT,FFT的运算复杂度大大降低了。实验步骤1 .若x(n)=cos(詈)是一个N二16的有限序列，利用MATLAB计算它的DFT并画出图形。若N不等于2的整数次方，假设N=12,利用MATLAB计算它的DFT并画出图形。2 .若x(n)=cos(詈)是一个N二16的有限序列，利用MATLAB计算它的FFT并画出图形。若N不等于2的整数次方，假设N=12,利用MATLAB计算它的FFT并画出图形。3 .根据上述两个步骤的实验结果，比较DF

17、T与FFT在算法与结果上的相同与差异。4 .通过对N二8的有限序列x(n)=2,4,T.0,-2.-43J进行FFT变换，在MATLAB中加入disp()函数以及序列存储单元，跟踪并分别展示：输入到各存储单元的数据倒序后各存储单元的数据各运算级次对应的该级运算后存储单元的数据输出各存储单元的数据通过数据分析FFT具体算法，更深入地理解FFTo5 .通过tic(),toe()函数，计算不同长度序列的DFT与FFT所用时间，得到二者算法差异。令N为不同长度的序列，分明算出DFT和FFT的耗时并进行比较： N二512点时 N二1024点时 N二2048点时 N二4096点时通过理论分析，找到DFT与

18、FFT的算法差异。实验结果1 .有限序列x(n)=cos(?)0>N=16利用MATLAB计算它的DFT：Xk=Columns1through52 .36604.4923+4.961Oi-0.6124-3.3371i0.3021-1.4336i0.5000-0.8660iColumns6through100.5766-0.5549i0.6124-0.3371i0.6290-0.1604i0.6340-0.0000i0.6290+0.1604iColumns11through150.6124+0.3371i0.5766+0.5549i0.5000+0.8660i0.3021+1.4336i

19、-0.6124+3.3371iCoIumn164.4923-4.961Oi输入有限序列x(n)：stem(xn)图5.116点DFT输入序列x(n)序列的DFT：stem(abs(Xk)7图5.216点输入序列x(n)的DFT»N二12利用MATLAB计算它的DFT：Xk=Columns1through4-0.00006.0000+0.0000i- 0.0000-0.0000i- 0.0000-0.0000iCoIumns5through8- 0.0000-0.0000i- 0.0000-0.0000i- 0.0000-0.0000i0.0000-0.0000iColumns9thr

20、ough120.0000-0.0000i0.0000-0.0000i0.0000-0.0000i6.0000+0.0000i输入有限序列x(n)：stem(xn)图6.112点DFT输入序列x(n)序列的DFT：stem(abs(Xk)0L0120©0©G0G246810图6.212点输入序列x(n)的DFT2.有限序列x(n)=cos(?)0>N=16利用MATLAB计算它的FFT：Xk二CoIumns12.3660Columns3-0.6124-CoIumns50.5000-Columns70.6124-CoIumns90.6340throughthrough3.

21、3371ithrough0.8660ithrough0.3371ithrough24.4923+4.961Oi40.3021-1.4336i60.5766-0.5549i80.6290-0.1604i100.6290+0.1604iColumns11through0.6124+0.3371iCoIumns13through0.5000+0.8660iColumns15through-0.6124+3.3371i120.5766+0.5549i140.3021+1.4336i164.4923-4.961Oi有限序列x(n)：stem(xn)0246810121416图7.116点FFT输入序列x

22、(n)序列的DFT：stem(abs(Xk)图7.216点输入序列x(n)的FFT>N=12利用MATLAB计算它的DFT：ans=Columns1through5-0.00006.0000-0.0000i-0.0000+0.0C00i-0.0000-0.0000i0.0000-0.0000iColumns6through100.0000+0.0000i0.00000.0000-0.0000i0.0000+0.0000i-0.000C+0.0000iCoIumns11through12-0.0000-0.0000i6.0000+0.0000i有限序列x(n)：stem(xn)图7.116

23、点FFT输入序列x(n)序列的DFT：stem(abs(Xk)6i91111>5-4-3-2-1-o1geeo$©$o-ee024681012图7.216点输入序列X(n)的FFT3 .根据上述两个步骤的实验结果，比较DFT与FFT在算法与结果上的相同与差异。> DFT与FFT在结果上相同。Xk相同，产生的图像也相同。> FFT是DFT的一种快速算法，对于16点的FFT,我们需要M=log?16=4阶运算,每一阶有八个蝴蝶因子，在每个蝴蝶因子中需要1次复数乘法与两次复数加法,因为每个蝴蝶因子都一样，所以FFT通过重复运算大大简化了算法与复杂度。> 对于16点的

24、DFT,需要复数乘法162=256次,FFT需要复数乘法510g2N=32次。> 对于12点的FFT,因为N不是2的整数次嘉，所以补4位0,得到N=16的FFT,并计算结果，对于12点的DFT与FFT,结果同样相同4 .通过对N二8的有限序列乂(向二2,4,-1,0,-2,-431进行FFT变换，在MATLAB中加入disp()函数以及序列存储单元，跟踪并分别展示：>输入到各存储单元的数据24-10-2-431>倒序后各存储单元的数据2-2-134-401>各运算级次对应的该级运算后存储单元的数据L=1042-4081-1L=2Columns1through224+4i

25、Columns3through4-24-41Columns5through618+1iColumns7through8-18-1iL=3Columns1through23.000010.3640-0.9497iColumns3through4-2.0000+1.0000i-2.3640-8.9497iColumns5through61.0000-2.3640+8.9497iColumns7through8-2.00001.0000i10.3640+0.9497i>输出各存储单元的数据Xk=Columns1through23.000010.3640-0.9497iColumns3throu

26、gh4-2.0000+1.0000i-2.3640-8.9497iCoIumnsthrough1.0000-2.3640+8.9497iColumnsthrough-2.0000-1.0000i10.3640+0.9497i有限序列X(n)：stem(xn)12345678图8.1x(n)输入序列序列的DFT：stem(abs(Xk)12图8.2序列x(n)的FFTXk模值表达pIot(Xk)图8.3序列x(n)的FFT复数表达5.利用tic函数，令N为不同长度的序列，分明算出DFT和FFT的耗时并进行比较,设x(n)=cos(毛)，N=16：N二512点时：dftcost_time=0.30

27、58fft_cost_time=0.0658 NF024点时：dft_cost_time=fft_cost_time= N二2048点时：dftcost_time=fft_cost_time二 N二4096点时：dft_cost_time=fftcosttime=0.93630.14463.73180.348916.75500.7791由上面数据表明，同样长度的信号，DFT耗时要比FFT耗时要少，表明FFT算法的有效性。实验分析1.DFT与FFT算法不同，FFT是通过DFT运算中存在对称性和周期性而做的化简，因而大大提高了变化速度。N点序列x(n)的DFT为：N-1X(k)=2x(n)W群n=

28、00<k<N-1由于系数=e-是一个周期函数:(N-k)且实对称的:W+N/2=_wnk快速傅里叶变换算法正是基于这样的基本思想而发展起来的，以DIT-FFT算法为例进行说明。首先将序列x(n)分解为两组，偶数项为一组，奇数项为一组：(x(2r)=xi(r)n(X(2r+1)=x2(r)ru，L，5(4)将Xi(r)和X2(r)分别进行N/2点的DFT得XI(k)和X?®,且:X(k)=Xi(k)+WX2(k)nxg+k)=x1(k)-wx2(k)r=0，l，门可以发现对于公式(5)的前k项与2/N+k项都包含相同的DFTXi(k)和X2(k),利用这一性质，我们可以简化

29、DFT算法。重复这一过程，我们可以将xi(r)和X2(。分别继续分成偶数项一组，奇数项一组，如图8所示：x(2r)=xI(r)<x(n)x(2r+l>x2(r)<X(21f(叫.,x1(2I+1)=x4(I)-bx,(21)=xf?(l)j工(21+A/(M*")=/(/)/=%(21+1)=5(04X】(A)=Xj(A)+“：相(外*(4+2=1(左)一73*曾马因=$(6+。，工“)W».Y")一L4V2=0.1-14V£=0-1从而得到如下算法:T5图9根据步骤4,通过对N=8的有限序列x(n闫2,4,-1。-2,-4,3,1进行

30、FFT变换,因为1(l-7)/V27-(l+7)/V2有理论计算结果:9-7j,V2-LU3640-0.9497J与如下MATLAB计算结果相同：输入到各存储单元的数据24-10-2-431倒序后各存储单元的数据2-2-134-401各运算级次对应的该级运算后存储单元的数据L=1042-4081Columns1through224+4iColumns3through4-24-41Columns5through618+1iColumns7through8-18-1iL=3Columns1through23.000010.3640-0.9497iColumns3through-2.0000+1.0

31、000iColumns5through1.0000Columns7through-2.0000-1.0000i4-2.3640-8.9497i6-2.3640+8.9497i810.3640+0.9497i因而，我们通过理论分析，精练了DFT算法，得到了FFT算法，它们两本质上是一样的,FFT是通过DFT运算中存在对称性和周期性而做的化简。2.DFT与FFT算法复杂度比较DFTN-10<k<N-1X(k)=Wn=0DFT算法复杂度的计算：复数乘法复数加法一次X(k)NN-1N次X(k)(N点DFT)IfN(N-1)对于N=512、1024和8192的DFT,分别需要复数乘法(次)N2=5122=2'8=262144N2=10242=220=1048576N2=81922=226=67108864这是一个巨大的数字，很难在实际中应用。X(k)=X1(k)+WX2(k)xG+k)=X1(k)-W/jX2(k)FFT算法复杂度的计算：> 对于N=2”点的FFT,需要M=log2N阶的重复运算;> 每一阶