fpu PD

fpu.c

#include <math.h>
#include <stdio.h>
#include "fpu.h"

static const double sr_default = 48000.0;
// time step of the simulation
static const double dt = 0.000001;
// iterations per sample
static const int ips = 32;

//---------------------
//---- static functions

// update forces
static void fpu_update_f(struct fpu* fpu) {
  // start and end positions are fixed at zero
  double f;
  double *x = fpu->x;
  double d0, d1;
  for(int i=1; i<FPU_NUM_MASSES-1; i++) {	
	d1 = x[i+1] - x[i];
	d0 = x[i] - x[i-1];
	f = d1 - d0 + fpu->epsilon * ((d1*d1*d1) - (d0*d0*d0));
	f *= fpu->beta;
	f -= fpu->rho * fpu->v[i];
	fpu->f[i] = f;
  } 
}

// update velocities
static void fpu_update_v(struct fpu* fpu) {
  // start and end positions are fixed at zero
  for(int i=1; i<FPU_NUM_MASSES-1; i++) {
	fpu->v[i] += fpu->f[i] * fpu->dt;
  }
}

// update positions
static void fpu_update_x(struct fpu* fpu) {
  // start and end positions are fixed at zero
  for(int i=1; i<FPU_NUM_MASSES-1; i++) {
	fpu->x[i] += fpu->v[i];
  }
}

//---------------------
//---- extern functions

void fpu_init(struct fpu* fpu) {
  fpu->dt = dt;
  fpu->sr = sr_default;
  fpu_clear_state(fpu);
}

void fpu_clear_state(struct fpu* fpu) {
  fpu->phase = 0.0;
  for(int i=0; i<FPU_NUM_MASSES; i++) {
	fpu->x[i] = 0.0;
	fpu->v[i] = 0.0;
	fpu->f[i] = 0.0;
  }
}

void fpu_set_sr(struct fpu* fpu, float val) {
  fpu->sr = val;
  // FIXME: might want to adjust dt and beta when SR changes...?
}

void fpu_set_pos(struct fpu* fpu, int pos, double val) {
  if(pos > 1 && pos < (FPU_NUM_MASSES-2)) {
	if(val > 1.0) { val = 1.0; }
	if(val < -1.0) { val = -1.0; }
	fpu->x[pos] = val;
  }
}

void fpu_set_beta(struct fpu* fpu, double val) {
  fpu->beta = val;
}

void fpu_set_rho(struct fpu* fpu, double val) {
  fpu->rho = val;
}

void fpu_set_epsilon(struct fpu* fpu, double val) {
  fpu->epsilon = val;
}

void fpu_update(struct fpu* fpu, double input) {
  for(int i=0; i<ips; i++) { 
	fpu_update_f(fpu);
	fpu_update_v(fpu);
	fpu_update_x(fpu);
  }
}

double fpu_get_output(struct fpu* fpu, double pos) {
  int pi0 =  (int)pos;
  double pf = pos - (double)pi0;
  pi0 = pi0 % (FPU_NUM_MASSES);
  int pi1 = (pi0 + 1) % FPU_NUM_MASSES;
  return fpu->x[pi0] + pf*(fpu->x[pi1] - fpu->x[pi0]);
}

fpu.h

#pragma once
/*
  fermi-pasta-ulam resonator

  this is a "physical" model of a nonlinear string.

  the string is a number of connected masses

  the force acting on a given mass is a function of its displacement.
  this function has a linear term and a variable nonlinear (cubic) term

  d1 = upper displacement
  d0 = lower displacement

  f = b* [ (d1 - d0) + e*(d1^3 - d0^3) ]

*/


#define FPU_NUM_MASSES 8

struct fpu {
  //--- synthesis parameters
  // sample rate
  float sr;
  // iterations per sample (function of samplerate)
  double ips;
  // phasor / counter to decouple time step from sample rate
  double phase;
  // output position (normalized)
  double pickup;
  // output value
  double out;
  //----- model parameters
  // "physical" time step (in seconds)
  double dt;
  // tension parameter (controls pitch)
  double beta;
  // restoring force (controls duration)
  double rho;
  // nonlinear parameter
  double epsilon;  
  // array of forces
  double f[FPU_NUM_MASSES];
  // array of velocities
  double v[FPU_NUM_MASSES];
  // array of positions
  double x[FPU_NUM_MASSES];
  
  
};

// initialize
extern void fpu_init(struct fpu* fpu);
// clear the state of the model
extern void fpu_clear_state(struct fpu* fpu);
// set sample rate
extern void fpu_set_sr(struct fpu* fpu, float val);
// set the position of a given mass
extern void fpu_set_pos(struct fpu* fpu, int pos, double val);
// set the beta coefficient (restoring force / tuning)
extern void fpu_set_beta(struct fpu* fpu, double val);
// set the rho coefficient (damping / decay)
extern void fpu_set_rho(struct fpu* fpu, double val);
// set epsilon coefficient (nonlinearity)
extern void fpu_set_epsilon(struct fpu* fpu, double val);
// set the time step of the model
extern void fpu_set_dt(struct fpu* fpu, double val);
// set the pickup position (applied as offset to all outputs)
// extern void fpu_set_pickup(struct fpu* fpu, double val);
// compute a single time step
extern void fpu_update(struct fpu* fpu, double input);
// get output value for given pickup position (interpolated)
extern double fpu_get_output(struct fpu* fpu, double pos);

