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

timblech at users.sourceforge.net timblech at users.sourceforge.net
Tue Mar 2 12:04:52 CET 2004


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





More information about the Pd-cvs mailing list