BeatTracker.h

// DSPark -- Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi -- MIT License
 
#pragma once
 
#include "../Core/DspMath.h"
#include "../Core/AudioBuffer.h"
#include "../Core/AudioSpec.h"
#include "../Core/detail/AccountedAllocator.h"
#include "OnsetDetector.h"
 
#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <span>
#include <vector>
 
namespace dspark {
 
namespace detail { template <FloatType T> struct OfflineBeatEngine; }
 
template <FloatType T>
class BeatTracker final
{
    template <typename U>
    using WorkVector = std::vector<U, detail::AccountedAllocator<U>>;
    friend struct detail::OfflineBeatEngine<T>;
 
public:
    BeatTracker()
    {
        const detail::AccountedAllocator<double> doubles(&allocationAccount_);
        for (auto* values : {&bankPeriod_, &bankCosCoef_, &bankSinCoef_, &bankDecay_,
                             &bankPrior_, &zRe_, &zIm_, &zMass_, &bankCoherence_,
                             &bankLevelRatio_, &env_, &envBal_, &local_, &cumScore_,
                             &acf_, &scratchA_, &kernel_, &localPeriod_, &sweepA_,
                             &sweepB_, &candScore_, &candExplains_})
            *values = WorkVector<double>(doubles);
        for (auto& values : envReg_) values = WorkVector<double>(doubles);
        envRef_ = WorkVector<int64_t>(detail::AccountedAllocator<int64_t>(&allocationAccount_));
        backlink_ = WorkVector<int>(detail::AccountedAllocator<int>(&allocationAccount_));
        candidates_ = WorkVector<int>(detail::AccountedAllocator<int>(&allocationAccount_));
    }
 
    struct Result
    {
        T tempoBpm = T(0);           
        T confidence = T(0);         
        T secondaryTempoBpm = T(0);  
        std::vector<int64_t> beatSamples; 
    };
 
    // -- Lifecycle -----------------------------------------------------------
 
    void prepare(const AudioSpec& spec)
    {
        if (!(spec.sampleRate > 0.0) || !std::isfinite(spec.sampleRate))
            return;
 
        prepared_.store(false, std::memory_order_relaxed);
 
        sampleRate_ = spec.sampleRate;
        onset_.prepare(spec);
        hop_ = onset_.getHopSize();
        if (hop_ < 1) hop_ = 1;
        frameRate_ = sampleRate_ / static_cast<double>(hop_);
        baselineCoef_ = 1.0 - std::exp(-1.0 / std::max(1.0, kBaselineSeconds * frameRate_));
 
        latencySamples_.store(static_cast<int64_t>(onset_.getEnvelopeLatencySamples())
                                  + static_cast<int64_t>(hop_),
                              std::memory_order_relaxed);
 
        buildBank();
        setTempoRange(T(kDefaultMinBpm), T(kDefaultMaxBpm));
        tightness_.store(kDefaultTightness, std::memory_order_relaxed);
 
        // Offline scratch that does not depend on the signal length.
        candidates_.reserve(64);
 
        resetState();
        prepared_.store(true, std::memory_order_relaxed);
    }
 
    void setTempoRange(T minBpm, T maxBpm) noexcept
    {
        if (!std::isfinite(minBpm) || !std::isfinite(maxBpm)) return;
        if (!(minBpm > T(0)) || !(maxBpm > T(0)) || !(minBpm <= maxBpm)) return;
        if (bankSize_ < 2) return;
 
        const double lo = std::clamp(static_cast<double>(minBpm),
                                     kBankMinBpm, kBankMaxBpm);
        const double hi = std::clamp(static_cast<double>(maxBpm),
                                     kBankMinBpm, kBankMaxBpm);
 
        // Bank index runs with PERIOD, so the fastest tempo is the lowest
        // index: the high BPM bound selects the low index bound.
        int iLo = bpmToIndex(hi);
        int iHi = bpmToIndex(lo);
        iLo = std::clamp(iLo, 0, bankSize_ - 1);
        iHi = std::clamp(iHi, 0, bankSize_ - 1);
        if (iLo > iHi) std::swap(iLo, iHi);
 
        // The two indices are a set with an invariant (lo <= hi), so they are
        // published as one word: half of one range combined with half of the
        // next is a range nobody asked for, and it would be searched.
        activeRange_.store(packRange(iLo, iHi), std::memory_order_relaxed);
    }
 
    void setTightness(T alpha) noexcept
    {
        if (!std::isfinite(alpha)) return;
        tightness_.store(std::max(kMinTightness, static_cast<double>(alpha)),
                         std::memory_order_relaxed);
    }
 
    [[nodiscard]] T getTightness() const noexcept
    {
        return static_cast<T>(tightness_.load(std::memory_order_relaxed));
    }
 
    // -- Offline -------------------------------------------------------------
 
    Result analyze(AudioBufferView<const T> whole)
    {
        if (!prepared_.load(std::memory_order_relaxed)) return {};
        if (whole.getNumChannels() < 1 || whole.getNumSamples() <= 0) return {};
        const int n = whole.getNumSamples();
        beginOffline(n);
        pushOffline(std::span<const T>(whole.getChannel(0), static_cast<size_t>(n)));
        return finishOffline();
    }
 
    void beginOffline(int64_t expectedSamples = 0)
    {
        env_.clear();
        envRef_.clear();
        for (auto& r : envReg_) r.clear();
        offlineOpen_ = false;
        if (!prepared_.load(std::memory_order_relaxed)) return;
 
        resetState();
        if (expectedSamples > 0)
        {
            const size_t frames = static_cast<size_t>(expectedSamples / std::max(1, hop_) + 2);
            env_.reserve(frames);
            envRef_.reserve(frames);
            for (auto& r : envReg_) r.reserve(frames);
        }
        offlinePushed_ = 0;
        offlineOpen_ = true;
    }
 
    void pushOffline(std::span<const T> samples)
    {
        if (!offlineOpen_) return;
        size_t off = 0;
        while (off < samples.size())
        {
            // Pieces end on the front end's frame boundaries, so each frame's
            // envelope value is read as the frame is computed.
            const int64_t since = offlinePushed_ % static_cast<int64_t>(hop_);
            const size_t take = std::min(samples.size() - off,
                                         static_cast<size_t>(static_cast<int64_t>(hop_) - since));
            onset_.pushSamples(samples.subspan(off, take));
            offlinePushed_ += static_cast<int64_t>(take);
            off += take;
            if (offlinePushed_ % static_cast<int64_t>(hop_) == 0) appendEnvelopeFrame();
        }
    }
 
    Result finishOffline()
    {
        Result out;
        if (!offlineOpen_) return out;
        offlineOpen_ = false;
        onset_.reset();
 
        return finishEnvelope();
    }
 
    // -- Audio path (causal, RT-safe) ---------------------------------------
 
    void processBlock(AudioBufferView<const T> in) noexcept
    {
        if (!prepared_.load(std::memory_order_relaxed) || in.getNumChannels() < 1)
            return;
        const T* ch0 = in.getChannel(0);
        pushSamples(std::span<const T>(ch0, static_cast<size_t>(in.getNumSamples())));
    }
 
    void pushSamples(std::span<const T> samples) noexcept
    {
        if (!prepared_.load(std::memory_order_relaxed)) return;
 
        bool fired = false;
        size_t off = 0;
        while (off < samples.size())
        {
            const int64_t sinceFrame = samplesPushed_ % static_cast<int64_t>(hop_);
            const size_t toBoundary = static_cast<size_t>(static_cast<int64_t>(hop_)
                                                          - sinceFrame);
            const size_t take = std::min(samples.size() - off, toBoundary);
 
            onset_.pushSamples(samples.subspan(off, take));
            samplesPushed_ += static_cast<int64_t>(take);
            off += take;
 
            if (static_cast<size_t>(take) == toBoundary)
            {
                if (advanceFrame()) fired = true;
            }
        }
 
        // Release, and the matching acquire is in beatNow(). The latch and the
        // beat's position are two separate words written in this same call,
        // and a caller that sees the latch set must not then read a position
        // from an earlier call: on a weakly ordered machine two relaxed loads
        // of distinct objects may be observed in either order, so without this
        // pairing a reader could be handed the PREVIOUS beat's position -- an
        // error of one whole beat period -- while being told a beat just
        // happened. The store is last on purpose: everything the frame wrote,
        // the position included, is ordered before it.
        beatLatched_.store(fired, std::memory_order_release);
    }
 
    // -- Readout (lock-free) -------------------------------------------------
 
    [[nodiscard]] T getRunningTempoBpm() const noexcept
    {
        float bpm = 0.0f, conf = 0.0f;
        unpackTempo(tempoAndConfidence_.load(std::memory_order_relaxed), bpm, conf);
        return static_cast<T>(bpm);
    }
 
    [[nodiscard]] T getConfidence() const noexcept
    {
        float bpm = 0.0f, conf = 0.0f;
        unpackTempo(tempoAndConfidence_.load(std::memory_order_relaxed), bpm, conf);
        return static_cast<T>(conf);
    }
 
    void getTempoAndConfidence(T& bpmOut, T& confidenceOut) const noexcept
    {
        float bpm = 0.0f, conf = 0.0f;
        unpackTempo(tempoAndConfidence_.load(std::memory_order_relaxed), bpm, conf);
        bpmOut = static_cast<T>(bpm);
        confidenceOut = static_cast<T>(conf);
    }
 
    [[nodiscard]] bool beatNow() const noexcept
    {
        return beatLatched_.load(std::memory_order_acquire);
    }
 
    [[nodiscard]] int64_t getLastBeatSample() const noexcept
    {
        return lastBeatSample_.load(std::memory_order_relaxed);
    }
 
    [[nodiscard]] int getLatencySamples() const noexcept
    {
        return static_cast<int>(latencySamples_.load(std::memory_order_relaxed));
    }
 
    [[nodiscard]] double getFrameRate() const noexcept { return frameRate_; }
 