fpu_pd.c

#include "../m_pd.h"
#include "fpu.h"

static t_class *fpu_tilde_class;

typedef struct _fpu_tilde {
  t_object x_obj;
  t_outlet *x_out[FPU_NUM_MASSES];
  t_sample *x_outvec[FPU_NUM_MASSES];
  float f;
  struct fpu fpu;
  double pickup_pos;
} t_fpu_tilde;


void* fpu_tilde_new(void) {
  t_fpu_tilde *x = (t_fpu_tilde *)pd_new(fpu_tilde_class);
  for (int i=0; i<FPU_NUM_MASSES; i++) {
	x->x_out[i] = outlet_new(&x->x_obj, &s_signal);
  }
  fpu_init(&x->fpu);
  return (void*)x;
}

void fpu_tilde_free(t_fpu_tilde *x) {
  for (int i=0; i<FPU_NUM_MASSES; i++) {
	outlet_free(x->x_out[i]);
  }
}

t_int* fpu_tilde_perform(t_int *w) {
  t_fpu_tilde* x = (t_fpu_tilde*)(w[1]);
  t_sample *in = (t_sample*)(w[2]);
  t_sample *out[FPU_NUM_MASSES];
  int n = (int)(w[3]);
  
  for(int i=0; i<FPU_NUM_MASSES; i++) {
	out[i] = x->x_outvec[i];
  }

  while(n--) {
	fpu_update(&x->fpu, *in++);
	for(int i=FPU_NUM_MASSES - 1; i >= 0; i--) {
	  *(out[i]) = fpu_get_output( &x->fpu, x->pickup_pos + i);
	  (out[i])++;
	}
  }
  return (w + 4);
}

void fpu_tilde_dsp(t_fpu_tilde *x, t_signal **sp) {
  
  fpu_set_sr(&x->fpu, sp[0]->s_sr);
  for(int i=0; i<FPU_NUM_MASSES; i++) {
	x->x_outvec[i] = sp[i+1]->s_vec;
  }
  
  dsp_add(fpu_tilde_perform, 3, x,
		  // input
		  sp[0]->s_vec,
		  sp[0]->s_n );
}

void fpu_tilde_beta(t_fpu_tilde *x, t_floatarg val) {
  fpu_set_beta(&x->fpu, val);
}

void fpu_tilde_pos(t_fpu_tilde *x, t_floatarg pos, t_floatarg val) {
  fpu_set_pos(&x->fpu, (int)pos, val);
}

void fpu_tilde_rho(t_fpu_tilde *x, t_floatarg val) {
  fpu_set_rho(&x->fpu, val);
}

void fpu_tilde_epsilon(t_fpu_tilde *x, t_floatarg val) {
  fpu_set_epsilon(&x->fpu, val);
}

void fpu_tilde_pickup(t_fpu_tilde *x, t_floatarg val) {
  x->pickup_pos = (double)val;
}

void fpu_tilde_clear(t_fpu_tilde *x, t_floatarg val) {
  fpu_clear_state(&x->fpu);
}
void fpu_tilde_setup(void) {
  fpu_tilde_class = class_new(gensym("fpu~"),
							  (t_newmethod)fpu_tilde_new,
							  0, sizeof(t_fpu_tilde),
							  CLASS_DEFAULT, 0);
  
  class_addmethod(fpu_tilde_class,  
				  (t_method)fpu_tilde_dsp, gensym("dsp"), 0);
    
  class_addmethod(fpu_tilde_class,
				  (t_method)fpu_tilde_pos, gensym("pos"), A_DEFFLOAT, A_DEFFLOAT, 0 );
  class_addmethod(fpu_tilde_class,
				  (t_method)fpu_tilde_beta, gensym("beta"), A_DEFFLOAT, 0 );
  
  class_addmethod(fpu_tilde_class,
				  (t_method)fpu_tilde_rho, gensym("rho"), A_DEFFLOAT, 0 );

  class_addmethod(fpu_tilde_class,
				  (t_method)fpu_tilde_epsilon, gensym("epsilon"), A_DEFFLOAT, 0 );
  
  class_addmethod(fpu_tilde_class,
				  (t_method)fpu_tilde_pickup, gensym("pickup"), A_DEFFLOAT, 0 );

  class_addmethod(fpu_tilde_class,
				  (t_method)fpu_tilde_clear, gensym("clear"), 0);

  CLASS_MAINSIGNALIN(fpu_tilde_class, t_fpu_tilde, f);  
}

Makefile

NAME = fpu~

CFLAGS = -fPIC -DLINUX

LDFLAGS = -export-dynamic -shared\
	-nostartfiles -nodefaultlibs -nostdlib

LIBS = -shared -Wl,--export-dynamic -lm

INC += ../

SRC = fpu_pd.c \
	fpu.c

OBJ = $(patsubst %.c, %.o, $(SRC))

TARGET = $(NAME).pd_linux

#STRIP = strip --strip-unneeded $(TARGET)

$(TARGET) : $(OBJ)
	gcc $(LDFLAGS) -o $(TARGET) $(OBJ) $(LIBS)
	$(STRIP)

clean:
	rm $(TARGET)
	rm *.o

install: $(TARGET)
	cp $(TARGET) ~/pd-externals