PhaseVocoderEngine.h

// DSPark - Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi - MIT License
 
#pragma once
 
#include "../../Core/DspMath.h"
#include "../../Core/FFT.h"
#include "../../Core/WindowFunctions.h"
 
#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <memory>
#include <numbers>
#include <vector>
 
namespace dspark {
namespace detail {
 
template <FloatType T>
class PhaseVocoderEngine
{
public:
    // The staged parameter words are single machine words on every supported
    // target; the audio path must never take a lock to read them.
    static_assert(std::atomic<double>::is_always_lock_free,
                  "PhaseVocoderEngine requires lock-free std::atomic<double>");
    static_assert(std::atomic<unsigned>::is_always_lock_free,
                  "PhaseVocoderEngine requires lock-free std::atomic<unsigned>");
    static_assert(std::atomic<bool>::is_always_lock_free,
                  "PhaseVocoderEngine requires lock-free std::atomic<bool>");
 
    struct Params
    {
        double targetSemitones  = 0.0;    
        bool   transientPreserve = true;  
        bool   formantPreserve   = false; 
        bool   phaseLock         = true;  
        bool   percussiveSplit   = false; 
    };
 
    // -- Lifecycle (setup thread) ---------------------------------------------
 
    bool prepare(double sampleRate, int numChannels, int fftSize,
                 bool resampleCompensation, bool separationSupport = false,
                 bool transientLockedHop = false,
                 bool spectralFluxDetector = false)
    {
        if (!(sampleRate > 0.0) || !std::isfinite(sampleRate) || numChannels < 1
            || (fftSize & (fftSize - 1)) != 0 || fftSize < 256 || fftSize > (1 << 20))
            return false;
 
        prepared_ = false;
 
        sampleRate_  = sampleRate;
        numChannels_ = numChannels;
        fftSize_     = fftSize;
        numBins_     = fftSize_ / 2 + 1;
        synthHop_    = fftSize_ / 4;
        ringMask_    = fftSize_ - 1;
        resampleCompensation_ = resampleCompensation;
 
        accumSize_ = fftSize_ * 4;          // power of two: frame + read margin
        accumMask_ = accumSize_ - 1;
 
        fft_ = std::make_unique<FFTReal<T>>(static_cast<size_t>(fftSize_));
 
        window_.resize(static_cast<size_t>(fftSize_));
        WindowFunctions<T>::hann(window_.data(), fftSize_, true);
        for (auto& w : window_) w = std::sqrt(w);   // sqrt-Hann analysis+synthesis
 
        const int nCh = numChannels_;
        inputRing_.assign(static_cast<size_t>(nCh), {});
        accum_.assign(static_cast<size_t>(nCh), {});
        for (int ch = 0; ch < nCh; ++ch)
        {
            inputRing_[static_cast<size_t>(ch)].assign(static_cast<size_t>(fftSize_), T(0));
            accum_[static_cast<size_t>(ch)].assign(static_cast<size_t>(accumSize_), T(0));
        }
 
        fftIn_.resize(static_cast<size_t>(fftSize_));
        spec_.resize(static_cast<size_t>(fftSize_ + 2));
        fftResult_.resize(static_cast<size_t>(fftSize_));
        prevAnalysis_.resize(static_cast<size_t>(fftSize_ + 2));
        prevSynth_.resize(static_cast<size_t>(fftSize_ + 2));
 
        // Formant preservation scratch (cepstral envelope).
        cepsTime_.resize(static_cast<size_t>(fftSize_));
        cepsSpec_.resize(static_cast<size_t>(fftSize_ + 2));
        envLog_.resize(static_cast<size_t>(numBins_));
        formantGain_.resize(static_cast<size_t>(numBins_));
        mag_.resize(static_cast<size_t>(numBins_));
        prevMag_.resize(static_cast<size_t>(numBins_));
        rotRe_.resize(static_cast<size_t>(numBins_));
        rotIm_.resize(static_cast<size_t>(numBins_));
        peakBin_.resize(static_cast<size_t>(numBins_ / 2 + 2));
 
        // Onset detection front end, sized only for the owner that asked for
        // it. Everything the flux detector needs per frame is allocated
        // here: the filterbank tables, the two band frames it differences,
        // and the flux history the running median is taken over. An owner on
        // the frame-energy test allocates none of it and runs the same
        // detector, over the same magnitudes, as it did before this
        // selection existed.
        fluxDetectorEnabled_ = spectralFluxDetector;
        if (fluxDetectorEnabled_)
        {
            buildOnsetFilterBank();
            bandCur_.assign(static_cast<size_t>(numBands_), T(0));
            bandPrev_.assign(static_cast<size_t>(numBands_), T(0));
            bandMaxPrev_.assign(static_cast<size_t>(numBands_), T(0));
            fluxHist_.assign(static_cast<size_t>(kFluxWindow), 0.0);
            fluxScratch_.assign(static_cast<size_t>(kFluxWindow), 0.0);
            // Un-normalised |X| grows linearly with the frame length, so the
            // band accumulation is referred to a 2048-sample frame before
            // the log compression: the growth cancels and one firing
            // threshold means the same sensitivity at every frame size.
            // fftSize is a power of two, so the factor is exact and
            // introduces no rounding.
            odfScale_ = static_cast<T>(kOdfRefFrame / static_cast<double>(fftSize_));
        }
        else
        {
            numBands_ = 0;
            fbStart_.clear();
            fbOffset_.clear();
            fbCount_.clear();
            fbWeights_.clear();
            bandCur_.clear();
            bandPrev_.clear();
            bandMaxPrev_.clear();
            fluxHist_.clear();
            fluxScratch_.clear();
            odfScale_ = T(1);
        }
        // One lock spans every analysis frame whose window can contain the
        // strike; that count is a property of the window geometry, not of
        // any signal. An owner that did not ask for the locked hop gets a
        // lock length of zero, which is the schedule with no lock in it.
        lockEnabled_ = transientLockedHop;
        lockFrames_ = lockEnabled_ ? (fftSize_ + synthHop_ - 1) / synthHop_ : 0;
 
        // Harmonic/percussive median filtering. Both window lengths are
        // chosen from a PHYSICAL span, so they follow the sample rate and the
        // frame size instead of being fixed bin/frame counts:
        //   - along time, 0.2 s of history, the span over which a musical
        //     partial is steady but a strike is not;
        //   - along frequency, 500 Hz, wide enough to bridge the harmonic
        //     comb of any pitched note in range and narrow enough to leave a
        //     broadband strike untouched.
        // Both are forced odd so the median is a sample of the window rather
        // than an average of two, and both are clamped so extreme frame sizes
        // cannot ask for a window longer than the data.
        separationEnabled_ = separationSupport;
        if (separationEnabled_)
        {
            const double hopSeconds = static_cast<double>(synthHop_) / sampleRate_;
            const double binHz = sampleRate_ / static_cast<double>(fftSize_);
            timeMedian_ = makeOdd(static_cast<int>(std::lround(0.2 / hopSeconds)), 3, 31);
            freqMedian_ = makeOdd(static_cast<int>(std::lround(500.0 / binHz)), 3, 63);
            freqMedian_ = std::min(freqMedian_, makeOdd(numBins_, 3, 63));
 
            magHist_.assign(static_cast<size_t>(timeMedian_) * static_cast<size_t>(numBins_), T(0));
            medianScratch_.resize(static_cast<size_t>(std::max(timeMedian_, freqMedian_)));
            maskHarm_.resize(static_cast<size_t>(numBins_));
            maskPerc_.resize(static_cast<size_t>(numBins_));
            specRaw_.resize(static_cast<size_t>(fftSize_ + 2));
        }
        else
        {
            timeMedian_ = 0;
            freqMedian_ = 0;
            magHist_.clear();
            medianScratch_.clear();
            maskHarm_.clear();
            maskPerc_.clear();
            specRaw_.clear();
        }
 
        prepared_ = true;
        return true;
    }
 
