aboutsummaryrefslogtreecommitdiff
path: root/src/mtx_rifft.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/mtx_rifft.c')
-rw-r--r--src/mtx_rifft.c113
1 files changed, 98 insertions, 15 deletions
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)