additive_bank.hpp

Go to the documentation of this file.

#pragma once

 
#include <pulp/signal/denormal.hpp>
#include <pulp/signal/additive_spectral_envelope.hpp>
#include <pulp/signal/envelope.hpp>
#include <pulp/signal/fast_math.hpp>
#include <pulp/signal/osc/phase.hpp>
#include <pulp/signal/rng.hpp>
#include <pulp/signal/units.hpp>
 
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <type_traits>
#include <vector>
 
namespace pulp::signal {

struct VoicePartial {
    double ratio = 1.0;      
    double amp = 1.0;        
    double phase01 = 0.0;    
    double decay_ms = 0.0;   
};

struct VoiceTable {
    static constexpr int kMaxPartials = 256;

    bool harmonic = true;
    int count = 0;
    std::array<VoicePartial, kMaxPartials> partials{};
 
    void clear() { count = 0; }
    bool add(const VoicePartial& p) {
        if (count >= kMaxPartials) return false;
        partials[static_cast<std::size_t>(count++)] = p;
        return true;
    }
};

inline VoiceTable make_organ_voice(int count = 64) {
    static constexpr double kDrawbarRatio[9] = {0.5, 1.0, 1.5, 2.0, 3.0,
                                                4.0, 5.0, 6.0, 8.0};
    static constexpr double kDrawbarAmp[9] = {0.50, 1.00, 0.35, 0.70, 0.45,
                                              0.55, 0.25, 0.20, 0.30};
    constexpr double kTailLevel = 0.30;
 
    VoiceTable v;
    v.harmonic = true;
    const int n = std::clamp(count, 1, VoiceTable::kMaxPartials);
    for (int p = 0; p < n; ++p) {
        VoicePartial row;
        if (p < 9) {
            row.ratio = kDrawbarRatio[p];
            row.amp = kDrawbarAmp[p];
        } else {
            // −12 dB/oct is amplitude ∝ 1/ratio², which is what the exponent 2
            // here means; writing it as a power of the ratio keeps the tail
            // scale-invariant rather than tied to a partial index.
            row.ratio = static_cast<double>(p + 1);
            row.amp = kTailLevel / (row.ratio * row.ratio);
        }
        row.phase01 = 0.0;
        row.decay_ms = 0.0;   // sustained
        v.add(row);
    }
    return v;
}

inline VoiceTable make_bell_voice(int count = 64) {
    static constexpr double kModeRatio[10] = {0.50, 1.00, 1.19, 1.50, 2.00,
                                              2.50, 2.66, 3.00, 4.00, 5.00};
    static constexpr double kModeAmp[10] = {0.60, 1.00, 0.70, 0.55, 0.80,
                                            0.35, 0.30, 0.40, 0.30, 0.20};
    static constexpr double kModeDecayMs[10] = {8000, 6000, 4000, 3000, 2500,
                                                1500, 1400, 1200, 800,  500};
    constexpr double kClusterRatioStep = 1.14;
    constexpr double kClusterLevel = 0.20;
    constexpr double kClusterAmpExponent = 1.5;
    constexpr double kClusterDecayScaleMs = 3000.0;
 
    VoiceTable v;
    v.harmonic = false;   // modal: the ratios are the measurement
    const int n = std::clamp(count, 1, VoiceTable::kMaxPartials);
    for (int p = 0; p < n; ++p) {
        VoicePartial row;
        if (p < 10) {
            row.ratio = kModeRatio[p];
            row.amp = kModeAmp[p];
            row.decay_ms = kModeDecayMs[p];
        } else {
            row.ratio = 5.0 * std::pow(kClusterRatioStep,
                                       static_cast<double>(p - 9));
            row.amp = kClusterLevel * std::pow(1.0 / row.ratio,
                                               kClusterAmpExponent);
            row.decay_ms = kClusterDecayScaleMs / row.ratio;
        }
        row.phase01 = 0.0;
        v.add(row);
    }
    return v;
}

template <typename SampleType = float>
class AdditiveBankT {
public:
    enum class RetrigPhase : std::uint8_t {
        reset_stored,    
        free_run,        
        seeded_random,   
    };

    enum class EnvelopeMode : std::uint8_t {
        shared_ar,           
        per_partial_decay,   
    };
 
