diff --git a/HANDOFF.md b/HANDOFF.md index b6b730b..f3394b3 100644 --- a/HANDOFF.md +++ b/HANDOFF.md @@ -136,16 +136,16 @@ carries the measured numbers; this is the map: carry the ABI tag: define it inside `namespace tap::mu::inline TAP_DSP_FFT_ABI`, or take the engine-dependent member as a template argument (as `aec_chain` does), and never forward-declare it in plain - `tap::mu`. That is transitive: it covers the five classes that hold a + `tap::mu`. That is transitive: it covers the six classes that hold a `basic_real_fft` by value (`partitioned_fdaf`, `partitioned_fdkf`, - `pem_afc`, `residual_suppressor`, `nn_suppressor`) and anything that + `pem_afc`, `residual_suppressor`, `nn_suppressor`, `spectral_reverb`) and anything that holds one of them at ``, or an `aec_chain`, by value — directly or through `optional` / `array` / `variant` — without naming it as a template argument. Double-only holders are exempt (the tag keys on the float default; `double` always runs the srdif engine, so the layout never changes): the C ABI handles, the ITU dump, and MuTap-Max's externals (`` behind `unique_ptr`), so MuTap-Max has nothing to - tag on its bump. `tests/test_fft_engine_contract.cpp` pins the five; it + tag on its bump. `tests/test_fft_engine_contract.cpp` pins the six; it cannot see a wrapper added later — review has to. (b) Pass every FFT size derived from configuration through `fft_detail::checked_fft_size` — CMSIS-DSP on the M55 supports only 32 … 4096 and an unchecked size @@ -343,38 +343,81 @@ carries the measured numbers; this is the map: - **Single-mic numbers here are not karaoke-suite numbers**: the chain's 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 - scratch harness (cabin band-limited, held note, wet 0.15, RT60 1 s, - PEM + FD-Kalman, like-for-like ASG, five seeds; not reproducible from - this repo until the class returns): - - **Shaped-from-F̂ loses to flat everywhere.** 40 s probes: S1 dry - +11.66 / flat +2.81 / shaped −27.66 (three seeds at the −30 bracket - floor); S3 +18.34 / +9.20 / +7.09; at the branch's 5.3 ms +5.45 / - +6.45 / +6.27. At S1 in studio / rehearsal / mt5 / mt9 (20 s) shaped - reads +0.92 / −1.60 / +6.21 / +3.28 against flat +7.77 / +6.48 / - +12.77 / +14.36. The branch's "shaped ≥ flat" came from 0.8 s probes - (at S1 then: shaped +12.29, flat +7.48, dry +10.41). - - **Reverb-with-canceller bisection had not converged by 40 s:** the - 20 s -> 40 s step still moved the medians by -25.73 (shaped, S1) to - +0.24 dB (flat, 5.3 ms), so no in-loop reverb claim was assertable. - - **Flat is unsafe open loop:** −18.05 to −24.61 dB against the dry - room (cabin S1 −24.61 at a 160 s probe; studio −18.05, rehearsal - −20.62, mt5 −20.86, mt9 −21.62 at 80 s). - - **Hypothesis, untested:** the shaping's per-bin allowance (up to - `shape_max` = 4, i.e. wet gain ×4 and decay time ×4 in the bins where - |F̂| is weakest) builds high-Q resonances where the canceller's - estimate is poorest, which the loop then finds. - - **Fix when it returns:** it holds a `basic_real_fft` by value - but is not defined inside `tap::mu::inline TAP_DSP_FFT_ABI` and does - not pass its FFT size through `fft_detail::checked_fft_size` (working - note 6 (a) and (b)); pin it in `test_fft_engine_contract.cpp`. +11. **`mutap::spectral_reverb`: RESOLVED, a measured negative result + (2026-10-08; phase 2 item 7, part 2; [`docs/reverb-afc.md`](docs/reverb-afc.md), + "Spectral reverb").** The per-bin reverb from the karaoke branch + (`8abe958`) is back as `include/mutap/spectral_reverb.h`, legal under + working note 6: defined inside `tap::mu::inline TAP_DSP_FFT_ABI`, its + 2 × block FFT size through `fft_detail::checked_fft_size`, pinned as the + sixth embedder in `test_fft_engine_contract.cpp`; an `afc_stage` with + `latency()` one block, a flat mode (`spectral_shaping::flat`) and an + explicit `shape_max`, shaped from F̂ (`reshape_from_impulse_response()`) + or the coherent bus sum |Σ F̂_m| (`reshape_from_impulse_responses()`). + Float plumbing `test_spectral_reverb.cpp` (both emulated selections; the + on-target count is now 118 in the comment, 120 selected by the filter on + the host build); host rows `test_spectral_reverb_host.cpp`; the + `MUTAP_SLOW` sweep `test_spectral_reverb_sweep.cpp` (ProbeConvergence, + Grid); `karaoke_ramp_dump --reverb spectral` and `reverb_audible.py + --plate spectral`. `decorrelated_loop`'s forward stage now takes a + stage's `latency()` out of the forward delay (0 for the plates: their + rows reproduce part 1's). Measured (macOS x86_64, AppleClang, Release, + five seed sets, held note unless stated): + - **Verdict: it loses to the plate everywhere measured.** Behind the + canceller (six rooms, S1 and S3, wet 0.15 / 0.30, rt60 1 s and the + plate's 1.1846 s T30, flat and shape_max 1 / 2 / 4, 40 s probes, NOT + converged and mostly understating the cost) it costs +2.83 to +30.57 + dB (median of room medians) where the as-shipped plate at decay 0.5 + costs −1.37 to +0.88; per seed the plate holds more in every room's + median for 22 of 24 rows, in 5 of 6 for the other two. At equal wet it + is 4.30 to 14.75 dB louder than the plate's mix; at matched loudness it + still costs more except shape_max 1 at S1, w 0.15 (−0.51 against the + plate's +1.44). Speech envelope (cabin and mt5, S1, w 0.30, 80 s): it + costs +4.30 to +17.09 dB where the plate raised the limit. + - **Convergence:** unlike the plates (10 s), only the cabin's flat, + shape_max 1 and 2 held rows converged (80 s); shape_max 4 and every mt5 + chain row had not by 160 s (cabin shape_max 4: −2.01 at 10 s, −19.02 at + 160 s). Gated rows use the converged ones only. + - **Flat is unsafe open loop, confirmed:** 17.15 to 31.80 dB below the + dry room (room medians, 160 s, both wets, S1 and S3); shaped rows + 12.56 to 29.84 dB below. Gated: on converged 80 s probes behind + `MUTAP_SLOW` (cabin / mt5) flat −30.51 / −22.77, shape_max 1 −19.61 / + −14.69, and the cost rows (cabin) +17.19 / +8.44; the default run + carries 20 s two-set proxies of both (conservative: short probes + understate), 104.5 s together on 3 threads with the level row. The + sanitizer legs could not afford the 80 s rows (697 s and 189 s). + - **Shaped vs flat:** shape_max 1 beats flat at S1 (+1.17 to +4.30 dB, + 3 to 5 of 6 rooms) and every ceiling does at S3, w 0.30 (+0.98 to + +3.71); shape_max 4 loses 7.91 to 12.79 dB at S1 in every room. The + branch's "shaped ≥ flat" (0.8 s probes) does not hold. + - **The hypothesis, measured:** at the moment of shaping, Spearman's ρ + between a bin's decay time and its misalignment |F − F̂| is negative in + every row (median −0.336 to −0.688): the long decays sit where the + absolute misalignment is SMALL, and the worst quarter of the bins gets + 0.022 to 0.178 of the tail energy (flat 0.246). Against the RELATIVE + misalignment |F − F̂|/|F| ρ is positive (+0.161 to +0.526). A + residual-loop indicator (max_k |F − F̂|·g_k against flat) rises with + the ceiling: −0.848 to −3.250 dB at shape_max 1, +5.487 to +9.717 at + shape_max 4. Per run, none of these predicts the run's shaped − flat + limit (ρ over 30 runs −0.351 to +0.699, no consistent sign). + - **Audible (cabin, ramp + criterion, PROTOCOL §7.3):** every + spectral row is audible earlier than the plate at the same wet (rt60 + 1.1846 s, five seeds; best shape_max 1 at S1, w 0.15, +11.90 against + the plate's +16.53), 14 of 16 below the dry loop's own audible limit on + the chain output (12 of 16 on c), and the ramp's runaway falls to + −6.94 … +18.31 where the plate's does not (+21.46 … +23.31). At S3 the + chain output is flagged near the ramp's start (−12.77 to −19.38) where + c reads −9.24 to +8.34; not examined: an open item of the same family + as part 1's S3 lines (the plates' ‡). 80 runs, 1380.5 s on 3 jobs, + WAVs deleted per run. + - **Kept, not removed:** the class stays as a measured negative result + (its header says so) because the gated rows and the sweep that hold + the verdict need it. Removing it again is the alternative; recommended + against for that reason. + - **Not run (the budget):** the six-room speech grid (speech is the two + ProbeConvergence rooms), converged probes for the grid (160 s would + have cost about four times the 40 s grid's 9.8 h), the backing track, + float32 in the loop, two microphones (the bus-sum shaping is tested for + plumbing only), audible rows beyond the cabin held note. 12. **The safety layer (2026-09-30; [`docs/howl-guard.md`](docs/howl-guard.md)).** PR A (tap/MuTap#77) adds `howl_detector`; PR B adds `include/mutap/howl_guard.h` (`tap::mu::howl_guard`: per-mic ARMING / @@ -434,7 +477,10 @@ carries the measured numbers; this is the map: 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). + item 7, part 1. Part 2 (2026-10-08, item 11) measured the spectral + reverb on this harness against these plates: it loses everywhere + measured and stays only as a negative result, so the plate is the + reverb this chain has; what follows is part 1's record. 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 @@ -492,6 +538,10 @@ carries the measured numbers; this is the map: 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. + - Part 2 changed one thing under part 1's rows: `decorrelated_loop`'s + forward stage now takes a stage's `latency()` out of the forward + delay. The plates report none, and their grid cost medians at S1 + repeat to the printed digit on part 2's 40 s probes. ## The next effort (Rev 4): AEC objects + echo chapter diff --git a/docs/reverb-afc.md b/docs/reverb-afc.md index c77ae39..a464f09 100644 --- a/docs/reverb-afc.md +++ b/docs/reverb-afc.md @@ -5,18 +5,20 @@ 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.* +number below comes from `tests/test_reverb_stage.cpp`, +`tests/test_spectral_reverb_host.cpp`, their `MUTAP_SLOW` sweeps +`tests/test_reverb_stage_sweep.cpp` and `tests/test_spectral_reverb_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). +probe lengths shown to converge first. The second half puts the spectral +reverb shaped from the canceller's estimate (`tap::mu::spectral_reverb`) +on the same harness beside the plate ([Spectral reverb](#spectral-reverb); +HANDOFF item 11): it loses. ## Headline @@ -42,6 +44,13 @@ on the same harness next (HANDOFF item 11). here. - **Probes: the held note converges at 10 s, decay 0.85 included**; the speech envelope drifts and is reported, not gated. +- **The spectral reverb loses to the plate everywhere measured** + ([Spectral reverb](#spectral-reverb)): behind the canceller it costs + +2.83 to +30.57 dB of runaway gain (six rooms, S1 and S3, wet 0.15 and + 0.30, flat and shaped; 40 s probes, not converged) where the plate costs + −1.37 to +0.88, and without a canceller the flat reverb sits 17 to 32 dB + below the dry room's limit. Shaping from F̂ helps only at shape_max 1. It + stays in the repository as a measured negative result. - **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 @@ -122,7 +131,10 @@ The karaoke suites' loop and protocol ([karaoke-afc.md](karaoke-afc.md), `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. + from silence, as it starts the loop's own history from zero. A stage's + own bulk latency is budgeted inside the forward delay, so the loop's total + delay stays S1's or S3's whatever the stage (the plates report none; the + spectral reverb's one block comes out of it, keeping it like for like). `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 @@ -518,6 +530,382 @@ decay and a wet: by the exact crossing from 0.73 dB below to 0.14 above. The sweep prints its full table. +## Spectral reverb + +Phase 2 item 7, part 2 (HANDOFF item 11): `tap::mu::spectral_reverb`, the +per-bin reverb from the karaoke branch (commit `8abe958`), back under the +ABI rules and measured on the harness above, beside the as-shipped plate. + +**Verdict: it loses to the plate, everywhere measured, and it stays in the +repository only as a measured negative result** (the header says so). Behind +the canceller it costs +2.83 to +30.57 dB of runaway gain (median over six +rooms, held note, 40 s probes) where the plate at the same wet costs −1.37 +to +0.88. Paired by seed, no room's median favours the spectral reverb +over the plate for 22 of the 24 spectral rows compared, and 1 of 6 rooms +does for the other two. Shaping from F̂ helps a flat spectral reverb at +shape_max 1 (and at S3, w 0.30, at every ceiling) but never closes the gap +to the plate, and shape_max 4 loses up to 12.79 dB more than flat at S1. +Without a canceller the flat reverb's limit is 17 to 32 dB below the dry +room's (room medians, 160 s probes), which confirms HANDOFF item 11's −18 +to −25. On the audible ramp every spectral row is audible earlier than the +plate at the same wet. The alternative, removing the class again, is +recommended against (HANDOFF item 11): the gated rows that hold the +negative result need it. + +### The spectral stage + +`include/mutap/spectral_reverb.h`: one complex one-pole per bin on +weighted-overlap-add spectra (Hann window of 2 × block, hop one block, no +synthesis window), S_k ← a_k S_k + X_k, Y_k = (1 − w) X_k + w_k S_k. The +whole output, dry path included, is the input one block late (`latency()`, +64 samples, 1.33 ms at 48 kHz). In the loop that block comes out of the +forward delay (`decorrelated_loop`'s forward stage now subtracts a stage's +`latency()`; 0 for the plates, whose rows are unchanged), so every +spectral row has the plates' total loop delay, 480 samples at S1 and 960 at +S3, as a chain's output delay line would absorb it. + +- **Flat** (`spectral_shaping::flat`): every bin at the configured rt60 and + wet; `reshape*()` is ignored. +- **Shaped** (`spectral_shaping::from_path`): per bin the allowance + (median|P| / |P_k|), clamped to [0.05, `shape_max`], scales the wet gain + and the decay time up where the path P is weak and down where it is + strong; the wet gains are renormalised to the flat shape's Σ w_k², the + decay times are not. P is the canceller's F̂ (`copy_impulse_response()`, + `reshape_from_impulse_response()`), or for several microphones the + coherent bus sum |Σ F̂_m| (`reshape_from_impulse_responses()`). +- **The rules:** defined inside `tap::mu::inline TAP_DSP_FFT_ABI` (it holds + a `basic_real_fft` by value) and pinned as the sixth embedder in + `tests/test_fft_engine_contract.cpp`; the transform size 2 × block goes + through `fft_detail::checked_fft_size` before any buffer is sized (block + 16 to 2048 on the Cortex-M55's CMSIS-DSP; an unchecked size there is a + release-mode HardFault). An `afc_stage` (`process_block(in, out, n)` + `noexcept`, n a multiple of the block), allocation only at construction, + float and double; `tests/test_spectral_reverb.cpp` runs its plumbing on + every target, both emulated selections included. + +What it is beside the plate (`SpectralReverbHost.DecayAndLevelBesideThePlate`): +its all-wet impulse response's T30 equals its rt60 (1.0000 s at rt60 1 s, +1.1846 s at 1.1846 s, the plate's decay-0.5 T30), and it is much louder +than the plate's mix at the same wet. For white input the flat stage reads ++3.40 / +7.60 dB at wet 0.15 / 0.30 (rt60 1 s) and +3.81 / +8.22 dB (rt60 +1.1846 s), against the plate's −0.956 / −1.687 dB: each bin's one-pole sums +about rt60 × 750 hops of its input. Shaped, the level rises with +`shape_max` (the decay times are not renormalised): the grid's medians run +from +3.34 (shape_max 1, w 0.15) to +13.06 dB (shape_max 4, w 0.30). + +### How the spectral rows are measured + +The plates' protocol ([How the loop is measured](#how-the-loop-is-measured)) +with one step the shaping needs (`tests/support/spectral_rig.h`): the +canceller converges with the reverb in the loop, flat (nothing to shape +from yet); the reverb is then shaped once from the converged F̂ and the +shape is held through every probe of the bisection, each probe starting the +tail from silence. That is the branch's protocol ("converge, then shape from +what the canceller has identified"). Shaped rows run without a canceller +(the open loop) are shaped from the F̂ a canceller converged to in the same +loop, then bisected with no canceller. The audible rows refresh the shape +from the live canceller every 256 blocks instead, since that canceller +starts cold. + +### Spectral probe convergence + +`SpectralReverbSweep.ProbeConvergence`: the cabin and mt5 at S1, wet 0.30, +rt60 1 s, flat and shaped (shape_max 1, 2, 4), held note to 160 s and +speech envelope to 80 s (the budget), each chain row with the open loop of +the same reverb beside it. Converged as for the plates. Median ASG, dB re +the dry open loop: + +| Room, material | Row (wet 0.30, rt60 1 s) | 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 | flat | −3.87 | −4.67 | −5.64 | −6.62 | −6.82 | 80 s | +| cabin, held | shape_max 1 | +3.36 | +2.27 | +1.48 | −0.37 | −0.27 | 80 s | +| cabin, held | shape_max 2 | +5.41 | +3.63 | +3.63 | +2.27 | +2.36 | 80 s | +| cabin, held | shape_max 4 | −2.01 | −3.30 | −5.74 | −17.46 | −19.02 | **not converged** | +| cabin, held | flat, no canceller | −29.14 | −29.75 | −30.14 | −30.33 | −30.43 | 80 s | +| cabin, held | shape_max 1, no canceller | −18.30 | −19.00 | −19.39 | −19.59 | −19.69 | 80 s | +| cabin, held | shape_max 2, no canceller | −18.50 | −19.30 | −19.79 | −20.08 | −20.18 | 80 s | +| cabin, held | shape_max 4, no canceller | −22.70 | −25.06 | −26.52 | −27.50 | −27.99 | **not converged** | +| mt5, held | canceller alone | +11.05 | +11.05 | +11.05 | +11.05 | | 10 s | +| mt5, held | flat | +1.21 | +1.00 | +0.41 | −0.27 | −0.57 | **not converged** | +| mt5, held | shape_max 1 | +10.57 | +8.03 | +9.49 | +8.03 | +8.42 | **not converged** | +| mt5, held | shape_max 2 | +8.52 | +5.00 | +5.98 | +6.86 | +6.46 | **not converged** | +| mt5, held | shape_max 4 | +5.98 | +2.75 | −13.85 | +2.07 | −20.00 | **not converged** | +| mt5, held | flat, no canceller | −21.91 | −22.23 | −22.52 | −22.81 | −22.91 | 80 s | +| mt5, held | shape_max 1, no canceller | −13.71 | −14.22 | −14.51 | −14.61 | −14.71 | 40 s | +| mt5, held | shape_max 2, no canceller | −16.35 | −17.73 | −18.61 | −19.10 | −19.39 | **not converged** | +| mt5, held | shape_max 4, no canceller | −19.86 | −22.03 | −23.40 | −24.18 | −24.57 | **not converged** | +| cabin, speech | canceller alone | +21.02 | +20.82 | +20.23 | +19.75 | | **not converged** | +| cabin, speech | flat | +3.26 | +3.05 | +2.95 | +2.66 | | **not converged** | +| cabin, speech | shape_max 1 | +15.47 | +15.55 | +15.45 | +15.45 | | 10 s | +| cabin, speech | shape_max 2 | +13.81 | +13.89 | +12.81 | +12.81 | | 40 s | +| cabin, speech | shape_max 4 | +12.42 | +10.76 | +11.25 | +10.86 | | **not converged** | +| mt5, speech | canceller alone | +23.26 | +23.16 | +22.87 | +22.68 | | 40 s | +| mt5, speech | flat | +10.18 | +11.25 | +9.30 | +8.81 | | **not converged** | +| mt5, speech | shape_max 1 | +14.28 | +14.38 | +13.59 | +13.40 | | 40 s | +| mt5, speech | shape_max 2 | +11.25 | +10.96 | +10.66 | +10.37 | | **not converged** | +| mt5, speech | shape_max 4 | +10.86 | +10.47 | +10.47 | +10.47 | | 20 s | + +- **Unlike the plates, the spectral rows do not converge at 10 s.** On the + held note the cabin's flat, shape_max 1 and 2 rows converged at 80 s; + shape_max 4 and every mt5 chain row had not converged by 160 s. The + medians mostly fall with probe length (more cost): the cabin's shape_max + 4 row goes from −2.01 at 10 s to −19.02 at 160 s, two of its five seeds + at the bracket floor (−20.00). This is HANDOFF item 11's "not converged by + 40 s", now measured to 160 s. +- **The open loop converges by 80 s** for flat and shape_max 1 (both rooms) + and shape_max 2 in the cabin; shape_max 4 drifts down to 160 s. +- **On speech** four of the eight spectral chain rows converged by 80 s + (cabin shape_max 1 at 10 s, shape_max 2 at 40 s; mt5 shape_max 1 at 40 s, + shape_max 4 at 20 s); the canceller alone in the cabin had not, as in the + plates' table. + +The grid below therefore runs the held note on 40 s probes (the budget: a +160 s row costs about four times a 40 s one), which is **not converged** and +is reported, not gated. The 40 s median ASG read at or above the 160 s one +in 7 of the 8 held rows above (mt5 shape_max 2: +5.98 against +6.46), so +the grid mostly understates the spectral reverb's cost. The open-loop rows +run on 160 s probes. The gated rows that use converged probes (80 s, the +cabin's flat and shape_max 1 chain rows, the cabin's and mt5's open loops) +run behind `MUTAP_SLOW`: at 697 s and 189 s on 3 threads they do not fit +the default run (the sanitizer legs). The default run carries a short proxy +for each, 20 s probes and two seed sets, and asserts the same directions. +Short probes understate both the spectral reverb's cost and the open +loop's danger (the table: every such median falls with probe length), so a +pass on 20 s is conservative. + +### Behind the canceller: the spectral reverb against the plate (runaway) + +`SpectralReverbSweep.Grid`, six rooms, five seed sets, held note, 40 s +probes. Each cell is the median over the six rooms of each room's five-seed +median, [the range of the room medians] and {the per-seed range over all +30 runs}. Cost = the canceller-alone limit − the limit with the reverb, per +seed; level = the stage's level change for white input (median over runs); +matched = cost − level. The plate rows are the as-shipped plate at decay +0.5 (T30 1.1846 s; rt60 1.1846 s is matched to it). Their cost medians +reproduce the plates' 20 s grid +([above](#what-the-plate-costs-behind-the-canceller-runaway)) to the +printed digit at S1 and at S3, w 0.15 (ASG medians within 0.08 dB); at S3, +w 0.30 the cost reads −1.37 against −1.17 there, the S3 canceller-alone +row moving between 20 and 40 s probes (+18.38 here against +18.40). + +**S1 (10 ms).** The canceller alone: ASG +11.74 [+11.05, +12.42] {+7.73, ++13.98}. + +| Reverb, wet | ASG | Cost [rooms] {seeds} | Level | Matched | +|---|---|---|---|---| +| plate 0.5, w 0.15 | +11.05 | +0.49 [−0.39, +2.05] {−4.00, +5.76} | −0.96 | +1.44 | +| spectral flat, rt60 1.18, w 0.15 | +7.34 | +5.66 [−1.76, +10.35] {−3.71, +13.48} | +3.81 | +1.85 | +| shape_max 1, rt60 1.18, w 0.15 | +8.61 | +2.83 [−2.15, +6.35] {−5.27, +10.55} | +3.34 | −0.51 | +| shape_max 2, rt60 1.18, w 0.15 | +6.46 | +7.03 [+0.98, +11.62] {−2.93, +15.43} | +4.91 | +2.12 | +| shape_max 4, rt60 1.18, w 0.15 | +0.21 | +14.26 [+3.61, +21.39] {−2.44, +31.74} | +7.25 | +7.01 | +| plate 0.5, w 0.30 | +11.54 | +0.88 [−0.78, +2.44] {−4.30, +4.20} | −1.69 | +2.57 | +| spectral flat, rt60 1, w 0.30 | +0.31 | +12.70 [+6.05, +17.38] {+4.39, +19.14} | +7.60 | +5.09 | +| shape_max 1, rt60 1, w 0.30 | +1.48 | +10.35 [+1.56, +12.01] {−2.54, +13.96} | +6.92 | +3.43 | +| shape_max 2, rt60 1, w 0.30 | +2.85 | +12.89 [+6.45, +14.45] {+2.73, +31.74} | +9.15 | +3.74 | +| shape_max 4, rt60 1, w 0.30 | −13.65 | +27.34 [+17.48, +32.42] {+5.37, +33.01} | +11.96 | +15.38 | +| spectral flat, rt60 1.18, w 0.30 | −0.76 | +13.09 [+9.57, +17.68] {+5.66, +21.97} | +8.22 | +4.87 | +| shape_max 1, rt60 1.18, w 0.30 | +4.22 | +10.74 [+3.91, +12.89] {+1.07, +18.46} | +7.61 | +3.13 | +| shape_max 2, rt60 1.18, w 0.30 | +0.02 | +13.48 [+5.08, +15.72] {+1.56, +30.66} | +9.94 | +3.53 | +| shape_max 4, rt60 1.18, w 0.30 | −18.34 | +30.57 [+13.38, +32.42] {+3.42, +33.98} | +12.22 | +18.34 | + +**S3 (20 ms).** The canceller alone: ASG +18.38 [+16.93, +18.59] {+16.52, ++19.36}. + +| Reverb, wet | ASG | Cost [rooms] {seeds} | Level | Matched | +|---|---|---|---|---| +| plate 0.5, w 0.15 | +18.87 | −0.59 [−1.76, +1.17] {−3.22, +3.71} | −0.96 | +0.37 | +| spectral flat, rt60 1.18, w 0.15 | +8.61 | +10.16 [+7.32, +11.91] {+5.86, +12.99} | +3.81 | +6.34 | +| shape_max 1, rt60 1.18, w 0.15 | +9.10 | +10.94 [+3.42, +12.79] {+1.46, +13.18} | +3.59 | +7.35 | +| shape_max 2, rt60 1.18, w 0.15 | +8.22 | +11.82 [+6.35, +14.84] {+2.44, +18.75} | +5.22 | +6.60 | +| shape_max 4, rt60 1.18, w 0.15 | +9.49 | +15.04 [+6.84, +17.77] {+1.46, +20.90} | +8.09 | +6.95 | +| plate 0.5, w 0.30 | +18.87 | −1.37 [−2.15, +0.49] {−4.20, +1.56} | −1.69 | +0.32 | +| spectral flat, rt60 1, w 0.30 | +0.61 | +18.07 [+13.48, +18.36] {+12.89, +20.02} | +7.60 | +10.46 | +| shape_max 1, rt60 1, w 0.30 | +3.14 | +16.99 [+10.35, +17.87] {+7.62, +20.61} | +6.94 | +10.05 | +| shape_max 2, rt60 1, w 0.30 | +1.68 | +17.77 [+12.01, +20.61] {+8.89, +23.44} | +9.04 | +8.74 | +| shape_max 4, rt60 1, w 0.30 | +1.78 | +18.07 [+12.50, +24.12] {+9.38, +27.15} | +12.19 | +5.87 | +| spectral flat, rt60 1.18, w 0.30 | −0.57 | +18.75 [+15.14, +19.92] {+14.65, +21.39} | +8.22 | +10.53 | +| shape_max 1, rt60 1.18, w 0.30 | +3.44 | +15.14 [+10.16, +17.77] {+7.23, +20.41} | +7.49 | +7.65 | +| shape_max 2, rt60 1.18, w 0.30 | +1.78 | +16.70 [+11.91, +18.65] {+7.81, +24.80} | +9.85 | +6.85 | +| shape_max 4, rt60 1.18, w 0.30 | +3.63 | +15.23 [+11.33, +17.09] {+3.91, +30.08} | +13.06 | +2.17 | + +The same runs paired by seed: the spectral chain's limit minus the plate's +at the same wet, and shaped minus flat at the same rt60 and wet; the +median over the six rooms of each room's per-seed median, [range], (rooms +where it is > 0): + +| | S1, w 0.15 (rt60 1.18) | S1, w 0.30 (rt60 1) | S1, w 0.30 (rt60 1.18) | S3, w 0.15 (rt60 1.18) | S3, w 0.30 (rt60 1) | S3, w 0.30 (rt60 1.18) | +|---|---|---|---|---|---|---| +| flat − plate | −3.52 (1 of 6) | −11.52 (0) | −11.72 (0) | −10.25 (0) | −17.77 (0) | −19.43 (0) | +| shape_max 1 − plate | −2.73 (1 of 6) | −8.50 (0) | −7.03 (0) | −8.20 (0) | −16.21 (0) | −15.43 (0) | +| shape_max 2 − plate | −4.59 (0) | −7.81 (0) | −10.45 (0) | −9.08 (0) | −16.89 (0) | −17.19 (0) | +| shape_max 4 − plate | −10.64 (0) | −24.51 (0) | −29.30 (0) | −8.98 (0) | −17.58 (0) | −15.72 (0) | +| shape_max 1 − flat | +1.17 (3 of 6) | +1.86 (4 of 6) | +4.30 (5 of 6) | +1.46 (3 of 6) | +1.95 (6 of 6) | +3.71 (6 of 6) | +| shape_max 2 − flat | −1.17 (2 of 6) | −0.49 (2 of 6) | −0.78 (2 of 6) | +0.59 (3 of 6) | +1.37 (4 of 6) | +2.64 (6 of 6) | +| shape_max 4 − flat | −7.91 (0) | −12.79 (0) | −12.60 (0) | +0.39 (3 of 6) | +0.98 (3 of 6) | +3.71 (6 of 6) | + +- **The plate holds more runaway gain than every spectral row:** the + median per-seed difference is −2.73 to −29.30 dB, and in 22 of the 24 + comparisons no room's median favours the spectral reverb (1 of 6 rooms + for flat and shape_max 1 at S1, w 0.15). Its cost, +2.83 to +30.57 dB, + starts above the plate's largest (+0.88 on this grid; +1.46 on the + plates' own). +- **Loudness explains part of it, not all.** The spectral stage at equal + wet is 4.30 to 14.75 dB louder than the plate's mix. At matched loudness the + spectral rows still cost +2.12 to +18.34 at S1 except shape_max 1 at + w 0.15 (−0.51; the plate +1.44) and flat at w 0.15 (+1.85; plate +1.44), + and +2.17 to +10.53 at S3 against the plate's +0.32 / +0.37. Matched is + the cost had the stage been turned down to the dry bus's loudness; the + measurement ran it as configured. +- **Shaping helps only at shape_max 1** at S1 (+1.17 to +4.30 over flat, 3 + to 5 of 6 rooms); shape_max 2 does not separate from flat there and + shape_max 4 loses 7.91 to 12.79 dB in every room. At S3 every ceiling + reads at or above flat at w 0.30 (+0.98 to +3.71). Where the branch had + "shaped ≥ flat" at shape_max 4 (0.8 s probes), the converged and 40 s + probes say the opposite at S1. +- **Speech envelope** (S1, wet 0.30, rt60 1 s; ProbeConvergence, cabin and + mt5 only, 80 s; the six-room speech grid was not run): the spectral rows + cost +4.30 to +17.09 (cabin) and +9.28 to +14.06 (mt5) against the + canceller alone, where the plate at decay 0.5 raised the limit (cabin + +21.31 against the canceller's +19.75, mt5 +24.24 against +22.68; the + plates' convergence table). + +### Without a canceller (open loop) + +The open loop with the reverb in it against the dry room (no canceller), +160 s probes, six rooms (the grid; rt60 1.1846 s), median of room medians +[range]: + +| | S1, w 0.15 | S1, w 0.30 | S3, w 0.15 | S3, w 0.30 | +|---|---|---|---|---| +| flat | −19.49 [−26.04, −17.15] | −25.45 [−31.80, −23.40] | −19.10 [−24.57, −17.83] | −25.16 [−30.33, −23.98] | +| shape_max 1 | −18.32 [−19.20, −14.71] | −18.52 [−22.13, −15.49] | −16.56 [−19.98, −12.56] | −17.44 [−23.30, −16.46] | +| shape_max 2 | −18.42 [−25.94, −15.68] | −20.96 [−25.84, −15.78] | −18.81 [−21.05, −13.34] | −19.79 [−24.77, −17.54] | +| shape_max 4 | −22.42 [−29.84, −18.22] | −23.69 [−28.87, −17.54] | −17.64 [−19.98, −13.73] | −17.44 [−25.74, −13.83] | + +The plate costs a bare loop at most +5.79 dB by the bound (decay 0.85, w +0.50; [The bare loop](#the-bare-loop-what-the-plate-costs-without-a-canceller)). +HANDOFF item 11's −18.05 to −24.61 for the flat reverb holds with converged +probes: −17.15 to −31.80 over the six rooms. Shaping does not make it safe: +every shaped row is 12.56 to 29.84 dB below the dry room's limit. Gated +(cabin / mt5, S1, w 0.30, rt60 1 s): on converged 80 s probes behind +`MUTAP_SLOW` (`SpectralReverbHost.UnsafeWithoutTheCanceller`) flat +−30.51 / −22.77, shape_max 1 −19.61 / −14.69 (five-set medians); in the +default run on 20 s probes (`UnsafeWithoutTheCancellerShortProbe`, two +sets) −29.80 / −22.42 and −19.79 / −13.98 (means). Behind the canceller, +`CostsMoreThanThePlateShortProbe` (default, cabin, 20 s, two sets) reads the +plate's limit above the spectral reverb's by +14.66 (flat) and +9.78 +(shape_max 1) dB, and `CostsMoreThanThePlateBehindTheCanceller` +(`MUTAP_SLOW`, 80 s) costs the spectral reverb +17.19 and +8.44 dB +against the canceller alone. + +### The hypothesis: does shaping put the reverb where the estimate is worst? + +HANDOFF item 11's hypothesis was that the allowance up to shape_max builds +high-Q resonances where the canceller's estimate is poorest. Measured at +the moment of shaping in every shaped grid run (`spectral_rig.h`'s +`per_bin`; 30 runs a row): per bin of the reverb's 65-bin grid, the +misalignment m_k = |F_k − F̂_k| (F the true band-limited path, both sampled +exactly on the grid), the bin's decay time T_k, its tail energy E_k = +w_k²/(1 − a_k²) and peak gain g_k = (1 − w) + w_k/(1 − a_k); per run, +Spearman's ρ over the bins, the share of Σ E_k in the quarter of the bins +with the largest m_k (a flat reverb's share is 16/65 = 0.246), and a +residual-loop indicator, 20 log10 of max_k m_k g_k over the flat reverb's. +Medians over the 30 runs [min, max]: + +| Row | ρ(T, m) | ρ(T, m / \|F\|) | E share, worst quarter | Peak m·g re flat, dB | +|---|---|---|---|---| +| S1, shape_max 1, w 0.15 | −0.535 [−0.865, −0.200] | +0.161 [−0.235, +0.644] | 0.134 [0.025, 0.217] | −2.456 [−6.914, +0.247] | +| S1, shape_max 2, w 0.15 | −0.589 [−0.876, −0.258] | +0.217 [−0.201, +0.689] | 0.082 [0.007, 0.238] | +2.375 [−4.481, +7.970] | +| S1, shape_max 4, w 0.15 | −0.601 [−0.881, −0.263] | +0.239 [−0.199, +0.680] | 0.036 [0.003, 0.312] | +8.251 [−3.879, +16.717] | +| S1, shape_max 1, w 0.30 | −0.336 [−0.784, −0.033] | +0.495 [−0.204, +0.625] | 0.175 [0.063, 0.230] | −1.015 [−6.522, +1.722] | +| S1, shape_max 2, w 0.30 | −0.437 [−0.843, −0.108] | +0.516 [−0.161, +0.657] | 0.091 [0.013, 0.217] | +3.706 [−4.829, +9.853] | +| S1, shape_max 4, w 0.30 | −0.458 [−0.840, −0.131] | +0.526 [−0.143, +0.640] | 0.050 [0.003, 0.189] | +8.408 [−5.239, +17.698] | +| S3, shape_max 1, w 0.15 | −0.625 [−0.877, −0.243] | +0.168 [−0.240, +0.537] | 0.110 [0.029, 0.230] | −3.250 [−6.873, +1.261] | +| S3, shape_max 2, w 0.15 | −0.683 [−0.942, −0.277] | +0.205 [−0.208, +0.588] | 0.058 [0.007, 0.177] | +0.438 [−5.248, +7.515] | +| S3, shape_max 4, w 0.15 | −0.688 [−0.945, −0.285] | +0.187 [−0.205, +0.585] | 0.022 [0.002, 0.211] | +5.487 [−2.705, +14.620] | +| S3, shape_max 1, w 0.30 | −0.336 [−0.687, −0.109] | +0.404 [−0.125, +0.596] | 0.177 [0.093, 0.230] | −0.848 [−5.504, +1.783] | +| S3, shape_max 2, w 0.30 | −0.381 [−0.755, −0.135] | +0.429 [−0.211, +0.655] | 0.110 [0.030, 0.314] | +4.981 [−2.803, +9.010] | +| S3, shape_max 4, w 0.30 | −0.402 [−0.754, −0.135] | +0.440 [−0.203, +0.654] | 0.059 [0.010, 0.363] | +9.717 [−3.728, +17.130] | + +(rt60 1.1846 s; the rt60 1 s rows at w 0.30 read within 0.11 of these in +ρ and within 0.017 in share, and the sweep prints them.) ρ(E, ·) and ρ(g, ·) equal ρ(T, ·): +all three rise with the allowance. + +- **In absolute terms, no: shaping puts the long decays where the + misalignment is smallest.** ρ(T, m) is negative in every row (median + −0.336 to −0.688), and the shaped reverb puts less of its tail energy in + the worst quarter of the bins than a flat one does (0.022 to 0.178 over + all 18 rows, against 0.246). +- **In relative terms, yes, weakly:** ρ(T, m / |F|) is positive in every + row (median +0.161 to +0.526). The bins where |F̂| is weak are where the + estimate is worst relative to the path. +- **What tracks the loss is the ceiling.** The residual-loop indicator + rises with shape_max: shape_max 1 lowers the worst bin's m·g against flat + (−0.848 to −3.250 dB), shape_max 4 raises it by +5.487 to +9.717 dB, + which matches shape_max 1 being the only ceiling that beats flat at S1 and + shape_max 4 losing there. Across runs, the indicator and the energy share + barely predict a run's own shaped − flat limit: Spearman's ρ over the 30 + runs of a row is −0.239 to +0.699 (indicator) and −0.351 to +0.434 + (share), with no consistent sign. In the measured form: a large ceiling + raises the largest residual m_k·g_k on the grid (+5.487 to +9.717 dB at + shape_max 4) although its long decays sit in bins of below-median + misalignment, so the hypothesis as worded (energy concentrated where the + estimate is worst) is not what the runs show; and none of the per-bin + statistics predicts a single run's outcome well. + +### Spectral audible limits + +The plates' audible protocol ([Audible limits](#audible-limits-the-product-chain-candidates)), +`reverb_audible.py --plate spectral`: the cabin, held note, the ramp from +20 dB under the dry loop's `exact_msg_db` after a 30 s warm-up, the live +canceller from a cold start, the shape refreshed from its F̂ every 256 +blocks (flat rows: never). PROTOCOL.md §7.3's reverb-only flags: the +chain output analysed (after the reverb, before the gain), c as the +cross-check, `--rt60 0.0328` and `--chain-rt60 1.1846` (the spectral +reverb's own T30 at rt60 1.1846 s, flat and shaped rows alike). rt60 +1.1846 s, five seeds, medians, dB over `exact_msg_db`; the first four rows +are the plates' table's (same protocol, this host): + +| Chain | S1 audible (c) | S1 runaway (ramp) | S3 audible (c) | S3 runaway (ramp) | +|---|---|---|---|---| +| dry (no canceller) | +0.18 | +1.08 | +0.21 | +1.62 | +| canceller | +16.57 | +20.70 | +12.37 | +21.37 | +| + 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 | +| + spectral flat, w 0.15 | −7.42 (−6.74) | +4.73 | −17.78 (−2.05) | +11.66 | +| + spectral flat, w 0.30 | −6.30 (−6.42) | −2.82 | −12.77 (−9.24) | +3.56 | +| + shape_max 1, w 0.15 | +11.90 (+11.90) | +18.31 | −13.14 (+4.98) | +17.98 | +| + shape_max 1, w 0.30 | −4.32 (−4.33) | +7.83 | −17.06 (−1.14) | +9.41 | +| + shape_max 2, w 0.15 | +4.52 (−1.39) | +13.52 | −18.66 (+6.70) | +17.07 | +| + shape_max 2, w 0.30 | −16.47 (−10.79) | +2.27 | −19.38 (−6.32) | +8.40 | +| + shape_max 4, w 0.15 | −9.40 (−7.99) | +1.16 | −17.23 (+8.34) | +14.11 | +| + shape_max 4, w 0.30 | −18.18 (−17.49) § | −6.94 | −18.83 (−6.15) | +7.84 | + +§ Four seeds: seed 62 ran away at the ramp's start gain (runaway −20.00, +the floor) and the criterion found no audible event. Over the 80 runs the +criterion raised 23 warnings on the chain output (flat 4 and 5 at S1, w +0.15 / 0.30, and 2 and 5 at S3; shape_max 2 at S1, w 0.30, 2; shape_max 4 +at S1, 2 and 3) and 2 on c; the driver's JSON keeps them. + +- **Every spectral row is audible earlier than the plate at the same wet.** + The best, shape_max 1 at S1, w 0.15, is audible at +11.90 against the + plate's +16.53. 14 of the 16 rows are audible below the dry loop's own + audible limit (+0.18 / +0.21) on the chain output, and 12 of 16 on c. +- **The ramp's runaway falls too**, unlike the plate's: −6.94 to +18.31 + against the plate's +21.46 to +23.31 and the canceller's +20.70 / + +21.37. On the ramp shape_max 1 and 2 run away later than flat at S1 + (+18.31 / +7.83 and +13.52 / +2.27 against +4.73 / −2.82; shape_max 4 + earlier, +1.16 / −6.94), and every shaped row later than flat at S3. +- **At S3 the chain output is flagged near the ramp's start** (−12.77 to + −19.38; the ramp starts at −20 after the warm-up) where c reads −9.24 to + +8.34. What the criterion flags in the chain output there was not + examined. On c as on the chain output, every spectral row sits below the + plate's row at the same wet (at S3, w 0.30, the plate's c cross-check + reads −0.18, its chain output −4.77 ‡). +- 80 runs on 3 jobs, 1380.5 s wall. Every WAV was deleted after its + analysis. + ## What did not converge, separate or get measured - **The speech envelope did not converge** in three of eight convergence @@ -532,6 +920,17 @@ decay and a wet: - **The reverb-only rows at S3, wet 0.30, are audible below the dry loop's limit in some seeds** (the audible table, ‡): real, feedback-dependent, bounded, at the plate's own modes; mechanism not identified. +- **The spectral reverb's S3 audible rows are flagged near the ramp's + start on the chain output** (−12.77 to −19.38, where c reads −9.24 to + +8.34; [Spectral audible limits](#spectral-audible-limits)): the same + family as the plate's S3 lines (‡ in the audible table), not examined + here. +- **The spectral reverb's in-loop rows did not converge at 10 s** and most + had not by 160 s ([Spectral reverb](#spectral-reverb)); its grid ran on + 40 s probes and is reported, not gated. Not run for it: the six-room + speech grid, S3 on speech, converged grid probes, float32 in the loop, + two microphones (the bus-sum shaping has plumbing tests only), and + audible rows beyond the cabin held note at rt60 1.1846 s. - **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 @@ -552,7 +951,8 @@ decay and a wet: 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; +work's; three for the spectral reverb's sweeps and audible runs). 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 @@ -600,6 +1000,45 @@ MUTAP_SLOW=1 build/tests/mutap_tests --gtest_filter='ReverbStageSweep.BareLoopCo | `Shift`, S3 | 4813.9 s (cpu-sum 19196.2 s) | 4 | 420 | | `BareLoopCost` | 53.5 s | 1 | (analytic) | +The spectral reverb ([Spectral reverb](#spectral-reverb)): its plumbing on +every target, the gated host rows (the default run's, then the two +converged rows behind `MUTAP_SLOW`), then the sweep (`SPECTRAL_SWEEP_*` +variables select subsets; the file comment lists them): + +```sh +build/tests/mutap_tests --gtest_filter='spectral_reverb_test*:SpectralReverbConfig*:SpectralReverbRt*' +MUTAP_SLOW_THREADS=3 build/tests/mutap_tests --gtest_filter='SpectralReverbHost.*' +MUTAP_SLOW=1 MUTAP_SLOW_THREADS=3 build/tests/mutap_tests --gtest_filter='SpectralReverbHost.*' +MUTAP_SLOW=1 MUTAP_SLOW_THREADS=3 SPECTRAL_SWEEP_MATS=held build/tests/mutap_tests --gtest_filter='SpectralReverbSweep.ProbeConvergence' +MUTAP_SLOW=1 MUTAP_SLOW_THREADS=3 SPECTRAL_SWEEP_MATS=speech SPECTRAL_SWEEP_PROBES=10,20,40,80 build/tests/mutap_tests --gtest_filter='SpectralReverbSweep.ProbeConvergence' +MUTAP_SLOW=1 MUTAP_SLOW_THREADS=3 SPECTRAL_SWEEP_MATS=held build/tests/mutap_tests --gtest_filter='SpectralReverbSweep.Grid' +``` + +| Test | Wall time | Threads | Jobs | +|---|---|---|---| +| `SpectralReverbHost.DecayAndLevelBesideThePlate` (default) | 0.25 s | 1 | (analytic) | +| `SpectralReverbHost.UnsafeWithoutTheCancellerShortProbe` (default) | 27.03 s | 3 | 12 open-loop bisections (4 after a convergence) | +| `SpectralReverbHost.CostsMoreThanThePlateShortProbe` (default) | 77.25 s | 3 | 6 chain bisections | +| `SpectralReverbHost.UnsafeWithoutTheCanceller` (`MUTAP_SLOW`) | 189.41 s | 3 | 30 open-loop bisections (10 after a convergence) | +| `SpectralReverbHost.CostsMoreThanThePlateBehindTheCanceller` (`MUTAP_SLOW`) | 696.95 s | 3 | 15 chain + 5 open-loop bisections | +| `ProbeConvergence`, held | 11349.4 s (cpu-sum 34003.1 s) | 3 | 290 | +| `ProbeConvergence`, speech | 6849.3 s (cpu-sum 20503.6 s) | 3 | 240 | +| `Grid`, held | 35418.6 s (cpu-sum 106197.4 s) | 3 | 1500 | + +A spectral job bisects the chain or the open loop with the reverb, after +converging the canceller where the shape needs it. The audible rows: + +```sh +python3 tools/notebook/reverb_audible.py --dump build/tools/notebook/karaoke_ramp_dump --work runs --jobs 3 \ + --plate spectral --decays 1.1846 --wets 0.15,0.3 --shape-maxes 0,1,2,4 \ + --configs spec0_1.1846_0.15,spec0_1.1846_0.3,spec1_1.1846_0.15,spec1_1.1846_0.3,spec2_1.1846_0.15,spec2_1.1846_0.3,spec4_1.1846_0.15,spec4_1.1846_0.3 \ + --json audible_spectral.json +``` + +80 runs on 3 jobs: 1380.5 s wall (per-run sum 4071.5 s). Each run's WAVs +are deleted after its analysis; the scratch directory held 139 MB of JSON +and gain logs at the end. + 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): diff --git a/include/mutap/fft.h b/include/mutap/fft.h index 04a4104..1ccd834 100644 --- a/include/mutap/fft.h +++ b/include/mutap/fft.h @@ -39,14 +39,15 @@ namespace tap::mu { // float engine the build selected, and so does the layout of every MuTap // class that holds a basic_real_fft by value: partitioned_fdaf // (fdaf.h), partitioned_fdkf (fd_kalman.h), pem_afc (pem_afc.h), - // residual_suppressor (postfilter.h) and nn_suppressor (nn_suppressor.h). - // Their member functions are weak symbols in every image that - // instantiates them, and a dynamic loader may coalesce weak definitions - // across images (macOS dyld binds them to one definition; on ELF a - // plugin resolves against the executable's exports first), so two + // residual_suppressor (postfilter.h), nn_suppressor (nn_suppressor.h) and + // spectral_reverb (spectral_reverb.h). Their member functions are weak + // symbols in every image that instantiates them, and a dynamic loader may + // coalesce weak definitions across images (macOS dyld binds them to one + // definition; on ELF a plugin resolves against the executable's exports + // first), so two // images built with different engines, loaded into one process, could // run one image's code over the other's layout. - // Each of the five is therefore defined inside + // Each of the six is therefore defined inside // // namespace tap::mu::inline TAP_DSP_FFT_ABI { ... } // @@ -63,7 +64,7 @@ namespace tap::mu { // // THE RULE IS TRANSITIVE. The hazard belongs to any class whose object // layout depends on the float engine and whose mangled name does not - // carry it: one that holds a basic_real_fft, any of the five + // carry it: one that holds a basic_real_fft, any of the six // above at , or an aec_chain by value — directly or // through std::optional / std::array / std::variant / a member struct — // without naming it as a template argument. Such a class must itself be @@ -79,7 +80,7 @@ namespace tap::mu { // harmless: they live in one test image with one engine. // // tests/test_fft_engine_contract.cpp pins the tag on every leg for the - // five (by qualified name at compile time, by typeid at run time); it + // six (by qualified name at compile time, by typeid at run time); it // cannot see a wrapper someone adds later, so the transitive rule is // enforced by review. Two more rules follow: // - never forward-declare a tagged class in plain tap::mu @@ -154,9 +155,9 @@ namespace tap::mu::inline TAP_DSP_FFT_ABI { /// of two): srdif (double always; float by default) 4 ... 2^30; /// vDSP (float, TAP_DSP_FFT_ACCELERATE) 4 ... 2^20; CMSIS-DSP (float, /// TAP_DSP_FFT_CMSIS, the Cortex-M55 default) 32 ... 4096. What a - /// MuTap block size maps to: N = 2 * block_size for the cancellers, - /// analysis_blocks * block_size for the residual suppressor, and - /// 2 * the trained hop for the learned suppressor. + /// MuTap block size maps to: N = 2 * block_size for the cancellers and + /// the spectral reverb, analysis_blocks * block_size for the residual + /// suppressor, and 2 * the trained hop for the learned suppressor. /// /// The message is a string literal from the call site, naming the /// class and how n is derived, followed by MUTAP_FFT_SIZE_RANGES: diff --git a/include/mutap/mutap.h b/include/mutap/mutap.h index b3b3977..bda07e9 100644 --- a/include/mutap/mutap.h +++ b/include/mutap/mutap.h @@ -25,3 +25,4 @@ #include "mutap/pem_afc.h" #include "mutap/postfilter.h" #include "mutap/reverb_stage.h" +#include "mutap/spectral_reverb.h" diff --git a/include/mutap/spectral_reverb.h b/include/mutap/spectral_reverb.h new file mode 100644 index 0000000..85a0fad --- /dev/null +++ b/include/mutap/spectral_reverb.h @@ -0,0 +1,384 @@ +/// @file spectral_reverb.h +/// @brief A per-bin block-frequency-domain reverb for afc_chain's reverb +/// slot, flat or with its decay shaped by the identified feedback path. +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors +// +// WHAT. One complex one-pole per frequency bin on weighted-overlap-add +// spectra (Hann analysis window of 2 * block_size, hop block_size, +// COLA-exact, no synthesis window): +// +// S_k[m] = a_k S_k[m - 1] + X_k[m] +// Y_k[m] = (1 - w) X_k[m] + w_k S_k[m] +// +// so the dry path is the input delayed by exactly one block (latency()) and +// each bin rings down with its own pole a_k. Flat, every a_k gives the +// configured rt60 and every w_k is the wet level w. SHAPED (from_path), the +// per-bin allowance (median|P| / |P_k|)^shape, clamped to [shape_min, +// shape_max], scales w_k up and the decay TIME up in the bins where the +// shaping input |P| is weak, and down where it is strong; the w_k are then +// renormalised so that sum_k w_k^2 equals the flat shape's (the decay times +// are not renormalised, so a shaped reverb's output energy is not the flat +// one's). P is the feedback path the canceller identified, F_hat +// (copy_impulse_response()), or for several microphones on one bus the +// coherent sum of their F_hat (reshape_from_impulse_responses()). +// +// A MEASURED NEGATIVE RESULT. Kept so the measurements that reject it can +// be re-run, not as a product stage: docs/reverb-afc.md, "Spectral reverb" +// (simulation, one host: macOS x86_64, AppleClang 17, Release). Behind PEM + +// FD-Kalman, six band-limited rooms, five seed sets, the held note, it costs +// +2.83 to +30.57 dB of runaway gain (median over the rooms, flat and shaped +// with shape_max 1, 2, 4, wet 0.15 and 0.30, S1 and S3) where the vendored +// Dattorro plate at the same wet costs -1.37 to +0.88; per seed, the +// plate's chain holds more in every room's median for 22 of the 24 spectral +// rows compared, and in 5 of 6 rooms for the other two. Shaping helps only at +// shape_max 1 and never reaches the plate. Without a canceller the flat +// reverb's limit is 17 to 32 dB below the dry room's (a bank of high-Q +// resonators), and shaping does not make it safe. At equal wet it is also +// 4.3 to 14.8 dB louder than the plate's mix. Use the plate +// (reverb_stage.h's reverb_mix over a mono-in reverb). +// +// LATENCY. latency() is block_size: the whole output (dry included) is the +// input one block late. On a chain the block comes out of the latency +// budget (afc_chain's output delay line), which is how the measurements +// above budgeted it: the loop's total delay is the plates' rows'. +// +// THE ABI TAG. The class holds a basic_real_fft by value, so its +// layout follows the build's float FFT engine: it is defined inside +// tap::mu::inline TAP_DSP_FFT_ABI (mutap/fft.h; HANDOFF working note 6 (a)), +// and tests/test_fft_engine_contract.cpp pins it as the sixth embedder. +// Never forward-declare it in plain tap::mu. +// +// THE SIZE GATE. The transform size 2 * block_size comes from +// configuration, so it passes through fft_detail::checked_fft_size before any +// buffer is sized from it: CMSIS-DSP on the Cortex-M55 supports 32 ... 4096 +// only, and an unchecked size there is a release-mode HardFault (working +// note 6 (b)). block_size 16 ... 2048 on that target. +// +// Real-time contract (as fd_kalman.h): the constructor validates its config +// (std::invalid_argument), allocates every buffer and may throw. +// process_block() and every other post-construction entry point are +// noexcept and allocation-free, reshape*() included (std::nth_element and +// std::pow on preallocated storage). Arithmetic is in Sample throughout: +// float32 is first-class and no double runs in the float hot path. + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "mutap/fft.h" + +namespace tap::mu { + + /// How a spectral_reverb's per-bin decay and wet gain are set. + enum class spectral_shaping { + flat, ///< every bin at the configured rt60 and wet; reshape*() is ignored + from_path ///< reshape*() shapes the bins from a feedback-path magnitude + }; + +} // namespace tap::mu + +// Inside the tag: the class holds a basic_real_fft by value. +namespace tap::mu::inline TAP_DSP_FFT_ABI { + + /// Weighted-overlap-add reverb with an independent one-pole decay per + /// frequency bin (see the file comment). An afc_stage: process_block(in, + /// out, n) with n a multiple of block_size(). + template + class spectral_reverb { + public: + struct config { + size_t block_size = 64; ///< hop; a power of two >= 4 with 2 * block_size an FFT size + Sample sample_rate = Sample(48000); ///< Hz, for the rt60 conversion + Sample rt60 = Sample(1); ///< seconds to -60 dB per bin, before shaping + Sample wet = Sample(0.15); ///< w in [0, 1]; the dry gain is 1 - w + spectral_shaping shaping = spectral_shaping::flat; + /// Shaping exponent (from_path): allowance = (median|P| / |P_k|)^shape. + Sample shape = Sample(1); + /// Floor and ceiling on the per-bin allowance: a bin's wet gain + /// scales by it (before the renormalisation) and its decay time + /// by it, so shape_max bounds how much longer and louder the + /// quietest bins of P ring than the flat reverb. + Sample shape_min = Sample(0.05); + Sample shape_max = Sample(4); + }; + + /// @throws std::invalid_argument on a bad config, the FFT size gate + /// included (before any buffer is sized from it). + explicit spectral_reverb(const config& cfg) + : m_cfg(validated(cfg)) + , m_window(fft_detail::checked_fft_size(2 * cfg.block_size, + "spectral_reverb: 2 * block_size" MUTAP_FFT_SIZE_RANGES)) + , m_bins(cfg.block_size + 1) + , m_fft(m_window) + , m_analysis(m_window) + , m_frame(m_window, Sample(0)) + , m_spectrum(m_window, Sample(0)) + , m_time(m_window, Sample(0)) + , m_tail(m_window, Sample(0)) + , m_state(m_window, Sample(0)) + , m_decay(m_bins, Sample(0)) + , m_wet(m_bins, Sample(0)) + , m_scratch(m_bins, Sample(0)) { + for (size_t i = 0; i < m_window; ++i) { + m_analysis[i] = Sample(0.5) + - Sample(0.5) + * std::cos(Sample(2) * std::numbers::pi_v + * static_cast(i) / static_cast(m_window)); + } + reshape_flat(); + } + + size_t block_size() const noexcept { return m_cfg.block_size; } + /// The transform size, 2 * block_size(). + size_t window_size() const noexcept { return m_window; } + /// Delay of the whole output, dry path included, in samples. + size_t latency() const noexcept { return m_cfg.block_size; } + /// Shaping bins: DC, block_size - 1 complex bins, Nyquist; bin k at + /// k * sample_rate / window_size() Hz. + size_t bins() const noexcept { return m_bins; } + Sample wet() const noexcept { return m_cfg.wet; } + Sample rt60() const noexcept { return m_cfg.rt60; } + spectral_shaping shaping() const noexcept { return m_cfg.shaping; } + Sample shape_max() const noexcept { return m_cfg.shape_max; } + /// Bin k's pole per hop and wet gain as currently shaped. @pre k < bins() + Sample bin_decay(size_t k) const noexcept { return m_decay[k]; } + Sample bin_wet(size_t k) const noexcept { return m_wet[k]; } + + /// Wet level in [0, 1], keeping the shape (the per-bin wet gains + /// scale with it) and the tail. Returns false (and changes nothing) + /// for a value outside [0, 1] or not finite. + bool set_wet(Sample wet) noexcept { + if (!std::isfinite(wet) || wet < Sample(0) || wet > Sample(1)) { + return false; + } + const Sample old = m_cfg.wet; + m_cfg.wet = wet; + if (old > Sample(0)) { + const Sample scale = wet / old; + for (auto& w : m_wet) { + w *= scale; + } + } + else { + reshape_flat(); // no shape survives a zero wet level + } + return true; + } + + /// Silence the tail and the analysis history; keeps the shape. + void reset() noexcept { + std::fill(m_frame.begin(), m_frame.end(), Sample(0)); + std::fill(m_tail.begin(), m_tail.end(), Sample(0)); + std::fill(m_state.begin(), m_state.end(), Sample(0)); + } + + /// Every bin at the configured rt60 and wet level. + void reshape_flat() noexcept { + const Sample a = flat_pole(); + for (size_t k = 0; k < m_bins; ++k) { + m_decay[k] = a; + m_wet[k] = m_cfg.wet; + } + } + + /// Shape from a per-bin magnitude of the path: `magnitude` holds + /// `bins` values, bins() of them. Ignored (false) in flat mode, for a + /// null pointer or a wrong count. + bool reshape(const Sample* magnitude, size_t bins) noexcept { + if (m_cfg.shaping != spectral_shaping::from_path || magnitude == nullptr || bins != m_bins) { + return false; + } + std::copy(magnitude, magnitude + m_bins, m_scratch.begin()); + apply_shape(); + return true; + } + + /// Shape from the canceller's F_hat as copy_impulse_response() hands + /// it over (`taps` samples). The response is folded onto the + /// window_size() grid, which samples its transfer function exactly + /// at the bins() frequencies. Ignored (false) in flat mode or while + /// the response is still all zero. + bool reshape_from_impulse_response(const Sample* ir, size_t taps) noexcept { + const Sample* const irs[1] = {ir}; + return reshape_from_impulse_responses(irs, 1, taps); + } + + /// Shape from the coherent bus sum |sum_m F_hat_m| of several + /// microphones' paths (`paths` responses of `taps` samples each): + /// the path the bus sees when every mic is summed into it. + bool reshape_from_impulse_responses(const Sample* const* irs, size_t paths, size_t taps) noexcept { + if (m_cfg.shaping != spectral_shaping::from_path || irs == nullptr || paths == 0 || taps == 0) { + return false; + } + for (size_t m = 0; m < paths; ++m) { + if (irs[m] == nullptr) { + return false; + } + } + std::fill(m_spectrum.begin(), m_spectrum.end(), Sample(0)); + for (size_t m = 0; m < paths; ++m) { + for (size_t i = 0; i < taps; ++i) { + m_spectrum[i % m_window] += irs[m][i]; + } + } + Sample energy = Sample(0); + for (const Sample x : m_spectrum) { + energy += x * x; + } + if (!(energy > k_tiny)) { + return false; // nothing identified yet: keep the current shape + } + m_fft.forward_inplace(m_spectrum.data()); + m_scratch[0] = std::abs(m_spectrum[0]); + m_scratch[m_bins - 1] = std::abs(m_spectrum[1]); + for (size_t k = 1; k + 1 < m_bins; ++k) { + m_scratch[k] = std::hypot(m_spectrum[2 * k], m_spectrum[2 * k + 1]); + } + apply_shape(); + return true; + } + + /// The afc_stage process: `n` samples in, `n` out, delayed by + /// latency(). @pre n is a multiple of block_size() (afc_chain passes + /// exactly block_size()); a remainder is written as silence. In + /// place (in == out) is allowed. + void process_block(const Sample* in, Sample* out, size_t n) noexcept { + const size_t b = m_cfg.block_size; + size_t i = 0; + for (; i + b <= n; i += b) { + hop(in + i, out + i); + } + for (; i < n; ++i) { + out[i] = Sample(0); + } + } + + private: + static constexpr Sample k_tiny = Sample(1e-12); + + static config validated(const config& cfg) { + if (cfg.block_size < 4 || (cfg.block_size & (cfg.block_size - 1)) != 0) { + throw std::invalid_argument("spectral_reverb: block_size must be a power of two >= 4"); + } + if (!std::isfinite(cfg.sample_rate) || !(cfg.sample_rate > Sample(0))) { + throw std::invalid_argument("spectral_reverb: sample_rate must be finite and > 0"); + } + if (!std::isfinite(cfg.rt60) || !(cfg.rt60 > Sample(0))) { + throw std::invalid_argument("spectral_reverb: rt60 must be finite and > 0"); + } + if (!std::isfinite(cfg.wet) || cfg.wet < Sample(0) || cfg.wet > Sample(1)) { + throw std::invalid_argument("spectral_reverb: wet must be finite and in [0, 1]"); + } + if (cfg.shaping != spectral_shaping::flat && cfg.shaping != spectral_shaping::from_path) { + throw std::invalid_argument("spectral_reverb: unknown shaping mode"); + } + if (!std::isfinite(cfg.shape) || cfg.shape < Sample(0)) { + throw std::invalid_argument("spectral_reverb: shape must be finite and >= 0"); + } + if (!std::isfinite(cfg.shape_max) || !(cfg.shape_min > Sample(0)) || !(cfg.shape_max >= cfg.shape_min)) { + throw std::invalid_argument("spectral_reverb: need 0 < shape_min <= shape_max, finite"); + } + return cfg; + } + + /// Per-hop pole giving the configured rt60 (-60 dB) before shaping. + Sample flat_pole() const noexcept { + const Sample hops = m_cfg.rt60 * m_cfg.sample_rate / static_cast(m_cfg.block_size); + return std::pow(Sample(10), Sample(-3) / std::max(hops, Sample(1))); + } + + /// The allowance per bin from the magnitudes in m_scratch (which it + /// overwrites while taking their median). + void apply_shape() noexcept { + // m_spectrum is free here: keep the magnitudes there while + // m_scratch is reordered for the median. + std::copy(m_scratch.begin(), m_scratch.end(), m_spectrum.begin()); + const size_t mid = m_bins / 2; + std::nth_element(m_scratch.begin(), m_scratch.begin() + static_cast(mid), m_scratch.end()); + const Sample median = std::max(m_scratch[mid], k_tiny); + + const Sample a = flat_pole(); + Sample sum = Sample(0); + for (size_t k = 0; k < m_bins; ++k) { + const Sample ratio = median / std::max(m_spectrum[k], k_tiny); + const Sample allowance = std::clamp(std::pow(ratio, m_cfg.shape), m_cfg.shape_min, m_cfg.shape_max); + m_wet[k] = m_cfg.wet * allowance; + m_decay[k] = std::pow(a, Sample(1) / allowance); + sum += m_wet[k] * m_wet[k]; + } + // Renormalise to the flat shape's sum of squared wet gains. + const Sample target = static_cast(m_bins) * m_cfg.wet * m_cfg.wet; + const Sample norm = (sum > k_tiny) ? std::sqrt(target / sum) : Sample(1); + for (auto& w : m_wet) { + w *= norm; + } + } + + /// One hop: block_size samples in, block_size out. + void hop(const Sample* in, Sample* out) noexcept { + const size_t b = m_cfg.block_size; + for (size_t i = 0; i < m_window - b; ++i) { + m_frame[i] = m_frame[i + b]; + } + for (size_t i = 0; i < b; ++i) { + m_frame[m_window - b + i] = in[i]; + } + for (size_t i = 0; i < m_window; ++i) { + m_spectrum[i] = m_frame[i] * m_analysis[i]; + } + m_fft.forward_inplace(m_spectrum.data()); + + const Sample dry = Sample(1) - m_cfg.wet; + // DC and Nyquist are real, in slots 0 and 1 (the packed layout). + m_state[0] = m_decay[0] * m_state[0] + m_spectrum[0]; + m_state[1] = m_decay[m_bins - 1] * m_state[1] + m_spectrum[1]; + m_spectrum[0] = dry * m_spectrum[0] + m_wet[0] * m_state[0]; + m_spectrum[1] = dry * m_spectrum[1] + m_wet[m_bins - 1] * m_state[1]; + for (size_t k = 1; k + 1 < m_bins; ++k) { + const Sample a = m_decay[k]; + const Sample w = m_wet[k]; + const size_t re = 2 * k; + const size_t im = re + 1; + m_state[re] = a * m_state[re] + m_spectrum[re]; + m_state[im] = a * m_state[im] + m_spectrum[im]; + m_spectrum[re] = dry * m_spectrum[re] + w * m_state[re]; + m_spectrum[im] = dry * m_spectrum[im] + w * m_state[im]; + } + + m_fft.inverse(m_spectrum.data(), m_time.data()); + for (size_t i = 0; i < m_window; ++i) { + m_tail[i] += m_time[i]; + } + for (size_t i = 0; i < b; ++i) { + out[i] = m_tail[i]; + } + for (size_t i = 0; i < m_window - b; ++i) { + m_tail[i] = m_tail[i + b]; + } + for (size_t i = m_window - b; i < m_window; ++i) { + m_tail[i] = Sample(0); + } + } + + config m_cfg; + size_t m_window; ///< 2 * block_size, through the size gate + size_t m_bins; + basic_real_fft m_fft; + std::vector m_analysis; ///< Hann, m_window + std::vector m_frame; ///< sliding time-domain window + std::vector m_spectrum; ///< analysis / synthesis scratch + std::vector m_time; ///< inverse-FFT scratch + std::vector m_tail; ///< overlap-add tail + std::vector m_state; ///< one complex one-pole per bin (packed) + std::vector m_decay; ///< per-bin pole per hop + std::vector m_wet; ///< per-bin wet gain + std::vector m_scratch; ///< magnitudes / median workspace + }; + +} // namespace tap::mu::inline TAP_DSP_FFT_ABI diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 0d76d43..38d43cf 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -60,7 +60,8 @@ add_executable(mutap_tests test_nn_suppressor.cpp test_fft_engine_contract.cpp test_portable_random.cpp - test_reverb_mix.cpp) + test_reverb_mix.cpp + test_spectral_reverb.cpp) target_link_libraries(mutap_tests PRIVATE MuTap::MuTap mutap_warnings) @@ -92,6 +93,9 @@ if(NOT MUTAP_BARE_METAL AND NOT CMAKE_CROSSCOMPILING) # 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) + # The spectral reverb in the same slot (docs/reverb-afc.md, "Spectral + # reverb"): its in-loop rows and sweep; its plumbing is test_spectral_reverb.cpp. + target_sources(mutap_tests PRIVATE test_spectral_reverb_host.cpp test_spectral_reverb_sweep.cpp) target_include_directories(mutap_tests PRIVATE ${PROJECT_SOURCE_DIR}/third_party/faust) endif() @@ -129,7 +133,7 @@ else() 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.*:reverb_mix_test/0.*:ReverbMixConfigValidation.*:ReverbMixRtContract.*") + "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.*:spectral_reverb_test/0.*:SpectralReverbConfigValidation.*:SpectralReverbRtContract.*") 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 0fa8201..6d605e8 100644 --- a/tests/bare_metal_main.cpp +++ b/tests/bare_metal_main.cpp @@ -55,18 +55,20 @@ int main() { "FrequencyShifterRtContract.*:" "howl_detector_test/0.*:HowlDetectorCrossPrecision.*:HowlDetectorConfigValidation.*:HowlDetectorRtContract.*:" "howl_guard_test/0.*:HowlGuardConfigValidation.*:HowlGuardRtContract.*:" - "reverb_mix_test/0.*:ReverbMixConfigValidation.*:ReverbMixRtContract.*"; + "reverb_mix_test/0.*:ReverbMixConfigValidation.*:ReverbMixRtContract.*:" + "spectral_reverb_test/0.*:SpectralReverbConfigValidation.*:SpectralReverbRtContract.*"; ::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 108 tests (the FFT suites left with the FFT; FftEngineContract, the + // is 118 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; the reverb - // stage's 7 with the stage); 30 leaves headroom for + // stage's 7 with the stage; the spectral reverb's 9, and its + // FftEngineContract row, with the spectral reverb); 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 8107a53..821304f 100644 --- a/tests/support/decorrelated_loop.h +++ b/tests/support/decorrelated_loop.h @@ -100,18 +100,22 @@ namespace mutap_test { }; /// 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. + /// (support for test_reverb_stage*.cpp and test_spectral_reverb*.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. A stage with a bulk latency (latency(), in + /// samples; 0 without one) has it taken out of the loop's forward_delay, + /// as a chain's output delay line would absorb it, so the loop's total + /// delay stays forward_delay. 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(). + /// Any object with process_block(in, out, n) (in != out) and reset(), + /// and optionally latency(). template static forward_stage of(Stage* stage) { forward_stage f; @@ -120,17 +124,23 @@ namespace mutap_test { static_cast(o)->process_block(in, out, n); }; f.m_reset = [](void* o) { static_cast(o)->reset(); }; + if constexpr (requires { stage->latency(); }) { + f.m_latency = static_cast(stage->latency()); + } 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); } + /// The stage's bulk latency in samples (0 for a stage without latency()). + size_t latency() const { return m_latency; } private: void* m_object = nullptr; void (*m_process)(void*, const Sample*, Sample*, size_t) = nullptr; void (*m_reset)(void*) = nullptr; + size_t m_latency = 0; }; /// Loop with a decorrelator and an auxiliary speaker feed. Same contract @@ -181,9 +191,11 @@ namespace mutap_test { // The IIR shifter is causal and needs no lookahead. const size_t lookahead = (m_cfg.mode == forward_mode::delay_modulation) ? static_cast(0.5 * m_cfg.depth + 1.0) : 0; - if (m_cfg.forward_delay < m_cfg.block_size + lookahead) { + const size_t stage_latency = m_cfg.stage ? m_cfg.stage.latency() : 0; + if (m_cfg.forward_delay < m_cfg.block_size + lookahead + stage_latency) { throw std::invalid_argument("decorrelated_loop: forward_delay too short for this mode"); } + m_read_delay = m_cfg.forward_delay - stage_latency; set_forward_gain_db(m_cfg.forward_gain_db); if (m_cfg.stage) { m_cfg.stage.reset(); @@ -297,7 +309,7 @@ namespace mutap_test { /// The loop signal as it reaches the loudspeaker at absolute time t /// (called once per sample, in order: the shifter carries state). double forward_sample(long long t) { - const double d = static_cast(m_cfg.forward_delay); + const double d = static_cast(m_read_delay); switch (m_cfg.mode) { case forward_mode::delay_modulation: { // Mean delay stays d; the wobble is +-depth/2 around it. @@ -319,7 +331,8 @@ namespace mutap_test { } config m_cfg; - double m_gain = 1.0; + double m_gain = 1.0; + size_t m_read_delay = 0; ///< forward_delay minus the stage's latency (the stage adds it back) /// Absolute sample clock, started mid-buffer so the first blocks' /// history reads land on the zero-filled past rather than wrapping. long long m_t = static_cast(k_history) / 2; diff --git a/tests/support/spectral_rig.h b/tests/support/spectral_rig.h new file mode 100644 index 0000000..cec28bb --- /dev/null +++ b/tests/support/spectral_rig.h @@ -0,0 +1,296 @@ +/// @file spectral_rig.h +/// @brief tap::mu::spectral_reverb as decorrelated_loop's forward stage, on +/// the reverb suites' protocol, and the per-bin measurements that +/// test the shaping hypothesis. +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors +// +// Shared by test_spectral_reverb_host.cpp (the gated rows), +// test_spectral_reverb_sweep.cpp (the MUTAP_SLOW sweep) and karaoke_ramp_dump +// (the audible rows), so a number in +// docs/reverb-afc.md's "Spectral reverb" section means the same measurement +// wherever it came from. Host-only, double precision (a property of the +// algorithm, not of a target's arithmetic). +// +// THE MEASUREMENT is support/karaoke_asg.h's, as the Dattorro rows ran it +// (test_reverb_stage_sweep.cpp), with one addition the shaping needs: the +// canceller converges with the reverb in the loop (flat: a shaped reverb has +// nothing to shape from yet), then the reverb is shaped ONCE from the +// converged canceller's impulse response F_hat (copy_impulse_response(), +// reshape_from_impulse_response(); a no-op in flat mode) and the shape is +// held through every probe of the bisection (each probe starts the reverb's +// tail from silence; reset() keeps the shape). That is the branch's +// protocol (commit 8abe958, "converge, then shape from what the canceller +// has identified"). The reverb's one block of latency comes out of the +// forward delay (decorrelated_loop's forward_stage), so the loop's total +// delay is S1's 480 or S3's 960 samples, as in the plates' rows. +// +// THE HYPOTHESIS (HANDOFF item 11): shaping up to shape_max puts long, +// loud bins where |F_hat| is weakest, which may be where the estimate is +// worst. At the moment of shaping, per bin k of the reverb's grid (65 bins, +// 375 Hz apart at block 64; folding a response onto the 128-point grid +// samples its transfer function exactly there): +// misalignment m_k = |F_k - F_hat_k| (F the true band-limited path) +// relative m_k / |F_k| +// decay time T_k = rt60 * ln(a_flat) / ln(a_k) (the bin's rt60, s) +// tail energy E_k = w_k^2 / (1 - a_k^2) (the bin's one-pole power gain per hop) +// peak gain g_k = (1 - w) + w_k / (1 - a_k) (the bin's gain at its centre) +// and per run: Spearman's rho of T_k against m_k and against m_k / |F_k| +// (E_k and g_k rank as T_k does: all three rise with the allowance), the +// share of sum_k E_k in the quarter of the bins with the largest m_k (16 of +// 65; a flat reverb's share is 16 / 65 = 0.246), and a residual-loop +// indicator: 20 log10 of max_k m_k g_k over the flat reverb's max_k m_k g +// (same rt60 and wet), positive when the shape makes the worst bin of the +// residual loop hotter. An indicator on the hop grid, not the loop's +// frequency response. +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "closed_loop.h" +#include "decorrelated_loop.h" +#include "karaoke_asg.h" +#include "mutap/fft.h" +#include "mutap/spectral_reverb.h" + +namespace mutap_test::spectral { + + namespace kk = karaoke; + + using reverb = tap::mu::spectral_reverb; + + /// One spectral reverb configuration. shape_max 0 is the flat mode. + struct params { + double rt60 = 1.0; + double wet = 0.30; + double shape_max = 0.0; ///< 0: flat (spectral_shaping::flat); else from_path with this ceiling + }; + + inline reverb::config config_of(const params& p) { + reverb::config c; + c.block_size = kk::k_block; + c.sample_rate = kk::k_fs; + c.rt60 = p.rt60; + c.wet = p.wet; + if (p.shape_max > 0.0) { + c.shaping = tap::mu::spectral_shaping::from_path; + c.shape_max = p.shape_max; + c.shape_min = std::min(c.shape_min, p.shape_max); + } + return c; + } + + /// Spearman's rank correlation (average ranks for ties). + inline double spearman(const std::vector& x, const std::vector& y) { + auto ranks = [](const std::vector& v) { + std::vector idx(v.size()); + std::iota(idx.begin(), idx.end(), size_t{0}); + std::sort(idx.begin(), idx.end(), [&](size_t a, size_t b) { return v[a] < v[b]; }); + std::vector r(v.size()); + for (size_t i = 0; i < idx.size();) { + size_t j = i; + while (j + 1 < idx.size() && v[idx[j + 1]] == v[idx[i]]) { + ++j; + } + const double avg = 0.5 * static_cast(i + j); + for (size_t k = i; k <= j; ++k) { + r[idx[k]] = avg; + } + i = j + 1; + } + return r; + }; + const auto rx = ranks(x); + const auto ry = ranks(y); + const auto n = static_cast(x.size()); + const auto mx = std::accumulate(rx.begin(), rx.end(), 0.0) / n; + const auto my = std::accumulate(ry.begin(), ry.end(), 0.0) / n; + double sxy = 0.0; + double sxx = 0.0; + double syy = 0.0; + for (size_t i = 0; i < rx.size(); ++i) { + sxy += (rx[i] - mx) * (ry[i] - my); + sxx += (rx[i] - mx) * (rx[i] - mx); + syy += (ry[i] - my) * (ry[i] - my); + } + return (sxx > 0.0 && syy > 0.0) ? sxy / std::sqrt(sxx * syy) : std::nan(""); + } + + /// A response folded onto an n-point grid and transformed: the transfer + /// function at k * 2 pi / n, k = 0 .. n / 2, as (re, im) pairs. + inline std::vector> folded_spectrum(const std::vector& h, size_t n) { + std::vector buf(n, 0.0); + for (size_t i = 0; i < h.size(); ++i) { + buf[i % n] += h[i]; + } + tap::mu::real_fft fft(n); + fft.forward_inplace(buf.data()); + std::vector> out(n / 2 + 1); + out[0] = {buf[0], 0.0}; + out[n / 2] = {buf[1], 0.0}; + for (size_t k = 1; k < n / 2; ++k) { + out[k] = {buf[2 * k], buf[2 * k + 1]}; + } + return out; + } + + /// The per-bin picture at the moment of shaping (see the file comment). + struct bin_stats { + std::vector mis; ///< m_k + std::vector rel; ///< m_k / |F_k| + std::vector t60; ///< T_k, s + std::vector e; ///< E_k + double rho_t_mis = 0.0; + double rho_t_rel = 0.0; + double share_top = 0.0; ///< share of sum E_k in the top quarter of m_k + double peak_rel_db = 0.0; ///< max m_k g_k, shaped over flat, dB + double misalign_db = 0.0; ///< 10 log10(sum m_k^2 / sum |F_k|^2) on the grid + }; + + inline bin_stats per_bin(const std::vector& path, const std::vector& f_hat, const reverb& rev) { + const size_t n = rev.window_size(); + const auto F = folded_spectrum(path, n); + const auto Fh = folded_spectrum(f_hat, n); + bin_stats s; + const double a_flat = [&] { + const double hops = rev.rt60() * kk::k_fs / static_cast(rev.block_size()); + return std::pow(10.0, -3.0 / std::max(hops, 1.0)); + }(); + const double w_flat = rev.wet(); + const double g_flat = (1.0 - w_flat) + w_flat / (1.0 - a_flat); + double num = 0.0; + double den = 0.0; + double peak_s = 0.0; + double peak_f = 0.0; + for (size_t k = 0; k < rev.bins(); ++k) { + const double dr = F[k].first - Fh[k].first; + const double di = F[k].second - Fh[k].second; + const double m = std::hypot(dr, di); + const double f = std::hypot(F[k].first, F[k].second); + const double a = rev.bin_decay(k); + const double w = rev.bin_wet(k); + s.mis.push_back(m); + s.rel.push_back(m / std::max(f, 1e-300)); + s.t60.push_back(rev.rt60() * std::log(a_flat) / std::log(a)); + s.e.push_back(w * w / (1.0 - a * a)); + peak_s = std::max(peak_s, m * ((1.0 - w_flat) + w / (1.0 - a))); + peak_f = std::max(peak_f, m * g_flat); + num += m * m; + den += f * f; + } + s.misalign_db = 10.0 * std::log10(num / den); + s.peak_rel_db = 20.0 * std::log10(peak_s / peak_f); + s.rho_t_mis = spearman(s.t60, s.mis); + s.rho_t_rel = spearman(s.t60, s.rel); + std::vector idx(s.mis.size()); + std::iota(idx.begin(), idx.end(), size_t{0}); + std::sort(idx.begin(), idx.end(), [&](size_t a, size_t b) { return s.mis[a] > s.mis[b]; }); + const size_t top = s.mis.size() / 4; // 16 of 65 + double e_top = 0.0; + for (size_t i = 0; i < top; ++i) { + e_top += s.e[idx[i]]; + } + s.share_top = e_top / std::accumulate(s.e.begin(), s.e.end(), 0.0); + return s; + } + + /// The stage's broadband level change for white input, dB: output over + /// input energy of `seconds` of white noise (the first second skipped), + /// through a copy of `rev` (its shape as it stands). + inline double level_db(const reverb& rev, double seconds = 20.0) { + reverb r = rev; + const size_t b = r.block_size(); + r.reset(); + const auto x = white_near_end(static_cast(seconds * kk::k_fs) / b * b, 4242); + std::vector y(x.size()); + for (size_t i = 0; i < x.size(); i += b) { + r.process_block(&x[i], &y[i], b); + } + const auto skip = static_cast(kk::k_fs); + double ex = 0.0; + double ey = 0.0; + for (size_t i = skip; i < x.size(); ++i) { + ex += x[i - r.latency()] * x[i - r.latency()]; + ey += y[i] * y[i]; + } + return 10.0 * std::log10(ey / ex); + } + + /// The reverb's all-wet impulse response (wet 1, flat): its return r, + /// `seconds` long, the impulse at the first sample of a hop. + inline std::vector wet_impulse_response(double rt60, double seconds) { + params p; + p.rt60 = rt60; + p.wet = 1.0; + reverb r(config_of(p)); + const size_t b = r.block_size(); + const size_t n = static_cast(seconds * kk::k_fs) / b * b; + std::vector x(n, 0.0); + std::vector y(n, 0.0); + x[0] = 1.0; + for (size_t i = 0; i < n; i += b) { + r.process_block(&x[i], &y[i], b); + } + return {y.begin() + static_cast(r.latency()), y.end()}; + } + + /// One run: converge, shape, bisect. + struct run { + double exact_db = 0.0; + double chain_db = 0.0; ///< with the canceller + double open_rev_db = std::numeric_limits::quiet_NaN(); ///< no canceller, the same reverb + double level = 0.0; ///< level_db() of the reverb as shaped + bool shaped = false; ///< reshape took effect + bin_stats bins; + }; + + /// Bisected limits on `path` (see the file comment): with `chain`, of the + /// canceller + spectral reverb `prm`; with `open_reverb`, of the open loop + /// (no canceller) with the same reverb, as shaped, over [exact + open_lo, + /// exact + open_hi]. The canceller converges (and the reverb is shaped) + /// either way. + inline run measure(const std::vector& path, kk::material mat, const kk::protocol& p, unsigned seed, + const params& prm, bool chain, bool open_reverb, double open_lo = -40.0, double open_hi = 10.0) { + run r; + r.exact_db = exact_msg_db(path, p.delay); + const auto v_conv = kk::near_end(mat, p.converge_blocks, seed); + const auto v_probe = kk::near_end(mat, p.probe_blocks, seed + 10); + reverb rev(config_of(prm)); + kk::setup s; + s.mat = mat; + s.stage = forward_stage::of(&rev); + const auto cfg = kk::loop_config(path, s, p.delay, nullptr); + + kk::kalman_afc::config c; + c.fdaf.block_size = kk::k_block; + c.fdaf.partitions = kk::k_parts; + kk::kalman_afc afc(c); + std::vector ir(afc.filter_length(), 0.0); + // The open loop alone of a flat reverb needs no canceller (nothing to + // shape from): F_hat stays zero and the per-bin misalignment is |F|. + if (chain || prm.shape_max > 0.0) { + kk::converge(afc, cfg, r.exact_db, p, v_conv); + afc.copy_impulse_response(ir.data()); + } + r.shaped = rev.reshape_from_impulse_response(ir.data(), ir.size()); + r.bins = per_bin(path, ir, rev); + r.level = level_db(rev); + + r.chain_db = std::numeric_limits::quiet_NaN(); + if (chain) { + r.chain_db = + decorrelated_msg_db(cfg, &afc, v_probe, r.exact_db + p.chain_lo, r.exact_db + p.chain_hi, p.tol_db); + } + if (open_reverb) { + r.open_rev_db = decorrelated_msg_db(cfg, static_cast(nullptr), v_probe, + r.exact_db + open_lo, r.exact_db + open_hi, p.tol_db); + } + return r; + } + +} // namespace mutap_test::spectral diff --git a/tests/test_fft_engine_contract.cpp b/tests/test_fft_engine_contract.cpp index 0a13c9d..51fd164 100644 --- a/tests/test_fft_engine_contract.cpp +++ b/tests/test_fft_engine_contract.cpp @@ -6,9 +6,10 @@ // // THE ABI TAG. Every MuTap class that holds a basic_real_fft by value // (partitioned_fdaf, partitioned_fdkf, pem_afc, residual_suppressor, -// nn_suppressor) is defined inside tap::mu::inline TAP_DSP_FFT_ABI, so the -// build's float FFT engine is part of its mangled name (mutap/fft.h has the -// why). Pinned two ways, on every leg that runs the battery: +// nn_suppressor and, the sixth, spectral_reverb) is defined inside +// tap::mu::inline TAP_DSP_FFT_ABI, so the build's float FFT engine is part of +// its mangled name (mutap/fft.h has the why). Pinned two ways, on every leg +// that runs the battery: // - at compile time, by naming each class through the tag namespace // explicitly (tap::mu::TAP_DSP_FFT_ABI::partitioned_fdaf): qualified // lookup into a namespace finds only what is declared there, so a class @@ -41,6 +42,7 @@ // checked before the core is built) // residual_suppressor N = analysis_blocks * block_size // nn_suppressor N = 2 * the weights' hop +// spectral_reverb N = 2 * block_size #include #include @@ -63,6 +65,7 @@ #include "mutap/nn_suppressor.h" #include "mutap/pem_afc.h" #include "mutap/postfilter.h" +#include "mutap/spectral_reverb.h" namespace { @@ -76,7 +79,8 @@ namespace { && std::is_same_v, tagged::partitioned_fdkf> && std::is_same_v, tagged::pem_afc> && std::is_same_v, tagged::residual_suppressor> - && std::is_same_v, tagged::nn_suppressor>; + && std::is_same_v, tagged::nn_suppressor> + && std::is_same_v, tagged::spectral_reverb>; static_assert(k_declared_in_tag, "every float FFT embedder is defined inside the ABI tag"); static_assert(k_declared_in_tag, "every double FFT embedder is defined inside the ABI tag"); @@ -100,6 +104,7 @@ namespace { tap::mu::pem_afc, tap::mu::partitioned_fdkf>>(); expect_only_this_builds_tag>(); expect_only_this_builds_tag>(); + expect_only_this_builds_tag>(); expect_only_this_builds_tag>(); expect_only_this_builds_tag>(); expect_only_this_builds_tag>(); @@ -230,6 +235,23 @@ namespace { EXPECT_FALSE(suppressor_accepts(k_top_power, 4)); } + template + bool spectral_reverb_accepts(std::size_t block) { + typename tap::mu::spectral_reverb::config cfg; + cfg.block_size = block; + return accepted([&] { return tap::mu::spectral_reverb(cfg); }); + } + + template + void expect_spectral_reverb_follows_the_range() { + for (const std::size_t b : block_sweep()) { + EXPECT_EQ(spectral_reverb_accepts(b), fft::supports_size(2 * b)) + << "spectral_reverb block " << b; + } + EXPECT_FALSE(spectral_reverb_accepts(std::size_t{1} << 30)); + EXPECT_FALSE(spectral_reverb_accepts(k_top_power)); + } + template void expect_nn_follows_the_range() { for (std::size_t hop = 16; hop <= 4096; hop *= 2) { @@ -267,6 +289,11 @@ TEST(FftEngineContract, LearnedSuppressorHopIsGatedBySupportsSize) { expect_nn_follows_the_range(); } +TEST(FftEngineContract, SpectralReverbBlockSizeIsGatedBySupportsSize) { + expect_spectral_reverb_follows_the_range(); + expect_spectral_reverb_follows_the_range(); +} + // The plan's rows, stated per build rather than read from the predicate: the // double profile always runs the srdif engine (4 ... 2^30); the float // profile runs CMSIS-DSP (32 ... 4096) on the Cortex-M55 leg, vDSP @@ -277,18 +304,26 @@ TEST(FftEngineContract, ConfiguredSizesThisBuildRejects) { EXPECT_TRUE(fdaf_accepts(4096)); // N = 8192 #if defined(TAP_DSP_FFT_CMSIS) static_assert(fft::k_min_size == 32 && fft::k_max_size == 4096); - EXPECT_FALSE(fdaf_accepts(8)); // N = 16 - EXPECT_FALSE(fdaf_accepts(4096)); // N = 8192 - EXPECT_TRUE(fdaf_accepts(16)); // N = 32 - EXPECT_TRUE(fdaf_accepts(2048)); // N = 4096 - EXPECT_FALSE(suppressor_accepts(4, 4)); // N = 16 - EXPECT_FALSE(nn_accepts(4096)); // N = 8192 + EXPECT_FALSE(fdaf_accepts(8)); // N = 16 + EXPECT_FALSE(fdaf_accepts(4096)); // N = 8192 + EXPECT_TRUE(fdaf_accepts(16)); // N = 32 + EXPECT_TRUE(fdaf_accepts(2048)); // N = 4096 + EXPECT_FALSE(suppressor_accepts(4, 4)); // N = 16 + EXPECT_FALSE(nn_accepts(4096)); // N = 8192 + EXPECT_FALSE(spectral_reverb_accepts(8)); // N = 16 + EXPECT_FALSE(spectral_reverb_accepts(4096)); // N = 8192 + EXPECT_TRUE(spectral_reverb_accepts(16)); // N = 32 + EXPECT_TRUE(spectral_reverb_accepts(2048)); // N = 4096 #else EXPECT_TRUE(fdaf_accepts(8)); // N = 16 EXPECT_TRUE(fdaf_accepts(4096)); // N = 8192 EXPECT_TRUE(suppressor_accepts(4, 4)); EXPECT_TRUE(nn_accepts(4096)); + EXPECT_TRUE(spectral_reverb_accepts(8)); // N = 16 + EXPECT_TRUE(spectral_reverb_accepts(4096)); // N = 8192 #endif + EXPECT_TRUE(spectral_reverb_accepts(8)); // N = 16 + EXPECT_TRUE(spectral_reverb_accepts(4096)); // N = 8192 } // pem_afc states its own gate before its core is constructed, so the diff --git a/tests/test_spectral_reverb.cpp b/tests/test_spectral_reverb.cpp new file mode 100644 index 0000000..c3f4a8e --- /dev/null +++ b/tests/test_spectral_reverb.cpp @@ -0,0 +1,336 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors +// +// tap::mu::spectral_reverb's plumbing (mutap/spectral_reverb.h): the dry path +// one block late, the flat decay against its rt60, the flat mode, the shape's +// bounds and renormalisation, the multi-mic bus sum, the afc_stage shape, +// config validation and the real-time contract. Float and double; this file +// builds on every target and both emulated selections run its float suite +// (tests/bare_metal_main.cpp, TEST_FILTER in tests/CMakeLists.txt). What the +// reverb does inside the loop is test_spectral_reverb_host.cpp (host-only) +// and docs/reverb-afc.md. + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "mutap/afc_chain.h" +#include "mutap/spectral_reverb.h" + +namespace { + + using tap::mu::spectral_shaping; + + template + class spectral_reverb_test : public ::testing::Test {}; + using sample_types = ::testing::Types; + TYPED_TEST_SUITE(spectral_reverb_test, sample_types); + + 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; + } + + template + typename tap::mu::spectral_reverb::config base_config() { + typename tap::mu::spectral_reverb::config c; + c.block_size = 64; + c.sample_rate = Sample(48000); + c.rt60 = Sample(1); + c.wet = Sample(0.3); + return c; + } + + /// A decaying, sign-changing test response of `taps` samples (no libm in + /// the values beyond one pow per tap). + template + std::vector test_response(size_t taps, unsigned salt) { + std::vector h(taps); + for (size_t i = 0; i < taps; ++i) { + const double s = ((i * 7 + salt * 13) % 11 < 5) ? -1.0 : 1.0; + h[i] = static_cast(s * std::pow(0.995, static_cast(i)) + * (1.0 + 0.1 * static_cast((i + salt) % 5))); + } + return h; + } + + // Wet 0 is the input exactly one block late (latency()), from the first + // sample, across calls of one and of three blocks, and in place. + TYPED_TEST(spectral_reverb_test, DryPathIsTheInputOneBlockLate) { + using sample = TypeParam; + auto c = base_config(); + c.wet = sample(0); + tap::mu::spectral_reverb rev(c); + const size_t b = rev.block_size(); + ASSERT_EQ(rev.latency(), b); + ASSERT_EQ(rev.window_size(), 2 * b); + ASSERT_EQ(rev.bins(), b + 1); + const size_t n = 40 * b; + const auto x = ramp_signal(n); + auto y = x; // in place + for (size_t i = 0; i < n;) { + const size_t len = ((i / b) % 2 == 0) ? b : 3 * b; + const size_t m = std::min(len, n - i); + rev.process_block(y.data() + i, y.data() + i, m); + i += m; + } + // Measured worst error (macOS x86_64, AppleClang 17, Release): float + // 2.98e-07, double 6.66e-16. + const sample tol = std::is_same_v ? sample(5e-6) : sample(1e-13); + size_t bad = 0; + sample worst = sample(0); + for (size_t i = 0; i < n; ++i) { + const sample expect = (i < b) ? sample(0) : x[i - b]; + worst = std::max(worst, std::abs(y[i] - expect)); + if (!(std::abs(y[i] - expect) <= tol)) { + ++bad; + } + } + std::printf("[ measured ] dry path, worst |error| %.3g\n", static_cast(worst)); + EXPECT_EQ(bad, 0U) << "worst " << worst; + } + + // The flat decay means its rt60: the all-wet response's per-block energy + // falls 60 dB in rt60 seconds (least-squares slope from -5 to -35 dB of + // the energy decay curve, Schroeder's T30), the impulse at the start of a + // hop. Measured (macOS x86_64, AppleClang 17, Release): rt60 0.5 s reads + // T30 0.5000 s and rt60 1.0 s reads 1.0000 s, float and double; the 2 % + // allowance is for other targets' FFT engines. + TYPED_TEST(spectral_reverb_test, FlatDecayMatchesItsRt60) { + using sample = TypeParam; + for (const sample rt60 : {sample(0.5), sample(1)}) { + auto c = base_config(); + c.wet = sample(1); + c.rt60 = rt60; + tap::mu::spectral_reverb rev(c); + const size_t b = rev.block_size(); + const size_t n = static_cast(1.5 * static_cast(rt60) * 48000.0) / b * b; + std::vector x(n, sample(0)); + std::vector y(n, sample(0)); + x[0] = sample(1); + for (size_t i = 0; i < n; i += b) { + rev.process_block(&x[i], &y[i], b); + } + // Schroeder integration in double (a test-side measurement). + std::vector edc(n); + double acc = 0.0; + for (size_t i = n; i-- > 0;) { + acc += static_cast(y[i]) * static_cast(y[i]); + edc[i] = acc; + } + double sx = 0.0; + double sy = 0.0; + double sxx = 0.0; + double sxy = 0.0; + double m = 0.0; + for (size_t i = 0; i < n; ++i) { + const double db = 10.0 * std::log10(edc[i] / edc[0]); + if (db <= -5.0 && db >= -35.0) { + const double t = static_cast(i) / 48000.0; + sx += t; + sy += db; + sxx += t * t; + sxy += t * db; + m += 1.0; + } + } + ASSERT_GT(m, 100.0); + const double t30 = -60.0 / ((m * sxy - sx * sy) / (m * sxx - sx * sx)); + std::printf("[ measured ] rt60 %.2f s: T30 %.4f s\n", static_cast(rt60), t30); + EXPECT_NEAR(t30, static_cast(rt60), 0.02 * static_cast(rt60)) << "rt60 " << rt60; + } + } + + // Flat mode ignores every reshape; from_path mode with an all-zero + // response keeps its shape. + TYPED_TEST(spectral_reverb_test, FlatModeAndEmptyResponsesKeepTheShape) { + using sample = TypeParam; + auto c = base_config(); + const auto h = test_response(1024, 1); + { + tap::mu::spectral_reverb flat(c); + EXPECT_EQ(flat.shaping(), spectral_shaping::flat); + EXPECT_FALSE(flat.reshape_from_impulse_response(h.data(), h.size())); + std::vector mag(flat.bins(), sample(1)); + mag[3] = sample(1e-3); + EXPECT_FALSE(flat.reshape(mag.data(), mag.size())); + for (size_t k = 0; k < flat.bins(); ++k) { + EXPECT_EQ(flat.bin_decay(k), flat.bin_decay(0)); + EXPECT_EQ(flat.bin_wet(k), c.wet); + } + } + c.shaping = spectral_shaping::from_path; + tap::mu::spectral_reverb rev(c); + const sample a0 = rev.bin_decay(0); + const std::vector zeros(1024, sample(0)); + EXPECT_FALSE(rev.reshape_from_impulse_response(zeros.data(), zeros.size())); + EXPECT_FALSE(rev.reshape_from_impulse_response(nullptr, 1024)); + std::vector short_mag(rev.bins() - 1, sample(1)); + EXPECT_FALSE(rev.reshape(short_mag.data(), short_mag.size())); + for (size_t k = 0; k < rev.bins(); ++k) { + EXPECT_EQ(rev.bin_decay(k), a0); + } + } + + // The shape: the allowance (median|P| / |P_k|), clamped to [shape_min, + // shape_max], lengthens the decay time of the weak bins by at most + // shape_max and shortens the strong ones by at most shape_min; the wet + // gains keep the flat shape's sum of squares; ordering follows |P|. + TYPED_TEST(spectral_reverb_test, ShapeIsBoundedAndRenormalised) { + using sample = TypeParam; + for (const sample cap : {sample(1), sample(2), sample(4)}) { + auto c = base_config(); + c.shaping = spectral_shaping::from_path; + c.shape_max = cap; + tap::mu::spectral_reverb rev(c); + const sample a_flat = rev.bin_decay(0); + std::vector mag(rev.bins()); + for (size_t k = 0; k < mag.size(); ++k) { + mag[k] = static_cast(std::pow(10.0, -2.0 + 4.0 * static_cast(k) / 64.0)); // 1e-2 .. 1e2 + } + ASSERT_TRUE(rev.reshape(mag.data(), mag.size())); + double sum_w2 = 0.0; + for (size_t k = 0; k < rev.bins(); ++k) { + // Decay time relative to flat: ln(a_flat) / ln(a_k). + const double ratio = + std::log(static_cast(a_flat)) / std::log(static_cast(rev.bin_decay(k))); + EXPECT_LE(ratio, static_cast(cap) * (1.0 + 1e-4)) << "bin " << k; + EXPECT_GE(ratio, static_cast(c.shape_min) * (1.0 - 1e-4)) << "bin " << k; + if (k > 0) { + EXPECT_LE(rev.bin_decay(k), rev.bin_decay(k - 1)) << "bin " << k; // stronger |P|, shorter decay + } + sum_w2 += static_cast(rev.bin_wet(k)) * static_cast(rev.bin_wet(k)); + } + const double target = static_cast(rev.bins()) * static_cast(c.wet * c.wet); + EXPECT_NEAR(sum_w2 / target, 1.0, 1e-4) << "cap " << cap; + // The weakest bin rings exactly cap times longer than flat. + EXPECT_NEAR(std::log(static_cast(a_flat)) / std::log(static_cast(rev.bin_decay(0))), + static_cast(cap), 1e-3 * static_cast(cap)); + } + } + + // Several microphones on one bus: shaping from their responses is + // shaping from the coherent sum of them. + TYPED_TEST(spectral_reverb_test, BusShapeIsTheCoherentSum) { + using sample = TypeParam; + auto c = base_config(); + c.shaping = spectral_shaping::from_path; + const auto h1 = test_response(1024, 1); + const auto h2 = test_response(1024, 4); + std::vector sum(1024); + for (size_t i = 0; i < sum.size(); ++i) { + sum[i] = h1[i] + h2[i]; + } + tap::mu::spectral_reverb a(c); + tap::mu::spectral_reverb b(c); + const sample* const irs[2] = {h1.data(), h2.data()}; + ASSERT_TRUE(a.reshape_from_impulse_responses(irs, 2, 1024)); + ASSERT_TRUE(b.reshape_from_impulse_response(sum.data(), sum.size())); + const sample tol = std::is_same_v ? sample(1e-5) : sample(1e-12); + for (size_t k = 0; k < a.bins(); ++k) { + EXPECT_NEAR(a.bin_decay(k), b.bin_decay(k), tol) << k; + EXPECT_NEAR(a.bin_wet(k), b.bin_wet(k), tol) << k; + } + const sample* const bad[2] = {h1.data(), nullptr}; + EXPECT_FALSE(a.reshape_from_impulse_responses(bad, 2, 1024)); + } + + // A call that is not a whole number of blocks: the whole blocks are + // processed and the remainder is silence. reset() silences the tail and + // keeps the shape; set_wet() scales the shaped gains. + TYPED_TEST(spectral_reverb_test, RemainderResetAndWet) { + using sample = TypeParam; + auto c = base_config(); + c.shaping = spectral_shaping::from_path; + tap::mu::spectral_reverb rev(c); + const auto h = test_response(1024, 2); + ASSERT_TRUE(rev.reshape_from_impulse_response(h.data(), h.size())); + const size_t b = rev.block_size(); + const auto x = ramp_signal(2 * b + 5); + std::vector y(x.size(), sample(7)); + rev.process_block(x.data(), y.data(), x.size()); + for (size_t i = 2 * b; i < y.size(); ++i) { + EXPECT_EQ(y[i], sample(0)) << i; + } + const sample w5 = rev.bin_wet(5); + const sample a5 = rev.bin_decay(5); + rev.reset(); + std::vector zeros(b, sample(0)); + std::vector out(b, sample(1)); + rev.process_block(zeros.data(), out.data(), b); + for (const sample v : out) { + EXPECT_EQ(v, sample(0)); + } + EXPECT_EQ(rev.bin_decay(5), a5); + EXPECT_TRUE(rev.set_wet(sample(0.15))); + EXPECT_NEAR(rev.bin_wet(5), w5 * sample(0.5), std::abs(w5) * sample(1e-5)); + EXPECT_EQ(rev.wet(), sample(0.15)); + EXPECT_FALSE(rev.set_wet(sample(1.5))); + EXPECT_FALSE(rev.set_wet(std::numeric_limits::quiet_NaN())); + EXPECT_EQ(rev.wet(), sample(0.15)); + } + + // The class fits afc_chain's reverb slot. + TYPED_TEST(spectral_reverb_test, FitsTheChainSlot) { + using sample = TypeParam; + static_assert(tap::mu::afc_stage, sample>); + tap::mu::spectral_reverb rev(base_config()); + tap::mu::afc_stage_ref ref(&rev); + EXPECT_TRUE(static_cast(ref)); + } + + TEST(SpectralReverbConfigValidation, RejectsBadConfigs) { + using rv = tap::mu::spectral_reverb; + auto make = [](auto edit) { + rv::config c; + edit(c); + return rv(c); + }; + EXPECT_THROW(make([](rv::config& c) { c.block_size = 100; }), std::invalid_argument); + EXPECT_THROW(make([](rv::config& c) { c.block_size = 2; }), std::invalid_argument); + EXPECT_THROW(make([](rv::config& c) { c.sample_rate = 0.0F; }), std::invalid_argument); + EXPECT_THROW(make([](rv::config& c) { c.rt60 = 0.0F; }), std::invalid_argument); + EXPECT_THROW(make([](rv::config& c) { c.rt60 = std::numeric_limits::infinity(); }), + std::invalid_argument); + EXPECT_THROW(make([](rv::config& c) { c.wet = -0.1F; }), std::invalid_argument); + EXPECT_THROW(make([](rv::config& c) { c.wet = 1.5F; }), std::invalid_argument); + EXPECT_THROW(make([](rv::config& c) { c.wet = std::nanf(""); }), std::invalid_argument); + EXPECT_THROW(make([](rv::config& c) { c.shape = -1.0F; }), std::invalid_argument); + EXPECT_THROW(make([](rv::config& c) { c.shape_min = 0.0F; }), std::invalid_argument); + EXPECT_THROW(make([](rv::config& c) { c.shape_max = 0.01F; }), std::invalid_argument); // < shape_min + EXPECT_THROW(make([](rv::config& c) { c.shaping = static_cast(7); }), std::invalid_argument); + EXPECT_NO_THROW(make([](rv::config& c) { c.wet = 1.0F; })); + EXPECT_NO_THROW(make([](rv::config& c) { + c.shaping = spectral_shaping::from_path; + c.shape_max = 1.0F; + })); + } + + TEST(SpectralReverbRtContract, PostConstructionEntryPointsAreNoexcept) { + using rv = tap::mu::spectral_reverb; + static_assert(noexcept(std::declval().process_block(nullptr, nullptr, 0))); + static_assert(noexcept(std::declval().reshape(nullptr, 0))); + static_assert(noexcept(std::declval().reshape_from_impulse_response(nullptr, 0))); + static_assert(noexcept(std::declval().reshape_from_impulse_responses(nullptr, 0, 0))); + static_assert(noexcept(std::declval().reshape_flat())); + static_assert(noexcept(std::declval().set_wet(0.1F))); + static_assert(noexcept(std::declval().reset())); + static_assert(noexcept(std::declval().latency())); + static_assert(noexcept(std::declval().bin_decay(0))); + SUCCEED(); + } + +} // namespace diff --git a/tests/test_spectral_reverb_host.cpp b/tests/test_spectral_reverb_host.cpp new file mode 100644 index 0000000..f3f3157 --- /dev/null +++ b/tests/test_spectral_reverb_host.cpp @@ -0,0 +1,388 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors +// +// The spectral reverb behind the canceller (docs/reverb-afc.md, "Spectral +// reverb"): what the reverb is (T30 against the plate's, its level) and the +// gated in-loop rows, on support/spectral_rig.h's measurement. Host-only, +// double precision, like every closed-loop claim (tests/CMakeLists.txt); the +// class's plumbing on every target is test_spectral_reverb.cpp, and the +// MUTAP_SLOW sweep behind these rows is test_spectral_reverb_sweep.cpp. +// +// Measured on macOS 15.7 x86_64 (i9-8950HK), AppleClang 17, Release; each +// test's comment has its numbers. + +#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" +#include "support/spectral_rig.h" + +namespace { + + namespace kk = mutap_test::karaoke; + namespace rv = mutap_test::reverb; + namespace sp = mutap_test::spectral; + using mutap_test::median; + using mutap_test::seed_in_set; + + /// The as-shipped plate's decay-0.5 T30, L (ReverbStage.PlateDecayAndMixLevel). + constexpr double k_plate_t30 = 1.1846; + + /// 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(); + } + } + + std::string per_seed(const std::vector& v) { + std::string s; + for (const double x : v) { + char b[16]; + std::snprintf(b, sizeof b, " %+.2f", x); + s += b; + } + return s; + } + + /// The held note on a band-limited room at S1, wet 0.30, rt60 1 s + /// (ProbeConvergence's configuration), five seed sets. Per seed: the dry + /// open loop's limit, and for each spectral row (shape_max 0 = flat) the + /// limit of the open loop with the reverb (no canceller) or of the chain; + /// absolute forward gains, dB; bisection to 0.5 dB over the given + /// brackets (dB re exact_msg_db) on `probe_s` probes. + struct spectral_rows { + std::vector open; + std::vector> value; ///< [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(value[row][s] - open[s]); + } + return a; + } + }; + + spectral_rows run_spectral(const std::string& room, double probe_s, const std::vector& shapes, bool chain, + double lo, double hi, unsigned sets = mutap_test::k_claim_seed_sets) { + const auto path = kk::room(room); + spectral_rows r; + r.open.assign(sets, 0.0); + r.value.assign(shapes.size(), std::vector(sets, 0.0)); + std::vector> jobs; + 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 = kk::k_s1; + p.probe_blocks = kk::probe_blocks(probe_s); + r.open[set] = + kk::open_loop_db(path, kk::k_s1, kk::near_end(kk::material::held, p.probe_blocks, seed + 10), p); + }); + for (size_t k = 0; k < shapes.size(); ++k) { + jobs.emplace_back([&, set, seed, k] { + kk::protocol p; + p.delay = kk::k_s1; + p.probe_blocks = kk::probe_blocks(probe_s); + p.chain_lo = lo; + p.chain_hi = hi; + sp::params prm; + prm.rt60 = 1.0; + prm.wet = 0.30; + prm.shape_max = shapes[k]; + const auto run = sp::measure(path, kk::material::held, p, seed, prm, chain, !chain, lo, hi); + r.value[k][set] = chain ? run.chain_db : run.open_rev_db; + }); + } + } + run_parallel(jobs); + return r; + } + + double mean(const std::vector& v) { + double sum = 0.0; + for (const double x : v) { + sum += x; + } + return v.empty() ? 0.0 : sum / static_cast(v.size()); + } + + /// The rows the default run cannot afford (ASan legs): MUTAP_SLOW=1. + bool slow_enabled() { + const char* v = std::getenv("MUTAP_SLOW"); + return v != nullptr && std::string(v) == "1"; + } + + /// The canceller alone's limit (karaoke_asg.h's measure, absolute dB) on + /// the same room and material, per seed set, bisected to 0.5 dB. + std::vector canceller_alone(const std::string& room, double probe_s) { + const auto path = kk::room(room); + std::vector a(mutap_test::k_claim_seed_sets, 0.0); + std::vector> jobs; + for (unsigned set = 0; set < a.size(); ++set) { + jobs.emplace_back([&, set] { + kk::protocol p; + p.delay = kk::k_s1; + p.probe_blocks = kk::probe_blocks(probe_s); + p.chain_lo = 2.0; + p.chain_hi = 20.0; + kk::setup s; + s.mat = kk::material::held; + a[set] = kk::measure(path, s, p, seed_in_set(2, set), false).chain_db; + }); + } + run_parallel(jobs); + return a; + } + +} // namespace + +// What the spectral reverb is, beside the plate it is compared with: its +// all-wet flat impulse response's T30 (Schroeder, -5 to -35 dB; reverb_rig.h's +// t30_seconds) at rt60 1 s and at the plate's 1.1846 s, and the stage's level +// change for white input (spectral_rig.h's level_db), flat, at wet 0.15 and +// 0.30. Measured (macOS x86_64, AppleClang 17, Release): T30 1.0000 s and +// 1.1846 s; level +3.40 / +7.60 dB at rt60 1 s and +3.81 / +8.22 dB at +// 1.1846 s (wet 0.15 / 0.30), where the plate's mix at decay 0.5 reads +// -0.956 / -1.687 dB: at equal wet the spectral reverb is the louder by 4.4 +// to 9.9 dB. Deterministic; the allowances are for other hosts' libm. +TEST(SpectralReverbHost, DecayAndLevelBesideThePlate) { + for (const double rt : {1.0, k_plate_t30}) { + const auto h = sp::wet_impulse_response(rt, 4.0 * rt); + const double t30 = rv::t30_seconds(h); + std::printf("spectral reverb rt60 %.4f s: T30 %.4f s; level", rt, t30); + EXPECT_NEAR(t30, rt, 0.005) << "rt60 " << rt; + for (const double w : {0.15, 0.30}) { + sp::params p; + p.rt60 = rt; + p.wet = w; + const double lvl = sp::level_db(sp::reverb(sp::config_of(p))); + std::printf(" w%.2f %+.2f dB", w, lvl); + EXPECT_GT(lvl, 0.0) << "the flat spectral reverb is louder than the dry bus"; + } + std::printf("\n"); + } + const auto plate = rv::plate_ir(rv::variant::shipped, 0.5, 0.0005); + EXPECT_NEAR(rv::t30_seconds(plate.l), k_plate_t30, 0.01); +} + +// WITH NO CANCELLER THE SPECTRAL REVERB IS UNSAFE (HANDOFF item 11's +// claim, re-measured with converged probes): the open loop with the reverb +// in it, against the dry open loop, held note, S1, wet 0.30, rt60 1 s, 80 s +// probes (SpectralReverbSweep.ProbeConvergence: the flat and shape_max-1 +// open rows converged by 80 s in the cabin and mt5). Shaped rows are shaped +// from the F_hat a canceller converged to in the same loop (spectral_rig.h) +// and then run without it. Measured (macOS x86_64, AppleClang 17, Release; +// per seed set 2 / 22 / 42 / 62 / 82, open-loop ASG, dB): +// +// cabin flat -30.51 -30.51 -30.51 -30.51 -30.51 median -30.51 +// cabin shape_max 1 -20.66 -19.61 -16.45 -20.31 -16.80 median -19.61 +// mt5 flat -22.77 -22.77 -22.77 -22.77 -22.77 median -22.77 +// mt5 shape_max 1 -13.98 -14.69 -15.04 -15.39 -14.34 median -14.69 +// +// (The flat rows repeat to the digit across seed sets: the held note's +// period is the same in every set and the open loop has no canceller to +// differ by.) Gates: flat median < -15 dB (largest measured -22.77), +// shape_max 1 median < -8 dB (largest -14.69). The plate costs a bare loop at +// most +5.79 dB by the magnitude bound (docs/reverb-afc.md). 189.41 s on 3 +// threads of the Intel Mac, so MUTAP_SLOW only; the default run carries +// UnsafeWithoutTheCancellerShortProbe below. +TEST(SpectralReverbHost, UnsafeWithoutTheCanceller) { + if (!slow_enabled()) { + GTEST_SKIP() << "set MUTAP_SLOW=1 (189 s on 3 threads); UnsafeWithoutTheCancellerShortProbe runs by default"; + } + for (const std::string room : {"cabin", "mt5"}) { + const auto r = run_spectral(room, 80.0, {0.0, 1.0}, false, -40.0, 5.0); + for (size_t k = 0; k < 2; ++k) { + const auto a = r.asg(k); + std::printf("%s open loop, %s: ASG%s median %+.2f\n", room.c_str(), k == 0 ? "flat" : "shape_max 1", + per_seed(a).c_str(), median(a)); + } + EXPECT_LT(median(r.asg(0)), -15.0) << room << " flat (measured -30.51 cabin, -22.77 mt5)"; + EXPECT_LT(median(r.asg(1)), -8.0) << room << " shape_max 1 (measured -19.61 cabin, -14.69 mt5)"; + } +} + +// BEHIND THE CANCELLER THE SPECTRAL REVERB COSTS MORE THAN THE PLATE +// (runaway, held note, cabin, S1, wet 0.30, rt60 1 s). The cost is the +// canceller-alone limit minus the limit with the reverb, per seed set; the +// spectral rows bisect on 80 s probes (ProbeConvergence: the cabin's flat, +// shape_max 1 and 2 chain rows converged at 80 s; shape_max 4 had not by +// 160 s and is not gated), the canceller alone on 20 s (converged at 10 s, +// its 20 to 80 s medians identical). ReverbStage.CancellerHoldsBehindThePlates +// bounds the plate's cost at +3.5 dB (measured <= +0.56 there, the as-shipped +// plate at decay 0.5 and 0.7, w 0.30, cabin S1). Measured (macOS x86_64, +// AppleClang 17, Release; per seed set, dB): +// +// ASG cost +// flat -5.94 -7.81 -5.62 -6.56 -5.00 +13.56 +21.62 +17.19 +18.69 +15.44 median +17.19 +// shape_max 1 -0.62 -0.62 +3.12 -0.94 +4.38 +8.25 +14.44 +8.44 +13.06 +6.06 median +8.44 +// +// against the canceller alone's median ASG of +11.74 there (ProbeConvergence, +// 20 to 80 s; its limits here, absolute dB: +0.73 +6.92 +4.67 +5.23 +3.54 +// at 20 s). Gates: the +// median cost of flat > +8 dB and of shape_max 1 > +4 dB, both above the +// plate's +3.5 dB ceiling (smallest measured +8.44; smallest single seed +// +6.06). 696.95 s on 3 threads of the Intel Mac, so MUTAP_SLOW only; the +// default run carries CostsMoreThanThePlateShortProbe below. +TEST(SpectralReverbHost, CostsMoreThanThePlateBehindTheCanceller) { + if (!slow_enabled()) { + GTEST_SKIP() << "set MUTAP_SLOW=1 (697 s on 3 threads); CostsMoreThanThePlateShortProbe runs by default"; + } + const auto canc = canceller_alone("cabin", 20.0); + const auto r = run_spectral("cabin", 80.0, {0.0, 1.0}, true, -20.0, 20.0); + std::printf("cabin canceller alone: limit%s\n", per_seed(canc).c_str()); + for (size_t k = 0; k < 2; ++k) { + std::vector cost; + for (size_t s = 0; s < canc.size(); ++s) { + cost.push_back(canc[s] - r.value[k][s]); + } + std::printf("cabin chain, %s: ASG%s median %+.2f | cost%s median %+.2f\n", k == 0 ? "flat" : "shape_max 1", + per_seed(r.asg(k)).c_str(), median(r.asg(k)), per_seed(cost).c_str(), median(cost)); + EXPECT_GT(median(cost), k == 0 ? 8.0 : 4.0) + << (k == 0 ? "flat (measured +17.19)" : "shape_max 1 (measured +8.44)"); + } +} + +// THE DEFAULT RUN'S PROXY FOR UnsafeWithoutTheCanceller (MUTAP_SLOW): the same +// open loops on 20 s probes and two seed sets (2, 22), cabin and mt5. Short +// probes UNDERSTATE how unsafe the open loop is (ProbeConvergence: every +// open row's median falls with probe length, e.g. cabin flat -29.75 at 20 s +// against -30.43 at 160 s), so a pass here is conservative. The gate is the +// mean of the two seed sets. Measured (macOS x86_64, AppleClang 17, Release; +// open-loop ASG, seed sets 2 / 22, dB; the 80 s medians from the slow row): +// +// 20 s, sets 2 / 22 mean 80 s median (5 sets) +// cabin flat -29.80 -29.80 -29.80 -30.51 +// cabin shape_max 1 -20.31 -19.26 -19.79 -19.61 +// mt5 flat -22.42 -22.42 -22.42 -22.77 +// mt5 shape_max 1 -13.63 -14.34 -13.98 -14.69 +// +// Gates, the slow row's: flat mean < -15 dB (largest measured -22.42), +// shape_max 1 mean < -8 dB (largest -13.98). 27.03 s on 3 threads of the +// Intel Mac. (Other hosts: not yet recorded; CI's legs run it.) +TEST(SpectralReverbHost, UnsafeWithoutTheCancellerShortProbe) { + for (const std::string room : {"cabin", "mt5"}) { + const auto r = run_spectral(room, 20.0, {0.0, 1.0}, false, -40.0, 5.0, 2); + for (size_t k = 0; k < 2; ++k) { + const auto a = r.asg(k); + std::printf("%s open loop, 20 s, %s: ASG%s mean %+.2f\n", room.c_str(), k == 0 ? "flat" : "shape_max 1", + per_seed(a).c_str(), mean(a)); + } + EXPECT_LT(mean(r.asg(0)), -15.0) << room << " flat (measured -29.80 cabin, -22.42 mt5)"; + EXPECT_LT(mean(r.asg(1)), -8.0) << room << " shape_max 1 (measured -19.79 cabin, -13.98 mt5)"; + } +} + +// THE DEFAULT RUN'S PROXY FOR CostsMoreThanThePlateBehindTheCanceller +// (MUTAP_SLOW): direction only, cabin, S1, wet 0.30, two seed sets (2, 22), +// 20 s probes. Per seed set, the limit of the canceller behind the +// as-shipped plate (decay 0.5) minus its limit behind the spectral reverb +// (rt60 1 s; flat and shape_max 1): positive = the spectral reverb costs more +// than the plate (the canceller-alone limit cancels out of the difference). +// Short probes UNDERSTATE the spectral reverb's cost (ProbeConvergence: the +// cabin's flat and shape_max 1 chain medians fall from 20 to 160 s, -4.67 +// to -6.82 and +2.27 to -0.27, where the plate's converged at 10 s), so a +// pass is conservative. The gate is the mean of the two seed sets. +// Measured (macOS x86_64, AppleClang 17, Release; dB): +// +// limits re exact_msg_db, plate - spectral, +// sets 2 / 22 sets 2 / 22 mean +// plate 0.5 +11.88 +9.69 +// flat -3.88 -3.88 +15.75 +13.56 +14.66 +// shape_max 1 +0.25 +1.75 +11.62 +7.94 +9.78 +// +// beside the slow row's 80 s costs against the canceller alone, +17.19 +// (flat) and +8.44 (shape_max 1), medians of five sets. Gates: the mean of +// plate - spectral > +6 dB for flat and > +3 dB for shape_max 1 (measured +// +14.66 and +9.78; smallest single set +7.94). 77.25 s on 3 threads of the +// Intel Mac. (Other hosts: not yet recorded; CI's legs run it.) +TEST(SpectralReverbHost, CostsMoreThanThePlateShortProbe) { + // One pool for the six bisections (plate, flat, shape_max 1 x two seed + // sets), each to 0.5 dB over a bracket around the 20 s limits measured + // (dB re exact_msg_db): the plate over [+4, +18], the spectral rows over + // [-14, +10]. + const auto path = kk::room("cabin"); + constexpr unsigned k_sets = 2; + std::vector plate(k_sets, 0.0); + std::vector> spec(2, std::vector(k_sets, 0.0)); + std::vector> jobs; + auto proto = [] { + kk::protocol p; + p.delay = kk::k_s1; + p.probe_blocks = kk::probe_blocks(20.0); + return p; + }; + for (unsigned set = 0; set < k_sets; ++set) { + const unsigned seed = seed_in_set(2, set); + jobs.emplace_back([&, set, seed] { + kk::protocol p = proto(); + p.chain_lo = 4.0; + p.chain_hi = 18.0; + rv::params pp; + pp.decay = 0.5; + pp.wet = 0.30; + rv::rig rig(pp); + kk::setup st; + st.mat = kk::material::held; + st.stage = mutap_test::forward_stage::of(&rig); + plate[set] = kk::measure(path, st, p, seed, false).chain_db; + }); + for (size_t k = 0; k < 2; ++k) { + jobs.emplace_back([&, set, seed, k] { + kk::protocol p = proto(); + p.chain_lo = -14.0; + p.chain_hi = 10.0; + sp::params prm; + prm.rt60 = 1.0; + prm.wet = 0.30; + prm.shape_max = (k == 0) ? 0.0 : 1.0; + spec[k][set] = sp::measure(path, kk::material::held, p, seed, prm, true, false).chain_db; + }); + } + } + run_parallel(jobs); + const double exact = mutap_test::exact_msg_db(path, kk::k_s1); + std::printf("cabin plate 0.5 w0.30, 20 s: limit re exact_msg_db %+.2f %+.2f\n", plate[0] - exact, plate[1] - exact); + for (size_t k = 0; k < 2; ++k) { + std::vector diff; + for (unsigned set = 0; set < k_sets; ++set) { + diff.push_back(plate[set] - spec[k][set]); + } + std::printf("cabin chain, 20 s, %s: limit re exact_msg_db %+.2f %+.2f | plate - spectral%s mean %+.2f\n", + k == 0 ? "flat" : "shape_max 1", spec[k][0] - exact, spec[k][1] - exact, per_seed(diff).c_str(), + mean(diff)); + EXPECT_GT(mean(diff), k == 0 ? 6.0 : 3.0) + << (k == 0 ? "flat (measured +14.66)" : "shape_max 1 (measured +9.78)"); + } +} diff --git a/tests/test_spectral_reverb_sweep.cpp b/tests/test_spectral_reverb_sweep.cpp new file mode 100644 index 0000000..d91b664 --- /dev/null +++ b/tests/test_spectral_reverb_sweep.cpp @@ -0,0 +1,893 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors +// +// The spectral reverb behind the canceller (docs/reverb-afc.md, "Spectral +// reverb"; HANDOFF item 11). 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; each finished bisection +// prints a "job" line (the checkpoint of a run that takes hours), every table +// the doc quotes is printed, and each row's median is recorded with +// RecordProperty. +// +// MUTAP_SLOW=1 MUTAP_SLOW_THREADS=3 build/tests/mutap_tests --gtest_filter='SpectralReverbSweep.*' +// +// Measurement: support/spectral_rig.h (the reverb suites' protocol, +// support/karaoke_asg.h, with tap::mu::spectral_reverb as decorrelated_loop's +// forward stage: converged 1500 blocks at exact_msg_db - 6 with the reverb +// flat in the loop, shaped once from the converged F_hat, the shape held +// through the bisection; its one block of latency inside the forward delay). +// Bisection to 0.1 dB over [exact - 20, exact + 30] with the canceller and +// [exact - 40, exact + 10] without; seed sets 2 / 22 / 42 / 62 / 82; six +// band-limited rooms of both generator families. Like for like with the +// plates (test_reverb_stage_sweep.cpp): +// +// ASG = the chain's limit - the dry open loop's (no canceller, no +// reverb), same probe material and length; +// cost = the canceller-alone limit - the chain's, per seed; +// open cost = the dry open loop's limit - the open loop's with the +// reverb (no canceller), per seed: HANDOFF item 11's +// "flat is unsafe open loop"; +// level = the stage's broadband level change for white input (the +// reverb as shaped in that run); matched = cost - level. +// +// ProbeConvergence cabin and mt5, S1, wet 0.30, rt60 1 s, flat and shaped +// (shape_max 1, 2, 4), held note and speech envelope, +// probes 10 / 20 / 40 / 80 / 160 s, with the open-loop +// limit of the same reverb; the canceller alone to 80 s. +// "Converged" as in the plates' sweep: the shortest probe +// from which every longer probe's five-seed median is +// within 0.25 dB of the longest's. +// Grid six rooms, S1 and S3: the held note, flat and shaped +// (1, 2, 4) at wet 0.15 / 0.30 and rt60 1 s and 1.1846 s +// (the as-shipped plate's decay-0.5 T30, L), beside the +// as-shipped plate at decay 0.5 and the same wets; the +// speech envelope at wet 0.30, rt60 1.1846 s. Probes from +// ProbeConvergence (grid_probe). Prints the cost, open +// cost and the hypothesis statistics (spectral_rig.h). +// +// Environment: SPECTRAL_SWEEP_PROBES ("10,20"), SPECTRAL_SWEEP_MATS ("held", +// "speech" or both), SPECTRAL_SWEEP_SHAPES ("0,1,2,4"; 0 = flat), +// SPECTRAL_SWEEP_DELAYS ("480,960"), SPECTRAL_SWEEP_RT60S, SPECTRAL_SWEEP_WETS, +// SPECTRAL_SWEEP_ROOMS, SPECTRAL_SWEEP_PROBE_HELD / _SPEECH, +// SPECTRAL_SWEEP_PLATES=0 (skip the plate rows). + +#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" +#include "support/spectral_rig.h" + +namespace { + + using mutap_test::median; + using mutap_test::seed_in_set; + namespace kk = mutap_test::karaoke; + namespace rv = mutap_test::reverb; + namespace sp = mutap_test::spectral; + + /// The as-shipped plate's decay-0.5 T30 (L return, damping 0.0005; + /// ReverbStage.PlateDecayAndMixLevel): the matched rt60. + constexpr double k_plate_t30 = 1.1846; + + bool slow_enabled() { + const char* v = std::getenv("MUTAP_SLOW"); + return v != nullptr && std::string(v) == "1"; + } + + 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) { + if (std::isnan(v)) { + return " n/a"; + } + char b[16]; + std::snprintf(b, sizeof b, "%+6.2f", v); + return b; + } + + std::string fmt3(double v) { + if (std::isnan(v)) { + return " n/a"; + } + char b[16]; + std::snprintf(b, sizeof b, "%+.3f", v); + return b; + } + + const char* mat_name(kk::material m) { + return m == kk::material::held ? "held" : "speech"; + } + + /// Comma-separated numbers from the environment, or the default. + std::vector list_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; + } + + std::vector mats_from_env(std::vector dflt) { + const char* v = std::getenv("SPECTRAL_SWEEP_MATS"); + if (v == nullptr) { + return dflt; + } + const std::string s = v; + std::vector out; + if (s.find("held") != std::string::npos) { + out.push_back(kk::material::held); + } + if (s.find("speech") != std::string::npos) { + out.push_back(kk::material::speech); + } + return out; + } + + std::vector rooms_from_env(std::vector dflt) { + const char* v = std::getenv("SPECTRAL_SWEEP_ROOMS"); + 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(s.substr(a, b == std::string::npos ? std::string::npos : b - a)); + if (b == std::string::npos) { + break; + } + a = b + 1; + } + return out; + } + + std::string spectral_label(const sp::params& p) { + char b[96]; + if (p.shape_max > 0.0) { + std::snprintf(b, sizeof b, "spectral shaped %g rt%.4g w%.2f", p.shape_max, p.rt60, p.wet); + } + else { + std::snprintf(b, sizeof b, "spectral flat rt%.4g w%.2f", p.rt60, p.wet); + } + return b; + } + + std::string plate_label(const rv::params& p) { + char b[96]; + std::snprintf(b, sizeof b, "plate %s d%.2f w%.2f", rv::name(p.plate), p.decay, p.wet); + return b; + } + + enum class kind { open, canceller, plate, spectral }; + + /// One bisection (a spectral job is the chain bisection and, with + /// open_reverb, the open loop with the reverb). + 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 plate; + sp::params spec; + bool chain = true; ///< spectral: bisect the chain + bool open_reverb = true; ///< spectral: bisect the open loop with the reverb + // results + double value = 0.0; ///< the bisected chain (or open) limit, dB + double exact = 0.0; + sp::run run; + double seconds = 0.0; + }; + + std::string job_label(const job& j) { + switch (j.what) { + case kind::open: + return "open"; + case kind::canceller: + return "canceller"; + case kind::plate: + return plate_label(j.plate); + case kind::spectral: + default: + return spectral_label(j.spec); + } + } + + 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; + switch (j.what) { + case kind::open: + j.value = kk::open_loop_db(path, j.delay, kk::near_end(j.mat, p.probe_blocks, j.seed + 10), p); + break; + case kind::canceller: + j.value = kk::measure(path, s, p, j.seed, false).chain_db; + break; + case kind::plate: { + rv::rig rig(j.plate); + s.stage = mutap_test::forward_stage::of(&rig); + j.value = kk::measure(path, s, p, j.seed, false).chain_db; + break; + } + case kind::spectral: + default: + j.run = sp::measure(path, j.mat, p, j.seed, j.spec, j.chain, j.open_reverb); + j.value = j.chain ? j.run.chain_db : j.run.open_rev_db; + break; + } + 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", j.room.c_str(), j.delay, mat_name(j.mat), + j.probe_s, j.seed, job_label(j).c_str(), j.value - j.exact); + if (j.what == kind::spectral) { + std::printf(" open %+.2f level %+.2f rhoT %+.3f share %.3f mis %+.2f", + j.run.open_rev_db - j.exact, j.run.level, j.run.bins.rho_t_mis, + j.run.bins.share_top, j.run.bins.misalign_db); + } + std::printf(" (%.0f s)\n", 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 jobs, %.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 jobs 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; + using seed_key = std::tuple; + + struct results { + std::map> asg; ///< per seed set, seed order + std::map> cost; ///< canceller - this + std::map> open_asg; ///< open with reverb - dry open (spectral) + std::map> level; ///< spectral + std::map> runs; ///< spectral jobs, seed order + std::map open; + std::map canc; + std::map, double> chain; ///< (row, seed) -> limit + }; + + /// Enqueue the references (dry open loop, canceller alone) once per + /// (room, delay, material, probe). + using ref_map = std::map, unsigned>; + + void enqueue_refs(std::vector& jobs, ref_map& refs, const std::string& room, size_t delay, kk::material mat, + double probe_s, bool canceller) { + const auto ref = std::make_tuple(room, delay, static_cast(mat), probe_s); + const unsigned have = refs[ref]; + const unsigned want = canceller ? 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 ((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; + } + + void enqueue_row(std::vector& jobs, const std::string& room, size_t delay, kk::material mat, double probe_s, + const job& proto) { + for (unsigned set = 0; set < mutap_test::k_claim_seed_sets; ++set) { + job j = proto; + j.room = room; + j.delay = delay; + j.mat = mat; + j.probe_s = probe_s; + j.seed = seed_in_set(2, set); + jobs.push_back(j); + } + } + + results collect(const std::vector& jobs) { + results r; + for (const auto& j : jobs) { + const seed_key k{j.room, j.delay, static_cast(j.mat), j.probe_s, j.seed}; + if (j.what == kind::open) { + r.open[k] = j.value; + } + else if (j.what == kind::canceller) { + r.canc[k] = j.value; + } + } + 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 seed_key s{j->room, j->delay, static_cast(j->mat), j->probe_s, j->seed}; + const auto o = r.open.find(s); + if (o == r.open.end()) { + continue; + } + if (j->what == kind::spectral && !j->chain) { // an open-loop-only job + r.open_asg[k].push_back(j->run.open_rev_db - o->second); + continue; // its shape repeats the chain job's (same seed, same convergence) + } + r.asg[k].push_back(j->value - o->second); + r.chain[{k, j->seed}] = j->value; + const auto c = r.canc.find(s); + if (j->what != kind::canceller && c != r.canc.end()) { + r.cost[k].push_back(c->second - j->value); + } + if (j->what == kind::spectral) { + r.runs[k].push_back(j); + r.level[k].push_back(j->run.level); + if (!std::isnan(j->run.open_rev_db)) { + r.open_asg[k].push_back(j->run.open_rev_db - o->second); + } + } + } + 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)); + } + + 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 %-40s 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)); + } + const auto o = r.open_asg.find(k); + if (o != r.open_asg.end()) { + std::printf(" | open ASG%s | median %s", per_seed(o->second).c_str(), fmt(median(o->second)).c_str()); + record(suite, k, "open_asg", median(o->second)); + } + const auto l = r.level.find(k); + if (l != r.level.end()) { + std::printf(" | level median %s", fmt(median(l->second)).c_str()); + } + std::printf("\n"); + } + } + + /// 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 = std::nan(""); + double room_lo = std::nan(""); + double room_hi = std::nan(""); + double seed_lo = std::nan(""); + double seed_hi = std::nan(""); + }; + + across summarize(const std::map>& m, const std::vector& rooms, + size_t delay, kk::material mat, double probe, const std::string& l) { + std::vector meds; + std::vector all; + for (const auto& room : rooms) { + const auto it = m.find(row_key{room, delay, static_cast(mat), probe, l}); + if (it == m.end() || it->second.empty()) { + 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; + } + + std::string across_str(const across& a) { + return fmt(a.median_of_medians) + " [" + fmt(a.room_lo) + ", " + fmt(a.room_hi) + "] {" + fmt(a.seed_lo) + ", " + + fmt(a.seed_hi) + "}"; + } + + /// The converged-probe summary: each row's median by probe and the + /// shortest converged probe (0.25 dB of the longest). + void convergence_table(const std::map>& m, const std::string& room, kk::material mat, + const std::string& l, const std::vector& probes, const char* what) { + std::vector med; + std::vector used; + for (const double p : probes) { + const auto a = m.find(row_key{room, kk::k_s1, static_cast(mat), p, l}); + if (a == m.end() || a->second.empty()) { + continue; + } + med.push_back(median(a->second)); + used.push_back(p); + } + if (med.empty()) { + return; + } + size_t conv = med.size(); + for (size_t i = med.size(); i-- > 0;) { + if (std::abs(med[i] - med.back()) <= 0.25) { + conv = i; + } + else { + break; + } + } + const bool ok = conv + 1 < med.size(); + std::printf("%-6s %-6s %-9s %-40s", room.c_str(), mat_name(mat), what, l.c_str()); + for (size_t i = 0; i < probes.size(); ++i) { + std::printf(" %8s", i < med.size() ? fmt(med[i]).c_str() : ""); + } + if (ok) { + std::printf(" %5.0f s\n", used[conv]); + } + else { + std::printf(" not converged\n"); + } + char name[192]; + std::snprintf(name, sizeof name, "conv_%s_%s_%s_%s_converged_s", room.c_str(), mat_name(mat), what, 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(used[conv]) : std::string("none")); + } + + /// The plate mix's level change for white input (reverb_rig.h), dB. + double plate_level_db(const rv::params& p) { + const auto h = rv::plate_ir(p.plate, p.decay, p.damping); + return rv::mix_energy_db(rv::mono_return(h, p.ret), p.wet); + } + + /// The probe the open-loop rows (no canceller, the reverb in the loop) + /// run with (ProbeConvergence; SPECTRAL_SWEEP_PROBE_OPEN overrides). + double open_probe(kk::material /*mat*/) { + if (const char* e = std::getenv("SPECTRAL_SWEEP_PROBE_OPEN")) { + return std::stod(e); + } + return 160.0; + } + + /// The probe each material's grid rows run with (ProbeConvergence; + /// SPECTRAL_SWEEP_PROBE_HELD / _SPEECH override). On the held note the + /// spectral chain rows converged at 80 s where they converged at all + /// (cabin flat, shape_max 1 and 2; the rest not by 160 s), and the 40 s + /// median ASG read at or above the 160 s one in 7 of the 8 cabin and mt5 + /// rows (mt5 shape_max 2: +5.98 against +6.46); the plate and canceller + /// rows converged at 10 s. 40 s is the budget: the held grid is not + /// converged, mostly understates the spectral reverb's cost, and is + /// reported, not gated. + double grid_probe(kk::material mat) { + if (mat == kk::material::held) { + if (const char* e = std::getenv("SPECTRAL_SWEEP_PROBE_HELD")) { + return std::stod(e); + } + return 40.0; + } + if (const char* e = std::getenv("SPECTRAL_SWEEP_PROBE_SPEECH")) { + return std::stod(e); + } + return 80.0; + } + +} // namespace + +// Probe-length convergence of the spectral rows, the chain's and the open +// loop's: how long a probe the claims need. +TEST(SpectralReverbSweep, ProbeConvergence) { + if (!slow_enabled()) { + GTEST_SKIP() << "set MUTAP_SLOW=1 for the spectral reverb sweep"; + } + const auto probes = list_from_env("SPECTRAL_SWEEP_PROBES", {10.0, 20.0, 40.0, 80.0, 160.0}); + const auto shapes = list_from_env("SPECTRAL_SWEEP_SHAPES", {0.0, 1.0, 2.0, 4.0}); + const auto mats = mats_from_env({kk::material::held, kk::material::speech}); + const auto rooms = rooms_from_env({"cabin", "mt5"}); + const double canceller_longest = 80.0; + + std::vector jobs; + ref_map refs; + for (const auto& room : rooms) { + for (const auto mat : mats) { + for (const double probe : probes) { + enqueue_refs(jobs, refs, room, kk::k_s1, mat, probe, probe <= canceller_longest); + for (const double sm : shapes) { + job proto; + proto.what = kind::spectral; + proto.spec.rt60 = 1.0; + proto.spec.wet = 0.30; + proto.spec.shape_max = sm; + enqueue_row(jobs, room, kk::k_s1, mat, probe, proto); + } + } + } + } + run_jobs(jobs); + const auto r = collect(jobs); + std::printf("\nProbe convergence, S1, wet 0.30, rt60 1 s; per seed set (2 22 42 62 82), dB\n"); + print_rows(r, "sconv"); + + std::printf("\nMedian ASG by probe (dB); chain = with the canceller, open = the open loop with the reverb (no " + "canceller), both re the dry open loop; converged = the shortest probe from which every longer " + "probe's median is within 0.25 dB of the longest's\n"); + std::printf("%-6s %-6s %-9s %-40s", "room", "mat", "", "row"); + for (const double p : probes) { + std::printf(" %6.0f s", p); + } + std::printf(" converged\n"); + for (const auto& room : rooms) { + for (const auto mat : mats) { + convergence_table(r.asg, room, mat, "canceller", probes, "chain"); + for (const double sm : shapes) { + sp::params p; + p.rt60 = 1.0; + p.wet = 0.30; + p.shape_max = sm; + convergence_table(r.asg, room, mat, spectral_label(p), probes, "chain"); + convergence_table(r.open_asg, room, mat, spectral_label(p), probes, "open"); + convergence_table(r.cost, room, mat, spectral_label(p), probes, "cost"); + } + } + } + for (const auto& [k, v] : r.asg) { + for (const double x : v) { + EXPECT_TRUE(std::isfinite(x)); + } + } +} + +// The grid: the spectral reverb against the plate behind the canceller and +// with no canceller, and the per-bin hypothesis, six rooms. +TEST(SpectralReverbSweep, Grid) { + if (!slow_enabled()) { + GTEST_SKIP() << "set MUTAP_SLOW=1 for the spectral reverb sweep"; + } + const auto rooms = rooms_from_env({"cabin", "studio", "rehearsal", "hall", "mt5", "mt9"}); + const auto shapes = list_from_env("SPECTRAL_SWEEP_SHAPES", {0.0, 1.0, 2.0, 4.0}); + const auto wets = list_from_env("SPECTRAL_SWEEP_WETS", {0.15, 0.30}); + const auto rt60s = list_from_env("SPECTRAL_SWEEP_RT60S", {1.0, k_plate_t30}); + const auto delays = list_from_env("SPECTRAL_SWEEP_DELAYS", {480.0, 960.0}); + const auto mats = mats_from_env({kk::material::held, kk::material::speech}); + const char* pe = std::getenv("SPECTRAL_SWEEP_PLATES"); + const bool plates = pe == nullptr || std::string(pe) != "0"; + + // The rows (the budget: ProbeConvergence put the spectral rows at 80 s + // and beyond, against the plates' 10 s). Held: every shape at both wets + // at the matched rt60 (the plate's decay-0.5 T30), and at rt60 1 s at + // wet 0.30 (ProbeConvergence's configuration), the plate at each wet; + // the open loop (no canceller) with each spectral row at the matched + // rt60. Speech, S1 only: wet 0.30, rt60 1 s (ProbeConvergence's), and the + // plate there; no open loop. + struct row { + job proto; + kk::material mat; + }; + std::vector rows; + for (const auto mat : mats) { + const bool held = mat == kk::material::held; + for (const double w : wets) { + if (!held && w != 0.30) { + continue; + } + if (plates) { + job pj; + pj.what = kind::plate; + pj.plate.plate = rv::variant::shipped; + pj.plate.decay = 0.5; + pj.plate.wet = w; + rows.push_back({pj, mat}); + } + for (const double rt : rt60s) { + if (held ? (rt != k_plate_t30 && w != 0.30) : rt != 1.0) { + continue; + } + for (const double sm : shapes) { + job sj; + sj.what = kind::spectral; + sj.spec.rt60 = rt; + sj.spec.wet = w; + sj.spec.shape_max = sm; + rows.push_back({sj, mat}); + } + } + } + } + std::vector jobs; + ref_map refs; + for (const auto& room : rooms) { + for (const double dd : delays) { + const auto d = static_cast(dd); + for (const auto& rw : rows) { + if (rw.mat == kk::material::speech && d != kk::k_s1) { + continue; + } + enqueue_refs(jobs, refs, room, d, rw.mat, grid_probe(rw.mat), true); + job chain_job = rw.proto; + chain_job.open_reverb = false; + enqueue_row(jobs, room, d, rw.mat, grid_probe(rw.mat), chain_job); + if (rw.proto.what == kind::spectral && rw.mat == kk::material::held + && rw.proto.spec.rt60 == k_plate_t30) { + // The open loop with the same reverb, on its own probe. + enqueue_refs(jobs, refs, room, d, rw.mat, open_probe(rw.mat), false); + job open_job = rw.proto; + open_job.chain = false; + enqueue_row(jobs, room, d, rw.mat, open_probe(rw.mat), open_job); + } + } + } + } + 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, "sgrid"); + + std::printf("\nGrid summary: median over the rooms of each room's five-seed median [per-room range] {per-seed " + "range}; cost = canceller alone - chain; open = the open loop with the reverb (no canceller) re the " + "dry open loop; level = the stage's level change for white input (median over runs); matched = " + "cost - level\n"); + for (const double dd : delays) { + const auto d = static_cast(dd); + for (const auto mat : mats) { + const double probe = grid_probe(mat); + const auto c = summarize(r.asg, rooms, d, mat, probe, "canceller"); + if (!c.found) { + continue; + } + std::printf("\n%s, %s, %g s probes. Canceller alone: ASG %s\n", d == kk::k_s1 ? "S1 (10 ms)" : "S3 (20 ms)", + mat_name(mat), probe, across_str(c).c_str()); + std::printf(" %-40s %-34s %-34s %-34s %7s %8s\n", "row", "ASG", "cost", "open ASG", "level", "matched"); + for (const auto& rw : rows) { + if (rw.mat != mat) { + continue; + } + const std::string l = job_label(rw.proto); + const auto a = summarize(r.asg, rooms, d, mat, probe, l); + const auto co = summarize(r.cost, rooms, d, mat, probe, l); + const auto o = summarize(r.open_asg, rooms, d, mat, open_probe(mat), l); + auto lv = summarize(r.level, rooms, d, mat, probe, l); + if (rw.proto.what == kind::plate) { + lv.median_of_medians = plate_level_db(rw.proto.plate); + } + if (!a.found) { + continue; + } + std::printf(" %-40s %-34s %-34s %-34s %7s %8s\n", l.c_str(), across_str(a).c_str(), + across_str(co).c_str(), o.found ? across_str(o).c_str() : "", + fmt(lv.median_of_medians).c_str(), + fmt(co.median_of_medians - lv.median_of_medians).c_str()); + } + } + } + + // Spectral against the plate at the same wet, per seed (chain limits), + // and shaped against flat. + std::printf("\nPer-seed differences of chain limits, median over the rooms of each room's median [per-room " + "range] (count of rooms where the room median > 0)\n"); + auto paired = [&](size_t d, kk::material mat, const std::string& a_label, const std::string& b_label) { + const double probe = grid_probe(mat); + std::vector meds; + for (const auto& room : rooms) { + const row_key ka{room, d, static_cast(mat), probe, a_label}; + const row_key kb{room, d, static_cast(mat), probe, b_label}; + std::vector diff; + for (unsigned set = 0; set < mutap_test::k_claim_seed_sets; ++set) { + const unsigned seed = seed_in_set(2, set); + const auto ia = r.chain.find({ka, seed}); + const auto ib = r.chain.find({kb, seed}); + if (ia != r.chain.end() && ib != r.chain.end()) { + diff.push_back(ia->second - ib->second); + } + } + if (!diff.empty()) { + meds.push_back(median(diff)); + } + } + if (meds.empty()) { + return; + } + const auto pos = std::count_if(meds.begin(), meds.end(), [](double x) { return x > 0.0; }); + std::printf(" %s %-6s %-40s - %-30s %s [%s, %s] (%ld of %zu)\n", d == kk::k_s1 ? "S1" : "S3", mat_name(mat), + a_label.c_str(), b_label.c_str(), fmt(median(meds)).c_str(), + fmt(*std::min_element(meds.begin(), meds.end())).c_str(), + fmt(*std::max_element(meds.begin(), meds.end())).c_str(), static_cast(pos), meds.size()); + }; + for (const double dd : delays) { + const auto d = static_cast(dd); + for (const auto& rw : rows) { + if (rw.proto.what != kind::spectral) { + continue; + } + const std::string l = job_label(rw.proto); + if (plates) { + rv::params pp; + pp.plate = rv::variant::shipped; + pp.decay = 0.5; + pp.wet = rw.proto.spec.wet; + paired(d, rw.mat, l, plate_label(pp)); + } + if (rw.proto.spec.shape_max > 0.0) { + sp::params flat = rw.proto.spec; + flat.shape_max = 0.0; + paired(d, rw.mat, l, spectral_label(flat)); + } + } + } + + // The hypothesis, per shaped row: medians over every run (rooms x seeds) + // of the per-run rank correlations and the worst-quarter energy share, + // and across runs, Spearman's rho of the share against the run's shaped + // - flat chain limit (same room, delay, seed). + std::printf("\nShaping hypothesis (spectral_rig.h): per run at the moment of shaping, rho = Spearman over the 65 " + "bins; medians over the runs [min, max]; flat share = 16/65 = 0.246\n"); + std::printf(" %-9s %-40s %-26s %-26s %-26s %-26s %-26s %-14s %-14s\n", "", "row", "rho(T60, |F-F^|)", + "rho(T60, |F-F^|/|F|)", "share(E, top quarter)", "peak m*g re flat, dB", "misalignment dB", + "rho(share, s-f)", "rho(peak, s-f)"); + for (const double dd : delays) { + const auto d = static_cast(dd); + for (const auto& rw : rows) { + if (rw.proto.what != kind::spectral || rw.proto.spec.shape_max <= 0.0) { + continue; + } + const std::string l = job_label(rw.proto); + const double probe = grid_probe(rw.mat); + sp::params flat = rw.proto.spec; + flat.shape_max = 0.0; + std::vector rt; + std::vector pk; + std::vector rr; + std::vector peak_x; + std::vector sh; + std::vector mi; + std::vector share_x; + std::vector gain_y; + for (const auto& room : rooms) { + const row_key k{room, d, static_cast(rw.mat), probe, l}; + const row_key kf{room, d, static_cast(rw.mat), probe, spectral_label(flat)}; + const auto it = r.runs.find(k); + if (it == r.runs.end()) { + continue; + } + for (const job* j : it->second) { + rt.push_back(j->run.bins.rho_t_mis); + pk.push_back(j->run.bins.peak_rel_db); + rr.push_back(j->run.bins.rho_t_rel); + sh.push_back(j->run.bins.share_top); + mi.push_back(j->run.bins.misalign_db); + const auto f = r.chain.find({kf, j->seed}); + if (f != r.chain.end()) { + share_x.push_back(j->run.bins.share_top); + peak_x.push_back(j->run.bins.peak_rel_db); + gain_y.push_back(j->value - f->second); + } + } + } + if (rt.empty()) { + continue; + } + auto mm = [](const std::vector& v) { + return fmt3(median(v)) + " [" + fmt3(*std::min_element(v.begin(), v.end())) + ", " + + fmt3(*std::max_element(v.begin(), v.end())) + "]"; + }; + const double rho_sg = share_x.size() > 2 ? sp::spearman(share_x, gain_y) : std::nan(""); + const double rho_pg = peak_x.size() > 2 ? sp::spearman(peak_x, gain_y) : std::nan(""); + std::printf(" %s %-6s %-40s %-26s %-26s %-26s %-26s %-26s %-14s %s (n %zu)\n", d == kk::k_s1 ? "S1" : "S3", + mat_name(rw.mat), l.c_str(), mm(rt).c_str(), mm(rr).c_str(), mm(sh).c_str(), mm(pk).c_str(), + mm(mi).c_str(), fmt3(rho_sg).c_str(), fmt3(rho_pg).c_str(), share_x.size()); + const row_key rk{"all", d, static_cast(rw.mat), probe, l}; + record("hyp", rk, "rho_t_mis", median(rt)); + record("hyp", rk, "share_top", median(sh)); + record("hyp", rk, "rho_share_gain", rho_sg); + record("hyp", rk, "peak_rel_db", median(pk)); + record("hyp", rk, "rho_peak_gain", rho_pg); + } + } + for (const auto& [k, v] : r.asg) { + for (const double x : v) { + EXPECT_TRUE(std::isfinite(x)); + } + } +} diff --git a/tools/notebook/karaoke_ramp_dump.cpp b/tools/notebook/karaoke_ramp_dump.cpp index f6c3f28..087a01a 100644 --- a/tools/notebook/karaoke_ramp_dump.cpp +++ b/tools/notebook/karaoke_ramp_dump.cpp @@ -42,6 +42,16 @@ // 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). // +// THE SPECTRAL REVERB (docs/reverb-afc.md, "Spectral reverb"; tests/support/ +// spectral_rig.h): `--reverb spectral` puts tap::mu::spectral_reverb in the +// same slot at `--wet` and `--rt60`, flat (`--shape-max 0`, the default) or +// shaped from the live canceller's F_hat with that ceiling, reshaped every +// `--reshape-blocks` blocks (256, 0.34 s: the branch measured the cadence +// not to matter) once F_hat is nonzero (the canceller starts cold, so the +// shape starts flat). Its one block of latency comes out of `--delay`, as in +// the suites. `--reverb-ir-only` writes its all-wet flat impulse response +// (twice, as L and R). +// // Outputs, with `--out PREFIX`: // PREFIX.wav float32, 48 kHz, 4 channels: c (the canceller output), // the voice stem (the near end v as it enters the mic), @@ -81,10 +91,12 @@ #include "fixtures/rir_studio.h" #include "mutap/fd_kalman.h" #include "mutap/pem_afc.h" +#include "mutap/spectral_reverb.h" #include "support/closed_loop.h" #include "support/decorrelated_loop.h" #include "support/reverb_rig.h" #include "support/rooms.h" +#include "support/spectral_rig.h" namespace { @@ -113,7 +125,7 @@ namespace { 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 - std::string reverb = "none"; ///< none | shipped | paper + std::string reverb = "none"; ///< none | shipped | paper | spectral double wet = 0.30; double decay = 0.5; double damping = 0.0005; @@ -121,6 +133,9 @@ namespace { 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 + double rt60 = 1.0; ///< spectral + double shape_max = 0.0; ///< spectral: 0 flat, else shaped from F_hat with this ceiling + size_t reshape = 256; ///< spectral: blocks between reshapes std::string out; }; @@ -132,7 +147,8 @@ namespace { " [--aux-db DB] [--seed 2] [--warmup 30] [--start-below 20] [--rate 0.5]\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", + " [--return L|M] [--topology bus|dry] [--lib-shift] [--reverb-ir-only]\n" + " [--reverb spectral [--rt60 1] [--shape-max 0] [--reshape-blocks 256]]\n", msg); std::exit(2); } @@ -214,6 +230,15 @@ namespace { else if (a == "--lib-shift") { o.lib_shift = true; } + else if (a == "--rt60") { + o.rt60 = std::stod(next()); + } + else if (a == "--shape-max") { + o.shape_max = std::stod(next()); + } + else if (a == "--reshape-blocks") { + o.reshape = static_cast(std::stoul(next())); + } else if (a == "--reverb-ir-only") { o.reverb_ir = true; } @@ -230,8 +255,14 @@ 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.reverb != "none" && o.reverb != "shipped" && o.reverb != "paper" && o.reverb != "spectral") { + usage("--reverb must be none, shipped, paper or spectral"); + } + if (o.reverb == "spectral" && (o.shift_hz != 0.0 || o.lib_shift)) { + usage("--reverb spectral takes no shift"); + } + if (o.reverb == "spectral" && o.reshape == 0) { + usage("--reshape-blocks must be >= 1"); } if (o.ret != "L" && o.ret != "M") { usage("--return must be L or M"); @@ -368,6 +399,19 @@ int main(int argc, char** argv) { } } + if (o.reverb_ir && o.reverb == "spectral") { + const auto h = mutap_test::spectral::wet_impulse_response(o.rt60, 4.0 * o.rt60); + std::vector> ch(2); + for (const double x : h) { + ch[0].push_back(static_cast(x)); + ch[1].push_back(static_cast(x)); + } + 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; + } if (o.reverb_ir) { const auto h = mutap_test::reverb::plate_ir(rig_params(o).plate, o.decay, o.damping); std::vector> ch(2); @@ -399,13 +443,23 @@ 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; - const bool staged = o.reverb != "none" || o.lib_shift; - std::unique_ptr rig; - if (staged) { + lc.feedback_path = path; + lc.block_size = k_block; + lc.forward_delay = o.delay; + lc.forward_gain_db = start_db; + const bool spectral = o.reverb == "spectral"; + const bool staged = o.reverb != "none" || o.lib_shift; + std::unique_ptr rig; + std::unique_ptr spec; + if (spectral) { + mutap_test::spectral::params sp; + sp.rt60 = o.rt60; + sp.wet = o.wet; + sp.shape_max = o.shape_max; + spec = std::make_unique(mutap_test::spectral::config_of(sp)); + lc.stage = mutap_test::forward_stage::of(spec.get()); + } + else if (staged) { rig = std::make_unique(rig_params(o)); lc.stage = mutap_test::forward_stage::of(rig.get()); } @@ -430,6 +484,9 @@ int main(int argc, char** argv) { } std::string csv = "time_s,gain_db\n"; + std::vector f_hat(afc.filter_length(), 0.0); + size_t reshapes = 0; + double runaway_db = NAN; double runaway_t = NAN; size_t stop_at = max_blocks; @@ -444,6 +501,10 @@ int main(int argc, char** argv) { std::snprintf(row, sizeof row, "%.6f,%.6f\n", t, g); csv += row; + if (spectral && o.canceller && o.shape_max > 0.0 && blk > 0 && blk % o.reshape == 0) { + afc.copy_impulse_response(f_hat.data()); + reshapes += spec->reshape_from_impulse_response(f_hat.data(), f_hat.size()) ? 1 : 0; + } 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(); @@ -475,12 +536,13 @@ int main(int argc, char** argv) { "\"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, \"reverb\": \"%s\", \"wet\": %g, " - "\"decay\": %g, \"damping\": %g, \"return\": \"%s\", \"topology\": \"%s\", \"staged\": %s}\n", + "\"decay\": %g, \"damping\": %g, \"return\": \"%s\", \"topology\": \"%s\", \"staged\": %s, " + "\"rt60\": %g, \"shape_max\": %g, \"reshape_blocks\": %zu, \"reshapes\": %zu}\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, o.reverb.c_str(), o.wet, o.decay, o.damping, o.ret.c_str(), - o.topology.c_str(), staged ? "true" : "false"); + o.topology.c_str(), staged ? "true" : "false", o.rt60, o.shape_max, o.reshape, reshapes); return 0; } diff --git a/tools/notebook/reverb_audible.py b/tools/notebook/reverb_audible.py index 8b33285..85c8e5e 100644 --- a/tools/notebook/reverb_audible.py +++ b/tools/notebook/reverb_audible.py @@ -35,6 +35,16 @@ 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. +With --plate spectral the reverb is tap::mu::spectral_reverb (docs/reverb-afc.md, +"Spectral reverb"; the dump's --reverb spectral) and the rows are + specS_R_W + the spectral reverb at rt60 R, wet W, flat (S = 0) or shaped + from the live canceller's F_hat with shape_max S (reshaped every + 256 blocks) +for S in --shape-maxes, R in --decays (read as rt60 values, seconds) and W in +--wets, beside dry and plain; no shift rows. Their --chain-rt60 is the +spectral reverb's own T30 at rt60 R (its all-wet flat impulse response through +schroeder()), for flat and shaped rows alike. + 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 @@ -85,14 +95,21 @@ import measure_rir # noqa: E402 -def configs(decays: list[float], wets: list[float], plate: str) -> dict: +def configs(decays: list[float], wets: list[float], plate: str, shape_maxes: list[float] | None = None) -> 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"), } + if plate == "spectral": + for sm in shape_maxes or [0.0]: + for rt in decays: + for w in wets: + out[f"spec{sm:g}_{rt:g}_{w:g}"] = (["--reverb", "spectral", "--rt60", f"{rt:g}", "--wet", f"{w:g}", + "--shape-max", f"{sm:g}"], "reverb") + return out + out["shift2"] = (["--lib-shift", "--shift-hz", "2"], "shift") + out["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}"] @@ -103,6 +120,22 @@ def configs(decays: list[float], wets: list[float], plate: str) -> dict: return out +def spectral_t30(dump: str, work: pathlib.Path, rt60: float) -> float: + """The spectral reverb's own T30 at rt60 (its all-wet flat impulse response).""" + prefix = work / f"spectral_rt{rt60:g}" + subprocess.run([dump, "--out", str(prefix), "--reverb", "spectral", "--rt60", f"{rt60:g}", "--reverb-ir-only"], + check=True) + fs, ir = wavfile.read(str(prefix) + ".reverb.wav") + return float(measure_rir.schroeder(ir[:, 0].astype(np.float64), fs)["T30"]) + + +def chain_key(meta: dict): + """The --chain-rt60 table's key for a dump's JSON line.""" + if meta["reverb"] == "spectral": + return ("spectral", float(meta["rt60"])) + return (meta["reverb"], float(meta["decay"]), float(meta["damping"]), meta["return"]) + + 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", @@ -167,7 +200,7 @@ def one_run(dump: str, work: pathlib.Path, name: str, spec: tuple, delay: int, s 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"])] + t_chain = chain_t30[chain_key(meta)] 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" @@ -217,7 +250,7 @@ def recover(dump: str, work: pathlib.Path, name: str, spec: tuple, delay: int, s 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"])] + t_chain = chain_t30[chain_key(meta)] 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) @@ -272,6 +305,8 @@ def print_table(runs: list, order: list) -> None: + f" median {fmt(statistics.median(gap))})") # The two shift topologies against each other: per seed, dry-only shift - whole-bus shift. + if not any(n.startswith("bus") for n in order): + return # the spectral rows have no shift print("\ndry-only shift - whole-bus shift, per seed (same seed, delay, shift, decay and wet), dB") pairs = 0 higher = 0 @@ -307,7 +342,8 @@ def main(argv=None) -> None: 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("--plate", default="shipped", choices=["shipped", "paper", "spectral"]) + ap.add_argument("--shape-maxes", default="0,1,2,4", help="with --plate spectral: 0 = flat, else the shape_max") 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") @@ -322,7 +358,7 @@ def main(argv=None) -> None: decays = [float(x) for x in args.decays.split(",")] wets = [float(x) for x in args.wets.split(",")] - cfgs = configs(decays, wets, args.plate) + cfgs = configs(decays, wets, args.plate, [float(x) for x in args.shape_maxes.split(",")]) order = list(cfgs) if args.merge: by_key = {} @@ -342,6 +378,11 @@ def main(argv=None) -> None: print(f"band-limited cabin (4096 taps) T30 = {room_t30:.4f} s (--rt60)") chain_t30 = {} for d in decays: + if args.plate == "spectral": + key = ("spectral", d) + chain_t30[key] = spectral_t30(args.dump, work, d) + print(f"spectral reverb, rt60 {d:g}: T30 = {chain_t30[key]:.4f} s (--chain-rt60)") + continue 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)")