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（）

解决我的问题。也许它也可以帮助你。