New checkin after sf upgrade

svn path=/trunk/externals/postlude/; revision=5152

5 files changed, 147 insertions, 64 deletions

diff --git a/flib/src/cc~.c b/flib/src/cc~.c
index 2d06d9b..da2d50b 100644
--- a/flib/src/cc~.c
+++ b/flib/src/cc~.c
@@ -17,29 +17,36 @@
    */
 
 
-/*Calculate the cross correlation of two signal vectors*/
+/*Calculate the cc correlation of two signal vectors*/
 
 /*The time domain implementation is based on code by Phil Bourke 
  * the frequency domain version is based on code by Charles Henry 
  * 
- * Specify a time delay as an argument for the time domain implemenation, for example an argument of 32 will give the correlation coefficients  for delays from -32 to 32 samples between the two input vectors
+ * Specify a time delay as an argument for the time domain implemenation, for example an argument of 32 will give the 
+ * correlation coefficients  for delays from -32 to 32 samples between the two input vectors
  *
- * Specify an argument of 'f' for the frequency domain implementation*/
+ * Specify an argument of 'f' for the frequency domain implementation,   
+ * or 'r' for the running cross covariance (not normalized)              instead of the numerical delay argument
+ * these two methods both have got positive delays on 0,N/2-1 and the negative delays (-N/2, -1) are indexed on N/2,N-1 */
 
 
 #include "flib.h"
 #define SQ(a) (a * a)
 
-static t_class *cross_class;
+static t_class *cc_class;
 
-typedef struct _cross {
+typedef struct _cc {
     t_object  x_obj;
     t_float f;
     t_int delay;
     t_int is_freq_domain;
-} t_cross;
+    t_float *buffer2Nsig1, *buffer2Nsig2;
+    t_float *output_prev_block;
+    t_int is_new;
+    t_int n;
+} t_cc;
 
-static t_int *cross_perform_time_domain(t_int *w)
+static t_int *cc_perform_time_domain(t_int *w)
 {
     t_sample *x = (t_sample *)(w[1]);
     t_sample *y = (t_sample *)(w[2]);
@@ -51,7 +58,7 @@ static t_int *cross_perform_time_domain(t_int *w)
   
    if(maxdelay > N * .5){
 	maxdelay = N * .5;
-	post("cross~: invalid maxdelay, must be <= blocksize/2");
+	post("cc~: invalid maxdelay, must be <= blocksize/2");
    }
  
 /* Calculate the mean of the two series x[], y[] */
@@ -97,103 +104,179 @@ static t_int *cross_perform_time_domain(t_int *w)
 }
 
 
