fet_compressor.hpp

Go to the documentation of this file.

#pragma once

 
#include <pulp/signal/biquad.hpp>
#include <pulp/signal/denormal.hpp>
#include <pulp/signal/dynamics_core.hpp>
#include <pulp/signal/dynamics_contract.hpp>
#include <pulp/signal/units.hpp>
#include <pulp/signal/windowed_sinc_design.hpp>
 
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <cstdint>
 
namespace pulp::signal {

enum class FetRatio : std::uint8_t {
    r4_1,
    r8_1,
    r12_1,
    r20_1,
    all_buttons_in,
};

template <typename SampleType = float>
class FetCompressorT {
public:
    using Ratio = FetRatio;
 
    // ── Resampler geometry: the latency contract's inputs ──────────────────

    static constexpr int kOversamplingFactor = 4;

    static constexpr int kResamplerTaps = 65;

    static constexpr double kResamplerKaiserBeta = 8.0;

    static constexpr double kResamplerCutoffFraction = 0.96;

    static constexpr int kLatencySamples =
        2 * ((kResamplerTaps - 1) / 2) / kOversamplingFactor;
 
    // ── Design parameters: the divider ────────────────────────────────────

    static constexpr double kFetOnResistance = 120.0;

    static constexpr double kFetOffResistance = 1.0e6;

    static constexpr double kGainReductionMaxDb = 40.0;
 
    // ── Design parameters: the coloration term ────────────────────────────

    static constexpr double kColorationAlphaMax = 0.08;

    static constexpr double kColorationAlphaCeiling = 0.16;

    static constexpr double kColorationHeadroomClamp = 1.2;
 
    // ── Design parameters: all-buttons-in ─────────────────────────────────
    //
    // **Honest gap, stated plainly.** Moore (2012) documents the QUALITATIVE
    // phenomenon — heavier, more distorted, less predictable compression — from
    // listening and usage analysis, not circuit-level component values. No
    // citable source gives numbers for any of the five constants below. They
    // are original engineering choices; re-tuning by ear must stay inside the
    // stated ranges and must preserve the DIRECTION of each deviation.

    static constexpr double kAbiRatio = 20.0;

    static constexpr double kAbiKneeWidenMult = 2.5;

    static constexpr double kAbiBiasShiftDb = 1.5;

    static constexpr double kAbiAttackScale = 0.6;

    static constexpr double kAbiReleaseScale = 0.75;

    static constexpr double kAbiAlphaExtraMult = 1.4;
 
    // ── Design parameters: the transformer tilt ───────────────────────────

    static constexpr double kTransformerCornerHz = 60.0;

    static constexpr double kTransformerTiltDbMax = 3.0;

    static constexpr double kTransformerQ = 0.70710678118654752;
 
    // ── Design parameters: everything else ────────────────────────────────

    static constexpr double kThresholdDbfs = -18.0;

    static constexpr double kLevelEpsilon = dynamics::kLevelEpsilon;
 
    // Control ranges. These mirror the catalog node's parameter table; the
    // table is the canonical declaration and these are the same numbers at the
    // call site, not a second independent one.
    static constexpr double kInputGainDbMin = -20.0;
    static constexpr double kInputGainDbMax = 40.0;
    static constexpr double kOutputGainDbMin = -20.0;
    static constexpr double kOutputGainDbMax = 20.0;
    static constexpr double kAttackUsMin = 20.0;
    static constexpr double kAttackUsMax = 800.0;
    static constexpr double kReleaseMsMin = 50.0;
    static constexpr double kReleaseMsMax = 1100.0;
    static constexpr double kKneeDbMin = 0.0;
    static constexpr double kKneeDbMax = 6.0;
 
    // ── Lifecycle ─────────────────────────────────────────────────────────
 
    FetCompressorT() { prepare(44100.0); }

    void prepare(double sample_rate) {
        sample_rate_ = sample_rate > 0.0 ? sample_rate : sample_rate_;
        oversampled_rate_ = sample_rate_ * kOversamplingFactor;
 
        design_resampler();
        solve_divider();
        update_ratio_dependents();
        update_coefficients();
        update_transformer();
        reset();
    }

