blip25-mbe 0.1.0

//! MBE baseline codec — TIA-102.BABA-A §1.10–§1.13 (synthesis pipeline)
//! and §1.8 (dequantize-side analysis is shared with the upstream pipeline in
//! [`crate::imbe_wire::dequantize`]).
//!
//! This module hosts the synthesizer that turns [`MbeParams`] into
//! 8 kHz 16-bit PCM, plus a (currently unimplemented) forward analysis
//! encoder at [`analysis`] that would take 8 kHz PCM back to [`MbeParams`].
//!
//! Synthesis pipeline (1993 P25 vocoder specification):
//!
//! - **§1.10 enhancement** — `M̃_l → M̄_l` via the W_l psycho-acoustic
//!   weighting.
//! - **§1.11 smoothing** — frame repeat / mute / V/UV smoothing.
//! - **§1.12.1 unvoiced** — LCG noise → DFT → per-band spectral shaping
//!   → IDFT → weighted overlap-add.
//! - **§1.12.2 voiced** — sinusoidal overlap-add with phase tracking
//!   across the four V/UV transition cases.
//!
//! Per TIA-102.BABG this baseline generation explicitly cannot pass the
//! enhanced-vocoder performance tests. It is implemented here for spec
//! completeness and as the reference point against which AMBE/+/+2
//! generations' improvements are measured.
//!
//! [`MbeParams`]: crate::mbe_params::MbeParams

use core::f64::consts::PI as PI64;

use crate::mbe_params::{L_MAX, MbeParams};

pub mod analysis;
pub mod phase_regen;

include!(concat!(env!("OUT_DIR"), "/annex_i_synth_window.rs"));

/// Lookup `wS(n)` for `n ∈ [-105, 105]`. Returns 0.0 outside that range
/// (matching the spec — `wS` is treated as zero outside its support).
#[inline]
pub fn synth_window(n: i32) -> f32 {
    if !(-105..=105).contains(&n) {
        return 0.0;
    }
    IMBE_SYNTH_WINDOW[(n + 105) as usize]
}

// ---------------------------------------------------------------------------
// §1.10 — Spectral amplitude enhancement (Eq. 105–111)
// ---------------------------------------------------------------------------

/// Initial value for the local-energy state `S_E` per BABA-A §10
/// Annex A. Used when the synthesizer is cold-started.
pub const INIT_S_E: f64 = 75000.0;

/// Lower bound on `S_E` per Eq. 111. Floors out the recurrence so very
/// quiet frames don't drive the predictor toward zero.
pub const S_E_FLOOR: f64 = 10000.0;

/// Apply BABA-A §1.10 spectral amplitude enhancement to the
/// reconstructed magnitudes `M̃_l`, producing the synthesizer-side
/// magnitudes `M̄_l` and the updated local-energy state `S_E(0)`.
///
/// Per §1.10:
/// 1. Compute energy moments `R_M0`, `R_M1` (Eq. 105–106).
/// 2. Per-harmonic weight `W_l` (Eq. 107) with `8·l ≤ L̃` bypass and
///    `[0.5, 1.2]` clamp on `W_l` (Eq. 108).
/// 3. Energy-preserving rescale by `γ = sqrt(R_M0 / Σ M̄_l²)` (Eq. 109–110).
/// 4. Update `S_E(0) = max(0.95·S_E(−1) + 0.05·R_M0, 10000)` (Eq. 111).
///
/// All internal arithmetic is `f64` per the spec's reference C
/// implementation. Inputs and outputs at the `MbeParams` boundary
/// stay in `f32`.
///
/// `m_tilde` must have length `L̃ ≥ 1`; the returned array's entries
/// `[0..L̃]` are populated, the rest are zero.
///
/// **Edge cases.** If `R_M0` is zero (silence frame) the enhancement
/// reduces to a passthrough — `W_l` is undefined (denominator = 0)
/// and `γ` is also undefined; we return `M̄_l = M̃_l` and only update
/// `S_E`. Similarly for `R_M0 == R_M1` (denominator zero by another
/// path).
pub fn enhance_spectral_amplitudes(
    m_tilde: &[f32],
    omega_0: f32,
    s_e_prev: f64,
) -> ([f32; L_MAX as usize], f64) {
    let l = m_tilde.len();
    debug_assert!(l > 0 && l <= L_MAX as usize);
    let omega_0 = f64::from(omega_0);

    // Eq. 105–106: energy moments.
    let mut r_m0 = 0.0f64;
    let mut r_m1 = 0.0f64;
    for (i, &m) in m_tilde.iter().enumerate() {
        let l_one = (i + 1) as f64;
        let m2 = f64::from(m) * f64::from(m);
        r_m0 += m2;
        r_m1 += m2 * (omega_0 * l_one).cos();
    }

    let mut m_bar_f64 = [0f64; L_MAX as usize];
    let denom = omega_0 * r_m0 * (r_m0 - r_m1);

    if r_m0 <= 0.0 || denom.abs() < 1e-30 {
        // Silence-frame degenerate path — bypass enhancement.
        for (i, &m) in m_tilde.iter().enumerate() {
            m_bar_f64[i] = f64::from(m);
        }
    } else {
        // Eq. 107–108: per-harmonic weight with bypass + clamp.
        for (i, &m) in m_tilde.iter().enumerate() {
            let l_one = (i + 1) as f64;
            let bar = if 8 * (i + 1) <= l {
                f64::from(m) // first branch — bypass for the lowest ⌊L̃/8⌋
            } else {
                let num = r_m0 * r_m0 + r_m1 * r_m1
                    - 2.0 * r_m0 * r_m1 * (omega_0 * l_one).cos();
                if num <= 0.0 {
                    f64::from(m) // numerator pathology → passthrough
                } else {
                    let w_l = f64::from(m).sqrt() * (0.96 * num / denom).powf(0.25);
                    if w_l > 1.2 {
                        1.2 * f64::from(m)
                    } else if w_l < 0.5 {
                        0.5 * f64::from(m)
                    } else {
                        w_l * f64::from(m)
                    }
                }
            };
            m_bar_f64[i] = bar;
        }

        // Eq. 109–110: energy-preserving rescale.
        let mut sum_sq = 0.0f64;
        for &v in m_bar_f64.iter().take(l) {
            sum_sq += v * v;
        }
        if sum_sq > 1e-30 {
            let gamma = (r_m0 / sum_sq).sqrt();
            for v in m_bar_f64.iter_mut().take(l) {
                *v *= gamma;
            }
        }
    }

    // Eq. 111: S_E recurrence with floor.
    let s_e = (0.95 * s_e_prev + 0.05 * r_m0).max(S_E_FLOOR);

    let mut m_bar = [0f32; L_MAX as usize];
    for (i, &v) in m_bar_f64.iter().take(l).enumerate() {
        m_bar[i] = v as f32;
    }
    (m_bar, s_e)
}

// ---------------------------------------------------------------------------
// §1.11 — Adaptive smoothing, frame repeat / mute, V/UV smoothing
// ---------------------------------------------------------------------------

/// Initial value for the amplitude-smoothing state `τ_M` per BABA-A
/// §10 Annex A.
pub const INIT_TAU_M: f64 = 20480.0;

/// Smoothed-error-rate threshold above which the frame is muted
/// (full-rate IMBE §1.11.2 mute condition).
pub const MUTE_EPSILON_R_THRESHOLD: f64 = 0.0875;

/// Smoothed-error-rate threshold above which a half-rate AMBE+2 frame
/// is muted (BABA-A §2.8.3 Eq. 199). Different from full-rate's 0.0875
/// because the half-rate ε_R recurrence (Eq. 197, weight 0.001064)
/// scales differently than the full-rate recurrence.
pub const MUTE_EPSILON_R_THRESHOLD_HALFRATE: f64 = 0.096;

/// Gain applied to the unvoiced LCG output when synthesizing the mute
/// frame's comfort-noise output. §1.11.2 calls for "random
/// small-amplitude noise"; we reuse the unvoiced synth's existing
/// noise window (centered, zero-mean after subtracting `LCG_MEAN`).
/// Calibrated so the i16 RMS sits around 100 (≈ 0.003 × 32768) — the
/// same low level JMBE / SDRTrunk emit on mute.
const MUTE_NOISE_GAIN: f64 = 0.0065;
/// Mean of the §1.12.1 LCG output range `[0, 53124]`. Subtracted from
/// each LCG sample to produce zero-mean noise.
const LCG_MEAN: f64 = 26562.0;

/// Default fundamental frequency used as the substitute when the
/// beyond-spec consecutive-repeat reset (`SynthState::set_repeat_reset_after`)
/// fires. Matches JMBE / SDRTrunk's `IMBEFundamentalFrequency.DEFAULT`
/// (b̂₀ = 134, ω₀ ≈ 0.2985·π ≈ 0.937 rad/sample, F0 ≈ 119 Hz). Not a
/// spec value — gap 0022 resolution confirms BABA-A §7.7 has no
/// reset behavior.
const DEFAULT_OMEGA_0: f32 = 0.937_544_4;

/// Per-frame error context derived from FEC decoding (§1.5) plus the
/// pitch-validity check (§1.3.1). Drives the §1.11 smoothing
/// decisions.
#[derive(Clone, Copy, Debug, Default)]
pub struct FrameErrorContext {
    /// Bit errors corrected in û₀ (Golay).
    pub epsilon_0: u8,
    /// Bit errors corrected in û₄ specifically (used in V_M branch 2).
    pub epsilon_4: u8,
    /// Total bit errors across all 7 FEC-protected vectors `û₀..û₆`.
    pub epsilon_t: u8,
    /// Set when the decoded pitch index `b̂₀` is in the reserved range
    /// `[208, 255]` (or otherwise invalid) — forces a frame repeat
    /// regardless of the error counts.
    pub bad_pitch: bool,
}

/// Action the synthesizer should take for the current frame, per
/// BABA-A §1.11.1 / §1.11.2.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub enum FrameDisposition {
    /// Synthesize from the current frame's parameters.
    Use,
    /// Substitute the previous frame's parameters (Eq. 99–104) and
    /// then synthesize.
    Repeat,
    /// Bypass the synthesizer and emit silence / low-amplitude noise.
    Mute,
}

/// Update the smoothed error rate `ε_R(0) = 0.95·ε_R(−1) + 0.05·(ε_T/144)`
/// and decide the frame disposition (full-rate IMBE).
///
/// Per §1.11.1 the repeat trigger is:
/// `b̂₀ invalid` **OR** (`ε₀ ≥ 2` AND `ε_T ≥ 10 + 40·ε_R(0)`).
///
/// Per §1.11.2 mute supersedes repeat when `ε_R(0) > 0.0875`.
pub fn frame_disposition(
    err: &FrameErrorContext,
    epsilon_r_prev: f64,
) -> (FrameDisposition, f64) {
    let epsilon_r =
        0.95 * epsilon_r_prev + 0.05 * (f64::from(err.epsilon_t) / 144.0);
    let disp = if epsilon_r > MUTE_EPSILON_R_THRESHOLD {
        FrameDisposition::Mute
    } else if err.bad_pitch
        || (err.epsilon_0 >= 2
            && f64::from(err.epsilon_t) >= 10.0 + 40.0 * epsilon_r)
    {
        FrameDisposition::Repeat
    } else {
        FrameDisposition::Use
    };
    (disp, epsilon_r)
}

/// Half-rate AMBE+2 counterpart of [`frame_disposition`] per BABA-A
/// §2.8.1–§2.8.3 (Eq. 196–199).
///
/// Differs from full-rate in four places:
///
/// - `ε_T = ε₀ + ε₁` (Eq. 196) — only the two Golay codewords contribute,
///   not all four cosets. The caller is responsible for passing
///   `err.epsilon_t = ε₀ + ε₁`; this function does not re-derive it.
/// - ε_R recurrence weight is `0.001064` (Eq. 197) rather than `0.05/144`.
///   Per §2.8.1 commentary this matches `0.05 / 47 ≈ 0.001064` scaled by
///   the half-rate frame size.
/// - Repeat trigger is `(ε₀ ≥ 4) OR (ε₀ ≥ 2 AND ε_T ≥ 6)` (Eq. 198-199) —
///   constant 6 threshold, not `10 + 40·ε_R`.
/// - Mute threshold is `ε_R(0) > 0.096` (§2.8.3) rather than 0.0875.
pub fn frame_disposition_halfrate(
    err: &FrameErrorContext,
    epsilon_r_prev: f64,
) -> (FrameDisposition, f64) {
    // Eq. 197: ε_R(0) = 0.95·ε_R(−1) + 0.001064·ε_T (no /144 divisor).
    let epsilon_r = 0.95 * epsilon_r_prev + 0.001064 * f64::from(err.epsilon_t);
    let disp = if epsilon_r > MUTE_EPSILON_R_THRESHOLD_HALFRATE {
        FrameDisposition::Mute
    } else if err.bad_pitch
        || err.epsilon_0 >= 4
        || (err.epsilon_0 >= 2 && err.epsilon_t >= 6)
    {
        FrameDisposition::Repeat
    } else {
        FrameDisposition::Use
    };
    (disp, epsilon_r)
}

/// Output of the §1.11.3 V/UV-and-amplitude smoothing step.
#[derive(Clone, Debug)]
pub struct SmoothedFrame {
    /// Smoothed per-harmonic spectral amplitudes (`M̄_l` after γ_M scaling).
    /// Entries `[0..L]` are populated.
    pub m_bar: [f32; L_MAX as usize],
    /// Smoothed per-harmonic V/UV decisions (`v̄_l`). Entries `[0..L]`
    /// are populated; higher indices are `false`.
    pub v_bar: [bool; L_MAX as usize],
    /// Updated amplitude-smoothing state `τ_M(0)` for the next frame.
    pub tau_m: f64,
}

/// Apply the §1.11.3 V/UV-and-amplitude smoothing pipeline.
///
/// Inputs:
/// - `m_bar`: enhanced spectral amplitudes from §1.10 (`M̄_l`).
/// - `v_tilde`: per-harmonic V/UV from the §1.3.2 expansion (`ṽ_l`).
/// - `s_e`: current-frame `S_E(0)` from §1.10.
/// - `epsilon_r`: smoothed error rate `ε_R(0)` from
///   [`frame_disposition`].
/// - `epsilon_t`, `epsilon_4`: per-frame error counts (raw, not smoothed).
/// - `tau_m_prev`: previous frame's `τ_M(−1)`.
///
/// Returns the smoothed `M̄_l`, the smoothed per-harmonic V/UV `v̄_l`,
/// and the updated `τ_M(0)`.
pub fn apply_smoothing(
    m_bar: &[f32],
    v_tilde: &[bool],
    s_e: f64,
    epsilon_r: f64,
    epsilon_t: u8,
    epsilon_4: u8,
    tau_m_prev: f64,
) -> SmoothedFrame {
    let l = m_bar.len();
    debug_assert_eq!(v_tilde.len(), l);

    // Eq. 112: V_M threshold (3-branch).
    let v_m: f64 = if epsilon_r <= 0.005 && epsilon_t <= 4 {
        f64::INFINITY
    } else if epsilon_r <= 0.0125 && epsilon_4 == 0 {
        45.255 * s_e.powf(0.375) * (-277.26 * epsilon_r).exp()
    } else {
        1.414 * s_e.powf(0.375)
    };

    // Eq. 113: per-harmonic V/UV override (force voiced if M̄_l > V_M).
    let mut v_bar = [false; L_MAX as usize];
    for (i, &m) in m_bar.iter().enumerate() {
        v_bar[i] = if f64::from(m) > v_m { true } else { v_tilde[i] };
    }

    // Eq. 114: total amplitude.
    let a_m: f64 = m_bar.iter().take(l).map(|&v| f64::from(v)).sum();

    // Eq. 115: τ_M recurrence.
    let tau_m = if epsilon_r <= 0.005 && epsilon_t <= 6 {
        20480.0
    } else {
        6000.0 - 300.0 * f64::from(epsilon_t) + tau_m_prev
    };

    // Eq. 116: γ_M scaling.
    let gamma_m: f64 = if tau_m > a_m { 1.0 } else if a_m > 0.0 { tau_m / a_m } else { 1.0 };

    let mut m_bar_out = [0f32; L_MAX as usize];
    for (i, &m) in m_bar.iter().enumerate() {
        m_bar_out[i] = ((f64::from(m)) * gamma_m) as f32;
    }

    SmoothedFrame { m_bar: m_bar_out, v_bar, tau_m }
}

