Matlab FFT и FFTW

Question

Я пытаюсь иметь тот же БПФ с FFTW и Matlab. Я использую файлы MEX, чтобы проверить, хорош ли FFTW. Я думаю, что я все правильно, но:

1. Я получаю абсурдные значения от FFTW,
2. Я не получить те же результаты при выполнении кода FFTW несколько раз на том же входном сигнале.

Может кто-нибудь помочь мне получить FFTW правильно?

-

EDIT 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 для выполнения своих преобразований Фурье, на моей платформе (Mac OS) это приводит к проблемам с компоновщиком, поскольку mex заменяет нужную библиотеку версией fftw Matlab. Чтобы избежать этой статической ссылки на библиотеку с помощью 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 является эрмитовым. Таким образом, чтобы получить полное использование выходной:

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

В случае, если "р = fftw_plan_dft_r2c_1d (NFFT, s, из, FFTW_ESTIMATE);"

be "p = fftw_plan_dft_r2c_1d (nfft, in, out, FFTW_ESTIMATE);".

'in' - 16-байтовый выровненный, но 's' не может.

Я не уверен, что это вызовет проблему. У меня есть аналогичный код о FFTW, он иногда дает мне правильный результат и другие времена NaN. Кроме того, я попытался проверить свой код с ctypes python, и на самом деле он имеет такое же поведение.

Наконец-то я нашел это сообщение Checking fftw3 with valgrind, и это помогает мне. Для меня проблема заключается в резервировании хранилища кучи в FFTW, которые не освобождаются даже после завершения программы.

fftw_cleanup()

решает мою проблему. Может быть, это тоже может вам помочь.