True peak after encoding: 17 chart tracks

"""Measure a week's chart through the codecs the platforms actually serve.

Seventeen tracks bought as lossless downloads (Qobuz, Japanese store, 2026-08-10),
one per artist, taken in chart order from Spotify Global for the week of 6 August.
Nothing here is chosen by genre, year or sound: the selection rule is written down
in docs/MEASUREMENT_STUDY.md and was fixed before anything was measured.

What this answers that the existing chart measurements do not: those measure the
file the label delivered. This measures what comes back out of the encoder the
listener is actually served, which is a different file.

Run it in the container: Homebrew's ffmpeg ships without libvorbis, and the
container's is the build the product actually serves.

    docker compose run --rm --no-deps \\
      -v "/Volumes/Extreme SSD/003_music/Download:/audio:ro" \\
      -v "$PWD/docs:/app/docs" \\
      web uv run --directory /app python /app/scripts/chart_codec_study.py \\
      --root /audio

Mounting an external drive that way is slow to start: a couple of minutes
can pass before the first line appears, on a cold disk. It is not stuck. On 2026-08-10
that wait got read as a failure, the audio was copied under media/ to work
around it, and the copy was never needed.

Writes docs/data/chart_codec_study.json. Only numbers are written; the audio
stays where it is and no encode is kept. Measuring a recording and reporting
the numbers is one thing, redistributing it is another.

The null-test primitives are imported from scripts/codec_null_test.py, which is
the script published with the guide, so these figures come out of the code we
already showed people rather than a second private copy of it.
"""

import argparse
import json
import os
import shutil
import subprocess
import sys
import tempfile
import time

import numpy as np
import soundfile as sf

HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)
sys.path.insert(0, os.path.join(HERE, '..', 'src'))

from codec_null_test import (                                    # noqa: E402
    align, band_residual, check_alignment, db, rms, specs, to_wav,
    true_peak_ebur128,
)
from study_metrics import (                                     # noqa: E402
    analyze_dynamic_range, analyze_loudness, generate_track_peaks,
)

DEFAULT_ROOT = '.'
OUT_DIR = os.path.join(HERE, '..', 'docs', 'data')

# Titles are the seller's, verbatim. Where a shop's metadata is odd - Thriller's
# sixth track is sold as "(Single Version)" while its neighbours say "(Album
# Version)" - the oddity is carried through rather than tidied, because the table
# has to describe what was bought.
MANIFEST = [
    # key            artist                    title as sold                         path
    ('ariana',      'Ariana Grande',          'hate that i made you love me',
     'アリアナ・グランデ/hate that i made you love me/01 - アリアナ・グランデ - hate that i made you love me.wav'),
    ('shakira',     'Shakira, Burna Boy',     'Dai Dai',
     'シャキーラ, バーナ・ボーイ/Dai Dai/01 - シャキーラ, バーナ・ボーイ - Dai Dai.wav'),
    ('katseye',     'KATSEYE',                'Animal',
     'KATSEYE/Animal/01 - KATSEYE - Animal.wav'),
    ('tempercity',  'Temper City',            'Self Aware',
     'Temper City/Self Aware/01 - Temper City - Self Aware.wav'),
    ('sombr',       'sombr',                  'back to friends',
     'Sombr/back to friends/01 - Sombr - back to friends.wav'),
    ('langley',     'Ella Langley',           "Choosin' Texas",
     "ELLA LANGLEY/Choosin' Texas/01 - ELLA LANGLEY - Choosin' Texas.wav"),
    ('djo',         'Djo',                    'End of Beginning',
     'Djo/End of Beginning/01 - Djo - End of Beginning.wav'),
    ('oliviadean',  'Olivia Dean',            'Man I Need',
     'オリヴィア・ディーン/Man I Need/01 - オリヴィア・ディーン - マン・アイ・ニード.wav'),
    ('bts',         'BTS',                    'SWIM',
     'BTS (防弾少年団)/ARIRANG/07 - BTS (防弾少年団) - SWIM.wav'),
    ('bieber',      'Justin Bieber',          'Beauty And A Beat (Album Version)',
     'Justin Bieber/Believe/10 - Justin Bieber - ビューティー・アンド・ア・ビート (Album Version).wav'),
    ('malcolmtodd', 'Malcolm Todd',           'Earrings',
     'Malcolm Todd/Sweet Boy/01 - Malcolm Todd - Earrings.wav'),
    ('mj',          'Michael Jackson',        'Billie Jean (Single Version)',
     'Michael Jackson/Thriller/06 - Michael Jackson - ビリー・ジーン (Single Version).wav'),
    ('fike',        'Dominic Fike',           'Babydoll',
     "ドミニク・ファイク/Don't Forget About Me, Demos/03 - ドミニク・ファイク - Babydoll.wav"),
    ('oasis',       'Oasis',                  'Wonderwall (Remastered)',
     'オアシス/(What’s The Story) Morning Glory (20th Anniversary Standard Edition)/03 - オアシス - Wonderwall (Remastered).wav'),
    ('katyperry',   'Katy Perry',             'The One That Got Away',
     'ケイティ・ペリー/Teenage Dream/07 - ケイティ・ペリー - ワン・ザット・ゴット・アウェイ.wav'),
    ('rodrigo',     'Olivia Rodrigo',         'stupid song',
     'オリヴィア・ロドリゴ/you seem pretty sad for a girl so in love/02 - オリヴィア・ロドリゴ - stupid song.wav'),
    ('iris',        'Goo Goo Dolls',          'Iris',
     'Goo Goo Dolls/Dizzy up the Girl/11 - Goo Goo Dolls - Iris.wav'),
]

