the shepard tone script

#!/usr/bin/env python3
"""
the risset eternal accelerando: build it, then measure whether it actually accelerates.

    ==========================================================================
     hi, i'm jonas. i build audio plugins on my own, in denmark.

       the plugins    https://kernaudio.io
                      no ilok, no subscription, no account, and they keep
                      working offline forever. CHECK is free.

       the videos     https://youtube.com/@kernaudioio
                      subscribe if you want the next one.

       field notes    https://kernaudio.io/field-notes
                      one email a week. mixing and dsp, with the working
                      out shown. no filler.
    ==========================================================================

this is the script that builds every demo in the video, and it is the one i promised at the end.
run it, change the numbers, break it. that is the point of putting it here.

the third act of the film rests on one claim: the thing does not speed up forever, because it is
a loop. that is my own measurement rather than something i read, so this file has to take the
measurement instead of asserting it. `--measure` is where that happens.

the construction follows dan stowell, "scheduling and composing with risset eternal accelerando
rhythms", icmc 2011 pp.474-477. parallel streams of one pattern at tempos an octave apart, each
stream's amplitude set by a cosine window in log rate, so the top stream fades out as a new slow
one fades in underneath. his observation is the reason it does not turn to mush at the top: "the
density in a risset rhythm stays approximately constant irrespective of accelerando or
decelerando."

one figure you will see quoted elsewhere and will not find here: 120 to 15,360 bpm. it traces to
a press release about a single seven minute track, it is 120 x 2^7 done on a calculator, and it
appears nowhere in the peer reviewed paper. it is arithmetic, not a measurement, so i do not
print it.

stdlib only. no numpy, no scipy, nothing to install.

    python3 shepard-tone.py --selftest      prove the onset counter works before trusting it
    python3 shepard-tone.py --measure       the measurement the third act rests on
    python3 shepard-tone.py --render        write the demo wavs

(in this repo the file is `stimuli.py`; it is published on kernaudio.io as `shepard-tone.py`,
which is the name anyone who downloads it will have, so the usage lines use that one.)
"""
from __future__ import annotations

import argparse
import json
import array
import math
import os
import sys
import wave

SR = 48000
OUT = os.path.dirname(os.path.abspath(__file__))

# one cycle. after TAU seconds every stream has doubled and taken its neighbour's place, so the
# state is identical to the start. this is the loop length, and it is the whole point of act three.
TAU = 12.0
STREAM_LO = -5          # streams BELOW the window, waiting to fade in as they accelerate
STREAM_HI = 6
BASE_RATE = 2.0          # events/sec of the slowest stream at t=0
CENTRE_LOG_RATE = math.log2(4.0)   # the loud middle of the tempo window
WINDOW_OCT = 2.2         # half-width of the cosine amplitude window, in octaves of rate


def stream_rate(v: int, t: float) -> float:
    """Instantaneous event rate of stream v at time t. Doubles every TAU seconds."""
    return BASE_RATE * (2.0 ** (v + t / TAU))


def stream_amp(rate: float) -> float:
    """Cosine window in LOG-rate. Silent at both ends, loud in the middle."""
    d = (math.log2(rate) - CENTRE_LOG_RATE) / WINDOW_OCT
    return 0.0 if abs(d) >= 1.0 else 0.5 + 0.5 * math.cos(math.pi * d)


def click(n: int, f0: float = 1800.0, decay: float = 260.0) -> list:
    return [math.sin(2 * math.pi * f0 * i / SR) * math.exp(-decay * i / SR) for i in range(n)]


def render(dur: float) -> list:
    """The full accelerando. Streams are scheduled by integrating rate over time."""
    n = int(SR * dur)
    out = [0.0] * n
    ck = click(int(SR * 0.05))
    # v has to extend below the window, not just above it. my first version ran v=0..5 and the
    # measurement caught the mistake: the rate drifted up across cycles instead of repeating,
    # because once the original streams accelerated past the window there was nothing slower left
    # to fade in underneath. a risset rhythm is periodic or it is just an accelerando. worth saying
    # plainly that the measurement found this, not me. i would have shipped it.
    for v in range(STREAM_LO, STREAM_HI):
        # phase advances at the stream's instantaneous rate; a click on each whole phase
        phase, t, i = 0.0, 0.0, 0
        dt = 1.0 / SR
        while i < n:
            r = stream_rate(v, t)
            phase += r * dt
            if phase >= 1.0:
                phase -= 1.0
                a = stream_amp(r)
                if a > 1e-4:
                    for k, s in enumerate(ck):
                        if i + k < n:
                            out[i + k] += a * s
            t += dt
            i += 1
    peak = max(abs(x) for x in out) or 1.0
    return [0.7 * x / peak for x in out]


# ---------------------------------------------------------------- the SHEPARD half (D1, D2)
#
# the loop is sample exact rather than close enough, and that is a decision with a reason behind
# it. the cold open plays this with no voice over it, and the third act's whole argument is that
# the thing is a loop. a demo that gives away its own loop point in the first twenty seconds kills
# the ending before the beginning is over.
#
# Partial k has instantaneous frequency f0 * 2^(k + t/TAU_SHEP), so its phase is
#     phi_k(t) = 2*pi * f0 * 2^k * (TAU_SHEP/ln2) * (2^(t/TAU_SHEP) - 1)
# At t = TAU_SHEP that is 2*pi * 2^k * (f0*TAU_SHEP/ln2). So if
#     M = f0 * TAU_SHEP / ln2
# is an integer, every partial closes on a whole number of cycles at once, and partial k lands
# exactly where partial k+1 started: same frequency, same amplitude, same phase. x(t+TAU) == x(t)
# to floating point. so i pick M first and derive f0 from it, which puts f0 at 27.4948 Hz, which is
# A0 to within a twentieth of a cent. TAU_SHEP = 12.0 s also makes the loop a whole 576,000
# samples, which is tidy.
SHEP_M = 476                                   # integer cycles of the lowest partial per loop
TAU_SHEP = 12.0                                # one octave per loop; 576,000 samples at 48 kHz
SHEP_F0 = SHEP_M * math.log(2) / TAU_SHEP      # 27.4948 Hz, derived rather than typed
SHEP_SPAN = 10                                 # bell is zero at p=0 and p=SHEP_SPAN
SHEP_N = 10                                    # ten copies of it, which is what i say in the video
NYQ_P = math.log2(0.45 * SR / SHEP_F0)         # hard-zero above this, so nothing ever aliases


