25 May 2006


      Update of /cvsroot/pure-data/externals/creb/modules
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv29110
Added Files:
    abs~.c bdiag~.c bfft~.c bitsplit~.c blocknorm~.c cheby~.c 
    cmath~.c diag~.c dist~.c dwt~.c dynwav~.c eadsr~.c ead~.c 
    ear~.c eblosc~.c fdn~.c junction~.c lattice~.c matrix~.c 
    permut~.c qmult~.c qnorm~.c ramp~.c resofilt~.c sbosc~.c 
    scrollgrid1D~.c statwav~.c tabreadmix~.c window~.c xfm~.c 
Removed Files:
    abs.c bdiag.c bfft.c bitsplit.c blocknorm.c cheby.c cmath.c 
    diag.c dist.c dwt.c dynwav.c ead.c eadsr.c ear.c eblosc.c 
    fdn.c junction.c lattice.c matrix.c permut.c qmult.c qnorm.c 
    ramp.c resofilt.c sbosc.c scrollgrid1D.c statwav.c 
    tabreadmix.c window.c xfm.c 
Log Message:
renamed files to match their class names
--- NEW FILE: ear~.c ---
/*
 *   ear.c  -  exponential attack release envelope 
 *   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 "extlib_util.h"
typedef struct earctl
{
  t_float c_attack;
  t_float c_release;
  t_float c_state;
  t_float c_target;
} t_earctl;
typedef struct ear
{
  t_object x_obj;
  t_float x_sr;
  t_earctl x_ctl;
} t_ear;
void ear_attack(t_ear *x, t_floatarg f)
{
  x->x_ctl.c_attack = milliseconds_2_one_minus_realpole(f);
}
void ear_release(t_ear *x, t_floatarg f)
{
  x->x_ctl.c_release = milliseconds_2_one_minus_realpole(f);
}
void ear_start(t_ear *x)
{
    x->x_ctl.c_target = 1;
    x->x_ctl.c_state = 0.0f;
}
void ear_stop(t_ear *x)
{
    x->x_ctl.c_target = 0;
}
void ear_float(t_ear *x, t_floatarg f)
{
    if (f == 0.0f) ear_stop(x);
    else ear_start(x);
}
static t_int *ear_perform(t_int *w)
{
    t_float *out   = (float *)(w[3]);
    t_earctl *ctl = (t_earctl *)(w[1]);
    t_float attack = ctl->c_attack;
    t_float release  = ctl->c_release;
    t_float state  = ctl->c_state;
    t_float target = ctl->c_target;
    t_int n = (t_int)(w[2]);
t_int i;
if (target) /* attack phase */
      for (i = 0; i < n; i++)
    {
      *out++ = state;
      state += attack*(1 - state);
    } 
    else /* release phase */
      for (i = 0; i < n; i++)
    {
      *out++ = state;
      state -= release*state;
    }
ctl->c_state = IS_DENORMAL(state) ? 0 : state;
    return (w+4);
}
static void ear_dsp(t_ear *x, t_signal **sp)
{
    x->x_sr = sp[0]->s_sr;
    dsp_add(ear_perform, 3, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec);
}                                  
void ear_free(void)
{
}
t_class *ear_class;  /* attack - release */
void *ear_new(t_floatarg attack, t_floatarg release)
{
    t_ear *x = (t_ear *)pd_new(ear_class);
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("attack"));  
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("release"));
    outlet_new(&x->x_obj, gensym("signal"));
ear_attack(x,attack);
    ear_release(x,release);
    x->x_ctl.c_state  = 0;
    x->x_ctl.c_target = 0;
return (void *)x;
}
void ear_tilde_setup(void)
{
    //post("ear~ v0.1");
    ear_class = class_new(gensym("ear~"), (t_newmethod)ear_new,
    	(t_method)ear_free, sizeof(t_ear), 0, A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(ear_class, (t_method)ear_float, gensym("float"), A_FLOAT, 0);
    class_addmethod(ear_class, (t_method)ear_start, gensym("start"), 0);
    class_addmethod(ear_class, (t_method)ear_start, gensym("bang"), 0);
    class_addmethod(ear_class, (t_method)ear_stop, gensym("stop"), 0);
    class_addmethod(ear_class, (t_method)ear_dsp, gensym("dsp"), 0); 
    class_addmethod(ear_class, 
    	    (t_method)ear_attack, gensym("attack"), A_FLOAT, 0);
    class_addmethod(ear_class, 
    	    (t_method)ear_release, gensym("release"), A_FLOAT, 0);
}
--- NEW FILE: diag~.c ---
/*
 *   diag.c  - 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>
#ifndef IS_DENORMAL
#include "extlib_util.h"
#endif
#define MAXORDER 64
typedef struct diagctl
{
  t_float *c_state;
  t_float *c_eigen;
  t_int c_order;
} t_diagctl;
typedef struct diag
{
  t_object x_obj;
  t_float x_f;
  t_diagctl x_ctl;
} t_diag;
static float randfloat(void){
  float r = rand ();
  r /= (RAND_MAX/2);
  r -= 1;
  return r;
}
static void diag_eigen(t_diag *x, t_floatarg index, t_floatarg val)
{
  int i = (int)index;
  if (i<0) return;
  if (i>=x->x_ctl.c_order) return;
  x->x_ctl.c_eigen[i] = val;
}
/* set decay time of pole at index */
static void diag_time(t_diag *x, t_floatarg index, t_floatarg time)
{
  float r;
/* time in ms */
  time *= 0.001;
if (time < 0.0f) time = 0.0f;
  r = pow(0.001f, (float)x->x_ctl.c_order / (time * sys_getsr()));
  if (r < 0.0f) r = 0.0f;
  if (r > 1.0f) r = 1.0f;
diag_eigen(x, index, r);
}
static void diag_reset(t_diag *x)
{
  int i;
for (i=0; i<x->x_ctl.c_order; i++)
    { 
      x->x_ctl.c_state[i] = 0;
    }
}
static void diag_random(t_diag *x)
{
  int i;
for (i=0; i<x->x_ctl.c_order; i++)
    { 
      x->x_ctl.c_state[i] = randfloat();
    }
}
static t_int *diag_perform(t_int *w)
{
t_float *in     = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_diagctl *ctl  = (t_diagctl *)(w[1]);
t_float *eigen = ctl->c_eigen;
  t_float *state = ctl->c_state;
  t_int n = (t_int)(w[2]);
  t_float newstate;
int i;
for (i=0; i<n; i++)
    {
      newstate = (*state + *in) * (*eigen);
      newstate = IS_DENORMAL(newstate) ? 0 : newstate;
      *state   = newstate;
      *out     = newstate;
in++;
      out++;
      state++;
      eigen++;
    }
return (w+5);
}
static void diag_dsp(t_diag *x, t_signal **sp)
{
int n = sp[0]->s_n;
  int i;
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(diag_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
}
static void diag_free(t_diag *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 *diag_class;
static void *diag_new(t_floatarg permute)
{
    t_diag *x = (t_diag *)pd_new(diag_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 diag_tilde_setup(void)
{
  //post("diag~ v0.1");
    diag_class = class_new(gensym("diag~"), (t_newmethod)diag_new,
    	(t_method)diag_free, sizeof(t_diag), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(diag_class, t_diag, x_f);
    class_addmethod(diag_class, (t_method)diag_dsp, gensym("dsp"), 0); 
    class_addmethod(diag_class, (t_method)diag_reset, gensym("reset"), 0); 
    class_addmethod(diag_class, (t_method)diag_random, gensym("random"), 0); 
    class_addmethod(diag_class, (t_method)diag_random, gensym("bang"), 0); 
    class_addmethod(diag_class, (t_method)diag_eigen, gensym("eigen"), A_DEFFLOAT, A_DEFFLOAT, 0); 
    class_addmethod(diag_class, (t_method)diag_time, gensym("time"), A_DEFFLOAT, A_DEFFLOAT, 0);
}
--- statwav.c DELETED ---
--- NEW FILE: bdiag~.c ---
/*
 *   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>
#ifndef IS_DENORMAL
#include "extlib_util.h"
#endif
#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;
s1new = IS_DENORMAL(s1new) ? 0 : s1new; /* clear denormals */
      s2new = IS_DENORMAL(s2new) ? 0 : s2new;
*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); 
}
--- tabreadmix.c DELETED ---
--- NEW FILE: resofilt~.c ---
/*
 *   resofilt.c  -  some high quality time variant filters
 *   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.
 */
/* some (expensive) time variant iir filters, 
   i.e. moog 4-pole transfer function using the impulse
   invariant transform with self osc and
   signal freq and reso inputs */
#include "m_pd.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>  
#include "filters.h"
typedef struct resofiltctl
{
    t_float c_state[4]; /* filter state */
    t_float c_f_prev;
    t_float c_r_prev;
} t_resofiltctl;
typedef struct resofilt
{
    t_object x_obj;
    t_float x_f;
    t_resofiltctl x_ctl;
    t_perfroutine x_dsp;
} t_resofilt;
static inline void _sat_state(t_float *x)
{
    const float norm_squared_max = 1.0f;
    float norm_squared = x[0] * x[0] + x[1] * x[1];
if (norm_squared > norm_squared_max){
    float scale = 1.0f / sqrt(norm_squared);
    x[0] *= scale;
    x[1] *= scale;
    }
}
/* the moog vcf discretized with an impulse invariant transform */
static t_int *resofilt_perform_fourpole(t_int *w)
{
t_resofiltctl *ctl = (t_resofiltctl *)(w[1]);
  t_int n            = (t_int)(w[2]);
  t_float *in        = (float *)(w[3]);
  t_float *freq      = (float *)(w[4]);
  t_float *reso      = (float *)(w[5]);
  t_float *out       = (float *)(w[6]);
int i;
  t_float inv_n = 1.0f / ((float)n);
  t_float inv_sr = 1.0f / sys_getsr();
t_float phasor[2], phasor_rot[2];
  t_float radior[2], radior_rot[2];
  t_float scalor, scalor_inc;
t_float f,r,freq_rms,reso_rms;
  t_float f_prev = ctl->c_f_prev;
  t_float r_prev = ctl->c_r_prev;
/* use rms of input to drive freq and reso
     this is such that connecting a dsp signal to the inlets has a reasonable result,
     even if the inputs are oscillatory. the rms values will be interpolated linearly
     (that is, linearly in the "analog" domain, so exponentially in the z-domain)  */
reso_rms = freq_rms = 0.0f;
  for (i=0; i<n; i++){
      t_float _freq = *freq++; /* first input is the reso frequency (absolute) */
      t_float _reso = *reso++; /* second input is the resonnance (0->1), >1 == self osc */
      freq_rms += _freq * _freq;
      reso_rms += _reso * _reso;
  }
  freq_rms = sqrt(freq_rms * inv_n) * inv_sr;
  reso_rms = sqrt(reso_rms * inv_n);
  f = (freq_rms > 0.5f) ? 0.5f : freq_rms;
  r = sqrt(sqrt(reso_rms));
/* calculate the new pole locations
     we use an impulse invariant transform: the moog vcf poles are located at
     s_p = (-1 +- r \sqrt{+- j}, with r = (k/4)^(1/4) \in [0,1]
the poles are always complex, so we can use an orthogonal implementation
     both conj pole pairs have the same angle, so we can use one phasor and 2 radii
  */
/* compute phasor, radius and update eqs
     since these are at k-rate, the expense is justifyable */
  phasor[0] = cos(2.0 * M_PI * r_prev * f_prev);
  phasor[1] = sin(2.0 * M_PI * r_prev * f_prev);
  phasor_rot[0] = cos(2.0 * M_PI * (r*f - r_prev*f_prev) * inv_n);
  phasor_rot[1] = sin(2.0 * M_PI * (r*f - r_prev*f_prev) * inv_n);
radior[0] = exp(f_prev * (-1.0 + r_prev)); /* dominant radius */
  radior[1] = exp(f_prev * (-1.0 - r_prev));
  radior_rot[0] = exp((f*(-1.0f + r) - f_prev * (-1.0 + r_prev)) * inv_n);
  radior_rot[1] = exp((f*(-1.0f - r) - f_prev * (-1.0 - r_prev)) * inv_n);
/* moog like scaling */
  if (1){
      float r2 = r_prev * r_prev;
      float r4 = r2 * r2;
      scalor = 1.0f + (1.0f + 4.0f * r4);
      r2 = r * r;
      r4 = r2 * r2;
      scalor_inc = ((1.0f + (1.0f + 4.0f * r4)) - scalor) * inv_n;
  }
/* no scaling */
  else{
      scalor = 1.0f;
      scalor_inc = 0.0f;
  }
ctl->c_f_prev = f;
  ctl->c_r_prev = r;
/* perform filtering */
  for (i=0; i<n; i++){
      float poleA[2], poleB[2];
      float *stateA = ctl->c_state;
      float *stateB = ctl->c_state+2;
float x;
/* compute poles */
      poleA[0] = radior[0] * phasor[0];
      poleA[1] = radior[0] * phasor[1];
      poleB[0] = radior[1] * phasor[0];
      poleB[1] = radior[1] * phasor[1];
/* perform filtering: use 2 DC normalized complex one-poles:
     y[k] = x[k] + a(y[k-1] - x[k]) or y(z) = (1-a)/(1-az^{-1}) x(z) */
x = *in++ * scalor;
/* nondominant pole first */
      stateB[0] -= x;
      vcmul2(stateB, poleB);
      x = stateB[0] += x;
/* dominant pole second */
      stateA[0] -= x;
      vcmul2(stateA, poleA);
      x = stateA[0] += x;
*out++ = x;
/* saturate (normalize if pow > 1) state to prevent explosion and to allow self-osc */
      _sat_state(stateA);
      _sat_state(stateB);
/* interpolate filter coefficients */
      vcmul2(phasor, phasor_rot);
      radior[0] *= radior_rot[0];
      radior[1] *= radior_rot[1];
      scalor += scalor_inc;
}
return (w+7);
}
/* 303-style modified moog vcf (3-pole) */
static t_int *resofilt_perform_threepole(t_int *w)
{
t_resofiltctl *ctl = (t_resofiltctl *)(w[1]);
  t_int n            = (t_int)(w[2]);
  t_float *in        = (float *)(w[3]);
  t_float *freq      = (float *)(w[4]);
  t_float *reso      = (float *)(w[5]);
  t_float *out       = (float *)(w[6]);
int i;
  t_float inv_n = 1.0f / ((float)n);
  t_float inv_sr = 1.0f / sys_getsr();
t_float phasor[2], phasor_rot[2];
  t_float radior[2], radior_rot[2];
  t_float scalor, scalor_inc;
t_float f,r,freq_rms,reso_rms;
  t_float f_prev = ctl->c_f_prev;
  t_float r_prev = ctl->c_r_prev;
t_float sqrt5 = sqrtf(5.0f);
/* use rms of input to drive freq and reso */
  reso_rms = freq_rms = 0.0f;
  for (i=0; i<n; i++){
      t_float _freq = *freq++; /* first input is the reso frequency (absolute) */
      t_float _reso = *reso++; /* second input is the resonnance (0->1), >1 == self osc */
      freq_rms += _freq * _freq;
      reso_rms += _reso * _reso;
  }
  freq_rms = sqrt(freq_rms * inv_n) * inv_sr;
  reso_rms = sqrt(reso_rms * inv_n);
  f = (freq_rms > 0.5f) ? 0.25f : freq_rms * 0.5f;
  r = cbrt(reso_rms);
/* calculate the new pole locations
     we use an impulse invariant transform: the 303 vcf poles are located at
     s_p = omega(-1 + r sqrt(5) e^{pi/3(1+2p)})
real pole: omega * (-1 -r)
     complex pole:
       real = omega (-1 + r)
       imag = omega (+- 2 r)
this is a strange beast. legend goes it was "invented" by taking the moog vcf
     and moving one cap up, such that the not-so controllable 3-pole that emerged
     would avoid the moog patent..
so, the sound is not so much the locations of the poles, but how the filter
     reacts to time varying controls. i.e. the pole movement with constant reso,
     used in the tb-303.
*/
/* compute phasor, radius and update eqs */
  phasor[0] = cos(2.0 * M_PI * r_prev * f_prev * 2.0f);
  phasor[1] = sin(2.0 * M_PI * r_prev * f_prev * 2.0f);
  phasor_rot[0] = cos(2.0 * M_PI * (r*f - r_prev*f_prev) * 2.0f * inv_n);
  phasor_rot[1] = sin(2.0 * M_PI * (r*f - r_prev*f_prev) * 2.0f * inv_n);
radior[0] = exp(f_prev * (-1.0 + r_prev)); /* dominant radius */
  radior[1] = exp(f_prev * (-1.0 - sqrt5 * r_prev));
  radior_rot[0] = exp((f*(-1.0f + r) - f_prev * (-1.0 + r_prev)) * inv_n);
  radior_rot[1] = exp((f*(-1.0f - r) - f_prev * (-1.0 - sqrt5 * r_prev)) * inv_n);
/* 303 like scaling */
  if (1){
      float r3 = r_prev * r_prev * r_prev;
      scalor = 1.0f + (1.0f + 3.0f * r_prev);
      r3 = r * r * r;
      scalor_inc = ((1.0f + (1.0f + 3.0f * r3)) - scalor) * inv_n;
  }
/* no scaling */
  else{
      scalor = 1.0f;
      scalor_inc = 0.0f;
  }
ctl->c_f_prev = f;
  ctl->c_r_prev = r;
ctl->c_state[3] = 0.0f;
  /* perform filtering */
  for (i=0; i<n; i++){
      float poleA[2], poleB[2];
      float *stateA = ctl->c_state;
      float *stateB = ctl->c_state+2;
float x;
/* compute poles */
      poleA[0] = radior[0] * phasor[0];
      poleA[1] = radior[0] * phasor[1];
poleB[0] = radior[1];
/* perform filtering: use (real part of) 2 DC normalized complex one-poles:
     y[k] = x[k] + a(y[k-1] - x[k]) or y(z) = (1-a)/(1-az^{-1}) x(z) */
x = *in++ * scalor;
/* nondominant pole first */
      stateB[0] -= x;
      stateB[0] *= poleB[0];
      x = stateB[0] += x;
/* dominant pole second */
      stateA[0] -= x;
      vcmul2(stateA, poleA);
      x = stateA[0] += x;
*out++ = x;
/* saturate (normalize if pow > 1) state to prevent explosion and to allow self-osc */
      _sat_state(stateA);
      _sat_state(stateB);
/* interpolate filter coefficients */
      vcmul2(phasor, phasor_rot);
      radior[0] *= radior_rot[0];
      radior[1] *= radior_rot[1];
      scalor += scalor_inc;
}
return (w+7);
}
static void resofilt_dsp(t_resofilt *x, t_signal **sp)
{
  int n = sp[0]->s_n;
dsp_add(x->x_dsp,
      6, 
      &x->x_ctl, 
      sp[0]->s_n, 
      sp[0]->s_vec, 
      sp[1]->s_vec, 
      sp[2]->s_vec, 
      sp[3]->s_vec);
}                                  
static void resofilt_free(t_resofilt *x)
{
}
t_class *resofilt_class;
static void *resofilt_new(t_floatarg algotype)
{
    t_resofilt *x = (t_resofilt *)pd_new(resofilt_class);
/* in */
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));  
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));
/* out */
    outlet_new(&x->x_obj, gensym("signal"));
/* set algo type */
    if (algotype == 3.0f){
    post("resofilt~: 3-pole lowpass ladder vcf");
    x->x_dsp = resofilt_perform_threepole;
    }
    else {
    post("resofilt~: 4-pole lowpass ladder vcf");
    x->x_dsp = resofilt_perform_fourpole;
    }
return (void *)x;
}
void resofilt_tilde_setup(void)
{
    resofilt_class = class_new(gensym("resofilt~"), (t_newmethod)resofilt_new,
    	(t_method)resofilt_free, sizeof(t_resofilt), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(resofilt_class, t_resofilt, x_f);
    class_addmethod(resofilt_class, (t_method)resofilt_dsp, gensym("dsp"), 0); 
}
--- NEW FILE: ramp~.c ---
/*
 *   ramp.c  - retriggerable counter for dsp signals 
 *   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>
typedef struct rampctl
{
    t_float c_offset;
    t_int c_blockscale;
} t_rampctl;
typedef struct ramp
{
    t_object x_obj;
    t_float x_f;
    t_rampctl x_ctl;
} t_ramp;
void ramp_offset(t_ramp *x,  t_floatarg f)
{
    x->x_ctl.c_offset = f;
}
void ramp_bang(t_ramp *x)
{
    ramp_offset(x, 0);
}
static t_int *ramp_perform(t_int *w)
{
    t_float *out    = (float *)(w[3]);
    t_rampctl *ctl  = (t_rampctl *)(w[1]);
    t_int i;
    t_int n = (t_int)(w[2]);
    t_float x;
t_float scale = ctl->c_blockscale ? 1.0f / (float)n : 1.0f;
x = ctl->c_offset;
for (i = 0; i < n; i++)
    {
        *out++ = ((float)x++) * scale;
    }
ctl->c_offset = x; /* save state */
return (w+4);
}
static void ramp_dsp(t_ramp *x, t_signal **sp)
{
    dsp_add(ramp_perform, 3, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec);
}                                  
void ramp_free(void)
{
}
t_class *ramp_class;
void *ramp_new(void)
{
    t_ramp *x = (t_ramp *)pd_new(ramp_class);
    outlet_new(&x->x_obj, gensym("signal")); 
    x->x_ctl.c_blockscale = 0;
    ramp_bang(x);
    return (void *)x;
}
void *blockramp_new(void)
{
    t_ramp *x = (t_ramp *)ramp_new();
    x->x_ctl.c_blockscale = 1;
    return (void *)x;
}
void ramp_tilde_setup(void)
{
  //post("ramp~ v0.1");
    ramp_class = class_new(gensym("ramp~"), (t_newmethod)ramp_new,
    	(t_method)ramp_free, sizeof(t_ramp), 0, 0);
class_addcreator((t_newmethod)blockramp_new, gensym("blockramp~"), A_NULL);
class_addmethod(ramp_class, (t_method)ramp_bang, gensym("bang"), 0);
    class_addmethod(ramp_class, (t_method)ramp_dsp, gensym("dsp"), 0); 
    class_addfloat(ramp_class, (t_method)ramp_offset); 
}
--- qnorm.c DELETED ---
--- cmath.c DELETED ---
--- ear.c DELETED ---
--- resofilt.c DELETED ---
--- NEW FILE: qmult~.c ---
/*
 *   qmult.c  - quaternion multiplication dsp object 
 *   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 <stdlib.h>
#include <stdio.h>
#include <math.h>
typedef struct qmultctl
{
  t_float *c_inputleft[4];
  t_float *c_inputright[4];
  t_float *c_output[4];
} t_qmultctl;
typedef struct qmult
{
  t_object x_obj;
  t_float x_f;
  t_qmultctl x_ctl;
} t_qmult;
static t_int *qmult_perform(t_int *word)
{
t_qmultctl *ctl     = (t_qmultctl *)(word[1]);
  t_int n             = (t_int)(word[2]);
  t_int i;
t_float *in0l        = ctl->c_inputleft[0];
  t_float *in1l        = ctl->c_inputleft[1];
  t_float *in2l        = ctl->c_inputleft[2];
  t_float *in3l        = ctl->c_inputleft[3];
t_float *in0r        = ctl->c_inputright[0];
  t_float *in1r        = ctl->c_inputright[1];
  t_float *in2r        = ctl->c_inputright[2];
  t_float *in3r        = ctl->c_inputright[3];
t_float *out0       = ctl->c_output[0];
  t_float *out1       = ctl->c_output[1];
  t_float *out2       = ctl->c_output[2];
  t_float *out3       = ctl->c_output[3];
t_float wl, xl, yl, zl;
  t_float wr, xr, yr, zr;
  t_float w, x, y, z;
for (i=0;i<n;i++)
    {
/* read input quaternions */
    wl = *in0l++;
    xl = *in1l++;
    yl = *in2l++;
    zl = *in3l++;
