From 587d8bdf5bb7af9d13978539d6d5dc7c926697a5 Mon Sep 17 00:00:00 2001 From: vangel76 Date: Sun, 25 Oct 2026 12:02:37 +0100 Subject: [PATCH 1/3] 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}"); + } } From fd84350122c851d234ff17975348348cf47bcb46 Mon Sep 17 00:00:00 2001 From: vangel76 Date: Sun, 25 Oct 2026 12:03:57 +0100 Subject: [PATCH 2/3] opusenc: halve the a2nlsf bisection steps Each root was refined with 40 bisection steps, which shrinks the bracket to about 3e-15 rad. A Q15 NLSF step is 9.6e-5 rad, so 20 steps (3e-9 rad) are far more than the quantiser can see. The saved cos calls are the remaining bulk of the root search after the grid table. The coded packet streams are byte-identical to the 40-step search on 24 configurations of 90 s: speech-like and music-like input, mono and stereo, SILK NB 8k, MB 12k, WB 16k and 32k, hybrid SWB 32k and FB 64k. Co-Authored-By: Claude Sonnet 5.5 --- crates/opusenc/src/silk/analysis.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/crates/opusenc/src/silk/analysis.rs b/crates/opusenc/src/silk/analysis.rs index 2808b64b1..c84619222 100644 --- a/crates/opusenc/src/silk/analysis.rs +++ b/crates/opusenc/src/silk/analysis.rs @@ -88,8 +88,8 @@ fn cos_series(p: &[f64], w: f64) -> f64 { 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; +/// Bisection steps per root: the bracket shrinks to 3e-9 rad, far below one Q15 NLSF step (9.6e-5 rad). +const BISECT_STEPS: usize = 20; /// `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. From 9a0e19aae73e0d20c03c253fd6cf667d58cf5e37 Mon Sep 17 00:00:00 2001 From: vangel76 Date: Sun, 25 Oct 2026 12:05:19 +0100 Subject: [PATCH 3/3] opusenc: memoise the pitch contour correlations The contour search tries 5 base lags times 34 codebook contours times 4 subframes, 680 normalised correlations, but only about 4 x 30 distinct (subframe, lag) pairs occur. Cache them per call. The values are the same f64 results as before, so the bitstream does not change (cmp on 24 SILK and hybrid configurations of 90 s, mono and stereo, speech and music). Co-Authored-By: Claude Sonnet 5.5 --- crates/opusenc/src/silk/analysis.rs | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/crates/opusenc/src/silk/analysis.rs b/crates/opusenc/src/silk/analysis.rs index c84619222..4f7c7bfab 100644 --- a/crates/opusenc/src/silk/analysis.rs +++ b/crates/opusenc/src/silk/analysis.rs @@ -229,9 +229,20 @@ pub fn pitch_search(x: &[f32], start: usize, len: usize, fs_khz: usize) -> Pitch let score = full[(lag - min_lag) as usize]; // Contour: per-subframe lag offsets around the chosen lag. let cb: &[[i32; 4]] = if fs_khz == 8 { &CB_NB_20 } else { &CB_WB_20 }; + // The contour search revisits the same (subframe, lag) pairs for neighbouring base lags and + // codebook entries, so each correlation is computed once. + let span = (18 * fs - min_lag + 1) as usize; + let memo = std::cell::RefCell::new(vec![f64::NAN; 4 * span]); let sub_score = |k: usize, l: i32| -> f64 { let l = l.clamp(min_lag, 18 * fs); - ncorr(x, start + k * sub..start + (k + 1) * sub, l as usize) + let slot = k * span + (l - min_lag) as usize; + let cached = memo.borrow()[slot]; + if !cached.is_nan() { + return cached; + } + let v = ncorr(x, start + k * sub..start + (k + 1) * sub, l as usize); + memo.borrow_mut()[slot] = v; + v }; let mut best = (f64::MIN, lag, 0usize); for l0 in (lag - 2).max(min_lag)..=(lag + 2).min(max_lag) {