def shep_amp(p: float) -> float:
    """the stationary bell over log frequency. zero at both ends.

    p is position in octaves above f0. the window does not move. the content moves through it,
    and that is the sentence the whole video is built to make unmissable.
    """
    if p <= 0.0 or p >= SHEP_SPAN or p >= NYQ_P:
        return 0.0
    d = (p - SHEP_SPAN / 2.0) / (SHEP_SPAN / 2.0)
    return 0.5 + 0.5 * math.cos(math.pi * d)


def shepard(dur: float, only: int | None = None, norm: float | None = None) -> list:
    """the endless rise. `only` solos a single partial, which is how the second demo exposes it.

    `norm` exists because peak normalising two arms separately destroys the relationship between
    them, and that relationship is the demo. the second demo plays one copy and then all ten.
    normalised independently, the solo is a pure sine with a crest factor around 1.41, so it came
    out 2.4 dB louder than the ten partial stack, which has a higher crest and therefore loses rms
    at the same peak. the demo was saying one layer is louder than the sum that contains it. that
    is false, and it is the opposite of what "the layers separated" is supposed to show.

    pass the same divisor for both arms and the relationship survives. returns its own peak when
    `norm` is None, so the caller can capture it and reuse it.""" 
    n = int(SR * dur)
    out = [0.0] * n
    c = TAU_SHEP / math.log(2.0)
    # render one partial below the window. my first version ran k=0..9 and the loop did not close,
    # worst delta 2.6e-2. at t+TAU the set {0..9} has become {1..10}, so partial 0 leaves and
    # nothing arrives underneath it, and partial 0 is not silent on its way out: it is still at
    # amplitude 0.095 by p=1. extending down to k=-1 gives the departing partial a replacement that
    # is genuinely inaudible, because p in [-1,0) means amplitude exactly zero.
    # k=-1 only closes its phase if M is even, and the selftest asserts that.
    # this is the same mistake as the one in the rhythm half above: a periodic construction needs
    # streams below the window, not just above it. i made it twice, in two different halves.
    for k in range(-1, SHEP_N):
        if only is not None and k != only:
            continue
        base = 2.0 * math.pi * SHEP_F0 * (2.0 ** k) * c
        for i in range(n):
            u = 2.0 ** ((i / SR) / TAU_SHEP)
            a = shep_amp(k + (i / SR) / TAU_SHEP)
            if a > 1e-6:
                out[i] += a * math.sin(base * (u - 1.0))
    peak = norm if norm is not None else (max((abs(x) for x in out), default=0.0) or 1.0)
    return [0.7 * x / peak for x in out]


def goertzel(x: list, f: float, lo: int, hi: int) -> float:
    """amplitude at one frequency over x[lo:hi]. stdlib, and cheap enough to run per partial.

    this is a goertzel rather than an fft because every claim in the video is about named
    frequencies, the partials themselves, not about the shape of a spectrum. and because i wanted
    this file to run with nothing installed.
    """
    n = hi - lo
    if n <= 0:
        return 0.0
    w = 2.0 * math.pi * f / SR
    cw, sw = math.cos(w), math.sin(w)
    coeff = 2.0 * cw
    s1 = s2 = 0.0
    wsum = 0.0
    for i in range(lo, hi):
        # Hann, so a partial that is sliding does not smear across the whole analysis band
        win = 0.5 - 0.5 * math.cos(2.0 * math.pi * (i - lo) / n)
        wsum += win
        s0 = x[i] * win + coeff * s1 - s2
        s2, s1 = s1, s0
    re, im = s1 - s2 * cw, s2 * sw
    # normalise by the window's coherent gain, not by n. dividing by n reports a hann windowed
    # 0.400 tone as 0.200, because the window sums to n/2, so every amplitude came out exactly
    # half. the selftest planted a known amplitude and caught it. by eye the numbers looked
    # completely plausible, which is the problem with checking by eye.
    return 2.0 * math.hypot(re, im) / (wsum or 1.0)


def partial_levels(x: list, t: float, win: float = 0.35) -> list:
    """Measured amplitude of each of the ten partials at time t, read off the DELIVERED samples."""
    lo = max(0, int(SR * (t - win / 2)))
    hi = min(len(x), int(SR * (t + win / 2)))
    p0 = t / TAU_SHEP
    out = []
    for k in range(SHEP_N):
        f = SHEP_F0 * (2.0 ** (k + p0))
        out.append(goertzel(x, f, lo, hi) if f < 0.45 * SR else 0.0)
    return out


def centroid(x: list, t: float) -> float:
    """spectral centroid in octaves above f0, measured from the samples rather than from the model.

    in the video i say the ladder never runs out. the measurable form of that is that the centre
    of mass comes back to where it started after one loop instead of climbing away with the notes.
    """
    lv = partial_levels(x, t)
    p0 = t / TAU_SHEP
    num = sum(a * (k + p0) for k, a in enumerate(lv))
    den = sum(lv)
    return num / den if den > 1e-12 else 0.0


# ---------------------------------------------------------------- the DROPPED-SUBDIVISION half (D4)

