diff options
Diffstat (limited to 'src/mtx_rifft.c')
| -rw-r--r-- | src/mtx_rifft.c | 113 |
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) |