TubePreamp.h

Go to the documentation of this file.

// DSPark - Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi - MIT License
 
#pragma once
 
#include "../Core/AudioBuffer.h"
#include "../Core/AudioSpec.h"
#include "../Core/Biquad.h"
#include "../Core/DenormalGuard.h"
#include "../Core/DspMath.h"
#include "../Core/Oversampling.h"
#include "../Core/StateBlob.h"
#include "../Core/WDF.h"
 
#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <memory>
#include <numbers>
#include <vector>
 
namespace dspark {
 
template <FloatType T>
class TubePreamp
{
public:
    // -- Lifecycle ---------------------------------------------------------------
 
    void prepare(const AudioSpec& spec)
    {
        if (!spec.isValid()) return;
        prepared_.store(false, std::memory_order_relaxed);
        spec_ = spec;
        sampleRate_ = spec.sampleRate;
        // Internal processing rate = active oversampling factor x base rate
        // (RF-009/ADR-011: factor is configurable, 1 = off).
        fs2_ = static_cast<double>(osFactor_) * sampleRate_;
        numChannels_ = spec.numChannels;
        maxBlock_ = std::max(spec.maxBlockSize, 1);
 
        if (osFactor_ > 1)
        {
            oversampler_ = std::make_unique<Oversampling<T>>(
                osFactor_, Oversampling<T>::Quality::High);
            oversampler_->prepare(spec);
        }
        else
        {
            oversampler_.reset();
        }
 
        channels_.clear();
        channels_.resize(static_cast<size_t>(numChannels_));
        for (auto& ch : channels_)
            ch = std::make_unique<ChannelState>(fs2_);
 
        latency_ = oversampler_ ? oversampler_->getLatency() : 0;
        drySize_ = 1;
        while (drySize_ < latency_ + maxBlock_ + 1) drySize_ <<= 1;
        dryRing_.assign(static_cast<size_t>(numChannels_),
                        std::vector<T>(static_cast<size_t>(drySize_), T(0)));
        dryPos_ = 0;
 
        calibrateReference();
 
        prepared_.store(true, std::memory_order_relaxed);
        dirty_.store(true, std::memory_order_release);
        reset();
    }
 
    void reset() noexcept
    {
        if (!prepared_.load(std::memory_order_relaxed)) return;
        const double sagR = static_cast<double>(sag_.load(std::memory_order_relaxed)) * 40e3;
        const int stages = stages_.load(std::memory_order_relaxed);
        for (auto& ch : channels_)
            ch->reset(sagR, stages);
        for (auto& d : dryRing_)
            std::fill(d.begin(), d.end(), T(0));
        dryPos_ = 0;
        if (oversampler_) oversampler_->reset();
        // Seed the anti-zipper ramps at their targets: no fade-in on start.
        hScaleSm_ = -1.0;
        outGainSm_ = -1.0;
        currentMix_ = mix_.load(std::memory_order_relaxed);
    }
 
    // -- Parameters (thread-safe) ---------------------------------------------------
 
    void setDrive(T db) noexcept
    {
        if (!std::isfinite(db)) return;
        driveDb_.store(std::clamp(db, T(-12), T(36)), std::memory_order_relaxed);
        dirty_.store(true, std::memory_order_release);
    }
 
    void setTreble(T treble) noexcept
    {
        if (!std::isfinite(treble)) return;
        treble_.store(std::clamp(treble, T(0), T(1)), std::memory_order_relaxed);
        dirty_.store(true, std::memory_order_release);
    }
 
    void setBass(T bass) noexcept
    {
        if (!std::isfinite(bass)) return;
        bass_.store(std::clamp(bass, T(0), T(1)), std::memory_order_relaxed);
        dirty_.store(true, std::memory_order_release);
    }
 
    void setMiddle(T middle) noexcept
    {
        if (!std::isfinite(middle)) return;
        middle_.store(std::clamp(middle, T(0), T(1)), std::memory_order_relaxed);
        dirty_.store(true, std::memory_order_release);
    }
 