# Two tracks were bought twice, in both tiers the shop sold. Iris is the question
# (a 1998 catalog title: is the high-res a different master?) and Djo is the
# control (a recent release, where the two tiers should be the same master). One
# pair on its own cannot tell a real difference from the error in the comparison,
# which is what the second pair is for.
PAIRS = [
    ('iris', 'Goo Goo Dolls', 'Iris (CD tier)',
     'Goo Goo Dolls/Dizzy up the Girl/01-11 Iris CD.flac'),
    ('djo', 'Djo', 'End of Beginning (CD tier)',
     'Djo/End of Beginning/01-01 End of Beginning.flac'),
]

BAND_WINDOW_S = 60      # centered excerpt for the per-band FFT, see band_window()

# Spotify publishes three Vorbis bitrates and most listeners are not on the
# lowest one. Measuring only 96 kbps would be measuring the worst case and
# calling it the platform, which is the first thing anyone would say about it,
# and they would be right. 96 stays in as the floor.
EXTRA_SPECS = {
    'spotify_160': {'encoder': 'libvorbis', 'bitrate': '160k', 'ext': 'ogg',
                    'label': 'Ogg Vorbis 160 kbps'},
    'spotify_320': {'encoder': 'libvorbis', 'bitrate': '320k', 'ext': 'ogg',
                    'label': 'Ogg Vorbis 320 kbps'},
}


def all_specs():
    """The product's three, plus the two extra Vorbis bitrates."""
    out = dict(specs())
    out.update(EXTRA_SPECS)
    return out


def toolchain():
    """Exactly what did the encoding, so the run can be repeated or disputed.

    None of these are the encoders the platforms run. Apple uses its own AAC
    implementation and Spotify builds its own Vorbis; a page that reports what
    an encoder did has to say which encoder.
    """
    import platform
    out = {'platform': platform.platform()}
    try:
        r = subprocess.run(['ffmpeg', '-version'], capture_output=True, text=True)
        lines = (r.stdout or '').splitlines()
        out['ffmpeg'] = lines[0] if lines else ''
        for ln in lines:
            if ln.startswith('configuration:'):
                out['ffmpeg_configuration'] = ln[len('configuration:'):].strip()
    except OSError:
        out['ffmpeg'] = 'unavailable'
    for lib in ('libvorbis', 'libopus', 'aac', 'wavpack'):
        r = subprocess.run(['ffmpeg', '-hide_banner', '-h', 'encoder=' + lib],
                           capture_output=True, text=True)
        head = (r.stdout or '').splitlines()
        out.setdefault('encoders', {})[lib] = head[0].strip() if head else 'missing'
    return out


def platforms():
    """Read the platform table out of the shipped tool rather than retyping it.

    Four copies of these numbers exist in the product and a test keeps them in
    step; a fifth copy living in a study script is exactly how an article ends up
    quoting a figure the product stopped using.
    """
    import re
    path = os.path.join(HERE, '..', 'src', 'theme', 'static', 'js',
                        'sn-loudness-tool.js')
    src = open(path, encoding='utf-8').read()
    m = re.search(r'var PLATFORMS = \[(.*?)\n    \];', src, re.DOTALL)
    if not m:
        sys.exit('could not find PLATFORMS in sn-loudness-tool.js')
    out = []
    for line in m.group(1).splitlines():
        line = line.strip()
        if not line.startswith('{'):
            continue
        key = re.search(r"key:\s*'([^']+)'", line)
        name = re.search(r"name:\s*'([^']+)'", line)
        lufs = re.search(r'lufs:\s*(-?[\d.]+|null)', line)
        tp = re.search(r'tp:\s*(-?[\d.]+)', line)
        if not (key and name and lufs):
            continue
        out.append({
            'key': key.group(1),
            'name': name.group(1),
            'lufs': None if lufs.group(1) == 'null' else float(lufs.group(1)),
            'tp': float(tp.group(1)) if tp else None,
            'down_only': 'down: true' in line,
            'normalizes': 'norm: false' not in line,
        })
    return out


def platform_gain(master_lufs, p):
    """The gain a platform applies. Arithmetic, not a model: target minus master,
    floored at zero where the platform only turns things down."""
    if not p['normalizes'] or p['lufs'] is None or master_lufs is None:
        return 0.0
    delta = p['lufs'] - master_lufs
    return min(0.0, delta) if p['down_only'] else delta


def band_window(x, sr):
    """A centered excerpt for the per-band FFT.

    The residual itself is computed on the whole track. The band breakdown is
    not, because a 192 kHz five-minute file is a 55-million-point transform and
    the same excerpt length on every track is more comparable than a full
    transform whose resolution changes with the length of the song.
    """
    n = int(BAND_WINDOW_S * sr)
    if len(x) <= n:
        return 0, len(x)
    start = (len(x) - n) // 2
    return start, start + n


