Port Cortex-M4 resonant filter to a pure vtable-free core/instrument/engine/filter module with 0.1-10 Q and log cutoff

CMakeLists.txt

@@ -995,6 +995,18 @@ target_link_libraries(curve_popup PUBLIC editor_geometry)
 add_library(master_gain STATIC src/core/instrument/engine/master_gain.cpp)
 target_include_directories(master_gain PUBLIC src)
 
+# filter — the per-voice 2-pole resonant low/high-pass, ported from Daniel's Cortex-M4 filter
+# with its virtual FilterBase/Filter/Biquad hierarchy flattened away (that hierarchy dispatched
+# virtually per channel per sample, which the per-voice per-sample path forbids). Control
+# mapping, RBJ coefficient math, feedback saturation, and the filter type each get their own
+# file; VoiceFilter::process is header-inline so the biquad kernel still inlines at the call
+# site. Standard library only. NEITHER SDK.
+add_library(filter STATIC
+    src/core/instrument/engine/filter/filter_params.cpp
+    src/core/instrument/engine/filter/filter_coeffs.cpp
+    src/core/instrument/engine/filter/voice_filter.cpp)
+target_include_directories(filter PUBLIC src)
+
 # sample_bands: the band-stack allocator's vertical inventory, asserted as pure geometry
 # (chrome / two-lane waveform / deck row) independent of any paint call — the contract the
 # band owners downstream read.
@@ -1095,6 +1107,13 @@ add_executable(master_gain_tests tests/test_master_gain.cpp)
 target_link_libraries(master_gain_tests PRIVATE master_gain)
 add_test(NAME master_gain_tests COMMAND master_gain_tests)
 
+# filter: the per-voice resonant filter. Pins the RBJ coefficients against an independent
+# textbook cos/sin derivation, asserts the cutoff/Q control mappings at their anchors, and
+# measures the resonant peak both analytically and by driving real sines. NEITHER SDK.
+add_executable(filter_tests tests/test_filter.cpp)
+target_link_libraries(filter_tests PRIVATE filter)
+add_test(NAME filter_tests COMMAND filter_tests)
+
 # ---------------------------------------------------------------------------
 # 4) The REAPER extension — a loadable module (dlopen'd by REAPER, not linked).
 # ---------------------------------------------------------------------------

src/core/instrument/engine/filter/CLAUDE.md

@@ -0,0 +1,92 @@
+# src/core/instrument/engine/filter — the per-voice resonant filter
+
+## Scope
+
+The pure 2-pole resonant low/high-pass a sounding voice runs. No REAPER, no VST3, no
+allocation, no I/O. Four files, one responsibility each:
+
+- `filter_params` — the control domain: `FilterMode`, normalized [0,1] knob position →
+  cutoff Hz and Q, and the exact inverses.
+- `filter_coeffs` — the DSP domain: `BiquadCoeffs` and the RBJ coefficient computation
+  from (mode, cutoff Hz, Q, sample rate).
+- `filter_saturate` — the high-pass feedback saturator (`tanhSaturate` /
+  `feedbackSaturate`). Header-only inline; it sits on the per-sample path.
+- `voice_filter` — `FilterSettings` and `VoiceFilter`, the concrete per-voice type.
+  `process()` is defined in the header.
+
+## Invariants
+
+### No vtable on the per-sample path
+
+This is a **port, not a relocation**. The Cortex-M4 source was a virtual hierarchy
+(`FilterBase` → `Filter` → `Biquad` → `{BiquadHP, BiquadLP}`) whose base class routed the
+channel loop through pure-virtual `process_channel_frame` / `filter` / `update_feedback`
+so a `FilterDecorator` chain could wrap it. **None of that came across, and none of it may
+come back.** `VoiceFilter` is concrete: mode is a member branch inside an inlined
+`process()`, predicted perfectly because it cannot change within a note. There is no
+`IFilter`, no decorator seam, no virtual `tick()`, and no allocation in `process()` — root
+`CLAUDE.md`'s structural heuristic 3 names this class of dispatch blowout directly.
+
+A non-type template parameter for the mode was considered and rejected: mode is a
+runtime-settable user parameter, so templating would only relocate the same branch to the
+call site and force the voice to hold two instances or switch over them.
+
+### Two modes, and only two
+
+2-pole high-pass and 2-pole low-pass. The source's `Biquad1PoleLP` is struck and was not
+ported. Further modes are deferred — **do not build a mode-extension framework** for them.
+
+### The cutoff control is sample-rate-free; the clamp is not
+
+`filterCutoffHzFromNorm` sweeps a fixed 20 Hz – 20 kHz (three exact decades, so norm 1/3
+is 200 Hz and 2/3 is 2 kHz) and takes no sample rate. The persisted value is the
+normalized knob position, so a rate-derived endpoint would make one preset sound different
+at 44.1k and 96k. The Nyquist clamp (`kFilterNyquistFraction`, 0.48) is a property of the
+bilinear transform — `tan(pi*fc/sr)` diverges at Nyquist — so it lives in `biquadCoeffs`
+where the rate is already a parameter. 20 kHz is under 0.48·sr at every supported rate, so
+the clamp never eats live knob travel; the source's hardcoded 23 kHz endpoint did exactly
+that at 44.1k.
+
+`biquadCoeffs` with a non-positive sample rate returns pass-through coefficients. It does
+**not** fall back to 44100 — that would breach the standing no-hardcoded-sample-rates
+ruling.
+
+### Q spans 0.1 → 10 with √2 at the center
+
+Settled by Daniel. The source's `Q = M_SQRT1_2 + resonance` mapping (floored at 0.707, no
+center anchor) was **rewritten, not ported**. The curve is quadratic in log Q through the
+three anchors rather than two spliced log segments — same anchors either way, but no slope
+kink at the center detent. The quadratic term is nonzero only because √2 is not the
+geometric mean of 0.1 and 10; `filterNormFromQ` divides by it.
+
+### The high-pass input feedback is load-bearing
+
+`kHighPassFeedbackShare` (0.24) times the raw **normalized** resonance, not Q — Q reaches
+10 and scaling the feedback by it would push loop gain past unity. The high-pass numerator
+collapses toward zero as cutoff falls, taking the resonance with it; the saturated
+feedback restores the character down there. Ported behavior; the constant is the tuning
+knob if the feel needs adjusting. `audio_saturate` and `H()` from the source were unused
+by the biquads and were not ported.
+
+### Denormal flushing
+
+`process()` flushes the **y** history to exact zero below `kFilterDenormalFloor` (1e-30).
+Only the recursive half needs it: a denormal in `y` self-sustains and stalls the FPU for
+thousands of samples on a ringing-out voice, while the `x` history is an FIR tail that
+shifts out within two samples. `isSilent()` reports the flushed state and is the honest
+signal that a voice's filter can no longer contribute output.
+
+## Gotchas
+
+- **The tan pre-warp is not a different filter.** By the half-angle identity
+  `cos(w0) = (1-w²)/(1+w²)` and `sin(w0) = 2w/(1+w²)` with `w = tan(pi*fc/sr)`, these are
+  the textbook RBJ cos/sin coefficients exactly — just computed in a form that stays
+  conditioned at low cutoff where `cos(w0) → 1`. `tests/test_filter.cpp` asserts the
+  equivalence against an independent derivation. Don't "simplify" it back to `std::cos`.
+- **`prepare()` deliberately does not clear history** — a live parameter move must glide,
+  not click. Call `reset()` at note-on.
+- **`a1`/`a2` are stored for a subtracting difference equation** (`y = ... - a1*y1 -
+  a2*y2`), so the transfer denominator is `1 + a1*z^-1 + a2*z^-2`. A sign convention slip
+  here inverts the poles.
+- **No call site yet.** Wiring the filter into the voice path is a separate track; nothing
+  in `sampler_core` references this module today.

src/core/instrument/engine/filter/filter_coeffs.cpp

