w1n5t0n/memlnaut-nisps
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

feat(nisps/dsp): lean DSP primitives ported from maximilian (meml-1v6)

Browse source 
Header-only, heap-free, sample-rate-aware DSP primitives for stream 3:
- Biquad (LPF/HPF/BPF/Notch/Peak/LowShelf/HighShelf, denormal flush)
- Delay<N> + DynamicDelay<N> (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<NCh> + EnvelopeFollower
- ADSR envelope generator
- SineOsc/SawOsc/SquareOsc, PAFOperator (port of maxiPAFOperator with
  static gauss/cauchy tables), FMOp single-operator FM building block
- PitchShifter<N> 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.
This commit is contained in:
w1n5t0n 2026-04-29 16:08:49 +03:00
parent 825ed6ad33
commit 973455b158
12 changed files with 1281 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

182
nisps/dsp/biquad.hpp Normal file
View file

@ -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 <cmath>
#include <cstddef>
#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

34
nisps/dsp/dc_blocker.hpp Normal file
View file

@ -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

94
nisps/dsp/delay.hpp Normal file
View file

@ -0,0 +1,94 @@
// nisps/dsp/delay.hpp — static-buffer delay lines.
//
// Delay<N> — fixed-length feedback delay (maximilian's maxiDelayline)
// DynamicDelay<N> — 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 <array>
#include <cstddef>
#include <cstring>
#include "../core/perf.hpp"
namespace nisps {
// Fixed-length feedback delay. `play()` reads the head, mixes feedback into the
// slot, and advances. Behavior matches maxiDelayline<N>::play().
template <std::size_t N>
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<float, N> 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<N>`.
template <std::size_t N>
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<float>(N)) {}
void clear() noexcept {
std::memset(line_.data(), 0, N * sizeof(float));
write_index_ = 0u;
smoothed_size_ = static_cast<float>(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<float>(write_index_) - smoothed_size_;
if (read_pos < 0.f) read_pos += static_cast<float>(N);
const std::size_t i1 = static_cast<std::size_t>(read_pos);
const float frac = read_pos - static_cast<float>(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<float, N> line_{};
std::size_t write_index_ = 0u;
float smoothed_size_;
float smooth_coeff_ = 0.997f;
};
} // namespace nisps

118
nisps/dsp/env.hpp Normal file
View file

@ -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 <cmath>
#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

149
nisps/dsp/filter.hpp Normal file
View file

@ -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<NCh>
// 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 <array>
#include <cmath>
#include <cstddef>
#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 <std::size_t NCh>
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<float, NCh> 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

279
nisps/dsp/osc.hpp Normal file
View file

@ -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 <array>
#include <cmath>
#include <cstddef>
#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<float>(kTabRange);
static const float kTabRangeRcpr = 1.f / static_cast<float>(kTabRange);
static const float kTabScale = static_cast<float>(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<float>(static_cast<int>(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<int>(halfsine_scaled);
const float tabfrac = halfsine_scaled - static_cast<float>(table_index);
if (table_index < 0) table_index = 0;
if (table_index > static_cast<int>(kTabSize) - 2) table_index = static_cast<int>(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<float>(kTabRange) * static_cast<float>(kTabRange));
const float cauchy_slope = (-2.f * static_cast<float>(kTabRange)) * cauchy_val * cauchy_val;
const float addsq = -cauchy_slope / (2.f * static_cast<float>(kTabRange));
const float fake_at3 = cauchy_val + addsq * static_cast<float>(kTabRange) * static_cast<float>(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<float>(i) * (static_cast<float>(kTabRange) / static_cast<float>(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

102
nisps/dsp/pitch_shift.hpp Normal file
View file

@ -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 <array>
#include <cmath>
#include <cstddef>
#include <cstring>
#include "../core/perf.hpp"
namespace nisps {
template <std::size_t N = 8192>
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<float>(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<std::size_t>(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<float>(N))
- std::floor(rp / static_cast<float>(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<float>(N)) rp -= static_cast<float>(N);
if (rp < 0.f) rp += static_cast<float>(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<float, N> buffer_{};
std::size_t write_idx_ = 0u;
float read_pos_[2] = {0.f, static_cast<float>(N) * 0.5f};
};
} // namespace nisps

122
nisps/dsp/reverb.hpp Normal file
View file

@ -0,0 +1,122 @@
// nisps/dsp/reverb.hpp — Schroeder-Moorer reverb building blocks.
//
// AllPass<N> — feedback all-pass section
// Comb<N> — simple feedback comb
// LpComb<N> — 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 <array>
#include <cstddef>
#include <cstring>
#include "../core/perf.hpp"
namespace nisps {
template <std::size_t N>
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<float, N> line_{};
std::size_t idx_ = 0u;
};
template <std::size_t N>
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<float, N> 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 <std::size_t N>
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<float, N> line_{};
std::size_t idx_ = 0u;
float lp_y_ = 0.f;
};
} // namespace nisps

68
tests/cpp/test_dsp_biquad.cpp Normal file
View file

@ -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 <cmath>
#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));
}
}

50
tests/cpp/test_dsp_delay.cpp Normal file
View file

@ -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<float>(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<float>(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);
}
}

46
tests/cpp/test_dsp_pitch_shift.cpp Normal file
View file

@ -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 <cmath>
#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);
}

37
tests/cpp/test_dsp_reverb.cpp Normal file
View file

@ -0,0 +1,37 @@
// tests/cpp/test_dsp_reverb.cpp — bounded-output sanity for reverb sections.
#include <cmath>
#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<float>(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));
}
}