Diff of /mex/comp_dct.c [3e30ee] .. [ba1662]  Maximize  Restore

Switch to side-by-side view

--- a/mex/comp_dct.c
+++ b/mex/comp_dct.c
@@ -13,19 +13,6 @@
 #define COMPLEXARGS
 #define NOCOMPLEXFMTCHANGE
 
-
-static fftw_plan* p_old = 0;
-
-static void fftrealAtExit(void)
-{
-  if(p_old!=0)
-  {
-     fftw_destroy_plan(*p_old);
-     free(p_old);
-  }
-}
-
-
 #endif /* _LTFAT_MEX_FILE */
 
 #define MEX_FILE __BASE_FILE__
@@ -36,204 +23,92 @@
 #include "config.h"
 
 
-// Calling convention:
-//  comp_dct(f,type);
-/*
-FFTW_REDFT00 computes an REDFT00 transform, i.e. a DCT-I. (Logical N=2*(n-1), inverse is FFTW_REDFT00.)
-FFTW_REDFT10 computes an REDFT10 transform, i.e. a DCT-II (sometimes called “the” DCT). (Logical N=2*n, inverse is FFTW_REDFT01.)
-FFTW_REDFT01 computes an REDFT01 transform, i.e. a DCT-III (sometimes called “the” IDCT, being the inverse of DCT-II). (Logical N=2*n, inverse is FFTW_REDFT=10.)
-FFTW_REDFT11 computes an REDFT11 transform, i.e. a DCT-IV. (Logical N=2*n, inverse is FFTW_REDFT11.)
-*/
+static LTFAT_FFTW(plan) LTFAT_NAME(p_old) = 0;
+
+static void LTFAT_NAME(dctMexAtExitFnc)()
+{
+    if(LTFAT_NAME(p_old)!=0)
+    {
+        LTFAT_FFTW(destroy_plan)(LTFAT_NAME(p_old));
+    }
+}
 
 
-void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
+void
+LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],
+                         int nrhs, const mxArray *prhs[] )
 {
-  // Register exit function only once
-  #ifdef LTFAT_DOUBLE
-  if(p_old==0)
-  {
-      mexAtExit(fftrealAtExit);
-  }
-  #endif
+    // Register exit function only once
+    static int atExitFncRegistered = 0;
+    if(!atExitFncRegistered)
+    {
+        LTFAT_NAME(ltfatMexAtExit)(LTFAT_NAME(dctMexAtExitFnc));
+        atExitFncRegistered = 1;
+    }
 
-  mwIndex L, W, N, type;
-  LTFAT_REAL *f_r, *f_i;
-  LTFAT_FFTW(iodim) dims[1], howmanydims[1];
-  LTFAT_FFTW(plan) p;
-  LTFAT_FFTW(r2r_kind) kind[1];
+    LTFAT_REAL *c_r, *c_i;
+    const LTFAT_REAL *f_r, *f_i;
+    dct_kind kind;
 
-  L = mxGetM(prhs[0]);
-  W = mxGetN(prhs[0]);
-  N = 2*L;
-  type = (mwIndex) mxGetScalar(prhs[1]);
+    mwIndex L = mxGetM(prhs[0]);
+    mwIndex W = mxGetN(prhs[0]);
+    mwIndex type = (mwIndex) mxGetScalar(prhs[1]);
 
-  // Copy inputs and get pointers
-  if( mxIsComplex(prhs[0]))
-  {
-     plhs[0] = ltfatCreateMatrix(L, W, LTFAT_MX_CLASSID, mxCOMPLEX);
-     f_i = (LTFAT_REAL*) mxGetPi(plhs[0]);
-     memcpy(f_i,mxGetPi(prhs[0]),W*L*sizeof(LTFAT_REAL));
-  }
-  else
-  {
-     plhs[0] = ltfatCreateMatrix(L, W, LTFAT_MX_CLASSID, mxREAL);
-  }
-  f_r = (LTFAT_REAL*) mxGetPr(plhs[0]);
-  memcpy(f_r,mxGetPr(prhs[0]),W*L*sizeof(LTFAT_REAL));
+// Copy inputs and get pointers
+    if( mxIsComplex(prhs[0]))
+    {
+        f_i =  mxGetImagData(prhs[0]);
+        plhs[0] = ltfatCreateMatrix(L, W, LTFAT_MX_CLASSID, mxCOMPLEX);
+        c_i = mxGetImagData(plhs[0]);
+    }
+    else
+    {
+        plhs[0] = ltfatCreateMatrix(L, W, LTFAT_MX_CLASSID, mxREAL);
+    }
 
-  // Create plan. Copy data from f to cout.
-  dims[0].n = L;
-  dims[0].is = 1;
-  dims[0].os = 1;
+    f_r = mxGetData(prhs[0]);
+    c_r = mxGetData(plhs[0]);
 
-  howmanydims[0].n = W;
-  howmanydims[0].is = L;
-  howmanydims[0].os = L;
-
-  LTFAT_REAL sqrt2 = (LTFAT_REAL) sqrt(2.0);
-  LTFAT_REAL postScale = (LTFAT_REAL) 1.0/sqrt2;
-  LTFAT_REAL scale = (LTFAT_REAL) sqrt2*(1.0/(double)N)*sqrt((double)L);
-
-  // Re-allocate and prescale input
-  if(type==1||type==3)
-  {
-     for(mwIndex ii=0;ii<W;ii++)
-     {
-	    f_r[ii*L] *= sqrt2;
-     }
-     if(mxIsComplex(prhs[0]))
-     {
-        for(mwIndex ii=0;ii<W;ii++)
-        {
-           f_i[ii*L] *= sqrt2;
-        }
-     }
-  }
-
-  switch(type)
-  {
-	 case 1:
-		N -= 2;
-		for(mwIndex ii=0;ii<W;ii++)
-        {
-	       f_r[(ii+1)*L-1] *= sqrt2;
-        }
-
-        if(mxIsComplex(prhs[0]))
-        {
-           for(mwIndex ii=0;ii<W;ii++)
-           {
-              f_i[(ii+1)*L-1] *= sqrt2;
-           }
-        }
-
-		scale = (LTFAT_REAL) sqrt2*(1.0/((double)N))*sqrt((double)L-1);
-        kind[0] = FFTW_REDFT00;
-     break;
-	 case 2:
-        kind[0] = FFTW_REDFT10;
-     break;
-	 case 3:
-        kind[0] = FFTW_REDFT01;
-     break;
-	 case 4:
-        kind[0] = FFTW_REDFT11;
-     break;
-     default:
-		 mexErrMsgTxt("Unknown type.");
-  }
+    switch(type)
+    {
+    case 1:
+        kind = DCTI;
+        break;
+    case 2:
+        kind = DCTII;
+        break;
+    case 3:
+        kind = DCTIII;
+        break;
+    case 4:
+        kind = DCTIV;
+        break;
+    default:
+        mexErrMsgTxt("Unknown type.");
+    }
 
 
 
-  // The calling prototype
-  //fftw_plan fftw_plan_guru_split_dft_r2c(
-  //        int rank, const fftw_iodim *dims,
-  //        int howmany_rank, const fftw_iodim *howmany_dims,
-  //        double *in, double *ro, double *io,
-  //        unsigned flags);
+    LTFAT_FFTW(plan) p = LTFAT_NAME(dct_init)( c_r, L, W, kind);
+    /*
+    The old plan is freed after the new one is cretaed.
+    According to the FFTW doc. creating new plan is quick as long as there 
+    already exists a plan for the same length.
+    */
 
 
-  p = LTFAT_FFTW(plan_guru_r2r)(1, dims,
-				   1, howmanydims,
-				   f_r, f_r,
-				   kind,
-				   FFTW_OPTITYPE);
-  /*
-  FFTW documentation qote http://www.fftw.org/fftw3_doc/New_002darray-Execute-Functions.html#New_002darray-Execute-Functions:
-  ...
-  creating a new plan is quick once one exists for a given size
-  ...
-  so why not to store the old plan..
-  */
+    LTFAT_NAME(dctMexAtExitFnc)();
+    LTFAT_NAME(p_old) = p;
+
+    LTFAT_NAME(dct_plan)(f_r,L,W,kind,c_r,p);
+
+    if( mxIsComplex(prhs[0]))
+    {
+        LTFAT_NAME(dct_plan)(f_i,L,W,kind,c_i,p);
+    }
 
 
-  if(p_old!=0)
-  {
-    fftw_destroy_plan(*p_old);
-    free(p_old);
-  }
-  p_old = malloc(sizeof(p));
-  memcpy(p_old,&p,sizeof(p));
-
-
-  // Real FFT.
-  LTFAT_FFTW(execute)(p);
-
-
-  // Do the normalization
-  for(int ii=0;ii<L*W;ii++)
-  {
-	  f_r[ii] *= scale;
-  }
-
-  if(type==1||type==2)
-  {
-     // Scale DC component
-     for(int ii=0;ii<W;ii++)
-     {
-	    f_r[ii*L] *= postScale;
-     }
-  }
-
-  if(type==1)
-  {
-     // Scale AC component
-     for(int ii=0;ii<W;ii++)
-     {
-	    f_r[(ii+1)*L-1] *= postScale;
-     }
-  }
-
-  // If the input is complex, process the imaginary part
-  if(mxIsComplex(prhs[0]))
-  {
-	  LTFAT_FFTW(execute_r2r)(p,f_i,f_i);
-	    // Do the normalization
-      for(int ii=0;ii<L*W;ii++)
-      {
-	     f_i[ii] *= scale;
-      }
-      if(type==1||type==2)
-      {
-         // Scale DC component
-         for(int ii=0;ii<W;ii++)
-         {
-	        f_i[ii*L] *= postScale;
-         }
-      }
-      if(type==1)
-      {
-         // Scale AC component
-         for(int ii=0;ii<W;ii++)
-         {
-	        f_i[(ii+1)*L-1] *= postScale;
-         }
-      }
-  }
-
-  // LTFAT_FFTW(destroy_plan)(p);
-
-  return;
+    return;
 }
 #endif
 

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks