Add the shared FIR substrate carried from SampleRateTap

CLAUDE.md

@@ -8,19 +8,27 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 portable C++20 (standard library only, no frameworks), consumed as a git submodule by the
 individual libraries (TapTools pins it as `submodules/dsptap`; the AmbiTap/MuTap lineage is where
 the FFT came from). Four primitives today: the real FFT (`fft.h`), the YIN pitch detector
-(`yin.h`), and two pitch shifters (`psola.h`, `pvoc.h`). See `README.md` for each primitive's
-contract summary.
+(`yin.h`), and two pitch shifters (`psola.h`, `pvoc.h`) — plus the FIR substrate carried from
+SampleRateTap for the two rate converters (SampleRateTap, RatioTap): Kaiser prototype design
+(`kaiser.h`), the sample-format traits (`sample_traits.h`: float/Q15/Q31), the FIR dot kernels
+(`fir_kernels.h`), row-sum-preserving quantization (`quantize.h`), and the measurement
+instruments (`analysis/`). See `README.md` for each asset's contract summary.
 
 ## The design discipline (load-bearing — every primitive follows it)
 
 - **Fixed numeric contracts.** Each header documents its packing, conventions, normalization, and
   latency as *numbers*, and the tests pin them. Changing a documented contract point is a breaking
   change for every consumer.
-- **Double is the golden model; float32 is the embedded profile.** `basic_*<Sample>` templates
-  with `using x = basic_x<double>` / `x32 = basic_x<float>` aliases. The double path never changes
-  for speed; accelerated float backends (vDSP, CMSIS-Helium for the FFT) must re-present the
-  *exact* golden contract, so the double test battery stays a valid oracle. Cross-precision
-  agreement is pinned by tests.
+- **Double is the golden model; float32 is the embedded profile; Q15/Q31 are format-limited
+  embedded profiles.** `basic_*<Sample>` templates with `using x = basic_x<double>` /
+  `x32 = basic_x<float>` aliases. The double path never changes for speed; accelerated float
+  backends (vDSP, CMSIS-Helium for the FFT) must re-present the *exact* golden contract, so the
+  double test battery stays a valid oracle. Cross-precision agreement is pinned by tests. The
+  fixed-point profiles (`sample_traits.h`) carry their contracts as numbers the same way — Q
+  formats, the single rounding point, saturation — and exist for M33/M55-class targets
+  (Bluetooth-adjacent converters, eurorack/pedal deployments) where double or any float is
+  unaffordable. Per-primitive fixed-point adoption is opt-in and is a documented Q-format design
+  each time, via traits over raw sample types, never wrapper classes.
 - **Real-time safe by construction.** Geometry fixed at construction, every buffer allocated
   there; processing is `noexcept` and allocation-free. Numerically fragile recursions (e.g. the
   order-48 Levinson–Durbin inside `pvoc`) run in double even in the float profile — documented

README.md

@@ -4,7 +4,9 @@ Shared DSP primitives for the **Tap** family of audio libraries. Header-only,
 plain portable C++ (C++20, standard library only), no Max/Min or framework
 dependency — consumed as a git submodule by the individual libraries.
 
