// 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::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;

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

    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. The settle and measure windows are wall-clock, not sample
// counts, so every rate integrates the same amount of signal.
static double measuredRms(FilterMode mode, float cutoffNorm, float resNorm, double freqHz,
                          double sr, double amp = 1.0) {
    VoiceFilter f;
    f.prepare({mode, cutoffNorm, resNorm}, sr);
    f.reset();

    const int settle = static_cast<int>(sr * 0.15);
    const int measure = static_cast<int>(sr * 0.10);
    double sumSq = 0.0;
    for (int i = 0; i < settle + measure; ++i) {
        const float x = static_cast<float>(amp * 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);
    }
}

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

// The rates the invariance claim is made over. 88.2k is deliberately included: it is the rate
// whose calibrated feedback delay lands between two whole taps, so it is the one the
// interpolating read has to earn.
static const double kRates[] = {44100.0, 48000.0, 88200.0, 96000.0, 192000.0};
static constexpr int kRateCount = 5;
static constexpr int kRef48k = 1;  // index of the reference rate within kRates

// Resonant emphasis: level at the cutoff over the passband level. Measured at the requested
// cutoff rather than at the scanned peak so no frequency-grid quantization leaks into the
// comparison. The passband reference is the same frequency at every rate, or the ratio would
// compare a different measurement at each rate -- and it must stay well clear of the LOWEST
// Nyquist tested, since a high-pass reference near 44.1k's band edge measures the bilinear
// warping rather than the resonance.
static double emphasisAtCutoff(FilterMode mode, double fcHz, float resNorm, double sr) {
    const float cn = filterNormFromCutoffHz(static_cast<float>(fcHz));
    const double refHz = (mode == FilterMode::LowPass) ? fcHz / 8.0 : fcHz * 8.0;
    return measuredRms(mode, cn, resNorm, fcHz, sr, 0.25) /
           measuredRms(mode, cn, resNorm, refHz, sr, 0.25);
}

// The feedback loop's contribution alone: the measured closed-loop level at a frequency over the
// level the bare coefficients predict there. Dividing the coefficient response out removes the
// bilinear discretization difference between rates -- which is real, correct, and not something
// a feedback fix can or should touch -- leaving exactly the loop under audit. In low-pass mode
// there is no loop, so this is identically 1 at every rate.
static double feedbackContribution(FilterMode mode, double fcHz, float resNorm, double sr) {
    const float cn = filterNormFromCutoffHz(static_cast<float>(fcHz));
    VoiceFilter f;
    f.prepare({mode, cn, resNorm}, sr);
    const double openLoopRms = magnitudeAt(f.coeffs(), fcHz, sr) * 0.25 / std::sqrt(2.0);
    return measuredRms(mode, cn, resNorm, fcHz, sr, 0.25) / openLoopRms;
}

// Where the response actually peaks, as a multiple of the requested cutoff.
static double peakOverCutoff(FilterMode mode, double fcHz, float resNorm, double sr) {
    const float cn = filterNormFromCutoffHz(static_cast<float>(fcHz));
    double peak = 0.0, peakF = 0.0;
    for (int i = 0; i <= 12; ++i) {
        const double f = fcHz * std::pow(2.0, -0.5 + i / 12.0);
        const double r = measuredRms(mode, cn, resNorm, f, sr, 0.25);
        if (r > peak) { peak = r; peakF = f; }
    }
    return peakF / fcHz;
}

// The defect these pin: the high-pass feedback loop closes once per sample, so while its tap was
// the immediately previous output the loop's phase at the cutoff -- and with it the resonant
// emphasis -- scaled with the sample rate. Against that one-sample tap, emphasisAtCutoff for
// fc=1 kHz, res=1.0 measured 5.46 at 48k rising monotonically to 6.04 at 192k (10.5%), and
// feedbackContribution for fc=4 kHz, res=1.0 ran 0.443 at 48k against 0.506 at 192k (14.4%).
// Both now sit inside the bounds below.
//
// The two tolerances split on the reference rate, and the split is load-bearing rather than
// convenient. At or above 48k the calibrated interval is at least one sample, so the tap
// reproduces it and only the bilinear discretization difference remains. Below it -- 44.1k --
// one sample is ALREADY longer than the interval, so the delay cannot be shortened to match
// without a sub-sample delay the loop cannot contain; 44.1k is left exactly where it has always
// been, which is up to 6% off 48k at the top of the cutoff range.
static constexpr double kAtOrAboveReferenceTolerance = 0.02;
static constexpr double kBelowReferenceTolerance = 0.08;

