Saturday, May 06, 2006

FFTW - 計算 FFT 的函式庫 (C 語言)

FFTW 是用 C 語言寫成的函式庫,用以計算 DFT 轉換。輸入的訊號可以是任意長度,可以是實數型式或是複數型式的訊號。當然,從它的名稱知道它是以 FFT 來實作。本來我已經有修改自 NUMERICAL RECIPES 的 fft 轉換程式,不過 NUMERICAL RECIPES 的授權限制滿大的,相較之下 FFTW 則是採用 GNU General Public License 授權,比較可以放心使用。

要在 VC 下使用 FFTW ,可以從這裡下載 zip 檔

安裝︰

解壓縮之後,開啟 命令提示字元 ,轉換目錄到解壓縮的目錄,再執行以下三個指令︰
lib /machine:i386 /def:libfftw3-3.def
lib /machine:i386 /def:libfftw3f-3.def
lib /machine:i386 /def:libfftw3l-3.def
這會在目錄下建立三個 lib 檔。

將 libfftw3l-3.dll, libfftw3f-3.dll, libfftw3-3.dll 這3個 dll 檔複製到 system32 目錄。
在 VC 專案中指定 libfftw3l-3.lib, libfftw3f-3.lib, libfftw3-3.lib 這3個lib檔及 fftw3.h 檔所在的目錄。

在 VC 專案中加入 libfftw3l-3.lib, libfftw3f-3.lib, libfftw3-3.lib 這3個 lib。


使用︰
首先要引入 fftw3.h 這個標頭檔,如果要做的是一維 DFT 轉換,可參考以下的程式碼(取自fftw的線上說明文件)

#include <fftw3.h>
...
{
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);
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
...
fftw_execute(p); /* repeat as needed */
...
fftw_destroy_plan(p);
fftw_free(in); fftw_free(out);
}

fftw_complex 是 FFTW 自訂的複數型態,是一個2個元素的浮點陣列 double[2],若 in 為 fftw_complex 所宣告的一個陣列
fftw_complex * in = new fftw_complex[10];
,則第一個元素的實部為 in[0][0], 虛部為 in[0][1]。
fftw_malloc 是類似於 malloc 的函數,之間的差別我還不太清楚,不過文件上有提到,使用其他開記憶體的方式也是可以的(例如: new),使用 fftw_malloc 的好處是,開啟的陣列元素(包括 fftw_complex 的實虛部)會初始化為0(並非如此)。

fftw_plan_dft_1d - 從函式名稱大概可以知道,這個函式是用來產生一個準備的工作,它並不會真的執行 DFT 轉換,真正的 DFT 轉換要再用 fftw_exectue 來執行。

