From 78c263a238ee7d6b5ce76e952bd78a166f5c9a46 Mon Sep 17 00:00:00 2001 From: Brandon Thomas Date: Sat, 3 Oct 2026 04:11:27 -0400 Subject: [PATCH] M9.15: AC-3 decoder from ATSC A/52 (all modes, coupling, rematrixing, block switching), wired into TS/PS/MP4/MKV Co-Authored-By: Claude Opus 5.5 --- Cargo.lock | 9 + Cargo.toml | 1 + crates/ac3/Cargo.toml | 16 + crates/ac3/src/lib.rs | 1054 ++++++++++++++++++++++++++++ crates/ac3/src/tables.rs | 69 ++ crates/ac3/src/tests.rs | 98 +++ crates/ac3/tests/oracle.rs | 147 ++++ crates/codecs/Cargo.toml | 1 + crates/codecs/src/audio.rs | 26 +- crates/codecs/src/mkv.rs | 1 + crates/codecs/tests/mpeg_oracle.rs | 7 +- 11 files changed, 1425 insertions(+), 4 deletions(-) create mode 100644 crates/ac3/Cargo.toml create mode 100644 crates/ac3/src/lib.rs create mode 100644 crates/ac3/src/tables.rs create mode 100644 crates/ac3/src/tests.rs create mode 100644 crates/ac3/tests/oracle.rs diff --git a/Cargo.lock b/Cargo.lock index 5d0e90bd3..53aba9696 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1684,6 +1684,14 @@ dependencies = [ "thiserror 2.0.21", ] +[[package]] +name = "filmcraft-ac3" +version = "0.1.0" +dependencies = [ + "filmcraft-testkit", + "thiserror 2.0.21", +] + [[package]] name = "filmcraft-audio-dsp" version = "0.1.0" @@ -1759,6 +1767,7 @@ name = "filmcraft-codecs" version = "0.1.0" dependencies = [ "filmcraft-aac", + "filmcraft-ac3", "filmcraft-av1", "filmcraft-bitstream", "filmcraft-color", diff --git a/Cargo.toml b/Cargo.toml index f9b43019c..02dfdf88b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -21,6 +21,7 @@ filmcraft-matroska = { path = "crates/matroska" } filmcraft-mxf = { path = "crates/mxf" } filmcraft-mpegts = { path = "crates/mpegts" } filmcraft-mpeg2v = { path = "crates/mpeg2v" } +filmcraft-ac3 = { path = "crates/ac3" } filmcraft-ogg = { path = "crates/ogg" } filmcraft-h264 = { path = "crates/h264" } filmcraft-hevc = { path = "crates/hevc" } diff --git a/crates/ac3/Cargo.toml b/crates/ac3/Cargo.toml new file mode 100644 index 000000000..e45965968 --- /dev/null +++ b/crates/ac3/Cargo.toml @@ -0,0 +1,16 @@ +[package] +name = "filmcraft-ac3" +description = "Clean-room AC-3 (Dolby Digital, ATSC A/52) audio decoder" +version.workspace = true +edition.workspace = true +license.workspace = true +rust-version.workspace = true + +[lints] +workspace = true + +[dependencies] +thiserror = { workspace = true } + +[dev-dependencies] +filmcraft-testkit = { workspace = true } diff --git a/crates/ac3/src/lib.rs b/crates/ac3/src/lib.rs new file mode 100644 index 000000000..ff06ea301 --- /dev/null +++ b/crates/ac3/src/lib.rs @@ -0,0 +1,1054 @@ +//! Clean-room AC-3 (Dolby Digital) audio decoder, written from ATSC A/52:2012 ("Digital Audio +//! Compression Standard", §5 bit stream syntax, §6-7 decoding). +//! +//! - All audio coding modes (1+1, 1/0 … 3/2) with or without LFE, 32 / 44.1 / 48 kHz, every +//! frame size (Table 5.18). +//! - Exponents (D15 / D25 / D45, reuse across blocks), the parametric bit allocation exactly as +//! specified (§7.2, fixed-point integer steps and Tables 7.6-7.16) including delta bit allocation, +//! grouped and asymmetric mantissas, dither for zero-bit mantissas (§7.3.4), channel coupling +//! with phase flags, rematrixing, dynamic range control (dynrng / dynrng2), block switching +//! (512- and 256-sample IMDCTs, §7.9.4) with the Kaiser-Bessel derived window. +//! - Output: planar f32 at full scale ±1, channels in WAV / SMPTE order (L R C LFE Ls Rs); no +//! downmix. +//! +//! E-AC-3 (bsid 11-16) is recognised and refused. + +mod tables; + +use std::sync::OnceLock; + +#[derive(Debug, thiserror::Error, PartialEq, Eq)] +pub enum Error { + #[error("not an AC-3 frame")] + NoSync, + #[error("unsupported AC-3: {0}")] + Unsupported(String), + #[error("invalid AC-3 frame: {0}")] + Invalid(&'static str), +} + +pub type Result = std::result::Result; + +/// Data rate in kb/s by frmsizecod / 2 (Table 5.18). +const BITRATES: [u32; 19] = [32, 40, 48, 56, 64, 80, 96, 112, 128, 160, 192, 224, 256, 320, 384, 448, 512, 576, 640]; + +/// Full-bandwidth channels per acmod (Table 5.8). +const NFCHANS: [usize; 8] = [2, 1, 2, 3, 3, 4, 4, 5]; + +/// syncinfo and the start of bsi. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub struct Header { + pub sample_rate: u32, + pub bitrate_kbps: u32, + /// Bytes in the syncframe. + pub frame_bytes: usize, + pub bsid: u8, + pub acmod: u8, + pub lfeon: bool, +} + +impl Header { + /// Output channels (full bandwidth + LFE). + pub fn channels(&self) -> usize { + NFCHANS[self.acmod as usize] + self.lfeon as usize + } +} + +/// Parse the frame header at the start of `b`. +pub fn parse_header(b: &[u8]) -> Result
{ + if b.len() < 7 || b[0] != 0x0B || b[1] != 0x77 { + return Err(Error::NoSync); + } + let fscod = b[4] >> 6; + let frmsizecod = b[4] & 0x3F; + let bsid = b[5] >> 3; + if bsid > 10 { + return Err(Error::Unsupported(format!("bsid {bsid} (E-AC-3)"))); + } + if fscod == 3 || frmsizecod >= 38 { + return Err(Error::Invalid("reserved fscod / frmsizecod")); + } + let kbps = BITRATES[(frmsizecod / 2) as usize]; + let words = match fscod { + 0 => kbps * 2, + 1 => kbps * 96_000 / 44_100 + (frmsizecod & 1) as u32, + _ => kbps * 3, + }; + // bsid 9 / 10: half / quarter sample rate (A/52 Annex note on "bsid 9 and 10") + let shift = bsid.saturating_sub(8) as u32; + let sample_rate = [48_000, 44_100, 32_000][fscod as usize] >> shift; + let mut r = Bits::new(&b[6..]); + let acmod = r.read(3) as u8; + if acmod & 1 != 0 && acmod != 1 { + r.skip(2); + } + if acmod & 4 != 0 { + r.skip(2); + } + if acmod == 2 { + r.skip(2); + } + let lfeon = r.read(1) == 1; + Ok(Header { sample_rate, bitrate_kbps: kbps >> shift, frame_bytes: words as usize * 2, bsid, acmod, lfeon }) +} + +/// MSB-first bit reader; reads past the end give zeros. +struct Bits<'a> { + d: &'a [u8], + pos: usize, +} + +impl<'a> Bits<'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 + } +} + +/// Index of the LFE and coupling channels in the per-channel arrays. +const LFE: usize = 5; +const CPL: usize = 6; + +/// Bit allocation parameter tables (Tables 7.6-7.11). +const SLOWDEC: [i32; 4] = [0x0f, 0x11, 0x13, 0x15]; +const FASTDEC: [i32; 4] = [0x3f, 0x53, 0x67, 0x7b]; +const SLOWGAIN: [i32; 4] = [0x540, 0x4d8, 0x478, 0x410]; +const DBPBTAB: [i32; 4] = [0x000, 0x700, 0x900, 0xb00]; +const FLOORTAB: [i32; 8] = [0x2f0, 0x2b0, 0x270, 0x230, 0x1f0, 0x170, 0x0f0, -0x800]; +const FASTGAIN: [i32; 8] = [0x080, 0x100, 0x180, 0x200, 0x280, 0x300, 0x380, 0x400]; +/// Mantissa bits of the asymmetric quantizers by bap (Table 7.18). +const QNTZTAB: [u32; 16] = [0, 0, 0, 3, 0, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 16]; + +/// The Kaiser-Bessel derived window of §7.9.4 (α = 5, 512 points), whose first half Table 7.33 +/// lists to five decimals. +fn window() -> &'static [f32; 256] { + static W: OnceLock<[f32; 256]> = OnceLock::new(); + W.get_or_init(|| { + let i0 = |x: f64| { + let (mut s, mut t) = (1.0f64, 1.0f64); + for k in 1..50 { + t *= (x / 2.0) / k as f64; + s += t * t; + } + s + }; + let alpha = 5.0 * std::f64::consts::PI; + let k: Vec = (0..=256).map(|j| i0(alpha * (1.0 - ((j as f64 - 128.0) / 128.0).powi(2)).max(0.0).sqrt())).collect(); + let total: f64 = k.iter().sum(); + let mut acc = 0.0; + let mut w = [0f32; 256]; + for n in 0..256 { + acc += k[n]; + w[n] = (acc / total).sqrt() as f32; + } + w + }) +} + +/// Twiddles for the IMDCTs: (xcos1, xsin1) for N = 512 and (xcos2, xsin2) for the 256-sample +/// transforms, and the IFFT kernels. +struct Twiddles { + c1: [f32; 128], + s1: [f32; 128], + c2: [f32; 64], + s2: [f32; 64], + /// e^{j 2π m / 128}, m = 0..128. + e: [(f32, f32); 128], +} + +fn twiddles() -> &'static Twiddles { + static T: OnceLock = OnceLock::new(); + T.get_or_init(|| { + use std::f64::consts::PI; + let n = 512.0; + let mut t = Twiddles { c1: [0.0; 128], s1: [0.0; 128], c2: [0.0; 64], s2: [0.0; 64], e: [(0.0, 0.0); 128] }; + for k in 0..128 { + let a = 2.0 * PI * (8 * k + 1) as f64 / (8.0 * n); + t.c1[k] = -a.cos() as f32; + t.s1[k] = -a.sin() as f32; + let m = 2.0 * PI * k as f64 / 128.0; + t.e[k] = (m.cos() as f32, m.sin() as f32); + } + for k in 0..64 { + let a = 2.0 * PI * (8 * k + 1) as f64 / (4.0 * n); + t.c2[k] = -a.cos() as f32; + t.s2[k] = -a.sin() as f32; + } + t + }) +} + +/// In-place complex inverse DFT of size `n` (128 or 64) by radix-2 decimation in time: +/// z[n] = Σ Z[k] e^{+j 2π k n / size}. +fn ifft(re: &mut [f32], im: &mut [f32]) { + let n = re.len(); + let e = &twiddles().e; + // bit reversal + let bits = n.trailing_zeros(); + for i in 0..n { + let j = i.reverse_bits() >> (usize::BITS - bits); + if j > i { + re.swap(i, j); + im.swap(i, j); + } + } + let mut len = 2; + while len <= n { + let step = 128 / len; + for s in (0..n).step_by(len) { + for k in 0..len / 2 { + let (wr, wi) = e[k * step]; + let (a, b) = (s + k, s + k + len / 2); + let tr = re[b] * wr - im[b] * wi; + let ti = re[b] * wi + im[b] * wr; + re[b] = re[a] - tr; + im[b] = im[a] - ti; + re[a] += tr; + im[a] += ti; + } + } + len *= 2; + } +} + +/// Inverse transform of one block of 256 coefficients into 512 windowed samples (§7.9.4). +fn imdct(x: &[f32; 256], short: bool, out: &mut [f32; 512]) { + let w = window(); + let t = twiddles(); + const N: usize = 512; + if !short { + let (mut zr, mut zi) = ([0f32; 128], [0f32; 128]); + for k in 0..N / 4 { + let (a, b) = (x[N / 2 - 2 * k - 1], x[2 * k]); + zr[k] = a * t.c1[k] - b * t.s1[k]; + zi[k] = b * t.c1[k] + a * t.s1[k]; + } + ifft(&mut zr, &mut zi); + let (mut yr, mut yi) = ([0f32; 128], [0f32; 128]); + for n in 0..N / 4 { + yr[n] = zr[n] * t.c1[n] - zi[n] * t.s1[n]; + yi[n] = zi[n] * t.c1[n] + zr[n] * t.s1[n]; + } + for n in 0..N / 8 { + out[2 * n] = -yi[N / 8 + n] * w[2 * n]; + out[2 * n + 1] = yr[N / 8 - n - 1] * w[2 * n + 1]; + out[N / 4 + 2 * n] = -yr[n] * w[N / 4 + 2 * n]; + out[N / 4 + 2 * n + 1] = yi[N / 4 - n - 1] * w[N / 4 + 2 * n + 1]; + out[N / 2 + 2 * n] = -yr[N / 8 + n] * w[N / 2 - 2 * n - 1]; + out[N / 2 + 2 * n + 1] = yi[N / 8 - n - 1] * w[N / 2 - 2 * n - 2]; + out[3 * N / 4 + 2 * n] = yi[n] * w[N / 4 - 2 * n - 1]; + out[3 * N / 4 + 2 * n + 1] = -yr[N / 4 - n - 1] * w[N / 4 - 2 * n - 2]; + } + return; + } + let mut y = [([0f32; 64], [0f32; 64]); 2]; + for (h, yh) in y.iter_mut().enumerate() { + let (mut zr, mut zi) = ([0f32; 64], [0f32; 64]); + for k in 0..N / 8 { + let (a, b) = (x[2 * (N / 4 - 2 * k - 1) + h], x[2 * (2 * k) + h]); + zr[k] = a * t.c2[k] - b * t.s2[k]; + zi[k] = b * t.c2[k] + a * t.s2[k]; + } + ifft(&mut zr, &mut zi); + for n in 0..N / 8 { + yh.0[n] = zr[n] * t.c2[n] - zi[n] * t.s2[n]; + yh.1[n] = zi[n] * t.c2[n] + zr[n] * t.s2[n]; + } + } + let ((yr1, yi1), (yr2, yi2)) = (y[0], y[1]); + for n in 0..N / 8 { + out[2 * n] = -yi1[n] * w[2 * n]; + out[2 * n + 1] = yr1[N / 8 - n - 1] * w[2 * n + 1]; + out[N / 4 + 2 * n] = -yr1[n] * w[N / 4 + 2 * n]; + out[N / 4 + 2 * n + 1] = yi1[N / 8 - n - 1] * w[N / 4 + 2 * n + 1]; + out[N / 2 + 2 * n] = -yr2[n] * w[N / 2 - 2 * n - 1]; + out[N / 2 + 2 * n + 1] = yi2[N / 8 - n - 1] * w[N / 2 - 2 * n - 2]; + out[3 * N / 4 + 2 * n] = yi2[n] * w[N / 4 - 2 * n - 1]; + out[3 * N / 4 + 2 * n + 1] = -yr2[N / 8 - n - 1] * w[N / 4 - 2 * n - 2]; + } +} + +fn logadd(a: i32, b: i32) -> i32 { + let c = a - b; + let addr = ((c.abs() >> 1) as usize).min(255); + if c >= 0 { a + tables::LATAB[addr] as i32 } else { b + tables::LATAB[addr] as i32 } +} + +fn calc_lowcomp(a: i32, b0: i32, b1: i32, bin: usize) -> i32 { + if bin < 7 { + if b0 + 256 == b1 { + 384 + } else if b0 > b1 { + (a - 64).max(0) + } else { + a + } + } else if bin < 20 { + if b0 + 256 == b1 { + 320 + } else if b0 > b1 { + (a - 64).max(0) + } else { + a + } + } else { + (a - 128).max(0) + } +} + +/// Inputs to the bit allocation of one channel. +struct AllocParams { + start: usize, + end: usize, + fgain: i32, + snroffset: i32, + /// Coupling channel leak initialisation. + leak: Option<(i32, i32)>, + + /// Delta bit allocation segments (offset, length, ba) when in use. + dba: Option>, +} + +/// The parametric bit allocation of §7.2.2 for one exponent set. +fn bit_allocation(exp: &[u8; 256], p: &AllocParams, g: &Globals, bap: &mut [u8; 256]) { + let (start, end) = (p.start, p.end); + if start >= end { + return; + } + let (sdecay, fdecay, sgain, dbknee, floor) = (SLOWDEC[g.sdcycod], FASTDEC[g.fdcycod], SLOWGAIN[g.sgaincod], DBPBTAB[g.dbpbcod], FLOORTAB[g.floorcod]); + let mut psd = [0i32; 256]; + for bin in start..end { + psd[bin] = 3072 - ((exp[bin] as i32) << 7); + } + // PSD integration + let mut bndpsd = [0i32; 50]; + let mut j = start; + let mut k = tables::MASKTAB[start] as usize; + loop { + let lastbin = (tables::BNDTAB[k] as usize + tables::BNDSZ[k] as usize).min(end); + bndpsd[k] = psd[j]; + j += 1; + while j < lastbin { + bndpsd[k] = logadd(bndpsd[k], psd[j]); + j += 1; + } + k += 1; + if end <= lastbin || k >= 50 { + break; + } + } + // excitation + let bndstrt = tables::MASKTAB[start] as usize; + let bndend = tables::MASKTAB[end - 1] as usize + 1; + let mut excite = [0i32; 50]; + let (mut fastleak, mut slowleak) = p.leak.unwrap_or((0, 0)); + let begin; + if bndstrt == 0 { + // the last band of the LFE channel (bndend 7) has no band above it + let lfe_last = |bin: usize| !(bndend == 7 && bin == 6); + let mut lowcomp = 0; + lowcomp = calc_lowcomp(lowcomp, bndpsd[0], bndpsd[1], 0); + excite[0] = bndpsd[0] - p.fgain - lowcomp; + lowcomp = calc_lowcomp(lowcomp, bndpsd[1], bndpsd[2], 1); + excite[1] = bndpsd[1] - p.fgain - lowcomp; + let mut b = 7; + for bin in 2..7 { + if lfe_last(bin) { + lowcomp = calc_lowcomp(lowcomp, bndpsd[bin], bndpsd[bin + 1], bin); + } + fastleak = bndpsd[bin] - p.fgain; + slowleak = bndpsd[bin] - sgain; + excite[bin] = fastleak - lowcomp; + if lfe_last(bin) && bndpsd[bin] <= bndpsd[bin + 1] { + b = bin + 1; + break; + } + } + for bin in b..bndend.min(22) { + if lfe_last(bin) { + lowcomp = calc_lowcomp(lowcomp, bndpsd[bin], bndpsd[bin + 1], bin); + } + fastleak -= fdecay; + fastleak = fastleak.max(bndpsd[bin] - p.fgain); + slowleak -= sdecay; + slowleak = slowleak.max(bndpsd[bin] - sgain); + excite[bin] = (fastleak - lowcomp).max(slowleak); + } + begin = 22; + } else { + begin = bndstrt; + } + for bin in begin..bndend { + fastleak -= fdecay; + fastleak = fastleak.max(bndpsd[bin] - p.fgain); + slowleak -= sdecay; + slowleak = slowleak.max(bndpsd[bin] - sgain); + excite[bin] = fastleak.max(slowleak); + } + // masking curve + let mut mask = [0i32; 50]; + for bin in bndstrt..bndend { + if bndpsd[bin] < dbknee { + excite[bin] += (dbknee - bndpsd[bin]) >> 2; + } + mask[bin] = excite[bin].max(tables::HTH[g.fscod][bin] as i32); + } + // delta bit allocation + if let Some(segs) = &p.dba { + let mut band = 0usize; + for &(off, len, ba) in segs { + band += off as usize; + let delta = if ba >= 4 { (ba as i32 - 3) << 7 } else { (ba as i32 - 4) << 7 }; + for _ in 0..len { + if band < 50 { + mask[band] += delta; + } + band += 1; + } + } + } + // bit allocation pointers + let mut i = start; + let mut j = tables::MASKTAB[start] as usize; + loop { + let lastbin = (tables::BNDTAB[j] as usize + tables::BNDSZ[j] as usize).min(end); + let mut m = mask[j] - p.snroffset - floor; + if m < 0 { + m = 0; + } + m &= 0x1fe0; + m += floor; + while i < lastbin { + let addr = ((psd[i] - m) >> 5).clamp(0, 63) as usize; + bap[i] = tables::BAPTAB[addr]; + i += 1; + } + j += 1; + if end <= lastbin || j >= 50 { + break; + } + } +} + +/// Bit allocation parameters common to all channels of a block. +#[derive(Clone, Copy, Default)] +struct Globals { + fscod: usize, + sdcycod: usize, + fdcycod: usize, + sgaincod: usize, + dbpbcod: usize, + floorcod: usize, +} + +/// State carried from block to block within a syncframe. +struct Frame { + nfchans: usize, + blksw: [bool; 5], + dithflag: [bool; 5], + dynrng: [f32; 2], + cplinu: bool, + chincpl: [bool; 5], + phsflginu: bool, + cplbegf: usize, + cplendf: usize, + /// Coupling band of each coupling sub-band (relative to cplbegf). + sub_to_band: [usize; 18], + ncplbnd: usize, + cplco: [[f32; 18]; 5], + phsflg: [bool; 18], + rematflg: [bool; 4], + expstr: [u8; 7], + endmant: [usize; 5], + exps: [[u8; 256]; 7], + g: Globals, + csnroffst: i32, + fsnroffst: [i32; 7], + fgaincod: [usize; 7], + cplleak: (i32, i32), + deltbae: [u8; 7], + deltsegs: [Vec<(u8, u8, u8)>; 7], + bap: [[u8; 256]; 7], +} + +impl Frame { + fn new(nfchans: usize, fscod: usize) -> Frame { + Frame { + nfchans, + blksw: [false; 5], + dithflag: [false; 5], + dynrng: [1.0; 2], + cplinu: false, + chincpl: [false; 5], + phsflginu: false, + cplbegf: 0, + cplendf: 0, + sub_to_band: [0; 18], + ncplbnd: 0, + cplco: [[0.0; 18]; 5], + phsflg: [false; 18], + rematflg: [false; 4], + expstr: [0; 7], + endmant: [0; 5], + exps: [[0; 256]; 7], + g: Globals { fscod, ..Default::default() }, + csnroffst: 0, + fsnroffst: [0; 7], + fgaincod: [0; 7], + cplleak: (0, 0), + deltbae: [2; 7], + deltsegs: Default::default(), + bap: [[0; 256]; 7], + } + } + fn cplstrtmant(&self) -> usize { + 37 + 12 * self.cplbegf + } + fn cplendmant(&self) -> usize { + 37 + 12 * (self.cplendf + 3) + } +} + +/// dynrng gain (§7.7.1.2): X (3-bit signed) gives 6.02 dB steps, Y a linear factor in [1/2, 1). +fn dynrng_gain(v: u32) -> f32 { + let x = ((v as i32) << 24) >> 29; + let y = (v & 0x1F) as f32; + 2f32.powi(x + 1) * (32.0 + y) / 64.0 +} + +/// Decode exponents of one set (§7.1.3) into `exp[first..]`; `ngrps` groups of 3 mapped values, +/// each expanded to `grpsize` bins after the absolute exponent. +fn decode_exponents(r: &mut Bits, absexp: i32, ngrps: usize, grpsize: usize, exp: &mut [u8; 256], first: usize, skip_abs: bool) -> Result<()> { + let mut prev = absexp; + let mut at = first; + if !skip_abs { + exp[at] = absexp.clamp(0, 24) as u8; + at += 1; + } + for _ in 0..ngrps { + let g = r.read(7) as i32; + if g > 124 { + return Err(Error::Invalid("exponent group out of range")); + } + for d in [g / 25, (g % 25) / 5, g % 5] { + prev += d - 2; + if !(0..=24).contains(&prev) { + return Err(Error::Invalid("exponent out of range")); + } + for _ in 0..grpsize { + if at < 256 { + exp[at] = prev as u8; + } + at += 1; + } + } + } + Ok(()) +} + +/// Grouped-mantissa buffers shared by all exponent sets of a block (§7.3.5). +#[derive(Default)] +struct Groups { + b1: ([f32; 3], usize), + b2: ([f32; 3], usize), + b4: ([f32; 2], usize), +} + +const fn sym(levels: i32, code: i32) -> f32 { + (2 * code - (levels - 1)) as f32 / levels as f32 +} + +/// The mantissa value (before the exponent shift) of a bin with allocation `bap`. +fn mantissa(r: &mut Bits, bap: u8, gr: &mut Groups) -> f32 { + match bap { + 1 => { + if gr.b1.1 == 0 { + let g = r.read(5) as i32; + gr.b1.0 = [sym(3, g / 9), sym(3, (g % 9) / 3), sym(3, g % 3)]; + gr.b1.1 = 3; + } + gr.b1.1 -= 1; + gr.b1.0[2 - gr.b1.1] + } + 2 => { + if gr.b2.1 == 0 { + let g = r.read(7) as i32; + gr.b2.0 = [sym(5, g / 25), sym(5, (g % 25) / 5), sym(5, g % 5)]; + gr.b2.1 = 3; + } + gr.b2.1 -= 1; + gr.b2.0[2 - gr.b2.1] + } + 3 => sym(7, r.read(3) as i32), + 4 => { + if gr.b4.1 == 0 { + let g = r.read(7) as i32; + gr.b4.0 = [sym(11, g / 11), sym(11, g % 11)]; + gr.b4.1 = 2; + } + gr.b4.1 -= 1; + gr.b4.0[1 - gr.b4.1] + } + 5 => sym(15, r.read(4) as i32), + _ => { + let bits = QNTZTAB[bap as usize & 15]; + let v = r.read(bits) as i32; + let v = (v << (32 - bits)) >> (32 - bits); + v as f32 / (1u32 << (bits - 1)) as f32 + } + } +} + +/// A decoder: keeps the overlap of each channel between frames. +pub struct Decoder { + delay: [[f32; 256]; 6], + dither: u32, + seed: u32, +} + +impl Default for Decoder { + fn default() -> Self { + Self::new() + } +} + +/// A decoded syncframe: 1536 samples per channel, WAV channel order. +#[derive(Clone, Debug)] +pub struct Decoded { + pub header: Header, + pub channels: Vec>, +} + +impl Decoder { + pub fn new() -> Decoder { + Decoder { delay: [[0.0; 256]; 6], dither: 0, seed: 0x1234_5678 } + } + + /// Seed the dither generator (dither is decoder-specific noise; tests compare two seeds). + pub fn set_dither_seed(&mut self, seed: u32) { + self.seed = seed; + } + + /// Forget the overlap (before decoding from another position). + pub fn reset(&mut self) { + self.delay = [[0.0; 256]; 6]; + } + + fn dither(&mut self) -> f32 { + // uniform in [-0.707, 0.707] (§7.3.4; any reasonably random sequence) + self.dither = self.dither.wrapping_mul(1_664_525).wrapping_add(1_013_904_223); + ((self.dither >> 8) as f32 / (1u32 << 24) as f32 * 2.0 - 1.0) * 0.707 + } + + /// Decode one syncframe. + pub fn decode(&mut self, frame: &[u8]) -> Result { + let h = parse_header(frame)?; + if frame.len() < h.frame_bytes { + return Err(Error::Invalid("truncated syncframe")); + } + let fscod = (frame[4] >> 6) as usize; + // Dither is reseeded per syncframe from its CRC words (plus the configured seed), so a + // frame decodes to the same samples however it was reached (seeking, render caches). + let crcs = u32::from_be_bytes([frame[2], frame[3], frame[h.frame_bytes - 2], frame[h.frame_bytes - 1]]); + self.dither = self.seed ^ crcs.wrapping_mul(0x9E37_79B9); + let mut r = Bits::new(&frame[..h.frame_bytes]); + r.skip(40); // syncinfo + // bsi + r.skip(8); // bsid, bsmod + let acmod = r.read(3) as usize; + if acmod & 1 != 0 && acmod != 1 { + r.skip(2); + } + if acmod & 4 != 0 { + r.skip(2); + } + if acmod == 2 { + r.skip(2); + } + let lfeon = r.bit(); + for k in 0..if acmod == 0 { 2 } else { 1 } { + let _ = k; + r.skip(5); // dialnorm + if r.bit() { + r.skip(8); // compr + } + if r.bit() { + r.skip(8); // langcod + } + if r.bit() { + r.skip(7); // mixlevel, roomtyp + } + } + r.skip(2); // copyrightb, origbs + if r.bit() { + r.skip(14); + } + if r.bit() { + r.skip(14); + } + if r.bit() { + let l = r.read(6) as usize; + r.skip((l + 1) * 8); + } + let nfchans = NFCHANS[acmod]; + let mut f = Frame::new(nfchans, fscod); + let nout = nfchans + lfeon as usize; + let mut pcm: Vec> = vec![Vec::with_capacity(1536); nout]; + for blk in 0..6 { + let mut coefs = [[0f32; 256]; 6]; + self.audio_block(&mut r, &mut f, acmod, lfeon, blk, &mut coefs)?; + if r.overrun() { + return Err(Error::Invalid("audio block runs past the frame")); + } + // inverse transform, overlap-add + let mut buf = [0f32; 512]; + for ch in 0..nout { + let src = if ch < nfchans { ch } else { LFE }; + let short = ch < nfchans && f.blksw[ch]; + imdct(&coefs[src], short, &mut buf); + let d = &mut self.delay[src]; + for n in 0..256 { + pcm[ch].push(2.0 * (buf[n] + d[n])); + d[n] = buf[256 + n]; + } + } + } + Ok(Decoded { header: h, channels: wav_order(pcm, acmod, lfeon) }) + } + + fn audio_block(&mut self, r: &mut Bits, f: &mut Frame, acmod: usize, lfeon: bool, blk: usize, coefs: &mut [[f32; 256]; 6]) -> Result<()> { + let nf = f.nfchans; + for ch in 0..nf { + f.blksw[ch] = r.bit(); + } + for ch in 0..nf { + f.dithflag[ch] = r.bit(); + } + for k in 0..if acmod == 0 { 2 } else { 1 } { + if r.bit() { + f.dynrng[k] = dynrng_gain(r.read(8)); + } else if blk == 0 { + f.dynrng[k] = 1.0; + } + } + // coupling strategy + if r.bit() { + f.cplinu = r.bit(); + if f.cplinu { + for ch in 0..nf { + f.chincpl[ch] = r.bit(); + } + if acmod == 2 { + f.phsflginu = r.bit(); + } + f.cplbegf = r.read(4) as usize; + f.cplendf = r.read(4) as usize; + if f.cplendf + 3 <= f.cplbegf { + return Err(Error::Invalid("coupling end below start")); + } + let nsub = 3 + f.cplendf - f.cplbegf; + let mut band = 0; + f.sub_to_band[0] = 0; + for s in 1..nsub { + if !r.bit() { + band += 1; + } + f.sub_to_band[s] = band; + } + f.ncplbnd = band + 1; + } else { + f.chincpl = [false; 5]; + } + } else if blk == 0 { + return Err(Error::Invalid("no coupling strategy in block 0")); + } + // coupling coordinates + if f.cplinu { + let mut any = false; + for ch in 0..nf { + if f.chincpl[ch] && r.bit() { + any = true; + let mstr = r.read(2) as i32; + for bnd in 0..f.ncplbnd { + let e = r.read(4) as i32; + let m = r.read(4) as f32; + let v = if e == 15 { m / 16.0 } else { (m + 16.0) / 32.0 }; + let co = v * 2f32.powi(-(e + 3 * mstr)); + for s in 0..(3 + f.cplendf - f.cplbegf) { + if f.sub_to_band[s] == bnd { + f.cplco[ch][s] = co; + } + } + } + } + } + if acmod == 2 && f.phsflginu && any { + for bnd in 0..f.ncplbnd { + let flag = r.bit(); + for s in 0..(3 + f.cplendf - f.cplbegf) { + if f.sub_to_band[s] == bnd { + f.phsflg[s] = flag; + } + } + } + } + } + // rematrixing + if acmod == 2 && r.bit() { + let n = if !f.cplinu || f.cplbegf > 2 { + 4 + } else if f.cplbegf > 0 { + 3 + } else { + 2 + }; + f.rematflg = [false; 4]; + for k in 0..n { + f.rematflg[k] = r.bit(); + } + } + // exponent strategies + if f.cplinu { + f.expstr[CPL] = r.read(2) as u8; + } + for ch in 0..nf { + f.expstr[ch] = r.read(2) as u8; + } + if lfeon { + f.expstr[LFE] = r.read(1) as u8; + } + for ch in 0..nf { + if f.expstr[ch] != 0 { + if f.chincpl[ch] && f.cplinu { + f.endmant[ch] = f.cplstrtmant(); + } else { + let bw = r.read(6) as usize; + if bw > 60 { + return Err(Error::Invalid("chbwcod > 60")); + } + f.endmant[ch] = (bw + 12) * 3 + 37; + } + } + } + // exponents + if f.cplinu && f.expstr[CPL] != 0 { + let abs = (r.read(4) as i32) << 1; + let grpsize = 1 << (f.expstr[CPL] - 1); + let ngrps = (f.cplendmant() - f.cplstrtmant()) / (3 * grpsize); + let first = f.cplstrtmant(); + decode_exponents(r, abs, ngrps, grpsize, &mut f.exps[CPL], first, true)?; + } + for ch in 0..nf { + if f.expstr[ch] != 0 { + let abs = r.read(4) as i32; + let grpsize = 1usize << (f.expstr[ch] - 1); + let ngrps = (f.endmant[ch] - 1 + 3 * grpsize - 3) / (3 * grpsize); + decode_exponents(r, abs, ngrps, grpsize, &mut f.exps[ch], 0, false)?; + r.skip(2); // gainrng + } + } + if lfeon && f.expstr[LFE] != 0 { + let abs = r.read(4) as i32; + decode_exponents(r, abs, 2, 1, &mut f.exps[LFE], 0, false)?; + } + // bit allocation parametric information + if r.bit() { + f.g.sdcycod = r.read(2) as usize; + f.g.fdcycod = r.read(2) as usize; + f.g.sgaincod = r.read(2) as usize; + f.g.dbpbcod = r.read(2) as usize; + f.g.floorcod = r.read(3) as usize; + } else if blk == 0 { + return Err(Error::Invalid("no bit allocation information in block 0")); + } + if r.bit() { + f.csnroffst = r.read(6) as i32; + if f.cplinu { + f.fsnroffst[CPL] = r.read(4) as i32; + f.fgaincod[CPL] = r.read(3) as usize; + } + for ch in 0..nf { + f.fsnroffst[ch] = r.read(4) as i32; + f.fgaincod[ch] = r.read(3) as usize; + } + if lfeon { + f.fsnroffst[LFE] = r.read(4) as i32; + f.fgaincod[LFE] = r.read(3) as usize; + } + } else if blk == 0 { + return Err(Error::Invalid("no SNR offsets in block 0")); + } + if f.cplinu && r.bit() { + f.cplleak = (((r.read(3) as i32) << 8) + 768, ((r.read(3) as i32) << 8) + 768); + } + // delta bit allocation + if blk == 0 { + f.deltbae = [2; 7]; + } + if r.bit() { + if f.cplinu { + f.deltbae[CPL] = r.read(2) as u8; + } + for ch in 0..nf { + f.deltbae[ch] = r.read(2) as u8; + } + let order: Vec = f.cplinu.then_some(CPL).into_iter().chain(0..nf).collect(); + for ch in order { + if f.deltbae[ch] == 1 { + let n = r.read(3) as usize + 1; + f.deltsegs[ch] = (0..n).map(|_| (r.read(5) as u8, r.read(4) as u8, r.read(3) as u8)).collect(); + } else if f.deltbae[ch] == 3 { + return Err(Error::Invalid("reserved delta bit allocation state")); + } + } + } + // skip field + if r.bit() { + let l = r.read(9) as usize; + r.skip(l * 8); + } + // bit allocation for every exponent set + let all_zero = f.csnroffst == 0 && (0..7).all(|c| f.fsnroffst[c] == 0); + let snr = |fine: i32| (((f.csnroffst - 15) << 4) + fine) << 2; + let mut sets: Vec<(usize, AllocParams)> = Vec::new(); + for ch in 0..nf { + let dba = (f.deltbae[ch] <= 1).then(|| f.deltsegs[ch].clone()); + sets.push((ch, AllocParams { start: 0, end: f.endmant[ch], fgain: FASTGAIN[f.fgaincod[ch]], snroffset: snr(f.fsnroffst[ch]), leak: None, dba })); + } + if f.cplinu { + let dba = (f.deltbae[CPL] <= 1).then(|| f.deltsegs[CPL].clone()); + sets.push(( + CPL, + AllocParams { + start: f.cplstrtmant(), + end: f.cplendmant(), + fgain: FASTGAIN[f.fgaincod[CPL]], + snroffset: snr(f.fsnroffst[CPL]), + leak: Some(f.cplleak), + dba, + }, + )); + } + if lfeon { + sets.push((LFE, AllocParams { start: 0, end: 7, fgain: FASTGAIN[f.fgaincod[LFE]], snroffset: snr(f.fsnroffst[LFE]), leak: None, dba: None })); + } + for (ch, p) in &sets { + f.bap[*ch] = [0; 256]; + if !all_zero { + let exps = f.exps[*ch]; + bit_allocation(&exps, p, &f.g, &mut f.bap[*ch]); + } + } + // mantissas + let mut gr = Groups::default(); + let mut cpl = [0f32; 256]; + let mut got_cpl = false; + for ch in 0..nf { + for bin in 0..f.endmant[ch] { + let b = f.bap[ch][bin]; + let m = if b == 0 { if f.dithflag[ch] { self.dither() } else { 0.0 } } else { mantissa(r, b, &mut gr) }; + coefs[ch][bin] = m * 2f32.powi(-(f.exps[ch][bin] as i32)); + } + if f.cplinu && f.chincpl[ch] && !got_cpl { + got_cpl = true; + for bin in f.cplstrtmant()..f.cplendmant() { + let b = f.bap[CPL][bin]; + // zero-bit coupling mantissas are dithered per channel when decoupling + cpl[bin] = if b == 0 { f32::NAN } else { mantissa(r, b, &mut gr) }; + } + } + } + if lfeon { + for bin in 0..7 { + let b = f.bap[LFE][bin]; + let m = if b == 0 { 0.0 } else { mantissa(r, b, &mut gr) }; + coefs[LFE][bin] = m * 2f32.powi(-(f.exps[LFE][bin] as i32)); + } + } + // decoupling + if f.cplinu { + for ch in 0..nf { + if !f.chincpl[ch] { + continue; + } + for bin in f.cplstrtmant()..f.cplendmant() { + let s = (bin - 37) / 12 - f.cplbegf; + let m = if cpl[bin].is_nan() { if f.dithflag[ch] { self.dither() } else { 0.0 } } else { cpl[bin] }; + let mut v = m * 2f32.powi(-(f.exps[CPL][bin] as i32)) * f.cplco[ch][s] * 8.0; + if acmod == 2 && ch == 1 && f.phsflginu && f.phsflg[s] { + v = -v; + } + coefs[ch][bin] = v; + } + } + } + // rematrixing + if acmod == 2 { + let end = if f.cplinu { f.cplstrtmant() } else { f.endmant[0].min(f.endmant[1]) }; + const B: [usize; 5] = [13, 25, 37, 61, 253]; + for k in 0..4 { + if f.rematflg[k] { + for bin in B[k]..B[k + 1].min(end) { + let (l, rr) = (coefs[0][bin], coefs[1][bin]); + coefs[0][bin] = l + rr; + coefs[1][bin] = l - rr; + } + } + } + } + // dynamic range + for ch in 0..nf { + let g = if acmod == 0 && ch == 1 { f.dynrng[1] } else { f.dynrng[0] }; + if g != 1.0 { + coefs[ch].iter_mut().for_each(|c| *c *= g); + } + } + if lfeon && f.dynrng[0] != 1.0 { + coefs[LFE].iter_mut().for_each(|c| *c *= f.dynrng[0]); + } + Ok(()) + } +} + +/// Coded channel order (Table 5.8) → WAV / SMPTE order: fronts L R C, LFE, surrounds. +fn wav_order(mut pcm: Vec>, acmod: usize, lfeon: bool) -> Vec> { + let nf = NFCHANS[acmod]; + let lfe = lfeon.then(|| pcm.pop().expect("lfe")); + let mut fbw = pcm; + // acmod 3, 5, 7 code L C R …: move C after R + if acmod & 1 == 1 && acmod != 1 { + let c = fbw.remove(1); + fbw.insert(2, c); + } + let fronts = match acmod { + 1 => 1, + 3 | 5 | 7 => 3, + _ => 2.min(nf), + }; + let mut out: Vec> = fbw.drain(..fronts).collect(); + if let Some(l) = lfe { + out.push(l); + } + out.extend(fbw); + out +} + +#[cfg(test)] +mod tests; diff --git a/crates/ac3/src/tables.rs b/crates/ac3/src/tables.rs new file mode 100644 index 000000000..29ec3cfb4 --- /dev/null +++ b/crates/ac3/src/tables.rs @@ -0,0 +1,69 @@ +// Generated from the tables of ATSC A/52:2012 (Tables 7.12-7.16, 7.33). +pub(crate) const MASKTAB: [u8; 256] = [ + 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 28, 28, 29, 29, 29, 30, 30, 30, 31, 31, 31, 32, + 32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 39, 40, + 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, + 43, 43, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 46, + 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, + 47, 47, 47, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 49, 49, + 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 49, 0, 0, 0, +]; + +pub(crate) const LATAB: [u16; 256] = [ + 64, 63, 62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, 52, 51, 50, 49, 48, 47, 47, 46, 45, 44, 44, 43, 42, 41, 41, 40, 39, 38, 38, 37, 36, 36, 35, 35, 34, 33, + 33, 32, 32, 31, 30, 30, 29, 29, 28, 28, 27, 27, 26, 26, 25, 25, 24, 24, 23, 23, 22, 22, 21, 21, 21, 20, 20, 19, 19, 19, 18, 18, 18, 17, 17, 17, 16, 16, 16, + 15, 15, 15, 14, 14, 14, 13, 13, 13, 13, 12, 12, 12, 12, 11, 11, 11, 11, 10, 10, 10, 10, 10, 9, 9, 9, 9, 9, 8, 8, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, + 6, 6, 6, 6, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, +]; + +pub(crate) const BAPTAB: [u8; 64] = [ + 0, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, + 13, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15, 15, 15, +]; + +pub(crate) const BNDTAB: [u8; 50] = [ + 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 31, 34, 37, 40, 43, 46, 49, 55, 61, 67, 73, 79, + 85, 97, 109, 121, 133, 157, 181, 205, 229, +]; + +pub(crate) const BNDSZ: [u8; 50] = [ + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6, 12, 12, 12, 12, 24, 24, 24, 24, + 24, +]; + +pub(crate) const HTH: [[u16; 50]; 3] = [ + [ + 1232, 1232, 1088, 1024, 992, 960, 944, 944, 928, 928, 928, 928, 928, 912, 912, 912, 896, 896, 880, 880, 864, 864, 848, 848, 832, 832, 816, 800, 784, + 768, 752, 752, 752, 752, 768, 784, 832, 912, 992, 1056, 1120, 1168, 1184, 1120, 1088, 1088, 1312, 2048, 2112, 2112, + ], + [ + 1264, 1264, 1120, 1040, 992, 976, 960, 944, 944, 928, 928, 928, 928, 928, 912, 912, 912, 896, 896, 896, 880, 880, 864, 864, 848, 848, 832, 832, 800, + 784, 768, 752, 752, 752, 752, 768, 800, 848, 912, 992, 1056, 1104, 1184, 1168, 1120, 1088, 1152, 1584, 2112, 2112, + ], + [ + 1408, 1408, 1200, 1104, 1056, 1008, 992, 976, 960, 944, 944, 944, 928, 928, 928, 928, 928, 928, 928, 928, 912, 912, 912, 912, 896, 896, 896, 880, 864, + 848, 832, 816, 800, 784, 768, 752, 752, 752, 768, 784, 816, 848, 960, 1040, 1136, 1184, 1120, 1088, 1104, 1248, + ], +]; +/// Table 7.33 as printed (5 decimals), for checking the computed window. +#[cfg(test)] +#[allow(clippy::approx_constant)] +pub(crate) const WINDOW_TABLE: [f32; 256] = [ + 0.00014, 0.00024, 0.00037, 0.00051, 0.00067, 0.00086, 0.00107, 0.0013, 0.00157, 0.00187, 0.0022, 0.00256, 0.00297, 0.00341, 0.0039, 0.00443, 0.00501, + 0.00564, 0.00632, 0.00706, 0.00785, 0.00871, 0.00962, 0.01061, 0.01166, 0.01279, 0.01399, 0.01526, 0.01662, 0.01806, 0.01959, 0.02121, 0.02292, 0.02472, + 0.02662, 0.02863, 0.03073, 0.03294, 0.03527, 0.0377, 0.04025, 0.04292, 0.04571, 0.04862, 0.05165, 0.05481, 0.0581, 0.06153, 0.06508, 0.06878, 0.07261, + 0.07658, 0.08069, 0.08495, 0.08935, 0.09389, 0.09859, 0.10343, 0.10842, 0.11356, 0.11885, 0.12429, 0.12988, 0.13563, 0.14152, 0.14757, 0.15376, 0.16011, + 0.16661, 0.17325, 0.18005, 0.18699, 0.19407, 0.2013, 0.20867, 0.21618, 0.22382, 0.23161, 0.23952, 0.24757, 0.25574, 0.26404, 0.27246, 0.281, 0.28965, + 0.29841, 0.30729, 0.31626, 0.32533, 0.3345, 0.34376, 0.35311, 0.36253, 0.37204, 0.38161, 0.39126, 0.40096, 0.41072, 0.42054, 0.4304, 0.4403, 0.45023, + 0.4602, 0.47019, 0.4802, 0.49022, 0.50025, 0.51028, 0.52031, 0.53033, 0.54033, 0.55031, 0.56026, 0.57019, 0.58007, 0.58991, 0.5997, 0.60944, 0.61912, + 0.62873, 0.63827, 0.64774, 0.65713, 0.66643, 0.67564, 0.68476, 0.69377, 0.70269, 0.7115, 0.72019, 0.72877, 0.73723, 0.74557, 0.75378, 0.76186, 0.76981, + 0.77762, 0.7853, 0.79283, 0.80022, 0.80747, 0.81457, 0.82151, 0.82831, 0.83496, 0.84145, 0.84779, 0.85398, 0.86001, 0.86588, 0.8716, 0.87716, 0.88257, + 0.88782, 0.89291, 0.89785, 0.90264, 0.90728, 0.91176, 0.9161, 0.92028, 0.92432, 0.92822, 0.93197, 0.93558, 0.93906, 0.9424, 0.9456, 0.94867, 0.95162, + 0.95444, 0.95713, 0.95971, 0.96217, 0.96451, 0.96674, 0.96887, 0.97089, 0.97281, 0.97463, 0.97635, 0.97799, 0.97953, 0.98099, 0.98236, 0.98366, 0.98488, + 0.98602, 0.9871, 0.98811, 0.98905, 0.98994, 0.99076, 0.99153, 0.99225, 0.99291, 0.99353, 0.99411, 0.99464, 0.99513, 0.99558, 0.996, 0.99639, 0.99674, + 0.99706, 0.99736, 0.99763, 0.99788, 0.99811, 0.99831, 0.9985, 0.99867, 0.99882, 0.99895, 0.99908, 0.99919, 0.99929, 0.99938, 0.99946, 0.99953, 0.99959, + 0.99965, 0.99969, 0.99974, 0.99978, 0.99981, 0.99984, 0.99986, 0.99988, 0.9999, 0.99992, 0.99993, 0.99994, 0.99995, 0.99996, 0.99997, 0.99998, 0.99998, + 0.99998, 0.99999, 0.99999, 0.99999, 0.99999, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, +]; diff --git a/crates/ac3/src/tests.rs b/crates/ac3/src/tests.rs new file mode 100644 index 000000000..d015acbe4 --- /dev/null +++ b/crates/ac3/src/tests.rs @@ -0,0 +1,98 @@ +use super::*; + +#[test] +fn computed_window_matches_table_7_33() { + let w = window(); + for (n, (&a, &b)) in w.iter().zip(tables::WINDOW_TABLE.iter()).enumerate() { + assert!((a - b).abs() <= 6e-6, "w[{n}] = {a}, table {b}"); + } + // princen-bradley: w[n]² + w[255-n]²... (the second half mirrors): w[n]² + w[511-n]² = 1 + for n in 0..256 { + let m = w[255 - n]; + assert!((w[n] * w[n] + m * m - 1.0).abs() < 1e-5); + } +} + +#[test] +fn tables_are_consistent() { + // masktab maps every bin of a band to that band + for b in 0..50 { + for bin in tables::BNDTAB[b] as usize..(tables::BNDTAB[b] + tables::BNDSZ[b]) as usize { + if bin < 253 { + assert_eq!(tables::MASKTAB[bin] as usize, b, "bin {bin}"); + } + } + } + assert_eq!(tables::BNDTAB[49] as usize + tables::BNDSZ[49] as usize, 253); + // latab is non-increasing from 64 to 0 + assert_eq!(tables::LATAB[0], 64); + assert!(tables::LATAB.windows(2).all(|w| w[0] >= w[1])); + assert!(tables::BAPTAB.windows(2).all(|w| w[0] <= w[1])); + assert_eq!((tables::BAPTAB[0], tables::BAPTAB[63]), (0, 15)); +} + +#[test] +fn inverse_fft_matches_the_definition() { + for n in [64usize, 128] { + let re0: Vec = (0..n).map(|i| ((i * 7 + 3) % 11) as f32 - 5.0).collect(); + let im0: Vec = (0..n).map(|i| ((i * 5 + 1) % 13) as f32 - 6.0).collect(); + let (mut re, mut im) = (re0.clone(), im0.clone()); + ifft(&mut re, &mut im); + for t in 0..n { + let (mut sr, mut si) = (0f64, 0f64); + for k in 0..n { + let a = 2.0 * std::f64::consts::PI * (k * t) as f64 / n as f64; + sr += re0[k] as f64 * a.cos() - im0[k] as f64 * a.sin(); + si += re0[k] as f64 * a.sin() + im0[k] as f64 * a.cos(); + } + assert!((sr - re[t] as f64).abs() < 1e-3 && (si - im[t] as f64).abs() < 1e-3, "n {n} t {t}"); + } + } +} + +#[test] +fn headers_and_gains() { + // 48 kHz, 448 kb/s, 3/2 + LFE: bsid 8, bsmod 0, acmod 7 (cmixlev, surmixlev), lfeon + let b = [0x0B, 0x77, 0, 0, 30, 8 << 3, 0b1110_0001, 0b0100_0000]; + let h = parse_header(&b).unwrap(); + assert_eq!((h.sample_rate, h.bitrate_kbps, h.frame_bytes, h.acmod, h.lfeon, h.channels()), (48_000, 448, 1792, 7, true, 6)); + assert_eq!(parse_header(&[0x0B, 0x77, 0, 0, 0, 16 << 3, 0, 0]), Err(Error::Unsupported("bsid 16 (E-AC-3)".into()))); + assert_eq!(dynrng_gain(0), 1.0); + assert!((dynrng_gain(0b1110_0000) - 0.5).abs() < 1e-7); + assert!((dynrng_gain(0b0111_1111) - 16.0 * 63.0 / 64.0).abs() < 1e-5); +} + +#[test] +fn channel_order() { + let ch = |n: usize| (0..n).map(|i| vec![i as f32]).collect::>(); + let ids = |v: Vec>| v.iter().map(|c| c[0] as usize).collect::>(); + // coded L C R Ls Rs LFE → L R C LFE Ls Rs + assert_eq!(ids(wav_order(ch(6), 7, true)), vec![0, 2, 1, 5, 3, 4]); + assert_eq!(ids(wav_order(ch(3), 2, true)), vec![0, 1, 2]); + assert_eq!(ids(wav_order(ch(1), 1, false)), vec![0]); +} + +#[test] +fn garbage_never_panics() { + let mut d = Decoder::new(); + let mut x = 0xACu32; + for len in [0usize, 5, 8, 100, 2000] { + for _ in 0..200 { + let mut f: Vec = (0..len) + .map(|_| { + x ^= x << 13; + x ^= x >> 17; + x ^= x << 5; + x as u8 + }) + .collect(); + if len >= 6 { + f[0] = 0x0B; + f[1] = 0x77; + f[4] &= 0x3F; + f[5] = (f[5] & 7) | (8 << 3); + } + let _ = d.decode(&f); + } + } +} diff --git a/crates/ac3/tests/oracle.rs b/crates/ac3/tests/oracle.rs new file mode 100644 index 000000000..9c2ea4222 --- /dev/null +++ b/crates/ac3/tests/oracle.rs @@ -0,0 +1,147 @@ +//! AC-3 against ffmpeg: ffmpeg's AC-3 encoder writes raw syncframes from generated tones and +//! noise; ffmpeg's decode of the same file is the reference. AC-3 output is not bit-exact across +//! decoders (zero-bit mantissas are filled with decoder-specific dither, §7.3.4; float vs fixed +//! arithmetic), so the criterion is SNR against ffmpeg's decode. + +use std::path::{Path, PathBuf}; +use std::process::{Command, Stdio}; + +struct Spec { + name: &'static str, + input: &'static str, + rate: u32, + channels: u32, + enc: &'static [&'static str], + /// Minimum SNR (dB) against ffmpeg. + min_snr: f64, +} + +const SPECS: &[Spec] = &[ + Spec { + name: "ac3_stereo_192k", + input: "aevalsrc=exprs=0.3*sin(2*PI*441*t)+0.1*sin(2*PI*3100*t)|0.25*sin(2*PI*660*t)+0.1*sin(2*PI*(500+4000*t)*t):s=48000", + rate: 48_000, + channels: 2, + enc: &["-b:a", "192k"], + min_snr: 60.0, + }, + // low bitrate stereo: coupling and rematrixing + Spec { + name: "ac3_stereo_96k_coupling", + input: "aevalsrc=exprs=0.3*sin(2*PI*441*t)+0.05*sin(2*PI*9100*t)|0.3*sin(2*PI*441*t)+0.05*sin(2*PI*11000*t):s=44100", + rate: 44_100, + channels: 2, + enc: &["-b:a", "96k", "-channel_coupling", "1"], + min_snr: 40.0, + }, + Spec { name: "ac3_mono_32k", input: "sine=frequency=1000:sample_rate=32000", rate: 32_000, channels: 1, enc: &["-b:a", "96k"], min_snr: 60.0 }, + // 5.1 with LFE + Spec { + name: "ac3_51_448k", + input: "aevalsrc=exprs=0.3*sin(2*PI*300*t)|0.3*sin(2*PI*400*t)|0.3*sin(2*PI*500*t)|0.3*sin(2*PI*60*t)|0.2*sin(2*PI*700*t)|0.2*sin(2*PI*800*t):s=48000:c=5.1", + rate: 48_000, + channels: 6, + enc: &["-b:a", "448k"], + min_snr: 55.0, + }, + // noise (wideband, many zero-bit mantissas): the SNR is limited by dither + Spec { + name: "ac3_noise_stereo", + input: "anoisesrc=color=pink:sample_rate=48000:amplitude=0.3:seed=3,pan=stereo|c0=c0|c1=-0.5*c0", + rate: 48_000, + channels: 2, + enc: &["-b:a", "128k"], + min_snr: 18.0, + }, +]; + +fn make(ff: &Path, s: &Spec) -> Option { + let out = filmcraft_testkit::fixtures_dir("ac3").join(format!("{}.ac3", s.name)); + filmcraft_testkit::fixtures::generate(&out, |tmp| { + let st = Command::new(ff) + .args(["-hide_banner", "-v", "error", "-y", "-f", "lavfi", "-i", s.input, "-t", "1"]) + .args(["-c:a", "ac3"]) + .args(s.enc) + .args(["-f", "ac3"]) + .arg(tmp) + .stdin(Stdio::null()) + .status(); + matches!(st, Ok(x) if x.success()) + }) +} + +fn frames(data: &[u8]) -> Vec<&[u8]> { + let mut out = Vec::new(); + let mut p = 0; + while let Ok(h) = filmcraft_ac3::parse_header(&data[p..]) { + if p + h.frame_bytes > data.len() { + break; + } + out.push(&data[p..p + h.frame_bytes]); + p += h.frame_bytes; + } + out +} + +#[test] +fn ffmpeg_encoded_streams_decode() { + let ff = filmcraft_testkit::require_ffmpeg!(); + for s in SPECS { + let f = make(&ff, s).unwrap_or_else(|| panic!("could not generate {}", s.name)); + let data = std::fs::read(&f).unwrap(); + let mut dec = filmcraft_ac3::Decoder::new(); + let mut ours: Vec> = vec![Vec::new(); s.channels as usize]; + for fr in frames(&data) { + let d = dec.decode(fr).unwrap(); + assert_eq!((d.header.sample_rate, d.channels.len() as u32), (s.rate, s.channels), "{}", s.name); + for (o, c) in ours.iter_mut().zip(d.channels) { + o.extend(c); + } + } + let o = Command::new(&ff).args(["-v", "error", "-i"]).arg(&f).args(["-f", "f32le", "-"]).output().unwrap(); + let want: Vec = o.stdout.chunks_exact(4).map(|b| f32::from_le_bytes([b[0], b[1], b[2], b[3]])).collect(); + let ch = s.channels as usize; + let n = want.len() / ch; + assert_eq!(n, ours[0].len(), "{}: sample count", s.name); + let mut worst = f64::INFINITY; + for c in 0..ch { + let (mut sig, mut err) = (0f64, 0f64); + for i in 0..n { + let (a, b) = (ours[c][i] as f64, want[i * ch + c] as f64); + sig += b * b; + err += (a - b) * (a - b); + } + let snr = 10.0 * (sig / err.max(1e-30)).log10(); + worst = worst.min(snr); + } + // the same stream decoded with another dither sequence: how far dither alone moves the + // output (the difference to ffmpeg must be of that order) + let mut dec2 = filmcraft_ac3::Decoder::new(); + dec2.set_dither_seed(0xDEAD_BEEF); + let other: Vec = frames(&data).iter().flat_map(|fr| dec2.decode(fr).unwrap().channels.swap_remove(0)).collect(); + let (mut sig, mut err) = (0f64, 0f64); + for i in 0..n { + sig += (ours[0][i] as f64).powi(2); + err += (ours[0][i] as f64 - other[i] as f64).powi(2); + } + let dither_snr = 10.0 * (sig / err.max(1e-30)).log10(); + println!( + "{}: {} Hz, {} ch, {} frames, worst channel SNR vs ffmpeg {worst:.1} dB (two dither seeds: {dither_snr:.1} dB)", + s.name, + s.rate, + ch, + frames(&data).len() + ); + assert!(worst >= s.min_snr, "{}: SNR {worst:.1} dB < {}", s.name, s.min_snr); + } +} + +#[test] +#[ignore] +fn generate_fixtures() { + let Some(ff) = filmcraft_testkit::ffmpeg_or_skip("ac3 fixtures") else { return }; + for s in SPECS { + let out = filmcraft_testkit::fixtures_dir("ac3").join(format!("{}.ac3", s.name)); + filmcraft_testkit::fixtures::generate_and_report(&format!("ac3/{}", s.name), &[out], || make(&ff, s)); + } +} diff --git a/crates/codecs/Cargo.toml b/crates/codecs/Cargo.toml index 432ddb0b1..90e1c9ca8 100644 --- a/crates/codecs/Cargo.toml +++ b/crates/codecs/Cargo.toml @@ -19,6 +19,7 @@ filmcraft-matroska = { workspace = true } filmcraft-mxf = { workspace = true } filmcraft-mpegts = { workspace = true } filmcraft-mpeg2v = { workspace = true } +filmcraft-ac3 = { workspace = true } filmcraft-ogg = { workspace = true } filmcraft-bitstream = { workspace = true } filmcraft-h264 = { workspace = true } diff --git a/crates/codecs/src/audio.rs b/crates/codecs/src/audio.rs index 53fa12627..3cf7f400f 100644 --- a/crates/codecs/src/audio.rs +++ b/crates/codecs/src/audio.rs @@ -24,6 +24,8 @@ enum Inner { Opus { dec: Box, order: Option<&'static [usize]> }, /// Bootstrap decoders (MP3, ALAC, FLAC, HE-AAC…) via symphonia. Symphonia(Box), + /// Our AC-3 decoder (ATSC A/52). + Ac3(Box), } /// Opus always decodes at 48 kHz (RFC 7845 §5.1: the input rate in the header is informational). @@ -122,7 +124,7 @@ impl PacketDecoder { } /// AC-3 (ATSC A/52). pub fn ac3() -> Result { - Err(CodecError::Unsupported("AC-3 audio".into())) + Ok(Self { inner: Inner::Ac3(Box::new(filmcraft_ac3::Decoder::new())), channels: 0 }) } /// Whether this decodes Opus (which needs [`OPUS_PRE_ROLL`] of pre-roll after a seek). pub fn is_opus(&self) -> bool { @@ -133,6 +135,7 @@ impl PacketDecoder { match c { C::Aac(a) => Self::aac(&a.asc, if a.sample_rate > 0 { a.sample_rate } else { rate }), C::Mp3 => Self::new(CODEC_TYPE_MP3, rate, None), + C::Ac3 { .. } => Self::ac3(), C::Alac { cookie } => Self::new(CODEC_TYPE_ALAC, rate, Some(cookie.clone())), C::Flac(_) => Self::new(CODEC_TYPE_FLAC, rate, None), C::Opus(o) => Self::opus(filmcraft_opus::OpusHead::from_dops(&o.to_bytes()).map_err(|e| CodecError::Unsupported(format!("Opus: {e}")))?), @@ -161,6 +164,24 @@ impl PacketDecoder { } Inner::Symphonia(d) => d, Inner::AacPending => return Err(CodecError::Decode("AAC: no configuration yet".into())), + Inner::Ac3(d) => { + // a packet may hold several syncframes (MP4 / Matroska) + let mut out: Vec> = Vec::new(); + let mut p = 0; + while p < data.len() { + let f = d.decode(&data[p..]).map_err(|e| CodecError::Decode(e.to_string()))?; + p += f.header.frame_bytes; + if out.is_empty() { + out = f.channels; + } else { + for (o, c) in out.iter_mut().zip(f.channels) { + o.extend(c); + } + } + } + self.channels = out.len(); + return Ok(out); + } }; let pkt = Packet::new_from_slice(0, ts, 0, data); let buf = dec.decode(&pkt).map_err(|e| CodecError::Decode(e.to_string()))?; @@ -188,13 +209,14 @@ impl PacketDecoder { } Inner::Opus { dec, .. } => dec.reset(), Inner::Symphonia(d) => d.reset(), + Inner::Ac3(d) => d.reset(), Inner::AacPending => {} } } } /// Whether AC-3 audio decodes (our ATSC A/52 decoder). -pub const AC3_DECODER: bool = false; +pub const AC3_DECODER: bool = true; /// An AudioSpecificConfig for AAC with these parameters (ISO/IEC 14496-3 §1.6.2.1, GA specific /// config with no extension flags). diff --git a/crates/codecs/src/mkv.rs b/crates/codecs/src/mkv.rs index ff817a343..08196e0b8 100644 --- a/crates/codecs/src/mkv.rs +++ b/crates/codecs/src/mkv.rs @@ -335,6 +335,7 @@ impl MkvSource { Codec::Opus { head } => PacketDecoder::opus(filmcraft_opus::OpusHead::parse(head).map_err(|e| CodecError::Unsupported(format!("Opus: {e}")))?), Codec::Mp3 => PacketDecoder::new(CODEC_TYPE_MP3, rate, None), Codec::Mp2 => PacketDecoder::mpeg_audio(2, rate), + Codec::Ac3 => PacketDecoder::ac3(), // symphonia wants the STREAMINFO block body: skip `fLaC` + the 4-byte block header. Codec::Flac { private } => PacketDecoder::new(CODEC_TYPE_FLAC, rate, private.get(8..42).map(<[u8]>::to_vec)), Codec::Vorbis { headers } if headers.len() == 3 => { diff --git a/crates/codecs/tests/mpeg_oracle.rs b/crates/codecs/tests/mpeg_oracle.rs index 2f01d9d3d..5cdc7b92c 100644 --- a/crates/codecs/tests/mpeg_oracle.rs +++ b/crates/codecs/tests/mpeg_oracle.rs @@ -350,7 +350,8 @@ fn avchd_ac3_audio() { let a = src.info().audio.clone().unwrap(); if filmcraft_codecs::audio::AC3_DECODER { assert_eq!((a.codec.as_str(), a.channels), ("AC-3", 2)); - let e = check_audio(&ff, &f, 1e-3); + // AC-3 zero-bit mantissas carry decoder-specific dither (A/52 §7.3.4) + let e = check_audio(&ff, &f, 5e-3); println!("m2ts_h264_ac3: AC-3 max error {e:.2e}"); } else { assert_eq!(a.codec, "AC-3 (unsupported)"); @@ -396,7 +397,9 @@ fn program_streams() { // AC-3 in private stream 1 let f = make(&ff, "ps_mpeg2_ac3.vob"); let src = open(&f); - assert!(src.info().audio.as_ref().unwrap().codec.starts_with("AC-3")); + assert_eq!(src.info().audio.as_ref().unwrap().codec, "AC-3"); + let e = check_audio(&ff, &f, 5e-3); + println!("ps_mpeg2_ac3: AC-3 max error {e:.2e}"); // MPEG-1 system stream let f = make(&ff, "mpeg1_system.mpg"); let src = open(&f);