    [[nodiscard]] const OnsetDetector<T>& getOnsetDetector() const noexcept
    {
        return onset_;
    }
 
    void reset() noexcept
    {
        offlineOpen_ = false;
        resetState();
    }
 
private:
    // Shared complete-envelope engine. The legacy sample frontend and the
    // pooled offline feature path both execute this exact arithmetic.
    Result finishEnvelope(void* beatContext = nullptr,
                          void (*receiveBeats)(void*, std::span<const int64_t>) = nullptr)
    {
        Result out;
        int iLo = 0, iHi = 0;
        unpackRange(activeRange_.load(std::memory_order_relaxed), iLo, iHi);
 
        const int n = static_cast<int>(env_.size());
        const double maxPeriod = bankPeriod_[static_cast<size_t>(iHi)];
        if (n < static_cast<int>(kMinAnalysisBeats * maxPeriod))
        {
            resetState();
            return out;
        }
 
        conditionEnvelope();
 
        // Two candidate periods from the tempo-prior-weighted autocorrelation,
        // at least a quarter octave apart so that the runner-up is a different
        // metrical reading and not the same peak one lag over.
        double tau1 = 0.0, tau2 = 0.0;
        const bool balanced = (envBal_.size() == env_.size());
        if (balanced) env_.swap(envBal_);
        pickCandidates(iLo, iHi, tau1, tau2);
        if (balanced) env_.swap(envBal_);
        if (!(tau1 > 0.0)) { resetState(); return out; }
 
        // The ranking above proposes a reading; which multiple of it is the
        // beat is settled by the level model (see chooseMetricalLevel()). Only
        // where the signal has a pulse of its own in the range: otherwise the
        // correlation decided on its own, and so it stays.
        if (fundamentalFound_) chooseMetricalLevel(iLo, iHi, tau1, tau2);
 
        // The metrical level is decided against one period for the whole
        // signal, which is the stable model, and only then is the grid laid
        // down against a period that is allowed to move.
        // A period that moves is only worth estimating where there is a pulse
        // to estimate it from. When the winner does not clear the fundamental
        // floor -- the caller restricted the range past the real pulse -- a
        // per-frame estimate over that range is tracking noise, and it drags
        // the delivered grid with it: measured at 96 BPM on 160 BPM material
        // searched over 60 to 100, where the answer a listener would give is
        // the 80 BPM the correlation points at.
        //
        // Both grids are built - one against the moving period, one against
        // the steady one - and the one that explains the envelope better is
        // delivered. A per-frame period is the right model for a tempo that
        // moves, and the wrong one for a steady tempo under a syncopated
        // part: there the resonator sweep is drawn, bar by bar, towards the
        // syncopation's own period, and the grid it drives falls off the beat
        // and onto the off-accents. Measured on a limited pop master at 120
        // BPM over a dotted-eighth guitar, the moving-period grid put 18 of 74
        // intervals at three quarters or five quarters of a beat; the steady
        // grid has none, and explains more of the envelope, so it is the one
        // delivered. On a tempo ramp the steady grid drifts off the beats and
        // the moving one wins by the same measure.
        WorkVector<int64_t> win(detail::AccountedAllocator<int64_t>{&allocationAccount_});
        localPeriod_.clear();
        WorkVector<int64_t> moving(detail::AccountedAllocator<int64_t>{&allocationAccount_});
        bool haveMoving = false;
        if (coherence(tau1) >= kMinFundamentalCoherence)
        {
            computeLocalPeriods(tau1);
            haveMoving = buildGrid(tau1, moving) && moving.size() >= 2;
        }
        localPeriod_.clear();
        bool haveSteady = buildGrid(tau1, win) && win.size() >= 2;
        // A third reading: the steady period held kStrictTightnessFactor times
        // more tightly. The tightness in force is tuned for rubato, where it
        // must let an interval stretch; under a steady beat with a strong
        // syncopated part the same slack lets single intervals jump onto the
        // off-accents. Whichever grid explains the envelope best is kept, so
        // rubato keeps its loose grid and a steady groove gets a steady one.
        {
            WorkVector<int64_t> strict(detail::AccountedAllocator<int64_t>{&allocationAccount_});
            alphaScale_ = kStrictTightnessFactor;
            const bool haveStrict = buildGrid(tau1, strict) && strict.size() >= 2;
            alphaScale_ = 1.0;
            if (haveStrict && (!haveSteady || gridCoherence(strict) > gridCoherence(win)))
            {
                win.swap(strict);
                haveSteady = true;
            }
        }
        if (haveMoving && (!haveSteady
                           || gridCoherence(moving) > kMovingGridMargin * gridCoherence(win)))
            win.swap(moving);
        else if (!haveSteady)
        {
            resetState();
            return out;
        }
 
        // The delivered tempo is fitted to the delivered grid, not read off
        // the autocorrelation lag: the grid spans the whole signal, so a
        // straight-line fit through it resolves the period far below the
        // frame spacing the lag is quantised to.
        const double slope = fitBeatSlope(win);
        if (receiveBeats)
            receiveBeats(beatContext, {win.data(), win.size()});
        else
            out.beatSamples.assign(win.begin(), win.end());
        out.tempoBpm = static_cast<T>((slope > 0.0) ? 60.0 * sampleRate_ / slope
                                                    : periodToBpm(tau1));
        out.secondaryTempoBpm = static_cast<T>((tau2 > 0.0) ? periodToBpm(tau2) : 0.0);
 
        // Confidence is the coherence of the delivered grid, reduced only if
        // the alternative explained the signal nearly as well. Coherence alone
        // answers "does this grid fit"; it cannot answer "and is it the only
        // grid that does", and on material that genuinely carries two pulses
        // -- three against two, most obviously -- a tracker that reports one
        // of them at full confidence is asserting something it has no evidence
        // for. See kAmbiguityFloor for why the reduction has a floor under it
        // instead of being applied in proportion.
        const double base = coherence((slope > 0.0) ? slope / static_cast<double>(hop_)
                                                    : tau1);
        out.confidence = static_cast<T>(std::clamp(base * ambiguityDiscount(secondaryShare_),
                                                   0.0, 1.0));
 
        resetState();
        return out;
    }
    // -- Constants -----------------------------------------------------------
 
    static constexpr double kBankMinBpm = 20.0;
    static constexpr double kBankMaxBpm = 480.0;
    static constexpr double kHarmonicReach = 5.0;
    static constexpr int kBinsPerOctave = 128; 
    static constexpr double kDefaultMinBpm = 40.0;
    static constexpr double kDefaultMaxBpm = 240.0;
 
    static constexpr int kNumHarmonics = 3;
    static constexpr double kHarmonicDivisors[kNumHarmonics] = { 2.0, 3.0, 5.0 };
 
    static constexpr double kIntegrationBeats = 8.0;
 
    static constexpr double kPriorCentreSeconds = 0.5;
    static constexpr double kPriorWidthOctaves = 1.4;
 
    static constexpr double kTactusFastOctaves = 0.45;
    static constexpr double kTactusSlowOctaves = kPriorWidthOctaves;
 
    static constexpr int kLevelCandidates = 3;
    static constexpr int kLevelMiddle = 1;
    static constexpr size_t kLevelRawFeatures = 18;
    static constexpr size_t kLevelShareAtPeriod = 1;   
    static constexpr size_t kLevelBalancedShareAtPeriod = 2; 
    static constexpr size_t kLevelShareAtDouble = 7;   
    static constexpr size_t kLevelShareAtTriple = 14;  
    static constexpr std::array<double, kLevelCandidates> kLevelMultiples = { 0.5, 1.0, 2.0 };
    static constexpr std::array<double, kLevelCandidates> kLevelBias = {
        -0.975336646, +0.092188351, +0.883148294 };
    static constexpr std::array<double, 2 * kLevelRawFeatures + 1> kLevelWeights = {
        // Each reading on its own:
        -0.542395096, -0.584902910, +0.498236961,
        +2.173179030, +0.290865194, +0.064019696,
        -0.056958253, -1.551139582, +2.419398334,
        +0.999335380, +1.771192530, +1.715518634,
        -0.269329648, -0.242521186, +0.512883972,
        +2.417326254, -4.177049390, -0.658688841,
        // The same, relative to the proposed reading:
        -0.740402303, -0.614477620, +0.536878776,
        +2.283943101, +0.321349606, +0.070567309,
        -0.060601959, -1.661774996, +2.331723324,
        +1.031193638, +1.895745003, +1.908437878,
        -0.571387423, -0.456177537, +0.860604058,
        +4.361899112, -4.360371230, -0.993521862,
        // Log tempo squared: the tapping preference.
        -2.979577030 };
    static constexpr std::array<double, 4> kTernaryMultiples = {
        1.0 / 3.0, 2.0 / 3.0, 1.5, 3.0 };
    static constexpr std::array<double, 4> kTernaryBias = {
        0.53602663631392855, -0.44218443943636193, -0.98844100123510192,
        -4.3780218173667675 };
    static constexpr std::array<double, 2 * kLevelRawFeatures + 1> kTernaryWeights = {
        0.135530296867787, 1.516476106187614, 0.50512741958750595,
        -3.2226909770906635, 0.76647827698580318, -0.5179245237669392,
        1.354793380662362, -2.8513739655680381, -2.7957191619709669,
        -0.93155463697868401, -0.59124401355810241, 0.57664680014957181,
        -0.83331592757939144, -0.92868430654211009, 0.72147135283011565,
        1.342529710852844, -0.78142066296951229, 2.2830394447846869,
        -0.71810354059356818, 0.89986068541218478, 0.41744186276874551,
        1.3685809837441467, 1.8839328472219159, -0.48244306625392142,
        0.89516613728223393, -0.60386251174497119, 1.8096407974850166,
        2.9213215771564656, 4.4113655210180074, 3.2066836727020585,
        2.0230309593047608, 2.7743701263252531, -2.8072203231223085,
        -0.60235479209693188, 0.10118126378289685, -2.5817564678948961,
        -0.67204847147376856 };
 
