diode_bridge_compressor.hpp

Go to the documentation of this file.

#pragma once

 
#include <pulp/signal/ballistics_filter.hpp>
#include <pulp/signal/dc_blocker.hpp>
#include <pulp/signal/denormal.hpp>
#include <pulp/signal/dynamics_core.hpp>
#include <pulp/signal/dynamics_contract.hpp>
#include <pulp/signal/junction.hpp>
#include <pulp/signal/tpt_filter.hpp>
#include <pulp/signal/units.hpp>
 
#include <algorithm>
#include <cmath>
#include <cstdint>
 
namespace pulp::signal {
 
// ── The current-controlled variable attenuator ────────────────────────────

template <typename SampleType = float>
class DiodeBridgeGainT {
public:
    // ── Design parameters (the complete roster) ───────────────────────────

    static constexpr double kIdeality = 1.9;

    static constexpr double kSaturationCurrent = 1e-14;

    static constexpr double kSeriesResistance = 1000.0;

    static constexpr double kBaseCurvature = 0.5;

    static constexpr double kDriveCurvature = 0.6;

    static constexpr double kDriveBridgeMax = 0.20;

    static constexpr double kDriveBridgeFloor = 0.3;

    static constexpr double kMaxControlDrive = 30.62;

    static constexpr double kMonotonicMargin = 0.98;

    static constexpr double kAdaaEpsilon = 1e-5;
 
    // ── Lifecycle ─────────────────────────────────────────────────────────
 
    DiodeBridgeGainT() { configure_junction(); }

    void prepare(double /*sample_rate*/) { reset(); }
 
    void reset() { previous_shaped_ = 0.0; }
 
    // ── Controls ──────────────────────────────────────────────────────────

    void set_character(double character01) {
        if (!std::isfinite(character01)) return;
        const double c = std::clamp(character01, 0.0, 1.0);
        drive_ = kDriveBridgeMax * (kDriveBridgeFloor + (1.0 - kDriveBridgeFloor) * c);
    }

    void set_adaa(bool on) { adaa_ = on; }

    double drive() const { return drive_; }
 
    // ── The control law ───────────────────────────────────────────────────

    static double gain_for_control_drive(double control_drive) {
        return 1.0 / (1.0 + std::max(0.0, control_drive));
    }

    static double control_drive_for_gain_db(double gain_reduction_db) {
        return units::db_to_linear(-std::min(0.0, gain_reduction_db)) - 1.0;
    }

    double control_drive_for_current(double amperes) const {
        const double bias = network_.knee_voltage(std::max(amperes, 0.0));
        return kSeriesResistance * network_.conductance(bias) - series_offset_;
    }

    double dynamic_resistance(double amperes) const {
        return 1.0 / network_.conductance(network_.knee_voltage(std::max(amperes, 0.0)));
    }
 
    // ── The colour ────────────────────────────────────────────────────────

    static double curvature(double control_drive) {
        return kBaseCurvature * (1.0 + kDriveCurvature * std::max(0.0, control_drive));
    }

    static double shape(double s, double beta) { return s - beta * s * s * s / 3.0; }

    static double shape_antiderivative(double s, double beta) {
        const double s2 = s * s;
        return 0.5 * s2 - beta * s2 * s2 / 12.0;
    }

    static double max_operating_amplitude(double beta) {
        return beta > 0.0 ? kMonotonicMargin / std::sqrt(beta) : 1e30;
    }

    static double third_harmonic_ratio(double s, double beta) {
        const double fundamental = 1.0 - beta * s * s / 4.0;
        if (!(fundamental > 0.0)) return 0.0;
        return (beta * s * s / 12.0) / fundamental;
    }
 
    // ── Processing ────────────────────────────────────────────────────────