def dropped(dur: float, tau: float = TAU_SHEP) -> list:
    """one pulse train whose rate doubles over `tau`, with every second hit fading out as it goes.

    this is the rhythm version of the same trick. by the time the rate has doubled, half the hits
    have gone, so the events per second land back where they started.

    "dropped subdivisions" is jake lizzio's phrase for this, not established terminology, so i
    credit him for it rather than using it as though it were a standard term.
    """
    n = int(SR * dur)
    out = [0.0] * n
    ck = click(int(SR * 0.05))
    rate0 = 4.0
    phase, i, t = 0.0, 0, 0.0
    dt = 1.0 / SR
    idx = 0
    while i < n:
        # both the rate and the fade run on the cycle phase, not on absolute time. my first version
        # accelerated on absolute t while the fade reset at each cycle boundary, so at t=TAU the
        # grid was at double rate with every hit back at full amplitude. it measured 4.50/s going
        # to 8.50/s, which is plainly not landing where it began. within one cycle the grid
        # doubles and the odd hits fade out, so the events per second come out where they started
        # and the next cycle inherits them as its full set.
        tp = t % tau
        r = rate0 * (2.0 ** (tp / tau))
        phase += r * dt
        if phase >= 1.0:
            phase -= 1.0
            # odd-numbered hits fade linearly to silence across one doubling
            a = 1.0 if idx % 2 == 0 else max(0.0, 1.0 - tp / tau)
            if a > 1e-4:
                for j, s in enumerate(ck):
                    if i + j < n:
                        out[i + j] += a * s
            idx += 1
        t += dt
        i += 1
    peak = max((abs(v) for v in out), default=0.0) or 1.0
    return [0.7 * v / peak for v in out]


# ---------------------------------------------------------------- the HEMIOLA (D5)
#
# one arm, played once, and this is the file that decides it rather than leaving it to whatever
# the picture happens to draw.
#
# in the video i say: same six hits, nothing changes in the audio, and the music flips between the
# two feels. if i played two arms here the audio would change, and that sentence would be false as
# spoken. so it is one bar of six subdivisions carrying both feels at once, repeated unchanged,
# and the picture does the re-pointing. the ear does the rest.
HEM_SUB = 0.25          # seconds per subdivision
HEM_BARS = 6


def hemiola(bars: int = HEM_BARS) -> list:
    """six subdivisions per bar. a low voice on the two feel, a high voice on the three feel, both
    at once."""
    n = int(SR * HEM_SUB * 6 * bars)
    out = [0.0] * n
    low = click(int(SR * 0.09), f0=220.0, decay=90.0)
    high = click(int(SR * 0.05), f0=880.0, decay=200.0)
    for b in range(bars):
        for s in range(6):
            at = int(SR * HEM_SUB * (6 * b + s))
            if s % 3 == 0:                      # two groups of three
                for j, v in enumerate(low):
                    if at + j < n:
                        out[at + j] += 0.9 * v
            if s % 2 == 0:                      # three groups of two
                for j, v in enumerate(high):
                    if at + j < n:
                        out[at + j] += 0.6 * v
    peak = max((abs(v) for v in out), default=0.0) or 1.0
    return [0.7 * v / peak for v in out]


# ---------------------------------------------------------------- measurement

def onsets(x: list, lo: int, hi: int, thresh: float = 0.12, refractory: float = 0.010) -> int:
    """count clicks in x[lo:hi] by rectified envelope threshold crossing."""
    a = math.exp(-2 * math.pi * 400.0 / SR)
    y = 0.0
    env = []
    for s in x[lo:hi]:
        y = (1 - a) * abs(s) + a * y
        env.append(y)
    peak = max(env) or 1.0
    gap = int(SR * refractory)
    count, last = 0, -gap
    for i, v in enumerate(env):
        if v > thresh * peak and i - last >= gap:
            count += 1
            last = i
    return count


def rate_at(x: list, t: float, win: float = 2.0) -> float:
    """measured events per second in a window centred on t."""
    lo = max(0, int(SR * (t - win / 2)))
    hi = min(len(x), int(SR * (t + win / 2)))
    return onsets(x, lo, hi) / ((hi - lo) / SR)


def phase_aligned_pair(rates: list, tol: float = 1e-6):
    """pick the widest pair of sampled points that are a whole number of loop cycles apart.

    this function exists because the summary it replaced was wrong, and wrong in an expensive way:
    it printed an instruction to whoever ran it. the old code took `rates[1]` and `rates[-1]`,
    which with the default sample points is t=6s and t=33s. that is 27 seconds against a 12 second
    loop, so they are different phases of the same cycle. it then announced that the later point
    was faster by 0.50/s, that the sign was the other way round, and that i should say so on
    camera.

    so it was not a stale comment. it was a measuring device telling me my own climax was
    backwards. in the video i compare t=3s and t=27s, which is 24 seconds, exactly two whole
    cycles, and both measure 9.50/s. the take was right and this file's headline disagreed with it.

    the reason is that the instantaneous rate oscillates between about 6.5 and 9.5 inside every
    cycle. that is not noise and it is not a fault. it is what a windowed onset count does to a
    construction whose streams fade in and out. so a rate difference between two arbitrary points
    tells you where in the cycle you looked, not whether the thing accelerated. only a whole
    number of cycles apart is a comparison at all.

    returns (a, b) as (t, rate) tuples, or (None, None) if no aligned pair exists, in which case
    the caller has to refuse to state a sign rather than pick the nearest thing and guess.
    """
    best = (None, None)
    widest = 0.0
    for i, a in enumerate(rates):
        for b in rates[i + 1:]:
            dt = b[0] - a[0]
            if abs(dt - round(dt / TAU) * TAU) <= tol and round(dt / TAU) >= 1 and dt > widest:
                widest, best = dt, (a, b)
    return best


