fuzz_pair.hpp

Go to the documentation of this file.

#pragma once

 
#include <pulp/signal/dc_blocker.hpp>
#include <pulp/signal/junction.hpp>
#include <pulp/signal/denormal.hpp>
#include <pulp/signal/oversampling.hpp>
#include <pulp/signal/rng.hpp>
#include <pulp/signal/units.hpp>
 
#include <algorithm>
#include <cmath>
#include <cstdint>
 
namespace pulp::signal {

enum class FuzzDevice : std::uint8_t {
    germanium,  
    silicon,    
};

template <typename SampleType = float>
class FuzzPairT {
public:
    using Device = FuzzDevice;
 
    // ── Device rows ───────────────────────────────────────────────────────
    // Ideality factor sets knee SOFTNESS; saturation current sets where the
    // knee falls. Both are honest [design parameter]s standing for a device
    // family's published order of magnitude, never a measured part.

    static constexpr double kGermaniumIdeality = 1.15;
    static constexpr double kGermaniumSaturation = 5e-6;
    static constexpr double kSiliconIdeality = 1.03;
    static constexpr double kSiliconSaturation = 5e-12;

    static constexpr double kThermalVoltage = junction::kThermalVoltage;
 
    // ── Circuit constants ─────────────────────────────────────────────────

    static constexpr double kNominalCollectorCurrent = 1.0;

    static constexpr double kQuiescentCollector = 0.5;

    static constexpr double kEmitterDegeneration = 0.18;

    static constexpr double kStarveExponent = 1.6;

    static constexpr double kMinAvailableCurrent = 1e-3;

    static constexpr double kInputImpedanceKohm = 68.0;

    static constexpr double kFeedbackConductanceMin = 0.05;
    static constexpr double kFeedbackConductanceMax = 1.0;

    static constexpr double kCouplingGain = 0.85;

    static constexpr double kLoopGainCeiling = 0.94;

    static constexpr int kSolverIterations = 24;

    static constexpr double kNewtonStepThermal = 2.0;

    static constexpr double kResidualTolerance = 1e-6;

    static constexpr double kInputDriveRails = 0.4;

    static constexpr double kGermaniumDriftOctaves = 0.6;
    static constexpr double kSiliconDriftOctaves = 0.05;

    static constexpr double kDcBlockerCornerHz = 12.0;

    static constexpr int kOversampleFactor = 4;
 
    // ── Lifecycle ─────────────────────────────────────────────────────────
 
    FuzzPairT() { update(); }

    void prepare(double sample_rate) {
        sample_rate_ = sample_rate > 0.0 ? sample_rate : sample_rate_;
        using Os = OversamplerT<SampleType>;
        oversampler_.set_kind(Os::Kind::linear_phase_fir);
        oversampler_.set_quality(Os::Quality::standard);
        oversampler_.set_factor(Os::Factor::x4);
        oversampler_.set_sample_rate(static_cast<SampleType>(sample_rate_));
        latency_ = oversampler_.latency_samples();
 
        // The blocker runs at the OVERSAMPLED rate, so its pole is set for
        // that rate: a corner low enough to pass the bass but high enough to
        // clear the DC step a starved operating point introduces.
        blocker_.set_pole(static_cast<SampleType>(
            1.0 - kDcBlockerCornerHz * 6.283185307179586 /
                      (sample_rate_ * kOversampleFactor)));
        drift_.prepare(sample_rate_ * kOversampleFactor);
        update();
        reset();
    }
 
    void reset() {
        collector1_ = quiescent1_;
        collector2_ = quiescent2_;
        blocker_.reset();
        drift_.reset();
        oversampler_.reset();
        worst_residual_ = 0.0;
    }
 
    // ── Controls ──────────────────────────────────────────────────────────
 
    void set_device(Device device) {
        device_ = device;
        update();
    }
 
    Device device() const { return device_; }