    void setSag(T sag) noexcept
    {
        if (!std::isfinite(sag)) return;
        sag_.store(std::clamp(sag, T(0), T(1)), std::memory_order_relaxed);
        dirty_.store(true, std::memory_order_release);
    }
 
    void setStages(int stages) noexcept
    {
        stages_.store(std::clamp(stages, 1, 2), std::memory_order_relaxed);
        dirty_.store(true, std::memory_order_release);
    }
 
    void setOversampling(int factor)
    {
        if (factor < 1 || factor > 16 || (factor & (factor - 1)) != 0) return;
        if (factor == osFactor_) return;
        osFactor_ = factor;
        // Rebuild the whole chain at the new internal rate if already prepared
        // (fs2_, oversampler, per-channel circuits and calibration all depend
        // on the factor). Same setup-thread cost as prepare().
        if (prepared_.load(std::memory_order_relaxed))
            prepare(spec_);
    }
 
    [[nodiscard]] int getOversamplingFactor() const noexcept { return osFactor_; }
 
    void setOutput(T db) noexcept
    {
        if (!std::isfinite(db)) return;
        outputDb_.store(std::clamp(db, T(-24), T(12)), std::memory_order_relaxed);
    }
 
    void setMix(T mix) noexcept
    {
        if (!std::isfinite(mix)) return;
        mix_.store(std::clamp(mix, T(0), T(1)), std::memory_order_relaxed);
    }
 
    [[nodiscard]] T getDrive() const noexcept { return driveDb_.load(std::memory_order_relaxed); }
    [[nodiscard]] T getTreble() const noexcept { return treble_.load(std::memory_order_relaxed); }
    [[nodiscard]] T getBass() const noexcept { return bass_.load(std::memory_order_relaxed); }
    [[nodiscard]] T getMiddle() const noexcept { return middle_.load(std::memory_order_relaxed); }
    [[nodiscard]] T getSag() const noexcept { return sag_.load(std::memory_order_relaxed); }
    [[nodiscard]] int getStages() const noexcept { return stages_.load(std::memory_order_relaxed); }
    [[nodiscard]] T getOutput() const noexcept { return outputDb_.load(std::memory_order_relaxed); }
    [[nodiscard]] T getMix() const noexcept { return mix_.load(std::memory_order_relaxed); }
 
    [[nodiscard]] int getLatency() const noexcept { return latency_; }
 
    [[nodiscard]] int getLatencySamples() const noexcept { return latency_; }
 
    [[nodiscard]] T getSupplyVoltage() const noexcept
    {
        return supplyNow_.load(std::memory_order_relaxed);
    }
 
    [[nodiscard]] std::vector<uint8_t> getState() const
    {
        StateWriter w(stateId("TUBE"), 1);
        // Explicit float casts: the blob stores float, and with T = double the
        // unqualified write(key, double) would be ambiguous (float/int32/bool).
        w.write("drive", static_cast<float>(driveDb_.load(std::memory_order_relaxed)));
        w.write("treble", static_cast<float>(treble_.load(std::memory_order_relaxed)));
        w.write("bass", static_cast<float>(bass_.load(std::memory_order_relaxed)));
        w.write("middle", static_cast<float>(middle_.load(std::memory_order_relaxed)));
        w.write("sag", static_cast<float>(sag_.load(std::memory_order_relaxed)));
        w.write("stages", stages_.load(std::memory_order_relaxed));
        w.write("output", static_cast<float>(outputDb_.load(std::memory_order_relaxed)));
        w.write("mix", static_cast<float>(mix_.load(std::memory_order_relaxed)));
        w.write("oversampling", osFactor_);
        return w.blob();
    }
 
    bool setState(const uint8_t* data, size_t size)
    {
        StateReader r(data, size);
        if (!r.isValid() || r.processorId() != stateId("TUBE")) return false;
        setDrive(static_cast<T>(r.read("drive", 0.0f)));
        setTreble(static_cast<T>(r.read("treble", 0.5f)));
        setBass(static_cast<T>(r.read("bass", 0.5f)));
        setMiddle(static_cast<T>(r.read("middle", 0.5f)));
        setSag(static_cast<T>(r.read("sag", 0.3f)));
        setStages(r.read("stages", 1));
        setOutput(static_cast<T>(r.read("output", 0.0f)));
        setMix(static_cast<T>(r.read("mix", 1.0f)));
        // Default 2 = the historical fixed factor (blobs written before RF-009
        // restore the 2x behaviour they were captured with).
        setOversampling(r.read("oversampling", 2));
        return true;
    }
 