def our_true_peak(path):
    """Our 4x FFT-upsampled true peak, plus the sample peak it sat above.

    Both this and ffmpeg's ebur128 are recorded for every file. Across the
    seventeen they agree within 0.11 dB at every rate from 44.1 to 192 kHz,
    which is the answer to anyone who would rather not trust a filter we wrote.

    An earlier version of this comment claimed ebur128 stops oversampling at
    high rates. It does not. That was inferred from one 192 kHz file where the
    true peak sits on the sample peak - which is what a densely sampled file
    looks like, not what a broken meter looks like.
    """
    try:
        _peaks, sample_peak, tp, _tp_data = generate_track_peaks(path)
        return tp, sample_peak
    except Exception as e:                                   # pragma: no cover
        print('    true peak failed: %s' % e)
        return None, None


def source_metrics(path):
    loud = analyze_loudness(path) or {}
    tp, sample_peak = our_true_peak(path)
    st_max = loud.get('short_term_max')
    info = sf.info(path)
    return {
        'sample_rate': info.samplerate,
        'channels': info.channels,
        'subtype': info.subtype,
        'duration_s': round(info.frames / info.samplerate, 3),
        'lufs': loud.get('lufs'),
        'lra': loud.get('lra'),
        'short_term_max': st_max,
        'true_peak_4x': tp,
        'sample_peak': sample_peak,
        'true_peak_ebur128': loud.get('true_peak'),
        'dynamic_range': analyze_dynamic_range(path),
        'psr': (tp - st_max) if (tp is not None and st_max is not None) else None,
    }


def delivery_rate(spec):
    """The rate the encoder is actually fed.

    Nobody hands Spotify a 192 kHz file and gets 192 kHz Ogg back: the platform
    resamples first, and the streams are 44.1 kHz for Vorbis and AAC and 48 kHz
    for Opus. Feeding the encoder the file's own rate would measure something no
    listener receives - and libvorbis refuses 96 kHz at these bitrates anyway,
    which is how this got noticed.

    Nulling happens at this rate too, against a reference resampled to it, so
    the residual is the codec and not the sample-rate conversion in front of it.
    For a 44.1 kHz source that makes this identical to the published script.
    """
    return 48000 if spec['encoder'] == 'libopus' else 44100


def encode_spec(src, dst, spec):
    """As the published script's encode(), plus an explicit sample format.

    The references here are 32-bit float so that resampling a master sitting at
    full scale cannot clip on the way in. FLAC has no float mode, so ffmpeg
    quietly picks 16-bit for the lossless control and the control then reports
    -97 dB - which is 16-bit quantisation noise wearing the costume of a broken
    method. The three real codecs take float happily and are unaffected.
    """
    cmd = ['ffmpeg', '-y', '-hide_banner', '-loglevel', 'error', '-i', src,
           '-map', '0:a:0', '-c:a', spec['encoder']]
    spec['_argv'] = cmd  # filled in below, recorded so the run can be repeated
    if spec.get('bitrate') and spec['bitrate'] != '0':
        cmd += ['-b:a', spec['bitrate']]
    if spec['encoder'] == 'aac':
        cmd += ['-profile:a', 'aac_low']
    if spec.get('sample_fmt'):
        cmd += ['-sample_fmt', spec['sample_fmt']]
    full = cmd + [dst]
    spec['_argv'] = ' '.join(full[:-1] + ['out.' + spec['ext']])
    run_ff(full)


def run_ff(cmd):
    r = subprocess.run(cmd, capture_output=True, text=True)
    if r.returncode != 0:
        raise RuntimeError((r.stderr or 'ffmpeg failed')[:400])


def codec_pass(ref_path, spec, tmp, sr):
    """Encode, decode, align, subtract. Same steps as the published script.

    ``ref_path`` is the master already resampled to the encoder's delivery rate
    and ``sr`` is that rate.
    """
    enc = os.path.join(tmp, 'e.' + spec['ext'])
    dec = os.path.join(tmp, 'd.wav')
    for p in (enc, dec):
        if os.path.exists(p):
            os.remove(p)
    encode_spec(ref_path, enc, spec)
    to_wav(enc, dec, sr)

    ref, _ = sf.read(ref_path, always_2d=True, dtype='float64')
    ref = ref[:, 0]
    test, _ = sf.read(dec, always_2d=True, dtype='float64')
    test = test[:, 0]

    lag = align(ref, test)
    if lag >= 0:
        a, b = ref[:len(ref) - lag], test[lag:lag + len(ref) - lag]
    else:
        a, b = ref[-lag:], test[:len(ref) + lag]
    n = min(len(a), len(b))
    a, b = a[:n], b[:n]

    # A codec can hand back a fraction of a dB off level, and uncorrected that
    # reads as damage. Both figures are reported so neither has to be trusted.
    g = float(np.dot(a, b) / np.dot(b, b)) if np.dot(b, b) else 1.0
    residual_raw = db(rms(a - b) / rms(a))
    diff_gain = a - b * g
    residual = db(rms(diff_gain) / rms(a))

    lo, hi = band_window(a, sr)
    bands = band_residual(a[lo:hi], diff_gain[lo:hi], sr)

    enc_tp_4x, enc_sample_peak = our_true_peak(dec)
    out = {
        'label': spec['label'],
        'lag_samples': int(lag),
        'gain_correction_db': db(abs(g)),
        'residual_db': None if residual == -np.inf else round(residual, 2),
        'residual_uncorrected_db': (None if residual_raw == -np.inf
                                    else round(residual_raw, 2)),
        'bands': [{'lo': lo_, 'hi': hi_,
                   'rel_db': None if r == -np.inf else round(r, 1)}
                  for (lo_, hi_), r in bands],
        'true_peak_4x': enc_tp_4x,
        'sample_peak': enc_sample_peak,
        'true_peak_ebur128': true_peak_ebur128(dec),
        'encoded_bytes': os.path.getsize(enc),
        'command': spec.get('_argv'),
    }
    del ref, test, a, b, diff_gain
    return out


