speaker_cabinet.hpp

Go to the documentation of this file.

#pragma once

 
#include <pulp/signal/biquad.hpp>
#include <pulp/signal/dc_blocker.hpp>
#include <pulp/signal/denormal.hpp>
#include <pulp/signal/smoothed_value.hpp>
#include <pulp/signal/tpt_filter.hpp>
#include <pulp/signal/units.hpp>
 
#include <algorithm>
#include <cmath>
#include <complex>
 
namespace pulp::signal {

enum class SpeakerBoxType { sealed, open_back };

struct SpeakerDriverArchetype {
    const char* name;
    double fs_hz;           
    double qts;             
    double vas_litres;      
    double f_breakup_hz;    
    double baffle_width_m;  
};
 
namespace detail {

template <typename SampleType>
class DualTapResonatorT {
public:
    struct Taps {
        SampleType highpass;
        SampleType lowpass;
    };
 
    void set_coefficients(double cutoff_hz, double q, double sample_rate) {
        constexpr double pi = 3.14159265358979323846;
        const double nyquist_limit = 0.49 * sample_rate;
        const double fc = std::clamp(cutoff_hz, 1.0, nyquist_limit);
        const double safe_q = std::max(q, 1.0e-3);
        const double g = std::tan(pi * fc / sample_rate);
        const double k = 1.0 / safe_q;
        const double a1 = 1.0 / (1.0 + g * (g + k));
        g_ = static_cast<SampleType>(g);
        k_ = static_cast<SampleType>(k);
        a1_ = static_cast<SampleType>(a1);
        a2_ = static_cast<SampleType>(g * a1);
        a3_ = static_cast<SampleType>(g * g * a1);
    }
 
    Taps tick(SampleType input) {
        const SampleType v3 = input - ic2_;
        const SampleType v1 = a1_ * ic1_ + a2_ * v3;
        const SampleType v2 = ic2_ + a2_ * ic1_ + a3_ * v3;
        // Snap both integrator states: the TPT feedback path has no FTZ guard
        // of its own, exactly as SvfT documents.
        ic1_ = snap_to_zero(SampleType{2} * v1 - ic1_);
        ic2_ = snap_to_zero(SampleType{2} * v2 - ic2_);
        return {input - k_ * v1 - v2, v2};
    }
 
    void reset() {
        ic1_ = SampleType{0};
        ic2_ = SampleType{0};
    }
 
private:
    SampleType g_ = SampleType{0};
    SampleType k_ = SampleType{1};
    SampleType a1_ = SampleType{1};
    SampleType a2_ = SampleType{0};
    SampleType a3_ = SampleType{0};
    SampleType ic1_ = SampleType{0};
    SampleType ic2_ = SampleType{0};
};
 
}  // namespace detail

template <typename SampleType = float>
class SpeakerModelT {
public:
    // ── Driver archetypes ────────────────────────────────────────────────────
 
    static constexpr int kArchetypeCount = 5;

    static const SpeakerDriverArchetype& archetype(int index) {
        static constexpr SpeakerDriverArchetype kTable[kArchetypeCount] = {
            // name              fs     Qts   Vas(L)  f_breakup  baffle W(m)
            {"Brit-12 ceramic",  75.0,  0.50, 200.0,  1900.0,    0.50},
            {"Amer-12 ceramic",  95.0,  0.60, 150.0,  2200.0,    0.50},
            {"Alnico-12",        85.0,  0.45, 170.0,  2000.0,    0.48},
            {"Brit-10",         100.0,  0.55,  90.0,  2600.0,    0.42},
            {"Bass-15",          45.0,  0.35, 300.0,  1200.0,    0.60},
        };
        return kTable[std::clamp(index, 0, kArchetypeCount - 1)];
    }
 
    // ── Design parameters ────────────────────────────────────────────────────
    // Every value here is ours to tune, declared with default + range.

    static constexpr double kSpeedOfSoundMs = 343.0;

