From 57cce940b64021905914a66be0e47e6300f52791 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 12:42:45 +0000 Subject: [PATCH] Adopt the DspTap FIR substrate via submodule MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit SampleRateTap now consumes the shared design/format/kernel/measurement layer from DspTap (submodules/dsptap, pinned at the substrate-extraction merge) instead of carrying it privately — the M0 adoption step of the RatioTap plan, following the MuTap fft.h re-export pattern throughout: - submodules/dsptap + tap::dsp linked into the INTERFACE target; every workflow checkout now fetches submodules recursively. - srt/detail/kaiser.h -> re-export shim (tap::samplerate::detail names unchanged for existing includes). - srt/sample_traits.h -> refinement: each specialization derives from the tap::dsp format core and adds the ASRC-specific blend stratum (make_blend_factor, blend_factor_from_q64, blend); the sample_type concept states the full combined contract. detail::round_sat/ clamp_sat re-exported. - srt/polyphase_filter.h -> dot kernels consumed from tap::dsp::fir_kernels (SRT_ macro names forwarded to the TAP_DSP_ gates); the bank constructor's inline largest-remainder correction replaced by tap::dsp::quantize_row_preserving_sum. - tests/support analysis headers -> shims over tap::dsp::analysis. - Book: 21 anchor includes repointed to the submodule paths (kaiser design, dot kernels, row-sum quantization, moved sample-traits regions, multitone/sine instruments); every remaining include verified to resolve. Behavior change: none. Full suite green under GCC and clang -Werror (73/73 both), clang-tidy clean; the icount CI gates (+-3% on M33/M55/ Hexagon) are the final proof the relocation costs nothing on the hot path. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Ldeq57sBySx2nFTQsGcQB6 --- .github/workflows/book-pages.yml | 2 + .github/workflows/ci-arm64.yml | 2 + .github/workflows/ci.yml | 20 ++ .github/workflows/compare.yml | 2 + .github/workflows/style.yml | 2 + .gitmodules | 3 + CMakeLists.txt | 7 + book/src/epilogue/letter.md | 10 +- book/src/part1/fractional-resampler.md | 4 +- book/src/part1/kaiser.md | 12 +- book/src/part1/sample-traits.md | 10 +- book/src/part3/c3-c5.md | 2 +- book/src/part3/c6.md | 4 +- include/srt/detail/kaiser.h | 371 ++----------------------- include/srt/polyphase_filter.h | 208 +++----------- include/srt/sample_traits.h | 147 +++------- submodules/dsptap | 1 + tests/support/multitone_analysis.h | 229 +-------------- tests/support/sine_analysis.h | 108 +------ 19 files changed, 162 insertions(+), 982 deletions(-) create mode 100644 .gitmodules create mode 160000 submodules/dsptap diff --git a/.github/workflows/book-pages.yml b/.github/workflows/book-pages.yml index 42a6695..a30d14b 100644 --- a/.github/workflows/book-pages.yml +++ b/.github/workflows/book-pages.yml @@ -43,6 +43,8 @@ jobs: url: ${{ steps.deployment.outputs.page_url }} steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Install Doxygen + Graphviz run: sudo apt-get update && sudo apt-get install -y doxygen graphviz diff --git a/.github/workflows/ci-arm64.yml b/.github/workflows/ci-arm64.yml index f9b6684..e284d17 100644 --- a/.github/workflows/ci-arm64.yml +++ b/.github/workflows/ci-arm64.yml @@ -26,6 +26,8 @@ jobs: timeout-minutes: 30 steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Configure (Release) run: cmake -B build -DCMAKE_BUILD_TYPE=Release -DSRT_WERROR=ON diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 037713f..d4dce5b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -39,6 +39,8 @@ jobs: werror: OFF steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Configure env: @@ -68,6 +70,8 @@ jobs: flags: -fsanitize=thread steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Configure env: @@ -109,6 +113,8 @@ jobs: HEXAGON_TOOLCHAIN_SHA256: "55b41922318f6331590ab7baa7f5dbdd99c109327a9c44a52c5e9878fab148c1" steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Cache toolchain id: cache @@ -193,6 +199,8 @@ jobs: timeout-minutes: 30 steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Install toolchain and QEMU run: > @@ -222,6 +230,8 @@ jobs: timeout-minutes: 30 steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Install toolchain and QEMU run: > @@ -294,6 +304,8 @@ jobs: HEXAGON_TOOLCHAIN_SHA256: "55b41922318f6331590ab7baa7f5dbdd99c109327a9c44a52c5e9878fab148c1" steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Install toolchains and QEMU run: > @@ -440,6 +452,8 @@ jobs: timeout-minutes: 15 steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Configure run: > @@ -464,6 +478,8 @@ jobs: timeout-minutes: 20 steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Install dependencies run: > @@ -493,6 +509,8 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive # Run the TapHouse pre-commit hook — the same pinned clang-format # developers get from `pre-commit install` (.pre-commit-config.yaml), so # local and CI share one version AND one file scope (all C/C++ minus @@ -515,6 +533,8 @@ jobs: MDBOOK_SHA256: "9ef07fd288ba58ff3b99d1c94e6d414d431c9a61fdb20348e5beb74b823d546b" steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Install mdBook (pinned) run: | diff --git a/.github/workflows/compare.yml b/.github/workflows/compare.yml index 691c92f..01409cd 100644 --- a/.github/workflows/compare.yml +++ b/.github/workflows/compare.yml @@ -26,6 +26,8 @@ jobs: HEXAGON_TOOLCHAIN_SHA256: "55b41922318f6331590ab7baa7f5dbdd99c109327a9c44a52c5e9878fab148c1" steps: - uses: actions/checkout@df4cb1c069e1874edd31b4311f1884172cec0e10 # v6 + with: + submodules: recursive - name: Install toolchains and QEMU run: > diff --git a/.github/workflows/style.yml b/.github/workflows/style.yml index 9178d13..f6c44ee 100644 --- a/.github/workflows/style.yml +++ b/.github/workflows/style.yml @@ -16,6 +16,8 @@ jobs: timeout-minutes: 30 steps: - uses: actions/checkout@v4 + with: + submodules: recursive - name: Install tools run: sudo apt-get update && sudo apt-get install -y clang-tidy-18 libgtest-dev cmake python3 - name: Configure (compile database) diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..312d995 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "submodules/dsptap"] + path = submodules/dsptap + url = https://github.com/tap/dsptap diff --git a/CMakeLists.txt b/CMakeLists.txt index c84cbcc..ae3bdac 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,12 +1,19 @@ cmake_minimum_required(VERSION 3.24) project(SampleRateTap VERSION 0.1.0 LANGUAGES CXX) +# Shared Tap-family substrate (tap::dsp): kaiser design math, sample-format +# traits, FIR dot kernels, row-sum quantization, measurement instruments. +# Pinned as a submodule per the DspTap release flow (changes land there first, +# consumers bump the pin). +add_subdirectory(submodules/dsptap) + add_library(SampleRateTap INTERFACE) add_library(SampleRateTap::SampleRateTap ALIAS SampleRateTap) target_include_directories(SampleRateTap INTERFACE $ $) target_compile_features(SampleRateTap INTERFACE cxx_std_20) +target_link_libraries(SampleRateTap INTERFACE tap::dsp) if(PROJECT_IS_TOP_LEVEL) option(SRT_BUILD_TESTS "Build SampleRateTap tests" ON) diff --git a/book/src/epilogue/letter.md b/book/src/epilogue/letter.md index 6b1bdde..190c5c5 100644 --- a/book/src/epilogue/letter.md +++ b/book/src/epilogue/letter.md @@ -87,7 +87,7 @@ design target by 1/sinc(f/fs) so the composite comes out flat *and* zeroed: ```cpp -{{#include ../../../include/srt/detail/kaiser.h:pw_comp_design}} +{{#include ../../../submodules/dsptap/include/tap/dsp/kaiser.h:pw_comp_design}} ``` Two implementation details earned their scars. The cosine-series trick is @@ -98,7 +98,7 @@ rect is applied as a running sum *after* windowing, which makes the zeros exact in the discrete composite: ```cpp -{{#include ../../../include/srt/detail/kaiser.h:pw_comp_rect}} +{{#include ../../../submodules/dsptap/include/tap/dsp/kaiser.h:pw_comp_rect}} ``` One bug in between is worth its paragraph, because the test suite caught @@ -158,14 +158,14 @@ log-spaced 60 Hz–16 kHz, amplitudes pink (equal energy per octave — the long-run average of real program), phases chosen for bounded crest: ```cpp -{{#include ../../../tests/support/multitone_analysis.h:pw_comb}} +{{#include ../../../submodules/dsptap/include/tap/dsp/analysis/multitone_analysis.h:pw_comb}} ``` The measurement fits every tone jointly and calls everything left over — aliases, servo FM, noise — the residual: ```cpp -{{#include ../../../tests/support/multitone_analysis.h:pw_metric}} +{{#include ../../../submodules/dsptap/include/tap/dsp/analysis/multitone_analysis.h:pw_metric}} ``` The instrument fought back twice, and both fights are preserved in its @@ -249,7 +249,7 @@ over with a widened tolerance earlier in this very chapter's story. The tricky thing is row-sum-preserving quantization: ```cpp -{{#include ../../../include/srt/polyphase_filter.h:pw_row_sum}} +{{#include ../../../submodules/dsptap/include/tap/dsp/quantize.h:pw_row_sum}} ``` After it, the measured result is the one his claim demanded: a DC input diff --git a/book/src/part1/fractional-resampler.md b/book/src/part1/fractional-resampler.md index bb79b00..55dd26d 100644 --- a/book/src/part1/fractional-resampler.md +++ b/book/src/part1/fractional-resampler.md @@ -72,7 +72,7 @@ The fix is to split the kernel at its natural seam: blend once per frame into a scratch row, then run a plain dot product per channel. ```cpp -{{#include ../../../include/srt/polyphase_filter.h:rs_dot_row}} +{{#include ../../../submodules/dsptap/include/tap/dsp/fir_kernels.h:rs_dot_row}} ``` Two things about this function beyond its arithmetic. First, the comment @@ -326,7 +326,7 @@ history to deliver all channels of tap t contiguously — the frame-major layout — and a register-blocked kernel: ```cpp -{{#include ../../../include/srt/polyphase_filter.h:rs_dot_rows_frame_major}} +{{#include ../../../submodules/dsptap/include/tap/dsp/fir_kernels.h:rs_dot_rows_frame_major}} ``` The measured C6 results, condensed (the full campaign, including the diff --git a/book/src/part1/kaiser.md b/book/src/part1/kaiser.md index f079c3a..072c76a 100644 --- a/book/src/part1/kaiser.md +++ b/book/src/part1/kaiser.md @@ -35,7 +35,7 @@ interpolates exactly. The `sinc` in this file is that function, with the one hazard a numeric programmer would expect handled explicitly: ```cpp -{{#include ../../../include/srt/detail/kaiser.h:kai_sinc}} +{{#include ../../../submodules/dsptap/include/tap/dsp/kaiser.h:kai_sinc}} ``` (The 0/0 at x = 0 is a *removable* singularity — the limit is 1 — but IEEE @@ -114,7 +114,7 @@ I₀(x) = Σₖ [ (x/2)ᵏ / k! ]² which converges for every finite x: ```cpp -{{#include ../../../include/srt/detail/kaiser.h:kai_besseli0}} +{{#include ../../../submodules/dsptap/include/tap/dsp/kaiser.h:kai_besseli0}} ``` Three details carry all the engineering. @@ -155,7 +155,7 @@ constant *relative* accuracy, which is what the window formula's ratio ## `kaiser_beta`: an empirical fit, taken as published ```cpp -{{#include ../../../include/srt/detail/kaiser.h:kai_beta}} +{{#include ../../../submodules/dsptap/include/tap/dsp/kaiser.h:kai_beta}} ``` This is Kaiser's published fit, digit for digit — `0.1102`, `0.5842`, @@ -182,7 +182,7 @@ harris) says taps scale linearly with attenuation and inversely with transition width: ```cpp -{{#include ../../../include/srt/detail/kaiser.h:kai_estimate}} +{{#include ../../../submodules/dsptap/include/tap/dsp/kaiser.h:kai_estimate}} ``` Note what the signature normalizes to: transition width *as a fraction of @@ -222,7 +222,7 @@ UB, however silly the input. ## `design_prototype`: where all of it lands ```cpp -{{#include ../../../include/srt/detail/kaiser.h:kai_prototype}} +{{#include ../../../submodules/dsptap/include/tap/dsp/kaiser.h:kai_prototype}} ``` One pass, one output array, but four decisions are packed into these lines. @@ -292,7 +292,7 @@ was rejected, and since the reasoning is a design artifact it is kept where refactors will trip over it: ```cpp -{{#include ../../../include/srt/detail/kaiser.h:kai_design_note}} +{{#include ../../../submodules/dsptap/include/tap/dsp/kaiser.h:kai_design_note}} ``` Present the alternative fairly, because it *almost* works: diff --git a/book/src/part1/sample-traits.md b/book/src/part1/sample-traits.md index 074b512..929ba22 100644 --- a/book/src/part1/sample-traits.md +++ b/book/src/part1/sample-traits.md @@ -156,7 +156,7 @@ precisely at the row where the response matters most. So the coefficients trade one precision bit for one headroom bit: ```cpp -{{#include ../../../include/srt/sample_traits.h:st_q15_coeff}} +{{#include ../../../submodules/dsptap/include/tap/dsp/sample_traits.h:st_q15_coeff}} ``` with the conversion doing round-half-away-from-zero and saturating at the @@ -164,7 +164,7 @@ integer limits (the *design* is checked separately; saturation here is a belt against future filter specs, not an expected event): ```cpp -{{#include ../../../include/srt/sample_traits.h:st_roundsat}} +{{#include ../../../submodules/dsptap/include/tap/dsp/sample_traits.h:st_roundsat}} ``` What did the traded bit cost? Quantizing coefficients to Q14 puts the @@ -190,7 +190,7 @@ and that is their virtue: `Q15::make_coeff(1.0) == 16384` is the sentence Here is the Q15 multiply-accumulate: ```cpp -{{#include ../../../include/srt/sample_traits.h:st_q15_mac}} +{{#include ../../../submodules/dsptap/include/tap/dsp/sample_traits.h:st_q15_mac}} ``` Two things are chosen here. The product is computed in `int32_t` — a @@ -218,7 +218,7 @@ blended-row rewrite could both be verified *bit-exact* rather than All of the rounding budget is spent in one place: ```cpp -{{#include ../../../include/srt/sample_traits.h:st_q15_finalize}} +{{#include ../../../submodules/dsptap/include/tap/dsp/sample_traits.h:st_q15_finalize}} ``` The accumulator holds a Q29 value (Q0.15 sample × Q1.14 coefficient); the @@ -256,7 +256,7 @@ accumulator worth having on the targets this path exists for. So each product gives up 16 bits *before* joining the sum: ```cpp -{{#include ../../../include/srt/sample_traits.h:st_q31_mac}} +{{#include ../../../submodules/dsptap/include/tap/dsp/sample_traits.h:st_q31_mac}} ``` Now redo the bound: Q45 products have worst-case magnitude 2⁴⁵, and diff --git a/book/src/part3/c3-c5.md b/book/src/part3/c3-c5.md index c084a08..07c4dcb 100644 --- a/book/src/part3/c3-c5.md +++ b/book/src/part3/c3-c5.md @@ -172,7 +172,7 @@ The subtle part of C4 is not the intrinsic; it is the **gate**. Here is the actual block from `include/srt/polyphase_filter.h`, pulled in live: ```cpp -{{#include ../../../include/srt/polyphase_filter.h:opt_smlald_gate}} +{{#include ../../../submodules/dsptap/include/tap/dsp/fir_kernels.h:opt_smlald_gate}} ``` Read the condition: DSP extension present *and MVE absent*. The naive gate diff --git a/book/src/part3/c6.md b/book/src/part3/c6.md index 70ef2b2..e390eef 100644 --- a/book/src/part3/c6.md +++ b/book/src/part3/c6.md @@ -145,7 +145,7 @@ demonstrably keeps in registers. The library's tile is a template on K, taken live from `include/srt/polyphase_filter.h`: ```cpp -{{#include ../../../include/srt/polyphase_filter.h:opt_dot_tile}} +{{#include ../../../submodules/dsptap/include/tap/dsp/fir_kernels.h:opt_dot_tile}} ``` The doc comment carries the 2.8× number in the header itself — the trap @@ -181,7 +181,7 @@ cascade of tiles — 8, then 4, then 2, then 1 — each a specialization of the same template: ```cpp -{{#include ../../../include/srt/polyphase_filter.h:opt_dot_rows}} +{{#include ../../../submodules/dsptap/include/tap/dsp/fir_kernels.h:opt_dot_rows}} ``` Twelve channels is 8+4. Sixteen is 8+8. The deployment shapes exercise diff --git a/include/srt/detail/kaiser.h b/include/srt/detail/kaiser.h index 2ad5891..805fb80 100644 --- a/include/srt/detail/kaiser.h +++ b/include/srt/detail/kaiser.h @@ -1,363 +1,26 @@ +/// \file kaiser.h +/// \brief Re-export of the shared tap::dsp Kaiser prototype designer. // SPDX-License-Identifier: MIT // Copyright 2026 SampleRateTap contributors -// ANCHOR: kai_design_note -/// \file kaiser.h -/// \brief Kaiser-window FIR prototype design for the polyphase interpolation bank. -/// -/// Design note — runtime vs constexpr: the prototype tables run 12K-33K taps and -/// each tap needs sin/sqrt plus a ~50-term Bessel I0 series. Constexpr evaluation -/// is interpreted (roughly 1e3-1e4x slower than native), would need hand-rolled -/// constexpr transcendentals before C++26, and would cost tens of seconds to -/// minutes of compile time in every including translation unit. Runtime design -/// takes well under 10 ms, runs once in a constructor, and is off the audio path, -/// so all design math here is plain runtime double precision. -// ANCHOR_END: kai_design_note +// +// The Kaiser-window FIR design math (bessel_i0, kaiser_beta, estimate_taps, +// sinc, design_prototype, design_prototype_compensated, solve_dense) used to +// live here. It now lives in DspTap (tap::dsp), consumed via the +// submodules/dsptap submodule — one implementation shared with RatioTap and +// the rest of the family. This header keeps the historical include path and +// the unqualified names working inside tap::samplerate::detail. #pragma once -#include -#include -#include -#include -#include +#include "tap/dsp/kaiser.h" namespace tap::samplerate::detail { - // ANCHOR: kai_besseli0 - /// Modified Bessel function of the first kind, order zero, by power series. - /// Converges for all practical Kaiser betas (|x| < ~40); terms are added until - /// they no longer contribute at double precision. - inline double bessel_i0(double x) noexcept { - const double half_x = 0.5 * x; - double term = 1.0; - double sum = 1.0; - for (int k = 1; k < 1000; ++k) { - const double r = half_x / static_cast(k); - term *= r * r; - sum += term; - if (term < 1e-21 * sum) { - break; - } - } - return sum; - } - // ANCHOR_END: kai_besseli0 - - // ANCHOR: kai_beta - /// Kaiser window shape parameter for a given stopband attenuation in dB - /// (Kaiser's published empirical fit). - inline double kaiser_beta(double atten_db) noexcept { - if (atten_db > 50.0) { - return 0.1102 * (atten_db - 8.7); - } - if (atten_db > 21.0) { - return 0.5842 * std::pow(atten_db - 21.0, 0.4) + 0.07886 * (atten_db - 21.0); - } - return 0.0; - } - // ANCHOR_END: kai_beta - - // ANCHOR: kai_estimate - /// Kaiser/harris FIR length estimate, expressed per polyphase branch. - /// - /// \param attenDb target stopband attenuation in dB - /// \param transWidthNorm transition width normalized to the *input* sample rate - /// (e.g. 8 kHz transition at 48 kHz -> 8000/48000) - /// \return estimated taps per polyphase phase: N = (A - 8) / (2.285 * 2*pi * df) - inline std::size_t estimate_taps(double atten_db, double trans_width_norm) noexcept { - // Clamp pathological inputs (attenDb < 8, non-positive width): the raw - // formula goes negative/infinite there and casting that to size_t is UB. - if (!(trans_width_norm > 0.0)) { - return 4; - } - const double n = (atten_db - 8.0) / (2.285 * 2.0 * std::numbers::pi * trans_width_norm); - return n > 4.0 ? static_cast(std::ceil(n)) : 4; - } - // ANCHOR_END: kai_estimate - - // ANCHOR: kai_sinc - /// sin(pi x)/(pi x) with the removable singularity handled. - inline double sinc(double x) noexcept { - if (std::abs(x) < 1e-12) { - return 1.0; - } - const double px = std::numbers::pi * x; - return std::sin(px) / px; - } - // ANCHOR_END: kai_sinc - - // ANCHOR: kai_prototype - /// Designs the Kaiser-windowed sinc prototype lowpass for an L-phase - /// interpolation bank. - /// - /// \param h output, length L*T (the full oversampled prototype) - /// \param numPhases L; the prototype is sampled on a grid of 1/L input samples - /// \param cutoffNorm cutoff normalized to the input Nyquist, i.e. 2*fc/fs in - /// (0, 1]; for a near-unity interpolator centered between a - /// 20 kHz passband and 28 kHz stopband at 48 kHz this is - /// (20000+28000)/48000 = 1.0 (cutoff at input Nyquist) - /// \param beta Kaiser shape parameter (see kaiserBeta) - /// - /// The result is normalized so that sum(h) == L, giving each polyphase branch a - /// DC gain of ~1 (deviation bounded by the stopband leakage). - inline void design_prototype(std::span h, std::size_t num_phases, double cutoff_norm, - double beta) noexcept { - const std::size_t n = h.size(); - const double center = 0.5 * static_cast(n - 1); - const double i0_beta = bessel_i0(beta); - double sum = 0.0; - for (std::size_t i = 0; i < n; ++i) { - const double t = (static_cast(i) - center) / static_cast(num_phases); - const double u = (static_cast(i) - center) / center; // window argument, [-1, 1] - const double w = bessel_i0(beta * std::sqrt(std::max(0.0, 1.0 - u * u))) / i0_beta; - h[i] = cutoff_norm * sinc(cutoff_norm * t) * w; - sum += h[i]; - } - const double gain = static_cast(num_phases) / sum; - for (auto& v : h) { - v *= gain; - } - } - // ANCHOR_END: kai_prototype - - /// Solves the dense n x n system m * out = rhs in place (Gaussian elimination - /// with partial pivoting; row-major m). Small systems only — the compensated - /// design below solves at most 15 unknowns. - inline void solve_dense(std::span m, std::span rhs, std::span out, std::size_t n) noexcept { - std::vector order(n); - for (std::size_t i = 0; i < n; ++i) { - order[i] = i; - } - for (std::size_t col = 0; col < n; ++col) { - std::size_t piv = col; - for (std::size_t r = col + 1; r < n; ++r) { - if (std::abs(m[order[r] * n + col]) > std::abs(m[order[piv] * n + col])) { - piv = r; - } - } - std::swap(order[col], order[piv]); - const std::size_t p = order[col]; - for (std::size_t r = col + 1; r < n; ++r) { - const std::size_t rr = order[r]; - const double f = m[rr * n + col] / m[p * n + col]; - for (std::size_t q = col; q < n; ++q) { - m[rr * n + q] -= f * m[p * n + q]; - } - rhs[rr] -= f * rhs[p]; - } - } - for (std::size_t col = n; col-- > 0;) { - const std::size_t p = order[col]; - double v = rhs[p]; - for (std::size_t q = col + 1; q < n; ++q) { - v -= m[p * n + q] * out[q]; - } - out[col] = v / m[p * n + col]; - } - } - - // ANCHOR: pw_comp_design - /// Designs a prototype with transmission zeros at every integer multiple of - /// the sample rate, passband droop pre-compensated (the music-dsp thread's - /// suggestion, done inside the passband spec — see the epilogue chapter of - /// the book and notebooks/asrc_rbj_analysis.ipynb). - /// - /// Construction: the zeros come from convolving with a one-input-sample rect - /// (multiplies the response by sinc(f/fs), zero at every k*fs — exactly where - /// the images of low-frequency program energy sit). The rect's passband droop - /// (-2.64 dB at 20 kHz) is cancelled by tilting the design target by - /// 1/sinc(f/fs), expressed as a short cosine series in f — which in time is a - /// weighted sum of the brickwall kernel at small integer shifts, so the whole - /// design stays closed-form: no FFT, no dependency. One correction pass - /// (measure the built passband by direct DFT, fold the deviation back into - /// the tilt) holds every preset's ripple within +/-0.005 dB, a >=2x margin - /// on the spec. (A second pass buys ~1 dB of margin for another kernel - /// build plus probe sweep — real money on soft-double targets, where every - /// design flop is a libcall; the pass count and probe count below are - /// sized by the M33 instruction-count ledger in docs/PERFORMANCE.md.) - /// - /// \param h output, length L*T for the TOTAL taps per phase T; the - /// sinc design uses T-1 taps and the rect supplies the - /// +1 (composite length L*(T-1)+L-1, one zero of padding) - /// \param numPhases L - /// \param cutoffNorm as designPrototype - /// \param beta Kaiser shape parameter for the T-1-tap base design - /// \param passbandNorm passband edge / sample rate (flatness is corrected and - /// verified up to here) - /// - /// Costs a few ms more than designPrototype (three kernel builds plus ~100 - /// direct-DFT probes); still constructor-only, off the audio path. Allocates - /// workspace; may throw std::bad_alloc. - // ANCHOR_END: pw_comp_design - inline void design_prototype_compensated(std::span h, std::size_t num_phases, double cutoff_norm, - double beta, double passband_norm) { - const std::size_t L = num_phases; - const std::size_t total = h.size() / L; // total taps per phase (with rect) - const std::size_t td = total - 1; // sinc-design taps per phase - const std::size_t n = L * td; // fine-grid design length - const std::size_t nc = n + L - 1; // composite length after rect - // Compensator order: enough cosine terms to hold the passband tilt to - // ~1e-4, capped so the shifted kernels stay well inside short windows. - const std::size_t M = std::min(14, (td - 1) / 5); - - constexpr std::size_t k_grid = 1001; // fit grid over f/fs in [0, 0.5] - constexpr std::size_t k_probe = 24; // passband correction probes - std::vector target(k_grid), a(M + 1), fine(n), probe(k_probe); - for (std::size_t g = 0; g < k_grid; ++g) { - const double f = 0.5 * static_cast(g) / static_cast(k_grid - 1); - const double pf = std::numbers::pi * f; - target[g] = f < 1e-9 ? 1.0 : pf / std::sin(pf); // 1/sinc(f/fs) - } - - const auto fit_cosine_series = [&] { - // Weighted LS of target on cos(2*pi*m*f): exact where flatness is - // specified (the passband, heavy weight), merely tracked above it. - // Basis by Chebyshev recurrence: cos(m x) from the two previous - // orders, one real cosine per grid point. - std::vector nm((M + 1) * (M + 1), 0.0), rhs(M + 1, 0.0), basis(M + 1); - for (std::size_t g = 0; g < k_grid; ++g) { - const double f = 0.5 * static_cast(g) / static_cast(k_grid - 1); - const double w2 = f <= passband_norm + 0.02 ? 1e8 : 1.0; // (weight 1e4)^2 - const double c1 = std::cos(2.0 * std::numbers::pi * f); - basis[0] = 1.0; - if (M >= 1) { - basis[1] = c1; - } - for (std::size_t m = 2; m <= M; ++m) { - basis[m] = 2.0 * c1 * basis[m - 1] - basis[m - 2]; - } - for (std::size_t r = 0; r <= M; ++r) { - for (std::size_t q = 0; q <= M; ++q) { - nm[r * (M + 1) + q] += w2 * basis[r] * basis[q]; - } - rhs[r] += w2 * basis[r] * target[g]; - } - } - solve_dense(nm, rhs, a, M + 1); - }; - - const auto build = [&] { - // Tilted ideal kernel: sum of the brickwall sinc at integer shifts. - // Centered at n/2 (not (n-1)/2): the even-length rect below shifts - // the composite by (L-1)/2, and n/2 + (L-1)/2 == (L*total - 1)/2 — - // the exact center designPrototype uses, so the bank's phase/delay - // convention is identical for both designs. (Getting this wrong is a - // half-fine-sample delay error: ~-72 dB at 1 kHz, worse by 6 dB per - // octave — the fractional-delay accuracy tests catch it.) - // - // Transcendental budget: the naive form of this loop calls libm sin - // once per tap per compensator term (~2M calls across the design; a - // measured +225M constructor instructions on Cortex-M55, and worse - // where doubles are soft). Instead: sin(pi*c*(t -+ m)) expands by - // angle addition over precomputed sin/cos(pi*c*m), and sin/cos of - // the per-tap angle advance by a unit rotator, re-synced with real - // libm calls every 4096 taps to bound drift far below the design's - // own accuracy floor. - const double center = 0.5 * static_cast(n); - const double i0_beta = bessel_i0(beta); - std::vector cs(M + 1), sn(M + 1); - for (std::size_t m = 0; m <= M; ++m) { - cs[m] = std::cos(std::numbers::pi * cutoff_norm * static_cast(m)); - sn[m] = std::sin(std::numbers::pi * cutoff_norm * static_cast(m)); - } - const double step = std::numbers::pi * cutoff_norm / static_cast(L); - const double step_c = std::cos(step), step_s = std::sin(step); - double ang_s = 0.0, ang_c = 1.0; // sin/cos(pi*c*t_i), re-synced below - for (std::size_t i = 0; i < n; ++i) { - const double t = (static_cast(i) - center) / static_cast(L); - if (i % 4096 == 0) { - ang_s = std::sin(std::numbers::pi * cutoff_norm * t); - ang_c = std::cos(std::numbers::pi * cutoff_norm * t); - } - const auto shifted_sinc = [&](double dm, double sin_shift, double cos_shift) { - const double x = cutoff_norm * (t - dm); // dm may be negative - if (std::abs(x) < 1e-12) { - return 1.0; - } - // sin(pi*c*(t - dm)) = sin(pi*c*t)cos(pi*c*dm) - cos(..)sin(..) - return (ang_s * cos_shift - ang_c * sin_shift) / (std::numbers::pi * x); - }; - double v = a[0] * cutoff_norm * shifted_sinc(0.0, 0.0, 1.0); - for (std::size_t m = 1; m <= M; ++m) { - const double dm = static_cast(m); - v += 0.5 * a[m] * cutoff_norm * (shifted_sinc(dm, sn[m], cs[m]) + shifted_sinc(-dm, -sn[m], cs[m])); - } - const double u = (static_cast(i) - center) / center; - fine[i] = v * bessel_i0(beta * std::sqrt(std::max(0.0, 1.0 - u * u))) / i0_beta; - const double next_s = ang_s * step_c + ang_c * step_s; - ang_c = ang_c * step_c - ang_s * step_s; - ang_s = next_s; - } - // ANCHOR: pw_comp_rect - // Rect convolution as a running sum: exact zeros at every k*fs. - double run = 0.0; - for (std::size_t i = 0; i < nc; ++i) { - run += i < n ? fine[i] : 0.0; - if (i >= L) { - run -= fine[i - L]; - } - h[i] = run / static_cast(L); - } - for (std::size_t i = nc; i < h.size(); ++i) { - h[i] = 0.0; - } - // ANCHOR_END: pw_comp_rect - double sum = 0.0; - for (std::size_t i = 0; i < nc; ++i) { - sum += h[i]; - } - const double gain = static_cast(L) / sum; - for (std::size_t i = 0; i < nc; ++i) { - h[i] *= gain; - } - }; - - for (int pass = 0; pass < 1; ++pass) { - fit_cosine_series(); - build(); - // Probe the built passband by direct DFT (cos projection about the - // composite's symmetry center, (L*total - 1)/2 == nc/2) and fold the - // deviation into the tilt. One rotator per probe frequency: two libm - // calls each instead of one per tap (rotator drift over ~2^14 steps - // is ~1e-12, five orders below the ripple being measured). - const double center = 0.5 * static_cast(nc); - for (std::size_t j = 0; j < k_probe; ++j) { - const double f = passband_norm * static_cast(j + 1) / k_probe; - const double th = 2.0 * std::numbers::pi * f / static_cast(L); - const double th_c = std::cos(th), th_s = std::sin(th); - double rc = std::cos(th * -center), rs = std::sin(th * -center); - double acc = 0.0; - for (std::size_t i = 0; i < nc; ++i) { - acc += h[i] * rc; - const double nrc = rc * th_c - rs * th_s; - rs = rs * th_c + rc * th_s; - rc = nrc; - } - probe[j] = std::abs(acc) / static_cast(L); - } - for (std::size_t g = 0; g < k_grid; ++g) { - const double f = 0.5 * static_cast(g) / static_cast(k_grid - 1); - if (f > passband_norm) { - continue; - } - // probe[j] sits at f = passbandNorm*(j+1)/kProbe, i.e. x = j+1 - const double x = f / passband_norm * k_probe - 1.0; - double d; - if (x <= 0.0) { - d = probe[0]; - } - else if (x >= static_cast(k_probe - 1)) { - d = probe[k_probe - 1]; - } - else { - const auto j = static_cast(x); - const double fr = x - static_cast(j); - d = probe[j] * (1.0 - fr) + probe[j + 1] * fr; - } - target[g] /= std::max(d, 0.5); - } - } - fit_cosine_series(); - build(); - } + using tap::dsp::bessel_i0; + using tap::dsp::design_prototype; + using tap::dsp::design_prototype_compensated; + using tap::dsp::estimate_taps; + using tap::dsp::kaiser_beta; + using tap::dsp::sinc; + using tap::dsp::solve_dense; } // namespace tap::samplerate::detail diff --git a/include/srt/polyphase_filter.h b/include/srt/polyphase_filter.h index bc11b00..69f18e2 100644 --- a/include/srt/polyphase_filter.h +++ b/include/srt/polyphase_filter.h @@ -8,60 +8,27 @@ #include #include #include +#include #include #include #include #include "srt/detail/kaiser.h" #include "srt/sample_traits.h" - -// No-alias qualifier for the kernel hot loops: without it the compiler -// versions the blend loop behind a runtime aliasing check (verified with -// -fopt-info-vec; see docs/PERFORMANCE.md, hypothesis 2). -#if defined(_MSC_VER) -#define SRT_RESTRICT __restrict -#else -#define SRT_RESTRICT __restrict__ -#endif - -// ANCHOR: opt_smlald_gate -// Dual 16x16 MAC (SMLALD) for the Q15 dot product on Arm cores that have -// the DSP extension but no Helium — the Cortex-M33/M4/M7 class (e.g. -// Raspberry Pi Pico 2). Gated off when MVE is present: on M55 the compiler -// already auto-vectorizes the scalar loop with Helium and the intrinsic -// path would replace vectors with dual-MACs (see docs/PERFORMANCE.md, -// hypothesis 4). Bit-exactness: each 16x16 product is exact in int32 and -// the int64 accumulation is associative, so pairing changes no output bit. -#if defined(__ARM_FEATURE_DSP) && !defined(__ARM_FEATURE_MVE) -#include -#define SRT_Q15_SMLALD 1 -#else -#define SRT_Q15_SMLALD 0 -#endif -// ANCHOR_END: opt_smlald_gate - -// Channel-parallel dot product for high channel counts (hypothesis C6, -// docs/PERFORMANCE.md): history stored frame-major so the per-tap inner -// loop runs across channels — contiguous loads, one accumulator lane per -// channel, coefficient broadcast. Bit-exact because each channel's -// accumulation order over taps is unchanged (lanes are channels, not -// taps), which is what lets the FLOAT path vectorize at all: its strict -// per-channel double accumulation forbids tap-axis SIMD (PERFORMANCE.md -// hypothesis 5), but the channel axis is free. Float-only by measurement: -// fixed-point planar dots already auto-vectorize over taps on hosts -// (integer reduction is exactly reassociable) and measured ~1.5x FASTER -// than the channel-parallel form. Host-only: the embedded targets keep -// their proven planar codegen (Helium on M55, SMLALD on M33-class, -// Hexagon's measured scalar floor — see PERFORMANCE.md C4/C5). -#if !defined(__ARM_FEATURE_MVE) && !defined(__ARM_FEATURE_DSP) && !defined(__hexagon__) -#define SRT_CHANNEL_PARALLEL 1 -#else -#define SRT_CHANNEL_PARALLEL 0 -#endif -// Minimum channel count for the frame-major path (overridable for A/B -// measurements; the blend-share planar path stays better at low counts). +#include "tap/dsp/fir_kernels.h" +#include "tap/dsp/quantize.h" + +// The kernel hot loops (dot_row, the channel-parallel tiles) and their target +// gates (SMLALD dual-MAC on DSP-extension Arm, the channel-parallel layout +// selection) moved to DspTap's fir_kernels.h with the rest of the shared FIR +// substrate; the measurements that sized them are recorded there and in +// docs/PERFORMANCE.md (hypotheses 2/4/5/C4-C6). The historical SRT_ macro +// names keep working, forwarded to the TAP_DSP_ gates. +#define SRT_RESTRICT TAP_DSP_RESTRICT +#define SRT_Q15_SMLALD TAP_DSP_Q15_SMLALD +#define SRT_CHANNEL_PARALLEL TAP_DSP_CHANNEL_PARALLEL #ifndef SRT_CP_MIN_CHANNELS -#define SRT_CP_MIN_CHANNELS 4 +#define SRT_CP_MIN_CHANNELS TAP_DSP_CP_MIN_CHANNELS #endif namespace tap::samplerate { @@ -191,48 +158,19 @@ namespace tap::samplerate { } m_table.resize((m_phases + 1) * m_taps); - std::vector remainder(m_taps); + // Row-sum-preserving quantization ("the coefficients of every + // phase must add to one" — R. Bristow-Johnson, music-dsp): the + // largest-remainder correction now lives in DspTap's quantize.h, + // shared with RatioTap's fixed-ratio tables. Gather each branch + // into storage (tap-reversed) order in double, then quantize the + // row as a unit so its DC sum survives fixed point exactly. + std::vector row_d(m_taps); for (std::size_t p = 0; p < m_phases; ++p) { - coeff* row = m_table.data() + p * m_taps; - double exact_sum = 0.0; - std::int64_t quant_sum = 0; for (std::size_t t = 0; t < m_taps; ++t) { - const std::size_t m = t * m_phases + p; // prototype index of (branch p, tap t) - const double v = (m < n) ? proto[m] : 0.0; - const double scaled = v * sample_traits::k_coeff_scale; - const coeff q = sample_traits::make_coeff(v); - row[m_taps - 1 - t] = q; - remainder[m_taps - 1 - t] = scaled - static_cast(q); - exact_sum += scaled; - quant_sum += static_cast(q); - } - if constexpr (!std::is_floating_point_v) { - // ANCHOR: pw_row_sum - // Row-sum-preserving quantization ("the coefficients of every - // phase must add to one" — R. Bristow-Johnson, music-dsp). In - // double the image_zeros designs make every branch's DC sum - // identical to machine epsilon (the k*fs zeros ARE branch-DC - // uniformity, stated in frequency), but independent per-tap - // rounding rebroke it by several LSB. Distribute each row's - // total rounding residual to the taps that were rounded - // furthest from it (largest-remainder method): every row then - // sums to llround(exact * scale), so DC gain deviates by at - // most one coefficient LSB across all mu. - std::int64_t residual = static_cast(std::llround(exact_sum)) - quant_sum; - while (residual != 0) { - const double sgn = residual > 0 ? 1.0 : -1.0; - std::size_t best = 0; - for (std::size_t u = 1; u < m_taps; ++u) { - if (sgn * remainder[u] > sgn * remainder[best]) { - best = u; - } - } - row[best] = static_cast(row[best] + (residual > 0 ? 1 : -1)); - remainder[best] -= sgn; - residual -= residual > 0 ? 1 : -1; - } - // ANCHOR_END: pw_row_sum + const std::size_t m = t * m_phases + p; // prototype index of (branch p, tap t) + row_d[m_taps - 1 - t] = (m < n) ? proto[m] : 0.0; } + tap::dsp::quantize_row_preserving_sum(row_d, std::span(m_table.data() + p * m_taps, m_taps)); } // The extra row L is row 0 advanced one input sample; copy it from // the QUANTIZED row 0 so the mu-wrap continuity invariant holds @@ -356,95 +294,15 @@ namespace tap::samplerate { } // ANCHOR_END: rs_interpolate_phase - // ANCHOR: rs_dot_row - /// Dot product of a pre-blended coefficient row against a history window. - /// Identical arithmetic to interpolate() given the same mu: blend then mac, - /// per tap, in the same order — outputs are bit-exact either way. - template - inline S dot_row(const typename sample_traits::coeff* SRT_RESTRICT row, const S* SRT_RESTRICT hist, - std::size_t taps) noexcept { - using tr = sample_traits; -#if SRT_Q15_SMLALD - if constexpr (std::is_same_v) { - std::int64_t acc = 0; - std::size_t t = 0; - for (; t + 1 < taps; t += 2) { - // memcpy keeps the 16-bit pair loads alignment-safe; both - // compile to a single 32-bit load (little-endian packing - // matches SMLALD's lo/hi lanes). - std::uint32_t h; - std::uint32_t r; - std::memcpy(&h, hist + t, sizeof h); - std::memcpy(&r, row + t, sizeof r); - acc = __smlald(static_cast(h), static_cast(r), acc); - } - for (; t < taps; ++t) // odd-tap tail; every preset is even - acc = tr::mac(acc, hist[t], row[t]); - return tr::finalize(acc); - } -#endif - typename tr::accum acc{}; - for (std::size_t t = 0; t < taps; ++t) { - acc = tr::mac(acc, hist[t], row[t]); - } - return tr::finalize(acc); - } - // ANCHOR_END: rs_dot_row - - // ANCHOR: opt_dot_tile - /// One K-channel tile of the channel-parallel dot (hypothesis C6): K - /// accumulators live in a constexpr-size local array — registers, not - /// memory — while the tap loop walks the frame-major window with stride - /// `stride` samples per frame. K is the register-blocking factor; a naive - /// channels-inner loop with accumulators in memory measures ~2.8x SLOWER - /// than planar (each mac round-trips its accumulator through the stack). - template - inline void dot_tile_frame_major(const typename sample_traits::coeff* SRT_RESTRICT row, const S* SRT_RESTRICT x, - std::size_t taps, std::size_t stride, S* SRT_RESTRICT out) noexcept { - using tr = sample_traits; - typename tr::accum acc[K]{}; - for (std::size_t t = 0; t < taps; ++t) { - const auto coeff = row[t]; - const S* SRT_RESTRICT frame = x + t * stride; - for (std::size_t k = 0; k < K; ++k) { - acc[k] = tr::mac(acc[k], frame[k], coeff); - } - } - for (std::size_t k = 0; k < K; ++k) { - out[k] = tr::finalize(acc[k]); - } - } - // ANCHOR_END: opt_dot_tile - - // ANCHOR: rs_dot_rows_frame_major - // ANCHOR: opt_dot_rows - /// Channel-parallel dot products over a frame-major history block: all - /// channels' outputs for one frame in register-blocked tiles of 8/4/2/1. - /// Per channel the accumulation order over taps equals dotRow's, so the - /// outputs are bit-exact vs the planar path for every sample type — float - /// included, since each channel's double accumulator still sums the taps - /// in the same order (lanes are channels, not taps). - template - inline void dot_rows_frame_major(const typename sample_traits::coeff* SRT_RESTRICT row, const S* SRT_RESTRICT x, - std::size_t taps, std::size_t channels, S* SRT_RESTRICT out) noexcept { - std::size_t c = 0; - for (; c + 8 <= channels; c += 8) { - dot_tile_frame_major(row, x + c, taps, channels, out + c); - } - if (c + 4 <= channels) { - dot_tile_frame_major(row, x + c, taps, channels, out + c); - c += 4; - } - if (c + 2 <= channels) { - dot_tile_frame_major(row, x + c, taps, channels, out + c); - c += 2; - } - if (c < channels) { - dot_tile_frame_major(row, x + c, taps, channels, out + c); - } - } - // ANCHOR_END: rs_dot_rows_frame_major - // ANCHOR_END: opt_dot_rows + // The dot-product kernels (dot_row with its SMLALD path, the + // channel-parallel dot_tile_frame_major / dot_rows_frame_major) moved to + // DspTap's fir_kernels.h with the rest of the shared FIR substrate — + // identical code, and the tap::dsp core traits they run on are the base + // of this library's traits, so outputs are bit-exact. Re-exported here so + // existing unqualified uses keep working. + using tap::dsp::dot_row; + using tap::dsp::dot_rows_frame_major; + using tap::dsp::dot_tile_frame_major; // ANCHOR: rs_class_doc /// Streaming fractional-delay engine for one converter instance. diff --git a/include/srt/sample_traits.h b/include/srt/sample_traits.h index fd14b4d..e45eee1 100644 --- a/include/srt/sample_traits.h +++ b/include/srt/sample_traits.h @@ -14,6 +14,14 @@ /// - std::int32_t : Q31 samples, Q1.30 coefficients, int64 accumulation, /// saturating output /// +/// The format core — coefficient Q formats, mac/finalize, rounding and +/// saturation — lives in DspTap (tap::dsp::sample_traits, consumed via the +/// submodules/dsptap submodule) and is shared with RatioTap and the rest of +/// the family. This header layers the ASRC-specific stratum on top: the +/// inter-phase coefficient *blend* machinery (mu interpolation between +/// adjacent polyphase rows, including the Q0.64 phase-accumulator entry +/// points), which a fixed-ratio converter has no use for. +/// /// The clock servo and the filter design always run in double regardless of /// sample type (control path and one-time init, not the audio path), so the /// fixed-point datapaths contain no floating-point inner loops. @@ -25,59 +33,33 @@ #include #include +#include "tap/dsp/sample_traits.h" + namespace tap::samplerate { namespace detail { - // ANCHOR: st_roundsat - /// Round-and-saturate a double to a signed integer coefficient/sample type. - template - constexpr I round_sat(double v) noexcept { - constexpr double lo = static_cast(std::numeric_limits::min()); - constexpr double hi = static_cast(std::numeric_limits::max()); - const double r = v < 0.0 ? v - 0.5 : v + 0.5; // round half away from zero - if (r <= lo) { - return std::numeric_limits::min(); - } - if (r >= hi) { - return std::numeric_limits::max(); - } - return static_cast(r); - } - // ANCHOR_END: st_roundsat - - /// Saturate a 64-bit accumulator result to a narrower signed integer. - template - constexpr I clamp_sat(std::int64_t v) noexcept { - constexpr auto lo = static_cast(std::numeric_limits::min()); - constexpr auto hi = static_cast(std::numeric_limits::max()); - return static_cast(v < lo ? lo : (v > hi ? hi : v)); - } + // Historical homes of the round/saturate helpers; the implementations + // moved to DspTap with the format core. + using tap::dsp::detail::clamp_sat; + using tap::dsp::detail::round_sat; } // namespace detail // ANCHOR: st_primary /// Primary template intentionally undefined; specialize per sample type. + /// Each specialization derives from the shared tap::dsp format core and + /// adds the blend stratum. template struct sample_traits; // ANCHOR_END: st_primary // ANCHOR: st_float - /// Float datapath: float samples and coefficients, double accumulation. - /// The double accumulator keeps the dot-product noise floor far below the - /// 120 dB transparency target; float coefficient storage quantizes the - /// filter at roughly -150 dB, negligible against the same target. + /// Float datapath: the tap::dsp float core (double accumulation) plus a + /// single-precision blend stratum. template <> - struct sample_traits { - using coeff = float; ///< stored filter coefficient type - using accum = double; ///< dot-product accumulator type - using blend_factor = float; ///< per-sample fractional blend representation - - /// Convert a double-precision designed coefficient to storage form. - static coeff make_coeff(double c) noexcept { return static_cast(c); } - /// Coefficient units per 1.0 (used by the bank's row-sum-preserving - /// quantization; unity for floating storage, where no correction runs). - static constexpr double k_coeff_scale = 1.0; + struct sample_traits : tap::dsp::sample_traits { + using blend_factor = float; ///< per-sample fractional blend representation /// Convert the intra-phase fraction (in [0,1)) once per output sample. static blend_factor make_blend_factor(double fr) noexcept { return static_cast(fr); } @@ -92,47 +74,20 @@ namespace tap::samplerate { } // ANCHOR_END: st_blend_q64_float - /// acc + x * c, in the accumulator domain. - static accum mac(accum acc, float x, coeff c) noexcept { - return acc + static_cast(x) * static_cast(c); - } - /// Linear blend between two adjacent-phase coefficients. static coeff blend(coeff a, coeff b, blend_factor fr) noexcept { return a + fr * (b - a); } - - /// Convert the accumulator to an output sample (saturates for fixed point). - static float finalize(accum acc) noexcept { return static_cast(acc); } - - /// The zero/silence sample value. - static float silence() noexcept { return 0.0f; } }; // ANCHOR_END: st_float // ANCHOR: st_q15_header - /// Q15 fixed-point datapath (samples are int16_t in Q0.15). - /// - /// Coefficients are stored in Q1.14: the prototype's peak tap reaches ~1.0 - /// (per-phase DC gain is 1), which does not fit Q0.15, so one headroom bit - /// is traded for one precision bit. Products are Q0.15 x Q1.14 = Q29 and are - /// summed exactly in int64 (48-80 taps add ~6-7 bits — no overflow, no - /// intermediate rounding). The single rounding happens in finalize(): - /// Q29 -> Q15 with round-half-up and saturation. Coefficient quantization - /// (Q14, ~-86 dB) and output quantization (Q15) set the noise floor — both - /// at the format's own limit, so the converter is Q15-transparent. + /// Q15 fixed-point datapath: the tap::dsp Q15 core (Q1.14 coefficients, + /// exact int64 accumulation, single Q29 -> Q15 rounding) plus an + /// integer-only blend stratum with the fraction carried in Q15. template <> - struct sample_traits { - using coeff = std::int16_t; - using accum = std::int64_t; + struct sample_traits : tap::dsp::sample_traits { using blend_factor = std::int32_t; ///< fraction in Q15 // ANCHOR_END: st_q15_header - // ANCHOR: st_q15_coeff - static coeff make_coeff(double c) noexcept { - return detail::round_sat(c * 16384.0); // Q1.14 - } - static constexpr double k_coeff_scale = 16384.0; // Q1.14 units per 1.0 - // ANCHOR_END: st_q15_coeff - static blend_factor make_blend_factor(double fr) noexcept { return static_cast(fr * 32768.0); // Q15 } @@ -145,12 +100,6 @@ namespace tap::samplerate { } // ANCHOR_END: st_q15_q64 - // ANCHOR: st_q15_mac - static accum mac(accum acc, std::int16_t x, coeff c) noexcept { - return acc + static_cast(static_cast(x) * static_cast(c)); - } - // ANCHOR_END: st_q15_mac - // ANCHOR: st_q15_blend static coeff blend(coeff a, coeff b, blend_factor fr) noexcept { // Q14 + (Q15 * Q14) >> 15, in int64: the worst-case int32 product @@ -162,39 +111,17 @@ namespace tap::samplerate { return static_cast(a + ((fr * diff) >> 15)); } // ANCHOR_END: st_q15_blend - - // ANCHOR: st_q15_finalize - static std::int16_t finalize(accum acc) noexcept { - // Round-half-up, not half-even: the bias is a fraction of one - // sub-LSB rounding step, far below the Q15 noise floor. - return detail::clamp_sat((acc + (1 << 13)) >> 14); // Q29 -> Q15 - } - // ANCHOR_END: st_q15_finalize - - static std::int16_t silence() noexcept { return 0; } }; // ANCHOR: st_q31 - /// Q31 fixed-point datapath (samples are int32_t in Q0.31). - /// - /// Coefficients are stored in Q1.30 (one headroom bit for the ~1.0 peak - /// tap). A full-precision product would be Q0.31 x Q1.30 = 62 bits, which - /// overflows int64 once ~48 of them are summed, so each product is - /// pre-shifted down 16 bits (Q45) before accumulation; the discarded bits - /// sit 14 bits below the final Q31 LSB, far beneath the format's noise - /// floor. finalize() rounds Q45 -> Q31 with saturation. The blend fraction - /// uses Q20 (the int64 blend path makes the extra precision free). + /// Q31 fixed-point datapath: the tap::dsp Q31 core (Q1.30 coefficients, + /// products pre-shifted to Q45) plus an integer-only blend stratum. The + /// blend fraction uses Q20 (the int64 blend path makes the extra + /// precision free). template <> - struct sample_traits { - using coeff = std::int32_t; - using accum = std::int64_t; + struct sample_traits : tap::dsp::sample_traits { using blend_factor = std::int32_t; ///< fraction in Q20 - static coeff make_coeff(double c) noexcept { - return detail::round_sat(c * 1073741824.0); // Q1.30 - } - static constexpr double k_coeff_scale = 1073741824.0; // Q1.30 units per 1.0 - static blend_factor make_blend_factor(double fr) noexcept { return static_cast(fr * 1048576.0); // Q20 } @@ -204,28 +131,18 @@ namespace tap::samplerate { return static_cast(frac >> 44); // Q20 } - // ANCHOR: st_q31_mac - static accum mac(accum acc, std::int32_t x, coeff c) noexcept { - return acc + ((static_cast(x) * c) >> 16); // Q61 -> Q45 - } - // ANCHOR_END: st_q31_mac - static coeff blend(coeff a, coeff b, blend_factor fr) noexcept { const std::int64_t diff = static_cast(b) - a; return static_cast(a + ((fr * diff) >> 20)); } - - static std::int32_t finalize(accum acc) noexcept { - return detail::clamp_sat((acc + (1 << 13)) >> 14); // Q45 -> Q31 - } - - static std::int32_t silence() noexcept { return 0; } }; // ANCHOR_END: st_q31 // ANCHOR: st_concept /// Satisfied by any type with a complete, well-formed sample_traits - /// specialization. + /// specialization: the tap::dsp format core (make_coeff/mac/finalize/ + /// silence — see tap::dsp::sample_type) refined with this library's blend + /// contract. template concept sample_type = requires(T x, double d, typename sample_traits::accum a, typename sample_traits::coeff c, typename sample_traits::blend_factor f) { diff --git a/submodules/dsptap b/submodules/dsptap new file mode 160000 index 0000000..98abb04 --- /dev/null +++ b/submodules/dsptap @@ -0,0 +1 @@ +Subproject commit 98abb04733c3f171276d07059bb3fe25a1fe7474 diff --git a/tests/support/multitone_analysis.h b/tests/support/multitone_analysis.h index bd31ca9..055573c 100644 --- a/tests/support/multitone_analysis.h +++ b/tests/support/multitone_analysis.h @@ -1,226 +1,21 @@ // SPDX-License-Identifier: MIT // Copyright 2026 SampleRateTap contributors -// Program-weighted multitone metric: a deterministic pink-weighted tone -// comb and the residual measurement that scores a converter against it. -// -// Why this exists: every other quality number in this suite is a single -// sine — the worst-case metric. A filter with transmission zeros at k*fs -// (filter_spec::imageZeros) is deliberately optimized for a different -// promise: alias rejection weighted by where real program energy lives -// (predominantly the bottom octaves). That promise is unverifiable with -// single sines, so this header supplies the instrument: K log-spaced tones -// with pink (equal-energy-per-octave) amplitudes, and a fit-subtract -// residual over the converted tail. See the book's epilogue chapter. +// Re-export of the shared tap::dsp::analysis program-weighted multitone +// instrument (pink tone comb + joint LS fit). The implementation moved to +// DspTap with the measurement harness; this header keeps the historical +// include path and the srt_test names working. The design history is the +// book's epilogue chapter. #pragma once -#include -#include -#include -#include -#include -#include - -#include "srt/detail/kaiser.h" // detail::solveDense for the joint LS -#include "support/sine_analysis.h" +#include "tap/dsp/analysis/multitone_analysis.h" namespace srt_test { - // ANCHOR: pw_comb - struct tone_comb { - std::vector freq_hz; // log-spaced - std::vector amplitude; // pink: a ~ 1/sqrt(f), scaled to peakSum - std::vector phase; // golden-angle sequence: bounded crest - - /// K tones from fLo to fHi; sum of amplitudes == peakSum, so the summed - /// signal can never exceed peakSum even in the worst phase alignment. - static tone_comb pink(std::size_t k, double f_lo, double f_hi, double peak_sum) { - tone_comb c; - double sum = 0.0; - for (std::size_t i = 0; i < k; ++i) { - const double f = f_lo * std::pow(f_hi / f_lo, static_cast(i) / static_cast(k - 1)); - c.freq_hz.push_back(f); - c.amplitude.push_back(1.0 / std::sqrt(f / f_lo)); - c.phase.push_back(2.0 * std::numbers::pi * 0.6180339887498949 * static_cast(i * i)); - sum += c.amplitude.back(); - } - for (auto& a : c.amplitude) { - a *= peak_sum / sum; - } - return c; - } - - /// Sample of the comb at input sample index i (rate fs). - double sample_at(std::uint64_t i, double fs) const { - double v = 0.0; - for (std::size_t k = 0; k < freq_hz.size(); ++k) { - v += amplitude[k] - * std::sin(2.0 * std::numbers::pi * freq_hz[k] / fs * static_cast(i) + phase[k]); - } - return v; - } - }; - // ANCHOR_END: pw_comb - - /// Fits a*sin + b*cos at a fixed normalized frequency (no DC term; the - /// joint fit models DC), returning the fitted component's power. - struct tone_fit { - double a = 0.0, b = 0.0; - double power() const { return 0.5 * (a * a + b * b); } - }; - - inline tone_fit fit_tone_fixed(std::span x, double freq_norm) { - const double w = 2.0 * std::numbers::pi * freq_norm; - double ss = 0.0, sc = 0.0, cc = 0.0, rs = 0.0, rc = 0.0; - for (std::size_t i = 0; i < x.size(); ++i) { - const double s = std::sin(w * static_cast(i)); - const double c = std::cos(w * static_cast(i)); - ss += s * s; - sc += s * c; - cc += c * c; - rs += s * x[i]; - rc += c * x[i]; - } - const double det = ss * cc - sc * sc; - tone_fit f; - f.a = (rs * cc - rc * sc) / det; - f.b = (rc * ss - rs * sc) / det; - return f; - } - - /// Refines a tone's frequency by comparing the fitted phase of the two - /// window halves (fitSineTracked's method, on the double work buffer and - /// without a DC term), returning the refined normalized frequency. - inline double track_tone_freq(std::span x, double freq_norm) { - double f = freq_norm; - const std::size_t half = x.size() / 2; - for (int iter = 0; iter < 4; ++iter) { - const tone_fit a = fit_tone_fixed(x.first(half), f); - const tone_fit b = fit_tone_fixed(x.subspan(half), f); - const double two_pi = 2.0 * std::numbers::pi; - const double predicted = std::atan2(a.b, a.a) + two_pi * f * static_cast(half); - const double dphi = std::remainder(std::atan2(b.b, b.a) - predicted, two_pi); - f += dphi / (two_pi * static_cast(half)); - } - return f; - } - - // ANCHOR: pw_metric - /// Joint least-squares fit of all tones at once (2K unknowns via normal - /// equations), writing per-tone fits and returning the residual out-of-model - /// power. Sequential fit-subtract is NOT enough here: 24 tones on a - /// rectangular window leak into each other far above the -120 dB floors - /// being measured, and Gauss-Seidel over that coupling converges too slowly - /// to be an instrument (measured: it floors near 48 dB on exact synthetic - /// tones; the joint solve reaches the float quantization floor). - inline double joint_fit_residual_power(std::span x, std::span nus, - std::span fits) { - const std::size_t k = nus.size(); - const std::size_t n2 = 2 * k + 1; // +1: a DC column. Subtracting the - // sample mean beforehand is WRONG: a finite window of pure tones has a - // legitimate nonzero mean (partial cycles of the low tones), and - // removing it injects a constant the sine basis cannot absorb — a - // measured -48 dB instrument floor. Modeled jointly, DC costs nothing. - std::vector ata(n2 * n2, 0.0), aty(n2, 0.0), basis(n2), sol(n2); - for (std::size_t i = 0; i < x.size(); ++i) { - for (std::size_t t = 0; t < k; ++t) { - const double w = 2.0 * std::numbers::pi * nus[t] * static_cast(i); - basis[2 * t] = std::sin(w); - basis[2 * t + 1] = std::cos(w); - } - basis[n2 - 1] = 1.0; - for (std::size_t r = 0; r < n2; ++r) { - for (std::size_t q = r; q < n2; ++q) { - ata[r * n2 + q] += basis[r] * basis[q]; - } - aty[r] += basis[r] * x[i]; - } - } - for (std::size_t r = 0; r < n2; ++r) { - for (std::size_t q = 0; q < r; ++q) { - ata[r * n2 + q] = ata[q * n2 + r]; - } - } - tap::samplerate::detail::solve_dense(ata, aty, sol, n2); - for (std::size_t t = 0; t < k; ++t) { - fits[t] = tone_fit{sol[2 * t], sol[2 * t + 1]}; - } - double resid = 0.0; - for (std::size_t i = 0; i < x.size(); ++i) { - double model = sol[n2 - 1]; // DC: modeled out, in neither bucket - for (std::size_t t = 0; t < k; ++t) { - const double w = 2.0 * std::numbers::pi * nus[t] * static_cast(i); - model += sol[2 * t] * std::sin(w) + sol[2 * t + 1] * std::cos(w); - } - const double r = x[i] - model; - resid += r * r; - } - return resid / static_cast(x.size()); - } - - /// Program-weighted SNR: total fitted tone power over the power of what is - /// left after subtracting every tone — aliases, servo FM, noise, all of it. - /// - /// The comb is generated at physical Hz, so the converted tones sit at - /// freqHz/fsOut regardless of the clock offset. At these SNR levels a fit - /// frequency must be exact to ~1e-9 relative, so the servo's sub-ppm - /// settling residual is estimated from the data: joint-fit at the nominal - /// ratio, re-track each tone against (residual + that tone), pool the - /// implied ratios amplitude-weighted — every tone rides the SAME physical - /// clock ratio, so pooling averages the tracking noise down — then - /// joint-fit once more at the pooled ratio. - inline double program_weighted_snr_db(std::span tail, const tone_comb& comb, double /*fsIn*/, - double fs_out) { - std::vector work(tail.begin(), tail.end()); - const std::size_t k = comb.freq_hz.size(); - std::vector nus(k); - for (std::size_t t = 0; t < k; ++t) { - nus[t] = comb.freq_hz[t] / fs_out; - } - std::vector fits(k); - joint_fit_residual_power(work, nus, fits); - - // Ratio refinement against the joint residual, two rounds. Pooling - // weight is (amplitude * frequency)^2 — Fisher weighting: a tone's - // phase drift over the window scales with its frequency, so the high - // tones carry nearly all the ratio information even though pink - // weighting makes them the quietest. (Amplitude-only weighting leaves - // the ratio unresolved and floors the whole instrument at -48 dB.) - std::vector lone(work.size()); - double resid_power = 0.0; - for (int round = 0; round < 2; ++round) { - std::vector resid(work); - for (std::size_t i = 0; i < work.size(); ++i) { - double model = 0.0; - for (std::size_t t = 0; t < k; ++t) { - const double w = 2.0 * std::numbers::pi * nus[t] * static_cast(i); - model += fits[t].a * std::sin(w) + fits[t].b * std::cos(w); - } - resid[i] -= model; - } - double rho_num = 0.0, rho_den = 0.0; - for (std::size_t t = 0; t < k; ++t) { - const double w = 2.0 * std::numbers::pi * nus[t]; - for (std::size_t i = 0; i < lone.size(); ++i) { - lone[i] = resid[i] + fits[t].a * std::sin(w * static_cast(i)) - + fits[t].b * std::cos(w * static_cast(i)); - } - const double refined = track_tone_freq(lone, nus[t]); - const double wt = comb.amplitude[t] * nus[t]; - rho_num += wt * wt * (refined / nus[t]); - rho_den += wt * wt; - } - const double rho = rho_num / rho_den; - for (std::size_t t = 0; t < k; ++t) { - nus[t] *= rho; - } - resid_power = joint_fit_residual_power(work, nus, fits); - } - double signal = 0.0; - for (const auto& f : fits) { - signal += f.power(); - } - return 10.0 * std::log10(signal / resid_power); - } - // ANCHOR_END: pw_metric + using tap::dsp::analysis::fit_tone_fixed; + using tap::dsp::analysis::joint_fit_residual_power; + using tap::dsp::analysis::program_weighted_snr_db; + using tap::dsp::analysis::tone_comb; + using tap::dsp::analysis::tone_fit; + using tap::dsp::analysis::track_tone_freq; } // namespace srt_test diff --git a/tests/support/sine_analysis.h b/tests/support/sine_analysis.h index 691a3a8..7bbed87 100644 --- a/tests/support/sine_analysis.h +++ b/tests/support/sine_analysis.h @@ -1,109 +1,17 @@ // SPDX-License-Identifier: MIT // Copyright 2026 SampleRateTap contributors -// Least-squares sine fitting for THD+N-style quality measurements. +// Re-export of the shared tap::dsp::analysis sine-fit instrument. The +// implementation moved to DspTap with the measurement harness; this header +// keeps the historical include path and the srt_test names working. #pragma once -#include -#include -#include -#include +#include "tap/dsp/analysis/sine_analysis.h" namespace srt_test { - struct sine_fit { - double amplitude = 0.0; - double phase = 0.0; - double dc = 0.0; - double residual_rms = 0.0; - double freq_norm = 0.0; - }; - - /// Fits x[i] ~ a*sin(w i) + b*cos(w i) + c by least squares (3x3 normal - /// equations) at the known normalized frequency, then measures the residual. - inline sine_fit fit_sine(std::span x, double freq_norm) { - const double w = 2.0 * std::numbers::pi * freq_norm; - double m[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}; - double rhs[3] = {0, 0, 0}; - for (std::size_t i = 0; i < x.size(); ++i) { - const double s = std::sin(w * static_cast(i)); - const double c = std::cos(w * static_cast(i)); - const double basis[3] = {s, c, 1.0}; - for (int r = 0; r < 3; ++r) { - for (int q = 0; q < 3; ++q) { - m[r][q] += basis[r] * basis[q]; - } - rhs[r] += basis[r] * static_cast(x[i]); - } - } - // Gaussian elimination with partial pivoting. - int order[3] = {0, 1, 2}; - for (int col = 0; col < 3; ++col) { - int piv = col; - for (int r = col + 1; r < 3; ++r) { - if (std::abs(m[order[r]][col]) > std::abs(m[order[piv]][col])) { - piv = r; - } - } - std::swap(order[col], order[piv]); - const int p = order[col]; - for (int r = col + 1; r < 3; ++r) { - const int rr = order[r]; - const double f = m[rr][col] / m[p][col]; - for (int q = col; q < 3; ++q) { - m[rr][q] -= f * m[p][q]; - } - rhs[rr] -= f * rhs[p]; - } - } - double sol[3]; - for (int col = 2; col >= 0; --col) { - const int p = order[col]; - double v = rhs[p]; - for (int q = col + 1; q < 3; ++q) { - v -= m[p][q] * sol[q]; - } - sol[col] = v / m[p][col]; - } - sine_fit fit; - fit.amplitude = std::hypot(sol[0], sol[1]); - fit.phase = std::atan2(sol[1], sol[0]); - fit.dc = sol[2]; - double sq = 0.0; - for (std::size_t i = 0; i < x.size(); ++i) { - const double s = std::sin(w * static_cast(i)); - const double c = std::cos(w * static_cast(i)); - const double r = static_cast(x[i]) - (sol[0] * s + sol[1] * c + sol[2]); - sq += r * r; - } - fit.residual_rms = std::sqrt(sq / static_cast(x.size())); - fit.freq_norm = freq_norm; - return fit; - } - - /// Like fitSine, but refines the frequency first (a few iterations comparing - /// the fitted phase of the two window halves). An ASRC's rate estimate - /// converges asymptotically, so the tail of a run can sit a fraction of a ppm - /// off the nominal ratio; a rigid fixed-frequency fit would book that - /// (inaudible) offset as residual. Tracking the fundamental is standard - /// THD-analyzer practice. - inline sine_fit fit_sine_tracked(std::span x, double freq_norm_guess) { - double f = freq_norm_guess; - const std::size_t half = x.size() / 2; - for (int iter = 0; iter < 4; ++iter) { - const sine_fit a = fit_sine(x.first(half), f); - const sine_fit b = fit_sine(x.subspan(half), f); - // b.phase is relative to the second half's start; predict it from a. - const double two_pi = 2.0 * std::numbers::pi; - const double predicted = a.phase + two_pi * f * static_cast(half); - const double dphi = std::remainder(b.phase - predicted, two_pi); - f += dphi / (two_pi * static_cast(half)); - } - return fit_sine(x, f); - } - - /// Signal-to-(residual) ratio in dB for a fitted sine. - inline double snr_db(const sine_fit& f) { - return 10.0 * std::log10((f.amplitude * f.amplitude * 0.5) / (f.residual_rms * f.residual_rms)); - } + using tap::dsp::analysis::fit_sine; + using tap::dsp::analysis::fit_sine_tracked; + using tap::dsp::analysis::sine_fit; + using tap::dsp::analysis::snr_db; } // namespace srt_test