diff --git a/HANDOFF.md b/HANDOFF.md index 88e1713..f51a889 100644 --- a/HANDOFF.md +++ b/HANDOFF.md @@ -293,6 +293,24 @@ residual-echo estimate, spectral suppression, comfort noise matched to the near-end noise floor; plus the integrated chain type. House workflow: scratch-measure first, thresholds with margin, rooms from both generator families. RT contract as everywhere. +*DONE — `mutap/postfilter.h` (`residual_suppressor` + `aec_chain`) + +`tests/test_postfilter.cpp` (12 tests, measured-first; every Stage 2 +deliverable inside its margin target — the measured table and the four +design decisions the numbers forced live in the matrix's "Stage 2 +delivered" section). The headline discovery: the AEC chain must NOT +use PEM — open-loop AEC has an exogenous far end and the predictor +refit floors misalignment near -20 dB where the raw FD-Kalman core +(transition 0.9998, initial uncertainty 10 — the measured AEC sweet +spot) reaches -75 dBm0(A) bare with double-talk immunity from its own +noise-PSD tracker. The suppressor correlates the MIC, not E, against +the echo estimate (orthogonality principle: coh(E,Yhat) saturates at +0.34 while adapting), and takes suppression DEPTH from a +coherence-gated leakage estimate so double-talk transparency is the +shape of the rule, not a detector. Two instrument fixes en route: +A-weighting poles above Nyquist now prewarp-clamped (16 kHz tests went +NaN), and pem_afc grew echo_estimate_block(). mutap.h includes the new +header; emulated-target selections unchanged (double-only ITU suites +stay host-side).* **Stage 3 — Compliance suite** (`tests/test_itu_*.cpp`). One gtest per matrix row asserting requirement + margin policy. Swept across fixture diff --git a/README.md b/README.md index f5b0ac1..45e3165 100644 --- a/README.md +++ b/README.md @@ -187,6 +187,23 @@ algorithmically complete. What exists today: double-talk the warped predictor beat the speech cascade's suppression in **all 18 room × seed pairs** measured (per-room medians 17.3…19.6 vs 15.2…16.7 dB, Kalman core). +- **Residual-echo post-filter + comfort noise** + ([`mutap/postfilter.h`](include/mutap/postfilter.h)) — the Stage 2 + deliverable of the [ITU compliance plan](docs/itu-compliance.md): + a per-bin Wiener suppressor driven by mic-vs-echo-estimate coherence + gating a learned leakage estimate, comfort noise matched to the + near-end floor by two-window minimum statistics, and `mutap::aec_chain` + composing it with a linear canceller (default: the **raw** FD-Kalman + core — open-loop AEC has an exogenous far end, so PEM's decorrelation + buys nothing and its predictor refit floors misalignment near −20 dB + where the raw core reaches −75 dBm0(A) bare). Measured on the ITU + battery ([`tests/test_postfilter.cpp`](tests/test_postfilter.cpp)): + single-talk residual **−79.9/−88.9 dBm0(A)** (cabin/studio; the P.1120 + clause wants < −58, our margin target < −64), double-talk near-end + attenuation **1.05 dB** (clause ≤ 3, target ≤ 1.5), double-talk echo + loss **≥ 34.1 dB in every band** (clause ≥ 27, target ≥ 33), comfort + noise matched **−1.2 dB / ≤ 1.7 dB per band** (clause +2/−5, half-mask + spectrum), noise pumping 3.3 dB, near-end build-up 20.9 ms. Next up (see [HANDOFF.md](HANDOFF.md) "What's next"): in-Max listening in a real room and the default-engine decision, then the M55 performance diff --git a/docs/itu-compliance.md b/docs/itu-compliance.md index e16573c..1644bc8 100644 --- a/docs/itu-compliance.md +++ b/docs/itu-compliance.md @@ -333,3 +333,53 @@ analogues), Hoth-spectrum noise for the P.340 rows. - Switching dynamics inside 25 ms build-up at the margin targets. - No noise pumping beyond 5 dB around speech bursts. - The chain's ERL trajectory above the convergence masks at 2x speed. + +### Stage 2 delivered (mutap/postfilter.h + tests/test_postfilter.cpp) + +All of the above, measured against this file's margin targets on the +golden model (double, block 256, 2048-tap unit-energy paths, 48 kHz; +one required-rate check at 16 kHz). Chain = RAW partitioned FD-Kalman +(transition 0.9998, initial uncertainty 10) + residual suppressor at +defaults: + +| Deliverable | Margin target | Measured | +|---|---|---| +| ST residual, cabin / studio | < -64 dBm0(A) | **-79.9 / -88.9** | +| ST residual, NLMS core / 16 kHz | < -64 dBm0(A) | -65.0 / -69.1 | +| DT send attenuation, cabin / studio | <= 1.5 dB | **1.05 / 0.93** | +| DT echo loss, worst band 200-6950 Hz | >= 33 dB | **38.0 / 34.1** | +| ERL by 600 / 1200 ms | >= 40 / >= 46 dB | 43.2 / 46.9 | +| Comfort-noise level match | +1 / -2.5 dB | -1.20 | +| Comfort-noise spectrum, worst band | half-mask (+-3..6) | 1.69 dB | +| Noise pumping | <= 5 dB | 3.3 | +| Near-end build-up at DT onset | <= 25 ms | 20.9 ms | + +Design decisions the numbers forced (full derivations in +postfilter.h's comments, rejected designs kept in git history): + +1. **The AEC chain does not use PEM.** Open-loop AEC has an exogenous + far end — no closed-loop bias to remove — and the predictor's + block-by-block refit injects gradient noise that floors misalignment + near -20 dB. The raw Kalman core measures -75.6 dBm0(A) bare on the + same scenario where PEM-Kalman plateaus at -28.9, and its per-bin + noise-PSD tracking IS the double-talk defense (post-DT residual + -70.9). PEM remains the right structure for the closed loop (AFC), + and aec_chain still composes with pem_afc. +2. **The suppressor's discriminator correlates the MIC (E + Yhat), + never E, against Yhat**: a converged adaptive filter keeps E + orthogonal to its reference, so coh(E, Yhat) saturates at 0.34 while + adapting; coh(D, Yhat) measures 0.99. +3. **Suppression depth comes from a leakage estimate gated by that + coherence, not from the coherence itself** — the AM-FM plans + interleave 20 Hz apart at the low end, no realizable gain filter + notches that selectively, and a coherence-proportional gain measured + 5.5 dB of near-end attenuation there. The Wiener-on-E form rides to + unity wherever near-end energy inflates |E|^2. +4. **Comfort noise fills to a two-window minimum-statistics floor + (bias x4)** — asymmetric-rate one-pole trackers measured either + -8 dB undershoot (raw-minimum bias) or tens-of-seconds acquisition. + +Still Stage 3's to prove: the full per-row multi-rate suite (both +required rates x both cores x all rooms), the spectral echo mask, the +switching/activation battery, stability sweep, G.168-adapted rows, and +the TCL / time-variant-path rows. diff --git a/include/mutap/mutap.h b/include/mutap/mutap.h index d161c80..7095b6a 100644 --- a/include/mutap/mutap.h +++ b/include/mutap/mutap.h @@ -21,3 +21,4 @@ #include "mutap/fft.h" #include "mutap/lpc.h" #include "mutap/pem_afc.h" +#include "mutap/postfilter.h" diff --git a/include/mutap/pem_afc.h b/include/mutap/pem_afc.h index e7c238e..cf517e2 100644 --- a/include/mutap/pem_afc.h +++ b/include/mutap/pem_afc.h @@ -178,6 +178,12 @@ namespace mutap { /// Current feedback-path estimate as filter_length() time-domain taps. void copy_impulse_response(Sample* dest) noexcept { m_fdaf.copy_impulse_response(dest); } + /// The echo estimate of the LAST processed block — block_size() + /// time-domain samples of what the cancellation subtracted (the + /// residual suppressor's reference, postfilter.h). Valid until + /// the next process_block()/reset(). + const Sample* echo_estimate_block() const noexcept { return m_time.data() + block_size(); } + private: static config validated(const config& cfg) { if (cfg.analysis_window < 2 * cfg.fdaf.block_size) { diff --git a/include/mutap/postfilter.h b/include/mutap/postfilter.h new file mode 100644 index 0000000..6f6bc18 --- /dev/null +++ b/include/mutap/postfilter.h @@ -0,0 +1,594 @@ +/// @file postfilter.h +/// @brief Residual-echo suppressor + comfort noise (the AEC post-filter), +/// and aec_chain — the linear canceller and the post-filter composed +/// into the unit the ITU compliance matrix measures. +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "mutap/fd_kalman.h" +#include "mutap/fft.h" +#include "mutap/pem_afc.h" + +namespace mutap { + + /// Residual-echo suppressor with matched comfort noise + /// (docs/itu-compliance.md, Stage 2). + /// + /// Linear cancellation measures ~20-25 dB of echo suppression; the + /// single-talk clauses of the automotive recs want a residual below + /// -58 dBm0(A) — every compliant product closes that gap with a + /// post-filter, and the double-talk clauses bound what it may cost: + /// no more than 3 dB of near-end attenuation (our margin target: + /// 1.5 dB). Both sides fall out of one rule: + /// + /// coh(k) = |S_dy(k)|^2 / (S_dd(k) S_yy(k)) MIC-vs-Yhat coherence + /// lam(k) <- S_ee(k) / S_yy(k) learned ONLY while coh(k) > gate + /// G(k) = max(g_min, 1 - beta lam(k) S_yy(k) / S_ee(k)) + /// + /// where Yhat is the linear canceller's own echo-estimate spectrum + /// and D = E + Yhat reconstructs the microphone. Two estimators + /// split the two jobs a residual suppressor has: + /// + /// - coh(D, Yhat) answers "is this bin echo right now?". The + /// discriminator deliberately correlates the MIC — not the + /// canceller output E — against Yhat: a converged adaptive + /// filter keeps E orthogonal to its reference by construction, + /// and its weight jitter makes the tiny residual's phase wander, + /// so coh(E, Yhat) saturates well below 1 (measured 0.34 + /// adapting / 0.53 frozen — the freeze experiment that located + /// this). The mic instead carries the FULL echo H x, so per bin + /// D/Yhat ~ H/Hhat ~ 1: insensitive to adaptation jitter and to + /// the convolution-tail (frame-multiplicative) error, since both + /// signals traverse nearly the same filter. + /// + /// - lam(k), the LEAKAGE (residual power per unit estimated-echo + /// power ~ the canceller's per-bin misalignment), answers "how + /// much of E does echo explain?". It is learned only while the + /// coherence certifies the bin echo-dominant, and held through + /// double talk. The applied gain is the Wiener gain ON E: during + /// single talk S_ee ~ lam S_yy so the gain floors; during double + /// talk near-end energy inflates S_ee while lam S_yy holds, so + /// the gain rides to 1 — even in bins where the analysis + /// resolution cannot separate the near-end line from an echo + /// line 20 Hz away (P.501's AM-FM plans interleave that tightly + /// at the low end; a gain keyed to coherence alone measurably + /// took 5.5 dB off the near end there, hitting both signals + /// 1-for-1). Transparency is intrinsic: not a detector decision + /// but the shape of the rule. + /// + /// Suppressed bins are refilled with COMFORT NOISE shaped to the + /// tracked near-end noise floor N(k) (two-window minimum + /// statistics): out = G E + sqrt(max(0, N - G^2 |E|^2)) . unit-noise, + /// so the output sits at the true noise floor instead of pumping + /// between "noise" and "digital silence" (P.111x comfort-noise + /// clauses: level within +2/-5 dB, spectrum within mask). + /// + /// Framing matches the canceller: overlap-save on [previous block, + /// current block], per-bin gains, last half kept — zero added + /// latency beyond the canceller's own block. Gains are time-smoothed + /// (fast attack toward suppression, slower release) and floor-capped + /// to bound musical noise. + /// + /// Real-time contract: constructor allocates and may throw; every + /// post-construction entry point is noexcept and allocation-free + /// (the comfort-noise source is a xorshift PRNG, not ). + template + class residual_suppressor { + public: + struct config { + size_t block_size = 256; ///< must match the canceller's block size + /// Analysis window = this many blocks (power of 2). The + /// suppressor's frequency resolution is what separates echo + /// from near-end structure — at block 256 / 48 kHz one block + /// is 93.75 Hz per bin, too coarse for the AM-FM double-talk + /// combs (~90-160 Hz interleave); 8 blocks (23.4 Hz bins, + /// Hann-windowed estimation) resolves them. + size_t analysis_blocks = 8; + Sample max_suppression_db = 40; ///< g_min = -this in dB (floor per bin) + Sample over_subtraction = Sample(1.2); ///< beta on the coherence (>= 1) + /// One-pole smoothing of the coherence/power accumulators, + /// [0, 1) (~100 ms at block 256 / 48 kHz). + Sample leakage_smoothing = Sample(0.95); + /// Coherence above which a bin is certified echo-dominant + /// and its leakage may be re-learned, [0, 1). + Sample leakage_gate = Sample(0.9); + /// Per-bin gain smoothing toward MORE suppression (attack) + /// and toward LESS (release), [0, 1). Attack is the SLOW one + /// (~265 ms at block 256 / 48 kHz): estimator variance at + /// bins the analysis cannot fully resolve makes an + /// instantaneous attack chatter, and the duty-cycle loss + /// measured 3 dB off the 250 Hz near-end comb line during + /// double talk (1.6 dB integrated; 1.0 dB at 0.98). Echo + /// containment doesn't need the postfilter to be fast — the + /// canceller's Kalman gain is the fast defense. Release is + /// the near-end onset path — the switching clauses' build-up + /// budget (T_r <= 25 ms target): 0.9 measured 26.0 ms to + /// within 3 dB of the settled double-talk send level, 0.8 + /// measures 15.7 ms and slower release saturates there. + Sample gain_attack = Sample(0.98); + Sample gain_release = Sample(0.85); + /// Noise-floor tracker (two-window minimum statistics on + /// smoothed |E|^2): power smoothing [0, 1), window length in + /// blocks (floor reacts within 1-2 windows; speech bursts + /// shorter than a window cannot drag it up), and the bias + /// factor compensating the minimum's undershoot of the mean + /// (calibrated against the P.111x comfort-noise level + /// clause; see the measured table in test_postfilter.cpp). + Sample floor_smoothing = Sample(0.9); + size_t floor_window = 128; + Sample floor_bias = Sample(4); + bool comfort_noise = true; ///< fill suppressed bins to the noise floor + }; + + explicit residual_suppressor(const config& cfg) + : m_cfg(validated(cfg)) + , m_n(cfg.analysis_blocks * cfg.block_size) + , m_fft(m_n) + , m_g_min(std::pow(Sample(10), -cfg.max_suppression_db / Sample(20))) + , m_input(m_n) + , m_ywin(m_n) + , m_yspec(m_n) + , m_dspec(m_n) + , m_window(m_n) + , m_spec(m_n) + , m_time(m_n) + , m_sdy_re(m_n / 2 + 1) + , m_sdy_im(m_n / 2 + 1) + , m_sdd(m_n / 2 + 1) + , m_syy(m_n / 2 + 1) + , m_see(m_n / 2 + 1) + , m_leak(m_n / 2 + 1) + , m_gain(m_n / 2 + 1) + , m_psm(m_n / 2 + 1) + , m_min_cur(m_n / 2 + 1) + , m_min_prev(m_n / 2 + 1) + , m_gspec(m_n) + , m_gtime(m_n) { + reset(); + } + + size_t block_size() const noexcept { return m_cfg.block_size; } + size_t analysis_bins() const noexcept { return m_n / 2 + 1; } + /// Introspection for tests/diagnostics: the smoothed per-bin gain + /// and tracked leakage (analysis-bin index). + Sample gain_at(size_t k) const noexcept { return m_gain[k]; } + /// The mic-vs-echo-estimate coherence gating the leakage learner. + Sample coherence_at(size_t k) const noexcept { + const Sample c = + (m_sdy_re[k] * m_sdy_re[k] + m_sdy_im[k] * m_sdy_im[k]) / (m_sdd[k] * m_syy[k] + Sample(1e-20)); + return c > Sample(1) ? Sample(1) : c; + } + /// The learned per-bin leakage lam (the canceller's effective + /// per-bin misalignment; 1 until first certified echo-dominant). + Sample leakage_at(size_t k) const noexcept { return m_leak[k]; } + + void reset() noexcept { + for (size_t i = 0; i < m_n; ++i) { // Hann for the ESTIMATION ffts + m_window[i] = Sample(0.5) + - Sample(0.5) + * std::cos(Sample(2) * static_cast(std::numbers::pi) + * static_cast(i) / static_cast(m_n - 1)); + } + for (auto& v : m_input) { + v = Sample(0); + } + for (auto& v : m_ywin) { + v = Sample(0); + } + for (auto& v : m_sdy_re) { + v = Sample(0); + } + for (auto& v : m_sdy_im) { + v = Sample(0); + } + for (auto& v : m_sdd) { + v = Sample(0); + } + for (auto& v : m_see) { + v = Sample(0); + } + for (auto& v : m_leak) { + v = Sample(1); // pessimistic: unconverged canceller leaks everything + } + for (auto& v : m_syy) { + v = Sample(0); + } + for (auto& v : m_gain) { + v = Sample(1); + } + for (auto& v : m_psm) { + v = Sample(0); + } + // The PREVIOUS-window minimum starts at 0: the floor — and + // with it the comfort fill — stays OFF until one full window + // of real signal has been observed. + for (auto& v : m_min_cur) { + v = std::numeric_limits::infinity(); + } + for (auto& v : m_min_prev) { + v = Sample(0); + } + m_min_count = 0; + m_rng = 0x2545F491U; + } + + /// Process one block: e is the linear canceller's output, + /// yhat_block the echo estimate for the SAME block + /// (pem_afc::echo_estimate_block(), block_size time samples — + /// windowed HERE at analysis resolution: the canceller's own + /// 2*block bins cannot resolve the AM-FM double-talk combs, and + /// a coarse reference bleeds echo credit onto near-end bins), + /// out receives the suppressed block. e and out may alias. The + /// constrained gain filter is linear-phase with block_size + /// samples of delay — the postfilter's only added latency. + void process_block(const Sample* e, const Sample* yhat_block, Sample* out) noexcept { + const size_t b = m_cfg.block_size; + const size_t bins = m_n / 2 + 1; + const Sample eps = Sample(1e-20); + + // Slide both analysis windows by one block. + for (size_t i = 0; i + b < m_n; ++i) { + m_input[i] = m_input[i + b]; + m_ywin[i] = m_ywin[i + b]; + } + for (size_t i = 0; i < b; ++i) { + m_input[m_n - b + i] = e[i]; + m_ywin[m_n - b + i] = yhat_block[i]; + } + // Signal-path FFT stays RECTANGULAR (exact overlap-save + // filtering); the ESTIMATION FFTs are Hann-windowed — + // rectangular sidelobes (-13 dB) bleed echo energy into + // near-end bins and were a measured double-talk + // transparency failure. + for (size_t i = 0; i < m_n; ++i) { + m_spec[i] = m_input[i]; + // d = e + yhat reconstructs the MICROPHONE frame — the + // discriminator's left-hand signal (see class comment). + m_dspec[i] = m_window[i] * (m_input[i] + m_ywin[i]); + m_yspec[i] = m_window[i] * m_ywin[i]; + } + m_fft.forward_inplace(m_spec.data()); + m_fft.forward_inplace(m_dspec.data()); + m_fft.forward_inplace(m_yspec.data()); + + // Pass 1 per analysis bin: gated power-domain leakage -> + // Wiener gain -> time smoothing (attack 0 = same-block + // suppression, so talk onsets cannot burst through). + for (size_t k = 0; k < bins; ++k) { + Sample d_re; + Sample d_im; + Sample y_re; + Sample y_im; + packed_bin(m_dspec.data(), k, m_n, d_re, d_im); + packed_bin(m_yspec.data(), k, m_n, y_re, y_im); + + const Sample pd = d_re * d_re + d_im * d_im + eps; + const Sample py = y_re * y_re + y_im * y_im + eps; + // E = D - Yhat by linearity of the windowed transform. + const Sample e_re = d_re - y_re; + const Sample e_im = d_im - y_im; + const Sample pe = e_re * e_re + e_im * e_im + eps; + + // THE DISCRIMINATOR: magnitude-squared coherence between + // the MIC (d = e + yhat) and the echo estimate, from + // smoothed COMPLEX cross-spectra. In single talk + // D/Yhat ~ H/Hhat ~ 1 per bin — stable against the + // canceller's weight jitter — so coh -> 1 and the gain + // collapses. Near-end energy is incoherent with Yhat + // and DILUTES coh: g = 1 - coh is the Wiener near-end- + // preserving gain, so double-talk transparency is + // intrinsic, not a detector decision. (Three rejected + // designs are in git history: power-leakage ratios hug + // their noise floor or re-learn near-end energy during + // sustained double talk; rules keyed to INSTANTANEOUS + // |E|^2 excursions read loud residual as near-end and + // pass exactly the blocks that matter; and coh(E, Yhat) + // saturates at 0.34 while adapting — the update keeps E + // orthogonal to its reference, so the canceller output + // is the one signal the echo estimate cannot explain.) + const Sample a_r = m_cfg.leakage_smoothing; + m_sdy_re[k] += (Sample(1) - a_r) * ((d_re * y_re + d_im * y_im) - m_sdy_re[k]); + m_sdy_im[k] += (Sample(1) - a_r) * ((d_im * y_re - d_re * y_im) - m_sdy_im[k]); + m_sdd[k] += (Sample(1) - a_r) * (pd - m_sdd[k]); + m_syy[k] += (Sample(1) - a_r) * (py - m_syy[k]); + m_see[k] += (Sample(1) - a_r) * (pe - m_see[k]); + + Sample coh = (m_sdy_re[k] * m_sdy_re[k] + m_sdy_im[k] * m_sdy_im[k]) / (m_sdd[k] * m_syy[k] + eps); + coh = coh > Sample(1) ? Sample(1) : coh; + if (m_syy[k] <= eps * Sample(100)) { + coh = Sample(0); // no estimated echo: nothing to suppress + } + + // Re-learn the leakage only while the coherence certifies + // the bin echo-dominant; hold it through double talk. + if (coh > m_cfg.leakage_gate) { + const Sample lam = m_see[k] / (m_syy[k] + eps); + m_leak[k] += (Sample(1) - a_r) * (lam - m_leak[k]); + } + + // Wiener gain ON E, with the estimated residual PSD + // lam S_yy in the numerator. The denominator takes the + // INSTANTANEOUS |E|^2 when it exceeds the smoothed power + // — a near-end onset inside the smoother's lag must ride + // through, not get clipped by yesterday's S_ee. + const Sample se = m_see[k] > pe ? m_see[k] : pe; + Sample r = m_leak[k] * m_syy[k] / (se + eps); + r = r > Sample(1) ? Sample(1) : r; + Sample g = Sample(1) - m_cfg.over_subtraction * r; + g = g < m_g_min ? m_g_min : (g > Sample(1) ? Sample(1) : g); + const Sample a_g = g < m_gain[k] ? m_cfg.gain_attack : m_cfg.gain_release; + m_gain[k] += (Sample(1) - a_g) * (g - m_gain[k]); + } + + // Pass 2: CONSTRAIN the gain response to a causal block_size- + // tap linear-phase filter (the postfilter's version of the + // FDAF's gradient constraint). Unconstrained per-bin gains + // act as a CIRCULAR filter on the frame and smear suppression + // into bins that must stay transparent — measured: 4-8 dB of + // near-end attenuation during double talk from this alone. + // Constrained and causalized, the last block of the frame is + // a true linear convolution; the price is the filter's + // block_size/2-sample linear-phase delay. + m_gspec[0] = m_gain[0]; + m_gspec[1] = m_gain[m_n / 2]; + for (size_t k = 1; k < m_n / 2; ++k) { + m_gspec[2 * k] = m_gain[k]; + m_gspec[2 * k + 1] = Sample(0); + } + m_fft.inverse(m_gspec.data(), m_gtime.data()); + // Causal window: taps [0, 2b) of h_c[tau] = g_ir[(tau - b) mod n] + // — a 2*block-tap filter resolves gain structure down to + // fs/(2*block) (needed for the double-talk combs; a + // block-tap filter measurably smeared suppression onto the + // near-end lines). Its linear phase is the postfilter's + // latency: one block. + for (size_t tau = 0; tau < m_n; ++tau) { + m_gspec[tau] = Sample(0); + } + for (size_t tau = 0; tau < 2 * b; ++tau) { + // Rectangular cut, deliberately: a Hann taper here was + // measured WORSE for double-talk transparency (2.25 dB + // near-end attenuation vs 1.60) — the taper's doubled + // mainlobe smears deep echo notches onto near-end comb + // lines more than the rectangular sidelobes do. + m_gspec[tau] = m_gtime[(tau + m_n - b) % m_n]; + } + m_fft.forward_inplace(m_gspec.data()); + + // Pass 3: apply the constrained gains + comfort-noise fill. + for (size_t k = 0; k < bins; ++k) { + Sample e_re; + Sample e_im; + Sample g_re; + Sample g_im; + packed_bin(m_spec.data(), k, m_n, e_re, e_im); + packed_bin(m_gspec.data(), k, m_n, g_re, g_im); + + Sample out_re = g_re * e_re - g_im * e_im; + Sample out_im = g_re * e_im + g_im * e_re; + + // Near-end noise floor by TWO-WINDOW MINIMUM STATISTICS + // on the smoothed PRE-gain power |E|^2: the minimum of + // the last one-to-two windows of smoothed power is the + // stationary floor — speech and echo bursts shorter than + // a window cannot drag it up, and the floor is acquired + // within one window of the first pause. (Two rejected + // trackers, both measured: asymmetric fast-fall/slow-rise + // one-poles either chase every downward fluctuation of + // the exponentially-distributed bin powers — the raw- + // minimum bias, -8 dB comfort-noise undershoot — or, + // rise-limited, take tens of seconds to acquire the + // floor at all.) The bias factor compensates what is + // left of the minimum's undershoot of the mean. Fill + // the gap between floor and suppressed power with + // random-phase noise so suppression never gates below + // the true near-end floor. + const Sample p_in = e_re * e_re + e_im * e_im; + m_psm[k] += (Sample(1) - m_cfg.floor_smoothing) * (p_in - m_psm[k]); + if (m_psm[k] < m_min_cur[k]) { + m_min_cur[k] = m_psm[k]; + } + if (m_cfg.comfort_noise) { + const Sample fl = m_cfg.floor_bias * (m_min_cur[k] < m_min_prev[k] ? m_min_cur[k] : m_min_prev[k]); + const Sample po = out_re * out_re + out_im * out_im; + const Sample deficit = fl - po; + if (deficit > Sample(0)) { + const Sample amp = std::sqrt(deficit / Sample(2)); + out_re += amp * next_noise(); + out_im += amp * next_noise(); + } + } + store_bin(m_spec.data(), k, m_n, out_re, out_im); + } + + // Rotate the minimum-statistics windows. + if (++m_min_count >= m_cfg.floor_window) { + m_min_count = 0; + for (size_t k = 0; k < bins; ++k) { + m_min_prev[k] = m_min_cur[k]; + m_min_cur[k] = std::numeric_limits::infinity(); + } + } + + m_fft.inverse(m_spec.data(), m_time.data()); + for (size_t i = 0; i < b; ++i) { + out[i] = m_time[m_n - b + i]; + } + } + + private: + static config validated(const config& cfg) { + if (cfg.block_size < 4 || (cfg.block_size & (cfg.block_size - 1)) != 0) { + throw std::invalid_argument("residual_suppressor: block_size must be a power of 2 >= 4"); + } + if (cfg.analysis_blocks < 4 || (cfg.analysis_blocks & (cfg.analysis_blocks - 1)) != 0) { + throw std::invalid_argument("residual_suppressor: analysis_blocks must be a power of 2 >= 4"); + } + if (cfg.max_suppression_db <= Sample(0)) { + throw std::invalid_argument("residual_suppressor: max_suppression_db must be positive"); + } + if (cfg.over_subtraction < Sample(1)) { + throw std::invalid_argument("residual_suppressor: over_subtraction must be >= 1"); + } + auto unit = [](Sample v) { return v >= Sample(0) && v < Sample(1); }; + if (!unit(cfg.leakage_smoothing) || !unit(cfg.gain_attack) || !unit(cfg.gain_release) + || !unit(cfg.floor_smoothing)) { + throw std::invalid_argument("residual_suppressor: smoothing constants must be in [0, 1)"); + } + if (!unit(cfg.leakage_gate)) { + throw std::invalid_argument("residual_suppressor: leakage_gate must be in [0, 1)"); + } + if (cfg.floor_window == 0 || cfg.floor_bias <= Sample(0)) { + throw std::invalid_argument("residual_suppressor: floor_window/floor_bias must be positive"); + } + return cfg; + } + + /// Packed-format bin access (Ooura layout: [DC, Nyquist, re, im, ...]). + static void packed_bin(const Sample* s, size_t k, size_t n, Sample& re, Sample& im) noexcept { + if (k == 0) { + re = s[0]; + im = Sample(0); + } + else if (k == n / 2) { + re = s[1]; + im = Sample(0); + } + else { + re = s[2 * k]; + im = s[2 * k + 1]; + } + } + + static void store_bin(Sample* s, size_t k, size_t n, Sample re, Sample im) noexcept { + if (k == 0) { + s[0] = re; + } + else if (k == n / 2) { + s[1] = re; + } + else { + s[2 * k] = re; + s[2 * k + 1] = im; + } + } + + /// Allocation-free noise source in [-1, 1] (xorshift32). + Sample next_noise() noexcept { + m_rng ^= m_rng << 13; + m_rng ^= m_rng >> 17; + m_rng ^= m_rng << 5; + return Sample(2) * (static_cast(m_rng) / Sample(4294967296.0)) - Sample(1); + } + + config m_cfg; + size_t m_n; ///< FFT size, 2 * block_size + basic_real_fft m_fft; + Sample m_g_min; + std::vector m_input; ///< sliding analysis window of e + std::vector m_ywin; ///< sliding analysis window of the echo estimate + std::vector m_yspec; ///< Hann-windowed estimation spectrum of the echo estimate + std::vector m_dspec; ///< Hann-windowed estimation spectrum of the mic d = e + yhat + std::vector m_window; + std::vector m_spec; + std::vector m_time; + std::vector m_sdy_re; ///< smoothed complex cross-spectrum D . conj(Yhat) + std::vector m_sdy_im; + std::vector m_sdd; ///< smoothed |D|^2 + std::vector m_syy; ///< smoothed |Yhat|^2 + std::vector m_see; ///< smoothed |E|^2 + std::vector m_leak; ///< per-bin leakage lam (residual per unit Yhat power) + std::vector m_gain; ///< smoothed per-bin gain + std::vector m_psm; ///< smoothed |E|^2 for the floor tracker + std::vector m_min_cur; ///< current-window minimum of m_psm + std::vector m_min_prev; ///< previous-window minimum + size_t m_min_count = 0; + std::vector m_gspec; ///< constrained gain spectrum (packed) + std::vector m_gtime; ///< gain impulse response workspace + std::uint32_t m_rng = 0x2545F491U; + }; + + /// The unit the compliance matrix measures: a linear canceller + /// followed by the residual suppressor, sharing one block size. + /// Same process_block(x, y, e) signature and real-time contract as + /// the cores themselves. + /// + /// The canceller is pluggable, and the default is the RAW + /// frequency-domain Kalman core — deliberately not pem_afc. PEM + /// prewhitening earns its keep in the CLOSED loop (feedback), where + /// the loudspeaker signal is correlated with the near end and a + /// naive update is biased; open-loop AEC has an exogenous far end, + /// no bias to fix, and the predictor's block-by-block refit just + /// injects gradient noise that floors the misalignment near -20 dB + /// (measured: raw fdkf single-talk residual -73 dBm0(A) where + /// PEM-Kalman plateaus at -29 on the same scenario, and double talk + /// barely moves the raw core — its per-bin noise-PSD tracker IS the + /// double-talk defense). Any type with the 4-argument raw-core + /// surface process_block(x, y, e, yhat) works, as does pem_afc's + /// 3-argument surface with echo_estimate_block(). + template > + class aec_chain { + public: + using canceller_type = Canceller; + + struct config { + typename Canceller::config canceller; + typename residual_suppressor::config postfilter; + }; + + explicit aec_chain(const config& cfg) + : m_afc(cfg.canceller) + , m_post(matched(cfg.postfilter, m_afc.block_size())) + , m_mid(m_afc.block_size()) + , m_yhat(m_afc.block_size()) {} + + size_t block_size() const noexcept { return m_afc.block_size(); } + + Canceller& canceller() noexcept { return m_afc; } + const Canceller& canceller() const noexcept { return m_afc; } + residual_suppressor& postfilter() noexcept { return m_post; } + const residual_suppressor& postfilter() const noexcept { return m_post; } + + void reset() noexcept { + m_afc.reset(); + m_post.reset(); + } + + void set_adaptation(bool enabled) noexcept { m_afc.set_adaptation(enabled); } + + void process_block(const Sample* x, const Sample* y, Sample* e) noexcept { + if constexpr (requires(Canceller& c) { c.process_block(x, y, e, e); }) { + m_afc.process_block(x, y, m_mid.data(), m_yhat.data()); + m_post.process_block(m_mid.data(), m_yhat.data(), e); + } + else { + m_afc.process_block(x, y, m_mid.data()); + m_post.process_block(m_mid.data(), m_afc.echo_estimate_block(), e); + } + } + + private: + static typename residual_suppressor::config matched(typename residual_suppressor::config pf, + size_t block_size) { + pf.block_size = block_size; // one block size for the chain + return pf; + } + + Canceller m_afc; + residual_suppressor m_post; + std::vector m_mid; + std::vector m_yhat; + }; + +} // namespace mutap diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3134be9..e716ff2 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -38,7 +38,8 @@ add_executable(mutap_tests test_closed_loop.cpp test_rir_fixtures.cpp test_aec.cpp - test_itu_signals.cpp) + test_itu_signals.cpp + test_postfilter.cpp) target_link_libraries(mutap_tests PRIVATE MuTap::MuTap mutap_warnings) diff --git a/tests/support/itu_levels.h b/tests/support/itu_levels.h index ebcb6b9..8d4c5d8 100644 --- a/tests/support/itu_levels.h +++ b/tests/support/itu_levels.h @@ -99,9 +99,14 @@ namespace mutap_test::itu { // from 1/(s + w), bilinear with PER-POLE frequency prewarping // (k_w = w / tan(w / 2fs)) so each pole lands at its analog // frequency — the plain transform shades the 12.2 kHz pole - // pair by ~1.5 dB at 10 kHz at audio rates. + // pair by ~1.5 dB at 10 kHz at audio rates. Prewarping is + // only defined below Nyquist; a pole above it (the 12.2 kHz + // pair at fs = 16 kHz) is clamped to 0.45 fs — it shapes + // nothing in band there, but tan() past pi/2 flips sign and + // the filter output went NaN (measured, 16 kHz chain tests). auto pole = [&](double w) { - const double kw = w / std::tan(w / (2.0 * fs)); + const double wc = std::min(w, 2.0 * pi * 0.45 * fs); + const double kw = wc / std::tan(wc / (2.0 * fs)); return std::pair{kw + w, w - kw}; // a0, a1 }; const auto [p1a0, p1a1] = pole(w1); diff --git a/tests/test_postfilter.cpp b/tests/test_postfilter.cpp new file mode 100644 index 0000000..dad9da0 --- /dev/null +++ b/tests/test_postfilter.cpp @@ -0,0 +1,576 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 MuTap contributors +// +// Stage 2 of the ITU compliance plan (docs/itu-compliance.md): the residual +// suppressor + comfort noise (postfilter.h) and the aec_chain it composes +// with the linear canceller. Every threshold here was measured first in the +// scratch harness and is asserted with margin; the matrix rows each test +// pre-proves are named in its comment (the full multi-rate compliance suite +// is Stage 3 — these pin the chain's headline behavior on the golden model). +// +// Chain under test: RAW partitioned frequency-domain Kalman canceller +// (transition 0.9998 and initial_uncertainty 10 — the AEC sweet spot +// measured in the Stage 2 sweeps: deeper steady state than the 0.9995 AFC +// default lifts the worst double-talk echo-loss band ~5 dB while 0.99995 +// costs double-talk transparency and tracking, and the larger P(0) buys +// ~4 dB of early convergence, 46.9 vs 45.2 dB ERL by 1200 ms) followed by +// the residual suppressor at its defaults. Note the chain deliberately +// does NOT use pem_afc: open-loop AEC has an exogenous far end, and PEM's +// predictor refit floors misalignment near -20 dB where the raw Kalman +// core reaches -75 dBm0(A) bare (measured; the full story is in +// postfilter.h's class comment). +// +// Measured (this file's exact protocols and chain config, double +// precision, block 256, 2048-tap unit-energy paths at 48 kHz unless +// stated): +// +// single-talk max residual, A-weighted 35 ms meter, CSS at -16 dBm0: +// kalman+pf cabin -79.9 dBm0(A) kalman+pf studio -88.9 +// nlms+pf cabin -65.0 kalman 16 kHz cabin -69.1 +// double-talk (AM-FM orthogonal pair, P.501 Table 7-6): +// send atten (send -1.7 dBPa): cabin 1.05 dB studio 0.93 +// echo loss worst band 200-6950 Hz (send -25.7 dBPa): +// cabin 38.0 dB (6660 Hz) studio 34.1 (270 Hz) +// convergence from cold start (CSS -16 dBm0, 35 ms meters): +// ERL 43.2 dB by 600 ms, 46.9 by 1200, 55.3 by 2000, 60.6 by 5000 +// comfort noise (Hoth at -46 dBm0 near end, CSS -16 dBm0 far end): +// level delta -1.20 dB; band spectrum deviation <= 1.69 dB; +// noise pumping 3.3 dB; DT-onset build-up 20.9 ms (gain_release +// 0.85; 0.9 measured 26.0 ms against the switching rows' 25 ms +// margin target, which is what set the default) +// +// Matrix margin targets asserted: ST residual < -64 dBm0(A) +// (ITU_EchoLevel), DT send atten <= 1.5 dB (ITU_DtSendAtten / +// ITU_DtSentSpeech), DT echo loss >= 33 dB per band (ITU_DtEchoLoss), +// ERL >= 40 dB by 600 ms and >= 46 dB by 1200 ms (ITU_ConvergenceQuiet), +// comfort noise +1/-2.5 dB level and half-mask spectrum +// (ITU_ComfortNoiseLevel/Spectrum), pumping <= 5 dB (ITU_NoisePump*), +// build-up <= 25 ms (switching rows). + +#include +#include +#include +#include + +#include + +#include "fixtures/rir_cabin.h" +#include "fixtures/rir_studio.h" +#include "mutap/fd_kalman.h" +#include "mutap/postfilter.h" +#include "support/echo_scenario.h" +#include "support/itu_levels.h" +#include "support/itu_signals.h" + +namespace { + + using namespace mutap_test; + + using chain_kalman = mutap::aec_chain; + using chain_nlms = mutap::aec_chain>; + + constexpr size_t k_block = 256; + constexpr size_t k_taps = 2048; + constexpr double k_fs = 48000.0; + + std::vector unit_energy(std::vector f) { + double e = 0.0; + for (double v : f) { + e += v * v; + } + for (auto& v : f) { + v /= std::sqrt(e); + } + return f; + } + + std::vector cabin_path() { + return unit_energy({fixtures::k_rir_cabin, fixtures::k_rir_cabin + k_taps}); + } + std::vector studio_path() { + return unit_energy({fixtures::k_rir_studio, fixtures::k_rir_studio + k_taps}); + } + + template + typename Chain::config aec_config(size_t block = k_block, size_t taps = k_taps) { + typename Chain::config cfg; + cfg.canceller.block_size = block; + cfg.canceller.partitions = taps / block; + if constexpr (requires { cfg.canceller.transition; }) { + cfg.canceller.transition = 0.9998; + cfg.canceller.initial_uncertainty = 10; + } + return cfg; + } + + /// Max of the A-weighted 35 ms level trace over [from, to). + double max_level_dbm0a(const std::vector& x, double fs, size_t from, size_t to) { + itu::a_weighting aw(fs); + auto xa = aw.apply(std::vector(x.begin() + static_cast(from), x.begin() + static_cast(to))); + itu::exp_level_meter m(fs, 0.035); + const auto tr = m.trace_dbm0(xa); + return *std::max_element(tr.begin() + static_cast(0.1 * fs), tr.end()); + } + + /// Single-talk protocol: CSS at -16 dBm0 through the path, no near + /// end; returns max A-weighted residual level over the last third. + template + double single_talk_residual(Proc& p, const std::vector& path, double fs, size_t block) { + typename echo_sim::config sc; + sc.echo_path = path; + sc.block_size = block; + echo_sim sim(sc); + + itu::css_config cc; + cc.periods = 60; + cc.shaped = true; + auto x = itu::make_css_at(cc, fs); + itu::set_level_dbm0(x, -16.0); + + std::vector out; + out.reserve(x.size()); + for (size_t blk = 0; blk + 1 <= x.size() / block; ++blk) { + sim.step(&x[blk * block], static_cast(nullptr), &p); + const auto& e = sim.error_block(); + out.insert(out.end(), e.begin(), e.end()); + } + return max_level_dbm0a(out, fs, out.size() * 2 / 3, out.size()); + } + + /// Double-talk protocol: converge on CSS, then AM-FM orthogonal pair — + /// far end on the receive plan at -16 dBm0 through the path, near end + /// on the send plan at send_dbpa. Returns the settled second half. + struct dt_run { + std::vector out; + std::vector x; + std::vector v; + size_t n = 0; + }; + + template + dt_run run_double_talk(Proc& p, const std::vector& path, double fs, double send_dbpa) { + typename echo_sim::config sc; + sc.echo_path = path; + sc.block_size = k_block; + echo_sim sim(sc); + + itu::css_config cc; + cc.periods = 30; + cc.shaped = true; + auto xc = itu::make_css_at(cc, fs); + itu::set_level_dbm0(xc, -16.0); + for (size_t blk = 0; blk + 1 <= xc.size() / k_block; ++blk) { + sim.step(&xc[blk * k_block], static_cast(nullptr), &p); + } + + dt_run r; + r.n = static_cast(8 * fs); + r.x = itu::make_amfm(itu::amfm_receive_plan(), r.n, fs); + itu::set_level_dbm0(r.x, -16.0); + r.v = itu::make_amfm(itu::amfm_send_plan(), r.n, fs); + const double vg = itu::dbpa_to_rms(send_dbpa) / itu::rms_of(r.v.data(), r.v.size()); + for (auto& s : r.v) { + s *= vg; + } + r.out.reserve(r.n); + for (size_t blk = 0; blk + 1 <= r.n / k_block; ++blk) { + sim.step(&r.x[blk * k_block], &r.v[blk * k_block], &p); + const auto& e = sim.error_block(); + r.out.insert(r.out.end(), e.begin(), e.end()); + } + return r; + } + + std::vector settled_half(const std::vector& s, size_t n) { + return {s.begin() + static_cast(n / 2), s.begin() + static_cast(n)}; + } + + // ---------------------------------------------------------------- tests + + // ITU_EchoLevel margin target: < -64 dBm0(A) single-talk residual. + // Measured: cabin -79.9, studio -88.9 dBm0(A). + TEST(ItuChain, SingleTalkResidualKalman) { + { + chain_kalman c(aec_config()); + EXPECT_LT(single_talk_residual(c, cabin_path(), k_fs, k_block), -72.0); + } + { + chain_kalman c(aec_config()); + EXPECT_LT(single_talk_residual(c, studio_path(), k_fs, k_block), -80.0); + } + } + + // The NLMS-core chain also meets the margin target (measured -65.0). + TEST(ItuChain, SingleTalkResidualNlmsCore) { + chain_nlms c(aec_config()); + EXPECT_LT(single_talk_residual(c, cabin_path(), k_fs, k_block), -64.0); + } + + // 16 kHz is a REQUIRED operating rate (matrix sample-rate policy). + // Cabin path resampled with the NOTE 2 resampler; measured -69.1. + TEST(ItuChain, SingleTalkResidualAt16k) { + const double fs16 = 16000.0; + std::vector cab48(fixtures::k_rir_cabin, fixtures::k_rir_cabin + 4096); + auto cab16 = itu::resample(cab48, static_cast(fixtures::k_rir_cabin_fs), fs16, 0); + cab16.resize(1024); // 64 ms covers the cabin RT60 + chain_kalman c(aec_config(128, 1024)); + EXPECT_LT(single_talk_residual(c, unit_energy(cab16), fs16, 128), -64.0); + } + + // ITU_DtSendAtten / ITU_DtSentSpeech margin target: <= 1.5 dB near-end + // attenuation during double talk, send at -1.7 dBPa. Measured: cabin + // 1.05, studio 0.93 dB (bare canceller -0.41: the postfilter's whole + // cost is ~1.4 dB, spent at the comb bins the analysis cannot split). + TEST(ItuChain, DoubleTalkSendAttenuation) { + const auto sp = itu::amfm_send_plan(); + for (const auto& path : {cabin_path(), studio_path()}) { + chain_kalman c(aec_config()); + const auto r = run_double_talk(c, path, k_fs, -1.7); + const auto vh = settled_half(r.v, r.n); + const auto oh = settled_half(r.out, r.n); + const double atten = + itu::comb_band_level_db(vh, sp, 0.0, k_fs) - itu::comb_band_level_db(oh, sp, 0.0, k_fs); + EXPECT_LE(atten, 1.5); + } + } + + // ITU_DtEchoLoss margin target: >= 33 dB echo loss during double talk + // in EACH receive band 200-6950 Hz, send at -25.7 dBPa. Measured worst + // bands: cabin 38.0 dB (6660 Hz), studio 34.1 dB (the 270 Hz line the + // gain constraint cannot notch selectively — the canceller carries it). + TEST(ItuChain, DoubleTalkEchoLossPerBand) { + const auto rp = itu::amfm_receive_plan(); + const double worst_expected[] = {35.0, 33.0}; // cabin, studio + size_t room = 0; + for (const auto& path : {cabin_path(), studio_path()}) { + chain_kalman c(aec_config()); + const auto r = run_double_talk(c, path, k_fs, -25.7); + const auto xh = settled_half(r.x, r.n); + const auto oh = settled_half(r.out, r.n); + double worst = 1e9; + for (size_t b = 0; b < rp.f0.size(); ++b) { + if (rp.f0[b] < 200.0 || rp.f0[b] > 6950.0) { + continue; + } + itu::amfm_plan one; + one.f0.push_back(rp.f0[b]); + one.df.push_back(rp.df[b]); + const double loss = + itu::comb_band_level_db(xh, one, 0.0, k_fs) - itu::comb_band_level_db(oh, one, 0.0, k_fs); + worst = std::min(worst, loss); + } + EXPECT_GE(worst, worst_expected[room++]); + } + } + + // ITU_ConvergenceQuiet margin targets: ERL >= 40 dB by 600 ms and + // >= 46 dB by 1200 ms from cold start. Measured 43.2 / 46.9 dB. + TEST(ItuChain, ConvergenceFromColdStart) { + chain_kalman c(aec_config()); + typename echo_sim::config sc; + sc.echo_path = cabin_path(); + sc.block_size = k_block; + echo_sim sim(sc); + + itu::css_config cc; + cc.periods = 15; + cc.shaped = true; + auto x = itu::make_css_at(cc, k_fs); + itu::set_level_dbm0(x, -16.0); + + std::vector out; + std::vector mic; + for (size_t blk = 0; blk + 1 <= x.size() / k_block; ++blk) { + sim.step(&x[blk * k_block], static_cast(nullptr), &c); + const auto& e = sim.error_block(); + out.insert(out.end(), e.begin(), e.end()); + const auto& y = sim.echo_block(); + mic.insert(mic.end(), y.begin(), y.end()); + } + itu::exp_level_meter m_y(k_fs, 0.035); + itu::exp_level_meter m_e(k_fs, 0.035); + const auto tr_y = m_y.trace_dbm0(mic); + const auto tr_e = m_e.trace_dbm0(out); + // ERL by t: best reading within the preceding CSS period (350 ms) + // — an instantaneous read lands wherever the CSS pause happens to + // fall and measures meter decay, not echo. + const auto erl_by = [&](double t) { + const size_t i0 = static_cast((t - 0.35) * k_fs); + const size_t i1 = static_cast(t * k_fs); + double best = -1e9; + for (size_t i = i0; i < i1; ++i) { + best = std::max(best, tr_y[i] - tr_e[i]); + } + return best; + }; + EXPECT_GE(erl_by(0.6), 40.0); + EXPECT_GE(erl_by(1.2), 46.0); + } + + // ITU_ComfortNoiseLevel (+1/-2.5 dB target), ITU_ComfortNoiseSpectrum + // (half-mask: +-6 dB 200-800 Hz, +-5 to 2 kHz, +-3 above) and + // ITU_NoisePump* (<= 5 dB target): Hoth noise at -46 dBm0 in the near + // end, far-end CSS bursts. Measured: level delta -1.20 dB, worst band + // deviation 1.69 dB, pumping 3.3 dB. + TEST(ItuChain, ComfortNoiseMatchesFloor) { + chain_kalman c(aec_config()); + typename echo_sim::config sc; + sc.echo_path = cabin_path(); + sc.block_size = k_block; + echo_sim sim(sc); + + const size_t n_half = static_cast(10 * k_fs); + auto noise = itu::make_hoth_noise(2 * n_half, 7, k_fs); + itu::set_level_dbm0(noise, -46.0); + itu::css_config cc; + cc.periods = 28; + cc.shaped = true; + auto x = itu::make_css_at(cc, k_fs); + itu::set_level_dbm0(x, -16.0); + x.resize(2 * n_half, 0.0); // second half: far end silent + + std::vector out; + for (size_t blk = 0; blk + 1 <= x.size() / k_block; ++blk) { + sim.step(&x[blk * k_block], &noise[blk * k_block], &c); + const auto& e = sim.error_block(); + out.insert(out.end(), e.begin(), e.end()); + } + const auto seg = [&](double t0, double t1) { + return std::vector(out.begin() + static_cast(t0 * k_fs), + out.begin() + static_cast(t1 * k_fs)); + }; + auto talk = seg(6.0, 9.5); // far end active: suppressed + comfort fill + auto quiet = seg(16.0, 19.5); // far end silent: the true noise floor + + // Level match, A-weighted. + itu::a_weighting aw(k_fs); + auto ta = aw.apply(talk); + auto qa = aw.apply(quiet); + const double lt = itu::level_dbov(ta.data(), ta.size()) - itu::k_dbov_per_dbm0; + const double lq = itu::level_dbov(qa.data(), qa.size()) - itu::k_dbov_per_dbm0; + EXPECT_LE(lt - lq, 1.0); + EXPECT_GE(lt - lq, -2.5); + + // Spectral match: Welch band levels, half-mask bounds. + const auto band_levels = [&](const std::vector& s) { + const size_t n = 8192; + mutap::basic_real_fft fft(n); + std::vector psd(n / 2 + 1, 0.0); + std::vector buf(n); + std::vector win(n); + for (size_t i = 0; i < n; ++i) { + win[i] = + 0.5 - 0.5 * std::cos(2.0 * std::numbers::pi * static_cast(i) / static_cast(n - 1)); + } + size_t frames = 0; + for (size_t off = 0; off + n <= s.size(); off += n / 2, ++frames) { + for (size_t i = 0; i < n; ++i) { + buf[i] = win[i] * s[off + i]; + } + fft.forward_inplace(buf.data()); + psd[0] += buf[0] * buf[0]; + psd[n / 2] += buf[1] * buf[1]; + for (size_t k = 1; k < n / 2; ++k) { + psd[k] += buf[2 * k] * buf[2 * k] + buf[2 * k + 1] * buf[2 * k + 1]; + } + } + const double edges[] = {200, 400, 800, 1600, 3150, 6300, 8000}; + std::vector bands; + for (size_t b = 0; b + 1 < 7; ++b) { + double sum = 0.0; + const size_t k0 = static_cast(edges[b] * n / k_fs); + const size_t k1 = static_cast(edges[b + 1] * n / k_fs); + for (size_t k = k0; k < k1; ++k) { + sum += psd[k]; + } + bands.push_back(10.0 * std::log10(sum / static_cast(frames) + 1e-30)); + } + return bands; + }; + const auto bt = band_levels(talk); + const auto bq = band_levels(quiet); + const double mask[] = {6.0, 6.0, 5.0, 3.0, 3.0, 3.0}; + for (size_t b = 0; b < 6; ++b) { + EXPECT_LE(std::abs(bt[b] - bq[b]), mask[b]) << "band " << b; + } + + // Noise pumping over the far-end-active segment. + itu::exp_level_meter m(k_fs, 0.035); + const auto tr = m.trace_dbm0(ta); + const double vmax = *std::max_element(tr.begin() + 1680, tr.end()); + const double vmin = *std::min_element(tr.begin() + 1680, tr.end()); + EXPECT_LE(vmax - vmin, 5.0); + } + + // Switching build-up target <= 25 ms: near-end onset mid-double-talk + // reaches within 3 dB of its settled send level in 20.9 ms (measured). + TEST(ItuChain, NearEndBuildUpTime) { + chain_kalman c(aec_config()); + typename echo_sim::config sc; + sc.echo_path = cabin_path(); + sc.block_size = k_block; + echo_sim sim(sc); + + itu::css_config cc; + cc.periods = 30; + cc.shaped = true; + auto x = itu::make_css_at(cc, k_fs); + itu::set_level_dbm0(x, -16.0); + + itu::css_config cd; + cd.periods = 30; + cd.kind = itu::css_kind::double_talk; + cd.shaped = true; + auto v = itu::make_css_at(cd, k_fs); + v.resize(x.size(), 0.0); + const double vg = itu::dbpa_to_rms(-1.7) / itu::rms_of(v.data(), v.size() / 2); + // Onset at 2/3 of the run — the settled-level window below needs + // 2 s of trace after it. (An onset at 5/6 left only 1.75 s and the + // window read past the trace: ASan aborted, MSVC compared against + // a garbage median. The other platforms' passes were luck.) + const size_t t0 = (x.size() / k_block * 2 / 3) * k_block; + std::vector vv(x.size(), 0.0); + for (size_t i = t0; i < x.size(); ++i) { + vv[i] = vg * v[i - t0]; + } + std::vector out; + for (size_t blk = 0; blk + 1 <= x.size() / k_block; ++blk) { + sim.step(&x[blk * k_block], &vv[blk * k_block], &c); + const auto& e = sim.error_block(); + out.insert(out.end(), e.begin(), e.end()); + } + itu::a_weighting aw(k_fs); + auto seg = aw.apply(std::vector(out.begin() + static_cast(t0) - 4800, out.end())); + itu::exp_level_meter m(k_fs, 0.005); + const auto tr = m.trace_dbm0(seg); + ASSERT_GE(tr.size(), 4800 + static_cast(2 * k_fs)); + std::vector settled(tr.begin() + 4800 + static_cast(k_fs), + tr.begin() + 4800 + static_cast(2 * k_fs)); + std::nth_element(settled.begin(), settled.begin() + static_cast(settled.size() / 2), settled.end()); + const double target = settled[settled.size() / 2]; + size_t t_hit = tr.size(); + for (size_t i = 4800; i < tr.size(); ++i) { + if (tr[i] >= target - 3.0) { + t_hit = i - 4800; + break; + } + } + EXPECT_LE(1000.0 * static_cast(t_hit) / k_fs, 25.0); + } + + // The pem_afc 3-argument surface still composes (the closed-loop AFC + // heritage path through aec_chain's echo_estimate_block() branch). + TEST(ItuChain, PemAfcCancellerStillComposes) { + using chain_pem = mutap::aec_chain>; + chain_pem::config cfg; + cfg.canceller.fdaf.block_size = k_block; + cfg.canceller.fdaf.partitions = k_taps / k_block; + chain_pem c(cfg); + // Composition smoke test, not a depth claim: the PEM structure's + // predictor refit floors misalignment near -20 dB (the reason the + // AEC chain defaults to the raw core). Measured 16.5 dB below the + // raw echo, converged. + typename echo_sim::config sc; + sc.echo_path = cabin_path(); + sc.block_size = k_block; + echo_sim sim(sc); + itu::css_config cc; + cc.periods = 20; + cc.shaped = true; + auto x = itu::make_css_at(cc, k_fs); + itu::set_level_dbm0(x, -16.0); + std::vector out; + std::vector mic; + for (size_t blk = 0; blk + 1 <= x.size() / k_block; ++blk) { + sim.step(&x[blk * k_block], static_cast(nullptr), &c); + const auto& e = sim.error_block(); + out.insert(out.end(), e.begin(), e.end()); + const auto& y = sim.echo_block(); + mic.insert(mic.end(), y.begin(), y.end()); + } + const size_t from = out.size() * 2 / 3; + const double le = itu::level_dbov(out.data() + from, out.size() - from); + const double ly = itu::level_dbov(mic.data() + from, mic.size() - from); + EXPECT_LT(le, ly - 12.0); + } + + // reset() restores the initial state exactly, comfort-noise PRNG + // included: two identical runs produce identical output. + TEST(ItuChain, ResetRestoresDeterminism) { + chain_kalman c(aec_config()); + auto run_once = [&] { + typename echo_sim::config sc; + sc.echo_path = cabin_path(); + sc.block_size = k_block; + echo_sim sim(sc); + itu::css_config cc; + cc.periods = 4; + cc.shaped = true; + auto x = itu::make_css_at(cc, k_fs); + itu::set_level_dbm0(x, -16.0); + std::vector out; + for (size_t blk = 0; blk + 1 <= x.size() / k_block; ++blk) { + sim.step(&x[blk * k_block], static_cast(nullptr), &c); + const auto& e = sim.error_block(); + out.insert(out.end(), e.begin(), e.end()); + } + return out; + }; + const auto first = run_once(); + c.reset(); + const auto second = run_once(); + ASSERT_EQ(first.size(), second.size()); + for (size_t i = 0; i < first.size(); ++i) { + ASSERT_EQ(first[i], second[i]) << "diverged at sample " << i; + } + } + + TEST(PostFilterConfigValidation, RejectsBadConfigs) { + using pf = mutap::residual_suppressor; + pf::config good; + EXPECT_NO_THROW(pf{good}); + { + auto c = good; + c.block_size = 100; // not a power of two + EXPECT_THROW(pf{c}, std::invalid_argument); + } + { + auto c = good; + c.analysis_blocks = 2; // < 4 + EXPECT_THROW(pf{c}, std::invalid_argument); + } + { + auto c = good; + c.max_suppression_db = 0.0; + EXPECT_THROW(pf{c}, std::invalid_argument); + } + { + auto c = good; + c.over_subtraction = 0.5; // < 1 + EXPECT_THROW(pf{c}, std::invalid_argument); + } + { + auto c = good; + c.leakage_smoothing = 1.0; // not in [0, 1) + EXPECT_THROW(pf{c}, std::invalid_argument); + } + { + auto c = good; + c.floor_window = 0; + EXPECT_THROW(pf{c}, std::invalid_argument); + } + } + + TEST(PostFilterRtContract, ProcessingPathIsNoexcept) { + using pf = mutap::residual_suppressor; + static_assert(noexcept(std::declval().process_block(nullptr, nullptr, nullptr))); + static_assert(noexcept(std::declval().reset())); + static_assert(noexcept(std::declval().gain_at(0))); + static_assert(noexcept(std::declval().coherence_at(0))); + static_assert(noexcept(std::declval().leakage_at(0))); + static_assert(noexcept(std::declval().process_block(nullptr, nullptr, nullptr))); + static_assert(noexcept(std::declval().reset())); + SUCCEED(); + } + +} // namespace