-t_int *cross_perform_freq_domain(t_int *w)
+static t_int *cc_perform_freq_domain(t_int *w)
 {
-  t_cross *x = (t_cross *)(w[1]);
+  t_cc *x = (t_cc *)(w[1]);
 
   t_sample *sig1 = (t_sample *)(w[2]);
   t_sample *sig2 = (t_sample *)(w[3]);
   t_sample *out  = (t_sample *)(w[4]);
-  long int size  = (long int) w[5];
-  long int k     = size/2;
-  float *expsig1 = NULL;
-  float *revsig2 = NULL;
-  float temp, temp2;
-  long int i=0;
-  int well_defined=1;
-  int qtr, thrqtr;
-
-// The two signals are created, nonzero on 0 to N/4 and 3N/4 to N
-// This will be revised
-
-  expsig1=(float *) alloca(size*sizeof(float));
-  revsig2=(float *) alloca(size*sizeof(float));
-  qtr = size/4;
-  thrqtr = 3*size/4;
-  for (i=0; i < qtr ; i++)
+  t_int size  = (int) w[5];
+  x->n = size;
+  t_int size2 = size*2;
+  t_int half = size/2;
+  t_float *expsig1 = NULL;
+  t_float *revsig2 = NULL;
+  t_float temp, temp2;
+  t_int i=0;
+  t_int thrhalf;
+
+// This stuff here sets up two buffers to hold the previous N samples
+// To get the usual overlapping block (2) design on each input
+
+  if (x->is_new)
   {
-    expsig1[i]=sig1[i];
+    x->buffer2Nsig1=getbytes(size*sizeof(t_float));    
+    x->buffer2Nsig2=getbytes(size*sizeof(t_float));
+    x->output_prev_block=getbytes(size*sizeof(t_float));        
+    x->is_new=0;
+  }
+//  Here we set the buffers for the next round
+  for(i=half; i < size; i++)
+ {
+    x->buffer2Nsig1[i]=sig1[i];
+    x->buffer2Nsig2[i]=sig2[i];
+  }
+// The two signals are created, nonzero on 0 to 1/4 and 3/4 to 1
+// Using a block size of 2N, --size2
+
+  expsig1=(float *) getbytes(size2*sizeof(float));
+  revsig2=(float *) getbytes(size2*sizeof(float));
+
+// Loops for assignment of old values in buffer + new block
+  thrhalf = 3*half;
+  for (i=0; i < half ; i++)
+  {
+    expsig1[i]=x->buffer2Nsig1[i];
     revsig2[i]=0;
   }
-  for (i=qtr; i < thrqtr ; i++)
+  for (i=half; i < size ; i++)
+  {
+    expsig1[i]=x->buffer2Nsig1[i];
+    revsig2[i]=sig2[size-i];            /// Needs revision here, not too clear
+  }
+  expsig1[size]=sig1[0];
+  revsig2[size]=sig2[0];
+  for (i=size+1; i < thrhalf ; i++)
   {
-    expsig1[i]=sig1[i];
-    revsig2[i]=sig2[size-i];
+    expsig1[i]=sig1[i-size];
+    revsig2[i]=x->buffer2Nsig2[size2-i];
   }
-  for (i=thrqtr; i < size ; i++)
+  for (i=thrhalf; i < size2 ; i++)
   {
-    expsig1[i]=sig1[i];
+    expsig1[i]=sig1[i-size];
     revsig2[i]=0;
   }
 
-  mayer_realfft(size, expsig1);
-  mayer_realfft(size, revsig2);
+//  fft the two blocks and multiply them
+  mayer_realfft(size2, expsig1);
+  mayer_realfft(size2, revsig2);
+
   expsig1[0]*=revsig2[0];
-  expsig1[k]*=revsig2[k];
-  for(i=1; i < k; i++)
+  expsig1[size]*=revsig2[size];
+  for(i=1; i < size2; i++)  
   {
     temp=expsig1[i];
-    temp2=expsig1[size-i];
-    expsig1[i]=temp*revsig2[i]-temp2*revsig2[size-i];
-    expsig1[size-i]=temp*revsig2[size-i]+temp2*revsig2[i];
+    temp2=expsig1[size2-i];
+    expsig1[i]=temp*revsig2[i]-temp2*revsig2[size2-i];
+    expsig1[size2-i]=temp*revsig2[size2-i]+temp2*revsig2[i];
   }
   
-  mayer_realifft(size, expsig1);
-  for(i=0; i < size; i++)
+//  ifft
+  mayer_realifft(size2, expsig1);
+
+//  format the output:  this section formats the ouptut either as
+//  a simple cc or as a running cc
+if (x->is_freq_domain == 1)
+{
+  for(i=0; i < half; i++)
+  {
+    out[i]=expsig1[i]/size2;
+    out[half + i]=expsig1[half + i]/size2;
+  }
+} else {
+  for(i=0; i < half; i++)
   {
-    out[i]=expsig1[i];
+    out[i]=x->output_prev_block[i] + expsig1[i]/size2;
+    out[half + i]=x->output_prev_block[half + i] + expsig1[half + i]/size2;
+    x->output_prev_block[i] = out[i];
+    x->output_prev_block[half + i] = out[half + i];
   }
+}
 
+freebytes(expsig1, size2*sizeof(float));
+freebytes(revsig2, size2*sizeof(float));
   return(w+6);
 
 }
 