@@ -0,0 +1,43 @@
+#include "core/instrument/engine/filter/filter_coeffs.h"
+
+#include <cmath>
+
+namespace reasampler::instrument::engine {
+namespace {
+
+// M_PI is not standard C++ and is absent on MSVC without _USE_MATH_DEFINES.
+constexpr double kPi = 3.14159265358979323846;
+
+double clampd(double v, double lo, double hi) { return v < lo ? lo : (v > hi ? hi : v); }
+
+}  // namespace
+
+BiquadCoeffs biquadCoeffs(FilterMode mode, float cutoffHz, float q, double sampleRate) {
+    if (!(sampleRate > 0.0)) return BiquadCoeffs{};
+
+    const double nyquistCeiling = kFilterNyquistFraction * sampleRate;
+    const double fc = clampd(cutoffHz, kFilterCutoffMinHz, nyquistCeiling);
+    const double qq = clampd(q, kFilterQMin, kFilterQMax);
+
+    const double w = std::tan(kPi * fc / sampleRate);
+    const double w2 = w * w;
+    const double cosw = (1.0 - w2) / (1.0 + w2);
+    const double sinw = 2.0 * w / (1.0 + w2);
+    const double alpha = sinw / (2.0 * qq);
+    const double norm = 1.0 / (1.0 + alpha);
+
+    // Both modes share the denominator; only the numerator's sign on cosw differs, and b1 is
+    // always +/-2*b0 — folding that in keeps the two branches from drifting apart.
+    const double b0 = (mode == FilterMode::HighPass ? (1.0 + cosw) : (1.0 - cosw)) * 0.5 * norm;
+    const double b1 = (mode == FilterMode::HighPass ? -2.0 : 2.0) * b0;
+
+    BiquadCoeffs c;
+    c.b0 = static_cast<float>(b0);
+    c.b1 = static_cast<float>(b1);
+    c.b2 = static_cast<float>(b0);
+    c.a1 = static_cast<float>(-2.0 * cosw * norm);
+    c.a2 = static_cast<float>((1.0 - alpha) * norm);
+    return c;
+}
+
+}  // namespace reasampler::instrument::engine

src/core/instrument/engine/filter/filter_coeffs.h

@@ -0,0 +1,31 @@
+// filter_coeffs.h — RBJ Audio EQ Cookbook Direct Form I biquad coefficients for the 2-pole
+// low/high-pass. Computed via the tan half-angle substitution w = tan(pi*fc/sr): by the
+// identity cos(w0) = (1-w^2)/(1+w^2), sin(w0) = 2w/(1+w^2) these ARE the textbook cos/sin
+// coefficients, in a form that stays conditioned at low cutoff where cos(w0) -> 1.
+
+#pragma once
+
+#include "core/instrument/engine/filter/filter_params.h"
+
+namespace reasampler::instrument::engine {
+
+// Already normalized by a0. The denominator is 1 + a1*z^-1 + a2*z^-2, so the difference
+// equation SUBTRACTS the a terms: y = b0*x + b1*x1 + b2*x2 - a1*y1 - a2*y2.
+struct BiquadCoeffs {
+    float b0 = 1.0f;
+    float b1 = 0.0f;
+    float b2 = 0.0f;
+    float a1 = 0.0f;
+    float a2 = 0.0f;
+};
+
+// Highest fraction of the sample rate the pre-warp stays well-conditioned at: tan() diverges
+// as fc approaches sr/2. Ported unchanged from the firmware, where it was already the ceiling.
+inline constexpr double kFilterNyquistFraction = 0.48;
+
+// cutoffHz is clamped into [kFilterCutoffMinHz, kFilterNyquistFraction*sampleRate] and q into
+// [kFilterQMin, kFilterQMax]. A non-positive sampleRate yields pass-through coefficients — the
+// no-hardcoded-sample-rates ruling means we refuse to invent a rate rather than assume 44.1k.
+BiquadCoeffs biquadCoeffs(FilterMode mode, float cutoffHz, float q, double sampleRate);
+
+}  // namespace reasampler::instrument::engine

src/core/instrument/engine/filter/filter_params.cpp

@@ -0,0 +1,59 @@
+#include "core/instrument/engine/filter/filter_params.h"
+
+#include <cmath>
+
+namespace reasampler::instrument::engine {
+namespace {
+
+double clamp01(double v) { return v < 0.0 ? 0.0 : (v > 1.0 ? 1.0 : v); }
+
+// log Q = A + B*n + C*n^2, solved from the three anchor points. C is nonzero precisely
+// because the center anchor sqrt(2) is not the geometric mean of the endpoints (which is 1);
+// were they equal the curve would degenerate to a plain log sweep and the inverse below
+// would divide by zero.
+struct QCurve {
+    double a, b, c;
+};
+
+QCurve qCurve() {
+    const double lo = std::log(static_cast<double>(kFilterQMin));
+    const double mid = std::log(static_cast<double>(kFilterQCenter));
+    const double hi = std::log(static_cast<double>(kFilterQMax));
+    return {lo, 4.0 * mid - 3.0 * lo - hi, 2.0 * lo + 2.0 * hi - 4.0 * mid};
+}
+
+}  // namespace
+
+float filterCutoffHzFromNorm(float norm) {
+    const double lo = std::log(static_cast<double>(kFilterCutoffMinHz));
+    const double hi = std::log(static_cast<double>(kFilterCutoffMaxHz));
+    return static_cast<float>(std::exp(lo + clamp01(norm) * (hi - lo)));
+}
+
+float filterNormFromCutoffHz(float hz) {
+    if (!(hz > 0.0f)) return 0.0f;
+    const double lo = std::log(static_cast<double>(kFilterCutoffMinHz));
+    const double hi = std::log(static_cast<double>(kFilterCutoffMaxHz));
+    return static_cast<float>(clamp01((std::log(static_cast<double>(hz)) - lo) / (hi - lo)));
+}
+
+float filterQFromNorm(float norm) {
+    const QCurve k = qCurve();
+    const double n = clamp01(norm);
+    return static_cast<float>(std::exp(k.a + n * (k.b + k.c * n)));
+}
+
+float filterNormFromQ(float q) {
+    if (!(q > kFilterQMin)) return 0.0f;
+    if (q >= kFilterQMax) return 1.0f;
+    // Clamping first is load-bearing, not just tidy: the parabola peaks at log Q well below
+    // an arbitrarily large q, so an unclamped out-of-range value has no real root at all.
+    const QCurve k = qCurve();
+    const double d = k.b * k.b - 4.0 * k.c * (k.a - std::log(static_cast<double>(q)));
+    if (!(d >= 0.0)) return 0.0f;
+    // Of the two roots only this one lies on the rising branch inside [0,1]; the parabola's
+    // vertex sits well above 1 for the settled anchors.
+    return static_cast<float>(clamp01((-k.b + std::sqrt(d)) / (2.0 * k.c)));
+}
+
+}  // namespace reasampler::instrument::engine

src/core/instrument/engine/filter/filter_params.h

@@ -0,0 +1,39 @@
+// filter_params.h — control-domain mapping for the voice filter: normalized [0,1] knob
+// positions to cutoff Hz and Q, plus the two-mode enum. Deliberately sample-rate-free —
+// the Nyquist clamp is a property of the bilinear transform and lives in filter_coeffs,
+// so the persisted normalized cutoff means the same frequency at every project rate.
+
+#pragma once
+
+namespace reasampler::instrument::engine {
+
+enum class FilterMode { LowPass, HighPass };
+
+// The audio band the cutoff control sweeps: three exact decades, so norm 1/3 is 200 Hz and
+// norm 2/3 is 2 kHz. NOT derived from the sample rate — a rate-dependent endpoint would make
+// one saved preset sound different at 44.1k and 96k, and at 44.1k the top of the travel would
+// be dead against the Nyquist clamp (the ported firmware's 23 kHz endpoint had exactly that
+// defect). 20 kHz sits under 0.48*sr at every rate we support, so the whole knob stays live.
+inline constexpr float kFilterCutoffMinHz = 20.0f;
+inline constexpr float kFilterCutoffMaxHz = 20000.0f;
+
+// Q spans the full range with Butterworth (sqrt(2)) at the control's center detent.
+inline constexpr float kFilterQMin = 0.1f;
+inline constexpr float kFilterQMax = 10.0f;
+inline constexpr float kFilterQCenter = 1.41421356f;
+
+// Out-of-range norm clamps to the endpoints.
+float filterCutoffHzFromNorm(float norm);
+
+// Exact inverse of filterCutoffHzFromNorm over the band; out-of-band Hz clamps to 0 or 1.
+float filterNormFromCutoffHz(float hz);
+
+// A single smooth curve — quadratic in log Q — through (0, kFilterQMin),
+// (0.5, kFilterQCenter), (1, kFilterQMax), rather than two spliced log segments. Same three
+// anchors either way, but the single curve has no slope kink at the center detent.
+float filterQFromNorm(float norm);
+
+// Exact inverse of filterQFromNorm; out-of-range Q clamps to 0 or 1.
+float filterNormFromQ(float q);
+
+}  // namespace reasampler::instrument::engine

