added fftw3 (uncompiled) to the mtx_rfft.c and mtx_rifft.c sources.

Unless DHAVE_FFTW3 is set, everything works like before.

svn path=/trunk/externals/iem/iemmatrix/; revision=10202

6 files changed, 213 insertions, 36 deletions

diff --git a/doc/mtx_rfft-help.pd b/doc/mtx_rfft-help.pd
index 2d20bc0..6548a01 100644
--- a/doc/mtx_rfft-help.pd
+++ b/doc/mtx_rfft-help.pd
@@ -18,17 +18,19 @@ has to have 2^k columns \; a N/2+1 columns matrix is the result \;
 #X obj 172 255 / 8;
 #X obj 49 306 t a a;
 #X obj 92 306 mtx_print original;
-#X obj 49 394 mtx_print rfft_realpart;
-#X obj 121 372 mtx_print rfft_imagpart;
+#X obj 49 414 mtx_print rfft_realpart;
+#X obj 121 392 mtx_print rfft_imagpart;
 #X text 266 419 see also:;
 #X obj 49 280 mtx_cos;
 #X obj 103 255 * 3.14159;
 #X text 109 280 cosine wave;
-#X obj 121 352 mtx_int;
-#X obj 49 373 mtx_int;
+#X obj 121 372 mtx_int;
+#X obj 49 393 mtx_int;
 #X text 89 15 [mtx_rfft];
 #X obj 49 330 mtx_rfft;
 #X obj 261 438 mtx_rifft;
+#X obj 121 353 mtx_+ 0.5;
+#X obj 49 373 mtx_+ 0.5;
 #X connect 4 0 8 0;
 #X connect 7 0 10 0;
 #X connect 8 0 9 0;
@@ -43,5 +45,7 @@ has to have 2^k columns \; a N/2+1 columns matrix is the result \;
 #X connect 21 0 11 1;
 #X connect 23 0 18 0;
 #X connect 24 0 17 0;
-#X connect 26 0 24 0;
-#X connect 26 1 23 0;
+#X connect 26 0 29 0;
+#X connect 26 1 28 0;
+#X connect 28 0 23 0;
+#X connect 29 0 24 0;
diff --git a/doc/mtx_rifft-help.pd b/doc/mtx_rifft-help.pd
index 94ea8ed..be8f908 100644
--- a/doc/mtx_rifft-help.pd
+++ b/doc/mtx_rifft-help.pd
@@ -20,12 +20,13 @@ the result \;;
 #X obj 49 307 mtx_cos;
 #X obj 104 307 mtx_sin;
 #X text 217 232 <-- scroll here to select delay;
-#X obj 49 381 mtx_print rifft;
-#X obj 49 359 mtx_int;
+#X obj 49 411 mtx_print rifft;
+#X obj 49 389 mtx_int;
 #X obj 103 255 * -3.14159;
 #X text 89 15 [mtx_rifft];
 #X obj 262 463 mtx_rfft;
 #X obj 49 336 mtx_rifft;
+#X obj 49 362 mtx_+ 0.5;
 #X connect 4 0 15 0;
 #X connect 5 0 7 0;
 #X connect 6 0 5 0;
@@ -40,4 +41,5 @@ the result \;;
 #X connect 17 0 24 1;
 #X connect 20 0 19 0;
 #X connect 21 0 8 1;
-#X connect 24 0 20 0;
+#X connect 24 0 25 0;
+#X connect 25 0 20 0;
diff --git a/src/iemmatrix_sources.c b/src/iemmatrix_sources.c
index 2369ab4..2d77635 100644
--- a/src/iemmatrix_sources.c
+++ b/src/iemmatrix_sources.c
@@ -37,11 +37,13 @@ void iemmatrix_sources_setup(void)
 	iemtx_eq_setup(); /* mtx_eq.c */
 	iemtx_exp_setup(); /* mtx_exp.c */
 	iemtx_eye_setup(); /* mtx_eye.c */
