// Standalone tests for the RUNNING per-voice TPT/SVF 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 derivation
// that shares no code with the implementation, and the responses against the analog 2-pole
// prototype evaluated at the bilinear-warped frequency. Sibling targets own the neighbouring
// domains: test_filter_params.cpp the control mappings, test_filter_morph.cpp the pure morph-weight
// algebra, test_filter_state.cpp the numerical/state behaviour. This file owns the analytic
// reference and the steady-state gain measurement, and everything here uses them.

#include "../src/core/instrument/engine/filter/filter_coeffs.h"
#include "../src/core/instrument/engine/filter/filter_morph.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 <cmath>
#include <cstdio>
#include <initializer_list>

using namespace reasampler::instrument::engine::filter;

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;

// Morph positions. The endpoints are the same pure taps under both laws; only the centre differs
// — a band-pass under HighBandLow, a notch under HighNotchLow.
static constexpr float kHighPass = 0.0f;
static constexpr float kBandPass = 0.5f;
static constexpr float kCentre = 0.5f;
static constexpr float kLowPass = 1.0f;

static const MorphLaw kBothLaws[] = {MorphLaw::HighBandLow, MorphLaw::HighNotchLow};
static const char* lawName(MorphLaw law) {
    return law == MorphLaw::HighBandLow ? "HP-BP-LP" : "HP-notch-LP";
}

// The rates the invariance claims are made over.
static const double kRates[] = {44100.0, 48000.0, 88200.0, 96000.0, 192000.0};
static constexpr int kRateCount = 5;

// The measurement pass's bar, and the bar the rewrite exists to hold: peak and passband agree
// with the analytic target to better than this at every rate, level, and morph position.
static constexpr double kAgreement = 0.004;

// ---------------------------------------------------------------------------
// Independent references
// ---------------------------------------------------------------------------

// The analog 2-pole prototype |H(jW)| evaluated at the bilinear-warped frequency. The TPT maps
// the digital frequency onto the prototype EXACTLY at the prewarped corner, so this is the exact
// digital magnitude — derived from the continuous-time prototype and the transform rather than
// from anything filter_coeffs computes.
static double analyticMag(float morph, double freq, double fc, double q, double sr) {
    const double w = std::tan(kPi * freq / sr) / std::tan(kPi * fc / sr);
    const double dRe = 1.0 - w * w, dIm = w / q;
    const double den = std::sqrt(dRe * dRe + dIm * dIm);
    if (morph == kHighPass) return w * w / den;
    if (morph == kBandPass) return w / den;
    return 1.0 / den;
}

// Steady-state gain of the running filter at one frequency. Windows are wall-clock, not sample
// counts, so every rate integrates the same amount of signal.
static double measuredGain(const FilterSettings& fs, double sr, double freq, double amp = 0.25,
                           double settleSec = 0.15, double measureSec = 0.10) {
    VoiceFilter f;
    f.prepare(fs, sr);
    f.reset();
    const int settle = static_cast<int>(sr * settleSec);
    const int measure = static_cast<int>(sr * measureSec);
    double sumSq = 0.0;
    for (int i = 0; i < settle + measure; ++i) {
        const float y = f.process(0, static_cast<float>(amp * std::sin(2.0 * kPi * freq * i / sr)));
        if (i >= settle) sumSq += static_cast<double>(y) * y;
    }
    return std::sqrt(sumSq / measure) / (amp / std::sqrt(2.0));
}

static FilterSettings at(double fcHz, float res, float morph, float drive = 0.0f,
                         MorphLaw law = MorphLaw::HighBandLow) {
    return {filterNormFromCutoffHz(static_cast<float>(fcHz)), res, morph, drive, law};
}

// ---------------------------------------------------------------------------
// SVF coefficients — pinned literals plus an independent derivation
// ---------------------------------------------------------------------------

