1、大神大學課程設計說明書題目：DFT與FFT計算速度比較分析系 別：工業自動化儀表年級專業：22級儀表1班學 號：2學生：儀表小子指導教師：林大神大神教師職稱：大神 大神大神大學課程設計（論文）任務書院（系）：大神工程學院基層教學單位：自動化儀表系學號2學生儀表小子專業（班級）22級儀表1班設計題目DFT與FFT計算速度比較分析設 計 技 術 參 數用MATLAB現DFT及FFT對任意長度的序列進行傅里葉變換DFT與FFT的運算時間比較設 計 要 求利用Matlab或者C語言設計DFT和FFT程序，比較兩種頻譜分析方法的 計算速度，并與理論值進行比較。工 作 量先對兩種算法進行介紹， 包括推導過

2、程及運算性質，然后用MATLA限現兩種算法，再分別對兩種算法進行運算時間對比，并分析時間長短的原因。工 作 計 劃A周第二周即-接受任務并查閱資料周二到周五上午學習相關知識下午編寫程序上機調試程序到周四上午學習相關知識下午編寫程序上機調試程序編寫任務書參 考 資 料1 .平、王娜、林洪斌等，信號處理原理及應用。機械工業，2008.102 .王宏，MATLAB 6.5 及其在信號處理中的應用。清華大學，2004.103 . Sanjit K.Mitra著洪、余翔宇等譯，數字信號處理實驗指導書。電子工業2005.1指導教師簽字林大神J大神？基層教學單位主任簽字大仙說明：此表一式四份，學生、指導教師

3、、基層教學單位、系部各一份。2011年7月13 日大神工程學院課程設計評審意見表指導教師評語： 認真工作態度 較認真不認真平時成績：理論分析止確完善較為合理較不指導教師簽字：軟件設計完善 合理較不1月日圖向及其它成績：答辯小組評語：清晰基本掌握原理了解不清楚優化設計止確基本止確 不止確答辯成績：組長簽字:1月日課程設計綜合成績：答辯小組成員簽字：年月日摘要時域分析方法和頻域分析方法是信號和系統的分析的兩種方法，本文介紹離散信號和系統的頻域分析方法，它和連續信號和系統的頻域分析方法有所不同，但也有相似之處。本說明書主要是在介紹兩種用于信號處理的傅里葉變換算法一一DFT （離散傅里葉變換）和FFT

4、 （快速傅里葉變換），分別介紹了這兩種運算的推導過程，并且對這兩種變換作了簡要的介紹，分析了各自的性質。然后通過MATLAB 分別實現了這兩種傅里葉變換，并對這兩種變換進行了運算時間的比較一一分別對同一函數進行 DFT和FFT計算兩者的運行時間，并作圖比較。本說明書的程序部分都是在MATLAB 環境下進行的運算。MATLAB 是矩陣實驗室（Matrix Laboratory）的簡稱，是美國MathWorks公司出品的商業數學軟件，用于算法開發、數據可視化、數據分析以及數值計算的高級技術計算語言和交互式環境，主要包括MATLAB 和 Simulink 兩大部分。MATLAB 的基本數據單位是矩陣

5、，它的指令表達式與數學、工程中常用的形式十分相似，故用MATLAB 來解算問題要比用C， FORTRAN 等語言完成相同的事情簡捷得多。在新的版本中也加入了對C， FORTRAN， C+ ， JAVA的支持，可以直接調用,用戶也可以將自己編寫的實用程序導入到MATLAB函數庫中方便自己以后調用。本文介紹了DFT 與 FFT 的原理與Matlab 實現程序，以及 DFT 與 FFT 的關鍵詞：計算速度的比較。并用guide函數親自編寫了一個界面。DFT、FFT、Matlab、運算速度、guide目錄摘 要 1第一章 DFT 原理與 Matlab 實現 31.1 DFT 的原理 31.2 DFT

6、的Matlab 實現 4第二章 FFT 的原理與Matlab 實現 62.1 FFT 的原理 62.1.1 FFT 的基本思想 62.1.2 基 2 FFT 算法 72.2 FFT 的Matlab 實現 9第三章 DFT 與FFT 計算速度比較分析 123.1 FFT 與直接計算DFT 的比較 123.2 FFT 與 DFT 運算時間Matlab 程序 133.2.1 隨機序列的DFT 計算時間程序133.2.2 分析兩者運算時間的差異：16第四章 心得體會17參考文獻：18第一章DFT原理與Matlab實現1.1 DFT的原理傅里葉變換就是在以時間為自變量的蓿號”與以頻率為自變量的 頻譜”函