    void set_fuzz(double amount01) {
        if (!std::isfinite(amount01)) return;
        fuzz_ = std::clamp(amount01, 0.0, 1.0);
        update();
    }

    void set_bias_starve(double starve01) {
        if (!std::isfinite(starve01)) return;
        starve_ = std::clamp(starve01, 0.0, 1.0);
        update();
    }

    void set_source_impedance_kohm(double kohm) {
        if (!std::isfinite(kohm)) return;
        source_kohm_ = std::clamp(kohm, 0.1, 1000.0);
        update();
    }
 
    void set_output_level_db(double db) {
        if (!std::isfinite(db)) return;
        output_gain_ = units::db_to_linear(std::clamp(db, -24.0, 24.0));
    }
 
    void set_mix(double mix01) {
        if (!std::isfinite(mix01)) return;
        mix_ = std::clamp(mix01, 0.0, 1.0);
    }

    void set_oversampling_enabled(bool on) { oversample_ = on; }
 
    void set_seed(std::uint32_t seed) { drift_.set_seed(seed); }

    void set_drift_enabled(bool on) { drift_enabled_ = on; }
 
    // ── Observables (for metering and for tests) ──────────────────────────

    double loading_factor() const { return loading_; }

    double bias_voltage() const { return knee2_; }

    double base_bias_voltage() const { return base_bias2_; }

    double quiescent_collector() const { return quiescent2_; }

    double stage_gain() const { return stage_gain_; }

    double input_scale_volts() const { return input_scale_; }

    double available_current() const { return available_current_; }

    double loop_gain() const {
        const double stage1 = collector_load_ * (kNominalCollectorCurrent /
                                                 (ideality_ * kThermalVoltage)) /
                              degeneration1_;
        const double stage2 =
            collector_load_ * (available_current_ / (ideality_ * kThermalVoltage)) / degeneration2_;
        return feedback_conductance_ * feedback_norm_ * kCouplingGain * kCouplingGain * stage1 *
               stage2 * loading_;
    }

    double worst_residual() const { return worst_residual_; }

    int latency_samples() const noexcept { return oversample_ ? latency_ : 0; }
 
    // ── Processing ────────────────────────────────────────────────────────
 
    SampleType process(SampleType input) {
        // Every stage below is recursive. Treat a non-finite sample as a
        // recovery boundary before it reaches either the oversampler or the
        // coupled operating-point state.
        if (!std::isfinite(static_cast<double>(input))) {
            reset();
            return SampleType{0};
        }
        const double dry = static_cast<double>(input);
        const auto stage = [this](SampleType x) {
            return static_cast<SampleType>(solve_stage(static_cast<double>(x)));
        };
        const double wet = oversample_ ? static_cast<double>(oversampler_.process(input, stage))
                                       : static_cast<double>(stage(input));
        const double output = (dry + mix_ * (wet - dry)) * output_gain_;
        // Also recover if an already-poisoned composed primitive or state ever
        // reaches this boundary. `reset()` is RT-safe and clears every recursive
        // state owned by the pair.
        if (!std::isfinite(output)) {
            reset();
            return SampleType{0};
        }
        return static_cast<SampleType>(snap_to_zero(output));
    }
 
    void process_block(const SampleType* in, SampleType* out, int n) {
        for (int i = 0; i < n; ++i) out[i] = process(in[i]);
    }
 
private:
    double junction_current(double v_be) const { return network_.current(v_be); }

    double transconductance(double v_be) const { return network_.conductance(v_be); }

    double collector(double current) const {
        return std::clamp(1.0 - collector_load_ * current, 0.0, 1.0);
    }
 
    struct PairState {
        double first;
        double second;
    };
 
    struct Residual {
        double first;
        double second;
 
        double infinity_norm() const {
            return std::max(std::abs(first), std::abs(second));
        }
        double euclidean_norm() const { return std::hypot(first, second); }
    };
 
    struct PairEvaluation {
        PairState current;
        PairState collector_voltage;
        Residual residual;
    };
 