    static constexpr double kBoxVolumeLDefault = 28.0;
    static constexpr double kBoxVolumeLMin = 5.0;
    static constexpr double kBoxVolumeLMax = 150.0;

    static constexpr double kResonanceTrimSemitonesMin = -12.0;
    static constexpr double kResonanceTrimSemitonesMax = 12.0;

    static constexpr double kQResonanceMin = 0.4;
    static constexpr double kQResonanceMax = 2.0;

    static constexpr double kConeBreakupAmountDefault = 100.0;

    static constexpr double kTrebleRolloffHzDefault = 4000.0;
    static constexpr double kTrebleRolloffHzMin = 1500.0;
    static constexpr double kTrebleRolloffHzMax = 8000.0;

    static constexpr double kSemiInductanceDefault = 0.25;
    static constexpr double kSemiInductanceMin = 0.0;
    static constexpr double kSemiInductanceMax = 1.0;
    static constexpr double kSemiInductanceShelfDb = -12.0;
    static constexpr double kSemiInductanceShelfDbMin = -18.0;
    static constexpr double kSemiInductanceShelfDbMax = -6.0;

    static constexpr double kDriveDbDefault = 0.0;
    static constexpr double kDriveDbMin = -12.0;
    static constexpr double kDriveDbMax = 24.0;

    static constexpr double kCompressionAmountDefault = 60.0;

    static constexpr double kExcursionReference = 0.5;
    static constexpr double kExcursionReferenceMin = 0.1;
    static constexpr double kExcursionReferenceMax = 2.0;

    static constexpr double kBlFalloff = 1.5;
    static constexpr double kBlFalloffMin = 0.0;
    static constexpr double kBlFalloffMax = 4.0;

    static constexpr double kCmsStiffening = 0.35;
    static constexpr double kCmsStiffeningMin = 0.0;
    static constexpr double kCmsStiffeningMax = 1.0;

    static constexpr double kMicDistanceCmDefault = 3.0;
    static constexpr double kMicDistanceCmMin = 1.0;
    static constexpr double kMicDistanceCmMax = 60.0;

    static constexpr double kProximityGainK = 12.0;
    static constexpr double kProximityGainKMin = 0.0;
    static constexpr double kProximityGainKMax = 24.0;
    static constexpr double kProximityReferenceCm = 30.0;
    static constexpr double kProximityReferenceCmMin = 20.0;
    static constexpr double kProximityReferenceCmMax = 50.0;
    static constexpr double kProximityCeilingDb = 8.0;
    static constexpr double kProximityCeilingDbMin = 6.0;
    static constexpr double kProximityCeilingDbMax = 12.0;
    static constexpr double kProximityCornerHz = 250.0;
    static constexpr double kProximityCornerHzMin = 150.0;
    static constexpr double kProximityCornerHzMax = 400.0;

    static constexpr double kAirLossDbPerCm = 0.05;
    static constexpr double kAirLossDbPerCmMin = 0.02;
    static constexpr double kAirLossDbPerCmMax = 0.10;
    static constexpr double kAirLossCornerHz = 6000.0;
    static constexpr double kAirLossCornerHzMin = 4000.0;
    static constexpr double kAirLossCornerHzMax = 8000.0;

    static constexpr double kPresenceCapDb = 3.0;
    static constexpr double kPresenceEdgeDb = -8.0;
    static constexpr double kPresenceDbMax = 6.0;
    static constexpr double kPresenceDbMin = -12.0;
    static constexpr double kPresenceCornerHz = 2500.0;
    static constexpr double kPresenceCornerHzMin = 1500.0;
    static constexpr double kPresenceCornerHzMax = 3500.0;
    static constexpr double kMicPositionPctDefault = 30.0;

    static constexpr double kOffAxisOnAxisHz = 8000.0;
    static constexpr double kOffAxisOnAxisHzMin = 5000.0;
    static constexpr double kOffAxisOnAxisHzMax = 12000.0;
    static constexpr double kOffAxisFactor = 0.7;
    static constexpr double kOffAxisFactorMin = 0.0;
    static constexpr double kOffAxisFactorMax = 0.95;
    static constexpr double kMicAxisDegMin = 0.0;
    static constexpr double kMicAxisDegMax = 90.0;

