# SPDX-License-Identifier: MIT
# © 2026 DTMFRAME
"""Synthetic-signal demonstration for the article reverb-time-decay-predelay.

Questions it answers with numbers (no product, no licence, no recorded material):
  D1  Does the Schroeder backward-integration measurement used here recover a known RT60?
      (exponentially decaying noise with a set RT60 of 0.5, 1, 2 and 4 s)
  D2  Does pre-delay change the reverberation time? (same tail, pre-delay 0, 20, 80 ms)
  D3  If a "decay" control is a feedback amount rather than a time, does the reverberation
      time follow the room size? (toy 8-line feedback delay network, size x0.5, x1, x2;
      mode A sets each line's gain from a target T60, mode B uses one fixed gain)
  D4  Is the reverberation time one number across frequency? (the same network with
      one-pole absorption filters; octave-band T30 against the design curve)
  D5  Where does pre-delay sit relative to the early reflections? (two constructions of
      the same impulse response, pre-delay before everything vs before the late tail only)

It also writes the article's concept diagram, out/concept_timeline.svg (figures.py; no data,
only the numbers above drawn to scale).

Run:  ../.venv/bin/python demo.py      (writes out/ and out/record.json)
"""
from __future__ import annotations

import sys
from pathlib import Path

import numpy as np
from scipy import signal

HERE = Path(__file__).resolve().parent
sys.path.insert(0, str(HERE.parent / "common"))
import demo_common as dc  # noqa: E402
import figures  # noqa: E402
import matplotlib.pyplot as plt  # noqa: E402

FS = dc.FS
SEED = 20261007
OUT = HERE / "out"
EVIDENCE_ID = "ev_msr_synth_reverb_decay_predelay_1"


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

def edc_db(h: np.ndarray) -> np.ndarray:
    """Schroeder backward integration: energy remaining after t, in dB re the total."""
    e = np.cumsum((h.astype(np.float64) ** 2)[::-1])[::-1]
    return 10.0 * np.log10(np.maximum(e / e[0], 1e-300))


def decay_time(h: np.ndarray, hi_db: float, lo_db: float, fs: int = FS) -> dict:
    """Least-squares line through the EDC between hi_db and lo_db, extrapolated to -60 dB."""
    edc = edc_db(h)
    idx = np.where((edc <= hi_db) & (edc >= lo_db))[0]
    t = idx / fs
    slope, intercept = np.polyfit(t, edc[idx], 1)  # dB per second
    resid = edc[idx] - (slope * t + intercept)
    return {"seconds": float(-60.0 / slope), "slope_db_per_s": float(slope),
            "fit_start_s": float(t[0]), "fit_end_s": float(t[-1]),
            "fit_rms_error_db": float(np.sqrt(np.mean(resid ** 2)))}


def t20(h):
    return decay_time(h, -5.0, -25.0)


def t30(h):
    return decay_time(h, -5.0, -35.0)


def edt(h):
    return decay_time(h, 0.0, -10.0)


def time_to_db(h: np.ndarray, level_db: float, fs: int = FS) -> float:
    edc = edc_db(h)
    return float(np.argmax(edc <= level_db) / fs)


def octave_band(h: np.ndarray, fc: float, fs: int = FS) -> np.ndarray:
    sos = signal.butter(4, [fc / np.sqrt(2), fc * np.sqrt(2)], btype="bandpass", fs=fs, output="sos")
    return signal.sosfilt(sos, h)


# ---------------------------------------------------------------- generators

def exp_noise_tail(rt60: float, length_s: float, rng: np.random.Generator, predelay_s: float = 0.0) -> np.ndarray:
    """Gaussian noise whose amplitude falls 60 dB in rt60 seconds, starting after predelay_s."""
    n = int(round(length_s * FS))
    t = np.arange(n) / FS
    env = 10.0 ** (-3.0 * t / rt60)  # amplitude: -60 dB at t = rt60
    tail = rng.standard_normal(n) * env
    d = int(round(predelay_s * FS))
    return np.concatenate([np.zeros(d), tail])