static void checkInvariant(const char* what, FilterMode mode, double fcHz, float resNorm,
                           double (*measure)(FilterMode, double, float, double)) {
    const double reference = measure(mode, fcHz, resNorm, kRates[kRef48k]);
    for (int r = 0; r < kRateCount; ++r) {
        const double v = measure(mode, fcHz, resNorm, kRates[r]);
        const double deviation = std::fabs(v - reference) / reference;
        const double tolerance = kRates[r] >= kRates[kRef48k] ? kAtOrAboveReferenceTolerance
                                                              : kBelowReferenceTolerance;
        if (!(deviation <= tolerance)) {
            std::printf("FAIL line %d: %s %s fc=%.0f res=%.2f at %.0f Hz: %.5f vs 48k %.5f "
                        "(%.2f%% > %.2f%%)\n",
                        __LINE__, what, mode == FilterMode::LowPass ? "LP" : "HP", fcHz, resNorm,
                        kRates[r], v, reference, deviation * 100.0, tolerance * 100.0);
            ++g_fail;
        }
    }
}

// End-to-end: the emphasis a listener hears, coefficients and feedback together. Held to cutoffs
// whose passband reference (8x the cutoff) stays well below 44.1k's band edge -- higher cutoffs
// are covered by the isolated test below, which does not need a passband reference at all.
static void testHighPassResonanceIsRateInvariant() {
    for (float res : {0.2f, 0.5f, 1.0f}) {
        checkInvariant("emphasis", FilterMode::HighPass, 250.0, res, emphasisAtCutoff);
        checkInvariant("emphasis", FilterMode::HighPass, 1000.0, res, emphasisAtCutoff);
    }
}

// The low-pass has no feedback path, so it was already invariant. Pinning it is the control: it
// proves the measurement detects what it claims to, and it keeps a future feedback path on the
// low-pass from acquiring the same defect unnoticed.
static void testLowPassResonanceIsRateInvariant() {
    for (float res : {0.2f, 0.5f, 1.0f}) {
        checkInvariant("emphasis", FilterMode::LowPass, 250.0, res, emphasisAtCutoff);
        checkInvariant("emphasis", FilterMode::LowPass, 1000.0, res, emphasisAtCutoff);
        checkInvariant("emphasis", FilterMode::LowPass, 4000.0, res, emphasisAtCutoff);
    }
}

// The precise form of the same claim, with the discretization difference divided out, so it also
// holds at the top of the cutoff range where a passband reference cannot sit clear of 44.1k's
// band edge.
static void testFeedbackLoopContributionIsRateInvariant() {
    for (float res : {0.2f, 0.5f, 1.0f}) {
        for (double fc : {250.0, 1000.0, 4000.0}) {
            checkInvariant("loop", FilterMode::HighPass, fc, res, feedbackContribution);
            checkInvariant("loop", FilterMode::LowPass, fc, res, feedbackContribution);
        }
    }
}

static void testResonantPeakTracksCutoffAtEveryRate() {
    for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) {
        for (double fc : {250.0, 1000.0, 4000.0}) {
            for (int r = 0; r < kRateCount; ++r) {
                // At full resonance there is a real peak to find; a quarter octave either side
                // of the requested cutoff is the same window the 48k-only test uses.
                const double ratio = peakOverCutoff(mode, fc, 1.0f, kRates[r]);
                if (!(ratio > 1.0 / 1.19 && ratio < 1.19)) {
                    std::printf("FAIL line %d: %s peak at %.3f x fc (fc=%.0f, sr=%.0f)\n",
                                __LINE__, mode == FilterMode::LowPass ? "LP" : "HP", ratio, fc,
                                kRates[r]);
                    ++g_fail;
                }
            }
        }
    }
}