// ---------------------------------------------------------------------------
// §1.12.1 — Unvoiced synthesis (Eq. 117–126)
// ---------------------------------------------------------------------------

/// Samples per 20 ms frame at 8 kHz (`N` in §1.12).
pub const FRAME_SAMPLES: usize = 160;

/// Spec-conformant initial state of the white-noise LCG (`u(−105) = 3147`
/// per BABA-A §10 Annex A). Used when [`UnvoicedNoiseGen::SpecLcg`]
/// is selected. The DVSI hardware does NOT use the spec generator on
/// the half-rate (AMBE+2) path — see [`UnvoicedNoiseGen::ChipLcg`].
pub const NOISE_INIT: u32 = 3147;

/// Unvoiced noise generator selector.
///
/// BABA-A §1.12.1 Eq. 117 specifies the LCG
/// `u(n+1) = (171·u(n) + 11213) mod 53125`, seed `u(−105) = 3147`. The
/// IMBE full-rate path on the DVSI chip matches this generator.
///
/// On the **half-rate (AMBE+2) path** the chip uses a different LCG —
/// `(173·u(n) + 13849) mod 65536` with an internal seed equivalent to
/// 60 584 at our `t=0 sample 0` — verified by Probe 1 of the chip-vs-spec
/// noise investigation (gap report `0025_chip_uses_pn_lcg_for_halfrate`).
/// That is the recurrence BABA-A §1.5 Eq. 84-85 specifies for **FEC
/// masking**, repurposed by DVSI for the AMBE+2 unvoiced synth.
///
/// Default: [`UnvoicedNoiseGen::ChipLcg`] — match the de-facto industry
/// reference. Switch to [`UnvoicedNoiseGen::SpecLcg`] for BABA-A §1.12.1
/// spec-literal output (useful for
/// conformance testing against the spec rather than the chip).
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub enum UnvoicedNoiseGen {
    /// `(171, 11213, 53125)` LCG with seed `3147`. BABA-A §1.12.1 Eq. 117.
    /// JMBE 1.0.9 uses this. Matches chip's IMBE full-rate path.
    SpecLcg,
    /// `(173, 13849, 65536)` LCG with seed `60584`. Matches the DVSI
    /// AMBE-3000R chip on the P25 half-rate path. Empirically derived
    /// by Probe 1 (correlation 0.945 vs chip on an all-UV probe stream).
    ChipLcg,
}

impl UnvoicedNoiseGen {
    /// LCG params: (multiplier, addend, modulus, initial state).
    pub const fn params(self) -> (u32, u32, u32, u32) {
        match self {
            UnvoicedNoiseGen::SpecLcg => (171, 11213, 53125, 3147),
            UnvoicedNoiseGen::ChipLcg => (173, 13849, 65536, 60584),
        }
    }
}

/// Default LCG generator when [`UnvoicedSynthState::new`] is called
/// without an explicit selection. Matches the DVSI chip's half-rate
/// behaviour (Probe 1 finding). Override via env var `BLIP25_LCG=spec`
/// for the BABA-A §1.12.1 spec-literal LCG.
///
/// Per-rate paths construct their state via [`UnvoicedSynthState::with_gen`]
/// to pick the right LCG explicitly:
///   - IMBE full-rate (Phase 1) → [`UnvoicedNoiseGen::SpecLcg`]
///   - AMBE+2 half-rate (Phase 2) → [`UnvoicedNoiseGen::ChipLcg`]
///
/// The env override only affects bare `new()` calls; rate-aware paths
/// continue to dispatch explicitly so existing test behaviour is
/// preserved.
fn default_lcg_gen() -> UnvoicedNoiseGen {
    match std::env::var("BLIP25_LCG").ok().as_deref() {
        Some("spec") | Some("Spec") => UnvoicedNoiseGen::SpecLcg,
        _ => UnvoicedNoiseGen::ChipLcg,
    }
}

/// Unvoiced-synthesis spectral scale `γ_w` per BABA-A §1.12.1 Eq. 121:
///
/// ```text
/// γ_w = [Σ_{n=−110}^{110} wR(n)]
///       · sqrt( Σ_{n=−104}^{104} wS²(n) / Σ_{n=−110}^{110} wR²(n) )
/// ```
///
/// Pre-evaluated from the committed Annex I (synthesis window) and
/// Annex C (pitch refinement window) CSVs. Decoder-init constant.
pub const GAMMA_W: f64 = 146.643269;


/// Cross-frame state for the unvoiced synthesizer (§1.12.1).
#[derive(Clone, Debug)]
pub struct UnvoicedSynthState {
    /// LCG multiplier `a`, addend `c`, modulus `m` and live state. The
    /// generator runs `lcg ← (a·lcg + c) mod m`. Default (spec LCG) is
    /// `(171, 11213, 53125)` with seed `3147`; the chip LCG used by
    /// the AMBE-3000R half-rate path is `(173, 13849, 65536)` with
    /// seed `60584` — selected via [`UnvoicedSynthState::with_gen`].
    lcg_a: u32,
    lcg_c: u32,
    lcg_m: u32,
    lcg: u32,
    /// Most recent 209 noise samples covering the synthesis window
    /// `wS` over `n = −104..104`.
    noise_window: [f64; 209],
    /// Previous frame's IDFT output `ũ_w(n, −1)` for `n = −128..127`.
    /// Indexed `[n + 128]`.
    prev_idft: [f64; 256],
    /// `false` until the first frame has populated `noise_window`.
    initialized: bool,
}

impl UnvoicedSynthState {
    /// Cold-start state per §1.13 / Annex A: empty noise window, LCG
    /// seeded per `default_lcg_gen` (chip generator by default, spec
    /// when `BLIP25_LCG=spec`), prev IDFT all zero.
    pub fn new() -> Self {
        Self::with_gen(default_lcg_gen())
    }

    /// Construct with an explicit LCG generator. Rate-aware code should
    /// always use this (e.g. IMBE full-rate → [`UnvoicedNoiseGen::SpecLcg`],
    /// AMBE+2 half-rate → [`UnvoicedNoiseGen::ChipLcg`]) so the env
    /// override doesn't accidentally cross-contaminate generations.
    pub fn with_gen(gen: UnvoicedNoiseGen) -> Self {
        let (a, c, m, seed) = gen.params();
        Self {
            lcg_a: a,
            lcg_c: c,
            lcg_m: m,
            lcg: seed,
            noise_window: [0.0; 209],
            prev_idft: [0.0; 256],
            initialized: false,
        }
    }

    /// Advance the LCG one step using the per-instance params:
    /// `lcg ← (a·lcg + c) mod m`. Defaults follow [`UnvoicedNoiseGen`].
    #[inline]
    fn next_noise(&mut self) -> f64 {
        self.lcg = self.lcg_a
            .wrapping_mul(self.lcg)
            .wrapping_add(self.lcg_c)
            % self.lcg_m;
        f64::from(self.lcg)
    }

    /// Per-frame: shift the noise window forward by N=160 samples.
    ///
    /// Frame-0 init scheme is controlled by env `BLIP25_LCG_INIT`:
    ///   - `burn` (default): fill all 209 positions from LCG up front.
    ///   - `zero`: start with all-zero buffer (JMBE-style; first frame
    ///     uses no LCG samples on the synth side).
    ///   - `burn_160` / `burn_96` / `burn_64` / `burn_49`: partial burn
    ///     for chip-init alignment probes.
    /// Subsequent frames always reuse the trailing 49 samples and
    /// generate 160 new ones (this matches the geometry of our
    /// 209-sample window).
    fn advance_window(&mut self) {
        if !self.initialized {
            let burn = match std::env::var("BLIP25_LCG_INIT").ok().as_deref() {
                Some("zero") | Some("Zero") => 0,
                Some("burn_64") => 64,
                Some("burn_96") => 96,
                Some("burn_128") => 128,
                Some("burn_160") => 160,
                Some("burn_49") => 49,
                _ => 209,
            };
            // Fill the latest `burn` positions of the window from LCG;
            // earlier positions stay at their default 0.0. When burn=209
            // (default) the whole window is fresh LCG; when burn=0 the
            // whole window is zero (JMBE-style); intermediate values
            // emulate a partial chip-side init.
            let start = 209usize.saturating_sub(burn);
            for i in start..209 {
                self.noise_window[i] = self.next_noise();
            }
            self.initialized = true;
        } else {
            self.noise_window.copy_within(160..209, 0);
            for i in 49..209 {
                self.noise_window[i] = self.next_noise();
            }
        }
    }
}

impl Default for UnvoicedSynthState {
    fn default() -> Self {
        Self::new()
    }
}

/// 256-point complex DFT of a windowed real input over `n = −104..=104`.
/// Returns `(re, im)` arrays of length 256 indexed `[m + 128]` for
/// `m = −128..127`.
///
/// Backed by [`rustfft`] — `O(N log N)` per frame. Input is mapped
/// from logical n ∈ [−104, 104] into DFT-natural index n_nat = n mod
/// 256; output is mapped from DFT-natural k into the centered m_idx
/// = m + 128 layout via the same 128-entry rotate. FFT is general
/// DSP (not P25 IP), so the swap is allowed under the clean-room rule
/// — same justification as the encode-side `signal_spectrum`.
fn dft_256_windowed(input: &[f64; 209]) -> ([f64; 256], [f64; 256]) {
    use num_complex::Complex;
    use rustfft::{Fft, FftPlanner};
    use std::sync::OnceLock;

    static FFT: OnceLock<std::sync::Arc<dyn Fft<f64>>> = OnceLock::new();
    let fft = FFT.get_or_init(|| {
        let mut planner = FftPlanner::<f64>::new();
        planner.plan_fft_forward(256)
    });

    // input[i] = signal at logical n = i - 104. Pack into natural
    // index n_nat = n.rem_euclid(256). For n ∈ [0, 104], n_nat = n;
    // for n ∈ [−104, −1], n_nat = n + 256 ∈ [152, 255].
    let mut buf = [Complex::<f64>::new(0.0, 0.0); 256];
    for (i, &x) in input.iter().enumerate() {
        let n = i as i32 - 104;
        let n_nat = n.rem_euclid(256) as usize;
        buf[n_nat].re = x;
    }
    fft.process(&mut buf);

    // Output rotation: re[m_idx] = buf[(m_idx + 128) % 256] places
    // m = 0 at m_idx = 128 and m = −128 at m_idx = 0.
    let mut re = [0f64; 256];
    let mut im = [0f64; 256];
    for m_idx in 0..256 {
        let k = (m_idx + 128) & 255;
        re[m_idx] = buf[k].re;
        im[m_idx] = buf[k].im;
    }
    (re, im)
}

/// 256-point inverse DFT producing real-valued output (Eq. 125).
/// The 1/256 normalization is included.
///
/// Backed by [`rustfft`]. Input is the centered (re, im) packing
/// produced by [`dft_256_windowed`] / [`shape_spectrum`]; output is
/// indexed by `n_idx = n + 128` for `n ∈ [−128, 127]`. The IDFT of
/// a hermitian-symmetric spectrum is real; we take the real part to
/// drop the round-off-only imaginary residue.
fn idft_256(re: &[f64; 256], im: &[f64; 256]) -> [f64; 256] {
    use num_complex::Complex;
    use rustfft::{Fft, FftPlanner};
    use std::sync::OnceLock;

    static IFFT: OnceLock<std::sync::Arc<dyn Fft<f64>>> = OnceLock::new();
    let fft = IFFT.get_or_init(|| {
        let mut planner = FftPlanner::<f64>::new();
        planner.plan_fft_inverse(256)
    });

    // Undo the centered (m + 128) rotation: freq[k] = packed[(k + 128) mod 256].
    let mut buf = [Complex::<f64>::new(0.0, 0.0); 256];
    for k in 0..256 {
        let m_idx = (k + 128) & 255;
        buf[k] = Complex::new(re[m_idx], im[m_idx]);
    }
    fft.process(&mut buf);

    // rustfft's inverse is unnormalized — divide by N. Output index
    // mirrors the input rotation: out[n_idx] = ifft_buf[(n_idx + 128) % 256].
    let mut out = [0f64; 256];
    for n_idx in 0..256 {
        let k = (n_idx + 128) & 255;
        out[n_idx] = buf[k].re / 256.0;
    }
    out
}

/// Apply the per-band spectral shaping (Eq. 119–124) in-place on the
/// DFT outputs.
fn shape_spectrum(
    re: &mut [f64; 256],
    im: &mut [f64; 256],
    omega_0: f64,
    m_bar: &[f32],
    v_bar: &[bool],
    gamma_w: f64,
) {
    let l_count = m_bar.len();
    let scale = 256.0 / (2.0 * PI64);

    // First, zero everything outside band 1..L̃ (Eq. 124). Band edges:
    //   ⌈ã_1⌉ at the low end, ⌈b̃_L̃⌉ at the high end.
    let a1 = (scale * 0.5 * omega_0).ceil() as i32;
    let b_last = (scale * (l_count as f64 + 0.5) * omega_0).ceil() as i32;
    for m_idx in 0..256 {
        let m = m_idx as i32 - 128;
        if m.unsigned_abs() < a1 as u32 || (m.unsigned_abs() as i32) >= b_last {
            re[m_idx] = 0.0;
            im[m_idx] = 0.0;
        }
    }

    // Pre-compute the unmodified spectrum power for each band's norm
    // sum, since we'll be overwriting `re`/`im` during the sweep.
    let mut band_norm: Vec<f64> = Vec::with_capacity(l_count);
    let mut band_edges: Vec<(i32, i32)> = Vec::with_capacity(l_count);
    for l in 1..=l_count as i32 {
        let l_f = f64::from(l);
        let a_l = (scale * (l_f - 0.5) * omega_0).ceil() as i32;
        let b_l = (scale * (l_f + 0.5) * omega_0).ceil() as i32;
        band_edges.push((a_l, b_l));
        // Norm sum: η in [⌈ã_l⌉, ⌈b̃_l⌉) (half-open, count = b - a).
        let mut norm_sum = 0f64;
        let count = (b_l - a_l).max(0) as usize;
        for eta in a_l..b_l {
            // Both +η and −η contribute (the spec's "|m|" means both signs).
            for &sign in &[1i32, -1] {
                let m_idx = (sign * eta + 128) as usize;
                if m_idx < 256 {
                    norm_sum += re[m_idx] * re[m_idx] + im[m_idx] * im[m_idx];
                }
            }
        }
        // Average over (count) magnitudes counted twice (positive + negative).
        let norm = if count > 0 && norm_sum > 0.0 {
            (norm_sum / (2.0 * count as f64)).sqrt()
        } else {
            1.0
        };
        band_norm.push(norm);
    }

    // Per-band assignment. Half-open `[a_l, b_l)` so adjacent bands
    // tile cleanly: `b_l == a_{l+1}` by construction (Eq. 122/123),
    // so the bin at `m = b_l` belongs to band l+1, not band l. Using
    // an inclusive range `[a_l, b_l]` here (an earlier port) caused
    // that shared bin to (a) get scaled by the latter band's factor
    // but (b) not contribute to the *current* band's norm sum, which
    // mis-balanced unvoiced energy on dense-band frames (audible on
    // fricatives like the /ks/ in "axe"). JMBE / the DVSI chip both
    // use `[a, b)` consistently across norm and application.
    for (l_idx, &(a_l, b_l)) in band_edges.iter().enumerate() {
        let voiced = v_bar[l_idx];
        if voiced {
            // Eq. 119: zero voiced bands.
            for m_abs in a_l..b_l {
                for &sign in &[1i32, -1] {
                    let m = sign * m_abs;
                    let m_idx = (m + 128) as usize;
                    if m_idx < 256 {
                        re[m_idx] = 0.0;
                        im[m_idx] = 0.0;
                    }
                }
            }
        } else {
            // Eq. 120: scale by γ_w · M̄_l / norm.
            let m_bar_l = f64::from(m_bar[l_idx]);
            let norm = band_norm[l_idx];
            let factor = if norm > 0.0 {
                gamma_w * m_bar_l / norm
            } else {
                0.0
            };
            for m_abs in a_l..b_l {
                for &sign in &[1i32, -1] {
                    let m = sign * m_abs;
                    let m_idx = (m + 128) as usize;
                    if m_idx < 256 {
                        re[m_idx] *= factor;
                        im[m_idx] *= factor;
                    }
                }
            }
        }
    }
}