    // -- Processing -------------------------------------------------------------------
 
    void processBlock(AudioBufferView<T> buffer) noexcept
    {
        if (!prepared_.load(std::memory_order_relaxed)) return;
        DenormalGuard guard;
 
        const int nCh = std::min(buffer.getNumChannels(), numChannels_);
        const int nS = buffer.getNumSamples();
        if (nCh == 0 || nS == 0) return;
 
        // Front-door non-finite guard: a single NaN/Inf input sample would
        // poison the recursive triode Newton-Raphson / WDF tone stack / supply
        // sag / output DC-blocker / flatten-EQ state PERMANENTLY (only reset()
        // clears it, not clean input). Replace bad samples with silence before
        // they reach any state, so a transient upstream glitch cannot corrupt
        // the channel for the rest of the stream (M-005 C1).
        for (int ch = 0; ch < nCh; ++ch)
        {
            T* d = buffer.getChannel(ch);
            for (int i = 0; i < nS; ++i)
                if (!std::isfinite(d[i])) d[i] = T(0);
        }
 
        // Acquire pairs with the setters' release stores so the recompute
        // always sees the values published before the flag.
        if (dirty_.load(std::memory_order_relaxed)
            && dirty_.exchange(false, std::memory_order_acquire))
            recompute();
 
        // Linear per-block mix ramp with exact landing (settled: step == 0
        // and the per-sample value reduces to the constant, bit-identically).
        const T mixTarget = mix_.load(std::memory_order_relaxed);
        const T mixStart  = currentMix_;
        const T mixStep   = (mixTarget - mixStart) / static_cast<T>(nS);
        const double outGain = mScale_
            * std::pow(10.0, static_cast<double>(outputDb_.load(std::memory_order_relaxed)) / 20.0);
 
        // Dry snapshot.
        for (int ch = 0; ch < nCh; ++ch)
        {
            const T* in = buffer.getChannel(ch);
            auto& dry = dryRing_[static_cast<size_t>(ch)];
            int dp = dryPos_;
            for (int i = 0; i < nS; ++i)
            {
                dry[static_cast<size_t>(dp)] = in[i];
                dp = (dp + 1) & (drySize_ - 1);
            }
        }
 
        // Nonlinear circuit at the active oversampling factor (1x = process the
        // base-rate buffer in place, no resampling; RF-009/ADR-011).
        {
            const bool osOn = (oversampler_ != nullptr);
            auto osView = osOn ? oversampler_->upsample(buffer) : buffer;
            const int osN = osView.getNumSamples();
 
            // Anti-zipper: GEOMETRIC in-block ramps toward the current
            // targets (~30 ms across blocks), shared by all channels. A flat
            // per-block step clicked while dragging drive; and since
            // (hScale, outGain) is a compensation pair spanning orders of
            // magnitude (outGain ~ 1/drive), linear interpolation transits
            // through over-gained states - power-law interpolation keeps
            // the pair on the calibration curve.
            if (outGainSm_ <= 0.0) outGainSm_ = outGain;
            if (hScaleSm_ <= 0.0)  hScaleSm_ = hScale_;
            const double kSm = 1.0 - std::exp(-static_cast<double>(osN) / (0.030 * fs2_));
            const double outEnd = outGainSm_ * std::pow(outGain / outGainSm_, kSm);
            const double hEnd   = hScaleSm_ * std::pow(hScale_ / hScaleSm_, kSm);
            const double outRat = std::pow(outEnd / outGainSm_, 1.0 / static_cast<double>(osN));
            const double hRat   = std::pow(hEnd / hScaleSm_, 1.0 / static_cast<double>(osN));
 
            for (int ch = 0; ch < nCh; ++ch)
            {
                T* d = osView.getChannel(ch);
                auto& state = *channels_[static_cast<size_t>(ch)];
                double g = outGainSm_, h = hScaleSm_;
                for (int i = 0; i < osN; ++i)
                {
                    g *= outRat;
                    h *= hRat;
                    d[i] = static_cast<T>(g
                        * state.processSample(h * static_cast<double>(d[i]),
                                              numStagesActive_, sagR_));
                }
            }
            outGainSm_ = outEnd;
            hScaleSm_ = hEnd;
            if (osOn) oversampler_->downsample(buffer);
            supplyNow_.store(static_cast<T>(kBplus - sagR_ * channels_[0]->ipLP),
                             std::memory_order_relaxed);
        }
 
        // Latency-compensated mix (ramped: a hard flip on the distorted wet
        // stream clicked at 1.6x the steady-state sample delta).
        for (int ch = 0; ch < nCh; ++ch)
        {
            T* d = buffer.getChannel(ch);
            const auto& dry = dryRing_[static_cast<size_t>(ch)];
            for (int i = 0; i < nS; ++i)
            {
                const int idx = (dryPos_ + i - latency_) & (drySize_ - 1);
                const T drySample = dry[static_cast<size_t>(idx)];
                const T mixVal = mixStart + mixStep * static_cast<T>(i);
                d[i] = drySample + (d[i] - drySample) * mixVal;
            }
        }
        currentMix_ = mixTarget;   // exact landing
        dryPos_ = (dryPos_ + nS) & (drySize_ - 1);
    }
 
private:
    // -- Circuit constants (classic 12AX7 common-cathode stage) -------------------
    static constexpr double kMu = 100.0, kEx = 1.4, kKg1 = 1060.0;
    static constexpr double kKp = 600.0, kKvb = 300.0;
    static constexpr double kBplus = 300.0;
    static constexpr double kRL = 100e3;
    static constexpr double kRk = 1.5e3;
    static constexpr double kCk = 22e-6;
    static constexpr double kInterstage = 0.12;   
 
