distortion.hpp

Go to the documentation of this file.

#pragma once

 
#include <pulp/signal/denormal.hpp>
#include <pulp/signal/junction.hpp>
#include <pulp/signal/tpt_filter.hpp>
#include <pulp/signal/units.hpp>
 
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <utility>
 
namespace pulp::signal {

enum class DiodeModel : std::uint8_t {
    silicon,    
    germanium,  
    led,        
};

enum class ClipperTopology : std::uint8_t {
    to_ground,  
    in_loop,    
};
 
namespace detail {

inline void apply_diode_model(junction::JunctionPair& network, DiodeModel model) {
    switch (model) {
        case DiodeModel::germanium:
            network.saturation_current = 1e-6;
            network.ideality = 1.1;
            break;
        case DiodeModel::led:
            network.saturation_current = 1e-18;
            network.ideality = 2.0;
            break;
        case DiodeModel::silicon:
        default:
            network.saturation_current = 1e-14;
            network.ideality = 1.0;
            break;
    }
}

inline void apply_symmetry(junction::JunctionPair& network, double symmetry) {
    if (!std::isfinite(symmetry)) return;
    const double s = std::clamp(symmetry, -1.0, 1.0);
    network.leg_a = 1.0 / (1.0 + std::max(0.0, s));
    network.leg_b = 1.0 + std::min(0.0, s);
}
 
}  // namespace detail

struct ClipperSolverConfig;
 
namespace detail {

template <typename Residual, typename Derivative>
inline double solve_safeguarded_newton(Residual&& f, Derivative&& df, double lo, double hi,
                                       double initial, double alternate, double tolerance,
                                       int max_iterations, int& iterations_out) {
    // Order the bracket and pad it, so a root sitting exactly on an endpoint is
    // still strictly interior.
    if (lo > hi) std::swap(lo, hi);
    const double pad = 1e-9 + 0.01 * (hi - lo);
    lo -= pad;
    hi += pad;
 
    double f_lo = f(lo);
    double f_hi = f(hi);
    // If the physical bracket does not actually straddle the root, widen a
    // bounded number of times rather than failing. Bounded because an unbounded
    // widening loop is exactly the RT hazard the iteration cap exists to
    // prevent.
    for (int widen = 0; widen < 8 && f_lo * f_hi > 0.0; ++widen) {
        const double span = std::max(hi - lo, 1e-3);
        lo -= span;
        hi += span;
        f_lo = f(lo);
        f_hi = f(hi);
    }
 
    // Two candidate starts: where the solve landed last sample (excellent when
    // the signal is moving slowly, which is the common case) and the analytic
    // conduction estimate (excellent on a cold start or a fast transient).
    // Taking the better of the two by residual costs two evaluations and saves
    // the dozens of iterations a bad start would need.
    double v = std::clamp(initial, lo, hi);
    const double candidate = std::clamp(alternate, lo, hi);
    if (std::abs(f(candidate)) < std::abs(f(v))) v = candidate;
    int iterations = 0;
    for (; iterations < max_iterations; ++iterations) {
        const double fv = f(v);
        if (std::abs(fv) < tolerance) break;
 
        // Maintain the bracket from every evaluation, so the fallback keeps
        // getting tighter even while Newton is doing the work.
        if (f_lo * fv < 0.0) {
            hi = v;
            f_hi = fv;
        } else {
            lo = v;
            f_lo = fv;
        }
 
        const double dfv = df(v);
        double next = dfv > 0.0 ? v - fv / dfv : v;
        // Outside the bracket, or non-finite: bisect instead. This is what
        // makes the iteration globally convergent rather than merely usually
        // convergent.
        if (!std::isfinite(next) || next <= lo || next >= hi) next = 0.5 * (lo + hi);
        v = next;
    }
    iterations_out = iterations;
    return std::isfinite(v) ? v : 0.0;
}
 
}  // namespace detail
 
struct ClipperSolverConfig {
    static constexpr double kNewtonTolerance = 1e-9;

    static constexpr int kMaxNewtonIterations = 8;

    static constexpr int kMaxNewtonIterationsPerSample = 2 * kMaxNewtonIterations;

