ffmpeg -i input.wav -af \
"loudnorm=I=-23:TP=-3:LRA=7, \
highpass=f=85, \
equalizer=f=250:t=q:w=1.5:g=-5, \
equalizer=f=450:t=q:w=1.2:g=-3, \
equalizer=f=4500:t=q:w=1:g=4, \
equalizer=f=8000:t=q:w=1:g=2, \
acompressor=threshold=-22dB:ratio=2.5:attack=8:release=120:makeup=3, \
loudnorm=I=-16:TP=-1.5:LRA=11" \
output.wav

uv run analyze_audio.py episode.wav

#!/usr/bin/env python3
# /// script
# requires-python = ">=3.9"
# dependencies = [
#     "numpy",
#     "scipy",
# ]
# ///
"""
analyze_audio.py — 教堂讲道录音质量诊断脚本

检测内容：
  1. 削波(clipping)检测 —— 找连续平顶采样，比单看峰值电平准确
  2. 频段能量分布 —— 判断是否存在"空洞/发闷"（150-1000Hz 占比过高，
     3000Hz+ 占比过低）
  3. 混响衰减粗估(RT60-ish) —— 基于语音瞬态后的能量衰减斜率
  4. 立体声 L/R 相关性 —— 排除梳状滤波/相位抵消问题
  5. 响度统计(RMS/峰值)

依赖：只有 numpy + scipy。用 scipy.io.wavfile 读 wav（纯 Python 解析文件头，
不需要编译好的二进制库），刻意不用 soundfile/librosa —— 前者依赖 libsndfile
这个原生二进制，在较新/较冷门的 Python 版本上容易出现"装得上但打不开文件"
的兼容性问题；后者依赖 numba/LLVM，在树莓派等 ARM 设备上装起来容易踩坑。
限制：只能读 wav（这个诊断场景下够用——DR-05X 录出来的本来就是 wav）。

用法：
    python3 analyze_audio.py episode1.wav [episode2.wav ...]
    python3 analyze_audio.py episode1.wav --json          # 输出JSON，便于pipeline解析
    uv run analyze_audio.py episode1.wav                  # 免装依赖，见脚本头部声明

作为库使用：
    from analyze_audio import analyze_audio
    report = analyze_audio("episode1.wav")
"""

import sys
import json
import argparse
import itertools

import numpy as np
from scipy.io import wavfile
from scipy.signal import stft as scipy_stft


def _read_wav_normalized(path):
    """读取 wav 并归一化到 [-1, 1] 浮点区间，返回 (data, sr)。
    data 形状为 (samples,) 或 (samples, channels)。
    """
    sr, data = wavfile.read(path)
    if np.issubdtype(data.dtype, np.floating):
        return data.astype(np.float64), sr
    if data.dtype == np.uint8:
        # 8-bit wav 是无符号的，中心在128
        return (data.astype(np.float64) - 128.0) / 128.0, sr
    # 有符号整数类型（int16/int32等）：按该dtype的满量程归一化
    max_val = float(np.iinfo(data.dtype).max) + 1.0
    return data.astype(np.float64) / max_val, sr


BANDS = [
    (20, 150, "20-150Hz"),
    (150, 400, "150-400Hz"),
    (400, 1000, "400-1kHz"),
    (1000, 3000, "1k-3kHz"),
    (3000, 6000, "3k-6kHz"),
    (6000, 12000, "6k-12kHz"),
    (12000, 20000, "12k-20kHz"),
]

# 判断"空洞/发闷"的经验阈值，来自对这个讲道系列前两条样本的实测
# （150-1000Hz 占比 > 65% 且 3000Hz+ 占比 < 3% 时，主观听感基本必发闷）
BOXY_LOWMID_THRESHOLD = 0.65
THIN_HIGH_THRESHOLD = 0.03


def _detect_clipping(y, thresh=0.999, min_run=3):
    """检测连续平顶采样（真实数字削波的特征），而不是只看单个峰值。
    返回 (削波采样总数, 最长连续平顶采样数, 削波片段数)。
    """
    if y.ndim == 1:
        y = y[:, None]
    total_clipped = 0
    longest_run = 0
    n_events = 0
    for ch in range(y.shape[1]):
        flat = np.abs(y[:, ch]) >= thresh
        runs = [len(list(g)) for k, g in itertools.groupby(flat) if k]
        real_runs = [r for r in runs if r >= min_run]
        total_clipped += sum(real_runs)
        n_events += len(real_runs)
        if real_runs:
            longest_run = max(longest_run, max(real_runs))
    return total_clipped, longest_run, n_events