    void reset() {
        base_history_.fill(0.0);
        os_history_.fill(0.0);
        dry_delay_.fill(0.0);
        base_write_ = 0;
        os_write_ = 0;
        dry_write_ = 0;
        detector_db_ = 0.0;
        control_ = 0.0;
        feedback_ = 0.0;
        transformer_.reset();
    }
 
    // ── Controls (real units throughout) ──────────────────────────────────
    //
    // Taken as `double` rather than `SampleType`. A `double` parameter accepts
    // everything a `SampleType` one would and rounds nothing on the way in, and
    // it matches `feedforward_compressor.hpp` so the two members of the family
    // present the same call shape.

    void set_input_gain_db(double db) {
        input_gain_db_ = std::clamp(dynamics::retain_finite(db, input_gain_db_),
                                    kInputGainDbMin, kInputGainDbMax);
        input_gain_ = units::db_to_linear(input_gain_db_);
    }

    void set_output_gain_db(double db) {
        output_gain_db_ = std::clamp(dynamics::retain_finite(db, output_gain_db_),
                                     kOutputGainDbMin, kOutputGainDbMax);
        output_gain_ = units::db_to_linear(output_gain_db_);
    }
 
    void set_ratio(Ratio r) {
        ratio_button_ = r;
        update_ratio_dependents();
        update_coefficients();
    }
 
    Ratio ratio() const { return ratio_button_; }
 
    void set_attack_us(double us) {
        attack_us_ = std::clamp(dynamics::retain_finite(us, attack_us_),
                                kAttackUsMin, kAttackUsMax);
        update_coefficients();
    }
 
    void set_release_ms(double ms) {
        release_ms_ = std::clamp(dynamics::retain_finite(ms, release_ms_),
                                 kReleaseMsMin, kReleaseMsMax);
        update_coefficients();
    }

    void set_knee_db(double db) {
        knee_request_db_ = std::clamp(dynamics::retain_finite(db, knee_request_db_),
                                      kKneeDbMin, kKneeDbMax);
        update_ratio_dependents();
    }

    void set_transformer_amount(double pct01) {
        transformer_amount_ = std::clamp(
            dynamics::retain_finite(pct01, transformer_amount_), 0.0, 1.0);
        update_transformer();
    }

    void set_mix(double pct01) {
        mix_ = std::clamp(dynamics::retain_finite(pct01, mix_), 0.0, 1.0);
    }
 
    // ── Reported quantities ───────────────────────────────────────────────

    int latency_samples() const noexcept { return kLatencySamples; }
 
    double sample_rate() const noexcept { return sample_rate_; }
    double oversampled_rate() const noexcept { return oversampled_rate_; }

    double static_curve_db(double level_db) const {
        return dynamics::soft_knee_output_db(level_db, kThresholdDbfs, knee_db_, ratio_value_);
    }

    double gain_computer_db(double level_db) const {
        return level_db - static_curve_db(level_db);
    }

    double measured_static_curve_db(double input_db) const {
        const double m = loop_slope();
        const double w = knee_db_;
        const double over = input_db - bias_shift_db_ - kThresholdDbfs;
        if (2.0 * over < -w) return kThresholdDbfs + over;
        if (w > 0.0 && m > 0.0 && 2.0 * over <= w * (1.0 + m)) {
            // Quadratic branch: (m/2W)·s² + s − (over + W/2) = 0, s = u + W/2.
            const double s =
                (w / m) * (std::sqrt(1.0 + (2.0 * m / w) * (over + 0.5 * w)) - 1.0);
            return kThresholdDbfs + s - 0.5 * w;
        }
        return kThresholdDbfs + over / (1.0 + m);
    }

    double measured_gain_reduction_db(double input_db) const {
        return input_db - measured_static_curve_db(input_db);
    }

    double loop_slope() const { return 1.0 - 1.0 / ratio_value_; }

    double measured_ratio() const { return 1.0 + loop_slope(); }

