From 1ebc34b0cf73d21001b846c00a5febd4368f9afc Mon Sep 17 00:00:00 2001 From: "Dylan M. Taylor" Date: Wed, 11 Nov 2026 17:26:39 -0500 Subject: [PATCH] craft-audio: MP3 (Layer III) encoder from FilmCraft (feature `mp3-encode`) (#869) Moves filmcraft-mp3enc (storytold/filmcraft@4bec22b) into craft_audio::mp3, on craft-media's BitWriter. `encode` turns the feature on, mp3::Error converts into AudioError::Encode, and a test covers every sample rate and bit rate with extreme signals. --- README.md | 2 +- crates/audio/Cargo.toml | 8 +- crates/audio/README.md | 35 +- crates/audio/src/lib.rs | 14 +- crates/audio/src/mp3/filterbank.rs | 111 ++++ crates/audio/src/mp3/huffman.rs | 228 ++++++++ crates/audio/src/mp3/huffman_tables.rs | 341 ++++++++++++ crates/audio/src/mp3/mod.rs | 502 ++++++++++++++++++ crates/audio/src/mp3/quantize.rs | 370 +++++++++++++ crates/audio/src/mp3/tables.rs | 125 +++++ crates/audio/src/mp3/tests.rs | 490 +++++++++++++++++ crates/audio/tests/codecs/main.rs | 1 + .../audio/tests/codecs/mp3_encoder_oracle.rs | 162 ++++++ 13 files changed, 2380 insertions(+), 9 deletions(-) create mode 100644 crates/audio/src/mp3/filterbank.rs create mode 100644 crates/audio/src/mp3/huffman.rs create mode 100644 crates/audio/src/mp3/huffman_tables.rs create mode 100644 crates/audio/src/mp3/mod.rs create mode 100644 crates/audio/src/mp3/quantize.rs create mode 100644 crates/audio/src/mp3/tables.rs create mode 100644 crates/audio/src/mp3/tests.rs create mode 100644 crates/audio/tests/codecs/mp3_encoder_oracle.rs diff --git a/README.md b/README.md index 4fc478c6..d592ad3d 100644 --- a/README.md +++ b/README.md @@ -131,7 +131,7 @@ craft- = { git = "https://github.com/storytold/craft-libs", rev = " for AudioError { } } +#[cfg(feature = "mp3-encode")] +impl From for AudioError { + fn from(e: mp3::Error) -> Self { + AudioError::Encode(e.to_string()) + } +} + #[cfg(feature = "flac-encode")] impl From for AudioError { fn from(e: flac::Error) -> Self { diff --git a/crates/audio/src/mp3/filterbank.rs b/crates/audio/src/mp3/filterbank.rs new file mode 100644 index 00000000..da8262aa --- /dev/null +++ b/crates/audio/src/mp3/filterbank.rs @@ -0,0 +1,111 @@ +//! The Layer III hybrid analysis filterbank: the polyphase analysis subband filter of +//! ISO/IEC 11172-3 Annex C.1.3 (32 subbands), then per subband an 18-line MDCT with the normal +//! (block type 0) window over two granules, and the encoder side of the alias-reduction +//! butterflies (the inverse of §2.4.3.4.10.4's decoder butterflies). + +use std::sync::OnceLock; + +use super::tables::{alias_coefficients, analysis_matrix, analysis_window}; + +/// Frequency lines per granule. +pub(crate) const GRANULE: usize = 576; + +/// The polyphase analysis filter of one channel: a 512-sample FIFO, newest sample first. +#[derive(Clone)] +pub(crate) struct Polyphase { + x: [f32; 512], +} + +impl Polyphase { + pub fn new() -> Self { + Polyphase { x: [0.0; 512] } + } + + /// Shift in 32 samples (oldest first) and compute the 32 subband samples (Annex C.1.3). + pub fn filter(&mut self, input: &[f32; 32], out: &mut [f32; 32]) { + self.x.copy_within(0..480, 32); + for (i, s) in input.iter().enumerate() { + // X[31] receives the first of the new samples, X[0] the last + self.x[31 - i] = *s; + } + let c = analysis_window(); + let mut y = [0f32; 64]; + for (i, yi) in y.iter_mut().enumerate() { + let mut acc = 0f32; + for j in 0..8 { + acc += c[i + 64 * j] * self.x[i + 64 * j]; + } + *yi = acc; + } + for (o, row) in out.iter_mut().zip(analysis_matrix()) { + *o = row.iter().zip(&y).map(|(m, y)| m * y).sum(); + } + } +} + +/// The MDCT kernel for normal blocks: `w[i] cos(π/72 (2i + 1 + 18)(2k + 1)) / 9`, `k` in `0..18`, +/// `i` in `0..36`, with the sine window `w[i] = sin(π/36 (i + 1/2))`. The 1/9 (= 2/18) makes it the +/// exact inverse of the decoder's unscaled IMDCT with windowed overlap-add. +fn mdct_kernel() -> &'static [[f32; 36]; 18] { + static K: OnceLock<[[f32; 36]; 18]> = OnceLock::new(); + K.get_or_init(|| { + let mut k = [[0f32; 36]; 18]; + for (kk, row) in k.iter_mut().enumerate() { + for (i, c) in row.iter_mut().enumerate() { + let w = (std::f64::consts::PI / 36.0 * (i as f64 + 0.5)).sin(); + let a = std::f64::consts::PI / 72.0 * (2 * i + 1 + 18) as f64 * (2 * kk + 1) as f64; + *c = (w * a.cos() / 9.0) as f32; + } + } + k + }) +} + +/// The hybrid filterbank of one channel: the polyphase filter and the previous granule's subband +/// samples for the MDCT overlap. +#[derive(Clone)] +pub(crate) struct Analysis { + poly: Polyphase, + prev: [[f32; 18]; 32], +} + +impl Analysis { + pub fn new() -> Self { + Analysis { poly: Polyphase::new(), prev: [[0.0; 18]; 32] } + } + + /// Transform one granule of 576 samples into 576 frequency lines (`xr[18 * sb + k]`). + pub fn granule(&mut self, pcm: &[f32; GRANULE], xr: &mut [f32; GRANULE]) { + let mut cur = [[0f32; 18]; 32]; + let mut slot = [0f32; 32]; + let mut out = [0f32; 32]; + for (t, chunk) in pcm.as_chunks::<32>().0.iter().enumerate().take(18) { + slot.copy_from_slice(chunk); + self.poly.filter(&slot, &mut out); + for (sb, v) in out.iter().enumerate() { + // frequency inversion: odd subbands, odd time slots (§2.4.3.4.10.5 mirrored) + cur[sb][t] = if sb % 2 == 1 && t % 2 == 1 { -v } else { *v }; + } + } + let kernel = mdct_kernel(); + for sb in 0..32 { + let mut z = [0f32; 36]; + z[..18].copy_from_slice(&self.prev[sb]); + z[18..].copy_from_slice(&cur[sb]); + for (k, row) in kernel.iter().enumerate() { + xr[18 * sb + k] = row.iter().zip(&z).map(|(a, b)| a * b).sum(); + } + } + self.prev = cur; + // alias reduction: the inverse of the decoder's butterflies between adjacent subbands + for sb in 1..32 { + for (i, (cs, ca)) in alias_coefficients().iter().enumerate() { + let up = 18 * sb - 1 - i; + let dn = 18 * sb + i; + let (a, b) = (xr[up], xr[dn]); + xr[up] = a * cs + b * ca; + xr[dn] = b * cs - a * ca; + } + } + } +} diff --git a/crates/audio/src/mp3/huffman.rs b/crates/audio/src/mp3/huffman.rs new file mode 100644 index 00000000..15bc559b --- /dev/null +++ b/crates/audio/src/mp3/huffman.rs @@ -0,0 +1,228 @@ +//! Huffman coding of a granule's quantized lines (ISO/IEC 11172-3 §2.4.2.7 `huffmancodebits`, +//! §2.4.3.4.6): the big-values region in pairs, split into three regions at scalefactor band +//! boundaries with a table each, then the count1 region in quadruples of magnitude ≤ 1, then the +//! zeros. The encoder picks the region split and tables that need the fewest bits. + +use craft_media::bitstream::BitWriter; + +use super::filterbank::GRANULE; +use super::huffman_tables::{LINBITS, QUAD_A_CODES, QUAD_A_LENS, QUAD_B_CODES, pair_table}; + +/// The tables without `linbits`, smallest first. +const SMALL: [usize; 13] = [1, 2, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 15]; +/// Escape tables: 16-23 share table 16's codes, 24-31 table 24's. +const ESC16: usize = 16; +const ESC24: usize = 24; + +/// How a granule's lines are coded, as signalled in the side information. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub(crate) struct Coding { + /// Pairs in the big-values region. + pub big_values: usize, + /// Quadruples in the count1 region. + pub count1: usize, + pub table_select: [usize; 3], + pub region0_count: usize, + pub region1_count: usize, + /// 0: table A, 1: table B. + pub count1_table: usize, + /// Huffman-coded bits (part 3 of `part2_3_length`). + pub bits: usize, +} + +/// Cumulative bits per table over the scalefactor bands of the big-values region: entry `b` sums +/// bands `0..b`. +struct Prefix { + small: [[u32; 23]; SMALL.len()], + esc16: [u32; 23], + esc24: [u32; 23], + escapes: [u32; 23], + max: [u16; 22], +} + +/// Bits needed for `v` (0 for 0). +fn bit_width(v: u32) -> u32 { + 32 - v.leading_zeros() +} + +impl Prefix { + fn new(ix: &[u16], bands: &[u16; 23], end: usize) -> (Prefix, usize) { + let mut p = Prefix { small: [[0; 23]; SMALL.len()], esc16: [0; 23], esc24: [0; 23], escapes: [0; 23], max: [0; 22] }; + let tables: Vec<_> = SMALL.iter().map(|&t| pair_table(t)).collect(); + let (t16, t24) = (pair_table(ESC16), pair_table(ESC24)); + let mut nb = 0; + for (b, w) in bands.windows(2).enumerate() { + let (lo, hi) = (usize::from(w[0]), usize::from(w[1]).min(end)); + if lo >= end || b >= 22 { + break; + } + let lines = ix.get(lo..hi).unwrap_or(&[]); + let max = lines.iter().copied().max().unwrap_or(0); + p.max[b] = max; + for (n, t) in tables.iter().enumerate() { + let bits = match t.as_ref().filter(|t| usize::from(max) < t.size) { + Some(t) => lines + .as_chunks::<2>() + .0 + .iter() + .map(|q| { + let (x, y) = (usize::from(q[0]), usize::from(q[1])); + u32::from(t.lens.get(x * t.size + y).copied().unwrap_or(32)) + u32::from(x != 0) + u32::from(y != 0) + }) + .sum(), + None => u32::MAX / 32, + }; + p.small[n][b + 1] = p.small[n][b].saturating_add(bits); + } + let (mut e16, mut e24, mut esc) = (0u32, 0u32, 0u32); + for q in lines.as_chunks::<2>().0 { + let (x, y) = (usize::from(q[0]), usize::from(q[1])); + let k = x.min(15) * 16 + y.min(15); + let signs = u32::from(x != 0) + u32::from(y != 0); + e16 += u32::from(t16.as_ref().and_then(|t| t.lens.get(k)).copied().unwrap_or(32)) + signs; + e24 += u32::from(t24.as_ref().and_then(|t| t.lens.get(k)).copied().unwrap_or(32)) + signs; + esc += u32::from(x >= 15) + u32::from(y >= 15); + } + p.esc16[b + 1] = p.esc16[b] + e16; + p.esc24[b + 1] = p.esc24[b] + e24; + p.escapes[b + 1] = p.escapes[b] + esc; + nb = b + 1; + } + (p, nb) + } + + /// The cheapest table for bands `a..b`: (table_select, bits). + fn best(&self, a: usize, b: usize) -> (usize, u32) { + let max = self.max.get(a..b).and_then(|m| m.iter().copied().max()).unwrap_or(0); + if max == 0 { + return (0, 0); + } + let mut best = (0usize, u32::MAX); + for (n, &t) in SMALL.iter().enumerate() { + let sum = self.small[n][b].saturating_sub(self.small[n][a]); + if usize::from(max) < pair_table(t).map_or(0, |p| p.size) && sum < best.1 { + best = (t, sum); + } + } + // the escape tables: the smallest linbits that hold the largest value + let need = if max >= 15 { bit_width(u32::from(max) - 15) } else { 0 }; + let escapes = self.escapes[b] - self.escapes[a]; + for (family, base) in [(ESC16, self.esc16[b] - self.esc16[a]), (ESC24, self.esc24[b] - self.esc24[a])] { + if let Some(t) = (family..family + 8).find(|&t| LINBITS.get(t).is_some_and(|&l| l >= need)) { + let bits = base.saturating_add(escapes.saturating_mul(LINBITS.get(t).copied().unwrap_or(0))); + if bits < best.1 { + best = (t, bits); + } + } + } + best + } +} + +/// Choose the regions and tables for `ix` (magnitudes ≤ [`super::tables::MAX_QUANT`]) with +/// scalefactor band boundaries `bands`. +pub(crate) fn choose(ix: &[u16; GRANULE], bands: &[u16; 23]) -> Coding { + let mut coding = Coding::default(); + let Some(last) = ix.iter().rposition(|&v| v != 0) else { + return coding; + }; + let n = (last + 2) & !1; // round up to a whole pair + let mut i = n; + while i >= 4 && ix[i - 4..i].iter().all(|&v| v <= 1) { + i -= 4; + } + // count1 region [i, n) + let (mut bits_a, mut bits_b, mut signs) = (0u32, 0u32, 0u32); + for q in ix[i..n].as_chunks::<4>().0 { + let idx = usize::from(q[0] * 8 + q[1] * 4 + q[2] * 2 + q[3]); + bits_a += u32::from(QUAD_A_LENS.get(idx).copied().unwrap_or(6)); + bits_b += 4; + signs += q.iter().map(|&v| u32::from(v)).sum::(); + } + coding.count1 = (n - i) / 4; + coding.count1_table = usize::from(bits_b < bits_a); + let count1_bits = bits_a.min(bits_b) + signs; + coding.big_values = i / 2; + // big-values region [0, i): costs per scalefactor band, clipped at i + let (prefix, nb) = Prefix::new(ix, bands, i); + let mut best: Option<(u32, usize, usize, [usize; 3])> = None; + if nb > 0 { + for r0 in 0..16usize { + for r1 in 0..8usize { + if r0 + r1 + 2 > 22 { + continue; + } + let s1 = (r0 + 1).min(nb); + let s2 = (r0 + r1 + 2).min(nb); + let (t0, b0) = prefix.best(0, s1); + let (t1, b1) = prefix.best(s1, s2); + let (t2, b2) = prefix.best(s2, nb); + let total = b0.saturating_add(b1).saturating_add(b2); + if best.is_none_or(|b| total < b.0) { + best = Some((total, r0, r1, [t0, t1, t2])); + } + if s2 >= nb && r1 > 0 { + break; // longer region 1 changes nothing + } + } + if r0 + 1 >= nb { + break; + } + } + } + let (big_bits, r0, r1, tables) = best.unwrap_or((0, 0, 0, [0; 3])); + coding.region0_count = r0; + coding.region1_count = r1; + coding.table_select = tables; + coding.bits = big_bits as usize + count1_bits as usize; + coding +} + +/// Write the Huffman-coded lines; `negative[i]` is the sign of line `i`. +pub(crate) fn write(w: &mut BitWriter, ix: &[u16; GRANULE], negative: &[bool; GRANULE], c: &Coding, bands: &[u16; 23]) { + let r1 = bands.get(c.region0_count + 1).map_or(GRANULE, |&b| usize::from(b)); + let r2 = bands.get(c.region0_count + c.region1_count + 2).map_or(GRANULE, |&b| usize::from(b)); + let end = (2 * c.big_values).min(GRANULE); + for i in (0..end).step_by(2) { + let select = c.table_select[if i < r1 { + 0 + } else if i < r2 { + 1 + } else { + 2 + }]; + let Some(t) = pair_table(select) else { continue }; + let linbits = LINBITS.get(select).copied().unwrap_or(0); + let (x, y) = (usize::from(ix[i]), usize::from(ix[i + 1])); + let (xe, ye) = if linbits > 0 { (x.min(15), y.min(15)) } else { (x, y) }; + let k = xe * t.size + ye; + w.write_bits(u32::from(t.codes.get(k).copied().unwrap_or(0)), u32::from(t.lens.get(k).copied().unwrap_or(0))); + if linbits > 0 && xe == 15 { + w.write_bits((x - 15) as u32, linbits); + } + if x != 0 { + w.write_bit(negative[i]); + } + if linbits > 0 && ye == 15 { + w.write_bits((y - 15) as u32, linbits); + } + if y != 0 { + w.write_bit(negative[i + 1]); + } + } + let c1_end = (end + 4 * c.count1).min(GRANULE); + for i in (end..c1_end).step_by(4) { + let Some(q) = ix.get(i..i + 4) else { break }; + let idx = usize::from(q[0] * 8 + q[1] * 4 + q[2] * 2 + q[3]) & 15; + if c.count1_table == 0 { + w.write_bits(u32::from(QUAD_A_CODES[idx]), u32::from(QUAD_A_LENS[idx])); + } else { + w.write_bits(u32::from(QUAD_B_CODES[idx]), 4); + } + for (k, &v) in q.iter().enumerate() { + if v != 0 { + w.write_bit(negative[i + k]); + } + } + } +} diff --git a/crates/audio/src/mp3/huffman_tables.rs b/crates/audio/src/mp3/huffman_tables.rs new file mode 100644 index 00000000..b98290e0 --- /dev/null +++ b/crates/audio/src/mp3/huffman_tables.rs @@ -0,0 +1,341 @@ +//! The Huffman code tables of Layer III (ISO/IEC 11172-3 Annex B, Table 3-B.7): `hcod` and +//! `hlen` per value pair `(x, y)` at index `x * size + y`, and the two count1 quadruple tables. +//! Tables 4 and 14 are not used by the standard; tables 16-23 share the codes of table 16 and +//! 24-31 those of table 24, each with its own number of `linbits`. + +/// One pair table: `size × size` codes for values `0..size` (the escape value 15 of tables 16 and +/// 24 is followed by `linbits` bits of `value - 15`). +pub(crate) struct PairTable { + pub size: usize, + pub codes: &'static [u16], + pub lens: &'static [u8], +} + +const T1C: [u16; 4] = [1, 1, 1, 0]; +const T1L: [u8; 4] = [1, 3, 2, 3]; + +const T2C: [u16; 9] = [1, 2, 1, 3, 1, 1, 3, 2, 0]; +const T2L: [u8; 9] = [1, 3, 6, 3, 3, 5, 5, 5, 6]; + +const T3C: [u16; 9] = [3, 2, 1, 1, 1, 1, 3, 2, 0]; +const T3L: [u8; 9] = [2, 2, 6, 3, 2, 5, 5, 5, 6]; + +const T5C: [u16; 16] = [1, 2, 6, 5, 3, 1, 4, 4, 7, 5, 7, 1, 6, 1, 1, 0]; +const T5L: [u8; 16] = [1, 3, 6, 7, 3, 3, 6, 7, 6, 6, 7, 8, 7, 6, 7, 8]; + +const T6C: [u16; 16] = [7, 3, 5, 1, 6, 2, 3, 2, 5, 4, 4, 1, 3, 3, 2, 0]; +const T6L: [u8; 16] = [3, 3, 5, 7, 3, 2, 4, 5, 4, 4, 5, 6, 6, 5, 6, 7]; + +#[rustfmt::skip] +const T7C: [u16; 36] = [ + 1, 2, 10, 19, 16, 10, + 3, 3, 7, 10, 5, 3, + 11, 4, 13, 17, 8, 4, + 12, 11, 18, 15, 11, 2, + 7, 6, 9, 14, 3, 1, + 6, 4, 5, 3, 2, 0, +]; +#[rustfmt::skip] +const T7L: [u8; 36] = [ + 1, 3, 6, 8, 8, 9, + 3, 4, 6, 7, 7, 8, + 6, 5, 7, 8, 8, 9, + 7, 7, 8, 9, 9, 9, + 7, 7, 8, 9, 9, 10, + 8, 8, 9, 10, 10, 10, +]; + +#[rustfmt::skip] +const T8C: [u16; 36] = [ + 3, 4, 6, 18, 12, 5, + 5, 1, 2, 16, 9, 3, + 7, 3, 5, 14, 7, 3, + 19, 17, 15, 13, 10, 4, + 13, 5, 8, 11, 5, 1, + 12, 4, 4, 1, 1, 0, +]; +#[rustfmt::skip] +const T8L: [u8; 36] = [ + 2, 3, 6, 8, 8, 9, + 3, 2, 4, 8, 8, 8, + 6, 4, 6, 8, 8, 9, + 8, 8, 8, 9, 9, 10, + 8, 7, 8, 9, 10, 10, + 9, 8, 9, 9, 11, 11, +]; + +#[rustfmt::skip] +const T9C: [u16; 36] = [ + 7, 5, 9, 14, 15, 7, + 6, 4, 5, 5, 6, 7, + 7, 6, 8, 8, 8, 5, + 15, 6, 9, 10, 5, 1, + 11, 7, 9, 6, 4, 1, + 14, 4, 6, 2, 6, 0, +]; +#[rustfmt::skip] +const T9L: [u8; 36] = [ + 3, 3, 5, 6, 8, 9, + 3, 3, 4, 5, 6, 8, + 4, 4, 5, 6, 7, 8, + 6, 5, 6, 7, 7, 8, + 7, 6, 7, 7, 8, 9, + 8, 7, 8, 8, 9, 9, +]; + +#[rustfmt::skip] +const T10C: [u16; 64] = [ + 1, 2, 10, 23, 35, 30, 12, 17, + 3, 3, 8, 12, 18, 21, 12, 7, + 11, 9, 15, 21, 32, 40, 19, 6, + 14, 13, 22, 34, 46, 23, 18, 7, + 20, 19, 33, 47, 27, 22, 9, 3, + 31, 22, 41, 26, 21, 20, 5, 3, + 14, 13, 10, 11, 16, 6, 5, 1, + 9, 8, 7, 8, 4, 4, 2, 0, +]; +#[rustfmt::skip] +const T10L: [u8; 64] = [ + 1, 3, 6, 8, 9, 9, 9, 10, + 3, 4, 6, 7, 8, 9, 8, 8, + 6, 6, 7, 8, 9, 10, 9, 9, + 7, 7, 8, 9, 10, 10, 9, 10, + 8, 8, 9, 10, 10, 10, 10, 10, + 9, 9, 10, 10, 11, 11, 10, 11, + 8, 8, 9, 10, 10, 10, 11, 11, + 9, 8, 9, 10, 10, 11, 11, 11, +]; + +#[rustfmt::skip] +const T11C: [u16; 64] = [ + 3, 4, 10, 24, 34, 33, 21, 15, + 5, 3, 4, 10, 32, 17, 11, 10, + 11, 7, 13, 18, 30, 31, 20, 5, + 25, 11, 19, 59, 27, 18, 12, 5, + 35, 33, 31, 58, 30, 16, 7, 5, + 28, 26, 32, 19, 17, 15, 8, 14, + 14, 12, 9, 13, 14, 9, 4, 1, + 11, 4, 6, 6, 6, 3, 2, 0, +]; +#[rustfmt::skip] +const T11L: [u8; 64] = [ + 2, 3, 5, 7, 8, 9, 8, 9, + 3, 3, 4, 6, 8, 8, 7, 8, + 5, 5, 6, 7, 8, 9, 8, 8, + 7, 6, 7, 9, 8, 10, 8, 9, + 8, 8, 8, 9, 9, 10, 9, 10, + 8, 8, 9, 10, 10, 11, 10, 11, + 8, 7, 7, 8, 9, 10, 10, 10, + 8, 7, 8, 9, 10, 10, 10, 10, +]; + +#[rustfmt::skip] +const T12C: [u16; 64] = [ + 9, 6, 16, 33, 41, 39, 38, 26, + 7, 5, 6, 9, 23, 16, 26, 11, + 17, 7, 11, 14, 21, 30, 10, 7, + 17, 10, 15, 12, 18, 28, 14, 5, + 32, 13, 22, 19, 18, 16, 9, 5, + 40, 17, 31, 29, 17, 13, 4, 2, + 27, 12, 11, 15, 10, 7, 4, 1, + 27, 12, 8, 12, 6, 3, 1, 0, +]; +#[rustfmt::skip] +const T12L: [u8; 64] = [ + 4, 3, 5, 7, 8, 9, 9, 9, + 3, 3, 4, 5, 7, 7, 8, 8, + 5, 4, 5, 6, 7, 8, 7, 8, + 6, 5, 6, 6, 7, 8, 8, 8, + 7, 6, 7, 7, 8, 8, 8, 9, + 8, 7, 8, 8, 8, 9, 8, 9, + 8, 7, 7, 8, 8, 9, 9, 10, + 9, 8, 8, 9, 9, 9, 9, 10, +]; + +#[rustfmt::skip] +const T13C: [u16; 256] = [ + 1, 5, 14, 21, 34, 51, 46, 71, 42, 52, 68, 52, 67, 44, 43, 19, + 3, 4, 12, 19, 31, 26, 44, 33, 31, 24, 32, 24, 31, 35, 22, 14, + 15, 13, 23, 36, 59, 49, 77, 65, 29, 40, 30, 40, 27, 33, 42, 16, + 22, 20, 37, 61, 56, 79, 73, 64, 43, 76, 56, 37, 26, 31, 25, 14, + 35, 16, 60, 57, 97, 75, 114, 91, 54, 73, 55, 41, 48, 53, 23, 24, + 58, 27, 50, 96, 76, 70, 93, 84, 77, 58, 79, 29, 74, 49, 41, 17, + 47, 45, 78, 74, 115, 94, 90, 79, 69, 83, 71, 50, 59, 38, 36, 15, + 72, 34, 56, 95, 92, 85, 91, 90, 86, 73, 77, 65, 51, 44, 43, 42, + 43, 20, 30, 44, 55, 78, 72, 87, 78, 61, 46, 54, 37, 30, 20, 16, + 53, 25, 41, 37, 44, 59, 54, 81, 66, 76, 57, 54, 37, 18, 39, 11, + 35, 33, 31, 57, 42, 82, 72, 80, 47, 58, 55, 21, 22, 26, 38, 22, + 53, 25, 23, 38, 70, 60, 51, 36, 55, 26, 34, 23, 27, 14, 9, 7, + 34, 32, 28, 39, 49, 75, 30, 52, 48, 40, 52, 28, 18, 17, 9, 5, + 45, 21, 34, 64, 56, 50, 49, 45, 31, 19, 12, 15, 10, 7, 6, 3, + 48, 23, 20, 39, 36, 35, 53, 21, 16, 23, 13, 10, 6, 1, 4, 2, + 16, 15, 17, 27, 25, 20, 29, 11, 17, 12, 16, 8, 1, 1, 0, 1, +]; +#[rustfmt::skip] +const T13L: [u8; 256] = [ + 1, 4, 6, 7, 8, 9, 9, 10, 9, 10, 11, 11, 12, 12, 13, 13, + 3, 4, 6, 7, 8, 8, 9, 9, 9, 9, 10, 10, 11, 12, 12, 12, + 6, 6, 7, 8, 9, 9, 10, 10, 9, 10, 10, 11, 11, 12, 13, 13, + 7, 7, 8, 9, 9, 10, 10, 10, 10, 11, 11, 11, 11, 12, 13, 13, + 8, 7, 9, 9, 10, 10, 11, 11, 10, 11, 11, 12, 12, 13, 13, 14, + 9, 8, 9, 10, 10, 10, 11, 11, 11, 11, 12, 11, 13, 13, 14, 14, + 9, 9, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 13, 13, 14, 14, + 10, 9, 10, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 14, 16, 16, + 9, 8, 9, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 14, 15, 15, + 10, 9, 10, 10, 11, 11, 11, 13, 12, 13, 13, 14, 14, 14, 16, 15, + 10, 10, 10, 11, 11, 12, 12, 13, 12, 13, 14, 13, 14, 15, 16, 17, + 11, 10, 10, 11, 12, 12, 12, 12, 13, 13, 13, 14, 15, 15, 15, 16, + 11, 11, 11, 12, 12, 13, 12, 13, 14, 14, 15, 15, 15, 16, 16, 16, + 12, 11, 12, 13, 13, 13, 14, 14, 14, 14, 14, 15, 16, 15, 16, 16, + 13, 12, 12, 13, 13, 13, 15, 14, 14, 17, 15, 15, 15, 17, 16, 16, + 12, 12, 13, 14, 14, 14, 15, 14, 15, 15, 16, 16, 19, 18, 19, 16, +]; + +#[rustfmt::skip] +const T15C: [u16; 256] = [ + 7, 12, 18, 53, 47, 76, 124, 108, 89, 123, 108, 119, 107, 81, 122, 63, + 13, 5, 16, 27, 46, 36, 61, 51, 42, 70, 52, 83, 65, 41, 59, 36, + 19, 17, 15, 24, 41, 34, 59, 48, 40, 64, 50, 78, 62, 80, 56, 33, + 29, 28, 25, 43, 39, 63, 55, 93, 76, 59, 93, 72, 54, 75, 50, 29, + 52, 22, 42, 40, 67, 57, 95, 79, 72, 57, 89, 69, 49, 66, 46, 27, + 77, 37, 35, 66, 58, 52, 91, 74, 62, 48, 79, 63, 90, 62, 40, 38, + 125, 32, 60, 56, 50, 92, 78, 65, 55, 87, 71, 51, 73, 51, 70, 30, + 109, 53, 49, 94, 88, 75, 66, 122, 91, 73, 56, 42, 64, 44, 21, 25, + 90, 43, 41, 77, 73, 63, 56, 92, 77, 66, 47, 67, 48, 53, 36, 20, + 71, 34, 67, 60, 58, 49, 88, 76, 67, 106, 71, 54, 38, 39, 23, 15, + 109, 53, 51, 47, 90, 82, 58, 57, 48, 72, 57, 41, 23, 27, 62, 9, + 86, 42, 40, 37, 70, 64, 52, 43, 70, 55, 42, 25, 29, 18, 11, 11, + 118, 68, 30, 55, 50, 46, 74, 65, 49, 39, 24, 16, 22, 13, 14, 7, + 91, 44, 39, 38, 34, 63, 52, 45, 31, 52, 28, 19, 14, 8, 9, 3, + 123, 60, 58, 53, 47, 43, 32, 22, 37, 24, 17, 12, 15, 10, 2, 1, + 71, 37, 34, 30, 28, 20, 17, 26, 21, 16, 10, 6, 8, 6, 2, 0, +]; +#[rustfmt::skip] +const T15L: [u8; 256] = [ + 3, 4, 5, 7, 7, 8, 9, 9, 9, 10, 10, 11, 11, 11, 12, 13, + 4, 3, 5, 6, 7, 7, 8, 8, 8, 9, 9, 10, 10, 10, 11, 11, + 5, 5, 5, 6, 7, 7, 8, 8, 8, 9, 9, 10, 10, 11, 11, 11, + 6, 6, 6, 7, 7, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, + 7, 6, 7, 7, 8, 8, 9, 9, 9, 9, 10, 10, 10, 11, 11, 11, + 8, 7, 7, 8, 8, 8, 9, 9, 9, 9, 10, 10, 11, 11, 11, 12, + 9, 7, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 11, 11, 12, 12, + 9, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 12, + 9, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 12, 12, 12, + 9, 8, 9, 9, 9, 9, 10, 10, 10, 11, 11, 11, 11, 12, 12, 12, + 10, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 11, 12, 13, 12, + 10, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 11, 12, 12, 12, 13, + 11, 10, 9, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 13, 13, + 11, 10, 10, 10, 10, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, + 12, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 12, 13, + 12, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, +]; + +#[rustfmt::skip] +const T16C: [u16; 256] = [ + 1, 5, 14, 44, 74, 63, 110, 93, 172, 149, 138, 242, 225, 195, 376, 17, + 3, 4, 12, 20, 35, 62, 53, 47, 83, 75, 68, 119, 201, 107, 207, 9, + 15, 13, 23, 38, 67, 58, 103, 90, 161, 72, 127, 117, 110, 209, 206, 16, + 45, 21, 39, 69, 64, 114, 99, 87, 158, 140, 252, 212, 199, 387, 365, 26, + 75, 36, 68, 65, 115, 101, 179, 164, 155, 264, 246, 226, 395, 382, 362, 9, + 66, 30, 59, 56, 102, 185, 173, 265, 142, 253, 232, 400, 388, 378, 445, 16, + 111, 54, 52, 100, 184, 178, 160, 133, 257, 244, 228, 217, 385, 366, 715, 10, + 98, 48, 91, 88, 165, 157, 148, 261, 248, 407, 397, 372, 380, 889, 884, 8, + 85, 84, 81, 159, 156, 143, 260, 249, 427, 401, 392, 383, 727, 713, 708, 7, + 154, 76, 73, 141, 131, 256, 245, 426, 406, 394, 384, 735, 359, 710, 352, 11, + 139, 129, 67, 125, 247, 233, 229, 219, 393, 743, 737, 720, 885, 882, 439, 4, + 243, 120, 118, 115, 227, 223, 396, 746, 742, 736, 721, 712, 706, 223, 436, 6, + 202, 224, 222, 218, 216, 389, 386, 381, 364, 888, 443, 707, 440, 437, 1728, 4, + 747, 211, 210, 208, 370, 379, 734, 723, 714, 1735, 883, 877, 876, 3459, 865, 2, + 377, 369, 102, 187, 726, 722, 358, 711, 709, 866, 1734, 871, 3458, 870, 434, 0, + 12, 10, 7, 11, 10, 17, 11, 9, 13, 12, 10, 7, 5, 3, 1, 3, +]; +#[rustfmt::skip] +const T16L: [u8; 256] = [ + 1, 4, 6, 8, 9, 9, 10, 10, 11, 11, 11, 12, 12, 12, 13, 9, + 3, 4, 6, 7, 8, 9, 9, 9, 10, 10, 10, 11, 12, 11, 12, 8, + 6, 6, 7, 8, 9, 9, 10, 10, 11, 10, 11, 11, 11, 12, 12, 9, + 8, 7, 8, 9, 9, 10, 10, 10, 11, 11, 12, 12, 12, 13, 13, 10, + 9, 8, 9, 9, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 9, + 9, 8, 9, 9, 10, 11, 11, 12, 11, 12, 12, 13, 13, 13, 14, 10, + 10, 9, 9, 10, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 14, 10, + 10, 9, 10, 10, 11, 11, 11, 12, 12, 13, 13, 13, 13, 15, 15, 10, + 10, 10, 10, 11, 11, 11, 12, 12, 13, 13, 13, 13, 14, 14, 14, 10, + 11, 10, 10, 11, 11, 12, 12, 13, 13, 13, 13, 14, 13, 14, 13, 11, + 11, 11, 10, 11, 12, 12, 12, 12, 13, 14, 14, 14, 15, 15, 14, 10, + 12, 11, 11, 11, 12, 12, 13, 14, 14, 14, 14, 14, 14, 13, 14, 11, + 12, 12, 12, 12, 12, 13, 13, 13, 13, 15, 14, 14, 14, 14, 16, 11, + 14, 12, 12, 12, 13, 13, 14, 14, 14, 16, 15, 15, 15, 17, 15, 11, + 13, 13, 11, 12, 14, 14, 13, 14, 14, 15, 16, 15, 17, 15, 14, 11, + 9, 8, 8, 9, 9, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 8, +]; + +#[rustfmt::skip] +const T24C: [u16; 256] = [ + 15, 13, 46, 80, 146, 262, 248, 434, 426, 669, 653, 649, 621, 517, 1032, 88, + 14, 12, 21, 38, 71, 130, 122, 216, 209, 198, 327, 345, 319, 297, 279, 42, + 47, 22, 41, 74, 68, 128, 120, 221, 207, 194, 182, 340, 315, 295, 541, 18, + 81, 39, 75, 70, 134, 125, 116, 220, 204, 190, 178, 325, 311, 293, 271, 16, + 147, 72, 69, 135, 127, 118, 112, 210, 200, 188, 352, 323, 306, 285, 540, 14, + 263, 66, 129, 126, 119, 114, 214, 202, 192, 180, 341, 317, 301, 281, 262, 12, + 249, 123, 121, 117, 113, 215, 206, 195, 185, 347, 330, 308, 291, 272, 520, 10, + 435, 115, 111, 109, 211, 203, 196, 187, 353, 332, 313, 298, 283, 531, 381, 17, + 427, 212, 208, 205, 201, 193, 186, 177, 169, 320, 303, 286, 268, 514, 377, 16, + 335, 199, 197, 191, 189, 181, 174, 333, 321, 305, 289, 275, 521, 379, 371, 11, + 668, 184, 183, 179, 175, 344, 331, 314, 304, 290, 277, 530, 383, 373, 366, 10, + 652, 346, 171, 168, 164, 318, 309, 299, 287, 276, 263, 513, 375, 368, 362, 6, + 648, 322, 316, 312, 307, 302, 292, 284, 269, 261, 512, 376, 370, 364, 359, 4, + 620, 300, 296, 294, 288, 282, 273, 266, 515, 380, 374, 369, 365, 361, 357, 2, + 1033, 280, 278, 274, 267, 264, 259, 382, 378, 372, 367, 363, 360, 358, 356, 0, + 43, 20, 19, 17, 15, 13, 11, 9, 7, 6, 4, 7, 5, 3, 1, 3, +]; +#[rustfmt::skip] +const T24L: [u8; 256] = [ + 4, 4, 6, 7, 8, 9, 9, 10, 10, 11, 11, 11, 11, 11, 12, 9, + 4, 4, 5, 6, 7, 8, 8, 9, 9, 9, 10, 10, 10, 10, 10, 8, + 6, 5, 6, 7, 7, 8, 8, 9, 9, 9, 9, 10, 10, 10, 11, 7, + 7, 6, 7, 7, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 7, + 8, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 10, 11, 7, + 9, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 10, 7, + 9, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 7, + 10, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 8, + 10, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 8, + 10, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 8, + 11, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 8, + 11, 10, 9, 9, 9, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 8, + 11, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 8, + 11, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 8, + 12, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 8, + 8, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 4, +]; + +/// The pair tables by `table_select` (0-31); `None` for table 0 (all zeros) and the unused 4 and 14. +pub(crate) fn pair_table(select: usize) -> Option { + let (size, codes, lens): (usize, &'static [u16], &'static [u8]) = match select { + 1 => (2, &T1C, &T1L), + 2 => (3, &T2C, &T2L), + 3 => (3, &T3C, &T3L), + 5 => (4, &T5C, &T5L), + 6 => (4, &T6C, &T6L), + 7 => (6, &T7C, &T7L), + 8 => (6, &T8C, &T8L), + 9 => (6, &T9C, &T9L), + 10 => (8, &T10C, &T10L), + 11 => (8, &T11C, &T11L), + 12 => (8, &T12C, &T12L), + 13 => (16, &T13C, &T13L), + 15 => (16, &T15C, &T15L), + 16..=23 => (16, &T16C, &T16L), + 24..=31 => (16, &T24C, &T24L), + _ => return None, + }; + Some(PairTable { size, codes, lens }) +} + +/// `linbits` of each `table_select`. +pub(crate) const LINBITS: [u32; 32] = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 6, 8, 10, 13, 4, 5, 6, 7, 8, 9, 11, 13]; + +/// Count1 table A (`count1table_select` 0): codes and lengths of the quadruple `8v + 4w + 2x + y`. +pub(crate) const QUAD_A_CODES: [u8; 16] = [1, 5, 4, 5, 6, 5, 4, 4, 7, 3, 6, 0, 7, 2, 3, 1]; +pub(crate) const QUAD_A_LENS: [u8; 16] = [1, 4, 4, 5, 4, 6, 5, 6, 4, 5, 5, 6, 5, 6, 6, 6]; +/// Count1 table B: four bits, the inverted quadruple. +pub(crate) const QUAD_B_CODES: [u8; 16] = [15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]; diff --git a/crates/audio/src/mp3/mod.rs b/crates/audio/src/mp3/mod.rs new file mode 100644 index 00000000..dd7591d6 --- /dev/null +++ b/crates/audio/src/mp3/mod.rs @@ -0,0 +1,502 @@ +//! Clean-room MPEG-1 / MPEG-2 Audio Layer III (MP3) encoder, written from ISO/IEC 11172-3:1993 +//! (MPEG-1 Audio: §2.4.1 bit stream syntax, §2.4.2 semantics, Annex B tables, Annex C encoding +//! process) and ISO/IEC 13818-3:1998 (the lower sampling frequencies). +//! +//! - Constant bit rate: 32-320 kbit/s at 32 / 44.1 / 48 kHz (MPEG-1, two granules per frame) and +//! 8-160 kbit/s at 16 / 22.05 / 24 kHz (MPEG-2 LSF, one granule per frame). +//! - Mono, stereo, or joint stereo with frame-by-frame mid/side coding. +//! - The polyphase analysis filterbank of Annex C.1.3, the normal-block MDCT and alias reduction; +//! long blocks only (no block switching yet). +//! - A simple psychoacoustic model (spread band energies against the threshold in quiet), the +//! inner / outer quantization loops with scalefactors, Huffman table and region selection, and +//! the bit reservoir (`main_data_begin`). +//! - Output: a plain stream of frames, no ID3 or Xing / Info tag, no CRC. The frames written one +//! after another are a playable `.mp3` file. +//! +//! No decoder here: [`mpa`](crate::mpa) decodes Layers I and II, and files with Layer III are read +//! through `symphonia` (feature `decode`). +//! +//! ``` +//! use craft_audio::mp3::{Config, DELAY, Encoder}; +//! +//! let left: Vec = (0..48_000).map(|i| (i as f32 * 0.05).sin() * 0.5).collect(); +//! let right = left.clone(); +//! let mut enc = Encoder::new(Config::cbr(48_000, 2, 192))?; +//! let mut mp3 = enc.encode(&[&left, &right]); // whole frames as they complete +//! mp3.extend(enc.flush()); // pads with silence to the last input sample +//! assert_eq!(mp3.len() % 576, 0); // 192 kbit/s at 48 kHz: 576 bytes per frame +//! assert_eq!(DELAY, 1057); // a decoder's output lags the input by this many samples +//! # Ok::<(), craft_audio::mp3::Error>(()) +//! ``` +//! +//! From FilmCraft's `filmcraft-mp3enc` (storytold/filmcraft@4bec22b). Public items keep their +//! names, so FilmCraft can re-export this module. + +// Index loops follow the specification's pseudo-code (FilmCraft allows this lint workspace-wide). +#![allow(clippy::needless_range_loop)] + +mod filterbank; +mod huffman; +mod huffman_tables; +mod quantize; +mod tables; + +use std::collections::VecDeque; + +use craft_media::bitstream::BitWriter; + +use filterbank::{Analysis, GRANULE}; +use quantize::{Granule, Model}; + +/// Settings an [`Encoder`] refuses. +#[non_exhaustive] +#[derive(Debug, thiserror::Error, PartialEq, Eq)] +pub enum Error { + /// A sample rate outside [`SAMPLE_RATES`]. + #[error("MP3: unsupported sample rate {0} Hz (32, 44.1, 48 kHz, or 16, 22.05, 24 kHz)")] + SampleRate(u32), + /// A bit rate that is not one of [`bitrates`] for the sample rate. + #[error("MP3: unsupported bit rate {kbps} kbit/s at {rate} Hz")] + Bitrate { + /// The bit rate asked for, kbit/s. + kbps: u32, + /// The sample rate, Hz. + rate: u32, + }, + /// More than two channels, or none. + #[error("MP3: unsupported channel count {0} (1 or 2)")] + Channels(usize), +} + +/// Result with this module's [`Error`]. +pub type Result = std::result::Result; + +/// How two channels are coded. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub enum StereoMode { + /// Left and right coded independently (`mode` 00). + Stereo, + /// Joint stereo (`mode` 01): each frame chooses left / right or mid / side. + #[default] + JointStereo, +} + +/// Encoder settings. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub struct Config { + /// One of [`SAMPLE_RATES`]. + pub sample_rate: u32, + /// 1 or 2. + pub channels: usize, + /// One of [`bitrates`] for the sample rate, kbit/s. + pub bitrate_kbps: u32, + /// How two channels are coded (ignored for mono). + pub stereo: StereoMode, +} + +impl Config { + /// `channels` at `sample_rate`, constant `bitrate_kbps`, joint stereo. + pub fn cbr(sample_rate: u32, channels: usize, bitrate_kbps: u32) -> Config { + Config { sample_rate, channels, bitrate_kbps, stereo: StereoMode::JointStereo } + } +} + +/// Whether the sample rate needs MPEG-2's lower sampling frequencies. +fn is_lsf(rate: u32) -> Option { + if tables::RATES_MPEG1.contains(&rate) { + Some(false) + } else if tables::RATES_LSF.contains(&rate) { + Some(true) + } else { + None + } +} + +/// The sample rates the encoder accepts. +pub const SAMPLE_RATES: [u32; 6] = [32_000, 44_100, 48_000, 16_000, 22_050, 24_000]; + +/// The bit rates (kbit/s) a stream at `sample_rate` can use; empty for an unsupported rate. +pub fn bitrates(sample_rate: u32) -> &'static [u32] { + match is_lsf(sample_rate) { + Some(false) => &tables::BITRATES_MPEG1, + Some(true) => &tables::BITRATES_LSF, + None => &[], + } +} + +/// The highest bit rate at `sample_rate` not above `kbps` (the lowest one when `kbps` is below +/// all of them); `None` for an unsupported rate. +pub fn nearest_bitrate(sample_rate: u32, kbps: u32) -> Option { + let rates = bitrates(sample_rate); + rates.iter().rev().find(|&&r| r <= kbps).or(rates.first()).copied() +} + +/// The supported sample rate an encoder should use for audio at `rate`: `rate` itself when +/// supported, else the closest supported one at or above it (48 kHz above 48 kHz). +pub fn nearest_sample_rate(rate: u32) -> u32 { + let mut sorted = SAMPLE_RATES; + sorted.sort_unstable(); + sorted.iter().copied().find(|&r| r >= rate).unwrap_or(48_000) +} + +/// Energy of the frequency lines of a full-scale sine (amplitude 1) over one granule. +const FULL_SCALE_TONE_ENERGY: f32 = 1.0; +/// Playback level of a full-scale sine assumed by the threshold in quiet, dB SPL. +const FULL_SCALE_SPL: f64 = 96.0; +/// Samples by which a standard decoder's output lags the input: the polyphase analysis and +/// synthesis filters (481) and the MDCT overlap of one granule (576). +pub const DELAY: u32 = 1057; + +/// A Layer III encoder. Feed planar `f32` samples (full scale ±1) to [`Encoder::encode`]; it +/// returns whole frames as they are completed. [`Encoder::flush`] pads the end with silence and +/// returns the rest. +pub struct Encoder { + cfg: Config, + lsf: bool, + sr_index: u32, + br_index: u32, + granules: usize, + analysis: Vec, + model: Model, + /// Buffered input per channel, and how much of it is consumed. + input: Vec>, + consumed: usize, + /// Frames whose main-data slot is not yet filled: header + side information, slot bytes. + pending: VecDeque<(Vec, usize)>, + /// Main data from the first pending frame's slot on. + main: Vec, + /// Slot bytes of the pending frames. + slot_total: usize, + /// Largest `main_data_begin`, bytes. + max_reservoir: usize, + /// Remainder accumulator for the padding slot (44.1 / 22.05 kHz). + pad_acc: u64, + samples_in: u64, + frames: u64, +} + +impl Encoder { + /// An encoder for `cfg`, or the setting it cannot carry. + pub fn new(cfg: Config) -> Result { + let lsf = is_lsf(cfg.sample_rate).ok_or(Error::SampleRate(cfg.sample_rate))?; + if !(1..=2).contains(&cfg.channels) { + return Err(Error::Channels(cfg.channels)); + } + let rates = if lsf { &tables::RATES_LSF } else { &tables::RATES_MPEG1 }; + let sr_index = rates.iter().position(|&r| r == cfg.sample_rate).ok_or(Error::SampleRate(cfg.sample_rate))?; + let br_index = bitrates(cfg.sample_rate) + .iter() + .position(|&b| b == cfg.bitrate_kbps) + .ok_or(Error::Bitrate { kbps: cfg.bitrate_kbps, rate: cfg.sample_rate })?; + let bands = if lsf { tables::SFB_LSF[sr_index] } else { tables::SFB_MPEG1[sr_index] }; + let line_hz = cfg.sample_rate as f64 / 2.0 / GRANULE as f64; + let mut ath = [0f32; 22]; + for (b, a) in ath.iter_mut().enumerate() { + let (lo, hi) = (usize::from(bands[b]), usize::from(bands[b + 1])); + let db = (lo..hi).map(|i| tables::ath_db((i as f64 + 0.5) * line_hz)).fold(f64::MAX, f64::min); + *a = (f64::from(FULL_SCALE_TONE_ENERGY) * 10f64.powf((db - FULL_SCALE_SPL) / 10.0)) as f32; + } + let cutoff_hz = cutoff_hz(cfg.bitrate_kbps, cfg.channels).min(cfg.sample_rate as f64 * 0.5); + let cutoff = ((cutoff_hz / line_hz).ceil() as usize).min(GRANULE); + let frame_bits = u64::from(cfg.bitrate_kbps) * 1000 * if lsf { 576 } else { 1152 } / u64::from(cfg.sample_rate); + // reservoir and frame must fit the decoder's 7680-bit input buffer + let max_reservoir = (7680u64.saturating_sub(frame_bits) / 8).min(if lsf { 255 } else { 511 }) as usize; + Ok(Encoder { + cfg, + lsf, + sr_index: sr_index as u32, + br_index: br_index as u32 + 1, + granules: if lsf { 1 } else { 2 }, + analysis: vec![Analysis::new(); cfg.channels], + model: Model::new(lsf, bands, ath, cutoff), + input: vec![Vec::new(); cfg.channels], + consumed: 0, + pending: VecDeque::new(), + main: Vec::new(), + slot_total: 0, + max_reservoir, + pad_acc: 0, + samples_in: 0, + frames: 0, + }) + } + + /// The settings. + pub fn config(&self) -> &Config { + &self.cfg + } + + /// Samples per channel in a frame: 1152 (MPEG-1) or 576 (MPEG-2 LSF). + pub fn frame_samples(&self) -> usize { + GRANULE * self.granules + } + + /// Samples by which a standard decoder's output lags the input. + pub fn delay(&self) -> u32 { + DELAY + } + + /// Frames written so far. + pub fn frames(&self) -> u64 { + self.frames + } + + /// Encode planar samples (one slice per channel; a missing channel is silence, extra ones + /// are ignored, every channel takes the first channel's length). Returns the completed frames. + pub fn encode(&mut self, planar: &[&[f32]]) -> Vec { + let n = planar.first().map_or(0, |p| p.len()); + for (c, buf) in self.input.iter_mut().enumerate() { + let src = planar.get(c).copied().unwrap_or(&[]); + buf.extend((0..n).map(|i| src.get(i).copied().filter(|v| v.is_finite()).unwrap_or(0.0).clamp(-4.0, 4.0))); + } + self.samples_in = self.samples_in.saturating_add(n as u64); + let mut out = Vec::new(); + self.drain(&mut out); + out + } + + /// Encode the buffered rest padded with silence, plus enough silence that a decoder outputs + /// every input sample, and return the remaining frames. + pub fn flush(&mut self) -> Vec { + let fs = self.frame_samples() as u64; + let needed = if self.samples_in == 0 { self.frames } else { self.samples_in.saturating_add(u64::from(DELAY)).div_ceil(fs) }; + let have = self.input.first().map_or(0, |b| b.len().saturating_sub(self.consumed)); + let want = (needed.saturating_sub(self.frames) as usize).saturating_mul(fs as usize); + if want > have { + for buf in &mut self.input { + buf.resize(self.consumed + want, 0.0); + } + } + let mut out = Vec::new(); + self.drain(&mut out); + // fill the last frames' slots and write them out + self.main.resize(self.slot_total.max(self.main.len()), 0); + self.emit(&mut out); + self.input.iter_mut().for_each(Vec::clear); + self.consumed = 0; + out + } + + /// Encode every whole frame of buffered input. + fn drain(&mut self, out: &mut Vec) { + let fs = self.frame_samples(); + let mut pcm = vec![[0f32; GRANULE]; self.cfg.channels * self.granules]; + while self.input.first().is_some_and(|b| b.len() - self.consumed >= fs) { + for (c, buf) in self.input.iter().enumerate() { + for gr in 0..self.granules { + let start = self.consumed + gr * GRANULE; + if let (Some(dst), Some(src)) = (pcm.get_mut(gr * self.cfg.channels + c), buf.get(start..start + GRANULE)) { + dst.copy_from_slice(src); + } + } + } + self.consumed += fs; + self.encode_frame(&pcm); + self.emit(out); + } + if self.consumed > 0 { + let k = self.consumed; + self.input.iter_mut().for_each(|b| { + b.drain(..k.min(b.len())); + }); + self.consumed = 0; + } + } + + /// Bytes of the next frame (with its padding slot). + fn next_frame_len(&mut self) -> usize { + let per = if self.lsf { 72_000 } else { 144_000 }; + let num = per * u64::from(self.cfg.bitrate_kbps); + let sr = u64::from(self.cfg.sample_rate); + self.pad_acc += num % sr; + let pad = self.pad_acc >= sr; + if pad { + self.pad_acc -= sr; + } + (num / sr) as usize + usize::from(pad) + } + + fn side_info_len(&self) -> usize { + match (self.lsf, self.cfg.channels) { + (false, 1) => 17, + (false, _) => 32, + (true, 1) => 9, + (true, _) => 17, + } + } + + /// Analyse, quantize and assemble one frame (`pcm[gr * channels + ch]`). + fn encode_frame(&mut self, pcm: &[[f32; GRANULE]]) { + let nch = self.cfg.channels; + let ngr = self.granules; + let mut xr = vec![[0f32; GRANULE]; nch * ngr]; + for gr in 0..ngr { + for ch in 0..nch { + let k = gr * nch + ch; + if let (Some(a), Some(p), Some(x)) = (self.analysis.get_mut(ch), pcm.get(k), xr.get_mut(k)) { + a.granule(p, x); + x.iter_mut().skip(self.model.cutoff).for_each(|v| *v = 0.0); + } + } + } + let mut xmin: Vec<[f32; 22]> = xr.iter().map(|x| self.model.thresholds(x)).collect(); + let ms = nch == 2 && self.cfg.stereo == StereoMode::JointStereo && use_mid_side(&xr); + if ms { + for gr in 0..ngr { + let (l, r) = xr.split_at_mut(gr * 2 + 1); + if let (Some(l), Some(r)) = (l.last_mut(), r.first_mut()) { + for (a, b) in l.iter_mut().zip(r.iter_mut()) { + let (m, s) = ((*a + *b) * std::f32::consts::FRAC_1_SQRT_2, (*a - *b) * std::f32::consts::FRAC_1_SQRT_2); + *a = m; + *b = s; + } + } + let t: [f32; 22] = std::array::from_fn(|b| xmin[gr * 2][b].min(xmin[gr * 2 + 1][b])); + xmin[gr * 2] = t; + xmin[gr * 2 + 1] = t; + } + } + // the reservoir: free bytes in the pending frames' slots, stuffed down to the maximum + let free = self.slot_total.saturating_sub(self.main.len()); + if free > self.max_reservoir { + self.main.resize(self.slot_total - self.max_reservoir, 0); + } + let main_data_begin = free.min(self.max_reservoir); + let frame_len = self.next_frame_len(); + let slot = frame_len.saturating_sub(4 + self.side_info_len()); + let budget = 8 * (main_data_begin + slot); + let n = nch * ngr; + let mean = 8 * slot / n; + let max_reservoir_bits = 8 * self.max_reservoir; + let mut used = 0usize; + let mut coded: Vec = Vec::with_capacity(n); + for (k, (x, t)) in xr.iter().zip(&xmin).enumerate() { + let left = n - k; + let remaining = budget.saturating_sub(used); + let surplus = remaining.saturating_sub(mean * left); + let max_bits = (mean + surplus * 3 / 5 / left).min(remaining).min(4095); + let room = max_reservoir_bits.saturating_sub(surplus); + let min_bits = mean.saturating_sub(room / left); + let g = self.model.encode(x, t, max_bits, min_bits); + used += g.bits(); + coded.push(g); + } + let mut md = BitWriter::new(); + for g in &coded { + g.write(&mut md, &self.model); + } + let mut main = md.finish(); + if main.len() > main_data_begin + slot || coded.iter().any(|g| g.bits() > 4095) { + // never expected: fall back to silence rather than write a broken frame + coded = (0..n).map(|_| self.model.encode(&[0.0; GRANULE], &[f32::MAX; 22], 0, 0)).collect(); + main = Vec::new(); + } + let head = self.header_and_side_info(frame_len, ms, main_data_begin, &coded); + self.pending.push_back((head, slot)); + self.slot_total += slot; + self.main.extend_from_slice(&main); + self.frames += 1; + } + + /// Write out the pending frames whose slots are full. + fn emit(&mut self, out: &mut Vec) { + while let Some((_, slot)) = self.pending.front() { + let slot = *slot; + if self.main.len() < slot { + break; + } + if let Some((head, _)) = self.pending.pop_front() { + out.extend_from_slice(&head); + out.extend(self.main.drain(..slot)); + self.slot_total -= slot; + } + } + } + + fn header_and_side_info(&self, frame_len: usize, ms: bool, main_data_begin: usize, coded: &[Granule]) -> Vec { + let nch = self.cfg.channels; + let padding = frame_len > (if self.lsf { 72_000 } else { 144_000 } * self.cfg.bitrate_kbps / self.cfg.sample_rate) as usize; + let mode = match (nch, self.cfg.stereo) { + (1, _) => 3, + (_, StereoMode::Stereo) => 0, + (_, StereoMode::JointStereo) => 1, + }; + let mut w = BitWriter::new(); + w.write_bits(0x7FF, 11); + w.write_bits(if self.lsf { 2 } else { 3 }, 2); // ID: MPEG-1 = 11, MPEG-2 = 10 + w.write_bits(1, 2); // Layer III + w.write_bit(true); // protection_bit: no CRC + w.write_bits(self.br_index, 4); + w.write_bits(self.sr_index, 2); + w.write_bit(padding); + w.write_bit(false); // private_bit + w.write_bits(mode, 2); + w.write_bits(if ms { 2 } else { 0 }, 2); // mode_extension: mid/side, no intensity + w.write_bit(false); // copyright + w.write_bit(true); // original + w.write_bits(0, 2); // emphasis + if self.lsf { + w.write_bits(main_data_begin as u32, 8); + w.write_bits(0, if nch == 1 { 1 } else { 2 }); + } else { + w.write_bits(main_data_begin as u32, 9); + w.write_bits(0, if nch == 1 { 5 } else { 3 }); + w.write_bits(0, 4 * nch as u32); // scfsi + } + for g in coded { + let c = &g.coding; + w.write_bits(g.bits() as u32, 12); + w.write_bits(c.big_values as u32, 9); + w.write_bits(g.global_gain, 8); + w.write_bits(g.scalefac_compress, if self.lsf { 9 } else { 4 }); + w.write_bit(false); // window_switching_flag + for t in c.table_select { + w.write_bits(t as u32, 5); + } + w.write_bits(c.region0_count as u32, 4); + w.write_bits(c.region1_count as u32, 3); + if !self.lsf { + w.write_bit(false); // preflag + } + w.write_bit(g.scalefac_scale); + w.write_bits(c.count1_table as u32, 1); + } + w.finish() + } +} + +/// Mid/side for a frame when the side signal is much weaker than the mid signal. +fn use_mid_side(xr: &[[f32; GRANULE]]) -> bool { + let (mut m, mut s) = (0f64, 0f64); + for pair in xr.as_chunks::<2>().0 { + for (l, r) in pair[0].iter().zip(pair[1].iter()) { + let (a, b) = (f64::from(*l), f64::from(*r)); + m += (a + b) * (a + b); + s += (a - b) * (a - b); + } + } + s < 0.25 * m +} + +/// The low-pass for a bit rate: the band the bits can carry (joint stereo counted as 1.6 +/// channels' worth). +fn cutoff_hz(kbps: u32, channels: usize) -> f64 { + let per = f64::from(kbps) / if channels > 1 { 1.6 } else { 1.0 }; + const STEPS: [(f64, f64); 11] = [ + (140.0, 20_000.0), + (120.0, 19_000.0), + (100.0, 18_000.0), + (80.0, 16_500.0), + (64.0, 15_000.0), + (56.0, 13_500.0), + (48.0, 12_000.0), + (40.0, 10_000.0), + (32.0, 8_000.0), + (24.0, 6_500.0), + (16.0, 5_000.0), + ]; + STEPS.iter().find(|(k, _)| per >= *k).map_or(3_500.0, |s| s.1) +} + +#[cfg(test)] +mod tests; diff --git a/crates/audio/src/mp3/quantize.rs b/crates/audio/src/mp3/quantize.rs new file mode 100644 index 00000000..24d07eeb --- /dev/null +++ b/crates/audio/src/mp3/quantize.rs @@ -0,0 +1,370 @@ +//! A simple psychoacoustic model and the quantization loops (after ISO/IEC 11172-3 Annex C.1.5.4). +//! +//! The model gives every scalefactor band an allowed noise energy: the band energies spread to +//! their neighbours (10 dB per band upwards, 25 dB per band downwards) lowered by a fixed +//! signal-to-mask offset, but never below the absolute threshold of hearing. The inner loop picks +//! the finest global gain whose Huffman code fits the granule's bit budget; the outer loop +//! amplifies (raises the scalefactors of) the bands whose quantization noise exceeds the allowed +//! noise and keeps the result with the lowest worst noise-to-mask ratio. When every band is under +//! its threshold and the bit reservoir can take the savings, the global gain is coarsened again +//! as far as the thresholds allow. + +use craft_media::bitstream::BitWriter; + +use super::filterbank::GRANULE; +use super::huffman::{self, Coding}; +use super::tables::{LSF_PARTITIONS, MAX_QUANT, SLEN, pow43}; + +/// Outer-loop iterations per granule (each raises at least one scalefactor). +const MAX_OUTER: usize = 24; +/// Signal-to-mask offset of the spread band energies, dB. +const MASK_OFFSET_DB: f32 = 12.0; + +/// The fixed per-stream parameters of the model and the loops. +pub(crate) struct Model { + pub lsf: bool, + pub bands: [u16; 23], + /// Allowed noise energy from the threshold in quiet, per band. + pub ath: [f32; 22], + /// Lines at and above this are not coded (the low-pass). + pub cutoff: usize, + /// `spread[b][j]`: how much of band `j`'s energy masks band `b`. + spread: [[f32; 22]; 22], +} + +/// One coded granule of one channel. +#[derive(Clone)] +pub(crate) struct Granule { + pub ix: [u16; GRANULE], + pub negative: [bool; GRANULE], + pub global_gain: u32, + pub scalefac: [u8; 22], + pub scalefac_scale: bool, + pub scalefac_compress: u32, + /// Bits per scalefactor of each band group (MPEG-1: bands 0-10, 11-20; MPEG-2: four partitions). + pub slen: [u32; 4], + pub part2_bits: usize, + pub coding: Coding, +} + +impl Granule { + fn silent() -> Granule { + Granule { + ix: [0; GRANULE], + negative: [false; GRANULE], + global_gain: 210, + scalefac: [0; 22], + scalefac_scale: false, + scalefac_compress: 0, + slen: [0; 4], + part2_bits: 0, + coding: Coding::default(), + } + } + + /// `part2_3_length`: scalefactor and Huffman bits. + pub fn bits(&self) -> usize { + self.part2_bits + self.coding.bits + } + + /// Write the scalefactors (no scfsi reuse) and the Huffman code. + pub fn write(&self, w: &mut BitWriter, model: &Model) { + let groups: [(usize, usize, u32); 4] = if model.lsf { + let mut start = 0; + let mut g = [(0, 0, 0); 4]; + for (p, (&n, gi)) in LSF_PARTITIONS.iter().zip(g.iter_mut()).enumerate() { + *gi = (start, start + n, self.slen[p]); + start += n; + } + g + } else { + [(0, 11, self.slen[0]), (11, 21, self.slen[1]), (21, 21, 0), (21, 21, 0)] + }; + for (a, b, bits) in groups { + for &sf in self.scalefac.get(a..b).unwrap_or(&[]) { + w.write_bits(u32::from(sf), bits); + } + } + huffman::write(w, &self.ix, &self.negative, &self.coding, &model.bands); + } +} + +/// Bits needed for `v` (0 for 0). +fn bit_width(v: u32) -> u32 { + 32 - v.leading_zeros() +} + +impl Model { + pub fn new(lsf: bool, bands: [u16; 23], ath: [f32; 22], cutoff: usize) -> Model { + let mut spread = [[0f32; 22]; 22]; + for (b, row) in spread.iter_mut().enumerate() { + for (j, s) in row.iter_mut().enumerate() { + let db = if b >= j { (b - j) as f32 * 10.0 } else { (j - b) as f32 * 25.0 }; + *s = 10f32.powf(-(db + MASK_OFFSET_DB) / 10.0); + } + } + Model { lsf, bands, ath, cutoff: cutoff.min(GRANULE), spread } + } + + fn band(&self, b: usize) -> std::ops::Range { + let lo = self.bands.get(b).map_or(GRANULE, |&v| usize::from(v)); + let hi = self.bands.get(b + 1).map_or(GRANULE, |&v| usize::from(v)); + lo..hi.max(lo) + } + + /// The largest scalefactor of band `b` the side information can carry. + fn sf_limit(&self, b: usize) -> u8 { + // MPEG-1 slen1 / MPEG-2 slen1-2 reach 4 bits, slen2 / slen3-4 3 bits; band 21 has none + match b { + 0..=10 => 15, + 11..=20 => 7, + _ => 0, + } + } + + /// Allowed noise energy per band for the lines `xr` (zero above the cutoff). + pub fn thresholds(&self, xr: &[f32; GRANULE]) -> [f32; 22] { + let mut en = [0f32; 22]; + for (b, e) in en.iter_mut().enumerate() { + *e = xr.get(self.band(b)).unwrap_or(&[]).iter().map(|v| v * v).sum(); + } + let mut xmin = [0f32; 22]; + for (b, x) in xmin.iter_mut().enumerate() { + let r = self.band(b); + if r.start >= self.cutoff { + *x = f32::MAX; + continue; + } + let mask: f32 = self.spread[b].iter().zip(&en).map(|(s, e)| s * e).sum(); + *x = mask.max(self.ath[b]); + } + xmin + } + + /// Scalefactor bits and `scalefac_compress` (and the slen per group) for `sf`. + fn part2(&self, sf: &[u8; 22]) -> (usize, u32, [u32; 4]) { + let max_in = |a: usize, b: usize| bit_width(u32::from(sf.get(a..b).and_then(|s| s.iter().copied().max()).unwrap_or(0))); + if self.lsf { + let mut slen = [0u32; 4]; + let mut start = 0; + let mut bits = 0usize; + for (p, &n) in LSF_PARTITIONS.iter().enumerate() { + slen[p] = max_in(start, start + n); + bits += n * slen[p] as usize; + start += n; + } + // scalefac_compress < 400: slen1, slen2 ≤ 4; slen3, slen4 ≤ 3 (13818-3 §2.4.3.2) + let sfc = ((slen[0] * 5 + slen[1]) << 4) + (slen[2] << 2) + slen[3]; + (bits, sfc, slen) + } else { + let (n1, n2) = (max_in(0, 11), max_in(11, 21)); + let mut best = (usize::MAX, 15u32, [4, 3, 0, 0]); + for (i, &(s1, s2)) in SLEN.iter().enumerate() { + let bits = 11 * s1 as usize + 10 * s2 as usize; + if s1 >= n1 && s2 >= n2 && bits < best.0 { + best = (bits, i as u32, [s1, s2, 0, 0]); + } + } + best + } + } + + /// Code one granule's lines `xr` with at most `max_bits` and, when the thresholds `xmin` are + /// met, preferably no fewer than `min_bits`. + pub fn encode(&self, xr: &[f32; GRANULE], xmin: &[f32; 22], max_bits: usize, min_bits: usize) -> Granule { + let mut q = Quantizer { model: self, xr, xr34: [0.0; GRANULE] }; + let mut peak = 0f32; + for (o, v) in q.xr34.iter_mut().zip(xr.iter()).take(self.cutoff) { + let a = v.abs(); + *o = a.sqrt() * a.sqrt().sqrt(); + peak = peak.max(a); + } + let mut silent = Granule::silent(); + if peak <= 1e-9 || max_bits == 0 { + return silent; + } + silent.negative.iter_mut().zip(xr.iter()).for_each(|(n, v)| *n = *v < 0.0); + let mut sf = [0u8; 22]; + let mut scale = false; + let mut hint: Option = None; + let mut best: Option<(f32, Granule)> = None; + for _ in 0..MAX_OUTER { + let (p2, sfc, slen) = self.part2(&sf); + if p2 >= max_bits { + break; + } + let Some((gg, ix, coding)) = q.inner(&sf, scale, max_bits - p2, hint) else { break }; + hint = Some(gg); + let nmr = q.nmr(&ix, gg, &sf, scale, xmin); + let worst = nmr.iter().copied().fold(0f32, f32::max); + if best.as_ref().is_none_or(|b| worst < b.0) { + let mut g = silent.clone(); + g.ix = ix; + g.global_gain = gg; + g.scalefac = sf; + g.scalefac_scale = scale; + g.scalefac_compress = sfc; + g.slen = slen; + g.part2_bits = p2; + g.coding = coding; + best = Some((worst, g)); + } + // amplify the bands over their threshold + let over: Vec = (0..21).filter(|&b| nmr[b] > 1.0).collect(); + if over.is_empty() || over.len() == 21 { + break; + } + for &b in &over { + sf[b] = sf[b].saturating_add(1); + } + if (0..21).any(|b| sf[b] > self.sf_limit(b)) { + if scale { + break; + } + scale = true; + sf.iter_mut().for_each(|s| *s = s.div_ceil(2)); + } + } + let Some((worst, mut g)) = best else { return silent }; + if worst <= 1.0 && g.bits() > min_bits { + q.coarsen(&mut g, xmin, min_bits); + } + g + } +} + +struct Quantizer<'a> { + model: &'a Model, + xr: &'a [f32; GRANULE], + /// `|xr|^(3/4)`. + xr34: [f32; GRANULE], +} + +impl Quantizer<'_> { + /// Quantize with global gain `gg` and the scalefactors; `None` when a line exceeds the + /// largest escape value. + fn quantize(&self, gg: u32, sf: &[u8; 22], scale: bool) -> Option<[u16; GRANULE]> { + let mut ix = [0u16; GRANULE]; + let sfm = if scale { 1.0 } else { 0.5 }; + for b in 0..22 { + let r = self.model.band(b); + let e = -0.1875 * (gg as f32 - 210.0) + 0.75 * sfm * f32::from(sf[b]); + let f = 2f32.powf(e); + for i in r.start..r.end.min(self.model.cutoff) { + let v = self.xr34[i] * f + 0.4054; + if v.is_nan() || v >= (MAX_QUANT + 1) as f32 { + return None; + } + ix[i] = v as u16; + } + } + Some(ix) + } + + /// Noise-to-mask ratio per band. + fn nmr(&self, ix: &[u16; GRANULE], gg: u32, sf: &[u8; 22], scale: bool, xmin: &[f32; 22]) -> [f32; 22] { + let p43 = pow43(); + let sfm = if scale { 1.0 } else { 0.5 }; + let mut out = [0f32; 22]; + for (b, o) in out.iter_mut().enumerate() { + let r = self.model.band(b); + let step = 2f32.powf(0.25 * (gg as f32 - 210.0) - sfm * f32::from(sf[b])); + let mut noise = 0f32; + for i in r.start..r.end.min(self.model.cutoff) { + let rec = p43.get(usize::from(ix[i])).copied().unwrap_or(0.0) * step; + let d = self.xr[i].abs() - rec; + noise += d * d; + } + *o = if xmin[b] > 0.0 { + noise / xmin[b] + } else if noise > 0.0 { + f32::MAX + } else { + 0.0 + }; + } + out + } + + /// The code of `gg` if it fits `budget` bits. + fn fits(&self, gg: u32, sf: &[u8; 22], scale: bool, budget: usize) -> Option<([u16; GRANULE], Coding)> { + let ix = self.quantize(gg, sf, scale)?; + let c = huffman::choose(&ix, &self.model.bands); + (c.bits <= budget).then_some((ix, c)) + } + + /// The inner loop: the smallest global gain whose code fits `budget` bits. + fn inner(&self, sf: &[u8; 22], scale: bool, budget: usize, hint: Option) -> Option<(u32, [u16; GRANULE], Coding)> { + match hint { + None => { + let mut hi = 255u32; + let mut found = self.fits(hi, sf, scale, budget)?; + let mut lo = 0u32; + while lo < hi { + let mid = (lo + hi) / 2; + match self.fits(mid, sf, scale, budget) { + Some(r) => { + hi = mid; + found = r; + } + None => lo = mid + 1, + } + } + Some((hi, found.0, found.1)) + } + Some(start) => { + let mut g = start.min(255); + if let Some(mut found) = self.fits(g, sf, scale, budget) { + while g > 0 { + match self.fits(g - 1, sf, scale, budget) { + Some(r) => { + g -= 1; + found = r; + } + None => break, + } + } + Some((g, found.0, found.1)) + } else { + while g < 255 { + g += 1; + if let Some(r) = self.fits(g, sf, scale, budget) { + return Some((g, r.0, r.1)); + } + } + None + } + } + } + } + + /// Raise the global gain of `g` as far as every band stays under its threshold and the code + /// keeps at least `min_bits`. + fn coarsen(&self, g: &mut Granule, xmin: &[f32; 22], min_bits: usize) { + let ok = |gg: u32| -> Option<([u16; GRANULE], Coding)> { + let ix = self.quantize(gg, &g.scalefac, g.scalefac_scale)?; + let c = huffman::choose(&ix, &self.model.bands); + let nmr = self.nmr(&ix, gg, &g.scalefac, g.scalefac_scale, xmin); + (nmr.iter().all(|&v| v <= 1.0) && g.part2_bits + c.bits >= min_bits).then_some((ix, c)) + }; + let (mut lo, mut hi) = (g.global_gain, 255u32); + let mut found = None; + while lo < hi { + let mid = (lo + hi).div_ceil(2); + match ok(mid) { + Some(r) => { + lo = mid; + found = Some(r); + } + None => hi = mid - 1, + } + } + if let Some((ix, c)) = found + && lo > g.global_gain + { + g.global_gain = lo; + g.ix = ix; + g.coding = c; + } + } +} diff --git a/crates/audio/src/mp3/tables.rs b/crates/audio/src/mp3/tables.rs new file mode 100644 index 00000000..054df513 --- /dev/null +++ b/crates/audio/src/mp3/tables.rs @@ -0,0 +1,125 @@ +//! Tables from ISO/IEC 11172-3 (MPEG-1) and ISO/IEC 13818-3 (MPEG-2 lower sampling frequencies). + +use std::sync::OnceLock; + +/// Layer III bit rates in kbit/s by `bitrate_index` 1-14 (§2.4.2.3); index 0 (free format) and +/// 15 (forbidden) are never written. +pub(crate) const BITRATES_MPEG1: [u32; 14] = [32, 40, 48, 56, 64, 80, 96, 112, 128, 160, 192, 224, 256, 320]; +/// Layer III bit rates of the lower sampling frequencies (ISO/IEC 13818-3 §2.4.2.3). +pub(crate) const BITRATES_LSF: [u32; 14] = [8, 16, 24, 32, 40, 48, 56, 64, 80, 96, 112, 128, 144, 160]; + +/// Sampling frequencies by `sampling_frequency` (0-2) of MPEG-1 and of MPEG-2 LSF. +pub(crate) const RATES_MPEG1: [u32; 3] = [44_100, 48_000, 32_000]; +pub(crate) const RATES_LSF: [u32; 3] = [22_050, 24_000, 16_000]; + +/// Long-block scalefactor band boundaries (Table 3-B.8 of 11172-3; Table B.2 of 13818-3), by +/// `sampling_frequency`: 22 bands over the 576 frequency lines of a granule. +pub(crate) const SFB_MPEG1: [[u16; 23]; 3] = [ + [0, 4, 8, 12, 16, 20, 24, 30, 36, 44, 52, 62, 74, 90, 110, 134, 162, 196, 238, 288, 342, 418, 576], + [0, 4, 8, 12, 16, 20, 24, 30, 36, 42, 50, 60, 72, 88, 106, 128, 156, 190, 230, 276, 330, 384, 576], + [0, 4, 8, 12, 16, 20, 24, 30, 36, 44, 54, 66, 82, 102, 126, 156, 194, 240, 296, 364, 448, 550, 576], +]; +pub(crate) const SFB_LSF: [[u16; 23]; 3] = [ + [0, 6, 12, 18, 24, 30, 36, 44, 54, 66, 80, 96, 116, 140, 168, 200, 238, 284, 336, 396, 464, 522, 576], + [0, 6, 12, 18, 24, 30, 36, 44, 54, 66, 80, 96, 114, 136, 162, 194, 232, 278, 332, 394, 464, 540, 576], + [0, 6, 12, 18, 24, 30, 36, 44, 54, 66, 80, 96, 116, 140, 168, 200, 238, 284, 336, 396, 464, 522, 576], +]; + +/// `slen1` (bands 0-10) and `slen2` (bands 11-20) by MPEG-1 `scalefac_compress` (§2.4.2.7). +pub(crate) const SLEN: [(u32, u32); 16] = + [(0, 0), (0, 1), (0, 2), (0, 3), (3, 0), (1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3), (3, 1), (3, 2), (3, 3), (4, 2), (4, 3)]; + +/// Scalefactor bands per partition for long blocks without intensity stereo when MPEG-2 +/// `scalefac_compress` < 400 (13818-3 §2.4.3.2, `nr_of_sfb_block[0][0]`). +pub(crate) const LSF_PARTITIONS: [usize; 4] = [6, 5, 5, 5]; + +/// Coefficients `D[0..=256]` of the synthesis window (Table 3-B.3), in units of 2⁻¹⁶: the table's +/// nine-decimal values are exact multiples of 1/65536. The analysis window of Table 3-C.1 is the +/// same prototype divided by 32: `C[i] = D[i] / 32`. +const WINDOW_HALF: [i32; 257] = [ + 0, -1, -1, -1, -1, -1, -1, -2, -2, -2, -2, -3, -3, -4, -4, -5, // + -5, -6, -7, -7, -8, -9, -10, -11, -13, -14, -16, -17, -19, -21, -24, -26, // + -29, -31, -35, -38, -41, -45, -49, -53, -58, -63, -68, -73, -79, -85, -91, -97, // + -104, -111, -117, -125, -132, -139, -147, -154, -161, -169, -176, -183, -190, -196, -202, -208, // + 213, 218, 222, 225, 227, 228, 228, 227, 224, 221, 215, 208, 200, 189, 177, 163, // + 146, 127, 106, 83, 57, 29, -2, -36, -72, -111, -153, -197, -244, -294, -347, -401, // + -459, -519, -581, -645, -711, -779, -848, -919, -991, -1064, -1137, -1210, -1283, -1356, -1428, -1498, // + -1567, -1634, -1698, -1759, -1817, -1870, -1919, -1962, -2001, -2032, -2057, -2075, -2085, -2087, -2080, -2063, // + 2037, 2000, 1952, 1893, 1822, 1739, 1644, 1535, 1414, 1280, 1131, 970, 794, 605, 402, 185, // + -45, -288, -545, -814, -1095, -1388, -1692, -2006, -2330, -2663, -3004, -3351, -3705, -4063, -4425, -4788, // + -5153, -5517, -5879, -6237, -6589, -6935, -7271, -7597, -7910, -8209, -8491, -8755, -8998, -9219, -9416, -9585, // + -9727, -9838, -9916, -9959, -9966, -9935, -9863, -9750, -9592, -9389, -9139, -8840, -8492, -8092, -7640, -7134, // + 6574, 5959, 5288, 4561, 3776, 2935, 2037, 1082, 70, -998, -2122, -3300, -4533, -5818, -7154, -8540, // + -9975, -11455, -12980, -14548, -16155, -17799, -19478, -21189, -22929, -24694, -26482, -28289, -30112, -31947, -33791, -35640, // + -37489, -39336, -41176, -43006, -44821, -46617, -48390, -50137, -51853, -53534, -55178, -56778, -58333, -59838, -61289, -62684, // + -64019, -65290, -66494, -67629, -68692, -69679, -70590, -71420, -72169, -72835, -73415, -73908, -74313, -74630, -74856, -74992, // + 75038, +]; + +/// The 512 synthesis window coefficients `D[i]` (Table 3-B.3). The table stores the symmetric +/// prototype `h[i] = h[512 - i]` with every other block of 64 negated, so for `i > 256` +/// `D[i] = -D[512 - i]`, except at multiples of 64 where both halves share the block parity. +pub(crate) fn synthesis_window() -> &'static [f32; 512] { + static W: OnceLock<[f32; 512]> = OnceLock::new(); + W.get_or_init(|| { + let mut w = [0f32; 512]; + for (i, d) in w.iter_mut().enumerate() { + let (j, sign) = if i <= 256 { (i, 1.0) } else { (512 - i, if i % 64 == 0 { 1.0 } else { -1.0 }) }; + *d = sign * WINDOW_HALF.get(j).copied().unwrap_or(0) as f32 / 65_536.0; + } + w + }) +} + +/// The analysis window `C[i]` (Table 3-C.1): `D[i] / 32`. +pub(crate) fn analysis_window() -> &'static [f32; 512] { + static C: OnceLock<[f32; 512]> = OnceLock::new(); + C.get_or_init(|| { + let mut c = [0f32; 512]; + c.iter_mut().zip(synthesis_window()).for_each(|(c, d)| *c = d / 32.0); + c + }) +} + +/// The analysis matrixing coefficients `M[i][k] = cos((2i + 1)(k - 16)π / 64)` (Annex C.1.3), +/// `i` in `0..32`, `k` in `0..64`. +pub(crate) fn analysis_matrix() -> &'static [[f32; 64]; 32] { + static M: OnceLock<[[f32; 64]; 32]> = OnceLock::new(); + M.get_or_init(|| { + let mut m = [[0f32; 64]; 32]; + for (i, row) in m.iter_mut().enumerate() { + for (k, c) in row.iter_mut().enumerate() { + *c = ((2 * i + 1) as f64 * (k as f64 - 16.0) * std::f64::consts::PI / 64.0).cos() as f32; + } + } + m + }) +} + +/// Alias-reduction butterfly coefficients `(cs[i], ca[i])` from `c[i]` (Table 3-B.9). +pub(crate) fn alias_coefficients() -> &'static [(f32, f32); 8] { + static A: OnceLock<[(f32, f32); 8]> = OnceLock::new(); + A.get_or_init(|| { + const C: [f64; 8] = [-0.6, -0.535, -0.33, -0.185, -0.095, -0.041, -0.0142, -0.0037]; + let mut a = [(0f32, 0f32); 8]; + for (o, c) in a.iter_mut().zip(C) { + let n = (1.0 + c * c).sqrt(); + *o = ((1.0 / n) as f32, (c / n) as f32); + } + a + }) +} + +/// `|i|^(4/3)` for every quantized magnitude a big-values escape can carry (`15 + 2^13 - 1`). +pub(crate) const MAX_QUANT: usize = 8206; +pub(crate) fn pow43() -> &'static [f32] { + static P: OnceLock> = OnceLock::new(); + P.get_or_init(|| (0..=MAX_QUANT).map(|i| (i as f64).powf(4.0 / 3.0) as f32).collect()) +} + +/// The absolute threshold of hearing in dB SPL at `f` Hz (Terhardt's approximation, as in the +/// psychoacoustic models of 11172-3 Annex D). +pub(crate) fn ath_db(f: f64) -> f64 { + let k = (f.max(20.0)) / 1000.0; + 3.64 * k.powf(-0.8) - 6.5 * (-0.6 * (k - 3.3) * (k - 3.3)).exp() + 0.001 * k.powi(4) +} diff --git a/crates/audio/src/mp3/tests.rs b/crates/audio/src/mp3/tests.rs new file mode 100644 index 00000000..a4160f0a --- /dev/null +++ b/crates/audio/src/mp3/tests.rs @@ -0,0 +1,490 @@ +use super::filterbank::Analysis; +use super::tables::{alias_coefficients, synthesis_window}; +use super::*; + +/// The decoder side of the hybrid filterbank (§2.4.3.4.10: alias butterflies, IMDCT with the +/// normal window and overlap-add, frequency inversion, polyphase synthesis), for checking the +/// analysis side. +struct Synth { + overlap: [[f32; 18]; 32], + v: [f32; 1024], +} + +impl Synth { + fn new() -> Synth { + Synth { overlap: [[0.0; 18]; 32], v: [0.0; 1024] } + } + + fn granule(&mut self, xr: &[f32; GRANULE]) -> Vec { + let mut x = *xr; + for sb in 1..32 { + for (i, (cs, ca)) in alias_coefficients().iter().enumerate() { + let (up, dn) = (18 * sb - 1 - i, 18 * sb + i); + let (bu, bd) = (x[up], x[dn]); + x[up] = bu * cs - bd * ca; + x[dn] = bd * cs + bu * ca; + } + } + let mut sub = [[0f32; 18]; 32]; + for sb in 0..32 { + let mut y = [0f32; 36]; + for (i, yi) in y.iter_mut().enumerate() { + let w = (std::f64::consts::PI / 36.0 * (i as f64 + 0.5)).sin(); + let mut acc = 0f64; + for k in 0..18 { + acc += f64::from(x[18 * sb + k]) * (std::f64::consts::PI / 72.0 * (2 * i + 1 + 18) as f64 * (2 * k + 1) as f64).cos(); + } + *yi = (acc * w) as f32; + } + for t in 0..18 { + let mut v = y[t] + self.overlap[sb][t]; + if sb % 2 == 1 && t % 2 == 1 { + v = -v; + } + sub[sb][t] = v; + self.overlap[sb][t] = y[t + 18]; + } + } + let d = synthesis_window(); + let mut out = Vec::with_capacity(GRANULE); + for t in 0..18 { + self.v.copy_within(0..960, 64); + for i in 0..64 { + self.v[i] = (0..32).map(|k| (((16 + i) * (2 * k + 1)) as f64 * std::f64::consts::PI / 64.0).cos() as f32 * sub[k][t]).sum(); + } + for j in 0..32 { + let mut acc = 0f32; + for i in 0..8 { + acc += self.v[i * 128 + j] * d[i * 64 + j] + self.v[i * 128 + 96 + j] * d[i * 64 + 32 + j]; + } + out.push(acc); + } + } + out + } +} + +fn noise(n: usize, amp: f32, seed: u64) -> Vec { + let mut s = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + (0..n) + .map(|_| { + s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + ((s >> 40) as f32 / (1u64 << 24) as f32 * 2.0 - 1.0) * amp + }) + .collect() +} + +fn sine(n: usize, rate: f32, hz: f32, amp: f32) -> Vec { + (0..n).map(|i| amp * (2.0 * std::f32::consts::PI * hz * i as f32 / rate).sin()).collect() +} + +fn analyse(pcm: &[f32]) -> Vec<[f32; GRANULE]> { + let mut a = Analysis::new(); + pcm.as_chunks::() + .0 + .iter() + .map(|c| { + let mut x = [0f32; GRANULE]; + a.granule(c, &mut x); + x + }) + .collect() +} + +#[test] +fn hybrid_filterbank_reconstructs_with_the_documented_delay() { + let input = noise(GRANULE * 12, 0.5, 3); + let mut s = Synth::new(); + let out: Vec = analyse(&input).iter().flat_map(|x| s.granule(x)).collect(); + let d = DELAY as usize; + let err = (0..input.len() - d).map(|i| (out[i + d] - input[i]).abs()).fold(0f32, f32::max); + assert!(err < 2e-4, "max reconstruction error {err} at delay {d}"); +} + +#[test] +fn full_scale_tone_energy_matches_the_threshold_scale() { + // the MDCT is not shift invariant: average a few granules + let x = analyse(&sine(GRANULE * 10, 44_100.0, 1000.0, 1.0)); + let e = x[2..].iter().flatten().map(|v| v * v).sum::() / 8.0; + assert!((e / FULL_SCALE_TONE_ENERGY - 1.0).abs() < 0.1, "tone energy {e}"); +} + +/// Frames (MPEG-1, 44.1 kHz mono, 320 kbit/s, no reservoir) carrying the given granules. +fn frames_of(granules: &[Granule]) -> Vec { + let mut enc = Encoder::new(Config::cbr(44_100, 1, 320)).unwrap(); + let mut out = Vec::new(); + for pair in granules.chunks(2) { + let mut coded = pair.to_vec(); + if coded.len() == 1 { + coded.push(enc.model.encode(&[0.0; GRANULE], &[f32::MAX; 22], 0, 0)); + } + let len = enc.next_frame_len(); + let slot = len - 4 - enc.side_info_len(); + let mut md = BitWriter::new(); + coded.iter().for_each(|g| g.write(&mut md, &enc.model)); + let mut main = md.finish(); + assert!(main.len() <= slot, "{} main-data bytes in a {slot}-byte slot", main.len()); + main.resize(slot, 0); + out.extend(enc.header_and_side_info(len, false, 0, &coded)); + out.extend(main); + } + out +} + +/// A granule coding `pairs` with table `table` in all three regions, then `quads` (magnitudes +/// ≤ 1) with count1 table `c1`; every third nonzero line negative. +fn forced_granule(model: &Model, pairs: &[(u16, u16)], quads: &[[u16; 4]], table: usize, c1: usize, gain: u32) -> Granule { + let mut g = model.encode(&[0.0; GRANULE], &[f32::MAX; 22], 0, 0); + for (k, (x, y)) in pairs.iter().enumerate() { + g.ix[2 * k] = *x; + g.ix[2 * k + 1] = *y; + } + for (k, q) in quads.iter().enumerate() { + g.ix[2 * pairs.len() + 4 * k..2 * pairs.len() + 4 * k + 4].copy_from_slice(q); + } + for (i, n) in g.negative.iter_mut().enumerate() { + *n = i % 3 == 1; + } + g.global_gain = gain; + g.coding = huffman::Coding { + big_values: pairs.len(), + count1: quads.len(), + table_select: [table; 3], + region0_count: 7, + region1_count: 7, + count1_table: c1, + bits: 0, + }; + let mut w = BitWriter::new(); + huffman::write(&mut w, &g.ix, &g.negative, &g.coding, &model.bands); + g.coding.bits = w.bit_len(); + assert!(g.bits() <= 4095, "table {table}: {} bits", g.bits()); + g +} + +/// What a standard decoder outputs for the granules (mono, no scalefactors). +fn reference_decode(granules: &[Granule]) -> Vec { + let mut s = Synth::new(); + let p43 = super::tables::pow43(); + let mut out = Vec::new(); + for g in granules { + let step = 2f32.powf(0.25 * (g.global_gain as f32 - 210.0)); + let mut xr = [0f32; GRANULE]; + for (i, x) in xr.iter_mut().enumerate() { + let v = p43[usize::from(g.ix[i])] * step; + *x = if g.negative[i] { -v } else { v }; + } + out.extend(s.granule(&xr)); + } + out +} + +fn ffmpeg_decode_mono(bytes: &[u8], name: &str) -> Option<(Vec, String)> { + let ff = craft_media_testkit::ffmpeg_or_skip("mp3 Huffman tables")?; + let path = craft_media_testkit::fixtures_dir("mp3enc").join(name); + std::fs::write(&path, bytes).unwrap(); + let out = + std::process::Command::new(ff).args(["-v", "error", "-i"]).arg(&path).args(["-f", "f32le", "-acodec", "pcm_f32le", "-"]).output().unwrap(); + let pcm = out.stdout.as_chunks::<4>().0.iter().map(|b| f32::from_le_bytes(*b)).collect(); + Some((pcm, String::from_utf8_lossy(&out.stderr).into_owned())) +} + +/// Every code of every pair table (escapes with several linbits values) and of both count1 +/// tables, as ffmpeg decodes them, against the reference decode of the intended lines. +#[test] +fn every_huffman_code_decodes_as_intended_in_ffmpeg() { + let enc = Encoder::new(Config::cbr(44_100, 1, 320)).unwrap(); + let mut granules = Vec::new(); + let mut owners = Vec::new(); + for table in (1..32).filter(|t| *t != 4 && *t != 14) { + let size = huffman_tables::pair_table(table).unwrap().size as u16; + let linbits = huffman_tables::LINBITS[table]; + // every code, with small escape values so that any decoding error stands out ... + let esc = |v: u16, k: u16, big: bool| match (linbits > 0 && v == 15, big) { + (false, _) => v, + (true, false) => 15 + k % 8 % (1 << linbits), + (true, true) => 15 + ((1u32 << linbits) - 1).saturating_sub(u32::from(k % 3)) as u16, + }; + let mut pairs = Vec::new(); + for x in 0..size { + for y in 0..size { + pairs.push((esc(x, x + 3 * y, false), esc(y, 7 * x + y, false))); + } + } + for chunk in pairs.chunks(48) { + granules.push(forced_granule(&enc.model, chunk, &[], table, 0, 190)); + owners.push(format!("table {table}")); + } + // ... and the largest escapes + if linbits > 0 { + let big: Vec<(u16, u16)> = (0..16).map(|k| (esc(15, k, true), esc(if k % 2 == 0 { 15 } else { k }, k + 1, true))).collect(); + granules.push(forced_granule(&enc.model, &big, &[], table, 0, 120)); + owners.push(format!("table {table}, large escapes")); + } + } + let quads: Vec<[u16; 4]> = (0..16u16).map(|q| [q >> 3, (q >> 2) & 1, (q >> 1) & 1, q & 1]).collect(); + for c1 in 0..2 { + granules.push(forced_granule(&enc.model, &[(1, 2), (0, 1)], &quads, 2, c1, 190)); + owners.push(format!("count1 table {}", if c1 == 0 { "A" } else { "B" })); + } + let bytes = frames_of(&granules); + let Some((ff, err)) = ffmpeg_decode_mono(&bytes, "huffman_tables.mp3") else { return }; + assert!(err.trim().is_empty(), "ffmpeg: {err}"); + let ours = reference_decode(&granules); + assert!(ff.len() >= ours.len(), "{} vs {}", ff.len(), ours.len()); + // a wrong code desynchronises the rest of its granule: report the first granule that differs + for (k, owner) in owners.iter().enumerate() { + let r = k * GRANULE..(k + 1) * GRANULE; + let peak = ours[r.clone()].iter().fold(1e-3f32, |m, v| m.max(v.abs())); + let diff = ours[r.clone()].iter().zip(&ff[r]).map(|(a, b)| (a - b).abs()).fold(0f32, f32::max); + assert!(diff < 1e-5 * peak, "{owner} (granule {k}): differs by {diff} (peak {peak})"); + } +} + +/// A frame's header and side-information fields we check. +#[derive(Debug)] +struct Parsed { + lsf: bool, + bitrate_index: u32, + rate_index: u32, + padding: bool, + mode: u32, + mode_extension: u32, + len: usize, + main_data_begin: usize, + part2_3: Vec, + global_gain: Vec, +} + +/// Walk the frames of a stream, checking sync and lengths. +fn parse(bytes: &[u8]) -> Vec { + let mut out = Vec::new(); + let mut pos = 0; + while pos < bytes.len() { + let h = u32::from_be_bytes(bytes[pos..pos + 4].try_into().unwrap()); + assert_eq!(h >> 21, 0x7FF, "sync at {pos}"); + let lsf = (h >> 19) & 3 == 2; + assert!(lsf || (h >> 19) & 3 == 3, "MPEG-1 or MPEG-2"); + assert_eq!((h >> 17) & 3, 1, "Layer III"); + assert_eq!((h >> 16) & 1, 1, "no CRC"); + let bitrate_index = (h >> 12) & 15; + let rate_index = (h >> 10) & 3; + let padding = (h >> 9) & 1 == 1; + let mode = (h >> 6) & 3; + let kbps = if lsf { tables::BITRATES_LSF } else { tables::BITRATES_MPEG1 }[bitrate_index as usize - 1]; + let rate = if lsf { tables::RATES_LSF } else { tables::RATES_MPEG1 }[rate_index as usize]; + let len = (if lsf { 72_000 } else { 144_000 } * kbps / rate) as usize + usize::from(padding); + let nch = if mode == 3 { 1 } else { 2 }; + let mut r = craft_media::bitstream::BitReader::new(&bytes[pos + 4..pos + len]); + let main_data_begin = r.read_bits(if lsf { 8 } else { 9 }).unwrap() as usize; + r.skip(if lsf { nch } else { [5, 3][nch - 1] + 4 * nch }).unwrap(); + let (mut part2_3, mut global_gain) = (Vec::new(), Vec::new()); + for _ in 0..(if lsf { 1 } else { 2 }) * nch { + part2_3.push(r.read_bits(12).unwrap() as usize); + let big_values = r.read_bits(9).unwrap(); + assert!(big_values <= 288); + global_gain.push(r.read_bits(8).unwrap()); + r.skip(if lsf { 9 } else { 4 }).unwrap(); + assert!(!r.read_flag().unwrap(), "long blocks only"); + r.skip(15 + 7 + if lsf { 2 } else { 3 }).unwrap(); + } + out.push(Parsed { lsf, bitrate_index, rate_index, padding, mode, mode_extension: (h >> 4) & 3, len, main_data_begin, part2_3, global_gain }); + pos += len; + } + assert_eq!(pos, bytes.len(), "whole frames"); + out +} + +fn encode_all(cfg: Config, input: &[Vec]) -> Vec { + let mut enc = Encoder::new(cfg).unwrap(); + let planes: Vec<&[f32]> = input.iter().map(Vec::as_slice).collect(); + let mut out = enc.encode(&planes); + out.extend(enc.flush()); + out +} + +#[test] +fn headers_carry_the_configuration_and_frames_have_their_nominal_sizes() { + let n = 44_100; + let input = vec![sine(n, 44_100.0, 440.0, 0.5), sine(n, 44_100.0, 660.0, 0.4)]; + let frames = parse(&encode_all(Config::cbr(44_100, 2, 128), &input)); + // the decoder needs (input + delay) samples: whole frames of 1152 + assert_eq!(frames.len(), (n + DELAY as usize).div_ceil(1152)); + let f = &frames[0]; + assert_eq!((f.lsf, f.bitrate_index, f.rate_index, f.mode), (false, 9, 0, 1)); + // 44.1 kHz: 417 or 418 bytes, padded so that the mean is exactly 128 kbit/s + let total: usize = frames.iter().map(|f| f.len).sum(); + let exact = 144_000.0 * 128.0 / 44_100.0 * frames.len() as f64; + assert!((total as f64 - exact).abs() < 1.0, "{total} vs {exact}"); + assert!(frames.iter().any(|f| f.padding) && frames.iter().any(|f| !f.padding)); + // the same tone in both channels at nearly the same level: mid/side + let ms = parse(&encode_all(Config::cbr(44_100, 2, 128), &[input[0].clone(), input[0].iter().map(|v| v * 0.9).collect()])); + assert!(ms.iter().skip(1).all(|f| f.mode_extension == 2), "mid/side expected"); + // simple stereo never uses mid/side + let st = parse(&encode_all(Config { stereo: StereoMode::Stereo, ..Config::cbr(48_000, 2, 192) }, &input)); + assert!(st.iter().all(|f| f.mode == 0 && f.mode_extension == 0)); + // MPEG-2 LSF mono + let lsf = parse(&encode_all(Config::cbr(22_050, 1, 32), &input[..1])); + assert_eq!((lsf[0].lsf, lsf[0].rate_index, lsf[0].mode, lsf[0].part2_3.len()), (true, 0, 3, 1)); + assert_eq!(lsf.len(), (n + DELAY as usize).div_ceil(576)); +} + +#[test] +fn the_reservoir_lends_bits_within_the_decoder_buffer() { + // quiet passages save bits that a loud burst of noise then borrows + let mut x = vec![0f32; 44_100]; + let burst = noise(4410, 0.7, 9); + x[22_050..26_460].copy_from_slice(&burst); + x.iter_mut().take(22_050).enumerate().for_each(|(i, v)| *v = 0.01 * (i as f32 * 0.05).sin()); + for (kbps, max_begin) in [(64u32, 511usize), (128, 511), (256, 125), (320, 0)] { + let frames = parse(&encode_all(Config::cbr(44_100, 1, kbps), &[x.clone()])); + let mut slots = 0usize; + let mut lent = 0; + for f in &frames { + let slot = f.len - 4 - 17; + assert!(f.main_data_begin <= max_begin, "{kbps} kbps: main_data_begin {}", f.main_data_begin); + assert!(f.main_data_begin <= slots, "{kbps} kbps: points before the stream"); + let bits: usize = f.part2_3.iter().sum(); + assert!(bits <= 8 * (f.main_data_begin + slot), "{kbps} kbps: {bits} bits of main data"); + assert!(f.part2_3.iter().all(|&b| b <= 4095)); + slots += slot; + lent = lent.max(f.main_data_begin); + } + if max_begin > 0 { + assert!(lent > 0, "{kbps} kbps: the reservoir was never used"); + } + } +} + +#[test] +fn hostile_settings_are_refused() { + for rate in [0, 8_000, 11_025, 44_101, 96_000, u32::MAX] { + assert_eq!(Encoder::new(Config::cbr(rate, 2, 128)).err(), Some(Error::SampleRate(rate))); + } + for kbps in [0, 8, 33, 144, 384, u32::MAX] { + assert!(matches!(Encoder::new(Config::cbr(44_100, 2, kbps)), Err(Error::Bitrate { .. })), "{kbps}"); + } + assert!(Encoder::new(Config::cbr(24_000, 1, 8)).is_ok()); + assert!(matches!(Encoder::new(Config::cbr(24_000, 1, 192)), Err(Error::Bitrate { .. }))); + for ch in [0, 3, 6, usize::MAX] { + assert_eq!(Encoder::new(Config::cbr(48_000, ch, 128)).err(), Some(Error::Channels(ch))); + } + assert_eq!(nearest_bitrate(44_100, 1000), Some(320)); + assert_eq!(nearest_bitrate(44_100, 130), Some(128)); + assert_eq!(nearest_bitrate(44_100, 0), Some(32)); + assert_eq!(nearest_bitrate(22_050, 320), Some(160)); + assert_eq!(nearest_bitrate(96_000, 128), None); + assert_eq!( + (nearest_sample_rate(44_100), nearest_sample_rate(96_000), nearest_sample_rate(8_000), nearest_sample_rate(0)), + (44_100, 48_000, 16_000, 16_000) + ); + assert_eq!(nearest_sample_rate(u32::MAX), 48_000); + assert!(bitrates(0).is_empty()); +} + +#[test] +fn hostile_samples_still_make_valid_frames() { + let mut enc = Encoder::new(Config::cbr(48_000, 2, 96)).unwrap(); + let mut out = Vec::new(); + let wild: Vec = (0..5000) + .map(|i| match i % 5 { + 0 => f32::NAN, + 1 => f32::INFINITY, + 2 => -1e30, + 3 => 1e-40, + _ => 3.0, + }) + .collect(); + out.extend(enc.encode(&[&wild, &wild[..10]])); // ragged: the second channel is padded with silence + out.extend(enc.encode(&[&wild[..7]])); // a missing channel is silence + out.extend(enc.encode(&[])); + out.extend(enc.encode(&[&[], &[], &[]])); + out.extend(enc.flush()); + out.extend(enc.flush()); // a second flush adds nothing + let frames = parse(&out); + assert_eq!(frames.len(), (5007 + DELAY as usize).div_ceil(1152)); + assert!(frames.iter().all(|f| f.global_gain.iter().all(|&g| g <= 255))); + // no input, no frames + let mut empty = Encoder::new(Config::cbr(44_100, 1, 64)).unwrap(); + assert!(empty.flush().is_empty()); +} + +#[test] +fn chunking_does_not_change_the_stream() { + let input = vec![sine(20_000, 32_000.0, 300.0, 0.3), noise(20_000, 0.1, 5)]; + let whole = encode_all(Config::cbr(32_000, 2, 96), &input); + let mut enc = Encoder::new(Config::cbr(32_000, 2, 96)).unwrap(); + let mut parts = Vec::new(); + for k in (0..20_000).step_by(777) { + let e = (k + 777).min(20_000); + parts.extend(enc.encode(&[&input[0][k..e], &input[1][k..e]])); + } + parts.extend(enc.flush()); + assert!(whole == parts, "chunked encode differs"); +} + +/// Every sample rate and bit rate in mono, and the lowest and highest bit rate of each sample +/// rate in both stereo modes, fed the signals that stress the +/// quantizer (full-scale and clipped noise, an impulse train, DC, the Nyquist tone, denormals and +/// silence): whole frames of the nominal size, every granule within its 12-bit length and the +/// reservoir within the frames before it. +#[test] +fn every_setting_makes_valid_frames_from_extreme_signals() { + let signals: [fn(usize, u64) -> Vec; 7] = [ + |n, s| noise(n, 1.0, s), + |n, s| noise(n, 40.0, s), // clamped to ±4 by the encoder + |n, _| (0..n).map(|i| if i % 97 == 0 { 1.0 } else { 0.0 }).collect(), + |n, _| vec![1.0; n], + |n, _| (0..n).map(|i| if i % 2 == 0 { 1.0 } else { -1.0 }).collect(), + |n, _| vec![1e-42; n], + |n, _| vec![0.0; n], + ]; + let mut k = 0usize; + for rate in SAMPLE_RATES { + let rates = bitrates(rate); + for (i, &kbps) in rates.iter().enumerate() { + let edge = i == 0 || i + 1 == rates.len(); + for (channels, stereo) in [(1, StereoMode::Stereo), (2, StereoMode::Stereo), (2, StereoMode::JointStereo)] { + if channels == 2 && !edge { + continue; + } + k += 1; + // two or three frames per setting keep a debug build quick + let n = 600; + let input: Vec> = (0..channels).map(|c| signals[(k + 3 * c) % signals.len()](n, (k + c) as u64)).collect(); + let cfg = Config { sample_rate: rate, channels, bitrate_kbps: kbps, stereo }; + let what = format!("{rate} Hz {kbps} kbit/s {channels} ch {stereo:?}"); + let frames = parse(&encode_all(cfg, &input)); + let fs = if rate >= 32_000 { 1152 } else { 576 }; + assert_eq!(frames.len(), (n + DELAY as usize).div_ceil(fs), "{what}"); + let side = match (rate >= 32_000, channels) { + (true, 1) => 17, + (true, _) => 32, + (false, 1) => 9, + (false, _) => 17, + }; + let mut before = 0usize; + for f in &frames { + assert_eq!(f.lsf, rate < 32_000, "{what}"); + assert_eq!( + f.mode, + if channels == 1 { + 3 + } else if stereo == StereoMode::Stereo { + 0 + } else { + 1 + }, + "{what}" + ); + let slot = f.len - 4 - side; + assert!(f.main_data_begin <= before, "{what}: main_data_begin {} points before the stream", f.main_data_begin); + assert!(f.part2_3.iter().all(|&b| b <= 4095), "{what}: {:?}", f.part2_3); + assert!(f.part2_3.iter().sum::() <= 8 * (f.main_data_begin + slot), "{what}: main data overruns its frame"); + before += slot; + } + } + } + } + assert_eq!(k, 6 * 14 + 6 * 2 * 2); +} diff --git a/crates/audio/tests/codecs/main.rs b/crates/audio/tests/codecs/main.rs index 030f3deb..0fafa5dd 100644 --- a/crates/audio/tests/codecs/main.rs +++ b/crates/audio/tests/codecs/main.rs @@ -12,6 +12,7 @@ mod ac3_encoder_oracle; mod ac3_oracle; mod alac; mod concealment; +mod mp3_encoder_oracle; mod mpa_oracle; mod opus; mod routing_oracle; diff --git a/crates/audio/tests/codecs/mp3_encoder_oracle.rs b/crates/audio/tests/codecs/mp3_encoder_oracle.rs new file mode 100644 index 00000000..c9bbf24f --- /dev/null +++ b/crates/audio/tests/codecs/mp3_encoder_oracle.rs @@ -0,0 +1,162 @@ +//! The encoder against ffmpeg as an external decoder (skipped without ffmpeg): ffprobe must see +//! an MP3 stream with our rate, channels and duration, ffmpeg must decode every frame without a +//! complaint, and its decode, aligned by [`DELAY`], must match the input within a signal-to-noise +//! ratio per bit rate. White noise loses what lies above the bit rate's low-pass, so its bound is +//! low; the bounds sit about 6 dB under what the encoder reaches today. + +use std::path::{Path, PathBuf}; +use std::process::{Command, Stdio}; + +use craft_audio::mp3::{Config, DELAY, Encoder, StereoMode}; + +fn fixtures() -> PathBuf { + crate::testkit::fixtures_dir("mp3enc") +} + +fn chords(rate: u32, n: usize, detune: f32) -> Vec { + let f = [220.0, 277.18, 329.63, 440.0, 1318.5, 3520.0]; + let a = [0.2, 0.15, 0.12, 0.1, 0.05, 0.02]; + (0..n) + .map(|i| { + let t = i as f32 / rate as f32; + f.iter().zip(a).map(|(f, a)| a * (2.0 * std::f32::consts::PI * f * (1.0 + detune) * t).sin()).sum::() + }) + .collect() +} + +fn noise(n: usize, amp: f32, seed: u64) -> Vec { + let mut s = seed.wrapping_mul(0x9E37_79B9_7F4A_7C15) | 1; + (0..n) + .map(|_| { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + ((s >> 40) as f32 / (1u64 << 24) as f32 * 2.0 - 1.0) * amp + }) + .collect() +} + +fn sweep(rate: u32, n: usize) -> Vec { + let secs = n as f32 / rate as f32; + let (f0, f1) = (40.0f32, rate as f32 * 0.4); + let k = (f1 / f0).ln() / secs; + (0..n).map(|i| 0.4 * (2.0 * std::f32::consts::PI * f0 * (((i as f32 / rate as f32) * k).exp() - 1.0) / k).sin()).collect() +} + +fn encode(cfg: Config, input: &[Vec]) -> Vec { + let mut enc = Encoder::new(cfg).unwrap(); + let mut out = Vec::new(); + let n = input[0].len(); + let mut pos = 0; + while pos < n { + let end = (pos + 1000).min(n); + let planes: Vec<&[f32]> = input.iter().map(|c| &c[pos..end]).collect(); + out.extend(enc.encode(&planes)); + pos = end; + } + out.extend(enc.flush()); + out +} + +fn ffmpeg_decode(ff: &Path, path: &Path, channels: usize) -> (Vec>, String) { + let out = Command::new(ff) + .args(["-v", "error", "-i"]) + .arg(path) + .args(["-f", "f32le", "-acodec", "pcm_f32le", "-"]) + .stdin(Stdio::null()) + .output() + .unwrap(); + let samples: Vec = out.stdout.as_chunks::<4>().0.iter().map(|b| f32::from_le_bytes(*b)).collect(); + let planar = (0..channels).map(|c| samples.iter().skip(c).step_by(channels).copied().collect()).collect(); + (planar, String::from_utf8_lossy(&out.stderr).into_owned()) +} + +/// SNR in dB of `dec` (aligned) against `input`. +fn snr(input: &[f32], dec: &[f32]) -> f64 { + let (mut s, mut n) = (0f64, 0f64); + for (a, b) in input.iter().zip(dec) { + s += f64::from(*a) * f64::from(*a); + n += f64::from(a - b) * f64::from(a - b); + } + 10.0 * (s / n.max(1e-30)).log10() +} + +struct Case { + name: &'static str, + rate: u32, + channels: usize, + kbps: u32, + stereo: StereoMode, + signal: fn(u32, usize, usize) -> Vec, + min_snr: f64, +} + +fn sig_chords(rate: u32, n: usize, c: usize) -> Vec { + chords(rate, n, 0.003 * c as f32) +} +fn sig_noise(_rate: u32, n: usize, c: usize) -> Vec { + noise(n, 0.25, 11 + c as u64) +} +fn sig_sweep(rate: u32, n: usize, c: usize) -> Vec { + sweep(rate, n).iter().map(|v| v * (1.0 - 0.2 * c as f32)).collect() +} + +fn run(ff: &Path, c: &Case) -> f64 { + let n = c.rate as usize * 2 + 777; + let input: Vec> = (0..c.channels).map(|ch| (c.signal)(c.rate, n, ch)).collect(); + let cfg = Config { sample_rate: c.rate, channels: c.channels, bitrate_kbps: c.kbps, stereo: c.stereo }; + let bytes = encode(cfg, &input); + let path = fixtures().join(format!("{}_{}_{}ch_{}k.mp3", c.name, c.rate, c.channels, c.kbps)); + std::fs::write(&path, &bytes).unwrap(); + let (dec, err) = ffmpeg_decode(ff, &path, c.channels); + assert!(err.trim().is_empty(), "{}: ffmpeg reported: {err}", path.display()); + let d = DELAY as usize; + assert!(dec[0].len() >= n + d, "{}: decoded {} samples, need {}", c.name, dec[0].len(), n + d); + // the bit rate is exact: every frame has the nominal size + let frame_samples = if c.rate >= 32_000 { 1152 } else { 576 }; + let frames = dec[0].len() / frame_samples; + let nominal = c.kbps as f64 * 1000.0 / 8.0 * (frames * frame_samples) as f64 / c.rate as f64; + assert!((bytes.len() as f64 / nominal - 1.0).abs() < 0.01, "{}: {} bytes vs {nominal}", c.name, bytes.len()); + if let Some(probe) = crate::testkit::oracle::find_tool("ffprobe") { + let out = Command::new(probe) + .args(["-v", "error", "-show_entries", "stream=codec_name,sample_rate,channels:format=duration", "-of", "default=nw=1"]) + .arg(&path) + .output() + .unwrap(); + let text = String::from_utf8_lossy(&out.stdout); + let field = |k: &str| text.lines().find_map(|l| l.strip_prefix(&format!("{k}="))).unwrap_or("").to_string(); + assert_eq!(field("codec_name"), "mp3", "{text}"); + assert_eq!(field("sample_rate"), c.rate.to_string(), "{text}"); + assert_eq!(field("channels"), c.channels.to_string(), "{text}"); + let dur: f64 = field("duration").parse().unwrap(); + let expect = (frames * frame_samples) as f64 / c.rate as f64; + assert!((dur - expect).abs() < 2.0 * frame_samples as f64 / c.rate as f64, "{}: duration {dur} vs {expect}", c.name); + } + let s = (0..c.channels).map(|ch| snr(&input[ch], &dec[ch][d..d + n])).sum::() / c.channels as f64; + println!("{:<8} {:>5} Hz {}ch {:>3} kbps {:?}: SNR {s:.2} dB", c.name, c.rate, c.channels, c.kbps, c.stereo); + s +} + +#[test] +fn ffmpeg_decodes_our_streams_at_every_rate() { + let ff = require_ffmpeg!(); + use StereoMode::*; + let cases = [ + Case { name: "chords", rate: 44_100, channels: 2, kbps: 128, stereo: JointStereo, signal: sig_chords, min_snr: 42.0 }, + Case { name: "chords", rate: 44_100, channels: 2, kbps: 320, stereo: JointStereo, signal: sig_chords, min_snr: 60.0 }, + Case { name: "sweep", rate: 48_000, channels: 2, kbps: 192, stereo: Stereo, signal: sig_sweep, min_snr: 22.0 }, + Case { name: "noise", rate: 48_000, channels: 2, kbps: 128, stereo: JointStereo, signal: sig_noise, min_snr: 2.5 }, + Case { name: "chords", rate: 32_000, channels: 1, kbps: 64, stereo: JointStereo, signal: sig_chords, min_snr: 40.0 }, + Case { name: "chords", rate: 22_050, channels: 2, kbps: 64, stereo: JointStereo, signal: sig_chords, min_snr: 33.0 }, + Case { name: "sweep", rate: 24_000, channels: 1, kbps: 48, stereo: JointStereo, signal: sig_sweep, min_snr: 24.0 }, + Case { name: "chords", rate: 16_000, channels: 1, kbps: 32, stereo: JointStereo, signal: sig_chords, min_snr: 30.0 }, + ]; + let mut failed = Vec::new(); + for c in &cases { + let s = run(&ff, c); + if s < c.min_snr { + failed.push(format!("{} {} Hz {} kbps: {s:.2} < {}", c.name, c.rate, c.kbps, c.min_snr)); + } + } + assert!(failed.is_empty(), "{failed:#?}"); +}