wr = *in0r++;
    xr = *in1r++;
    yr = *in2r++;
    zr = *in3r++;
/* multiply quaternions */
    w = wl * wr;
    x = wl * xr;
    y = wl * yr;
    z = wl * zr;
w -= xl * xr;
    x += xl * wr;
    y -= xl * zr;
    z += xl * yr;
w -= yl * yr;
    x += yl * zr;
    y += yl * wr;
    z -= yl * xr;
w -= zl * zr;
    x -= zl * yr;
    y += zl * xr;
    z += zl * wr;
/* write output quaternion */
    *out0++ = w;
    *out1++ = x;
    *out2++ = y;
    *out3++ = z;
    }
return (word+3);
}
static void qmult_dsp(t_qmult *x, t_signal **sp)
{
int i;
  for (i=0;i<4;i++)
    {
      x->x_ctl.c_inputleft[i] = sp[i]->s_vec;
      x->x_ctl.c_inputright[i] = sp[i+4]->s_vec;
      x->x_ctl.c_output[i] = sp[i+8]->s_vec;
    }
dsp_add(qmult_perform, 2, &x->x_ctl, sp[0]->s_n);
}
static void qmult_free(t_qmult *x)
{
}
t_class *qmult_class;
static void *qmult_new(t_floatarg channels)
{
    int i;
    t_qmult *x = (t_qmult *)pd_new(qmult_class);
for (i=1;i<8;i++) inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal")); 
    for (i=0;i<4;i++) outlet_new(&x->x_obj, gensym("signal"));
return (void *)x;
}
void qmult_tilde_setup(void)
{
  //post("qmult~ v0.1");
    qmult_class = class_new(gensym("qmult~"), (t_newmethod)qmult_new,
    	(t_method)qmult_free, sizeof(t_qmult), 0, 0);
    CLASS_MAINSIGNALIN(qmult_class, t_qmult, x_f); 
    class_addmethod(qmult_class, (t_method)qmult_dsp, gensym("dsp"), 0);
}
--- xfm.c DELETED ---
--- NEW FILE: eadsr~.c ---
/*
 *   eadsr.c  -  exponential attack decay sustain release envelope 
 *   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 "extlib_util.h"
typedef struct eadsrctl
{
  t_float c_attack;
  t_float c_decay;
  t_float c_sustain;
  t_float c_release;
  t_float c_state;
  t_float c_target;
} t_eadsrctl;
typedef struct eadsr
{
  t_object x_obj;
  t_float x_sr;
  t_eadsrctl x_ctl;
} t_eadsr;
void eadsr_attack(t_eadsr *x, t_floatarg f)
{
  x->x_ctl.c_attack = milliseconds_2_one_minus_realpole(f);
}
void eadsr_decay(t_eadsr *x, t_floatarg f)
{
  x->x_ctl.c_decay = milliseconds_2_one_minus_realpole(f);
}
void eadsr_sustain(t_eadsr *x, t_floatarg f)
{
  if (f>ENVELOPE_MAX) f = ENVELOPE_MAX;
  if (f<ENVELOPE_MIN) f = ENVELOPE_MIN;
x->x_ctl.c_sustain = f;
}
void eadsr_release(t_eadsr *x, t_floatarg f)
{
  x->x_ctl.c_release = milliseconds_2_one_minus_realpole(f);
}
void eadsr_start(t_eadsr *x)
{
    x->x_ctl.c_target = 1;
    x->x_ctl.c_state = 0.0f;
}
void eadsr_stop(t_eadsr *x)
{
    x->x_ctl.c_target = 0;
}
void eadsr_float(t_eadsr *x, t_floatarg f)
{
    if (f == 0.0f) eadsr_stop(x);
    else eadsr_start(x);
}
static t_int *eadsr_perform(t_int *w)
{
    t_float *out    = (float *)(w[3]);
    t_eadsrctl *ctl  = (t_eadsrctl *)(w[1]);
    t_float attack  = ctl->c_attack;
    t_float decay   = ctl->c_decay;
    t_float sustain = ctl->c_sustain;
    t_float release = ctl->c_release;
    t_float state   = ctl->c_state;
    t_float target  = ctl->c_target;
    t_int n = (t_int)(w[2]);
t_int i;
for (i = 0; i < n; i++){
    if (target == 1.0f){
        /* attack */
        *out++ = state;
        state += attack*(1 - state);
        target = (state > ENVELOPE_MAX) ? sustain : 1.0f;
    }
    else if (target == 0.0f){
        /* release */
        *out++ = state;
        state -= release*state;
    }
    else{
        /* decay */
        *out++ = state;
        state -= decay*(state-sustain);
    }
    }
/* save state */
    ctl->c_state = IS_DENORMAL(state) ? 0 : state;
    ctl->c_target = target;
    return (w+4);
}
static void eadsr_dsp(t_eadsr *x, t_signal **sp)
{
    x->x_sr = sp[0]->s_sr;
    dsp_add(eadsr_perform, 3, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec);
}                                  
void eadsr_free(void)
{
}
t_class *eadsr_class;
void *eadsr_new(t_floatarg attack, t_floatarg decay, t_floatarg sustain, t_floatarg release)
{
    t_eadsr *x = (t_eadsr *)pd_new(eadsr_class);
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("attack"));  
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("decay"));  
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("sustain"));     
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("release"));
    outlet_new(&x->x_obj, gensym("signal"));
x->x_ctl.c_state  = 0;
    x->x_ctl.c_target = 0;
    eadsr_attack(x, attack);
    eadsr_decay(x, decay);
    eadsr_sustain(x, sustain);
    eadsr_release(x, release);
return (void *)x;
}
void eadsr_tilde_setup(void)
{
    //post("eadsr~ v0.1");
    eadsr_class = class_new(gensym("eadsr~"), (t_newmethod)eadsr_new,
    	(t_method)eadsr_free, sizeof(t_eadsr), 0,  A_DEFFLOAT, A_DEFFLOAT, A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(eadsr_class, (t_method)eadsr_float, gensym("float"), A_FLOAT, 0);
    class_addmethod(eadsr_class, (t_method)eadsr_start, gensym("start"), 0);
    class_addmethod(eadsr_class, (t_method)eadsr_start, gensym("bang"), 0);
    class_addmethod(eadsr_class, (t_method)eadsr_stop, gensym("stop"), 0);
    class_addmethod(eadsr_class, (t_method)eadsr_dsp, gensym("dsp"), 0); 
    class_addmethod(eadsr_class, (t_method)eadsr_attack, gensym("attack"), A_FLOAT, 0);
    class_addmethod(eadsr_class, (t_method)eadsr_decay, gensym("decay"), A_FLOAT, 0);
    class_addmethod(eadsr_class, (t_method)eadsr_sustain, gensym("sustain"), A_FLOAT, 0);
    class_addmethod(eadsr_class, (t_method)eadsr_release, gensym("release"), A_FLOAT, 0);
}
--- NEW FILE: statwav~.c ---
/*
 *   dynwav.c  - static wavetable oscillator (scale + tabread)
 *   data organization is in (real, imag) pairs
 *   the first 2 components are (DC, NY)
 *   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>
#include <string.h>
static t_class *statwav_tilde_class;
typedef struct _statwav_tilde
{
    t_object x_obj;
    int x_npoints;
    float *x_vec;
    t_symbol *x_arrayname;
    float x_f;
} t_statwav_tilde;
static void *statwav_tilde_new(t_symbol *s)
{
    t_statwav_tilde *x = (t_statwav_tilde *)pd_new(statwav_tilde_class);
    x->x_arrayname = s;
    x->x_vec = 0;
    outlet_new(&x->x_obj, gensym("signal"));
    x->x_f = 0;
    return (x);
}
static t_int *statwav_tilde_perform(t_int *w)
{
    t_statwav_tilde *x = (t_statwav_tilde *)(w[1]);
    t_float *in = (t_float *)(w[2]);
    t_float *out = (t_float *)(w[3]);
    int n = (int)(w[4]);
    float maxindex;
    int imaxindex;
    float *buf = x->x_vec, *fp;
    int i;
maxindex = x->x_npoints;
    imaxindex = x->x_npoints;
if (!buf) goto zero;
for (i = 0; i < n; i++)
    {
      float phase = *in++;
      float modphase = phase - (int)phase;
      float findex;
      int index;
      int ia, ib, ic, id;
      float frac,  a,  b,  c,  d, cminusb;
      static int count;
if (modphase < 0.0f) modphase += 1.0f;
      findex = modphase * maxindex;
      index = findex;
frac = findex - index;
      ia = (imaxindex+index-1) % imaxindex;
      ib = index;
      ic = (index+1) % imaxindex;
      id = (index+2) % imaxindex;
a = buf[ia];
      b = buf[ib];
      c = buf[ic];
      d = buf[id];
      /* if (!i && !(count++ & 1023))
     post("fp = %lx,  shit = %lx,  b = %f",  fp, buf->b_shit,  b); */
      cminusb = c-b;
*out++ = b + frac * (
            cminusb - 0.5f * (frac-1.) * (
                (a - d + 3.0f * cminusb) * frac + (b - a - cminusb)
            )
        );                                                         
    }
    
    return (w+5);
 zero:
    while (n--) *out++ = 0;
return (w+5);
}
static void statwav_tilde_set(t_statwav_tilde *x, t_symbol *s)
{
    t_garray *a;
x->x_arrayname = s;
    if (!(a = (t_garray *)pd_findbyclass(x->x_arrayname, garray_class)))
    {
        if (*s->s_name)
            error("statwav~: %s: no such array", x->x_arrayname->s_name);
        x->x_vec = 0;
    }
    else if (!garray_getfloatarray(a, &x->x_npoints, &x->x_vec))
    {
        error("%s: bad template for statwav~", x->x_arrayname->s_name);
        x->x_vec = 0;
    }
    else garray_usedindsp(a);
}
static void statwav_tilde_dsp(t_statwav_tilde *x, t_signal **sp)
{
    statwav_tilde_set(x, x->x_arrayname);
dsp_add(statwav_tilde_perform, 4, x,
        sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
}
static void statwav_tilde_free(t_statwav_tilde *x)
{
}
void statwav_tilde_setup(void)
{
  //post("statwav~ v0.1");
    statwav_tilde_class = class_new(gensym("statwav~"),
        (t_newmethod)statwav_tilde_new, (t_method)statwav_tilde_free,
        sizeof(t_statwav_tilde), 0, A_DEFSYM, 0);
    CLASS_MAINSIGNALIN(statwav_tilde_class, t_statwav_tilde, x_f);
    class_addmethod(statwav_tilde_class, (t_method)statwav_tilde_dsp,
        gensym("dsp"), 0);
    class_addmethod(statwav_tilde_class, (t_method)statwav_tilde_set,
        gensym("set"), A_SYMBOL, 0);
}
--- diag.c DELETED ---
--- qmult.c DELETED ---
--- NEW FILE: bitsplit~.c ---
/*
 *   bitsplit.c  - convert a signal to a binary vector 
 *   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 <stdlib.h>
#include <stdio.h>
#include <math.h>
#define MAXCHANNELS 24
typedef struct bitsplitctl
{
    t_int c_outputs;
    t_float *c_input;
    t_float **c_output;
} t_bitsplitctl;
typedef struct bitsplit
{
  t_object x_obj;
  t_float x_f;
  t_bitsplitctl x_ctl;
} t_bitsplit;
static t_int *bitsplit_perform(t_int *word)
{
t_bitsplitctl *ctl  = (t_bitsplitctl *)(word[1]);
    t_int n             = (t_int)(word[2]);
    t_float *in         = ctl->c_input;
    t_int outputs       = ctl->c_outputs;
    t_float **out       = ctl->c_output;
    t_int i,j;
for (i=0;i<n;i++){
    long word = (in[i] * (float)(0x7fffffff));
    for (j=0; j<outputs; j++){
        out[j][i] = (float)((word >> 31) & 1);
        word <<= 1;
    }
    }
return (word+3);
}
static void bitsplit_dsp(t_bitsplit *x, t_signal **sp)
{
int i;
    x->x_ctl.c_input = sp[0]->s_vec;
    for (i=0;i<x->x_ctl.c_outputs;i++){
    x->x_ctl.c_output[i] = sp[i+1]->s_vec;
    }
    dsp_add(bitsplit_perform, 2, &x->x_ctl, sp[0]->s_n);
}
static void bitsplit_free(t_bitsplit *x)
{
    free (x->x_ctl.c_output);
}
t_class *bitsplit_class;
static void *bitsplit_new(t_floatarg channels)
{
    int i = (int)channels;
    t_bitsplit *x = (t_bitsplit *)pd_new(bitsplit_class);
if (i<1) i = 1;
    if (i>MAXCHANNELS) i = MAXCHANNELS;
    x->x_ctl.c_outputs = i;
    x->x_ctl.c_output = malloc(sizeof(float)*i);
while (i--) outlet_new(&x->x_obj, gensym("signal"));
return (void *)x;
}
void bitsplit_tilde_setup(void)
{
    bitsplit_class = class_new(gensym("bitsplit~"), (t_newmethod)bitsplit_new,
    	(t_method)bitsplit_free, sizeof(t_bitsplit), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(bitsplit_class, t_bitsplit, x_f); 
    class_addmethod(bitsplit_class, (t_method)bitsplit_dsp, gensym("dsp"), 0); 
}
--- dwt.c DELETED ---
--- ramp.c DELETED ---
--- NEW FILE: xfm~.c ---
/*
 *   xfm.c  -  cross frequency modulation object
 *   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.
 */
/*  
coupled fm. osc state equations:
phasor for system i =
[  1  -phi ]
(1+phi^2)^(1/2)  *   [ phi   1  ]
with phi = 2*pi*(freq_base + freq_mod * out_other) / sr
ideal phasor would be
[ cos(phi)  - sin(phi) ]
[ sin(phi)    cos(phi) ]
this means frequencies are warped:
2*pi*f_real = atan(2*pi*f)
some (possible) enhancements:
    + add an integrator to get phase modulation
    + undo the frequency warping
*/
#include "m_pd.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define SINSAMPLES 512
#define MYPI 3.1415927
#define DISTORTED   0
#define NORMALIZED  1
typedef struct xfmctl
{
    //t_float c_sintab[SINSAMPLES + 1];
    t_float c_x1, c_y1; /* state osc 1 */
    t_float c_x2, c_y2; /* state osc 2 */
    t_int c_type; /* type of algo */
} t_xfmctl;
typedef struct xfm
{
    t_object x_obj;
    t_float x_f;
    t_xfmctl x_ctl;
} t_xfm;
void xfm_type(t_xfm *x, t_float f)
{
    int t = (int)f;
if (t == DISTORTED) x->x_ctl.c_type = t;
    if (t == NORMALIZED) x->x_ctl.c_type = t;
}
static inline t_float xfm_sat(t_float x)
{
  const float max =  1;
  const float min = -1;
x = (x > max) ? (max) : (x);
  x = (x < min) ? (min) : (x);
return(x);
}
static t_int *xfm_perform(t_int *w)
{
t_float *inA     = (float *)(w[3]);
  t_float *inB     = (float *)(w[4]);
  t_float *fbA     = (float *)(w[5]);
  t_float *fbB     = (float *)(w[6]);
  t_float *outA    = (float *)(w[7]);
  t_float *outB    = (float *)(w[8]);
  t_xfmctl *ctl    = (t_xfmctl *)(w[1]);
  t_int n          = (t_int)(w[2]);
  //t_float *tab     = ctl->c_sintab;
t_float x1 = ctl->c_x1, y1 = ctl->c_y1,   z1, dx1, dy1, inv_norm1;
  t_float x2 = ctl->c_x2, y2 = ctl->c_y2,   z2, dx2, dy2, inv_norm2;
t_float scale = 2 * M_PI / sys_getsr();
t_int i;
switch(ctl->c_type){
  default:
  case  DISTORTED:
/* this is a 4 degree of freedom hyperchaotic system */
      /* two coupled saturated unstable oscillators */
for (i=0; i<n; i++){
      /* osc 1 */
      z1 = scale * (x2 * (*fbA++) + (*inA++));
      
      dx1 = x1 - z1*y1;
      dy1 = y1 + z1*x1;
      
      x1 = xfm_sat(dx1);
      y1 = xfm_sat(dy1);
      
      /* osc 2*/
      z2 = scale * (x1 * (*fbB++) + (*inB++));
          
      dx2 = x2 - z2*y2;
      dy2 = y2 + z2*x2;
      
      x2 = xfm_sat(dx2);
      y2 = xfm_sat(dy2);
      
      /* output */
      (*outA++) = x1;
      (*outB++) = x2;
      
      }
      break;
case NORMALIZED:
/* this is a an effective 2 degree of freedom quasiperiodic system */
      /* two coupled stable oscillators */
for (i=0; i<n; i++){
/* osc 1 */
      z1 = scale * (x2 * (*fbA++) + (*inA++));
      
      dx1 = x1 - z1*y1;
      dy1 = y1 + z1*x1;
      inv_norm1 = 1.0f / hypot(dx1, dy1);
      
          
      /* osc 2*/
      z2 = scale * (x1 * (*fbB++) + (*inB++));
      
      dx2 = x2 - z2*y2;
      dy2 = y2 + z2*x2;
      inv_norm2 = 1.0f / hypot(dx2, dy2);
      
      /* renormalize */
      x1 = dx1 * inv_norm1;
      y1 = dy1 * inv_norm1;
      x2 = dx2 * inv_norm2;
      y2 = dy2 * inv_norm2;
      
      /* output */
      (*outA++) = x1;
      (*outB++) = x2;
      }
      break;
  }
ctl->c_x1 = x1;
  ctl->c_y1 = y1;
  ctl->c_x2 = x2;
  ctl->c_y2 = y2;
return (w+9);
}
static void xfm_dsp(t_xfm *x, t_signal **sp)
{
  int n = sp[0]->s_n;
  int k;
dsp_add(xfm_perform, 
      8, 
      &x->x_ctl, 
      sp[0]->s_n, 
      sp[0]->s_vec, 
      sp[1]->s_vec, 
      sp[2]->s_vec,
      sp[3]->s_vec,
      sp[4]->s_vec,
      sp[5]->s_vec);
}                                  
static void xfm_free(t_xfm *x)
{
}
static void xfm_reset(t_xfm *x)
{
    x->x_ctl.c_x1 = 1;
    x->x_ctl.c_y1 = 0;
    x->x_ctl.c_x2 = 1;
    x->x_ctl.c_y2 = 0;
}
t_class *xfm_class;
static void *xfm_new(t_floatarg algotype)
{
    t_xfm *x = (t_xfm *)pd_new(xfm_class);
/* ins */
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));  
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));  
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));
/* outs */
    outlet_new(&x->x_obj, gensym("signal")); 
    outlet_new(&x->x_obj, gensym("signal"));
/* init data */
    xfm_reset(x);
    xfm_type(x, algotype);