src/core/instrument/engine/filter/filter_saturate.h

@@ -0,0 +1,27 @@
+// filter_saturate.h — the high-pass feedback-path saturator, ported from the Cortex-M4
+// filter. Header-inline: it sits on the per-voice per-sample path, and a rational
+// approximation is here precisely to avoid a transcendental tanh() call there.
+
+#pragma once
+
+namespace reasampler::instrument::engine {
+
+// Rational tanh approximation inside +/-threshold, continued past it with a gentle 0.1 slope
+// anchored at the threshold value so the curve stays continuous rather than hard-clipping.
+inline float tanhSaturate(float x, float threshold, float a, float b) {
+    if (x > threshold) {
+        const float satAtThreshold = threshold * a / (a + b + threshold * threshold);
+        return satAtThreshold + (x - threshold) * 0.1f;
+    }
+    if (x < -threshold) {
+        const float satAtThreshold = -threshold * a / (a + b + threshold * threshold);
+        return satAtThreshold + (x + threshold) * 0.1f;
+    }
+    return x * a / (a + b + x * x);
+}
+
+// TB-303-style hard feedback saturation. Tuned for the large excursions a resonant feedback
+// path produces, not for audio-level signals — do not reuse it as a general waveshaper.
+inline float feedbackSaturate(float x) { return tanhSaturate(x, 2.0f, 27.0f, 9.0f); }
+
+}  // namespace reasampler::instrument::engine

src/core/instrument/engine/filter/voice_filter.cpp

@@ -0,0 +1,29 @@
+#include "core/instrument/engine/filter/voice_filter.h"
+
+namespace reasampler::instrument::engine {
+
+void VoiceFilter::prepare(const FilterSettings& settings, double sampleRate) {
+    mode_ = settings.mode;
+    const float cutoffHz = filterCutoffHzFromNorm(settings.cutoffNorm);
+    coeffs_ = biquadCoeffs(settings.mode, cutoffHz, filterQFromNorm(settings.resonanceNorm),
+                           sampleRate);
+    const float res = settings.resonanceNorm < 0.0f
+                          ? 0.0f
+                          : (settings.resonanceNorm > 1.0f ? 1.0f : settings.resonanceNorm);
+    fbAmount_ = res * kHighPassFeedbackShare;
+}
+
+void VoiceFilter::reset() {
+    for (State& s : state_) s = State{};
+}
+
+bool VoiceFilter::isSilent() const {
+    for (const State& s : state_) {
+        if (s.x1 != 0.0f || s.x2 != 0.0f || s.y1 != 0.0f || s.y2 != 0.0f || s.fb != 0.0f) {
+            return false;
+        }
+    }
+    return true;
+}
+
+}  // namespace reasampler::instrument::engine

src/core/instrument/engine/filter/voice_filter.h

@@ -0,0 +1,113 @@
+// voice_filter.h — per-voice 2-pole resonant low/high-pass. Concrete type, no vtable: this
+// sits on the per-voice per-sample path, so process() is header-inline and mode is a member
+// branch. No allocation, no virtual dispatch, no I/O anywhere in process().
+
+#pragma once
+
+#include <cassert>
+#include <type_traits>
+
+#include "core/instrument/engine/filter/filter_coeffs.h"
+#include "core/instrument/engine/filter/filter_params.h"
+#include "core/instrument/engine/filter/filter_saturate.h"
+
+namespace reasampler::instrument::engine {
+
+// Normalized control positions, as the editor moves them and the persisted state carries them.
+struct FilterSettings {
+    FilterMode mode = FilterMode::LowPass;
+    float cutoffNorm = 1.0f;
+    float resonanceNorm = 0.0f;
+};
+
+// Share of the last output fed back into the high-pass input at full resonance. Driven by the
+// raw control position rather than by Q: Q reaches 10, and scaling the feedback by it would
+// push the loop gain past unity at the top of the range.
+inline constexpr float kHighPassFeedbackShare = 0.24f;
+
+// Below this the recursion has decayed past -600 dB. Flushing keeps the history out of the
+// subnormal range, where a ringing-out voice would otherwise stall the FPU for thousands of
+// samples. Chosen well above FLT_MIN so a flushed state can never re-enter that range.
+inline constexpr float kFilterDenormalFloor = 1e-30f;
+
+class VoiceFilter {
+public:
+    // The instrument's output bus is permanently stereo; one history line per channel.
+    static constexpr int kMaxChannels = 2;
+
+    struct State {
+        float x1 = 0.0f;
+        float x2 = 0.0f;
+        float y1 = 0.0f;
+        float y2 = 0.0f;
+        float fb = 0.0f;  // last output; the high-pass input-feedback tap
+    };
+
+    // Recomputes coefficients from the control positions. History is deliberately preserved so
+    // a live parameter move glides instead of clicking; call reset() at note-on.
+    void prepare(const FilterSettings& settings, double sampleRate);
+
+    void reset();
+
+    // Hot path. `channel` must be in [0, kMaxChannels).
+    float process(int channel, float x) {
+        assert(channel >= 0 && channel < kMaxChannels);
+        State& s = state_[channel];
+
+        // The high-pass numerator collapses toward zero as cutoff falls, taking the resonance
+        // with it; feeding a saturated share of the last output back into the input restores
+        // the character the coefficients alone stop producing down there.
+        const float in = (mode_ == FilterMode::HighPass)
+                             ? x - fbAmount_ * feedbackSaturate(s.fb * 0.9f)
+                             : x;
+
+        const float y = coeffs_.b0 * in + coeffs_.b1 * s.x1 + coeffs_.b2 * s.x2
+                        - coeffs_.a1 * s.y1 - coeffs_.a2 * s.y2;
+
+        s.x2 = s.x1;
+        s.x1 = in;
+        s.y2 = s.y1;
+        s.y1 = y;
+
+        // Snap the WHOLE state once the recursion as a whole has decayed past -600 dB.
+        // Zeroing individual samples instead does not work: a resonator swings through zero
+        // twice a cycle, so a per-sample flush injects a step in phase with the resonance,
+        // which the resonance then amplifies — the filter limit-cycles at the floor forever
+        // rather than going quiet. Testing y1 AND y2 tests the envelope, not one sample.
+        if (s.y1 > -kFilterDenormalFloor && s.y1 < kFilterDenormalFloor &&
+            s.y2 > -kFilterDenormalFloor && s.y2 < kFilterDenormalFloor) {
+            s = State{};
+        }
+        s.fb = s.y1;
+        return y;
+    }
+
+    void processFrame(float* samples, int channelCount) {
+        assert(channelCount >= 0 && channelCount <= kMaxChannels);
+        for (int c = 0; c < channelCount; ++c) samples[c] = process(c, samples[c]);
+    }
+
+    // True once every history line has flushed to exact zero — the voice's filter has stopped
+    // ringing and cannot contribute further output.
+    bool isSilent() const;
+
+    const State& state(int channel) const {
+        assert(channel >= 0 && channel < kMaxChannels);
+        return state_[channel];
+    }
+    const BiquadCoeffs& coeffs() const { return coeffs_; }
+
+private:
+    BiquadCoeffs coeffs_{};
+    FilterMode mode_ = FilterMode::LowPass;
+    float fbAmount_ = 0.0f;
+    State state_[kMaxChannels]{};
+};
+
+// The port's whole point, enforced by the compiler rather than by review: the source was a
+// virtual hierarchy dispatching per channel per sample, and this type must never grow one
+// back. Trivially copyable also means nothing here is heap-owned.
+static_assert(!std::is_polymorphic_v<VoiceFilter>, "no vtable on the per-sample path");
+static_assert(std::is_trivially_copyable_v<VoiceFilter>, "state is plain values, never owned");
+
+}  // namespace reasampler::instrument::engine

