kithara dsp - zvuk/kithara GitHub Wiki

kithara-dsp

Documentation reviewed from source revision 4a470dd2c. This records the documented contract at that revision; it is not a new runtime validation. API and usage · All crates.

Vector DSP kernels over planar f32 slices (layout, sanitize, downmix, biquad cascades, window interpolation, the real FFT with bin magnitude and phase, autocorrelation and the sum of squares), plus the one import path for firewheel's fade, mix and smoothing helpers and for the biquad crate's filter designs.

Purpose

Layout and sanitize loops used to live inline in kithara-signal and kithara-decode, on the scalar path or on fast-interleave. kithara-dsp owns them as free functions whose backend the build target picks, so a hot loop is vectorised once and every consumer gets the same bits.

The offline analyzers followed the same pattern: kithara-beat and kithara-waveform each planned their own realfft transform and built their own Hann table, and kithara-waveform, kithara-analysis's producer and its beat analyzer each averaged channels in their own loop. They now take the spectrum from spectrum::Fft and the mono mix from downmix, so one length contract, one window and one downmix serve every analyzer.

Public surface

  • deinterleave_variable repeats the signature of fast_interleave's function for f32 and returns (). It writes the same bits (pinned against fast_interleave for 1–8 × 1–8 channels, with NaN payloads, infinities and subnormals in the data). Surplus output channels stay untouched, surplus input channels are ignored.
  • interleave_channel_major and deinterleave_channel_major address every plane of one strided channel-major slice (channel c fills slice[c * stride..(c + 1) * stride]), so any channel count runs without a slice of plane references and without allocating. They equal the fast_interleave functions over chunks_exact(stride); samples past the last whole plane are ignored or stay untouched.
  • Deliberate difference: where fast_interleave panics on a short slice, these handle the whole frames every used slice holds (the common prefix) and leave the rest untouched.
  • sanitize is the only function that rewrites values: NaN, ±infinity, subnormals and −0.0 become +0.0, bit-identical to kithara_signal::sanitize_sample (pinned by an exhaustive 2³² bit-pattern test in tests/crates/dsp). No other kernel sanitizes implicitly.
  • filter::Biquad runs biquad sections in cascade over planar channels. new allocates the state (and the vDSP setup on Apple); retune, process, settle and copy_state do not allocate. retune takes filter::Coefficients, rounds them to f32 and refuses a non-finite value or a pole outside the unit circle. A call of at least two frames whose input and output peaks stay at or below 1e-9 resets the state, so silence never decays into subnormals. The portable kernel also ends every call by zeroing each channel's section memory that has fallen below 1e-30, so a silent channel beside a playing one stops ringing before subnormals. settle(levels) puts every channel in the steady state for a constant input, so a consumer enters filtering without a step.
  • filter::{Coefficients, Type, Hertz, Errors} re-export the biquad crate's RBJ cookbook designs as the one import path; the crate designs, Biquad executes.
  • interp::interpolate(method, window, positions, output) reads a window at fractional positions with interp::Interpolation (linear, quadratic, Hermite, Watte). Interpolation::padding() names the frames a method reads before and after ⌊p⌋. A NaN position, or one outside before ≤ p < len − after, returns InterpError::OutOfWindow and leaves the output untouched.
  • interp::RateRamp is a read rate moving linearly from from to to over frames frames, then holding to. Offsets, the rest of a ramp and a block of positions come from the closed form, so a split block advances like one block and the ramp lands on to exactly.
  • downmix(interleaved, channels, mono) averages each whole interleaved frame into one mono sample and returns the frames written: the whole frames both slices hold. One channel copies, two run the pair kernel (l + r) · 0.5, and more sum the frame in channel order and multiply by 1/N; every path writes the bits of that scalar formula.
  • sum_squares(samples) reduces a slice to Σx²; an empty slice is +0.0, and the result stays within N·ε·Σx² of the f64 sum.
  • spectrum::FftLen is a length f·2ⁿ with f ∈ {1, 3, 5, 15} and n ≥ 4. FftLen::new is const and refuses any other length with SpectrumError::Length; bins() is N/2 + 1.
  • spectrum::Fft::new(len) plans the transform and a symmetric Hann table (denominator N − 1). spectrum(pools) builds a Spectrum — bins() real and imaginary parts and the frame buffer the transform works in — whose sample planes come from the caller's kithara-bufpool region, PoolError when they do not fit; off Apple it also allocates realfft's complex output and scratch. forward(frame, spectrum) windows the frame, zero-pads one shorter than N and writes bins 0..=N/2 unscaled on every backend; the DC and Nyquist bins have an imaginary part of exactly 0.0. A frame longer than N, or a Spectrum of another length, is SpectrumError::Shape. Fft is Send + Sync: one transform serves concurrent calls, each with its own Spectrum.
  • spectrum::magnitude(re, im, output) writes √(re² + im²) per bin, within 2 ε of hypot; unlike hypot it overflows once |re| or |im| exceeds 1.8e19, far above any audio bin. spectrum::phase writes atan2(im, re) within 1e-6 rad, and phase(+0, +0) is +0. Both handle the bins all three slices hold, return their count, and run one fearless_simd kernel on every platform, so a bin's value depends on the bin alone.
  • spectrum::Autocorrelation::new(len, pools) takes its zero padding (2N − 1 samples) once from the caller's region, PoolError when it does not fit. process(frame, output) writes lags 0..N of the first N samples of the frame, zero-extended, each divided by its N − lag products (the unbiased estimate), as many as output holds.
  • The FFT, the bin kernels and the autocorrelation compile only with the spectrum feature, off by default. Decode, effects, host, play, signal, resampler and warp depend on the crate and run no analysis, so they build no FFT. FftLen and SpectrumError are always available: kithara-waveform names the length in parameters it serves without its analyzer. kithara-beat and kithara-waveform turn the feature on with their dsp.
  • Nothing panics, and only filter::Biquad::new, spectrum::Fft::new and, off Apple, Fft::spectrum allocate; Fft::spectrum takes its sample planes, and spectrum::Autocorrelation::new its padding, from the caller's region. The crate is #![forbid(unsafe_code)] and denies clippy::indexing_slicing, clippy::arithmetic_side_effects and clippy::panic.