    static constexpr double kAdaaEpsilon = 1e-5;
};
 
namespace detail {

inline double solve_implicit_flow_stage(const junction::JunctionPair& network,
                                        double history_voltage, double history_flow,
                                        double source_current, double flow_weight,
                                        double c_over_step, double resistance,
                                        double initial, int& iterations_out) {
    const auto residual = [&](double v) {
        const double flow = source_current - v / resistance - network.current(v);
        return c_over_step * (v - history_voltage) -
               flow_weight * (history_flow + flow);
    };
    const auto derivative = [&](double v) {
        return c_over_step +
               flow_weight * (1.0 / resistance + network.conductance(v));
    };
 
    const double open_circuit = source_current * resistance;
    return solve_safeguarded_newton(
        residual, derivative,
        /*lo=*/std::min({0.0, history_voltage, open_circuit}),
        /*hi=*/std::max({0.0, history_voltage, open_circuit}),
        initial, network.conduction_estimate(open_circuit, resistance),
        ClipperSolverConfig::kNewtonTolerance / resistance,
        ClipperSolverConfig::kMaxNewtonIterations, iterations_out);
}

inline double solve_tr_bdf2(const junction::JunctionPair& network, double previous_voltage,
                            double previous_source_current, double source_current,
                            double capacitance, double timestep, double resistance,
                            int& iterations_out) {
    constexpr double gamma = 2.0 - 1.4142135623730950488016887242097;
    constexpr double bdf2_step = (1.0 - gamma) / (2.0 - gamma);
    constexpr double stage_weight = 1.0 / (gamma * (2.0 - gamma));
 
    const double source_at_stage =
        previous_source_current + gamma * (source_current - previous_source_current);
    const double previous_flow = previous_source_current - previous_voltage / resistance -
                                 network.current(previous_voltage);
    const double c_over_stage_t = capacitance / (gamma * timestep);
    int stage_iterations = 0;
    const double stage_voltage = solve_implicit_flow_stage(
        network, previous_voltage, previous_flow, source_at_stage,
        /*flow_weight=*/0.5, c_over_stage_t, resistance,
        /*initial=*/previous_voltage, stage_iterations);
 
    const double bdf2_history =
        stage_weight * stage_voltage + (1.0 - stage_weight) * previous_voltage;
    const double c_over_bdf2_t = capacitance / (bdf2_step * timestep);
    int final_iterations = 0;
    const double result = solve_implicit_flow_stage(
        network, bdf2_history, /*history_flow=*/0.0, source_current,
        /*flow_weight=*/1.0, c_over_bdf2_t, resistance,
        /*initial=*/stage_voltage, final_iterations);
    iterations_out = stage_iterations + final_iterations;
    return result;
}
 
}  // namespace detail

template <typename SampleType = float>
class DiodeClipperT {
public:
    using Model = DiodeModel;

    static constexpr double kResistanceDefault = 10000.0;

    static constexpr double kCapacitanceDefault = 51e-12;
 
    DiodeClipperT() {
        detail::apply_diode_model(network_, Model::silicon);
        detail::apply_symmetry(network_, 0.0);
    }
 
    void prepare(double sample_rate) {
        sample_rate_ = std::isfinite(sample_rate) && sample_rate > 0.0 ? sample_rate
                                                                       : sample_rate_;
        timestep_ = 1.0 / sample_rate_;
        reset();
    }
 
    void set_diode_model(Model model) { detail::apply_diode_model(network_, model); }

    void set_symmetry(double symmetry) { detail::apply_symmetry(network_, symmetry); }
 
    void set_resistance(double ohms) {
        if (std::isfinite(ohms)) resistance_ = std::max(ohms, 1.0);
    }

    void set_capacitance(double farads) {
        if (std::isfinite(farads)) capacitance_ = std::max(farads, 0.0);
    }

 
    void reset() {
        voltage_ = 0.0;
        previous_source_current_ = 0.0;
        previous_voltage_ = 0.0;
        last_iterations_ = 0;
    }

    int last_iteration_count() const { return last_iterations_; }
 
    SampleType process(SampleType input) {
        if (!std::isfinite(static_cast<double>(input))) {
            reset();
            return SampleType{0};
        }
        const double vin = static_cast<double>(input);
        const double v = solve(vin);
        if (!std::isfinite(v)) {
            reset();
            return SampleType{0};
        }
        previous_source_current_ = vin / resistance_;
        previous_voltage_ = v;
        voltage_ = v;
        return static_cast<SampleType>(snap_to_zero(v));
    }

    double voltage() const { return voltage_; }