return (void *)x;
}
void xfm_tilde_setup(void)
{
    //post("xfm~ v0.1");
    xfm_class = class_new(gensym("xfm~"), (t_newmethod)xfm_new,
    	(t_method)xfm_free, sizeof(t_xfm), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(xfm_class, t_xfm, x_f);
    class_addmethod(xfm_class, (t_method)xfm_type, gensym("type"), A_FLOAT, 0);
    class_addmethod(xfm_class, (t_method)xfm_dsp, gensym("dsp"), 0); 
    class_addmethod(xfm_class, (t_method)xfm_reset, gensym("reset"), 0);
}
--- abs.c DELETED ---
--- ead.c DELETED ---
--- eadsr.c DELETED ---
--- fdn.c DELETED ---
--- NEW FILE: dwt~.c ---
/*
 *   dwt.c  - code for discrete wavelet transform 
 *   (symmetric interpolating biorthogonal wavelets using the lifting transform) 
 *   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>
#include <string.h>
#define MAXORDER 64
typedef enum{DWT,IDWT,DWT16,IDWT16} t_dwttype;
typedef struct dwtctl
{
  t_float c_update[MAXORDER];
  t_float c_predict[MAXORDER];
  t_int c_nupdate;
  t_int c_npredict;
  t_int c_levels;
  t_int c_fakein;
  t_float c_fakeval;
  t_int c_mask;
  char c_name[16];
  t_int *c_clutter;
  t_int *c_unclutter;
  t_int c_permute;
  t_dwttype c_type;
} t_dwtctl;
typedef struct dwt
{
  t_object x_obj;
  t_float x_f;
  t_dwtctl x_ctl;
} t_dwt;
static void dwt_even(t_dwt *x, t_floatarg f)
{
  int k = (int)f;
  int i, j;
  float *p = x->x_ctl.c_predict;
  float *u = x->x_ctl.c_update;
  float  l, xi, xj;
if ((k>0) && (k<MAXORDER/2))
    {
      for (i=0; i<k; i++){
    l = 1;
    xi = 1+2*i;
    for (j=0; j<k; j++)
      {
        xj = 1+2*j;
        if (i != j)  l /= (1 - ((xi*xi) / (xj*xj)));
      }
    l *= .5;
p[k-i-1] = l;
    p[k+i] = l;
    u[k-i-1] = l/2;
    u[k+i] = l/2;
      }
x->x_ctl.c_npredict = 2*k;
      x->x_ctl.c_nupdate  = 2*k;
    }
}
static void dwt_wavelet(t_dwt *x, t_floatarg f)
{
  int k = (int)f;
  t_float *p = x->x_ctl.c_predict;
  t_float *u = x->x_ctl.c_update;
  t_int *np = &x->x_ctl.c_npredict;
  t_int *nu = &x->x_ctl.c_nupdate;
switch(k)
    {
    default:
    case 1: /* haar */
      *np = *nu = 2; /* actual order is one */
      p[0] = 1;
      p[1] = 0;
      u[0] = 0;
      u[1] = .5;
      break;
case 2: /* hat */
    case 3:
      *np = *nu = 2;
      p[0] = .5;
      p[1] = .5;
      u[0] = .25;
      u[1] = .25;
      break;
case 4: /* N = 4, N~ = 4 */
    case 5:
      *np = *nu = 4;
      p[0] = -0.0625;
      p[1] =  0.5625;
      p[2] =  0.5625;
      p[3] = -0.0625;
      u[0] = -0.03125;
      u[1] =  0.28125;
      u[2] =  0.28125;
      u[3] = -0.03125;
      break;
case 6:
    case 7:
      *np = *nu = 6;
      p[0] =    0.01171875000000;
      p[1] =   -0.09765625000000;
      p[2] =    0.58593750000000;
      p[3] =    0.58593750000000;
      p[4] =   -0.09765625000000;
      p[5] =    0.01171875000000;
      u[0] =    0.00585937500000;
      u[1] =   -0.04882812500000;
      u[2] =    0.29296875000000;
      u[3] =    0.29296875000000;
      u[4] =   -0.04882812500000;
      u[5] =    0.00585937500000;
      break;
    }
}
static inline void dwt_perform_permutation(t_float *S, int n, t_int *f)
{
  t_int k,l;
  t_float swap;
  for(k=0; k<n; k++)
    {
      l = f[k];
      while (l<k) l = f[l];
      swap = S[k];
      S[k] = S[l];
      S[l] = swap;
    }
}
static void dwt_permutation(t_dwt *x, t_int n){
t_dwtctl *ctl = &x->x_ctl;
  t_int k, L=0, l, start, power;
  t_int nsave = n;
while(nsave>>=1) L++;
if (ctl->c_clutter)   free(ctl->c_clutter);
  if (ctl->c_unclutter) free(ctl->c_unclutter);
ctl->c_clutter = (t_int *)malloc(n*sizeof(t_int));
  ctl->c_unclutter = (t_int *)malloc(n*sizeof(t_int));
for(l = L, start = n/2, power=1; l>0; l--, start /=2, power *=2)
    {
      for(k=0; k<start; k++)
    {
      ctl->c_unclutter[start+k] = (1 + 2*k) * power;
    }
    }
  ctl->c_unclutter[0] = 0;
for(k=0; k<n; k++)
    ctl->c_clutter[ctl->c_unclutter[k]] = k;
return;
/* debug */
  for(k=0; k<n; k++)
    printf("clutter[%ld] = %ld\n", k, ctl->c_clutter[k]);
  for(k=0; k<n; k++)
    printf("unclutter[%ld] = %ld\n", k, ctl->c_unclutter[k]);
exit(1);
}
static void idwt_coef(t_dwt *x, t_floatarg index, t_floatarg value)
{
  x->x_ctl.c_fakein = (int)index;
  x->x_ctl.c_fakeval = value;
}
static void dwt_print(t_dwt *x)
{
  int i;
printf("%s: predict: [ ", x->x_ctl.c_name);
  for (i=0; i<x->x_ctl.c_npredict; i++) printf("%f ", x->x_ctl.c_predict[i]);
  printf("], ");
printf("update: [ ");
  for (i=0; i<x->x_ctl.c_nupdate; i++) printf("%f ", x->x_ctl.c_update[i]);
  printf("]\n");
}
static void dwt_filter(t_dwt *x,  t_symbol *s, int argc, t_atom *argv)
{
  int invalid_argument = 0;
  int i;
char *name = x->x_ctl.c_name;
float *pfilter = x->x_ctl.c_predict; 
  float *ufilter = x->x_ctl.c_update; 
  float *mask = NULL;
t_int *length = NULL;
  float sum = 0;
if (s == gensym("predict"))
    {
      mask = pfilter;
      length = &(x->x_ctl.c_npredict);
    }
  else if (s == gensym("update"))
    {
      mask = ufilter;
      length = &(x->x_ctl.c_nupdate);
    }
  else if (s == gensym("mask"))
    {
      mask = NULL;
    }
  else
    {
      return;
    }
if (argc >= MAXORDER) post("%s: error, maximum order exceeded.",name);
  else if ((x->x_ctl.c_type == DWT16 || x->x_ctl.c_type == IDWT16 ) && (argc != 16))
    post("%s: error, need to have 16 coefficients.",name);
  else if (argc == 0) post("%s: no arguments given.",name);
  else if (argc & 1) post("%s: error, only an even number of coefficients is allowed.", name);
  else
    {
      for (i=0; i<argc; i++){
    if (argv[i].a_type != A_FLOAT )
      {
        invalid_argument = 1;
        break;
      }
      }
if (invalid_argument) post("%s: invalid argument, must be a number.", name);
      else
    {
      if (mask) /* one of update / predict */
        {
          for (i=0; i<argc; i++) mask[i] = argv[i].a_w.w_float;
          *length = argc;
        }
      else /* both + normalization */
        {
          for (i=0; i<argc; i++) sum +=  argv[i].a_w.w_float;
          for (i=0; i<argc; i++)
    	{
    	  pfilter[i] =  argv[i].a_w.w_float / sum;
    	  ufilter[i] =  argv[i].a_w.w_float / (sum*2);
    	}
          x->x_ctl.c_npredict = argc;
          x->x_ctl.c_nupdate = argc;
        }
    }
}
}
static inline void dwtloop(float *vector, 
    		   int source,
    		   int dest,
    	     int increment,
    	     int backup,
    	     int numcoef, 
    	     int mask,
    	     float *filter,
    	     int filtlength,
    	     float sign)
{
int k,m;
  float acc;
for (k = 0; k < numcoef; k++)
    {
      acc = 0;
      for (m = 0; m < filtlength; m++)
    {
acc += filter[m] * vector[source];
      source += increment;
      source &= mask;
    }
      vector[dest] += sign * acc;
      dest += increment;
      source -= backup;
      source &= mask;
    }
}
static inline void dwtloop16(float *vector, 
    	     int source,
    	     int dest,
    	     int increment,
    	     int backup,
    	     int numcoef, 
    	     int mask,
    	     float *filter,
    	     int filtlength, /* ignored, set to 16 */
    	     float sign)
{
int k,m;
  float acc;
for (k = 0; k < numcoef; k++)
    {
      acc = 0;
acc += filter[0] * vector[source];
      source += increment;
      source &= mask;
acc += filter[1] * vector[source];
      source += increment;
      source &= mask;
acc += filter[2] * vector[source];
      source += increment;
      source &= mask;
acc += filter[3] * vector[source];
      source += increment;
      source &= mask;
acc += filter[4] * vector[source];
      source += increment;
      source &= mask;
acc += filter[5] * vector[source];
      source += increment;
      source &= mask;
acc += filter[6] * vector[source];
      source += increment;
      source &= mask;
acc += filter[7] * vector[source];
      source += increment;
      source &= mask;
acc += filter[8] * vector[source];
      source += increment;
      source &= mask;
acc += filter[9] * vector[source];
      source += increment;
      source &= mask;
acc += filter[10] * vector[source];
      source += increment;
      source &= mask;
acc += filter[11] * vector[source];
      source += increment;
      source &= mask;
acc += filter[12] * vector[source];
      source += increment;
      source &= mask;
acc += filter[13] * vector[source];
      source += increment;
      source &= mask;
acc += filter[14] * vector[source];
      source += increment;
      source &= mask;
acc += filter[15] * vector[source];
      source += increment;
      source &= mask;
vector[dest] += sign * acc;
      dest += increment;
      source -= backup;
      source &= mask;
    }
}
static t_int *dwt_perform(t_int *w)
{
t_float *in     = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_dwtctl *ctl  = (t_dwtctl *)(w[1]);
t_int n = (t_int)(w[2]);
int i;
int numcoef = n/2;
  /*  int source_u = ((1 - ctl->c_nupdate)/2 - 1);
   *  int source_p = ((1 - ctl->c_npredict)/2);
   */
  int source_u = ((2 - ctl->c_nupdate) - 1);
  int source_p = ((2 - ctl->c_npredict));
  int increment = 2;
  int dest = 1;
  int backup_u = (ctl->c_nupdate-1)*2;
  int backup_p = (ctl->c_npredict-1)*2;
/* copy input to output */
  if (in != out)
    for (i=0; i<n; i++) out[i]=in[i];
/* fake input */
  /*    for (i=0; i<n; i++) out[i]=0; out[n/8]=1;*/
/* backward transform */
/* iterate over all levels */
  for (i=0; i < ctl->c_levels; i++){
/* foreward predict */
    dwtloop(out, (source_p & (n-1)), dest, increment, backup_p, numcoef, n-1, ctl->c_predict, ctl->c_npredict, -1);
/* foreward update */
    dwtloop(out, (source_u & (n-1)), 0,    increment, backup_u, numcoef, n-1, ctl->c_update,  ctl->c_nupdate,  +1);
/* update control parameters */
    numcoef /= 2;
    source_p *= 2;
    source_u *= 2;
    backup_p *= 2;
    backup_u *= 2;
    increment *= 2;
    dest *= 2;
  }
if (ctl->c_permute) 
    dwt_perform_permutation(out, n, ctl->c_unclutter);
return (w+5);
}
static t_int *idwt_perform(t_int *w)
{
t_float *in     = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_dwtctl *ctl  = (t_dwtctl *)(w[1]);
t_int n = (t_int)(w[2]);
int i;
int numcoef = 1;
  int source_u = ((2 - ctl->c_nupdate) - 1) * (n/2);
  int source_p = ((2 - ctl->c_npredict)) * (n/2);
  int increment = n;
  int dest = n/2;
  int backup_u = (ctl->c_nupdate-1)*n;
  int backup_p = (ctl->c_npredict-1)*n;
  int fake_in = ctl->c_fakein;
  float fake_val = ctl->c_fakeval;
/* copy input to output */
  if (in != out)
    for (i=0; i<n; i++) out[i]=in[i];
/* fake input */
if ((fake_in >= 0) && (fake_in<n)){
    for (i=0; i<n; i++) out[i]=0; 
    out[fake_in]=fake_val;
  }
if (ctl->c_permute) 
    dwt_perform_permutation(out, n, ctl->c_clutter);
/* backward transform */
/* iterate over all levels */
  for (i=0; i < ctl->c_levels; i++){
/* backward update */
    dwtloop(out, (source_u & (n-1)), 0,    increment, backup_u, numcoef, n-1, ctl->c_update,  ctl->c_nupdate,  -1);
/* backward predict */
    dwtloop(out, (source_p & (n-1)), dest, increment, backup_p, numcoef, n-1, ctl->c_predict, ctl->c_npredict, +1);
/* update control parameters */
    numcoef *= 2;
    source_p /= 2;
    source_u /= 2;
    backup_p /= 2;
    backup_u /= 2;
    increment /= 2;
    dest /= 2;
  }
return (w+5);
}
static t_int *dwt16_perform(t_int *w)
{
t_float *in     = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_dwtctl *ctl  = (t_dwtctl *)(w[1]);
t_int n = (t_int)(w[2]);
int i;
int numcoef = n/2;
  /*  int source_u = ((1 - ctl->c_nupdate)/2 - 1);
   *  int source_p = ((1 - ctl->c_npredict)/2);
   */
  int source_u = ((2 - ctl->c_nupdate) - 1);
  int source_p = ((2 - ctl->c_npredict));
  int increment = 2;
  int dest = 1;
  int backup_u = (ctl->c_nupdate-1)*2;
  int backup_p = (ctl->c_npredict-1)*2;
/* copy input to output */
  if (in != out)
    for (i=0; i<n; i++) out[i]=in[i];
/* fake input */
  /*    for (i=0; i<n; i++) out[i]=0; out[n/8]=1;*/
/* backward transform */
/* iterate over all levels */
  for (i=0; i < ctl->c_levels; i++){
/* foreward predict */
    dwtloop16(out, (source_p & (n-1)), dest, increment, backup_p, numcoef, n-1, ctl->c_predict, 16, -1);
/* foreward update */
    dwtloop16(out, (source_u & (n-1)), 0,    increment, backup_u, numcoef, n-1, ctl->c_update,  16,  +1);
/* update control parameters */
    numcoef /= 2;
    source_p *= 2;
    source_u *= 2;
    backup_p *= 2;
    backup_u *= 2;
    increment *= 2;
    dest *= 2;
  }
if (ctl->c_permute) 
    dwt_perform_permutation(out, n, ctl->c_unclutter);
return (w+5);
}
static t_int *idwt16_perform(t_int *w)
{
t_float *in     = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_dwtctl *ctl  = (t_dwtctl *)(w[1]);
t_int n = (t_int)(w[2]);
int i;
int numcoef = 1;
  int source_u = ((2 - ctl->c_nupdate) - 1) * (n/2);
  int source_p = ((2 - ctl->c_npredict)) * (n/2);
  int increment = n;
  int dest = n/2;
  int backup_u = (ctl->c_nupdate-1)*n;
  int backup_p = (ctl->c_npredict-1)*n;
  int fake_in = ctl->c_fakein;
  float fake_val = ctl->c_fakeval;
/* copy input to output */
  if (in != out)
    for (i=0; i<n; i++) out[i]=in[i];
/* fake input */
if ((fake_in >= 0) && (fake_in<n)){
    for (i=0; i<n; i++) out[i]=0; 
    out[fake_in]=fake_val;
  }
if (ctl->c_permute) 
    dwt_perform_permutation(out, n, ctl->c_clutter);
/* backward transform */
/* iterate over all levels */
  for (i=0; i < ctl->c_levels; i++){
/* backward update */
    dwtloop16(out, (source_u & (n-1)), 0,    increment, backup_u, numcoef, n-1, ctl->c_update,  16,  -1);
/* backward predict */
    dwtloop16(out, (source_p & (n-1)), dest, increment, backup_p, numcoef, n-1, ctl->c_predict, 16, +1);
/* update control parameters */
    numcoef *= 2;
    source_p /= 2;
    source_u /= 2;
    backup_p /= 2;
    backup_u /= 2;
    increment /= 2;
    dest /= 2;
  }
return (w+5);
}
static void dwt_dsp(t_dwt *x, t_signal **sp)
{
int n = sp[0]->s_n;
  int ln = 0;
dwt_permutation(x, n);
x->x_ctl.c_mask = n-1;
  while (n >>= 1) ln++;     
  x->x_ctl.c_levels = ln;
switch(x->x_ctl.c_type){
  case DWT:
    dsp_add(dwt_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
    break;
  case IDWT:
    dsp_add(idwt_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
    break;
  case DWT16:
    dsp_add(dwt16_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
    break;
  case IDWT16:
    dsp_add(idwt16_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
    break;
}
}
void dwt_free(t_dwt *x)
{
if (x->x_ctl.c_clutter)   free(x->x_ctl.c_clutter);
  if (x->x_ctl.c_unclutter) free(x->x_ctl.c_unclutter);
}
t_class *dwt_class, *idwt_class, *dwt16_class, *idwt16_class;
static void dwt_reset(t_dwt *x)
{
  bzero(x->x_ctl.c_update, 16*sizeof(t_float));
  bzero(x->x_ctl.c_predict, 16*sizeof(t_float));
x->x_ctl.c_update[7] = .25;
    x->x_ctl.c_update[8] = .25;
    x->x_ctl.c_nupdate = 16;
x->x_ctl.c_predict[7] = .5;
    x->x_ctl.c_predict[8] = .5;
    x->x_ctl.c_npredict = 16;
x->x_ctl.c_fakein = -1;
    x->x_ctl.c_fakeval = 0;
}
static void *dwt_new_common(t_floatarg permute)
{
    t_dwt *x = (t_dwt *)pd_new(dwt_class);
    int i;
outlet_new(&x->x_obj, gensym("signal"));
/* init data */
    dwt_reset(x);
x->x_ctl.c_clutter = NULL;
    x->x_ctl.c_unclutter = NULL;
    x->x_ctl.c_permute = (t_int) permute;
return (void *)x;
}
static void *dwt_new(t_floatarg permute)
{
    t_dwt *x = dwt_new_common(permute);
    sprintf(x->x_ctl.c_name,"dwt");
    x->x_ctl.c_type = DWT;
    return (void *)x;
}
static void *idwt_new(t_floatarg permute)
{
    t_dwt *x = dwt_new_common(permute);
    sprintf(x->x_ctl.c_name,"idwt");
    x->x_ctl.c_type = IDWT;
    return (void *)x;
}
static void *dwt16_new(t_floatarg permute)
{
    t_dwt *x = dwt_new_common(permute);
    sprintf(x->x_ctl.c_name,"dwt16");
    x->x_ctl.c_type = DWT16;
    return (void *)x;
}
static void *idwt16_new(t_floatarg permute)
{
    t_dwt *x = dwt_new_common(permute);
    sprintf(x->x_ctl.c_name,"idwt16");
    x->x_ctl.c_type = IDWT16;
    return (void *)x;
}
void dwt_tilde_setup(void)
{
  //post("dwt~ v0.1");
dwt_class = class_new(gensym("dwt~"), (t_newmethod)dwt_new,
    	(t_method)dwt_free, sizeof(t_dwt), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(dwt_class, t_dwt, x_f);
    class_addmethod(dwt_class, (t_method)dwt_print, gensym("print"), 0);
    class_addmethod(dwt_class, (t_method)dwt_reset, gensym("reset"), 0);
    class_addmethod(dwt_class, (t_method)dwt_dsp, gensym("dsp"), 0);
class_addmethod(dwt_class, (t_method)dwt_filter, gensym("predict"), A_GIMME, 0); 
    class_addmethod(dwt_class, (t_method)dwt_filter, gensym("update"), A_GIMME, 0); 
    class_addmethod(dwt_class, (t_method)dwt_filter, gensym("mask"), A_GIMME, 0);
class_addmethod(dwt_class, (t_method)dwt_even, gensym("even"), A_DEFFLOAT, 0);
    class_addmethod(dwt_class, (t_method)idwt_coef, gensym("coef"), A_DEFFLOAT, A_DEFFLOAT, 0);
/*class_addmethod(dwt_class, (t_method)dwt_wavelet, gensym("wavelet"), A_DEFFLOAT, 0); */
idwt_class = class_new(gensym("idwt~"), (t_newmethod)idwt_new,
    	(t_method)dwt_free, sizeof(t_dwt), 0, A_DEFFLOAT, 0);
CLASS_MAINSIGNALIN(idwt_class, t_dwt, x_f);
    class_addmethod(idwt_class, (t_method)dwt_print, gensym("print"), 0);
    class_addmethod(idwt_class, (t_method)dwt_dsp, gensym("dsp"), 0);
class_addmethod(idwt_class, (t_method)dwt_filter, gensym("predict"), A_GIMME, 0); 
    class_addmethod(idwt_class, (t_method)dwt_filter, gensym("update"), A_GIMME, 0); 
    class_addmethod(idwt_class, (t_method)dwt_filter, gensym("mask"), A_GIMME, 0);
class_addmethod(idwt_class, (t_method)idwt_coef, gensym("coef"), A_DEFFLOAT, A_DEFFLOAT, 0);
class_addmethod(idwt_class, (t_method)dwt_even, gensym("even"), A_DEFFLOAT, 0);
dwt16_class = class_new(gensym("dwt16~"), (t_newmethod)dwt16_new,
    	(t_method)dwt_free, sizeof(t_dwt), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(dwt16_class, t_dwt, x_f);
    class_addmethod(dwt16_class, (t_method)dwt_print, gensym("print"), 0);
    class_addmethod(dwt16_class, (t_method)dwt_reset, gensym("reset"), 0);
    class_addmethod(dwt16_class, (t_method)dwt_dsp, gensym("dsp"), 0);
class_addmethod(dwt16_class, (t_method)dwt_filter, gensym("predict"), A_GIMME, 0); 
    class_addmethod(dwt16_class, (t_method)dwt_filter, gensym("update"), A_GIMME, 0); 
    class_addmethod(dwt16_class, (t_method)dwt_filter, gensym("mask"), A_GIMME, 0);
idwt16_class = class_new(gensym("idwt16~"), (t_newmethod)idwt16_new,
    	(t_method)dwt_free, sizeof(t_dwt), 0, A_DEFFLOAT, 0);
CLASS_MAINSIGNALIN(idwt16_class, t_dwt, x_f);
    class_addmethod(idwt16_class, (t_method)dwt_print, gensym("print"), 0);
    class_addmethod(idwt16_class, (t_method)dwt_dsp, gensym("dsp"), 0);
class_addmethod(idwt16_class, (t_method)dwt_filter, gensym("predict"), A_GIMME, 0); 
    class_addmethod(idwt16_class, (t_method)dwt_filter, gensym("update"), A_GIMME, 0); 
    class_addmethod(idwt16_class, (t_method)dwt_filter, gensym("mask"), A_GIMME, 0);
class_addmethod(idwt16_class, (t_method)idwt_coef, gensym("coef"), A_DEFFLOAT, A_DEFFLOAT, 0);
}
--- blocknorm.c DELETED ---
--- junction.c DELETED ---
--- cheby.c DELETED ---
--- matrix.c DELETED ---
--- NEW FILE: tabreadmix~.c ---
/*
 *   tabreadmix.c  -  an overlap add tabread~ clone
 *   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 "extlib_util.h"
/******************** tabreadmix~ ***********************/
static t_class *tabreadmix_tilde_class;
typedef struct _tabreadmix_tilde
{
    t_object x_obj;
    int x_npoints;
    float *x_vec;
    t_symbol *x_arrayname;
    float x_f;
/* file position vars */
    int x_currpos;
    int x_prevpos;
/* cross fader state vars */
    int x_xfade_size;
    int x_xfade_phase;
    float x_xfade_cos;
    float x_xfade_sin;
    float x_xfade_state_c;
    float x_xfade_state_s;
} t_tabreadmix_tilde;
inline void tabreadmix_tilde_wrapindices(t_tabreadmix_tilde *x)
{
    int max;
/* modulo */
    x->x_currpos %= x->x_npoints;
    x->x_prevpos %= x->x_npoints;
/* make sure 0<=..<x->x_npoints */
    //if (x->x_currpos < 0) x->x_currpos += x->x_npoints;
    //if (x->x_prevpos < 0) x->x_prevpos += x->x_npoints;
    x->x_currpos += (x->x_currpos < 0) * x->x_npoints;
    x->x_prevpos += (x->x_prevpos < 0) * x->x_npoints;
}
#define min(x,y) ((x)<(y)?(x):(y))
static t_int *tabreadmix_tilde_perform(t_int *w)
{
    t_tabreadmix_tilde *x = (t_tabreadmix_tilde *)(w[1]);
    t_float *pos = (t_float *)(w[2]);
    t_float *out = (t_float *)(w[3]);
    int n = (int)(w[4]);    
    int maxxindex;
    float *buf = x->x_vec;
    int i;
    float currgain, prevgain;
    float c,s;
    int chunk;
    int leftover;
    int newpos = (int)*pos;
maxxindex = x->x_npoints;
    if (!buf) goto zero;
    if (maxxindex <= 0) goto zero;
while (n){
/* process as much data as possible */
    leftover = x->x_xfade_size - x->x_xfade_phase;
    chunk = min(n, leftover);
for (i = 0; i < chunk; i++){
        /* compute crossfade gains from oscillator state */
        currgain = 0.5f - x->x_xfade_state_c;
        prevgain = 0.5f + x->x_xfade_state_c;
        
        /* check indices & wrap */
        tabreadmix_tilde_wrapindices(x);
/* mix and write */
        newpos = (int)(*pos++);
        *out++ = currgain * buf[x->x_currpos++] + prevgain * buf[x->x_prevpos++];
        
        /* advance oscillator */
        c =  x->x_xfade_state_c * x->x_xfade_cos -  x->x_xfade_state_s * x->x_xfade_sin;
        s =  x->x_xfade_state_c * x->x_xfade_sin +  x->x_xfade_state_s * x->x_xfade_cos;
        x->x_xfade_state_c = c;
        x->x_xfade_state_s = s;
    }
/* update indices */
    x->x_xfade_phase += chunk;
    n -= chunk;
    //pos += chunk;
/* check if prev chunk is finished */
    if (x->x_xfade_size == x->x_xfade_phase){
        x->x_prevpos = x->x_currpos;
        x->x_currpos = newpos;
        x->x_xfade_state_c = 0.5f;
        x->x_xfade_state_s = 0.0f;
        x->x_xfade_phase = 0;
    }
}
/* return if we ran out of data */
    return (w+5);
zero:
    while (n--) *out++ = 0;
    return (w+5);
}
static void tabreadmix_tilde_blocksize(t_tabreadmix_tilde *x, t_float size)
{
    double prev_phase;
    int max;
    float fmax = (float)x->x_npoints * 0.5f;
if (size < 1.0) size = 1.0;
prev_phase = (double)x->x_xfade_phase;
    prev_phase *= size;
    prev_phase /= (double)x->x_xfade_size;
/* preserve the crossfader state */
    x->x_xfade_phase = (int)prev_phase;
    x->x_xfade_size = (int)size;
x->x_xfade_cos = cos(M_PI / (float)x->x_xfade_size);
    x->x_xfade_sin = sin(M_PI / (float)x->x_xfade_size);
/* make sure indices are inside array */
    if (x->x_npoints == 0){
    x->x_currpos = 0;
    x->x_prevpos = 0;
    }
//else tabreadmix_tilde_wrapindices(x);
}
void tabreadmix_tilde_pitch(t_tabreadmix_tilde *x, t_float f)
{
    if (f < 1) f = 1;
tabreadmix_tilde_blocksize(x, sys_getsr() / f);
}
void tabreadmix_tilde_chunks(t_tabreadmix_tilde *x, t_float f)
{
    if (f < 1.0f) f = 1.0f;
    tabreadmix_tilde_blocksize(x, (float)x->x_npoints / f);
}
void tabreadmix_tilde_bang(t_tabreadmix_tilde *x, t_float f)
{
    //trigger a chunk reset on next dsp call
    x->x_xfade_phase = x->x_xfade_size;
}
void tabreadmix_tilde_set(t_tabreadmix_tilde *x, t_symbol *s)
{
    t_garray *a;
x->x_arrayname = s;
    if (!(a = (t_garray *)pd_findbyclass(x->x_arrayname, garray_class)))
    {
        if (*s->s_name)
            error("tabreadmix~: %s: no such array", x->x_arrayname->s_name);
        x->x_vec = 0;
    }
    else if (!garray_getfloatarray(a, &x->x_npoints, &x->x_vec))
    {
        error("%s: bad template for tabreadmix~", x->x_arrayname->s_name);
        x->x_vec = 0;
    }
    else garray_usedindsp(a);
/* make sure indices are inside array */
    if (x->x_npoints == 0){
    x->x_currpos = 0;
    x->x_prevpos = 0;
    }
//else tabreadmix_tilde_wrapindices(x);
}
static void tabreadmix_tilde_dsp(t_tabreadmix_tilde *x, t_signal **sp)
{
    tabreadmix_tilde_set(x, x->x_arrayname);
dsp_add(tabreadmix_tilde_perform, 4, x,
        sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
}
static void tabreadmix_tilde_free(t_tabreadmix_tilde *x)
{
}
static void *tabreadmix_tilde_new(t_symbol *s)
{
    t_tabreadmix_tilde *x = (t_tabreadmix_tilde *)pd_new(tabreadmix_tilde_class);
    x->x_arrayname = s;
    x->x_vec = 0;
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("blocksize"));  
    outlet_new(&x->x_obj, gensym("signal"));
    x->x_f = 0;
    x->x_xfade_phase = 0;
    x->x_xfade_size = 1024;
    x->x_currpos = 0;
    x->x_prevpos = 0;
    x->x_xfade_state_c = 0.5f;
    x->x_xfade_state_s = 0.0f;
    tabreadmix_tilde_blocksize(x, 1024);
    return (x);
}
void tabreadmix_tilde_setup(void)
{
    tabreadmix_tilde_class = class_new(gensym("tabreadmix~"),
        (t_newmethod)tabreadmix_tilde_new, (t_method)tabreadmix_tilde_free,
        sizeof(t_tabreadmix_tilde), 0, A_DEFSYM, 0);
    CLASS_MAINSIGNALIN(tabreadmix_tilde_class, t_tabreadmix_tilde, x_f);
    class_addmethod(tabreadmix_tilde_class, (t_method)tabreadmix_tilde_dsp,
        gensym("dsp"), 0);
    class_addmethod(tabreadmix_tilde_class, (t_method)tabreadmix_tilde_set,
        gensym("set"), A_SYMBOL, 0);
    class_addmethod(tabreadmix_tilde_class, (t_method)tabreadmix_tilde_blocksize,
        gensym("blocksize"), A_FLOAT, 0);
    class_addmethod(tabreadmix_tilde_class, (t_method)tabreadmix_tilde_pitch,
        gensym("pitch"), A_FLOAT, 0);
    class_addmethod(tabreadmix_tilde_class, (t_method)tabreadmix_tilde_chunks,
        gensym("chunks"), A_FLOAT, 0);
    class_addmethod(tabreadmix_tilde_class, (t_method)tabreadmix_tilde_bang,
        gensym("bang"), 0);
}
--- NEW FILE: lattice~.c ---
/*
 *   lattice.c  - a lattice filter for pd
 *   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>
#define MAXORDER 1024
#define MAXREFCO 0.9999f
typedef struct latticesegment
{
    t_float delay; // delay element
    t_float rc;    // reflection coefficient
} t_latticesegment;
typedef struct latticectl
{
    t_latticesegment c_segment[MAXORDER]; // array of lattice segment data
    t_int c_segments;
} t_latticectl;
typedef struct lattice
{
  t_object x_obj;
  t_float x_f;
  t_latticectl x_ctl;
} t_lattice;
static t_int *lattice_perform(t_int *w)
{
t_float *in    = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_latticectl *ctl  = (t_latticectl *)(w[1]);
  t_int i,j;
  t_int n = (t_int)(w[2]);
  t_float x,rc,d;
t_latticesegment* seg = ctl->c_segment;
  t_int segments = ctl->c_segments;
  for (i=0; i<n; i++)
  {
      x = *in++;
      // section 0
      rc = seg[0].rc;
      x += seg[0].delay * rc;
// section 1 to segments-1
      for (j=1; j < segments; j++)
      {
      rc = seg[j].rc;
      d  = seg[j].delay;
      x += d * rc;
      seg[j-1].delay = d - rc * x; 
      }
      // stub
      seg[segments-1].delay = x;
*out++ = x;
  }
return (w+5);
}
static void lattice_dsp(t_lattice *x, t_signal **sp)
{
    dsp_add(lattice_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
}                                  
static void lattice_free(void)
{
}
t_class *lattice_class;
static void lattice_rc(t_lattice *x, t_float segment, t_float refco)
{
    t_int seg = (t_float)segment;
    if ((seg >= 0) && (seg < x->x_ctl.c_segments)){
    if (refco >= MAXREFCO) refco = MAXREFCO;
    if (refco <= -MAXREFCO) refco = -MAXREFCO;
    x->x_ctl.c_segment[seg].rc = refco;
    }
}
static void lattice_reset(t_lattice *x)
{
    t_float* buf = (t_float *)x->x_ctl.c_segment;
    t_int n = x->x_ctl.c_segments;
    t_int i;
    for (i=0; i<n; i++) buf[i]=0;
}
static void *lattice_new(t_floatarg segments)
{
    t_lattice *x = (t_lattice *)pd_new(lattice_class);
    t_int seg = (t_int)segments;
outlet_new(&x->x_obj, gensym("signal"));
if (seg < 1) seg = 1;
    if (seg > MAXORDER) seg = MAXORDER;
x->x_ctl.c_segments = seg;
lattice_reset(x);
return (void *)x;
}
void lattice_tilde_setup(void)
{
  //post("lattice~ v0.1");
    lattice_class = class_new(gensym("lattice~"), (t_newmethod)lattice_new,
    	(t_method)lattice_free, sizeof(t_lattice), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(lattice_class, t_lattice, x_f); 
    class_addmethod(lattice_class, (t_method)lattice_dsp, gensym("dsp"), 0); 
    class_addmethod(lattice_class, (t_method)lattice_reset, gensym("reset"), 0); 
    class_addmethod(lattice_class, (t_method)lattice_rc, gensym("rc"), A_FLOAT, A_FLOAT, 0);
}
--- NEW FILE: junction~.c ---
/*
 *   junction.c  - computes a lossless circulant junction 
 *   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 <stdlib.h>
#include <stdio.h>
#include <math.h>
typedef struct junctionctl
{
  t_int c_channels;
  t_float **c_in;
  t_float **c_out;
  t_float *c_buffer;
  t_float *c_coef;
  t_float c_norm;
} t_junctionctl;
typedef struct junction
{
  t_object x_obj;
  t_float x_f;
  t_junctionctl x_ctl;
} t_junction;
void junction_bang(t_junction *x)
{
  int i, n       = x->x_ctl.c_channels;
  t_float *coef  = x->x_ctl.c_coef;
  t_float r;
for (i=1; i<n/2; i++)
    {
      r = rand();
      r *= ((2 * M_PI)/RAND_MAX);
      coef[i]= cos(r);
      coef[n-i] = sin(r);
    }
coef[0] = (rand() & 1) ? 1 : -1;
  coef[n/2] = (rand() & 1) ? 1 : -1;
/*  mayer_realfft(n, coef); */
}
void junction_random(t_junction *x, t_floatarg f)
{
  srand((int)f);
  junction_bang(x);
}
static t_int *junction_perform(t_int *w)
{
t_junctionctl *ctl  = (t_junctionctl *)(w[1]);
  t_int n             = (t_int)(w[2]);
  t_int i,j;
  t_float x,y;
t_int c             = ctl->c_channels;
  t_float **in        = ctl->c_in;
  t_float **out       = ctl->c_out;
  t_float *buf        = ctl->c_buffer;
  t_float *coef       = ctl->c_coef;
t_float norm        = ctl->c_norm;
for (i=0;i<n;i++)
    {
/* read input */
      for (j=0; j<c; j++)
    {
      buf[j] = in[j][i];
    }
/* transform */
      mayer_realfft(c, buf);
      for (j=1; j<c/2; j++)
    {
      float x,y,a,b;
      x = buf[j];
      y = buf[c-j];
      a = coef[j];
      b = coef[c-j];
      buf[j]   = a * x - b * y;
      buf[c-j] = a * y + b * x;
    }
      buf[0] *= coef[0];
      buf[c/2] *= coef[c/2];
mayer_realifft(c, buf);
/* write output */
      for (j=0; j<c; j++) 
    { 
      out[j][i] = buf[j] * norm;
    }  
    }
return (w+3);
}
static void junction_dsp(t_junction *x, t_signal **sp)
{
  int i, c = x->x_ctl.c_channels;
  float norm;
for (i=0;i<c;i++)
    {
      x->x_ctl.c_in[i] = sp[i]->s_vec;
      x->x_ctl.c_out[i] = sp[i+c]->s_vec;
    }
norm = c;
  norm = 1. / (norm);
  x->x_ctl.c_norm =  norm;
dsp_add(junction_perform, 2, &x->x_ctl, sp[0]->s_n);
/*    dsp_add(junction_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);*/
}
void junction_free(t_junction *x)
{
if (x->x_ctl.c_in) free (x->x_ctl.c_in);
  if (x->x_ctl.c_out) free (x->x_ctl.c_out);
  if (x->x_ctl.c_buffer) free (x->x_ctl.c_buffer);
  if (x->x_ctl.c_coef) free (x->x_ctl.c_coef);
}
t_class *junction_class;
void *junction_new(t_floatarg channels)
{
int l = ilog2(channels);
  int i,n;
t_junction *x = (t_junction *)pd_new(junction_class);
if (l<2) l = 2;
  if (l>4) l = 4;
n=1;
  while (l--) n *= 2;
for (i=1;i<n;i++) inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal")); 
  for (i=0;i<n;i++) outlet_new(&x->x_obj, gensym("signal"));
x->x_ctl.c_in = (float **)malloc(n*sizeof(float *));
  x->x_ctl.c_out = (float **)malloc(n*sizeof(float *));
  x->x_ctl.c_buffer = (float *)malloc(n*sizeof(float));
  x->x_ctl.c_coef = (float *)malloc(n*sizeof(float));
  x->x_ctl.c_channels = n;
junction_bang(x);
return (void *)x;
}
void junction_tilde_setup(void)
{
  //post("junction~ v0.1");
    junction_class = class_new(gensym("junction~"), (t_newmethod)junction_new,
    	(t_method)junction_free, sizeof(t_junction), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(junction_class, t_junction, x_f); 
    class_addmethod(junction_class, (t_method)junction_bang, gensym("bang"), 0);
    class_addmethod(junction_class, (t_method)junction_random, gensym("random"), A_FLOAT, 0);
    class_addmethod(junction_class, (t_method)junction_dsp, gensym("dsp"), 0);
}
--- NEW FILE: abs~.c ---
// since this is present in a lot of libs, it is conditionally compiled
#ifdef HAVE_ABS_TILDE
#include "m_pd.h"
#include <math.h>
/* ------------------------- abs~ -------------------------- */
static t_class *abs_class;
typedef struct _abs
{
    t_object x_obj;
} t_abs;
static t_int *abs_perform(t_int *w)
{
    t_abs *x = (t_abs *)(w[1]);
    t_float *in = (t_float *)(w[2]);
    t_float *out = (t_float *)(w[3]);
    int n = (int)(w[4]);
    while (n--)
    {
    	float f = *in++;
    	if (f < 0) f = -f;
    	*out++ = f;
    }
    return (w+5);
}
static void abs_dsp(t_abs *x, t_signal **sp)
{
    dsp_add(abs_perform, 4, x, sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
}
static void *abs_new(void)
{
    t_abs *x = (t_abs *)pd_new(abs_class);
    outlet_new(&x->x_obj, &s_signal);
    return (x);
}
void abs_tilde_setup(void)
{
    abs_class = class_new(gensym("abs~"), (t_newmethod)abs_new, 0,
    	sizeof(t_abs), 0, A_NULL);
    class_addmethod(abs_class, (t_method)nullfn, &s_signal, A_NULL);
    class_addmethod(abs_class, (t_method)abs_dsp, gensym("dsp"), A_NULL);
}
#endif
--- NEW FILE: qnorm~.c ---
/*
 *   qnorm.c  - quaternion normalization dsp object 
 *   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 <stdlib.h>
#include <stdio.h>
#include <math.h>
typedef struct qnormctl
{
  t_float *c_input[4];
  t_float *c_output[4];
} t_qnormctl;
typedef struct qnorm
{
  t_object x_obj;
  t_float x_f;
  t_qnormctl x_ctl;
} t_qnorm;
static t_int *qnorm_perform(t_int *word)
{
t_qnormctl *ctl     = (t_qnormctl *)(word[1]);
  t_int n             = (t_int)(word[2]);
  t_int i;
t_float *in0        = ctl->c_input[0];
  t_float *in1        = ctl->c_input[1];
  t_float *in2        = ctl->c_input[2];
  t_float *in3        = ctl->c_input[3];
t_float *out0       = ctl->c_output[0];
  t_float *out1       = ctl->c_output[1];
  t_float *out2       = ctl->c_output[2];
  t_float *out3       = ctl->c_output[3];
t_float w, x, y, z;
  t_float norm;
  t_float inorm;
for (i=0;i<n;i++)
    {
/* read input */
    w = *in0++;
    x = *in1++;
    y = *in2++;
    z = *in3++;
/* transform */
    norm = w * w;
    norm += x * x;
    norm += y * y;
    norm += z * z;
    
    inorm = (norm == 0.0f) ? (0.0f) : (1.0f / sqrt(norm));
/* write output */
    *out0++ = w * inorm;
    *out1++ = x * inorm;
    *out2++ = y * inorm;
    *out3++ = z * inorm;
    }
return (word+3);
}
static void qnorm_dsp(t_qnorm *x, t_signal **sp)
{
int i;
  for (i=0;i<4;i++)
    {
      x->x_ctl.c_input[i] = sp[i]->s_vec;
      x->x_ctl.c_output[i] = sp[i+4]->s_vec;
    }
dsp_add(qnorm_perform, 2, &x->x_ctl, sp[0]->s_n);
}
static void qnorm_free(t_qnorm *x)
{
}
t_class *qnorm_class;
static void *qnorm_new(t_floatarg channels)
{
    int i;
    t_qnorm *x = (t_qnorm *)pd_new(qnorm_class);
for (i=1;i<4;i++) inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal")); 
    for (i=0;i<4;i++) outlet_new(&x->x_obj, gensym("signal"));
return (void *)x;
}
void qnorm_tilde_setup(void)
{
  //post("qnorm~ v0.1");
    qnorm_class = class_new(gensym("qnorm~"), (t_newmethod)qnorm_new,
    	(t_method)qnorm_free, sizeof(t_qnorm), 0, 0);
    CLASS_MAINSIGNALIN(qnorm_class, t_qnorm, x_f); 
    class_addmethod(qnorm_class, (t_method)qnorm_dsp, gensym("dsp"), 0);
}
--- lattice.c DELETED ---
--- NEW FILE: sbosc~.c ---
/*
 *   sbosc.c  - smallband oscillator. periodic, linear interpolated frequency center.
 *   data organization is in (real, imag) pairs
 *   the first 2 components are (DC, NY)
 *   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>
#include <string.h>
#define LOGTABSIZE 10
#define TABSIZE (1<<LOGTABSIZE)
#define MASKTABSIZE (TABSIZE-1)
#define SHIFTTABSIZE ((sizeof(unsigned int) * 8) - LOGTABSIZE)
#define FRACTABSIZE (1<<SHIFTTABSIZE)
#define INVFRACTABSIZE (1.0f / (float)(FRACTABSIZE))
#define MASKFRACTABSIZE (FRACTABSIZE-1)
#define PITCHLIMIT 20.0f
static float costable[TABSIZE];
static inline void _exp_j2pi(unsigned int t, float *real, float *imag)
{
    unsigned int i1 = t >> SHIFTTABSIZE;
    float f2 = (t & MASKFRACTABSIZE) * INVFRACTABSIZE;
    unsigned int i2 = (i1+1) & MASKTABSIZE;
    unsigned int i3 = (i1 - (TABSIZE>>2)) & MASKTABSIZE;
    unsigned int i4 = (i2 + 1 - (TABSIZE>>2)) & MASKTABSIZE;
    float f1 = 1.0f - f2;
    float a1 = f1 * costable[i1];
    float a2 = f2 * costable[i2];
    float b1 = f1 * costable[i3];
    float b2 = f2 * costable[i4];
    *real = a1 + a2;
    *imag = b1 + b2;
}
static t_class *sbosc_tilde_class;
typedef struct _sbosc_tilde
{
    t_object x_obj;
    float x_f;
/* state vars */
    unsigned int x_phase;          // phase of main pitch osc
    unsigned int x_phase_inc;      // frequency of main pitch osc
    unsigned int x_harmonic;       // first harmonic
    float x_frac;                  // fraction of first harmonic
} t_sbosc_tilde;
static void *sbosc_tilde_new(void)
{
    t_sbosc_tilde *x = (t_sbosc_tilde *)pd_new(sbosc_tilde_class);
    x->x_phase = 0;
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));  
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("phase"));  
    outlet_new(&x->x_obj, gensym("signal"));
    outlet_new(&x->x_obj, gensym("signal"));
    x->x_f = 0;
    return (x);
}
static t_int *sbosc_tilde_perform(t_int *w)
{
    t_sbosc_tilde *x = (t_sbosc_tilde *)(w[1]);
    t_float *pitch = (t_float *)(w[2]);
    t_float *center= (t_float *)(w[3]);
    t_float *out_real = (t_float *)(w[4]);
    t_float *out_imag = (t_float *)(w[5]);
    int n = (int)(w[6]);
    int i;
t_float pitch_to_phase =  4294967295.0f / sys_getsr();
for (i = 0; i < n; i++)
    {
    float p = *pitch++;
    float c = *center++;
      	float r1,r2,i1,i2;
/* compute harmonic mixture */
    unsigned int h1 = x->x_phase * x->x_harmonic;
    unsigned int h2 = h1 + x->x_phase;
    _exp_j2pi(h1, &r1, &i1);
    _exp_j2pi(h2, &r2, &i2);
    r1 *= x->x_frac;
    i1 *= x->x_frac;
    r2 *= 1.0f - x->x_frac;
    i2 *= 1.0f - x->x_frac;
*out_real++ = r1 + r2;
    *out_imag++ = i1 + i2;
x->x_phase += x->x_phase_inc;
/* check for phase wrap & update osc */
    if ((x->x_phase <= x->x_phase_inc)) 
        {
    	float p_plus = (p < 0.0f) ? -p : p;
    	float p_limit = (p_plus < PITCHLIMIT) ? PITCHLIMIT : p_plus;
    	float c_plus = (c < 0.0f) ? -c : c;
    	float harmonic = c_plus/p_limit;
    	x->x_phase_inc = pitch_to_phase * p_limit;
    	x->x_harmonic = harmonic;
    	x->x_frac = 1.0f - (harmonic - x->x_harmonic);
        }
}
    
    return (w+7);
}
static void sbosc_tilde_dsp(t_sbosc_tilde *x, t_signal **sp)
{
    dsp_add(sbosc_tilde_perform, 6, x,
        sp[0]->s_vec, sp[1]->s_vec, sp[2]->s_vec, sp[3]->s_vec, sp[0]->s_n);
}
static void sbosc_tilde_free(t_sbosc_tilde *x)
{
}
static void sbosc_tilde_phase(t_sbosc_tilde *x, t_floatarg f)
{
    x->x_phase = f * (1.0f / 4294967295.0f);
}
void sbosc_tilde_setup(void)
{
    int i;
// init tables
    for (i=0; i<TABSIZE; i++)
    costable[i] = cos(2.0 * M_PI * (double)i / (double)TABSIZE);
// class setup
    sbosc_tilde_class = class_new(gensym("sbosc~"),
        (t_newmethod)sbosc_tilde_new, (t_method)sbosc_tilde_free,
        sizeof(t_sbosc_tilde), 0, A_DEFSYM, 0);
    CLASS_MAINSIGNALIN(sbosc_tilde_class, t_sbosc_tilde, x_f);
    class_addmethod(sbosc_tilde_class, (t_method)sbosc_tilde_dsp,
        gensym("dsp"), 0);
    class_addmethod(sbosc_tilde_class, (t_method)sbosc_tilde_phase,
        gensym("phase"), A_FLOAT, 0);
}
--- NEW FILE: cmath~.c ---
/*
 *   cmath.c  - some complex math dsp objects
 *   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>
#define MINNORM 0.0000000001
typedef struct cmath
{
    t_object x_obj;
    t_float x_f;
    t_perfroutine x_perf;
} t_cmath;
static t_int *cmath_perform_clog(t_int *w)
{
    t_float *inx    = (float *)(w[2]);
    t_float *iny    = (float *)(w[3]);
    t_float *outx    = (float *)(w[4]);
    t_float *outy    = (float *)(w[5]);
    t_int i;
    t_int n = (t_int)(w[1]);
    t_float x;
while (n--){
    float x = *inx++;
    float y = *iny++;
    float norm = sqrt(x*x + y*y);
    float arg = atan2(y, x);
    if (norm < MINNORM){
        norm = MINNORM;
    }
    *outx++ = log(norm);
    *outy++ = arg;
    }
return (w+6);
}
static t_int *cmath_perform_cexp(t_int *w)
{
    t_float *inx    = (float *)(w[2]);
    t_float *iny    = (float *)(w[3]);
    t_float *outx    = (float *)(w[4]);
    t_float *outy    = (float *)(w[5]);
    t_int i;
    t_int n = (t_int)(w[1]);
    t_float x;
while (n--){
    float x = *inx++;
    float y = *iny++;
    float norm = exp(x);
    *outx++ = norm * cos(y);
    *outy++ = norm * sin(y);
    }
return (w+6);
}
static t_int *cmath_perform_nfft(t_int *w)
{
    t_float *inx    = (float *)(w[2]);
    t_float *iny    = (float *)(w[3]);
    t_float *outx    = (float *)(w[4]);
    t_float *outy    = (float *)(w[5]);
    t_int i;
    t_int n = (t_int)(w[1]);
    t_float x;
    t_float scale = 1.0f / (sqrt((float)n));
mayer_fft(n, inx, outx);
while (n--){
    float x = *inx++;
    float y = *iny++;
    *outx++ = scale * x;
    *outy++ = scale * y;
    }
return (w+6);
}
static t_int *cmath_perform_nifft(t_int *w)
{
    t_float *inx    = (float *)(w[2]);
    t_float *iny    = (float *)(w[3]);
    t_float *outx    = (float *)(w[4]);
    t_float *outy    = (float *)(w[5]);
    t_int i;
    t_int n = (t_int)(w[1]);
    t_float x;
    t_float scale = 1.0f / (sqrt((float)n));
mayer_ifft(n, inx, outx);
while (n--){
    float x = *inx++;
    float y = *iny++;
    *outx++ = scale * x;
    *outy++ = scale * y;
    }
return (w+6);
}
static void cmath_dsp(t_cmath *x, t_signal **sp)
{
    dsp_add(x->x_perf, 5, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec, sp[2]->s_vec, sp[3]->s_vec);
}                                  
void cmath_free(void)
{
}
t_class *cmath_class;
t_cmath *cmath_new_common(void)
{
    t_cmath *x = (t_cmath *)pd_new(cmath_class);
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));  
    outlet_new(&x->x_obj, gensym("signal")); 
    outlet_new(&x->x_obj, gensym("signal")); 
    return x;
}
#define DEFNEWCMATH(name, perfmethod)		\
void * name (void)				\
{						\
    t_cmath *x = cmath_new_common();		\
    x->x_perf = perfmethod ;			\
    return (void*)x;				\
}
DEFNEWCMATH(cmath_new_clog, cmath_perform_clog)
DEFNEWCMATH(cmath_new_cexp, cmath_perform_cexp)
DEFNEWCMATH(cmath_new_nfft, cmath_perform_nfft)
DEFNEWCMATH(cmath_new_nifft, cmath_perform_nifft)
void cmath_tilde_setup(void)
{
  //post("cmath~ v0.1");
    cmath_class = class_new(gensym("clog~"), (t_newmethod)cmath_new_clog,
    	(t_method)cmath_free, sizeof(t_cmath), 0, 0);
class_addcreator((t_newmethod)cmath_new_cexp, gensym("cexp~"), A_NULL);
    class_addcreator((t_newmethod)cmath_new_nfft, gensym("nfft~"), A_NULL);
    class_addcreator((t_newmethod)cmath_new_nifft, gensym("nifft~"), A_NULL);
CLASS_MAINSIGNALIN(cmath_class, t_cmath, x_f);
class_addmethod(cmath_class, (t_method)cmath_dsp, gensym("dsp"), 0); 
}
--- bfft.c DELETED ---
--- NEW FILE: ead~.c ---
/*
 *   ead.c  -  exponential attack decay envelope 
 *   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 "extlib_util.h"
/* pointer to */
t_class *ead_class;
/* state data fpr attack/decay dsp plugin */
typedef struct eadctl
{
  t_float c_attack;
  t_float c_decay;
  t_float c_state;
  t_int c_target;
} t_eadctl;
/* object data structure */
typedef struct ead
{
  t_object x_obj;  
  t_eadctl x_ctl;   
} t_ead;
static void ead_attack(t_ead *x, t_floatarg f)
{
  x->x_ctl.c_attack = milliseconds_2_one_minus_realpole(f);
}
static void ead_decay(t_ead *x, t_floatarg f)
{
  x->x_ctl.c_decay = milliseconds_2_one_minus_realpole(f);
}
static void ead_start(t_ead *x)
{
    /* reset state if necessary to prevent skipping */
// always reset, seems to be safest
    //if (x->x_ctl.c_target == 1)
x->x_ctl.c_state = 0.0f; 
    x->x_ctl.c_target = 1;
}
/* dsp callback function, not a method */
static t_int *ead_perform(t_int *w)
{
/* interprete arguments */
    t_float *out    = (float *)(w[3]);
    t_eadctl *ctl  = (t_eadctl *)(w[1]);
    t_float attack  = ctl->c_attack;
    t_float decay   = ctl->c_decay;
    t_float state   = ctl->c_state;
    t_int n = (t_int)(w[2]);
t_int i;
/* A/D code */
for (i = 0; i < n; i++){
    switch(ctl->c_target){
    case 1:
        /* attack phase */
        *out++ = state;
        state += attack*(1 - state);
        ctl->c_target = (state <= ENVELOPE_MAX);
        break;
    default:
        /* decay phase */
        *out++ = state;
        state -= decay*state;
        break;
    }
    
    }
/* save state */
    ctl->c_state = IS_DENORMAL(state) ? 0 : state;
return (w+4); 
}
static void ead_dsp(t_ead *x, t_signal **sp)
{
  dsp_add(ead_perform, 3, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec);
}
/* destructor */
static void ead_free(void)
{
}
/* constructor */
static void *ead_new(t_floatarg attack, t_floatarg decay, t_floatarg sustain, t_floatarg release)
{
    /* create instance */
    t_ead *x = (t_ead *)pd_new(ead_class);
    /* create new inlets, convert incoming message float to attack/decay */
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("attack"));  
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("decay"));  
    /* create a dsp outlet */
    outlet_new(&x->x_obj, gensym("signal"));