    SampleType process(SampleType input, double control_drive) {
        if (!std::isfinite(static_cast<double>(input)) || !std::isfinite(control_drive)) {
            reset();
            return SampleType{0};
        }
        const double x = std::max(0.0, control_drive);
        const double beta = curvature(x);
        const double limit = max_operating_amplitude(beta);
        const double s = std::clamp(static_cast<double>(input) * drive_, -limit, limit);
 
        double shaped;
        if (adaa_) {
            const double delta = s - previous_shaped_;
            shaped = std::abs(delta) > kAdaaEpsilon
                         ? (shape_antiderivative(s, beta) -
                            shape_antiderivative(previous_shaped_, beta)) /
                               delta
                         : shape(0.5 * (s + previous_shaped_), beta);
        } else {
            shaped = shape(s, beta);
        }
        previous_shaped_ = s;
 
        // Divide the drive back out: the shaper contributes harmonics, never
        // level (series law 1). See the class doc block.
        const double output = snap_to_zero(gain_for_control_drive(x) * shaped / drive_);
        if (!std::isfinite(output)) {
            reset();
            return SampleType{0};
        }
        return static_cast<SampleType>(output);
    }
 
private:
    void configure_junction() {
        network_.saturation_current = kSaturationCurrent;
        network_.ideality = kIdeality;
        // The audio sees the balanced bridge as an antiparallel MATCHED pair —
        // symmetric, hence the odd-only curvature above. One diode per leg.
        network_.leg_a = 1.0;
        network_.leg_b = 1.0;
        // `conductance` at zero bias is not zero (it is `2·Is/θ`), so the drive
        // law is offset to make `x(0) = 0` exactly. Without this a zero control
        // current would ask for a gain of 1/(1 + 1.5e-10) rather than unity —
        // inaudible, but it would make "no current means no reduction" an
        // approximation instead of an identity.
        series_offset_ = kSeriesResistance * network_.conductance(0.0);
    }
 
    junction::JunctionPair network_{};
    double series_offset_ = 0.0;
    double drive_ = kDriveBridgeMax * (kDriveBridgeFloor + (1.0 - kDriveBridgeFloor) * 0.35);
    double previous_shaped_ = 0.0;
    bool adaa_ = true;
};
 
using DiodeBridgeGain = DiodeBridgeGainT<float>;
using DiodeBridgeGain64 = DiodeBridgeGainT<double>;
 
// ── The transformer bracket ───────────────────────────────────────────────

template <typename SampleType = float>
class TransformerBracketT {
public:
    static constexpr double kLowCornerHz = 20.0;

    static constexpr double kHighCornerHz = 28000.0;

    static constexpr double kSaturationDepth = 0.15;

    static constexpr double kEvenAsymmetry = 0.5;

    static constexpr double kMonotonicMargin = 0.98;

    static constexpr double kAdaaEpsilon = 1e-5;

    static constexpr double kPeakGain = 1.0;
 
    void prepare(double sample_rate) {
        highpass_.prepare(static_cast<SampleType>(sample_rate));
        lowpass_.prepare(static_cast<SampleType>(sample_rate));
        highpass_.set_cutoff(static_cast<SampleType>(kLowCornerHz));
        lowpass_.set_cutoff(static_cast<SampleType>(kHighCornerHz));
        reset();
    }
 
    void reset() {
        highpass_.reset();
        lowpass_.reset();
        previous_input_ = 0.0;
    }

    void set_character(double character01) {
        if (!std::isfinite(character01)) return;
        const double a = kSaturationDepth * std::clamp(character01, 0.0, 1.0);
        positive_ = a * (1.0 + kEvenAsymmetry);
        negative_ = a * (1.0 - kEvenAsymmetry);
        symmetric_ = a;
        asymmetric_ = a * kEvenAsymmetry;
        const double steepest = std::max(positive_, negative_);
        limit_ = steepest > 0.0 ? kMonotonicMargin / steepest : 1e30;
    }
 
    void set_adaa(bool on) { adaa_ = on; }

    double saturate(double u) const {
        return u - 0.5 * symmetric_ * u * std::abs(u) - 0.5 * asymmetric_ * u * u;
    }

    double saturate_antiderivative(double u) const {
        const double u2 = u * u;
        return 0.5 * u2 - symmetric_ * u2 * std::abs(u) / 6.0 - asymmetric_ * u2 * u / 6.0;
    }
 
