/*******************************************************************************
 * r2fftMex.c   :   r2fft interface to Matlab
 * author       :   duweitao
 * usuage       :   compile this file and cfft.c cfft.h with 
 *                  "mex r2fftMex.c cfft.c cfft.h"
 *                  then run b = r2fftMex(a) , a must be a vector with length 
 *                  of 2^N , N >= 1 && N <= 15
 * date         :   2008-July-5th
 * version      :   0.1
 * modify       :   initial version
 ******************************************************************************/
#include "mex.h"
#include "cfft.h"

/* Input Arguments */
#define X_IN    (prhs[0])

/* Output Arguments */
#define Y_OUT   (plhs[0])

#if !defined(MAX)
#define	MAX(A, B)	((A) > (B) ? (A) : (B))
#endif

#if !defined(MIN)
#define	MIN(A, B)	((A) < (B) ? (A) : (B))
#endif


void r2fftMexTb(
  unsigned int  fftLen      ,
  double        *inReal     ,
  double        *inImag     , 
  double        *outReal    ,
  double        *outImag    )
{
  cfftdata_t cfftInReal, cfftInImag;
  cfftdata_t*   inCfftPtr;
  cfftdata_t*   outCfftPtr;

  unsigned int idx, err;
  handleCfft_t hCfft = NULL;
  inCfftPtr  = (cfftdata_t*) malloc(sizeof(cfftdata_t)*fftLen*2);
  memset(inCfftPtr, 0, sizeof(cfftdata_t)*fftLen*2);
  outCfftPtr = (cfftdata_t*) malloc(sizeof(cfftdata_t)*fftLen*2);
  memset(outCfftPtr, 0, sizeof(cfftdata_t)*fftLen*2);

  for(idx = 0; idx < fftLen; idx ++){
    cfftInReal = (cfftdata_t) inReal[idx];
    cfftInImag = (cfftdata_t) inImag[idx];
    inCfftPtr[2*idx]   = cfftInReal;
    inCfftPtr[2*idx+1] = cfftInImag;
  }

  err = CfftInit(fftLen, &hCfft);
  if(err){
    mexPrintf("%s\n", hCfft->errMsg);
    return;
  }
  err = CfftRun(inCfftPtr, outCfftPtr, hCfft);
  CfftClose(hCfft);


  for(idx = 0; idx < fftLen; idx ++){
    outReal[idx] = (double) outCfftPtr[2*idx];
    outImag[idx] = (double) outCfftPtr[2*idx + 1];
  }
  return;
} // r2fftMexTb()

void mexFunction( 
    int nlhs, mxArray *plhs[], 
    int nrhs, const mxArray*prhs[] )
     
{ 
    mxArray *InBuffer;
    double *xRealPtr, *xImagPtr, *yRealPtr, *yImagPtr, *inPtrR, *inPtrI;
    unsigned int m,n, isComplex, fftLen; 
    unsigned int idx; 
    
    /* Check for proper number of arguments */
    
    if (nrhs != 1) { 
    mexErrMsgTxt("ONE input arguments required."); 
    } else if (nlhs > 1) {
    mexErrMsgTxt("ONLY one output arguments."); 
    } 
    
    /* Check the dimensions of X_IN. */ 
    m = mxGetM(X_IN); 
    n = mxGetN(X_IN);
    if ((MIN(m,n) != 1)) { 
      mexErrMsgTxt("Requires that input is vector."); 
    }
    fftLen = MAX(m,n);
    isComplex = mxIsComplex(X_IN);
    if (!mxIsComplex(X_IN) && !mxIsDouble(X_IN)) {
      mexErrMsgTxt("Requires input double or complex data ."); 
    }
    
    
    InBuffer = mxCreateDoubleMatrix(fftLen, 1, mxCOMPLEX); 
    xRealPtr = mxGetPr(InBuffer);
    xImagPtr = mxGetPi(InBuffer);
    if(isComplex) {
    /* Assign pointers to input data */
      for(idx = 0; idx < fftLen; idx ++){
        xRealPtr[idx] = (mxGetPr(X_IN))[idx];
        xImagPtr[idx] = (mxGetPi(X_IN))[idx];
      }
    }
    else{
      for(idx = 0; idx < fftLen; idx ++){
        xRealPtr[idx] = (mxGetPr(X_IN))[idx];
        xImagPtr[idx] = 0;
      }
    }
    /* Create a matrix for the return data */ 
    Y_OUT = mxCreateDoubleMatrix(fftLen, 1, mxCOMPLEX); 

    yRealPtr = mxGetPr(Y_OUT); 
    yImagPtr = mxGetPi(Y_OUT); 
    r2fftMexTb(fftLen, xRealPtr, xImagPtr, yRealPtr, yImagPtr);
    return;
    
}