def study_track(key, artist, title, path, tmp, quick=False):
    print('  %s' % os.path.basename(path))
    t0 = time.time()
    row = {'key': key, 'artist': artist, 'title': title,
           'file': os.path.basename(path)}
    row['source'] = source_metrics(path)

    # One resampled reference per delivery rate, reused by every encoder that
    # feeds off it. Measured as well as encoded: the platform's own downsample
    # moves the peaks before the codec has touched anything, and separating
    # those two steps is most of what makes the codec figure mean something.
    refs = {}
    for rate in sorted({delivery_rate(s) for s in all_specs().values()} | {44100}):
        rp = os.path.join(tmp, 'ref%d.wav' % rate)
        to_wav(path, rp, rate)
        tp, sp = our_true_peak(rp)
        refs[rate] = rp
        row.setdefault('resampled', {})[str(rate)] = {
            'true_peak_4x': tp, 'sample_peak': sp,
        }

    # The control runs first and prints first. A pipeline that cannot return
    # "identical" when the answer is identical is not measuring anything, so if
    # this is not far below the codec rows the rest of the row is decoration.
    # WavPack rather than FLAC, because FLAC has no float mode and several of
    # these masters land above full scale once they are resampled to 44.1 kHz.
    # Writing those to a fixed-point file clips them, and the control then
    # reports the clipping as if the method were broken. WavPack round-trips
    # float bit-identically, so the control measures only what it is for.
    row['control'] = codec_pass(refs[44100],
                                {'encoder': 'wavpack', 'bitrate': '0',
                                 'ext': 'wv',
                                 'label': 'WavPack (lossless float)'},
                                tmp, 44100)
    print('    control %s' % (row['control']['residual_db'] if
                              row['control']['residual_db'] is not None
                              else 'perfect null'))

    row['codecs'] = {}
    if not quick:
        for pkey, spec in all_specs().items():
            rate = delivery_rate(spec)
            c = codec_pass(refs[rate], spec, tmp, rate)
            c['encoder_input_rate'] = rate
            row['codecs'][pkey] = c
            print('    %-18s residual %6s dB   TP %+.2f -> %+.2f (4x, %d Hz)'
                  % (spec['label'], c['residual_db'],
                     row['resampled'][str(rate)]['true_peak_4x'] or 0.0,
                     c['true_peak_4x'] or 0.0, rate))

    row['platform_gain_db'] = {
        p['key']: round(platform_gain(row['source']['lufs'], p), 2)
        for p in platforms()
    }
    row['seconds'] = round(time.time() - t0, 1)
    return row


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument('--root', default=DEFAULT_ROOT)
    ap.add_argument('--only', help='one manifest key, for iterating')
    ap.add_argument('--quick', action='store_true',
                    help='source metrics and control only, no codec passes')
    # Only ./src, ./media and ./scripts are mounted into the running web
    # container, so `docker compose exec` needs somewhere else for the results
    # to land. `docker compose run -v` (see the module docstring) can mount
    # docs/ directly and this can be left alone.
    ap.add_argument('--out', default=OUT_DIR)
    args = ap.parse_args()
    out_dir = args.out

    if not shutil.which('ffmpeg'):
        sys.exit('ffmpeg not found')

    # Every row in this study reports "aligned +0 samples", because ffmpeg
    # compensates encoder delay and Opus pre-skip when it decodes. That is also
    # precisely what a dead aligner looks like, so it has to prove itself on
    # known shifts before any of those zeros are believed.
    check_alignment()

    entries = [e for e in MANIFEST if not args.only or e[0] == args.only]
    pairs = [e for e in PAIRS if not args.only or e[0] == args.only]

    missing = [p for _, _, _, p in entries + pairs
               if not os.path.exists(os.path.join(args.root, p))]
    if missing:
        sys.exit('missing:\n  ' + '\n  '.join(missing))

    tmp = tempfile.mkdtemp(prefix='chartstudy-')
    result = {
        'measured_at': time.strftime('%Y-%m-%d'),
        'chart': 'Spotify Global, week of 2026-08-06',
        'selection': 'one track per artist, in chart order',
        'source': 'purchased lossless downloads (Qobuz, Japanese store)',
        'band_window_s': BAND_WINDOW_S,
        'toolchain': toolchain(),
        'tracks': [],
        'pairs': [],
    }

    print('tracks')
    for key, artist, title, rel in entries:
        result['tracks'].append(
            study_track(key, artist, title, os.path.join(args.root, rel),
                        tmp, args.quick))

    # The second tier of the two tracks sold in two. Source metrics only: the
    # question is whether the master differs, which loudness and range answer,
    # and running the codecs on both tiers would double the table for nothing.
    print('second tiers')
    for key, artist, title, rel in pairs:
        p = os.path.join(args.root, rel)
        print('  %s' % os.path.basename(p))
        result['pairs'].append({'key': key, 'artist': artist, 'title': title,
                                'file': os.path.basename(p),
                                'source': source_metrics(p)})

    os.makedirs(out_dir, exist_ok=True)
    out_json = os.path.join(out_dir, 'chart_codec_study.json')
    with open(out_json, 'w', encoding='utf-8') as f:
        json.dump(result, f, indent=2, ensure_ascii=False)
    print('\nwrote %s' % out_json)

    bad = [t['key'] for t in result['tracks']
           if t['control']['residual_db'] is not None
           and t['control']['residual_db'] > -120]
    if bad:
        print('WARNING: the lossless control did not null on: %s' % ', '.join(bad))
        print('Nothing else in this file means anything until it does.')


