shithub: leaf

ref: 4cb315f58a570ac0f640c6ea135b70fcbbf1f32f
dir: /LEAF/Src/leaf-filters.c/

View raw version
/*==============================================================================

    leaf-filter.c
    Created: 20 Jan 2017 12:01:10pm
    Author:  Michael R Mulshine

==============================================================================*/

#if _WIN32 || _WIN64

#include "..\Inc\leaf-filters.h"
#include "..\Inc\leaf-tables.h"
#include "..\leaf.h"

#else

#include "../Inc/leaf-filters.h"
#include "../Inc/leaf-tables.h"
#include "../leaf.h"
//#include "tim.h"
#endif

// ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ OnePole Filter ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ //
void    tAllpass_init(tAllpass* const ft, float initDelay, uint32_t maxDelay)
{
    _tAllpass* f = *ft = (_tAllpass*) leaf_alloc(sizeof(_tAllpass));
    
    f->gain = 0.7f;
    
    f->lastOut = 0.0f;
    
    tLinearDelay_init(&f->delay, initDelay, maxDelay);
}

void tAllpass_free(tAllpass* const ft)
{
    _tAllpass* f = *ft;
    
    tLinearDelay_free(&f->delay);
    leaf_free(f);
}

void    tAllpass_initToPool     (tAllpass* const ft, float initDelay, uint32_t maxDelay, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tAllpass* f = *ft = (_tAllpass*) mpool_alloc(sizeof(_tAllpass), m);
    
    f->gain = 0.7f;
    
    f->lastOut = 0.0f;
    
    tLinearDelay_initToPool(&f->delay, initDelay, maxDelay, mp);
}

void    tAllpass_freeFromPool   (tAllpass* const ft, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tAllpass* f = *ft;
    
    tLinearDelay_freeFromPool(&f->delay, mp);
    mpool_free(f, m);
}

void    tAllpass_setDelay(tAllpass* const ft, float delay)
{
    _tAllpass* f = *ft;
    
    tLinearDelay_setDelay(&f->delay, delay);
}

void    tAllpass_setGain(tAllpass* const ft, float gain)
{
    _tAllpass* f = *ft;
    
    f->gain = gain;
}

float   tAllpass_tick(tAllpass* const ft, float input)
{
    _tAllpass* f = *ft;
    
    float s1 = (-f->gain) * f->lastOut + input;
    
    float s2 = tLinearDelay_tick(&f->delay, s1) + (f->gain) * input;
    
    f->lastOut = s2;
    
    return f->lastOut;
}

// ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ OnePole Filter ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ //
void    tOnePole_init(tOnePole* const ft, float freq)
{
    _tOnePole* f = *ft = (_tOnePole*) leaf_alloc(sizeof(_tOnePole));
    
    f->gain = 1.0f;
    f->a0 = 1.0;
    
    tOnePole_setFreq(ft, freq);
    
    f->lastIn = 0.0f;
    f->lastOut = 0.0f;
}

void    tOnePole_free(tOnePole* const ft)
{
    _tOnePole* f = *ft;
    
    leaf_free(f);
}

void    tOnePole_initToPool     (tOnePole* const ft, float freq, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tOnePole* f = *ft = (_tOnePole*) mpool_alloc(sizeof(_tOnePole), m);
    
    f->gain = 1.0f;
    f->a0 = 1.0;
    
    tOnePole_setFreq(ft, freq);
    
    f->lastIn = 0.0f;
    f->lastOut = 0.0f;
}

void    tOnePole_freeFromPool   (tOnePole* const ft, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tOnePole* f = *ft;
    
    mpool_free(f, m);
}

void    tOnePole_setB0(tOnePole* const ft, float b0)
{
    _tOnePole* f = *ft;
    f->b0 = b0;
}

void    tOnePole_setA1(tOnePole* const ft, float a1)
{
    _tOnePole* f = *ft;
    if (a1 >= 1.0f)     a1 = 0.999999f;
    f->a1 = a1;
}

void    tOnePole_setPole(tOnePole* const ft, float thePole)
{
    _tOnePole* f = *ft;
    
    if (thePole >= 1.0f)    thePole = 0.999999f;
    
    // Normalize coefficients for peak unity gain.
    if (thePole > 0.0f)     f->b0 = (1.0f - thePole);
    else                    f->b0 = (1.0f + thePole);
    
    f->a1 = -thePole;
}

void        tOnePole_setFreq        (tOnePole* const ft, float freq)
{
    _tOnePole* f = *ft;
    f->b0 = freq * leaf.twoPiTimesInvSampleRate;
    f->b0 = LEAF_clip(0.0f, f->b0, 1.0f);
    f->a1 = 1.0f - f->b0;
}

void    tOnePole_setCoefficients(tOnePole* const ft, float b0, float a1)
{
    _tOnePole* f = *ft;
    if (a1 >= 1.0f)     a1 = 0.999999f;
    f->b0 = b0;
    f->a1 = a1;
}

void    tOnePole_setGain(tOnePole* const ft, float gain)
{
    _tOnePole* f = *ft;
    f->gain = gain;
}

float   tOnePole_tick(tOnePole* const ft, float input)
{
    _tOnePole* f = *ft;
    
    float in = input * f->gain;
    float out = (f->b0 * in) + (f->a1 * f->lastOut);
    
    f->lastIn = in;
    f->lastOut = out;
    
    return out;
}

// ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ TwoPole Filter ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ //
void    tTwoPole_init(tTwoPole* const ft)
{
    _tTwoPole* f = *ft = (_tTwoPole*) leaf_alloc(sizeof(_tTwoPole));
    
    f->gain = 1.0f;
    f->a0 = 1.0;
    f->b0 = 1.0;
    
    f->lastOut[0] = 0.0f;
    f->lastOut[1] = 0.0f;
}

void    tTwoPole_free(tTwoPole* const ft)
{
    _tTwoPole* f = *ft;
    
    leaf_free(f);
}

void    tTwoPole_initToPool     (tTwoPole* const ft, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tTwoPole* f = *ft = (_tTwoPole*) mpool_alloc(sizeof(_tTwoPole), m);
    
    f->gain = 1.0f;
    f->a0 = 1.0;
    f->b0 = 1.0;
    
    f->lastOut[0] = 0.0f;
    f->lastOut[1] = 0.0f;
}

void    tTwoPole_freeFromPool   (tTwoPole* const ft, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tTwoPole* f = *ft;
    
    mpool_free(f, m);
}

