MFTとMATLAB引数の問題を伴うFFTW

Question

MATLABからのFFT計算に使用されるスレッドの数を制御するために、FFTWライブラリを使用して次のC/MEXコードを書きました。このコードは、MATLABよりも遅いですが、計画者のFFTW_ESTIMATE引数で素晴らしい（複雑なFFTの順方向および逆方向）機能を果たします。しかし、FFTWプランナーを調整するためにFFTW_MEASURE引数に切り替えると、1つのFFTを前方に適用し、次に1つの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切り替え上記とは異なり、CreatePlanでpr2 <-> pi2 FFTWの文書に示唆されているように、計画の実行にあたっては、

私はそれぞれFFTW_MEASUREとFFTW_ESTIMATE引数付き

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

を実行した場合、私は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 quad core）でMATLAB 2011bを使用しています。

Answer 1

FFTWは正規化変換しませ計算し、ここを参照してください：あなたは、直接逆1に続いて変換を実行するときhttp://www.fftw.org/doc/What-FFTW-Really-Computes.html

は大雑把に言えば、あなたはの長さを乗じた バック入力（プラス丸め誤差）を取得しますあなたのデータ。

あなたはFFTW_ESTIMATE以外のフラグを使用して計画を作成すると、あなたの入力が上書きされます。 http://www.fftw.org/doc/Planner-Flags.html