    PairEvaluation evaluate_pair(PairState voltage, double scale, double drive,
                                 double feedback) const {
        const PairState current{scale * junction_current(voltage.first),
                                scale * junction_current(voltage.second)};
        const PairState collector_voltage{collector(current.first), collector(current.second)};
        const Residual residual{
            voltage.first + kEmitterDegeneration * current.first -
                (base_bias1_ + drive +
                 feedback * (collector_voltage.second - quiescent2_) * kCouplingGain),
            voltage.second + kEmitterDegeneration * current.second -
                (base_bias2_ + (collector_voltage.first - quiescent1_) * kCouplingGain)};
        return {current, collector_voltage, residual};
    }
 
    void accept_if_better(PairState candidate, double scale, double drive,
                          double feedback, PairState& best, double& best_merit) const {
        const double merit =
            evaluate_pair(candidate, scale, drive, feedback).residual.euclidean_norm();
        if (merit < best_merit) {
            best = candidate;
            best_merit = merit;
        }
    }

    bool solve_iteration(PairState& voltage, double scale, double drive,
                         double feedback, double step_limit) const {
        const PairEvaluation point = evaluate_pair(voltage, scale, drive, feedback);
        if (point.residual.infinity_norm() < kResidualTolerance) return true;
 
        const double gm1 = scale * transconductance(voltage.first);
        const double gm2 = scale * transconductance(voltage.second);
        // A clamped collector has zero coupling derivative; emitter
        // degeneration remains active on the diagonal.
        const double slope1 = point.collector_voltage.first > 0.0 &&
                                      point.collector_voltage.first < 1.0
                                  ? collector_load_ * gm1
                                  : 0.0;
        const double slope2 = point.collector_voltage.second > 0.0 &&
                                      point.collector_voltage.second < 1.0
                                  ? collector_load_ * gm2
                                  : 0.0;
        const double j11 = 1.0 + kEmitterDegeneration * gm1;
        const double j22 = 1.0 + kEmitterDegeneration * gm2;
        const double j12 = feedback * kCouplingGain * slope2;
        const double j21 = kCouplingGain * slope1;
        const double determinant = j11 * j22 - j12 * j21;
        if (std::abs(determinant) <= 1e-12) return false;
 
        const PairState newton_step{
            (j22 * point.residual.first - j12 * point.residual.second) / determinant,
            (j11 * point.residual.second - j21 * point.residual.first) / determinant};
        const double newton_scale =
            std::min(1.0, step_limit /
                              std::max({std::abs(newton_step.first),
                                        std::abs(newton_step.second), 1e-30}));
 
        const PairState gradient{
            j11 * point.residual.first + j21 * point.residual.second,
            j12 * point.residual.first + j22 * point.residual.second};
        const double gradient_scale =
            step_limit /
            std::max({std::abs(gradient.first), std::abs(gradient.second), 1e-30});
 
        PairState best = voltage;
        double best_merit = point.residual.euclidean_norm();
        const PairState full_newton{voltage.first - newton_scale * newton_step.first,
                                    voltage.second - newton_scale * newton_step.second};
        const double full_newton_merit =
            evaluate_pair(full_newton, scale, drive, feedback).residual.euclidean_norm();
        if (full_newton_merit < best_merit) {
            // Preserve Newton's quadratic path in its local basin.
            voltage = full_newton;
            return false;
        }
 
        // At a rail kink Newton's active-set Jacobian can point across the
        // kink. A fixed ladder along both Newton and residual-gradient
        // directions keeps the cost bounded while retaining a descent option.
        double damping = 1.0;
        for (int trial = 0; trial < 6; ++trial) {
            accept_if_better({voltage.first - damping * newton_scale * newton_step.first,
                              voltage.second - damping * newton_scale * newton_step.second},
                             scale, drive, feedback, best, best_merit);
            accept_if_better({voltage.first - damping * gradient_scale * gradient.first,
                              voltage.second - damping * gradient_scale * gradient.second},
                             scale, drive, feedback, best, best_merit);
            damping *= 0.5;
        }
        voltage = best;
        return false;
    }

