aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorFranz Zotter <fzotter@users.sourceforge.net>2008-07-22 13:20:18 +0000
committerFranz Zotter <fzotter@users.sourceforge.net>2008-07-22 13:20:18 +0000
commit1f9a0c772d2178f0df6c2a6ea16e3c0986abeed4 (patch)
tree40a6ff6061f8d6aaaf6d6cabd36e808cd59c2fe0
parent626388b61dd670db7c6390334922a6c18a5a458a (diff)
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
-rw-r--r--doc/mtx_rfft-help.pd16
-rw-r--r--doc/mtx_rifft-help.pd8
-rw-r--r--src/iemmatrix_sources.c8
-rw-r--r--src/iemmatrix_sources.h8
-rw-r--r--src/mtx_rfft.c96
-rw-r--r--src/mtx_rifft.c113
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)