float   tTwoPole_tick(tTwoPole* const ft, float input)
{
    _tTwoPole* f = *ft;
    
    float in = input * f->gain;
    float out = (f->b0 * in) - (f->a1 * f->lastOut[0]) - (f->a2 * f->lastOut[1]);
    
    f->lastOut[1] = f->lastOut[0];
    f->lastOut[0] = out;
    
    return out;
}

void    tTwoPole_setB0(tTwoPole* const ft, float b0)
{
    _tTwoPole* f = *ft;
    f->b0 = b0;
}

void    tTwoPole_setA1(tTwoPole* const ft, float a1)
{
    _tTwoPole* f = *ft;
    f->a1 = a1;
}

void    tTwoPole_setA2(tTwoPole* const ft, float a2)
{
    _tTwoPole* f = *ft;
    f->a2 = a2;
}


void    tTwoPole_setResonance(tTwoPole* const ft, float frequency, float radius, oBool normalize)
{
    _tTwoPole* f = *ft;
    
    if (frequency < 0.0f)   frequency = 0.0f;
    if (frequency > (leaf.sampleRate * 0.49f))   frequency = leaf.sampleRate * 0.49f;
    if (radius < 0.0f)      radius = 0.0f;
    if (radius >= 1.0f)     radius = 0.999999f;
    
    f->radius = radius;
    f->frequency = frequency;
    f->normalize = normalize;
    
    f->a2 = radius * radius;
    f->a1 =  -2.0f * radius * cosf(frequency * leaf.twoPiTimesInvSampleRate);
    
    if ( normalize )
    {
        // Normalize the filter gain ... not terribly efficient.
        float real = 1 - radius + (f->a2 - radius) * cosf(2 * frequency * leaf.twoPiTimesInvSampleRate);
        float imag = (f->a2 - radius) * sinf(2 * frequency * leaf.twoPiTimesInvSampleRate);
        f->b0 = sqrtf( powf(real, 2) + powf(imag, 2) );
    }
}

void    tTwoPole_setCoefficients(tTwoPole* const ft, float b0, float a1, float a2)
{
    _tTwoPole* f = *ft;
    f->b0 = b0;
    f->a1 = a1;
    f->a2 = a2;
}

void    tTwoPole_setGain(tTwoPole* const ft, float gain)
{
    _tTwoPole* f = *ft;
    f->gain = gain;
}

void     tTwoPoleSampleRateChanged (tTwoPole* const ft)
{
    _tTwoPole* f = *ft;
    
    f->a2 = f->radius * f->radius;
    f->a1 =  -2.0f * f->radius * cosf(f->frequency * leaf.twoPiTimesInvSampleRate);
    
    if ( f->normalize )
    {
        // Normalize the filter gain ... not terribly efficient.
        float real = 1 - f->radius + (f->a2 - f->radius) * cosf(2 * f->frequency * leaf.twoPiTimesInvSampleRate);
        float imag = (f->a2 - f->radius) * sinf(2 * f->frequency * leaf.twoPiTimesInvSampleRate);
        f->b0 = sqrtf( powf(real, 2) + powf(imag, 2) );
    }
}

// ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ OneZero Filter ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ //
void    tOneZero_init(tOneZero* const ft, float theZero)
{
    _tOneZero* f = *ft = (_tOneZero*) leaf_alloc(sizeof(_tOneZero));
    
    f->gain = 1.0f;
    f->lastIn = 0.0f;
    f->lastOut = 0.0f;
    tOneZero_setZero(ft, theZero);
}

void    tOneZero_free(tOneZero* const ft)
{
    _tOneZero* f = *ft;
    
    leaf_free(f);
}

void    tOneZero_initToPool     (tOneZero* const ft, float theZero, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tOneZero* f = *ft = (_tOneZero*) mpool_alloc(sizeof(_tOneZero), m);
    
    f->gain = 1.0f;
    f->lastIn = 0.0f;
    f->lastOut = 0.0f;
    tOneZero_setZero(ft, theZero);
}

void    tOneZero_freeFromPool   (tOneZero* const ft, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tOneZero* f = *ft;
    
    mpool_free(f, m);
}

float   tOneZero_tick(tOneZero* const ft, float input)
{
    _tOneZero* f = *ft;
    
    float in = input * f->gain;
    float out = f->b1 * f->lastIn + f->b0 * in;
    
    f->lastIn = in;
    
    return out;
}

void    tOneZero_setZero(tOneZero* const ft, float theZero)
{
    _tOneZero* f = *ft;
    
    if (theZero > 0.0f) f->b0 = 1.0f / (1.0f + theZero);
    else                f->b0 = 1.0f / (1.0f - theZero);
    
    f->b1 = -theZero * f->b0;
    
}

void    tOneZero_setB0(tOneZero* const ft, float b0)
{
    _tOneZero* f = *ft;
    
    f->b0 = b0;
}

void    tOneZero_setB1(tOneZero* const ft, float b1)
{
    _tOneZero* f = *ft;
    f->b1 = b1;
}

void    tOneZero_setCoefficients(tOneZero* const ft, float b0, float b1)
{
    _tOneZero* f = *ft;
    f->b0 = b0;
    f->b1 = b1;
}

void    tOneZero_setGain(tOneZero *ft, float gain)
{
    _tOneZero* f = *ft;
    f->gain = gain;
}

float   tOneZero_getPhaseDelay(tOneZero* const ft, float frequency )
{
    _tOneZero* f = *ft;
    
    if ( frequency <= 0.0f) frequency = 0.05f;
    
    f->frequency = frequency;
    
    float omegaT = 2 * PI * frequency * leaf.invSampleRate;
    float real = 0.0, imag = 0.0;
    
    real += f->b0;
    
    real += f->b1 * cosf(omegaT);
    imag -= f->b1 * sinf(omegaT);
    
    real *= f->gain;
    imag *= f->gain;
    
    float phase = atan2f( imag, real );
    
    real = 0.0; imag = 0.0;
    
    phase -= atan2f( imag, real );
    
    phase = fmodf( -phase, 2 * PI );
    
    return phase / omegaT;
}

// ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ TwoZero Filter ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ //
void    tTwoZero_init(tTwoZero* const ft)
{
    _tTwoZero* f = *ft = (_tTwoZero*) leaf_alloc(sizeof(_tTwoZero));
    
    f->gain = 1.0f;
    f->lastIn[0] = 0.0f;
    f->lastIn[1] = 0.0f;
}

