#!/usr/bin/env python3 """HOW LOUD IS THE DISC? ROADMAP P6, the item FINDINGS 66.3 reopened. python3 tools/analysis/35_audio_level.py [--streams 00000-00201] [--json out] FINDINGS 66 asked MAME's MSM6258 which decoder it is and got four axes back. The one with a bill attached is the CLAMP: the chip's accumulator is **10 bits** and it clamps INSIDE the recursion, so the reachable set of reconstructed samples is [-512, 511] in the 12-bit units everything in this project counts in -- a quarter of the 12-bit word `adpcm.py` used to clamp at. `pack.py` hands the encoder `s16 >> 4`, i.e. it maps the disc's full scale onto the 12-bit word, and 66.3 measured the Singe window peaking at **435 of 511**. That fit with 1.4 dB to spare, and it fit BY ACCIDENT: the window is a -13.4 dBFS passage. Any passage more than 1.4 dB louder does not merely distort at the top, it drives the predictor -- a clamped accumulator is a WRONG STATE that the next nibble is applied to, so the error outlives the loud sample. So the level cannot be chosen from the ten seconds this project gates on. It has to be chosen from the loudest thing the game will ever play, and this file measures that: every stream of the unique scene footage -- `00000`-`00201`, 1366.6 s, FINDINGS 32.1 -- through the SAME chain `extract_audio.py` uses (AC-3 5.1, ffmpeg's default downmix matrix, mono, 15,625 Hz), because a level measured through a different resampler is a level for a different encoder. Two statistics, and the difference between them is the whole argument: PEAK max |x| over the disc. What must fit under 511 for NOTHING to clamp. PASSAGE the loudest ~1 s window's peak and RMS. What the ear gets. A single sample 6 dB above everything else is a click and costs one clamp; a passage 6 dB above the gate window is where the recursion lives for fifteen thousand samples. It prints the attenuation each choice implies, in dB and as the shift `pack.py` would have to make, and it does NOT choose. Choosing needs the other half -- what attenuation costs at the quiet end, where the OKI step table's floor of 16 (12-bit units) does not scale with the signal -- and that is `--ladder`, which encodes real passages at real gains with `adpcm.CHIP` and reports the SNR. """ import argparse, getpass, json, os, subprocess, sys import numpy as np sys.path.insert(0, "tools/encoder") import adpcm BDROM = os.environ.get("DLX_BDROM") or f"/media/{getpass.getuser()}/BDROM" STREAM_DIR = f"{BDROM}/BDMV/STREAM" HZ = 15625 # the chip's rate, and the only one budgeted for FPS = 12 LUMP_FRAMES = 11 # FINDINGS 65.3's cadence: 11 frames of audio WIN = LUMP_FRAMES * HZ // FPS # 14,322 samples ~ 0.917 s -- one audio lump HOP = HZ // 4 # 0.25 s blocks; the window is 4 of them (rounded) CLAMP_LO, CLAMP_HI = adpcm.clamp_bounds(adpcm.CHIP["bits"]) # -512, 511 FULL12 = 2048 # what `s16 >> 4` maps full scale to def db(x, ref=FULL12): return -np.inf if x <= 0 else 20 * np.log10(x / ref) _PCM_CACHE = {} def pcm12(stream, start=None, dur=None): """One stream as 12-bit signed samples, through extract_audio.py's chain. Cached, because the scan, the census and the event walk are three passes over the same 20 million samples and the whole game is 40 MB of int16. """ ck = (stream, start, dur) if ck in _PCM_CACHE: return _PCM_CACHE[ck] cmd = ["ffmpeg", "-v", "error"] if start is not None: cmd += ["-ss", str(start)] if dur is not None: cmd += ["-t", str(dur)] cmd += ["-i", f"{STREAM_DIR}/{stream}.m2ts", "-vn", "-ac", "1", "-ar", str(HZ), "-f", "s16le", "-acodec", "pcm_s16le", "-"] p = subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE) if p.returncode: # 00176 is 3.0 s of mpeg2video with NO AUDIO TRACK AT ALL. That is a # fact about the disc, not a failure here, so it is reported rather # than swallowed -- but anything else is a real error. if b"does not contain any stream" not in p.stderr: raise SystemExit(f"ffmpeg failed on {stream}: " f"{p.stderr.decode(errors='replace')[:400]}") _PCM_CACHE[ck] = None return None x = np.frombuffer(p.stdout, "> 4 is -1. out = np.clip(x >> 4, -FULL12, FULL12 - 1).astype(np.int16) _PCM_CACHE[ck] = out return out def scan(streams): """Per-stream peak and loudest-passage statistics, in 12-bit units.""" rows, mute = [], [] for s in streams: x = pcm12(s) if x is None: mute.append(s) continue if x.size == 0: continue a = np.abs(x).astype(np.float64) nb = a.size // HOP if nb: bmax = a[:nb * HOP].reshape(nb, HOP).max(1) bsq = (a[:nb * HOP].reshape(nb, HOP) ** 2).sum(1) k = max(1, round(WIN / HOP)) if nb >= k: # sliding sum over k blocks == the ~1 s lump window cs = np.concatenate(([0.0], np.cumsum(bsq))) wrms = np.sqrt((cs[k:] - cs[:-k]) / (k * HOP)) wpk = np.array([bmax[i:i + k].max() for i in range(nb - k + 1)]) else: wrms = np.array([np.sqrt((a ** 2).mean())]) wpk = np.array([a.max()]) else: wrms = np.array([np.sqrt((a ** 2).mean())]) wpk = np.array([a.max()]) ipk = int(np.argmax(a)) irms = int(np.argmax(wrms)) rows.append(dict(stream=s, n=int(x.size), secs=x.size / HZ, peak=float(a.max()), peak_t=ipk / HZ, rms=float(np.sqrt((a ** 2).mean())), wpeak=float(wpk.max()), wrms=float(wrms.max()), wrms_t=irms * HOP / HZ)) return rows, mute def report(rows, mute): tot = sum(r["secs"] for r in rows) peak = max(rows, key=lambda r: r["peak"]) loud = max(rows, key=lambda r: r["wrms"]) disc_peak = peak["peak"] print(f"=== THE DISC, {len(rows)} streams, {tot:,.1f} s = {tot/60:.1f} min " f"(FINDINGS 32.1 says 1,366.6) ===\n") if mute: print(f" {len(mute)} stream(s) carry NO AUDIO TRACK: {', '.join(mute)}" f" -- a fact about the disc, and a case a shipping encoder has\n" f" to have an answer for (silence of the right length).\n") print(f'{"stream":>8}{"secs":>8}{"peak":>7}{"dBFS":>8}{"passage pk":>12}' f'{"passage rms":>13}{"dBFS":>8} at') for r in sorted(rows, key=lambda r: -r["wrms"])[:12]: print(f'{r["stream"]:>8}{r["secs"]:8.1f}{r["peak"]:7.0f}{db(r["peak"]):8.2f}' f'{r["wpeak"]:12.0f}{r["wrms"]:13.1f}{db(r["wrms"]):8.2f}' f' {r["wrms_t"]:6.2f} s') print(" (the twelve loudest PASSAGES; the table is sorted by the window " "RMS, not the peak)\n") print(f" DISC PEAK {disc_peak:.0f} of {FULL12} = {db(disc_peak):.2f} dBFS" f" ({peak['stream']} @ {peak['peak_t']:.2f} s)") print(f" LOUDEST PASSAGE peak {loud['wpeak']:.0f}, rms {loud['wrms']:.1f}" f" = {db(loud['wrms']):.2f} dBFS ({loud['stream']} @ {loud['wrms_t']:.2f} s)") print(f" THE CLAMP +{CLAMP_HI} / {CLAMP_LO} (it is not symmetric), " f"{db(CLAMP_HI):.2f} dBFS in the same units\n") need = disc_peak / CLAMP_HI print("=== WHAT THAT COSTS, AS A LEVEL ===\n") print(f" `s16 >> 4` is what pack.py does today and it puts the disc's own") print(f" peak at {disc_peak:.0f} against a clamp of {CLAMP_HI}: " f"{'OVER by' if need > 1 else 'under by'} {abs(20*np.log10(need)):.2f} dB.") print(f" Fitting the whole disc under the clamp with no sample clamped at") print(f" all needs a gain of {1/need:.4f} = {-20*np.log10(need):.2f} dB, i.e.") for sh in (4, 5, 6, 7): pk = disc_peak / (1 << (sh - 4)) mark = " <- fits" if pk <= CLAMP_HI else "" # ~ because a further right shift floors again and this halves; the # difference is one count and the column is a signpost, not a spec. print(f" s16 >> {sh} disc peak ~{pk:7.1f} " f"{'clamps' if pk > CLAMP_HI else 'clear':>6} by " f"{abs(20*np.log10(pk/CLAMP_HI)):5.2f} dB{mark}") print() # How much of the disc is actually above the clamp at today's level: the # number that decides whether this is a level question or a limiter question. return dict(rows=rows, disc_peak=disc_peak, peak_stream=peak["stream"], peak_t=peak["peak_t"], loud_stream=loud["stream"], loud_t=loud["wrms_t"], loud_rms=loud["wrms"], loud_wpeak=loud["wpeak"], clamp=CLAMP_HI, mute=mute) def clip_census(rows, streams, gains): """At each candidate gain, how many samples of the WHOLE DISC clamp? A peak is one number and this is the distribution behind it. A gain that clamps 12 samples in 22 minutes is a different object from one that clamps thousands, and the peak alone cannot tell them apart. """ print("=== THE CENSUS: how much of the disc is ABOVE the clamp, by gain ===\n") print(f'{"gain":>8}{"dB":>8}{"samples over":>14}{"of":>12}{"share":>10}' f'{"worst over":>12}') tot = 0 over = {g: 0 for g in gains} worst = {g: 0.0 for g in gains} for s in streams: x = pcm12(s) if x is None: continue a = np.abs(x).astype(np.float64) tot += a.size for g in gains: # ROUNDED, exactly as the ladder and pack.py requantise. Comparing # the float product instead makes 511/946 report one sample over # its own clamp, which is arithmetic about floats and not about # the disc. v = np.round(a * g) m = v > CLAMP_HI over[g] += int(m.sum()) if m.any(): worst[g] = max(worst[g], float(v.max() / CLAMP_HI)) for g in gains: w = f"{20*np.log10(worst[g]):.2f} dB" if worst[g] else "-" print(f'{g:8.4f}{20*np.log10(g):8.2f}{over[g]:14,}{tot:12,}' f'{100*over[g]/tot:9.4f}%{w:>12}') print() return dict(total=tot, over={f"{g:.4f}": over[g] for g in gains}) def clamp_events(streams, gain=1.0): """WHERE the over-clamp samples are, not just how many. 687 isolated samples in 22 minutes and one sustained 44 ms burst are the same census row and completely different sounds, and a clamp inside a recursion is not a clipped sample -- it is a wrong predictor state that the next nibble is applied to. So the run lengths are the statistic. """ runs = [] for st in streams: x = pcm12(st) if x is None: continue m = np.round(np.abs(x).astype(np.float64) * gain) > CLAMP_HI if not m.any(): continue d = np.diff(np.concatenate(([0], m.view(np.int8), [0]))) beg = np.where(d == 1)[0] end = np.where(d == -1)[0] for b, e in zip(beg, end): runs.append((int(e - b), st, b / HZ)) runs.sort(reverse=True) n = sum(r[0] for r in runs) print(f"=== WHERE THE CLAMPS ARE at gain {gain:.4f} " f"({len(runs)} events, {n:,} samples = {1000*n/HZ:.1f} ms) ===\n") print(f'{"run":>6}{"ms":>8} stream at') for r, st, t in runs[:10]: print(f'{r:6}{1000*r/HZ:8.2f} {st} {t:7.2f} s') if runs: print(f" longest run {runs[0][0]} samples = {1000*runs[0][0]/HZ:.2f} ms; " f"median run {sorted(r[0] for r in runs)[len(runs)//2]}") print() return dict(events=len(runs), samples=n, longest=runs[0][0] if runs else 0) def ladder(where, gains, dur, label): """Encode a real passage at each gain with adpcm.CHIP and report the SNR. This is the half a peak measurement cannot do. Attenuation buys headroom at the top and spends resolution at the bottom, because the OKI step table's floor is a constant 16 in 12-bit units and does not scale with the signal. The SNR is reported against the SCALED source, which is the honest comparison: the encoder's job is to reproduce what it was handed, and the listener's volume knob is not this project's problem. """ stream, start = where x = pcm12(stream, start, dur) print(f"=== THE LADDER: {label} -- {stream} @ {start:.2f} s, {dur:.2f} s, " f"{x.size:,} samples ===\n") print(f'{"gain":>8}{"dB":>8}{"src peak":>10}{"clamped":>9}{"SNR dB":>9}' f'{"vs 1.0":>8}') base = None out = [] for g in gains: src = np.clip(np.round(x * g), -FULL12, FULL12 - 1).astype(int).tolist() nib = adpcm.encode(src, variant=adpcm.CHIP["variant"], init=adpcm.CHIP["init"], bits=adpcm.CHIP["bits"]) rec = adpcm.decode(nib, variant=adpcm.CHIP["variant"], init=adpcm.CHIP["init"], bits=adpcm.CHIP["bits"]) s = np.array(src, dtype=np.float64) r = np.array(rec, dtype=np.float64) e = ((s - r) ** 2).sum() snr = 10 * np.log10((s ** 2).sum() / e) if e else np.inf nclamp = int((np.abs(s) > CLAMP_HI).sum()) if base is None: base = snr print(f'{g:8.4f}{20*np.log10(g):8.2f}{np.abs(s).max():10.0f}{nclamp:9,}' f'{snr:9.2f}{snr-base:+8.2f}') out.append(dict(gain=g, snr=snr, clamped=nclamp, peak=float(np.abs(s).max()))) print() return out def survey(rows, gains, n, dur, seed=20260825): """THE DISC, not three passages of it. Three hand-picked passages can be argued with; a sample cannot. `n` windows are drawn uniformly over the game's own timeline -- weighted by stream length, so a 24 s stream gets twenty times the draws of a 1.2 s one -- and every one is encoded at every gain with `adpcm.CHIP`. What is reported is the distribution: the mean SNR is what the level costs on average, and the WORST window is what it costs where it matters, because a level is chosen for the passage it fails on. """ rng = np.random.default_rng(seed) pool = [r for r in rows if r["secs"] >= dur] w = np.array([r["secs"] for r in pool], dtype=np.float64) w /= w.sum() picks = [] for _ in range(n): r = pool[int(rng.choice(len(pool), p=w))] t = float(rng.uniform(0, r["secs"] - dur)) picks.append((r["stream"], t)) print(f"=== THE SURVEY: {n} windows of {dur:.1f} s drawn over the whole " f"{sum(r['secs'] for r in rows)/60:.1f} min, encoded at every gain ===\n") src = [pcm12(st, t, dur) for st, t in picks] print(f'{"gain":>8}{"dB":>8}{"mean SNR":>10}{"median":>9}{"WORST":>8}' f'{"windows w/ clamp":>18}{"samples":>9}') out = [] for g in gains: snrs, nclamp, ncw = [], 0, 0 for x in src: v = np.clip(np.round(x * g), -FULL12, FULL12 - 1).astype(int) k = int((np.abs(v) > CLAMP_HI).sum()) nclamp += k ncw += 1 if k else 0 nib = adpcm.encode(v.tolist(), variant=adpcm.CHIP["variant"], init=adpcm.CHIP["init"], bits=adpcm.CHIP["bits"]) rec = np.array(adpcm.decode(nib, variant=adpcm.CHIP["variant"], init=adpcm.CHIP["init"], bits=adpcm.CHIP["bits"]), dtype=np.float64) f = v.astype(np.float64) e = ((f - rec) ** 2).sum() snrs.append(10 * np.log10((f ** 2).sum() / e) if e else np.inf) a = np.array(snrs) print(f'{g:8.4f}{20*np.log10(g):8.2f}{a.mean():10.2f}' f'{np.median(a):9.2f}{a.min():8.2f}{ncw:14} of {len(src)}{nclamp:9,}', flush=True) out.append(dict(gain=g, mean=float(a.mean()), median=float(np.median(a)), worst=float(a.min()), clamped=nclamp, windows=ncw)) print() return out def recover(stream, gain, control, K=64): """DOES A CLAMP OUTLIVE THE SAMPLE IT HAPPENS ON? 66.3 said it would. The worry was exact and it is the right worry for a recursive codec: a clamped accumulator is a WRONG STATE and the next nibble is applied to it, so the error should persist after the loud sample has gone. Measuring the error after a clamp run does show it elevated -- and that is not evidence, because the samples after a clamp run are LOUD samples, where the step is large and the error is large anyway. So the control is the same window at the gain that never clamps, rescaled to the same units and read at the SAME sample indices. What the ratio isolates is the clamp and nothing else. """ x = pcm12(stream).astype(np.float64) def enc(g): src = np.clip(np.round(x * g), -FULL12, FULL12 - 1).astype(int) nib = adpcm.encode(src.tolist(), variant=adpcm.CHIP["variant"], init=adpcm.CHIP["init"], bits=adpcm.CHIP["bits"]) rec = adpcm.decode(nib, variant=adpcm.CHIP["variant"], init=adpcm.CHIP["init"], bits=adpcm.CHIP["bits"]) return src.astype(np.float64), np.array(rec, dtype=np.float64) s1, r1 = enc(gain) s2, r2 = enc(control) e1 = np.abs(s1 - r1) e2 = np.abs(s2 - r2) / control * gain # the control, in gain's units m = np.abs(s1) > CLAMP_HI d = np.diff(np.concatenate(([0], m.view(np.int8), [0]))) ends = [e for e in np.where(d == -1)[0] if e + K <= e1.size] p1 = np.array([e1[e:e + K] for e in ends], dtype=np.float64).mean(0) p2 = np.array([e2[e:e + K] for e in ends], dtype=np.float64).mean(0) print(f"=== DOES THE CLAMP OUTLIVE THE SAMPLE? {stream}, gain {gain:g} " f"against a control at {control:g} ===\n") print(f" {int(m.sum())} samples clamp in {len(ends)} runs; the profile is " f"the mean |error| at each\n offset after a run ENDS, in 12-bit units, " f"against the same offsets of a\n window that never clamps at all.\n") print(f'{"after":>7}{"clamped":>10}{"control":>10}{"ratio":>8}') for i in (0, 1, 2, 4, 8, 16, 32, K - 1): print(f'{"+" + str(i):>7}{p1[i]:10.2f}{p2[i]:10.2f}{p1[i]/p2[i]:8.2f}') off = ~m print(f'\n off-clamp mean |err| {e1[off].mean():.2f} vs {e2[off].mean():.2f}') print(f' whole-window mean |err| {e1.mean():.2f} vs {e2.mean():.2f}') print(f' worst ratio over the {K} offsets: {(p1/p2).max():.2f}\n') return dict(stream=stream, gain=gain, control=control, runs=len(ends), clamped=int(m.sum()), worst_ratio=float((p1 / p2).max()), mean_err=float(e1.mean()), mean_err_control=float(e2.mean())) def main(): ap = argparse.ArgumentParser() ap.add_argument("--first", type=int, default=0) ap.add_argument("--last", type=int, default=201, help="the unique scene footage is 00000-00201 (FINDINGS 32.1); " "00215/00216/00223 are compilations of the same material") ap.add_argument("--ladder", action="store_true", help="also encode the loudest and a quiet passage at each gain") ap.add_argument("--dur", type=float, default=4.0, help="ladder passage seconds") ap.add_argument("--survey", type=int, default=0, help="encode N windows drawn over the whole game at each gain") ap.add_argument("--survey-dur", type=float, default=2.0) ap.add_argument("--recover", action="store_true", help="does a clamp outlive its sample? 66.3 said it would") ap.add_argument("--gate", action="store_true", help="assert FINDINGS 69's headline numbers, exit 1 if not") ap.add_argument("--json") a = ap.parse_args() streams = [f"{i:05d}" for i in range(a.first, a.last + 1)] streams = [s for s in streams if os.path.exists(f"{STREAM_DIR}/{s}.m2ts")] if not streams: sys.exit(f"no streams under {STREAM_DIR} -- is the Blu-ray mounted? " f"(DLX_BDROM)") rows, mute = scan(streams) summary = report(rows, mute) # The odd one is not a round number and is not meant to be: it is # CLAMP/disc peak, the gain at which the disc's own loudest sample lands # EXACTLY on the clamp, computed from the scan rather than typed in. exact = round(CLAMP_HI / summary["disc_peak"], 4) gains = sorted({1.0, 0.7071, exact, 0.5, 0.3536, 0.25}, reverse=True) summary["exact_gain"] = exact summary["census"] = clip_census(rows, streams, gains) summary["events"] = clamp_events(streams, 1.0) if a.ladder: # Three passages, because they answer three different questions. # PEAK what CLAMPING costs, since this is the only place on the disc # that clamps at today's level. # LOUD the loudest sustained window that is long enough to encode. # QUIET what ATTENUATION costs, which is the other end of the same # decision and the reason -15 dB is not free. long = [r for r in rows if r["secs"] >= a.dur] pk = max(rows, key=lambda r: r["peak"]) # the DISC peak, however short loud = max(long, key=lambda r: r["wrms"]) quiet = min(long, key=lambda r: r["wrms"]) at = lambda r, t: (r["stream"], min(max(0.0, t - a.dur / 2), max(0.0, r["secs"] - a.dur))) pkdur = min(a.dur, pk["secs"]) summary["ladder_peak"] = ladder( (pk["stream"], min(max(0.0, pk["peak_t"] - pkdur / 2), max(0.0, pk["secs"] - pkdur))), gains, pkdur, "THE DISC PEAK ITSELF") summary["ladder_loud"] = ladder(at(loud, loud["wrms_t"]), gains, a.dur, "THE LOUDEST SUSTAINED PASSAGE") summary["ladder_quiet"] = ladder(at(quiet, quiet["wrms_t"]), gains, a.dur, "A QUIET PASSAGE, for the other end") if a.survey: summary["survey"] = survey(rows, gains, a.survey, a.survey_dur) if a.recover: summary["recover"] = recover(summary["peak_stream"], 1.0, exact) if a.gate: expect = dict(disc_peak=946.0, peak_stream="00200", clamp=511, events=402, over=687) bad = [] for k, v in expect.items(): got = (summary["events"]["events"] if k == "events" else summary["events"]["samples"] if k == "over" else summary[k]) if got != v: bad.append(f"{k}: expected {v}, measured {got}") if bad: print("LEVEL GATE RED -- the disc does not measure as FINDINGS 69 " "recorded it:") for b in bad: print(" " + b) print(" (a different pressing is a legitimate cause; a different " "ffmpeg downmix is not)") sys.exit(1) print("LEVEL GATE GREEN: disc peak 946 of 2048 at 00200, 5.35 dB over " "the chip's 511,\n 687 samples in 402 events = 44.0 ms of the " "game's 21.5 min of audio.") if a.json: json.dump(summary, open(a.json, "w"), indent=1) print(f"-> {a.json}") if __name__ == "__main__": main()