FFT演算法實現與分析MATLAB

FFT演算法實現

厚

2.1實驗目的

I、加深對快速傅立葉變換的理解。
II、掌握 FFT 演算法及其程式的編寫。
III、掌握演算法效能評測的方法。
IV、熟悉MatLab程式設計。

2.2實驗原理

一個連續訊號Xa(t)的頻譜可以用它的傅立葉變換表示為：

如果對該訊號進行理想取樣，可以得到取樣序列：

同樣可以對該序列進行z變換，其中T為取樣週期：

當z=e^jω的時候，我們就得到了序列的傅立葉變換：

其中w稱為數字頻率，它和模擬域頻域的關係為：

其中fs是取樣頻率。上式說明數字頻率是模擬頻率對取樣率 fs 的歸一化。同模擬域的情況相似，數字頻率代表了序列值變化的速率，而序列的傅立葉變換稱為序列的頻譜。序列的傅立葉變換和對應的取樣訊號頻譜具有下式的對應關係。

即序列的頻譜是取樣訊號頻譜的週期延拓。從上式可以看出，只要分析取樣序列的頻譜，就可以得到相應的連續訊號的頻譜。
在各種訊號序列中，有限長序列在數字訊號處理中佔有很重要的地位。無限長的序列也往往可以用有限長序列來逼近。對於有限長的序列我們可以使用離散傅立葉（DFT），這一變換可以很好地反應序列的頻域特性，並且容易利用快速演算法在計算機上實現當序列的長度是 N 時，我們定義離散傅立葉變換為：

其中W_N=e^(-2Πj/N)，它的反變換定義為：

若直接計算DFT變換，整個DFT運算需要4N^2次實數相乘和2N(2N-1)次實數相加。
所以直接計算乘法次數與加法次數都和N^2成正比。例如N=10點的DFT，需要100次複數相乘，而N=1024時則需要1,048,576即一百多萬次複數乘法運算。這對於實時性要求很強的訊號處理來說，必將對計算速度有十分嚴苛的要求。為此，FFT作為對DFT的改進誕生。
快速傅立葉變換FFT並不是與DFT不相同的另一種變換，而是在DFT計算規律上建立的一種減少運算次數的快速演算法。常用FFT是以基-2的，長度N=2^M。運算效率高，程式簡單，使用方便。本實驗就使用以2為基實現FFT。演算法流程圖可以用蝶形演算法來表示，以8點的基2-FFT演算法為例：

每個蝶形運算為：

可以看到，每個蝶形運算都可原位運算。
當需要進行變換的序列長度不是2 的整數次方的時候，為了使用以 2 為基的 FFT，可以用末尾補零的方法，使其長度延長至 2 的整數次方。IFFT一般可以通過 FFT 程式來完成，只要對 X(k)取共軛，進行 FFT 運算，然後再取共軛，並乘以因子 1/N，就可以完成 IFFT。

2.3實驗內容

1、演算法實現
以對高斯序列進行FFT程式碼為例：

①、使用bitrevorder()函式對待變換訊號順序進行調整，存入x1中。例如M=8時
原順序：000(0)，001(1)，010(2)，011(3)，100(4)，101(5)，110(6)，111(7)

對二進位制翻轉之後為：

新順序：000(0)，100(4)，010(2)，110(6)，001(1)，101(5)，011(3)，111(7)

對比實驗原理中M=8的FFT輸入序列發現兩者順序一致。
②、n=1:m1表示一共有m1層運算，比如M=8時有3層。
k=1:M/(2^n)表示第n層分為k部分，比如M=8的第一層分為4個單獨蝶形運算，第二層分為2個兩兩蝶形運算，第三部分為1個四四蝶形運算。
m表示第k部分的第m個蝶形運算。
③、因為蝶形運算是原位運算，就不需要另外開闢空間，計算結果仍然存在原位置。
④、繪出編寫的FFT計算結果與MATLAB的FFT計算結果，以及兩者的差。
2、選取實驗 1 中的典型訊號序列驗證演算法的有效性
I、三角波

II、反三角波

III、高斯序列

IV、衰減正弦序列

V、單位脈衝序列

VI、矩形窗序列

VII、理想取樣序列

注意每個序列自己編寫的FFT計算結果和MATLAB的FFT計算結果，兩者相差數量級
在10^(-16)到10(-14)。說明編寫的FFT演算法正確性沒有問題。

3、對所編制的FFT演算法進行效能評估
演算法的評估首先是其正確性，是否能夠完成預期功能決定了該演算法是否有意義，
上一部分已經通過典型訊號序列驗證了編寫的FFT演算法的有效性。
這裡接下來主要從時間複雜度和空間複雜度兩方面來進行評價。空間複雜度是指該演算法執行過程需要佔用多少記憶體空間，隨著半導體產業的發展，記憶體空間變得越來越廉價，空間複雜度對演算法效能的影響也越來越小。人們往往根據時間複雜度來評價一個演算法的效能。而時間複雜度主要依賴於演算法的計算次數。已知基2-FFT演算法的複數乘法次數為1/2 〖Nlog〗_2 N。相比於加法，乘法在運算中要複雜得多，佔用資源也更多，所以時間複雜度主要依賴乘法次數。從理論複雜度來看，FFT顯然比DFT的N^2更優。在MARLAB中檢視程式運算時間有以下辦法：