多維度號的轉換只需要將 fftw_plan_dft_1d 轉為以下任一個函式。
fftw_plan fftw_plan_dft_2d(int nx, int ny,
fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
fftw_plan fftw_plan_dft_3d(int nx, int ny, int nz,
fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
fftw_plan fftw_plan_dft(int rank, const int *n,
fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
輸出訊號格式︰ fftw_plan_dft_1d : 若輸入訊號的長度為 32 ,則輸出的FFT訊號從第17個元素(從0起算)開始重複。 out[17][0] == out[15][0], out[17][1] == out[15][1] out[18][0] == out[14][0], out[18][1] == out[14][1] . . out[31][0] == out[2][0], out[31][1] == out[2][1] fftw_plan_dft_r2c_1d : 這個函式的傳入陣列為實數陣列,若輸入訊號長度為32,則輸出的FFT訊號會在 第17個元素之後均為0

28 comments:

AllenLi said...

您好,我用您的方式在vc++6上跑fftw,
但是在compile的時候有錯誤,
因為我對vc++不熟,也許在那邊出了錯,
不知道可否請您給我詳細的步驟,
越細越好,萬分感謝,
您的舉手之勞,將能讓我們節省精力與時間。

Chinson said...

Hi Allen,
何不請您寫一下您詳細的設定步驟,我再來看看可能哪邊有問題呢?

Allen said...

您好,首先先建立一個project,選擇win32 console applocation,
再選an empty project,然後新增c/c++ source file(取消add to project)
程式碼:
#include iostream.h(這邊不能用大於小於所以省略)
#include fftw3.h

int main()
{
// 資料宣告
int N = 256;
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);
// 製作輸入資料
for( int i=0; i N; i++){
in[i][0] = i;
in[i][1] = i+1;
}
// 產生最佳化的程式碼
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
// 執行這個最佳化的程式碼
fftw_execute(p); /* repeat as needed */
// 釋放空間
fftw_destroy_plan(p);
fftw_free(in); fftw_free(out);
cout << out[0][0];
}

然後選project中的add to project,將fftw3.h、libfftw3-3.lib、libfftw3l-3.lib、libfftw3f-3.lib加入,再compile,得到錯誤訊息:
--------------------Configuration: fftw - Win32 Debug--------------------
Linking...
LINK : warning LNK4001: no object files specified; libraries used
LINK : error LNK2001: unresolved external symbol _mainCRTStartup
Debug/fftw.exe : fatal error LNK1120: 1 unresolved externals
Error executing link.exe.

fftw.exe - 2 error(s), 1 warning(s)
請幫忙看一下,謝謝!

Chinson said...

"然後新增c/c++ source file(取消add to project)" -> 為何要取消? 要跟專案一起編譯的 cpp 檔應該要加進專案

"然後選project中的add to project,將fftw3.h、libfftw3-3.lib、libfftw3l-3.lib、libfftw3f-3.lib加入" -> 這邊 3 個 lib 檔並不是直接用 add to project 的,請在專案開啟時,從主選單的 "Project" - "Settings" 開啟設定視窗,並在 "Link" - "General" - "Object/Library modules" 內加入這3個檔案(以空白作分隔)

ps. 為了方便起見,fftw3.h、libfftw3-3.lib、libfftw3l-3.lib、libfftw3f-3.lib 這 4 個檔就直接放在專案目錄下,libfftw3l-3.dll, libfftw3f-3.dll, libfftw3-3.dll 這 3 個 dll 檔可以放在 C:\Windows\system32 目錄下。(也可以直接放在專案目錄下)

這是我看到的可能錯誤。

Allen said...

"然後新增c/c++ source file(取消add to project) -> 為何要取消? 要跟專案一起編譯的 cpp 檔應該要加進專案"
因為不取消沒辦法下一步

謝謝您幫我釐清有關lib跟dll該如何作用,我再繼續試試,再次謝謝!

Anonymous said...

请问:
FFTW怎样得到傅立叶级数的系数
f(x)=a[0]+sigma(a[n]cosnx+b[n]sinnx)

就是如何得到a[n]和b[n]啊?

Chinson said...

建議你看一下DFT、FFT的定義,教科書上都有詳細的說明,Wiki 上也有一些介紹

ayumyss said...

請問你知道Dev C的fftw怎麼安裝嗎.謝謝^^

Chinson said...

to ayumyss

抱歉喔, 我沒用過 DEV C 耶
不過國外有人幫fftw弄了 DEV C 的 package, 如果熟 DEV C 的話或許可以試試

http://devpaks.org/details.php?devpak=23

Chinson

立源 said...

hi
您好
我用vc6照您的方法作
我在解壓縮之後,開啟 命令提示字元 ,轉換目錄到解壓縮的目錄,再執行以下三個指令︰

lib /machine:i386 /def:libfftw3-3.def
lib /machine:i386 /def:libfftw3f-3.def
lib /machine:i386 /def:libfftw3l-3.def
這邊我在命令提示字元中無法執行您可否幫我一下

Chinson said...

To 立源
我猜測是系統找不到 lib.exe ,因為該檔案所在目錄可能不在系統的搜尋目錄內。
請用檔案總管的搜尋功能,找到 lib.exe 所在的目錄,並把該目錄加入到環境變數 "path" 中。
環境變數的修改請從 "控制台->效能及維護->系統->進階->環境變數->系統變數" 進行修改。
修改完後重新開啟新的 "命令提示字元" 視窗,應該就可以執行 lib 指令了。
如果執行時出現找不到某些 dll 的情況,也請找出該 dll 所在的目錄,並加入到 path 變數中。

ice said...

HI
我裝好了
但目前還不知如何使用
請問256及512長度的fft該如何使用
還有做完fft後輸出的資料是代表什麼
實部,虛部及每壹點代表什麼頻率
謝謝

Chinson said...

To ice,
內文中有如何操作長度 fftw 的方式,要進行長度 256 或 512 的操作應該沒什麼問題

fftw 的計算是屬於 DFT 的操作,如果要換算成頻率,要看你的來源訊號的取樣頻率等等,必須要了解 DTFT 與 DFT 之間的關係。

抱歉,因為太久沒用了,觀念不見得正確,不適合講太多,詳情還是請參考書上的說明比較好。

ps. Lathi 著作的 "Signal Processing & Linear Systems" 這本書在10.6節(p.641)有說明

ice said...

抱歉
我想直接內文引用的程式
但我看不太懂程式碼
輸入資料是 float data[256];
資料是實數
請問要改什麼地方??

Chinson said...

fftw_complex *in, *out;

將現有的 float data[256] 轉成 fftw_complex data_c[256]

ice said...

C:\>lib /machine:i386 /def:libfftw3-3.def
Microsoft (R) Library Manager Version 6.00.8168
Copyright (C) Microsoft Corp 1992-1998. All rights reserved.

LIB : fatal error LNK1104: cannot open file "libfftw3-3.def"

您好
之前可以執行的
但不知道為什麼他現在又出現錯誤

Chinson said...

libfftw3-3.def 這個檔案有放在 C:\ 目錄下嗎?

ice said...

#include "fftw3.h"
void main()
{
fftw_complex *in, *out;
fftw_plan p;
int N= 32;
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

for( int i=0; i<=N; i++)
{
in[i][0] = 1;
in[i][1] = 0;
}

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

fftw_execute(p); /* repeat as needed*/

fftw_destroy_plan(p);
fftw_free(in); fftw_free(out);
printf("%f",out[0][0]);

}

您好
這是我修改之後的程式碼
長度在10的時候有值
但是當我把N改成32
他就沒友值了
可以麻煩您幫我看一下嗎??

Chinson said...

這個程式中有2個錯誤,

1. for( int i=0; i<=N; i++)
 i 必須小於(<)N ,而不是小於或等於(<=)

2. ftw_free(in); fftw_free(out);
printf("%f",out[0][0]);

第一行已經把 out 這個陣列釋放掉了,無法在第2行再進行顯示

這2個問題是 C/C++ 基本語法的錯誤

皓量 said...

您好:
因為vc6會有compile的問題
所以我嘗試在vs2005使用fftw
但一直會有link error,類似下面的錯誤訊息

bilateral_filter.obj : error LNK2001: 無法解析的外部符號 "private: static unsigned int FFT::Support_3D::fftw_flags" (?fftw_flags@Support_3D@FFT@@0IA)

bilateral_filter.obj : error LNK2001: 無法解析的外部符號 "private: static bool FFT::Support_3D::wisdom_loaded" (?wisdom_loaded@Support_3D@FFT@@0_NA)


我想要執行的code是
http://people.csail.mit.edu/sparis/bf/#code
這裡的Fast bilateral filter • full kernel (FFT convolution)


不知您有沒有什麼建議可以解決此問題
謝謝您的幫忙!

Chinson said...

To 皓量 :
看起來像是fftw函式庫設定不正確,確認一下是否有照下面的步驟設定︰

在專案開啟時,從主選單的 "Project" - "Settings" 開啟設定視窗,並在 "Link" - "General" - "Object/Library modules" 內加入 libfftw3-3.lib、libfftw3l-3.lib、libfftw3f-3.lib 這3個檔案(以空白作分隔)

為了方便起見,fftw的相關檔案 .h, .lib, .dll 都先放在專案目錄下測試。

Bert said...

想請問一下 如果我的一個matrix a[x][y],我想做對每一列做fft,有沒辦法直接呼叫fftw 呢? 還是我得對一列都要轉成a1[y]然後呼叫fftw x 次呢。

Chinson said...

To Bert,
我想是的。

莊書瑜 said...

我在BCB上享用fft按照教學步驟做完跑出以下問題
[Linker Error] Unresolved external '_fftw_malloc' referenced from E:\BCB\PROJECTS\HOMEWORK1.OBJ
[Linker Error] Unresolved external '_fftw_plan_dft_2d' referenced from E:\BCB\PROJECTS\HOMEWORK1.OBJ
[Linker Error] Unresolved external '_fftw_execute' referenced from E:\BCB\PROJECTS\HOMEWORK1.OBJ
[Linker Error] Unresolved external '_fftw_destroy_plan' referenced from E:\BCB\PROJECTS\HOMEWORK1.OBJ
[Linker Error] Unresolved external '_fftw_free' referenced from E:\BCB\PROJECTS\HOMEWORK1.OBJ
不知道能否幫忙解決

莊書瑜 said...
This comment has been removed by the author.
Chinson Yeh said...

看起來是fftw的lib沒設好。

莊書瑜 said...

看起來是fftw的lib沒設好。

///////////////////////////////////////

不知道您所謂的lib檔沒設好是甚麼意思?
是我沒將lib檔放入專案嗎?

Project->Add to project ...
我是這樣放進專案中的不知道是否有錯?

Chinson Yeh said...

請參考文章內新增的圖示說明。