/// Synthesize the unvoiced component for one 20 ms frame per BABA-A
/// §1.12.1 Eq. 117–126.
///
/// Returns 160 PCM-domain samples (in `f64`; the final cast to `i16`
/// happens after voiced + unvoiced sum at the synthesis top level).
///
/// `gamma_w` is the spectral scale from Eq. 121. Pass
/// [`GAMMA_W`] until Annex C lands or a DVSI fixture
/// value is available.
pub fn synthesize_unvoiced(
    omega_0: f32,
    m_bar: &[f32],
    v_bar: &[bool],
    gamma_w: f64,
    state: &mut UnvoicedSynthState,
) -> [f64; FRAME_SAMPLES] {
    state.advance_window();

    // Window the noise: u(n) · wS(n) for n = −104..104.
    let mut windowed = [0f64; 209];
    for i in 0..209 {
        let n = i as i32 - 104;
        windowed[i] = state.noise_window[i] * f64::from(synth_window(n));
    }

    let (mut re, mut im) = dft_256_windowed(&windowed);
    shape_spectrum(&mut re, &mut im, f64::from(omega_0), m_bar, v_bar, gamma_w);
    let u_w = idft_256(&re, &im);

    // Eq. 126: weighted overlap-add with previous frame.
    let mut out = [0f64; FRAME_SAMPLES];
    for n in 0..FRAME_SAMPLES as i32 {
        let ws_n = f64::from(synth_window(n));
        let ws_n_minus_n = f64::from(synth_window(n - FRAME_SAMPLES as i32));
        let prev_term = if (0..=127).contains(&n) {
            ws_n * state.prev_idft[(n + 128) as usize]
        } else {
            0.0
        };
        let curr_term = if (32..=159).contains(&n) {
            ws_n_minus_n * u_w[(n - 32) as usize]
        } else {
            0.0
        };
        let denom = ws_n * ws_n + ws_n_minus_n * ws_n_minus_n;
        out[n as usize] = if denom > 1e-30 {
            (prev_term + curr_term) / denom
        } else {
            0.0
        };
    }

    state.prev_idft = u_w;
    out
}

/// Snapshot the noise samples that the voiced synthesizer needs for
/// phase randomization (`u(l)` for `l = 1..=56`, per Eq. 141).
///
/// Reads from the *current* (post-`advance_window`) noise window, so
/// the caller must invoke this after [`synthesize_unvoiced`] — at
/// which point both synthesizers see the same shifted sequence per
/// the §1.13 cross-frame state contract.
pub fn voiced_noise_samples(state: &UnvoicedSynthState) -> [f64; L_MAX as usize] {
    let mut out = [0f64; L_MAX as usize];
    for l in 1..=L_MAX as usize {
        // u(l) at global n = l, indexed in the window at l + 104.
        out[l - 1] = state.noise_window[l + 104];
    }
    out
}

// ---------------------------------------------------------------------------
// §1.12.2 — Voiced synthesis (Eq. 127–141)
// ---------------------------------------------------------------------------

/// Initial value for `ω̃₀(−1)` per BABA-A §10 Annex A.
pub const INIT_PREV_OMEGA_0: f64 = 0.02985 * PI64;

/// Cross-frame state for the voiced synthesizer (§1.12.2).
#[derive(Clone, Debug)]
pub struct VoicedSynthState {
    /// `φ_l(−1)` for `l = 1..=56`; index 0 unused. Init: 0.
    phi: [f64; L_MAX as usize + 1],
    /// `ψ_l(−1)` (auxiliary phase) for `l = 1..=56`; index 0 unused. Init: 0.
    psi: [f64; L_MAX as usize + 1],
    /// Previous frame's enhanced amplitudes `M̄_l(−1)`; index 0 unused. Init: 0.
    prev_m_bar: [f64; L_MAX as usize + 1],
    /// Previous frame's per-harmonic V/UV `v̄_l(−1)`; index 0 unused. Init: false.
    prev_v_bar: [bool; L_MAX as usize + 1],
    /// Previous frame's harmonic count `L̃(−1)`. Init: 30.
    prev_l: u8,
    /// Previous frame's fundamental frequency `ω̃₀(−1)`. Init: 0.02985·π.
    prev_omega_0: f64,
}

impl VoicedSynthState {
    /// Cold-start state per §1.13 / Annex A.
    pub fn new() -> Self {
        Self {
            phi: [0.0; L_MAX as usize + 1],
            psi: [0.0; L_MAX as usize + 1],
            prev_m_bar: [0.0; L_MAX as usize + 1],
            prev_v_bar: [false; L_MAX as usize + 1],
            prev_l: 30,
            prev_omega_0: INIT_PREV_OMEGA_0,
        }
    }
}

impl Default for VoicedSynthState {
    fn default() -> Self {
        Self::new()
    }
}

/// Wrap a phase value to the principal range `[−π, π)` per Eq. 138's
/// `Δφ − 2π·⌊(Δφ+π)/(2π)⌋` formulation.
#[inline]
fn wrap_phase(delta_phi: f64) -> f64 {
    delta_phi - 2.0 * PI64 * ((delta_phi + PI64) / (2.0 * PI64)).floor()
}

/// Phase computation mode for the voiced synthesizer.
///
/// - [`PhaseMode::Baseline`] follows BABA-A §1.12.2 Eq. 141 — noise-
///   perturbed phase driven by the LCG u(n) sequence. Used by full-rate
///   IMBE (P25 Phase 1).
/// - [`PhaseMode::AmbePlus`] replaces Eq. 141 with US5701390 phase
///   regeneration (Hilbert-transform of log-magnitude). Used by
///   half-rate AMBE+2 (P25 Phase 2). See
///   `~/blip25-specs/DVSI/AMBE-3000/AMBE-3000_Decoder_Implementation_Spec.md`
///   §5 and [`phase_regen`].
#[derive(Copy, Clone, Debug, Eq, PartialEq)]
pub enum PhaseMode {
    /// BABA-A Eq. 141 noise-perturbed phase. Full-rate IMBE default.
    Baseline,
    /// US5701390 §5 Hilbert-transform phase regeneration.
    /// AMBE+2 / P25 half-rate.
    AmbePlus,
}

/// Synthesize the voiced component for one 20 ms frame per BABA-A
/// §1.12.2 Eq. 127–141.
///
/// Returns 160 PCM-domain samples (in `f64`; the final `i16` cast
/// happens at the synthesis top level after summing with the unvoiced
/// component per Eq. 142).
///
/// `noise_samples` carries `u(l)` for `l = 1..=56` (index 0 = `u(1)`),
/// extracted from [`voiced_noise_samples`] after the unvoiced synth's
/// noise advance. Consumed only when `phase_mode == PhaseMode::Baseline`.
///
/// `phase_mode` selects between BABA-A Eq. 141 (baseline IMBE) and
/// US5701390 phase regeneration (AMBE+2 half-rate). See [`PhaseMode`].
pub fn synthesize_voiced(
    omega_0: f32,
    m_bar: &[f32],
    v_bar: &[bool],
    noise_samples: &[f64; L_MAX as usize],
    phase_mode: PhaseMode,
    state: &mut VoicedSynthState,
) -> [f64; FRAME_SAMPLES] {
    // Default to the spec LCG modulus for ρ_l (Eq. 141) — used by all
    // existing tests and by full-rate IMBE. AmbePlus phase mode ignores
    // the modulus (it bypasses Eq. 141 in favour of phase regen).
    synthesize_voiced_with_lcg_modulus(omega_0, m_bar, v_bar, noise_samples, phase_mode, 53125, state)
}

/// Variant of [`synthesize_voiced`] that lets the caller specify the
/// LCG modulus used in Eq. 141's `ρ_l = (2π/m)·u(l) − π`. Set this to
/// the same modulus used by the unvoiced [`UnvoicedSynthState`]'s
/// generator so ρ_l stays uniform in `[−π, π)` regardless of which
/// LCG generates the noise samples.
pub fn synthesize_voiced_with_lcg_modulus(
    omega_0: f32,
    m_bar: &[f32],
    v_bar: &[bool],
    noise_samples: &[f64; L_MAX as usize],
    phase_mode: PhaseMode,
    state_lcg_m: u32,
    state: &mut VoicedSynthState,
) -> [f64; FRAME_SAMPLES] {
    let l_curr = m_bar.len() as u8;
    let l_prev = state.prev_l;
    let max_l = l_curr.max(l_prev) as usize;

    let omega_curr = f64::from(omega_0);
    let omega_prev = state.prev_omega_0;
    let n_f = FRAME_SAMPLES as f64;
    let delta_omega = omega_curr - omega_prev;

    // ψ_l(0) = ψ_l(−1) + (ω̃₀(−1) + ω̃₀(0))·l·N/2  (Eq. 139, all l in 1..=56)
    let mut psi_curr = [0f64; L_MAX as usize + 1];
    for l in 1..=L_MAX as usize {
        let l_f = l as f64;
        psi_curr[l] = state.psi[l] + (omega_prev + omega_curr) * l_f * n_f / 2.0;
    }

    // φ_l(0) — Eq. 140 with branch dispatched by phase_mode.
    // Branch 1 for 1 ≤ l ≤ ⌊L̃/4⌋: φ = ψ  (both modes agree).
    // Branch 2 for ⌊L̃/4⌋ < l ≤ max(L̃(−1), L̃(0)):
    //   Baseline (Eq. 141): φ = ψ + L̃_uv · ρ_l / L̃ · l
    //     where ρ_l = (2π/53125)·u(l) − π.
    //   AmbePlus (US5701390 §5.3): φ = ψ + φ_regen,l
    //     where φ_regen,l = γ · Σ h(m)·log₂(M̄_{l+m}).
    // For l beyond max (up to 56): apply branch 2 too — these harmonics
    // are out-of-range (treated as unvoiced per Eq. 128–129) and their
    // phase still evolves for next-frame continuity.
    let l_quarter = (l_curr as usize) / 4;
    let mut phi_curr = [0f64; L_MAX as usize + 1];
    match phase_mode {
        PhaseMode::Baseline => {
            let l_uv: u8 = (0..l_curr as usize).filter(|&i| !v_bar[i]).count() as u8;
            // The 53125 in Eq. 141's ρ_l = (2π/53125)·u(l) − π is the
            // SPEC LCG modulus. For the chip LCG (modulus 65536), the
            // formula should use 65536 instead so ρ_l stays uniform
            // in [-π, π). Map the divisor to the live LCG modulus.
            let lcg_m = f64::from(state_lcg_m);
            for l in 1..=L_MAX as usize {
                if l <= l_quarter {
                    phi_curr[l] = psi_curr[l];
                } else {
                    let u_l = noise_samples[l - 1];
                    let rho_l = (2.0 * PI64 / lcg_m) * u_l - PI64; // Eq. 141
                    let l_curr_f = if l_curr > 0 { f64::from(l_curr) } else { 1.0 };
                    phi_curr[l] = psi_curr[l]
                        + f64::from(l_uv) * rho_l / l_curr_f * (l as f64);
                }
            }
        }
        PhaseMode::AmbePlus => {
            // Compute φ_regen,l for l = 1..=L̃ from the current frame's
            // M̄_l. noise_samples is ignored in this branch.
            let mut phi_regen = [0f64; L_MAX as usize + 1];
            phase_regen::ambe_phase_regen(m_bar, &mut phi_regen);
            for l in 1..=L_MAX as usize {
                if l <= l_quarter {
                    phi_curr[l] = psi_curr[l];
                } else {
                    // Harmonic indices past L̃(0) get φ_regen = 0 from
                    // the untouched-slot convention in ambe_phase_regen,
                    // so φ = ψ there (consistent with the spec's
                    // zero-padding of B_l beyond L̃).
                    phi_curr[l] = psi_curr[l] + phi_regen[l];
                }
            }
        }
    }

    // Pre-cache wS(n) and wS(n−N) over n = 0..N to avoid 56·160·2
    // bounds-checked lookups in the inner loop.
    let mut ws_n = [0f64; FRAME_SAMPLES];
    let mut ws_nm = [0f64; FRAME_SAMPLES];
    for n in 0..FRAME_SAMPLES {
        ws_n[n] = f64::from(synth_window(n as i32));
        ws_nm[n] = f64::from(synth_window(n as i32 - FRAME_SAMPLES as i32));
    }

    // Per-harmonic synthesis + accumulate.
    //
    // The three constant-frequency branches (VUv, UvV, VVSum) each
    // need cos(ω·l·n + φ) sampled over n = 0..N. We replace those
    // with a complex-phasor recursion: starting from z₀ = exp(jφ₀)
    // and stepping by w = exp(jω·l) per sample, z_{n+1} = z_n · w
    // and Re(z_n) gives the cosine. This costs 2 cos + 2 sin per
    // harmonic (init + step rotation) plus 4 mults + 2 adds per
    // sample, vs. one transcendental per sample in the direct form
    // — a ~30× math reduction in the inner loop. Magnitude drift
    // over 160 steps is ≪ 1e−13, well below audible.
    //
    // VVRamp's quadratic phase term (Eq. 134) doesn't reduce to a
    // constant rotation; keep it on `cos()`. It only fires for
    // low-l harmonics with small |Δω·l/ω| and is the rarer branch.
    enum Branch {
        UvUv,
        VUv,
        UvV,
        VVSum,
        VVRamp,
    }
    let mut s_v = [0f64; FRAME_SAMPLES];
    for l in 1..=max_l {
        let l_f = l as f64;
        let m_curr_l = if l <= l_curr as usize { f64::from(m_bar[l - 1]) } else { 0.0 };
        let m_prev_l = if l <= l_prev as usize { state.prev_m_bar[l] } else { 0.0 };
        let v_curr = l <= l_curr as usize && v_bar[l - 1];
        let v_prev = l <= l_prev as usize && state.prev_v_bar[l];

        let branch = match (v_prev, v_curr) {
            (false, false) => Branch::UvUv,
            (true, false) => Branch::VUv,
            (false, true) => Branch::UvV,
            (true, true) => {
                let pitch_change_ratio = if omega_curr.abs() > 1e-30 {
                    (delta_omega * l_f / omega_curr).abs()
                } else {
                    0.0
                };
                if l >= 8 || pitch_change_ratio >= 0.1 {
                    Branch::VVSum
                } else {
                    Branch::VVRamp
                }
            }
        };

        match branch {
            Branch::UvUv => {}
            Branch::VUv => {
                // Phase is omega_prev * n * l + state.phi[l], n = 0..N.
                let phi0 = state.phi[l];
                let (mut pre, mut pim) = (phi0.cos(), phi0.sin());
                let step = omega_prev * l_f;
                let (dc, ds) = (step.cos(), step.sin());
                let scale = 2.0 * m_prev_l;
                for n in 0..FRAME_SAMPLES {
                    s_v[n] += scale * ws_n[n] * pre;
                    let npre = pre * dc - pim * ds;
                    let npim = pre * ds + pim * dc;
                    pre = npre;
                    pim = npim;
                }
            }
            Branch::UvV => {
                // Phase is omega_curr * (n - N) * l + phi_curr[l].
                // Start phase at n = 0: omega_curr * (-N) * l + phi_curr[l].
                let phi0 = phi_curr[l] - omega_curr * n_f * l_f;
                let (mut cre, mut cim) = (phi0.cos(), phi0.sin());
                let step = omega_curr * l_f;
                let (dc, ds) = (step.cos(), step.sin());
                let scale = 2.0 * m_curr_l;
                for n in 0..FRAME_SAMPLES {
                    s_v[n] += scale * ws_nm[n] * cre;
                    let ncre = cre * dc - cim * ds;
                    let ncim = cre * ds + cim * dc;
                    cre = ncre;
                    cim = ncim;
                }
            }
            Branch::VVSum => {
                // Two phasors: prev frame (omega_prev) and current
                // frame (omega_curr) overlap-added with the cross
                // window pair.
                let (mut pre, mut pim) = (state.phi[l].cos(), state.phi[l].sin());
                let pstep = omega_prev * l_f;
                let (pdc, pds) = (pstep.cos(), pstep.sin());
                let phi_c0 = phi_curr[l] - omega_curr * n_f * l_f;
                let (mut cre, mut cim) = (phi_c0.cos(), phi_c0.sin());
                let cstep = omega_curr * l_f;
                let (cdc, cds) = (cstep.cos(), cstep.sin());
                let scale_p = 2.0 * m_prev_l;
                let scale_c = 2.0 * m_curr_l;
                for n in 0..FRAME_SAMPLES {
                    s_v[n] += scale_p * ws_n[n] * pre + scale_c * ws_nm[n] * cre;
                    let npre = pre * pdc - pim * pds;
                    let npim = pre * pds + pim * pdc;
                    pre = npre;
                    pim = npim;
                    let ncre = cre * cdc - cim * cds;
                    let ncim = cre * cds + cim * cdc;
                    cre = ncre;
                    cim = ncim;
                }
            }
            Branch::VVRamp => {
                // Quadratic phase — keep direct cos().
                let delta_phi = phi_curr[l] - state.phi[l]
                    - (omega_prev + omega_curr) * l_f * n_f / 2.0;
                let wrapped = wrap_phase(delta_phi);
                let delta_omega_l = wrapped / n_f;
                let theta0 = state.phi[l];
                let coef_n2 = (omega_curr * l_f - omega_prev * l_f) / (2.0 * n_f);
                for n in 0..FRAME_SAMPLES {
                    let nf = n as f64;
                    let a_l = m_prev_l + (nf / n_f) * (m_curr_l - m_prev_l);
                    let theta_l = theta0
                        + (omega_prev * l_f + delta_omega_l) * nf
                        + coef_n2 * nf * nf;
                    s_v[n] += 2.0 * a_l * theta_l.cos();
                }
            }
        }
    }

    // Update state for next frame.
    for l in 1..=L_MAX as usize {
        state.psi[l] = psi_curr[l];
        state.phi[l] = phi_curr[l];
    }
    for l in 1..=l_curr as usize {
        state.prev_m_bar[l] = f64::from(m_bar[l - 1]);
        state.prev_v_bar[l] = v_bar[l - 1];
    }
    // Zero out fields beyond the current L̃ so they don't leak.
    for l in (l_curr as usize + 1)..=L_MAX as usize {
        state.prev_m_bar[l] = 0.0;
        state.prev_v_bar[l] = false;
    }
    state.prev_l = l_curr;
    state.prev_omega_0 = omega_curr;

    s_v
}

