From 587d8bdf5bb7af9d13978539d6d5dc7c926697a5 Mon Sep 17 00:00:00 2001 From: vangel76 Date: Sun, 25 Oct 2026 12:02:37 +0100 Subject: [PATCH] opusenc: tabulate the cosine grid of the a2nlsf root search The LPC to LSF conversion scans a 1024-point grid for sign changes of two cosine series and called libm cos for every series term, about 21,000 calls per channel and frame (30-40 % of SILK and hybrid encode time). Compute cos(m * w) for every grid point and multiple once, with the same expression as before, so the values are bit-identical and the bitstream does not change. Only the bisection steps still call cos. A test compares the new root search against the previous one bit for bit on random symmetric polynomials; a2nlsf rejects orders beyond the table. Co-Authored-By: Claude Sonnet 5.5 --- crates/opusenc/src/silk/analysis.rs | 99 +++++++++++++++++++++++++++-- 1 file changed, 92 insertions(+), 7 deletions(-) diff --git a/crates/opusenc/src/silk/analysis.rs b/crates/opusenc/src/silk/analysis.rs index de510433e..2808b64b1 100644 --- a/crates/opusenc/src/silk/analysis.rs +++ b/crates/opusenc/src/silk/analysis.rs @@ -69,28 +69,55 @@ pub fn levinson(r: &[f64], order: usize) -> Vec { } /// Evaluates the symmetric polynomial `p` (length `d + 1`, `d` even) as the cosine series -/// `e^{jωd/2} P(e^{jω})`. -fn cos_series(p: &[f64], w: f64) -> f64 { +/// `e^{jωd/2} P(e^{jω})`, with `cos_m(m)` giving `cos(m ω)`. +#[inline] +fn cos_series_with(p: &[f64], cos_m: impl Fn(usize) -> f64) -> f64 { let h = (p.len() - 1) / 2; let mut s = p[h]; for k in 0..h { - s += 2.0 * p[k] * ((h - k) as f64 * w).cos(); + s += 2.0 * p[k] * cos_m(h - k); } s } +fn cos_series(p: &[f64], w: f64) -> f64 { + cos_series_with(p, |m| (m as f64 * w).cos()) +} + +/// Number of grid intervals the root search scans over `[0, π]`. +const GRID: usize = 1024; +/// Largest multiple of ω in the series (half the LPC order, at most 16 / 2). +const MAX_H: usize = 8; +/// Bisection steps per root. +const BISECT_STEPS: usize = 40; + +/// `cos(m · π g / GRID)` for every grid point `g` and `m` in `0..=MAX_H`, computed with the same +/// expression as [`cos_series`], so the table values are bit-identical to calling `cos` there. +fn grid_cos() -> &'static [[f64; MAX_H + 1]] { + static TABLE: std::sync::OnceLock> = std::sync::OnceLock::new(); + TABLE.get_or_init(|| { + (0..=GRID) + .map(|g| { + let w = std::f64::consts::PI * g as f64 / GRID as f64; + std::array::from_fn(|m| (m as f64 * w).cos()) + }) + .collect() + }) +} + fn roots(p: &[f64], out: &mut Vec) { - const GRID: usize = 1024; + let grid = grid_cos(); + let at_grid = |g: usize| cos_series_with(p, |m| grid[g][m]); let mut w0 = 0.0; - let mut f0 = cos_series(p, w0); + let mut f0 = at_grid(0); for g in 1..=GRID { let w1 = std::f64::consts::PI * g as f64 / GRID as f64; - let f1 = cos_series(p, w1); + let f1 = at_grid(g); if f0 == 0.0 { out.push(w0); } else if f0 * f1 < 0.0 { let (mut lo, mut hi, mut flo) = (w0, w1, f0); - for _ in 0..40 { + for _ in 0..BISECT_STEPS { let mid = 0.5 * (lo + hi); let fm = cos_series(p, mid); if fm * flo <= 0.0 { @@ -111,6 +138,9 @@ fn roots(p: &[f64], out: &mut Vec) { /// when the root search fails (unstable or degenerate filter). pub fn a2nlsf(a: &[f64]) -> Option> { let d = a.len(); + if d > 2 * MAX_H { + return None; + } let mut c = vec![0f64; d + 2]; c[0] = 1.0; for k in 0..d { @@ -319,4 +349,59 @@ mod tests { } } } + + /// The root search as it was before the grid table: `cos` called for every evaluation. + fn roots_reference(p: &[f64], out: &mut Vec) { + let mut w0 = 0.0; + let mut f0 = cos_series(p, w0); + for g in 1..=GRID { + let w1 = std::f64::consts::PI * g as f64 / GRID as f64; + let f1 = cos_series(p, w1); + if f0 == 0.0 { + out.push(w0); + } else if f0 * f1 < 0.0 { + let (mut lo, mut hi, mut flo) = (w0, w1, f0); + for _ in 0..BISECT_STEPS { + let mid = 0.5 * (lo + hi); + let fm = cos_series(p, mid); + if fm * flo <= 0.0 { + hi = mid; + } else { + lo = mid; + flo = fm; + } + } + out.push(0.5 * (lo + hi)); + } + w0 = w1; + f0 = f1; + } + } + + #[test] + fn grid_table_roots_are_bit_identical_to_direct_cos() { + let mut seed = 7u64; + let mut rnd = move || { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + (seed >> 11) as f64 / (1u64 << 53) as f64 * 2.0 - 1.0 + }; + let mut found = 0; + for len in [11usize, 17] { + for _ in 0..300 { + let mut p = vec![0f64; len]; + for k in 0..len / 2 + 1 { + p[k] = rnd(); + p[len - 1 - k] = p[k]; + } + let (mut a, mut b) = (Vec::new(), Vec::new()); + roots(&p, &mut a); + roots_reference(&p, &mut b); + assert_eq!(a.iter().map(|v| v.to_bits()).collect::>(), b.iter().map(|v| v.to_bits()).collect::>()); + found += a.len(); + } + } + assert!(found > 1000, "{found}"); + } }