static void testSvfCoefficientsMatchPinnedValues() {
    const double sr = 48000.0, fc = 1000.0, q = std::sqrt(2.0);
    const SvfCoeffs c = svfCoeffs(static_cast<float>(fc), static_cast<float>(q), sr);

    // Pinned literals: change the math and these fail.
    CHECK_NEAR(c.g, 0.0655434653, 2e-9);
    CHECK_NEAR(c.k, 0.7071067691, 2e-9);
    CHECK_NEAR(c.a1, 0.9517988563, 2e-9);
    CHECK_NEAR(c.a2, 0.0623841919, 2e-9);
    CHECK_NEAR(c.a3, 0.0040888758, 2e-9);

    // Independent derivation — proves the pins are the TPT solve and not just "what we emit".
    const double g = std::tan(kPi * fc / sr);
    const double k = 1.0 / q;
    const double denom = 1.0 + g * g + g * k;  // written out rather than factored as g*(g+k)
    CHECK_NEAR(c.g, g, 1e-7);
    CHECK_NEAR(c.k, k, 1e-7);
    CHECK_NEAR(c.a1, 1.0 / denom, 1e-7);
    CHECK_NEAR(c.a2, g / denom, 1e-7);
    CHECK_NEAR(c.a3, g * g / denom, 1e-7);
}

static void testTheSampleRateEntersOnlyThroughG() {
    // k and the cutoff mapping are rate-free; only g moves with the rate. A reference rate
    // creeping back into the module would break this.
    const SvfCoeffs a = svfCoeffs(1000.0f, 2.0f, 48000.0);
    const SvfCoeffs b = svfCoeffs(1000.0f, 2.0f, 96000.0);
    CHECK(a.k == b.k);
    CHECK(a.g != b.g);
    CHECK_NEAR(b.g, std::tan(kPi * 1000.0 / 96000.0), 1e-7);

    // Requesting above 0.48*sr clamps rather than diverging through tan().
    const SvfCoeffs clamped = svfCoeffs(20000.0f, 1.0f, 32000.0);
    CHECK_NEAR(clamped.g, std::tan(kPi * 0.48), 1e-5);
    CHECK(std::isfinite(clamped.a1) && std::isfinite(clamped.a3));

    // A non-positive rate yields g == 0 instead of inventing 44.1k.
    CHECK(svfCoeffs(1000.0f, 1.0f, 0.0).g == 0.0f);
    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}) {
        VoiceFilter f;
        f.prepare({0.5f, 0.5f, morph, 0.0f}, 0.0);
        f.reset();
        for (int i = 0; i < 64; ++i) {
            const float x = static_cast<float>(std::sin(0.1 * i));
            CHECK(f.process(0, x) == x);
        }
    }
}

// ---------------------------------------------------------------------------
// Morph — measured, under both laws
// ---------------------------------------------------------------------------

// The endpoints are exact 2-pole HP and LP under BOTH laws; only the centre is law-specific, so
// the centre is asserted here only for the law that has a pure tap there.
static void testMorphEndpointsMatchTheAnalyticTwoPoleTargets() {
    const double sr = 48000.0, fc = 1000.0;
    for (MorphLaw law : kBothLaws) {
        for (float res : {0.0f, 0.5f, 1.0f}) {
            const double q = filterQFromNorm(res);
            for (float morph : {kHighPass, kBandPass, kLowPass}) {
                if (morph == kBandPass && law != MorphLaw::HighBandLow) continue;
                for (double f : {125.0, 500.0, 1000.0, 2000.0, 8000.0}) {
                    const double got = measuredGain(at(fc, res, morph, 0.0f, law), sr, f);
                    const double want = analyticMag(morph, f, fc, q, sr);
                    if (!(std::fabs(got / want - 1.0) <= kAgreement)) {
                        std::printf("FAIL line %d: %s morph %.1f res %.1f at %.0f Hz: %.6f vs "
                                    "analytic %.6f (%.3f%%)\n",
                                    __LINE__, lawName(law), morph, res, f, got, want,
                                    (got / want - 1.0) * 100.0);
                        ++g_fail;
                    }
                }
            }
        }
    }
}