void    tTwoZero_free(tTwoZero* const ft)
{
    _tTwoZero* f = *ft;
    
    leaf_free(f);
}

void    tTwoZero_initToPool     (tTwoZero* const ft, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tTwoZero* f = *ft = (_tTwoZero*) mpool_alloc(sizeof(_tTwoZero), m);
    
    f->gain = 1.0f;
    f->lastIn[0] = 0.0f;
    f->lastIn[1] = 0.0f;
}

void    tTwoZero_freeFromPool   (tTwoZero* const ft, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tTwoZero* f = *ft;
    
    mpool_free(f, m);
}

float   tTwoZero_tick(tTwoZero* const ft, float input)
{
    _tTwoZero* f = *ft;
    
    float in = input * f->gain;
    float out = f->b2 * f->lastIn[1] + f->b1 * f->lastIn[0] + f->b0 * in;
    
    f->lastIn[1] = f->lastIn[0];
    f->lastIn[0] = in;
    
    return out;
}

void    tTwoZero_setNotch(tTwoZero* const ft, float freq, float radius)
{
    _tTwoZero* f = *ft;
    
    // Should also deal with frequency being > half sample rate / nyquist. See STK
    if (freq < 0.0f)    freq = 0.0f;
    if (radius < 0.0f)  radius = 0.0f;
    
    f->frequency = freq;
    f->radius = radius;
    
    f->b2 = radius * radius;
    f->b1 = -2.0f * radius * cosf(freq * leaf.twoPiTimesInvSampleRate); // OPTIMIZE with LOOKUP or APPROXIMATION
    
    // Normalize the filter gain. From STK.
    if ( f->b1 > 0.0f ) // Maximum at z = 0.
        f->b0 = 1.0f / ( 1.0f + f->b1 + f->b2 );
    else            // Maximum at z = -1.
        f->b0 = 1.0f / ( 1.0f - f->b1 + f->b2 );
    f->b1 *= f->b0;
    f->b2 *= f->b0;
    
}

void    tTwoZero_setB0(tTwoZero* const ft, float b0)
{
    _tTwoZero* f = *ft;
    f->b0 = b0;
}

void    tTwoZero_setB1(tTwoZero* const ft, float b1)
{
    _tTwoZero* f = *ft;
    f->b1 = b1;
}

void    tTwoZero_setCoefficients(tTwoZero* const ft, float b0, float b1, float b2)
{
    _tTwoZero* f = *ft;
    f->b0 = b0;
    f->b1 = b1;
    f->b2 = b2;
}

void    tTwoZero_setGain(tTwoZero* const ft, float gain)
{
    _tTwoZero* f = *ft;
    f->gain = gain;
}

void tTwoZeroSampleRateChanged(tTwoZero* const ft)
{
    _tTwoZero* f = *ft;
    
    tTwoZero_setNotch(ft, f->frequency, f->radius);
}

// ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ PoleZero Filter ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ //
void   tPoleZero_init(tPoleZero* const pzf)
{
    _tPoleZero* f = *pzf = (_tPoleZero*) leaf_alloc(sizeof(_tPoleZero));
    
    f->gain = 1.0f;
    f->b0 = 1.0;
    f->a0 = 1.0;
    
    f->lastIn = 0.0f;
    f->lastOut = 0.0f;
}

void   tPoleZero_free(tPoleZero* const pzf)
{
    _tPoleZero* f = *pzf;
    
    leaf_free(f);
}

void    tPoleZero_initToPool        (tPoleZero* const pzf, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tPoleZero* f = *pzf = (_tPoleZero*) mpool_alloc(sizeof(_tPoleZero), m);
    
    f->gain = 1.0f;
    f->b0 = 1.0;
    f->a0 = 1.0;
    
    f->lastIn = 0.0f;
    f->lastOut = 0.0f;
}

void    tPoleZero_freeFromPool      (tPoleZero* const pzf, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tPoleZero* f = *pzf;
    
    mpool_free(f, m);
}

void    tPoleZero_setB0(tPoleZero* const pzf, float b0)
{
    _tPoleZero* f = *pzf;
    f->b0 = b0;
}

void    tPoleZero_setB1(tPoleZero* const pzf, float b1)
{
    _tPoleZero* f = *pzf;
    f->b1 = b1;
}

void    tPoleZero_setA1(tPoleZero* const pzf, float a1)
{
    _tPoleZero* f = *pzf;
    
    if (a1 >= 1.0f) // a1 should be less than 1.0
    {
        a1 = 0.999999f;
    }
    
    f->a1 = a1;
}

void    tPoleZero_setCoefficients(tPoleZero* const pzf, float b0, float b1, float a1)
{
    _tPoleZero* f = *pzf;
    
    if (a1 >= 1.0f) // a1 should be less than 1.0
    {
        a1 = 0.999999f;
    }
    
    f->b0 = b0;
    f->b1 = b1;
    f->a1 = a1;
}

void    tPoleZero_setAllpass(tPoleZero* const pzf, float coeff)
{
    _tPoleZero* f = *pzf;
    
    if (coeff >= 1.0f) // allpass coefficient >= 1.0 makes filter unstable
    {
        coeff = 0.999999f;
    }
    
    f->b0 = coeff;
    f->b1 = 1.0f;
    f->a0 = 1.0f;
    f->a1 = coeff;
}

void    tPoleZero_setBlockZero(tPoleZero* const pzf, float thePole)
{
    _tPoleZero* f = *pzf;
    
    if (thePole >= 1.0f) // allpass coefficient >= 1.0 makes filter unstable
    {
        thePole = 0.999999f;
    }
    
    f->b0 = 1.0f;
    f->b1 = -1.0f;
    f->a0 = 1.0f;
    f->a1 = -thePole;
}

void    tPoleZero_setGain(tPoleZero* const pzf, float gain)
{
    _tPoleZero* f = *pzf;
    f->gain = gain;
}

float   tPoleZero_tick(tPoleZero* const pzf, float input)
{
    _tPoleZero* f = *pzf;
    
    float in = input * f->gain;
    float out = (f->b0 * in) + (f->b1 * f->lastIn) - (f->a1 * f->lastOut);
    
    f->lastIn = in;
    f->lastOut = out;
    
    return out;
}

// ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ BiQuad Filter ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ //
void    tBiQuad_init(tBiQuad* const ft)
{
    _tBiQuad* f = *ft = (_tBiQuad*) leaf_alloc(sizeof(_tBiQuad));
    
    f->gain = 1.0f;
    
    f->b0 = 0.0f;
    f->a0 = 0.0f;
    
    f->lastIn[0] = 0.0f;
    f->lastIn[1] = 0.0f;
    f->lastOut[0] = 0.0f;
    f->lastOut[1] = 0.0f;
}

void    tBiQuad_free(tBiQuad* const ft)
{
    _tBiQuad* f = *ft;
    
    leaf_free(f);
}

void    tBiQuad_initToPool     (tBiQuad* const ft, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tBiQuad* f = *ft = (_tBiQuad*) mpool_alloc(sizeof(_tBiQuad), m);
    
    f->gain = 1.0f;
    
    f->b0 = 0.0f;
    f->a0 = 0.0f;
    
    f->lastIn[0] = 0.0f;
    f->lastIn[1] = 0.0f;
    f->lastOut[0] = 0.0f;
    f->lastOut[1] = 0.0f;
}

void    tBiQuad_freeFromPool   (tBiQuad* const ft, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tBiQuad* f = *ft;
    
    mpool_free(f, m);
}

float   tBiQuad_tick(tBiQuad* const ft, float input)
{
    _tBiQuad* f = *ft;
    
    float in = input * f->gain;
    float out = f->b0 * in + f->b1 * f->lastIn[0] + f->b2 * f->lastIn[1];
    out -= f->a2 * f->lastOut[1] + f->a1 * f->lastOut[0];
    
    f->lastIn[1] = f->lastIn[0];
    f->lastIn[0] = in;
    
    f->lastOut[1] = f->lastOut[0];
    f->lastOut[0] = out;
    
    return out;
}

void    tBiQuad_setResonance(tBiQuad* const ft, float freq, float radius, oBool normalize)
{
    _tBiQuad* f = *ft;
    
    if (freq < 0.0f)    freq = 0.0f;
    if (freq > (leaf.sampleRate * 0.49f))   freq = leaf.sampleRate * 0.49f;
    if (radius < 0.0f)  radius = 0.0f;
    if (radius >= 1.0f)  radius = 1.0f;
    
    f->frequency = freq;
    f->radius = radius;
    f->normalize = normalize;
    
    f->a2 = radius * radius;
    f->a1 = -2.0f * radius * cosf(freq * leaf.twoPiTimesInvSampleRate);
    
    if (normalize)
    {
        f->b0 = 0.5f - 0.5f * f->a2;
        f->b1 = 0.0f;
        f->b2 = -f->b0;
    }
}

void    tBiQuad_setNotch(tBiQuad* const ft, float freq, float radius)
{
    _tBiQuad* f = *ft;
    
    if (freq < 0.0f)    freq = 0.0f;
    if (freq > (leaf.sampleRate * 0.49f))   freq = leaf.sampleRate * 0.49f;
    if (radius < 0.0f)  radius = 0.0f;
    
    f->b2 = radius * radius;
    f->b1 = -2.0f * radius * cosf(freq * leaf.twoPiTimesInvSampleRate); // OPTIMIZE with LOOKUP or APPROXIMATION
    
    // Does not attempt to normalize filter gain.
}

void tBiQuad_setEqualGainZeros(tBiQuad* const ft)
{
    _tBiQuad* f = *ft;
    f->b0 = 1.0f;
    f->b1 = 0.0f;
    f->b2 = -1.0f;
}

void    tBiQuad_setB0(tBiQuad* const ft, float b0)
{
    _tBiQuad* f = *ft;
    f->b0 = b0;
}

void    tBiQuad_setB1(tBiQuad* const ft, float b1)
{
    _tBiQuad* f = *ft;
    f->b1 = b1;
}

void    tBiQuad_setB2(tBiQuad* const ft, float b2)
{
    _tBiQuad* f = *ft;
    f->b2 = b2;
}

void    tBiQuad_setA1(tBiQuad* const ft, float a1)
{
    _tBiQuad* f = *ft;
    f->a1 = a1;
}

void    tBiQuad_setA2(tBiQuad* const ft, float a2)
{
    _tBiQuad* f = *ft;
    f->a2 = a2;
}

void    tBiQuad_setCoefficients(tBiQuad* const ft, float b0, float b1, float b2, float a1, float a2)
{
    _tBiQuad* f = *ft;
    f->b0 = b0;
    f->b1 = b1;
    f->b2 = b2;
    f->a1 = a1;
    f->a2 = a2;
}

void    tBiQuad_setGain(tBiQuad* const ft, float gain)
{
    _tBiQuad* f = *ft;
    f->gain = gain;
}

void    tBiQuadSampleRateChanged(tBiQuad* const ft)
{
    _tBiQuad* f = *ft;
    f->a2 = f->radius * f->radius;
    f->a1 = -2.0f * f->radius * cosf(f->frequency * leaf.twoPiTimesInvSampleRate);
    
    if (f->normalize)
    {
        f->b0 = 0.5f - 0.5f * f->a2;
        f->b1 = 0.0f;
        f->b2 = -f->b0;
    }
}

// Less efficient, more accurate version of SVF, in which cutoff frequency is taken as floating point Hz value and tanf
// is calculated when frequency changes.
void tSVF_init(tSVF* const svff, SVFType type, float freq, float Q)
{
    tSVF_initToPool     (svff, type, freq, Q, &leaf.mempool);

    // or maybe this?
    /*
     * hp=1 bp=A/Q (where A is 10^(G/40) and G is gain in decibels) and lp = 1
     */

}

void tSVF_free(tSVF* const svff)
{
    tSVF_freeFromPool     (svff, &leaf.mempool);
}