-Today it holds four primitives:
+Today it holds four primitives, plus the [FIR substrate](#the-fir-substrate) —
+the shared design-math / sample-format / kernel layer under SampleRateTap and
+RatioTap:
 
 ## `tap::dsp::real_fft` — real FFT with a fixed numeric contract
 
@@ -125,6 +127,85 @@ per analysis frame, with every relocated bin rescaled by
 At ratio 1 the correction is exactly unity, so the identity contract holds
 either way.
 
+## The FIR substrate
+
+Five headers carried from **SampleRateTap** (where they design and run the
+ASRC's polyphase datapath) and promoted here so **RatioTap**'s fixed-ratio
+44.1↔48 converter — and any future FIR consumer — shares one implementation.
+The performance-sensitive pieces are regression-gated in SampleRateTap's
+instruction-count CI (Cortex-M33/M55, Hexagon, ±3%); treat measured claims in
+the header comments as contracts.
+
+### `tap/dsp/kaiser.h` — FIR prototype design
+
+Kaiser-windowed sinc prototype design for L-phase polyphase banks: `bessel_i0`,
+`kaiser_beta` (Kaiser's empirical fit), `estimate_taps` (the harris length
+estimate), `design_prototype`, and `design_prototype_compensated` — the
+zeros-at-k·fs variant with passband droop pre-compensated (closed-form, no
+FFT), which turns branch-DC uniformity into exact transmission zeros at every
+multiple of the sample rate. Runtime design in double, deliberately not
+constexpr (the header's design note does the arithmetic); run it in a
+constructor, off the audio path. Also exports `solve_dense`, the small dense
+solver the compensated design and the analysis instruments share.
+
+### `tap/dsp/sample_traits.h` — sample formats: float, Q15, Q31
+
+The family's sample-format substrate: how each sample type stores
+coefficients, accumulates dot products, and rounds/saturates back to samples.
+
+| Type | Coefficients | Accumulation | Output |
+|---|---|---|---|
+| `float` | float | double | plain cast |
+| `std::int16_t` | Q1.14 | int64, exact | single Q29→Q15 round-half-up, saturating |
+| `std::int32_t` | Q1.30 | int64, products pre-shifted to Q45 | Q45→Q31, saturating |
+
+**Fixed point is a first-class embedded direction, not a legacy path.** The
+Q15/Q31 profiles exist for targets where double (sometimes any float) is
+unaffordable — SampleRateTap measured its float datapath at ~19× the
+instruction count of Q15 on a Cortex-M33 (soft-double accumulation). Expected
+deployments include Bluetooth-adjacent conversion (RatioTap) and M33/M55-class
+eurorack and pedal targets running TapTools primitives. Per-primitive adoption
+is opt-in, and each adoption is its own documented Q-format design: the ladder
+of headroom bits, pre-shifts, and the single rounding point is a per-datapath
+decision. That is also why these are *traits over raw sample types* rather
+than `q15`/`q31` wrapper classes — the arithmetic contract stays visible at
+the use site and pinnable by tests, buffers arrive from codecs and C ABIs as
+plain `int16_t`/`int32_t`, and the SMLALD kernel's paired loads stay legal.
+
+This header is the format core only. Engine-specific extensions (e.g.
+SampleRateTap's inter-phase coefficient blending) derive from these
+specializations and refine the `tap::dsp::sample_type` concept.
+
+### `tap/dsp/fir_kernels.h` — dot-product kernels
+
+The FIR hot loops, target-gated the way SampleRateTap's optimization campaign
+measured them: `dot_row` (planar; routes Q15 through a dual-MAC SMLALD loop on
+DSP-extension Arm cores without Helium — bit-exact by construction), and the
+channel-parallel pair `dot_tile_frame_major` / `dot_rows_frame_major`
+(register-blocked 8/4/2/1 tiles over frame-major storage, coefficient
+broadcast across channel lanes — bit-exact against the planar path for every
+sample type, float included, because lanes are channels, not taps). The
+`TAP_DSP_CHANNEL_PARALLEL` / `TAP_DSP_CP_MIN_CHANNELS` gates encode which
+targets prefer which layout.
+
+### `tap/dsp/quantize.h` — row-sum-preserving quantization
+
+`quantize_row_preserving_sum`: quantizes one polyphase branch to a fixed-point
+coefficient format while preserving the row's DC sum *exactly*
+(largest-remainder distribution of the rounding residual — "the coefficients
+of every phase must add to one", R. Bristow-Johnson, music-dsp). Plain
+conversion for float. Design-time code.
+
+### `tap/dsp/analysis/` — measurement instruments
+
+The quality-measurement harness the converter suites share: `sine_analysis.h`
+(least-squares single-tone fit, frequency-tracked variant, `snr_db`) and
+`multitone_analysis.h` (pink log-spaced `tone_comb`, joint least-squares
+multitone fit, `program_weighted_snr_db` — the program-weighted metric with
+Fisher-weighted ratio pooling). Instrument floors on exact synthetic signals
+are pinned by `tests/test_analysis.cpp`, so a consumer's quality gate never
+silently rests on a degraded instrument.
+
 ## Notebooks
 
 `notebooks/pitchshift.ipynb` measures the three pitch primitives — driving the
@@ -175,6 +256,13 @@ would silently miss the other. DspTap is the consolidation: one wrapper, one
 contract, one home for the next backend. The unified wrapper is MuTap's
 backend-capable `basic_real_fft`, generalized to the `tap::dsp` namespace.
 
+The FIR substrate is the second consolidation wave, moved here from
+**SampleRateTap** (its `srt/detail/kaiser.h`, `srt/sample_traits.h` format
+core, the dot kernels from `srt/polyphase_filter.h`, the row-sum quantization
+from its bank constructor, and the `tests/support/` measurement harness) at
+the moment **RatioTap** became the second consumer — the same
+extract-on-second-consumer rule that created this repo.
+
 See [`third_party/ooura/readme.txt`](third_party/ooura/readme.txt) and
 [`third_party/cmsis-dsp/VENDOR.md`](third_party/cmsis-dsp/VENDOR.md) for the
 vendored-code provenance and licenses.

include/tap/dsp/analysis/multitone_analysis.h

@@ -0,0 +1,229 @@
+/// @file multitone_analysis.h
+/// @brief Program-weighted multitone metric: pink tone comb + joint-fit residual.
+// SPDX-License-Identifier: MIT
+// Copyright 2026 Timothy Place and the DspTap contributors.
+//
+// Carried from SampleRateTap's test harness (tests/support/), promoted here
+// as a shared measurement instrument.
+//
+// Why this exists: single-sine SNR is the worst-case metric. A filter with
+// transmission zeros at k*fs (the compensated designs in kaiser.h) 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. The design
+// history is SampleRateTap's epilogue chapter.
+#pragma once
+
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <numbers>
+#include <span>
+#include <vector>
+
+#include "tap/dsp/analysis/sine_analysis.h"
+#include "tap/dsp/kaiser.h" // solve_dense for the joint LS
+
+namespace tap::dsp::analysis {
+
+    // ANCHOR: pw_comb
+    struct tone_comb {
+        std::vector<double> freq_hz;   // log-spaced
+        std::vector<double> amplitude; // pink: a ~ 1/sqrt(f), scaled to peakSum
+        std::vector<double> 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<double>(i) / static_cast<double>(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<double>(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<double>(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<const double> 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<double>(i));
+            const double c = std::cos(w * static_cast<double>(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<const double> 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<double>(half);
+            const double   dphi      = std::remainder(std::atan2(b.b, b.a) - predicted, two_pi);
+            f += dphi / (two_pi * static_cast<double>(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<const double> x, std::span<const double> nus,
+                                           std::span<tone_fit> 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<double> 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<double>(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];
+            }
+        }
+        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<double>(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<double>(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 any clock offset. At these SNR levels a fit
+    /// frequency must be exact to ~1e-9 relative, so a converter'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<const float> tail, const tone_comb& comb, double /*fsIn*/,
+                                          double fs_out) {
+        std::vector<double> work(tail.begin(), tail.end());
+        const std::size_t   k = comb.freq_hz.size();
+        std::vector<double> nus(k);
+        for (std::size_t t = 0; t < k; ++t) {
+            nus[t] = comb.freq_hz[t] / fs_out;
+        }
+        std::vector<tone_fit> 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<double> lone(work.size());
+        double              resid_power = 0.0;
+        for (int round = 0; round < 2; ++round) {
+            std::vector<double> 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<double>(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<double>(i))
+                              + fits[t].b * std::cos(w * static_cast<double>(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
+
+} // namespace tap::dsp::analysis