def measure() -> int:
    dur = TAU * 3 + 4.0
    x = render(dur)
    print(f"Risset accelerando: {dur:.1f}s, streams v={STREAM_LO}..{STREAM_HI-1}, loop TAU = {TAU}s\n")
    print(f"{'t (s)':>8}{'measured events/sec':>22}")
    pts = [3.0, 6.0, 9.0, 15.0, 18.0, 21.0, 27.0, 30.0, 33.0]
    rates = []
    for t in pts:
        if t + 1.0 > dur:
            continue
        r = rate_at(x, t)
        rates.append((t, r))
        print(f"{t:>8.1f}{r:>22.2f}")

    early = [r for t, r in rates if t < TAU]
    late = [r for t, r in rates if t > 2 * TAU]
    print(f"\nmean rate in the first cycle : {sum(early)/len(early):.2f} events/sec")
    print(f"mean rate in the third cycle : {sum(late)/len(late):.2f} events/sec")

    # the claim the third act rests on, stated by the measurement rather than by me
    a, b = phase_aligned_pair(rates)
    if a is None:
        # refusing is the correct output here. see the note on phase alignment above.
        print("\ncannot summarise: no two sampled points are a whole number of cycles apart, "
              f"so every available pair compares DIFFERENT PHASES of the {TAU:.0f}s loop.")
        print("     A rate difference between two phases says nothing about whether the thing "
              "accelerates. Sample points TAU apart and re-run.")
        return 1

    cycles = round((b[0] - a[0]) / TAU)
    print(f"\nt={a[0]:.0f}s and t={b[0]:.0f}s are {cycles} full cycle(s) apart "
          f"({cycles}x TAU={TAU:.0f}s), which is the same phase of the loop.")
    print(f"     t={a[0]:.0f}s measures {a[1]:.2f}/s; t={b[0]:.0f}s measures {b[1]:.2f}/s")
    if abs(b[1] - a[1]) <= 0.05:
        print("     the same, to the resolution of this counter. it did not accelerate at all.")
    elif b[1] < a[1]:
        print(f"     the later point is slower, by {a[1]-b[1]:.2f} events/sec.")
    else:
        print(f"     the later point is faster by {b[1]-a[1]:.2f}/s, so the sign is the other way. "
              "Say that on camera, not the assumption.")
    print("\nwhichever it is, that is what i say in the video. measure first, write second.")
    return 0


# ---------------------------------------------------------------- selftest

def selftest() -> int:
    global SHEP_F0          # the defect fixture below detunes this on purpose, then puts it back
    fails = []
    # 1. the counter has to recover a known click rate before i trust anything it says
    for want in (3.0, 8.0, 16.0):
        n = int(SR * 4.0)
        x = [0.0] * n
        ck = click(int(SR * 0.05))
        step = int(SR / want)
        for start in range(0, n, step):
            for k, s in enumerate(ck):
                if start + k < n:
                    x[start + k] += s
        got = rate_at(x, 2.0, win=2.0)
        if abs(got - want) > 0.6:
            fails.append(f"counter: planted {want}/s, measured {got:.2f}")
        else:
            print(f"  [ok] counter recovers a planted {want:.0f}/s train as {got:.2f}")

    # 2. silence has to measure zero. a counter that finds events in nothing will find them anywhere
    if rate_at([0.0] * int(SR * 3), 1.5) > 0.4:
        fails.append("the counter finds events in silence")
    else:
        print("  [ok] silence measures ~0 events/sec")

    # 3. the loop invariant. the construction has to repeat with period TAU
    a = stream_rate(0, 0.0) * 2
    b = stream_rate(0, TAU)
    if abs(a - b) > 1e-9:
        fails.append(f"stream rate does not double over TAU ({a} vs {b})")
    else:
        print(f"  [ok] every stream doubles over TAU={TAU}s, so the state repeats")

    # 4. the summary has to compare phase aligned points, or refuse. three checks, and the second
    #    one is a fixture built to fail the behaviour this file had until i fixed it.
    aligned = [(3.0, 9.50), (6.0, 6.50), (15.0, 9.50), (27.0, 9.50), (33.0, 7.00)]
    a, b = phase_aligned_pair(aligned)
    if a is None or (a[0], b[0]) != (3.0, 27.0):
        fails.append(f"phase picker chose {a}..{b}, wanted t=3 -> t=27 (the widest whole-cycle pair)")
    else:
        print("  [ok] phase picker chooses t=3 -> t=27, two whole cycles apart, the pair i quote")

    # the old code took rates[1] and rates[-1], which is t=6 and t=33. 27 seconds against a 12
    # second loop is not a whole number of cycles, and asserting it here is what stops that
    # coming back.
    dt_old = 33.0 - 6.0
    if abs(dt_old - round(dt_old / TAU) * TAU) <= 1e-6:
        fails.append(f"t=6 -> t=33 ({dt_old}s) is being treated as phase-aligned against TAU={TAU}")
    else:
        print(f"  [ok] the old pair t=6 -> t=33 ({dt_old:.0f}s) is correctly not phase-aligned")

    # and a set with no aligned pair at all has to return None, so the caller refuses rather than
    # picking the nearest thing and asserting a direction from it.
    if phase_aligned_pair([(1.0, 5.0), (4.0, 9.0), (20.0, 6.0)])[0] is not None:
        fails.append("phase picker returned a pair from a set with no whole-cycle spacing")
    else:
        print("  [ok] a set with no whole-cycle pair returns None, so the summary has to refuse")

    # 5. the pitch side instruments. same discipline as the onset counter above: each one gets a
    #    planted signal whose answer i already know, and then the half that actually matters, a
    #    planted defect it has to go red on. a checker that has never been seen to fail is not
    #    evidence that anything is right. i have shipped blind checkers before and it is a
    #    horrible way to find out.
    # the loop closes only if M is an integer AND even (the k=-1 partial closes on pi*M)
    if SHEP_M % 2 != 0:
        fails.append(f"SHEP_M={SHEP_M} is odd, so the partial below the window closes on a half cycle")
    elif abs(SHEP_F0 * TAU_SHEP / math.log(2) - SHEP_M) > 1e-9:
        fails.append("SHEP_F0 is not derived from SHEP_M, so no partial closes on a whole cycle")
    elif not float(SR * TAU_SHEP).is_integer():
        fails.append(f"a loop is {SR * TAU_SHEP} samples, not a whole number")
    else:
        print(f"  [ok] M={SHEP_M} even, f0 derived ({SHEP_F0:.4f} Hz), loop = {int(SR*TAU_SHEP)} samples")

    n = int(SR * 0.5)
    tone = [0.4 * math.sin(2 * math.pi * 1000.0 * i / SR) for i in range(n)]
    got = goertzel(tone, 1000.0, 0, n)
    if abs(got - 0.4) > 0.02:
        fails.append(f"goertzel: planted a 0.400 tone at 1 kHz, measured {got:.3f}")
    else:
        print(f"  [ok] goertzel recovers a planted 0.400 amplitude at 1 kHz as {got:.3f}")
    off = goertzel(tone, 3000.0, 0, n)
    if off > 0.02:
        fails.append(f"goertzel finds {off:.3f} at 3 kHz in a signal that only contains 1 kHz")
    else:
        print(f"  [ok] goertzel reads {off:.4f} at a frequency that is not there")

    # 5b. the defect fixtures. every one of these has to be rejected.
    stack = shepard(2.5)
    if not assert_d2(stack, stack):
        fails.append("the separation check accepted the full stack as a solo, so it is blind")
    else:
        print("  [ok] assert_d2 rejects a stack passed off as a solo")

    broken = hemiola(bars=3)
    bar = int(SR * HEM_SUB * 6)
    broken[bar + 1000] += 0.25          # one sample of difference, one bar in
    if not assert_d5(broken):
        fails.append("the hemiola check accepted bars that differ, so 'nothing changes' is unguarded")
    else:
        print("  [ok] assert_d5 rejects a bar with a single altered sample")

    # a shepard tone whose M is not an integer does not close, and the seamlessness check has to say so
    # (declared at the top of selftest)
    keep = SHEP_F0
    try:
        SHEP_F0 = keep * 1.004                      # M = 476 -> 477.9, phase no longer closes
        if not assert_d1(shepard(TAU_SHEP * 2)):
            fails.append("the seamlessness check accepted a loop whose phase does not close")
        else:
            print("  [ok] assert_d1 rejects a Shepard whose partials do not close on a whole cycle")
    finally:
        SHEP_F0 = keep

    if fails:
        print("\nselftest failed:")
        for f in fails:
            print("  -", f)
        return 1
    print("\nselftest passed")
    return 0