/* initialize */
    x->x_ctl.c_state  = 0;
    x->x_ctl.c_target = 0;
    ead_attack(x, attack);
    ead_decay(x, decay);
/* return instance */
    return (void *)x;
}
void ead_tilde_setup(void)
{
  //post("ead~ v0.1");
ead_class = class_new(gensym("ead~"), (t_newmethod)ead_new,
    	(t_method)ead_free, sizeof(t_ead), 0,  A_DEFFLOAT, A_DEFFLOAT, 0);
class_addmethod(ead_class, (t_method)ead_start, gensym("start"), 0);
    class_addmethod(ead_class, (t_method)ead_start, gensym("bang"), 0);
    class_addmethod(ead_class, (t_method)ead_dsp, gensym("dsp"), 0); 
    class_addmethod(ead_class, (t_method)ead_attack, gensym("attack"), A_FLOAT, 0);
    class_addmethod(ead_class, (t_method)ead_decay, gensym("decay"), A_FLOAT, 0);
}
--- sbosc.c DELETED ---
--- permut.c DELETED ---
--- NEW FILE: bfft~.c ---
/*
 *   bfft.c  - code for some fourrier transform variants and utility functions
 *   data organization is in (real, imag) pairs
 *   the first 2 components are (DC, NY)
 *   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 bfftctl
{
    t_int c_levels;
    char c_name[16];
    t_int *c_clutter;
    t_int *c_unclutter;
    t_int c_kill_DC;
    t_int c_kill_NY;
} t_bfftctl;
typedef struct bfft
{
    t_object x_obj;
    t_float x_f;
    t_bfftctl x_ctl;
} t_bfft;
t_class *bfft_class, *ibfft_class, *fht_class;
static inline void bfft_perform_permutation(t_float *S, int n, t_int *f)
{
  t_int k,l;
  t_float swap;
  for(k=0; k<n; k++)
    {
      l = f[k];
      while (l<k) l = f[l];
      swap = S[k];
      S[k] = S[l];
      S[l] = swap;
    }
}
static void bfft_permutation(t_bfft *x, t_int n){
t_bfftctl *ctl = &x->x_ctl;
  int i;
if (ctl->c_clutter)   free(ctl->c_clutter);
  if (ctl->c_unclutter) free(ctl->c_unclutter);
ctl->c_clutter = (t_int *)malloc(n*sizeof(t_int));
  ctl->c_unclutter = (t_int *)malloc(n*sizeof(t_int));
ctl->c_unclutter[0] = 0;
  ctl->c_unclutter[1] = n/2;
  for (i=1; i<n/2; i++){
    ctl->c_unclutter[2*i] = i;
    ctl->c_unclutter[2*i+1] = n-i;
  }
for(i=0; i<n; i++)
    ctl->c_clutter[ctl->c_unclutter[i]] = i;
return;
/* debug */
  /*  for(k=0; k<n; k++)
  **    printf("clutter[%d] = %d\n", k, ctl->c_clutter[k]);
  **  for(k=0; k<n; k++)
  **    printf("unclutter[%d] = %d\n", k, ctl->c_unclutter[k]);
  **
  ** exit(1);
  */
}
static t_int *bfft_perform(t_int *w)
{
t_float *in     = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_bfftctl *ctl  = (t_bfftctl *)(w[1]);
  t_int n = (t_int)(w[2]);
  t_float scale = sqrt(1.0f / (float)(n));
mayer_fht(out, n);
  bfft_perform_permutation(out, n, ctl->c_unclutter);
while (n--) *out++ *= scale;
return (w+5);
}
static t_int *ibfft_perform(t_int *w)
{
t_float *in     = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_bfftctl *ctl  = (t_bfftctl *)(w[1]);
  t_int n = (t_int)(w[2]);
  t_float scale = sqrt(1.0f / (float)(n));
if (ctl->c_kill_DC) {out[0] = 0.0f;}
  if (ctl->c_kill_NY) {out[1] = 0.0f;}
bfft_perform_permutation(out, n, ctl->c_clutter);
  mayer_fht(out, n);
while (n--) *out++ *= scale;
return (w+5);
}
static t_int *fht_perform(t_int *w)
{
t_float *in     = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_bfftctl *ctl  = (t_bfftctl *)(w[1]);
t_int n = (t_int)(w[2]);
mayer_fht(out, n);
return (w+5);
}
static void bfft_dsp(t_bfft *x, t_signal **sp)
{
int n = sp[0]->s_n;
  t_float *in = sp[0]->s_vec;
  t_float *out = sp[1]->s_vec;
bfft_permutation(x, n);
if (in != out)
    {
      dsp_add_copy(in,out,n);
      in = out;
    }
dsp_add(bfft_perform, 4, &x->x_ctl, n, in, out);
}
static void ibfft_dsp(t_bfft *x, t_signal **sp)
{
int n = sp[0]->s_n;
  t_float *in = sp[0]->s_vec;
  t_float *out = sp[1]->s_vec;
bfft_permutation(x, n);
if (in != out)
    {
      dsp_add_copy(in,out,n);
      in = out;
    }
dsp_add(ibfft_perform, 4, &x->x_ctl, n, in, out);
}
static void fht_dsp(t_bfft *x, t_signal **sp)
{
int n = sp[0]->s_n;
  t_float *in = sp[0]->s_vec;
  t_float *out = sp[1]->s_vec;
if (in != out)
    {
      dsp_add_copy(in,out,n);
      in = out;
    }
dsp_add(fht_perform, 4, &x->x_ctl, n, in, out);
}
static void bfft_free(t_bfft *x)
{
if (x->x_ctl.c_clutter)   free(x->x_ctl.c_clutter);
  if (x->x_ctl.c_unclutter) free(x->x_ctl.c_unclutter);
}
static void *bfft_new(void)
{
    t_bfft *x = (t_bfft *)pd_new(bfft_class);
    int i;
outlet_new(&x->x_obj, gensym("signal"));
sprintf(x->x_ctl.c_name,"bfft");
x->x_ctl.c_clutter = NULL;
    x->x_ctl.c_unclutter = NULL;
return (void *)x;
}
static void *ibfft_new(t_symbol *s)
{
    t_bfft *x = (t_bfft *)pd_new(ibfft_class);
    int i;
outlet_new(&x->x_obj, gensym("signal"));
if (s == gensym("killDCNY")){
    x->x_ctl.c_kill_DC = 1;
    x->x_ctl.c_kill_NY = 1;
    post("ibfft: removing DC and NY components.");
    }
    else{
    x->x_ctl.c_kill_DC = 0;
    x->x_ctl.c_kill_NY = 0;
    }
x->x_ctl.c_clutter = NULL;
    x->x_ctl.c_unclutter = NULL;
sprintf(x->x_ctl.c_name,"ibfft");
return (void *)x;
}
static void *fht_new(void)
{
    t_bfft *x = (t_bfft *)pd_new(fht_class);
    int i;
outlet_new(&x->x_obj, gensym("signal"));
x->x_ctl.c_clutter = NULL;
    x->x_ctl.c_unclutter = NULL;
sprintf(x->x_ctl.c_name,"fht");
return (void *)x;
}
void bfft_tilde_setup(void)
{
  //post("bfft~ v0.1");
    bfft_class = class_new(gensym("bfft~"), (t_newmethod)bfft_new,
    	(t_method)bfft_free, sizeof(t_bfft), 0, 0);
    CLASS_MAINSIGNALIN(bfft_class, t_bfft, x_f);
    class_addmethod(bfft_class, (t_method)bfft_dsp, gensym("dsp"), 0);
ibfft_class = class_new(gensym("ibfft~"), (t_newmethod)ibfft_new,
    	(t_method)bfft_free, sizeof(t_bfft), 0, A_DEFSYMBOL, A_NULL);
/* add the more logical bifft~ alias */
    class_addcreator((t_newmethod)ibfft_new, 
    	     gensym("bifft~"), 0, A_DEFSYMBOL, A_NULL);
CLASS_MAINSIGNALIN(ibfft_class, t_bfft, x_f);
    class_addmethod(ibfft_class, (t_method)ibfft_dsp, gensym("dsp"), 0);
fht_class = class_new(gensym("fht~"), (t_newmethod)fht_new,
    	(t_method)bfft_free, sizeof(t_bfft), 0, 0);
CLASS_MAINSIGNALIN(fht_class, t_bfft, x_f);
    class_addmethod(fht_class, (t_method)fht_dsp, gensym("dsp"), 0);
}
--- dynwav.c DELETED ---
--- NEW FILE: eblosc~.c ---
/*
 *   eblosc.c  - bandlimited oscillators with infinite support discontinuities 
 *   using minimum phase impulse, step & ramp
 *   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>
#include <string.h>
#include "filters.h"
typedef unsigned long long u64;
typedef unsigned long u32;
#define LPHASOR      (8*sizeof(u32)) // the phasor logsize
#define VOICES       8 // the number of oscillators
#define CUTOFF       0.8f // fraction of nyquist for impulse cutoff
typedef struct ebloscctl
{
    t_float c_pole[VOICES*2];      // complex poles
    t_float c_gain[VOICES*2];      // complex gains (waveform specific constants)
    t_float c_state[VOICES*2];     // complex state
u32 c_phase;                // phase of main oscillator
    u32 c_phase2;               // phase of secondairy oscillator
    t_float c_prev_amp;         // previous input of comparator
    t_float c_phase_inc_scale;
    t_float c_scale;
    t_float c_scale_update;
    t_symbol *c_waveform;
} t_ebloscctl;
typedef struct eblosc
{
  t_object x_obj;
  t_float x_f;
  t_ebloscctl x_ctl;
} t_eblosc;
/* phase converters */
static inline float _phase_to_float(u32 p){return ((float)p) * (1.0f / 4294967296.0f);}
static inline u32 _float_to_phase(float f){return (u32)(f * 4294967296.0f);}
/* get one sample from the oscillator bank and perform time tick */
static inline t_float _osc_tick(t_ebloscctl *ctl)
{
    float sum = 0.0f;
    int i;
    /* sum  all voices */
    for (i=0; i<VOICES*2; i+=2){
    /* rotate state */
    vcmul2(ctl->c_state+i, ctl->c_pole+i);
/* get real part and add to output */
    sum += ctl->c_state[0];
    }
return sum;
}
/* add shifted impulse */
static inline void _add_impulse(t_ebloscctl *ctl, t_float t0)
{
    int i;
    for (i=0; i<VOICES*2; i+=2){
    /* contribution is a_i z_i^t_0 */
float real = 1.0f;
    float imag = 0.0f;
        
    ctl->c_state[0] += real;
    ctl->c_state[1] += imag;
    }
}
/* add step */
static inline void _add_step(t_ebloscctl *ctl)
{
    int i;
    for (i=0; i<VOICES*2; i+=2){
    /* contribution is a_i (1 - z_i) */
float real = 1.0f;
    float imag = 0.0f;
        
    ctl->c_state[0] += real;
    ctl->c_state[1] += imag;
    }
}
/* add shifted step */
static inline void _add_shifted_step(t_ebloscctl *ctl, t_float t0)
{
    int i;
    for (i=0; i<VOICES*2; i+=2){
    /* contribution is a_i (1 - z_i^t_0) */
float real = 1.0f;
    float imag = 0.0f;
        
    ctl->c_state[0] += real;
    ctl->c_state[1] += imag;
    }
}
#if 0
/* update waveplayers on zero cross */
static void _bang_comparator(t_ebloscctl *ctl, float prev, float curr)
{
/* check for sign change */
    if ((prev * curr) < 0.0f){
int voice;
/* determine the location of the discontinuity (in oversampled coordiates
      using linear interpolation */
float f = (float)S * curr / (curr - prev);
/* get the offset in the oversample table */
u32 table_index = (u32)f;
/* determine the fractional part (in oversampled coordinates)
       for linear interpolation */
float table_frac_index = f - (float)table_index;
/* set state (+ or -) */
ctl->c_state =  (curr > 0.0f) ? 0.5f : -0.5f;
    
    /* steal the oldest voice */
voice = ctl->c_next_voice++;
    ctl->c_next_voice &= VOICES-1;
        
    /* initialize the new voice index and interpolation fraction */
ctl->c_index[voice] = table_index;
    ctl->c_frac[voice] = table_frac_index;
    ctl->c_vscale[voice] = -ctl->c_scale * 2.0f * ctl->c_state;
}
}
/* advance phasor and update waveplayers on phase wrap */
static void _bang_phasor(t_ebloscctl *ctl, float freq)
{
    u32 phase = ctl->c_phase;
    u32 phase_inc; 
    u32 oldphase;
    int voice;
    float scale = ctl->c_scale;
/* get increment */
    float inc = freq * ctl->c_phase_inc_scale;
/* calculate new phase
       the increment (and the phase) should be a multiple of S */
    if (inc < 0.0f) inc = -inc;
    phase_inc = ((u32)inc) & ~(S-1);
    oldphase = phase;
    phase += phase_inc;
/* check for phase wrap */
    if (phase < oldphase){
    u32 phase_inc_decimated = phase_inc >> LOVERSAMPLE;
    u32 table_index;
    u32 table_phase;
    
    /* steal the oldest voice if we have a phase wrap */
        
    voice = ctl->c_next_voice++;
    ctl->c_next_voice &= VOICES-1;
        
    /* determine the location of the discontinuity (in oversampled coordinates)
       which is S * (new phase) / (increment) */
        
    table_index = phase / phase_inc_decimated;
        
    /* determine the fractional part (in oversampled coordinates)
       for linear interpolation */
table_phase = phase - (table_index * phase_inc_decimated);
        
    /* use it to initialize the new voice index and interpolation fraction */
        
    ctl->c_index[voice] = table_index;
    ctl->c_frac[voice] = (float)table_phase / (float)phase_inc_decimated;
    ctl->c_vscale[voice] = scale;
    scale = scale * ctl->c_scale_update;
}
/* save state */
    ctl->c_phase = phase;
    ctl->c_scale = scale;
}
/* the 2 oscillator version:
   the second osc can reset the first osc's phase (hence it determines the pitch)
   the first osc determines the waveform */
