aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--flib/src/cc~.c199
-rw-r--r--flib/src/flib.h2
-rw-r--r--flib/src/irreg~.c2
-rw-r--r--flib/src/sfm~.c2
-rw-r--r--flib/src/ss~.c6
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);
generated by cgit v1.2.1 (git 2.18.0) at 2024-12-25 06:36:15 +0000