def write_wav(name: str, x: list) -> None:
    d = array.array("h", (int(max(-1.0, min(1.0, s)) * 32767) for s in x))
    with wave.open(os.path.join(OUT, name), "wb") as w:
        w.setnchannels(1)
        w.setsampwidth(2)
        w.setframerate(SR)
        w.writeframes(d.tobytes())
    print("wrote", name)


# ---------------------------------------------------------------- every demo is a claim
#
# every demo in the video is a claim, so every one of them gets measured against the sentence i
# say over it. each function below names the line it is defending and returns a list of failures.
# nothing is written to disk until all of those lists are empty.

MAX_TRAILING_SILENCE = 0.9      # a demo has to fill the window it was given, or the picture keeps moving over silence


def _trailing_silence(x: list, thresh: float = 1e-3) -> float:
    for i in range(len(x) - 1, -1, -1):
        if abs(x[i]) > thresh:
            return (len(x) - 1 - i) / SR
    return len(x) / SR


def assert_d1(x: list) -> list:
    """the cold open. three claims, three measurements."""
    f = []
    n_loop = int(SR * TAU_SHEP)

    # (1) seamless. the loop point has to be inaudible, and the surest way is for it not to exist.
    if len(x) >= 2 * n_loop:
        # compare every sample. my first version strided by 97, and a fixture proved why that is
        # not a check: a planted one sample defect at offset 1000 fell between the strides and the
        # harness reported clean. a hole in the sampling is a hole in the guarantee.
        a, b = x[:n_loop], x[n_loop:2 * n_loop]
        worst = 0.0 if a == b else max(abs(p - q) for p, q in zip(a, b))
        # a tolerance rather than equality, and the number is chosen against the medium. demanding
        # exact float equality failed the in memory buffer at 8.4e-11, which is accumulated rounding
        # in the phase integral, while the delivered int16 file compared bit identical, because one
        # int16 step is 3.05e-5, which is six orders of magnitude larger. a guarantee tighter than
        # the format can carry is a false alarm, not a stronger guarantee. 1e-9 is still about
        if worst > 1e-9:
            f.append(f"D1 loop is not sample-exact: worst |x[i]-x[i+TAU]| = {worst:.3e}")
        else:
            print(f"  [ok] D1 loops sample-exactly, every sample compared (worst {worst:.1e}, "
                  f"one int16 step is 3.1e-5), so the joint does not exist")
    else:
        f.append("D1 rendered shorter than two loops, so seamlessness was not tested")

    # (2) in the video: loud in the middle of your hearing, silent at both ends.
    #     that is a claim about the window, so i measure the window at its ends, where it is
    #     exactly zero by construction, then state how far down the outermost sounding partial is.
    #     my first version compared the outermost partials against the loudest and demanded
    #     -60 dB. that measured the wrong thing. at t=3 the lowest partial sits at p=0.25 and the
    #     highest at p=9.25, and neither of those is an end of the window. the window is silent at
    #     its ends. the partials near them are quiet, not silent. and -60 was a number i picked.
    if shep_amp(0.0) != 0.0 or shep_amp(float(SHEP_SPAN)) != 0.0:
        f.append("D1 window is not zero at its ends, so 'silent at both ends' is false")
    else:
        lv = partial_levels(x, 3.0)
        mid = max(lv)
        db = 20 * math.log10((max(lv[0], lv[-1]) + 1e-15) / (mid + 1e-15))
        print(f"  [ok] D1 window is exactly zero at both ends; outermost partial sits {db:.1f} dB down")

    # (3) the illusion itself, and i had this assertion backwards the first time.
    #     i demanded that the centroid move, and treated a still one as a fault. the still one is
    #     the whole point. in the video i say the ambiguity is not in the signal, it is in what
    #     you were tracking. every partial climbs, and the centre of mass of the spectrum does
    #     not go anywhere. that pair is the thesis, so that pair is what i measure.
    cs = [centroid(x, 1.0 + i * TAU_SHEP / 8) for i in range(9)]
    drift = max(cs) - min(cs)
    if drift > 0.15:
        f.append(f"D1 centroid drifts {drift:.3f} octaves across a loop, so it is rising, not circling")
    else:
        print(f"  [ok] D1 centroid holds within {drift:.3f} octaves across a full loop, so it goes nowhere")

    # and meanwhile a single partial demonstrably climbs. read off the delivered samples at the
    # frequency a rising partial should occupy, half a loop apart.
    a = goertzel(x, SHEP_F0 * 2 ** (5 + 1.0 / TAU_SHEP), int(SR * 0.8), int(SR * 1.2))
    b = goertzel(x, SHEP_F0 * 2 ** (5 + (1.0 + TAU_SHEP / 2) / TAU_SHEP),
                 int(SR * (0.8 + TAU_SHEP / 2)), int(SR * (1.2 + TAU_SHEP / 2)))
    if min(a, b) < 0.2 * max(partial_levels(x, 1.0)):
        f.append("D1 energy did not follow the rising partial, so nothing is actually climbing")
    else:
        print(f"  [ok] D1 partial 5 is still loud half an octave higher ({a:.3f} -> {b:.3f}): it climbs")
    return f