static void _bang_hardsync_phasor(t_ebloscctl *ctl, float freq, float freq2)
{
    u32 phase = ctl->c_phase;
    u32 phase2 = ctl->c_phase2;
    u32 phase_inc; 
    u32 phase_inc2; 
    u32 oldphase;
    u32 oldphase2;
    int voice;
    float scale = ctl->c_scale;
/* get increment */
    float inc = freq * ctl->c_phase_inc_scale;
    float inc2 = freq2 * ctl->c_phase_inc_scale;
/* calculate new phases
       the increment (and the phase) should be a multiple of S */
/* save previous phases */
    oldphase = phase;
    oldphase2 = phase2;
/* update second osc */
    if (inc2 < 0.0f) inc2 = -inc2;
    phase_inc2 = ((u32)inc2) & ~(S-1);
    phase2 += phase_inc2;
/* update first osc (freq should be >= freq of sync osc */
    if (inc < 0.0f) inc = -inc;
    phase_inc = ((u32)inc) & ~(S-1);
    if (phase_inc < phase_inc2) phase_inc = phase_inc2;
    phase += phase_inc;
/* check for sync discontinuity (osc 2) */
    if (phase2 < oldphase2) {
/* adjust phase depending on the location of the discontinuity in phase2:
       phase/phase_inc == phase2/phase_inc2 */
    
    u64 pi = phase_inc >> LOVERSAMPLE;
    u64 pi2 = phase_inc2 >> LOVERSAMPLE;
    u64 lphase = ((u64)phase2 * pi) / pi2;
    phase = lphase & ~(S-1);
    }
/* check for phase discontinuity (osc 1) */
    if (phase < oldphase){
    u32 phase_inc_decimated = phase_inc >> LOVERSAMPLE;
    u32 table_index;
    u32 table_phase;
    float stepsize;
    
    /* steal the oldest voice if we have a phase wrap */
        
    voice = ctl->c_next_voice++;
    ctl->c_next_voice &= VOICES-1;
        
    /* determine the location of the discontinuity (in oversampled coordinates)
       which is S * (new phase) / (increment) */
table_index = phase / phase_inc_decimated;
        
    /* determine the fractional part (in oversampled coordinates)
       for linear interpolation */
table_phase = phase - (table_index * phase_inc_decimated);
/* determine the step size
       as opposed to saw/impulse waveforms, the step is not always equal to one. it is:
           oldphase - phase + phase_inc 
       but for the unit step this will overflow to zero, so we
       reduce the bit depth to prevent overflow */
stepsize = _phase_to_float(((oldphase-phase) >> LOVERSAMPLE)
    			   + phase_inc_decimated) * (float)S;
        
    /* use it to initialize the new voice index and interpolation fraction */
        
    ctl->c_index[voice] = table_index;
    ctl->c_frac[voice] = (float)table_phase / (float)phase_inc_decimated;
    ctl->c_vscale[voice] = scale * stepsize;
    scale = scale * ctl->c_scale_update;
}
/* save state */
    ctl->c_phase = phase;
    ctl->c_phase2 = phase2;
    ctl->c_scale = scale;
}
static t_int *eblosc_perform_hardsync_saw(t_int *w)
{
    t_float *freq     = (float *)(w[3]);
    t_float *freq2     = (float *)(w[4]);
    t_float *out      = (float *)(w[5]);
    t_ebloscctl *ctl  = (t_ebloscctl *)(w[1]);
    t_int n           = (t_int)(w[2]);
    t_int i;
/* set postfilter cutoff */
    ctl->c_butter->setButterHP(0.85f * (*freq / sys_getsr()));
while (n--) {
    float frequency = *freq++;
    float frequency2 = *freq2++;
/* get the bandlimited discontinuity */
    float sample = _get_bandlimited_discontinuity(ctl, bls);
/* add aliased sawtooth wave */
    sample += _phase_to_float(ctl->c_phase) - 0.5f;
/* highpass filter output to remove DC offset and low frequency aliasing */
    ctl->c_butter->BangSmooth(sample, sample, 0.05f);
/* send to output */
    *out++ = sample;
/* advance phasor */
    _bang_hardsync_phasor(ctl, frequency2, frequency);
    
    }
return (w+6);
}
static t_int *eblosc_perform_saw(t_int *w)
{
    t_float *freq     = (float *)(w[3]);
    t_float *out      = (float *)(w[4]);
    t_ebloscctl *ctl  = (t_ebloscctl *)(w[1]);
    t_int n           = (t_int)(w[2]);
    t_int i;
while (n--) {
    float frequency = *freq++;
/* get the bandlimited discontinuity */
    float sample = _get_bandlimited_discontinuity(ctl, bls);
/* add aliased sawtooth wave */
    sample += _phase_to_float(ctl->c_phase) - 0.5f;
/* send to output */
    *out++ = sample;
/* advance phasor */
    _bang_phasor(ctl, frequency);
    
    }
return (w+5);
}
static t_int *eblosc_perform_pulse(t_int *w)
{
    t_float *freq     = (float *)(w[3]);
    t_float *out      = (float *)(w[4]);
    t_ebloscctl *ctl  = (t_ebloscctl *)(w[1]);
    t_int n           = (t_int)(w[2]);
    t_int i;
/* set postfilter cutoff */
    ctl->c_butter->setButterHP(0.85f * (*freq / sys_getsr()));
while (n--) {
    float frequency = *freq++;
/* get the bandlimited discontinuity */
    float sample = _get_bandlimited_discontinuity(ctl, bli);
/* highpass filter output to remove DC offset and low frequency aliasing */
    ctl->c_butter->BangSmooth(sample, sample, 0.05f);
/* send to output */
    *out++ = sample;
/* advance phasor */
    _bang_phasor(ctl, frequency);
    
    }
return (w+5);
}
static t_int *eblosc_perform_comparator(t_int *w)
{
    t_float *amp      = (float *)(w[3]);
    t_float *out      = (float *)(w[4]);
    t_ebloscctl *ctl  = (t_ebloscctl *)(w[1]);
    t_int n           = (t_int)(w[2]);
    t_int i;
    t_float prev_amp = ctl->c_prev_amp;
while (n--) {
    float curr_amp = *amp++;
/* exact zero won't work for zero detection (sic) */
    if (curr_amp == 0.0f) curr_amp = 0.0000001f;
/* get the bandlimited discontinuity */
    float sample = _get_bandlimited_discontinuity(ctl, bls);
/* add the block wave state */
    sample += ctl->c_state;
/* send to output */
    *out++ = sample;
/* advance phasor */
    _bang_comparator(ctl, prev_amp, curr_amp);
prev_amp = curr_amp;
    
    }
ctl->c_prev_amp = prev_amp;
return (w+5);
}
static void eblosc_phase(t_eblosc *x, t_float f)
{
    x->x_ctl.c_phase = _float_to_phase(f);
    x->x_ctl.c_phase2 = _float_to_phase(f);
}
static void eblosc_phase1(t_eblosc *x, t_float f)
{
    x->x_ctl.c_phase = _float_to_phase(f);
}
static void eblosc_phase2(t_eblosc *x, t_float f)
{
    x->x_ctl.c_phase2 = _float_to_phase(f);
}
static void eblosc_dsp(t_eblosc *x, t_signal **sp)
{
  int n = sp[0]->s_n;
/* set sampling rate scaling for phasors */
  x->x_ctl.c_phase_inc_scale = 4.0f * (float)(1<<(LPHASOR-2)) / sys_getsr();
/* setup & register the correct process routine depending on the waveform */
/* 2 osc */
  if (x->x_ctl.c_waveform == gensym("syncsaw")){
      x->x_ctl.c_scale = 1.0f;
      x->x_ctl.c_scale_update = 1.0f;
      dsp_add(eblosc_perform_hardsync_saw, 5, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec, sp[2]->s_vec);
  }
/* 1 osc */
  else if (x->x_ctl.c_waveform == gensym("pulse")){
      x->x_ctl.c_scale = 1.0f;
      x->x_ctl.c_scale_update = 1.0f;
      dsp_add(eblosc_perform_pulse, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
  }
  else if (x->x_ctl.c_waveform == gensym("pulse2")){
      x->x_ctl.c_phase_inc_scale *= 2;
      x->x_ctl.c_scale = 1.0f;
      x->x_ctl.c_scale_update = -1.0f;
      dsp_add(eblosc_perform_pulse, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
  }
  else if (x->x_ctl.c_waveform == gensym("comparator")){
      x->x_ctl.c_scale = 1.0f;
      x->x_ctl.c_scale_update = 1.0f;
      dsp_add(eblosc_perform_comparator, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
  }
  else{
       x->x_ctl.c_scale = 1.0f;
      x->x_ctl.c_scale_update = 1.0f;
      dsp_add(eblosc_perform_saw, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
  }
}                                  
static void eblosc_free(t_eblosc *x)
{
    delete x->x_ctl.c_butter;
}
t_class *eblosc_class;
static void *eblosc_new(t_symbol *s)
{
    t_eblosc *x = (t_eblosc *)pd_new(eblosc_class);
    int i;
/* out 1 */
    outlet_new(&x->x_obj, gensym("signal"));
/* optional signal inlets */
    if (s == gensym("syncsaw")){
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));  
    }
/* optional phase inlet */
    if (s != gensym("comparator")){
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("phase"));  
    }
/* create the postfilter */
    x->x_ctl.c_butter = new DSPIfilterSeries(3);
/* init oscillators */
    for (i=0; i<VOICES; i++) {
      x->x_ctl.c_index[i] = N-2;
      x->x_ctl.c_frac[i] = 0.0f;
    }
/* init rest of state data */
    eblosc_phase(x, 0);
    eblosc_phase2(x, 0);
    x->x_ctl.c_state = 0.0;
    x->x_ctl.c_prev_amp = 0.0;
    x->x_ctl.c_next_voice = 0;
    x->x_ctl.c_scale = 1.0f;
    x->x_ctl.c_scale_update = 1.0f;
    x->x_ctl.c_waveform = s;
return (void *)x;
}
extern "C"
{
    void eblosc_tilde_setup(void)
    {
    //post("eblosc~ v0.1");
    
    build_tables();
    
    eblosc_class = class_new(gensym("eblosc~"), (t_newmethod)eblosc_new,
    			(t_method)eblosc_free, sizeof(t_eblosc), 0, A_DEFSYMBOL, A_NULL);
    CLASS_MAINSIGNALIN(eblosc_class, t_eblosc, x_f);
    class_addmethod(eblosc_class, (t_method)eblosc_dsp, gensym("dsp"), A_NULL); 
    class_addmethod(eblosc_class, (t_method)eblosc_phase, gensym("phase"), A_FLOAT, A_NULL); 
    class_addmethod(eblosc_class, (t_method)eblosc_phase2, gensym("phase2"), A_FLOAT, A_NULL);
}
}
#endif
--- NEW FILE: scrollgrid1D~.c ---
/*
 *   scrollgrid1D.c  -  1D scroll grid attractor
 *   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.
 */
/*  1D scroll grid attractor
    for more information see:
Yalcin M., Ozoguz S., Suykens J.A.K., Vandewalle J., 
    ``Families of Scroll Grid Attractors'', 
    International Journal of Bifurcation and Chaos, vol. 12, no. 1, Jan. 2002, pp. 23-41.
this file implements a digital variant of the method introduced in the paper,
    so that it can be used as a parametrizable, bounded chatotic oscillator.
    in short it is a switched linear system, with some added hard limiting to 
    convert unstable oscillations into stable ones.
*/
#include "m_pd.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>  
#include "filters.h"
typedef struct scrollgrid1Dctl
{
    t_float c_x, c_y, c_z; /* state */
} t_scrollgrid1Dctl;
typedef struct scrollgrid1D
{
    t_object x_obj;
    t_float x_f;
    t_scrollgrid1Dctl x_ctl;
} t_scrollgrid1D;
static inline float _fixedpoint(float x, int n)
{
  int ix = (x + 0.5f);
  if (ix < 0) ix = 0;
  else if (ix >= n) ix = n-1;
  return (float)ix;
}
static inline float _sat(float x, float upper)
{
    float lower = -1.0f;
    if (x < lower) x = lower;
    else if (x > upper) x = upper;
    return x;
}
static t_int *scrollgrid1D_perform(t_int *w)
{
t_float *freq    = (float *)(w[3]);
  t_float *t1      = (float *)(w[4]);
  t_float *t2      = (float *)(w[5]);
  t_float *order   = (float *)(w[6]);
  t_float *outx    = (float *)(w[7]);
  t_float *outy    = (float *)(w[8]);
  t_float *outz    = (float *)(w[9]);
  t_scrollgrid1Dctl *ctl    = (t_scrollgrid1Dctl *)(w[1]);
  t_int n          = (t_int)(w[2]);
t_int i;
  t_float inv_sr = 1.0f /sys_getsr();
  t_float state[3] = {ctl->c_x, ctl->c_y, ctl->c_z};
  t_float c,f;
  t_float pole[2], r1, r2;
  t_int o;
  t_float x,y,z;
for (i=0; i<n; i++){
/* get params */
      r1 = exp(1000.0f * inv_sr / (0.01f + fabs(*t1++)));
      r2 = exp(-1000.0f * inv_sr / (0.01f + fabs(*t2++)));
      f = *freq++;
      o = (int)(*order++);
      if (o < 2) o = 2;
      pole[0] = r1 * cos(2.0f * M_PI * inv_sr * f);
      pole[1] = r1 * sin(2.0f * M_PI * inv_sr * f);
/* debug */
      //post("%f", r1);
/* base transform + clipping to prevent blowup */
      x = _sat(0.5f * (state[0] - state[2]), (float)o); /* projection onto axis containing fixed */
      y = _sat(0.5f * state[1], 1.0f);                         /* the "pure" oscillation axis */
      z = _sat(0.5f * (state[0] + state[2]), 1.0f);     /* orthogonal complement of x */
/* output */
      *outx++ = x;
      *outy++ = y;
      *outz++ = z;
/* calculate fixed point location (c, 0, -c) */
      c = _fixedpoint(x, o);
/* inverse base transform */
      state[0] = x + z;
      state[1] = 2.0f * y;
      state[2] = -x + z;
/* update transformed linear system around unstable fixed point */
      state[0] -= c;
      state[2] += c;
      vcmul2(state, pole);
      state[2] *= r2;
      state[0] += c;
      state[2] -= c;
}
ctl->c_x = state[0];
  ctl->c_y = state[1];
  ctl->c_z = state[2];
return (w+10);
}
static void scrollgrid1D_dsp(t_scrollgrid1D *x, t_signal **sp)
{
  int n = sp[0]->s_n;
  int k;
dsp_add(scrollgrid1D_perform, 
      9, 
      &x->x_ctl, 
      sp[0]->s_n, 
      sp[0]->s_vec, 
      sp[1]->s_vec, 
      sp[2]->s_vec,
      sp[3]->s_vec,
      sp[4]->s_vec,
      sp[5]->s_vec,
      sp[6]->s_vec);
}                                  
static void scrollgrid1D_free(t_scrollgrid1D *x)
{
}
static void scrollgrid1D_reset(t_scrollgrid1D *x)
{
    x->x_ctl.c_x = 1;
    x->x_ctl.c_y = 1;
    x->x_ctl.c_z = 1;
}
t_class *scrollgrid1D_class;
static void *scrollgrid1D_new(t_floatarg algotype)
{
    t_scrollgrid1D *x = (t_scrollgrid1D *)pd_new(scrollgrid1D_class);
/* ins */
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));  
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));  
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));
/* outs */
    outlet_new(&x->x_obj, gensym("signal")); 
    outlet_new(&x->x_obj, gensym("signal")); 
    outlet_new(&x->x_obj, gensym("signal"));
