Commit [f1ff16] Maximize Restore History

Closing #22 FEATURE: Complete C implementation of Wilson and MDCT transform.

Rewrite of the dgt comp. routines.

prusa prusa 2013-11-29

<< < 1 2 3 4 5 > >> (Page 4 of 5)
added comp/comp_isepdgt.m
added comp/comp_isepdgtreal.m
added comp/comp_sepdgtreal.m
added src/idgt.c
changed .gitignore
changed comp
changed comp/comp_dgt.m
changed comp/comp_dgtreal.m
changed comp/comp_dwilt.m
changed comp/comp_dwiltiii.m
changed comp/comp_idgt.m
changed comp/comp_idgt_fac.m
changed comp/comp_idgt_fb.m
changed comp/comp_idgtreal.m
changed comp/comp_idgtreal_fac.m
changed comp/comp_idgtreal_fb.m
changed comp/comp_idwilt.m
changed comp/comp_idwiltiii.m
changed mex
changed mex/filedefs.mk
changed mex/ltfat_mex_template_helper.h
changed oct
changed oct/Makefile_mingwoct
changed oct/comp_atrousfilterbank_td.cc
changed oct/comp_col2diag.cc
changed oct/comp_dct.cc
changed oct/comp_fftreal.cc
changed oct/comp_ifftreal.cc
changed oct/comp_wfac.cc
changed oct/config.h
changed oct/ltfat_oct_template_helper.h
changed src
changed src/c-safe-memalloc.c
changed src/dgt.c
changed src/dgt_fac.c
changed src/dgt_fb.c
changed src/dgt_long.h
changed src/dgt_multi.c
changed src/dgt_ola.c
changed src/dgt_shear.c
changed src/dgt_shearola.c
changed src/dgt_walnut.c
changed src/dgtreal_fac.c
changed src/dwilt.c
changed src/filedefs.mk
changed src/filterbank.c
changed src/gabdual.c
changed src/gabtight.c
changed src/gabtight_fac.c
changed src/goertzel.c
changed src/idgt_fac.c
changed src/idgt_fb.c
changed src/iwfac.c
changed src/ltfat.h
changed src/ltfat_typecomplexindependent.h
changed src/ltfat_typeindependent.h
changed src/reassign.c
changed src/wavelets.c
changed testing
changed testing/test_dgt.m
changed testing/test_dwilt.m
changed testing/test_wmdct.m
copied comp/comp_dwilt_fb.m -> comp/comp_idgtreal_long.m
copied comp/comp_dwilt_long.m -> comp/comp_idgt_long.m
copied mex/comp_dgt_fb.c -> oct/comp_sepdgtreal.cc
copied mex/comp_dgt_long.c -> oct/comp_isepdgtreal.cc
copied mex/comp_dgtreal_fb.c -> mex/comp_isepdgt.c
copied mex/comp_dgtreal_long.c -> mex/comp_sepdgtreal.c
copied mex/comp_dwilt_long.c -> mex/comp_dwilt.c
copied mex/comp_idgt_fac.c -> oct/comp_isepdgt.cc
copied mex/comp_idgt_fb.c -> src/idwilt.c
copied mex/comp_idgtreal_fac.c -> mex/comp_sepdgt.c
copied mex/comp_idgtreal_fb.c -> mex/comp_isepdgtreal.c
copied oct/comp_dgt_fb.cc -> oct/comp_dwilt.cc
copied oct/comp_dgt_long.cc -> oct/comp_idwilt.cc
copied oct/comp_dgtreal_fb.cc -> oct/comp_idwiltiii.cc
copied oct/comp_dgtreal_long.cc -> oct/comp_sepdgt.cc
copied oct/comp_dwilt_long.cc -> src/wmdct.c
copied oct/comp_idgt_fac.cc -> src/iwmdct.c
copied oct/comp_idgt_fb.cc -> mex/comp_dwiltiii.c
copied oct/comp_idgtreal_fac.cc -> mex/comp_idwilt.c
copied oct/comp_idgtreal_fb.cc -> mex/comp_idwiltiii.c
copied src/tfutil.c -> comp/comp_sepdgt.m
comp/comp_isepdgt.m Diff Switch to side-by-side view
Loading...
comp/comp_isepdgtreal.m Diff Switch to side-by-side view
Loading...
comp/comp_sepdgtreal.m Diff Switch to side-by-side view
Loading...
src/idgt.c Diff Switch to side-by-side view
Loading...
.gitignore Diff Switch to side-by-side view
Loading...
comp
Directory.
comp/comp_dgt.m Diff Switch to side-by-side view
Loading...
comp/comp_dgtreal.m Diff Switch to side-by-side view
Loading...
comp/comp_dwilt.m Diff Switch to side-by-side view
Loading...
comp/comp_dwiltiii.m Diff Switch to side-by-side view
Loading...
comp/comp_idgt.m Diff Switch to side-by-side view
Loading...
comp/comp_idgt_fac.m Diff Switch to side-by-side view
Loading...
comp/comp_idgt_fb.m Diff Switch to side-by-side view
Loading...
comp/comp_idgtreal.m Diff Switch to side-by-side view
Loading...
comp/comp_idgtreal_fac.m Diff Switch to side-by-side view
Loading...
comp/comp_idgtreal_fb.m Diff Switch to side-by-side view
Loading...
comp/comp_idwilt.m Diff Switch to side-by-side view
Loading...
comp/comp_idwiltiii.m Diff Switch to side-by-side view
Loading...
mex
Directory.
mex/filedefs.mk Diff Switch to side-by-side view
Loading...
mex/ltfat_mex_template_helper.h Diff Switch to side-by-side view
Loading...
oct
Directory.
oct/Makefile_mingwoct Diff Switch to side-by-side view
Loading...
oct/comp_atrousfilterbank_td.cc Diff Switch to side-by-side view
Loading...
oct/comp_col2diag.cc Diff Switch to side-by-side view
Loading...
oct/comp_dct.cc Diff Switch to side-by-side view
Loading...
oct/comp_fftreal.cc Diff Switch to side-by-side view
Loading...
oct/comp_ifftreal.cc Diff Switch to side-by-side view
Loading...
oct/comp_wfac.cc Diff Switch to side-by-side view
Loading...
oct/config.h Diff Switch to side-by-side view
Loading...
oct/ltfat_oct_template_helper.h Diff Switch to side-by-side view
Loading...
src
Directory.
src/c-safe-memalloc.c Diff Switch to side-by-side view
Loading...
src/dgt.c Diff Switch to side-by-side view
Loading...
src/dgt_fac.c Diff Switch to side-by-side view
Loading...
src/dgt_fb.c Diff Switch to side-by-side view
Loading...
src/dgt_long.h Diff Switch to side-by-side view
Loading...
src/dgt_multi.c Diff Switch to side-by-side view
Loading...
src/dgt_ola.c Diff Switch to side-by-side view
Loading...
src/dgt_shear.c Diff Switch to side-by-side view
Loading...
src/dgt_shearola.c Diff Switch to side-by-side view
Loading...
src/dgt_walnut.c Diff Switch to side-by-side view
Loading...
src/dgtreal_fac.c Diff Switch to side-by-side view
Loading...
src/dwilt.c Diff Switch to side-by-side view
Loading...
src/filedefs.mk Diff Switch to side-by-side view
Loading...
src/filterbank.c Diff Switch to side-by-side view
Loading...
src/gabdual.c Diff Switch to side-by-side view
Loading...
src/gabtight.c Diff Switch to side-by-side view
Loading...
src/gabtight_fac.c Diff Switch to side-by-side view
Loading...
src/goertzel.c Diff Switch to side-by-side view
Loading...
src/idgt_fac.c Diff Switch to side-by-side view
Loading...
src/idgt_fb.c Diff Switch to side-by-side view
Loading...
src/iwfac.c Diff Switch to side-by-side view
Loading...
src/ltfat.h Diff Switch to side-by-side view
Loading...
src/ltfat_typecomplexindependent.h Diff Switch to side-by-side view
Loading...
src/ltfat_typeindependent.h Diff Switch to side-by-side view
Loading...
src/reassign.c Diff Switch to side-by-side view
Loading...
src/wavelets.c Diff Switch to side-by-side view
Loading...
testing
Directory.
testing/test_dgt.m Diff Switch to side-by-side view
Loading...
testing/test_dwilt.m Diff Switch to side-by-side view
Loading...
testing/test_wmdct.m Diff Switch to side-by-side view
Loading...
comp/comp_dwilt_fb.m to comp/comp_idgtreal_long.m
--- a/comp/comp_dwilt_fb.m
+++ b/comp/comp_idgtreal_long.m
@@ -1,68 +1,35 @@
-function [coef]=comp_dwilt_fb(f,g,M,L)
-%COMP_DWILT_FB  Compute Discrete Wilson transform.
-%   
+function f=comp_idgtreal_long(coef,g,L,a,M)
+%COMP_IDGTREAL_FAC  Full-window factorization of a Gabor matrix assuming.
+%   Usage:  f=comp_idgtreal_long(c,g,L,a,M)
+%
+%   Input parameters:
+%         c     : M x N x W array of coefficients.
+%         g     : window (from facgabm).
+%         a     : Length of time shift.
+%         M     : Number of frequency shifts.
+%   Output parameters:
+%         f     : Reconstructed signal.
+%
+%   Do not call this function directly, use IDGT.
+%   This function does not check input parameters!
+%
+%   If input is a matrix, the transformation is applied to
+%   each column.
+%
+%   This function does not handle multidimensional data, take care before
+%   you call it.
+%
+%   References: so07-2 st98-8
 
-a=M;
-N=L/a;
-W=size(f,2);
+%   AUTHOR : Peter L. Søndergaard.
+%   TESTING: OK
+%   REFERENCE: OK
 
-
-coef=zeros(2*M,N/2,W,assert_classname(f,g));
-
-if (isreal(f) && isreal(g))
-
-  coef2=comp_dgt_fb(f,g,a,2*M);
-
-  % If the input coefficients are real, the calculations can be
-  % be simplified. The complex case code also works for the real case.
-
-  % Unmodulated case.
-  coef(1,:,:)=coef2(1,1:2:N,:);
-
-  % cosine, first column.
-  coef(3:2:M,:,:)=sqrt(2)*real(coef2(3:2:M,1:2:N,:));
-  
-  % sine, second column
-  coef(M+3:2:2*M,:,:)=-sqrt(2)*imag(coef2(3:2:M,2:2:N,:));
-  
-  % sine, first column.
-  coef(2:2:M,:,:)=-sqrt(2)*imag(coef2(2:2:M,1:2:N,:));
-  
-  % cosine, second column
-  coef(M+2:2:2*M,:,:)=sqrt(2)*real(coef2(2:2:M,2:2:N,:));
-
-  % Nyquest case
-  if mod(M,2)==0
-    coef(M+1,:,:) = coef2(M+1,1:2:N,:);
-  else
-    coef(M+1,:,:) = coef2(M+1,2:2:N,:);
-  end;
-
-
-else
-  % Complex valued case
-  
-  coef2=comp_dgt_fb(f,g,a,2*M);
-  
-  % Unmodulated case.
-  coef(1,:,:)=coef2(1,1:2:N,:);
-
-  % odd value of m
-  coef(2:2:M,:,:)=1/sqrt(2)*i*(coef2(2:2:M,1:2:N,:)-coef2(2*M:-2:M+2,1:2:N,:));
-  coef(M+2:2:2*M,:,:)=1/sqrt(2)*(coef2(2:2:M,2:2:N,:)+coef2(2*M:-2:M+2,2:2:N,:));
-
-  % even value of m
-  coef(3:2:M,:,:)=1/sqrt(2)*(coef2(3:2:M,1:2:N,:)+coef2(2*M-1:-2:M+2,1:2:N,:));
-  coef(M+3:2:2*M,:,:)=1/sqrt(2)*i*(coef2(3:2:M,2:2:N,:)-coef2(2*M-1:-2:M+2,2:2:N,:));
-
-  % Nyquest case
-  if mod(M,2)==0
-    coef(M+1,:,:) = coef2(M+1,1:2:N,:);
-  else
-    coef(M+1,:,:) = coef2(M+1,2:2:N,:);
-  end;
-
-end;
+% Get the factorization of the window.
+gf = comp_wfac(g,a,M);      
+        
+% Call the computational subroutine.
+f = comp_idgtreal_fac(coef,gf,L,a,M);
 
 
 
comp/comp_dwilt_long.m to comp/comp_idgt_long.m
--- a/comp/comp_dwilt_long.m
+++ b/comp/comp_idgt_long.m
@@ -1,70 +1,30 @@
-function [coef]=comp_dwilt_long(f,g,M,L)
-%COMP_DWILT_LONG  Compute Discrete Wilson transform.
-%   
+function f=comp_idgt_long(coef,g,L,a,M)
+%COMP_IDGT_FAC  Full-window factorization of a Gabor matrix.
+%   Usage:  f=comp_idgt_long(c,g,L,a,M)
+%
+%   Input parameters:
+%         c     : M x N x W array of coefficients.
+%         g     : Window.
+%         a     : Length of time shift.
+%         M     : Number of frequency shifts.
+%   Output parameters:
+%         f     : Reconstructed signal.
+%
+%   Do not call this function directly, use IDGT.
+%   This function does not check input parameters!
+%
+%   If input is a matrix, the transformation is applied to
+%   each column.
+%
+%   This function does not handle multidimensional data, take care before
+%   you call it.
+%
+%   References: so07-2 st98-8
 
-a=M;
-N=L/a;
-W=size(f,2);
+%   AUTHOR : Peter L. Søndergaard.
 
+% Get the factorization of the window.
+gf = comp_wfac(g,a,M);      
 
-coef=zeros(2*M,N/2,W,assert_classname(f,g));
-
-if (isreal(f) && isreal(g))
-
-  coef2=comp_dgt_long(f,g,a,2*M);
-
-  % If the input coefficients are real, the calculations can be
-  % be simplified. The complex case code also works for the real case.
-
-  % Unmodulated case.
-  coef(1,:,:)=coef2(1,1:2:N,:);
-
-  % cosine, first column.
-  coef(3:2:M,:,:)=sqrt(2)*real(coef2(3:2:M,1:2:N,:));
-  
-  % sine, second column
-  coef(M+3:2:2*M,:,:)=-sqrt(2)*imag(coef2(3:2:M,2:2:N,:));
-  
-  % sine, first column.
-  coef(2:2:M,:,:)=-sqrt(2)*imag(coef2(2:2:M,1:2:N,:));
-  
-  % cosine, second column
-  coef(M+2:2:2*M,:,:)=sqrt(2)*real(coef2(2:2:M,2:2:N,:));
-
-  % Nyquest case
-  if mod(M,2)==0
-    coef(M+1,:,:) = coef2(M+1,1:2:N,:);
-  else
-    coef(M+1,:,:) = coef2(M+1,2:2:N,:);
-  end;
-
-
-else
-  % Complex valued case
-  
-  coef2=comp_dgt_long(f,g,a,2*M);
-  
-  % Unmodulated case.
-  coef(1,:,:)=coef2(1,1:2:N,:);
-
-  % odd value of m
-  coef(2:2:M,:,:)=1/sqrt(2)*i*(coef2(2:2:M,1:2:N,:)-coef2(2*M:-2:M+2,1:2:N,:));
-  coef(M+2:2:2*M,:,:)=1/sqrt(2)*(coef2(2:2:M,2:2:N,:)+coef2(2*M:-2:M+2,2:2:N,:));
-
-  % even value of m
-  coef(3:2:M,:,:)=1/sqrt(2)*(coef2(3:2:M,1:2:N,:)+coef2(2*M-1:-2:M+2,1:2:N,:));
-  coef(M+3:2:2*M,:,:)=1/sqrt(2)*i*(coef2(3:2:M,2:2:N,:)-coef2(2*M-1:-2:M+2,2:2:N,:));
-
-  % Nyquest case
-  if mod(M,2)==0
-    coef(M+1,:,:) = coef2(M+1,1:2:N,:);
-  else
-    coef(M+1,:,:) = coef2(M+1,2:2:N,:);
-  end;
-
-
- 
-end;
-
-
-
+% Call the computational subroutine.
+f  = comp_idgt_fac(coef,gf,L,a,M);
mex/comp_dgt_fb.c to oct/comp_sepdgtreal.cc
--- a/mex/comp_dgt_fb.c
+++ b/oct/comp_sepdgtreal.cc
@@ -1,125 +1,78 @@
-#ifndef _LTFAT_MEX_FILE
-#define _LTFAT_MEX_FILE
-
-#define ISNARGINEQ 4
 #define TYPEDEPARGS 0, 1
 #define SINGLEARGS
