TPT diode ladder

diode_ladder.h

#include <cmath>
#include <numbers>

// See: https://www.discodsp.net/VAFilterDesign_2.1.2.pdf (p. 170)

template<typename X>
struct Diode_ladder {
    
    struct Params {
        X fc = 1000;
        X k = 0; // 0...16
    };
    
    auto reset(float sr) -> void
    {
        sp = 1 / sr;
        _s1 = _s2 = _s3 = _s4 = 0;
        prepare({});
    }
    
    auto prepare(const Params& params) -> void
    {
        this->params = params;
        const auto g = std::tan(std::numbers::pi_v<X> * params.fc * sp);
        _G = g / (1 + g);
    }
    
    auto process(X x) -> X
    {
        const auto k = params.k;
        const auto G = _G;
        
        // scale states
        const auto gp = 1 - G; // 1 / (1 + g)
        const auto s1 = _s1*gp;
        const auto s2 = _s2*gp;
        const auto s3 = _s3*gp;
        const auto s4 = _s4*gp;
        
        const auto gh = 0.5f * G;
        const auto g34 = gh;
        const auto g23 = gh / (1 - gh*gh);
        const auto g12 = gh / (1 - gh*g23);
        const auto g01 = 2*gh / (1 - 2*gh*g12);
                          
        const auto s34 = s4;
        const auto s23 = (gh*s4 + s3) / (1 - gh*gh);
        const auto s12 = (gh*s23 + s2) / (1 - gh*g23);
        const auto s01 = (2*gh*s12 + s1) / (1 - 2*gh*g12);
        
        const auto g04 = g01*g12*g23*g34;
        const auto s04 = g12*g23*g34*s01 + g23*g34*s12 + g34*s23 + s34;
        const auto u = (x - k*s04) / (1 + k*g04);
        
        const auto y1 = g01*u + s01;
        const auto y2 = g12*y1 + s12;
        const auto y3 = g23*y2 + s23;
        const auto y4 = g34*y3 + s34;
        
        _s1 = 2*y1 - _s1;
        _s2 = 2*y2 - _s2;
        _s3 = 2*y3 - _s3;
        _s4 = 2*y4 - _s4;
        
        return y4;
    }
    
private:
    
    Params params{};
    
    X sp = 1 / 48000.f;
    X _G{};
    
    X _s1{};
    X _s2{};
    X _s3{};
    X _s4{};

};