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 seven primitives, plus the FIR substrate — the shared design-math / sample-format / kernel layer under SampleRateTap and RatioTap:
include/tap/dsp/fft.h wraps the vendored Ooura split-radix real
FFT behind a small, well-specified interface, with
optional, mutually-exclusive float32 backends that re-present the exact
same numeric contract for speed on specific hardware:
| Backend | Build option | Target | Notes |
|---|---|---|---|
| Ooura (default) | — | everywhere | the golden model; double is always Ooura |
| CMSIS-DSP Helium | TAP_DSP_FFT_CMSIS |
bare-metal Cortex-M55 (MVE) | ~3× fewer instructions/transform |
| Apple vDSP | TAP_DSP_FFT_ACCELERATE |
macOS / Apple Silicon | ~3× faster/transform |
The two float32 backends conjugate imaginary bins and rescale so every
intermediate spectrum matches the Ooura build to single-precision rounding —
so the whole double-precision test battery stays a valid oracle for the
accelerated float paths. tests/test_fft_backend.cpp pins each backend to
Ooura's rdft_f bin-for-bin at the certified geometries (N = 512, 2048).
#include "tap/dsp/fft.h"
tap::dsp::real_fft fft(1024); // double, the desktop/golden profile
tap::dsp::real_fft32 fft32(1024); // float, the embedded / accelerated profile
std::vector<double> x(1024, 0.0);
fft.forward_inplace(x.data()); // packed spectrum, W = exp(+2πi/N)
fft.inverse(x.data(), x.data()); // out-of-place inverse, normalized (2/N)Key contract points (full detail in the header docstring):
- Packing (N/2 + 1 bins):
data[0]= DC real,data[1]= Nyquist real,data[2k]/data[2k+1]= bin k real/imag for1 ≤ k < N/2. - Sign convention
W = exp(+2πi/N)— imaginary parts are conjugated relative to the engineering-convention DFT. Consistent across every operand, so spectral products (fast convolution, adaptive-filter regressors) are unaffected; conjugate only when importing spectra computed elsewhere. - Normalization:
*_inplaceinverse is unnormalized (multiply by2/N); the out-of-placeinverse()applies the2/Nfor you. - Transforms are
noexceptand allocation-free after construction — real-time safe. Size must be a power of two,≥ 4, fixed at construction.
include/tap/dsp/yin.h implements the time-domain YIN estimator (de Cheveigné
& Kawahara 2002, steps 1–5): squared-difference function, cumulative-mean
normalization, absolute threshold with local-minimum descent, and parabolic
interpolation for sub-sample period precision. Header-only, allocation-free
after construction, noexcept on the analysis path.
#include "tap/dsp/yin.h"
tap::dsp::yin det(800, 20, 800); // window, tau_min, tau_max — double golden model
tap::dsp::yin32 det32(800, 20, 800); // float, the embedded profile
const auto r = det.analyze(frame); // frame_size() == window + tau_max samples
if (r.voiced()) {
const double freq = sample_rate / r.period; // fractional-sample period
}Key contract points (full detail in the header docstring):
- Geometry fixed at construction: integration
window, searched lag range[tau_min, tau_max](bound from your frequency range asτ = sr / f), withwindow ≥ tau_max;analyze()readswindow + tau_maxsamples, oldest first. - Result: fractional period in samples (0 = unvoiced) plus the normalized aperiodicity at the chosen lag (global minimum when unvoiced).
- Threshold: the paper's absolute threshold on the normalized difference, default 0.1, settable at runtime.
- Both precisions run the identical algorithm; the test battery pins sine/ sawtooth accuracy (sub-cent in double), octave robustness, unvoiced rejection, and float/double agreement. The difference-function inner loop is the designated Helium-MVE / HVX backend candidate behind this same contract, mirroring the FFT's golden-model-plus-backends pattern.
include/tap/dsp/psola.h is the real-time TD-PSOLA resynthesis stage:
Hann grains two source periods long, extracted at period-spaced analysis marks
and overlap-added at period/ratio-spaced synthesis marks with sub-sample
(Hermite) placement. Detection-agnostic — the caller supplies the period (from
tap::dsp::yin or any tracker); fixed latency of 2 * max_period + 2 samples.
tap::dsp::psola shifter(900); // deepest period it will be given
double y = shifter.process(x, period, ratio); // per sample; ratio 2 = octave upKnow what PSOLA is: it resamples the source's spectral envelope at the new harmonic spacing — which is why it preserves formants on voice, and why a pure tone shifted far from any new harmonic thins toward silence. Feed it harmonic-rich material; both behaviors are pinned by the tests.
include/tap/dsp/pvoc.h is an STFT pitch shifter (Hann, 4× overlap, built on
tap::dsp::real_fft so the float profile rides the vDSP/CMSIS backends) using
Laroche–Dolson-style peak-region shifting: each spectral peak's region is
translated rigidly by an integer bin offset and rotated by one accumulated
residual phase, so phase relationships across the peak stay intact. At
ratio 1 the output reconstructs the input's waveform delayed by exactly one
FFT frame (pinned by the tests). Latency = the FFT size (1024 default).
tap::dsp::pvoc shifter(1024);
shifter.set_formant(true); // optional LPC formant preservation
double y = shifter.process(x, ratio); // per sample; ratio sampled per hopOptional formant preservation (set_formant) uses the classic source-filter
method: an LPC spectral envelope (autocorrelation + Levinson–Durbin, order 48)
per analysis frame, with every relocated bin rescaled by
envelope(target)/envelope(source) — the excitation moves, the envelope stays.
At ratio 1 the correction is exactly unity, so the identity contract holds
either way.
include/tap/dsp/log_mel.h is the streaming analysis front end of a keyword
spotter: windowed real FFT (on tap::dsp::real_fft, so the float profile
rides the vDSP/CMSIS backends), power spectrum, triangular mel filterbank,
then a floored affine log10 or per-channel energy normalization (PCEN,
Wang et al. 2017). The header owns the formula-level contract as
numbers — HTK mel, unit-peak triangles, periodic Hann or sqrt-Hann, FFT
zero-padded at the end of the frame, frame t ending at sample
(t+1)*hop - 1, DC excluded, the PCEN recursion and its first-frame
priming — and stamps it with k_contract_version. Everything a trainer
tunes (band count, fmin/fmax, log floor/shift/scale, every PCEN parameter,
pre-emphasis) is runtime geometry in log_mel_geometry, carried by a trained
model, so retraining never touches the header. Reference geometry: 16 kHz,
400 / 160 / 512, 40 bands 20–7600 Hz. Latency = the frame length.
tap::dsp::log_mel_geometry g; // the reference geometry
g.pcen.enabled = true; // or leave the plain-log path
tap::dsp::log_mel32 fe(g); // float embedded profile; log_mel is the double golden model
std::vector<float> feats(fe.frames_for(n) * fe.bands());
size_t frames = fe.process(x, n, feats.data(), fe.frames_for(n)); // any chunking, same featuresPinned by tests/test_log_mel.cpp against the committed numpy restatement
(tools/reference/make_frontend_reference.py → tests/reference/frontend_vectors.h
at the reference geometry and frontend_vectors_tuned.h at a geometry with
every runtime field off its default; the script's Geometry mirrors
log_mel_geometry field for field and is the only numpy copy of these
formulas in the family — MuTap's KWS feature module imports it through the
submodule for its parity self-check, and trains on the C++ front end through
the C ABI): both paths sample-for-sample at both geometries, chunking
invariance, alignment and latency, the filterbank formulas, PCEN's gain
tracking on a level step and its reset semantics, and float/double agreement
as a measured number (6.7e-7 log, 5.3e-6 PCEN, pinned at 2×). No double
arithmetic on the float path: the RP2350's Cortex-M33 has no FP64.
include/tap/dsp/decimate.h is the host-rate stage in front of the 16 kHz
front end: basic_decimator<Sample, M> for M = 2, 3, 6 (32 / 48 / 96 kHz in),
in RatioTap's pattern — ratio as a type, Kaiser-windowed sinc from
kaiser.h with the cutoff at the output Nyquist and DC gain exactly 1,
fir_kernels.h's dot_row over the sample_traits.h formats (double golden,
float embedded, Q15 / Q31 with row-sum-preserving quantization). It is deliberately not
RatioTap, whose charter is 44.1 ↔ 48 only; a 44.1 kHz host composes RatioTap's
44.1 → 48 in front of the by-3 stage. Odd tap counts, integer group delay
(taps - 1) / 2, one output as input k*M arrives.
| profile | stopband | passband | taps by 2 / 3 / 6 |
|---|---|---|---|
| economy | 70 dB | 7000 Hz | 81 / 121 / 239 |
| transparent | 100 dB | 7600 Hz | 259 / 389 / 773 |
tap::dsp::decimate_by_3 dec; // 48 kHz -> 16 kHz, economy, float
std::vector<float> y(dec.outputs_for(n));
dec.process(x, n, y.data()); // noexcept, allocation-free, chunking-invariantPinned by tests/test_decimate.cpp: the tap counts against the searched
minima, the float output sample-for-sample against the numpy reference,
the passband and stopband numbers measured from the shipped coefficients,
unity DC (exact in Q15), the group delay, and Q15 tracking float within the
format's floor.
include/tap/dsp/nn.h is the arithmetic of a small recurrent gain network:
basic_dense<Sample> (y = act(W x + b), row-major [out x in], linear /
tanh / sigmoid) and basic_gru<Sample> (Cho et al. 2014 in PyTorch's
nn.GRU convention, gates ordered r, z, n along the 3*hidden axis). Lifted
from MuTap's learned residual suppressor at its wake-word plan's M3; the
keyword spotter is a different head on the same layers. Weights are stored
as float32 whatever the profile, the storage precision of a trained model's
file, and converted to Sample at the point of use; every dot product
accumulates in Sample bias-first in ascending input order, so the float
profile contains no double arithmetic. Weight vectors are moved in and
owned (a copied layer is a deep copy); apply() / step() are noexcept and
allocation-free. k_contract_version stamps the formulas.
using namespace tap::dsp::nn;
dense32 din(w_in, b_in, 64, 56, activation::tanh); // float embedded profile; dense/gru are the double golden model
gru32 cell(w_ih, w_hh, b_ih, b_hh, 96, 64);
dense32 dout(w_out, b_out, 26, 96, activation::sigmoid);
din.apply(features, hidden); cell.step(hidden); dout.apply(cell.state(), gains);Pinned by tests/test_nn.cpp: the layout on hand-computed numbers, the
activation forms, the gate order by isolating each block through its
biases, the GRU formula against an independent long-double restatement
that sums in the opposite order, reset and copy semantics, noexcept, and
float-tracks-double as a measured number at the suppressor's geometry. The
end-to-end oracle stays in MuTap: its Python parity CI job and its
suppressor's cross-precision pin must be unchanged by the promotion.
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,
plus the FFT's butterfly arithmetic trait (fft/fft_arith.h), which is built
over the same sample formats and documented here until the Stage 3b README
rewrite moves it to the FFT section's profiles table.
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.
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.
The family's sample-format substrate: how each sample type stores
coefficients, accumulates dot products, and rounds/saturates back to samples.
double is a sample format because it is the golden model of every
primitive: a traits-based primitive instantiates its reference profile
through the same substrate as its embedded profiles, and the cross-precision
pins measure float/Q15/Q31 against it.
| Type | Coefficients | Accumulation | Output |
|---|---|---|---|
double |
double | double | identity (the golden model) |
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 | single Q45→Q31 round-half-up, saturating |
The Q ladder is spelled as named constants on each fixed-point
specialization (k_sample_frac_bits, k_coeff_frac_bits,
k_accum_pre_shift) with the accumulator format and the single 14-bit
finalize shift derived from them and static_asserted; k_coeff_scale is
2^k_coeff_frac_bits, and k_is_fixed_point is what selects the
fixed-point algorithm in quantize.h. Every member is constexpr, and the
sample_type concept requires a value-initialized accumulator to be the
additive identity.
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.
The sibling trait the fixed-point real FFT is written against (its docstrings
are that kernel's specification — shift-before-butterfly, the magnitude bound,
the rounding count per complex product, the BFP headroom rule, the Q31 input
pre-shift, the twiddle generator; tests/test_fft_arith.cpp pins every
number). Documented here until the Stage 3b README rewrite; it moves to the
FFT section's profiles table then. One int32 kernel serves both fixed profiles: fft_arith<std::int32_t>
carries mul_coeff (int32 × Q1.30 → int64, >> 30 with one round-half-up,
saturating), saturating add / sub, shr_round (round-half-up), and
headroom_bits (the block's shared redundant sign bits, via
std::countl_zero); fft_arith<std::int16_t> is the Q15 I/O width —
widen (<< 14, two guard bits) and narrow (round-half-up, saturating) —
and names the int32 trait as its work. Twiddles are
sample_traits<std::int32_t>::coeff (Q1.30) for both fixed profiles, and 1.0
is representable. The float / double specializations are the same names
over plain arithmetic. Everything is constexpr and noexcept.
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.
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). Selected by
the trait's k_is_fixed_point: plain conversion for double and float. A tap
at the format's rail is never wrapped: each step goes to the largest-remainder
tap that can still move in that direction, so the sum is preserved whenever
one exists. Design-time code.
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.
The notebooks drive the actual shipping C++ through the C ABI in
tools/capi/ (ctypes bridge: notebooks/dsptap_py.py, which builds
build_capi/ on first import; the bridge also exposes LogMel and
Decimator for MuTap's keyword-spotting notebooks). They are committed
executed; re-execute with
jupyter nbconvert --to notebook --execute --inplace notebooks/<name>.ipynb
when a primitive's behavior changes.
notebooks/pitchshift.ipynb measures the three pitch primitives. It
documents the two findings from the primitives' development: PSOLA's
envelope-resampling nature (why it preserves formants and why a pure tone
shifted an octave thins out), and the measured level collapse of naive
phase-vocoder bin remapping vs the shipping peak-locked design — plus the LPC
formant-preservation demo.
notebooks/fft.ipynb measures the real FFT through dsptap_py.RealFFT
(profiles "double" and "float"; unpack/pack convert the header's
packed exp(+i) layout to and from numpy.fft.rfft's convention): the packing
and sign contract against numpy, the float32 profile's per-bin error against
the double golden model vs N — on the profile's own arithmetic via the raw
in-place entry points, not a double round trip — and the round-trip error of
both profiles. Its fixed-point section (Q15/Q31 noise floors under Welch's
model) is a designed placeholder until Stage 3b/3c of
docs/audit-fft-and-code-smells.md lands.
Standalone (builds the Ooura path, plus vDSP on macOS, and runs the tests):
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build
ctest --test-dir build --output-on-failureCI also runs an emulation-sized selection of the battery
(tests/bare_metal_main.cpp) bare-metal under QEMU on four Cortex-M legs:
cortex-m4-softfp and cortex-m4f (cmake/arm-cortex-m4-mps2.cmake, the
FPU flavour selected by -DTAP_DSP_M4_FPU=ON), cortex-m33
(cmake/arm-cortex-m33-mps2.cmake) and cortex-m55
(cmake/arm-cortex-m55-mps3.cmake, where the CMSIS-DSP Helium FFT backend is
ON and its parity suite runs against Ooura). Every suite compiled into a test
executable runs on the target unless excluded by name in
tests/CMakeLists.txt (a negative filter; each exclusion is a budget note).
To run one locally, with arm-none-eabi-g++ and qemu-system-arm on PATH:
cmake -S . -B build-m33 -DCMAKE_BUILD_TYPE=MinSizeRel \
-DCMAKE_TOOLCHAIN_FILE=cmake/arm-cortex-m33-mps2.cmake
cmake --build build-m33
ctest --test-dir build-m33 --output-on-failureBenchmarks and the per-target instruction-count ratchet build with
-DTAP_DSP_BUILD_BENCH=ON; policy and workflow in bench/README.md.
add_subdirectory(submodules/dsptap) # or however it is pinned
target_link_libraries(my_dsp PRIVATE tap::dsp)tap::dsp is an INTERFACE target (the headers + the compiled Ooura
static lib tap::dsp_fft); it does not build the tests when added as a
subdirectory (TAP_DSP_BUILD_TESTS defaults OFF unless top-level). The
per-platform float32 backend defaults follow the target: vDSP on Apple, CMSIS
on the bare-metal M55 profile, Ooura elsewhere — override with
-DTAP_DSP_FFT_ACCELERATE=OFF etc.
This code was carried, with textually identical vendored Ooura sources, inside
both MuTap (adaptive filtering) and AmbiTap (ambisonics / binaural
convolution). The C++ wrappers had begun to diverge — MuTap grew the templated
basic_real_fft<Sample> and the CMSIS/vDSP backends; AmbiTap kept an older
double-engine wrapper with no backends — so a bug fix or a new backend in one
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 and
third_party/cmsis-dsp/VENDOR.md for the
vendored-code provenance and licenses. The design note for the planned C++20
port of the same split-radix transform (and the Q15 / Q31 profiles) is
docs/fft-design.md, filled in as each stage lands; the
plan of record, docs/audit-fft-and-code-smells.md, lands with the plan PR
(#25). Nothing from either has shipped yet.
DspTap's own code is MIT (LICENSE). Vendored third-party code keeps its own
license. The Ooura FFT is under its author's own terms, stated in
third_party/ooura/readme.txt: "You may use, copy, modify this code for any
purpose and without fee. You may distribute this ORIGINAL package." — a grant
of use, copying and modification, and of distribution of the original package.
DspTap ships that file today with the notice attached, and the planned C++
port is a derivative work whose redistribution relies on the modification
grant (SPDX LicenseRef-Ooura AND MIT for the port header, with the notice
text in LICENSES/LicenseRef-Ooura.txt; the readme stays at
third_party/ooura/readme.txt permanently). CMSIS-DSP / CMSIS-Core are
Apache-2.0 with SPDX headers retained in every file. The canonical statement,
and the maintainer's reading of what it covers — a judgement call, not legal
advice — is NOTICE.md.