if __name__ == '__main__':
    main()

"""Add an Apple-encoder AAC column to the chart study, measured on macOS.

The main run happens in the container, where ffmpeg has libvorbis. Its AAC
column is ffmpeg's own encoder, which is not what Apple Music is served with,
and saying so in the copy is honest but weaker than measuring the real one.
macOS ships Apple's AAC encoder as afconvert, so this pass runs on the host and
merges a second AAC column into the same JSON.

    scripts/chart_codec_apple.py [--root DIR]

Requires afconvert (macOS), ffmpeg for decoding, and numpy/soundfile. It reads
docs/data/chart_codec_study.json, adds codecs['apple_native'] to every track and
a second toolchain entry, and writes the file back. Run it after the main study,
never before: it will not create rows it does not find.
"""

import argparse
import json
import os
import platform
import shutil
import subprocess
import sys
import tempfile

import numpy as np
import soundfile as sf

HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)
sys.path.insert(0, os.path.join(HERE, '..', 'src'))

from codec_null_test import align, band_residual, db, rms, to_wav   # noqa: E402
from chart_codec_study import (                                     # noqa: E402
    BAND_WINDOW_S, DEFAULT_ROOT, MANIFEST, band_window, our_true_peak,
    true_peak_ebur128,
)

DATA = os.path.join(HERE, '..', 'docs', 'data', 'chart_codec_study.json')

# Constrained VBR at 256 kbps with the codec's quality dial at maximum. Apple
# Music is served as 256 kbps VBR AAC; this is that bitrate through Apple's
# encoder rather than a guess at Apple's exact internal settings, and the
# command line goes into the data so the choice is visible.
AFCONVERT = ['-f', 'm4af', '-d', 'aac', '-b', '256000', '-q', '127', '-s', '2']


def run(cmd):
    r = subprocess.run(cmd, capture_output=True, text=True)
    if r.returncode != 0:
        raise RuntimeError(' '.join(cmd[:4]) + ': ' + (r.stderr or '')[:300])


def apple_version():
    r = subprocess.run(['afconvert', '--help'], capture_output=True, text=True)
    first = [ln.strip() for ln in (r.stdout or '').splitlines() if ln.strip()]
    return {
        'tool': 'afconvert (Apple CoreAudio)',
        'version': next((ln for ln in first if ln.startswith('Version')), ''),
        'host': platform.platform(),
    }