    SampleType process(SampleType input) {
        if (!std::isfinite(static_cast<double>(input))) {
            reset();
            return SampleType{0};
        }
        const double u = std::clamp(static_cast<double>(input), -limit_, limit_);
 
        double shaped;
        if (adaa_) {
            const double delta = u - previous_input_;
            shaped = std::abs(delta) > kAdaaEpsilon
                         ? (saturate_antiderivative(u) - saturate_antiderivative(previous_input_)) /
                               delta
                         : saturate(0.5 * (u + previous_input_));
        } else {
            shaped = saturate(u);
        }
        previous_input_ = u;
 
        const SampleType blocked = highpass_.process_highpass(static_cast<SampleType>(shaped));
        const SampleType output = lowpass_.process_lowpass(blocked);
        if (!std::isfinite(static_cast<double>(output))) {
            reset();
            return SampleType{0};
        }
        return output;
    }
 
private:
    TptFilterT<SampleType> highpass_{};
    TptFilterT<SampleType> lowpass_{};
    double symmetric_ = kSaturationDepth * 0.35;
    double asymmetric_ = kSaturationDepth * 0.35 * kEvenAsymmetry;
    double positive_ = kSaturationDepth * 0.35 * (1.0 + kEvenAsymmetry);
    double negative_ = kSaturationDepth * 0.35 * (1.0 - kEvenAsymmetry);
    double limit_ = kMonotonicMargin / (kSaturationDepth * 0.35 * (1.0 + kEvenAsymmetry));
    double previous_input_ = 0.0;
    bool adaa_ = true;
};
 
using TransformerBracket = TransformerBracketT<float>;
using TransformerBracket64 = TransformerBracketT<double>;
 
// ── The complete device ───────────────────────────────────────────────────

template <typename SampleType = float>
class DiodeBridgeCompressorT {
public:
    // ── Control ranges (mirrored by the catalog node's table) ─────────────
    static constexpr double kThresholdDbMin = -40.0;
    static constexpr double kThresholdDbMax = 6.0;
    static constexpr double kRatioMin = 1.5;
    static constexpr double kRatioMax = 20.0;
    static constexpr double kAttackMsMin = 0.5;
    static constexpr double kAttackMsMax = 100.0;
    static constexpr double kReleaseMsMin = 50.0;
    static constexpr double kReleaseMsMax = 2000.0;
    static constexpr double kKneeDbMin = 0.0;
    static constexpr double kKneeDbMax = 18.0;
    static constexpr double kMakeupDbMax = 24.0;
    static constexpr double kScHpfHzMin = 20.0;
    static constexpr double kScHpfHzMax = 400.0;
 
    // ── Design parameters ─────────────────────────────────────────────────

    static constexpr double kLimitRatio = 20.0;

    static constexpr double kCtrlSmoothMs = 1.0;

    static constexpr double kAutoSlowFactor = 6.0;

    static constexpr double kSustainTauMs = 200.0;

    static constexpr double kSustainThresholdDb = 3.0;

    static constexpr double kEnvelopeEpsilon = 1e-9;
 
    // ── Lifecycle ─────────────────────────────────────────────────────────

    void prepare(double sample_rate) {
        sample_rate_ = sample_rate > 0.0 ? sample_rate : sample_rate_;
        input_bracket_.prepare(sample_rate_);
        output_bracket_.prepare(sample_rate_);
        bridge_.prepare(sample_rate_);
        sidechain_hpf_.prepare(static_cast<SampleType>(sample_rate_));
        sidechain_hpf_.set_cutoff(static_cast<SampleType>(sc_hpf_hz_));
        fast_follower_.prepare(static_cast<SampleType>(sample_rate_));
        slow_follower_.prepare(static_cast<SampleType>(sample_rate_));
        fast_follower_.set_mode(BallisticsFilterT<SampleType>::Mode::peak);
        slow_follower_.set_mode(BallisticsFilterT<SampleType>::Mode::peak);
        update_ballistics();
        update_smoothing();
        reset();
    }
 
    void reset() {
        input_bracket_.reset();
        output_bracket_.reset();
        bridge_.reset();
        sidechain_dc_.reset();
        sidechain_hpf_.reset();
        fast_follower_.reset();
        slow_follower_.reset();
        smoothed_reduction_db_ = 0.0;
        sustain_ = 0.0;
        feedback_tap_ = 0.0;
    }
 
    // ── Controls (real units throughout) ──────────────────────────────────
 
    void set_threshold_db(double db) {
        threshold_db_ = std::clamp(dynamics::retain_finite(db, threshold_db_),
                                   kThresholdDbMin, kThresholdDbMax);
    }
 
