pauljd/craft-libs
1
0
Fork 0
mirror of https://github.com/storytold/craft-libs.git synced 2026-11-12 08:52:32 +00:00
Code Issues Packages Projects Releases Wiki Activity

craft-audio, craft-media: the changes FilmCraft and SoundCraft made to their copies since extraction (#807)

Browse Source 
* craft-audio: FilmCraft's FLAC encoder and MP1/MP2 decoder speed-ups

FLAC: per-order predictor kernels, merged Rice partition sums, the LPC
window kept per block length, and batches consumed by offset
(filmcraft 5eaff794, d9dad4a5). The output is unchanged: 64 golden cases
hash the same before and after.

MP1/MP2: the faster synthesis filterbank (filmcraft 51c755e5), with its
bit-exact test.

* craft-audio: FilmCraft's AAC and Opus decoder speed-ups and refactors

AAC: table-driven Huffman decoding, per-stage FFT twiddles and fused
stages, reused per-channel work buffers, and the encoder's memoised
per-band work (filmcraft f9627280, d587ceed, 6121635c, dc124ab7,
7a532eca).

Opus: radix-2/3/4 FFT kernels in the CELT IMDCT, the PVQ table built
once, one de-emphasis loop for both channels, the SILK resampler
without per-sample divisions and with eight partial sums, concealment
and FEC capped at 120 ms (filmcraft bf7ac420, 12e32026, c4dc64d1,
eb352ec8, 87b8e22d, 014b5e56, 78d77b72).

The per-frame mid-only flag this crate fixed in the SILK decoder is
kept in the refactored decode.

* craft-media: Ogg CRC-32 by slice-by-16

FilmCraft's 6763f9bf: sixteen independent lookups per 16 input bytes, so
opening a file is several times faster. Checked against a bit-at-a-time
reference on every short length and alignment, and across every split.

* craft-audio: SoundCraft's CAF cookie-size and FLAC vendor-length checks

A CAF magic cookie size above the signed range is refused instead of
passed to symphonia (soundcraft 3bd3965), and a Vorbis comment block
whose vendor string claims more than the block holds is refused
(soundcraft ddfe991). Each with its regression test.

* craft-audio: AC-3 reads bits through the shared BitReader

FilmCraft's 1c9c2e82: the AC-3 decoder's private bit-at-a-time reader is
replaced by craft-media's BitReader, with its two tests: the decoded
output of 2,000 synthetic and 12 real syncframes (and their damaged
variants) hashes the same as before, and the reader agrees with the old
one on random reads, skips and overruns. The three AC-3 features now
turn on craft-media's bitstream module.

* craft-audio: the quantiser rounds without a libm call

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