    double resistive_residual(double vin, double v) const {
        return (vin - v) / resistance_ - network_.current(v);
    }
 
private:
    double solve(double vin) {
        // With no capacitor the ODE is not an ODE: `C·dv/dt` vanishes and the
        // circuit equation is purely algebraic, `(vin − v)/R = i_D(v)`. Pushing
        // that through trapezoidal rule would NOT solve it — TR at `C = 0`
        // enforces `g(vₙ) = −g(vₙ₋₁)`, so the residual alternates sign every
        // sample instead of vanishing. That is the textbook non-L-stable
        // behaviour of TR at a stiff limit, not a bug in the solve, and it is
        // why the degenerate case gets the algebraic equation directly.
        const bool memoryless = capacitance_ <= 0.0;
        if (!memoryless) {
            return detail::solve_tr_bdf2(network_, previous_voltage_, previous_source_current_,
                                         vin / resistance_, capacitance_, timestep_, resistance_,
                                         last_iterations_);
        }
 
        const auto residual = [&](double v) {
            return network_.current(v) - (vin - v) / resistance_;
        };
        // Every term is strictly positive, so the residual is strictly
        // increasing in v — which is what makes the bracketed fallback safe.
        const auto derivative = [&](double v) {
            return network_.conductance(v) + 1.0 / resistance_;
        };
 
        return detail::solve_safeguarded_newton(
            residual, derivative,
            /*lo=*/std::min({0.0, previous_voltage_, vin}),
            /*hi=*/std::max({0.0, previous_voltage_, vin}),
            /*initial=*/previous_voltage_,
            /*alternate=*/network_.conduction_estimate(vin, resistance_),
            ClipperSolverConfig::kNewtonTolerance / resistance_,
            ClipperSolverConfig::kMaxNewtonIterations, last_iterations_);
    }
 
    junction::JunctionPair network_{};
    double sample_rate_ = 44100.0;
    double timestep_ = 1.0 / 44100.0;
    double resistance_ = kResistanceDefault;
    double capacitance_ = kCapacitanceDefault;
 
    double voltage_ = 0.0;
    double previous_source_current_ = 0.0;
    double previous_voltage_ = 0.0;
    int last_iterations_ = 0;
};
 
using DiodeClipper = DiodeClipperT<float>;
using DiodeClipper64 = DiodeClipperT<double>;

template <typename SampleType = float>
class FeedbackClipperT {
public:
    using Model = DiodeModel;
    using Topology = ClipperTopology;

    static constexpr double kFeedbackResistanceDefault = 51000.0;

    static constexpr double kInputResistanceDefault = 4700.0;

    static constexpr double kKneeCornerHzDefault = 720.0;
 
    FeedbackClipperT() {
        detail::apply_diode_model(network_, Model::silicon);
        detail::apply_symmetry(network_, 0.0);
        update_capacitance();
    }
 
    void prepare(double sample_rate) {
        sample_rate_ = std::isfinite(sample_rate) && sample_rate > 0.0 ? sample_rate
                                                                       : sample_rate_;
        timestep_ = 1.0 / sample_rate_;
        ground_stage_.prepare(sample_rate_);
        reset();
    }
 
    void set_topology(Topology topology) { topology_ = topology; }
    Topology topology() const { return topology_; }
 
    void set_diode_model(Model model) {
        detail::apply_diode_model(network_, model);
        ground_stage_.set_diode_model(model);
    }
 
    void set_symmetry(double symmetry) {
        detail::apply_symmetry(network_, symmetry);
        ground_stage_.set_symmetry(symmetry);
    }
 
    void set_feedback_resistance(double ohms) {
        if (!std::isfinite(ohms)) return;
        feedback_resistance_ = std::max(ohms, 1.0);
        update_capacitance();
    }
 
    void set_input_resistance(double ohms) {
        if (std::isfinite(ohms)) input_resistance_ = std::max(ohms, 1.0);
    }
 
    void set_knee_corner_hz(double hz) {
        if (!std::isfinite(hz)) return;
        knee_corner_hz_ = std::clamp(hz, 1.0, 20000.0);
        update_capacitance();
    }
 

    double linear_gain() const { return feedback_resistance_ / input_resistance_; }
 
    void reset() {
        previous_source_current_ = 0.0;
        previous_voltage_ = 0.0;
        last_iterations_ = 0;
        ground_stage_.reset();
    }
 
    int last_iteration_count() const {
        return topology_ == Topology::to_ground ? ground_stage_.last_iteration_count()
                                                : last_iterations_;
    }
 