    static constexpr double kBaffleStepDb = 6.0;

    static constexpr double kEdgeDistanceM = 0.15;
    static constexpr double kEdgeDistanceMMin = 0.08;
    static constexpr double kEdgeDistanceMMax = 0.30;
    static constexpr double kRipple1Q = 2.0;
    static constexpr double kRipple1Db = 2.5;
    static constexpr double kRipple2Q = 2.5;
    static constexpr double kRipple2Db = -2.0;
    static constexpr double kRippleQMin = 1.5;
    static constexpr double kRippleQMax = 4.0;
    static constexpr double kRippleDbMax = 4.0;
    static constexpr double kDiffractionAmountDefault = 50.0;

    static constexpr double kDipolePathM = 0.30;
    static constexpr double kDipolePathMMin = 0.15;
    static constexpr double kDipolePathMMax = 0.50;
    static constexpr double kDipoleNotchRatio = 1.8;
    static constexpr double kDipoleNotchRatioMin = 1.4;
    static constexpr double kDipoleNotchRatioMax = 2.2;
    static constexpr double kDipoleNotchQ = 1.5;
    static constexpr double kDipoleNotchQMin = 1.0;
    static constexpr double kDipoleNotchQMax = 3.0;
    static constexpr double kDipoleNotchDb = -4.0;
    static constexpr double kDipoleNotchDbMin = -8.0;
    static constexpr double kDipoleNotchDbMax = -2.0;

    static constexpr double kOutputTrimDbMin = -24.0;
    static constexpr double kOutputTrimDbMax = 24.0;

    static constexpr double kDcBlockerHz = 5.0;
    static constexpr double kDcBlockerHzMin = 2.0;
    static constexpr double kDcBlockerHzMax = 20.0;

    static constexpr double kSmoothingMs = 15.0;
    static constexpr double kSmoothingMsMin = 5.0;
    static constexpr double kSmoothingMsMax = 30.0;

    static constexpr int kControlBlock = 32;
    static constexpr int kControlBlockMin = 16;
    static constexpr int kControlBlockMax = 64;

    static constexpr double kWorstCaseGain = 10.0;
 
    // ── Cone-breakup mode shape ─────────────────────────────────────────────
 
    static constexpr int kBreakupModeCount = 4;

    struct BreakupMode {
        double ratio;
        double q;
        double gain_db;
    };

    static const BreakupMode& breakup_mode(int index) {
        static constexpr BreakupMode kShape[kBreakupModeCount] = {
            {1.00, 3.0, 4.0},
            {1.55, 4.0, 3.0},
            {2.20, 5.0, -4.0},
            {3.10, 6.0, 2.5},
        };
        return kShape[std::clamp(index, 0, kBreakupModeCount - 1)];
    }
    static constexpr double kBreakupRatioMin = 1.0;
    static constexpr double kBreakupRatioMax = 4.0;
    static constexpr double kBreakupQMin = 1.5;
    static constexpr double kBreakupQMax = 8.0;
    static constexpr double kBreakupGainDbMax = 6.0;
 
    // ── Lifecycle ────────────────────────────────────────────────────────────