    static constexpr double kLevelContrastFloor = 0.1;
    static constexpr double kLevelWindowBeats = 8.0;
 
    static constexpr double kDefaultTightness = 25.0;
    static constexpr double kMinTightness = 1e-3;
 
    static constexpr double kMinFundamentalCoherence = 0.1;
 
    static constexpr double kAmbiguityFloor = 0.5;
 
    static constexpr double kCandidateSeparationOctaves = 0.25;
 
    static constexpr double kCandidatePeakFloor = 0.25;
 
    static constexpr double kStrictTightnessFactor = 4.0;
 
    static constexpr double kMovingGridMargin = 1.5;
 
    static constexpr int kLocalPeriodSpan = kBinsPerOctave / 6;
 
    static constexpr double kMinAnalysisBeats = 4.0;
 
    static constexpr double kCoherenceWindowBeats = 8.0;
 
    static constexpr double kBaselineSeconds = 1.0;
 
    static constexpr double kLocalScoreSigmaPeriods = 1.0 / 32.0;
 
    static constexpr double kSearchLowFactor = 0.5;
    static constexpr double kSearchHighFactor = 2.0;
 
    static constexpr double kSnapFraction = 0.125;
 
    static constexpr double kTrimFraction = 0.5;
 
    static_assert(std::atomic<std::uint64_t>::is_always_lock_free,
                  "audio-thread stores must not lock");
    static_assert(std::atomic<std::uint32_t>::is_always_lock_free,
                  "audio-thread stores must not lock");
    static_assert(std::atomic<std::int64_t>::is_always_lock_free,
                  "audio-thread stores must not lock");
    static_assert(std::atomic<bool>::is_always_lock_free,
                  "audio-thread stores must not lock");
    static_assert(std::atomic<double>::is_always_lock_free,
                  "audio-thread stores must not lock");
 
    // -- Packing -------------------------------------------------------------
 
    [[nodiscard]] static std::uint64_t packTempo(float bpm, float conf) noexcept
    {
        std::uint32_t a = 0, b = 0;
        std::memcpy(&a, &bpm, sizeof(a));
        std::memcpy(&b, &conf, sizeof(b));
        return (static_cast<std::uint64_t>(a) << 32) | static_cast<std::uint64_t>(b);
    }
 
    static void unpackTempo(std::uint64_t v, float& bpm, float& conf) noexcept
    {
        const std::uint32_t a = static_cast<std::uint32_t>(v >> 32);
        const std::uint32_t b = static_cast<std::uint32_t>(v & 0xFFFFFFFFu);
        std::memcpy(&bpm, &a, sizeof(bpm));
        std::memcpy(&conf, &b, sizeof(conf));
    }
 
    [[nodiscard]] static std::uint32_t packRange(int lo, int hi) noexcept
    {
        return (static_cast<std::uint32_t>(lo & 0xFFFF) << 16)
             | static_cast<std::uint32_t>(hi & 0xFFFF);
    }
 
    static void unpackRange(std::uint32_t v, int& lo, int& hi) noexcept
    {
        lo = static_cast<int>((v >> 16) & 0xFFFFu);
        hi = static_cast<int>(v & 0xFFFFu);
    }
 
    // -- Bank ----------------------------------------------------------------
 
    [[nodiscard]] double periodToBpm(double periodFrames) const noexcept
    {
        return (periodFrames > 0.0) ? 60.0 * frameRate_ / periodFrames : 0.0;
    }
 
    [[nodiscard]] double bpmToPeriod(double bpm) const noexcept
    {
        return (bpm > 0.0) ? 60.0 * frameRate_ / bpm : 0.0;
    }
 
    [[nodiscard]] int bpmToIndex(double bpm) const noexcept
    {
        const double p = bpmToPeriod(bpm);
        if (!(p > 0.0) || bankSize_ < 1) return 0;
        const double x = std::log2(p / bankPeriod_[0])
                       * static_cast<double>(kBinsPerOctave);
        return std::clamp(static_cast<int>(std::lround(x)), 0, bankSize_ - 1);
    }
 
    void buildBank()
    {
        const double pMin = bpmToPeriod(kBankMaxBpm);
        const double pMax = bpmToPeriod(kBankMinBpm) * kHarmonicReach;
        const double octaves = std::log2(pMax / pMin);
        bankSize_ = static_cast<int>(std::lround(octaves
                                                 * static_cast<double>(kBinsPerOctave))) + 1;
        bankSize_ = std::max(bankSize_, 2 * kBinsPerOctave + 2);
 
        bankPeriod_.assign(static_cast<size_t>(bankSize_), 0.0);
        bankCosCoef_.assign(static_cast<size_t>(bankSize_), 0.0);
        bankSinCoef_.assign(static_cast<size_t>(bankSize_), 0.0);
        bankDecay_.assign(static_cast<size_t>(bankSize_), 0.0);
        bankPrior_.assign(static_cast<size_t>(bankSize_), 0.0);
        zRe_.assign(static_cast<size_t>(bankSize_), 0.0);
        zIm_.assign(static_cast<size_t>(bankSize_), 0.0);
        zMass_.assign(static_cast<size_t>(bankSize_), 0.0);
        bankCoherence_.assign(static_cast<size_t>(bankSize_), 0.0);
 
        const double centre = kPriorCentreSeconds * frameRate_;
        for (int i = 0; i < bankSize_; ++i)
        {
            const double p = pMin * std::pow(2.0, static_cast<double>(i)
                                                  / static_cast<double>(kBinsPerOctave));
            bankPeriod_[static_cast<size_t>(i)] = p;
 
            const double w = 2.0 * static_cast<double>(pi<double>) / p;
            const double lambda = std::exp(-1.0 / (kIntegrationBeats * p));
            bankDecay_[static_cast<size_t>(i)] = lambda;
            bankCosCoef_[static_cast<size_t>(i)] = lambda * std::cos(w);
            bankSinCoef_[static_cast<size_t>(i)] = lambda * std::sin(w);
 
            const double lg = std::log2(p / centre) / kPriorWidthOctaves;
            bankPrior_[static_cast<size_t>(i)] = std::exp(-0.5 * lg * lg);
        }
 
        // Index distance to each subharmonic. The bank is geometric, so a
        // ratio in period is a constant offset in index; rounding it to the
        // nearest bin is worth nothing at all here, because a resonator with a
        // memory of eight periods is a hundred times broader than one bin.
        for (int k = 0; k < kNumHarmonics; ++k)
            harmonicOffset_[k] = static_cast<int>(std::lround(
                std::log2(kHarmonicDivisors[k]) * static_cast<double>(kBinsPerOctave)));
 
        // The tactus preference of each subharmonic against its own candidate,
        // tabulated here so the per-frame score stays a few multiplies.
        bankLevelRatio_.assign(static_cast<size_t>(kNumHarmonics)
                                   * static_cast<size_t>(bankSize_), 1.0);
        for (int k = 0; k < kNumHarmonics; ++k)
            for (int i = 0; i < bankSize_; ++i)
                bankLevelRatio_[static_cast<size_t>(k) * static_cast<size_t>(bankSize_)
                                + static_cast<size_t>(i)]
                    = levelRatio(bankPeriod_[static_cast<size_t>(i)], kHarmonicDivisors[k]);
    }
 
    // -- Causal path ---------------------------------------------------------
 