def next_prime(n: int) -> int:
    def is_prime(k):
        if k < 2:
            return False
        r = int(k ** 0.5)
        return all(k % p for p in range(2, r + 1))
    while not is_prime(n):
        n += 1
    return n


BASE_DELAYS = np.array([1153, 1327, 1559, 1753, 1999, 2203, 2423, 2687])  # samples at 48 kHz (24-56 ms)


def fdn_delays(size: float) -> np.ndarray:
    return np.array([next_prime(int(round(m * size))) for m in BASE_DELAYS])


def hadamard8() -> np.ndarray:
    h2 = np.array([[1.0, 1.0], [1.0, -1.0]])
    return np.kron(np.kron(h2, h2), h2) / np.sqrt(8.0)


def fdn_ir(delays: np.ndarray, k: np.ndarray, b: np.ndarray, length_s: float) -> np.ndarray:
    """Impulse response of an 8-line FDN (orthogonal Hadamard feedback).

    Each line i: delay m_i, then the absorption filter H_i(z) = k_i (1 - b_i) / (1 - b_i z^-1)
    (b_i = 0 gives a frequency-independent gain k_i). Processed in blocks no longer than the
    shortest delay, so every block reads only samples already written.
    """
    n_total = int(round(length_s * FS))
    N = len(delays)
    A = hadamard8()
    cin = np.ones(N)
    cout = np.array([1, -1, 1, -1, 1, -1, 1, -1], dtype=float) / np.sqrt(N)
    x = np.zeros(n_total)
    x[0] = 1.0
    line_in = np.zeros((N, n_total))
    out = np.zeros(n_total)
    zi = np.zeros((N, 1))
    L = int(delays.min())
    for t0 in range(0, n_total, L):
        t1 = min(t0 + L, n_total)
        y = np.zeros((N, t1 - t0))
        for i, m in enumerate(delays):
            s0, s1 = t0 - m, t1 - m
            if s1 <= 0:
                continue
            a0 = max(s0, 0)
            y[i, a0 - s0:] = line_in[i, a0:s1]
        s = np.empty_like(y)
        for i in range(N):
            s[i], zi[i] = signal.lfilter([k[i] * (1.0 - b[i])], [1.0, -b[i]], y[i], zi=zi[i])
        out[t0:t1] = cout @ s
        line_in[:, t0:t1] = A @ s + np.outer(cin, x[t0:t1])
    return out


def gains_for_t60(delays: np.ndarray, t60: float) -> np.ndarray:
    """Per-line gain so that each pass loses 60 dB per t60 seconds: k = 10^(-3 m / (fs T60))."""
    return 10.0 ** (-3.0 * delays / (FS * t60))


def absorption_for(delays: np.ndarray, t60_dc: float, t60_nyq: float) -> tuple[np.ndarray, np.ndarray]:
    """One-pole absorption filters matched exactly at DC and at Nyquist.

    H(z) = k (1 - b) / (1 - b z^-1) has gain k at DC and k (1 - b) / (1 + b) at Nyquist.
    Set k for T60 at DC and solve (1 - b) / (1 + b) = r for the Nyquist gain:
    b = (1 - r) / (1 + r), r = g_nyq / k.
    """
    k = gains_for_t60(delays, t60_dc)
    g_nyq = gains_for_t60(delays, t60_nyq)
    r = g_nyq / k
    b = (1.0 - r) / (1.0 + r)
    return k, b


def per_line_t60(delays: np.ndarray, k: np.ndarray, b: np.ndarray, freqs: np.ndarray) -> np.ndarray:
    """T60 each line's filter implies at each frequency: -60 m / (fs * 20 log10 |H(f)|). Shape (lines, freqs)."""
    w = 2 * np.pi * freqs / FS
    rows = []
    for m, kk, bb in zip(delays, k, b):
        mag = np.abs(kk * (1 - bb) / (1 - bb * np.exp(-1j * w)))
        rows.append(-60.0 * m / (FS * 20 * np.log10(mag)))
    return np.array(rows)


