daniel/reasampler
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Γ-W1-T7: make Preserve's splices pitch-synchronous — the jump is a whole number of the source's own period, detected once at load

Browse Source 
30 Hz out-of-band energy 15.45% -> 0.00%; the 29 Hz rate-2.0 detune -133 -> +0 cents.
An unknown period keeps the fixed-window geometry bit for bit. The detector cannot
reach process(): sampler_core does not link it.
This commit is contained in:
Daniel Harvey
2026-08-01 23:19:16 -04:00
parent 163ab11e05
commit 93230208ff
17 changed files with 1085 additions and 137 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/energy_outside_fundamental.h
-50
View File
@@ -1,50 +0,0 @@
#pragma once
// Out-of-band spectral energy metric: the same period-grid, Hann-windowed direct-evaluation
// approach as test_preserve_low_frequency.cpp's reportSpectrum. Chosen over zero-crossing
// counting because splice debris adds spurious crossings that make that estimator
// anti-correlated with severity (a render can read a badly wrong PERIOD while this metric
// shows it is mostly clean, or vice versa). Grid/segment sizes are smaller than the hand-run
// harness's — this one runs inside the gated suite.
#include <cmath>
#include <cstddef>
#include <vector>
namespace reasampler::test_support {
// Percentage (0..100) of the segment [from, from+len)'s spectral energy that falls outside
// +/- 6% of `wantPeriod` (frames). 0 = a clean single tone at that period; higher values mean
// harmonics, splice-cadence sidebands, or crossfade cancellation debris are present.
inline double energyOutsideFundamentalPercent(const std::vector<double>& v, std::size_t from,
std::size_t len, double wantPeriod) {
constexpr double kPi = 3.14159265358979323846;
constexpr int kGrid = 400;
const double pLo = 30.0, pHi = 8000.0;
std::vector<double> mag(static_cast<std::size_t>(kGrid));
std::vector<double> per(static_cast<std::size_t>(kGrid));
for (int g = 0; g < kGrid; ++g) {
// Geometric grid: constant relative resolution across the swept period range.
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 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;
}
return eTotal > 0.0 ? 100.0 * (1.0 - eFund / eTotal) : 0.0;
}
} // namespace reasampler::test_support
tests/test_period_detect.cpp
+266
View File
@@ -0,0 +1,266 @@
// Standalone tests for reasampler::instrument::engine::detectPeriod — the offline source-period
// estimate behind Preserve's pitch-synchronous splices. No VST3, no REAPER, no test framework.
//
// Covers:
// 1. accuracy on pure tones across the searched band, at 44.1k and 48k, including the
// non-integer periods every real capture actually has — the splice jump is n periods, so
// a fractional-frame error lands multiplied by n.
// 2. the fundamental, not a harmonic: a sawtooth and a missing-fundamental stack must both
// report the repeat period, which is what a splice has to align on.
// 3. graceful degradation — noise, silence, and a source whose period changes mid-sample all
// return NONE. That is the contract the shifter's fixed-window fallback rests on: an
// estimate that is merely wrong would misalign every splice, which is worse than none.
// 4. the band edges and the short-sample path.
// 5. what the load pays, and that it does not grow with the sample length.
#include "../src/core/instrument/engine/period_detect.h"
#include <chrono>
#include <cmath>
#include <cstdint>
#include <cstdio>
#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;
static std::vector<AudioSample> sineOfPeriod(std::size_t frames, double period,
double phase = 0.0) {
std::vector<AudioSample> s(frames);
for (std::size_t i = 0; i < frames; ++i) {
s[i] = static_cast<AudioSample>(
std::sin(2.0 * kPi * static_cast<double>(i) / period + phase));
}
return s;
}
// --- 1. Accuracy on pure tones -------------------------------------------------------------
static void testPureTonePeriodIsFoundToBetterThanATenthOfAFrame() {
// Deliberately non-integer periods: an integer-only estimator passes an integer-period
// sweep and still misaligns every real capture.
const double periods[] = {23.7, 50.0, 100.25, 200.45, 441.0, 999.9, 1470.0, 2000.3, 2756.0};
for (double p : periods) {
const std::vector<AudioSample> src = sineOfPeriod(120000, p);
const PeriodEstimate est = detectPeriod(src, 44100);
CHECK(est.valid());
if (!est.valid()) {
std::printf(" period %.2f: NOT DETECTED\n", p);
continue;
}
const double errFrames = std::fabs(est.frames - p);
std::printf(" period %8.2f -> %8.4f (err %.4f fr, conf %.3f)\n", p, est.frames,
errFrames, est.confidence);
CHECK(errFrames < 0.1);
CHECK(est.confidence > 0.8);
}
}
static void testTheEstimateIsInSourceFramesSoTheRateOnlyMovesTheBand() {
// The same 30 Hz tone at two rates: the answer is frames, so it must track the rate. This
// is what lets the shifter compare it against a window that is also in frames.
for (int rate : {44100, 48000}) {
const double p = static_cast<double>(rate) / 30.0;
const std::vector<AudioSample> src = sineOfPeriod(160000, p);
const PeriodEstimate est = detectPeriod(src, rate);
CHECK(est.valid());
if (est.valid()) {
std::printf(" 30 Hz @ %d: %.3f fr (want %.3f)\n", rate, est.frames, p);
CHECK(std::fabs(est.frames - p) < 0.5);
}
}
}
// --- 2. The fundamental, not a harmonic ----------------------------------------------------
static void testHarmonicRichSourceReportsTheRepeatPeriodNotAPartial() {
// A sawtooth's strongest correlation dips at EVERY multiple of its period; a global-minimum
// estimator picks 2P or 3P about as often as P. YIN's first-dip rule is what makes this
// pass, and a jump quantized to 2P would splice a whole cycle out of phase half the time.
const double p = 512.0;
std::vector<AudioSample> src(120000);
for (std::size_t i = 0; i < src.size(); ++i) {
double v = 0.0;
for (int h = 1; h <= 12; ++h) {
v += std::sin(2.0 * kPi * h * static_cast<double>(i) / p) / h;
}
src[i] = static_cast<AudioSample>(0.5 * v);
}
const PeriodEstimate est = detectPeriod(src, 44100);
CHECK(est.valid());
if (est.valid()) {
std::printf(" sawtooth P=512 -> %.3f\n", est.frames);
CHECK(std::fabs(est.frames - p) < 1.0);
}
}
static void testMissingFundamentalStillReportsTheRepeatPeriod() {
// Partials 2..6 of a 700-frame period: there is no energy AT the fundamental, but the
// waveform still repeats every 700 frames — and repetition, not spectral content, is what
// a splice has to land on.
const double p = 700.0;
std::vector<AudioSample> src(120000);
for (std::size_t i = 0; i < src.size(); ++i) {
double v = 0.0;
for (int h = 2; h <= 6; ++h) {
v += std::sin(2.0 * kPi * h * static_cast<double>(i) / p);
}
src[i] = static_cast<AudioSample>(0.2 * v);
}
const PeriodEstimate est = detectPeriod(src, 44100);
CHECK(est.valid());
if (est.valid()) {
std::printf(" missing fundamental P=700 -> %.3f\n", est.frames);
CHECK(std::fabs(est.frames - p) < 2.0);
}
}
// --- 3. Graceful degradation ---------------------------------------------------------------
static void testNoiseSilenceAndAPeriodChangeAllReportNone() {
// White noise: no dip below the absolute threshold anywhere.
{
std::vector<AudioSample> src(120000);
std::uint32_t rng = 22222u;
for (auto& x : src) {
rng = rng * 1664525u + 1013904223u;
x = static_cast<AudioSample>((static_cast<double>(rng >> 8) / 8388608.0) - 1.0);
}
const PeriodEstimate est = detectPeriod(src, 44100);
std::printf(" white noise -> %s (%.3f)\n", est.valid() ? "DETECTED" : "none",
est.frames);
CHECK(!est.valid());
}
// Digital silence: the difference function is degenerate, not merely inconclusive.
{
const std::vector<AudioSample> src(120000, 0.0f);
CHECK(!detectPeriod(src, 44100).valid());
}
// Two halves at genuinely different periods: the probes disagree, so there is no ONE
// period, and reporting either half's would misalign every splice in the other half.
{
std::vector<AudioSample> src(160000);
double phase = 0.0;
for (std::size_t i = 0; i < src.size(); ++i) {
phase += 2.0 * kPi / (i < 80000 ? 300.0 : 700.0);
src[i] = static_cast<AudioSample>(std::sin(phase));
}
const PeriodEstimate est = detectPeriod(src, 44100);
std::printf(" period change 300->700 -> %s (%.3f)\n", est.valid() ? "DETECTED" : "none",
est.frames);
CHECK(!est.valid());
}
// Degenerate inputs.
CHECK(!detectPeriod({}, 44100).valid());
CHECK(!detectPeriod(sineOfPeriod(120000, 441.0), 0).valid());
}
static void testAPercussiveDecayIsNotForcedIntoAPeriod() {
// Filtered noise with a fast decay — the shape of a one-shot drum hit. Nothing repeats, so
// the answer must be none rather than whatever the envelope's own length looks like.
std::vector<AudioSample> src(120000);
std::uint32_t rng = 909090u;
double lp = 0.0;
for (std::size_t i = 0; i < src.size(); ++i) {
rng = rng * 1664525u + 1013904223u;
const double n = (static_cast<double>(rng >> 8) / 8388608.0) - 1.0;
lp += 0.25 * (n - lp);
const double env = std::exp(-static_cast<double>(i % 22050) / 2000.0);
src[i] = static_cast<AudioSample>(lp * env);
}
const PeriodEstimate est = detectPeriod(src, 44100);
std::printf(" percussive decay -> %s (%.3f)\n", est.valid() ? "DETECTED" : "none",
est.frames);
CHECK(!est.valid());
}
// --- 4. Band edges and short sources -------------------------------------------------------
static void testBelowTheBandReportsNoneAndAboveItReportsAWholeMultiple() {
// Below kPeriodDetectMinHz: none. This is the load-bearing edge — such a period cannot fit
// the splice jump anyway, so an answer here would only be one the shifter must reject.
const std::vector<AudioSample> low = sineOfPeriod(200000, 44100.0 / 8.0); // 8 Hz
CHECK(!detectPeriod(low, 44100).valid());
// Above kPeriodDetectMaxHz the search floor sits well above the true period, so what comes
// back is a WHOLE MULTIPLE of it — which is still an exactly aligned splice target, since
// every multiple of a period is a period. That is why the high edge needs no special
// handling: being outside the band costs nothing, because alignment was never in question
// for a tone this short-period.
const double p = 44100.0 / 6000.0; // 7.35 frames
const PeriodEstimate high = detectPeriod(sineOfPeriod(120000, p), 44100);
std::printf(" 6 kHz (P=%.3f) -> %s (%.3f, = %.3f periods)\n", p,
high.valid() ? "detected" : "none", high.frames, high.frames / p);
if (high.valid()) {
const double n = high.frames / p;
CHECK(std::fabs(n - std::floor(n + 0.5)) < 0.02);
}
}
static void testAShortSourceShortensTheSearchRatherThanRefusing() {
// A 12000-frame one-shot cannot host a full-band probe; the search band shortens to fit and
// a 200-frame period is still found. Below that the answer is none, not a guess.
const std::vector<AudioSample> shortSrc = sineOfPeriod(12000, 200.0);
const PeriodEstimate est = detectPeriod(shortSrc, 44100);
std::printf(" 12000-frame source, P=200 -> %s (%.3f)\n", est.valid() ? "detected" : "none",
est.frames);
CHECK(est.valid());
if (est.valid()) CHECK(std::fabs(est.frames - 200.0) < 0.5);
// Too short for even the minimum lag: none.
CHECK(!detectPeriod(sineOfPeriod(40, 20.0), 44100).valid());
}
// --- 5. What the load pays ------------------------------------------------------------------
// The whole reason a detector is affordable in a sampler is that it runs ONCE, off the audio
// thread, on a source that is already fully known. This prints what that once costs, and
// asserts the property that makes it safe: the cost does NOT grow with the sample length —
// a fixed number of fixed-size probes is analysed however long the capture is. Meaningful
// only in a Release build; asserted as a RATIO so it holds at either optimization level.
static void testDetectionCostIsBoundedRegardlessOfSampleLength() {
double shortMs = 0.0, longMs = 0.0;
for (std::size_t frames : {std::size_t{220500}, std::size_t{4410000}}) { // 5 s and 100 s
const std::vector<AudioSample> src = sineOfPeriod(frames, 441.0);
const int reps = 5;
const auto t0 = std::chrono::steady_clock::now();
double guard = 0.0;
for (int r = 0; r < reps; ++r) guard += detectPeriod(src, 44100).frames;
const double ms =
1000.0 * std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count()
/ reps;
CHECK(guard > 0.0);
std::printf(" [measure] detectPeriod over %7.1f s of source: %.3f ms\n",
static_cast<double>(frames) / 44100.0, ms);
(frames == 220500 ? shortMs : longMs) = ms;
}
// 20x the source for well under 2x the cost — the probes are fixed-size and fixed in
// number, so the only length dependence left is the cache behaviour of reaching further
// into the buffer.
CHECK(longMs < shortMs * 2.0 + 0.5);
}
int main() {
testPureTonePeriodIsFoundToBetterThanATenthOfAFrame();
testTheEstimateIsInSourceFramesSoTheRateOnlyMovesTheBand();
testHarmonicRichSourceReportsTheRepeatPeriodNotAPartial();
testMissingFundamentalStillReportsTheRepeatPeriod();
testNoiseSilenceAndAPeriodChangeAllReportNone();
testAPercussiveDecayIsNotForcedIntoAPeriod();
testBelowTheBandReportsNoneAndAboveItReportsAWholeMultiple();
testAShortSourceShortensTheSearchRatherThanRefusing();
testDetectionCostIsBoundedRegardlessOfSampleLength();
if (g_fail == 0) {
std::printf("all period_detect tests passed\n");
return 0;
}
std::printf("%d period_detect check(s) failed\n", g_fail);
return 1;
}
tests/test_pitch_shift.cpp
+228 -8
View File
@@ -30,11 +30,17 @@
// 8. stereo linked lag (Q-W0 T1-01) — a follower channel driven via processLinked() mirrors
// the master's splice decision (jump/lag/frac/fadeLen AND firing frame) exactly, on
// decorrelated stereo content where an independent per-channel search provably diverges.
// 10. pitch-synchronous splices — the nominal jump snapped to a whole number of source
// periods: the jump law itself, the bit-identical unknown-period fallback, 30 Hz and
// 29 Hz (the two symptoms of the unalignable gap), and the cadence corner, which this
// leaves where it found it.
#include "../src/core/instrument/engine/pitch_shift.h"
#include "energy_outside_fundamental.h"
#include "tone_metrics.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <vector>
@@ -595,14 +601,17 @@ static void testStereoLinkedLagSharedSchedule() {
// source frame due on an output frame go through writeFrame (no output), the last through
// process(); an output frame with none due takes processNoInput(). Returns the output plus,
// via `consumed`, how much source it ate.
// `sourcePeriod` > 0 puts the shifter on the pitch-synchronous jump the loader would have
// given it; 0 (the default) is the fixed-window fallback every pre-PSOLA call here exercises.
static std::vector<double> runStretch(const std::vector<AudioSample>& src, std::int64_t w,
double feedRate, double shift, std::size_t outFrames,
std::size_t* consumed) {
std::size_t* consumed, double sourcePeriod = 0.0) {
PitchShifter ps;
ps.configure(w);
ps.prime(src.data(), w);
ps.setShiftRatio(shift);
ps.setFeedRate(feedRate);
ps.setSourcePeriod(sourcePeriod);
std::size_t pos = static_cast<std::size_t>(w);
double debt = 0.0;
std::vector<double> out(outFrames);
@@ -746,12 +755,23 @@ static void testStretchAndShiftComposeSafely() {
// observable here: an investigation (test_preserve_low_frequency.cpp) found the P=500
// render's FUNDAMENTAL within 0.03% of target by autocorrelation and spectral peak alike,
// while the zero-crossing estimator read 23% flat — splice debris adds spurious crossings
// the count cannot tell from a real detune. Energy outside the fundamental tracks the actual
// damage instead: measured here (same rate/shift/source, this file's own metric parameters)
// at 7.31% / 14.41% / 21.22% for P=500/600/700, against 0.10% on an alignable control (P=200,
// below the safe floor) at the same rate and shift — so that is what this asserts: a known,
// characterised property of the range, not a pass/fail on a period estimate. A failure on
// either bound below is a finding — report it, don't retune the thresholds to hide it.
// the count cannot tell from a real detune. Energy outside the fundamental is measured here
// (same rate/shift/source, this file's own metric parameters) at 7.31% / 14.41% / 21.22% for
// P=500/600/700, against 0.10% on an alignable control (P=200, below the safe floor) at the
// same rate and shift. Those readings are stable and are what the bounds below hold.
//
// **CORRECTED — what those three readings MEAN.** They were once read as the corner's damage.
// They are almost entirely the metric's own floor: an ideal tone at the same want-period,
// measured identically, reads 7.08% / 13.90% / 20.92% (idealToneFloorPercent, below), because
// a long period under a 32768-frame segment leaks part of its own mainlobe outside the +/-6%
// band. The corner's real EXCESS over that floor is 0.23% / 0.51% / 0.30% — small, real, and
// nothing like the headline numbers. The alignable control's 0.10% is genuinely near-zero only
// because its want-period is short enough to have almost no floor. The bounds below are kept
// as a stable regression tripwire on the raw readings; read the excess, not the reading.
//
// This measures the FIXED-WINDOW path — no source period is set, which is what a capture with
// no single period (percussive, polyphonic, noise) gets. What the same corner does once the
// splice is pitch-synchronous is the test immediately after this one.
static void testStretchCadenceCornerArtifactEnergyAtRate2ShiftQuarter() {
using reasampler::test_support::energyOutsideFundamentalPercent;
const std::int64_t w = 2205;
@@ -804,6 +824,201 @@ static void testStretchCadenceCornerArtifactEnergyAtRate2ShiftQuarter() {
}
}
// --- 10. Pitch-synchronous splices: the nominal jump is a whole number of SOURCE periods. ---
// The out-of-band metric's OWN floor at a given period: a mathematically perfect tone, measured
// with exactly the parameters a render is. A long period under a fixed segment leaks part of
// its own mainlobe outside the +/-6% band, and that leakage grows steeply with the period — so
// a raw reading at period 2800 is not comparable with one at period 800, and neither is
// comparable with zero. The EXCESS over this floor is the honest "how much of this render is
// not the tone" number.
static double idealToneFloorPercent(double wantPeriod, std::size_t from, std::size_t len) {
std::vector<double> v(from + len + 2);
for (std::size_t i = 0; i < v.size(); ++i) {
v[i] = std::sin(2.0 * kPi * static_cast<double>(i) / wantPeriod);
}
return reasampler::test_support::energyOutsideFundamentalPercent(v, from, len, wantPeriod);
}
// A splice can only phase-align on a landing point a whole number of source periods away, and
// the search only reaches [0.75, 1.25] windows. periodAlignedJump is what puts an aligned point
// at the CENTRE of that interval instead of hoping one falls inside it.
static void testPeriodAlignedJumpSnapsToWholePeriodsWithinTheReachableBound() {
const std::int64_t w = 2205; // the product window at 44.1k
const std::int64_t maxJump = 2756; // 1.25 * w, the shifter's own jumpMax_
// Unknown period, and a period too long for even ONE whole one to fit: the fixed window,
// unchanged. Both are the documented fallback, and both must be EXACTLY today's geometry.
CHECK(periodAlignedJump(w, maxJump, 0.0) == w);
CHECK(periodAlignedJump(w, maxJump, -5.0) == w);
CHECK(periodAlignedJump(w, maxJump, 3000.0) == w);
// 30 Hz at 44.1k (P = 1470): two periods overshoot the bound, so it takes ONE — which is
// the case the whole track exists for. The pre-PSOLA geometry could reach neither 1470 nor
// 2940 from a 2205 nominal, since the search only spans [1654, 2756].
CHECK(periodAlignedJump(w, maxJump, 1470.0) == 1470);
CHECK(2 * 1470 > maxJump); // the witness that one period is forced, not merely chosen
// 220 Hz (P = 200.4545): eleven periods land within a frame of the window itself, so the
// splice cadence is essentially untouched while every landing is aligned.
CHECK(periodAlignedJump(w, maxJump, 44100.0 / 220.0) == 2205);
// A period just under the bound is taken whole; the result is never over the bound, at any
// period in the band. Sweeping is what proves the shrink loop terminates correctly rather
// than one hand-picked value doing so.
for (double p = 20.0; p < 3200.0; p += 0.37) {
const std::int64_t j = periodAlignedJump(w, maxJump, p);
CHECK(j >= 1);
if (p > static_cast<double>(maxJump)) {
CHECK(j == w); // out of reach -> fallback
} else {
CHECK(j <= maxJump);
// Aligned: the jump is a whole number of periods, to within the rounding to frames.
const double n = static_cast<double>(j) / p;
CHECK(std::fabs(n - std::floor(n + 0.5)) * p < 0.51);
}
}
}
// setSourcePeriod(0) and never calling it are the same state, not merely similar ones — the
// cheap half of the fallback claim. The EXPENSIVE half, that the fallback still matches the
// engine as it shipped, is testPreserveUnityRateIsBitIdenticalToTheShippedRead in
// test_sampler_core.cpp: it hashes four rendered streams (including transposed ones that
// really splice) against a baseline captured from commit 0a7778b, and it passes unmodified.
static void testAnUnknownPeriodIsBitIdenticalToTheFixedWindowGeometry() {
const std::int64_t w = 2205;
const std::vector<AudioSample> src = sine(400000, 400000.0 / 196.37);
const std::vector<double> never = runStretch(src, w, 1.0, 1.5, 40000, nullptr);
const std::vector<double> zeroed = runStretch(src, w, 1.0, 1.5, 40000, nullptr, 0.0);
bool same = true;
for (std::size_t i = 0; i < never.size(); ++i) if (never[i] != zeroed[i]) same = false;
CHECK(same);
}
// THE case this track exists for. 30 Hz sits in the only unalignable gap above 16 Hz at the
// product's 50 ms window: its nearest whole multiple misses the reachable interval by 184
// frames (45 degrees of phase), and the investigation measured the resulting sidebands at
// 3.57% out-of-band at +2 st / rate 1.0 and 15.45% at rate 2.0, against 0.00% on alignable
// controls. Here the same two conditions run with and without the source period, against a
// 34 Hz control that was alignable all along.
//
// The absolute numbers are NOT the harness's: a 1470-frame period under a 32768-frame segment
// leaks part of its own mainlobe outside the +/-6% band, so every reading here carries the same
// floor. That is exactly why the control is measured at the same length — the assertion is that
// 30 Hz reaches the control's floor, not that it reaches zero.
static void testThirtyHertzSplicesAlignOnceTheSourcePeriodIsKnown() {
using reasampler::test_support::energyOutsideFundamentalPercent;
const std::int64_t w = 2205;
const std::size_t srcLen = 400000, outFrames = 60000, from = 20000, len = 32768;
struct Row { const char* label; double freq; double rate; double semis; };
const Row rows[] = {
{"30 Hz +2 st, rate 1.0", 30.0, 1.0, 2.0},
{"30 Hz rate 2.0", 30.0, 2.0, 0.0},
{"34 Hz +2 st, rate 1.0", 34.0, 1.0, 2.0}, // control: alignable without a period
{"34 Hz rate 2.0", 34.0, 2.0, 0.0},
};
double controlWorst = 0.0, subjectWorst = 0.0;
for (const Row& r : rows) {
const double period = 44100.0 / r.freq;
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));
}
const double shift = std::pow(2.0, r.semis / 12.0);
const double want = period / shift;
const std::vector<double> off = runStretch(src, w, r.rate, shift, outFrames, nullptr);
const std::vector<double> on =
runStretch(src, w, r.rate, shift, outFrames, nullptr, period);
for (double v : on) CHECK(std::isfinite(v));
const double floor = idealToneFloorPercent(want, from, len);
const double pctOff = energyOutsideFundamentalPercent(off, from, len, want) - floor;
const double pctOn = energyOutsideFundamentalPercent(on, from, len, want) - floor;
std::printf(" [30 Hz] %-24s (want %6.1f fr, metric floor %.2f%%) excess energy: "
"fixed window %6.2f%% -> pitch-synchronous %6.2f%%\n", r.label, want, floor,
pctOff, pctOn);
if (r.freq == 34.0) controlWorst = std::max(controlWorst, pctOn);
else subjectWorst = std::max(subjectWorst, pctOn);
}
// 30 Hz stops being a special case: with the period known its excess over the metric's own
// floor is no worse than the alignable neighbour's, measured identically. Against the
// control rather than against a fixed number, so the assertion cannot be satisfied by a
// change that merely raised the floor everywhere.
std::printf(" [30 Hz] worst subject excess %.2f%% vs worst control excess %.2f%%\n",
subjectWorst, controlWorst);
CHECK(subjectWorst < 0.10);
CHECK(subjectWorst <= controlWorst + 0.05); // 0.05 absorbs the floor subtraction's sign noise
}
// The sharpest single symptom of the geometry: at 29 Hz the nearest multiple misses the
// reachable interval ONE-SIDED rather than straddling, so the per-splice phase steps stop
// cancelling and accumulate into a real detune — the investigation measured -133 cents at
// rate 2.0 with NO transposition at all. Rate moves duration; it must not move pitch.
static void testTwentyNineHertzAtRateTwoKeepsItsPitch() {
using reasampler::test_support::autocorrelationPeriod;
const std::int64_t w = 2205;
const double period = 44100.0 / 29.0; // 1520.7 frames
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));
}
auto centsOf = [&](double sourcePeriod) {
const std::vector<double> out =
runStretch(src, w, /*rate=*/2.0, /*shift=*/1.0, 60000, nullptr, sourcePeriod);
const double got = autocorrelationPeriod(out, 20000, 20000,
static_cast<std::int64_t>(period * 0.5),
static_cast<std::int64_t>(period * 1.7));
return 1200.0 * std::log2(got / period);
};
const double centsOff = centsOf(0.0);
const double centsOn = centsOf(period);
std::printf(" [29 Hz] rate 2.0, no transposition: fixed window %+.1f cents -> "
"pitch-synchronous %+.1f cents\n", centsOff, centsOn);
CHECK(std::fabs(centsOn) < 10.0);
// The fixed-window reading is asserted too, and that is what makes the pair non-vacuous: a
// setSourcePeriod that silently did nothing would render both identically and fail here.
CHECK(std::fabs(centsOff) > 50.0);
}
// The cadence corner (rate 2.0, -24 st, source periods above the 315-frame safe floor) is the
// OTHER mechanism — a splice landing inside a single perceived cycle. Measured against the
// metric's own floor, PSOLA moves it by nothing: 0.23/0.51/0.30% excess becomes 0.24/0.49/0.30%.
//
// That is not a shortfall, it is what the corner turned out to be. Correcting the previous
// test's reading (see its comment) shrank the corner from a 7-21% headline to a sub-1% excess,
// which leaves PSOLA nothing to recover there — a pitch-synchronous jump makes each splice
// land in phase, and these splices already did; what it cannot do is make them less frequent.
// So this asserts NO REGRESSION, not an improvement, and says so rather than claiming one.
static void testCadenceCornerIsUnmovedByAPitchSynchronousSplice() {
using reasampler::test_support::energyOutsideFundamentalPercent;
const std::int64_t w = 2205;
const double shift = std::pow(2.0, -24.0 / 12.0);
const std::size_t outFrames = 60000, from = 20000, len = 32768;
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));
}
const std::vector<double> off = runStretch(src, w, 2.0, shift, outFrames, nullptr);
const std::vector<double> on =
runStretch(src, w, 2.0, shift, outFrames, nullptr, period);
for (double v : on) CHECK(std::isfinite(v));
const double want = period / shift;
const double floor = idealToneFloorPercent(want, from, len);
const double pctOff = energyOutsideFundamentalPercent(off, from, len, want);
const double pctOn = energyOutsideFundamentalPercent(on, from, len, want);
std::printf(" [cadence corner, PSOLA] period %.0f (want %.0f, metric floor %.2f%%): "
"excess %.2f%% -> %.2f%%\n", period, want, floor, pctOff - floor,
pctOn - floor);
CHECK(pctOn - floor < 1.0); // the corner's real excess, PSOLA or not
CHECK(pctOn < pctOff + 0.05); // and PSOLA costs it nothing
}
}
// The two new entry points on a shifter that was never configured (a Varispeed voice's) —
// neither may touch the empty ring.
static void testStretchEntryPointsOnPassThrough() {
@@ -826,6 +1041,11 @@ int main() {
testStretchMovesDurationNotPitch();
testStretchAndShiftComposeSafely();
testStretchCadenceCornerArtifactEnergyAtRate2ShiftQuarter();
testPeriodAlignedJumpSnapsToWholePeriodsWithinTheReachableBound();
testAnUnknownPeriodIsBitIdenticalToTheFixedWindowGeometry();
testThirtyHertzSplicesAlignOnceTheSourcePeriodIsKnown();
testTwentyNineHertzAtRateTwoKeepsItsPitch();
testCadenceCornerIsUnmovedByAPitchSynchronousSplice();
testStretchEntryPointsOnPassThrough();
if (g_fail == 0) {
tests/test_preserve_low_frequency.cpp
+72 -48
View File
@@ -14,6 +14,11 @@
// 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.
@@ -24,9 +29,11 @@
// 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>
@@ -44,12 +51,20 @@ static int g_fail = 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`).
// 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;
@@ -61,6 +76,7 @@ static SampleData sineSample(double freqHz, int sampleRate, std::size_t frames,
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;
}
@@ -162,37 +178,10 @@ static double medianResidual(const std::vector<double>& v, std::size_t from, std
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<double>& 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<double> score(static_cast<std::size_t>(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<std::size_t>(lag) < v.size();
++k) {
const double b = v[from + k + static_cast<std::size_t>(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<std::size_t>(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<std::size_t>(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<double>(bestLag) + frac;
}
// 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
@@ -264,24 +253,28 @@ struct SpliceStats {
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)
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) {
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;
@@ -308,7 +301,7 @@ static SpliceStats spliceGeometry(const std::vector<AudioSample>& src, std::int6
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 (std::llabs(ev.jump) != st.nominalJump) st.jumpAlwaysNominal = false;
if (lastSpliceAt != 0) { intervalSum += static_cast<double>(i - lastSpliceAt); ++intervals; }
lastSpliceAt = i;
}
@@ -503,24 +496,29 @@ static void measureRow(const char* label, double freqHz, int sr, std::int64_t wi
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);
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 | 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,
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 reachable interval");
reach ? "" : " <-- no whole period in the fixed-window reachable interval");
CHECK(finite);
}
@@ -536,7 +534,7 @@ static void deepDive(const char* label, double freqHz, int sr, std::int64_t wind
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,
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));
@@ -617,10 +615,11 @@ static void reportFloorProbeMechanism() {
std::sin(2.0 * kPi * static_cast<double>(i) / period));
}
std::vector<double> out;
const SpliceStats st = spliceGeometry(src, w, rate, shift, 60000, &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 = autocorrPeriod(out, 20000, 20000,
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%%) "
@@ -681,14 +680,39 @@ static void reportAlignableControls() {
deepDive("220 Hz rate 2.0", 220.0, 44100, 2205, 60, 2.0);
}
int main() {
reportReachableInterval();
reportReachabilityByFrequency();
// 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) {
tests/test_sampler_core.cpp
+27 -6
View File
@@ -3214,7 +3214,10 @@ static void testPreserveStretchThirtyTwoVoicesHoldUp() {
const std::size_t blockFrames = 44100; // one second of audio
const std::size_t voiceCount = 32;
const int kWarmupReps = 2;
const int kTimedReps = 7;
// 5 rather than 7: the source-period rows below doubled the row count, and run-to-run
// spread on this machine is ~5% either way, so the extra reps bought precision the number
// does not carry while costing the Debug gate real seconds.
const int kTimedReps = 5;
SampleData s = stretchProbeSample(200000, true);
s.loop.hasLoop = true; // held notes: all 32 sound for the whole run
s.loop.start = 40000;
@@ -3223,7 +3226,24 @@ static void testPreserveStretchThirtyTwoVoicesHoldUp() {
// 1.0 is the reference: it is the cost the shipped Preserve read already carries, so the
// two stretched rows are read as a delta against it rather than in isolation.
for (double rate : {1.0, 0.5, 2.0}) {
//
// The last two rows carry a SOURCE PERIOD, which is where the pitch-synchronous splice can
// cost something: the jump becomes a whole number of periods, and when that is shorter than
// the window the splice cadence rises in proportion — more correlation searches per second.
// 1470 (30 Hz at 44.1k) is the worst realistic case in the audible band, forcing a jump of
// 2/3 the window and therefore 1.5x the searches. 220.5 (200 Hz) is the typical one: ten
// periods land exactly on the window, so the cadence is unchanged and the row should read
// as the no-period one — which is the measurement that separates "the mechanism costs
// something" from "a shorter jump costs something". There is no per-frame cost either way:
// the jump is resolved once at note-on.
struct Row { double rate; double period; };
const Row rows[] = {
{1.0, 0.0}, {0.5, 0.0}, {2.0, 0.0},
{1.0, 220.5}, {1.0, 1470.0}, {2.0, 1470.0},
};
for (const Row& row : rows) {
const double rate = row.rate;
s.sourcePeriodFrames = row.period;
std::vector<double> nsPerVoiceFrame;
nsPerVoiceFrame.reserve(kTimedReps);
for (int rep = 0; rep < kWarmupReps + kTimedReps; ++rep) {
@@ -3263,10 +3283,11 @@ static void testPreserveStretchThirtyTwoVoicesHoldUp() {
const double medianNs = nsPerVoiceFrame[nsPerVoiceFrame.size() / 2];
const double secsAtMedian =
medianNs * static_cast<double>(blockFrames) * static_cast<double>(voiceCount) / 1e9;
std::printf(" [measure] 32 stereo Preserve voices @ rate %.2f: median %.1f ns/voice/"
"frame [%.1f .. %.1f] over %d reps (%.1f%% of realtime at the median)\n",
rate, medianNs, nsPerVoiceFrame.front(), nsPerVoiceFrame.back(), kTimedReps,
100.0 * secsAtMedian);
std::printf(" [measure] 32 stereo Preserve voices @ rate %.2f, source period %6.1f: "
"median %.1f ns/voice/frame [%.1f .. %.1f] over %d reps (%.1f%% of realtime "
"at the median)\n",
rate, row.period, medianNs, nsPerVoiceFrame.front(), nsPerVoiceFrame.back(),
kTimedReps, 100.0 * secsAtMedian);
}
}
tests/tone_metrics.h
+93
View File
@@ -0,0 +1,93 @@
#pragma once
// The two tone metrics the Preserve tests and the hand-run low-frequency harness share, so
// there is one of each rather than a copy per file.
//
// Zero-crossing counting is deliberately NOT among them: splice debris adds spurious crossings
// that make that estimator anti-correlated with severity (a render can read a badly wrong
// PERIOD while it is spectrally clean, or vice versa).
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <vector>
namespace reasampler::test_support {
// Percentage (0..100) of the segment [from, from+len)'s spectral energy that falls outside
// +/- 6% of `wantPeriod` (frames), evaluated directly on a geometric period grid (no FFT-bin
// quantization). 0 = a clean single tone at that period; higher values mean harmonics,
// splice-cadence sidebands, or crossfade cancellation debris are present.
//
// `grid` trades resolution for cost: the gated tests run the default, the hand-run harness
// raises it. The value is NOT comparable across grid sizes or segment lengths — a long period
// under a short segment leaks part of its own mainlobe outside the +/-6% band, so a reading is
// only meaningful against a control measured at the SAME len and grid.
inline double energyOutsideFundamentalPercent(const std::vector<double>& v, std::size_t from,
std::size_t len, double wantPeriod,
int grid = 400) {
constexpr double kPi = 3.14159265358979323846;
const int kGrid = grid;
const double pLo = 30.0, pHi = 8000.0;
std::vector<double> mag(static_cast<std::size_t>(kGrid));
std::vector<double> per(static_cast<std::size_t>(kGrid));
for (int g = 0; g < kGrid; ++g) {
// Geometric grid: constant relative resolution across the swept period range.
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 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;
}
return eTotal > 0.0 ? 100.0 * (1.0 - eFund / eTotal) : 0.0;
}
// Period (frames) of the highest normalized-autocorrelation peak over [minLag, maxLag], with a
// parabolic refinement so the answer is not quantized to whole frames. Bracket the caller's
// range to roughly [0.5, 1.7] x the expected period: a pure tone autocorrelates equally at
// EVERY multiple of its period, so an unbounded search reports 2P about half the time.
inline double autocorrelationPeriod(const std::vector<double>& v, std::size_t from,
std::size_t len, std::int64_t minLag, std::int64_t maxLag) {
if (maxLag <= minLag) return 0.0;
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;
std::vector<double> score(static_cast<std::size_t>(maxLag - minLag + 1), 0.0);
double best = -1e18;
std::int64_t bestLag = minLag;
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<std::size_t>(lag) < v.size();
++k) {
const double b = v[from + k + static_cast<std::size_t>(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<std::size_t>(lag - minLag)] = r;
if (r > best) { best = r; bestLag = lag; }
}
const std::size_t i = static_cast<std::size_t>(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<double>(bestLag) + frac;
}
} // namespace reasampler::test_support