-#define COMPLEXINDEPENDENT
-
-#endif // _LTFAT_MEX_FILE - INCLUDED ONCE
-
-#define MEX_FILE __BASE_FILE__
-#include "ltfat_mex_template_helper.h"
-
-#if defined(LTFAT_SINGLE) || defined(LTFAT_DOUBLE)
-#include "ltfat_types.h"
-
-// Calling convention:
-//  comp_dgt_fb(f,g,a,M);
-
-void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
-{
-   int L  = mxGetM(prhs[0]);
-   int W  = mxGetN(prhs[0]);
-   int gl = mxGetM(prhs[1]);
-
-   int a=(int)mxGetScalar(prhs[2]);
-   int M=(int)mxGetScalar(prhs[3]);
-
-   int N=L/a;
-
-   mwSize ndim=3;
-   if (W==1)
-   {
-      ndim=2;
-   }
-
-   mwSize dims[3] = {M, N, W};
-   plhs[0] = ltfatCreateNdimArray(ndim,dims,LTFAT_MX_CLASSID,mxCOMPLEX);
-   const LTFAT_TYPE* f_combined = (const LTFAT_TYPE*) mxGetData(prhs[0]);
-   const LTFAT_TYPE* g_combined = (const LTFAT_TYPE*) mxGetData(prhs[1]);
-   LTFAT_COMPLEX* out_combined = (LTFAT_COMPLEX*) mxGetData(plhs[0]);
-
-   LTFAT_NAME(dgt_fb)(f_combined, g_combined,L,gl,W,a,M,out_combined);
-}
-#endif
+#define REALARGS
+#define OCTFILENAME comp_sepdgtreal // change to filename
+#define OCTFILEHELP "This function calls the C-library\n c=comp_dgtreal_fb(f,g,a,M,boundary);\n Yeah."
 
 
+#include "ltfat_oct_template_helper.h"
+// octave_idx_type is 32 or 64 bit signed integer
 /*
-#include "mex.h"
-#include "config.h"
-#include "ltfat.h"
-#include "ltfat-mex-helper.h"
+   dgtreal_fb forwarders
+   */
 
-// Calling convention:
-// comp_dgt_fb(f,g,a,M);
+static inline void
+fwd_dgtreal_fb(const double *f, const double *g,
+               const octave_idx_type L, const octave_idx_type gl,
+               const octave_idx_type W, const octave_idx_type a,
+               const octave_idx_type M, Complex *cout)
+{
+   dgtreal_fb_d(f,g,L,gl,W,a,M,reinterpret_cast<double _Complex*>(cout));
+}
 
-void mexFunction( int nlhs, mxArray *plhs[],
-		  int nrhs, const mxArray *prhs[] )
+static inline void
+fwd_dgtreal_fb(const float *f, const float *g,
+               const octave_idx_type L, const octave_idx_type gl,
+               const octave_idx_type W, const octave_idx_type a,
+               const octave_idx_type M, FloatComplex *cout)
+{
+   dgtreal_fb_s(f,g,L,gl,W,a,M,reinterpret_cast<float _Complex*>(cout));
+}
 
+static inline void
+fwd_dgtreal_long(const double *f, const double *g, const octave_idx_type L,
+                 const octave_idx_type W, const octave_idx_type a,
+                const octave_idx_type M, Complex *cout)
 {
-   int L, gl,W, a, M, N;
-   mwSize ndim;
-   mwSize dims[3];
+   dgtreal_long_d(f,g,L,W,a,M,reinterpret_cast<double _Complex*>(cout));
+}
 
-   ltfat_complex *f_combined, *g_combined;
-   ltfat_complex *out_combined;
+static inline void
+fwd_dgtreal_long(const float *f, const float *g, const octave_idx_type L,
+                 const octave_idx_type W, const octave_idx_type a,
+                 const octave_idx_type M, FloatComplex *cout)
+{
+   dgtreal_long_s(f,g,L,W,a,M,reinterpret_cast<float _Complex*>(cout));
+}
+template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
+octave_value_list octFunction(const octave_value_list& args, int nargout)
+{
+   DEBUGINFO;
+   const octave_idx_type a = args(2).int_value();
+   const octave_idx_type M = args(3).int_value();
+   const octave_idx_type M2 = M/2+1;
 
-   // Get matrix dimensions.
-   L  = mxGetM(prhs[0]);
-   W  = mxGetN(prhs[0]);
-   gl = mxGetM(prhs[1]);
+   MArray<LTFAT_TYPE> f = ltfatOctArray<LTFAT_TYPE>(args(0));
+   MArray<LTFAT_TYPE> g = ltfatOctArray<LTFAT_TYPE>(args(1));
+   const octave_idx_type L  = f.rows();
+   const octave_idx_type W  = f.columns();
+   const octave_idx_type gl = g.rows();
+   const octave_idx_type N = L/a;
 
-   a=(int)mxGetScalar(prhs[2]);
-   M=(int)mxGetScalar(prhs[3]);
+   dim_vector dims_out(M2,N,W);  
+   dims_out.chop_trailing_singletons();
 
-   N=L/a;
+   MArray<LTFAT_COMPLEX> cout(dims_out); 
+   cout.fill(0);
 
-   dims[0]=M;
-   dims[1]=N;
-   dims[2]=W;
-   ndim=3;
-   if (W==1)
+   if(gl<L)
    {
-      ndim=2;
-   }
-
-   out_combined=mxMalloc(M*N*W*sizeof(ltfat_complex));
-
-   // Copy the data.
-
-   if (mxIsComplex(prhs[0]) || mxIsComplex(prhs[1]))
-   {
-      f_combined=mxMalloc(L*W*sizeof(ltfat_complex));
-      split2combined(L*W, prhs[0], f_combined);
-
-      g_combined=mxMalloc(L*2*sizeof(ltfat_complex));
-      split2combined(gl,  prhs[1], g_combined);
-
-      dgt_fb((const ltfat_complex*)f_combined,(const ltfat_complex*)g_combined,
-	     L,gl,W,a,M,
-	     (ltfat_complex*)out_combined);
-
-      mxFree(f_combined);
-      mxFree(g_combined);
-
+      fwd_dgtreal_fb(f.data(),g.data(),L,gl,W,a,M,cout.fortran_vec());
    }
    else
    {
-      dgt_fb_r((const double*)mxGetPr(prhs[0]),
-	       (const double*)mxGetPr(prhs[1]),
-	       L,gl,W,a,M,
-	       (ltfat_complex*)out_combined);
+      fwd_dgtreal_long(f.data(),g.data(),L,W,a,M,cout.fortran_vec());
    }
 
-   plhs[0] = mxCreateNumericArray(ndim,dims,mxDOUBLE_CLASS,mxCOMPLEX);
-
-   combined2split(M*N*W, (const ltfat_complex*)out_combined, plhs[0]);
-
-   mxFree(out_combined);
-
-   return;
-
+   return octave_value(cout);
 }
-
-*/
mex/comp_dgt_long.c to oct/comp_isepdgtreal.cc
--- a/mex/comp_dgt_long.c
+++ b/oct/comp_isepdgtreal.cc
@@ -1,125 +1,81 @@
-#ifndef _LTFAT_MEX_FILE
-#define _LTFAT_MEX_FILE
-
-#define ISNARGINEQ 4
-#define TYPEDEPARGS 0, 1
+#define TYPEDEPARGS 0
 #define SINGLEARGS
+#define MATCHEDARGS 1
 #define COMPLEXARGS
-
-#endif // _LTFAT_MEX_FILE - INCLUDED ONCE
-
-#define MEX_FILE __BASE_FILE__
-#include "ltfat_mex_template_helper.h"
-
-#if defined(LTFAT_SINGLE) || defined(LTFAT_DOUBLE)
-#include "ltfat_types.h"
-
-// Calling convention:
-// comp_dgt_long(f,g,a,M);
-
-void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
-{
-   int L, W, a, M, N;
-   mwSize ndim;
-   mwSize dims[3];
-
-   // Get matrix dimensions.
-   L = mxGetM(prhs[0]);
-   W = mxGetN(prhs[0]);
-
-   a=(int)mxGetScalar(prhs[2]);
-   M=(int)mxGetScalar(prhs[3]);
-
-   N=L/a;
-
-   dims[0]=M;
-   dims[1]=N;
-   dims[2]=W;
-   ndim=3;
-   if (W==1)
-   {
-      ndim=2;
-   }
-
-   plhs[0] = ltfatCreateNdimArray(ndim,dims,LTFAT_MX_CLASSID,mxCOMPLEX);
-   const LTFAT_COMPLEX* f_combined = (const LTFAT_COMPLEX*) mxGetData(prhs[0]);
-   const LTFAT_COMPLEX* g_combined = (const LTFAT_COMPLEX*) mxGetData(prhs[1]);
-   LTFAT_COMPLEX* out_combined = (LTFAT_COMPLEX*) mxGetData(plhs[0]);
-
-   LTFAT_NAME(dgt_long)(f_combined,g_combined,L,W,a,M,out_combined);
-
-   return;
-}
-#endif
+#define OCTFILENAME comp_isepdgtreal // change to filename
+#define OCTFILEHELP "This function calls the C-library\n\
+                    c=comp_idgtereal_fb(c,g,L,a,M);\n\
+                    Yeah.\n"
 
 
-/*
-#include "mex.h"
-#include "config.h"
-#include "ltfat.h"
-#include "ltfat-mex-helper.h"
+#include "ltfat_oct_template_helper.h"
+// octave_idx_type is 32 or 64 bit signed integer
 
-// Calling convention:
-// comp_dgt_long(f,g,a,M);
-
-
-
-void dothedouble(int nlhs, mxArray *plhs[],
-		 int nrhs, const mxArray *prhs[],
-		 int L, int W, int a, int M, int N)
+static inline void
+fwd_idgtreal_fb(const Complex *coef, const double *gf,
+                const octave_idx_type L, const octave_idx_type gl,
+                const octave_idx_type W, const octave_idx_type a,
+                const octave_idx_type M, double *f)
 {
-
-
+   idgtreal_fb_d(reinterpret_cast<const double _Complex*>(coef),
+                 gf,L,gl,W,a,M,f);
 }
 
-void mexFunction( int nlhs, mxArray *plhs[],
-		  int nrhs, const mxArray *prhs[] )
+static inline void
+fwd_idgtreal_fb(const FloatComplex *coef, const float *gf,
+                const octave_idx_type L, const octave_idx_type gl,
+                const octave_idx_type W, const octave_idx_type a,
+                const octave_idx_type M, float *f)
+{
+   idgtreal_fb_s(reinterpret_cast<const float _Complex*>(coef),
+                 gf,L,gl,W,a,M,f);
+}
 
+static inline void
+fwd_idgtreal_long(const Complex *coef, const double *gf,
+                  const octave_idx_type L, const octave_idx_type W,
+                  const octave_idx_type a, const octave_idx_type M,
+                  double *f)
 {
-   int L, W, a, M, N;
-   mwSize ndim;
-   mwSize dims[3];
+   idgtreal_long_d(reinterpret_cast<const double _Complex*>(coef),
+                   gf,L,W,a,M,f);
+}
 
-   ltfat_complex *f_combined, *g_combined, *out_combined;
+static inline void
+fwd_idgtreal_long(const FloatComplex *coef, const float *gf,
+                  const octave_idx_type L, const octave_idx_type W, 
+                  const octave_idx_type a, const octave_idx_type M,
+                  float *f)
+{
+   idgtreal_long_s(reinterpret_cast<const float _Complex*>(coef),
+                   gf,L,W,a,M,f);
+}
 
-   // Get matrix dimensions.
-   L = mxGetM(prhs[0]);
-   W = mxGetN(prhs[0]);
+template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
+octave_value_list octFunction(const octave_value_list& args, int nargout)
+{
+   DEBUGINFO;
 
-   a=(int)mxGetScalar(prhs[2]);
-   M=(int)mxGetScalar(prhs[3]);
+   MArray<LTFAT_TYPE> coef = ltfatOctArray<LTFAT_TYPE>(args(0));
+   MArray<LTFAT_REAL> gf = ltfatOctArray<LTFAT_REAL>(args(1));
+   const octave_idx_type L = args(2).int_value();
+   const octave_idx_type a = args(3).int_value();
+   const octave_idx_type M = args(4).int_value();
+   const octave_idx_type N = L/a;
+   const octave_idx_type M2 = M/2 +1;
+   const octave_idx_type W = coef.nelem()/(N*M2);
+   const octave_idx_type gl = gf.rows();
 
-   N=L/a;
-
-   dims[0]=M;
-   dims[1]=N;
-   dims[2]=W;
-   ndim=3;
-   if (W==1)
+   MArray<LTFAT_REAL> f(dim_vector(L,W)); 
+   f.fill(0);
+   
+   if(gl<L)
    {
-      ndim=2;
+      fwd_idgtreal_fb(coef.data(),gf.data(),L,gl,W,a,M,f.fortran_vec());
    }
-
-   f_combined=(ltfat_complex*)mxMalloc(L*W*sizeof(ltfat_complex));
-   g_combined=mxMalloc(L*sizeof(ltfat_complex));
-   out_combined=(ltfat_complex*)mxMalloc(M*N*W*sizeof(ltfat_complex));
-
-   split2combined(L*W, prhs[0], f_combined);
-   split2combined(L, prhs[1], g_combined);
-
-   dgt_long((const ltfat_complex*)f_combined,(const ltfat_complex*)g_combined,
-      L,W,a,M,out_combined);
-
-   mxFree(f_combined);
-   mxFree(g_combined);
-
-   plhs[0] = mxCreateNumericArray(ndim,dims,mxDOUBLE_CLASS,mxCOMPLEX);
-
-   combined2split(M*N*W, (const ltfat_complex*)out_combined, plhs[0]);
-
-   mxFree(out_combined);
-   return;
-
+   else
+   {
+      fwd_idgtreal_long(coef.data(),gf.data(),L,W,a,M,f.fortran_vec());
+   }
+   return octave_value(f);
 }
-*/
-
mex/comp_dgtreal_fb.c to mex/comp_isepdgt.c
--- a/mex/comp_dgtreal_fb.c
+++ b/mex/comp_isepdgt.c
@@ -1,10 +1,10 @@
 #ifndef _LTFAT_MEX_FILE
 #define _LTFAT_MEX_FILE
 
-#define ISNARGINEQ 4
+#define ISNARGINEQ 5
 #define TYPEDEPARGS 0, 1
 #define SINGLEARGS
-#define REALARGS
+#define COMPLEXARGS
 
 #endif // _LTFAT_MEX_FILE - INCLUDED ONCE
 
@@ -15,82 +15,41 @@
 #include "ltfat_types.h"
 
 // Calling convention:
-// comp_dgtreal_fb(f,g,a,M);
+//  comp_isepdgt(coef,g,L,a,M);
 
 void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
 {
-   int L, gl,W, a, M, N, M2;
-
+   int L, W, a, M, N, gl;
    // Get matrix dimensions.
-   L  = mxGetM(prhs[0]);
-   W  = mxGetN(prhs[0]);
-   gl = mxGetM(prhs[1]);
-
-   a=(int)mxGetScalar(prhs[2]);
-   M=(int)mxGetScalar(prhs[3]);
-   M2=M/2+1;
-
+   L=(int)mxGetScalar(prhs[2]);
+   a=(int)mxGetScalar(prhs[3]);
+   M=(int)mxGetScalar(prhs[4]);
    N=L/a;
 
-   plhs[0] = ltfatCreateMatrix(M2, N*W,LTFAT_MX_CLASSID,mxCOMPLEX);
-   const LTFAT_REAL * f = (const LTFAT_REAL *) mxGetData(prhs[0]);
-   const LTFAT_REAL * g = (const LTFAT_REAL *) mxGetData(prhs[1]);
-   LTFAT_COMPLEX* out_combined = (LTFAT_COMPLEX*) mxGetData(plhs[0]);
+   gl = mxGetM(prhs[1]);
+   W  = mxGetNumberOfElements(prhs[0])/(M*N);
 
-   LTFAT_NAME(dgtreal_fb)(f,g,L,gl,W,a,M,out_combined);
+   plhs[0] = ltfatCreateMatrix(L, W,LTFAT_MX_CLASSID, mxCOMPLEX);
+   const LTFAT_COMPLEX* c_combined = (const LTFAT_COMPLEX*) mxGetData(prhs[0]);
+   const LTFAT_COMPLEX* g_combined = (const LTFAT_COMPLEX*) mxGetData(prhs[1]);
+   LTFAT_COMPLEX* f_combined = (LTFAT_COMPLEX*) mxGetData(plhs[0]);
 
+   if(gl<L)
+   {
+      LTFAT_NAME(idgt_fb)(c_combined,g_combined, L,gl,W,a,M, f_combined);
+   }
+   else
+   {
+      LTFAT_NAME(idgt_long)(c_combined,g_combined, L,W,a,M, f_combined); 
+   }
+  /* #else
+   NOT CALLING idgt_fb_r:
+   TO DO: Do it better.
+   LTFAT_NAME(idgt_fb_r)((const LTFAT_REAL (*)[2])c_combined,
+                         g_combined,
+                         L,gl,W,a,M,(LTFAT_REAL (*)[2]) f_combined);
+   #endif
+   */
    return;
 }
 #endif