def measure(src_path, tmp):
    ref_path = os.path.join(tmp, 'ref.wav')
    enc = os.path.join(tmp, 'a.m4a')
    dec = os.path.join(tmp, 'a.wav')
    for p in (ref_path, enc, dec):
        if os.path.exists(p):
            os.remove(p)

    to_wav(src_path, ref_path, 44100)
    cmd = ['afconvert'] + AFCONVERT + [ref_path, enc]
    run(cmd)
    to_wav(enc, dec, 44100)

    ref, _ = sf.read(ref_path, always_2d=True, dtype='float64')
    ref = ref[:, 0]
    test, _ = sf.read(dec, always_2d=True, dtype='float64')
    test = test[:, 0]

    lag = align(ref, test)
    if lag >= 0:
        a, b = ref[:len(ref) - lag], test[lag:lag + len(ref) - lag]
    else:
        a, b = ref[-lag:], test[:len(ref) + lag]
    n = min(len(a), len(b))
    a, b = a[:n], b[:n]

    g = float(np.dot(a, b) / np.dot(b, b)) if np.dot(b, b) else 1.0
    diff = a - b * g
    residual = db(rms(diff) / rms(a))
    lo, hi = band_window(a, 44100)
    tp, sample_peak = our_true_peak(dec)
    return {
        'label': 'AAC-LC 256 kbps (Apple encoder)',
        'lag_samples': int(lag),
        'gain_correction_db': db(abs(g)),
        'residual_db': None if residual == -np.inf else round(residual, 2),
        'bands': [{'lo': l_, 'hi': h_,
                   'rel_db': None if r == -np.inf else round(r, 1)}
                  for (l_, h_), r in band_residual(a[lo:hi], diff[lo:hi], 44100)],
        'true_peak_4x': tp,
        'sample_peak': sample_peak,
        'true_peak_ebur128': true_peak_ebur128(dec),
        'encoded_bytes': os.path.getsize(enc),
        'encoder_input_rate': 44100,
        'command': ' '.join(['afconvert'] + AFCONVERT + ['in.wav', 'out.m4a']),
    }


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument('--root', default=DEFAULT_ROOT)
    args = ap.parse_args()
    for tool in ('afconvert', 'ffmpeg'):
        if not shutil.which(tool):
            sys.exit('%s not found; this pass only runs on macOS' % tool)

    d = json.load(open(DATA, encoding='utf-8'))
    by_key = {t['key']: t for t in d['tracks']}
    paths = {k: os.path.join(args.root, rel) for k, _a, _t, rel in MANIFEST}

    tmp = tempfile.mkdtemp(prefix='apple-aac-')
    for key, track in by_key.items():
        src = paths.get(key)
        if not src or not os.path.exists(src):
            sys.exit('missing audio for %s' % key)
        r = measure(src, tmp)
        track['codecs']['apple_native'] = r
        ffm = track['codecs']['apple']
        print('%-22s ffmpeg %6.1f dB / TP %+.2f    Apple %6.1f dB / TP %+.2f'
              % (track['artist'], ffm['residual_db'], ffm['true_peak_4x'],
                 r['residual_db'], r['true_peak_4x']))

    d.setdefault('toolchain', {})['apple'] = apple_version()
    with open(DATA, 'w', encoding='utf-8') as f:
        json.dump(d, f, indent=2, ensure_ascii=False)
    print('\nmerged apple_native into %s' % DATA)


if __name__ == '__main__':
    main()

"""Standalone measurement functions for the published chart-codec study.

Snapshot from projects/audio.py, 2026-09-26. The application upload-size guard
is omitted so these functions require no Django settings. Calculations are
unchanged. Keep this file beside the other three downloaded scripts.
Requires ffmpeg, numpy and soundfile. All input audio remains local.
"""
import os
import logging
from contextlib import contextmanager


def _resample_loudness_series(pairs, n=2048, max_gap_s=1.0):
    """Bucket (time, short-term-LUFS) pairs into ``n`` evenly-timed buckets
    (averaging within each), aligned to the waveform's time axis. Buckets with
    no data; silent stretches, the initial 3 s gating ramp (values < -70 are
    dropped); are ``None`` so the front-end breaks the line there rather than
    diving to -120. Returns a list of ``n`` floats/None, or None if no usable
    data.

    ffmpeg emits a short-term value only ~every 0.1 s, so a short track has far
    fewer source points than ``n`` buckets; naive bucketing then leaves ~half
    the buckets empty, interspersed, and the front-end breaks the line at every
    hole → the curve renders as a dotted line. So we linearly interpolate across
    SHORT interior gaps (resolution artefacts, < ``max_gap_s``) and keep only
    long gaps (real silence) as ``None`` for the line to break there."""
    if not pairs:
        return None
    dur = pairs[-1][0]
    if not dur or dur <= 0:
        return None
    sums = [0.0] * n
    counts = [0] * n
    for t, s in pairs:
        if s <= -70:           # gating ramp / silence placeholder; skip
            continue
        bi = int(t / dur * n)
        if bi < 0:
            bi = 0
        elif bi >= n:
            bi = n - 1
        sums[bi] += s
        counts[bi] += 1
    if not any(counts):
        return None
    out = [round(sums[i] / counts[i], 1) if counts[i] else None for i in range(n)]

    # Bridge short interior gaps so the contour is continuous through dense
    # audio; leave long gaps (real silence) broken.
    bucket_s = dur / n
    max_run = max(1, int(round(max_gap_s / bucket_s))) if bucket_s > 0 else n
    i = 0
    while i < n:
        if out[i] is not None:
            i += 1
            continue
        j = i
        while j < n and out[j] is None:
            j += 1
        # out[i..j-1] is a None run bounded by out[i-1] and out[j]; bridge it
        # only when it sits between two real values and is short enough.
        if i > 0 and j < n and (j - i) <= max_run:
            lo, hi = out[i - 1], out[j]
            span = j - (i - 1)
            for k in range(i, j):
                f = (k - (i - 1)) / span
                out[k] = round(lo * (1 - f) + hi * f, 1)
        i = j
    return out