src/temp-cortex/biquad.hpp

@@ -1,73 +0,0 @@
-#pragma once
-
-#include "util.hpp"
-#include "filter.hpp"
-#include "filter_params.hpp"
-
-template <int k_channels, typename TUIParams>
-class Biqaud : public Filter<k_channels, FeedbackLine, NormalCoefficients, TUIParams, FilterParameters>
-{
-    public:
-        Biqaud(const uint32_t& sample_rate, FilterParameters *params);
-
-        void prepare_parameters(const TUIParams& params) override;
-        
-    protected:
-        uint32_t sample_rate;
-
-        void process_channel_frame(FeedbackLine& state, 
-                                   const NormalCoefficients& coeff, 
-                                   const float& x, 
-                                   float& y) override;
-
-        void filter(FeedbackLine& state, 
-                    const NormalCoefficients& coeff, 
-                    const float& x, 
-                    float& y) override;
-
-        void update_feedback(FeedbackLine& state, 
-                             const NormalCoefficients& coeff, 
-                             const float& x, 
-                             float& y) override;
-};
-
-template <int k_channels, typename TUIParams>
-class BiquadHP : public Biqaud<k_channels, TUIParams>
-{
-    public:
-        BiquadHP(const uint32_t& sample_rate, FilterParameters *params);
-
-        NormalCoefficients prepare_coefficients() override;
-        
-    protected:
-        void process_channel_frame(FeedbackLine& state, 
-                                   const NormalCoefficients& coeff, 
-                                   const float& x, 
-                                   float& y) override;
-};
-
-template <int k_channels, typename TUIParams>
-class BiquadLP : public Biqaud<k_channels, TUIParams>
-{
-    public:
-        BiquadLP(const uint32_t& sample_rate, FilterParameters *params);
-
-        NormalCoefficients prepare_coefficients() override;
-};
-
-template <int k_channels, typename TUIParams>
-class Biquad1PoleLP : public BiquadLP<k_channels, TUIParams>
-{
-    public:
-        Biquad1PoleLP(const uint32_t& sample_rate, FilterParameters *params);
-
-        NormalCoefficients prepare_coefficients() override;
-        
-    protected:
-        void process_channel_frame(FeedbackLine& state, 
-                                   const NormalCoefficients& coeff, 
-                                   const float& x, 
-                                   float& y) override;
-};
-
-#include "biquad.tpp"

src/temp-cortex/biquad.tpp

@@ -1,195 +0,0 @@
-#pragma once
-
-#include "basicmaths.h"
-#include "biquad.hpp"
-
-template <int k_channels, typename TUIParams>
-Biqaud<k_channels, TUIParams>::Biqaud(const uint32_t& p_sample_rate, FilterParameters *p_params)
-: Filter<k_channels, FeedbackLine, NormalCoefficients, TUIParams, FilterParameters>(p_params), sample_rate(p_sample_rate) 
-{ 
-    // Initialize filter state to zero to prevent random behavior
-    for (int i = 0; i < k_channels; i++) {
-        this->state[i].x[0] = 0.0f;
-        this->state[i].x[1] = 0.0f;
-        this->state[i].y[0] = 0.0f;
-        this->state[i].y[1] = 0.0f;
-        this->state[i].fb = 0.0f;
-    }
-}
-
-template <int k_channels, typename TUIParams>
-void Biqaud<k_channels, TUIParams>::prepare_parameters(const TUIParams& params)
-{
-    // Direct logarithmic interpolation for smooth frequency scaling using standard math
-    const float min_freq = 10.f;
-    const float max_freq = 23000.f;  
-    float log_freq = logf(min_freq) + params.p_cutoff * (logf(max_freq) - logf(min_freq));
-    float raw_cutoff = expf(log_freq);
-    this->params->cutoff = fminf(raw_cutoff, 0.48f * this->sample_rate);  // Allow closer to Nyquist
-
-    // Resonance response
-    this->params->res = params.p_resonance;
-
-    // Base Q of 0.707 plus resonance
-    this->params->Q = M_SQRT1_2 + this->params->res;
-}
-
-template <int k_channels, typename TUIParams>
-void Biqaud<k_channels, TUIParams>::process_channel_frame(FeedbackLine& state,
-                                                               const NormalCoefficients& coeff,
-                                                               const float& x, 
-                                                               float& y)
-{
-    this->filter(state, coeff, x, y);
-    this->update_feedback(state, coeff, x, y);
-}
-
-template <int k_channels, typename TUIParams>
-void Biqaud<k_channels, TUIParams>::filter(FeedbackLine &state, const NormalCoefficients &coeff, const float &x, float &y)
-{
-    // debugMessage("Biqaud::filter");
-    // Direct Form I biquad - matches Audio EQ Cookbook exactly
-    y = coeff.b0 * x + coeff.b1 * state.x[0] + coeff.b2 * state.x[1] 
-        - coeff.a1 * state.y[0] - coeff.a2 * state.y[1];
-}
-
-template <int k_channels, typename TUIParams>
-void Biqaud<k_channels, TUIParams>::update_feedback(FeedbackLine& state, const NormalCoefficients& coeff, const float& x, float& y)
-{
-    // debugMessage("State x[0], x[1], y[0]: ", state.x[0], state.x[1], state.y[0]);
-    // Update feedback state
-    state.x[1] = state.x[0];
-    state.x[0] = x;
-    state.y[1] = state.y[0];
-    state.y[0] = y;
-    state.fb = y;
-}
-
-template <int k_channels, typename TUIParams>
-BiquadHP<k_channels, TUIParams>::BiquadHP(const uint32_t& p_sample_rate, FilterParameters *p_params)
-: Biqaud<k_channels, TUIParams>(p_sample_rate, p_params) {}
-
-template <int k_channels, typename TUIParams>
-void BiquadHP<k_channels, TUIParams>::process_channel_frame(FeedbackLine& state,
-                                                                 const NormalCoefficients& coeff,
-                                                                 const float& x, 
-                                                                 float& y)
-{
-    // CRITICAL: Highpass filters require input feedback to work properly
-    // This compensates for coefficient collapse at low frequencies
-    const float fb_amount = this->params->res * 0.24f;
-    float input = x - fb_amount * feedback_saturate(state.fb * 0.9f);
-    
-    Biqaud<k_channels, TUIParams>::process_channel_frame(state, coeff, input, y);
-}
-
-template <int k_channels, typename TUIParams>
-NormalCoefficients BiquadHP<k_channels, TUIParams>::prepare_coefficients()
-{
-    // Pre-warped bilinear transform - same topology as lowpass but for highpass
-    const float w = tanf(M_PI * this->params->cutoff / this->sample_rate);
-    const float w2 = w * w;
-    const float cosw = (1.0f - w2) / (1.0f + w2);
-    const float sinw = 2.0f * w / (1.0f + w2);
-    const float alpha = sinw / (2.0f * this->params->Q);
-    
-    // Standard RBJ highpass with pre-warped frequency
-    const float norm = 1.0f / (1.0f + alpha);
-    const float b0 = (1.0f + cosw) * 0.5f * norm;
-    const float b1 = -(1.0f + cosw) * norm;
-    const float b2 = (1.0f + cosw) * 0.5f * norm;
-    const float a1 = -2.0f * cosw * norm;
-    const float a2 = (1.0f - alpha) * norm;
-    
-    NormalCoefficients coeff = {
-        .a1 = a1,
-        .a2 = a2,
-        .b0 = b0,
-        .b1 = b1,
-        .b2 = b2
-    };
-
-    return coeff;
-}
-
-template <int k_channels, typename TUIParams>
-BiquadLP<k_channels, TUIParams>::BiquadLP(const uint32_t& p_sample_rate, FilterParameters *p_params)
-: Biqaud<k_channels, TUIParams>(p_sample_rate, p_params) {}
-
-template <int k_channels, typename TUIParams>
-NormalCoefficients BiquadLP<k_channels, TUIParams>::prepare_coefficients()
-{
-    // Pre-warped bilinear transform - correct implementation
-    const float w = tanf(M_PI * this->params->cutoff / this->sample_rate);
-    const float w2 = w * w;
-    const float cosw = (1.0f - w2) / (1.0f + w2);
-    const float sinw = 2.0f * w / (1.0f + w2);
-    const float alpha = sinw / (2.0f * this->params->Q);
-    
-    // Standard RBJ lowpass with pre-warped frequency
-    const float norm = 1.0f / (1.0f + alpha);
-    const float b0 = (1.0f - cosw) * 0.5f * norm;
-    const float b1 = (1.0f - cosw) * norm;
-    const float b2 = (1.0f - cosw) * 0.5f * norm;
-    const float a1 = -2.0f * cosw * norm;
-    const float a2 = (1.0f - alpha) * norm;
-    
-    NormalCoefficients coeff = {
-        .a1 = a1,
-        .a2 = a2,
-        .b0 = b0,
-        .b1 = b1,
-        .b2 = b2
-    };
-
-    return coeff;
-}
-
-template <int k_channels, typename TUIParams>
-Biquad1PoleLP<k_channels, TUIParams>::Biquad1PoleLP(const uint32_t& p_sample_rate, FilterParameters *p_params)
-: BiquadLP<k_channels, TUIParams>(p_sample_rate, p_params) {}
-
-template <int k_channels, typename TUIParams>
-void Biquad1PoleLP<k_channels, TUIParams>::process_channel_frame(FeedbackLine& state,
-                                                                      const NormalCoefficients& coeff,
-                                                                      const float& x, 
-                                                                      float& y)
-{
-    // Stable Moog-style feedback with conservative limits
-    // Much more conservative k values for single-pole stability
-    const float k_max = 3.8f;  // Much lower max for stability
-    const float k = fminf(k_max, fmaxf(0.0f, (this->params->Q - M_SQRT1_2)));  // Conservative Q mapping
-    
-    // Conservative gain compensation
-    const float makeup_gain = 1.0f + k * 0.5f;  // Gentler compensation
-    
-    // Stable global feedback with limiting
-    float resonant_input = (x - k * feedback_saturate(state.fb * 0.8f)) * makeup_gain;
-    
-    // Process with stable resonant input
-    Biqaud<k_channels, TUIParams>::process_channel_frame(state, coeff, resonant_input, y);
-}
-
-template <int k_channels, typename TUIParams>
-NormalCoefficients Biquad1PoleLP<k_channels, TUIParams>::prepare_coefficients()
-{
-    // Correct 1-pole lowpass using bilinear transform
-    // H(s) = wc/(s + wc) -> H(z) = b0*(1+z^-1)/(1 + a1*z^-1)
-    const float w = tanf(M_PI * this->params->cutoff / this->sample_rate);
-    
-    // Bilinear transform gives both b0 and b1 coefficients
-    const float norm = 1.0f / (1.0f + w);
-    const float b0 = w * norm;     // Coefficient for x[n]
-    const float b1 = w * norm;     // Coefficient for x[n-1] (same as b0)
-    const float a1 = (w - 1.0f) * norm;  // Pole coefficient
-    
-    NormalCoefficients coeff = {
-        .a1 = a1,     // Pole coefficient  
-        .a2 = 0.0f,   // 1-pole has no second pole
-        .b0 = b0,     // Current input coefficient
-        .b1 = b1,     // Previous input coefficient  
-        .b2 = 0.0f    // 1-pole has no z^-2 numerator
-    };
-    
-    return coeff;
-}