+	iemtx_fft_setup(); /* mtx_fft.c */
 	iemtx_fill_setup(); /* mtx_fill.c */
 	iemtx_find_setup(); /* mtx_find.c */
 	iemtx_gauss_setup(); /* mtx_gauss.c */
 	iemtx_ge_setup(); /* mtx_ge.c */
 	iemtx_gt_setup(); /* mtx_gt.c */
+	iemtx_ifft_setup(); /* mtx_ifft.c */
 	iemtx_index_setup(); /* mtx_index.c */
 	iemtx_int_setup(); /* mtx_int.c */
 	iemtx_inverse_setup(); /* mtx_inverse.c */
@@ -53,8 +55,8 @@ void iemmatrix_sources_setup(void)
 	iemtx_mean_setup(); /* mtx_mean.c */
 	iemtx_min2_setup(); /* mtx_min2.c */
 	iemtx_minmax_setup(); /* mtx_minmax.c */
-	iemtx_mul__setup(); /* mtx_mul~.c */
 	iemtx_mul_setup(); /* mtx_mul.c */
+	iemtx_mul__setup(); /* mtx_mul~.c */
 	iemtx_neq_setup(); /* mtx_neq.c */
 	iemtx_not_setup(); /* mtx_not.c */
 	iemtx_ones_setup(); /* mtx_ones.c */
@@ -68,11 +70,11 @@ void iemmatrix_sources_setup(void)
 	iemtx_repmat_setup(); /* mtx_repmat.c */
 	iemtx_resize_setup(); /* mtx_resize.c */
 	iemtx_reverse_setup(); /* mtx_reverse.c */
+	iemtx_rfft_setup(); /* mtx_rfft.c */
+	iemtx_rifft_setup(); /* mtx_rifft.c */
 	iemtx_rmstodb_setup(); /* mtx_rmstodb.c */
 	iemtx_roll_setup(); /* mtx_roll.c */
 	iemtx_row_setup(); /* mtx_row.c */
-	iemtx_rowrfft_setup(); /* mtx_rowrfft.c */
-	iemtx_rowrifft_setup(); /* mtx_rowrifft.c */
 	iemtx_scroll_setup(); /* mtx_scroll.c */
 	iemtx_sin_setup(); /* mtx_sin.c */
 	iemtx_size_setup(); /* mtx_size.c */
diff --git a/src/iemmatrix_sources.h b/src/iemmatrix_sources.h
index 61738c5..fcd09d8 100644
--- a/src/iemmatrix_sources.h
+++ b/src/iemmatrix_sources.h
@@ -35,11 +35,13 @@ void iemtx_element_setup(void); /* mtx_element.c */
 void iemtx_eq_setup(void); /* mtx_eq.c */
 void iemtx_exp_setup(void); /* mtx_exp.c */
 void iemtx_eye_setup(void); /* mtx_eye.c */
+void iemtx_fft_setup(void); /* mtx_fft.c */
 void iemtx_fill_setup(void); /* mtx_fill.c */
 void iemtx_find_setup(void); /* mtx_find.c */
 void iemtx_gauss_setup(void); /* mtx_gauss.c */
 void iemtx_ge_setup(void); /* mtx_ge.c */
 void iemtx_gt_setup(void); /* mtx_gt.c */
+void iemtx_ifft_setup(void); /* mtx_ifft.c */
 void iemtx_index_setup(void); /* mtx_index.c */
 void iemtx_int_setup(void); /* mtx_int.c */
 void iemtx_inverse_setup(void); /* mtx_inverse.c */
@@ -51,8 +53,8 @@ void iemtx_max2_setup(void); /* mtx_max2.c */
 void iemtx_mean_setup(void); /* mtx_mean.c */
 void iemtx_min2_setup(void); /* mtx_min2.c */
 void iemtx_minmax_setup(void); /* mtx_minmax.c */
-void iemtx_mul__setup(void); /* mtx_mul~.c */
 void iemtx_mul_setup(void); /* mtx_mul.c */
+void iemtx_mul__setup(void); /* mtx_mul~.c */
 void iemtx_neq_setup(void); /* mtx_neq.c */
 void iemtx_not_setup(void); /* mtx_not.c */
 void iemtx_ones_setup(void); /* mtx_ones.c */