/* init data */
    scrollgrid1D_reset(x);
return (void *)x;
}
void scrollgrid1D_tilde_setup(void)
{
    //post("scrollgrid1D~ v0.1");
    scrollgrid1D_class = class_new(gensym("scrollgrid1D~"), (t_newmethod)scrollgrid1D_new,
    	(t_method)scrollgrid1D_free, sizeof(t_scrollgrid1D), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(scrollgrid1D_class, t_scrollgrid1D, x_f);
    class_addmethod(scrollgrid1D_class, (t_method)scrollgrid1D_dsp, gensym("dsp"), 0); 
    class_addmethod(scrollgrid1D_class, (t_method)scrollgrid1D_reset, gensym("reset"), 0);
}
--- NEW FILE: blocknorm~.c ---
/*
 *   blocknorm.c  - normalize an array of dsp blocks (for spectral processing) 
 *   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 <stdlib.h>
#include <stdio.h>
#include <math.h>
#define MAXCHANNELS 32
typedef struct blocknormctl
{
    t_int c_channels;
    t_float **c_input;
    t_float **c_output;
} t_blocknormctl;
typedef struct blocknorm
{
  t_object x_obj;
  t_float x_f;
  t_blocknormctl x_ctl;
} t_blocknorm;
static t_int *blocknorm_perform(t_int *word)
{
t_blocknormctl *ctl  = (t_blocknormctl *)(word[1]);
    t_int n             = (t_int)(word[2]);
    t_float **in        = ctl->c_input;
    t_float **out       = ctl->c_output;
    t_int c             = ctl->c_channels;
    t_int i,j;
t_float p = 0.0f;
    t_float x, s;
/* get power */
    for (j=0;j<c;j++){
    for (i=0;i<n;i++){
        x = in[j][i];
        p += x*x;
    }
    }
/* compute normalization */
    if (p == 0.0f) s = 1.0f;
    else s =sqrt(((float)(c * n)) /  p);
/* normalize */
    for (j=0;j<c;j++){
    for (i=0;i<n;i++){
        out[j][i] *= s;
    }
    }
return (word+3);
}
static void blocknorm_dsp(t_blocknorm *x, t_signal **sp)
{
int i;
    int c = x->x_ctl.c_channels;
    for (i=0;i<c;i++){
    x->x_ctl.c_input[i] = sp[i]->s_vec;
    x->x_ctl.c_output[i] = sp[c+i]->s_vec;
    }
    dsp_add(blocknorm_perform, 2, &x->x_ctl, sp[0]->s_n);
}
static void blocknorm_free(t_blocknorm *x)
{
    free (x->x_ctl.c_output);
    free (x->x_ctl.c_input);
}
t_class *blocknorm_class;
static void *blocknorm_new(t_floatarg channels)
{
    int i = (int)channels;
    int j;
    t_blocknorm *x = (t_blocknorm *)pd_new(blocknorm_class);
if (i<1) i = 1;
    if (i>MAXCHANNELS) i = MAXCHANNELS;
    x->x_ctl.c_channels = i;
    x->x_ctl.c_input = malloc(sizeof(float)*i);
    x->x_ctl.c_output = malloc(sizeof(float)*i);
j = i;
    while (--j) inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));  
    while (i--) outlet_new(&x->x_obj, gensym("signal"));