    static void koren(double vpk, double vgk, double& ip,
                      double& dIpdVpk, double& dIpdVgk) noexcept
    {
        vpk = std::max(vpk, 0.0);
        const double s = std::sqrt(kKvb + vpk * vpk);
        const double u = kKp * (1.0 / kMu + vgk / s);
 
        double sp = 0.0, sig = 0.0;                  // softplus(u), logistic(u)
        if (u > 30.0)       { sp = u; sig = 1.0; }
        else if (u < -30.0) { sp = std::exp(u); sig = sp; }
        else
        {
            const double eu = std::exp(u);
            sp = std::log1p(eu);
            sig = eu / (1.0 + eu);
        }
 
        const double e1 = (vpk / kKp) * sp;
        if (e1 <= 0.0)
        {
            ip = 0.0;
            dIpdVpk = 0.0;
            dIpdVgk = 0.0;
            return;
        }
        const double e1ex1 = std::pow(e1, kEx - 1.0);
        ip = 2.0 * e1ex1 * e1 / kKg1;
        const double dIpdE1 = 2.0 * kEx * e1ex1 / kKg1;
 
        const double dUdVpk = -kKp * vgk * vpk / (s * s * s);
        const double dE1dVpk = sp / kKp + (vpk / kKp) * sig * dUdVpk;
        const double dE1dVgk = vpk * sig / s;
        dIpdVpk = dIpdE1 * dE1dVpk;
        dIpdVgk = dIpdE1 * dE1dVgk;
    }
 
    struct TriodeStage
    {
        double fs2 = 96000.0;
        double ip = 8e-4;          
        double vk = 1.2;           
        double fPrev = 0.0;        
        double vpDC = 200.0;       
 