def analyze_loudness(input_path, timeout=180):
    """EBU R128 loudness + short-term max via ffmpeg ``ebur128`` filter.

    Returns ``{'lufs', 'true_peak', 'lra', 'short_term_max'}`` or None.

    We deliberately do NOT chain ``loudnorm`` here; when its dynamic-norm
    output is piped into ebur128 in the same graph, the short-term LUFS
    stream reflects the *normalized* signal, which understates the master's
    real short-term loudness by ~15 dB. ebur128 alone gives us the original
    integrated I, range LRA, true peak (oversampled) AND the per-time S
    stream we scan for short_term_max.
    """
    import subprocess
    import shutil
    import math
    import re
    log = logging.getLogger('soneam.loudness')

    if not shutil.which('ffmpeg'):
        return None
    try:
        result = subprocess.run([
            'ffmpeg', '-hide_banner', '-nostats',
            '-i', input_path,
            '-af', 'ebur128=peak=true',
            '-f', 'null', '-',
        ], capture_output=True, text=True, timeout=timeout)
        err = result.stderr or ''

        # ── Summary block at end of stderr (multi-line, use DOTALL) ────────
        def _grab(pattern):
            m = re.search(pattern, err, re.DOTALL)
            if not m:
                return None
            try:
                v = float(m.group(1))
            except (TypeError, ValueError):
                return None
            return v if math.isfinite(v) else None

        lufs      = _grab(r'Integrated loudness:.*?I:\s*(-?[\d.]+)\s*LUFS')
        lra       = _grab(r'Loudness range:.*?LRA:\s*(-?[\d.]+)\s*LU')
        # ebur128 prints "True peak: Peak: X dBFS"; that *is* the oversampled
        # true peak; we keep it stored under true_peak_dbtp for consistency
        # with the rest of the app even though ffmpeg labels it dBFS.
        true_peak = _grab(r'True peak:.*?Peak:\s*(-?[\d.]+)\s*dBFS')

        # ── Per-time short-term LUFS (the loudness-over-time curve) ─────────
        # ebur128 logs a line per ~0.1s carrying both t: (seconds) and S:
        # (3-second short-term LUFS). Collect (t, S), then bucket into a fixed
        # series aligned to the waveform so the front-end can draw it as a
        # polyline over the wave (like the TP envelope).
        st_pairs = []
        for m in re.finditer(r't:\s*([\d.]+)[^\n]*?S:\s*(-?[\d.]+)', err):
            try:
                t = float(m.group(1)); s = float(m.group(2))
            except (TypeError, ValueError):
                continue
            if math.isfinite(t) and math.isfinite(s):
                st_pairs.append((t, s))

        # Loudest moment (ignoring the gating ramp / silence padding near -120).
        st_vals = [s for (_, s) in st_pairs if s > -70]
        st_max = max(st_vals) if st_vals else None
        st_series = _resample_loudness_series(st_pairs, n=2048)

        return {
            'lufs':           lufs,
            'true_peak':      true_peak,
            'lra':            lra,
            'short_term_max': st_max,
            'st_series':      st_series,
        }
    except subprocess.TimeoutExpired:
        log.warning('ebur128 timed out for %s', os.path.basename(input_path))
        return None
    except Exception as e:
        log.warning('ebur128 failed for %s: %s', os.path.basename(input_path), e)
        return None


@contextmanager
def _analysis_soundfile(path):
    """Read WAVs whose ancillary chunks libsndfile rejects, without touching
    the original or converting its PCM. ffmpeg rewrites only the container.
    The temporary copy is used solely for metrics and removed on every exit.
    """
    import soundfile as sf
    import subprocess
    import tempfile
    try:
        audio = sf.SoundFile(path)
    except RuntimeError:
        with open(path, 'rb') as source:
            header = source.read(12)
        if header[:4] not in (b'RIFF', b'RF64') or header[8:12] != b'WAVE':
            raise
        with tempfile.TemporaryDirectory(prefix='soneam-analysis-') as directory:
            normalized = os.path.join(directory, 'analysis.wav')
            result = subprocess.run([
                'ffmpeg', '-nostdin', '-v', 'error', '-y', '-i', str(path),
                '-map', '0:a:0', '-c:a', 'copy', '-map_metadata', '-1',
                '-rf64', 'auto', normalized,
            ], capture_output=True, timeout=120)
            if result.returncode:
                raise RuntimeError('WAV analysis container rewrite failed')
            with sf.SoundFile(normalized) as audio:
                yield audio
        return
    with audio:
        yield audio


