ref: a90cbc5b8eac9a63d3e0bbe8f326af7adf44b395
parent: 635257887ab374a68acf386814df0c59aeec41a1
author: spiricom <jeff@snyderphonics.com>
date: Mon Apr 27 17:17:28 EDT 2020
wavefolder edits
binary files a/.DS_Store b/.DS_Store differ
--- a/LEAF/Inc/leaf-distortion.h
+++ b/LEAF/Inc/leaf-distortion.h
@@ -80,7 +80,7 @@
{
double Ln1;
double Fn1;
- float xn1;
+ double xn1;
double RL;
double R;
@@ -92,9 +92,16 @@
double d;
// Antialiasing error threshold
- double thresh;
+ double AAthresh;
double half_a;
double longthing;
+
+ double LambertThresh, w, expw, p, r, s, myerr;
+ double l;
+ double u, Ln, Fn;
+ double tempsDenom;
+ double tempErrDenom;
+ double tempOutDenom;
} _tLockhartWavefolder;
--- a/LEAF/Inc/leaf-filters.h
+++ b/LEAF/Inc/leaf-filters.h
@@ -417,9 +417,13 @@
float fs; // sample-rate
float fc; // characteristic frequency
float G; // gain
+ float invG; //1/gain
float B; // bandwidth (in octaves)
float m; // morph parameter (0...1)
+ float sr; //local sampling rate of filter (may be different from leaf sr if oversampled)
+ float inv_sr;
+
} _tVZFilter;
typedef _tVZFilter* tVZFilter;
@@ -428,7 +432,7 @@
void tVZFilter_free (tVZFilter* const);
void tVZFilter_initToPool (tVZFilter* const, VZFilterType type, float freq, float Q, tMempool* const);
void tVZFilter_freeFromPool (tVZFilter* const, tMempool* const);
-
+ void tVZFilter_setSampleRate (tVZFilter* const, float sampleRate);
float tVZFilter_tick (tVZFilter* const, float input);
void tVZFilter_calcCoeffs (tVZFilter* const);
void tVZFilter_setBandwidth (tVZFilter* const, float bandWidth);
@@ -436,6 +440,37 @@
void tVZFilter_setGain (tVZFilter* const, float gain);
void tVZFilter_setType (tVZFilter* const, VZFilterType type);
float tVZFilter_BandwidthToR (tVZFilter* const vf, float B);
+
+
+
+
+ //diode ladder filter
+ typedef struct _tDiodeFilter
+ {
+
+ float f;
+ float r;
+ float Vt;
+ float n;
+ float gamma;
+ float zi;
+ float g0inv;
+ float g1inv;
+ float g2inv;
+ float s0, s1, s2, s3;
+
+ } _tDiodeFilter;
+
+ typedef _tDiodeFilter* tDiodeFilter;
+
+ void tDiodeFilter_init (tDiodeFilter* const, float freq, float Q);
+ void tDiodeFilter_free (tDiodeFilter* const);
+ void tDiodeFilter_initToPool (tDiodeFilter* const, float freq, float Q, tMempool* const);
+ void tDiodeFilter_freeFromPool (tDiodeFilter* const, tMempool* const);
+
+ float tDiodeFilter_tick (tDiodeFilter* const, float input);
+ void tDiodeFilter_setFreq (tDiodeFilter* const vf, float cutoff);
+ void tDiodeFilter_setQ (tDiodeFilter* const vf, float resonance);
#ifdef __cplusplus
--- a/LEAF/Inc/leaf-math.h
+++ b/LEAF/Inc/leaf-math.h
@@ -145,10 +145,21 @@
void LEAF_crossfade(float fade, float* volumes);
+ float LEAF_tanh(float x);
+ void LEAF_generate_sine(float* buffer, int size);
+ void LEAF_generate_sawtooth(float* buffer, float basefreq, int size);
+ void LEAF_generate_triangle(float* buffer, float basefreq, int size);
+ void LEAF_generate_square(float* buffer, float basefreq, int size);
+ float LEAF_midiToFrequency(float f);
+
float fast_mtof(float f);
+ double fastexp(double x);
+
float fastexpf(float x);
+
+ double fasterexp(double x);
float fasterexpf(float x);
--- a/LEAF/Src/leaf-distortion.c
+++ b/LEAF/Src/leaf-distortion.c
@@ -17,7 +17,7 @@
#include "../Inc/leaf-tables.h"
//testing
-//#include "gpio.h"
+#include "gpio.h"
#endif
@@ -353,26 +353,7 @@
void tLockhartWavefolder_init(tLockhartWavefolder* const wf)
{- _tLockhartWavefolder* w = *wf = (_tLockhartWavefolder*) leaf_alloc(sizeof(_tLockhartWavefolder));
-
- w->Ln1 = 0.0;
- w->Fn1 = 0.0;
- w->xn1 = 0.0f;
-
- w->RL = 7.5e3;
- w->R = 15e3;
- w->VT = 26e-3;
- w->Is = 10e-16;
-
- w->a = 2.0*w->RL/w->R;
- w->b = (w->R+2.0*w->RL)/(w->VT*w->R);
- w->d = (w->RL*w->Is)/w->VT;
- w->half_a = 0.5 * w->a;
- w->longthing = (0.5*w->VT/w->b);
-
-
- // Antialiasing error threshold
- w->thresh = 10e-10;
+ tLockhartWavefolder_initToPool (wf, &leaf.mempool);
}
void tLockhartWavefolder_free(tLockhartWavefolder* const wf)
@@ -389,7 +370,7 @@
w->Ln1 = 0.0;
w->Fn1 = 0.0;
- w->xn1 = 0.0f;
+ w->xn1 = 0.0;
w->RL = 7.5e3;
w->R = 15e3;
@@ -404,7 +385,26 @@
// Antialiasing error threshold
- w->thresh = 10e-10;
+ w->AAthresh = 10e-10; //10
+
+ w->LambertThresh = 10e-10; //12 //was 8
+
+
+ w->w = 0.0f;
+ w->expw = 0.0f;
+ w->p = 0.0f;
+ w->r = 0.0f;
+ w->s= 0.0f;
+ w->myerr = 0.0f;
+ w->l = 0.0f;
+ w->u = 0.0f;
+ w->Ln = 0.0f;
+ w->Fn = 0.0f;
+ w->tempsDenom = 0.0f;
+ w->tempErrDenom = 0.0f;
+ w->tempOutDenom = 0.0f;
+
+
}
void tLockhartWavefolder_freeFromPool (tLockhartWavefolder* const wf, tMempool* const mp)
@@ -415,71 +415,182 @@
mpool_free(w, m);
}
-double tLockhartWavefolderLambert(double x, double ln)
+
+
+double tLockhartWavefolderLambert(tLockhartWavefolder* const wf, double x, double ln)
{- double thresh, w, expw, p, r, s, err;
- // Error threshold
- thresh = 10e-12;
+ _tLockhartWavefolder* mwf = *wf;
+
+
// Initial guess (use previous value)
- w = ln;
+ mwf->w = ln;
// Haley's method (Sec. 4.2 of the paper)
- for(int i=0; i<1000; i+=1) {+ for(int i=0; i<3000; i+=1) { //1000- expw = exp(w);
-
- p = w*expw - x;
- r = (w+1.0)*expw;
- s = (w+2.0)/(2.0*(w+1.0));
- err = (p/(r-(p*s)));
-
- if (fabs(err)<thresh) {+ mwf->expw = exp(mwf->w);
+ /*
+ if (isinf(mwf->expw) || isnan(mwf->expw))
+ {+ mwf->expw = 10e-5;
+ LEAF_error();
+ }
+ */
+ mwf->p = mwf->w*mwf->expw - x;
+ /*
+ if (isinf(mwf->p) || isnan(mwf->p))
+ {+ mwf->p = 10e-5;
+ LEAF_error();
+ }
+ */
+ mwf->r = (mwf->w+1.0)*mwf->expw;
+ /*
+ if (isinf(mwf->r) || isnan(mwf->r))
+ {+ mwf->r = 10e-5;
+ LEAF_error();
+ }
+ */
+ mwf->tempsDenom = (2.0*(mwf->w+1.0));
+ /*
+ if ((mwf->tempsDenom == 0.0) || isinf(mwf->tempsDenom) || isnan(mwf->tempsDenom))
+ {+ mwf->tempsDenom = 10e-5;
+ LEAF_error();
+ }
+ */
+ mwf->s = (mwf->w+2.0)/mwf->tempsDenom;
+ /*
+ if (isnan(mwf->s) || isinf(mwf->s))
+ {+ mwf->s = 10e-5;
+ LEAF_error();
+ }
+ */
+ mwf->tempErrDenom = (mwf->r-(mwf->p*mwf->s));
+ /*
+ if ((mwf->tempErrDenom == 0.0) || isinf(mwf->tempErrDenom) || isnan(mwf->tempErrDenom))
+ {+ mwf->tempErrDenom = 10e-5;
+ LEAF_error();
+ }
+ */
+ mwf->myerr = (mwf->p/mwf->tempErrDenom);
+ /*
+ if (isnan(mwf->myerr) || isinf(mwf->myerr))
+ {+ mwf->myerr = 10e-5;
+ LEAF_error();
+ }
+ */
+ if ((fabs(mwf->myerr))<mwf->LambertThresh) {+
break;
}
-
- w = w - err;
- if (i == 999)
+
+ mwf->w = mwf->w - mwf->myerr;
+
+ /*
+ if (isinf(mwf->w) || isnan(mwf->w))
{- //HAL_GPIO_WritePin(GPIOG, GPIO_PIN_7, GPIO_PIN_SET);
+ mwf->w = 10e-5;
+ LEAF_error();
}
+*/
}
- return w;
+ return mwf->w;
}
-float tLockhartWavefolder_tick(tLockhartWavefolder* const wf, float samp)
+float tLockhartWavefolder_tick(tLockhartWavefolder* const wf, float in)
{_tLockhartWavefolder* w = *wf;
-
+
float out = 0.0f;
// Compute Antiderivative
- int l = (samp > 0.0f) - (samp < 0.0f);
- double u = w->d*exp(l*w->b*samp);
- double Ln = tLockhartWavefolderLambert(u,w->Ln1);
- double Fn = w->longthing*(Ln*(Ln + 2.0)) - w->half_a*samp*samp;
-
+ w->l = (in > 0.0) - (in < 0.0);
+ w->u = w->d*exp(w->l*w->b*in);
+ /*
+ if ((w->u == 0.0) || isinf(w->u) || isnan(w->u))
+ {+ w->u = 10e-5;
+ LEAF_error();
+ }
+ */
+
+ w->Ln = tLockhartWavefolderLambert(wf,w->u,w->Ln1);
+ /*
+ if ((w->Ln == 0.0) || isinf(w->Ln) || isnan(w->Ln))
+ {+ w->Ln = 10e-5;
+ LEAF_error();
+ }
+*/
+ w->Fn = (w->longthing*(w->Ln*(w->Ln + 2.0))) - (w->half_a*in*in);
+ /*
+ if ((w->Fn == 0.0) || isinf(w->Fn) || isnan(w->Fn))
+ {+ w->Fn = 10e-5;
+ LEAF_error();
+ }
+ */
// Check for ill-conditioning
- if (fabs(samp-w->xn1)<w->thresh) {+
+ if (fabs(in-w->xn1)<w->AAthresh)
+ {// Compute Averaged Wavefolder Output
- float xn = 0.5f*(samp+w->xn1);
- u = w->d*exp(l*w->b*xn);
- Ln = tLockhartWavefolderLambert(u,w->Ln1);
- out = (float) (l*w->VT*Ln - w->a*xn);
+ double xn = 0.5*(in+w->xn1);
+ w->u = w->d*exp(w->l*w->b*xn);
+ /*
+ if ((w->u == 0.0) || isinf(w->u) || isnan(w->u))
+ {+ w->u = 10e-5;
+ LEAF_error();
+ }
+ */
+ w->Ln = tLockhartWavefolderLambert(wf,w->u,w->Ln1);
+ /*
+ if ((w->Ln == 0.0) || isinf(w->Ln) || isnan(w->Ln))
+ {+ w->Ln = 10e-5;
+ LEAF_error();
+ }
+ */
+ out = (float)((w->l*w->VT*w->Ln) - (w->a*xn));
+
}
- else {-
+ else
+ {+
// Apply AA Form
- out = (float) ((Fn-w->Fn1)/(samp-w->xn1));
+ w->tempOutDenom = (in-w->xn1);
+ /*
+ if ((w->tempOutDenom == 0.0) || isinf(w->tempOutDenom))
+ {+ w->tempOutDenom = 10e-5;
+ LEAF_error();
+ }
+ */
+ out = ((w->Fn-w->Fn1)/w->tempOutDenom);
+ /*
+ if (isinf(out) || isnan(out))
+ {+ out = 10e-5;
+ LEAF_error();
+ }
+ */
+
}
-
+
// Update States
- w->Ln1 = Ln;
- w->Fn1 = Fn;
- w->xn1 = samp;
+ w->Ln1 = w->Ln;
+ w->Fn1 = w->Fn;
+ w->xn1 = (double)in;
return out;
}
--- a/LEAF/Src/leaf-filters.c
+++ b/LEAF/Src/leaf-filters.c
@@ -17,7 +17,7 @@
#include "../Inc/leaf-filters.h"
#include "../Inc/leaf-tables.h"
#include "../leaf.h"
-
+#include "tim.h"
#endif
// ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ OnePole Filter ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ //
@@ -1398,6 +1398,8 @@
f->m = 0.0f;
f->s1 = 0.0f;
f->s2 = 0.0f;
+ f->sr = leaf.sampleRate;
+ f->inv_sr = leaf.invSampleRate;
tVZFilter_calcCoeffs(vf);
@@ -1409,6 +1411,13 @@
mpool_free(f, m);
}
+void tVZFilter_setSampleRate (tVZFilter* const vf, float sampleRate)
+{
+ _tVZFilter* f = *vf;
+ f->sr = sampleRate;
+ f->inv_sr = 1.0f/sampleRate;
+}
+
float tVZFilter_tick (tVZFilter* const vf, float in)
{
_tVZFilter* f = *vf;
@@ -1444,13 +1453,13 @@
{
_tVZFilter* f = *vf;
- f->g = tanf(PI * f->fc * leaf.invSampleRate); // embedded integrator gain (Fig 3.11)
+ f->g = tanf(PI * f->fc * f->inv_sr); // embedded integrator gain (Fig 3.11)
switch( f->type )
{
case Bypass:
{
- f->R2 = 1.0f / f->G; // can we use an arbitrary value here, for example R2 = 1?
+ f->R2 = f->invG; // can we use an arbitrary value here, for example R2 = 1?
f->cL = 1.0f;
f->cB = f->R2;
f->cH = 1.0f;
@@ -1458,19 +1467,19 @@
break;
case Lowpass:
{
- f->R2 = 1.0f / f->G;
+ f->R2 = f->invG;
f->cL = 1.0f; f->cB = 0.0f; f->cH = 0.0f;
}
break;
case Highpass:
{
- f->R2 = 1.0f / f->G;
+ f->R2 = f->invG;
f->cL = 0.0f; f->cB = 0.0f; f->cH = 1.0f;
}
break;
case BandpassSkirt:
{
- f->R2 = 1.0f / f->G;
+ f->R2 = f->invG;
f->cL = 0.0f; f->cB = 1.0f; f->cH = 0.0f;
}
break;
@@ -1489,7 +1498,7 @@
case Bell:
{
float fl = f->fc*powf(2.0f, (-f->B)*0.5f); // lower bandedge frequency (in Hz)
- float wl = tanf(PI*fl*leaf.invSampleRate); // warped radian lower bandedge frequency /(2*fs)
+ float wl = tanf(PI*fl*f->inv_sr); // warped radian lower bandedge frequency /(2*fs)
float r = f->g/wl;
r *= r; // warped frequency ratio wu/wl == (wc/wl)^2 where wu is the
// warped upper bandedge, wc the center
@@ -1523,7 +1532,7 @@
// experimental - maybe we must find better curves for cL, cB, cH:
case Morph:
{
- f->R2 = 1.0f / f->G;
+ f->R2 = f->invG;
float x = 2.0f*f->m-1.0f;
f->cL = maximum(-x, 0.0f); /*cL *= cL;*/
@@ -1562,6 +1571,7 @@
{
_tVZFilter* f = *vf;
f->G = LEAF_clip(0.000001f, gain, 100.0f);
+ f->invG = 1.0f/f->G;
tVZFilter_calcCoeffs(vf);
}
@@ -1583,9 +1593,148 @@
{
_tVZFilter* f = *vf;
float fl = f->fc*powf(2.0f, -B*0.5f); // lower bandedge frequency (in Hz)
- float gl = tanf(PI*fl*leaf.invSampleRate); // warped radian lower bandedge frequency /(2*fs)
+ float gl = tanf(PI*fl*f->inv_sr); // warped radian lower bandedge frequency /(2*fs)
float r = gl/f->g; // ratio between warped lower bandedge- and center-frequencies
// unwarped: r = pow(2, -B/2) -> approximation for low
// center-frequencies
return sqrtf((1.0f-r*r)*(1.0f-r*r)/(4.0f*r*r));
}
+
+
+
+void tDiodeFilter_init (tDiodeFilter* const vf, float cutoff, float resonance)
+{
+ tDiodeFilter_initToPool(vf, cutoff, resonance, &leaf.mempool);
+}
+
+void tDiodeFilter_free (tDiodeFilter* const vf)
+{
+ tDiodeFilter_freeFromPool(vf, &leaf.mempool);
+}
+void tDiodeFilter_initToPool (tDiodeFilter* const vf, float cutoff, float resonance, tMempool* const mp)
+{
+
+ _tMempool* m = *mp;
+ _tDiodeFilter* f = *vf = (_tDiodeFilter*) mpool_alloc(sizeof(_tDiodeFilter), m);
+ // initialization (the resonance factor is between 0 and 8 according to the article)
+ f->f = tan(PI * cutoff/leaf.sampleRate);
+ f->r = (7.f * resonance + 0.5f);
+ f->Vt = 0.5f;
+ f->n = 1.836f;
+ f->zi = 0.0f; //previous input value
+ f->gamma = f->Vt*f->n;
+ f->s0 = 0.01f;
+ f->s1 = 0.02f;
+ f->s2 = 0.03f;
+ f->s3 = 0.04f;
+ f->g0inv = 1.f/(2.f*f->Vt);
+ f->g1inv = 1.f/(2.f*f->gamma);
+ f->g2inv = 1.f/(6.f*f->gamma);
+
+
+}
+void tDiodeFilter_freeFromPool (tDiodeFilter* const vf, tMempool* const mp)
+{
+ _tMempool* m = *mp;
+ _tDiodeFilter* f = *vf = (_tDiodeFilter*) mpool_alloc(sizeof(_tDiodeFilter), m);
+ mpool_free(f, m);
+}
+
+float tanhXdX(float x)
+{
+ float a = x*x;
+ // IIRC I got this as Pade-approx for tanh(sqrt(x))/sqrt(x)
+ return ((a + 105.0f)*a + 945.0f) / ((15.0f*a + 420.0f)*a + 945.0f);
+}
+
+float tDiodeFilter_tick (tDiodeFilter* const vf, float in)
+{
+ _tDiodeFilter* f = *vf;
+
+ // the input x[n+1] is given by 'in', and x[n] by zi
+ // input with half delay
+ float ih = 0.5f * (in + f->zi);
+
+ // evaluate the non-linear factors
+ float t0 = f->f*tanhXdX((ih - f->r * f->s3)*f->g0inv)*f->g0inv;
+ float t1 = f->f*tanhXdX((f->s1-f->s0)*f->g1inv)*f->g1inv;
+ float t2 = f->f*tanhXdX((f->s2-f->s1)*f->g1inv)*f->g1inv;
+ float t3 = f->f*tanhXdX((f->s3-f->s2)*f->g1inv)*f->g1inv;
+ float t4 = f->f*tanhXdX((f->s3)*f->g2inv)*f->g2inv;
+
+
+
+ // This formula gives the result for y3 thanks to MATLAB
+ float y3 = (f->s2 + f->s3 + t2*(f->s1 + f->s2 + f->s3 + t1*(f->s0 + f->s1 + f->s2 + f->s3 + t0*in)) + t1*(2.0f*f->s2 + 2.0f*f->s3))*t3 + f->s3 + 2.0f*f->s3*t1 + t2*(2.0f*f->s3 + 3.0f*f->s3*t1);
+ if (isnan(y3))
+ {
+ __HAL_TIM_SET_COMPARE(&htim3, TIM_CHANNEL_2, 400);
+ }
+ float tempy3denom = (t4 + t1*(2.0f*t4 + 4.0f) + t2*(t4 + t1*(t4 + f->r*t0 + 4.0f) + 3.0f) + 2.0f)*t3 + t4 + t1*(2.0f*t4 + 2.0f) + t2*(2.0f*t4 + t1*(3.0f*t4 + 3.0f) + 2.0f) + 1.0f;
+ if (isnan(tempy3denom))
+ {
+ __HAL_TIM_SET_COMPARE(&htim3, TIM_CHANNEL_2, 400);
+ }
+ if (tempy3denom == 0.0f)
+ {
+ tempy3denom = 0.000001f;
+ }
+ y3 = y3 / tempy3denom;
+ if (isnan(y3))
+ {
+ __HAL_TIM_SET_COMPARE(&htim3, TIM_CHANNEL_2, 400);
+ }
+ if (t1 == 0.0f)
+ {
+ t1 = 0.000001f;
+ }
+ if (t2 == 0.0f)
+ {
+ t2 = 0.000001f;
+ }
+ if (t3 == 0.0f)
+ {
+ t3 = 0.000001f;
+ }
+ // Other outputs
+ float y2 = (f->s3 - (1+t4+t3)*y3) / (-t3);
+ float y1 = (f->s2 - (1+t3+t2)*y2 + t3*y3) / (-t2);
+ float y0 = (f->s1 - (1+t2+t1)*y1 + t2*y2) / (-t1);
+ float xx = (in - f->r*y3);
+
+ // update state
+ f->s0 += 2.0f * (t0*xx + t1*(y1-y0));
+ if (isnan(f->s0))
+ {
+ __HAL_TIM_SET_COMPARE(&htim3, TIM_CHANNEL_2, 400);
+ }
+
+ if (isinf(f->s0))
+ {
+ __HAL_TIM_SET_COMPARE(&htim3, TIM_CHANNEL_2, 400);
+ }
+ f->s1 += 2.0f * (t2*(y2-y1) - t1*(y1-y0));
+ f->s2 += 2.0f * (t3*(y3-y2) - t2*(y2-y1));
+ f->s3 += 2.0f * (-t4*(y3) - t3*(y3-y2));
+
+ f->zi = in;
+ return y3*f->r;
+
+}
+
+
+
+void tDiodeFilter_setFreq (tDiodeFilter* const vf, float cutoff)
+{
+ _tDiodeFilter* f = *vf;
+ f->f = tanf(PI * LEAF_clip(10.0f, cutoff, 20000.0f)*leaf.invSampleRate);
+}
+
+
+void tDiodeFilter_setQ (tDiodeFilter* const vf, float resonance)
+{
+ _tDiodeFilter* f = *vf;
+ f->r = LEAF_clip(0.5, (7.f * resonance + 0.5f), 8.0f);
+
+}
+
--- a/LEAF/Src/leaf-math.c
+++ b/LEAF/Src/leaf-math.c
@@ -90,6 +90,13 @@
return alias.f;
}
+double fastexp(double x) {
+ x = 1.0 + (x * 0.0009765625);
+ x *= x; x *= x; x *= x; x *= x;
+ x *= x; x *= x; x *= x; x *= x;
+ x *= x; x *= x;
+ return x;
+}
float fastexpf(float x) {
x = 1.0f + (x * 0.0009765625f);
@@ -99,8 +106,15 @@
return x;
}
+double fasterexp(double x) {
+ x = 1.0 + (x * 0.00390625);
+ x *= x; x *= x; x *= x; x *= x;
+ x *= x; x *= x; x *= x; x *= x;
+ return x;
+}
+
float fasterexpf(float x) {
- x = 1.0 + (x * 0.00390625f);
+ x = 1.0f + (x * 0.00390625f);
x *= x; x *= x; x *= x; x *= x;
x *= x; x *= x; x *= x; x *= x;
return x;