        void settleDC(double bplus) noexcept
        {
            // Static operating point: Vk = Ip*Rk (capacitor fully charged).
            double i = 8e-4;
            for (int it = 0; it < 60; ++it)
            {
                const double vkS = i * kRk;
                const double vpk = bplus - i * kRL - vkS;
                double ipK = 0.0, dVpk = 0.0, dVgk = 0.0;
                koren(vpk, -vkS, ipK, dVpk, dVgk);
                const double f = i - ipK;
                const double fp = 1.0 - (dVpk * (-(kRL + kRk)) + dVgk * (-kRk));
                const double di = f / fp;
                i -= di;
                i = std::clamp(i, 0.0, bplus / (kRL + kRk));
                if (std::abs(di) < 1e-15) break;
            }
            ip = i;
            vk = i * kRk;
            fPrev = 0.0;
            vpDC = bplus - i * kRL;
        }
 
        [[nodiscard]] double processSample(double vg, double bplusEff) noexcept
        {
            // Soft grid-conduction clamp toward +0.7 V.
            if (vg > 0.0)
                vg = 0.7 * std::tanh(vg / 0.7);
 
            // Trapezoidal cathode bypass: Ck dVk/dt = Ip - Vk/Rk, with fPrev
            // holding the previous NET CURRENT (Ip - Vk/Rk), so the update is
            //   Vk_n = Vk_{n-1} + (T/2Ck)(I_n + I_{n-1}) = kA + kB * Ip_n.
            const double h2c = 0.5 / (fs2 * kCk);
            const double denom = 1.0 + h2c / kRk;
            const double kA = (vk + h2c * fPrev) / denom;
            const double kB = h2c / denom;
 
            const double iMax = bplusEff / kRL + 1e-3;
            double i = std::clamp(ip, 0.0, iMax);
 
            for (int it = 0; it < 8; ++it)
            {
                const double vkN = kA + kB * i;
                const double vpk = bplusEff - i * kRL - vkN;
                const double vgk = vg - vkN;
                double ipK = 0.0, dVpk = 0.0, dVgk = 0.0;
                koren(vpk, vgk, ipK, dVpk, dVgk);
 
                const double f = i - ipK;
                const double fp = 1.0 - (dVpk * (-kRL - kB) + dVgk * (-kB));
                const double di = f / std::max(fp, 1e-6);
                i -= di;
                i = std::clamp(i, 0.0, iMax);
                if (std::abs(di) < 1e-12) break;
            }
 
            const double vkN = kA + kB * i;
            fPrev = i - vkN / kRk;      // net capacitor current for the trapezoid
            vk = vkN;
            ip = i;
            return (bplusEff - i * kRL) - vpDC;     // AC component
        }
    };
 
    struct ChannelState
    {
        explicit ChannelState(double fs2In)
            : fmv(38e3, 1e6)   // driven by the stage's real output impedance
        {
            stage1.fs2 = fs2In;
            stage2.fs2 = fs2In;
            fs2 = fs2In;
            fmv.prepare(fs2In);
        }
 
        void reset(double sagR, int numStages) noexcept
        {
            // Fixed point over the ACTIVE stage count: processSample only
            // draws current from the stages in use, so seeding ipLP with both
            // stages' current at 1-stage settings left a ~70 ms sag transient
            // (audible activation drift, and it polluted the noise floor).
            double bp = kBplus;
            double iTotal = 0.0;
            for (int it = 0; it < 8; ++it)
            {
                stage1.settleDC(bp);
                stage2.settleDC(bp);
                iTotal = stage1.ip + (numStages > 1 ? stage2.ip : 0.0);
                const double bpNew = kBplus - sagR * iTotal;
                if (std::abs(bpNew - bp) < 1e-9) break;
                bp = bpNew;
            }
            ipLP = iTotal;
            outHpX = outHpY = 0.0;
            fmv.reset();
            for (auto& set : flatten)
                for (auto& f : set)
                    f.reset();
        }
 
        void setFlattenCoeffs(int stageCount,
                              const std::array<BiquadCoeffs, 3>& c) noexcept
        {
            auto& set = flatten[static_cast<size_t>(stageCount - 1)];
            for (int k = 0; k < 3; ++k)
                set[static_cast<size_t>(k)].setCoeffs(c[static_cast<size_t>(k)]);
        }
 