def _band_energy(y_mono, sr, n_fft=4096):
    """返回各频段的平均幅度及占比。"""
    freqs, times, Zxx = scipy_stft(y_mono, fs=sr, nperseg=n_fft)
    mag = np.abs(Zxx)
    avg_spec = np.mean(mag, axis=1)

    results = []
    for lo, hi, label in BANDS:
        idx = np.where((freqs >= lo) & (freqs < hi))[0]
        e = float(np.mean(avg_spec[idx])) if len(idx) else 0.0
        results.append({"label": label, "lo": lo, "hi": hi, "energy": e})

    total = sum(b["energy"] for b in results) or 1e-12
    for b in results:
        b["pct"] = 100.0 * b["energy"] / total
    return results


def _spectral_centroid(y_mono, sr, n_fft=4096):
    freqs, times, Zxx = scipy_stft(y_mono, fs=sr, nperseg=n_fft)
    mag = np.abs(Zxx)
    avg_spec = np.mean(mag, axis=1)
    if avg_spec.sum() == 0:
        return 0.0
    return float(np.sum(freqs * avg_spec) / np.sum(avg_spec))


def _rt60_estimate(y_mono, sr, frame=512, hop=128, drop_db=15, max_events=300, rel_thresh_db=10):
    """粗估混响衰减：找语音瞬态峰值，测量衰减 drop_db 所需时间，外推到60dB。
    注意：这是基于能量包络的经验估算，不是声学测量级的RT60，仅用于趋势判断。

    峰值检测阈值是相对该文件自身最大电平动态确定的（max_db - rel_thresh_db），
    而不是写死的绝对dB值——不然遇到整体电平偏低的文件（比如某些降噪工具处理后
    的输出），会因为所有峰值都低于固定阈值而一个瞬态都找不到，导致估算结果
    误报"None"（不是真的没有混响，是检测器没找到东西可测）。
    """
    n_frames = 1 + (len(y_mono) - frame) // hop
    if n_frames <= 0:
        return None, 0
    rms = np.array([
        np.sqrt(np.mean(y_mono[i * hop: i * hop + frame] ** 2) + 1e-12)
        for i in range(n_frames)
    ])
    db = 20 * np.log10(rms + 1e-9)
    times = np.arange(n_frames) * hop / sr

    thresh = db.max() - rel_thresh_db
    peaks = []
    for i in range(2, len(db) - 2):
        if (db[i] > thresh and db[i] >= db[i - 1] and db[i] >= db[i + 1]
                and db[i] > db[i - 2] and db[i] > db[i + 2]):
            peaks.append(i)

    decays = []
    max_lookahead = int(0.4 * sr / hop)
    for p in peaks[:max_events]:
        peak_val = db[p]
        target = peak_val - drop_db
        for j in range(p, min(p + max_lookahead, len(db))):
            if db[j] < target:
                t_decay = times[j] - times[p]
                if t_decay > 0.02:
                    decays.append(60 / (drop_db / t_decay))
                break

    decays = np.array(decays)
    decays = decays[(decays > 0.05) & (decays < 3)]
    if len(decays) == 0:
        return None, 0
    return float(np.median(decays)), len(decays)