    double measured_knee_db() const { return knee_db_ * (2.0 + loop_slope()) * 0.5; }

    double nominal_ratio() const { return ratio_value_; }

    double effective_knee_db() const { return knee_db_; }

    double bias_shift_db() const { return bias_shift_db_; }

    double coloration_depth() const { return alpha_max_; }

    double attack_coefficient() const { return attack_coef_; }
    double release_coefficient() const { return release_coef_; }

    double gain_reduction_db() const { return detector_db_; }
    GainReduction gain_reduction() const noexcept {
        return GainReduction::from_magnitude_db(detector_db_);
    }

    double control_voltage() const { return control_; }
 
    // ── The divider law ───────────────────────────────────────────────────

    double divider_conductance(double c) const {
        return g_min_ + std::clamp(c, 0.0, 1.0) * (g_max_ - g_min_);
    }

    double divider_small_signal_gain(double c) const {
        return trim_ / (1.0 + series_resistance_ * divider_conductance(c));
    }

    double divider_gain(double c, double v) const {
        return trim_ / (1.0 + series_resistance_ *
                                  divider_conductance(c) * coloration_multiplier(c, v));
    }

    double coloration_multiplier(double c, double v) const {
        const double clamped = std::clamp(v, -kColorationHeadroomClamp, kColorationHeadroomClamp);
        return 1.0 + alpha_max_ * std::clamp(c, 0.0, 1.0) * clamped;
    }

    double coloration_multiplier_bound() const {
        return 1.0 + alpha_max_ * kColorationHeadroomClamp;
    }

    double control_for_reduction_db(double db) const {
        if (!(db > 0.0)) return 0.0;
        const double rest = 1.0 + series_resistance_ * g_min_;
        const double target = rest * std::pow(10.0, db / 20.0);
        const double c = (target - rest) / (series_resistance_ * (g_max_ - g_min_));
        return std::clamp(c, 0.0, 1.0);
    }

    static constexpr double kDividerGainSupremum = 1.0;
 
    bool divider_supremum_is_provable() const {
        const double a = kColorationAlphaCeiling * kColorationHeadroomClamp;
        return (g_max_ - g_min_) * (1.0 - a) > a * g_min_;
    }

    double resampler_peak_gain_bound() const { return resampler_peak_bound_; }

    double worst_case_gain() const {
        const double wet = input_gain_ * kDividerGainSupremum * resampler_peak_bound_ *
                           output_gain_ * kTransformerGainSupremum;
        return (1.0 - mix_) + mix_ * wet;
    }

    static constexpr double kTransformerGainSupremum = 1.0;
 
    // ── Processing ────────────────────────────────────────────────────────

    SampleType process(SampleType input) {
        // The feedback tap, both resampler histories, the dry delay, and the
        // transformer are recursive. Reject a non-finite sample before any of
        // them sees it, then restart from the class's documented fresh state.
        if (!std::isfinite(static_cast<double>(input))) {
            reset();
            return SampleType{0};
        }
        const double x = static_cast<double>(input);
        const double dry = push_dry(x);
        const double driven = x * input_gain_;
 
        // Interpolate: one base sample in, `kOversamplingFactor` out.
        base_history_[base_write_] = driven;
 
        double decimated = 0.0;
        for (int phase = 0; phase < kOversamplingFactor; ++phase) {
            const double v = interpolate_phase(phase);
 
            // ── the loop ──────────────────────────────────────────────────
            // The detector reads the PREVIOUS oversampled divider output. That
            // one line is the entire topological difference from the
            // feedforward design.
            const double level_db = dynamics::amplitude_db(feedback_, kLevelEpsilon);
            const double requested = gain_computer_db(level_db) + bias_shift_db_;
            detector_db_ = snap_to_zero(
                requested > detector_db_
                    ? attack_coef_ * detector_db_ + (1.0 - attack_coef_) * requested
                    : release_coef_ * detector_db_ + (1.0 - release_coef_) * requested);
            control_ = control_for_reduction_db(detector_db_);
 
            const double y = v * divider_gain(control_, v);
            feedback_ = y;
 
            os_history_[os_write_] = y;
            // The first oversample is the decimator's output grid point; taking
            // it here keeps the composite delay an exact integer.
            if (phase == 0) decimated = decimate();
            os_write_ = os_write_ + 1 < kResamplerTaps ? os_write_ + 1 : 0;
        }
        base_write_ = base_write_ + 1 < kBaseHistory ? base_write_ + 1 : 0;
 
        const double wet = transformer_.process(decimated * output_gain_);
        return static_cast<SampleType>(snap_to_zero(dry + mix_ * (wet - dry)));
    }
 
