diff --git a/src/core/instrument/engine/filter/CLAUDE.md b/src/core/instrument/engine/filter/CLAUDE.md index 3c09ce2..577ff38 100644 --- a/src/core/instrument/engine/filter/CLAUDE.md +++ b/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. diff --git a/src/core/instrument/engine/filter/filter_coeffs.h b/src/core/instrument/engine/filter/filter_coeffs.h index 0a483d7..8db6869 100644 --- a/src/core/instrument/engine/filter/filter_coeffs.h +++ b/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 diff --git a/src/core/instrument/engine/filter/filter_morph.h b/src/core/instrument/engine/filter/filter_morph.h index d622a30..15b484f 100644 --- a/src/core/instrument/engine/filter/filter_morph.h +++ b/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 diff --git a/src/core/instrument/engine/filter/voice_filter.cpp b/src/core/instrument/engine/filter/voice_filter.cpp index 69088ec..2b0f369 100644 --- a/src/core/instrument/engine/filter/voice_filter.cpp +++ b/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() { diff --git a/src/core/instrument/engine/filter/voice_filter.h b/src/core/instrument/engine/filter/voice_filter.h index 5c91474..079df92 100644 --- a/src/core/instrument/engine/filter/voice_filter.h +++ b/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]{}; }; diff --git a/tests/test_filter.cpp b/tests/test_filter.cpp index 839525d..390bf7b 100644 --- a/tests/test_filter.cpp +++ b/tests/test_filter.cpp @@ -14,6 +14,7 @@ #include #include #include +#include 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(kFilterQMin)); + const double mid = std::log(static_cast(kFilterQCenter)); + const double hi = std::log(static_cast(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(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(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(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::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(static_cast(rng >> 9) - (1 << 22)) / + static_cast(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(static_cast(rng >> 9) - (1 << 22)) / - static_cast(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();