    // ── Calibration constants ─────────────────────────────────────────────
 
    static constexpr int kMaxPartialsCeiling = VoiceTable::kMaxPartials;
    static constexpr int kMaxPartialsDefault = 128;
    static constexpr int kPartialCountDefault = 64;
 
    static constexpr double kFundamentalMinHz = 20.0;
    static constexpr double kFundamentalMaxHz = 8000.0;
    static constexpr double kFundamentalDefaultHz = 220.0;
 
    static constexpr double kInharmonicityDefault = 0.0;
    static constexpr double kInharmonicityCitedMax = 0.02;
    static constexpr double kInharmonicityMax = 0.05;
    static constexpr double kBPianoTenor = 0.0002;
    static constexpr double kBPianoMid = 0.0008;
    static constexpr double kBUpperTreble = 0.004;
    static constexpr double kBElectricPianoTine = 0.02;
 
    static constexpr double kSpectralTiltMinDbOct = -24.0;
    static constexpr double kSpectralTiltMaxDbOct = 6.0;
    static constexpr double kSpectralTiltDefaultDbOct = -6.0;
 
    static constexpr double kMasterGainMinDb = -60.0;
    static constexpr double kMasterGainMaxDb = 6.0;
    static constexpr double kMasterGainDefaultDb = -6.0;
 
    static constexpr double kAttackDefaultMs = 5.0;
    static constexpr double kAttackMinMs = 0.1;
    static constexpr double kAttackMaxMs = 2000.0;
    static constexpr double kReleaseDefaultMs = 400.0;
    static constexpr double kReleaseMinMs = 1.0;
    static constexpr double kReleaseMaxMs = 20000.0;
 
    static constexpr double kDetuneDefaultCents = 0.0;
    static constexpr double kDetuneMaxCents = 50.0;
    static constexpr double kDoubletJitterSpread = 0.5;

    static constexpr int kAmpUpdateCadenceSamples = 32;

    static constexpr double kNyquistGuardBandHz = 1000.0;

    static constexpr double kEnvelopeCurve = 0.5;
 
    static constexpr double kTwoPi = 6.283185307179586476925286766559;
 
    // ── Static laws (shipped, so tests compute expectations from them) ────

    static double partial_frequency_hz(double f0_hz, double ratio,
                                       double inharmonicity_b) {
        return f0_hz * ratio *
               std::sqrt(1.0 + inharmonicity_b * ratio * ratio);
    }

    static double nyquist_guard_gain(double f_hz, double sample_rate) {
        const double nyquist = 0.5 * sample_rate;
        const double taper_start = nyquist - kNyquistGuardBandHz;
        if (f_hz <= taper_start) return 1.0;
        if (f_hz >= nyquist) return 0.0;
        const double t = (f_hz - taper_start) / kNyquistGuardBandHz;
        return 0.5 * (1.0 + std::cos(kPi() * t));
    }

    static double worst_case_gain() {
        return units::db_to_linear(kMasterGainMaxDb);
    }

    static constexpr int latency_samples() { return 0; }
 
    // ── Lifecycle ─────────────────────────────────────────────────────────
 
    AdditiveBankT() { prepare(44100.0, kMaxPartialsDefault); }

    void prepare(double sample_rate, int max_partials = kMaxPartialsDefault) {
        if (std::isfinite(sample_rate) && sample_rate > 0.0) sample_rate_ = sample_rate;
        max_partials_ = std::clamp(max_partials, 1, kMaxPartialsCeiling);
        partial_count_ = std::min(partial_count_, max_partials_);
 
        // Two slots per partial: a doublet is a genuine pair of oscillators,
        // not a modulation trick, so the storage is sized for it up front and
        // `prepare` stays the only allocation point whether or not the user
        // later engages the detune.
        const auto slots = static_cast<std::size_t>(2 * max_partials_);
        phase_.resize(slots);
        gain_.assign(slots, 0.0);
        gain_step_.assign(slots, 0.0);
        increment_.assign(slots, 0.0);
        increment_step_.assign(slots, 0.0);
        decay_env_.assign(slots, 1.0);
        decay_coef_.assign(slots, 1.0);
        slot_detune_.assign(slots, 1.0);
        render_slots_.assign(slots, 0);
 
        onset_.prepare(sample_rate_);
        onset_.set_curve(kEnvelopeCurve);
        refresh_onset();
        reset();
    }