    bool advanceFrame() noexcept
    {
        const typename OnsetDetector<T>::OdfFrame f = onset_.getLastOdfFrame();
        ++frameIndex_;
 
        const double raw = (std::isfinite(f.value) && f.value > T(0))
                               ? static_cast<double>(f.value) : 0.0;
        lastFrameRef_ = f.referenceSample;
 
        // Same conditioning the offline path applies, in the only form a
        // causal path can have it: a trailing mean instead of a centred one.
        // Without it the log-compressed envelope's positive floor counts as
        // onset mass that is spread evenly in phase, so the coherence ratio
        // reads the noise floor of the material rather than its pulse -- on a
        // bed at a twentieth of the click amplitude that alone took the
        // published confidence from 0.98 to 0.25 while the tempo stayed
        // correct to a tenth of a percent, which is an under-reading a caller
        // gating on confidence would act on.
        baseline_ += (raw - baseline_) * baselineCoef_;
        const double o = std::max(0.0, raw - baseline_);
 
        for (int i = 0; i < bankSize_; ++i)
        {
            const size_t u = static_cast<size_t>(i);
            const double cr = bankCosCoef_[u];
            const double ci = bankSinCoef_[u];
            const double re = zRe_[u];
            const double im = zIm_[u];
            zRe_[u] = cr * re - ci * im + o;
            zIm_[u] = cr * im + ci * re;
            zMass_[u] = bankDecay_[u] * zMass_[u] + o;
 
            const double mass = zMass_[u];
            bankCoherence_[u] = (mass > kMassFloor)
                                    ? std::min(1.0, std::sqrt(zRe_[u] * zRe_[u]
                                                              + zIm_[u] * zIm_[u]) / mass)
                                    : 0.0;
        }
 
        int iLo = 0, iHi = 0;
        unpackRange(activeRange_.load(std::memory_order_relaxed), iLo, iHi);
        iLo = std::clamp(iLo, 0, bankSize_ - 1);
        iHi = std::clamp(iHi, iLo, bankSize_ - 1);
 
        int best = -1;
        double bestScore = 0.0;
        for (int i = iLo; i <= iHi; ++i)
        {
            const double s = scoreAt(i);
            if (s > bestScore) { bestScore = s; best = i; }
        }
 
        if (best < 0)
        {
            // Nothing in the searched range explains the envelope yet. Keep
            // the previous published pair rather than publishing a tempo of
            // zero and a confidence of zero, which a caller cannot tell from
            // a real reading of silence.
            return false;
        }
 
        // The runner-up at a different metrical level, and how nearly it
        // explained the envelope as well as the winner did. Same discount the
        // offline path applies, for the same reason and on the same scale: a
        // grid that fits is not the same claim as a grid that is the only one
        // that fits, and a caller gating on one number must be told the
        // difference.
        double runnerUp = 0.0;
        const int apart = std::max(1, static_cast<int>(kCandidateSeparationOctaves
                                                       * kBinsPerOctave));
        for (int i = iLo; i <= iHi; ++i)
        {
            if (std::abs(i - best) < apart) continue;
            runnerUp = std::max(runnerUp, scoreAt(i, false));
        }
        const double bestExplains = scoreAt(best, false);
        const double share = (bestExplains > 0.0)
            ? std::clamp(runnerUp / bestExplains, 0.0, 1.0) : 0.0;
 
        // The level is adjudicated with a preference that changes from bin to
        // bin, so the bin it selects can sit one bin off the coherence peak of
        // the very level it selected. The period and the phase are read from
        // the peak itself, which is where the signal actually is: two steps
        // are enough, because the two functions never differ by more than the
        // preference's slope over one bin, and it costs -0.27% -> -0.14% on
        // the worst reported tempo at the top of the default range.
        for (int guard = 0; guard < 2; ++guard)
        {
            if (best > 0 && bankCoherence_[static_cast<size_t>(best - 1)]
                            > bankCoherence_[static_cast<size_t>(best)]) --best;
            else if (best + 1 < bankSize_
                     && bankCoherence_[static_cast<size_t>(best + 1)]
                        > bankCoherence_[static_cast<size_t>(best)]) ++best;
            else break;
        }
 
        const size_t ub = static_cast<size_t>(best);
        const double period = refinePeriod(best);
        const double conf = std::clamp(bankCoherence_[ub] * ambiguityDiscount(share),
                                       0.0, 1.0);
 
        tempoAndConfidence_.store(packTempo(static_cast<float>(periodToBpm(period)),
                                            static_cast<float>(conf)),
                                  std::memory_order_relaxed);
 
        // Phase. arg(z) advances by one turn per period, so the frames since
        // the last pulse are arg(z) scaled by the period; the turn completing
        // -- that quantity falling back through zero -- is the beat.
        const double w = 2.0 * static_cast<double>(pi<double>) / bankPeriod_[ub];
        double sinceBeat = std::atan2(zIm_[ub], zRe_[ub]) / w;
        const double p = bankPeriod_[ub];
        while (sinceBeat < 0.0) sinceBeat += p;
        while (sinceBeat >= p) sinceBeat -= p;
 
        ++framesSinceBeat_;
        bool beat = false;
        // The turn completes when the frames-since-pulse quantity falls back
        // through zero. The refractory is not a smoother: the winning bin can
        // step between neighbours from one frame to the next, and the two
        // neighbours do not carry exactly the same phase, so without it a
        // one-bin step near a turn reports the same beat twice.
        if (havePhase_ && (prevSinceBeat_ - sinceBeat) > 0.5 * p
            && static_cast<double>(framesSinceBeat_) >= kBeatRefractory * p)
        {
            const int64_t backFrames = static_cast<int64_t>(std::lround(sinceBeat));
            const int64_t ref = lastFrameRef_ - backFrames * static_cast<int64_t>(hop_);
            lastBeatSample_.store(ref, std::memory_order_relaxed);
            framesSinceBeat_ = 0;
            beat = true;
        }
        prevSinceBeat_ = sinceBeat;
        lastBestIndex_ = best;
        havePhase_ = true;
 
        return beat;
    }
 
    [[nodiscard]] double refinePeriod(int i) const noexcept
    {
        if (i <= 0 || i >= bankSize_ - 1) return bankPeriod_[static_cast<size_t>(i)];
        const double a = scoreAt(i - 1, false);
        const double b = scoreAt(i, false);
        const double c = scoreAt(i + 1, false);
        const double den = a - 2.0 * b + c;
        if (!(std::abs(den) > 0.0)) return bankPeriod_[static_cast<size_t>(i)];
        double d = 0.5 * (a - c) / den;
        d = std::clamp(d, -0.5, 0.5);
        return bankPeriod_[static_cast<size_t>(i)]
             * std::pow(2.0, d / static_cast<double>(kBinsPerOctave));
    }
 
    [[nodiscard]] double scoreAt(int i, bool withPreference = true) const noexcept
    {
        const size_t u = static_cast<size_t>(i);
        double s = bankPrior_[u] * bankCoherence_[u];
        for (int k = 0; k < kNumHarmonics; ++k)
        {
            const int sub = i + harmonicOffset_[k];
            if (sub >= bankSize_) continue;
            const double r = withPreference
                ? bankLevelRatio_[static_cast<size_t>(k) * static_cast<size_t>(bankSize_) + u]
                : 1.0;
            s *= (1.0 - std::min(1.0, r * bankCoherence_[static_cast<size_t>(sub)]));
        }
        return s;
    }
 
    [[nodiscard]] static double levelPreference(double periodSeconds) noexcept
    {
        const double lg = std::log2(periodSeconds / kPriorCentreSeconds);
        const double w = (lg >= 0.0) ? kTactusSlowOctaves : kTactusFastOctaves;
        const double x = lg / w;
        return std::exp(-0.5 * x * x);
    }
 
    [[nodiscard]] double levelRatio(double periodFrames, double divisor) const noexcept
    {
        if (!(frameRate_ > 0.0) || !(periodFrames > 0.0)) return 1.0;
        const double p = periodFrames / frameRate_;
        const double here = levelPreference(p);
        if (!(here > 0.0)) return 1.0;
        return std::max(1.0, levelPreference(p * divisor) / here);
    }
 
    // -- Offline stages ------------------------------------------------------
 
    void appendEnvelopeFrame()
    {
        const typename OnsetDetector<T>::OdfFrame f = onset_.getLastOdfFrame();
        const double v = std::isfinite(f.value) ? static_cast<double>(f.value) : 0.0;
        env_.push_back(std::max(0.0, v));
        envRef_.push_back(f.referenceSample);
        for (int g = 0; g < kRegisters; ++g)
        {
            const double r = static_cast<double>(f.registers[static_cast<size_t>(g)]);
            envReg_[static_cast<size_t>(g)].push_back(std::isfinite(r) ? std::max(0.0, r) : 0.0);
        }
    }
 
    void conditionEnvelope()
    {
        const int n = static_cast<int>(env_.size());
        if (n < 1) return;
 
        // Register balance. The onset function is a mean over quarter-tone
        // bands, and the upper registers own most of them: a kick lights a
        // handful of bands below 200 Hz while hats, strums and consonants light
        // a hundred. On a dense master the pulse the kick states then barely
        // shows in the mean, and the backbeat and syllable rhythm decide the
        // metrical level instead. So the level is decided on a second,
        // register-balanced envelope: each register conditioned on its own -
        // floor removed, scaled to unit spread - and the four summed, so that
        // a register carries evidence in proportion to how regular it is, not
        // to how many bands it spans. The grid itself is still laid on the
        // plain envelope, where every onset keeps its own weight. Measured on a limited pop master at 135
        // BPM (kick on every beat, snare on two and four, strummed guitar and
        // a syllabic vocal on top), the coherence at half the tempo falls from
        // 0.88 of the coherence at the beat to 0.38.
        envBal_.clear();
        bool haveRegisters = false;
        for (const auto& r : envReg_)
            if (static_cast<int>(r.size()) == n)
                for (const double v : r)
                    if (v > 0.0) { haveRegisters = true; break; }
        if (haveRegisters)
        {
            // A register counts in proportion to how periodic it is: its
            // normalised autocorrelation peak over the searched lags. A
            // register holding one isolated event - a tone switching on - or
            // unpatterned material has no such peak and adds nothing, where
            // plain unit-spread scaling would have amplified it to equal voice.
            int iLo = 0, iHi = 0;
            unpackRange(activeRange_.load(std::memory_order_relaxed), iLo, iHi);
            const int lagLo = std::max(1, static_cast<int>(std::floor(bankPeriod_[static_cast<size_t>(iLo)])));
            const int lagHi = std::min(n / 2, static_cast<int>(std::ceil(bankPeriod_[static_cast<size_t>(iHi)])));
            envBal_.assign(static_cast<size_t>(n), 0.0);
            for (auto& r : envReg_)
            {
                conditionSeries(r, true);
                double zero = 0.0;
                for (const double v : r) zero += v * v;
                zero /= static_cast<double>(n);
                double peak = 0.0;
                if (zero > 0.0)
                    for (int lag = lagLo; lag <= lagHi; ++lag)
                    {
                        checkOfflineWork();
                        double acc = 0.0;
                        for (int i = 0; i + lag < n; ++i)
                            acc += r[static_cast<size_t>(i)] * r[static_cast<size_t>(i + lag)];
                        peak = std::max(peak, acc / static_cast<double>(n - lag) / zero);
                    }
                const double w = std::clamp(peak, 0.0, 1.0);
                for (int i = 0; i < n; ++i) envBal_[static_cast<size_t>(i)] += w * r[static_cast<size_t>(i)];
            }
            conditionSeries(envBal_, false);   // already floor-free: rescale only
        }
        conditionSeries(env_, true);
    }
 
    void conditionSeries(WorkVector<double>& e, bool removeFloor)
    {
        const int n = static_cast<int>(e.size());
        if (n < 1) return;
 
        const int half = std::max(1, static_cast<int>(std::lround(0.5 * kBaselineSeconds
                                                                  * frameRate_)));
        scratchA_.assign(static_cast<size_t>(n) + 1, 0.0);
        for (int i = 0; i < n; ++i)
            scratchA_[static_cast<size_t>(i) + 1] = scratchA_[static_cast<size_t>(i)]
                                                  + e[static_cast<size_t>(i)];
 
        double sumSq = 0.0;
        for (int i = 0; i < n; ++i)
        {
            if ((i & 255) == 0) checkOfflineWork();
            const int lo = std::max(0, i - half);
            const int hi = std::min(n, i + half + 1);
            const double mean = (scratchA_[static_cast<size_t>(hi)]
                                 - scratchA_[static_cast<size_t>(lo)])
                              / static_cast<double>(hi - lo);
            const double v = removeFloor ? std::max(0.0, e[static_cast<size_t>(i)] - mean)
                                         : e[static_cast<size_t>(i)];
            e[static_cast<size_t>(i)] = v;
            sumSq += v * v;
        }
 
        const double rms = std::sqrt(sumSq / static_cast<double>(n));
        if (rms > 0.0)
        {
            const double g = 1.0 / rms;
            for (int i = 0; i < n; ++i) e[static_cast<size_t>(i)] *= g;
        }
    }
 