// LAW-SPECIFIC, and deliberately not generalized: this guarantee belongs to HighBandLow alone.
// At the corner the three taps are HP = jQ, BP = Q, LP = -jQ — ADJACENT taps in exact quadrature
// — so a cos/sin pair holds the corner magnitude at exactly Q the whole way across. A linear
// crossfade would sag to Q/sqrt(2) mid-leg, a 3 dB hole that would read as a defect rather than
// as character. HighNotchLow deliberately violates this (its corner magnitude goes to zero at the
// centre); weakening this assertion to accommodate that law would throw the guarantee away.
static void testCornerMagnitudeIsFlatAtQAcrossTheHighBandLowSweep() {
    const double sr = 48000.0, fc = 1000.0;
    for (float res : {0.0f, 0.5f, 1.0f}) {
        const double q = filterQFromNorm(res);
        for (int i = 0; i <= 16; ++i) {
            const float m = static_cast<float>(i) / 16.0f;
            const double got = measuredGain(at(fc, res, m, 0.0f, MorphLaw::HighBandLow), sr, fc);
            if (!(std::fabs(got / q - 1.0) <= kAgreement)) {
                std::printf("FAIL line %d: morph %.4f res %.1f corner gain %.6f, expected Q "
                            "%.6f (%.3f%%)\n",
                            __LINE__, m, res, got, q, (got / q - 1.0) * 100.0);
                ++g_fail;
            }
        }
    }
}

// The SEM's centre is a genuine null, not merely a dip: the corner magnitude falls to the float
// noise floor because HP and LP sit at exactly +90 and -90 degrees there, so equal weights cancel
// by construction. Measured worst case across this whole grid is 3.8e-05 (-88 dB); the typical
// figure is -110 to -145 dB. The settle window has to clear the resonator's ring-down before the
// residual means anything — at 0.15 s and Q=10 the leftover transient alone reads as -52 dB and
// would be mistaken for the floor.
static void testHighNotchLowCentreIsATrueNullAtTheCorner() {
    for (int r = 0; r < kRateCount; ++r) {
        for (double fc : {250.0, 1000.0, 4000.0}) {
            for (float res : {0.0f, 0.5f, 1.0f}) {
                const double got = measuredGain(at(fc, res, kCentre, 0.0f, MorphLaw::HighNotchLow),
                                                kRates[r], fc, 0.25, 2.0, 0.5);
                if (!(got < 2e-4)) {
                    std::printf("FAIL line %d: SEM notch at sr %.0f fc %.0f res %.1f is %.3e "
                                "(%.1f dB) — not a null\n",
                                __LINE__, kRates[r], fc, res, got,
                                20.0 * std::log10(got + 1e-300));
                    ++g_fail;
                }
            }
        }
    }
}

// The null sits AT the cutoff, not merely somewhere nearby: the response falls monotonically into
// fc from both sides and is orders of magnitude below its own immediate neighbours. At fc=1 kHz,
// +/-5% off the notch already reads -20 dB while the notch itself reads -127 dB.
static void testHighNotchLowNullIsLocatedAtTheCutoff() {
    const double sr = 48000.0, fc = 1000.0;
    for (float res : {0.0f, 0.5f, 1.0f}) {
        const FilterSettings fs = at(fc, res, kCentre, 0.0f, MorphLaw::HighNotchLow);
        const double below[] = {0.5, 0.8, 0.95};
        double prev = 1e30;
        for (double ratio : below) {
            const double got = measuredGain(fs, sr, fc * ratio);
            CHECK(got < prev);
            prev = got;
        }
        const double atCorner = measuredGain(fs, sr, fc, 0.25, 2.0, 0.5);
        CHECK(atCorner < prev);

        prev = atCorner;
        for (double ratio : {1.05, 1.25, 2.0}) {
            const double got = measuredGain(fs, sr, fc * ratio);
            CHECK(got > prev);
            prev = got;
        }
        // Against its own immediate neighbours, so this is a null rather than a broad scoop.
        CHECK(atCorner < 1e-3 * measuredGain(fs, sr, fc * 0.95));
    }
}