    void reset() noexcept
    {
        if (!prepared_) return;
        for (auto& r : inputRing_) std::fill(r.begin(), r.end(), T(0));
        for (auto& a : accum_)     std::fill(a.begin(), a.end(), T(0));
        std::fill(prevAnalysis_.begin(), prevAnalysis_.end(), T(0));
        std::fill(prevSynth_.begin(), prevSynth_.end(), T(0));
        std::fill(prevMag_.begin(), prevMag_.end(), T(0));
        std::fill(magHist_.begin(), magHist_.end(), T(0));
        histPos_ = 0;
        histFilled_ = 0;
 
        std::fill(bandPrev_.begin(), bandPrev_.end(), T(0));
        std::fill(bandMaxPrev_.begin(), bandMaxPrev_.end(), T(0));
        std::fill(fluxHist_.begin(), fluxHist_.end(), 0.0);
        fluxPos_ = 0;
        fluxFilled_ = 0;
        onsetEnv_ = 0.0;
        lockLeft_ = 0;
        debt_ = 0.0;
 
        inputPos_ = 0;
        writeHead_ = static_cast<int64_t>(accumSize_);       // keep indices positive
 
        adoptParamsIfDirty();
        stActive_ = std::clamp(targetSemitones_, -12.0, 12.0);
        ratioActive_ = std::exp2(stActive_ / 12.0);
        hopCarry_ = 0.0;
        inputSinceHop_ = 0;
        analysisHop_ = computeAnalysisHop();
        firstFrame_ = true;
    }
 
    // -- Parameters (control thread) ------------------------------------------
 
    void publishParams(const Params& p) noexcept
    {
        seq_.fetch_add(1, std::memory_order_acq_rel);   // -> odd
        // Release fence: pairs with the reader's acquire fence via
        // fence-fence synchronization ([atomics.fences]/2). Without it the
        // relaxed data stores below are not ordered against the counter and
        // a torn copy could pass validation on weakly ordered targets.
        std::atomic_thread_fence(std::memory_order_release);
        stgSemitones_.store(p.targetSemitones, std::memory_order_relaxed);
        stgTransient_.store(p.transientPreserve, std::memory_order_relaxed);
        stgFormant_.store(p.formantPreserve, std::memory_order_relaxed);
        stgPhaseLock_.store(p.phaseLock, std::memory_order_relaxed);
        stgSplit_.store(p.percussiveSplit, std::memory_order_relaxed);
        seq_.fetch_add(1, std::memory_order_release);   // -> even
        dirty_.store(true, std::memory_order_release);
    }
 
    // -- Streaming (audio thread, stream owner) -------------------------------
 
    [[nodiscard]] int samplesToNextHop() const noexcept
    {
        return analysisHop_ - inputSinceHop_;
    }
 
    void pushInput(int ch, const T* src, int count) noexcept
    {
        auto& ring = inputRing_[static_cast<size_t>(ch)];
        int wp = inputPos_;
        for (int k = 0; k < count; ++k)
        {
            ring[static_cast<size_t>(wp)] = src[k];
            wp = (wp + 1) & ringMask_;
        }
    }
 