        void setToneControls(double t, double b, double m) noexcept
        {
            fmv.setControls(t, b, m);   // rebuilds the R-type scattering
        }
 
        [[nodiscard]] double processSample(double vgIn, int numStages, double sagR) noexcept
        {
            // Supply sag: B+ droops with smoothed total plate current.
            const double sagAlpha = 1.0 - std::exp(-1.0 / (0.07 * fs2));
            const double iTotal = stage1.ip + (numStages > 1 ? stage2.ip : 0.0);
            ipLP += sagAlpha * (iTotal - ipLP);
            const double bplusEff = kBplus - sagR * ipLP;
 
            // Stage 1 -> FMV tone stack -> (stage 2) -> output high-pass.
            double v = stage1.processSample(vgIn, bplusEff);
            v = fmv.processSample(v);
            if (numStages > 1)
                v = stage2.processSample(v * kInterstage, bplusEff);
 
            // Output coupling high-pass (~8 Hz, removes residual sag drift).
            const double a = 1.0 - 2.0 * std::numbers::pi * 8.0 / fs2;
            const double y = a * (outHpY + v - outHpX);
            outHpX = v;
            outHpY = y;
            // Restore absolute polarity for single-stage use.
            double out = (numStages > 1) ? y : -y;
 
            // Reference-flattening EQ: undoes the FMV stack's fixed envelope
            // at the neutral tone setting (designed in calibrateReference
            // from the measured response), so neutral knobs sound neutral
            // and the tone controls act RELATIVE to flat. The triode's
            // harmonic character is untouched (this stage is linear).
            auto& fl = flatten[numStages > 1 ? 1 : 0];
            out = fl[0].processSample(out, 0);
            out = fl[1].processSample(out, 0);
            out = fl[2].processSample(out, 0);
            return out;
        }
 
        double fs2 = 96000.0;
        TriodeStage stage1, stage2;
        double ipLP = 1.6e-3;
        double outHpX = 0.0, outHpY = 0.0;
        std::array<std::array<Biquad<double, 1>, 3>, 2> flatten;   
 
        wdf::ToneStackFMV<double> fmv;   
    };
 
    void recompute() noexcept
    {
        const double drive = std::pow(10.0, static_cast<double>(
            driveDb_.load(std::memory_order_relaxed)) / 20.0);
        const double t = static_cast<double>(treble_.load(std::memory_order_relaxed));
        const double b = static_cast<double>(bass_.load(std::memory_order_relaxed));
        const double m = static_cast<double>(middle_.load(std::memory_order_relaxed));
        const double sag = static_cast<double>(sag_.load(std::memory_order_relaxed));
        numStagesActive_ = stages_.load(std::memory_order_relaxed);
 
        hScale_ = drive;                       // 0 dBFS -> 1 V grid at drive 0
        // Note on physics: a class-A preamp draws near-constant average
        // current, so supply sag shifts the operating point rather than
        // pumping like a push-pull power amp. The audible touch response of
        // this model comes from the cathode-bypass bias shift (modelled in
        // TriodeStage); the sag control changes voicing, not loudness.
        sagR_ = sag * 40e3;
        for (auto& ch : channels_)
            ch->setToneControls(t, b, m);
 
        // Loudness: divide out the circuit's program gain at the REFERENCE
        // tone setting (measured once in prepare on a settled channel) and
        // MOST of the drive factor. The link is partial (drive^0.75, i.e. a
        // residual +0.25 dB/dB slope): an exact 1/drive link leaves the knob
        // audibly dead below 0 dB (the circuit is still clean there) and
        // turns it into a pure attenuator above (compression eats level
        // faster than the link returns it). With the residual slope, -12 dB
        // drive sits ~3 dB lower and clean, high drive holds level while the
        // density grows. The tone knobs stay fully audible: their deviation
        // from the 0.5/0.5/0.5 reference is part of the tone, not the level.
        // Cheap by construction (no scratch processing on the audio thread).
        //
        // The divisor is the circuit's MEASURED program gain at this drive
        // (prepare-time LUT, log-interpolated) - same contract as
        // TapeMachine/TransformerModel, whose per-drive scratch calibration
        // is what kept their knobs healthy. The previous analytic
        // 1/drive^0.75 broke at high drive: once the triode pins at its
        // ceiling, output stops growing with drive while the divisor keeps
        // rising, so the level FELL hard instead of holding.
        const double driveDbNow = static_cast<double>(driveDb_.load(std::memory_order_relaxed));
        mScale_ = std::pow(drive, 0.25)
                / std::max(programGainAt(driveDbNow, numStagesActive_), 1e-9);
    }
 