    SampleType process(SampleType input) {
        if (!std::isfinite(static_cast<double>(input))) {
            reset();
            return SampleType{0};
        }
        const double vin = static_cast<double>(input);
        if (topology_ == Topology::to_ground) {
            // A plain linear gain stage with the network shunting its output —
            // which IS the diode clipper fed by a pre-gain, so it composes one
            // rather than re-deriving the same circuit with the wrong
            // capacitor. That distinction is not cosmetic: the feedback cap
            // that sets the in-loop knee sits at 720 Hz, and reusing it here
            // would put a shunt clipper behind an audio-band filter, which
            // measurably tilts the supposedly level-independent plateau this
            // topology exists to provide.
            return ground_stage_.process(static_cast<SampleType>(vin * linear_gain()));
        }
 
        // In-loop: the same ODE with v ≡ −vout, R ≡ Rf, C ≡ Cf, and the input
        // current vin/Rin as an extra forcing term.
        const double forcing = vin / input_resistance_;
        const double v = solve(forcing);
        if (!std::isfinite(v)) {
            reset();
            return SampleType{0};
        }
        previous_source_current_ = forcing;
        previous_voltage_ = v;
        // v is −vout by the substitution, so the stage's output is its negation
        // — an inverting stage, which is what the topology is.
        return static_cast<SampleType>(snap_to_zero(-v));
    }
 
private:
    void update_capacitance() {
        constexpr double two_pi = 6.283185307179586476925286766559;
        capacitance_ = 1.0 / (two_pi * feedback_resistance_ * knee_corner_hz_);
    }
 
    double solve(double source_current) {
        return detail::solve_tr_bdf2(network_, previous_voltage_, previous_source_current_,
                                     source_current, capacitance_, timestep_,
                                     feedback_resistance_, last_iterations_);
    }
 
    junction::JunctionPair network_{};
    Topology topology_ = Topology::in_loop;
    double sample_rate_ = 44100.0;
    double timestep_ = 1.0 / 44100.0;
    double feedback_resistance_ = kFeedbackResistanceDefault;
    double input_resistance_ = kInputResistanceDefault;
    double knee_corner_hz_ = kKneeCornerHzDefault;
    double capacitance_ = 0.0;
 
    double previous_source_current_ = 0.0;
    double previous_voltage_ = 0.0;
    int last_iterations_ = 0;

    DiodeClipperT<SampleType> ground_stage_{};
};
 
using FeedbackClipper = FeedbackClipperT<float>;
using FeedbackClipper64 = FeedbackClipperT<double>;

template <typename SampleType = float>
class ToneStackT {
public:
    static constexpr double kPreGainDbMax = 12.0;
 
    void prepare(double sample_rate) {
        pre_.prepare(static_cast<SampleType>(sample_rate));
        post_.prepare(static_cast<SampleType>(sample_rate));
        reset();
    }
 
    void set_pre_tone_hz(double hz) {
        if (!std::isfinite(hz)) return;
        pre_hz_ = std::clamp(hz, 20.0, 20000.0);
        pre_.set_cutoff(static_cast<SampleType>(pre_hz_));
    }
 
    void set_post_tone_hz(double hz) {
        if (!std::isfinite(hz)) return;
        post_hz_ = std::clamp(hz, 20.0, 20000.0);
        post_.set_cutoff(static_cast<SampleType>(post_hz_));
    }

    void set_pre_gain_db(double db) {
        if (!std::isfinite(db)) return;
        pre_gain_ = units::db_to_linear(std::clamp(db, -kPreGainDbMax, kPreGainDbMax));
    }

    void set_tone_mix(double mix01) {
        if (std::isfinite(mix01)) mix_ = std::clamp(mix01, 0.0, 1.0);
    }
 
    void reset() {
        pre_.reset();
        post_.reset();
    }

    SampleType process_pre(SampleType input) {
        if (!std::isfinite(static_cast<double>(input))) {
            reset();
            return SampleType{0};
        }
        const auto out = pre_.process(input);
        const double shaped =
            static_cast<double>(out.lowpass) + pre_gain_ * static_cast<double>(out.highpass);
        const double result = static_cast<double>(input) +
                              mix_ * (shaped - static_cast<double>(input));
        if (!std::isfinite(result)) {
            reset();
            return SampleType{0};
        }
        return static_cast<SampleType>(result);
    }

    SampleType process_post(SampleType input) {
        if (!std::isfinite(static_cast<double>(input))) {
            reset();
            return SampleType{0};
        }
        const double filtered = static_cast<double>(post_.process_lowpass(input));
        const double result = static_cast<double>(input) +
                              mix_ * (filtered - static_cast<double>(input));
        if (!std::isfinite(result)) {
            reset();
            return SampleType{0};
        }
        return static_cast<SampleType>(result);
    }
 
private:
    TptFilterT<SampleType> pre_{};
    TptFilterT<SampleType> post_{};
    double pre_hz_ = 720.0;
    double post_hz_ = 4000.0;
    double pre_gain_ = 1.0;
    double mix_ = 1.0;
};
 
using ToneStack = ToneStackT<float>;
using ToneStack64 = ToneStackT<double>;
 
}  // namespace pulp::signal