1

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

我的代码如下：

Matlab函数：

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 向后：

正如 FFTW 文档中所建议的，此 MEX 函数与上述函数的不同之处仅在于切换指针pr <-> pi、函数和计划pr2 <-> pi2的CreatePlan执行。

如果我跑

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 四核）上使用 MATLAB 2011b。

4

1 回答 1

4

FFTW 计算非标准化变换，请参见此处： http ://www.fftw.org/doc/What-FFTW-Really-Computes.html

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

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

于 2012-11-07T16:55:16.027 回答