-
-/*
-
-#include "mex.h"
-#include "config.h"
-#include "ltfat.h"
-#include "ltfat-mex-helper.h"
-
-// Calling convention:
-// comp_dgtreal_fb(f,g,a,M);
-
-
- void mexFunction( int nlhs, mxArray *plhs[],
-		  int nrhs, const mxArray *prhs[] )
-
-{
-   int L, gl,W, a, M, N, M2;
-   ltfat_complex *out_combined;
-   double *f, *g;
-
-   // Get matrix dimensions.
-   L  = mxGetM(prhs[0]);
-   W  = mxGetN(prhs[0]);
-   gl = mxGetM(prhs[1]);
-
-   a=(int)mxGetScalar(prhs[2]);
-   M=(int)mxGetScalar(prhs[3]);
-   M2=M/2+1;
-
-   N=L/a;
-
-   // Create temporary matrices to convert to correct complex layout.
-
-   out_combined=mxMalloc(M2*N*W*sizeof(ltfat_complex));
-
-   f=mxGetPr(prhs[0]);
-   g=mxGetPr(prhs[1]);
-
-   dgtreal_fb((const double*)f,(const double*)g,L,gl,W,a,M,
-              out_combined);
-
-   plhs[0] = mxCreateDoubleMatrix(M2, N*W, mxCOMPLEX);
-
-   combined2split(M2*N*W, (const ltfat_complex*)out_combined, plhs[0]);
-
-
-   mxFree(out_combined);
-   return;
-
-}
-
-*/
mex/comp_dgtreal_long.c to mex/comp_sepdgtreal.c
--- a/mex/comp_dgtreal_long.c
+++ b/mex/comp_sepdgtreal.c
@@ -15,77 +15,43 @@
 #include "ltfat_types.h"
 
 // Calling convention:
-// comp_dgt_long(f,g,a,M);
+// comp_sepdgtreal(f,g,a,M);
 
 void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
 {
-   int L, W, a, M, N, M2;
+   int L, gl,W, a, M, N, M2;
 
    // Get matrix dimensions.
-   L = mxGetM(prhs[0]);
-   W = mxGetN(prhs[0]);
+   L  = mxGetM(prhs[0]);
+   W  = mxGetN(prhs[0]);
+   gl = mxGetM(prhs[1]);
 
    a=(int)mxGetScalar(prhs[2]);
    M=(int)mxGetScalar(prhs[3]);
+   M2=M/2+1;
+   N=L/a;
+   
+   mwSize ndim = 3;
+   if(W==1)
+   {
+      ndim = 2;
+   }
 
-   N=L/a;
-   M2 = (M/2)+1;
-
-   plhs[0] = ltfatCreateMatrix(M2, N*W,LTFAT_MX_CLASSID,mxCOMPLEX);
+   mwSize dims[] = {M2, N, W};
+   plhs[0] = ltfatCreateNdimArray(ndim, dims,LTFAT_MX_CLASSID,mxCOMPLEX);
    const LTFAT_REAL * f = (const LTFAT_REAL *) mxGetData(prhs[0]);
    const LTFAT_REAL * g = (const LTFAT_REAL *) mxGetData(prhs[1]);
    LTFAT_COMPLEX* out_combined = (LTFAT_COMPLEX*) mxGetData(plhs[0]);
 
-   LTFAT_NAME(dgtreal_long)(f,g, L, W, a, M, out_combined);
-
+   if(gl<L)
+   {
+      LTFAT_NAME(dgtreal_fb)(f,g,L,gl,W,a,M,out_combined);
+   }
+   else
+   {
+      LTFAT_NAME(dgtreal_long)(f,g, L, W, a, M, out_combined);
+   }
    return;
 }
 #endif
 
-/*
-#include "mex.h"
-#include "config.h"
-#include "ltfat.h"
-#include "ltfat-mex-helper.h"
-
-// Calling convention:
-// comp_dgt_long(f,g,a,M);
-
-
-void mexFunction( int nlhs, mxArray *plhs[],
-		  int nrhs, const mxArray *prhs[] )
-
-{
-   int L, W, a, M, N, M2;
-
-   double *f, *g;
-   ltfat_complex *out_combined;
-
-   // Get matrix dimensions.
-   L = mxGetM(prhs[0]);
-   W = mxGetN(prhs[0]);
-
-   a=(int)mxGetScalar(prhs[2]);
-   M=(int)mxGetScalar(prhs[3]);
-
-   N=L/a;
-   M2 = (M/2)+1;
-
-   f=mxGetPr(prhs[0]);
-   g=mxGetPr(prhs[1]);
-
-   out_combined=(ltfat_complex*)mxMalloc(M2*N*W*sizeof(ltfat_complex));
-
-   dgtreal_long((const double*)f,(const double*)g,
-		L, W, a, M, out_combined);
-
-   plhs[0] = mxCreateDoubleMatrix(M2, N*W, mxCOMPLEX);
-
-   combined2split(M2*N*W, (const ltfat_complex*)out_combined, plhs[0]);
-
-   mxFree(out_combined);
-   return;
-
-}
-*/
-
mex/comp_dwilt_long.c to mex/comp_dwilt.c
--- a/mex/comp_dwilt_long.c
+++ b/mex/comp_dwilt.c
@@ -1,7 +1,7 @@
 #ifndef _LTFAT_MEX_FILE
 #define _LTFAT_MEX_FILE
 
-#define ISNARGINEQ 4
+#define ISNARGINEQ 3
 #define TYPEDEPARGS 0, 1
 #define SINGLEARGS
 #define COMPLEXINDEPENDENT
@@ -15,18 +15,19 @@
 #include "ltfat_types.h"
 
 // Calling convention:
-//  comp_dwilt_long(f,g,M,L);
+//  comp_dwilt(f,g,M);
 
 void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
 {
-   int M, N, L, W;
+   int M, N, L, gl, W;
 
    mwSize ndim;
    mwSize dims[3];
 
    // Get matrix dimensions.
    M=(int)mxGetScalar(prhs[2]);
-   L=(int)mxGetScalar(prhs[3]);
+   L=(int)mxGetM(prhs[0]);
+   gl=(int) mxGetM(prhs[1]);
    W = mxGetN(prhs[0]);
 
    N=L/M;
@@ -50,10 +51,13 @@
    const LTFAT_TYPE* g = (const LTFAT_TYPE*) mxGetData(prhs[1]);
    LTFAT_TYPE* cout = (LTFAT_TYPE*) mxGetData(plhs[0]);
 
-   LTFAT_NAME(dwilt_long)(f,g,L, W, M, cout);
-
-
-
-   return;
+   if(gl<L)
+   {
+      LTFAT_NAME(dwilt_fb)(f,g,L,gl, W, M, cout);
+   }
+   else
+   {
+      LTFAT_NAME(dwilt_long)(f,g,L,W,M,cout);
+   }
 }
 #endif
mex/comp_idgt_fac.c to oct/comp_isepdgt.cc
--- a/mex/comp_idgt_fac.c
+++ b/oct/comp_isepdgt.cc
@@ -1,94 +1,85 @@
-#ifndef _LTFAT_MEX_FILE
-#define _LTFAT_MEX_FILE
-
-#define ISNARGINEQ 5
 #define TYPEDEPARGS 0, 1
 #define SINGLEARGS
 #define COMPLEXARGS