    void reset() {
        rng_.reset();
        onset_.reset();
        for (auto& p : phase_) p.reset(0.0);
        std::fill(decay_env_.begin(), decay_env_.end(), 1.0);
        std::fill(gain_.begin(), gain_.end(), 0.0);
        std::fill(gain_step_.begin(), gain_step_.end(), 0.0);
        std::fill(increment_step_.begin(), increment_step_.end(), 0.0);
        draw_doublet_jitter();
        apply_stored_phases();
        cadence_countdown_ = 0;
        update_control_block(true);
    }
 
    // ── Global voice ──────────────────────────────────────────────────────
 
    void set_fundamental_hz(double f0) {
        if (!std::isfinite(f0)) return;
        fundamental_hz_ = std::clamp(f0, kFundamentalMinHz, kFundamentalMaxHz);
    }
    double fundamental_hz() const { return fundamental_hz_; }
 
    void set_partial_count(int n) {
        // `std::max(1, ...)` is not defensive padding: an empty voice table
        // would otherwise hand `std::clamp` a hi below its lo, which is UB.
        partial_count_ =
            std::clamp(n, 1, std::max(1, std::min(max_partials_, voice_.count)));
    }
    int partial_count() const { return partial_count_; }
    int max_partials() const { return max_partials_; }

    void set_inharmonicity_b(double b) {
        if (!std::isfinite(b)) return;
        inharmonicity_b_ = std::clamp(b, 0.0, kInharmonicityMax);
    }
    double inharmonicity_b() const { return inharmonicity_b_; }
 
    void set_spectral_tilt_db_oct(double s) {
        if (!std::isfinite(s)) return;
        tilt_db_oct_ = std::clamp(s, kSpectralTiltMinDbOct, kSpectralTiltMaxDbOct);
    }
    double spectral_tilt_db_oct() const { return tilt_db_oct_; }
 
    void set_master_gain_db(double g) {
        if (!std::isfinite(g)) return;
        master_gain_db_ = std::clamp(g, kMasterGainMinDb, kMasterGainMaxDb);
        master_linear_ = units::db_to_linear(master_gain_db_);
    }
    double master_gain_db() const { return master_gain_db_; }

    bool set_trig_profile(FastTrigProfile profile) noexcept {
        if (profile != FastTrigProfile::reference &&
            profile != FastTrigProfile::realtime_precise)
            return false;
        if constexpr (!std::is_same_v<SampleType, float>) {
            if (profile != FastTrigProfile::reference) return false;
        }
        trig_profile_ = profile;
        return true;
    }
    FastTrigProfile trig_profile() const noexcept { return trig_profile_; }
 
    static constexpr bool trig_profile_has_vector_path(FastTrigProfile profile) {
#if defined(__APPLE__) && defined(__clang__) && defined(__aarch64__)
        return std::is_same_v<SampleType, float> &&
               profile == FastTrigProfile::realtime_precise;
#else
        (void)profile;
        return false;
#endif
    }
 
    // ── The partial table ─────────────────────────────────────────────────
 
    void load_voice(const VoiceTable& v) {
        if (v.count < 0 || v.count > VoiceTable::kMaxPartials) return;
        for (int p = 0; p < v.count; ++p) {
            const auto& row = v.partials[static_cast<std::size_t>(p)];
            if (!std::isfinite(row.ratio) || !std::isfinite(row.amp) ||
                !std::isfinite(row.phase01) || !std::isfinite(row.decay_ms))
                return;
        }
        voice_ = v;
        partial_count_ = std::clamp(partial_count_, 1,
                                    std::max(1, std::min(max_partials_, voice_.count)));
    }
    const VoiceTable& voice() const { return voice_; }

    void set_partial(int p, double ratio, double amp_lin, double phase01,
                     double decay_ms) {
        if (p < 0 || p >= VoiceTable::kMaxPartials) return;
        if (!std::isfinite(ratio) || !std::isfinite(amp_lin) ||
            !std::isfinite(phase01) || !std::isfinite(decay_ms))
            return;
        voice_.partials[static_cast<std::size_t>(p)] = {ratio, amp_lin, phase01,
                                                        decay_ms};
        voice_.count = std::max(voice_.count, p + 1);
    }
 