    void prepare(double sample_rate) {
        if (std::isfinite(sample_rate) && sample_rate > 0.0) sample_rate_ = sample_rate;
        const auto smoothing = static_cast<SampleType>(kSmoothingMs / 1000.0);
        const auto fs = static_cast<SampleType>(sample_rate_);
        drive_lin_.set_ramp_time(smoothing, fs);
        out_trim_lin_.set_ramp_time(smoothing, fs);
        treble_hz_.set_ramp_time(smoothing, fs);
        q_smoothed_.set_ramp_time(smoothing, fs);
        presence_db_.set_ramp_time(smoothing, fs);
        proximity_db_.set_ramp_time(smoothing, fs);
        offaxis_hz_.set_ramp_time(smoothing, fs);
        inductance_lp_.prepare(fs);
        offaxis_lp_.prepare(fs);
        dipole_hp_.prepare(fs);
        // One-pole high-pass pole for the requested corner: p = 1 - 2*pi*fc/fs.
        constexpr double two_pi = 6.283185307179586476925286766559;
        dc_.set_pole(static_cast<SampleType>(
            std::clamp(1.0 - two_pi * kDcBlockerHz / sample_rate_, 0.0, 0.99999)));
        drive_lin_.set_immediate(static_cast<SampleType>(units::db_to_linear(drive_db_)));
        out_trim_lin_.set_immediate(static_cast<SampleType>(units::db_to_linear(out_trim_db_)));
        treble_hz_.set_immediate(static_cast<SampleType>(treble_hz_target_));
        q_smoothed_.set_immediate(static_cast<SampleType>(resonance_q()));
        presence_db_.set_immediate(static_cast<SampleType>(presence_shelf_db()));
        proximity_db_.set_immediate(static_cast<SampleType>(proximity_gain_db()));
        offaxis_hz_.set_immediate(static_cast<SampleType>(offaxis_corner_hz()));
        control_counter_ = 0;
        update_coefficients();
        reset();
    }
 
    // ── Parameters ───────────────────────────────────────────────────────────
 
    void set_driver_archetype(int index) {
        archetype_index_ = std::clamp(index, 0, kArchetypeCount - 1);
        retarget_smoothed();
    }
    void set_box_type(SpeakerBoxType type) {
        box_type_ = type;
        retarget_smoothed();
    }
    void set_box_volume_l(double litres) {
        if (!std::isfinite(litres)) return;
        box_volume_l_ = std::clamp(litres, kBoxVolumeLMin, kBoxVolumeLMax);
        retarget_smoothed();
    }
    void set_resonance_trim_semitones(double semitones) {
        if (!std::isfinite(semitones)) return;
        resonance_trim_ =
            std::clamp(semitones, kResonanceTrimSemitonesMin, kResonanceTrimSemitonesMax);
    }

    void set_q_resonance(double q) {
        if (!std::isfinite(q)) return;
        q_override_ = q > 0.0 ? std::clamp(q, kQResonanceMin, kQResonanceMax) : 0.0;
        retarget_smoothed();
    }
    void set_cone_breakup_amount(double percent) {
        if (!std::isfinite(percent)) return;
        breakup_amount_ = std::clamp(percent, 0.0, 100.0) / 100.0;
    }
    void set_treble_rolloff_hz(double hz) {
        if (!std::isfinite(hz)) return;
        treble_hz_target_ = std::clamp(hz, kTrebleRolloffHzMin, kTrebleRolloffHzMax);
        treble_hz_.set_target(static_cast<SampleType>(treble_hz_target_));
    }
    void set_drive_db(double db) {
        if (!std::isfinite(db)) return;
        drive_db_ = std::clamp(db, kDriveDbMin, kDriveDbMax);
        drive_lin_.set_target(static_cast<SampleType>(units::db_to_linear(drive_db_)));
    }
    void set_compression_amount(double percent) {
        if (!std::isfinite(percent)) return;
        compression_amount_ = std::clamp(percent, 0.0, 100.0) / 100.0;
    }
    void set_mic_distance_cm(double cm) {
        if (!std::isfinite(cm)) return;
        mic_distance_cm_ = std::clamp(cm, kMicDistanceCmMin, kMicDistanceCmMax);
        proximity_db_.set_target(static_cast<SampleType>(proximity_gain_db()));
    }
    void set_mic_position_pct(double percent) {
        if (!std::isfinite(percent)) return;
        mic_position_ = std::clamp(percent, 0.0, 100.0) / 100.0;
        presence_db_.set_target(static_cast<SampleType>(presence_shelf_db()));
    }
    void set_mic_axis_deg(double degrees) {
        if (!std::isfinite(degrees)) return;
        mic_axis_deg_ = std::clamp(degrees, kMicAxisDegMin, kMicAxisDegMax);
        offaxis_hz_.set_target(static_cast<SampleType>(offaxis_corner_hz()));
    }
    void set_diffraction_amount(double percent) {
        if (!std::isfinite(percent)) return;
        diffraction_amount_ = std::clamp(percent, 0.0, 100.0) / 100.0;
    }
    void set_output_trim_db(double db) {
        if (!std::isfinite(db)) return;
        out_trim_db_ = std::clamp(db, kOutputTrimDbMin, kOutputTrimDbMax);
        out_trim_lin_.set_target(static_cast<SampleType>(units::db_to_linear(out_trim_db_)));
    }
 