Backends

Backends are private modules chosen by cfg; callers never name one.

Channels macOS, iOS everything else wasm32
1 copy_from_slice copy_from_slice copy_from_slice
2 vDSP ztoc/ctoz through kithara-apple (the only unsafe) fearless_simd at the best level the CPU reports fearless_simd, +simd128 fixed at compile time
3+ scalar strided copy, four frames per iteration same same
sanitize fearless_simd fearless_simd fearless_simd
biquad cascade vDSP biquadm through kithara-apple fearless_simd direct form I, channels across lanes same
interpolation vDSP vlint/vqint (linear, quadratic); scalar kernel (Hermite, Watte) scalar kernel scalar kernel
peak (silence rule) vDSP maxmgv fearless_simd fearless_simd
real FFT vDSP DFT_zrop through kithara-apple, the frame split into even and odd samples in 64-byte-aligned pool planes realfft even-length transform planned with rustfft same
window vDSP vmul fearless_simd same
bin magnitude, bin phase fearless_simd (NEON); phase through a polynomial atan2 same same
autocorrelation vDSP conv fearless_simd dot product per lag same
sum of squares vDSP svesq fearless_simd same
downmix 2 channels: vDSP vasm; 1: copy; 3+: scalar sum × 1/N 2 channels: fearless_simd; 1 and 3+ as on Apple same

On x86 the SIMD level is detected once per process; elsewhere it is fixed at compile time. The crate's unit tests run every pair backend (portable at the native level, portable at Level::fallback(), Accelerate on Apple) against a scalar oracle for every size and offset, specials included.

Why the strided copy moves four frames per iteration: a loop that moves one sample per iteration runs at half speed whenever it straddles a 4096-byte page, and opt-level = "z" neither unrolls nor aligns it. Four samples per iteration amortize that fetch and roughly halve the cost on every target.

Why not BLAS for 3+ channels: Accelerate's cblas_scopy with a stride of 2 or more quiets a signaling NaN (sets bit 22), so it is not a bit-exact copy. At stride 1, and in vDSP ztoc/ctoz, bits are preserved. The strided path stays a plain scalar copy on every target.

Why the designs come from biquad but the execution does not: measured as the minimum of 3000 runs at opt-level = "z" on aarch64, vDSP biquadm runs 4–10× faster than the crate's direct form, and the portable cascade is faster from two channels on (stereo 4.0 against 6.0 ns per sample) because it puts channels across SIMD lanes; in mono it is 25 % slower. x86 AVX2 was not measured.