①、tic和toc命令組合
tic；
operation；
toc；
tic用來儲存當前時間，也就是operation開始執行時間，toc用來記錄程式完成時間。MATLAB會自動計算時間差並顯示（以秒為單位但能精確到小數點後6位，即us）。

②、etime(t1,t2)和clock配合
t1=clock；
operation；
t2=clock；
etime(t1,t2);
通過呼叫windows系統時鐘進行時間差計算得到執行時間，t1和t2之間的時間差。

③cputime函式
t1=cputime；
operation；
t2=cputime-t1；
使用cpu主頻計算執行時間差，得到程式執行時間。

tic/toc是MATLAB自身計數器，精度要高於後兩者。而且，如②呼叫系統時鐘計算時間差，這段時間中系統可能還有其他後臺程式。
MATLAB官方推薦使用tic/toc組合，When timing the duration of an event,use the tic and toc functions instead of clock or etime.所以接下來的評估程式執行時間，本實驗使用tic/toc命令。此外，程式執行時間和計算機本身的計算能力有著直接關係，以下資料都是在個人膝上型電腦測得。由於電腦屬於商務本，計算能力很有限，時間相對會稍長一些。
以上主要說明了FFT演算法評估方法和側重點，具體評估資料在下面的dofft與DFT、dofft與MATLAB-FFT的效能比較中給出。

2.4實驗報告要求

1、總結自己實現的FFT演算法時採用了哪些方法減少了運算量。
1）使用matlab的bitrevorder()函式實現二進位制翻轉，由於matlab的函式是基於更底層的的c語言編寫的，有很專業的優化，執行速度肯定更快。
2）儘量使用小迴圈套大迴圈，因為執行的跳轉原因，大迴圈單次執行時間優於小迴圈。
3）利用蝶形運算的原位性，使用同一個地址空間儲存變換前序列和變換後序列。
4）每相鄰計算的蝶形運算資料在地址上儘量連續，減少定址時間。
2、給出自己的FFT演算法與實驗1中的DFT演算法效能比較結果。
為避免運算時間過短不利於記錄，使用20次迴圈。dofft程式見2.3，其餘程式如下
DFT程式：

dofft測試程式：

DFT測試程式：

執行時間記錄如下：
N
演算法 8 32 256 512 1024 2048
dofft 1.524102 1.563602 1.820989 1.984634 2.814433 3.088558
DFT 0.367555 0.375953 0.945036 4.491770 22.746014 107.771205
作圖：

測試序列為理想取樣序列。為避免迴圈執行時，MATLAB程式在前一迴圈已在記憶體中開闢空間和留有資料，測試程式中使用了clear all命令來清除記憶體中的資料。這樣測得的時間更加準確。從記錄的執行時間中可以看到，當N比較小的時候，DFT執行時間比dofft時間更小，這是因為DFT演算法使用的是向量運算，而dofft中使用了迴圈。MATLAB本身對於向量計算的速度快於迴圈的計算速度。所以如果進一步優化dofft演算法，可以改用向量運算，避免迴圈。當N趨於更大時,DFT執行時間迅速上升，很快 和dofft執行時間不在一個數量級。這和DFT、FFT兩演算法的理論時間複雜度一致。
3、給出自己的FFT演算法和MATLAB中fft演算法效能比較結果。
採用與2中相同的測試方法，同樣使用20次迴圈。

fft程式：

fft測試程式：

執行時間記錄如下：
N
演算法 512 1024 2048 4096 8192
dofft 1.984634 2.814433 3.088558 3.579388 5.438975
FFT 0.349001 0.369433 0.365632 0.403413 0.334651
N
演算法 16384 32768 65536 131072 262144
dofft 13.307104 18.961649 28.585727 61.441558 125.271445
FFT 0.350225 0.431561 0.430717 0.445077 0.530932
作圖：

自己編寫的dofft在進行變換長度的橫座標下近似線性增長。而MATLAB本身的FFT基本上隨著N的增大，執行時間基本上沒有變化。顯然dofft效能比fft效能差。想必MATLAB中的fft函式進行了更多技巧和優化。

4、總結實驗中根據實驗現象得到的其他結論。

①實驗中測執行時間，當前後兩個程式執行時間相差不大時，可能時間大小有波動。比如MATLAB中的fft函式在做2048點計算時所用時間比1024點計算時間稍小，這與MATLAB當前佔CPU和電腦狀態相關。也會出現同一個程式在不同時刻測執行時間大小稍有差異，多測幾次時間取平均時間長。
②DFT在N較小時執行時間小於dofft，說明MATLAB更優於計算向量。所以編寫MATLAB程式時應儘量把for迴圈改為矩陣運算，儘量向量化。
③同樣的演算法，基於不同的程式設計，在執行速度上仍然會有很大的不同，比如在MATLAB中使用for迴圈編寫，MATLAB本身的FFT，和用C語言編寫的FFT在執行速度上都會有很大不同。所以提高程式設計技巧，瞭解程式具體流水線、地址開闢、迴圈巢狀等等如何對優化程式有很大意義。
④MATLAB帶有眾多功能強大，高優化水平的函式，在編寫程式時，儘可能查詢MATLAB有無相關函式，充分利用，以提高編寫的程式的執行效率，比如dofft中的bitreorder。
⑤測量執行時間的時候要把繪圖的函式註釋掉，否者會佔用程式大量執行時間。