FFTW с проблемами аргументов MEX и MATLAB

Я написал следующий код C/MEX, используя библиотеку FFTW, чтобы управлять количеством потоков, используемых для вычисления FFT из MATLAB. Код отлично работает (сложное БПФ вперед и назад) с аргументом FFTW_ESTIMATE в планировщике, хотя он медленнее, чем MATLAB. Но когда я переключаюсь на аргумент FFTW_MEASURE для настройки планировщика FFTW, оказывается, что применение одного БПФ вперед, а затем одного БПФ назад не возвращает исходное изображение. Вместо этого изображение масштабируется с коэффициентом. Использование FFTW_PATIENT дает мне еще худший результат с нулевыми матрицами.

Мой код выглядит следующим образом:

Функции Matlab:

БПФ вперед:

function Y = fftNmx(X,NumCPU)   

if nargin < 2
    NumCPU = maxNumCompThreads;
    disp('Warning: Use the max maxNumCompThreads');
end
Y = FFTN_mx(X,NumCPU)./numel(X);

БПФ назад:

function Y = ifftNmx(X,NumCPU)

if nargin < 2
    NumCPU = maxNumCompThreads;
    disp('Warning: Use the max maxNumCompThreads');
end

Y = iFFTN_mx(X,NumCPU);

Мекс-функции:

БПФ вперед:

# include <string.h>
# include <stdlib.h>
# include <stdio.h>
# include <mex.h>
# include <matrix.h>
# include <math.h>
# include </home/nicolas/Code/C/lib/include/fftw3.h>

char *Wisfile = NULL;
char *Wistemplate = "%s/.fftwis";
#define WISLEN 8

void set_wisfile(void)
{
    char *home;
    if (Wisfile) return;
    home = getenv("HOME");
    Wisfile = (char *)malloc(strlen(home) + WISLEN + 1);
    sprintf(Wisfile, Wistemplate, home);
}


fftw_plan CreatePlan(int NumDims, int N[], double *XReal, double *XImag, double *YReal, double *YImag)
{
  fftw_plan Plan;
  fftw_iodim Dim[NumDims];
  int k, NumEl;
  FILE *wisdom;

  for(k = 0, NumEl = 1; k < NumDims; k++)
  {
    Dim[NumDims - k - 1].n = N[k];
    Dim[NumDims - k - 1].is = Dim[NumDims - k - 1].os = (k == 0) ? 1 : (N[k-1] * Dim[NumDims-k].is);
    NumEl *= N[k];
  }

/* Import the wisdom. */
  set_wisfile();
  wisdom = fopen(Wisfile, "r");
  if (wisdom) {
    fftw_import_wisdom_from_file(wisdom);
    fclose(wisdom);
  }

  if(!(Plan = fftw_plan_guru_split_dft(NumDims, Dim, 0, NULL, XReal, XImag, YReal, YImag, FFTW_MEASURE *(or FFTW_ESTIMATE respectively)* )))
    mexErrMsgTxt("FFTW3 failed to create plan.");

/* Save the wisdom.  */
  wisdom = fopen(Wisfile, "w");
  if (wisdom) {
    fftw_export_wisdom_to_file(wisdom);
    fclose(wisdom);
  }

  return Plan;
}


void mexFunction( int nlhs, mxArray *plhs[],
              int nrhs, const mxArray *prhs[] )
{
  #define B_OUT     plhs[0]

  int k, numCPU, NumDims;
  const mwSize *N;
  double *pr, *pi, *pr2, *pi2;
  static long MatLeng = 0;
  fftw_iodim Dim[NumDims];
  fftw_plan PlanForward;
  int NumEl = 1;
  int *N2;

  if (nrhs != 2) {
      mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
                "Two input argument required.");
  }

  if (!mxIsDouble(prhs[0])) {
      mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
                "Array must be double");
  }

  numCPU = (int) mxGetScalar(prhs[1]);
  if (numCPU > 8) {
      mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
                "NumOfThreads < 8 requested");
  }

  if (!mxIsComplex(prhs[0])) {
      mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
                "Array must be complex");
  }


  NumDims = mxGetNumberOfDimensions(prhs[0]);
  N = mxGetDimensions(prhs[0]);
  N2 = (int*) mxMalloc( sizeof(int) * NumDims);
  for(k=0;k<NumDims;k++) {
    NumEl *= NumEl * N[k];
    N2[k] = N[k];
  }

  pr = (double *) mxGetPr(prhs[0]);
  pi = (double *) mxGetPi(prhs[0]);

  //B_OUT = mxCreateNumericArray(NumDims, N, mxDOUBLE_CLASS, mxCOMPLEX);
  B_OUT  = mxCreateNumericMatrix(0, 0, mxDOUBLE_CLASS, mxCOMPLEX);
  mxSetDimensions(B_OUT , N, NumDims);
  mxSetData(B_OUT , (double* ) mxMalloc( sizeof(double) * mxGetNumberOfElements(prhs[0]) ));
  mxSetImagData(B_OUT , (double* ) mxMalloc( sizeof(double) * mxGetNumberOfElements(prhs[0]) ));

  pr2 = (double* ) mxGetPr(B_OUT);
  pi2 = (double* ) mxGetPi(B_OUT);

  fftw_init_threads();
  fftw_plan_with_nthreads(numCPU);  
  PlanForward = CreatePlan(NumDims, N2, pr, pi, pr2, pi2);
  fftw_execute_split_dft(PlanForward, pr, pi, pr2, pi2);
  fftw_destroy_plan(PlanForward);
  fftw_cleanup_threads();

}

БПФ назад:

Эта MEX-функция отличается от приведенной выше только переключением указателей pr <-> pi и pr2 <-> pi2 в функции CreatePlan и выполнением плана, как это предлагается в документации FFTW.

Если я побегу

A = imread('cameraman.tif');
>> A = double(A) + i*double(A);
>> B = fftNmx(A,8);
>> C = ifftNmx(B,8);
>> figure,imagesc(real(C))

с аргументами FFTW_MEASURE и FFTW_ESTIMATE соответственно я получаю этот результат.

Интересно, это из-за ошибки в моем коде или в библиотеке? Я пробовал разные вещи вокруг мудрости, экономя, а не экономя. Использование мудрости, производимой автономным инструментом FFTW, для создания мудрости. Я не видел никаких улучшений. Кто-нибудь может подсказать, почему это происходит?

Дополнительная информация:

Я компилирую код MEX, используя статические библиотеки:

mex FFTN_Meas_mx.cpp /home/nicolas/Code/C/lib/lib/libfftw3.a /home/nicolas/Code/C/lib/lib/libfftw3_threads.a -lm

Библиотека FFTW не была скомпилирована с:

./configure  CFLAGS="-fPIC" --prefix=/home/nicolas/Code/C/lib --enable-sse2 --enable-threads --&& make && make install

Я пробовал разные флаги без успеха. Я использую MATLAB 2011b на 64-битной Linux-станции (четырехъядерный AMD Opteron).