src/temp-cortex/filter.hpp

@@ -1,62 +0,0 @@
-#pragma once
-
-#include "util.hpp"
-// #include "fdecorator.hpp"
-
-template <int k_channels, typename TFeedbackLine, typename TCoefficients, typename TUIParams, typename TFilterParams>
-class FilterBase
-{
-    public:        
-        FilterBase(TFilterParams *p) : params(p) {}
-        
-        /// @brief Prepare the filter channels to process all frames in this block
-        virtual void prepare_parameters(const TUIParams& params) = 0;
-
-        /// @brief Prepare the filter channels to process all frames in this block
-        virtual TCoefficients prepare_coefficients() = 0;
-                
-        /// @brief process the current frame samples for all channels
-        /// @param x inputs samples
-        /// @param y output samples
-        virtual void process_frame(const TCoefficients& coeff, const float x[k_channels], float y[k_channels])
-        {
-            // Handle channel iteration in the base class to ensure virtual dispatch through decorator chain
-            for (uint16_t channel = 0; channel < k_channels; channel++) {
-                this->process_channel_frame(this->state[channel], coeff, x[channel], y[channel]);
-            }
-        }
-
-    protected: 
-        /// @brief process the current frame sample for given channel
-        /// @param x inputs sample
-        /// @param y output sample
-        virtual void process_channel_frame(TFeedbackLine& state, const TCoefficients& coeff, const float& x, float& y) = 0;
-
-        /// @brief filter the current frame sample for given channel
-        /// @param state filter state
-        /// @param coeff filter coefficients
-        /// @param x input sample
-        /// @param y output sample
-        virtual void filter(TFeedbackLine& state, const TCoefficients& coeff, const float& x, float& y) = 0;
-
-        /// @brief update the feedback line for the next frame
-        /// @param state filter state
-        /// @param coeff filter coefficients
-        /// @param x input sample
-        /// @param y output sample
-        virtual void update_feedback(TFeedbackLine& state, const TCoefficients& coeff, const float& x, float& y) = 0;
-
-        TFilterParams* params;
-        TFeedbackLine state[k_channels];
-
-        template <int, typename, typename, typename, typename, typename, typename> 
-        friend class FilterDecorator;
-};
-
-template <int k_channels, typename TFeedbackLine, typename TCoefficients, typename TUIParams, typename TFilterParams>
-class Filter : public FilterBase<k_channels, TFeedbackLine, TCoefficients, TUIParams, TFilterParams>
-{
-    public:  
-        Filter(TFilterParams *p) 
-        : FilterBase<k_channels, TFeedbackLine, TCoefficients, TUIParams, TFilterParams>(p) {}
-};

src/temp-cortex/filter_params.hpp

@@ -1,22 +0,0 @@
-#pragma once
-
-typedef struct
-{
-    float cutoff;
-    float res;
-    float Q;
-} FilterParameters;
-
-typedef struct {
-    float a1;
-    float a2;
-    float b0;
-    float b1;
-    float b2;
-} NormalCoefficients;
-
-typedef struct {
-    float x[2];  // Previous inputs
-    float y[2];  // Previous outputs
-    float fb;  // Feedback value for resonance
-} FeedbackLine;

src/temp-cortex/util.hpp

@@ -1,46 +0,0 @@
-#pragma once
-#include "basicmaths.h"
-
-// Improved tanh approximation with proper continuity
-static inline float tanh_saturate(float x, float threshold, float a, float b)
-{
-    if (x > threshold) {
-        float excess = x - threshold;
-        float sat_val = threshold * a / (a + b + threshold * threshold);  // Value at threshold
-        return sat_val + excess * 0.1f;  // Gentle slope beyond threshold
-    }
-    if (x < -threshold) {
-        float excess = x + threshold; 
-        float sat_val = -threshold * a / (a + b + threshold * threshold);  // Value at -threshold
-        return sat_val + excess * 0.1f;  // Gentle slope beyond -threshold
-    }
-    const float x2 = x * x;
-    return x * a / (a + b + x2);
-}
-
-// TB-303 style feedback saturation
-// Hard saturation for filter feedback (handles large values)
-static inline float feedback_saturate(float x)
-{
-    // More aggressive saturation for feedback control
-    return tanh_saturate(x, 2.0f, 27.f, 9.f);
-}
-
-// Gentle saturation for audio signals (subtle, musical)
-static inline float audio_saturate(float x)
-{
-    // Adjusted parameters to maintain more volume at threshold
-    // At x=0.92: output ≈ 0.85 (much better than previous 0.57)
-    return tanh_saturate(x, 0.92f, 15.0f, 1.0f);
-}
-
-/// @brief Tunable logistic function (sigmoid)
-/// @param a slope
-/// @param b slope 2
-/// @param c offset
-/// @param z portion scalar
-/// @return H(x)
-static inline float H(float x, float a, float b, float c, float z)
-{
-    return z * a / (a + expf(b * (c - x))) - 0.02f;
-}