    void commitInput(int count, int numActiveChannels) noexcept
    {
        inputPos_ = (inputPos_ + count) & ringMask_;
        inputSinceHop_ += count;
        if (inputSinceHop_ >= analysisHop_)
        {
            inputSinceHop_ = 0;
            processHop(numActiveChannels);
        }
    }
 
    [[nodiscard]] double activeRatio() const noexcept { return ratioActive_; }
 
    [[nodiscard]] Params activeParams() const noexcept
    {
        Params p;
        p.targetSemitones   = targetSemitones_;
        p.transientPreserve = transientOn_;
        p.formantPreserve   = formantOn_;
        p.phaseLock         = phaseLockOn_;
        p.percussiveSplit   = splitOn_;
        return p;
    }
 
    [[nodiscard]] const T* olaData(int ch) const noexcept
    {
        return accum_[static_cast<size_t>(ch)].data();
    }
 
    [[nodiscard]] int64_t olaMask() const noexcept
    {
        return static_cast<int64_t>(accumMask_);
    }
 
    [[nodiscard]] int olaSize() const noexcept { return accumSize_; }
 
    [[nodiscard]] int64_t writeHead() const noexcept { return writeHead_; }
 
    [[nodiscard]] int fftSize() const noexcept { return prepared_ ? fftSize_ : 0; }
 
    [[nodiscard]] int synthHop() const noexcept { return synthHop_; }
 
private:
    static constexpr double kTwoPi = 2.0 * std::numbers::pi;
 
    static constexpr int kSeqlockMaxAttempts = 3;
 
    static constexpr double kSeparationFactor = 2.0;
 
    static constexpr double kOnsetFMin = 27.5;
    static constexpr double kOnsetFMax = 16000.0;
    static constexpr int kOnsetBandsPerOctave = 24;
    static constexpr double kOdfRefFrame = 2048.0;
    static constexpr int kFluxWindow = 12;
    static constexpr double kFluxDelta = 0.02;
 
    void adoptParamsIfDirty() noexcept
    {
        if (!dirty_.exchange(false, std::memory_order_acquire)) return;
 
        // Bounded seqlock read (kSeqlockMaxAttempts): never adopts a torn set
        // (the accept test is the unbounded form's) and on give-up keeps the
        // previously adopted private copy.
        bool adopted = false;
        double st = 0.0;
        bool tr = true, fo = false, pl = true, sp = false;
        for (int attempt = 0; attempt < kSeqlockMaxAttempts; ++attempt)
        {
            const unsigned s0 = seq_.load(std::memory_order_acquire);
            if ((s0 & 1u) != 0u) continue;   // writer mid-publish: do not copy
            st = stgSemitones_.load(std::memory_order_relaxed);
            tr = stgTransient_.load(std::memory_order_relaxed);
            fo = stgFormant_.load(std::memory_order_relaxed);
            pl = stgPhaseLock_.load(std::memory_order_relaxed);
            sp = stgSplit_.load(std::memory_order_relaxed);
            // The fence orders the copy above BEFORE the re-read below and
            // pairs with the writer's release fence ([atomics.fences]/2); a
            // plain load-acquire orders only LATER accesses, so without it
            // the copy could sink below the re-read on weakly ordered CPUs.
            std::atomic_thread_fence(std::memory_order_acquire);
            if (s0 == seq_.load(std::memory_order_relaxed)) { adopted = true; break; }
        }
        if (!adopted)
        {
            // Adopt nothing: the active copies do not move. Re-arm and pick
            // the update up at a later hop.
            dirty_.store(true, std::memory_order_release);
            return;
        }
        // Commit-on-validate: the active copies are written only here, after
        // the copy validated, so a give-up can never leave a torn set live.
        targetSemitones_ = st;
        transientOn_     = tr;
        formantOn_       = fo;
        phaseLockOn_     = pl;
        splitOn_         = sp && separationEnabled_;
    }
 
    [[nodiscard]] static double princArg(double x) noexcept
    {
        return x - kTwoPi * std::round(x / kTwoPi);
    }
 
    [[nodiscard]] int computeAnalysisHop() noexcept
    {
        double ideal;
        if (lockLeft_ > 0)
        {
            --lockLeft_;
            debt_ += static_cast<double>(synthHop_) / ratioActive_
                   - static_cast<double>(synthHop_);
            ideal = static_cast<double>(synthHop_) + hopCarry_;
        }
        else
        {
            const double cap = 0.5 * static_cast<double>(synthHop_);
            const double repay = std::clamp(debt_, -cap, cap);
            debt_ -= repay;
            ideal = static_cast<double>(synthHop_) / ratioActive_ + repay + hopCarry_;
        }
        int hop = static_cast<int>(ideal);
        hop = std::clamp(hop, 1, fftSize_);
        hopCarry_ = ideal - static_cast<double>(hop);
        return hop;
    }
 
    void armTransientLock(bool transient) noexcept
    {
        if (!lockEnabled_) return;
        if (!transient) return;
        if (firstFrame_) return;
        if (lockLeft_ > 0) return;
        if (std::abs(debt_) >= 1.0) return;
        lockLeft_ = lockFrames_;
    }
 
