Matlab FFT和FFTW [英] Matlab FFT and FFTW

问题描述

我正在尝试使用FFTW和Matlab进行相同的FFT.我使用MEX文件检查FFTW是否良好.我想我一切都正确，但:

I'm trying to have the same FFT with FFTW and Matlab. I use MEX files to check if FFTW is good. I think I have everything correct but :

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

有人可以帮助我正确设置FFTW吗?

Can someone help me get FFTW right?

-

我终于弄清楚出了什么问题，但是... FFTW非常不稳定:我在5次中得到1次正确的频谱！ 怎么会?另外，如果我做对了，它就没有对称性了(这不是一个很严重的问题，但这太糟糕了).

EDIT 1 : I finally figured out what was wrong, BUT... FFTW is very unstable : I get the right spectrum 1 time out of 5 ! How come? Plus, when I get it right, it doesn't have symmetry (which is not a very serious problem but that's too bad).

-

这是Matlab代码来比较两者:

Here is the Matlab code to compare both :

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;
}

推荐答案

Matlab在我的平台(Mac OS)上使用fftw库执行其ftft，这导致链接器出现问题，因为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而言，它不是true).输出的大小为nfft/2 + 1，因为实际fft的输出为Hermitian.因此，要获得完整的输出，请使用:

Matlab uses the fftw library to perform its ffts, on my platform (Mac OS) this leads to issues with the linker as mex replaces the desired library with Matlab's version of fftw. To avoid this static link to the library using mex "-I/usr/local/include /usr/local/lib/libfftw3.a mexfftw.cpp". The input of fftw_plan_dft_r2c_1d is not destroyed so you don't need to duplicate the input (note: not true for fftw_plan_dft_c2r_1d). The output has size nfft/2+1 as the output of a real fft is Hermitian. So to get the full output use:

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];
}

这篇关于Matlab FFT和FFTW的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！