diff --git a/bench/README.md b/bench/README.md index 8eccfbb..e740497 100644 --- a/bench/README.md +++ b/bench/README.md @@ -50,6 +50,50 @@ medians of 5, load average 4.11 — a shared machine, so pessimistic): 324621.6 / 324455.2 ns that `pem_afc.h` quotes for one canceller's `process_block` at the same block size on the same CPU: 0.94 % / 0.55 %. +The same file breaks one block into stages, float and double, so one run +localizes the cost on a host: `howl_bank` (the resonators and their +envelopes, one fused loop, through a detector whose tick never comes), +`howl_tick` (one `howl_detail::decide()` on envelope snapshots), +`howl_tick_logs` (its 66 `std::log10` calls alone), `howl_fit` (one band's +growth fit), `howl_readouts` (the harmonic readouts and the 32 +`band_level_db()` reads attribution makes), `bench_guard/policy_m1` (the +guard's own step, no detector) and `bench_guard/m1`, `m2` (the whole guard +per block, as `HowlGuardHost.CostPerBlock` times it). + +**The guard on Linux GCC (tap/MuTap#87).** `CostPerBlock` had read the +guard at ~4x its Intel Mac cost on the Linux GCC runner. The breakdown, on +the `Benchmark smoke` job (GCC 13.3, `4 X 3241 MHz` / `4 X 3244.61 MHz`; +medians of 5 x 0.3 s; ns per block, float / double): + +| stage | before (main) | after | +|---|---:|---:| +| `howl_bank` | 10291 / 10182 | 1288 / 2434 | +| `howl_tick` | 730 / 945 | 725 / 945 | +| `howl_tick_logs` | 541 / 759 | 517 / 767 | +| `howl_fit` | 34.9 / 36.7 | 34.5 / 36.4 | +| `howl_readouts` | 241 / 293 | 248 / 294 | +| `bench_guard/policy_m1` | 30.9 / 48.1 | 29.9 / 48.8 | +| `howl_detector/b64` | 10885 / 11135 | 2006 / 3337 | +| `bench_guard/m1` | 10919 / 11250 | 2059 / 3407 | +| `bench_guard/m2` | 18204 / 18475 | 4113 / 6754 | + +All of the excess was the bank. GCC left its band loop scalar +(`-fopt-info-vec`: "couldn't vectorize loop", "no vectype for stmt" on the +coefficient loads) and compiled the envelope's attack / release pick as a +`comiss` / `jbe` branch, which noise mispredicts; Clang vectorized the loop. +`__restrict` on the loop's local pointers changed nothing under GCC (10603 +ns, a replica); the same loop with its arrays as `__restrict` parameters +vectorized (1164 / 2273 ns), so the loop now lives in +`howl_detail::run_bank` with restrict parameters, its arithmetic unchanged. +Linux Clang on the same PR run reads 1320 / 2445 (`howl_bank`), 2000 / 3363 +(the detector) and 2061 / 3432 (`m1`): the two compilers now match. A dump of +every detector readout after every `process_block` call (4 configs, 6 +signals that drive all three trigger paths, irregular partitions; 3,440,844 +values) is bit-identical between main's header and the new one on Linux GCC, +Linux Clang and macOS arm64 CI and on the Intel Mac. The tick's +`std::log10` calls are now the largest part of what is left after the bank +(517 / 767 of the tick's 725 / 945 ns) and were left alone. + ## Scalar baselines (reference container, 2.8 GHz x86, GCC Release: -O3 -DNDEBUG, medians of 5, idle machine) | layer | 48 kHz f64 | 48 kHz f32 | 16 kHz f64 | 16 kHz f32 | diff --git a/bench/bench_howl_detector.cpp b/bench/bench_howl_detector.cpp index 03fc544..ce10537 100644 --- a/bench/bench_howl_detector.cpp +++ b/bench/bench_howl_detector.cpp @@ -2,13 +2,33 @@ // Copyright 2026 MuTap contributors // // Host benchmark for the safety layer's per-mic howl detector -// (mutap/howl_detector.h) at its defaults: 32 resonators, block 64 at -// 48 kHz, float32 and double. One iteration is one block, so the reported -// time is the per-block (per-mic) cost; items are samples, as in -// bench_aec.cpp. The input cycles through a precomputed white-noise corpus -// at unit RMS, so the bank, the envelopes and the growth ring run their -// steady path (no denormals, no trip). bench/README.md has the workflow and -// the recorded numbers. +// (mutap/howl_detector.h) and the guard over it (mutap/howl_guard.h), at the +// detector's defaults: 32 resonators, block 64 at 48 kHz, float32 and double. +// One iteration is one block, so the reported time is the per-block (per-mic) +// cost; items are samples, as in bench_aec.cpp. The input cycles through a +// precomputed white-noise corpus at unit RMS, so the bank, the envelopes and +// the growth ring run their steady path (no denormals). bench/README.md has +// the workflow and the recorded numbers. +// +// Besides the whole detector (howl_detector/b64_*) the file breaks one block +// into its stages, so one run localizes the cost on a host: +// * howl_bank the resonators and their envelopes over 64 samples (one +// fused loop in the detector: a detector whose tick never +// comes, so process_block runs the bank alone); +// * howl_tick one tick, howl_detail::decide (levels, prominence, line +// ring, per-band growth fit, level tracker, verdict) on +// envelope snapshots taken from a running detector; +// * howl_tick_logs the tick's transcendentals alone: the 32 + 32 + 2 +// std::log10 calls decide() makes on the same snapshots; +// * howl_fit one band's least-squares fit over a full ring (decide() +// runs it for the peak band and every prominent one); +// * howl_readouts the harmonic / subharmonic readouts and the 32 +// band_level_db() reads the guard's attribution makes; +// * guard_policy howl_guard::update() + apply() per mic, the detector not +// run (the guard's own step); +// * guard the whole guard per block, as HowlGuardHost.CostPerBlock +// times it: analyze() per mic, update(), apply() per mic. +#include #include #include #include @@ -17,27 +37,37 @@ #include #include "mutap/howl_detector.h" +#include "mutap/howl_guard.h" namespace { - template - void bench_howl_detector(benchmark::State& state, size_t block) { - typename tap::mu::howl_detector::config cfg; - cfg.sample_rate = 48000.0; - cfg.block_size = block; - tap::mu::howl_detector det(cfg); + constexpr double k_fs = 48000.0; + constexpr size_t k_block = 64; + constexpr size_t k_blocks = 512; ///< corpus length, blocks + constexpr size_t k_warm = 2000; - constexpr size_t k_blocks = 512; + template + std::vector white_corpus(size_t block) { std::mt19937 rng(0xB0B0); std::normal_distribution nd(0.0, 1.0); std::vector x(k_blocks * block); for (auto& v : x) { v = static_cast(nd(rng)); } - for (size_t i = 0; i < 2000; ++i) { + return x; + } + + template + void bench_howl_detector(benchmark::State& state, size_t block) { + typename tap::mu::howl_detector::config cfg; + cfg.sample_rate = k_fs; + cfg.block_size = block; + tap::mu::howl_detector det(cfg); + const std::vector x = white_corpus(block); + for (size_t i = 0; i < k_warm; ++i) { det.process_block(&x[(i % k_blocks) * block], block); } - size_t i = 2000; + size_t i = k_warm; for (auto _ : state) { det.process_block(&x[(i % k_blocks) * block], block); benchmark::DoNotOptimize(det.prominence_db()); @@ -46,7 +76,252 @@ namespace { state.SetItemsProcessed(static_cast(state.iterations()) * static_cast(block)); } + /// The resonators and envelopes alone: a detector whose tick is 2^40 + /// samples away runs only run_bank() in process_block (the growth + /// window is stretched to the 3 ticks validation needs). + template + void bench_howl_bank(benchmark::State& state) { + typename tap::mu::howl_detector::config cfg; + cfg.sample_rate = k_fs; + cfg.block_size = size_t{1} << 40; + cfg.growth_window_s = static_cast(4.0 * static_cast(cfg.block_size) / k_fs); + tap::mu::howl_detector det(cfg); + const std::vector x = white_corpus(k_block); + for (size_t i = 0; i < k_warm; ++i) { + det.process_block(&x[(i % k_blocks) * k_block], k_block); + } + size_t i = k_warm; + for (auto _ : state) { + det.process_block(&x[(i % k_blocks) * k_block], k_block); + benchmark::DoNotOptimize(det.band_level_db(0)); + ++i; + } + state.SetItemsProcessed(static_cast(state.iterations()) * static_cast(k_block)); + } + + /// k_blocks ticks of a running default detector on the corpus: each + /// tick's 32 band envelopes (from band_level_db) and block power. + template + struct snapshots { + size_t bands = 0; + std::vector env; + std::vector power; + snapshots() { + typename tap::mu::howl_detector::config cfg; + cfg.sample_rate = k_fs; + cfg.block_size = k_block; + tap::mu::howl_detector det(cfg); + const std::vector x = white_corpus(k_block); + bands = det.bands(); + for (size_t i = 0; i < k_warm; ++i) { + det.process_block(&x[(i % k_blocks) * k_block], k_block); + } + for (size_t k = 0; k < k_blocks; ++k) { + const S* in = &x[k * k_block]; + det.process_block(in, k_block); + for (size_t b = 0; b < bands; ++b) { + env.push_back(std::pow(S(10), det.band_level_db(b) / S(20))); + } + S sq = S(0); + for (size_t n = 0; n < k_block; ++n) { + sq += in[n] * in[n]; + } + power.push_back(sq / static_cast(k_block)); + } + } + }; + + /// howl_detail::decision_params as howl_detector's constructor derives + /// them at the defaults (block 64 at 48 kHz: 16 slots, growth 0.04 dB + /// per tick, ceiling -6 dB, hold 150 ticks). + template + tap::mu::howl_detail::decision_params default_params(size_t bands) { + const typename tap::mu::howl_detector::config cfg; + const double tick = static_cast(k_block) / k_fs; + tap::mu::howl_detail::decision_params p; + p.bands = bands; + p.slots = static_cast(std::lround(static_cast(cfg.growth_window_s) / tick)); + p.inv_others = static_cast(1.0 / static_cast(bands - 1)); + p.prominence = cfg.prominence_db; + p.growth_tick = + static_cast(static_cast(cfg.growth_db_per_pass) * tick / static_cast(cfg.loop_period_s)); + p.linearity = cfg.linearity_db; + p.rise = cfg.rise_db; + p.level_catch = cfg.level_catch; + p.ceiling_pow = static_cast(std::pow(10.0, static_cast(cfg.ceiling_db) / 10.0)); + p.level = cfg.level_db; + p.level_prom = cfg.level_prominence_db; + p.program_coef = static_cast(1.0 - std::exp(-tick / static_cast(cfg.program_s))); + p.hold_ticks = static_cast(std::lround(static_cast(cfg.hold_s) / tick)); + p.conf_att = static_cast(std::exp(-tick / static_cast(cfg.hold_s))); + p.conf_rel = static_cast(std::exp(-tick / (4.0 * static_cast(cfg.hold_s)))); + p.slope_w.resize(p.slots); + const double mean = 0.5 * static_cast(p.slots - 1); + double sxx = 0.0; + for (size_t k = 0; k < p.slots; ++k) { + sxx += (static_cast(k) - mean) * (static_cast(k) - mean); + } + for (size_t k = 0; k < p.slots; ++k) { + p.slope_w[k] = static_cast((static_cast(k) - mean) / sxx); + } + return p; + } + + template + void bench_howl_tick(benchmark::State& state) { + const snapshots snap; + const auto p = default_params(snap.bands); + tap::mu::howl_detail::decision_state s; + s.resize(p.bands, p.slots); + s.reset(); + for (size_t i = 0; i < k_warm; ++i) { + const size_t k = i % k_blocks; + tap::mu::howl_detail::decide(p, s, &snap.env[k * snap.bands], snap.power[k]); + } + size_t i = k_warm; + for (auto _ : state) { + const size_t k = i % k_blocks; + tap::mu::howl_detail::decide(p, s, &snap.env[k * snap.bands], snap.power[k]); + benchmark::DoNotOptimize(s.prominence); + ++i; + } + state.SetItemsProcessed(static_cast(state.iterations()) * static_cast(k_block)); + } + + /// decide()'s std::log10 calls on the same snapshots: 20 log10(env) and + /// 20 log10 of the prominence ratio per band, then the block power's two. + template + void bench_howl_tick_logs(benchmark::State& state) { + const snapshots snap; + const size_t nb = snap.bands; + std::vector level(nb); + std::vector prom(nb); + size_t i = 0; + for (auto _ : state) { + const size_t k = i % k_blocks; + const S* env = &snap.env[k * nb]; + for (size_t b = 0; b < nb; ++b) { + level[b] = S(20) * std::log10(env[b]); + } + for (size_t b = 0; b < nb; ++b) { + prom[b] = S(20) * std::log10((env[b] + S(1e-9)) / (env[(b + 1) % nb] + S(1e-9))); + } + const S pw = snap.power[k] + S(1e-30); + S a = S(10) * std::log10(pw); + S c = S(10) * std::log10(pw / S(0.9)); + benchmark::DoNotOptimize(level.data()); + benchmark::DoNotOptimize(prom.data()); + benchmark::DoNotOptimize(a); + benchmark::DoNotOptimize(c); + benchmark::ClobberMemory(); + ++i; + } + state.SetItemsProcessed(static_cast(state.iterations()) * static_cast(k_block)); + } + + template + void bench_howl_fit(benchmark::State& state) { + const snapshots snap; + const auto p = default_params(snap.bands); + tap::mu::howl_detail::decision_state s; + s.resize(p.bands, p.slots); + s.reset(); + for (size_t i = 0; i < k_warm; ++i) { + const size_t k = i % k_blocks; + tap::mu::howl_detail::decide(p, s, &snap.env[k * snap.bands], snap.power[k]); + } + size_t b = 0; + for (auto _ : state) { + S slope; + S residual; + tap::mu::howl_detail::fit(p, s, b, slope, residual); + benchmark::DoNotOptimize(slope); + benchmark::DoNotOptimize(residual); + b = (b + 1 == p.bands) ? 0 : b + 1; + } + } + + template + void bench_howl_readouts(benchmark::State& state) { + typename tap::mu::howl_detector::config cfg; + cfg.sample_rate = k_fs; + cfg.block_size = k_block; + tap::mu::howl_detector det(cfg); + const std::vector x = white_corpus(k_block); + for (size_t i = 0; i < k_warm; ++i) { + det.process_block(&x[(i % k_blocks) * k_block], k_block); + } + for (auto _ : state) { + benchmark::DoNotOptimize(det.harmonic_ratio_db()); + benchmark::DoNotOptimize(det.subharmonic_ratio_db()); + for (size_t b = 0; b < det.bands(); ++b) { + benchmark::DoNotOptimize(det.band_level_db(b)); + } + } + } + + /// The guard as tests/support/guard_loop.h builds it (S1's 10 ms loop + /// period, the harness's +30 dB ceiling: no trip on the corpus). + template + typename tap::mu::howl_guard::config guard_config(size_t mics) { + typename tap::mu::howl_guard::config c; + c.microphones = mics; + c.block_size = k_block; + c.sample_rate = k_fs; + c.detector.loop_period_s = S(0.010); + c.detector.ceiling_db = S(30); + return c; + } + + template + void bench_guard(benchmark::State& state, size_t mics, bool analyze) { + tap::mu::howl_guard g(guard_config(mics)); + const std::vector x = white_corpus(k_block); + std::vector out(k_block); + const auto step = [&](size_t i) { + const S* e = &x[(i % k_blocks) * k_block]; + if (analyze) { + for (size_t m = 0; m < mics; ++m) { + g.analyze(m, e, S(1e-3), S(0.5)); + } + } + g.update(); + for (size_t m = 0; m < mics; ++m) { + g.apply(m, e, out.data()); + } + }; + for (size_t m = 0; m < mics; ++m) { + g.analyze(m, x.data(), S(1e-3), S(0.5)); + } + for (size_t i = 0; i < k_warm; ++i) { + step(i); + } + size_t i = k_warm; + for (auto _ : state) { + step(i); + benchmark::DoNotOptimize(out.data()); + ++i; + } + state.SetItemsProcessed(static_cast(state.iterations()) * static_cast(k_block)); + } + } // namespace BENCHMARK_CAPTURE(bench_howl_detector, b64_f64, 64); BENCHMARK_CAPTURE(bench_howl_detector, b64_f32, 64); +BENCHMARK_TEMPLATE(bench_howl_bank, double); +BENCHMARK_TEMPLATE(bench_howl_bank, float); +BENCHMARK_TEMPLATE(bench_howl_tick, double); +BENCHMARK_TEMPLATE(bench_howl_tick, float); +BENCHMARK_TEMPLATE(bench_howl_tick_logs, double); +BENCHMARK_TEMPLATE(bench_howl_tick_logs, float); +BENCHMARK_TEMPLATE(bench_howl_fit, double); +BENCHMARK_TEMPLATE(bench_howl_fit, float); +BENCHMARK_TEMPLATE(bench_howl_readouts, double); +BENCHMARK_TEMPLATE(bench_howl_readouts, float); +BENCHMARK_CAPTURE(bench_guard, policy_m1_f64, 1, false); +BENCHMARK_CAPTURE(bench_guard, policy_m1_f32, 1, false); +BENCHMARK_CAPTURE(bench_guard, m1_f64, 1, true); +BENCHMARK_CAPTURE(bench_guard, m1_f32, 1, true); +BENCHMARK_CAPTURE(bench_guard, m2_f64, 2, true); +BENCHMARK_CAPTURE(bench_guard, m2_f32, 2, true); diff --git a/include/mutap/howl_detector.h b/include/mutap/howl_detector.h index 135d0a7..6a36ee1 100644 --- a/include/mutap/howl_detector.h +++ b/include/mutap/howl_detector.h @@ -140,7 +140,10 @@ // COST (bench/bench_howl_detector.cpp, defaults, block 64 at 48 kHz, the same // i9-8950HK at load 4.1): 1786 ns per block in float32, 3039 ns in double - // 0.55 % / 0.94 % of one pem_afc process_block (pem_afc.h: 324455.2 / -// 324621.6 ns). +// 324621.6 ns). On the Linux GCC CI runner (GCC 13.3): 2006 / 3337 ns, the +// same as Linux Clang there, since run_bank's arrays became __restrict +// parameters (10885 / 11135 ns before: GCC did not vectorize the bank; +// tap/MuTap#87, bench/README.md). // // Real-time contract (as frequency_shifter.h): the constructor validates its // config, may throw std::invalid_argument and allocates (the bank's state @@ -257,6 +260,42 @@ namespace tap::mu { } }; + /// The bank over len samples: nb resonators (w = x + bias - a1 w1 - + /// a2 w2, y = |g (w - w2)| + env_floor) and their attack / release + /// envelopes. Returns sq plus the sum of the inputs' squares. + /// + /// A free function for one reason: its arrays are __restrict + /// PARAMETERS (each is the detector's own, so none aliases another + /// or the input). The vectorizer runs the band loop across bands; + /// without the qualifiers GCC 13 left it scalar, with a branch on + /// y > env that noise mispredicts, while Clang vectorized it - 10103 + /// vs 1164 ns per block on the Linux GCC runner (float; double 10383 + /// vs 2273), the guard's whole Linux excess. The same qualifiers on + /// local pointers did not help GCC; on parameters it honours them. + /// Each band's arithmetic is unchanged, operation for operation + /// (every readout measured bit-identical to the previous loop's on + /// Linux GCC / Clang and macOS; bench/README.md, tap/MuTap#87). + template + Sample run_bank(const Sample* __restrict in, size_t len, size_t nb, const Sample* __restrict a1, + const Sample* __restrict a2, const Sample* __restrict g, Sample* __restrict w1, + Sample* __restrict w2, Sample* __restrict env, Sample att, Sample rel, Sample sq, Sample bias, + Sample env_floor) noexcept { + for (size_t i = 0; i < len; ++i) { + const Sample x = in[i]; + const Sample xb = x + bias; + sq += x * x; + for (size_t b = 0; b < nb; ++b) { + const Sample w = xb - a1[b] * w1[b] - a2[b] * w2[b]; + const Sample y = std::abs(g[b] * (w - w2[b])) + env_floor; + w2[b] = w1[b]; + w1[b] = w; + const Sample c = (y > env[b]) ? att : rel; + env[b] = y + c * (env[b] - y); + } + } + return sq; + } + /// Least-squares slope (dB per tick) and RMS fit residual (dB) of /// band b's line level over the ring, oldest slot first. template @@ -678,30 +717,9 @@ namespace tap::mu { /// The resonators and their envelopes over len samples. void run_bank(const Sample* in, size_t len) noexcept { - const size_t nb = m_cfg.bands; - const Sample* a1 = m_a1.data(); - const Sample* a2 = m_a2.data(); - const Sample* g = m_g.data(); - Sample* w1 = m_w1.data(); - Sample* w2 = m_w2.data(); - Sample* env = m_env.data(); - const Sample att = m_att; - const Sample rel = m_rel; - Sample sq = m_block_sq; - for (size_t i = 0; i < len; ++i) { - const Sample x = in[i]; - const Sample xb = x + k_bias; - sq += x * x; - for (size_t b = 0; b < nb; ++b) { - const Sample w = xb - a1[b] * w1[b] - a2[b] * w2[b]; - const Sample y = std::abs(g[b] * (w - w2[b])) + k_env_floor; - w2[b] = w1[b]; - w1[b] = w; - const Sample c = (y > env[b]) ? att : rel; - env[b] = y + c * (env[b] - y); - } - } - m_block_sq = sq; + m_block_sq = + howl_detail::run_bank(in, len, m_cfg.bands, m_a1.data(), m_a2.data(), m_g.data(), m_w1.data(), + m_w2.data(), m_env.data(), m_att, m_rel, m_block_sq, k_bias, k_env_floor); } config m_cfg; diff --git a/tests/test_howl_guard_host.cpp b/tests/test_howl_guard_host.cpp index da91c54..0307324 100644 --- a/tests/test_howl_guard_host.cpp +++ b/tests/test_howl_guard_host.cpp @@ -1651,7 +1651,11 @@ TEST(HowlGuardHost, CostPerBlock) { // 2246 / 3962 and 3180 / 6311 ns). // Linux GCC CI (job 110465810942): float 8575 / 12794 ns, 4.61 / // 6.86 %; double 8706 / 13571 ns, 4.59 / 7.15 % (the canceller - // there 186-190 us, the guard 4x Intel's). + // there 186-190 us, the guard 4x Intel's): GCC did not vectorize + // the detector's bank. Since tap/MuTap#87 (run_bank's arrays are + // __restrict parameters; bench/README.md), job 113795845819: + // float 1998 / 3986 ns, 0.62 / 1.24 %; double 3371 / 6793 ns, + // 1.03 / 2.07 % (that runner's canceller 321-327 us). // A wall-clock ratio on shared CI hardware moves by the factor // between those hosts, so it is not gated tightly: the bound only // catches a gross regression (25 %, 3.5x over the worst CI row).