    [[nodiscard]] double programGainAt(double driveDb, int stages) const noexcept
    {
        const auto& lut = gProgLut_[static_cast<size_t>(stages - 1)];
        // Ordered min/max instead of clamp: a NaN input resolves to 0 here
        // (defence in depth; the setter already rejects non-finite values, and
        // an unguarded NaN would UB-cast into a wild LUT index).
        const double pos = std::min(std::max(0.0, (driveDb - kDriveLutMinDb) / kDriveLutStepDb),
                                    static_cast<double>(kDriveLutN - 1));
        const int idx = std::min(static_cast<int>(pos), kDriveLutN - 2);
        const double frac = pos - static_cast<double>(idx);
        const double a = std::max(lut[static_cast<size_t>(idx)], 1e-9);
        const double b = std::max(lut[static_cast<size_t>(idx) + 1], 1e-9);
        return a * std::pow(b / a, frac);   // gains are smooth in dB
    }
 
    void calibrateReference() noexcept
    {
        constexpr double kRefTones[] = { 100.0, 200.0, 400.0, 800.0, 1600.0, 3200.0, 6400.0 };
        constexpr double kSagRef = 0.3 * 40e3;
        const int settle = static_cast<int>(0.060 * fs2_);
        const int meas   = static_cast<int>(0.050 * fs2_);
        constexpr double kAmpPerTone = 0.095 / 2.6457513;
 
        for (int st = 1; st <= 2; ++st)
        {
            ChannelState cal(fs2_);   // fresh: flattener is passthrough here
            cal.setToneControls(0.5, 0.5, 0.5);
            cal.reset(kSagRef, st);
 
            // Per-tone Goertzel over the measured tail.
            std::array<double, 7> gs1 {}, gs2 {}, gc {};
            for (int k = 0; k < 7; ++k)
                gc[static_cast<size_t>(k)] =
                    2.0 * std::cos(2.0 * std::numbers::pi * kRefTones[k] / fs2_);
            for (int i = 0; i < settle + meas; ++i)
            {
                double x = 0.0;
                for (int k = 0; k < 7; ++k)
                    x += std::sin(2.0 * std::numbers::pi * kRefTones[k] * i / fs2_ + k * 1.7);
                x *= kAmpPerTone;
                const double y = cal.processSample(x, st, kSagRef);
                if (i >= settle)
                    for (int k = 0; k < 7; ++k)
                    {
                        auto& s1 = gs1[static_cast<size_t>(k)];
                        auto& s2 = gs2[static_cast<size_t>(k)];
                        const double s0 = y + gc[static_cast<size_t>(k)] * s1 - s2;
                        s2 = s1; s1 = s0;
                    }
            }
            std::array<double, 7> gainDb {};
            for (int k = 0; k < 7; ++k)
            {
                const double c = gc[static_cast<size_t>(k)] * 0.5;
                const double s1 = gs1[static_cast<size_t>(k)], s2 = gs2[static_cast<size_t>(k)];
                const double mag = std::sqrt(std::max(s1 * s1 + s2 * s2 - 2.0 * c * s1 * s2, 0.0))
                                 * 2.0 / meas;
                gainDb[static_cast<size_t>(k)] =
                    20.0 * std::log10(std::max(mag / kAmpPerTone, 1e-9));
            }
 
            // Flattening EQ from the band means, relative to the overall mean.
            double mean = 0.0;
            for (double g : gainDb) mean += g / 7.0;
            const double lo = 0.5 * (gainDb[0] + gainDb[1]) - mean;
            const double mid = (gainDb[2] + gainDb[3] + gainDb[4]) / 3.0 - mean;
            const double hi = 0.5 * (gainDb[5] + gainDb[6]) - mean;
 
            std::array<BiquadCoeffs, 3> fc = {
                BiquadCoeffs::makeLowShelf(fs2_, 180.0, -lo),
                BiquadCoeffs::makePeak(fs2_, 800.0, 0.55, -mid),
                BiquadCoeffs::makeHighShelf(fs2_, 4500.0, -hi)
            };
            for (auto& ch : channels_)
                ch->setFlattenCoeffs(st, fc);
 
            // Phase 2: program gain vs DRIVE, on a chained scratch channel
            // with the flattener installed (the chain as it really sounds).
            // Each grid point gets a re-settle (bias/sag adapt to the new
            // level) before its measurement window. recompute() interpolates
            // this LUT, so the loudness link tracks the circuit's actual
            // compression - the fix for the level falling at high drive.
            ChannelState sweep(fs2_);
            sweep.setToneControls(0.5, 0.5, 0.5);
            sweep.setFlattenCoeffs(st, fc);
            sweep.reset(kSagRef, st);
            const int settle2 = static_cast<int>(0.030 * fs2_);
            const int meas2   = static_cast<int>(0.030 * fs2_);
            int t = 0;   // continuous tone phase across the whole sweep
            for (int kd = 0; kd < kDriveLutN; ++kd)
            {
                const double d = std::pow(10.0,
                    (kDriveLutMinDb + kd * kDriveLutStepDb) / 20.0);
                double inSq = 0.0, outSq = 0.0;
                for (int i = 0; i < settle2 + meas2; ++i, ++t)
                {
                    double x = 0.0;
                    for (int k = 0; k < 7; ++k)
                        x += std::sin(2.0 * std::numbers::pi * kRefTones[k] * t / fs2_ + k * 1.7);
                    x *= kAmpPerTone;
                    const double y = sweep.processSample(d * x, st, kSagRef);
                    if (i >= settle2)
                    {
                        inSq += x * x;
                        outSq += y * y;
                    }
                }
                gProgLut_[static_cast<size_t>(st - 1)][static_cast<size_t>(kd)] =
                    (outSq > 0.0 && inSq > 0.0) ? std::sqrt(outSq / inSq) : 1.0;
            }
        }
    }
 
