This commit was generated by cvs2svn to compensate for changes in r350,svn2git-root

which included commits to RCS files with non-trunk default branches.

svn path=/trunk/externals/creb/; revision=351

1 files changed, 272 insertions, 0 deletions

diff --git a/modules/bdiag.c b/modules/bdiag.c
new file mode 100644
index 0000000..37a0349
--- /dev/null
+++ b/modules/bdiag.c
@@ -0,0 +1,272 @@
+/*
+ *   bdiag.c  - block diagonal state space system
+ *   treats input dsp block as n parallel signals
+ *
+ *     s1 = (a * s1) + (b * s2) + u1;
+ *     s2 = (a * s2) - (b * s1) + u2;
+ *
+ *   Copyright (c) 2000-2003 by Tom Schouten
+ *
+ *   This program is free software; you can redistribute it and/or modify
+ *   it under the terms of the GNU General Public License as published by
+ *   the Free Software Foundation; either version 2 of the License, or
+ *   (at your option) any later version.
+ *
+ *   This program is distributed in the hope that it will be useful,
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *   GNU General Public License for more details.
+ *
+ *   You should have received a copy of the GNU General Public License
+ *   along with this program; if not, write to the Free Software
+ *   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ */
+
+
+#include "m_pd.h"
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define MAXORDER 64
+
+typedef struct bdiagctl
+{
+  t_float *c_state;
+  t_float *c_eigen;
+  t_int c_order;
+} t_bdiagctl;
+
+typedef struct bdiag
+{
+  t_object x_obj;
+  t_float x_f;
+  t_bdiagctl x_ctl;
+} t_bdiag;
+
+
+static float randfloat(void){
+  float r = rand ();
+  r /= (RAND_MAX/2);
+  r -= 1;
+  return r;
+
+}
+
+static void bdiag_random(t_bdiag *x)
+{
+  int i;
+
+  for (i=0; i<x->x_ctl.c_order; i++)
+    { 
+      x->x_ctl.c_state[i] = randfloat();
+    }
+  
+}
+
+
+static void bdiag_reset(t_bdiag *x)
+{
+  int i;
+
+  for (i=0; i<x->x_ctl.c_order; i++)
+    { 
+      x->x_ctl.c_state[i] = 0;
+    }
+  
+}
+
+
+
+
+
+
+static void bdiag_eigen(t_bdiag *x, t_floatarg index, t_floatarg aval, t_floatarg bval)
+{
+  int i = (int)index;
+  if (i<0) return;
+  if (i>=x->x_ctl.c_order/2) return;
+  x->x_ctl.c_eigen[2*i+0] = aval;
+  x->x_ctl.c_eigen[2*i+1] = bval;
+  
+}
+
+/* set decay time and frequency of pole at index */
+static void bdiag_timefreq(t_bdiag *x, t_floatarg index, t_floatarg time, t_floatarg freq)
+{
+  float r,a,b,n;
+  float sr = sys_getsr() / (float)x->x_ctl.c_order;
+
+  /* time in ms */
+  time *= 0.001;
+
+  if (time < 0.0f) time = 0.0f;
+  r = pow(0.001f, 1.0f / (time * sr));
+  if (r < 0.0f) r = 0.0f;
+  if (r > 1.0f) r = 1.0f;
+
+  a = cos(2*M_PI*freq/sr);
+  b = sin(2*M_PI*freq/sr);
+ 
+  /* normalize to be sure */
+  n = 1.0f / sqrt(a*a + b*b);
+  a *= n;
+  b *= n;
+
+  bdiag_eigen(x, index, r*a, r*b);
+}
+
+static void bdiag_preset(t_bdiag *x, t_floatarg preset)
+{
+  int p = preset;
+  int i;
+  float a, b, w, r;
+
+  switch(p){
+  case 0:
+    post("preset 0");
+    for (i=0; i<x->x_ctl.c_order/2; i++){
+      w = randfloat() * .001;
+      r = 1. - (((float)i + 1.)/1000.);
+      a = cos(w) * r;
+      b = sin(w) * r;
+      post("%f %f %f %f", w, r, a, b);
+      bdiag_eigen(x,i,a,b);
+    }
+    break;
+  case 1:
+  default:
+    break;
+  
+  }
+}
+
+static t_int *bdiag_perform(t_int *w)
+{
+
+
+  t_float *in     = (float *)(w[3]);
+  t_float *out    = (float *)(w[4]);
+  t_bdiagctl *ctl  = (t_bdiagctl *)(w[1]);
+
+  t_float *eigen = ctl->c_eigen;
+  t_float *state = ctl->c_state;
+  t_int n = (t_int)(w[2]);
+
+  t_float u1,u2,a,b,s1,s2,s1new,s2new;
+
+  int i;
+
+  for (i=0; i<n; i+=2)
+    {
+      u1 = *in++;
+      u2 = *in++;
+      a = *eigen++; /* real part */
+      b = *eigen++; /* imag part */
+      s1 = state[0];
+      s2 = state[1];
+
+
+      s1new = (a * s1) - (b * s2) + u1; /* update state */
+      s2new = (a * s2) + (b * s1) + u2;
+
+      *state++ = s1new; /* store state */
+      *state++ = s2new;
+
+      *out++ = s1new; /* output state */
+      *out++ = s2new;
+    }
+
+  return (w+5);
+}
+
+
+static void bdiag_dsp(t_bdiag *x, t_signal **sp)
+{
+
+  int n = sp[0]->s_n;
+  int i;
+
+  if (n == 1)
+    {
+      post("bdiag: doesnt work with blocksize == 1");
+      dsp_add_copy(sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
+    }
+  else
+    {
+      if (x->x_ctl.c_order != n)
+	{
+	  if (x->x_ctl.c_state) free(x->x_ctl.c_state);
+	  if (x->x_ctl.c_eigen) free(x->x_ctl.c_eigen);
+	  
+	  x->x_ctl.c_state = (t_float *)malloc(n*sizeof(t_float));
+	  x->x_ctl.c_eigen = (t_float *)malloc(n*sizeof(t_float));
+	  
+	  for(i=0;i<n;i++)
+	    {
+	      x->x_ctl.c_state[i] = 0;
+	      x->x_ctl.c_eigen[i] = 0;
+	    }
+	  
+	  x->x_ctl.c_order = n;
+	}
+      
+           
+      dsp_add(bdiag_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
+    }
+
+}
+                                  
+static void bdiag_free(t_bdiag *x)
+{
+
+  if (x->x_ctl.c_state) free(x->x_ctl.c_state);
+  if (x->x_ctl.c_eigen) free(x->x_ctl.c_eigen);
+
+
+}
+
+t_class *bdiag_class;
+
+
+static void *bdiag_new(t_floatarg permute)
+{
+    t_bdiag *x = (t_bdiag *)pd_new(bdiag_class);
+    int i, n=64;
+
+    outlet_new(&x->x_obj, gensym("signal")); 
+
+    x->x_ctl.c_state = (t_float *)malloc(n*sizeof(t_float));
+    x->x_ctl.c_eigen = (t_float *)malloc(n*sizeof(t_float));
+
+    for(i=0;i<n;i++)
+      {
+	x->x_ctl.c_state[i] = 0;
+	x->x_ctl.c_eigen[i] = 0;
+      }
+    
+    x->x_ctl.c_order = n;
+
+
+    return (void *)x;
+
+
+}
+
+
+void bdiag_tilde_setup(void)
+{
+  //post("bdiag~ v0.1");
+    bdiag_class = class_new(gensym("bdiag~"), (t_newmethod)bdiag_new,
+    	(t_method)bdiag_free, sizeof(t_bdiag), 0, A_DEFFLOAT, 0);
+    CLASS_MAINSIGNALIN(bdiag_class, t_bdiag, x_f);
+    class_addmethod(bdiag_class, (t_method)bdiag_random, gensym("random"), 0);
+    class_addmethod(bdiag_class, (t_method)bdiag_random, gensym("bang"), 0);
+    class_addmethod(bdiag_class, (t_method)bdiag_reset, gensym("reset"), 0);
+    class_addmethod(bdiag_class, (t_method)bdiag_dsp, gensym("dsp"), 0); 
+
+    class_addmethod(bdiag_class, (t_method)bdiag_eigen, gensym("eigen"), A_DEFFLOAT, A_DEFFLOAT, A_DEFFLOAT, 0); 
+    class_addmethod(bdiag_class, (t_method)bdiag_timefreq, gensym("timefreq"), A_DEFFLOAT, A_DEFFLOAT, A_DEFFLOAT, 0); 
+    class_addmethod(bdiag_class, (t_method)bdiag_preset, gensym("preset"), A_DEFFLOAT, 0); 
+}
+