    // ── Audio ────────────────────────────────────────────────────────────────
 
    SampleType process(SampleType input) {
        if (!std::isfinite(static_cast<double>(input))) {
            reset();
            return SampleType{0};
        }
        if (control_counter_ == 0) update_coefficients();
        control_counter_ = (control_counter_ + 1) % kControlBlock;
 
        const SampleType drive = drive_lin_.next();
        const SampleType trim = out_trim_lin_.next();
 
        // Thiele-Small resonance. HP tap is the radiated pressure; LP tap is
        // proportional to cone displacement and feeds the nonlinear stage.
        const auto taps = resonance_.tick(dc_.process(input) * drive);
 
        // BL(x): even, strictly compressive, exactly unity at rest.
        const SampleType x_norm = taps.lowpass * excursion_scale_;
        const SampleType g_bl = SampleType{1} / (SampleType{1} + beta_ * x_norm * x_norm);
        excursion_ = x_norm;
        // Peak-track within the control block: the Cms stiffening is an
        // envelope effect, and sampling the instantaneous excursion once per
        // block would alias the carrier's own phase into the coefficient.
        const SampleType magnitude = std::abs(x_norm);
        if (magnitude > excursion_peak_) excursion_peak_ = magnitude;
 
        SampleType s = taps.highpass * g_bl;
        s = inductance(s);
        for (int i = 0; i < kBreakupModeCount; ++i) s = breakup_[i].process(s);
        s = baffle_.process(s);
        s = ripple_[0].process(s);
        s = ripple_[1].process(s);
        if (box_type_ == SpeakerBoxType::open_back) {
            s = dipole_hp_.process_highpass(s);
            s = dipole_notch_.process(s);
        }
        s = proximity_.process(s);
        s = air_loss_.process(s);
        s = presence_.process(s);
        s = offaxis_lp_.process_lowpass(s);
        const SampleType output = s * trim;
        if (!std::isfinite(static_cast<double>(output))) {
            reset();
            return SampleType{0};
        }
        return output;
    }
 
    void process(const SampleType* in, SampleType* out, int n) {
        for (int i = 0; i < n; ++i) out[i] = process(in[i]);
    }

    void reset() {
        dc_.reset();
        resonance_.reset();
        inductance_lp_.reset();
        semi_shelf_.reset();
        for (auto& mode : breakup_) mode.reset();
        baffle_.reset();
        ripple_[0].reset();
        ripple_[1].reset();
        dipole_hp_.reset();
        dipole_notch_.reset();
        proximity_.reset();
        air_loss_.reset();
        presence_.reset();
        offaxis_lp_.reset();
        excursion_ = SampleType{0};
        excursion_peak_ = SampleType{0};
        fc_dynamic_ = resonance_fc_hz();
        control_counter_ = 0;
    }

    int latency_samples() const { return 0; }

    double worst_case_gain() const { return kWorstCaseGain; }
 
    // ── Derived physics (tests compute their expected values from these) ─────

    double compliance_ratio() const {
        if (box_type_ == SpeakerBoxType::open_back) return 0.0;
        return archetype(archetype_index_).vas_litres / box_volume_l_;
    }

    double resonance_fc_hz() const {
        const double base = archetype(archetype_index_).fs_hz *
                            std::sqrt(1.0 + compliance_ratio());
        return base * std::pow(2.0, resonance_trim_ / 12.0);
    }