    [[nodiscard]] double coherence(double periodFrames) const
    {
        return coherence(env_, periodFrames);
    }
 
    [[nodiscard]] double coherence(const WorkVector<double>& envelope,
                                    double periodFrames) const
    {
        const int n = static_cast<int>(envelope.size());
        if (n < 2 || !(periodFrames > 1.0)) return 0.0;
 
        const int win = std::max(4, static_cast<int>(std::lround(kCoherenceWindowBeats
                                                                 * periodFrames)));
        if (win > n) return coherenceWindow(envelope, 0, n, periodFrames);
 
        const int step = std::max(1, win / 2);
        double acc = 0.0;
        int count = 0;
        for (int start = 0; start + win <= n; start += step)
        {
            acc += coherenceWindow(envelope, start, start + win, periodFrames);
            ++count;
        }
        if (count == 0) return coherenceWindow(envelope, 0, n, periodFrames);
        return acc / static_cast<double>(count);
    }
 
    [[nodiscard]] double coherenceWindow(const WorkVector<double>& envelope,
                                          int from, int to, double periodFrames) const
    {
        const double w = 2.0 * static_cast<double>(pi<double>) / periodFrames;
        double re = 0.0, im = 0.0, mass = 0.0;
        for (int i = from; i < to; ++i)
        {
            if ((i & 4095) == 0) checkOfflineWork();
            const double v = envelope[static_cast<size_t>(i)];
            const double a = w * static_cast<double>(i);
            re += v * std::cos(a);
            im -= v * std::sin(a);
            mass += v;
        }
        if (!(mass > kMassFloor)) return 0.0;
        return std::min(1.0, std::sqrt(re * re + im * im) / mass);
    }
 
    [[nodiscard]] std::array<double, 2> offlineCandidateScores(
        const WorkVector<double>& envelope, double period) const
    {
        const double centre = kPriorCentreSeconds * frameRate_;
        const double lg = std::log2(period / centre) / kPriorWidthOctaves;
        const double base = std::exp(-0.5 * lg * lg) * coherence(envelope, period);
        double adjudicated = base, explains = base;
        for (const double multiple : kHarmonicDivisors)
        {
            const double share = coherence(envelope, multiple * period);
            adjudicated *= (1.0 - std::min(1.0, levelRatio(period, multiple) * share));
            explains *= (1.0 - std::min(1.0, share));
        }
        return { adjudicated, explains };
    }
 
    void pickCandidates(int iLo, int iHi, double& tau1, double& tau2)
    {
        tau1 = 0.0;
        tau2 = 0.0;
        secondaryShare_ = 0.0;
        fundamentalFound_ = false;
 
        const int n = static_cast<int>(env_.size());
        const int lagLo = std::max(1, static_cast<int>(std::floor(
                                          bankPeriod_[static_cast<size_t>(iLo)])));
        const int lagHi = std::min(n / 2, static_cast<int>(std::ceil(
                                              bankPeriod_[static_cast<size_t>(iHi)])));
        if (lagHi <= lagLo) return;
 
        // Score one guard lag on either side so a genuine maximum exactly at
        // the requested range endpoint can be tested against both neighbors.
        // Previously both endpoints were excluded from candidate generation.
        acf_.assign(static_cast<size_t>(lagHi + 2), 0.0);
        const double centre = kPriorCentreSeconds * frameRate_;
        for (int lag = std::max(1, lagLo - 1); lag <= lagHi + 1; ++lag)
        {
            checkOfflineWork();
            double s = 0.0;
            const int m = n - lag;
            for (int i = 0; i < m; ++i)
                s += env_[static_cast<size_t>(i)] * env_[static_cast<size_t>(i + lag)];
            s /= static_cast<double>(m);
            const double lg = std::log2(static_cast<double>(lag) / centre)
                            / kPriorWidthOctaves;
            acf_[static_cast<size_t>(lag)] = s * std::exp(-0.5 * lg * lg);
        }
 
        // Peaks, not merely the argmax: the runner-up must be a peak of its
        // own, or it is the same peak read one lag over.
        candidates_.clear();
        for (int lag = lagLo; lag <= lagHi; ++lag)
        {
            const double v = acf_[static_cast<size_t>(lag)];
            if (v > acf_[static_cast<size_t>(lag - 1)]
                && v >= acf_[static_cast<size_t>(lag + 1)] && v > 0.0)
            {
                // A guard lag helps locate the maximum; it must not introduce
                // an out-of-range metrical hypothesis after interpolation.
                const double refined = refineLag(lag);
                if (refined >= bankPeriod_[static_cast<size_t>(iLo)]
                    && refined <= bankPeriod_[static_cast<size_t>(iHi)])
                    candidates_.push_back(lag);
            }
        }
        if (candidates_.empty()) return;
 
        // A candidate is a lag at which the envelope actually repeats. Local
        // maxima of the correlation's noise floor between its real peaks - a
        // hundredth of their height - are lags at which it does not, and the
        // coherence ranking below cannot tell them apart from real pulses on
        // dense material, where every coherence is small: measured on a
        // limited 92 BPM master, a 141 BPM bump at 0.3% of the tallest peak
        // outranked the 92 BPM peak itself. Only peaks within a fixed ratio
        // of the tallest are ranked.
        {
            double top = 0.0;
            for (const int lag : candidates_) top = std::max(top, acf_[static_cast<size_t>(lag)]);
            size_t w = 0;
            for (const int lag : candidates_)
                if (acf_[static_cast<size_t>(lag)] >= kCandidatePeakFloor * top)
                    candidates_[w++] = lag;
            candidates_.resize(w);
        }
 
        // Rank by the octave discriminator, not by correlation height. Two
        // scores per candidate, for the two questions of scoreAt(): the
        // adjudicated one decides the level, the explanatory one is what the
        // confidence is discounted against.
        candScore_.assign(candidates_.size(), 0.0);
        candExplains_.assign(candidates_.size(), 0.0);
        for (size_t k = 0; k < candidates_.size(); ++k)
        {
            const double tau = refineLag(candidates_[k]);
            const auto scores = offlineCandidateScores(env_, tau);
            candScore_[k] = scores[0];
            candExplains_[k] = scores[1];
        }
 
        size_t k1 = 0;
        for (size_t k = 1; k < candidates_.size(); ++k)
            if (candScore_[k] > candScore_[k1]) k1 = k;
 
        // The discriminator can only choose between candidates that look like
        // a fundamental. When the caller has restricted the range so that the
        // signal's own period falls outside it, every candidate left inside is
        // one of that period's subharmonics -- and a subharmonic is exactly
        // what the coherence term is built to score at nothing, so the ranking
        // would be choosing between residues. That is not a failure of the
        // caller's request: a listener held to 60..100 BPM on a 160 BPM pulse
        // hears 80, and 80 is a peak of the correlation. So when no candidate
        // clears the floor, the correlation decides on its own.
        bool anyFundamental = false;
        for (size_t k = 0; k < candidates_.size(); ++k)
            if (coherence(refineLag(candidates_[k])) >= kMinFundamentalCoherence)
            { anyFundamental = true; break; }
        fundamentalFound_ = anyFundamental;
        if (!anyFundamental)
        {
            k1 = 0;
            for (size_t k = 1; k < candidates_.size(); ++k)
                if (acf_[static_cast<size_t>(candidates_[k])]
                    > acf_[static_cast<size_t>(candidates_[k1])]) k1 = k;
            for (size_t k = 0; k < candidates_.size(); ++k)
                candScore_[k] = acf_[static_cast<size_t>(candidates_[k])];
        }
        tau1 = refineLag(candidates_[k1]);
 
        // The reported alternative and the ambiguity it implies are two
        // different quantities and are taken from the two different scores.
        // What gets REPORTED is the level a listener could plausibly have
        // tapped instead, so it is the best adjudicated alternative. What
        // DISCOUNTS the confidence is the best explanation of the envelope at
        // any other level, tapping preference set aside -- on a polyrhythm the
        // most competitive explanation is the tatum, which is not a level
        // anyone taps, and a confidence blind to it would call the signal
        // unambiguous because the reported alternative happens to be modest.
        int k2 = -1;
        double bestOther = 0.0;
        for (size_t k = 0; k < candidates_.size(); ++k)
        {
            const double sep = std::abs(std::log2(static_cast<double>(candidates_[k])
                                                  / static_cast<double>(candidates_[k1])));
            if (sep < kCandidateSeparationOctaves) continue;
            if (k2 < 0 || candScore_[k] > candScore_[static_cast<size_t>(k2)])
                k2 = static_cast<int>(k);
            bestOther = std::max(bestOther, candExplains_[k]);
        }
        if (k2 >= 0)
        {
            tau2 = refineLag(candidates_[static_cast<size_t>(k2)]);
            const double top = candExplains_[k1];
            secondaryShare_ = (top > 0.0) ? std::clamp(bestOther / top, 0.0, 1.0) : 0.0;
        }
    }
 
    [[nodiscard]] double refineLag(int lag) const noexcept
    {
        if (lag <= 0 || lag + 1 >= static_cast<int>(acf_.size()))
            return static_cast<double>(lag);
        const double a = acf_[static_cast<size_t>(lag - 1)];
        const double b = acf_[static_cast<size_t>(lag)];
        const double c = acf_[static_cast<size_t>(lag + 1)];
        const double den = a - 2.0 * b + c;
        if (!(std::abs(den) > 0.0)) return static_cast<double>(lag);
        const double d = std::clamp(0.5 * (a - c) / den, -0.5, 0.5);
        return static_cast<double>(lag) + d;
    }
 
