93 lines
3.2 KiB
Python
93 lines
3.2 KiB
Python
"""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 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)
|