tests/test_filter.cpp

@@ -0,0 +1,450 @@
+// Standalone tests for the per-voice filter — no VST3, no REAPER, no framework. Same fast
+// assert loop as the sibling pure tests. The coefficient pins are literals so a refactor that
+// changes the DSP fails loudly; they are cross-checked in-test against a textbook RBJ
+// derivation (std::cos/std::sin) that shares no code with the implementation.
+
+#include "../src/core/instrument/engine/filter/filter_coeffs.h"
+#include "../src/core/instrument/engine/filter/filter_params.h"
+#include "../src/core/instrument/engine/filter/filter_saturate.h"
+#include "../src/core/instrument/engine/filter/voice_filter.h"
+
+#include <cfloat>
+#include <cmath>
+#include <cstdio>
+#include <initializer_list>
+
+using namespace reasampler::instrument::engine;
+
+static int g_fail = 0;
+#define CHECK(cond) do { if(!(cond)) { \
+    std::printf("FAIL line %d: %s\n", __LINE__, #cond); ++g_fail; } } while(0)
+#define CHECK_NEAR(a, b, eps) do { const double a_ = (a), b_ = (b); \
+    if (!(std::fabs(a_ - b_) <= (eps))) { \
+        std::printf("FAIL line %d: %s (%.10f) != %s (%.10f), delta %.3e\n", \
+                    __LINE__, #a, a_, #b, b_, std::fabs(a_ - b_)); ++g_fail; } } while(0)
+
+static constexpr double kPi = 3.14159265358979323846;
+
+// ---------------------------------------------------------------------------
+// Cutoff mapping
+// ---------------------------------------------------------------------------
+
+static void testCutoffMapsThreeDecadesLogarithmically() {
+    CHECK_NEAR(filterCutoffHzFromNorm(0.0f), 20.0, 1e-3);
+    CHECK_NEAR(filterCutoffHzFromNorm(1.0f), 20000.0, 1e-2);
+
+    // Exactly three decades, so the decade midpoints land on round numbers.
+    CHECK_NEAR(filterCutoffHzFromNorm(1.0f / 3.0f), 200.0, 1e-3);
+    CHECK_NEAR(filterCutoffHzFromNorm(2.0f / 3.0f), 2000.0, 1e-2);
+
+    // Half-decade steps confirm the sweep is log, not linear.
+    CHECK_NEAR(filterCutoffHzFromNorm(1.0f / 6.0f), 20.0 * std::sqrt(10.0), 1e-3);
+    CHECK_NEAR(filterCutoffHzFromNorm(0.5f), 20.0 * std::sqrt(1000.0), 1e-2);
+
+    // A linear sweep would put the midpoint at 10010 Hz; the log sweep is nowhere near it.
+    CHECK(filterCutoffHzFromNorm(0.5f) < 1000.0f);
+
+    CHECK_NEAR(filterCutoffHzFromNorm(-1.0f), 20.0, 1e-3);
+    CHECK_NEAR(filterCutoffHzFromNorm(2.0f), 20000.0, 1e-2);
+}
+
+static void testCutoffNormRoundTrips() {
+    for (int i = 0; i <= 20; ++i) {
+        const float n = static_cast<float>(i) / 20.0f;
+        CHECK_NEAR(filterNormFromCutoffHz(filterCutoffHzFromNorm(n)), n, 1e-6);
+    }
+    CHECK_NEAR(filterNormFromCutoffHz(200.0f), 1.0 / 3.0, 1e-6);
+    CHECK_NEAR(filterNormFromCutoffHz(2000.0f), 2.0 / 3.0, 1e-6);
+    CHECK(filterNormFromCutoffHz(1.0f) == 0.0f);
+    CHECK(filterNormFromCutoffHz(0.0f) == 0.0f);
+    CHECK(filterNormFromCutoffHz(48000.0f) == 1.0f);
+}
+
+// ---------------------------------------------------------------------------
+// Q mapping
+// ---------------------------------------------------------------------------
+
+static void testQSpansPointOneToTenWithRootTwoAtCenter() {
+    CHECK_NEAR(filterQFromNorm(0.0f), 0.1, 1e-6);
+    CHECK_NEAR(filterQFromNorm(0.5f), std::sqrt(2.0), 1e-5);
+    CHECK_NEAR(filterQFromNorm(1.0f), 10.0, 1e-4);
+
+    CHECK_NEAR(filterQFromNorm(-1.0f), 0.1, 1e-6);
+    CHECK_NEAR(filterQFromNorm(2.0f), 10.0, 1e-4);
+
+    // Strictly monotonic across the whole travel — no fold-back from the quadratic term.
+    float prev = -1.0f;
+    for (int i = 0; i <= 1000; ++i) {
+        const float q = filterQFromNorm(static_cast<float>(i) / 1000.0f);
+        CHECK(q > prev);
+        prev = q;
+    }
+}
+
+static void testQNormRoundTrips() {
+    for (int i = 0; i <= 20; ++i) {
+        const float n = static_cast<float>(i) / 20.0f;
+        CHECK_NEAR(filterNormFromQ(filterQFromNorm(n)), n, 1e-5);
+    }
+    CHECK_NEAR(filterNormFromQ(static_cast<float>(std::sqrt(2.0))), 0.5, 1e-5);
+    CHECK(filterNormFromQ(0.0f) == 0.0f);
+    CHECK(filterNormFromQ(1000.0f) == 1.0f);
+}
+
+// ---------------------------------------------------------------------------
+// Coefficients — pinned literals plus an independent textbook derivation
+// ---------------------------------------------------------------------------
+
+// Textbook RBJ Audio EQ Cookbook, computed straight from cos(w0)/sin(w0). Shares no code with
+// filter_coeffs, which reaches the same numbers through the tan half-angle substitution.
+static void rbjReference(bool highPass, double fc, double q, double sr, double out[5]) {
+    const double w0 = 2.0 * kPi * fc / sr;
+    const double c = std::cos(w0);
+    const double s = std::sin(w0);
+    const double alpha = s / (2.0 * q);
+    const double a0 = 1.0 + alpha;
+    const double n = highPass ? (1.0 + c) : (1.0 - c);
+    out[0] = n / 2.0 / a0;                        // b0
+    out[1] = (highPass ? -n : n) / a0;            // b1
+    out[2] = n / 2.0 / a0;                        // b2
+    out[3] = -2.0 * c / a0;                       // a1
+    out[4] = (1.0 - alpha) / a0;                  // a2
+}
+
+static void testCoefficientsMatchPinnedRbjValues() {
+    const double sr = 48000.0, fc = 1000.0, q = std::sqrt(2.0);
+
+    const BiquadCoeffs lp = biquadCoeffs(FilterMode::LowPass, static_cast<float>(fc),
+                                         static_cast<float>(q), sr);
+    const BiquadCoeffs hp = biquadCoeffs(FilterMode::HighPass, static_cast<float>(fc),
+                                         static_cast<float>(q), sr);
+
+    // Pinned literals: change the math and these fail.
+    CHECK_NEAR(lp.b0, 0.0040888771, 2e-6);
+    CHECK_NEAR(lp.b1, 0.0081777542, 2e-6);
+    CHECK_NEAR(lp.b2, 0.0040888771, 2e-6);
+    CHECK_NEAR(lp.a1, -1.8954199076, 2e-6);
+    CHECK_NEAR(lp.a2, 0.9117754318, 2e-6);
+
+    CHECK_NEAR(hp.b0, 0.9517988338, 2e-6);
+    CHECK_NEAR(hp.b1, -1.9035976676, 2e-6);
+    CHECK_NEAR(hp.b2, 0.9517988338, 2e-6);
+    CHECK_NEAR(hp.a1, -1.8954199076, 2e-6);
+    CHECK_NEAR(hp.a2, 0.9117754318, 2e-6);
+
+    // Independent derivation — proves the pinned literals are RBJ and not just "what we emit".
+    double ref[5];
+    rbjReference(false, fc, q, sr, ref);
+    CHECK_NEAR(lp.b0, ref[0], 1e-6);
+    CHECK_NEAR(lp.b1, ref[1], 1e-6);
+    CHECK_NEAR(lp.b2, ref[2], 1e-6);
+    CHECK_NEAR(lp.a1, ref[3], 1e-6);
+    CHECK_NEAR(lp.a2, ref[4], 1e-6);
+
+    rbjReference(true, fc, q, sr, ref);
+    CHECK_NEAR(hp.b0, ref[0], 1e-6);
+    CHECK_NEAR(hp.b1, ref[1], 1e-6);
+    CHECK_NEAR(hp.b2, ref[2], 1e-6);
+    CHECK_NEAR(hp.a1, ref[3], 1e-6);
+    CHECK_NEAR(hp.a2, ref[4], 1e-6);
+}
+
+static void testCoefficientsTrackSampleRateAndClampBelowNyquist() {
+    // Same fc at a different rate must give the RBJ answer for THAT rate, not a cached one.
+    double ref[5];
+    rbjReference(false, 1000.0, 2.0, 44100.0, ref);
+    const BiquadCoeffs at441 = biquadCoeffs(FilterMode::LowPass, 1000.0f, 2.0f, 44100.0);
+    CHECK_NEAR(at441.a1, ref[3], 1e-6);
+    CHECK_NEAR(at441.a2, ref[4], 1e-6);
+
+    // Requesting above 0.48*sr clamps rather than diverging through tan().
+    const BiquadCoeffs clamped = biquadCoeffs(FilterMode::LowPass, 20000.0f, 1.0f, 32000.0);
+    rbjReference(false, 0.48 * 32000.0, 1.0, 32000.0, ref);
+    CHECK_NEAR(clamped.b0, ref[0], 1e-6);
+    CHECK(std::isfinite(clamped.a1) && std::isfinite(clamped.a2));
+
+    // A non-positive rate passes through instead of inventing 44.1k.
+    const BiquadCoeffs bypass = biquadCoeffs(FilterMode::LowPass, 1000.0f, 1.0f, 0.0);
+    CHECK(bypass.b0 == 1.0f && bypass.b1 == 0.0f && bypass.b2 == 0.0f);
+    CHECK(bypass.a1 == 0.0f && bypass.a2 == 0.0f);
+}
+
+// DC gain of a lowpass and Nyquist gain of a highpass are both exactly unity — an independent
+// structural check on the coefficient set that a sign slip would break.
+static void testPassbandGainIsUnity() {
+    for (double q : {0.1, std::sqrt(2.0), 10.0}) {
+        const BiquadCoeffs lp =
+            biquadCoeffs(FilterMode::LowPass, 1000.0f, static_cast<float>(q), 48000.0);
+        CHECK_NEAR((lp.b0 + lp.b1 + lp.b2) / (1.0 + lp.a1 + lp.a2), 1.0, 1e-4);
+
+        const BiquadCoeffs hp =
+            biquadCoeffs(FilterMode::HighPass, 1000.0f, static_cast<float>(q), 48000.0);
+        CHECK_NEAR((hp.b0 - hp.b1 + hp.b2) / (1.0 - hp.a1 + hp.a2), 1.0, 1e-4);
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Resonance
+// ---------------------------------------------------------------------------
+
+// |H(e^jw)| for y = b0*x + b1*x1 + b2*x2 - a1*y1 - a2*y2.
+static double magnitudeAt(const BiquadCoeffs& c, double freqHz, double sr) {
+    const double w = 2.0 * kPi * freqHz / sr;
+    const double nRe = c.b0 + c.b1 * std::cos(w) + c.b2 * std::cos(2 * w);
+    const double nIm = -(c.b1 * std::sin(w) + c.b2 * std::sin(2 * w));
+    const double dRe = 1.0 + c.a1 * std::cos(w) + c.a2 * std::cos(2 * w);
+    const double dIm = -(c.a1 * std::sin(w) + c.a2 * std::sin(2 * w));
+    return std::sqrt(nRe * nRe + nIm * nIm) / std::sqrt(dRe * dRe + dIm * dIm);
+}
+
+static void testHighQPeaksAtCutoffInBothModes() {
+    const double sr = 48000.0, fc = 1000.0;
+    const float qHigh = filterQFromNorm(1.0f);   // 10
+    const float qLow = filterQFromNorm(0.0f);    // 0.1
+
+    for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) {
+        const BiquadCoeffs hi = biquadCoeffs(mode, static_cast<float>(fc), qHigh, sr);
+
+        // Scan a log grid and locate the maximum.
+        double peakMag = 0.0, peakFreq = 0.0;
+        for (int i = 0; i <= 600; ++i) {
+            const double f = 20.0 * std::pow(1000.0, static_cast<double>(i) / 600.0);
+            const double m = magnitudeAt(hi, f, sr);
+            if (m > peakMag) { peakMag = m; peakFreq = f; }
+        }
+        // The peak is at the cutoff, not at a band edge — within a quarter octave.
+        CHECK(peakFreq > fc / 1.19 && peakFreq < fc * 1.19);
+        // An RBJ 2-pole peaks at Q; assert most of that emphasis is really there.
+        CHECK(peakMag > 8.0);
+
+        // The emphasis is relative to the passband, not just a loud filter.
+        const double passband = magnitudeAt(hi, mode == FilterMode::LowPass ? 20.0 : 20000.0, sr);
+        CHECK_NEAR(passband, 1.0, 0.05);
+        CHECK(peakMag / passband > 8.0);
+
+        // At the bottom of the Q control there is no peak at all: the response is monotone
+        // over the band, so high Q is genuinely doing the work.
+        const BiquadCoeffs lo = biquadCoeffs(mode, static_cast<float>(fc), qLow, sr);
+        double prev = magnitudeAt(lo, 20.0, sr);
+        bool monotone = true;
+        for (int i = 1; i <= 600; ++i) {
+            const double f = 20.0 * std::pow(1000.0, static_cast<double>(i) / 600.0);
+            const double m = magnitudeAt(lo, f, sr);
+            if (mode == FilterMode::LowPass ? (m > prev + 1e-9) : (m < prev - 1e-9)) {
+                monotone = false;
+            }
+            prev = m;
+        }
+        CHECK(monotone);
+    }
+}
+
+// Drive real sines through VoiceFilter and measure steady-state RMS. Unlike the analytic
+// check above this also exercises the high-pass input-feedback path, which is outside the
+// coefficient transfer function.
+static double measuredRms(FilterMode mode, float cutoffNorm, float resNorm, double freqHz,
+                          double sr) {
+    VoiceFilter f;
+    f.prepare({mode, cutoffNorm, resNorm}, sr);
+    f.reset();
+
+    const int settle = 24000, measure = 24000;
+    double sumSq = 0.0;
+    for (int i = 0; i < settle + measure; ++i) {
+        const float x = static_cast<float>(std::sin(2.0 * kPi * freqHz * i / sr));
+        const float y = f.process(0, x);
+        if (i >= settle) sumSq += static_cast<double>(y) * y;
+    }
+    return std::sqrt(sumSq / measure);
+}
+
+static void testMeasuredResponsePeaksAtCutoffInBothModes() {
+    const double sr = 48000.0;
+    const float cutoffNorm = filterNormFromCutoffHz(1000.0f);
+
+    for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) {
+        double peakRms = 0.0, peakFreq = 0.0;
+        for (int i = 0; i <= 40; ++i) {
+            const double f = 100.0 * std::pow(100.0, static_cast<double>(i) / 40.0);
+            const double r = measuredRms(mode, cutoffNorm, 1.0f, f, sr);
+            if (r > peakRms) { peakRms = r; peakFreq = f; }
+        }
+        CHECK(peakFreq > 1000.0 / 1.3 && peakFreq < 1000.0 * 1.3);
+
+        const double passband =
+            measuredRms(mode, cutoffNorm, 1.0f, mode == FilterMode::LowPass ? 100.0 : 10000.0, sr);
+        CHECK(peakRms / passband > 3.0);
+
+        // Same measurement at the bottom of the resonance control shows no such emphasis.
+        const double flatAtCutoff = measuredRms(mode, cutoffNorm, 0.0f, 1000.0, sr);
+        const double flatPassband =
+            measuredRms(mode, cutoffNorm, 0.0f, mode == FilterMode::LowPass ? 100.0 : 10000.0, sr);
+        CHECK(flatAtCutoff / flatPassband < 1.0);
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Stability
+// ---------------------------------------------------------------------------
+
+static void testFullRangeCutoffSweepAtAudioRateStaysBounded() {
+    // Deterministic pseudo-noise; a fixed sine would miss the resonant frequency on most steps.
+    unsigned rng = 0x13579bdfu;
+    auto noise = [&rng]() {
+        rng = rng * 1664525u + 1013904223u;
+        return static_cast<float>(static_cast<int>(rng >> 9) - (1 << 22)) / static_cast<float>(1 << 22);
+    };
+
+    for (double sr : {44100.0, 48000.0, 96000.0}) {
+        for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) {
+            for (float res : {0.0f, 0.5f, 1.0f}) {
+                for (int direction = 0; direction < 2; ++direction) {
+                    VoiceFilter f;
+                    f.reset();
+                    const int n = 48000;
+                    for (int i = 0; i < n; ++i) {
+                        const float t = static_cast<float>(i) / static_cast<float>(n - 1);
+                        // Per-sample coefficient update across the whole cutoff travel.
+                        f.prepare({mode, direction == 0 ? t : 1.0f - t, res}, sr);
+                        const float y = f.process(0, noise());
+                        CHECK(std::isfinite(y));
+                        CHECK(std::fabs(y) < 100.0f);
+                        if (!std::isfinite(y)) return;  // stop before the log floods
+                    }
+                }
+            }
+        }
+    }
+}
+
+static void testStateFlushesToZeroWithoutStallingInDenormals() {
+    const double sr = 48000.0;
+    for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) {
+        VoiceFilter f;
+        f.prepare({mode, filterNormFromCutoffHz(1000.0f), 1.0f}, sr);
+        f.reset();
+
+        // Excite, then hard-cut to silence the way a released voice does.
+        for (int i = 0; i < 480; ++i) {
+            f.process(0, 0.5f * static_cast<float>(std::sin(2.0 * kPi * 1000.0 * i / sr)));
+        }
+
+        int subnormalSamples = 0;
+        int silentAt = -1;
+        for (int i = 0; i < 20000; ++i) {
+            f.process(0, 0.0f);
+            const VoiceFilter::State& s = f.state(0);
+            const float vals[5] = {s.x1, s.x2, s.y1, s.y2, s.fb};
+            for (float v : vals) {
+                if (v != 0.0f && std::fabs(v) < FLT_MIN) { ++subnormalSamples; break; }
+            }
+            if (silentAt < 0 && f.isSilent()) silentAt = i;
+        }
+        // Without the flush the state grinds down through the subnormal range for thousands
+        // of samples; a stray sample or two at a zero crossing is not a stall.
+        CHECK(subnormalSamples <= 2);
+        CHECK(silentAt >= 0);
+        CHECK(silentAt < 20000);
+        // And it stays silent — a flush that perturbs the feedback loop would re-excite it.
+        for (int i = 0; i < 1000; ++i) CHECK(f.process(0, 0.0f) == 0.0f);
+        CHECK(f.isSilent());
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Impulse / step sanity and saturation
+// ---------------------------------------------------------------------------
+
+static void testImpulseResponseMatchesDifferenceEquation() {
+    const double sr = 48000.0;
+    VoiceFilter f;
+    f.prepare({FilterMode::LowPass, filterNormFromCutoffHz(1000.0f), 0.5f}, sr);
+    f.reset();
+    const BiquadCoeffs c = f.coeffs();
+
+    // First three impulse-response taps follow directly from the coefficients.
+    const float h0 = f.process(0, 1.0f);
+    const float h1 = f.process(0, 0.0f);
+    const float h2 = f.process(0, 0.0f);
+    CHECK_NEAR(h0, c.b0, 1e-6);
+    CHECK_NEAR(h1, c.b1 - c.a1 * c.b0, 1e-6);
+    CHECK_NEAR(h2, c.b2 - c.a1 * h1 - c.a2 * h0, 1e-6);
+}
+
+static void testLowpassStepSettlesToUnity() {
+    const double sr = 48000.0;
+    VoiceFilter f;
+    f.prepare({FilterMode::LowPass, filterNormFromCutoffHz(1000.0f), 0.0f}, sr);
+    f.reset();
+    float y = 0.0f;
+    for (int i = 0; i < 48000; ++i) y = f.process(0, 1.0f);
+    CHECK_NEAR(y, 1.0, 1e-3);  // DC passes a lowpass at unity
+
+    VoiceFilter hp;
+    hp.prepare({FilterMode::HighPass, filterNormFromCutoffHz(1000.0f), 0.0f}, sr);
+    hp.reset();
+    for (int i = 0; i < 48000; ++i) y = hp.process(0, 1.0f);
+    CHECK_NEAR(y, 0.0, 1e-3);  // and is fully rejected by a highpass
+}
+
+static void testResetClearsHistoryButPrepareKeepsIt() {
+    VoiceFilter f;
+    f.prepare({FilterMode::LowPass, 0.5f, 0.5f}, 48000.0);
+    f.process(0, 1.0f);
+    CHECK(!f.isSilent());
+
+    // A live parameter move must not zero the history — that is what would click.
+    f.prepare({FilterMode::LowPass, 0.6f, 0.5f}, 48000.0);
+    CHECK(!f.isSilent());
+
+    f.reset();
+    CHECK(f.isSilent());
+}
+
+static void testChannelStateIsIndependent() {
+    VoiceFilter f;
+    f.prepare({FilterMode::LowPass, 0.5f, 0.5f}, 48000.0);
+    f.reset();
+    f.process(0, 1.0f);
+    CHECK(f.state(0).x1 == 1.0f);
+    CHECK(f.state(1).x1 == 0.0f);
+
+    float frame[2] = {1.0f, -1.0f};
+    f.processFrame(frame, 2);
+    CHECK(f.state(1).x1 == -1.0f);
+    CHECK(frame[0] != frame[1]);
+}
+
+static void testFeedbackSaturationIsContinuousAndBounded() {
+    CHECK_NEAR(feedbackSaturate(0.0f), 0.0, 1e-9);
+    // Odd symmetry.
+    CHECK_NEAR(feedbackSaturate(1.5f), -feedbackSaturate(-1.5f), 1e-6);
+    // Continuous across the threshold at +/-2.
+    CHECK_NEAR(feedbackSaturate(2.0f - 1e-4f), feedbackSaturate(2.0f + 1e-4f), 1e-4);
+    // Compresses hard: a 100x input does not give a 100x output.
+    CHECK(std::fabs(feedbackSaturate(100.0f)) < 12.0f);
+    CHECK(feedbackSaturate(100.0f) > feedbackSaturate(50.0f));
+}
+
+int main() {
+    testCutoffMapsThreeDecadesLogarithmically();
+    testCutoffNormRoundTrips();
+    testQSpansPointOneToTenWithRootTwoAtCenter();
+    testQNormRoundTrips();
+    testCoefficientsMatchPinnedRbjValues();
+    testCoefficientsTrackSampleRateAndClampBelowNyquist();
+    testPassbandGainIsUnity();
+    testHighQPeaksAtCutoffInBothModes();
+    testMeasuredResponsePeaksAtCutoffInBothModes();
+    testFullRangeCutoffSweepAtAudioRateStaysBounded();
+    testStateFlushesToZeroWithoutStallingInDenormals();
+    testImpulseResponseMatchesDifferenceEquation();
+    testLowpassStepSettlesToUnity();
+    testResetClearsHistoryButPrepareKeepsIt();
+    testChannelStateIsIndependent();
+    testFeedbackSaturationIsContinuousAndBounded();
+
+    if (g_fail == 0) std::printf("filter_tests: all passed\n");
+    else std::printf("filter_tests: %d FAILED\n", g_fail);
+    return g_fail == 0 ? 0 : 1;
+}