    // ── Spectral envelope and morph ───────────────────────────────────────
 
    void set_envelope_a(const SpectralEnvelope& e) { envelope_a_ = e; }
    void set_envelope_b(const SpectralEnvelope& e) { envelope_b_ = e; }

    void set_morph(float m) {
        if (!std::isfinite(static_cast<double>(m))) return;
        morph_ = std::clamp(static_cast<double>(m), 0.0, 1.0);
    }
    double morph() const { return morph_; }
 
    void set_spectral_domain(SpectralDomain d) { domain_ = d; }
    SpectralDomain spectral_domain() const { return domain_; }

    double envelope_db_at(double f_hz) const {
        const double abscissa = domain_ == SpectralDomain::absolute_hz
                                    ? f_hz
                                    : f_hz / fundamental_hz_;
        const double a = envelope_a_.gain_db_at(abscissa);
        const double b = envelope_b_.gain_db_at(abscissa);
        return (1.0 - morph_) * a + morph_ * b +
               tilt_db_oct_ * std::log2(f_hz / fundamental_hz_);
    }
 
    // ── Envelope and trajectory ───────────────────────────────────────────
 
    void set_envelope_mode(EnvelopeMode mode) { envelope_mode_ = mode; }
    EnvelopeMode envelope_mode() const { return envelope_mode_; }
 
    void set_attack_ms(double ms) {
        if (!std::isfinite(ms)) return;
        // Zero is the explicit instantaneous-envelope mode used by the raw
        // oscillator/test surface. Non-zero values obey the published catalog
        // range, including its reachable 0.1 ms lower endpoint.
        attack_ms_ = ms == 0.0 ? 0.0 : std::clamp(ms, kAttackMinMs, kAttackMaxMs);
        refresh_onset();
    }
    void set_release_ms(double ms) {
        if (!std::isfinite(ms)) return;
        // As above, exact zero is a deliberate instantaneous internal mode;
        // all non-zero requests obey the declared 1 ms minimum.
        release_ms_ = ms == 0.0 ? 0.0 : std::clamp(ms, kReleaseMinMs, kReleaseMaxMs);
        refresh_onset();
    }

    void set_detune_cents(double cents) {
        if (!std::isfinite(cents)) return;
        detune_cents_ = std::clamp(cents, 0.0, kDetuneMaxCents);
    }
    double detune_cents() const { return detune_cents_; }
    bool doublet_active() const { return detune_cents_ > 0.0; }

    void set_pitch_glide(double start_cents, double time_ms) {
        if (!std::isfinite(start_cents) || !std::isfinite(time_ms)) return;
        glide_start_cents_ = start_cents;
        glide_time_ms_ = std::max(time_ms, 0.0);
    }
 
    void set_retrig_phase(RetrigPhase p) { retrig_phase_ = p; }
    RetrigPhase retrig_phase() const { return retrig_phase_; }

    void set_seed(std::uint32_t seed) {
        rng_.set_seed(seed);
        draw_doublet_jitter();
    }

    void retrigger() {
        switch (retrig_phase_) {
            case RetrigPhase::reset_stored:
                apply_stored_phases();
                break;
            case RetrigPhase::free_run:
                break;
            case RetrigPhase::seeded_random:
                for (int p = 0; p < partial_count_; ++p)
                    for (int k = 0; k < 2; ++k)
                        phase_[slot(p, k)].reset(rng_.next_unit<double>());
                break;
        }
        std::fill(decay_env_.begin(), decay_env_.end(), 1.0);
        glide_position_ = 0.0;
        onset_.gate_on();
        update_control_block(true);
    }

    void release() { onset_.gate_off(); }
 
    bool active() const { return onset_.active(); }

    double partial_frequency(int p) const {
        if (p < 0 || p >= partial_count_) return 0.0;
        return increment_[slot(p, 0)] * sample_rate_;
    }
 
    // ── Render ────────────────────────────────────────────────────────────
 
    SampleType next() {
        SampleType out{};
        process(&out, 1);
        return out;
    }
 
