利用ＦＦＴＷ實現快速傅里葉變換

摘要： 本文對FFTW的工作原理及機制作了詳細的分析，并通過實例講解了FFTW實現FFT的具體方法及性能的測試。

關鍵詞： FFTW 工作原理 使用方法 性能評測

引言

快速傅里葉變換(Fast Fourier transform簡稱FFT)作為時域和頻域轉換的分析工具，廣泛應用在數字通信、音頻信號處理、圖像處理、 譜估計、 雷達信號處理、地震勘探等各領域中。在FFT的軟件實現中，比較有名的庫有fftpack，mkl及FFTW等。

一、FFT及FFTW

FFT是離散傅里葉變換的快速算法，也可用于計算離散傅里葉變換的逆變換。快速傅里葉變換有廣泛的應用，如數字信號處理、計算大整數乘法、求解偏微分方程，等等。

對于復數序列{x(n)}(n=0，1，2，...，N-1)，則該復數序列的離散傅里葉變換為［１］：

X(K)= x(n)W，K=0，1，2，...N-1

其中W =e 。

直接按照上述變換公式的計算復雜度是O（n ）；快速傅里葉變換（FFT）可以計算出與直接計算相同的結果，但只需要計算復雜度O(nlog(n))。而實現FFT的算法有AMD Core Math Library (ACML)，intel mkl，cwplib，dfftpack，kissfft，mixfft，monnier，jmffte，numutils，FFTW3等，FFTW以其在各個平臺上的優秀表現，應用相當的廣泛。

FFTW是Fastest Fourier Transform in the West(西方快速Fourier變換)的縮寫。它是計算離散Fourier變換（DFT）的快速C程序的一個完整集合，它可計算一維或多維、實和復數據以及任意規模的DFT。FFTW還包含對共享和分布式存儲系統的并行變換。FFTW由麻省理工學院計算機科學實驗室超級計算技術組開發。

二、FFTW的工作原理

FFTW在各種平臺上都能實現相當高的性能，主是由于以下兩個FFTW工作特點［２］。

1.將長度為N的序列分解成長度為N 和N 的短序列進行計算。這些短序列的計算由一系列經過高度優化的C代碼片段來完成，這些代碼片段稱為solvers。solvers由planner來管理與組合起來完成計算，其中一種組合稱為一個plan。

2.對于同一個長度為N，其分解方式有多種，不同的組合路徑對于運算的速度影響不同，在不同的機器與平臺下其表現也會有所不同，因此需要確定哪種組合的效果最佳。FFTW通過對各種組合的評估，選擇出其中執行速度最快的一個plan。將此plan記錄下來，以后每次做同樣長度的變換時，只要使用相當的plan即可，這樣使FFTW能適應各種運行環境下的最優。

FFTW的核心工作模塊主要由planner，solver和plan來組成。其中solver是解決變換的一種獨立的方法；plan是由solver生成的包含計算變換的路徑和所有輸入輸出數組的指針的一個數據結構，可以供后面的算法使用；planner則是管理solver生成plan，并通過一定的算法找出執行速度最快的plan，并將其返回給用戶。具體的執行情況如下：

首先planner 對一系列的solver進行初始化，然后對給出一個指定長度變換，planner依次調用每個solver來生成plan，當此solver可以解決的話就返回一個plan指針，不能解決時則生成一個空指針。Planner對生成的各個plan進行執行并測量，選擇執行速度最快的plan返回給用戶。而用戶則調用返回的plan來執行傅里葉變換。可以將此包含執行路徑及組合算法的plan保存下來，重復使用。

由于FFTW的 planner是一些獨立的解決特定問題的plan的模塊。故同樣的planner可以用于復數的DFT，實數的DFT及其它的變換。

三、FFTW的使用方法

在使用FFTW之前應該配置好其開發環境，本文以VC6.0為例來講述其環境的配置方法。其步驟如下：