    void processHop(int nCh) noexcept
    {
        // Adopt pending control publications at the hop boundary (cold by
        // construction: it runs only when the control thread published; the
        // audio thread never publishes to itself through the staged channel).
        adoptParamsIfDirty();
 
        // --- glide the active ratio toward the target (max 0.5 st per hop) ----
        const double stTarget = std::clamp(targetSemitones_, -12.0, 12.0);
        stActive_ += std::clamp(stTarget - stActive_, -0.5, 0.5);
        ratioActive_ = std::exp2(stActive_ / 12.0);
 
        // --- reference-channel analysis ---------------------------------------
        analyzeChannel(0);
 
        const bool transient = detectTransient();
        armTransientLock(transient);
        buildRotations(transient);
        // The magnitude history is kept whenever the owner allocated it, not
        // only while the split is engaged: turning the split on mid-stream
        // must find a full window behind it, not an empty one.
        if (separationEnabled_)
        {
            pushMagnitudeHistory();
            if (splitOn_) updateSeparationMasks();
        }
 
        // Snapshot the reference analysis spectrum for the next frame's
        // heterodyned phase increments (after decisions, before modification).
        std::copy(spec_.begin(), spec_.end(), prevAnalysis_.begin());
        std::copy(mag_.begin(), mag_.end(), prevMag_.begin());
 
        // --- synthesis: reference channel first, then the rest -----------------
        synthesizeChannel(0, true);
        for (int ch = 1; ch < nCh; ++ch)
        {
            analyzeChannel(ch);
            synthesizeChannel(ch, false);
        }
 
        analysisHop_ = computeAnalysisHop();
        writeHead_ += synthHop_;
        firstFrame_ = false;
    }
 
    void analyzeChannel(int ch) noexcept
    {
        const auto& ring = inputRing_[static_cast<size_t>(ch)];
        const int readPos = inputPos_;   // oldest sample (ring size == fftSize)
        for (int k = 0; k < fftSize_; ++k)
        {
            const int idx = (readPos + k) & ringMask_;
            fftIn_[static_cast<size_t>(k)] = ring[static_cast<size_t>(idx)]
                                           * window_[static_cast<size_t>(k)];
        }
        fft_->forward(fftIn_.data(), spec_.data());
 
        if (ch == 0)
        {
            for (int k = 0; k < numBins_; ++k)
            {
                const T re = spec_[static_cast<size_t>(2 * k)];
                const T im = spec_[static_cast<size_t>(2 * k + 1)];
                mag_[static_cast<size_t>(k)] = std::sqrt(re * re + im * im);
            }
        }
    }
 
    [[nodiscard]] bool detectTransient() noexcept
    {
        return fluxDetectorEnabled_ ? detectTransientByFlux()
                                    : detectTransientByEnergy();
    }
 
    [[nodiscard]] bool detectTransientByEnergy() noexcept
    {
        double energy = 0.0;
        for (int k = 0; k < numBins_; ++k)
        {
            const double m = static_cast<double>(mag_[static_cast<size_t>(k)]);
            energy += m * m;
        }
        const double prevEnv = onsetEnv_;
        onsetEnv_ = std::max(energy, onsetEnv_ * 0.7);
        if (!transientOn_)
            return firstFrame_;
        return firstFrame_ || (energy > 4.0 * prevEnv && energy > 1e-12);
    }
 
    [[nodiscard]] bool detectTransientByFlux() noexcept
    {
        const double flux = computeSpectralFlux();
 
        // Median of the trailing window, taken before this frame is pushed:
        // a frame is compared with the material that preceded it, never with
        // itself.
        bool rise = false;
        if (fluxFilled_ >= kFluxWindow)
        {
            std::copy(fluxHist_.begin(), fluxHist_.end(), fluxScratch_.begin());
            const auto mid = fluxScratch_.begin() + kFluxWindow / 2;
            std::nth_element(fluxScratch_.begin(), mid, fluxScratch_.end());
            rise = (flux - *mid) >= kFluxDelta;
        }
 
        fluxHist_[static_cast<size_t>(fluxPos_)] = flux;
        fluxPos_ = (fluxPos_ + 1) % kFluxWindow;
        if (fluxFilled_ < kFluxWindow) ++fluxFilled_;
 
        if (!transientOn_)
            return firstFrame_;
        return firstFrame_ || rise;
    }
 
    [[nodiscard]] double computeSpectralFlux() noexcept
    {
        filterLogBands();
 
        double flux = 0.0;
        for (int b = 0; b < numBands_; ++b)
        {
            const double d = static_cast<double>(bandCur_[static_cast<size_t>(b)])
                           - static_cast<double>(bandMaxPrev_[static_cast<size_t>(b)]);
            if (d > 0.0) flux += d;
        }
        flux /= static_cast<double>(std::max(1, numBands_));
 
        // Rotate: previous <- current, and rebuild the maximum-filtered
        // reference the next frame will be differenced against.
        bandPrev_ = bandCur_;
        maxFilterBands();
        return flux;
    }
 
    void filterLogBands() noexcept
    {
        for (int b = 0; b < numBands_; ++b)
        {
            const int start = fbStart_[static_cast<size_t>(b)];
            const int off   = fbOffset_[static_cast<size_t>(b)];
            const int cnt   = fbCount_[static_cast<size_t>(b)];
            T acc = T(0);
            for (int i = 0; i < cnt; ++i)
            {
                const int k = start + i;
                if (k >= 0 && k < numBins_)
                    acc += mag_[static_cast<size_t>(k)]
                         * fbWeights_[static_cast<size_t>(off + i)];
            }
            bandCur_[static_cast<size_t>(b)] = std::log10(acc * odfScale_ + T(1));
        }
    }
 
    void maxFilterBands() noexcept
    {
        for (int b = 0; b < numBands_; ++b)
        {
            T mx = bandPrev_[static_cast<size_t>(b)];
            if (b > 0) mx = std::max(mx, bandPrev_[static_cast<size_t>(b - 1)]);
            if (b + 1 < numBands_) mx = std::max(mx, bandPrev_[static_cast<size_t>(b + 1)]);
            bandMaxPrev_[static_cast<size_t>(b)] = mx;
        }
    }
 