@@ -66,11 +68,11 @@ void iemtx_rand_setup(void); /* mtx_rand.c */
 void iemtx_repmat_setup(void); /* mtx_repmat.c */
 void iemtx_resize_setup(void); /* mtx_resize.c */
 void iemtx_reverse_setup(void); /* mtx_reverse.c */
+void iemtx_rfft_setup(void); /* mtx_rfft.c */
+void iemtx_rifft_setup(void); /* mtx_rifft.c */
 void iemtx_rmstodb_setup(void); /* mtx_rmstodb.c */
 void iemtx_roll_setup(void); /* mtx_roll.c */
 void iemtx_row_setup(void); /* mtx_row.c */
-void iemtx_rowrfft_setup(void); /* mtx_rowrfft.c */
-void iemtx_rowrifft_setup(void); /* mtx_rowrifft.c */
 void iemtx_scroll_setup(void); /* mtx_scroll.c */
 void iemtx_sin_setup(void); /* mtx_sin.c */
 void iemtx_size_setup(void); /* mtx_size.c */
diff --git a/src/mtx_rfft.c b/src/mtx_rfft.c
index e97f338..bad5ef2 100644
--- a/src/mtx_rfft.c
+++ b/src/mtx_rfft.c
@@ -15,17 +15,32 @@
 #include "iemmatrix.h"
 #include <stdlib.h>
 
+#ifdef DHAVE_FFTW3
+#include <fftw3.h>
+#endif
+
 static t_class *mtx_rfft_class;
 
