898 lines
35 KiB
Python
898 lines
35 KiB
Python
"""DSP 核心：线性相位滤波、电平计量、真峰值、ITU-R BS.1770 响度、离线处理。
|
||
|
||
设计原则是"录音链路尽量透明"：
|
||
|
||
* 除用户显式开启的功能（低切、门限、归一化）外，不做任何隐式处理；
|
||
* 实时链路上的滤波器一律使用**线性相位 FIR**，不引入相位失真；
|
||
* 所有电平/响度计量按国际标准实现（真峰值 4 倍过采样、BS.1770 K 加权）；
|
||
* 不依赖 scipy，只用 numpy，便于打包分发。
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import math
|
||
from collections import deque
|
||
from dataclasses import dataclass, field
|
||
|
||
import numpy as np
|
||
|
||
__all__ = [
|
||
"db_to_lin", "lin_to_db", "rms", "peak",
|
||
"design_highpass_fir", "design_lowpass_fir", "FIRFilter",
|
||
"true_peak", "true_peak_dbfs", "dbfs",
|
||
"kweighting_coeffs", "kweighting_response_sq",
|
||
"LoudnessMeter", "LevelMeter", "analyze_data", "analyze_file",
|
||
"remove_dc", "trim_silence", "normalize_peak", "normalize_loudness",
|
||
"fade_edges", "noise_gate", "mixdown_mono", "highpass_offline",
|
||
"SILENCE_DBFS_FLOOR",
|
||
]
|
||
|
||
# 低于该电平即可视为数字静音（16 位理论本底约 -96 dBFS）
|
||
SILENCE_DBFS_FLOOR = -120.0
|
||
|
||
|
||
# --------------------------------------------------------------------- 基础
|
||
def db_to_lin(db: float) -> float:
|
||
return float(10.0 ** (float(db) / 20.0))
|
||
|
||
|
||
def lin_to_db(x: float, floor: float = SILENCE_DBFS_FLOOR) -> float:
|
||
x = float(x)
|
||
if x <= 0 or not math.isfinite(x):
|
||
return float(floor)
|
||
v = 20.0 * math.log10(x)
|
||
return v if v > floor else float(floor)
|
||
|
||
|
||
def dbfs(x: np.ndarray | float, floor: float = SILENCE_DBFS_FLOOR) -> np.ndarray | float:
|
||
"""线性幅度 -> dBFS（数组安全）。"""
|
||
arr = np.asarray(x, dtype=np.float64)
|
||
with np.errstate(divide="ignore", invalid="ignore"):
|
||
out = 20.0 * np.log10(np.abs(arr))
|
||
out = np.where(np.isfinite(out) & (out > floor), out, floor)
|
||
return out if isinstance(x, np.ndarray) else float(out)
|
||
|
||
|
||
def rms(x: np.ndarray, axis=None) -> np.ndarray | float:
|
||
arr = np.asarray(x, dtype=np.float64)
|
||
if arr.size == 0:
|
||
return 0.0
|
||
return np.sqrt(np.mean(np.square(arr), axis=axis))
|
||
|
||
|
||
def peak(x: np.ndarray, axis=None) -> np.ndarray | float:
|
||
arr = np.asarray(x, dtype=np.float64)
|
||
if arr.size == 0:
|
||
return 0.0
|
||
return np.max(np.abs(arr), axis=axis)
|
||
|
||
|
||
# ------------------------------------------------------------------- 滤波器
|
||
def _windowed_sinc(numtaps: int, cutoff: float, *, highpass: bool,
|
||
beta: float = 8.6) -> np.ndarray:
|
||
"""窗函数法设计线性相位 FIR。``cutoff`` 为归一化频率（cycles/sample）。"""
|
||
if numtaps % 2 == 0:
|
||
numtaps += 1 # 奇数长度 -> 整数群延迟，且可表示真正的直流零点
|
||
n = np.arange(numtaps) - (numtaps - 1) / 2.0
|
||
h = 2.0 * cutoff * np.sinc(2.0 * cutoff * n)
|
||
if highpass:
|
||
h = -h
|
||
h[(numtaps - 1) // 2] += 1.0
|
||
h *= np.kaiser(numtaps, beta)
|
||
if highpass:
|
||
# 强制直流增益为严格 0：彻底消除直流漂移，同时保持线性相位（对称加权）
|
||
h -= h.mean()
|
||
return h
|
||
|
||
|
||
def design_highpass_fir(cutoff_hz: float, samplerate: int, *,
|
||
numtaps: int | None = None,
|
||
beta: float = 8.0,
|
||
max_taps: int = 8193) -> np.ndarray:
|
||
"""设计低切（高通）线性相位 FIR，通带增益归一到 0 dB。
|
||
|
||
抽头数按 Kaiser 公式反推，保证阻带衰减约 80 dB 且过渡带足够窄
|
||
（过渡带宽度取截止频率的一半），这样 100 Hz 以上的语音基频几乎不受影响。
|
||
"""
|
||
if cutoff_hz <= 0:
|
||
raise ValueError("截至频率必须为正")
|
||
nyq = samplerate / 2.0
|
||
if cutoff_hz >= nyq * 0.9:
|
||
raise ValueError("截至频率过高")
|
||
if numtaps is None:
|
||
atten = max(40.0, 9.0 * beta) # beta 与阻带衰减的经验对应
|
||
trans = max(cutoff_hz * 0.5, 4.0) # 过渡带宽度 (Hz)
|
||
d_omega = 2.0 * math.pi * trans / samplerate
|
||
numtaps = int(math.ceil((atten - 7.95) / (2.285 * d_omega)))
|
||
numtaps = int(np.clip(numtaps, 129, max_taps))
|
||
h = _windowed_sinc(int(numtaps), cutoff_hz / samplerate,
|
||
highpass=True, beta=beta)
|
||
# 用 1 kHz（或奈奎斯特的 10%）处的响应把通带增益归一到 1.0
|
||
ref = min(1000.0, nyq * 0.1)
|
||
w = 2.0 * math.pi * ref / samplerate
|
||
k = np.arange(h.size)
|
||
gain = np.abs(np.sum(h * np.exp(-1j * w * k)))
|
||
if gain > 1e-12:
|
||
h = h / gain
|
||
return h
|
||
|
||
|
||
def design_lowpass_fir(cutoff_hz: float, samplerate: int,
|
||
numtaps: int = 129, beta: float = 8.6) -> np.ndarray:
|
||
h = _windowed_sinc(int(numtaps), cutoff_hz / samplerate,
|
||
highpass=False, beta=beta)
|
||
s = h.sum()
|
||
return h / s if abs(s) > 1e-12 else h
|
||
|
||
|
||
class FIRFilter:
|
||
"""多通道流式 FIR（重叠相加），保持跨块状态。"""
|
||
|
||
def __init__(self, coeffs: np.ndarray, channels: int):
|
||
self.h = np.asarray(coeffs, dtype=np.float64)
|
||
self.channels = int(channels)
|
||
self._tail = np.zeros((self.h.size - 1, self.channels), dtype=np.float64)
|
||
|
||
@property
|
||
def latency_samples(self) -> int:
|
||
return (self.h.size - 1) // 2
|
||
|
||
def reset(self) -> None:
|
||
self._tail[:] = 0.0
|
||
|
||
def process(self, block: np.ndarray) -> np.ndarray:
|
||
x = np.asarray(block, dtype=np.float64)
|
||
if x.ndim == 1:
|
||
x = x[:, None]
|
||
n = x.shape[0]
|
||
out = np.empty((n, self.channels), dtype=np.float64)
|
||
for c in range(self.channels):
|
||
y = np.convolve(x[:, c], self.h)
|
||
y[:self.h.size - 1] += self._tail[:, c]
|
||
out[:, c] = y[:n]
|
||
self._tail[:, c] = y[n:n + self.h.size - 1]
|
||
return out
|
||
|
||
|
||
# ------------------------------------------------------------------ 真峰值
|
||
_OVERSAMPLE_CACHE: dict[tuple[int, int], list[np.ndarray]] = {}
|
||
|
||
|
||
def _polyphase_bank(factor: int = 4, taps_per_phase: int = 12) -> list[np.ndarray]:
|
||
"""L 倍过采样的多相插值滤波器组。
|
||
|
||
``h_p[j] = g[p + jL]``，其中 ``g`` 是以 1/(2L) 为截止频率的窗函数 sinc
|
||
（在过采样后的采样率下）。零相位支路 ``h_0`` 即近似恒等，
|
||
因此插值结果绝不会低于原始采样峰值。
|
||
"""
|
||
key = (factor, taps_per_phase)
|
||
if key in _OVERSAMPLE_CACHE:
|
||
return _OVERSAMPLE_CACHE[key]
|
||
numtaps = factor * taps_per_phase + 1
|
||
n = np.arange(numtaps) - (numtaps - 1) / 2.0
|
||
# 关键：插值核峰值必须为 1（不能用 sum 归一化），这样零相位支路恰好是
|
||
# 恒等滤波，插值结果的峰值才不会低于原始采样峰值。
|
||
g = np.sinc(n / factor) * np.kaiser(numtaps, 9.0)
|
||
bank = [g[p::factor].copy() for p in range(factor)]
|
||
_OVERSAMPLE_CACHE[key] = bank
|
||
return bank
|
||
|
||
|
||
def true_peak(x: np.ndarray, *, factor: int = 4) -> np.ndarray:
|
||
"""按 ITU-R BS.1770 附录 2 计算真峰值（返回每声道线性幅度）。
|
||
|
||
做法是多相插值到 ``factor`` 倍采样率后取最大绝对值，
|
||
可以捕捉到采样点之间的过冲（inter-sample peak）。
|
||
"""
|
||
arr = np.asarray(x, dtype=np.float64)
|
||
if arr.ndim == 1:
|
||
arr = arr[:, None]
|
||
if arr.size == 0:
|
||
return np.zeros(arr.shape[1])
|
||
bank = _polyphase_bank(factor)
|
||
out = np.zeros(arr.shape[1])
|
||
for c in range(arr.shape[1]):
|
||
col = arr[:, c]
|
||
best = 0.0
|
||
for h in bank:
|
||
y = np.convolve(col, h)
|
||
m = float(np.max(np.abs(y))) if y.size else 0.0
|
||
if m > best:
|
||
best = m
|
||
out[c] = best
|
||
return out
|
||
|
||
|
||
def true_peak_dbfs(x: np.ndarray, *, factor: int = 4) -> np.ndarray:
|
||
return dbfs(true_peak(x, factor=factor))
|
||
|
||
|
||
# ------------------------------------------------------- BS.1770 K 加权响度
|
||
# ITU-R BS.1770-4 在 48 kHz 下给出的参考系数（用于自检比对）
|
||
_KW_48K_SHELF = ([1.53512485958697, -2.69169618940638, 1.19839281085285],
|
||
[1.0, -1.69065929318241, 0.73248077421585])
|
||
_KW_48K_HPF = ([1.0, -2.0, 1.0],
|
||
[1.0, -1.99004745483398, 0.99007225036621])
|
||
|
||
# 等效设计参数（与 ITU 系数在 48 kHz 下逐位吻合，见自检用例）
|
||
_SHELF_GAIN_DB = 3.999843853973347
|
||
_SHELF_FC = 1681.974450955533
|
||
_SHELF_Q = 0.7071752369554196
|
||
_SHELF_VB_EXP = 0.4996667741545416
|
||
_HPF_FC = 38.13547087602444
|
||
_HPF_Q = 0.5003270373238773
|
||
|
||
|
||
def _itu_highshelf(samplerate: int):
|
||
"""BS.1770 第一级：高频搁架（+4 dB @ 高频段）。
|
||
|
||
采用 ITU 参考实现所用的参数化公式，在 48 kHz 下与标准给定系数
|
||
完全一致（误差 < 1e-15），其它采样率下即为标准的双线性变换推广。
|
||
"""
|
||
k = math.tan(math.pi * _SHELF_FC / samplerate)
|
||
vh = 10.0 ** (_SHELF_GAIN_DB / 20.0)
|
||
vb = vh ** _SHELF_VB_EXP
|
||
a0 = 1.0 + k / _SHELF_Q + k * k
|
||
b = np.array([(vh + vb * k / _SHELF_Q + k * k) / a0,
|
||
2.0 * (k * k - vh) / a0,
|
||
(vh - vb * k / _SHELF_Q + k * k) / a0])
|
||
a = np.array([1.0,
|
||
2.0 * (k * k - 1.0) / a0,
|
||
(1.0 - k / _SHELF_Q + k * k) / a0])
|
||
return b, a
|
||
|
||
|
||
def _itu_highpass(samplerate: int):
|
||
"""BS.1770 第二级：约 38 Hz 的二阶高通。"""
|
||
k = math.tan(math.pi * _HPF_FC / samplerate)
|
||
a0 = 1.0 + k / _HPF_Q + k * k
|
||
b = np.array([1.0, -2.0, 1.0])
|
||
a = np.array([1.0,
|
||
2.0 * (k * k - 1.0) / a0,
|
||
(1.0 - k / _HPF_Q + k * k) / a0])
|
||
return b, a
|
||
|
||
|
||
def kweighting_coeffs(samplerate: int):
|
||
"""返回 K 加权的两级 biquad 系数 ``((b1,a1),(b2,a2))``。"""
|
||
if int(samplerate) == 48000:
|
||
# 48 kHz 直接使用标准原文系数，保证逐位一致
|
||
b1, a1 = _KW_48K_SHELF
|
||
b2, a2 = _KW_48K_HPF
|
||
return (np.array(b1, dtype=np.float64), np.array(a1, dtype=np.float64)), \
|
||
(np.array(b2, dtype=np.float64), np.array(a2, dtype=np.float64))
|
||
return _itu_highshelf(int(samplerate)), _itu_highpass(int(samplerate))
|
||
|
||
|
||
def biquad_response_sq(b: np.ndarray, a: np.ndarray,
|
||
w: np.ndarray) -> np.ndarray:
|
||
"""|H(e^{jw})|^2，``w`` 为归一化角频率数组。"""
|
||
z1 = np.exp(-1j * w)
|
||
z2 = z1 * z1
|
||
num = b[0] + b[1] * z1 + b[2] * z2
|
||
den = a[0] + a[1] * z1 + a[2] * z2
|
||
return np.abs(num / den) ** 2
|
||
|
||
|
||
def kweighting_response_sq(samplerate: int, fft_size: int) -> np.ndarray:
|
||
"""K 加权滤波器在 rfft 各 bin 上的功率响应（长度 ``fft_size//2+1``）。"""
|
||
(b1, a1), (b2, a2) = kweighting_coeffs(samplerate)
|
||
w = 2.0 * math.pi * np.arange(fft_size // 2 + 1) / fft_size
|
||
return biquad_response_sq(b1, a1, w) * biquad_response_sq(b2, a2, w)
|
||
|
||
|
||
# BS.1770 的绝对门限与相对门限
|
||
_LUFS_ABS_GATE = -70.0
|
||
_LUFS_REL_GATE = -10.0
|
||
_LUFS_OFFSET = -0.691
|
||
|
||
|
||
def _lufs_from_ms(mean_square: float) -> float:
|
||
if mean_square <= 0:
|
||
return float("-inf")
|
||
return _LUFS_OFFSET + 10.0 * math.log10(mean_square)
|
||
|
||
|
||
class LoudnessMeter:
|
||
"""流式响度计（ITU-R BS.1770-4 门限算法）。
|
||
|
||
以 100 ms 为 hop、400 ms 为瞬时块（即 75% 重叠）统计能量，
|
||
再按绝对门限 -70 LUFS、相对门限 -10 LU 计算整体响度。
|
||
为实现高效，K 加权在频域按 |H(f)|² 施加——对"能量"计量而言与
|
||
时域滤波完全等价（K 加权滤波器冲激响应仅数毫秒，块边界残留可忽略）。
|
||
"""
|
||
|
||
HOP_SECONDS = 0.1
|
||
MOMENTARY_HOPS = 4 # 400 ms
|
||
SHORT_TERM_HOPS = 30 # 3 s
|
||
|
||
def __init__(self, samplerate: int, channels: int, *,
|
||
keep_hops: bool = True):
|
||
self.samplerate = int(samplerate)
|
||
self.channels = int(channels)
|
||
self.hop = max(1, int(round(self.HOP_SECONDS * self.samplerate)))
|
||
self._resp = kweighting_response_sq(self.samplerate, self.hop)
|
||
# rfft 的单边能量补正（直流与奈奎斯特 bin 不翻倍）
|
||
self._fold = np.full(self._resp.size, 2.0)
|
||
self._fold[0] = 1.0
|
||
if self.hop % 2 == 0:
|
||
self._fold[-1] = 1.0
|
||
self._weight = self._resp * self._fold / float(self.hop * self.hop)
|
||
|
||
self._acc = np.zeros((0, self.channels), dtype=np.float64)
|
||
self._hops: deque[float] = deque(maxlen=self.SHORT_TERM_HOPS)
|
||
self._all_hops: list[float] = [] if keep_hops else []
|
||
self._keep_hops = keep_hops
|
||
self.processed_frames = 0
|
||
|
||
# ------------------------------------------------------------------ push
|
||
def push(self, block: np.ndarray) -> None:
|
||
"""送入一块音频（float，(n, channels)）。"""
|
||
x = np.asarray(block, dtype=np.float64)
|
||
if x.ndim == 1:
|
||
x = x[:, None]
|
||
if x.shape[1] != self.channels:
|
||
raise ValueError("通道数与响度计配置不符")
|
||
self.processed_frames += x.shape[0]
|
||
if self._acc.size:
|
||
x = np.concatenate((self._acc, x), axis=0)
|
||
while x.shape[0] >= self.hop:
|
||
self._hop_energy(x[:self.hop])
|
||
x = x[self.hop:]
|
||
self._acc = x
|
||
|
||
def _hop_energy(self, hop_data: np.ndarray) -> None:
|
||
"""计算一个 hop 的加权均方（所有声道求和，声道权重取 1.0）。"""
|
||
spec = np.fft.rfft(hop_data, axis=0)
|
||
power = np.abs(spec) ** 2 # (bins, channels)
|
||
weighted = np.einsum("kc,k->", power, self._weight)
|
||
ms = float(weighted)
|
||
self._hops.append(ms)
|
||
if self._keep_hops:
|
||
self._all_hops.append(ms)
|
||
|
||
# ------------------------------------------------------------- readings
|
||
def _window_ms(self, hops: int) -> float:
|
||
n = min(hops, len(self._hops))
|
||
if n == 0:
|
||
return 0.0
|
||
vals = list(self._hops)[-n:]
|
||
return float(np.mean(vals))
|
||
|
||
@property
|
||
def momentary(self) -> float:
|
||
"""瞬时响度（400 ms 窗），LUFS。"""
|
||
return _lufs_from_ms(self._window_ms(self.MOMENTARY_HOPS))
|
||
|
||
@property
|
||
def short_term(self) -> float:
|
||
"""短时响度（3 s 窗），LUFS。"""
|
||
return _lufs_from_ms(self._window_ms(self.SHORT_TERM_HOPS))
|
||
|
||
@property
|
||
def integrated(self) -> float:
|
||
"""整体响度（带门限），LUFS。"""
|
||
if not self._keep_hops or len(self._all_hops) < self.MOMENTARY_HOPS:
|
||
return float("-inf")
|
||
# 400 ms 块（100 ms 步进 = 75% 重叠）的块能量序列
|
||
arr = np.asarray(self._all_hops, dtype=np.float64)
|
||
if arr.size < self.MOMENTARY_HOPS:
|
||
return float("-inf")
|
||
kernel = np.ones(self.MOMENTARY_HOPS) / self.MOMENTARY_HOPS
|
||
blocks = np.convolve(arr, kernel, mode="valid")
|
||
|
||
abs_gate = 10.0 ** ((_LUFS_ABS_GATE - _LUFS_OFFSET) / 10.0)
|
||
keep = blocks > abs_gate
|
||
if not np.any(keep):
|
||
return float("-inf")
|
||
mean1 = float(np.mean(blocks[keep]))
|
||
rel_gate = mean1 * (10.0 ** (_LUFS_REL_GATE / 10.0))
|
||
keep2 = blocks > max(rel_gate, abs_gate)
|
||
if not np.any(keep2):
|
||
return _lufs_from_ms(mean1)
|
||
return _lufs_from_ms(float(np.mean(blocks[keep2])))
|
||
|
||
def loudness_range_estimate(self) -> float:
|
||
"""LRA 粗估（10%~95% 分位的短时响度差）。
|
||
|
||
真正符合 EBU Tech 3342 的 LRA 需要 3 s 窗 + 10 s 分片，这里给出
|
||
一个足够稳定的工程近似值，用于提示动态范围。
|
||
"""
|
||
if not self._keep_hops:
|
||
return 0.0
|
||
arr = np.asarray(self._all_hops, dtype=np.float64)
|
||
if arr.size < self.SHORT_TERM_HOPS:
|
||
return 0.0
|
||
kernel = np.ones(self.SHORT_TERM_HOPS) / self.SHORT_TERM_HOPS
|
||
st = np.convolve(arr, kernel, mode="valid")
|
||
st = st[st > 10.0 ** ((_LUFS_ABS_GATE - _LUFS_OFFSET) / 10.0)]
|
||
if st.size < 2:
|
||
return 0.0
|
||
loud = _LUFS_OFFSET + 10.0 * np.log10(st)
|
||
p10, p95 = np.percentile(loud, 10.0), np.percentile(loud, 95.0)
|
||
return float(p95 - p10)
|
||
|
||
def reset(self) -> None:
|
||
self._acc = np.zeros((0, self.channels), dtype=np.float64)
|
||
self._hops.clear()
|
||
self._all_hops.clear()
|
||
self.processed_frames = 0
|
||
|
||
def flush(self) -> None:
|
||
"""把不足一个 hop 的尾部补零送算，保证短录音也有响度读数。"""
|
||
if self._acc.size:
|
||
tail = self._acc
|
||
pad = np.zeros((self.hop - tail.shape[0], self.channels))
|
||
self._hop_energy(np.concatenate((tail, pad), axis=0))
|
||
self._acc = np.zeros((0, self.channels))
|
||
|
||
|
||
# --------------------------------------------------------------- 实时电平表
|
||
@dataclass
|
||
class MeterSnapshot:
|
||
rms_db: list[float] = field(default_factory=list)
|
||
peak_db: list[float] = field(default_factory=list)
|
||
hold_db: list[float] = field(default_factory=list)
|
||
true_peak_db: list[float] = field(default_factory=list)
|
||
clipped: list[bool] = field(default_factory=list)
|
||
clip_count: int = 0
|
||
momentary_lufs: float = float("-inf")
|
||
short_term_lufs: float = float("-inf")
|
||
integrated_lufs: float = float("-inf")
|
||
|
||
|
||
class LevelMeter:
|
||
"""实时电平表：每块 RMS/峰值 + 峰值保持 + 削波锁存 + 真峰值 + 响度。"""
|
||
|
||
def __init__(self, samplerate: int, channels: int, *,
|
||
hold_decay_db_per_s: float = 12.0,
|
||
true_peak_enabled: bool = True,
|
||
loudness_enabled: bool = True):
|
||
self.samplerate = int(samplerate)
|
||
self.channels = int(channels)
|
||
self.hold_decay = float(hold_decay_db_per_s)
|
||
self.true_peak_enabled = bool(true_peak_enabled)
|
||
self.loudness_enabled = bool(loudness_enabled)
|
||
|
||
self._hold = np.full(self.channels, SILENCE_DBFS_FLOOR)
|
||
self._hold_timer = np.zeros(self.channels)
|
||
self._tp = np.zeros(self.channels)
|
||
self._clip = np.zeros(self.channels, dtype=bool)
|
||
self.clip_count = 0
|
||
self._last_block_time: float | None = None
|
||
self._tp_hop = max(1, int(round(0.1 * self.samplerate)))
|
||
self._acc: list[np.ndarray] = []
|
||
self._acc_frames = 0
|
||
self._preroll = np.zeros((0, self.channels))
|
||
|
||
self.loudness = LoudnessMeter(self.samplerate, self.channels) \
|
||
if self.loudness_enabled else None
|
||
self._rms_db = [SILENCE_DBFS_FLOOR] * self.channels
|
||
self._peak_db = [SILENCE_DBFS_FLOOR] * self.channels
|
||
|
||
# ------------------------------------------------------------------ push
|
||
def process(self, block: np.ndarray, *, now: float | None = None) -> None:
|
||
x = np.asarray(block, dtype=np.float64)
|
||
if x.ndim == 1:
|
||
x = x[:, None]
|
||
if x.size == 0:
|
||
return
|
||
|
||
# 1) 逐块 RMS / 峰值（更新快，界面刷新率高）
|
||
self._rms_db = [lin_to_db(v) for v in rms(x, axis=0)]
|
||
self._peak_db = [lin_to_db(v) for v in peak(x, axis=0)]
|
||
|
||
# 2) 峰值保持（带自动回落，避免"钉住"不动）
|
||
if now is None:
|
||
import time
|
||
now = time.monotonic()
|
||
dt = 0.0 if self._last_block_time is None else max(0.0, now - self._last_block_time)
|
||
self._last_block_time = now
|
||
self._hold = np.maximum(self._hold - self.hold_decay * dt,
|
||
np.asarray(self._peak_db, dtype=np.float64))
|
||
|
||
# 3) 削波锁存
|
||
over = np.any(np.abs(x) >= 0.99999, axis=0)
|
||
self._clip |= over
|
||
self.clip_count += int(np.count_nonzero(np.abs(x) >= 0.99999))
|
||
|
||
# 4) 响度
|
||
if self.loudness is not None:
|
||
self.loudness.push(x)
|
||
|
||
# 5) 真峰值：按 100 ms 累积（含 12 样本前导以避免块边界低估）
|
||
if self.true_peak_enabled:
|
||
self._acc.append(x)
|
||
self._acc_frames += x.shape[0]
|
||
while self._acc_frames >= self._tp_hop:
|
||
merged = np.concatenate(self._acc, axis=0)
|
||
window = merged[:self._tp_hop]
|
||
rest = merged[self._tp_hop:]
|
||
seg = np.concatenate((self._preroll, window), axis=0) \
|
||
if self._preroll.size else window
|
||
tp = true_peak(seg)
|
||
# 只统计新增部分的真峰值（近似：整段最大值）
|
||
self._tp = np.maximum(self._tp, tp)
|
||
self._preroll = window[-12:] if window.shape[0] >= 12 else window
|
||
self._acc = [rest] if rest.size else []
|
||
self._acc_frames = rest.shape[0]
|
||
|
||
# ------------------------------------------------------------- readings
|
||
def snapshot(self) -> MeterSnapshot:
|
||
return MeterSnapshot(
|
||
rms_db=[float(v) for v in self._rms_db],
|
||
peak_db=list(self._peak_db),
|
||
hold_db=[float(v) for v in self._hold],
|
||
true_peak_db=[lin_to_db(v) for v in self._tp],
|
||
clipped=[bool(v) for v in self._clip],
|
||
clip_count=self.clip_count,
|
||
momentary_lufs=self.loudness.momentary if self.loudness else float("-inf"),
|
||
short_term_lufs=self.loudness.short_term if self.loudness else float("-inf"),
|
||
integrated_lufs=self.loudness.integrated if self.loudness else float("-inf"),
|
||
)
|
||
|
||
def reset_clip(self) -> None:
|
||
self._clip[:] = False
|
||
|
||
def reset(self) -> None:
|
||
self._hold[:] = SILENCE_DBFS_FLOOR
|
||
self._tp[:] = 0.0
|
||
self._clip[:] = False
|
||
self.clip_count = 0
|
||
self._acc = []
|
||
self._acc_frames = 0
|
||
self._preroll = np.zeros((0, self.channels))
|
||
self._last_block_time = None
|
||
if self.loudness is not None:
|
||
self.loudness.reset()
|
||
|
||
|
||
# --------------------------------------------------------------- 离线分析
|
||
def _noise_floor(block_rms: np.ndarray) -> tuple[list[float], bool]:
|
||
"""从 50 ms 块 RMS 序列估计本底噪声。
|
||
|
||
取安静段的 10% 分位作为本底，但要求整体动态范围（90% 分位 / 10% 分位）
|
||
至少 10 dB——否则说明录音里根本没有"安静段"（例如全程持续发声的
|
||
测试信号），此时返回不可测，而不是编造一个数字。
|
||
"""
|
||
if block_rms.size < 4:
|
||
return [], False
|
||
p10 = np.percentile(block_rms, 10.0, axis=0)
|
||
p90 = np.percentile(block_rms, 90.0, axis=0)
|
||
if np.any(p10 <= 0) or np.any(p90 <= 0):
|
||
return [], False
|
||
if np.any(20.0 * np.log10(p90 / p10) < 10.0):
|
||
return [], False
|
||
return [round(float(v), 2) for v in np.atleast_1d(dbfs(p10))], True
|
||
|
||
|
||
def analyze_data(data: np.ndarray, samplerate: int, *,
|
||
true_peak_factor: int = 4) -> dict:
|
||
"""对一段音频做全面体检，返回可直接写入元数据的字典。"""
|
||
x = np.asarray(data, dtype=np.float64)
|
||
if x.ndim == 1:
|
||
x = x[:, None]
|
||
n, ch = x.shape
|
||
duration = n / float(samplerate) if samplerate else 0.0
|
||
|
||
pk = peak(x, axis=0)
|
||
pk_db = dbfs(pk)
|
||
tp = true_peak(x, factor=true_peak_factor) if n else np.zeros(ch)
|
||
tp_db = dbfs(tp)
|
||
rms_db = dbfs(rms(x, axis=0))
|
||
|
||
# 直流偏移
|
||
dc = np.mean(x, axis=0) if n else np.zeros(ch)
|
||
|
||
# 噪声本底：只在检测到安静段时才给出
|
||
block = max(1, int(0.05 * samplerate))
|
||
usable = (n // block) * block
|
||
if usable >= block:
|
||
blocks = x[:usable].reshape(-1, block, ch)
|
||
blk_rms = np.sqrt(np.mean(np.square(blocks), axis=1)) # (nblocks, ch)
|
||
floor_db, floor_ok = _noise_floor(blk_rms)
|
||
else:
|
||
floor_db, floor_ok = [], False
|
||
|
||
# 削波统计
|
||
clip_counts = np.count_nonzero(np.abs(x) >= 0.99999, axis=0)
|
||
|
||
# 静音占比（低于 -60 dBFS 的采样比例）
|
||
quiet = float(np.mean(np.abs(x) < db_to_lin(-60.0))) if n else 1.0
|
||
|
||
# 响度
|
||
lm = LoudnessMeter(samplerate, ch)
|
||
step = 1 << 16
|
||
for i in range(0, n, step):
|
||
lm.push(x[i:i + step])
|
||
lm.flush()
|
||
|
||
return {
|
||
"samplerate": int(samplerate),
|
||
"channels": int(ch),
|
||
"frames": int(n),
|
||
"duration": duration,
|
||
"peak_dbfs": [round(float(v), 3) for v in np.atleast_1d(pk_db)],
|
||
"true_peak_dbtp": [round(float(v), 3) for v in np.atleast_1d(tp_db)],
|
||
"rms_dbfs": [round(float(v), 3) for v in np.atleast_1d(rms_db)],
|
||
"integrated_lufs": round(float(lm.integrated), 3)
|
||
if math.isfinite(lm.integrated) else None,
|
||
"short_term_max_lufs": None,
|
||
"loudness_range_lu": round(float(lm.loudness_range_estimate()), 3),
|
||
"dc_offset": [round(float(v), 8) for v in np.atleast_1d(dc)],
|
||
"noise_floor_dbfs": floor_db,
|
||
"noise_floor_available": floor_ok,
|
||
"clipped_samples": [int(v) for v in np.atleast_1d(clip_counts)],
|
||
"clipped_total": int(np.sum(clip_counts)),
|
||
"silence_ratio": round(quiet, 4),
|
||
}
|
||
|
||
|
||
def analyze_file(path: str, *, max_seconds: float | None = None) -> dict:
|
||
"""对 WAV 文件做体检（分块读取，适合超长录音）。"""
|
||
from .wavfile import WavReader
|
||
|
||
with WavReader(path) as r:
|
||
info = r.info()
|
||
limit = r.frames if max_seconds is None else min(
|
||
r.frames, int(max_seconds * r.samplerate))
|
||
lm = LoudnessMeter(r.samplerate, r.channels)
|
||
pk = np.zeros(r.channels)
|
||
tp = np.zeros(r.channels)
|
||
rms_acc = np.zeros(r.channels, dtype=np.float64)
|
||
dc_acc = np.zeros(r.channels, dtype=np.float64)
|
||
clips = np.zeros(r.channels, dtype=np.int64)
|
||
count = 0
|
||
blk_rms: list[np.ndarray] = []
|
||
preroll = np.zeros((0, r.channels), dtype=np.float64)
|
||
while count < limit:
|
||
block = r.read(count, min(1 << 16, limit - count))
|
||
if block.size == 0:
|
||
break
|
||
x64 = block # 已是 float64
|
||
pk = np.maximum(pk, np.max(np.abs(x64), axis=0))
|
||
# 真峰值：带 32 样本前导，避免块边界漏掉采样点之间的过冲
|
||
seg = np.concatenate((preroll, x64), axis=0) if preroll.size else x64
|
||
if seg.shape[0] > 128:
|
||
tp = np.maximum(tp, true_peak(seg))
|
||
else:
|
||
tp = np.maximum(tp, np.max(np.abs(seg), axis=0))
|
||
preroll = x64[-32:] if x64.shape[0] >= 32 else x64
|
||
rms_acc += np.sum(np.square(x64), axis=0)
|
||
dc_acc += np.sum(x64, axis=0)
|
||
clips += np.count_nonzero(np.abs(x64) >= 0.99999, axis=0)
|
||
lm.push(x64)
|
||
sub = max(1, int(0.05 * r.samplerate))
|
||
usable = (x64.shape[0] // sub) * sub
|
||
if usable:
|
||
blk_rms.append(np.sqrt(np.mean(
|
||
np.square(x64[:usable].reshape(-1, sub, r.channels)), axis=1)))
|
||
count += block.shape[0]
|
||
lm.flush()
|
||
|
||
duration = count / r.samplerate if r.samplerate else 0.0
|
||
floor_db, floor_ok = _noise_floor(np.concatenate(blk_rms, axis=0)) \
|
||
if blk_rms else ([], False)
|
||
out = dict(info)
|
||
out.update({
|
||
"analyzed_frames": count,
|
||
"analyzed_duration": duration,
|
||
"peak_dbfs": [round(float(v), 3) for v in np.atleast_1d(dbfs(pk))],
|
||
"true_peak_dbtp": [round(float(v), 3) for v in np.atleast_1d(dbfs(tp))],
|
||
"rms_dbfs": [round(float(v), 3) for v in np.atleast_1d(
|
||
dbfs(np.sqrt(rms_acc / max(1, count))))],
|
||
"integrated_lufs": round(float(lm.integrated), 3)
|
||
if math.isfinite(lm.integrated) else None,
|
||
"loudness_range_lu": round(float(lm.loudness_range_estimate()), 3),
|
||
"dc_offset": [round(float(v), 8) for v in dc_acc / max(1, count)],
|
||
"noise_floor_dbfs": floor_db,
|
||
"noise_floor_available": floor_ok,
|
||
"clipped_samples": [int(v) for v in clips],
|
||
"clipped_total": int(clips.sum()),
|
||
})
|
||
return out
|
||
|
||
|
||
# ------------------------------------------------------------- 离线处理链
|
||
def remove_dc(x: np.ndarray, samplerate: int | None = None,
|
||
window_seconds: float = 1.0) -> np.ndarray:
|
||
"""去除直流与极低频漂移。
|
||
|
||
给出 ``samplerate`` 时减去**滑动平均**（默认 1 秒窗）而不是全局均值：
|
||
录音中途如果直流电平发生漂移（设备热漂移、某些 USB 声卡的固有偏移），
|
||
全局均值法只会把漂移原样留成低频噪声，滑动平均才能真正跟掉它。
|
||
对 1 秒窗而言，2 Hz 以上的频响基本平坦，不会动到有用信号。
|
||
"""
|
||
arr = np.asarray(x, dtype=np.float64)
|
||
if arr.size == 0:
|
||
return arr
|
||
if arr.ndim == 1:
|
||
arr = arr[:, None]
|
||
n, ch = arr.shape
|
||
w = int(window_seconds * samplerate) if samplerate else 0
|
||
if w < 3 or w >= n:
|
||
return arr - np.mean(arr, axis=0, keepdims=True)
|
||
cum = np.cumsum(np.concatenate([np.zeros((1, ch)), arr], axis=0), axis=0)
|
||
pad = w // 2
|
||
idx = np.arange(n)
|
||
lo = np.clip(idx - pad, 0, n)
|
||
hi = np.clip(idx - pad + w, 0, n)
|
||
count = np.maximum(1, hi - lo)
|
||
avg = (cum[hi] - cum[lo]) / count[:, None]
|
||
return arr - avg
|
||
|
||
|
||
def highpass_offline(x: np.ndarray, samplerate: int, cutoff_hz: float) -> np.ndarray:
|
||
h = design_highpass_fir(cutoff_hz, samplerate)
|
||
f = FIRFilter(h, x.shape[1] if x.ndim > 1 else 1)
|
||
y = f.process(x)
|
||
# 补偿滤波群延迟，保持时间轴对齐
|
||
d = f.latency_samples
|
||
return y[d:] if d else y
|
||
|
||
|
||
def trim_silence(x: np.ndarray, samplerate: int, *,
|
||
threshold_dbfs: float = -50.0,
|
||
min_silence: float = 0.35,
|
||
pad: float = 0.15) -> tuple[np.ndarray, dict]:
|
||
"""裁掉首尾静音，返回 ``(新数据, 裁剪信息)``。
|
||
|
||
只有当首/尾静音长于 ``min_silence`` 时才裁剪，避免把短促的气口误伤；
|
||
裁剪后仍保留 ``pad`` 秒的呼吸空间，防止削掉起音与尾音的自然衰减。
|
||
"""
|
||
if x.size == 0:
|
||
return x, {"trimmed_seconds": 0.0}
|
||
thr = db_to_lin(threshold_dbfs)
|
||
win = max(1, int(0.01 * samplerate))
|
||
usable = (x.shape[0] // win) * win
|
||
if usable < win:
|
||
return x, {"trimmed_seconds": 0.0}
|
||
env = np.max(np.abs(x[:usable]), axis=1).reshape(-1, win).max(axis=1)
|
||
idx = np.flatnonzero(env >= thr)
|
||
if idx.size == 0:
|
||
return x, {"trimmed_seconds": 0.0, "all_silent": True}
|
||
|
||
pad_blocks = int(pad * samplerate / win)
|
||
min_blocks = max(1, int(min_silence * samplerate / win))
|
||
start_block = int(idx[0])
|
||
end_block = int(idx[-1]) + 1
|
||
# 首部静音不够长就不动它
|
||
if start_block < min_blocks:
|
||
start_block = 0
|
||
else:
|
||
start_block = max(0, start_block - pad_blocks)
|
||
# 尾部同理
|
||
if env.size - end_block < min_blocks:
|
||
end_block = env.size
|
||
else:
|
||
end_block = min(env.size, end_block + pad_blocks)
|
||
|
||
start = start_block * win
|
||
end = min(x.shape[0], end_block * win)
|
||
trimmed = x[start:end]
|
||
if trimmed.size == 0:
|
||
return x, {"trimmed_seconds": 0.0}
|
||
return trimmed, {
|
||
"trimmed_seconds": round((x.shape[0] - trimmed.shape[0]) / samplerate, 3),
|
||
"lead_seconds": round(start / samplerate, 3),
|
||
"tail_seconds": round((x.shape[0] - end) / samplerate, 3),
|
||
"min_silence": min_silence,
|
||
"threshold_dbfs": threshold_dbfs,
|
||
}
|
||
|
||
|
||
def fade_edges(x: np.ndarray, samplerate: int, *,
|
||
fade_in: float = 0.005, fade_out: float = 0.02) -> np.ndarray:
|
||
y = np.array(x, copy=True)
|
||
n = y.shape[0]
|
||
for dur, head in ((fade_in, True), (fade_out, False)):
|
||
k = min(n, int(dur * samplerate))
|
||
if k <= 1:
|
||
continue
|
||
ramp = np.linspace(0.0, 1.0, k, dtype=np.float64)
|
||
if head:
|
||
y[:k] *= ramp[:, None]
|
||
else:
|
||
y[n - k:] *= ramp[::-1, None]
|
||
return y
|
||
|
||
|
||
def normalize_peak(x: np.ndarray, target_dbfs: float = -1.0,
|
||
*, true_peak_mode: bool = True,
|
||
samplerate: int | None = None) -> tuple[np.ndarray, float]:
|
||
"""按峰值归一化（默认按真峰值，避免转码后过载）。"""
|
||
if x.size == 0:
|
||
return x, 0.0
|
||
cur = float(np.max(true_peak(x))) if (true_peak_mode and x.shape[0] > 64) \
|
||
else float(peak(x))
|
||
if cur <= 1e-12:
|
||
return x, 0.0
|
||
gain = db_to_lin(target_dbfs) / cur
|
||
return x * gain, float(20.0 * math.log10(gain))
|
||
|
||
|
||
def normalize_loudness(x: np.ndarray, samplerate: int, target_lufs: float,
|
||
*, max_gain_db: float = 30.0,
|
||
true_peak_ceiling_dbfs: float = -1.0,
|
||
) -> tuple[np.ndarray, dict]:
|
||
"""按整体响度（BS.1770）归一化，并保证真峰值不超过上限。"""
|
||
if x.size == 0:
|
||
return x, {"gain_db": 0.0}
|
||
ch = x.shape[1] if x.ndim > 1 else 1
|
||
lm = LoudnessMeter(samplerate, ch)
|
||
step = 1 << 16
|
||
for i in range(0, x.shape[0], step):
|
||
lm.push(x[i:i + step])
|
||
lm.flush()
|
||
current = float(lm.integrated)
|
||
if not math.isfinite(current):
|
||
return x, {"gain_db": 0.0, "note": "无法测定响度，未做处理"}
|
||
gain_db = target_lufs - current
|
||
limited = False
|
||
if gain_db > max_gain_db:
|
||
gain_db, limited = max_gain_db, True
|
||
if gain_db < -max_gain_db:
|
||
gain_db, limited = -max_gain_db, True
|
||
|
||
y = x * db_to_lin(gain_db)
|
||
# 峰值保护：若真峰值越界则回退增益
|
||
tp = float(np.max(true_peak(y))) if y.shape[0] > 64 else float(peak(y))
|
||
ceiling = db_to_lin(true_peak_ceiling_dbfs)
|
||
if tp > ceiling and tp > 0:
|
||
extra = ceiling / tp
|
||
y = y * extra
|
||
gain_db += 20.0 * math.log10(extra)
|
||
return y, {
|
||
"gain_db": round(float(gain_db), 3),
|
||
"measured_lufs": round(current, 3),
|
||
"target_lufs": float(target_lufs),
|
||
"gain_limited": limited,
|
||
}
|
||
|
||
|
||
def noise_gate(x: np.ndarray, samplerate: int, *, threshold_dbfs: float = -60.0,
|
||
attack: float = 0.005, release: float = 0.12,
|
||
hold: float = 0.1) -> tuple[np.ndarray, dict]:
|
||
"""带包络跟随的软门限（仅在用户显式开启时使用）。"""
|
||
if x.size == 0:
|
||
return x, {"gated_ratio": 0.0}
|
||
thr = db_to_lin(threshold_dbfs)
|
||
win = max(1, int(0.005 * samplerate))
|
||
usable = (x.shape[0] // win) * win
|
||
env = np.max(np.abs(x[:usable]), axis=1).reshape(-1, win).max(axis=1)
|
||
open_mask = env >= thr
|
||
# 保持（hold）+ 攻放（attack/release）平滑，避免产生"抽气"感
|
||
hop = win / samplerate
|
||
hold_frames = max(1, int(hold / hop))
|
||
atk = max(1, int(attack / hop))
|
||
rel = max(1, int(release / hop))
|
||
gain = np.zeros(env.size)
|
||
g = 0.0
|
||
counter = 0
|
||
for i in range(env.size):
|
||
if open_mask[i]:
|
||
counter = hold_frames
|
||
target = 1.0
|
||
elif counter > 0:
|
||
counter -= 1
|
||
target = 1.0
|
||
else:
|
||
target = 0.0
|
||
step = 1.0 / (atk if target > g else rel)
|
||
g += np.clip(target - g, -step, step)
|
||
gain[i] = g
|
||
# 上采样到采样级并做线性插值，避免阶梯噪声
|
||
up = np.interp(np.arange(x.shape[0]),
|
||
np.arange(env.size) * win, gain)
|
||
y = x * up[:, None]
|
||
return y, {
|
||
"threshold_dbfs": threshold_dbfs,
|
||
"gated_ratio": round(float(np.mean(up < 0.5)), 4),
|
||
}
|
||
|
||
|
||
def mixdown_mono(x: np.ndarray) -> np.ndarray:
|
||
if x.ndim == 1 or x.shape[1] == 1:
|
||
return x
|
||
return np.mean(x, axis=1, keepdims=True)
|