    void buildOnsetFilterBank()
    {
        fbStart_.clear();
        fbOffset_.clear();
        fbCount_.clear();
        fbWeights_.clear();
        numBands_ = 0;
 
        const double binHz = sampleRate_ / static_cast<double>(fftSize_);
        const double fMax = std::min(kOnsetFMax, sampleRate_ * 0.5 * 0.999);
 
        std::vector<int> centres;
        for (int i = 0; ; ++i)
        {
            const double f = kOnsetFMin * std::pow(2.0, static_cast<double>(i)
                                 / static_cast<double>(kOnsetBandsPerOctave));
            if (f > fMax) break;
            int bin = static_cast<int>(std::lround(f / binHz));
            bin = std::clamp(bin, 0, numBins_ - 1);
            if (centres.empty() || bin > centres.back())
                centres.push_back(bin);
        }
 
        for (size_t j = 1; j + 1 < centres.size(); ++j)
        {
            const int lo = centres[j - 1];
            const int ce = centres[j];
            const int hi = centres[j + 1];
            if (!(lo < ce && ce < hi)) continue;
 
            fbStart_.push_back(lo);
            fbOffset_.push_back(static_cast<int>(fbWeights_.size()));
            for (int k = lo; k <= hi; ++k)
            {
                const T wv = (k <= ce)
                    ? static_cast<T>(static_cast<double>(k - lo)
                                     / static_cast<double>(ce - lo))
                    : static_cast<T>(static_cast<double>(hi - k)
                                     / static_cast<double>(hi - ce));
                fbWeights_.push_back(wv);
            }
            ++numBands_;
        }
        for (int b = 0; b < numBands_; ++b)
        {
            const int off = fbOffset_[static_cast<size_t>(b)];
            const int nextOff = (b + 1 < numBands_)
                                ? fbOffset_[static_cast<size_t>(b + 1)]
                                : static_cast<int>(fbWeights_.size());
            fbCount_.push_back(nextOff - off);
        }
        // A frame size too small to carry a single triangle leaves the flux
        // at zero, which never fires: the detector goes quiet rather than
        // reading out of an empty table.
        if (numBands_ < 1) numBands_ = 0;
    }
 
    void buildRotations(bool transient) noexcept
    {
        if (transient)
        {
            std::fill(rotRe_.begin(), rotRe_.end(), T(1));
            std::fill(rotIm_.begin(), rotIm_.end(), T(0));
            return;
        }
 
        if (!phaseLockOn_)
        {
            buildPerBinRotations();
            return;
        }
 
        // --- find spectral peaks (local maxima over +-2 bins, above floor) -----
        T maxMag = T(0);
        for (int k = 0; k < numBins_; ++k)
            maxMag = std::max(maxMag, mag_[static_cast<size_t>(k)]);
 
        numPeaks_ = 0;
        if (maxMag > T(1e-9))
        {
            const T floorMag = maxMag * T(1e-4);   // -80 dB relative floor
            // Candidates stop at numBins-3 so the +-2-bin neighbour tests
            // stay inside mag_ (the old bound of numBins-2 read mag_[k+2]
            // one float past the end whenever the bin below Nyquist was a
            // candidate above the floor - broadband material at small frame
            // sizes reached it routinely).
            const int last = numBins_ - 3;
            for (int k = 2; k <= last; ++k)
            {
                const T m = mag_[static_cast<size_t>(k)];
                if (m < floorMag) continue;
                if (m > mag_[static_cast<size_t>(k - 1)] && m >= mag_[static_cast<size_t>(k + 1)]
                    && m > mag_[static_cast<size_t>(k - 2)] && m >= mag_[static_cast<size_t>(k + 2)])
                {
                    peakBin_[static_cast<size_t>(numPeaks_++)] = k;
                    k += 2;   // a neighbour cannot also be a peak
                }
            }
        }
 
        if (numPeaks_ == 0)
        {
            std::fill(rotRe_.begin(), rotRe_.end(), T(1));
            std::fill(rotIm_.begin(), rotIm_.end(), T(0));
            return;
        }
 
        // --- per-peak heterodyned phase propagation ----------------------------
        // analysisHop_ still holds the hop just consumed (it is recomputed
        // after this call), which is exactly the Ra of this frame interval.
        const double Ra = static_cast<double>(analysisHop_);
        const double Rs = static_cast<double>(synthHop_);
        const double binW = kTwoPi / static_cast<double>(fftSize_);
 
        int regionStart = 0;
        for (int p = 0; p < numPeaks_; ++p)
        {
            const int bin = peakBin_[static_cast<size_t>(p)];
            const double re  = static_cast<double>(spec_[static_cast<size_t>(2 * bin)]);
            const double im  = static_cast<double>(spec_[static_cast<size_t>(2 * bin + 1)]);
            const double pre = static_cast<double>(prevAnalysis_[static_cast<size_t>(2 * bin)]);
            const double pim = static_cast<double>(prevAnalysis_[static_cast<size_t>(2 * bin + 1)]);
 
            double rotR = 1.0, rotI = 0.0;
            // A peak with no spectral history is a fresh partial: keep its
            // analysis phase instead of propagating from noise.
            if (prevMag_[static_cast<size_t>(bin)] > T(0.1) * mag_[static_cast<size_t>(bin)])
            {
                // Measured inter-frame phase advance, inherently wrapped.
                const double deltaPhi = std::atan2(im * pre - re * pim, re * pre + im * pim);
                const double omegaK = binW * static_cast<double>(bin);
                const double omegaInst = omegaK + princArg(deltaPhi - omegaK * Ra) / Ra;
 
                const double psiPrev = std::atan2(
                    static_cast<double>(prevSynth_[static_cast<size_t>(2 * bin + 1)]),
                    static_cast<double>(prevSynth_[static_cast<size_t>(2 * bin)]));
                const double phi = std::atan2(im, re);
                const double theta = psiPrev + Rs * omegaInst - phi;
                rotR = std::cos(theta);
                rotI = std::sin(theta);
            }
 
            // Region of influence: up to the magnitude valley before the next peak.
            int regionEnd = numBins_;   // exclusive
            if (p + 1 < numPeaks_)
            {
                const int nextBin = peakBin_[static_cast<size_t>(p + 1)];
                int valley = bin + 1;
                T valleyMag = mag_[static_cast<size_t>(valley)];
                for (int k = bin + 2; k < nextBin; ++k)
                {
                    if (mag_[static_cast<size_t>(k)] < valleyMag)
                    {
                        valleyMag = mag_[static_cast<size_t>(k)];
                        valley = k;
                    }
                }
                regionEnd = valley + 1;
            }
 
            for (int k = regionStart; k < regionEnd; ++k)
            {
                rotRe_[static_cast<size_t>(k)] = static_cast<T>(rotR);
                rotIm_[static_cast<size_t>(k)] = static_cast<T>(rotI);
            }
            regionStart = regionEnd;
        }
 
        // DC and Nyquist are real-valued in the packed spectrum: never rotate.
        rotRe_[0] = T(1);                                     rotIm_[0] = T(0);
        rotRe_[static_cast<size_t>(numBins_ - 1)] = T(1);     rotIm_[static_cast<size_t>(numBins_ - 1)] = T(0);
    }
 