// The SEM's zero is AT the notch frequency, not a broadband level sag: away from the corner the
// two taps are still an equal-power pair, so the sweep holds constant power on its legs. Measured
// deep in each tap's own passband — 50 Hz for the low tap, 20 kHz for the high tap, both far from
// a 1 kHz corner — and divided by that tap's OWN analytic response there, so what is left is the
// weight the law applied. That normalization is load-bearing, not cosmetic: at Q = 0.1 a 2-pole
// approaches its passband so slowly that the pure low tap still reads 0.896 at 50 Hz, and a raw
// reading would report a 20% "sag" that is the Q, not the morph. A LINEAR crossfade would give
// 0.5 at the centre instead of 1.0, so this tolerance discriminates equal-power from linear
// decisively rather than merely confirming a plausible shape.
static void testHighNotchLowLegsHoldConstantPowerAwayFromTheNotch() {
    const double sr = 48000.0, fc = 1000.0;
    for (float res : {0.0f, 0.5f, 1.0f}) {
        const double q = filterQFromNorm(res);
        const double lowRef = analyticMag(kLowPass, 50.0, fc, q, sr);
        const double highRef = analyticMag(kHighPass, 20000.0, fc, q, sr);
        for (int i = 0; i <= 8; ++i) {
            const float m = static_cast<float>(i) / 8.0f;
            const FilterSettings fs = at(fc, res, m, 0.0f, MorphLaw::HighNotchLow);
            const double low = measuredGain(fs, sr, 50.0) / lowRef;
            const double high = measuredGain(fs, sr, 20000.0) / highRef;
            const double power = low * low + high * high;
            if (!(std::fabs(power - 1.0) <= 0.02)) {
                std::printf("FAIL line %d: SEM morph %.3f res %.1f leg power %.6f (low %.6f, "
                            "high %.6f) — expected 1.0\n",
                            __LINE__, m, res, power, low, high);
                ++g_fail;
            }
        }
    }
}

// Continuity as a control, not just at the corner: no step between adjacent morph positions at
// any fixed frequency, under either law. A coefficient switch at the centre — the thing an enum
// over TOPOLOGIES would have forced — shows up here as a jump. Measured off the SEM's notch
// frequency, since the null itself is a legitimate near-step in the response.
static void testMorphSweepHasNoDiscontinuity() {
    const double sr = 48000.0, fc = 1000.0;
    constexpr int kSteps = 40;
    for (MorphLaw law : kBothLaws) {
        for (float res : {0.0f, 0.5f, 1.0f}) {
            for (double f : {250.0, 1000.0, 4000.0}) {
                if (f == fc && law == MorphLaw::HighNotchLow) continue;
                double prev = -1.0;
                for (int i = 0; i <= kSteps; ++i) {
                    const float m = static_cast<float>(i) / kSteps;
                    const double got = measuredGain(at(fc, res, m, 0.0f, law), sr, f);
                    if (prev >= 0.0) {
                        // Scaled by the response's own magnitude at this setting — the passband is
                        // unity and the corner is Q, so below Q=1 the passband is what a step has
                        // to be small against, not Q.
                        const double scale = std::fmax(1.0, filterQFromNorm(res));
                        // One step is 1/40 of the travel; the steepest leg moves well under a
                        // tenth of that scale over one step (measured worst case is 0.03).
                        const double jump = std::fabs(got - prev) / scale;
                        if (!(jump < 0.1)) {
                            std::printf("FAIL line %d: %s morph %.4f res %.1f at %.0f Hz jumps "
                                        "%.4f\n",
                                        __LINE__, lawName(law), m, res, f, jump);
                            ++g_fail;
                        }
                    }
                    prev = got;
                }
            }
        }
    }
}

// The law selects a MIX, computed once per prepare(); it must not reach the coefficient solve at
// all. Asserted bit-exactly rather than by tolerance — the cutoff, the damping term, and the
// zero-delay-loop solution are the same floats under either law, so no cutoff/Q/rate behaviour
// can differ between them by construction.
static void testMorphLawDoesNotDisturbTheCoefficients() {
    for (int r = 0; r < kRateCount; ++r) {
        for (int ci = 0; ci <= 8; ++ci) {
            for (float res : {0.0f, 0.5f, 1.0f}) {
                for (int mi = 0; mi <= 4; ++mi) {
                    VoiceFilter band, sem;
                    const float m = mi / 4.0f;
                    band.prepare({ci / 8.0f, res, m, 0.5f, MorphLaw::HighBandLow}, kRates[r]);
                    sem.prepare({ci / 8.0f, res, m, 0.5f, MorphLaw::HighNotchLow}, kRates[r]);
                    const SvfCoeffs& a = band.coeffs();
                    const SvfCoeffs& b = sem.coeffs();
                    CHECK(a.g == b.g && a.k == b.k);
                    CHECK(a.a1 == b.a1 && a.a2 == b.a2 && a.a3 == b.a3);
                }
            }
        }
    }
}