def predicted_t60_curve(delays: np.ndarray, k: np.ndarray, b: np.ndarray, freqs: np.ndarray) -> np.ndarray:
    """Design T60 per frequency: the mean over the lines (they agree within a few % with this design)."""
    return per_line_t60(delays, k, b, freqs).mean(axis=0)


def predicted_band_t30(fc: float, delays: np.ndarray, k: np.ndarray, b: np.ndarray, length_s: float) -> float:
    """T30 of the decay curve the design predicts inside one octave band.

    Model: flat initial energy spectrum; each frequency decays at its own design T60, so its
    remaining energy at t is proportional to T60(f) 10^(-6 t / T60(f)). Weight by the band
    filter's |H(f)|^2, sum over f, and fit -5..-35 dB like the measurement.
    """
    sos = signal.butter(4, [fc / np.sqrt(2), fc * np.sqrt(2)], btype="bandpass", fs=FS, output="sos")
    f = np.linspace(5.0, FS / 2 - 5.0, 24000)
    _, H = signal.sosfreqz(sos, worN=f, fs=FS)
    w2 = np.abs(H) ** 2
    t60f = predicted_t60_curve(delays, k, b, f)
    t = np.arange(0.0, length_s, 0.001)
    edc = (w2 * t60f)[None, :] * 10.0 ** (-6.0 * t[:, None] / t60f[None, :])
    e = edc.sum(axis=1)
    e_db = 10 * np.log10(e / e[0])
    idx = np.where((e_db <= -5) & (e_db >= -35))[0]
    slope, _ = np.polyfit(t[idx], e_db[idx], 1)
    return float(-60.0 / slope)


# ---------------------------------------------------------------- sound examples

def pluck(length_s: float = 0.5) -> np.ndarray:
    """A synthetic plucked tone: 330 Hz plus two harmonics, 1 ms fade-in, 60 ms decay constant."""
    n = int(round(length_s * FS))
    t = np.arange(n) / FS
    x = (np.sin(2 * np.pi * 330 * t) + 0.5 * np.sin(2 * np.pi * 660 * t) + 0.25 * np.sin(2 * np.pi * 990 * t))
    x *= np.exp(-t / 0.06)
    fade = int(0.001 * FS)
    x[:fade] *= np.linspace(0.0, 1.0, fade)
    return x / np.max(np.abs(x))


def render(dry: np.ndarray, ir: np.ndarray, wet_db: float, total_s: float) -> np.ndarray:
    n = int(round(total_s * FS))
    ir = ir / np.sqrt(np.sum(ir ** 2))
    wet = signal.fftconvolve(dry, ir)[:n]
    y = np.zeros(n)
    y[: len(dry)] += dry
    y[: len(wet)] += 10 ** (wet_db / 20) * wet
    return y


# ---------------------------------------------------------------- main