    [[nodiscard]] std::array<double, kLevelRawFeatures> levelFeatures(
        double period, const WorkVector<double>& balanced) const noexcept
    {
        std::array<double, kLevelRawFeatures> features {};
        size_t k = 0;
        features[k++] = std::log2(periodToBpm(period) / 120.0);
        features[k++] = pulseShare(env_, period);
        features[k++] = pulseShare(balanced, period);
        for (const auto& r : envReg_) features[k++] = pulseShare(r, period);
        features[k++] = pulseShare(env_, 2.0 * period);
        for (const auto& r : envReg_) features[k++] = pulseShare(r, 2.0 * period);
        features[k++] = repetition(env_, period);
        features[k++] = repetition(balanced, period);
        features[k++] = pulseShare(env_, 3.0 * period);
        features[k++] = pulseShare(envReg_[0], 3.0 * period);
        features[k++] = pulseShare(balanced, 2.0 * period);
        features[k++] = pulseShare(balanced, 3.0 * period);
        return features;
    }
 
    void chooseMetricalLevel(int iLo, int iHi, double& tau1, double& tau2)
    {
        const double pMin = bankPeriod_[static_cast<size_t>(iLo)];
        const double pMax = bankPeriod_[static_cast<size_t>(iHi)];
        const WorkVector<double>& balanced = (envBal_.size() == env_.size()) ? envBal_ : env_;
        std::array<std::array<double, kLevelRawFeatures>, kLevelCandidates> raw {};
        for (int c = 0; c < kLevelCandidates; ++c)
            raw[static_cast<size_t>(c)] = levelFeatures(
                tau1 * kLevelMultiples[static_cast<size_t>(c)], balanced);
 
        const auto& proposed = raw[static_cast<size_t>(kLevelMiddle)];
        if (!(proposed[kLevelShareAtPeriod] > 0.0)) return;
        const double contrastFloor = kLevelContrastFloor * proposed[kLevelShareAtPeriod];
        int best = kLevelMiddle;
        if (proposed[kLevelShareAtDouble] >= contrastFloor)
        {
            best = -1;
            double bestScore = 0.0;
            for (int c = 0; c < kLevelCandidates; ++c)
            {
                const double period = tau1 * kLevelMultiples[static_cast<size_t>(c)];
                if (c != kLevelMiddle && (period < pMin || period > pMax)) continue;
                const auto& f = raw[static_cast<size_t>(c)];
                double score = kLevelBias[static_cast<size_t>(c)];
                for (size_t j = 0; j < kLevelRawFeatures; ++j)
                    score += kLevelWeights[j] * f[j]
                           + kLevelWeights[kLevelRawFeatures + j] * (f[j] - proposed[j]);
                score += kLevelWeights[2 * kLevelRawFeatures] * f[0] * f[0];
                if (best < 0 || score > bestScore) { best = c; bestScore = score; }
            }
        }
 
        const double binaryPeriod = tau1 * kLevelMultiples[static_cast<size_t>(best)];
        double chosenPeriod = binaryPeriod;
        bool ternaryChosen = false;
        std::array<double, kTernaryMultiples.size()> supportedPeriods {};
        size_t supportedCount = 0;
        if (std::max(proposed[kLevelShareAtDouble], proposed[kLevelShareAtTriple]) >= contrastFloor)
        {
            const auto& reference = raw[static_cast<size_t>(best)];
            double bestGain = 0.0;
            for (size_t c = 0; c < kTernaryMultiples.size(); ++c)
            {
                const double period = tau1 * kTernaryMultiples[c];
                if (period < pMin || period > pMax) continue;
                const auto f = levelFeatures(period, balanced);
                if (f[kLevelBalancedShareAtPeriod] < kMinFundamentalCoherence) continue;
                supportedPeriods[supportedCount++] = period;
                double gain = kTernaryBias[c];
                for (size_t j = 0; j < kLevelRawFeatures; ++j)
                    gain += kTernaryWeights[j] * f[j]
                          + kTernaryWeights[kLevelRawFeatures + j] * (f[j] - reference[j]);
                gain += kTernaryWeights[2 * kLevelRawFeatures] * f[0] * f[0];
                if (gain > bestGain)
                {
                    bestGain = gain;
                    chosenPeriod = period;
                    ternaryChosen = true;
                }
            }
        }
        if (ternaryChosen)
        {
            tau2 = binaryPeriod;
            tau1 = chosenPeriod;
        }
        else if (best != kLevelMiddle)
        {
            tau2 = tau1;
            tau1 = binaryPeriod;
        }
        if (ternaryChosen || best != kLevelMiddle)
            refreshMetricalAmbiguity(balanced, tau1, tau2,
                                    { supportedPeriods.data(), supportedCount });
    }
 
    void refreshMetricalAmbiguity(const WorkVector<double>& envelope,
                                  double chosen, double displaced,
                                  std::span<const double> addedPeriods)
    {
        const double top = offlineCandidateScores(envelope, chosen)[1];
        double bestOther = 0.0;
        for (size_t k = 0; k < candidates_.size(); ++k)
            if (std::abs(std::log2(refineLag(candidates_[k]) / chosen))
                    >= kCandidateSeparationOctaves)
                bestOther = std::max(bestOther, candExplains_[k]);
        const auto consider = [&](double period) {
            if (period > 0.0 && std::abs(std::log2(period / chosen))
                                  >= kCandidateSeparationOctaves)
                bestOther = std::max(bestOther, offlineCandidateScores(envelope, period)[1]);
        };
        consider(displaced);
        for (const double period : addedPeriods) consider(period);
        secondaryShare_ = top > 0.0 ? std::clamp(bestOther / top, 0.0, 1.0)
                                   : (bestOther > 0.0 ? 1.0 : 0.0);
    }
 
    [[nodiscard]] static double pulseShare(const WorkVector<double>& e,
                                           double periodFrames) noexcept
    {
        const int n = static_cast<int>(e.size());
        const int win = static_cast<int>(kLevelWindowBeats * periodFrames);
        if (win < 4) return 0.0;
        const double w = 2.0 * static_cast<double>(pi<double>) / periodFrames;
        const double dc = std::cos(w), ds = std::sin(w);
        double acc = 0.0, mass = 0.0;
        for (int start = 0; start + win <= n; start += win / 2)
        {
            // The phasor is seeded exactly at each window start and rotated
            // within it; over one window the rotation's rounding is far
            // below anything the features resolve.
            double c = std::cos(w * static_cast<double>(start));
            double s = std::sin(w * static_cast<double>(start));
            double re = 0.0, im = 0.0, m = 0.0;
            for (int i = start; i < start + win; ++i)
            {
                const double v = std::max(0.0, e[static_cast<size_t>(i)]);
                re += v * c;
                im += v * s;
                m += v;
                const double cn = c * dc - s * ds;
                s = s * dc + c * ds;
                c = cn;
            }
            if (m > 0.0) { acc += std::sqrt(re * re + im * im); mass += m; }
        }
        return (mass > 0.0) ? acc / mass : 0.0;
    }
 
    [[nodiscard]] static double repetition(const WorkVector<double>& e,
                                           double periodFrames) noexcept
    {
        const int n = static_cast<int>(e.size());
        const int lag = static_cast<int>(std::lround(periodFrames));
        if (lag <= 0 || lag >= n) return 0.0;
        double s0 = 0.0, s = 0.0;
        for (int i = 0; i < n; ++i) s0 += e[static_cast<size_t>(i)] * e[static_cast<size_t>(i)];
        for (int i = 0; i + lag < n; ++i)
            s += e[static_cast<size_t>(i)] * e[static_cast<size_t>(i + lag)];
        return (s0 > 0.0) ? s / s0 * static_cast<double>(n) / static_cast<double>(n - lag)
                          : 0.0;
    }
 
    void computeLocalPeriods(double globalTau)
    {
        const int n = static_cast<int>(env_.size());
        localPeriod_.assign(static_cast<size_t>(n), globalTau);
        if (n < 8 || bankSize_ < 2) return;
 
        // Confined twice over: to kLocalPeriodSpan either side of the period
        // the global analysis settled on, AND to the range the caller asked to
        // be searched. The first keeps a syncopated part from moving the grid
        // off the beat; the second is the caller's contract, and a local
        // estimate that escaped it would deliver a tempo outside the range
        // that was asked for -- measured at 110 BPM on 160 BPM material
        // searched over 60 to 100 before this clamp existed.
        int aLo = 0, aHi = 0;
        unpackRange(activeRange_.load(std::memory_order_relaxed), aLo, aHi);
        aLo = std::clamp(aLo, 0, bankSize_ - 1);
        aHi = std::clamp(aHi, aLo, bankSize_ - 1);
 
        const int gi = periodToIndex(globalTau);
        int lo = std::clamp(gi - kLocalPeriodSpan, aLo, aHi);
        int hi = std::clamp(gi + kLocalPeriodSpan, aLo, aHi);
        if (hi < lo) std::swap(lo, hi);
        if (hi <= lo) return;
 
        sweepA_.assign(static_cast<size_t>(n), globalTau);
        sweepB_.assign(static_cast<size_t>(n), globalTau);
        sweepPeriods(lo, hi, false, sweepA_);
        sweepPeriods(lo, hi, true, sweepB_);
 
        for (int i = 0; i < n; ++i)
        {
            if ((i & 255) == 0) checkOfflineWork();
            const double a = sweepA_[static_cast<size_t>(i)];
            const double b = sweepB_[static_cast<size_t>(i)];
            localPeriod_[static_cast<size_t>(i)] =
                (a > 0.0 && b > 0.0) ? std::sqrt(a * b) : globalTau;
        }
 
        // Smooth over a couple of beats: a per-frame winner steps by whole
        // bank bins, and a target period that steps is a target period the
        // interval cost can be made to chase.
        const int half = std::max(1, static_cast<int>(std::lround(globalTau)));
        scratchA_.assign(static_cast<size_t>(n) + 1, 0.0);
        for (int i = 0; i < n; ++i)
            scratchA_[static_cast<size_t>(i) + 1] = scratchA_[static_cast<size_t>(i)]
                + std::log(localPeriod_[static_cast<size_t>(i)]);
        for (int i = 0; i < n; ++i)
        {
            if ((i & 255) == 0) checkOfflineWork();
            const int a = std::max(0, i - half);
            const int b = std::min(n, i + half + 1);
            const double m = (scratchA_[static_cast<size_t>(b)]
                              - scratchA_[static_cast<size_t>(a)])
                           / static_cast<double>(b - a);
            localPeriod_[static_cast<size_t>(i)] = std::exp(m);
        }
    }
 
