Index: src/ltfat.h
===================================================================
--- src/ltfat.h	(revision 879)
+++ src/ltfat.h	(working copy)
@@ -8,8 +8,11 @@
 /* FFTW_ESTIMATE is the only one to work, as the other types destroys input */
 #define FFTW_OPTITYPE FFTW_ESTIMATE
 
+/* Allow using arbitrary complex types. Assumes identical memory layout. */
+#ifndef LTFAT_USER_COMPLEX
 typedef fftw_complex  ltfat_complex;
 typedef fftwf_complex ltfat_scomplex;
+#endif
 
 /* BEGIN_C_DECLS */
 
Index: oct/comp_dgt_fac.cc
===================================================================
--- oct/comp_dgt_fac.cc	(revision 879)
+++ oct/comp_dgt_fac.cc	(working copy)
@@ -2,57 +2,64 @@
 #include "config.h"
 #include "ltfat.h"
 
+template <class MT1, class MT2>
+static octave_value 
+do_comp_dgt_fac (const octave_value_list& args,
+                 const int a, const int M,
+                 void (*lib_func) (typename MT1::element_type *,
+                                   typename MT2::element_type *,
+                                   const int, const int,
+                                   const int, const int,
+                                   const int,
+                                   typename MT2::element_type *))
+{
+ typedef typename MT1::element_type T1;
+ typedef typename MT2::element_type T2;
+
+ const MT1 f = octave_value_extract<MT1> args(0);
+ const MT2 gf = octave_value_extract<MT1> args(1);
+
+ const int L = f.rows();
+ const int W = f.columns();
+ const int R = gf.rows()*gf.columns()/L;
+
+ const int N = L/a;
+
+ MT2 cout(M,N*W*R);
+
+ // const_casts are necessary because LTFAT doesn't specify const on input data.
+ lib_func (const_cast<T1 *>(f.fortran_vec ()), const_cast<T2 *>(gf.fortran_vec ()),
+           L, W, R, a, M, cout.fortran_vec ());
+
+ return cout;
+}
+
+template <class MT1, class MT2>
 DEFUN_DLD (comp_dgt_fac, args, ,
   "This function calls the C-library\n\
   c=comp_dgt_fac(f,gf,a,M);\n\
   Yeah.")
 {
+  octave_value_list retval;
 
   const bool f_is_complex  = args(0).is_complex_type();
 
   const int a = args(2).int_value();
   const int M = args(3).int_value();
 
-  if (f_is_complex)
-  {
-
-     const ComplexMatrix f = args(0).complex_matrix_value();
-     const ComplexMatrix gf = args(1).complex_matrix_value();
-     
-     const int L = f.rows();
-     const int W = f.columns();
-     const int R = gf.rows()*gf.columns()/L;
-     
-     const int N = L/a;
-     
-     ComplexMatrix cout(M,N*W*R);  
-     
-     dgt_fac((ltfat_complex*)f.data(),(ltfat_complex*)gf.data(),
-     L, W, R, a, M,(ltfat_complex*)cout.data());
-
-     return octave_value (cout);
-     
-
-  }
+  if (args(0).is_single_type ())
+    {
+      if (f_is_complex)
+        do_comp_dgt_fac (args, a, M, sdgt_fac);
+      else
+        do_comp_dgt_fac (args, a, M, sdgt_fac_r);
+    }
   else
-  {
+    {
+      if (f_is_complex)
+        do_comp_dgt_fac (args, a, M, dgt_fac);
+      else
+        do_comp_dgt_fac (args, a, M, dgt_fac_r);
+    }
 
-     const Matrix f = args(0).matrix_value();
-     const ComplexMatrix gf = args(1).complex_matrix_value();
-     
-     const int L = f.rows();
-     const int W = f.columns();
-     const int R = gf.rows()*gf.columns()/L;
-     
-     const int N = L/a;
-     
-     ComplexMatrix cout(M,N*W*R);  
-     
-     dgt_fac_r((double*)f.data(),(ltfat_complex*)gf.data(),
-     L, W, R, a, M,(ltfat_complex*)cout.data());
-
-     return octave_value (cout);
-  
-  }
-
 }
Index: oct/config.h
===================================================================
--- oct/config.h	(revision 879)
+++ oct/config.h	(working copy)
@@ -5,13 +5,15 @@
 #define CONFIG_H 1
 
 #include "fftw3.h"
+
+#define LTFAT_USER_COMPLEX
+typedef Complex  ltfat_complex;
+typedef FloatComplex ltfat_scomplex;
+
 #include "ltfat.h"
 
 #define FFTW_OPTITYPE FFTW_ESTIMATE
 
-typedef fftw_complex  ltfat_complex;
-typedef fftwf_complex ltfat_scomplex;
-
 void* ltfat_malloc (size_t n)
 {
   return fftw_malloc(n);