void    tSVF_initToPool     (tSVF* const svff, SVFType type, float freq, float Q, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tSVF* svf = *svff = (_tSVF*) mpool_alloc(sizeof(_tSVF), m);

    svf->type = type;

    svf->ic1eq = 0;
    svf->ic2eq = 0;
    svf->Q = Q;
    svf->cutoff = freq;
    svf->g = tanf(PI * freq * leaf.invSampleRate);
    svf->k = 1.0f/Q;
    svf->a1 = 1.0f/(1.0f + svf->g * (svf->g + svf->k));
    svf->a2 = svf->g*svf->a1;
    svf->a3 = svf->g*svf->a2;
    svf->cH = 0.0f;
    svf->cB = 0.0f;
    svf->cL = 1.0f;

    if (type == SVFTypeLowpass)
    {
        svf->cH = 0.0f;
        svf->cB = 0.0f;
        svf->cBK = 0.0f;
        svf->cL = 1.0f;
    }
    else if (type == SVFTypeBandpass)
    {
        svf->cH = 0.0f;
        svf->cB = 1.0f;
        svf->cBK = 0.0f;
        svf->cL = 0.0f;
    }

    else if (type == SVFTypeHighpass)
    {
        svf->cH = 1.0f;
        svf->cB = 0.0f;
        svf->cBK = -1.0f;
        svf->cL = -1.0f;
    }

    else if (type == SVFTypeNotch)
    {
        svf->cH = 1.0f;
        svf->cB = 0.0f;
        svf->cBK = -1.0f;
        svf->cL = 0.0f;
    }


    else if (type == SVFTypePeak)
    {
        svf->cH = 1.0f;
        svf->cB = 0.0f;
        svf->cBK = -1.0f;
        svf->cL = -2.0f;
    }
}

void    tSVF_freeFromPool   (tSVF* const svff, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tSVF* svf = *svff;
    
    mpool_free(svf, m);
}

float   tSVF_tick(tSVF* const svff, float v0)
{
    _tSVF* svf = *svff;
    
    float v1,v2,v3;
    v3 = v0 - svf->ic2eq;
    v1 = (svf->a1 * svf->ic1eq) + (svf->a2 * v3);
    v2 = svf->ic2eq + (svf->a2 * svf->ic1eq) + (svf->a3 * v3);
    svf->ic1eq = (2.0f * v1) - svf->ic1eq;
    svf->ic2eq = (2.0f * v2) - svf->ic2eq;
    
    return (v0 * svf->cH) + (v1 * svf->cB) + (svf->k * v1 * svf->cBK) + (v2 * svf->cL);
}

void     tSVF_setFreq(tSVF* const svff, float freq)
{
    _tSVF* svf = *svff;
    svf->cutoff = freq;
    svf->g = tanf(PI * freq * leaf.invSampleRate);
    svf->a1 = 1.0f/(1.0f + svf->g * (svf->g + svf->k));
    svf->a2 = svf->g * svf->a1;
    svf->a3 = svf->g * svf->a2;
}

void     tSVF_setQ(tSVF* const svff, float Q)
{
    _tSVF* svf = *svff;
    svf->Q = Q;
    svf->k = 1.0f/Q;

    svf->a1 = 1.0f/(1.0f + svf->g * (svf->g + svf->k));
    svf->a2 = svf->g * svf->a1;
    svf->a3 = svf->g * svf->a2;
}

void    tSVF_setFreqAndQ(tSVF* const svff, float freq, float Q)
{
    _tSVF* svf = *svff;
    svf->k = 1.0f/Q;
    svf->g = tanf(PI * freq * leaf.invSampleRate);
    svf->a1 = 1.0f/(1.0f + svf->g * (svf->g + svf->k));
    svf->a2 = svf->g * svf->a1;
    svf->a3 = svf->g * svf->a2;
}

// Efficient version of tSVF where frequency is set based on 12-bit integer input for lookup in tanh wavetable.
void   tEfficientSVF_init(tEfficientSVF* const svff, SVFType type, uint16_t input, float Q)
{
    _tEfficientSVF* svf = *svff = (_tEfficientSVF*) leaf_alloc(sizeof(_tEfficientSVF));
    
    svf->type = type;
    
    svf->ic1eq = 0;
    svf->ic2eq = 0;
    
    svf->g = __leaf_table_filtertan[input];
    svf->k = 1.0f/Q;
    svf->a1 = 1.0f/(1.0f+svf->g*(svf->g+svf->k));
    svf->a2 = svf->g*svf->a1;
    svf->a3 = svf->g*svf->a2;
}

void tEfficientSVF_free(tEfficientSVF* const svff)
{
    _tEfficientSVF* svf = *svff;
    
    leaf_free(svf);
}

void    tEfficientSVF_initToPool    (tEfficientSVF* const svff, SVFType type, uint16_t input, float Q, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tEfficientSVF* svf = *svff = (_tEfficientSVF*) mpool_alloc(sizeof(_tEfficientSVF), m);
    
    svf->type = type;
    
    svf->ic1eq = 0;
    svf->ic2eq = 0;
    
    svf->g = __leaf_table_filtertan[input];
    svf->k = 1.0f/Q;
    svf->a1 = 1.0f/(1.0f+svf->g*(svf->g+svf->k));
    svf->a2 = svf->g*svf->a1;
    svf->a3 = svf->g*svf->a2;
}

void    tEfficientSVF_freeFromPool  (tEfficientSVF* const svff, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tEfficientSVF* svf = *svff;
    
    mpool_free(svf, m);
}

float   tEfficientSVF_tick(tEfficientSVF* const svff, float v0)
{
    _tEfficientSVF* svf = *svff;
    
    float v1,v2,v3;
    v3 = v0 - svf->ic2eq;
    v1 = (svf->a1 * svf->ic1eq) + (svf->a2 * v3);
    v2 = svf->ic2eq + (svf->a2 * svf->ic1eq) + (svf->a3 * v3);
    svf->ic1eq = (2.0f * v1) - svf->ic1eq;
    svf->ic2eq = (2.0f * v2) - svf->ic2eq;
    
    if (svf->type == SVFTypeLowpass)        return v2;
    else if (svf->type == SVFTypeBandpass)  return v1;
    else if (svf->type == SVFTypeHighpass)  return v0 - (svf->k * v1) - v2;
    else if (svf->type == SVFTypeNotch)     return v0 - (svf->k * v1);
    else if (svf->type == SVFTypePeak)      return v0 - (svf->k * v1) - (2.0f * v2);
    else                                    return 0.0f;
    
}

void     tEfficientSVF_setFreq(tEfficientSVF* const svff, uint16_t input)
{
    _tEfficientSVF* svf = *svff;
    
    svf->g = __leaf_table_filtertan[input];
    svf->a1 = 1.0f/(1.0f + svf->g * (svf->g + svf->k));
    svf->a2 = svf->g * svf->a1;
    svf->a3 = svf->g * svf->a2;
}

void     tEfficientSVF_setQ(tEfficientSVF* const svff, float Q)
{
    _tEfficientSVF* svf = *svff;
    
    svf->k = 1.0f/Q;
    svf->a1 = 1.0f/(1.0f + svf->g * (svf->g + svf->k));
    svf->a2 = svf->g * svf->a1;
    svf->a3 = svf->g * svf->a2;
}

