Matlab FFT和FFTW

Question

我正在嘗試使用FFTW和Matlab進行相同的FFT。我使用MEX文件來檢查FFTW是否良好。我想我擁有了一切正確的，但：

1. 我從FFTW得到荒謬的價值觀，
2. 相同的輸入信號運行FFTW碼好幾次，當我沒有得到相同的結果。

有人可以幫我搞定FFTW嗎？

-

編輯1：我終於想通了，什麼是錯的，但是... FFTW很不穩定：我得到正確的光譜1的時間超過了5！ 怎麼回事？另外，當我做對了，它沒有對稱性（這不是一個非常嚴重的問題，但那太糟糕了）。

-

這裏是Matlab代碼來比較兩種：

fs = 2000;     % sampling rate 
T = 1/fs;      % sampling period 
t = (0:T:0.1);    % time vector 

f1 = 50;      % frequency in Hertz 
omega1 = 2*pi*f1;    % angular frequency in radians 

phi = 2*pi*0.25;    % arbitrary phase offset = 3/4 cycle 
x1 = cos(omega1*t + phi);  % sinusoidal signal, amplitude = 1 

%% 

mex -I/usr/local/include -L/usr/local/lib/ -lfftw3 mexfftw.cpp 

N=256; 
S1=mexfftw(x1,N); 
S2=fft(x1,N); 
plot(abs(S1)),hold,plot(abs(S2),'r'), legend('FFTW','Matlab') 

這裏是MEX文件：

/********************************************************************* 
* mex -I/usr/local/include -L/usr/local/lib/ -lfftw3 mexfftw.cpp 
* Use above to compile ! 
* 
********************************************************************/ 
#include <matrix.h> 
#include <mex.h> 

#include "fftw3.h" 

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 

//declare variables 
mxArray *sig_v, *fft_v; 
int nfft; 

const mwSize *dims; 
double *s, *fr, *fi; 
int dimx, dimy, numdims; 

//associate inputs 
sig_v = mxDuplicateArray(prhs[0]); 
nfft = static_cast<int>(mxGetScalar(prhs[1])); 

//figure out dimensions 
dims = mxGetDimensions(prhs[0]); 
numdims = mxGetNumberOfDimensions(prhs[0]); 
dimy = (int)dims[0]; dimx = (int)dims[1]; 

//associate outputs 
fft_v = plhs[0] = mxCreateDoubleMatrix(nfft, 1, mxCOMPLEX); 

//associate pointers 
s = mxGetPr(sig_v); 
fr = mxGetPr(fft_v); 
fi = mxGetPi(fft_v); 

//do something 
double *in; 
fftw_complex *out; 
fftw_plan p;   

in = (double*) fftw_malloc(sizeof(double) * dimy); 
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * nfft); 
p = fftw_plan_dft_r2c_1d(nfft, s, out, FFTW_ESTIMATE); 

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

for (int i=0; i<nfft; i++) { 
    fr[i] = out[i][0]; 
    fi[i] = out[i][1]; 
} 

fftw_destroy_plan(p); 
fftw_free(in); 
fftw_free(out); 

return; 
} 

Answer 1

Matlab的使用FFTW庫，以執行其的FFT，在我的平臺上（Mac OS），這會導致鏈接器出現問題，因爲mex會用Matlab的fftw版本替換所需的庫。使用mex「-I/usr/local/include /usr/local/lib/libfftw3.a mexfftw.cpp」來避免此庫的靜態鏈接。 fftw_plan_dft_r2c_1d的輸入不會被破壞，因此您不需要重複輸入（注意：對於fftw_plan_dft_c2r_1d不是這樣）。輸出的大小爲nfft/2 + 1，因爲真實fft的輸出是Hermitian。因此，要得到完整的輸出使用：

for (i=0; i<nfft/2+1; i++) { 
    fr[i] = out[i][0]; 
    fi[i] = out[i][1]; 
} 
for (i=1; i<nfft/2+1; i++) { 
    fr[nfft-i] = out[i][0]; 
    fi[nfft-i] = out[i][1]; 
} 

Answer 2

應該 「P = fftw_plan_dft_r2c_1d（N FFT，S，出，FFTW_ESTIMATE）;」

是 「p = fftw_plan_dft_r2c_1d（nfft，in，out，FFTW_ESTIMATE）;」。

'in'是16字節對齊的，但's'可能不是。

我不確定是否會導致問題。我有一個關於FFTW的類似代碼， 它有時會給我正確的結果和其他時間NaN。此外，我試圖用python ctypes測試 我的代碼，實際上它具有相同的行爲。

最後，我發現這個職位Checking fftw3 with valgrind，它可以幫助我。 對於我來說，問題是即使在程序終止後，仍然沒有釋放出 的FFTW中的堆存儲空間。

fftw_cleanup（）

解決我的問題。也許它也可以幫助你。