使用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