1.下載FFTW for windwos版本（ftp://ftp.fftw.org/pub/fftw/fftw-3.2-dll.zip）。

2.解壓fftw-3.2-dll.zip到一指定目錄，如c:＼fftw3。在此目錄下有編譯所需要的動態鏈接庫，庫的定義文件等，但缺少靜態鏈接的lib文件，可利用def文件和VC的命令來生成。

3.使用VC下的命令lib來生成VC靜態鏈接時需要的lib文件。如下命令將生成libfftw3-3.lib文件。

lib/def:libfftw3-3.def

4.在VC的包含文件路徑中加入fftw3頭文件的路徑，如c:＼fftw3。

5.在新建的工程中添加上生成的libfftw3-3.lib。

6.然后在工程就可以使用FFTW庫。

FFTW的使用過程要使用到的數據結構為fftw_complex 。它的定義如下：

typedef double fftw_complex;

其是一個復數，其數據的存儲格式為：部分存儲實數部分，部分存儲虛數部分。此數據結構由fftw_malloc和fftw_free函數分配與釋放。其函數原型為：

#include

void *fftw_malloc(size_t n);

void fftw_free(void *p);

在使用過程中使用到的函數有fftw_plan_dft_*(*號代表各種方式的變換)和fftw_execute函數。其中fftw_plan_dft_*可實現生成一維及多維的及時變換和復雜變換的plan。fftw_execute則根據生成的plan來執行變換。其函數原型為：

void fftw_execute(const fftw_plan plan);

FFTW函數的使用步驟如下：

首先，使用fttw_malloc分配變換中要使用的輸入輸出數組。

第二步，利用它生成一個plan，用來包含FFTW計算FFT是要包含的所有的數據。

當輸入輸出數組使用同一個數組時，將執行原位運算。

第三步，利用先前生成的plan作為參數，調用fftw_execute執行傅里葉變換。當需要多次執行時，可以重復調用fftw_execute來完成任務。

最后，當不使用時，可以調用fftw_destroy_plan，fftw_free函數將相應的plan和輸入輸出數組銷毀。

具體的使用實例如下所示：

/* 定義輸入輸出數組及fftw_plan*/

fftw_complex *in， *out;

fftw_plan p;

/* 為輸入輸出數組分配空間 */

in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

/*調用planner生成最優plan */

p = fftw_plan_dft_1d(N， in， out， FFTW_FORWARD， FFTW_ESTIMATE);

/*給輸入數組in添加數據*/

...

/*執行變換*/

fftw_execute(p); /* 根據需要可反復調用 */

/*執行完后在輸出數組out中取出變換后的輸出數據 */

...

/*不使用該plan，將其和輸入輸出數組銷毀 */

fftw_destroy_plan(p);

fftw_free(in); fftw_free(out);

四、FFTW的性能評測

在2.2 GHz 32位Dual Core AMD Opteron Processor 275處理器上，分別對AMD Core Math Library (ACML)，intel mkl，cwplib，dfftpack，kissfft，mixfft，monnier，jmffte，numutils，FFTW3對雙精度一維的不同長度的DFT進行測試。其測試的速度如圖1所示［３］。

從上圖可以看出其性能明顯要優于其它的FFT算法的執行速度。在實踐中使用FFTW將會顯著縮短算法的執行時間，提高整個系統的實時性能。

本人在做地震實時相關處理項目中使用FFTW做相關運算的處理，取得了顯著的性能提升。由此可見在實踐中使用FFTW實現FFT是切實可行的，并且能夠高效地執行。

參考文獻：

［1］劉益成，孫祥娥.數字信號處理.北京：電子工業出版社，2004.

［2］王剛.遙感圖像快速傅立葉變換新算法的研究與實現.中科院遙感應用研究所，2003.

［3］FFTW User Manual.http：//www.fftw.org.

注：“本文中所涉及到的圖表、注解、公式等內容請以PDF格式閱讀原文。”