    void process_block(SampleType* io, int n) {
        for (int i = 0; i < n; ++i) io[i] = process(io[i]);
    }
 
private:
    static constexpr int kBaseHistory =
        (kResamplerTaps + kOversamplingFactor - 1) / kOversamplingFactor;
 
    // ── Resampler ─────────────────────────────────────────────────────────
 
    void design_resampler() {
        // Cutoff expressed in cycles per OVERSAMPLED sample: base Nyquist is
        // `0.5/kOversamplingFactor` there.
        const double cutoff =
            kResamplerCutoffFraction * 0.5 / static_cast<double>(kOversamplingFactor);
        const auto prototype = design_windowed_sinc(static_cast<std::size_t>(kResamplerTaps),
                                                    cutoff, kResamplerKaiserBeta);
        double decimation_l1 = 0.0;
        for (int i = 0; i < kResamplerTaps; ++i) {
            decimation_[static_cast<std::size_t>(i)] = prototype[static_cast<std::size_t>(i)];
            decimation_l1 += std::abs(prototype[static_cast<std::size_t>(i)]);
        }
 
        // Polyphase decomposition of the interpolator. Zero-stuffing by L drops
        // the DC gain to 1/L, so each branch carries the compensating factor L.
        double interpolation_l1 = 0.0;
        for (int phase = 0; phase < kOversamplingFactor; ++phase) {
            int count = 0;
            double branch_l1 = 0.0;
            for (int tap = phase; tap < kResamplerTaps; tap += kOversamplingFactor) {
                const double coefficient =
                    kOversamplingFactor * prototype[static_cast<std::size_t>(tap)];
                interpolation_[static_cast<std::size_t>(phase)][static_cast<std::size_t>(count)] =
                    coefficient;
                branch_l1 += std::abs(coefficient);
                ++count;
            }
            interpolation_length_[static_cast<std::size_t>(phase)] = count;
            interpolation_l1 = std::max(interpolation_l1, branch_l1);
        }
        resampler_peak_bound_ = interpolation_l1 * decimation_l1;
    }

    double interpolate_phase(int phase) const {
        const auto p = static_cast<std::size_t>(phase);
        const int count = interpolation_length_[p];
        double sum = 0.0;
        int index = base_write_;
        for (int k = 0; k < count; ++k) {
            sum += interpolation_[p][static_cast<std::size_t>(k)] *
                   base_history_[static_cast<std::size_t>(index)];
            if (--index < 0) index = kBaseHistory - 1;
        }
        return sum;
    }

    double decimate() const {
        double sum = 0.0;
        int index = os_write_;
        for (int k = 0; k < kResamplerTaps; ++k) {
            sum += decimation_[static_cast<std::size_t>(k)] *
                   os_history_[static_cast<std::size_t>(index)];
            if (--index < 0) index = kResamplerTaps - 1;
        }
        return sum;
    }

    double push_dry(double x) {
        constexpr std::size_t size = kLatencySamples + 1;
        const std::size_t read = (dry_write_ + size - kLatencySamples) % size;
        const double out = dry_delay_[read];
        dry_delay_[dry_write_] = x;
        dry_write_ = (dry_write_ + 1) % size;
        return out;
    }
 
    // ── Derived quantities ────────────────────────────────────────────────