    double resonance_q() const {
        if (q_override_ > 0.0) return q_override_;
        return archetype(archetype_index_).qts * std::sqrt(1.0 + compliance_ratio());
    }

    double resonance_peak_db() const {
        const double q = resonance_q();
        const double inner = 1.0 - 1.0 / (4.0 * q * q);
        if (inner <= 0.0) return 0.0;
        return units::linear_to_db(q / std::sqrt(inner));
    }

    double resonance_peak_hz() const {
        const double q = resonance_q();
        const double inner = 1.0 - 1.0 / (2.0 * q * q);
        if (inner <= 0.0) return resonance_fc_hz();
        return resonance_fc_hz() / std::sqrt(inner);
    }

    double baffle_step_hz() const {
        return kSpeedOfSoundMs / (3.14159265358979323846 * archetype(archetype_index_).baffle_width_m);
    }

    double ripple_hz() const { return kSpeedOfSoundMs / (2.0 * kEdgeDistanceM); }

    double dipole_hz() const { return kSpeedOfSoundMs / (2.0 * kDipolePathM); }

    double breakup_mode_hz(int index) const {
        return archetype(archetype_index_).f_breakup_hz * breakup_mode(index).ratio;
    }

    double offaxis_corner_hz() const {
        constexpr double deg_to_rad = 3.14159265358979323846 / 180.0;
        const double reduction = kOffAxisFactor * std::sin(mic_axis_deg_ * deg_to_rad);
        return kOffAxisOnAxisHz * std::max(0.05, 1.0 - reduction);
    }

    double proximity_gain_db() const {
        const double raw = kProximityGainK * (1.0 / mic_distance_cm_ - 1.0 / kProximityReferenceCm);
        return std::clamp(raw, 0.0, kProximityCeilingDb);
    }

    double presence_shelf_db() const {
        return kPresenceCapDb + mic_position_ * (kPresenceEdgeDb - kPresenceCapDb);
    }

    double air_loss_db() const { return -kAirLossDbPerCm * mic_distance_cm_; }

    double inductance_magnitude_db(double hz) const {
        constexpr double pi = 3.14159265358979323846;
        const double w = 2.0 * pi * hz / sample_rate_;
        const std::complex<double> z_inv = std::polar(1.0, -w);
 
        // TPT one-pole lowpass: H(z) = g(1 + z^-1) / (1 - (1 - 2g) z^-1),
        // with the same bilinear g the filter itself computes from its corner.
        const double fc = static_cast<double>(inductance_lp_.cutoff());
        const double wa = 2.0 * sample_rate_ * std::tan(pi * fc / sample_rate_);
        const double g = wa / (2.0 * sample_rate_ + wa);
        const std::complex<double> lp =
            g * (1.0 + z_inv) / (1.0 - (1.0 - 2.0 * g) * z_inv);
 
        // The shelf's live coefficients, evaluated as H(e^jw).
        const auto c = semi_shelf_.coefficients();
        const std::complex<double> z2 = z_inv * z_inv;
        const std::complex<double> shelf =
            (static_cast<double>(c.b0) + static_cast<double>(c.b1) * z_inv +
             static_cast<double>(c.b2) * z2) /
            (1.0 + static_cast<double>(c.a1) * z_inv + static_cast<double>(c.a2) * z2);
 
        const double s = static_cast<double>(semi_);
        return units::linear_to_db(std::abs((1.0 - s) * lp + s * shelf));
    }

    double bl_beta() const { return compression_amount_ * kBlFalloff; }

    double cms_gamma() const { return compression_amount_ * kCmsStiffening; }

    double excursion() const { return static_cast<double>(excursion_); }

    double dynamic_fc_hz() const { return fc_dynamic_; }
 