// ---------------------------------------------------------------------------
// Top-level synthesis (Eq. 142)
// ---------------------------------------------------------------------------

/// Bundled cross-frame state for the full §1.10–§1.12 synthesis
/// pipeline. Owns the unvoiced and voiced sub-states plus the scalar
/// state from §1.10 / §1.11 (S_E, τ_M, ε_R). Also carries the
/// per-frame [`FrameErrorContext`] and calibration scale `γ_w` so the
/// codec-generation entry points can hide those knobs from consumers
/// that don't need to tune them.
#[derive(Clone, Debug)]
pub struct SynthState {
    /// §1.10 local-energy state.
    pub s_e: f64,
    /// §1.11 amplitude-smoothing state.
    pub tau_m: f64,
    /// §1.11 smoothed error rate.
    pub epsilon_r: f64,
    /// §1.12.1 unvoiced sub-state.
    pub unvoiced: UnvoicedSynthState,
    /// §1.12.2 voiced sub-state.
    pub voiced: VoicedSynthState,
    /// Per-frame error context read by the codec-generation entry
    /// points ([`crate::codecs::ambe::synthesize_frame`] and siblings)
    /// on each call. Defaults to "error-free" on
    /// [`SynthState::new`]; consumers with live FEC error counts
    /// (e.g. from a [`crate::imbe_wire::frame::Frame::errors`]
    /// total) should update it before each synth call.
    pub err: FrameErrorContext,
    /// §1.12.1 spectral calibration scale. Defaults to [`GAMMA_W`];
    /// override for calibration studies.
    pub gamma_w: f64,
    /// Snapshot of the last successfully-synthesized frame's parameters
    /// (for Eq. 99–104 substitution on Repeat/Mute). `None` until the
    /// first valid frame.
    last_good: Option<LastGoodFrame>,
    /// Count of consecutive Repeat/Mute frames since the last `Use`.
    /// Drives the optional `repeat_reset_after` heuristic.
    repeat_count: u32,
    /// Beyond-spec consecutive-repeat reset threshold. When `Some(n)`,
    /// after `n` consecutive Repeat/Mute frames the substituted
    /// parameters fall back to a default-fundamental + amps=1.0 frame
    /// instead of the prior `last_good` snapshot. `None` (default) is
    /// the spec-faithful path: literal Eq. 99–104 indefinitely (gap
    /// 0022 resolution — PDF §7.7 prescribes plain assignment with no
    /// per-repeat attenuation; this knob exists only for chip-stream
    /// interop quality, not for spec conformance).
    repeat_reset_after: Option<u32>,
    /// Beyond-spec error-rate freeze on Repeat. When `true`, a frame
    /// whose disposition is `Repeat` does not advance `ε_R`; the
    /// previous frame's `ε_R` is carried into the next frame
    /// unchanged. Mirrors JMBE's `IMBEModelParameters.copy()` calling
    /// `setErrorRate(previous.getErrorRate())` (gap 0021 resolution).
    /// Default `false` keeps the spec-faithful path that always
    /// smooths `ε_R = 0.95·prev + 0.05·(ε_T/144)` regardless of
    /// disposition. Enable for chip-encoded-bits decode quality (the
    /// chip emits frames with deterministic 13-error patterns that
    /// would otherwise drive `ε_R` past the 0.0875 Mute threshold and
    /// silence ~99% of frames).
    chip_compat: bool,
    /// Beyond-spec spectral-discontinuity clamp (gap 0026). When
    /// `true`, frames whose harmonic count `L` jumps by more than
    /// [`SPECTRAL_CLAMP_L_THRESHOLD`] from the previous frame have
    /// their post-smoothing amplitudes scaled by
    /// [`SPECTRAL_CLAMP_GAMMA`] for one frame. Mirrors observed
    /// DVSI AMBE-3000R behavior on pitch/L jumps (chip RMS ~73 % of
    /// ours at the jump frame, transient single-frame). BABA-A
    /// §1.11.3 does not describe this clamp; gap 0026 documents the
    /// probe data. Default `false` keeps the spec-faithful path.
    ///
    /// **Tracks a sweep characterisation, not a chip model.** Gap
    /// 0029 falsified the per-bin amplitude lowpass mechanism at
    /// the M̄ level — the chip ≈ spec-faithful `√(Σ M̄²)` energy
    /// on average, with a structured ~8% residual concentrated in
    /// the (low_b̂₀ ∈ [10, 40], high_b̂₀ ∈ [60, 119]) corner of the
    /// (b̂₀_lo, b̂₀_hi) plane. The |ΔL|-keyed heuristic implemented
    /// here is a correlation, not a mechanism; spec-author
    /// recommends strict-spec (`γ_M = 1.0`) as the default and
    /// flags this knob as "likely to be revised pending the 4%
    /// PCM-vs-`spec_atten` offset investigation on the implementer
    /// side." Keep the default `false`; only enable for downstream
    /// demos that already opted into chip-correlation behaviour.
    chip_compat_spectral_clamp: bool,
    /// Previous frame's harmonic count `L`, for
    /// [`Self::chip_compat_spectral_clamp`] discontinuity detection.
    /// 0 means no prior frame yet (no clamp on first frame).
    prev_l: u8,
    /// The disposition (`Use` / `Repeat` / `Mute`) selected by the
    /// most recent [`synthesize_frame`] / [`synthesize_frame_ambe_plus`]
    /// call. `None` until at least one frame has been synthesized.
    /// Updated **inside** the synth before any state advancement so
    /// callers reading it after the synth returns get the current
    /// frame's decision, not the prior one.
    last_disposition: Option<FrameDisposition>,
}

/// Minimum |ΔL| between consecutive frames that triggers the beyond-spec
/// spectral-discontinuity clamp ([`SynthState::chip_compat_spectral_clamp`]).
/// Calibrated from the 2-D b̂₀ sweep (2026-05-14, 25 probes spanning
/// ΔL ∈ [−37, 40] on a steady-state-with-single-jump-frame schedule):
///
/// | \|ΔL\| range | chip/ours ratio at jump frame |
/// |--------------|-------------------------------|
/// |   1 – 5      | 0.99 – 1.04 (within probe noise, no clamp) |
/// |   6 – 25     | 0.94 – 1.17 (mostly chip *louder* — overlap-add redistribution) |
/// |  26 – 40     | 0.74 – 0.82 (chip clamps consistently) |
/// |  35 – 37 down| 0.81 – 0.83 (clamp fires on big down-jumps too) |
///
/// So the clamp fires symmetrically for `|ΔL| > 25`. A previous, narrower
/// `> 5` threshold over-fired on natural speech transitions where the chip
/// does not clamp.
pub const SPECTRAL_CLAMP_L_THRESHOLD: i16 = 25;

/// Per-frame amplitude scale factor applied to `M̄_l` when the spectral-
/// discontinuity clamp fires. Mean chip/ours ratio across the
/// `|ΔL| > 25` regime in the 2-D sweep is ~0.78; on the noisy #3537 case
/// (FEC-error-driven trigger) the chip's RMS is ~0.72 of ours. We pick
/// the average so both regimes land within a few percent of chip behaviour.
pub const SPECTRAL_CLAMP_GAMMA: f64 = 0.75;

#[derive(Clone, Debug)]
struct LastGoodFrame {
    omega_0: f32,
    l: u8,
    voiced: [bool; L_MAX as usize],
    m_tilde: [f32; L_MAX as usize],
}

impl SynthState {
    /// Cold-start synthesis state per §1.13 / Annex A. Uses the
    /// default LCG generator (chip LCG unless `BLIP25_LCG=spec`).
    /// Rate-aware paths should prefer [`SynthState::with_unvoiced_gen`]
    /// to make the LCG choice explicit.
    pub fn new() -> Self {
        Self::with_unvoiced_gen(default_lcg_gen())
    }

    /// Cold-start with an explicit unvoiced LCG generator. Used by
    /// rate-aware codec wrappers:
    ///   - IMBE full-rate → [`UnvoicedNoiseGen::SpecLcg`]
    ///   - AMBE+2 half-rate → [`UnvoicedNoiseGen::ChipLcg`]
    pub fn with_unvoiced_gen(gen: UnvoicedNoiseGen) -> Self {
        Self {
            s_e: INIT_S_E,
            tau_m: INIT_TAU_M,
            epsilon_r: 0.0,
            unvoiced: UnvoicedSynthState::with_gen(gen),
            voiced: VoicedSynthState::new(),
            err: FrameErrorContext::default(),
            gamma_w: GAMMA_W,
            last_good: None,
            repeat_count: 0,
            repeat_reset_after: None,
            chip_compat: false,
            chip_compat_spectral_clamp: false,
            prev_l: 0,
            last_disposition: None,
        }
    }

    /// The disposition (`Use` / `Repeat` / `Mute`) selected by the
    /// most recent synth call. `None` until at least one frame has
    /// been synthesized; reset by [`Self::new`].
    #[inline]
    pub fn last_disposition(&self) -> Option<FrameDisposition> {
        self.last_disposition
    }

    /// Enable the beyond-spec consecutive-repeat reset heuristic. After
    /// `n` consecutive Repeat/Mute frames, fall back to a default
    /// fundamental + amps=1.0 frame instead of replaying the prior
    /// snapshot. `None` (default) keeps the literal Eq. 99–104 path —
    /// the spec-faithful behavior per gap 0022 resolution.
    ///
    /// JMBE's `IMBEModelParameters.copy()` uses `n = 3`; SDRTrunk
    /// matches. Use the same as a pragmatic chip-interop default when
    /// enabling the knob.
    pub fn set_repeat_reset_after(&mut self, n: Option<u32>) {
        self.repeat_reset_after = n;
    }

    /// Current consecutive-repeat reset threshold (`None` = disabled,
    /// spec-faithful).
    #[inline]
    pub fn repeat_reset_after(&self) -> Option<u32> {
        self.repeat_reset_after
    }

    /// Enable JMBE-style error-rate freeze on Repeat (beyond-spec).
    /// When enabled, frames whose disposition is `Repeat` do not
    /// advance `ε_R`; the previous frame's value is carried forward.
    /// This prevents runs of high-error chip-encoded frames from
    /// driving `ε_R` past the Mute threshold. `false` (default) is
    /// the spec-faithful path. See gap 0021.
    pub fn set_chip_compat(&mut self, on: bool) {
        self.chip_compat = on;
    }

    /// Current chip-compat (error-rate freeze on Repeat) setting.
    #[inline]
    pub fn chip_compat(&self) -> bool {
        self.chip_compat
    }

    /// Force the beyond-spec spectral-discontinuity clamp on
    /// independently of [`Self::set_chip_compat`]. Most consumers
    /// should leave this off and rely on the umbrella `chip_compat`
    /// flag, which auto-enables the clamp alongside the gap-0021
    /// ε_R freeze on Repeat. This standalone toggle is for the
    /// narrow case of wanting the clamp without the ε_R freeze.
    pub fn set_chip_compat_spectral_clamp(&mut self, on: bool) {
        self.chip_compat_spectral_clamp = on;
    }

