diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index 0e85749f..dfe5e003 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -77,3 +77,4 @@ and the exact transfer revision are recorded in each vendored crate's SOURCE.jso | `hayro-interpret 0.7.0 / assets/CGATS001Compat-v2-micro.icc` | Clinton Ingram (saucecontrol/Compact-ICC-Profiles) | https://github.com/storytold/pdfcraft/blob/7437bde378f43856eada40af25b983812e615e8a/vendor/hayro-interpret/assets/CGATS001Compat-v2-micro.icc | CC0-1.0 | `73e1ba37d2bad5bab2a964f40a9eed96209666efc067c3322626214bbef234a0` | Default CMYK colour conversion (pinned hayro-interpret) | | `hayro-cmap 0.1.0 / assets/cmaps.brotli` | Adobe | https://github.com/storytold/pdfcraft/blob/a2ec27b4775f53c40cad3456936d9406611e3a8d/vendor/hayro-cmap/assets/cmaps.brotli | BSD-3-Clause | `ddf37af521d904d5646c2fc37b2b7e912dc1b56f0d210190b61f598154758e46` | Character-code mapping tables required to read CJK PDFs. Non-visual dependency data; existing embedding unchanged | | `hayro-interpret 0.7.0 / src/font/generated/metrics.rs` | Adobe (Core 14 AFM metrics, PDF specification encodings); tables by The Hayro Authors | https://github.com/storytold/pdfcraft/blob/7437bde378f43856eada40af25b983812e615e8a/vendor/hayro-interpret/src/font/generated/metrics.rs | MIT OR Apache-2.0 | `2ec9e608573eda459de1e1d9bda3bb494315e36bee885acc88912a11fe0caa41` | Glyph widths and encodings for non-embedded standard fonts. Non-visual dependency data; existing embedding unchanged | +| `crates/audio/src/ac3/golden_frames.bin` | FilmCraft contributors (original work: 12 AC-3 and E-AC-3 syncframes that ffmpeg's `ac3` and `eac3` encoders made from tones and seeded noise synthesised by ffmpeg's lavfi generators; no third-party or recorded audio) | [FilmCraft `crates/ac3`](https://github.com/storytold/filmcraft/tree/1c9c2e82/crates/ac3/src) | MIT OR Apache-2.0 | diff --git a/crates/audio/Cargo.toml b/crates/audio/Cargo.toml index 9c6c1d22..53949299 100644 --- a/crates/audio/Cargo.toml +++ b/crates/audio/Cargo.toml @@ -34,13 +34,13 @@ flac-encode = ["craft-media/bitstream"] aac-decode = ["craft-media/bitstream"] aac-encode = ["craft-media/bitstream"] opus-decode = [] -ac3-decode = [] +ac3-decode = ["craft-media/bitstream"] mpa-decode = ["craft-media/bitstream"] # ALAC packet decoder and encoder (`alac`), ported from Apple's reference implementation. alac-decode = [] alac-encode = [] -ac3-encode = [] -eac3-encode = [] +ac3-encode = ["craft-media/bitstream"] +eac3-encode = ["craft-media/bitstream"] # Multi-resolution waveform overviews (`peaks`). peaks = [] # Encode AAC channels in parallel (with `aac-encode`). diff --git a/crates/audio/README.md b/crates/audio/README.md index 1b67e4ea..e05ae3e7 100644 --- a/crates/audio/README.md +++ b/crates/audio/README.md @@ -36,8 +36,8 @@ No default features. The shared types (`AudioBuffer`, `AudioInfo`, `ChannelLayou | `alac-decode` | `alac::Decoder` | none | | `alac-encode` | `alac::Encoder` | none | | `opus-decode` | `opus::Decoder`/`StreamDecoder` | none | -| `ac3-decode` | `ac3::Decoder`, `ac3::ProgrammeDecoder` | none | -| `ac3-encode`, `eac3-encode` | `ac3::encoder::Encoder` | none | +| `ac3-decode` | `ac3::Decoder`, `ac3::ProgrammeDecoder` | craft-media `bitstream` | +| `ac3-encode`, `eac3-encode` | `ac3::encoder::Encoder` | craft-media `bitstream` | | `mpa-decode` | `mpa::Decoder` | craft-media `bitstream` | | `peaks` | `peaks::Peaks` | none | | `rayon` | parallel AAC encoding (with `aac-encode`) | `rayon` | diff --git a/crates/audio/src/aac/decoder.rs b/crates/audio/src/aac/decoder.rs index d072796e..8c4d4130 100644 --- a/crates/audio/src/aac/decoder.rs +++ b/crates/audio/src/aac/decoder.rs @@ -5,8 +5,8 @@ use craft_media::bitstream::BitReader; use crate::aac::config::{AudioSpecificConfig, ElementType, ProgramConfig}; use crate::aac::huffman::{self, INTENSITY_HCB, INTENSITY_HCB2, NOISE_HCB, ZERO_HCB}; use crate::aac::ics::{IcsInfo, WindowSequence}; -use crate::aac::mdct::{WindowShape, mdct_long, mdct_short, window}; -use crate::aac::tables::{pow43, sf_gain}; +use crate::aac::mdct::{MdctScratch, WindowShape, mdct_long, mdct_short, window}; +use crate::aac::tables::{pow43, pow43_table, sf_gain}; use crate::aac::tns::{self, TnsData}; use crate::aac::{Error, Result}; @@ -45,6 +45,24 @@ struct ChannelState { prev_shape: WindowShape, } +/// Buffers reused by every frame so the steady state does not allocate per channel. +#[derive(Clone)] +struct Scratch { + /// Quantised values of the channel being parsed. + q: Vec, + /// Windowed output of the current block (2048) and the raw inverse transform (2048 / 256). + buf: Vec, + y: Vec, + long: MdctScratch, + short: MdctScratch, +} + +impl Scratch { + fn new() -> Scratch { + Scratch { q: vec![0; 1024], buf: vec![0.0; 2048], y: vec![0.0; 2048], long: mdct_long().scratch(), short: mdct_short().scratch() } + } +} + /// AAC-LC decoder. #[derive(Clone)] pub struct Decoder { @@ -54,6 +72,7 @@ pub struct Decoder { offsets: Vec, channels: Vec, noise_state: u32, + scratch: Scratch, /// A fill element carried SBR data (implicitly signalled HE-AAC). sbr_seen: bool, } @@ -83,7 +102,15 @@ impl Decoder { if layout.is_empty() && asc.channel_config != 0 { return Err(Error::InvalidConfig("unsupported channel configuration")); } - let mut d = Decoder { asc, layout: Vec::new(), offsets: Vec::new(), channels: Vec::new(), noise_state: NOISE_SEED, sbr_seen: false }; + let mut d = Decoder { + asc, + layout: Vec::new(), + offsets: Vec::new(), + channels: Vec::new(), + noise_state: NOISE_SEED, + scratch: Scratch::new(), + sbr_seen: false, + }; d.set_layout(layout); Ok(d) } @@ -406,7 +433,9 @@ impl Decoder { } // spectral_data - let mut q = vec![0i32; 1024]; + let mut q = std::mem::take(&mut self.scratch.q); + q.clear(); + q.resize(1024, 0); let mut vals = [0i32; 4]; for g in 0..info.num_groups { let w0 = info.group_start(g); @@ -438,6 +467,7 @@ impl Decoder { } // dequantise + let t43 = pow43_table(); for g in 0..info.num_groups { let w0 = info.group_start(g); for sfb in 0..info.max_sfb { @@ -447,14 +477,16 @@ impl Decoder { } let gain = sf_gain(ics.sf[g][sfb]); for w in w0..w0 + info.group_len[g] as usize { - for k in swb[sfb] as usize..swb[sfb + 1] as usize { - let v = q[w * 128 + k]; - let m = pow43(v.unsigned_abs()) * gain; - ics.spec[w * 128 + k] = if v < 0 { -m } else { m }; + let r = w * 128 + swb[sfb] as usize..w * 128 + swb[sfb + 1] as usize; + for (o, &v) in ics.spec[r.clone()].iter_mut().zip(&q[r]) { + let a = v.unsigned_abs(); + let m = t43.get(a as usize).copied().unwrap_or_else(|| pow43(a)) * gain; + *o = if v < 0 { -m } else { m }; } } } } + self.scratch.q = q; Ok(ics) } @@ -504,23 +536,24 @@ impl Decoder { if let Some(t) = &ics.tns { tns::apply_decoder(&mut ics.spec, &ics.info, self.sf_index(), t); } - synth(&mut self.channels[ch], &ics.info, &ics.spec) + synth(&mut self.channels[ch], &mut self.scratch, &ics.info, &ics.spec) } } /// Inverse MDCT, windowing and overlap-add (§4.6.18), output scaled to ±1.0 full scale. -fn synth(state: &mut ChannelState, info: &IcsInfo, spec: &[f32]) -> Vec { - let mut buf = vec![0f32; 2048]; +fn synth(state: &mut ChannelState, sc: &mut Scratch, info: &IcsInfo, spec: &[f32]) -> Vec { + let Scratch { buf, y, long, short, .. } = sc; + buf.fill(0.0); let prev = state.prev_shape; let cur = info.window_shape; let norm = 1.0 / 32768.0; match info.window_sequence { WindowSequence::EightShort => { let t = mdct_short(); - let mut y = vec![0f32; 256]; + let y = &mut y[..256]; let wc = window(cur, false); for w in 0..8 { - t.inverse(&spec[w * 128..(w + 1) * 128], &mut y); + t.inverse_with(&spec[w * 128..(w + 1) * 128], y, short); let wl = window(if w == 0 { prev } else { cur }, false); let off = 448 + 128 * w; let s = norm / 128.0; @@ -531,8 +564,7 @@ fn synth(state: &mut ChannelState, info: &IcsInfo, spec: &[f32]) -> Vec { } } seq => { - let mut y = vec![0f32; 2048]; - mdct_long().inverse(&spec[..1024], &mut y); + mdct_long().inverse_with(&spec[..1024], y, long); let s = norm / 1024.0; match seq { WindowSequence::LongStop => { @@ -711,7 +743,7 @@ mod tests { fn non_finite_spectra_do_not_reach_the_output_or_the_overlap() { let mut state = ChannelState { overlap: vec![0.0; 1024], prev_shape: WindowShape::Sine }; let info = IcsInfo::parse(&mut BitReader::new(&[0; 2]), 4).unwrap(); - let out = synth(&mut state, &info, &vec![f32::INFINITY; 1024]); + let out = synth(&mut state, &mut Scratch::new(), &info, &vec![f32::INFINITY; 1024]); assert!(out.iter().chain(&state.overlap).all(|s| s.is_finite())); } } diff --git a/crates/audio/src/aac/encoder/mod.rs b/crates/audio/src/aac/encoder/mod.rs index 5230df5b..5a0ea5ac 100644 --- a/crates/audio/src/aac/encoder/mod.rs +++ b/crates/audio/src/aac/encoder/mod.rs @@ -25,7 +25,7 @@ use crate::aac::tables::{sample_rate_index, swb_offsets_long, swb_offsets_short, use crate::aac::tns::{self, TnsData, TnsFilter}; use crate::aac::{Error, FRAME_LEN, Result}; use psy::Psy; -use quant::{ChannelInput, QuantChannel, quantise_channel, write_ics}; +use quant::{ChannelInput, QuantCache, QuantChannel, quantise_channel, write_ics}; /// Rate control mode. #[derive(Debug, Clone, Copy, PartialEq)] @@ -153,6 +153,8 @@ struct ElementState { struct ElementFrame { info: IcsInfo, chans: Vec, + /// Rate-loop memo per channel, filled while the frame is quantised. + caches: Vec, ms_mask: Option<(u32, Vec)>, } @@ -658,14 +660,15 @@ impl Encoder { }; ms_mask = Some((present, used)); } - let chans = specs.iter().zip(thrs.iter()).zip(tnss).map(|((s, th), tn)| ChannelInput::new(info, swb, s, th, tn)).collect(); + let chans: Vec = specs.iter().zip(thrs.iter()).zip(tnss).map(|((s, th), tn)| ChannelInput::new(info, swb, s, th, tn)).collect(); + let caches = chans.iter().map(QuantCache::new).collect(); let e = &mut self.elements[ei]; e.prev_seq = seq; e.short_next = Some(short_next); if short { self.stats.short_blocks += 1; } - (ElementFrame { info, chans, ms_mask }, pe) + (ElementFrame { info, chans, caches, ms_mask }, pe) } fn element_overhead(&self, ei: usize, f: &ElementFrame) -> u32 { @@ -682,15 +685,21 @@ impl Encoder { } } - fn quantise_all(frames: &[ElementFrame], lambda: f32) -> Vec> { + fn quantise_all(frames: &mut [ElementFrame], lambda: f32) -> Vec> { #[cfg(feature = "rayon")] { use rayon::prelude::*; - frames.par_iter().map(|f| f.chans.par_iter().map(|c| quantise_channel(c, lambda)).collect()).collect() + frames + .par_iter_mut() + .map(|f| f.chans.par_iter().zip(f.caches.par_iter_mut()).map(|(c, cache)| quantise_channel(c, cache, lambda)).collect()) + .collect() } #[cfg(not(feature = "rayon"))] { - frames.iter().map(|f| f.chans.iter().map(|c| quantise_channel(c, lambda)).collect()).collect() + frames + .iter_mut() + .map(|f| f.chans.iter().zip(f.caches.iter_mut()).map(|(c, cache)| quantise_channel(c, cache, lambda)).collect()) + .collect() } } @@ -743,16 +752,16 @@ impl Encoder { BitrateMode::Vbr(q) => { let d = q.clamp(0.5, 6.0) - 1.0; let mut lambda = 10f32.powf(-0.2 - 1.2 * d + 0.09 * d * d); - let mut res = Self::quantise_all(&frames, lambda); + let mut res = Self::quantise_all(&mut frames, lambda); while total_bits(&res) > target && lambda < 1e9 { lambda *= 2.0; - res = Self::quantise_all(&frames, lambda); + res = Self::quantise_all(&mut frames, lambda); } res } BitrateMode::Cbr(_) => { let mut lam = self.lambda.clamp(1e-6, 1e9); - let mut cur = Self::quantise_all(&frames, lam); + let mut cur = Self::quantise_all(&mut frames, lam); let (mut lo, mut hi, mut best); if total_bits(&cur) <= target { hi = lam; @@ -763,7 +772,7 @@ impl Encoder { if lo < 1e-6 { break; } - cur = Self::quantise_all(&frames, lo); + cur = Self::quantise_all(&mut frames, lo); if total_bits(&cur) > target { break; } @@ -774,7 +783,7 @@ impl Encoder { lo = lam; loop { lam *= 4.0; - cur = Self::quantise_all(&frames, lam); + cur = Self::quantise_all(&mut frames, lam); if total_bits(&cur) <= target || lam > 1e9 { break; } @@ -788,7 +797,7 @@ impl Encoder { break; } let mid = (lo * hi).sqrt(); - let r = Self::quantise_all(&frames, mid); + let r = Self::quantise_all(&mut frames, mid); if total_bits(&r) <= target { hi = mid; best = r; diff --git a/crates/audio/src/aac/encoder/quant.rs b/crates/audio/src/aac/encoder/quant.rs index 5ca742d2..d19b68f9 100644 --- a/crates/audio/src/aac/encoder/quant.rs +++ b/crates/audio/src/aac/encoder/quant.rs @@ -8,7 +8,7 @@ use craft_media::bitstream::BitWriter; use crate::aac::huffman::{self, ZERO_HCB}; use crate::aac::ics::IcsInfo; -use crate::aac::tables::pow43; +use crate::aac::tables::pow43_table; use crate::aac::tns::TnsData; /// Magic rounding offset for AAC quantisation (`0.5 - 0.0946`), minimising the error for the @@ -67,6 +67,42 @@ impl ChannelInput { } } +/// Scalefactor-keyed results of one band, valid until the band is quantised at another scalefactor. +#[derive(Clone, Copy)] +struct BandMemo { + /// Scalefactor the cached values belong to (`NO_SF` when empty). + sf: i32, + maxq: i32, + /// Spectral bits per codebook plus one for an all-zero band (see `band_costs`), once computed. + cost: Option<[u32; 12]>, +} + +const NO_SF: i32 = i32::MIN; + +/// Work that depends only on (band, scalefactor), not on λ, kept across the rate-loop iterations of +/// one frame. One cache per channel; callers hold it mutably, so parallel channels share nothing. +pub struct QuantCache { + /// Distortion per (band, scalefactor), NaN while unknown. + dist: Vec, + bands: Vec, + /// Quantised values of every band at its memoised scalefactor. + q: Vec, + /// Cost row of a band that quantises to zeros. + zero_cost: Vec>, +} + +impl QuantCache { + pub fn new(ch: &ChannelInput) -> QuantCache { + let nb = ch.bands.len(); + QuantCache { + dist: vec![f32::NAN; nb * 256], + bands: vec![BandMemo { sf: NO_SF, maxq: 0, cost: None }; nb], + q: vec![0; ch.x.len()], + zero_cost: vec![None; nb], + } + } +} + /// A quantised channel ready to be written. #[derive(Clone, Default)] pub struct QuantChannel { @@ -90,46 +126,54 @@ fn gain(sf: i32) -> f32 { (0.25 * (sf as f32 - 100.0)).exp2() } -/// Quantise a slice at `sf`, returning the distortion; writes the (signed) values to `out` if given. -fn quantise(x: &[f32], xabs: &[f32], x34: &[f32], sf: i32, out: Option<&mut [i32]>) -> f32 { +/// Distortion of a slice quantised at `sf`. Kept out of line so the accumulator stays in a register. +#[inline(never)] +fn distortion(xabs: &[f32], x34: &[f32], sf: i32) -> f32 { let g34 = gain34(sf); let g = gain(sf); + let t = pow43_table(); let mut d = 0f32; - match out { - Some(out) => { - for i in 0..x.len() { - let q = ((x34[i] * g34 + ROUND) as i32).min(MAX_Q); - let e = xabs[i] - pow43(q as u32) * g; - d += e * e; - out[i] = if x[i] < 0.0 { -q } else { q }; - } - } - None => { - for i in 0..x.len() { - let q = ((x34[i] * g34 + ROUND) as i32).min(MAX_Q); - let e = xabs[i] - pow43(q as u32) * g; - d += e * e; - } - } + for (&a, &b) in xabs.iter().zip(x34) { + let q = ((b * g34 + ROUND) as i32).min(MAX_Q); + let e = a - t[(q as usize) & 8191] * g; + d += e * e; } d } +/// Quantise a slice at `sf`, writing the (signed) values to `out`. +fn quantise(x: &[f32], sf: i32, x34: &[f32], out: &mut [i32]) { + let g34 = gain34(sf); + for ((o, &xv), &b) in out.iter_mut().zip(x).zip(x34) { + let q = ((b * g34 + ROUND) as i32).min(MAX_Q); + *o = if xv < 0.0 { -q } else { q }; + } +} + /// Largest scalefactor whose distortion stays within `allowed`. -fn search_sf(ch: &ChannelInput, b: &Band, allowed: f32) -> i32 { +fn search_sf(ch: &ChannelInput, bi: usize, dist: &mut [f32], allowed: f32) -> i32 { + let b = &ch.bands[bi]; let r = b.start..b.end; - let (x, xa, x34) = (&ch.x[r.clone()], &ch.xabs[r.clone()], &ch.x34[r]); + let (xa, x34) = (&ch.xabs[r.clone()], &ch.x34[r]); + let dist = &mut dist[bi * 256..bi * 256 + 256]; + let mut dist_at = |sf: i32| { + let d = &mut dist[sf as usize]; + if d.is_nan() { + *d = distortion(xa, x34, sf); + } + *d + }; let l2 = b.max34.log2(); let lo = ((100.0 - ((MAX_Q as f32 + 0.5).log2() - l2) / 0.1875).ceil() as i32).clamp(0, 255); let zero = ((100.0 - ((1.0 - ROUND).log2() - l2) / 0.1875).floor() as i32 + 1).clamp(0, 256); let mut a = lo; let mut hi = (zero - 1).clamp(lo, 255); - if quantise(x, xa, x34, a, None) > allowed { + if dist_at(a) > allowed { return a; } while a < hi { let mid = (a + hi + 1) / 2; - if quantise(x, xa, x34, mid, None) <= allowed { + if dist_at(mid) <= allowed { a = mid; } else { hi = mid - 1; @@ -151,9 +195,193 @@ fn candidates(m: i32) -> &'static [u8] { } } +/// Bring the scalefactor differences between coded bands (`sf >= 0`) within ±60 by lowering the +/// offenders, sweeping up and down until stable. +fn limit_sf_steps(sf: &mut [i32]) { + let coded: Vec = (0..sf.len()).filter(|&i| sf[i] >= 0).collect(); + for _ in 0..8 { + let mut changed = false; + for w in coded.windows(2) { + if sf[w[1]] > sf[w[0]] + 60 { + sf[w[1]] = sf[w[0]] + 60; + changed = true; + } + } + for w in coded.windows(2).rev() { + if sf[w[0]] > sf[w[1]] + 60 { + sf[w[0]] = sf[w[1]] + 60; + changed = true; + } + } + if !changed { + break; + } + } +} + +/// Quantise every band at its scalefactor into `q`, recording each band's largest magnitude in +/// `maxq`; a band already quantised at that scalefactor is taken from `cache`. A band that quantises +/// to all zeros leaves the chain (`sf = -1`); returns whether any did. +fn quantise_bands(ch: &ChannelInput, cache: &mut QuantCache, sf: &mut [i32], q: &mut [i32], maxq: &mut [i32]) -> bool { + let mut removed = false; + for (i, b) in ch.bands.iter().enumerate() { + let r = b.start..b.end; + if sf[i] < 0 { + q[r].fill(0); + maxq[i] = 0; + continue; + } + let memo = &mut cache.bands[i]; + if memo.sf != sf[i] { + quantise(&ch.x[r.clone()], sf[i], &ch.x34[r.clone()], &mut cache.q[r.clone()]); + *memo = BandMemo { sf: sf[i], maxq: cache.q[r.clone()].iter().map(|v| v.abs()).max().unwrap_or(0), cost: None }; + } + q[r.clone()].copy_from_slice(&cache.q[r]); + maxq[i] = memo.maxq; + if maxq[i] == 0 { + sf[i] = -1; + removed = true; + } + } + removed +} + +/// Cost sentinel for a codebook that cannot code a band. +const INF: u32 = u32::MAX / 4; + +/// Bits per candidate codebook for a band holding `vals` (largest magnitude `maxq`); `INF` +/// marks codebooks that cannot code it. An all-zero band inside a coded section also needs a (zero) +/// scalefactor difference, hence the extra bit. +fn band_costs(vals: &[i32], maxq: i32) -> [u32; 12] { + let mut cost = [INF; 12]; + let extra = (maxq == 0) as u32; + for &c in candidates(maxq) { + match c { + 0 => cost[0] = 0, + // odd codebook below 11 and its successor share everything but the length table + c if c % 2 == 1 && c < huffman::ESC_HCB => { + let (a, b) = huffman::count_bits_pair(c, vals); + cost[c as usize] = a + extra; + cost[c as usize + 1] = b + extra; + } + c if c % 2 == 0 => {} // filled in by the pair above + c => cost[c as usize] = huffman::count_bits(c, vals) + extra, + } + } + cost +} + +/// Bits to code each band of a group with each codebook (`INF` when the codebook cannot), from the +/// per-band memo. +fn codebook_costs(ch: &ChannelInput, cache: &mut QuantCache, sf: &[i32], q: &[i32], maxq: &[i32], base: usize, n: usize) -> Vec<[u32; 12]> { + let mut cost = vec![[INF; 12]; n]; + for i in 0..n { + let bi = base + i; + let b = &ch.bands[bi]; + let vals = &q[b.start..b.end]; + // A band that ends up uncoded (`sf < 0`) holds zeros, whatever its memo says. + let row = if sf[bi] >= 0 { &mut cache.bands[bi].cost } else { &mut cache.zero_cost[bi] }; + cost[i] = *row.get_or_insert_with(|| band_costs(vals, maxq[bi])); + } + cost +} + +/// Dynamic programming over the codebook of each band of a group, charging `new_sec` bits for every +/// section start. Returns the chosen codebook per band. +fn best_codebooks(cost: &[[u32; 12]], new_sec: u32) -> Vec { + let n = cost.len(); + let mut best = [INF; 12]; + let mut back = vec![[0u8; 12]; n]; + for c in 0..12 { + if cost[0][c] < INF { + best[c] = new_sec + cost[0][c]; + } + } + for i in 1..n { + let (mut min_c, mut min_v) = (0usize, INF); + for c in 0..12 { + if best[c] < min_v { + min_v = best[c]; + min_c = c; + } + } + let mut nb_ = [INF; 12]; + for c in 0..12 { + if cost[i][c] >= INF { + continue; + } + let stay = best[c]; + let switch = min_v + new_sec; + if stay <= switch { + nb_[c] = stay + cost[i][c]; + back[i][c] = c as u8; + } else { + nb_[c] = switch + cost[i][c]; + back[i][c] = min_c as u8; + } + } + best = nb_; + } + let mut cbs = vec![0u8; n]; + let mut c = (0..12).min_by_key(|&c| best[c]).unwrap_or(0); + for i in (0..n).rev() { + cbs[i] = c as u8; + c = back[i][c] as usize; + } + cbs +} + +/// Run-length encode a group's codebooks into `(codebook, length)` sections, and sum the spectral +/// bits of its non-zero sections from the cost rows. `base` is the index of the group's first band. +fn group_sections(cost: &[[u32; 12]], maxq: &[i32], cbs: &[u8], base: usize) -> (Vec<(u8, usize)>, u32) { + let mut secs: Vec<(u8, usize)> = Vec::new(); + let mut spectral_bits = 0u32; + for (i, &c) in cbs.iter().enumerate() { + match secs.last_mut() { + Some((lc, len)) if *lc == c => *len += 1, + _ => secs.push((c, 1)), + } + if c != ZERO_HCB { + // the cost row carries one extra bit for an all-zero band (its scalefactor difference) + spectral_bits += cost[i][c as usize] - (maxq[base + i] == 0) as u32; + } + } + (secs, spectral_bits) +} + +/// The sections of each window group: `(codebook, bands)` runs. +type GroupSections = Vec>; + +/// Section selection: a codebook per band, the sections of each group and the bits they cost +/// `(codebooks, sections, spectral_bits, section_bits)`. +fn select_sections(ch: &ChannelInput, cache: &mut QuantCache, sf: &[i32], q: &[i32], maxq: &[i32]) -> (Vec, GroupSections, u32, u32) { + let info = &ch.info; + let (sect_bits, esc) = if info.is_short() { (3u32, 7usize) } else { (5u32, 31usize) }; + let mut cb = vec![ZERO_HCB; ch.bands.len()]; + let mut sections = Vec::with_capacity(info.num_groups); + let mut spectral_bits = 0u32; + let mut section_bits = 0u32; + let n = info.max_sfb; + for g in 0..info.num_groups { + let base = g * n; + if n == 0 { + sections.push(Vec::new()); + continue; + } + let cost = codebook_costs(ch, cache, sf, q, maxq, base, n); + cb[base..base + n].copy_from_slice(&best_codebooks(&cost, 4 + sect_bits)); + let (secs, bits) = group_sections(&cost, maxq, &cb[base..base + n], base); + spectral_bits += bits; + for &(_, len) in &secs { + section_bits += 4 + sect_bits * (len / esc + 1) as u32; + } + sections.push(secs); + } + (cb, sections, spectral_bits, section_bits) +} + /// Quantise one channel at rate-distortion multiplier `lambda` and select sections. -#[allow(clippy::cognitive_complexity, reason = "the AAC rate/distortion loop, ported as-is from FilmCraft")] -pub fn quantise_channel(ch: &ChannelInput, lambda: f32) -> QuantChannel { +pub fn quantise_channel(ch: &ChannelInput, cache: &mut QuantCache, lambda: f32) -> QuantChannel { let nb = ch.bands.len(); let mut sf = vec![-1i32; nb]; for (i, b) in ch.bands.iter().enumerate() { @@ -161,139 +389,20 @@ pub fn quantise_channel(ch: &ChannelInput, lambda: f32) -> QuantChannel { if b.max34 <= 0.0 || b.en <= allowed { continue; } - sf[i] = search_sf(ch, b, allowed); + sf[i] = search_sf(ch, i, &mut cache.dist, allowed); } // Scalefactor differences between coded bands must stay within ±60: lower (refine) offenders, // then quantise; repeat if a band unexpectedly quantised to all zeros (it leaves the chain). let mut q = vec![0i32; ch.x.len()]; let mut maxq = vec![0i32; nb]; loop { - for _ in 0..8 { - let mut changed = false; - let coded: Vec = (0..nb).filter(|&i| sf[i] >= 0).collect(); - for w in coded.windows(2) { - if sf[w[1]] > sf[w[0]] + 60 { - sf[w[1]] = sf[w[0]] + 60; - changed = true; - } - } - for w in coded.windows(2).rev() { - if sf[w[0]] > sf[w[1]] + 60 { - sf[w[0]] = sf[w[1]] + 60; - changed = true; - } - } - if !changed { - break; - } - } - let mut removed = false; - for (i, b) in ch.bands.iter().enumerate() { - let r = b.start..b.end; - if sf[i] < 0 { - q[r].fill(0); - maxq[i] = 0; - continue; - } - quantise(&ch.x[r.clone()], &ch.xabs[r.clone()], &ch.x34[r.clone()], sf[i], Some(&mut q[r.clone()])); - maxq[i] = q[r].iter().map(|v| v.abs()).max().unwrap_or(0); - if maxq[i] == 0 { - sf[i] = -1; - removed = true; - } - } - if !removed { + limit_sf_steps(&mut sf); + if !quantise_bands(ch, cache, &mut sf, &mut q, &mut maxq) { break; } } - // Section selection: dynamic programming per group over codebooks. - let info = &ch.info; - let (sect_bits, esc) = if info.is_short() { (3u32, 7usize) } else { (5u32, 31usize) }; - let mut cb = vec![ZERO_HCB; nb]; - let mut sections = Vec::with_capacity(info.num_groups); - let mut spectral_bits = 0u32; - let mut section_bits = 0u32; - let per_group = info.max_sfb; - for g in 0..info.num_groups { - let base = g * per_group; - let n = per_group; - if n == 0 { - sections.push(Vec::new()); - continue; - } - // cost[i][c] - const INF: u32 = u32::MAX / 4; - let mut cost = vec![[INF; 12]; n]; - for i in 0..n { - let bi = base + i; - let b = &ch.bands[bi]; - let vals = &q[b.start..b.end]; - for &c in candidates(maxq[bi]) { - cost[i][c as usize] = if c == 0 { - 0 - } else { - // all-zero band inside a coded section also needs a (zero) scalefactor difference - huffman::count_bits(c, vals) + if maxq[bi] == 0 { 1 } else { 0 } - }; - } - } - let new_sec = 4 + sect_bits; - let mut best = [INF; 12]; - let mut back = vec![[0u8; 12]; n]; - for c in 0..12 { - if cost[0][c] < INF { - best[c] = new_sec + cost[0][c]; - } - } - for i in 1..n { - let (mut min_c, mut min_v) = (0usize, INF); - for c in 0..12 { - if best[c] < min_v { - min_v = best[c]; - min_c = c; - } - } - let mut nb_ = [INF; 12]; - for c in 0..12 { - if cost[i][c] >= INF { - continue; - } - let stay = best[c]; - let switch = min_v + new_sec; - if stay <= switch { - nb_[c] = stay + cost[i][c]; - back[i][c] = c as u8; - } else { - nb_[c] = switch + cost[i][c]; - back[i][c] = min_c as u8; - } - } - best = nb_; - } - let mut c = (0..12).min_by_key(|&c| best[c]).unwrap_or(0); - for i in (0..n).rev() { - cb[base + i] = c as u8; - c = back[i][c] as usize; - } - // run-length sections and exact costs - let mut secs: Vec<(u8, usize)> = Vec::new(); - for i in 0..n { - let c = cb[base + i]; - match secs.last_mut() { - Some((lc, len)) if *lc == c => *len += 1, - _ => secs.push((c, 1)), - } - if c != ZERO_HCB { - let b = &ch.bands[base + i]; - spectral_bits += huffman::count_bits(c, &q[b.start..b.end]); - } - } - for &(_, len) in &secs { - section_bits += 4 + sect_bits * (len / esc + 1) as u32; - } - sections.push(secs); - } + let (cb, sections, spectral_bits, section_bits) = select_sections(ch, cache, &sf, &q, &maxq); // Scalefactors: global gain is the first coded band's value; bands coded with a non-zero // codebook but all-zero values reuse the previous value (difference 0). @@ -310,7 +419,7 @@ pub fn quantise_channel(ch: &ChannelInput, lambda: f32) -> QuantChannel { sf_bits += huffman::sf_bits(sf[i] - prev); prev = sf[i]; } - let tns_bits = 1 + ch.tns.as_ref().map(|t| t.bit_count(info)).unwrap_or(0); + let tns_bits = 1 + ch.tns.as_ref().map(|t| t.bit_count(&ch.info)).unwrap_or(0); let bits = 8 + section_bits + sf_bits + 1 + tns_bits + 1 + spectral_bits; QuantChannel { q, sf, cb, sections, global_gain: first, bits } } @@ -374,10 +483,31 @@ mod tests { let thr = vec![vec![1000.0f32; 64]]; let ch = ChannelInput::new(info, swb, &spec, &thr, None); for lambda in [0.01f32, 1.0, 100.0, 1e5] { - let qc = quantise_channel(&ch, lambda); + let mut cache = QuantCache::new(&ch); + let qc = quantise_channel(&ch, &mut cache, lambda); let mut bw = BitWriter::new(); write_ics(&mut bw, &ch, &qc, false); assert_eq!(bw.bit_len() as u32, qc.bits, "lambda {lambda}"); } } + + #[test] + fn reused_cache_gives_the_same_result_as_a_fresh_one() { + let swb = swb_offsets_long(4); + let info = IcsInfo { window_sequence: WindowSequence::OnlyLong, window_shape: WindowShape::Sine, max_sfb: 40, ..Default::default() }; + let spec: Vec = (0..1024).map(|i| ((i * 7919 % 1000) as f32 - 500.0) * (1000.0 / (i as f32 + 10.0))).collect(); + let thr = vec![vec![1000.0f32; 64]]; + let ch = ChannelInput::new(info, swb, &spec, &thr, None); + let mut shared = QuantCache::new(&ch); + // up and down, like the rate loop, with repeats + for lambda in [1.0f32, 4.0, 16.0, 0.25, 1.0, 100.0, 3.0, 3.0, 1e5, 0.01, 1.0] { + let a = quantise_channel(&ch, &mut shared, lambda); + let b = quantise_channel(&ch, &mut QuantCache::new(&ch), lambda); + assert_eq!( + (&a.q, &a.sf, &a.cb, &a.sections, a.global_gain, a.bits), + (&b.q, &b.sf, &b.cb, &b.sections, b.global_gain, b.bits), + "lambda {lambda}" + ); + } + } } diff --git a/crates/audio/src/aac/huffman.rs b/crates/audio/src/aac/huffman.rs index 1de2bb19..7d41a2c7 100644 --- a/crates/audio/src/aac/huffman.rs +++ b/crates/audio/src/aac/huffman.rs @@ -53,60 +53,92 @@ pub fn book(cb: u8) -> &'static Book { &BOOKS[(cb as usize).clamp(1, 11) - 1] } -const LEAF: u32 = 1 << 31; +/// Bits resolved by the first table level. +const L1: u32 = 9; +/// Flag of a first-level entry that points to a sub-table. +const SUB: u32 = 1 << 31; -/// Binary decoding tree for one codebook. -struct Tree { - nodes: Vec<[u32; 2]>, +/// Two-level decoding table for one prefix code. A leaf entry is `len << 24 | symbol`; a first-level +/// entry with `SUB` set is `SUB | offset << 5 | width` and selects a sub-table of `1 << width` +/// leaves indexed by the next `width` bits. An entry of 0 is a bit pattern that no codeword starts +/// with (the codebooks are complete, so this does not occur). +struct Lut { + t: Vec, } -impl Tree { - fn build(codes: &[(u32, u8)]) -> Tree { - let mut nodes: Vec<[u32; 2]> = vec![[0, 0]]; - for (sym, &(code, len)) in codes.iter().enumerate() { - let mut node = 0usize; - for i in (0..len).rev() { - let bit = ((code >> i) & 1) as usize; - if i == 0 { - nodes[node][bit] = LEAF | sym as u32; - } else { - if nodes[node][bit] == 0 { - nodes.push([0, 0]); - nodes[node][bit] = (nodes.len() - 1) as u32; - } - node = nodes[node][bit] as usize; - } +impl Lut { + fn build(codes: &[(u32, u8)]) -> Lut { + let mut t = vec![0u32; 1 << L1]; + // sub-table width per first-level prefix + let mut width = vec![0u32; 1 << L1]; + for &(code, len) in codes { + let len = len as u32; + if len > L1 { + let p = (code >> (len - L1)) as usize; + width[p] = width[p].max(len - L1); } } - Tree { nodes } + for (p, &w) in width.iter().enumerate() { + if w > 0 { + t[p] = SUB | (t.len() as u32) << 5 | w; + t.resize(t.len() + (1 << w), 0); + } + } + for (sym, &(code, len)) in codes.iter().enumerate() { + let len = len as u32; + let leaf = len << 24 | sym as u32; + if len <= L1 { + let first = (code << (L1 - len)) as usize; + t[first..first + (1 << (L1 - len))].fill(leaf); + } else { + let e = t[(code >> (len - L1)) as usize]; + let (off, w) = (((e & !SUB) >> 5) as usize, e & 31); + let rest = code & ((1 << (len - L1)) - 1); + let first = off + (rest << (w - (len - L1))) as usize; + t[first..first + (1 << (w - (len - L1)))].fill(leaf); + } + } + Lut { t } } + /// The symbol of the codeword at the reader position; `Eof` if the codeword is cut off. #[inline] fn decode(&self, br: &mut BitReader) -> Result { - let mut node = 0usize; - loop { - let bit = br.read_bits(1)? as usize; - let next = self.nodes[node][bit]; - if next & LEAF != 0 { - return Ok(next & !LEAF); - } - if next == 0 { - return Err(Error::Bitstream("invalid Huffman codeword")); - } - node = next as usize; + let bits = peek32(br); + let mut e = self.t[(bits >> (32 - L1)) as usize]; + if e & SUB != 0 { + let w = e & 31; + e = self.t[((e & !SUB) >> 5) as usize + ((bits << L1) >> (32 - w)) as usize]; } + let len = e >> 24; + if len == 0 { + return Err(Error::Bitstream("invalid Huffman codeword")); + } + br.skip(len as usize).map_err(|_| Error::Eof)?; + Ok(e & 0xFF_FFFF) } } -fn trees() -> &'static [Tree; 12] { - static TREES: OnceLock<[Tree; 12]> = OnceLock::new(); - TREES.get_or_init(|| std::array::from_fn(|i| if i == 0 { Tree::build(&SCALEFACTOR) } else { Tree::build(BOOKS[i - 1].codes) })) +/// The next 32 bits (zero padded past the end), from one 8-byte load when possible. +#[inline] +fn peek32(br: &BitReader) -> u32 { + let pos = br.position(); + let byte = pos >> 3; + match br.data().get(byte..byte.saturating_add(8)).and_then(|b| <[u8; 8]>::try_from(b).ok()) { + Some(b) => ((u64::from_be_bytes(b) << (pos & 7)) >> 32) as u32, + None => br.peek(32), + } +} + +fn luts() -> &'static [Lut; 12] { + static LUTS: OnceLock<[Lut; 12]> = OnceLock::new(); + LUTS.get_or_init(|| std::array::from_fn(|i| if i == 0 { Lut::build(&SCALEFACTOR) } else { Lut::build(BOOKS[i - 1].codes) })) } /// Decode one scalefactor difference (`-60..=60`). #[inline] pub fn decode_sf(br: &mut BitReader) -> Result { - Ok(trees()[0].decode(br)? as i32 - 60) + Ok(luts()[0].decode(br)? as i32 - 60) } /// Write one scalefactor difference (`-60..=60`). @@ -125,7 +157,7 @@ pub fn sf_bits(diff: i32) -> u32 { /// Decode one codeword of codebook `cb` (1..=11) into `out[..dim]`, including sign bits and escapes. pub fn decode_spectral(br: &mut BitReader, cb: u8, out: &mut [i32]) -> Result<()> { let b = book(cb); - let mut idx = trees()[cb as usize].decode(br)?; + let mut idx = luts()[cb as usize].decode(br)?; for i in (0..b.dim).rev() { out[i] = (idx % b.modulo) as i32; idx /= b.modulo; @@ -198,6 +230,27 @@ pub fn count_bits(cb: u8, vals: &[i32]) -> u32 { bits } +/// Bits for codebooks `cb` and `cb + 1` in one pass. `cb` must be odd and below 11: each such pair +/// shares dimension, modulus and signedness, so the index is computed once and only the length table +/// differs. +pub fn count_bits_pair(cb: u8, vals: &[i32]) -> (u32, u32) { + debug_assert!(cb % 2 == 1 && cb < ESC_HCB); + let (b, b2) = (book(cb), book(cb + 1)); + let (mut x, mut y) = (0u32, 0u32); + for chunk in vals.chunks_exact(b.dim) { + let i = index_of(b, chunk); + x += b.codes[i].1 as u32; + y += b2.codes[i].1 as u32; + if !b.signed { + // sign bits (no escapes below codebook 11) + let nz = chunk.iter().filter(|&&v| v != 0).count() as u32; + x += nz; + y += nz; + } + } + (x, y) +} + /// Write `vals` with codebook `cb` (1..=11). pub fn write_spectral(bw: &mut BitWriter, cb: u8, vals: &[i32]) { let b = book(cb); @@ -259,6 +312,53 @@ mod tests { } } + /// Every codeword, followed by arbitrary bits, decodes to its own symbol and length. + #[test] + fn lut_decodes_every_codeword() { + let mut seed = 7u32; + for (i, lut) in luts().iter().enumerate() { + let codes: &[(u32, u8)] = if i == 0 { &SCALEFACTOR } else { BOOKS[i - 1].codes }; + for (sym, &(code, len)) in codes.iter().enumerate() { + for _ in 0..3 { + seed = seed.wrapping_mul(1664525).wrapping_add(1013904223); + let mut bw = BitWriter::new(); + bw.write_bits(code, len as u32); + bw.write_bits(seed >> 8, 24); + let data = bw.finish(); + let mut br = BitReader::new(&data); + assert_eq!(lut.decode(&mut br).unwrap(), sym as u32, "book {i}"); + assert_eq!(br.position(), len as usize); + } + } + } + } + + /// A codeword cut off by the end of the data is `Eof`, also at an unaligned position. + #[test] + fn truncated_codeword_is_eof() { + for (i, lut) in luts().iter().enumerate() { + let codes: &[(u32, u8)] = if i == 0 { &SCALEFACTOR } else { BOOKS[i - 1].codes }; + let &(code, len) = codes.iter().max_by_key(|c| c.1).unwrap(); + for lead in 0..8u32 { + let mut bw = BitWriter::new(); + bw.write_bits(0, lead); + bw.write_bits(code, len as u32); + let total = lead as usize + len as usize; + let data = bw.finish(); + // drop the last byte if that cuts the codeword, else the codeword is whole + let short = &data[..data.len() - 1]; + let mut br = BitReader::new(short); + br.skip(lead as usize).unwrap(); + let r = lut.decode(&mut br); + if short.len() * 8 < total { + assert!(matches!(r, Err(Error::Eof)), "book {i} lead {lead}"); + } else { + assert!(r.is_ok()); + } + } + } + } + #[test] fn spectral_roundtrip_every_book() { for cb in 1..=11u8 { @@ -291,6 +391,25 @@ mod tests { } } + #[test] + fn pair_count_matches_single_count() { + let mut seed = 99u32; + for cb in (1..11u8).step_by(2) { + let b = book(cb); + for round in 0..40 { + let mut vals = Vec::new(); + for _ in 0..(8 * b.dim) { + seed = seed.wrapping_mul(1664525).wrapping_add(1013904223); + let m = ((seed >> 8) % (b.lav + 1)) as i32; + let m = if round % 4 == 0 { 0 } else { m }; + let neg = if b.signed { seed & 0x8000_0000 != 0 } else { seed & 1 == 1 }; + vals.push(if neg { -m } else { m }); + } + assert_eq!(count_bits_pair(cb, &vals), (count_bits(cb, &vals), count_bits(cb + 1, &vals)), "codebook {cb}"); + } + } + } + #[test] fn scalefactor_roundtrip() { let mut bw = BitWriter::new(); diff --git a/crates/audio/src/aac/huffman_tables.rs b/crates/audio/src/aac/huffman_tables.rs index c65c7df6..e48f5fad 100644 --- a/crates/audio/src/aac/huffman_tables.rs +++ b/crates/audio/src/aac/huffman_tables.rs @@ -5,8 +5,6 @@ //! (every table is a complete prefix code — Kraft sum exactly 1 — covering every index once; see //! `tests::codebooks_are_complete_prefix_codes`). -#![allow(clippy::unreadable_literal)] - pub(crate) static SCALEFACTOR: [(u32, u8); 121] = [ (0x3ffe8, 18), (0x3ffe6, 18), diff --git a/crates/audio/src/aac/mdct.rs b/crates/audio/src/aac/mdct.rs index b829030d..4549cb7d 100644 --- a/crates/audio/src/aac/mdct.rs +++ b/crates/audio/src/aac/mdct.rs @@ -7,7 +7,11 @@ use std::sync::OnceLock; /// Radix-2 complex FFT (forward, `e^{-2πi nk/P}`). struct Fft { n: usize, - twiddle: Vec<[f32; 2]>, + /// Twiddle of the length-2 stage, `e^{-2πi 0/n}`. + w0: [f32; 2], + /// Twiddles of the stages of length 4, 8, ..., n, one contiguous run of `len / 2` per stage. + stage_re: Vec, + stage_im: Vec, rev: Vec, } @@ -16,13 +20,23 @@ impl Fft { assert!(n.is_power_of_two()); let bits = n.trailing_zeros(); let rev = (0..n as u32).map(|i| if bits == 0 { 0 } else { i.reverse_bits() >> (32 - bits) }).collect(); - let twiddle = (0..n / 2) + let twiddle: Vec<[f32; 2]> = (0..n / 2) .map(|k| { let a = -2.0 * PI * k as f64 / n as f64; [a.cos() as f32, a.sin() as f32] }) .collect(); - Fft { n, twiddle, rev } + let (mut stage_re, mut stage_im) = (Vec::new(), Vec::new()); + let mut len = 4; + while len <= n { + let step = n / len; + for k in 0..len / 2 { + stage_re.push(twiddle[k * step][0]); + stage_im.push(twiddle[k * step][1]); + } + len <<= 1; + } + Fft { n, w0: twiddle.first().copied().unwrap_or([1.0, 0.0]), stage_re, stage_im, rev } } fn run(&self, re: &mut [f32], im: &mut [f32]) { @@ -34,24 +48,80 @@ impl Fft { im.swap(i, j); } } - let mut len = 2; - while len <= n { + if n >= 2 { + let [wr, wi] = self.w0; + let ((re2, _), (im2, _)) = (re.as_chunks_mut::<2>(), im.as_chunks_mut::<2>()); + for (r, i) in re2.iter_mut().zip(im2) { + let xr = r[1] * wr - i[1] * wi; + let xi = r[1] * wi + i[1] * wr; + r[1] = r[0] - xr; + i[1] = i[0] - xi; + r[0] += xr; + i[0] += xi; + } + } + let mut len = 4; + let mut off = 0; + // Two radix-2 stages (len, 2 * len) per pass over the data: the same butterflies in the same + // order per element, so the result is bit-identical to running the stages one by one. + while len * 2 <= n { + let h = len / 2; + let (t1r, t1i) = (&self.stage_re[off..off + h], &self.stage_im[off..off + h]); + let (t2r, t2i) = (&self.stage_re[off + h..off + 3 * h], &self.stage_im[off + h..off + 3 * h]); + off += 3 * h; + for (rc, ic) in re.chunks_exact_mut(4 * h).zip(im.chunks_exact_mut(4 * h)) { + let (ra, rest) = rc.split_at_mut(h); + let (rb, rest) = rest.split_at_mut(h); + let (rc2, rd) = rest.split_at_mut(h); + let (ia, rest) = ic.split_at_mut(h); + let (ib, rest) = rest.split_at_mut(h); + let (ic2, id) = rest.split_at_mut(h); + for k in 0..h { + let (w1r, w1i) = (t1r[k], t1i[k]); + // stage `len`: (a, b) and (c, d) + let xr = rb[k] * w1r - ib[k] * w1i; + let xi = rb[k] * w1i + ib[k] * w1r; + let (br, bi) = (ra[k] - xr, ia[k] - xi); + let (ar, ai) = (ra[k] + xr, ia[k] + xi); + let xr = rd[k] * w1r - id[k] * w1i; + let xi = rd[k] * w1i + id[k] * w1r; + let (dr, di) = (rc2[k] - xr, ic2[k] - xi); + let (cr, ci) = (rc2[k] + xr, ic2[k] + xi); + // stage `2 * len`: (a, c) with w2[k], (b, d) with w2[k + h] + let (w2r, w2i) = (t2r[k], t2i[k]); + let xr = cr * w2r - ci * w2i; + let xi = cr * w2i + ci * w2r; + rc2[k] = ar - xr; + ic2[k] = ai - xi; + ra[k] = ar + xr; + ia[k] = ai + xi; + let (w2r, w2i) = (t2r[k + h], t2i[k + h]); + let xr = dr * w2r - di * w2i; + let xi = dr * w2i + di * w2r; + rd[k] = br - xr; + id[k] = bi - xi; + rb[k] = br + xr; + ib[k] = bi + xi; + } + } + len <<= 2; + } + if len <= n { let half = len / 2; - let step = n / len; - for start in (0..n).step_by(len) { + let (twr, twi) = (&self.stage_re[off..off + half], &self.stage_im[off..off + half]); + for (rc, ic) in re.chunks_exact_mut(len).zip(im.chunks_exact_mut(len)) { + let (ra, rb) = rc.split_at_mut(half); + let (ia, ib) = ic.split_at_mut(half); for k in 0..half { - let [wr, wi] = self.twiddle[k * step]; - let a = start + k; - let b = a + half; - let xr = re[b] * wr - im[b] * wi; - let xi = re[b] * wi + im[b] * wr; - re[b] = re[a] - xr; - im[b] = im[a] - xi; - re[a] += xr; - im[a] += xi; + let (wr, wi) = (twr[k], twi[k]); + let xr = rb[k] * wr - ib[k] * wi; + let xi = rb[k] * wi + ib[k] * wr; + rb[k] = ra[k] - xr; + ib[k] = ia[k] - xi; + ra[k] += xr; + ia[k] += xi; } } - len <<= 1; } } } @@ -83,19 +153,22 @@ impl Mdct { Mdct { m, fft: Fft::new(q), pre, post } } - /// DCT-IV: `out[k] = Σ u[n] cos(π/M (n+½)(k+½))`. - fn dct4(&self, u: &[f32], out: &mut [f32]) { + /// Working buffers for this transform size. + pub fn scratch(&self) -> MdctScratch { + MdctScratch { re: vec![0.0; self.m / 2], im: vec![0.0; self.m / 2], u: vec![0.0; self.m] } + } + + /// DCT-IV: `out[k] = Σ u[n] cos(π/M (n+½)(k+½))`; `re` and `im` (length M/2) are overwritten. + fn dct4(&self, u: &[f32], out: &mut [f32], re: &mut [f32], im: &mut [f32]) { let m = self.m; let q = m / 2; - let mut re = vec![0f32; q]; - let mut im = vec![0f32; q]; for n in 0..q { let (a, b) = (u[2 * n], u[m - 1 - 2 * n]); let [c, s] = self.pre[n]; re[n] = a * c - b * s; im[n] = a * s + b * c; } - self.fft.run(&mut re, &mut im); + self.fft.run(re, im); for k in 0..q { let [c, s] = self.post[k]; let wr = re[k] * c - im[k] * s; @@ -109,21 +182,28 @@ impl Mdct { pub fn forward(&self, x: &[f32], out: &mut [f32]) { let m = self.m; let h = m / 2; - let mut u = vec![0f32; m]; + let MdctScratch { re, im, mut u } = self.scratch(); for n in 0..h { // quarters a, b, c, d of length M/2 u[n] = -x[m + h - 1 - n] - x[m + h + n]; u[h + n] = x[n] - x[m - 1 - n]; } - self.dct4(&u, out); + let (mut re, mut im) = (re, im); + self.dct4(&u, out, &mut re, &mut im); } /// Inverse MDCT (unscaled): `y[n] = Σ X[k] cos(π/M (n + ½ + M/2)(k + ½))`, `y` has length 2M. + #[cfg(test)] pub fn inverse(&self, spec: &[f32], y: &mut [f32]) { + self.inverse_with(spec, y, &mut self.scratch()); + } + + /// [`Self::inverse`] with caller-provided working buffers from [`Self::scratch`]. + pub fn inverse_with(&self, spec: &[f32], y: &mut [f32], sc: &mut MdctScratch) { let m = self.m; let h = m / 2; - let mut u = vec![0f32; m]; - self.dct4(spec, &mut u); + let MdctScratch { re, im, u } = sc; + self.dct4(spec, u, re, im); for n in 0..h { y[n] = u[h + n]; y[m - 1 - n] = -u[h + n]; @@ -133,6 +213,14 @@ impl Mdct { } } +/// Working buffers of one [`Mdct`] size, reusable across transforms. +#[derive(Clone)] +pub struct MdctScratch { + re: Vec, + im: Vec, + u: Vec, +} + /// Cached transforms for 1024 and 128 lines. pub fn mdct_long() -> &'static Mdct { static M: OnceLock = OnceLock::new(); diff --git a/crates/audio/src/aac/tables.rs b/crates/audio/src/aac/tables.rs index 6dc6fb6f..87fb55d0 100644 --- a/crates/audio/src/aac/tables.rs +++ b/crates/audio/src/aac/tables.rs @@ -101,12 +101,22 @@ pub fn tns_max_order(short: bool) -> usize { if short { 7 } else { 12 } } -/// `x^(4/3)` for the integer quantised magnitudes 0..=8191. -pub fn pow43(q: u32) -> f32 { +/// The `x^(4/3)` table for the integer quantised magnitudes 0..=8191. +pub fn pow43_table() -> &'static [f32; 8192] { use std::sync::OnceLock; - static TABLE: OnceLock> = OnceLock::new(); - let t = TABLE.get_or_init(|| (0..8192u32).map(|i| (i as f64).powf(4.0 / 3.0) as f32).collect()); - match t.get(q as usize) { + static TABLE: OnceLock> = OnceLock::new(); + TABLE.get_or_init(|| { + let mut t = Box::new([0f32; 8192]); + for (i, v) in t.iter_mut().enumerate() { + *v = (i as f64).powf(4.0 / 3.0) as f32; + } + t + }) +} + +/// `x^(4/3)` for an integer quantised magnitude (table lookup up to 8191). +pub fn pow43(q: u32) -> f32 { + match pow43_table().get(q as usize) { Some(v) => *v, None => (q as f64).powf(4.0 / 3.0) as f32, } @@ -115,6 +125,16 @@ pub fn pow43(q: u32) -> f32 { /// `2^(0.25 * (sf - 100))`, the dequantisation gain of scalefactor `sf`. #[inline] pub fn sf_gain(sf: i32) -> f32 { + use std::sync::OnceLock; + static TABLE: OnceLock<[f32; 256]> = OnceLock::new(); + let t = TABLE.get_or_init(|| std::array::from_fn(|i| sf_gain_exact(i as i32))); + match usize::try_from(sf).ok().and_then(|i| t.get(i)) { + Some(g) => *g, + None => sf_gain_exact(sf), + } +} + +fn sf_gain_exact(sf: i32) -> f32 { (0.25 * (sf as f64 - 100.0)).exp2() as f32 } diff --git a/crates/audio/src/ac3/golden_frames.bin b/crates/audio/src/ac3/golden_frames.bin new file mode 100644 index 00000000..4f8a589a Binary files /dev/null and b/crates/audio/src/ac3/golden_frames.bin differ diff --git a/crates/audio/src/ac3/mod.rs b/crates/audio/src/ac3/mod.rs index e04142b8..a8e8920d 100644 --- a/crates/audio/src/ac3/mod.rs +++ b/crates/audio/src/ac3/mod.rs @@ -38,6 +38,7 @@ pub use program::{ProgrammeDecoder, ProgrammeFrame}; pub mod encoder; mod tables; +use craft_media::bitstream::BitReader; use std::sync::OnceLock; /// Errors from header parsing and decoding. @@ -165,33 +166,38 @@ pub fn parse_header(b: &[u8]) -> Result
{ }) } -/// MSB-first bit reader; reads past the end give zeros. +/// MSB-first bit reader over the shared [`BitReader`]. Reads past the end give zeros, and +/// `skip` may move past the end: a damaged frame decodes to something, and [`Bits::overrun`] +/// tells afterwards whether it ran out of data. struct Bits<'a> { - d: &'a [u8], - pos: usize, + r: BitReader<'a>, } impl<'a> Bits<'a> { fn new(d: &'a [u8]) -> Self { - Self { d, pos: 0 } + Self { r: BitReader::new(d) } } + #[inline] fn read(&mut self, n: u32) -> u32 { - let mut v = 0u32; - for _ in 0..n { - let byte = self.d.get(self.pos >> 3).copied().unwrap_or(0); - v = (v << 1) | ((byte >> (7 - (self.pos & 7))) & 1) as u32; - self.pos += 1; - } + let v = self.r.peek(n); + self.skip(n as usize); v } + #[inline] fn bit(&mut self) -> bool { self.read(1) == 1 } + #[inline] fn skip(&mut self, n: usize) { - self.pos += n; + self.r.seek_bits(self.r.position() + n); } fn overrun(&self) -> bool { - self.pos > self.d.len() * 8 + self.r.position() > self.r.data().len() * 8 + } + /// The position in bits. + #[cfg(test)] + fn pos(&self) -> usize { + self.r.position() } } @@ -1276,7 +1282,7 @@ impl Decoder { } // mantissas #[cfg(test)] - let mantissa_start = r.pos; + let mantissa_start = r.pos(); let mut gr = Groups::default(); let mut cpl = [0f32; 256]; let mut got_cpl = false; @@ -1303,7 +1309,7 @@ impl Decoder { } } #[cfg(test)] - self.block_bit_positions.push((mantissa_start, r.pos)); + self.block_bit_positions.push((mantissa_start, r.pos())); // decoupling if f.cplinu { for ch in 0..nf { diff --git a/crates/audio/src/ac3/tests.rs b/crates/audio/src/ac3/tests.rs index 830347a4..a9a83c10 100644 --- a/crates/audio/src/ac3/tests.rs +++ b/crates/audio/src/ac3/tests.rs @@ -272,3 +272,244 @@ fn eac3_garbage_never_panics() { } } } + +/// Golden hashes of the decoder output, so that a change of the bit reader cannot change a +/// decoded sample or an error. The values were made with the private bit-by-bit reader, before it +/// was replaced by the shared one: every frame below is decoded with a +/// fresh decoder, and every output sample bit (as `f32::to_bits`, channel after channel) and every +/// error (its `Debug` text) goes into an FNV-1a 64 hash. The frames are deterministic: valid +/// AC-3 and E-AC-3 syncframe headers (all fscod, frmsizecod, bsid, acmod and lfeon values in +/// the tables below) with random, zero-heavy or all-ones payload bits, plus truncated and +/// bit-flipped variants of the same frames. +struct Fnv(u64); + +impl Fnv { + fn new() -> Self { + Fnv(0xcbf2_9ce4_8422_2325) + } + fn bytes(&mut self, b: &[u8]) { + for &x in b { + self.0 = (self.0 ^ u64::from(x)).wrapping_mul(0x0000_0100_0000_01b3); + } + } +} + +struct Rng(u32); + +impl Rng { + fn next(&mut self) -> u32 { + self.0 ^= self.0 << 13; + self.0 ^= self.0 >> 17; + self.0 ^= self.0 << 5; + self.0 + } +} + +/// One payload byte for the given flavour: random, zero-heavy, or ones-heavy. +fn payload_byte(rng: &mut Rng, flavour: usize) -> u8 { + let v = rng.next(); + match flavour % 3 { + 0 => (v >> 7) as u8, + 1 => ((v >> 7) as u8) & ((v >> 19) as u8) & ((v >> 3) as u8), + _ => ((v >> 7) as u8) | ((v >> 19) as u8) | ((v >> 3) as u8), + } +} + +/// See `synthetic_frames`: 12 syncframes, back to back. +const REAL_FRAMES: &[u8] = include_bytes!("golden_frames.bin"); + +fn synthetic_frames() -> Vec> { + let mut rng = Rng(0x00C0_FFEE); + let mut frames = Vec::new(); + // AC-3: bsid 8 and the half / quarter rate variants 9 and 10, every fscod and a spread of + // frmsizecod, acmod 0 to 7, lfeon off and on. + for n in 0..1200usize { + let fscod = n % 3; + let frmsizecod = (n * 7 + n / 3) % 38; + let bsid = [8u8, 8, 8, 9, 10, 6][n % 6]; + let acmod = (n / 3) % 8; + let lfeon = (n / 24) % 2; + let Ok(h) = parse_header(&[0x0B, 0x77, 0, 0, ((fscod << 6) | frmsizecod) as u8, bsid << 3, (acmod as u8) << 5, 0]) else { continue }; + let mut f = vec![0u8; h.frame_bytes]; + for b in f.iter_mut() { + *b = payload_byte(&mut rng, n); + } + f[0] = 0x0B; + f[1] = 0x77; + f[4] = ((fscod << 6) | frmsizecod) as u8; + f[5] = bsid << 3; + // acmod in the top 3 bits; the mix-level bits that follow depend on acmod and are random + f[6] = (f[6] & 0x1F) | ((acmod as u8) << 5); + // lfeon follows the optional 2-bit fields (cmixlev, surmixlev, dsurmod) + let mut at = 3; + if acmod & 1 != 0 && acmod != 1 { + at += 2; + } + if acmod & 4 != 0 { + at += 2; + } + if acmod == 2 { + at += 2; + } + let bit = 48 + at; + f[bit / 8] = (f[bit / 8] & !(0x80 >> (bit % 8))) | ((lfeon as u8) * (0x80 >> (bit % 8))); + frames.push(f); + } + // E-AC-3: strmtyp 0 / 1 / 2, fscod 0..=3 (3 = reduced rate), numblkscod, acmod, lfeon + for n in 0..800usize { + let words = 24 + (n * 13) % 400; + let fscod = n % 4; + let blk = (n / 4) % 4; + let acmod = (n / 16) % 8; + let lfeon = (n / 128) % 2; + let strmtyp = [0u32, 0, 0, 1, 2][n % 5]; + let bsid = 11 + (n % 6) as u32; + let mut f = vec![0u8; words * 2]; + for b in f.iter_mut() { + *b = payload_byte(&mut rng, n); + } + let head: u32 = (strmtyp << 30) + | (((n % 7) as u32) << 27) + | (((words - 1) as u32) << 16) + | ((fscod as u32) << 14) + | ((if fscod == 3 { (n / 7 % 3) as u32 } else { blk as u32 }) << 12) + | ((acmod as u32) << 9) + | ((lfeon as u32) << 8) + | (bsid << 3) + | (f[5] as u32 & 7); + f[0] = 0x0B; + f[1] = 0x77; + f[2..6].copy_from_slice(&head.to_be_bytes()); + frames.push(f); + } + // real frames (two of each stream, made with ffmpeg's ac3 / eac3 encoders from tones and noise: + // mono 64k 32 kHz, coupled stereo 96k 44.1 kHz, noise stereo 96k, 5.1 192k; E-AC-3 noise + // stereo 96k and 5.1 192k) reach exponents, bit allocation and mantissas of every block, which + // random payload bits almost never do. Each one is decoded as is, with 60 bit-flipped copies + // and cut at every 16th byte. + let mut p = 0; + while let Ok(h) = parse_header(&REAL_FRAMES[p..]) { + let f = &REAL_FRAMES[p..p + h.frame_bytes]; + p += h.frame_bytes; + frames.push(f.to_vec()); + for _ in 0..60 { + let mut m = f.to_vec(); + for _ in 0..1 + rng.next() % 3 { + let bit = 64 + rng.next() as usize % (m.len() * 8 - 64); + m[bit / 8] ^= 0x80 >> (bit % 8); + } + frames.push(m); + } + for cut in (8..f.len()).step_by(16) { + frames.push(f[..cut].to_vec()); + } + } + // truncated and bit-flipped variants of every 5th frame + let base = frames.clone(); + for (i, f) in base.iter().enumerate().step_by(5) { + let mut t = f.clone(); + t.truncate(rng.next() as usize % f.len()); + frames.push(t); + let mut m = f.clone(); + for _ in 0..1 + rng.next() % 6 { + let bit = rng.next() as usize % (m.len() * 8); + m[bit / 8] ^= 0x80 >> (bit % 8); + } + frames.push(m); + // cut the payload to just a few bytes after the header, so that the reads run past the end + let mut s = f.clone(); + s.truncate((8 + i % 40).min(f.len())); + frames.push(s); + } + frames +} + +#[test] +fn decoded_output_matches_golden_hashes() { + let (mut all, mut ok, mut err, mut audio) = (Fnv::new(), 0usize, 0usize, 0usize); + for (i, f) in synthetic_frames().iter().enumerate() { + let mut dec = Decoder::new(); + // decode each frame twice with the same decoder: the second pass also covers the + // E-AC-3 carry-over and the overlap state + for _ in 0..2 { + match dec.decode(f) { + Ok(d) => { + ok += 1; + audio += usize::from(!d.channels.is_empty()); + all.bytes(&[1]); + for ch in &d.channels { + for s in ch { + all.bytes(&s.to_bits().to_le_bytes()); + } + } + } + Err(e) => { + err += 1; + all.bytes(format!("{i}:{e:?}").as_bytes()); + } + } + } + } + eprintln!("golden: ok {ok} ({audio} with audio) err {err} hash {:#018x}", all.0); + assert_eq!((ok, err, all.0), GOLDEN); +} + +const GOLDEN: (usize, usize, u64) = (2960, 7024, 0x1dc9_31d6_0e89_e8c8); + +/// The private reader this crate used before it moved to the shared `BitReader`, verbatim. +struct OldBits<'a> { + d: &'a [u8], + pos: usize, +} + +impl<'a> OldBits<'a> { + fn new(d: &'a [u8]) -> Self { + Self { d, pos: 0 } + } + fn read(&mut self, n: u32) -> u32 { + let mut v = 0u32; + for _ in 0..n { + let byte = self.d.get(self.pos >> 3).copied().unwrap_or(0); + v = (v << 1) | ((byte >> (7 - (self.pos & 7))) & 1) as u32; + self.pos += 1; + } + v + } + fn bit(&mut self) -> bool { + self.read(1) == 1 + } + fn skip(&mut self, n: usize) { + self.pos += n; + } + fn overrun(&self) -> bool { + self.pos > self.d.len() * 8 + } +} + +#[test] +fn bits_match_the_old_bit_by_bit_reader() { + let mut rng = Rng(0x5EED_0001); + for round in 0..3000usize { + let len = (rng.next() % 40) as usize; + let data: Vec = (0..len).map(|_| payload_byte(&mut rng, round)).collect(); + let (mut old, mut new) = (OldBits::new(&data), Bits::new(&data)); + for op in 0..80 { + match rng.next() % 5 { + // widths 0 to 40 include more than 32 bits, where only the low 32 bits are kept + 0 | 1 => { + let n = rng.next() % 41; + assert_eq!(old.read(n), new.read(n), "round {round} op {op}: read({n})"); + } + 2 => assert_eq!(old.bit(), new.bit(), "round {round} op {op}: bit"), + 3 => { + let n = (rng.next() % 70) as usize; + old.skip(n); + new.skip(n); + } + _ => {} + } + assert_eq!(old.overrun(), new.overrun(), "round {round} op {op}: overrun"); + assert_eq!(old.pos, new.pos(), "round {round} op {op}: position"); + } + } +} diff --git a/crates/audio/src/mpa/mod.rs b/crates/audio/src/mpa/mod.rs index a8cd8e9b..7ba224d7 100644 --- a/crates/audio/src/mpa/mod.rs +++ b/crates/audio/src/mpa/mod.rs @@ -190,6 +190,9 @@ type Subbands = [[[f32; tables::SUBBANDS]; 36]; 2]; pub struct Decoder { /// The synthesis filter's `V` vector per channel (§2.4.3.2, Figure A.2). v: [[f32; 1024]; 2], + /// Where the newest 64 values of each `V` start: `V` is a ring, logical index `l` sits at + /// `(pos + l) % 1024`. + pos: [usize; 2], sb: Box, } @@ -202,12 +205,13 @@ impl Default for Decoder { impl Decoder { /// A decoder with empty synthesis filter state. pub fn new() -> Self { - Self { v: [[0.0; 1024]; 2], sb: Box::new([[[0.0; 32]; 36]; 2]) } + Self { v: [[0.0; 1024]; 2], pos: [0; 2], sb: Box::new([[[0.0; 32]; 36]; 2]) } } /// Forget the filter history. pub fn reset(&mut self) { self.v = [[0.0; 1024]; 2]; + self.pos = [0; 2]; } /// Decode the frame at the start of `frame` (which may run on into the next frames; only @@ -237,12 +241,12 @@ impl Decoder { 36 }; let window = tables::window(); - let matrix = tables::matrix(); + let matrix = tables::matrix_t(); let mut channels = vec![vec![0f32; slots * 32]; nch]; for (ch, out) in channels.iter_mut().enumerate() { - let (Some(v), Some(sb)) = (self.v.get_mut(ch), self.sb.get(ch)) else { continue }; + let (Some(v), Some(pos), Some(sb)) = (self.v.get_mut(ch), self.pos.get_mut(ch), self.sb.get(ch)) else { continue }; for (slot, o) in sb.iter().take(slots).zip(out.as_chunks_mut::<32>().0) { - synthesize(v, matrix, window, slot, o); + synthesize(v, pos, matrix, window, slot, o); } } Ok(Frame { header, channels }) @@ -414,18 +418,43 @@ fn layer2(r: &mut BitReader, h: &Header, out: &mut Subbands) -> Result<()> { /// The synthesis subband filter (§2.4.3.2, Figure A.2) for one time slot: 32 subband samples in, /// 32 PCM samples out. -fn synthesize(v: &mut [f32; 1024], n: &[[f32; 32]; 64], d: &[f32; 512], s: &[f32; 32], out: &mut [f32; 32]) { - v.copy_within(0..960, 64); - for (vi, row) in v.iter_mut().zip(n.iter()) { - *vi = row.iter().zip(s.iter()).map(|(a, b)| a * b).sum(); - } - for (j, o) in out.iter_mut().enumerate() { - let mut acc = 0f32; - for i in 0..8 { - acc += v[i * 128 + j] * d[i * 64 + j] + v[i * 128 + 96 + j] * d[i * 64 + 32 + j]; +fn synthesize(v: &mut [f32; 1024], pos: &mut usize, n: &[[f32; 36]; 32], d: &[f32; 512], s: &[f32; 32], out: &mut [f32; 32]) { + // shift V by 64: move the ring's start back instead of copying 960 values + *pos = (*pos + 1024 - 64) % 1024; + let p = *pos; + // Only 33 of the 64 rows of N are distinct (N[32 - i] = -N[i], N[96 - i] = N[i]): compute the + // sums of rows 0..=16 and 48..64 (padded to 36 for the vector width), accumulated column by + // column, each sum adding its terms in k order from -0.0 like `Iterator::sum` + let mut acc = [-0.0f32; 36]; + for (col, &sk) in n.iter().zip(s.iter()) { + for (a, &c) in acc.iter_mut().zip(col.iter()) { + *a += c * sk; } - *o = acc; } + if let Some(head) = v.get_mut(p..).and_then(|t| t.first_chunk_mut::<64>()) { + for i in 0..16 { + head[i] = acc[i]; + head[32 - i] = -acc[i]; + } + head[16] = acc[16]; + for j in 0..16 { + head[48 + j] = acc[17 + j]; + if j > 0 { + head[48 - j] = acc[17 + j]; + } + } + } + let mut o = [0f32; 32]; + // p is a multiple of 64, so each 32-value run below is contiguous in the ring + let run = |l: usize| v.get((p + l) % 1024..).and_then(|t| t.first_chunk::<32>()); + for (i, db) in d.as_chunks::<64>().0.iter().enumerate() { + let (Some(va), Some(vc)) = (run(i * 128), run(i * 128 + 96)) else { continue }; + let (Some(da), Some(dc)) = (db.first_chunk::<32>(), db.last_chunk::<32>()) else { continue }; + for ((((o, &a), &b), &c), &e) in o.iter_mut().zip(va).zip(da).zip(vc).zip(dc) { + *o += a * b + c * e; + } + } + *out = o; } #[cfg(test)] diff --git a/crates/audio/src/mpa/tables.rs b/crates/audio/src/mpa/tables.rs index 11d56c5a..11a60c7d 100644 --- a/crates/audio/src/mpa/tables.rs +++ b/crates/audio/src/mpa/tables.rs @@ -72,6 +72,23 @@ pub(crate) fn matrix() -> &'static [[f32; 32]; 64] { }) } +/// The 33 distinct rows of [`matrix`] (0..=16 and 48..64), transposed to `[k][row]` and padded +/// with zeros to 36. +pub(crate) fn matrix_t() -> &'static [[f32; 36]; 32] { + static NT: OnceLock<[[f32; 36]; 32]> = OnceLock::new(); + NT.get_or_init(|| { + let n = matrix(); + let mut t = [[0f32; 36]; 32]; + for (k, col) in t.iter_mut().enumerate() { + for (r, c) in col.iter_mut().take(33).enumerate() { + let i = if r <= 16 { r } else { 48 + r - 17 }; + *c = n.get(i).and_then(|row| row.get(k)).copied().unwrap_or(0.0); + } + } + t + }) +} + /// Scalefactors (Table 3-B.1): `2^(1 - index / 3)`. Index 63 is not used by encoders; it continues /// the progression rather than failing the frame. pub(crate) fn scalefactor(index: u32) -> f32 { diff --git a/crates/audio/src/mpa/tests.rs b/crates/audio/src/mpa/tests.rs index 490d7259..941509bc 100644 --- a/crates/audio/src/mpa/tests.rs +++ b/crates/audio/src/mpa/tests.rs @@ -152,3 +152,99 @@ fn random_frames_never_panic() { } } } + +/// Output of the decoder, bit for bit, on deterministic synthetic streams. +/// +/// Each stream is a run of valid headers (mono, stereo, joint stereo at several bounds and dual +/// channel; Layer I and II; MPEG-1 and MPEG-2; several bit rates and sampling frequencies) +/// followed by xorshift noise, decoded by one `Decoder` so that the synthesis filter's history +/// carries from frame to frame. Frames the decoder rejects are counted and skipped. The hash is +/// FNV-1a over the bits of every decoded sample. +/// +/// The golden values were captured from the decoder before the synthesis filterbank was +/// rewritten (`synthesize` with the shifting `V` buffer and the 64-row matrix): check out that +/// commit, run `cargo test -p filmcraft-mpa decoded_samples_are_bit_exact -- --nocapture`, and +/// copy the printed `(ok, hash)` pairs. A change of the filter's arithmetic (summation order, +/// grouping, fused multiply-add) shows up here; a legitimate one needs new goldens and a reason. +#[test] +fn decoded_samples_are_bit_exact() { + // (lsf, layer, bitrate_index, rate_index, mode, mode_extension) + const STREAMS: &[(bool, u8, u8, u8, u8, u8)] = &[ + (false, 2, 10, 1, 0, 0), + (false, 2, 14, 1, 3, 0), + (false, 2, 6, 0, 1, 0), + (false, 2, 9, 0, 1, 3), + (false, 2, 5, 2, 2, 0), + (false, 2, 2, 1, 0, 0), + (true, 2, 8, 0, 1, 2), + (true, 2, 4, 2, 3, 0), + (false, 1, 12, 0, 1, 2), + (false, 1, 14, 1, 0, 0), + (false, 1, 8, 2, 3, 0), + (true, 1, 10, 0, 2, 0), + ]; + let mut got = vec![]; + for (n, &(lsf, layer, br, rate, mode, ext)) in STREAMS.iter().enumerate() { + let mut x = 0x9E37_79B9u32 ^ (n as u32).wrapping_mul(0x85EB_CA6B); + let mut rnd = move || { + x ^= x << 13; + x ^= x >> 17; + x ^= x << 5; + x + }; + let mut d = Decoder::new(); + let (mut ok, mut hash) = (0u32, 0xCBF2_9CE4_8422_2325u64); + for _ in 0..60 { + let h = header(lsf, layer, br, rate, mode, ext); + let len = parse_header(&h).map_or(0, |p| p.frame_bytes); + let mut f = h.to_vec(); + // sparse noise: most allocations are zero, so many frames are valid and quiet + f.extend((0..len.saturating_sub(4)).map(|_| (rnd() & rnd() & rnd()) as u8)); + if let Ok(fr) = d.decode(&f) { + ok += 1; + for s in fr.channels.iter().flatten() { + hash = (hash ^ u64::from(s.to_bits())).wrapping_mul(0x0000_0100_0000_01B3); + } + } + } + got.push((ok, hash)); + } + println!("{got:?}"); + const GOLDEN: &[(u32, u64)] = &[ + (59, 0xea8448e85fe3df1b), + (60, 0xc6ecc7f674c186d4), + (60, 0xb38a4d7cba7bfbde), + (59, 0xf37f8a5abff50363), + (60, 0x1e489c563c40c852), + (49, 0x52ab5c8f9db13cc1), + (60, 0x401b8039ec831783), + (60, 0x3555c5d49d12ecca), + (59, 0xb063c5f8c79edaef), + (58, 0xcedb49aa85eb2e29), + (58, 0x3ee32a2c6b9da56b), + (57, 0xbb913530021557b1), + ]; + assert_eq!(got, GOLDEN); +} + +/// `reset` forgets the filter history, including where the ring buffer of `V` starts. +#[test] +fn reset_matches_a_fresh_decoder() { + let mut x = 0x2545_F491u32; + let mut f = header(false, 2, 10, 1, 0, 0).to_vec(); + f.extend((0..572).map(|_| { + x ^= x << 13; + x ^= x >> 17; + x ^= x << 5; + (x & (x >> 7) & (x >> 13)) as u8 + })); + let first = Decoder::new().decode(&f).unwrap(); + let mut d = Decoder::new(); + for _ in 0..3 { + d.decode(&f).unwrap(); + } + d.reset(); + let again = d.decode(&f).unwrap(); + assert!(first.channels.iter().flatten().any(|&s| s != 0.0)); + assert_eq!(first.channels, again.channels); +} diff --git a/crates/audio/src/opus/decode/celt/bands.rs b/crates/audio/src/opus/decode/celt/bands.rs index bbd041c6..5544ef32 100644 --- a/crates/audio/src/opus/decode/celt/bands.rs +++ b/crates/audio/src/opus/decode/celt/bands.rs @@ -13,23 +13,40 @@ pub fn lcg_rand(seed: u32) -> u32 { seed.wrapping_mul(1_664_525).wrapping_add(1_013_904_223) } -/// Decodes a PVQ codeword of `k` pulses in `n` dimensions; returns the sum of squares. -fn decode_pulses(y: &mut [i32], n: usize, k: usize, dec: &mut RangeDecoder, v: &mut Vec) -> f32 { - // Table of V(nn, kk) for nn in 0..=n, kk in 0..=k (saturating). - let w = k + 1; - v.clear(); - v.resize((n + 1) * w, 0); - v[0] = 1; - for nn in 1..=n { - v[nn * w] = 1; - for kk in 1..=k { - let a = v[(nn - 1) * w + kk]; - let b = v[nn * w + kk - 1]; - let c = v[(nn - 1) * w + kk - 1]; - v[nn * w + kk] = a.saturating_add(b).saturating_add(c); +/// Largest band partition (n) and pulse count (k) that can be coded: the last band of a 20 ms +/// frame has 176 bins, and `get_pulses` yields at most 128 pulses. +const PVQ_MAX_N: usize = 176; +const PVQ_MAX_K: usize = 128; +const PVQ_W: usize = PVQ_MAX_K + 1; + +/// V(n, k) for n <= PVQ_MAX_N, k <= PVQ_MAX_K, saturated at `u32::MAX`, built once. +/// +/// Every (n, k) the bitstream can code has V(n, k) < 2^32 (see `rate::fits_in32`), and the +/// decoder only reads entries V(n', k') <= V(n, k), so the saturation is never observable. +fn pvq_table() -> &'static [u32] { + static TABLE: std::sync::OnceLock> = std::sync::OnceLock::new(); + TABLE.get_or_init(|| { + let mut v = vec![0u32; (PVQ_MAX_N + 1) * PVQ_W]; + v[0] = 1; + for n in 1..=PVQ_MAX_N { + v[n * PVQ_W] = 1; + for k in 1..PVQ_W { + let a = v[(n - 1) * PVQ_W + k]; + let b = v[n * PVQ_W + k - 1]; + let c = v[(n - 1) * PVQ_W + k - 1]; + v[n * PVQ_W + k] = a.saturating_add(b).saturating_add(c); + } } - } - let total = v[n * w + k]; + v + }) +} + +/// Decodes a PVQ codeword of `k` pulses in `n` dimensions; returns the sum of squares. +fn decode_pulses(y: &mut [i32], n: usize, k: usize, dec: &mut RangeDecoder) -> f32 { + debug_assert!(n <= PVQ_MAX_N && k <= PVQ_MAX_K); + let table = pvq_table(); + let vv = |nn: usize, kx: usize| table[nn * PVQ_W + kx] as u64; + let total = vv(n, k); let ft = total.min(u32::MAX as u64) as u32; let mut i = if ft >= 2 { dec.uint(ft) as u64 } else { 0 }; let mut kk = k; @@ -40,7 +57,6 @@ fn decode_pulses(y: &mut [i32], n: usize, k: usize, dec: &mut RangeDecoder, v: & y[j] = 0; continue; } - let vv = |nn: usize, kx: usize| v[nn * w + kx]; let mut p = (vv(rem - 1, kk) + vv(rem, kk)) / 2; let neg = i >= p; if neg { @@ -71,7 +87,7 @@ fn exp_rotation1(x: &mut [f32], len: usize, stride: usize, c: f32, s: f32) { x[i + stride] = c * x2 + s * x1; x[i] = c * x1 + ms * x2; } - if len >= 2 * stride + 1 { + if len > 2 * stride { for i in (0..=len - 2 * stride - 1).rev() { let x1 = x[i]; let x2 = x[i + stride]; @@ -255,7 +271,6 @@ pub struct BandCtx<'m, 'd, 'b> { pub seed: u32, pub disable_inv: bool, pub avoid_split_noise: bool, - pub pvq_v: Vec, pub iy: Vec, pub tmp: Vec, } @@ -269,8 +284,48 @@ struct Split { qalloc: i32, } +#[derive(Clone, Copy)] +struct ThetaBand { + bb: usize, + b0: usize, + lm: i32, + stereo: bool, +} + +struct QuantPartitionInput<'a> { + samples: usize, + bits: i32, + blocks: usize, + lowband: Option<&'a [f32]>, + lm: i32, + gain: f32, + fill: u32, +} + +struct QuantBandInput<'a> { + samples: usize, + bits: i32, + blocks: usize, + lowband: Option<&'a mut [f32]>, + lm: i32, + lowband_out: Option<&'a mut [f32]>, + gain: f32, + fill: u32, +} + +struct QuantBandStereoInput<'a> { + samples: usize, + bits: i32, + blocks: usize, + lowband: Option<&'a mut [f32]>, + lm: i32, + lowband_out: Option<&'a mut [f32]>, + fill: u32, +} + impl BandCtx<'_, '_, '_> { - fn compute_theta(&mut self, n: i32, b: &mut i32, bb: usize, b0: usize, lm: i32, stereo: bool, fill: &mut u32) -> Split { + fn compute_theta(&mut self, n: i32, b: &mut i32, band: ThetaBand, fill: &mut u32) -> Split { + let ThetaBand { bb, b0, lm, stereo } = band; let i = self.i; let pulse_cap = self.m.log_n[i] + lm * (1 << BITRES); let offset = (pulse_cap >> 1) - if stereo && n == 2 { QTHETA_OFFSET_TWOPHASE } else { QTHETA_OFFSET }; @@ -301,13 +356,13 @@ impl BandCtx<'_, '_, '_> { if fm < (half * (half + 1)) >> 1 { let it = (isqrt32(8 * fm + 1) - 1) >> 1; fs = it + 1; - fl = it * (it + 1) >> 1; + fl = (it * (it + 1)) >> 1; itheta = it as i32; } else { let q = qn as u32; let it = (2 * (q + 1) - isqrt32(8 * (ft - fm - 1) + 1)) >> 1; fs = q + 1 - it; - fl = ft - ((q + 1 - it) * (q + 2 - it) >> 1); + fl = ft - (((q + 1 - it) * (q + 2 - it)) >> 1); itheta = it as i32; } self.dec.update(fl, fl + fs, ft); @@ -362,18 +417,8 @@ impl BandCtx<'_, '_, '_> { 1 } - #[allow(clippy::too_many_arguments)] - fn quant_partition( - &mut self, - x: &mut [f32], - n: usize, - mut b: i32, - mut bb: usize, - lowband: Option<&[f32]>, - mut lm: i32, - gain: f32, - mut fill: u32, - ) -> u32 { + fn quant_partition(&mut self, x: &mut [f32], input: QuantPartitionInput<'_>) -> u32 { + let QuantPartitionInput { samples: n, bits: mut b, blocks: mut bb, lowband, mut lm, gain, mut fill } = input; let b0 = bb; let i = self.i; let cache = self.m.cache(lm, i); @@ -386,7 +431,7 @@ impl BandCtx<'_, '_, '_> { fill = (fill & 1) | (fill << 1); } bb = (bb + 1) >> 1; - let sp = self.compute_theta(n as i32, &mut b, bb, b0, lm, false, &mut fill); + let sp = self.compute_theta(n as i32, &mut b, ThetaBand { bb, b0, lm, stereo: false }, &mut fill); let mid = sp.imid as f32 / 32768.0; let side = sp.iside as f32 / 32768.0; let mut delta = sp.delta; @@ -404,19 +449,25 @@ impl BandCtx<'_, '_, '_> { let mut rebalance = self.remaining_bits; let mut cm; if mbits >= sbits { - cm = self.quant_partition(xa, n, mbits, bb, lowband, lm, gain * mid, fill); + cm = self.quant_partition(xa, QuantPartitionInput { samples: n, bits: mbits, blocks: bb, lowband, lm, gain: gain * mid, fill }); rebalance = mbits - (rebalance - self.remaining_bits); if rebalance > 3 << BITRES && sp.itheta != 0 { sbits += rebalance - (3 << BITRES); } - cm |= self.quant_partition(ya, n, sbits, bb, next_lowband2, lm, gain * side, fill >> bb) << (b0 >> 1); + cm |= self.quant_partition( + ya, + QuantPartitionInput { samples: n, bits: sbits, blocks: bb, lowband: next_lowband2, lm, gain: gain * side, fill: fill >> bb }, + ) << (b0 >> 1); } else { - cm = self.quant_partition(ya, n, sbits, bb, next_lowband2, lm, gain * side, fill >> bb) << (b0 >> 1); + cm = self.quant_partition( + ya, + QuantPartitionInput { samples: n, bits: sbits, blocks: bb, lowband: next_lowband2, lm, gain: gain * side, fill: fill >> bb }, + ) << (b0 >> 1); rebalance = sbits - (rebalance - self.remaining_bits); if rebalance > 3 << BITRES && sp.itheta != 16384 { mbits += rebalance - (3 << BITRES); } - cm |= self.quant_partition(xa, n, mbits, bb, lowband, lm, gain * mid, fill); + cm |= self.quant_partition(xa, QuantPartitionInput { samples: n, bits: mbits, blocks: bb, lowband, lm, gain: gain * mid, fill }); } return cm; } @@ -436,7 +487,7 @@ impl BandCtx<'_, '_, '_> { let mut iy = std::mem::take(&mut self.iy); iy.clear(); iy.resize(n, 0); - let ryy = decode_pulses(&mut iy, n, k, self.dec, &mut self.pvq_v); + let ryy = decode_pulses(&mut iy, n, k, self.dec); let g = gain / ryy.sqrt(); for (xv, &yv) in x[..n].iter_mut().zip(iy.iter()) { *xv = g * yv as f32; @@ -475,19 +526,8 @@ impl BandCtx<'_, '_, '_> { } } - #[allow(clippy::too_many_arguments)] - pub fn quant_band( - &mut self, - x: &mut [f32], - n: usize, - b: i32, - mut bb: usize, - mut lowband: Option<&mut [f32]>, - lm: i32, - lowband_out: Option<&mut [f32]>, - gain: f32, - mut fill: u32, - ) -> u32 { + fn quant_band(&mut self, x: &mut [f32], input: QuantBandInput<'_>) -> u32 { + let QuantBandInput { samples: n, bits: b, blocks: mut bb, mut lowband, lm, lowband_out, gain, mut fill } = input; let n0 = n; let mut n_b = n; let b0_orig = bb; @@ -529,7 +569,7 @@ impl BandCtx<'_, '_, '_> { { deinterleave_hadamard(lb, n_b >> recombine, b0 << recombine, long_blocks, &mut tmp); } - let mut cm = self.quant_partition(x, n, b, bb, lowband.as_deref(), lm, gain, fill); + let mut cm = self.quant_partition(x, QuantPartitionInput { samples: n, bits: b, blocks: bb, lowband: lowband.as_deref(), lm, gain, fill }); if b0 > 1 { interleave_hadamard(x, n_b >> recombine, b0 << recombine, long_blocks, &mut tmp); } @@ -557,24 +597,13 @@ impl BandCtx<'_, '_, '_> { cm & ((1u32 << bb) - 1) } - #[allow(clippy::too_many_arguments)] - pub fn quant_band_stereo( - &mut self, - x: &mut [f32], - y: &mut [f32], - n: usize, - mut b: i32, - bb: usize, - lowband: Option<&mut [f32]>, - lm: i32, - lowband_out: Option<&mut [f32]>, - mut fill: u32, - ) -> u32 { + fn quant_band_stereo(&mut self, x: &mut [f32], y: &mut [f32], input: QuantBandStereoInput<'_>) -> u32 { + let QuantBandStereoInput { samples: n, bits: mut b, blocks: bb, lowband, lm, lowband_out, mut fill } = input; let orig_fill = fill; if n == 1 { return self.quant_band_n1(x, Some(y), lowband_out); } - let sp = self.compute_theta(n as i32, &mut b, bb, bb, lm, true, &mut fill); + let sp = self.compute_theta(n as i32, &mut b, ThetaBand { bb, b0: bb, lm, stereo: true }, &mut fill); let mid = sp.imid as f32 / 32768.0; let side = sp.iside as f32 / 32768.0; let mut cm; @@ -593,11 +622,11 @@ impl BandCtx<'_, '_, '_> { } let sign = (1 - 2 * sign) as f32; if c { - cm = self.quant_band(y, n, mbits, bb, lowband, lm, lowband_out, 1.0, orig_fill); + cm = self.quant_band(y, QuantBandInput { samples: n, bits: mbits, blocks: bb, lowband, lm, lowband_out, gain: 1.0, fill: orig_fill }); x[0] = -sign * y[1]; x[1] = sign * y[0]; } else { - cm = self.quant_band(x, n, mbits, bb, lowband, lm, lowband_out, 1.0, orig_fill); + cm = self.quant_band(x, QuantBandInput { samples: n, bits: mbits, blocks: bb, lowband, lm, lowband_out, gain: 1.0, fill: orig_fill }); y[0] = -sign * x[1]; y[1] = sign * x[0]; } @@ -617,19 +646,25 @@ impl BandCtx<'_, '_, '_> { self.remaining_bits -= sp.qalloc; let mut rebalance = self.remaining_bits; if mbits >= sbits { - cm = self.quant_band(x, n, mbits, bb, lowband, lm, lowband_out, 1.0, fill); + cm = self.quant_band(x, QuantBandInput { samples: n, bits: mbits, blocks: bb, lowband, lm, lowband_out, gain: 1.0, fill }); rebalance = mbits - (rebalance - self.remaining_bits); if rebalance > 3 << BITRES && sp.itheta != 0 { sbits += rebalance - (3 << BITRES); } - cm |= self.quant_band(y, n, sbits, bb, None, lm, None, side, fill >> bb); + cm |= self.quant_band( + y, + QuantBandInput { samples: n, bits: sbits, blocks: bb, lowband: None, lm, lowband_out: None, gain: side, fill: fill >> bb }, + ); } else { - cm = self.quant_band(y, n, sbits, bb, None, lm, None, side, fill >> bb); + cm = self.quant_band( + y, + QuantBandInput { samples: n, bits: sbits, blocks: bb, lowband: None, lm, lowband_out: None, gain: side, fill: fill >> bb }, + ); rebalance = sbits - (rebalance - self.remaining_bits); if rebalance > 3 << BITRES && sp.itheta != 16384 { mbits += rebalance - (3 << BITRES); } - cm |= self.quant_band(x, n, mbits, bb, lowband, lm, lowband_out, 1.0, fill); + cm |= self.quant_band(x, QuantBandInput { samples: n, bits: mbits, blocks: bb, lowband, lm, lowband_out, gain: 1.0, fill }); } } if n != 2 { @@ -670,28 +705,50 @@ fn stereo_merge(x: &mut [f32], y: &mut [f32], mid: f32) { } /// Decodes all band shapes (`quant_all_bands`, decoder side). -#[allow(clippy::too_many_arguments)] -pub fn quant_all_bands( - m: &Mode, - start: usize, - end: usize, - x_buf: &mut [f32], - y_buf: Option<&mut [f32]>, - collapse_masks: &mut [u32], - pulses: &[i32; NB_EBANDS], - short_blocks: bool, - spread: usize, - mut dual_stereo: bool, - intensity: usize, - tf_res: &[i32; NB_EBANDS], - total_bits: i32, - mut balance: i32, - dec: &mut RangeDecoder, - lm: usize, - coded_bands: usize, - seed: &mut u32, - disable_inv: bool, -) { +pub struct QuantAllBandsInput<'a, 'd> { + pub mode: &'a Mode, + pub start: usize, + pub end: usize, + pub x_buf: &'a mut [f32], + pub y_buf: Option<&'a mut [f32]>, + pub collapse_masks: &'a mut [u32], + pub pulses: &'a [i32; NB_EBANDS], + pub short_blocks: bool, + pub spread: usize, + pub dual_stereo: bool, + pub intensity: usize, + pub tf_res: &'a [i32; NB_EBANDS], + pub total_bits: i32, + pub balance: i32, + pub decoder: &'a mut RangeDecoder<'d>, + pub lm: usize, + pub coded_bands: usize, + pub seed: &'a mut u32, + pub disable_inv: bool, +} + +pub fn quant_all_bands(input: QuantAllBandsInput<'_, '_>) { + let QuantAllBandsInput { + mode: m, + start, + end, + x_buf, + y_buf, + collapse_masks, + pulses, + short_blocks, + spread, + mut dual_stereo, + intensity, + tf_res, + total_bits, + mut balance, + decoder: dec, + lm, + coded_bands, + seed, + disable_inv, + } = input; let mm = 1usize << lm; let bb = if short_blocks { mm } else { 1 }; let norm_offset = mm * EBANDS[start] as usize; @@ -715,7 +772,6 @@ pub fn quant_all_bands( seed: *seed, disable_inv, avoid_split_noise: bb > 1, - pvq_v: Vec::new(), iy: Vec::new(), tmp: Vec::new(), }; @@ -800,13 +856,19 @@ pub fn quant_all_bands( &mut scratch[..n] }); let lo = if last { None } else { Some(&mut norm[out_off..out_off + n]) }; - x_cm = ctx.quant_band(x, n, b / 2, bb, lb, lm as i32, lo, 1.0, x_cm); + x_cm = ctx.quant_band( + x, + QuantBandInput { samples: n, bits: b / 2, blocks: bb, lowband: lb, lm: lm as i32, lowband_out: lo, gain: 1.0, fill: x_cm }, + ); let lb2 = effective_lowband.map(|e| { scratch2[..n].copy_from_slice(&norm2[e..e + n]); &mut scratch2[..n] }); let lo2 = if last { None } else { Some(&mut norm2[out_off..out_off + n]) }; - y_cm = ctx.quant_band(y, n, b / 2, bb, lb2, lm as i32, lo2, 1.0, y_cm); + y_cm = ctx.quant_band( + y, + QuantBandInput { samples: n, bits: b / 2, blocks: bb, lowband: lb2, lm: lm as i32, lowband_out: lo2, gain: 1.0, fill: y_cm }, + ); } else { let lb = effective_lowband.map(|e| { scratch[..n].copy_from_slice(&norm[e..e + n]); @@ -815,9 +877,16 @@ pub fn quant_all_bands( let lo = if last { None } else { Some(&mut norm[out_off..out_off + n]) }; if let Some(yb) = y_buf.as_deref_mut() { let y = &mut yb[xo..xo + n]; - x_cm = ctx.quant_band_stereo(x, y, n, b, bb, lb, lm as i32, lo, x_cm | y_cm); + x_cm = ctx.quant_band_stereo( + x, + y, + QuantBandStereoInput { samples: n, bits: b, blocks: bb, lowband: lb, lm: lm as i32, lowband_out: lo, fill: x_cm | y_cm }, + ); } else { - x_cm = ctx.quant_band(x, n, b, bb, lb, lm as i32, lo, 1.0, x_cm | y_cm); + x_cm = ctx.quant_band( + x, + QuantBandInput { samples: n, bits: b, blocks: bb, lowband: lb, lm: lm as i32, lowband_out: lo, gain: 1.0, fill: x_cm | y_cm }, + ); } y_cm = x_cm; } @@ -831,21 +900,23 @@ pub fn quant_all_bands( } /// Anti-collapse noise injection (RFC 6716 §4.3.5). -#[allow(clippy::too_many_arguments)] -pub fn anti_collapse( - x_buf: &mut [f32], - collapse_masks: &[u32], - lm: usize, - c: usize, - size: usize, - start: usize, - end: usize, - log_e: &[f32], - prev1: &[f32], - prev2: &[f32], - pulses: &[i32; NB_EBANDS], - mut seed: u32, -) { +pub struct AntiCollapseInput<'a> { + pub x_buf: &'a mut [f32], + pub collapse_masks: &'a [u32], + pub lm: usize, + pub channels: usize, + pub size: usize, + pub start: usize, + pub end: usize, + pub log_e: &'a [f32], + pub prev1: &'a [f32], + pub prev2: &'a [f32], + pub pulses: &'a [i32; NB_EBANDS], + pub seed: u32, +} + +pub fn anti_collapse(input: AntiCollapseInput<'_>) { + let AntiCollapseInput { x_buf, collapse_masks, lm, channels: c, size, start, end, log_e, prev1, prev2, pulses, mut seed } = input; for i in start..end { let n0 = (EBANDS[i + 1] - EBANDS[i]) as usize; let depth = ((1 + pulses[i]) as u32 / n0 as u32) >> lm; @@ -885,17 +956,19 @@ pub fn anti_collapse( } /// Multiplies the normalised shapes by the band energies (`denormalise_bands`). -#[allow(clippy::too_many_arguments)] -pub fn denormalise_bands( - x: &[f32], - freq: &mut [f32], - band_log_e: &[f32], - mut start: usize, - mut end: usize, - mm: usize, - downsample: usize, - silence: bool, -) { +pub struct DenormaliseBandsInput<'a> { + pub x: &'a [f32], + pub freq: &'a mut [f32], + pub band_log_e: &'a [f32], + pub start: usize, + pub end: usize, + pub mm: usize, + pub downsample: usize, + pub silence: bool, +} + +pub fn denormalise_bands(input: DenormaliseBandsInput<'_>) { + let DenormaliseBandsInput { x, freq, band_log_e, mut start, mut end, mm, downsample, silence } = input; let n = mm * SHORT_MDCT; let mut bound = mm * EBANDS[end] as usize; if downsample != 1 { @@ -944,6 +1017,18 @@ mod tests { } } + #[test] + fn pvq_table_matches_reference_counts() { + let table = pvq_table(); + for n in 0..=PVQ_MAX_N { + for k in 0..=PVQ_MAX_K { + let want = super::super::rate::pvq_v(n, k).min(u32::MAX as u128) as u32; + // V(0, k > 0) = 0 and V(n, 0) = 1 in both. + assert_eq!(table[n * PVQ_W + k], want, "n={n} k={k}"); + } + } + } + #[test] fn pvq_roundtrip_all_codewords_small() { use crate::opus::decode::range::enc::RangeEncoder; @@ -958,8 +1043,7 @@ mod tests { let buf = e.done(); let mut d = RangeDecoder::new(&buf); let mut y = vec![0i32; n]; - let mut v = Vec::new(); - let ryy = decode_pulses(&mut y, n, k, &mut d, &mut v); + let ryy = decode_pulses(&mut y, n, k, &mut d); assert_eq!(y.iter().map(|v| v.unsigned_abs() as usize).sum::(), k); assert_eq!(ryy as i32, y.iter().map(|v| v * v).sum::()); assert!(seen.insert(y.clone()), "duplicate codeword n={n} k={k} idx={idx}"); diff --git a/crates/audio/src/opus/decode/celt/energy.rs b/crates/audio/src/opus/decode/celt/energy.rs index f3275614..fccb767e 100644 --- a/crates/audio/src/opus/decode/celt/energy.rs +++ b/crates/audio/src/opus/decode/celt/energy.rs @@ -44,7 +44,6 @@ pub fn laplace_decode(dec: &mut RangeDecoder, mut fs: u32, decay: u32) -> i32 { } /// `unquant_coarse_energy`. -#[allow(clippy::too_many_arguments)] pub fn unquant_coarse_energy(start: usize, end: usize, old_e: &mut [f32], intra: bool, dec: &mut RangeDecoder, c: usize, lm: usize) { let prob = &E_PROB_MODEL[lm][intra as usize]; let (coef, beta) = if intra { (0.0, BETA_INTRA) } else { (PRED_COEF[lm], BETA_COEF[lm]) }; @@ -88,17 +87,19 @@ pub fn unquant_fine_energy(start: usize, end: usize, old_e: &mut [f32], fine_qua } } -#[allow(clippy::too_many_arguments)] -pub fn unquant_energy_finalise( - start: usize, - end: usize, - old_e: &mut [f32], - fine_quant: &[i32; NB_EBANDS], - fine_priority: &[i32; NB_EBANDS], - mut bits_left: i32, - dec: &mut RangeDecoder, - c: usize, -) { +pub struct EnergyFinalise<'a, 'd> { + pub start: usize, + pub end: usize, + pub old_e: &'a mut [f32], + pub fine_quant: &'a [i32; NB_EBANDS], + pub fine_priority: &'a [i32; NB_EBANDS], + pub bits_left: i32, + pub decoder: &'a mut RangeDecoder<'d>, + pub channels: usize, +} + +pub fn unquant_energy_finalise(input: EnergyFinalise<'_, '_>) { + let EnergyFinalise { start, end, old_e, fine_quant, fine_priority, mut bits_left, decoder: dec, channels: c } = input; for prio in 0..2 { let mut i = start; while i < end && bits_left >= c as i32 { diff --git a/crates/audio/src/opus/decode/celt/mdct.rs b/crates/audio/src/opus/decode/celt/mdct.rs index 13dc4a36..b3bfcf0f 100644 --- a/crates/audio/src/opus/decode/celt/mdct.rs +++ b/crates/audio/src/opus/decode/celt/mdct.rs @@ -19,9 +19,14 @@ impl Cpx { fn add(self, o: Cpx) -> Cpx { Cpx { re: self.re + o.re, im: self.im + o.im } } + #[inline] + fn sub(self, o: Cpx) -> Cpx { + Cpx { re: self.re - o.re, im: self.im - o.im } + } } -/// Mixed-radix (2, 3, 4, 5) forward complex FFT, unscaled. +/// Mixed-radix (2, 3, 4, 5) forward complex FFT, unscaled. Radix 2, 3 and 4 use dedicated +/// butterflies; radix 5 uses a generic one. pub struct Fft { n: usize, factors: Vec, @@ -66,20 +71,70 @@ impl Fft { self.rec(&mut out[q * m..(q + 1) * m], &inp[q * in_stride..], in_stride * p, m, fi + 1, tw_stride * p); } } - let big = self.n; - let mut t = [Cpx::default(); 5]; - for k in 0..m { - for q in 0..p { - let w = self.twiddle[(q * k * tw_stride) % big]; - t[q] = out[q * m + k].mul(w); + let tw = &self.twiddle; + match p { + 2 => { + for k in 0..m { + let a0 = out[k]; + let a1 = out[m + k].mul(tw[k * tw_stride]); + out[k] = a0.add(a1); + out[m + k] = a0.sub(a1); + } } - for q2 in 0..p { - let mut acc = Cpx::default(); - for (q, tq) in t.iter().enumerate().take(p) { - let w = self.twiddle[((q * q2) % p) * (big / p)]; - acc = acc.add(tq.mul(w)); + 4 => { + for k in 0..m { + let t0 = out[k]; + let t1 = out[m + k].mul(tw[k * tw_stride]); + let t2 = out[2 * m + k].mul(tw[2 * k * tw_stride]); + let t3 = out[3 * m + k].mul(tw[3 * k * tw_stride]); + let s0 = t0.add(t2); + let d0 = t0.sub(t2); + let s1 = t1.add(t3); + let d1 = t1.sub(t3); + out[k] = s0.add(s1); + out[2 * m + k] = s0.sub(s1); + // -j * d1 = (d1.im, -d1.re) + out[m + k] = Cpx { re: d0.re + d1.im, im: d0.im - d1.re }; + out[3 * m + k] = Cpx { re: d0.re - d1.im, im: d0.im + d1.re }; + } + } + 3 => { + // sin(-2*pi/3) = -sqrt(3)/2 + let h = tw[self.n / 3].im; + for k in 0..m { + let t0 = out[k]; + let t1 = out[m + k].mul(tw[k * tw_stride]); + let t2 = out[2 * m + k].mul(tw[2 * k * tw_stride]); + let s = t1.add(t2); + let d = t1.sub(t2); + let a = Cpx { re: t0.re - 0.5 * s.re, im: t0.im - 0.5 * s.im }; + out[k] = t0.add(s); + // j * h * d = (-h * d.im, h * d.re) + out[m + k] = Cpx { re: a.re - h * d.im, im: a.im + h * d.re }; + out[2 * m + k] = Cpx { re: a.re + h * d.im, im: a.im - h * d.re }; + } + } + _ => { + // Radix 5 (at most once per transform): generic butterfly with the five roots of + // unity looked up once. + let mut w5 = [Cpx::default(); 5]; + for (q, w) in w5.iter_mut().enumerate() { + *w = tw[q * (self.n / p)]; + } + let mut t = [Cpx::default(); 5]; + for k in 0..m { + t[0] = out[k]; + for q in 1..p { + t[q] = out[q * m + k].mul(tw[q * k * tw_stride]); + } + for q2 in 0..p { + let mut acc = t[0]; + for q in 1..p { + acc = acc.add(t[q].mul(w5[(q * q2) % p])); + } + out[q2 * m + k] = acc; + } } - out[q2 * m + k] = acc; } } } @@ -159,7 +214,7 @@ mod tests { #[test] fn fft_matches_dft() { - for &n in &[60usize, 120, 240, 480] { + for &n in &[15usize, 30, 60, 120, 240, 480, 960] { let fft = Fft::new(n); let x: Vec = (0..n).map(|i| Cpx { re: ((i * 7) % 13) as f32 - 6.0, im: ((i * 3) % 5) as f32 - 2.0 }).collect(); let mut y = x.clone(); @@ -176,6 +231,63 @@ mod tests { } } + /// The previous generic implementation (every butterfly multiplies by every DFT root), kept + /// to measure how far the specialised kernels move the output. + fn reference_fft(n: usize, data: &mut [Cpx]) { + let fft = Fft::new(n); + let inp = data.to_vec(); + rec_ref(&fft, data, &inp, 1, n, 0, 1); + } + + fn rec_ref(f: &Fft, out: &mut [Cpx], inp: &[Cpx], in_stride: usize, n: usize, fi: usize, tw_stride: usize) { + let p = f.factors[fi]; + let m = n / p; + if m == 1 { + for (i, o) in out.iter_mut().enumerate().take(p) { + *o = inp[i * in_stride]; + } + } else { + for q in 0..p { + rec_ref(f, &mut out[q * m..(q + 1) * m], &inp[q * in_stride..], in_stride * p, m, fi + 1, tw_stride * p); + } + } + let big = f.n; + let mut t = [Cpx::default(); 5]; + for k in 0..m { + for q in 0..p { + let w = f.twiddle[(q * k * tw_stride) % big]; + t[q] = out[q * m + k].mul(w); + } + for q2 in 0..p { + let mut acc = Cpx::default(); + for (q, tq) in t.iter().enumerate().take(p) { + let w = f.twiddle[((q * q2) % p) * (big / p)]; + acc = acc.add(tq.mul(w)); + } + out[q2 * m + k] = acc; + } + } + } + + /// Every size the CELT frame sizes use (N/4 for N = 2 * {120, 240, 480, 960}, plus the + /// short-block sizes) agrees with the previous implementation to f32 rounding. + #[test] + fn fft_matches_reference_implementation() { + for &n in &[15usize, 30, 60, 120, 240, 480, 960] { + let fft = Fft::new(n); + let x: Vec = (0..n) + .map(|i| Cpx { re: (((i * 7919) % 2003) as f32 - 1000.0) / 500.0, im: (((i * 104729) % 1009) as f32 - 500.0) / 250.0 }) + .collect(); + let mut new = x.clone(); + fft.process(&mut new, &mut Vec::new()); + let mut old = x.clone(); + reference_fft(n, &mut old); + let peak = old.iter().map(|c| c.re.abs().max(c.im.abs())).fold(0f32, f32::max); + let max_diff = new.iter().zip(&old).map(|(a, b)| (a.re - b.re).abs().max((a.im - b.im).abs())).fold(0f32, f32::max); + assert!(max_diff <= 1e-5 * peak, "n={n} max_diff={max_diff} peak={peak}"); + } + } + /// The folded output must equal the direct IMDCT formula folded by TDAC symmetry. #[test] fn imdct_matches_direct_formula() { diff --git a/crates/audio/src/opus/decode/celt/mod.rs b/crates/audio/src/opus/decode/celt/mod.rs index fa002646..51b04512 100644 --- a/crates/audio/src/opus/decode/celt/mod.rs +++ b/crates/audio/src/opus/decode/celt/mod.rs @@ -8,10 +8,10 @@ pub(crate) mod tables; use std::sync::OnceLock; -use bands::{anti_collapse, denormalise_bands, quant_all_bands, renormalise_vector}; -use energy::{unquant_coarse_energy, unquant_energy_finalise, unquant_fine_energy}; +use bands::{AntiCollapseInput, DenormaliseBandsInput, QuantAllBandsInput, anti_collapse, denormalise_bands, quant_all_bands, renormalise_vector}; +use energy::{EnergyFinalise, unquant_coarse_energy, unquant_energy_finalise, unquant_fine_energy}; use mdct::ImdctScratch; -use rate::{Mode, compute_allocation}; +use rate::{AllocationInput, Mode, compute_allocation}; use tables::{COMB_GAINS, EBANDS, MAX_LM, NB_EBANDS, OVERLAP, PREEMPH, SHORT_MDCT, SPREAD_ICDF, TAPSET_ICDF, TF_SELECT_TABLE, TRIM_ICDF}; use crate::opus::decode::range::{BITRES, RangeDecoder}; @@ -21,11 +21,58 @@ const COMBFILTER_MINPERIOD: i32 = 15; const PLC_PITCH_LAG_MAX: usize = 720; const PLC_PITCH_LAG_MIN: usize = 100; +struct SynthesisFrame<'a> { + bands: &'a [f32], + samples: usize, + start_band: usize, + end_band: usize, + coded_channels: usize, + transient: bool, + lm: usize, + silence: bool, +} + pub(crate) fn mode() -> &'static Mode { static MODE: OnceLock = OnceLock::new(); MODE.get_or_init(Mode::new) } +fn packet_silence(dec: &mut RangeDecoder, total_bits: i32) -> (bool, i32) { + let mut tell = dec.tell(); + let silence = if tell >= total_bits { + true + } else if tell == 1 { + dec.bit_logp(15) + } else { + false + }; + if silence { + // Pretend we've read all the remaining bits. + dec.skip_to_end(total_bits); + tell = total_bits; + } + (silence, tell) +} + +fn read_postfilter(dec: &mut RangeDecoder, start: usize, mut tell: i32, total_bits: i32) -> (f32, i32, usize, i32) { + let mut gain = 0.0; + let mut pitch = 0; + let mut tapset = 0; + if start == 0 && tell + 16 <= total_bits { + if dec.bit_logp(1) { + let octave = dec.uint(6); + pitch = ((16 << octave) + dec.bits(4 + octave)) as i32 - 1; + let qg = dec.bits(3); + if dec.tell() + 2 <= total_bits { + tapset = dec.icdf(&TAPSET_ICDF, 2); + } + gain = 0.09375 * (qg + 1) as f32; + } + tell = dec.tell(); + } + (gain, pitch, tapset, tell) +} + /// CELT decoder state for one (mono or stereo) stream. pub struct CeltDecoder { m: &'static Mode, @@ -181,34 +228,8 @@ impl CeltDecoder { } } let total_bits_raw = len as i32 * 8; - let mut tell = dec.tell(); - let silence = if tell >= total_bits_raw { - true - } else if tell == 1 { - dec.bit_logp(15) - } else { - false - }; - if silence { - // Pretend we've read all the remaining bits. - dec.skip_to_end(total_bits_raw); - tell = total_bits_raw; - } - let mut postfilter_gain = 0.0f32; - let mut postfilter_pitch = 0i32; - let mut postfilter_tapset = 0usize; - if start == 0 && tell + 16 <= total_bits_raw { - if dec.bit_logp(1) { - let octave = dec.uint(6); - postfilter_pitch = ((16 << octave) + dec.bits(4 + octave)) as i32 - 1; - let qg = dec.bits(3); - if dec.tell() + 2 <= total_bits_raw { - postfilter_tapset = dec.icdf(&TAPSET_ICDF, 2); - } - postfilter_gain = 0.09375 * (qg + 1) as f32; - } - tell = dec.tell(); - } + let (silence, tell) = packet_silence(dec, total_bits_raw); + let (postfilter_gain, postfilter_pitch, postfilter_tapset, mut tell) = read_postfilter(dec, start, tell, total_bits_raw); let is_transient = if lm > 0 && tell + 3 <= total_bits_raw { let t = dec.bit_logp(3); tell = dec.tell(); @@ -255,7 +276,18 @@ impl CeltDecoder { let mut bits = ((len as i32 * 8) << BITRES) - dec.tell_frac() as i32 - 1; let anti_collapse_rsv = if is_transient && lm >= 2 && bits >= ((lm as i32 + 2) << BITRES) { 1 << BITRES } else { 0 }; bits -= anti_collapse_rsv; - let alloc = compute_allocation(m, start, end, &offsets, &cap, alloc_trim, bits, c, lm, dec); + let alloc = compute_allocation(AllocationInput { + mode: m, + start, + end, + offsets: &offsets, + cap: &cap, + trim: alloc_trim, + total: bits, + channels: c, + lm, + decoder: dec, + }); unquant_fine_energy(start, end, &mut self.old_band_e, &alloc.fine_quant, dec, c); // Shift the decode memory by one frame. @@ -273,27 +305,27 @@ impl CeltDecoder { let (xa, xb) = x.split_at_mut(n); let yb = if c == 2 { Some(xb) } else { None }; let mut seed = self.rng; - quant_all_bands( - m, + quant_all_bands(QuantAllBandsInput { + mode: m, start, end, - xa, - yb, - &mut collapse_masks, - &alloc.pulses, + x_buf: xa, + y_buf: yb, + collapse_masks: &mut collapse_masks, + pulses: &alloc.pulses, short_blocks, spread, - alloc.dual_stereo, - alloc.intensity, - &tf_res, - ((len as i32) * (8 << BITRES)) - anti_collapse_rsv, - alloc.balance, - dec, + dual_stereo: alloc.dual_stereo, + intensity: alloc.intensity, + tf_res: &tf_res, + total_bits: ((len as i32) * (8 << BITRES)) - anti_collapse_rsv, + balance: alloc.balance, + decoder: dec, lm, - alloc.coded_bands, - &mut seed, - self.disable_inv, - ); + coded_bands: alloc.coded_bands, + seed: &mut seed, + disable_inv: self.disable_inv, + }); self.rng = seed; } let mut anti_collapse_on = false; @@ -301,29 +333,47 @@ impl CeltDecoder { anti_collapse_on = dec.bits(1) != 0; } let bits_left = len as i32 * 8 - dec.tell(); - unquant_energy_finalise(start, end, &mut self.old_band_e, &alloc.fine_quant, &alloc.fine_priority, bits_left, dec, c); + unquant_energy_finalise(EnergyFinalise { + start, + end, + old_e: &mut self.old_band_e, + fine_quant: &alloc.fine_quant, + fine_priority: &alloc.fine_priority, + bits_left, + decoder: dec, + channels: c, + }); if anti_collapse_on { - anti_collapse( - &mut x, - &collapse_masks, + anti_collapse(AntiCollapseInput { + x_buf: &mut x, + collapse_masks: &collapse_masks, lm, - c, - n, + channels: c, + size: n, start, end, - &self.old_band_e, - &self.old_log_e, - &self.old_log_e2, - &alloc.pulses, - self.rng, - ); + log_e: &self.old_band_e, + prev1: &self.old_log_e, + prev2: &self.old_log_e2, + pulses: &alloc.pulses, + seed: self.rng, + }); } if silence { for v in self.old_band_e[..c * NB_EBANDS].iter_mut() { *v = -28.0; } } - self.synthesis(&x, n, start, eff_end, c, is_transient, lm, silence); + self.synthesis(SynthesisFrame { + bands: &x, + samples: n, + start_band: start, + end_band: eff_end, + coded_channels: c, + transient: is_transient, + lm, + silence, + }); self.x = x; // Post-filter. @@ -334,28 +384,32 @@ impl CeltDecoder { let off = DECODE_BUFFER_SIZE - n; comb_filter( mem, - off, - self.postfilter_period_old, - self.postfilter_period, - SHORT_MDCT, - self.postfilter_gain_old, - self.postfilter_gain, - self.postfilter_tapset_old, - self.postfilter_tapset, - &m.window, + CombFilterParams { + offset: off, + old_period: self.postfilter_period_old, + new_period: self.postfilter_period, + samples: SHORT_MDCT, + old_gain: self.postfilter_gain_old, + new_gain: self.postfilter_gain, + old_tapset: self.postfilter_tapset_old, + new_tapset: self.postfilter_tapset, + window: &m.window, + }, ); if lm != 0 { comb_filter( mem, - off + SHORT_MDCT, - self.postfilter_period, - postfilter_pitch, - n - SHORT_MDCT, - self.postfilter_gain, - postfilter_gain, - self.postfilter_tapset, - postfilter_tapset, - &m.window, + CombFilterParams { + offset: off + SHORT_MDCT, + old_period: self.postfilter_period, + new_period: postfilter_pitch, + samples: n - SHORT_MDCT, + old_gain: self.postfilter_gain, + new_gain: postfilter_gain, + old_tapset: self.postfilter_tapset, + new_tapset: postfilter_tapset, + window: &m.window, + }, ); } } @@ -398,8 +452,9 @@ impl CeltDecoder { } /// Denormalisation + IMDCT into the decode memory (`celt_synthesis`). - #[allow(clippy::too_many_arguments)] - fn synthesis(&mut self, x: &[f32], n: usize, start: usize, eff_end: usize, c: usize, is_transient: bool, lm: usize, silence: bool) { + fn synthesis(&mut self, frame: SynthesisFrame<'_>) { + let SynthesisFrame { bands: x, samples: n, start_band: start, end_band: eff_end, coded_channels: c, transient: is_transient, lm, silence } = + frame; let cc = self.channels; let m = self.m; let mm = 1usize << lm; @@ -411,7 +466,16 @@ impl CeltDecoder { let (f0, f1) = freq.split_at_mut(n); let out_off = DECODE_BUFFER_SIZE - n; if cc == 2 && c == 1 { - denormalise_bands(x, f0, &self.old_band_e, start, eff_end, mm, self.downsample, silence); + denormalise_bands(DenormaliseBandsInput { + x, + freq: f0, + band_log_e: &self.old_band_e, + start, + end: eff_end, + mm, + downsample: self.downsample, + silence, + }); f1.copy_from_slice(f0); for (ch, f) in [(0usize, &*f1), (1usize, &*f0)] { let mem = &mut self.decode_mem[ch * ml..(ch + 1) * ml]; @@ -420,8 +484,26 @@ impl CeltDecoder { } } } else if cc == 1 && c == 2 { - denormalise_bands(x, f0, &self.old_band_e, start, eff_end, mm, self.downsample, silence); - denormalise_bands(&x[n..], f1, &self.old_band_e[NB_EBANDS..], start, eff_end, mm, self.downsample, silence); + denormalise_bands(DenormaliseBandsInput { + x, + freq: f0, + band_log_e: &self.old_band_e, + start, + end: eff_end, + mm, + downsample: self.downsample, + silence, + }); + denormalise_bands(DenormaliseBandsInput { + x: &x[n..], + freq: f1, + band_log_e: &self.old_band_e[NB_EBANDS..], + start, + end: eff_end, + mm, + downsample: self.downsample, + silence, + }); for i in 0..n { f0[i] = 0.5 * (f0[i] + f1[i]); } @@ -431,7 +513,16 @@ impl CeltDecoder { } } else { for ch in 0..cc { - denormalise_bands(&x[ch * n..], f0, &self.old_band_e[ch * NB_EBANDS..], start, eff_end, mm, self.downsample, silence); + denormalise_bands(DenormaliseBandsInput { + x: &x[ch * n..], + freq: f0, + band_log_e: &self.old_band_e[ch * NB_EBANDS..], + start, + end: eff_end, + mm, + downsample: self.downsample, + silence, + }); let mem = &mut self.decode_mem[ch * ml..(ch + 1) * ml]; for b in 0..bcount { m.mdct[shift].backward(&f0[b..], bcount, &mut mem[out_off + nb * b..], &m.window, OVERLAP, &mut self.imdct); @@ -446,6 +537,43 @@ impl CeltDecoder { let ml = self.mem_len(); let ds = self.downsample; let out_off = DECODE_BUFFER_SIZE - n; + if ds == 1 && cc == 2 { + // Both channels in one loop: two independent IIR chains overlap their latency. + let (m0, m1) = self.decode_mem.split_at(ml); + let left = &m0[out_off..out_off + n]; + let right = &m1[out_off..out_off + n]; + let [mut s0, mut s1] = self.preemph_mem; + for ((o, &v0), &v1) in pcm[..2 * n].as_chunks_mut::<2>().0.iter_mut().zip(left).zip(right) { + let t0 = v0 + 1e-30 + s0; + let t1 = v1 + 1e-30 + s1; + s0 = PREEMPH * t0; + s1 = PREEMPH * t1; + if accum { + o[0] += t0 * (1.0 / 32768.0); + o[1] += t1 * (1.0 / 32768.0); + } else { + o[0] = t0 * (1.0 / 32768.0); + o[1] = t1 * (1.0 / 32768.0); + } + } + self.preemph_mem = [s0, s1]; + return; + } + if ds == 1 && cc == 1 { + let mem = &self.decode_mem[out_off..out_off + n]; + let mut state = self.preemph_mem[0]; + for (o, &v) in pcm[..n].iter_mut().zip(mem) { + let t = v + 1e-30 + state; + state = PREEMPH * t; + if accum { + *o += t * (1.0 / 32768.0); + } else { + *o = t * (1.0 / 32768.0); + } + } + self.preemph_mem[0] = state; + return; + } let mut scratch = std::mem::take(&mut self.pcm_scratch); for ch in 0..cc { let mem = &self.decode_mem[ch * ml + out_off..ch * ml + out_off + n]; @@ -508,7 +636,16 @@ impl CeltDecoder { } let saved_sc = self.stream_channels; self.stream_channels = cc; - self.synthesis(&x, n, self.start, eff_end, cc, false, lm, false); + self.synthesis(SynthesisFrame { + bands: &x, + samples: n, + start_band: self.start, + end_band: eff_end, + coded_channels: cc, + transient: false, + lm, + silence: false, + }); self.stream_channels = saved_sc; self.x = x; // As libopus: after noise, the history holds no pitch to extend until two good frames @@ -586,7 +723,7 @@ fn tf_decode(start: usize, end: usize, is_transient: bool, tf_res: &mut [i32; NB let mut budget = dec.storage() as i32 * 8; let mut tell = dec.tell(); let mut logp: u32 = if is_transient { 2 } else { 4 }; - let tf_select_rsv = lm > 0 && tell + logp as i32 + 1 <= budget; + let tf_select_rsv = lm > 0 && tell + (logp as i32) < budget; budget -= tf_select_rsv as i32; let mut tf_changed = 0; let mut curr = 0; @@ -610,8 +747,30 @@ fn tf_decode(start: usize, end: usize, is_transient: bool, tf_res: &mut [i32; NB } /// Pitch post-filter (in place on `buf[off .. off + n]`, reading history before `off`). -#[allow(clippy::too_many_arguments)] -fn comb_filter(buf: &mut [f32], off: usize, t0: i32, t1: i32, n: usize, g0: f32, g1: f32, tapset0: usize, tapset1: usize, window: &[f32; OVERLAP]) { +struct CombFilterParams<'a> { + offset: usize, + old_period: i32, + new_period: i32, + samples: usize, + old_gain: f32, + new_gain: f32, + old_tapset: usize, + new_tapset: usize, + window: &'a [f32; OVERLAP], +} + +fn comb_filter(buf: &mut [f32], params: CombFilterParams<'_>) { + let CombFilterParams { + offset: off, + old_period: t0, + new_period: t1, + samples: n, + old_gain: g0, + new_gain: g1, + old_tapset: tapset0, + new_tapset: tapset1, + window, + } = params; if g0 == 0.0 && g1 == 0.0 { return; } diff --git a/crates/audio/src/opus/decode/celt/rate.rs b/crates/audio/src/opus/decode/celt/rate.rs index cbbca119..782d68d5 100644 --- a/crates/audio/src/opus/decode/celt/rate.rs +++ b/crates/audio/src/opus/decode/celt/rate.rs @@ -142,7 +142,7 @@ impl Mode { let mut cap = [0i32; NB_EBANDS]; for (i, cp) in cap.iter_mut().enumerate() { let n = (EBANDS[i + 1] - EBANDS[i]) << lm; - *cp = (self.caps[NB_EBANDS * (2 * lm + c - 1) + i] as i32 + 64) * c as i32 * n >> 2; + *cp = ((self.caps[NB_EBANDS * (2 * lm + c - 1) + i] as i32 + 64) * c as i32 * n) >> 2; } cap } @@ -217,7 +217,7 @@ fn compute_caps(m: &Mode) -> Vec { for j in 0..nb { let mut n0 = EBANDS[j + 1] - EBANDS[j]; let max_bits = if n0 << i == 1 { - c * (1 + MAX_FINE_BITS) << BITRES + (c * (1 + MAX_FINE_BITS)) << BITRES } else { let mut lm0 = 0; if n0 > 2 { @@ -257,7 +257,7 @@ fn compute_caps(m: &Mode) -> Vec { let num = mb + ndof * offset; let den = (ndof - 1) << BITRES; let qb = ((num + (den >> 1)) / den).min(MAX_FINE_BITS); - mb += c * qb << BITRES; + mb += (c * qb) << BITRES; mb }; let w = c * ((EBANDS[j + 1] - EBANDS[j]) << i); @@ -281,19 +281,21 @@ pub struct Allocation { } /// `clt_compute_allocation` (decoder side). -#[allow(clippy::too_many_arguments)] -pub fn compute_allocation( - m: &Mode, - start: usize, - end: usize, - offsets: &[i32; NB_EBANDS], - cap: &[i32; NB_EBANDS], - alloc_trim: i32, - total: i32, - c: usize, - lm: usize, - dec: &mut RangeDecoder, -) -> Allocation { +pub struct AllocationInput<'a, 'd> { + pub mode: &'a Mode, + pub start: usize, + pub end: usize, + pub offsets: &'a [i32; NB_EBANDS], + pub cap: &'a [i32; NB_EBANDS], + pub trim: i32, + pub total: i32, + pub channels: usize, + pub lm: usize, + pub decoder: &'a mut RangeDecoder<'d>, +} + +pub fn compute_allocation(input: AllocationInput<'_, '_>) -> Allocation { + let AllocationInput { mode: m, start, end, offsets, cap, trim: alloc_trim, total, channels: c, lm, decoder: dec } = input; let ci = c as i32; let lmi = lm as i32; let mut total = total.max(0); @@ -318,8 +320,8 @@ pub fn compute_allocation( let mut trim_offset = [0i32; NB_EBANDS]; for j in start..end { let w = EBANDS[j + 1] - EBANDS[j]; - thresh[j] = (ci << BITRES).max((3 * w << lmi << BITRES) >> 4); - trim_offset[j] = ci * w * (alloc_trim - 5 - lmi) * (end as i32 - j as i32 - 1) * (1 << (lmi + BITRES as i32)) >> 6; + thresh[j] = (ci << BITRES).max(((3 * w) << lmi << BITRES) >> 4); + trim_offset[j] = (ci * w * (alloc_trim - 5 - lmi) * (end as i32 - j as i32 - 1) * (1 << (lmi + BITRES as i32))) >> 6; if w << lmi == 1 { trim_offset[j] -= ci << BITRES; } @@ -333,7 +335,7 @@ pub fn compute_allocation( let mid = (lo + hi) >> 1; for j in (start..end).rev() { let n = EBANDS[j + 1] - EBANDS[j]; - let mut bitsj = ci * n * (ALLOC[mid as usize][j] as i32) << lmi >> 2; + let mut bitsj = (ci * n * (ALLOC[mid as usize][j] as i32)) << lmi >> 2; if bitsj > 0 { bitsj = (bitsj + trim_offset[j]).max(0); } @@ -358,8 +360,8 @@ pub fn compute_allocation( lo -= 1; for j in start..end { let n = EBANDS[j + 1] - EBANDS[j]; - let mut bits1j = ci * n * (ALLOC[lo as usize][j] as i32) << lmi >> 2; - let mut bits2j = if hi >= nb_alloc { cap[j] } else { ci * n * (ALLOC[hi as usize][j] as i32) << lmi >> 2 }; + let mut bits1j = (ci * n * (ALLOC[lo as usize][j] as i32)) << lmi >> 2; + let mut bits2j = if hi >= nb_alloc { cap[j] } else { (ci * n * (ALLOC[hi as usize][j] as i32)) << lmi >> 2 }; if bits1j > 0 { bits1j = (bits1j + trim_offset[j]).max(0); } @@ -377,27 +379,61 @@ pub fn compute_allocation( bits1[j] = bits1j; bits2[j] = bits2j; } - interp_bits2pulses(m, start, end, skip_start, &bits1, &bits2, &thresh, cap, total, skip_rsv, intensity_rsv, dual_stereo_rsv, c, lm, dec) + interp_bits2pulses(InterpPulsesInput { + mode: m, + start, + end, + skip_start, + bits1: &bits1, + bits2: &bits2, + threshold: &thresh, + cap, + total, + skip_reservation: skip_rsv, + intensity_reservation: intensity_rsv, + stereo_reservation: dual_stereo_rsv, + channels: c, + lm, + decoder: dec, + }) } -#[allow(clippy::too_many_arguments)] -fn interp_bits2pulses( - m: &Mode, +struct InterpPulsesInput<'a, 'd> { + mode: &'a Mode, start: usize, end: usize, skip_start: usize, - bits1: &[i32; NB_EBANDS], - bits2: &[i32; NB_EBANDS], - thresh: &[i32; NB_EBANDS], - cap: &[i32; NB_EBANDS], - mut total: i32, - skip_rsv: i32, - mut intensity_rsv: i32, - mut dual_stereo_rsv: i32, - c: usize, + bits1: &'a [i32; NB_EBANDS], + bits2: &'a [i32; NB_EBANDS], + threshold: &'a [i32; NB_EBANDS], + cap: &'a [i32; NB_EBANDS], + total: i32, + skip_reservation: i32, + intensity_reservation: i32, + stereo_reservation: i32, + channels: usize, lm: usize, - dec: &mut RangeDecoder, -) -> Allocation { + decoder: &'a mut RangeDecoder<'d>, +} + +fn interp_bits2pulses(input: InterpPulsesInput<'_, '_>) -> Allocation { + let InterpPulsesInput { + mode: m, + start, + end, + skip_start, + bits1, + bits2, + threshold: thresh, + cap, + mut total, + skip_reservation: skip_rsv, + intensity_reservation: mut intensity_rsv, + stereo_reservation: mut dual_stereo_rsv, + channels: c, + lm, + decoder: dec, + } = input; let ci = c as i32; let alloc_floor = ci << BITRES; let stereo = (c > 1) as i32; @@ -409,7 +445,7 @@ fn interp_bits2pulses( let mut psum = 0; let mut done = false; for j in (start..end).rev() { - let tmp = bits1[j] + (mid * bits2[j] >> ALLOC_STEPS); + let tmp = bits1[j] + ((mid * bits2[j]) >> ALLOC_STEPS); if tmp >= thresh[j] || done { done = true; psum += tmp.min(cap[j]); @@ -427,7 +463,7 @@ fn interp_bits2pulses( let mut psum = 0; let mut done = false; for j in (start..end).rev() { - let mut tmp = bits1[j] + (lo * bits2[j] >> ALLOC_STEPS); + let mut tmp = bits1[j] + ((lo * bits2[j]) >> ALLOC_STEPS); if tmp < thresh[j] && !done { tmp = if tmp >= alloc_floor { alloc_floor } else { 0 }; } else { @@ -507,9 +543,9 @@ fn interp_bits2pulses( if n == 2 { offset += (den << BITRES) >> 2; } - if bits[j] + offset < den * 2 << BITRES { + if bits[j] + offset < (den * 2) << BITRES { offset += nclogn >> 2; - } else if bits[j] + offset < den * 3 << BITRES { + } else if bits[j] + offset < (den * 3) << BITRES { offset += nclogn >> 3; } let mut eb = (bits[j] + offset + (den << (BITRES - 1))).max(0); @@ -520,7 +556,7 @@ fn interp_bits2pulses( eb = eb.min(MAX_FINE_BITS); ebits[j] = eb; fine_priority[j] = (eb * (den << BITRES) >= bits[j] + offset) as i32; - bits[j] -= ci * eb << BITRES; + bits[j] -= (ci * eb) << BITRES; } else { excess = (bit - (ci << BITRES)).max(0); bits[j] = bit - excess; @@ -530,7 +566,7 @@ fn interp_bits2pulses( if excess > 0 { let extra_fine = (excess >> (stereo + BITRES as i32)).min(MAX_FINE_BITS - ebits[j]); ebits[j] += extra_fine; - let extra_bits = extra_fine * ci << BITRES; + let extra_bits = (extra_fine * ci) << BITRES; fine_priority[j] = (extra_bits >= excess - balance) as i32; excess -= extra_bits; } diff --git a/crates/audio/src/opus/decode/decoder.rs b/crates/audio/src/opus/decode/decoder.rs index 0307b62b..46678f39 100644 --- a/crates/audio/src/opus/decode/decoder.rs +++ b/crates/audio/src/opus/decode/decoder.rs @@ -29,6 +29,34 @@ pub struct StreamDecoder { gain_q8: i16, } +struct SilkFrame { + mode: Mode, + bandwidth: Option, + has_packet: bool, + fec: bool, + frame_size: usize, + audio_size: usize, + rate: usize, + channels: usize, +} + +struct Redundancy { + present: bool, + celt_to_silk: bool, + bytes: usize, + packet_len: usize, +} + +fn apply_output_gain(pcm: &mut [f32], samples: usize, channels: usize, gain_q8: i16) { + if gain_q8 == 0 { + return; + } + let gain = (6.488_141e-4f32 * gain_q8 as f32).exp2(); + for sample in pcm.iter_mut().take(samples * channels) { + *sample *= gain; + } +} + impl StreamDecoder { /// Creates a decoder producing `channels` (1 or 2) at `rate` (8/12/16/24/48 kHz). pub fn new(rate: u32, channels: usize) -> Result { @@ -111,12 +139,20 @@ impl StreamDecoder { } } - /// Conceals `n` samples (a multiple of 2.5 ms). + /// Longest concealment or FEC span accepted per call: 120 ms, the longest Opus packet. + fn max_span(&self) -> usize { + self.rate as usize * 3 / 25 + } + + /// Conceals `n` samples (a multiple of 2.5 ms, at most 120 ms). pub fn conceal(&mut self, n: usize, out: &mut Vec) -> Result { let f2_5 = self.rate as usize / 400; if n == 0 || !n.is_multiple_of(f2_5) { return Err(Error::InvalidArgument("concealment length must be a multiple of 2.5 ms")); } + if n > self.max_span() { + return Err(Error::InvalidArgument("concealment length exceeds 120 ms")); + } let start = out.len(); out.resize(start + n * self.channels, 0.0); let mut done = 0; @@ -129,7 +165,11 @@ impl StreamDecoder { } /// Decodes the in-band FEC (LBRR) of `packet` for the previous (lost) packet of `n` samples. + /// `n` is at most 120 ms. pub fn decode_fec(&mut self, packet: &[u8], n: usize, out: &mut Vec) -> Result { + if n > self.max_span() { + return Err(Error::InvalidArgument("FEC length exceeds 120 ms")); + } let p = Packet::parse(packet)?; let pfs = p.toc.frame_samples(self.rate); if n < pfs || p.toc.mode == Mode::CeltOnly || self.mode == Some(Mode::CeltOnly) { @@ -170,6 +210,87 @@ impl StreamDecoder { Ok(nb) } + fn decode_silk_frame(&mut self, frame: SilkFrame, dec: &mut Option<&mut RangeDecoder>) -> Result> { + let SilkFrame { mode, bandwidth, has_packet, fec, frame_size, audio_size, rate, channels } = frame; + let mut pcm = Vec::new(); + if mode == Mode::CeltOnly { + return Ok(pcm); + } + pcm.resize(frame_size.max(rate / 100) * channels, 0.0); + if self.prev_mode == Some(Mode::CeltOnly) { + self.silk.reset(); + } + let internal_rate = if mode == Mode::SilkOnly { + match bandwidth { + Some(Bandwidth::Narrowband) => 8000, + Some(Bandwidth::Mediumband) => 12000, + _ => 16000, + } + } else { + 16000 + }; + if has_packet { + self.silk.set_stream(self.stream_channels, internal_rate); + } + let lost = if !has_packet { + LostFlag::Lost + } else if fec { + LostFlag::Lbrr + } else { + LostFlag::Normal + }; + let payload_ms = (1000 * audio_size / rate).max(10); + let mut decoded = 0; + while decoded < frame_size { + let first = decoded == 0; + let count = match self.silk.decode(dec.as_deref_mut(), lost, first, payload_ms, &mut pcm[decoded * channels..]) { + Ok(count) => count, + Err(_error) if lost != LostFlag::Normal => { + let remaining = frame_size - decoded; + pcm[decoded * channels..(decoded + remaining) * channels].fill(0.0); + remaining + } + Err(error) => return Err(error), + }; + if count == 0 { + break; + } + decoded += count; + } + Ok(pcm) + } + + fn read_redundancy(&self, mode: Mode, fec: bool, packet_len: usize, dec: &mut Option<&mut RangeDecoder>) -> Redundancy { + let mut state = Redundancy { present: false, celt_to_silk: false, bytes: 0, packet_len }; + let Some(decoder) = dec.as_deref_mut() else { return state }; + let hybrid = self.mode == Some(Mode::Hybrid); + if fec || mode == Mode::CeltOnly || decoder.tell() + 17 + 20 * i32::from(hybrid) > 8 * state.packet_len as i32 { + return state; + } + state.present = if mode == Mode::Hybrid { decoder.bit_logp(12) } else { true }; + if !state.present { + return state; + } + state.celt_to_silk = decoder.bit_logp(1); + state.bytes = if mode == Mode::Hybrid { decoder.uint(256) as usize + 2 } else { state.packet_len - ((decoder.tell() as usize + 7) >> 3) }; + if state.bytes > state.packet_len || (state.packet_len - state.bytes) * 8 < decoder.tell() as usize { + state.packet_len = 0; + state.bytes = 0; + state.present = false; + } else { + state.packet_len -= state.bytes; + } + decoder.shrink_storage(state.bytes); + state + } + + fn is_mode_transition(&self, has_packet: bool, mode: Mode) -> bool { + has_packet + && self.prev_mode.is_some() + && ((mode == Mode::CeltOnly && self.prev_mode != Some(Mode::CeltOnly) && !self.prev_redundancy) + || (mode != Mode::CeltOnly && self.prev_mode == Some(Mode::CeltOnly))) + } + /// `opus_decode_frame`: decodes one frame (or conceals when `data` is `None`/≤1 byte). #[allow(clippy::cognitive_complexity, reason = "the RFC 6716 frame decoding sequence, ported as-is from FilmCraft")] fn decode_frame(&mut self, data: Option<&[u8]>, pcm: &mut [f32], frame_size: usize, fec: bool) -> Result { @@ -229,14 +350,7 @@ impl StreamDecoder { } None => None, }; - let mut transition = false; - if data.is_some() - && self.prev_mode.is_some() - && ((mode == Mode::CeltOnly && self.prev_mode != Some(Mode::CeltOnly) && !self.prev_redundancy) - || (mode != Mode::CeltOnly && self.prev_mode == Some(Mode::CeltOnly))) - { - transition = true; - } + let mut transition = self.is_mode_transition(data.is_some(), mode); let mut pcm_transition = Vec::new(); if transition && mode == Mode::CeltOnly { pcm_transition.resize(f5 * ch, 0.0); @@ -249,78 +363,14 @@ impl StreamDecoder { frame_size = audiosize; // SILK. - let mut pcm_silk: Vec = Vec::new(); - if mode != Mode::CeltOnly { - pcm_silk.resize(frame_size.max(f10) * ch, 0.0); - if self.prev_mode == Some(Mode::CeltOnly) { - self.silk.reset(); - } - let payload_ms = (1000 * audiosize / fs).max(10); - let internal_rate = if mode == Mode::SilkOnly { - match bandwidth { - Some(Bandwidth::Narrowband) => 8000, - Some(Bandwidth::Mediumband) => 12000, - _ => 16000, - } - } else { - 16000 - }; - if data.is_some() { - self.silk.set_stream(self.stream_channels, internal_rate); - } - let lost = if data.is_none() { - LostFlag::Lost - } else if fec { - LostFlag::Lbrr - } else { - LostFlag::Normal - }; - let mut decoded = 0; - while decoded < frame_size { - let first = decoded == 0; - let n = match self.silk.decode(dec.as_deref_mut(), lost, first, payload_ms, &mut pcm_silk[decoded * ch..]) { - Ok(n) => n, - Err(e) => { - if lost != LostFlag::Normal { - let n = frame_size - decoded; - pcm_silk[decoded * ch..(decoded + n) * ch].fill(0.0); - n - } else { - return Err(e); - } - } - }; - if n == 0 { - break; - } - decoded += n; - } - } + let pcm_silk = self.decode_silk_frame( + SilkFrame { mode, bandwidth, has_packet: data.is_some(), fec, frame_size, audio_size: audiosize, rate: fs, channels: ch }, + &mut dec, + )?; + let Redundancy { present: redundancy, celt_to_silk, bytes: redundancy_bytes, packet_len: len } = + self.read_redundancy(mode, fec, data_len, &mut dec); let mut start_band = 0; - let mut redundancy = false; - let mut celt_to_silk = false; - let mut redundancy_bytes = 0usize; - let mut len = data_len; - if !fec - && mode != Mode::CeltOnly - && let Some(d) = dec.as_deref_mut() - && d.tell() + 17 + 20 * (self.mode == Some(Mode::Hybrid)) as i32 <= 8 * len as i32 - { - redundancy = if mode == Mode::Hybrid { d.bit_logp(12) } else { true }; - if redundancy { - celt_to_silk = d.bit_logp(1); - redundancy_bytes = if mode == Mode::Hybrid { d.uint(256) as usize + 2 } else { len - ((d.tell() as usize + 7) >> 3) }; - if redundancy_bytes > len || (len - redundancy_bytes) * 8 < d.tell() as usize { - len = 0; - redundancy_bytes = 0; - redundancy = false; - } else { - len -= redundancy_bytes; - } - d.shrink_storage(redundancy_bytes); - } - } if mode != Mode::CeltOnly { start_band = 17; } @@ -392,12 +442,7 @@ impl StreamDecoder { smooth_fade(&pcm_transition, &cur, pcm, f2_5, ch, window, self.rate); } } - if self.gain_q8 != 0 { - let g = (6.488_141e-4f32 * self.gain_q8 as f32).exp2(); - for v in pcm[..frame_size * ch].iter_mut() { - *v *= g; - } - } + apply_output_gain(pcm, frame_size, ch, self.gain_q8); self.range_final = match dec.as_deref() { Some(d) if data_len > 1 => d.rng() ^ redundant_rng, _ => 0, diff --git a/crates/audio/src/opus/decode/silk/mod.rs b/crates/audio/src/opus/decode/silk/mod.rs index 62886fb9..338b700a 100644 --- a/crates/audio/src/opus/decode/silk/mod.rs +++ b/crates/audio/src/opus/decode/silk/mod.rs @@ -694,8 +694,94 @@ impl SilkDecoder { self.internal_rate = internal_rate; } + fn prepare_frame(&mut self, new_packet: bool, payload_ms: usize) -> Result { + let nci = self.channels_internal; + if new_packet { + for channel in &mut self.ch { + channel.n_frames_decoded = 0; + } + } + if nci > self.channels_internal_prev && self.channels_internal_prev != 0 { + self.ch[1] = Channel::new(); + } + let stereo_to_mono = nci == 1 && self.channels_internal_prev == 2 && self.internal_rate == 1000 * self.ch[0].fs_khz as u32; + if self.ch[0].n_frames_decoded == 0 { + let (frames, subframes) = match payload_ms { + 10 => (1, 2), + 20 => (1, 4), + 40 => (2, 4), + 60 => (3, 4), + _ => return Err(Error::InvalidPacket("bad SILK frame size")), + }; + let fs_khz = ((self.internal_rate >> 10) + 1) as usize; + for channel in self.ch.iter_mut().take(nci) { + channel.n_frames_per_packet = frames; + channel.nb_subfr = subframes; + channel.set_fs(fs_khz, self.api_rate); + } + } + if self.api_channels == 2 && nci == 2 && self.channels_internal_prev == 1 { + self.stereo.pred_prev_q13 = [0; 2]; + self.stereo.s_side = [0; 2]; + self.ch[1].resampler = self.ch[0].resampler.clone(); + } + self.channels_internal_prev = nci; + Ok(stereo_to_mono) + } + + fn decode_packet_header(&mut self, dec: &mut Option<&mut RangeDecoder>, lost: LostFlag) -> Result<()> { + let nci = self.channels_internal; + if lost == LostFlag::Lost || self.ch[0].n_frames_decoded != 0 { + return Ok(()); + } + let decoder = dec.as_deref_mut().ok_or(Error::InvalidPacket("missing SILK data"))?; + for channel in self.ch.iter_mut().take(nci) { + for index in 0..channel.n_frames_per_packet { + channel.vad_flags[index] = decoder.bit_logp(1); + } + channel.lbrr_flag = decoder.bit_logp(1); + } + for channel in self.ch.iter_mut().take(nci) { + channel.lbrr_flags = [false; 3]; + if channel.lbrr_flag { + if channel.n_frames_per_packet == 1 { + channel.lbrr_flags[0] = true; + } else { + let sym = decoder.icdf(LBRR_FLAGS[channel.n_frames_per_packet - 2], 8) + 1; + for index in 0..channel.n_frames_per_packet { + channel.lbrr_flags[index] = (sym >> index) & 1 != 0; + } + } + } + } + if lost == LostFlag::Normal { + self.skip_lbrr_data(decoder, nci); + } + Ok(()) + } + + fn skip_lbrr_data(&mut self, decoder: &mut RangeDecoder, nci: usize) { + for frame in 0..self.ch[0].n_frames_per_packet { + for channel in 0..nci { + if !self.ch[channel].lbrr_flags[frame] { + continue; + } + if nci == 2 && channel == 0 { + stereo_decode_pred(decoder); + if !self.ch[1].lbrr_flags[frame] { + decoder.icdf(&MID_ONLY, 8); + } + } + let condition = if frame > 0 && self.ch[channel].lbrr_flags[frame - 1] { Cond::Conditionally } else { Cond::Independently }; + let state = &mut self.ch[channel]; + state.decode_indices(decoder, frame, true, condition); + let mut pulses = [0i32; MAX_FRAME + 16]; + decode_pulses(decoder, &mut pulses, state.indices.signal_type, state.indices.quant_offset_type, state.frame_length); + } + } + } + /// Decodes one SILK frame (10 or 20 ms) into interleaved `out` at the API rate. - #[allow(clippy::cognitive_complexity, reason = "the RFC 6716 SILK frame decode, ported as-is from FilmCraft")] pub fn decode( &mut self, mut dec: Option<&mut RangeDecoder>, @@ -705,79 +791,8 @@ impl SilkDecoder { out: &mut [f32], ) -> Result { let nci = self.channels_internal; - if new_packet { - for c in self.ch.iter_mut() { - c.n_frames_decoded = 0; - } - } - if nci > self.channels_internal_prev && self.channels_internal_prev != 0 { - self.ch[1] = Channel::new(); - } - let stereo_to_mono = nci == 1 && self.channels_internal_prev == 2 && self.internal_rate == 1000 * self.ch[0].fs_khz as u32; - if self.ch[0].n_frames_decoded == 0 { - for c in self.ch.iter_mut().take(nci) { - let (nfpp, nb) = match payload_ms { - 10 => (1, 2), - 20 => (1, 4), - 40 => (2, 4), - 60 => (3, 4), - _ => return Err(Error::InvalidPacket("bad SILK frame size")), - }; - c.n_frames_per_packet = nfpp; - c.nb_subfr = nb; - let fs_khz = ((self.internal_rate >> 10) + 1) as usize; - c.set_fs(fs_khz, self.api_rate); - } - } - if self.api_channels == 2 && nci == 2 && self.channels_internal_prev == 1 { - self.stereo.pred_prev_q13 = [0; 2]; - self.stereo.s_side = [0; 2]; - self.ch[1].resampler = self.ch[0].resampler.clone(); - } - self.channels_internal_prev = nci; - - if lost != LostFlag::Lost && self.ch[0].n_frames_decoded == 0 { - let d = dec.as_deref_mut().ok_or(Error::InvalidPacket("missing SILK data"))?; - for c in self.ch.iter_mut().take(nci) { - for i in 0..c.n_frames_per_packet { - c.vad_flags[i] = d.bit_logp(1); - } - c.lbrr_flag = d.bit_logp(1); - } - for c in self.ch.iter_mut().take(nci) { - c.lbrr_flags = [false; 3]; - if c.lbrr_flag { - if c.n_frames_per_packet == 1 { - c.lbrr_flags[0] = true; - } else { - let sym = d.icdf(LBRR_FLAGS[c.n_frames_per_packet - 2], 8) + 1; - for i in 0..c.n_frames_per_packet { - c.lbrr_flags[i] = (sym >> i) & 1 != 0; - } - } - } - } - if lost == LostFlag::Normal { - // Skip LBRR data. - for i in 0..self.ch[0].n_frames_per_packet { - for n in 0..nci { - if self.ch[n].lbrr_flags[i] { - if nci == 2 && n == 0 { - stereo_decode_pred(d); - if !self.ch[1].lbrr_flags[i] { - d.icdf(&MID_ONLY, 8); - } - } - let cond = if i > 0 && self.ch[n].lbrr_flags[i - 1] { Cond::Conditionally } else { Cond::Independently }; - let c = &mut self.ch[n]; - c.decode_indices(d, i, true, cond); - let mut pulses = [0i32; MAX_FRAME + 16]; - decode_pulses(d, &mut pulses, c.indices.signal_type, c.indices.quant_offset_type, c.frame_length); - } - } - } - } - } + let stereo_to_mono = self.prepare_frame(new_packet, payload_ms)?; + self.decode_packet_header(&mut dec, lost)?; // Stereo prediction. Each frame starts out coding the side channel, as in libopus. let mut decode_only_middle = false; if nci == 2 { diff --git a/crates/audio/src/opus/decode/silk/resampler.rs b/crates/audio/src/opus/decode/silk/resampler.rs index b8ae6f5b..0671778c 100644 --- a/crates/audio/src/opus/decode/silk/resampler.rs +++ b/crates/audio/src/opus/decode/silk/resampler.rs @@ -4,6 +4,26 @@ //! RFC 6716 Table 54 (and the integer delays used when no rate conversion is needed), so that the //! SILK output lines up with the CELT layer in Hybrid mode. +/// Dot product with eight independent partial sums, so the additions pipeline and the loop +/// vectorises instead of waiting on one 4-cycle add chain. Rounds differently from a +/// sequential sum (the resampler is not normative). +#[inline] +fn dot8(a: &[f32], b: &[f32]) -> f32 { + let mut s = [0f32; 8]; + let (ca, ra) = a.as_chunks::<8>(); + let (cb, rb) = b.as_chunks::<8>(); + for (x, y) in ca.iter().zip(cb) { + for l in 0..8 { + s[l] += x[l] * y[l]; + } + } + let mut acc = ((s[0] + s[4]) + (s[1] + s[5])) + ((s[2] + s[6]) + (s[3] + s[7])); + for (x, y) in ra.iter().zip(rb) { + acc += x * y; + } + acc +} + #[derive(Clone)] pub struct Resampler { fin: u32, @@ -142,17 +162,32 @@ impl Resampler { self.hist.extend(input.iter().map(|&v| v as f32)); self.n_in += input.len() as i64; let total_out = (self.n_in as u64 * self.fout as u64) / self.fin as u64; + // phases[ph] repeats every `period` outputs, advancing the input position by `step` + // (= fin * period / fout, exact because period = fout / gcd). + let step = self.fin as i64 * self.period as i64 / self.fout as i64; + let mut cycle = (self.n_out / self.period) as i64; + let mut ph = (self.n_out % self.period) as usize; while self.n_out < total_out { - let j = self.n_out; - let cycle = (j / self.period) as i64; - let (base0, ref taps) = self.phases[(j % self.period) as usize]; - let base = base0 + cycle * self.fin as i64 * self.period as i64 / self.fout as i64; + let (base0, ref taps) = self.phases[ph]; + let base = base0 + cycle * step; + ph += 1; + if ph as u64 == self.period { + ph = 0; + cycle += 1; + } let first = base - self.half + 1; + let i0 = first - self.hist_start; let mut acc = 0f32; - for (m, &t) in taps.iter().enumerate() { - let idx = first + m as i64 - self.hist_start; - if idx >= 0 && (idx as usize) < self.hist.len() { - acc += self.hist[idx as usize] * t; + if i0 >= 0 && i0 as usize + taps.len() <= self.hist.len() { + // Window fully inside the history (always, after the first call). + let win = &self.hist[i0 as usize..i0 as usize + taps.len()]; + acc = dot8(win, taps); + } else { + for (m, &t) in taps.iter().enumerate() { + let idx = first + m as i64 - self.hist_start; + if idx >= 0 && (idx as usize) < self.hist.len() { + acc += self.hist[idx as usize] * t; + } } } out.push(acc); @@ -202,4 +237,14 @@ mod tests { assert!(snr > 35.0, "{fin}->{fout}: {snr} dB"); } } + + #[test] + fn dot8_matches_sequential_sum() { + for len in [1usize, 7, 8, 9, 35, 36, 71] { + let a: Vec = (0..len).map(|i| ((i * 37 % 19) as f32 - 9.0) * 1000.0).collect(); + let b: Vec = (0..len).map(|i| ((i * 11 % 7) as f32 - 3.0) / 7.0).collect(); + let seq: f64 = a.iter().zip(&b).map(|(&x, &y)| x as f64 * y as f64).sum(); + assert!((dot8(&a, &b) as f64 - seq).abs() < 1e-2, "len {len}"); + } + } } diff --git a/crates/audio/src/read/caf.rs b/crates/audio/src/read/caf.rs index 94ae7f04..d1c6f69c 100644 --- a/crates/audio/src/read/caf.rs +++ b/crates/audio/src/read/caf.rs @@ -20,7 +20,7 @@ fn chunks(b: &[u8], mut pos: usize) -> Result<()> { // Whether the packet table lists each packet's size and frame count (variable per packet). let (mut lists_sizes, mut lists_frames, mut described) = (false, false, false); while let (Some(kind), Some(size)) = (b.get(pos..pos + 4), be_u64(b, pos + 4)) { - let (body, size) = (pos + 12, size as i64); + let body = pos + 12; let used = match kind { b"desc" if size == 32 => { lists_sizes = be_u32(b, body + 16) == Some(0); @@ -28,12 +28,12 @@ fn chunks(b: &[u8], mut pos: usize) -> Result<()> { described = true; Some(32) } - b"data" if size == -1 => Some(4), + b"data" if size == u64::MAX => Some(4), b"chan" if size >= 12 => be_u32(b, body + 8).map(|n| 12 + u64::from(n) * 20), b"pakt" if size >= 24 && described => packet_table(b, body, lists_sizes, lists_frames)?, - b"kuki" if size >= 0 && size as u64 > b.len().saturating_sub(body) as u64 => return Err(malformed("magic cookie")), + b"kuki" if size > b.len().saturating_sub(body) as u64 => return Err(malformed("magic cookie")), b"desc" | b"pakt" | b"chan" => None, - _ => u64::try_from(size).ok(), + _ => Some(size), }; match used.filter(|_| described).and_then(|n| usize::try_from(n).ok()).and_then(|n| body.checked_add(n)) { Some(next) if next < b.len() => pos = next, @@ -85,3 +85,13 @@ fn varint(b: &[u8], at: usize) -> Option<(u64, usize)> { } None } + +#[cfg(test)] +mod tests { + #[test] + fn rejects_cookie_size_above_signed_range() { + let mut bytes = b"caff\0\x01\0\0kuki".to_vec(); + bytes.extend_from_slice(&(1u64 << 63).to_be_bytes()); + assert!(super::check(&bytes, 0).is_err()); + } +} diff --git a/crates/audio/src/read/flac.rs b/crates/audio/src/read/flac.rs index a16e7bdf..e0828dc4 100644 --- a/crates/audio/src/read/flac.rs +++ b/crates/audio/src/read/flac.rs @@ -154,10 +154,13 @@ fn comments_mask(block: &[u8]) -> Option { None } -/// Whether a comment in a (little-endian) Vorbis comment block claims more than the block holds. -/// The vendor string is skipped unread. +/// Whether a vendor string or comment in a (little-endian) Vorbis comment block claims more than +/// the block holds. fn comments_overrun(b: &[u8], body: usize, end: usize) -> bool { - let Some(vendor) = le_u32(b, body) else { return false }; + let Some(vendor) = le_u32(b, body).filter(|_| body.saturating_add(4) <= end) else { return false }; + if vendor as usize > end - (body + 4) { + return true; + } let at = (body + 4).saturating_add(vendor as usize); let Some(count) = le_u32(b, at).filter(|_| at.saturating_add(4) <= end) else { return false }; overruns(b, at + 4, end, le_u32, std::iter::repeat_n(0, count as usize)) diff --git a/crates/audio/src/write/flac/encoder.rs b/crates/audio/src/write/flac/encoder.rs index 2e01c62d..02a22a0c 100644 --- a/crates/audio/src/write/flac/encoder.rs +++ b/crates/audio/src/write/flac/encoder.rs @@ -56,6 +56,8 @@ pub struct Encoder { min_frame: u32, max_frame: u32, md5: super::md5::Md5, + /// The LPC analysis window for the current block length. + window: Vec, } impl Encoder { @@ -84,6 +86,7 @@ impl Encoder { // STREAMINFO holds an MD5 signature of the unencoded audio (RFC 9639 §8.2); this is an // integrity checksum, not a password hash or a security boundary. md5: super::md5::Md5::default(), // NOSONAR + window: Vec::new(), }) } @@ -157,9 +160,16 @@ impl Encoder { } let bs = usize::from(self.cfg.block_size); let mut out = Vec::new(); - while self.pending.first().is_some_and(|p| p.len() >= bs) { - let block: Vec> = self.pending.iter_mut().map(|p| p.drain(..bs).collect()).collect(); + // Walk the buffered samples by offset and compact once: draining each block from the + // front moved the whole remainder per frame, quadratic in the batch length. + let mut at = 0; + while self.pending.first().is_some_and(|p| p.len() - at >= bs) { + let block: Vec> = self.pending.iter().map(|p| p.get(at..at + bs).unwrap_or_default().to_vec()).collect(); out.push(self.frame(&block)); + at += bs; + } + for p in &mut self.pending { + p.drain(..at); } Ok(out) } @@ -187,7 +197,11 @@ impl Encoder { self.md5.update(&raw); let effort = self.effort; - let code = |x: &[i32], b: u8| super::predict::choose(x, b, effort); + if self.window.len() != n { + self.window = super::predict::window(n); + } + let win = &self.window; + let code = |x: &[i32], b: u8| super::predict::choose(x, b, effort, win); let data = match block { // the side channel takes one bit more than the samples: at most 32 (i32) [l, r] if effort.stereo && bps < 32 => { diff --git a/crates/audio/src/write/flac/frame.rs b/crates/audio/src/write/flac/frame.rs index 34756ca0..06a1bf5f 100644 --- a/crates/audio/src/write/flac/frame.rs +++ b/crates/audio/src/write/flac/frame.rs @@ -95,7 +95,8 @@ pub struct FrameSpec<'a> { pub channels: &'a [(&'a [i32], Subframe, u8)], } -fn write_residual(w: &mut BitWriter, residual: &Residual, r: &[i64], block: usize, order: usize) { +fn write_residual(w: &mut BitWriter, residual: &Residual, block: usize, order: usize) { + let r = &residual.folded; let wide = residual.wide(); w.write_bits(u32::from(wide), 2); w.write_bits(u32::from(residual.partition_order), 4); @@ -105,8 +106,8 @@ fn write_residual(w: &mut BitWriter, residual: &Residual, r: &[i64], block: usiz for (p, &k) in residual.params.iter().enumerate() { w.write_bits(u32::from(k), if wide { 5 } else { 4 }); let n = if p == 0 { len.saturating_sub(order) } else { len }; - for &v in r.get(start..start + n).unwrap_or_default() { - let u = super::predict::fold(v); + for &u in r.get(start..start + n).unwrap_or_default() { + let u = u64::from(u); let mut q = u >> k; while q >= 32 { w.write_bits(0, 32); @@ -145,8 +146,7 @@ fn write_subframe(w: &mut BitWriter, x: &[i32], sub: &Subframe, bps: u8, wasted: header(w, 0b001000 | u32::from(*order)); let o = usize::from(*order); x.iter().take(o).for_each(|&v| sample(w, v)); - let r = super::predict::fixed_residual(x, o).unwrap_or_default(); - write_residual(w, residual, &r, x.len(), o); + write_residual(w, residual, x.len(), o); } Subframe::Lpc { coefs, precision, shift, residual } => { let o = coefs.len(); @@ -155,8 +155,7 @@ fn write_subframe(w: &mut BitWriter, x: &[i32], sub: &Subframe, bps: u8, wasted: w.write_bits(u32::from(*precision) - 1, 4); w.write_bits(u32::from(*shift), 5); coefs.iter().for_each(|&c| w.write_bits(c as u32, u32::from(*precision))); - let r = super::predict::lpc_residual(x, coefs, *shift).unwrap_or_default(); - write_residual(w, residual, &r, x.len(), o); + write_residual(w, residual, x.len(), o); } Subframe::Wasted { bits, coded } => { let shifted: Vec = x.iter().map(|&v| v >> bits).collect(); diff --git a/crates/audio/src/write/flac/golden_tests.rs b/crates/audio/src/write/flac/golden_tests.rs new file mode 100644 index 00000000..ec5ef562 --- /dev/null +++ b/crates/audio/src/write/flac/golden_tests.rs @@ -0,0 +1,249 @@ +//! Golden hashes of the encoder output for a fixed set of synthetic signals. +//! +//! Encoder speed-ups must not change a single output byte, so this pins the encoded frames and the +//! final STREAMINFO (which holds the MD5 of the input) of 64 deterministic cases: noise, sines, silence, +//! full-scale squares, wasted bits and a music-like mix, at 8/12/16/24/32 bits, 1 to 4 channels, +//! several block sizes (with a short final frame), and levels 0 to 8. +//! +//! The hashes were captured from FilmCraft's encoder (68919ed8), which this one came from, before the per-order +//! predictor kernels and merged partition sums were added. Six cases differ from FilmCraft's table +//! on purpose and were captured here, identical before and after those changes: the five `wasted` +//! cases (this encoder codes wasted bits) and `click-32b` (32-bit stereo is coded as independent +//! channels). Signals use integer arithmetic only. The +//! encoder itself uses `f64` (window, autocorrelation), so a platform whose libm rounds `cos` +//! differently could in principle differ; the hashes were made on x86-64 Linux. +//! +//! An intentional change of the output (for example a different LPC order choice) regenerates the +//! table with `cargo test -p craft-audio --features flac-encode golden -- --ignored --nocapture print`. + +use super::encoder::{Encoder, EncoderConfig}; + +struct Rng(u64); +impl Rng { + fn next(&mut self) -> u64 { + self.0 ^= self.0 << 13; + self.0 ^= self.0 >> 7; + self.0 ^= self.0 << 17; + self.0 + } +} + +/// FNV-1a, 64 bit. +fn fnv(bytes: &[u8]) -> u64 { + bytes.iter().fold(0xcbf2_9ce4_8422_2325, |h, &b| (h ^ u64::from(b)).wrapping_mul(0x0000_0100_0000_01b3)) +} + +/// A sine-like wave from an integer resonator: `y[n] = 2cos(w) y[n-1] - y[n-2]` in Q14, amplitude `amp`. +fn resonator(n: usize, coef_q14: i64, amp: i64) -> Vec { + let (mut y1, mut y2) = (amp, amp * 9 / 10); + (0..n) + .map(|_| { + let y = ((coef_q14 * y1) >> 14) - y2; + y2 = y1; + y1 = y.clamp(-amp, amp); + y1 + }) + .collect() +} + +/// One channel of `kind` at `bps` bits, values always inside the sample range. +fn signal(kind: &str, n: usize, bps: u8, seed: u64) -> Vec { + let max = (1i64 << (bps - 1)) - 1; + let mut r = Rng(seed | 1); + let sine = resonator(n, 32_000 - (seed as i64 % 7) * 900, max / 2); + let sine2 = resonator(n, 27_000, max / 4); + (0..n) + .map(|i| { + let noise = |r: &mut Rng, amp: i64| (r.next() >> 11) as i64 % (2 * amp + 1) - amp; + let v = match kind { + "noise" => noise(&mut r, max), + "sine" => sine[i], + "music" => sine[i] + sine2[i] + noise(&mut r, max / 64), + "square" => { + if (i / 37) % 2 == 0 { + max + } else { + -max - 1 + } + } + // low 3 bits always zero (wasted bits) on a quiet sine + "wasted" => (sine[i] / 8) * 8, + "click" => { + if i == n / 2 { + max + } else { + 0 + } + } + // a silent stretch between two noisy ones + "gap" => { + if (n / 3..2 * n / 3).contains(&i) { + 0 + } else { + noise(&mut r, max / 3) + } + } + _ => 0, + }; + v.clamp(-max - 1, max) as i32 + }) + .collect() +} + +struct Case { + kind: &'static str, + bps: u8, + channels: u8, + block: u16, + level: u8, + samples: usize, + chunk: usize, +} + +impl Case { + fn name(&self) -> String { + format!("{}-{}b-{}ch-bs{}-L{}-n{}-c{}", self.kind, self.bps, self.channels, self.block, self.level, self.samples, self.chunk) + } + + /// Hash of every frame and the final STREAMINFO (the file header carries the crate version). + fn hash(&self) -> u64 { + let planar: Vec> = + (0..self.channels).map(|c| signal(self.kind, self.samples, self.bps, 0x9E37_79B9 + u64::from(c) * 7919)).collect(); + let mut cfg = EncoderConfig::new(44_100, self.channels, self.bps); + cfg.block_size = self.block; + cfg.level = self.level; + let mut enc = Encoder::new(cfg).expect("valid config"); + let mut out = Vec::new(); + let mut at = 0; + while at < self.samples { + let end = (at + self.chunk).min(self.samples); + let parts: Vec<&[i32]> = planar.iter().map(|c| &c[at..end]).collect(); + out.extend(enc.encode(&parts).expect("samples in range")); + at = end; + } + out.extend(enc.finish()); + out.extend_from_slice(&enc.streaminfo()); + fnv(&out) + } +} + +fn cases() -> Vec { + const KINDS: [&str; 8] = ["noise", "sine", "music", "square", "wasted", "click", "gap", "silence"]; + const BPS: [u8; 5] = [8, 12, 16, 24, 32]; + const BLOCKS: [u16; 5] = [16, 576, 1000, 4096, 4608]; + let mut v = Vec::new(); + // every kind at every bit depth; channels, block size and level rotate through their ranges + for (k, kind) in KINDS.iter().enumerate() { + for (b, &bps) in BPS.iter().enumerate() { + let block = BLOCKS[(k + 2 * b) % BLOCKS.len()]; + let samples = 2 * usize::from(block) + 37 + (k * 101 + b * 13) % 700; + v.push(Case { + kind, + bps, + channels: 1 + ((k + b) % 4) as u8, + block, + level: ((k * 5 + b * 3) % 9) as u8, + samples, + chunk: 1 + (k * 977 + b * 311) % 5000, + }); + } + } + // all levels on stereo music at 16 and 24 bits, and the 4096 blocks real files use + for level in 0..=8 { + for bps in [16u8, 24] { + v.push(Case { kind: "music", bps, channels: 2, block: 4096, level, samples: 9000, chunk: 4096 }); + } + } + // 4 channels, short blocks and a block that is not a power of two + for level in [0u8, 5, 8] { + v.push(Case { kind: "noise", bps: 16, channels: 4, block: 1000, level, samples: 2500, chunk: 777 }); + v.push(Case { kind: "gap", bps: 16, channels: 2, block: 576, level, samples: 2000, chunk: 2000 }); + } + v +} + +/// (case name, hash), captured on FilmCraft 68919ed8. +const GOLDEN: &[(&str, u64)] = &[ + ("noise-8b-1ch-bs16-L0-n69-c1", 0x69cba522cd83247a), + ("noise-12b-2ch-bs1000-L3-n2050-c312", 0xbd8aa93b40b5b078), + ("noise-16b-3ch-bs4608-L6-n9279-c623", 0xe1969205fcfb7d74), + ("noise-24b-4ch-bs576-L0-n1228-c934", 0xd444404f16d70169), + ("noise-32b-1ch-bs4096-L3-n8281-c1245", 0x45ca92f380216922), + ("sine-8b-2ch-bs576-L5-n1290-c978", 0x135ff72f3b74e284), + ("sine-12b-3ch-bs4096-L8-n8343-c1289", 0x95c4d0720e59d826), + ("sine-16b-4ch-bs16-L2-n196-c1600", 0x84bf0b50ef111b77), + ("sine-24b-1ch-bs1000-L5-n2177-c1911", 0xf4cb2b792c100e1c), + ("sine-32b-2ch-bs4608-L8-n9406-c2222", 0x45aaf5d6091ecc85), + ("music-8b-3ch-bs1000-L1-n2239-c1955", 0x3540af737b613e03), + ("music-12b-4ch-bs4608-L4-n9468-c2266", 0x07d8c76385536b70), + ("music-16b-1ch-bs576-L7-n1417-c2577", 0x1b5300d81d3a8eac), + ("music-24b-2ch-bs4096-L1-n8470-c2888", 0x570f7821161e7932), + ("music-32b-3ch-bs16-L4-n323-c3199", 0xb720fd6b9fb8d7c9), + ("square-8b-4ch-bs4096-L6-n8532-c2932", 0xf1e702077432c346), + ("square-12b-1ch-bs16-L0-n385-c3243", 0x2820a440a7465bb2), + ("square-16b-2ch-bs1000-L3-n2366-c3554", 0x827ee17d45c447bc), + ("square-24b-3ch-bs4608-L6-n9595-c3865", 0x88b6af280896f131), + ("square-32b-4ch-bs576-L0-n1544-c4176", 0x8a72918d0834a2f1), + ("wasted-8b-1ch-bs4608-L2-n9657-c3909", 0xa7516cdcf761110e), + ("wasted-12b-2ch-bs576-L5-n1606-c4220", 0xf1894eef730d7546), + ("wasted-16b-3ch-bs4096-L8-n8659-c4531", 0x39d34c6ea3569b7c), + ("wasted-24b-4ch-bs16-L2-n512-c4842", 0xfa3786f5b0f1004d), + ("wasted-32b-1ch-bs1000-L5-n2493-c153", 0x00934303beb10016), + ("click-8b-2ch-bs16-L7-n574-c4886", 0x54b956cd16590b80), + ("click-12b-3ch-bs1000-L1-n2555-c197", 0x5b5c00ba683e16fb), + ("click-16b-4ch-bs4608-L4-n9784-c508", 0x906ae4960be0eb16), + ("click-24b-1ch-bs576-L7-n1733-c819", 0x815a9e001b296ae1), + ("click-32b-2ch-bs4096-L1-n8786-c1130", 0x9d4424021be06705), + ("gap-8b-3ch-bs576-L3-n1795-c863", 0xa372ac31d477b4de), + ("gap-12b-4ch-bs4096-L6-n8848-c1174", 0x0f98e2cbcff69b73), + ("gap-16b-1ch-bs16-L0-n701-c1485", 0xf0de17bb815bc8eb), + ("gap-24b-2ch-bs1000-L3-n2682-c1796", 0x09a7b31c70740cc1), + ("gap-32b-3ch-bs4608-L6-n9911-c2107", 0x8ba6d94ee8c85b0c), + ("silence-8b-4ch-bs1000-L8-n2044-c1840", 0x845724c286fa7300), + ("silence-12b-1ch-bs4608-L2-n9273-c2151", 0xe09802e7a78bf137), + ("silence-16b-2ch-bs576-L5-n1222-c2462", 0x9bfd29a65f5d4a8c), + ("silence-24b-3ch-bs4096-L8-n8275-c2773", 0xfcf34cb0553f5b39), + ("silence-32b-4ch-bs16-L2-n128-c3084", 0xc4c8f1f5d73208e0), + ("music-16b-2ch-bs4096-L0-n9000-c4096", 0x7ee21e0646af3021), + ("music-24b-2ch-bs4096-L0-n9000-c4096", 0x384bb9626dd0eeb3), + ("music-16b-2ch-bs4096-L1-n9000-c4096", 0x7ee21e0646af3021), + ("music-24b-2ch-bs4096-L1-n9000-c4096", 0x384bb9626dd0eeb3), + ("music-16b-2ch-bs4096-L2-n9000-c4096", 0x7ee21e0646af3021), + ("music-24b-2ch-bs4096-L2-n9000-c4096", 0x384bb9626dd0eeb3), + ("music-16b-2ch-bs4096-L3-n9000-c4096", 0xf948f408d6ea5f90), + ("music-24b-2ch-bs4096-L3-n9000-c4096", 0xc448a8514df4434d), + ("music-16b-2ch-bs4096-L4-n9000-c4096", 0xd6524def929fdb8e), + ("music-24b-2ch-bs4096-L4-n9000-c4096", 0xf8f9ddf3036880be), + ("music-16b-2ch-bs4096-L5-n9000-c4096", 0xd6524def929fdb8e), + ("music-24b-2ch-bs4096-L5-n9000-c4096", 0xf8f9ddf3036880be), + ("music-16b-2ch-bs4096-L6-n9000-c4096", 0xd6524def929fdb8e), + ("music-24b-2ch-bs4096-L6-n9000-c4096", 0xf8f9ddf3036880be), + ("music-16b-2ch-bs4096-L7-n9000-c4096", 0x5be70a6777d3303e), + ("music-24b-2ch-bs4096-L7-n9000-c4096", 0x70773b947683c189), + ("music-16b-2ch-bs4096-L8-n9000-c4096", 0x5be70a6777d3303e), + ("music-24b-2ch-bs4096-L8-n9000-c4096", 0x70773b947683c189), + ("noise-16b-4ch-bs1000-L0-n2500-c777", 0x5d0bbda30208b8e9), + ("gap-16b-2ch-bs576-L0-n2000-c2000", 0x294b47f7663f3861), + ("noise-16b-4ch-bs1000-L5-n2500-c777", 0x5d0bbda30208b8e9), + ("gap-16b-2ch-bs576-L5-n2000-c2000", 0xb8a5a6d18a43f0da), + ("noise-16b-4ch-bs1000-L8-n2500-c777", 0x5d0bbda30208b8e9), + ("gap-16b-2ch-bs576-L8-n2000-c2000", 0xb8a5a6d18a43f0da), +]; + +#[test] +fn encoder_output_matches_golden_hashes() { + let all = cases(); + assert_eq!(all.len(), GOLDEN.len(), "golden table out of date: regenerate it"); + for (case, &(name, want)) in all.iter().zip(GOLDEN) { + assert_eq!(case.name(), name, "case list changed"); + assert_eq!(case.hash(), want, "output changed for {name}"); + } +} + +#[test] +#[ignore = "prints the golden table"] +fn print() { + for c in cases() { + println!(" (\"{}\", {:#018x}),", c.name(), c.hash()); + } +} diff --git a/crates/audio/src/write/flac/mod.rs b/crates/audio/src/write/flac/mod.rs index 720b59b4..f2674d42 100644 --- a/crates/audio/src/write/flac/mod.rs +++ b/crates/audio/src/write/flac/mod.rs @@ -10,6 +10,8 @@ mod predict; #[cfg(test)] use craft_media_testkit as common; #[cfg(test)] +mod golden_tests; +#[cfg(test)] mod oracle_tests; use super::pcm::{Quantizer, validate_buffer, validate_layout}; diff --git a/crates/audio/src/write/flac/predict.rs b/crates/audio/src/write/flac/predict.rs index 7a99464e..6ab7e123 100644 --- a/crates/audio/src/write/flac/predict.rs +++ b/crates/audio/src/write/flac/predict.rs @@ -13,6 +13,8 @@ pub struct Residual { pub partition_order: u8, /// Total bits of the residual section, coding method included. pub bits: u64, + /// The folded residual after the warm-up samples. + pub folded: Vec, } impl Residual { @@ -56,16 +58,14 @@ fn rice_estimate(sum: u64, n: u64, k: u32) -> u64 { (sum >> k) + n * (u64::from(k) + 1) } -fn rice_bits_exact(u: &[u64], k: u32) -> u64 { - u.iter().map(|&v| (v >> k) + u64::from(k) + 1).sum() +fn rice_bits_exact(u: &[u32], k: u32) -> u64 { + u.iter().map(|&v| u64::from(v >> k)).sum::() + u.len() as u64 * (u64::from(k) + 1) } -/// The best parameter for a partition of folded values. -fn best_param(u: &[u64]) -> (u8, u64) { - let n = u.len() as u64; - let sum: u64 = u.iter().sum(); +/// The parameter for a partition of `n` folded values summing to `sum`. +fn param_for(sum: u64, n: u64) -> u32 { if n == 0 { - return (0, 0); + return 0; } let mean = sum / n; let guess = if mean == 0 { 0 } else { 63 - mean.leading_zeros() }; @@ -76,85 +76,147 @@ fn best_param(u: &[u64]) -> (u8, u64) { best = (k, b); } } - (best.0 as u8, rice_bits_exact(u, best.0)) + best.0 +} + +/// One partition while choosing the partition order. +#[derive(Clone, Copy)] +struct Part { + sum: u64, + n: u64, + k: u32, + bits: u64, } /// Choose the partition order and parameters for the residual of a block of `block` samples whose /// first `order` samples are warm-up (`u` holds the folded residual after them). -fn code_residual(u: &[u64], block: usize, order: usize, max_partition_order: u8) -> Residual { - let mut best: Option = None; +/// +/// The partition sums of the finest usable order are added pairwise for the coarser orders, and a +/// partition whose two halves chose the same parameter reuses their exact bit counts. +fn code_residual(u: Vec, block: usize, order: usize, max_partition_order: u8) -> Residual { + let mut top = None; for po in 0..=max_partition_order { let parts = 1usize << po; - if !block.is_multiple_of(parts) { + if !block.is_multiple_of(parts) || block / parts <= order { break; } - let len = block / parts; - if len <= order { - break; - } - let mut params = Vec::with_capacity(parts); - let mut bits = 0u64; + top = Some(po); + } + let Some(top) = top else { return Residual { params: vec![0], partition_order: 0, bits: u64::MAX, folded: u } }; + let len = block >> top; + let mut finest = Vec::with_capacity(1 << top); + let mut start = 0usize; + for p in 0..1usize << top { + let n = if p == 0 { len - order } else { len }; + let chunk = u.get(start..start + n).unwrap_or_default(); + let sum: u64 = chunk.iter().map(|&v| u64::from(v)).sum(); + let k = param_for(sum, n as u64); + finest.push(Part { sum, n: n as u64, k, bits: rice_bits_exact(chunk, k) }); + start += n; + } + let mut levels = vec![finest]; + for po in (0..top).rev() { + let len = block >> po; + let Some(finer) = levels.last() else { break }; + let mut level = Vec::with_capacity(finer.len() / 2); let mut start = 0usize; - for p in 0..parts { - let n = if p == 0 { len - order } else { len }; - let Some(chunk) = u.get(start..start + n) else { break }; - let (k, b) = best_param(chunk); - params.push(k); - bits += b; - start += n; + for (p, &[a, b]) in finer.as_chunks::<2>().0.iter().enumerate() { + // a and b: the two halves of partition p, one order finer + let (sum, n) = (a.sum + b.sum, a.n + b.n); + let k = param_for(sum, n); + let bits = if a.k == k && b.k == k { a.bits + b.bits } else { rice_bits_exact(u.get(start..start + n as usize).unwrap_or_default(), k) }; + level.push(Part { sum, n, k, bits }); + start += if p == 0 { len - order } else { len }; } - let wide = params.iter().any(|&k| k > 14); - bits += 2 + 4 + parts as u64 * if wide { 5 } else { 4 }; - if best.as_ref().is_none_or(|b| bits < b.bits) { - best = Some(Residual { params, partition_order: po, bits }); + levels.push(level); + } + levels.reverse(); + let mut best: Option<(usize, u64)> = None; + for (po, level) in levels.iter().enumerate() { + let wide = level.iter().any(|p| p.k > 14); + let bits = level.iter().map(|p| p.bits).sum::() + 2 + 4 + level.len() as u64 * if wide { 5 } else { 4 }; + if best.is_none_or(|b| bits < b.1) { + best = Some((po, bits)); } } - best.unwrap_or(Residual { params: vec![0], partition_order: 0, bits: u64::MAX }) + let (po, bits) = best.unwrap_or((0, u64::MAX)); + let params = levels.get(po).map_or_else(Vec::new, |l| l.iter().map(|p| p.k as u8).collect()); + Residual { params, partition_order: po as u8, bits, folded: u } } -/// The residual of fixed predictor `order` (0–4), or None when it does not fit. -pub fn fixed_residual(x: &[i32], order: usize) -> Option> { - let x: Vec = x.iter().map(|&v| i64::from(v)).collect(); - let mut out = Vec::with_capacity(x.len().saturating_sub(order)); - for i in order..x.len() { - let s = |j: usize| x.get(i - j).copied().unwrap_or(0); - let r = match order { - 0 => s(0), - 1 => s(0) - s(1), - 2 => s(0) - 2 * s(1) + s(2), - 3 => s(0) - 3 * s(1) + 3 * s(2) - s(3), - _ => s(0) - 4 * s(1) + 6 * s(2) - 4 * s(3) + s(4), - }; +/// The folded residual of an order-`N` integer predictor (`c[0]` weighs the previous sample), or +/// None when a residual does not fit. Monomorphised per order so the dot product is unrolled. +fn predict_fold(x: &[i32], c: &[i32], shift: u8, lpc: bool) -> Option> { + let c: [i32; N] = c.try_into().ok()?; + let mut out = Vec::with_capacity(x.len().saturating_sub(N)); + for w in x.windows(N + 1) { + let mut p = 0i64; + for j in 0..N { + p += i64::from(c[j]) * i64::from(w[N - 1 - j]); + } + let r = i64::from(w[N]) - prediction(p >> shift, lpc)?; if r.abs() > MAX_RESIDUAL { return None; } - out.push(r); + out.push(fold(r) as u32); } Some(out) } -/// The residual of quantised linear prediction, or None when it does not fit. -pub fn lpc_residual(x: &[i32], coefs: &[i32], shift: u8) -> Option> { - let order = coefs.len(); - let mut out = Vec::with_capacity(x.len().saturating_sub(order)); - for i in order..x.len() { - let hist = x.get(i - order..i)?; - // coefs[0] weighs the previous sample - let pred: i64 = coefs.iter().zip(hist.iter().rev()).map(|(&c, &s)| i64::from(c) * i64::from(s)).sum(); - // Valid past the i32 range too, but symphonia adds it to the residual as an i32, which - // panics in a debug build: its overflow checks are on in apps without our profile override. - let pred = i32::try_from(pred >> shift).ok()?; - let r = i64::from(*x.get(i)?) - i64::from(pred); +/// A prediction to subtract. A linear prediction (`lpc`) must fit `i32`: larger ones are valid +/// FLAC, but symphonia adds the prediction to the residual as an `i32`, which panics in a debug +/// build (its overflow checks are on in apps without our profile override). +fn prediction(p: i64, lpc: bool) -> Option { + if lpc { i32::try_from(p).ok().map(i64::from) } else { Some(p) } +} + +fn predict_fold_any(x: &[i32], c: &[i32], shift: u8, lpc: bool) -> Option> { + let mut out = Vec::with_capacity(x.len().saturating_sub(c.len())); + for w in x.windows(c.len() + 1) { + let p: i64 = c.iter().zip(w.iter().rev().skip(1)).map(|(&c, &s)| i64::from(c) * i64::from(s)).sum(); + let r = i64::from(*w.last()?) - prediction(p >> shift, lpc)?; if r.abs() > MAX_RESIDUAL { return None; } - out.push(r); + out.push(fold(r) as u32); } Some(out) } +/// The folded residual of linear prediction with `c` (any order), or None when it does not fit. +pub fn lpc_fold(x: &[i32], c: &[i32], shift: u8) -> Option> { + fold_with(x, c, shift, true) +} + +fn fold_with(x: &[i32], c: &[i32], shift: u8, lpc: bool) -> Option> { + match c.len() { + 0 => predict_fold::<0>(x, c, shift, lpc), + 1 => predict_fold::<1>(x, c, shift, lpc), + 2 => predict_fold::<2>(x, c, shift, lpc), + 3 => predict_fold::<3>(x, c, shift, lpc), + 4 => predict_fold::<4>(x, c, shift, lpc), + 5 => predict_fold::<5>(x, c, shift, lpc), + 6 => predict_fold::<6>(x, c, shift, lpc), + 7 => predict_fold::<7>(x, c, shift, lpc), + 8 => predict_fold::<8>(x, c, shift, lpc), + 9 => predict_fold::<9>(x, c, shift, lpc), + 10 => predict_fold::<10>(x, c, shift, lpc), + 11 => predict_fold::<11>(x, c, shift, lpc), + 12 => predict_fold::<12>(x, c, shift, lpc), + _ => predict_fold_any(x, c, shift, lpc), + } +} + +/// The fixed polynomial predictors of RFC 9639 §9.2.5 as integer coefficients. +const FIXED: [&[i32]; 5] = [&[], &[1], &[2, -1], &[3, -3, 1], &[4, -6, 4, -1]]; + +/// The folded residual of fixed predictor `order` (0–4), or None when it does not fit. +pub fn fixed_fold(x: &[i32], order: usize) -> Option> { + fold_with(x, FIXED.get(order)?, 0, false) +} + /// Tukey(0.5) window: flat in the middle half, raised-cosine tapers. -fn window(n: usize) -> Vec { +pub fn window(n: usize) -> Vec { let taper = (n / 4).max(1); (0..n) .map(|i| { @@ -166,13 +228,12 @@ fn window(n: usize) -> Vec { /// Linear prediction coefficients of every order 1..=`max` (Levinson–Durbin on the windowed /// autocorrelation); empty for a silent or degenerate block. -fn lpc_orders(x: &[i32], max: usize) -> Vec> { +fn lpc_orders(x: &[i32], max: usize, w: &[f64]) -> Vec> { let n = x.len(); - if n <= max { + if n <= max || w.len() != n { return Vec::new(); } - let w = window(n); - let s: Vec = x.iter().zip(&w).map(|(&v, &w)| f64::from(v) * w).collect(); + let s: Vec = x.iter().zip(w).map(|(&v, &w)| f64::from(v) * w).collect(); let r: Vec = (0..=max).map(|lag| s.iter().zip(s.iter().skip(lag)).map(|(a, b)| a * b).sum()).collect(); let Some(&r0) = r.first() else { return Vec::new() }; if r0 <= 0.0 { @@ -251,7 +312,7 @@ impl Effort { } /// The cheapest coding of one channel of `bps`-bit samples, and its size in bits. -pub fn choose(x: &[i32], bps: u8, effort: Effort) -> (Subframe, u64) { +pub fn choose(x: &[i32], bps: u8, effort: Effort, win: &[f64]) -> (Subframe, u64) { let n = x.len(); let b = u64::from(bps); let head = 8; // zero bit, type, wasted-bits flag @@ -263,7 +324,7 @@ pub fn choose(x: &[i32], bps: u8, effort: Effort) -> (Subframe, u64) { let wasted = x.iter().fold(0, |bits, &v| bits | v).trailing_zeros(); if wasted > 0 { let shifted: Vec = x.iter().map(|&v| v >> wasted).collect(); - let (coded, bits) = choose(&shifted, bps - wasted as u8, effort); + let (coded, bits) = choose(&shifted, bps - wasted as u8, effort, win); return (Subframe::Wasted { bits: wasted as u8, coded: Box::new(coded) }, bits + u64::from(wasted)); } let mut best = (Subframe::Verbatim, head + n as u64 * b); @@ -272,20 +333,19 @@ pub fn choose(x: &[i32], bps: u8, effort: Effort) -> (Subframe, u64) { best = (sub, bits); } }; - let folded = |r: &[i64]| r.iter().map(|&v| fold(v)).collect::>(); for order in 0..=4usize.min(n.saturating_sub(1)) { - let Some(r) = fixed_residual(x, order) else { continue }; - let residual = code_residual(&folded(&r), n, order, effort.max_partition_order); + let Some(r) = fixed_fold(x, order) else { continue }; + let residual = code_residual(r, n, order, effort.max_partition_order); let bits = head + order as u64 * b + residual.bits; consider(Subframe::Fixed { order: order as u8, residual }, bits); } // 15-bit coefficients, fewer for short blocks (their side information weighs more) let precision: u8 = if n <= 1152 { 12 } else { 15 }; - for coefs in lpc_orders(x, effort.max_lpc_order) { + for coefs in lpc_orders(x, effort.max_lpc_order, win) { let Some((q, shift)) = quantise(&coefs, precision) else { continue }; let order = q.len(); - let Some(r) = lpc_residual(x, &q, shift) else { continue }; - let residual = code_residual(&folded(&r), n, order, effort.max_partition_order); + let Some(r) = lpc_fold(x, &q, shift) else { continue }; + let residual = code_residual(r, n, order, effort.max_partition_order); let bits = head + order as u64 * (b + u64::from(precision)) + 4 + 5 + residual.bits; consider(Subframe::Lpc { coefs: q, precision, shift, residual }, bits); } @@ -305,34 +365,34 @@ mod tests { fn predictable_signals_beat_verbatim() { // a ramp is exactly predicted by the order-2 fixed predictor let ramp: Vec = (0..4096).map(|i| i * 3 - 6000).collect(); - let (sub, bits) = choose(&ramp, 16, Effort::level(5)); + let (sub, bits) = choose(&ramp, 16, Effort::level(5), &window(ramp.len())); assert!(matches!(sub, Subframe::Fixed { order: 2, .. } | Subframe::Lpc { .. }), "{sub:?}"); assert!(bits < 4096 * 2, "{bits}"); // a sine is well predicted by LPC let sine: Vec = (0..4096).map(|i| ((i as f64 * 0.05).sin() * 20000.0) as i32).collect(); - let (_, bits) = choose(&sine, 16, Effort::level(5)); + let (_, bits) = choose(&sine, 16, Effort::level(5), &window(sine.len())); assert!(bits < 4096 * 6, "{bits}"); - assert_eq!(choose(&[7; 100], 16, Effort::level(5)), (Subframe::Constant(7), 24)); + assert_eq!(choose(&[7; 100], 16, Effort::level(5), &window(100)), (Subframe::Constant(7), 24)); } #[test] fn residuals_may_reach_the_limit_of_a_32_bit_integer() { - assert!(fixed_residual(&[i32::MAX, -i32::MAX], 0).is_some()); - assert!(fixed_residual(&[0, i32::MAX], 1).is_some()); + assert!(fixed_fold(&[i32::MAX, -i32::MAX], 0).is_some()); + assert!(fixed_fold(&[0, i32::MAX], 1).is_some()); // the most negative value is out of range (§9.2.7.3) - assert!(fixed_residual(&[i32::MIN], 0).is_none()); - assert!(fixed_residual(&[0, 1 << 30, -(1 << 30)], 1).is_none()); - assert!(lpc_residual(&[0, i32::MAX], &[1], 0).is_some()); + assert!(fixed_fold(&[i32::MIN], 0).is_none()); + assert!(fixed_fold(&[0, 1 << 30, -(1 << 30)], 1).is_none()); + assert!(lpc_fold(&[0, i32::MAX], &[1], 0).is_some()); } #[test] fn linear_predictions_stay_in_the_range_of_a_32_bit_integer() { let x = [1 << 30, 1 << 30]; - assert!(lpc_residual(&x, &[1], 0).is_some()); + assert!(lpc_fold(&x, &[1], 0).is_some()); // predicts 2^31, and the residual would be small - assert!(lpc_residual(&x, &[2], 0).is_none()); - assert!(lpc_residual(&x, &[4], 1).is_none()); - assert!(lpc_residual(&x, &[3], 1).is_some()); + assert!(lpc_fold(&x, &[2], 0).is_none()); + assert!(lpc_fold(&x, &[4], 1).is_none()); + assert!(lpc_fold(&x, &[3], 1).is_some()); } #[test] @@ -340,13 +400,13 @@ mod tests { let x = tones(4096); let padded: Vec = x.iter().map(|&s| s << 8).collect(); let effort = Effort::level(5); - let (sub, bits) = choose(&x, 16, effort); + let (sub, bits) = choose(&x, 16, effort, &window(x.len())); // the same coding, plus the 8 bits that say how many are wasted - assert_eq!(choose(&padded, 24, effort), (Subframe::Wasted { bits: 8, coded: Box::new(sub) }, bits + 8)); + assert_eq!(choose(&padded, 24, effort, &window(padded.len())), (Subframe::Wasted { bits: 8, coded: Box::new(sub) }, bits + 8)); // at most what leaves one bit per sample, and nothing for silence, a constant keeps its bits let lowest = [0, -(1 << 23)]; - assert!(matches!(choose(&lowest, 24, effort).0, Subframe::Wasted { bits: 23, .. })); - assert_eq!(choose(&[256; 10], 24, effort).0, Subframe::Constant(256)); + assert!(matches!(choose(&lowest, 24, effort, &window(lowest.len())).0, Subframe::Wasted { bits: 23, .. })); + assert_eq!(choose(&[256; 10], 24, effort, &window([256; 10].len())).0, Subframe::Constant(256)); } /// Two tones and a little noise, from a fixed seed. @@ -366,8 +426,8 @@ mod tests { fn linear_prediction_beats_the_fixed_predictors() { let x = tones(4096); let effort = Effort::level(5); - let (_, lpc) = choose(&x, 16, effort); - let (_, fixed) = choose(&x, 16, Effort { max_lpc_order: 0, ..effort }); + let (_, lpc) = choose(&x, 16, effort, &window(x.len())); + let (_, fixed) = choose(&x, 16, Effort { max_lpc_order: 0, ..effort }, &window(x.len())); assert!(lpc < fixed * 4 / 5, "lpc {lpc} fixed {fixed}"); } @@ -376,17 +436,18 @@ mod tests { // quiet, then loud: one Rice parameter cannot suit both let x: Vec = tones(4096).iter().enumerate().map(|(i, &s)| if i < 2048 { s / 256 } else { s }).collect(); let effort = Effort::level(5); - let (_, partitioned) = choose(&x, 16, effort); - let (_, whole) = choose(&x, 16, Effort { max_partition_order: 0, ..effort }); + let (_, partitioned) = choose(&x, 16, effort, &window(x.len())); + let (_, whole) = choose(&x, 16, Effort { max_partition_order: 0, ..effort }, &window(x.len())); assert!(partitioned < whole * 19 / 20, "partitioned {partitioned} whole {whole}"); } #[test] fn lpc_residual_inverts() { let x: Vec = (0..500).map(|i| ((i as f64 * 0.1).sin() * 1000.0) as i32).collect(); - let coefs = lpc_orders(&x, 4).pop().unwrap(); + let coefs = lpc_orders(&x, 4, &window(x.len())).pop().unwrap(); let (q, shift) = quantise(&coefs, 15).unwrap(); - let r = lpc_residual(&x, &q, shift).unwrap(); + let r: Vec = + lpc_fold(&x, &q, shift).unwrap().iter().map(|&u| if u & 1 == 0 { i64::from(u >> 1) } else { -i64::from(u >> 1) - 1 }).collect(); // rebuild the signal from the warm-up and the residual, as a decoder would let mut y: Vec = x[..q.len()].to_vec(); for (i, &res) in r.iter().enumerate() { @@ -405,12 +466,12 @@ mod tests { (seed >> 8) as i32 - (1 << 23) }) .collect(); - let (_, bits) = choose(&noise, 24, Effort::level(8)); + let (_, bits) = choose(&noise, 24, Effort::level(8), &window(noise.len())); assert!(bits <= 8 + 4096 * 24); let square: Vec = (0..4096).map(|i| if i % 2 == 0 { (1 << 23) - 1 } else { -(1 << 23) }).collect(); - choose(&square, 24, Effort::level(8)); - choose(&[], 16, Effort::level(5)); - choose(&[1], 16, Effort::level(5)); - choose(&[1, 2], 16, Effort::level(8)); + choose(&square, 24, Effort::level(8), &window(square.len())); + choose(&[], 16, Effort::level(5), &window(0)); + choose(&[1], 16, Effort::level(5), &window(1)); + choose(&[1, 2], 16, Effort::level(8), &window(2)); } } diff --git a/crates/audio/src/write/pcm.rs b/crates/audio/src/write/pcm.rs index 8374f3a0..614d5a4a 100644 --- a/crates/audio/src/write/pcm.rs +++ b/crates/audio/src/write/pcm.rs @@ -43,10 +43,18 @@ impl Quantizer { v += self.uniform() - self.uniform(); } let max = scale - 1.0; - v.round().clamp(-scale, max) as i32 + round_half_away(v.clamp(-scale, max)) as i32 } } +/// `v.round()` for values in the `i32` range, without the libm call: the bounce's inner loop. +#[inline] +fn round_half_away(v: f64) -> f64 { + let t = v as i64; + let f = v - t as f64; + (t + i64::from(f >= 0.5) - i64::from(f <= -0.5)) as f64 +} + /// Bits per sample at `depth`, and whether they are float. #[cfg(feature = "pcm-encode")] pub(crate) fn sample_format(depth: BitDepth) -> (u16, bool) { @@ -122,3 +130,38 @@ pub(crate) fn validate_layout(channels: usize, sample_rate: u32) -> Result<()> { } Ok(()) } + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn round_half_away_matches_f64_round_on_the_quantizer_range() { + let edges = [ + 0.0, + -0.0, + 0.5, + -0.5, + 1.5, + -1.5, + 2.5, + 0.499_999_999_999_999_94, + -0.499_999_999_999_999_94, + 8_388_607.5, + -8_388_608.0, + 2_147_483_647.0, + -2_147_483_648.0, + ]; + let mut x = 0x2545_f491_4f6c_dd1du64; + let spread = (0..10_000).map(|_| { + x ^= x << 13; + x ^= x >> 7; + x ^= x << 17; + (x >> 11) as f64 / (1u64 << 53) as f64 * 4_294_967_295.0 - 2_147_483_648.0 + }); + for v in edges.into_iter().chain(spread) { + assert!(round_half_away(v) == v.round(), "{v}"); + assert_eq!(round_half_away(v) as i32, v.round() as i32, "{v}"); + } + } +} diff --git a/crates/audio/tests/audio/malformed.rs b/crates/audio/tests/audio/malformed.rs index a4feb37d..0c48968f 100644 --- a/crates/audio/tests/audio/malformed.rs +++ b/crates/audio/tests/audio/malformed.rs @@ -321,7 +321,8 @@ fn flac_metadata_fields_beyond_their_block_are_refused() { let picture = |data: u32| [&[0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0][..], &[0; 16], &data.to_be_bytes(), &[1, 2, 3]].concat(); assert_eq!(decode(&with_block(6, &picture(3)), None).unwrap().1.frames(), 300); let comment = [&0u32.to_le_bytes()[..], &1u32.to_le_bytes(), &u32::MAX.to_le_bytes()].concat(); - for f in [with_block(6, &picture(0x7fff_fff0)), with_block(4, &comment)] { + let vendor_overflow = [&0x7fff_fff0u32.to_le_bytes()[..], &0u32.to_le_bytes()].concat(); + for f in [with_block(6, &picture(0x7fff_fff0)), with_block(4, &comment), with_block(4, &vendor_overflow)] { assert!(matches!(decode(&f, None), Err(AudioError::Malformed(m)) if m.contains("larger than its block"))); } } diff --git a/crates/media/src/ogg/mod.rs b/crates/media/src/ogg/mod.rs index 1ab5f5ad..6311335f 100644 --- a/crates/media/src/ogg/mod.rs +++ b/crates/media/src/ogg/mod.rs @@ -66,19 +66,63 @@ pub fn sniff(head: &[u8]) -> bool { /// CRC-32 of RFC 3533 §6: polynomial 0x04C11DB7, initial value 0, no reflection, no final XOR. pub fn crc32(data: &[u8]) -> u32 { - static TABLE: std::sync::OnceLock<[u32; 256]> = std::sync::OnceLock::new(); - let t = TABLE.get_or_init(|| { - let mut t = [0u32; 256]; - for (i, e) in t.iter_mut().enumerate() { - let mut r = (i as u32) << 24; - for _ in 0..8 { - r = if r & 0x8000_0000 != 0 { (r << 1) ^ 0x04C1_1DB7 } else { r << 1 }; - } - *e = r; + crc32_update(0, data) +} + +/// `CRC_TABLES[k][b]`: the CRC register after feeding byte `b` followed by `k` zero bytes +/// (slice-by-16: sixteen independent lookups per 16 input bytes instead of a serial chain of 16). +static CRC_TABLES: [[u32; 256]; 16] = crc_tables(); + +const fn crc_tables() -> [[u32; 256]; 16] { + let mut t = [[0u32; 256]; 16]; + let mut i = 0; + while i < 256 { + let mut r = (i as u32) << 24; + let mut b = 0; + while b < 8 { + r = if r & 0x8000_0000 != 0 { (r << 1) ^ 0x04C1_1DB7 } else { r << 1 }; + b += 1; } - t - }); - data.iter().fold(0u32, |c, &b| (c << 8) ^ t[((c >> 24) as u8 ^ b) as usize]) + t[0][i] = r; + i += 1; + } + let mut k = 1; + while k < 16 { + let mut i = 0; + while i < 256 { + let p = t[k - 1][i]; + t[k][i] = (p << 8) ^ t[0][(p >> 24) as usize]; + i += 1; + } + k += 1; + } + t +} + +/// Continue an RFC 3533 CRC-32 over `data` from register `crc`. +pub fn crc32_update(mut crc: u32, data: &[u8]) -> u32 { + let t = &CRC_TABLES; + let (chunks, rest) = data.as_chunks::<16>(); + for c in chunks { + let w = crc ^ u32::from_be_bytes([c[0], c[1], c[2], c[3]]); + crc = t[15][(w >> 24) as usize] + ^ t[14][((w >> 16) & 0xFF) as usize] + ^ t[13][((w >> 8) & 0xFF) as usize] + ^ t[12][(w & 0xFF) as usize] + ^ t[11][c[4] as usize] + ^ t[10][c[5] as usize] + ^ t[9][c[6] as usize] + ^ t[8][c[7] as usize] + ^ t[7][c[8] as usize] + ^ t[6][c[9] as usize] + ^ t[5][c[10] as usize] + ^ t[4][c[11] as usize] + ^ t[3][c[12] as usize] + ^ t[2][c[13] as usize] + ^ t[1][c[14] as usize] + ^ t[0][c[15] as usize]; + } + rest.iter().fold(crc, |c, &b| (c << 8) ^ t[0][((c >> 24) as u8 ^ b) as usize]) } /// The codec of a logical bitstream (from its first packet). diff --git a/crates/media/src/ogg/tests.rs b/crates/media/src/ogg/tests.rs index ce906a7d..bc47bfc0 100644 --- a/crates/media/src/ogg/tests.rs +++ b/crates/media/src/ogg/tests.rs @@ -647,3 +647,52 @@ fn theora_codec_params_carry_the_picture_size() { let p = open(&f).unwrap().streams[0].codec_params(); assert_eq!((p.width, p.height), (Some(1920), Some(1080))); } + +/// Bit-at-a-time CRC of RFC 3533 §6 (polynomial 0x04C11DB7, initial value 0, no reflection, no final XOR). +fn crc32_bitwise(mut crc: u32, data: &[u8]) -> u32 { + for &b in data { + crc ^= u32::from(b) << 24; + for _ in 0..8 { + crc = if crc & 0x8000_0000 != 0 { (crc << 1) ^ 0x04C1_1DB7 } else { crc << 1 }; + } + } + crc +} + +#[test] +fn crc_matches_bitwise_reference_on_every_length_and_alignment() { + let mut seed = 7u64; + let buf: Vec = (0..70_100) + .map(|_| { + seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + (seed >> 33) as u8 + }) + .collect(); + // every short length (covers the 16-byte main loop and every tail length), at several start offsets + for start in 0..4 { + for len in 0..=64 { + let d = &buf[start..start + len]; + assert_eq!(crc32(d), crc32_bitwise(0, d), "start {start} len {len}"); + } + } + for len in [65, 255, 4096, 65_307, 70_000] { + for start in [0, 1, 3] { + let d = &buf[start..start + len]; + assert_eq!(crc32(d), crc32_bitwise(0, d), "start {start} len {len}"); + } + } + // all-zero and all-ones input + assert_eq!(crc32(&[0u8; 100]), 0); + assert_eq!(crc32(&[0xFFu8; 100]), crc32_bitwise(0, &[0xFFu8; 100])); +} + +#[test] +fn crc_update_continues_across_any_split() { + let data: Vec = (0..200u32).map(|i| (i * 31 + 5) as u8).collect(); + let whole = crc32(&data); + for cut in 0..=data.len() { + let (a, b) = data.split_at(cut); + assert_eq!(crc32_update(crc32_update(0, a), b), whole, "cut {cut}"); + assert_eq!(crc32_update(crc32_bitwise(0, a), b), whole, "cut {cut} (reference head)"); + } +}