Matlab FFT和FFTW

Question

我试着用FFTW和Matlab实现同样的FFT。我使用MEX文件来检查FFTW是否良好。我认为一切都是正确的，但是：

1. 我从自由工联那里得到荒谬的价值观，
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库来执行其fftw，在我的平台(Mac )上，这会导致链接器的问题，因为mex用Matlab的fftw版本替换了所需的库。为了避免这个静态链接使用mex "-I/usr/local/include /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(nfft，s，out，FFTW_ESTIMATE);“

be "p = fftw_plan_dft_r2c_1d(nfft，in，out，FFTW_ESTIMATE);“

“‘in”是16字节对齐的，但“s”可能不对。

我不知道这是否会引起问题。我有一个类似的关于FFTW的代码，它有时给我正确的结果，有时给我NaN。此外，我尝试用python ctype测试我的代码，实际上它也有相同的行为。

最后，我找到了这个帖子Checking fftw3 with valgrind，它帮助了我。对我来说，问题是在FFTW中保留堆存储，即使在程序终止之后，堆存储也不会被释放。

fftw_cleanup()

解决了我的问题。也许它也能帮到你。