2

我正在尝试使用 FFTW 和 Matlab 进行相同的 FFT。我使用 MEX 文件来检查 FFTW 是否良好。我认为我的一切都是正确的，但是：

  1. 我从 FFTW 得到荒谬的价值观，
  2. 在同一输入信号上多次运行 FFTW 代码时，我没有得到相同的结果。

有人可以帮我把 FFTW 弄对吗？

--

编辑 1：我终于弄清楚出了什么问题，但是...... FFTW 非常不稳定：我在 5 次中得到了正确的频谱 1 次！怎么来的？另外，当我做对时，它没有对称性（这不是一个非常严重的问题，但这太糟糕了）。

--

这是比较两者的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;
}
4

2 回答 2

1

Matlab 使用 fftw 库来执行其 ffts，在我的平台（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];
}
于 2013-06-16T02:27:37.677 回答
0

应该“p = fftw_plan_dft_r2c_1d(nfft, s, out, 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()

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

于 2014-07-07T22:01:00.107 回答