使用MEX和MATLAB参数问题的FFTW

Question

我使用FFTW库编写了以下C/MEX代码来控制用于从MATLAB计算FFT的线程数。该代码在规划器中的作用很好（复杂的FFT向前和向后）与FFTW_ESTIMATE参数，虽然它比MATLAB慢。但是，当我切换到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向后：

此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参数我得到this result。

我不知道这是由于我的代码或库中的错误。我围绕智慧尝试了不同的东西，节省不了。使用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。

Answer 1

FFTW计算未归一化改造，在这里看到： http://www.fftw.org/doc/What-FFTW-Really-Computes.html

粗略地说，当你执行直接变换，然后通过反向一个，你会得到 回输入（加上舍入误差）乘以长度你的数据。

当您使用比FFTW_ESTIMATE其他标志计划，您的输入将被覆盖： http://www.fftw.org/doc/Planner-Flags.html