-
-#endif // _LTFAT_MEX_FILE - INCLUDED ONCE
-
-#define MEX_FILE __BASE_FILE__
-#include "ltfat_mex_template_helper.h"
-
-#if defined(LTFAT_SINGLE) || defined(LTFAT_DOUBLE)
-#include "ltfat_types.h"
-
-// Calling convention:
-//  comp_idgt_fac(coef,gf,L,a,M);
-
-void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
-{
-   int L, W, a, M, N;
-
-   // Get matrix dimensions.
-   L=(int)mxGetScalar(prhs[2]);
-   a=(int)mxGetScalar(prhs[3]);
-   M=(int)mxGetScalar(prhs[4]);
-   N=L/a;
-   W = mxGetM(prhs[0])*mxGetN(prhs[0])/(M*N);
+#define OCTFILENAME comp_isepdgt // change to filename
+#define OCTFILEHELP "This function calls the C-library\n\
+                    c=comp_idgt_fb(coef,g,L,a,M);\n\
+                    Yeah."
 
 
-   plhs[0] = ltfatCreateMatrix(L, W,LTFAT_MX_CLASSID, mxCOMPLEX);
-   const LTFAT_COMPLEX* c_combined = (const LTFAT_COMPLEX*) mxGetData(prhs[0]);
-   const LTFAT_COMPLEX* gf_combined = (const LTFAT_COMPLEX*) mxGetData(prhs[1]);
-   LTFAT_COMPLEX* f_combined = (LTFAT_COMPLEX*) mxGetData(plhs[0]);
+#include "ltfat_oct_template_helper.h"
+// octave_idx_type is 32 or 64 bit signed integer
 
-   LTFAT_NAME(idgt_fac)(c_combined,gf_combined,L,W,a,M,f_combined);
+static inline void
+fwd_idgt_fb(const Complex *coef, const Complex *gf,
+            const octave_idx_type L,const octave_idx_type gl,
+            const octave_idx_type W,const octave_idx_type a,
+            const octave_idx_type M,Complex *f)
+{
+   idgt_fb_d(reinterpret_cast<const double _Complex*>(coef),
+             reinterpret_cast<const double _Complex*>(gf),
+             L,gl,W,a,M,reinterpret_cast<double _Complex*>(f));
+}
 
+static inline void
+fwd_idgt_fb(const FloatComplex *coef, const FloatComplex *gf,
+            const octave_idx_type L,const octave_idx_type gl, 
+            const octave_idx_type W,const octave_idx_type a,
+            const octave_idx_type M,FloatComplex *f)
+{
+   idgt_fb_s(reinterpret_cast<const float _Complex*>(coef),
+             reinterpret_cast<const float _Complex*>(gf),
+             L,gl,W,a,M,reinterpret_cast<float _Complex*>(f));
+}
 
+static inline void
+fwd_idgt_long(const Complex *coef, const Complex *gf,
+              const octave_idx_type L, const octave_idx_type W,
+              const octave_idx_type a, const octave_idx_type M,
+              Complex *f)
+{
+   idgt_long_d(reinterpret_cast<const double _Complex*>(coef),
+               reinterpret_cast<const double _Complex*>(gf),
+               L,W,a,M,
+               reinterpret_cast<double _Complex*>(f));
+}
 
-   return;
+static inline void 
+fwd_idgt_long(const FloatComplex *coef, const FloatComplex *gf,
+              const octave_idx_type L, const octave_idx_type W, 
+              const octave_idx_type a, const octave_idx_type M,
+              FloatComplex *f)
+{
+   idgt_long_s(reinterpret_cast<const float _Complex*>(coef),
+               reinterpret_cast<const float _Complex*>(gf),
+               L,W,a,M,
+               reinterpret_cast<float _Complex*>(f));
 }
-#endif
 
-/*
-#include "mex.h"
-#include "config.h"
-#include "ltfat.h"
-#include "ltfat-mex-helper.h"
+template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
+octave_value_list octFunction(const octave_value_list& args, int nargout)
+{
+   DEBUGINFO;
 
-// Calling convention:
-//  comp_idgt_fac(coef,gf,L,a,M);
+   MArray<LTFAT_TYPE> coef = ltfatOctArray<LTFAT_TYPE>(args(0));
+   MArray<LTFAT_TYPE> gf = ltfatOctArray<LTFAT_TYPE>(args(1));
+   const octave_idx_type L = args(2).int_value();
+   const octave_idx_type a = args(3).int_value();
+   const octave_idx_type M = args(4).int_value();
+   const octave_idx_type N = L/a;
+   const octave_idx_type gl = gf.rows();
 
+   const octave_idx_type W = coef.nelem()/(M*N);
 
-void mexFunction( int nlhs, mxArray *plhs[],
-		  int nrhs, const mxArray *prhs[] )
+   MArray<LTFAT_COMPLEX> f(dim_vector(L,W)); 
 
-{
-   int L, W, a, M, N;
-   ltfat_complex *f_combined, *gf_combined, *c_combined;
-
-   // Get matrix dimensions.
-   L=(int)mxGetScalar(prhs[2]);
-   a=(int)mxGetScalar(prhs[3]);
-   M=(int)mxGetScalar(prhs[4]);
-   N=L/a;
-   W = mxGetM(prhs[0])*mxGetN(prhs[0])/(M*N);
-
-   // Create temporary matrices to convert to correct complex layout.
-
-   c_combined=mxMalloc(M*N*W*sizeof(ltfat_complex));
-   gf_combined=mxMalloc(L*sizeof(ltfat_complex));
-   f_combined=mxMalloc(L*W*sizeof(ltfat_complex));
-
-   split2combined(M*N*W, prhs[0], c_combined);
-   split2combined(L, prhs[1], gf_combined);
-
-   idgt_fac((const ltfat_complex*)c_combined,(const ltfat_complex*)gf_combined,L,W,a,M,
-	    f_combined);
-
-   mxFree(c_combined);
-   mxFree(gf_combined);
-
-   plhs[0] = mxCreateDoubleMatrix(L, W, mxCOMPLEX);
-
-   combined2split(L*W, (const ltfat_complex*)f_combined, plhs[0]);
-
-   mxFree(f_combined);
-
-   return;
-
+   if(gl<L)
+   {
+      fwd_idgt_fb(coef.data(),gf.data(),L,gl,W,a,M,f.fortran_vec());
+   }
+   else
+   {
+      fwd_idgt_long(coef.data(),gf.data(),L,W,a,M,f.fortran_vec());
+   }
+   return octave_value(f);
 }
-*/
-
mex/comp_idgt_fb.c to src/idwilt.c
--- a/mex/comp_idgt_fb.c
+++ b/src/idwilt.c
@@ -1,116 +1,154 @@
-#ifndef _LTFAT_MEX_FILE
-#define _LTFAT_MEX_FILE
-
-#define ISNARGINEQ 5
-#define TYPEDEPARGS 0, 1
-#define SINGLEARGS
-#define COMPLEXARGS
-
-#endif // _LTFAT_MEX_FILE - INCLUDED ONCE
-
-#define MEX_FILE __BASE_FILE__
-#include "ltfat_mex_template_helper.h"
-
-#if defined(LTFAT_SINGLE) || defined(LTFAT_DOUBLE)
+#include "config.h"
+#include "ltfat.h"
 #include "ltfat_types.h"
 
-// Calling convention:
-//  comp_idgt_fb(coef,g,L,a,M);
+#define CH(name) LTFAT_COMPLEXH_NAME(name)
 
-void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
-{
-   int L, W, a, M, N, gl;
-   // Get matrix dimensions.
-   L=(int)mxGetScalar(prhs[2]);
-   a=(int)mxGetScalar(prhs[3]);
-   M=(int)mxGetScalar(prhs[4]);
-   N=L/a;
+#define PREPROC_REAL \
+  for (int n=0;n<N*W;n+=2) \
+  { \
+     pcoef2[0]=pcoef[0]; \
+\
+     for (int m=1;m<M;m+=2) \
+     { \
+        pcoef2[m] = -I*scalconst*(pcoef[m]); \
+        pcoef2[m+coef2_ld] = scalconst*(pcoef[m+M]); \
+     } \
+ \
+     for (int m=2;m<M;m+=2) \
+     { \
+        pcoef2[m] = scalconst*(pcoef[m]); \
+        pcoef2[m+coef2_ld] = -I*scalconst*(pcoef[m+M]); \
+     } \
+ \
+     pcoef2[M+nyquestadd] = pcoef[M]; \
+     pcoef+=2*M; \
+     pcoef2+=2*coef2_ld; \
+  }
 
-   gl = mxGetM(prhs[1]);
-   W  = mxGetM(prhs[0])*mxGetN(prhs[0])/(M*N);
-
-   plhs[0] = ltfatCreateMatrix(L, W,LTFAT_MX_CLASSID, mxCOMPLEX);
-   const LTFAT_COMPLEX* c_combined = (const LTFAT_COMPLEX*) mxGetData(prhs[0]);
-   const LTFAT_COMPLEX* g_combined = (const LTFAT_COMPLEX*) mxGetData(prhs[1]);
-   LTFAT_COMPLEX* f_combined = (LTFAT_COMPLEX*) mxGetData(plhs[0]);
-
-   //#ifdef LTFAT_COMPLEXTYPE
-   LTFAT_NAME(idgt_fb)(c_combined,g_combined, L,gl,W,a,M, f_combined);
-
-  /* #else
-   NOT CALLING idgt_fb_r:
-   TO DO: Do it better.
-   LTFAT_NAME(idgt_fb_r)((const LTFAT_REAL (*)[2])c_combined,
-                         g_combined,
-                         L,gl,W,a,M,(LTFAT_REAL (*)[2]) f_combined);
-   #endif
-   */
-   return;
-}
-#endif
-
-/*
-#include "mex.h"
-#include "config.h"
-#include "ltfat.h"
-#include "ltfat-mex-helper.h"
+#define PREPROC_COMPLEX \
+  for (int n=0;n<N*W;n+=2) \
+  { \
+     pcoef2[0] = pcoef[0]; \
+ \
+     for (int m=1;m<M;m+=2) \
+     { \
+        pcoef2[m] = -I*scalconst*pcoef[m]; \
+        pcoef2[M2-m] = I*scalconst*pcoef[m]; \
+        pcoef2[M2+m] =  scalconst*pcoef[m+M]; \
+        pcoef2[M4-m] = scalconst*pcoef[m+M]; \
+     } \
+ \
+     for (int m=2;m<M;m+=2) \
+     { \
+        pcoef2[m] = scalconst*pcoef[m]; \
+        pcoef2[M2-m] = scalconst*pcoef[m]; \
+        pcoef2[M2+m] =  -I*scalconst*pcoef[m+M]; \
+        pcoef2[M4-m] = I*scalconst*pcoef[m+M]; \
+     } \
+ \
+     pcoef2[M+nyquestadd] = pcoef[M]; \
+     pcoef+=M2; \
+     pcoef2+=M4; \
+  }
 
 
-// Calling convention:
-//  comp_idgt_fb(coef,g,L,a,M);
+LTFAT_EXTERN void
+LTFAT_NAME_COMPLEX(idwilt_long)(const LTFAT_COMPLEX *c, const LTFAT_COMPLEX *g,
+			   const int L, const int W, const int M,
+			   LTFAT_COMPLEX *f)
+{
+   const int N=L/M;
+   const int M2=2*M;
+   const int M4=4*M;
+   const LTFAT_REAL scalconst = 1.0/sqrt(2.0);
 
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_calloc(2*M*N*W,sizeof(LTFAT_COMPLEX));
 
-void mexFunction( int nlhs, mxArray *plhs[],
-		  int nrhs, const mxArray *prhs[] )
+  const int nyquestadd = (M%2)*M2;
 
-{
-   int L, W, a, M, N, gl;
-   ltfat_complex *f_combined, *g_combined, *c_combined;
+  LTFAT_COMPLEX *pcoef  = c;
+  LTFAT_COMPLEX *pcoef2 = coef2;
 
-   // Get matrix dimensions.
-   L=(int)mxGetScalar(prhs[2]);
-   a=(int)mxGetScalar(prhs[3]);
-   M=(int)mxGetScalar(prhs[4]);
-   N=L/a;
+  PREPROC_COMPLEX
 
-   gl = mxGetM(prhs[1]);
-   W  = mxGetM(prhs[0])*mxGetN(prhs[0])/(M*N);
+  LTFAT_NAME(idgt_long)(coef2, g, L, W, M, 2*M, f);
 
-   // Create temporary matrices to convert to correct complex layout.
-
-   c_combined=mxMalloc(M*N*W*sizeof(ltfat_complex));
-   f_combined=mxMalloc(L*W*sizeof(ltfat_complex));
-
-   split2combined(M*N*W, prhs[0], c_combined);
-
-   if (mxIsComplex(prhs[1]))
-   {
-
-      g_combined=mxMalloc(gl*sizeof(ltfat_complex));
-      split2combined(   gl, prhs[1], g_combined);
-
-      idgt_fb((const ltfat_complex*)c_combined,(const ltfat_complex*)g_combined,L,gl,W,a,M,
-	      f_combined);
-
-      mxFree(g_combined);
-   }
-   else
-   {
-      idgt_fb_r((const ltfat_complex*)c_combined,
-		(const double*)mxGetPr(prhs[1]),L,gl,W,a,M,
-		f_combined);
-   }
-   mxFree(c_combined);
-
-
-   plhs[0] = mxCreateDoubleMatrix(L, W, mxCOMPLEX);
-
-   combined2split(L*W, (const ltfat_complex*)f_combined, plhs[0]);
-
-   mxFree(f_combined);
-
-   return;
+  ltfat_free(coef2);
 
 }
-*/
 
+LTFAT_EXTERN void
+LTFAT_NAME_REAL(idwilt_long)(const LTFAT_REAL *c, const LTFAT_REAL *g,
+			   const int L, const int W, const int M,
+			   LTFAT_REAL *f)
+{
+   const int N=L/M;
+   const int coef2_ld = M+1;
+   const LTFAT_REAL scalconst = 1.0/sqrt(2.0);
+   const int nyquestadd = (M%2)*coef2_ld;
+
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_calloc((M+1)*N*W,sizeof(LTFAT_COMPLEX));
+
+  LTFAT_REAL *pcoef  = c;
+  LTFAT_COMPLEX *pcoef2 = coef2;
+
+  PREPROC_REAL
+
+  LTFAT_NAME(idgtreal_long)(coef2, g, L, W, M, 2*M, f);
+
+  ltfat_free(coef2);
+
+}
+
+LTFAT_EXTERN void
+LTFAT_NAME_COMPLEX(idwilt_fb)(const LTFAT_COMPLEX *c, const LTFAT_COMPLEX *g,
+			    const Lint L, const Lint gl, const Lint W, const Lint M,
+			   LTFAT_COMPLEX *f)
+{
+   const int N=L/M;
+   const int M2=2*M;
+   const int M4=4*M;
+   const LTFAT_REAL scalconst = 1.0/sqrt(2.0);
+
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_calloc(2*M*N*W,sizeof(LTFAT_COMPLEX));
+
+  const int nyquestadd = (M%2)*M2;
+
+  LTFAT_COMPLEX *pcoef  = c;
+  LTFAT_COMPLEX *pcoef2 = coef2;
+
+  PREPROC_COMPLEX
+
+  LTFAT_NAME(idgt_fb)(coef2, g, L, gl, W, M, 2*M, f);
+
+  ltfat_free(coef2);
+
+}
+
+LTFAT_EXTERN void
+LTFAT_NAME_REAL(idwilt_fb)(const LTFAT_REAL *c, const LTFAT_REAL *g,
+			   const Lint L, const Lint gl, const Lint W, const Lint M,
+			   LTFAT_REAL *f)
+{
+   const Lint N = L/M;
+   const Lint coef2_ld = M + 1;
+   const Lint nyquestadd = (M%2)*coef2_ld;
+   const LTFAT_REAL scalconst = (LTFAT_REAL) 1.0/sqrt(2.0);
+
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_calloc((M+1)*N*W,sizeof(LTFAT_COMPLEX));
+
+   LTFAT_REAL* pcoef  = c;
+   LTFAT_COMPLEX* pcoef2 = coef2;
+
+   PREPROC_REAL
+
+   LTFAT_NAME(idgtreal_fb)(coef2, g, L, gl, W, M, 2*M, f);
+
+   ltfat_free(coef2);
+}
+
+#undef CH
+#undef PREPROC_REAL
+#undef PREPROC_COMPLEX
+
mex/comp_idgtreal_fac.c to mex/comp_sepdgt.c
--- a/mex/comp_idgtreal_fac.c
+++ b/mex/comp_sepdgt.c
@@ -1,10 +1,10 @@
 #ifndef _LTFAT_MEX_FILE
 #define _LTFAT_MEX_FILE
 
-#define ISNARGINEQ 5
+#define ISNARGINEQ 4
 #define TYPEDEPARGS 0, 1
 #define SINGLEARGS
-#define COMPLEXARGS
+#define COMPLEXINDEPENDENT
 
 #endif // _LTFAT_MEX_FILE - INCLUDED ONCE
 
@@ -15,82 +15,39 @@
 #include "ltfat_types.h"
 
 // Calling convention:
-// comp_idgtreal_fac(coef,gf,L,a,M);
+//  comp_dgt_fb(f,g,a,M);
 
 void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
 {
-   int L, W, a, M, N;
+   int L  = mxGetM(prhs[0]);
+   int W  = mxGetN(prhs[0]);
+   int gl = mxGetNumberOfElements(prhs[1]);
 
-   // Get matrix dimensions.
-   L=(int)mxGetScalar(prhs[2]);
-   a=(int)mxGetScalar(prhs[3]);
-   M=(int)mxGetScalar(prhs[4]);
-   N=L/a;
-   W = mxGetN(prhs[0])/N;
+   int a=(int)mxGetScalar(prhs[2]);
+   int M=(int)mxGetScalar(prhs[3]);
 
-   plhs[0] = ltfatCreateMatrix(L, W,LTFAT_MX_CLASSID,mxREAL);
-   const LTFAT_COMPLEX* c_combined = (const LTFAT_COMPLEX*) mxGetData(prhs[0]);
-   const LTFAT_COMPLEX* gf_combined = (const LTFAT_COMPLEX*) mxGetData(prhs[1]);
-   LTFAT_REAL* f_r = (LTFAT_REAL*) mxGetData(plhs[0]);
+   int N=L/a;
 
-   LTFAT_NAME(idgtreal_fac)(c_combined,gf_combined,L,W,a,M,f_r);
+   mwSize ndim=3;
+   if (W==1)
+   {
+      ndim=2;
+   }
 
+   mwSize dims[3] = {M, N, W};
+   plhs[0] = ltfatCreateNdimArray(ndim,dims,LTFAT_MX_CLASSID,mxCOMPLEX);
+   const LTFAT_TYPE* f_combined = (const LTFAT_TYPE*) mxGetData(prhs[0]);
+   const LTFAT_TYPE* g_combined = (const LTFAT_TYPE*) mxGetData(prhs[1]);
+   LTFAT_COMPLEX* out_combined = (LTFAT_COMPLEX*) mxGetData(plhs[0]);
 
-   return;
+   if(gl<L)
+   {
+      LTFAT_NAME(dgt_fb)(f_combined, g_combined,L,gl,W,a,M,out_combined);
+   }
+   else
+   {
+      LTFAT_NAME(dgt_long)(f_combined, g_combined,L,W,a,M,out_combined);
+   }
 }
 #endif
 
-/*
-#include "mex.h"
-#include "config.h"
-#include "ltfat.h"
-#include "ltfat-mex-helper.h"
-
-
-// Calling convention:
-// comp_idgtreal_fac(coef,gf,L,a,M);
-
-void mexFunction( int nlhs, mxArray *plhs[],
-		  int nrhs, const mxArray *prhs[] )
-
-{
-   int L, W, a, M, N, M2;
-   ltfat_complex *gf_combined, *c_combined;
-   double *f_r,*gf_r,*gf_i,*c_r,*c_i;
-
-   int ii;
-
-   // This is a floor operation.
-
-   // Get matrix dimensions.
-   L=(int)mxGetScalar(prhs[2]);
-   a=(int)mxGetScalar(prhs[3]);
-   M=(int)mxGetScalar(prhs[4]);
-   N=L/a;
-   W = mxGetN(prhs[0])/N;
-
-   M2= M/2+1;
-
-   // Create temporary matrices to convert to correct complex layout.
-   c_combined=mxMalloc(M2*N*W*sizeof(ltfat_complex));
-   gf_combined=mxMalloc(L*sizeof(ltfat_complex));
-
-   split2combined(M2*N*W, prhs[0], c_combined);
-   split2combined(L, prhs[1], gf_combined);
-
-   plhs[0] = mxCreateDoubleMatrix(L, W, mxREAL);
-
-   f_r=mxGetPr(plhs[0]);
-
-   LTFAT_NAME_DOUBLE(idgtreal_fac)((const ltfat_complex*)c_combined,
-		(const ltfat_complex*)gf_combined,L,W,a,M,
-		(double*)f_r);
-
-   mxFree(c_combined);
-   mxFree(gf_combined);
-
-   return;
-
-}
-*/
-
mex/comp_idgtreal_fb.c to mex/comp_isepdgtreal.c
--- a/mex/comp_idgtreal_fb.c
+++ b/mex/comp_isepdgtreal.c
@@ -17,79 +17,37 @@
 #include "ltfat_types.h"
 
 // Calling convention:
-//  comp_idgtreal_fb(coef,g,L,a,M);
+//  comp_isepdgtreal(coef,g,L,a,M);
 
 void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
 {
-   int L, W, a, M, N, gl;
+   int L, W, a, M, N, gl, M2;
 
    // Get matrix dimensions.
    L=(int)mxGetScalar(prhs[2]);
    a=(int)mxGetScalar(prhs[3]);
    M=(int)mxGetScalar(prhs[4]);
    N=L/a;
-   W = mxGetN(prhs[0])/N;
+   M2 = M/2+1;
+   W = mxGetNumberOfElements(prhs[0])/(N*M2);
    gl = mxGetM(prhs[1]);
 
 
    plhs[0] = ltfatCreateMatrix(L, W,LTFAT_MX_CLASSID,mxREAL);
    const LTFAT_COMPLEX* c_combined = (const LTFAT_COMPLEX*) mxGetData(prhs[0]);
-   const LTFAT_REAL * gf = (const LTFAT_REAL *) mxGetData(prhs[1]);
+   const LTFAT_REAL * g = (const LTFAT_REAL *) mxGetData(prhs[1]);
    LTFAT_REAL* f_r = (LTFAT_REAL*) mxGetData(plhs[0]);
 
-   LTFAT_NAME(idgtreal_fb)(c_combined,gf,L,gl,W,a,M,f_r);
+   if(gl<L)
+   {    
+      LTFAT_NAME(idgtreal_fb)(c_combined,g,L,gl,W,a,M,f_r);
+   }
+   else
+   {
+      LTFAT_NAME(idgtreal_long)(c_combined,g,L,W,a,M,f_r);
+   }
 
 
    return;
 }
 #endif
-
-/*
-#include "mex.h"
-#include "config.h"
-#include "ltfat.h"
-#include "ltfat-mex-helper.h"
-
-
-// Calling convention:
-//  comp_idgtreal_fb(coef,g,L,a,M);
-
-
-void mexFunction( int nlhs, mxArray *plhs[],
-		  int nrhs, const mxArray *prhs[] )
-
-{
-   int L, W, a, M, N, M2, gl;
-   ltfat_complex *c_combined;
-
-   int ii;
-
-   // Get matrix dimensions.
-   L=(int)mxGetScalar(prhs[2]);
-   a=(int)mxGetScalar(prhs[3]);
-   M=(int)mxGetScalar(prhs[4]);
-   N=L/a;
-   W = mxGetN(prhs[0])/N;
-   gl = mxGetM(prhs[1]);
-
-
-   M2= M/2+1;
-
-   // Create temporary matrices to convert to correct complex layout.
-   c_combined=mxMalloc(M2*N*W*sizeof(ltfat_complex));
-
-   split2combined(M2*N*W, prhs[0], c_combined);
-
-   plhs[0] = mxCreateDoubleMatrix(L, W, mxREAL);
-
-   LTFAT_NAME_DOUBLE(idgtreal_fb)((const ltfat_complex*)c_combined,
-		(const double*)mxGetPr(prhs[1]),L,gl,W,a,M,
-		(double*)mxGetPr(plhs[0]));
-
-   mxFree(c_combined);
-
-   return;
-
-}
-*/
-
oct/comp_dgt_fb.cc to oct/comp_dwilt.cc
--- a/oct/comp_dgt_fb.cc
+++ b/oct/comp_dwilt.cc
@@ -1,145 +1,124 @@
 #define TYPEDEPARGS 0, 1
 #define SINGLEARGS
 #define COMPLEXINDEPENDENT
-#define OCTFILENAME comp_dgt_fb // change to filename
-#define OCTFILEHELP "This function calls the C-library\n c=comp_dgt_fb(f,g,a,M,boundary);\n Yeah."
+#define OCTFILENAME comp_dwilt // change to filename
+#define OCTFILEHELP "This function calls the C-library\n\
+                     coef=comp_dwilt_fb(f,g,M,L);\n\
+                     Yeah."
+
 
 #include "ltfat_oct_template_helper.h"
-// octave_idx_type 32 or 64 bit signed integer
+// octave_idx_type is 32 or 64 bit signed integer
 /*
-  dgt_fb forwarders
+  dgtreal_ola forwarders
 */
 
-static inline void fwd_dgt_fb(const Complex *f, const Complex *g,
-                              const octave_idx_type L, const octave_idx_type gl,
-		                      const octave_idx_type W, const octave_idx_type a,
-							  const octave_idx_type M, Complex *cout)
+static inline void
+fwd_dwilt_fb(const Complex *f, const Complex *g,
+             const octave_idx_type L, const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             Complex *cout)
 {
-   dgt_fb_cd(reinterpret_cast<const double _Complex*>(f),
-             reinterpret_cast<const double _Complex*>(g),
-             L,gl,W,a,M,
-			 reinterpret_cast<double _Complex*>(cout));
+   dwilt_fb_cd(reinterpret_cast<const double _Complex*>(f),
+               reinterpret_cast<const double _Complex*>(g),
+               L, gl, W,M, reinterpret_cast<double _Complex*>(cout));
 }
 
-static inline void fwd_dgt_fb(const FloatComplex *f, const FloatComplex *g,
-                              const octave_idx_type L, const octave_idx_type gl,
-		                      const octave_idx_type W, const octave_idx_type a,
-							  const octave_idx_type M, FloatComplex *cout)
+static inline void 
+fwd_dwilt_fb(const FloatComplex *f, const FloatComplex *g,
+             const octave_idx_type L,  const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             FloatComplex *cout)
 {
-   dgt_fb_cs(reinterpret_cast<const float _Complex*>(f),
-             reinterpret_cast<const float _Complex*>(g),
-             L,gl,W,a,M,
-			 reinterpret_cast<float _Complex*>(cout));
+   dwilt_fb_cs(reinterpret_cast<const float _Complex*>(f),
+               reinterpret_cast<const float _Complex*>(g),
+               L, gl,W,M,
+               reinterpret_cast<float _Complex*>(cout));
 }
 
-static inline void fwd_dgt_fb(const double *f, const double *g,
-                              const octave_idx_type L, const octave_idx_type gl,
-		                      const octave_idx_type W, const octave_idx_type a,
-							  const octave_idx_type M, Complex *cout)
+static inline void
+fwd_dwilt_fb(const double *f, const double *g,
+             const octave_idx_type L,  const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             double *cout)
 {
-   dgt_fb_d(reinterpret_cast<const double*>(f),
-            reinterpret_cast<const double*>(g),
-            L,gl,W,a,M,
-			reinterpret_cast<double _Complex*>(cout));
+   dwilt_fb_d(f,g,L,gl,W,M,cout);
 }
 
-static inline void fwd_dgt_fb(const float *f, const float *g,
-                              const octave_idx_type L, const octave_idx_type gl,
-		                      const octave_idx_type W, const octave_idx_type a,
-							  const octave_idx_type M, FloatComplex *cout)
+static inline void
+fwd_dwilt_fb(const float *f, const float *g,
+             const octave_idx_type L,  const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             float *cout)
 {
-   dgt_fb_s(reinterpret_cast<const float*>(f),
-            reinterpret_cast<const float*>(g),
-            L,gl,W,a,M,
-			reinterpret_cast<float _Complex*>(cout));
+   dwilt_fb_s(f,g,L,gl,W,M,cout);
+}
+
+static inline void
+fwd_dwilt_long(const Complex *f, const Complex *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, Complex *cout)
+{
+   dwilt_long_cd(reinterpret_cast<const double _Complex*>(f),
+                 reinterpret_cast<const double _Complex*>(g),
+                 L,W,M, reinterpret_cast<double _Complex*>(cout));
+}
+
+static inline void 
+fwd_dwilt_long(const FloatComplex *f, const FloatComplex *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, FloatComplex *cout)
+{
+   dwilt_long_cs(reinterpret_cast<const float _Complex*>(f),
+                 reinterpret_cast<const float _Complex*>(g),
+                 L,W,M, reinterpret_cast<float _Complex*>(cout));
+}
+
+static inline void
+fwd_dwilt_long(const double *f, const double *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, double *cout)
+{
+   dwilt_long_d(f,g,L,W,M,cout);
+}
+
+static inline void
+fwd_dwilt_long(const float *f, const float *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, float *cout)
+{
+   dwilt_long_s(f,g,L,W,M,cout);
 }
 
 template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
 octave_value_list octFunction(const octave_value_list& args, int nargout)
 {
-     DEBUGINFO;
-	 const octave_idx_type a = args(2).int_value();
-     const octave_idx_type M = args(3).int_value();
-	 
-	 MArray<LTFAT_TYPE> f = ltfatOctArray<LTFAT_TYPE>(args(0));
-	 MArray<LTFAT_TYPE> g = ltfatOctArray<LTFAT_TYPE>(args(1));
-	 const octave_idx_type L  = f.rows();
-     const octave_idx_type W  = f.columns();
-     const octave_idx_type gl = g.rows();
-     const octave_idx_type N = L/a;
-	 
-	 dim_vector dims_out(M,N,W);  
-     dims_out.chop_trailing_singletons();
+   DEBUGINFO;
+   const octave_idx_type M = args(2).int_value();
 
-     MArray<LTFAT_COMPLEX> cout(dims_out); 
-     cout.fill(0);
-	 
-	 fwd_dgt_fb(f.data(),g.data(),L,gl,W,a,M,cout.fortran_vec());
-	 
-     return octave_value(cout);
-}
-/* OLD CODE
-#include <octave/oct.h>
-#include "config.h"
-#include "ltfat.h"
-
-DEFUN_DLD (comp_dgt_fb, args, ,
-  "This function calls the C-library\n\
-  c=comp_dgt_fb(f,g,a,M,boundary);\n\
-  Yeah.")
-{
-
-  const bool use_complex  = args(0).is_complex_type() || args(1).is_complex_type();
-
-  const int a = args(2).int_value();
-  const int M = args(3).int_value();
+   MArray<LTFAT_TYPE> f = ltfatOctArray<LTFAT_TYPE>(args(0));
+   MArray<LTFAT_TYPE> g = ltfatOctArray<LTFAT_TYPE>(args(1));
+   const octave_idx_type gl = g.nelem();
+   const octave_idx_type W = f.columns();
+   const octave_idx_type L = f.rows();
+   const octave_idx_type N = L/M;
 
 
-  if (use_complex)
-  {  
-     const ComplexMatrix f = args(0).complex_matrix_value();
-     const ComplexMatrix g = args(1).complex_matrix_value();
-     
-     const int L  = f.rows();
-     const int W  = f.columns();
-     const int gl = g.rows();
-     const int N = L/a;
+   dim_vector dims_out(2*M,N/2,W);  
+   dims_out.chop_trailing_singletons();
 
-     dim_vector dims_out(M,N,W);  
-     dims_out.chop_trailing_singletons();
-     
-     ComplexNDArray cout(dims_out);  
-     
-     dgt_fb((ltfat_complex*)f.data(),(ltfat_complex*)g.data(),L,gl,W,a,M,
-     (ltfat_complex*)cout.data());  
+   MArray<LTFAT_TYPE> cout(dims_out); 
+   cout.fill(0);
 
-     return octave_value (cout);
-  
-  }
-  else
-  {
+   if(gl<L)
+   {     
+      fwd_dwilt_fb(f.data(),g.data(),L, gl, W, M,cout.fortran_vec());
+   }
+   else
+   {
+      fwd_dwilt_long(f.data(),g.data(),L, W, M,cout.fortran_vec());
+   }
 
-     const Matrix f = args(0).matrix_value();
-     const Matrix g = args(1).matrix_value();
-     
-     const int L  = f.rows();
-     const int W  = f.columns();
-     const int gl = g.rows();
-     const int N = L/a;
-     
-     dim_vector dims_out(M,N,W);  
-     dims_out.chop_trailing_singletons();
+   return octave_value(cout);
+}
 
-     ComplexNDArray cout(dims_out);  
-     
-     dgt_fb_r((double*)f.data(),(double*)g.data(),L,gl,W,a,M,
-     (ltfat_complex*)cout.data());  
-
-     return octave_value (cout);
-
-
-  }
-
-
-}
-*/
oct/comp_dgt_long.cc to oct/comp_idwilt.cc
--- a/oct/comp_dgt_long.cc
+++ b/oct/comp_idwilt.cc
@@ -1,138 +1,120 @@
 #define TYPEDEPARGS 0, 1
 #define SINGLEARGS
-#define COMPLEXARGS
-#define OCTFILENAME comp_dgt_long // change to filename
-#define OCTFILEHELP "This function calls the C-library\n c=comp_dgt_long(f,g,a,M);\n Yeah."
+#define COMPLEXINDEPENDENT
+#define OCTFILENAME comp_idwilt // change to filename
+#define OCTFILEHELP "This function calls the C-library\n\
+                     coef=comp_dwilt_fb(f,g,M,L);\n\
+                     Yeah."
 
 
 #include "ltfat_oct_template_helper.h"
 // octave_idx_type is 32 or 64 bit signed integer
 /*
-  dgt_long forwarders
+  dgtreal_ola forwarders
 */
 
-static inline void fwd_dgt_long(const Complex *f, const Complex *g,
-                              const octave_idx_type L, const octave_idx_type W,
-							  const octave_idx_type a, const octave_idx_type M,
-							  Complex *cout)
+static inline void
+fwd_idwilt_fb(const Complex *c, const Complex *g,
+             const octave_idx_type L, const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             Complex *f)
 {
-   dgt_long_d(reinterpret_cast<const double _Complex*>(f),
-              reinterpret_cast<const double _Complex*>(g),
-              L,W,a,M,
-			  reinterpret_cast<double _Complex*>(cout));
+   idwilt_fb_cd(reinterpret_cast<const double _Complex*>(c),
+               reinterpret_cast<const double _Complex*>(g),
+               L, gl, W,M, reinterpret_cast<double _Complex*>(f));
 }
 
-static inline void fwd_dgt_long(const FloatComplex *f, const FloatComplex *g,
-                              const octave_idx_type L, const octave_idx_type W, 
-							  const octave_idx_type a, const octave_idx_type M,
-							  FloatComplex *cout)
+static inline void 
+fwd_idwilt_fb(const FloatComplex *c, const FloatComplex *g,
+             const octave_idx_type L,  const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             FloatComplex *f)
 {
-   dgt_long_s(reinterpret_cast<const float _Complex*>(f),
-              reinterpret_cast<const float _Complex*>(g),
-              L,W,a,M,
-			  reinterpret_cast<float _Complex*>(cout));
+   idwilt_fb_cs(reinterpret_cast<const float _Complex*>(c),
+               reinterpret_cast<const float _Complex*>(g),
+               L, gl,W,M,
+               reinterpret_cast<float _Complex*>(f));
+}
+
+static inline void
+fwd_idwilt_fb(const double *c, const double *g,
+             const octave_idx_type L,  const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             double *f)
+{
+   idwilt_fb_d(c,g,L,gl,W,M,f);
+}
+
+static inline void
+fwd_idwilt_fb(const float *c, const float *g,
+             const octave_idx_type L,  const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             float *f)
+{
+   idwilt_fb_s(c,g,L,gl,W,M,f);
+}
+
+static inline void
+fwd_idwilt_long(const Complex *c, const Complex *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, Complex *f)
+{
+   idwilt_long_cd(reinterpret_cast<const double _Complex*>(c),
+                 reinterpret_cast<const double _Complex*>(g),
+                 L,W,M, reinterpret_cast<double _Complex*>(f));
+}
+
+static inline void 
+fwd_idwilt_long(const FloatComplex *c, const FloatComplex *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, FloatComplex *f)
+{
+   idwilt_long_cs(reinterpret_cast<const float _Complex*>(c),
+                 reinterpret_cast<const float _Complex*>(g),
+                 L,W,M, reinterpret_cast<float _Complex*>(f));
+}
+
+static inline void
+fwd_idwilt_long(const double *c, const double *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, double *f)
+{
+   idwilt_long_d(c,g,L,W,M,f);
+}
+
+static inline void
+fwd_idwilt_long(const float *c, const float *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, float *f)
+{
+   idwilt_long_s(c,g,L,W,M,f);
 }
 
 template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
 octave_value_list octFunction(const octave_value_list& args, int nargout)
 {
-     DEBUGINFO;
-	 const octave_idx_type a = args(2).int_value();
-     const octave_idx_type M = args(3).int_value();
-	 
-	 MArray<LTFAT_TYPE> f = ltfatOctArray<LTFAT_TYPE>(args(0));
-	 MArray<LTFAT_TYPE> g = ltfatOctArray<LTFAT_TYPE>(args(1));
-	 const octave_idx_type L  = f.rows();
-     const octave_idx_type W  = f.columns();
-     const octave_idx_type N = L/a;
-	 
-	 dim_vector dims_out(M,N,W);  
-     dims_out.chop_trailing_singletons();
+   DEBUGINFO;
 
-     MArray<LTFAT_TYPE> cout(dims_out); 
-     cout.fill(0);
-	 
-	 fwd_dgt_long(f.data(),g.data(),L,W,a,M,cout.fortran_vec());
-	 
-     return octave_value(cout);
+   MArray<LTFAT_TYPE> c = ltfatOctArray<LTFAT_TYPE>(args(0));
+   MArray<LTFAT_TYPE> g = ltfatOctArray<LTFAT_TYPE>(args(1));
+   
+   const octave_idx_type M = c.rows()/2;
+   const octave_idx_type gl = g.nelem();
+   const octave_idx_type N = 2*c.columns();
+   const octave_idx_type L = N*M;
+   const octave_idx_type W = c.nelem()/L;
+
+   MArray<LTFAT_TYPE> f(dim_vector(L,W)); 
+
+   if(gl<L)
+   {     
+      fwd_idwilt_fb(c.data(),g.data(),L, gl, W, M,f.fortran_vec());
+   }
+   else
+   {
+      fwd_idwilt_long(c.data(),g.data(),L, W, M,f.fortran_vec());
+   }
+
+   return octave_value(f);
 }
 
-/* OLD CODE
-#include <octave/oct.h>
-#include "config.h"
-#include "ltfat.h"
-
-
-DEFUN_DLD (comp_dgt_long, args, ,
-  "This function calls the C-library\n\
-  c=comp_dgt_long(f,g,a,M);\n\
-  Yeah.")
-{
-
-  //const bool f_is_complex  = args(0).is_complex_type();
-
-  const int a = args(2).int_value();
-  const int M = args(3).int_value();
-
-  const ComplexMatrix f = args(0).complex_matrix_value();
-  const ComplexMatrix g = args(1).complex_matrix_value();
-  
-  const int L = f.rows();
-  const int W = f.columns();
-  const int N = L/a;
-
-  dim_vector dims_out(M,N,W);  
-  dims_out.chop_trailing_singletons();
-
-  ComplexNDArray cout(dims_out);  
-  
-  dgt_long((ltfat_complex*)f.data(),(ltfat_complex*)g.data(),
-	   L, W, a, M,(ltfat_complex*)cout.data());
-  
-  return octave_value (cout);
-
-
-
-  if (f_is_complex)
-  {
-
-
-  }
-  else
-  {
-
-    const Matrix f = args(0).matrix_value();
-    ltfat_complex *gf;
-
-    const int L = f.rows();
-    const int W = f.columns();
-    const int N = L/a;
-
-    gf = (ltfat_complex*)ltfat_malloc(L*sizeof(ltfat_complex));
-    
-    if (args(1).is_complex_type())
-    {
-       const ComplexMatrix g = args(1).complex_matrix_value();
-       wfac((ltfat_complex*)g.data(),L,1,a,M,gf);
-    }
-    else
-    {
-       const Matrix g = args(1).matrix_value();
-       wfac_r((double*)g.data(),L,1,a,M,gf);
-    }
-    
-    ComplexMatrix cout(M,N*W);  
-    
-    dgt_fac_r((double*)f.data(),gf,L, W, a, M,
-	      (ltfat_complex*)cout.data());
-    
-    ltfat_free(gf);
-
-    return octave_value (cout);
-
-  }
-
-  
-
-}
-*/
oct/comp_dgtreal_fb.cc to oct/comp_idwiltiii.cc
--- a/oct/comp_dgtreal_fb.cc
+++ b/oct/comp_idwiltiii.cc
@@ -1,87 +1,121 @@
 #define TYPEDEPARGS 0, 1
 #define SINGLEARGS
-#define REALARGS
-#define OCTFILENAME comp_dgtreal_fb // change to filename
-#define OCTFILEHELP "This function calls the C-library\n c=comp_dgtreal_fb(f,g,a,M,boundary);\n Yeah."
+#define COMPLEXINDEPENDENT
+#define OCTFILENAME comp_idwiltiii // change to filename
+#define OCTFILEHELP "This function calls the C-library\n\
+                     coef=comp_dwilt_fb(f,g,M,L);\n\
+                     Yeah."
 
 
 #include "ltfat_oct_template_helper.h"
 // octave_idx_type is 32 or 64 bit signed integer
 /*
-  dgtreal_fb forwarders
+  dgtreal_ola forwarders
 */
 
-static inline void fwd_dgtreal_fb(const double *f, const double *g,
-                                  const octave_idx_type L, const octave_idx_type gl,
-							      const octave_idx_type W, const octave_idx_type a,
-							      const octave_idx_type M, Complex *cout)
+static inline void
+fwd_idwiltiii_fb(const Complex *c, const Complex *g,
+             const octave_idx_type L, const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             Complex *f)
 {
-   dgtreal_fb_d(f,g,L,gl,W,a,M,reinterpret_cast<double _Complex*>(cout));
+   idwiltiii_fb_cd(reinterpret_cast<const double _Complex*>(c),
+               reinterpret_cast<const double _Complex*>(g),
+               L, gl, W,M, reinterpret_cast<double _Complex*>(f));
 }
 
-static inline void fwd_dgtreal_fb(const float *f, const float *g,
-                                const octave_idx_type L, const octave_idx_type gl,
-							    const octave_idx_type W, const octave_idx_type a,
-							    const octave_idx_type M, FloatComplex *cout)
+static inline void 
+fwd_idwiltiii_fb(const FloatComplex *c, const FloatComplex *g,
+             const octave_idx_type L,  const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             FloatComplex *f)
 {
-   dgtreal_fb_s(f,g,L,gl,W,a,M,reinterpret_cast<float _Complex*>(cout));
+   idwiltiii_fb_cs(reinterpret_cast<const float _Complex*>(c),
+               reinterpret_cast<const float _Complex*>(g),
+               L, gl,W,M,
+               reinterpret_cast<float _Complex*>(f));
+}
+
+static inline void
+fwd_idwiltiii_fb(const double *c, const double *g,
+             const octave_idx_type L,  const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             double *f)
+{
+   idwiltiii_fb_d(c,g,L,gl,W,M,f);
+}
+
+static inline void
+fwd_idwiltiii_fb(const float *c, const float *g,
+             const octave_idx_type L,  const octave_idx_type gl,
+             const octave_idx_type W, const octave_idx_type M, 
+             float *f)
+{
+   idwiltiii_fb_s(c,g,L,gl,W,M,f);
+}
+
+static inline void
+fwd_idwiltiii_long(const Complex *c, const Complex *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, Complex *f)
+{
+   idwiltiii_long_cd(reinterpret_cast<const double _Complex*>(c),
+                 reinterpret_cast<const double _Complex*>(g),
+                 L,W,M, reinterpret_cast<double _Complex*>(f));
+}
+
+static inline void 
+fwd_idwiltiii_long(const FloatComplex *c, const FloatComplex *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, FloatComplex *f)
+{
+   idwiltiii_long_cs(reinterpret_cast<const float _Complex*>(c),
+                 reinterpret_cast<const float _Complex*>(g),
+                 L,W,M, reinterpret_cast<float _Complex*>(f));
+}
+
+static inline void
+fwd_idwiltiii_long(const double *c, const double *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, double *f)
+{
+   idwiltiii_long_d(c,g,L,W,M,f);
+}
+
+static inline void
+fwd_idwiltiii_long(const float *c, const float *g,
+               const octave_idx_type L, const octave_idx_type W, 
+               const octave_idx_type M, float *f)
+{
+   idwiltiii_long_s(c,g,L,W,M,f);
 }
 
 template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
 octave_value_list octFunction(const octave_value_list& args, int nargout)
 {
-     DEBUGINFO;
-	 const octave_idx_type a = args(2).int_value();
-     const octave_idx_type M = args(3).int_value();
-	 const octave_idx_type M2 = M/2+1;
-	 
-	 MArray<LTFAT_TYPE> f = ltfatOctArray<LTFAT_TYPE>(args(0));
-	 MArray<LTFAT_TYPE> g = ltfatOctArray<LTFAT_TYPE>(args(1));
-	 const octave_idx_type L  = f.rows();
-     const octave_idx_type W  = f.columns();
-	 const octave_idx_type gl = g.rows();
-     const octave_idx_type N = L/a;
-	 
-	 dim_vector dims_out(M2,N*W);  
-     dims_out.chop_trailing_singletons();
+   DEBUGINFO;
 
-     MArray<LTFAT_COMPLEX> cout(dims_out); 
-     cout.fill(0);
-	 
-	 fwd_dgtreal_fb(f.data(),g.data(),L,gl,W,a,M,cout.fortran_vec());
-	 
-     return octave_value(cout);
+   MArray<LTFAT_TYPE> c = ltfatOctArray<LTFAT_TYPE>(args(0));
+   MArray<LTFAT_TYPE> g = ltfatOctArray<LTFAT_TYPE>(args(1));
+   
+   const octave_idx_type M = c.rows();
+   const octave_idx_type gl = g.nelem();
+   const octave_idx_type N = c.columns();
+   const octave_idx_type L = N*M;
+   const octave_idx_type W = c.nelem()/L;
+
+
+   MArray<LTFAT_TYPE> f(dim_vector(L,W)); 
+
+   if(gl<L)
+   {     
+      fwd_idwiltiii_fb(c.data(),g.data(),L, gl, W, M,f.fortran_vec());
+   }
+   else
+   {
+      fwd_idwiltiii_long(c.data(),g.data(),L, W, M,f.fortran_vec());
+   }
+
+   return octave_value(f);
 }
 
-
-/* OLD CODE
-#include <octave/oct.h>
-#include "config.h"
-#include "ltfat.h"
-
-DEFUN_DLD (comp_dgtreal_fb, args, ,
-  "This function calls the C-library\n\
-  c=comp_dgtreal_fb(f,g,a,M,boundary);\n\
-  Yeah.")
-{
-  const int a = args(2).int_value();
-  const int M = args(3).int_value();
-  const int M2 = M/2+1;
-
-  const Matrix f = args(0).matrix_value();
-  const Matrix g = args(1).matrix_value();
-  
-  const int L  = f.rows();
-  const int W  = f.columns();
-  const int gl = g.rows();
-  const int N = L/a;
-  
-  ComplexMatrix cout(M2,N*W);  
-  
-  dgtreal_fb((double*)f.data(),(double*)g.data(),L,gl,W,a,M,
-	     (ltfat_complex*)cout.data());  
-  
-  return octave_value (cout);
-
-}
-*/
oct/comp_dgtreal_long.cc to oct/comp_sepdgt.cc
--- a/oct/comp_dgtreal_long.cc
+++ b/oct/comp_sepdgt.cc
@@ -1,91 +1,118 @@
 #define TYPEDEPARGS 0, 1
 #define SINGLEARGS
-#define REALARGS
-#define OCTFILENAME comp_dgtreal_long // change to filename
-#define OCTFILEHELP "This function calls the C-library\n\
-                     c=comp_dgtreal_long(f,gf,a,M);\n\
-                     Yeah."
-
+#define COMPLEXINDEPENDENT
+#define OCTFILENAME comp_sepdgt // change to filename
+#define OCTFILEHELP "This function calls the C-library\n c=comp_dgt_fb(f,g,a,M) or comp_dgt_long(f,g,a,M);\n Yeah."
 
 #include "ltfat_oct_template_helper.h"
-// octave_idx_type is 32 or 64 bit signed integer
+// octave_idx_type 32 or 64 bit signed integer
 /*
-  dgtreal_fb forwarders
+  dgt_fb forwarders
 */
 
-static inline void fwd_dgtreal_long(const double *f, const double *g,
-                                  const octave_idx_type L,
-							      const octave_idx_type W, const octave_idx_type a,
-							      const octave_idx_type M, Complex *cout)
+static inline void fwd_dgt_fb(const Complex *f, const Complex *g,
+                              const octave_idx_type L, const octave_idx_type gl,
+                              const octave_idx_type W, const octave_idx_type a,
+                              const octave_idx_type M, Complex *cout)
 {
-   dgtreal_long_d(f,g,L,W,a,M,reinterpret_cast<double _Complex*>(cout));
+   dgt_fb_cd(reinterpret_cast<const double _Complex*>(f),
+             reinterpret_cast<const double _Complex*>(g),
+             L,gl,W,a,M,reinterpret_cast<double _Complex*>(cout));
 }
 
-static inline void fwd_dgtreal_long(const float *f, const float *g,
-                                const octave_idx_type L,
-							    const octave_idx_type W, const octave_idx_type a,
-							    const octave_idx_type M, FloatComplex *cout)
+static inline void fwd_dgt_fb(const FloatComplex *f, const FloatComplex *g,
+                              const octave_idx_type L, const octave_idx_type gl,
+                              const octave_idx_type W, const octave_idx_type a,
+                              const octave_idx_type M, FloatComplex *cout)
 {
-   dgtreal_long_s(f,g,L,W,a,M,reinterpret_cast<float _Complex*>(cout));
+   dgt_fb_cs(reinterpret_cast<const float _Complex*>(f),
+             reinterpret_cast<const float _Complex*>(g),
+             L,gl,W,a,M,reinterpret_cast<float _Complex*>(cout));
 }
 
+static inline void fwd_dgt_fb(const double *f, const double *g,
+                              const octave_idx_type L, const octave_idx_type gl,
+                              const octave_idx_type W, const octave_idx_type a,
+                              const octave_idx_type M, Complex *cout)
+{
+   dgt_fb_d(reinterpret_cast<const double*>(f),
+            reinterpret_cast<const double*>(g),
+            L,gl,W,a,M,reinterpret_cast<double _Complex*>(cout));
+}
+
+static inline void fwd_dgt_fb(const float *f, const float *g,
+                              const octave_idx_type L, const octave_idx_type gl,
+                              const octave_idx_type W, const octave_idx_type a,
+                              const octave_idx_type M, FloatComplex *cout)
+{
+   dgt_fb_s(reinterpret_cast<const float*>(f),
+            reinterpret_cast<const float*>(g),
+            L,gl,W,a,M,reinterpret_cast<float _Complex*>(cout));
+}
+
+static inline void fwd_dgt_long(const Complex *f, const Complex *g,
+                                const octave_idx_type L, const octave_idx_type W,
+                                const octave_idx_type a, const octave_idx_type M,
+                                Complex *cout)
+{
+   dgt_long_cd(reinterpret_cast<const double _Complex*>(f),
+              reinterpret_cast<const double _Complex*>(g),
+              L,W,a,M,reinterpret_cast<double _Complex*>(cout));
+}
+
+static inline void fwd_dgt_long(const FloatComplex *f, const FloatComplex *g,
+                                const octave_idx_type L, const octave_idx_type W, 
+                                const octave_idx_type a, const octave_idx_type M,
+                                FloatComplex *cout)
+{
+   dgt_long_cs(reinterpret_cast<const float _Complex*>(f),
+              reinterpret_cast<const float _Complex*>(g),
+              L,W,a,M,reinterpret_cast<float _Complex*>(cout));
+}
+
+static inline void fwd_dgt_long(const double *f, const double *g,
+                                const octave_idx_type L, const octave_idx_type W,
+                                const octave_idx_type a, const octave_idx_type M,
+                                Complex *cout)
+{
+   dgt_long_d(f,g,L,W,a,M,reinterpret_cast<double _Complex*>(cout));
+}
+
+static inline void fwd_dgt_long(const float *f, const float *g,
+                                const octave_idx_type L, const octave_idx_type W, 
+                                const octave_idx_type a, const octave_idx_type M,
+                                FloatComplex *cout)
+{
+   dgt_long_s(f,g,L,W,a,M,reinterpret_cast<float _Complex*>(cout));
+}
 template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
 octave_value_list octFunction(const octave_value_list& args, int nargout)
 {
-     DEBUGINFO;
-	 const octave_idx_type a = args(2).int_value();
-     const octave_idx_type M = args(3).int_value();
+   DEBUGINFO;
+   const octave_idx_type a = args(2).int_value();
+   const octave_idx_type M = args(3).int_value();
+         
+   MArray<LTFAT_TYPE> f = ltfatOctArray<LTFAT_TYPE>(args(0));
+   MArray<LTFAT_TYPE> g = ltfatOctArray<LTFAT_TYPE>(args(1));
+   const octave_idx_type L  = f.rows();
+   const octave_idx_type W  = f.columns();
+   const octave_idx_type gl = g.rows();
+   const octave_idx_type N = L/a;
+         
+   dim_vector dims_out(M,N,W);  
+   dims_out.chop_trailing_singletons();
 
-	 MArray<LTFAT_TYPE> f = ltfatOctArray<LTFAT_TYPE>(args(0));
-	 MArray<LTFAT_TYPE> g = ltfatOctArray<LTFAT_TYPE>(args(1));
-  
-     const octave_idx_type L = f.rows();
-     const octave_idx_type W = f.columns();
-     const octave_idx_type N = L/a;
-
-     const octave_idx_type M2 = (M/2)+1;
-	 
-	 dim_vector dims_out(M2,N*W);  
-     dims_out.chop_trailing_singletons();
-
-     MArray<LTFAT_COMPLEX> cout(dims_out); 
-     cout.fill(0);
-	 
-	 fwd_dgtreal_long(f.data(),g.data(),L,W,a,M,cout.fortran_vec());
-	 
-     return octave_value(cout);
+   MArray<LTFAT_COMPLEX> cout(dims_out); 
+   cout.fill(0);
+   
+   if(gl<L)
+   {  
+      fwd_dgt_fb(f.data(),g.data(),L,gl,W,a,M,cout.fortran_vec());
+   }
+   else
+   {
+      fwd_dgt_long(f.data(),g.data(),L,W,a,M,cout.fortran_vec());
+   }
+         
+   return octave_value(cout);
 }
-
-
-/* OLD CODE
-#include <octave/oct.h>
-#include "config.h"
-#include "ltfat.h"
-
-DEFUN_DLD (comp_dgtreal_long, args, ,
-  "This function calls the C-library\n\
-  c=comp_dgtreal_long(f,gf,a,M);\n\
-  Yeah.")
-{
-
-  const int a = args(2).int_value();
-  const int M = args(3).int_value();
-
-  const Matrix f = args(0).matrix_value();
-  const Matrix g = args(1).matrix_value();
-  
-  const int L = f.rows();
-  const int W = f.columns();
-  const int N = L/a;
-
-  const int M2 = (M/2)+1;
-  
-  ComplexMatrix cout(M2,N*W);  
-  
-  dgtreal_long((double*)f.data(),(double*)g.data(),
-	    L, W, a, M,(ltfat_complex*)cout.data());
-  
-  return octave_value (cout);
-
-}
-*/
oct/comp_dwilt_long.cc to src/wmdct.c
--- a/oct/comp_dwilt_long.cc
+++ b/src/wmdct.c
@@ -1,138 +1,170 @@
-#define TYPEDEPARGS 0, 1
-#define SINGLEARGS
-#define COMPLEXINDEPENDENT
-#define OCTFILENAME comp_dwilt_long // change to filename
-#define OCTFILEHELP "This function calls the C-library\n\
-                     coef=comp_dwilt_long(f,g,M,L);\n\
-                     Yeah."
+#include "config.h"
+#include "ltfat.h"
+#include "ltfat_types.h"
+
+#define CH(name) LTFAT_COMPLEXH_NAME(name)
+
+#define POSTPROC_REAL \
+  for (Lint n=0;n<N*W;n+=2) \
+  { \
+     for (Lint m=0;m<M;m+=2) \
+     { \
+       pcoef[m]=CH(creal)(pcoef2[m])+CH(cimag)(pcoef2[m]); \
+       pcoef[m+M]=CH(creal)(pcoef2[m+M2])-CH(cimag)(pcoef2[m+M2]); \
+     } \
+     \
+     for (Lint m=1;m<M;m+=2) \
+     { \
+       pcoef[m]=CH(creal)(pcoef2[m])-CH(cimag)(pcoef2[m]); \
+       pcoef[m+M]=CH(creal)(pcoef2[m+M2])+CH(cimag)(pcoef2[m+M2]); \
+     } \
+ \
+     pcoef+=M2; \
+     pcoef2+=M4; \
+  }
+
+#define POSTPROC_COMPLEX \
+  for (Lint n=0;n<N*W;n+=2) \
+  { \
+     for (Lint m=0;m<M;m+=2) \
+     { \
+         pcoef[m] =   scalconst*(emipi4*pcoef2[m]  +eipi4*pcoef2[M2-1-m]); \
+         pcoef[m+M] = scalconst*(eipi4*pcoef2[m+M2]+emipi4*pcoef2[M4-1-m]); \
+     } \
+ \
+     for (Lint m=1;m<M;m+=2) \
+     { \
+       pcoef[m] = scalconst*(eipi4*pcoef2[m]    +emipi4*pcoef2[M2-1-m]); \
+       pcoef[m+M]=scalconst*(emipi4*pcoef2[m+M2]+eipi4*pcoef2[M4-1-m]); \
+     } \
+ \
+     pcoef+=M2; \
+     pcoef2+=M4; \
+  }
+
+#define PREPROC \
+   for(Lint n=0;n<L;n++) \
+      f2[n] = (LTFAT_COMPLEX) cexp(-PI*I*n/(2.0*M)); \
+   for(Lint w=W-1;w>=0;w--) \
+      for(Lint n=0;n<L;n++) \
+         f2[n+w*L] = f2[n]*f[n+w*L];
 
 
-#include "ltfat_oct_template_helper.h"
-// octave_idx_type is 32 or 64 bit signed integer
-/*
-  dgtreal_ola forwarders
-*/
+LTFAT_EXTERN void
+LTFAT_NAME_COMPLEX(dwiltiii_long)(const LTFAT_COMPLEX *f, const LTFAT_COMPLEX *g,
+			   const Lint L, const Lint W, const Lint M,
+			   LTFAT_COMPLEX *cout)
+{
+   const Lint N=L/M;
+   const Lint M2=2*M;
+   const Lint M4=4*M;
+   const LTFAT_REAL scalconst = 1.0/sqrt(2.0);
+   const LTFAT_COMPLEX eipi4 = cexp(I*PI/4.0);
+   const LTFAT_COMPLEX emipi4 = cexp(-I*PI/4.0);
 
-static inline void fwd_dwilt_long(const Complex *f, const Complex *g,
-                                  const octave_idx_type L, 
-							      const octave_idx_type W, 
-							      const octave_idx_type M, 
-								  Complex *cout)
-{
-   dwilt_long_cd(reinterpret_cast<const double _Complex*>(f),
-                reinterpret_cast<const double _Complex*>(g),
-				L,W,M,
-				reinterpret_cast<double _Complex*>(cout));
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_malloc(2*M*N*W*sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *f2 = (LTFAT_COMPLEX*)ltfat_malloc(L*W*sizeof(LTFAT_COMPLEX));
+
+   PREPROC
+
+   LTFAT_NAME_COMPLEX(dgt_long)(f2, g, L, W, M, 2*M, coef2);
+
+   LTFAT_COMPLEX *pcoef  = cout;
+   LTFAT_COMPLEX *pcoef2 = coef2;
+
+   POSTPROC_COMPLEX
+
+   LTFAT_SAFEFREEALL(coef2,f2);
 }
 
-static inline void fwd_dwilt_long(const FloatComplex *f, 
-                                  const FloatComplex *g,
-                                  const octave_idx_type L, 
-							      const octave_idx_type W, 
-							      const octave_idx_type M, 
-								  FloatComplex *cout)
+LTFAT_EXTERN void
+LTFAT_NAME_REAL(dwiltiii_long)(const LTFAT_REAL *f, const LTFAT_REAL *g,
+			   const Lint L, const Lint W, const Lint M,
+			   LTFAT_REAL *cout)
 {
-   dwilt_long_cs(reinterpret_cast<const float _Complex*>(f),
-                reinterpret_cast<const float _Complex*>(g),
-				L,W,M,
-				reinterpret_cast<float _Complex*>(cout));
+   const Lint N=L/M;
+   const Lint M2 = 2*M;
+   const Lint M4=4*M;
+
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_malloc(2*M*N*W*sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *f2 = (LTFAT_COMPLEX*)ltfat_malloc(L*W*sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *g2 = (LTFAT_COMPLEX*)ltfat_malloc(L*sizeof(LTFAT_COMPLEX));
+   for(Lint ii=0;ii<L;ii++)
+       g2[ii]=g[ii];
+
+   PREPROC
+
+   LTFAT_NAME_COMPLEX(dgt_long)(f2, g2, L, W, M, 2*M, coef2);
+
+
+   LTFAT_REAL *pcoef  = cout;
+   LTFAT_COMPLEX *pcoef2 = coef2;
+
+   POSTPROC_REAL
+
+   LTFAT_SAFEFREEALL(coef2,f2,g2);
+
 }
 
-static inline void fwd_dwilt_long(const double *f, const double *g,
-                                  const octave_idx_type L, 
-							      const octave_idx_type W, 
-							      const octave_idx_type M, 
-								  double *cout)
+LTFAT_EXTERN void
+LTFAT_NAME_COMPLEX(dwiltiii_fb)(const LTFAT_COMPLEX *f, const LTFAT_COMPLEX *g,
+			    const Lint L, const Lint gl, const Lint W, const Lint M,
+			   LTFAT_COMPLEX *cout)
 {
-   dwilt_long_d(f,g,L,W,M,cout);
+   const int N=L/M;
+   const int M2=2*M;
+   const int M4=4*M;
+   const LTFAT_REAL scalconst = 1.0/sqrt(2.0);
+   const LTFAT_COMPLEX eipi4 = cexp(I*PI/4.0);
+   const LTFAT_COMPLEX emipi4 = cexp(-I*PI/4.0);
+
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_malloc(2*M*N*W*sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *f2 = (LTFAT_COMPLEX*)ltfat_malloc(L*W*sizeof(LTFAT_COMPLEX));
+
+  PREPROC
+
+  /* coef2=comp_dgt(f,g,a,2*M,L); */
+  LTFAT_NAME_COMPLEX(dgt_fb)(f2, g, L, gl, W, M, 2*M, coef2);
+
+
+  LTFAT_COMPLEX *pcoef  = cout;
+  LTFAT_COMPLEX *pcoef2 = coef2;
+
+  POSTPROC_COMPLEX
+
+  LTFAT_SAFEFREEALL(coef2,f2);
+
 }
 
-static inline void fwd_dwilt_long(const float *f, const float *g,
-                                  const octave_idx_type L, 
-							      const octave_idx_type W, 
-							      const octave_idx_type M, 
-								  float *cout)
+LTFAT_EXTERN void
+LTFAT_NAME_REAL(dwiltiii_fb)(const LTFAT_REAL *f, const LTFAT_REAL *g,
+			   const Lint L, const Lint gl, const Lint W, const Lint M,
+			   LTFAT_REAL *cout)
 {
-   dwilt_long_s(f,g,L,W,M,cout);
+   const Lint N = L/M;
+   const Lint M2 = 2*M;
+   const Lint M4=4*M;
+
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_malloc(2*M*N*W*sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *f2 = (LTFAT_COMPLEX*)ltfat_malloc(L*W*sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *g2 = (LTFAT_COMPLEX*)ltfat_malloc(gl*sizeof(LTFAT_COMPLEX));
+   for(Lint ii=0;ii<gl;ii++)
+       g2[ii]=g[ii];
+
+   PREPROC
+
+   LTFAT_NAME_COMPLEX(dgt_fb)(f2, g2, L, gl, W, M, 2*M, coef2);
+
+   LTFAT_REAL* pcoef  = cout;
+   LTFAT_COMPLEX* pcoef2 = coef2;
+
+   POSTPROC_REAL
+
+   LTFAT_SAFEFREEALL(coef2,f2,g2);
 }
 
-template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
-octave_value_list octFunction(const octave_value_list& args, int nargout)
-{
-     DEBUGINFO;
-     const octave_idx_type M = args(2).int_value();
-     const octave_idx_type L = args(3).int_value();
-     const octave_idx_type N = L/M;
+#undef CH
+#undef POSTPROC_REAL
+#undef POSTPROC_COMPLEX
+#undef PREPROC
 
-	 MArray<LTFAT_TYPE> f = ltfatOctArray<LTFAT_TYPE>(args(0));
-	 MArray<LTFAT_TYPE> g = ltfatOctArray<LTFAT_TYPE>(args(1));
-  
-     const int W = f.columns();
-
-     dim_vector dims_out(2*M,N/2,W);  
-     dims_out.chop_trailing_singletons();
-
-     MArray<LTFAT_TYPE> cout(dims_out); 
-     cout.fill(0);
-	 
-	 fwd_dwilt_long(f.data(),g.data(),L, W, M,cout.fortran_vec());
-	 
-     return octave_value(cout);
-}
-
-/*
-#include <octave/oct.h>
-#include "config.h"
-#include "ltfat.h"
-
-
-DEFUN_DLD (comp_dwilt_long, args, ,
-  "This function calls the C-library\n\
-  coef=comp_dwilt_long(f,g,M,L);\n\
-  Yeah.")
-{
-  const int M = args(2).int_value();
-  const int L = args(3).int_value();
-  const int N = L/M;
-
-  if (args(0).is_complex_type())
-  {
-
-    const ComplexMatrix f = args(0).complex_matrix_value();
-    const ComplexMatrix g = args(1).complex_matrix_value();
-    
-    const int W = f.columns();
-
-    dim_vector dims_out(2*M,N/2,W);  
-    dims_out.chop_trailing_singletons();
-    
-    ComplexNDArray cout(dims_out);  
-    
-    dwilt_long((ltfat_complex*)f.data(),(ltfat_complex*)g.data(), L, W, M,
-	       (ltfat_complex*)cout.data());
-
-    return octave_value (cout);
-
-  }
-  else
-  {
-
-    const Matrix f = args(0).matrix_value();
-    const Matrix g = args(1).matrix_value();
-    
-    const int W = f.columns();
-
-    dim_vector dims_out(2*M,N/2,W);  
-    dims_out.chop_trailing_singletons();
-    
-    NDArray cout(dims_out);  
-    
-    dwiltreal_long((double*)f.data(),(double*)g.data(), L, W, M,
-		  (double*)cout.data());
-	
-    return octave_value (cout);
-
-  }
-
-}
-*/
oct/comp_idgt_fac.cc to src/iwmdct.c
--- a/oct/comp_idgt_fac.cc
+++ b/src/iwmdct.c
@@ -1,82 +1,164 @@
-#define TYPEDEPARGS 0, 1
-#define SINGLEARGS
-#define COMPLEXARGS
-#define OCTFILENAME comp_idgt_fac // change to filename
-#define OCTFILEHELP "This function calls the C-library\n\
-                     c=comp_idgt_fac(c,gf,L,a,M);\n\
-                     Yeah.\n"
+#include "config.h"
+#include "ltfat.h"
+#include "ltfat_types.h"
+
+#define CH(name) LTFAT_COMPLEXH_NAME(name)
+
+#define PREPROC_COMPLEX \
+  for (Lint n=0;n<N*W;n+=2) \
+  { \
+     for (Lint m=0;m<M;m+=2) \
+     { \
+        pcoef2[m] = eipi4*pcoef[m]; \
+        pcoef2[M2-1-m] = emipi4*pcoef[m]; \
+        pcoef2[m+M2] = emipi4*pcoef[m+M]; \
+        pcoef2[M4-1-m] = eipi4*pcoef[m+M]; \
+     } \
+ \
+     for (Lint m=1;m<M;m+=2) \
+     { \
+        pcoef2[m] = emipi4*pcoef[m]; \
+        pcoef2[M2-1-m] = eipi4*pcoef[m]; \
+        pcoef2[m+M2] = eipi4*pcoef[m+M];  \
+        pcoef2[M4-1-m] = emipi4*pcoef[m+M]; \
+     } \
+ \
+     pcoef+=M2; \
+     pcoef2+=M4; \
+  }
+
+#define POSTPROC_REAL \
+   for(Lint w=0;w<W;w++) \
+      for(Lint n=0;n<L;n++) \
+         f[n+w*L] = scalconst*CH(creal)(f2[n+w*L]*CH(cexp)(I*PI*n/(2.0*M)));
+
+#define POSTPROC_COMPLEX \
+   for(Lint w=0;w<W;w++) \
+      for(Lint n=0;n<L;n++) \
+         f[n+w*L] = scalconst*f2[n+w*L]*CH(cexp)(I*PI*n/(2.0*M));
+
+LTFAT_EXTERN void
+LTFAT_NAME_COMPLEX(idwiltiii_long)(const LTFAT_COMPLEX *c, const LTFAT_COMPLEX *g,
+			   const Lint L, const Lint W, const Lint M,
+			   LTFAT_COMPLEX *f)
+{
+   const Lint N=L/M;
+   const Lint M2=2*M;
+   const Lint M4=4*M;
+   const LTFAT_REAL scalconst = 1.0/sqrt(2.0);
+   const LTFAT_COMPLEX eipi4 = cexp(I*PI/4.0);
+   const LTFAT_COMPLEX emipi4 = cexp(-I*PI/4.0);
+
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_calloc(2*M*N*W,sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *f2 = (LTFAT_COMPLEX*)ltfat_malloc(L*W*sizeof(LTFAT_COMPLEX));
 
 
-#include "ltfat_oct_template_helper.h"
-// octave_idx_type is 32 or 64 bit signed integer
+   LTFAT_COMPLEX *pcoef  = c;
+   LTFAT_COMPLEX *pcoef2 = coef2;
 
-static inline void fwd_idgt_fac(const Complex *coef, const Complex *gf,
-                              const octave_idx_type L, const octave_idx_type W,
-							  const octave_idx_type a, const octave_idx_type M,
-							  Complex *f)
-{
-   idgt_fac_d(reinterpret_cast<const double _Complex*>(coef),
-              reinterpret_cast<const double _Complex*>(gf),
-              L,W,a,M,
-			  reinterpret_cast<double _Complex*>(f));
+   PREPROC_COMPLEX
+
+   LTFAT_NAME(idgt_long)(coef2, g, L, W, M, 2*M, f2);
+   
+   POSTPROC_COMPLEX
+
+   LTFAT_SAFEFREEALL(coef2,f2);
 }
 
-static inline void fwd_idgt_fac(const FloatComplex *coef, const FloatComplex *gf,
-                              const octave_idx_type L, const octave_idx_type W, 
-							  const octave_idx_type a, const octave_idx_type M,
-							  FloatComplex *f)
+LTFAT_EXTERN void
+LTFAT_NAME_REAL(idwiltiii_long)(const LTFAT_REAL *c, const LTFAT_REAL *g,
+			   const Lint L, const Lint W, const Lint M,
+			   LTFAT_REAL *f)
 {
-   idgt_fac_s(reinterpret_cast<const float _Complex*>(coef),
-              reinterpret_cast<const float _Complex*>(gf),
-              L,W,a,M,
-			  reinterpret_cast<float _Complex*>(f));
+   const Lint N=L/M;
+   const Lint M2 = 2*M;
+   const Lint M4=4*M;
+   const LTFAT_REAL scalconst = 1.0/sqrt(2.0);
+   const LTFAT_COMPLEX eipi4 = cexp(I*PI/4.0);
+   const LTFAT_COMPLEX emipi4 = cexp(-I*PI/4.0);
+
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_calloc(2*M*N*W,sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *f2 = (LTFAT_COMPLEX*)ltfat_malloc(L*W*sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *g2 = (LTFAT_COMPLEX*)ltfat_malloc(L*sizeof(LTFAT_COMPLEX));
+   for(Lint ii=0;ii<L;ii++)
+       g2[ii]=g[ii];
+
+
+   LTFAT_REAL *pcoef  = c;
+   LTFAT_COMPLEX *pcoef2 = coef2;
+
+   PREPROC_COMPLEX
+
+   LTFAT_NAME(idgt_long)(coef2, g2, L, W, M, 2*M, f2);
+
+   POSTPROC_REAL
+
+   LTFAT_SAFEFREEALL(coef2,f2,g2);
+
 }
 
-template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
-octave_value_list octFunction(const octave_value_list& args, int nargout)
+LTFAT_EXTERN void
+LTFAT_NAME_COMPLEX(idwiltiii_fb)(const LTFAT_COMPLEX *c, const LTFAT_COMPLEX *g,
+			    const Lint L, const Lint gl, const Lint W, const Lint M,
+			   LTFAT_COMPLEX *f)
 {
-     DEBUGINFO;
-	 
-	 MArray<LTFAT_TYPE> coef = ltfatOctArray<LTFAT_TYPE>(args(0));
-	 MArray<LTFAT_TYPE> gf = ltfatOctArray<LTFAT_TYPE>(args(1));
-	 const octave_idx_type L = args(2).int_value();
-     const octave_idx_type a = args(3).int_value();
-     const octave_idx_type M = args(4).int_value();
-     const octave_idx_type N = L/a;
-     const octave_idx_type W = coef.rows()*coef.columns()/(M*N);
+   const int N=L/M;
+   const int M2=2*M;
+   const int M4=4*M;
+   const LTFAT_REAL scalconst = 1.0/sqrt(2.0);
+   const LTFAT_COMPLEX eipi4 = cexp(I*PI/4.0);
+   const LTFAT_COMPLEX emipi4 = cexp(-I*PI/4.0);
+
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_calloc(2*M*N*W,sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *f2 = (LTFAT_COMPLEX*)ltfat_malloc(L*W*sizeof(LTFAT_COMPLEX));
 
 
-     MArray<LTFAT_COMPLEX> f(dim_vector(L,W)); 
-     f.fill(0);
-	 
-	 fwd_idgt_fac(coef.data(),gf.data(),L,W,a,M,f.fortran_vec());
-	 
-     return octave_value(f);
+   LTFAT_COMPLEX *pcoef  = c;
+   LTFAT_COMPLEX *pcoef2 = coef2;
+
+   PREPROC_COMPLEX
+   
+   LTFAT_NAME(idgt_fb)(coef2, g, L, gl, W, M, 2*M, f2);
+   
+   POSTPROC_COMPLEX
+
+   LTFAT_SAFEFREEALL(coef2,f2);
+
 }
-/*
-#include <octave/oct.h>
-#include "config.h"
-#include "ltfat.h"
 
-DEFUN_DLD (comp_idgt_fac, args, ,
-  "This function calls the C-library\n\
-  c=comp_idgt_fac(c,gf,L,a,M);\n\
-  Yeah.\n")
+LTFAT_EXTERN void
+LTFAT_NAME_REAL(idwiltiii_fb)(const LTFAT_REAL *c, const LTFAT_REAL *g,
+			   const Lint L, const Lint gl, const Lint W, const Lint M,
+			   LTFAT_REAL *f)
 {
+   const Lint N = L/M;
+   const Lint M2 = 2*M;
+   const Lint M4=4*M;
+   const LTFAT_REAL scalconst = 1.0/sqrt(2.0);
+   const LTFAT_COMPLEX eipi4 = cexp(I*PI/4.0);
+   const LTFAT_COMPLEX emipi4 = cexp(-I*PI/4.0);
 
-  const ComplexMatrix coef = args(0).complex_matrix_value();
-  const ComplexMatrix gf = args(1).complex_matrix_value();
-  const int L = args(2).int_value();
-  const int a = args(3).int_value();
-  const int M = args(4).int_value();
-  const int N = L/a;
-  const int W = coef.rows()*coef.columns()/(M*N);
+   LTFAT_COMPLEX *coef2 = (LTFAT_COMPLEX*)ltfat_calloc(2*M*N*W,sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *f2 = (LTFAT_COMPLEX*)ltfat_malloc(L*W*sizeof(LTFAT_COMPLEX));
+   LTFAT_COMPLEX *g2 = (LTFAT_COMPLEX*)ltfat_malloc(gl*sizeof(LTFAT_COMPLEX));
+   for(Lint ii=0;ii<gl;ii++)
+       g2[ii]=g[ii];
 
-  ComplexMatrix f(L,W);  
-  
-  idgt_fac((ltfat_complex*)coef.data(),(ltfat_complex*)gf.data(),L,W,a,M,
-	  (ltfat_complex*)f.data());
 
-  return octave_value (f);
+   LTFAT_REAL* pcoef  = c;
+   LTFAT_COMPLEX* pcoef2 = coef2;
+
+   PREPROC_COMPLEX
+   
+   LTFAT_NAME(idgt_fb)(coef2, g2, L, gl, W, M, 2*M, f2);
+   
+   POSTPROC_REAL
+   
+   LTFAT_SAFEFREEALL(coef2,f2,g2);
 }
-*/
+
+#undef CH
+#undef PREPROC_COMPLEX
+#undef POSTPROC_REAL
+#undef POSTPROC_COMPLEX
+
oct/comp_idgt_fb.cc to mex/comp_dwiltiii.c
--- a/oct/comp_idgt_fb.cc
+++ b/mex/comp_dwiltiii.c
@@ -1,102 +1,67 @@
+#ifndef _LTFAT_MEX_FILE
+#define _LTFAT_MEX_FILE
+
+#define ISNARGINEQ 3
 #define TYPEDEPARGS 0, 1
 #define SINGLEARGS
-#define COMPLEXARGS
-#define OCTFILENAME comp_idgt_fb // change to filename
-#define OCTFILEHELP "This function calls the C-library\n\
-                    c=comp_idgt_fb(coef,g,L,a,M);\n\
-                    Yeah."
+#define COMPLEXINDEPENDENT
 
+#endif // _LTFAT_MEX_FILE - INCLUDED ONCE
 
-#include "ltfat_oct_template_helper.h"
-// octave_idx_type is 32 or 64 bit signed integer
+#define MEX_FILE __BASE_FILE__
+#include "ltfat_mex_template_helper.h"
 
-static inline void fwd_idgt_fb(const Complex *coef, const Complex *gf,
-                               const octave_idx_type L,
-                               const octave_idx_type gl, const octave_idx_type W,
-							   const octave_idx_type a, const octave_idx_type M,
-							   Complex *f)
+#if defined(LTFAT_SINGLE) || defined(LTFAT_DOUBLE)
+#include "ltfat_types.h"
+
+// Calling convention:
+//  comp_dwiltiii(f,g,M);
+
+void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
 {
-   idgt_fb_d(reinterpret_cast<const double _Complex*>(coef),
-              reinterpret_cast<const double _Complex*>(gf),
-              L,gl,W,a,M,
-			  reinterpret_cast<double _Complex*>(f));
-}
+   int M, N, L, gl, W;
 
-static inline void fwd_idgt_fb(const FloatComplex *coef, const FloatComplex *gf,
-                               const octave_idx_type L,
-                               const octave_idx_type gl, const octave_idx_type W, 
-							   const octave_idx_type a, const octave_idx_type M,
-							   FloatComplex *f)
-{
-   idgt_fb_s(reinterpret_cast<const float _Complex*>(coef),
-              reinterpret_cast<const float _Complex*>(gf),
-              L,gl,W,a,M,
-			  reinterpret_cast<float _Complex*>(f));
-}
+   mwSize ndim;
+   mwSize dims[3];
 
-template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
-octave_value_list octFunction(const octave_value_list& args, int nargout)
-{
-     DEBUGINFO;
-	 
-	 MArray<LTFAT_TYPE> coef = ltfatOctArray<LTFAT_TYPE>(args(0));
-	 MArray<LTFAT_TYPE> gf = ltfatOctArray<LTFAT_TYPE>(args(1));
-	 const octave_idx_type L = args(2).int_value();
-     const octave_idx_type a = args(3).int_value();
-     const octave_idx_type M = args(4).int_value();
-     const octave_idx_type N = L/a;
-	 const octave_idx_type gl = gf.rows();
+   // Get matrix dimensions.
+   M=(int)mxGetScalar(prhs[2]);
+   L=(int)mxGetM(prhs[0]);
+   gl=(int) mxGetM(prhs[1]);
+   W = mxGetN(prhs[0]);
 
-     const octave_idx_type W = coef.rows()*coef.columns()/(M*N);
+   N=L/M;
 
+   dims[0]=M;
+   dims[1]=N;
 
-     MArray<LTFAT_COMPLEX> f(dim_vector(L,W)); 
-     f.fill(0);
-	 
-	 fwd_idgt_fb(coef.data(),gf.data(),L,gl,W,a,M,f.fortran_vec());
-	 
-     return octave_value(f);
-}
-/*
-#include <octave/oct.h>
-#include "config.h"
-#include "ltfat.h"
-
-DEFUN_DLD (comp_idgt_fb, args, ,
-  "This function calls the C-library\n\
-  c=comp_idgt_fb(coef,g,L,a,M);\n\
-  Yeah.")
-{   
-   const ComplexMatrix coef = args(0).complex_matrix_value();
-
-   const int L = args(2).int_value();
-   const int a = args(3).int_value();
-   const int M = args(4).int_value();
-   const int N = L/a;
-
-   const int W = coef.rows()*coef.columns()/(M*N);
-   
-   ComplexMatrix f(L,W);  
-   
-   if (args(1).is_complex_type())
+   if (W==1)
    {
-      const ComplexMatrix g = args(1).complex_matrix_value();
-      const int gl = g.rows();
-      
-      idgt_fb((const ltfat_complex*)coef.data(),(const ltfat_complex*)g.data(),
-	      L,gl,W,a,M,
-	      (ltfat_complex*)f.fortran_vec());        
+      ndim=2;
    }
    else
    {
-      const Matrix g = args(1).matrix_value();
-      const int gl = g.rows();
-      
-      idgt_fb_r((const ltfat_complex*)coef.data(),g.data(),
-		L,gl,W,a,M,
-		(ltfat_complex*)f.fortran_vec());      
+      ndim=3;
+      dims[2]=W;
    }
 
-   return octave_value (f);               
+   plhs[0] = ltfatCreateNdimArray(ndim,dims,LTFAT_MX_CLASSID,LTFAT_MX_COMPLEXITY);
+
+   const LTFAT_TYPE* f = (const LTFAT_TYPE*) mxGetData(prhs[0]);
+   const LTFAT_TYPE* g = (const LTFAT_TYPE*) mxGetData(prhs[1]);
+   LTFAT_TYPE* cout = (LTFAT_TYPE*) mxGetData(plhs[0]);
+
+   if(gl<L)
+   {
+      LTFAT_NAME(dwiltiii_fb)(f,g,L,gl, W, M, cout);
+   }
+   else
+   {
+      LTFAT_NAME(dwiltiii_long)(f,g,L,W,M,cout);
+   }
+
+
+
+   return;
 }
-*/
+#endif
oct/comp_idgtreal_fac.cc to mex/comp_idwilt.c
--- a/oct/comp_idgtreal_fac.cc
+++ b/mex/comp_idwilt.c
@@ -1,80 +1,48 @@
+#ifndef _LTFAT_MEX_FILE
+#define _LTFAT_MEX_FILE
+
+#define ISNARGINEQ 2
 #define TYPEDEPARGS 0, 1
 #define SINGLEARGS
-#define COMPLEXARGS
-#define OCTFILENAME comp_idgtreal_fac // change to filename
-#define OCTFILEHELP "This function calls the C-library\n\
-                    c=comp_idgtreal_fac(c,gf,L,W,a,M);\n\
-                    Yeah.\n"
+#define COMPLEXINDEPENDENT
 
+#endif // _LTFAT_MEX_FILE - INCLUDED ONCE
 
-#include "ltfat_oct_template_helper.h"
-// octave_idx_type is 32 or 64 bit signed integer
+#define MEX_FILE __BASE_FILE__
+#include "ltfat_mex_template_helper.h"
 
-static inline void fwd_idgtreal_fac(const Complex *coef, const Complex *gf,
-                               const octave_idx_type L, const octave_idx_type W,
-							   const octave_idx_type a, const octave_idx_type M,
-							   double *f)
+#if defined(LTFAT_SINGLE) || defined(LTFAT_DOUBLE)
+#include "ltfat_types.h"
+
+// Calling convention:
+//  comp_idwilt(c,g);
+
+void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
 {
-   idgtreal_fac_d(reinterpret_cast<const double _Complex*>(coef),
-              reinterpret_cast<const double _Complex*>(gf),
-              L,W,a,M,f);
+   int M, N, L, gl, W;
+
+   // Get matrix dimensions.
+   M = mxGetM(prhs[0])/2;
+   // This is woraround to het number of columns
+   N = 2*ltfatGetN(prhs[0]);
+   gl= mxGetNumberOfElements(prhs[1]);
+   W = mxGetNumberOfElements(prhs[0])/(M*N);
+
+   L=N*M;
+
+   plhs[0] = ltfatCreateMatrix(L,W,LTFAT_MX_CLASSID,LTFAT_MX_COMPLEXITY);
+
+   const LTFAT_TYPE* c = (const LTFAT_TYPE*) mxGetData(prhs[0]);
+   const LTFAT_TYPE* g = (const LTFAT_TYPE*) mxGetData(prhs[1]);
+   LTFAT_TYPE* f = (LTFAT_TYPE*) mxGetData(plhs[0]);
+
+   if(gl<L)
+   {
+      LTFAT_NAME(idwilt_fb)(c,g,L,gl,W,M,f);
+   }
+   else
+   {
+      LTFAT_NAME(idwilt_long)(c,g,L,W,M,f);
+   }
 }
-
-static inline void fwd_idgtreal_fac(const FloatComplex *coef, const FloatComplex *gf,
-                               const octave_idx_type L, const octave_idx_type W, 
-							   const octave_idx_type a, const octave_idx_type M,
-							   float *f)
-{
-   idgtreal_fac_s(reinterpret_cast<const float _Complex*>(coef),
-              reinterpret_cast<const float _Complex*>(gf),
-              L,W,a,M,f);
-}
-
-template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
-octave_value_list octFunction(const octave_value_list& args, int nargout)
-{
-     DEBUGINFO;
-	 
-	 MArray<LTFAT_TYPE> coef = ltfatOctArray<LTFAT_TYPE>(args(0));
-	 MArray<LTFAT_TYPE> gf = ltfatOctArray<LTFAT_TYPE>(args(1));
-	 const octave_idx_type L = args(2).int_value();
-     const octave_idx_type a = args(3).int_value();
-     const octave_idx_type M = args(4).int_value();
-     const octave_idx_type N = L/a;
-     const octave_idx_type W = coef.columns()/N;
-
-     MArray<LTFAT_REAL> f(dim_vector(L,W)); 
-     f.fill(0);
-	 
-	 fwd_idgtreal_fac(coef.data(),gf.data(),L,W,a,M,f.fortran_vec());
-	 
-     return octave_value(f);
-}
-
-/*
-#include <octave/oct.h>
-#include "config.h"
-#include "ltfat.h"
-
-DEFUN_DLD (comp_idgtreal_fac, args, ,
-  "This function calls the C-library\n\
-  c=comp_idgt_fac(c,gf,L,a,M);\n\
-  Yeah.\n")
-{
-
-  const ComplexMatrix coef = args(0).complex_matrix_value();
-  const ComplexMatrix gf = args(1).complex_matrix_value();
-  const int L = args(2).int_value();
-  const int a = args(3).int_value();
-  const int M = args(4).int_value();
-  const int N = L/a;
-  const int W = coef.columns()/N;
-
-  Matrix f(L,W);  
-  
-  idgtreal_fac((ltfat_complex*)coef.data(),(ltfat_complex*)gf.data(),L,W,a,M,
-	  (double*)f.data());
-
-  return octave_value (f);
-}
-*/
+#endif
oct/comp_idgtreal_fb.cc to mex/comp_idwiltiii.c
--- a/oct/comp_idgtreal_fb.cc
+++ b/mex/comp_idwiltiii.c
@@ -1,83 +1,47 @@
-#define TYPEDEPARGS 0
+#ifndef _LTFAT_MEX_FILE
+#define _LTFAT_MEX_FILE
+
+#define ISNARGINEQ 2
+#define TYPEDEPARGS 0, 1
 #define SINGLEARGS
-#define MATCHEDARGS 1
-#define COMPLEXARGS
-#define OCTFILENAME comp_idgtreal_fb // change to filename
-#define OCTFILEHELP "This function calls the C-library\n\
-                    c=comp_idgtereal_fb(c,g,L,a,M);\n\
-                    Yeah.\n"
+#define COMPLEXINDEPENDENT
 
+#endif // _LTFAT_MEX_FILE - INCLUDED ONCE
 
-#include "ltfat_oct_template_helper.h"
-// octave_idx_type is 32 or 64 bit signed integer
+#define MEX_FILE __BASE_FILE__
+#include "ltfat_mex_template_helper.h"
 
-static inline void fwd_idgtreal_fb(const Complex *coef, const double *gf,
-                               const octave_idx_type L,
-                               const octave_idx_type gl, const octave_idx_type W,
-							   const octave_idx_type a, const octave_idx_type M,
-							   double *f)
+#if defined(LTFAT_SINGLE) || defined(LTFAT_DOUBLE)
+#include "ltfat_types.h"
+
+// Calling convention:
+//  comp_idwilt(c,g);
+
+void LTFAT_NAME(ltfatMexFnc)( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
 {
-   idgtreal_fb_d(reinterpret_cast<const double _Complex*>(coef),
-                 gf,L,gl,W,a,M,f);
+   int M, N, L, gl, W;
+
+   // Get matrix dimensions.
+   M= mxGetM(prhs[0]);
+   N= ltfatGetN(prhs[0]);
+   gl= mxGetNumberOfElements(prhs[1]);
+   W = mxGetNumberOfElements(prhs[0])/(M*N);
+
+   L=N*M;
+
+   plhs[0] = ltfatCreateMatrix(L,W,LTFAT_MX_CLASSID,LTFAT_MX_COMPLEXITY);
+
+   const LTFAT_TYPE* c = (const LTFAT_TYPE*) mxGetData(prhs[0]);
+   const LTFAT_TYPE* g = (const LTFAT_TYPE*) mxGetData(prhs[1]);
+   LTFAT_TYPE* f = (LTFAT_TYPE*) mxGetData(plhs[0]);
+
+   if(gl<L)
+   {
+      LTFAT_NAME(idwiltiii_fb)(c,g,L,gl,W,M,f);
+   }
+   else
+   {
+      LTFAT_NAME(idwiltiii_long)(c,g,L,W,M,f);
+   }
 }
-
-static inline void fwd_idgtreal_fb(const FloatComplex *coef, const float *gf,
-                               const octave_idx_type L,
-                               const octave_idx_type gl, const octave_idx_type W, 
-							   const octave_idx_type a, const octave_idx_type M,
-							   float *f)
-{
-   idgtreal_fb_s(reinterpret_cast<const float _Complex*>(coef),
-                 gf,L,gl,W,a,M,f);
-}
-
-template <class LTFAT_TYPE, class LTFAT_REAL, class LTFAT_COMPLEX>
-octave_value_list octFunction(const octave_value_list& args, int nargout)
-{
-     DEBUGINFO;
-	 
-	 MArray<LTFAT_TYPE> coef = ltfatOctArray<LTFAT_TYPE>(args(0));
-	 MArray<LTFAT_REAL> gf = ltfatOctArray<LTFAT_REAL>(args(1));
-	 const octave_idx_type L = args(2).int_value();
-     const octave_idx_type a = args(3).int_value();
-     const octave_idx_type M = args(4).int_value();
-     const octave_idx_type N = L/a;
-     const octave_idx_type W = coef.columns()/N;
-     const octave_idx_type gl = gf.rows();
-
-     MArray<LTFAT_REAL> f(dim_vector(L,W)); 
-     f.fill(0);
-	 
-	 fwd_idgtreal_fb(coef.data(),gf.data(),L,gl,W,a,M,f.fortran_vec());
-	 
-     return octave_value(f);
-}
-
-/*
-#include <octave/oct.h>
-#include "config.h"
-#include "ltfat.h"
-
-DEFUN_DLD (comp_idgtreal_fb, args, ,
-  "This function calls the C-library\n\
-  c=comp_idgt_fb(c,g,L,a,M);\n\
-  Yeah.\n")
-{
-
-  const ComplexMatrix coef = args(0).complex_matrix_value();
-  const Matrix g = args(1).matrix_value();
-  const int L = args(2).int_value();
-  const int a = args(3).int_value();
-  const int M = args(4).int_value();
-  const int N = L/a;
-  const int W = coef.columns()/N;
-  const int gl = g.rows();
-
-  Matrix f(L,W);  
-  
-  idgtreal_fb((ltfat_complex*)coef.data(),g.data(),L,gl,W,a,M,
-	  f.fortran_vec());
-
-  return octave_value (f);
-}
-*/
+#endif
src/tfutil.c to comp/comp_sepdgt.m
--- a/src/tfutil.c
+++ b/comp/comp_sepdgt.m
@@ -1,35 +1,20 @@
-#include "ltfat.h"
+function [coef]=comp_sepdgt(f,g,a,M)
+%COMP_SEPDGT  Separable DGT
+%   Usage:  c=comp_sepdgt(f,g,a,M);
+%  
+%   This is a computational routine. Do not call it directly.
+%
+%   See help on DGT.
 
-/*
-int LTFAT_NAME(complexprod)(LTFAT_COMPLEX *c, const LTFAT_COMPLEX a, const LTFAT_COMPLEX b)
-{
-#ifdef HAVE_COMPLEX_H
-  (*c)=a*b;
-#else
+%   AUTHOR : Peter L. Søndergaard.
 
-  (*c)[0]=a[0]*b[0]-a[1]*b[1];
-  (*c)[1]=a[1]*b[0]+a[0]*b[1];
+L=size(f,1);
+Lwindow=size(g,1);
 
-#endif
-  return (0);
-}
-*/
-
-/* --- usefull for debugging --------
-void print_z(const int N, LTFAT_COMPLEX *p)
-{
-  int i;
-
-  for(i=0;i<N;i++)
-  {
-#ifdef HAVE_COMPLEX_H
-    printf("%.16lf  %.16lf\n",creal(p[i]),cimag(p[i]));
-#else
-    printf("%.16lf  %.16lf\n",p[i][0],p[i][1]);
-#endif
-  }
-  fflush(NULL);
-
-}
-
-*/
+if Lwindow<L
+   % Do the filter bank algorithm
+   coef=comp_dgt_fb(f,g,a,M);
+else
+   % Do the factorization algorithm
+   coef=comp_dgt_long(f,g,a,M);
+end;
<< < 1 2 3 4 5 > >> (Page 4 of 5)