+#ifdef DHAVE_FFTW3
+enum ComplexPart { REALPART=0,  IMAGPART=1};
+#endif
+
 typedef struct _MTXRfft_ MTXRfft;
 struct _MTXRfft_
 {
   t_object x_obj;
   int size;
   int size2;
-
+#ifdef DHAVE_FFTW3
+  int fftn;
+  int rows;
+  fftw_plan *fftplan;
+  fftw_complex *f_out;
+  double *f_in;
+#else
   t_float *f_re;
   t_float *f_im;
+#endif
 
   t_outlet *list_re_out;
   t_outlet *list_im_out;
@@ -36,10 +51,23 @@ struct _MTXRfft_
 
 static void deleteMTXRfft (MTXRfft *x) 
 {
+#if DHAVE_FFTW3
+  int n;
+  if (x->fftplan) {
+     for (n=0; n<x->rows; n++) 
+	fftw_destroy_plan(x->fftplan[n]);
+     free(x->fftplan);
+  }
+  if (x->f_out)
+     free(x->f_out);
+  if (x->f_in)
+     free(x->f_in);
+#else
   if (x->f_re)
      free (x->f_re);
   if (x->f_im) 
      free (x->f_im);
+#endif
   if (x->list_re)
      free (x->list_re);
   if (x->list_im)
@@ -51,9 +79,16 @@ static void *newMTXRfft (t_symbol *s, int argc, t_atom *argv)
   MTXRfft *x = (MTXRfft *) pd_new (mtx_rfft_class);
   x->list_re_out = outlet_new (&x->x_obj, gensym("matrix"));
   x->list_im_out = outlet_new (&x->x_obj, gensym("matrix"));
-
   x->size=x->size2=0;
+#ifdef DHAVE_FFTW3
+  x->fftn=0;
+  x->rows=0;
+  x->f_in=0;
+  x->f_out=0;
+  x->fftplan=0;
+#else
   x->f_re=x->f_im=0;
+#endif
   x->list_re=x->list_im=0;
   
   return ((void *) x);
@@ -94,12 +129,26 @@ static void writeFloatIntoList (int n, t_atom *l, t_float *f)
   for (;n--;f++, l++) 
     SETFLOAT (l, *f);
 }
+
 static void readFloatFromList (int n, t_atom *l, t_float *f) 
 {
   while (n--) 
     *f++ = atom_getfloat (l++);
 }
 
+#ifdef DHAVE_FFTW3
+static void writeFFTWComplexPartIntoList (int n, t_atom *l, fftw_complex *f, enum ComplexPart p) 
+{
+  for (;n--;f++, l++) 
+    SETFLOAT (l, ((t_float)*f[p]));
+}
+static void readDoubleFromList (int n, t_atom *l, double *f) 
+{
+  while (n--) 
+    *f++ = (double)atom_getfloat (l++);
+}
+#endif
+
 static void mTXRfftMatrix (MTXRfft *x, t_symbol *s, 
 			      int argc, t_atom *argv)
 {
@@ -112,8 +161,13 @@ static void mTXRfftMatrix (MTXRfft *x, t_symbol *s,
   int fft_count;
   t_atom *list_re = x->list_re;
   t_atom *list_im = x->list_im;
+#ifdef DHAVE_FFTW3
+  fftw_complex *f_out = x->f_out;
+  double *f_in = x->f_in;
+#else
   t_float *f_re = x->f_re;
   t_float *f_im = x->f_im;
+#endif
 
   /* fftsize check */
   if (!size)
@@ -127,8 +181,30 @@ static void mTXRfftMatrix (MTXRfft *x, t_symbol *s,
     /* ok, do the FFT! */
 
     /* memory things */
+#ifdef DHAVE_FFTW3
+    if ((x->rows!=rows)||(columns!=x->fftn)){
+       f_out=(fftw_complex*)realloc(f_out, sizeof(fftw_complex)*(size2-2));
+       f_in=(double*)realloc(f_in, sizeof(double)*size);
+       x->f_in = f_in;
+       x->f_out = f_out;
+       for (fft_count=0; fft_count<x->rows; fft_count++) {
+          fftw_destroy_plan(x->fftplan[fft_count]);
+       }
+       x->fftplan = (fftw_plan*)realloc(x->fftplan, sizeof(fftplan)*rows);
+       for (fft_count=0; fft_count<rows; fft_count++, f_in+=columns, f_out+=columns_re) {
+	  x->fftplan[fft_count] = fftw_plan_dft_r2c_1d (columns,f_in,f_out,FFTW_ESTIMATE);
+       }
+       x->fftn=columns;
+       x->rows=rows;
+       f_in=x->f_in;
+       f_out=x->f_out;
+    }
+#else
     f_re=(t_float*)realloc(f_re, sizeof(t_float)*size);
     f_im=(t_float*)realloc(f_im, sizeof(t_float)*size);
+    x->f_re = f_re;
+    x->f_im = f_im;
+#endif
     list_re=(t_atom*)realloc(list_re, sizeof(t_atom)*size2);
     list_im=(t_atom*)realloc(list_im, sizeof(t_atom)*size2);
 
@@ -136,22 +212,30 @@ static void mTXRfftMatrix (MTXRfft *x, t_symbol *s,
     x->size2 = size2;
     x->list_im = list_im;
     x->list_re = list_re;
-    x->f_re = f_re;
-    x->f_im = f_im;
 
     /* main part */
+#ifdef DHAVE_FFTW3
+    readDoubleFromList (size, argv, f_in);
+#else
     readFloatFromList (size, argv, f_re);
+#endif
 
-    fft_count = rows;
     list_re += 2;
     list_im += 2;
-    while (fft_count--){ 
+    for (fft_count=0;fft_count<rows;fft_count++){ 
+#if DHAVE_FFTW3
+      fftw_execute(x->fftplan[fft_count]);
+      writeFFTWComplexPartIntoList(columns_re,list_re,f_out,REALPART);
+      writeFFTWComplexPartIntoList(columns_re,list_im,f_out,IMAGPART);
+      f_out+=columns_re;
+#else
       mayer_realfft (columns, f_re);
       fftRestoreImag (columns, f_re, f_im);
       writeFloatIntoList (columns_re, list_re, f_re);
       writeFloatIntoList (columns_re, list_im, f_im);
       f_im += columns;
       f_re += columns;
+#endif
       list_re += columns_re;
       list_im += columns_re;
     }
diff --git a/src/mtx_rifft.c b/src/mtx_rifft.c
index 7a1142f..d10fc28 100644
--- a/src/mtx_rifft.c
+++ b/src/mtx_rifft.c
@@ -15,8 +15,16 @@
 #include "iemmatrix.h"
 #include <stdlib.h>
 
+#ifdef DHAVE_FFTW3
+#include <fftw3.h>
+#endif
+
 static t_class *mtx_rifft_class;
 
+#ifdef DHAVE_FFTW3
+enum ComplexPart { REALPART=0,  IMAGPART=1};
+#endif
+
 typedef struct _MTXRifft_
 {
   t_object x_obj;
@@ -25,10 +33,16 @@ typedef struct _MTXRifft_
   int columns_re;
   int size;
   int size2;
+#ifdef DHAVE_FFTW3  
+  fftw_plan *fftplan;
+  fftw_complex *f_in;
+  double *f_out;
+#else
   t_float renorm_fac;
 
   t_float *f_re;
   t_float *f_im;
+#endif
 
   t_outlet *list_re_out;
   t_outlet *list_im_out;
@@ -77,6 +91,18 @@ static void ifftPrepareReal (int n, t_float *re, t_float *im)
     *++re = -*--im;
 }
 
+#ifdef DHAVE_FFTW3
+static void readFFTWComplexPartFromList (int n, t_atom *l, fftw_complex *f, enum ComplexPart p) 
+{
+  for (;n--;f++, l++) 
+    *f[p] = (double) atom_getfloat (l);
+}
+static void writeDoubleIntoList (int n, t_atom *l, double *f) 
+{
+  while (n--) 
+    SETFLOAT (l++,((t_float)(*f++)));
+}
+#endif
 
 static void *newMTXRifft (t_symbol *s, int argc, t_atom *argv)
 {
@@ -98,8 +124,13 @@ static void mTXRifftMatrixCold (MTXRifft *x, t_symbol *s,
   int size = rows * columns;
   int ifft_count;
   t_atom *list_re = x->list_re;
+#ifdef DHAVE_FFTW3
+  fftw_complex *f_in = x->f_in;
+  double *f_out = x->f_out;
+#else
   t_float *f_re = x->f_re;
   t_float *f_im = x->f_im;
+#endif
 
   /* ifftsize check */
   if (columns_re < 3)
@@ -113,26 +144,53 @@ static void mTXRifftMatrixCold (MTXRifft *x, t_symbol *s,
   else if (columns == (1 << ilog2(columns))) {
 
     /* memory things */
+#ifdef DHAVE_FFTW3
+    if ((x->rows!=rows)||(columns!=x->columns)){
+       for (ifft_count=0;ifft_count<rows;ifft_count++) {
+	  fftw_destroy_plan(x->fftplan[ifft_count]);
+       }
+       x->fftplan=(fftw_plan*)realloc(x->fftplan,sizeof(fftw_plan)*rows);
+       f_in=(fftw_complex*)realloc(f_in,sizeof(fftw_complex)*size2);
+       f_out=(double*)realloc(f_out,sizeof(double)*size);
+       x->f_out = f_out;
+       x->f_in = f_in;
+       for (ifft_count=0;ifft_count<rows;ifft_count++) {
+	  x->fftplan[ifft_count]=fftw_plan_dft_c2r_1d(columns,f_in,f_out,FFTW_ESTIMATE);
+	  f_out+=columns;
+	  f_in+=columns_re;
+       }
+       f_in=x->f_in;
+       f_out=x->f_out;
+ }
+#else
     f_re=(t_float*)realloc(f_re, sizeof(t_float)*size);
     f_im=(t_float*)realloc(f_im, sizeof(t_float)*size);
-    list_re=(t_atom*)realloc(list_re, sizeof(t_atom)*(size+2));
+    x->f_re = f_re;
+    x->f_im = f_im;
+#endif
 
+    list_re=(t_atom*)realloc(list_re, sizeof(t_atom)*(size+2));
     x->size = size;
     x->size2 = size2;
     x->rows = rows;
     x->columns = columns;
     x->columns_re = columns_re;
     x->list_re = list_re;
-    x->f_re = f_re;
-    x->f_im = f_im;
       
     /* main part: reading imaginary part */
     ifft_count = rows;
+#ifndef DHAVE_FFTW3
     x->renorm_fac = 1.0f / columns;
+#endif
     while (ifft_count--) {
+#ifdef DHAVE_FFTW3
+      readFFTWComplexPartFromList(columns_re, argv, f_in, IMAGPART);
+      f_in += columns_re;
+#else
       readFloatFromList (columns_re, argv, f_im);
-      argv += columns_re;
       f_im += columns;
+#endif
+      argv += columns_re;
     }
     /* do nothing else! */
   }
@@ -150,10 +208,14 @@ static void mTXRifftMatrixHot (MTXRifft *x, t_symbol *s,
   int in_size = argc-2;
   int size2 = x->size2;
   int ifft_count;
-  t_atom *ptr_re = x->list_re;
+#ifdef DHAVE_FFTW3
+  fftw_complex *f_in = x->f_in;
+  double *f_out = x->f_out;
+#else
   t_float *f_re = x->f_re;
   t_float *f_im = x->f_im;
   t_float renorm_fac = x->renorm_fac;
+#endif
 
   /* ifftsize check */
   if ((rows != x->rows) || 
@@ -164,27 +226,35 @@ static void mTXRifftMatrixHot (MTXRifft *x, t_symbol *s,
   else if (!x->size2)
     post("mtx_rifft: invalid right side matrix");
   else { /* main part */
-    ifft_count = rows;
-    ptr_re += 2;
-    while (ifft_count--){ 
+    for (ifft_count=0;ifft_count<rows;ifft_count++){ 
+#ifdef DHAVE_FFTW3
+      readFFTWComplexPartFromList(columns_re,argv,f_in,REALPART);
+      fftw_execute(x->fftw_execute[ifft_count]);
+      f_in+=columns_re;
+#else
       readFloatFromList (columns_re, argv, f_re);
       ifftPrepareReal (columns, f_re, f_im);
       mayer_realifft (columns, f_re);
       multiplyVector (columns, f_re, renorm_fac);
       f_im += columns;
       f_re += columns;
-      ptr_re += columns;
+#endif
       argv += columns_re;
     }
-    ptr_re = x->list_re;
+#ifndef DHAVE_FFTW3
     f_re = x->f_re;
+#endif
     size2 = x->size2;
 
-    SETSYMBOL(ptr_re, gensym("matrix"));
-    SETFLOAT(ptr_re, rows);
-    SETFLOAT(&ptr_re[1], x->columns);
-    writeFloatIntoList (size, ptr_re+2, f_re);
-    outlet_anything(x->list_re_out, gensym("matrix"), size+2, ptr_re);
+    SETSYMBOL(x->list_re, gensym("matrix"));
+    SETFLOAT(x->list_re, rows);
+    SETFLOAT(x->list_re+1, x->columns);
+#ifdef DHAVE_FFTW3
+    writeDoubleIntoList (size, x->list_re+2, f_out);
+#else
+    writeFloatIntoList (size, x->list_re+2, f_re);
+#endif
+    outlet_anything(x->list_re_out, gensym("matrix"), size+2, x->list_re);
   }
 }
 
@@ -198,10 +268,23 @@ static void mTXRifftBang (MTXRifft *x)
 
 static void deleteMTXRifft (MTXRifft *x) 
 {
+#ifdef DHAVE_FFTW3
+  int n;
+  if (x->fftplan) {
+     for (n=0; n<x->rows; n++) 
+	fftw_destroy_plan(x->fftplan[n]);
+     free(x->fftplan);
+  }
+  if (x->f_out)
+     free(x->f_out);
+  if (x->f_in)
+     free(x->f_in);
+#else
   if (x->f_re)
      free(x->f_re);
   if (x->f_im)
      free(x->f_im);
+#endif
   if (x->list_re)
      free(x->list_re);
   if (x->list_im)