// 48k is the rate the feedback constants were voiced at, and the rate Daniel's ear judgments
// were made against, so making the other rates match it must not move it. These literals were
// captured from the build BEFORE the fixed-time feedback tap landed; the tap resolves to
// exactly one sample at 48k, so they must reproduce bit-for-bit rather than merely closely.
static void testFortyEightKilohertzBehaviorIsUnchanged() {
    struct Pin {
        FilterMode mode;
        double y1, y7, y31, y127, energy, sineRms;
    };
    const Pin pins[2] = {
        {FilterMode::LowPass, 0.016871979, 0.098936319, -0.084338546, -0.044524558, 0.652648822,
         1.767755710},
        {FilterMode::HighPass, -0.184770823, -0.082661532, 0.057580549, -0.008336116, 1.427662234,
         0.895141269},
    };

    for (const Pin& p : pins) {
        VoiceFilter f;
        f.prepare({p.mode, filterNormFromCutoffHz(1000.0f), 1.0f}, 48000.0);
        f.reset();
        double energy = 0.0;
        for (int i = 0; i < 4096; ++i) {
            const float y = f.process(0, i == 0 ? 1.0f : 0.0f);
            energy += static_cast<double>(y) * y;
            if (i == 1) CHECK_NEAR(y, p.y1, 1e-7);
            if (i == 7) CHECK_NEAR(y, p.y7, 1e-7);
            if (i == 31) CHECK_NEAR(y, p.y31, 1e-7);
            if (i == 127) CHECK_NEAR(y, p.y127, 1e-7);
        }
        CHECK_NEAR(energy, p.energy, 1e-7);

        VoiceFilter g;
        g.prepare({p.mode, filterNormFromCutoffHz(1000.0f), 1.0f}, 48000.0);
        g.reset();
        double sumSq = 0.0;
        for (int i = 0; i < 28800; ++i) {
            const float x = static_cast<float>(0.25 * std::sin(2.0 * kPi * 1000.0 * i / 48000.0));
            const float y = g.process(0, x);
            if (i >= 14400) sumSq += static_cast<double>(y) * y;
        }
        CHECK_NEAR(std::sqrt(sumSq / 14400.0), p.sineRms, 1e-7);
    }
}

// The tap is a fixed INTERVAL, so the sample offset it resolves to scales with the rate. Read
// out of the filter's behavior, not its internals: run an impulse through the high-pass and
// alongside it the bare difference equation on the SAME coefficients with no feedback at all.
// The tap reads y[n-D], and every earlier history slot is zero, so the first sample at which the
// two can possibly diverge is exactly D. Against the pre-fix one-sample tap this reports 1 at
// every rate; it must now report 1, 1, 1, 2, 4.
static void testFeedbackTapOffsetScalesWithSampleRate() {
    const int expected[kRateCount] = {1, 1, 1, 2, 4};
    for (int r = 0; r < kRateCount; ++r) {
        VoiceFilter f;
        f.prepare({FilterMode::HighPass, filterNormFromCutoffHz(1000.0f), 1.0f}, kRates[r]);
        f.reset();
        const BiquadCoeffs c = f.coeffs();

        float x1 = 0.0f, x2 = 0.0f, y1 = 0.0f, y2 = 0.0f;
        int firstDivergence = -1;
        for (int i = 0; i < 64 && firstDivergence < 0; ++i) {
            const float x = (i == 0) ? 1.0f : 0.0f;
            const float actual = f.process(0, x);
            const float noFeedback = c.b0 * x + c.b1 * x1 + c.b2 * x2 - c.a1 * y1 - c.a2 * y2;
            x2 = x1;
            x1 = x;
            y2 = y1;
            y1 = noFeedback;
            if (actual != noFeedback) firstDivergence = i;
        }
        if (firstDivergence != expected[r]) {
            std::printf("FAIL line %d: sr=%.0f feedback first reaches the output at sample %d, "
                        "expected %d\n",
                        __LINE__, kRates[r], firstDivergence, expected[r]);
            ++g_fail;
        }
    }
}

// The floor is load-bearing, not defensive: below 48k one sample is ALREADY longer than the
// calibrated interval, so the offset cannot shrink to match without a sub-sample delay the loop
// cannot contain -- it would be algebraic and uncomputable. A rate at or below the reference
// therefore keeps the firmware's single tap, and a non-positive rate lands on the same floor
// rather than on an invented rate.
static void testFeedbackTapNeverFallsBelowOneSample() {
    for (double sr : {-48000.0, 0.0, 1000.0, 22050.0, 44100.0, 48000.0}) {
        VoiceFilter f;
        f.prepare({FilterMode::HighPass, filterNormFromCutoffHz(1000.0f), 1.0f}, sr);
        f.reset();
        for (int i = 0; i < 512; ++i) CHECK(std::isfinite(f.process(0, i == 0 ? 1.0f : 0.0f)));
    }
}

// ---------------------------------------------------------------------------
// 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 (int r = 0; r < kRateCount; ++r) {
        const double sr = kRates[r];
        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();
                    // A fixed WALL-CLOCK sweep: the same cutoff travel per second at every rate,
                    // so the per-sample coefficient step gets no gentler as the rate rises.
                    const int n = static_cast<int>(sr * 0.25);
                    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
                    }
                }
            }
        }
    }
}

