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 b438235..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. @@ -245,16 +248,30 @@ 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 | +| `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`, +`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 @@ -273,6 +290,24 @@ 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 — 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 +`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 +325,11 @@ 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. 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 d6c971e..79a86f4 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,16 @@ // 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 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. // // Construction designs the filter (runtime double, off the audio path) and // allocates; process() and reset() are noexcept and allocation-free. @@ -185,7 +191,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..bc95af1 --- /dev/null +++ b/include/tap/dsp/fft/fft_arith.h @@ -0,0 +1,328 @@ +/// @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 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 +// 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 + +#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 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); } + + 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). + /// @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}; + } + 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). @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 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 { + 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 + /// 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 + 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 { + assert(bits >= 0 && bits <= 31); + 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 — 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`. + /// - 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_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); } + + 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: 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>; + 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); + 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..cb28121 100644 --- a/include/tap/dsp/quantize.h +++ b/include/tap/dsp/quantize.h @@ -11,8 +11,8 @@ #include #include +#include #include -#include #include #include "tap/dsp/sample_traits.h" @@ -35,18 +35,31 @@ 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). 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() + /// @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]); } @@ -65,16 +78,22 @@ 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; } } - dst[best] = static_cast(dst[best] + (residual > 0 ? 1 : -1)); + if (best == n) { + break; // every tap is at this rail: the residual is unabsorbable without a wrap + } + 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 -= residual > 0 ? 1 : -1; + residual -= step; } } } diff --git a/include/tap/dsp/sample_traits.h b/include/tap/dsp/sample_traits.h index 7ad396e..e124e1c 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,88 @@ 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{}) == 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{}) == 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{}) + == 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{}) + == 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); // 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..d617e44 100644 --- a/tests/test_decimate.cpp +++ b/tests/test_decimate.cpp @@ -2,18 +2,24 @@ // 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, 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 #include @@ -26,6 +32,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 +77,54 @@ 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. 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) { + 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 +149,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 +203,115 @@ 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"; + } + + // 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) { 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..a8f3b1d 100644 --- a/tests/test_quantize.cpp +++ b/tests/test_quantize.cpp @@ -64,6 +64,71 @@ 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 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) { + // 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, 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(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{-2.0001, -0.4}; + std::vector n15(neg.size()); + quantize_row_preserving_sum(neg, n15); + EXPECT_EQ(n15[0], -32768); + 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) { 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