def analyze_dynamic_range(input_path):
    """Approximate TT-style Dynamic Range from peak/RMS of the loudest blocks.

    Splits the file into ~3-second blocks, computes peak and RMS per block,
    sorts by RMS, takes the top 20%, and returns 20·log10(mean(peak)/mean(RMS)) for
    that subset. The result is in dB; typical commercial masters land in
    the 4–10 dB range; dynamic mixes show 12–20 dB.

    Not a strict TT-DR (which has additional smoothing & 3 s vs 20 s rules
    that vary by spec version), but close enough for engineer guidance and
    free since we already read the file for peaks generation.
    """
    import math
    log = logging.getLogger('soneam.dr')
    try:
        import soundfile as sf
        import numpy as np
    except ImportError:
        return None

    try:
        with _analysis_soundfile(input_path) as f:
            sr = f.samplerate
            channels = f.channels
            if sr <= 0 or channels <= 0:
                return None
            block = max(1, int(sr * 3))                # 3 s blocks
            peaks, rmses = [], []
            while True:
                chunk = f.read(block, dtype='float32', always_2d=True)
                if chunk.shape[0] == 0:
                    break
                # Combine channels: max(abs) for peak, sqrt(mean(square)) for RMS.
                p = float(np.max(np.abs(chunk)))
                r = float(np.sqrt(np.mean(chunk ** 2)))
                if p > 0 and r > 0:
                    peaks.append(p); rmses.append(r)
                if chunk.shape[0] < block:
                    break

            if not rmses:
                return None
            # Top 20% by RMS (the loudest stretches drive perceived dynamics).
            n = max(1, len(rmses) // 5)
            order = np.argsort(rmses)[::-1][:n]
            top_p = float(np.mean([peaks[i] for i in order]))
            top_r = float(np.mean([rmses[i] for i in order]))
            if top_r <= 0 or top_p <= 0:
                return None
            return 20.0 * math.log10(top_p / top_r)
    except Exception as e:
        log.warning('DR analysis failed for %s: %s', os.path.basename(input_path), e)
        return None


def generate_track_peaks(file_path, n=2048):
    """Return waveform buckets, sample peak and a 4× FFT true-peak estimate.

    Snapshot of the estimator used for the study. ffmpeg ebur128 is recorded
    independently; neither estimator is treated as evidence of audibility.
    """
    import math
    log = logging.getLogger('soneam.peaks')

    try:
        import soundfile as sf
        import numpy as np
    except ImportError:
        log.warning('soundfile/numpy not installed; cannot generate peaks')
        return None, None, None, None

    try:
        with _analysis_soundfile(file_path) as f:
            channels = f.channels
            total_frames = f.frames
            sr = f.samplerate
            if total_frames == 0 or sr <= 0:
                return None, None, None, None
            bucket_size = max(1, total_frames // n)
            actual_n = total_frames // bucket_size
            n_ch = min(channels, 2)

            peaks = [[] for _ in range(n_ch)]
            tp_bucket = []           # max true peak per bucket (linear, max across L/R)
            global_peak = 0.0
            global_tp = 0.0
            UP = 4  # 4× oversample for ITU-R BS.1770-style true peak.
            # Context samples read on each side of a bucket and upsampled with it,
            # then DISCARDED. A per-bucket FFT zero-pad upsample treats the bucket
            # as periodic; the wrap-around discontinuity turns into Gibbs ringing
            # that inflated the true peak by ~2 dB. Upsampling with margins and
            # keeping only the bucket interior removes that artifact.
            MARGIN = 128

            for i in range(actual_n):
                start = i * bucket_size
                rstart = max(0, start - MARGIN)
                pre = start - rstart                       # left context actually available
                f.seek(rstart)
                ext = f.read(bucket_size + pre + MARGIN, dtype='float32', always_2d=True)
                if ext.shape[0] == 0:
                    break
                blen = min(bucket_size, ext.shape[0] - pre)  # real bucket length within ext
                if blen <= 0:
                    break
                real = ext[pre:pre + blen]

                max_abs = np.max(np.abs(real), axis=0)
                for ch_idx in range(n_ch):
                    peaks[ch_idx].append(float(max_abs[ch_idx]))
                local_peak = float(np.max(max_abs[:n_ch]))
                if local_peak > global_peak:
                    global_peak = local_peak

                # True peak per bucket: upsample the EXTENDED block, read the peak
                # only from the real-bucket interior (margins discarded). Floor at
                # the bucket's sample peak (TP can never be below it).
                local_tp = local_peak
                if ext.shape[0] >= 8:
                    lo, hi = pre * UP, (pre + blen) * UP
                    for ch in range(n_ch):
                        col = ext[:, ch].astype(np.float64, copy=False)
                        spec = np.fft.rfft(col)
                        padded = np.zeros(col.size * UP // 2 + 1, dtype=complex)
                        padded[:spec.shape[0]] = spec
                        up = np.fft.irfft(padded, n=col.size * UP) * UP
                        p = float(np.max(np.abs(up[lo:hi])))
                        if p > local_tp:
                            local_tp = p
                    if local_tp > global_tp:
                        global_tp = local_tp
                tp_bucket.append(round(local_tp, 4))

            global_tp = max(global_tp, global_peak)   # TP ≥ sample peak, always

            result = {
                'l': peaks[0],
                'r': peaks[1] if channels >= 2 else None,
                'n': len(peaks[0]),
            }
            tp_result = {'t': tp_bucket, 'n': len(tp_bucket)} if tp_bucket else None
            sample_peak_dbfs = 20.0 * math.log10(global_peak) if global_peak > 0 else -120.0
            true_peak_dbtp   = 20.0 * math.log10(global_tp)   if global_tp   > 0 else -120.0
            return result, sample_peak_dbfs, true_peak_dbtp, tp_result
    except Exception as e:
        log.warning('peaks generation failed for %s: %s', os.path.basename(file_path), e)
        return None, None, None, None