7、數之間的某種變換關系。隨時間自變量形式的不同，其傅里葉變換的形式也有不同：周期序列的離散傅里葉級數（DFS）和非周期序列的傅里葉變換 （DTFT）,其表示式分別為：泌 k DFS %n1j nk%ne Nn 0X(ej ) DTFT x(n) x(n)en(1.1.1)(1.1.2)在實際工作中，當用數字計算機對信號進行頻譜分析時，要求信號必須以有限長度 的離散值作為輸入，而計算所得的頻譜值自然也是有限、離散的。上述兩種形式的傅里 葉變換中，DFS變換滿足時、頻域自變量的離散化，但其時間變量和頻率變量又同時具 有周期性；DTFT變換滿足時間自變量的有限長度（非周期能量有限信號），但其頻率變量為

8、連續形式。可見，這兩種變換都難以實際應用。考慮到 DFS變換的時、頻域形式雖是 周期序列，但每個周期卻只有 N個獨立的復值，知道其一個周期的容即可得到其它的容。 因此，若從DFS變換的時、頻域各取出一個周期，即可構造出時間和頻率自變量皆為離散、有限長度的傅氏變換，這就是離散傅里葉變換（DFT）的引出思想，下面進行具體推導。設xn是一個長度為M的有限長序列，由周期序列與有限長序列的本質聯系，可以N（ N M ）為周期將x n展開為無重疊的周期序列，即周期延拓為x n rNr再利用式（1.1.1）對%nXk的主值序列（k 0,1,.,N進彳T DFS變換，得到周期離散的頻譜(1.1.3)X%k (

9、k (,萬，取1），代入DFS反變換公式（4-3b）,即分析可見：2j -n-nke N (n 0,1,.,N一個周期的值(kN 1i 2 kn%n IDFS )%k制k eN k 0(1.1.4)在 DFS正變換中，只要把一個周期的xn乘以對應的1),即可得任意k下的X k ；同理，在DFS反變換中，僅用 X k的0,1,.,N 1),即可得到任意n下的x n 。如果同時限制(1.1.1)式中的n和(1.1.4)式中的k,使其都只在區間取值,就得到了一個周期的X n和一個周期的X k間的對應關系1N-k1x01X0WNknknW N,2j-式中，WN e " , N為DFT變換區間

10、的長度,(1.1.5)(1.1.6)上兩式即稱為有限長序列的離散傅里葉變換對。(1.1.5)式稱為離散傅里葉變換，簡稱DFT; (1.1.6)式稱為離散傅里葉逆變換(Inverse Discrete Fourier Transform), 簡稱 IDFT。1.2 DFT的Matlab實現程序DFT函數:functionxk=dft(xn,N)%計算離散傅里葉變換%Xk=dft(xn,N)%Xk=在 0<=k<=N-1 間的%xn=N點有限長度序列%N=DFT的長度n=0:1:N-1;k=0:1:N-1;WN=exp(-1i*2*pi/N);nk=n'*k;WNnk=WN.A

11、nk;q=xn*WNnk;DFT系數數組%n的行向量%k的行向量%Wn因子%產生一個含nk值白N乘N維矩陣%DFT矩陣%DFT系數的行向量對一個單位抽樣序列的DFT變換(N=64):clear all;N=64;x=zeros(1,N);x(1)=1;xn=0:N-1;subplot(121)stem(xn,x);axis(-1 33 0 1.1);XK=dft(xn,N);subplot(122);stem(abs(XK);DFT Matlab處理結果：Ficiire- I=_ 二區 I第二章 FFT 的原理與Matlab 實現2.1 FFT 的原理DFT 在數字信號處理中有很重要的作用，如

12、頻譜分析、線性卷積等，但由于直接計算 DFT 的計算量與變換區間長度N 的平方成正比，當N 較大時，計算量太大，所以在快速傅里葉變換(Fast Fourier Transform,簡稱FFT)出現前，直接用 DFT進行譜分析和信 號的實時處理是不切實際的。直到1965年庫利(J. W Cooley)和圖基(J. W. Tukey)提出了 DFT 運算的一種快速方法以后，情況才發生了根本的變化。多年來，人們不斷改進和完善，形成了一系列新型FFT 算法，如基2FFT 算法、分裂基FFT 算法、 N 為復合數的 FFT 算法等。必須強調指出：FFT 并不是與DFT 不同的另外一種變換，而是為減少DF

13、T 計算次數的一種快速算法。為了解FFT 高效算法的重要及實現思路，先介紹DFT 的運算特點，再具體討論高效算法。2.1.1 FFT 的基本思想從哪些方面能改進DFT 的運算以減少運算工作量呢? 如前所述，DFT 的運算量是與2N 成正比的。顯然，如果一個大點數N 的 DFT 能分解為若干小點數DFT 的組合，則可達到減少運算工作量的效果。充分利用系數WNnk 的下列固有周期性和對稱性，使DFT運算中的有些項可以合并，從而使長序列的DFT分解為更小點數的 DFT,提高運算效率。nkk N nkn knWNnk的對稱性WNWN WNnknk knNkN nWNnk 的周期性WNnkWNknNWN

14、kN n快速傅里葉變換算是基于上述基本思想而發展起來的。下面介紹最常用的基2 FFT算法 ( N 2m ，庫利一圖基算法)。2.1.2 基2 FFT算法基2 FFT算法主要包括兩類：按時間抽取(Decimation-in-time ,簡稱DIT)法和按頻率且N 2m，其中m為整數。如果不抽取(Decmation-in-Frequence ,簡稱 DIF)法，本文只介紹 DIT算法。設xn是列長為N n 0，1，N 1的輸入序列,滿足這個條件，可以人為地加入若干零點來達到。將x n按n的奇偶分成兩個子序列則(1.1.5)式可化為X k DFTN 2 1xr 0N 2 1Xir 0x 2rx1 r

15、x 2r 1 x2 rN 1x n x nn偶數N 2 1 2rk2r Wnxr 0N2 rk kr Wn WnN，r 0,1, 1(2.1.1)N 1W:x n wN1kn奇數2r 1 k2r 1 Wn2 12 rkX2 r Wn0.2 j .2N0 j 2由于WN e N e 2WN ,故上式又可表示為2N 2 1X kX1r WNkr 022N 2 1X2 r WNkX1 kW：X2 kr 0-2k 0,1,N 1(2.1.2)其中X1 k和X2 k分別是x1 r及x2 r的N/2點的DFTN 21(2.1.3)(2.1.4)X 1 kx1 r WNrk2DFT x1 r r 0N 2

16、1X2 kx2 r WNrk2DFT x2 r r 0(2.1.2)式表明了一個N點的DFT被分解為兩個N'，2點的DFT。但是這里有一個問題，即x1 r , x2 r的列長為N2 ,它們的DFT X1 k , X2 k的點數也是N2 ,即, N ,k 0,1,12 ，而Xk卻有N個點，所以按(2.1.2)式計算得到的只是X k(k 0，1， N 1)的前一半項數的結果，要用Xi kX2 k來表達全部的X k值還必須應用W系數的周期性，即W N*2這樣可得NX1萬同理可得另外又考慮到W；的對稱性因此，將上述公式代入X”X10rWnX2WnN ；2N22.N - k221X1 r WNk

17、2X10X2Nwn2WNkwNk(2.1.2)中，X k又可表達為X1 kW：X2 kk0,1,L N 12(2.1.5)X NX1萬N .5 k NWn2X2 k2X1kWNX2k k。儲1(2.1.6)由上分析可見，只要求出N N Y 0 12 區間各個整數k值所對應的X1 k和X2 k值，即可求出0N 1區間的全部Xk值，這一點恰恰是FFT能大量節省計算的關鍵所在。WkX2 k式(2.1.5)、(2.1.6)的運算可用圖2.1.1的信號流圖符號表示，根據其形狀稱之為蝶形運算符號。依此類推得8點的DIT-FFT的運算流圖如下:x(0)x(1)2.2 FFT 的 Matlab 實現程序FFT

18、 函數：function xn=myfft(x)N=length(x);M=log2(N);xtmp=zeros(1,N);value=zeros(1,M);for i=0:N-1repr=i;for t=1:1:Mrepr=bitshift(i,1-t);value(t)=bitand(repr,1);endpos=0;for k=1:1:Mpos=pos+value(k)*2A(M-k);endxtmp(pos+1)=x(i+1);end for i=1:Mdeepth=2A(i-1);width=2A(M-i);for t=1:2Ai：Nfor k=1:deepthtmp=xtmp(t+

19、k-1);wn=width*(k-1);xtmp(t+k-1)=tmp+exp(-j*2*pi*wn/N)*xtmp(t+k+deepth-1);xtmp(t+k+deepth-1)=tmp-exp(-j*2*pi*wn/N)*xtmp(t+k+deepth-1); endendendxn=xtmp;對一個單位抽樣序列的FFT 變換(N=64)：clear all;N=64;x=zeros(1,N);x(1)=1;xn=0:N-1;subplot(121),stem(xn,x);axis(-1 33 0 1.1);XK=fft(xn);subplot(122);stem(abs(XK);FFT

20、 Matlab 處理結果第三章 DFT與FFT計算速度比較分析3.1 FFT與直接計算 DFT的比較對于N點的DFT,需要計算的復數乘法和復數加法次數分別是N2和N(N-1)。由前面介紹的按時間抽取的FFT運算流圖可見，每一級都由 N2個蝶形運算構成。因此每一級運算都需要 N/2次復乘和N次復加(每個結加減各一次)。這樣m級運算總共需要NN ,復乘數mFm log 2 N22復加數年N m NlogNFFT計算量與 DFT計算量比較NDFTFFTDFT的乘法 次數與FFT 的乘法次數 之比DFT的加法 次數與FFT 加法次數之 比復數乘法次 數復數加法次 數復數乘法次 數復數加法次 數2421

21、241416124841.58645612245.32.31625624032648.03.753210249928016012.86.2644096403219238421.310.5128163841625644889636.618.125665536652801024204864.031.951226214426163223044608113.856.8102410485761047552512010240204.8102.32048419430441922561126422528372.4186.1409616777216167331202457649152682.7341.38192

22、6710886467100672532481064961260.3630.03.2 FFT 與 DFT 運算時間Matlab 程序3.2.1 隨機序列的DFT 計算時間程序Guide 程序下的主函數set(handles.TITLE,'string','calculating');Nmax=str2num(get(handles.N,'string');axes(handles.FFT_axes)cla;axes(handles.DFT_axes)cla;fft_time=zeros(1,Nmax);for n=1:1:Nmaxx=rand(1,

23、n);t=clock;my_fft(x);fft_time(n)=etime(clock,t);endn=1:1:Nmax;axes(handles.FFT_axes)plot(n,fft_time,'.');xlabel('N');ylabel(' 時間 ( 單位 : 秒 )')title('FFT 執行時間')my_dft_time=zeros(1,Nmax);for n=1:1:Nmaxx=rand(1,n);t=clock;my_dft(x,n);my_dft_time(n)=etime(clock,t);endn=1:1

24、:Nmax;axes(handles.DFT_axes)plot(n,my_dft_time,'.');xlabel('N');ylabel(' 時間 ( 單位 : 秒 ) ')title(' DFT 執行時間')set(handles.TITLE,'string','DFT_AND_FFT_TIME');3.2.2 Matlab實驗結果:Guide程序界面運行結果Nmax =1284-J-wtqNmax =512irn4行口可thNFFlIMtltMN=1024工彈o5 4 3-31 曲aQ 口

25、£a*lrl-_teN=2048ijfrFrihri+ia3.2.2 分析兩者運算時間的差異:由前面介紹的按時間抽取的FFT運算流圖可見，每一級都由N,2個蝶形運算構成。因此每一級運算都需要N:2次復乘和N次復加(每個結加減各一次)。這樣m級運算總共復乘數mFN ,m log2 N(4.3.45)復加數aF實際計算量和這個數字稍有出入，因為l 1l 1n l1 2 4222 1 Nn 0m Nlog2N(4.3.46)wN0 i 由運算流圖可見，這種情況共有1次,WN21, Wn4 j ,這幾個系數是都不用乘法運算的，但這種情況在直接計算DFT中也是有的，且當 N較大時，這些影響也較小。所以為了統一作比較我們不考慮以上特例。綜上所述，可以得出如下結論;NlogzN成正比的，而直接計算按時間抽取法所需的復乘數和復加數都是與DFT時所需的復乘數為 N2,復數加為N(N-1)次。當N>1 時,DIT-FFT算法比直接計算 DFT的運算次數大大