def main() -> None:
    dc.clear_outputs(HERE, OUT)
    dc.setup_style()
    rng = np.random.default_rng(SEED)
    results: dict = {}

    # D1: recovering a known RT60
    d1 = []
    edcs = {}
    for rt in (0.5, 1.0, 2.0, 4.0):
        h = exp_noise_tail(rt, length_s=rt * 1.6, rng=rng)
        row = {"set_rt60_s": rt, "t20": t20(h), "t30": t30(h), "edt": edt(h)}
        d1.append(row)
        edcs[rt] = edc_db(h)
    devs = [abs(r[k]["seconds"] - r["set_rt60_s"]) / r["set_rt60_s"] * 100 for r in d1 for k in ("t20", "t30", "edt")]
    results["D1_known_rt60"] = {"rows": d1, "max_abs_deviation_percent": max(devs)}

    fig, ax = plt.subplots(figsize=(7.2, 3.8))
    for i, rt in enumerate((0.5, 1.0, 2.0, 4.0)):
        e = edcs[rt]
        t = np.arange(len(e)) / FS
        keep = e > -80
        ax.plot(t[keep], e[keep], color=dc.SERIES[i], label=f"RT60 set {rt:g} s")
        ax.annotate(f"{rt:g} s", xy=(rt, -60), xytext=(rt, -66), ha="center", fontsize=8, color=dc.TEXT_2)
    ax.axhline(-60, color=dc.TEXT_2, lw=0.8, ls="--")
    ax.set_xlim(0, 4.6)
    ax.set_ylim(-80, 2)
    ax.set_xlabel("time (s)")
    ax.set_ylabel("energy remaining (dB)")
    ax.set_title("D1  Schroeder decay curves of tails with a set RT60")
    ax.legend(loc="upper right")
    fig.tight_layout()
    dc.save_figure(fig, OUT / "d1_edc_known_rt60.png")
    dc.write_curve_csv(OUT / "d1_known_rt60.csv", "curve_csv",
                       ["設定した残響時間 [s]", "T20 [s]", "T30 [s]", "EDT [s]"],
                       [np.array([r["set_rt60_s"] for r in d1])] + [np.array([r[k]["seconds"] for r in d1]) for k in ("t20", "t30", "edt")],
                       [3, 6, 6, 6])

    # D2: pre-delay
    d2 = []
    base = exp_noise_tail(1.5, length_s=2.6, rng=rng)
    edc_pd = {}
    for pd_ms in (0, 20, 80):
        h = np.concatenate([np.zeros(int(round(pd_ms * FS / 1000))), base])
        r = t30(h)
        d2.append({"predelay_ms": pd_ms, "t30_s": r["seconds"], "slope_db_per_s": r["slope_db_per_s"],
                   "time_from_start_to_minus30db_s": time_to_db(h, -30.0),
                   "first_nonzero_sample_s": float(np.argmax(np.abs(h) > 0) / FS)})
        edc_pd[pd_ms] = edc_db(h)
    results["D2_predelay"] = {"set_rt60_s": 1.5, "rows": d2}

    fig, ax = plt.subplots(figsize=(7.2, 3.6))
    for i, pd_ms in enumerate((0, 20, 80)):
        e = edc_pd[pd_ms]
        t = np.arange(len(e)) / FS
        keep = (e > -50) & (t < 0.6)
        ax.plot(t[keep] * 1000, e[keep], color=dc.SERIES[i], label=f"pre-delay {pd_ms} ms")
    ax.set_xlim(0, 600)
    ax.set_ylim(-25, 2)
    ax.set_xlabel("time from the dry signal (ms)")
    ax.set_ylabel("energy remaining (dB)")
    ax.set_title("D2  Same tail (RT60 1.5 s), three pre-delays: shifted, same slope")
    ax.legend(loc="upper right")
    fig.tight_layout()
    dc.save_figure(fig, OUT / "d2_predelay_edc.png")
    t_ms = np.arange(0, 601, 1.0)
    dc.write_curve_csv(OUT / "d2_predelay_edc.csv", "decay_curve_csv",
                       ["時間 [ms]", "プリディレイ0ms [dB]", "プリディレイ20ms [dB]", "プリディレイ80ms [dB]"],
                       [t_ms] + [edc_pd[p_][(t_ms * FS / 1000).astype(int)] for p_ in (0, 20, 80)], [0, 3, 3, 3])
    dc.write_curve_csv(OUT / "d2_predelay_t30.csv", "curve_csv",
                       ["プリディレイ [ms]", "T30 [s]", "-30dBに達する時刻 [s]"],
                       [np.array([r["predelay_ms"] for r in d2]), np.array([r["t30_s"] for r in d2]),
                        np.array([r["time_from_start_to_minus30db_s"] for r in d2])], [0, 6, 6])

    # D3: decay as a time (mode A) versus decay as a feedback amount (mode B)
    d3 = []
    sizes = (0.5, 1.0, 2.0)
    g_fixed = 0.93
    target = 2.0
    for size in sizes:
        m = fdn_delays(size)
        hA = fdn_ir(m, gains_for_t60(m, target), np.zeros(8), length_s=target * 1.6 + 0.3)
        hB = fdn_ir(m, np.full(8, g_fixed), np.zeros(8), length_s=12.0 * size + 0.3)
        d3.append({"size": size, "delays_samples": m.tolist(), "mean_delay_ms": float(m.mean() / FS * 1000),
                   "mode_A_target_t60_s": target, "mode_A_t30": t30(hA),
                   "mode_B_gain": g_fixed, "mode_B_t30": t30(hB)})
    for r in d3:
        r["mode_A_t30_ratio_to_size1"] = r["mode_A_t30"]["seconds"] / d3[1]["mode_A_t30"]["seconds"]
        r["mode_B_t30_ratio_to_size1"] = r["mode_B_t30"]["seconds"] / d3[1]["mode_B_t30"]["seconds"]
    results["D3_decay_semantics"] = d3

    fig, ax = plt.subplots(figsize=(6.4, 3.6))
    xs = np.array(sizes)
    ya = np.array([r["mode_A_t30"]["seconds"] for r in d3])
    yb = np.array([r["mode_B_t30"]["seconds"] for r in d3])
    ax.plot(xs, ya, color=dc.SERIES[0], marker="o", ms=7, label="A: decay set as a time (T60 2 s)")
    ax.plot(xs, yb, color=dc.SERIES[1], marker="o", ms=7, label=f"B: decay set as feedback gain {g_fixed}")
    for x, y in zip(xs, ya):
        ax.annotate(f"{y:.2f} s", (x, y), textcoords="offset points", xytext=(0, 8), ha="center", fontsize=8, color=dc.TEXT_2)
    for x, y in zip(xs, yb):
        ax.annotate(f"{y:.2f} s", (x, y), textcoords="offset points", xytext=(0, -14), ha="center", fontsize=8, color=dc.TEXT_2)
    ax.set_xscale("log", base=2)
    ax.set_xticks(xs, [f"x{s:g}" for s in sizes])
    ax.set_ylim(0, max(yb.max(), ya.max()) * 1.25)
    ax.set_xlabel("size (all delay lengths scaled)")
    ax.set_ylabel("measured T30 (s)")
    ax.set_title("D3  Toy FDN: what 'decay' means changes what size does")
    ax.legend(loc="upper left")
    fig.tight_layout()
    dc.save_figure(fig, OUT / "d3_size_vs_decay.png")
    dc.write_curve_csv(OUT / "d3_size_vs_decay.csv", "curve_csv",
                       ["大きさ [倍]", "A 時間として扱う [s]", "B 帰還の量として扱う [s]"], [xs, ya, yb], [2, 6, 6])

    # D4: frequency-dependent decay
    m = fdn_delays(1.0)
    t60_dc, t60_nyq = 2.5, 0.2
    k, b = absorption_for(m, t60_dc, t60_nyq)
    h_damped = fdn_ir(m, k, b, length_s=t60_dc * 1.6 + 0.3)
    h_flat = fdn_ir(m, gains_for_t60(m, t60_dc), np.zeros(8), length_s=t60_dc * 1.6 + 0.3)
    bands = np.array([125, 250, 500, 1000, 2000, 4000, 8000], dtype=float)
    pred = predicted_t60_curve(m, k, b, bands)
    d4 = []
    for fc, pr in zip(bands, pred):
        d4.append({"band_hz": float(fc), "predicted_t60_s": float(pr),
                   "predicted_band_t30_s": predicted_band_t30(fc, m, k, b, t60_dc * 1.6),
                   "damped_t30": t30(octave_band(h_damped, fc)),
                   "flat_t30": t30(octave_band(h_flat, fc))})
    spread = per_line_t60(m, k, b, bands)
    for r, lo, hi in zip(d4, spread.min(axis=0), spread.max(axis=0)):
        r["per_line_design_t60_min_s"], r["per_line_design_t60_max_s"] = float(lo), float(hi)
    results["D4_frequency_dependent"] = {"design_t60_dc_s": t60_dc, "design_t60_nyquist_s": t60_nyq,
                                         "broadband_t30_damped": t30(h_damped), "broadband_t30_flat": t30(h_flat),
                                         "rows": d4}

    fig, ax = plt.subplots(figsize=(6.8, 3.8))
    fgrid = np.geomspace(60, 16000, 200)
    ax.plot(fgrid, predicted_t60_curve(m, k, b, fgrid), color=dc.SERIES[1], lw=1.5, ls="--", label="damped: design T60 per frequency")
    ax.plot(bands, [r["predicted_band_t30_s"] for r in d4], color=dc.SERIES[1], marker="s", ms=8, lw=0, mfc="none", mew=1.5, label="damped: design, octave-band T30")
    ax.plot(bands, [r["damped_t30"]["seconds"] for r in d4], color=dc.SERIES[1], marker="o", ms=5, lw=0, label="damped: measured octave-band T30")
    ax.plot(bands, [r["flat_t30"]["seconds"] for r in d4], color=dc.SERIES[0], marker="o", ms=7, label="no damping: measured octave-band T30")
    ax.set_xscale("log")
    ax.set_xticks(bands, ["125", "250", "500", "1k", "2k", "4k", "8k"])
    ax.minorticks_off()
    ax.set_ylim(0, 3.0)
    ax.set_xlabel("octave band centre (Hz)")
    ax.set_ylabel("reverberation time (s)")
    ax.set_title("D4  Same network, T60 2.5 s at low frequency: with and without HF damping")
    ax.legend(loc="lower left")
    fig.tight_layout()
    dc.save_figure(fig, OUT / "d4_band_rt.png")
    dc.write_curve_csv(OUT / "d4_band_rt.csv", "curve_csv",
                       ["帯域の中心 [Hz]", "ダンピングなし [s]", "ダンピングあり [s]", "ダンピングあり 設計の帯域平均 [s]"],
                       [bands, np.array([r["flat_t30"]["seconds"] for r in d4]), np.array([r["damped_t30"]["seconds"] for r in d4]),
                        np.array([r["predicted_band_t30_s"] for r in d4])], [0, 6, 6, 6])

    # D5: pre-delay before everything vs before the late tail only
    er_times_ms = np.array([7.0, 11.0, 17.0, 23.0, 31.0, 41.0])
    er_gains = np.array([0.8, -0.65, 0.55, -0.45, 0.4, -0.3])
    late_start_ms = 45.0
    late = exp_noise_tail(1.2, length_s=2.0, rng=rng) * 0.12
    pd_ms = 30.0

    def build(er_shift_ms: float, late_shift_ms: float) -> np.ndarray:
        n = int(round((late_start_ms + late_shift_ms) * FS / 1000)) + len(late) + FS // 10
        h = np.zeros(n)
        for tm, g in zip(er_times_ms + er_shift_ms, er_gains):
            h[int(round(tm * FS / 1000))] += g
        s = int(round((late_start_ms + late_shift_ms) * FS / 1000))
        h[s:s + len(late)] += late
        return h

    h_p0 = build(0.0, 0.0)
    h_p1 = build(pd_ms, pd_ms)   # pre-delay before the early reflections and the tail
    h_p2 = build(0.0, pd_ms)     # pre-delay before the late tail only
    d5 = {"predelay_ms": pd_ms, "early_reflection_times_ms": er_times_ms.tolist(), "late_start_ms_without_predelay": late_start_ms}
    for name, h in (("none", h_p0), ("before_all", h_p1), ("before_late_only", h_p2)):
        d5[name] = {"first_reflection_ms": float(np.argmax(np.abs(h) > 0) / FS * 1000),
                    "late_start_ms": float((late_start_ms + (0 if name == "none" else pd_ms))),
                    "t30": t30(h)}
    results["D5_predelay_placement"] = d5

    fig, axes = plt.subplots(3, 1, figsize=(7.2, 5.2), sharex=True)
    for ax, (label, h), col in zip(axes, (("no pre-delay", h_p0), ("30 ms before early reflections and tail", h_p1),
                                          ("30 ms before the late tail only", h_p2)), (dc.SERIES[0], dc.SERIES[1], dc.SERIES[2])):
        t_ms = np.arange(len(h)) / FS * 1000
        keep = t_ms < 160
        ax.vlines(t_ms[keep][h[keep] != 0], 0, np.abs(h[keep][h[keep] != 0]), color=col, lw=1.0)
        ax.set_ylim(0, 0.9)
        ax.set_yticks([])
        ax.text(0.99, 0.85, label, transform=ax.transAxes, ha="right", va="top", fontsize=9, color=dc.TEXT)
    axes[-1].set_xlabel("time from the dry signal (ms)")
    axes[0].set_title("D5  Where the pre-delay goes (|amplitude| of the first 160 ms)")
    fig.tight_layout()
    dc.save_figure(fig, OUT / "d5_predelay_placement.png")
    bin_n = int(0.0005 * FS)                      # 0.5 ms bins, peak |amplitude| in each
    n_bins = int(0.160 * FS) // bin_n

    def binned(h):
        return np.abs(h[: n_bins * bin_n]).reshape(n_bins, bin_n).max(axis=1)

    dc.write_curve_csv(OUT / "d5_predelay_placement.csv", "curve_csv",
                       ["時間 [ms]", "プリディレイなし [-]", "全体の前に30ms [-]", "後部残響の前だけに30ms [-]"],
                       [np.arange(n_bins) * 0.5, binned(h_p0), binned(h_p1), binned(h_p2)], [1, 6, 6, 6])

    # Sound examples (dry pluck plus wet at -12 dB re the dry; one common gain for all files)
    dry = pluck()
    m1 = fdn_delays(1.0)
    ir_1s = fdn_ir(m1, gains_for_t60(m1, 1.0), np.zeros(8), length_s=2.0)
    ir_2p5_flat = h_flat
    ir_2p5_damped = h_damped
    total = 4.0
    clips = {
        "a1_dry.wav": np.pad(dry, (0, int(total * FS) - len(dry))),
        "a2_t60_1s_predelay_0ms.wav": render(dry, ir_1s, -12.0, total),
        "a3_t60_1s_predelay_80ms.wav": render(dry, np.concatenate([np.zeros(int(0.08 * FS)), ir_1s]), -12.0, total),
        "a4_t60_2p5s_no_damping.wav": render(dry, ir_2p5_flat, -12.0, total),
        "a5_t60_2p5s_hf_damping.wav": render(dry, ir_2p5_damped, -12.0, total),
    }
    peak = max(np.max(np.abs(x)) for x in clips.values())
    common_gain = 10 ** (-1.0 / 20) / peak
    for name, x in clips.items():
        dc.write_wav(OUT / name, x * common_gain)
    results["sound_examples"] = {"source": "synthetic pluck 330/660/990 Hz, 60 ms decay constant",
                                 "mix": "impulse response normalised to unit energy, convolved, added to the dry at -12 dB",
                                 "wet_gain_db": -12.0, "common_gain_db": float(20 * np.log10(common_gain)),
                                 "files": list(clips), "committed": False}

    detail = {
        "evidence_id": EVIDENCE_ID,
        "article_slug": "reverb-time-decay-predelay",
        "signal_seed": SEED,
        "sample_rate_hz": FS,
        "common_module_sha256": dc.sha256_file(Path(dc.__file__)),
        "figures_module_sha256": dc.sha256_file(HERE / "figures.py"),
        "svg_module_sha256": dc.sha256_file(HERE.parent / "common" / "svgkit.py"),
        "method": {
            "decay_curve": "Schroeder backward integration of the squared impulse response, dB re total energy",
            "t20": "least-squares line on the decay curve from -5 to -25 dB, extrapolated to -60 dB",
            "t30": "least-squares line on the decay curve from -5 to -35 dB, extrapolated to -60 dB",
            "edt": "least-squares line on the decay curve from 0 to -10 dB, extrapolated to -60 dB",
            "predelay": "pure delay: zeros prepended to the response",
            "octave_bands": "4th-order Butterworth band-pass (scipy.signal.butter, sosfilt), edges fc/sqrt2 to fc*sqrt2",
            "fdn": "8 delay lines, Hadamard/sqrt(8) feedback, per-line gain k=10^(-3m/(fs T60)); damping: one-pole k(1-b)/(1-b z^-1) matched to T60 at DC and at Nyquist, b=(1-r)/(1+r), r=g_nyq/k",
        },
        "results": results,
    }
    ref = next(r for r in d1 if r["set_rt60_s"] == 2.0)
    ledger = {
        "summary_editorial": "合成した応答で、残響時間の測り方と、プリディレイ、帰還の量、ダンピングの影響を確かめた実証の記録です。",
        "subject_ja": "残響時間を後方積分で測る方法と、プリディレイ・帰還の量・高域のダンピングが残響時間に与える影響を、合成した応答で確かめる実証です。",
        "material_description_ja": "スクリプトが乱数と式から生成したインパルス応答（指数減衰の雑音と、8本の遅延線の帰還遅延網）と、合成した撥弦の音です。",
        "settings": [
            {"name": "標本化周波数", "value": "48000 Hz"},
            {"name": "残響時間の計算", "value": "後方積分の曲線の-5〜-35 dBに直線を当て、-60 dBまで延ばす（T30）"},
            {"name": "結果の基準の尾", "value": "指数減衰の雑音、設定した残響時間 2.0 s"},
            {"name": "プリディレイ", "value": "0 ms、20 ms、80 ms（純粋な遅延）"},
            {"name": "帰還遅延網", "value": "遅延線8本、Hadamard 行列、1倍のとき 1153〜2687 サンプル"},
            {"name": "帰還の量（モードB）", "value": "0.93"},
            {"name": "ダンピング", "value": "1次フィルター、0 Hzで2.5 s、ナイキスト周波数で0.2 s"},
            {"name": "帯域", "value": "1オクターブ、4次のバターワース帯域通過"},
        ],
        "comparison_conditions_ja": ["同じ seed で生成した同じ尾と同じ遅延網を、1つの条件だけ変えて比べました。"],
        "procedure_ja": [
            "乱数の seed を固定し、残響時間を決めた指数減衰の雑音と、8本の遅延線の帰還遅延網のインパルス応答を生成しました。",
            "応答を2乗して後ろから積分し、減衰の曲線をdBで求めました。",
            "曲線の-5 dBから-35 dBの区間に直線を当てはめ、60 dB下がる時間に延ばした値をT30としました。",
            "同じ尾に0 ms、20 ms、80 msのプリディレイを付け、T30と-30 dBに達する時刻を比べました。",
            "帰還の量を目標の残響時間から計算する場合と0.93に固定する場合で、遅延線の長さを変えてT30を比べました。",
            "遅延線ごとに1次のダンピングフィルターを入れ、1オクターブの帯域ごとにT30を求めました。",
            "初期反射と後部残響を分けた応答で、プリディレイを全体の前に置く場合と後部残響の前だけに置く場合を比べました。",
        ],
        "results": {"rt60_seconds": {"value": round(ref["t30"]["seconds"], 3), "decimals": 3}},
        "signal_seed": SEED,
    }
    csv_kinds = {"d1_known_rt60.csv": "curve_csv", "d2_predelay_edc.csv": "decay_curve_csv", "d2_predelay_t30.csv": "curve_csv",
                 "d3_size_vs_decay.csv": "curve_csv", "d4_band_rt.csv": "curve_csv", "d5_predelay_placement.csv": "curve_csv"}
    figures.write_all(OUT)
    rec = dc.write_outputs(HERE, OUT, EVIDENCE_ID, csv_kinds, detail, ledger)
    print(rec)

if __name__ == "__main__":
    main()