Why the interpolation kernel is scalar: the taps of each position are a gather, and moving them lane by lane into vectors costs more than the formula saves. The kernel evaluates in f64 and rounds to f32 once: evaluated in f32, Hermite on a window whose neighbours change sign overshot the peak by more than 4 ulp.

Why FftLen is vDSP's length set on every target: vDSP_DFT_zrop runs only f·2ⁿ. A length only the portable path accepted would build and pass on Linux and fail at construction on Apple, so both backends take the same set. Every such length is even, so the portable path always takes realfft's even-length transform.

Why the portable FFT is planned through rustfft::FftPlanner: RealFftPlanner hands out std::sync::Arc, which the arch check platform_layer_hygiene denies outside kithara-platform. A RealToComplexEven built from an FftPlanner is held by value.

Why the spectrum's sample planes come from the caller's region: analysis runs inside the caller's kithara-bufpool budget, and the perf.prefer-primitive-pool rule steers f32 scratch into the pool. The pool holds no Complex<f32>, so off Apple realfft's complex output and scratch are allocated with the Spectrum, once per spectrum. vDSP's real DFT also refuses planes off a 64-byte boundary (RealDft::ALIGN in kithara-apple), so on Apple the pool planes carry slack, vDSP works in their aligned windows, and the bins are then copied into the Spectrum: the result is a function of the frame alone, wherever the planes sit.

Why the Apple FFT is vDSP: on an Apple M3 Pro at opt-level = "z" the kernels bench measures Fft::forward at 1.31/2.91/4.00 µs for 1024/3072/4096 points, against 1.89/5.54/7.23 µs for the same windowed frame through realfft.

Why magnitude and phase run one portable kernel on Apple too: vDSP zvabs rounds its last elements by output address, so the same bins in another pool buffer gave other bits, and zvphas's vector path breaks the atan2 contract at zero-magnitude bins, whose phase feeds beat's novelty. The SDK promises placement independence for no vDSP routine; the portable √(mul_add) and polynomial atan2 have it by construction and give the same bits on every target. On the same machine they take 1.21 µs and 1.90 µs per 1024 bins.

Why a downmix above two channels stays scalar: a vector sum adds the channels in another order and moves the low bits. Offline analysis gains nothing there, and the scalar formula keeps the output bit-identical to the code it replaced.

wasm +simd128 minimum browsers

Safari / iOS 16.4, Chrome 91, Firefox 89. .cargo/config.toml sets +simd128 for every wasm32 build; the kithara-ffi web build passes --enable-simd to wasm-opt.

Integration

kithara-signal interleaves and deinterleaves its buffers through the layout functions; its pooled planar buffer goes through the channel-major ones, so no channel count allocates. kithara-decode sanitizes resampled planes.

kithara-resampler's Glide backend filters through filter::Biquad, interpolates through interp::interpolate and places its positions with interp::RateRamp.

kithara-beat's spectral detector reads its novelty from spectrum::Fft, magnitude and phase and its tempo period from spectrum::Autocorrelation. kithara-waveform sums band energy with sum_squares over the bins of an Fft; it, kithara-analysis's producer and its beat analyzer downmix through downmix.

Facade over firewheel

kithara_dsp::fade::FadeCurve and kithara_dsp::param::* (SmoothedParam, SmootherConfig, SmoothingFilter, SmoothingFilterCoeff, Mix, MixDSP and their constants) re-export firewheel's types as the one import path. The arch check firewheel_dsp_facade denies importing these modules from firewheel directly; node, event and buffer APIs stay direct firewheel imports. A re-export can later become a local type of the same name without touching consumers.

Crossfade gains in kithara-host (crossfader_gain) and kithara-play (CrossfadeSettings::gains) speak FadeCurve, so the host crossfader and the play crossfade share firewheel's curve, including its exact 0/1 ends below 1e-5 and above 0.99999.

no_panic spike (outcome)

#[no_panic] on the portable kernels was evaluated and rejected: the attribute moves the kernel body into a closure that does not inherit the fearless_simd target-feature context, so on x86_64 the AVX2 arm stops inlining its intrinsics. The panic guarantee is carried by the restriction lints above instead.

Build notes

  • [profile.dev.package.kithara-dsp] opt-level = 2: the kernels sit on kithara-decode's dev path, and opt-0 fearless_simd turns every lane op into a call.
  • fearless_simd is in hakari final-excludes so the dev-only force_support_fallback feature does not leak into release builds.
⚠️ **GitHub.com Fallback** ⚠️