#!/bin/bash
# 响度校准: HRTF 卷积会把能量摊薄, 算准该补多少 dB 再补, 别拍脑袋
set -eu
cd "/home/lou/桌面/工作区/实验/空间音频"
SOFA=/usr/share/libmysofa/MIT_KEMAR_normal_pinna.sofa

echo "===== 1. 量现有两版的响度 ====="
python3 - <<'PY'
import math
import struct
import wave


def stat(path: str) -> tuple[float, float]:
    """返回 (峰值, 均方根), 都归一化到 1.0。"""
    with wave.open(path) as w:
        channels, frames = w.getnchannels(), w.getnframes()
        raw = w.readframes(frames)
    data = struct.unpack("<%dh" % (len(raw) // 2), raw)
    peak = max(abs(x) for x in data) / 32768.0
    rms = (sum(x * x for x in data) / len(data)) ** 0.5 / 32768.0
    return peak, rms


peak_a, rms_a = stat("A-普通下混.wav")
peak_b, rms_b = stat("B-HRTF环绕.wav")
need = 20 * math.log10(rms_a / rms_b)
print("  A 普通下混   峰值 %.3f (%.1f dBFS)   均方根 %.4f (%.1f dBFS)"
      % (peak_a, 20 * math.log10(peak_a), rms_a, 20 * math.log10(rms_a)))
print("  B HRTF       峰值 %.3f (%.1f dBFS)   均方根 %.4f (%.1f dBFS)"
      % (peak_b, 20 * math.log10(peak_b), rms_b, 20 * math.log10(rms_b)))
print("  => 需补 %.1f dB; 补后 B 峰值约 %.1f dBFS" % (need, 20 * math.log10(peak_b) + need))
PY

echo
echo "===== 2. 按算出的量重渲染 B 版 ====="
GAIN=$(python3 -c "
import math, struct, wave
def rms(p):
    with wave.open(p) as w:
        c, n = w.getnchannels(), w.getnframes()
        raw = w.readframes(n)
    d = struct.unpack('<%dh' % (len(raw)//2), raw)
    return (sum(x*x for x in d)/len(d))**0.5
print('%.1f' % (20*math.log10(rms('A-普通下混.wav')/rms('B-HRTF环绕.wav'))))
")
echo "  gain = ${GAIN} dB"
ffmpeg -y -hide_banner -loglevel warning -i 51-test.wav \
  -af "sofalizer=sofa=$SOFA:type=freq:radius=1.2:gain=${GAIN}" \
  -c:a pcm_s16le "B-HRTF环绕.wav"

echo
echo "===== 3. 复测(峰值别顶到 0dBFS, 否则削波) ====="
python3 - <<'PY'
import math
import struct
import wave

for f in ("A-普通下混.wav", "B-HRTF环绕.wav"):
    with wave.open(f) as w:
        channels, frames = w.getnchannels(), w.getnframes()
        raw = w.readframes(frames)
    data = struct.unpack("<%dh" % (len(raw) // 2), raw)
    peak = max(abs(x) for x in data) / 32768.0
    rms = (sum(x * x for x in data) / len(data)) ** 0.5 / 32768.0
    clip = sum(1 for x in data if abs(x) >= 32767)
    print("  %-18s 峰值 %.1f dBFS   均方根 %.1f dBFS   削波点 %d"
          % (f, 20 * math.log10(peak), 20 * math.log10(rms), clip))
PY