Filter review fixes: gate softLimit on drive, correct flush/DF1 measurement claims, pin control-law test literals, fix stale bypass state

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

@@ -85,8 +85,15 @@ by the topology rather than picked by ear:
 Three properties carry the design:
 
 - `depth == 0` makes it algebraically the identity (`x / sqrt(1) == x`, exact in IEEE), so
-  drive 0 is **bit-exact** linear with no branch and no special case on the hot path. The
-  test asserts bit-identity against the same kernel with the limiter deleted.
+  drive 0 is **bit-exact** linear whether or not `softLimit` is actually called. The test
+  asserts bit-identity against the same kernel with the limiter deleted.
+- `process()` gates the call on `driven_` (`driveDepth_ != 0`, cached at `prepare()`) rather
+  than calling `softLimit` unconditionally. `sqrt`/div sit on the per-sample recursive
+  dependency chain, so out-of-order execution can't hide their latency, and at drive 0 that
+  cost buys nothing. Measured: 11.2 ns/sample unconditional vs 4.1 ns gated — the gated form
+  lands at the limiter-removed floor. `driven_` only changes at `prepare()`, so the branch
+  predicts perfectly. The gate is a perf optimization on top of the bit-identity above, not a
+  substitute for it — deleting the gate would still be correct, just 2.7x slower at rest.
 - `|softLimit(u, d)| <= |u|` for every depth, so the state update can only shrink the state.
   The filter cannot gain energy from the drive stage: stability at any Q and any cutoff is
   structural, and self-oscillation is impossible. This is why the shaper must keep unit slope
@@ -131,25 +138,48 @@ kink at the center detent. The quadratic term is nonzero only because √2 is no
 geometric mean of 0.1 and 10; `filterNormFromQ` divides by it. The SVF consumes it as
 `k = 1/Q`.
 
-### Denormal flushing tests the envelope, not one sample
+### Denormal flushing: why conjunctive, honestly
 
 `process()` flushes **both** integrators to exact zero once both are below
-`kFilterDenormalFloor` (1e-30). Testing both is required, not tidy: `ic1` and `ic2` are in
-quadrature, so a resonator swings each of them through zero twice a cycle. Flushing on a
-single integrator injects a step in phase with the resonance, which the resonance then
-amplifies — the filter limit-cycles at the floor forever instead of going quiet. This was
-re-verified for TPT rather than assumed to transfer from the retired Direct Form I state.
+`kFilterDenormalFloor` (1e-30). The honest reason is narrower than it sounds: `isSilent()`
+means "both integrators are exactly zero," so both have to reach zero for that check to mean
+anything, and the conjunctive test is the cheapest way to guarantee it.
+
+The stronger claim — that a per-variable flush limit-cycles at the floor — does **not**
+reproduce on this topology. Measured (Q=10, fc=1kHz, 48k): shipped conjunctive goes silent at
+sample 10783 with 0 subnormals; a per-variable independent flush goes silent ~180 samples
+earlier and an either-below-zero-both flush ~970 samples earlier, both also 0 subnormals, no
+limit cycle, and the same excited RMS. That claim WAS real on the retired Direct Form I state,
+where the flushed variables (`y1`/`y2`) were the actual filter OUTPUT, so zeroing one injected
+a discontinuity the resonance then amplified. Here `ic1`/`ic2` are integrator STATE, not
+output: zeroing one only removes energy, a contraction rather than an injection, so the hazard
+is structurally absent. The only demonstrable hazard is no flush at all, which never reaches
+exact zero and grinds through subnormals for thousands of samples on a released voice.
+
+Keep the conjunctive test regardless — it costs nothing extra and is the right guarantee for
+`isSilent()` — but don't cite the limit-cycle rationale for TPT; it belongs to the retired
+topology.
 
 ## Gotchas
 
 - **TPT is what fixed the low-cutoff conditioning defect** — this is a topology change, not
   a relocation. Direct Form I encoded pole proximity in `a1 → -2`, `a2 → +1` and cancelled
