From 76e2b6ab47b8b15b1ef461acb55a344f7ab4998e Mon Sep 17 00:00:00 2001 From: Timothy Place Date: Thu, 17 Sep 2026 18:18:08 +0000 Subject: [PATCH 1/2] Stage 3a: sample_traits, named Q ladder, fft_arith trait, decimate re-pin The sample-format substrate now covers the golden model. sample_traits (coeff = double, accum = double, k_coeff_scale = 1.0, identity conversions) lands, so the sample_type concept covers four types and a traits-based primitive can instantiate its double profile; test_sample_traits.cpp's static_assert(!sample_type) flips. The concept now requires coeff, accum, k_coeff_scale, a new k_is_fixed_point, and (documented, static_asserted per type) that accum{} is the additive identity. Every trait member is constexpr. The Q ladder is spelled as named constants (k_sample_frac_bits 15/31, k_coeff_frac_bits 14/30, k_accum_pre_shift 0/16) with the accumulator format and the 14-bit finalize shift derived from them and static_asserted; k_coeff_scale is 2^k_coeff_frac_bits. Bit-identical: the Q15/Q31 decimate outputs and quantized rows were dumped before and after and match byte for byte, and every existing test passes unchanged. fft/fft_arith.h is the butterfly arithmetic trait the Stage 3b fixed-point kernel is written against: fft_arith carries mul_coeff (int32 x Q1.30 -> int64, >> 30 with one round-half-up, saturating), saturating add/sub, shr_round and headroom_bits; fft_arith is the Q15 I/O width (widen << 14 with two guard bits, saturating round-half-up narrow) over the int32 trait as its work type; float/double are the same names over plain arithmetic. Twiddles are Q1.30 for both fixed profiles and 1.0 is representable. tests/test_fft_arith.cpp pins every contract number. quantize.h selects its algorithm by k_is_fixed_point and routes the +/-1 correction through a saturating step so a tap at the rail cannot wrap. decimate: basic_decimator instantiates; the float/double contract tests are typed; the Q15 cross-precision pin is re-pointed at double (measured 1.57e-4, pinned 3.2e-4) with a Q31-vs-double sibling (2.45e-9, pinned 4.9e-9) and a double-vs-numpy pin at the reference's own float32 rounding (5.75e-8, pinned 1.15e-7); the float-vs-numpy pin is unchanged. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019ZPTzNxo5Fe4EtpXXKf7Sy --- README.md | 37 ++++- include/tap/dsp/decimate.h | 16 +- include/tap/dsp/fft/fft_arith.h | 251 ++++++++++++++++++++++++++++++++ include/tap/dsp/quantize.h | 27 +++- include/tap/dsp/sample_traits.h | 171 ++++++++++++++++++---- tests/CMakeLists.txt | 1 + tests/test_decimate.cpp | 154 ++++++++++++++------ tests/test_fft_arith.cpp | 249 +++++++++++++++++++++++++++++++ tests/test_fir_kernels.cpp | 9 +- tests/test_quantize.cpp | 32 ++++ tests/test_sample_traits.cpp | 65 ++++++++- 11 files changed, 912 insertions(+), 100 deletions(-) create mode 100644 include/tap/dsp/fft/fft_arith.h create mode 100644 tests/test_fft_arith.cpp diff --git a/README.md b/README.md index b438235..44aa32a 100644 --- a/README.md +++ b/README.md @@ -245,17 +245,31 @@ constexpr (the header's design note does the arithmetic); run it in a constructor, off the audio path. Also exports `solve_dense`, the small dense solver the compensated design and the analysis instruments share. -### `tap/dsp/sample_traits.h` — sample formats: float, Q15, Q31 +### `tap/dsp/sample_traits.h` — sample formats: double, float, Q15, Q31 The family's sample-format substrate: how each sample type stores coefficients, accumulates dot products, and rounds/saturates back to samples. +`double` is a sample format because it is the golden model of every +primitive: a traits-based primitive instantiates its reference profile +through the same substrate as its embedded profiles, and the cross-precision +pins measure float/Q15/Q31 against it. | Type | Coefficients | Accumulation | Output | |---|---|---|---| +| `double` | double | double | identity (the golden model) | | `float` | float | double | plain cast | | `std::int16_t` | Q1.14 | int64, exact | single Q29→Q15 round-half-up, saturating | | `std::int32_t` | Q1.30 | int64, products pre-shifted to Q45 | Q45→Q31, saturating | +The Q ladder is spelled as named constants on each fixed-point +specialization (`k_sample_frac_bits`, `k_coeff_frac_bits`, +`k_accum_pre_shift`) with the accumulator format and the single 14-bit +finalize shift derived from them and `static_assert`ed; `k_coeff_scale` is +`2^k_coeff_frac_bits`, and `k_is_fixed_point` is what selects the +fixed-point algorithm in `quantize.h`. Every member is `constexpr`, and the +`sample_type` concept requires a value-initialized accumulator to be the +additive identity. + **Fixed point is a first-class embedded direction, not a legacy path.** The Q15/Q31 profiles exist for targets where double (sometimes any float) is unaffordable — SampleRateTap measured its float datapath at ~19× the @@ -273,6 +287,21 @@ This header is the format core only. Engine-specific extensions (e.g. SampleRateTap's inter-phase coefficient blending) derive from these specializations and refine the `tap::dsp::sample_type` concept. +### `tap/dsp/fft/fft_arith.h` — butterfly arithmetic for the FFT profiles + +The sibling trait the fixed-point real FFT is written against (its docstrings +are that kernel's specification; `tests/test_fft_arith.cpp` pins every +number). One int32 kernel serves both fixed profiles: `fft_arith` +carries `mul_coeff` (int32 × Q1.30 → int64, `>> 30` with one round-half-up, +saturating), saturating `add` / `sub`, `shr_round` (round-half-up), and +`headroom_bits` (the block's shared redundant sign bits, via +`std::countl_zero`); `fft_arith` is the Q15 I/O width — +`widen` (`<< 14`, two guard bits) and `narrow` (round-half-up, saturating) — +and names the int32 trait as its `work`. Twiddles are +`sample_traits::coeff` (Q1.30) for both fixed profiles, and 1.0 +is representable. The `float` / `double` specializations are the same names +over plain arithmetic. Everything is `constexpr` and `noexcept`. + ### `tap/dsp/fir_kernels.h` — dot-product kernels The FIR hot loops, target-gated the way SampleRateTap's optimization campaign @@ -290,8 +319,10 @@ targets prefer which layout. `quantize_row_preserving_sum`: quantizes one polyphase branch to a fixed-point coefficient format while preserving the row's DC sum *exactly* (largest-remainder distribution of the rounding residual — "the coefficients -of every phase must add to one", R. Bristow-Johnson, music-dsp). Plain -conversion for float. Design-time code. +of every phase must add to one", R. Bristow-Johnson, music-dsp). Selected by +the trait's `k_is_fixed_point`: plain conversion for double and float. Each +correction step saturates, so a tap at the format's rail is never wrapped. +Design-time code. ### `tap/dsp/analysis/` — measurement instruments diff --git a/include/tap/dsp/decimate.h b/include/tap/dsp/decimate.h index d6c971e..5bb2025 100644 --- a/include/tap/dsp/decimate.h +++ b/include/tap/dsp/decimate.h @@ -7,9 +7,9 @@ // 96 kHz in, 16 kHz out. Built in RatioTap's pattern — the ratio is a // compile-time type, the prototype is a Kaiser-windowed sinc from kaiser.h, // the hot loop is fir_kernels.h's dot_row over the sample_traits.h formats -// (float golden, Q15/Q31 fixed point) — but deliberately NOT RatioTap: that -// library's charter is 44.1 <-> 48 only. A 44.1 kHz host composes RatioTap's -// 44.1 -> 48 in front of the by-3 stage here. +// (double golden, float embedded, Q15/Q31 fixed point) — but deliberately +// NOT RatioTap: that library's charter is 44.1 <-> 48 only. A 44.1 kHz host +// composes RatioTap's 44.1 -> 48 in front of the by-3 stage here. // // Contract, as numbers: // - Ratios: 2, 3 and 6 (input 32 / 48 / 96 kHz for a 16 kHz output). The @@ -37,10 +37,14 @@ // bands stop at 7.6 kHz and its features carry little above 7 kHz, and // economy costs 121 MACs per 16 kHz output (by 3) — under 20k MACs per // 10 ms hop. transparent exists for offline use. -// - Sample formats: float (double accumulation, the golden model, pinned +// - Sample formats: double (the golden model, as for every primitive), +// float (double accumulation, the embedded profile, pinned // sample-for-sample against a committed numpy reference), Q15 and Q31 // through sample_traits.h with row-sum-preserving quantization so DC -// gain stays exactly 1. Mono: the consumer is single-channel by charter. +// gain stays exactly 1; the Q15 profile is pinned against double as a +// measured number. The decimate_by_* aliases stay float: they name the +// 16 kHz front end's deployed profile. Mono: the consumer is +// single-channel by charter. // // Construction designs the filter (runtime double, off the audio path) and // allocates; process() and reset() are noexcept and allocation-free. @@ -185,7 +189,7 @@ namespace tap::dsp { std::size_t m_phase = 0; ///< inputs since the last emitted output, in [0, M) }; - using decimate_by_2 = basic_decimator; ///< 32 kHz -> 16 kHz, float golden + using decimate_by_2 = basic_decimator; ///< 32 kHz -> 16 kHz, float embedded profile using decimate_by_3 = basic_decimator; ///< 48 kHz -> 16 kHz using decimate_by_6 = basic_decimator; ///< 96 kHz -> 16 kHz diff --git a/include/tap/dsp/fft/fft_arith.h b/include/tap/dsp/fft/fft_arith.h new file mode 100644 index 0000000..b91af46 --- /dev/null +++ b/include/tap/dsp/fft/fft_arith.h @@ -0,0 +1,251 @@ +/// @file fft_arith.h +/// @brief Butterfly arithmetic trait for the real FFT profiles: double, float, Q15 I/O, Q31. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// The sibling of sample_traits.h for the FFT. The FIR trait's mac()/finalize() +// shape (a long exact accumulation with one rounding at the end) is the wrong +// one for a butterfly, which multiplies a sample by a unit-magnitude twiddle +// and adds two samples, stage after stage, with the data staying in its own +// format throughout. This header is the arithmetic the fixed-point kernel +// (fft/fixed_point.h, Stage 3b of docs/audit-fft-and-code-smells.md) is +// written against; the docstrings below are that kernel's specification, and +// tests/test_fft_arith.cpp pins every number in them. +// +// Design (audit Part 7, from Welch 1969 and Oppenheim & Weinstein 1972 on +// finite-register FFT error, and the substrate's rule that a Q ladder is +// visible at the use site): +// +// - One int32 kernel for both fixed profiles. Q31 data is int32 as is. The +// Q15 profile widens int16 to int32 on input (<< 14, two guard bits) and +// narrows on output (round-half-up, saturating, the substrate's finalize +// rule); fft_arith is exactly that pair of I/O-width +// conversions plus the name of the int32 trait that does the work. +// - Twiddles are sample_traits::coeff, i.e. Q1.30 int32, for +// BOTH fixed profiles (fft_arith::coeff is the same type). +// 1.0 is representable (k_coeff_one = 2^30 < 2^31 - 1), so DC and Nyquist +// twiddles need no special case; every twiddle is generated in double and +// rounded once by make_coeff (round_sat: half away from zero, saturating). +// - The multiply is int32 x Q1.30 -> int64, then >> 30 with ONE round-half-up. +// That is the single documented rounding point per product; Armv7E-M / +// Armv8-M's SMMULR (32x32 high half with rounding) is the designated seam +// behind the same contract, as SMLALD is for the FIR kernels. +// - Growth is handled by the kernel with shr_round (fixed scaling) or +// headroom_bits + shr_round (block floating point); the trait only supplies +// the operations, never a policy. +// +// Every operation is constexpr and noexcept. The int32 operations preserve +// the data's Q format whatever it is (Q0.31 for the Q31 profile, Q2.29 for +// the widened Q15 profile): only the twiddle's format (Q1.30) enters the +// arithmetic. The floating specializations are the same names over plain +// float/double arithmetic so a profile-generic wrapper can be written once. +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "tap/dsp/sample_traits.h" + +namespace tap::dsp { + + /// Primary template intentionally undefined; specialize per sample type. + template + struct fft_arith; + + /// Floating profiles (double is the golden model, float the embedded + /// profile): every operation is the plain IEEE one, nothing rounds + /// beyond the format, nothing saturates, and there are no guard bits. + /// wide == sample, widen/narrow are identities, headroom_bits is 0 + /// because a floating block never needs block scaling. + template + struct fft_arith { + using sample = F; + using wide = F; ///< the kernel's data width: the sample itself + using coeff = F; ///< twiddle type + using work = fft_arith; + + static constexpr bool k_is_fixed_point = false; + static constexpr int k_guard_bits = 0; + static constexpr int k_coeff_frac_bits = 0; + static constexpr F k_coeff_one = F{1}; + + static constexpr coeff make_coeff(double w) noexcept { return static_cast(w); } + + static constexpr wide widen(sample x) noexcept { return x; } + static constexpr sample narrow(wide x) noexcept { return x; } + + static constexpr F mul_coeff(F x, coeff w) noexcept { return x * w; } + static constexpr F add(F a, F b) noexcept { return a + b; } + static constexpr F sub(F a, F b) noexcept { return a - b; } + + /// x * 2^-bits, exact (a power-of-two scale never rounds in binary + /// floating point short of underflow). + static constexpr F shr_round(F x, int bits) noexcept { + F s = F{1}; + for (int i = 0; i < bits; ++i) { + s *= F{0.5}; + } + return x * s; + } + + /// Floating data carries its exponent per element: no block headroom. + static constexpr int headroom_bits(const F* /*x*/, std::size_t /*n*/) noexcept { return 0; } + }; + + // ANCHOR: fa_q31 + /// The int32 kernel's arithmetic: Q31 data, and the Q15 profile's widened + /// data. Contracts, per operation (format in -> format out, rounding, + /// saturation), each pinned by the named test in tests/test_fft_arith.cpp: + /// + /// - make_coeff(double) -> coeff : Q1.30, round half away from zero, + /// saturating (sample_traits::make_coeff). 1.0 -> 2^30 + /// exactly. `TwiddleUnityIsRepresentable`. + /// - mul_coeff(x, w) -> sample : Qm.n x Q1.30 -> Qm.n. Exact int64 + /// product, then >> 30 with round-half-up (add 2^29, arithmetic shift): + /// the single rounding point. Saturates on the way back to int32; with + /// |w| <= 1.0 the only input that saturates is INT32_MIN x (-1.0) + /// (+2^31 -> INT32_MAX). `MulCoeffIsExactOnRepresentableProducts`, + /// `MulCoeffRoundsHalfUp`, `MulCoeffSaturatesAtTheRail`. + /// - add(a, b), sub(a, b) -> sample : Qm.n +/- Qm.n -> Qm.n, computed in + /// int64 and saturated to [INT32_MIN, INT32_MAX]. Never wraps. + /// `AddSubSaturateAtTheRails`. + /// - shr_round(x, bits) -> sample : x * 2^-bits, round-half-up in the + /// discarded bits (add 2^(bits-1), arithmetic shift, in int64 so the + /// rounding add cannot overflow). bits == 0 is the identity; the + /// result always fits for 0 <= bits <= 31. `ShrRoundRoundsHalfUp`. + /// - headroom_bits(x, n) -> int : the number of redundant sign bits + /// shared by the whole block — the largest s such that every x[i] << s + /// still fits in int32 (Arm's CLS, computed as countl_zero of the OR of + /// the one's-complement magnitudes x ^ (x >> 31), minus the sign bit). + /// 31 for an all-zero or empty block, 30 for {1}, 0 for INT32_MAX or + /// INT32_MIN. `HeadroomBitsOnZeroOneAndFullScale`. + /// - widen / narrow : identities (wide == sample). + template <> + struct fft_arith { + using sample = std::int32_t; + using wide = std::int32_t; ///< the kernel's data width + using coeff = sample_traits::coeff; ///< Q1.30 twiddle + using work = fft_arith; + /// Product width: int32 x Q1.30 is exact in int64. + using product = std::int64_t; + + static constexpr bool k_is_fixed_point = true; + static constexpr int k_guard_bits = 0; ///< Q0.31 data has no spare width + /// Q1.30: the twiddle's fraction bits, shared with the FIR Q31 coefficient. + static constexpr int k_coeff_frac_bits = sample_traits::k_coeff_frac_bits; + static constexpr coeff k_coeff_one = coeff{1} << k_coeff_frac_bits; ///< 1.0 in Q1.30 + static_assert(k_coeff_frac_bits == 30 && k_coeff_one == 1073741824, "twiddles are Q1.30; 1.0 is representable"); + static constexpr int k_sample_bits = 32; + + static constexpr coeff make_coeff(double w) noexcept { return sample_traits::make_coeff(w); } + + static constexpr wide widen(sample x) noexcept { return x; } + static constexpr sample narrow(wide x) noexcept { return x; } + + static constexpr sample mul_coeff(sample x, coeff w) noexcept { + const product p = static_cast(x) * static_cast(w); + return detail::clamp_sat((p + (product{1} << (k_coeff_frac_bits - 1))) >> k_coeff_frac_bits); + } + + static constexpr sample add(sample a, sample b) noexcept { + return detail::clamp_sat(static_cast(a) + static_cast(b)); + } + + static constexpr sample sub(sample a, sample b) noexcept { + return detail::clamp_sat(static_cast(a) - static_cast(b)); + } + + static constexpr sample shr_round(sample x, int bits) noexcept { + if (bits <= 0) { + return x; + } + return static_cast((static_cast(x) + (std::int64_t{1} << (bits - 1))) >> bits); + } + + static constexpr int headroom_bits(const sample* x, std::size_t n) noexcept { + std::uint32_t magnitudes = 0; + for (std::size_t i = 0; i < n; ++i) { + // One's-complement magnitude: x for x >= 0, ~x for x < 0, so + // INT32_MIN maps to 0x7fffffff (zero headroom) and -1 to 0. + magnitudes |= static_cast(x[i] ^ (x[i] >> (k_sample_bits - 1))); + } + return std::countl_zero(magnitudes) - 1; + } + }; + // ANCHOR_END: fa_q31 + + // ANCHOR: fa_q15 + /// The Q15 profile's I/O width. The kernel does not compute in int16: it + /// widens each Q0.15 sample to int32 on the way in and narrows on the way + /// out, and every butterfly runs through fft_arith (named + /// here as `work`). Contracts: + /// + /// - widen(x) -> wide : Q0.15 -> Q2.29, x << 14. Exact. The two guard + /// bits above full scale are what let worst-case radix-4 growth run + /// under fixed scaling without an input-level contract (audit Part 7). + /// INT16_MIN -> -2^29. `WidenPlacesTwoGuardBits`. + /// - narrow(y) -> sample : Q2.29 -> Q0.15, (y + 2^13) >> 14 in int64, + /// round-half-up, saturating to [INT16_MIN, INT16_MAX] — the + /// substrate's finalize rule. `NarrowRoundsHalfUpAndSaturates`. + /// - narrow(widen(x)) == x for every int16 value. `WidenNarrowRoundTripsEveryInt16`. + /// - coeff, make_coeff, k_coeff_one : the int32 trait's Q1.30 twiddles. + /// The Q15 profile does NOT use Q1.14 twiddles. + template <> + struct fft_arith { + using sample = std::int16_t; + using wide = std::int32_t; ///< the kernel's data width + using coeff = sample_traits::coeff; ///< Q1.30 twiddle, same as Q31 + using work = fft_arith; ///< the arithmetic the kernel runs on `wide` + + static constexpr bool k_is_fixed_point = true; + static constexpr int k_guard_bits = 2; ///< spare bits above full scale after widen() + static constexpr int k_widen_shift = 32 - 16 - k_guard_bits; ///< 14: Q0.15 -> Q2.29 + static constexpr int k_coeff_frac_bits = work::k_coeff_frac_bits; + static constexpr coeff k_coeff_one = work::k_coeff_one; + static_assert(k_widen_shift == 14, "Q15 I/O: two guard bits means a 14-bit widen"); + + static constexpr coeff make_coeff(double w) noexcept { return work::make_coeff(w); } + + static constexpr wide widen(sample x) noexcept { return static_cast(x) << k_widen_shift; } + + static constexpr sample narrow(wide x) noexcept { + return detail::clamp_sat((static_cast(x) + (std::int64_t{1} << (k_widen_shift - 1))) + >> k_widen_shift); + } + }; + // ANCHOR_END: fa_q15 + + /// Satisfied by the four FFT profiles' sample types: the I/O-width pair, + /// the twiddle type and its generator, and a `work` trait that carries the + /// butterfly operations over `wide`. + template + concept fft_sample_type = requires(T x, typename fft_arith::wide y, typename fft_arith::coeff w, double d) { + typename fft_arith::wide; + typename fft_arith::coeff; + typename fft_arith::work; + requires std::is_same_v::k_is_fixed_point), const bool>; + { fft_arith::make_coeff(d) } -> std::same_as::coeff>; + { fft_arith::widen(x) } -> std::same_as::wide>; + { fft_arith::narrow(y) } -> std::same_as; + { fft_arith::work::mul_coeff(y, w) } -> std::same_as::wide>; + { fft_arith::work::add(y, y) } -> std::same_as::wide>; + { fft_arith::work::sub(y, y) } -> std::same_as::wide>; + { fft_arith::work::shr_round(y, 1) } -> std::same_as::wide>; + { fft_arith::work::headroom_bits(&y, std::size_t{1}) } -> std::same_as; + }; + + static_assert(fft_sample_type); + static_assert(fft_sample_type); + static_assert(fft_sample_type); + static_assert(fft_sample_type); + + // Both fixed profiles share one twiddle type and 1.0 is representable in it. + static_assert(std::is_same_v::coeff, fft_arith::coeff>); + static_assert(fft_arith::make_coeff(1.0) == fft_arith::k_coeff_one); + static_assert(fft_arith::k_coeff_one < std::numeric_limits::max()); + +} // namespace tap::dsp diff --git a/include/tap/dsp/quantize.h b/include/tap/dsp/quantize.h index ff46f28..e7d9a36 100644 --- a/include/tap/dsp/quantize.h +++ b/include/tap/dsp/quantize.h @@ -12,7 +12,6 @@ #include #include #include -#include #include #include "tap/dsp/sample_traits.h" @@ -35,18 +34,25 @@ namespace tap::dsp { /// llround(exact_sum * k_coeff_scale), so a table built row-by-row holds /// DC gain within one coefficient LSB across all phases. /// - /// For floating-point coefficient types this is a plain conversion (no - /// correction runs; k_coeff_scale is 1). Design-time code: allocates a - /// scratch vector for integer formats, so keep it off the audio path. + /// The algorithm is selected by the trait's k_is_fixed_point property: + /// for the floating formats (double, float) this is a plain conversion + /// (no correction runs; k_coeff_scale is 1). Each +/-1 correction step is + /// saturating: a tap already at the coefficient format's rail is never + /// wrapped. Such a tap cannot absorb a step, so the correction stops + /// there and the row keeps its saturated sum — a row whose design exceeds + /// the format (|c| >= 2.0 in Q1.14 / Q1.30) had no representable DC gain + /// to preserve, and every tap of it is still the nearest representable + /// value. Design-time code: allocates a scratch vector for integer + /// formats, so keep it off the audio path. /// - /// \pre dst.size() == src.size() + /// @pre dst.size() == src.size() template inline void quantize_row_preserving_sum(std::span src, std::span::coeff> dst) { using tr = sample_traits; using coeff = typename tr::coeff; const std::size_t n = src.size(); - if constexpr (std::is_floating_point_v) { + if constexpr (!tr::k_is_fixed_point) { for (std::size_t t = 0; t < n; ++t) { dst[t] = tr::make_coeff(src[t]); } @@ -72,9 +78,14 @@ namespace tap::dsp { best = u; } } - dst[best] = static_cast(dst[best] + (residual > 0 ? 1 : -1)); + const std::int64_t step = residual > 0 ? 1 : -1; + const coeff stepped = detail::clamp_sat(static_cast(dst[best]) + step); + if (stepped == dst[best]) { + break; // at the rail: saturate, never wrap, and the residual is unabsorbable + } + dst[best] = stepped; remainder[best] -= sgn; - residual -= residual > 0 ? 1 : -1; + residual -= step; } } } diff --git a/include/tap/dsp/sample_traits.h b/include/tap/dsp/sample_traits.h index 7ad396e..c4d4bb6 100644 --- a/include/tap/dsp/sample_traits.h +++ b/include/tap/dsp/sample_traits.h @@ -1,5 +1,5 @@ /// @file sample_traits.h -/// @brief Sample-format customization point for FIR datapaths: float, Q15, Q31. +/// @brief Sample-format customization point for FIR datapaths: double, float, Q15, Q31. // SPDX-License-Identifier: MIT // Copyright 2026 Timothy Place and the DspTap contributors. // @@ -11,10 +11,18 @@ // dot products, and rounds/saturates back to samples. Engine-specific // extensions (e.g. SampleRateTap's inter-phase coefficient blending) layer on // top by deriving from these specializations and refining the concept. +// Butterfly arithmetic for the fixed-point FFT (multiply by a Q1.30 twiddle, +// saturating add/sub, block headroom) is the sibling trait in +// fft/fft_arith.h, built over the same constants. // -// Three sample types are provided: +// Four sample types are provided: // -// - float : float I/O and coefficients, double accumulation +// - double : double I/O, coefficients and accumulation — the golden +// model of every primitive (CLAUDE.md), so it is a sample +// format here too and a traits-based primitive can +// instantiate its golden profile +// - float : float I/O and coefficients, double accumulation — the +// embedded floating profile // - std::int16_t : Q15 samples, Q1.14 coefficients, int64 accumulation, // saturating output // - std::int32_t : Q31 samples, Q1.30 coefficients, int64 accumulation, @@ -35,12 +43,31 @@ // int16_t/int32_t also match what codecs and C ABIs actually deliver, and // keep memcpy/SIMD access patterns (e.g. the SMLALD pair loads in // fir_kernels.h) legal. +// +// The Q ladder is spelled as named constants on each fixed-point +// specialization, and the shifts are derived from them: +// +// k_sample_frac_bits fraction bits of the sample (Q0.15 -> 15, Q0.31 -> 31) +// k_coeff_frac_bits fraction bits of the stored coefficient (Q1.14 -> 14, +// Q1.30 -> 30): one headroom bit for the ~1.0 peak tap +// k_accum_pre_shift bits each product drops before it is accumulated +// (0 for Q15, 16 for Q31) +// k_accum_frac_bits = k_sample_frac_bits + k_coeff_frac_bits - k_accum_pre_shift +// (Q29 for Q15, Q45 for Q31): the accumulator's format +// k_finalize_shift = k_accum_frac_bits - k_sample_frac_bits +// = k_coeff_frac_bits - k_accum_pre_shift +// (14 for both profiles): the single rounding point, +// round-half-up then saturate +// +// Every member is constexpr and noexcept; the concept requires accum{} to be +// the additive identity, which is what lets the kernels start from a +// value-initialized accumulator. #pragma once -#include #include #include #include +#include namespace tap::dsp { @@ -77,6 +104,33 @@ namespace tap::dsp { template struct sample_traits; + /// Double datapath: double samples, coefficients and accumulation. The + /// golden model — every primitive's double profile is the reference the + /// float and fixed-point profiles are measured against, so double is a + /// sample format of the substrate and not only its design/accumulator + /// domain. Every conversion is the identity; nothing rounds or saturates. + template <> + struct sample_traits { + using coeff = double; ///< stored filter coefficient type + using accum = double; ///< dot-product accumulator type + + /// Coefficient units per 1.0: unity, no correction runs (quantize.h). + static constexpr double k_coeff_scale = 1.0; + static constexpr bool k_is_fixed_point = false; + + /// Convert a double-precision designed coefficient to storage form. + static constexpr coeff make_coeff(double c) noexcept { return c; } + + /// acc + x * c, in the accumulator domain. + static constexpr accum mac(accum acc, double x, coeff c) noexcept { return acc + x * c; } + + /// Convert the accumulator to an output sample. + static constexpr double finalize(accum acc) noexcept { return acc; } + + /// The zero/silence sample value. + static constexpr double silence() noexcept { return 0.0; } + }; + /// Float datapath: float samples and coefficients, double accumulation. /// The double accumulator keeps the dot-product noise floor far below a /// 120 dB transparency target; float coefficient storage quantizes the @@ -86,22 +140,24 @@ namespace tap::dsp { using coeff = float; ///< stored filter coefficient type using accum = double; ///< dot-product accumulator type - /// Convert a double-precision designed coefficient to storage form. - static coeff make_coeff(double c) noexcept { return static_cast(c); } /// Coefficient units per 1.0 (used by row-sum-preserving quantization; /// unity for floating storage, where no correction runs). - static constexpr double k_coeff_scale = 1.0; + static constexpr double k_coeff_scale = 1.0; + static constexpr bool k_is_fixed_point = false; + + /// Convert a double-precision designed coefficient to storage form. + static constexpr coeff make_coeff(double c) noexcept { return static_cast(c); } /// acc + x * c, in the accumulator domain. - static accum mac(accum acc, float x, coeff c) noexcept { + static constexpr accum mac(accum acc, float x, coeff c) noexcept { return acc + static_cast(x) * static_cast(c); } /// Convert the accumulator to an output sample (saturates for fixed point). - static float finalize(accum acc) noexcept { return static_cast(acc); } + static constexpr float finalize(accum acc) noexcept { return static_cast(acc); } /// The zero/silence sample value. - static float silence() noexcept { return 0.0f; } + static constexpr float silence() noexcept { return 0.0f; } }; // ANCHOR: st_q15_core @@ -121,28 +177,42 @@ namespace tap::dsp { using coeff = std::int16_t; using accum = std::int64_t; + static constexpr int k_sample_frac_bits = 15; ///< Q0.15 samples + static constexpr int k_coeff_frac_bits = 14; ///< Q1.14 coefficients + static constexpr int k_accum_pre_shift = 0; ///< products accumulate exactly + /// Q0.15 x Q1.14 = Q29, accumulated as is. + static constexpr int k_accum_frac_bits = k_sample_frac_bits + k_coeff_frac_bits - k_accum_pre_shift; + /// Q29 -> Q15: the single rounding, 14 bits. + static constexpr int k_finalize_shift = k_accum_frac_bits - k_sample_frac_bits; + static_assert(k_accum_frac_bits == 29 && k_finalize_shift == 14, + "Q15 ladder: Q0.15 x Q1.14 = Q29 -> Q15 is a 14-bit shift"); + // ANCHOR: st_q15_coeff - static coeff make_coeff(double c) noexcept { - return detail::round_sat(c * 16384.0); // Q1.14 + static constexpr double k_coeff_scale = + static_cast(std::int64_t{1} << k_coeff_frac_bits); ///< Q1.14 units per 1.0 + static constexpr bool k_is_fixed_point = true; + + static constexpr coeff make_coeff(double c) noexcept { + return detail::round_sat(c * k_coeff_scale); // Q1.14 } - static constexpr double k_coeff_scale = 16384.0; // Q1.14 units per 1.0 // ANCHOR_END: st_q15_coeff // ANCHOR: st_q15_mac - static accum mac(accum acc, std::int16_t x, coeff c) noexcept { + static constexpr accum mac(accum acc, std::int16_t x, coeff c) noexcept { return acc + static_cast(static_cast(x) * static_cast(c)); } // ANCHOR_END: st_q15_mac // ANCHOR: st_q15_finalize - static std::int16_t finalize(accum acc) noexcept { + static constexpr std::int16_t finalize(accum acc) noexcept { // Round-half-up, not half-even: the bias is a fraction of one // sub-LSB rounding step, far below the Q15 noise floor. - return detail::clamp_sat((acc + (1 << 13)) >> 14); // Q29 -> Q15 + return detail::clamp_sat((acc + (std::int64_t{1} << (k_finalize_shift - 1))) + >> k_finalize_shift); // Q29 -> Q15 } // ANCHOR_END: st_q15_finalize - static std::int16_t silence() noexcept { return 0; } + static constexpr std::int16_t silence() noexcept { return 0; } }; // ANCHOR_END: st_q15_core @@ -160,43 +230,86 @@ namespace tap::dsp { using coeff = std::int32_t; using accum = std::int64_t; - static coeff make_coeff(double c) noexcept { - return detail::round_sat(c * 1073741824.0); // Q1.30 + static constexpr int k_sample_frac_bits = 31; ///< Q0.31 samples + static constexpr int k_coeff_frac_bits = 30; ///< Q1.30 coefficients + static constexpr int k_accum_pre_shift = 16; ///< each Q61 product drops 16 bits before the add + /// Q0.31 x Q1.30 = Q61, pre-shifted to Q45 for accumulation. + static constexpr int k_accum_frac_bits = k_sample_frac_bits + k_coeff_frac_bits - k_accum_pre_shift; + /// Q45 -> Q31: the single rounding, 14 bits — the same shift as Q15's + /// Q29 -> Q15 because the pre-shift (16) equals the extra coefficient + /// precision (30 - 14). + static constexpr int k_finalize_shift = k_accum_frac_bits - k_sample_frac_bits; + static_assert(k_accum_frac_bits == 45 && k_finalize_shift == 14, + "Q31 ladder: Q0.31 x Q1.30 = Q61 >> 16 = Q45 -> Q31 is a 14-bit shift"); + + static constexpr double k_coeff_scale = + static_cast(std::int64_t{1} << k_coeff_frac_bits); ///< Q1.30 units per 1.0 + static constexpr bool k_is_fixed_point = true; + + static constexpr coeff make_coeff(double c) noexcept { + return detail::round_sat(c * k_coeff_scale); // Q1.30 } - static constexpr double k_coeff_scale = 1073741824.0; // Q1.30 units per 1.0 // ANCHOR: st_q31_mac - static accum mac(accum acc, std::int32_t x, coeff c) noexcept { - return acc + ((static_cast(x) * c) >> 16); // Q61 -> Q45 + static constexpr accum mac(accum acc, std::int32_t x, coeff c) noexcept { + return acc + ((static_cast(x) * c) >> k_accum_pre_shift); // Q61 -> Q45 } // ANCHOR_END: st_q31_mac - static std::int32_t finalize(accum acc) noexcept { - return detail::clamp_sat((acc + (1 << 13)) >> 14); // Q45 -> Q31 + static constexpr std::int32_t finalize(accum acc) noexcept { + return detail::clamp_sat((acc + (std::int64_t{1} << (k_finalize_shift - 1))) + >> k_finalize_shift); // Q45 -> Q31 } - static std::int32_t silence() noexcept { return 0; } + static constexpr std::int32_t silence() noexcept { return 0; } }; // ANCHOR_END: st_q31_core // ANCHOR: st_core_concept /// Satisfied by any type with a complete format-core sample_traits - /// specialization — everything the FIR kernels (fir_kernels.h) and table - /// builders require. Consumers with richer datapaths refine this concept - /// over their own trait extensions (SampleRateTap's sample_type adds the - /// inter-phase blend contract on top). + /// specialization — everything the FIR kernels (fir_kernels.h), the table + /// builders and quantize.h require: the coeff/accum types, the + /// k_coeff_scale / k_is_fixed_point properties, and the four operations. + /// Semantic requirement the kernels rely on and the static_asserts below + /// pin for the four shipping formats: a value-initialized accum{} is the + /// additive identity, so mac(accum{}, x, c) is the product x * c alone + /// and finalize(accum{}) is silence(). Consumers with richer datapaths + /// refine this concept over their own trait extensions (SampleRateTap's + /// sample_type adds the inter-phase blend contract on top). template concept sample_type = requires(T x, double d, typename sample_traits::accum a, typename sample_traits::coeff c) { + typename sample_traits::coeff; + typename sample_traits::accum; + requires std::is_same_v::k_coeff_scale), const double>; + requires std::is_same_v::k_is_fixed_point), const bool>; { sample_traits::make_coeff(d) } -> std::same_as::coeff>; { sample_traits::mac(a, x, c) } -> std::same_as::accum>; { sample_traits::finalize(a) } -> std::same_as; { sample_traits::silence() } -> std::same_as; }; + static_assert(sample_type); static_assert(sample_type); static_assert(sample_type); static_assert(sample_type); + + // accum{} is the additive identity: finalize(accum{}) is silence, and a + // single mac from accum{} yields the bare product in the accumulator + // domain (Q29 for Q15, Q45 for Q31, the exact double product for the + // floating formats). + static_assert(sample_traits::finalize(sample_traits::accum{}) == 0.0); + static_assert(sample_traits::mac(sample_traits::accum{}, 0.5, 0.5) == 0.25); + static_assert(sample_traits::finalize(sample_traits::accum{}) == 0.0f); + static_assert(sample_traits::mac(sample_traits::accum{}, 0.5f, 0.5f) == 0.25); + static_assert(sample_traits::finalize(sample_traits::accum{}) == 0); + static_assert(sample_traits::mac(sample_traits::accum{}, std::int16_t{32767}, + std::int16_t{16384}) + == std::int64_t{32767} * 16384); + static_assert(sample_traits::finalize(sample_traits::accum{}) == 0); + static_assert(sample_traits::mac(sample_traits::accum{}, std::int32_t{1} << 30, + std::int32_t{1} << 30) + == std::int64_t{1} << 44); // ANCHOR_END: st_core_concept } // namespace tap::dsp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index eaaa14b..ec83dd8 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -125,6 +125,7 @@ tap_dsp_add_gtest_executable(tap_dsp_tests test_analysis.cpp test_decimate.cpp test_fft.cpp + test_fft_arith.cpp test_fft_backend.cpp test_fir_kernels.cpp test_kaiser.cpp diff --git a/tests/test_decimate.cpp b/tests/test_decimate.cpp index 141bb40..03860e6 100644 --- a/tests/test_decimate.cpp +++ b/tests/test_decimate.cpp @@ -2,18 +2,21 @@ // Copyright 2026 Timothy Place and the DspTap contributors. // // Locks down the tap::dsp::basic_decimator contract: the pinned tap counts, -// sample-for-sample agreement of the float golden model with the committed -// numpy reference for every ratio and profile, chunking invariance, exact -// output counting, unity DC gain (exact in Q15 through row-sum-preserving +// sample-for-sample agreement of the float profile with the committed numpy +// reference for every ratio and profile (and of the double golden model to +// the reference's own float32 rounding), chunking invariance, exact output +// counting, unity DC gain (exact in Q15 through row-sum-preserving // quantization), the stated passband/stopband numbers measured from the -// shipped coefficients, the integer group delay, and the Q15 profile -// tracking float within the format's own floor. +// shipped coefficients, the integer group delay, and the Q15 / Q31 profiles +// tracking the double golden model within each format's own floor. The +// float/double contract tests are typed over both profiles. #include #include #include #include #include +#include #include #include @@ -26,6 +29,11 @@ namespace { using tap::dsp::basic_decimator; using tap::dsp::decimate_profile; + template + class decimate_test : public ::testing::Test {}; + using floating_types = ::testing::Types; + TYPED_TEST_SUITE(decimate_test, floating_types, ); + // Float coefficient storage against float64 reference coefficients: the // same -90 dB floor RatioTap pins at 3e-5 on a 0.9 peak. constexpr float k_reference_tolerance = 3e-5f; @@ -66,19 +74,50 @@ namespace { k_dec_transparent_6.size()); } - TEST(Decimate, ChunkingIsBitIdentical) { - const auto& x = frontend_ref::k_dec_input; - basic_decimator whole; - std::vector ref(whole.outputs_for(x.size())); - whole.process(x.data(), x.size(), ref.data()); + // The double golden model against the numpy float64 reference, which is + // committed rounded to float32: the disagreement is the reference's own + // rounding. Measured 2026-09 at 5.75e-8 over every ratio and profile + // (half a float ulp at the signal's level); pinned at 2x. + TEST(Decimate, DoubleMatchesNumpyReferenceToItsFloatRounding) { + using namespace frontend_ref; + const auto check = [](auto ratio, const decimate_profile& p, const float* ref, std::size_t n_ref) { + constexpr std::size_t M = decltype(ratio)::value; + basic_decimator dec(p); + const std::vector x(k_dec_input.begin(), k_dec_input.end()); + std::vector y(n_ref); + ASSERT_EQ(dec.process(x.data(), x.size(), y.data()), n_ref); + double worst = 0.0; + for (std::size_t i = 0; i < n_ref; ++i) { + worst = std::max(worst, std::abs(y[i] - static_cast(ref[i]))); + } + EXPECT_LT(worst, 1.15e-7) << "ratio " << M << " taps " << dec.taps(); + }; + using two = std::integral_constant; + using three = std::integral_constant; + using six = std::integral_constant; + check(two{}, decimate_profile::economy(), k_dec_economy_2.data(), k_dec_economy_2.size()); + check(three{}, decimate_profile::economy(), k_dec_economy_3.data(), k_dec_economy_3.size()); + check(six{}, decimate_profile::economy(), k_dec_economy_6.data(), k_dec_economy_6.size()); + check(two{}, decimate_profile::transparent(), k_dec_transparent_2.data(), k_dec_transparent_2.size()); + check(three{}, decimate_profile::transparent(), k_dec_transparent_3.data(), k_dec_transparent_3.size()); + check(six{}, decimate_profile::transparent(), k_dec_transparent_6.data(), k_dec_transparent_6.size()); + } + + TYPED_TEST(decimate_test, ChunkingIsBitIdentical) { + using sample = TypeParam; + const auto& x = frontend_ref::k_dec_input; + const std::vector xs(x.begin(), x.end()); + basic_decimator whole; + std::vector ref(whole.outputs_for(xs.size())); + whole.process(xs.data(), xs.size(), ref.data()); for (const std::size_t chunk : {std::size_t{1}, std::size_t{2}, std::size_t{7}, std::size_t{255}, std::size_t{1000}}) { - basic_decimator dec; - std::vector y; - for (std::size_t i = 0; i < x.size(); i += chunk) { - const std::size_t n = std::min(chunk, x.size() - i); - std::vector part(dec.outputs_for(n)); - ASSERT_EQ(dec.process(x.data() + i, n, part.data()), part.size()); + basic_decimator dec; + std::vector y; + for (std::size_t i = 0; i < xs.size(); i += chunk) { + const std::size_t n = std::min(chunk, xs.size() - i); + std::vector part(dec.outputs_for(n)); + ASSERT_EQ(dec.process(xs.data() + i, n, part.data()), part.size()); y.insert(y.end(), part.begin(), part.end()); } ASSERT_EQ(y.size(), ref.size()) << "chunk " << chunk; @@ -103,13 +142,16 @@ namespace { EXPECT_EQ(dec.outputs_for(6), 2U); // x[9] too } - TEST(Decimate, DcGainIsUnity) { - basic_decimator dec; - std::vector x(dec.taps() * 2, 0.5f); - std::vector y(dec.outputs_for(x.size())); + TYPED_TEST(decimate_test, DcGainIsUnity) { + using sample = TypeParam; + basic_decimator dec; + std::vector x(dec.taps() * 2, sample{0.5}); + std::vector y(dec.outputs_for(x.size())); dec.process(x.data(), x.size(), y.data()); - EXPECT_NEAR(y.back(), 0.5f, 1e-6f); + EXPECT_NEAR(y.back(), 0.5, 1e-6); + } + TEST(Decimate, DcGainIsExactInQ15) { basic_decimator q15; std::vector xq(q15.taps() * 2, std::int16_t{16384}); std::vector yq(q15.outputs_for(xq.size())); @@ -154,42 +196,60 @@ namespace { expect_profile_met<6>(decimate_profile::transparent()); } - TEST(Decimate, GroupDelayIsHalfTheTaps) { - basic_decimator dec; + TYPED_TEST(decimate_test, GroupDelayIsHalfTheTaps) { + using sample = TypeParam; + basic_decimator dec; EXPECT_EQ(dec.latency_input_samples(), (dec.taps() - 1) / 2); - std::vector x(dec.taps() + 3, 0.0f); - x[0] = 1.0f; - std::vector y(dec.outputs_for(x.size())); + std::vector x(dec.taps() + 3, sample{0}); + x[0] = sample{1}; + std::vector y(dec.outputs_for(x.size())); dec.process(x.data(), x.size(), y.data()); const auto peak = static_cast(std::max_element(y.begin(), y.end()) - y.begin()); // y[k] = h[3k]; the centre tap (index 60) lands on output 20 exactly. EXPECT_EQ(peak * 3, dec.latency_input_samples()); - EXPECT_FLOAT_EQ(y[peak], dec.coefficients()[dec.latency_input_samples()]); - } - - TEST(Decimate, Q15TracksFloatWithinTheFormatFloor) { - // The input is scaled to half full scale: at 0.9 the float output - // overshoots 1.0 on this noise and Q15 saturates, which is the format's - // contract, not tracking error. Measured 2026-09 at 1.6e-4 (Q1.14 - // coefficient rounding over 121 taps, plus one output LSB); pinned at 2x. - const auto& x = frontend_ref::k_dec_input; - basic_decimator ff; - basic_decimator fq; - std::vector xh(x.size()); - std::vector xq(x.size()); + EXPECT_EQ(y[peak], dec.coefficients()[dec.latency_input_samples()]); + } + + // The fixed-point profiles against the double golden model (the house + // rule: cross-precision compares to double, never to a sibling profile). + // The input is scaled to half full scale: at 0.9 the output overshoots + // 1.0 on this noise and the fixed formats saturate, which is their + // contract, not tracking error. Returns the worst deviation in full-scale + // units. + template + double worst_deviation_from_double(double full_scale) { + const auto& x = frontend_ref::k_dec_input; + basic_decimator fd; + basic_decimator fq; + std::vector xd(x.size()); + std::vector xq(x.size()); for (std::size_t i = 0; i < x.size(); ++i) { - xh[i] = 0.5f * x[i]; - xq[i] = static_cast(std::lround(xh[i] * 32768.0f)); + xd[i] = 0.5 * static_cast(x[i]); + xq[i] = static_cast(std::llround(xd[i] * full_scale)); } - std::vector yf(ff.outputs_for(xh.size())); - std::vector yq(fq.outputs_for(xq.size())); - ff.process(xh.data(), xh.size(), yf.data()); + std::vector yd(fd.outputs_for(xd.size())); + std::vector yq(fq.outputs_for(xq.size())); + fd.process(xd.data(), xd.size(), yd.data()); fq.process(xq.data(), xq.size(), yq.data()); double worst = 0.0; - for (std::size_t i = 0; i < yf.size(); ++i) { - worst = std::max(worst, std::abs(static_cast(yq[i]) / 32768.0 - static_cast(yf[i]))); + for (std::size_t i = 0; i < yd.size(); ++i) { + worst = std::max(worst, std::abs(static_cast(yq[i]) / full_scale - yd[i])); } - EXPECT_LT(worst, 3.2e-4) << "Q15 vs float, in full-scale units"; + return worst; + } + + TEST(Decimate, Q15TracksDoubleWithinTheFormatFloor) { + // Measured 2026-09 at 1.57e-4 (Q1.14 coefficient rounding over 121 + // taps, plus one output LSB); pinned at 2x. The same number as the + // former Q15-vs-float pin: float sits 2.5e-8 from double here. + EXPECT_LT(worst_deviation_from_double(32768.0), 3.2e-4) << "Q15 vs double, in full-scale units"; + } + + TEST(Decimate, Q31TracksDoubleWithinTheFormatFloor) { + // Measured 2026-09 at 2.45e-9 (Q1.30 coefficient rounding over 121 + // taps, the 16-bit product pre-shift, one output LSB); pinned at 2x. + EXPECT_LT(worst_deviation_from_double(2147483648.0), 4.9e-9) + << "Q31 vs double, in full-scale units"; } TEST(Decimate, Q15SaturatesInsteadOfWrapping) { diff --git a/tests/test_fft_arith.cpp b/tests/test_fft_arith.cpp new file mode 100644 index 0000000..717d5ea --- /dev/null +++ b/tests/test_fft_arith.cpp @@ -0,0 +1,249 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Contract battery for the butterfly arithmetic trait (fft/fft_arith.h), +// the specification the fixed-point FFT kernel is written against. Every +// number here is a documented contract point: the Q1.30 twiddle format, the +// single rounding point of the multiply, saturation at the int32 rails, the +// Q15 profile's two guard bits, and the block-headroom definition. + +#include +#include +#include +#include +#include + +#include + +#include "tap/dsp/fft/fft_arith.h" + +namespace { + + using q31 = tap::dsp::fft_arith; + using q15 = tap::dsp::fft_arith; + using f32 = tap::dsp::fft_arith; + using f64 = tap::dsp::fft_arith; + + constexpr std::int32_t k_i32_max = std::numeric_limits::max(); + constexpr std::int32_t k_i32_min = std::numeric_limits::min(); + constexpr std::int16_t k_i16_max = std::numeric_limits::max(); + constexpr std::int16_t k_i16_min = std::numeric_limits::min(); + + TEST(FftArith, TwiddleUnityIsRepresentable) { + // Q1.30 for both fixed profiles, one type, 1.0 = 2^30 exactly. + static_assert(std::is_same_v); + static_assert(std::is_same_v); + static_assert(q31::k_coeff_frac_bits == 30 && q15::k_coeff_frac_bits == 30); + EXPECT_EQ(q31::make_coeff(1.0), 1073741824); + EXPECT_EQ(q15::make_coeff(1.0), 1073741824); + EXPECT_EQ(q31::make_coeff(-1.0), -1073741824); + EXPECT_EQ(q31::make_coeff(0.5), 536870912); + EXPECT_EQ(q31::k_coeff_one, 1073741824); + // Rounded once, half away from zero, saturating (round_sat). + EXPECT_EQ(q31::make_coeff(0x1p-31), 1); // 0.5 LSB rounds up + EXPECT_EQ(q31::make_coeff(-0x1p-31), -1); + EXPECT_EQ(q31::make_coeff(3.0), k_i32_max); + // A unity twiddle is the identity on any sample. + EXPECT_EQ(q31::mul_coeff(123456789, q31::k_coeff_one), 123456789); + EXPECT_EQ(q31::mul_coeff(k_i32_max, q31::k_coeff_one), k_i32_max); + EXPECT_EQ(q31::mul_coeff(k_i32_min, q31::k_coeff_one), k_i32_min); + } + + TEST(FftArith, MulCoeffIsExactOnRepresentableProducts) { + // x * w in Qm.n x Q1.30: when the exact product has no bits below the + // shift, the result is the exact product (hand-computed). + EXPECT_EQ(q31::mul_coeff(1 << 20, 1 << 29), 1 << 19); // 2^20 * 0.5 + EXPECT_EQ(q31::mul_coeff(-(1 << 20), 1 << 29), -(1 << 19)); + EXPECT_EQ(q31::mul_coeff(1 << 20, -(1 << 29)), -(1 << 19)); + EXPECT_EQ(q31::mul_coeff(3 << 10, 1 << 28), 3 << 8); // 3072 * 0.25 = 768 + EXPECT_EQ(q31::mul_coeff(1000, 0), 0); + EXPECT_EQ(q31::mul_coeff(0, 1 << 30), 0); + EXPECT_EQ(q31::mul_coeff(1 << 30, 1 << 30), 1 << 30); // 0.5 * 1.0 in Q0.31 + EXPECT_EQ(q31::mul_coeff(1 << 30, 3 << 28), 3 << 28); // 0.5 * 0.75 = 0.375 + // The float and double profiles are the plain product. + EXPECT_EQ(f32::mul_coeff(0.5f, 0.75f), 0.375f); + EXPECT_EQ(f64::mul_coeff(0.5, 0.75), 0.375); + } + + TEST(FftArith, MulCoeffRoundsHalfUp) { + // 1 * 0.5 in Q1.30 = 2^29 in the product: exactly half an LSB, which + // rounds up to 1 for the positive product and to 0 for the negative + // product (half-up, not half-away). + EXPECT_EQ(q31::mul_coeff(1, 1 << 29), 1); + EXPECT_EQ(q31::mul_coeff(-1, 1 << 29), 0); + EXPECT_EQ(q31::mul_coeff(1, -(1 << 29)), 0); + // Just below half rounds down; just above rounds up, on both sides. + EXPECT_EQ(q31::mul_coeff(1, (1 << 29) - 1), 0); + EXPECT_EQ(q31::mul_coeff(1, (1 << 29) + 1), 1); + EXPECT_EQ(q31::mul_coeff(-1, (1 << 29) - 1), 0); + EXPECT_EQ(q31::mul_coeff(-1, (1 << 29) + 1), -1); + // 3 * 0.5 = 1.5 -> 2; -3 * 0.5 = -1.5 -> -1. + EXPECT_EQ(q31::mul_coeff(3, 1 << 29), 2); + EXPECT_EQ(q31::mul_coeff(-3, 1 << 29), -1); + // One rounding, not two: 2^30 - 1 (0.999999999) times 2^30 - 1. + // Exact product = (2^30 - 1)^2 = 2^60 - 2^31 + 1; >> 30 with half-up: + // floor((2^60 - 2^31 + 1 + 2^29) / 2^30) = 2^30 - 2 + floor((1 + 2^29) / 2^30) = 2^30 - 2. + EXPECT_EQ(q31::mul_coeff((1 << 30) - 1, (1 << 30) - 1), (1 << 30) - 2); + } + + TEST(FftArith, MulCoeffSaturatesAtTheRail) { + // With |w| <= 1.0 the only saturating input is INT32_MIN * -1.0, whose + // exact product is +2^31. + EXPECT_EQ(q31::mul_coeff(k_i32_min, -q31::k_coeff_one), k_i32_max); + EXPECT_EQ(q31::mul_coeff(k_i32_min, -q31::k_coeff_one + 1), k_i32_max - 1); // exact 2^31 - 2 + 2^-30: fits + EXPECT_EQ(q31::mul_coeff(k_i32_max, -q31::k_coeff_one), k_i32_min + 1); + EXPECT_EQ(q31::mul_coeff(k_i32_min, q31::k_coeff_one), k_i32_min); + // Out-of-contract twiddles (|w| > 1.0) saturate rather than wrap. + EXPECT_EQ(q31::mul_coeff(k_i32_max, k_i32_max), k_i32_max); + EXPECT_EQ(q31::mul_coeff(k_i32_min, k_i32_max), k_i32_min); + EXPECT_EQ(q31::mul_coeff(k_i32_min, k_i32_min), k_i32_max); + } + + TEST(FftArith, AddSubSaturateAtTheRails) { + EXPECT_EQ(q31::add(1, 2), 3); + EXPECT_EQ(q31::sub(1, 2), -1); + EXPECT_EQ(q31::add(k_i32_max, 1), k_i32_max); + EXPECT_EQ(q31::add(k_i32_max, k_i32_max), k_i32_max); + EXPECT_EQ(q31::add(k_i32_min, -1), k_i32_min); + EXPECT_EQ(q31::add(k_i32_min, k_i32_min), k_i32_min); + EXPECT_EQ(q31::sub(k_i32_min, 1), k_i32_min); + EXPECT_EQ(q31::sub(k_i32_max, -1), k_i32_max); + EXPECT_EQ(q31::sub(0, k_i32_min), k_i32_max); // -INT32_MIN does not exist + EXPECT_EQ(q31::sub(k_i32_min, k_i32_max), k_i32_min); + // In range, the plain results. + EXPECT_EQ(q31::add(k_i32_max, k_i32_min), -1); + EXPECT_EQ(q31::sub(k_i32_min, k_i32_min), 0); + EXPECT_EQ(f64::add(0.5, 0.25), 0.75); + EXPECT_EQ(f32::sub(0.5f, 0.25f), 0.25f); + } + + TEST(FftArith, ShrRoundRoundsHalfUp) { + EXPECT_EQ(q31::shr_round(8, 2), 2); + EXPECT_EQ(q31::shr_round(2, 2), 1); // 0.5 -> 1 + EXPECT_EQ(q31::shr_round(-2, 2), 0); // -0.5 -> 0 (half-up) + EXPECT_EQ(q31::shr_round(1, 2), 0); // 0.25 -> 0 + EXPECT_EQ(q31::shr_round(-1, 2), 0); // -0.25 -> 0 + EXPECT_EQ(q31::shr_round(3, 2), 1); // 0.75 -> 1 + EXPECT_EQ(q31::shr_round(-3, 2), -1); // -0.75 -> -1 + EXPECT_EQ(q31::shr_round(-6, 2), -1); // -1.5 -> -1 + EXPECT_EQ(q31::shr_round(6, 2), 2); // 1.5 -> 2 + EXPECT_EQ(q31::shr_round(7, 0), 7); // identity + // The rounding add cannot overflow: full scale by one bit. + EXPECT_EQ(q31::shr_round(k_i32_max, 1), 1 << 30); // (2^31 - 1 + 1) / 2 + EXPECT_EQ(q31::shr_round(k_i32_min, 1), -(1 << 30)); + EXPECT_EQ(q31::shr_round(k_i32_min, 31), -1); + EXPECT_EQ(q31::shr_round(k_i32_max, 31), 1); // 0.99999 -> 1 + EXPECT_EQ(q31::shr_round(k_i32_min + 1, 31), -1); // -0.99999 -> -1 + // Floating: exact power-of-two scale. + EXPECT_EQ(f64::shr_round(3.0, 2), 0.75); + EXPECT_EQ(f32::shr_round(-3.0f, 1), -1.5f); + EXPECT_EQ(f64::shr_round(1.0, 0), 1.0); + } + + TEST(FftArith, HeadroomBitsOnZeroOneAndFullScale) { + const std::array zeros{0, 0, 0, 0}; + EXPECT_EQ(q31::headroom_bits(zeros.data(), zeros.size()), 31); + EXPECT_EQ(q31::headroom_bits(zeros.data(), 0), 31); // empty block + const std::array one{1}; + EXPECT_EQ(q31::headroom_bits(one.data(), one.size()), 30); + const std::array minus_one{-1}; + EXPECT_EQ(q31::headroom_bits(minus_one.data(), 1), 31); // -1 << 31 still fits + const std::array minus_two{-2}; + EXPECT_EQ(q31::headroom_bits(minus_two.data(), 1), 30); + const std::array full{0, k_i32_max, 0}; + EXPECT_EQ(q31::headroom_bits(full.data(), full.size()), 0); + const std::array full_neg{0, k_i32_min, 0}; + EXPECT_EQ(q31::headroom_bits(full_neg.data(), full_neg.size()), 0); + const std::array half{1 << 30, -(1 << 30)}; + EXPECT_EQ(q31::headroom_bits(half.data(), half.size()), 0); // 2^30 << 1 = 2^31 does not fit + const std::array quarter{(1 << 29) + 5, -(1 << 29)}; + EXPECT_EQ(q31::headroom_bits(quarter.data(), quarter.size()), 1); + // The definition, checked on the whole bit ladder: x = 1 << (30 - s) + // has headroom s, and -x has s + 1 (two's complement holds one more + // negative power of two: -2^k << (31 - k) is INT32_MIN and fits). + for (int s = 0; s <= 30; ++s) { + const std::int32_t x = std::int32_t{1} << (30 - s); + EXPECT_EQ(q31::headroom_bits(&x, 1), s) << "s=" << s; + const std::int32_t nx = -x; + EXPECT_EQ(q31::headroom_bits(&nx, 1), s + 1) << "s=" << s; + } + // Floating blocks report no headroom: no block scaling. + const std::array fd{0.5, 1e-9}; + EXPECT_EQ(f64::headroom_bits(fd.data(), fd.size()), 0); + } + + TEST(FftArith, WidenPlacesTwoGuardBits) { + static_assert(q15::k_guard_bits == 2 && q15::k_widen_shift == 14); + static_assert(std::is_same_v); + static_assert(std::is_same_v); + EXPECT_EQ(q15::widen(1), 1 << 14); + EXPECT_EQ(q15::widen(k_i16_max), 32767 << 14); + EXPECT_EQ(q15::widen(k_i16_min), -(1 << 29)); + // Full scale sits two bits below the int32 sign: headroom exactly 2. + const std::int32_t fs = q15::widen(k_i16_max); + EXPECT_EQ(q31::headroom_bits(&fs, 1), 2); + const std::int32_t nfs = q15::widen(k_i16_min); + EXPECT_EQ(q31::headroom_bits(&nfs, 1), 2); + // 4x growth from full scale fits (the fixed-scaling argument). + EXPECT_EQ(q31::add(q31::add(nfs, nfs), q31::add(nfs, nfs)), k_i32_min); + EXPECT_EQ(q31::add(q31::add(fs, fs), q31::add(fs, fs)), 4 * (32767 << 14)); + } + + TEST(FftArith, NarrowRoundsHalfUpAndSaturates) { + EXPECT_EQ(q15::narrow(1 << 14), 1); + EXPECT_EQ(q15::narrow(1 << 13), 1); // exactly half rounds up + EXPECT_EQ(q15::narrow((1 << 13) - 1), 0); + EXPECT_EQ(q15::narrow(-(1 << 13)), 0); // -0.5 LSB -> 0 (half-up) + EXPECT_EQ(q15::narrow(-(1 << 13) - 1), -1); + EXPECT_EQ(q15::narrow(k_i32_max), k_i16_max); + EXPECT_EQ(q15::narrow(k_i32_min), k_i16_min); + EXPECT_EQ(q15::narrow(32768 << 14), k_i16_max); // one LSB over full scale saturates + EXPECT_EQ(q15::narrow((32767 << 14) + (1 << 13)), k_i16_max); + EXPECT_EQ(q15::narrow(-(32769 << 14)), k_i16_min); + EXPECT_EQ(q15::narrow(-(1 << 29)), k_i16_min); + } + + TEST(FftArith, WidenNarrowRoundTripsEveryInt16) { + for (int v = k_i16_min; v <= k_i16_max; ++v) { + const auto x = static_cast(v); + ASSERT_EQ(q15::narrow(q15::widen(x)), x) << v; + } + // And the identities of the other profiles. + EXPECT_EQ(q31::narrow(q31::widen(k_i32_min)), k_i32_min); + EXPECT_EQ(f32::narrow(f32::widen(0.5f)), 0.5f); + EXPECT_EQ(f64::narrow(f64::widen(-0.25)), -0.25); + } + + TEST(FftArith, EveryOperationIsConstexprAndNoexcept) { + static_assert(noexcept(q31::mul_coeff(1, 1)) && noexcept(q31::add(1, 1)) && noexcept(q31::sub(1, 1))); + static_assert(noexcept(q31::shr_round(1, 1)) && noexcept(q31::headroom_bits(nullptr, 0))); + static_assert(noexcept(q15::widen(1)) && noexcept(q15::narrow(1)) && noexcept(q31::make_coeff(1.0))); + static_assert(q31::mul_coeff(1 << 20, 1 << 29) == 1 << 19); + static_assert(q31::add(k_i32_max, 1) == k_i32_max); + static_assert(q31::shr_round(-2, 2) == 0); + static_assert(q15::narrow(q15::widen(-12345)) == -12345); + constexpr std::array block{1 << 20, -(1 << 24)}; + static_assert(q31::headroom_bits(block.data(), block.size()) == 7); // -2^24 << 7 = INT32_MIN + static_assert(f64::shr_round(1.0, 3) == 0.125); + } + + // A butterfly on the widened Q15 profile, hand-computed end to end: the + // arithmetic the kernel will compose, checked once here so the pieces + // are known to fit together (formats, rounding point, narrowing). + TEST(FftArith, WidenedButterflyHandComputed) { + // x0 = 0.5, x1 = 0.25 in Q0.15; w = cos(60 deg) = 0.5 in Q1.30. + const std::int32_t a = q15::widen(16384); + const std::int32_t b = q15::widen(8192); + const auto w = q15::make_coeff(0.5); + const std::int32_t t = q31::mul_coeff(b, w); // 0.125 in Q2.29 + EXPECT_EQ(t, 4096 << 14); + const std::int32_t y0 = q31::add(a, t); // 0.625 + const std::int32_t y1 = q31::sub(a, t); // 0.375 + EXPECT_EQ(q15::narrow(y0), 20480); + EXPECT_EQ(q15::narrow(y1), 12288); + // Scaled by one bit before narrowing (fixed scaling): half those. + EXPECT_EQ(q15::narrow(q31::shr_round(y0, 1)), 10240); + EXPECT_EQ(q15::narrow(q31::shr_round(y1, 1)), 6144); + } + +} // namespace diff --git a/tests/test_fir_kernels.cpp b/tests/test_fir_kernels.cpp index 4f13c01..baace59 100644 --- a/tests/test_fir_kernels.cpp +++ b/tests/test_fir_kernels.cpp @@ -5,7 +5,8 @@ // extracted from SampleRateTap's multichannel suite, where it gates the // frame-major fast path — is bit-exactness: the channel-parallel kernel must // produce the identical bits as the planar dot for every sample type, because -// consumers switch between the layouts by channel count and target. +// consumers switch between the layouts by channel count and target. Typed +// over the four substrate formats; double is the golden model. #include #include @@ -33,6 +34,10 @@ namespace { template S sample_from(std::uint32_t r); template <> + double sample_from(std::uint32_t r) { + return (static_cast(r % 65536) - 32768.0) / 65536.0; + } + template <> float sample_from(std::uint32_t r) { return (static_cast(r % 65536) - 32768.0f) / 65536.0f; } @@ -53,7 +58,7 @@ namespace { template class fir_kernels_test : public ::testing::Test {}; - using sample_types = ::testing::Types; + using sample_types = ::testing::Types; TYPED_TEST_SUITE(fir_kernels_test, sample_types, ); // dot_row must equal the reference accumulation: mac per tap in order, diff --git a/tests/test_quantize.cpp b/tests/test_quantize.cpp index df3a9d3..9055921 100644 --- a/tests/test_quantize.cpp +++ b/tests/test_quantize.cpp @@ -64,6 +64,38 @@ namespace { check_rows_sum_exact(); } + TEST(Quantize, DoubleIsPlainConversion) { + // The golden model's coefficients are the designed values themselves. + const std::vector row{0.25, -0.125, 1.0, -0.9999, 0.0, 1e-300}; + std::vector q(row.size()); + quantize_row_preserving_sum(row, q); + EXPECT_EQ(q, row); + } + + // The +/-1 correction steps saturate. A row whose design exceeds the + // coefficient format has a tap at the rail with a huge positive remainder, + // so the largest-remainder pick lands on it: the step must leave it at + // the rail, never wrap it to the opposite sign, and the routine returns. + TEST(Quantize, CorrectionNeverWrapsATapAtTheRail) { + const std::vector row{3.0, 0.4, -0.25}; + std::vector q15(row.size()); + quantize_row_preserving_sum(row, q15); + EXPECT_EQ(q15[0], 32767); // saturated, not wrapped + EXPECT_EQ(q15[1], 6554); // 0.4 * 16384 = 6553.6, its own rounding, untouched + EXPECT_EQ(q15[2], -4096); + std::vector q31(row.size()); + quantize_row_preserving_sum(row, q31); + EXPECT_EQ(q31[0], 2147483647); + EXPECT_EQ(q31[1], 429496730); // 0.4 * 2^30 = 429496729.6 + EXPECT_EQ(q31[2], -268435456); + // The negative rail too. + const std::vector neg{-3.0, -0.4}; + std::vector n15(neg.size()); + quantize_row_preserving_sum(neg, n15); + EXPECT_EQ(n15[0], -32768); + EXPECT_EQ(n15[1], -6554); + } + TEST(Quantize, FloatIsPlainConversion) { const std::vector row{0.25, -0.125, 1.0, -0.9999, 0.0}; std::vector q(row.size()); diff --git a/tests/test_sample_traits.cpp b/tests/test_sample_traits.cpp index 57aff55..de0b239 100644 --- a/tests/test_sample_traits.cpp +++ b/tests/test_sample_traits.cpp @@ -8,6 +8,7 @@ #include #include +#include #include @@ -18,6 +19,7 @@ namespace { using q15 = tap::dsp::sample_traits; using q31 = tap::dsp::sample_traits; using f32 = tap::dsp::sample_traits; + using f64 = tap::dsp::sample_traits; TEST(SampleTraits, CoefficientConversionRoundsAndSaturates) { EXPECT_EQ(q15::make_coeff(0.0), 0); @@ -28,6 +30,33 @@ namespace { EXPECT_EQ(q31::make_coeff(1.0), 1073741824); // Q1.30 EXPECT_EQ(q31::make_coeff(10.0), 2147483647); // saturates EXPECT_FLOAT_EQ(f32::make_coeff(0.25), 0.25f); // unity scale, plain cast + EXPECT_EQ(f64::make_coeff(0.1), 0.1); // identity + } + + TEST(SampleTraits, CoefficientScaleIsTheFormatsUnity) { + // k_coeff_scale is derived from the named fraction-bit constants and + // must be the number make_coeff maps 1.0 to; quantize.h multiplies by + // one and compares against the other. + static_assert(q15::k_coeff_frac_bits == 14); + static_assert(q31::k_coeff_frac_bits == 30); + EXPECT_EQ(q15::k_coeff_scale, 16384.0); + EXPECT_EQ(q31::k_coeff_scale, 1073741824.0); + EXPECT_EQ(static_cast(q15::make_coeff(1.0)), q15::k_coeff_scale); + EXPECT_EQ(static_cast(q31::make_coeff(1.0)), q31::k_coeff_scale); + EXPECT_EQ(f32::k_coeff_scale, 1.0); + EXPECT_EQ(f64::k_coeff_scale, 1.0); + } + + TEST(SampleTraits, QLadderClosesFromTheNamedConstants) { + // The shifts are derived, not written: accumulator format = sample + + // coefficient fraction bits - pre-shift; finalize = accumulator -> + // sample. Both profiles land on the documented 14-bit rounding. + static_assert(q15::k_sample_frac_bits == 15 && q15::k_accum_pre_shift == 0); + static_assert(q15::k_accum_frac_bits == 29 && q15::k_finalize_shift == 14); + static_assert(q31::k_sample_frac_bits == 31 && q31::k_accum_pre_shift == 16); + static_assert(q31::k_accum_frac_bits == 45 && q31::k_finalize_shift == 14); + static_assert(q15::k_is_fixed_point && q31::k_is_fixed_point); + static_assert(!f32::k_is_fixed_point && !f64::k_is_fixed_point); } TEST(SampleTraits, RoundSatRoundsHalfAwayFromZero) { @@ -56,9 +85,11 @@ namespace { EXPECT_EQ(q15::finalize((std::int64_t{1} << 13) - 1), 0); EXPECT_EQ(q15::finalize(-(std::int64_t{1} << 13)), 0); EXPECT_EQ(q15::finalize(-(std::int64_t{1} << 13) - 1), -1); - // Q45 -> Q31 uses the identical constant. + // Q45 -> Q31 uses the identical constant, on both sides of zero. EXPECT_EQ(q31::finalize((std::int64_t{1} << 13)), 1); EXPECT_EQ(q31::finalize((std::int64_t{1} << 13) - 1), 0); + EXPECT_EQ(q31::finalize(-(std::int64_t{1} << 13)), 0); + EXPECT_EQ(q31::finalize(-(std::int64_t{1} << 13) - 1), -1); } TEST(SampleTraits, Q15MacIsExactInInt64) { @@ -94,18 +125,42 @@ namespace { EXPECT_DOUBLE_EQ(acc, 1.0 + 0x1p-30); } + TEST(SampleTraits, DoubleIsTheIdentityDatapath) { + static_assert(std::is_same_v); + static_assert(std::is_same_v); + // Nothing rounds, nothing saturates: the golden model passes values + // through untouched, so a traits-based primitive's double profile is + // the same arithmetic as its hand-written double sibling. + EXPECT_EQ(f64::make_coeff(1e300), 1e300); + EXPECT_EQ(f64::finalize(1e300), 1e300); + EXPECT_EQ(f64::mac(1.0, 0x1p-40, 1.0), 1.0 + 0x1p-40); + EXPECT_EQ(f64::silence(), 0.0); + } + TEST(SampleTraits, SilenceIsZero) { EXPECT_EQ(q15::silence(), 0); EXPECT_EQ(q31::silence(), 0); EXPECT_EQ(f32::silence(), 0.0f); + EXPECT_EQ(f64::silence(), 0.0); } - // The core concept is the kernels' requirement set; all three shipping - // formats satisfy it (also statically asserted in the header). + // The core concept is the kernels' requirement set; all four shipping + // formats satisfy it (also statically asserted in the header). double is + // a sample format because it is the golden model of every primitive + // (CLAUDE.md): a traits-based primitive instantiates its reference + // profile through the same substrate its embedded profiles use, and the + // cross-precision tests measure Q15/Q31/float against it. + static_assert(tap::dsp::sample_type); static_assert(tap::dsp::sample_type); static_assert(tap::dsp::sample_type); static_assert(tap::dsp::sample_type); - static_assert(!tap::dsp::sample_type); // no specialization on purpose: - // double is the golden-model/accumulator domain, not an I/O sample format. + static_assert(!tap::dsp::sample_type); // no specialization: not a format + static_assert(!tap::dsp::sample_type); + + // Every operation is constexpr: the concept's additive-identity + // requirement is checked at compile time in the header, and here for the + // value the kernels start from. + static_assert(q15::mac(q15::accum{}, std::int16_t{-32768}, std::int16_t{-16384}) == std::int64_t{32768} * 16384); + static_assert(f64::mac(f64::accum{}, 3.0, -2.0) == -6.0); } // namespace From 8d3d213d60e3e689f950b837b86356e2317793dc Mon Sep 17 00:00:00 2001 From: Timothy Place Date: Thu, 17 Sep 2026 20:05:46 +0000 Subject: [PATCH 2/2] Stage 3a review fixes: rail-aware quantize, fft_arith spec, committed table bit-pin Addresses the two hostile reviews on #22 (plan Part 13). quantize.h: the +/-1 correction now picks the largest remainder among taps that can still move in the step's direction and breaks only when none exists, so a row with a saturated tap still preserves its sum through the movable taps (a saturated tap has the largest remainder by construction and used to take, and abandon, the step). Docstring threshold corrected to 2 - 2^-15 / 2 - 2^-31 and the loop cost stated. New test NearRailRowStillPreservesItsSum ({1.99997, 0.3, 0.3} -> 42598; Q31 -> 2791728742); the rail test now covers the break path (every tap at the rail). fft_arith.h, as the Stage 3b kernel's spec: shift-before-butterfly (2 bits per radix-4 stage, 1 per radix-2, 1 before the real post-pass) and the resulting magnitude bound (sqrt(2) x full scale: 2^29.5 for widened Q15, 1.5 bits below the rail) as the condition under which the two guard bits suffice; two roundings per complex-product component, with the int64-accumulate variant named as a different, non-bit-identical design; SMMULR is not a bit-exact seam (rounds at bit 32, not 30) and would be a separately pinned Arm profile; a BFP kernel never consumes the full reported headroom (-2^20 << 11 is INT32_MIN) and 31 means "no information" ({0} and {-1} both); k_fixed_scaling_input_pre_shift (1 for Q31, 0 otherwise) with its rationale; twiddle generation source and the host-libm identity caveat. Concept fft_sample_type now requires k_guard_bits, k_coeff_frac_bits, k_coeff_one, k_fixed_scaling_input_pre_shift and noexcept on every operation. shr_round gains @pre 0 <= bits <= 31, asserted. sample_traits.h header static_asserts compare finalize(accum{}) to silence(), not literals. Decimate.FixedPointTablesAreBitPinned: row sum (2^14 / 2^30) and a byte-wise FNV-1a-64 of every Q15/Q31 coefficient table (both profiles, M = 2/3/6), measured on this substrate and identical on the pre-3a substrate, so the bit-identity evidence lives in the repo. Docs: CLAUDE.md sample_traits line lists double; README decimate paragraph follows the house rule (double golden, float embedded); "Five headers" plus the fft_arith trait with a forward note to the 3b README rewrite; the Q31 row says round-half-up; decimate.h notes Q31 is pinned against double too. test_fft_arith.cpp is listed alphabetically in tests/CMakeLists.txt. Bit identity re-verified: the pre-change Q15/Q31 dump matches byte for byte. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019ZPTzNxo5Fe4EtpXXKf7Sy --- CLAUDE.md | 2 +- README.md | 25 ++++-- include/tap/dsp/decimate.h | 6 +- include/tap/dsp/fft/fft_arith.h | 143 ++++++++++++++++++++++++-------- include/tap/dsp/quantize.h | 42 ++++++---- include/tap/dsp/sample_traits.h | 10 ++- tests/test_decimate.cpp | 70 +++++++++++++++- tests/test_quantize.cpp | 57 ++++++++++--- 8 files changed, 273 insertions(+), 82 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index e06ce30..5a97065 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -12,7 +12,7 @@ the FFT came from). Seven primitives today: the real FFT (`fft.h`), the YIN pitc the log-mel/PCEN feature extractor (`log_mel.h`) and the fixed-ratio decimators to 16 kHz (`decimate.h`) — and the dense/GRU inference kernels (`nn.h`) — plus the FIR substrate carried from SampleRateTap for the two rate converters (SampleRateTap, RatioTap): Kaiser prototype design -(`kaiser.h`), the sample-format traits (`sample_traits.h`: float/Q15/Q31), the FIR dot kernels +(`kaiser.h`), the sample-format traits (`sample_traits.h`: double/float/Q15/Q31), the FIR dot kernels (`fir_kernels.h`), row-sum-preserving quantization (`quantize.h`), and the measurement instruments (`analysis/`). See `README.md` for each asset's contract summary. diff --git a/README.md b/README.md index 44aa32a..ec1de08 100644 --- a/README.md +++ b/README.md @@ -170,8 +170,8 @@ arithmetic on the float path: the RP2350's Cortex-M33 has no FP64. front end: `basic_decimator` for M = 2, 3, 6 (32 / 48 / 96 kHz in), in RatioTap's pattern — ratio as a type, Kaiser-windowed sinc from `kaiser.h` with the cutoff at the output Nyquist and DC gain exactly 1, -`fir_kernels.h`'s `dot_row` over the `sample_traits.h` formats (float golden, -Q15 / Q31 with row-sum-preserving quantization). It is deliberately not +`fir_kernels.h`'s `dot_row` over the `sample_traits.h` formats (double golden, +float embedded, Q15 / Q31 with row-sum-preserving quantization). It is deliberately not RatioTap, whose charter is 44.1 ↔ 48 only; a 44.1 kHz host composes RatioTap's 44.1 → 48 in front of the by-3 stage. Odd tap counts, integer group delay `(taps - 1) / 2`, one output as input `k*M` arrives. @@ -228,7 +228,10 @@ suppressor's cross-precision pin must be unchanged by the promotion. Five headers carried from **SampleRateTap** (where they design and run the ASRC's polyphase datapath) and promoted here so **RatioTap**'s fixed-ratio -44.1↔48 converter — and any future FIR consumer — shares one implementation. +44.1↔48 converter — and any future FIR consumer — shares one implementation, +plus the FFT's butterfly arithmetic trait (`fft/fft_arith.h`), which is built +over the same sample formats and documented here until the Stage 3b README +rewrite moves it to the FFT section's profiles table. The performance-sensitive pieces are regression-gated in SampleRateTap's instruction-count CI (Cortex-M33/M55, Hexagon, ±3%); treat measured claims in the header comments as contracts. @@ -259,7 +262,7 @@ pins measure float/Q15/Q31 against it. | `double` | double | double | identity (the golden model) | | `float` | float | double | plain cast | | `std::int16_t` | Q1.14 | int64, exact | single Q29→Q15 round-half-up, saturating | -| `std::int32_t` | Q1.30 | int64, products pre-shifted to Q45 | Q45→Q31, saturating | +| `std::int32_t` | Q1.30 | int64, products pre-shifted to Q45 | single Q45→Q31 round-half-up, saturating | The Q ladder is spelled as named constants on each fixed-point specialization (`k_sample_frac_bits`, `k_coeff_frac_bits`, @@ -290,8 +293,11 @@ specializations and refine the `tap::dsp::sample_type` concept. ### `tap/dsp/fft/fft_arith.h` — butterfly arithmetic for the FFT profiles The sibling trait the fixed-point real FFT is written against (its docstrings -are that kernel's specification; `tests/test_fft_arith.cpp` pins every -number). One int32 kernel serves both fixed profiles: `fft_arith` +are that kernel's specification — shift-before-butterfly, the magnitude bound, +the rounding count per complex product, the BFP headroom rule, the Q31 input +pre-shift, the twiddle generator; `tests/test_fft_arith.cpp` pins every +number). Documented here until the Stage 3b README rewrite; it moves to the +FFT section's profiles table then. One int32 kernel serves both fixed profiles: `fft_arith` carries `mul_coeff` (int32 × Q1.30 → int64, `>> 30` with one round-half-up, saturating), saturating `add` / `sub`, `shr_round` (round-half-up), and `headroom_bits` (the block's shared redundant sign bits, via @@ -320,9 +326,10 @@ targets prefer which layout. coefficient format while preserving the row's DC sum *exactly* (largest-remainder distribution of the rounding residual — "the coefficients of every phase must add to one", R. Bristow-Johnson, music-dsp). Selected by -the trait's `k_is_fixed_point`: plain conversion for double and float. Each -correction step saturates, so a tap at the format's rail is never wrapped. -Design-time code. +the trait's `k_is_fixed_point`: plain conversion for double and float. A tap +at the format's rail is never wrapped: each step goes to the largest-remainder +tap that can still move in that direction, so the sum is preserved whenever +one exists. Design-time code. ### `tap/dsp/analysis/` — measurement instruments diff --git a/include/tap/dsp/decimate.h b/include/tap/dsp/decimate.h index 5bb2025..79a86f4 100644 --- a/include/tap/dsp/decimate.h +++ b/include/tap/dsp/decimate.h @@ -41,8 +41,10 @@ // float (double accumulation, the embedded profile, pinned // sample-for-sample against a committed numpy reference), Q15 and Q31 // through sample_traits.h with row-sum-preserving quantization so DC -// gain stays exactly 1; the Q15 profile is pinned against double as a -// measured number. The decimate_by_* aliases stay float: they name the +// gain stays exactly 1; the Q15 and Q31 profiles are pinned against +// double as measured numbers, and their coefficient tables are bit-pinned +// (row sum and FNV-1a-64 per ratio and profile). The decimate_by_* +// aliases stay float: they name the // 16 kHz front end's deployed profile. Mono: the consumer is // single-channel by charter. // diff --git a/include/tap/dsp/fft/fft_arith.h b/include/tap/dsp/fft/fft_arith.h index b91af46..bc95af1 100644 --- a/include/tap/dsp/fft/fft_arith.h +++ b/include/tap/dsp/fft/fft_arith.h @@ -27,12 +27,63 @@ // twiddles need no special case; every twiddle is generated in double and // rounded once by make_coeff (round_sat: half away from zero, saturating). // - The multiply is int32 x Q1.30 -> int64, then >> 30 with ONE round-half-up. -// That is the single documented rounding point per product; Armv7E-M / -// Armv8-M's SMMULR (32x32 high half with rounding) is the designated seam -// behind the same contract, as SMLALD is for the FIR kernels. -// - Growth is handled by the kernel with shr_round (fixed scaling) or -// headroom_bits + shr_round (block floating point); the trait only supplies -// the operations, never a policy. +// That is the single documented rounding point per real product. A +// complex twiddle multiply is therefore TWO roundings per output +// component: y_r = sub(mul_coeff(x_r, w_r), mul_coeff(x_i, w_i)), +// y_i = add(mul_coeff(x_r, w_i), mul_coeff(x_i, w_r)). An +// accumulate-both-products-in-int64-then-shift-once variant is a +// different design with different (lower) noise and is NOT bit-identical +// to this one; the kernel and its Welch-model pins are written against +// the two-rounding form, and the trait deliberately offers no fused op. +// - Arm's SMMULR/SMMLAR (32x32, high word, rounded) is NOT a bit-exact seam +// for mul_coeff: it rounds at bit 32 of the product where mul_coeff rounds +// at bit 30, and no arrangement reproduces one from the other (shifting +// the twiddle up two bits makes 1.0 unrepresentable, which is why Q1.30 +// was chosen; shr_round(mul_coeff(x, w), 2) is a double rounding that +// differs from (x*w + 2^31) >> 32 on about an eighth of all products). +// An SMMULR build is therefore a separately pinned Arm profile with its +// own noise numbers, decided only when an on-target measurement justifies +// it — unlike SMLALD for the FIR kernels, which is exact. No second +// multiply op is defined here so that the portable arithmetic stays the +// one contract every host and every QEMU leg runs bit-identically. +// - Scaling is the kernel's policy; the trait supplies the operations and +// states the arithmetic conditions under which they suffice: +// * Fixed scaling is SHIFT-BEFORE-BUTTERFLY: each stage's inputs are +// shr_round'ed by 2 bits (radix-4), 1 bit (radix-2), and 1 bit before +// the real post-pass, before the butterfly is computed. Under that +// ordering the complex magnitude never grows (four inputs each scaled +// by 1/4, unit-magnitude twiddles), so the magnitude bound is set by +// the input alone: at most sqrt(2) x full scale for a packed real +// pair (both components at full scale), plus at most half an LSB per +// rounding. Shift-AFTER-butterfly is not covered: four rotated +// full-scale Q2.29 inputs at a 45 degree twiddle sum to 2^31 before +// the shift and saturate. +// * Q15 (Q2.29 after widen, full scale 2^29): the bound is 2^29.5, 1.5 +// bits below the int32 rail; the two guard bits suffice with no +// input-level contract and k_fixed_scaling_input_pre_shift is 0. +// * Q31 (Q0.31, full scale 2^31 - 1): the bound 2^31.5 exceeds the rail, +// so fixed scaling takes k_fixed_scaling_input_pre_shift = 1 bit +// (shr_round, one rounding, -6 dB, one bit of 31) before the first +// stage, giving a bound of 2^30.5, half a bit below the rail. +// * Block floating point takes no input pre-shift. Before each stage the +// kernel reads headroom_bits over the block and right-shifts by +// max(0, growth bits for the stage - headroom), accumulating the +// exponent; the trait has no left shift, so BFP never normalizes a +// quiet block upward, and a shift of 0 with headroom to spare simply +// keeps the extra precision. The kernel must never let a stage +// consume the full reported headroom: a block value of -2^20 reports +// 11, and -2^20 << 11 is INT32_MIN, at which sub(0, x) saturates and +// mul_coeff(x, -1.0) is the one saturating product. Leave at least +// one bit unused. headroom_bits is 31 for an all-zero block AND for +// {-1}: 31 means "no information", not "silence". +// - Twiddles are generated at construction in double, w_k = cos/sin of +// 2*pi*k/N through std::cos/std::sin, then rounded once by make_coeff. +// Host libm last-bit differences (glibc, newlib, UCRT, Apple) can move a +// double that lies within 2^-31 of a Q1.30 rounding boundary onto the +// other side, so fixed-point outputs are host-identical only if the table +// is; the 3b battery pins the table's checksum for each certified N so a +// libm difference is detected rather than silently absorbed. (The +// quantization bound |w_q - w| <= 0.5 LSB holds on every host.) // // Every operation is constexpr and noexcept. The int32 operations preserve // the data's Q format whatever it is (Q0.31 for the Q31 profile, Q2.29 for @@ -42,6 +93,7 @@ #pragma once #include +#include #include #include #include @@ -68,10 +120,11 @@ namespace tap::dsp { using coeff = F; ///< twiddle type using work = fft_arith; - static constexpr bool k_is_fixed_point = false; - static constexpr int k_guard_bits = 0; - static constexpr int k_coeff_frac_bits = 0; - static constexpr F k_coeff_one = F{1}; + static constexpr bool k_is_fixed_point = false; + static constexpr int k_guard_bits = 0; + static constexpr int k_coeff_frac_bits = 0; + static constexpr int k_fixed_scaling_input_pre_shift = 0; ///< floating data has an exponent + static constexpr F k_coeff_one = F{1}; static constexpr coeff make_coeff(double w) noexcept { return static_cast(w); } @@ -84,7 +137,9 @@ namespace tap::dsp { /// x * 2^-bits, exact (a power-of-two scale never rounds in binary /// floating point short of underflow). + /// @pre 0 <= bits <= 31 static constexpr F shr_round(F x, int bits) noexcept { + assert(bits >= 0 && bits <= 31); F s = F{1}; for (int i = 0; i < bits; ++i) { s *= F{0.5}; @@ -115,14 +170,21 @@ namespace tap::dsp { /// `AddSubSaturateAtTheRails`. /// - shr_round(x, bits) -> sample : x * 2^-bits, round-half-up in the /// discarded bits (add 2^(bits-1), arithmetic shift, in int64 so the - /// rounding add cannot overflow). bits == 0 is the identity; the - /// result always fits for 0 <= bits <= 31. `ShrRoundRoundsHalfUp`. + /// rounding add cannot overflow). @pre 0 <= bits <= 31 (asserted); + /// bits == 0 is the identity and the result always fits. + /// `ShrRoundRoundsHalfUp`. /// - headroom_bits(x, n) -> int : the number of redundant sign bits /// shared by the whole block — the largest s such that every x[i] << s /// still fits in int32 (Arm's CLS, computed as countl_zero of the OR of /// the one's-complement magnitudes x ^ (x >> 31), minus the sign bit). - /// 31 for an all-zero or empty block, 30 for {1}, 0 for INT32_MAX or - /// INT32_MIN. `HeadroomBitsOnZeroOneAndFullScale`. + /// 31 for an all-zero or empty block and for {-1} (31 is "no + /// information", not a silence signal), 30 for {1}, 0 for INT32_MAX or + /// INT32_MIN. A BFP kernel never consumes the full amount (file + /// header). `HeadroomBitsOnZeroOneAndFullScale`. + /// - k_fixed_scaling_input_pre_shift = 1 : under fixed scaling the Q31 + /// profile shr_rounds its input by one bit before the first stage + /// (file header: the bound 2^31.5 of a full-scale packed pair exceeds + /// the rail; with the pre-shift it is 2^30.5). BFP does not use it. /// - widen / narrow : identities (wide == sample). template <> struct fft_arith { @@ -135,6 +197,9 @@ namespace tap::dsp { static constexpr bool k_is_fixed_point = true; static constexpr int k_guard_bits = 0; ///< Q0.31 data has no spare width + /// One bit (-6 dB) before the first stage under fixed scaling: the + /// price of no guard bits. Not applied under block floating point. + static constexpr int k_fixed_scaling_input_pre_shift = 1; /// Q1.30: the twiddle's fraction bits, shared with the FIR Q31 coefficient. static constexpr int k_coeff_frac_bits = sample_traits::k_coeff_frac_bits; static constexpr coeff k_coeff_one = coeff{1} << k_coeff_frac_bits; ///< 1.0 in Q1.30 @@ -160,7 +225,8 @@ namespace tap::dsp { } static constexpr sample shr_round(sample x, int bits) noexcept { - if (bits <= 0) { + assert(bits >= 0 && bits <= 31); + if (bits == 0) { return x; } return static_cast((static_cast(x) + (std::int64_t{1} << (bits - 1))) >> bits); @@ -186,8 +252,13 @@ namespace tap::dsp { /// /// - widen(x) -> wide : Q0.15 -> Q2.29, x << 14. Exact. The two guard /// bits above full scale are what let worst-case radix-4 growth run - /// under fixed scaling without an input-level contract (audit Part 7). - /// INT16_MIN -> -2^29. `WidenPlacesTwoGuardBits`. + /// under fixed scaling without an input-level contract — ON THE + /// CONDITION that the kernel shifts before each butterfly (2 bits per + /// radix-4 stage, 1 per radix-2, 1 before the real post-pass), which + /// bounds the complex magnitude by sqrt(2) x 2^29 = 2^29.5, 1.5 bits + /// below the int32 rail (file header). No input pre-shift: + /// k_fixed_scaling_input_pre_shift = 0. INT16_MIN -> -2^29. + /// `WidenPlacesTwoGuardBits`. /// - narrow(y) -> sample : Q2.29 -> Q0.15, (y + 2^13) >> 14 in int64, /// round-half-up, saturating to [INT16_MIN, INT16_MAX] — the /// substrate's finalize rule. `NarrowRoundsHalfUpAndSaturates`. @@ -201,11 +272,12 @@ namespace tap::dsp { using coeff = sample_traits::coeff; ///< Q1.30 twiddle, same as Q31 using work = fft_arith; ///< the arithmetic the kernel runs on `wide` - static constexpr bool k_is_fixed_point = true; - static constexpr int k_guard_bits = 2; ///< spare bits above full scale after widen() - static constexpr int k_widen_shift = 32 - 16 - k_guard_bits; ///< 14: Q0.15 -> Q2.29 - static constexpr int k_coeff_frac_bits = work::k_coeff_frac_bits; - static constexpr coeff k_coeff_one = work::k_coeff_one; + static constexpr bool k_is_fixed_point = true; + static constexpr int k_guard_bits = 2; ///< spare bits above full scale after widen() + static constexpr int k_fixed_scaling_input_pre_shift = 0; ///< the guard bits cover the growth bound + static constexpr int k_widen_shift = 32 - 16 - k_guard_bits; ///< 14: Q0.15 -> Q2.29 + static constexpr int k_coeff_frac_bits = work::k_coeff_frac_bits; + static constexpr coeff k_coeff_one = work::k_coeff_one; static_assert(k_widen_shift == 14, "Q15 I/O: two guard bits means a 14-bit widen"); static constexpr coeff make_coeff(double w) noexcept { return work::make_coeff(w); } @@ -219,23 +291,28 @@ namespace tap::dsp { }; // ANCHOR_END: fa_q15 - /// Satisfied by the four FFT profiles' sample types: the I/O-width pair, - /// the twiddle type and its generator, and a `work` trait that carries the - /// butterfly operations over `wide`. + /// Satisfied by the four FFT profiles' sample types: everything a + /// profile-generic kernel reads — the I/O-width pair, the twiddle type, + /// its generator and its constants, the scaling constants, and a `work` + /// trait that carries the butterfly operations over `wide`, all noexcept. template concept fft_sample_type = requires(T x, typename fft_arith::wide y, typename fft_arith::coeff w, double d) { typename fft_arith::wide; typename fft_arith::coeff; typename fft_arith::work; requires std::is_same_v::k_is_fixed_point), const bool>; - { fft_arith::make_coeff(d) } -> std::same_as::coeff>; - { fft_arith::widen(x) } -> std::same_as::wide>; - { fft_arith::narrow(y) } -> std::same_as; - { fft_arith::work::mul_coeff(y, w) } -> std::same_as::wide>; - { fft_arith::work::add(y, y) } -> std::same_as::wide>; - { fft_arith::work::sub(y, y) } -> std::same_as::wide>; - { fft_arith::work::shr_round(y, 1) } -> std::same_as::wide>; - { fft_arith::work::headroom_bits(&y, std::size_t{1}) } -> std::same_as; + requires std::is_same_v::k_guard_bits), const int>; + requires std::is_same_v::k_coeff_frac_bits), const int>; + requires std::is_same_v::k_fixed_scaling_input_pre_shift), const int>; + requires std::is_same_v::k_coeff_one), const typename fft_arith::coeff>; + { fft_arith::make_coeff(d) } noexcept -> std::same_as::coeff>; + { fft_arith::widen(x) } noexcept -> std::same_as::wide>; + { fft_arith::narrow(y) } noexcept -> std::same_as; + { fft_arith::work::mul_coeff(y, w) } noexcept -> std::same_as::wide>; + { fft_arith::work::add(y, y) } noexcept -> std::same_as::wide>; + { fft_arith::work::sub(y, y) } noexcept -> std::same_as::wide>; + { fft_arith::work::shr_round(y, 1) } noexcept -> std::same_as::wide>; + { fft_arith::work::headroom_bits(&y, std::size_t{1}) } noexcept -> std::same_as; }; static_assert(fft_sample_type); diff --git a/include/tap/dsp/quantize.h b/include/tap/dsp/quantize.h index e7d9a36..cb28121 100644 --- a/include/tap/dsp/quantize.h +++ b/include/tap/dsp/quantize.h @@ -11,6 +11,7 @@ #include #include +#include #include #include @@ -36,14 +37,20 @@ namespace tap::dsp { /// /// The algorithm is selected by the trait's k_is_fixed_point property: /// for the floating formats (double, float) this is a plain conversion - /// (no correction runs; k_coeff_scale is 1). Each +/-1 correction step is - /// saturating: a tap already at the coefficient format's rail is never - /// wrapped. Such a tap cannot absorb a step, so the correction stops - /// there and the row keeps its saturated sum — a row whose design exceeds - /// the format (|c| >= 2.0 in Q1.14 / Q1.30) had no representable DC gain - /// to preserve, and every tap of it is still the nearest representable - /// value. Design-time code: allocates a scratch vector for integer - /// formats, so keep it off the audio path. + /// (no correction runs; k_coeff_scale is 1). The +/-1 correction steps + /// never wrap: each step is taken by the largest-remainder tap among + /// those that can still move in the step's direction (not already at + /// that rail), so the sum is preserved whenever at least one such tap + /// exists, and only a row whose every tap sits at the needed rail keeps + /// its saturated sum. A tap saturates in make_coeff from + /// |c| >= 2 - 2^-15 (Q1.14) / 2 - 2^-31 (Q1.30); such a tap has the + /// row's largest remainder by construction, so without this rule it + /// would take, and wrap on, every step. Cost: one pass over the row per + /// residual LSB. A row inside the format has |residual| <= taps / 2 + /// (rounding only); a tap saturated by k LSB adds k, so a grossly + /// out-of-format row (a design bug) is slow as well as wrong. + /// Design-time code: allocates a scratch vector for integer formats, so + /// keep it off the audio path. /// /// @pre dst.size() == src.size() template @@ -71,19 +78,20 @@ namespace tap::dsp { } std::int64_t residual = static_cast(std::llround(exact_sum)) - quant_sum; while (residual != 0) { - const double sgn = residual > 0 ? 1.0 : -1.0; - std::size_t best = 0; - for (std::size_t u = 1; u < n; ++u) { - if (sgn * remainder[u] > sgn * remainder[best]) { + const double sgn = residual > 0 ? 1.0 : -1.0; + const coeff rail = residual > 0 ? std::numeric_limits::max() : std::numeric_limits::min(); + // Largest remainder among the taps that can still move this way. + std::size_t best = n; + for (std::size_t u = 0; u < n; ++u) { + if (dst[u] != rail && (best == n || sgn * remainder[u] > sgn * remainder[best])) { best = u; } } - const std::int64_t step = residual > 0 ? 1 : -1; - const coeff stepped = detail::clamp_sat(static_cast(dst[best]) + step); - if (stepped == dst[best]) { - break; // at the rail: saturate, never wrap, and the residual is unabsorbable + if (best == n) { + break; // every tap is at this rail: the residual is unabsorbable without a wrap } - dst[best] = stepped; + const std::int64_t step = residual > 0 ? 1 : -1; + dst[best] = static_cast(dst[best] + step); // off the rail: cannot overflow remainder[best] -= sgn; residual -= step; } diff --git a/include/tap/dsp/sample_traits.h b/include/tap/dsp/sample_traits.h index c4d4bb6..e124e1c 100644 --- a/include/tap/dsp/sample_traits.h +++ b/include/tap/dsp/sample_traits.h @@ -298,15 +298,17 @@ namespace tap::dsp { // single mac from accum{} yields the bare product in the accumulator // domain (Q29 for Q15, Q45 for Q31, the exact double product for the // floating formats). - static_assert(sample_traits::finalize(sample_traits::accum{}) == 0.0); + static_assert(sample_traits::finalize(sample_traits::accum{}) == sample_traits::silence()); static_assert(sample_traits::mac(sample_traits::accum{}, 0.5, 0.5) == 0.25); - static_assert(sample_traits::finalize(sample_traits::accum{}) == 0.0f); + static_assert(sample_traits::finalize(sample_traits::accum{}) == sample_traits::silence()); static_assert(sample_traits::mac(sample_traits::accum{}, 0.5f, 0.5f) == 0.25); - static_assert(sample_traits::finalize(sample_traits::accum{}) == 0); + static_assert(sample_traits::finalize(sample_traits::accum{}) + == sample_traits::silence()); static_assert(sample_traits::mac(sample_traits::accum{}, std::int16_t{32767}, std::int16_t{16384}) == std::int64_t{32767} * 16384); - static_assert(sample_traits::finalize(sample_traits::accum{}) == 0); + static_assert(sample_traits::finalize(sample_traits::accum{}) + == sample_traits::silence()); static_assert(sample_traits::mac(sample_traits::accum{}, std::int32_t{1} << 30, std::int32_t{1} << 30) == std::int64_t{1} << 44); diff --git a/tests/test_decimate.cpp b/tests/test_decimate.cpp index 03860e6..d617e44 100644 --- a/tests/test_decimate.cpp +++ b/tests/test_decimate.cpp @@ -7,15 +7,18 @@ // the reference's own float32 rounding), chunking invariance, exact output // counting, unity DC gain (exact in Q15 through row-sum-preserving // quantization), the stated passband/stopband numbers measured from the -// shipped coefficients, the integer group delay, and the Q15 / Q31 profiles -// tracking the double golden model within each format's own floor. The -// float/double contract tests are typed over both profiles. +// shipped coefficients, the integer group delay, the Q15 / Q31 profiles +// tracking the double golden model within each format's own floor, and a +// byte-wise bit pin of every fixed-point coefficient table. The float/double +// contract tests are typed over both profiles. #include #include #include #include +#include #include +#include #include #include @@ -77,7 +80,11 @@ namespace { // The double golden model against the numpy float64 reference, which is // committed rounded to float32: the disagreement is the reference's own // rounding. Measured 2026-09 at 5.75e-8 over every ratio and profile - // (half a float ulp at the signal's level); pinned at 2x. + // (half a float ulp at the signal's level); pinned at 2x. This bounds + // the double path to one float ulp of numpy, no tighter: a double + // reference (the generator already emits k_mel_log as double) would pin + // it near 1e-12 — a recorded follow-up for the reference-script owner, + // since the numpy pin was frozen for this change. TEST(Decimate, DoubleMatchesNumpyReferenceToItsFloatRounding) { using namespace frontend_ref; const auto check = [](auto ratio, const decimate_profile& p, const float* ref, std::size_t n_ref) { @@ -252,6 +259,61 @@ namespace { << "Q31 vs double, in full-scale units"; } + // FNV-1a-64 over the coefficient bytes in memory order (little-endian on + // every CI host and target; the pin is byte-wise so a one-LSB move in + // any tap changes it). + template + std::uint64_t fnv1a64(std::span c) { + std::uint64_t h = 0xcbf29ce484222325ULL; + for (const C v : c) { + unsigned char bytes[sizeof(C)]; + std::memcpy(bytes, &v, sizeof bytes); + for (const unsigned char b : bytes) { + h ^= b; + h *= 0x100000001b3ULL; + } + } + return h; + } + + template + void expect_table_pinned(const decimate_profile& p, std::size_t taps, std::int64_t sum, std::uint64_t fnv) { + basic_decimator dec(p); + const auto c = dec.coefficients(); + ASSERT_EQ(c.size(), taps); + std::int64_t s = 0; + for (const auto v : c) { + s += v; + } + EXPECT_EQ(s, sum) << "DC gain exactly 1 in the format's unity"; + EXPECT_EQ(fnv1a64(c), fnv) << "format " << sizeof(S) * 8 << " ratio " << M << " taps " << taps; + } + + // The committed bit pin for the fixed-point tables: every Q15 and Q31 + // coefficient table, both profiles, all three ratios — row sum (the + // "DC gain exactly 1" contract, 2^14 / 2^30) and a byte-wise FNV-1a-64. + // Measured 2026-09 on the Stage 3a substrate and identical on the + // substrate before it (bare Q-ladder literals): this is the evidence + // that the named constants changed no bit, kept in-tree so a later + // kernel or backend change that moves a coefficient LSB inside the + // tolerance pins is caught here. + TEST(Decimate, FixedPointTablesAreBitPinned) { + const auto eco = decimate_profile::economy(); + const auto tra = decimate_profile::transparent(); + expect_table_pinned(eco, 81, 16384, 0x9c824c3cc6602f9dULL); + expect_table_pinned(tra, 259, 16384, 0x57fcb157f4ee1c25ULL); + expect_table_pinned(eco, 121, 16384, 0x6b81f5469b8d019cULL); + expect_table_pinned(tra, 389, 16384, 0xbb59dcc474a57d61ULL); + expect_table_pinned(eco, 239, 16384, 0x8810d24884cf9bc5ULL); + expect_table_pinned(tra, 773, 16384, 0xa69d8469e463b919ULL); + expect_table_pinned(eco, 81, 1073741824, 0x8ba26dffc4b03946ULL); + expect_table_pinned(tra, 259, 1073741824, 0x980e87e92485f881ULL); + expect_table_pinned(eco, 121, 1073741824, 0xa8406e283f4cef8eULL); + expect_table_pinned(tra, 389, 1073741824, 0x3550a72e43996bd7ULL); + expect_table_pinned(eco, 239, 1073741824, 0xe20c46ad95d31338ULL); + expect_table_pinned(tra, 773, 1073741824, 0x4cd54a0678ad8bcaULL); + } + TEST(Decimate, Q15SaturatesInsteadOfWrapping) { // A full-scale step: the linear-phase lowpass pre-rings below zero by up // to the Gibbs undershoot (about 9 % of the step; measured 6.8 % here), diff --git a/tests/test_quantize.cpp b/tests/test_quantize.cpp index 9055921..a8f3b1d 100644 --- a/tests/test_quantize.cpp +++ b/tests/test_quantize.cpp @@ -72,28 +72,61 @@ namespace { EXPECT_EQ(q, row); } - // The +/-1 correction steps saturate. A row whose design exceeds the - // coefficient format has a tap at the rail with a huge positive remainder, - // so the largest-remainder pick lands on it: the step must leave it at - // the rail, never wrap it to the opposite sign, and the routine returns. + // The +/-1 correction never wraps a tap at the rail. A tap that saturated + // in make_coeff has the row's largest positive remainder, so a naive + // largest-remainder pick would step it: instead the step goes to the + // largest remainder among taps that can still move, and the sum is + // preserved through them. When every tap sits at the needed rail the + // correction stops and the row keeps its saturated sum. TEST(Quantize, CorrectionNeverWrapsATapAtTheRail) { - const std::vector row{3.0, 0.4, -0.25}; + // Just over the rail: 2.0001 * 16384 = 32769.6 -> 32767 (3 LSB short). + const std::vector row{2.0001, 0.4, -0.25}; std::vector q15(row.size()); quantize_row_preserving_sum(row, q15); - EXPECT_EQ(q15[0], 32767); // saturated, not wrapped - EXPECT_EQ(q15[1], 6554); // 0.4 * 16384 = 6553.6, its own rounding, untouched - EXPECT_EQ(q15[2], -4096); + EXPECT_EQ(q15[0], 32767); // saturated, not wrapped, untouched + EXPECT_EQ(std::int64_t{q15[0]} + q15[1] + q15[2], std::llround((2.0001 + 0.4 - 0.25) * 16384)); std::vector q31(row.size()); quantize_row_preserving_sum(row, q31); EXPECT_EQ(q31[0], 2147483647); - EXPECT_EQ(q31[1], 429496730); // 0.4 * 2^30 = 429496729.6 - EXPECT_EQ(q31[2], -268435456); + EXPECT_EQ(std::int64_t{q31[0]} + q31[1] + q31[2], std::llround((2.0001 + 0.4 - 0.25) * 1073741824.0)); // The negative rail too. - const std::vector neg{-3.0, -0.4}; + const std::vector neg{-2.0001, -0.4}; std::vector n15(neg.size()); quantize_row_preserving_sum(neg, n15); EXPECT_EQ(n15[0], -32768); - EXPECT_EQ(n15[1], -6554); + EXPECT_EQ(std::int64_t{n15[0]} + n15[1], std::llround((-2.0001 - 0.4) * 16384)); + // Every tap at the rail: nothing can move, the loop stops, no wrap. + const std::vector all{2.5, 3.0}; + std::vector a15(all.size()); + quantize_row_preserving_sum(all, a15); + EXPECT_EQ(a15[0], 32767); + EXPECT_EQ(a15[1], 32767); + const std::vector alln{-2.5, -3.0}; + quantize_row_preserving_sum(alln, a15); + EXPECT_EQ(a15[0], -32768); + EXPECT_EQ(a15[1], -32768); + } + + // The rule that makes the rail case honest: a saturated tap has the + // row's largest positive remainder, so a naive largest-remainder pick + // would land on it and either wrap (main before this) or give up on a + // sum that the other taps can still absorb. 1.99997 is the first value + // that saturates in Q1.14 (32767.5 rounds to the rail); the row's sum + // 42597.9 -> 42598 is representable and must be reached through taps 1 + // and 2. + TEST(Quantize, NearRailRowStillPreservesItsSum) { + const std::vector row{1.99997, 0.3, 0.3}; + std::vector q15(row.size()); + quantize_row_preserving_sum(row, q15); + EXPECT_EQ(q15[0], 32767); // at the rail, untouched + EXPECT_EQ(q15[1] + q15[2], 4915 + 4915 + 1); // the step went to a movable tap + EXPECT_EQ(std::int64_t{q15[0]} + q15[1] + q15[2], std::llround(1.99997 * 16384 + 0.6 * 16384)); + const std::vector row31{1.9999999995, 0.3, 0.3}; + std::vector q31(row31.size()); + quantize_row_preserving_sum(row31, q31); + EXPECT_EQ(q31[0], 2147483647); + EXPECT_EQ(std::int64_t{q31[0]} + q31[1] + q31[2], std::llround((1.9999999995 + 0.6) * 1073741824.0)); + EXPECT_EQ(std::int64_t{q31[0]} + q31[1] + q31[2], 2791728742); } TEST(Quantize, FloatIsPlainConversion) {