diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index eaaa14b..93b23b3 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -140,3 +140,124 @@ tap_dsp_add_gtest_executable(tap_dsp_tests target_link_libraries(tap_dsp_tests PRIVATE tap::dsp tap_dsp_warnings) + +# ============================================================================== +# BEGIN Stage 2a additions (docs/audit-fft-and-code-smells.md, Part 9): the +# independent oracle and the real-time guard join the main executable; the +# Ooura bit-identity gate gets its own targets so it can own its compiler flags. +# Registration goes through tap_dsp_add_gtest_executable above, so the two +# targets are hosted (discovered) or bare-metal one-shot exactly like +# tap_dsp_tests. Keep this block self-contained and at the end of the file. +# ============================================================================== + +# Largest transform the FFT sweeps run. Hosted, the default covers every +# consumer geometry (2^20: five 8 MB buffers plus tables). Cross-compiled, the +# emulated Cortex-M legs (Part 10) get 4096, what fits the MPS2/MPS3 data +# regions; the toolchain files set the same value in the cache themselves and +# a plain cache set here never overrides theirs or a -D on the command line. +# The parity sweep, its 2^20 test (compiled out below the cap: a GTEST_SKIP is +# rc=1 on the target) and the oracle's closed-form sweep all read it. +if(CMAKE_CROSSCOMPILING) + set(_tap_dsp_parity_max_n_default 4096) +else() + set(_tap_dsp_parity_max_n_default 1048576) +endif() +set(TAP_DSP_PARITY_MAX_N ${_tap_dsp_parity_max_n_default} CACHE STRING + "Largest FFT size the Ooura parity and oracle sweeps run (power of two >= 4; 4096 on emulated targets)") + +target_sources(tap_dsp_tests PRIVATE + test_fft_oracle.cpp + test_fft_rt.cpp) +target_compile_definitions(tap_dsp_tests PRIVATE TAP_DSP_PARITY_MAX_N=${TAP_DSP_PARITY_MAX_N}) + +# ------------------------------------------------------------------------------ +# Two private reference builds of the vendored C. Neither is the tap_dsp_fft +# the library links: that target disappears at Stage 2c (it survives only under +# TAP_DSP_FFT_CMSIS, where it also carries CMSIS objects), and the gate's whole +# point is to compare against the C independently of how the library builds it. +# The parity executables link one of these and NOT tap::dsp, so their +# rdft/rdft_f resolve here and no backend define (vDSP on Apple, CMSIS on the +# M55) reaches them. +# +# Stage 2c note: fftsg.c moves to tests/reference/ooura/ (D6); this gate needs +# fftsg_float.c to move with it, or its float side becomes port-vs-port. +# +# ..._nocontract fp-contraction OFF. Bit identity between the C and its C++ +# transliteration holds only if both compile to the same +# sequence of IEEE operations, and CMake gives neither side an +# fp-contract setting: gcc contracts C in gnu17 mode, g++ +# contracts C++ in every mode, clang contracts everywhere +# (Part 4, "fp-contraction policy"). MSVC: nothing is passed; +# its default /fp:precise has not contracted since VS 2022 +# 17.0 made /fp:contract opt-in (earlier x64 /arch:AVX2 and +# ARM64 builds did contract under /fp:precise; windows-latest +# is well past 17.0). +# ..._default the same two files at default flags, for the informational +# target, so it measures "port at default flags vs C at +# default flags" and nothing else. +# ------------------------------------------------------------------------------ +set(_tap_dsp_nocontract_c "$<$:-ffp-contract=off>") +set(_tap_dsp_nocontract_cxx "$<$:-ffp-contract=off>") +set(_tap_dsp_reference_sources + ${PROJECT_SOURCE_DIR}/third_party/ooura/fftsg.c + ${PROJECT_SOURCE_DIR}/third_party/ooura/fftsg_float.c) + +add_library(tap_dsp_fft_reference_nocontract STATIC ${_tap_dsp_reference_sources}) +target_compile_options(tap_dsp_fft_reference_nocontract PRIVATE ${_tap_dsp_nocontract_c}) + +add_library(tap_dsp_fft_reference_default STATIC ${_tap_dsp_reference_sources}) + +# THE GATE: -ffp-contract=off on both sides, memcmp bit identity. On the QEMU +# legs it runs whole (no MAIN_FILTER): at N <= 4096 the sweep is milliseconds. +tap_dsp_add_gtest_executable(tap_dsp_fft_parity + SOURCES test_fft_parity_ooura.cpp + LABELS parity) +target_include_directories(tap_dsp_fft_parity PRIVATE ${PROJECT_SOURCE_DIR}/include) +target_compile_features(tap_dsp_fft_parity PRIVATE cxx_std_20) +target_compile_options(tap_dsp_fft_parity PRIVATE ${_tap_dsp_nocontract_cxx}) +target_compile_definitions(tap_dsp_fft_parity PRIVATE TAP_DSP_PARITY_MAX_N=${TAP_DSP_PARITY_MAX_N}) +target_link_libraries(tap_dsp_fft_parity PRIVATE + tap_dsp_fft_reference_nocontract + tap_dsp_warnings) + +# INFORMATIONAL: the same source at default flags on both sides, measuring the +# max-ulp deviation per N and never failing. It is not a pass/fail gate with an +# assumed bound because a different fusion choice per butterfly stage +# accumulates over log2 N stages; the number goes on the record per platform +# instead (Part 6, item N3). Two channels carry it, because ctest hides the +# stdout of a passing test under --output-on-failure: (1) CI runs the parity +# label with -V as its own step (ci.yml, Part 13); (2) the JUnit XML below. +# The report is ONE test covering both precisions so a single invocation +# produces the complete XML. +tap_dsp_add_gtest_executable(tap_dsp_fft_parity_default_flags + SOURCES test_fft_parity_ooura.cpp + LABELS parity) +target_include_directories(tap_dsp_fft_parity_default_flags PRIVATE ${PROJECT_SOURCE_DIR}/include) +target_compile_features(tap_dsp_fft_parity_default_flags PRIVATE cxx_std_20) +target_compile_definitions(tap_dsp_fft_parity_default_flags PRIVATE + TAP_DSP_PARITY_INFORMATIONAL + TAP_DSP_PARITY_MAX_N=${TAP_DSP_PARITY_MAX_N}) +target_link_libraries(tap_dsp_fft_parity_default_flags PRIVATE + tap_dsp_fft_reference_default + tap_dsp_warnings) + +# JUnit XML for the informational report (hosted only: bare metal has no file +# system). gtest reads GTEST_OUTPUT from the environment at start-up, and the +# helper exposes no per-test properties or extra arguments, so the ENVIRONMENT +# property is attached the way gtest_discover_tests itself attaches its +# registrations: a script on the directory's TEST_INCLUDE_FILES, which ctest +# runs after the discovery script that defines the test (if(TEST ...) is not +# available in that context, so the call is unconditional: a discovery failure +# then errors here as well as in its own _NOT_BUILT entry). Should the helper +# grow an EXTRA_ARGS/PROPERTIES pass-through, this collapses to one line there. +if(NOT TAP_DSP_BARE_METAL) + set(_tap_dsp_parity_xml_script ${CMAKE_CURRENT_BINARY_DIR}/tap_dsp_fft_parity_default_flags_env.cmake) + file(WRITE ${_tap_dsp_parity_xml_script} + "set_tests_properties(fft_parity_ooura_default_flags.ReportsMaxUlpVersusOoura PROPERTIES\n" + " ENVIRONMENT \"GTEST_OUTPUT=xml:${CMAKE_BINARY_DIR}/parity-ulp.xml\")\n") + set_property(DIRECTORY APPEND PROPERTY TEST_INCLUDE_FILES ${_tap_dsp_parity_xml_script}) +endif() + +# ============================================================================== +# END Stage 2a additions +# ============================================================================== diff --git a/tests/support/signals.h b/tests/support/signals.h new file mode 100644 index 0000000..0a6474d --- /dev/null +++ b/tests/support/signals.h @@ -0,0 +1,78 @@ +/// @file signals.h +/// @brief Shared test-signal generators for the DspTap test battery. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// One xorshift32, one random_signal and one tone synthesizer, so a new test +// file does not grow its own copy (three already exist in test_fft.cpp, +// test_fft_backend.cpp and test_nn.cpp; migrating them, and adding the dB +// 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. + +#pragma once + +#include +#include +#include +#include +#include + +namespace tap::dsp::test { + + /// 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. + class xorshift32 { + public: + explicit xorshift32(std::uint32_t seed) noexcept + : m_s(seed != 0u ? seed : 0x9E3779B9u) {} + + /// Next raw 32-bit state. + std::uint32_t next_u32() noexcept { + m_s ^= m_s << 13; + m_s ^= m_s >> 17; + m_s ^= m_s << 5; + return m_s; + } + + /// Uniform in [-1, 1), computed in double so float and double signals + /// drawn from the same seed are the same values up to the final cast. + double next_unit() noexcept { return static_cast(next_u32()) / 2147483648.0 - 1.0; } + + private: + std::uint32_t m_s; + }; + + /// n samples uniform in [-amplitude, amplitude), from a fixed seed. + 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()); + } + return x; + } + + /// amplitude * cos(2*pi*bin*j/n + phase) for j in [0, n). + /// + /// The angle is formed as 2*pi * fmod(bin*j, n) / n, not as (2*pi*bin/n) * j: + /// with a rounded per-sample increment the phase error grows like j*ulp, + /// and at n = 65536 a supposedly on-bin cosine is ~1e-11 off, a thousand + /// times the engine's own error — which is what a closed-form oracle would + /// then wrongly report. fmod is exact, bin*j is exact for integer bins below + /// 2^53, and the one division by a power-of-two n is exact, so an integer-bin + /// tone here is the sampled cosine to one rounding of cos() per sample. + template + std::vector tone(std::size_t n, double bin, double amplitude, double phase = 0.0) { + std::vector x(n); + 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)); + } + return x; + } + +} // namespace tap::dsp::test diff --git a/tests/test_fft_oracle.cpp b/tests/test_fft_oracle.cpp new file mode 100644 index 0000000..76845d8 --- /dev/null +++ b/tests/test_fft_oracle.cpp @@ -0,0 +1,755 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// THE INDEPENDENT ORACLE for tap::dsp::basic_real_fft (Stage 2a of +// docs/audit-fft-and-code-smells.md; Part 9 names this file). +// +// test_fft_parity_ooura.cpp proves the port IS Ooura, bit for bit. It cannot +// prove Ooura is a Fourier transform: if the reference C and the port shared a +// defect, parity would be green. This file answers that with two references +// that share nothing with fftsg.c: +// +// 1. Closed-form vectors whose transforms are known exactly — an impulse +// (flat spectrum), DC (N in slot 0), the Nyquist alternation (N in slot +// 1), an on-bin cosine and sine at bin k (+N/2 in the real, resp. IMAG, +// part of bin k — the plus sign in the imaginary part is the documented +// W = exp(+2*pi*i/N) convention, the conjugate of the engineering DFT), +// and a two-tone superposition — and their unnormalized inverses, whose +// gain is N/2 per the contract in fft.h. Checked to a tolerance DERIVED +// from the sample type's epsilon and log2 N (see "Tolerance" below). +// +// 2. A compensated-summation DFT for N <= 256: the O(N^2) definition with +// double-double accumulation (TwoSum / TwoProd, Dekker 1971 and Knuth +// TAOCP vol. 2 4.2.2; Shewchuk 1997 for the error-free transformations) +// and double-double twiddles from a Taylor series after exact octant +// reduction on the integer bin*sample index. Not `long double`: that is +// `double` on MSVC and on Apple arm64 (Part 6, item N19), so it would be +// no oracle at all on two of the three hosted CI legs. +// +// Relation to test_fft.cpp. That file's contract battery already pins the +// impulse, DC/Nyquist packing and the +i sign convention at one size each +// (N = 64 and 128) with a round tolerance; Part 9 lists the closed forms in +// both files on purpose. What this file adds is the sweep to N = 65536, the +// inverse closed forms, the derived (not round) tolerance, and the DFT +// 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. + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "support/signals.h" +#include "tap/dsp/fft.h" + +#ifndef TAP_DSP_PARITY_MAX_N +#define TAP_DSP_PARITY_MAX_N (1 << 20) +#endif + +namespace { + + // ------------------------------------------------------------------------ + // Tolerance. + // + // Higham, Accuracy and Stability of Numerical Algorithms (2nd ed., 2002), + // Theorem 24.2: a radix-2 FFT of length N = 2^t whose twiddles carry + // relative error at most mu computes y_hat with + // + // || y_hat - y ||_2 <= t * eta / (1 - t * eta) * || y ||_2, + // eta = mu + gamma_4 (1 + mu), gamma_4 ~ 4u, + // + // where u is the unit roundoff (epsilon / 2). mu is NOT one rounding here: + // Ooura's makewt takes cos/sin from libm for a quarter of the table and + // derives the rest arithmetically (w[2] = 0.5 / cos(2 delta), the halving + // recurrences 0.5 / wk1r, ...), and the float instantiation forms + // delta * j in float, so mu is a small multiple of u. With mu ~ 2-3u, + // eta ~ 6-7u and, for t * eta << 1, the 2-norm error is below about + // 7 * u * log2(N) * ||y||_2. The largest single-element error cannot exceed + // the 2-norm, and ||y||_2 comes free with the exact answer. Casting the + // closed-form input to Sample adds at most u * ||y||_2 (Parseval), and + // Ooura is split-radix rather than the theorem's radix-2, which changes the + // constant by a factor of order one. The bound is taken as + // + // k_higham_constant * epsilon * log2(N) * ||y||_2, k_higham_constant = 4 + // + // (epsilon = 2u, so 8u: the ~7u above plus the input cast), i.e. derived + // with slack, not fitted. Measured against it on x86-64 Linux (GCC 13, -O3, + // glibc, Ooura as the engine): the largest |error| / (epsilon * log2(N) * + // ||y||_2) over every test in this file is 0.25 for double and 0.17 for + // float, 16x and 23x inside the constant. + // + // What the bound cannot see. A 2-norm used per element is loose by up to + // sqrt(N): at float, N = 65536, DC input, the tolerance is ~0.5 on a + // spectrum whose one live bin is 65536, so a single empty bin that is off + // by 0.4 would pass. That is acceptable for the closed forms (their job is + // the sweep over sizes and the sign convention) and is exactly why the DFT + // comparison is confined to N <= 256, where sqrt(N) <= 16 and the reference + // is exact to ~1e-30: there the check is tight enough to catch a single + // wrong bin. Precision pins as measured numbers live in test_fft.cpp and + // test_fft_backend.cpp; this file is an oracle for correctness. + // ------------------------------------------------------------------------ + constexpr double k_higham_constant = 4.0; + + double higham_tolerance(double epsilon, std::size_t n, double spectrum_norm2) { + const double log2n = std::log2(static_cast(n)); + return k_higham_constant * epsilon * log2n * spectrum_norm2; + } + + // ------------------------------------------------------------------------ + // Per-profile traits — THE EXTENSION POINT FOR THE FIXED-POINT STAGE. + // + // 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. + // ------------------------------------------------------------------------ + template + struct profile; + + // Size range a profile's engine accepts. fft.h promises any power of two + // >= 4 for Ooura, and the double profile is always Ooura. The float + // profile is whatever backend the build selected: under TAP_DSP_FFT_CMSIS + // (the M55 leg) CMSIS-DSP's arm_rfft_fast_init_f32 accepts 32..4096 only + // (arm_rfft_fast_init_f32.c, the switch at the end) and fft.h's wrapper + // does not check its return status, so a size outside that range is + // undefined behaviour — a HardFault at N = 4 on the QEMU M55 leg is how + // this was found. Until fft.h rejects or falls back on those sizes (a + // finding for the fft.h owner, not this file), the float sweeps on that + // backend run over the range the backend supports. vDSP (macOS) takes the + // full range. + constexpr std::size_t k_ooura_min_n = 4; + constexpr std::size_t k_ooura_max_n = std::size_t{1} << 20; +#if defined(TAP_DSP_FFT_CMSIS) + constexpr std::size_t k_float_backend_min_n = 32; + constexpr std::size_t k_float_backend_max_n = 4096; +#else + constexpr std::size_t k_float_backend_min_n = k_ooura_min_n; + constexpr std::size_t k_float_backend_max_n = k_ooura_max_n; +#endif + + template <> + struct profile { + static constexpr double k_epsilon = std::numeric_limits::epsilon(); + static constexpr double k_full_scale = 1.0; ///< closed-form drive amplitude + static constexpr std::size_t k_min_n = k_ooura_min_n; + static constexpr std::size_t k_max_n = k_ooura_max_n; + static double to_double(double v) { return v; } + static double from_double(double v) { return v; } + static double forward_scale(std::size_t) { return 1.0; } ///< engine forward = scale * DFT + static double inverse_scale(std::size_t) { return 1.0; } ///< engine inverse = scale * unnormalized + static double tolerance(std::size_t n, double norm2) { return higham_tolerance(k_epsilon, n, norm2); } + }; + + template <> + struct profile { + static constexpr double k_epsilon = std::numeric_limits::epsilon(); + static constexpr double k_full_scale = 1.0; + static constexpr std::size_t k_min_n = k_float_backend_min_n; + static constexpr std::size_t k_max_n = k_float_backend_max_n; + static double to_double(float v) { return static_cast(v); } + static float from_double(double v) { return static_cast(v); } + static double forward_scale(std::size_t) { return 1.0; } + static double inverse_scale(std::size_t) { return 1.0; } + static double tolerance(std::size_t n, double norm2) { return higham_tolerance(k_epsilon, n, norm2); } + }; + + 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) { + double e = 0.0; + for (const double v : x) { + e += v * v; + } + return std::sqrt(static_cast(x.size()) * e); + } + + std::vector scale_by(std::vector x, double factor) { + for (double& v : x) { + v *= factor; + } + return x; + } + + template + std::vector to_doubles(const std::vector& x) { + std::vector d(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) { + d[i] = profile::to_double(x[i]); + } + return d; + } + + template + std::vector from_doubles(const std::vector& x) { + std::vector s(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) { + s[i] = profile::from_double(x[i]); + } + return s; + } + + // ------------------------------------------------------------------------ + // Double-double arithmetic: an unevaluated sum hi + lo with |lo| <= ulp(hi)/2, + // ~106 bits of significand. Only what the DFT below needs. + // ------------------------------------------------------------------------ + struct dd { + double hi; + double lo; + }; + + /// Knuth TwoSum: s + e == a + b exactly, no ordering assumption. + inline dd two_sum(double a, double b) { + const double s = a + b; + const double bb = s - a; + const double e = (a - (s - bb)) + (b - bb); + return {s, e}; + } + + /// Dekker FastTwoSum, valid when |a| >= |b|. + inline dd quick_two_sum(double a, double b) { + const double s = a + b; + const double e = b - (s - a); + return {s, e}; + } + + /// TwoProd via a correctly rounded fma: p + e == a * b exactly. std::fma + /// is required to round once, whatever the hardware. + inline dd two_prod(double a, double b) { + const double p = a * b; + const double e = std::fma(a, b, -p); + return {p, e}; + } + + inline dd dd_add(dd a, dd b) { + dd s = two_sum(a.hi, b.hi); + dd t = two_sum(a.lo, b.lo); + s.lo += t.hi; + s = quick_two_sum(s.hi, s.lo); + s.lo += t.lo; + return quick_two_sum(s.hi, s.lo); + } + + inline dd dd_neg(dd a) { + return {-a.hi, -a.lo}; + } + + inline dd dd_mul(dd a, dd b) { + dd p = two_prod(a.hi, b.hi); + p.lo += a.hi * b.lo + a.lo * b.hi; + return quick_two_sum(p.hi, p.lo); + } + + inline dd dd_mul_d(dd a, double b) { + dd p = two_prod(a.hi, b); + p.lo += a.lo * b; + return quick_two_sum(p.hi, p.lo); + } + + inline dd dd_div_d(dd a, double b) { + const double q1 = a.hi / b; + const dd r = dd_add(a, dd_neg(dd_mul_d({b, 0.0}, q1))); + const double q2 = r.hi / b; + return quick_two_sum(q1, q2); + } + + inline double dd_round(dd a) { + return a.hi + a.lo; + } + + // pi as a double-double: the standard 3.141592653589793 + 1.2246467991473532e-16. + constexpr dd k_pi_dd{3.141592653589793, 1.2246467991473532e-16}; + + /// sin and cos of |theta| <= pi/4 by Taylor series in double-double. Terms + /// fall below 1e-34 relative by k = 15 for |theta| <= pi/4; 20 terms are + /// summed for margin, largest first, which double-double handles exactly + /// enough for a 1e-30 result. The oracle needs ~1e-20, so this is ample. + struct dd_sincos { + dd sin; + dd cos; + }; + + dd_sincos sincos_small(dd theta) { + const dd theta2 = dd_mul(theta, theta); + dd s_term = theta; + dd c_term = {1.0, 0.0}; + dd s = s_term; + dd c = c_term; + for (int k = 1; k <= 20; ++k) { + // sin: term_k = -term_{k-1} * theta^2 / ((2k)(2k+1)) + // cos: term_k = -term_{k-1} * theta^2 / ((2k-1)(2k)) + s_term = dd_neg(dd_div_d(dd_mul(s_term, theta2), static_cast((2 * k) * (2 * k + 1)))); + c_term = dd_neg(dd_div_d(dd_mul(c_term, theta2), static_cast((2 * k - 1) * (2 * k)))); + s = dd_add(s, s_term); + c = dd_add(c, c_term); + } + return {s, c}; + } + + /// (cos, sin)(2*pi*m/n) for integer m in [0, n) and n a power of two >= 4, + /// by exact octant reduction on the integers (m and n are exact, m/n is + /// exact because n is a power of two) followed by the small-angle series. + dd_sincos twiddle(std::size_t m, std::size_t n) { + const std::size_t half = n / 2; + const std::size_t quarter = n / 4; + bool negate = false; + bool rotate = false; // (cos, sin) -> (-sin, cos) + m %= n; + if (m >= half) { + m -= half; + negate = true; + } + if (m >= quarter) { + m -= quarter; + rotate = true; + } + // m in [0, quarter): theta in [0, pi/2). Fold to [0, pi/4] via + // cos(theta) = sin(pi/2 - theta). + bool swap = false; + if (2 * m > quarter) { + m = quarter - m; + swap = true; + } + const double fraction = static_cast(m) / static_cast(n); // exact + const dd theta = dd_mul_d({2.0 * k_pi_dd.hi, 2.0 * k_pi_dd.lo}, fraction); + const dd_sincos sc = sincos_small(theta); + dd c = swap ? sc.sin : sc.cos; + dd s = swap ? sc.cos : sc.sin; + if (rotate) { + const dd t = c; + c = dd_neg(s); + s = t; + } + if (negate) { + c = dd_neg(c); + s = dd_neg(s); + } + return {s, c}; + } + + // ------------------------------------------------------------------------ + // The compensated DFT, in Ooura's packing and sign convention (fft.h): + // forward: a[0] = Re X[0], a[1] = Re X[N/2], a[2k] = Re X[k], a[2k+1] = Im X[k], + // X[k] = sum_j x[j] * exp(+2*pi*i*j*k/N) + // inverse (unnormalized, gain N/2 on a round trip): + // x[k] = (R[0] + R[N/2] * (-1)^k) / 2 + // + sum_{j=1}^{N/2-1} (R[j] cos(2*pi*j*k/N) + I[j] sin(2*pi*j*k/N)) + // The inverse formula is fftsg.c's own statement of what rdft(n, -1, ...) + // computes; both are restated from the definition, not from the code path. + // ------------------------------------------------------------------------ + class compensated_dft { + public: + explicit compensated_dft(std::size_t n) + : m_n(n) + , m_twiddles(n) { + for (std::size_t m = 0; m < n; ++m) { + m_twiddles[m] = twiddle(m, n); + } + } + + std::vector forward(const std::vector& x) const { + std::vector a(m_n, 0.0); + for (std::size_t k = 0; k <= m_n / 2; ++k) { + dd re{0.0, 0.0}; + dd im{0.0, 0.0}; + for (std::size_t j = 0; j < m_n; ++j) { + const dd_sincos& w = m_twiddles[(j * k) % m_n]; + re = dd_add(re, dd_mul_d(w.cos, x[j])); + im = dd_add(im, dd_mul_d(w.sin, x[j])); + } + if (k == 0) { + a[0] = dd_round(re); + } + else if (k == m_n / 2) { + a[1] = dd_round(re); + } + else { + a[2 * k] = dd_round(re); + a[2 * k + 1] = dd_round(im); + } + } + return a; + } + + std::vector inverse_unnormalized(const std::vector& a) const { + std::vector x(m_n, 0.0); + for (std::size_t k = 0; k < m_n; ++k) { + // (R[0] +- R[N/2]) / 2: the sum error-free via TwoSum, the + // halving exact. + dd acc = dd_mul_d(two_sum(a[0], (k % 2 == 0) ? a[1] : -a[1]), 0.5); + for (std::size_t j = 1; j < m_n / 2; ++j) { + const dd_sincos& w = m_twiddles[(j * k) % m_n]; + acc = dd_add(acc, dd_mul_d(w.cos, a[2 * j])); + acc = dd_add(acc, dd_mul_d(w.sin, a[2 * j + 1])); + } + x[k] = dd_round(acc); + } + return x; + } + + const dd_sincos& twiddle_at(std::size_t m) const { return m_twiddles[m % m_n]; } + + private: + std::size_t m_n; + std::vector m_twiddles; + }; + + // ------------------------------------------------------------------------ + // Shared checkers. + // ------------------------------------------------------------------------ + template + void expect_close(const std::vector& got, const std::vector& expected, double tol, const char* what, + std::size_t n) { + ASSERT_EQ(got.size(), expected.size()); + for (std::size_t i = 0; i < got.size(); ++i) { + const double g = profile::to_double(got[i]); + ASSERT_NEAR(g, expected[i], tol) << what << " N=" << n << " index " << i; + } + } + + /// Powers of two from the profile's minimum up to min(limit, + /// TAP_DSP_PARITY_MAX_N, the profile's maximum). The cap is the same knob + /// the parity gate uses (tests/CMakeLists.txt); the QEMU legs set 4096 + /// because a 65536-point double sweep needs several 512 KB buffers that + /// the MPS2 data region does not have. + template + std::vector sizes_up_to(std::size_t limit) { + const std::size_t top = + std::min({limit, static_cast(TAP_DSP_PARITY_MAX_N), profile::k_max_n}); + std::vector s; + for (std::size_t n = profile::k_min_n; n <= top; n *= 2) { + s.push_back(n); + } + return s; + } + + /// The closed-form sweep: up to 65536. + template + std::vector closed_form_sizes() { + return sizes_up_to(65536); + } + + /// The compensated-DFT sizes: up to 256 (see "What the bound cannot see"). + template + std::vector dft_sizes() { + return sizes_up_to(256); + } + + /// Bin used for the on-bin materials: n/8, or 1 below n = 8. Always in + /// [1, n/2), i.e. a genuinely complex bin, never DC or Nyquist. + std::size_t tone_bin(std::size_t n) { + return std::max(1, n / 8); + } + + // Runs the engine forward on x (given in the double domain, cast to the + // profile) and checks it against the exact spectrum `expected`, with the + // profile's forward scale applied. + template + void check_forward(const std::vector& x, const std::vector& expected, const char* what) { + const std::size_t n = x.size(); + tap::dsp::basic_real_fft fft(n); + std::vector buf = from_doubles(x); + fft.forward_inplace(buf.data()); + const std::vector scaled = scale_by(expected, profile::forward_scale(n)); + expect_close(buf, scaled, profile::tolerance(n, spectrum_norm2(x)), what, n); + } + + // Runs the engine's UNNORMALIZED inverse on the packed spectrum `a` and + // checks it against `expected` (which already carries the N/2 gain), with + // the profile's inverse scale applied. + template + void check_inverse(const std::vector& a, const std::vector& expected, const char* what) { + const std::size_t n = a.size(); + tap::dsp::basic_real_fft fft(n); + std::vector buf = from_doubles(a); + fft.inverse_inplace(buf.data()); + const std::vector scaled = scale_by(expected, profile::inverse_scale(n)); + // ||expected||_2 plays the role of ||y||_2 for the inverse direction; + // the packed spectrum's 2-norm is within sqrt(2) of the true one, and + // the N/2 gain means ||expected||_2 >= ||a||_2 * sqrt(N)/2 for these + // vectors, so the bound stays of the same form. + double e = 0.0; + for (const double v : expected) { + e += v * v; + } + expect_close(buf, scaled, profile::tolerance(n, std::sqrt(e)), what, n); + } + + template + class fft_oracle_test : public ::testing::Test {}; + TYPED_TEST_SUITE(fft_oracle_test, oracle_types); + + // ======================================================================== + // Closed forms, forward. Every one is driven at the profile's full-scale + // amplitude (1.0 for float and double); the exact answers scale with it. + // ======================================================================== + + TYPED_TEST(fft_oracle_test, ImpulseIsFlatAtEverySize) { + const double a = profile::k_full_scale; + for (const std::size_t n : closed_form_sizes()) { + std::vector x(n, 0.0); + x[0] = a; + std::vector expected(n, 0.0); + expected[0] = a; // DC + expected[1] = a; // Nyquist + for (std::size_t k = 1; k < n / 2; ++k) { + expected[2 * k] = a; + } + check_forward(x, expected, "impulse"); + } + } + + TYPED_TEST(fft_oracle_test, DcLandsAsNInSlotZero) { + const double a = profile::k_full_scale; + for (const std::size_t n : closed_form_sizes()) { + std::vector x(n, a); + std::vector expected(n, 0.0); + expected[0] = a * static_cast(n); + check_forward(x, expected, "dc"); + } + } + + TYPED_TEST(fft_oracle_test, NyquistAlternationLandsAsNInSlotOne) { + const double a = profile::k_full_scale; + for (const std::size_t n : closed_form_sizes()) { + std::vector x(n); + for (std::size_t j = 0; j < n; ++j) { + x[j] = (j % 2 == 0) ? a : -a; + } + std::vector expected(n, 0.0); + expected[1] = a * static_cast(n); + check_forward(x, expected, "nyquist"); + } + } + + TYPED_TEST(fft_oracle_test, OnBinCosineIsPlusHalfNInTheRealSlot) { + const double a = profile::k_full_scale; + for (const std::size_t n : closed_form_sizes()) { + const std::size_t k = tone_bin(n); + const auto x = tap::dsp::test::tone(n, static_cast(k), a, 0.0); + std::vector expected(n, 0.0); + expected[2 * k] = a * static_cast(n) / 2.0; + check_forward(x, expected, "cosine"); + } + } + + // The sign convention: a sine at bin k lands at +N/2 in the imaginary slot + // (it would be -N/2 in the engineering convention exp(-2*pi*i/N)). + TYPED_TEST(fft_oracle_test, OnBinSineIsPlusHalfNInTheImaginarySlot) { + const double a = profile::k_full_scale; + for (const std::size_t n : closed_form_sizes()) { + const std::size_t k = tone_bin(n); + // sin(w j) = cos(w j - pi/2) + const auto x = tap::dsp::test::tone(n, static_cast(k), a, -std::numbers::pi / 2.0); + std::vector expected(n, 0.0); + expected[2 * k + 1] = a * static_cast(n) / 2.0; + check_forward(x, expected, "sine"); + } + } + + TYPED_TEST(fft_oracle_test, TwoToneSuperposesLinearly) { + const double a = profile::k_full_scale; + for (const std::size_t n : closed_form_sizes()) { + if (n < 8) { + continue; // needs two distinct complex bins + } + const std::size_t k1 = tone_bin(n); + const std::size_t k2 = k1 + 1; + const double a1 = 0.75 * a; + const double a2 = -0.25 * a; + const auto c = tap::dsp::test::tone(n, static_cast(k1), a1, 0.0); + const auto s = tap::dsp::test::tone(n, static_cast(k2), a2, -std::numbers::pi / 2.0); + std::vector x(n); + for (std::size_t j = 0; j < n; ++j) { + x[j] = c[j] + s[j]; + } + std::vector expected(n, 0.0); + expected[2 * k1] = a1 * static_cast(n) / 2.0; + expected[2 * k2 + 1] = a2 * static_cast(n) / 2.0; + check_forward(x, expected, "two-tone"); + } + } + + // ======================================================================== + // Closed forms, inverse (unnormalized: gain N/2). + // ======================================================================== + + TYPED_TEST(fft_oracle_test, InverseOfFlatSpectrumIsHalfNImpulse) { + const double amp = profile::k_full_scale; + for (const std::size_t n : closed_form_sizes()) { + std::vector a(n, 0.0); + a[0] = amp; + a[1] = amp; + for (std::size_t k = 1; k < n / 2; ++k) { + a[2 * k] = amp; + } + std::vector expected(n, 0.0); + expected[0] = amp * static_cast(n) / 2.0; + check_inverse(a, expected, "inverse impulse"); + } + } + + TYPED_TEST(fft_oracle_test, InverseOfDcSpectrumIsHalfNConstant) { + const double amp = profile::k_full_scale; + for (const std::size_t n : closed_form_sizes()) { + std::vector a(n, 0.0); + a[0] = amp; + std::vector expected(n, 0.5 * amp); + check_inverse(a, expected, "inverse dc"); + } + } + + TYPED_TEST(fft_oracle_test, InverseOfNyquistSpectrumIsHalfNAlternation) { + const double amp = profile::k_full_scale; + for (const std::size_t n : closed_form_sizes()) { + std::vector a(n, 0.0); + a[1] = amp; + std::vector expected(n); + for (std::size_t k = 0; k < n; ++k) { + expected[k] = (k % 2 == 0) ? 0.5 * amp : -0.5 * amp; + } + check_inverse(a, expected, "inverse nyquist"); + } + } + + TYPED_TEST(fft_oracle_test, InverseOfOnBinSpectrumIsHalfNTone) { + const double amp = profile::k_full_scale; + for (const std::size_t n : closed_form_sizes()) { + const std::size_t k = tone_bin(n); + std::vector a(n, 0.0); + a[2 * k] = amp; // cosine part + a[2 * k + 1] = amp; // sine part, +i convention + // x[j] = amp * (cos(w j) + sin(w j)): the unnormalized inverse has + // gain N/2 on the round trip, so a bin of amp comes back with + // amplitude amp. + const auto c = tap::dsp::test::tone(n, static_cast(k), amp, 0.0); + const auto s = tap::dsp::test::tone(n, static_cast(k), amp, -std::numbers::pi / 2.0); + std::vector expected(n); + for (std::size_t j = 0; j < n; ++j) { + expected[j] = c[j] + s[j]; + } + check_inverse(a, expected, "inverse tone"); + } + } + + // ======================================================================== + // The compensated DFT, N <= 256. + // ======================================================================== + + // The oracle checks itself first: its twiddles agree with libm to double + // rounding, and its own forward followed by its own inverse is the + // identity to double rounding, so a disagreement below is the engine's. + // The two constants are derived bounds, not fitted numbers. Twiddles: the + // comparison forms theta = 2*pi*m/n in double, and rounding theta (at most + // half an ulp of 2*pi, 4.4e-16 = 2 eps) moves cos/sin by up to that much; + // libm adds at most 1 ulp (1 eps for values in [0.5, 1]); dd_round adds + // half an ulp. Bound 3.5 eps, constant 4 eps; measured maximum 6.4e-16 = + // 2.9 eps (x86-64 Linux, glibc 2.39), i.e. the dd twiddles are as good as + // libm and the argument rounding dominates. Round trip: forward then + // inverse re-rounds through two dd_round steps plus the final 2/N + // multiply (exact, power of two), so ~1.5 eps relative on unit-scale data; + // constant 8 eps; measured maximum 1.1e-16 = 0.5 eps. + TEST(fft_oracle_self_check, TwiddlesMatchLibmToDoubleRounding) { + for (const std::size_t n : dft_sizes()) { + const compensated_dft oracle(n); + for (std::size_t m = 0; m < n; ++m) { + const double theta = 2.0 * std::numbers::pi * static_cast(m) / static_cast(n); + const auto& w = oracle.twiddle_at(m); + EXPECT_NEAR(dd_round(w.cos), std::cos(theta), 4.0 * std::numeric_limits::epsilon()) + << "N=" << n << " m=" << m; + EXPECT_NEAR(dd_round(w.sin), std::sin(theta), 4.0 * std::numeric_limits::epsilon()) + << "N=" << n << " m=" << m; + // The low word is a genuine correction, not noise: it is at + // most half an ulp of the high word. + EXPECT_LE(std::fabs(w.cos.lo), std::fabs(w.cos.hi) * std::numeric_limits::epsilon() + 1e-300); + } + } + } + + TEST(fft_oracle_self_check, RoundTripIsIdentityToDoubleRounding) { + for (const std::size_t n : dft_sizes()) { + const compensated_dft oracle(n); + const std::vector x = tap::dsp::test::random_signal(n, 0xC0FFEEu); + const std::vector a = oracle.forward(x); + const std::vector back = oracle.inverse_unnormalized(a); + const double gain = 2.0 / static_cast(n); + for (std::size_t i = 0; i < n; ++i) { + EXPECT_NEAR(back[i] * gain, x[i], 8.0 * std::numeric_limits::epsilon()) + << "N=" << n << " " << i; + } + } + } + + template + void expect_forward_matches_dft(std::size_t n, const std::vector& x, const char* what) { + const compensated_dft oracle(n); + const std::vector expected = oracle.forward(x); + check_forward(x, expected, what); + } + + template + void expect_inverse_matches_dft(std::size_t n, const std::vector& a, const char* what) { + const compensated_dft oracle(n); + const std::vector expected = oracle.inverse_unnormalized(a); + check_inverse(a, expected, what); + } + + TYPED_TEST(fft_oracle_test, ForwardMatchesCompensatedDftOnBroadband) { + for (const std::size_t n : dft_sizes()) { + // Drawn in the profile first so the oracle sees exactly the bits + // the engine sees; from_doubles is then the identity. + const auto x = to_doubles(tap::dsp::test::random_signal(n, 0x2545F491u)); + expect_forward_matches_dft(n, x, "dft forward broadband"); + } + } + + TYPED_TEST(fft_oracle_test, ForwardMatchesCompensatedDftOnOffBinTone) { + for (const std::size_t n : dft_sizes()) { + // Off-bin (k + 0.37) so every bin carries leakage and none is + // numerically empty: the opposite regime from the closed forms. + const auto x = + to_doubles(tap::dsp::test::tone(n, static_cast(tone_bin(n)) + 0.37, 0.8, 1.1)); + expect_forward_matches_dft(n, x, "dft forward off-bin tone"); + } + } + + TYPED_TEST(fft_oracle_test, InverseMatchesCompensatedDftOnBroadbandSpectrum) { + for (const std::size_t n : dft_sizes()) { + // An arbitrary packed spectrum, not one produced by a forward, so + // the inverse is judged on its own definition. + const auto a = to_doubles(tap::dsp::test::random_signal(n, 0x1D872B41u)); + expect_inverse_matches_dft(n, a, "dft inverse broadband"); + } + } + + TYPED_TEST(fft_oracle_test, InverseMatchesCompensatedDftOnToneSpectrum) { + for (const std::size_t n : dft_sizes()) { + const compensated_dft oracle(n); + const auto x = + to_doubles(tap::dsp::test::tone(n, static_cast(tone_bin(n)) + 0.37, 0.8, 1.1)); + // The spectrum is rounded to the profile before both sides see it. + const auto a = to_doubles(from_doubles(oracle.forward(x))); + expect_inverse_matches_dft(n, a, "dft inverse tone spectrum"); + } + } + +} // namespace diff --git a/tests/test_fft_parity_ooura.cpp b/tests/test_fft_parity_ooura.cpp new file mode 100644 index 0000000..410f11e --- /dev/null +++ b/tests/test_fft_parity_ooura.cpp @@ -0,0 +1,426 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// THE BIT-IDENTITY GATE FOR THE C++20 PORT OF OOURA'S rdft (Stage 2a of +// docs/audit-fft-and-code-smells.md, Part 3; the policy behind it is Part 4's +// "fp-contraction policy" section and Part 6 items N3/N4). +// +// Two sides, each re-pointable in ONE place below: +// +// ooura_ref the raw rdft / rdft_f of Takuya Ooura's fftsg.c, +// from the REFERENCE copy of the C that this target's +// CMake block compiles (tests/CMakeLists.txt). +// engine_under_test whatever the library routes the real FFT to. Today +// that is basic_real_fft, i.e. the same C, so this +// suite is Ooura-vs-Ooura and trivially green. That +// is the point: the port agent re-points this alias +// at detail::split_radix_rdft, watches the suite go +// red, and works it green statement by statement. +// +// The comparison is memcmp over the raw output bytes: not EXPECT_EQ (which +// calls +0.0 and -0.0 equal and any NaN unequal to itself), not a tolerance. +// Forward and inverse are both taken in place and unnormalized, exactly as the +// C exposes them, for every power of two from 4 to 65536 plus one run at 2^20 +// (the largest size any consumer uses; item N20), on five materials chosen to +// reach different arithmetic: broadband xorshift noise (every butterfly busy), +// an on-bin tone (one live bin, the rest numerically empty), an impulse (flat +// spectrum, no cancellation), DC (maximal cancellation in every non-zero bin) +// and a full-scale alternating +1/-1 (all energy at Nyquist). +// +// Why this target owns its compiler flags. Bit identity between a C function +// and a C++ transliteration of it holds only if both are compiled to the same +// sequence of IEEE operations. Measured for Part 4: gcc -std=c17 does not +// contract a*b+c into an FMA; gcc -std=gnu17 does; g++ contracts in both c++20 +// and gnu++20; clang contracts in every mode, statement-scoped. CMake sets none +// of this, so on any ISA with FMA (Apple arm64, Cortex-M55, x86 with -march) +// the two sides fuse differently by default and the identity is false for +// reasons that have nothing to do with the port. The gate therefore compiles +// BOTH the reference C and this TU with -ffp-contract=off (MSVC: /fp:precise +// does not contract, so nothing is passed there). Same binary, same libm: +// cross-platform identity is not claimed and does not hold today either, +// because libm's cos/sin differ in the last bit between glibc, newlib, UCRT +// and Apple (Part 4, "The float twiddles"). +// +// To see that the flag is load-bearing rather than take it on faith, count +// the fused instructions in the reference object on an FMA-capable target: +// +// cc -O3 -march=haswell [-ffp-contract=off] -c third_party/ooura/fftsg.c -o probe.o +// objdump -d probe.o | grep -ciE 'vfmadd|vfmsub|vfnmadd|vfnmsub' +// +// Measured: gcc 13.3 200 / clang 18.1 171 fused instructions by default, 0 +// with the flag, for both compilers. On plain x86-64 (no -march) the count is +// 0 either way because the ISA has no FMA, which is why the hosted CI legs +// cannot see a missing flag; the objdump check is the one that can. +// +// The same source builds a SECOND, informational target at default flags +// (TAP_DSP_PARITY_INFORMATIONAL). It measures the max-ulp deviation per N and +// never fails. It is not a pass/fail gate with an assumed bound because a +// different fusion choice per butterfly stage accumulates over log2 N stages +// and "1 ulp" would be a guess, not a measurement (item N3); its job is to +// put the number on the record for each platform so the policy fft.h states +// about exporting -ffp-contract=off can be decided on data. Where the number +// lands: ctest hides the stdout of a passing test under --output-on-failure, +// so the table is visible only when the label is run with -V (CI does that +// as its own step, Part 13) and, as RecordProperty values, in the JUnit XML +// that the CMake block requests with --gtest_output=xml. +// +// TAP_DSP_PARITY_MAX_N caps the sizes run, so the emulated QEMU legs can take +// the suite at N <= 4096 (Part 9, Part 10) where 2^20 does not fit in RAM; +// the CMake block and the toolchain files default it to 4096 when +// cross-compiling, and the 2^20 test is compiled out (not skipped) below it. +// +// Stage 2c note. fftsg.c moves to tests/reference/ooura/ (Decision D6), and +// the 2c text deletes fftsg_float.c. This gate needs BOTH files: the float +// side compares against rdft_f, and without fftsg_float.c it degenerates to +// port-vs-port. 2c must move fftsg_float.c alongside fftsg.c, or delete the +// float half of this file and say so in the 2c PR. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "support/signals.h" +#include "tap/dsp/fft.h" + +#ifndef TAP_DSP_PARITY_MAX_N +#define TAP_DSP_PARITY_MAX_N (1 << 20) +#endif +static_assert(TAP_DSP_PARITY_MAX_N >= 4 && (TAP_DSP_PARITY_MAX_N & (TAP_DSP_PARITY_MAX_N - 1)) == 0, + "TAP_DSP_PARITY_MAX_N must be a power of two >= 4, or every sweep below is vacuously green"); + +namespace { + + // ------------------------------------------------------------------------ + // The reference side. Re-point here, nowhere else. + // + // Raw Ooura, called exactly as fftsg.c documents: ip[0] = 0 requests table + // initialization on the first call, and the workspace geometry is the one + // readme.txt prescribes (ip: 2 + sqrt(n/2), w: n/2). rdft / rdft_f come + // from the reference C library this target links (see the CMake block), + // not from the library's tap_dsp_fft, so the reference stays the C even + // after the library stops linking it (Stage 2c). + // ------------------------------------------------------------------------ + inline void ooura_rdft_ref(int n, int isgn, double* a, int* ip, double* w) { + rdft(n, isgn, a, ip, w); + } + inline void ooura_rdft_ref(int n, int isgn, float* a, int* ip, float* w) { + rdft_f(n, isgn, a, ip, w); + } + + template + class ooura_ref { + public: + explicit ooura_ref(std::size_t n) + : m_n(static_cast(n)) + , m_ip(2 + static_cast(std::sqrt(static_cast(n) / 2.0)) + 1, 0) + , m_w(n / 2, Sample(0)) { + m_ip[0] = 0; + } + void forward_inplace(Sample* a) { ooura_rdft_ref(m_n, 1, a, m_ip.data(), m_w.data()); } + void inverse_inplace(Sample* a) { ooura_rdft_ref(m_n, -1, a, m_ip.data(), m_w.data()); } + + private: + int m_n; + std::vector m_ip; + std::vector m_w; + }; + + // ------------------------------------------------------------------------ + // The side under test. Re-point here, nowhere else. + // + // Stage 2a (the port lands beside the C, nothing routed): re-point at + // tap::dsp::detail::split_radix_rdft + // Stage 2b (routing flipped): back to basic_real_fft, which then IS the + // port, and the suite guards the flip. + // Both expose the constructor-from-size / forward_inplace / inverse_inplace + // surface Part 4 fixes, so nothing else in this file changes. + // ------------------------------------------------------------------------ + template + using engine_under_test = tap::dsp::basic_real_fft; + + // ------------------------------------------------------------------------ + // Materials. + // ------------------------------------------------------------------------ + enum class material { broadband, on_bin_tone, impulse, dc, alternating }; + + constexpr material k_materials[] = {material::broadband, material::on_bin_tone, material::impulse, material::dc, + material::alternating}; + + const char* name_of(material m) { + switch (m) { + case material::broadband: + return "broadband"; + case material::on_bin_tone: + return "on-bin tone"; + case material::impulse: + return "impulse"; + case material::dc: + return "dc"; + case material::alternating: + return "alternating +-1"; + } + return "?"; + } + + template + std::vector make_material(material m, std::size_t n) { + switch (m) { + case material::broadband: + return tap::dsp::test::random_signal(n, 0x9E3779B9u); + case material::on_bin_tone: + // Bin n/8 (bin 1 below n = 8): far from DC and Nyquist so the + // butterflies that combine it see a real twiddle, not +-1 or +-i. + return tap::dsp::test::tone(n, static_cast(std::max(1, n / 8)), 0.5, 0.3); + case material::impulse: { + std::vector x(n, Sample(0)); + x[0] = Sample(1); + return x; + } + case material::dc: + return std::vector(n, Sample(1)); + case material::alternating: { + std::vector x(n); + for (std::size_t i = 0; i < n; ++i) { + x[i] = (i % 2 == 0) ? Sample(1) : Sample(-1); + } + return x; + } + } + return {}; + } + + // ------------------------------------------------------------------------ + // Measurement. ulp distance over the ordered-integer image of the bits + // (sign-magnitude folded to two's complement), so +0/-0 count as 1 apart, + // adjacent finite values as 1, and a NaN as "far" rather than as an + // exception. Used only to make a mismatch report readable and to feed the + // informational target; the gate itself is memcmp. + // ------------------------------------------------------------------------ + template + std::uint64_t ulp_distance(Sample a, Sample b) { + using bits_t = std::conditional_t; + static_assert(sizeof(bits_t) == sizeof(Sample)); + bits_t ia = 0; + bits_t ib = 0; + std::memcpy(&ia, &a, sizeof(Sample)); + std::memcpy(&ib, &b, sizeof(Sample)); + if (ia < 0) { + ia = std::numeric_limits::min() - ia; + } + if (ib < 0) { + ib = std::numeric_limits::min() - ib; + } + const std::uint64_t ua = static_cast(ia); + const std::uint64_t ub = static_cast(ib); + return ua > ub ? ua - ub : ub - ua; + } + + struct comparison { + bool identical = true; + std::size_t first_diff = 0; + std::uint64_t max_ulp = 0; + }; + + template + comparison compare(const std::vector& got, const std::vector& ref) { + comparison c; + c.identical = std::memcmp(got.data(), ref.data(), got.size() * sizeof(Sample)) == 0; + if (c.identical) { + return c; + } + bool seen = false; + for (std::size_t i = 0; i < got.size(); ++i) { + const std::uint64_t d = ulp_distance(got[i], ref[i]); + if (d != 0 && !seen) { + c.first_diff = i; + seen = true; + } + c.max_ulp = std::max(c.max_ulp, d); + } + return c; + } + + // Runs one (N, material) pair in both directions. The inverse input is the + // reference forward's output, so both sides receive identical bits and the + // inverse is judged on its own, not through the forward. + template + struct pair_result { + comparison forward; + comparison inverse; + }; + + template + pair_result run_pair(std::size_t n, material m) { + const std::vector x = make_material(m, n); + + ooura_ref ref(n); + engine_under_test dut(n); + std::vector ref_spec = x; + std::vector dut_spec = x; + ref.forward_inplace(ref_spec.data()); + dut.forward_inplace(dut_spec.data()); + + std::vector ref_time = ref_spec; + std::vector dut_time = ref_spec; + ref.inverse_inplace(ref_time.data()); + dut.inverse_inplace(dut_time.data()); + + return {compare(dut_spec, ref_spec), compare(dut_time, ref_time)}; + } + + template + const char* precision_name() { + return sizeof(Sample) == 4 ? "float" : "double"; + } + + /// The sweep: every power of two from 4 to min(65536, TAP_DSP_PARITY_MAX_N). + std::vector sweep_sizes() { + std::vector sizes; + for (std::size_t n = 4; n <= 65536 && n <= static_cast(TAP_DSP_PARITY_MAX_N); n *= 2) { + sizes.push_back(n); + } + return sizes; + } + + // Used by the gate only when it is compiled in (see the #if below). + [[maybe_unused]] constexpr std::size_t k_large_n = std::size_t{1} << 20; + +#if !defined(TAP_DSP_PARITY_INFORMATIONAL) + + // ======================================================================== + // THE GATE: -ffp-contract=off on both sides, memcmp identity. + // ======================================================================== + + template + void expect_identical(std::size_t n, material m, bool forward) { + const pair_result r = run_pair(n, m); + const comparison& c = forward ? r.forward : r.inverse; + ASSERT_TRUE(c.identical) << (forward ? "forward" : "inverse") << " " << precision_name() << " N=" << n + << " material=" << name_of(m) << ": first difference at index " << c.first_diff + << ", max " << c.max_ulp << " ulp. The engine under test is not the same sequence " + << "of IEEE operations as Ooura's rdft; see the statement-fidelity rule in Part 4."; + } + + template + void expect_identical_all_sizes(bool forward) { + for (const std::size_t n : sweep_sizes()) { + for (const material m : k_materials) { + expect_identical(n, m, forward); + if (::testing::Test::HasFatalFailure()) { + return; + } + } + } + } + + TEST(fft_parity_ooura, ForwardIsBitIdenticalToOouraDouble) { + expect_identical_all_sizes(true); + } + + TEST(fft_parity_ooura, ForwardIsBitIdenticalToOouraFloat) { + expect_identical_all_sizes(true); + } + + TEST(fft_parity_ooura, InverseIsBitIdenticalToOouraDouble) { + expect_identical_all_sizes(false); + } + + TEST(fft_parity_ooura, InverseIsBitIdenticalToOouraFloat) { + expect_identical_all_sizes(false); + } + + // One run at 2^20 — the largest geometry any consumer uses (AmbiTap's + // long-partition convolution) — kept as its own test so its runtime shows + // separately. COMPILED OUT, not skipped, when TAP_DSP_PARITY_MAX_N is + // below 2^20: the bare-metal one-shot main counts a GTEST_SKIP as a failed + // gate (tests/bare_metal_main.cpp), and the QEMU legs run at 4096. +#if TAP_DSP_PARITY_MAX_N >= (1 << 20) + template + void expect_identical_large() { + for (const material m : k_materials) { + expect_identical(k_large_n, m, true); + expect_identical(k_large_n, m, false); + if (::testing::Test::HasFatalFailure()) { + return; + } + } + } + + TEST(fft_parity_ooura, LargeTransformIsBitIdenticalToOouraDouble) { + expect_identical_large(); + } + + TEST(fft_parity_ooura, LargeTransformIsBitIdenticalToOouraFloat) { + expect_identical_large(); + } +#endif // TAP_DSP_PARITY_MAX_N >= (1 << 20) + +#else // TAP_DSP_PARITY_INFORMATIONAL + + // ======================================================================== + // INFORMATIONAL: default flags on both sides, measured max ulp per N, + // never a failure. Printed to stdout (visible under `ctest -V`, NOT under + // --output-on-failure, which hides a passing test's output) and recorded + // as gtest properties for the JUnit XML the CMake block requests. One + // test for both precisions, so one process writes the whole XML. + // ======================================================================== + + template + void report_max_ulp() { + std::vector sizes = sweep_sizes(); + if (k_large_n <= static_cast(TAP_DSP_PARITY_MAX_N)) { + sizes.push_back(k_large_n); + } + std::uint64_t overall = 0; + std::printf("[ ulp ] %s: engine under test vs Ooura, default compiler flags (informational)\n", + precision_name()); + std::printf("[ ulp ] %10s %14s %14s %s\n", "N", "forward", "inverse", "worst material (fwd / inv)"); + for (const std::size_t n : sizes) { + std::uint64_t fwd = 0; + std::uint64_t inv = 0; + const char* fwd_worst = name_of(k_materials[0]); + const char* inv_worst = name_of(k_materials[0]); + for (const material m : k_materials) { + const pair_result r = run_pair(n, m); + if (r.forward.max_ulp > fwd) { + fwd = r.forward.max_ulp; + fwd_worst = name_of(m); + } + if (r.inverse.max_ulp > inv) { + inv = r.inverse.max_ulp; + inv_worst = name_of(m); + } + } + overall = std::max({overall, fwd, inv}); + // %lu, not %zu: newlib-nano's printf on the QEMU legs prints "zu". + std::printf("[ ulp ] %10lu %14llu %14llu %s / %s\n", static_cast(n), + static_cast(fwd), static_cast(inv), fwd_worst, + inv_worst); + ::testing::Test::RecordProperty(std::string(precision_name()) + "_N" + std::to_string(n) + "_fwd", + std::to_string(fwd)); + ::testing::Test::RecordProperty(std::string(precision_name()) + "_N" + std::to_string(n) + "_inv", + std::to_string(inv)); + } + std::printf("[ ulp ] %s: max over all N and materials = %llu ulp\n", precision_name(), + static_cast(overall)); + std::fflush(stdout); + SUCCEED() << "informational only: " << overall << " ulp max"; + } + + TEST(fft_parity_ooura_default_flags, ReportsMaxUlpVersusOoura) { + report_max_ulp(); + report_max_ulp(); + } + +#endif // TAP_DSP_PARITY_INFORMATIONAL + +} // namespace diff --git a/tests/test_fft_rt.cpp b/tests/test_fft_rt.cpp new file mode 100644 index 0000000..3258bb8 --- /dev/null +++ b/tests/test_fft_rt.cpp @@ -0,0 +1,405 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// THE REAL-TIME GUARD for tap::dsp::basic_real_fft (Part 9 of +// docs/audit-fft-and-code-smells.md names this file). +// +// fft.h promises that the transforms are noexcept and allocation-free once +// the object is constructed, so they are safe on an audio thread. This file +// 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); +// - "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 +// is program-wide (that is how replaceable allocation functions work), +// so every window below is kept tight: nothing between the two reads of +// the counter but the call under test, and every buffer allocated before +// the first read. +// +// WHAT THE GUARD SEES AND WHAT IT DOES NOT. It counts C++ allocation +// functions only. An engine that allocates through malloc/calloc directly, or +// inside a vendor library (Apple's vDSP on the macOS leg, CMSIS on the M55), +// is invisible to it. For today's Ooura path the claim is complete: makewt +// and makect write into the caller's preallocated ip/w tables (fftsg.c +// 655-756) and the C never calls an allocator. For the backends the guard +// covers the wrapper's own code and nothing more; a malloc interposer +// (glibc's __libc_malloc, or DYLD_INTERPOSE) is the tool for the vendor +// layer and is out of scope here. The self-test CountsAVectorAllocation +// keeps the counter honest: if the replacement ever stopped being picked up, +// that test fails first. +// +// The FIRST call after construction is covered as well as a steady-state one: +// today Ooura builds its trig and bit-reversal tables lazily on the first +// transform (Part 1, item F6), and that initialization must be allocation-free +// too, because the first transform a consumer runs is very often on the audio +// thread already. The port moves table construction into the constructor +// (Part 4); this test is indifferent to where it happens, only to what it +// allocates. +// +// NOT covered, deliberately: the float-I/O convenience overloads on the +// double engine (forward(const float*, float*) / inverse(const float*, +// float*)). fft.h documents them as a setup-time path that allocates a +// staging buffer per call (Part 1, item F7) and Decision D5 deprecates them +// for one consumer cycle; guarding them would pin a behaviour the plan is +// 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 +// 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 +// tables travel with the object, including through operator=. +// +// Stage 4 extension points, per Part 9: the std::span overloads equal the +// pointer overloads, and Engine::is_shareable is what the header says. Neither +// exists yet, so neither is asserted here. + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "support/signals.h" +#include "tap/dsp/fft.h" + +// ---------------------------------------------------------------------------- +// Counting replacements for the replaceable global allocation functions +// ([new.delete.single] and [new.delete.array]). All eight allocating forms +// (throwing / nothrow x plain / aligned x single / array) go through one +// counter; every deallocating form is matched so nothing falls back to the +// library default. Aligned forms are served by over-allocating and stashing +// the raw pointer just below the aligned block, which needs no platform +// aligned-malloc. +// ---------------------------------------------------------------------------- +namespace { + + std::atomic allocation_count{0}; + + void* counted_malloc(std::size_t n) { + allocation_count.fetch_add(1, std::memory_order_relaxed); + return std::malloc(n == 0 ? 1 : n); + } + + void* counted_aligned_malloc(std::size_t n, std::size_t alignment) { + allocation_count.fetch_add(1, std::memory_order_relaxed); + const std::size_t slack = alignment + sizeof(void*); + void* const raw = std::malloc(n + slack); + if (raw == nullptr) { + return nullptr; + } + const std::uintptr_t base = reinterpret_cast(raw) + sizeof(void*); + const std::uintptr_t aligned = (base + alignment - 1) & ~(static_cast(alignment) - 1); + void* const p = reinterpret_cast(aligned); + std::memcpy(reinterpret_cast(aligned) - 1, &raw, sizeof(void*)); + return p; + } + + void counted_aligned_free(void* p) noexcept { + if (p == nullptr) { + return; + } + void* raw = nullptr; + std::memcpy(&raw, static_cast(p) - 1, sizeof(void*)); + std::free(raw); + } + +} // namespace + +// GCC pairs the `new` expressions it inlines in this TU with the replaced +// operator delete below, sees std::free on the result of an "allocation +// function" it treats as opaque, and reports -Wmismatched-new-delete once per +// inlining site (GCC bug 101480: replaced allocation functions that forward to +// malloc/free). The pairing is correct by construction here — every allocating +// form returns std::malloc's result and every deallocating form frees it — so +// the diagnostic is silenced around the replacement functions only. +#if defined(__GNUC__) && !defined(__clang__) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wmismatched-new-delete" +#endif + +void* operator new(std::size_t n) { + void* p = counted_malloc(n); + if (p == nullptr) { + throw std::bad_alloc(); + } + return p; +} +void* operator new[](std::size_t n) { + void* p = counted_malloc(n); + if (p == nullptr) { + throw std::bad_alloc(); + } + return p; +} +void* operator new(std::size_t n, const std::nothrow_t&) noexcept { + return counted_malloc(n); +} +void* operator new[](std::size_t n, const std::nothrow_t&) noexcept { + return counted_malloc(n); +} +void* operator new(std::size_t n, std::align_val_t al) { + void* p = counted_aligned_malloc(n, static_cast(al)); + if (p == nullptr) { + throw std::bad_alloc(); + } + return p; +} +void* operator new[](std::size_t n, std::align_val_t al) { + void* p = counted_aligned_malloc(n, static_cast(al)); + if (p == nullptr) { + throw std::bad_alloc(); + } + return p; +} +void* operator new(std::size_t n, std::align_val_t al, const std::nothrow_t&) noexcept { + return counted_aligned_malloc(n, static_cast(al)); +} +void* operator new[](std::size_t n, std::align_val_t al, const std::nothrow_t&) noexcept { + return counted_aligned_malloc(n, static_cast(al)); +} + +void operator delete(void* p) noexcept { + std::free(p); +} +void operator delete[](void* p) noexcept { + std::free(p); +} +void operator delete(void* p, std::size_t) noexcept { + std::free(p); +} +void operator delete[](void* p, std::size_t) noexcept { + std::free(p); +} +void operator delete(void* p, const std::nothrow_t&) noexcept { + std::free(p); +} +void operator delete[](void* p, const std::nothrow_t&) noexcept { + std::free(p); +} +void operator delete(void* p, std::align_val_t) noexcept { + counted_aligned_free(p); +} +void operator delete[](void* p, std::align_val_t) noexcept { + counted_aligned_free(p); +} +void operator delete(void* p, std::size_t, std::align_val_t) noexcept { + counted_aligned_free(p); +} +void operator delete[](void* p, std::size_t, std::align_val_t) noexcept { + counted_aligned_free(p); +} +void operator delete(void* p, std::align_val_t, const std::nothrow_t&) noexcept { + counted_aligned_free(p); +} +void operator delete[](void* p, std::align_val_t, const std::nothrow_t&) noexcept { + counted_aligned_free(p); +} + +#if defined(__GNUC__) && !defined(__clang__) +#pragma GCC diagnostic pop +#endif + +namespace { + + template + using fft_t = tap::dsp::basic_real_fft; + + // ------------------------------------------------------------------------ + // noexcept, as a compile-time fact on both instantiations. + // ------------------------------------------------------------------------ + 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()))); + 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()); + + // ------------------------------------------------------------------------ + // The allocation guard. + // ------------------------------------------------------------------------ + class allocation_guard { + public: + allocation_guard() noexcept + : m_start(allocation_count.load(std::memory_order_relaxed)) {} + std::size_t allocations_since() const noexcept { + return allocation_count.load(std::memory_order_relaxed) - m_start; + } + + private: + std::size_t m_start; + }; + + constexpr std::size_t k_guarded_sizes[] = {512, 4096}; + + template + class fft_rt_test : public ::testing::Test {}; + + using sample_types = ::testing::Types; + TYPED_TEST_SUITE(fft_rt_test, sample_types); + + // Already proved by the namespace-scope static_asserts above; this exists + // so the promise has a row in the test listing per profile. + TYPED_TEST(fft_rt_test, TransformsAreNoexcept) { + EXPECT_TRUE(transforms_are_noexcept()); + } + + // The guard itself must see allocations, or every test below passes for + // the wrong reason. + TEST(fft_rt_guard, CountsAVectorAllocation) { + allocation_guard guard; + { + std::vector v(64); + EXPECT_NE(v.data(), nullptr); + } + EXPECT_GE(guard.allocations_since(), 1u); + } + + 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; + + // First call after construction (lazy table init today, F6). + { + allocation_guard guard; + call(fft, in.data(), out.data()); + const std::size_t count = guard.allocations_since(); + EXPECT_EQ(count, 0u) << what << " N=" << n << " allocated " << count + << " time(s) on the FIRST call after construction"; + } + // Steady state. + { + allocation_guard guard; + call(fft, in.data(), out.data()); + const std::size_t count = guard.allocations_since(); + EXPECT_EQ(count, 0u) << what << " N=" << n << " allocated " << count << " time(s) in steady state"; + } + } + + TYPED_TEST(fft_rt_test, ForwardInplaceAllocatesNothing) { + 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); }); + } + } + + TYPED_TEST(fft_rt_test, InverseInplaceAllocatesNothing) { + 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); }); + } + } + + TYPED_TEST(fft_rt_test, ForwardOutOfPlaceAllocatesNothing) { + 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); }); + } + } + + TYPED_TEST(fft_rt_test, InverseOutOfPlaceAllocatesNothing) { + 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); }); + } + } + + // ------------------------------------------------------------------------ + // Value semantics: a copy computes exactly what its source computes. + // ------------------------------------------------------------------------ + template + void expect_bit_identical(const std::vector& a, const std::vector& b, const char* what) { + ASSERT_EQ(a.size(), b.size()); + for (std::size_t i = 0; i < a.size(); ++i) { + // memcmp, not ==: the comparison is over the exact bits. + ASSERT_EQ(std::memcmp(&a[i], &b[i], sizeof(Sample)), 0) + << what << " differs from its source at index " << i; + } + } + + template + struct both_directions { + std::vector spectrum; + std::vector time; + }; + + 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()); + return r; + } + + TYPED_TEST(fft_rt_test, CopyProducesBitIdenticalOutput) { + constexpr std::size_t n = 512; + const auto x = tap::dsp::test::random_signal(n, 0x2545F491u); + + fft_t 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"); + + // 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"); + } + + TYPED_TEST(fft_rt_test, CopyOfUnwarmedSourceProducesBitIdenticalOutput) { + constexpr std::size_t n = 512; + 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"); + } + + TYPED_TEST(fft_rt_test, CopyAssignmentProducesBitIdenticalOutput) { + constexpr std::size_t n = 512; + const auto x = tap::dsp::test::random_signal(n, 0x2545F491u); + + fft_t 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); + 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"); + } + +} // namespace