/* Highpass */
void    tHighpass_init(tHighpass* const ft, float freq)
{
    _tHighpass* f = *ft = (_tHighpass*) leaf_alloc(sizeof(_tHighpass));
    
    f->R = (1.0f - (freq * leaf.twoPiTimesInvSampleRate));
    f->ys = 0.0f;
    f->xs = 0.0f;
    
    f->frequency = freq;
}

void    tHighpass_free(tHighpass* const ft)
{
    _tHighpass* f = *ft;
    
    leaf_free(f);
}

void    tHighpass_initToPool    (tHighpass* const ft, float freq, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tHighpass* f = *ft = (_tHighpass*) mpool_calloc(sizeof(_tHighpass), m);
    
    f->R = (1.0f - (freq * leaf.twoPiTimesInvSampleRate));
    f->ys = 0.0f;
    f->xs = 0.0f;
    
    f->frequency = freq;
}

void    tHighpass_freeFromPool  (tHighpass* const ft, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tHighpass* f = *ft;
    
    mpool_free(f, m);
}

void     tHighpass_setFreq(tHighpass* const ft, float freq)
{
    _tHighpass* f = *ft;
    f->frequency = freq;
    f->R = (1.0f - (freq * leaf.twoPiTimesInvSampleRate));
    
}

float     tHighpass_getFreq(tHighpass* const ft)
{
    _tHighpass* f = *ft;
    return f->frequency;
}

// From JOS DC Blocker
float   tHighpass_tick(tHighpass* const ft, float x)
{
    _tHighpass* f = *ft;
    f->ys = x - f->xs + f->R * f->ys;
    f->xs = x;
    return f->ys;
}

void tHighpassSampleRateChanged(tHighpass* const ft)
{
    _tHighpass* f = *ft;
    f->R = (1.0f-((f->frequency * 2.0f * 3.14f) * leaf.invSampleRate));
}

void tButterworth_init(tButterworth* const ft, int N, float f1, float f2)
{
    _tButterworth* f = *ft = (_tButterworth*) leaf_alloc(sizeof(_tButterworth));
    
    f->f1 = f1;
    f->f2 = f2;
    f->gain = 1.0f;
    
    f->N = N;
    
    if (f->N > NUM_SVF_BW) f->N = NUM_SVF_BW;
    
    for(int i = 0; i < N/2; ++i)
    {
        tSVF_init(&f->low[i], SVFTypeHighpass, f1, 0.5f/cosf((1.0f+2.0f*i)*PI/(2*N)));
        tSVF_init(&f->high[i], SVFTypeLowpass, f2, 0.5f/cosf((1.0f+2.0f*i)*PI/(2*N)));
    }
}

void tButterworth_free(tButterworth* const ft)
{
    _tButterworth* f = *ft;
    
    for(int i = 0; i < f->N/2; ++i)
    {
        tSVF_free(&f->low[i]);
        tSVF_free(&f->high[i]);
    }
    
    leaf_free(f);
}

void    tButterworth_initToPool     (tButterworth* const ft, int N, float f1, float f2, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tButterworth* f = *ft = (_tButterworth*) mpool_alloc(sizeof(_tButterworth), m);
    
    f->f1 = f1;
    f->f2 = f2;
    f->gain = 1.0f;
    
    f->N = N;
    
    if (f->N > NUM_SVF_BW) f->N = NUM_SVF_BW;
    
    for(int i = 0; i < N/2; ++i)
    {
        tSVF_initToPool(&f->low[i], SVFTypeHighpass, f1, 0.5f/cosf((1.0f+2.0f*i)*PI/(2*N)), mp);
        tSVF_initToPool(&f->high[i], SVFTypeLowpass, f2, 0.5f/cosf((1.0f+2.0f*i)*PI/(2*N)), mp);
    }
}

void    tButterworth_freeFromPool   (tButterworth* const ft, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tButterworth* f = *ft;
    
    for(int i = 0; i < f->N/2; ++i)
    {
        tSVF_freeFromPool(&f->low[i], mp);
        tSVF_freeFromPool(&f->high[i], mp);
    }
    
    mpool_free(f, m);
}

float tButterworth_tick(tButterworth* const ft, float samp)
{
    _tButterworth* f = *ft;
    
    for(int i = 0; i < ((f->N)/2); ++i)
    {
        samp = tSVF_tick(&f->low[i],samp);
        samp = tSVF_tick(&f->high[i],samp);
    }
    return samp;
}

void tButterworth_setF1(tButterworth* const ft, float f1)
{
    _tButterworth* f = *ft;
    
    f->f1 = f1;
    for(int i = 0; i < ((f->N)/2); ++i)        tSVF_setFreq(&f->low[i], f1);
}

void tButterworth_setF2(tButterworth* const ft, float f2)
{
    _tButterworth* f = *ft;
    
    f->f2 = f2;
    for(int i = 0; i < ((f->N)/2); ++i)        tSVF_setFreq(&f->high[i], f2);
}

void tButterworth_setFreqs(tButterworth* const ft, float f1, float f2)
{
    _tButterworth* f = *ft;
    
    f->f1 = f1;
    f->f2 = f2;
    for(int i = 0; i < ((f->N)/2); ++i)
    {
        tSVF_setFreq(&f->low[i], f1);
        tSVF_setFreq(&f->high[i], f2);
    }
}

void	tFIR_init(tFIR* const firf, float* coeffs, int numTaps)
{
    _tFIR* fir = *firf = (_tFIR*) leaf_alloc(sizeof(_tFIR));
    
    fir->numTaps = numTaps;
    fir->coeff = coeffs;
    fir->past = (float*)leaf_alloc(sizeof(float) * fir->numTaps);
    for (int i = 0; i < fir->numTaps; ++i) fir->past[i] = 0.0f;
}

void    tFIR_free(tFIR* const firf)
{
    _tFIR* fir = *firf;
    
    leaf_free(fir->past);
    leaf_free(fir);
}

void    tFIR_initToPool     (tFIR* const firf, float* coeffs, int numTaps, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tFIR* fir = *firf = (_tFIR*) mpool_alloc(sizeof(_tFIR), m);
    
    fir->numTaps = numTaps;
    fir->coeff = coeffs;
    fir->past = (float*) mpool_alloc(sizeof(float) * fir->numTaps, m);
    for (int i = 0; i < fir->numTaps; ++i) fir->past[i] = 0.0f;
}