// The default is the reviewed-and-measured law, not the SEM leg. The editor and any persisted-
// state codec read this default, so a preset saved before the selector existed must still sound
// exactly as it did — asserted on the folded mix, which is the only thing the kernel sees.
static void testFilterSettingsDefaultsToTheHighBandLowLaw() {
    CHECK(FilterSettings{}.morphLaw == MorphLaw::HighBandLow);

    VoiceFilter defaulted, explicitLaw;
    defaulted.prepare({0.5f, 0.5f, kCentre, 0.0f}, 48000.0);
    explicitLaw.prepare({0.5f, 0.5f, kCentre, 0.0f, MorphLaw::HighBandLow}, 48000.0);
    CHECK(defaulted.mix().m0 == explicitLaw.mix().m0);
    CHECK(defaulted.mix().m1 == explicitLaw.mix().m1);
    CHECK(defaulted.mix().m2 == explicitLaw.mix().m2);
}

// ---------------------------------------------------------------------------
// Drive
// ---------------------------------------------------------------------------

// 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. 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 checkDriveZeroBitIdentity(float morph, MorphLaw law) {
    VoiceFilter f;
    f.prepare(at(1000.0, 1.0f, morph, 0.0f, law), 48000.0);
    f.reset();
    LinearKernelRef ref0{f.coeffs(), f.mix()};
    LinearKernelRef ref1{f.coeffs(), f.mix()};

    unsigned rng0 = 0x13579bdfu;
    for (int i = 0; i < 4096; ++i) {
        const float x = nextNoise(rng0);
        CHECK(f.process(0, x) == ref0.step(x));
    }

    // 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);
    }
}

static void testDriveZeroIsBitIdenticalToTheLinearKernel() {
    for (MorphLaw law : kBothLaws) {
        for (float morph : {kHighPass, kBandPass, kLowPass}) checkDriveZeroBitIdentity(morph, law);
    }
}

// The complaint the rewrite answers: resonance must not track how hard the sample hits the
// filter unless the user asked for it. At drive 0 the response is identical over a 1000:1 level
// range; the tap this replaced moved by 14% over the same span.
static void testDriveZeroResponseIsLevelInvariant() {
    const double sr = 48000.0, fc = 1000.0;
    for (float morph : {kHighPass, kBandPass, kLowPass}) {
        const double q = filterQFromNorm(1.0f);
        const double want = analyticMag(morph, fc, fc, q, sr);
        for (double amp : {0.001, 0.01, 0.1, 1.0}) {
            const double got = measuredGain(at(fc, 1.0f, morph), sr, fc, amp);
            if (!(std::fabs(got / want - 1.0) <= kAgreement)) {
                std::printf("FAIL line %d: morph %.1f amp %g gain %.6f vs analytic %.6f "
                            "(%.3f%%)\n",
                            __LINE__, morph, amp, got, want, (got / want - 1.0) * 100.0);
                ++g_fail;
            }
        }
    }
}

// Drive is bounded by construction, not by tuning: softLimit is a contraction, so the state
// update can only ever shrink the state and the filter cannot gain energy from it. This sweeps
// the corners that would expose a tuned margin instead.
static void testFullDriveStaysBoundedAtEveryCutoffResonanceAndRate() {
    unsigned rng = 0x2468aceu;
    auto noise = [&rng]() {
        rng = rng * 1664525u + 1013904223u;
        return static_cast<float>(static_cast<int>(rng >> 9) - (1 << 22)) /
               static_cast<float>(1 << 22);
    };

    for (int r = 0; r < kRateCount; ++r) {
        const double sr = kRates[r];
        for (MorphLaw law : kBothLaws) {
            for (int ci = 0; ci <= 8; ++ci) {
                for (int mi = 0; mi <= 4; ++mi) {
                    for (float res : {0.0f, 0.5f, 1.0f}) {
                        VoiceFilter f;
                        f.prepare({ci / 8.0f, res, mi / 4.0f, 1.0f, law}, sr);
                        f.reset();
                        for (int i = 0; i < 4000; ++i) {
                            const float y = f.process(0, noise());
                            if (!std::isfinite(y) || std::fabs(y) > 8.0f) {
                                std::printf("FAIL line %d: %s sr=%.0f cutoff=%.2f morph=%.2f "
                                            "res=%.1f full drive produced %g\n",
                                            __LINE__, lawName(law), sr, ci / 8.0, mi / 4.0, res, y);
                                ++g_fail;
                                return;
                            }
                        }
                    }
                }
            }
        }
    }
}