    void set_ratio(double r) {
        ratio_ = std::clamp(dynamics::retain_finite(r, ratio_), kRatioMin, kRatioMax);
    }
 
    void set_knee_db(double w) {
        knee_db_ = std::clamp(dynamics::retain_finite(w, knee_db_),
                              kKneeDbMin, kKneeDbMax);
    }
 
    void set_attack_ms(double ms) {
        attack_ms_ = std::clamp(dynamics::retain_finite(ms, attack_ms_),
                                kAttackMsMin, kAttackMsMax);
        update_ballistics();
    }
 
    void set_release_ms(double ms) {
        release_ms_ = std::clamp(dynamics::retain_finite(ms, release_ms_),
                                 kReleaseMsMin, kReleaseMsMax);
        update_ballistics();
    }
 
    void set_makeup_db(double db) {
        makeup_db_ = std::clamp(dynamics::retain_finite(db, makeup_db_),
                                0.0, kMakeupDbMax);
        makeup_linear_ = units::db_to_linear(makeup_db_);
    }

    void set_character(double character01) {
        character_ = std::clamp(dynamics::retain_finite(character01, character_),
                                0.0, 1.0);
        bridge_.set_character(character_);
        input_bracket_.set_character(character_);
        output_bracket_.set_character(character_);
    }

    void set_mix_percent(double percent) {
        const double current_percent = mix_ * 100.0;
        mix_ = std::clamp(dynamics::retain_finite(percent, current_percent),
                          0.0, 100.0) / 100.0;
    }
 
    void set_sc_hpf_hz(double hz) {
        sc_hpf_hz_ = std::clamp(dynamics::retain_finite(hz, sc_hpf_hz_),
                                kScHpfHzMin, kScHpfHzMax);
        sidechain_hpf_.set_cutoff(static_cast<SampleType>(sc_hpf_hz_));
    }

    void set_auto_release(bool on) { auto_release_ = on; }

    void set_feedback(bool on) { feedback_ = on; }
 
    void set_adaa(bool on) {
        bridge_.set_adaa(on);
        input_bracket_.set_adaa(on);
        output_bracket_.set_adaa(on);
    }
 
    // ── Queries ───────────────────────────────────────────────────────────

    int latency_samples() const noexcept { return 0; }

    double gain_reduction_db() const { return smoothed_reduction_db_; }
    GainReduction gain_reduction() const noexcept {
        return GainReduction::from_signed_db(smoothed_reduction_db_);
    }

    double control_drive() const {
        return DiodeBridgeGainT<SampleType>::control_drive_for_gain_db(smoothed_reduction_db_);
    }

    double static_curve_db(double input_db) const {
        const double slope = ratio_ >= kLimitRatio ? 0.0 : 1.0 / ratio_;
        return dynamics::soft_knee_reduction_db(input_db - threshold_db_, knee_db_, slope);
    }

    double static_curve_feedback_db(double input_db) const {
        double gr = 0.0;
        for (int i = 0; i < 96; ++i) {
            const double next = 0.5 * (gr + static_curve_db(input_db + gr));
            if (next == gr) break;
            gr = next;
        }
        return gr;
    }

    static double worst_case_gain() {
        return units::db_to_linear(kMakeupDbMax) * TransformerBracketT<SampleType>::kPeakGain *
               TransformerBracketT<SampleType>::kPeakGain;
    }
 
    // ── Processing ────────────────────────────────────────────────────────
 
    SampleType process(SampleType input) {
        // Every stage below owns recursive state, including the feedback tap.
        // A non-finite input is therefore a recovery boundary: do not let it
        // enter the loop, and make the next finite sample start validly.
        if (!std::isfinite(static_cast<double>(input))) {
            reset();
            return SampleType{0};
        }
        const double dry = static_cast<double>(input);
 
        // Detection. The feedback tap is last sample's pre-makeup output: the
        // loop's only delay, and it lives in the CONTROL path, not the audio.
        const double detector_input = feedback_ ? feedback_tap_ : dry;
        const double envelope = detect(detector_input);
        const double envelope_db =
            units::linear_to_db(std::abs(envelope) + kEnvelopeEpsilon);
 
        const double target_db = static_curve_db(envelope_db);
        smoothed_reduction_db_ = snap_to_zero(smoothed_reduction_db_ +
                                              smooth_coef_ * (target_db - smoothed_reduction_db_));
 
        // Gain and drive are derived from ONE smoothed quantity, so they cannot
        // disagree about how much reduction is being applied.
        const double gain = units::db_to_linear(smoothed_reduction_db_);
        const double drive = 1.0 / gain - 1.0;
 
        const SampleType bracketed = input_bracket_.process(input);
        const SampleType attenuated = bridge_.process(bracketed, drive);
        const SampleType wet = output_bracket_.process(attenuated);
 
        feedback_tap_ = static_cast<double>(wet);
        return static_cast<SampleType>(dry * (1.0 - mix_) +
                                       static_cast<double>(wet) * makeup_linear_ * mix_);
    }
 
