"""EBU R128 loudness + true peak (D1).

- integrated_lufs: K-weighting (BS.1770) + 400 ms windows, 100 ms hop,
  absolute (-70 LUFS) + relative (-10 LU) gating -> LUFS.
- true_peak_db: 4x oversample (Kaiser-16) -> max |x| in dBTP.

Stereo-only (channel weight 1.0 L/R); mono broadcast to 2 ch = same result.
"""
import numpy as np
from scipy.signal import lfilter, resample_poly


def _k_weighting_coeffs(sr: float):
    """Two biquad stages per ITU-R BS.1770-4 (RBJ cookbook form)."""
    # Stage 1: high shelf +4 dB, f0=1681.97 Hz, Q=0.7072
    G = 3.999843853973347
    f0 = 1681.974450955533
    Q = 0.7071752369554196
    A = 10 ** (G / 40.0)
    w0 = 2 * np.pi * f0 / sr
    alpha = np.sin(w0) / (2 * Q)
    cosw = np.cos(w0)
    b = np.array([
        A * ((A + 1) + (A - 1) * cosw + 2 * np.sqrt(A) * alpha),
        -2 * A * ((A - 1) + (A + 1) * cosw),
        A * ((A + 1) + (A - 1) * cosw - 2 * np.sqrt(A) * alpha),
    ])
    a = np.array([
        (A + 1) - (A - 1) * cosw + 2 * np.sqrt(A) * alpha,
        2 * ((A - 1) - (A + 1) * cosw),
        (A + 1) - (A - 1) * cosw - 2 * np.sqrt(A) * alpha,
    ])
    # Stage 2: high-pass f0=38.14 Hz, Q=0.5003
    f1 = 38.13547087602444
    Q1 = 0.5003270373238773
    w1 = 2 * np.pi * f1 / sr
    alpha1 = np.sin(w1) / (2 * Q1)
    cosw1 = np.cos(w1)
    hp_b = np.array([(1 + cosw1) / 2, -(1 + cosw1), (1 + cosw1) / 2])
    hp_a = np.array([1 + alpha1, -2 * cosw1, 1 - alpha1])
    return b / a[0], a / a[0], hp_b, hp_a


def _k_weighted(audio: np.ndarray, sr: float) -> np.ndarray:
    b1, a1, b2, a2 = _k_weighting_coeffs(sr)
    out = np.empty_like(audio)
    for c in range(audio.shape[0]):
        out[c] = lfilter(b2, a2, lfilter(b1, a1, audio[c]))
    return out


def integrated_lufs(audio: np.ndarray, sr: float) -> float:
    """Integrated loudness in LUFS (EBU R128 / BS.1770-4). -inf for silence."""
    audio = np.asarray(audio, dtype=np.float64)
    if audio.ndim == 1:
        audio = audio[np.newaxis, :]
    if audio.shape[1] == 0:
        return -np.inf
    k = _k_weighted(audio, sr)
    block = int(round(0.4 * sr))
    hop = int(round(0.1 * sr))
    n = max(1, (k.shape[1] - block) // hop + 1)
    powers = np.empty(n)
    for i in range(n):
        seg = k[:, i * hop:i * hop + block]
        powers[i] = np.mean(seg ** 2)
    # Absolute gate: -70 LUFS -> power = 10^((-70 + 0.691)/10)
    abs_gate = 10 ** ((-70.0 + 0.691) / 10.0)
    above_abs = powers[powers > abs_gate]
    if above_abs.size == 0:
        return -np.inf
    abs_loud = -0.691 + 10.0 * np.log10(np.mean(above_abs))
    # Relative gate: abs_loud - 10 LU
    rel_gate = 10 ** ((abs_loud - 10.0 + 0.691) / 10.0)
    gated = above_abs[above_abs > rel_gate]
    if gated.size == 0:
        gated = above_abs
    return float(-0.691 + 10.0 * np.log10(np.mean(gated)))


def mono_side_ratio(audio: np.ndarray) -> float:
    """Side/Mid RMS ratio (D3 mono-compat). 0 = mono-safe (L==R); >1 = side
    louder than mid — mono sum mất nhiều âm thanh."""
    audio = np.asarray(audio, dtype=np.float64)
    if audio.ndim == 1:
        audio = audio[np.newaxis, :]
    if audio.shape[0] < 2 or audio.shape[1] == 0:
        return 0.0
    mid = (audio[0] + audio[1]) * 0.5
    side = (audio[0] - audio[1]) * 0.5
    m = float(np.sqrt(np.mean(mid ** 2))) + 1e-12
    return float(np.sqrt(np.mean(side ** 2))) / m


def true_peak_db(audio: np.ndarray, sr: float, oversample: int = 4) -> float:
    """True peak in dBTP (4x oversample). -inf for digital silence."""
    audio = np.asarray(audio, dtype=np.float64)
    if audio.ndim == 1:
        audio = audio[np.newaxis, :]
    if audio.shape[1] == 0:
        return -np.inf
    up = resample_poly(audio, oversample, 1, axis=1, window=("kaiser", 16))
    peak = float(np.max(np.abs(up)))
    if peak <= 1e-12:
        return -np.inf
    return 20.0 * np.log10(peak)