// 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.
static double measuredRms(FilterMode mode, float cutoffNorm, float resNorm, double freqHz,
                          double sr) {
    VoiceFilter f;
    f.prepare({mode, cutoffNorm, resNorm}, sr);
    f.reset();

    const int settle = 24000, measure = 24000;
    double sumSq = 0.0;
    for (int i = 0; i < settle + measure; ++i) {
        const float x = static_cast<float>(std::sin(2.0 * kPi * freqHz * i / sr));
        const float y = f.process(0, x);
        if (i >= settle) sumSq += static_cast<double>(y) * y;
    }
    return std::sqrt(sumSq / measure);
}

static void testMeasuredResponsePeaksAtCutoffInBothModes() {
    const double sr = 48000.0;
    const float cutoffNorm = filterNormFromCutoffHz(1000.0f);

    for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) {
        double peakRms = 0.0, peakFreq = 0.0;
        for (int i = 0; i <= 40; ++i) {
            const double f = 100.0 * std::pow(100.0, static_cast<double>(i) / 40.0);
            const double r = measuredRms(mode, cutoffNorm, 1.0f, f, sr);
            if (r > peakRms) { peakRms = r; peakFreq = f; }
        }
        CHECK(peakFreq > 1000.0 / 1.3 && peakFreq < 1000.0 * 1.3);

        const double passband =
            measuredRms(mode, cutoffNorm, 1.0f, mode == FilterMode::LowPass ? 100.0 : 10000.0, sr);
        CHECK(peakRms / passband > 3.0);

        // Same measurement at the bottom of the resonance control shows no such emphasis.
        const double flatAtCutoff = measuredRms(mode, cutoffNorm, 0.0f, 1000.0, sr);
        const double flatPassband =
            measuredRms(mode, cutoffNorm, 0.0f, mode == FilterMode::LowPass ? 100.0 : 10000.0, sr);
        CHECK(flatAtCutoff / flatPassband < 1.0);
    }
}

// ---------------------------------------------------------------------------
// Stability
// ---------------------------------------------------------------------------

static void testFullRangeCutoffSweepAtAudioRateStaysBounded() {
    // Deterministic pseudo-noise; a fixed sine would miss the resonant frequency on most steps.
    unsigned rng = 0x13579bdfu;
    auto noise = [&rng]() {
        rng = rng * 1664525u + 1013904223u;
        return static_cast<float>(static_cast<int>(rng >> 9) - (1 << 22)) / static_cast<float>(1 << 22);
    };

    for (double sr : {44100.0, 48000.0, 96000.0}) {
        for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) {
            for (float res : {0.0f, 0.5f, 1.0f}) {
                for (int direction = 0; direction < 2; ++direction) {
                    VoiceFilter f;
                    f.reset();
                    const int n = 48000;
                    for (int i = 0; i < n; ++i) {
                        const float t = static_cast<float>(i) / static_cast<float>(n - 1);
                        // Per-sample coefficient update across the whole cutoff travel.
                        f.prepare({mode, direction == 0 ? t : 1.0f - t, res}, sr);
                        const float y = f.process(0, noise());
                        CHECK(std::isfinite(y));
                        CHECK(std::fabs(y) < 100.0f);
                        if (!std::isfinite(y)) return;  // stop before the log floods
                    }
                }
            }
        }
    }
}

static void testStateFlushesToZeroWithoutStallingInDenormals() {
    const double sr = 48000.0;
    for (FilterMode mode : {FilterMode::LowPass, FilterMode::HighPass}) {
        VoiceFilter f;
        f.prepare({mode, filterNormFromCutoffHz(1000.0f), 1.0f}, sr);
        f.reset();

        // Excite, then hard-cut to silence the way a released voice does.
        for (int i = 0; i < 480; ++i) {
            f.process(0, 0.5f * static_cast<float>(std::sin(2.0 * kPi * 1000.0 * i / sr)));
        }

        int subnormalSamples = 0;
        int silentAt = -1;
        for (int i = 0; i < 20000; ++i) {
            f.process(0, 0.0f);
            const VoiceFilter::State& s = f.state(0);
            const float vals[5] = {s.x1, s.x2, s.y1, s.y2, s.fb};
            for (float v : vals) {
                if (v != 0.0f && std::fabs(v) < FLT_MIN) { ++subnormalSamples; break; }
            }
            if (silentAt < 0 && f.isSilent()) silentAt = i;
        }
        // Without the flush the state grinds down through the subnormal range for thousands
        // of samples; a stray sample or two at a zero crossing is not a stall.
        CHECK(subnormalSamples <= 2);
        CHECK(silentAt >= 0);
        CHECK(silentAt < 20000);
        // And it stays silent — a flush that perturbs the feedback loop would re-excite it.
        for (int i = 0; i < 1000; ++i) CHECK(f.process(0, 0.0f) == 0.0f);
        CHECK(f.isSilent());
    }
}

// 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.
static void testHighPassSustainedDCDoesNotReRing() {
    const double sr = 48000.0;
    VoiceFilter f;
    f.prepare({FilterMode::HighPass, filterNormFromCutoffHz(1000.0f), 0.0f}, sr);
    f.reset();

    const int settle = 1000;
    float worstAfterSettle = 0.0f;
    for (int i = 0; i < 20000; ++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 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();
    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;
}