void    tFIR_freeFromPool   (tFIR* const firf, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tFIR* fir = *firf;
    
    mpool_free(fir->past, m);
    mpool_free(fir, m);
}

float	tFIR_tick(tFIR* const firf, float input)
{
    _tFIR* fir = *firf;
    
	fir->past[0] = input;
	float y = 0.0f;
	for (int i = 0; i < fir->numTaps; ++i) y += fir->past[i]*fir->coeff[i];
	for (int i = fir->numTaps-1; i > 0; --i) fir->past[i] = fir->past[i-1];
	return y;
}

//---------------------------------------------
////
/// Median filter implemented based on James McCartney's median filter in Supercollider,
/// translated from a Gen~ port of the Supercollider code that I believe was made by Rodrigo Costanzo and which I got from PA Tremblay - JS


void    tMedianFilter_init           (tMedianFilter* const f, int size)
{
	tMedianFilter_initToPool(f, size, &leaf.mempool);
}
void    tMedianFilter_free           (tMedianFilter* const f)
{
	tMedianFilter_freeFromPool(f, &leaf.mempool);
}
void    tMedianFilter_initToPool     (tMedianFilter* const mf, int size, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tMedianFilter* f = *mf = (_tMedianFilter*) mpool_alloc(sizeof(_tMedianFilter), m);

    f->size = size;
    f->middlePosition = size / 2;
	f->last = size - 1;
	f->pos = -1;
    f->val = (float*) mpool_alloc(sizeof(float) * size, m);
    f->age = (int*) mpool_alloc(sizeof(int) * size, m);
    for (int i = 0; i < f->size; ++i)
	{
    	f->val[i] = 0.0f;
    	f->age[i] = i;
	}

}
void    tMedianFilter_freeFromPool   (tMedianFilter* const mf, tMempool* const mp)
{
    _tMempool* m = *mp;
    _tMedianFilter* f = *mf;

    mpool_free(f->val, m);
    mpool_free(f->age, m);
    mpool_free(f, m);
}

float   tMedianFilter_tick           (tMedianFilter* const mf, float input)
{
	_tMedianFilter* f = *mf;

	for(int i=0; i<f->size; i++) {
		int thisAge = f->age[i];
		if(thisAge == f->last) {
			f->pos = i;
		}
		else {
			thisAge++;
			f->age[i] = thisAge;
		}
	}

	while( f->pos!=0 ) {
		float test = f->val[f->pos-1];
		if(input < test) {
			f->val[f->pos]=test;
			f->age[f->pos]=f->age[f->pos-1];
			f->pos -= 1;
		} else {break;}
	}

	while(f->pos != f->last) {
		float test = f->val[f->pos+1];
		if( input > test) {
			f->val[f->pos] = test;
			f->age[f->pos] = f->age[f->pos+1];
			f->pos += 1;
		} else {break;}
	}

	f->val[f->pos] = input;
	f->age[f->pos] = 0;

	return  f->val[f->middlePosition];
}

/////

void    tVZFilter_init           (tVZFilter* const vf, VZFilterType type, float freq, float bandWidth)
{
	tVZFilter_initToPool(vf, type, freq, bandWidth, &leaf.mempool);
}

void    tVZFilter_free           (tVZFilter* const vf)
{
	tVZFilter_freeFromPool(vf, &leaf.mempool);
}
void    tVZFilter_initToPool     (tVZFilter* const vf, VZFilterType type, float freq, float bandWidth, tMempool* const mp)
{

	 _tMempool* m = *mp;
	 _tVZFilter* f = *vf = (_tVZFilter*) mpool_alloc(sizeof(_tVZFilter), m);
	f->fc   = freq;
	f->type = type;
	f->G    = ONE_OVER_SQRT2;
	f->invG    = 1.0f/ONE_OVER_SQRT2;
	f->B    = bandWidth;
	f->m    = 0.0f;
	f->s1 = 0.0f;
	f->s2 = 0.0f;
	f->sr = leaf.sampleRate;
	f->inv_sr = leaf.invSampleRate;
	tVZFilter_calcCoeffs(vf);


}
void    tVZFilter_freeFromPool   (tVZFilter* const vf, tMempool* const mp)
{
	 _tMempool* m = *mp;
		 _tVZFilter* f = *vf = (_tVZFilter*) mpool_alloc(sizeof(_tVZFilter), m);
		 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;

	float yL, yB, yH;

	// compute highpass output via Eq. 5.1:
	yH = (in - f->R2*f->s1 - f->g*f->s1 - f->s2) * f->h;

	// compute bandpass output by applying 1st integrator to highpass output:
	yB = tanhf(f->g*yH) + f->s1;
	f->s1 = f->g*yH + yB; // state update in 1st integrator

	// compute lowpass output by applying 2nd integrator to bandpass output:
	yL = tanhf(f->g*yB) + f->s2;
	f->s2 = f->g*yB + yL; // state update in 2nd integrator

	//according to the Vadim paper, we could add saturation to this model by adding a tanh in the integration stage.
	//
	//seems like that might look like this:
    // y = tanh(g*x) + s; // output computation
    // s = g*x + y; // state update

	//instead of this:
    // y = g*x + s; // output computation
    // s = g*x + y; // state update

	return f->cL*yL + f->cB*yB + f->cH*yH;

}

float   tVZFilter_tickEfficient           	(tVZFilter* const vf, float in)
{
	_tVZFilter* f = *vf;

	float yL, yB, yH;

	// compute highpass output via Eq. 5.1:
	yH = (in - f->R2*f->s1 - f->g*f->s1 - f->s2) * f->h;

	// compute bandpass output by applying 1st integrator to highpass output:
	yB = (f->g*yH) + f->s1;
	f->s1 = f->g*yH + yB; // state update in 1st integrator

	// compute lowpass output by applying 2nd integrator to bandpass output:
	yL = (f->g*yB) + f->s2;
	f->s2 = f->g*yB + yL; // state update in 2nd integrator

	//according to the Vadim paper, we could add saturation to this model by adding a tanh in the integration stage.
	//
	//seems like that might look like this:
    // y = tanh(g*x) + s; // output computation
    // s = g*x + y; // state update

	//instead of this:
    // y = g*x + s; // output computation
    // s = g*x + y; // state update

	return f->cL*yL + f->cB*yB + f->cH*yH;

}