return (void *)x;
}
void blocknorm_tilde_setup(void)
{
    blocknorm_class = class_new(gensym("blocknorm~"), (t_newmethod)blocknorm_new,
    	(t_method)blocknorm_free, sizeof(t_blocknorm), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(blocknorm_class, t_blocknorm, x_f); 
    class_addmethod(blocknorm_class, (t_method)blocknorm_dsp, gensym("dsp"), 0); 
}
--- NEW FILE: fdn~.c ---
/*
 *   fdn.c  - a feedback delay network (reverb tail)
 *   using a housholder reflection feedback matrix (In - 2/n 11T)
 *   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.
 */
/* TODO: CLEAN UP THIS MESS
add delay time generation code
add prime calculation routine (for prime delay line lengths)
add more diffuse feedback matrix (hadamard)
check filtering code
*/
#include <m_pd.h>
#include <math.h>
#include "extlib_util.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#define FDN_MIN_DECAY_TIME .01f
/*
#define NBPRIMES
int prime[NBPRIMES];
static int isprime(int n)
{
  int i=1;
  int d,m,p;
  while(1){
    p = prime[i++];
    m = n % p;
    if (m == 0) return 0; // it is a prime
    d = n / p;
    if (d < p) return 1;  // it is not a prime
  }
}
static int initprimes(void)
{
  int i, curprime;
  prime[0] = 1;
  prime[1] = 2;
  curprime = 3;
for(i=2; i<NBPRIMES; i++){
    while (!isprime(curprime)) curprime++;
    prime[i] = curprime;
    //printf("%d, ", curprime);
    curprime++;
  }
  printf("\n");
  return 0;
}
*/
//static int find_nearest_prime(int n){ return n;}
typedef struct fdnctl
{
  t_int c_order;      /* veelvoud van 4 */
  t_int c_maxorder;
  t_float c_leak;
  t_float c_input;
  t_float c_output;
  t_float *c_buf;
  t_float *c_gain_in;
  t_float *c_gain_state;
  t_float c_timehigh;
  t_float c_timelow;
  t_int *c_tap;       /* cirular feed: N+1 pointers: 1 read, (N-1)r/w, 1 write */ 
  t_float *c_length;  /* delay lengths in ms */
  t_int c_bufsize;
  t_float c_fsample;
  t_float *c_vector[2];
  t_float *c_vectorbuffer;
  t_int c_curvector;
} t_fdnctl;
typedef struct fdn
{
  t_object x_obj;
  t_float x_f;
  t_fdnctl x_ctl;
} t_fdn;
static void fdn_order(t_fdn *x, t_int order){
  if (order > x->x_ctl.c_maxorder) {
    post("fdn: this should not happen (panic!) order %d is larger than maxorder %d:", 
     order, x->x_ctl.c_maxorder ); 
    exit(1);
  }
  x->x_ctl.c_order = order;
  x->x_ctl.c_leak = -2./ order;
  x->x_ctl.c_input = 1./ sqrt(order); //????????????????
}
static void fdn_print(t_fdn *x)
{
  int i;
  fprintf(stderr, "fdn: delay coefficients (ms)\n");
  for (i=0;i<x->x_ctl.c_order;i++) {
    fprintf(stderr, "%f ", x->x_ctl.c_length[i]);
  }
  fprintf(stderr, "\n");
}
static void fdn_reset(t_fdn *x)
{
  int i;
  if (x->x_ctl.c_buf) memset(x->x_ctl.c_buf, 0, x->x_ctl.c_bufsize * sizeof(float));
  if (x->x_ctl.c_vectorbuffer) memset(x->x_ctl.c_vectorbuffer, 0, x->x_ctl.c_maxorder * 2 * sizeof(float));
}
static t_int *fdn_perform(t_int *w)
{
t_float *in      = (float *)(w[3]);
  t_float *outr    = (float *)(w[4]);
  t_float *outl    = (float *)(w[5]);
  t_fdnctl *ctl    = (t_fdnctl *)(w[1]);
  t_int n          = (t_int)(w[2]);
  t_float input    = ctl->c_input;
  t_float output   = ctl->c_output;
  t_float *gain_in    = ctl->c_gain_in;
  t_float *gain_state = ctl->c_gain_state;
  t_int order      = ctl->c_order;
  t_int *tap       = ctl->c_tap;
  t_float *buf     = ctl->c_buf;
  t_int mask       = ctl->c_bufsize - 1;
t_int i,j;
  t_float x,y,v,left,right,z;
  t_float filt_in, filt_last;
t_float *cvec, *lvec;
t_float save;
for(i=0;i<n;i++){
x = *in++;
    y = 0;
    left = 0;
    right = 0;
/* get temporary vector buffers */
    cvec = ctl->c_vector[ctl->c_curvector];
    lvec = ctl->c_vector[ctl->c_curvector ^ 1];
    ctl->c_curvector ^= 1;
/* read input vector + get sum and left/right output*/
    for(j=0;j<order;) 
      {
    z = buf[tap[j]];
    cvec[j] = z;
    y     += z;
    left  += z;
    right += z;
    j++;
z = buf[tap[j]];
    cvec[j] = z;
    y     += z;
    left  -= z;
    right += z;
    j++;
z = buf[tap[j]];
    cvec[j] = z;
    y     += z;
    left  += z;
    right -= z;
    j++;
z = buf[tap[j]];
    cvec[j] = z;
    y     += z;
    left  -= z;
    right -= z;
    j++;
      }
/* write output */
    *outl++ = left;
    *outr++ = right;
/* y == leak to all inputs */
    y *=  ctl->c_leak;
/* perform feedback */
    /* todo: decouple feedback & permutation */
    save = cvec[0];
    for (j=0; j<order-1; j++){
      cvec[j] = cvec[j+1] + y + x; 
    }
    cvec[order-1] = save + y + x;
/* apply gain + store result vector in delay lines + increment taps*/
    tap[0] = (tap[0]+1)&mask;
    for(j=0;j<order;j++) {
    save = gain_in[j] * cvec[j] + gain_state[j] * lvec[j];
    save = IS_DENORMAL(save) ? 0 : save;
    cvec[j] = save;
    buf[tap[j+1]] = save;
    tap[j+1] = (tap[j+1] + 1) & mask;
    }	
  }
return (w+6);
}
static void fdn_dsp(t_fdn *x, t_signal **sp)
{
x->x_ctl.c_fsample = sp[0]->s_sr;
  dsp_add(fdn_perform, 
      5, 
      &x->x_ctl, 
      sp[0]->s_n, 
      sp[0]->s_vec, 
      sp[1]->s_vec, 
      sp[2]->s_vec);
}                                  
static void fdn_free(t_fdn *x)
{
  if ( x->x_ctl.c_tap) free( x->x_ctl.c_tap);
  if ( x->x_ctl.c_length) free( x->x_ctl.c_length);
  if ( x->x_ctl.c_gain_in) free( x->x_ctl.c_gain_in);
  if ( x->x_ctl.c_gain_state) free( x->x_ctl.c_gain_state);
  if ( x->x_ctl.c_buf) free ( x->x_ctl.c_buf);
  if ( x->x_ctl.c_vectorbuffer) free ( x->x_ctl.c_vectorbuffer );
}
/*
each delay line is filtered with a first order iir filter:
(gl: dc gain, gh: ny gain)
H(z) = 2 gl gh / (gl + gh - z^-1 (gl - gh))
this results in the difference equation
yk = (2 gl gh ) / (gl + gh) x + (gl - gh) / (gl + gh) yk-1
*/
static void fdn_time(t_fdn *x, t_float timelow, t_float timehigh){
  t_float elow, ehigh;
  t_int i;
  t_float gainlow, gainhigh, gainscale;
if (timelow < FDN_MIN_DECAY_TIME) timelow = FDN_MIN_DECAY_TIME;
  if (timehigh < FDN_MIN_DECAY_TIME) timehigh = FDN_MIN_DECAY_TIME;
elow = -.003 / (timelow);
  ehigh = -.003 / (timehigh);
/* setup gains */
  for(i=0;i<x->x_ctl.c_order;i++){
    gainlow = pow(10, elow * (x->x_ctl.c_length[i]));
    gainhigh = pow(10, ehigh * (x->x_ctl.c_length[i]));
    gainscale = 1.0f / (gainlow + gainhigh);
    x->x_ctl.c_gain_in[i]    = 2.0f * gainlow * gainhigh * gainscale;
    x->x_ctl.c_gain_state[i] = (gainlow - gainhigh) * gainscale;
  }
  x->x_ctl.c_timehigh = timehigh;
  x->x_ctl.c_timelow = timelow;
}
static void fdn_updatedamping(t_fdn *x)
{
  fdn_time(x, x->x_ctl.c_timelow, x->x_ctl.c_timehigh);
}
static void fdn_timelow(t_fdn *x, t_float f){
  x->x_ctl.c_timelow = fabs(f);
  fdn_updatedamping(x);
}
static void fdn_timehigh(t_fdn *x, t_float f){
  x->x_ctl.c_timehigh = fabs(f);
  fdn_updatedamping(x);
}
static void fdn_setupdelayline(t_fdn *x){
  int sum, t, n;
  int mask = x->x_ctl.c_bufsize - 1;
  int start =  x->x_ctl.c_tap[0];
  t_int *tap = x->x_ctl.c_tap;
  float *length = x->x_ctl.c_length;
  float scale = sys_getsr() * .001f;
sum = 0;
  tap[0] = (start & mask);
  for (t=1; t<= x->x_ctl.c_order; t++){
    sum += (int)(length[t-1] * scale);
    tap[t]=(start+sum)&mask;
  }
if (sum > mask){ 
    post("fdn: warning: not enough delay memory, behaviour is undefined (this could lead to instability...)");
  }
}
static void fdn_list (t_fdn *x,  t_symbol *s, int argc, t_atom *argv){
  int i;
  float l;
  int sum=0;
int order = argc & 0xfffffffc;
if (order < 4) return;
  if (order > x->x_ctl.c_maxorder) return;
fdn_order(x, order);
  for(i=0; i<order; i++)
      if (argv[i].a_type == A_FLOAT) x->x_ctl.c_length[i] = argv[i].a_w.w_float;
fdn_setupdelayline(x);
  fdn_updatedamping(x);
}
static void fdn_linear(t_fdn *x, t_float forder, t_float min, t_float max)
{
    t_int order = ((int)forder) & 0xfffffffc;
    t_float length, inc;
    t_int i;
if (order < 4) return;
    if (order > x->x_ctl.c_maxorder) return;
    if (min <= 0) return;
    if (max <= 0) return;
inc = (max - min) / (float)(order - 1);
    length = min;
for (i=0; i<order; i++){
    x->x_ctl.c_length[i] = length;
    length += inc;
    }
fdn_order(x, order);
    fdn_setupdelayline(x);
    fdn_updatedamping(x);
}
static void fdn_exponential(t_fdn *x, t_float forder, t_float min, t_float max)
{
    t_int order = ((int)forder) & 0xfffffffc;
    t_float length, inc;
    t_int i;
if (order < 4) return;
    if (order > x->x_ctl.c_maxorder) return;
    if (min <= 0) return;
    if (max <= 0) return;
inc = pow (max / min, 1.0f / ((float)(order - 1)));
    length = min;
for (i=0; i<order; i++){
    x->x_ctl.c_length[i] = length;
    length *= inc;
    }
fdn_order(x, order);
    fdn_setupdelayline(x);
    fdn_updatedamping(x);
}
t_class *fdn_class;
static void *fdn_new(t_floatarg maxiorder, t_floatarg maxibufsize)
{
  t_int order = maxiorder;
  t_int bufround;
  t_fdn *x = (t_fdn *)pd_new(fdn_class);
  t_float scale = sys_getsr() * .001f;
  t_int bufsize = (t_int)(scale * maxibufsize);
inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("timelow"));  
  inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("timehigh"));  
  outlet_new(&x->x_obj, gensym("signal")); 
  outlet_new(&x->x_obj, gensym("signal"));
/* init data */
  if (order < 4) order = 8;
  if (bufsize < 64) bufsize = 65536;
bufround = 1;
  while (bufround < bufsize) bufround *= 2;
  bufsize = bufround;
post("fdn: maximum nb of delay lines %d, total buffer size %d samples (%f seconds)", 
       order, bufsize, ((float)bufsize) / sys_getsr());
x->x_ctl.c_maxorder = order;
  x->x_ctl.c_buf = (float *)malloc(sizeof(float) * bufsize);
  x->x_ctl.c_bufsize = bufsize;
  x->x_ctl.c_fsample = sys_getsr();
  x->x_ctl.c_tap = (t_int *)malloc((order + 1) * sizeof(t_int)); 
  x->x_ctl.c_length = (t_float *)malloc(order * sizeof(t_int)); 
  x->x_ctl.c_gain_in = (t_float *)malloc(order * sizeof(t_float));
  x->x_ctl.c_gain_state = (t_float *)malloc(order * sizeof(t_float));
  x->x_ctl.c_vectorbuffer = (t_float *)malloc(order * 2 * sizeof(float));
  memset(x->x_ctl.c_vectorbuffer, 0, order * 2 * sizeof(float));
  x->x_ctl.c_curvector = 0;
  x->x_ctl.c_vector[0] = &x->x_ctl.c_vectorbuffer[0];
  x->x_ctl.c_vector[1] = &x->x_ctl.c_vectorbuffer[order];
/* preset */
  fdn_order(x,8);
  x->x_ctl.c_length[0]= 29.0f;
  x->x_ctl.c_length[1]= 31.0f;
  x->x_ctl.c_length[2]= 37.0f;
  x->x_ctl.c_length[3]= 67.0f;
  x->x_ctl.c_length[4]= 82.0f;
  x->x_ctl.c_length[5]= 110.0f;
  x->x_ctl.c_length[6]= 172.0f;
  x->x_ctl.c_length[7]= 211.0f;
  fdn_setupdelayline(x);
  fdn_time(x, 4, 1);
/* reset delay memory to zero */
  fdn_reset(x);
return (void *)x;
}
void fdn_tilde_setup(void)
{
  //post("fdn~ v0.1");
    fdn_class = class_new(gensym("fdn~"), (t_newmethod)fdn_new,
    	(t_method)fdn_free, sizeof(t_fdn), 0, A_DEFFLOAT, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(fdn_class, t_fdn, x_f); 
    class_addmethod(fdn_class, (t_method)fdn_print, gensym("print"), 0);
    class_addmethod(fdn_class, (t_method)fdn_reset, gensym("reset"), 0);
    class_addmethod(fdn_class, (t_method)fdn_timehigh, gensym("timehigh"), A_DEFFLOAT, 0);
    class_addmethod(fdn_class, (t_method)fdn_timelow, gensym("timelow"), A_DEFFLOAT, 0);
    class_addmethod(fdn_class, (t_method)fdn_list, gensym("lines"), A_GIMME, 0);
    class_addmethod(fdn_class, (t_method)fdn_dsp, gensym("dsp"), 0); 
    class_addmethod(fdn_class, (t_method)fdn_linear, gensym("linear"), A_FLOAT, A_FLOAT, A_FLOAT, 0); 
    class_addmethod(fdn_class, (t_method)fdn_exponential, gensym("exponential"), A_FLOAT, A_FLOAT, A_FLOAT, 0);
}
--- NEW FILE: dynwav~.c ---
/*
 *   blosc.c  - bandlimited oscillators
 *   data organization is in (real, imag) pairs
 *   the first 2 components are (DC, NY)
 *   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>
#include <string.h>
#define MAXORDER 1024
typedef struct dynwavctl
{
  t_float *c_buf1; /* current */
  t_float *c_buf2; /* old */
  t_int c_order;
} t_dynwavctl;
typedef struct dynwav
{
  t_object x_obj;
  t_float x_f;
  t_dynwavctl x_ctl;
} t_dynwav;
static t_int *dynwav_perform(t_int *w)
{
t_float *wave     = (float *)(w[3]);
  t_float *freq     = (float *)(w[4]);
  t_float *out      = (float *)(w[5]);
  t_dynwavctl *ctl  = (t_dynwavctl *)(w[1]);
  t_int n           = (t_int)(w[2]);
t_float *buf, *dbuf, *swap;
int i;
  int mask = n-1;
/* swap buffer pointers */
  swap = ctl->c_buf1;               /* this is the last one stored */
  buf  = ctl->c_buf1 = ctl->c_buf2; /* put oldest in newest to overwrite */
  dbuf = ctl->c_buf2 = swap;        /* put last one in oldest */
if (buf && dbuf)
    {
/* store input wavetable in buffer */
      memcpy(buf, wave, n*sizeof(t_float));
for (i = 0; i < n; i++)
    {
      float findex = *freq++ * (t_float)n;
      int index = findex;
      float frac,  a,  b,  c,  d, cminusb, q, r;
      int ia, ib, ic, id;
      
      frac = findex - index;
ia = (index-1) & mask;
      ib = (index  ) & mask;
      ic = (index+1) & mask;
      id = (index+2) & mask;
q = i+1;
      q /= n;
r = n-1-i;
      r /= n;
      
      /* get 4 points, wrap index */
      a = q * buf[ia] + r * dbuf[ia];
      b = q * buf[ib] + r * dbuf[ib];
      c = q * buf[ic] + r * dbuf[ic];
      d = q * buf[id] + r * dbuf[id];
      
      cminusb = c-b;
      *out++ = b + frac * (cminusb - 0.5f * (frac-1.) * 
    		       ((a - d + 3.0f * cminusb) * frac + 
    			(b - a - cminusb)));
    }
}
  return (w+6);
}
static t_int *dynwav_perform_8point(t_int *w) /* werkt nog nie tegoei */
{
t_float *wave     = (float *)(w[3]);
  t_float *freq     = (float *)(w[4]);
  t_float *out      = (float *)(w[5]);
  t_dynwavctl *ctl  = (t_dynwavctl *)(w[1]);
  t_int n           = (t_int)(w[2]);
t_float *buf, *dbuf, *swap;
int i;
  int mask = n-1;
/* swap buffer pointers */
  swap = ctl->c_buf1;               /* this is the last one stored */
  buf  = ctl->c_buf1 = ctl->c_buf2; /* put oldest in newest to overwrite */
  dbuf = ctl->c_buf2 = swap;        /* put last one in oldest */
if (buf && dbuf)
    {
/* const float N1 = 1 / (   2      * (1-(1/9))  * (1-(1/25))  * (1-(1/49))  );
      ** const float N2 = 1 / ( (1-(9))  *      2     * (1-(9/25))  * (1-(9/49))  );
      ** const float N3 = 1 / ( (1-(25)) * (1-(25/9)) *      2      * (1-(25/49)) );
      ** const float N4 = 1 / ( (1-(49)) * (1-(49/9)) * (1-(49/25)) *    2        );
      */
const float N1 =   0.59814453125;
      const float N2 =  -0.11962890625;
      const float N3 =   0.02392578125;
      const float N4 =  -0.00244140625;
/* store input wavetable in buffer */
      memcpy(buf, wave, n*sizeof(t_float));
for (i = 0; i < n; i++)
    {
      float findex = *freq++ * (t_float)n;
      int index = findex;
      float frac, q, r, fm, fp, fe, fo;
      float x1, x2, x3, x4;
      float g1, g2, g3, g4;
      float gg, g2g3g4, g1g3g4, g1g2g4, g1g2g3;
      float acc;
      int im, ip;
      
      frac = 2 *(findex - index) - 1;
x1 = frac;
      x2 = frac/3;
      x3 = frac/5;
      x4 = frac/7;
g1 = 1 - x1*x1;
      g2 = 1 - x2*x2;
      g3 = 1 - x3*x3;
      g4 = 1 - x4*x4;
gg       = g3 * g4;
      g2g3g4   = g2 * gg;       /* 1 */
      g1g3g4   = g1 * gg;       /* 2 */
      gg       = g1 * g2;
      g1g2g4   = g4 * gg;       /* 3 */
      g1g2g3   = g3 * gg;       /* 4 */
/* triangle interpolation between current and past wavetable*/
      q = i+1;
      q /= n;
r = n-1-i;
      r /= n;
/* 1, -1*/
      im =  (index  ) & mask;
      ip =  (index+1) & mask;
      fm =  q * buf[im] + r * dbuf[im];
      fp =  q * buf[ip] + r * dbuf[ip];
      fe = fp + fm;
      fo = fp - fm;
acc = N1 * g2g3g4 * (fe + x1*fo);
/* 2, -2 */
      im =  (index-1) & mask;
      ip =  (index+2) & mask;
      fm =  q * buf[im] + r * dbuf[im];
      fp =  q * buf[ip] + r * dbuf[ip];
      fe = fp + fm;
      fo = fp - fm;
acc += N2 * g1g3g4 * (fe + x2*fo);
/* 3, -3 */
      im =  (index-2) & mask;
      ip =  (index+3) & mask;
      fm =  q * buf[im] + r * dbuf[im];
      fp =  q * buf[ip] + r * dbuf[ip];
      fe = fp + fm;
      fo = fp - fm;
acc += N3 * g1g2g4 * (fe + x3*fo);
/* 4, -4 */
      im =  (index-3) & mask;
      ip =  (index+4) & mask;
      fm =  q * buf[im] + r * dbuf[im];
      fp =  q * buf[ip] + r * dbuf[ip];
      fe = fp + fm;
      fo = fp - fm;
acc += N4 * g1g2g3 * (fe + x4*fo);
*out++ = acc;
}
}
  return (w+6);
}
static void dynwav_dsp(t_dynwav *x, t_signal **sp)
{
  int n = sp[0]->s_n;
  int k;
if (x->x_ctl.c_order != n)
    {
      if (x->x_ctl.c_buf1) free (x->x_ctl.c_buf1);
      if (x->x_ctl.c_buf2) free (x->x_ctl.c_buf2);
x->x_ctl.c_buf1 = (t_float *)malloc(n*sizeof(t_float));
      x->x_ctl.c_buf2 = (t_float *)malloc(n*sizeof(t_float));
for(k=0; k<n; k++)
    {
      x->x_ctl.c_buf1[k] = 0;
      x->x_ctl.c_buf2[k] = 0;
    }
x->x_ctl.c_order = n;
    }
dsp_add(dynwav_perform_8point, 5, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec, sp[2]->s_vec);
}                                  
static void dynwav_free(t_dynwav *x)
{
if (x->x_ctl.c_buf1) free (x->x_ctl.c_buf1);
  if (x->x_ctl.c_buf2) free (x->x_ctl.c_buf2);
}
t_class *dynwav_class;
static void *dynwav_new(t_floatarg order)
{
    t_dynwav *x = (t_dynwav *)pd_new(dynwav_class);
    int iorder = (int)order;
    int i, n=64, k;
/* in 2 */
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("signal"), gensym("signal"));
/* out 1 */
    outlet_new(&x->x_obj, gensym("signal"));
