#include "detect.hpp"
#include <cmath>
#include <cstring>
#include <algorithm>

Detector::Detector(size_t nfft, float sample_rate)
    : nfft_(nfft), sample_rate_(sample_rate),
      sharpness_(1.0f), selectivity_(0.5f), depth_(0.0f) {
    envelope_.resize(nfft, 0.0f);
    prev_mask_.resize(nfft, 1.0f);
    smooth_buf_.resize(nfft, 0.0f);
}

void Detector::setParams(float sharpness, float selectivity, float depth) {
    sharpness_ = sharpness;
    selectivity_ = selectivity;
    depth_ = depth;
}

// Bell curve weight: H_q(f; fc, Qeff)
// H(f) = 1/sqrt(1 + (Qeff * A)^2) where A = f/fc - fc/f
static float bell_curve(float freq, float fc, float qeff) {
    if (freq <= 0.0f || fc <= 0.0f) return 0.0f;
    float a = freq / fc - fc / freq;
    float qa = qeff * a;
    return 1.0f / std::sqrt(1.0f + qa * qa);
}

// Floor function: base reduction depending on input level (dBFS)
// From measured data: floor(L) ≈ 1.972 + 0.2584*(L+24)
static float floor_func(float level_db) {
    return 1.972f + 0.2584f * (level_db + 24.0f);
}

void Detector::processFrame(const std::complex<double>* spectrum, float* mask) {
    size_t half = nfft_ / 2;
    float bin_hz = sample_rate_ / static_cast<float>(nfft_);

    // Compute magnitude per bin
    std::vector<float> mag(half + 1);
    for (size_t i = 0; i <= half; i++) {
        mag[i] = static_cast<float>(std::sqrt(
            spectrum[i].real() * spectrum[i].real() +
            spectrum[i].imag() * spectrum[i].imag()));
    }

    // Regional level: RMS over a moving window of bins, calibrated to dBFS
    // Empirically: peek level -8 dBFS → floor 6.1; below ~-32 dBFS → no cut
    const size_t reg = 16;
    std::vector<float> reg_db(half + 1);
    for (size_t i = 0; i <= half; i++) {
        size_t lo = (i >= reg) ? i - reg : 0;
        size_t hi = std::min(half, i + reg);
        float e = 0.0f;
        for (size_t k = lo; k <= hi; k++) e += mag[k] * mag[k];
        float rms = std::sqrt(e / (hi - lo + 1u));
        reg_db[i] = 20.0f * std::log10(std::max(rms, 1e-10f)) - 26.1f;
    }

    float qeff = 1.54f * std::pow(std::max(sharpness_, 0.5f), 1.33f);
    float sens_weight = 6.02f * std::min(1.0f, 12.0f / 12.0f);  // sens=12 → full

    for (size_t i = 0; i <= half; i++) {
        // Level-dependent floor: louder bins cut deeper (matches soothe mask shape)
        float floor_red = std::max(0.0f, floor_func(reg_db[i]));

        // Total reduction in dB (no additive band bell — tone boost is implicit
        // via higher regional level at the tone frequency)
        float total_red = 1.041f * depth_ * floor_red;

        // Convert to linear mask: mask = 10^(-total_red/20)
        mask[i] = std::pow(10.0f, -total_red / 20.0f);
    }

    // Mirror for negative frequencies
    for (size_t i = half + 1; i < nfft_; i++) {
        mask[i] = mask[nfft_ - i];
    }

    // Temporal smoothing
    const float smooth_alpha = 0.3f;
    for (size_t i = 0; i < nfft_; i++) {
        mask[i] = prev_mask_[i] + smooth_alpha * (mask[i] - prev_mask_[i]);
        prev_mask_[i] = mask[i];
    }
}