    void process(SampleType* out, int n) {
        switch (trig_profile_) {
            case FastTrigProfile::realtime_precise:
                process_with_trig<FastTrigProfile::realtime_precise>(out, n);
                return;
            case FastTrigProfile::reference:
            default:
                process_with_trig<FastTrigProfile::reference>(out, n);
                return;
        }
    }
 
private:
    template <FastTrigProfile Profile>
    void process_with_trig(SampleType* out, int n) {
        for (int i = 0; i < n; ++i) {
            if (cadence_countdown_ <= 0) update_control_block(false);
            --cadence_countdown_;
 
            double sum = 0.0;
#if defined(__APPLE__) && defined(__clang__) && defined(__aarch64__)
            if constexpr (std::is_same_v<SampleType, float> &&
                          Profile == FastTrigProfile::realtime_precise)
                sum = render_precise_simd_sample();
#endif
            constexpr bool rendered_simd =
#if defined(__APPLE__) && defined(__clang__) && defined(__aarch64__)
                std::is_same_v<SampleType, float> &&
                Profile == FastTrigProfile::realtime_precise;
#else
                false;
#endif
            for (int r = rendered_simd ? render_count_ : 0; r < render_count_; ++r) {
                const auto s = static_cast<std::size_t>(render_slots_[
                    static_cast<std::size_t>(r)]);
                // Evaluate at the CURRENT phase, then advance — the catalog's
                // oscillator convention (`osc::VaOscillator::next` reads its
                // entry phase). Advancing first would put the first sample one
                // increment past the stored phase, so a partial stored at 0.25
                // would start at 0.9983 rather than at 1, and a bank would not
                // agree sample-for-sample with a catalog sine osc set to the
                // same phase.
                const double sine = [&] {
                    if constexpr (Profile == FastTrigProfile::reference)
                        return std::sin(kTwoPi * phase_[s].phase());
                    else
                        return static_cast<double>(
                            FastMath::sin_cycles<Profile>(
                                static_cast<float>(phase_[s].phase())));
                }();
                sum += gain_[s] * decay_env_[s] * sine;
                phase_[s].advance(increment_[s]);
                gain_[s] += gain_step_[s];
                increment_[s] += increment_step_[s];
                decay_env_[s] *= decay_coef_[s];
            }
 
            const double env = static_cast<double>(onset_.next());
            // The accumulator is the one recursive-looking quantity here; the
            // per-partial decays are geometric and would otherwise drizzle
            // subnormals into the sum for minutes after a long tail.
            out[i] = static_cast<SampleType>(
                snap_to_zero(sum * env * master_linear_));
        }
    }
    static constexpr double kPi() { return 3.14159265358979323846; }

    static constexpr std::size_t slot(int p, int k) {
        return static_cast<std::size_t>(2 * p + k);
    }
 
#if defined(__APPLE__) && defined(__clang__) && defined(__aarch64__)
    double render_precise_simd_sample() {
        simd_float4 accumulated = simd_make_float4(0.0f);
        int r = 0;
        for (; r + 3 < render_count_; r += 4) {
            const auto s0 = static_cast<std::size_t>(render_slots_[r]);
            const auto s1 = static_cast<std::size_t>(render_slots_[r + 1]);
            const auto s2 = static_cast<std::size_t>(render_slots_[r + 2]);
            const auto s3 = static_cast<std::size_t>(render_slots_[r + 3]);
            const simd_float4 phase = simd_make_float4(
                static_cast<float>(phase_[s0].phase()),
                static_cast<float>(phase_[s1].phase()),
                static_cast<float>(phase_[s2].phase()),
                static_cast<float>(phase_[s3].phase()));
            const simd_float4 sine = FastMath::sin_cycles_precise(phase);
            const simd_float4 amplitude = simd_make_float4(
                static_cast<float>(gain_[s0] * decay_env_[s0]),
                static_cast<float>(gain_[s1] * decay_env_[s1]),
                static_cast<float>(gain_[s2] * decay_env_[s2]),
                static_cast<float>(gain_[s3] * decay_env_[s3]));
            accumulated += amplitude * sine;
 
            const std::array<std::size_t, 4> slots{s0, s1, s2, s3};
            for (const auto s : slots) {
                phase_[s].advance(increment_[s]);
                gain_[s] += gain_step_[s];
                increment_[s] += increment_step_[s];
                decay_env_[s] *= decay_coef_[s];
            }
        }
        double sum = static_cast<double>(simd_reduce_add(accumulated));
        for (; r < render_count_; ++r) {
            const auto s = static_cast<std::size_t>(render_slots_[r]);
            sum += gain_[s] * decay_env_[s] *
                   static_cast<double>(FastMath::sin_cycles<
                                       FastTrigProfile::realtime_precise>(
                       static_cast<float>(phase_[s].phase())));
            phase_[s].advance(increment_[s]);
            gain_[s] += gain_step_[s];
            increment_[s] += increment_step_[s];
            decay_env_[s] *= decay_coef_[s];
        }
        return sum;
    }
#endif