    int archetype_index() const { return archetype_index_; }
    SpeakerBoxType box_type() const { return box_type_; }
    double sample_rate() const { return sample_rate_; }
 
private:
    // ── Voice-coil inductance: first-order rolloff blended with a shelf ──────
    // Leach's semi-inductive motor measures a gentler slope than a clean
    // -6 dB/oct, because eddy-current losses make the coil's impedance rise
    // sub-linearly. Blending a fraction of a high-shelved copy of the input
    // back in floors the decay instead of letting it fall forever.
    SampleType inductance(SampleType x) {
        const SampleType rolled = inductance_lp_.process_lowpass(x);
        const SampleType shelved = semi_shelf_.process(x);
        return (SampleType{1} - semi_) * rolled + semi_ * shelved;
    }
 
    void retarget_smoothed() {
        q_smoothed_.set_target(static_cast<SampleType>(resonance_q()));
        presence_db_.set_target(static_cast<SampleType>(presence_shelf_db()));
        proximity_db_.set_target(static_cast<SampleType>(proximity_gain_db()));
        offaxis_hz_.set_target(static_cast<SampleType>(offaxis_corner_hz()));
    }

    void update_coefficients() {
        const auto fs = static_cast<SampleType>(sample_rate_);
 
        // Advance the smoothers by one control block so their ramps track
        // real time even though coefficients only move at this cadence.
        const double q = static_cast<double>(q_smoothed_.current());
        const double treble = static_cast<double>(treble_hz_.current());
        const double presence = static_cast<double>(presence_db_.current());
        const double proximity = static_cast<double>(proximity_db_.current());
        const double offaxis = static_cast<double>(offaxis_hz_.current());
        q_smoothed_.skip(kControlBlock);
        treble_hz_.skip(kControlBlock);
        presence_db_.skip(kControlBlock);
        proximity_db_.skip(kControlBlock);
        offaxis_hz_.skip(kControlBlock);
 
        // Cms(x) stiffening from the PREVIOUS block's peak excursion: raising
        // fc at high excursion reduces the displacement that produced it, so
        // the loop is compressive negative feedback and cannot run away.
        const double peak = static_cast<double>(excursion_peak_);
        fc_dynamic_ = resonance_fc_hz() * (1.0 + cms_gamma() * peak * peak);
        excursion_peak_ = SampleType{0};
        resonance_.set_coefficients(fc_dynamic_, q, sample_rate_);
 
        beta_ = static_cast<SampleType>(bl_beta());
        excursion_scale_ = static_cast<SampleType>(1.0 / kExcursionReference);
 
        inductance_lp_.set_cutoff(static_cast<SampleType>(treble));
        semi_ = static_cast<SampleType>(semi_inductance_);
        semi_shelf_.set_coefficients(Biquad::Type::high_shelf, static_cast<SampleType>(treble),
                                     kShelfQ, fs,
                                     static_cast<SampleType>(kSemiInductanceShelfDb));
 
        for (int i = 0; i < kBreakupModeCount; ++i) {
            const auto& mode = breakup_mode(i);
            breakup_[i].set_coefficients(
                Biquad::Type::peaking, static_cast<SampleType>(breakup_mode_hz(i)),
                static_cast<SampleType>(mode.q), fs,
                static_cast<SampleType>(mode.gain_db * breakup_amount_));
        }
 
        baffle_.set_coefficients(Biquad::Type::high_shelf,
                                 static_cast<SampleType>(baffle_step_hz()), kShelfQ, fs,
                                 static_cast<SampleType>(kBaffleStepDb * diffraction_amount_));
        const double ripple1 = ripple_hz();
        ripple_[0].set_coefficients(Biquad::Type::peaking, static_cast<SampleType>(ripple1),
                                    static_cast<SampleType>(kRipple1Q), fs,
                                    static_cast<SampleType>(kRipple1Db * diffraction_amount_));
        ripple_[1].set_coefficients(Biquad::Type::peaking, static_cast<SampleType>(2.0 * ripple1),
                                    static_cast<SampleType>(kRipple2Q), fs,
                                    static_cast<SampleType>(kRipple2Db * diffraction_amount_));
 
        dipole_hp_.set_cutoff(static_cast<SampleType>(dipole_hz()));
        dipole_notch_.set_coefficients(
            Biquad::Type::peaking, static_cast<SampleType>(dipole_hz() * kDipoleNotchRatio),
            static_cast<SampleType>(kDipoleNotchQ), fs,
            static_cast<SampleType>(kDipoleNotchDb));
 
        proximity_.set_coefficients(Biquad::Type::low_shelf,
                                    static_cast<SampleType>(kProximityCornerHz), kShelfQ, fs,
                                    static_cast<SampleType>(proximity));
        air_loss_.set_coefficients(Biquad::Type::high_shelf,
                                   static_cast<SampleType>(kAirLossCornerHz), kShelfQ, fs,
                                   static_cast<SampleType>(air_loss_db()));
        presence_.set_coefficients(Biquad::Type::high_shelf,
                                   static_cast<SampleType>(kPresenceCornerHz), kShelfQ, fs,
                                   static_cast<SampleType>(presence));
        offaxis_lp_.set_cutoff(static_cast<SampleType>(offaxis));
    }

