FFTW 与 MEX 和 MATLAB 参数问题

Question

我使用 FFTW 库编写了以下 C/MEX 代码来控制用于 MATLAB FFT 计算的线程数。代码在规划器中使用FFTW_ESTIMATE参数时效果很好(向前和向后复杂 FFT(，尽管它比 MATLAB 慢。但是，当我切换到FFTW_MEASURE参数来调整 FFTW 规划器时，事实证明，向前应用一个 FFT 然后向后应用一个 FFT 不会返回初始图像。相反，图像按一个因子缩放。使用FFTW_PATIENT会给我带来更糟糕的空矩阵结果。

我的代码如下：

矩阵函数：

FFT 转发：

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

FFT 向后：

function Y = ifftNmx(X,NumCPU)
if nargin < 2
    NumCPU = maxNumCompThreads;
    disp('Warning: Use the max maxNumCompThreads');
end
Y = iFFTN_mx(X,NumCPU);

墨西哥函数：

FFT 转发：

# 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();
}

FFT 向后：

此 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

我尝试了不同的标志但没有成功。我在 Linux 64 位工作站(AMD opteron quad core(上使用 MATLAB 2011b。

Answer 1

FFTW 计算非规范化转换，请参阅此处：http://www.fftw.org/doc/What-FFTW-Really-Computes.html

粗略地说，当您执行直接变换后跟逆变换时，您会得到将输入(加上舍入误差(乘以数据的长度。

使用 FFTW_ESTIMATE 以外的标志创建计划时，您的输入将被覆盖：http://www.fftw.org/doc/Planner-Flags.html