    void process_block(SampleType* io, int n) {
        for (int i = 0; i < n; ++i) io[i] = process(io[i]);
    }
 
private:
    double detect(double detector_input) {
        const SampleType blocked = sidechain_dc_.process(static_cast<SampleType>(detector_input));
        const SampleType filtered = sidechain_hpf_.process_highpass(blocked);
 
        const double fast = static_cast<double>(fast_follower_.process(filtered));
        // Keep the slow branch synchronized even while its readout is disabled.
        // Otherwise enabling Auto mid-release blends `fast` against a stale
        // zero state and briefly recovers gain faster than the manual follower.
        const double slow = static_cast<double>(slow_follower_.process(filtered));
        if (!auto_release_) return fast;
 
        // Program-dependent dual slope. The spec says `max(fast, slow)`, which
        // is identically the SLOW follower — during release both decay from the
        // same value and the slow one is always the larger, so the maximum
        // never selects the fast branch and the mode degenerates into a single
        // long release with nothing program-dependent about it. Blending on a
        // sustain indicator is what makes it dual-slope: fast recovery right
        // after a transient, slow hold once compression has been sustained.
        // Equals `max(fast, slow)` in the sustained limit.
        const double sustained = -smoothed_reduction_db_ > kSustainThresholdDb ? 1.0 : 0.0;
        sustain_ = snap_to_zero(sustain_ + sustain_coef_ * (sustained - sustain_));
        return fast + sustain_ * (slow - fast);
    }
 
    void update_ballistics() {
        const auto attack = static_cast<SampleType>(attack_ms_);
        fast_follower_.set_attack_ms(attack);
        slow_follower_.set_attack_ms(attack);
        fast_follower_.set_release_ms(static_cast<SampleType>(release_ms_));
        slow_follower_.set_release_ms(static_cast<SampleType>(release_ms_ * kAutoSlowFactor));
    }
 
    void update_smoothing() {
        smooth_coef_ = units::ms_to_onepole_coef(kCtrlSmoothMs, sample_rate_);
        sustain_coef_ = units::ms_to_onepole_coef(kSustainTauMs, sample_rate_);
    }
 
    double sample_rate_ = 48000.0;
 
    double threshold_db_ = -12.0;
    double ratio_ = 4.0;
    double knee_db_ = 6.0;
    double attack_ms_ = 3.0;
    double release_ms_ = 400.0;
    double makeup_db_ = 0.0;
    double character_ = 0.35;
    double mix_ = 1.0;
    double sc_hpf_hz_ = 100.0;
    bool auto_release_ = false;
    bool feedback_ = true;
 
    double makeup_linear_ = 1.0;
    double smooth_coef_ = 1.0;
    double sustain_coef_ = 1.0;
 
    double smoothed_reduction_db_ = 0.0;
    double sustain_ = 0.0;
    double feedback_tap_ = 0.0;
 
    TransformerBracketT<SampleType> input_bracket_{};
    TransformerBracketT<SampleType> output_bracket_{};
    DiodeBridgeGainT<SampleType> bridge_{};
    DcBlocker<SampleType> sidechain_dc_{};
    TptFilterT<SampleType> sidechain_hpf_{};
    BallisticsFilterT<SampleType> fast_follower_{};
    BallisticsFilterT<SampleType> slow_follower_{};
};
 
using DiodeBridgeCompressor = DiodeBridgeCompressorT<float>;
using DiodeBridgeCompressor64 = DiodeBridgeCompressorT<double>;
 
}  // namespace pulp::signal