diff --git a/HANDOFF.md b/HANDOFF.md index 2276789..b6b730b 100644 --- a/HANDOFF.md +++ b/HANDOFF.md @@ -344,6 +344,10 @@ carries the measured numbers; this is the map: reference lags one block, so the cabin S1 held note reads +9.71 against the karaoke suite's +11.63. 11. **Deferred to phase 2's reverb-integration item: `mutap::spectral_reverb`.** + **Its return is the next PR (phase 2 item 7, part 2)**, on the harness + part 1 built (item 13: `reverb_mix`, the loop's forward stage, the + probe-convergence sweep, the audible driver; `docs/reverb-afc.md`), so + the spectral reverb is measured exactly as the Dattorro plates were. The per-bin spectral reverb shaped from F̂ is commit `8abe958` on the local branch `karaoke-afc-decorrelation-reverb` (from the karaoke bundle); it did not land. Measured during the karaoke landing with a @@ -428,6 +432,66 @@ carries the measured numbers; this is the map: the canceller's limit + 6 instead. - The latch never engaged in a gated row; it is covered by the state machine suite only. +13. **The reverb behind the canceller, part 1: the Dattorro plates + (2026-10-02; [`docs/reverb-afc.md`](docs/reverb-afc.md)).** Phase 2 + item 7, part 1; the spectral reverb's return is part 2 (item 11). + Landed: `include/mutap/reverb_stage.h` (`tap::mu::reverb_mix`, y = (1 − w)·x + w·r with r the reverb's L or (L + R)/2, and + `shifted_dry_mix`, the dry-only-shift + topology; any reverb with the FAUST shim's `process(in, out, n)`; + header-only, no FFT, no ABI tag, `latency()` 0); the float plumbing + suite `test_reverb_mix.cpp` in both emulated selections (7 tests; the + bare-metal count comment is now 108); `decorrelated_loop`'s forward + stage (reset per probe; bit-identical without one, checked on the + karaoke measurement against origin/main); `tests/support/reverb_rig.h`; + the gated host rows `test_reverb_stage.cpp`; the `MUTAP_SLOW` sweep + `test_reverb_stage_sweep.cpp` (ProbeConvergence, Grid, Shift, + BareLoopCost); `karaoke_ramp_dump --reverb ...` and + `tools/notebook/reverb_audible.py`. Measured (macOS x86_64, AppleClang, + Release, five seed sets): + - **Probe convergence:** every held-note row (both plates, decay 0.5 to + 0.85, the 5 Hz shift rows; cabin and mt5, S1) converged at 10 s, its + 20 to 160 s medians identical; phase 0's 384 s was the dry loop with + the plate and does not carry over to the canceller-held loop. The + speech envelope drifts down with probe length; 3 of its 8 rows had + not converged by their longest probe (cabin canceller to 80 s, cabin + decay 0.5 and mt5 decay 0.7 to 160 s); none is gated. + - **Runaway cost behind the canceller (Grid, six rooms, held note, 20 s + probes):** the as-shipped plate at decay 0.5 / 0.7 / 0.85 × wet 0.15 / + 0.30 / 0.50 costs at most +1.46 dB (median of room medians, S1; +2.69 + at matched loudness) and nothing at S3; decay 0.85 at w 0.50 RAISES + the S1 limit by 4.39 dB. The bare loop charges the same plate up to + +5.79 dB (bound). Damping 0.5, the (L + R)/2 return and the paper + plate (decay 0.5) move the cost by at most 1.07 dB, inside the rows' + room spreads. Speech envelope (S1, 80 s, wet 0.30): −1.17 / −0.88 / + −0.29 dB raw. + - **Audible (cabin, ramp + criterion, PROTOCOL §7.3 flags):** the + plate lowers the audible limit (S1: canceller +16.57; plate 0.5 / 0.7 + at w 0.30 +12.89 / +10.46) while the ramp's runaway does not fall. + With a 2 or 5 Hz shift the dry-only topology (`shifted_dry_mix`) is + audible later than the whole-bus shift in 13 of 16 pairs; bisected + runaway says the opposite at S1 (whole bus higher in 22 of 24 room + medians; gated in the cabin, +3.09 dB) and at S3 at 2 Hz; at S3, 5 Hz + the two do not separate. + - **dattorro_paper bare-loop table** (phase 0's method; reproduces the + as-shipped table to the digit): at equal decay it costs more than as + shipped, at equal T30 less (paper 0.5 +0.30 against as-shipped 0.7 + +1.11 dB bound at w 0.30). T30 at decay 0.3 / 0.5 / 0.7 / 0.85: 1.7252 + / 2.1186 / 3.7882 / 8.2457 s (L). + Gated (`ReverbStage.*`, host-only): the mix over both plates, the T30s + and mix levels, three bare-loop cells, the cabin S1 cost ceiling (+3.5 + dB on medians measured ≤ +0.56) and the topology direction. Every + number is one host's; no second host ran the sweep. Open: + - **The S3 reverb-only rows at wet 0.30 are audible below the dry + loop's limit in some seeds** (−13.55 to +1.89 at decay 0.5); what the + criterion flags there was not examined. + - **Runaway and audible disagree on the shift topology**; phase 1's + rooms decide. + - **Not run:** speech at S3, the paper plate in the loop at decay 0.7 / + 0.85, bisected shift rows at wet 0.15 / 0.50, audible rows beyond the + cabin held note, float32 in the loop, the backing track, two mics. + - The audible driver's first invocation filled the scratch disk (about + 180 MB of WAV a run); it now deletes WAVs and checkpoints per run. ## The next effort (Rev 4): AEC objects + echo chapter diff --git a/README.md b/README.md index 2bca29b..15fbcc6 100644 --- a/README.md +++ b/README.md @@ -299,6 +299,14 @@ algorithmically complete. What exists today: `howl_detector` on each residual, attached to `afc_chain` with `set_guard()`, with a soundcheck sampler that sets each mic's verdict thresholds; measured live in [`docs/howl-guard.md`](docs/howl-guard.md). +- **The reverb stage** — `tap::mu::reverb_mix` and `tap::mu::shifted_dry_mix` + ([`include/mutap/reverb_stage.h`](include/mutap/reverb_stage.h)): any + mono-in reverb engine in `afc_chain`'s reverb slot as + y = (1 − w)·x + w·r, and the topology that frequency-shifts the dry path + only so a reverb tail never recirculates through the shifter; + header-only, no FFT. The vendored FAUST Dattorro plates measured behind + the canceller, with converged probes, in + [`docs/reverb-afc.md`](docs/reverb-afc.md). Next up (see [HANDOFF.md](HANDOFF.md) "What's next"): in-Max listening in a real room and the default-engine decision, then the M55 performance diff --git a/docs/karaoke-afc.md b/docs/karaoke-afc.md index 111d5c9..a192336 100644 --- a/docs/karaoke-afc.md +++ b/docs/karaoke-afc.md @@ -314,7 +314,10 @@ does not land with this work, for three reasons: - The class also needs the FFT ABI tag. It returns with phase 2's reverb-integration item. HANDOFF records what was -measured and the code's location. +measured and the code's location. That item's first part put the FAUST +Dattorro plates behind the canceller and measured them with converged +probes ([A reverb behind the canceller](reverb-afc.md)); the spectral +reverb comes back on the same harness in the next part. ## Reading these numbers diff --git a/docs/reverb-afc.md b/docs/reverb-afc.md new file mode 100644 index 0000000..b1985bf --- /dev/null +++ b/docs/reverb-afc.md @@ -0,0 +1,599 @@ +# A reverb behind the canceller + +*Measured October 2026 in simulation: the cabin RIR fixture and five other +rooms from both generator families, band-limited through a loudspeaker +model, and phase 0's four room fixtures at full length for the bare loop. +Host: macOS 15.7 x86_64 (i9-8950HK), AppleClang 17, Release, on a shared +machine ([Hosts](#hosts)). Nothing here has been heard in a real room. Every +number below comes from `tests/test_reverb_stage.cpp`, its `MUTAP_SLOW` +sweep `tests/test_reverb_stage_sweep.cpp`, or the +`tools/notebook/reverb_audible.py` driver; [Provenance](#provenance) gives +the commands.* + +A reverb in the voice bus sits inside the feedback loop: whatever it does +to the bus, it does to the howl. This document puts FAUST's Dattorro plate, +as shipped and with the paper's delay lengths, in `afc_chain`'s reverb slot +behind the PEM + FD-Kalman canceller, and measures what it costs, with the +probe lengths shown to converge first. It is the first half of the reverb +work; the spectral reverb shaped from the canceller's estimate comes back +on the same harness next (HANDOFF item 11). + +## Headline + +- **Behind the canceller the plate costs almost no runaway gain.** Held + note, six rooms, the as-shipped plate at decay 0.5 to 0.85 and wet 0.15 + to 0.50: the largest median cost is +1.46 dB at S1 and none at S3 (at + most +2.69 dB at matched loudness). Long decays at high wet even raise the + runaway limit (S1, decay 0.85, w 0.50: −4.39 dB of cost). The bare loop, + with no canceller, charges that same setting +5.79 dB by the magnitude + bound (+6.25 exact, four fixtures at d = 480). +- **The audible limit is where the plate costs.** In the cabin at S1 the + canceller alone is audible at +16.57 dB; with the plate at decay 0.5 / + 0.7 and wet 0.30, at +12.89 / +10.46, while the ramp's runaway does not + move (+20.70 alone, +22.07 / +21.76 with the plate). At S3, wet 0.30, + some seeds are audible below the dry loop's limit (not examined). +- **With a frequency shift, shifting the dry path only is audible later** + than shifting the whole bus ahead of the plate: per-seed medians favour + it in 13 of 16 decay × wet × shift × delay pairs, by up to +5.11 dB + (5 Hz, decay 0.7, w 0.30, S3). Bisected runaway says the opposite (the + whole-bus shift holds more in 22 of 24 room medians at S1; at S3 at 2 Hz + only), so the topology is a listening and room question, not settled + here. +- **Probes: the held note converges at 10 s, decay 0.85 included**; the + speech envelope drifts and is reported, not gated. +- **The paper plate**, new in the bare-loop table: at equal decay it costs + more than the plate as shipped (it rings longer); at equal T30 (paper 0.5 + against as-shipped 0.7, 2.12 s) it costs less by the bound, +0.30 against + +1.11 dB at w 0.30. Behind the canceller, at decay 0.5, the two do not + separate. + +Every runaway and audible number is the held note's or the speech +envelope's on simulated rooms, double precision, on one host; read the +held-note shift rows as the karaoke document does, as a worst case for the +shift. + +## The stage + +`include/mutap/reverb_stage.h` puts any mono-in reverb engine in +`afc_chain`'s reverb slot (decorrelator → reverb → safety): + +- **`tap::mu::reverb_mix`**: y = (1 − w)·x + w·r(x), where + r is the reverb's mono return, output L or (L + R)/2 + (`reverb_return::left` / `mid`), and w the wet level in [0, 1]. +- **`tap::mu::shifted_dry_mix`**: y = (1 − w)·s(x) + + w·r(x). The shifter sits on the dry path inside the stage and the reverb + is fed the unshifted bus, so the reverb's tail never passes through the + shifter. It goes in the reverb slot with the decorrelator slot empty. + +`Reverb` is any type with `process(const Sample* const* in, Sample* const* +out, int n) noexcept`, one input and N ≥ 1 outputs: the shape of the FAUST +shim's `mutap_faust::faust_block`. The FAUST-generated plates stay in +`third_party/faust/`; the header knows only that signature. The stages do +not own the reverb or the shifter. Both are header-only, hold no FFT (no ABI +tag), report `latency()` 0, allocate their return buffers (channels × block) +at construction, and are `noexcept` and allocation-free afterwards, float32 +included. `set_wet()` and `set_return_mode()` change the mix between blocks. +`tests/test_reverb_mix.cpp` checks the plumbing on every target, both +emulated selections included, with a stand-in reverb. + +The mix is a crossfade, so it lowers the bus level. For white input the +level change is 10·log10 of the energy of (1 − w)·δ + w·r. For the vendored +plates (L return, damping 0.0005; `ReverbStage.PlateDecayAndMixLevel`): + +| decay | w 0.15 | w 0.30 | w 0.50 | +|---|---|---|---| +| 0.3 | −0.963 | −1.721 | −2.135 | +| 0.5 | −0.956 | −1.687 | −2.033 | +| 0.7 | −0.925 | −1.546 | −1.622 | +| 0.85 | −0.833 | −1.134 | −0.547 | + +Both plates give the same values to the printed digit. A loop measured at +equal forward gain with and without the mix therefore includes this level +change; the tables below give the cost both raw and at matched loudness +(raw cost minus the level change). + +## The two plates + +Both are vendored in `third_party/faust/` (its README has the details): +`dattorro` is FAUST's `re.dattorro_rev` as shipped (the source of the +`dattorro~` Max external), and `dattorro_paper` uses the paper's delay +lengths with the tank bug fixed. At equal `decay` the paper plate rings +longer. T30 by Schroeder integration, −5 to −35 dB, L / (L + R)/2, damping +0.0005 (`ReverbStage.PlateDecayAndMixLevel`): + +| decay | as shipped | paper lengths | +|---|---|---| +| 0.3 | 1.0144 / 0.8620 s | 1.7252 / 1.5975 s | +| 0.5 | 1.1846 / 1.3082 s | 2.1186 / 2.3626 s | +| 0.7 | 2.1173 / 2.1530 s | 3.7882 / 3.8598 s | +| 0.85 | 4.6220 / 4.5825 s | 8.2457 / 8.1939 s | + +The paper plate at decay 0.5 has the as-shipped plate's decay-0.7 T30 (L +2.1186 against 2.1173 s). + +## How the loop is measured + +The karaoke suites' loop and protocol ([karaoke-afc.md](karaoke-afc.md), +"How it is measured"), with the reverb in the chain's slot: + +- **Loop.** `tests/support/decorrelated_loop.h`: mic → canceller → forward + delay → stage → gain → speaker → room. The stage is + `decorrelated_loop`'s new forward stage, run on each block after the + forward delay and before the gain: the chain's reverb slot. Every loop + built with a stage resets it, so each bisection probe starts the plate + from silence, as it starts the loop's own history from zero. + `tests/support/reverb_rig.h` builds the stage: `reverb_mix` over a + vendored plate, the library's `frequency_shifter` then `reverb_mix` (the + whole-bus shift), `shifted_dry_mix` (the dry-only shift), or the shifter + alone. Every shift row in this document uses the library's shifter (its + recursive oscillator); the karaoke suites' shift rows use the loop's + exact-ramp oscillator, with the same Hilbert pair. +- **Rooms.** Six band-limited 1024-tap rooms: the image-source fixtures + cabin, studio, rehearsal and hall, and the `random_decaying_rir` rooms mt5 + and mt9 (both generator families). S1 is a 480-sample forward delay + (10 ms), S3 960 (20 ms). +- **Materials.** The held 300 Hz note (`voiced_near_end`, 160-sample + period) and the speech-envelope material (`ar_near_end`). +- **Bisection.** `tests/support/karaoke_asg.h`: PEM + FD-Kalman (speech + cascade, default configuration) converges 1500 blocks at `exact_msg_db` + − 6 dB **with the reverb in the loop**; then the forward gain is bisected + to 0.1 dB over [exact − 20, exact + 30] with the 40 dB runaway rule on a + second seed of the material, a fresh copy of the converged canceller and + a silent plate per probe. Five seed sets (2, 22, 42, 62, 82). +- **Like for like.** + - ASG = the chain's limit − the dry open loop's limit (no canceller, no + reverb), bisected on the same probe material and length. + - Cost = the canceller-alone limit − the canceller-plus-reverb limit, per + seed; positive means the reverb costs stable gain. + - Matched = cost − the mix's level change (the table in + [The stage](#the-stage)): the cost had the mix been made as loud as the + dry bus. The mix is quieter, so matched is the larger. +- **Medians.** Per room, the median of the five seed sets; across rooms, + the median of the six per-room medians, with the per-room range and the + per-seed range over all 30 runs beside it. The loop is chaotic, so single + seeds move by dB; read directions, not single magnitudes. + +## Probe convergence first + +Phase 0 needed 384 s probes before a dry loop with the plate in it +converged (4.6 s tail), and on the karaoke branch the spectral reverb's +in-loop bisection was still moving at 40 s. So before any in-loop number, +`ReverbStageSweep.ProbeConvergence` bisects the cabin and mt5 at S1, wet +0.30, with probes of 10, 20, 40, 80 and 160 s. **Converged** means the +shortest probe from which every longer probe's five-seed median ASG lies +within 0.25 dB of the row's longest probe. Median ASG, dB: + +| Room, material | Row (wet 0.30, damping 0.0005, L) | 10 s | 20 s | 40 s | 80 s | 160 s | Converged | +|---|---|---|---|---|---|---|---| +| cabin, held | canceller alone | +11.66 | +11.74 | +11.74 | +11.74 | | 10 s | +| cabin, held | as shipped, decay 0.5 | +10.98 | +11.05 | +11.05 | +11.05 | +11.05 | 10 s | +| cabin, held | as shipped, decay 0.7 | +10.98 | +11.05 | +11.05 | +11.05 | +11.05 | 10 s | +| cabin, held | as shipped, decay 0.85 | +13.71 | +13.79 | +13.79 | +13.79 | +13.79 | 10 s | +| cabin, held | paper, decay 0.5 | +11.17 | +11.25 | +11.25 | +11.25 | +11.25 | 10 s | +| cabin, held | decay 0.7 + 5 Hz, whole bus | +18.59 | +18.67 | +18.67 | +18.67 | | 10 s | +| cabin, held | decay 0.7 + 5 Hz, dry only | +17.03 | +17.11 | +17.11 | +17.11 | | 10 s | +| mt5, held | canceller alone | +11.05 | +11.05 | +11.05 | +11.05 | | 10 s | +| mt5, held | as shipped, decay 0.5 | +11.56 | +11.64 | +11.64 | +11.64 | +11.64 | 10 s | +| mt5, held | as shipped, decay 0.7 | +12.64 | +12.71 | +12.71 | +12.71 | +12.71 | 10 s | +| mt5, held | as shipped, decay 0.85 | +12.44 | +12.52 | +12.52 | +12.52 | +12.52 | 10 s | +| mt5, held | paper, decay 0.5 | +11.35 | +11.35 | +11.35 | +11.35 | +11.35 | 10 s | +| mt5, held | decay 0.7 + 5 Hz, whole bus | +18.50 | +18.57 | +18.57 | +18.57 | | 10 s | +| mt5, held | decay 0.7 + 5 Hz, dry only | +17.79 | +17.79 | +17.79 | +17.79 | | 10 s | +| cabin, speech | canceller alone | +21.02 | +20.82 | +20.23 | +19.75 | | **not converged** | +| cabin, speech | as shipped, decay 0.5 | +21.70 | +21.80 | +21.31 | +21.31 | +20.82 | **not converged** | +| cabin, speech | as shipped, decay 0.7 | +21.41 | +21.60 | +21.21 | +20.72 | +20.82 | 80 s | +| cabin, speech | as shipped, decay 0.85 | +20.53 | +20.53 | +20.33 | +20.04 | +19.94 | 80 s | +| mt5, speech | canceller alone | +23.26 | +23.16 | +22.87 | +22.68 | | 40 s | +| mt5, speech | as shipped, decay 0.5 | +24.53 | +24.73 | +24.14 | +24.24 | +24.14 | 40 s | +| mt5, speech | as shipped, decay 0.7 | +24.73 | +24.73 | +23.95 | +24.14 | +23.85 | **not converged** | +| mt5, speech | as shipped, decay 0.85 | +24.04 | +23.65 | +23.65 | +23.55 | +23.36 | 80 s | + +The canceller-alone rows and the shift rows ran to 80 s, the rest to 160 s. + +- **The held note converges at 10 s, decay 0.85 included.** Every held-note + row reads the same at 20, 40, 80 and 160 s, seed by seed. With the + canceller in the loop, the held note either runs away within the first + 20 s of a probe or not at all. Phase 0's 384 s was the dry loop with the + plate, no canceller; it does not carry over. +- **The speech envelope drifts down with probe length** and three of its + eight rows had not converged by their longest probe: the canceller alone + in the cabin moves 1.27 dB from 10 to 80 s. The karaoke suites saw the + same drift on this material. No speech-envelope row is gated; the grid + runs it on 80 s probes and reports it. +- The plate's cost converges faster than either limit: on speech it reads + between −0.10 and −1.86 dB at every probe, the plate raising the limit. + +The grid and the shift rows therefore run the held note on 20 s probes, +twice the converged length, and the gated rows on 10 s. + +## What the plate costs behind the canceller (runaway) + +`ReverbStageSweep.Grid`, six rooms, five seed sets each. Each cell is the +median over the six rooms of each room's five-seed median, with [the +range of the room medians] and {the per-seed range over all 30 runs}. +Cost is positive when the plate lowers the stable gain; matched is cost +minus the mix's level change (from [The stage](#the-stage)). Damping +0.0005 and the L return unless the row says otherwise. + +**Held note, S1 (10 ms), 20 s probes.** The canceller alone: ASG +11.74 +[+11.05, +12.42] {+7.73, +13.98}. + +| Plate, decay, wet | ASG | Cost [rooms] {seeds} | Mix level | Matched | +|---|---|---|---|---| +| as shipped 0.5, w 0.15 | +10.98 | +0.49 [−0.39, +2.05] {−4.00, +5.76} | −0.96 | +1.44 | +| as shipped 0.5, w 0.30 | +11.46 | +0.88 [−0.78, +2.44] {−4.30, +4.20} | −1.69 | +2.57 | +| as shipped 0.5, w 0.50 | +12.62 | −1.07 [−1.66, −0.39] {−4.30, +2.15} | −2.03 | +0.96 | +| as shipped 0.7, w 0.15 | +10.68 | +1.46 [−0.20, +1.76] {−3.91, +3.32} | −0.93 | +2.39 | +| as shipped 0.7, w 0.30 | +11.45 | +0.39 [−0.29, +0.68] {−4.00, +3.81} | −1.55 | +1.94 | +| as shipped 0.7, w 0.50 | +14.49 | −2.15 [−3.91, −1.86] {−4.98, +0.00} | −1.62 | −0.53 | +| as shipped 0.85, w 0.15 | +12.03 | −0.20 [−1.17, +2.34] {−4.30, +3.42} | −0.83 | +0.64 | +| as shipped 0.85, w 0.30 | +13.32 | −1.46 [−2.05, −0.49] {−6.35, +3.03} | −1.13 | −0.33 | +| as shipped 0.85, w 0.50 | +16.33 | −4.39 [−6.15, −3.52] {−8.50, −1.27} | −0.55 | −3.85 | +| paper 0.5, w 0.15 | +11.35 | +0.59 [+0.10, +1.46] {−4.30, +3.61} | −0.96 | +1.54 | +| paper 0.5, w 0.30 | +11.35 | +0.29 [−0.29, +1.46] {−4.79, +3.71} | −1.69 | +1.98 | +| paper 0.5, w 0.50 | +13.79 | −1.56 [−2.73, −1.07] {−6.35, +0.29} | −2.03 | +0.47 | +| as shipped 0.5, w 0.30, damping 0.5 | +11.84 | +0.39 [−0.78, +0.59] {−3.52, +3.91} | −1.71 | +2.10 | +| as shipped 0.5, w 0.30, (L + R)/2 | +11.54 | +0.29 [−0.78, +1.07] {−4.30, +5.08} | −1.70 | +2.00 | +| as shipped 0.7, w 0.30, damping 0.5 | +11.54 | +0.29 [−1.07, +0.49] {−4.20, +3.42} | −1.67 | +1.96 | +| as shipped 0.7, w 0.30, (L + R)/2 | +11.46 | +1.07 [−0.29, +1.56] {−5.86, +3.61} | −1.62 | +2.69 | +| as shipped 0.85, w 0.30, damping 0.5 | +13.20 | −1.37 [−2.15, +0.49] {−6.05, +2.73} | −1.58 | +0.21 | +| as shipped 0.85, w 0.30, (L + R)/2 | +12.81 | −1.17 [−3.12, −0.49] {−5.66, +1.27} | −1.32 | +0.15 | + +**Held note, S3 (20 ms), 20 s probes.** The canceller alone: ASG +18.40 +[+17.25, +18.59] {+16.46, +19.28}. + +| Plate, decay, wet | ASG | Cost [rooms] {seeds} | Mix level | Matched | +|---|---|---|---|---| +| as shipped 0.5, w 0.15 | +18.79 | −0.59 [−1.76, +1.17] {−3.22, +3.71} | −0.96 | +0.37 | +| as shipped 0.5, w 0.30 | +18.87 | −1.17 [−2.15, +0.49] {−4.20, +1.56} | −1.69 | +0.51 | +| as shipped 0.5, w 0.50 | +19.84 | −1.76 [−3.22, −1.17] {−4.39, +0.29} | −2.03 | +0.28 | +| as shipped 0.7, w 0.15 | +18.79 | −0.59 [−0.98, +0.20] {−3.71, +2.73} | −0.93 | +0.34 | +| as shipped 0.7, w 0.30 | +19.36 | −0.39 [−2.15, +0.10] {−4.10, +1.37} | −1.55 | +1.16 | +| as shipped 0.7, w 0.50 | +20.92 | −1.66 [−4.10, −0.68] {−4.49, +0.29} | −1.62 | −0.04 | +| as shipped 0.85, w 0.15 | +18.98 | −0.59 [−1.76, +0.68] {−2.64, +3.03} | −0.83 | +0.25 | +| as shipped 0.85, w 0.30 | +19.36 | −0.98 [−2.05, −0.39] {−3.61, +1.66} | −1.13 | +0.16 | +| as shipped 0.85, w 0.50 | +19.38 | −0.39 [−3.03, +0.98] {−5.27, +1.66} | −0.55 | +0.16 | +| paper 0.5, w 0.15 | +18.79 | −0.68 [−1.07, −0.10] {−2.83, +1.66} | −0.96 | +0.27 | +| paper 0.5, w 0.30 | +18.87 | −0.39 [−1.46, +0.68] {−3.32, +3.12} | −1.69 | +1.30 | +| paper 0.5, w 0.50 | +19.59 | −1.56 [−2.44, −0.29] {−4.98, +1.46} | −2.03 | +0.47 | +| as shipped 0.5, w 0.30, damping 0.5 | +18.30 | −0.10 [−1.37, +1.17] {−3.03, +2.25} | −1.71 | +1.62 | +| as shipped 0.5, w 0.30, (L + R)/2 | +19.06 | −0.49 [−1.17, +0.00] {−3.61, +2.54} | −1.70 | +1.22 | +| as shipped 0.7, w 0.30, damping 0.5 | +19.55 | −1.27 [−2.83, +0.59] {−3.81, +2.73} | −1.67 | +0.40 | +| as shipped 0.7, w 0.30, (L + R)/2 | +18.87 | −0.20 [−1.95, +0.29] {−4.49, +3.22} | −1.62 | +1.42 | +| as shipped 0.85, w 0.30, damping 0.5 | +19.65 | −1.46 [−1.95, −0.29] {−4.20, +1.27} | −1.58 | +0.11 | +| as shipped 0.85, w 0.30, (L + R)/2 | +19.28 | −0.98 [−3.22, −0.29] {−4.20, +2.44} | −1.32 | +0.34 | + +**Speech envelope, S1, 80 s probes** (not converged in every room; read as +directions). The canceller alone: ASG +19.75 [+18.48, +22.68] {+17.70, ++22.87}. + +| Plate, decay, wet | ASG | Cost [rooms] {seeds} | Mix level | Matched | +|---|---|---|---|---| +| as shipped 0.5, w 0.30 | +21.31 | −1.17 [−1.86, −0.68] {−2.44, +0.20} | −1.69 | +0.51 | +| as shipped 0.7, w 0.30 | +21.21 | −0.88 [−1.66, −0.78] {−1.86, +0.20} | −1.55 | +0.67 | +| as shipped 0.85, w 0.30 | +20.53 | −0.29 [−0.78, +0.10] {−1.95, +0.78} | −1.13 | +0.84 | + +What the grid says: + +- **Behind the canceller the plate costs little runaway gain, at any decay + measured.** The largest raw cost, median over the rooms, is +1.46 dB + (S1, held, decay 0.7, w 0.15), and +2.69 dB at matched loudness (S1, + held, decay 0.7, w 0.30, (L + R)/2). At S3 no raw cost is positive. + Where the bare loop (no canceller, the four fixtures at 4096 taps) + charges decay 0.85 at w 0.30 +3.71 dB by the bound, the canceller-held + loop reads −1.46 (S1). +- **Long decays and high wet raise the runaway limit on the held note.** + At S1, decay 0.85, w 0.50 the chain holds 4.39 dB more than the + canceller alone (−3.85 at matched loudness), in every room ({−8.50, + −1.27} per seed). Runaway is not the whole story: the audible rows below + fall as wet and decay rise. +- **Damping, the return and the plate do not separate.** At w 0.30, + damping 0.5 against 0.0005 moves the cost by −0.49 / −0.10 / +0.09 dB at + S1 (decay 0.5 / 0.7 / 0.85) and +1.07 / −0.88 / −0.48 at S3; the + (L + R)/2 return against L by −0.59 / +0.68 / +0.29 at S1. The paper + plate at decay 0.5 sits within 0.78 dB of the as-shipped plate at every + wet and delay. Every one of these differences is smaller than the + per-room ranges of the two rows it compares. +- **Per-seed spread dominates:** single seeds run from −8.50 to +5.76 dB + of cost on the held note. The gated row + (`ReverbStage.CancellerHoldsBehindThePlates`, cabin, S1) bounds the + median cost at +3.5 dB. + +## Shift and reverb + +With a frequency shift in the decorrelator slot and the plate behind it, +the plate's tail goes round the loop through the shifter and climbs by the +shift on every pass. `shifted_dry_mix` is the alternative: the shift on the +dry path only, the plate fed the unshifted bus. Both are measured here +against the shift alone, all on the library's shifter. + +### Runaway (bisected), six rooms + +`ReverbStageSweep.Shift`: held note, 20 s probes, the as-shipped plate at +wet 0.30. ASG is the median of the six rooms' five-seed medians, [the range +of the room medians]; whole-bus − dry-only is the per-seed difference of +the two chain limits, its median per room, then the median of those, with +the count of rooms where it is negative: + +| S1 | 2 Hz | 5 Hz | +|---|---|---| +| canceller alone | +11.74 [+11.05, +12.42] | (same) | +| shift alone | +17.50 [+15.18, +19.06] | +17.23 [+15.37, +20.53] | +| plate alone, decay 0.5 | +11.46 [+9.98, +11.86] | (same) | +| whole bus + plate 0.5 | +17.52 [+16.82, +19.75] | +19.36 [+16.54, +19.65] | +| dry only + plate 0.5 | +14.10 [+13.01, +15.25] | +16.91 [+15.84, +18.28] | +| whole bus − dry only, decay 0.5 | +4.49 (0 of 6 rooms < 0) | +1.46 (0 of 6) | +| plate alone, decay 0.7 | +11.45 [+10.66, +12.71] | (same) | +| whole bus + plate 0.7 | +17.99 [+16.74, +19.45] | +18.57 [+14.88, +18.77] | +| dry only + plate 0.7 | +15.45 [+13.91, +16.23] | +17.11 [+16.15, +18.77] | +| whole bus − dry only, decay 0.7 | +3.42 (0 of 6) | +0.68 (1 of 6; −0.98) | + +| S3 | 2 Hz | 5 Hz | +|---|---|---| +| canceller alone | +18.40 [+17.25, +18.59] | (same) | +| shift alone | +19.38 [+16.74, +21.60] | +18.89 [+17.52, +21.31] | +| plate alone, decay 0.5 | +18.87 [+18.30, +20.25] | (same) | +| whole bus + plate 0.5 | +20.45 [+18.89, +22.58] | +20.33 [+19.00, +21.80] | +| dry only + plate 0.5 | +18.18 [+17.23, +19.08] | +20.16 [+18.98, +21.21] | +| whole bus − dry only, decay 0.5 | +1.95 (0 of 6) | +0.20 (3 of 6) | +| plate alone, decay 0.7 | +19.36 [+17.73, +19.75] | (same) | +| whole bus + plate 0.7 | +20.35 [+19.20, +22.38] | +20.25 [+19.00, +21.70] | +| dry only + plate 0.7 | +19.67 [+18.40, +20.04] | +20.45 [+18.32, +22.29] | +| whole bus − dry only, decay 0.7 | +1.17 (1 of 6) | +0.68 (3 of 6) | + +- **By runaway the whole-bus shift holds more than the dry-only shift at + S1**: its room medians are positive in 22 of 24, zero in one and + negative in one (−0.98, 5 Hz, decay 0.7). At S3 the same holds at 2 Hz + (11 of 12 positive); at 5 Hz the two do not separate (3 of 6 rooms each + way at both decays). The dry-only shift leaves the wet path + (gain 0.30) unshifted and gives back part of the shift's runaway gain: + at 2 Hz it reads +14.10 / +15.45 against the shift alone's +17.50. +- **By the audible limit it is the other way round** (next section): on + the cabin ramp the dry-only shift is audible later in 13 of 16 pairs. + The two measurements differ in protocol as well as criterion (a converged + canceller bisected against a canceller adapting from a cold start under + a ramp; karaoke-afc.md, "Ramp vs bisected runaway"), and the ramp's own + runaway gain also favours the dry-only shift (16 of 16). Which one a room + shows is phase 1's to measure. +- `ReverbStage.WholeBusShiftHoldsMoreRunawayThanDryOnly` gates the + direction in the cabin at 2 Hz, decay 0.5 (per-seed median +3.09 dB). + +### Audible limits: the product-chain candidates + +The audible limit on the canceller output c, by the offline howl criterion +under the protocol's ramp (30 s warm-up 20 dB under the open-loop limit, +then +1 dB per 2 s; the live canceller from a cold start), with +PROTOCOL.md §7.3's flags: `--rt60 0.0328` (the band-limited cabin) on every +row, `--chain-rt60` the plate's T30 (1.1846 s at decay 0.5, 2.1167 s at 0.7, +from `measure_rir.py`'s `schroeder()`) on every row with the plate, and +`--dechirp` on every row with a shift. The reverb-only rows are analysed on +the chain output (after the plate, before the gain), with c as the +cross-check. The cabin, held note, the as-shipped plate, the library's +shifter, five seeds; dB over the dry loop's `exact_msg_db` (S1 −6.8967, +S3 −6.8603); medians (`tools/notebook/reverb_audible.py`): + +| Chain | S1 audible | S1 runaway (ramp) | S3 audible | S3 runaway (ramp) | +|---|---|---|---|---| +| dry (no canceller) | +0.18 | +1.08 | +0.21 | +1.62 | +| canceller | +16.57 | +20.70 | +12.37 | +21.37 | +| + 2 Hz † | +4.03 | +20.13 | +6.50 | +18.73 | +| + 5 Hz † | +4.43 | +19.39 | +2.18 | +19.49 | +| + plate 0.5, w 0.15 | +16.53 | +21.46 | +10.38 | +23.01 | +| + plate 0.5, w 0.30 | +12.89 | +22.07 | −4.77 ‡ | +23.31 | +| + plate 0.7, w 0.15 | +13.53 | +20.89 | +11.44 | +23.18 | +| + plate 0.7, w 0.30 | +10.46 | +21.76 | +3.62 ‡ | +24.19 | +| + 2 Hz whole bus + plate 0.5, w 0.15 † | +5.92 | +20.71 | +6.77 | +20.68 | +| + 2 Hz dry only + plate 0.5, w 0.15 † | +8.55 | +21.31 | +7.90 | +21.63 | +| + 5 Hz whole bus + plate 0.5, w 0.15 † | +6.66 | +20.93 | +2.95 | +21.01 | +| + 5 Hz dry only + plate 0.5, w 0.15 † | +6.48 | +21.28 | +3.53 | +22.24 | +| + 2 Hz whole bus + plate 0.5, w 0.30 † | +7.53 | +21.70 | +7.74 | +22.62 | +| + 2 Hz dry only + plate 0.5, w 0.30 † | +8.51 | +22.24 | +8.98 | +22.62 | +| + 5 Hz whole bus + plate 0.5, w 0.30 † | +6.43 | +22.12 | +5.71 | +22.91 | +| + 5 Hz dry only + plate 0.5, w 0.30 † | +8.06 | +22.66 | +7.10 | +23.55 | +| + 2 Hz whole bus + plate 0.7, w 0.15 † | +6.62 | +20.85 | +7.50 | +18.99 | +| + 2 Hz dry only + plate 0.7, w 0.15 † | +8.71 | +20.85 | +8.68 | +22.12 | +| + 5 Hz whole bus + plate 0.7, w 0.15 † | +6.93 | +20.45 | +3.12 | +20.26 | +| + 5 Hz dry only + plate 0.7, w 0.15 † | +7.13 | +21.75 | +4.45 | +21.12 | +| + 2 Hz whole bus + plate 0.7, w 0.30 † | +8.87 | +20.47 | +6.74 | +21.31 | +| + 2 Hz dry only + plate 0.7, w 0.30 † | +7.69 | +21.47 | +6.19 | +22.67 | +| + 5 Hz whole bus + plate 0.7, w 0.30 † | +5.73 | +21.67 | +4.63 | +22.15 | +| + 5 Hz dry only + plate 0.7, w 0.30 † | +10.30 | +22.62 | +9.46 | +24.42 | + +† **Worst case for the shift**, as in the karaoke document: the held note is +a perfectly periodic synthetic note with a 40 dB floor between its +harmonics, so partials recirculating through the shifter stand out. + +‡ The reverb-only rows at S3, wet 0.30, are audible below the dry loop's +limit in some seeds: per seed −7.89 −0.17 +1.89 −4.77 −13.55 at decay 0.5 +(c cross-check −2.48 −0.18 +1.90 +6.39 −13.55) and +4.56 +6.94 −7.15 +3.62 +−2.56 at 0.7. Their runaway limits are the highest in the table. What the +criterion flags there was not examined. + +What the audible rows say: + +- **The dry-only shift does better than the whole-bus shift.** Per seed, + dry-only minus whole-bus, the median is positive in 13 of the 16 + shift + plate pairs for the audible limit and in 16 of 16 for the ramp's + runaway; the three exceptions are 2 Hz at decay 0.7, wet 0.30 (S1 −1.97, + S3 −0.41) and 5 Hz at decay 0.5, wet 0.15, S1 (−0.18). The largest gains + are at 5 Hz, decay 0.7, wet 0.30: +4.34 (S1) and +5.11 (S3). +- **The plate alone costs the canceller audible gain as wet and decay + rise:** at S1 +16.53 / +12.89 at decay 0.5 (w 0.15 / 0.30) and +13.53 / + +10.46 at 0.7, against the canceller's +16.57. Its runaway limit does not + fall (+20.89 to +22.07 against +20.70). +- **Adding the plate to the shift does not lower the audible limit.** The + shift + plate rows read +5.73 to +10.30 at S1 against +4.03 / +4.43 for + the shift alone, but they are judged with the plate's T30 in the + criterion's hold and the shift-alone rows are not. Re-analysed with the + same hold (`--control`), the shift alone reads +5.84 / +6.10 at 2 Hz and + +5.28 / +5.37 at 5 Hz (S1; 1.185 / 2.117 s in the hold), and at S3 +6.50 + and +2.18 (unchanged). +- Seed spread is large: with the canceller in, per-seed audible − runaway + ranges from −2.53 to −36.79 dB, and a single seed reads up to 21.60 dB + from its row's median (2 Hz alone, S1, seed 82: −17.57 against +4.03). + +## The bare loop: what the plate costs without a canceller + +Phase 0's method, reproduced for both plates (`ReverbStageSweep.BareLoopCost`): +the loop K·e^{−jωd}·H(ω)·F(ω) with H = (1 − w) + w·R the mix and F one of +the four room fixtures at its full 4096 taps (studio, rehearsal, hall, +cabin; not band-limited), d = 480 samples (S1). Each cell is the cost in +open-loop stable gain, dry limit minus mix limit, positive = less stable +gain, as **magnitude bound / exact Nyquist crossing**: the bound +(−20·log10 max|HF|) holds under any delay; the exact limit is the first K at +which K·L(ω) = 1 at d = 480. Median over the four rooms (the mean of the +middle two, as phase 0); L return. "Matched" is the w 0.30 bound minus the +mix's level change at w 0.30. Spectra on an rfft grid of 2N, N = +nextpow2(4·(IR length + 4095 + 480)); over all 144 cells (both plates, both +returns) the cost moves at most 0.0309 dB (bound) and 0.0693 dB (exact) +from N to 2N. + +As shipped (reproduces phase 0's table to the printed digit): + +| decay (T30 L) | damping | w 0.15 | w 0.30 | w 0.50 | matched (w 0.30) | +|---|---|---|---|---|---| +| 0.3 (1.01 s) | 0.0005 | −0.16 / −0.56 | −0.21 / −0.83 | −0.16 / −0.66 | +1.51 | +| 0.5 (1.18 s) | 0.0005 | +0.01 / −0.62 | +0.12 / −0.60 | +0.42 / −0.09 | +1.81 | +| 0.5 | 0.5 | −0.09 / −0.67 | −0.09 / −0.83 | +0.03 / −0.54 | +1.63 | +| 0.7 (2.12 s) | 0.0005 | +0.46 / −0.28 | +1.11 / +0.45 | +2.17 / +1.85 | +2.65 | +| 0.7 | 0.5 | +0.03 / −0.74 | +0.18 / −0.67 | +0.96 / −0.07 | +1.85 | +| 0.85 (4.62 s) | 0.0005 | +1.78 / +0.65 | +3.71 / +2.76 | +5.79 / +6.25 | +4.84 | +| 0.85 | 0.5 | +0.50 / −0.55 | +2.22 / +1.18 | +4.26 / +3.12 | +3.79 | + +**Paper lengths (`dattorro_paper`), new:** + +| decay (T30 L) | damping | w 0.15 | w 0.30 | w 0.50 | matched (w 0.30) | +|---|---|---|---|---|---| +| 0.3 (1.73 s) | 0.0005 | −0.21 / −0.50 | −0.27 / −0.51 | −0.20 / −0.41 | +1.45 | +| 0.5 (2.12 s) | 0.0005 | +0.08 / −0.31 | +0.30 / +0.06 | +0.74 / +0.43 | +1.99 | +| 0.5 | 0.2 | −0.03 / −0.34 | +0.09 / −0.01 | +0.37 / +0.38 | +1.79 | +| 0.5 | 0.5 | −0.16 / −0.37 | −0.18 / −0.16 | −0.00 / +0.15 | +1.53 | +| 0.7 (3.79 s) | 0.0005 | +0.76 / +0.38 | +1.61 / +1.14 | +2.77 / +2.42 | +3.16 | +| 0.7 | 0.2 | +0.34 / +0.22 | +0.80 / +0.87 | +1.53 / +1.99 | +2.42 | +| 0.7 | 0.5 | −0.02 / −0.04 | +0.14 / +0.37 | +0.98 / +0.97 | +1.81 | +| 0.85 (8.25 s) | 0.0005 | +2.30 / +1.43 | +4.32 / +4.00 | +6.59 / +6.47 | +5.45 | +| 0.85 | 0.2 | +0.91 / +0.88 | +2.78 / +3.11 | +4.92 / +5.35 | +4.19 | +| 0.85 | 0.5 | +0.33 / +0.13 | +2.04 / +1.42 | +3.92 / +4.12 | +3.62 | + +What the bare-loop table says, for a listening pass that has to pick a +decay and a wet: + +- **At equal decay the paper plate costs more**, because it rings longer: + at w 0.30, damping 0.0005, the bound reads +0.30 against +0.12 at decay + 0.5, +1.61 against +1.11 at 0.7, +4.32 against +3.71 at 0.85. +- **At equal T30 the paper plate costs less.** At decay 0.5 (T30 2.12 s) + it reads +0.08 / +0.30 / +0.74 by the bound at w 0.15 / 0.30 / 0.50; the + as-shipped plate at decay 0.7, the same T30, reads +0.46 / +1.11 / +2.17. + Why was not measured (phase 0 measured the as-shipped plate's |R| peak + only); nor does the bare loop say how the two sound. +- **Damping trims the long decays, more on the paper plate.** At w 0.30 + the bound falls from +3.71 to +2.22 dB (as shipped, decay 0.85, damping + 0.0005 → 0.5) and from +4.32 to +2.04 (paper, decay 0.85); at decay 0.5 + it moves the as-shipped plate by 0.21 dB and the paper plate by 0.48. +- **The (L + R)/2 return never costs more by the bound:** at w 0.30 it + reads 0.00 to 0.39 dB below L in all 24 rows (plate × decay × damping); + by the exact crossing from 0.73 dB below to 0.14 above. The sweep prints + its full table. + +## What did not converge, separate or get measured + +- **The speech envelope did not converge** in three of eight convergence + rows by the longest probe (cabin canceller to 80 s; cabin decay 0.5 and + mt5 decay 0.7 to 160 s). Its grid rows (80 s) are directions, not gated. +- **Damping, the return and the plate do not separate** behind the + canceller (the grid); only the bare loop separates them. +- **The two shift topologies disagree by measurement:** the whole-bus + shift holds more runaway gain by bisection, the dry-only shift is audible + later on the ramp. Neither protocol is wrong; which a room shows is a + phase 1 measurement. +- **The reverb-only rows at S3, wet 0.30, are audible below the dry loop's + limit in some seeds** (the audible table, ‡). Not examined. +- **Not run:** the speech envelope at S3; the paper plate in the loop at + decay 0.7 and 0.85 (the bare loop covers them); the shift rows at wet + 0.15 and 0.50 by bisection (the audible driver covers 0.15); the + audible rows for the paper plate, for other rooms than the cabin and for + the speech envelope; float32 in the loop (every in-loop number is + double; `test_reverb_mix.cpp` runs the stage's float plumbing on every + target); the backing track; two microphones. +- **The audible table's canceller row is not the karaoke document's.** + karaoke-afc.md's headline analysed every row with `--dechirp`; here + PROTOCOL.md §7.3's flags apply (no `--dechirp` without a shift), and the + canceller alone reads +16.57 at S1 (karaoke-afc.md: +14.96 on this host). + The loop itself is unchanged: `decorrelated_loop` with no stage gives the + karaoke measurement's limits bit for bit against origin/main (cabin and + mt5, held and speech, with and without the shift, seed 2, 5 s probes). + +## Hosts + +Every number here was measured on **macOS 15.7.9 x86_64 (i9-8950HK), +AppleClang 17.0.0, Release**, double precision in the loop, on a machine +shared with another job (at most four of its twelve threads were this +work's). The analytic tables (T30, mix level, bare loop) are deterministic; +the loop is chaotic, which is why the in-loop claims are medians over five +seed sets and six rooms, and the gated rows bound medians with several dB +of margin. No second host has run the sweep or the audible driver. CI runs +the gated rows on its macOS arm64 and Linux legs (pass/fail only). The +audible driver ran on Python 3.10.7, numpy 2.2.6, scipy 1.15.3. + +## Provenance + +```sh +cmake -S . -B build -G Ninja -DCMAKE_BUILD_TYPE=Release -DMUTAP_WERROR=ON -DMUTAP_BUILD_KARAOKE_DUMP=ON +cmake --build build +``` + +The stage's plumbing (every target; both emulated selections run the float +half) and the gated host rows: + +```sh +build/tests/mutap_tests --gtest_filter='reverb_mix_test*:ReverbMix*' +build/tests/mutap_tests --gtest_filter='ReverbStage.*' +``` + +Wall times in the full `ctest` run of this PR (396 of 396 passed, 3233.67 s, +serial; each in-loop row on 4 threads): `CancellerHoldsBehindThePlates` +133.29 s, `WholeBusShiftHoldsMoreRunawayThanDryOnly` 63.73 s, +`PaperPlateBareLoopCost` 1.35 s, `PlateDecayAndMixLevel` 0.34 s, the rest +0.01 s each. + +The sweep, behind `MUTAP_SLOW=1` (at most 4 threads by default; +`MUTAP_SLOW_THREADS` overrides; every finished bisection prints a `job` +line): + +```sh +MUTAP_SLOW=1 MUTAP_SLOW_THREADS=3 build/tests/mutap_tests --gtest_filter='ReverbStageSweep.ProbeConvergence' +MUTAP_SLOW=1 MUTAP_SLOW_THREADS=3 build/tests/mutap_tests --gtest_filter='ReverbStageSweep.Grid' +MUTAP_SLOW=1 MUTAP_SLOW_THREADS=1 REVERB_SWEEP_DELAYS=480 build/tests/mutap_tests --gtest_filter='ReverbStageSweep.Shift' +MUTAP_SLOW=1 MUTAP_SLOW_THREADS=4 REVERB_SWEEP_DELAYS=960 build/tests/mutap_tests --gtest_filter='ReverbStageSweep.Shift' +MUTAP_SLOW=1 build/tests/mutap_tests --gtest_filter='ReverbStageSweep.BareLoopCost' +``` + +| Test | Wall time | Threads | Bisections | +|---|---|---|---| +| `ProbeConvergence` | 22607.7 s (cpu-sum 67810.6 s) | 3 | 610 | +| `Grid` | 24684.3 s (cpu-sum 73998.3 s) | 3 | 1350 | +| `Shift`, S1 | 17611.7 s (cpu-sum 17611.6 s) | 1 | 420 | +| `Shift`, S3 | 4813.9 s (cpu-sum 19196.2 s) | 4 | 420 | +| `BareLoopCost` | 53.5 s | 1 | (analytic) | + +The audible table needs numpy and scipy; each run is about 2 minutes of +dump and one to three analyses. One invocation with `--control` reproduces +all of it (the runs are deterministic): + +```sh +python3 tools/notebook/reverb_audible.py --dump build/tools/notebook/karaoke_ramp_dump --work runs --jobs 1 \ + --control --json audible.json +python3 tools/notebook/reverb_audible.py --merge audible.json --work runs # reprint +``` + +How this table was actually produced: a first invocation without +`--control` ran 230 of its 240 runs on one job in about 5 h 10 min and +then died on a write: it kept every run's WAVs (about 180 MB a run) and +the scratch disk filled. The driver now deletes each run's WAVs after the +analysis (`--keep-wav` keeps them) and writes each run's record as it +finishes. `--recover` rebuilt the 230 records from the criterion's JSON +and the gain logs (the runaway gain by the dump's own rule; before the WAVs +were deleted, a scratch check found it equal to the WAVs' first 40 dB block +in 230 of 230 runs); a resumed invocation ran the last 10 runs (786.2 s); +the 20 shift-only runs were re-run with `--control` (3104.9 s), and their +limits repeat the first run's to the printed digit. diff --git a/include/mutap/mutap.h b/include/mutap/mutap.h index e78b945..b3b3977 100644 --- a/include/mutap/mutap.h +++ b/include/mutap/mutap.h @@ -24,3 +24,4 @@ #include "mutap/lpc.h" #include "mutap/pem_afc.h" #include "mutap/postfilter.h" +#include "mutap/reverb_stage.h" diff --git a/include/mutap/reverb_stage.h b/include/mutap/reverb_stage.h new file mode 100644 index 0000000..80ea8ef --- /dev/null +++ b/include/mutap/reverb_stage.h @@ -0,0 +1,296 @@ +/// @file reverb_stage.h +/// @brief Reverb-slot stages for afc_chain: a dry/wet mix around any mono-in +/// reverb, and the variant that frequency-shifts the dry path only. +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors +// +// WHAT. afc_chain's reverb slot (decorrelator -> reverb -> safety, see +// afc_chain.h) takes any noexcept process_block(in, out, n). A reverb +// engine does not have that shape: it is mono in, N >= 1 channels out, and +// says nothing about how much of it reaches the loudspeaker. The two +// classes here are that glue, with the mix written down: +// +// reverb_mix: y = (1 - w) x + w r(x) +// shifted_dry_mix: y = (1 - w) shift(x) + w r(x) +// +// where r is the reverb's mono return (output channel 0, "L", or the +// average of channels 0 and 1, "(L + R) / 2"; reverb_return) and w the wet +// level in [0, 1]. The mix is a crossfade, so it changes the bus level: for +// white input its power gain is 10 log10 of the sum of the squared impulse +// response of (1 - w) delta + w r: for the vendored Dattorro plates, +// -0.547 to -2.135 dB over wet 0.15 to 0.50 and decay 0.3 to 0.85 (L +// return, damping 0.0005; docs/reverb-afc.md). A loop measurement that +// compares a mix against the dry chain at equal forward gain therefore +// includes that level change; docs/reverb-afc.md reports both. +// +// THE DRY-ONLY SHIFT. With a frequency shifter in the decorrelator slot +// AND a reverb behind it, the reverb's tail goes round the loop through the +// shifter and climbs by the shift on every pass. shifted_dry_mix is the +// alternative topology for that case: the shifter sits on the dry +// path inside this stage (leave the chain's decorrelator slot empty), and +// the reverb is fed the UNshifted bus, so the wet path never passes through +// the shifter. What it costs and buys in the loop is measured, not assumed +// (docs/reverb-afc.md, "Shift and reverb"). +// +// THE REVERB TYPE. Anything with +// +// void process(const Sample* const* in, Sample* const* out, int n) noexcept +// +// taking ONE input channel and writing `channels` output channels: the shape +// of the FAUST shim's mutap_faust::faust_block (third_party/faust/), which +// the tests instantiate with the vendored Dattorro plates. FAUST-generated +// code never enters include/mutap: these classes know only that signature. +// If the type also has num_outputs() (faust_block does), the constructor +// checks it against config::channels. The stages do not own the reverb (as +// afc_chain does not own its stages): the caller keeps it alive, sets its +// parameters (decay, damping, ...) through its own interface and clears it +// with its own reset when the chain is reset. +// +// LATENCY. latency() is 0 for both: the mix adds no delay of its own. A +// reverb with a pre-delay delays only the wet path, and a shifter's +// frequency-dependent group delay is not a bulk latency +// (frequency_shifter::latency()). +// +// Real-time contract (as fd_kalman.h): the constructor validates its config, +// allocates the per-channel return buffers (channels x block_size samples, +// plus block_size for the shifted dry path) and may throw +// std::invalid_argument. process_block() and every other post-construction +// entry point are noexcept and allocation-free. process_block() accepts any +// n (it works in pieces of at most block_size) and may run in place +// (in == out). Arithmetic is in Sample throughout: float32 is first-class, +// and no double runs in the float hot path. + +#pragma once + +#include +#include +#include +#include + +// No ABI tag: nothing here holds an FFT, so the layout does not depend on the +// build's FFT engine and the classes sit in plain tap::mu (as +// frequency_shifter.h's do). Header-only; the Reverb and Shifter are +// template arguments. +namespace tap::mu { + + /// Which of the reverb's outputs the mix returns to the mono bus. + enum class reverb_return { + left, ///< output channel 0 + mid ///< (channel 0 + channel 1) / 2; needs >= 2 channels + }; + + /// A mono-in reverb engine the stages can drive (see the file comment). + template + concept mono_reverb = requires(Reverb& reverb, const Sample* const* in, Sample* const* out, int n) { + { reverb.process(in, out, n) } noexcept; + }; + + /// A block process the dry path of shifted_dry_mix can run (the shape of + /// afc_chain's afc_stage: tap::mu::frequency_shifter, for one). + template + concept dry_path_stage = requires(Shifter& stage, const Sample* in, Sample* out, size_t n) { + { stage.process_block(in, out, n) } noexcept; + }; + + namespace reverb_detail { + + /// The return half both stages share: the reverb's output buffers, + /// the wet level and the return mode. + template + class reverb_return_core { + public: + reverb_return_core(Reverb& reverb, size_t block_size, size_t channels, Sample wet, reverb_return mode) + : m_reverb(&reverb) + , m_block(validated_block(block_size)) + , m_channels(validated_channels(reverb, channels)) + , m_buffer(channels * block_size, Sample(0)) + , m_outs(channels, nullptr) { + for (size_t c = 0; c < channels; ++c) { + m_outs[c] = m_buffer.data() + c * block_size; + } + if (!set_wet(wet)) { + throw std::invalid_argument("reverb_mix: wet must be finite and in [0, 1]"); + } + if (!set_return_mode(mode)) { + throw std::invalid_argument("reverb_mix: the (L + R) / 2 return needs >= 2 reverb channels"); + } + } + + bool set_wet(Sample wet) noexcept { + if (!std::isfinite(wet) || wet < Sample(0) || wet > Sample(1)) { + return false; + } + m_wet = wet; + m_dry = Sample(1) - wet; + return true; + } + + bool set_return_mode(reverb_return mode) noexcept { + if (mode == reverb_return::mid && m_channels < 2) { + return false; + } + m_mode = mode; + return true; + } + + Sample wet() const noexcept { return m_wet; } + Sample dry_gain() const noexcept { return m_dry; } + reverb_return return_mode() const noexcept { return m_mode; } + size_t block_size() const noexcept { return m_block; } + size_t channels() const noexcept { return m_channels; } + + /// Run the reverb on `n` <= block_size() samples of `in`; its + /// outputs land in the channel buffers. + void run_reverb(const Sample* in, size_t n) noexcept { + const Sample* const ins[1] = {in}; + m_reverb->process(ins, m_outs.data(), static_cast(n)); + } + + /// out[i] = dry_gain() * dry[i] + wet() * r[i] for the last + /// run_reverb() block; dry may equal out. + void mix(const Sample* dry, Sample* out, size_t n) const noexcept { + const Sample* l = m_outs[0]; + if (m_mode == reverb_return::mid) { + const Sample* r = m_outs[1]; + const Sample half = m_wet * Sample(0.5); + for (size_t i = 0; i < n; ++i) { + out[i] = m_dry * dry[i] + half * (l[i] + r[i]); + } + } + else { + for (size_t i = 0; i < n; ++i) { + out[i] = m_dry * dry[i] + m_wet * l[i]; + } + } + } + + private: + static size_t validated_block(size_t block_size) { + if (block_size == 0) { + throw std::invalid_argument("reverb_mix: block_size must be >= 1"); + } + return block_size; + } + + static size_t validated_channels(Reverb& reverb, size_t channels) { + if (channels == 0) { + throw std::invalid_argument("reverb_mix: channels must be >= 1"); + } + if constexpr (requires(const Reverb& r) { r.num_outputs(); }) { + if (static_cast(reverb.num_outputs()) != channels) { + throw std::invalid_argument("reverb_mix: channels differs from the reverb's num_outputs()"); + } + } + if constexpr (requires(const Reverb& r) { r.num_inputs(); }) { + if (reverb.num_inputs() != 1) { + throw std::invalid_argument("reverb_mix: the reverb must take exactly one input"); + } + } + return channels; + } + + Reverb* m_reverb; + size_t m_block; + size_t m_channels; + std::vector m_buffer; ///< channel-major, block_size each + std::vector m_outs; ///< one pointer per channel into m_buffer + Sample m_wet = Sample(0); + Sample m_dry = Sample(1); + reverb_return m_mode = reverb_return::left; + }; + + } // namespace reverb_detail + + /// The reverb slot's mix: y = (1 - w) x + w r(x), r the reverb's mono + /// return. An afc_stage (afc_chain::set_reverb(&mix)). Does not own the + /// reverb. See the file comment for the contract. + template + requires mono_reverb + class reverb_mix { + public: + struct config { + size_t block_size = 64; ///< the largest piece process_block() hands the reverb; >= 1 + size_t channels = 2; ///< the reverb's output channels; >= 1 + Sample wet = Sample(0.3); ///< w, finite, in [0, 1] + reverb_return mode = reverb_return::left; + }; + + /// @throws std::invalid_argument on a bad config, or when the + /// reverb reports a channel count other than cfg.channels + reverb_mix(Reverb& reverb, const config& cfg) + : m_core(reverb, cfg.block_size, cfg.channels, cfg.wet, cfg.mode) {} + + /// Set w; a non-finite w, or one outside [0, 1], is rejected (false, + /// nothing changes). Takes effect from the next block. + bool set_wet(Sample wet) noexcept { return m_core.set_wet(wet); } + /// Choose the return; mid with a one-channel reverb is rejected. + bool set_return_mode(reverb_return mode) noexcept { return m_core.set_return_mode(mode); } + + Sample wet() const noexcept { return m_core.wet(); } + reverb_return return_mode() const noexcept { return m_core.return_mode(); } + size_t block_size() const noexcept { return m_core.block_size(); } + size_t channels() const noexcept { return m_core.channels(); } + + /// 0 samples: the mix adds no delay of its own (file comment). + static constexpr size_t latency() noexcept { return 0; } + + /// Mix n samples. in == out is allowed. + void process_block(const Sample* in, Sample* out, size_t n) noexcept { + const size_t b = m_core.block_size(); + for (size_t off = 0; off < n; off += b) { + const size_t len = (n - off < b) ? n - off : b; + m_core.run_reverb(in + off, len); + m_core.mix(in + off, out + off, len); + } + } + + private: + reverb_detail::reverb_return_core m_core; + }; + + /// The dry-only-shift topology: y = (1 - w) shift(x) + w r(x), the + /// reverb fed the unshifted input so its tail never recirculates through + /// the shifter. Put it in the chain's REVERB slot and leave the + /// decorrelator slot empty. An afc_stage. Owns neither the shifter nor + /// the reverb. See the file comment for the contract. + template + requires dry_path_stage && mono_reverb + class shifted_dry_mix { + public: + using config = typename reverb_mix::config; + + /// @throws std::invalid_argument as reverb_mix's constructor + shifted_dry_mix(Shifter& shifter, Reverb& reverb, const config& cfg) + : m_shifter(&shifter) + , m_core(reverb, cfg.block_size, cfg.channels, cfg.wet, cfg.mode) + , m_dry(cfg.block_size, Sample(0)) {} + + bool set_wet(Sample wet) noexcept { return m_core.set_wet(wet); } + bool set_return_mode(reverb_return mode) noexcept { return m_core.set_return_mode(mode); } + + Sample wet() const noexcept { return m_core.wet(); } + reverb_return return_mode() const noexcept { return m_core.return_mode(); } + size_t block_size() const noexcept { return m_core.block_size(); } + size_t channels() const noexcept { return m_core.channels(); } + + /// 0 samples (file comment; the shifter's group delay is its own). + static constexpr size_t latency() noexcept { return 0; } + + /// Mix n samples. in == out is allowed. + void process_block(const Sample* in, Sample* out, size_t n) noexcept { + const size_t b = m_core.block_size(); + for (size_t off = 0; off < n; off += b) { + const size_t len = (n - off < b) ? n - off : b; + m_core.run_reverb(in + off, len); + m_shifter->process_block(in + off, m_dry.data(), len); + m_core.mix(m_dry.data(), out + off, len); + } + } + + private: + Shifter* m_shifter; + reverb_detail::reverb_return_core m_core; + std::vector m_dry; ///< the shifted dry path, one piece + }; + +} // namespace tap::mu diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 61e7451..0d76d43 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -59,7 +59,8 @@ add_executable(mutap_tests test_float32.cpp test_nn_suppressor.cpp test_fft_engine_contract.cpp - test_portable_random.cpp) + test_portable_random.cpp + test_reverb_mix.cpp) target_link_libraries(mutap_tests PRIVATE MuTap::MuTap mutap_warnings) @@ -88,6 +89,9 @@ if(NOT MUTAP_BARE_METAL AND NOT CMAKE_CROSSCOMPILING) # The guard's live claims (tests/support/guard_runs.h): host-only for the # same reasons, and one row runs the vendored Dattorro in the reverb slot. target_sources(mutap_tests PRIVATE test_howl_guard_host.cpp) + # The reverb behind the canceller (docs/reverb-afc.md): host-only for the + # same reasons, and it drives the vendored Dattorro plates. + target_sources(mutap_tests PRIVATE test_reverb_stage.cpp test_reverb_stage_sweep.cpp) target_include_directories(mutap_tests PRIVATE ${PROJECT_SOURCE_DIR}/third_party/faust) endif() @@ -117,13 +121,15 @@ else() # validation/contract tests, the frequency shifter's float suite (not # its 10-minute drift test) and the howl detector's float suite with its # canary (its acoustic claims are test_howl_detector_host.cpp, host-only), - # and the guard's float state-machine suite (its live claims are - # test_howl_guard_host.cpp, host-only). + # the guard's float state-machine suite (its live claims are + # test_howl_guard_host.cpp, host-only), and the reverb stage's float + # plumbing (its Dattorro and in-loop claims are test_reverb_stage.cpp, + # host-only). # (The FFT suites left with the FFT: they run in DspTap's own CI.) set(mutap_emulated_selection "") if(CMAKE_CROSSCOMPILING) set(mutap_emulated_selection TEST_FILTER - "fdaf_test/0.*:FdafCrossPrecision.*:FdafConfigValidation.*:FdafRtContract.*:fd_kalman_test/0.*:fd_kalman_test/1.*:kalman_loop_test/0.*:FdKalmanConfigValidation.*:FdKalmanRtContract.*:Levinson.*:LpcPredictor.*:SpeechPredictor.*:WarpedLpcPredictor.*:PredictorConfigValidation.*:pem_afc_test/0.*:PemAfcConfigValidation.*:PemAfcRtContract.*:closed_loop_test/0.*:burst_test/0.*:aec_test/0.*:AdaptationControlConfigValidation.*:nn_suppressor_test/0.*:NnSuppressorCrossPrecision.*:NnChainFloat32.*:FftEngineContract.*:PortableRandom.*:frequency_shifter_test/0.*:FrequencyShifterCrossPrecision.*:FrequencyShifterConfigValidation.*:FrequencyShifterRtContract.*:howl_detector_test/0.*:HowlDetectorCrossPrecision.*:HowlDetectorConfigValidation.*:HowlDetectorRtContract.*:howl_guard_test/0.*:HowlGuardConfigValidation.*:HowlGuardRtContract.*") + "fdaf_test/0.*:FdafCrossPrecision.*:FdafConfigValidation.*:FdafRtContract.*:fd_kalman_test/0.*:fd_kalman_test/1.*:kalman_loop_test/0.*:FdKalmanConfigValidation.*:FdKalmanRtContract.*:Levinson.*:LpcPredictor.*:SpeechPredictor.*:WarpedLpcPredictor.*:PredictorConfigValidation.*:pem_afc_test/0.*:PemAfcConfigValidation.*:PemAfcRtContract.*:closed_loop_test/0.*:burst_test/0.*:aec_test/0.*:AdaptationControlConfigValidation.*:nn_suppressor_test/0.*:NnSuppressorCrossPrecision.*:NnChainFloat32.*:FftEngineContract.*:PortableRandom.*:frequency_shifter_test/0.*:FrequencyShifterCrossPrecision.*:FrequencyShifterConfigValidation.*:FrequencyShifterRtContract.*:howl_detector_test/0.*:HowlDetectorCrossPrecision.*:HowlDetectorConfigValidation.*:HowlDetectorRtContract.*:howl_guard_test/0.*:HowlGuardConfigValidation.*:HowlGuardRtContract.*:reverb_mix_test/0.*:ReverbMixConfigValidation.*:ReverbMixRtContract.*") endif() gtest_discover_tests(mutap_tests DISCOVERY_TIMEOUT 120 ${mutap_emulated_selection} PROPERTIES TIMEOUT 900) diff --git a/tests/bare_metal_main.cpp b/tests/bare_metal_main.cpp index b31514b..0fa8201 100644 --- a/tests/bare_metal_main.cpp +++ b/tests/bare_metal_main.cpp @@ -54,17 +54,19 @@ int main() { "frequency_shifter_test/0.*:FrequencyShifterCrossPrecision.*:FrequencyShifterConfigValidation.*:" "FrequencyShifterRtContract.*:" "howl_detector_test/0.*:HowlDetectorCrossPrecision.*:HowlDetectorConfigValidation.*:HowlDetectorRtContract.*:" - "howl_guard_test/0.*:HowlGuardConfigValidation.*:HowlGuardRtContract.*"; + "howl_guard_test/0.*:HowlGuardConfigValidation.*:HowlGuardRtContract.*:" + "reverb_mix_test/0.*:ReverbMixConfigValidation.*:ReverbMixRtContract.*"; ::testing::InitGoogleTest(); const int rc = RUN_ALL_TESTS(); // A filter typo selects zero tests and RUN_ALL_TESTS() returns 0 — an // empty run must not pass green. Checked after the run because gtest // only applies the filter inside RUN_ALL_TESTS. The on-target selection - // is 101 tests (the FFT suites left with the FFT; FftEngineContract, the + // is 108 tests (the FFT suites left with the FFT; FftEngineContract, the // engine-contract rows MuTap owns, joined at the DspTap 0db95b6 bump; // PortableRandom, the fixtures' variates, when they replaced 's; // the frequency shifter's 7 with the shifter; the howl detector's 12 - // with the detector; the howl guard's 19 with the guard); 30 leaves headroom for + // with the detector; the howl guard's 19 with the guard; the reverb + // stage's 7 with the stage); 30 leaves headroom for // legitimate removals without masking a typo. const int selected = ::testing::UnitTest::GetInstance()->test_to_run_count(); if (selected < 30) { diff --git a/tests/support/decorrelated_loop.h b/tests/support/decorrelated_loop.h index 8d37b0d..8107a53 100644 --- a/tests/support/decorrelated_loop.h +++ b/tests/support/decorrelated_loop.h @@ -99,6 +99,40 @@ namespace mutap_test { tap::mu::allpass_hilbert m_pair; }; + /// A forward-path stage after the decorrelator: the chain's reverb slot + /// (support for test_reverb_stage*.cpp). Non-owning and type-erased, like + /// tap::mu::afc_stage_ref, plus a reset: every loop built with one calls + /// reset() in its constructor, so each bisection probe (a fresh loop with + /// a zero error history) also starts the stage from silence. The caller + /// owns the object and keeps it alive while any loop config refers to it; + /// one object per thread. + template + class forward_stage { + public: + forward_stage() = default; + + /// Any object with process_block(in, out, n) (in != out) and reset(). + template + static forward_stage of(Stage* stage) { + forward_stage f; + f.m_object = stage; + f.m_process = [](void* o, const Sample* in, Sample* out, size_t n) { + static_cast(o)->process_block(in, out, n); + }; + f.m_reset = [](void* o) { static_cast(o)->reset(); }; + return f; + } + + explicit operator bool() const { return m_process != nullptr; } + void process_block(const Sample* in, Sample* out, size_t n) const { m_process(m_object, in, out, n); } + void reset() const { m_reset(m_object); } + + private: + void* m_object = nullptr; + void (*m_process)(void*, const Sample*, Sample*, size_t) = nullptr; + void (*m_reset)(void*) = nullptr; + }; + /// Loop with a decorrelator and an auxiliary speaker feed. Same contract /// as closed_loop_sim: step() advances one block and returns the RMS of /// the loop output (+inf once the loop has gone non-finite). @@ -123,6 +157,10 @@ namespace mutap_test { /// short buffer loops. Null for none. const std::vector* aux = nullptr; double aux_gain = 1.0; + + /// Run on each block of the forward signal after the decorrelator + /// and before the gain (the chain's reverb slot); empty for none. + forward_stage stage; }; explicit decorrelated_loop(config cfg) @@ -132,7 +170,9 @@ namespace mutap_test { , m_u(m_cfg.block_size) , m_y(m_cfg.block_size) , m_e(m_cfg.block_size) - , m_aux(m_cfg.block_size, Sample(0)) { + , m_aux(m_cfg.block_size, Sample(0)) + , m_fwd(m_cfg.block_size, Sample(0)) + , m_staged(m_cfg.block_size, Sample(0)) { if (m_cfg.feedback_path.empty()) { throw std::invalid_argument("decorrelated_loop: empty feedback path"); } @@ -145,6 +185,9 @@ namespace mutap_test { throw std::invalid_argument("decorrelated_loop: forward_delay too short for this mode"); } set_forward_gain_db(m_cfg.forward_gain_db); + if (m_cfg.stage) { + m_cfg.stage.reset(); + } } void set_forward_gain_db(double k_db) { @@ -162,12 +205,22 @@ namespace mutap_test { const size_t lf = m_cfg.feedback_path.size(); const double lim = m_cfg.speaker_limit; + // The forward signal is rounded to Sample before the gain, as + // the loudspeaker path of a Sample-typed chain would be. Every + // sample it reads is at least one block old (forward_delay >= + // block_size), so the whole block exists before the stage runs. + for (size_t i = 0; i < b; ++i) { + m_fwd[i] = static_cast(forward_sample(m_t + static_cast(i))); + } + const Sample* fwd_block = m_fwd.data(); + if (m_cfg.stage) { + m_cfg.stage.process_block(m_fwd.data(), m_staged.data(), b); + fwd_block = m_staged.data(); + } for (size_t i = 0; i < b; ++i) { - const long long t = m_t + static_cast(i); - // The forward signal is rounded to Sample before the gain, - // as the loudspeaker path of a Sample-typed chain would be. - const auto fwd = static_cast(forward_sample(t)); - double u = m_gain * static_cast(fwd); + const long long t = m_t + static_cast(i); + const Sample fwd = fwd_block[i]; + double u = m_gain * static_cast(fwd); if (m_cfg.aux != nullptr && !m_cfg.aux->empty()) { const double a = m_cfg.aux_gain * static_cast((*m_cfg.aux)[static_cast(t) % m_cfg.aux->size()]); @@ -222,6 +275,10 @@ namespace mutap_test { const std::vector& mic_block() const { return m_y; } /// The aux feed's contribution to the last speaker block (zeros without one). const std::vector& aux_block() const { return m_aux; } + /// The forward signal of the last speaker block before the gain: after + /// the decorrelator and the stage (the chain's output, PROTOCOL.md + /// §7.3's "chain output"). + const std::vector& forward_block() const { return m_cfg.stage ? m_staged : m_fwd; } private: static constexpr size_t k_history = 8192; ///< power of two; >> forward_delay @@ -271,7 +328,9 @@ namespace mutap_test { std::vector m_u; std::vector m_y; std::vector m_e; - std::vector m_aux; ///< the aux feed as played, this block + std::vector m_aux; ///< the aux feed as played, this block + std::vector m_fwd; ///< the forward signal after the decorrelator, this block + std::vector m_staged; ///< ... and after the stage iir_ssb_shifter m_shifter; }; diff --git a/tests/support/karaoke_asg.h b/tests/support/karaoke_asg.h index 78a361e..4c52345 100644 --- a/tests/support/karaoke_asg.h +++ b/tests/support/karaoke_asg.h @@ -119,6 +119,10 @@ namespace mutap_test::karaoke { engine core = engine::kalman; material mat = material::held; bool aux = false; ///< white aux feed at the singer's level (0 dB) + /// After the decorrelator, before the gain: the chain's reverb slot + /// (decorrelated_loop.h's forward_stage; caller-owned, one per + /// thread). Empty for none. + forward_stage stage; }; /// How it is measured. @@ -147,6 +151,7 @@ namespace mutap_test::karaoke { cfg.forward_delay = delay; cfg.mode = s.mode; cfg.shift_hz = s.shift_hz; + cfg.stage = s.stage; if (s.aux && aux != nullptr) { cfg.aux = aux; cfg.aux_gain = 1.0; diff --git a/tests/support/reverb_rig.h b/tests/support/reverb_rig.h new file mode 100644 index 0000000..bd52e88 --- /dev/null +++ b/tests/support/reverb_rig.h @@ -0,0 +1,420 @@ +/// @file reverb_rig.h +/// @brief The vendored Dattorro plates behind tap::mu::reverb_mix / +/// shifted_dry_mix as a forward-path stage of decorrelated_loop, and +/// the reverb's own measurements (impulse response, T30, mix level). +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors +// +// Shared by test_reverb_stage.cpp (the gated rows) and +// test_reverb_stage_sweep.cpp (the MUTAP_SLOW sweep), so a number in +// docs/reverb-afc.md means the same measurement wherever it came from. +// Host-only: it includes the FAUST-generated classes (third_party/faust/). +// +// A rig is one forward path behind the canceller, in the chain's order +// (decorrelator -> reverb): +// +// reverb y = (1 - w) x + w r(x) reverb_mix +// bus_shift y = (1 - w) s(x) + w r(s(x)) frequency_shifter, then reverb_mix +// dry_shift y = (1 - w) s(x) + w r(x) shifted_dry_mix +// shift_only y = s(x) frequency_shifter +// +// r the plate's mono return (L or (L + R) / 2), s the library's SSB +// shifter (tap::mu::frequency_shifter). Every shift row of the +// reverb suites uses the library's shifter, shift_only included, so the +// three shift rows differ only in where the reverb sits. (The karaoke suites' +// shift rows use the loop's exact-ramp oscillator instead; the Hilbert pair +// is the same.) +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../fixtures/rir_cabin.h" +#include "../fixtures/rir_hall.h" +#include "../fixtures/rir_rehearsal.h" +#include "../fixtures/rir_studio.h" +#include "faust_generated.h" +#include "mutap/fft.h" +#include "mutap/frequency_shifter.h" +#include "mutap/reverb_stage.h" +#include "tap/dsp/math.h" + +namespace mutap_test::reverb { + + inline constexpr int k_fs = 48000; + inline constexpr size_t k_block = 64; + + /// Which vendored plate. + enum class variant { + shipped, ///< re.dattorro_rev as FAUST 2.88.0 ships it (dattorro.dsp) + paper ///< the paper's delay lengths, tank bug fixed (dattorro_paper.dsp) + }; + + enum class topology { reverb, bus_shift, dry_shift, shift_only }; + + inline const char* name(variant v) { + return v == variant::shipped ? "shipped" : "paper"; + } + inline const char* name(topology t) { + switch (t) { + case topology::bus_shift: + return "bus shift + reverb"; + case topology::dry_shift: + return "dry-only shift"; + case topology::shift_only: + return "shift only"; + case topology::reverb: + default: + return "reverb"; + } + } + + /// One reverb configuration. + struct params { + variant plate = variant::shipped; + double decay = 0.5; + double damping = 0.0005; + double wet = 0.30; + tap::mu::reverb_return ret = tap::mu::reverb_return::left; + topology topo = topology::reverb; + double shift_hz = 0.0; + }; + + /// FAUST's init() writes class-static state on some classes; construct + /// one plate at a time (faust_shim.h). + inline std::mutex& construction_mutex() { + static std::mutex m; + return m; + } + + /// A rig over one plate class. + class rig { + public: + explicit rig(const params& p) + : m_p(p) { + tap::mu::frequency_shifter::config sc; + sc.sample_rate = k_fs; + sc.shift_hz = p.shift_hz; + m_shifter = std::make_unique>(sc); + m_tmp.resize(k_block); + if (p.topo != topology::shift_only) { + const std::lock_guard lock(construction_mutex()); + if (p.plate == variant::shipped) { + m_impl = std::make_unique>(p, *m_shifter); + } + else { + m_impl = std::make_unique>(p, *m_shifter); + } + } + } + + void process_block(const double* in, double* out, size_t n) { + switch (m_p.topo) { + case topology::shift_only: + m_shifter->process_block(in, out, n); + return; + case topology::bus_shift: + m_shifter->process_block(in, m_tmp.data(), n); + m_impl->mix(m_tmp.data(), out, n); + return; + case topology::dry_shift: + m_impl->dry_mix(in, out, n); + return; + case topology::reverb: + default: + m_impl->mix(in, out, n); + return; + } + } + + /// Silence the plate and restart the shifter (phase 0). + void reset() { + if (m_impl) { + m_impl->clear(); + } + m_shifter->reset(); + } + + const params& parameters() const { return m_p; } + + private: + struct base { + base() = default; + base(const base&) = delete; + base& operator=(const base&) = delete; + virtual ~base() = default; + virtual void mix(const double* in, double* out, size_t n) = 0; + virtual void dry_mix(const double* in, double* out, size_t n) = 0; + virtual void clear() = 0; + }; + + template + struct impl final : base { + using block_t = mutap_faust::faust_block; + using mix_t = tap::mu::reverb_mix; + using dry_t = tap::mu::shifted_dry_mix, block_t>; + + impl(const params& p, tap::mu::frequency_shifter& shifter) + : plate(k_fs) + , cfg(make_config(p)) + , m(plate, cfg) + , d(shifter, plate, cfg) { + plate.set("decay", p.decay); + plate.set("damping", p.damping); + } + static typename mix_t::config make_config(const params& p) { + typename mix_t::config c; + c.block_size = k_block; + c.channels = 2; + c.wet = p.wet; + c.mode = p.ret; + return c; + } + void mix(const double* in, double* out, size_t n) override { m.process_block(in, out, n); } + void dry_mix(const double* in, double* out, size_t n) override { d.process_block(in, out, n); } + void clear() override { plate.clear(); } + + block_t plate; + typename mix_t::config cfg; + mix_t m; + dry_t d; + }; + + params m_p; + std::unique_ptr m_impl; + std::unique_ptr> m_shifter; + std::vector m_tmp; + }; + + /// The plate's impulse response, L and R, at 48 kHz: run for up to + /// `cap_s` seconds and truncate one past the last sample where max(|L|, + /// |R|) >= peak * 10^(-90/20) (phase 0's dattorro_ir.cpp rule). + struct impulse_response { + std::vector l; + std::vector r; + bool capped = false; ///< the tail had not fallen 90 dB by cap_s + }; + + inline impulse_response plate_ir(variant v, double decay, double damping, double cap_s = 20.0) { + auto run = [&](auto& plate) { + plate.set("decay", decay); + plate.set("damping", damping); + const auto n = static_cast(cap_s * k_fs); + impulse_response h; + h.l.assign(n, 0.0); + h.r.assign(n, 0.0); + std::vector in(k_block, 0.0); + for (size_t i = 0; i < n; i += k_block) { + const size_t b = std::min(k_block, n - i); + std::fill(in.begin(), in.end(), 0.0); + if (i == 0) { + in[0] = 1.0; + } + double* outs[2] = {&h.l[i], &h.r[i]}; + plate.process_mono(in.data(), outs, static_cast(b)); + } + double peak = 0.0; + for (size_t i = 0; i < n; ++i) { + peak = std::max({peak, std::abs(h.l[i]), std::abs(h.r[i])}); + } + const double thr = peak * std::pow(10.0, -90.0 / 20.0); + size_t len = 0; + for (size_t i = n; i-- > 0;) { + if (std::abs(h.l[i]) >= thr || std::abs(h.r[i]) >= thr) { + len = i + 1; + break; + } + } + h.capped = len == n; + h.l.resize(len); + h.r.resize(len); + return h; + }; + std::unique_lock lock(construction_mutex()); + if (v == variant::shipped) { + mutap_faust::faust_block plate(k_fs); + lock.unlock(); + return run(plate); + } + mutap_faust::faust_block plate(k_fs); + lock.unlock(); + return run(plate); + } + + /// The mono return of an impulse response: L or (L + R) / 2. + inline std::vector mono_return(const impulse_response& h, tap::mu::reverb_return ret) { + if (ret == tap::mu::reverb_return::left) { + return h.l; + } + std::vector m(h.l.size()); + for (size_t i = 0; i < m.size(); ++i) { + m[i] = 0.5 * (h.l[i] + h.r[i]); + } + return m; + } + + /// T30 by Schroeder backward integration: a least-squares line through + /// the energy decay curve from its first -5 dB to its first -35 dB + /// sample (test_faust_vendored.cpp's method). + inline double t30_seconds(const std::vector& h) { + std::vector edc(h.size()); + double acc = 0.0; + for (size_t i = h.size(); i-- > 0;) { + acc += h[i] * h[i]; + edc[i] = acc; + } + size_t lo = 0; + while (lo < h.size() && tap::dsp::power_db(edc[lo] / edc[0]) > -5.0) { + ++lo; + } + size_t hi = lo; + while (hi < h.size() && tap::dsp::power_db(edc[hi] / edc[0]) > -35.0) { + ++hi; + } + double sx = 0.0; + double sy = 0.0; + double sxx = 0.0; + double sxy = 0.0; + for (size_t i = lo; i <= hi && i < h.size(); ++i) { + const double t = static_cast(i) / k_fs; + const double y = tap::dsp::power_db(edc[i] / edc[0]); + sx += t; + sy += y; + sxx += t * t; + sxy += t * y; + } + const auto m = static_cast(hi - lo + 1); + const double slope = (m * sxy - sx * sy) / (m * sxx - sx * sx); // dB per second + return -60.0 / slope; + } + + /// The mix's broadband level change for white input, dB: 10 log10 of + /// the energy of (1 - w) delta + w r. Negative = the mix is quieter + /// than the dry bus (the crossfade gives up more level than the tail + /// adds back). + inline double mix_energy_db(const std::vector& r, double wet) { + double e = 0.0; + for (size_t i = 0; i < r.size(); ++i) { + double h = wet * r[i]; + if (i == 0) { + h += 1.0 - wet; + } + e += h * h; + } + return tap::dsp::power_db(e); + } + + // ------------------------------------------------- the bare-loop analysis + // + // Phase 0's method (_afc_poc/phase0-data/dattorro/dattorro_msg.py, + // reproduced): the loop L(w) = K e^{-jwd} H(w) F(w) with no canceller, + // H = (1 - w) + w R the mix, F a room fixture. Two limits per loop: + // bound = -20 log10 max |H F| (magnitude only; any delay) + // exact = -20 log10 of the largest positive real value of + // e^{-jwd} H F where its phase crosses 0 mod 2 pi (DC and + // Nyquist included when positive): the first K at which + // K L = 1, the Nyquist crossing (closed_loop.h's + // exact_msg_db, generalized to a mix) + // and the cost of the mix = dry limit - mix limit (positive = less + // stable gain). Spectra on an rfft grid with linear interpolation of + // the crossings between bins. + + /// X(w_k) = sum_n x[n] e^{-j w_k n}, k = 0 .. n/2, of x zero-padded to n + /// (a power of two). + inline std::vector> spectrum(const std::vector& x, size_t n) { + tap::mu::real_fft fft(n); + std::vector buf(n, 0.0); + for (size_t i = 0; i < x.size() && i < n; ++i) { + buf[i] = x[i]; + } + fft.forward_inplace(buf.data()); + // Ooura packing with the exp(+i) sign (closed_loop.h): X(w_k) = + // buf[2k] - j buf[2k+1]; DC and Nyquist are real, in buf[0], buf[1]. + std::vector> X(n / 2 + 1); + X[0] = {buf[0], 0.0}; + X[n / 2] = {buf[1], 0.0}; + for (size_t k = 1; k < n / 2; ++k) { + X[k] = {buf[2 * k], -buf[2 * k + 1]}; + } + return X; + } + + struct loop_limits { + double bound_db = 0.0; + double exact_db = 0.0; + }; + + /// The two limits of e^{-jwd} H F with H = (1 - wet) + wet R (R null: + /// H = 1, the dry loop), every spectrum on the same n-point grid. + inline loop_limits analytic_limits(const std::vector>& F, + const std::vector>* R, double wet, size_t d) { + const size_t half = F.size() - 1; + const size_t n = 2 * half; + auto at = [&](size_t k) { + std::complex h = F[k]; + if (R != nullptr) { + h *= (1.0 - wet) + wet * (*R)[k]; + } + const auto m = static_cast((static_cast(k) * d) % n); + const double ph = -2.0 * std::numbers::pi * static_cast(m) / static_cast(n); + return std::pair, double>{h * std::polar(1.0, ph), std::abs(h)}; + }; + double peak = 0.0; + auto [l0, a0] = at(0); + peak = a0; + double best = std::max(l0.real(), 0.0); + for (size_t k = 0; k < half; ++k) { + auto [l1, a1] = at(k + 1); + peak = std::max(peak, a1); + if (std::signbit(l0.imag()) != std::signbit(l1.imag())) { + const double t = l0.imag() / (l0.imag() - l1.imag()); + best = std::max(best, l0.real() + t * (l1.real() - l0.real())); + } + l0 = l1; + } + best = std::max(best, l0.real()); // Nyquist + return {-20.0 * std::log10(peak), -20.0 * std::log10(best)}; + } + + /// A room fixture as phase 0's bare-loop analysis loaded it: all 4096 + /// taps, as stored (not band-limited). "studio", "rehearsal", "hall"; + /// anything else is the cabin. + inline std::vector fixture_4096(const std::string& name) { + const float* f = fixtures::k_rir_cabin; + if (name == "studio") { + f = fixtures::k_rir_studio; + } + else if (name == "rehearsal") { + f = fixtures::k_rir_rehearsal; + } + else if (name == "hall") { + f = fixtures::k_rir_hall; + } + return {f, f + 4096}; + } + + /// Phase 0's median of the four fixtures: numpy's, the mean of the + /// middle two. + inline double median_of_four(std::vector v) { + std::sort(v.begin(), v.end()); + return 0.5 * (v[1] + v[2]); + } + + /// The smallest power of two >= x. + inline size_t next_pow2(size_t x) { + size_t n = 1; + while (n < x) { + n <<= 1U; + } + return n; + } + +} // namespace mutap_test::reverb diff --git a/tests/test_reverb_mix.cpp b/tests/test_reverb_mix.cpp new file mode 100644 index 0000000..2878cf2 --- /dev/null +++ b/tests/test_reverb_mix.cpp @@ -0,0 +1,310 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors +// +// The reverb-slot stages' plumbing (mutap/reverb_stage.h): the mix formula, +// both returns, the dry-only-shift topology, any block length, in place, +// config validation and the real-time contract. A small stateful stand-in +// replaces the reverb, so this file builds on every target and the emulated +// selections run its float suite (tests/bare_metal_main.cpp, TEST_FILTER in +// tests/CMakeLists.txt). The vendored Dattorro plates behind the same +// classes, and what they cost inside the loop, are test_reverb_stage.cpp +// (host-only). + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "mutap/afc_chain.h" +#include "mutap/frequency_shifter.h" +#include "mutap/reverb_stage.h" + +namespace { + + using tap::mu::reverb_return; + + /// A stand-in reverb with the FAUST shim's shape: mono in, `Channels` + /// out, state carried across calls. Channel c is a one-pole recursion + /// y = x + a_c y[n-1] with a different a per channel, so L and R differ + /// and a block split that lost state would show. + template + class stand_in_reverb { + public: + void process(const Sample* const* in, Sample* const* out, int n) noexcept { + for (int c = 0; c < Channels; ++c) { + const Sample a = Sample(0.5) + Sample(0.25) * static_cast(c); + for (int i = 0; i < n; ++i) { + auto& s = m_state[static_cast(c)]; + s = in[0][i] + a * s; + out[c][i] = s; + } + } + } + int num_inputs() const noexcept { return 1; } + int num_outputs() const noexcept { return Channels; } + void clear() noexcept { + for (auto& s : m_state) { + s = Sample(0); + } + } + + private: + std::array(Channels)> m_state{}; + }; + + /// The same reverb without num_outputs() / num_inputs(): the stages + /// must not require them. + template + struct bare_reverb { + void process(const Sample* const* in, Sample* const* out, int n) noexcept { + for (int i = 0; i < n; ++i) { + out[0][i] = Sample(2) * in[0][i]; + } + } + }; + + /// A dry-path stand-in: y = 3 x + the previous input (stateful). + template + struct stand_in_shifter { + void process_block(const Sample* in, Sample* out, size_t n) noexcept { + for (size_t i = 0; i < n; ++i) { + const Sample x = in[i]; + out[i] = Sample(3) * x + m_prev; + m_prev = x; + } + } + Sample m_prev = Sample(0); + }; + + template + std::vector ramp_signal(size_t n) { + std::vector x(n); + for (size_t i = 0; i < n; ++i) { + // Deterministic, sign-changing, no libm: a sawtooth of period 37. + x[i] = static_cast(static_cast(i % 37) / 18.0 - 1.0); + } + return x; + } + + /// Run a reference reverb over x in one call and return its channels. + template + std::vector> reference_channels(const std::vector& x) { + stand_in_reverb rev; + std::vector> ch(Channels, std::vector(x.size())); + std::vector outs(Channels); + for (int c = 0; c < Channels; ++c) { + outs[static_cast(c)] = ch[static_cast(c)].data(); + } + const Sample* const ins[1] = {x.data()}; + rev.process(ins, outs.data(), static_cast(x.size())); + return ch; + } + + template + class reverb_mix_test : public ::testing::Test {}; + using sample_types = ::testing::Types; + TYPED_TEST_SUITE(reverb_mix_test, sample_types); + + template + Sample tolerance() { + return std::is_same_v ? Sample(1e-5) : Sample(1e-12); + } + + // y = (1 - w) x + w L, and with the mid return w (L + R) / 2, across + // calls longer and shorter than the block (the stage works in pieces). + TYPED_TEST(reverb_mix_test, MixIsTheCrossfadeOfDryAndReturn) { + using sample = TypeParam; + const size_t n = 1000; + const auto x = ramp_signal(n); + const auto ref = reference_channels(x); + const sample w = sample(0.3); + const sample tol = tolerance(); + using reverb_t = stand_in_reverb; + using mix_t = tap::mu::reverb_mix; + const size_t calls[] = {64, 200, 7, 64, 1, 300, 364}; + for (const auto mode : {reverb_return::left, reverb_return::mid}) { + reverb_t rev; + typename mix_t::config cfg; + cfg.block_size = 64; + cfg.wet = w; + cfg.mode = mode; + mix_t mix(rev, cfg); + std::vector y(n); + size_t off = 0; + for (const size_t len : calls) { + mix.process_block(x.data() + off, y.data() + off, len); + off += len; + } + ASSERT_EQ(off, n); + size_t bad = 0; + for (size_t i = 0; i < n; ++i) { + const sample r = (mode == reverb_return::left) ? ref[0][i] : sample(0.5) * (ref[0][i] + ref[1][i]); + const sample expect = (sample(1) - w) * x[i] + w * r; + if (!(std::abs(y[i] - expect) <= tol * (sample(1) + std::abs(expect)))) { + ++bad; + } + } + EXPECT_EQ(bad, 0U) << (mode == reverb_return::left ? "L" : "(L + R) / 2"); + } + } + + // Wet 0 is the dry signal and wet 1 the return, sample for sample; a + // live set_wet() takes effect from the next block; in place works. + TYPED_TEST(reverb_mix_test, WetEndpointsAndLiveChanges) { + using sample = TypeParam; + using reverb_t = stand_in_reverb; + using mix_t = tap::mu::reverb_mix; + const size_t n = 256; + const auto x = ramp_signal(n); + const auto ref = reference_channels(x); + + reverb_t rev; + typename mix_t::config cfg; + cfg.wet = sample(0); + mix_t mix(rev, cfg); + std::vector y = x; // in place + mix.process_block(y.data(), y.data(), 64); + for (size_t i = 0; i < 64; ++i) { + EXPECT_EQ(y[i], x[i]) << i; + } + ASSERT_TRUE(mix.set_wet(sample(1))); + mix.process_block(y.data() + 64, y.data() + 64, 64); + for (size_t i = 64; i < 128; ++i) { + EXPECT_EQ(y[i], ref[0][i]) << i; + } + EXPECT_EQ(mix.wet(), sample(1)); + EXPECT_EQ(mix.channels(), 2U); + EXPECT_EQ(mix.block_size(), 64U); + EXPECT_EQ(mix_t::latency(), 0U); + } + + // The dry-only shift: y = (1 - w) s(x) + w r(x), the reverb fed the + // UNshifted input. + TYPED_TEST(reverb_mix_test, DryOnlyShiftFeedsTheReverbTheUnshiftedInput) { + using sample = TypeParam; + using reverb_t = stand_in_reverb; + using shift_t = stand_in_shifter; + using dry_t = tap::mu::shifted_dry_mix; + const size_t n = 500; + const auto x = ramp_signal(n); + const auto ref = reference_channels(x); + shift_t ref_shift; + std::vector s(n); + ref_shift.process_block(x.data(), s.data(), n); + + reverb_t rev; + shift_t sh; + typename dry_t::config cfg; + cfg.wet = sample(0.25); + cfg.mode = reverb_return::mid; + dry_t mix(sh, rev, cfg); + std::vector y(n); + mix.process_block(x.data(), y.data(), 130); + mix.process_block(x.data() + 130, y.data() + 130, n - 130); + const sample tol = tolerance(); + size_t bad = 0; + for (size_t i = 0; i < n; ++i) { + const sample r = sample(0.5) * (ref[0][i] + ref[1][i]); + const sample expect = sample(0.75) * s[i] + sample(0.25) * r; + if (!(std::abs(y[i] - expect) <= tol * (sample(1) + std::abs(expect)))) { + ++bad; + } + } + EXPECT_EQ(bad, 0U); + EXPECT_EQ(dry_t::latency(), 0U); + } + + // With the library's shifter at 0 Hz (its sample-exact bypass) and wet 0, + // the dry-only stage is the identity. + TYPED_TEST(reverb_mix_test, DryOnlyShiftAtZeroHertzAndZeroWetIsTheIdentity) { + using sample = TypeParam; + using reverb_t = stand_in_reverb; + using shift_t = tap::mu::frequency_shifter; + using dry_t = tap::mu::shifted_dry_mix; + typename shift_t::config sc; + sc.shift_hz = sample(0); + shift_t sh(sc); + reverb_t rev; + typename dry_t::config cfg; + cfg.wet = sample(0); + dry_t mix(sh, rev, cfg); + const auto x = ramp_signal(300); + std::vector y(x.size()); + mix.process_block(x.data(), y.data(), x.size()); + EXPECT_EQ(y, x); + } + + // Both stages fit afc_chain's reverb slot. + TYPED_TEST(reverb_mix_test, StagesFitTheChainSlot) { + using sample = TypeParam; + using reverb_t = stand_in_reverb; + static_assert(tap::mu::afc_stage, sample>); + static_assert( + tap::mu::afc_stage, reverb_t>, sample>); + static_assert(tap::mu::afc_stage>, sample>); + reverb_t rev; + tap::mu::reverb_mix mix(rev, {}); + tap::mu::afc_stage_ref ref(&mix); + EXPECT_TRUE(static_cast(ref)); + } + + TEST(ReverbMixConfigValidation, RejectsBadConfigs) { + using reverb_t = stand_in_reverb; + using mix_t = tap::mu::reverb_mix; + reverb_t rev; + auto make = [&](auto edit) { + mix_t::config c; + edit(c); + return mix_t(rev, c); + }; + EXPECT_THROW(make([](mix_t::config& c) { c.block_size = 0; }), std::invalid_argument); + EXPECT_THROW(make([](mix_t::config& c) { c.channels = 0; }), std::invalid_argument); + EXPECT_THROW(make([](mix_t::config& c) { c.channels = 3; }), std::invalid_argument); // reverb reports 2 + EXPECT_THROW(make([](mix_t::config& c) { c.wet = -0.1F; }), std::invalid_argument); + EXPECT_THROW(make([](mix_t::config& c) { c.wet = 1.5F; }), std::invalid_argument); + EXPECT_THROW(make([](mix_t::config& c) { c.wet = std::nanf(""); }), std::invalid_argument); + EXPECT_NO_THROW(make([](mix_t::config& c) { c.wet = 1.0F; })); + + // One channel: the (L + R) / 2 return is refused, at construction and live. + using mono_t = tap::mu::reverb_mix>; + stand_in_reverb mono; + mono_t::config mc; + mc.channels = 1; + mc.mode = reverb_return::mid; + EXPECT_THROW(mono_t(mono, mc), std::invalid_argument); + mc.mode = reverb_return::left; + mono_t m(mono, mc); + EXPECT_FALSE(m.set_return_mode(reverb_return::mid)); + EXPECT_EQ(m.return_mode(), reverb_return::left); + + // set_wet refuses what the constructor refuses, changing nothing. + EXPECT_FALSE(m.set_wet(-0.01F)); + EXPECT_FALSE(m.set_wet(1.01F)); + EXPECT_FALSE(m.set_wet(std::numeric_limits::infinity())); + EXPECT_EQ(m.wet(), 0.3F); + } + + TEST(ReverbMixRtContract, PostConstructionEntryPointsAreNoexcept) { + using rf = stand_in_reverb; + using mixf = tap::mu::reverb_mix; + using dryf = tap::mu::shifted_dry_mix, rf>; + static_assert(noexcept(std::declval().process_block(nullptr, nullptr, 0))); + static_assert(noexcept(std::declval().set_wet(0.0F))); + static_assert(noexcept(std::declval().set_return_mode(reverb_return::left))); + static_assert(noexcept(std::declval().wet())); + static_assert(noexcept(mixf::latency())); + static_assert(mixf::latency() == 0); + static_assert(noexcept(std::declval().process_block(nullptr, nullptr, 0))); + static_assert(noexcept(std::declval().set_wet(0.0F))); + static_assert(noexcept(std::declval().set_return_mode(reverb_return::mid))); + static_assert(dryf::latency() == 0); + SUCCEED(); + } + +} // namespace diff --git a/tests/test_reverb_stage.cpp b/tests/test_reverb_stage.cpp new file mode 100644 index 0000000..1ab4786 --- /dev/null +++ b/tests/test_reverb_stage.cpp @@ -0,0 +1,443 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors +// +// The reverb behind the canceller (docs/reverb-afc.md): mutap/reverb_stage.h +// over the vendored Dattorro plates (third_party/faust/), what the plates +// are (T30 per decay, the mix's level), what they cost a bare loop, and the +// gated in-loop rows. Host-only, like every FAUST and closed-loop claim +// (tests/CMakeLists.txt); the stages' plumbing on every target is +// test_reverb_mix.cpp. +// +// Measured on macOS 15.7 x86_64 (i9-8950HK), AppleClang 17, Release; each +// test's comment has its numbers. The in-loop rows bisect on 10 s probes, +// the length every held-note row of the MUTAP_SLOW sweep's +// ProbeConvergence converged at (test_reverb_stage_sweep.cpp), and bound +// medians over the five seed sets with several dB of margin for other +// hosts; single seeds move by up to 6.5 dB here. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "faust_generated.h" +#include "mutap/frequency_shifter.h" +#include "mutap/reverb_stage.h" +#include "support/karaoke_asg.h" +#include "support/reverb_rig.h" +#include "support/rooms.h" + +namespace { + + namespace rv = mutap_test::reverb; + namespace kk = mutap_test::karaoke; + using mutap_test::median; + using mutap_test::seed_in_set; + using tap::mu::reverb_return; + + /// At most 4 workers (the suites share their machines); + /// MUTAP_SLOW_THREADS overrides. + unsigned worker_count() { + if (const char* v = std::getenv("MUTAP_SLOW_THREADS")) { + return std::max(1U, static_cast(std::strtoul(v, nullptr, 10))); + } + const unsigned hw = std::thread::hardware_concurrency(); + return std::clamp(hw > 1 ? hw - 1 : 1U, 1U, 4U); + } + + void run_parallel(const std::vector>& jobs) { + std::atomic next{0}; + std::vector pool; + for (unsigned t = 0; t < worker_count(); ++t) { + pool.emplace_back([&] { + for (size_t j = next++; j < jobs.size(); j = next++) { + jobs[j](); + } + }); + } + for (auto& t : pool) { + t.join(); + } + } + + /// One in-loop configuration: the canceller alone (no plate) or with a + /// rig, and the chain bracket it is bisected over (dB re exact_msg_db). + struct loop_row { + const char* name; + bool with_rig = false; + rv::params p; + double lo = -15.0; + double hi = 30.0; + }; + + /// Per seed set: the dry open loop and each row's chain limit, on the + /// band-limited `room` at `delay`, `mat`, a `probe_s` probe; bisection + /// to 0.5 dB (support/karaoke_asg.h). The rows' ASG and cost per seed + /// are printed. + struct loop_result { + std::vector open; + std::vector> chain; ///< [row][seed set] + std::vector asg(size_t row) const { + std::vector a; + for (size_t s = 0; s < open.size(); ++s) { + a.push_back(chain[row][s] - open[s]); + } + return a; + } + /// chain[0] - chain[row]: the plate's cost against row 0 (the + /// canceller alone), per seed. + std::vector cost(size_t row) const { + std::vector c; + for (size_t s = 0; s < open.size(); ++s) { + c.push_back(chain[0][s] - chain[row][s]); + } + return c; + } + }; + + loop_result run_loop_rows(const std::string& room, size_t delay, kk::material mat, double probe_s, + const std::vector& rows) { + const auto path = kk::room(room); + const unsigned sets = mutap_test::k_claim_seed_sets; + loop_result r; + std::vector> jobs; + r.open.assign(sets, 0.0); + r.chain.assign(rows.size(), std::vector(sets, 0.0)); + for (unsigned set = 0; set < sets; ++set) { + const unsigned seed = seed_in_set(2, set); + jobs.emplace_back([&, set, seed] { + kk::protocol p; + p.delay = delay; + p.probe_blocks = kk::probe_blocks(probe_s); + r.open[set] = kk::open_loop_db(path, delay, kk::near_end(mat, p.probe_blocks, seed + 10), p); + }); + for (size_t k = 0; k < rows.size(); ++k) { + jobs.emplace_back([&, set, seed, k] { + kk::protocol p; + p.delay = delay; + p.probe_blocks = kk::probe_blocks(probe_s); + p.chain_lo = rows[k].lo; + p.chain_hi = rows[k].hi; + kk::setup s; + s.mat = mat; + if (rows[k].with_rig) { + rv::rig rig(rows[k].p); + s.stage = mutap_test::forward_stage::of(&rig); + r.chain[k][set] = kk::measure(path, s, p, seed, false).chain_db; + } + else { + r.chain[k][set] = kk::measure(path, s, p, seed, false).chain_db; + } + }); + } + } + run_parallel(jobs); + for (size_t k = 0; k < rows.size(); ++k) { + std::string line = std::string(rows[k].name) + " ASG"; + for (const double a : r.asg(k)) { + char b[16]; + std::snprintf(b, sizeof b, " %+.2f", a); + line += b; + } + char m[48]; + std::snprintf(m, sizeof m, " median %+.2f", median(r.asg(k))); + line += m; + if (k > 0) { + line += " | cost"; + for (const double c : r.cost(k)) { + char b[16]; + std::snprintf(b, sizeof b, " %+.2f", c); + line += b; + } + std::snprintf(m, sizeof m, " median %+.2f", median(r.cost(k))); + line += m; + } + std::printf("%s %s d=%zu: %s\n", room.c_str(), mat == kk::material::held ? "held" : "speech", delay, + line.c_str()); + } + return r; + } + + using shipped_t = mutap_faust::faust_block; + using paper_t = mutap_faust::faust_block; + + std::vector test_signal(size_t n) { + std::vector x(n, 0.0); + for (size_t i = 0; i < n; ++i) { + x[i] = std::sin(0.013 * static_cast(i)) * ((i % 480) < 240 ? 1.0 : 0.25); + } + return x; + } + + /// The plate's own outputs for x, run in 64-sample blocks. + template + std::vector> plate_outputs(Block& plate, const std::vector& x) { + std::vector> y(2, std::vector(x.size())); + for (size_t i = 0; i < x.size(); i += rv::k_block) { + double* outs[2] = {&y[0][i], &y[1][i]}; + plate.process_mono(&x[i], outs, static_cast(rv::k_block)); + } + return y; + } + +} // namespace + +// reverb_mix over the vendored plate is the crossfade of the input and the +// plate's own return, both returns, both plates. +TEST(ReverbStage, MixOverTheVendoredPlates) { + const size_t n = 64 * 200; + const auto x = test_signal(n); + for (const auto mode : {reverb_return::left, reverb_return::mid}) { + shipped_t ref_plate(rv::k_fs); + ref_plate.set("decay", 0.7); + const auto ref = plate_outputs(ref_plate, x); + + shipped_t plate(rv::k_fs); + plate.set("decay", 0.7); + tap::mu::reverb_mix::config cfg; + cfg.wet = 0.3; + cfg.mode = mode; + tap::mu::reverb_mix mix(plate, cfg); + std::vector y(n); + for (size_t i = 0; i < n; i += 128) { // two plate blocks per call + mix.process_block(&x[i], &y[i], 128); + } + double worst = 0.0; + for (size_t i = 0; i < n; ++i) { + const double r = (mode == reverb_return::left) ? ref[0][i] : 0.5 * (ref[0][i] + ref[1][i]); + worst = std::max(worst, std::abs(y[i] - (0.7 * x[i] + 0.3 * r))); + } + EXPECT_LT(worst, 1e-12) << (mode == reverb_return::left ? "L" : "M"); + } + // The paper plate fits the same stage, and a channel count it does not + // have is refused. + paper_t paper(rv::k_fs); + tap::mu::reverb_mix::config cfg; + EXPECT_NO_THROW((tap::mu::reverb_mix(paper, cfg))); + cfg.channels = 1; + EXPECT_THROW((tap::mu::reverb_mix(paper, cfg)), std::invalid_argument); +} + +// The dry-only shift feeds the plate the unshifted input: its wet part is +// the plate's return of x, not of shift(x). +TEST(ReverbStage, DryOnlyShiftKeepsTheTailOutOfTheShifter) { + const size_t n = 64 * 200; + const auto x = test_signal(n); + shipped_t ref_plate(rv::k_fs); + const auto ref = plate_outputs(ref_plate, x); + tap::mu::frequency_shifter::config sc; + sc.shift_hz = 5.0; + tap::mu::frequency_shifter ref_shift(sc); + std::vector s(n); + ref_shift.process_block(x.data(), s.data(), n); + + shipped_t plate(rv::k_fs); + tap::mu::frequency_shifter shift(sc); + using dry_t = tap::mu::shifted_dry_mix, shipped_t>; + dry_t::config cfg; + cfg.wet = 0.3; + dry_t mix(shift, plate, cfg); + std::vector y(n); + for (size_t i = 0; i < n; i += 64) { + mix.process_block(&x[i], &y[i], 64); + } + double worst = 0.0; + for (size_t i = 0; i < n; ++i) { + worst = std::max(worst, std::abs(y[i] - (0.7 * s[i] + 0.3 * ref[0][i]))); + } + EXPECT_LT(worst, 1e-12); +} + +// What the two plates are: T30 per decay (Schroeder, -5 to -35 dB) and the +// mix's broadband level for white input at each wet, damping 0.0005. +// Measured (macOS x86_64, AppleClang 17, Release; the impulse response is +// deterministic, so the margins absorb only libm / FMA differences): +// +// T30, L / (L + R) / 2, s decay 0.3 0.5 0.7 0.85 +// as shipped 1.0144 / 0.8620 1.1846 / 1.3082 2.1173 / 2.1530 4.6220 / 4.5825 +// paper lengths 1.7252 / 1.5975 2.1186 / 2.3626 3.7882 / 3.8598 8.2457 / 8.1939 +// +// mix level, L, dB w 0.15 w 0.30 w 0.50 +// decay 0.5 (both plates) -0.956 -1.687 -2.033 +// decay 0.7 (both plates) -0.925 -1.546 -1.622 +// decay 0.85 (both plates) -0.833 -1.134 -0.547 +// +// The two plates' mix levels agree to the printed digit at every decay +// (measured; not derived here). +TEST(ReverbStage, PlateDecayAndMixLevel) { + struct expect { + rv::variant v; + double decay; + double t30_l; + double t30_m; + }; + const expect rows[] = { + {rv::variant::shipped, 0.3, 1.0144, 0.8620}, {rv::variant::shipped, 0.5, 1.1846, 1.3082}, + {rv::variant::shipped, 0.7, 2.1173, 2.1530}, {rv::variant::shipped, 0.85, 4.6220, 4.5825}, + {rv::variant::paper, 0.3, 1.7252, 1.5975}, {rv::variant::paper, 0.5, 2.1186, 2.3626}, + {rv::variant::paper, 0.7, 3.7882, 3.8598}, {rv::variant::paper, 0.85, 8.2457, 8.1939}, + }; + for (const auto& e : rows) { + const auto h = rv::plate_ir(e.v, e.decay, 0.0005); + const double t_l = rv::t30_seconds(h.l); + const double t_m = rv::t30_seconds(rv::mono_return(h, reverb_return::mid)); + std::printf("%-7s decay %.2f: T30 L %.4f s, (L+R)/2 %.4f s, IR %zu samples%s; mix level L", rv::name(e.v), + e.decay, t_l, t_m, h.l.size(), h.capped ? " (capped)" : ""); + for (const double w : {0.15, 0.30, 0.50}) { + std::printf(" w%.2f %+.3f", w, rv::mix_energy_db(h.l, w)); + } + std::printf(" dB\n"); + EXPECT_FALSE(h.capped); + EXPECT_NEAR(t_l, e.t30_l, 0.01) << rv::name(e.v) << " decay " << e.decay; + EXPECT_NEAR(t_m, e.t30_m, 0.01) << rv::name(e.v) << " decay " << e.decay; + if (e.decay == 0.5) { + EXPECT_NEAR(rv::mix_energy_db(h.l, 0.30), -1.687, 0.005); + } + if (e.decay == 0.85) { + EXPECT_NEAR(rv::mix_energy_db(h.l, 0.30), -1.134, 0.005); + } + } +} + +// The bare-loop table for the paper plate (no canceller; the method of the +// anti-howl PoC's phase 0, test_reverb_stage_sweep.cpp's BareLoopCost): +// the mix's cost in open-loop stable gain by the magnitude bound and by the +// exact Nyquist crossing at d = 480, median over the four fixtures at 4096 +// taps (mean of the middle two), L return, damping 0.0005, w 0.30. The +// as-shipped row reproduces phase 0's published +0.12 / -0.60 dB, the +// method check. Measured (bound / exact, dB): +// +// as shipped, decay 0.5 (T30 1.18 s) +0.1245 / -0.6004 +// paper, decay 0.5 (T30 2.12 s) +0.2990 / +0.0625 +// paper, decay 0.7 (T30 3.79 s) +1.6114 / +1.1364 +// +// Deterministic: the impulse response and the FFTs repeat exactly on a +// host, so the 0.02 dB allowance is for other hosts' libm and FMA, not +// noise. (The grid itself is not that fine: over the sweep's 144 cells the +// cost moves up to 0.0309 dB by the bound and 0.0693 dB exact from the N +// grid to the 2N grid used here.) +TEST(ReverbStage, PaperPlateBareLoopCost) { + struct expect { + rv::variant v; + double decay; + double bound; + double exact; + }; + const expect rows[] = { + {rv::variant::shipped, 0.5, 0.1245, -0.6004}, + {rv::variant::paper, 0.5, 0.2990, 0.0625}, + {rv::variant::paper, 0.7, 1.6114, 1.1364}, + }; + const std::vector rooms = {"studio", "rehearsal", "hall", "cabin"}; + constexpr size_t k_d = 480; + for (const auto& e : rows) { + const auto h = rv::plate_ir(e.v, e.decay, 0.0005); + const size_t n = 2 * rv::next_pow2(4 * (h.l.size() + 4095 + k_d)); + const auto R = rv::spectrum(h.l, n); + std::vector bounds; + std::vector exacts; + for (const auto& room : rooms) { + const auto F = rv::spectrum(rv::fixture_4096(room), n); + const auto dry = rv::analytic_limits(F, nullptr, 0.0, k_d); + const auto mix = rv::analytic_limits(F, &R, 0.30, k_d); + bounds.push_back(dry.bound_db - mix.bound_db); + exacts.push_back(dry.exact_db - mix.exact_db); + } + const double b = rv::median_of_four(bounds); + const double x = rv::median_of_four(exacts); + std::printf("%-7s decay %.2f, w 0.30, L: bare-loop cost bound %+.4f dB, exact %+.4f dB\n", rv::name(e.v), + e.decay, b, x); + EXPECT_NEAR(b, e.bound, 0.02) << rv::name(e.v) << " decay " << e.decay; + EXPECT_NEAR(x, e.exact, 0.02) << rv::name(e.v) << " decay " << e.decay; + } +} + +// THE CANCELLER HOLDS BEHIND THE PLATES (runaway, held note, cabin, S1). +// Bisected like the sweep (support/karaoke_asg.h, the plate in the loop +// through tap::mu::reverb_mix) but to 0.5 dB over a narrowed bracket, on +// 10 s probes: every held-note row of ReverbStageSweep.ProbeConvergence +// converged at 10 s (its 20 to 160 s medians identical). The claim is a +// ceiling on the plate's cost, the canceller-alone limit minus the limit +// with the plate, as a median over the five seed sets, and a floor on the +// chain's ASG. Measured (macOS x86_64, AppleClang 17, Release; per seed set +// 2 / 22 / 42 / 62 / 82, dB): +// +// ASG cost +// canceller alone +7.62 +13.81 +11.56 +12.12 +10.44 (median +11.56) +// shipped d0.5 w0.15 +11.56 +8.19 +8.75 +12.41 +10.72 -3.94 +5.62 +2.81 -0.28 -0.28 median -0.28 +// shipped d0.5 w0.30 +11.84 +9.59 +11.00 +12.41 +9.59 -4.22 +4.22 +0.56 -0.28 +0.84 median +0.56 +// shipped d0.7 w0.30 +11.00 +10.72 +9.59 +11.84 +11.28 -3.37 +3.09 +1.97 +0.28 -0.84 median +0.28 +// shipped d0.85 w0.30 +14.09 +12.41 +10.72 +13.81 +13.81 -6.47 +1.41 +0.84 -1.69 -3.38 median -1.69 +// paper d0.5 w0.30 +11.00 +10.16 +11.56 +10.16 +12.97 -3.37 +3.66 +0.00 +1.97 -2.53 median +0.00 +// +// Gates: every plate's median cost < +3.5 dB (largest measured +0.56) and +// every row's median ASG > +7 dB (smallest +10.72). Single seeds swing by +// up to 6.5 dB either way, so no single-seed or direction claim is made +// (the sweep's six rooms carry the rest; docs/reverb-afc.md). 133.29 s on +// 4 threads of the Intel Mac in a full ctest run. +TEST(ReverbStage, CancellerHoldsBehindThePlates) { + std::vector rows = {{"canceller", false, {}, 2.0, 20.0}}; + rows.push_back({"shipped d0.5 w0.15", true, {}, 2.0, 20.0}); + rows.back().p.wet = 0.15; + rows.push_back({"shipped d0.5 w0.30", true, {}, 2.0, 20.0}); + rows.push_back({"shipped d0.7 w0.30", true, {}, 2.0, 20.0}); + rows.back().p.decay = 0.7; + rows.push_back({"shipped d0.85 w0.30", true, {}, 2.0, 20.0}); + rows.back().p.decay = 0.85; + rows.push_back({"paper d0.5 w0.30", true, {}, 2.0, 20.0}); + rows.back().p.plate = rv::variant::paper; + const auto r = run_loop_rows("cabin", kk::k_s1, kk::material::held, 10.0, rows); + for (size_t k = 0; k < rows.size(); ++k) { + EXPECT_GT(median(r.asg(k)), 7.0) << rows[k].name << " (measured >= +10.72)"; + if (k > 0) { + EXPECT_LT(median(r.cost(k)), 3.5) << rows[k].name << " (measured <= +0.56)"; + } + } +} + +// THE TWO SHIFT TOPOLOGIES BY RUNAWAY (held note, cabin, S1, 10 s probes, +// 0.5 dB): with the as-shipped plate at decay 0.5, wet 0.30 and a 2 Hz +// shift, the whole-bus shift (shifter in the decorrelator slot, the plate +// behind it) against the dry-only shift (shifted_dry_mix). Measured (macOS +// x86_64, AppleClang 17, Release; seed sets 2 / 22 / 42 / 62 / 82): +// +// whole-bus ASG +20.38 +18.12 +19.81 +17.00 +15.88 median +18.12 +// dry-only ASG +16.44 +15.03 +12.22 +16.16 +13.34 median +15.03 +// whole-bus - dry-only, per seed: +3.94 +3.09 +7.59 +0.84 +2.53, median +3.09 +// +// Gate: the per-seed median of whole-bus - dry-only > 0 (margin 3.09 dB). +// In the sweep (ReverbStageSweep.Shift, 20 s, six rooms, 2 and 5 Hz, decay +// 0.5 and 0.7) the per-room median is >= 0 in 23 of 24 at S1 (one -0.98). +// THIS IS RUNAWAY ONLY: on the ramp with the criterion the dry-only shift +// is the one audible later, in 13 of 16 pairs (docs/reverb-afc.md). 63.73 s +// on 4 threads of the Intel Mac in a full ctest run. +TEST(ReverbStage, WholeBusShiftHoldsMoreRunawayThanDryOnly) { + std::vector rows = {{"canceller", false, {}, 2.0, 20.0}}; + rows.push_back({"bus-shift 2 Hz d0.5 w0.30", true, {}, 8.0, 26.0}); + rows.back().p.topo = rv::topology::bus_shift; + rows.back().p.shift_hz = 2.0; + rows.push_back({"dry-shift 2 Hz d0.5 w0.30", true, {}, 8.0, 26.0}); + rows.back().p.topo = rv::topology::dry_shift; + rows.back().p.shift_hz = 2.0; + const auto r = run_loop_rows("cabin", kk::k_s1, kk::material::held, 10.0, rows); + std::vector diff; + for (size_t s = 0; s < r.open.size(); ++s) { + diff.push_back(r.chain[1][s] - r.chain[2][s]); + } + std::string line = "whole-bus - dry-only, per seed:"; + for (const double d : diff) { + char b[16]; + std::snprintf(b, sizeof b, " %+.2f", d); + line += b; + } + std::printf("%s median %+.2f\n", line.c_str(), median(diff)); + EXPECT_GT(median(diff), 0.0) << "measured +3.09"; +} diff --git a/tests/test_reverb_stage_sweep.cpp b/tests/test_reverb_stage_sweep.cpp new file mode 100644 index 0000000..db48287 --- /dev/null +++ b/tests/test_reverb_stage_sweep.cpp @@ -0,0 +1,909 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors +// +// The reverb-behind-the-canceller sweep behind test_reverb_stage.cpp and +// docs/reverb-afc.md. SLOW: skipped unless MUTAP_SLOW=1 (the repo's C++ +// slow-test convention; see tests/CMakeLists.txt). Jobs run on at most 4 +// threads unless MUTAP_SLOW_THREADS says otherwise; every table the doc +// quotes is printed, each finished bisection prints a "job" line (the +// checkpoint of a run that takes hours), and each row's median is recorded +// with RecordProperty (see --gtest_output=xml). +// +// MUTAP_SLOW=1 MUTAP_SLOW_THREADS=4 build/tests/mutap_tests --gtest_filter='ReverbStageSweep.*' +// +// Measurement (support/karaoke_asg.h, the karaoke suites' protocol, with +// the reverb as a forward-path stage of decorrelated_loop, support/ +// reverb_rig.h): band-limited 1024-tap rooms; PEM + FD-Kalman (speech +// cascade) converged 1500 blocks at exact_msg_db - 6 WITH the reverb in the +// loop; the forward gain bisected to 0.1 dB over [exact - 20, exact + 30] +// with the 40 dB runaway rule on a second seed's material (seed + 10); a +// fresh canceller copy and a silent reverb per probe; seeds 2 / 22 / 42 / +// 62 / 82 (rooms.h's five seed sets). Like for like: +// +// ASG = the chain's bisected limit - the dry open loop's (no canceller, no +// reverb), on the same probe material and length; +// cost = canceller-alone limit - canceller-plus-reverb limit, per seed +// (positive = the reverb costs stable gain); +// mix level = the mix's broadband level change for white input (10 log10 +// of the energy of (1 - w) delta + w r; reverb_rig.h), and +// matched-loudness cost = cost - mix level: the cost had the mix been +// made as loud as the dry bus. +// +// ProbeConvergence cabin and mt5 at S1, wet 0.30, probes 10 / 20 / 40 / +// 80 / 160 s: the as-shipped plate at decay 0.5 / 0.7 +// / 0.85 on the held note and the speech envelope; the +// paper plate at decay 0.5 on the held note; the 5 Hz +// shift with the as-shipped plate at decay 0.7, whole +// bus and dry-only, held note, to 80 s; the canceller +// alone to 80 s. "Converged" = the shortest probe from +// which every longer probe's five-seed median ASG is +// within 0.25 dB of the row's longest probe's. +// Grid six rooms (cabin, studio, rehearsal, hall; mt5, mt9). +// Held note, S1 and S3, 20 s probes: wet 0.15 / 0.30 / +// 0.50 x decay 0.5 / 0.7 / 0.85 on the as-shipped plate +// and at decay 0.5 on the paper plate (damping 0.0005, +// L return), damping 0.5 and the (L + R) / 2 return at +// wet 0.30 on the as-shipped plate. Speech envelope, S1, +// 80 s probes: the as-shipped plate at wet 0.30. Probes +// from ProbeConvergence (grid_probe). +// Shift the same rooms, held note, 20 s probes; 2 and 5 Hz: +// the shift alone, and with the as-shipped plate at +// decay 0.5 / 0.7, wet 0.30, the plate alone, the +// whole-bus shift + plate and the dry-only shift. +// Runaway limits. REVERB_SWEEP_DELAYS picks S1 / S3. +// BareLoopCost phase 0's analytic bare-loop table for both plates +// (no canceller): four fixtures at 4096 taps, d = 480, +// magnitude bound / exact Nyquist cost, and each +// plate's T30 per decay. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "support/karaoke_asg.h" +#include "support/reverb_rig.h" +#include "support/rooms.h" + +namespace { + + using mutap_test::median; + using mutap_test::seed_in_set; + namespace kk = mutap_test::karaoke; + namespace rv = mutap_test::reverb; + + bool slow_enabled() { + const char* v = std::getenv("MUTAP_SLOW"); + return v != nullptr && std::string(v) == "1"; + } + + /// At most 4 workers by default (the suites share their machines); + /// MUTAP_SLOW_THREADS overrides. + unsigned worker_count() { + if (const char* v = std::getenv("MUTAP_SLOW_THREADS")) { + return std::max(1U, static_cast(std::strtoul(v, nullptr, 10))); + } + const unsigned hw = std::thread::hardware_concurrency(); + return std::clamp(hw > 1 ? hw - 1 : 1U, 1U, 4U); + } + + std::string fmt(double v) { + char b[16]; + std::snprintf(b, sizeof b, "%+6.2f", v); + return b; + } + + const char* mat_name(kk::material m) { + return m == kk::material::held ? "held" : "speech"; + } + + /// A row's label: "canceller", or the reverb configuration. + std::string label(const rv::params* p) { + if (p == nullptr) { + return "canceller"; + } + char b[128]; + if (p->topo == rv::topology::shift_only) { + std::snprintf(b, sizeof b, "shift %g Hz", p->shift_hz); + return b; + } + std::snprintf(b, sizeof b, "%s d%.2f w%.2f m%g %s%s", rv::name(p->plate), p->decay, p->wet, p->damping, + p->ret == tap::mu::reverb_return::left ? "L" : "M", ""); + std::string s = b; + if (p->topo == rv::topology::bus_shift || p->topo == rv::topology::dry_shift) { + char t[64]; + std::snprintf(t, sizeof t, " %s %g Hz", p->topo == rv::topology::bus_shift ? "bus-shift" : "dry-shift", + p->shift_hz); + s += t; + } + return s; + } + + enum class kind { open, canceller, rig }; + + /// One bisection. + struct job { + std::string room; + size_t delay = kk::k_s1; + kk::material mat = kk::material::held; + double probe_s = 20.0; + unsigned seed = 2; + kind what = kind::canceller; + rv::params p; + double value = 0.0; ///< the bisected limit, dB + double exact = 0.0; + double seconds = 0.0; + }; + + std::string job_label(const job& j) { + return j.what == kind::rig ? label(&j.p) : (j.what == kind::open ? "open" : "canceller"); + } + + kk::protocol protocol_for(const job& j) { + kk::protocol p; + p.delay = j.delay; + p.probe_blocks = kk::probe_blocks(j.probe_s); + p.tol_db = 0.1; + p.chain_lo = -20.0; + p.chain_hi = 30.0; + return p; + } + + void run_one(job& j) { + const auto t0 = std::chrono::steady_clock::now(); + const auto path = kk::room(j.room); + kk::protocol p = protocol_for(j); + j.exact = mutap_test::exact_msg_db(path, j.delay); + kk::setup s; + s.mat = j.mat; + if (j.what == kind::open) { + const auto v = kk::near_end(j.mat, p.probe_blocks, j.seed + 10); + j.value = kk::open_loop_db(path, j.delay, v, p); + } + else if (j.what == kind::canceller) { + j.value = kk::measure(path, s, p, j.seed, false).chain_db; + } + else { + rv::rig rig(j.p); + s.stage = mutap_test::forward_stage::of(&rig); + j.value = kk::measure(path, s, p, j.seed, false).chain_db; + } + j.seconds = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); + } + + void run_jobs(std::vector& jobs) { + // Longest probes first, so the pool drains evenly. + std::stable_sort(jobs.begin(), jobs.end(), [](const job& a, const job& b) { return a.probe_s > b.probe_s; }); + std::atomic next{0}; + std::vector pool; + std::mutex io; + size_t done = 0; + const auto t0 = std::chrono::steady_clock::now(); + for (unsigned t = 0; t < worker_count(); ++t) { + pool.emplace_back([&] { + for (size_t i = next++; i < jobs.size(); i = next++) { + run_one(jobs[i]); + const std::lock_guard lock(io); + const job& j = jobs[i]; + std::printf(" job %s d=%zu %s %g s seed %u %s: %+.2f (%.0f s)\n", j.room.c_str(), j.delay, + mat_name(j.mat), j.probe_s, j.seed, job_label(j).c_str(), j.value - j.exact, j.seconds); + std::fflush(stdout); // the job lines are the checkpoint of a run that takes hours + if (++done % 50 == 0 || done == jobs.size()) { + std::printf(" ... %zu / %zu bisections, %.0f s\n", done, jobs.size(), + std::chrono::duration(std::chrono::steady_clock::now() - t0).count()); + std::fflush(stdout); + } + } + }); + } + for (auto& t : pool) { + t.join(); + } + double cpu = 0.0; + for (const auto& j : jobs) { + cpu += j.seconds; + } + std::printf("%zu bisections on %u threads: wall %.1f s, cpu-sum %.1f s\n", jobs.size(), worker_count(), + std::chrono::duration(std::chrono::steady_clock::now() - t0).count(), cpu); + } + + /// Row key: room, delay, material, probe, label. + using row_key = std::tuple; + + struct results { + std::map> chain; ///< per seed set, seed order + std::map> asg; + std::map> cost; ///< canceller - this, per seed (rig rows) + std::map exact; + }; + + /// Enqueue `what` for every seed set at (room, delay, mat, probe), plus + /// the open loop and the canceller alone there (deduplicated). + using ref_map = std::map, unsigned>; ///< 1 open, 2 canceller + + /// Enqueue `p` (null: nothing but the references) for every seed set at + /// (room, delay, mat, probe), plus the open loop there and, with + /// `canceller_ref`, the canceller alone (each reference once). + void enqueue(std::vector& jobs, ref_map& refs, const std::string& room, size_t delay, kk::material mat, + double probe_s, const rv::params* p, bool canceller_ref = true) { + const auto ref = std::make_tuple(room, delay, static_cast(mat), probe_s); + const unsigned have = refs[ref]; + const unsigned want = canceller_ref ? 3U : 1U; + for (unsigned set = 0; set < mutap_test::k_claim_seed_sets; ++set) { + job j; + j.room = room; + j.delay = delay; + j.mat = mat; + j.probe_s = probe_s; + j.seed = seed_in_set(2, set); + if (p != nullptr) { + j.what = kind::rig; + j.p = *p; + jobs.push_back(j); + } + if ((want & 1U) != 0 && (have & 1U) == 0) { + job o = j; + o.what = kind::open; + jobs.push_back(o); + } + if ((want & 2U) != 0 && (have & 2U) == 0) { + job c = j; + c.what = kind::canceller; + jobs.push_back(c); + } + } + refs[ref] = have | want; + } + + results collect(const std::vector& jobs) { + std::map, double> open; + std::map, double> canc; + for (const auto& j : jobs) { + const auto k = std::make_tuple(j.room, j.delay, static_cast(j.mat), j.probe_s, j.seed); + if (j.what == kind::open) { + open[k] = j.value; + } + else if (j.what == kind::canceller) { + canc[k] = j.value; + } + } + results r; + // Seed order: jobs were sorted by probe only (stable), so within a + // row the seed sets stay in enqueue order; sort explicitly anyway. + std::vector rows; + for (const auto& j : jobs) { + if (j.what != kind::open) { + rows.push_back(&j); + } + } + std::stable_sort(rows.begin(), rows.end(), [](const job* a, const job* b) { return a->seed < b->seed; }); + for (const job* j : rows) { + const row_key k{j->room, j->delay, static_cast(j->mat), j->probe_s, job_label(*j)}; + const auto s = std::make_tuple(j->room, j->delay, static_cast(j->mat), j->probe_s, j->seed); + r.chain[k].push_back(j->value); + r.asg[k].push_back(j->value - open.at(s)); + const auto c = canc.find(s); + if (j->what == kind::rig && c != canc.end()) { + r.cost[k].push_back(c->second - j->value); + } + r.exact[k] = j->exact; + } + return r; + } + + std::string per_seed(const std::vector& v) { + std::string s; + for (const double x : v) { + s += " " + fmt(x); + } + return s; + } + + void record(const char* suite, const row_key& k, const char* what, double v) { + char name[256]; + std::snprintf(name, sizeof name, "%s_%s_d%zu_%s_%gs_%s_%s", suite, std::get<0>(k).c_str(), std::get<1>(k), + mat_name(static_cast(std::get<2>(k))), std::get<3>(k), std::get<4>(k).c_str(), + what); + std::string key = name; + std::replace_if(key.begin(), key.end(), [](char c) { return c == ' ' || c == '+' || c == '.'; }, '_'); + ::testing::Test::RecordProperty(key, fmt(v)); + } + + /// Print every row: ASG per seed and median, and the cost per seed and + /// median for reverb rows. + void print_rows(const results& r, const char* suite) { + for (const auto& [k, v] : r.asg) { + std::printf("%-9s d=%-4zu %-6s %5.0f s %-44s ASG%s | median %s", std::get<0>(k).c_str(), std::get<1>(k), + mat_name(static_cast(std::get<2>(k))), std::get<3>(k), std::get<4>(k).c_str(), + per_seed(v).c_str(), fmt(median(v)).c_str()); + record(suite, k, "asg", median(v)); + const auto c = r.cost.find(k); + if (c != r.cost.end()) { + std::printf(" | cost%s | median %s", per_seed(c->second).c_str(), fmt(median(c->second)).c_str()); + record(suite, k, "cost", median(c->second)); + } + std::printf("\n"); + } + } + + rv::params plate(rv::variant v, double decay, double wet, double damping = 0.0005, + tap::mu::reverb_return ret = tap::mu::reverb_return::left) { + rv::params p; + p.plate = v; + p.decay = decay; + p.wet = wet; + p.damping = damping; + p.ret = ret; + return p; + } + + /// Probes from an environment override ("10,20,40") or the default. + std::vector probes_from_env(const char* var, std::vector dflt) { + const char* v = std::getenv(var); + if (v == nullptr) { + return dflt; + } + std::vector out; + std::string s = v; + size_t a = 0; + while (a < s.size()) { + const size_t b = s.find(',', a); + out.push_back(std::stod(s.substr(a, b == std::string::npos ? std::string::npos : b - a))); + if (b == std::string::npos) { + break; + } + a = b + 1; + } + return out; + } + +} // namespace + +// Probe-length convergence of the reverb rows: how long a probe the +// in-loop claims need. +TEST(ReverbStageSweep, ProbeConvergence) { + if (!slow_enabled()) { + GTEST_SKIP() << "set MUTAP_SLOW=1 for the reverb sweep"; + } + const auto probes = probes_from_env("REVERB_SWEEP_PROBES", {10.0, 20.0, 40.0, 80.0, 160.0}); + const auto both = std::vector{kk::material::held, kk::material::speech}; + const auto held = std::vector{kk::material::held}; + + /// One convergence row: a configuration, the materials and the longest + /// probe it runs to. + struct row { + rv::params p; + std::vector mats; + double longest = 160.0; + }; + // The as-shipped plate at the three decays carries the full probe set + // on both materials; the paper plate at decay 0.5 (the as-shipped + // decay-0.7 T30) on the held note only; the 5 Hz shift with the plate, + // whole-bus and dry-only, to 80 s on the held note (the budget). + std::vector rows; + for (const double decay : {0.5, 0.7, 0.85}) { + rows.push_back({plate(rv::variant::shipped, decay, 0.30), both, 160.0}); + } + rows.push_back({plate(rv::variant::paper, 0.5, 0.30), held, 160.0}); + for (const auto t : {rv::topology::bus_shift, rv::topology::dry_shift}) { + rv::params p = plate(rv::variant::shipped, 0.7, 0.30); + p.topo = t; + p.shift_hz = 5.0; + rows.push_back({p, held, 80.0}); + } + if (const char* only = std::getenv("REVERB_SWEEP_ROWS")) { // e.g. "0,1" (indices into rows) + const auto keep = probes_from_env("REVERB_SWEEP_ROWS", {}); + std::vector kept; + for (const double k : keep) { + kept.push_back(rows.at(static_cast(k))); + } + rows = kept; + (void)only; + } + const double canceller_longest = 80.0; // the canceller alone (the cost reference) up to 80 s + + std::vector jobs; + ref_map refs; + for (const std::string room : {"cabin", "mt5"}) { + for (const auto& rw : rows) { + for (const auto mat : rw.mats) { + for (const double probe : probes) { + if (probe <= rw.longest) { + enqueue(jobs, refs, room, kk::k_s1, mat, probe, &rw.p, probe <= canceller_longest); + } + } + } + } + } + run_jobs(jobs); + const auto r = collect(jobs); + std::printf("\nProbe convergence, S1, wet 0.30, damping 0.0005, L return; per seed set (2 22 42 62 82), dB\n"); + print_rows(r, "conv"); + + // The summary: each row's median ASG per probe, and the converged probe. + std::printf("\nMedian ASG by probe (dB); converged = the shortest probe from which every longer probe's median " + "is within 0.25 dB of the row's longest probe's\n"); + std::printf("%-9s %-6s %-44s", "room", "mat", "row"); + for (const double p : probes) { + std::printf(" %7.0f s", p); + } + std::printf(" converged | cost by probe\n"); + for (const std::string room : {"cabin", "mt5"}) { + for (const auto mat : both) { + std::vector> labels = {{"canceller", canceller_longest}}; + for (const auto& rw : rows) { + if (std::find(rw.mats.begin(), rw.mats.end(), mat) != rw.mats.end()) { + labels.emplace_back(label(&rw.p), rw.longest); + } + } + for (const auto& [l, longest] : labels) { + std::vector med; + std::vector cost; + for (const double p : probes) { + const row_key k{room, kk::k_s1, static_cast(mat), p, l}; + const auto a = r.asg.find(k); + if (p > longest || a == r.asg.end()) { + continue; + } + med.push_back(median(a->second)); + const auto c = r.cost.find(k); + cost.push_back(c == r.cost.end() ? std::nan("") : median(c->second)); + } + if (med.empty()) { + continue; + } + size_t conv = med.size(); // none + for (size_t i = med.size(); i-- > 0;) { + if (std::abs(med[i] - med.back()) <= 0.25) { + conv = i; + } + else { + break; + } + } + // The longest probe alone agrees with itself: converged + // only if the next-shorter one agrees too. + const bool ok = conv + 1 < med.size(); + std::printf("%-9s %-6s %-44s", room.c_str(), mat_name(mat), l.c_str()); + for (size_t i = 0; i < probes.size(); ++i) { + std::printf(" %9s", i < med.size() ? fmt(med[i]).c_str() : ""); + } + if (ok) { + std::printf(" %7.0f s |", probes[conv]); + } + else { + std::printf(" not conv. |"); + } + if (l != "canceller") { + for (const double c : cost) { + std::printf(" %s", std::isnan(c) ? " n/a" : fmt(c).c_str()); + } + } + std::printf("\n"); + char name[160]; + std::snprintf(name, sizeof name, "conv_%s_%s_%s_converged_s", room.c_str(), mat_name(mat), l.c_str()); + std::string key = name; + std::replace_if(key.begin(), key.end(), [](char ch) { return ch == ' ' || ch == '.'; }, '_'); + ::testing::Test::RecordProperty(key, ok ? fmt(probes[conv]) : std::string("none")); + } + } + } + for (const auto& [k, v] : r.asg) { + for (const double x : v) { + EXPECT_TRUE(std::isfinite(x)); + } + } +} + +namespace { + + const std::vector& sweep_rooms() { + static const std::vector k_rooms = {"cabin", "studio", "rehearsal", "hall", "mt5", "mt9"}; + return k_rooms; + } + + /// The mix level (dB, white input) of one plate configuration, cached. + double mix_level_db(const rv::params& p) { + static std::map, double> cache; + static std::mutex m; + const auto key = std::make_tuple(static_cast(p.plate), p.decay, p.damping, static_cast(p.ret), p.wet); + { + const std::lock_guard lock(m); + const auto it = cache.find(key); + if (it != cache.end()) { + return it->second; + } + } + const auto h = rv::plate_ir(p.plate, p.decay, p.damping); + const double v = rv::mix_energy_db(rv::mono_return(h, p.ret), p.wet); + const std::lock_guard lock(m); + cache[key] = v; + return v; + } + + /// Across rooms: the median of the per-room medians, their range, and + /// the per-seed range over every room. + struct across { + bool found = false; + double median_of_medians = 0.0; + double room_lo = 0.0; + double room_hi = 0.0; + double seed_lo = 0.0; + double seed_hi = 0.0; + }; + + across summarize(const std::map>& m, size_t delay, kk::material mat, double probe, + const std::string& l) { + std::vector meds; + std::vector all; + for (const auto& room : sweep_rooms()) { + const auto it = m.find(row_key{room, delay, static_cast(mat), probe, l}); + if (it == m.end()) { + continue; + } + meds.push_back(median(it->second)); + all.insert(all.end(), it->second.begin(), it->second.end()); + } + across a; + if (meds.empty()) { + return a; + } + a.found = true; + a.median_of_medians = median(meds); + a.room_lo = *std::min_element(meds.begin(), meds.end()); + a.room_hi = *std::max_element(meds.begin(), meds.end()); + a.seed_lo = *std::min_element(all.begin(), all.end()); + a.seed_hi = *std::max_element(all.begin(), all.end()); + return a; + } + + /// The probe a grid or shift row on `mat` runs with, from + /// ProbeConvergence (macOS x86_64): every held-note row there (both + /// plates, decay 0.5 to 0.85, the 5 Hz shift rows; cabin and mt5) + /// converged at 10 s, its 20 to 160 s medians identical, so the held + /// rows run 20 s, twice that; the speech rows converged at 40 or 80 s + /// where they converged at all (three of eight did not by 160 s), so + /// they run 80 s and are reported, not gated. + /// REVERB_SWEEP_PROBE overrides every row. + double grid_probe(kk::material mat) { + if (const char* e = std::getenv("REVERB_SWEEP_PROBE")) { + return std::stod(e); + } + return mat == kk::material::held ? 20.0 : 80.0; + } + + /// The decays the grid runs (REVERB_SWEEP_DECAYS, e.g. "0.5,0.7"). + std::vector grid_decays() { + return probes_from_env("REVERB_SWEEP_DECAYS", {0.5, 0.7, 0.85}); + } + + /// A delay and a material. + struct condition { + size_t delay; + kk::material mat; + }; + +} // namespace + +// ASG and the reverb's cost versus wet, decay, damping and return, behind +// the canceller, six rooms of both generator families. +TEST(ReverbStageSweep, Grid) { + if (!slow_enabled()) { + GTEST_SKIP() << "set MUTAP_SLOW=1 for the reverb sweep"; + } + const condition s1_held{kk::k_s1, kk::material::held}; + const condition s1_speech{kk::k_s1, kk::material::speech}; + const condition s3_held{kk::k_s3, kk::material::held}; + const condition s3_speech{kk::k_s3, kk::material::speech}; + const auto held_conds = std::vector{s1_held, s3_held}; + const auto all_conds = std::vector{s1_held, s3_held, s1_speech}; + + // The held note at S1 and S3: wet 0.15 / 0.30 / 0.50 x decay 0.5 / 0.7 + // / 0.85 on the as-shipped plate and at decay 0.5 on the paper plate + // (its convergence was measured there only), damping 0.0005, L return; + // damping 0.5 and the (L + R) / 2 return at wet 0.30 on the as-shipped + // plate. The speech envelope at S1, wet 0.30, the as-shipped plate, 80 s + // probes (grid_probe). Not run (the budget): the speech envelope at S3, + // and the paper plate at decay 0.7 and 0.85 (the bare loop covers them). + std::vector>> rows; + for (const auto v : {rv::variant::shipped, rv::variant::paper}) { + for (const double decay : grid_decays()) { + if (v == rv::variant::paper && decay > 0.6) { + continue; + } + for (const double wet : {0.15, 0.30, 0.50}) { + const bool speech = v == rv::variant::shipped && wet == 0.30; + rows.emplace_back(plate(v, decay, wet), speech ? all_conds : held_conds); + } + if (v == rv::variant::shipped) { + rows.emplace_back(plate(v, decay, 0.30, 0.5), held_conds); + rows.emplace_back(plate(v, decay, 0.30, 0.0005, tap::mu::reverb_return::mid), held_conds); + } + } + } + std::vector jobs; + ref_map refs; + for (const auto& room : sweep_rooms()) { + for (const auto& [c, conds] : rows) { + for (const auto& cond : conds) { + enqueue(jobs, refs, room, cond.delay, cond.mat, grid_probe(cond.mat), &c); + } + } + } + run_jobs(jobs); + const auto r = collect(jobs); + std::printf("\nGrid, band-limited rooms, per seed set (2 22 42 62 82), dB\n"); + print_rows(r, "grid"); + + std::printf("\nGrid summary: median over the six rooms of each room's five-seed median; [per-room median range]; " + "{per-seed range, 30 runs}; mix level = the mix's broadband level change (white input); matched = " + "cost - mix level\n"); + for (const auto& cond : {s1_held, s1_speech, s3_held, s3_speech}) { + if (!summarize(r.asg, cond.delay, cond.mat, grid_probe(cond.mat), "canceller").found) { + continue; // a condition this grid did not run + } + std::printf("\n%s, %s\n", cond.delay == kk::k_s1 ? "S1 (10 ms)" : "S3 (20 ms)", mat_name(cond.mat)); + for (const double probe : {10.0, 20.0, 40.0, 80.0, 160.0}) { + const auto c = summarize(r.asg, cond.delay, cond.mat, probe, "canceller"); + if (!c.found) { + continue; + } + std::printf(" canceller alone (%3.0f s probe): ASG %s [%s, %s] {%s, %s}\n", probe, + fmt(c.median_of_medians).c_str(), fmt(c.room_lo).c_str(), fmt(c.room_hi).c_str(), + fmt(c.seed_lo).c_str(), fmt(c.seed_hi).c_str()); + } + std::printf(" %-32s %5s %-26s %-40s %7s %8s\n", "row", "probe", "ASG [rooms]", "cost [rooms] {seeds}", + "mix lvl", "matched"); + for (const auto& [cfg, conds] : rows) { + const std::string l = label(&cfg); + const double probe = grid_probe(cond.mat); + const auto a = summarize(r.asg, cond.delay, cond.mat, probe, l); + const auto c = summarize(r.cost, cond.delay, cond.mat, probe, l); + if (!a.found) { + continue; + } + const double lvl = mix_level_db(cfg); + std::printf(" %-32s %3.0f s %s [%s, %s] %s [%s, %s] {%s, %s} %s %s\n", l.c_str(), probe, + fmt(a.median_of_medians).c_str(), fmt(a.room_lo).c_str(), fmt(a.room_hi).c_str(), + fmt(c.median_of_medians).c_str(), fmt(c.room_lo).c_str(), fmt(c.room_hi).c_str(), + fmt(c.seed_lo).c_str(), fmt(c.seed_hi).c_str(), fmt(lvl).c_str(), + fmt(c.median_of_medians - lvl).c_str()); + } + } + for (const auto& [k, v] : r.asg) { + for (const double x : v) { + EXPECT_TRUE(std::isfinite(x)); + } + } +} + +// The shift with the reverb: whole-bus shift + reverb (the chain's slot +// order, the tail recirculating through the shifter) against the dry-only +// shift (shifted_dry_mix), runaway limits. +TEST(ReverbStageSweep, Shift) { + if (!slow_enabled()) { + GTEST_SKIP() << "set MUTAP_SLOW=1 for the reverb sweep"; + } + const char* env = std::getenv("REVERB_SWEEP_SHIFT_PROBE"); + const double probe = env != nullptr ? std::stod(env) : grid_probe(kk::material::held); + std::vector configs; + for (const double hz : {2.0, 5.0}) { + rv::params s; + s.topo = rv::topology::shift_only; + s.shift_hz = hz; + configs.push_back(s); + } + // Wet 0.30 only (the budget; the audible driver covers 0.15 too). + for (const double decay : {0.5, 0.7}) { + for (const double wet : {0.30}) { + configs.push_back(plate(rv::variant::shipped, decay, wet)); + for (const double hz : {2.0, 5.0}) { + for (const auto t : {rv::topology::bus_shift, rv::topology::dry_shift}) { + rv::params p = plate(rv::variant::shipped, decay, wet); + p.topo = t; + p.shift_hz = hz; + configs.push_back(p); + } + } + } + } + // REVERB_SWEEP_DELAYS ("480", "960" or "480,960") runs one delay per + // invocation: the two halves' tables are independent. + std::vector delays; + for (const double d : probes_from_env("REVERB_SWEEP_DELAYS", {480.0, 960.0})) { + delays.push_back(static_cast(d)); + } + std::vector jobs; + ref_map refs; + for (const auto& room : sweep_rooms()) { + for (const size_t d : delays) { + for (const auto& c : configs) { + enqueue(jobs, refs, room, d, kk::material::held, probe, &c); + } + } + } + run_jobs(jobs); + const auto r = collect(jobs); + std::printf("\nShift, band-limited rooms, held note, %g s probes, per seed set (2 22 42 62 82), dB\n", probe); + print_rows(r, "shift"); + + std::printf("\nShift summary: median over the six rooms of each room's five-seed median [per-room range]; " + "bus - dry = per-seed (whole-bus shift + reverb) - (dry-only shift), chain limits\n"); + for (const size_t d : delays) { + std::printf("\n%s, held\n", d == kk::k_s1 ? "S1 (10 ms)" : "S3 (20 ms)"); + const auto c = summarize(r.asg, d, kk::material::held, probe, "canceller"); + std::printf(" %-44s ASG %s [%s, %s]\n", "canceller", fmt(c.median_of_medians).c_str(), fmt(c.room_lo).c_str(), + fmt(c.room_hi).c_str()); + for (const auto& cfg : configs) { + const std::string l = label(&cfg); + const auto a = summarize(r.asg, d, kk::material::held, probe, l); + std::printf(" %-44s ASG %s [%s, %s]", l.c_str(), fmt(a.median_of_medians).c_str(), fmt(a.room_lo).c_str(), + fmt(a.room_hi).c_str()); + if (cfg.topo == rv::topology::bus_shift) { + rv::params dry = cfg; + dry.topo = rv::topology::dry_shift; + std::vector room_meds; + for (const auto& room : sweep_rooms()) { + const auto& b = r.chain.at(row_key{room, d, 0, probe, l}); + const auto& o = r.chain.at(row_key{room, d, 0, probe, label(&dry)}); + std::vector diff; + for (size_t i = 0; i < b.size() && i < o.size(); ++i) { + diff.push_back(b[i] - o[i]); + } + room_meds.push_back(median(diff)); + } + std::printf(" bus - dry, per-room medians:"); + for (const double x : room_meds) { + std::printf(" %s", fmt(x).c_str()); + } + std::printf(" | median %s, %zu of %zu rooms < 0", fmt(median(room_meds)).c_str(), + static_cast( + std::count_if(room_meds.begin(), room_meds.end(), [](double x) { return x < 0.0; })), + room_meds.size()); + } + std::printf("\n"); + } + } + for (const auto& [k, v] : r.asg) { + for (const double x : v) { + EXPECT_TRUE(std::isfinite(x)); + } + } +} + +// Phase 0's bare-loop cost table (no canceller), for both plates: what the +// mix costs the open loop by the magnitude bound and by the exact Nyquist +// crossing at d = 480, four fixtures at 4096 taps, and each plate's T30. +TEST(ReverbStageSweep, BareLoopCost) { + if (!slow_enabled()) { + GTEST_SKIP() << "set MUTAP_SLOW=1 for the reverb sweep"; + } + const std::vector rooms = {"studio", "rehearsal", "hall", "cabin"}; + const std::vector decays = {0.3, 0.5, 0.7, 0.85}; + const std::vector damps = {0.0005, 0.2, 0.5}; + const std::vector wets = {0.15, 0.30, 0.50}; + constexpr size_t k_d = 480; + const auto t0 = std::chrono::steady_clock::now(); + + struct cell { + double bound = 0.0; + double exact = 0.0; + double dn_bound = 0.0; + double dn_exact = 0.0; + double level = 0.0; + }; + using cell_key = std::tuple; // plate, decay, damping, return, wet + std::map table; + std::map, double> t30; + std::map, size_t> ir_len; + + std::printf("\nBare-loop cost per room (dB, bound / exact at d = 480, on the 2N grid); rooms studio rehearsal " + "hall cabin\n"); + for (const auto v : {rv::variant::shipped, rv::variant::paper}) { + const int vi = static_cast(v); + for (const double decay : decays) { + for (const double damping : damps) { + const auto h = rv::plate_ir(v, decay, damping); + ir_len[{vi, decay, damping}] = h.l.size(); + const size_t n1 = rv::next_pow2(4 * (h.l.size() + 4095 + k_d)); + for (const auto mode : {tap::mu::reverb_return::left, tap::mu::reverb_return::mid}) { + const int mi = static_cast(mode); + const auto r = rv::mono_return(h, mode); + t30[{vi, decay, damping, mi}] = rv::t30_seconds(r); + // [grid][wet] -> per room (bound cost, exact cost) + std::map, std::vector>> costs; + for (const size_t n : {n1, 2 * n1}) { + const auto R = rv::spectrum(r, n); + for (const auto& room : rooms) { + const auto F = rv::spectrum(rv::fixture_4096(room), n); + const auto dry = rv::analytic_limits(F, nullptr, 0.0, k_d); + for (const double w : wets) { + const auto mix = rv::analytic_limits(F, &R, w, k_d); + costs[{n, w}].emplace_back(dry.bound_db - mix.bound_db, dry.exact_db - mix.exact_db); + } + } + } + for (const double w : wets) { + const auto& a = costs.at({n1, w}); + const auto& b = costs.at({2 * n1, w}); + std::vector bounds; + std::vector exacts; + cell c; + std::printf("%-7s decay %.2f damping %-6g %s w %.2f:", rv::name(v), decay, damping, + mode == tap::mu::reverb_return::left ? "L" : "M", w); + for (size_t i = 0; i < rooms.size(); ++i) { + bounds.push_back(b[i].first); + exacts.push_back(b[i].second); + c.dn_bound = std::max(c.dn_bound, std::abs(b[i].first - a[i].first)); + c.dn_exact = std::max(c.dn_exact, std::abs(b[i].second - a[i].second)); + std::printf(" %s / %s", fmt(b[i].first).c_str(), fmt(b[i].second).c_str()); + } + c.bound = rv::median_of_four(bounds); + c.exact = rv::median_of_four(exacts); + c.level = rv::mix_energy_db(r, w); + table[{vi, decay, damping, mi, w}] = c; + std::printf(" | median %s / %s, max |2N - N| %.4f / %.4f, mix level %s\n", + fmt(c.bound).c_str(), fmt(c.exact).c_str(), c.dn_bound, c.dn_exact, + fmt(c.level).c_str()); + } + } + } + } + } + + std::printf("\nT30 (s; Schroeder, -5 to -35 dB) and IR length to -90 dB re peak (samples)\n"); + std::printf("%-8s %6s %8s %10s %12s %10s\n", "plate", "decay", "damping", "T30 L", "T30 (L+R)/2", "IR len"); + for (const auto v : {rv::variant::shipped, rv::variant::paper}) { + const int vi = static_cast(v); + for (const double decay : decays) { + for (const double damping : damps) { + std::printf("%-8s %6.2f %8g %10.4f %12.4f %10zu\n", rv::name(v), decay, damping, + t30.at({vi, decay, damping, 0}), t30.at({vi, decay, damping, 1}), + ir_len.at({vi, decay, damping})); + } + } + } + for (const int mi : {0, 1}) { + std::printf("\nBare-loop cost, %s return: median over the four rooms (mean of the middle two, as phase 0), " + "bound / exact at d = 480; matched = the w 0.30 bound - the mix level at w 0.30\n", + mi == 0 ? "L" : "(L + R) / 2"); + std::printf("%-8s %6s %8s %7s | %-16s %-16s %-16s | %8s %8s\n", "plate", "decay", "damping", "T30", "w 0.15", + "w 0.30", "w 0.50", "lvl 0.30", "matched"); + for (const auto v : {rv::variant::shipped, rv::variant::paper}) { + const int vi = static_cast(v); + for (const double decay : decays) { + for (const double damping : damps) { + std::printf("%-8s %6.2f %8g %7.2f |", rv::name(v), decay, damping, + t30.at({vi, decay, damping, mi})); + for (const double w : wets) { + const auto& c = table.at({vi, decay, damping, mi, w}); + std::printf(" %s / %s ", fmt(c.bound).c_str(), fmt(c.exact).c_str()); + } + const auto& c = table.at({vi, decay, damping, mi, 0.30}); + std::printf("| %s %s\n", fmt(c.level).c_str(), fmt(c.bound - c.level).c_str()); + char name[96]; + std::snprintf(name, sizeof name, "bare_%s_decay%g_damping%g_%s_w030_bound", rv::name(v), decay, + damping, mi == 0 ? "L" : "M"); + std::string key = name; + std::replace(key.begin(), key.end(), '.', '_'); + ::testing::Test::RecordProperty(key, fmt(c.bound)); + } + } + } + } + std::printf("\nBareLoopCost wall %.1f s\n", + std::chrono::duration(std::chrono::steady_clock::now() - t0).count()); + for (const auto& [k, c] : table) { + EXPECT_TRUE(std::isfinite(c.bound) && std::isfinite(c.exact)); + } +} diff --git a/third_party/faust/README.md b/third_party/faust/README.md index d9d0991..4f06590 100644 --- a/third_party/faust/README.md +++ b/third_party/faust/README.md @@ -14,9 +14,15 @@ here: This is reference and test material. It is not part of MuTap's library: generated code never enters `include/mutap/`, nothing under `include/` -includes it, and it is compiled only by `tests/test_faust_vendored.cpp`, on -host builds (the cross and bare-metal test builds skip it; see -`tests/CMakeLists.txt`). Like every `third_party/` directory it is excluded +includes it, and it is compiled only by host-only tests +(`tests/test_faust_vendored.cpp`, `test_howl_guard_host.cpp`, +`test_reverb_stage.cpp` and its sweep) and the opt-in +`tools/notebook/karaoke_ramp_dump.cpp` (the cross and bare-metal test +builds skip it; see `tests/CMakeLists.txt`). The library drives the +Dattorro plates without knowing them: `tap::mu::reverb_mix` +(`include/mutap/reverb_stage.h`) takes any type with `faust_block`'s +`process(in, out, n)`; `docs/reverb-afc.md` measures them behind the +canceller. Like every `third_party/` directory it is excluded from the clang-format hook (`.pre-commit-config.yaml`: `^third_party/`) and from the clang-tidy sweep (`scripts/tidy.sh`: `'third_party' not in f`; `.clang-tidy`'s `HeaderFilterRegex` covers only `include/` and `tests/`). diff --git a/tools/notebook/CMakeLists.txt b/tools/notebook/CMakeLists.txt index 75eafb4..31d7d59 100644 --- a/tools/notebook/CMakeLists.txt +++ b/tools/notebook/CMakeLists.txt @@ -18,5 +18,8 @@ endif() if(MUTAP_BUILD_KARAOKE_DUMP) add_executable(karaoke_ramp_dump karaoke_ramp_dump.cpp) target_link_libraries(karaoke_ramp_dump PRIVATE MuTap::MuTap mutap_warnings) - target_include_directories(karaoke_ramp_dump PRIVATE ${CMAKE_SOURCE_DIR}/tests) + # tests/ for the suites' loop machinery; third_party/faust for the + # vendored Dattorro plates its --reverb rows put behind the canceller. + target_include_directories(karaoke_ramp_dump PRIVATE ${CMAKE_SOURCE_DIR}/tests + ${CMAKE_SOURCE_DIR}/third_party/faust) endif() diff --git a/tools/notebook/karaoke_ramp_dump.cpp b/tools/notebook/karaoke_ramp_dump.cpp index 30d1231..f6c3f28 100644 --- a/tools/notebook/karaoke_ramp_dump.cpp +++ b/tools/notebook/karaoke_ramp_dump.cpp @@ -29,11 +29,26 @@ // mic): the criterion's reference ramp. `--ir-only` writes PREFIX.ir.wav and // exits (the driver uses it at 4096 taps for the room's T30). // +// A REVERB BEHIND THE CANCELLER (docs/reverb-afc.md; tests/support/ +// reverb_rig.h, the reverb suites' forward path): `--reverb shipped|paper` +// puts a vendored Dattorro plate (FAUST, third_party/faust/) in the chain's +// reverb slot through tap::mu::reverb_mix, at `--wet`, `--decay`, +// `--damping` and `--return L|M` (L or (L + R) / 2). With a reverb, or with +// `--lib-shift`, the shift (`--shift-hz`) runs inside that stage on the +// library's frequency_shifter instead of the loop's exact-ramp oscillator: +// `--topology bus` shifts the whole bus before the reverb (the chain's slot +// order: the tail recirculates through the shifter), `--topology dry` uses +// tap::mu::shifted_dry_mix (only the dry path is shifted; the reverb is fed +// the unshifted bus). `--reverb-ir-only` writes the plate's impulse response +// (L, R) to PREFIX.reverb.wav and exits (the driver's --chain-rt60). +// // Outputs, with `--out PREFIX`: -// PREFIX.wav float32, 48 kHz, 3 channels: c (the canceller output), -// the voice stem (the near end v as it enters the mic) and +// PREFIX.wav float32, 48 kHz, 4 channels: c (the canceller output), +// the voice stem (the near end v as it enters the mic), // the track stem (the aux feed as the loudspeaker plays it; -// silence when there is none); +// silence when there is none) and the chain output (the +// forward signal after the shift and the reverb, before the +// gain: PROTOCOL.md 7.3's signal for a reverb-only chain); // PREFIX.gain.csv rows "time_s,gain_db", one per block (absolute bus gain); // PREFIX.ir.wav float32 mono, the band-limited room path in the loop (for // measure_rir.py's Schroeder T30, the criterion's --rt60); @@ -56,6 +71,7 @@ #include #include #include +#include #include #include @@ -67,6 +83,7 @@ #include "mutap/pem_afc.h" #include "support/closed_loop.h" #include "support/decorrelated_loop.h" +#include "support/reverb_rig.h" #include "support/rooms.h" namespace { @@ -94,8 +111,16 @@ namespace { double rate = 0.5; double max_over = 40.0; double tail_s = 2.0; - bool canceller = true; ///< false: the dry open loop (the reference ramp) - bool ir_only = false; ///< write PREFIX.ir.wav and stop + bool canceller = true; ///< false: the dry open loop (the reference ramp) + bool ir_only = false; ///< write PREFIX.ir.wav and stop + std::string reverb = "none"; ///< none | shipped | paper + double wet = 0.30; + double decay = 0.5; + double damping = 0.0005; + std::string ret = "L"; ///< L | M ((L + R) / 2) + std::string topology = "bus"; ///< bus | dry (with a shift and a reverb) + bool lib_shift = false; ///< shift inside the stage, no reverb + bool reverb_ir = false; ///< write PREFIX.reverb.wav and stop std::string out; }; @@ -105,7 +130,9 @@ namespace { "usage: karaoke_ramp_dump --out PREFIX [--room cabin|studio|rehearsal|hall|mtN]\n" " [--taps 1024] [--raw] [--delay 480] [--shift-hz 0] [--material held|ar]\n" " [--aux-db DB] [--seed 2] [--warmup 30] [--start-below 20] [--rate 0.5]\n" - " [--max-over 40] [--tail 2] [--no-canceller] [--ir-only]\n", + " [--max-over 40] [--tail 2] [--no-canceller] [--ir-only]\n" + " [--reverb none|shipped|paper] [--wet 0.3] [--decay 0.5] [--damping 0.0005]\n" + " [--return L|M] [--topology bus|dry] [--lib-shift] [--reverb-ir-only]\n", msg); std::exit(2); } @@ -166,6 +193,30 @@ namespace { else if (a == "--ir-only") { o.ir_only = true; } + else if (a == "--reverb") { + o.reverb = next(); + } + else if (a == "--wet") { + o.wet = std::stod(next()); + } + else if (a == "--decay") { + o.decay = std::stod(next()); + } + else if (a == "--damping") { + o.damping = std::stod(next()); + } + else if (a == "--return") { + o.ret = next(); + } + else if (a == "--topology") { + o.topology = next(); + } + else if (a == "--lib-shift") { + o.lib_shift = true; + } + else if (a == "--reverb-ir-only") { + o.reverb_ir = true; + } else if (a == "--out") { o.out = next(); } @@ -179,6 +230,18 @@ namespace { if (o.material != "held" && o.material != "ar") { usage("--material must be held or ar"); } + if (o.reverb != "none" && o.reverb != "shipped" && o.reverb != "paper") { + usage("--reverb must be none, shipped or paper"); + } + if (o.ret != "L" && o.ret != "M") { + usage("--return must be L or M"); + } + if (o.topology != "bus" && o.topology != "dry") { + usage("--topology must be bus or dry"); + } + if (o.reverb_ir && o.reverb == "none") { + usage("--reverb-ir-only needs --reverb"); + } return o; } @@ -263,6 +326,28 @@ namespace { return std::fclose(f) == 0; } + /// The reverb suites' forward path for these options (reverb_rig.h). + mutap_test::reverb::params rig_params(const options& o) { + namespace rv = mutap_test::reverb; + rv::params p; + p.plate = (o.reverb == "paper") ? rv::variant::paper : rv::variant::shipped; + p.decay = o.decay; + p.damping = o.damping; + p.wet = o.wet; + p.ret = (o.ret == "M") ? tap::mu::reverb_return::mid : tap::mu::reverb_return::left; + p.shift_hz = o.shift_hz; + if (o.reverb == "none") { + p.topo = rv::topology::shift_only; + } + else if (o.shift_hz != 0.0) { + p.topo = (o.topology == "dry") ? rv::topology::dry_shift : rv::topology::bus_shift; + } + else { + p.topo = rv::topology::reverb; + } + return p; + } + } // namespace int main(int argc, char** argv) { @@ -283,6 +368,20 @@ int main(int argc, char** argv) { } } + if (o.reverb_ir) { + const auto h = mutap_test::reverb::plate_ir(rig_params(o).plate, o.decay, o.damping); + std::vector> ch(2); + for (size_t i = 0; i < h.l.size(); ++i) { + ch[0].push_back(static_cast(h.l[i])); + ch[1].push_back(static_cast(h.r[i])); + } + if (!write_wav(o.out + ".reverb.wav", ch)) { + std::fprintf(stderr, "karaoke_ramp_dump: cannot write %s.reverb.wav\n", o.out.c_str()); + return 1; + } + return 0; + } + const double exact = mutap_test::exact_msg_db(path, o.delay); const double max_f = mutap_test::theoretical_msg_db(path); const double start_db = exact - o.start_below; @@ -300,11 +399,17 @@ int main(int argc, char** argv) { const std::vector aux = mutap_test::white_near_end(120000, o.seed + 777); decorrelated_loop::config lc; - lc.feedback_path = path; - lc.block_size = k_block; - lc.forward_delay = o.delay; - lc.forward_gain_db = start_db; - if (o.shift_hz != 0.0) { + lc.feedback_path = path; + lc.block_size = k_block; + lc.forward_delay = o.delay; + lc.forward_gain_db = start_db; + const bool staged = o.reverb != "none" || o.lib_shift; + std::unique_ptr rig; + if (staged) { + rig = std::make_unique(rig_params(o)); + lc.stage = mutap_test::forward_stage::of(rig.get()); + } + else if (o.shift_hz != 0.0) { lc.mode = forward_mode::shift; lc.shift_hz = o.shift_hz; } @@ -319,7 +424,7 @@ int main(int argc, char** argv) { pc.fdaf.partitions = o.taps / k_block; kalman_afc afc(pc); - std::vector> wav(3); + std::vector> wav(4); for (auto& ch : wav) { ch.reserve(n); } @@ -342,10 +447,12 @@ int main(int argc, char** argv) { const double rms = sim.step(&v[blk * k_block], o.canceller ? &afc : nullptr); const auto& e = sim.error_block(); const auto& a = sim.aux_block(); + const auto& fw = sim.forward_block(); for (size_t i = 0; i < k_block; ++i) { wav[0].push_back(static_cast(e[i])); wav[1].push_back(static_cast(v[blk * k_block + i])); wav[2].push_back(static_cast(a[i])); + wav[3].push_back(static_cast(fw[i])); } if (std::isnan(runaway_db) && rms >= 100.0) { runaway_db = g; @@ -367,11 +474,13 @@ int main(int argc, char** argv) { std::printf("{\"room\": \"%s\", \"taps\": %zu, \"banded\": %s, \"delay\": %zu, \"shift_hz\": %g, " "\"material\": \"%s\", \"aux_db\": %s, \"seed\": %u, \"canceller\": %s, \"exact_msg_db\": %.4f, " "\"theoretical_msg_db\": %.4f, \"start_db\": %.4f, \"warmup_s\": %g, \"rate_db_per_s\": %g, " - "\"runaway_db\": %s, \"runaway_t\": %s, \"seconds\": %.3f}\n", + "\"runaway_db\": %s, \"runaway_t\": %s, \"seconds\": %.3f, \"reverb\": \"%s\", \"wet\": %g, " + "\"decay\": %g, \"damping\": %g, \"return\": \"%s\", \"topology\": \"%s\", \"staged\": %s}\n", o.room.c_str(), o.taps, o.banded ? "true" : "false", o.delay, o.shift_hz, o.material.c_str(), o.aux ? std::to_string(o.aux_db).c_str() : "null", o.seed, o.canceller ? "true" : "false", exact, max_f, start_db, o.warmup_s, o.rate, std::isnan(runaway_db) ? "null" : std::to_string(runaway_db).c_str(), std::isnan(runaway_t) ? "null" : std::to_string(runaway_t).c_str(), - static_cast(blk * k_block) / k_fs); + static_cast(blk * k_block) / k_fs, o.reverb.c_str(), o.wet, o.decay, o.damping, o.ret.c_str(), + o.topology.c_str(), staged ? "true" : "false"); return 0; } diff --git a/tools/notebook/reverb_audible.py b/tools/notebook/reverb_audible.py new file mode 100644 index 0000000..8b33285 --- /dev/null +++ b/tools/notebook/reverb_audible.py @@ -0,0 +1,372 @@ +#!/usr/bin/env python3 +# SPDX-License-Identifier: MIT +# Copyright 2026 MuTap contributors +"""Audible limit of the simulated loop with a Dattorro plate behind the canceller. + +The reverb suites (tests/test_reverb_stage*.cpp) measure runaway limits by +bisection. This driver measures the AUDIBLE limit the way the karaoke one +(karaoke_audible.py) does, with the reverb in the chain: it runs +karaoke_ramp_dump with --reverb (the suites' forward path, tests/support/ +reverb_rig.h: tap::mu::reverb_mix over a vendored plate, the shift inside the +stage on the library's frequency_shifter) under the protocol's ramp (30 s +warm-up 20 dB under exact_msg_db, then 1 dB per 2 s, the live PEM + FD-Kalman +canceller adapting from a cold start), and applies tools/fixtures/ +howl_criterion.py by PROTOCOL.md section 7.3's table: + + dry, canceller c --rt60 ROOM + + shift c --rt60 ROOM --dechirp + + reverb only chain output --rt60 ROOM --chain-rt60 T (c as the cross-check) + + shift + reverb c --rt60 ROOM --chain-rt60 T --dechirp + +c is the canceller output (the dump's channel 0), the chain output its +channel 3 (after the shift and the reverb, before the gain). ROOM is the +band-limited cabin fixture's T30 at 4096 taps (karaoke_audible.py's +cabin_t30); T is the plate's own T30 for the row's decay, damping and return, +from the dump's --reverb-ir-only impulse response through measure_rir.py's +schroeder(). Analysing c with the plate's T30 in the hold is what counts a +reverb tail recirculating through the shifter: it is part of c's loop. + +Configurations (the product-chain candidates; --plate picks the plate): + dry no canceller (the audible open-loop reference) + plain the canceller alone + shift2/5 the canceller + a 2 / 5 Hz shift (library shifter, in the stage) + rev_D_W + the plate at decay D, wet W + bus2/5_D_W + shift on the whole bus, then the plate (the chain's slot order) + dry2/5_D_W + shifted_dry_mix: only the dry path shifted, the plate fed the unshifted bus +for D in --decays and W in --wets. + +Printed per configuration and delay: each seed's audible limit and the +ramp's runaway gain, both in dB over the dry loop's exact_msg_db, and their +medians; for reverb-only rows also the cross-check on c. --merge prints the +tables from earlier --json records. + +--control adds, for the shift-only rows, the same analysis of c with each +plate's T30 in the criterion's hold (--chain-rt60), the flags the product +rows get: without it the shift-only and the shift + plate rows are judged +by different holds. + +Disk and resumption: a run's dump is ~90 MB of float WAV (two minutes, four +channels) and its mono splits as much again, so each run's WAVs are deleted +once the criterion has read them (--keep-wav keeps them), and each run's +record is written beside it (PREFIX.record.json) as soon as it finishes; a +run whose record exists is not run again. --recover rebuilds the records of +runs that finished their analysis without one (the criterion's JSON and the +gain log are kept): exact_msg_db from a one-block dump with the run's +arguments, the runaway gain from the gain log by the dump's own rule (the +log stops 2 s after the runaway block, unless the ramp reached its end). + + cmake -S . -B build -DCMAKE_BUILD_TYPE=Release -DMUTAP_BUILD_KARAOKE_DUMP=ON + cmake --build build --target karaoke_ramp_dump + python3 tools/notebook/reverb_audible.py --dump build/tools/notebook/karaoke_ramp_dump --jobs 4 + +Requires numpy and scipy (the criterion's own requirements). +""" + +from __future__ import annotations + +import argparse +import concurrent.futures +import json +import pathlib +import statistics +import subprocess +import sys +import time + +import numpy as np +from scipy.io import wavfile + +HERE = pathlib.Path(__file__).resolve().parent +FIXTURES = HERE.parent / "fixtures" +sys.path.insert(0, str(FIXTURES)) +sys.path.insert(0, str(HERE)) +sys.dont_write_bytecode = True +import karaoke_audible # noqa: E402 (cabin_t30, plain_pass_limit, fmt) +import measure_rir # noqa: E402 + + +def configs(decays: list[float], wets: list[float], plate: str) -> dict: + """name -> (dump arguments, kind); kind picks section 7.3's row.""" + out = { + "dry": (["--no-canceller"], "plain"), + "plain": ([], "plain"), + "shift2": (["--lib-shift", "--shift-hz", "2"], "shift"), + "shift5": (["--lib-shift", "--shift-hz", "5"], "shift"), + } + for d in decays: + for w in wets: + rv = ["--reverb", plate, "--decay", f"{d:g}", "--wet", f"{w:g}"] + out[f"rev_{d:g}_{w:g}"] = (rv, "reverb") + for hz in (2, 5): + out[f"bus{hz}_{d:g}_{w:g}"] = (rv + ["--shift-hz", str(hz), "--topology", "bus"], "product") + out[f"dry{hz}_{d:g}_{w:g}"] = (rv + ["--shift-hz", str(hz), "--topology", "dry"], "product") + return out + + +def plate_t30(dump: str, work: pathlib.Path, plate: str, decay: float, damping: float, ret: str) -> float: + prefix = work / f"plate_{plate}_d{decay:g}_m{damping:g}" + subprocess.run([dump, "--out", str(prefix), "--reverb", plate, "--decay", f"{decay:g}", "--damping", + f"{damping:g}", "--reverb-ir-only"], check=True) + fs, ir = wavfile.read(str(prefix) + ".reverb.wav") + ir = ir.astype(np.float64) + mono = ir[:, 0] if ret == "L" else 0.5 * (ir[:, 0] + ir[:, 1]) + return float(measure_rir.schroeder(mono, fs)["T30"]) + + +def criterion(prefix: pathlib.Path, wav: str, programs: list[str], room_t30: float, chain_t30: float | None, + dechirp: bool, tag: str) -> dict: + cmd = [sys.executable, str(FIXTURES / "howl_criterion.py"), "analyze", wav, "--channel", "0"] + for p in programs: + cmd += ["--program", p] + cmd += ["--gain-log", str(prefix) + ".gain.csv", "--rt60", f"{room_t30:.6f}"] + if chain_t30 is not None: + cmd += ["--chain-rt60", f"{chain_t30:.6f}"] + if dechirp: + cmd += ["--dechirp"] + out = str(prefix) + f".{tag}.howl.json" + cmd += ["--json", out] + subprocess.run(cmd, check=True, capture_output=True, text=True) + return json.loads(pathlib.Path(out).read_text()) + + +# The dump's ramp defaults (karaoke_ramp_dump.cpp): the gain log has one row per block. +FS = 48000 +BLOCK = 64 +WARMUP_S, START_BELOW, MAX_OVER, RATE, TAIL_S = 30.0, 20.0, 40.0, 0.5, 2.0 + + +def run_args(name: str, spec: tuple, delay: int, seed: int) -> list: + return ["--delay", str(delay), "--seed", str(seed), "--material", "held", *spec[0]] + + +def prefix_of(work: pathlib.Path, name: str, delay: int, seed: int) -> pathlib.Path: + return work / f"held_{name}_d{delay}_s{seed}" + + +def runaway_from_gain_log(csv: str): + """The dump stops tail_s after the runaway block; a log that reaches max_blocks never ran away.""" + rows = np.loadtxt(csv, delimiter=",", skiprows=1) + max_blocks = int(WARMUP_S * FS / BLOCK) + int((START_BELOW + MAX_OVER) / RATE * FS / BLOCK) + 1 + tail = int(TAIL_S * FS / BLOCK) + return None if rows.shape[0] >= max_blocks else float(rows[rows.shape[0] - 1 - tail, 1]) + + +def one_run(dump: str, work: pathlib.Path, name: str, spec: tuple, delay: int, seed: int, room_t30: float, + chain_t30: dict, keep_wav: bool = False, control: bool = False) -> dict: + extra, kind = spec + t0 = time.perf_counter() + prefix = prefix_of(work, name, delay, seed) + record = pathlib.Path(str(prefix) + ".record.json") + if record.exists(): + return json.loads(record.read_text()) + out = subprocess.run([dump, "--out", str(prefix), *run_args(name, spec, delay, seed)], check=True, + capture_output=True, text=True) + meta = json.loads(out.stdout) + fs, x = wavfile.read(str(prefix) + ".wav") + voice = str(prefix) + ".voice.wav" + wavfile.write(voice, fs, x[:, 1].copy()) + t_chain = None + if meta.get("reverb", "none") != "none": + t_chain = chain_t30[(meta["reverb"], float(meta["decay"]), float(meta["damping"]), meta["return"])] + rec = dict(name=name, kind=kind, delay=delay, seed=seed, exact=meta["exact_msg_db"], + max_f=meta["theoretical_msg_db"], runaway=meta["runaway_db"], chain_t30=t_chain) + c_wav = str(prefix) + ".c.wav" + wavfile.write(c_wav, fs, x[:, 0].copy()) + if kind == "reverb": + chain_wav = str(prefix) + ".chain.wav" + wavfile.write(chain_wav, fs, x[:, 3].copy()) + h = criterion(prefix, chain_wav, [voice], room_t30, t_chain, False, "chain") + hc = criterion(prefix, c_wav, [voice], room_t30, t_chain, False, "c") + rec.update(limit=h.get("limit_db"), limit_c=hc.get("limit_db"), warnings=h.get("warnings", []), + warnings_c=hc.get("warnings", [])) + else: + h = criterion(prefix, c_wav, [voice], room_t30, t_chain, kind in ("shift", "product"), "c") + rec.update(limit=h.get("limit_db"), limit_plain=karaoke_audible.plain_pass_limit(h), + warnings=h.get("warnings", [])) + if kind == "shift" and control: + # The control: the same c analysed with each plate T30 in the hold, as the + # product rows are, so shift alone and shift + plate meet one criterion. + rec["limit_ctl"] = {} + for t in sorted(set(chain_t30.values())): + hc = criterion(prefix, c_wav, [voice], room_t30, t, True, f"ctl{t:.3f}") + rec["limit_ctl"][f"{t:.4f}"] = hc.get("limit_db") + rec["seconds"] = time.perf_counter() - t0 + if not keep_wav: + for w in pathlib.Path(work).glob(prefix.name + ".*wav"): + if not w.name.endswith(".ir.wav"): + w.unlink() + record.write_text(json.dumps(rec)) + return rec + + +def recover(dump: str, work: pathlib.Path, name: str, spec: tuple, delay: int, seed: int, chain_t30: dict): + """Rebuild the record of a run whose analysis finished without one (see the module docstring).""" + extra, kind = spec + prefix = prefix_of(work, name, delay, seed) + record = pathlib.Path(str(prefix) + ".record.json") + c_json = pathlib.Path(str(prefix) + ".c.howl.json") + if record.exists() or not c_json.exists() or (kind == "reverb" and + not pathlib.Path(str(prefix) + ".chain.howl.json").exists()): + return None + meta_prefix = work / f"meta_{name}_d{delay}_s{seed}" + out = subprocess.run([dump, "--out", str(meta_prefix), *run_args(name, spec, delay, seed), "--warmup", "0", + "--start-below", "0", "--max-over", "0", "--tail", "0"], check=True, capture_output=True, + text=True) + meta = json.loads(out.stdout) + for f in pathlib.Path(work).glob(meta_prefix.name + ".*"): + f.unlink() + t_chain = None + if meta.get("reverb", "none") != "none": + t_chain = chain_t30[(meta["reverb"], float(meta["decay"]), float(meta["damping"]), meta["return"])] + rec = dict(name=name, kind=kind, delay=delay, seed=seed, exact=meta["exact_msg_db"], + max_f=meta["theoretical_msg_db"], runaway=runaway_from_gain_log(str(prefix) + ".gain.csv"), + chain_t30=t_chain, recovered=True) + hc = json.loads(c_json.read_text()) + if kind == "reverb": + h = json.loads(pathlib.Path(str(prefix) + ".chain.howl.json").read_text()) + rec.update(limit=h.get("limit_db"), limit_c=hc.get("limit_db"), warnings=h.get("warnings", []), + warnings_c=hc.get("warnings", [])) + else: + rec.update(limit=hc.get("limit_db"), limit_plain=karaoke_audible.plain_pass_limit(hc), + warnings=hc.get("warnings", [])) + record.write_text(json.dumps(rec)) + return rec + + +def med(xs): + xs = [x for x in xs if x is not None] + return statistics.median(xs) if xs else None + + +def print_table(runs: list, order: list) -> None: + fmt = karaoke_audible.fmt + print("\nheld note, cabin (band-limited); dB over exact_msg_db (dry loop, phase-exact); audible = howl criterion " + "(PROTOCOL 7.3), runaway = the ramp's 40 dB rule; per seed in seed order") + for d in sorted({r["delay"] for r in runs}): + for n in order: + rs = sorted((r for r in runs if r["name"] == n and r["delay"] == d), key=lambda r: r["seed"]) + if not rs: + continue + aud = [None if r["limit"] is None else r["limit"] - r["exact"] for r in rs] + run = [None if r["runaway"] is None else r["runaway"] - r["exact"] for r in rs] + warn = sum(len(r["warnings"]) for r in rs) + t = rs[0].get("chain_t30") + print(f"d={d:4d} {n:13s} audible " + " ".join(fmt(a) for a in aud) + f" median {fmt(med(aud))} | " + "runaway " + " ".join(fmt(a) for a in run) + f" median {fmt(med(run))} | warnings {warn}" + + ("" if t is None else f" | chain T30 {t:.3f} s")) + if any("limit_ctl" in r for r in rs): + for t in sorted({t for r in rs for t in r.get("limit_ctl", {})}): + cc = [None if r.get("limit_ctl", {}).get(t) is None else r["limit_ctl"][t] - r["exact"] + for r in rs] + print(f"{'':20s}(control, c with --chain-rt60 {float(t):.3f}: " + " ".join(fmt(a) for a in cc) + + f" median {fmt(med(cc))})") + if rs[0]["kind"] == "reverb": + cc = [None if r.get("limit_c") is None else r["limit_c"] - r["exact"] for r in rs] + wc = sum(len(r.get("warnings_c", [])) for r in rs) + print(f"{'':20s}(c cross-check: " + " ".join(fmt(a) for a in cc) + + f" median {fmt(med(cc))}, warnings {wc})") + # audible - runaway per seed, where both exist + gap = [a - b for a, b in zip(aud, run) if a is not None and b is not None] + if gap: + print(f"{'':20s}(audible - runaway per seed: " + " ".join(fmt(g) for g in gap) + + f" median {fmt(statistics.median(gap))})") + + # The two shift topologies against each other: per seed, dry-only shift - whole-bus shift. + print("\ndry-only shift - whole-bus shift, per seed (same seed, delay, shift, decay and wet), dB") + pairs = 0 + higher = 0 + higher_run = 0 + for d in sorted({r["delay"] for r in runs}): + for n in order: + if not n.startswith("bus"): + continue + other = "dry" + n[3:] + b = {r["seed"]: r for r in runs if r["name"] == n and r["delay"] == d} + o = {r["seed"]: r for r in runs if r["name"] == other and r["delay"] == d} + seeds = sorted(set(b) & set(o)) + if not seeds: + continue + row = {} + for key in ("limit", "runaway"): + row[key] = [o[s][key] - b[s][key] for s in seeds if o[s][key] is not None and b[s][key] is not None] + print(f"d={d:4d} {other:13s} - {n:13s} audible " + " ".join(fmt(x) for x in row["limit"]) + + f" median {fmt(med(row['limit']))} | runaway " + " ".join(fmt(x) for x in row["runaway"]) + + f" median {fmt(med(row['runaway']))}") + if row["limit"] and row["runaway"]: + pairs += 1 + higher += med(row["limit"]) > 0 + higher_run += med(row["runaway"]) > 0 + print(f"dry-only shift higher (median over seeds): audible in {higher} of {pairs} pairs, runaway in " + f"{higher_run} of {pairs}") + + +def main(argv=None) -> None: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--dump", help="path to the karaoke_ramp_dump executable") + ap.add_argument("--work", default="reverb_audible_runs", help="directory for the runs (WAVs, JSON)") + ap.add_argument("--jobs", type=int, default=4) + ap.add_argument("--seeds", default="2,22,42,62,82") + ap.add_argument("--delays", default="480,960") + ap.add_argument("--plate", default="shipped", choices=["shipped", "paper"]) + ap.add_argument("--decays", default="0.5,0.7") + ap.add_argument("--wets", default="0.15,0.3") + ap.add_argument("--configs", default=None, help="comma-separated subset of the configuration names") + ap.add_argument("--json", default=None, help="write every run's record here") + ap.add_argument("--merge", nargs="+", default=None, metavar="JSON", help="print the tables from earlier records") + ap.add_argument("--keep-wav", action="store_true", help="keep each run's WAVs (about 180 MB a run)") + ap.add_argument("--control", action="store_true", + help="also analyse the shift-only rows' c with each plate T30 in the hold (--chain-rt60)") + ap.add_argument("--recover", action="store_true", + help="rebuild missing records from the work directory's criterion JSON and gain logs, then exit") + args = ap.parse_args(argv) + + decays = [float(x) for x in args.decays.split(",")] + wets = [float(x) for x in args.wets.split(",")] + cfgs = configs(decays, wets, args.plate) + order = list(cfgs) + if args.merge: + by_key = {} + for path in args.merge: + for r in json.loads(pathlib.Path(path).read_text()): + by_key[(r["name"], r["delay"], r["seed"])] = r + print_table(list(by_key.values()), order) + return + if not args.dump: + ap.error("--dump is required unless --merge is given") + names = args.configs.split(",") if args.configs else order + + work = pathlib.Path(args.work) + work.mkdir(parents=True, exist_ok=True) + t0 = time.perf_counter() + room_t30 = karaoke_audible.cabin_t30(args.dump, work) + print(f"band-limited cabin (4096 taps) T30 = {room_t30:.4f} s (--rt60)") + chain_t30 = {} + for d in decays: + key = (args.plate, d, 0.0005, "L") + chain_t30[key] = plate_t30(args.dump, work, args.plate, d, 0.0005, "L") + print(f"{args.plate} plate, decay {d:g}, damping 0.0005, L: T30 = {chain_t30[key]:.4f} s (--chain-rt60)") + + seeds = [int(s) for s in args.seeds.split(",")] + delays = [int(d) for d in args.delays.split(",")] + if args.recover: + got = [recover(args.dump, work, n, cfgs[n], d, s, chain_t30) for n in names for d in delays for s in seeds] + print(f"recovered {sum(r is not None for r in got)} records") + return + runs = [] + with concurrent.futures.ThreadPoolExecutor(max_workers=args.jobs) as pool: + futs = [pool.submit(one_run, args.dump, work, n, cfgs[n], d, s, room_t30, chain_t30, args.keep_wav, + args.control) + for n in names for d in delays for s in seeds] + for f in concurrent.futures.as_completed(futs): + runs.append(f.result()) + print_table(runs, order) + timed = [r for r in runs if not r.get("recovered")] + print(f"\n{len(runs)} runs ({len(runs) - len(timed)} from recovered records) on {args.jobs} jobs: wall " + f"{time.perf_counter() - t0:.1f} s, per-run sum of this invocation's and earlier timed runs " + f"{sum(r['seconds'] for r in timed):.1f} s") + if args.json: + pathlib.Path(args.json).write_text(json.dumps(runs, indent=1)) + + +if __name__ == "__main__": + main()