// Measurement harness for the Preserve engine's behaviour on LOW-FREQUENCY material.
// Reports numbers; it renders no perceptual verdict and changes no DSP.
//
// The question it answers: a splice relocates the read tap by a nominal `window` refined by a
// correlation search over +/- maxLag, so the reachable relocation distances form ONE bounded
// interval. Phase-aligning a splice needs a WHOLE NUMBER OF SOURCE PERIODS inside that
// interval, and for some periods none exists — a geometric limit, separate from the
// splice-CADENCE inequality. Both now sit in time_stretch.h; this is what measured them.
//
// Two findings shaped the sections below and are worth knowing before reading the output:
// whether an unreachable multiple accumulates into a DETUNE or only wobbles the phase depends
// on whether the nearest multiple misses on one side or straddles (a straddle cancels in the
// mean); and PITCH IS THE WRONG THING TO MEASURE HERE — the fundamental usually survives, so
// the load-bearing metric is energy outside it. Zero-crossing counting in particular reports a
// wrong period on renders whose fundamental is provably correct, which is section C.
//
// It now runs every frequency-dependent section TWICE — once with splices falling back to the
// fixed window (the behaviour every number above was measured on) and once pitch-synchronous,
// with the period detected from the PCM exactly as the loader would. The two columns differ in
// that one thing, so the comparison needs no second binary and no remembered baseline.
//
// Measures, at both 44.1k and 48k geometry:
//   A. the reachable relocation interval, observed rather than derived (jump/lag/frac off
//      every SpliceEvent), and the alignment-reachability predicate over frequency.
//   B. Voice-level renders at 30 Hz: the shipped default path, then transposition at rate 1.0
//      and rate 0.5/2.0 with none, then a 20-200 Hz sweep to locate the turnover.
//   C. the P=500-frame (~88 Hz) case the pitch_shift floor probe fails on, measured with three
//      independent pitch estimators to separate the two mechanisms from a measurement artifact.
//   D. a window sweep at 30 Hz — what a larger window would buy, and what it would cost.
//   E. alignable frequencies under identical conditions, without which D and B have no scale.

#include "../src/core/instrument/engine/period_detect.h"
#include "../src/core/instrument/engine/pitch_shift.h"
#include "../src/core/instrument/engine/time_stretch.h"
#include "../src/core/instrument/engine/voice.h"
#include "tone_metrics.h"

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <utility>
#include <vector>

using namespace reasampler;
using namespace reasampler::instrument::engine;