float   tVZFilter_tickEfficientBP           	(tVZFilter* const vf, float in)
{
	_tVZFilter* f = *vf;

	float yL, yB, yH;

	// compute highpass output via Eq. 5.1:
	yH = (in - f->R2*f->s1 - f->g*f->s1 - f->s2) * f->h;

	// compute bandpass output by applying 1st integrator to highpass output:
	yB = (f->g*yH) + f->s1;
	f->s1 = f->g*yH + yB; // state update in 1st integrator

	// compute lowpass output by applying 2nd integrator to bandpass output:
	yL = (f->g*yB) + f->s2;
	f->s2 = f->g*yB + yL; // state update in 2nd integrator

	//according to the Vadim paper, we could add saturation to this model by adding a tanh in the integration stage.
	//
	//seems like that might look like this:
    // y = tanh(g*x) + s; // output computation
    // s = g*x + y; // state update

	//instead of this:
    // y = g*x + s; // output computation
    // s = g*x + y; // state update

	return f->cL*yL + f->cB*yB + f->cH*yH;

}


void   tVZFilter_calcCoeffs           (tVZFilter* const vf)
{

	_tVZFilter* f = *vf;
	f->g = tanf(PI * f->fc * f->inv_sr);  // embedded integrator gain (Fig 3.11)

	  switch( f->type )
	  {
	  case Bypass:
		{
		  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;
		}
		break;
	  case Lowpass:
		{
			f->R2 = f->invG;
			f->cL = 1.0f; f->cB = 0.0f; f->cH = 0.0f;
		}
		break;
	  case Highpass:
		{
			f->R2 = f->invG;
			f->cL = 0.0f; f->cB = 0.0f; f->cH = 1.0f;
		}
		break;
	  case BandpassSkirt:
		{
			f->R2 = f->invG;
			f->cL = 0.0f; f->cB = 1.0f; f->cH = 0.0f;
		}
		break;
	  case BandpassPeak:
		{
			f->R2 = 2.0f*tVZFilter_BandwidthToR(vf, f->B);
			f->cL = 0.0f; f->cB = f->R2; f->cH = 0.0f;
		}
		break;
	  case BandReject:
		{
			f->R2 = 2.0f*tVZFilter_BandwidthToR(vf, f->B);
			f->cL = 1.0f; f->cB = 0.0f; f->cH = 1.0f;
		}
		break;
	  case Bell:
		{
			float fl = f->fc*powf(2.0f, (-f->B)*0.5f); // lower bandedge frequency (in Hz)
			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
			f->R2 = 2.0f*sqrtf(((r*r+1.0f)/r-2.0f)/(4.0f*f->G));
			f->cL = 1.0f; f->cB = f->R2*f->G; f->cH = 1.0f;
		}
		break;
	  case Lowshelf:
		{
			float A = sqrtf(f->G);
		  f->g /= sqrtf(A);               // scale SVF-cutoff frequency for shelvers
		  f->R2 = 2*sinhf(f->B*logf(2.0f)*0.5f);
		  f->cL = f->G; f->cB = f->R2*A; f->cH = 1.0f;
		}
		break;
	  case Highshelf:
		{
		  float A = sqrtf(f->G);
		  f->g *= sqrtf(A);               // scale SVF-cutoff frequency for shelvers
		  f->R2 = 2.0f*sinhf(f->B*logf(2.0f)*0.5f);
		  f->cL = 1.0f; f->cB = f->R2*A; f->cH = f->G;
		}
		break;
	  case Allpass:
		{
			f->R2 = 2.0f*tVZFilter_BandwidthToR(vf, f->B);
			f->cL = 1.0f; f->cB = -f->R2; f->cH = 1.0f;
		}
		break;

		// experimental - maybe we must find better curves for cL, cB, cH:
	  case Morph:
		{
			f->R2 = f->invG;
		  float x  = 2.0f*f->m-1.0f;

		  f->cL = maximum(-x, 0.0f); /*cL *= cL;*/
		  f->cH = minimum( x, 0.0f); /*cH *= cH;*/
		  f->cB = 1.0f-x*x;

			// bottom line: we need to test different versions for how they feel when tweaking the
			// morph parameter

		  // this scaling ensures constant magnitude at the cutoff point (we divide the coefficients by
		  // the magnitude response value at the cutoff frequency and scale back by the gain):
		  float s = f->G * sqrtf((f->R2*f->R2) / (f->cL*f->cL + f->cB*f->cB + f->cH*f->cH - 2.0f*f->cL*f->cH));
		  f->cL *= s; f->cB *= s; f->cH *= s;
		}
		break;

	  }

	  f->h = 1.0f / (1.0f + f->R2*f->g + f->g*f->g);  // factor for feedback precomputation
}


void   tVZFilter_setBandwidth           	(tVZFilter* const vf, float B)
{
	_tVZFilter* f = *vf;
	f->B = LEAF_clip(0.0f, B, 100.0f);
	tVZFilter_calcCoeffs(vf);
}
void   tVZFilter_setFreq           (tVZFilter* const vf, float freq)
{
	_tVZFilter* f = *vf;
	f->fc = LEAF_clip(0.0f, freq, 0.5f*leaf.sampleRate);
	tVZFilter_calcCoeffs(vf);
}
void   tVZFilter_setFreqAndBandwidth           (tVZFilter* const vf, float freq, float bw)
{
	_tVZFilter* f = *vf;
	f->B = LEAF_clip(0.0f,bw, 100.0f);
	f->fc = LEAF_clip(0.0f, freq, 0.5f*leaf.sampleRate);
	tVZFilter_calcCoeffs(vf);
}

void   tVZFilter_setGain          		(tVZFilter* const vf, float gain)
{
	_tVZFilter* f = *vf;
	f->G = LEAF_clip(0.000001f, gain, 100.0f);
	f->invG = 1.0f/f->G;
	tVZFilter_calcCoeffs(vf);
}

void   tVZFilter_setMorph          		(tVZFilter* const vf, float morph)
{
	_tVZFilter* f = *vf;
	f->m = LEAF_clip(0.0f, morph, 1.0f);
	tVZFilter_calcCoeffs(vf);
}

void tVZFilter_setType         		(tVZFilter* const vf, VZFilterType type)
{
	_tVZFilter* f = *vf;
	f->type = type;
	tVZFilter_calcCoeffs(vf);
}

float tVZFilter_BandwidthToR(tVZFilter* const vf, float B)
{
	_tVZFilter* f = *vf;
  float fl = f->fc*powf(2.0f, -B*0.5f); // lower bandedge frequency (in Hz)
  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);

}