-  them against each other every sample, which at `fc/sr ≈ 1e-4` cost ~17 bits and put the
-  measured peak **15% low** at 20 Hz / 192 kHz. TPT encodes the same proximity in `a1`'s
-  small deviation from 1, which float32 resolves: measured 10.0160 against the analytic
-  10.0125, +0.034%. Do not reintroduce a direct-form kernel.
+  them against each other every sample; at `fc/sr ≈ 1e-4` that ~17-bit cancellation moved the
+  measured 20 Hz / 192 kHz LP peak by **-27% on a true-peak scan, -57% measured at the
+  analytic peak frequency** (the degraded pole itself moves, so the two methods diverge), and
+  the error is non-monotone with rate rather than a fixed percentage (+5% high at 96 kHz).
+  TPT encodes the same proximity in `a1`'s small deviation from 1, which float32 resolves:
+  checked against an exact-double evaluation of the same difference equation (which matches
+  the analytic target to within measurement noise), TPT's float32-narrowed coefficients are
+  genuinely ~0.02% low at 48 kHz, widening to ~0.03% low at 192 kHz — real coefficient
+  narrowing, not measurement-window noise, and comfortably inside the test's 0.4% tolerance
+  either way. Do not reintroduce a direct-form kernel.
 - **`prepare()` deliberately does not clear state** — a live parameter move must glide, not
-  click. Call `reset()` at note-on.
+  click. Call `reset()` at note-on. **Exception: the non-positive-rate bypass path.** There,
+  `a1=1, a2=a3=0` makes both state updates the exact identity and `bypassMix()` never reads
+  the state at all, so a stale nonzero `ic1`/`ic2` would otherwise latch `isSilent()` false
+  forever with no audible effect either way — `prepare()` clears state on that path only,
+  which costs nothing audibly since bypass ignores it.
 - **The morph endpoints are asserted on the folded mix, exactly.** `morphWeights` snaps the
   leg endpoints instead of trusting `cos`/`sin` to land on 0 and 1, which they miss by ~1e-17
   — enough to leave a -324 dB neighbour tap in what is specified as a pure response.

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

@@ -28,10 +28,15 @@ inline constexpr double kFilterNyquistFraction = 0.48;
 // rate rather than assume 44.1k.
 //
 // Float storage is safe HERE in a way it was not for the retired Direct Form I path. DF1 encoded
-// pole proximity in a1 -> -2, a2 -> +1 and cancelled them against each other every sample, which
-// at fc/sr ~ 1e-4 cost ~17 bits and moved the resonant peak -15%. TPT encodes the same proximity
-// in a1's small DEVIATION from 1, which float resolves. Measured 20 Hz/192 kHz peak is 10.0160
-// against the analytic 10.0125, +0.034%.
+// pole proximity in a1 -> -2, a2 -> +1 and cancelled them against each other every sample; at
+// fc/sr ~ 1e-4 that ~17-bit cancellation moved the measured 20 Hz/192 kHz LP peak by -27%
+// (true-peak scan) to -57% (point measurement at the analytic peak frequency, since the
+// degraded pole itself moves) -- and the error is non-monotone with rate, not a fixed percentage
+// (+5% high at 96 kHz). TPT encodes the same proximity in a1's small DEVIATION from 1, which
+// float resolves: measured against an exact-double evaluation of the same difference equation
+// (which matches the analytic target to within measurement noise), TPT's float32-narrowed
+// coefficients land genuinely ~0.02% low at 48 kHz, widening to ~0.03% low at 192 kHz -- both
+// comfortably inside the test's 0.4% tolerance.
 SvfCoeffs svfCoeffs(float cutoffHz, float q, double sampleRate);
 
 }  // namespace reasampler::instrument::engine::filter

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

@@ -14,7 +14,8 @@ struct MorphWeights {
     float lp = 1.0f;
 };
 