def assert_d2(solo: list, full: list) -> list:
    """the layers separated, and audible.

    i have shipped a masked beep as "audible" before, because a dB figure in the code meant
    something other than what i thought. so this measures the delivered samples and states the
separation as a number instead of trusting the word.
    """
    f = []
    s = partial_levels(solo, 2.0)
    top = max(s)
    others = sorted(s, reverse=True)[1]
    sep = 20 * math.log10((others + 1e-15) / (top + 1e-15))
    if sep > -20.0:
        f.append(f"D2 solo is not solo: next partial only {sep:.1f} dB down, wanted <= -20")
    else:
        print(f"  [ok] D2 soloed layer stands {abs(sep):.1f} dB clear of every other partial")

    lv = partial_levels(full, 2.0)
    loud = max(lv)
    audible = sum(1 for a in lv if a > loud * 0.1)      # within 20 dB of the loudest
    if audible < 3:
        f.append(f"D2 stack has only {audible} partial(s) within 20 dB, so the layers are not audible")
    else:
        print(f"  [ok] D2 stack carries {audible} partials within 20 dB of the loudest")
    return f


def assert_d4(x: list) -> list:
    """the dropped subdivisions. in the video: the number of events per second stays roughly
constant."""
    f = []
    a, b = rate_at(x, 1.5), rate_at(x, 1.5 + TAU_SHEP)
    if abs(b - a) > 1.2:
        f.append(f"D4 did not land where it began: {a:.2f}/s -> {b:.2f}/s one doubling later")
    else:
        print(f"  [ok] D4 measures {a:.2f}/s and {b:.2f}/s one full doubling apart, so it lands home")
    return f


def assert_d5(x: list) -> list:
    """the hemiola. in the video: same six hits, nothing changes in the audio."""
    f = []
    bar = int(SR * HEM_SUB * 6)
    if len(x) >= 3 * bar:
        # exact, every sample. see the note in the cold open check above. a strided version of this
        # passed a bar carrying a planted one sample difference.
        a, b = x[bar:2 * bar], x[2 * bar:3 * bar]
        if a != b:
            worst = max(abs(p - q) for p, q in zip(a, b))
            f.append(f"D5 bars are not identical (worst delta {worst:.3e}), so 'nothing changes' is false")
        else:
            print("  [ok] D5 bar 2 is sample-identical to bar 3, so nothing changes in the audio")
    # both feels have to actually be present, or the sentence i say over it is undemonstrated
    lo_hits = onsets([v if abs(v) > 0 else 0.0 for v in x], 0, bar * 2)
    if lo_hits < 6:
        f.append(f"D5 bar carries only {lo_hits} onsets across two bars, wanted >= 6")
    else:
        print(f"  [ok] D5 carries {lo_hits} onsets across two bars, so both feels are sounding")
    return f


def assert_d6(x: list) -> list:
    """the payoff. the two arms have to be the two moments i name out loud, and they have to
measure alike.

    this is the check that would have caught the old splice, and nothing else could have. the file
    was the right length, it filled its window, it had no trailing silence, and it passed every
    other check i own. the fault only existed in the relationship between the audio and the
sentence read over it.
    """
    f = []
    a, b = rate_at(x, 3.0), rate_at(x, 9.0)   # centre of arm one, centre of arm two (6 s in, which is t=27 in the source)
    if abs(a - b) > 0.05:
        f.append(f"D6 arms do not match: {a:.2f}/s and {b:.2f}/s, but i say both are the same")
    elif abs(a - 9.5) > 0.05:
        f.append(f"D6 arms both measure {a:.2f}/s, but i say nine and a half in the video")
    else:
        print(f"  [ok] D6 both arms measure {a:.2f}/s, the two moments i name, and they agree")
    return f