    // -- Members --------------------------------------------------------------------
    AudioSpec spec_ {};                
    double sampleRate_ = 48000.0;
    double fs2_ = 96000.0;
    int numChannels_ = 0;
    int maxBlock_ = 0;
    std::atomic<bool> prepared_ { false };
    int latency_ = 0;
    int drySize_ = 1;
    int osFactor_ = 2;                  
 
    std::unique_ptr<Oversampling<T>> oversampler_;
    std::vector<std::unique_ptr<ChannelState>> channels_;
 
    std::vector<std::vector<T>> dryRing_;
    int dryPos_ = 0;
 
    static constexpr int    kDriveLutN = 13;       
    static constexpr double kDriveLutMinDb = -12.0;
    static constexpr double kDriveLutStepDb = 4.0;
 
    double hScale_ = 1.0;
    double mScale_ = 1.0;
    double sagR_ = 0.0;
    int numStagesActive_ = 1;
    std::array<std::array<double, kDriveLutN>, 2> gProgLut_ {
        { { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 },
          { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 } } };
    double hScaleSm_ = -1.0;                    
    double outGainSm_ = -1.0;                   
    T currentMix_ = T(1);                       
 
    std::atomic<T> driveDb_ { T(0) };
    std::atomic<T> treble_ { T(0.5) };
    std::atomic<T> bass_ { T(0.5) };
    std::atomic<T> middle_ { T(0.5) };
    std::atomic<T> sag_ { T(0.3) };
    std::atomic<int> stages_ { 1 };
    std::atomic<T> outputDb_ { T(0) };
    std::atomic<T> mix_ { T(1) };
    std::atomic<T> supplyNow_ { T(300) };
    std::atomic<bool> dirty_ { true };
};
 
} // namespace dspark