// The decay to the floor is a fixed WALL-CLOCK time (~0.21 s at these settings), not a fixed
// sample count -- so the budget has to scale with the rate. A fixed 20000-sample budget is itself
// a rate assumption: it is ample at 48k and expires mid-decay at 96k and above.
static void testStateFlushesToZeroWithoutStallingInDenormals() {
    for (int r = 0; r < kRateCount; ++r) {
        const double sr = kRates[r];
        const int budget = static_cast<int>(sr * 0.5);
        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.
            const int excite = static_cast<int>(sr * 0.01);
            for (int i = 0; i < excite; ++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 < budget; ++i) {
                f.process(0, 0.0f);
                const VoiceFilter::State& s = f.state(0);
                bool subnormal = false;
                for (float v : {s.x1, s.x2, s.y1, s.y2}) {
                    if (v != 0.0f && std::fabs(v) < FLT_MIN) subnormal = true;
                }
                for (float v : s.fb) {
                    if (v != 0.0f && std::fabs(v) < FLT_MIN) subnormal = true;
                }
                if (subnormal) ++subnormalSamples;
                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. The feedback
            // tap line holds copies of the flushed y, so it drains behind it rather than feeding
            // subnormals back into the loop.
            CHECK(subnormalSamples <= 2);
            CHECK(silentAt >= 0);
            CHECK(silentAt < budget);
            // 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());
        }
    }
}

// A high-pass has an exact DC null (b1 == -2*b0 bit-exactly), so under sustained DC the
// recursive y decays to zero while x1/x2 sit pinned at the DC level -- the case the zero-input
// test above cannot see, since there x1/x2 are zero anyway. A flush that clears x1/x2 along
// with y1/y2 discards that pinned history; the next sample then recomputes a full-amplitude
// step from b0*in alone, which re-rings and repeats forever (a click train). This must fail
// against a flush that also clears x1/x2.
// Run at full resonance as well as none: at res=0 the feedback share is zero and the tap line is
// inert, so that case alone would never notice the tap line failing to drain behind a flush.
static void testHighPassSustainedDCDoesNotReRing() {
    for (int r = 0; r < kRateCount; ++r) {
        for (float res : {0.0f, 1.0f}) {
            const double sr = kRates[r];
            VoiceFilter f;
            f.prepare({FilterMode::HighPass, filterNormFromCutoffHz(1000.0f), res}, sr);
            f.reset();

            const int settle = static_cast<int>(sr * 0.05);
            float worstAfterSettle = 0.0f;
            for (int i = 0; i < static_cast<int>(sr * 0.5); ++i) {
                const float y = f.process(0, 1.0f);
                if (i >= settle) {
                    const float a = std::fabs(y);
                    if (a > worstAfterSettle) worstAfterSettle = a;
                }
            }
            // A correct flush leaves the settled output pinned near zero. The click train this
            // regresses against recurs every ~4760 samples at 48k at a magnitude around 0.6 --
            // nowhere near this tolerance.
            CHECK(worstAfterSettle < 1e-3f);
        }
    }
}

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

    // A DC step through a highpass should settle to (and STAY AT) zero. Sampling only the
    // final value is not enough to prove that: a resonator swings through zero twice a cycle,
    // so a single late sample can land near zero while the envelope is still ringing well
    // above it elsewhere in the same run -- track the worst case over the settled region.
    VoiceFilter hp;
    hp.prepare({FilterMode::HighPass, filterNormFromCutoffHz(1000.0f), 0.0f}, sr);
    hp.reset();
    const int settle = 200;
    float worstAfterSettle = 0.0f;
    for (int i = 0; i < 48000; ++i) {
        y = hp.process(0, 1.0f);
        if (i >= settle) {
            const float a = std::fabs(y);
            if (a > worstAfterSettle) worstAfterSettle = a;
        }
    }
    CHECK(worstAfterSettle < 1e-3f);  // fully rejected by a highpass, not just at one instant
}

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 testFeedbackSaturationIsContinuousWithGentleLinearTail() {
    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);
    // Past the threshold the curve continues on a 0.1 slope rather than hard-clipping -- it is
    // NOT bounded, so this pins the linear continuation's shallow slope, not a ceiling.
    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();
    testHighPassResonanceIsRateInvariant();
    testLowPassResonanceIsRateInvariant();
    testFeedbackLoopContributionIsRateInvariant();
    testResonantPeakTracksCutoffAtEveryRate();
    testFortyEightKilohertzBehaviorIsUnchanged();
    testFeedbackTapOffsetScalesWithSampleRate();
    testFeedbackTapNeverFallsBelowOneSample();
    testFullRangeCutoffSweepAtAudioRateStaysBounded();
    testStateFlushesToZeroWithoutStallingInDenormals();
    testHighPassSustainedDCDoesNotReRing();
    testImpulseResponseMatchesDifferenceEquation();
    testLowpassStepSettlesToUnity();
    testResetClearsHistoryButPrepareKeepsIt();
    testChannelStateIsIndependent();
    testFeedbackSaturationIsContinuousWithGentleLinearTail();

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