#!/usr/bin/env python3
"""soothe2 поведенческий симулятор.
Модель: out = x − amount(t)·bp(x; fc, Q)
  bp — 2nd-order bandpass, пик gain 1.0 на fc; Q из selectivity.
  amount(t) = amt_total·env(t), env — one-pole со скоростью attack (в рост) / release (в спад).
  env управляется энергией полосы (детектор): bp-loudness > порога → целимся в amt, иначе в 0.
Кривые из fit_curves.py / измерения.
"""
import numpy as np
import wave

def _lut_sat(xs, ys, x):
    return np.interp(x, xs, ys, left=ys[0], right=ys[-1])

# --- измеренные кривые как LUT (максимальная точность) ---
_DEPTH_LUT = (np.array([0.5, 0.864, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0]),
              np.array([0.5605, 0.5879, 0.5977, 0.6673, 0.7274, 0.8199, 0.9387, 0.9893]))
_SENS_LUT = (np.array([0.0, 6.0, 12.0, 24.0, 48.0]),
             np.array([0.4814, 0.7597, 1.0, 1.0, 1.0]))
_SHARP_LUT = (np.array([1.0, 3.0, 5.0, 10.0, 20.0]),
              np.array([0.2041, 0.5913, 0.8386, 1.0, 1.0]))

def amount_total(depth=0.864, sens=12.0, sharp=10.0, mode=1):
    A = _lut_sat(*_DEPTH_LUT, depth)
    S = _lut_sat(*_SENS_LUT, sens)
    H = _lut_sat(*_SHARP_LUT, sharp)
    M = 0.745 if mode == 0 else 1.0
    return A * S * H * M

def read_wav(path):
    w = wave.open(path, 'rb')
    sw, nc, n = w.getsampwidth(), w.getnchannels(), w.getnframes()
    d = np.frombuffer(w.readframes(n), dtype=np.uint8).reshape(n, nc, sw)
    w.close()
    ch = d[:, 0, :]
    if sw == 3:
        v = (ch[:, 0].astype(np.int64) | (ch[:, 1].astype(np.int64) << 8) |
             (ch[:, 2].astype(np.int64) << 16))
        v = (v ^ (1 << 23)) - (1 << 23)
        return v.astype(np.float64) / (1 << 23)
    elif sw == 2:
        v = (ch[:, 0].astype(np.int64) | (ch[:, 1].astype(np.int64) << 8))
        v = (v ^ (1 << 15)) - (1 << 15)
        return v.astype(np.float64) / (1 << 15)
    else:
        v = ch[:, 0].astype(np.float64)
        return (v - 128) / 128.0

def write_wav(path, x, sr=44100):
    x = np.clip(x, -1, 1)
    frames = (x.astype(np.float64) * (1 << 23)).astype(np.int64)
    out = np.empty((frames.size, 1, 3), dtype=np.uint8)
    out[:, 0, 0] = (frames & 0xFF).astype(np.uint8)
    out[:, 0, 1] = ((frames >> 8) & 0xFF).astype(np.uint8)
    out[:, 0, 2] = ((frames >> 16) & 0xFF).astype(np.uint8)
    w = wave.open(path, 'wb')
    w.setnchannels(1); w.setsampwidth(3); w.setframerate(sr)
    w.writeframes(np.ascontiguousarray(out).tobytes()); w.close()

def tau_attack(attack):
    return 0.020 * np.exp(attack / 1.955)     # a=0→20ms, a=5→258ms, a=10→3.3s

def tau_release(release):
    # измеренные tau (экспоненциальный фит спада к floor 0.2dB)
    xs = np.array([0.0, 1.0, 2.0, 5.0, 10.0, 100.0])
    ys = np.array([0.027, 0.08, 0.11, 0.14, 15.0, 15.0])
    return float(np.interp(release, xs, ys))

def q_from_sel(sel):
    # лучший фит ширины (probes-зонд): sel1→3.5, sel8→6.0 (2-порядка полоса)
    return 3.0 + 0.36 * min(sel, 10.0)

# --- bandpass biquad: peak gain = 1.0 на fc (RBJ normalized) ---
def bp_coeffs(fc, q, sr):
    w0 = 2 * np.pi * fc / sr
    alpha = np.sin(w0) / (2 * q)
    b = np.array([alpha, 0.0, -alpha]) / alpha   # [1, 0, -1]
    a = np.array([1 + alpha, -2 * np.cos(w0), 1 - alpha])
    an = a / a[0]   # normalize denominator: a[0]=1
    # численное сканирование → норм на пик=1 (fn scale)
    w = np.linspace(w0 * 0.4, w0 * 1.6, 8000)
    H = np.abs(b[0] / a[0] * (1 - np.exp(-2j * w)) /
               (1 + an[1] * np.exp(-1j * w) + an[2] * np.exp(-2j * w)))
    G = H.max()
    bn = b / a[0] / G
    return bn, an

def apply_bp(x, fc, q, sr):
    b, a = bp_coeffs(fc, q, sr)
    b0, b1, b2 = b
    a1, a2 = a[1], a[2]
    y = np.zeros_like(x)
    z1 = z2 = 0.0
    for i in range(x.size):
        v = x[i]
        out = b0 * v + z1
        z1 = b1 * v - a1 * out + z2
        z2 = b2 * v - a2 * out
        y[i] = out
    return y

def simulate(x, sr=44100, fc=500.0, depth=0.864, sens=12.0, sharp=10.0,
             sel=10.0, mode=1, attack=0.0, release=0.0, mix=100.0, thr=0.010):
    amt = amount_total(depth, sens, sharp, mode)
    q = q_from_sel(sel)
    ta, tr = tau_attack(attack), tau_release(release)
    a_alpha = 1 - np.exp(-1.0 / (ta * sr))
    r_alpha = 1 - np.exp(-1.0 / (tr * sr))
    bp = apply_bp(x, fc, q, sr)
    # детектор: мгновенная амплитуда полосы → smoothed loudness
    if len(x) > 0:
        loud = np.abs(bp)
        N = int(0.005 * sr)
        win = np.ones(N) / N
        loud = np.convolve(loud, win, mode='same')
    env = 0.0
    out = np.zeros_like(x)
    lp = 0.0
    for i in range(x.size):
        tar = amt if loud[i] > thr else 0.0
        env += (a_alpha if tar > env else r_alpha) * (tar - env)
        out[i] = x[i] - env * bp[i]
    if mix < 100:
        out = (100 - mix) / 100.0 * x + mix / 100.0 * out
    return out