-static void cross_dsp(t_cross *x, t_signal **sp)
+static void cc_dsp(t_cc *x, t_signal **sp)
 {
     if(!x->is_freq_domain)
-	dsp_add(cross_perform_time_domain, 5,
+	dsp_add(cc_perform_time_domain, 5,
 		sp[0]->s_vec, sp[1]->s_vec, sp[2]->s_vec, sp[0]->s_n, x->delay);
     else
-	dsp_add(cross_perform_freq_domain, 5, x,
+	dsp_add(cc_perform_freq_domain, 5, x,
 		sp[0]->s_vec, sp[1]->s_vec, sp[2]->s_vec, sp[0]->s_n);
 }
 
-static void *cross_new(t_symbol *s, t_int argc, t_atom *argv)
+
+//  For using with running calculation, send a bang to clear the buffer
+//  and start over with calculations
+
+static void cc_bang(t_cc *x)
 {
-    t_cross *x = (t_cross *)pd_new(cross_class);
+  int i;
+  for(i=0;i<x->n;i++)
+    x->output_prev_block[i]=0;
+}
+
 
+static void *cc_new(t_symbol *s, t_int argc, t_atom *argv)
+{
+    t_cc *x = (t_cc *)pd_new(cc_class);
     if(atom_getsymbol(argv) == gensym("f")){
 	    x->is_freq_domain = 1;
-	    post("flib: cross: Frequency domain selected");
+	    post("flib: cc: Frequency domain selected");
+    }
+    else if(atom_getsymbol(argv) == gensym("r")){
+	    x->is_freq_domain = 2;
+	    post("flib: cc: Running frequency domain selected");
     }
     else {
 	x->delay = atom_getfloat(argv);
-	post("flib: cross: Time domain selected");
+	post("flib: cc: Time domain selected");
     } 
     inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_signal, &s_signal);
     outlet_new(&x->x_obj, &s_signal);
+    x->is_new=1;
+    x->buffer2Nsig1=NULL;
+    x->buffer2Nsig2=NULL;
+    x->output_prev_block=NULL;
     return (void *)x;
 }
 
+static void cc_free(t_cc *x)
+{
+  if     (x->buffer2Nsig1 != NULL)
+    freebytes(x->buffer2Nsig1, x->n*sizeof(float));
+  if     (x->buffer2Nsig2 != NULL)
+    freebytes(x->buffer2Nsig2, x->n*sizeof(float));
+  if     (x->output_prev_block != NULL)
+    freebytes(x->output_prev_block, x->n*sizeof(float));
+}
 
-void cross_tilde_setup(void) {
-    cross_class = class_new(gensym("cross~"),
-	    (t_newmethod)cross_new,
-	    0, sizeof(t_cross),
+
+void cc_tilde_setup(void) {
+    cc_class = class_new(gensym("cc~"),
+	    (t_newmethod)cc_new,
+	    (t_method)cc_free, sizeof(t_cc),
 	    CLASS_DEFAULT, A_GIMME, 0);
 
-    class_addmethod(cross_class,
-	    (t_method)cross_dsp, gensym("dsp"), 0);
-    CLASS_MAINSIGNALIN(cross_class, t_cross,f);
-    class_sethelpsymbol(cross_class, gensym("help-flib"));
+    class_addbang(cc_class, (t_method)cc_bang);
+    class_addmethod(cc_class,
+	    (t_method)cc_dsp, gensym("dsp"), 0);
+    CLASS_MAINSIGNALIN(cc_class, t_cc,f);
+    class_sethelpsymbol(cc_class, gensym("help-flib"));
 }
+
diff --git a/flib/src/flib.h b/flib/src/flib.h
index c43f63b..4e93d2e 100644
--- a/flib/src/flib.h
+++ b/flib/src/flib.h
@@ -21,7 +21,7 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 #include <stdlib.h>
 #include <string.h>
 
-#define VERSION "0.82"
+#define VERSION "0.83"
 
 void sc_tilde_setup(void);
 void ss_tilde_setup(void);
diff --git a/flib/src/irreg~.c b/flib/src/irreg~.c
index 7616c63..2b17828 100644
--- a/flib/src/irreg~.c
+++ b/flib/src/irreg~.c
@@ -18,7 +18,7 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
 
 
-/* calculates the spectral irreg of one frame according to Jensen et al 1999*/
+/* calculates the spectral irreg of one frame according to Jensen  1999. Original formula by Krimphoff et al 1994*/
 
 #include "flib.h"
 static t_class *irreg_class;
diff --git a/flib/src/sfm~.c b/flib/src/sfm~.c
index ee745c2..f2fbfdd 100644
--- a/flib/src/sfm~.c
+++ b/flib/src/sfm~.c
@@ -18,7 +18,7 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
 
 
-/* calculates spectral flatness measure as described by Tae Hong Park*/
+/* calculates spectral flatness measure: Geometric Mean/ Arithemtic Mean (In this case converted to a DB scale */
 
 #include "m_pd.h"
 #include <math.h>
diff --git a/flib/src/ss~.c b/flib/src/ss~.c
index 8ad1caa..83baf82 100644
--- a/flib/src/ss~.c
+++ b/flib/src/ss~.c
@@ -18,7 +18,7 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
 
 
-/* calculates the spectral smoothness of one frame according to Krimphoff et al. 1994*/
+/* calculates the spectral smoothness of one frame according to McAdams 1999. Takes magnitude spectrum as input*/
 
 #include "flib.h"
 static t_class *ss_class;
@@ -35,9 +35,9 @@ static t_int *ss_perform(t_int *w)
   t_int N = (t_int)(w[3]),M = N >> 1,n;
   t_float I = 0, buf[M];
   for(n = 0; n < M; buf[n++] = *in++);
-  for(n = 2; n < M - 1; n++){
+  for(n = 2; n < M - 1; n++){ 
 	  if(buf[n] != 0 && buf[n-1] != 0 && buf[n+1] != 0)
-		  I += (buf[n] - (buf[n-1] + buf[n] + buf[n+1]) / 3);
+		  I += (20 * log10(buf[n]) - (20 * log10(buf[n-1]) + 20 * log10(buf[n]) + 20 * log10(buf[n+1])) / 3);
   }
 
   outlet_float(x->x_obj.ob_outlet, I);