/* init data */
    x->x_ctl.c_buf1 = (t_float *)malloc(n*sizeof(t_float));
    x->x_ctl.c_buf2 = (t_float *)malloc(n*sizeof(t_float));
for(k=0; k<n; k++)
      {
    x->x_ctl.c_buf1[k] = 0;
    x->x_ctl.c_buf2[k] = 0;
      }
x->x_ctl.c_order = n;
return (void *)x;
}
void dynwav_tilde_setup(void)
{
  //post("dynwav~ v0.1");
    dynwav_class = class_new(gensym("dynwav~"), (t_newmethod)dynwav_new,
    	(t_method)dynwav_free, sizeof(t_dynwav), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(dynwav_class, t_dynwav, x_f);
    class_addmethod(dynwav_class, (t_method)dynwav_dsp, gensym("dsp"), 0);
}
--- NEW FILE: cheby~.c ---
/*
 *   cheby.c  - chebyshev polynomial evaluation 
 *   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>
#define MAX_ORDER 1024
#define DEFAULT_ORDER 4
typedef struct chebyctl
{
  t_float c_gain[MAX_ORDER];
  t_int c_order;
} t_chebyctl;
typedef struct cheby
{
  t_object x_obj;
  t_float x_f;
  t_chebyctl x_ctl;
} t_cheby;
static void cheby_bang(t_cheby *x)
{
}
static void cheby_coef(t_cheby *x, t_floatarg coef, t_floatarg f)
{
  int i = (int)coef;
  if ((i > 0) && (i < x->x_ctl.c_order + 1)){
    x->x_ctl.c_gain[i-1] = f;
  /*   post("cheby: harmonic %d set to %f", i, f); */
  }
}
static t_int *cheby_perform(t_int *w)
{
t_float *in    = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_chebyctl *ctl  = (t_chebyctl *)(w[1]);
  t_float *gain  = ctl->c_gain;
  t_int i;
  t_int n = (t_int)(w[2]), k;
  t_float x,y,t1,t2,t,acc;
for (i = 0; i < n; i++) 
    {
      x = *in++;
gain = ctl->c_gain;
      t2 = 1;                      /* T_0 */
      t1 = x;                      /* T_1 */
acc = *gain++ * x;           /* a_1 T_1 */
      for (k=2; k<=ctl->c_order; k++){
    t = 2*x*t1 - t2;           /* T_k = 2 x T_{k-1} - T_{k-2} */
    acc += *gain++ * t;        /* a_k T_k */
    t2 = t1;
    t1 = t;
      }
*out++ = acc;
}
return (w+5);
}
static void cheby_dsp(t_cheby *x, t_signal **sp)
{
    dsp_add(cheby_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
}                                  
static void cheby_free(void)
{
}
t_class *cheby_class;
static void *cheby_new(t_floatarg order_f)
{
  int i;
  int order = (int)order_f;
t_cheby *x = (t_cheby *)pd_new(cheby_class);
    outlet_new(&x->x_obj, gensym("signal"));
if (order < 1)         order = DEFAULT_ORDER;  /* default */
    if (order > MAX_ORDER) order = MAX_ORDER;      /* maximum */
//post("cheby: order = %d", order);
x->x_ctl.c_order = order;
    cheby_coef(x, 1, 1);
    for (i=2; i<order+1; i++){
      cheby_coef(x, 0, i);
    }
return (void *)x;
}
void cheby_tilde_setup(void)
{
  //post("cheby~ v0.1");
    cheby_class = class_new(gensym("cheby~"), (t_newmethod)cheby_new,
    	(t_method)cheby_free, sizeof(t_cheby), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(cheby_class, t_cheby, x_f); 
    class_addmethod(cheby_class, (t_method)cheby_bang, gensym("bang"), 0);
    class_addmethod(cheby_class, (t_method)cheby_dsp, gensym("dsp"), 0); 
    class_addmethod(cheby_class, (t_method)cheby_coef, gensym("coef"), A_DEFFLOAT, A_DEFFLOAT, 0);
}
--- eblosc.c DELETED ---
--- NEW FILE: matrix~.c ---
/*
 *   matrix.c  - applies a matrix transform to a signal block
 *   intended for spectral processing, dynwav
 *
 *   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>
#include <string.h>
#define MAXORDER 1024
typedef struct matrixctl
{
  t_float *c_A; /* matrix */
  t_float *c_x; /* vector */
  t_int c_order;
} t_matrixctl;
typedef struct matrix
{
  t_object x_obj;
  t_float x_f;
  t_matrixctl x_ctl;
} t_matrix;
static void matrix_load(t_matrix *x, t_symbol *s)
{
  FILE *matrix;
if(s && s->s_name)
    {
      post("matrix: loading %s",s->s_name);
      if(matrix = fopen(s->s_name, "r"))
    {
      int n = x->x_ctl.c_order;
      fread(x->x_ctl.c_A, sizeof(float), n*n, matrix);
    }
      else post("matrix: error, cant open file.");
    }
}
static t_int *matrix_perform(t_int *w)
{
t_float *in       = (float *)(w[3]);
  t_float *out      = (float *)(w[4]);
  t_matrixctl *ctl  = (t_matrixctl *)(w[1]);
  t_int n           = (t_int)(w[2]);
t_int i,j;
  t_float *A = ctl->c_A;
  t_float *x = ctl->c_x;
if (in == out) /* store input if ness. */
    {
      memcpy(x, in, sizeof(t_float)*n);
      in = x;
    }
  bzero(out, sizeof(t_float)*n); /* init output */
for (j=0; j<n; j++)
      for (i=0; i<n; i++) 
    out[i] += A[i+j*n] * in[j];
return (w+5);
}
static void matrix_dsp(t_matrix *x, t_signal **sp)
{
  int n = sp[0]->s_n;
  int k,i;
if (x->x_ctl.c_order != n)
    {
      if (x->x_ctl.c_A) free (x->x_ctl.c_A);
x->x_ctl.c_A = (t_float *)calloc(n*n,sizeof(t_float));
      x->x_ctl.c_x = (t_float *)calloc(n,sizeof(t_float));
      x->x_ctl.c_order = n;
    }
for (i=0;i<n;i++) x->x_ctl.c_A[i] = 1;
dsp_add(matrix_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
}                                  
static void matrix_free(t_matrix *x)
{
if (x->x_ctl.c_A) free (x->x_ctl.c_A);
  if (x->x_ctl.c_x) free (x->x_ctl.c_x);
}
t_class *matrix_class;
static void *matrix_new(t_floatarg order)
{
    t_matrix *x = (t_matrix *)pd_new(matrix_class);
    int iorder = (int)order;
    int i, n=64, k;
/* out 1 */
    outlet_new(&x->x_obj, gensym("signal"));
/* init data */
x->x_ctl.c_A = (t_float *)calloc(n*n,sizeof(t_float));
    x->x_ctl.c_x = (t_float *)calloc(n,sizeof(t_float));
for (i=0;i<n;i++) x->x_ctl.c_A[i] = 1;
    x->x_ctl.c_order = n;
return (void *)x;
}
void bmatrix_tilde_setup(void)
{
  //post("matrix~ v0.1");
    matrix_class = class_new(gensym("bmatrix~"), (t_newmethod)matrix_new,
    	(t_method)matrix_free, sizeof(t_matrix), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(matrix_class, t_matrix, x_f);
    class_addmethod(matrix_class, (t_method)matrix_dsp, gensym("dsp"), 0); 
    class_addmethod(matrix_class, (t_method)matrix_load, gensym("load"), A_SYMBOL,0);
}
--- NEW FILE: permut~.c ---
/*
 *   permut.c  - applies a (random) permutation on a signal block
 *   intended for spectral processing, dynwav
 *
 *   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 "extlib_util.h"
#include <stdlib.h>
typedef struct permutctl
{
  char c_type;
  t_int *c_permutationtable;
  int c_blocksize;
} t_permutctl;
typedef struct permut
{
  t_object x_obj;
  t_float x_f;
  t_permutctl x_ctl;
} t_permut;
static inline void permut_perform_permutation(t_float *S, int n, t_int *f)
{
  t_int k,l;
  t_float swap;
  for(k=0; k<n; k++)
    {
      l = f[k];
      while (l<k) l = f[l];
      swap = S[k];
      S[k] = S[l];
      S[l] = swap;
    }
}
static void permut_random(t_permut *x, t_floatarg seed)
{
  int i,j;
  int N = x->x_ctl.c_blocksize;
  int mask = N-1;
  t_int *p = x->x_ctl.c_permutationtable;
  int r, last = 0;
//srand(* ((unsigned int *)(&seed)));
  srand (((t_flint)seed).i);
if(p)
    {
      p[0] = rand() & mask;
      for (i=1;i<N;i++)
    {
      r = rand() & mask;
      j = 0;
      while(j<i)
        {
          if (p[j] == r)
    	{
    	  r = (r + 1) & mask;
    	  j = 0;
    	}
          else j++;
        }
      p[i] = r;
      
    }
    }
}
static void permut_bang(t_permut *x)
{
    unsigned int r = rand();
    //permut_random(x, *((float *)(&r)));
    permut_random(x, ((t_flint)r).f);
}
static void permut_resize_table(t_permut *x, int size)
{
  if (x->x_ctl.c_blocksize != size)
    {
      if (x->x_ctl.c_permutationtable)
    free(x->x_ctl.c_permutationtable);
      x->x_ctl.c_permutationtable = (t_int *)malloc(sizeof(int)*size);
      x->x_ctl.c_blocksize = size;
/* make sure it's initialized */
      permut_bang(x);
} 
}
static t_int *permut_perform(t_int *w)
{
t_float *in    = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_permutctl *ctl  = (t_permutctl *)(w[1]);
  t_int i;
  t_int n = (t_int)(w[2]);
  t_float x,y;
  t_int *p =  ctl->c_permutationtable;
if (in != out)
    for (i=0; i<n; i++) out[i] = in[i];
permut_perform_permutation(out, n, p);
return (w+5);
}
static void permut_dsp(t_permut *x, t_signal **sp)
{
  permut_resize_table(x, sp[0]->s_n);
  dsp_add(permut_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
}
static void permut_free(t_permut *x)
{
  if (x->x_ctl.c_permutationtable) free(x->x_ctl.c_permutationtable);
}
t_class *permut_class;
static void *permut_new(void)
{
    t_permut *x = (t_permut *)pd_new(permut_class);
    outlet_new(&x->x_obj, gensym("signal"));
x->x_ctl.c_permutationtable = 0;
    x->x_ctl.c_blocksize = 0;
    permut_resize_table(x, 64);
    permut_random(x, 0);
return (void *)x;
}
void permut_tilde_setup(void)
{
  //post("permut~ v0.1");
    permut_class = class_new(gensym("permut~"), (t_newmethod)permut_new,
    	(t_method)permut_free, sizeof(t_permut), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(permut_class, t_permut, x_f); 
    class_addmethod(permut_class, (t_method)permut_random, gensym("random"), A_FLOAT, 0);
    class_addmethod(permut_class, (t_method)permut_bang, gensym("bang"), 0);
    class_addmethod(permut_class, (t_method)permut_dsp, gensym("dsp"), 0);
}
--- NEW FILE: dist~.c ---
/*
 *   dist.c  -  wave shaping extern 
 *   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>
#define     CLIP          0
#define     INVERSE       1
#define     INVERSESQ     2
#define     INVERSECUB    3
#define     RAT1          4
#define     RAT2          5
#define     FULLRECT      6
#define     HALFRECT      7
#define     PULSE         8
#define     NEWTON1       9
#define     UPPERCLIP    10
#define     COMPARATOR   11
typedef struct distctl
{
  t_float c_gain;
  t_float c_delay;
  char c_type;
} t_distctl;
typedef struct dist
{
  t_object x_obj;
  t_float x_f;
  t_distctl x_ctl;
} t_dist;
void dist_bang(t_dist *x)
{
}
void dist_gain(t_dist *x,  t_floatarg f)
{
  x->x_ctl.c_gain = f;
}
static t_int *dist_perform(t_int *w)
{
t_float *in    = (float *)(w[3]);
  t_float *out    = (float *)(w[4]);
  t_distctl *ctl  = (t_distctl *)(w[1]);
  t_float gain  = ctl->c_gain;
  t_int i;
  t_int n = (t_int)(w[2]);
  t_float x,y,v;
  t_float z = ctl->c_delay;
switch(ctl->c_type){
  case CLIP:	
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
    x = (x >  1) ? ( 1.) : x;
    x = (x < -1) ? (-1.) : x;
    *out++ = 0.9999 * x;
}
    break;
case INVERSE:	
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
    x = (x > 1) ? (2. - 1/x) : x;
    x = (x < -1) ? (-2. - 1/x) : x;
    *out++ = x/2.0001;
}
    break;
case INVERSESQ:	
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
    x = (x > 1) ? (2. - 1/x) : x;
    x = (x < -1) ? (-2. - 1/x) : x;
    x /= 2;
    *out++ = 1.999*x*x-1;
}
    break;
case INVERSECUB:	
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
    x = (x > 1) ? (2. - 1/x) : x;
    x = (x < -1) ? (-2. - 1/x) : x;
    x /= 2;
    *out++ = .9999 * x*x*x;
}
    break;
case RAT1: /*(2*d./((1+(d).^2)))*/   
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
    y = (1. + x*x);
    x = 1.9999*x/y;
    *out++ = x;
      }
    break;
case RAT2: /*(2*d./((1+(d).^16)))*/   
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
    y = x*x;
    y *= y;
    y *= y;
    y *= y;
    y = (1. + y);
    x = 1.2*x/y;
    *out++ = x;
      }
    break;
case FULLRECT:
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
    x = (x>0) ? x : -x;
    x = (x>1) ? 1 : x;
    *out++ = 1.9999*(x-.5);
      }
    break;
case HALFRECT:
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
    x = (x>0) ? x : 0;
    x = (x>1) ? 1 : x;
    *out++ = 1.9999*(x-.5);
      }
    break;
case PULSE:
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
    y = (x>0) ? (1):(-1);
    x = (z*y > 0) ? (0) : (y);
    *out++ = .9999 * x;
    z = x;
    
      }
    ctl->c_delay = z;
    break;
case NEWTON1:
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
    y = 1./(1.+x*x);
z = .5;
    z = .5*(y/z + z);
    z = .5*(y/z + z);
    z = .5*(y/z + z);
/*	z = .5*(y/z + z);
     *	z = .5*(y/z + z);
     *	z = .5*(y/z + z);
     */
*out++ = x * z; 
    
      }
    ctl->c_delay = z;
    break;
case UPPERCLIP:
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
x = (x < 0.0f) ? 0.0f : x;
    x = (x > 0.9999f) ? 0.9999f : x;
*out++ = x; 
    
      }
    break;
case COMPARATOR:
    for (i = 0; i < n; i++)
      {
    x = *in++ * gain;
x = (x > 0.0f) ? 1.0f : -1.0f;
*out++ = x; 
    
      }
    break;
default:
for (i = 0; i < n; i++) *out++ = *in++;
    break;
}
return (w+5);
}
static void dist_dsp(t_dist *x, t_signal **sp)
{
    dsp_add(dist_perform, 4, &x->x_ctl, sp[0]->s_n, sp[0]->s_vec, sp[1]->s_vec);
}                                  
void dist_free(void)
{
}
t_class *dist_class;
void *dist_new(t_floatarg type)
{
    t_dist *x = (t_dist *)pd_new(dist_class);
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("gain"));  
    outlet_new(&x->x_obj, gensym("signal"));
dist_gain(x, 1);
    x->x_ctl.c_type = (char)type;
    x->x_ctl.c_delay = 0;
return (void *)x;
}
void dist_tilde_setup(void)
{
  //post("dist~ v0.1");
    dist_class = class_new(gensym("dist~"), (t_newmethod)dist_new,
    	(t_method)dist_free, sizeof(t_dist), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(dist_class, t_dist, x_f); 
    class_addmethod(dist_class, (t_method)dist_bang, gensym("bang"), 0);
    class_addmethod(dist_class, (t_method)dist_dsp, gensym("dsp"), 0); 
    class_addmethod(dist_class, (t_method)dist_gain, gensym("gain"), A_FLOAT, 0);
}
--- dist.c DELETED ---
--- bdiag.c DELETED ---
--- NEW FILE: window~.c ---
/*
 *   window.c  - window generation abstraction
 *   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 <stdlib.h>
static t_class *window_class;
typedef struct _window
{
    t_object x_obj;
    t_float x_f;
    t_float *x_window;
    t_int x_size;
    t_symbol *x_type;
    t_float x_typearg;
} t_window;
static t_int *window_perform(t_int *w)
{
    t_window *x = (t_window *)(w[1]);
    t_float *in = (t_float *)(w[2]);
    t_float *out = (t_float *)(w[3]);
    t_float *window = x->x_window;
    int n = (int)(w[4]);
    while (n--)
    {
    	*out++ = (*in++) * (*window++);
    }
    return (w+5);
}
static void window_size(t_window *x, t_int n)
{
    if (x->x_size != n){
    if (x->x_window) free(x->x_window);
    x->x_window = malloc(sizeof(float)*n);
    x->x_size = n;
    }
}
static void window_type(t_window *x, t_symbol *s, t_float f)
{
    int i;
    float a = 0;
    float a_inc = 2 * M_PI / (float)(x->x_size);
    if (!s) s = gensym("hamming");
    if (s == gensym("hamming")){
    for (i=0; i<x->x_size; i++){
        float c = cos(a);
        x->x_window[i] = (0.54 - 0.46 * c);
        a += a_inc;
    }
    }
    else if (s == gensym("hann")){
    for (i=0; i<x->x_size; i++){
        float c = cos(a);
        x->x_window[i] = (0.5 - 0.5 * c);
        a += a_inc;
    }
    }
    else if (s == gensym("hann/hamming")){
    for (i=0; i<x->x_size; i++) {
        float c = cos(a);
        x->x_window[i] = (0.5 - 0.5 * c) / (0.54 - 0.46 * c);
        a += a_inc;
    }
    }
    else if (s == gensym("bfft_pink")){
    x->x_window[0] = 1.0f; //DC
    x->x_window[1] = 1.0f / sqrt((double)(x->x_size>>1)); //NY
    for (i=2; i<x->x_size; i+=2) {
        double freq = (double)(i>>1);
        float amp = sqrt(1.0 / freq);
        x->x_window[i] = amp;
        x->x_window[i+1] = amp;
    }
    }   
    else if (s == gensym("bfft_blue")){
    x->x_window[0] = 1.0f; //DC
    x->x_window[1] = sqrt((double)(x->x_size>>1)); //NY
    for (i=2; i<x->x_size; i+=2) {
        double freq = (double)(i>>1);
        float amp = sqrt(freq);
        x->x_window[i] = amp;
        x->x_window[i+1] = amp;
    }
    }   
    else if (s == gensym("bfft_db/octave")){
    float power = f/6.0;
    x->x_window[0] = 1.0f; //DC
    x->x_window[1] = pow((double)(x->x_size>>1), power); //NY
    for (i=2; i<x->x_size; i+=2) {
        double freq = (double)(i>>1);
        float amp = pow(freq, power);
        x->x_window[i] = amp;
        x->x_window[i+1] = amp;
    }
    }
/* default is no window */
    else{
    post("bwin~: unknown window type, using rectangular");
    for (i=0; i<x->x_size; i++) x->x_window[i] = 1.0f;
    }
x->x_type = s;
    x->x_typearg = f;
}
static void window_dsp(t_window *x, t_signal **sp)
{
    int n =  sp[0]->s_n;
    if (x->x_size != n){
    window_size(x, n);
    window_type(x, x->x_type, x->x_typearg);
    }
dsp_add(window_perform, 4, x, sp[0]->s_vec, sp[1]->s_vec, n);
}
static void window_free(t_window *x)
{
    free(x->x_window);
}
static void *window_new(t_symbol *s)
{
    t_window *x = (t_window *)pd_new(window_class);
    outlet_new(&x->x_obj, &s_signal);
    x->x_window = 0;
    window_size(x, 64);
    window_type(x, s, 0);
    return (x);
}
void window_tilde_setup(void)
{
    window_class = class_new(gensym("bwin~"), (t_newmethod)window_new, (t_method)window_free,
    	sizeof(t_window), 0, A_DEFSYMBOL, A_NULL);
    CLASS_MAINSIGNALIN(window_class, t_window, x_f);
    class_addmethod(window_class, (t_method)window_dsp, gensym("dsp"), A_NULL);
    class_addmethod(window_class, (t_method)window_type, gensym("type"), A_SYMBOL, A_DEFFLOAT, A_NULL);
}
--- bitsplit.c DELETED ---
--- scrollgrid1D.c DELETED ---
--- window.c DELETED ---