    /// Current standalone spectral-discontinuity clamp setting.
    /// Note: the clamp also fires whenever
    /// [`Self::chip_compat`] is `true`.
    #[inline]
    pub fn chip_compat_spectral_clamp(&self) -> bool {
        self.chip_compat_spectral_clamp
    }
}

impl Default for SynthState {
    fn default() -> Self {
        Self::new()
    }
}

/// Convert one frame of `f64` synthesizer output into 16-bit signed
/// PCM with rounding and clamp at the i16 limits.
#[inline]
pub fn pcm_from_f64(samples: &[f64; FRAME_SAMPLES]) -> [i16; FRAME_SAMPLES] {
    let mut out = [0i16; FRAME_SAMPLES];
    for (i, &s) in samples.iter().enumerate() {
        let clamped = s.clamp(f64::from(i16::MIN), f64::from(i16::MAX));
        out[i] = clamped.round() as i16;
    }
    out
}

/// Synthesize one 20 ms frame (160 PCM samples) from `MbeParams` per
/// the §1.10–§1.12 pipeline, using baseline BABA-A Eq. 141 phase
/// (full-rate IMBE / P25 Phase 1).
///
/// `gamma_w` is the §1.12.1 spectral scale; pass [`GAMMA_W`]
/// until a calibrated value is available.
///
/// For AMBE+2 / half-rate P25 Phase 2 use
/// [`synthesize_frame_ambe_plus`] instead, which replaces Eq. 141
/// with US5701390 phase regeneration.
pub fn synthesize_frame(
    params: &MbeParams,
    err: &FrameErrorContext,
    gamma_w: f64,
    state: &mut SynthState,
) -> [i16; FRAME_SAMPLES] {
    synthesize_frame_with_mode(params, err, gamma_w, PhaseMode::Baseline, state)
}

/// AMBE+2 variant of [`synthesize_frame`] — applies US5701390 §5
/// phase regeneration in place of BABA-A Eq. 141. Used for P25
/// half-rate (Phase 2 TDMA) decoding.
pub fn synthesize_frame_ambe_plus(
    params: &MbeParams,
    err: &FrameErrorContext,
    gamma_w: f64,
    state: &mut SynthState,
) -> [i16; FRAME_SAMPLES] {
    synthesize_frame_with_mode(params, err, gamma_w, PhaseMode::AmbePlus, state)
}

/// Synthesize a repeated frame for an erasure, using `state.last_good`
/// as the source parameters. Cold-start (no prior good frame) emits
/// silence.
///
/// This is the decode-side counterpart to the
/// [`crate::ambe_plus2_wire::dequantize::Decoded::Erasure`] /
/// [`crate::imbe_wire::dequantize::Decoded::Erasure`] variants —
/// consumers call it when the wire layer signals an erasure without
/// constructing any `MbeParams` of their own. `phase_mode` selects
/// baseline vs. AMBE+ phase handling; the codec-generation modules
/// pre-pick the right value for their generation.
///
/// The underlying synth still advances per-frame state (`s_e`, `τ_M`,
/// `ε_R`, phase/noise substates), so re-acquisition remains smooth
/// after the erasure clears.
pub fn synthesize_repeat_with_mode(
    phase_mode: PhaseMode,
    state: &mut SynthState,
) -> [i16; FRAME_SAMPLES] {
    // Force Repeat disposition by setting bad_pitch on a local err
    // copy; do not touch state.err so consumer-driven error accounting
    // is preserved.
    let err = FrameErrorContext { bad_pitch: true, ..state.err };
    let params = MbeParams::silence();
    synthesize_frame_with_mode(&params, &err, state.gamma_w, phase_mode, state)
}

/// Baseline-phase convenience wrapper around
/// [`synthesize_repeat_with_mode`]. Used by
/// [`crate::codecs::ambe::synthesize_repeat`].
pub fn synthesize_repeat(state: &mut SynthState) -> [i16; FRAME_SAMPLES] {
    synthesize_repeat_with_mode(PhaseMode::Baseline, state)
}

/// AMBE+ phase-regen wrapper around [`synthesize_repeat_with_mode`].
/// Used by [`crate::codecs::ambe_plus::synthesize_repeat`] and
/// [`crate::codecs::ambe_plus2::synthesize_repeat`].
pub fn synthesize_repeat_ambe_plus(state: &mut SynthState) -> [i16; FRAME_SAMPLES] {
    synthesize_repeat_with_mode(PhaseMode::AmbePlus, state)
}

