Added frequency domain correlation to cross~

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

3 files changed, 158 insertions, 15 deletions

diff --git a/flib/makefile b/flib/makefile
index 0974eb0..babda2f 100644
--- a/flib/makefile
+++ b/flib/makefile
@@ -1,6 +1,6 @@
-# Makefile hacked together from various others, mostly the old zexy makefile I think. All credit to IOhannes Zmoelnig and anyone else who understands make better than me.
+# Makefile hacked together from various others, mostly the old zexy makefile and pique. All credit to IOhannes Zmoelnig, Miller Puckette and anyone else who understands make better than me.
 #
-# Should work with Linux and OS X. No idea about Windows.
+# Should work with Linux and OS X. No idea about Windows and Irix.
 
 NAME=flib
 CSYM=flib
@@ -152,6 +152,71 @@ LINUXINCLUDE =  -I/usr/include
 	ld -export_dynamic  -shared -o $(NAME).pd_linux $(NAME).o $(LINUXOBJECTS) -lc -lm
 	strip --strip-unneeded $(NAME).pd_linux
 
+# ----------------------- IRIX 5.x -----------------------
+
+pd_irix5: src/$(NAME).pd_irix5
+
+.SUFFIXES: .pd_irix5
+
+SGIOBJECTS = $(TARGETS:%=%.o)
+
+SGICFLAGS5 = -o32 -DPD -DUNIX -DIRIX -O2
+
+SGIINCLUDE =  -I/usr/include/ -I/usr/local/include/
+
+.c.pd_irix5:
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/sc~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/ss~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/irreg~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/mspec~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/peak~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/pspec~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/sfm~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/trist~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/pp~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/bmax~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/melf~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/clean~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/wdv~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/hca~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/cross~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c $*.c
+	ld -elf -shared -rdata_shared -o $(NAME).pd_irix5 $(NAME).o $(LINUXOBJECTS) -lm -lc
+	rm $*.o
+
+# ----------------------- IRIX 6.x -----------------------
+
+pd_irix6: src/$(NAME).pd_irix6
+
+.SUFFIXES: .pd_irix6
+
+SGICFLAGS6 = -n32 -DPD -DUNIX -DIRIX -DN32 -woff 1080,1064,1185 \
+	-OPT:roundoff=3 -OPT:IEEE_arithmetic=3 -OPT:cray_ivdep=true \
+	-Ofast=ip32
+
+SGIOBJECTS = $(TARGETS:%=%.o)
+
+SGIINCLUDE =  -I/usr/include/ -I/usr/local/include/
+
+.c.pd_irix6:
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/ss~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/irreg~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/mspec~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/peak~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/pspec~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/sfm~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/trist~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/pp~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/bmax~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/melf~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/clean~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/wdv~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/hca~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c src/cross~.c
+	$(CC) $(SGICFLAGS) $(SGIINCLUDE) -c $*.c
+	ld -n32 -IPA -shared -rdata_shared -o $*.pd_irix6 $(NAME).o $(LINUXOBJECTS) -lm -lc
+	rm $*.o
+
 # ----------------------------------------------------------
 
 install:
diff --git a/flib/src/cross~.c b/flib/src/cross~.c
index d7ea4bb..2d06d9b 100644
--- a/flib/src/cross~.c
+++ b/flib/src/cross~.c
@@ -17,8 +17,14 @@
    */
 
 
-/*Calculate the (non optimized) cross correlation of two signal vectors*/
-/*Based on code by Phil Bourke */
+/*Calculate the cross 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 an argument of 'f' for the frequency domain implementation*/
 
 
 #include "flib.h"
@@ -30,9 +36,10 @@ typedef struct _cross {
     t_object  x_obj;
     t_float f;
     t_int delay;
+    t_int is_freq_domain;
 } t_cross;
 
-static t_int *cross_perform(t_int *w)
+static t_int *cross_perform_time_domain(t_int *w)
 {
     t_sample *x = (t_sample *)(w[1]);
     t_sample *y = (t_sample *)(w[2]);
@@ -45,12 +52,9 @@ static t_int *cross_perform(t_int *w)
    if(maxdelay > N * .5){
 	maxdelay = N * .5;
 	post("cross~: invalid maxdelay, must be <= blocksize/2");
-   } 
-   else if (!maxdelay){
-       maxdelay = N * .5;
    }
  
-   /* Calculate the mean of the two series x[], y[] */
+/* Calculate the mean of the two series x[], y[] */
    mx = 0;
    my = 0;   
    for (i=0;i<N;i++) {
@@ -92,16 +96,90 @@ static t_int *cross_perform(t_int *w)
 
 }
 
+
+t_int *cross_perform_freq_domain(t_int *w)
+{
+  t_cross *x = (t_cross *)(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++)
+  {
+    expsig1[i]=sig1[i];
+    revsig2[i]=0;
+  }
+  for (i=qtr; i < thrqtr ; i++)
+  {
+    expsig1[i]=sig1[i];
+    revsig2[i]=sig2[size-i];
+  }
+  for (i=thrqtr; i < size ; i++)
+  {
+    expsig1[i]=sig1[i];
+    revsig2[i]=0;
+  }
+
+  mayer_realfft(size, expsig1);
+  mayer_realfft(size, revsig2);
+  expsig1[0]*=revsig2[0];
+  expsig1[k]*=revsig2[k];
+  for(i=1; i < k; 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];
+  }
+  
+  mayer_realifft(size, expsig1);
+  for(i=0; i < size; i++)
+  {
+    out[i]=expsig1[i];
+  }
+
+  return(w+6);
+
+}
+
 static void cross_dsp(t_cross *x, t_signal **sp)
 {
-    dsp_add(cross_perform, 5,
-	    sp[0]->s_vec, sp[1]->s_vec, sp[2]->s_vec, sp[0]->s_n, x->delay);
+    if(!x->is_freq_domain)
+	dsp_add(cross_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,
+		sp[0]->s_vec, sp[1]->s_vec, sp[2]->s_vec, sp[0]->s_n);
 }
 
-static void *cross_new(t_floatarg f)
+static void *cross_new(t_symbol *s, t_int argc, t_atom *argv)
 {
     t_cross *x = (t_cross *)pd_new(cross_class);
-    x->delay = (t_int)f;
+
+    if(atom_getsymbol(argv) == gensym("f")){
+	    x->is_freq_domain = 1;
+	    post("flib: cross: Frequency domain selected");
+    }
+    else {
+	x->delay = atom_getfloat(argv);
+	post("flib: cross: Time domain selected");
+    } 
     inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_signal, &s_signal);
     outlet_new(&x->x_obj, &s_signal);
     return (void *)x;
@@ -112,7 +190,7 @@ void cross_tilde_setup(void) {
     cross_class = class_new(gensym("cross~"),
 	    (t_newmethod)cross_new,
 	    0, sizeof(t_cross),
-	    CLASS_DEFAULT, A_DEFFLOAT, 0);
+	    CLASS_DEFAULT, A_GIMME, 0);
 
     class_addmethod(cross_class,
 	    (t_method)cross_dsp, gensym("dsp"), 0);
diff --git a/flib/src/flib.h b/flib/src/flib.h
index fba78f0..868f683 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.8"
+#define VERSION "0.82"
 
 void sc_tilde_setup(void);
 void ss_tilde_setup(void);