diff options
| -rw-r--r-- | doc/mtx_rfft-help.pd | 16 | ||||
| -rw-r--r-- | doc/mtx_rifft-help.pd | 8 | ||||
| -rw-r--r-- | src/iemmatrix_sources.c | 8 | ||||
| -rw-r--r-- | src/iemmatrix_sources.h | 8 | ||||
| -rw-r--r-- | src/mtx_rfft.c | 96 | ||||
| -rw-r--r-- | src/mtx_rifft.c | 113 |
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) |