/// Core of [`synthesize_frame`] / [`synthesize_frame_ambe_plus`]
/// parameterized by [`PhaseMode`].
fn synthesize_frame_with_mode(
    params: &MbeParams,
    err: &FrameErrorContext,
    gamma_w: f64,
    phase_mode: PhaseMode,
    state: &mut SynthState,
) -> [i16; FRAME_SAMPLES] {
    // §1.11 (full-rate) vs §2.8 (half-rate) use different repeat/mute
    // thresholds and ε_R recurrences. In this codebase `PhaseMode::AmbePlus`
    // is uniquely used by the half-rate Phase 2 path (AMBE+/AMBE+2 codecs),
    // and `PhaseMode::Baseline` is uniquely used by full-rate IMBE.
    let (disp, epsilon_r) = match phase_mode {
        PhaseMode::Baseline => frame_disposition(err, state.epsilon_r),
        PhaseMode::AmbePlus => frame_disposition_halfrate(err, state.epsilon_r),
    };
    // Beyond-spec error-rate freeze on Repeat (gap 0021, opt-in via
    // `chip_compat`). When the disposition is Repeat, JMBE's
    // `IMBEModelParameters.copy()` carries `previous.errorRate` into
    // the next frame instead of advancing the smoothed value. This
    // prevents `ε_R` from climbing past the 0.0875 Mute threshold
    // during runs of high-error chip-encoded frames. Mute is checked
    // first by `frame_disposition`, so the freeze only applies once
    // Mute has not already been triggered.
    state.epsilon_r = if state.chip_compat && disp == FrameDisposition::Repeat {
        state.epsilon_r
    } else {
        epsilon_r
    };
    state.last_disposition = Some(disp);

    // Track consecutive Repeat/Mute frames for the optional reset.
    let next_repeat_count = match disp {
        FrameDisposition::Use => 0,
        FrameDisposition::Repeat | FrameDisposition::Mute => state.repeat_count + 1,
    };

    // Beyond-spec: when `repeat_reset_after` is set and the threshold
    // is exceeded, fall back to a default-fundamental + amps=1.0 frame
    // for substitution (matches JMBE's `MAX_HEADROOM_THRESHOLD = 3`
    // behavior). Spec-faithful path leaves this `None` and uses
    // literal Eq. 99–104 indefinitely (gap 0022 resolution).
    let force_default = matches!(disp, FrameDisposition::Repeat | FrameDisposition::Mute)
        && state
            .repeat_reset_after
            .is_some_and(|n| state.repeat_count >= n);
    // JMBE-cyclic behavior: when force_default fires, the substituted
    // params become this frame's `last_good`, and `repeat_count` resets
    // to 0 so subsequent frames inherit from the substituted defaults
    // (rather than continuously substituting on every frame). This
    // reproduces JMBE's `IMBEModelParameters.copy()` path where the
    // new instance's `repeatCount` is default-init at 0 when defaults
    // are substituted. Without this reset, our output settles into a
    // continuous default-substitute drone instead of the cyclical
    // recovery pattern JMBE produces.
    let next_repeat_count = if force_default { 0 } else { next_repeat_count };

    // Choose which set of parameters drives this frame's synthesis.
    // On Mute we still want to advance every cross-frame piece of
    // state (so re-acquisition works after the mute clears) but we
    // emit comfort noise — Eq. 99–104 still apply to substitute the
    // last good frame's values into the synth pipeline so its state
    // evolves smoothly.
    let (omega_0, l, voiced_arr, m_tilde_arr) = if force_default {
        // Default comfort-tone substitute: ω₀ ≈ 0.2985π (≈ 119 Hz F0,
        // ~3.4 ms period — JMBE / SDRTrunk's IMBEFundamentalFrequency.
        // DEFAULT, b̂₀ = 134), L = 30, **all bands UNVOICED**, amps = 1.0.
        // Matches JMBE's `IMBEModelParameters.copy()` default-substitute
        // branch: `setVoicingDecisions(new boolean[lplus1])` is
        // boolean-default-init = all `false`. Synthesizing all-voiced
        // produces a buzzy 30-harmonic tone; all-unvoiced produces soft
        // noise that's perceptually closer to JMBE's recovery output and
        // PESQ-scores materially higher on chip-encoded P25 traffic.
        let l = 30u8;
        let voiced = [false; L_MAX as usize];
        let mut m_tilde = [0f32; L_MAX as usize];
        for i in 0..l as usize {
            // voiced[i] stays false — JMBE-style all-unvoiced default.
            m_tilde[i] = 1.0;
        }
        (DEFAULT_OMEGA_0, l, voiced, m_tilde)
    } else {
        match (disp, &state.last_good) {
            (FrameDisposition::Use, _) => {
                let l = params.harmonic_count();
                let mut voiced = [false; L_MAX as usize];
                let mut m_tilde = [0f32; L_MAX as usize];
                voiced[..l as usize].copy_from_slice(params.voiced_slice());
                m_tilde[..l as usize].copy_from_slice(params.amplitudes_slice());
                (params.omega_0(), l, voiced, m_tilde)
            }
            (FrameDisposition::Repeat | FrameDisposition::Mute, Some(prev)) => {
                (prev.omega_0, prev.l, prev.voiced, prev.m_tilde)
            }
            (FrameDisposition::Repeat | FrameDisposition::Mute, None) => {
                // First-frame edge case: nothing to repeat. Fall back
                // to current frame's params (treats the disposition
                // as Use).
                let l = params.harmonic_count();
                let mut voiced = [false; L_MAX as usize];
                let mut m_tilde = [0f32; L_MAX as usize];
                voiced[..l as usize].copy_from_slice(params.voiced_slice());
                m_tilde[..l as usize].copy_from_slice(params.amplitudes_slice());
                (params.omega_0(), l, voiced, m_tilde)
            }
        }
    };
    state.repeat_count = next_repeat_count;

    // §1.10 enhancement: M̃_l → M̄_l, update S_E.
    let (m_bar_full, s_e_new) = analysis::profile::time(
        analysis::profile::Stage::Enhance,
        || enhance_spectral_amplitudes(&m_tilde_arr[..l as usize], omega_0, state.s_e),
    );
    state.s_e = s_e_new;

    // §1.11.3 V/UV-and-amplitude smoothing.
    let mut smoothed = analysis::profile::time(analysis::profile::Stage::Smoothing, || {
        apply_smoothing(
            &m_bar_full[..l as usize],
            &voiced_arr[..l as usize],
            state.s_e,
            state.epsilon_r,
            err.epsilon_t,
            err.epsilon_4,
            state.tau_m,
        )
    });
    state.tau_m = smoothed.tau_m;

    // Beyond-spec spectral-discontinuity clamp (gap 0026). DVSI
    // AMBE-3000R attenuates frames where the spectral envelope grows
    // abruptly relative to the previous frame, *or* where FEC reported
    // ε₀ ≥ 2 (the boundary where Golay correction confidence drops).
    // §1.11.3's τ_M/γ_M is magnitude-based (Eq. 115-116) and doesn't
    // fire when A_M < τ_M, missing both cases. Triggers, per the 2-D
    // b̂₀ sweep (2026-05-14, see SPECTRAL_CLAMP_L_THRESHOLD):
    //   - |L_curr − L_prev| > SPECTRAL_CLAMP_L_THRESHOLD (bidirectional).
    //     Small |ΔL| (≤ 5) shows no chip clamp; medium |ΔL| (6–25) has
    //     the chip running *louder* than us (overlap-add redistribution
    //     across the boundary); only |ΔL| > 25 consistently shows the
    //     chip attenuating. Down-jumps clamp at the same magnitude.
    //   - err.epsilon_0 ≥ 2: covers Repeat frames (the spec-faithful
    //     Repeat path replays last_good's params, so ΔL = 0 and the
    //     jump detector misses #3537 f=773-style cases). Real-world
    //     noisy field audio hits this branch on every Golay-corrected
    //     2-3-error frame, which is why it's the dominant useful
    //     trigger on noisy traffic.
    // Enabled either by the umbrella `chip_compat` flag (paired with
    // the gap-0021 ε_R freeze for full chip-interop) or by the
    // standalone `chip_compat_spectral_clamp` flag (for consumers who
    // want only the clamp). Both default off, so spec-faithful mode is
    // unchanged.
    //
    // **Half-rate (PhaseMode::AmbePlus) only.** The empirical formula
    // and constants are fit against the AMBE-3000R chip at
    // PKT_RATET 33; that chip is an AMBE+2 codec and is not a true
    // IMBE oracle. Applying this heuristic to full-rate IMBE
    // (PhaseMode::Baseline) is speculative — we have no IMBE chip
    // oracle to justify it. Strict-spec IMBE behavior remains the
    // default for full-rate regardless of chip_compat.
    if matches!(phase_mode, PhaseMode::AmbePlus)
        && (state.chip_compat || state.chip_compat_spectral_clamp)
    {
        let l_jump = state.prev_l > 0
            && (i16::from(l) - i16::from(state.prev_l)).abs() > SPECTRAL_CLAMP_L_THRESHOLD;
        let err_jump = err.epsilon_0 >= 2;
        if l_jump || err_jump {
            for v in smoothed.m_bar.iter_mut().take(l as usize) {
                *v = (f64::from(*v) * SPECTRAL_CLAMP_GAMMA) as f32;
            }
        }
    }
    state.prev_l = l;

    // §1.12.1 unvoiced + §1.12.2 voiced + Eq. 142 sum.
    let s_uv = analysis::profile::time(analysis::profile::Stage::SynthUnvoiced, || {
        synthesize_unvoiced(
            omega_0,
            &smoothed.m_bar[..l as usize],
            &smoothed.v_bar[..l as usize],
            gamma_w,
            &mut state.unvoiced,
        )
    });
    let noise_samples = voiced_noise_samples(&state.unvoiced);
    let lcg_m = state.unvoiced.lcg_m;
    let s_v = analysis::profile::time(analysis::profile::Stage::SynthVoiced, || {
        synthesize_voiced_with_lcg_modulus(
            omega_0,
            &smoothed.m_bar[..l as usize],
            &smoothed.v_bar[..l as usize],
            &noise_samples,
            phase_mode,
            lcg_m,
            &mut state.voiced,
        )
    });

    analysis::profile::time(analysis::profile::Stage::SynthMix, || {
        let mut s = [0f64; FRAME_SAMPLES];
        if disp == FrameDisposition::Mute {
            // §1.11.2: bypass synthesis output, emit small-amplitude
            // noise (the spec's primary recommendation; "true silence" is
            // the alternative). State above already advanced normally.
            // Reuse the unvoiced LCG's noise window (already advanced by
            // `synthesize_unvoiced`) so we don't burn extra entropy or
            // need a separate RNG. The noise window stores 209 samples
            // covering n=−104..104; take the central 160 (offsets 24..184
            // → n=−80..79). Center via LCG_MEAN and scale to a comfort
            // level matching JMBE / SDRTrunk.
            for i in 0..FRAME_SAMPLES {
                let raw = state.unvoiced.noise_window[24 + i];
                s[i] = (raw - LCG_MEAN) * MUTE_NOISE_GAIN;
            }
        } else {
            for i in 0..FRAME_SAMPLES {
                s[i] = s_uv[i] + s_v[i];
            }
        }

        // Snapshot for next-frame substitution. Storing (ω₀, L, voiced,
        // M̃) is enough — rerunning §1.10 enhancement on M̃ with the
        // current S_E reproduces M̄ exactly, so no separate M̄ snapshot
        // is needed.
        let mut snap_voiced = [false; L_MAX as usize];
        let mut snap_m_tilde = [0f32; L_MAX as usize];
        snap_voiced[..l as usize].copy_from_slice(&voiced_arr[..l as usize]);
        snap_m_tilde[..l as usize].copy_from_slice(&m_tilde_arr[..l as usize]);
        state.last_good = Some(LastGoodFrame {
            omega_0,
            l,
            voiced: snap_voiced,
            m_tilde: snap_m_tilde,
        });

        pcm_from_f64(&s)
    })
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn synth_window_endpoints_are_zero() {
        // Per Annex I header, wS(±105) = 0.
        assert_eq!(synth_window(-105), 0.0);
        assert_eq!(synth_window(105), 0.0);
    }

    #[test]
    fn synth_window_symmetric() {
        for n in 1..=105 {
            assert_eq!(
                synth_window(-n),
                synth_window(n),
                "wS(-{n}) != wS({n})"
            );
        }
    }

    #[test]
    fn synth_window_peaks_near_zero() {
        // The window is a tapered shape that peaks near n=0. Verify
        // the centre is the maximum (or among the max set).
        let centre = synth_window(0);
        for n in -105..=105 {
            assert!(
                synth_window(n) <= centre + 1e-6,
                "wS({n}) = {} exceeds wS(0) = {centre}",
                synth_window(n)
            );
        }
    }

    #[test]
    fn synth_window_out_of_range_is_zero() {
        // Spec treats wS outside [-105, 105] as zero (used in the OLA
        // formula in §1.12).
        assert_eq!(synth_window(-106), 0.0);
        assert_eq!(synth_window(106), 0.0);
        assert_eq!(synth_window(1000), 0.0);
        assert_eq!(synth_window(-1000), 0.0);
    }

    #[test]
    fn synth_window_table_length_matches_constant() {
        assert_eq!(IMBE_SYNTH_WINDOW.len(), SYNTH_WINDOW_LEN);
        assert_eq!(SYNTH_WINDOW_LEN, 211);
    }

    // ---- §1.10 enhancement -----------------------------------------------

    fn sum_sq(arr: &[f32]) -> f64 {
        arr.iter().map(|&v| f64::from(v) * f64::from(v)).sum()
    }

    #[test]
    fn enhance_silence_passes_through_and_floors_s_e() {
        // All-zero M̃_l → R_M0 = 0 → silence-frame degenerate path.
        let m_tilde = vec![0.0f32; 30];
        let (m_bar, s_e) = enhance_spectral_amplitudes(&m_tilde, 0.1, INIT_S_E);
        for i in 0..30 {
            assert_eq!(m_bar[i], 0.0, "M̄_{i} = {}", m_bar[i]);
        }
        // S_E recurrence with R_M0 = 0: 0.95·75000 = 71250 (above floor).
        assert!((s_e - 71250.0).abs() < 1e-6, "S_E = {s_e}");
    }

    #[test]
    fn enhance_preserves_energy_within_a_few_percent() {
        // The Eq. 109–110 rescale is *exact* energy-preserving in
        // double precision before the f64 → f32 cast; we just check
        // the round-tripped energy matches within a small tolerance.
        let l = 24usize;
        let m_tilde: Vec<f32> = (1..=l).map(|i| (i as f32 * 0.3).sin().abs() + 0.1).collect();
        let omega_0 = std::f32::consts::PI / 12.0;
        let (m_bar, _) = enhance_spectral_amplitudes(&m_tilde, omega_0, INIT_S_E);
        let energy_in = sum_sq(&m_tilde);
        let energy_out = sum_sq(&m_bar[..l]);
        let rel_err = (energy_out - energy_in).abs() / energy_in;
        assert!(rel_err < 1e-4, "energy drift {rel_err} (in={energy_in}, out={energy_out})");
    }

    #[test]
    fn enhance_s_e_floor_engages() {
        // After many quiet frames, S_E should converge to the floor.
        let m_tilde = vec![0.0f32; 30];
        let mut s_e = INIT_S_E;
        for _ in 0..200 {
            let (_, new_s_e) = enhance_spectral_amplitudes(&m_tilde, 0.1, s_e);
            s_e = new_s_e;
        }
        assert!((s_e - S_E_FLOOR).abs() < 1e-6, "S_E = {s_e}");
    }

    #[test]
    fn enhance_s_e_recurrence_responds_to_input_energy() {
        // Pumping in higher-energy frames should drive S_E up from
        // its initial value (with the 0.05 mixing factor).
        let m_tilde: Vec<f32> = (1..=20).map(|_| 100.0).collect(); // R_M0 ≈ 200000
        let (_, s_e1) = enhance_spectral_amplitudes(&m_tilde, 0.1, INIT_S_E);
        // S_E = 0.95·75000 + 0.05·200000 = 71250 + 10000 = 81250
        assert!((s_e1 - 81250.0).abs() < 1e-2, "S_E = {s_e1}");
    }

    #[test]
    fn enhance_low_harmonic_bypass_holds_per_section_1_10() {
        // For l with 8·l ≤ L̃, the W_l weighting is bypassed (Eq. 108
        // first branch). The bypass gives M̄_l = M̃_l *before* the
        // global γ rescale. Verify the per-harmonic ratio M̄_l / M̃_l
        // is identical across all bypassed l (i.e., they all get the
        // same γ scaling and no per-l W weighting).
        let l = 56usize; // many harmonics so 8·l ≤ L̃ for several values
        let m_tilde: Vec<f32> = (1..=l).map(|i| 1.0 + (i as f32 * 0.05).sin()).collect();
        let omega_0 = std::f32::consts::PI / 28.0;
        let (m_bar, _) = enhance_spectral_amplitudes(&m_tilde, omega_0, INIT_S_E);
        // Bypassed harmonics: l = 1..=⌊L/8⌋ = 1..=7. They all share
        // the same final γ multiplier.
        let r0 = m_bar[0] / m_tilde[0];
        for i in 1..7 {
            let r = m_bar[i] / m_tilde[i];
            assert!(
                (r - r0).abs() < 1e-4,
                "bypassed harmonic l={}: ratio {r} vs {r0}",
                i + 1
            );
        }
    }

    // ---- §1.11 smoothing / disposition -----------------------------------

    fn err(epsilon_0: u8, epsilon_t: u8, epsilon_4: u8) -> FrameErrorContext {
        FrameErrorContext { epsilon_0, epsilon_4, epsilon_t, bad_pitch: false }
    }

    #[test]
    fn frame_disposition_clean_frame_uses() {
        let (d, er) = frame_disposition(&err(0, 0, 0), 0.0);
        assert_eq!(d, FrameDisposition::Use);
        assert_eq!(er, 0.0);
    }

    #[test]
    fn frame_disposition_bad_pitch_repeats() {
        let mut e = err(0, 0, 0);
        e.bad_pitch = true;
        let (d, _) = frame_disposition(&e, 0.0);
        assert_eq!(d, FrameDisposition::Repeat);
    }

    #[test]
    fn frame_disposition_joint_error_threshold_repeats() {
        // ε₀ ≥ 2 AND ε_T ≥ 10 + 40·ε_R(0). With ε_R(prev) = 0 and ε_T = 12,
        // ε_R(0) = 0.05·12/144 ≈ 0.00417. Threshold = 10 + 40·0.00417 ≈ 10.17.
        // ε_T = 12 ≥ 10.17 → repeat.
        let (d, _) = frame_disposition(&err(2, 12, 0), 0.0);
        assert_eq!(d, FrameDisposition::Repeat);
        // ε₀ = 1 below threshold → use.
        let (d, _) = frame_disposition(&err(1, 12, 0), 0.0);
        assert_eq!(d, FrameDisposition::Use);
    }

    #[test]
    fn frame_disposition_high_error_rate_mutes() {
        // ε_R prev high enough that the smoothed value crosses 0.0875.
        let (d, _) = frame_disposition(&err(0, 50, 0), 0.5);
        assert_eq!(d, FrameDisposition::Mute);
    }

    /// Spec-faithful default: a run of high-error frames advances ε_R
    /// every frame and eventually trips Mute. Reproduces the chip.bit
    /// failure mode (gap 0021).
    #[test]
    fn synth_baseline_drives_chip_pattern_into_mute() {
        let mut state = SynthState::new();
        // Default: chip_compat = false. ε_R advances every frame.
        assert!(!state.chip_compat());
        // Chip.bit's deterministic 13-error pattern, ε_0 = 3 (≥ 2).
        let err_ctx = err(3, 13, 1);
        let voiced = vec![true; crate::mbe_params::L_MIN as usize];
        let amps = vec![1.0f32; crate::mbe_params::L_MIN as usize];
        let params = MbeParams::new(
            2.0 * core::f32::consts::PI / 50.0,
            crate::mbe_params::L_MIN,
            &voiced,
            &amps,
        )
        .unwrap();
        // Run enough frames for ε_R to climb past 0.0875. Steady-state
        // is 0.05·13/144 / (1−0.95) = 0.0903, just above the 0.0875
        // Mute threshold, but it takes ≈70 frames to converge close
        // enough (ε_R = 0.0903·(1−0.95^n), need n ≳ 70 for ε_R > 0.0875).
        for _ in 0..120 {
            let _ = synthesize_frame(&params, &err_ctx, GAMMA_W, &mut state);
        }
        assert_eq!(
            state.last_disposition(),
            Some(FrameDisposition::Mute),
            "spec-faithful path should mute under sustained chip.bit errors"
        );
    }

    /// chip_compat = true: same sustained 13-error pattern stays in
    /// Repeat indefinitely. ε_R never advances past the first frame,
    /// Mute never trips. Mirrors JMBE's `IMBEModelParameters.copy()`
    /// `setErrorRate(previous.getErrorRate())` semantics.
    #[test]
    fn synth_chip_compat_freezes_error_rate_on_repeat() {
        let mut state = SynthState::new();
        state.set_chip_compat(true);
        let err_ctx = err(3, 13, 1);
        let voiced = vec![true; crate::mbe_params::L_MIN as usize];
        let amps = vec![1.0f32; crate::mbe_params::L_MIN as usize];
        let params = MbeParams::new(
            2.0 * core::f32::consts::PI / 50.0,
            crate::mbe_params::L_MIN,
            &voiced,
            &amps,
        )
        .unwrap();
        for _ in 0..120 {
            let _ = synthesize_frame(&params, &err_ctx, GAMMA_W, &mut state);
        }
        assert_eq!(
            state.last_disposition(),
            Some(FrameDisposition::Repeat),
            "chip_compat path should hold Repeat (ε_R frozen below Mute)"
        );
        // ε_R should be near the cold-start single-step value
        // (0.05·13/144 ≈ 0.00451), not the climbed steady-state.
        assert!(
            state.epsilon_r < 0.0875,
            "frozen ε_R should stay below Mute threshold (got {})",
            state.epsilon_r
        );
    }

    #[test]
    fn smoothing_low_error_uses_infinite_v_m() {
        // Branch 1: ε_R ≤ 0.005 AND ε_T ≤ 4 → V_M = ∞. Even very tall
        // M̄_l should not flip an unvoiced harmonic to voiced.
        let m_bar = vec![1e6f32; 9];
        let v_tilde = vec![false; 9];
        let s = apply_smoothing(&m_bar, &v_tilde, 75000.0, 0.001, 2, 0, INIT_TAU_M);
        for i in 0..9 {
            assert!(!s.v_bar[i], "v̄_{i} flipped to voiced under V_M = ∞");
        }
    }

    #[test]
    fn smoothing_v_uv_override_when_amplitude_exceeds_threshold() {
        // Force branch 3 (V_M = 1.414 · S_E^0.375). With S_E small, V_M
        // is small; large M̄ should flip to voiced.
        let s_e = 100.0; // V_M = 1.414 · 100^0.375 ≈ 1.414 · 5.84 ≈ 8.27
        let m_bar = vec![1.0f32, 100.0, 1.0, 100.0, 1.0, 100.0, 1.0, 100.0, 1.0];
        let v_tilde = vec![false; 9];
        let s = apply_smoothing(&m_bar, &v_tilde, s_e, 0.05, 10, 1, INIT_TAU_M);
        for i in 0..9 {
            let expected = m_bar[i] > 8.27;
            assert_eq!(s.v_bar[i], expected, "v̄_{i}");
        }
    }

    #[test]
    fn smoothing_amplitude_low_error_resets_tau_m() {
        // Branch 1 of Eq. 115: ε_R ≤ 0.005 AND ε_T ≤ 6 → τ_M = 20480.
        let m_bar = vec![10.0f32; 9];
        let s = apply_smoothing(&m_bar, &vec![true; 9], 1000.0, 0.001, 3, 0, 5000.0);
        assert_eq!(s.tau_m, 20480.0);
    }

    #[test]
    fn smoothing_amplitude_recurrence_under_errors() {
        // Branch 2: τ_M(0) = 6000 − 300·ε_T + τ_M(−1)
        let m_bar = vec![10.0f32; 9];
        let s = apply_smoothing(&m_bar, &vec![true; 9], 1000.0, 0.05, 8, 1, 10000.0);
        // expected τ_M = 6000 − 300·8 + 10000 = 13600
        assert!((s.tau_m - 13600.0).abs() < 1e-6);
    }

    #[test]
    fn smoothing_gamma_m_is_one_when_tau_m_exceeds_amplitude_total() {
        // A_M = Σ M̄_l. Tiny amplitudes → A_M < τ_M → γ_M = 1 (no scaling).
        let m_bar = vec![0.01f32; 9];
        let s = apply_smoothing(&m_bar, &vec![true; 9], 1000.0, 0.001, 0, 0, INIT_TAU_M);
        for i in 0..9 {
            assert!((s.m_bar[i] - 0.01).abs() < 1e-6);
        }
    }

    // ---- §1.12.1 unvoiced synthesis --------------------------------------

    #[test]
    fn lcg_first_few_values_match_recurrence() {
        // The cold-start `UnvoicedSynthState::new()` now defaults to the
        // chip's LCG (probed: 173/13849/65536, seed 60584) rather than
        // BABA-A §1.12.1's spec LCG. Verify both recurrences via their
        // first advance from the seed:
        //
        // Chip LCG: (173·60584 + 13849) mod 65536
        //         = (10481032 + 13849) mod 65536 = 10494881 mod 65536
        //         = 10494881 - 160·65536 = 10494881 - 10485760 = 9121
        let mut chip_state = UnvoicedSynthState::new();
        assert_eq!(chip_state.next_noise() as u32, 9121);

        // Spec LCG: u(n+1) = (171·u(n) + 11213) mod 53125, u(−105) = 3147.
        // First call: (171·3147 + 11213) mod 53125
        //           = (538137 + 11213) mod 53125 = 549350 mod 53125 = 18100
        let (a, c, m, seed) = UnvoicedNoiseGen::SpecLcg.params();
        assert_eq!((a.wrapping_mul(seed).wrapping_add(c)) % m, 18100);
        // Second: (171·18100 + 11213) mod 53125 = 25063
        assert_eq!((a.wrapping_mul(18100).wrapping_add(c)) % m, 25063);
    }

    #[test]
    fn dft_idft_roundtrip_recovers_input() {
        // A real input that fits in n=−104..=104, surrounded by zero.
        let mut input = [0f64; 209];
        for i in 0..209 {
            let n = i as i32 - 104;
            input[i] = (f64::from(n) * 0.05).sin() * 100.0;
        }
        let (re, im) = dft_256_windowed(&input);
        let recovered = idft_256(&re, &im);
        // The IDFT spans n=−128..127 (256 values), our input spans n=-104..=104.
        // Recovered values for n in [-104, 104] should match the input
        // at the same n-position. recovered[(n+128)] = input[(n+104)].
        for n in -104..=104i32 {
            let in_val = input[(n + 104) as usize];
            let out_val = recovered[(n + 128) as usize];
            assert!(
                (out_val - in_val).abs() < 1e-6,
                "n={n}: in={in_val}, out={out_val}"
            );
        }
    }

    #[test]
    fn unvoiced_silence_input_produces_silence_output() {
        let m_bar = vec![0.0f32; 12];
        let v_bar = vec![false; 12];
        let mut state = UnvoicedSynthState::new();
        let pcm = synthesize_unvoiced(0.2, &m_bar, &v_bar, GAMMA_W, &mut state);
        for (i, &s) in pcm.iter().enumerate() {
            assert!(s.abs() < 1e-6, "n={i}: {s}");
        }
    }

    #[test]
    fn unvoiced_all_voiced_first_frame_produces_silence_output() {
        // Voiced bands → Ũ_w = 0 in those bands; with no unvoiced bands
        // there's no signal. Previous IDFT is also zero (init), so OLA
        // output is silence.
        let m_bar = vec![10.0f32; 12];
        let v_bar = vec![true; 12];
        let mut state = UnvoicedSynthState::new();
        let pcm = synthesize_unvoiced(0.2, &m_bar, &v_bar, GAMMA_W, &mut state);
        for (i, &s) in pcm.iter().enumerate() {
            assert!(s.abs() < 1e-6, "n={i}: {s}");
        }
    }

    #[test]
    fn unvoiced_nonzero_input_produces_nonzero_output() {
        let m_bar = vec![10.0f32; 12];
        let v_bar = vec![false; 12];
        let mut state = UnvoicedSynthState::new();
        // First frame's OLA blends prev IDFT (zero) with current IDFT.
        // So we run two frames and check the second has non-trivial energy.
        let _ = synthesize_unvoiced(0.2, &m_bar, &v_bar, GAMMA_W, &mut state);
        let pcm = synthesize_unvoiced(0.2, &m_bar, &v_bar, GAMMA_W, &mut state);
        let energy: f64 = pcm.iter().map(|&s| s * s).sum();
        assert!(energy > 1e-3, "energy too low: {energy}");
    }

    #[test]
    fn unvoiced_state_advances_between_frames() {
        let m_bar = vec![1.0f32; 9];
        let v_bar = vec![false; 9];
        let mut state = UnvoicedSynthState::new();
        let lcg_init = state.lcg;
        let _ = synthesize_unvoiced(0.2, &m_bar, &v_bar, GAMMA_W, &mut state);
        let lcg_after_1 = state.lcg;
        assert_ne!(lcg_after_1, lcg_init);
        let _ = synthesize_unvoiced(0.2, &m_bar, &v_bar, GAMMA_W, &mut state);
        let lcg_after_2 = state.lcg;
        assert_ne!(lcg_after_2, lcg_after_1);
    }

    #[test]
    fn unvoiced_output_length_is_one_frame() {
        let m_bar = vec![1.0f32; 9];
        let v_bar = vec![false; 9];
        let mut state = UnvoicedSynthState::new();
        let pcm = synthesize_unvoiced(0.2, &m_bar, &v_bar, GAMMA_W, &mut state);
        assert_eq!(pcm.len(), FRAME_SAMPLES);
        assert_eq!(FRAME_SAMPLES, 160);
    }

    #[test]
    fn smoothing_gamma_m_clamps_loud_frames() {
        // Loud amplitudes → A_M > τ_M → γ_M = τ_M / A_M < 1, all M̄_l scaled
        // by the same γ_M.
        let m_bar = vec![10000.0f32; 9];
        let v_tilde = vec![true; 9];
        let s = apply_smoothing(&m_bar, &v_tilde, 1000.0, 0.001, 0, 0, INIT_TAU_M);
        // A_M = 90000, τ_M = 20480, γ_M ≈ 0.2276.
        let expected_gamma = 20480.0 / 90000.0;
        for i in 0..9 {
            let ratio = s.m_bar[i] / 10000.0;
            assert!(
                (ratio - expected_gamma as f32).abs() < 1e-3,
                "γ_M ratio at l={}: {ratio} vs {expected_gamma}",
                i + 1
            );
        }
    }

    // ---- §1.12.2 voiced synthesis ----------------------------------------

    #[test]
    fn voiced_state_init_matches_annex_a() {
        let s = VoicedSynthState::new();
        assert_eq!(s.prev_l, 30);
        assert!((s.prev_omega_0 - 0.02985 * PI64).abs() < 1e-12);
        for l in 1..=L_MAX as usize {
            assert_eq!(s.phi[l], 0.0);
            assert_eq!(s.psi[l], 0.0);
            assert_eq!(s.prev_m_bar[l], 0.0);
            assert!(!s.prev_v_bar[l]);
        }
    }

    #[test]
    fn wrap_phase_into_principal_range() {
        let cases = [
            (0.0f64, 0.0),
            (PI64, -PI64),
            (-PI64, -PI64),
            (3.0 * PI64, -PI64),
            (-3.0 * PI64, -PI64),
        ];
        for (input, expected) in cases {
            let w = wrap_phase(input);
            assert!(
                (w - expected).abs() < 1e-12,
                "wrap({input}) = {w}, expected {expected}"
            );
        }
        // Verify range invariant for arbitrary inputs.
        for k in -10..=10 {
            let raw = (k as f64) * 0.7 * PI64;
            let w = wrap_phase(raw);
            assert!((-PI64..PI64).contains(&w), "wrap({raw}) = {w}");
        }
    }

    #[test]
    fn voiced_all_unvoiced_produces_silence() {
        // (UV, UV) for every harmonic → branch UvUv → s̃_{v,l} = 0.
        let m_bar = vec![10.0f32; 12];
        let v_bar = vec![false; 12];
        let noise = [0f64; L_MAX as usize];
        let mut state = VoicedSynthState::new();
        // First call: prev is also init (all UV by Annex A).
        let pcm = synthesize_voiced(0.2, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state);
        for &s in pcm.iter() {
            assert!(s.abs() < 1e-9);
        }
    }

    #[test]
    fn voiced_silence_amplitudes_produces_silence() {
        let m_bar = vec![0f32; 12];
        let v_bar = vec![true; 12];
        let noise = [0f64; L_MAX as usize];
        let mut state = VoicedSynthState::new();
        // First frame: prev is all UV (Annex A init), so harmonics use
        // Branch::UvV with M̄_curr = 0 → output is zero.
        let pcm = synthesize_voiced(0.2, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state);
        for &s in pcm.iter() {
            assert!(s.abs() < 1e-9, "{s}");
        }
    }

    #[test]
    fn voiced_state_advances_phase_for_all_56_harmonics() {
        // ψ_l(0) = ψ_l(−1) + (ω_prev + ω_curr)·l·N/2 — for ALL l in 1..=56,
        // regardless of L̃ or V/UV.
        let m_bar = vec![1.0f32; 9]; // L = 9
        let v_bar = vec![false; 9];
        let noise = [0f64; L_MAX as usize];
        let mut state = VoicedSynthState::new();
        let psi_init = state.psi;
        let _ = synthesize_voiced(0.2, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state);
        // Every harmonic's ψ should have advanced.
        for l in 1..=L_MAX as usize {
            assert!(
                state.psi[l] > psi_init[l] - 1e-12,
                "ψ_{l} did not advance: {} → {}",
                psi_init[l],
                state.psi[l]
            );
            assert!(state.psi[l] != psi_init[l], "ψ_{l} unchanged");
        }
    }

    #[test]
    fn voiced_steady_state_pure_tone() {
        // Steady-state sinusoid: identical M̄, ω₀, V/UV across two
        // frames should produce a near-pure-tone-like waveform on the
        // second frame (after warm-up).
        let omega_0 = 0.2f32;
        let mut m_bar = [0f32; 12];
        m_bar[0] = 100.0; // fundamental only, harmonic 1 voiced
        let mut v_bar = [false; 12];
        v_bar[0] = true;
        let noise = [0f64; L_MAX as usize];
        let mut state = VoicedSynthState::new();
        let _ = synthesize_voiced(omega_0, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state);
        let pcm = synthesize_voiced(omega_0, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state);

        // Output should have measurable energy concentrated around
        // ω₀. We just sanity-check non-zero variance.
        let mean: f64 = pcm.iter().sum::<f64>() / pcm.len() as f64;
        let var: f64 = pcm.iter().map(|&s| (s - mean).powi(2)).sum::<f64>()
            / pcm.len() as f64;
        assert!(var > 100.0, "voiced steady-state variance too low: {var}");
    }

    #[test]
    fn voiced_phase_carries_across_frames_for_continuity() {
        // The φ_l update should keep the synthesized waveform roughly
        // continuous between consecutive steady-state frames. We
        // compare the last sample of frame N with the first sample of
        // frame N+1 — they should be close (within a few % of peak).
        let omega_0 = 0.15f32;
        let mut m_bar = [0f32; 12];
        m_bar[0] = 100.0;
        m_bar[1] = 60.0;
        let mut v_bar = [false; 12];
        v_bar[0] = true;
        v_bar[1] = true;
        let noise = [0f64; L_MAX as usize];
        let mut state = VoicedSynthState::new();
        // Warm up.
        let _ = synthesize_voiced(omega_0, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state);
        let _ = synthesize_voiced(omega_0, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state);
        let pcm_a = synthesize_voiced(omega_0, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state);
        let pcm_b = synthesize_voiced(omega_0, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state);
        let last_a = pcm_a[FRAME_SAMPLES - 1];
        let first_b = pcm_b[0];
        // Approximate continuity: difference well below the typical
        // sample magnitude.
        let typical_mag = pcm_a.iter().map(|s| s.abs()).fold(0f64, f64::max);
        assert!(
            (last_a - first_b).abs() < typical_mag,
            "phase discontinuity: last={last_a}, first={first_b}, max={typical_mag}"
        );
    }

    #[test]
    fn voiced_output_length_is_one_frame() {
        let m_bar = [0f32; 9];
        let v_bar = [false; 9];
        let noise = [0f64; L_MAX as usize];
        let mut state = VoicedSynthState::new();
        let pcm = synthesize_voiced(0.2, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state);
        assert_eq!(pcm.len(), FRAME_SAMPLES);
    }

    #[test]
    fn voiced_ambe_plus_mode_diverges_from_baseline_on_high_l() {
        // Same MbeParams and cold-start state run through Baseline
        // (Eq. 141 noise-perturbed phase) and AmbePlus (Hilbert regen)
        // must produce different PCM on high-l harmonics because the
        // phase computation differs there. Low-l region (l ≤ L/4) is
        // common to both and should agree.
        let omega_0 = 0.15f32;
        // L = 12, so L/4 = 3. Harmonics 1..=3 use the common ψ branch;
        // 4..=12 differ between modes.
        let mut m_bar = [0f32; 12];
        for (l, m) in m_bar.iter_mut().enumerate() {
            *m = 100.0 * ((l as f32) + 1.0).sqrt();
        }
        let v_bar = [true; 12];
        let noise = [100f64; L_MAX as usize]; // non-trivial noise so Eq. 141 fires

        let mut state_base = VoicedSynthState::new();
        let _ = synthesize_voiced(
            omega_0, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state_base,
        );
        let pcm_base = synthesize_voiced(
            omega_0, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state_base,
        );

        let mut state_amb = VoicedSynthState::new();
        let _ = synthesize_voiced(
            omega_0, &m_bar, &v_bar, &noise, PhaseMode::AmbePlus, &mut state_amb,
        );
        let pcm_amb = synthesize_voiced(
            omega_0, &m_bar, &v_bar, &noise, PhaseMode::AmbePlus, &mut state_amb,
        );

        // Overall waveforms should differ.
        let diff_rms: f64 = pcm_base
            .iter()
            .zip(pcm_amb.iter())
            .map(|(a, b)| (a - b).powi(2))
            .sum::<f64>()
            / FRAME_SAMPLES as f64;
        assert!(
            diff_rms.sqrt() > 1.0,
            "Baseline and AmbePlus produced identical voiced output — phase mode not wired through"
        );
    }

    #[test]
    fn voiced_ambe_plus_matches_baseline_on_all_unvoiced() {
        // When no harmonics are voiced, both phase modes produce
        // identical (zero) output — the phase branch never activates.
        let m_bar = vec![100f32; 9];
        let v_bar = vec![false; 9];
        let noise = [50f64; L_MAX as usize];

        let mut state_base = VoicedSynthState::new();
        let pcm_base = synthesize_voiced(
            0.2, &m_bar, &v_bar, &noise, PhaseMode::Baseline, &mut state_base,
        );
        let mut state_amb = VoicedSynthState::new();
        let pcm_amb = synthesize_voiced(
            0.2, &m_bar, &v_bar, &noise, PhaseMode::AmbePlus, &mut state_amb,
        );
        for (i, (a, b)) in pcm_base.iter().zip(pcm_amb.iter()).enumerate() {
            assert!(
                (a - b).abs() < 1e-9,
                "i={i}: Baseline {a} vs AmbePlus {b} — should match on all-UV frame"
            );
        }
    }

    #[test]
    fn voiced_noise_samples_extracts_from_window() {
        let mut us = UnvoicedSynthState::new();
        us.advance_window();
        let n = voiced_noise_samples(&us);
        for l in 1..=L_MAX as usize {
            assert_eq!(n[l - 1], us.noise_window[l + 104]);
        }
    }

    // ---- Top-level synthesize_frame --------------------------------------

    fn build_params(omega_0: f32, voiced: &[bool], amplitudes: &[f32]) -> MbeParams {
        let l = voiced.len() as u8;
        MbeParams::new(omega_0, l, voiced, amplitudes).expect("valid params")
    }

    #[test]
    fn synthesize_frame_silence_input_produces_silence_output() {
        let p = build_params(0.2, &[false; 9], &[0f32; 9]);
        let err = FrameErrorContext::default();
        let mut state = SynthState::new();
        let pcm = synthesize_frame(&p, &err, GAMMA_W, &mut state);
        for s in pcm.iter() {
            assert_eq!(*s, 0);
        }
    }

    #[test]
    fn synthesize_frame_advances_all_substates() {
        let voiced = vec![true; 9];
        let amps = vec![100f32; 9];
        let p = build_params(0.2, &voiced, &amps);
        let err = FrameErrorContext::default();
        let mut state = SynthState::new();
        let s_e_init = state.s_e;
        let lcg_init = state.unvoiced.lcg;
        let psi_init = state.voiced.psi;
        let _ = synthesize_frame(&p, &err, GAMMA_W, &mut state);
        assert_ne!(state.s_e, s_e_init, "S_E unchanged");
        assert_ne!(state.unvoiced.lcg, lcg_init, "LCG unchanged");
        // ψ for at least one harmonic should have advanced.
        let any_changed = (1..=L_MAX as usize).any(|l| state.voiced.psi[l] != psi_init[l]);
        assert!(any_changed, "no ψ advanced");
        assert!(state.last_good.is_some(), "no snapshot stored");
    }

    #[test]
    fn synthesize_frame_mute_emits_low_amplitude_noise_but_advances_state() {
        // Force a Mute by pre-loading a high ε_R into state. §1.11.2's
        // primary recommendation is "random small-amplitude noise"; we
        // follow that (silence was the prior alternative).
        let voiced = vec![true; 9];
        let amps = vec![100f32; 9];
        let p = build_params(0.2, &voiced, &amps);
        let err = FrameErrorContext { epsilon_t: 100, ..Default::default() };
        let mut state = SynthState::new();
        state.epsilon_r = 0.5; // > MUTE_EPSILON_R_THRESHOLD
        let s_e_init = state.s_e;
        let pcm = synthesize_frame(&p, &err, GAMMA_W, &mut state);
        // Output is comfort noise: peak under 1500 i16 (~5% full
        // scale), and at least some samples non-zero.
        let peak = pcm.iter().map(|s| s.unsigned_abs()).max().unwrap_or(0);
        let nonzero = pcm.iter().filter(|&&s| s != 0).count();
        assert!(peak < 1500, "Mute noise peak too loud: {peak}");
        assert!(nonzero > FRAME_SAMPLES / 2, "Mute noise too sparse: {nonzero}");
        // State should still have advanced (so re-acquisition works).
        assert_ne!(state.s_e, s_e_init);
    }

    #[test]
    fn pcm_from_f64_clamps_and_rounds() {
        let mut input = [0f64; FRAME_SAMPLES];
        input[0] = 0.0;
        input[1] = 100.4; // → 100
        input[2] = 100.6; // → 101
        input[3] = -100.6; // → -101
        input[4] = 50000.0; // → clamp to 32767
        input[5] = -50000.0; // → clamp to -32768
        let out = pcm_from_f64(&input);
        assert_eq!(out[0], 0);
        assert_eq!(out[1], 100);
        assert_eq!(out[2], 101);
        assert_eq!(out[3], -101);
        assert_eq!(out[4], 32767);
        assert_eq!(out[5], -32768);
    }

    #[test]
    fn synthesize_frame_repeat_reuses_last_good() {
        // Run a clean frame to populate last_good, then a Repeat frame
        // using completely different params — output should match what
        // the clean frame would have produced for the SAME index in
        // its OWN second call (since substitution + state evolution is
        // deterministic).
        let voiced = vec![true; 9];
        let amps = vec![50f32; 9];
        let p_good = build_params(0.2, &voiced, &amps);
        let p_repeat = build_params(0.1, &vec![false; 12], &vec![1.0; 12]);
        let err_use = FrameErrorContext::default();
        let err_repeat = FrameErrorContext { bad_pitch: true, ..Default::default() };

        let mut state_a = SynthState::new();
        let _ = synthesize_frame(&p_good, &err_use, GAMMA_W, &mut state_a);
        let pcm_repeat = synthesize_frame(&p_repeat, &err_repeat, GAMMA_W, &mut state_a);

        let mut state_b = SynthState::new();
        let _ = synthesize_frame(&p_good, &err_use, GAMMA_W, &mut state_b);
        let pcm_use_again = synthesize_frame(&p_good, &err_use, GAMMA_W, &mut state_b);

        // The Repeat path used last_good (= p_good), so its output
        // should equal the second p_good frame, modulo the change in
        // ε_R (Repeat carries non-zero ε_T = 0 by default so ε_R is
        // unchanged in this test).
        for i in 0..FRAME_SAMPLES {
            assert_eq!(
                pcm_repeat[i], pcm_use_again[i],
                "n={i}: repeat={} use={}",
                pcm_repeat[i], pcm_use_again[i]
            );
        }
    }

    #[test]
    fn spectral_clamp_does_not_fire_on_full_rate_imbe() {
        // We have no IMBE chip oracle (AMBE-3000R is AMBE+2, not IMBE),
        // so the empirically-fit clamp is gated to PhaseMode::AmbePlus
        // (half-rate). Full-rate IMBE goes through PhaseMode::Baseline
        // and must NOT trigger the clamp regardless of chip_compat.
        let p_low = build_params(0.5, &vec![true; 9], &vec![50f32; 9]);
        let p_high = build_params(0.1, &vec![true; 30], &vec![50f32; 30]);
        let err = FrameErrorContext::default();

        let mut state_off = SynthState::new();
        let _ = synthesize_frame(&p_low, &err, GAMMA_W, &mut state_off);
        let pcm_off = synthesize_frame(&p_high, &err, GAMMA_W, &mut state_off);

        let mut state_on = SynthState::new();
        state_on.set_chip_compat(true);
        state_on.set_chip_compat_spectral_clamp(true);
        let _ = synthesize_frame(&p_low, &err, GAMMA_W, &mut state_on);
        let pcm_on = synthesize_frame(&p_high, &err, GAMMA_W, &mut state_on);

        assert_eq!(
            &pcm_off[..], &pcm_on[..],
            "full-rate IMBE (PhaseMode::Baseline) must be unaffected by chip_compat flags"
        );
    }

    #[test]
    fn spectral_clamp_auto_enables_under_chip_compat() {
        // The umbrella `chip_compat` flag should enable the spectral
        // clamp without needing the standalone toggle. Auto-on under
        // chip_compat is the recommended consumer profile. Uses a
        // large L jump (9→45, |ΔL|=36) to exceed the threshold of 25.
        let p_low = build_params(0.5, &vec![true; 9], &vec![50f32; 9]);
        let p_high = build_params(0.05, &vec![true; 45], &vec![50f32; 45]);
        let err = FrameErrorContext::default();

        let rms = |pcm: &[i16]| {
            (pcm.iter().map(|&s| f64::from(s).powi(2)).sum::<f64>() / pcm.len() as f64).sqrt()
        };

        let mut state_off = SynthState::new();
        let _ = synthesize_frame_ambe_plus(&p_low, &err, GAMMA_W, &mut state_off);
        let pcm_off = synthesize_frame_ambe_plus(&p_high, &err, GAMMA_W, &mut state_off);

        let mut state_on = SynthState::new();
        state_on.set_chip_compat(true);
        assert!(!state_on.chip_compat_spectral_clamp(),
            "standalone flag should stay off — clamp is enabled via chip_compat");
        let _ = synthesize_frame_ambe_plus(&p_low, &err, GAMMA_W, &mut state_on);
        let pcm_on = synthesize_frame_ambe_plus(&p_high, &err, GAMMA_W, &mut state_on);

        let rms_off = rms(&pcm_off);
        let rms_on = rms(&pcm_on);
        assert!(rms_on < rms_off,
            "chip_compat should auto-enable the clamp on large L-jump: off={rms_off:.1} on={rms_on:.1}");
    }

    #[test]
    fn spectral_clamp_default_off_leaves_synth_unchanged() {
        // With the flag off (default), the synth path must match the
        // spec-faithful behavior bit-for-bit.
        let p_low = build_params(0.5, &vec![true; 9], &vec![50f32; 9]);
        let p_high = build_params(0.1, &vec![true; 30], &vec![50f32; 30]); // L=9 → 30
        let err = FrameErrorContext::default();

        let mut state_a = SynthState::new();
        let _ = synthesize_frame_ambe_plus(&p_low, &err, GAMMA_W, &mut state_a);
        let pcm_a = synthesize_frame_ambe_plus(&p_high, &err, GAMMA_W, &mut state_a);

        let mut state_b = SynthState::new();
        assert!(!state_b.chip_compat_spectral_clamp(), "default must be off");
        let _ = synthesize_frame_ambe_plus(&p_low, &err, GAMMA_W, &mut state_b);
        let pcm_b = synthesize_frame_ambe_plus(&p_high, &err, GAMMA_W, &mut state_b);

        assert_eq!(&pcm_a[..], &pcm_b[..], "default-off must not alter synth output");
    }

    #[test]
    fn spectral_clamp_fires_on_large_l_jump_when_enabled() {
        // L jumps 9 → 45 (Δ = +36, well above the threshold of 25).
        // Clamp-on should produce smaller RMS than clamp-off.
        let p_low = build_params(0.5, &vec![true; 9], &vec![50f32; 9]);
        let p_high = build_params(0.05, &vec![true; 45], &vec![50f32; 45]);
        let err = FrameErrorContext::default();

        let rms = |pcm: &[i16]| {
            (pcm.iter().map(|&s| f64::from(s).powi(2)).sum::<f64>() / pcm.len() as f64).sqrt()
        };

        let mut state_off = SynthState::new();
        let _ = synthesize_frame_ambe_plus(&p_low, &err, GAMMA_W, &mut state_off);
        let pcm_off = synthesize_frame_ambe_plus(&p_high, &err, GAMMA_W, &mut state_off);

        let mut state_on = SynthState::new();
        state_on.set_chip_compat_spectral_clamp(true);
        let _ = synthesize_frame_ambe_plus(&p_low, &err, GAMMA_W, &mut state_on);
        let pcm_on = synthesize_frame_ambe_plus(&p_high, &err, GAMMA_W, &mut state_on);

        let rms_off = rms(&pcm_off);
        let rms_on = rms(&pcm_on);
        assert!(rms_on < rms_off, "clamp must reduce RMS on large L-jump: off={rms_off:.1} on={rms_on:.1}");
        // Roughly 0.75× per sweep data — allow generous bounds.
        let ratio = rms_on / rms_off;
        assert!(ratio < 0.95, "expected ratio < 0.95, got {ratio:.3}");
        assert!(ratio > 0.5, "expected ratio > 0.5 (sanity), got {ratio:.3}");
    }

    #[test]
    fn spectral_clamp_fires_on_large_l_down_jump() {
        // L drops 45 → 9 (Δ = −36). 2-D sweep (2026-05-14) shows the
        // chip clamps both directions when |ΔL| > 25 — the previous
        // "down-jumps don't clamp" hypothesis was based on a different
        // probe schedule (single-frame return-to-baseline) that
        // confounded the direction with the recovery transient.
        let p_high = build_params(0.05, &vec![true; 45], &vec![50f32; 45]);
        let p_low = build_params(0.5, &vec![true; 9], &vec![50f32; 9]);
        let err = FrameErrorContext::default();

        let rms = |pcm: &[i16]| {
            (pcm.iter().map(|&s| f64::from(s).powi(2)).sum::<f64>() / pcm.len() as f64).sqrt()
        };

        let mut state_off = SynthState::new();
        let _ = synthesize_frame_ambe_plus(&p_high, &err, GAMMA_W, &mut state_off);
        let pcm_off = synthesize_frame_ambe_plus(&p_low, &err, GAMMA_W, &mut state_off);

        let mut state_on = SynthState::new();
        state_on.set_chip_compat_spectral_clamp(true);
        let _ = synthesize_frame_ambe_plus(&p_high, &err, GAMMA_W, &mut state_on);
        let pcm_on = synthesize_frame_ambe_plus(&p_low, &err, GAMMA_W, &mut state_on);

        let rms_off = rms(&pcm_off);
        let rms_on = rms(&pcm_on);
        assert!(rms_on < rms_off,
            "clamp must reduce RMS on large L-down-jump: off={rms_off:.1} on={rms_on:.1}");
    }

    #[test]
    fn spectral_clamp_does_not_fire_on_small_l_change() {
        // L jumps 18 → 25 (Δ = +7). Under the new |ΔL| > 25 threshold,
        // this is below the chip-clamp regime and the clamp must NOT
        // fire. Verifies we don't over-clamp natural voice transitions.
        let p_a = build_params(0.2, &vec![true; 18], &vec![50f32; 18]);
        let p_b = build_params(0.15, &vec![true; 25], &vec![50f32; 25]);
        let err = FrameErrorContext::default();

        let mut state_off = SynthState::new();
        let _ = synthesize_frame_ambe_plus(&p_a, &err, GAMMA_W, &mut state_off);
        let pcm_off = synthesize_frame_ambe_plus(&p_b, &err, GAMMA_W, &mut state_off);

        let mut state_on = SynthState::new();
        state_on.set_chip_compat_spectral_clamp(true);
        let _ = synthesize_frame_ambe_plus(&p_a, &err, GAMMA_W, &mut state_on);
        let pcm_on = synthesize_frame_ambe_plus(&p_b, &err, GAMMA_W, &mut state_on);

        assert_eq!(&pcm_off[..], &pcm_on[..],
            "small L change (|ΔL|=7, below threshold of 25) must not trigger clamp");
    }

    #[test]
    fn spectral_clamp_fires_on_high_epsilon_0() {
        // No L change, but ε₀ ≥ 2 → clamp must fire (covers the Repeat
        // case where last_good replay yields ΔL = 0 yet the chip
        // attenuates). Peak is dominated by overlap-add from the prior
        // frame's window tail (which the clamp doesn't reach), so
        // assert on RMS instead — it integrates the whole frame.
        let p = build_params(0.3, &vec![true; 20], &vec![50f32; 20]);
        let err_clean = FrameErrorContext::default();
        let err_e2 = FrameErrorContext { epsilon_0: 2, epsilon_t: 2, ..Default::default() };

        // Warm up state with a clean frame so ΔL=0 on the next call.
        let mut state_off = SynthState::new();
        let _ = synthesize_frame_ambe_plus(&p, &err_clean, GAMMA_W, &mut state_off);
        let pcm_off = synthesize_frame_ambe_plus(&p, &err_e2, GAMMA_W, &mut state_off);

        let mut state_on = SynthState::new();
        state_on.set_chip_compat_spectral_clamp(true);
        let _ = synthesize_frame_ambe_plus(&p, &err_clean, GAMMA_W, &mut state_on);
        let pcm_on = synthesize_frame_ambe_plus(&p, &err_e2, GAMMA_W, &mut state_on);

        let rms = |pcm: &[i16]| {
            (pcm.iter().map(|&s| f64::from(s).powi(2)).sum::<f64>() / pcm.len() as f64).sqrt()
        };
        let rms_off = rms(&pcm_off);
        let rms_on = rms(&pcm_on);
        assert!(rms_on < rms_off, "ε₀≥2 must trigger clamp: off={rms_off:.1} on={rms_on:.1}");
    }

    #[test]
    fn enhance_clamps_w_l_within_brackets() {
        // For an extreme spectral shape we expect some W_l values to
        // hit the clamp (1.2 high, 0.5 low). Non-bypassed M̄_l / M̃_l
        // ratios must lie in [0.5·γ, 1.2·γ] for some shared γ.
        let l = 12usize; // small enough that 8·l > L̃ for all l ≥ 2
        let m_tilde: Vec<f32> = vec![1e-3, 10.0, 0.01, 5.0, 0.05, 3.0, 0.1, 2.0, 0.2, 1.5, 0.3, 1.0];
        let omega_0 = std::f32::consts::PI / 6.0;
        let (m_bar, _) = enhance_spectral_amplitudes(&m_tilde, omega_0, INIT_S_E);
        // All non-bypassed (l ≥ 2 since L=12 means only l=1 is bypassed)
        // ratios M̄/M̃ before γ would be in [0.5, 1.2]. After γ scaling
        // they all share the same γ. So the ratio of any two should be
        // bounded by (1.2/0.5) = 2.4.
        let mut ratios: Vec<f32> = (1..l).map(|i| m_bar[i] / m_tilde[i]).collect();
        ratios.sort_by(|a, b| a.partial_cmp(b).unwrap());
        let span = ratios[ratios.len() - 1] / ratios[0];
        assert!(span <= 2.4 + 1e-4, "ratio span {span} exceeds 2.4");
    }
}