    void buildPerBinRotations() noexcept
    {
        const double Ra = static_cast<double>(analysisHop_);
        const double Rs = static_cast<double>(synthHop_);
        const double binW = kTwoPi / static_cast<double>(fftSize_);
 
        // DC and Nyquist are real-valued in the packed spectrum: never rotate.
        rotRe_[0] = T(1);                                   rotIm_[0] = T(0);
        rotRe_[static_cast<size_t>(numBins_ - 1)] = T(1);
        rotIm_[static_cast<size_t>(numBins_ - 1)] = T(0);
 
        for (int k = 1; k < numBins_ - 1; ++k)
        {
            double rotR = 1.0, rotI = 0.0;
            // A bin with no magnitude history carries no phase to propagate.
            if (prevMag_[static_cast<size_t>(k)] > T(0.1) * mag_[static_cast<size_t>(k)])
            {
                const double re  = static_cast<double>(spec_[static_cast<size_t>(2 * k)]);
                const double im  = static_cast<double>(spec_[static_cast<size_t>(2 * k + 1)]);
                const double pre = static_cast<double>(prevAnalysis_[static_cast<size_t>(2 * k)]);
                const double pim = static_cast<double>(prevAnalysis_[static_cast<size_t>(2 * k + 1)]);
 
                const double deltaPhi = std::atan2(im * pre - re * pim, re * pre + im * pim);
                const double omegaK = binW * static_cast<double>(k);
                const double omegaInst = omegaK + princArg(deltaPhi - omegaK * Ra) / Ra;
 
                const double psiPrev = std::atan2(
                    static_cast<double>(prevSynth_[static_cast<size_t>(2 * k + 1)]),
                    static_cast<double>(prevSynth_[static_cast<size_t>(2 * k)]));
                const double phi = std::atan2(im, re);
                const double theta = psiPrev + Rs * omegaInst - phi;
                rotR = std::cos(theta);
                rotI = std::sin(theta);
            }
            rotRe_[static_cast<size_t>(k)] = static_cast<T>(rotR);
            rotIm_[static_cast<size_t>(k)] = static_cast<T>(rotI);
        }
    }
 
    [[nodiscard]] static int makeOdd(int v, int lo, int hi) noexcept
    {
        v = std::clamp(v, lo, hi);
        if ((v & 1) == 0) --v;
        return std::max(v, lo | 1);
    }
 
    void pushMagnitudeHistory() noexcept
    {
        T* slot = magHist_.data() + static_cast<size_t>(histPos_) * static_cast<size_t>(numBins_);
        std::copy(mag_.begin(), mag_.end(), slot);
        histPos_ = (histPos_ + 1) % timeMedian_;
        if (histFilled_ < timeMedian_) ++histFilled_;
    }
 
    void updateSeparationMasks() noexcept
    {
        const int nHist = histFilled_;
        const int half = freqMedian_ / 2;
 
        for (int k = 0; k < numBins_; ++k)
        {
            // Median along time at this bin, over the frames seen so far.
            for (int f = 0; f < nHist; ++f)
                medianScratch_[static_cast<size_t>(f)] =
                    magHist_[static_cast<size_t>(f) * static_cast<size_t>(numBins_)
                             + static_cast<size_t>(k)];
            auto mid = medianScratch_.begin() + nHist / 2;
            std::nth_element(medianScratch_.begin(), mid,
                             medianScratch_.begin() + nHist);
            const double h = static_cast<double>(*mid);
 
            // Median along frequency in this frame, clamped at both edges so
            // the window never reads outside the spectrum.
            const int lo = std::max(0, k - half);
            const int hi = std::min(numBins_ - 1, k + half);
            const int n = hi - lo + 1;
            for (int j = 0; j < n; ++j)
                medianScratch_[static_cast<size_t>(j)] = mag_[static_cast<size_t>(lo + j)];
            auto midF = medianScratch_.begin() + n / 2;
            std::nth_element(medianScratch_.begin(), midF, medianScratch_.begin() + n);
            const double p = static_cast<double>(*midF);
 
            // A bin goes down the percussive path only when the frequency
            // median dominates the time median by the separation factor.
            // Everything else, including everything ambiguous, stays
            // harmonic: the phase-propagated path is the one that is right
            // for most material, so the split must earn each bin it takes
            // rather than divide every bin in proportion. A proportional
            // (Wiener) mask measured 1.5 dB of log-spectral distance on a
            // purely harmonic bed, where the correct answer is that nothing
            // is percussive at all.
            const bool percussive = (p > kSeparationFactor * h);
            maskHarm_[static_cast<size_t>(k)] = percussive ? T(0) : T(1);
            maskPerc_[static_cast<size_t>(k)] = percussive ? T(1) : T(0);
        }
    }
 
