From ac653aa581237736b04e47f9a3099abf071b9428 Mon Sep 17 00:00:00 2001 From: daniel-c-harvey Date: Sat, 1 Aug 2026 21:11:08 -0400 Subject: [PATCH] Measure Preserve's splice-alignment geometry on low-frequency material MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A splice can only relocate by [0.75, 1.25]*window, so periods with no multiple in that interval never phase-align — at 50 ms, f < 16 Hz and 26.7-32 Hz. Harness runs by hand; too slow to gate. --- src/core/instrument/engine/CMakeLists.txt | 10 + src/core/instrument/engine/time_stretch.h | 15 + tests/test_preserve_low_frequency.cpp | 700 ++++++++++++++++++++++ 3 files changed, 725 insertions(+) create mode 100644 tests/test_preserve_low_frequency.cpp diff --git a/src/core/instrument/engine/CMakeLists.txt b/src/core/instrument/engine/CMakeLists.txt index 46690f9..864b9b9 100644 --- a/src/core/instrument/engine/CMakeLists.txt +++ b/src/core/instrument/engine/CMakeLists.txt @@ -50,6 +50,16 @@ reasampler_test(live_delivery LINK sampler_core) # both mode shapes share, and the Trigger tail's terminal behaviour. reasampler_test(staged_envelopes LINK sampler_core) +# Measurement harness for Preserve on low-frequency material: how the splice search's +# reachable relocation interval interacts with a long source period. Written longhand and +# deliberately NOT add_test()'d — it sweeps frequencies, windows and spectra and takes ~2m40s +# in Debug, which does not belong in a gate whose other targets run in seconds. It still +# builds with everything else, so it cannot rot into non-compilation. Run it by hand, in +# Release, when the question is what Preserve does to a given frequency. +add_executable(preserve_low_frequency_tests + ${REASAMPLER_TESTS_DIR}/test_preserve_low_frequency.cpp) +target_link_libraries(preserve_low_frequency_tests PRIVATE sampler_core) + # The Preserve read's source-feed schedule — the TIME half beside pitch_shift's PITCH half. # Header-only (it sits on the per-sample feed), hence INTERFACE. add_library(time_stretch INTERFACE) diff --git a/src/core/instrument/engine/time_stretch.h b/src/core/instrument/engine/time_stretch.h index d13ad5a..d2dd645 100644 --- a/src/core/instrument/engine/time_stretch.h +++ b/src/core/instrument/engine/time_stretch.h @@ -27,6 +27,21 @@ namespace reasampler::instrument::engine { // -24 st is reachable from the Pitch knob alone. (The pre-stretch rate-1.0 engine's floor by // the same inequality is P > 735, ~60 Hz — what this range raises the floor from, not what // it removes.) +// +// A SECOND, INDEPENDENT limit binds the same material, and no rate bound touches it. A splice +// relocates the tap by the nominal window refined by a search over +/- window/4, so the +// reachable relocation distances are exactly [0.75, 1.25] * window; a phase-aligned splice +// needs a WHOLE NUMBER of source periods inside that one interval. The interval is 0.5*window +// wide, so any period <= window/2 always has a multiple in it — but above that, coverage +// breaks into disjoint bands (n=1 covers periods [0.75, 1.25]*window, n=2 covers +// [0.375, 0.625]*window) and the gap between them is reachable by nothing. Because both the +// interval and the period scale with the sample rate, the unalignable set is fixed in Hz by +// the window's MILLISECONDS: at 50 ms that is f < 16 Hz and 26.7 Hz < f < 32 Hz. Measured +// (Release, 44.1k and 48k) at 30 Hz: the rendered pitch stays correct, but energy outside the +// fundamental is 3.6% at +2 st / rate 1.0 and 15.5% at rate 2.0, against 0.00% at 34 Hz under +// identical conditions; at 29 Hz / rate 2.0 the tone itself lands 7.4% flat. Unlike the +// cadence inequality above, this one is not about how OFTEN a splice fires — a window of at +// least two source periods removes it outright, and nothing else does. inline constexpr double kStretchRateMin = 0.5; inline constexpr double kStretchRateMax = 2.0; inline constexpr int kMaxFeedPerFrame = 2; // ceil(kStretchRateMax) diff --git a/tests/test_preserve_low_frequency.cpp b/tests/test_preserve_low_frequency.cpp new file mode 100644 index 0000000..7c02540 --- /dev/null +++ b/tests/test_preserve_low_frequency.cpp @@ -0,0 +1,700 @@ +// 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. +// +// 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/pitch_shift.h" +#include "../src/core/instrument/engine/time_stretch.h" +#include "../src/core/instrument/engine/voice.h" + +#include +#include +#include +#include +#include +#include + +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; + +// --------------------------------------------------------------------------------------- +// 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`). +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(sampleRate); + for (std::size_t i = 0; i < frames; ++i) { + s.frames[i] = static_cast(std::sin(w * static_cast(i) + phase)); + } + s.sampleRate = sampleRate; + s.rootNote = 60; + s.play.pitchEngine = engine; // Gate, no loop, default (fully open) AHDSR + 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 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 out(outFrames, 0.0); + for (std::size_t i = 0; i < outFrames; ++i) { + out[i] = static_cast(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& 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(i - prev); + prev = i; + ++count; + } + } + return count > 1 ? sum / static_cast(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& 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(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& v, std::size_t from, std::size_t to, + double cyclesPerFrame, double* worstStepCycles) { + const std::size_t blk = static_cast(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& v, std::size_t from, std::size_t to, + double cyclesPerFrame) { + const std::size_t blk = static_cast(1.0 / cyclesPerFrame); + std::vector 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(mid), r.end()); + return r[mid]; +} + +// Period of the highest normalized-autocorrelation peak over [minLag, maxLag] — 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. +static double autocorrPeriod(const std::vector& v, std::size_t from, std::size_t len, + std::int64_t minLag, std::int64_t maxLag) { + double e0 = 0.0; + for (std::size_t k = 0; k < len && from + k < v.size(); ++k) e0 += v[from + k] * v[from + k]; + if (e0 <= 0.0) return 0.0; + double best = -1e18; std::int64_t bestLag = 0; + std::vector score(static_cast(maxLag - minLag + 1), 0.0); + for (std::int64_t lag = minLag; lag <= maxLag; ++lag) { + double s = 0.0, e = 0.0; + for (std::size_t k = 0; k < len && from + k + static_cast(lag) < v.size(); + ++k) { + const double b = v[from + k + static_cast(lag)]; + s += v[from + k] * b; + e += b * b; + } + const double r = e > 0.0 ? s / std::sqrt(e0 * e) : 0.0; + score[static_cast(lag - minLag)] = r; + if (r > best) { best = r; bestLag = lag; } + } + // Parabolic refinement so the estimate isn't quantized to whole frames. + const std::size_t i = static_cast(bestLag - minLag); + double frac = 0.0; + if (i > 0 && i + 1 < score.size()) { + const double den = score[i - 1] - 2.0 * score[i] + score[i + 1]; + if (den < 0.0) frac = 0.5 * (score[i - 1] - score[i + 1]) / den; + } + return static_cast(bestLag) + frac; +} + +// 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& 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 mag(static_cast(kGrid), 0.0); + std::vector per(static_cast(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(g) / (kGrid - 1)); + per[static_cast(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(k) / + static_cast(len))); + const double x = v[from + k] * hann; + re += x * std::cos(w * static_cast(k)); + im += x * std::sin(w * static_cast(k)); + } + mag[static_cast(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> peaks; // (magnitude, period) + for (int g = 1; g + 1 < kGrid; ++g) { + const std::size_t i = static_cast(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& a, const std::pair& 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(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| == window on every splice (steady state) +}; + +// 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& src, std::int64_t window, + double rate, double shift, std::size_t outFrames, + std::vector* audio = nullptr) { + PitchShifter ps; + ps.configure(window); + ps.prime(src.data(), window); + ps.setShiftRatio(shift); + ps.setFeedRate(rate); + StretchCursor cur; + cur.start(window); + loop::ResolvedLoop lp{}; // inactive: the source is long enough to run straight through + + SpliceStats st; + 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(q) < src.size()) + ? src[static_cast(q)] : 0.0f; + fed = true; + } + const AudioSample o = fed ? ps.process(last) : ps.processNoInput(); + if (audio != nullptr) (*audio)[i] = static_cast(o); + const SpliceEvent& ev = ps.lastSplice(); + if (!ev.fired) continue; + ++st.count; + const double reloc = std::fabs(static_cast(ev.jump) - + static_cast(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) != window) st.jumpAlwaysNominal = false; + if (lastSpliceAt != 0) { intervalSum += static_cast(i - lastSpliceAt); ++intervals; } + lastSpliceAt = i; + } + st.meanInterval = intervals > 0 ? intervalSum / static_cast(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{44100, 2205}, + std::pair{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 noise(600000); + std::uint32_t rng = 12345u; + for (auto& x : noise) { + rng = rng * 1664525u + 1013904223u; + x = static_cast((static_cast(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(w), + 1000.0 * static_cast(w) / sr); + std::printf(" up-shift splices %lld, lag [%lld, %lld], reloc [%.2f, %.2f]\n", + up.count, static_cast(up.minLag), + static_cast(up.maxLag), up.minReloc, up.maxReloc); + std::printf(" down-shift splices %lld, lag [%lld, %lld], reloc [%.2f, %.2f]\n", + dn.count, static_cast(dn.minLag), + static_cast(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(w - w / 4), static_cast(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{44100, 2205}, + std::pair{48000, 2400}}) { + const int sr = g.first; + const std::int64_t w = g.second; + const double lo = static_cast(w - w / 4), hi = static_cast(w + w / 4); + std::printf(" %d Hz / window %lld, interval [%.0f, %.0f] frames:\n", sr, + static_cast(w), lo, hi); + for (double f : freqs) { + const double P = static_cast(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 p = renderVoice(pres, 60, 1.0, w, outFrames); + const std::vector 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 p48 = renderVoice(pres48, 60, 1.0, 2400, outFrames); + const std::vector 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 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 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(big), up[0], + static_cast(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(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 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(shortSample.frames[0])); + CHECK(shortOut[0] == static_cast(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 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(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 src(s.frames.begin(), s.frames.end()); + const SpliceStats st = spliceGeometry(src, window, rate, shift, outFrames); + const double lo = static_cast(window - window / 4); + const double hi = static_cast(window + window / 4); + int n = 0; + 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(st.count) : 0.0; + std::printf(" %-26s f=%6.1f Hz shift=%.4f rate=%.2f | 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 | " + "align %s%s\n", + label, freqHz, shift, rate, 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 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 out = renderVoice(s, note, rate, window, 300000); + const double wantPeriod = (static_cast(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 = autocorrPeriod(out, 60000, 60000, + std::max(40, + static_cast(wantPeriod * 0.5)), + static_cast(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(w - w / 4), hi = static_cast(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(sr) / period; + // The cadence inequality from time_stretch.h, evaluated for this row. + const double cadence = static_cast(w) / 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 src(srcLen); + for (std::size_t i = 0; i < srcLen; ++i) { + src[i] = static_cast( + std::sin(2.0 * kPi * static_cast(i) / period)); + } + std::vector out; + const SpliceStats st = spliceGeometry(src, w, rate, shift, 60000, &out); + const double want = period / shift; + const double zc = periodIn(out, 20000, 50000); + const double ac = autocorrPeriod(out, 20000, 20000, + static_cast(want * 0.5), + static_cast(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(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(w - w / 4), hi = static_cast(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(w), 1000.0 * static_cast(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(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); +} + +int main() { + reportReachableInterval(); + reportReachabilityByFrequency(); + testRootRateUnityIsBitIdenticalToTheDirectRead(); + reportTransposedAt30Hz(); + reportFrequencySweep(); + reportFloorProbeMechanism(); + reportAlignableControls(); + 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; +}