    double solve_stage(double input) {
        // Thermal drift rides on the saturation current, which is the
        // temperature-sensitive quantity — germanium's far more so than
        // silicon's.
        double saturation = saturation_;
        if (drift_enabled_) {
            const double octaves = device_ == Device::germanium ? kGermaniumDriftOctaves
                                                                : kSiliconDriftOctaves;
            saturation *= std::pow(2.0, octaves * (static_cast<double>(drift_.next()) - 1.0));
        }
        const double scale = saturation / saturation_;
 
        const double drive = input * input_scale_ * loading_;
 
        // The feedback path carries the source-impedance loading for the same
        // reason the input does: it is the same node. One ratio, two
        // consequences — and because it multiplies the loop here, the ceiling
        // `loop_gain()` reports is the one the solver actually runs at.
        const double beta = feedback_conductance_ * feedback_norm_ * loading_;
 
        // Warm start. Both stages' base voltages are explicit given the PREVIOUS
        // sample's collectors, so the per-stage equation `v + R_e·i(v) = base`
        // can be pre-solved approximately before the coupled iteration begins.
        // Without it the solve starts wherever the last sample left it, which on
        // a transient is far enough away that the fixed iteration count does not
        // recover: a junction Newton descends by only about one thermal voltage
        // per step while it is deep in conduction, so an initial guess 400 mV
        // high needs a dozen iterations to walk back.
        //
        // TWO estimates are taken and the smaller kept, because each is a bound
        // that is tight in the regime the other is loose in, and both sit ABOVE
        // the root (the base-voltage function is convex, so Newton from above
        // descends monotonically — starting below it can overshoot into the
        // steep region instead):
        //
        //   * the linearised inverse — divide the base excursion by `1 + R_e·gm`
        //     at the bias point — which is exact for small excursions and
        //     increasingly conservative for large ones;
        //   * the conduction estimate `θ·ln(base/(R_e·Is))`, which assumes the
        //     whole base voltage is dropped across the degeneration. That is the
        //     right asymptote once the stage is driven hard, and it is what
        //     rescues the deeply starved settings: there the bias sits near
        //     cutoff, `1 + R_e·gm` is barely above 1, and the linear estimate
        //     alone lands hundreds of millivolts too high.
        //
        // The estimate is asked for at `R_e·scale` rather than `R_e` so the
        // drift-perturbed saturation current is the one it inverts — the ratio
        // it forms is `base/(R_e·Is·scale)`, which is the quantity that matters.
        const double base1 =
            base_bias1_ + drive + beta * (collector2_ - quiescent2_) * kCouplingGain;
        const double base2 = base_bias2_ + (collector1_ - quiescent1_) * kCouplingGain;
        const double degeneration_scale = kEmitterDegeneration * scale;
        PairState voltage{
            std::min(knee1_ + (base1 - base_bias1_) / degeneration1_,
                     network_.conduction_estimate(base1, degeneration_scale)),
            std::min(knee2_ + (base2 - base_bias2_) / degeneration2_,
                     network_.conduction_estimate(base2, degeneration_scale))};
 
        // Per-iteration step limit, in thermal voltages. The standard remedy
        // for Newton on a junction: unlimited, a step taken where the
        // exponential is steep lands somewhere it can take many iterations to
        // climb back from.
        const double step_limit = kNewtonStepThermal * ideality_ * kThermalVoltage;
        bool converged = false;
        for (int iter = 0; iter < kSolverIterations; ++iter) {
            // The loop bound remains compile-time fixed for RT predictability;
            // a converged sample simply skips further transcendental work.
            if (!converged)
                converged = solve_iteration(voltage, scale, drive, beta, step_limit);
        }
 
        const PairEvaluation solved = evaluate_pair(voltage, scale, drive, beta);
 
        // The only state carried to the next sample. The junction voltages are
        // NOT kept: the warm start re-derives them from these collectors, so
        // holding them too would be a second copy of the same operating point
        // that could disagree with this one after a parameter change.
        collector1_ = solved.collector_voltage.first;
        collector2_ = solved.collector_voltage.second;
        worst_residual_ = std::max(worst_residual_, solved.residual.infinity_norm());
 
        // The output is Q2's collector deviation — bounded to the rail by
        // construction, which is also where the hard fuzz clipping comes from.
        return static_cast<double>(
            blocker_.process(static_cast<SampleType>(collector2_ - quiescent2_)));
    }
 