    void synthesizeChannel(int ch, bool isReference) noexcept
    {
        // The percussive path re-uses this channel's UNROTATED spectrum, so
        // it has to be kept before the rotation overwrites it.
        if (splitOn_)
            std::copy(spec_.begin(), spec_.end(), specRaw_.begin());
 
        // Rigid per-region rotation preserves intra-region (and inter-channel)
        // relative phases exactly.
        for (int k = 1; k < numBins_ - 1; ++k)
        {
            const T re = spec_[static_cast<size_t>(2 * k)];
            const T im = spec_[static_cast<size_t>(2 * k + 1)];
            const T rr = rotRe_[static_cast<size_t>(k)];
            const T ri = rotIm_[static_cast<size_t>(k)];
            spec_[static_cast<size_t>(2 * k)]     = re * rr - im * ri;
            spec_[static_cast<size_t>(2 * k + 1)] = re * ri + im * rr;
        }
 
        // Formant preservation: pre-warp synthesis magnitudes by
        // env(k*ratio)/env(k) so the owner's output resampler puts the
        // spectral envelope back where the input had it. The gain table is
        // computed once on the reference channel and shared, keeping the
        // stereo image coherent.
        if (formantOn_ && std::abs(ratioActive_ - 1.0) > 1e-6)
        {
            if (isReference)
                computeFormantGains();
            for (int k = 1; k < numBins_ - 1; ++k)
            {
                const T g = formantGain_[static_cast<size_t>(k)];
                spec_[static_cast<size_t>(2 * k)]     *= g;
                spec_[static_cast<size_t>(2 * k + 1)] *= g;
            }
        }
 
        // Anti-alias for upward shifts (resample-back owners): the owner's
        // resampler multiplies frequencies by `ratio`, so taper bins that
        // would land above Nyquist.
        if (resampleCompensation_ && ratioActive_ > 1.0)
        {
            const int cut = static_cast<int>(static_cast<double>(fftSize_ / 2) / ratioActive_);
            const int taperStart = std::max(1, cut - 4);
            for (int k = taperStart; k < numBins_; ++k)
            {
                T g = T(0);
                if (k <= cut)
                    g = static_cast<T>(cut - k + 1) / static_cast<T>(cut - taperStart + 1);
                spec_[static_cast<size_t>(2 * k)]     *= g;
                spec_[static_cast<size_t>(2 * k + 1)] *= g;
            }
        }
 
        if (isReference)
            std::copy(spec_.begin(), spec_.end(), prevSynth_.begin());
 
        // Harmonic/percussive recombination. The phase recursion above is fed
        // the fully propagated spectrum (the line just past), never the mix:
        // the percussive share deliberately breaks the phase recursion, and
        // letting it back in would poison the next frame's propagation.
        if (splitOn_)
        {
            for (int k = 0; k < numBins_; ++k)
            {
                const T mh = maskHarm_[static_cast<size_t>(k)];
                const T mp = maskPerc_[static_cast<size_t>(k)];
                const auto re = static_cast<size_t>(2 * k);
                const auto im = static_cast<size_t>(2 * k + 1);
                spec_[re] = spec_[re] * mh + specRaw_[re] * mp;
                spec_[im] = spec_[im] * mh + specRaw_[im] * mp;
            }
        }
 
        fft_->inverse(spec_.data(), fftResult_.data());
 
        synthesizeTail(ch);
    }
 
    void computeFormantGains() noexcept
    {
        // Even-symmetric log magnitude.
        for (int k = 0; k < numBins_; ++k)
        {
            const T re = spec_[static_cast<size_t>(2 * k)];
            const T im = spec_[static_cast<size_t>(2 * k + 1)];
            cepsTime_[static_cast<size_t>(k)] =
                std::log(std::sqrt(re * re + im * im) + T(1e-9));
        }
        for (int k = numBins_; k < fftSize_; ++k)
            cepsTime_[static_cast<size_t>(k)] =
                cepsTime_[static_cast<size_t>(fftSize_ - k)];
 
        fft_->forward(cepsTime_.data(), cepsSpec_.data());
 
        // Lifter: keep quefrencies below ~1 ms (the smooth envelope), zero
        // the rest (the harmonic comb). Cepstral index is quefrency in
        // samples: 1 ms = 0.001 * fs.
        const int qCut = std::max(8, static_cast<int>(0.001 * sampleRate_));
        const int keep = std::min(qCut, numBins_ - 1);
        for (int k = keep + 1; k < numBins_; ++k)
        {
            cepsSpec_[static_cast<size_t>(2 * k)]     = T(0);
            cepsSpec_[static_cast<size_t>(2 * k + 1)] = T(0);
        }
 
        fft_->inverse(cepsSpec_.data(), cepsTime_.data());
        for (int k = 0; k < numBins_; ++k)
            envLog_[static_cast<size_t>(k)] = cepsTime_[static_cast<size_t>(k)];
 
        // Gains with linear interpolation at k*ratio, clamped to +-40 dB.
        for (int k = 0; k < numBins_; ++k)
        {
            const double pos = std::min(static_cast<double>(k) * ratioActive_,
                                        static_cast<double>(numBins_ - 1));
            const auto i0 = static_cast<int>(pos);
            const auto frac = static_cast<T>(pos - i0);
            const T target = envLog_[static_cast<size_t>(i0)]
                + (envLog_[static_cast<size_t>(std::min(i0 + 1, numBins_ - 1))]
                   - envLog_[static_cast<size_t>(i0)]) * frac;
            const T delta = std::clamp(target - envLog_[static_cast<size_t>(k)],
                                       T(-4.6), T(4.6));
            formantGain_[static_cast<size_t>(k)] = std::exp(delta);
        }
    }
 