// Full drive at full resonance with no input must still go quiet. A nonlinearity in the loop is
// exactly where a self-oscillator would hide, and softLimit's sub-unit slope is what forbids it.
static void testFullDriveDoesNotSelfOscillate() {
    for (int r = 0; r < kRateCount; ++r) {
        const double sr = kRates[r];
        for (MorphLaw law : kBothLaws) {
            for (float morph : {kHighPass, kBandPass, kLowPass}) {
                VoiceFilter f;
                f.prepare(at(1000.0, 1.0f, morph, 1.0f, law), sr);
                f.reset();
                const int excite = static_cast<int>(sr * 0.01);
                for (int i = 0; i < excite; ++i) {
                    f.process(0, static_cast<float>(std::sin(2.0 * kPi * 1000.0 * i / sr)));
                }
                for (int i = 0; i < static_cast<int>(sr * 0.5); ++i) f.process(0, 0.0f);
                CHECK(f.isSilent());
            }
        }
    }
}

// Drive has to actually do something at the top of its travel, and do it monotonically — the
// brief's "extreme, not politely warm". Measured at the corner, where the resonance state is
// what the limiter sees.
static void testDriveCompressesTheResonantPeakMonotonically() {
    const double sr = 48000.0, fc = 1000.0;
    double prev = 1e30;
    for (int i = 0; i <= 8; ++i) {
        const double got = measuredGain(at(fc, 1.0f, kLowPass, i / 8.0f), sr, fc, 1.0);
        CHECK(got < prev);
        prev = got;
    }
    // Full drive against no drive: a large, unmistakable reduction of the resonant peak.
    CHECK(prev < 0.5 * filterQFromNorm(1.0f));

    // And the passband is left alone at every drive setting — drive colours the resonance, it
    // is not a distortion box in series with the signal.
    for (int i = 0; i <= 4; ++i) {
        CHECK_NEAR(measuredGain(at(fc, 1.0f, kLowPass, i / 4.0f), sr, 100.0, 1.0), 1.0, 0.05);
    }
}

static void testSoftLimitIsOddMonotoneBoundedAndExactAtZeroDepth() {
    for (double x : {-3.0, -0.5, 0.0, 1e-9, 0.25, 7.0}) {
        // Depth 0 is the identity by algebra, so drive 0 needs no special case on the hot path.
        CHECK(softLimit(static_cast<float>(x), 0.0f) == static_cast<float>(x));
    }
    CHECK_NEAR(softLimit(1.5f, 2.0f), -softLimit(-1.5f, 2.0f), 1e-9);

    for (float depth : {0.5f, 4.0f, 64.0f}) {
        // The two properties the stability argument rests on, over the whole excursion range a
        // resonating state can reach. Monotonicity is NOT asserted here: far past the knee the
        // curve is asymptotically flat, so the true increment between adjacent samples falls
        // below float epsilon and rounding can walk it backwards by an ulp.
        for (int i = -400; i <= 400; ++i) {
            const float x = static_cast<float>(i) * 0.05f;
            const float y = softLimit(x, depth);
            CHECK(std::fabs(y) <= std::fabs(x));        // a contraction — the stability argument
            CHECK(std::fabs(y) < 1.0f / depth + 1e-6f); // bounded by the knee
        }
        // Strictly increasing across the knee, which is where the shaping actually happens.
        const float knee = 1.0f / depth;
        float prev = -1e30f;
        for (int i = -20; i <= 20; ++i) {
            const float y = softLimit(static_cast<float>(i) * 0.1f * knee, depth);
            CHECK(y > prev);
            prev = y;
        }
    }
}

// ---------------------------------------------------------------------------
// Sample-rate invariance
// ---------------------------------------------------------------------------