    void update() {
        if (device_ == Device::germanium) {
            ideality_ = kGermaniumIdeality;
            saturation_ = kGermaniumSaturation;
        } else {
            ideality_ = kSiliconIdeality;
            saturation_ = kSiliconSaturation;
        }
        network_.ideality = ideality_;
        network_.saturation_current = saturation_;
        network_.leg_a = 1.0;
        // A transistor's base-emitter junction conducts one way only.
        network_.leg_b = 0.0;
 
        available_current_ =
            std::max(kNominalCollectorCurrent * std::pow(1.0 - starve_, kStarveExponent),
                     kMinAvailableCurrent);
 
        // The collector load puts the UNSTARVED quiescent collector mid-rail.
        collector_load_ = (1.0 - kQuiescentCollector) / kNominalCollectorCurrent;
 
        // One ratio, two consequences: attenuation and loop-gain loading.
        loading_ = kInputImpedanceKohm / (source_kohm_ + kInputImpedanceKohm);
 
        // ── The two operating points ──────────────────────────────────────
        //
        // Starvation is applied to STAGE 2 ONLY, because that is what it is: the
        // supply feeding Q2's collector node runs down. Stage 1 keeps its
        // healthy bias throughout.
        //
        // Starving both stages together is the tempting simplification and it
        // quietly removes two of the three behaviours the starve control exists
        // for. With both stages down, the loop gain falls as the SQUARE of the
        // current, so the feedback path is gone long before the operating point
        // is asymmetric enough to fold — yet the feedback path amplifying the
        // asymmetry is precisely the stated mechanism for octave-up. And a
        // uniformly quiet pair is not a gate: what gates is a healthy first
        // stage still swinging its full rail into a second stage that can no
        // longer follow it.
        //
        // The bias of each stage is its own conduction knee at its own current,
        // INVERTED FROM THE SAME LAW the forward path uses rather than written
        // out again — so the documented 363 mV germanium / 693 mV silicon
        // figures are consequences of the shipped `(n, Is)` rows and cannot
        // drift from them. The base sits one emitter drop above that knee.
        knee1_ = network_.knee_voltage(kNominalCollectorCurrent);
        knee2_ = network_.knee_voltage(available_current_);
        base_bias1_ = knee1_ + kEmitterDegeneration * kNominalCollectorCurrent;
        base_bias2_ = knee2_ + kEmitterDegeneration * available_current_;
 
        const double gm1 = kNominalCollectorCurrent / (ideality_ * kThermalVoltage);
        const double gm2 = available_current_ / (ideality_ * kThermalVoltage);
        degeneration1_ = 1.0 + gm1 * kEmitterDegeneration;
        degeneration2_ = 1.0 + gm2 * kEmitterDegeneration;
 
        // Where each collector actually rests. Starving raises Stage 2's toward
        // the supply rail, because less current means less drop across the load
        // — that asymmetric headroom is what folds the waveform and produces the
        // octave, so it must not be assumed away.
        //
        // Every coupling in the solve is taken as a deviation from THESE, not
        // from `kQuiescentCollector`. Referring them to the unstarved mid-rail
        // instead injects a DC step of `R_c·(I_nominal − I_available)` — nearly
        // half a rail at full starvation — into the next stage's base, which
        // drives the starved stage hard ON. That is backwards from what
        // starvation means, and it is also what makes the solve intractable
        // there: the junction is pushed hundreds of millivolts past its knee
        // every sample, where a fixed iteration count cannot follow.
        quiescent1_ = kQuiescentCollector;
        quiescent2_ = std::clamp(1.0 - collector_load_ * available_current_, 0.0, 1.0);
 
        // The stage gain, and with it the input scale, come from STAGE 1 at its
        // nominal operating point — the stage the input actually drives, and one
        // whose bias the starve control does not touch. Referring the drive to
        // the starved stage instead would raise it by exactly the factor
        // starvation lowered that stage's gain, and the starve control would do
        // almost nothing.
        stage_gain_ = collector_load_ * gm1 / degeneration1_;
        const double rail_to_rail = 1.0 / std::max(stage_gain_, 1e-12);
        input_scale_ = kInputDriveRails * rail_to_rail;
 
        feedback_conductance_ =
            kFeedbackConductanceMin + fuzz_ * (kFeedbackConductanceMax - kFeedbackConductanceMin);
 
        // The ceiling is enforced STRUCTURALLY rather than hoped for: the
        // feedback path is normalised so that the worst case over the whole
        // parameter space — maximum feedback, healthy bias, zero source
        // impedance — lands exactly on `kLoopGainCeiling`. The specification
        // describes this as "clamping feedback_conductance at the top of its
        // travel such that the measured maximum never reaches 1.0"; deriving
        // the normaliser is that clamp, with the arithmetic done once instead
        // of a magic number chosen by trial.
        //
        // The worst case uses the SMALLER ideality factor, because a smaller
        // one means a steeper exponential and therefore more transconductance
        // at the same current — silicon is the binding device, not germanium.
        // Degeneration narrows that margin (it divides out most of a `gm`
        // difference) without reversing it, so the binding device is unchanged.
        constexpr double worst_ideality = std::min(kGermaniumIdeality, kSiliconIdeality);
        const double worst_gm = kNominalCollectorCurrent / (worst_ideality * kThermalVoltage);
        const double worst_gm_eff = worst_gm / (1.0 + worst_gm * kEmitterDegeneration);
        const double worst_stage = collector_load_ * worst_gm_eff;
        const double worst_open_loop =
            kCouplingGain * kCouplingGain * worst_stage * worst_stage;
        feedback_norm_ =
            kLoopGainCeiling / std::max(worst_open_loop * kFeedbackConductanceMax, 1e-12);
    }
 