    void synthesizeTail(int ch) noexcept
    {
        // Overlap-add at the fixed synthesis hop. sqrt-Hann analysis+synthesis
        // with Rs = N/4 overlap-adds to a constant 2.0, hence the 0.5 norm.
        auto& acc = accum_[static_cast<size_t>(ch)];
        const auto m = static_cast<int64_t>(accumMask_);
 
        // The tail region [w + N - Rs, w + N) is entered for the first time by
        // this frame: clear the stale ring content before accumulating.
        for (int k = fftSize_ - synthHop_; k < fftSize_; ++k)
            acc[static_cast<size_t>((writeHead_ + k) & m)] = T(0);
 
        constexpr T kNorm = T(0.5);
        for (int k = 0; k < fftSize_; ++k)
        {
            const auto idx = static_cast<size_t>((writeHead_ + k) & m);
            acc[idx] += fftResult_[static_cast<size_t>(k)]
                      * window_[static_cast<size_t>(k)] * kNorm;
        }
    }
 
    // -- Members -----------------------------------------------------------------
    double sampleRate_ = 48000.0;
    int numChannels_ = 0;
    bool prepared_ = false;
    bool resampleCompensation_ = false;
 
    int fftSize_ = 2048;
    int numBins_ = 1025;
    int synthHop_ = 512;
    int ringMask_ = 2047;
    int accumSize_ = 8192;
    int accumMask_ = 8191;
 
    std::unique_ptr<FFTReal<T>> fft_;
    std::vector<T> window_;
 
    std::vector<std::vector<T>> inputRing_;   
    std::vector<std::vector<T>> accum_;       
 
    std::vector<T> fftIn_, spec_, fftResult_;
    std::vector<T> prevAnalysis_, prevSynth_; 
    std::vector<T> mag_, prevMag_;            
    std::vector<T> cepsTime_, cepsSpec_;      
    std::vector<T> envLog_, formantGain_;     
    std::vector<T> rotRe_, rotIm_;            
    std::vector<int> peakBin_;
    int numPeaks_ = 0;
 
    // Onset detection. The frame-energy test needs one scalar; the flux
    // detector needs the log-frequency filterbank tables, the two band
    // frames the flux is taken between, and the flux history the firing
    // threshold's running median is taken over. Which of the two runs is the
    // owner's choice, fixed in prepare(), and only the selected one is
    // allocated for.
    bool fluxDetectorEnabled_ = false;        
    double onsetEnv_ = 0.0;                   
    std::vector<int> fbStart_, fbOffset_, fbCount_;
    std::vector<T> fbWeights_;
    int numBands_ = 0;
    T odfScale_ = T(1);
    std::vector<T> bandCur_, bandPrev_, bandMaxPrev_;
    std::vector<double> fluxHist_;            
    std::vector<double> fluxScratch_;         
    int fluxPos_ = 0;                         
    int fluxFilled_ = 0;                      
 
    // Transient-locked hop.
    bool lockEnabled_ = false;                
    int lockFrames_ = 0;                      
    int lockLeft_ = 0;                        
    double debt_ = 0.0;                       
 
    // Harmonic/percussive split (allocated only when the owner asked for it).
    bool separationEnabled_ = false;
    int timeMedian_ = 0;                      
    int freqMedian_ = 0;                      
    std::vector<T> magHist_;                  
    std::vector<T> medianScratch_;            
    std::vector<T> maskHarm_, maskPerc_;      
    std::vector<T> specRaw_;                  
    int histPos_ = 0;                         
    int histFilled_ = 0;                      
 
    int inputPos_ = 0;
    int64_t writeHead_ = 0;
 
    // Audio-thread-private active state: the stream owner writes these as
    // plain members and never publishes them through the staged channel.
    double stActive_ = 0.0;
    double ratioActive_ = 1.0;
    double hopCarry_ = 0.0;
    int analysisHop_ = 512;
    int inputSinceHop_ = 0;
    bool firstFrame_ = true;
 
    // Audio-thread-private active copies of the published parameter set
    // (committed only by a validated bounded read; defaults == Params{}).
    double targetSemitones_ = 0.0;
    bool transientOn_ = true;
    bool formantOn_ = false;
    bool phaseLockOn_ = true;
    bool splitOn_ = false;
 
    // Control -> audio staged parameter channel (canonical seqlock; the ONLY
    // cross-thread handoff in this class).
    std::atomic<double> stgSemitones_ { 0.0 };
    std::atomic<bool> stgTransient_ { true };
    std::atomic<bool> stgFormant_ { false };
    std::atomic<bool> stgPhaseLock_ { true };
    std::atomic<bool> stgSplit_ { false };
    std::atomic<unsigned> seq_ { 0 };        
    std::atomic<bool> dirty_ { false };      
};
 
} // namespace detail
} // namespace dspark