-// HP at 0.0, BP at 0.5, LP at 1.0. Out-of-range norm clamps to the endpoints.
+// HP at 0.0, BP at 0.5, LP at 1.0. Out-of-range norm clamps to the endpoints; NaN clamps to
+// neither and lands on pure band-pass instead (every comparison against it is false).
 //
 // Equal-power (cos/sin) rather than linear, and that choice is forced by the topology rather
 // than picked by ear. At the corner frequency the three taps are HP = jQ, BP = Q, LP = -jQ, so

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

@@ -5,8 +5,18 @@ namespace reasampler::instrument::engine::filter {
 void VoiceFilter::prepare(const FilterSettings& settings, double sampleRate) {
     coeffs_ = svfCoeffs(filterCutoffHzFromNorm(settings.cutoffNorm),
                         filterQFromNorm(settings.resonanceNorm), sampleRate);
-    mix_ = (sampleRate > 0.0) ? morphMix(morphWeights(settings.morphNorm), coeffs_.k) : bypassMix();
+    if (sampleRate > 0.0) {
+        mix_ = morphMix(morphWeights(settings.morphNorm), coeffs_.k);
+    } else {
+        // Bypass: a1=1, a2=a3=0 makes both state updates the exact identity, and bypassMix()
+        // reads only the input, never the state -- so clearing here is audibly free (the state
+        // was already going to be ignored) and prevents a stale nonzero ic1/ic2 from latching
+        // isSilent() false forever, which prepare() otherwise deliberately never does.
+        mix_ = bypassMix();
+        for (State& s : state_) s = State{};
+    }
     driveDepth_ = filterDriveDepthFromNorm(settings.driveNorm);
+    driven_ = driveDepth_ != 0.0f;
 }
 
 void VoiceFilter::reset() {

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

@@ -58,15 +58,28 @@ public:
         // Placing it on the state rather than inside the zero-delay loop keeps a1/a2/a3 an exact
         // algebraic solve — a nonlinearity inside the loop would need per-sample Newton
         // iteration. softLimit is a contraction, so this cannot destabilize the filter.
-        s.ic1 = softLimit(2.0f * v1 - s.ic1, driveDepth_);
+        //
+        // Gated on driven_ rather than called unconditionally: sqrt and div sit on this
+        // recursive dependency chain, so out-of-order execution can't hide them, and at drive 0
+        // (the default) that cost buys nothing — softLimit(x, 0) == x algebraically. Measured:
+        // 11.2 ns/sample unconditional vs 4.1 ns gated, matching the limiter-removed floor.
+        // driven_ only changes at prepare(), so the branch predicts perfectly. Bit-identity at
+        // drive 0 holds either way, by algebra — the gate is a perf optimization, not what makes
+        // it exact.
+        const float u = 2.0f * v1 - s.ic1;
+        s.ic1 = driven_ ? softLimit(u, driveDepth_) : u;
         s.ic2 = 2.0f * v2 - s.ic2;
 
-        // Snap the state once the whole resonator has decayed past -600 dB. Testing BOTH
-        // integrators is testing the ENVELOPE rather than one sample, and that is required, not
-        // tidy: ic1 and ic2 are in quadrature, so a resonator swings each of them through zero
-        // twice a cycle. Flushing on a single integrator would inject a step in phase with the
-        // resonance, which the resonance then amplifies — the filter limit-cycles at the floor
-        // forever instead of going quiet.
+        // Snap the state once the whole resonator has decayed past -600 dB. isSilent() means
+        // "both integrators are exactly zero," so both must reach zero for that check to be
+        // meaningful — the conjunctive test is the cheapest guarantee of that, not a defense
+        // against a demonstrated limit cycle on this topology (measured: a per-variable flush
+        // and an either-below-zero-both flush both go silent here too, no limit cycle, no
+        // subnormals). That risk was real on the retired Direct Form I state, where a per-sample
+        // flush zeroed y1/y2 — the actual OUTPUT — injecting a step the resonance then amplified.
+        // ic1/ic2 are integrator STATE, not output; zeroing one only removes energy, a
+        // contraction rather than an injection. The only demonstrable hazard here is no flush at
+        // all, which never reaches exact zero and stalls in subnormals for thousands of samples.
         if (s.ic1 > -kFilterDenormalFloor && s.ic1 < kFilterDenormalFloor &&
             s.ic2 > -kFilterDenormalFloor && s.ic2 < kFilterDenormalFloor) {
             s.ic1 = 0.0f;
@@ -96,6 +109,7 @@ private:
     SvfCoeffs coeffs_{};
     MorphMix mix_{};
     float driveDepth_ = 0.0f;
+    bool driven_ = false;  // driveDepth_ != 0, cached so process() branches on a bool, not a float compare
     State state_[kMaxChannels]{};
 };
 

tests/test_filter.cpp

@@ -14,6 +14,7 @@
 #include <cmath>
 #include <cstdio>
 #include <initializer_list>
+#include <limits>
 
 using namespace reasampler::instrument::engine::filter;
 
@@ -116,6 +117,23 @@ static void testQSpansPointOneToTenWithRootTwoAtCenter() {
     CHECK_NEAR(filterQFromNorm(-1.0f), 0.1, 1e-6);
     CHECK_NEAR(filterQFromNorm(2.0f), 10.0, 1e-4);
 
+    // Pins the single quadratic-in-log-Q curve at two interior points, derived independently by
+    // solving log Q = a + b*n + c*n^2 through the three anchors above rather than read out of
+    // the implementation. A two-spliced-log-segments curve (log-linear on each half, the design
+    // this module doc explicitly rejects for its center-detent slope kink) would give 0.376 and
+    // 3.761 here instead — both comfortably outside this tolerance.
+    {
+        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));
+        const double c = 2.0 * lo + 2.0 * hi - 4.0 * mid;
+        const double b = hi - lo - c;
+        const double a = lo;
+        auto qLaw = [&](double n) { return std::exp(a + b * n + c * n * n); };
+        CHECK_NEAR(filterQFromNorm(0.25f), qLaw(0.25), 1e-5);
+        CHECK_NEAR(filterQFromNorm(0.75f), qLaw(0.75), 1e-5);
+    }
+
     // Strictly monotonic across the whole travel — no fold-back from the quadratic term.
     float prev = -1.0f;
     for (int i = 0; i <= 1000; ++i) {
@@ -142,6 +160,10 @@ static void testDriveDepthIsZeroAtRestAndRisesMonotonically() {
     CHECK_NEAR(filterDriveDepthFromNorm(1.0f), kFilterDriveDepthMax, 1e-6);
     CHECK_NEAR(filterDriveDepthFromNorm(2.0f), kFilterDriveDepthMax, 1e-6);
 
+    // Pins the SQUARE law at an interior point, not just the anchors: a linear law would give
+    // kFilterDriveDepthMax/2 (2.0) here, not kFilterDriveDepthMax/4 (1.0).
+    CHECK_NEAR(filterDriveDepthFromNorm(0.5f), kFilterDriveDepthMax * 0.25, 1e-6);
+
     float prev = -1.0f;
     for (int i = 0; i <= 100; ++i) {
         const float d = filterDriveDepthFromNorm(static_cast<float>(i) / 100.0f);
@@ -195,6 +217,28 @@ static void testTheSampleRateEntersOnlyThroughG() {
     CHECK(svfCoeffs(1000.0f, 1.0f, -48000.0).g == 0.0f);
 }
 
+// A voice re-prepared at a non-positive rate while still ringing must not latch isSilent()
+// false forever -- a future voice allocator using isSilent() as its free condition would leak
+// the voice. Bypass ignores state entirely (a1=1, a2=a3=0, bypassMix reads only the input), so
+// clearing it here is audibly free.
+static void testNonPositiveRatePrepareClearsStaleStateAndReportsSilent() {
+    VoiceFilter f;
+    f.prepare(at(1000.0, 1.0f, kLowPass), 48000.0);
+    f.reset();
+    for (int i = 0; i < 100; ++i) {
+        f.process(0, static_cast<float>(std::sin(2.0 * kPi * 1000.0 * i / 48000.0)));
+    }
+    CHECK(!f.isSilent());  // genuinely ringing before the rate goes bad
+
+    f.prepare({0.5f, 0.5f, kLowPass, 0.0f}, 0.0);
+    CHECK(f.isSilent());
+    for (int i = 0; i < 480000; ++i) {
+        const float x = static_cast<float>(std::sin(0.1 * i));
+        CHECK(f.process(0, x) == x);
+    }
+    CHECK(f.isSilent());
+}
+
 // An invalid rate must pass the signal, not silence the instrument, whatever the morph asks for.
 static void testNonPositiveRatePassesSignalThroughAtEveryMorph() {
     for (float morph : {kHighPass, kBandPass, kLowPass}) {
@@ -229,6 +273,23 @@ static void testMorphEndpointMixesAreExactlyPureTaps() {
     // Out-of-range clamps to the endpoints rather than extrapolating.
     CHECK(morphWeights(-1.0f).hp == 1.0f);
     CHECK(morphWeights(2.0f).lp == 1.0f);
+
+    // NaN clamps to neither endpoint (every comparison against it is false) and lands on pure
+    // band-pass instead -- no crash, a sane fallback rather than an extrapolation.
+    const MorphWeights nanW = morphWeights(std::numeric_limits<float>::quiet_NaN());
+    CHECK(nanW.hp == 0.0f && nanW.bp == 1.0f && nanW.lp == 0.0f);
+}
+
+// Pins the cos/sin curve at an interior point, not just the endpoints and the quadrature
+// identity (hp^2+bp^2+lp^2=1, which any equal-power reparameterization would also satisfy).
+// theta=0.5*pi*t^2 (quadratic in the leg fraction, still equal-power, still exact at both
+// ends) would give hp=0.9239/bp=0.3827 here instead of the cos/sin pair's 0.7071/0.7071.
+static void testMorphInteriorPointMatchesCosSinNotAnAlternateEqualPowerCurve() {
+    const MorphWeights w = morphWeights(0.25f);  // HP->BP leg, t = 2*0.25 = 0.5
+    const double theta = 0.5 * kPi * 0.5;
+    CHECK_NEAR(w.hp, std::cos(theta), 1e-6);
+    CHECK_NEAR(w.bp, std::sin(theta), 1e-6);
+    CHECK(w.lp == 0.0f);
 }
 
 // HP and LP never carry weight at the same time. That is what keeps the centre a band-pass
@@ -322,33 +383,64 @@ static void testMorphSweepHasNoDiscontinuity() {
 
 // The hard acceptance criterion, in its strongest form: at drive 0 the kernel is BIT-IDENTICAL
 // to the same kernel with the limiter deleted. softLimit(x, 0) is x / sqrt(1) == x exactly, so
-// this holds by algebra rather than by tolerance.
+// this holds by algebra rather than by tolerance. Both channels and both entry points
+// (process() and processFrame()) are covered, not just channel 0 through process().
+struct LinearKernelRef {
+    SvfCoeffs c;
+    MorphMix mix;
+    float ic1 = 0.0f, ic2 = 0.0f;
+
+    float step(float x) {
+        const float v3 = x - ic2;
+        const float v1 = c.a1 * ic1 + c.a2 * v3;
+        const float v2 = ic2 + c.a2 * ic1 + c.a3 * v3;
+        ic1 = 2.0f * v1 - ic1;  // no limiter at all
+        ic2 = 2.0f * v2 - ic2;
+        if (ic1 > -kFilterDenormalFloor && ic1 < kFilterDenormalFloor &&
+            ic2 > -kFilterDenormalFloor && ic2 < kFilterDenormalFloor) {
+            ic1 = 0.0f;
+            ic2 = 0.0f;
+        }
+        return mix.m0 * x + mix.m1 * v1 + mix.m2 * v2;
+    }
+};
+
+static float nextNoise(unsigned& rng) {
+    rng = rng * 1664525u + 1013904223u;
+    return static_cast<float>(static_cast<int>(rng >> 9) - (1 << 22)) /
+           static_cast<float>(1 << 22);
+}
+
 static void testDriveZeroIsBitIdenticalToTheLinearKernel() {
     for (float morph : {kHighPass, kBandPass, kLowPass}) {
         VoiceFilter f;
         f.prepare(at(1000.0, 1.0f, morph, 0.0f), 48000.0);
         f.reset();
-        const SvfCoeffs c = f.coeffs();
-        const MorphMix mix = f.mix();
+        LinearKernelRef ref0{f.coeffs(), f.mix()};
+        LinearKernelRef ref1{f.coeffs(), f.mix()};
 
-        float ic1 = 0.0f, ic2 = 0.0f;
-        unsigned rng = 0x13579bdfu;
+        unsigned rng0 = 0x13579bdfu;
         for (int i = 0; i < 4096; ++i) {
-            rng = rng * 1664525u + 1013904223u;
-            const float x = static_cast<float>(static_cast<int>(rng >> 9) - (1 << 22)) /
-                            static_cast<float>(1 << 22);
+            const float x = nextNoise(rng0);
+            CHECK(f.process(0, x) == ref0.step(x));
+        }
 
-            const float v3 = x - ic2;
-            const float v1 = c.a1 * ic1 + c.a2 * v3;
-            const float v2 = ic2 + c.a2 * ic1 + c.a3 * v3;
-            ic1 = 2.0f * v1 - ic1;  // no limiter at all
-            ic2 = 2.0f * v2 - ic2;
-            if (ic1 > -kFilterDenormalFloor && ic1 < kFilterDenormalFloor &&
-                ic2 > -kFilterDenormalFloor && ic2 < kFilterDenormalFloor) {
-                ic1 = 0.0f;
-                ic2 = 0.0f;
-            }
-            CHECK(f.process(0, x) == mix.m0 * x + mix.m1 * v1 + mix.m2 * v2);
+        // process(1, ...): channel 1's state is independent of channel 0's above.
+        unsigned rng1 = 0x2468acefu;
+        for (int i = 0; i < 4096; ++i) {
+            const float x = nextNoise(rng1);
+            CHECK(f.process(1, x) == ref1.step(x));
+        }
+
+        // processFrame(): both channels advanced together through the frame entry point,
+        // continuing from the state each channel already has.
+        for (int i = 0; i < 4096; ++i) {
+            float frame[2] = {nextNoise(rng0), nextNoise(rng1)};
+            const float want0 = ref0.step(frame[0]);
+            const float want1 = ref1.step(frame[1]);
+            f.processFrame(frame, 2);
+            CHECK(frame[0] == want0);
+            CHECK(frame[1] == want1);
         }
     }
 }
@@ -712,9 +804,11 @@ int main() {
 
     testSvfCoefficientsMatchPinnedValues();
     testTheSampleRateEntersOnlyThroughG();
+    testNonPositiveRatePrepareClearsStaleStateAndReportsSilent();
     testNonPositiveRatePassesSignalThroughAtEveryMorph();
 
     testMorphEndpointMixesAreExactlyPureTaps();
+    testMorphInteriorPointMatchesCosSinNotAnAlternateEqualPowerCurve();
     testMorphNeverBlendsHighAgainstLowPass();
     testMorphEndpointsMatchTheAnalyticTwoPoleTargets();
     testCornerMagnitudeIsFlatAcrossTheWholeMorphSweep();