[PD-cvs] pd/src d_fft.c,1.1.1.1,1.1.1.1.16.1

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 ! }