    void refresh_onset() {
        onset_.set_attack_ms(attack_ms_);
        onset_.set_release_ms(release_ms_);
    }
 
    void apply_stored_phases() {
        for (int p = 0; p < std::max(partial_count_, 1); ++p) {
            const double stored =
                voice_.partials[static_cast<std::size_t>(p)].phase01;
            for (int k = 0; k < 2; ++k) {
                const auto s = slot(p, k);
                if (s < phase_.size()) phase_[s].reset(stored);
            }
        }
    }

    void draw_doublet_jitter() {
        for (std::size_t s = 0; s < slot_detune_.size(); ++s)
            slot_detune_[s] = 1.0 + kDoubletJitterSpread *
                                        rng_.next_bipolar<double>();
    }

    void update_control_block(bool immediate) {
        cadence_countdown_ = kAmpUpdateCadenceSamples;
 
        const double glide = glide_ratio();
        advance_glide();
 
        const bool pair = doublet_active();
        const int sub_slots = pair ? 2 : 1;
        const double pair_scale = pair ? 0.5 : 1.0;
 
        // Pass one: frequencies, guard gains, envelope gains, and the sum the
        // normaliser needs. The normaliser has to see every factor that reaches
        // the output, or the bound it computes is a bound on a different signal.
        double magnitude_sum = 0.0;
        render_count_ = 0;
        for (int p = 0; p < partial_count_; ++p) {
            const auto& row = voice_.partials[static_cast<std::size_t>(p)];
            const double base_hz =
                voice_.harmonic
                    ? partial_frequency_hz(fundamental_hz_, row.ratio,
                                           inharmonicity_b_)
                    : fundamental_hz_ * row.ratio;
 
            for (int k = 0; k < sub_slots; ++k) {
                const auto s = slot(p, k);
                // The jitter is indexed by PARTIAL, not by slot, so the pair
                // straddles its nominal symmetrically. Drawing it per slot
                // would detune the two members by different amounts, which
                // shifts the mode's centre frequency instead of splitting it —
                // a mistuned bell rather than a beating one.
                const double cents =
                    pair ? (k == 0 ? -0.5 : 0.5) * detune_cents_ *
                               slot_detune_[slot(p, 0)]
                         : 0.0;
                const double f_hz =
                    base_hz * glide * units::cents_to_ratio(cents);
 
                const double guard = nyquist_guard_gain(f_hz, sample_rate_);
                const double env_db = f_hz > 0.0 ? envelope_db_at(f_hz) : 0.0;
                const double linear = units::db_to_linear(env_db);
                const double signed_gain =
                    row.amp * pair_scale * linear * guard;
 
                target_gain_[s] = signed_gain;
                target_increment_[s] = f_hz / sample_rate_;
                magnitude_sum += std::abs(signed_gain);
 
                // `shared_ar` means exactly that: the onset is the whole
                // amplitude story and the table's per-partial τ is ignored, so
                // every partial holds while gated. Reading `decay_ms` in this
                // mode would make the organ table's sustained rows and the
                // bell's 500 ms rows behave differently under a control the
                // user thinks selects between them.
                decay_coef_[s] =
                    envelope_mode_ == EnvelopeMode::per_partial_decay
                        ? decay_coefficient(row.decay_ms)
                        : 1.0;
                render_slots_[static_cast<std::size_t>(render_count_++)] =
                    static_cast<int>(s);
            }
        }
 
        // Pass two: the normaliser, then the ramps. Applying `g_norm` to the
        // per-partial gains rather than to the summed output is what makes the
        // bound survive the ramp — a post-sum normaliser would be a different
        // number on each side of a control-block edge.
        normaliser_ = 1.0 / std::max(1.0, magnitude_sum);
        const double step = 1.0 / static_cast<double>(kAmpUpdateCadenceSamples);
        for (int r = 0; r < render_count_; ++r) {
            const auto s = static_cast<std::size_t>(
                render_slots_[static_cast<std::size_t>(r)]);
            const double g = target_gain_[s] * normaliser_;
            if (immediate) {
                gain_[s] = g;
                increment_[s] = target_increment_[s];
                gain_step_[s] = 0.0;
                increment_step_[s] = 0.0;
            } else {
                gain_step_[s] = (g - gain_[s]) * step;
                increment_step_[s] = (target_increment_[s] - increment_[s]) * step;
            }
        }
    }