    static constexpr SampleType kShelfQ = static_cast<SampleType>(0.70710678118654752440);
 
    using Biquad = BiquadT<SampleType>;
 
    // ── Configuration ────────────────────────────────────────────────────────
    double sample_rate_ = 48000.0;
    int archetype_index_ = 0;
    SpeakerBoxType box_type_ = SpeakerBoxType::open_back;
    double box_volume_l_ = kBoxVolumeLDefault;
    double resonance_trim_ = 0.0;
    double q_override_ = 0.0;
    double breakup_amount_ = kConeBreakupAmountDefault / 100.0;
    double treble_hz_target_ = kTrebleRolloffHzDefault;
    double semi_inductance_ = kSemiInductanceDefault;
    double drive_db_ = kDriveDbDefault;
    double compression_amount_ = kCompressionAmountDefault / 100.0;
    double mic_distance_cm_ = kMicDistanceCmDefault;
    double mic_position_ = kMicPositionPctDefault / 100.0;
    double mic_axis_deg_ = 0.0;
    double diffraction_amount_ = kDiffractionAmountDefault / 100.0;
    double out_trim_db_ = 0.0;
 
    // ── Smoothed control values ──────────────────────────────────────────────
    SmoothedValue<SampleType> drive_lin_{SampleType{1}};
    SmoothedValue<SampleType> out_trim_lin_{SampleType{1}};
    SmoothedValue<SampleType> treble_hz_{static_cast<SampleType>(kTrebleRolloffHzDefault)};
    SmoothedValue<SampleType> q_smoothed_{SampleType{1}};
    SmoothedValue<SampleType> presence_db_{SampleType{0}};
    SmoothedValue<SampleType> proximity_db_{SampleType{0}};
    SmoothedValue<SampleType> offaxis_hz_{static_cast<SampleType>(kOffAxisOnAxisHz)};
 
    // ── Audio state ──────────────────────────────────────────────────────────
    DcBlocker<SampleType> dc_;
    detail::DualTapResonatorT<SampleType> resonance_;
    TptFilterT<SampleType> inductance_lp_;
    Biquad semi_shelf_;
    Biquad breakup_[kBreakupModeCount];
    Biquad baffle_;
    Biquad ripple_[2];
    TptFilterT<SampleType> dipole_hp_;
    Biquad dipole_notch_;
    Biquad proximity_;
    Biquad air_loss_;
    Biquad presence_;
    TptFilterT<SampleType> offaxis_lp_;
 
    SampleType beta_ = SampleType{0};
    SampleType semi_ = static_cast<SampleType>(kSemiInductanceDefault);
    SampleType excursion_scale_ = static_cast<SampleType>(1.0 / kExcursionReference);
    SampleType excursion_ = SampleType{0};
    SampleType excursion_peak_ = SampleType{0};
    double fc_dynamic_ = 0.0;
    int control_counter_ = 0;
};
 
using SpeakerModel = SpeakerModelT<float>;
using SpeakerModel64 = SpeakerModelT<double>;

 
}  // namespace pulp::signal