diff --git a/include/tap/dsp/psola.h b/include/tap/dsp/psola.h index c208262..bd22d0f 100644 --- a/include/tap/dsp/psola.h +++ b/include/tap/dsp/psola.h @@ -31,6 +31,7 @@ #include #include #include +#include #include #include @@ -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 class basic_psola { static_assert(std::is_same_v || std::is_same_v, @@ -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(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(ring * std::max(1, static_cast(k_clock_span) / ring)); clear(); } @@ -84,6 +100,11 @@ namespace tap::dsp { size_t max_period() const noexcept { return static_cast(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(m_wrap); } + /// Zero all running state (buffers, marks, counters). void clear() noexcept { std::fill(m_input.begin(), m_input.end(), Sample(0)); @@ -102,7 +123,7 @@ namespace tap::dsp { const double r = std::clamp(static_cast(ratio), static_cast(k_min_ratio), static_cast(k_max_ratio)); - m_input[static_cast(m_n % static_cast(m_input.size()))] = in; + m_input[static_cast(m_n % static_cast(m_input.size()))] = in; // Analysis marks: free-running, one per source period. const double now = static_cast(m_n); @@ -134,17 +155,50 @@ namespace tap::dsp { // Emit, then release the slot for reuse. Sample y = Sample(0); - if (m_n >= static_cast(m_latency)) { - const size_t slot = - static_cast((m_n - static_cast(m_latency)) % static_cast(m_accum.size())); - y = m_accum[slot]; - m_accum[slot] = Sample(0); + if (m_n >= static_cast(m_latency)) { + const size_t slot = static_cast((m_n - static_cast(m_latency)) + % static_cast(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(m_wrap); + const std::uint64_t target = static_cast(m_n) + (samples / ring) * ring; + const std::uint64_t folded = (target < wrap) ? target : wrap + (target - wrap) % wrap; + shift_clock(static_cast(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(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 { @@ -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 m_input; std::vector 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}; diff --git a/include/tap/dsp/pvoc.h b/include/tap/dsp/pvoc.h index 2cc9307..13861de 100644 --- a/include/tap/dsp/pvoc.h +++ b/include/tap/dsp/pvoc.h @@ -35,6 +35,7 @@ #include #include #include +#include #include #include @@ -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). @@ -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(fft_size)) , m_hop(static_cast(fft_size) / k_overlap) , m_bins(static_cast(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(m_n_size), Sample(0)); for (int i = 0; i < m_n_size; ++i) { @@ -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)); @@ -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(m_input.size()); - const long an = static_cast(m_accum.size()); + const std::int32_t in_size = static_cast(m_input.size()); + const std::int32_t an = static_cast(m_accum.size()); m_input[static_cast(m_n % in_size)] = in; - if ((m_n + 1) % m_hop == 0 && m_n + 1 >= static_cast(m_n_size)) { + if ((m_n + 1) % m_hop == 0 && m_n + 1 >= m_n_size) { const double r = std::clamp(static_cast(ratio), static_cast(k_min_ratio), static_cast(k_max_ratio)); run_frame(r); } Sample y = Sample(0); - if (m_n >= static_cast(m_n_size)) { - const size_t slot = static_cast((m_n - static_cast(m_n_size)) % an); + if (m_n >= m_n_size) { + const size_t slot = static_cast((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(m_n) + (samples / wrap) * wrap; + const std::uint64_t folded = (target < wrap) ? target : wrap + (target - wrap) % wrap; + m_n = static_cast(folded); + } + private: void run_frame(double r) noexcept { - const long in_size = static_cast(m_input.size()); - const long start = m_n + 1 - static_cast(m_n_size); + const std::int32_t in_size = static_cast(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) { @@ -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(m_accum.size()); - const Sample norm = m_cola_norm * (Sample(2) / static_cast(m_n_size)); + const std::int32_t an = static_cast(m_accum.size()); + const Sample norm = m_cola_norm * (Sample(2) / static_cast(m_n_size)); for (int i = 0; i < m_n_size; ++i) { const size_t slot = static_cast((start + i) % an); m_accum[slot] += m_synth[static_cast(i)] * m_window[static_cast(i)] * norm; @@ -362,7 +395,7 @@ namespace tap::dsp { std::vector m_lpc_tmp; std::vector 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. diff --git a/tests/test_psola.cpp b/tests/test_psola.cpp index 33dbf6c..33bc0cb 100644 --- a/tests/test_psola.cpp +++ b/tests/test_psola.cpp @@ -8,6 +8,7 @@ // oracle is tap::dsp::yin, certified by its own battery. #include +#include #include #include @@ -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; + 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(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(x / 3.6); + const double a = static_cast(ref.process(xs, period, TypeParam(1.5))); + const double b = static_cast(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); diff --git a/tests/test_pvoc.cpp b/tests/test_pvoc.cpp index 9f37483..3f42730 100644 --- a/tests/test_pvoc.cpp +++ b/tests/test_pvoc.cpp @@ -8,6 +8,7 @@ // oracle is tap::dsp::yin, certified by its own battery. #include +#include #include #include @@ -199,6 +200,50 @@ namespace { } } + TYPED_TEST(pvoc_test, ClockWrapIsBitExactAndCountersDoNotOverflow) { + // The run-time contract is a bound, not an observable: the 32-bit sample + // clock stays below 2 * clock_wrap(), the overlap-add ring (3 * fft_size, + // a multiple of the input ring and the hop), so the wrap changes nothing + // at all. Every 1 s run above already wraps 14 times; this test makes + // the pin explicit. Reference: a shifter run from zero, whose identity + // reconstruction must hold over the WHOLE run (every wrap included), not + // just the tail. 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); its output is + // bit-identical to the reference's. The seam is applied at a warm-up + // BELOW clock_wrap(): the reference clock then reads 2048 while the + // seeded one folds to 3072 + 2048 = 5120, and the two differ until both + // wrap to 3072 1024 samples later, which is what ASSERT_EQ pins (a clock + // value and that value plus clock_wrap() are indistinguishable). Past + // clock_wrap() the seam would fold to the reference's own value and the + // comparison would be vacuous. + tap::dsp::basic_pvoc ref(1024); + tap::dsp::basic_pvoc sub(1024); + ref.set_formant(true); // the LPC path runs in the frame too + sub.set_formant(true); + ASSERT_EQ(sub.clock_wrap(), 3u * 1024u); + + const double tolerance = std::is_same_v ? 1e-8 : 2e-3; + const int latency = static_cast(ref.latency()); + const int warm = 2048; // past the warm-up guard (fft_size), below clock_wrap() + const int run = 48000; + 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) + 777); + } + const TypeParam x = static_cast(std::sin(2.0 * k_pi * 440.0 * t / k_sr)); + const TypeParam a = ref.process(x, TypeParam(1)); + const TypeParam b = sub.process(x, TypeParam(1)); + ASSERT_EQ(a, b) << "t " << t; + if (t >= 2 * latency) { // steady state: four frames overlap from here on + const double expected = std::sin(2.0 * k_pi * 440.0 * (t - latency) / k_sr); + worst = std::max(worst, std::abs(static_cast(a) - expected)); + } + } + EXPECT_LT(worst, tolerance); + } + TEST(pvoc_cross_precision, FloatAgreesWithDoubleGoldenModel) { tap::dsp::pvoc gold(1024); tap::dsp::pvoc32 fast(1024);