static int g_fail = 0;
#define CHECK(cond) do { if(!(cond)) { \
    std::printf("FAIL line %d: %s\n", __LINE__, #cond); ++g_fail; } } while(0)

constexpr double kPi = 3.14159265358979323846;

// Whether a source carries its detected period into the shifter — i.e. whether splices are
// pitch-synchronous or fall back to the fixed window. Every section below runs under whichever
// is set, so main() can drive the SAME measurements both ways from one binary and the two
// columns are comparable by construction.
static bool g_pitchSynchronous = true;

// ---------------------------------------------------------------------------------------
// Source + render helpers
// ---------------------------------------------------------------------------------------

// A pure sine at `freqHz`, phase-continuous, long enough that a rate-2.0 render never
// exhausts it (the caller sizes `frames`). The period is DETECTED rather than computed from
// freqHz on purpose: that is the number the loader would actually hand the engine, so the
// measurement includes any detector error rather than assuming it away.
static SampleData sineSample(double freqHz, int sampleRate, std::size_t frames,
                             PitchEngine engine, double phase = 0.0) {
    SampleData s;
    s.frames.resize(frames);
    const double w = 2.0 * kPi * freqHz / static_cast<double>(sampleRate);
    for (std::size_t i = 0; i < frames; ++i) {
        s.frames[i] = static_cast<AudioSample>(std::sin(w * static_cast<double>(i) + phase));
    }
    s.sampleRate = sampleRate;
    s.rootNote = 60;
    s.play.pitchEngine = engine;  // Gate, no loop, default (fully open) AHDSR
    if (g_pitchSynchronous) s.sourcePeriodFrames = detectPeriod(s.frames, sampleRate).frames;
    return s;
}

// One note through the REAL Voice: presize (off-thread step), start with the stretch rate,
// then pull `outFrames` mono frames. `note - 60` is the transposition in semitones.
static std::vector<double> renderVoice(const SampleData& s, int note, double stretchRate,
                                       std::int64_t window, std::size_t outFrames) {
    Voice v;
    v.presizePreserveShifters(window);
    v.start(note, 127, s, /*declickTakeover=*/false, stretchRate);
    std::vector<double> out(outFrames, 0.0);
    for (std::size_t i = 0; i < outFrames; ++i) {
        out[i] = static_cast<double>(v.renderFrame());
    }
    return out;
}

// ---------------------------------------------------------------------------------------
// Metrics
// ---------------------------------------------------------------------------------------

// Mean spacing between positive-going zero crossings over [from, to) — the same estimator
// test_pitch_shift.cpp uses, kept identical so the two files' numbers are comparable.
static double periodIn(const std::vector<double>& v, std::size_t from, std::size_t to) {
    double sum = 0.0;
    std::size_t prev = 0, count = 0;
    for (std::size_t i = from + 1; i < to && i < v.size(); ++i) {
        if (v[i - 1] <= 0.0 && v[i] > 0.0) {
            if (count > 0) sum += static_cast<double>(i - prev);
            prev = i;
            ++count;
        }
    }
    return count > 1 ? sum / static_cast<double>(count - 1) : 0.0;
}

// Least-squares fit of a single tone at `cyclesPerFrame` over [from, from+len): returns the
// fitted phase and writes the residual energy fraction (1 - explained), which is the
// single-tone-purity metric — 0 = a perfect sine at that frequency, 1 = none of the energy
// is there. Robust to amplitude but NOT to phase drift within the block, which is why the
// caller keeps blocks near one period.
static double toneFit(const std::vector<double>& v, std::size_t from, std::size_t len,
                      double cyclesPerFrame, double* residFraction) {
    double sc = 0.0, ss = 0.0, cc = 0.0, s2 = 0.0, cs = 0.0, e = 0.0;
    for (std::size_t k = 0; k < len && from + k < v.size(); ++k) {
        const double t = 2.0 * kPi * cyclesPerFrame * static_cast<double>(k);
        const double c = std::cos(t), s = std::sin(t);
        const double x = v[from + k];
        sc += x * c; ss += x * s; cc += c * c; s2 += s * s; cs += c * s; e += x * x;
    }
    const double det = cc * s2 - cs * cs;
    double a = 0.0, b = 0.0;
    if (std::fabs(det) > 1e-12) {
        a = (sc * s2 - ss * cs) / det;
        b = (ss * cc - sc * cs) / det;
    }
    const double explained = a * sc + b * ss;  // energy captured by the fit
    if (residFraction != nullptr) *residFraction = e > 0.0 ? 1.0 - explained / e : 0.0;
    return std::atan2(b, a);
}

// Total unwrapped phase drift (in CYCLES) of the render relative to an ideal tone at
// `cyclesPerFrame`, measured across [from, to) in one-period blocks. This is the direct
// observable behind "the rendered pitch is wrong": a nonzero drift IS a frequency error.
static double phaseDriftCycles(const std::vector<double>& v, std::size_t from, std::size_t to,
                               double cyclesPerFrame, double* worstStepCycles) {
    const std::size_t blk = static_cast<std::size_t>(1.0 / cyclesPerFrame);
    double total = 0.0, prev = 0.0, worst = 0.0;
    bool first = true;
    for (std::size_t p = from; p + blk <= to && p + blk < v.size(); p += blk) {
        const double ph = toneFit(v, p, blk, cyclesPerFrame, nullptr);
        if (!first) {
            double d = ph - prev;
            while (d > kPi) d -= 2.0 * kPi;
            while (d < -kPi) d += 2.0 * kPi;
            total += d / (2.0 * kPi);
            if (std::fabs(d) / (2.0 * kPi) > worst) worst = std::fabs(d) / (2.0 * kPi);
        }
        prev = ph;
        first = false;
    }
    if (worstStepCycles != nullptr) *worstStepCycles = worst;
    return total;
}

// Median single-tone residual fraction over the render, in one-period blocks.
static double medianResidual(const std::vector<double>& v, std::size_t from, std::size_t to,
                             double cyclesPerFrame) {
    const std::size_t blk = static_cast<std::size_t>(1.0 / cyclesPerFrame);
    std::vector<double> r;
    for (std::size_t p = from; p + blk <= to && p + blk < v.size(); p += blk) {
        double resid = 0.0;
        toneFit(v, p, blk, cyclesPerFrame, &resid);
        r.push_back(resid);
    }
    if (r.empty()) return 0.0;
    std::size_t mid = r.size() / 2;
    std::nth_element(r.begin(), r.begin() + static_cast<std::ptrdiff_t>(mid), r.end());
    return r[mid];
}

// A pitch estimator that, unlike zero-crossing counting, is not fooled by a low-level fast
// component adding spurious crossings. The two disagreeing is itself the diagnosis. Shared with
// the gated tests (tone_metrics.h) so there is one estimator and not two.
using reasampler::test_support::autocorrelationPeriod;

// The five strongest spectral peaks over a Hann-windowed segment, scanned on a fine period
// grid (Goertzel-style direct evaluation, no FFT-bin quantization). Prints period in frames
// and magnitude relative to the strongest — the decisive "is the rendered pitch wrong, or is
// there a second component fooling the zero-crossing count" measurement.
static void reportSpectrum(const char* label, const std::vector<double>& v, std::size_t from,
                           std::size_t len, double wantPeriod) {
    const int kGrid = 2000;
    const double pLo = 30.0, pHi = 8000.0;
    std::vector<double> mag(static_cast<std::size_t>(kGrid), 0.0);
    std::vector<double> per(static_cast<std::size_t>(kGrid), 0.0);
    for (int g = 0; g < kGrid; ++g) {
        // Geometric grid: constant relative resolution across three octaves of period.
        const double p = pLo * std::pow(pHi / pLo, static_cast<double>(g) / (kGrid - 1));
        per[static_cast<std::size_t>(g)] = p;
        double re = 0.0, im = 0.0;
        const double w = 2.0 * kPi / p;
        for (std::size_t k = 0; k < len && from + k < v.size(); ++k) {
            const double hann = 0.5 * (1.0 - std::cos(2.0 * kPi * static_cast<double>(k) /
                                                      static_cast<double>(len)));
            const double x = v[from + k] * hann;
            re += x * std::cos(w * static_cast<double>(k));
            im += x * std::sin(w * static_cast<double>(k));
        }
        mag[static_cast<std::size_t>(g)] = std::sqrt(re * re + im * im);
    }
    double top = 0.0;
    for (double m : mag) top = std::max(top, m);
    // Local maxima, ranked by MAGNITUDE (not by grid order) so the fundamental cannot be
    // pushed off the list by low-level debris at a shorter period.
    std::vector<std::pair<double, double>> peaks;  // (magnitude, period)
    for (int g = 1; g + 1 < kGrid; ++g) {
        const std::size_t i = static_cast<std::size_t>(g);
        if (mag[i] <= mag[i - 1] || mag[i] < mag[i + 1]) continue;
        if (mag[i] < 0.02 * top) continue;
        peaks.emplace_back(mag[i], per[i]);
    }
    std::sort(peaks.begin(), peaks.end(),
              [](const std::pair<double, double>& a, const std::pair<double, double>& b) {
                  return a.first > b.first;
              });
    // Energy fraction OUTSIDE the fundamental's mainlobe — the honest "how much of this render
    // is not the wanted tone" number, since a peak list alone can hide broadband debris.
    double eTotal = 0.0, eFund = 0.0;
    for (int g = 0; g < kGrid; ++g) {
        const std::size_t i = static_cast<std::size_t>(g);
        const double e = mag[i] * mag[i];
        eTotal += e;
        if (std::fabs(per[i] - wantPeriod) / wantPeriod < 0.06) eFund += e;
    }
    std::printf("    %s spectrum (want period %.1f fr); strongest peaks >2%% of max:\n", label,
                wantPeriod);
    for (std::size_t k = 0; k < peaks.size() && k < 8; ++k) {
        std::printf("      period %8.1f fr  rel %.4f%s\n", peaks[k].second,
                    peaks[k].first / top,
                    std::fabs(peaks[k].second - wantPeriod) / wantPeriod < 0.03
                        ? "   <-- the wanted tone" : "");
    }
    std::printf("      energy outside the wanted tone's mainlobe: %.2f%%\n",
                eTotal > 0.0 ? 100.0 * (1.0 - eFund / eTotal) : 0.0);
}

// ---------------------------------------------------------------------------------------
// A. Splice geometry, observed off the shifter's own SpliceEvent stream
// ---------------------------------------------------------------------------------------

struct SpliceStats {
    long long count = 0;
    double minReloc = 1e18, maxReloc = -1e18;
    std::int64_t minLag = 1LL << 40, maxLag = -(1LL << 40);
    double meanInterval = 0.0;
    bool jumpAlwaysNominal = true;  // |jump| == the nominal on every splice (steady state)
    std::int64_t nominalJump = 0;   // what the shifter itself resolved the nominal to
};

// Drives a bare PitchShifter over the same feed schedule Voice uses, recording every splice.
// The audio is not kept — this measures the DECISIONS, not the sound.
static SpliceStats spliceGeometry(const std::vector<AudioSample>& src, std::int64_t window,
                                  double rate, double shift, std::size_t outFrames,
                                  std::vector<double>* audio = nullptr,
                                  double sourcePeriod = 0.0) {
    PitchShifter ps;
    ps.configure(window);
    ps.prime(src.data(), window);
    ps.setShiftRatio(shift);
    ps.setFeedRate(rate);
    ps.setSourcePeriod(sourcePeriod);
    StretchCursor cur;
    cur.start(window);
    loop::ResolvedLoop lp{};  // inactive: the source is long enough to run straight through

    SpliceStats st;
    st.nominalJump = ps.spliceJump();
    if (audio != nullptr) audio->assign(outFrames, 0.0);
    std::size_t lastSpliceAt = 0;
    double intervalSum = 0.0;
    long long intervals = 0;
    for (std::size_t i = 0; i < outFrames; ++i) {
        const std::int64_t due = cur.due(rate);
        AudioSample last = 0.0f;
        bool fed = false;
        for (std::int64_t k = 0; k < due; ++k) {
            if (fed) ps.writeFrame(last);
            const std::int64_t q = cur.next(lp);
            last = (q >= 0 && static_cast<std::size_t>(q) < src.size())
                       ? src[static_cast<std::size_t>(q)] : 0.0f;
            fed = true;
        }
        const AudioSample o = fed ? ps.process(last) : ps.processNoInput();
        if (audio != nullptr) (*audio)[i] = static_cast<double>(o);
        const SpliceEvent& ev = ps.lastSplice();
        if (!ev.fired) continue;
        ++st.count;
        const double reloc = std::fabs(static_cast<double>(ev.jump) -
                                       static_cast<double>(ev.lag) - ev.frac);
        if (reloc < st.minReloc) st.minReloc = reloc;
        if (reloc > st.maxReloc) st.maxReloc = reloc;
        if (ev.lag < st.minLag) st.minLag = ev.lag;
        if (ev.lag > st.maxLag) st.maxLag = ev.lag;
        if (std::llabs(ev.jump) != st.nominalJump) st.jumpAlwaysNominal = false;
        if (lastSpliceAt != 0) { intervalSum += static_cast<double>(i - lastSpliceAt); ++intervals; }
        lastSpliceAt = i;
    }
    st.meanInterval = intervals > 0 ? intervalSum / static_cast<double>(intervals) : 0.0;
    return st;
}

// Is there a whole number of source periods inside the reachable relocation interval?
static bool alignmentReachable(double periodFrames, double lo, double hi, int* whichN) {
    for (int n = 1; n <= 64; ++n) {
        const double m = periodFrames * n;
        if (m > hi) break;
        if (m >= lo) { if (whichN != nullptr) *whichN = n; return true; }
    }
    if (whichN != nullptr) *whichN = 0;
    return false;
}

// ---------------------------------------------------------------------------------------
// 1. The reachable relocation interval, measured
// ---------------------------------------------------------------------------------------

static void reportReachableInterval() {
    std::printf("\n=== A. Reachable splice relocation interval (measured) ===\n");
    for (const auto& g : {std::pair<int, std::int64_t>{44100, 2205},
                          std::pair<int, std::int64_t>{48000, 2400}}) {
        const int sr = g.first;
        const std::int64_t w = g.second;
        // Broadband noise: every lag is a plausible candidate, so the search's own limits —
        // not the source's periodicity — set the observed extremes.
        std::vector<AudioSample> noise(600000);
        std::uint32_t rng = 12345u;
        for (auto& x : noise) {
            rng = rng * 1664525u + 1013904223u;
            x = static_cast<AudioSample>((static_cast<double>(rng >> 8) / 8388608.0) - 1.0);
        }
        SpliceStats up = spliceGeometry(noise, w, 1.0, 1.5, 120000);   // up-shift: jump +w
        SpliceStats dn = spliceGeometry(noise, w, 1.0, 0.7, 120000);   // down-shift: jump -w
        const double lo = std::min(up.minReloc, dn.minReloc);
        const double hi = std::max(up.maxReloc, dn.maxReloc);
        std::printf("  %d Hz, window %lld frames (%.1f ms):\n", sr, static_cast<long long>(w),
                    1000.0 * static_cast<double>(w) / sr);
        std::printf("    up-shift   splices %lld, lag [%lld, %lld], reloc [%.2f, %.2f]\n",
                    up.count, static_cast<long long>(up.minLag),
                    static_cast<long long>(up.maxLag), up.minReloc, up.maxReloc);
        std::printf("    down-shift splices %lld, lag [%lld, %lld], reloc [%.2f, %.2f]\n",
                    dn.count, static_cast<long long>(dn.minLag),
                    static_cast<long long>(dn.maxLag), dn.minReloc, dn.maxReloc);
        std::printf("    observed reachable relocation interval: [%.2f, %.2f] frames "
                    "= [%.2f, %.2f] ms\n", lo, hi, 1000.0 * lo / sr, 1000.0 * hi / sr);
        std::printf("    structural bound (window +/- window/4): [%lld, %lld]\n",
                    static_cast<long long>(w - w / 4), static_cast<long long>(w + w / 4));
        // The jump is nominal and the lag is inside +/- window/4 — the two facts the
        // reachable interval is derived from.
        CHECK(up.jumpAlwaysNominal && dn.jumpAlwaysNominal);
        CHECK(up.minLag >= -(w / 4) && up.maxLag <= w / 4);
        CHECK(dn.minLag >= -(w / 4) && dn.maxLag <= w / 4);
    }
}

// The reachability predicate over frequency, at both geometries. Pure arithmetic over the
// interval measured above — no render, stated as such.
static void reportReachabilityByFrequency() {
    std::printf("\n=== A2. Alignment reachability by frequency (arithmetic, not rendered) ===\n");
    const double freqs[] = {12, 14, 16, 18, 20, 22, 24, 26, 26.6, 28, 30, 31, 31.9, 32,
                            34, 36, 40, 50, 60, 80, 88.2, 100, 140, 200};
    for (const auto& g : {std::pair<int, std::int64_t>{44100, 2205},
                          std::pair<int, std::int64_t>{48000, 2400}}) {
        const int sr = g.first;
        const std::int64_t w = g.second;
        const double lo = static_cast<double>(w - w / 4), hi = static_cast<double>(w + w / 4);
        std::printf("  %d Hz / window %lld, interval [%.0f, %.0f] frames:\n", sr,
                    static_cast<long long>(w), lo, hi);
        for (double f : freqs) {
            const double P = static_cast<double>(sr) / f;
            int n = 0;
            const bool ok = alignmentReachable(P, lo, hi, &n);
            if (ok) {
                std::printf("    %6.1f Hz  P=%8.1f  ALIGNABLE (n=%d, n*P=%.1f)\n", f, P, n,
                            n * P);
            } else {
                // How far the nearest multiple sits outside the interval, and the phase error
                // that residual forces at every splice.
                double best = 1e18; double bestM = 0.0;
                for (int k = 1; k <= 64; ++k) {
                    const double m = P * k;
                    const double d = m < lo ? lo - m : (m > hi ? m - hi : 0.0);
                    if (d < best) { best = d; bestM = m; }
                }
                std::printf("    %6.1f Hz  P=%8.1f  UNALIGNABLE (nearest n*P=%.1f, off by "
                            "%.1f frames = %.1f deg of phase)\n",
                            f, P, bestM, best, 360.0 * best / P);
            }
        }
    }
}

// ---------------------------------------------------------------------------------------
// 2. Voice-level renders at 30 Hz
// ---------------------------------------------------------------------------------------

// The reassurance case: the SHIPPED default path. 30 Hz played at its root, rate 1.0, no
// transposition. The shift is exactly 1.0, so the tap's delay never drifts and no splice can
// fire; a primed shifter at unity is a bit-exact pass-through. Baseline is the SAME source
// through Varispeed at the root, which is a straight readPos_ += 1.0 read of the PCM — i.e.
// the unprocessed sample. This is NOT a comparison against a pre-change binary; it is the
// stronger claim that the path is transparent.
static void testRootRateUnityIsBitIdenticalToTheDirectRead() {
    std::printf("\n=== B1. 30 Hz, root note, rate 1.0, no transposition ===\n");
    const int sr = 44100;
    const std::int64_t w = 2205;
    const std::size_t frames = 300000, outFrames = 250000;
    const SampleData pres = sineSample(30.0, sr, frames, PitchEngine::Preserve);
    const SampleData vari = sineSample(30.0, sr, frames, PitchEngine::Varispeed);
    const std::vector<double> p = renderVoice(pres, 60, 1.0, w, outFrames);
    const std::vector<double> v = renderVoice(vari, 60, 1.0, w, outFrames);
    std::size_t firstDiff = outFrames;
    for (std::size_t i = 0; i < outFrames; ++i) {
        if (p[i] != v[i]) { firstDiff = i; break; }
    }
    std::printf("  Preserve vs Varispeed at root, %zu frames: %s\n", outFrames,
                firstDiff == outFrames ? "BIT-IDENTICAL"
                                       : "differ (first at frame ?)");
    if (firstDiff != outFrames) {
        std::printf("  first difference at frame %zu (%.9f vs %.9f)\n", firstDiff, p[firstDiff],
                    v[firstDiff]);
    }
    CHECK(firstDiff == outFrames);
    // And the same claim at the 48k geometry.
    const SampleData pres48 = sineSample(30.0, 48000, frames, PitchEngine::Preserve);
    const SampleData vari48 = sineSample(30.0, 48000, frames, PitchEngine::Varispeed);
    const std::vector<double> p48 = renderVoice(pres48, 60, 1.0, 2400, outFrames);
    const std::vector<double> v48 = renderVoice(vari48, 60, 1.0, 2400, outFrames);
    bool same48 = true;
    for (std::size_t i = 0; i < outFrames && same48; ++i) if (p48[i] != v48[i]) same48 = false;
    std::printf("  same at 48k / window 2400: %s\n", same48 ? "BIT-IDENTICAL" : "DIFFER");
    CHECK(same48);
    // Splice count on the same conditions, read off the shifter directly.
    std::vector<AudioSample> src(pres.frames.begin(), pres.frames.end());
    const SpliceStats st = spliceGeometry(src, w, 1.0, 1.0, outFrames);
    std::printf("  splices fired over %zu frames at shift 1.0, rate 1.0: %lld\n", outFrames,
                st.count);
    CHECK(st.count == 0);

    // Onset at a MUCH larger window — the cost side of any window-resize option. prime()
    // parks the tap on src[0] whatever the window, so frame 0 must still be source frame 0.
    // Started at quarter-phase so src[0] is FULL SCALE, not the zero a sine would give: a
    // frame-0 match against 0.0 would also pass on a voice that produced silence.
    const SampleData cosPhase =
        sineSample(30.0, sr, 300000, PitchEngine::Preserve, kPi / 2.0);
    CHECK(cosPhase.frames[0] == 1.0f);
    for (std::int64_t big : {std::int64_t{2205}, std::int64_t{8820}}) {
        const std::vector<double> up = renderVoice(cosPhase, 67, 1.0, big, 64);  // +7 st
        std::printf("  window %5lld, +7 st: out[0]=%.9f (src[0]=%.9f), |out| over frames 1..63 "
                    "min %.6f\n", static_cast<long long>(big), up[0],
                    static_cast<double>(cosPhase.frames[0]),
                    *std::min_element(up.begin() + 1, up.end(),
                        [](double a, double b) { return std::fabs(a) < std::fabs(b); }));
        CHECK(up[0] == static_cast<double>(cosPhase.frames[0]));  // zero added latency
    }
    // A sample SHORTER than the window: start() primes the whole playable span and freezes
    // the writer immediately. A larger window moves that threshold, so check it still speaks
    // on frame 0 at the largest window swept below.
    const SampleData shortSample =
        sineSample(30.0, sr, 3000, PitchEngine::Preserve, kPi / 2.0);
    const std::vector<double> shortOut = renderVoice(shortSample, 67, 1.0, 8820, 64);
    std::printf("  3000-frame sample under a 8820-frame window, +7 st: out[0]=%.9f src[0]=%.9f\n",
                shortOut[0], static_cast<double>(shortSample.frames[0]));
    CHECK(shortOut[0] == static_cast<double>(shortSample.frames[0]));
}

// One measured row: render through the Voice and report every metric for that condition.
static void measureRow(const char* label, double freqHz, int sr, std::int64_t window,
                       int note, double rate) {
    const std::size_t frames = 900000;
    const std::size_t outFrames = 300000;
    const SampleData s = sineSample(freqHz, sr, frames, PitchEngine::Preserve);
    const double shift = std::pow(2.0, (note - 60) / 12.0);
    const std::vector<double> out = renderVoice(s, note, rate, window, outFrames);

    bool finite = true;
    double peak = 0.0;
    for (double x : out) { if (!std::isfinite(x)) finite = false; peak = std::max(peak, std::fabs(x)); }

    const double srcPeriod = static_cast<double>(sr) / freqHz;
    const double wantPeriod = srcPeriod / shift;           // pitch is the TAP's, not the feed's
    const double wantCpf = 1.0 / wantPeriod;
    const std::size_t from = 40000, to = 280000;
    const double gotPeriod = periodIn(out, from, to);
    double worstStep = 0.0;
    const double drift = phaseDriftCycles(out, from, to, wantCpf, &worstStep);
    const double resid = medianResidual(out, from, to, wantCpf);

    std::vector<AudioSample> src(s.frames.begin(), s.frames.end());
    const SpliceStats st =
        spliceGeometry(src, window, rate, shift, outFrames, nullptr, s.sourcePeriodFrames);
    const double lo = static_cast<double>(window - window / 4);
    const double hi = static_cast<double>(window + window / 4);
    int n = 0;
    // Reachability of the FIXED-window interval. With a source period known this is no longer
    // the binding question — the nominal jump is a multiple of the period by construction —
    // but it stays reported because it is what the "NO" rows below were diagnosed by.
    const bool reach = alignmentReachable(srcPeriod, lo, hi, &n);

    // Effective frequency error implied by the drift, and the phase step it works out to per
    // splice — the number that says whether a splice is stepping the phase or not.
    const double driftPerSplice = st.count > 0 ? drift / static_cast<double>(st.count) : 0.0;
    std::printf("  %-26s f=%6.1f Hz shift=%.4f rate=%.2f | P det %8.2f jump %5lld | "
                "period got %8.2f want %8.2f (%+.2f%%) | splices %4lld every %7.0f fr | "
                "phase drift %+8.3f cyc (%+7.1f deg/splice, worst step %.1f deg) | "
                "resid %.4f | peak %.3f | fixed-window align %s%s\n",
                label, freqHz, shift, rate, s.sourcePeriodFrames,
                static_cast<long long>(st.nominalJump), gotPeriod, wantPeriod,
                wantPeriod > 0.0 ? 100.0 * (gotPeriod - wantPeriod) / wantPeriod : 0.0,
                st.count, st.meanInterval, drift, 360.0 * driftPerSplice, 360.0 * worstStep,
                resid, peak, reach ? "YES" : "NO",
                reach ? "" : " <-- no whole period in the fixed-window reachable interval");
    CHECK(finite);
}

// Three independent pitch estimators plus the spectrum, on one condition. Where the
// zero-crossing count and the autocorrelation disagree, the render is not simply detuned —
// something else is crossing zero.
static void deepDive(const char* label, double freqHz, int sr, std::int64_t window, int note,
                     double rate) {
    const SampleData s = sineSample(freqHz, sr, 900000, PitchEngine::Preserve);
    const double shift = std::pow(2.0, (note - 60) / 12.0);
    const std::vector<double> out = renderVoice(s, note, rate, window, 300000);
    const double wantPeriod = (static_cast<double>(sr) / freqHz) / shift;
    const double zc = periodIn(out, 40000, 280000);
    // Search bounded to [0.5, 1.7] x the wanted period: a pure sine autocorrelates equally at
    // EVERY multiple of its period, so an unbounded search reports 2P about half the time.
    const double ac = autocorrelationPeriod(out, 60000, 60000,
                                     std::max<std::int64_t>(40,
                                         static_cast<std::int64_t>(wantPeriod * 0.5)),
                                     static_cast<std::int64_t>(wantPeriod * 1.7));
    std::printf("  %s (f=%.1f Hz, shift %.4f, rate %.2f, want period %.1f fr):\n", label, freqHz,
                shift, rate, wantPeriod);
    std::printf("    zero-crossing period %.2f (%+.2f%%) | autocorrelation period %.2f "
                "(%+.2f%%)\n", zc, 100.0 * (zc - wantPeriod) / wantPeriod, ac,
                100.0 * (ac - wantPeriod) / wantPeriod);
    reportSpectrum(label, out, 60000, 131072, wantPeriod);
}

static void reportTransposedAt30Hz() {
    std::printf("\n=== B2. 30 Hz transposed, rate 1.0 (44.1k / window 2205) ===\n");
    measureRow("30 Hz +2 st", 30.0, 44100, 2205, 62, 1.0);
    measureRow("30 Hz +7 st", 30.0, 44100, 2205, 67, 1.0);
    measureRow("30 Hz -7 st", 30.0, 44100, 2205, 53, 1.0);
    std::printf("\n=== B3. 30 Hz stretched, no transposition (44.1k / window 2205) ===\n");
    measureRow("30 Hz rate 0.5", 30.0, 44100, 2205, 60, 0.5);
    measureRow("30 Hz rate 2.0", 30.0, 44100, 2205, 60, 2.0);
    std::printf("\n=== B4. the same six at 48k / window 2400 ===\n");
    measureRow("30 Hz +2 st @48k", 30.0, 48000, 2400, 62, 1.0);
    measureRow("30 Hz +7 st @48k", 30.0, 48000, 2400, 67, 1.0);
    measureRow("30 Hz rate 2.0 @48k", 30.0, 48000, 2400, 60, 2.0);
}

// Where does the behaviour actually turn over? Swept at a fixed, modest transposition so the
// only thing changing is the source period against the reachable interval.
static void reportFrequencySweep() {
    std::printf("\n=== B5. Frequency sweep, +2 st, rate 1.0 (44.1k / window 2205) ===\n");
    const double freqs[] = {20, 22, 24, 25, 26, 26.5, 27, 28, 29, 30, 31, 31.5, 32, 33,
                            34, 36, 40, 45, 50, 60, 70, 80, 88.2, 100, 120, 140, 170, 200};
    for (double f : freqs) measureRow("sweep +2 st", f, 44100, 2205, 62, 1.0);
    std::printf("\n=== B6. Same sweep at rate 2.0, NO transposition ===\n");
    for (double f : freqs) measureRow("sweep rate 2.0", f, 44100, 2205, 60, 2.0);
}

// ---------------------------------------------------------------------------------------
// 3. The P=500 case: geometric, or cadence?
// ---------------------------------------------------------------------------------------

// pitch_shift_tests' floor probe fails at source period 500 frames, rate 2.0, shift 0.25 and
// holds at 600/700. Under the geometric claim, alignment needs a whole number of source
// periods in [0.75w, 1.25w] = [1654, 2756]. This reports whether that is satisfied for each of
// the three periods — separating "no aligned landing point exists" (geometric) from "an
// aligned landing point exists but the cadence is too fast to use it" (the inequality already
// in time_stretch.h).
static void reportFloorProbeMechanism() {
    std::printf("\n=== C. The P=500/600/700 floor probe: which mechanism? ===\n");
    const std::int64_t w = 2205;
    const int sr = 44100;
    const double rate = 2.0;
    const double shift = std::pow(2.0, -24.0 / 12.0);  // 0.25
    const double lo = static_cast<double>(w - w / 4), hi = static_cast<double>(w + w / 4);
    for (double period : {500.0, 600.0, 700.0}) {
        int n = 0;
        const bool reach = alignmentReachable(period, lo, hi, &n);
        const double freq = static_cast<double>(sr) / period;
        // The cadence inequality from time_stretch.h, evaluated for this row against the
        // ACTUAL nominal jump this geometry resolves to — under g_pitchSynchronous that is
        // periodAlignedJump's answer, not always the fixed window, so the two passes of this
        // report (fixed-window / pitch-synchronous) must not print the same number.
        PitchShifter jumpProbe;
        jumpProbe.configure(w);
        jumpProbe.setSourcePeriod(g_pitchSynchronous ? period : 0.0);
        const double cadence =
            static_cast<double>(jumpProbe.spliceJump()) / std::fabs(rate - shift);
        const double outPeriod = period / shift;
        std::printf("  P=%5.0f (%.1f Hz): alignable in [%.0f,%.0f]? %s%s | cadence %.0f fr vs "
                    "output period %.0f fr -> %s\n",
                    period, freq, lo, hi, reach ? "YES" : "NO",
                    reach ? "" : " (geometric failure)",
                    cadence, outPeriod,
                    cadence < outPeriod ? "SPLICE INSIDE A CYCLE (cadence failure)" : "ok");
        measureRow("floor probe", freq, sr, w, 60 - 24, rate);
    }
    // The probe's OWN signal, reproduced exactly: a bare PitchShifter fed by the same
    // schedule test_pitch_shift.cpp's runStretch uses, not the Voice. Its zero-crossing
    // number is the one that is currently RED, so it is the one that has to be explained.
    std::printf("\n  --- the probe's exact signal (bare PitchShifter, runStretch schedule) ---\n");
    for (double period : {500.0, 600.0, 700.0}) {
        const std::size_t srcLen = 400000;
        std::vector<AudioSample> src(srcLen);
        for (std::size_t i = 0; i < srcLen; ++i) {
            src[i] = static_cast<AudioSample>(
                std::sin(2.0 * kPi * static_cast<double>(i) / period));
        }
        std::vector<double> out;
        const SpliceStats st = spliceGeometry(src, w, rate, shift, 60000, &out,
                                              g_pitchSynchronous ? period : 0.0);
        const double want = period / shift;
        const double zc = periodIn(out, 20000, 50000);
        const double ac = autocorrelationPeriod(out, 20000, 20000,
                                         static_cast<std::int64_t>(want * 0.5),
                                         static_cast<std::int64_t>(want * 1.7));
        std::printf("  P=%.0f: zero-crossing %.2f (%+.2f%%) | autocorrelation %.2f (%+.2f%%) "
                    "| splices %lld every %.0f fr\n", period, zc, 100.0 * (zc - want) / want,
                    ac, 100.0 * (ac - want) / want, st.count, st.meanInterval);
        reportSpectrum("probe", out, 20000, 32768, want);
    }

    std::printf("\n  --- independent pitch estimators on the same three (through the Voice) ---\n");
    deepDive("P=500", 44100.0 / 500.0, sr, w, 36, rate);
    deepDive("P=600", 44100.0 / 600.0, sr, w, 36, rate);
    deepDive("P=700", 44100.0 / 700.0, sr, w, 36, rate);
    std::printf("\n  --- and on the geometric cases, for contrast ---\n");
    deepDive("30 Hz +2 st rate 1.0", 30.0, sr, w, 62, 1.0);
    deepDive("30 Hz rate 2.0", 30.0, sr, w, 60, 2.0);
    deepDive("29 Hz rate 2.0", 29.0, sr, w, 60, 2.0);
}

// What would a bigger window buy? The reachable interval is [0.75w, 1.25w], so it contains a
// whole number of source periods for EVERY period P <= 0.5w — i.e. a window of at least TWO
// source periods makes alignment reachable unconditionally. This sweeps 30 Hz across windows
// spanning that threshold (1470 * 2 = 2940 frames = 66.7 ms at 44.1k) and reports what
// actually changes. The window is an ARGUMENT to configure(); nothing shipped is altered.
static void reportWindowSweep() {
    std::printf("\n=== D. Window sweep at 30 Hz — what a larger window would buy ===\n");
    const int sr = 44100;
    const double P = static_cast<double>(sr) / 30.0;
    std::printf("  source period %.1f frames; alignment is unconditional once window >= 2P = "
                "%.0f frames (%.1f ms)\n", P, 2.0 * P, 2000.0 * P / sr);
    for (std::int64_t w : {std::int64_t{2205}, std::int64_t{2646}, std::int64_t{2940},
                           std::int64_t{3528}, std::int64_t{4410}, std::int64_t{8820}}) {
        const double lo = static_cast<double>(w - w / 4), hi = static_cast<double>(w + w / 4);
        int n = 0;
        const bool reach = alignmentReachable(P, lo, hi, &n);
        std::printf("\n  window %lld fr (%.1f ms), interval [%.0f, %.0f]: %s\n",
                    static_cast<long long>(w), 1000.0 * static_cast<double>(w) / sr, lo, hi,
                    reach ? "ALIGNABLE" : "unalignable");
        // Per-voice Preserve state: two shifter rings of 2*window floats (L/R) plus the
        // window-sized prime scratch = 5*window floats (voice.cpp presizePreserveShifters).
        const double bytes = 5.0 * static_cast<double>(w) * 4.0;
        std::printf("    per-voice Preserve state %.1f KB; at the 32-voice ceiling %.2f MB\n",
                    bytes / 1024.0, 32.0 * bytes / (1024.0 * 1024.0));
        measureRow("  30 Hz +2 st", 30.0, sr, w, 62, 1.0);
        measureRow("  30 Hz rate 2.0", 30.0, sr, w, 60, 2.0);
        deepDive("    +2 st", 30.0, sr, w, 62, 1.0);
        deepDive("    rate 2.0", 30.0, sr, w, 60, 2.0);
    }
}

// Alignable neighbours under identical conditions — without these the out-of-band-energy
// numbers above have no scale.
static void reportAlignableControls() {
    std::printf("\n=== E. Alignable controls (same conditions, a frequency that CAN align) ===\n");
    deepDive("34 Hz +2 st rate 1.0", 34.0, 44100, 2205, 62, 1.0);
    deepDive("34 Hz rate 2.0", 34.0, 44100, 2205, 60, 2.0);
    deepDive("20 Hz +2 st rate 1.0", 20.0, 44100, 2205, 62, 1.0);
    deepDive("220 Hz +2 st rate 1.0", 220.0, 44100, 2205, 62, 1.0);
    deepDive("220 Hz rate 2.0", 220.0, 44100, 2205, 60, 2.0);
}

// The frequency-dependent sections, run under whichever splice geometry is set. Everything
// that can differ between the two is in here; section A (the reachable interval, measured on
// noise) and the reachability arithmetic are properties of the fixed-window search alone and
// run once.
static void runFrequencySections() {
    testRootRateUnityIsBitIdenticalToTheDirectRead();
    reportTransposedAt30Hz();
    reportFrequencySweep();
    reportFloorProbeMechanism();
    reportAlignableControls();
}

int main() {
    reportReachableInterval();
    reportReachabilityByFrequency();

    // The same measurements twice, from one binary, so the two columns differ in exactly one
    // thing. The FIXED-WINDOW pass reproduces the pre-PSOLA engine — it is the baseline every
    // number in the investigation was taken against.
    g_pitchSynchronous = false;
    std::printf("\n\n##################################################################\n");
    std::printf("### FIXED-WINDOW splices (no source period) — the prior behaviour ###\n");
    std::printf("##################################################################\n");
    runFrequencySections();

    g_pitchSynchronous = true;
    std::printf("\n\n##################################################################\n");
    std::printf("### PITCH-SYNCHRONOUS splices (detected source period)           ###\n");
    std::printf("##################################################################\n");
    runFrequencySections();

    // The window sweep is about what a LARGER WINDOW would buy, which was the alternative to
    // this track. Run under the shipped geometry only.
    reportWindowSweep();

    if (g_fail == 0) {
        std::printf("\nall preserve_low_frequency measurements completed\n");
        return 0;
    }
    std::printf("\n%d preserve_low_frequency check(s) failed\n", g_fail);
    return 1;
}