def guard_delivered(name: str, check) -> list:
    """re-open the written file and re-assert on the bytes that actually ship.

    this exists because i once had a verify pass while the delivered audio was 2 dB wrong. every
    check had measured intermediates rather than the file itself. this script post-processes an arm
    after its checks have run, since the payoff demo is a splice of two slices of another one,
    which is exactly the hole that lets it happen again.
    """
    with wave.open(os.path.join(OUT, name), "rb") as w:
        raw = w.readframes(w.getnframes())
    d = array.array("h")
    d.frombytes(raw)
    x = [v / 32767.0 for v in d]
    fails = [f"{name}: {m}" for m in check(x)]
    tail = _trailing_silence(x)
    if tail > MAX_TRAILING_SILENCE:
        fails.append(f"{name}: {tail:.2f}s of trailing silence, over the {MAX_TRAILING_SILENCE}s budget")
    return fails


def do_render() -> int:
    fails = []

    print("D3 / D6: the risset accelerando, and the splice that breaks it")
    x = render(TAU * 3 + 4.0)
    write_wav("DEMO_D3_RISSET.wav", x)
    # the splice is 0-6 and 24-30, and it used to be 0-6 and 30-36. that was a real fault on
    # THE FILM'S MOST IMPORTANT BEAT, and it is the JOIN class: the audio and the take were each
    # fine alone and disagreed with each other.
    # in the video i quote t=3 s and t=27 s, which is 24 s apart, two whole cycles, and that is
    # are comparable at all (`phase_aligned_pair` exists to enforce that). The old splice put the
    # second arm's centre at t=33, which is 30 s after t=3, so two and a half cycles, anti-phase. it
    # measures 7.00/s. So the viewer would have HEARD a slower second arm while the voice said both
    # measured nine and a half, on the beat the whole film is built to deliver.
    # Now each arm is centred on a moment the take actually names.
    d6 = x[:int(SR * 6)] + x[int(SR * 24):int(SR * 30)]
    fails += assert_d6(d6)
    write_wav("DEMO_D6_BREAK.wav", d6)

    print("\nD1: the endless rise (the cold open)")
    d1 = shepard(TAU_SHEP * 3)
    fails += assert_d1(d1)
    write_wav("DEMO_D1_SHEPARD.wav", d1)

    print("\nD2: the same tone with one layer soloed, then the stack")
    # one divisor for both arms, see `shepard(norm=)`. the stack sets it, because the stack is the
    # thing the solo is a part of.
    _raw_full = shepard(4.0, norm=1.0)
    _pk = max(abs(v) for v in _raw_full) / 0.7
    solo = shepard(4.0, only=5, norm=_pk)
    full = shepard(4.0, norm=_pk)
    fails += assert_d2(solo, full)
    write_wav("DEMO_D2_OCTAVES.wav", solo + full)

    print("\nD4: dropped subdivisions")
    d4 = dropped(TAU_SHEP + 3.0)
    fails += assert_d4(d4)
    write_wav("DEMO_D4_DROPPED.wav", d4)

    print("\nD5: the hemiola, one arm")
    d5 = hemiola()
    fails += assert_d5(d5)
    write_wav("DEMO_D5_HEMIOLA.wav", d5)

    # ---- and now the only measurement that counts: the bytes on disk
    print("\ndelivered-file guard: re-opening every wav that ships")
    fails += guard_delivered("DEMO_D1_SHEPARD.wav", assert_d1)
    fails += guard_delivered("DEMO_D4_DROPPED.wav", assert_d4)
    fails += guard_delivered("DEMO_D5_HEMIOLA.wav", assert_d5)
    fails += guard_delivered("DEMO_D6_BREAK.wav", assert_d6)
    for nm in ("DEMO_D2_OCTAVES.wav", "DEMO_D3_RISSET.wav"):
        fails += guard_delivered(nm, lambda _x: [])

    if fails:
        print("\nrender failed its own claims:")
        for m in fails:
            print("  -", m)
        return 1
    print("\nall six demos written, and every one measured against the sentence it has to prove")
    return 0


# ---------------------------------------------------------------- the listening windows
#
# a window's wav is its duration. the picture reads these files, so a demo cannot under-fill the
# window it was given. i have had that go wrong before: three demos stopped seconds before their
# reserved window ended, and the playhead kept sweeping over silence.
#
# the dropped subdivisions window is 12 s rather than the 8 s i first budgeted. its whole claim
# is that it lands where it began, and where it began is one full doubling away. an 8 s window
# would cut the cycle two thirds of the way through and show an accelerando that never returns,
# which is the opposite of the point. a longer demo window is teaching, not padding.
WINDOWS = [
    # 12.0 s, not the 20 s i originally wanted, and the number is not a compromise.
    # 20 s of tone with no voice put the first word at 21.6 s. the steepest drop on any video is the
    # first 30 seconds, so that is the most expensive stretch in the whole piece to spend on
    # silence.
    # TAU_SHEP is 12.0 s, so 12.0 s is exactly one loop. the cold open plays the whole cycle and
    # ends on the sample it started from, which is the thesis performed before a word is spoken.
    # it is also the shortest window that still contains a complete statement of the illusion, and
    # it moves the first word to 13.6 s. shorter would cut the loop mid cycle and give away
    # nothing. longer buys no new information at all, because after one cycle the tone repeats
    # exactly. nobody is any less fooled by one cycle than by two, so the extra 8 seconds of
    # silence bought nothing.
    ("W1_RISE",     "DEMO_D1_SHEPARD.wav",  0.0, 12.0, 0.40),
    ("W2_LAYERS",   "DEMO_D2_OCTAVES.wav",  0.0,  8.0, 0.02),
    ("W3_TIME",     "DEMO_D3_RISSET.wav",   0.0,  8.0, 0.02),
    ("W4_DROPPED",  "DEMO_D4_DROPPED.wav",  0.0, 12.0, 0.02),
    ("W5_HEMIOLA",  "DEMO_D5_HEMIOLA.wav",  0.0,  9.0, 0.02),
    ("W6_BREAK",    "DEMO_D6_BREAK.wav",    0.0, 12.0, 0.02),
]

