From 79233dd4a5b0540c4897be7577e28289d63e8cf2 Mon Sep 17 00:00:00 2001 From: Timothy Place Date: Fri, 18 Sep 2026 01:27:35 +0000 Subject: [PATCH 1/4] tests: fixed-point FFT battery scaffold, written against the Stage 3b contract Adds tests/test_fft_fixed.cpp (Part 9's fixed-point battery: exponents, exact fixed scale, saturation sweep, round trip per policy, BFP against fixed, rounding symmetry, Q15 vs Q31, the Welch noise-floor model and the twiddle-table pins) typed over {Q15, Q31} x {fixed, block_floating}, and tests/support/sample_scale so the shared generators land in the fixed profiles. Every measured pin is 0.0 (red) until the real kernel supplies it. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019ZPTzNxo5Fe4EtpXXKf7Sy --- tests/CMakeLists.txt | 1 + tests/support/signals.h | 57 +- tests/test_fft_fixed.cpp | 1212 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 1267 insertions(+), 3 deletions(-) create mode 100644 tests/test_fft_fixed.cpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index da9020f..a68cc58 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -127,6 +127,7 @@ tap_dsp_add_gtest_executable(tap_dsp_tests test_fft.cpp test_fft_arith.cpp test_fft_backend.cpp + test_fft_fixed.cpp test_fir_kernels.cpp test_kaiser.cpp test_log_mel.cpp diff --git a/tests/support/signals.h b/tests/support/signals.h index 0a6474d..e5b54b0 100644 --- a/tests/support/signals.h +++ b/tests/support/signals.h @@ -9,17 +9,67 @@ // helper Part 9 lists once a caller exists, is the Stage 6 hygiene item in // docs/audit-fft-and-code-smells.md, not this file's job). Everything here is deterministic — fixed seeds, no wall // clock, no filesystem — so it can run unchanged on the bare-metal QEMU legs. +// +// sample_scale is how a generator lands in a profile's sample type: +// the identity cast for float and double, and round-to-nearest with saturation +// into Q0.15 / Q0.31 for the fixed-point profiles (Stage 3b), so the same +// seed produces the same signal in every profile up to the profile's own +// quantisation. Four files share it (test_fft.cpp, test_fft_fixed.cpp, +// test_fft_oracle.cpp, test_fft_rt.cpp), which is the bar for a helper here. #pragma once #include +#include #include #include +#include #include +#include #include +#include "tap/dsp/sample_traits.h" + namespace tap::dsp::test { + /// How a value in the double domain (fractions of full scale, 1.0 = full + /// scale) crosses into and out of a profile's sample type. + /// - float, double : the plain cast; k_lsb is 0 and k_full_scale is 1. + /// - int16_t (Q0.15), int32_t (Q0.31): from_double rounds half away + /// from zero and saturates (the substrate's round_sat, the same + /// rounding the coefficient generators use); k_lsb is 2^-15 / 2^-31 + /// and k_full_scale the largest representable positive value, + /// 1 - k_lsb. Full-scale negative (-1.0) is INT_MIN exactly. + template + struct sample_scale; + + template + struct sample_scale { + static constexpr int k_frac_bits = 0; + static constexpr double k_lsb = 0.0; + static constexpr double k_full_scale = 1.0; + static constexpr double to_double(F v) noexcept { return static_cast(v); } + static constexpr F from_double(double v) noexcept { return static_cast(v); } + }; + + template + struct sample_scale { + static_assert(std::is_same_v || std::is_same_v, + "the fixed-point profiles are Q0.15 (int16_t) and Q0.31 (int32_t)"); + static constexpr int k_frac_bits = std::numeric_limits::digits; // 15 or 31 + static constexpr double k_scale = static_cast(std::int64_t{1} << k_frac_bits); + static constexpr double k_lsb = 1.0 / k_scale; + static constexpr double k_full_scale = 1.0 - k_lsb; + static constexpr double to_double(I v) noexcept { return static_cast(v) / k_scale; } + static constexpr I from_double(double v) noexcept { return tap::dsp::detail::round_sat(v * k_scale); } + }; + + static_assert(sample_scale::from_double(1.0) == std::numeric_limits::max()); + static_assert(sample_scale::from_double(-1.0) == std::numeric_limits::min()); + static_assert(sample_scale::from_double(0.5) == std::int32_t{1} << 30); + static_assert(sample_scale::to_double(std::int32_t{1} << 30) == 0.5); + static_assert(sample_scale::from_double(0.25) == 0.25f); + /// Marsaglia's xorshift32 (Journal of Statistical Software 8(14), 2003), /// the generator every existing copy in this repo already uses. The state /// must be non-zero; a zero seed is replaced by a fixed non-zero constant. @@ -44,13 +94,14 @@ namespace tap::dsp::test { std::uint32_t m_s; }; - /// n samples uniform in [-amplitude, amplitude), from a fixed seed. + /// n samples uniform in [-amplitude, amplitude), from a fixed seed, landed + /// in the profile through sample_scale (a plain cast for float/double). template std::vector random_signal(std::size_t n, std::uint32_t seed, double amplitude = 1.0) { xorshift32 rng(seed); std::vector x(n); for (auto& v : x) { - v = static_cast(amplitude * rng.next_unit()); + v = sample_scale::from_double(amplitude * rng.next_unit()); } return x; } @@ -70,7 +121,7 @@ namespace tap::dsp::test { const double period = static_cast(n); for (std::size_t j = 0; j < n; ++j) { const double turns = std::fmod(bin * static_cast(j), period) / period; - x[j] = static_cast(amplitude * std::cos(2.0 * std::numbers::pi * turns + phase)); + x[j] = sample_scale::from_double(amplitude * std::cos(2.0 * std::numbers::pi * turns + phase)); } return x; } diff --git a/tests/test_fft_fixed.cpp b/tests/test_fft_fixed.cpp new file mode 100644 index 0000000..41156fe --- /dev/null +++ b/tests/test_fft_fixed.cpp @@ -0,0 +1,1212 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// THE FIXED-POINT CONTRACT BATTERY for tap::dsp::basic_real_fft over the Q15 +// (std::int16_t) and Q31 (std::int32_t) profiles under both scaling policies +// (Stage 3b of docs/audit-fft-and-code-smells.md; Part 9 names this file and +// its promises). Written test-first against the API contract; the kernel +// (fft/fixed_point.h) is what has to pass it, and every number below is a +// contract point the kernel's header cites by test name. +// +// The contract, as the tests read it (fft.h and fft/fft_arith.h are the +// specification; this is the summary the assertions are written against): +// +// - forward_inplace / inverse_inplace / forward / inverse return an +// exponent e. Read the buffer as fractions of full scale (Q0.15 / Q0.31). +// With G the double golden model (basic_real_fft) on the SAME +// quantised input read as fractions: +// forward: G.forward_inplace(x) == out * 2^e +// inverse: G.inverse_inplace(a) (UNNORMALISED) == out * 2^e +// up to the kernel's rounding noise, same packing, same exp(+i) sign. +// The fixed-point inverse applies NO 2/N: the exponent carries the scale. +// - scaling::fixed: e is the constant fixed_scaling_exponent(n) = +// log2(n) + fft_arith::k_fixed_scaling_input_pre_shift in both +// directions: log2 n for Q15 (forward output exactly X/N), log2 n + 1 for +// Q31 (the one-bit input pre-shift that buys the missing guard bits). +// - scaling::block_floating: 0 <= e <= that constant, data-dependent; the +// kernel shifts only when growth requires it. +// - round trip, both policies: x == out * 2^(e_fwd + e_inv + 1 - log2 n) +// up to rounding noise (Ooura's unnormalised inverse has gain N/2). +// - saturation-free for every input under both policies; no alignment +// requirement; allocation at construction only (test_fft_rt.cpp). +// +// Measured numbers. Every tolerance and ratio here is a number measured on the +// real kernel and pinned at 2x (the log_mel pattern), never a round number; +// the table `pins` carries them all in one place with the host, compiler and +// date they were taken on. Fixed seeds, no wall clock, no filesystem, no +// : the battery runs unchanged on the four QEMU legs (Part 10), where +// N <= 2048 fixed point is cheap and the double golden model is the expensive +// part, so sizes are kept modest and TAP_DSP_PARITY_MAX_N caps the sweeps. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "support/signals.h" +#include "tap/dsp/fft.h" +#include "tap/dsp/fft/fft_arith.h" +#include "tap/dsp/fft/tables.h" + +#ifndef TAP_DSP_PARITY_MAX_N +#define TAP_DSP_PARITY_MAX_N 65536 +#endif + +namespace { + + namespace scaling = tap::dsp::scaling; + using tap::dsp::test::sample_scale; + + // ------------------------------------------------------------------------ + // The four configurations under test. + // ------------------------------------------------------------------------ + template + struct config { + using sample = Sample; + using policy = Scaling; + using fft = tap::dsp::basic_real_fft; + using arith = tap::dsp::fft_arith; + + static constexpr bool k_is_q15 = std::is_same_v; + static constexpr bool k_is_bfp = std::is_same_v; + static constexpr double k_lsb = sample_scale::k_lsb; + static constexpr double k_full = sample_scale::k_full_scale; + static constexpr Sample k_max = std::numeric_limits::max(); + static constexpr Sample k_min = std::numeric_limits::min(); + static constexpr int k_pre = arith::k_fixed_scaling_input_pre_shift; + static constexpr int k_int_frac = k_is_q15 ? 29 : 31; ///< fraction bits of the int32 work format + + static const char* name() { + if constexpr (k_is_q15) { + return k_is_bfp ? "Q15/bfp " : "Q15/fixed"; + } + else { + return k_is_bfp ? "Q31/bfp " : "Q31/fixed"; + } + } + }; + + using q15_fixed = config; + using q31_fixed = config; + using q15_bfp = config; + using q31_bfp = config; + + using all_configs = ::testing::Types; + using fixed_configs = ::testing::Types; + using bfp_configs = ::testing::Types; + + // The aliases are the documented spellings of the four configurations. + static_assert(std::is_same_v); + static_assert(std::is_same_v); + static_assert(std::is_same_v); + static_assert(std::is_same_v); + // The default policy is fixed scaling. + static_assert(std::is_same_v, q15_fixed::fft>); + + // ------------------------------------------------------------------------ + // THE PINS. Measured on the real kernel and pinned at 2x; see each test + // for what the number bounds. Host, compiler and date beside each block. + // A pin of 0.0 means "not yet measured" and fails the test on purpose. + // ------------------------------------------------------------------------ + struct pin_table { + double saturation_max_lsb; ///< SaturationFreeWorstCaseDoesNotWrap: max |out - G/2^e| in output LSB + double round_trip_k; ///< RoundTripReconstructsInputPerPolicy: max error in reconstructed LSB + double bfp_vs_fixed_lsb; ///< BfpMatchesFixedAfterShift: max |shifted bfp - fixed| in output LSB + double negation_sum_max_lsb; ///< RoundingBiasOnNegatedInputIsBounded: max |F(x) + F(-x)| in LSB + double negation_bias_lsb; ///< ... and the mean of F(x) + F(-x) in LSB + double noise_ratio_noise; ///< NoiseFloorTracksWelchModel: max measured/model, white noise + double noise_ratio_tone; ///< ... on-bin tone + double q15_vs_q31_lsb; ///< Q15AndQ31AgreeToTheQ15Floor: max |v15 - v31| in Q15 output LSB + }; + + // Not yet measured: the kernel branch has not landed. Every pin is 0.0 + // so the battery is red until the numbers are taken from the real kernel. + constexpr pin_table k_pins_q15_fixed{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + constexpr pin_table k_pins_q31_fixed{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + constexpr pin_table k_pins_q15_bfp{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + constexpr pin_table k_pins_q31_bfp{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + + template + constexpr const pin_table& pins() { + if constexpr (std::is_same_v) { + return k_pins_q15_fixed; + } + else if constexpr (std::is_same_v) { + return k_pins_q31_fixed; + } + else if constexpr (std::is_same_v) { + return k_pins_q15_bfp; + } + else { + return k_pins_q31_bfp; + } + } + + // ------------------------------------------------------------------------ + // Small helpers. + // ------------------------------------------------------------------------ + constexpr std::size_t k_max_n = std::min(65536, TAP_DSP_PARITY_MAX_N); + + constexpr int log2_size(std::size_t n) { + int l = 0; + while ((std::size_t{1} << l) < n) { + ++l; + } + return l; + } + static_assert(log2_size(4) == 2 && log2_size(65536) == 16); + + double pow2(int e) { + return std::ldexp(1.0, e); + } + + double db(double power) { + return 10.0 * std::log10(std::max(power, 1e-300)); + } + + template + std::vector fractions(const std::vector& x) { + std::vector d(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) { + d[i] = sample_scale::to_double(x[i]); + } + return d; + } + + template + std::vector quantise(const std::vector& x) { + std::vector s(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) { + s[i] = sample_scale::from_double(x[i]); + } + return s; + } + + /// The double golden model's unnormalised forward, on fractions. + std::vector golden_forward(std::vector x) { + tap::dsp::basic_real_fft g(x.size()); + g.forward_inplace(x.data()); + return x; + } + + /// The double golden model's UNNORMALISED inverse, on fractions. + std::vector golden_inverse(std::vector a) { + tap::dsp::basic_real_fft g(a.size()); + g.inverse_inplace(a.data()); + return a; + } + + template + struct transform_result { + std::vector out; + int exponent; + }; + + template + transform_result run_forward(const std::vector& x) { + typename Cfg::fft fft(x.size()); + transform_result r{x, 0}; + r.exponent = fft.forward_inplace(r.out.data()); + return r; + } + + template + transform_result run_inverse(const std::vector& a) { + typename Cfg::fft fft(a.size()); + transform_result r{a, 0}; + r.exponent = fft.inverse_inplace(r.out.data()); + return r; + } + + /// |out - golden / 2^e| over the whole buffer, in output LSBs. + struct deviation { + double max_lsb; + double rms_lsb; + std::size_t argmax; + }; + + template + deviation deviation_from_golden(const std::vector& out, const std::vector& golden, + int e) { + deviation d{0.0, 0.0, 0}; + const double scale = pow2(-e) / Cfg::k_lsb; // golden -> output LSB units + for (std::size_t i = 0; i < out.size(); ++i) { + const double err = static_cast(out[i]) - golden[i] * scale; + if (std::fabs(err) > d.max_lsb) { + d.max_lsb = std::fabs(err); + d.argmax = i; + } + d.rms_lsb += err * err; + } + d.rms_lsb = std::sqrt(d.rms_lsb / static_cast(out.size())); + return d; + } + + /// Named test vectors. + template + struct pattern { + const char* name; + std::vector x; + }; + + /// Full-scale adversarial time-domain patterns (fft_arith.h, "Scaling"): + /// every one has |x| at a rail in every sample, so it drives the + /// magnitude bound; the rotated packed pairs and the square-wave + /// exponentials are the ones that meet a 45 degree twiddle (Part 6, N6). + template + std::vector> adversarial_time_patterns(std::size_t n) { + constexpr Sample hi = std::numeric_limits::max(); + constexpr Sample lo = std::numeric_limits::min(); + + std::vector> p; + p.push_back({"dc +full", std::vector(n, hi)}); + p.push_back({"dc INT_MIN", std::vector(n, lo)}); + { + std::vector x(n); + for (std::size_t j = 0; j < n; ++j) { + x[j] = (j % 2 == 0) ? hi : lo; + } + p.push_back({"nyquist +-", x}); + for (std::size_t j = 0; j < n; ++j) { + x[j] = (j % 2 == 0) ? lo : hi; + } + p.push_back({"nyquist -+", x}); + } + { + std::vector x(n, Sample{0}); + x[0] = hi; + p.push_back({"impulse +full", x}); + x[0] = lo; + p.push_back({"impulse INT_MIN", x}); + } + // Packed pairs z_j = x[2j] + i x[2j+1] at a rail on BOTH components + // (|z| = sqrt(2) full scale), rotating by +90 / -90 degrees per pair: + // a complex exponential at bin M/4 of the complex sequence whose whole + // energy lands in one bin, coherent through every butterfly. + { + const std::array, 4> cw{{{hi, hi}, {lo, hi}, {lo, lo}, {hi, lo}}}; + const std::array, 4> ccw{{{hi, hi}, {hi, lo}, {lo, lo}, {lo, hi}}}; + const std::array, 4> cw_neg{{{lo, lo}, {hi, lo}, {hi, hi}, {lo, hi}}}; + std::vector x(n); + for (std::size_t j = 0; j < n / 2; ++j) { + x[2 * j] = cw[j % 4][0]; + x[2 * j + 1] = cw[j % 4][1]; + } + p.push_back({"pair rotation +90", x}); + for (std::size_t j = 0; j < n / 2; ++j) { + x[2 * j] = ccw[j % 4][0]; + x[2 * j + 1] = ccw[j % 4][1]; + } + p.push_back({"pair rotation -90", x}); + for (std::size_t j = 0; j < n / 2; ++j) { + x[2 * j] = cw_neg[j % 4][0]; + x[2 * j + 1] = cw_neg[j % 4][1]; + } + p.push_back({"pair rotation +90 from INT_MIN", x}); + } + // Square-wave complex exponentials: z_j = sgn cos(theta_j + phi) + + // i sgn sin(theta_j + phi), theta_j = 2 pi j k / M. Every component at + // a rail; the fundamental sits at complex bin k where the stage + // twiddles are 45 degrees (k = M/8), 22.5 degrees (M/16) and 135 + // degrees (3M/8). phi keeps the samples off the zero crossings. + for (const double k_over_m : {1.0 / 8.0, 1.0 / 16.0, 3.0 / 8.0}) { + std::vector x(n); + for (std::size_t j = 0; j < n / 2; ++j) { + const double theta = + 2.0 * std::numbers::pi * k_over_m * static_cast(j) + std::numbers::pi / 8.0; + x[2 * j] = std::cos(theta) >= 0.0 ? hi : lo; + x[2 * j + 1] = std::sin(theta) >= 0.0 ? hi : lo; + } + p.push_back({"square exponential", x}); + } + // Real square waves at bins N/8 and N/4 + 1. + for (const double k_over_n : {1.0 / 8.0, 1.0 / 4.0 + 1.0 / static_cast(n)}) { + std::vector x(n); + for (std::size_t j = 0; j < n; ++j) { + const double theta = 2.0 * std::numbers::pi * k_over_n * static_cast(j) + 0.3; + x[j] = std::cos(theta) >= 0.0 ? hi : lo; + } + p.push_back({"square wave", x}); + } + // Full-scale binary noise: the broadband worst case (fixed seeds). + for (const std::uint32_t seed : {0x2545F491u, 0x9E3779B9u, 0x1D872B41u, 0xC0FFEE01u}) { + tap::dsp::test::xorshift32 rng(seed); + std::vector x(n); + for (std::size_t j = 0; j < n; ++j) { + x[j] = (rng.next_u32() & 0x10000u) != 0u ? hi : lo; + } + p.push_back({"binary noise", x}); + } + return p; + } + + /// Full-scale adversarial packed spectra for the inverse direction. + template + std::vector> adversarial_spectrum_patterns(std::size_t n) { + constexpr Sample hi = std::numeric_limits::max(); + constexpr Sample lo = std::numeric_limits::min(); + + // The time-domain set is a fine set of spectra too (all at a rail), + // plus the spectral shapes whose inverse concentrates: one full-scale + // complex bin (a tone at amplitude sqrt(2)), and DC + Nyquist only. + std::vector> p = adversarial_time_patterns(n); + { + std::vector a(n, Sample{0}); + const std::size_t k = std::max(1, n / 8 + 1); + a[2 * k] = hi; + a[2 * k + 1] = hi; + p.push_back({"one full complex bin", a}); + a.assign(n, Sample{0}); + a[0] = hi; + a[1] = hi; + p.push_back({"dc and nyquist full", a}); + a[0] = lo; + a[1] = lo; + p.push_back({"dc and nyquist INT_MIN", a}); + } + return p; + } + + /// Sizes for the per-pattern sweeps: the two smallest geometries (M = 2 + /// and M = 4: no radix-4 stage at all, one radix-2), an odd and an even + /// log2 M, and the two log_mel / pvoc geometries. + constexpr std::array k_sweep_sizes{4, 8, 16, 64, 512, 1024, 2048}; + + // ======================================================================== + // Exponents. + // ======================================================================== + + template + class fft_fixed_scaling_test : public ::testing::Test {}; + TYPED_TEST_SUITE(fft_fixed_scaling_test, fixed_configs); + + template + class fft_bfp_test : public ::testing::Test {}; + TYPED_TEST_SUITE(fft_bfp_test, bfp_configs); + + template + class fft_fixed_point_test : public ::testing::Test {}; + TYPED_TEST_SUITE(fft_fixed_point_test, all_configs); + + // scaling::fixed: e = log2(n) + k_fixed_scaling_input_pre_shift, a + // compile-time constant, in both directions, from every entry point. + TYPED_TEST(fft_fixed_scaling_test, FixedExponentIsTheStatedConstant) { + using cfg = TypeParam; + using fft = typename cfg::fft; + static_assert(fft::fixed_scaling_exponent(1024) == 10 + cfg::k_pre); + static_assert(fft::fixed_scaling_exponent(4) == 2 + cfg::k_pre); + static_assert(noexcept(fft::fixed_scaling_exponent(4))); + if constexpr (cfg::k_is_q15) { + static_assert(fft::fixed_scaling_exponent(512) == 9, "Q15: no pre-shift, output exactly X/N"); + } + else { + static_assert(fft::fixed_scaling_exponent(512) == 10, "Q31: the one-bit input pre-shift"); + } + for (std::size_t n = 4; n <= k_max_n; n *= 2) { + const int expected = log2_size(n) + cfg::k_pre; + EXPECT_EQ(fft::fixed_scaling_exponent(n), expected) << "n=" << n; + fft f(n); + const auto x = + tap::dsp::test::random_signal(n, 0x2545F491u ^ static_cast(n)); + auto a = x; + EXPECT_EQ(f.forward_inplace(a.data()), expected) << "forward_inplace n=" << n; + EXPECT_EQ(f.inverse_inplace(a.data()), expected) << "inverse_inplace n=" << n; + std::vector out(n); + EXPECT_EQ(f.forward(x.data(), out.data()), expected) << "forward n=" << n; + EXPECT_EQ(f.inverse(out.data(), out.data()), expected) << "inverse (aliased) n=" << n; + // And a silent block reports the same constant: the exponent is + // static, not a headroom measurement. + std::vector zeros(n, 0); + EXPECT_EQ(f.forward_inplace(zeros.data()), expected) << "forward of silence n=" << n; + } + } + + // scaling::block_floating: 0 <= e <= the fixed constant, in both + // directions, whatever the input; a louder input never needs a smaller + // exponent than silence does (0 is the floor). + TYPED_TEST(fft_bfp_test, BfpExponentIsWithinRange) { + using cfg = TypeParam; + using fft = typename cfg::fft; + using s = typename cfg::sample; + for (const std::size_t n : {std::size_t{4}, std::size_t{16}, std::size_t{512}, std::size_t{2048}}) { + const int top = fft::fixed_scaling_exponent(n); + fft f(n); + + std::vector> inputs; + inputs.push_back(std::vector(n, s{0})); + inputs.push_back(std::vector(n, s{0})); + inputs.back()[0] = s{1}; + inputs.push_back(tap::dsp::test::random_signal(n, 0x9E3779B9u, 1e-3)); + inputs.push_back(tap::dsp::test::random_signal(n, 0x9E3779B9u, 0.1)); + inputs.push_back(tap::dsp::test::random_signal(n, 0x9E3779B9u, 1.0)); + inputs.push_back(std::vector(n, cfg::k_max)); + inputs.push_back(std::vector(n, cfg::k_min)); + for (auto& p : adversarial_time_patterns(n)) { + inputs.push_back(std::move(p.x)); + } + + for (std::size_t i = 0; i < inputs.size(); ++i) { + auto a = inputs[i]; + const int ef = f.forward_inplace(a.data()); + EXPECT_GE(ef, 0) << "forward n=" << n << " input " << i; + EXPECT_LE(ef, top) << "forward n=" << n << " input " << i; + const int ei = f.inverse_inplace(a.data()); + EXPECT_GE(ei, 0) << "inverse n=" << n << " input " << i; + EXPECT_LE(ei, top) << "inverse n=" << n << " input " << i; + } + } + } + + // ======================================================================== + // Exact scale under fixed scaling. + // ======================================================================== + + // An integer-valued impulse of A output LSBs has X_k = A for every k; a + // constant of c LSBs has X_0 = N c and nothing else. Under fixed scaling + // the output is X * 2^-e with e = log2 n + pre: when that is an integer + // number of LSBs no rounding can occur anywhere in the kernel (every + // stage shift discards zero bits: the impulse's value returns to A / 4^s + // exactly and the constant's to c at every stage), so the result is + // bit-exact. A = m * 2^e and c a multiple of 2^(pre + 2) (the Q31 + // pre-shift plus one radix-4 shift of the constant path; Q15's widen + // supplies fourteen zero bits so any int16 works there). + TYPED_TEST(fft_fixed_scaling_test, FixedForwardScaleIsExactlyXOverN) { + using cfg = TypeParam; + using s = typename cfg::sample; + for (const std::size_t n : k_sweep_sizes) { + const int e = cfg::fft::fixed_scaling_exponent(n); + // Impulses: A = m * 2^e for several m, the largest that fits. + const std::int64_t unit = std::int64_t{1} << e; + const std::int64_t top_m = static_cast(cfg::k_max) / unit; + for (const std::int64_t m : {std::int64_t{1}, std::int64_t{3}, top_m, -top_m}) { + if (m == 0) { + continue; + } + std::vector x(n, s{0}); + x[0] = static_cast(m * unit); + const auto r = run_forward(x); + ASSERT_EQ(r.exponent, e); + for (std::size_t i = 0; i < n; ++i) { + ASSERT_EQ(static_cast(r.out[i]), m) + << "impulse A=" << m * unit << " n=" << n << " index " << i; + } + } + // Constants: c a multiple of 2^(pre + 2); the extreme values too. + const std::int64_t step = std::int64_t{1} << (cfg::k_pre + 2); + std::vector constants{step, -step, 5 * step, + (static_cast(cfg::k_max) / step) * step, + static_cast(cfg::k_min)}; + if constexpr (cfg::k_is_q15) { + constants.push_back(cfg::k_max); // 32767: odd, exact through the 14-bit widen + constants.push_back(1); + constants.push_back(-1); + } + for (const std::int64_t c : constants) { + const std::vector x(n, static_cast(c)); + const auto r = run_forward(x); + ASSERT_EQ(r.exponent, e); + // X_0 = N c; out_0 = N c / 2^e = c / 2^pre. + const std::int64_t expected_dc = c / (std::int64_t{1} << cfg::k_pre); + ASSERT_EQ(static_cast(r.out[0]), expected_dc) << "dc c=" << c << " n=" << n; + for (std::size_t i = 1; i < n; ++i) { + ASSERT_EQ(r.out[i], s{0}) << "dc c=" << c << " n=" << n << " index " << i; + } + } + } + } + + // The Q31 profile's fixed forward is X / (2N): the extra bit is the input + // pre-shift (fft_arith::k_fixed_scaling_input_pre_shift = 1), + // stated here in the profile's own numbers. + TEST(fft_fixed_q31, FixedForwardScaleIsExactlyXOverTwoN) { + using fft = tap::dsp::real_fft_q31; + static_assert(tap::dsp::fft_arith::k_fixed_scaling_input_pre_shift == 1); + static_assert(tap::dsp::fft_arith::k_fixed_scaling_input_pre_shift == 0); + static_assert(fft::fixed_scaling_exponent(512) == 10); + constexpr std::size_t n = 512; + fft f(n); + // A constant of 2^30 (0.5): X_0 = 512 * 2^30, out_0 = 2^29 (0.25). + std::vector dc(n, std::int32_t{1} << 30); + EXPECT_EQ(f.forward_inplace(dc.data()), 10); + EXPECT_EQ(dc[0], std::int32_t{1} << 29); + // An impulse of 2^20: out = 2^20 / 1024 = 2^10 in every slot. + std::vector imp(n, 0); + imp[0] = std::int32_t{1} << 20; + EXPECT_EQ(f.forward_inplace(imp.data()), 10); + for (std::size_t i = 0; i < n; ++i) { + ASSERT_EQ(imp[i], std::int32_t{1} << 10) << i; + } + // The Q15 profile, for contrast, is X / N: the same 0.5 constant + // comes back as 0.5. + tap::dsp::real_fft_q15 f15(n); + std::vector dc15(n, std::int16_t{1} << 14); + EXPECT_EQ(f15.forward_inplace(dc15.data()), 9); + EXPECT_EQ(dc15[0], std::int16_t{1} << 14); + } + + // The inverse under fixed scaling carries the same constant and is + // exact on the same kind of input: the unnormalised inverse of a + // DC-only spectrum a[0] = c is c/2 everywhere (fft.h's packing), so the + // output is c / 2^(e+1); of a flat spectrum (impulse response) it is + // N c / 2 at k = 0 and 0 elsewhere, so out_0 = c / 2^(pre + 1). + TYPED_TEST(fft_fixed_scaling_test, FixedInverseCarriesTheSameExponent) { + using cfg = TypeParam; + using s = typename cfg::sample; + for (const std::size_t n : k_sweep_sizes) { + const int e = cfg::fft::fixed_scaling_exponent(n); + const std::int64_t unit = std::int64_t{1} << (e + 1); + for (const std::int64_t m : + {std::int64_t{1}, std::int64_t{-3}, static_cast(cfg::k_max) / unit}) { + if (m == 0) { + continue; + } + std::vector a(n, s{0}); + a[0] = static_cast(m * unit); + const auto r = run_inverse(a); + ASSERT_EQ(r.exponent, e); + for (std::size_t i = 0; i < n; ++i) { + ASSERT_EQ(static_cast(r.out[i]), m) << "dc spectrum n=" << n << " index " << i; + } + } + // The flat spectrum: a[0] = a[1] = a[2k] = c, a[2k+1] = 0. + const std::int64_t step = std::int64_t{1} << (cfg::k_pre + 2); + for (const std::int64_t c : {step, -step, (static_cast(cfg::k_max) / step) * step}) { + std::vector a(n, s{0}); + a[0] = static_cast(c); + a[1] = static_cast(c); + for (std::size_t k = 1; k < n / 2; ++k) { + a[2 * k] = static_cast(c); + } + const auto r = run_inverse(a); + ASSERT_EQ(r.exponent, e); + // x_0 = N c / 2, out_0 = N c / 2^(e + 1) = c / 2^(pre + 1). + ASSERT_EQ(static_cast(r.out[0]), c / (std::int64_t{1} << (cfg::k_pre + 1))) + << "flat spectrum c=" << c << " n=" << n; + for (std::size_t i = 1; i < n; ++i) { + ASSERT_EQ(r.out[i], s{0}) << "flat spectrum c=" << c << " n=" << n << " index " << i; + } + } + } + } + + // ======================================================================== + // Saturation. + // ======================================================================== + + // Every adversarial full-scale pattern, both directions, both policies: + // the output tracks the golden model at the returned exponent to within + // the pinned number of LSBs (a wrap is 2^15 / 2^31 LSBs off, a clamp at + // least the amount clamped), and no output sample sits on a rail unless + // the golden model puts it within the pin of that rail. + TYPED_TEST(fft_fixed_point_test, SaturationFreeWorstCaseDoesNotWrap) { + using cfg = TypeParam; + using s = typename cfg::sample; + const auto pin = pins().saturation_max_lsb; + ASSERT_GT(pin, 0.0) << "unmeasured pin"; + double worst = 0.0; + for (const std::size_t n : k_sweep_sizes) { + for (const auto& p : adversarial_time_patterns(n)) { + const auto r = run_forward(p.x); + const auto golden = golden_forward(fractions(p.x)); + const auto d = deviation_from_golden(r.out, golden, r.exponent); + worst = std::max(worst, d.max_lsb); + EXPECT_LE(d.max_lsb, pin) << cfg::name() << " forward " << p.name << " n=" << n << " e=" << r.exponent + << " at index " << d.argmax; + for (std::size_t i = 0; i < n; ++i) { + if (r.out[i] == cfg::k_max || r.out[i] == cfg::k_min) { + const double g = std::fabs(golden[i]) * pow2(-r.exponent) / cfg::k_lsb; + EXPECT_GE(g, static_cast(cfg::k_max) - pin) + << cfg::name() << " forward " << p.name << " n=" << n << ": sample " << i + << " sits on a rail the golden model does not reach"; + } + } + } + for (const auto& p : adversarial_spectrum_patterns(n)) { + const auto r = run_inverse(p.x); + const auto golden = golden_inverse(fractions(p.x)); + const auto d = deviation_from_golden(r.out, golden, r.exponent); + worst = std::max(worst, d.max_lsb); + EXPECT_LE(d.max_lsb, pin) << cfg::name() << " inverse " << p.name << " n=" << n << " e=" << r.exponent + << " at index " << d.argmax; + for (std::size_t i = 0; i < n; ++i) { + if (r.out[i] == cfg::k_max || r.out[i] == cfg::k_min) { + const double g = std::fabs(golden[i]) * pow2(-r.exponent) / cfg::k_lsb; + EXPECT_GE(g, static_cast(cfg::k_max) - pin) + << cfg::name() << " inverse " << p.name << " n=" << n << ": sample " << i + << " sits on a rail the golden model does not reach"; + } + } + } + } + std::printf("[ measured ] %s saturation sweep: max |out - G/2^e| = %.3f LSB (pin %.3f)\n", cfg::name(), worst, + pin); + } + + // A silent block is silent out, in both directions. + TYPED_TEST(fft_fixed_point_test, SilenceIsSilence) { + using cfg = TypeParam; + using s = typename cfg::sample; + for (const std::size_t n : k_sweep_sizes) { + const std::vector zeros(n, s{0}); + const auto f = run_forward(zeros); + const auto i = run_inverse(zeros); + for (std::size_t j = 0; j < n; ++j) { + ASSERT_EQ(f.out[j], s{0}) << "forward n=" << n << " index " << j; + ASSERT_EQ(i.out[j], s{0}) << "inverse n=" << n << " index " << j; + } + EXPECT_GE(f.exponent, 0); + EXPECT_LE(f.exponent, cfg::fft::fixed_scaling_exponent(n)); + } + } + + // ======================================================================== + // Round trip. + // ======================================================================== + + // x == out * 2^(e_fwd + e_inv + 1 - log2 n) up to rounding noise. The + // error is measured in units of the reconstructed LSB, i.e. the output + // LSB times that same power of two: under fixed scaling the round trip + // discards log2(N) bits twice and the unit is coarse (Q15, N = 1024: 2^-4 + // of full scale, an honest number the header states), under block + // floating point the exponents are smaller and the unit finer. The pin + // is the largest error in that unit over the sizes and levels below. + TYPED_TEST(fft_fixed_point_test, RoundTripReconstructsInputPerPolicy) { + using cfg = TypeParam; + using s = typename cfg::sample; + const auto pin = pins().round_trip_k; + ASSERT_GT(pin, 0.0) << "unmeasured pin"; + double worst = 0.0; + for (const std::size_t n : {std::size_t{16}, std::size_t{64}, std::size_t{256}, std::size_t{1024}}) { + for (const double amplitude : {1.0, 0.25, 0.01}) { + const auto x = + tap::dsp::test::random_signal(n, 0xC0FFEE01u ^ static_cast(n), amplitude); + const auto xf = fractions(x); + const auto fwd = run_forward(x); + const auto inv = run_inverse(fwd.out); + const int shift = fwd.exponent + inv.exponent + 1 - log2_size(n); + const double unit = cfg::k_lsb * pow2(shift); + double max_err = 0.0; + for (std::size_t i = 0; i < n; ++i) { + const double back = sample_scale::to_double(inv.out[i]) * pow2(shift); + max_err = std::max(max_err, std::fabs(back - xf[i]) / unit); + } + worst = std::max(worst, max_err); + EXPECT_LE(max_err, pin) << cfg::name() << " n=" << n << " amplitude " << amplitude + << " e_fwd=" << fwd.exponent << " e_inv=" << inv.exponent << " unit=" << unit; + } + } + std::printf("[ measured ] %s round trip: max error %.3f reconstructed LSB (pin %.3f)\n", cfg::name(), worst, + pin); + } + + // forward()/inverse() are copy-then-in-place: same bits, same exponent, + // with and without aliasing; and the fixed-point inverse() applies no + // 2/N (the exponent carries the scale). + TYPED_TEST(fft_fixed_point_test, OutOfPlaceIsCopyThenInPlace) { + using cfg = TypeParam; + using s = typename cfg::sample; + for (const std::size_t n : {std::size_t{8}, std::size_t{512}}) { + typename cfg::fft fft(n); + const auto x = tap::dsp::test::random_signal(n, 0x1D872B41u, 0.9); + + auto inplace = x; + const int e_inplace = fft.forward_inplace(inplace.data()); + std::vector out(n, s{0}); + const int e_out = fft.forward(x.data(), out.data()); + auto alias = x; + const int e_alias = fft.forward(alias.data(), alias.data()); + EXPECT_EQ(e_out, e_inplace); + EXPECT_EQ(e_alias, e_inplace); + EXPECT_EQ(std::memcmp(out.data(), inplace.data(), n * sizeof(s)), 0) << "forward out-of-place n=" << n; + EXPECT_EQ(std::memcmp(alias.data(), inplace.data(), n * sizeof(s)), 0) << "forward aliased n=" << n; + + auto inv_inplace = inplace; + const int ei_inplace = fft.inverse_inplace(inv_inplace.data()); + std::vector inv_out(n, s{0}); + const int ei_out = fft.inverse(inplace.data(), inv_out.data()); + auto inv_alias = inplace; + const int ei_alias = fft.inverse(inv_alias.data(), inv_alias.data()); + EXPECT_EQ(ei_out, ei_inplace); + EXPECT_EQ(ei_alias, ei_inplace); + EXPECT_EQ(std::memcmp(inv_out.data(), inv_inplace.data(), n * sizeof(s)), 0) + << "inverse out-of-place n=" << n; + EXPECT_EQ(std::memcmp(inv_alias.data(), inv_inplace.data(), n * sizeof(s)), 0) << "inverse aliased n=" << n; + } + } + + // ======================================================================== + // Block floating point against fixed scaling. + // ======================================================================== + + // The BFP output, shifted right (round-half-up) by e_fixed - e_bfp, + // agrees with the fixed output to within the pinned number of LSBs (the + // difference is the fixed path's extra rounding noise plus one rounding + // of the shift); and on an input where BFP reports the full constant it + // shifted like fixed at every stage, so the two are bit-identical. + TYPED_TEST(fft_bfp_test, BfpMatchesFixedAfterShift) { + using cfg = TypeParam; + using s = typename cfg::sample; + using fixed = config; + const auto pin = pins().bfp_vs_fixed_lsb; + ASSERT_GT(pin, 0.0) << "unmeasured pin"; + double worst = 0.0; + for (const std::size_t n : k_sweep_sizes) { + std::vector> inputs; + inputs.push_back({"noise 0 dBFS", tap::dsp::test::random_signal(n, 0x2545F491u, 1.0)}); + inputs.push_back({"noise -20 dBFS", tap::dsp::test::random_signal(n, 0x2545F491u, 0.1)}); + inputs.push_back({"noise -40 dBFS", tap::dsp::test::random_signal(n, 0x2545F491u, 0.01)}); + inputs.push_back({"tone -6 dBFS", tap::dsp::test::tone(n, static_cast(n / 8 + 1), 0.5, 0.3)}); + inputs.push_back({"dc +full", std::vector(n, cfg::k_max)}); + for (const auto& p : inputs) { + const auto f = run_forward(p.x); + const auto b = run_forward(p.x); + ASSERT_LE(b.exponent, f.exponent) << p.name << " n=" << n; + const int shift = f.exponent - b.exponent; + for (std::size_t i = 0; i < n; ++i) { + const std::int64_t bv = static_cast(b.out[i]); + const std::int64_t shifted = + shift == 0 ? bv : ((bv + (std::int64_t{1} << (shift - 1))) >> shift); // round-half-up + const double diff = std::fabs(static_cast(shifted - static_cast(f.out[i]))); + worst = std::max(worst, diff); + EXPECT_LE(diff, pin) << cfg::name() << " " << p.name << " n=" << n << " e_fixed=" << f.exponent + << " e_bfp=" << b.exponent << " index " << i; + } + } + } + std::printf("[ measured ] %s vs fixed after shift: max %.1f LSB (pin %.1f)\n", cfg::name(), worst, pin); + } + + TYPED_TEST(fft_bfp_test, BfpAtTheFullExponentIsBitIdenticalToFixed) { + using cfg = TypeParam; + using s = typename cfg::sample; + using fixed = config; + std::size_t attained = 0; + std::size_t tried = 0; + for (const std::size_t n : k_sweep_sizes) { + for (const auto& p : adversarial_time_patterns(n)) { + ++tried; + const auto f = run_forward(p.x); + const auto b = run_forward(p.x); + if (b.exponent != f.exponent) { + continue; + } + ++attained; + EXPECT_EQ(std::memcmp(b.out.data(), f.out.data(), n * sizeof(s)), 0) + << cfg::name() << " " << p.name << " n=" << n << ": BFP reported the full exponent " << f.exponent + << " but its output is not the fixed output"; + } + } + // The contract's upper bound is attainable: some full-scale input + // makes BFP shift like fixed at every stage. + EXPECT_GT(attained, 0u) << cfg::name() << ": no full-scale pattern out of " << tried + << " reached the fixed exponent"; + std::printf("[ measured ] %s: %zu of %zu full-scale patterns reach the fixed exponent\n", cfg::name(), attained, + tried); + } + + // ======================================================================== + // Rounding symmetry. + // ======================================================================== + + // Every rounding in the kernel is round-half-up (fft_arith.h), which is + // not odd-symmetric: F(-x) is not exactly -F(x). The sum F(x) + F(-x) is + // the asymmetry, at most one LSB per rounding that lands on a tie and + // reaches the output; its maximum and its mean (the bias) are pinned. + TYPED_TEST(fft_fixed_point_test, RoundingBiasOnNegatedInputIsBounded) { + using cfg = TypeParam; + using s = typename cfg::sample; + const auto max_pin = pins().negation_sum_max_lsb; + const auto bias_pin = pins().negation_bias_lsb; + ASSERT_GT(max_pin, 0.0) << "unmeasured pin"; + ASSERT_GT(bias_pin, 0.0) << "unmeasured pin"; + double worst_max = 0.0; + double worst_bias = 0.0; + for (const std::size_t n : k_sweep_sizes) { + for (const double amplitude : {0.999, 0.1}) { + // 0.999 keeps -x representable (-INT_MIN would saturate). + const auto x = + tap::dsp::test::random_signal(n, 0x9E3779B9u ^ static_cast(n), amplitude); + std::vector neg(n); + for (std::size_t i = 0; i < n; ++i) { + neg[i] = static_cast(-x[i]); + } + for (const bool inverse : {false, true}) { + const auto a = inverse ? run_inverse(x) : run_forward(x); + const auto b = inverse ? run_inverse(neg) : run_forward(neg); + ASSERT_EQ(a.exponent, b.exponent) << "negation changed the exponent, n=" << n; + double max_sum = 0.0; + double mean = 0.0; + for (std::size_t i = 0; i < n; ++i) { + const double sum = static_cast(a.out[i]) + static_cast(b.out[i]); + max_sum = std::max(max_sum, std::fabs(sum)); + mean += sum; + } + mean = std::fabs(mean / static_cast(n)); + worst_max = std::max(worst_max, max_sum); + worst_bias = std::max(worst_bias, mean); + EXPECT_LE(max_sum, max_pin) << cfg::name() << (inverse ? " inverse" : " forward") << " n=" << n + << " amplitude " << amplitude; + EXPECT_LE(mean, bias_pin) << cfg::name() << (inverse ? " inverse" : " forward") << " n=" << n + << " amplitude " << amplitude; + } + } + } + std::printf("[ measured ] %s F(x)+F(-x): max %.2f LSB (pin %.2f), bias %.4f LSB (pin %.4f)\n", cfg::name(), + worst_max, max_pin, worst_bias, bias_pin); + } + + // ======================================================================== + // Q15 against Q31: the one sibling-profile comparison (Part 9, Rules). + // ======================================================================== + + // The same signal (a Q15 vector, widened to Q31 exactly) through both + // profiles under the same policy: brought to a common exponent, the two + // spectra differ by the Q15 output quantisation plus the Q15 profile's + // (coarser) internal noise, pinned in Q15 output LSBs. + template + double q15_vs_q31(const char* what, double pin) { + using c15 = config; + using c31 = config; + double worst = 0.0; + for (const std::size_t n : k_sweep_sizes) { + for (const double amplitude : {1.0, 0.05}) { + const auto x15 = tap::dsp::test::random_signal( + n, 0xC0FFEE01u ^ static_cast(n), amplitude); + std::vector x31(n); + for (std::size_t i = 0; i < n; ++i) { + x31[i] = static_cast(x15[i]) << 16; // exact + } + const auto r15 = run_forward(x15); + const auto r31 = run_forward(x31); + // Compare in the unnormalised frame, in Q15 output LSBs at e15. + const double unit = c15::k_lsb * pow2(r15.exponent); + for (std::size_t i = 0; i < n; ++i) { + const double v15 = sample_scale::to_double(r15.out[i]) * pow2(r15.exponent); + const double v31 = sample_scale::to_double(r31.out[i]) * pow2(r31.exponent); + const double diff = std::fabs(v15 - v31) / unit; + worst = std::max(worst, diff); + EXPECT_LE(diff, pin) << what << " n=" << n << " amplitude " << amplitude << " index " << i + << " e15=" << r15.exponent << " e31=" << r31.exponent; + } + } + } + std::printf("[ measured ] Q15 vs Q31 (%s): max %.3f Q15 LSB (pin %.3f)\n", what, worst, pin); + return worst; + } + + TEST(fft_fixed_profiles, Q15AndQ31AgreeToTheQ15Floor) { + ASSERT_GT(k_pins_q15_fixed.q15_vs_q31_lsb, 0.0) << "unmeasured pin"; + ASSERT_GT(k_pins_q15_bfp.q15_vs_q31_lsb, 0.0) << "unmeasured pin"; + q15_vs_q31("fixed", k_pins_q15_fixed.q15_vs_q31_lsb); + q15_vs_q31("bfp", k_pins_q15_bfp.q15_vs_q31_lsb); + } + + // ======================================================================== + // The noise floor against Welch's model. + // ======================================================================== + // + // THE MODEL (Welch 1969, "A fixed-point fast Fourier transform error + // analysis"; Oppenheim & Weinstein 1972), specialised to this kernel's + // arithmetic as fft_arith.h states it: + // + // - Internal LSB q: 2^-31 for Q31, 2^-29 for the widened Q15 (Q2.29), + // both as fractions of full scale. Every rounding (shr_round, + // mul_coeff) is additive noise uniform in (-q/2, q/2], variance + // q^2/12, independent of every other rounding. + // - Structure (Part 7): a radix-4 DIF complex FFT of length M = N/2 + // (floor(log2 M / 2) radix-4 stages, one radix-2 stage when log2 M is + // odd), then Ooura's real post-pass, which pairs bins k and M-k and + // applies one complex product per pair. Under fixed scaling the Q31 + // profile pre-shifts its input by one bit (one rounding, no sum). + // - Per stage, per output component. Shift-before-butterfly: each of + // the stage's sum_gain inputs (4, 2, 2 for radix-4, radix-2, the + // post-pass) is rounded once per component and reaches every output + // component through the butterfly with unit weight, so the shift + // injects sum_gain * q^2/12 -- when the stage shifts at all (a BFP + // stage shifting by 0 rounds nothing: shr_round(x, 0) is exact). The + // twiddle product is the two-rounding complex multiply: 2 q^2/12 per + // component. The twiddle's own quantisation, |dw| <= 2^-31 per + // component (0.5 LSB of Q1.30), uniform, contributes 2 * P * 2^-62 / 3 + // where P is the signal variance per component entering the product. + // - Propagation. Noise present after stage t reaches the output through + // each later stage u with gain sum_gain_u * 4^-shift_u: a sum of + // sum_gain_u uncorrelated terms through unit-magnitude twiddles, then + // the shift. Signal: a white input of variance s^2 per sample enters + // as M complex values of variance s^2 per component and follows the + // same gain, without the injections. + // - The Q15 narrowing adds (2^-15)^2 / 12 per output value. + // - Shift schedule. scaling::fixed: every stage shifts its full amount + // (2, 1, 1; plus the Q31 pre-shift) and the sum is the fixed exponent. + // scaling::block_floating: the kernel returns only the total e, and + // the schedule is not observable, so the model takes the EARLY + // schedule (e assigned to the first stages, each up to its fixed + // amount), which is the worst case for noise: a bit shifted early + // attenuates nothing injected after it. The late schedule (what an + // on-demand scaler does for a quiet input) is the best case, and the + // table prints both so the measured floor can be read against them. + // + // The model's number is the noise variance per output component in the + // output's own units (fractions of full scale at exponent e). The + // measurement is the mean over the interior bins (indices 2 .. N-1; DC + // and Nyquist take a different post-pass path) of |out - G(x)/2^e|^2 per + // component, with G the double golden model on the SAME quantised input, + // so only the kernel's roundings and its twiddle quantisation are in it. + // The pin is the largest measured/model ratio over the sizes and levels, + // per configuration and material, at 2x; the table below is printed for + // a -V run to record (Part 13). + struct stage_spec { + int sum_gain; ///< inputs summed into each output component + int max_shift; ///< the fixed-scaling shift, and the BFP maximum + int product_roundings; ///< roundings per output component in the twiddle product + }; + + std::vector kernel_stages(std::size_t n, bool with_pre_shift) { + std::vector stages; + if (with_pre_shift) { + stages.push_back({1, 1, 0}); + } + const int log2_m = log2_size(n / 2); + for (int s = 0; s < log2_m / 2; ++s) { + stages.push_back({4, 2, 2}); + } + if (log2_m % 2 == 1) { + stages.push_back({2, 1, 2}); + } + stages.push_back({2, 1, 2}); // the real post-pass + return stages; + } + + enum class schedule { fixed, early, late }; + + std::vector shifts_for(const std::vector& stages, int e, schedule which) { + std::vector shifts(stages.size(), 0); + if (which == schedule::fixed) { + for (std::size_t t = 0; t < stages.size(); ++t) { + shifts[t] = stages[t].max_shift; + } + return shifts; + } + int remaining = e; + if (which == schedule::early) { + for (std::size_t t = 0; t < stages.size() && remaining > 0; ++t) { + shifts[t] = std::min(stages[t].max_shift, remaining); + remaining -= shifts[t]; + } + } + else { + for (std::size_t t = stages.size(); t-- > 0 && remaining > 0;) { + shifts[t] = std::min(stages[t].max_shift, remaining); + remaining -= shifts[t]; + } + } + return shifts; + } + + struct welch_prediction { + double noise; ///< variance per output component, output units + double signal; ///< the model's signal variance per output component + }; + + welch_prediction welch_model(const std::vector& stages, const std::vector& shifts, double q_int, + double q_out_extra, double input_variance) { + constexpr double twiddle_lsb = 0x1p-31; // 0.5 LSB of Q1.30 + const double rounding = q_int * q_int / 12.0; + double noise = 0.0; + double signal = input_variance; + for (std::size_t t = 0; t < stages.size(); ++t) { + const stage_spec& st = stages[t]; + const double shift_gain = pow2(-2 * shifts[t]); + const double stage_gain = static_cast(st.sum_gain) * shift_gain; + noise = noise * stage_gain; + signal = signal * stage_gain; + if (shifts[t] > 0) { + noise += static_cast(st.sum_gain) * rounding; + } + if (st.product_roundings > 0) { + noise += static_cast(st.product_roundings) * rounding; + noise += 2.0 * signal * twiddle_lsb * twiddle_lsb / 3.0; + } + } + noise += q_out_extra * q_out_extra / 12.0; + return {noise, signal}; + } + + struct noise_row { + const char* material; + std::size_t n; + double level_db; + int exponent; + double signal_db; + double floor_db; + double model_db; + double model_late_db; + double ratio; + }; + + template + noise_row measure_noise_floor(const char* material, std::size_t n, double level_db, + const std::vector& x, double input_variance) { + const auto r = run_forward(x); + const auto golden = golden_forward(fractions(x)); + const double scale = pow2(-r.exponent); + + double noise = 0.0; + double signal = 0.0; + for (std::size_t i = 2; i < n; ++i) { + const double g = golden[i] * scale; + const double err = sample_scale::to_double(r.out[i]) - g; + noise += err * err; + signal += g * g; + } + const auto count = static_cast(n - 2); + noise /= count; + signal /= count; + + const double q_int = pow2(-Cfg::k_int_frac); + const double q_out_extra = Cfg::k_is_q15 ? Cfg::k_lsb : 0.0; + const auto stages = kernel_stages(n, !Cfg::k_is_bfp && Cfg::k_pre > 0); + const auto model_shifts = shifts_for(stages, r.exponent, Cfg::k_is_bfp ? schedule::early : schedule::fixed); + const auto late_shifts = shifts_for(stages, r.exponent, Cfg::k_is_bfp ? schedule::late : schedule::fixed); + const auto model = welch_model(stages, model_shifts, q_int, q_out_extra, input_variance); + const auto late = welch_model(stages, late_shifts, q_int, q_out_extra, input_variance); + + noise_row row{material, n, level_db, r.exponent, + db(signal), db(noise), db(model.noise), db(late.noise), + noise / model.noise}; + std::printf("[ floor ] %s %-6s N=%5zu %4.0f dBFS e=%2d signal %7.2f floor %7.2f model %7.2f (late %7.2f) " + "dBFS/component ratio %.3f snr %6.2f dB\n", + Cfg::name(), material, n, level_db, r.exponent, row.signal_db, row.floor_db, row.model_db, + row.model_late_db, row.ratio, row.signal_db - row.floor_db); + return row; + } + + TYPED_TEST(fft_fixed_point_test, NoiseFloorTracksWelchModel) { + using cfg = TypeParam; + using s = typename cfg::sample; + const auto noise_pin = pins().noise_ratio_noise; + const auto tone_pin = pins().noise_ratio_tone; + ASSERT_GT(noise_pin, 0.0) << "unmeasured pin"; + ASSERT_GT(tone_pin, 0.0) << "unmeasured pin"; + double worst_noise = 0.0; + double worst_tone = 0.0; + for (const std::size_t n : {std::size_t{256}, std::size_t{512}, std::size_t{2048}}) { + for (const double level_db : {0.0, -20.0, -40.0, -60.0}) { + const double amplitude = std::pow(10.0, level_db / 20.0); + // White noise, uniform in [-A, A): variance A^2 / 3 per sample. + const auto noise = + tap::dsp::test::random_signal(n, 0x2545F491u ^ static_cast(n), amplitude); + const auto nrow = measure_noise_floor("noise", n, level_db, noise, amplitude * amplitude / 3.0); + worst_noise = std::max(worst_noise, nrow.ratio); + // On-bin tone at bin N/8 + 3, peak A * full scale: variance A^2 / 2. + const double a = amplitude * cfg::k_full; + const auto tone = tap::dsp::test::tone(n, static_cast(n / 8 + 3), a, 0.3); + const auto trow = measure_noise_floor("tone", n, level_db, tone, a * a / 2.0); + worst_tone = std::max(worst_tone, trow.ratio); + } + } + EXPECT_LE(worst_noise, noise_pin) << cfg::name() << " white-noise floor above the pinned ratio to the model"; + EXPECT_LE(worst_tone, tone_pin) << cfg::name() << " on-bin tone floor above the pinned ratio to the model"; + // A floor far BELOW the model would mean a different arithmetic (a + // fused complex multiply, an exact shift) or a broken measurement, + // not a better kernel: the two-rounding form is the contract. + EXPECT_GE(worst_noise, noise_pin / 8.0) << cfg::name() << " white-noise floor implausibly far below the model"; + std::printf("[ measured ] %s noise/model ratio: noise %.3f (pin %.3f), tone %.3f (pin %.3f)\n", cfg::name(), + worst_noise, noise_pin, worst_tone, tone_pin); + } + + // ======================================================================== + // The twiddle table (fft/tables.h). + // ======================================================================== + + /// FNV-1a 64 over the table as little-endian int32 bytes: the fold is + /// over bit patterns, host-endianness independent, and one differing + /// last bit anywhere changes it (Part 13: an integer fold, never a float). + std::uint64_t fnv1a64(const std::vector& table) { + std::uint64_t h = 0xcbf29ce484222325ull; + for (const std::int32_t v : table) { + const auto u = static_cast(v); + for (int b = 0; b < 4; ++b) { + h ^= static_cast((u >> (8 * b)) & 0xffu); + h *= 0x100000001b3ull; + } + } + return h; + } + + /// The kernel's Q1.30 table for a transform of size n, as tables.h lays + /// it out: entry 2k is cos(2 pi k / n), entry 2k + 1 is sin(2 pi k / n). + std::vector twiddle_table(std::size_t n) { + return tap::dsp::detail::make_twiddle_table_q30(n); + } + + // |w_q - w| <= 0.5 LSB of Q1.30 on every host: the table is cos/sin from + // libm rounded once (make_coeff, half away from zero), never a recurrence. + // The reference is libm too, evaluated on the exact turn fraction; two + // double evaluations of the same angle differ by ~1e-16, far inside the + // 2^-20 LSB slack, so a last-bit libm difference passes here (the + // checksum below is where it is detected) while a recurrence-drifted or + // Q1.14-derived table fails. + TEST(fft_fixed_tables, TwiddleTableIsWithinHalfLsb) { + constexpr double one = 0x1p30; + for (std::size_t n = 4; n <= k_max_n; n *= 2) { + const auto table = twiddle_table(n); + ASSERT_EQ(table.size(), 2 * n) << "n=" << n; + double worst = 0.0; + for (std::size_t k = 0; k < n; ++k) { + const double turns = static_cast(k) / static_cast(n); // exact + const double c = std::cos(2.0 * std::numbers::pi * turns) * one; + const double s = std::sin(2.0 * std::numbers::pi * turns) * one; + const double ec = std::fabs(static_cast(table[2 * k]) - c); + const double es = std::fabs(static_cast(table[2 * k + 1]) - s); + worst = std::max({worst, ec, es}); + ASSERT_LE(ec, 0.5 + 0x1p-20) << "cos n=" << n << " k=" << k; + ASSERT_LE(es, 0.5 + 0x1p-20) << "sin n=" << n << " k=" << k; + } + // 1.0 and 0 are exact; the table is exactly symmetric where the + // angles are (k and n - k: same cos, opposite sin; n/4: (0, 1)). + EXPECT_EQ(table[0], std::int32_t{1} << 30); + EXPECT_EQ(table[1], 0); + EXPECT_EQ(table[2 * (n / 4)], 0); + EXPECT_EQ(table[2 * (n / 4) + 1], std::int32_t{1} << 30); + for (std::size_t k = 1; k < n / 2; ++k) { + ASSERT_EQ(table[2 * k], table[2 * (n - k)]) << "cos symmetry n=" << n << " k=" << k; + ASSERT_EQ(table[2 * k + 1], -table[2 * (n - k) + 1]) << "sin symmetry n=" << n << " k=" << k; + } + if (n <= 2048) { + std::printf("[ measured ] twiddle table n=%zu: max |w_q - w| = %.6f LSB\n", n, worst); + } + } + } + + // Host libm last-bit differences can move a double lying within 2^-31 of + // a Q1.30 rounding boundary to the other side (fft_arith.h, "Twiddles"), + // and then fixed-point outputs differ between hosts by design. The + // checksum makes that visible instead of absorbed: a disagreement between + // the CI hosts is a finding to record (which host, which N), not a skip. + TEST(fft_fixed_tables, TwiddleTableChecksumIsPinned) { + struct pinned { + std::size_t n; + std::uint64_t fnv1a64; + }; + // Not yet measured: filled from the real kernel branch; the host that + // produced each value is recorded here. + constexpr std::array expected{{{256, 0}, {512, 0}, {2048, 0}}}; + for (const auto& p : expected) { + const auto got = fnv1a64(twiddle_table(p.n)); + std::printf("[ checksum ] twiddle table n=%zu fnv1a64=%016llx\n", p.n, + static_cast(got)); + ASSERT_NE(p.fnv1a64, 0u) << "unmeasured checksum n=" << p.n; + EXPECT_EQ(got, p.fnv1a64) << "n=" << p.n << ": this host's libm produced a different Q1.30 table"; + } + } + +} // namespace From 665a15bf803d3f7cf223d46f85b7ae1e6eb56d0b Mon Sep 17 00:00:00 2001 From: Timothy Place Date: Fri, 18 Sep 2026 01:35:44 +0000 Subject: [PATCH 2/4] tests: widen test_fft.cpp, the oracle and the rt guard to the fixed profiles test_fft.cpp's typed contract suite runs on float, double, Q15 and Q31 through a per-profile trait (input conversion, forward scale from the returned exponent, tolerance), with Q15TracksDouble / Q31TracksDouble beside FloatTracksDouble; the float/double comparisons are unchanged. test_fft_oracle.cpp gains profile / profile with a tolerance derived from fft_arith.h's rounding count. test_fft_rt.cpp is typed over the six FFT classes (both scaling policies) for noexcept, the allocation guard and copy bit-identity including the exponent. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019ZPTzNxo5Fe4EtpXXKf7Sy --- tests/test_fft.cpp | 341 +++++++++++++++++++++++++++++++------- tests/test_fft_fixed.cpp | 111 ++++++++----- tests/test_fft_oracle.cpp | 95 +++++++++-- tests/test_fft_rt.cpp | 156 +++++++++++------ 4 files changed, 530 insertions(+), 173 deletions(-) diff --git a/tests/test_fft.cpp b/tests/test_fft.cpp index 91a629a..0516bf6 100644 --- a/tests/test_fft.cpp +++ b/tests/test_fft.cpp @@ -4,44 +4,184 @@ // Locks down the tap::dsp::basic_real_fft contract: round-trip identity, the // packed-spectrum layout (DC / Nyquist in slots 0 and 1), the documented // W = exp(+2*pi*i/N) sign convention, Parseval energy conservation, and -// float/double cross-precision agreement. FDAF correctness later depends on -// every one of these staying exactly as documented in fft.h. +// cross-precision agreement against the double golden model. FDAF +// correctness later depends on every one of these staying exactly as +// documented in fft.h. +// +// The typed suite runs on all four profiles (Part 9 of +// docs/audit-fft-and-code-smells.md): float and double, and the Q15 / Q31 +// fixed-point profiles (Stage 3b) through the per-profile trait `profile` +// below, which supplies the input conversion, the forward scale read from +// the returned exponent, and the tolerance. The float and double checks are +// numerically what they were before the widening — same signals, same +// tolerances, same comparisons; the trait is the identity for them. #include +#include #include #include +#include #include #include +#include "support/signals.h" #include "tap/dsp/fft.h" namespace { + using tap::dsp::test::sample_scale; + template std::vector random_signal(size_t n, unsigned seed) { std::mt19937 gen(seed); std::uniform_real_distribution dist(-1.0, 1.0); std::vector x(n); for (auto& v : x) { - v = static_cast(dist(gen)); + v = sample_scale::from_double(dist(gen)); } return x; } + // Absolute tolerance in the sample's own units (fractions of full scale + // for the fixed profiles). Float and double: the golden battery's + // numbers, unchanged. Q15 / Q31: measured against the closed forms below + // and pinned at 2x; see profile<> for the exponent handling. + // MEASURED: not yet (the kernel branch has not landed); 0.0 fails on purpose. template constexpr double k_tolerance = 0.0; template <> constexpr double k_tolerance = 1e-12; template <> constexpr double k_tolerance = 2e-5; + template <> + constexpr double k_tolerance = 0.0; + template <> + constexpr double k_tolerance = 0.0; + + // ------------------------------------------------------------------------ + // Per-profile trait. A forward transform's result is read as + // to_double(data[i]) * 2^e, with e the exponent the fixed-point profiles + // return (0 for float and double, whose transforms return void); the + // inverse used here is the UNNORMALISED one in every profile, so a round + // trip is x == back * 2^(e_fwd + e_inv + 1 - log2 n) in all four (for the + // floating profiles that is the 2/n their inverse() applies). + // ------------------------------------------------------------------------ + template + struct profile { + using fft = tap::dsp::basic_real_fft; + + static constexpr bool k_fixed_point = std::is_integral_v; + + static double to_double(Sample v) { return sample_scale::to_double(v); } + static Sample from_double(double v) { return sample_scale::from_double(v); } + + static int forward(fft& f, Sample* data) { + if constexpr (k_fixed_point) { + return f.forward_inplace(data); + } + else { + f.forward_inplace(data); + return 0; + } + } + static int inverse_unnormalised(fft& f, Sample* data) { + if constexpr (k_fixed_point) { + return f.inverse_inplace(data); + } + else { + f.inverse_inplace(data); + return 0; + } + } + /// Tolerance on a value the golden battery compared at k_tolerance * n: + /// a spectrum value of magnitude ~n in the floating profiles, which the + /// fixed profiles hold scaled by 2^-e, so their tolerance does not grow. + static double tolerance_at_n(size_t n) { + if constexpr (k_fixed_point) { + return k_tolerance; + } + else { + return k_tolerance * static_cast(n); + } + } + }; template class real_fft_test : public ::testing::Test {}; - using sample_types = ::testing::Types; + using sample_types = ::testing::Types; TYPED_TEST_SUITE(real_fft_test, sample_types); + constexpr int log2_of(size_t n) { + int l = 0; + while ((size_t{1} << l) < n) { + ++l; + } + return l; + } + + /// forward() then the unnormalised inverse, reconstructed to the input's + /// scale: back[i] * 2^(e_fwd + e_inv + 1 - log2 n), as doubles. + template + std::vector round_trip(size_t n, const std::vector& x, bool in_place_aliased) { + using p = profile; + typename p::fft fft(n); + std::vector buf = x; + int e_fwd = 0; + int e_inv = 0; + if (in_place_aliased) { + // forward()/inverse() explicitly permit output aliasing input. + if constexpr (p::k_fixed_point) { + e_fwd = fft.forward(buf.data(), buf.data()); + e_inv = fft.inverse(buf.data(), buf.data()); + } + else { + fft.forward(buf.data(), buf.data()); + fft.inverse(buf.data(), buf.data()); // the 2/n-normalised inverse, as before + } + } + else { + std::vector spectrum(n); + if constexpr (p::k_fixed_point) { + e_fwd = fft.forward(x.data(), spectrum.data()); + e_inv = fft.inverse(spectrum.data(), buf.data()); + } + else { + fft.forward(x.data(), spectrum.data()); + std::vector back(n); + fft.inverse(spectrum.data(), back.data()); // the 2/n-normalised inverse, as before + buf = back; + } + } + std::vector out(n); + const double scale = p::k_fixed_point ? std::ldexp(1.0, e_fwd + e_inv + 1 - log2_of(n)) : 1.0; + for (size_t i = 0; i < n; ++i) { + out[i] = p::to_double(buf[i]) * scale; + } + return out; + } + + /// The round-trip tolerance: the profile's k_tolerance for float and + /// double (as before); for the fixed profiles the reconstruction unit is + /// the output LSB times the same power of two (under fixed scaling the + /// trip discards log2 n bits twice: Q15 at n = 1024 reconstructs in + /// steps of 2^-4, the honest number fft.h states), times the pinned + /// number of those units. MEASURED: not yet; 0.0 fails on purpose. + template + constexpr double k_round_trip_units = 0.0; + + template + double round_trip_tolerance(size_t n) { + if constexpr (profile::k_fixed_point) { + const int e = profile::fft::fixed_scaling_exponent(n); + return k_round_trip_units * sample_scale::k_lsb * std::ldexp(1.0, 2 * e + 1 - log2_of(n)); + } + else { + return k_tolerance; + } + } + TYPED_TEST(real_fft_test, SizeAndBinCount) { tap::dsp::basic_real_fft fft(1024); EXPECT_EQ(fft.size(), 1024u); @@ -51,71 +191,72 @@ namespace { TYPED_TEST(real_fft_test, RoundTripReproducesInput) { constexpr size_t n = 1024; - tap::dsp::basic_real_fft fft(n); - const auto x = random_signal(n, 42); - - std::vector spectrum(n); - std::vector back(n); - fft.forward(x.data(), spectrum.data()); - fft.inverse(spectrum.data(), back.data()); + const auto x = random_signal(n, 42); + const auto back = round_trip(n, x, false); + const auto tol = round_trip_tolerance(n); + ASSERT_GT(tol, 0.0) << "unmeasured pin"; for (size_t i = 0; i < n; ++i) { - EXPECT_NEAR(back[i], x[i], k_tolerance) << "sample " << i; + EXPECT_NEAR(back[i], profile::to_double(x[i]), tol) << "sample " << i; } } TYPED_TEST(real_fft_test, RoundTripInPlaceAndAliased) { constexpr size_t n = 256; - tap::dsp::basic_real_fft fft(n); - const auto x = random_signal(n, 7); - - // forward()/inverse() explicitly permit output aliasing input. - auto buf = x; - fft.forward(buf.data(), buf.data()); - fft.inverse(buf.data(), buf.data()); + const auto x = random_signal(n, 7); + const auto back = round_trip(n, x, true); + const auto tol = round_trip_tolerance(n); + ASSERT_GT(tol, 0.0) << "unmeasured pin"; for (size_t i = 0; i < n; ++i) { - EXPECT_NEAR(buf[i], x[i], k_tolerance) << "sample " << i; + EXPECT_NEAR(back[i], profile::to_double(x[i]), tol) << "sample " << i; } } TYPED_TEST(real_fft_test, ImpulseHasFlatSpectrum) { + using p = profile; constexpr size_t n = 64; - tap::dsp::basic_real_fft fft(n); - std::vector x(n, TypeParam(0)); - x[0] = TypeParam(1); + typename p::fft fft(n); + std::vector x(n, TypeParam(0)); + x[0] = p::from_double(1.0); - fft.forward_inplace(x.data()); + const int e = p::forward(fft, x.data()); + const double scale = std::ldexp(1.0, -e); + const double tol = k_tolerance; + ASSERT_GT(tol, 0.0) << "unmeasured pin"; - EXPECT_NEAR(x[0], 1.0, k_tolerance); // DC - EXPECT_NEAR(x[1], 1.0, k_tolerance); // Nyquist + EXPECT_NEAR(p::to_double(x[0]), 1.0 * scale, tol); // DC + EXPECT_NEAR(p::to_double(x[1]), 1.0 * scale, tol); // Nyquist for (size_t k = 1; k < n / 2; ++k) { - EXPECT_NEAR(x[2 * k], 1.0, k_tolerance) << "bin " << k << " real"; - EXPECT_NEAR(x[2 * k + 1], 0.0, k_tolerance) << "bin " << k << " imag"; + EXPECT_NEAR(p::to_double(x[2 * k]), 1.0 * scale, tol) << "bin " << k << " real"; + EXPECT_NEAR(p::to_double(x[2 * k + 1]), 0.0, tol) << "bin " << k << " imag"; } } TYPED_TEST(real_fft_test, DcAndNyquistPacking) { + using p = profile; constexpr size_t n = 64; - tap::dsp::basic_real_fft fft(n); + typename p::fft fft(n); + const double tol = p::tolerance_at_n(n); + ASSERT_GT(tol, 0.0) << "unmeasured pin"; // Constant input: all energy in DC = data[0]. - std::vector dc(n, TypeParam(1)); - fft.forward_inplace(dc.data()); - EXPECT_NEAR(dc[0], static_cast(n), k_tolerance * n); - EXPECT_NEAR(dc[1], 0.0, k_tolerance * n); + std::vector dc(n, p::from_double(1.0)); + const int e_dc = p::forward(fft, dc.data()); + EXPECT_NEAR(p::to_double(dc[0]), static_cast(n) * std::ldexp(1.0, -e_dc), tol); + EXPECT_NEAR(p::to_double(dc[1]), 0.0, tol); // Alternating +1/-1: all energy in Nyquist = data[1]. std::vector nyq(n); for (size_t i = 0; i < n; ++i) { - nyq[i] = (i % 2 == 0) ? TypeParam(1) : TypeParam(-1); + nyq[i] = (i % 2 == 0) ? p::from_double(1.0) : p::from_double(-1.0); } - fft.forward_inplace(nyq.data()); - EXPECT_NEAR(nyq[0], 0.0, k_tolerance * n); - EXPECT_NEAR(nyq[1], static_cast(n), k_tolerance * n); + const int e_nyq = p::forward(fft, nyq.data()); + EXPECT_NEAR(p::to_double(nyq[0]), 0.0, tol); + EXPECT_NEAR(p::to_double(nyq[1]), static_cast(n) * std::ldexp(1.0, -e_nyq), tol); } // The documented sign convention (W = exp(+2*pi*i/N)): a pure cosine at @@ -123,61 +264,115 @@ namespace { // in the bin's IMAG part with a PLUS sign (+N/2) — the conjugate of the // engineering-convention DFT, where it would be -N/2. TYPED_TEST(real_fft_test, SignConventionIsPlusI) { + using p = profile; constexpr size_t n = 128; constexpr size_t k = 5; - tap::dsp::basic_real_fft fft(n); - const double w = 2.0 * std::numbers::pi * static_cast(k) / static_cast(n); + typename p::fft fft(n); + const double w = 2.0 * std::numbers::pi * static_cast(k) / static_cast(n); std::vector cosine(n); std::vector sine(n); for (size_t j = 0; j < n; ++j) { - cosine[j] = static_cast(std::cos(w * static_cast(j))); - sine[j] = static_cast(std::sin(w * static_cast(j))); + cosine[j] = p::from_double(std::cos(w * static_cast(j))); + sine[j] = p::from_double(std::sin(w * static_cast(j))); } - fft.forward_inplace(cosine.data()); - fft.forward_inplace(sine.data()); - - const double half = static_cast(n) / 2.0; - const double tol = k_tolerance * static_cast(n); - EXPECT_NEAR(cosine[2 * k], half, tol); - EXPECT_NEAR(cosine[2 * k + 1], 0.0, tol); - EXPECT_NEAR(sine[2 * k], 0.0, tol); - EXPECT_NEAR(sine[2 * k + 1], half, tol); // +N/2, not -N/2 + const int e_cos = p::forward(fft, cosine.data()); + const int e_sin = p::forward(fft, sine.data()); + + const double half_cos = static_cast(n) / 2.0 * std::ldexp(1.0, -e_cos); + const double half_sin = static_cast(n) / 2.0 * std::ldexp(1.0, -e_sin); + const double tol = p::tolerance_at_n(n); + ASSERT_GT(tol, 0.0) << "unmeasured pin"; + EXPECT_NEAR(p::to_double(cosine[2 * k]), half_cos, tol); + EXPECT_NEAR(p::to_double(cosine[2 * k + 1]), 0.0, tol); + EXPECT_NEAR(p::to_double(sine[2 * k]), 0.0, tol); + EXPECT_NEAR(p::to_double(sine[2 * k + 1]), half_sin, tol); // +N/2, not -N/2 // And nothing leaks into any other bin. for (size_t bin = 1; bin < n / 2; ++bin) { if (bin == k) { continue; } - EXPECT_NEAR(cosine[2 * bin], 0.0, tol) << "cos leak, bin " << bin; - EXPECT_NEAR(sine[2 * bin], 0.0, tol) << "sin leak, bin " << bin; + EXPECT_NEAR(p::to_double(cosine[2 * bin]), 0.0, tol) << "cos leak, bin " << bin; + EXPECT_NEAR(p::to_double(sine[2 * bin]), 0.0, tol) << "sin leak, bin " << bin; } } + // Parseval's relative tolerance for the fixed profiles: the rounding + // noise adds power, so the relative error is of the order of the per-bin + // noise-to-signal ratio; measured and pinned at 2x. + // MEASURED: not yet; 0.0 fails on purpose. + template + constexpr double k_parseval_relative = 0.0; + TYPED_TEST(real_fft_test, ParsevalEnergyConservation) { + using p = profile; constexpr size_t n = 512; - tap::dsp::basic_real_fft fft(n); - const auto x = random_signal(n, 1234); + typename p::fft fft(n); + const auto x = random_signal(n, 1234); double time_energy = 0.0; for (size_t i = 0; i < n; ++i) { - time_energy += static_cast(x[i]) * static_cast(x[i]); + time_energy += p::to_double(x[i]) * p::to_double(x[i]); } - auto spectrum = x; - fft.forward_inplace(spectrum.data()); - double freq_energy = static_cast(spectrum[0]) * static_cast(spectrum[0]) - + static_cast(spectrum[1]) * static_cast(spectrum[1]); + auto spectrum = x; + const int e = p::forward(fft, spectrum.data()); + const double scale = std::ldexp(1.0, e); // spectrum values are X * 2^-e + double freq_energy = p::to_double(spectrum[0]) * scale * (p::to_double(spectrum[0]) * scale) + + p::to_double(spectrum[1]) * scale * (p::to_double(spectrum[1]) * scale); for (size_t k = 1; k < n / 2; ++k) { - const double re = static_cast(spectrum[2 * k]); - const double im = static_cast(spectrum[2 * k + 1]); + const double re = p::to_double(spectrum[2 * k]) * scale; + const double im = p::to_double(spectrum[2 * k + 1]) * scale; freq_energy += 2.0 * (re * re + im * im); } freq_energy /= static_cast(n); - EXPECT_NEAR(freq_energy, time_energy, k_tolerance * time_energy * static_cast(n)); + if constexpr (p::k_fixed_point) { + ASSERT_GT(k_parseval_relative, 0.0) << "unmeasured pin"; + EXPECT_NEAR(freq_energy, time_energy, k_parseval_relative * time_energy); + } + else { + EXPECT_NEAR(freq_energy, time_energy, k_tolerance * time_energy * static_cast(n)); + } + } + + // ------------------------------------------------------------------------ + // Cross-precision: every profile against the double golden model. + // ------------------------------------------------------------------------ + + /// The relative 2-norm error of a profile's forward transform against the + /// double golden model on the same signal: sqrt(sum |s - g * 2^-e|^2 / + /// sum |g * 2^-e|^2) over the packed spectrum, with g the double result + /// on the profile's own quantised input. + template + double relative_error_vs_double(size_t n, unsigned seed) { + using p = profile; + const auto x = random_signal(n, seed); + std::vector xd(n); + for (size_t i = 0; i < n; ++i) { + xd[i] = p::to_double(x[i]); + } + + tap::dsp::real_fft fft64(n); + typename p::fft fftp(n); + auto sd = xd; + auto sp = x; + fft64.forward_inplace(sd.data()); + const int e = p::forward(fftp, sp.data()); + const double scale = std::ldexp(1.0, -e); + + double err = 0.0; + double ref = 0.0; + for (size_t i = 0; i < n; ++i) { + const double g = sd[i] * scale; + const double d = p::to_double(sp[i]) - g; + err += d * d; + ref += g * g; + } + return std::sqrt(err / ref); } // The two precisions are the same algorithm on the same data: float must @@ -209,6 +404,28 @@ namespace { EXPECT_LT(std::sqrt(err / ref), 1e-6); } + // The fixed-point profiles against the same golden model: the relative + // 2-norm error of the fixed-scaling forward at n = 1024 on full-scale + // uniform noise, a measured number pinned at 2x (the log_mel pattern). + // Q15's error is the output narrowing (2^-15 / sqrt(12) per value on a + // spectrum whose per-component RMS is sqrt(1/3 / 2n)); Q31's is the + // int32 kernel's rounding noise, orders of magnitude lower. + // MEASURED: not yet; 0.0 fails on purpose. + constexpr double k_q15_tracks_double = 0.0; + constexpr double k_q31_tracks_double = 0.0; + + TEST(RealFftCrossPrecision, Q15TracksDouble) { + ASSERT_GT(k_q15_tracks_double, 0.0) << "unmeasured pin"; + const double err = relative_error_vs_double(1024, 99); + EXPECT_LT(err, k_q15_tracks_double) << "measured " << err; + } + + TEST(RealFftCrossPrecision, Q31TracksDouble) { + ASSERT_GT(k_q31_tracks_double, 0.0) << "unmeasured pin"; + const double err = relative_error_vs_double(1024, 99); + EXPECT_LT(err, k_q31_tracks_double) << "measured " << err; + } + // The double-engine float-I/O convenience overloads (used by AmbiTap's HRTF // analysis): a float-buffer round trip through the double FFT reproduces the // input to float precision. diff --git a/tests/test_fft_fixed.cpp b/tests/test_fft_fixed.cpp index 41156fe..e282558 100644 --- a/tests/test_fft_fixed.cpp +++ b/tests/test_fft_fixed.cpp @@ -1124,12 +1124,14 @@ namespace { } // ======================================================================== - // The twiddle table (fft/tables.h). + // The tables (fft/tables.h). // ======================================================================== - /// FNV-1a 64 over the table as little-endian int32 bytes: the fold is - /// over bit patterns, host-endianness independent, and one differing - /// last bit anywhere changes it (Part 13: an integer fold, never a float). + /// FNV-1a 64 over a table as little-endian int32 bytes, written here + /// independently of detail::table_checksum so that helper is itself + /// pinned: the fold is over bit patterns, host-endianness independent, + /// and one differing last bit anywhere changes it (Part 13: an integer + /// fold, never a float). std::uint64_t fnv1a64(const std::vector& table) { std::uint64_t h = 0xcbf29ce484222325ull; for (const std::int32_t v : table) { @@ -1142,47 +1144,61 @@ namespace { return h; } - /// The kernel's Q1.30 table for a transform of size n, as tables.h lays - /// it out: entry 2k is cos(2 pi k / n), entry 2k + 1 is sin(2 pi k / n). - std::vector twiddle_table(std::size_t n) { - return tap::dsp::detail::make_twiddle_table_q30(n); - } - - // |w_q - w| <= 0.5 LSB of Q1.30 on every host: the table is cos/sin from - // libm rounded once (make_coeff, half away from zero), never a recurrence. - // The reference is libm too, evaluated on the exact turn fraction; two - // double evaluations of the same angle differ by ~1e-16, far inside the - // 2^-20 LSB slack, so a last-bit libm difference passes here (the - // checksum below is where it is detected) while a recurrence-drifted or - // Q1.14-derived table fails. + // |w_q - w| <= 0.5 LSB of Q1.30 on every host, for the kernel twiddles + // (W_M^k, M = N/2, interleaved cos/sin) and the real post-pass pairs + // (0.5 - 0.5 sin, 0.5 cos over N): each is cos/sin from libm rounded + // once by make_coeff (half away from zero), never a recurrence. The + // reference is libm too, on the exact turn fraction; two double + // evaluations of the same angle differ by ~1e-16, far inside the 2^-20 + // LSB slack, so a last-bit libm difference passes here (the checksum + // below is where it is detected) while a recurrence-drifted or + // Q1.14-derived table fails. The exact entries (1, 0, i) and the + // k <-> M-k symmetry are checked bit-exactly. TEST(fft_fixed_tables, TwiddleTableIsWithinHalfLsb) { - constexpr double one = 0x1p30; + constexpr double one = 0x1p30; + constexpr double slack = 0.5 + 0x1p-20; for (std::size_t n = 4; n <= k_max_n; n *= 2) { - const auto table = twiddle_table(n); - ASSERT_EQ(table.size(), 2 * n) << "n=" << n; + const std::size_t m = n / 2; + const auto table = tap::dsp::detail::make_twiddle_table(m); + ASSERT_EQ(table.size(), 2 * m) << "n=" << n; double worst = 0.0; - for (std::size_t k = 0; k < n; ++k) { - const double turns = static_cast(k) / static_cast(n); // exact + for (std::size_t k = 0; k < m; ++k) { + const double turns = static_cast(k) / static_cast(m); // exact const double c = std::cos(2.0 * std::numbers::pi * turns) * one; const double s = std::sin(2.0 * std::numbers::pi * turns) * one; const double ec = std::fabs(static_cast(table[2 * k]) - c); const double es = std::fabs(static_cast(table[2 * k + 1]) - s); worst = std::max({worst, ec, es}); - ASSERT_LE(ec, 0.5 + 0x1p-20) << "cos n=" << n << " k=" << k; - ASSERT_LE(es, 0.5 + 0x1p-20) << "sin n=" << n << " k=" << k; + ASSERT_LE(ec, slack) << "cos m=" << m << " k=" << k; + ASSERT_LE(es, slack) << "sin m=" << m << " k=" << k; + } + EXPECT_EQ(table[0], std::int32_t{1} << 30) << "m=" << m; + EXPECT_EQ(table[1], 0) << "m=" << m; + if (m >= 4) { + EXPECT_EQ(table[2 * (m / 4)], 0) << "m=" << m; + EXPECT_EQ(table[2 * (m / 4) + 1], std::int32_t{1} << 30) << "m=" << m; } - // 1.0 and 0 are exact; the table is exactly symmetric where the - // angles are (k and n - k: same cos, opposite sin; n/4: (0, 1)). - EXPECT_EQ(table[0], std::int32_t{1} << 30); - EXPECT_EQ(table[1], 0); - EXPECT_EQ(table[2 * (n / 4)], 0); - EXPECT_EQ(table[2 * (n / 4) + 1], std::int32_t{1} << 30); - for (std::size_t k = 1; k < n / 2; ++k) { - ASSERT_EQ(table[2 * k], table[2 * (n - k)]) << "cos symmetry n=" << n << " k=" << k; - ASSERT_EQ(table[2 * k + 1], -table[2 * (n - k) + 1]) << "sin symmetry n=" << n << " k=" << k; + for (std::size_t k = 1; k < m / 2; ++k) { + ASSERT_EQ(table[2 * k], table[2 * (m - k)]) << "cos symmetry m=" << m << " k=" << k; + ASSERT_EQ(table[2 * k + 1], -table[2 * (m - k) + 1]) << "sin symmetry m=" << m << " k=" << k; + } + + const auto post = tap::dsp::detail::make_real_post_pass_table(n); + ASSERT_EQ(post.size(), 2 * (n / 4)) << "n=" << n; + EXPECT_EQ(post[0], std::int32_t{1} << 29) << "n=" << n; // 0.5 - 0.5 sin 0 + EXPECT_EQ(post[1], std::int32_t{1} << 29) << "n=" << n; // 0.5 cos 0 + for (std::size_t k = 0; k < n / 4; ++k) { + const double turns = static_cast(k) / static_cast(n); + const double wkr = (0.5 - 0.5 * std::sin(2.0 * std::numbers::pi * turns)) * one; + const double wki = 0.5 * std::cos(2.0 * std::numbers::pi * turns) * one; + const double er = std::fabs(static_cast(post[2 * k]) - wkr); + const double ei = std::fabs(static_cast(post[2 * k + 1]) - wki); + worst = std::max({worst, er, ei}); + ASSERT_LE(er, slack) << "wkr n=" << n << " k=" << k; + ASSERT_LE(ei, slack) << "wki n=" << n << " k=" << k; } if (n <= 2048) { - std::printf("[ measured ] twiddle table n=%zu: max |w_q - w| = %.6f LSB\n", n, worst); + std::printf("[ measured ] tables n=%zu: max |w_q - w| = %.6f LSB\n", n, worst); } } } @@ -1192,20 +1208,33 @@ namespace { // and then fixed-point outputs differ between hosts by design. The // checksum makes that visible instead of absorbed: a disagreement between // the CI hosts is a finding to record (which host, which N), not a skip. + // Both tables a transform of N builds are pinned, for the three + // certified geometries. TEST(fft_fixed_tables, TwiddleTableChecksumIsPinned) { struct pinned { std::size_t n; - std::uint64_t fnv1a64; + std::uint64_t twiddles; ///< make_twiddle_table(n / 2) + std::uint64_t post_pass; ///< make_real_post_pass_table(n) }; // Not yet measured: filled from the real kernel branch; the host that // produced each value is recorded here. - constexpr std::array expected{{{256, 0}, {512, 0}, {2048, 0}}}; + constexpr std::array expected{{{256, 0, 0}, {512, 0, 0}, {2048, 0, 0}}}; for (const auto& p : expected) { - const auto got = fnv1a64(twiddle_table(p.n)); - std::printf("[ checksum ] twiddle table n=%zu fnv1a64=%016llx\n", p.n, - static_cast(got)); - ASSERT_NE(p.fnv1a64, 0u) << "unmeasured checksum n=" << p.n; - EXPECT_EQ(got, p.fnv1a64) << "n=" << p.n << ": this host's libm produced a different Q1.30 table"; + const auto twiddles = tap::dsp::detail::make_twiddle_table(p.n / 2); + const auto post = tap::dsp::detail::make_real_post_pass_table(p.n); + const auto tw_sum = fnv1a64(twiddles); + const auto post_sum = fnv1a64(post); + // Printed before any assertion so a -V log carries every host's + // values whether or not they match. + std::printf("[ checksum ] n=%zu twiddles fnv1a64=%016llx post-pass fnv1a64=%016llx\n", p.n, + static_cast(tw_sum), static_cast(post_sum)); + // The kernel's own checksum helper is the same fold. + EXPECT_EQ(tap::dsp::detail::table_checksum(twiddles.data(), twiddles.size()), tw_sum); + EXPECT_EQ(tap::dsp::detail::table_checksum(post.data(), post.size()), post_sum); + EXPECT_NE(p.twiddles, 0u) << "unmeasured checksum n=" << p.n; + EXPECT_EQ(tw_sum, p.twiddles) << "n=" << p.n << ": this host's libm produced a different twiddle table"; + EXPECT_EQ(post_sum, p.post_pass) + << "n=" << p.n << ": this host's libm produced a different post-pass table"; } } diff --git a/tests/test_fft_oracle.cpp b/tests/test_fft_oracle.cpp index 76845d8..3dcbf43 100644 --- a/tests/test_fft_oracle.cpp +++ b/tests/test_fft_oracle.cpp @@ -34,11 +34,14 @@ // reference; the test names are distinct so the ctest listing carries each // promise once. // -// The suite is typed over float and double now. The fixed-point stage -// (Stage 3, Part 7) extends `profile` below with int16_t and int32_t; -// the tests are written against that trait (scale as a function of N, the -// profile's own tolerance, a full-scale amplitude below saturation), so the -// change is confined to the trait and the type list. +// The suite is typed over all four profiles: float and double (the Ooura +// engines) and, since Stage 3b (Part 7), int16_t and int32_t, the Q15 and +// Q31 fixed-point profiles under scaling::fixed, through `profile` +// below (scale as a function of N from the fixed exponent, the profile's own +// tolerance derived from the arithmetic's rounding count, a full-scale +// amplitude below saturation). Block floating point is not an oracle +// question (its exponent is data-dependent); test_fft_fixed.cpp pins it +// against the fixed result. #include #include @@ -46,12 +49,14 @@ #include #include #include +#include #include #include #include "support/signals.h" #include "tap/dsp/fft.h" +#include "tap/dsp/fft/fft_arith.h" #ifndef TAP_DSP_PARITY_MAX_N #define TAP_DSP_PARITY_MAX_N (1 << 20) @@ -107,20 +112,19 @@ namespace { } // ------------------------------------------------------------------------ - // Per-profile traits — THE EXTENSION POINT FOR THE FIXED-POINT STAGE. + // Per-profile traits. // // A profile says how a sample crosses into the double domain, the amplitude // the closed forms are driven at, what scale the engine's forward and // unnormalized inverse carry relative to the mathematical DFT as a function // of N, and its own tolerance. For the float profiles: full scale 1.0, both - // scales 1, and the Higham bound above. Stage 3 (Part 7) adds - // template <> struct profile { ... }; - // template <> struct profile { ... }; - // with k_full_scale below 1 - 2^-15 (a full-scale 1.0 input whose X/N - // expectation is exactly 1.0 saturates in Q15), forward_scale(n) = 1/n - // under `fixed` scaling (BFP reports its exponent alongside), the profile's - // inverse scale, and a tolerance built from Part 7's quantization-noise - // numbers rather than from epsilon. Then append the types to oracle_types. + // scales 1, and the Higham bound above. For the fixed-point profiles + // (Stage 3b): k_full_scale is the largest representable value 1 - 2^-15 / + // 1 - 2^-31 (a 1.0 whose X/N expectation is exactly 1.0 does not exist in + // Q15), both scales are 2^-e with e = fixed_scaling_exponent(n) (Q15: + // 1/n; Q31: 1/2n, the input pre-shift), and the tolerance is the + // worst-case rounding bound below, built from fft_arith.h's rounding + // count rather than from an epsilon. // ------------------------------------------------------------------------ template struct profile; @@ -172,7 +176,68 @@ namespace { static double tolerance(std::size_t n, double norm2) { return higham_tolerance(k_epsilon, n, norm2); } }; - using oracle_types = ::testing::Types; + // ------------------------------------------------------------------------ + // Fixed-point tolerance: a worst-case (max-abs) bound, derived. + // + // fft_arith.h fixes the arithmetic: shift-before-butterfly with + // round-half-up, the two-rounding complex multiply, Q1.30 twiddles rounded + // once. Each rounding is off by at most half an internal LSB q (2^-29 for + // the widened Q15 data, 2^-31 for Q31, as fractions of full scale) and + // each twiddle by at most 2^-31 relative. Per kernel stage, per output + // component, with every error aligned (the worst case, not the RMS): + // - the shift roundings of the stage's four inputs, each <= q/2, reach + // the output through the butterfly with unit weight: <= 2 q; + // - the two product roundings: <= q; + // - the twiddle quantisation on a value of magnitude <= 2^30.5 q + // (the fixed-scaling bound): <= sqrt(2) * 2^30.5 * 2^-31 q < 1 q; + // so <= 4 q per stage. Under shift-before-butterfly a stage attenuates + // the noise it receives by the same factor the sum can amplify it, so the + // worst-case bound simply adds across the ceil(log2(N/2) / 2) kernel + // stages, plus 4 q for the real post-pass (one shift, one complex + // product) and q/2 for the Q31 input pre-shift. On top: the input + // quantisation (<= half an input LSB per sample, summed over N terms and + // scaled by 2^-e: <= q_io / 2), and the output narrowing (Q15: half an + // output LSB). The bound is doubled for margin -- derived with slack, not + // fitted -- and the pins in test_fft.cpp / test_fft_fixed.cpp carry the + // measured numbers. + // ------------------------------------------------------------------------ + template + struct fixed_profile { + using fft = tap::dsp::basic_real_fft; + static constexpr int k_io_bits = std::numeric_limits::digits; // 15 or 31 + static constexpr int k_internal_bits = std::is_same_v ? 29 : 31; + static constexpr double k_q_io = tap::dsp::test::sample_scale::k_lsb; + static constexpr double k_q_int = 1.0 / static_cast(std::int64_t{1} << k_internal_bits); + static constexpr double k_full_scale = tap::dsp::test::sample_scale::k_full_scale; + static constexpr std::size_t k_min_n = 4; + static constexpr std::size_t k_max_n = 65536; + static constexpr double k_roundings_per_stage = 4.0; + static constexpr double k_post_pass = 4.0; + static constexpr double k_margin = 2.0; + + static double to_double(I v) { return tap::dsp::test::sample_scale::to_double(v); } + static I from_double(double v) { return tap::dsp::test::sample_scale::from_double(v); } + static double forward_scale(std::size_t n) { return std::ldexp(1.0, -fft::fixed_scaling_exponent(n)); } + static double inverse_scale(std::size_t n) { return std::ldexp(1.0, -fft::fixed_scaling_exponent(n)); } + + static int kernel_stages(std::size_t n) { + const int log2_m = static_cast(std::lround(std::log2(static_cast(n / 2)))); + return (log2_m + 1) / 2; // ceil(log2 M / 2): radix-4 stages plus the odd radix-2 + } + static double tolerance(std::size_t n, double /*norm2*/) { + const double pre_shift = tap::dsp::fft_arith::k_fixed_scaling_input_pre_shift * 0.5; + const double internal = (k_roundings_per_stage * kernel_stages(n) + k_post_pass + pre_shift) * k_q_int; + const double io = 0.5 * k_q_io + (k_internal_bits == k_io_bits ? 0.0 : 0.5 * k_q_io); + return k_margin * (internal + io); + } + }; + + template <> + struct profile : fixed_profile {}; + template <> + struct profile : fixed_profile {}; + + using oracle_types = ::testing::Types; /// ||y||_2 of the full complex spectrum of x, from Parseval: sqrt(N) * ||x||_2. double spectrum_norm2(const std::vector& x) { diff --git a/tests/test_fft_rt.cpp b/tests/test_fft_rt.cpp index 3258bb8..0c5ad21 100644 --- a/tests/test_fft_rt.cpp +++ b/tests/test_fft_rt.cpp @@ -9,7 +9,8 @@ // makes both halves of that promise checkable: // // - noexcept is a static_assert on the four transform entry points, on -// both instantiations (the pattern test_nn.cpp uses); +// all six instantiations -- float, double, and the Q15 / Q31 fixed-point +// profiles under both scaling policies (the pattern test_nn.cpp uses); // - "allocation-free" is checked by REPLACING THE GLOBAL operator new / // operator delete in this translation unit with counting versions and // asserting the count does not move across a transform. The replacement @@ -46,7 +47,9 @@ // removing. They are also not noexcept, consistently with that. // // Copy and copy-assignment are pinned to produce bit-identical output to the -// source in both directions and both precisions. test_fft_backend.cpp's +// source in both directions and all six instantiations (for the fixed-point +// profiles that includes the returned exponent, and the Q15 profile's int32 +// work buffer travelling with the copy). test_fft_backend.cpp's // CopiesAgreeWithTheirSource covers the float FORWARD at the certified // geometries while walking the heap (its purpose is vDSP's alignment // dispatch); this one is the general value-semantics check that the engine's @@ -62,6 +65,7 @@ #include #include #include +#include #include #include @@ -209,25 +213,41 @@ void operator delete[](void* p, std::align_val_t, const std::nothrow_t&) noexcep namespace { - template - using fft_t = tap::dsp::basic_real_fft; + // The suite is typed over the FFT class itself, so the two scaling + // policies of each fixed-point profile are distinct rows. + template + struct sample_of; + template + struct sample_of> { + using type = Sample; + }; + template + using sample_of_t = typename sample_of::type; // ------------------------------------------------------------------------ - // noexcept, as a compile-time fact on both instantiations. + // noexcept, as a compile-time fact on every instantiation. // ------------------------------------------------------------------------ - template + template constexpr bool transforms_are_noexcept() { - using fft = fft_t; - static_assert(noexcept(std::declval().forward_inplace(std::declval()))); - static_assert(noexcept(std::declval().inverse_inplace(std::declval()))); - static_assert(noexcept(std::declval().forward(std::declval(), std::declval()))); - static_assert(noexcept(std::declval().inverse(std::declval(), std::declval()))); + using fft = Fft; + using sample = sample_of_t; + static_assert(noexcept(std::declval().forward_inplace(std::declval()))); + static_assert(noexcept(std::declval().inverse_inplace(std::declval()))); + static_assert(noexcept(std::declval().forward(std::declval(), std::declval()))); + static_assert(noexcept(std::declval().inverse(std::declval(), std::declval()))); static_assert(noexcept(std::declval().size())); static_assert(noexcept(std::declval().num_bins())); return true; } - static_assert(transforms_are_noexcept()); - static_assert(transforms_are_noexcept()); + static_assert(transforms_are_noexcept()); + static_assert(transforms_are_noexcept()); + static_assert(transforms_are_noexcept()); + static_assert(transforms_are_noexcept()); + static_assert(transforms_are_noexcept()); + static_assert(transforms_are_noexcept()); + // The fixed-point profiles' constant exponent is a noexcept constexpr too. + static_assert(noexcept(tap::dsp::real_fft_q15::fixed_scaling_exponent(4))); + static_assert(noexcept(tap::dsp::real_fft_q31_bfp::fixed_scaling_exponent(4))); // ------------------------------------------------------------------------ // The allocation guard. @@ -246,11 +266,12 @@ namespace { constexpr std::size_t k_guarded_sizes[] = {512, 4096}; - template + template class fft_rt_test : public ::testing::Test {}; - using sample_types = ::testing::Types; - TYPED_TEST_SUITE(fft_rt_test, sample_types); + using fft_types = ::testing::Types; + TYPED_TEST_SUITE(fft_rt_test, fft_types); // Already proved by the namespace-scope static_asserts above; this exists // so the promise has a row in the test listing per profile. @@ -269,11 +290,12 @@ namespace { EXPECT_GE(guard.allocations_since(), 1u); } - template + template void expect_no_allocation(std::size_t n, const char* what, Call&& call) { - fft_t fft(n); - std::vector in = tap::dsp::test::random_signal(n, 0x9E3779B9u); - std::vector out = in; + using sample = sample_of_t; + Fft fft(n); + std::vector in = tap::dsp::test::random_signal(n, 0x9E3779B9u); + std::vector out = in; // First call after construction (lazy table init today, F6). { @@ -293,32 +315,34 @@ namespace { } TYPED_TEST(fft_rt_test, ForwardInplaceAllocatesNothing) { + using sample = sample_of_t; for (const std::size_t n : k_guarded_sizes) { expect_no_allocation( - n, "forward_inplace", - [](fft_t& fft, TypeParam*, TypeParam* out) { fft.forward_inplace(out); }); + n, "forward_inplace", [](TypeParam& fft, sample*, sample* out) { (void)fft.forward_inplace(out); }); } } TYPED_TEST(fft_rt_test, InverseInplaceAllocatesNothing) { + using sample = sample_of_t; for (const std::size_t n : k_guarded_sizes) { expect_no_allocation( - n, "inverse_inplace", - [](fft_t& fft, TypeParam*, TypeParam* out) { fft.inverse_inplace(out); }); + n, "inverse_inplace", [](TypeParam& fft, sample*, sample* out) { (void)fft.inverse_inplace(out); }); } } TYPED_TEST(fft_rt_test, ForwardOutOfPlaceAllocatesNothing) { + using sample = sample_of_t; for (const std::size_t n : k_guarded_sizes) { expect_no_allocation( - n, "forward", [](fft_t& fft, const TypeParam* in, TypeParam* out) { fft.forward(in, out); }); + n, "forward", [](TypeParam& fft, const sample* in, sample* out) { (void)fft.forward(in, out); }); } } TYPED_TEST(fft_rt_test, InverseOutOfPlaceAllocatesNothing) { + using sample = sample_of_t; for (const std::size_t n : k_guarded_sizes) { expect_no_allocation( - n, "inverse", [](fft_t& fft, const TypeParam* in, TypeParam* out) { fft.inverse(in, out); }); + n, "inverse", [](TypeParam& fft, const sample* in, sample* out) { (void)fft.inverse(in, out); }); } } @@ -339,67 +363,89 @@ namespace { struct both_directions { std::vector spectrum; std::vector time; + int forward_exponent = 0; ///< fixed-point profiles; 0 for float/double + int inverse_exponent = 0; }; - template - both_directions run_both(fft_t& fft, const std::vector& x) { - both_directions r; - r.spectrum = x; - fft.forward_inplace(r.spectrum.data()); - r.time = r.spectrum; - fft.inverse_inplace(r.time.data()); + /// The returned exponent, or 0 where the transform returns void. + template + int exponent_of(Call&& call) { + if constexpr (std::is_void_v) { + call(); + return 0; + } + else { + return call(); + } + } + + template + both_directions> run_both(Fft& fft, const std::vector>& x) { + both_directions> r; + r.spectrum = x; + r.forward_exponent = exponent_of([&] { return fft.forward_inplace(r.spectrum.data()); }); + r.time = r.spectrum; + r.inverse_exponent = exponent_of([&] { return fft.inverse_inplace(r.time.data()); }); return r; } + template + void expect_same_result(const both_directions& a, const both_directions& b, const char* what) { + expect_bit_identical(a.spectrum, b.spectrum, what); + expect_bit_identical(a.time, b.time, what); + EXPECT_EQ(a.forward_exponent, b.forward_exponent) << what << ": forward exponent differs"; + EXPECT_EQ(a.inverse_exponent, b.inverse_exponent) << what << ": inverse exponent differs"; + } + TYPED_TEST(fft_rt_test, CopyProducesBitIdenticalOutput) { + using sample = sample_of_t; constexpr std::size_t n = 512; - const auto x = tap::dsp::test::random_signal(n, 0x2545F491u); + const auto x = tap::dsp::test::random_signal(n, 0x2545F491u); - fft_t original(n); + TypeParam original(n); // Source already warmed (tables built) so the copy carries built tables. const auto from_original = run_both(original, x); - fft_t copy(original); - const auto from_copy = run_both(copy, x); - expect_bit_identical(from_copy.spectrum, from_original.spectrum, "copy forward"); - expect_bit_identical(from_copy.time, from_original.time, "copy inverse"); + TypeParam copy(original); + const auto from_copy = run_both(copy, x); + expect_same_result(from_copy, from_original, "copy"); // And the source is unaffected by having been copied. const auto again = run_both(original, x); - expect_bit_identical(again.spectrum, from_original.spectrum, "source forward after copy"); - expect_bit_identical(again.time, from_original.time, "source inverse after copy"); + expect_same_result(again, from_original, "source after copy"); } TYPED_TEST(fft_rt_test, CopyOfUnwarmedSourceProducesBitIdenticalOutput) { + using sample = sample_of_t; constexpr std::size_t n = 512; - const auto x = tap::dsp::test::random_signal(n, 0x2545F491u); + const auto x = tap::dsp::test::random_signal(n, 0x2545F491u); // Copied BEFORE any transform: with lazy tables (F6) both objects // build their own; with constructor-built tables both carry the same. - fft_t original(n); - fft_t copy(original); - const auto from_copy = run_both(copy, x); - const auto from_original = run_both(original, x); - expect_bit_identical(from_copy.spectrum, from_original.spectrum, "unwarmed copy forward"); - expect_bit_identical(from_copy.time, from_original.time, "unwarmed copy inverse"); + TypeParam original(n); + TypeParam copy(original); + const auto from_copy = run_both(copy, x); + const auto from_original = run_both(original, x); + expect_same_result(from_copy, from_original, "unwarmed copy"); } TYPED_TEST(fft_rt_test, CopyAssignmentProducesBitIdenticalOutput) { + using sample = sample_of_t; constexpr std::size_t n = 512; - const auto x = tap::dsp::test::random_signal(n, 0x2545F491u); + const auto x = tap::dsp::test::random_signal(n, 0x2545F491u); - fft_t original(n); - const auto from_original = run_both(original, x); + TypeParam original(n); + const auto from_original = run_both(original, x); // Assigned over an engine of a DIFFERENT geometry, so the assignment - // has to replace every table, not just refresh one of the same size. - fft_t target(4096); + // has to replace every table (and the Q15 work buffer), not just + // refresh one of the same size. + TypeParam target(4096); target = original; EXPECT_EQ(target.size(), n); EXPECT_EQ(target.num_bins(), n / 2 + 1); const auto from_target = run_both(target, x); - expect_bit_identical(from_target.spectrum, from_original.spectrum, "assigned forward"); - expect_bit_identical(from_target.time, from_original.time, "assigned inverse"); + expect_same_result(from_target, from_original, "assigned"); } } // namespace From 8b86d48ebdd0ea5f9450c92cb4a49023f554d979 Mon Sep 17 00:00:00 2001 From: Timothy Place Date: Fri, 18 Sep 2026 01:48:22 +0000 Subject: [PATCH 3/4] tests: pin the fixed-point battery at the kernel's measured numbers Measured 2026-09-18 on x86-64 Linux (GCC 13.3.0 and clang 18.1.3 agree) against claude/wave2-stage3b-kernel a55a14f and pinned at 2x: the saturation-sweep deviation, round-trip error, BFP-vs-fixed shift agreement, negation asymmetry and bias, Welch-model ratios for white noise and an on-bin tone, Q15-vs-Q31 agreement, the Q15/Q31 contract tolerances, Parseval and TracksDouble numbers, and the Q1.30 table checksums for N = 256 / 512 / 2048. The Welch model is refined to the kernel's stated structure (no rotation on the L = 4 and radix-2 stages, q = 0 skips, post-pass propagation gain 1, folded Q31 pre-shift, exact shift variances); every pinned test prints its measurement before it asserts; the inverse-exactness premise carries pre + 3 zero bits. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019ZPTzNxo5Fe4EtpXXKf7Sy --- tests/test_fft.cpp | 145 +++++++++---- tests/test_fft_fixed.cpp | 446 +++++++++++++++++++++++---------------- 2 files changed, 368 insertions(+), 223 deletions(-) diff --git a/tests/test_fft.cpp b/tests/test_fft.cpp index 0516bf6..cc92d85 100644 --- a/tests/test_fft.cpp +++ b/tests/test_fft.cpp @@ -16,8 +16,10 @@ // numerically what they were before the widening — same signals, same // tolerances, same comparisons; the trait is the identity for them. +#include #include #include +#include #include #include #include @@ -45,9 +47,13 @@ namespace { // Absolute tolerance in the sample's own units (fractions of full scale // for the fixed profiles). Float and double: the golden battery's - // numbers, unchanged. Q15 / Q31: measured against the closed forms below - // and pinned at 2x; see profile<> for the exponent handling. - // MEASURED: not yet (the kernel branch has not landed); 0.0 fails on purpose. + // numbers, unchanged. Q15 / Q31: the largest |got - expected| over the + // closed-form tests below (ImpulseHasFlatSpectrum, DcAndNyquistPacking, + // SignConventionIsPlusI) measured 2026-09-18 on x86-64 Linux, GCC + // 13.3.0 -O3, kernel a55a14f: 1.0 LSB for Q15 (the DC constant is + // 1 - 2^-15, compared against an ideal 1.0), 1.0 LSB for Q31 (the + // on-bin sine's rounding); pinned at 2 LSB each. See profile<> for the + // exponent handling. template constexpr double k_tolerance = 0.0; template <> @@ -55,9 +61,9 @@ namespace { template <> constexpr double k_tolerance = 2e-5; template <> - constexpr double k_tolerance = 0.0; + constexpr double k_tolerance = 2.0 / 32768.0; // 2 LSB of Q0.15 template <> - constexpr double k_tolerance = 0.0; + constexpr double k_tolerance = 2.0 / 2147483648.0; // 2 LSB of Q0.31 // ------------------------------------------------------------------------ // Per-profile trait. A forward transform's result is read as @@ -113,6 +119,30 @@ namespace { using sample_types = ::testing::Types; TYPED_TEST_SUITE(real_fft_test, sample_types); + /// EXPECT_NEAR that also tracks the largest |got - expected| seen, so a + /// test can print the measured number its pin was taken from (in the + /// sample's units; for the fixed profiles the pin is stated in LSBs). + class near_tracker { + public: + void check(double got, double expected, double tol, const char* what, size_t i) { + m_max = std::max(m_max, std::fabs(got - expected)); + EXPECT_NEAR(got, expected, tol) << what << " " << i; + } + double max() const { return m_max; } + + template + void report(const char* test) const { + if constexpr (profile::k_fixed_point) { + std::printf("[ measured ] %s %s: max |got - expected| = %.4f LSB (tolerance %.4f LSB)\n", + sizeof(Sample) == 2 ? "Q15" : "Q31", test, m_max / sample_scale::k_lsb, + k_tolerance / sample_scale::k_lsb); + } + } + + private: + double m_max = 0.0; + }; + constexpr int log2_of(size_t n) { int l = 0; while ((size_t{1} << l) < n) { @@ -167,9 +197,17 @@ namespace { /// the output LSB times the same power of two (under fixed scaling the /// trip discards log2 n bits twice: Q15 at n = 1024 reconstructs in /// steps of 2^-4, the honest number fft.h states), times the pinned - /// number of those units. MEASURED: not yet; 0.0 fails on purpose. + /// number of those units. Measured 2026-09-18 (x86-64 Linux, GCC 13.3.0 + /// -O3, kernel a55a14f) over RoundTripReproducesInput (n = 1024) and + /// RoundTripInPlaceAndAliased (n = 256): Q15 0.5054 / 0.5098 (the + /// output narrowing's half LSB), Q31 3.342 / 2.199 (the kernel's rounding + /// noise over two transforms); pinned at 2x the larger. template constexpr double k_round_trip_units = 0.0; + template <> + constexpr double k_round_trip_units = 1.02; + template <> + constexpr double k_round_trip_units = 6.7; template double round_trip_tolerance(size_t n) { @@ -188,17 +226,31 @@ namespace { EXPECT_EQ(fft.num_bins(), 513u); } + /// The measured round-trip error in reconstructed LSBs (fixed profiles). + template + void report_round_trip(const char* test, size_t n, double max_error) { + if constexpr (profile::k_fixed_point) { + const int e = profile::fft::fixed_scaling_exponent(n); + const double unit = sample_scale::k_lsb * std::ldexp(1.0, 2 * e + 1 - log2_of(n)); + std::printf("[ measured ] %s %s n=%zu: max error %.4f reconstructed LSB (unit %.3g; pin %.3f)\n", + sizeof(Sample) == 2 ? "Q15" : "Q31", test, n, max_error / unit, unit, + k_round_trip_units); + } + } + TYPED_TEST(real_fft_test, RoundTripReproducesInput) { constexpr size_t n = 1024; const auto x = random_signal(n, 42); const auto back = round_trip(n, x, false); const auto tol = round_trip_tolerance(n); - ASSERT_GT(tol, 0.0) << "unmeasured pin"; + near_tracker t; for (size_t i = 0; i < n; ++i) { - EXPECT_NEAR(back[i], profile::to_double(x[i]), tol) << "sample " << i; + t.check(back[i], profile::to_double(x[i]), tol, "sample", i); } + report_round_trip("RoundTripReproducesInput", n, t.max()); + EXPECT_GT(tol, 0.0) << "unmeasured pin"; } TYPED_TEST(real_fft_test, RoundTripInPlaceAndAliased) { @@ -207,11 +259,13 @@ namespace { const auto x = random_signal(n, 7); const auto back = round_trip(n, x, true); const auto tol = round_trip_tolerance(n); - ASSERT_GT(tol, 0.0) << "unmeasured pin"; + near_tracker t; for (size_t i = 0; i < n; ++i) { - EXPECT_NEAR(back[i], profile::to_double(x[i]), tol) << "sample " << i; + t.check(back[i], profile::to_double(x[i]), tol, "sample", i); } + report_round_trip("RoundTripInPlaceAndAliased", n, t.max()); + EXPECT_GT(tol, 0.0) << "unmeasured pin"; } TYPED_TEST(real_fft_test, ImpulseHasFlatSpectrum) { @@ -225,14 +279,16 @@ namespace { const int e = p::forward(fft, x.data()); const double scale = std::ldexp(1.0, -e); const double tol = k_tolerance; - ASSERT_GT(tol, 0.0) << "unmeasured pin"; - EXPECT_NEAR(p::to_double(x[0]), 1.0 * scale, tol); // DC - EXPECT_NEAR(p::to_double(x[1]), 1.0 * scale, tol); // Nyquist + near_tracker t; + t.check(p::to_double(x[0]), 1.0 * scale, tol, "dc", 0); + t.check(p::to_double(x[1]), 1.0 * scale, tol, "nyquist", 1); for (size_t k = 1; k < n / 2; ++k) { - EXPECT_NEAR(p::to_double(x[2 * k]), 1.0 * scale, tol) << "bin " << k << " real"; - EXPECT_NEAR(p::to_double(x[2 * k + 1]), 0.0, tol) << "bin " << k << " imag"; + t.check(p::to_double(x[2 * k]), 1.0 * scale, tol, "bin real", k); + t.check(p::to_double(x[2 * k + 1]), 0.0, tol, "bin imag", k); } + t.report("ImpulseHasFlatSpectrum"); + EXPECT_GT(tol, 0.0) << "unmeasured pin"; } TYPED_TEST(real_fft_test, DcAndNyquistPacking) { @@ -241,13 +297,13 @@ namespace { typename p::fft fft(n); const double tol = p::tolerance_at_n(n); - ASSERT_GT(tol, 0.0) << "unmeasured pin"; + near_tracker t; // Constant input: all energy in DC = data[0]. std::vector dc(n, p::from_double(1.0)); const int e_dc = p::forward(fft, dc.data()); - EXPECT_NEAR(p::to_double(dc[0]), static_cast(n) * std::ldexp(1.0, -e_dc), tol); - EXPECT_NEAR(p::to_double(dc[1]), 0.0, tol); + t.check(p::to_double(dc[0]), static_cast(n) * std::ldexp(1.0, -e_dc), tol, "dc slot", 0); + t.check(p::to_double(dc[1]), 0.0, tol, "dc slot", 1); // Alternating +1/-1: all energy in Nyquist = data[1]. std::vector nyq(n); @@ -255,8 +311,10 @@ namespace { nyq[i] = (i % 2 == 0) ? p::from_double(1.0) : p::from_double(-1.0); } const int e_nyq = p::forward(fft, nyq.data()); - EXPECT_NEAR(p::to_double(nyq[0]), 0.0, tol); - EXPECT_NEAR(p::to_double(nyq[1]), static_cast(n) * std::ldexp(1.0, -e_nyq), tol); + t.check(p::to_double(nyq[0]), 0.0, tol, "nyquist slot", 0); + t.check(p::to_double(nyq[1]), static_cast(n) * std::ldexp(1.0, -e_nyq), tol, "nyquist slot", 1); + t.report("DcAndNyquistPacking"); + EXPECT_GT(tol, 0.0) << "unmeasured pin"; } // The documented sign convention (W = exp(+2*pi*i/N)): a pure cosine at @@ -283,28 +341,35 @@ namespace { const double half_cos = static_cast(n) / 2.0 * std::ldexp(1.0, -e_cos); const double half_sin = static_cast(n) / 2.0 * std::ldexp(1.0, -e_sin); const double tol = p::tolerance_at_n(n); - ASSERT_GT(tol, 0.0) << "unmeasured pin"; - EXPECT_NEAR(p::to_double(cosine[2 * k]), half_cos, tol); - EXPECT_NEAR(p::to_double(cosine[2 * k + 1]), 0.0, tol); - EXPECT_NEAR(p::to_double(sine[2 * k]), 0.0, tol); - EXPECT_NEAR(p::to_double(sine[2 * k + 1]), half_sin, tol); // +N/2, not -N/2 + near_tracker t; + t.check(p::to_double(cosine[2 * k]), half_cos, tol, "cos re", k); + t.check(p::to_double(cosine[2 * k + 1]), 0.0, tol, "cos im", k); + t.check(p::to_double(sine[2 * k]), 0.0, tol, "sin re", k); + t.check(p::to_double(sine[2 * k + 1]), half_sin, tol, "sin im (+N/2, not -N/2)", k); // And nothing leaks into any other bin. for (size_t bin = 1; bin < n / 2; ++bin) { if (bin == k) { continue; } - EXPECT_NEAR(p::to_double(cosine[2 * bin]), 0.0, tol) << "cos leak, bin " << bin; - EXPECT_NEAR(p::to_double(sine[2 * bin]), 0.0, tol) << "sin leak, bin " << bin; + t.check(p::to_double(cosine[2 * bin]), 0.0, tol, "cos leak, bin", bin); + t.check(p::to_double(sine[2 * bin]), 0.0, tol, "sin leak, bin", bin); } + t.report("SignConventionIsPlusI"); + EXPECT_GT(tol, 0.0) << "unmeasured pin"; } // Parseval's relative tolerance for the fixed profiles: the rounding // noise adds power, so the relative error is of the order of the per-bin - // noise-to-signal ratio; measured and pinned at 2x. - // MEASURED: not yet; 0.0 fails on purpose. + // noise-to-signal ratio. Measured 2026-09-18 (x86-64 Linux, GCC 13.3.0 + // -O3, kernel a55a14f) at n = 512 on full-scale uniform noise: Q15 + // 7.15e-6, Q31 1.83e-9; pinned at 2x. template constexpr double k_parseval_relative = 0.0; + template <> + constexpr double k_parseval_relative = 1.43e-5; + template <> + constexpr double k_parseval_relative = 3.7e-9; TYPED_TEST(real_fft_test, ParsevalEnergyConservation) { using p = profile; @@ -331,7 +396,10 @@ namespace { freq_energy /= static_cast(n); if constexpr (p::k_fixed_point) { - ASSERT_GT(k_parseval_relative, 0.0) << "unmeasured pin"; + std::printf("[ measured ] %s ParsevalEnergyConservation: relative error %.3e (pin %.3e)\n", + sizeof(TypeParam) == 2 ? "Q15" : "Q31", std::fabs(freq_energy - time_energy) / time_energy, + k_parseval_relative); + EXPECT_GT(k_parseval_relative, 0.0) << "unmeasured pin"; EXPECT_NEAR(freq_energy, time_energy, k_parseval_relative * time_energy); } else { @@ -408,21 +476,24 @@ namespace { // 2-norm error of the fixed-scaling forward at n = 1024 on full-scale // uniform noise, a measured number pinned at 2x (the log_mel pattern). // Q15's error is the output narrowing (2^-15 / sqrt(12) per value on a - // spectrum whose per-component RMS is sqrt(1/3 / 2n)); Q31's is the - // int32 kernel's rounding noise, orders of magnitude lower. - // MEASURED: not yet; 0.0 fails on purpose. - constexpr double k_q15_tracks_double = 0.0; - constexpr double k_q31_tracks_double = 0.0; + // spectrum whose per-component RMS is sqrt(1/3 / 2n) = 0.0128: predicted + // 6.9e-4); Q31's is the int32 kernel's rounding noise, orders of + // magnitude lower. Measured 2026-09-18 (x86-64 Linux, GCC 13.3.0 -O3, + // kernel a55a14f): Q15 6.949e-4, Q31 5.458e-8. + constexpr double k_q15_tracks_double = 1.4e-3; + constexpr double k_q31_tracks_double = 1.1e-7; TEST(RealFftCrossPrecision, Q15TracksDouble) { - ASSERT_GT(k_q15_tracks_double, 0.0) << "unmeasured pin"; const double err = relative_error_vs_double(1024, 99); + std::printf("[ measured ] Q15TracksDouble: relative 2-norm error %.3e (pin %.3e)\n", err, k_q15_tracks_double); + EXPECT_GT(k_q15_tracks_double, 0.0) << "unmeasured pin"; EXPECT_LT(err, k_q15_tracks_double) << "measured " << err; } TEST(RealFftCrossPrecision, Q31TracksDouble) { - ASSERT_GT(k_q31_tracks_double, 0.0) << "unmeasured pin"; const double err = relative_error_vs_double(1024, 99); + std::printf("[ measured ] Q31TracksDouble: relative 2-norm error %.3e (pin %.3e)\n", err, k_q31_tracks_double); + EXPECT_GT(k_q31_tracks_double, 0.0) << "unmeasured pin"; EXPECT_LT(err, k_q31_tracks_double) << "measured " << err; } diff --git a/tests/test_fft_fixed.cpp b/tests/test_fft_fixed.cpp index e282558..f7fa1d5 100644 --- a/tests/test_fft_fixed.cpp +++ b/tests/test_fft_fixed.cpp @@ -33,7 +33,8 @@ // Measured numbers. Every tolerance and ratio here is a number measured on the // real kernel and pinned at 2x (the log_mel pattern), never a round number; // the table `pins` carries them all in one place with the host, compiler and -// date they were taken on. Fixed seeds, no wall clock, no filesystem, no +// date they were taken on, and every pinned test prints what it measured so +// a -V run records the current value beside the pin. Fixed seeds, no wall clock, no filesystem, no // : the battery runs unchanged on the four QEMU legs (Part 10), where // N <= 2048 fixed point is cheap and the double golden model is the expensive // part, so sizes are kept modest and TAP_DSP_PARITY_MAX_N caps the sweeps. @@ -114,8 +115,8 @@ namespace { // ------------------------------------------------------------------------ // THE PINS. Measured on the real kernel and pinned at 2x; see each test - // for what the number bounds. Host, compiler and date beside each block. - // A pin of 0.0 means "not yet measured" and fails the test on purpose. + // for what the number bounds. A pin of 0.0 on a field a row reads means + // "not measured" and fails that test on purpose. // ------------------------------------------------------------------------ struct pin_table { double saturation_max_lsb; ///< SaturationFreeWorstCaseDoesNotWrap: max |out - G/2^e| in output LSB @@ -128,12 +129,23 @@ namespace { double q15_vs_q31_lsb; ///< Q15AndQ31AgreeToTheQ15Floor: max |v15 - v31| in Q15 output LSB }; - // Not yet measured: the kernel branch has not landed. Every pin is 0.0 - // so the battery is red until the numbers are taken from the real kernel. - constexpr pin_table k_pins_q15_fixed{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - constexpr pin_table k_pins_q31_fixed{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - constexpr pin_table k_pins_q15_bfp{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - constexpr pin_table k_pins_q31_bfp{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + // Measured 2026-09-18 on x86-64 Linux (Ubuntu 24.04, glibc 2.39), GCC + // 13.3.0 and clang 18.1.3 -O3 (identical: the kernel is integer + // arithmetic and the tables are the same libm), kernel at + // claude/wave2-stage3b-kernel a55a14f, and pinned at 2x. The measured + // values, in the pins' order: + // Q15/fixed 0.500 / 0.562 / - / 1.00 / 0.0005 / 1.055 / 0.022 / 0.500 + // Q31/fixed 4.250 / 3.382 / - / 6.00 / 0.8125 / 1.744 / 0.383 / - + // Q15/bfp 0.750 / 5.000 / 1.0 / 1.00 / 0.0156 / 1.041 / 0.535 / 0.500 + // Q31/bfp 15.988 / 41.50 / 4.0 / 62.0 / 1.0469 / 2.321 / 0.644 / - + // The Q31 block-floating maxima all sit at index 0 or 1: the DC/Nyquist + // path, where round-half-up biases add coherently (fixed_point.h, + // "Honest limit"). bfp_vs_fixed_lsb is read by the block-floating rows + // only and q15_vs_q31_lsb by the Q15 rows only; the others carry 0.0. + constexpr pin_table k_pins_q15_fixed{1.0, 1.13, 0.0, 2.0, 0.001, 2.11, 0.044, 1.0}; + constexpr pin_table k_pins_q31_fixed{8.5, 6.77, 0.0, 12.0, 1.63, 3.49, 0.77, 0.0}; + constexpr pin_table k_pins_q15_bfp{1.5, 10.0, 2.0, 2.0, 0.032, 2.09, 1.07, 1.0}; + constexpr pin_table k_pins_q31_bfp{32.0, 83.0, 8.0, 124.0, 2.1, 4.65, 1.29, 0.0}; template constexpr const pin_table& pins() { @@ -390,8 +402,8 @@ namespace { TYPED_TEST_SUITE(fft_fixed_scaling_test, fixed_configs); template - class fft_bfp_test : public ::testing::Test {}; - TYPED_TEST_SUITE(fft_bfp_test, bfp_configs); + class fft_fixed_bfp_test : public ::testing::Test {}; + TYPED_TEST_SUITE(fft_fixed_bfp_test, bfp_configs); template class fft_fixed_point_test : public ::testing::Test {}; @@ -433,7 +445,7 @@ namespace { // scaling::block_floating: 0 <= e <= the fixed constant, in both // directions, whatever the input; a louder input never needs a smaller // exponent than silence does (0 is the floor). - TYPED_TEST(fft_bfp_test, BfpExponentIsWithinRange) { + TYPED_TEST(fft_fixed_bfp_test, BfpExponentIsWithinRange) { using cfg = TypeParam; using fft = typename cfg::fft; using s = typename cfg::sample; @@ -495,8 +507,10 @@ namespace { x[0] = static_cast(m * unit); const auto r = run_forward(x); ASSERT_EQ(r.exponent, e); + // Flat: DC, Nyquist and every real part m; every imaginary part 0. for (std::size_t i = 0; i < n; ++i) { - ASSERT_EQ(static_cast(r.out[i]), m) + const std::int64_t expected = (i < 2 || i % 2 == 0) ? m : 0; + ASSERT_EQ(static_cast(r.out[i]), expected) << "impulse A=" << m * unit << " n=" << n << " index " << i; } } @@ -543,7 +557,7 @@ namespace { imp[0] = std::int32_t{1} << 20; EXPECT_EQ(f.forward_inplace(imp.data()), 10); for (std::size_t i = 0; i < n; ++i) { - ASSERT_EQ(imp[i], std::int32_t{1} << 10) << i; + ASSERT_EQ(imp[i], (i < 2 || i % 2 == 0) ? std::int32_t{1} << 10 : 0) << i; } // The Q15 profile, for contrast, is X / N: the same 0.5 constant // comes back as 0.5. @@ -556,7 +570,7 @@ namespace { // The inverse under fixed scaling carries the same constant and is // exact on the same kind of input: the unnormalised inverse of a // DC-only spectrum a[0] = c is c/2 everywhere (fft.h's packing), so the - // output is c / 2^(e+1); of a flat spectrum (impulse response) it is + // output is c / 2^(e+1); of a flat spectrum (the impulse response) it is // N c / 2 at k = 0 and 0 elsewhere, so out_0 = c / 2^(pre + 1). TYPED_TEST(fft_fixed_scaling_test, FixedInverseCarriesTheSameExponent) { using cfg = TypeParam; @@ -577,8 +591,13 @@ namespace { ASSERT_EQ(static_cast(r.out[i]), m) << "dc spectrum n=" << n << " index " << i; } } - // The flat spectrum: a[0] = a[1] = a[2k] = c, a[2k+1] = 0. - const std::int64_t step = std::int64_t{1} << (cfg::k_pre + 2); + // The flat spectrum: a[0] = a[1] = a[2k] = c, a[2k+1] = 0. The + // pre-pass takes pre + 1 bits, then the constant c / 2^(pre+1) + // rides the unrotated path through every stage's 2-bit shift, so + // c must carry pre + 3 zero bits for the trip to be exact (with + // fewer, round-half-up doubles a 2-LSB constant at every stage: + // that is the kernel's honest bias, not a scale error). + const std::int64_t step = std::int64_t{1} << (cfg::k_pre + 3); for (const std::int64_t c : {step, -step, (static_cast(cfg::k_max) / step) * step}) { std::vector a(n, s{0}); a[0] = static_cast(c); @@ -607,48 +626,60 @@ namespace { // the pinned number of LSBs (a wrap is 2^15 / 2^31 LSBs off, a clamp at // least the amount clamped), and no output sample sits on a rail unless // the golden model puts it within the pin of that rail. + struct worst_case { + double value = 0.0; + const char* what = ""; + const char* name = ""; + std::size_t n = 0; + int e = 0; + std::size_t index = 0; + void note(double v, const char* w, const char* nm, std::size_t size, int exponent, std::size_t i) { + if (v > value) { + value = v; + what = w; + name = nm; + n = size; + e = exponent; + index = i; + } + } + }; + TYPED_TEST(fft_fixed_point_test, SaturationFreeWorstCaseDoesNotWrap) { using cfg = TypeParam; using s = typename cfg::sample; const auto pin = pins().saturation_max_lsb; - ASSERT_GT(pin, 0.0) << "unmeasured pin"; - double worst = 0.0; + worst_case worst; + worst_case rail; // largest shortfall of the golden model below a rail the output sits on for (const std::size_t n : k_sweep_sizes) { - for (const auto& p : adversarial_time_patterns(n)) { - const auto r = run_forward(p.x); - const auto golden = golden_forward(fractions(p.x)); - const auto d = deviation_from_golden(r.out, golden, r.exponent); - worst = std::max(worst, d.max_lsb); - EXPECT_LE(d.max_lsb, pin) << cfg::name() << " forward " << p.name << " n=" << n << " e=" << r.exponent - << " at index " << d.argmax; - for (std::size_t i = 0; i < n; ++i) { - if (r.out[i] == cfg::k_max || r.out[i] == cfg::k_min) { - const double g = std::fabs(golden[i]) * pow2(-r.exponent) / cfg::k_lsb; - EXPECT_GE(g, static_cast(cfg::k_max) - pin) - << cfg::name() << " forward " << p.name << " n=" << n << ": sample " << i - << " sits on a rail the golden model does not reach"; - } - } - } - for (const auto& p : adversarial_spectrum_patterns(n)) { - const auto r = run_inverse(p.x); - const auto golden = golden_inverse(fractions(p.x)); - const auto d = deviation_from_golden(r.out, golden, r.exponent); - worst = std::max(worst, d.max_lsb); - EXPECT_LE(d.max_lsb, pin) << cfg::name() << " inverse " << p.name << " n=" << n << " e=" << r.exponent - << " at index " << d.argmax; - for (std::size_t i = 0; i < n; ++i) { - if (r.out[i] == cfg::k_max || r.out[i] == cfg::k_min) { - const double g = std::fabs(golden[i]) * pow2(-r.exponent) / cfg::k_lsb; - EXPECT_GE(g, static_cast(cfg::k_max) - pin) - << cfg::name() << " inverse " << p.name << " n=" << n << ": sample " << i - << " sits on a rail the golden model does not reach"; + for (const bool inverse : {false, true}) { + const auto patterns = inverse ? adversarial_spectrum_patterns(n) : adversarial_time_patterns(n); + for (const auto& p : patterns) { + const auto r = inverse ? run_inverse(p.x) : run_forward(p.x); + const auto golden = inverse ? golden_inverse(fractions(p.x)) : golden_forward(fractions(p.x)); + const auto d = deviation_from_golden(r.out, golden, r.exponent); + worst.note(d.max_lsb, inverse ? "inverse" : "forward", p.name, n, r.exponent, d.argmax); + for (std::size_t i = 0; i < n; ++i) { + if (r.out[i] == cfg::k_max || r.out[i] == cfg::k_min) { + const double g = std::fabs(golden[i]) * pow2(-r.exponent) / cfg::k_lsb; + rail.note(static_cast(cfg::k_max) - g, inverse ? "inverse" : "forward", p.name, n, + r.exponent, i); + } } } } } - std::printf("[ measured ] %s saturation sweep: max |out - G/2^e| = %.3f LSB (pin %.3f)\n", cfg::name(), worst, - pin); + std::printf("[ measured ] %s saturation sweep: max |out - G/2^e| = %.3f LSB (pin %.3f) at %s %s n=%zu e=%d " + "index %zu; largest rail shortfall %.3f LSB\n", + cfg::name(), worst.value, pin, worst.what, worst.name, worst.n, worst.e, worst.index, rail.value); + EXPECT_GT(pin, 0.0) << "unmeasured pin"; + EXPECT_LE(worst.value, pin) << cfg::name() << " " << worst.what << " " << worst.name << " n=" << worst.n + << " e=" << worst.e << " index " << worst.index; + // No output sits on a rail unless the golden model is within the pin + // of that rail: a clamp would show as a shortfall of at least what was + // clamped, a wrap as 2^15 / 2^31 LSB of deviation above. + EXPECT_LE(rail.value, pin) << cfg::name() << " " << rail.what << " " << rail.name << " n=" << rail.n + << ": sample " << rail.index << " sits on a rail the golden model does not reach"; } // A silent block is silent out, in both directions. @@ -683,29 +714,31 @@ namespace { using cfg = TypeParam; using s = typename cfg::sample; const auto pin = pins().round_trip_k; - ASSERT_GT(pin, 0.0) << "unmeasured pin"; - double worst = 0.0; + worst_case worst; for (const std::size_t n : {std::size_t{16}, std::size_t{64}, std::size_t{256}, std::size_t{1024}}) { for (const double amplitude : {1.0, 0.25, 0.01}) { const auto x = tap::dsp::test::random_signal(n, 0xC0FFEE01u ^ static_cast(n), amplitude); - const auto xf = fractions(x); - const auto fwd = run_forward(x); - const auto inv = run_inverse(fwd.out); - const int shift = fwd.exponent + inv.exponent + 1 - log2_size(n); - const double unit = cfg::k_lsb * pow2(shift); - double max_err = 0.0; + const auto xf = fractions(x); + const auto fwd = run_forward(x); + const auto inv = run_inverse(fwd.out); + const int shift = fwd.exponent + inv.exponent + 1 - log2_size(n); + const double unit = cfg::k_lsb * pow2(shift); for (std::size_t i = 0; i < n; ++i) { const double back = sample_scale::to_double(inv.out[i]) * pow2(shift); - max_err = std::max(max_err, std::fabs(back - xf[i]) / unit); + worst.note(std::fabs(back - xf[i]) / unit, + amplitude == 1.0 ? "0 dBFS" : (amplitude == 0.25 ? "-12 dBFS" : "-40 dBFS"), + "round trip", n, fwd.exponent * 100 + inv.exponent, i); } - worst = std::max(worst, max_err); - EXPECT_LE(max_err, pin) << cfg::name() << " n=" << n << " amplitude " << amplitude - << " e_fwd=" << fwd.exponent << " e_inv=" << inv.exponent << " unit=" << unit; } } - std::printf("[ measured ] %s round trip: max error %.3f reconstructed LSB (pin %.3f)\n", cfg::name(), worst, - pin); + std::printf("[ measured ] %s round trip: max error %.3f reconstructed LSB (pin %.3f) at %s n=%zu " + "e_fwd=%d e_inv=%d index %zu\n", + cfg::name(), worst.value, pin, worst.what, worst.n, worst.e / 100, worst.e % 100, worst.index); + EXPECT_GT(pin, 0.0) << "unmeasured pin"; + EXPECT_LE(worst.value, pin) << cfg::name() << " " << worst.what << " n=" << worst.n + << " e_fwd=" << worst.e / 100 << " e_inv=" << worst.e % 100 << " index " + << worst.index; } // forward()/inverse() are copy-then-in-place: same bits, same exponent, @@ -752,13 +785,12 @@ namespace { // difference is the fixed path's extra rounding noise plus one rounding // of the shift); and on an input where BFP reports the full constant it // shifted like fixed at every stage, so the two are bit-identical. - TYPED_TEST(fft_bfp_test, BfpMatchesFixedAfterShift) { + TYPED_TEST(fft_fixed_bfp_test, BfpMatchesFixedAfterShift) { using cfg = TypeParam; using s = typename cfg::sample; using fixed = config; const auto pin = pins().bfp_vs_fixed_lsb; - ASSERT_GT(pin, 0.0) << "unmeasured pin"; - double worst = 0.0; + worst_case worst; for (const std::size_t n : k_sweep_sizes) { std::vector> inputs; inputs.push_back({"noise 0 dBFS", tap::dsp::test::random_signal(n, 0x2545F491u, 1.0)}); @@ -767,25 +799,33 @@ namespace { inputs.push_back({"tone -6 dBFS", tap::dsp::test::tone(n, static_cast(n / 8 + 1), 0.5, 0.3)}); inputs.push_back({"dc +full", std::vector(n, cfg::k_max)}); for (const auto& p : inputs) { - const auto f = run_forward(p.x); - const auto b = run_forward(p.x); - ASSERT_LE(b.exponent, f.exponent) << p.name << " n=" << n; - const int shift = f.exponent - b.exponent; - for (std::size_t i = 0; i < n; ++i) { - const std::int64_t bv = static_cast(b.out[i]); - const std::int64_t shifted = - shift == 0 ? bv : ((bv + (std::int64_t{1} << (shift - 1))) >> shift); // round-half-up - const double diff = std::fabs(static_cast(shifted - static_cast(f.out[i]))); - worst = std::max(worst, diff); - EXPECT_LE(diff, pin) << cfg::name() << " " << p.name << " n=" << n << " e_fixed=" << f.exponent - << " e_bfp=" << b.exponent << " index " << i; + for (const bool inverse : {false, true}) { + const auto f = inverse ? run_inverse(p.x) : run_forward(p.x); + const auto b = inverse ? run_inverse(p.x) : run_forward(p.x); + ASSERT_LE(b.exponent, f.exponent) << p.name << " n=" << n; + const int shift = f.exponent - b.exponent; + for (std::size_t i = 0; i < n; ++i) { + const std::int64_t bv = static_cast(b.out[i]); + const std::int64_t shifted = + shift == 0 ? bv : ((bv + (std::int64_t{1} << (shift - 1))) >> shift); + const double diff = + std::fabs(static_cast(shifted - static_cast(f.out[i]))); + worst.note(diff, inverse ? "inverse" : "forward", p.name, n, f.exponent * 100 + b.exponent, i); + } } } } - std::printf("[ measured ] %s vs fixed after shift: max %.1f LSB (pin %.1f)\n", cfg::name(), worst, pin); + std::printf("[ measured ] %s vs fixed after shift: max %.1f LSB (pin %.1f) at %s %s n=%zu e_fixed=%d e_bfp=%d " + "index %zu\n", + cfg::name(), worst.value, pin, worst.what, worst.name, worst.n, worst.e / 100, worst.e % 100, + worst.index); + EXPECT_GT(pin, 0.0) << "unmeasured pin"; + EXPECT_LE(worst.value, pin) << cfg::name() << " " << worst.what << " " << worst.name << " n=" << worst.n + << " e_fixed=" << worst.e / 100 << " e_bfp=" << worst.e % 100 << " index " + << worst.index; } - TYPED_TEST(fft_bfp_test, BfpAtTheFullExponentIsBitIdenticalToFixed) { + TYPED_TEST(fft_fixed_bfp_test, BfpAtTheFullExponentIsBitIdenticalToFixed) { using cfg = TypeParam; using s = typename cfg::sample; using fixed = config; @@ -826,10 +866,8 @@ namespace { using s = typename cfg::sample; const auto max_pin = pins().negation_sum_max_lsb; const auto bias_pin = pins().negation_bias_lsb; - ASSERT_GT(max_pin, 0.0) << "unmeasured pin"; - ASSERT_GT(bias_pin, 0.0) << "unmeasured pin"; - double worst_max = 0.0; - double worst_bias = 0.0; + worst_case worst_max; + worst_case worst_bias; for (const std::size_t n : k_sweep_sizes) { for (const double amplitude : {0.999, 0.1}) { // 0.999 keeps -x representable (-INT_MIN would saturate). @@ -843,25 +881,32 @@ namespace { const auto a = inverse ? run_inverse(x) : run_forward(x); const auto b = inverse ? run_inverse(neg) : run_forward(neg); ASSERT_EQ(a.exponent, b.exponent) << "negation changed the exponent, n=" << n; - double max_sum = 0.0; - double mean = 0.0; + double mean = 0.0; for (std::size_t i = 0; i < n; ++i) { const double sum = static_cast(a.out[i]) + static_cast(b.out[i]); - max_sum = std::max(max_sum, std::fabs(sum)); + worst_max.note(std::fabs(sum), inverse ? "inverse" : "forward", + amplitude > 0.5 ? "-0 dBFS" : "-20 dBFS", n, a.exponent, i); mean += sum; } - mean = std::fabs(mean / static_cast(n)); - worst_max = std::max(worst_max, max_sum); - worst_bias = std::max(worst_bias, mean); - EXPECT_LE(max_sum, max_pin) << cfg::name() << (inverse ? " inverse" : " forward") << " n=" << n - << " amplitude " << amplitude; - EXPECT_LE(mean, bias_pin) << cfg::name() << (inverse ? " inverse" : " forward") << " n=" << n - << " amplitude " << amplitude; + // The bias is a mean over the block: read it where the + // block is large enough for a mean to say something. + if (n >= 64) { + worst_bias.note(std::fabs(mean / static_cast(n)), inverse ? "inverse" : "forward", + amplitude > 0.5 ? "-0 dBFS" : "-20 dBFS", n, a.exponent, 0); + } } } } - std::printf("[ measured ] %s F(x)+F(-x): max %.2f LSB (pin %.2f), bias %.4f LSB (pin %.4f)\n", cfg::name(), - worst_max, max_pin, worst_bias, bias_pin); + std::printf("[ measured ] %s F(x)+F(-x): max %.2f LSB (pin %.2f) at %s %s n=%zu index %zu; bias %.4f LSB " + "(pin %.4f) at %s %s n=%zu\n", + cfg::name(), worst_max.value, max_pin, worst_max.what, worst_max.name, worst_max.n, worst_max.index, + worst_bias.value, bias_pin, worst_bias.what, worst_bias.name, worst_bias.n); + EXPECT_GT(max_pin, 0.0) << "unmeasured pin"; + EXPECT_GT(bias_pin, 0.0) << "unmeasured pin"; + EXPECT_LE(worst_max.value, max_pin) << cfg::name() << " " << worst_max.what << " " << worst_max.name + << " n=" << worst_max.n << " index " << worst_max.index; + EXPECT_LE(worst_bias.value, bias_pin) + << cfg::name() << " " << worst_bias.what << " " << worst_bias.name << " n=" << worst_bias.n; } // ======================================================================== @@ -873,10 +918,10 @@ namespace { // spectra differ by the Q15 output quantisation plus the Q15 profile's // (coarser) internal noise, pinned in Q15 output LSBs. template - double q15_vs_q31(const char* what, double pin) { - using c15 = config; - using c31 = config; - double worst = 0.0; + void q15_vs_q31(const char* what, double pin) { + using c15 = config; + using c31 = config; + worst_case worst; for (const std::size_t n : k_sweep_sizes) { for (const double amplitude : {1.0, 0.05}) { const auto x15 = tap::dsp::test::random_signal( @@ -885,27 +930,29 @@ namespace { for (std::size_t i = 0; i < n; ++i) { x31[i] = static_cast(x15[i]) << 16; // exact } - const auto r15 = run_forward(x15); - const auto r31 = run_forward(x31); - // Compare in the unnormalised frame, in Q15 output LSBs at e15. - const double unit = c15::k_lsb * pow2(r15.exponent); - for (std::size_t i = 0; i < n; ++i) { - const double v15 = sample_scale::to_double(r15.out[i]) * pow2(r15.exponent); - const double v31 = sample_scale::to_double(r31.out[i]) * pow2(r31.exponent); - const double diff = std::fabs(v15 - v31) / unit; - worst = std::max(worst, diff); - EXPECT_LE(diff, pin) << what << " n=" << n << " amplitude " << amplitude << " index " << i - << " e15=" << r15.exponent << " e31=" << r31.exponent; + for (const bool inverse : {false, true}) { + const auto r15 = inverse ? run_inverse(x15) : run_forward(x15); + const auto r31 = inverse ? run_inverse(x31) : run_forward(x31); + // Compare in the unnormalised frame, in Q15 output LSBs at e15. + const double unit = c15::k_lsb * pow2(r15.exponent); + for (std::size_t i = 0; i < n; ++i) { + const double v15 = sample_scale::to_double(r15.out[i]) * pow2(r15.exponent); + const double v31 = sample_scale::to_double(r31.out[i]) * pow2(r31.exponent); + worst.note(std::fabs(v15 - v31) / unit, inverse ? "inverse" : "forward", + amplitude == 1.0 ? "0 dBFS" : "-26 dBFS", n, r15.exponent * 100 + r31.exponent, i); + } } } } - std::printf("[ measured ] Q15 vs Q31 (%s): max %.3f Q15 LSB (pin %.3f)\n", what, worst, pin); - return worst; + std::printf( + "[ measured ] Q15 vs Q31 (%s): max %.3f Q15 LSB (pin %.3f) at %s %s n=%zu e15=%d e31=%d index %zu\n", what, + worst.value, pin, worst.what, worst.name, worst.n, worst.e / 100, worst.e % 100, worst.index); + EXPECT_GT(pin, 0.0) << "unmeasured pin"; + EXPECT_LE(worst.value, pin) << what << " " << worst.what << " " << worst.name << " n=" << worst.n + << " e15=" << worst.e / 100 << " e31=" << worst.e % 100 << " index " << worst.index; } TEST(fft_fixed_profiles, Q15AndQ31AgreeToTheQ15Floor) { - ASSERT_GT(k_pins_q15_fixed.q15_vs_q31_lsb, 0.0) << "unmeasured pin"; - ASSERT_GT(k_pins_q15_bfp.q15_vs_q31_lsb, 0.0) << "unmeasured pin"; q15_vs_q31("fixed", k_pins_q15_fixed.q15_vs_q31_lsb); q15_vs_q31("bfp", k_pins_q15_bfp.q15_vs_q31_lsb); } @@ -916,72 +963,96 @@ namespace { // // THE MODEL (Welch 1969, "A fixed-point fast Fourier transform error // analysis"; Oppenheim & Weinstein 1972), specialised to this kernel's - // arithmetic as fft_arith.h states it: + // arithmetic as fft_arith.h states it and to the structure fixed_point.h + // documents (the model is written here, independently, from those two + // headers; the kernel's design note carries its own derivation): // // - Internal LSB q: 2^-31 for Q31, 2^-29 for the widened Q15 (Q2.29), - // both as fractions of full scale. Every rounding (shr_round, - // mul_coeff) is additive noise uniform in (-q/2, q/2], variance - // q^2/12, independent of every other rounding. - // - Structure (Part 7): a radix-4 DIF complex FFT of length M = N/2 - // (floor(log2 M / 2) radix-4 stages, one radix-2 stage when log2 M is - // odd), then Ooura's real post-pass, which pairs bins k and M-k and - // applies one complex product per pair. Under fixed scaling the Q31 - // profile pre-shifts its input by one bit (one rounding, no sum). + // both as fractions of full scale. Every rounding is additive noise, + // independent of every other: a mul_coeff rounding is uniform on a + // continuous interval, variance q^2/12; an s-bit shr_round of an + // integer takes one of 2^s equally likely residues, variance + // (1 - 4^-s) q^2/12 (3/4 of q^2/12 for one bit, 15/16 for two). + // - Structure: a radix-4 DIF complex FFT of length M = N/2 (one stage + // per factor of 4, spans M, M/4, ..., down to 4; a radix-2 stage when + // log2 M is odd), then Ooura's real post-pass pairing bins k and + // M - k with one complex product per pair (bins 1 .. N/4 - 1 and + // their mirrors: every interior bin but N/4). Under fixed scaling the + // Q31 profile folds its one-bit input pre-shift into the first + // stage's shift (a single 3-bit rounding). // - Per stage, per output component. Shift-before-butterfly: each of - // the stage's sum_gain inputs (4, 2, 2 for radix-4, radix-2, the - // post-pass) is rounded once per component and reaches every output - // component through the butterfly with unit weight, so the shift - // injects sum_gain * q^2/12 -- when the stage shifts at all (a BFP - // stage shifting by 0 rounds nothing: shr_round(x, 0) is exact). The - // twiddle product is the two-rounding complex multiply: 2 q^2/12 per - // component. The twiddle's own quantisation, |dw| <= 2^-31 per - // component (0.5 LSB of Q1.30), uniform, contributes 2 * P * 2^-62 / 3 - // where P is the signal variance per component entering the product. - // - Propagation. Noise present after stage t reaches the output through - // each later stage u with gain sum_gain_u * 4^-shift_u: a sum of + // the stage's sum_gain inputs (4 for radix-4, 2 for radix-2) is + // rounded once per component and reaches every output component with + // unit weight, so the shift injects sum_gain * v(s) -- when the stage + // shifts at all (shr_round(x, 0) is exact, so a BFP stage shifting by + // 0 injects nothing). The twiddle rotation is the two-rounding + // complex multiply, 2 q^2/12 per component, on 3 of the 4 outputs of + // every butterfly whose index q is not 0: averaged over a stage of + // span L that is 2 * (3/4) * (1 - 4/L) q^2/12 per component, zero for + // the last radix-4 stage (L = 4) and for the radix-2 stage. The + // twiddle's own quantisation, |dw| <= 2^-31 per component (half an + // LSB of Q1.30), uniform, adds 2 * P * 2^-62 / 3 on the rotated + // outputs, P being the signal variance per component entering the + // rotation. + // - Post-pass: its one-bit shift injects v(1) (the two paired values' + // roundings reach the output with weights |1 - w|^2 + |w|^2 = 1), + // its rotation 2 q^2/12; it carries earlier noise and the signal with + // gain 1 * 4^-shift. + // - Propagation: noise present after a stage reaches the output through + // each later stage u with power gain sum_gain_u * 4^-shift_u (a sum of // sum_gain_u uncorrelated terms through unit-magnitude twiddles, then - // the shift. Signal: a white input of variance s^2 per sample enters - // as M complex values of variance s^2 per component and follows the - // same gain, without the injections. - // - The Q15 narrowing adds (2^-15)^2 / 12 per output value. + // the shift). Signal: a white input of variance sigma^2 per sample + // enters as M complex values of variance sigma^2 per component and + // follows the same gains without the injections; through the + // post-pass it lands at sigma^2 / (2N) * 4^(log2 N - e), the DFT's + // N sigma^2 / 2 per component at exponent e. + // - The Q15 narrowing adds (2^-15)^2 / 12 per output value; under block + // floating point one more shift (growth 1) precedes it. // - Shift schedule. scaling::fixed: every stage shifts its full amount - // (2, 1, 1; plus the Q31 pre-shift) and the sum is the fixed exponent. - // scaling::block_floating: the kernel returns only the total e, and - // the schedule is not observable, so the model takes the EARLY - // schedule (e assigned to the first stages, each up to its fixed - // amount), which is the worst case for noise: a bit shifted early - // attenuates nothing injected after it. The late schedule (what an - // on-demand scaler does for a quiet input) is the best case, and the - // table prints both so the measured floor can be read against them. + // and the sum is fixed_scaling_exponent(N). scaling::block_floating: + // the kernel returns the total e and the per-stage schedule is not + // observable, so the model brackets it: the EARLY schedule (e + // assigned to the first stages, each up to its fixed amount) is the + // worst case for noise, since a bit shifted early attenuates nothing + // injected after it, and is what the pin is measured against; the + // LATE schedule is the best case and the table prints it beside. // // The model's number is the noise variance per output component in the // output's own units (fractions of full scale at exponent e). The // measurement is the mean over the interior bins (indices 2 .. N-1; DC - // and Nyquist take a different post-pass path) of |out - G(x)/2^e|^2 per - // component, with G the double golden model on the SAME quantised input, - // so only the kernel's roundings and its twiddle quantisation are in it. - // The pin is the largest measured/model ratio over the sizes and levels, - // per configuration and material, at 2x; the table below is printed for - // a -V run to record (Part 13). + // and Nyquist take the glue path) of |out - G(x)/2^e|^2 per component, + // with G the double golden model on the SAME quantised input, so only + // the kernel's roundings and its twiddle quantisation are in it. What + // the variance model leaves out is the round-half-up bias: an s-bit + // shr_round has mean +2^-(s+1) LSB, which the final shifts put on every + // bin and the earlier ones spread with rotated phases; it shows as a + // measured/model ratio above one and is pinned as such. The pin is the + // largest measured/model ratio over the sizes and levels, per + // configuration and material, at 2x; the table is printed for a -V run + // to record (Part 13). struct stage_spec { - int sum_gain; ///< inputs summed into each output component - int max_shift; ///< the fixed-scaling shift, and the BFP maximum - int product_roundings; ///< roundings per output component in the twiddle product + int sum_gain; ///< inputs summed into each output component + int max_shift; ///< the fixed-scaling shift, and the BFP maximum + double rotation_fraction; ///< fraction of output components that see a twiddle rotation }; - std::vector kernel_stages(std::size_t n, bool with_pre_shift) { + std::vector kernel_stages(std::size_t n, bool fold_pre_shift, bool q15_final_shift) { std::vector stages; - if (with_pre_shift) { - stages.push_back({1, 1, 0}); + const std::size_t m = n / 2; + for (std::size_t span = m; span >= 4; span /= 4) { + const double rotated = 0.75 * (1.0 - 4.0 / static_cast(span)); // 3 of 4 outputs, q != 0 + stages.push_back({4, 2, rotated}); } - const int log2_m = log2_size(n / 2); - for (int s = 0; s < log2_m / 2; ++s) { - stages.push_back({4, 2, 2}); + if (stages.empty() || (m >> (2 * stages.size())) == 2) { + stages.push_back({2, 1, 0.0}); // the radix-2 final stage (log2 M odd) } - if (log2_m % 2 == 1) { - stages.push_back({2, 1, 2}); + if (fold_pre_shift) { + stages.front().max_shift += 1; // Q31 fixed: pre-shift and first stage, one rounding + } + stages.push_back({1, 1, 1.0}); // the real post-pass + if (q15_final_shift) { + stages.push_back({1, 31, 0.0}); // Q15 BFP: the shift before the narrow, unbounded } - stages.push_back({2, 1, 2}); // the real post-pass return stages; } @@ -1018,7 +1089,7 @@ namespace { welch_prediction welch_model(const std::vector& stages, const std::vector& shifts, double q_int, double q_out_extra, double input_variance) { - constexpr double twiddle_lsb = 0x1p-31; // 0.5 LSB of Q1.30 + constexpr double twiddle_lsb = 0x1p-31; // half an LSB of Q1.30 const double rounding = q_int * q_int / 12.0; double noise = 0.0; double signal = input_variance; @@ -1029,11 +1100,10 @@ namespace { noise = noise * stage_gain; signal = signal * stage_gain; if (shifts[t] > 0) { - noise += static_cast(st.sum_gain) * rounding; + noise += static_cast(st.sum_gain) * (1.0 - pow2(-2 * shifts[t])) * rounding; } - if (st.product_roundings > 0) { - noise += static_cast(st.product_roundings) * rounding; - noise += 2.0 * signal * twiddle_lsb * twiddle_lsb / 3.0; + if (st.rotation_fraction > 0.0) { + noise += st.rotation_fraction * (2.0 * rounding + 2.0 * signal * twiddle_lsb * twiddle_lsb / 3.0); } } noise += q_out_extra * q_out_extra / 12.0; @@ -1073,7 +1143,7 @@ namespace { const double q_int = pow2(-Cfg::k_int_frac); const double q_out_extra = Cfg::k_is_q15 ? Cfg::k_lsb : 0.0; - const auto stages = kernel_stages(n, !Cfg::k_is_bfp && Cfg::k_pre > 0); + const auto stages = kernel_stages(n, !Cfg::k_is_bfp && Cfg::k_pre > 0, Cfg::k_is_bfp && Cfg::k_is_q15); const auto model_shifts = shifts_for(stages, r.exponent, Cfg::k_is_bfp ? schedule::early : schedule::fixed); const auto late_shifts = shifts_for(stages, r.exponent, Cfg::k_is_bfp ? schedule::late : schedule::fixed); const auto model = welch_model(stages, model_shifts, q_int, q_out_extra, input_variance); @@ -1082,7 +1152,7 @@ namespace { noise_row row{material, n, level_db, r.exponent, db(signal), db(noise), db(model.noise), db(late.noise), noise / model.noise}; - std::printf("[ floor ] %s %-6s N=%5zu %4.0f dBFS e=%2d signal %7.2f floor %7.2f model %7.2f (late %7.2f) " + std::printf("[ floor ] %s %-5s N=%5zu %4.0f dBFS e=%2d signal %7.2f floor %7.2f model %7.2f (late %7.2f) " "dBFS/component ratio %.3f snr %6.2f dB\n", Cfg::name(), material, n, level_db, r.exponent, row.signal_db, row.floor_db, row.model_db, row.model_late_db, row.ratio, row.signal_db - row.floor_db); @@ -1090,14 +1160,12 @@ namespace { } TYPED_TEST(fft_fixed_point_test, NoiseFloorTracksWelchModel) { - using cfg = TypeParam; - using s = typename cfg::sample; - const auto noise_pin = pins().noise_ratio_noise; - const auto tone_pin = pins().noise_ratio_tone; - ASSERT_GT(noise_pin, 0.0) << "unmeasured pin"; - ASSERT_GT(tone_pin, 0.0) << "unmeasured pin"; - double worst_noise = 0.0; - double worst_tone = 0.0; + using cfg = TypeParam; + using s = typename cfg::sample; + const auto noise_pin = pins().noise_ratio_noise; + const auto tone_pin = pins().noise_ratio_tone; + double worst_noise = 0.0; + double worst_tone = 0.0; for (const std::size_t n : {std::size_t{256}, std::size_t{512}, std::size_t{2048}}) { for (const double level_db : {0.0, -20.0, -40.0, -60.0}) { const double amplitude = std::pow(10.0, level_db / 20.0); @@ -1113,14 +1181,16 @@ namespace { worst_tone = std::max(worst_tone, trow.ratio); } } + std::printf("[ measured ] %s noise/model ratio: noise %.3f (pin %.3f), tone %.3f (pin %.3f)\n", cfg::name(), + worst_noise, noise_pin, worst_tone, tone_pin); + EXPECT_GT(noise_pin, 0.0) << "unmeasured pin"; + EXPECT_GT(tone_pin, 0.0) << "unmeasured pin"; EXPECT_LE(worst_noise, noise_pin) << cfg::name() << " white-noise floor above the pinned ratio to the model"; EXPECT_LE(worst_tone, tone_pin) << cfg::name() << " on-bin tone floor above the pinned ratio to the model"; // A floor far BELOW the model would mean a different arithmetic (a // fused complex multiply, an exact shift) or a broken measurement, // not a better kernel: the two-rounding form is the contract. EXPECT_GE(worst_noise, noise_pin / 8.0) << cfg::name() << " white-noise floor implausibly far below the model"; - std::printf("[ measured ] %s noise/model ratio: noise %.3f (pin %.3f), tone %.3f (pin %.3f)\n", cfg::name(), - worst_noise, noise_pin, worst_tone, tone_pin); } // ======================================================================== @@ -1216,9 +1286,13 @@ namespace { std::uint64_t twiddles; ///< make_twiddle_table(n / 2) std::uint64_t post_pass; ///< make_real_post_pass_table(n) }; - // Not yet measured: filled from the real kernel branch; the host that - // produced each value is recorded here. - constexpr std::array expected{{{256, 0, 0}, {512, 0, 0}, {2048, 0, 0}}}; + // Taken 2026-09-18 on x86-64 Linux, glibc 2.39 (GCC 13.3.0 and clang + // 18.1.3 agree), kernel a55a14f. The other CI hosts (macOS arm64, + // Windows UCRT, the newlib QEMU legs) either reproduce these or the + // difference is a recorded finding per host and N. + constexpr std::array expected{{{256, 0x95f5c68afe494835ull, 0x66c84a75861eafb6ull}, + {512, 0x6df6ff99a3ed3c85ull, 0x4b1374200abed27cull}, + {2048, 0xe42c528f3ae88b45ull, 0xa0f40e80bbf4efb9ull}}}; for (const auto& p : expected) { const auto twiddles = tap::dsp::detail::make_twiddle_table(p.n / 2); const auto post = tap::dsp::detail::make_real_post_pass_table(p.n); From b1b165ddba813b4c5562097e02360c94a94503b4 Mon Sep 17 00:00:00 2001 From: Timothy Place Date: Fri, 18 Sep 2026 01:52:00 +0000 Subject: [PATCH 4/4] tests: newlib-safe printf formats in the fixed-point battery newlib's printf on the QEMU legs has no %zu: it printed "zu", left the varargs desynchronised, and the next %s HardFaulted the M55 run inside RoundingBiasOnNegatedInputIsBounded. Sizes now go through %lu as unsigned long and the 64-bit table checksums print as two %08lx halves. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019ZPTzNxo5Fe4EtpXXKf7Sy --- tests/test_fft.cpp | 7 ++--- tests/test_fft_fixed.cpp | 58 ++++++++++++++++++++++++---------------- 2 files changed, 39 insertions(+), 26 deletions(-) diff --git a/tests/test_fft.cpp b/tests/test_fft.cpp index cc92d85..d1373c7 100644 --- a/tests/test_fft.cpp +++ b/tests/test_fft.cpp @@ -232,9 +232,10 @@ namespace { if constexpr (profile::k_fixed_point) { const int e = profile::fft::fixed_scaling_exponent(n); const double unit = sample_scale::k_lsb * std::ldexp(1.0, 2 * e + 1 - log2_of(n)); - std::printf("[ measured ] %s %s n=%zu: max error %.4f reconstructed LSB (unit %.3g; pin %.3f)\n", - sizeof(Sample) == 2 ? "Q15" : "Q31", test, n, max_error / unit, unit, - k_round_trip_units); + // %lu, not %zu: newlib's printf on the QEMU legs has no %zu. + std::printf("[ measured ] %s %s n=%lu: max error %.4f reconstructed LSB (unit %.3g; pin %.3f)\n", + sizeof(Sample) == 2 ? "Q15" : "Q31", test, static_cast(n), max_error / unit, + unit, k_round_trip_units); } } diff --git a/tests/test_fft_fixed.cpp b/tests/test_fft_fixed.cpp index f7fa1d5..88730cf 100644 --- a/tests/test_fft_fixed.cpp +++ b/tests/test_fft_fixed.cpp @@ -181,6 +181,13 @@ namespace { return std::ldexp(1.0, e); } + /// printf argument for a size: newlib's printf on the QEMU legs has no + /// %lu (it prints "zu" and desynchronises the varargs), so sizes go + /// through %lu as unsigned long everywhere in this file. + unsigned long ul(std::size_t v) { + return static_cast(v); + } + double db(double power) { return 10.0 * std::log10(std::max(power, 1e-300)); } @@ -669,9 +676,10 @@ namespace { } } } - std::printf("[ measured ] %s saturation sweep: max |out - G/2^e| = %.3f LSB (pin %.3f) at %s %s n=%zu e=%d " - "index %zu; largest rail shortfall %.3f LSB\n", - cfg::name(), worst.value, pin, worst.what, worst.name, worst.n, worst.e, worst.index, rail.value); + std::printf("[ measured ] %s saturation sweep: max |out - G/2^e| = %.3f LSB (pin %.3f) at %s %s n=%lu e=%d " + "index %lu; largest rail shortfall %.3f LSB\n", + cfg::name(), worst.value, pin, worst.what, worst.name, ul(worst.n), worst.e, ul(worst.index), + rail.value); EXPECT_GT(pin, 0.0) << "unmeasured pin"; EXPECT_LE(worst.value, pin) << cfg::name() << " " << worst.what << " " << worst.name << " n=" << worst.n << " e=" << worst.e << " index " << worst.index; @@ -732,9 +740,10 @@ namespace { } } } - std::printf("[ measured ] %s round trip: max error %.3f reconstructed LSB (pin %.3f) at %s n=%zu " - "e_fwd=%d e_inv=%d index %zu\n", - cfg::name(), worst.value, pin, worst.what, worst.n, worst.e / 100, worst.e % 100, worst.index); + std::printf("[ measured ] %s round trip: max error %.3f reconstructed LSB (pin %.3f) at %s n=%lu " + "e_fwd=%d e_inv=%d index %lu\n", + cfg::name(), worst.value, pin, worst.what, ul(worst.n), worst.e / 100, worst.e % 100, + ul(worst.index)); EXPECT_GT(pin, 0.0) << "unmeasured pin"; EXPECT_LE(worst.value, pin) << cfg::name() << " " << worst.what << " n=" << worst.n << " e_fwd=" << worst.e / 100 << " e_inv=" << worst.e % 100 << " index " @@ -815,10 +824,10 @@ namespace { } } } - std::printf("[ measured ] %s vs fixed after shift: max %.1f LSB (pin %.1f) at %s %s n=%zu e_fixed=%d e_bfp=%d " - "index %zu\n", - cfg::name(), worst.value, pin, worst.what, worst.name, worst.n, worst.e / 100, worst.e % 100, - worst.index); + std::printf("[ measured ] %s vs fixed after shift: max %.1f LSB (pin %.1f) at %s %s n=%lu e_fixed=%d e_bfp=%d " + "index %lu\n", + cfg::name(), worst.value, pin, worst.what, worst.name, ul(worst.n), worst.e / 100, worst.e % 100, + ul(worst.index)); EXPECT_GT(pin, 0.0) << "unmeasured pin"; EXPECT_LE(worst.value, pin) << cfg::name() << " " << worst.what << " " << worst.name << " n=" << worst.n << " e_fixed=" << worst.e / 100 << " e_bfp=" << worst.e % 100 << " index " @@ -849,8 +858,8 @@ namespace { // makes BFP shift like fixed at every stage. EXPECT_GT(attained, 0u) << cfg::name() << ": no full-scale pattern out of " << tried << " reached the fixed exponent"; - std::printf("[ measured ] %s: %zu of %zu full-scale patterns reach the fixed exponent\n", cfg::name(), attained, - tried); + std::printf("[ measured ] %s: %lu of %lu full-scale patterns reach the fixed exponent\n", cfg::name(), + ul(attained), ul(tried)); } // ======================================================================== @@ -897,10 +906,11 @@ namespace { } } } - std::printf("[ measured ] %s F(x)+F(-x): max %.2f LSB (pin %.2f) at %s %s n=%zu index %zu; bias %.4f LSB " - "(pin %.4f) at %s %s n=%zu\n", - cfg::name(), worst_max.value, max_pin, worst_max.what, worst_max.name, worst_max.n, worst_max.index, - worst_bias.value, bias_pin, worst_bias.what, worst_bias.name, worst_bias.n); + std::printf("[ measured ] %s F(x)+F(-x): max %.2f LSB (pin %.2f) at %s %s n=%lu index %lu; bias %.4f LSB " + "(pin %.4f) at %s %s n=%lu\n", + cfg::name(), worst_max.value, max_pin, worst_max.what, worst_max.name, ul(worst_max.n), + ul(worst_max.index), worst_bias.value, bias_pin, worst_bias.what, worst_bias.name, + ul(worst_bias.n)); EXPECT_GT(max_pin, 0.0) << "unmeasured pin"; EXPECT_GT(bias_pin, 0.0) << "unmeasured pin"; EXPECT_LE(worst_max.value, max_pin) << cfg::name() << " " << worst_max.what << " " << worst_max.name @@ -945,8 +955,8 @@ namespace { } } std::printf( - "[ measured ] Q15 vs Q31 (%s): max %.3f Q15 LSB (pin %.3f) at %s %s n=%zu e15=%d e31=%d index %zu\n", what, - worst.value, pin, worst.what, worst.name, worst.n, worst.e / 100, worst.e % 100, worst.index); + "[ measured ] Q15 vs Q31 (%s): max %.3f Q15 LSB (pin %.3f) at %s %s n=%lu e15=%d e31=%d index %lu\n", what, + worst.value, pin, worst.what, worst.name, ul(worst.n), worst.e / 100, worst.e % 100, ul(worst.index)); EXPECT_GT(pin, 0.0) << "unmeasured pin"; EXPECT_LE(worst.value, pin) << what << " " << worst.what << " " << worst.name << " n=" << worst.n << " e15=" << worst.e / 100 << " e31=" << worst.e % 100 << " index " << worst.index; @@ -1152,9 +1162,9 @@ namespace { noise_row row{material, n, level_db, r.exponent, db(signal), db(noise), db(model.noise), db(late.noise), noise / model.noise}; - std::printf("[ floor ] %s %-5s N=%5zu %4.0f dBFS e=%2d signal %7.2f floor %7.2f model %7.2f (late %7.2f) " + std::printf("[ floor ] %s %-5s N=%5lu %4.0f dBFS e=%2d signal %7.2f floor %7.2f model %7.2f (late %7.2f) " "dBFS/component ratio %.3f snr %6.2f dB\n", - Cfg::name(), material, n, level_db, r.exponent, row.signal_db, row.floor_db, row.model_db, + Cfg::name(), material, ul(n), level_db, r.exponent, row.signal_db, row.floor_db, row.model_db, row.model_late_db, row.ratio, row.signal_db - row.floor_db); return row; } @@ -1268,7 +1278,7 @@ namespace { ASSERT_LE(ei, slack) << "wki n=" << n << " k=" << k; } if (n <= 2048) { - std::printf("[ measured ] tables n=%zu: max |w_q - w| = %.6f LSB\n", n, worst); + std::printf("[ measured ] tables n=%lu: max |w_q - w| = %.6f LSB\n", ul(n), worst); } } } @@ -1300,8 +1310,10 @@ namespace { const auto post_sum = fnv1a64(post); // Printed before any assertion so a -V log carries every host's // values whether or not they match. - std::printf("[ checksum ] n=%zu twiddles fnv1a64=%016llx post-pass fnv1a64=%016llx\n", p.n, - static_cast(tw_sum), static_cast(post_sum)); + // Two 32-bit halves: newlib's printf has no %llx either. + std::printf("[ checksum ] n=%lu twiddles fnv1a64=%08lx%08lx post-pass fnv1a64=%08lx%08lx\n", ul(p.n), + static_cast(tw_sum >> 32), static_cast(tw_sum & 0xffffffffu), + static_cast(post_sum >> 32), static_cast(post_sum & 0xffffffffu)); // The kernel's own checksum helper is the same fold. EXPECT_EQ(tap::dsp::detail::table_checksum(twiddles.data(), twiddles.size()), tw_sum); EXPECT_EQ(tap::dsp::detail::table_checksum(post.data(), post.size()), post_sum);