    void sweepPeriods(int lo, int hi, bool reverse, WorkVector<double>& out)
    {
        const int n = static_cast<int>(env_.size());
        std::fill(zRe_.begin(), zRe_.end(), 0.0);
        std::fill(zIm_.begin(), zIm_.end(), 0.0);
        std::fill(zMass_.begin(), zMass_.end(), 0.0);
 
        for (int k = 0; k < n; ++k)
        {
            if ((k & 255) == 0) checkOfflineWork();
            const int i = reverse ? (n - 1 - k) : k;
            const double o = env_[static_cast<size_t>(i)];
 
            int best = -1;
            double bestValue = 0.0;
            for (int b = lo; b <= hi; ++b)
            {
                const size_t u = static_cast<size_t>(b);
                const double re = zRe_[u], im = zIm_[u];
                zRe_[u] = bankCosCoef_[u] * re - bankSinCoef_[u] * im + o;
                zIm_[u] = bankCosCoef_[u] * im + bankSinCoef_[u] * re;
                zMass_[u] = bankDecay_[u] * zMass_[u] + o;
                const double mass = zMass_[u];
                const double c = (mass > kMassFloor)
                    ? std::sqrt(zRe_[u] * zRe_[u] + zIm_[u] * zIm_[u]) / mass : 0.0;
                bankCoherence_[u] = c;
                const double s = bankPrior_[u] * c;
                if (s > bestValue) { bestValue = s; best = b; }
            }
            out[static_cast<size_t>(i)] = (best >= 0)
                ? bankPeriod_[static_cast<size_t>(best)] : 0.0;
        }
    }
 
    [[nodiscard]] static double ambiguityDiscount(double share) noexcept
    {
        if (!(share > kAmbiguityFloor)) return 1.0;
        return std::clamp(1.0 - (share - kAmbiguityFloor) / (1.0 - kAmbiguityFloor),
                          0.0, 1.0);
    }
 
    [[nodiscard]] int periodToIndex(double periodFrames) const noexcept
    {
        if (!(periodFrames > 0.0) || bankSize_ < 1) return 0;
        const double x = std::log2(periodFrames / bankPeriod_[0])
                       * static_cast<double>(kBinsPerOctave);
        return std::clamp(static_cast<int>(std::lround(x)), 0, bankSize_ - 1);
    }
 
    bool buildGrid(double periodFrames, WorkVector<int64_t>& gridOut)
    {
        gridOut.clear();
        WorkVector<int> beatFrames(detail::AccountedAllocator<int>{&allocationAccount_});
        if (!runDp(periodFrames, beatFrames) || beatFrames.size() < 2) return false;
        gridOut.reserve(beatFrames.size());
        for (const int f : beatFrames)
            gridOut.push_back(beatSample(f));
        return true;
    }
 
    bool runDp(double periodFrames, WorkVector<int>& beatFrames)
    {
        beatFrames.clear();
        const int n = static_cast<int>(env_.size());
        const bool varying = (localPeriod_.size() == env_.size());
        const int lo = std::max(1, static_cast<int>(std::lround(kSearchLowFactor
                                                                * periodFrames)));
        const int hi = std::max(lo + 1, static_cast<int>(std::lround(kSearchHighFactor
                                                                     * periodFrames)));
        if (n <= hi + 2) return false;
 
        buildLocalScore(periodFrames);
 
        const double alpha = tightness_.load(std::memory_order_relaxed) * alphaScale_;
        cumScore_.assign(static_cast<size_t>(n), 0.0);
        backlink_.assign(static_cast<size_t>(n), -1);
 
        for (int i = 0; i < n; ++i)
        {
            if ((i & 255) == 0) checkOfflineWork();
            const double target = varying ? localPeriod_[static_cast<size_t>(i)]
                                          : periodFrames;
            const int dLo = varying
                ? std::max(1, static_cast<int>(std::lround(kSearchLowFactor * target)))
                : lo;
            const int dHiBound = varying
                ? std::max(dLo + 1, static_cast<int>(std::lround(kSearchHighFactor
                                                                 * target)))
                : hi;
 
            double best = 0.0;
            int bestJ = -1;
            const int dHi = std::min(dHiBound, i);
            for (int d = dLo; d <= dHi; ++d)
            {
                const int j = i - d;
                const double lr = std::log(static_cast<double>(d) / target);
                const double s = cumScore_[static_cast<size_t>(j)] - alpha * lr * lr;
                if (bestJ < 0 || s > best) { best = s; bestJ = j; }
            }
            cumScore_[static_cast<size_t>(i)] = local_[static_cast<size_t>(i)]
                                              + ((bestJ >= 0) ? best : 0.0);
            backlink_[static_cast<size_t>(i)] = bestJ;
        }
 
        int end = 0;
        for (int i = 1; i < n; ++i)
            if (cumScore_[static_cast<size_t>(i)] > cumScore_[static_cast<size_t>(end)])
                end = i;
 
        for (int i = end; i >= 0; i = backlink_[static_cast<size_t>(i)])
        {
            beatFrames.push_back(i);
            if (backlink_[static_cast<size_t>(i)] < 0) break;
        }
        std::reverse(beatFrames.begin(), beatFrames.end());
 
        trimBeats(beatFrames);
        if (beatFrames.size() < 2) return false;
 
        snapBeats(beatFrames, periodFrames);
        extendGrid(beatFrames);
        snapBeats(beatFrames, periodFrames);
        return beatFrames.size() >= 2;
    }
 
    void snapBeats(WorkVector<int>& beatFrames, double periodFrames) const
    {
        const int n = static_cast<int>(env_.size());
        const int win = std::max(1, static_cast<int>(std::lround(periodFrames
                                                                 * kSnapFraction)));
        int previous = -1;
        for (size_t k = 0; k < beatFrames.size(); ++k)
        {
            const int f = beatFrames[k];
            int best = f;
            double bestValue = local_[static_cast<size_t>(f)];
            const int lo = std::max(previous + 1, f - win);
            const int hi = std::min(n - 1, f + win);
            for (int i = lo; i <= hi; ++i)
            {
                if (local_[static_cast<size_t>(i)] > bestValue)
                {
                    bestValue = local_[static_cast<size_t>(i)];
                    best = i;
                }
            }
            beatFrames[k] = best;
            previous = best;
        }
    }
 
    void extendGrid(WorkVector<int>& beatFrames) const
    {
        if (beatFrames.size() < 3) return;
        const int n = static_cast<int>(env_.size());
 
        WorkVector<int> ibi(detail::AccountedAllocator<int>{&allocationAccount_});
        ibi.reserve(beatFrames.size());
        for (size_t k = 1; k < beatFrames.size(); ++k)
            ibi.push_back(beatFrames[k] - beatFrames[k - 1]);
        std::sort(ibi.begin(), ibi.end());
        const int step = ibi[ibi.size() / 2];
        if (step < 2) return;
 
        double mean = 0.0;
        for (const int f : beatFrames) mean += local_[static_cast<size_t>(f)];
        mean /= static_cast<double>(beatFrames.size());
        const double floorValue = kTrimFraction * mean;
        const int win = std::max(1, static_cast<int>(std::lround(step * kSnapFraction)));
 
        WorkVector<int> front(detail::AccountedAllocator<int>{&allocationAccount_});
        for (int guess = beatFrames.front() - step; guess - win >= 0; guess -= step)
        {
            int best = -1;
            double bestValue = floorValue;
            for (int i = std::max(0, guess - win); i <= std::min(n - 1, guess + win); ++i)
                if (local_[static_cast<size_t>(i)] > bestValue)
                { bestValue = local_[static_cast<size_t>(i)]; best = i; }
            if (best < 0) break;
            front.push_back(best);
            guess = best;
        }
        std::reverse(front.begin(), front.end());
 
        WorkVector<int> back(detail::AccountedAllocator<int>{&allocationAccount_});
        for (int guess = beatFrames.back() + step; guess + win <= n - 1; guess += step)
        {
            int best = -1;
            double bestValue = floorValue;
            for (int i = std::max(0, guess - win); i <= std::min(n - 1, guess + win); ++i)
                if (local_[static_cast<size_t>(i)] > bestValue)
                { bestValue = local_[static_cast<size_t>(i)]; best = i; }
            if (best < 0) break;
            back.push_back(best);
            guess = best;
        }
 
        if (front.empty() && back.empty()) return;
        WorkVector<int> merged(detail::AccountedAllocator<int>{&allocationAccount_});
        merged.reserve(front.size() + beatFrames.size() + back.size());
        merged.insert(merged.end(), front.begin(), front.end());
        merged.insert(merged.end(), beatFrames.begin(), beatFrames.end());
        merged.insert(merged.end(), back.begin(), back.end());
        beatFrames.swap(merged);
    }
 