def analyze_audio(path, n_fft=4096):
    y, sr = _read_wav_normalized(path)
    is_stereo = y.ndim == 2 and y.shape[1] >= 2
    y_mono = y.mean(axis=1) if y.ndim == 2 else y

    duration = len(y_mono) / sr
    peak_dbfs = 20 * np.log10(np.max(np.abs(y)) + 1e-12)
    rms_db = 20 * np.log10(np.sqrt(np.mean(y_mono ** 2)) + 1e-12)

    clipped_samples, longest_run, n_clip_events = _detect_clipping(y)

    bands = _band_energy(y_mono, sr, n_fft=n_fft)
    centroid = _spectral_centroid(y_mono, sr, n_fft=n_fft)
    rt60, rt60_n = _rt60_estimate(y_mono, sr)

    lowmid_pct = sum(b["pct"] for b in bands if b["lo"] >= 150 and b["hi"] <= 1000) / 100
    high_pct = sum(b["pct"] for b in bands if b["lo"] >= 3000) / 100

    corr = None
    side_mid_ratio = None
    if is_stereo:
        L, R = y[:, 0], y[:, 1]
        corr = float(np.corrcoef(L, R)[0, 1])
        mid, side = (L + R) / 2, (L - R) / 2
        mid_e, side_e = np.mean(mid ** 2), np.mean(side ** 2)
        side_mid_ratio = float(side_e / mid_e) if mid_e > 0 else None

    warnings = []
    if clipped_samples > 0:
        warnings.append(
            f"检测到真实削波：{clipped_samples}个采样点，{n_clip_events}处片段，"
            f"最长连续{longest_run}个采样。这部分失真后期无法修复，建议检查录制增益。"
        )
    if lowmid_pct >= BOXY_LOWMID_THRESHOLD and high_pct <= THIN_HIGH_THRESHOLD:
        warnings.append(
            f"频谱失衡：150-1000Hz占比{lowmid_pct*100:.1f}%，3000Hz+占比仅{high_pct*100:.1f}%，"
            "符合'空洞/发闷'特征（麦克风距离过远/混响过重）。"
        )
    if rt60 is not None and rt60 > 0.8:
        warnings.append(f"混响衰减估算约{rt60:.2f}秒，偏长（语音清晰度理想值在0.3-0.6秒）。")
    if corr is not None and corr < 0.5:
        warnings.append(f"左右声道相关性较低({corr:.2f})，可能存在相位抵消/梳状滤波问题。")

    return {
        "file": path,
        "duration_sec": duration,
        "sample_rate": sr,
        "stereo": is_stereo,
        "peak_dbfs": peak_dbfs,
        "rms_dbfs": rms_db,
        "clipping": {
            "clipped_samples": clipped_samples,
            "longest_run": longest_run,
            "n_events": n_clip_events,
        },
        "spectral_centroid_hz": centroid,
        "bands": bands,
        "lowmid_pct": lowmid_pct * 100,
        "high_pct": high_pct * 100,
        "rt60_estimate_sec": rt60,
        "rt60_n_samples": rt60_n,
        "lr_correlation": corr,
        "side_mid_ratio": side_mid_ratio,
        "warnings": warnings,
    }


def print_report(report):
    print(f"\n=== {report['file']} ===")
    print(f"时长: {report['duration_sec']:.1f}s  采样率: {report['sample_rate']}Hz  "
          f"{'立体声' if report['stereo'] else '单声道'}")
    print(f"峰值: {report['peak_dbfs']:.1f} dBFS   平均电平: {report['rms_dbfs']:.1f} dBFS")
    print(f"频谱质心: {report['spectral_centroid_hz']:.0f} Hz")
    print("频段分布:")
    for b in report["bands"]:
        bar = "#" * int(b["pct"] / 2)
        print(f"  {b['label']:>10}: {b['pct']:5.1f}%  {bar}")
    if report["rt60_estimate_sec"] is not None:
        print(f"混响衰减估算: {report['rt60_estimate_sec']:.2f}s (基于{report['rt60_n_samples']}个瞬态样本)")
    if report["lr_correlation"] is not None:
        print(f"L/R相关性: {report['lr_correlation']:.3f}  side/mid能量比: {report['side_mid_ratio']:.3f}")
    c = report["clipping"]
    print(f"削波检测: {c['clipped_samples']}个采样点, {c['n_events']}处片段, 最长连续{c['longest_run']}")

    if report["warnings"]:
        print("\n⚠ 警告:")
        for w in report["warnings"]:
            print(f"  - {w}")
    else:
        print("\n✓ 未发现明显问题")


def main():
    parser = argparse.ArgumentParser(description="教堂讲道录音频谱诊断")
    parser.add_argument("files", nargs="+", help="待分析的音频文件（wav格式）")
    parser.add_argument("--json", action="store_true", help="以JSON输出，便于pipeline脚本解析")
    args = parser.parse_args()

    reports = [analyze_audio(f) for f in args.files]

    if args.json:
        print(json.dumps(reports, ensure_ascii=False, indent=2))
    else:
        for r in reports:
            print_report(r)


if __name__ == "__main__":
    main()