Skip to content

Bound the psola/pvoc sample clocks so 32-bit targets never overflow (Stage 0) #20

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
tap merged 2 commits into from
Sep 17, 2026
+203 −21
Merged

Bound the psola/pvoc sample clocks so 32-bit targets never overflow (Stage 0) #20

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
907cbf1
Bound the psola/pvoc sample clocks so 32-bit targets never overflow
tap Sep 17, 2026
dd7a8ef
Address the two hostile reviews on the sample-clock fix
tap Sep 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
84 changes: 74 additions & 10 deletions include/tap/dsp/psola.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <type_traits>
#include <vector>

Expand All @@ -55,6 +56,15 @@ namespace tap::dsp {
/// spacing, the Hann windows sum to exactly one, and the output is the
/// input delayed by latency() (plus interpolation error) — pinned by
/// the test battery.
/// - Run time is unbounded: the sample clock is a fixed-width 32-bit
/// count bounded below 2 * clock_wrap() (a multiple of the ring size),
/// so no counter overflows on any target. The contract is that bound,
/// not an observable. The wrap moves the fractional mark positions in
/// magnitude only, so from 2 * clock_wrap() samples on (at most 2^19,
/// 10.9 s at 48 kHz; 519,552 samples for max_period 900) the output
/// diverges from an unbounded-clock reference at the ~1e-8 relative
/// level (measured 3.8e-9 absolute in double and one float ulp in
/// float, against a 0.35 peak). Pinned by the test battery.
template <typename Sample>
class basic_psola {
static_assert(std::is_same_v<Sample, float> || std::is_same_v<Sample, double>,
Expand All @@ -65,17 +75,23 @@ namespace tap::dsp {
static constexpr Sample k_min_ratio = Sample(0.25);
static constexpr Sample k_max_ratio = Sample(4);

/// @pre max_period >= 16 — the deepest period process() will be given.
/// @pre 16 <= max_period < 2^26 — the deepest period process() will be
/// given; the upper bound keeps 2 * clock_wrap() and every grain position
/// inside the int32 sample clock.
explicit basic_psola(size_t max_period)
: m_max_period(static_cast<Sample>(max_period))
, m_latency(2 * max_period + 2) {
assert(max_period >= 16);
assert(max_period < (size_t{1} << 26)); // see @pre
// One ring size serves both buffers; the clock wrap relies on that.
// Input history: a grain reaches back to (mark - period) and marks lag the
// input cursor by up to two periods -> three periods of history plus slack.
m_input.assign(4 * max_period + 8, Sample(0));
// Output accumulator: emission lags by latency(); grains extend up to one
// period past their mark -> latency + period ahead of the emit cursor.
m_accum.assign(4 * max_period + 8, Sample(0));
const size_t ring = 4 * max_period + 8;
m_input.assign(ring, Sample(0));
m_accum.assign(ring, Sample(0));
m_wrap = static_cast<std::int32_t>(ring * std::max<size_t>(1, static_cast<size_t>(k_clock_span) / ring));
clear();
}

Expand All @@ -84,6 +100,11 @@ namespace tap::dsp {

size_t max_period() const noexcept { return static_cast<size_t>(m_max_period); }

/// Period of the sample clock's wrap, in samples: the largest multiple of
/// the ring size not above 2^18 (or one ring, if the ring is larger). The
/// clock first wraps at 2 * clock_wrap() and every clock_wrap() after.
size_t clock_wrap() const noexcept { return static_cast<size_t>(m_wrap); }

/// Zero all running state (buffers, marks, counters).
void clear() noexcept {
std::fill(m_input.begin(), m_input.end(), Sample(0));
Expand All @@ -102,7 +123,7 @@ namespace tap::dsp {
const double r = std::clamp(static_cast<double>(ratio), static_cast<double>(k_min_ratio),
static_cast<double>(k_max_ratio));

m_input[static_cast<size_t>(m_n % static_cast<long>(m_input.size()))] = in;
m_input[static_cast<size_t>(m_n % static_cast<std::int32_t>(m_input.size()))] = in;

// Analysis marks: free-running, one per source period.
const double now = static_cast<double>(m_n);
Expand Down Expand Up @@ -134,17 +155,50 @@ namespace tap::dsp {

// Emit, then release the slot for reuse.
Sample y = Sample(0);
if (m_n >= static_cast<long>(m_latency)) {
const size_t slot =
static_cast<size_t>((m_n - static_cast<long>(m_latency)) % static_cast<long>(m_accum.size()));
y = m_accum[slot];
m_accum[slot] = Sample(0);
if (m_n >= static_cast<std::int32_t>(m_latency)) {
const size_t slot = static_cast<size_t>((m_n - static_cast<std::int32_t>(m_latency))
% static_cast<std::int32_t>(m_accum.size()));
y = m_accum[slot];
m_accum[slot] = Sample(0);
}
++m_n;
if (m_n >= 2 * m_wrap) {
shift_clock(-m_wrap);
}
return y;
}

/// Testing seam: advance the sample clock by `samples` (rounded down to a
/// multiple of the ring size) as if that many samples had elapsed with the
/// ring contents unchanged, wrapping the clock exactly as process() does.
/// Lets a test cross the wrap, or an elapsed count past 2^31, in O(1)
/// instead of processing that many samples. Fractional mark positions that
/// are not multiples of the new magnitude's ulp may round by that ulp.
/// Not a contract point; may change or disappear without a version note.
void advance_clock_for_testing(std::uint64_t samples) noexcept {
const std::uint64_t ring = m_input.size();
const std::uint64_t wrap = static_cast<std::uint64_t>(m_wrap);
const std::uint64_t target = static_cast<std::uint64_t>(m_n) + (samples / ring) * ring;
const std::uint64_t folded = (target < wrap) ? target : wrap + (target - wrap) % wrap;
shift_clock(static_cast<std::int32_t>(folded) - m_n);
}

private:
/// Move the clock and every absolute position by `by` samples, a multiple of
/// the ring size, so every ring index is unchanged. process() calls it with
/// -m_wrap when the clock reaches 2 * m_wrap; m_wrap is at least one ring, so
/// the clock stays at or above latency() and the warm-up guard stays true.
/// That subtraction is exact in double: m_wrap is an integer, hence a
/// multiple of every position's ulp, and the result's magnitude does not
/// exceed the position's.
void shift_clock(std::int32_t by) noexcept {
const double d = static_cast<double>(by);
m_n += by;
m_next_mark += d;
m_prev_mark += d;
m_next_synth += d;
}

/// Overlap-add one Hann grain: output slots o in [s - t, s + t] receive the
/// source at m + (o - s), read with Hermite interpolation (s is fractional).
void place_grain(double s, double m, double t, Sample gain) noexcept {
Expand Down Expand Up @@ -181,12 +235,22 @@ namespace tap::dsp {

static constexpr double k_pi = 3.14159265358979323846;

/// Span of the sample clock: clock_wrap() is the largest multiple of the
/// ring size not above this (or one ring, if the ring is larger). The
/// value is arbitrary within (longest test run, 2^29): 2^18 keeps every
/// position below 2^19 + 5 * max_period, i.e. with at least 33 fractional
/// bits, and wraps every 5.4 s at 48 kHz (first at 10.9 s) so the wrap
/// path is exercised routinely rather than once a shift. Not a contract
/// point; the tests pin clock_wrap() <= 2^18 as a literal.
static constexpr std::int32_t k_clock_span = std::int32_t{1} << 18;

Sample m_max_period;
size_t m_latency;

std::vector<Sample> m_input;
std::vector<Sample> m_accum;
long m_n{0};
std::int32_t m_n{0}; // sample clock, in [0, 2 * m_wrap); see shift_clock()
std::int32_t m_wrap{0}; // clock wrap period: a multiple of the ring size
double m_next_mark{0.0};
double m_prev_mark{0.0};
bool m_have_mark{false};
Expand Down
55 changes: 44 additions & 11 deletions include/tap/dsp/pvoc.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <type_traits>
#include <vector>

Expand All @@ -48,6 +49,12 @@ namespace tap::dsp {
/// Geometry (FFT size, 4x overlap) is fixed at construction and every
/// buffer is allocated there; process() is noexcept and allocation-free,
/// safe on a real-time audio thread. Latency is exactly latency() samples.
/// Run time is unbounded: the sample clock is a fixed-width 32-bit count
/// bounded below 2 * clock_wrap() (the overlap-add ring, 3 * fft_size, a
/// multiple of the input ring and the hop), so no counter overflows on any
/// target; the contract is that bound, not an observable. The wrap changes
/// no index and no value, so the output is bit-identical across it (pinned
/// by the test battery).
/// At ratio == 1 every region's bin offset and residual are zero, so the
/// output reconstructs the input's waveform delayed by exactly one FFT
/// frame (pinned by the test battery).
Expand Down Expand Up @@ -75,13 +82,15 @@ namespace tap::dsp {
static constexpr double k_env_floor = 1e-6; // |A| guard when inverting the envelope
static constexpr double k_max_boost = 16.0; // per-bin formant-correction gain cap

/// @pre fft_size is a power of two, >= 64.
/// @pre fft_size is a power of two in [64, 2^28]; the upper bound keeps
/// 6 * fft_size inside the int32 sample clock.
explicit basic_pvoc(size_t fft_size = 1024)
: m_n_size(static_cast<int>(fft_size))
, m_hop(static_cast<int>(fft_size) / k_overlap)
, m_bins(static_cast<int>(fft_size) / 2 + 1)
, m_fft(fft_size) {
assert(fft_size >= 64 && (fft_size & (fft_size - 1)) == 0);
assert(fft_size <= (size_t{1} << 28)); // see @pre

m_window.assign(static_cast<size_t>(m_n_size), Sample(0));
for (int i = 0; i < m_n_size; ++i) {
Expand Down Expand Up @@ -130,6 +139,10 @@ namespace tap::dsp {
void set_formant(bool on) noexcept { m_formant = on; }
bool formant() const noexcept { return m_formant; }

/// Period of the sample clock's wrap, in samples: the overlap-add ring.
/// The clock first wraps at 2 * clock_wrap() and every clock_wrap() after.
size_t clock_wrap() const noexcept { return m_accum.size(); }

/// Zero all running state (buffers, phases, counters).
void clear() noexcept {
std::fill(m_input.begin(), m_input.end(), Sample(0));
Expand All @@ -141,31 +154,51 @@ namespace tap::dsp {

/// Consume one input sample; produce the output sample for time n - latency().
Sample process(Sample in, Sample ratio) noexcept {
const long in_size = static_cast<long>(m_input.size());
const long an = static_cast<long>(m_accum.size());
const std::int32_t in_size = static_cast<std::int32_t>(m_input.size());
const std::int32_t an = static_cast<std::int32_t>(m_accum.size());

m_input[static_cast<size_t>(m_n % in_size)] = in;

if ((m_n + 1) % m_hop == 0 && m_n + 1 >= static_cast<long>(m_n_size)) {
if ((m_n + 1) % m_hop == 0 && m_n + 1 >= m_n_size) {
const double r = std::clamp(static_cast<double>(ratio), static_cast<double>(k_min_ratio),
static_cast<double>(k_max_ratio));
run_frame(r);
}

Sample y = Sample(0);
if (m_n >= static_cast<long>(m_n_size)) {
const size_t slot = static_cast<size_t>((m_n - static_cast<long>(m_n_size)) % an);
if (m_n >= m_n_size) {
const size_t slot = static_cast<size_t>((m_n - m_n_size) % an);
y = m_accum[slot];
m_accum[slot] = Sample(0);
}
++m_n;
// Sample clock: runs in [0, 2 * an). Pulling it back by an — a multiple
// of the input ring and of the hop — leaves every ring index and the
// frame schedule unchanged, keeps it at or above the warm-up guards
// (an >= m_n_size), and bounds it forever without a 64-bit division.
if (m_n >= 2 * an) {
m_n -= an;
}
return y;
}

/// Testing seam: advance the sample clock by `samples` (rounded down to a
/// multiple of clock_wrap()) as if that many samples had elapsed with the
/// ring contents unchanged, wrapping the clock exactly as process() does.
/// Lets a test cross an elapsed count past 2^31 in O(1) instead of
/// processing that many samples. Not a contract point; may change or
/// disappear without a version note.
void advance_clock_for_testing(std::uint64_t samples) noexcept {
const std::uint64_t wrap = m_accum.size();
const std::uint64_t target = static_cast<std::uint64_t>(m_n) + (samples / wrap) * wrap;
const std::uint64_t folded = (target < wrap) ? target : wrap + (target - wrap) % wrap;
m_n = static_cast<std::int32_t>(folded);
}

private:
void run_frame(double r) noexcept {
const long in_size = static_cast<long>(m_input.size());
const long start = m_n + 1 - static_cast<long>(m_n_size);
const std::int32_t in_size = static_cast<std::int32_t>(m_input.size());
const std::int32_t start = m_n + 1 - m_n_size;

// analysis: window the newest N samples and transform
for (int i = 0; i < m_n_size; ++i) {
Expand Down Expand Up @@ -271,8 +304,8 @@ namespace tap::dsp {

// synthesis window + COLA-normalized overlap-add (fold in the raw
// inverse's 2/N normalization here, one multiply per sample)
const long an = static_cast<long>(m_accum.size());
const Sample norm = m_cola_norm * (Sample(2) / static_cast<Sample>(m_n_size));
const std::int32_t an = static_cast<std::int32_t>(m_accum.size());
const Sample norm = m_cola_norm * (Sample(2) / static_cast<Sample>(m_n_size));
for (int i = 0; i < m_n_size; ++i) {
const size_t slot = static_cast<size_t>((start + i) % an);
m_accum[slot] += m_synth[static_cast<size_t>(i)] * m_window[static_cast<size_t>(i)] * norm;
Expand Down Expand Up @@ -362,7 +395,7 @@ namespace tap::dsp {
std::vector<double> m_lpc_tmp;
std::vector<double> m_env;
bool m_formant{false};
long m_n{0};
std::int32_t m_n{0}; // sample clock, in [0, 2 * m_accum.size()); see process()
};

/// Double-precision shifter — the desktop/golden-model profile.
Expand Down
40 changes: 40 additions & 0 deletions tests/test_psola.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
// oracle is tap::dsp::yin, certified by its own battery.

#include <cmath>
#include <cstdint>
#include <vector>

#include <gtest/gtest.h>
Expand Down Expand Up @@ -146,6 +147,45 @@ namespace {
}
}

TYPED_TEST(psola_test, ClockWrapIsSeamlessAndCountersDoNotOverflow) {
// The run-time contract is a bound, not an observable: the 32-bit sample
// clock stays below 2 * clock_wrap(), a multiple of the ring size, and
// the wrap changes nothing but the magnitude of the fractional mark
// positions. Reference: a shifter run from zero. Subject: the same
// shifter with its clock advanced past 2^31 elapsed samples through the
// documented O(1) seam (the count a 48 kHz Cortex-M or Windows build
// reaches after 12.4 h), then run for more than one wrap period so the
// wrap fires inside the run. Ratio 1.5 makes the synthesis step t / r
// inexact, so the wrap's magnitude change is exercised; the outputs agree
// to the rounding of those positions (measured 2e-9 against a 0.35 peak),
// five orders below any dropped or doubled grain.
using shifter_t = tap::dsp::basic_psola<TypeParam>;
shifter_t ref(900);
shifter_t sub(900);
ASSERT_LE(sub.clock_wrap(), size_t{1} << 18); // the documented span
ASSERT_EQ(sub.clock_wrap() % (4 * 900 + 8), 0u); // a multiple of the ring size

const TypeParam period = static_cast<TypeParam>(k_sr / 150.0);
const int warm = 12000; // seam applied past the warm-up, mid-stream
const int run = (1 << 18) + 8192; // > clock_wrap(), so the wrap fires inside the run
double worst = 0.0;
for (int t = 0; t < warm + run; ++t) {
if (t == warm) {
sub.advance_clock_for_testing((std::uint64_t{1} << 31) + 12345);
}
double x = 0.0;
for (int h = 1; h <= 20; ++h) {
x += std::sin(2.0 * k_pi * 150.0 * h * t / k_sr) / h;
}
const TypeParam xs = static_cast<TypeParam>(x / 3.6);
const double a = static_cast<double>(ref.process(xs, period, TypeParam(1.5)));
const double b = static_cast<double>(sub.process(xs, period, TypeParam(1.5)));
ASSERT_TRUE(std::isfinite(b)) << "t " << t;
worst = std::max(worst, std::abs(a - b));
}
EXPECT_LT(worst, 1e-6);
}

TEST(psola_cross_precision, FloatAgreesWithDoubleGoldenModel) {
tap::dsp::psola gold(900);
tap::dsp::psola32 fast(900);
Expand Down
Loading
Loading