    double decay_coefficient(double decay_ms) const {
        if (!(decay_ms > 0.0)) return 1.0;   // sustained — the table's "∞"
        const double tau_samples = decay_ms * 0.001 * sample_rate_;
        return tau_samples > 0.0 ? std::exp(-1.0 / tau_samples) : 0.0;
    }
 
    double glide_ratio() const {
        if (glide_time_ms_ <= 0.0 || glide_start_cents_ == 0.0) return 1.0;
        const double t = std::min(glide_position_, 1.0);
        return units::cents_to_ratio(glide_start_cents_ * (1.0 - t));
    }
 
    void advance_glide() {
        if (glide_time_ms_ <= 0.0) {
            glide_position_ = 1.0;
            return;
        }
        const double block_ms = 1000.0 *
                                static_cast<double>(kAmpUpdateCadenceSamples) /
                                sample_rate_;
        glide_position_ =
            std::min(1.0, glide_position_ + block_ms / glide_time_ms_);
    }
 
    double sample_rate_ = 44100.0;
    int max_partials_ = kMaxPartialsDefault;
    int partial_count_ = kPartialCountDefault;
 
    double fundamental_hz_ = kFundamentalDefaultHz;
    double inharmonicity_b_ = kInharmonicityDefault;
    double tilt_db_oct_ = kSpectralTiltDefaultDbOct;
    double master_gain_db_ = kMasterGainDefaultDb;
    double master_linear_ = units::db_to_linear(kMasterGainDefaultDb);
    double morph_ = 0.0;
    double detune_cents_ = kDetuneDefaultCents;
    double attack_ms_ = kAttackDefaultMs;
    double release_ms_ = kReleaseDefaultMs;
    double glide_start_cents_ = 0.0;
    double glide_time_ms_ = 0.0;
    double glide_position_ = 1.0;
    double normaliser_ = 1.0;
 
    SpectralDomain domain_ = SpectralDomain::absolute_hz;
    EnvelopeMode envelope_mode_ = EnvelopeMode::per_partial_decay;
    RetrigPhase retrig_phase_ = RetrigPhase::reset_stored;
    FastTrigProfile trig_profile_ = FastTrigProfile::reference;
 
    VoiceTable voice_ = make_organ_voice(kPartialCountDefault);
    SpectralEnvelope envelope_a_{};
    SpectralEnvelope envelope_b_{};
 
    ArT<double> onset_{};
    Xorshift32 rng_{0x9E3779B9u};
 
    std::vector<osc::PhaseAccumulator> phase_{};
    std::vector<double> gain_{}, gain_step_{};
    std::vector<double> increment_{}, increment_step_{};
    std::vector<double> decay_env_{}, decay_coef_{};
    std::vector<double> slot_detune_{};
    std::vector<int> render_slots_{};
    std::array<double, 2 * VoiceTable::kMaxPartials> target_gain_{};
    std::array<double, 2 * VoiceTable::kMaxPartials> target_increment_{};
 
    int render_count_ = 0;
    int cadence_countdown_ = 0;
};
 
using AdditiveBank = AdditiveBankT<float>;
using AdditiveBank64 = AdditiveBankT<double>;
 
}  // namespace pulp::signal