    void buildLocalScore(double periodFrames)
    {
        const int n = static_cast<int>(env_.size());
        local_.assign(static_cast<size_t>(n), 0.0);
 
        const double sigma = std::max(0.5, kLocalScoreSigmaPeriods * periodFrames);
        const int half = std::max(1, static_cast<int>(std::lround(3.0 * sigma)));
        kernel_.assign(static_cast<size_t>(2 * half + 1), 0.0);
        double ksum = 0.0;
        for (int k = -half; k <= half; ++k)
        {
            const double x = static_cast<double>(k) / sigma;
            const double v = std::exp(-0.5 * x * x);
            kernel_[static_cast<size_t>(k + half)] = v;
            ksum += v;
        }
        if (ksum > 0.0)
            for (double& v : kernel_) v /= ksum;
 
        for (int i = 0; i < n; ++i)
        {
            if ((i & 255) == 0) checkOfflineWork();
            double acc = 0.0;
            const int kLo = std::max(-half, -i);
            const int kHi = std::min(half, n - 1 - i);
            for (int k = kLo; k <= kHi; ++k)
                acc += env_[static_cast<size_t>(i + k)]
                     * kernel_[static_cast<size_t>(k + half)];
            local_[static_cast<size_t>(i)] = acc;
        }
    }
 
    void trimBeats(WorkVector<int>& beatFrames) const
    {
        if (beatFrames.size() < 3) return;
 
        double mean = 0.0;
        for (const int f : beatFrames) mean += local_[static_cast<size_t>(f)];
        mean /= static_cast<double>(beatFrames.size());
        const double floorValue = kTrimFraction * mean;
 
        size_t first = 0;
        while (first + 2 < beatFrames.size()
               && local_[static_cast<size_t>(beatFrames[first])] < floorValue)
            ++first;
 
        size_t last = beatFrames.size() - 1;
        while (last > first + 1
               && local_[static_cast<size_t>(beatFrames[last])] < floorValue)
            --last;
 
        if (first > 0 || last + 1 < beatFrames.size())
            beatFrames = WorkVector<int>(beatFrames.begin()
                                              + static_cast<std::ptrdiff_t>(first),
                                          beatFrames.begin()
                                              + static_cast<std::ptrdiff_t>(last) + 1,
                                          beatFrames.get_allocator());
    }
 
    [[nodiscard]] int64_t beatSample(int f) const noexcept
    {
        const int n = static_cast<int>(env_.size());
        const int64_t base = envRef_[static_cast<size_t>(f)];
        if (f <= 0 || f >= n - 1) return base;
        const double a = env_[static_cast<size_t>(f - 1)];
        const double b = env_[static_cast<size_t>(f)];
        const double c = env_[static_cast<size_t>(f + 1)];
        const double den = a - 2.0 * b + c;
        if (!(den < 0.0)) return base;   // not a peak: nothing to interpolate
        const double d = std::clamp(0.5 * (a - c) / den, -0.5, 0.5);
        return base + static_cast<int64_t>(std::llround(d * static_cast<double>(hop_)));
    }
 
    [[nodiscard]] double gridCoherence(const WorkVector<int64_t>& beats) const noexcept
    {
        if (beats.size() < 2) return 0.0;
        // Judged on the register-balanced envelope when there is one, for the
        // reason the metrical level is: on a dense mix the plain envelope is
        // mostly strums and syllables, and a grid that chases them "explains"
        // it better than the grid on the kick does.
        const WorkVector<double>& e = (envBal_.size() == env_.size()) ? envBal_ : env_;
        double re = 0.0, mass = 0.0;
        size_t k = 0;
        for (size_t i = 0; i < e.size(); ++i)
        {
            const int64_t r = envRef_[i];
            if (r < beats.front()) continue;
            while (k + 1 < beats.size() && beats[k + 1] <= r) ++k;
            if (k + 1 >= beats.size()) break;
            const double span = static_cast<double>(beats[k + 1] - beats[k]);
            if (!(span > 0.0)) continue;
            const double phase = static_cast<double>(r - beats[k]) / span;
            const double o = e[i];
            re += o * std::cos(2.0 * pi<double> * phase);
            mass += o;
        }
        return (mass > kMassFloor) ? re / mass : 0.0;
    }
 
    [[nodiscard]] static double fitBeatSlope(const WorkVector<int64_t>& beats)
    {
        const size_t n = beats.size();
        if (n < 2) return 0.0;
        WorkVector<double> ibi(detail::AccountedAllocator<double>{beats.get_allocator().account()});
        ibi.reserve(n - 1);
        for (size_t k = 1; k < n; ++k) ibi.push_back(static_cast<double>(beats[k] - beats[k - 1]));
        WorkVector<double> sorted = ibi;
        std::nth_element(sorted.begin(), sorted.begin() + static_cast<std::ptrdiff_t>(sorted.size() / 2), sorted.end());
        const double med = sorted[sorted.size() / 2];
        if (!(med > 0.0)) return 0.0;
 
        // Reuse the median scratch for interval counts. The fractional
        // transition ends exactly on a half-beat; its internal indices keep
        // their relative timing. No segment of the grid is discarded.
        for (size_t j = 0; j < ibi.size(); ++j)
            sorted[j] = std::max(1.0, std::round(ibi[j] / med));
        const auto stable = [&](size_t j) {
            const double r = ibi[j] / med;
            return r >= 0.75 && std::abs(r - std::round(r)) < 0.125;
        };
        for (size_t j = 3; j + 3 < ibi.size(); ++j)
        {
            if (stable(j) || !stable(j - 1) || !stable(j - 2) || !stable(j - 3))
                continue;
            for (size_t width = 1; width <= 3; ++width)
            {
                const size_t end = j + width;
                if (end + 3 > ibi.size()
                    || !stable(end) || !stable(end + 1) || !stable(end + 2))
                    continue;
                double total = 0.0;
                for (size_t q = j; q < end; ++q) total += ibi[q] / med;
                const double half = std::floor(total) + 0.5;
                if (std::abs(total - half) >= 0.125) continue;
                for (size_t q = j; q < end; ++q)
                    sorted[q] = (ibi[q] / med) * half / total;
                j = end - 1;
                break;
            }
        }
 
        double sx = 0.0, sy = 0.0, sxx = 0.0, sxy = 0.0, index = 0.0;
        for (size_t k = 0; k < n; ++k)
        {
            if (k > 0) index += sorted[k - 1];
            const double y = static_cast<double>(beats[k]);
            sx += index; sy += y; sxx += index * index; sxy += index * y;
        }
        const double dn = static_cast<double>(n);
        const double den = dn * sxx - sx * sx;
        if (!(std::abs(den) > 0.0)) return 0.0;
        return (dn * sxy - sx * sy) / den;
    }
 
    // -- State ---------------------------------------------------------------
 
    void resetState() noexcept
    {
        onset_.reset();
        std::fill(zRe_.begin(), zRe_.end(), 0.0);
        std::fill(zIm_.begin(), zIm_.end(), 0.0);
        std::fill(zMass_.begin(), zMass_.end(), 0.0);
        std::fill(bankCoherence_.begin(), bankCoherence_.end(), 0.0);
 
        samplesPushed_ = 0;
        frameIndex_ = 0;
        lastFrameRef_ = 0;
        baseline_ = 0.0;
        prevSinceBeat_ = 0.0;
        framesSinceBeat_ = 0;
        lastBestIndex_ = -1;
        havePhase_ = false;
 
        tempoAndConfidence_.store(packTempo(0.0f, 0.0f), std::memory_order_relaxed);
        beatLatched_.store(false, std::memory_order_relaxed);
        lastBeatSample_.store(-1, std::memory_order_relaxed);
    }
 
    static constexpr double kBeatRefractory = 0.5;
 
    static constexpr double kMassFloor = 1e-12;
 
    // -- Members -------------------------------------------------------------
 
    detail::AllocationAccount allocationAccount_;
    void* offlineCheckpointContext_ = nullptr;
    void (*offlineCheckpoint_)(void*) = nullptr;
    void checkOfflineWork() const
    {
        if (offlineCheckpoint_) offlineCheckpoint_(offlineCheckpointContext_);
    }
    OnsetDetector<T> onset_;
 
    double sampleRate_ = 44100.0;
    double frameRate_ = 200.0;
    int hop_ = 221;
 
    // Resonator bank (built at prepare, read on the audio thread).
    int bankSize_ = 0;
    WorkVector<double> bankPeriod_;
    WorkVector<double> bankCosCoef_, bankSinCoef_, bankDecay_, bankPrior_;
    int harmonicOffset_[kNumHarmonics] = { 0, 0, 0 };
    WorkVector<double> zRe_, zIm_, zMass_, bankCoherence_;
    WorkVector<double> bankLevelRatio_;
 
    // Causal stream state (audio thread only).
    int64_t samplesPushed_ = 0;
    int64_t frameIndex_ = 0;
    int64_t lastFrameRef_ = 0;
    double baseline_ = 0.0;
    double baselineCoef_ = 0.005;
    double prevSinceBeat_ = 0.0;
    int64_t framesSinceBeat_ = 0;
    int lastBestIndex_ = -1;
    bool havePhase_ = false;
 
    // Offline scratch (analyze() and the offline session only).
    int64_t offlinePushed_ = 0;
    bool offlineOpen_ = false;
    WorkVector<double> env_;
    WorkVector<int64_t> envRef_;
    static constexpr int kRegisters = OnsetDetector<T>::kNumRegisters;
    static_assert(kRegisters == 4, "the level model (kLevelWeights) reads four register envelopes");
    std::array<WorkVector<double>, kRegisters> envReg_;   
    WorkVector<double> envBal_;   
    double alphaScale_ = 1.0;      
    WorkVector<double> local_, cumScore_, acf_, scratchA_, kernel_;
    WorkVector<double> localPeriod_, sweepA_, sweepB_;
    WorkVector<int> backlink_;
    WorkVector<int> candidates_;
    WorkVector<double> candScore_, candExplains_;
    double secondaryShare_ = 0.0;
    bool fundamentalFound_ = false;   
 
    // Published state.
    std::atomic<bool> prepared_ { false };
    std::atomic<std::uint64_t> tempoAndConfidence_ { 0 };
    std::atomic<std::uint32_t> activeRange_ { 0 };
    std::atomic<double> tightness_ { kDefaultTightness };
    std::atomic<bool> beatLatched_ { false };
    std::atomic<std::int64_t> lastBeatSample_ { -1 };
    std::atomic<std::int64_t> latencySamples_ { 549 };
};
 
} // namespace dspark