    void solve_divider() {
        g_min_ = 1.0 / kFetOffResistance;
        g_max_ = 1.0 / kFetOnResistance;
        const double k = std::pow(10.0, -kGainReductionMaxDb / 20.0);
        series_resistance_ = kFetOffResistance * kFetOnResistance * (k - 1.0) /
                             (kFetOnResistance - k * kFetOffResistance);
        trim_ = 1.0 + series_resistance_ * g_min_;
    }
 
    void update_ratio_dependents() {
        const bool abi = ratio_button_ == Ratio::all_buttons_in;
        switch (ratio_button_) {
            case Ratio::r4_1: ratio_value_ = 4.0; break;
            case Ratio::r8_1: ratio_value_ = 8.0; break;
            case Ratio::r12_1: ratio_value_ = 12.0; break;
            case Ratio::r20_1: ratio_value_ = 20.0; break;
            case Ratio::all_buttons_in: ratio_value_ = kAbiRatio; break;
        }
        knee_db_ = abi ? knee_request_db_ * kAbiKneeWidenMult : knee_request_db_;
        bias_shift_db_ = abi ? kAbiBiasShiftDb : 0.0;
        alpha_max_ = abi ? std::min(kColorationAlphaMax * kAbiAlphaExtraMult,
                                    kColorationAlphaCeiling)
                         : kColorationAlphaMax;
    }
 
    void update_coefficients() {
        const bool abi = ratio_button_ == Ratio::all_buttons_in;
        const double attack_s =
            attack_us_ * 1e-6 * (abi ? kAbiAttackScale : 1.0);
        const double release_s =
            release_ms_ * 1e-3 * (abi ? kAbiReleaseScale : 1.0);
        attack_coef_ = retain(attack_s);
        release_coef_ = retain(release_s);
    }

    double retain(double tau_seconds) const {
        return dynamics::one_pole_retain(tau_seconds, oversampled_rate_);
    }
 
    void update_transformer() {
        transformer_.set_coefficients(
            BiquadT<double>::Type::low_shelf, static_cast<double>(kTransformerCornerHz),
            kTransformerQ, sample_rate_, -transformer_amount_ * kTransformerTiltDbMax);
    }
 
    // ── State ─────────────────────────────────────────────────────────────
 
    double sample_rate_ = 44100.0;
    double oversampled_rate_ = 44100.0 * kOversamplingFactor;
 
    // Controls.
    double input_gain_db_ = 0.0;
    double output_gain_db_ = 0.0;
    double attack_us_ = 200.0;
    double release_ms_ = 300.0;
    double knee_request_db_ = 1.0;
    double transformer_amount_ = 0.6;
    double mix_ = 1.0;
    Ratio ratio_button_ = Ratio::r4_1;
 
    // Derived from controls.
    double input_gain_ = 1.0;
    double output_gain_ = 1.0;
    double ratio_value_ = 4.0;
    double knee_db_ = 1.0;
    double bias_shift_db_ = 0.0;
    double alpha_max_ = kColorationAlphaMax;
    double attack_coef_ = 0.0;
    double release_coef_ = 0.0;
 
    // Derived from the divider's design parameters.
    double g_min_ = 0.0;
    double g_max_ = 0.0;
    double series_resistance_ = 0.0;
    double trim_ = 1.0;
 
    // Resampler.
    std::array<double, kResamplerTaps> decimation_{};
    std::array<std::array<double, kBaseHistory>, kOversamplingFactor> interpolation_{};
    std::array<int, kOversamplingFactor> interpolation_length_{};
    double resampler_peak_bound_ = 1.0;
 
    // Running state.
    std::array<double, kBaseHistory> base_history_{};
    std::array<double, kResamplerTaps> os_history_{};
    std::array<double, kLatencySamples + 1> dry_delay_{};
    int base_write_ = 0;
    int os_write_ = 0;
    std::size_t dry_write_ = 0;
 
    double detector_db_ = 0.0;
    double control_ = 0.0;
    double feedback_ = 0.0;
    BiquadT<double> transformer_{};
};
 
using FetCompressor = FetCompressorT<float>;
using FetCompressor64 = FetCompressorT<double>;
 
}  // namespace pulp::signal