// The rate must enter only through g = tan(pi*fc/sr), so the response at a given cutoff and Q is
// the same filter at every rate. The retired feedback tap made this false: it closed the loop
// once per SAMPLE, so emphasis ran 5.02 at 48k against 8.52 at 192k.
static void testResponseIsRateInvariantAtEveryMorph() {
    for (float morph : {kHighPass, kBandPass, kLowPass}) {
        for (float res : {0.2f, 0.5f, 1.0f}) {
            const double q = filterQFromNorm(res);
            for (double fc : {250.0, 1000.0, 4000.0}) {
                for (int r = 0; r < kRateCount; ++r) {
                    const double got = measuredGain(at(fc, res, morph), kRates[r], fc);
                    const double want = analyticMag(morph, fc, fc, q, kRates[r]);
                    if (!(std::fabs(got / want - 1.0) <= kAgreement)) {
                        std::printf("FAIL line %d: morph %.1f res %.1f fc %.0f at %.0f Hz: %.6f "
                                    "vs analytic %.6f (%.3f%%)\n",
                                    __LINE__, morph, res, fc, kRates[r], got, want,
                                    (got / want - 1.0) * 100.0);
                        ++g_fail;
                    }
                }
            }
        }
    }
}

// The conditioning corner: fc/sr ~ 1e-4. Float32 Direct Form I encoded pole proximity in
// a1 -> -2, a2 -> +1 and cancelled them every sample, costing ~17 bits and putting the measured
// peak 15% LOW at 20 Hz / 192 kHz. TPT encodes the same proximity in a1's small deviation from
// 1, which float resolves; this pins that the defect is gone at every rate.
static void testLowCutoffHighRateCornerHoldsTheAnalyticPeak() {
    const double q = filterQFromNorm(1.0f);
    // A 2-pole low-pass peaks at W = sqrt(1 - 1/(2Q^2)), where |H| = Q / sqrt(1 - 1/(4Q^2)).
    const double wPeak = std::sqrt(1.0 - 1.0 / (2.0 * q * q));
    const double want = q / std::sqrt(1.0 - 1.0 / (4.0 * q * q));
    CHECK_NEAR(want, 10.012516, 1e-5);  // the figure the measurement pass quoted

    for (int r = 0; r < kRateCount; ++r) {
        const double sr = kRates[r];
        const double fPeak = sr / kPi * std::atan(wPeak * std::tan(kPi * 20.0 / sr));
        // Q=10 at 20 Hz rings for ~0.16 s, so the settle window has to be seconds, not samples.
        const double got = measuredGain(at(20.0, 1.0f, kLowPass), sr, fPeak, 0.25, 3.0, 1.0);
        if (!(std::fabs(got / want - 1.0) <= kAgreement)) {
            std::printf("FAIL line %d: 20 Hz peak at %.0f Hz is %.6f vs analytic %.6f (%.3f%%)\n",
                        __LINE__, sr, got, want, (got / want - 1.0) * 100.0);
            ++g_fail;
        }
    }
}

int main() {
    testSvfCoefficientsMatchPinnedValues();
    testTheSampleRateEntersOnlyThroughG();
    testNonPositiveRatePrepareClearsStaleStateAndReportsSilent();
    testNonPositiveRatePassesSignalThroughAtEveryMorph();

    testMorphEndpointsMatchTheAnalyticTwoPoleTargets();
    testCornerMagnitudeIsFlatAtQAcrossTheHighBandLowSweep();
    testHighNotchLowCentreIsATrueNullAtTheCorner();
    testHighNotchLowNullIsLocatedAtTheCutoff();
    testHighNotchLowLegsHoldConstantPowerAwayFromTheNotch();
    testMorphSweepHasNoDiscontinuity();
    testMorphLawDoesNotDisturbTheCoefficients();
    testFilterSettingsDefaultsToTheHighBandLowLaw();

    testDriveZeroIsBitIdenticalToTheLinearKernel();
    testDriveZeroResponseIsLevelInvariant();
    testFullDriveStaysBoundedAtEveryCutoffResonanceAndRate();
    testFullDriveDoesNotSelfOscillate();
    testDriveCompressesTheResonantPeakMonotonically();
    testSoftLimitIsOddMonotoneBoundedAndExactAtZeroDepth();

    testResponseIsRateInvariantAtEveryMorph();
    testLowCutoffHighRateCornerHoldsTheAnalyticPeak();

    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;
}