desc:hard clipping with less aliasing
slider1:0<-6, 48, 1>Drive
slider2:1<0, 1, 1>Oversampling
slider3:1<0, 1, 1>Fancy Nonlinearity

@init
function apply(x0, c*)
local(y0)
instance(x1, x2, y1, y2)
(
    y0 = x0 + c.b1 * x1 + c.b2 * x2
            - c.a1 * y1 - c.a2 * y2;
    x2 = x1; x1 = x0;
    y2 = y1; y1 = y0;
    y0 * c.g;
);

eq0.g = 0.00155912;
eq1.g = 1.0;
eq2.g = 1.0;
eq3.g = 1.0;
eq4.g = 1.0;
eq5.g = 1.0;
eq6.g = 1.0;
eq7.g = 1.0;

eq0.b1 = 2.08319328;
eq0.b2 = 1.24778165;
eq0.a1 = -1.15588622;
eq0.a2 = 0.43252556;
eq1.b1 = 1.15573368;
eq1.b2 = 0.95390880;
eq1.a1 = -1.18077287;
eq1.a2 = 0.57137601;
eq2.b1 = -2.58088534;
eq2.b2 = 1.74249472;
eq2.a1 = -1.50190306;
eq2.a2 = 0.58565543;
eq3.b1 = -2.12178667;
eq3.b2 = 1.74699392;
eq3.a1 = -0.50043552;
eq3.a2 = 0.59393911;
eq4.b1 = -1.26659076;
eq4.b2 = 1.74391445;
eq4.a1 = -0.79694471;
eq4.a2 = 0.60465343;
eq5.b1 = 0.63756680;
eq5.b2 = 0.94955405;
eq5.a1 = -0.23035542;
eq5.a2 = 0.73111164;
eq6.b1 = 0.39615867;
eq6.b2 = 0.98029175;
eq6.a1 = -0.07856381;
eq6.a2 = 0.90186786;
eq7.b1 = 0.30122276;
eq7.b2 = 0.99782300;
eq7.a1 = 0.04415921;
eq7.a2 = 0.98139113;

compensation.g = 1.60769517;
compensation.b1 = -0.00000001;
compensation.b2 = -0.00000007;
compensation.a1 = 0.53589834;
compensation.a2 = 0.07179671;

function iir(x)
instance(q0, q1, q2, q3, q4, q5, q6, q7)
(
    x = q0.apply(x, eq0);
    x = q1.apply(x, eq1);
    x = q2.apply(x, eq2);
    x = q3.apply(x, eq3);
    x = q4.apply(x, eq4);
    x = q5.apply(x, eq5);
    x = q6.apply(x, eq6);
    x = q7.apply(x, eq7);
    x
);

function f(x)
(
    min(max(x, -1), 1);
);

function F0(x)
(
    abs(x) <= 1 ? x * x * (1/2) : x * sign(x) - (1/2)
);

function F1(x)
(
    abs(x) <= 1 ? x * x * x * (1/3) : (x * x * (1/2) - (1/6)) * sign(x)
);

function distort(x0)
local(y, den0, den1, term0, term1)
instance(x1, x2)
(
    slider3 ? (
        den0 = x0 - x1; den0 *= den0;
        den1 = x2 - x1; den1 *= den1;
        term0 = den0 < 0.0001
              ? 0.5 * f((x0 + 2 * x1) * (1/3))
              : (x0 * (F0(x0) - F0(x1)) - (F1(x0) - F1(x1))) / den0;
        term1 = den1 < 0.0001
              ? 0.5 * f((x2 + 2 * x1) * (1/3))
              : (x2 * (F0(x2) - F0(x1)) - (F1(x2) - F1(x1))) / den1;
        y = term0 + term1;
    ) : (
        y = f(x0);
    );
    x2 = x1; x1 = x0;
    y
);

function process(x)
instance(ix, ox, dd, ddc, d, dc)
(
    x = x * drive;
    slider2 ? (
        xa = ix.iir(x) * 2;
        xb = ix.iir(0) * 2;
        xa = dd.distort(xa);
        xb = dd.distort(xb);
        slider3 ? (
            xa = dc.apply(xa, compensation);
            xb = dc.apply(xb, compensation);
            x = ox.iir(xa);
                ox.iir(xb);
        ) : (
                ox.iir(xa);
            x = ox.iir(xb);
        );
    ) : (
        x = d.distort(x);
        slider3 ? (
            x = dc.apply(x, compensation);
        );
    );
    x / drive
);

pdc_bot_ch = 0;
pdc_top_ch = 2;
pdc_delay = slider2 ? 16 : 0;

@slider
drive = pow(10, (min(max(slider1, -6), 48) / 20));

@sample
spl0 = s0.process(spl0);
spl1 = s1.process(spl1);