使用FFTW作FFT變換

通常FFT分析，都是求一個時間歷程的頻率成分，

用FFTW中的

fftw_plan_dft_r2c_1d

輸入是一個實數列向量in[N]，

輸出是一個複數空間二維矩陣out[N][2]，第1列是實數、第2列是複數。

 

求幅值時，要對複數向量的每一個，取向量模，乘以2，再除以點數N。

注意：

1. 若是採樣時間不到 1 s, 就設法補成 1 s。

一種方法是補零，（參考 彭真明 老師博文 link2 ）

還有一種方法，若是已知頻率大概範圍，折算出週期後，用正週期的時間歷程估計。

 

2. 若是有直流項（平均值大於零），須要把直流項濾掉。

 

fftw_complex *out = NULL;

fftw_plan p;

in  = (double *)fftw_malloc(sizeof(double) * len);

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

p = fftw_plan_dft_r2c_1d(len, in, out, FFTW_ESTIMATE);

fftw_execute(p);

for (int i = 0; i < len; i++)

{

amp  = sqrt(out[i][0]*out[i][0] + out[i][1]*out[i][1]);

amp  = 2.0* amp/len;

freq  =  i*1.0/t1;

}