Update of /cvsroot/pure-data/pd/src
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv29916
Modified Files:
Tag: devel_0_37
d_fft.c
Log Message:
added support for fftw
Index: d_fft.c
===================================================================
RCS file: /cvsroot/pure-data/pd/src/d_fft.c,v
retrieving revision 1.1.1.1
retrieving revision 1.1.1.1.16.1
diff -C2 -d -r1.1.1.1 -r1.1.1.1.16.1
*** d_fft.c 29 Jul 2002 17:05:56 -0000 1.1.1.1
--- d_fft.c 2 Mar 2004 11:04:50 -0000 1.1.1.1.16.1
***************
*** 5,8 ****
--- 5,14 ----
#include "m_pd.h"
+ #ifdef USE_FFTW
+ #include "fftw3.h"
+ #include <string.h>
+ #endif
+
+ #ifndef USE_FFTW
/* ------------------------ fft~ and ifft~ -------------------------------- */
static t_class *sigfft_class, *sigifft_class;
***************
*** 41,45 ****
for (;n--; in1++, in2++)
{
! float f = *in1;
*in1 = *in2;
*in2 = f;
--- 47,51 ----
for (;n--; in1++, in2++)
{
! float f = *in1;
*in1 = *in2;
*in2 = f;
***************
*** 141,144 ****
--- 147,151 ----
}
+
static t_int *sigrfft_perform(t_int *w)
{
***************
*** 185,190 ****
--- 192,200 ----
}
+
/* ----------------------- rifft~ -------------------------------- */
+
+
static t_class *sigrifft_class;
***************
*** 245,248 ****
--- 255,702 ----
}
+ #else
+
+ /* Support for fftw3 by Tim Blechmann */
+
+ /* ------------------------ fft~ and ifft~ -------------------------------- */
+
+ static t_class *sigfftw_class, *sigifftw_class;
+
+ typedef struct fftw
+ {
+ t_object x_obj;
+ float x_f;
+
+ fftwf_plan plan;
+ t_int bins;
+ fftwf_complex * incomplex;
+ fftwf_complex * outcomplex;
+
+ } t_sigfftw;
+
+ static void *sigfftw_new(void)
+ {
+ t_sigfftw *x = (t_sigfftw *)pd_new(sigfftw_class);
+ outlet_new(&x->x_obj, gensym("signal"));
+ outlet_new(&x->x_obj, gensym("signal"));
+ inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_signal, &s_signal);
+ x->x_f = 0;
+
+ //get ready for the default blocksize
+ x->bins=64;
+ x->incomplex = fftwf_malloc(sizeof(fftw_complex) * x->bins);
+ x->outcomplex = fftwf_malloc(sizeof(fftw_complex) * x->bins);
+ x->plan = fftwf_plan_dft_1d(x->bins,x->incomplex,x->outcomplex,
+ FFTW_FORWARD,FFTW_MEASURE);
+
+ return (x);
+ }
+
+ static void *sigifftw_new(void)
+ {
+ t_sigfftw *x = (t_sigfftw *)pd_new(sigifftw_class);
+ outlet_new(&x->x_obj, gensym("signal"));
+ outlet_new(&x->x_obj, gensym("signal"));
+ inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_signal, &s_signal);
+ x->x_f = 0;
+
+ //get ready for the default blocksize
+ x->bins=64;
+ x->incomplex = fftwf_malloc(sizeof(fftw_complex) * x->bins);
+ x->outcomplex = fftwf_malloc(sizeof(fftw_complex) * x->bins);
+ x->plan = fftwf_plan_dft_1d(x->bins,x->incomplex,x->outcomplex,
+ FFTW_BACKWARD,FFTW_MEASURE);
+ return (x);
+ }
+
+ static void sigfftw_free(t_sigfftw * x)
+ {
+ fftwf_free(x->incomplex);
+ fftwf_free(x->outcomplex);
+ fftwf_destroy_plan(x->plan);
+
+ }
+
+ static void sigifftw_free(t_sigfftw * x)
+ {
+ fftwf_free(x->incomplex);
+ fftwf_free(x->outcomplex);
+ fftwf_destroy_plan(x->plan);
+ }
+
+ static t_int *sigfftw_perform(t_int *w)
+ {
+ t_sigfftw * x =(t_sigfftw *) w[1];
+ t_sample * in1 =(t_sample *) w[2];
+ t_sample * in2 =(t_sample *) w[3];
+ t_sample * out1 =(t_sample *) w[4];
+ t_sample * out2 =(t_sample *) w[5];
+ t_int n =(t_int) w[6];
+
+ if (n != x->bins)
+ {
+ x->bins=n;
+
+ //re-allocate fft buffers
+ fftwf_free(x->incomplex);
+ x->incomplex = fftwf_malloc(sizeof(fftwf_complex) * x->bins);
+ fftwf_free(x->outcomplex);
+ x->outcomplex = fftwf_malloc(sizeof(fftwf_complex) * x->bins);
+
+ //set plan, this might take a few seconds
+ //but you don't have to do that on the fly...
+ fftwf_destroy_plan(x->plan);
+ x->plan = fftwf_plan_dft_1d(x->bins,x->incomplex,x->outcomplex,
+ FFTW_FORWARD,FFTW_MEASURE);
+ }
+
+ t_int i=n;
+ while(i!=0)
+ {
+ --i;
+ x->incomplex[i][0]=*(in1+i);
+ x->incomplex[i][1]=*(in2+i);
+ }
+ fftwf_execute(x->plan);
+
+ i=n;
+ while(i!=0)
+ {
+ --i;
+ *(out1+i)=x->outcomplex[i][0];
+ *(out2+i)=x->outcomplex[i][1];
+ }
+
+ return (w+7);
+ }
+
+ static t_int *sigifftw_perform(t_int *w)
+ {
+ t_sigfftw * x =(t_sigfftw *) w[1];
+ t_sample * in1 =(t_sample *) w[2];
+ t_sample * in2 =(t_sample *) w[3];
+ t_sample * out1 =(t_sample *) w[4];
+ t_sample * out2 =(t_sample *) w[5];
+ t_int n =(t_int) w[6];
+
+ if (n != x->bins)
+ {
+ x->bins=n;
+
+ //re-allocate fft buffers
+ fftwf_free(x->incomplex);
+ x->incomplex = fftwf_malloc(sizeof(fftwf_complex) * x->bins);
+ fftwf_free(x->outcomplex);
+ x->outcomplex = fftwf_malloc(sizeof(fftwf_complex) * x->bins);
+
+ //set plan, this might take a few seconds
+ //but you don't have to do that on the fly...
+ fftwf_destroy_plan(x->plan);
+ x->plan = fftwf_plan_dft_1d(x->bins,x->incomplex,x->outcomplex,
+ FFTW_BACKWARD,FFTW_MEASURE);
+ }
+
+ t_int i=n;
+ while(i!=0)
+ {
+ --i;
+ x->incomplex[i][0]=*(in1+i);
+ x->incomplex[i][1]=*(in2+i);
+ }
+ fftwf_execute(x->plan);
+
+ i=n;
+ while(i!=0)
+ {
+ --i;
+ *(out1+i)=x->outcomplex[i][0];
+ *(out2+i)=x->outcomplex[i][1];
+ }
+
+ return (w+7);
+
+ }
+
+ static void sigfftw_dspx(t_sigfftw *x, t_signal **sp, t_int *(*f)(t_int *w))
+ {
+ int n = sp[0]->s_n;
+ float *in1 = sp[0]->s_vec;
+ float *in2 = sp[1]->s_vec;
+ float *out1 = sp[2]->s_vec;
+ float *out2 = sp[3]->s_vec;
+ dsp_add(f, 6, x, in1, in2, out1, out2, n);
+ }
+
+ static void sigfftw_dsp(t_sigfftw *x, t_signal **sp)
+ {
+ sigfftw_dspx(x, sp, sigfftw_perform);
+ }
+
+ static void sigifftw_dsp(t_sigfftw *x, t_signal **sp)
+ {
+ sigfftw_dspx(x, sp, sigifftw_perform);
+ }
+
+ static void sigfftw_setup(void)
+ {
+ sigfftw_class = class_new(gensym("fft~"), sigfftw_new,
+ (t_method) sigfftw_free,
+ sizeof(t_sigfftw), 0, 0);
+ CLASS_MAINSIGNALIN(sigfftw_class, t_sigfftw, x_f);
+ class_addmethod(sigfftw_class, (t_method)sigfftw_dsp,
+ gensym("dsp"), 0);
+
+ sigifftw_class = class_new(gensym("ifft~"), sigifftw_new,
+ (t_method) sigifftw_free,
+ sizeof(t_sigfftw), 0, 0);
+ CLASS_MAINSIGNALIN(sigifftw_class, t_sigfftw, x_f);
+ class_addmethod(sigifftw_class, (t_method)sigifftw_dsp,
+ gensym("dsp"), 0);
+ class_sethelpsymbol(sigifftw_class, gensym("fft~"));
+ }
+
+
+ /* ----------------------- rfft~ --------------------------------- */
+
+ static t_class *sigrfftw_class;
+
+ typedef struct rfftw
+ {
+ t_object x_obj; //me
+
+ fftwf_plan plan; //fftw plan
+ t_int bins; //number of bins
+ t_sample * infft; //array fftw is working on
+ t_sample * outfft; //array fftw uses to output it's values
+
+ float x_f;
+
+ } t_sigrfftw;
+
+ //constructor
+ static void *sigrfftw_new(void)
+ {
+ t_sigrfftw *x = (t_sigrfftw *)pd_new(sigrfftw_class);
+ outlet_new(&x->x_obj, gensym("signal"));
+ outlet_new(&x->x_obj, gensym("signal"));
+ x->x_f = 0;
+
+ //get ready for the default blocksize
+ x->bins=64;
+ x->infft = fftwf_malloc(sizeof(t_sample) * x->bins);
+ x->outfft = fftwf_malloc(sizeof(t_sample) * x->bins);
+ //it seems that the halfcomplex transformation is a bit faster
+ x->plan = fftwf_plan_r2r_1d(x->bins,x->infft,x->outfft,
+ FFTW_FORWARD,FFTW_MEASURE);
+
+ return (x);
+ }
+
+ //destructor
+ static void sigrfftw_free(t_sigrfftw *x)
+ {
+ fftwf_free(x->infft);
+ fftwf_free(x->outfft);
+ fftwf_destroy_plan(x->plan);
+ }
+
+ static t_int *sigrfftw_perform(t_int *w)
+ {
+ t_sigrfftw * x =(t_sigrfftw *) w[1];
+ t_sample * in1 =(t_sample *) w[2];
+ t_sample * out1 =(t_sample *) w[3];
+ t_sample * out2 =(t_sample *) w[4];
+ t_int n =(t_int) w[5];
+
+ if (n != x->bins)
+ {
+ x->bins=n;
+
+ //re-allocate fft buffers
+ fftwf_free(x->infft);
+ x->infft = fftwf_malloc(sizeof(float) * x->bins);
+ fftwf_free(x->outfft);
+ x->outfft = fftwf_malloc(sizeof(float) * x->bins);
+
+ //set plan, this might take a few seconds
+ //but you don't have to do that on the fly...
+ fftwf_destroy_plan(x->plan);
+ x->plan = fftwf_plan_r2r_1d(x->bins,x->infft,x->outfft,
+ FFTW_FORWARD,FFTW_MEASURE);
+ }
+
+ memcpy(x->infft,in1,n*sizeof(t_sample));
+ fftwf_execute(x->plan);
+ memcpy(out1,x->outfft,n/2*(sizeof(t_sample)));
+
+ t_int i = n/2;
+ while (i)
+ {
+ out2[i] = -(x->outfft[n-i]);
+ --i;
+ }
+ return (w+6);
+ }
+
+ static void sigrfftw_dsp(t_sigrfftw *x, t_signal **sp)
+ {
+ int n = sp[0]->s_n, n2 = (n>>1);
+ float *in1 = sp[0]->s_vec;
+ float *out1 = sp[1]->s_vec;
+ float *out2 = sp[2]->s_vec;
+
+ if (n < 4)
+ {
+ error("fft: minimum 4 points");
+ return;
+ }
+ else
+ dsp_add(sigrfftw_perform,5,x,in1,out1,out2,n);
+
+ dsp_add_zero(out1 + n2, n2);
+ dsp_add_zero(out2 + n2, n2);
+ }
+
+ static void sigrfftw_setup(void)
+ {
+ sigrfftw_class = class_new(gensym("rfft~"), sigrfftw_new,
+ (t_method)sigrfftw_free,
+ sizeof(t_sigrfftw), 0, 0);
+ CLASS_MAINSIGNALIN(sigrfftw_class, t_sigrfftw, x_f);
+ class_addmethod(sigrfftw_class, (t_method)sigrfftw_dsp,
+ gensym("dsp"), 0);
+ class_sethelpsymbol(sigrfftw_class, gensym("fft~"));
+ }
+
+
+ /* ----------------------- rifft~ -------------------------------- */
+
+ static t_class *sigrifftw_class;
+
+ typedef struct rifftw
+ {
+ t_object x_obj; //me
+
+ fftwf_iodim iodim;
+
+ fftwf_plan plan; //fftw plan
+ t_int bins; //number of bins
+ t_sample * infft;
+ t_sample * inimag;
+ t_sample * outfft;
+
+ float x_f;
+
+ } t_sigrifftw;
+
+ //constructor
+ static void *sigrifftw_new(void)
+ {
+ t_sigrifftw *x = (t_sigrifftw *)pd_new(sigrifftw_class);
+ inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_signal, &s_signal);
+ outlet_new(&x->x_obj, gensym("signal"));
+ x->x_f = 0;
+
+ //get ready for the default blocksize
+ x->bins=64;
+ x->infft = fftwf_malloc(sizeof(t_sample) * (x->bins));
+ //x->inimag = fftwf_malloc(sizeof(t_sample) * (x->bins)/2);
+ x->outfft = fftwf_malloc(sizeof(t_sample) * x->bins);
+
+ // x->iodim->1;
+
+ x->plan = fftwf_plan_r2r_1d(x->bins,x->infft,x->outfft,
+ FFTW_BACKWARD,FFTW_MEASURE);
+
+
+ return (x);
+ }
+
+ //destructor
+ static void sigrifftw_free(t_sigrifftw *x)
+ {
+ fftwf_free(x->infft);
+ fftwf_free(x->inimag);
+ fftwf_free(x->outfft);
+ fftwf_destroy_plan(x->plan);
+ }
+
+ static t_int *sigrifftw_perform(t_int *w)
+ {
+ t_sigrifftw * x =(t_sigrifftw *) w[1];
+ t_sample * in1 =(t_sample *) w[2];
+ t_sample * in2 =(t_sample *) w[3];
+ t_sample * out1 =(t_sample *) w[4];
+ t_int n =(t_int) w[5];
+
+ if (n != x->bins)
+ {
+ x->bins=n;
+
+ //re-allocate fft buffers
+ fftwf_free(x->infft);
+ x->infft = fftwf_malloc(sizeof(float) * (x->bins));
+ //fftwf_free(x->inimag);
+ //x->inimag = fftwf_malloc(sizeof(float) * (x->bins)/2);
+ fftwf_free(x->outfft);
+ x->outfft = fftwf_malloc(sizeof(float) * x->bins);
+
+ //set plan, this might take a few seconds
+ //but you don't have to do that on the fly...
+ fftwf_destroy_plan(x->plan);
+ x->plan = fftwf_plan_r2r_1d(x->bins,x->infft,x->outfft,
+ FFTW_BACKWARD,FFTW_MEASURE);
+
+ //x->plan = fftwf_plan_guru_split_dft_c2r(1,iodim ,x->bins,x->infft,x->outfft,FFTW_BACKWARD,FFTW_MEASURE);
+ }
+
+ memcpy(x->infft,in1,n/2*sizeof(t_sample));
+
+ t_int i = n/2;
+ while (i!=n)
+ {
+ x->infft[i]=-in2[n-i];
+ ++i;
+ }
+
+ fftwf_execute(x->plan);
+ memcpy(out1,x->outfft,n*(sizeof(t_sample)));
+
+ return (w+6);
+ }
+
+ static void sigrifftw_dsp(t_sigrifftw *x, t_signal **sp)
+ {
+ int n = sp[0]->s_n, n2 = (n>>1);
+ float *in1 = sp[0]->s_vec;
+ float *in2 = sp[1]->s_vec;
+ float *out1 = sp[2]->s_vec;
+
+ if (n < 4)
+ {
+ error("fft: minimum 4 points");
+ return;
+ }
+
+ else
+ dsp_add(sigrifftw_perform,5,x,in1,in2,out1,n);
+
+ }
+
+ static void sigrifftw_setup(void)
+ {
+ sigrifftw_class = class_new(gensym("rifft~"), sigrifftw_new,
+ (t_method)sigrifftw_free,
+ sizeof(t_sigrifftw), 0, 0);
+ CLASS_MAINSIGNALIN(sigrifftw_class, t_sigrifftw, x_f);
+ class_addmethod(sigrifftw_class, (t_method)sigrifftw_dsp,
+ gensym("dsp"), 0);
+ class_sethelpsymbol(sigrifftw_class, gensym("fft~"));
+ }
+
+ #endif
+ /* end of FFTW support */
+
+
/* ----------------------- framp~ -------------------------------- */
***************
*** 336,342 ****
void d_fft_setup(void)
{
sigfft_setup();
sigrfft_setup();
sigrifft_setup();
! sigframp_setup();
}
--- 790,804 ----
void d_fft_setup(void)
{
+ sigframp_setup();
+
+ #ifndef USE_FFTW
sigfft_setup();
sigrfft_setup();
sigrifft_setup();
! #else
! sigfftw_setup(); /* added by Tim Blechmann to support fftw */
! sigrifftw_setup();
! sigrfftw_setup();
! #endif
!
}
