diff --git a/CLAUDE.md b/CLAUDE.md index 20a7d86..344176f 100644 --- a/CLAUDE.md +++ b/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_*` templates - with `using x = basic_x` / `x32 = basic_x` 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_*` templates with `using x = basic_x` / + `x32 = basic_x` 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 diff --git a/README.md b/README.md index 076dd82..a7c992f 100644 --- a/README.md +++ b/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. diff --git a/a.out b/a.out deleted file mode 100755 index 67cd4ba..0000000 Binary files a/a.out and /dev/null differ diff --git a/include/tap/dsp/analysis/multitone_analysis.h b/include/tap/dsp/analysis/multitone_analysis.h new file mode 100644 index 0000000..b60c037 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include + +#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 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]; + } + } + 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 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 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 + +} // namespace tap::dsp::analysis diff --git a/include/tap/dsp/analysis/sine_analysis.h b/include/tap/dsp/analysis/sine_analysis.h new file mode 100644 index 0000000..63d8de0 --- /dev/null +++ b/include/tap/dsp/analysis/sine_analysis.h @@ -0,0 +1,116 @@ +/// @file sine_analysis.h +/// @brief Least-squares sine fitting for THD+N-style quality measurements. +// 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: fit a sine of known (or tracked) +// frequency by least squares, subtract it, and score the residual. Used by +// converter quality suites (SampleRateTap, RatioTap) and available to any +// primitive that needs a single-tone SNR number. +#pragma once + +#include +#include +#include +#include + +namespace tap::dsp::analysis { + + 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). A converter'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)); + } + +} // namespace tap::dsp::analysis diff --git a/include/tap/dsp/fir_kernels.h b/include/tap/dsp/fir_kernels.h new file mode 100644 index 0000000..70069b0 --- /dev/null +++ b/include/tap/dsp/fir_kernels.h @@ -0,0 +1,168 @@ +/// @file fir_kernels.h +/// @brief FIR dot-product kernels: planar, SMLALD dual-MAC, and channel-parallel. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Carried from SampleRateTap (include/srt/polyphase_filter.h), where these +// kernels are the ASRC's hot loop; promoted here so RatioTap's fixed-ratio +// converter runs the identical, already-measured code. The performance claims +// cited below were measured in SampleRateTap's optimization campaign and are +// regression-gated there by instruction-count CI on Cortex-M33/M55 and +// Hexagon (see that repo's docs/PERFORMANCE.md); a change here lands in every +// consumer on its next submodule bump, so treat the measured comments as +// contracts. +#pragma once + +#include +#include +#include +#include + +#include "tap/dsp/sample_traits.h" + +// No-alias qualifier for the kernel hot loops: without it the compiler +// versions loops over distinct row/history spans behind a runtime aliasing +// check (verified with -fopt-info-vec; SampleRateTap docs/PERFORMANCE.md, +// hypothesis 2). +#if defined(_MSC_VER) +#define TAP_DSP_RESTRICT __restrict +#else +#define TAP_DSP_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 (SampleRateTap +// 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 TAP_DSP_Q15_SMLALD 1 +#else +#define TAP_DSP_Q15_SMLALD 0 +#endif +// ANCHOR_END: opt_smlald_gate + +// Channel-parallel dot product for high channel counts (SampleRateTap +// docs/PERFORMANCE.md, hypothesis C6): 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 +// (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 — hypotheses C4/C5). +#if !defined(__ARM_FEATURE_MVE) && !defined(__ARM_FEATURE_DSP) && !defined(__hexagon__) +#define TAP_DSP_CHANNEL_PARALLEL 1 +#else +#define TAP_DSP_CHANNEL_PARALLEL 0 +#endif +// Minimum channel count for the frame-major path (overridable for A/B +// measurements; a blend-share planar path stays better at low counts). +#ifndef TAP_DSP_CP_MIN_CHANNELS +#define TAP_DSP_CP_MIN_CHANNELS 4 +#endif + +namespace tap::dsp { + + // ANCHOR: rs_dot_row + /// Dot product of a coefficient row against a history window, in the + /// sample type's accumulator domain (see sample_traits.h): tap-order + /// accumulation, single rounding in finalize. + template + inline S dot_row(const typename sample_traits::coeff* TAP_DSP_RESTRICT row, const S* TAP_DSP_RESTRICT hist, + std::size_t taps) noexcept { + using tr = sample_traits; +#if TAP_DSP_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 + 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* TAP_DSP_RESTRICT row, + const S* TAP_DSP_RESTRICT x, std::size_t taps, std::size_t stride, + S* TAP_DSP_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* TAP_DSP_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 dot_row'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* TAP_DSP_RESTRICT row, + const S* TAP_DSP_RESTRICT x, std::size_t taps, std::size_t channels, + S* TAP_DSP_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 + +} // namespace tap::dsp diff --git a/include/tap/dsp/kaiser.h b/include/tap/dsp/kaiser.h new file mode 100644 index 0000000..99a3fd7 --- /dev/null +++ b/include/tap/dsp/kaiser.h @@ -0,0 +1,375 @@ +/// @file kaiser.h +/// @brief Kaiser-window FIR prototype design for polyphase filter banks. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Carried from SampleRateTap (include/srt/detail/kaiser.h), where it designs +// the ASRC's interpolation bank; promoted here so RatioTap's fixed-ratio +// tables and any future FIR design share one implementation. All methods are +// published literature: Kaiser window design (Kaiser 1974), band-limited +// interpolation (J. O. Smith, CCRMA), the harris length estimate, and the +// zeros-at-k*fs droop-compensated variant worked out in the June 2026 +// music-dsp thread (verified in SampleRateTap's asrc_rbj_analysis notebook). +// +// ANCHOR: kai_design_note +/// Design note — runtime vs constexpr: 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 +#pragma once + +#include +#include +#include +#include +#include + +namespace tap::dsp { + + // 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, and the analysis instruments' + /// joint tone fits stay under ~50. + 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 SampleRateTap's epilogue + /// chapter and its asrc_rbj_analysis notebook). + /// + /// 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 SampleRateTap's + /// 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 a consuming 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 — SampleRateTap's 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(); + } + +} // namespace tap::dsp diff --git a/include/tap/dsp/quantize.h b/include/tap/dsp/quantize.h new file mode 100644 index 0000000..ff46f28 --- /dev/null +++ b/include/tap/dsp/quantize.h @@ -0,0 +1,83 @@ +/// @file quantize.h +/// @brief Row-sum-preserving coefficient quantization for polyphase tables. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Carried from SampleRateTap, where this correction ran inline in the +// polyphase bank constructor; promoted here because every polyphase table +// with unity per-phase DC gain wants it — RatioTap's fixed-ratio tables +// exactly as much as the ASRC's interpolation bank. +#pragma once + +#include +#include +#include +#include +#include + +#include "tap/dsp/sample_traits.h" + +namespace tap::dsp { + + // ANCHOR: pw_row_sum + /// Quantizes one polyphase branch's coefficients (double, in units where + /// 1.0 is unity gain) to the sample type's coefficient format, preserving + /// the row's DC sum exactly. + /// + /// "Note that for DC, this should get you infinite S/N ratio... for every + /// phase or fractional delay, the FIR coefficients must add to 1." + /// -- R. Bristow-Johnson, music-dsp. In double, a well-designed bank + /// makes every branch's DC sum identical to machine epsilon (zeros at + /// k*fs ARE branch-DC uniformity, stated in frequency), but independent + /// per-tap rounding re-breaks it by several LSB. This distributes 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_sum * k_coeff_scale), so a table built row-by-row holds + /// DC gain within one coefficient LSB across all phases. + /// + /// For floating-point coefficient types this is a plain conversion (no + /// correction runs; k_coeff_scale is 1). Design-time code: allocates a + /// scratch vector for integer formats, so keep it off the audio path. + /// + /// \pre dst.size() == src.size() + template + inline void quantize_row_preserving_sum(std::span src, + std::span::coeff> dst) { + using tr = sample_traits; + using coeff = typename tr::coeff; + const std::size_t n = src.size(); + if constexpr (std::is_floating_point_v) { + for (std::size_t t = 0; t < n; ++t) { + dst[t] = tr::make_coeff(src[t]); + } + } + else { + std::vector remainder(n); + double exact_sum = 0.0; + std::int64_t quant_sum = 0; + for (std::size_t t = 0; t < n; ++t) { + const double scaled = src[t] * tr::k_coeff_scale; + const coeff q = tr::make_coeff(src[t]); + dst[t] = q; + remainder[t] = scaled - static_cast(q); + exact_sum += scaled; + quant_sum += static_cast(q); + } + std::int64_t residual = static_cast(std::llround(exact_sum)) - quant_sum; + while (residual != 0) { + const double sgn = residual > 0 ? 1.0 : -1.0; + std::size_t best = 0; + for (std::size_t u = 1; u < n; ++u) { + if (sgn * remainder[u] > sgn * remainder[best]) { + best = u; + } + } + dst[best] = static_cast(dst[best] + (residual > 0 ? 1 : -1)); + remainder[best] -= sgn; + residual -= residual > 0 ? 1 : -1; + } + } + } + // ANCHOR_END: pw_row_sum + +} // namespace tap::dsp diff --git a/include/tap/dsp/sample_traits.h b/include/tap/dsp/sample_traits.h new file mode 100644 index 0000000..7ad396e --- /dev/null +++ b/include/tap/dsp/sample_traits.h @@ -0,0 +1,202 @@ +/// @file sample_traits.h +/// @brief Sample-format customization point for FIR datapaths: float, Q15, Q31. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Carried from SampleRateTap (include/srt/sample_traits.h), where it +// parameterizes the ASRC datapath; promoted here as the family's sample-format +// substrate, shared with RatioTap's fixed-ratio hot loop and available to any +// primitive growing an embedded fixed-point profile. This header is the +// *format core* only — how each sample type stores coefficients, accumulates +// dot products, and rounds/saturates back to samples. Engine-specific +// extensions (e.g. SampleRateTap's inter-phase coefficient blending) layer on +// top by deriving from these specializations and refining the concept. +// +// Three sample types are provided: +// +// - float : float I/O and coefficients, double accumulation +// - std::int16_t : Q15 samples, Q1.14 coefficients, int64 accumulation, +// saturating output +// - std::int32_t : Q31 samples, Q1.30 coefficients, int64 accumulation, +// saturating output +// +// The fixed-point profiles exist for targets where double (and sometimes any +// float) is unaffordable: M33-class cores measured ~19x the instruction count +// running SampleRateTap's float datapath vs Q15 (soft-double accumulation). +// Expected deployments include Bluetooth-adjacent converters (RatioTap) and +// M33/M55-class eurorack/pedal targets (TapTools primitives, opt-in per +// primitive — each adoption is its own documented Q-format design). +// +// Deliberately raw sample types + traits, not wrapper classes with operator +// overloads: the Q-format ladder (headroom bits, accumulator pre-shifts, +// where the single rounding happens) is a per-datapath design decision that +// must stay visible at the use site and pinnable by tests — see the Q31 +// mac() below for an example a `q31` class could not express honestly. Raw +// int16_t/int32_t also match what codecs and C ABIs actually deliver, and +// keep memcpy/SIMD access patterns (e.g. the SMLALD pair loads in +// fir_kernels.h) legal. +#pragma once + +#include +#include +#include +#include + +namespace tap::dsp { + + 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)); + } + + } // namespace detail + + /// Primary template intentionally undefined; specialize per sample type. + template + struct sample_traits; + + /// Float datapath: float samples and coefficients, double accumulation. + /// The double accumulator keeps the dot-product noise floor far below a + /// 120 dB transparency target; float coefficient storage quantizes the + /// filter at roughly -150 dB, negligible against the same target. + template <> + struct sample_traits { + using coeff = float; ///< stored filter coefficient type + using accum = double; ///< dot-product accumulator type + + /// Convert a double-precision designed coefficient to storage form. + static coeff make_coeff(double c) noexcept { return static_cast(c); } + /// Coefficient units per 1.0 (used by row-sum-preserving quantization; + /// unity for floating storage, where no correction runs). + static constexpr double k_coeff_scale = 1.0; + + /// 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); + } + + /// 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: st_q15_core + /// Q15 fixed-point datapath (samples are int16_t in Q0.15). + /// + /// Coefficients are stored in Q1.14: a unity-DC prototype's peak tap + /// reaches ~1.0, 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 a converter built on this is + /// Q15-transparent. + template <> + struct sample_traits { + using coeff = std::int16_t; + using accum = std::int64_t; + + // 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 + + // 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_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_END: st_q15_core + + // ANCHOR: st_q31_core + /// 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. + template <> + struct sample_traits { + using coeff = std::int32_t; + using accum = std::int64_t; + + static coeff make_coeff(double c) noexcept { + return detail::round_sat(c * 1073741824.0); // Q1.30 + } + static constexpr double k_coeff_scale = 1073741824.0; // Q1.30 units per 1.0 + + // ANCHOR: st_q31_mac + static accum mac(accum acc, std::int32_t x, coeff c) noexcept { + return acc + ((static_cast(x) * c) >> 16); // Q61 -> Q45 + } + // ANCHOR_END: st_q31_mac + + 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_core + + // ANCHOR: st_core_concept + /// Satisfied by any type with a complete format-core sample_traits + /// specialization — everything the FIR kernels (fir_kernels.h) and table + /// builders require. Consumers with richer datapaths refine this concept + /// over their own trait extensions (SampleRateTap's sample_type adds the + /// inter-phase blend contract on top). + template + concept sample_type = + requires(T x, double d, typename sample_traits::accum a, typename sample_traits::coeff c) { + { sample_traits::make_coeff(d) } -> std::same_as::coeff>; + { sample_traits::mac(a, x, c) } -> std::same_as::accum>; + { sample_traits::finalize(a) } -> std::same_as; + { sample_traits::silence() } -> std::same_as; + }; + + static_assert(sample_type); + static_assert(sample_type); + static_assert(sample_type); + // ANCHOR_END: st_core_concept + +} // namespace tap::dsp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 75f134c..e3abd01 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -20,10 +20,15 @@ set(INSTALL_GTEST OFF CACHE BOOL "" FORCE) FetchContent_MakeAvailable(googletest) add_executable(tap_dsp_tests + test_analysis.cpp test_fft.cpp test_fft_backend.cpp + test_fir_kernels.cpp + test_kaiser.cpp test_psola.cpp test_pvoc.cpp + test_quantize.cpp + test_sample_traits.cpp test_yin.cpp) target_link_libraries(tap_dsp_tests PRIVATE tap::dsp diff --git a/tests/test_analysis.cpp b/tests/test_analysis.cpp new file mode 100644 index 0000000..235e1e8 --- /dev/null +++ b/tests/test_analysis.cpp @@ -0,0 +1,79 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Smoke battery for the measurement instruments. These headers earn their +// keep in consumer suites (SampleRateTap's quality tests, RatioTap's alias +// acceptance); here we pin that the instruments themselves reach their +// documented floors on exact synthetic signals — an instrument that cannot +// measure a clean signal cleanly would silently weaken every consumer's +// quality gate. + +#include +#include +#include + +#include + +#include "tap/dsp/analysis/multitone_analysis.h" +#include "tap/dsp/analysis/sine_analysis.h" + +namespace { + + namespace an = tap::dsp::analysis; + + TEST(SineAnalysis, FitRecoversExactTone) { + const double nu = 997.0 / 48000.0; + std::vector x(16384); + for (std::size_t i = 0; i < x.size(); ++i) { + x[i] = static_cast(0.5 * std::sin(2.0 * std::numbers::pi * nu * static_cast(i) + 0.3)); + } + const an::sine_fit fit = an::fit_sine(x, nu); + EXPECT_NEAR(fit.amplitude, 0.5, 1e-6); + EXPECT_NEAR(fit.dc, 0.0, 1e-6); + // Residual is the float storage quantization only: > 100 dB SNR. + EXPECT_GT(an::snr_db(fit), 100.0); + } + + TEST(SineAnalysis, TrackedFitAbsorbsSmallFrequencyOffset) { + // The tone sits 5 ppm off the guess — the situation the tracker exists + // for (a converter's rate estimate settling asymptotically). A rigid + // fit at the guess books the offset as residual; the tracked fit must + // recover the true frequency and the quantization-limited floor. + const double nu_true = (997.0 / 48000.0) * (1.0 + 5e-6); + std::vector x(65536); + for (std::size_t i = 0; i < x.size(); ++i) { + x[i] = static_cast(0.5 * std::sin(2.0 * std::numbers::pi * nu_true * static_cast(i))); + } + const an::sine_fit fit = an::fit_sine_tracked(x, 997.0 / 48000.0); + EXPECT_NEAR(fit.freq_norm, nu_true, 1e-10); + EXPECT_NEAR(fit.amplitude, 0.5, 1e-6); + EXPECT_GT(an::snr_db(fit), 100.0); + } + + TEST(MultitoneAnalysis, CombIsBoundedByPeakSum) { + const an::tone_comb comb = an::tone_comb::pink(24, 40.0, 20000.0, 0.5); + ASSERT_EQ(comb.freq_hz.size(), 24u); + double amp_sum = 0.0; + for (const double a : comb.amplitude) { + amp_sum += a; + } + EXPECT_NEAR(amp_sum, 0.5, 1e-12); + for (std::uint64_t i = 0; i < 48000; ++i) { + ASSERT_LE(std::abs(comb.sample_at(i, 48000.0)), 0.5 + 1e-9); + } + } + + TEST(MultitoneAnalysis, JointFitReachesQuantizationFloorOnExactTones) { + // The documented instrument floor claim: sequential fit-subtract + // floors near 48 dB on exact synthetic tones; the joint solve reaches + // the float storage quantization floor. Pin the latter. + const an::tone_comb comb = an::tone_comb::pink(12, 40.0, 20000.0, 0.5); + std::vector tail(32768); + for (std::size_t i = 0; i < tail.size(); ++i) { + tail[i] = static_cast(comb.sample_at(i, 48000.0)); + } + const double snr = an::program_weighted_snr_db(tail, comb, 48000.0, 48000.0); + EXPECT_GT(snr, 90.0); + } + +} // namespace diff --git a/tests/test_fir_kernels.cpp b/tests/test_fir_kernels.cpp new file mode 100644 index 0000000..cfc5cc5 --- /dev/null +++ b/tests/test_fir_kernels.cpp @@ -0,0 +1,130 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Contract battery for the FIR dot kernels. The load-bearing promise — +// extracted from SampleRateTap's multichannel suite, where it gates the +// frame-major fast path — is bit-exactness: the channel-parallel kernel must +// produce the identical bits as the planar dot for every sample type, because +// consumers switch between the layouts by channel count and target. + +#include +#include + +#include + +#include "tap/dsp/fir_kernels.h" + +namespace { + + using tap::dsp::dot_row; + using tap::dsp::dot_rows_frame_major; + using tap::dsp::sample_traits; + + // Deterministic pseudo-random generator (xorshift), mapped per sample type + // to well-inside-full-scale values so no test depends on saturation. + inline std::uint32_t next(std::uint32_t& s) { + s ^= s << 13; + s ^= s >> 17; + s ^= s << 5; + return s; + } + + template + S sample_from(std::uint32_t r); + template <> + float sample_from(std::uint32_t r) { + return (static_cast(r % 65536) - 32768.0f) / 65536.0f; + } + template <> + std::int16_t sample_from(std::uint32_t r) { + return static_cast(static_cast(r % 32768) - 16384); + } + template <> + std::int32_t sample_from(std::uint32_t r) { + return static_cast(r) / 4; + } + + template + typename sample_traits::coeff coeff_from(std::uint32_t r) { + // Coefficients in roughly [-0.5, 0.5) of the format's unity. + return sample_traits::make_coeff((static_cast(r % 4096) - 2048.0) / 4096.0); + } + + template + class fir_kernels_test : public ::testing::Test {}; + using sample_types = ::testing::Types; + TYPED_TEST_SUITE(fir_kernels_test, sample_types, ); + + // dot_row must equal the reference accumulation: mac per tap in order, + // one finalize. (On SMLALD targets this also pins the dual-MAC pairing; + // pairing is bit-free because each 16x16 product is exact in int32.) + TYPED_TEST(fir_kernels_test, DotRowMatchesReferenceAccumulation) { + using sample = TypeParam; + using tr = sample_traits; + constexpr std::size_t k_taps = 48; + std::uint32_t seed = 0x12345678u; + std::vector hist(k_taps); + std::vector row(k_taps); + for (std::size_t t = 0; t < k_taps; ++t) { + hist[t] = sample_from(next(seed)); + row[t] = coeff_from(next(seed)); + } + typename tr::accum acc{}; + for (std::size_t t = 0; t < k_taps; ++t) { + acc = tr::mac(acc, hist[t], row[t]); + } + EXPECT_EQ(dot_row(row.data(), hist.data(), k_taps), tr::finalize(acc)); + } + + // The frame-major channel-parallel kernel is bit-exact against the planar + // dot for every channel count that exercises the 8/4/2/1 tiling, float + // included (lanes are channels, not taps: each channel's double + // accumulation order is unchanged). + TYPED_TEST(fir_kernels_test, ChannelParallelMatchesPlanarBitExact) { + using sample = TypeParam; + using tr = sample_traits; + constexpr std::size_t k_taps = 48; + for (std::size_t channels : {1u, 2u, 3u, 4u, 7u, 8u, 11u, 12u, 16u}) { + std::uint32_t seed = 0x9e3779b9u + static_cast(channels); + std::vector x(k_taps * channels); // frame-major + std::vector row(k_taps); + for (auto& v : x) { + v = sample_from(next(seed)); + } + for (auto& c : row) { + c = coeff_from(next(seed)); + } + std::vector out(channels); + dot_rows_frame_major(row.data(), x.data(), k_taps, channels, out.data()); + for (std::size_t c = 0; c < channels; ++c) { + std::vector planar(k_taps); + for (std::size_t t = 0; t < k_taps; ++t) { + planar[t] = x[t * channels + c]; + } + EXPECT_EQ(out[c], dot_row(row.data(), planar.data(), k_taps)) + << "channels=" << channels << " c=" << c; + } + } + } + + // Odd tap counts exercise dot_row's scalar tail on SMLALD targets; on + // hosts this just pins the same contract at a second geometry. + TYPED_TEST(fir_kernels_test, OddTapCountMatchesReference) { + using sample = TypeParam; + using tr = sample_traits; + constexpr std::size_t k_taps = 33; + std::uint32_t seed = 0xdeadbeefu; + std::vector hist(k_taps); + std::vector row(k_taps); + for (std::size_t t = 0; t < k_taps; ++t) { + hist[t] = sample_from(next(seed)); + row[t] = coeff_from(next(seed)); + } + typename tr::accum acc{}; + for (std::size_t t = 0; t < k_taps; ++t) { + acc = tr::mac(acc, hist[t], row[t]); + } + EXPECT_EQ(dot_row(row.data(), hist.data(), k_taps), tr::finalize(acc)); + } + +} // namespace diff --git a/tests/test_kaiser.cpp b/tests/test_kaiser.cpp new file mode 100644 index 0000000..ee866dd --- /dev/null +++ b/tests/test_kaiser.cpp @@ -0,0 +1,196 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Contract battery for the Kaiser prototype designer, ported from +// SampleRateTap's test_kaiser.cpp. The spec structs below carry that +// library's shipping preset numbers (fast/balanced/transparent/economy) so +// the designs it depends on stay pinned at their source. + +#include +#include +#include +#include +#include + +#include + +#include "tap/dsp/kaiser.h" + +namespace { + + using namespace tap::dsp; + + TEST(Kaiser, BesselI0ReferenceValues) { + EXPECT_DOUBLE_EQ(bessel_i0(0.0), 1.0); + EXPECT_NEAR(bessel_i0(1.0), 1.2660658777520084, 1e-12); + EXPECT_NEAR(bessel_i0(5.0), 27.239871823604442, 1e-9); + EXPECT_NEAR(bessel_i0(12.0), 18948.925349296309, 1e-6); + } + + TEST(Kaiser, BetaReferenceValues) { + EXPECT_NEAR(kaiser_beta(120.0), 0.1102 * (120.0 - 8.7), 1e-12); + EXPECT_NEAR(kaiser_beta(40.0), 0.5842 * std::pow(19.0, 0.4) + 0.07886 * 19.0, 1e-12); + EXPECT_DOUBLE_EQ(kaiser_beta(15.0), 0.0); + } + + TEST(Kaiser, TapEstimateMatchesHarrisFormula) { + // 120 dB over a 20->28 kHz transition at 48 kHz: ~47 taps per phase. + const std::size_t taps = estimate_taps(120.0, 8000.0 / 48000.0); + EXPECT_GE(taps, 45u); + EXPECT_LE(taps, 49u); + } + + // A prototype specification in SampleRateTap's filter_spec vocabulary. + struct proto_spec { + std::size_t num_phases; + std::size_t taps_per_phase; + double passband_hz; + double stopband_hz; + double stopband_atten_db; + bool image_zeros; + + proto_spec scaled_to(double sample_rate_hz) const { + constexpr double k_design_rate_hz = 48000.0; + proto_spec s = *this; + s.passband_hz *= sample_rate_hz / k_design_rate_hz; + s.stopband_hz *= sample_rate_hz / k_design_rate_hz; + return s; + } + }; + + // SampleRateTap's shipping presets (filter_spec::fast/balanced/transparent/ + // economy), pinned here at the substrate so a design change is caught where + // the code now lives. + constexpr proto_spec k_fast{128, 32, 18000.0, 30000.0, 96.0, false}; + constexpr proto_spec k_balanced{256, 48, 20000.0, 28000.0, 120.0, true}; + constexpr proto_spec k_transparent{512, 80, 20000.0, 26000.0, 140.0, true}; + constexpr proto_spec k_economy{512, 32, 18000.0, 30000.0, 96.0, true}; + + // Direct DFT magnitude of the double-precision prototype, normalized so the + // passband sits at 0 dB. f is in Hz; the prototype rate is L * fs. + double response_db(const std::vector& h, std::size_t num_phases, double fs, double f) { + const double proto_rate = static_cast(num_phases) * fs; + std::complex acc{0.0, 0.0}; + for (std::size_t m = 0; m < h.size(); ++m) { + const double ang = -2.0 * std::numbers::pi * f * static_cast(m) / proto_rate; + acc += h[m] * std::polar(1.0, ang); + } + return 20.0 * std::log10(std::abs(acc) / static_cast(num_phases)); + } + + void check_prototype_meets_spec(const proto_spec& spec, double fs) { + const std::size_t phases = std::bit_ceil(spec.num_phases); + const std::size_t n = phases * spec.taps_per_phase; + std::vector h(n); + const double cutoff_norm = (spec.passband_hz + spec.stopband_hz) / fs; + if (spec.image_zeros) { + design_prototype_compensated(h, phases, cutoff_norm, kaiser_beta(spec.stopband_atten_db), + spec.passband_hz / fs); + } + else { + design_prototype(h, phases, cutoff_norm, kaiser_beta(spec.stopband_atten_db)); + } + + // Passband: flat within +/-0.01 dB up to the passband edge. For the + // compensated designs this is the claim the droop pre-compensation + // exists to defend (the raw rect would sag -2.64 dB at 20 kHz). + for (double f = 0.0; f <= spec.passband_hz; f += 500.0) { + EXPECT_NEAR(response_db(h, spec.num_phases, fs, f), 0.0, 0.01) << "passband deviation at " << f << " Hz"; + } + + // Stopband: at least the rated attenuation (1 dB grace) from the stopband + // edge out to well past the first few images. + for (double f = spec.stopband_hz; f <= 3.0 * fs; f += 250.0) { + EXPECT_LT(response_db(h, spec.num_phases, fs, f), -(spec.stopband_atten_db - 1.0)) + << "stopband leakage at " << f << " Hz"; + } + + // Transmission zeros at every k*fs: exact in exact arithmetic, so demand + // far below the rated stopband (double rounding measures ~-300 dB). + if (spec.image_zeros) { + for (int k = 1; k <= 3; ++k) { + EXPECT_LT(response_db(h, spec.num_phases, fs, static_cast(k) * fs), -150.0) + << "missing transmission zero at " << k << "*fs"; + } + } + } + + TEST(Kaiser, FastPrototypeMeetsSpec) { + check_prototype_meets_spec(k_fast, 48000.0); + } + + TEST(Kaiser, BalancedPrototypeMeetsSpec) { + check_prototype_meets_spec(k_balanced, 48000.0); + } + + TEST(Kaiser, TransparentPrototypeMeetsSpec) { + check_prototype_meets_spec(k_transparent, 48000.0); + } + + TEST(Kaiser, EconomyPrototypeMeetsSpec) { + check_prototype_meets_spec(k_economy, 48000.0); + } + + // The compensated presets must also hold their specs at scaled rates (the + // 16 kHz deployment path): normalized design, same numbers. + // The k*fs transmission zeros ARE branch-DC uniformity, stated in the + // frequency domain: with exact zeros, every polyphase branch's coefficient + // sum is identical (measured spread 1.8e-15 -- machine epsilon -- vs 4.7e-6 + // for the plain fast design, whose spread is its stopband leakage at fs). + TEST(Kaiser, CompensatedBranchSumsAreUniform) { + const proto_spec spec = k_balanced; + const std::size_t phases = std::bit_ceil(spec.num_phases); + std::vector h(phases * spec.taps_per_phase); + design_prototype_compensated(h, phases, (spec.passband_hz + spec.stopband_hz) / 48000.0, + kaiser_beta(spec.stopband_atten_db), spec.passband_hz / 48000.0); + double lo = 1e9; + double hi = -1e9; + for (std::size_t p = 0; p < phases; ++p) { + double sum = 0.0; + for (std::size_t t = 0; t < spec.taps_per_phase; ++t) { + sum += h[t * phases + p]; + } + lo = std::min(lo, sum); + hi = std::max(hi, sum); + } + EXPECT_LT(hi - lo, 1e-12); + EXPECT_NEAR(lo, 1.0, 1e-9); + } + + TEST(Kaiser, CompensatedSpecsHoldAt16k) { + check_prototype_meets_spec(k_balanced.scaled_to(16000.0), 16000.0); + check_prototype_meets_spec(k_economy.scaled_to(16000.0), 16000.0); + } + + // Non-power-of-two phase counts are first-class here (RatioTap designs at + // L = 147 and 160): the designer must meet spec without the bit_ceil + // rounding its historical consumer applied. 70 dB / 19->22.05 kHz is + // RatioTap's economy-profile shape for the 48->44.1 direction. + TEST(Kaiser, RationalPhaseCountMeetsSpec) { + const std::size_t phases = 147; + const double fs = 48000.0, pass = 19000.0, stop = 22050.0, atten = 70.0; + const std::size_t taps = estimate_taps(atten, (stop - pass) / fs); + std::vector h(phases * taps); + design_prototype(h, phases, (pass + stop) / fs, kaiser_beta(atten)); + for (double f = 0.0; f <= pass; f += 500.0) { + EXPECT_NEAR(response_db(h, phases, fs, f), 0.0, 0.05) << "passband deviation at " << f << " Hz"; + } + for (double f = stop; f <= 3.0 * fs; f += 250.0) { + EXPECT_LT(response_db(h, phases, fs, f), -(atten - 1.0)) << "stopband leakage at " << f << " Hz"; + } + } + + TEST(Kaiser, SolveDenseSolvesKnownSystem) { + // 3x3 that requires pivoting (zero leading entry); the solver works in + // place, so verify by residual against saved copies. + const std::vector m0{0.0, 2.0, 1.0, /**/ 1.0, 1.0, 1.0, /**/ 2.0, 0.0, 3.0}; + const std::vector rhs0{4.0, 6.0, 5.0}; + std::vector m = m0, rhs = rhs0, x(3); + solve_dense(m, rhs, x, 3); + for (std::size_t r = 0; r < 3; ++r) { + const double v = m0[3 * r] * x[0] + m0[3 * r + 1] * x[1] + m0[3 * r + 2] * x[2]; + EXPECT_NEAR(v, rhs0[r], 1e-12) << "row " << r; + } + } + +} // namespace diff --git a/tests/test_quantize.cpp b/tests/test_quantize.cpp new file mode 100644 index 0000000..df3a9d3 --- /dev/null +++ b/tests/test_quantize.cpp @@ -0,0 +1,101 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Contract battery for row-sum-preserving quantization, ported from the +// row-sum checks in SampleRateTap's test_fixed_point.cpp (there stated +// against the assembled polyphase bank; here against the utility itself, +// over rows of a genuinely designed prototype). + +#include +#include +#include + +#include + +#include "tap/dsp/kaiser.h" +#include "tap/dsp/quantize.h" + +namespace { + + using tap::dsp::quantize_row_preserving_sum; + using tap::dsp::sample_traits; + + // A realistic source: the RatioTap-shaped design (L = 160, plain Kaiser) + // whose branches this utility will quantize in production. + std::vector designed_prototype(std::size_t phases, std::size_t taps) { + std::vector h(phases * taps); + tap::dsp::design_prototype(h, phases, (19000.0 + 22050.0) / 48000.0, tap::dsp::kaiser_beta(70.0)); + return h; + } + + template + void check_rows_sum_exact() { + constexpr std::size_t k_phases = 160; + constexpr std::size_t k_taps = 32; + const auto proto = designed_prototype(k_phases, k_taps); + using coeff = typename sample_traits::coeff; + const double scale = sample_traits::k_coeff_scale; + + std::vector row(k_taps); + std::vector q(k_taps); + for (std::size_t p = 0; p < k_phases; ++p) { + double exact = 0.0; + for (std::size_t t = 0; t < k_taps; ++t) { + row[t] = proto[t * k_phases + p]; + exact += row[t] * scale; + } + quantize_row_preserving_sum(row, q); + std::int64_t sum = 0; + for (std::size_t t = 0; t < k_taps; ++t) { + sum += q[t]; + // Each tap stays within the rounding step plus at most the + // residual corrections: never grossly redistributed. + EXPECT_LE(std::abs(static_cast(q[t]) - row[t] * scale), 1.5) << "phase " << p << " tap " << t; + } + ASSERT_EQ(sum, std::llround(exact)) << "phase " << p; + } + } + + TEST(Quantize, RowSumsAreExactQ15) { + check_rows_sum_exact(); + } + + TEST(Quantize, RowSumsAreExactQ31) { + check_rows_sum_exact(); + } + + TEST(Quantize, FloatIsPlainConversion) { + const std::vector row{0.25, -0.125, 1.0, -0.9999, 0.0}; + std::vector q(row.size()); + quantize_row_preserving_sum(row, q); + for (std::size_t t = 0; t < row.size(); ++t) { + EXPECT_FLOAT_EQ(q[t], static_cast(row[t])); + } + } + + // The consequence the correction exists for: a unity-DC row quantizes to a + // row whose coefficient sum is exactly the format's unity, so DC gain + // through any phase deviates by at most one output LSB. + TEST(Quantize, UnityDcRowSumsToFormatUnity) { + constexpr std::size_t k_taps = 48; + // A smooth row normalized to sum exactly 1.0 in double. + std::vector row(k_taps); + double sum = 0.0; + for (std::size_t t = 0; t < k_taps; ++t) { + const double u = (static_cast(t) - 23.5) / 24.0; + row[t] = std::exp(-4.0 * u * u); + sum += row[t]; + } + for (auto& v : row) { + v /= sum; + } + std::vector q(k_taps); + quantize_row_preserving_sum(row, q); + std::int64_t qsum = 0; + for (const auto c : q) { + qsum += c; + } + EXPECT_EQ(qsum, 16384); // Q1.14 unity, exactly + } + +} // namespace diff --git a/tests/test_sample_traits.cpp b/tests/test_sample_traits.cpp new file mode 100644 index 0000000..57aff55 --- /dev/null +++ b/tests/test_sample_traits.cpp @@ -0,0 +1,111 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// Contract battery for the sample-format traits, ported and extended from +// SampleRateTap's test_fixed_point.cpp. Every number here is a documented +// contract point (Q formats, rounding mode, saturation, accumulator +// pre-shift); changing one is a breaking change for every consumer. + +#include +#include + +#include + +#include "tap/dsp/sample_traits.h" + +namespace { + + using q15 = tap::dsp::sample_traits; + using q31 = tap::dsp::sample_traits; + using f32 = tap::dsp::sample_traits; + + TEST(SampleTraits, CoefficientConversionRoundsAndSaturates) { + EXPECT_EQ(q15::make_coeff(0.0), 0); + EXPECT_EQ(q15::make_coeff(1.0), 16384); // Q1.14 + EXPECT_EQ(q15::make_coeff(-1.0), -16384); + EXPECT_EQ(q15::make_coeff(10.0), 32767); // saturates + EXPECT_EQ(q15::make_coeff(-10.0), -32768); + EXPECT_EQ(q31::make_coeff(1.0), 1073741824); // Q1.30 + EXPECT_EQ(q31::make_coeff(10.0), 2147483647); // saturates + EXPECT_FLOAT_EQ(f32::make_coeff(0.25), 0.25f); // unity scale, plain cast + } + + TEST(SampleTraits, RoundSatRoundsHalfAwayFromZero) { + using tap::dsp::detail::round_sat; + EXPECT_EQ(round_sat(0.5), 1); + EXPECT_EQ(round_sat(-0.5), -1); + EXPECT_EQ(round_sat(0.49), 0); + EXPECT_EQ(round_sat(-0.49), 0); + EXPECT_EQ(round_sat(1e9), std::numeric_limits::max()); + EXPECT_EQ(round_sat(-1e9), std::numeric_limits::min()); + } + + TEST(SampleTraits, FinalizeSaturates) { + // Far beyond full scale in the accumulator domain must clamp, not wrap. + EXPECT_EQ(q15::finalize(std::int64_t{1} << 40), 32767); + EXPECT_EQ(q15::finalize(-(std::int64_t{1} << 40)), -32768); + EXPECT_EQ(q31::finalize(std::int64_t{1} << 60), 2147483647); + EXPECT_EQ(q31::finalize(-(std::int64_t{1} << 60)), -2147483648LL); + } + + TEST(SampleTraits, FinalizeRoundsHalfUp) { + // Q29 -> Q15: the rounding constant is 2^13; exactly half rounds up, + // one below half rounds down, and the same holds on the negative side + // (round-half-up, not half-away: -0.5 LSB lands on 0). + EXPECT_EQ(q15::finalize((std::int64_t{1} << 13)), 1); + EXPECT_EQ(q15::finalize((std::int64_t{1} << 13) - 1), 0); + EXPECT_EQ(q15::finalize(-(std::int64_t{1} << 13)), 0); + EXPECT_EQ(q15::finalize(-(std::int64_t{1} << 13) - 1), -1); + // Q45 -> Q31 uses the identical constant. + EXPECT_EQ(q31::finalize((std::int64_t{1} << 13)), 1); + EXPECT_EQ(q31::finalize((std::int64_t{1} << 13) - 1), 0); + } + + TEST(SampleTraits, Q15MacIsExactInInt64) { + // Q0.15 x Q1.14 products are exact in int32 and summed without + // intermediate rounding: full-scale sample times unity coefficient. + const auto acc = q15::mac(0, std::int16_t{32767}, q15::make_coeff(1.0)); + EXPECT_EQ(acc, std::int64_t{32767} * 16384); + // finalize returns the sample: one rounding, at the end. + EXPECT_EQ(q15::finalize(acc), 32767); + } + + TEST(SampleTraits, Q31MacPreShiftsSixteenBits) { + // Q0.31 x Q1.30 = Q61 would overflow int64 after ~48 accumulations, so + // each product drops 16 bits before the add (Q45). Contract: exactly + // 16, no more (precision), no fewer (headroom). + const auto acc = q31::mac(0, std::int32_t{1} << 30, std::int32_t{1} << 30); + EXPECT_EQ(acc, std::int64_t{1} << 44); // (2^60) >> 16 + // Headroom check: ~80 full-scale accumulations stay inside int64. + std::int64_t big = 0; + for (int i = 0; i < 80; ++i) { + big = q31::mac(big, std::numeric_limits::max(), std::numeric_limits::max()); + } + EXPECT_GT(big, 0); // no wrap + EXPECT_EQ(q31::finalize(big), std::numeric_limits::max()); + } + + TEST(SampleTraits, FloatMacAccumulatesInDouble) { + static_assert(std::is_same_v); + // A double accumulator must hold a contribution a float one would drop: + // 1.0 + 2^-30 survives in double, vanishes in float. + const double acc = f32::mac(f32::mac(0.0, 1.0f, 1.0f), 0x1p-15f, 0x1p-15f); + EXPECT_GT(acc, 1.0); + EXPECT_DOUBLE_EQ(acc, 1.0 + 0x1p-30); + } + + TEST(SampleTraits, SilenceIsZero) { + EXPECT_EQ(q15::silence(), 0); + EXPECT_EQ(q31::silence(), 0); + EXPECT_EQ(f32::silence(), 0.0f); + } + + // The core concept is the kernels' requirement set; all three shipping + // formats satisfy it (also statically asserted in the header). + static_assert(tap::dsp::sample_type); + static_assert(tap::dsp::sample_type); + static_assert(tap::dsp::sample_type); + static_assert(!tap::dsp::sample_type); // no specialization on purpose: + // double is the golden-model/accumulator domain, not an I/O sample format. + +} // namespace