# the arms are a contract with the picture, and two of them are decisions rather than facts.
# the picture refuses to draw if a window's display has a different number of rows than the
# window has distinct arms. that guard exists because a clamp once made a video name the wrong
# sound for months, and every check i had passed it. so this table is where "how many things is
# the viewer hearing" gets decided once, in the file that makes the sound, instead of drifting
# between the audio and the picture.
#
#   the hemiola is one arm on purpose. in the video i say: same six hits, nothing changes in the
#   audio, and the music flips between the two feels. two arms would make that sentence false as
#   spoken. so it is one arm, played once, and the picture does the re-pointing.
#
#
#   the break is two arms, and they are not a comparison of two sounds. they are two slices of
#   one file, six seconds from the start and six from near the end. that is why the picture draws
#   two panels with a visible break between them rather than one axis with a jumping playhead. a
#   playhead sweeping 0 to 36 s would spend 24 seconds crossing material nobody is hearing.
ARMS = {
    "W1_RISE":    [("a tone that climbs and never arrives", 0.0, 12.0)],
    "W2_LAYERS":  [("one layer, alone", 0.0, 4.0), ("all ten together", 4.0, 8.0)],
    "W3_TIME":    [("three tempos, crossfaded", 0.0, 8.0)],
    "W4_DROPPED": [("every second hit fading out", 0.0, 12.0)],
    "W5_HEMIOLA": [("one bar of six, both feels at once", 0.0, 9.0)],
    # the source times are the label. in the video i name "three seconds in" and "twenty seven
    # seconds in", so the two panels have to say which part of the run each one is, or you cannot
    # tell they are two moments of one recording rather than two different recordings.
    "W6_BREAK":   [("t = 0-6 s", 0.0, 6.0), ("t = 24-30 s", 6.0, 12.0)],
}


def do_windows() -> int:
    wdir = os.path.join(OUT, "windows")
    os.makedirs(wdir, exist_ok=True)
    fails = []
    for wid, src, t0, dur, fade in WINDOWS:
        with wave.open(os.path.join(OUT, src), "rb") as w:
            raw, sr = w.readframes(w.getnframes()), w.getframerate()
        d = array.array("h")
        d.frombytes(raw)
        x = [v / 32767.0 for v in d][int(t0 * sr):int((t0 + dur) * sr)]
        if len(x) < int(dur * sr) - 1:
            fails.append(f"{wid}: {src} is too short for a {dur}s window")
            continue
        # a de-click at both edges. the cold open gets a real fade in, since it opens the video
        nf = int(fade * sr)
        for i in range(nf):
            g = i / nf
            x[i] *= g
            x[-1 - i] *= g
        tail = _trailing_silence(x)
        if tail > MAX_TRAILING_SILENCE:
            fails.append(f"{wid}: {tail:.2f}s of trailing silence in a {dur}s window")
        d2 = array.array("h", (int(max(-1.0, min(1.0, s)) * 32767) for s in x))
        with wave.open(os.path.join(wdir, wid + ".wav"), "wb") as w:
            w.setnchannels(1)
            w.setsampwidth(2)
            w.setframerate(sr)
            w.writeframes(d2.tobytes())
        print(f"  wrote windows/{wid}.wav  {dur:5.2f}s  from {src}")
    # the manifest the picture reads. derived from ARMS above and never hand typed, so the picture's
    # row count and the audio's arm count cannot drift apart.
    man = {}
    for wid, src, t0, dur, fade in WINDOWS:
        arms = ARMS[wid]
        if abs(arms[-1][2] - dur) > 1e-6 or abs(arms[0][1]) > 1e-6:
            fails.append(f"{wid}: ARMS span {arms[0][1]}..{arms[-1][2]} but the window is 0..{dur}")
        man[wid] = {
            "caption": arms[0][0] if len(arms) == 1 else " / ".join(a[0] for a in arms),
            "duration": dur,
            "segments": [{"name": n, "start": a, "end": b} for n, a, b in arms],
        }
    with open(os.path.join(wdir, "windows.json"), "w") as fh:
        json.dump(man, fh, indent=1)

    if fails:
        print("\nwindows failed:")
        for m in fails:
            print("  -", m)
        return 1
    print(f"\n{len(WINDOWS)} windows written; each one's duration is its file, so none can under-fill")
    print(f"windows.json: {sum(len(v['segments']) for v in man.values())} arms across {len(man)} "
          f"windows (W5 deliberately 1, W6 deliberately 2)")
    return 0


def main() -> int:
    ap = argparse.ArgumentParser()
    ap.add_argument("--selftest", action="store_true")
    ap.add_argument("--measure", action="store_true")
    ap.add_argument("--render", action="store_true")
    ap.add_argument("--windows", action="store_true")
    a = ap.parse_args()
    if a.selftest:
        return selftest()
    if a.render:
        return do_render()
    if a.windows:
        return do_windows()
    return measure()


if __name__ == "__main__":
    sys.exit(main())