相信內(nèi)行一看標(biāo)題就知道我要說誰了疑枯！沒錯，就是 FFTW 蛔六。這是一個求解快速傅里葉變換的開源庫荆永，就連大佬級別的Matlab也是用FFTW家的庫。當(dāng)然了国章，也有人知道這個庫的存在具钥，可能苦于這個庫還不足夠傻瓜，而長時間沒有親自體驗液兽。本文的目的就是一步一步道來FFTW的使用過程骂删。linux和mac系統(tǒng)下的使用太簡單了，問后稍微提一下就行抵碟。主要講解windows系統(tǒng)下的使用，我這里使用的環(huán)境是：
windows 10+Visual Studio 2017 Community版本
Windows
下載
獲取FFTW庫有兩種方式坏匪，一種是從 源碼編譯 拟逮；一種是直接使用 官方發(fā)布的dll，有 x63 和 win32 兩種版本可供選擇适滓。我推薦 第二種方式 敦迄。我測試的是x64版本，直接 下載 后解壓即可。解壓后的文件列表如下圖所示：
綠色方框內(nèi)的 .dll 文件和 .def 文件是兩個關(guān)鍵文件罚屋。還有類似的兩組苦囱，分別后面帶有 f 和 l 應(yīng)該表示數(shù)值精度的區(qū)別，同樣的用法脾猛。
生成lib
一般在VS中編寫C/C++程序調(diào)用dll時候撕彤，同時也需要lib文件作為連接庫文件。雖然FFTW下載文件里面沒有對應(yīng)的lib文件猛拴，但是沒關(guān)系羹铅，可以使用一個工具生成lib文件。
這個工具就是vs2017自帶的一個程序 lib.exe愉昆，位于 C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Tools\MSVC\14.15.26726\bin\Hostx64\x64 目錄（自己的系統(tǒng)下也可以按照類似的路徑去尋找）职员。
然后在命令行中切換到此目錄下，運行命令跛溉。
命令格式為：lib /out:保存lib的路徑\fftw3-3.lib /machine:x64或者win32 /def:.def文件的路徑\fftw3-3.def
運行命令： lib /out:D:\SoftwareProject\fftw-3.3.5-dll64\fftw3-3.lib /machine:x64 /def:D:\SoftwareProject\fftw-3.3.5-dll64\libfftw3-3.def 焊切。
或者直接在隨便一個目錄下新建一個.bat文件(比如createlib.bat)，里面寫入如下代碼芳室，然后雙擊運行即可
C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Tools\MSVC\14.15.26726\bin\Hostx64\x64\lib /out:D:\SoftwareProject\fftw-3.3.5-dll64\fftw3-3.lib /machine:x64 /def:D:\SoftwareProject\fftw-3.3.5-dll64\libfftw3-3.def
然后在代碼中填入的lib生成路徑(D:\SoftwareProject\fftw-3.3.5-dll64)下看到新生成的兩個文件: fftw3-3.lib和fftw3-3.exp∽ǚ荆現(xiàn)在就可以正式編寫代碼使用fftw進行快速傅里葉變換啦。
使用fftw
這里舉個一維傅里葉變換的例子渤愁，二維和三維的用法與一維的類似牵祟，細(xì)節(jié)可以參考 官方文檔 。
問題描述和模擬數(shù)據(jù)生成
（1） 問題： 提取一個疊加正弦信號的頻率抖格。
（2）模擬信號： 用代碼生成頻率為30 Hz和90Hz的兩個正弦信號的疊加作為模擬信號诺苹，可用下面這段matlab代碼生成，并保存到文件 signal_30_90.txt雹拄。
% create a 1-D signal with frequancy of 30Hz and 90Hz
f1=30; f2=90;
x=linspace(-3*pi,3*pi,2000);
y=sin(2*pi*f1*x)+sin(2*pi*f2*x);
figure(1);
plot(x,y);
fpout=fopen('signal_30_90.txt','w');
for i=1:length(x)
fprintf(fpout,'%.5E %.5E\n',x(i),y(i));
end
fclose(fpout);
編寫C++代碼使用fftw進行頻譜分析
前奏
（1）添加fftw頭文件
有兩種方式添加fftw的頭文件到自己的項目中：第一種方式是直接將解壓文件中的fftw3.h拷貝到項目目錄下收奔；第二種方式是在vs2017的項目屬性中添加包含目錄中加入fftw3.h所在的路徑。然后就可以在自己的代碼中用 #include "fftw3.h包含頭文件了滓玖，不然編譯器會提示錯誤：找不到頭文件坪哄。推薦第一種方式，因為畢竟文件很小很少势篡，而且這種方式方便移植程序給別人 翩肌。
（2）添加鏈接庫
與（1）類似，將libfftw3-3.lib拷貝到項目目錄下禁悠，然后在vs2017項目屬性中的連接庫下面添加libfftw3-3.lib念祭，操作過程如下圖所示：
在代碼中使用fftw
首先將fftw通過 #include "fftw3.h"將fftw的頭文件加入代碼中，然后就可以用fftw中的函數(shù)進行正碍侦、逆傅里葉變換粱坤。
下面代碼功能是從文件 signal_30_90.txt 中讀取信號隶糕，然后對其進行傅里葉變化求出頻譜并保存到文件 Spectrum_Singal.txt 中，結(jié)果見下圖
// TestFFT.cpp : 此文件包含 "main" 函數(shù)站玄。程序執(zhí)行將在此處開始并結(jié)束枚驻。
//
#define PI 3.141592653
#include <iostream>
#include <fstream>
using namespace std;
// 將fftw的使用進行一定的封裝，??使用過程更簡便
void FFT1d(vector<double> invector_r, fftw_complex* out_c);
//測試一維傅里葉變換
int testFFT1d();
int main()
{
//測試一維傅里葉變換：求解一個正弦疊加信號的頻譜
testFFT1d();
return 0;
}
void FFT1d(vector<double> invector_r, fftw_complex* out_c)
{
//1. 申明數(shù)據(jù)輸入輸出數(shù)組
fftw_complex *in;
//2. 申明fft執(zhí)行變量
fftw_plan p;//,q;
//3. 給出數(shù)據(jù)個數(shù)
int N = invector_r.size();
//4. 根據(jù)數(shù)據(jù)個數(shù)開辟動態(tài)數(shù)組空間
in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)* N);
//out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)* N);
//5. 給輸入數(shù)據(jù)數(shù)組賦值
for (int i = 0; i < N; i++)
{
in[i][0] = invector_r[i];
in[i][1] = 0.0;
//printf("%6.2f \n", in[i][0]);
}
//6. 調(diào)用fft執(zhí)行函數(shù)
p = fftw_plan_dft_1d(N, in, out_c, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p); /* repeat as needed*/
//8. 最后銷毀fft相關(guān)的變量
fftw_destroy_plan(p);
//fftw_destroy_plan(q);
fftw_free(in);
}
int testFFT1d()
{
string filename = "data\\signal_30_90.txt";
ifstream fin = ifstream(filename);
vector<double> t;
vector<double> amplitude;
vector<double> fvector;
vector<double> SpectrumVector;
double t0, amp0;
if (!fin)
{
cout << "打開文件失斨昕酢： " << endl;
return 0;
}
while (!fin.eof())
{
fin >> t0 >> amp0;
t.push_back(t0);
amplitude.push_back(amp0);
}
cout << t[0] << " " << amplitude[0] << " " << t.size() << endl;
//2. Calculate FFT
fftw_complex* fftout=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*t.size());
FFT1d(amplitude, fftout);
//3. calculate spectrum
GetSpectrum1d(t.size(), t[1]-t[0],fftout, fvector, SpectrumVector);
//write spectrum to file
ofstream fout = ofstream("data\\Spectrum_Singal.txt");
for (int i = 0; i < SpectrumVector.size(); i++)
{
fout<<fvector[i]<<"\t" << SpectrumVector[i] << endl;
}
fout.close();
//free pointer
fftw_free(fftout);
return 0;
}
