0

I tried to calculate fft in 2D using the following code, however, the results are incorrect or Maltab displays an error message and closes automatically. I couldn't find where is the mistake. Could you please help me ? Here is the code :

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "C:\exemple\fftw3.h"
#include <mex.h>


void DivideArray(double *Data, int NumEl, double Divisor)
{
  int n;
  for(n = 0; n < NumEl; n++)
    Data[n] /= Divisor;
}


void FFT2DIntervaled(int M ,int N, double *X, double *Y, int Sign)
{
  fftw_plan Plan;
  if(!(Plan = fftw_plan_dft_2d(M, N, (fftw_complex *)X, (fftw_complex *)Y, Sign, FFTW_ESTIMATE)))
    mexErrMsgTxt("FFTW3 failed to create plan");
  fftw_execute(Plan);
  fftw_destroy_plan(Plan);
  if(Sign == 1)
    DivideArray(Y, 2*M*N, M*N);
}

void mexFunction(int nlhs,mxArray *plhs[] ,int nrhs ,const mxArray *prhs[])
{
  plhs[0] = mxCreateDoubleMatrix(mxGetM(X_IN),mxGetN(X_IN), mxCOMPLEX);
  FFT2DIntervaled(mxGetM(prhs[0]),mxGetN(prhs[0]),(double*)mxGetPr(prhs[0]),mxGetPr(plhs[0]),1);
}

N.B: I dont think it's duplicate because in my case Matlab crashes and closes automatically.


EDIT Taking into account your remarks , here is the new version of the code :

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "C:\exemple\fftw3.h"
#include <mex.h>

void DivideArray(double *Data, int NumEl, double Divisor)
{
    int n;
    for(n = 0; n < NumEl; n++)
        Data[n] /= Divisor;
}


void FFT2DIntervaled(int M,int NN, double *Xr, double *Xi,double *Yr, double *Yi, int Sign){
    int k,NumEl;
    fftw_plan Plan;
    fftw_iodim Dim[2];
    int N[2]={M,NN};
    for(k = 0, NumEl = 1; k < 2; k++)
    {
        Dim[1-k].n = N[k];
        Dim[1-k].is = Dim[1-k].os = (k == 0) ? 1 : (N[k-1] * Dim[2-k].is);
        NumEl *= N[k];}
    if(!(Plan = fftw_plan_guru_split_dft(2, Dim, 0, NULL,
            Xr, Xi, Yr, Yi, FFTW_ESTIMATE)))
        mexErrMsgTxt("FFTW3 failed to create plan.");
    if(Sign == -1)
        fftw_execute_split_dft(Plan, Xr, Xi, Yr, Yi);
    else
    {fftw_execute_split_dft(Plan, Xi, Xr, Yi, Yr);
     DivideArray(Yr, NumEl, NumEl);
     DivideArray(Yi, NumEl, NumEl);
    }

    fftw_destroy_plan(Plan);
    return;

}

void mexFunction(int nlhs,mxArray *plhs[] ,int nrhs ,const mxArray *prhs[])

{

    if (!mxIsDouble(prhs[0])) {
        mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
                "Array must be double");
    }

    plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[0]),mxGetN(prhs[0]), mxCOMPLEX);


    mxSetData(plhs[0], mxMalloc( sizeof(double) * mxGetNumberOfElements(prhs[0]) ));
    mxSetImagData(plhs[0], mxMalloc( sizeof(double) * mxGetNumberOfElements(prhs[0]) ));


    FFT2DIntervaled(mxGetM(prhs[0]),mxGetN(prhs[0]),(double*)mxGetPr(prhs[0]),(double*)mxGetPi(prhs[0]),mxGetPr(plhs[0]),mxGetPi(plhs[0]),1);
}

It still doesnt work. Could you help me to find the mistake?

ransa
  • 71
  • 10
  • 2
    You're casting a double pointer to a complex pointer and wondering why your program crashes? Read the documentation to mxGetPr and mxGetPi. – Cris Luengo Apr 10 '17 at 00:21
  • @cris-luengo You are right , I didn't noticed that . I changed my code . The newest version is on my first post . It still doesn't work. Could you help me to find where is the mistake ? – ransa Apr 10 '17 at 08:36

0 Answers0