// limiter.cpp — see limiter.h.

#include "core/instrument/engine/limiter.h"

#include <algorithm>
#include <cmath>

namespace reasampler::instrument::engine {

namespace {

constexpr int kProtoLen = kLimiterOversample * kLimiterOsTaps + 1;  // 33: odd, so phase 0 is exact

double sincPi(double x) {
    if (x == 0.0) return 1.0;
    const double a = 3.14159265358979323846 * x;
    return std::sin(a) / a;
}

}  // namespace

double limiterCeilingLinear() { return std::pow(10.0, kLimiterCeilingDbTp / 20.0); }

int limiterLookaheadSamples(double sampleRate) {
    if (!(sampleRate > 0.0)) return 0;
    const int n = static_cast<int>(kLimiterLookaheadSeconds * sampleRate + 0.5);
    // One sample above the detector's group delay is the floor: the smoothing window must have
    // at least one entry of its own for the no-overshoot bound to say anything.
    return n > kLimiterOsDelay ? n : kLimiterOsDelay + 1;
}

void Limiter::prepare(double sampleRate) {
    latency_ = limiterLookaheadSamples(sampleRate);
    if (latency_ <= 0) latency_ = kLimiterOsDelay + 1;
    window_ = latency_ - kLimiterOsDelay + 1;
    ceiling_ = static_cast<float>(limiterCeilingLinear());
    const double rate = sampleRate > 0.0 ? sampleRate : 48000.0;
    releaseCoeff_ = static_cast<float>(1.0 - std::exp(-1.0 / (kLimiterReleaseSeconds * rate)));
    mixStep_ = static_cast<float>(1.0 / (kLimiterCrossfadeSeconds * rate));

    // Windowed-sinc polyphase interpolator, built here because it costs transcendentals.
    // Phase 0's taps all land on sinc zeros except the centre, so it is an exact delay and is
    // read straight out of the history instead of being convolved.
    for (int p = 0; p < kLimiterOversample; ++p) {
        for (int k = 0; k < kLimiterOsTaps; ++k) {
            const int i = kLimiterOversample * k + p;
            const double centred = static_cast<double>(i) - (kProtoLen - 1) / 2.0;
            const double hann =
                0.5 - 0.5 * std::cos(2.0 * 3.14159265358979323846 * i / (kProtoLen - 1));
            osTaps_[p][k] = static_cast<float>(sincPi(centred / kLimiterOversample) * hann);
        }
    }

    delayL_.assign(static_cast<std::size_t>(latency_), 0.f);
    delayR_.assign(static_cast<std::size_t>(latency_), 0.f);
    wedgeVal_.assign(static_cast<std::size_t>(window_), 1.f);
    wedgeIdx_.assign(static_cast<std::size_t>(window_), 0);
    avgRing_.assign(static_cast<std::size_t>(window_), 1.f);
    reset();
}

void Limiter::clearState() {
    std::fill(delayL_.begin(), delayL_.end(), 0.f);
    std::fill(delayR_.begin(), delayR_.end(), 0.f);
    delayPos_ = 0;
    for (int i = 0; i < kLimiterOsTaps; ++i) { histL_[i] = 0.f; histR_[i] = 0.f; }
    histPos_ = 0;
    wedgeHead_ = 0;
    wedgeCount_ = 0;
    pushIndex_ = 0;
    std::fill(avgRing_.begin(), avgRing_.end(), 1.f);
    avgSum_ = static_cast<double>(window_);
    avgPos_ = 0;
    releaseGain_ = 1.f;
}

void Limiter::reset() {
    clearState();
    active_ = target_.load(std::memory_order_relaxed);
    mix_ = active_ ? 1.f : 0.f;
    primeRemaining_ = 0;
}

void Limiter::setEnabled(bool on) { target_.store(on, std::memory_order_relaxed); }

float Limiter::detectTruePeak(float xl, float xr, bool stereo) {
    histPos_ = (histPos_ + 1) & (kLimiterOsTaps - 1);
    histL_[histPos_] = xl;
    if (stereo) histR_[histPos_] = xr;

    // Phase 0 is the exact delay, so the sample under test is read, not convolved.
    const int base = (histPos_ - kLimiterOsDelay + kLimiterOsTaps) & (kLimiterOsTaps - 1);
    float peak = std::fabs(histL_[base]);
    if (stereo) {
        const float r0 = std::fabs(histR_[base]);
        if (r0 > peak) peak = r0;
    }
    for (int p = 1; p < kLimiterOversample; ++p) {
        float accL = 0.f, accR = 0.f;
        for (int k = 0; k < kLimiterOsTaps; ++k) {
            const int idx = (histPos_ - k + kLimiterOsTaps) & (kLimiterOsTaps - 1);
            accL += osTaps_[p][k] * histL_[idx];
            if (stereo) accR += osTaps_[p][k] * histR_[idx];
        }
        const float al = std::fabs(accL);
        if (al > peak) peak = al;
        if (stereo) {
            const float ar = std::fabs(accR);
            if (ar > peak) peak = ar;
        }
    }
    return peak;
}

float Limiter::smoothGain(float target) {
    // Sliding minimum over `window_` via a monotonic wedge. Expiring the front BEFORE the push
    // is what bounds the wedge to `window_` entries — pushing first can lap the ring. Wraps by
    // compare-and-subtract, matching delayPos_/avgPos_: window_ is not a power of two, so `%`
    // would not strength-reduce on this per-sample path.
    while (wedgeCount_ > 0 &&
           wedgeIdx_[static_cast<std::size_t>(wedgeHead_)] <= pushIndex_ - window_) {
        wedgeHead_ = (wedgeHead_ + 1 == window_) ? 0 : wedgeHead_ + 1;
        --wedgeCount_;
    }
    while (wedgeCount_ > 0) {
        const int backSum = wedgeHead_ + wedgeCount_ - 1;
        const int back = (backSum >= window_) ? backSum - window_ : backSum;
        if (wedgeVal_[static_cast<std::size_t>(back)] < target) break;
        --wedgeCount_;
    }
    const int slotSum = wedgeHead_ + wedgeCount_;
    const int slot = (slotSum >= window_) ? slotSum - window_ : slotSum;
    wedgeVal_[static_cast<std::size_t>(slot)] = target;
    wedgeIdx_[static_cast<std::size_t>(slot)] = pushIndex_;
    ++wedgeCount_;
    ++pushIndex_;
    const float windowMin = wedgeVal_[static_cast<std::size_t>(wedgeHead_)];

    // Moving average of the same width over those minima.
    avgSum_ += static_cast<double>(windowMin) - static_cast<double>(avgRing_[static_cast<std::size_t>(avgPos_)]);
    avgRing_[static_cast<std::size_t>(avgPos_)] = windowMin;
    avgPos_ = (avgPos_ + 1 == window_) ? 0 : avgPos_ + 1;
    float smoothed = static_cast<float>(avgSum_ / window_);
    // Never above unity — the structural form of "no makeup gain, ever", and what makes the
    // at-rest gain land on EXACTLY 1.0f after the running sum has been added to and subtracted
    // from for hours.
    if (!(smoothed < 1.f)) smoothed = 1.f;

    // Release: falls with the smoother, rises no faster than the one-pole. Staying at or below
    // `smoothed` is what preserves the no-overshoot bound.
    if (smoothed < releaseGain_) releaseGain_ = smoothed;
    else releaseGain_ += (smoothed - releaseGain_) * releaseCoeff_;
    return releaseGain_;
}

float Limiter::process(float* left, float* right, int frames) {
    if (!left || frames <= 0 || latency_ <= 0) return 1.f;
    const bool want = target_.load(std::memory_order_relaxed);
    if (!want && !active_) return 1.f;  // settled bypass: not one sample read or written
    if (want && !active_) {
        // A live engage. Start dry, fill the delay line, then crossfade — so the wet path is
        // never silence weighted above zero.
        clearState();
        active_ = true;
        mix_ = 0.f;
        primeRemaining_ = latency_;
    }

    const bool stereo = (right != nullptr);
    float blockMin = 1.f;
    for (int i = 0; i < frames; ++i) {
        const float dryL = left[i];
        const float dryR = stereo ? right[i] : 0.f;

        const float peak = detectTruePeak(dryL, dryR, stereo);
        const float targetGain = peak > ceiling_ ? ceiling_ / peak : 1.f;
        const float gain = smoothGain(targetGain);

        const std::size_t slot = static_cast<std::size_t>(delayPos_);
        const float wetL = delayL_[slot] * gain;
        const float wetR = stereo ? delayR_[slot] * gain : 0.f;
        delayL_[slot] = dryL;
        if (stereo) delayR_[slot] = dryR;
        delayPos_ = (delayPos_ + 1 == latency_) ? 0 : delayPos_ + 1;

        // The endpoints are branches rather than blend arithmetic so a settled state is exact:
        // dry + (wet - dry) * 1.0f is not wet in floating point. At m <= 0 the buffer is left
        // untouched, which is the dry sample already in it.
        const float m = mix_;
        // The reported minimum is the gain actually reaching the output, not the limiter's raw
        // target — mid-crossfade only a fraction `m` of the reduction is audible, so the meter
        // (whose contract is "smallest gain APPLIED") must blend the same way the signal does:
        // unity at m=0, `gain` at m=1, linear between.
        const float effectiveGain = 1.f - m + m * gain;
        if (effectiveGain < blockMin) blockMin = effectiveGain;
        if (m >= 1.f) {
            left[i] = wetL;
            if (stereo) right[i] = wetR;
        } else if (m > 0.f) {
            left[i] = dryL + (wetL - dryL) * m;
            if (stereo) right[i] = dryR + (wetR - dryR) * m;
        }

        if (primeRemaining_ > 0) {
            --primeRemaining_;
        } else if (want) {
            mix_ = (mix_ + mixStep_ >= 1.f) ? 1.f : mix_ + mixStep_;
        } else {
            mix_ = (mix_ - mixStep_ <= 0.f) ? 0.f : mix_ - mixStep_;
        }
    }
    if (!want && mix_ <= 0.f && primeRemaining_ == 0) active_ = false;
    return blockMin;
}

}  // namespace reasampler::instrument::engine