    Device device_ = Device::germanium;
    double sample_rate_ = 44100.0;
    double fuzz_ = 0.65;
    double starve_ = 0.0;
    double source_kohm_ = 10.0;
    double mix_ = 1.0;
    double output_gain_ = 1.0;
    bool oversample_ = true;
    bool drift_enabled_ = false;
 
    junction::JunctionPair network_{};
    double ideality_ = kGermaniumIdeality;
    double saturation_ = kGermaniumSaturation;
    double available_current_ = kNominalCollectorCurrent;
    double knee1_ = 0.0;
    double knee2_ = 0.0;
    double base_bias1_ = 0.0;
    double base_bias2_ = 0.0;
    double degeneration1_ = 1.0;
    double degeneration2_ = 1.0;
    double stage_gain_ = 1.0;
    double collector_load_ = 0.5;
    double loading_ = 1.0;
    double input_scale_ = 1.0;
    double feedback_conductance_ = 0.0;
    double feedback_norm_ = 1.0;
    int latency_ = 0;
 
    double quiescent1_ = kQuiescentCollector;
    double quiescent2_ = kQuiescentCollector;
    double collector1_ = kQuiescentCollector;
    double collector2_ = kQuiescentCollector;
    double worst_residual_ = 0.0;
 
    DcBlocker<SampleType> blocker_{};
    DriftT<double> drift_{};
    OversamplerT<SampleType> oversampler_{};
};
 
using FuzzPair = FuzzPairT<float>;
using FuzzPair64 = FuzzPairT<double>;
 
}  // namespace pulp::signal