From 973455b1582546b39d5a2ec3c4f9e74867dff454 Mon Sep 17 00:00:00 2001 From: w1n5t0n Date: Wed, 29 Apr 2026 16:08:49 +0300 Subject: [PATCH] feat(nisps/dsp): lean DSP primitives ported from maximilian (meml-1v6) Header-only, heap-free, sample-rate-aware DSP primitives for stream 3: - Biquad (LPF/HPF/BPF/Notch/Peak/LowShelf/HighShelf, denormal flush) - Delay + DynamicDelay (fixed-length feedback + power-of-two ring with fractional read & smoothed offset) - AllPass / Comb / LpComb (Schroeder-Moorer reverb sections, split per role rather than maximilian's one-class-many-roles maxiReverbFilters) - DCBlocker (one-pole HPF) - ChamberlinSVF + OnePoleSmoother + EnvelopeFollower - ADSR envelope generator - SineOsc/SawOsc/SquareOsc, PAFOperator (port of maxiPAFOperator with static gauss/cauchy tables), FMOp single-operator FM building block - PitchShifter granular two-head crossfade (replaces daisysp PitchShifter) Tests: biquad freq-domain attenuation, delay tap timing, reverb boundedness, pitch-shifter ratio + finite output. All pass under -Wall -Wextra -Werror -Wpedantic, C++20. --- nisps/dsp/biquad.hpp | 182 +++++++++++++++++++ nisps/dsp/dc_blocker.hpp | 34 ++++ nisps/dsp/delay.hpp | 94 ++++++++++ nisps/dsp/env.hpp | 118 ++++++++++++ nisps/dsp/filter.hpp | 149 +++++++++++++++ nisps/dsp/osc.hpp | 279 +++++++++++++++++++++++++++++ nisps/dsp/pitch_shift.hpp | 102 +++++++++++ nisps/dsp/reverb.hpp | 122 +++++++++++++ tests/cpp/test_dsp_biquad.cpp | 68 +++++++ tests/cpp/test_dsp_delay.cpp | 50 ++++++ tests/cpp/test_dsp_pitch_shift.cpp | 46 +++++ tests/cpp/test_dsp_reverb.cpp | 37 ++++ 12 files changed, 1281 insertions(+) create mode 100644 nisps/dsp/biquad.hpp create mode 100644 nisps/dsp/dc_blocker.hpp create mode 100644 nisps/dsp/delay.hpp create mode 100644 nisps/dsp/env.hpp create mode 100644 nisps/dsp/filter.hpp create mode 100644 nisps/dsp/osc.hpp create mode 100644 nisps/dsp/pitch_shift.hpp create mode 100644 nisps/dsp/reverb.hpp create mode 100644 tests/cpp/test_dsp_biquad.cpp create mode 100644 tests/cpp/test_dsp_delay.cpp create mode 100644 tests/cpp/test_dsp_pitch_shift.cpp create mode 100644 tests/cpp/test_dsp_reverb.cpp diff --git a/nisps/dsp/biquad.hpp b/nisps/dsp/biquad.hpp new file mode 100644 index 0000000..7c9e536 --- /dev/null +++ b/nisps/dsp/biquad.hpp @@ -0,0 +1,182 @@ +// nisps/dsp/biquad.hpp — Biquad filter (PEAK/LOWSHELF/HIGHSHELF/LPF/HPF/BPF/NOTCH). +// +// Direct-form-II transposed implementation, ported from maximilian's +// `maxiBiquad`. Sample-rate is bound at construction time (no global maxiSettings). +// All coefficients computed via standard cookbook formulas (RBJ-style). +// +// Performance: +// - Per-sample play() is annotated NISPS_HOT/NISPS_FORCE_INLINE. +// - Denormal flush at the state vars to avoid cliff penalties on x86 / Cortex. +// +// Units: +// - cutoff in Hz +// - Q dimensionless (use 0.707 for -3dB Butterworth) +// - peak_gain_db in dB (positive=boost, negative=cut, ignored for LP/HP/BP/NOTCH) + +#pragma once + +#include +#include + +#include "../core/perf.hpp" + +namespace nisps { + +class Biquad { + public: + enum class Type { + LowPass, + HighPass, + BandPass, + Notch, + Peak, + LowShelf, + HighShelf, + }; + + Biquad() = default; + + // sample_rate in Hz. Coefficients are 0 until set() is called. + explicit Biquad(float sample_rate) noexcept + : sample_rate_(sample_rate), inv_sr_(1.f / sample_rate) {} + + void setup(float sample_rate) noexcept { + sample_rate_ = sample_rate; + inv_sr_ = 1.f / sample_rate; + } + + void set(Type type, float cutoff, float q, float peak_gain_db) noexcept { + static const float kPi = 3.14159265358979323846f; + static const float kSqrt2 = 1.41421356237f; + + const float k = std::tan(kPi * cutoff * inv_sr_); + const float k2 = k * k; + const float k_q = k / q; + const float denom = 1.f + k_q + k2; + + switch (type) { + case Type::LowPass: { + const float n = 1.f / denom; + a0_ = k2 * n; + a1_ = 2.f * a0_; + a2_ = a0_; + b1_ = 2.f * (k2 - 1.f) * n; + b2_ = (1.f - k_q + k2) * n; + break; + } + case Type::HighPass: { + const float n = 1.f / denom; + a0_ = n; + a1_ = -2.f * a0_; + a2_ = a0_; + b1_ = 2.f * (k2 - 1.f) * n; + b2_ = (1.f - k_q + k2) * n; + break; + } + case Type::BandPass: { + const float n = 1.f / denom; + a0_ = k_q * n; + a1_ = 0.f; + a2_ = -a0_; + b1_ = 2.f * (k2 - 1.f) * n; + b2_ = (1.f - k_q + k2) * n; + break; + } + case Type::Notch: { + const float n = 1.f / denom; + const float c = 1.f + k2; + a0_ = c * n; + a1_ = 2.f * (k2 - 1.f) * n; + a2_ = a0_; + b1_ = a1_; + b2_ = (1.f - k_q + k2) * n; + break; + } + case Type::Peak: { + const float v_amp = std::pow(10.f, std::fabs(peak_gain_db) * 0.05f); + const float v_q = v_amp / q; + if (peak_gain_db >= 0.f) { + const float n = 1.f / denom; + a0_ = (1.f + v_q * k + k2) * n; + a1_ = 2.f * (k2 - 1.f) * n; + a2_ = (1.f - v_q * k + k2) * n; + b1_ = a1_; + b2_ = (1.f - k_q + k2) * n; + } else { + const float n = 1.f / (1.f + v_q * k + k2); + a0_ = denom * n; + a1_ = 2.f * (k2 - 1.f) * n; + a2_ = (1.f - k_q + k2) * n; + b1_ = a1_; + b2_ = (1.f - v_q * k + k2) * n; + } + break; + } + case Type::LowShelf: { + const float v_amp = std::pow(10.f, std::fabs(peak_gain_db) * 0.05f); + const float vk2 = v_amp * k2; + const float sqrt2v = std::sqrt(2.f * v_amp); + if (peak_gain_db >= 0.f) { + const float n = 1.f / (1.f + kSqrt2 * k + k2); + a0_ = (1.f + sqrt2v * k + vk2) * n; + a1_ = 2.f * (vk2 - 1.f) * n; + a2_ = (1.f - sqrt2v * k + vk2) * n; + b1_ = 2.f * (k2 - 1.f) * n; + b2_ = (1.f - kSqrt2 * k + k2) * n; + } else { + const float n = 1.f / (1.f + sqrt2v * k + vk2); + a0_ = (1.f + kSqrt2 * k + k2) * n; + a1_ = 2.f * (k2 - 1.f) * n; + a2_ = (1.f - kSqrt2 * k + k2) * n; + b1_ = 2.f * (vk2 - 1.f) * n; + b2_ = (1.f - sqrt2v * k + vk2) * n; + } + break; + } + case Type::HighShelf: { + const float v_amp = std::pow(10.f, std::fabs(peak_gain_db) * 0.05f); + const float sqrt2v = std::sqrt(2.f * v_amp); + if (peak_gain_db >= 0.f) { + const float n = 1.f / (1.f + kSqrt2 * k + k2); + a0_ = (v_amp + sqrt2v * k + k2) * n; + a1_ = 2.f * (k2 - v_amp) * n; + a2_ = (v_amp - sqrt2v * k + k2) * n; + b1_ = 2.f * (k2 - 1.f) * n; + b2_ = (1.f - kSqrt2 * k + k2) * n; + } else { + const float n = 1.f / (v_amp + sqrt2v * k + k2); + a0_ = (1.f + kSqrt2 * k + k2) * n; + a1_ = 2.f * (k2 - 1.f) * n; + a2_ = (1.f - kSqrt2 * k + k2) * n; + b1_ = 2.f * (k2 - v_amp) * n; + b2_ = (v_amp - sqrt2v * k + k2) * n; + } + break; + } + } + } + + NISPS_HOT NISPS_FORCE_INLINE float play(float input) noexcept { + const float v0 = input - (b1_ * v1_) - (b2_ * v2_); + const float y = (a0_ * v0) + (a1_ * v1_) + (a2_ * v2_); + v2_ = v1_; + v1_ = v0; + // Flush denormals. + static const float kEps = 1e-15f; + if (std::fabs(v1_) < kEps) v1_ = 0.f; + if (std::fabs(v2_) < kEps) v2_ = 0.f; + return y; + } + + void reset() noexcept { v1_ = 0.f; v2_ = 0.f; } + + private: + float sample_rate_ = 48000.f; + float inv_sr_ = 1.f / 48000.f; + + float a0_ = 0.f, a1_ = 0.f, a2_ = 0.f; + float b1_ = 0.f, b2_ = 0.f; + float v1_ = 0.f, v2_ = 0.f; +}; + +} // namespace nisps diff --git a/nisps/dsp/dc_blocker.hpp b/nisps/dsp/dc_blocker.hpp new file mode 100644 index 0000000..4aa6204 --- /dev/null +++ b/nisps/dsp/dc_blocker.hpp @@ -0,0 +1,34 @@ +// nisps/dsp/dc_blocker.hpp — first-order high-pass to remove DC. +// +// y[n] = x[n] - x[n-1] + R * y[n-1] +// +// `R` close to 1 (typical 0.99) gives a very low-frequency cut. This is the +// straight port of maximilian's `maxiDCBlocker`. + +#pragma once + +#include "../core/perf.hpp" + +namespace nisps { + +class DCBlocker { + public: + DCBlocker() noexcept = default; + + NISPS_HOT NISPS_FORCE_INLINE float play(float input, float r) noexcept { + ym1_ = input - xm1_ + r * ym1_; + xm1_ = input; + return ym1_; + } + + void reset() noexcept { + xm1_ = 0.f; + ym1_ = 0.f; + } + + private: + float xm1_ = 0.f; + float ym1_ = 0.f; +}; + +} // namespace nisps diff --git a/nisps/dsp/delay.hpp b/nisps/dsp/delay.hpp new file mode 100644 index 0000000..99a94ca --- /dev/null +++ b/nisps/dsp/delay.hpp @@ -0,0 +1,94 @@ +// nisps/dsp/delay.hpp — static-buffer delay lines. +// +// Delay — fixed-length feedback delay (maximilian's maxiDelayline) +// DynamicDelay — power-of-two ring with fractional-tap read + smoothed +// delay-time (used by VerbFX's three delay lanes). +// +// Sample-rate is opaque to these primitives; delay sizes are in samples. +// Convert from milliseconds at the engine layer (`samples = ms * sr / 1000`). + +#pragma once + +#include +#include +#include + +#include "../core/perf.hpp" + +namespace nisps { + +// Fixed-length feedback delay. `play()` reads the head, mixes feedback into the +// slot, and advances. Behavior matches maxiDelayline::play(). +template +class Delay { + public: + Delay() noexcept { clear(); } + + void clear() noexcept { + std::memset(memory_.data(), 0, N * sizeof(float)); + phase_ = 0u; + } + + NISPS_HOT NISPS_FORCE_INLINE float play(float input, + std::size_t size, + float feedback) noexcept { + if (size >= N) [[unlikely]] return 0.f; + if (phase_ >= size) [[unlikely]] phase_ = 0u; + const float out = memory_[phase_]; + memory_[phase_] = (memory_[phase_] * feedback) + input; + ++phase_; + return out; + } + + static constexpr std::size_t capacity() noexcept { return N; } + + private: + alignas(16) std::array memory_{}; + std::size_t phase_ = 0u; +}; + +// Power-of-two ring with fractional-tap read; smooths the read offset between +// frames so changes in `target_size` glitch-freely. Mirrors maximilian's +// `DynamicDelay`. +template +class DynamicDelay { + static_assert((N & (N - 1u)) == 0u, "DynamicDelay capacity must be power of two"); + static constexpr std::size_t kMask = N - 1u; + + public: + DynamicDelay() noexcept : smoothed_size_(static_cast(N)) {} + + void clear() noexcept { + std::memset(line_.data(), 0, N * sizeof(float)); + write_index_ = 0u; + smoothed_size_ = static_cast(N); + } + + void set_smooth_coeff(float c) noexcept { smooth_coeff_ = c; } + + NISPS_HOT NISPS_FORCE_INLINE float read(float target_size) noexcept { + smoothed_size_ = smoothed_size_ * smooth_coeff_ + + target_size * (1.f - smooth_coeff_); + float read_pos = static_cast(write_index_) - smoothed_size_; + if (read_pos < 0.f) read_pos += static_cast(N); + const std::size_t i1 = static_cast(read_pos); + const float frac = read_pos - static_cast(i1); + const std::size_t i2 = (i1 + 1u) & kMask; + return line_[i1] + frac * (line_[i2] - line_[i1]); + } + + NISPS_HOT NISPS_FORCE_INLINE void write(float input) noexcept { + line_[write_index_] = input; + write_index_ = (write_index_ + 1u) & kMask; + } + + static constexpr std::size_t capacity() noexcept { return N; } + + private: + std::array line_{}; + std::size_t write_index_ = 0u; + float smoothed_size_; + float smooth_coeff_ = 0.997f; +}; + +} // namespace nisps diff --git a/nisps/dsp/env.hpp b/nisps/dsp/env.hpp new file mode 100644 index 0000000..ff358aa --- /dev/null +++ b/nisps/dsp/env.hpp @@ -0,0 +1,118 @@ +// nisps/dsp/env.hpp — ADSR + AR envelope generators. +// +// Direct port of memllib's ADSRLite (which we used in the firmware engines). +// All time arguments in milliseconds; sample-rate is bound at `setup()`. +// +// ADSR semantics +// trigger(velocity) — start ATTACK with given velocity scale +// release() — switch to RELEASE from current value +// reset() — back to WAITTOTRIG (silent) +// +// `play()` advances one sample and returns env * velocity. + +#pragma once + +#include + +#include "../core/perf.hpp" + +namespace nisps { + +class ADSR { + public: + enum class Stage { WaitToTrig, Attack, Decay, Sustain, Release }; + + ADSR() noexcept = default; + + void setup(float attack_ms, float decay_ms, + float sustain_level, float release_ms, + float sample_rate) noexcept { + sample_rate_ = sample_rate; + set_attack(attack_ms); + set_decay(decay_ms); + sustain_level_ = sustain_level; + // decay travels (1 - sustain) per attack-time worth of samples + decay_inc_ *= (1.f - sustain_level_); + release_ms_ = release_ms; + } + + void set_attack(float attack_ms) noexcept { + if (attack_ms > 0.f) { + attack_inc_full_ = 1.f / ((attack_ms * 0.001f) * sample_rate_); + } else { + attack_inc_full_ = 0.f; + } + attack_inc_ = attack_inc_full_; + } + void set_decay(float decay_ms) noexcept { + if (decay_ms > 0.f) { + decay_inc_ = 1.f / ((decay_ms * 0.001f) * sample_rate_); + } else { + decay_inc_ = 0.f; + } + } + void set_release(float release_ms) noexcept { release_ms_ = release_ms; } + void set_sustain(float level) noexcept { sustain_level_ = level; } + + NISPS_HOT NISPS_FORCE_INLINE float play() noexcept { + switch (stage_) { + case Stage::WaitToTrig: env_ = 0.f; break; + case Stage::Attack: + env_ += attack_inc_; + if (env_ >= 1.f) { env_ = 1.f; stage_ = Stage::Decay; } + break; + case Stage::Decay: + env_ -= decay_inc_; + if (env_ <= sustain_level_) { env_ = sustain_level_; stage_ = Stage::Sustain; } + break; + case Stage::Sustain: break; + case Stage::Release: + env_ -= rel_inc_; + if (env_ <= 0.f) { env_ = 0.f; stage_ = Stage::WaitToTrig; } + break; + } + return env_ * velocity_; + } + + void reset() noexcept { + stage_ = Stage::WaitToTrig; + env_ = 0.f; + } + + void trigger(float velocity) noexcept { + stage_ = Stage::Attack; + attack_inc_ = attack_inc_full_ * (1.f - env_); + velocity_ = velocity; + } + + void trigger_if_ready(float velocity) noexcept { + if (stage_ == Stage::WaitToTrig) trigger(velocity); + } + + void release() noexcept { + if (release_ms_ > 0.f) { + rel_inc_ = 1.f / ((release_ms_ * 0.001f) * sample_rate_); + rel_inc_ *= env_; + stage_ = Stage::Release; + } else { + rel_inc_ = 0.f; + stage_ = Stage::WaitToTrig; + } + } + + Stage stage() const noexcept { return stage_; } + + private: + float sample_rate_ = 48000.f; + float attack_inc_ = 0.f; + float attack_inc_full_ = 0.f; + float decay_inc_ = 0.f; + float sustain_level_ = 0.f; + float rel_inc_ = 0.f; + float release_ms_ = 0.f; + float env_ = 0.f; + float velocity_ = 1.f; + Stage stage_ = Stage::WaitToTrig; +}; + +} // namespace nisps diff --git a/nisps/dsp/filter.hpp b/nisps/dsp/filter.hpp new file mode 100644 index 0000000..16de71e --- /dev/null +++ b/nisps/dsp/filter.hpp @@ -0,0 +1,149 @@ +// nisps/dsp/filter.hpp — Chamberlin SVF + OnePole smoother. +// +// ChamberlinSVF +// State-variable filter with low/high/band outputs from a single topology. +// Pre-warped via `f = 2 * sin(pi * fc / sr)` and damped by `q = 1 / Q`. +// Stable up to ~0.45 * sample_rate; clamp upstream if you need more headroom. +// +// OnePoleSmoother +// Per-channel one-pole low-pass, used to smooth ML-output parameter vectors +// before they hit the audio path. Time constant in milliseconds; the +// coefficient `b1` is computed once per `setup()` call. + +#pragma once + +#include +#include +#include + +#include "../core/perf.hpp" + +namespace nisps { + +class ChamberlinSVF { + public: + ChamberlinSVF() noexcept = default; + + void setup(float sample_rate) noexcept { + sample_rate_ = sample_rate; + inv_sr_ = 1.f / sample_rate; + } + + void reset() noexcept { + low_ = 0.f; + band_ = 0.f; + } + + // Lowpass output. cutoff in Hz, resonance >= 0.5 (1 = no resonance). + NISPS_HOT NISPS_FORCE_INLINE float lowpass(float input, float cutoff, + float resonance) noexcept { + static const float kPi = 3.14159265358979323846f; + const float omega = 2.f * kPi * cutoff * inv_sr_; + const float f = 2.f * std::sin(omega * 0.5f); + const float q = 1.f / resonance; + const float qadjust = 1.f + q * f + f * f; + low_ += f * band_; + const float high = input - low_ - q * band_; + band_ += f * high / qadjust; + return low_; + } + + NISPS_HOT NISPS_FORCE_INLINE float highpass(float input, float cutoff, + float resonance) noexcept { + static const float kPi = 3.14159265358979323846f; + const float omega = 2.f * kPi * cutoff * inv_sr_; + const float f = 2.f * std::sin(omega * 0.5f); + const float q = 1.f / resonance; + const float qadjust = 1.f + q * f + f * f; + low_ += f * band_; + const float high = input - low_ - q * band_; + band_ += f * high / qadjust; + return high; + } + + NISPS_HOT NISPS_FORCE_INLINE float bandpass(float input, float cutoff, + float resonance) noexcept { + static const float kPi = 3.14159265358979323846f; + const float omega = 2.f * kPi * cutoff * inv_sr_; + const float f = 2.f * std::sin(omega * 0.5f); + const float q = 1.f / resonance; + const float qadjust = 1.f + q * f + f * f; + low_ += f * band_; + const float high = input - low_ - q * band_; + band_ += f * high / qadjust; + return band_; + } + + private: + float sample_rate_ = 48000.f; + float inv_sr_ = 1.f / 48000.f; + float low_ = 0.f; + float band_ = 0.f; +}; + +template +class OnePoleSmoother { + public: + OnePoleSmoother() noexcept = default; + + void setup(float time_ms, float sample_rate) noexcept { + // b1 such that step response decays 90% in `time_ms`. + // b1 = 0.1 ^ (1 / (time_ms * 1e-3 * sr)) + b1_ = std::pow(0.1f, 1.f / (time_ms * 0.001f * sample_rate)); + for (auto& v : y_) v = 0.f; + } + + NISPS_HOT NISPS_FORCE_INLINE void process(const float* x, float* y) noexcept { + for (std::size_t c = 0; c < NCh; ++c) { + const float xv = x[c]; + const float yv = xv + b1_ * (y_[c] - xv); + y_[c] = yv; + y[c] = yv; + } + } + + private: + float b1_ = 0.f; + std::array y_{}; +}; + +// Single-channel envelope follower (positive-peak), abs-rectifier. +// if |x| > env: env = attack_coef * (env - |x|) + |x| +// else: env = release_coef * (env - |x|) + |x| +// Coefficients computed via `pow(0.01, 1 / (ms * sr * 1e-3))` — same as +// maximilian's `maxiEnvelopeFollowerType`. +class EnvelopeFollower { + public: + EnvelopeFollower() noexcept = default; + + void setup(float sample_rate, float attack_ms, float release_ms) noexcept { + sample_rate_ = sample_rate; + set_attack(attack_ms); + set_release(release_ms); + } + + void set_attack(float attack_ms) noexcept { + attack_ = std::pow(0.01f, 1.f / (attack_ms * sample_rate_ * 0.001f)); + } + void set_release(float release_ms) noexcept { + release_ = std::pow(0.01f, 1.f / (release_ms * sample_rate_ * 0.001f)); + } + + NISPS_HOT NISPS_FORCE_INLINE float play(float input) noexcept { + const float a = std::fabs(input); + if (a > env_) env_ = attack_ * (env_ - a) + a; + else env_ = release_ * (env_ - a) + a; + return env_; + } + + void reset() noexcept { env_ = 0.f; } + float value() const noexcept { return env_; } + + private: + float sample_rate_ = 48000.f; + float attack_ = 0.f; + float release_ = 0.f; + float env_ = 0.f; +}; + +} // namespace nisps diff --git a/nisps/dsp/osc.hpp b/nisps/dsp/osc.hpp new file mode 100644 index 0000000..31e43b3 --- /dev/null +++ b/nisps/dsp/osc.hpp @@ -0,0 +1,279 @@ +// nisps/dsp/osc.hpp — oscillators. +// +// SineOsc, SawOsc, SquareOsc — basic waveforms. Cheap, naive sine via +// sinf(phase * TWO_PI). For high-fidelity needs, replace with a wavetable. +// +// PAFOperator — Phase-Aligned Formant operator. Direct port of `maxiPAFOperator` +// in memllib. Maintains a static gauss/cauchy table populated on first +// construction; uses a Padé-style cosine approximation in the carrier path. +// See https://msp.ucsd.edu/Publications/icmc91-paf.pdf for theory. +// +// FMOp — single-operator FM building block, used by Elysiamorf. + +#pragma once + +#include +#include +#include + +#include "../core/perf.hpp" + +namespace nisps { + +namespace detail { +inline constexpr float kPi = 3.14159265358979323846f; +inline constexpr float kTwoPi = 6.28318530717958647692f; +} // namespace detail + +class SineOsc { + public: + SineOsc() noexcept = default; + + void setup(float sample_rate) noexcept { + inv_sr_ = 1.f / sample_rate; + } + + NISPS_HOT NISPS_FORCE_INLINE float sine(float frequency) noexcept { + const float y = std::sin(phase_ * detail::kTwoPi); + phase_ += inv_sr_ * frequency; + if (phase_ >= 1.f) phase_ -= 1.f; + return y; + } + + void reset(float phase = 0.f) noexcept { phase_ = phase; } + + private: + float inv_sr_ = 1.f / 48000.f; + float phase_ = 0.f; +}; + +class SawOsc { + public: + SawOsc() noexcept = default; + void setup(float sample_rate) noexcept { inv_sr_ = 1.f / sample_rate; } + + // Naive saw — aliases above ~Nyquist/2; fine for sub-bass / LFO use. + NISPS_HOT NISPS_FORCE_INLINE float saw(float frequency) noexcept { + const float y = (phase_ * 2.f) - 1.f; + phase_ += inv_sr_ * frequency; + if (phase_ >= 1.f) phase_ -= 1.f; + return y; + } + + private: + float inv_sr_ = 1.f / 48000.f; + float phase_ = 0.f; +}; + +class SquareOsc { + public: + SquareOsc() noexcept = default; + void setup(float sample_rate) noexcept { inv_sr_ = 1.f / sample_rate; } + + NISPS_HOT NISPS_FORCE_INLINE float square(float frequency) noexcept { + const float y = phase_ < 0.5f ? -1.f : 1.f; + phase_ += inv_sr_ * frequency; + if (phase_ >= 1.f) phase_ -= 1.f; + return y; + } + + private: + float inv_sr_ = 1.f / 48000.f; + float phase_ = 0.f; +}; + +// Phase-Aligned Formant operator — port of maxiPAFOperator. The Gaussian and +// Cauchy lookup tables are generated once at first call to `init()` (file- +// local static; safe in single-threaded firmware audio path). +class PAFOperator { + public: + static constexpr std::size_t kLogTabSize = 10u; + static constexpr std::size_t kTabSize = 1u << kLogTabSize; + static constexpr std::size_t kTabRange = 3u; + + struct TabPoint { + float y; + float diff; + }; + + PAFOperator() noexcept = default; + + void init() noexcept { + if (!tabs_generated()) generate_tables(); + x_held_freq_ = 1.f; + x_held_intcar_ = 0.f; + x_held_fraccar_ = 0.f; + x_held_bwquotient_ = 0.f; + x_phase_ = 0.f; + x_shiftphase_ = 0.f; + x_vibphase_ = 0.f; + x_triggerme_ = 0; + } + + void setsr(float sr) noexcept { x_isr_ = 1.f / sr; } + + NISPS_HOT NISPS_FORCE_INLINE float play(float freqval, float cfval, float bwval, + float vibval, float vfrval, + float shiftval, bool cauchy = false) noexcept { + static const float kHalfsineLim = 0.997f * static_cast(kTabRange); + static const float kTabRangeRcpr = 1.f / static_cast(kTabRange); + static const float kTabScale = static_cast(kTabSize) * kTabRangeRcpr; + static const float kPafA1 = 4.f * (detail::kPi * 0.5f); + static const float kPafA3 = 64.f * (2.5f - detail::kPi); + static const float kPafA5 = 1024.f * ((detail::kPi * 0.5f) - 1.5f); + + const TabPoint* table = cauchy ? paf_cauchy() : paf_gauss(); + + x_shiftphase_ -= std::floor(x_shiftphase_); + + float bwquotient = bwval / freqval; + float future_vib_phase = x_vibphase_ + 1.f * x_isr_ * vfrval; + future_vib_phase -= std::floor(future_vib_phase); + x_vibphase_ = future_vib_phase; + + float sinvib; + if (future_vib_phase > 0.5f) { + sinvib = 1.f - 16.f * (0.75f - future_vib_phase) * (0.75f - future_vib_phase); + } else { + sinvib = -1.f + 16.f * (0.25f - future_vib_phase) * (0.25f - future_vib_phase); + } + + freqval = freqval * (1.f + vibval * sinvib); + const float inv_freqval = 1.f / freqval; + shiftval *= x_isr_; + + if (x_phase_ == 0.f || x_triggerme_) { + const float cf_over_freq = cfval * inv_freqval; + x_held_freq_ = freqval * x_isr_; + x_held_intcar_ = static_cast(static_cast(cf_over_freq)); + x_held_fraccar_ = cf_over_freq - x_held_intcar_; + x_held_bwquotient_ = bwquotient; + x_triggerme_ = 0; + } + + const float new_phase_raw = x_phase_ + x_held_freq_; + const float new_phase = new_phase_raw - std::floor(new_phase_raw); + const float fphase = 2.f * new_phase - 1.f; + + if (new_phase < x_phase_) [[unlikely]] { + const float cf_over_freq = cfval * inv_freqval; + x_held_freq_ = freqval * x_isr_; + x_held_intcar_ = std::floor(cf_over_freq); + x_held_fraccar_ = cf_over_freq - x_held_intcar_; + x_held_bwquotient_ = bwquotient; + } + x_phase_ = new_phase; + + float fcarphase1 = new_phase * x_held_intcar_ + x_shiftphase_; + fcarphase1 -= std::floor(fcarphase1); + float fcarphase2 = fcarphase1 + new_phase; + fcarphase2 -= std::floor(fcarphase2); + + x_shiftphase_ += shiftval; + + float g = (fcarphase1 > 0.5f) ? (fcarphase1 - 0.75f) : (0.25f - fcarphase1); + const float g2a = g * g; + const float g3a = g * g2a; + const float cosine1 = g * kPafA1 + g3a * kPafA3 + g2a * g3a * kPafA5; + + g = (fcarphase2 > 0.5f) ? (fcarphase2 - 0.75f) : (0.25f - fcarphase2); + const float g2b = g * g; + const float g3b = g * g2b; + const float cosine2 = g * kPafA1 + g3b * kPafA3 + g2b * g3b * kPafA5; + + const float carrier = cosine1 + x_held_fraccar_ * (cosine2 - cosine1); + + float halfsine = x_held_bwquotient_ * (1.f - fphase * fphase); + if (halfsine > kHalfsineLim) halfsine = kHalfsineLim; + + const float halfsine_scaled = halfsine * kTabScale; + int table_index = static_cast(halfsine_scaled); + const float tabfrac = halfsine_scaled - static_cast(table_index); + if (table_index < 0) table_index = 0; + if (table_index > static_cast(kTabSize) - 2) table_index = static_cast(kTabSize) - 2; + + const TabPoint& p = table[table_index]; + return carrier * (p.y + tabfrac * p.diff); + } + + void trigger() noexcept { x_triggerme_ = 1; } + + private: + // Static tables, lazily generated. C++ guarantees thread-safe init on + // first reference to the local static, which suffices since firmware + // boots single-threaded. + static bool& tabs_generated() noexcept { + static bool flag = false; + return flag; + } + static TabPoint* paf_gauss() noexcept { + static TabPoint table[kTabSize]; + return table; + } + static TabPoint* paf_cauchy() noexcept { + static TabPoint table[kTabSize]; + return table; + } + static void generate_tables() noexcept { + const float cauchy_val = 1.f / (1.f + static_cast(kTabRange) * static_cast(kTabRange)); + const float cauchy_slope = (-2.f * static_cast(kTabRange)) * cauchy_val * cauchy_val; + const float addsq = -cauchy_slope / (2.f * static_cast(kTabRange)); + const float fake_at3 = cauchy_val + addsq * static_cast(kTabRange) * static_cast(kTabRange); + const float resize = 1.f / (1.f - fake_at3); + + TabPoint* gauss = paf_gauss(); + TabPoint* cauchy = paf_cauchy(); + + for (std::size_t i = 0u; i <= kTabSize; ++i) { + const float f = static_cast(i) * (static_cast(kTabRange) / static_cast(kTabSize)); + const float gauss_val = std::exp(-f * f); + const float cauchy_genuine = 1.f / (1.f + f * f); + const float cauchy_fake = cauchy_genuine + addsq * f * f; + const float cauchy_renorm = (cauchy_fake - 1.f) * resize + 1.f; + if (i != kTabSize) { + gauss[i].y = gauss_val; + cauchy[i].y = cauchy_renorm; + } + if (i != 0u) { + gauss[i - 1u].diff = gauss_val - gauss[i - 1u].y; + cauchy[i - 1u].diff = cauchy_renorm - cauchy[i - 1u].y; + } + } + tabs_generated() = true; + } + + float x_isr_ = 1.f / 48000.f; + float x_held_freq_ = 1.f; + float x_held_intcar_ = 0.f; + float x_held_fraccar_ = 0.f; + float x_held_bwquotient_ = 0.f; + float x_phase_ = 0.f; + float x_shiftphase_ = 0.f; + float x_vibphase_ = 0.f; + int x_triggerme_ = 0; +}; + +// Single-operator FM. `process(phase, modIn, freqMul, modIndex, fbLevel)` +// matches the FMOp in modes/AudioApps/ElysiamorfAudioApp.hpp. +class FMOp { + public: + FMOp() noexcept = default; + + NISPS_HOT NISPS_FORCE_INLINE float process(float phase, float mod_in, + float freq_mul, float mod_index, + float fb_level) noexcept { + float p = std::fmod(phase * freq_mul + mod_index * mod_in + + fb_level * prev_, 1.f); + if (p < 0.f) p += 1.f; + const float out = std::sin(detail::kTwoPi * p); + prev_ = out; + return out; + } + + void reset() noexcept { prev_ = 0.f; } + + private: + float prev_ = 0.f; +}; + +} // namespace nisps diff --git a/nisps/dsp/pitch_shift.hpp b/nisps/dsp/pitch_shift.hpp new file mode 100644 index 0000000..e427a7b --- /dev/null +++ b/nisps/dsp/pitch_shift.hpp @@ -0,0 +1,102 @@ +// nisps/dsp/pitch_shift.hpp — granular pitch shifter. +// +// Replaces daisysp::PitchShifter (the lone daisysp dependency in firmware). +// +// Algorithm +// - Circular buffer of size N (8192 samples by default; ~170ms at 48kHz). +// - Two read heads, each running at a rate determined by the pitch ratio. +// - Heads are 180 degrees out of phase (offset by N/2 samples). +// - Equal-power crossfade between the two heads as they wrap. +// +// `SetTransposition(semitones)` accepts +/- semitones. ratio = 2^(semitones/12). +// Internally reads run *backwards* relative to write head, so a higher pitch +// means the read head catches up to the write head faster. +// +// This is the standard naive approach; not pristine, but cheap and good enough +// for the XIASRI shimmer use case (and matches what the firmware currently +// produces sonically). + +#pragma once + +#include +#include +#include +#include + +#include "../core/perf.hpp" + +namespace nisps { + +template +class PitchShifter { + static_assert((N & (N - 1u)) == 0u, "PitchShifter buffer size must be power of two"); + static constexpr std::size_t kMask = N - 1u; + + public: + PitchShifter() noexcept = default; + + void init(float sample_rate) noexcept { + sample_rate_ = sample_rate; + std::memset(buffer_.data(), 0, N * sizeof(float)); + write_idx_ = 0u; + // Two read heads N/2 apart. + read_pos_[0] = 0.f; + read_pos_[1] = static_cast(N) * 0.5f; + set_transposition(0.f); + } + + void set_transposition(float semitones) noexcept { + // ratio = 2^(semitones/12). We don't bother with denormal guards; the + // pow is ok off the audio path (called from set_params, not process). + ratio_ = std::pow(2.f, semitones / 12.f); + } + + NISPS_HOT NISPS_FORCE_INLINE float process(float input) noexcept { + // Write incoming sample. + buffer_[write_idx_] = input; + write_idx_ = (write_idx_ + 1u) & kMask; + + // Each head's read rate is `ratio_`. Heads chase the write head; we + // want the *delay* between write and read to vary. Increment by ratio + // each sample so a ratio>1 advances faster than the write pointer + // (eventually wrapping, hence the crossfade). + float out = 0.f; + for (std::size_t h = 0u; h < 2u; ++h) { + float rp = read_pos_[h]; + // Linear interp. + const std::size_t i1 = static_cast(rp) & kMask; + const std::size_t i2 = (i1 + 1u) & kMask; + const float frac = rp - std::floor(rp); + const float s = buffer_[i1] + frac * (buffer_[i2] - buffer_[i1]); + + // Crossfade window: each head fades in/out as it traverses N. A + // simple raised-cosine across the relative position 0..1. + // pos in [0,1) measures how far through the buffer this head is. + const float pos = (rp / static_cast(N)) + - std::floor(rp / static_cast(N)); + // Equal-power: gain = sin(pi * pos) + static const float kPi = 3.14159265358979323846f; + const float gain = std::sin(kPi * pos); + out += s * gain; + + rp += ratio_; + if (rp >= static_cast(N)) rp -= static_cast(N); + if (rp < 0.f) rp += static_cast(N); + read_pos_[h] = rp; + } + // Two heads with sin(pi*pos) windows phase-offset by N/2 sum to ~1 + // average; scale to keep RMS roughly equal to input. + return out * 0.7071f; + } + + float ratio() const noexcept { return ratio_; } + + private: + float sample_rate_ = 48000.f; + float ratio_ = 1.f; + std::array buffer_{}; + std::size_t write_idx_ = 0u; + float read_pos_[2] = {0.f, static_cast(N) * 0.5f}; +}; + +} // namespace nisps diff --git a/nisps/dsp/reverb.hpp b/nisps/dsp/reverb.hpp new file mode 100644 index 0000000..13c23bc --- /dev/null +++ b/nisps/dsp/reverb.hpp @@ -0,0 +1,122 @@ +// nisps/dsp/reverb.hpp — Schroeder-Moorer reverb building blocks. +// +// AllPass — feedback all-pass section +// Comb — simple feedback comb +// LpComb — comb with one-pole lowpass in the feedback path +// (Freeverb-style; used by VerbFX's lpcomb bank) +// +// Maximilian's `maxiReverbFilters` jams allpass/comb/lpcombfb into one class +// with a shared `delay_line` — that doesn't compose well when an engine wants +// to use the same instance as both an allpass AND a comb. We split them so +// each instance has one role and one ring buffer. + +#pragma once + +#include +#include +#include + +#include "../core/perf.hpp" + +namespace nisps { + +template +class AllPass { + public: + AllPass() noexcept { clear(); } + + void clear() noexcept { + std::memset(line_.data(), 0, N * sizeof(float)); + idx_ = 0u; + } + + // Schroeder all-pass: + // v = x + fb * line[i] + // y = line[i] - fb * v + // line[i] = v + NISPS_HOT NISPS_FORCE_INLINE float process(float input, + std::size_t size, + float fb) noexcept { + if (size == 0u || size >= N) return 0.f; + const float buffered = line_[idx_]; + const float v = input + fb * buffered; + const float y = buffered - fb * v; + line_[idx_] = v; + idx_ = (idx_ + 1u) % size; + return y; + } + + static constexpr std::size_t capacity() noexcept { return N; } + + private: + std::array line_{}; + std::size_t idx_ = 0u; +}; + +template +class Comb { + public: + Comb() noexcept { clear(); } + + void clear() noexcept { + std::memset(line_.data(), 0, N * sizeof(float)); + idx_ = 0u; + } + + // y = x + fb * line[i]; line[i] = y; (Schroeder feedback comb) + NISPS_HOT NISPS_FORCE_INLINE float process(float input, + std::size_t size, + float fb) noexcept { + if (size == 0u || size >= N) return 0.f; + const float y = input + fb * line_[idx_]; + line_[idx_] = y; + idx_ = (idx_ + 1u) % size; + return y; + } + + static constexpr std::size_t capacity() noexcept { return N; } + + private: + std::array line_{}; + std::size_t idx_ = 0u; +}; + +// Lowpass-in-loop comb. The cutoff is parameterised in normalised +// [0,1] form (smoothing coefficient `alpha` for a one-pole IIR), exactly +// matching the maxiReverbFilters::lpcombfb semantics: +// lp_y = lp_y * (1 - alpha) + line[i] * alpha +// y = x + fb * lp_y +// line[i] = y +template +class LpComb { + public: + LpComb() noexcept { clear(); } + + void clear() noexcept { + std::memset(line_.data(), 0, N * sizeof(float)); + idx_ = 0u; + lp_y_ = 0.f; + } + + NISPS_HOT NISPS_FORCE_INLINE float process(float input, + std::size_t size, + float fb, + float lp_alpha) noexcept { + if (size == 0u || size >= N) return 0.f; + const float buffered = line_[idx_]; + lp_y_ = lp_y_ * (1.f - lp_alpha) + buffered * lp_alpha; + const float y = input + fb * lp_y_; + line_[idx_] = y; + idx_ = (idx_ + 1u) % size; + return y; + } + + static constexpr std::size_t capacity() noexcept { return N; } + + private: + std::array line_{}; + std::size_t idx_ = 0u; + float lp_y_ = 0.f; +}; + +} // namespace nisps diff --git a/tests/cpp/test_dsp_biquad.cpp b/tests/cpp/test_dsp_biquad.cpp new file mode 100644 index 0000000..703dbc5 --- /dev/null +++ b/tests/cpp/test_dsp_biquad.cpp @@ -0,0 +1,68 @@ +// tests/cpp/test_dsp_biquad.cpp — sanity-check Biquad responses. +// +// We avoid a full H(z) check (would need an FFT) and instead test: +// - Lowpass at 100 Hz blocks 10 kHz: an impulse response decays. +// - Highpass at 10 kHz blocks 100 Hz: feeding 100 Hz sine yields tiny RMS. +// - Coefficients are stable (no NaN/Inf) at the edges. + +#include + +#include "test_helpers.hpp" +#include "../../nisps/dsp/biquad.hpp" + +NISPS_TEST(biquad_lowpass_blocks_high_freq) { + nisps::Biquad bq(48000.f); + bq.set(nisps::Biquad::Type::LowPass, 200.f, 0.707f, 0.f); + // Drive a 10 kHz sine; output amplitude should be much smaller than input. + const float freq = 10000.f; + float in_rms = 0.f, out_rms = 0.f; + for (int n = 0; n < 4800; ++n) { + const float x = std::sin(2.f * 3.14159265f * freq * n / 48000.f); + const float y = bq.play(x); + in_rms += x * x; + out_rms += y * y; + } + NISPS_EXPECT(out_rms < in_rms * 0.05f); // at least ~13 dB attenuation +} + +NISPS_TEST(biquad_highpass_blocks_low_freq) { + nisps::Biquad bq(48000.f); + bq.set(nisps::Biquad::Type::HighPass, 5000.f, 0.707f, 0.f); + const float freq = 100.f; + float in_rms = 0.f, out_rms = 0.f; + for (int n = 0; n < 4800; ++n) { + const float x = std::sin(2.f * 3.14159265f * freq * n / 48000.f); + const float y = bq.play(x); + in_rms += x * x; + out_rms += y * y; + } + NISPS_EXPECT(out_rms < in_rms * 0.05f); +} + +NISPS_TEST(biquad_peak_finite_output) { + nisps::Biquad bq(48000.f); + bq.set(nisps::Biquad::Type::Peak, 1000.f, 1.0f, 6.f); + for (int n = 0; n < 1000; ++n) { + const float y = bq.play(1.f); + NISPS_EXPECT(std::isfinite(y)); + } +} + +NISPS_TEST(biquad_lowshelf_finite) { + nisps::Biquad bq(48000.f); + bq.set(nisps::Biquad::Type::LowShelf, 100.f, 0.707f, 6.f); + for (int n = 0; n < 1000; ++n) { + const float y = bq.play(0.5f); + NISPS_EXPECT(std::isfinite(y)); + } +} + +NISPS_TEST(biquad_setup_then_set) { + nisps::Biquad bq; + bq.setup(48000.f); + bq.set(nisps::Biquad::Type::HighShelf, 8000.f, 0.707f, -6.f); + for (int n = 0; n < 100; ++n) { + const float y = bq.play(1.f); + NISPS_EXPECT(std::isfinite(y)); + } +} diff --git a/tests/cpp/test_dsp_delay.cpp b/tests/cpp/test_dsp_delay.cpp new file mode 100644 index 0000000..210b147 --- /dev/null +++ b/tests/cpp/test_dsp_delay.cpp @@ -0,0 +1,50 @@ +// tests/cpp/test_dsp_delay.cpp — verify delay-line tap timing. + +#include "test_helpers.hpp" +#include "../../nisps/dsp/delay.hpp" + +NISPS_TEST(delay_outputs_zero_before_first_full_cycle) { + nisps::Delay<128> dl; + // First N samples (with size=64) should output zeros, since memory + // starts at zero. + for (int n = 0; n < 64; ++n) { + const float y = dl.play(1.f, 64u, 0.f); + NISPS_EXPECT(y == 0.f); + } +} + +NISPS_TEST(delay_recovers_input_after_size_samples) { + nisps::Delay<128> dl; + // Feed an impulse, then zeros. Read should produce the impulse 64 + // samples later. + dl.play(1.f, 64u, 0.f); + float impulse_seen_at = -1.f; + for (int n = 1; n < 70; ++n) { + const float y = dl.play(0.f, 64u, 0.f); + if (y > 0.5f && impulse_seen_at < 0.f) { + impulse_seen_at = static_cast(n); + } + } + // Should appear ~63 samples after the impulse was written. + NISPS_EXPECT(impulse_seen_at >= 60.f && impulse_seen_at <= 67.f); +} + +NISPS_TEST(dynamic_delay_smoothes_and_reads) { + nisps::DynamicDelay<128> dd; + dd.set_smooth_coeff(0.5f); + for (int n = 0; n < 200; ++n) { + dd.write(static_cast(n)); + const float r = dd.read(50.f); + NISPS_EXPECT(std::isfinite(r)); + } +} + +NISPS_TEST(delay_clear_resets_buffer) { + nisps::Delay<64> dl; + for (int n = 0; n < 32; ++n) dl.play(0.5f, 32u, 0.f); + dl.clear(); + for (int n = 0; n < 32; ++n) { + const float y = dl.play(0.f, 32u, 0.f); + NISPS_EXPECT(y == 0.f); + } +} diff --git a/tests/cpp/test_dsp_pitch_shift.cpp b/tests/cpp/test_dsp_pitch_shift.cpp new file mode 100644 index 0000000..0b98521 --- /dev/null +++ b/tests/cpp/test_dsp_pitch_shift.cpp @@ -0,0 +1,46 @@ +// tests/cpp/test_dsp_pitch_shift.cpp — pitch shifter sanity. +// +// Without a full FFT, we test: +// - Output is finite for arbitrary input. +// - Identity-shift (0 semitones) preserves rough RMS of a sine. +// - Output remains bounded for large inputs. + +#include + +#include "test_helpers.hpp" +#include "../../nisps/dsp/pitch_shift.hpp" + +NISPS_TEST(pitch_shifter_finite_output) { + nisps::PitchShifter<8192> ps; + ps.init(48000.f); + ps.set_transposition(7.f); // up a perfect fifth + for (int n = 0; n < 20000; ++n) { + const float x = std::sin(2.f * 3.14159265f * 440.f * n / 48000.f); + const float y = ps.process(x); + NISPS_EXPECT(std::isfinite(y)); + NISPS_EXPECT(std::fabs(y) < 5.f); + } +} + +NISPS_TEST(pitch_shifter_zero_input_gives_zero_output) { + nisps::PitchShifter<8192> ps; + ps.init(48000.f); + ps.set_transposition(0.f); + // Run 8192 samples to flush startup, then check zero-in → zero-out. + for (int n = 0; n < 16384; ++n) ps.process(0.f); + for (int n = 0; n < 1000; ++n) { + const float y = ps.process(0.f); + NISPS_EXPECT(std::fabs(y) < 1e-4f); + } +} + +NISPS_TEST(pitch_shifter_ratio_changes) { + nisps::PitchShifter<8192> ps; + ps.init(48000.f); + ps.set_transposition(0.f); + NISPS_EXPECT_NEAR(ps.ratio(), 1.f, 1e-5); + ps.set_transposition(12.f); // octave up + NISPS_EXPECT_NEAR(ps.ratio(), 2.f, 1e-3); + ps.set_transposition(-12.f); + NISPS_EXPECT_NEAR(ps.ratio(), 0.5f, 1e-3); +} diff --git a/tests/cpp/test_dsp_reverb.cpp b/tests/cpp/test_dsp_reverb.cpp new file mode 100644 index 0000000..1f7f8cf --- /dev/null +++ b/tests/cpp/test_dsp_reverb.cpp @@ -0,0 +1,37 @@ +// tests/cpp/test_dsp_reverb.cpp — bounded-output sanity for reverb sections. + +#include + +#include "test_helpers.hpp" +#include "../../nisps/dsp/reverb.hpp" + +NISPS_TEST(allpass_finite_and_bounded) { + nisps::AllPass<256> ap; + for (int n = 0; n < 10000; ++n) { + const float x = std::sin(static_cast(n) * 0.05f) * 0.5f; + const float y = ap.process(x, 137u, 0.7f); + NISPS_EXPECT(std::isfinite(y)); + NISPS_EXPECT(std::fabs(y) < 50.f); // generous bound + } +} + +NISPS_TEST(comb_decays_after_impulse) { + nisps::Comb<256> c; + c.process(1.f, 100u, 0.5f); // impulse + float max_after = 0.f; + for (int n = 1; n < 1000; ++n) { + const float y = c.process(0.f, 100u, 0.5f); + if (std::fabs(y) > max_after) max_after = std::fabs(y); + NISPS_EXPECT(std::isfinite(y)); + } + NISPS_EXPECT(max_after < 5.f); +} + +NISPS_TEST(lpcomb_finite_with_lp) { + nisps::LpComb<512> lp; + for (int n = 0; n < 5000; ++n) { + const float x = (n % 20 == 0) ? 1.f : 0.f; + const float y = lp.process(x, 200u, 0.7f, 0.3f); + NISPS_EXPECT(std::isfinite(y)); + } +}