diff --git a/Cargo.lock b/Cargo.lock index 1d14bab9..866ae15f 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1221,6 +1221,7 @@ version = "0.1.0" dependencies = [ "effectcraft-color", "effectcraft-media", + "effectcraft-opusenc", "effectcraft-project", "effectcraft-raster", "effectcraft-render", @@ -1228,6 +1229,7 @@ dependencies = [ "filmcraft-aac", "filmcraft-h264enc", "filmcraft-isobmff", + "filmcraft-opus", "filmcraft-prores", "gif", "image", @@ -1347,6 +1349,10 @@ dependencies = [ "thiserror 2.0.21", ] +[[package]] +name = "effectcraft-opusenc" +version = "0.1.0" + [[package]] name = "effectcraft-path" version = "0.1.0" diff --git a/Cargo.toml b/Cargo.toml index 7aabc4af..08a7d19d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -34,6 +34,7 @@ effectcraft-engine = { path = "crates/engine" } effectcraft-host = { path = "crates/host" } effectcraft-ui-egui = { path = "crates/ui-egui" } effectcraft-automation = { path = "crates/automation" } +effectcraft-opusenc = { path = "crates/opusenc" } # foundation thiserror = "2" @@ -68,6 +69,8 @@ filmcraft-h264enc = { git = "https://github.com/storytold/filmcraft", rev = "89c filmcraft-prores = { git = "https://github.com/storytold/filmcraft", rev = "89cfbf24f698b0f0e538eb7795929f19bbfb1dbb" } filmcraft-aac = { git = "https://github.com/storytold/filmcraft", rev = "89cfbf24f698b0f0e538eb7795929f19bbfb1dbb" } filmcraft-isobmff = { git = "https://github.com/storytold/filmcraft", rev = "89cfbf24f698b0f0e538eb7795929f19bbfb1dbb" } +# decoder used as the round-trip oracle of effectcraft-opusenc (crates/export dev-dependency) +filmcraft-opus = { git = "https://github.com/storytold/filmcraft", rev = "89cfbf24f698b0f0e538eb7795929f19bbfb1dbb" } # gpu (the same wgpu eframe / egui-wgpu use) wgpu = { version = "30.0.1", features = ["fragile-send-sync-non-atomic-wasm"] } @@ -131,6 +134,9 @@ opt-level = 3 [profile.dev.package.effectcraft-media] opt-level = 3 +[profile.dev.package.effectcraft-opusenc] +opt-level = 3 + [profile.dev.package.effectcraft-export] opt-level = 3 diff --git a/crates/export/Cargo.toml b/crates/export/Cargo.toml index 34c83aab..5bdd0fee 100644 --- a/crates/export/Cargo.toml +++ b/crates/export/Cargo.toml @@ -18,6 +18,7 @@ filmcraft-h264enc = { workspace = true } filmcraft-prores = { workspace = true } filmcraft-aac = { workspace = true } filmcraft-isobmff = { workspace = true } +effectcraft-opusenc = { workspace = true } image = { workspace = true, features = ["png", "jpeg", "tiff", "exr"] } gif = { workspace = true } rayon = { workspace = true } @@ -28,3 +29,4 @@ log = { workspace = true } [dev-dependencies] effectcraft-media = { workspace = true } effectcraft-color = { workspace = true } +filmcraft-opus = { workspace = true } diff --git a/crates/export/tests/opus_roundtrip.rs b/crates/export/tests/opus_roundtrip.rs new file mode 100644 index 00000000..5ea099f0 --- /dev/null +++ b/crates/export/tests/opus_roundtrip.rs @@ -0,0 +1,364 @@ +//! Round trip of `effectcraft-opusenc` through FilmCraft's spec-derived Opus decoder, plus an +//! ffmpeg oracle (skipped when ffmpeg is not installed) reading an Ogg Opus file written here. + +use std::process::Command; + +use effectcraft_opusenc::OpusEncoder; + +const SR: f32 = 48_000.0; +const FRAME: usize = OpusEncoder::FRAME_SIZE; + +/// Encodes interleaved `pcm`, padding the tail so the encoder's look-ahead is flushed. +fn encode(pcm: &[f32], channels: u8, bitrate: u32) -> (OpusEncoder, Vec>) { + let c = channels as usize; + let mut enc = OpusEncoder::new(channels, bitrate); + let total = pcm.len() / c + enc.pre_skip() as usize; + let frames = total.div_ceil(FRAME); + let mut padded = pcm.to_vec(); + padded.resize(frames * FRAME * c, 0.0); + let packets = padded.chunks(FRAME * c).map(|f| enc.encode_float(f)).collect(); + (enc, packets) +} + +/// Decodes with filmcraft-opus, drops the pre-skip and returns planar channels of `len` samples. +fn decode(head: &[u8], packets: &[Vec], len: usize) -> Vec> { + let mut dec = filmcraft_opus::Decoder::new(head).expect("OpusHead"); + dec.set_trim_pre_skip(true); + let mut out: Vec> = vec![Vec::new(); dec.channels()]; + for p in packets { + let planes = dec.decode(Some(p)).expect("decode"); + for (o, p) in out.iter_mut().zip(planes) { + o.extend_from_slice(&p); + } + } + for o in &mut out { + assert!(o.len() >= len, "decoded {} < {len}", o.len()); + o.truncate(len); + } + out +} + +fn roundtrip(planar: &[Vec], bitrate: u32) -> (Vec>, Vec>) { + let c = planar.len(); + let len = planar[0].len(); + let mut pcm = vec![0f32; len * c]; + for (ch, p) in planar.iter().enumerate() { + for (i, &v) in p.iter().enumerate() { + pcm[i * c + ch] = v; + } + } + let (enc, packets) = encode(&pcm, c as u8, bitrate); + let decoded = decode(&enc.opus_head(), &packets, len); + (packets, decoded) +} + +/// SNR in dB over the steady part (the first and last 20 ms are excluded). +fn snr(orig: &[f32], dec: &[f32]) -> f64 { + let (a, b) = (FRAME, orig.len() - FRAME); + let mut s = 0f64; + let mut e = 0f64; + for i in a..b { + s += (orig[i] as f64).powi(2); + e += (orig[i] as f64 - dec[i] as f64).powi(2); + } + 10.0 * (s / e.max(1e-30)).log10() +} + +fn sine(freq: f32, amp: f32, len: usize) -> Vec { + (0..len).map(|i| amp * (2.0 * std::f32::consts::PI * freq * i as f32 / SR).sin()).collect() +} + +fn chirp(f0: f32, f1: f32, amp: f32, len: usize) -> Vec { + let dur = len as f32 / SR; + let k = (f1 - f0) / dur; + (0..len) + .map(|i| { + let t = i as f32 / SR; + amp * (2.0 * std::f32::consts::PI * (f0 * t + 0.5 * k * t * t)).sin() + }) + .collect() +} + +fn noise(amp: f32, len: usize, seed: u32) -> Vec { + let mut s = seed; + (0..len) + .map(|_| { + s = s.wrapping_mul(1_664_525).wrapping_add(1_013_904_223); + amp * ((s >> 8) as f32 / (1u32 << 24) as f32 * 2.0 - 1.0) + }) + .collect() +} + +/// Energy (dB) of `x` in each of a few octave bands up to 16 kHz, from a Hann-windowed FFT +/// averaged over 2048-sample blocks. +fn band_energies(x: &[f32]) -> Vec { + const N: usize = 2048; + let edges = [100.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0, 8000.0, 16000.0]; + let mut acc = vec![0f64; edges.len() - 1]; + let mut start = N; + while start + N <= x.len() - N { + let mut re: Vec = (0..N).map(|i| x[start + i] as f64 * (0.5 - 0.5 * (2.0 * std::f64::consts::PI * i as f64 / N as f64).cos())).collect(); + let mut im = vec![0f64; N]; + fft(&mut re, &mut im); + for k in 1..N / 2 { + let f = k as f64 * SR as f64 / N as f64; + if let Some(b) = edges.windows(2).position(|w| f >= w[0] && f < w[1]) { + acc[b] += re[k] * re[k] + im[k] * im[k]; + } + } + start += N; + } + acc.iter().map(|e| 10.0 * (e + 1e-20).log10()).collect() +} + +fn fft(re: &mut [f64], im: &mut [f64]) { + let n = re.len(); + let mut j = 0; + for i in 1..n { + let mut bit = n >> 1; + while j & bit != 0 { + j ^= bit; + bit >>= 1; + } + j |= bit; + if i < j { + re.swap(i, j); + im.swap(i, j); + } + } + let mut len = 2; + while len <= n { + let ang = -2.0 * std::f64::consts::PI / len as f64; + for s in (0..n).step_by(len) { + for k in 0..len / 2 { + let (wr, wi) = ((ang * k as f64).cos(), (ang * k as f64).sin()); + let (a, b) = (s + k, s + k + len / 2); + let tr = re[b] * wr - im[b] * wi; + let ti = re[b] * wi + im[b] * wr; + re[b] = re[a] - tr; + im[b] = im[a] - ti; + re[a] += tr; + im[a] += ti; + } + } + len <<= 1; + } +} + +#[test] +fn sine_440_mono() { + let x = sine(440.0, 0.5, 48_000); + for rate in [64_000, 128_000] { + let (packets, d) = roundtrip(std::slice::from_ref(&x), rate); + assert!(packets.iter().all(|p| p.len() == rate as usize / 400)); + let s = snr(&x, &d[0]); + eprintln!("sine 440 Hz mono @ {} kb/s: SNR {s:.1} dB", rate / 1000); + assert!(s > 20.0, "sine SNR {s:.1} dB at {rate}"); + } +} + +#[test] +fn chirp_mono() { + let x = chirp(100.0, 10_000.0, 0.4, 96_000); + for rate in [64_000, 128_000] { + let (_, d) = roundtrip(std::slice::from_ref(&x), rate); + let s = snr(&x, &d[0]); + eprintln!("chirp 100 Hz-10 kHz mono @ {} kb/s: SNR {s:.1} dB", rate / 1000); + assert!(s > 15.0, "chirp SNR {s:.1} dB at {rate}"); + } +} + +#[test] +fn white_noise_mono() { + let x = noise(0.3, 48_000, 12345); + for rate in [64_000, 128_000] { + let (_, d) = roundtrip(std::slice::from_ref(&x), rate); + let s = snr(&x, &d[0]); + let (eo, ed) = (band_energies(&x), band_energies(&d[0])); + let err: Vec = x.iter().zip(&d[0]).map(|(a, b)| a - b).collect(); + let ee = band_energies(&err); + let octave_snr: Vec = eo.iter().zip(&ee).map(|(a, b)| a - b).collect(); + let worst = eo.iter().zip(&ed).map(|(a, b)| (a - b).abs()).fold(0.0, f64::max); + let min_octave = octave_snr.iter().copied().fold(f64::MAX, f64::min); + eprintln!( + "white noise mono @ {} kb/s: full-band SNR {s:.1} dB, octave SNRs 100 Hz-16 kHz {:.1?} dB, worst octave energy error {worst:.2} dB", + rate / 1000, + octave_snr + ); + // Full-band noise loses everything above 20 kHz (CELT's last band ends there) and spends + // few bits per coefficient on 16-20 kHz, so the full-band waveform SNR is low; the band + // energies must match closely and the octaves below 16 kHz must be waveform-coded. + assert!(worst < 1.5, "band energy error {worst:.2} dB at {rate}: {eo:?} vs {ed:?}"); + assert!(s > 0.5, "noise SNR {s:.1} dB at {rate}"); + assert!(min_octave > if rate >= 128_000 { 10.0 } else { 3.0 }, "octave SNRs {octave_snr:?} at {rate}"); + } +} + +#[test] +fn lowpass_noise_mono() { + // White noise through a 4-pole low-pass at about 2 kHz: a dense spectrum the encoder can + // afford to code as a waveform. + let w = noise(0.8, 48_000, 999); + let a = (-2.0 * std::f32::consts::PI * 2_000.0 / SR).exp(); + let mut st = [0f32; 4]; + let x: Vec = w + .iter() + .map(|&v| { + let mut y = v; + for s in st.iter_mut() { + *s = (1.0 - a) * y + a * *s; + y = *s; + } + y + }) + .collect(); + for rate in [64_000, 128_000] { + let (_, d) = roundtrip(std::slice::from_ref(&x), rate); + let s = snr(&x, &d[0]); + eprintln!("low-pass noise mono @ {} kb/s: SNR {s:.1} dB", rate / 1000); + assert!(s > 15.0, "low-pass noise SNR {s:.1} dB at {rate}"); + } +} + +#[test] +fn silence_stays_silent() { + let x = vec![0f32; 48_000]; + let (packets, d) = roundtrip(&[x.clone(), x], 128_000); + assert!(packets.iter().all(|p| p.len() == 3), "silence frames are 3-byte packets"); + for ch in &d { + let peak = ch.iter().fold(0f32, |m, v| m.max(v.abs())); + assert!(peak < 1e-6, "peak {peak}"); + } +} + +#[test] +fn stereo_distinct_channels() { + let l = sine(440.0, 0.5, 48_000); + let r = chirp(200.0, 6_000.0, 0.3, 48_000); + let (packets, d) = roundtrip(&[l.clone(), r.clone()], 256_000); + assert!(packets.iter().all(|p| p[0] == 0xFC && p.len() == 640)); + let (sl, sr) = (snr(&l, &d[0]), snr(&r, &d[1])); + eprintln!("stereo @ 256 kb/s: L (sine) SNR {sl:.1} dB, R (chirp) SNR {sr:.1} dB"); + assert!(sl > 20.0 && sr > 15.0, "stereo SNR L {sl:.1} R {sr:.1}"); + // No leakage of the chirp into the left channel (2-4 kHz octave: chirp only). + let (el, er) = (band_energies(&d[0]), band_energies(&r)); + assert!(el[4] < er[4] - 30.0, "chirp energy leaked into L: {} vs {}", el[4], er[4]); +} + +#[test] +fn stereo_128k_total() { + let l = sine(1000.0, 0.4, 48_000); + let r = sine(330.0, 0.4, 48_000); + let (_, d) = roundtrip(&[l.clone(), r.clone()], 128_000); + let (sl, sr) = (snr(&l, &d[0]), snr(&r, &d[1])); + eprintln!("stereo @ 128 kb/s: L SNR {sl:.1} dB, R SNR {sr:.1} dB"); + assert!(sl > 15.0 && sr > 15.0, "stereo SNR L {sl:.1} R {sr:.1}"); +} + +#[test] +fn signal_after_silence_and_loud_transient() { + // Silence, then a burst at full scale, then silence again: energy prediction must recover. + let mut x = vec![0f32; 24_000]; + x.extend(sine(2_000.0, 0.99, 24_000)); + x.extend(vec![0f32; 24_000]); + let (_, d) = roundtrip(std::slice::from_ref(&x), 96_000); + let s = snr(&x[24_000 + 2 * FRAME..48_000], &d[0][24_000 + 2 * FRAME..48_000]); + eprintln!("burst after silence @ 96 kb/s: SNR {s:.1} dB"); + assert!(s > 15.0, "burst SNR {s:.1}"); + let tail_peak = d[0][48_000 + 2 * FRAME..].iter().fold(0f32, |m, v| m.max(v.abs())); + assert!(tail_peak < 1e-3, "tail peak {tail_peak}"); +} + +// ---- Ogg Opus (RFC 3533 + RFC 7845) for the ffmpeg oracle ---- + +fn ogg_crc(data: &[u8]) -> u32 { + let mut crc = 0u32; + for &b in data { + crc ^= (b as u32) << 24; + for _ in 0..8 { + crc = if crc & 0x8000_0000 != 0 { (crc << 1) ^ 0x04C1_1DB7 } else { crc << 1 }; + } + } + crc +} + +fn ogg_page(out: &mut Vec, packet: &[u8], header_type: u8, granule: u64, seq: u32) { + let mut lacing = vec![255u8; packet.len() / 255]; + lacing.push((packet.len() % 255) as u8); + let start = out.len(); + out.extend_from_slice(b"OggS"); + out.push(0); + out.push(header_type); + out.extend_from_slice(&granule.to_le_bytes()); + out.extend_from_slice(&0x4543_4F50u32.to_le_bytes()); // serial + out.extend_from_slice(&seq.to_le_bytes()); + out.extend_from_slice(&0u32.to_le_bytes()); // CRC placeholder + out.push(lacing.len() as u8); + out.extend_from_slice(&lacing); + out.extend_from_slice(packet); + let crc = ogg_crc(&out[start..]); + out[start + 22..start + 26].copy_from_slice(&crc.to_le_bytes()); +} + +fn ogg_opus(enc: &OpusEncoder, packets: &[Vec], len: usize) -> Vec { + let mut out = Vec::new(); + ogg_page(&mut out, &enc.opus_head(), 0x02, 0, 0); + let vendor = b"effectcraft-opusenc"; + let mut tags = b"OpusTags".to_vec(); + tags.extend_from_slice(&(vendor.len() as u32).to_le_bytes()); + tags.extend_from_slice(vendor); + tags.extend_from_slice(&0u32.to_le_bytes()); + ogg_page(&mut out, &tags, 0, 0, 1); + let pre_skip = enc.pre_skip() as u64; + for (i, p) in packets.iter().enumerate() { + let last = i + 1 == packets.len(); + let granule = if last { pre_skip + len as u64 } else { ((i + 1) * FRAME) as u64 }; + ogg_page(&mut out, p, if last { 0x04 } else { 0 }, granule, i as u32 + 2); + } + out +} + +#[test] +fn ogg_crc_known_value() { + // CRC-32 (poly 0x04C11DB7, no reflection, init 0) of "123456789". + assert_eq!(ogg_crc(b"123456789"), 0x89A1_897F); +} + +#[test] +fn ffmpeg_oracle_decodes_ogg_opus() { + if Command::new("ffmpeg").arg("-version").output().is_err() { + eprintln!("ffmpeg not installed; skipping the oracle test"); + return; + } + for (channels, rate) in [(1u8, 96_000u32), (2, 192_000)] { + let len = 48_000; + let planar: Vec> = + if channels == 1 { vec![chirp(150.0, 9_000.0, 0.4, len)] } else { vec![sine(440.0, 0.5, len), chirp(150.0, 9_000.0, 0.3, len)] }; + let c = channels as usize; + let mut pcm = vec![0f32; len * c]; + for (ch, p) in planar.iter().enumerate() { + for (i, &v) in p.iter().enumerate() { + pcm[i * c + ch] = v; + } + } + let (enc, packets) = encode(&pcm, channels, rate); + let path = std::env::temp_dir().join(format!("effectcraft-opusenc-{}-{channels}.opus", std::process::id())); + std::fs::write(&path, ogg_opus(&enc, &packets, len)).unwrap(); + let out = Command::new("ffmpeg") + .args(["-v", "error", "-i"]) + .arg(&path) + .args(["-f", "f32le", "-ac", &channels.to_string(), "-ar", "48000", "-"]) + .output() + .expect("run ffmpeg"); + let _ = std::fs::remove_file(&path); + assert!(out.status.success(), "ffmpeg failed: {}", String::from_utf8_lossy(&out.stderr)); + let dec: Vec = out.stdout.chunks_exact(4).map(|b| f32::from_le_bytes([b[0], b[1], b[2], b[3]])).collect(); + assert_eq!(dec.len(), len * c, "ffmpeg must honour pre-skip and the end granule"); + for ch in 0..c { + let d: Vec = (0..len).map(|i| dec[i * c + ch]).collect(); + let s = snr(&planar[ch], &d); + eprintln!("ffmpeg {channels} ch @ {} kb/s, channel {ch}: SNR {s:.1} dB", rate / 1000); + assert!(s > 15.0, "ffmpeg-decoded SNR {s:.1} dB (channel {ch})"); + } + } +} diff --git a/crates/opusenc/Cargo.toml b/crates/opusenc/Cargo.toml new file mode 100644 index 00000000..90b0fd02 --- /dev/null +++ b/crates/opusenc/Cargo.toml @@ -0,0 +1,12 @@ +[package] +name = "effectcraft-opusenc" +description = "Clean-room pure-Rust Opus encoder (RFC 6716 CELT-only, fullband 48 kHz, 20 ms; RFC 7845 OpusHead) for WebM export audio" +version.workspace = true +edition.workspace = true +license.workspace = true +rust-version.workspace = true + +[lints] +workspace = true + +[dependencies] diff --git a/crates/opusenc/README.md b/crates/opusenc/README.md new file mode 100644 index 00000000..0c76d49f --- /dev/null +++ b/crates/opusenc/README.md @@ -0,0 +1,33 @@ +# effectcraft-opusenc + +Clean-room, pure-Rust Opus encoder for EffectCraft's WebM export audio (layer L0, `std` only, +no `unsafe`). + +- **CELT-only** packets (TOC configuration 31: fullband 20 kHz, 20 ms, one frame per packet), + 48 kHz input, mono or stereo (stereo coded as dual L/R stereo), constant packet size from the + requested bitrate; all-zero frames become 3-byte silence packets. +- Forward MDCT with the low-overlap window, band energies with inter/intra coarse prediction + (cheaper of the two per frame), fine and final energy bits, the decoder's exact bit allocation, + spreading and recursive band splitting, greedy PVQ search and codeword indexing. +- Not implemented: SILK/hybrid modes, transient (short-block) frames, TF changes, pitch + post-filter, intensity stereo, dynamic allocation boosts, VBR rate control. +- `OpusEncoder::pre_skip()` is 120 samples (the MDCT overlap); `opus_head()` returns the RFC 7845 + identification header. + +```rust +let mut enc = effectcraft_opusenc::OpusEncoder::new(2, 128_000); +let packet = enc.encode_float(&vec![0.0; effectcraft_opusenc::OpusEncoder::FRAME_SIZE * 2]); +``` + +## Specifications + +- RFC 6716, *Definition of the Opus Audio Codec* (September 2012), with the RFC 8251 updates + (October 2017) as implemented by the decoder we round-trip against. +- RFC 7845, *Ogg Encapsulation for the Opus Audio Codec* (April 2016), §5.1 `OpusHead` and + §4.2 pre-skip. + +Implemented from the RFC text. Constant tables (band edges, `eMeans`, the coarse-energy Laplace +model, the allocation table, ICDFs) and the decoder-side allocation arithmetic are mirrored from +FilmCraft's first-party, spec-derived `filmcraft-opus` decoder (MIT OR Apache-2.0), which is also +the round-trip oracle in `crates/export/tests/opus_roundtrip.rs` (ffmpeg is a second, external +oracle when installed). No third-party encoder or decoder source was consulted. diff --git a/crates/opusenc/src/bands.rs b/crates/opusenc/src/bands.rs new file mode 100644 index 00000000..744de3fe --- /dev/null +++ b/crates/opusenc/src/bands.rs @@ -0,0 +1,497 @@ +//! Band shape coding: PVQ search and codeword indexing, spreading, recursive band splitting +//! (RFC 6716 §4.3.4, encoder direction). +//! +//! The encoder codes 20 ms long-block frames only (no transients, no TF changes), and stereo +//! as dual (independent L/R) stereo, so a band is always coded by `quant_partition` with a single +//! block. Every bit-budget computation mirrors the decoder so both sides take the same +//! split/pulse decisions. + +use crate::range::{BITRES, RangeEncoder, ilog}; +use crate::rate::{Mode, QTHETA_OFFSET, get_pulses}; +use crate::tables::*; + +pub const SPREAD_NORMAL: usize = 2; + +#[inline] +fn frac_mul16(a: i32, b: i32) -> i32 { + (16384 + (a as i16 as i32) * (b as i16 as i32)) >> 15 +} + +fn bitexact_cos(x: i32) -> i32 { + let tmp = (4096 + x * x) >> 13; + let x2 = tmp; + let x2 = (32767 - x2) + frac_mul16(x2, -7651 + frac_mul16(x2, 8277 + frac_mul16(-626, x2))); + 1 + x2 +} + +fn bitexact_log2tan(isin: i32, icos: i32) -> i32 { + let lc = ilog(icos as u32); + let ls = ilog(isin as u32); + let icos = icos << (15 - lc); + let isin = isin << (15 - ls); + (ls - lc) * (1 << 11) + frac_mul16(isin, frac_mul16(isin, -2597) + 7932) - frac_mul16(icos, frac_mul16(icos, -2597) + 7932) +} + +fn compute_qn(n: i32, b: i32, offset: i32, pulse_cap: i32, stereo: bool) -> i32 { + const EXP2_TABLE8: [i32; 8] = [16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048]; + let mut n2 = 2 * n - 1; + if stereo && n == 2 { + n2 -= 1; + } + let mut qb = (b + n2 * offset) / n2; + qb = qb.min(b - pulse_cap - (4 << BITRES)); + qb = qb.min(8 << BITRES); + if qb < (1 << BITRES >> 1) { + 1 + } else { + let qn = EXP2_TABLE8[(qb & 0x7) as usize] >> (14 - (qb >> BITRES)); + (qn + 1) >> 1 << 1 + } +} + +/// One elementary pass of the spreading rotation, inverted (the decoder applies the forward +/// pass; the encoder undoes it before the pulse search). +fn exp_rotation1_inv(x: &mut [f32], len: usize, stride: usize, c: f32, s: f32) { + if len <= stride { + return; + } + if len >= 2 * stride + 1 { + for i in 0..=len - 2 * stride - 1 { + let x1 = x[i]; + let x2 = x[i + stride]; + x[i + stride] = c * x2 - s * x1; + x[i] = c * x1 + s * x2; + } + } + for i in (0..len - stride).rev() { + let x1 = x[i]; + let x2 = x[i + stride]; + x[i + stride] = c * x2 - s * x1; + x[i] = c * x1 + s * x2; + } +} + +/// Inverse of the decoder's spreading rotation (§4.3.4.3) for one block of `len` samples. +pub fn exp_rotation_inv(x: &mut [f32], len: usize, k: usize, spread: usize) { + const SPREAD_FACTOR: [usize; 3] = [15, 10, 5]; + if 2 * k >= len || spread == 0 { + return; + } + let factor = SPREAD_FACTOR[spread - 1]; + let gain = len as f32 / (len + factor * k) as f32; + let theta = 0.5 * gain * gain; + let c = (0.5 * std::f32::consts::PI * theta).cos(); + let s = (0.5 * std::f32::consts::PI * (1.0 - theta)).cos(); + let mut stride2 = 0; + if len >= 8 { + stride2 = 1; + while stride2 * stride2 + stride2 < len { + stride2 += 1; + } + } + // The decoder applies (stride2: s, c) then (1: c, s); undo in reverse order. + exp_rotation1_inv(x, len, 1, c, s); + if stride2 != 0 { + exp_rotation1_inv(x, len, stride2, s, c); + } +} + +/// Greedy PVQ search: the integer vector `y` with `sum |y| = k` closest in direction to `x`. +pub fn pvq_search(x: &[f32], k: usize, y: &mut [i32]) { + let n = x.len(); + y[..n].fill(0); + let abs: Vec = x.iter().map(|v| v.abs()).collect(); + let sum: f32 = abs.iter().sum(); + if !sum.is_finite() || sum <= 1e-15 { + y[0] = k as i32; + return; + } + let mut placed = 0usize; + let mut xy = 0f32; + let mut yy = 0f32; + if k > n / 2 { + // Pre-project onto the pyramid (rounding down) and place the rest greedily. + let r = (k as f32 - 1.0) / sum; + for j in 0..n { + let p = (abs[j] * r).floor() as i32; + y[j] = p; + placed += p as usize; + xy += abs[j] * p as f32; + yy += (p * p) as f32; + } + } + while placed < k { + let mut best = 0usize; + let mut best_num = -1f32; + let mut best_den = 1f32; + for j in 0..n { + let num = xy + abs[j]; + let num = num * num; + let den = yy + 2.0 * y[j] as f32 + 1.0; + if num * best_den > best_num * den { + best_num = num; + best_den = den; + best = j; + } + } + xy += abs[best]; + yy += 2.0 * y[best] as f32 + 1.0; + y[best] += 1; + placed += 1; + } + for j in 0..n { + if x[j] < 0.0 { + y[j] = -y[j]; + } + } +} + +/// Encodes the PVQ codeword `y` (`sum |y| = k`) as the index the RFC 6716 §4.3.4.2 decoder maps +/// back to `y`. +pub fn encode_pulses(y: &[i32], k: usize, enc: &mut RangeEncoder, v: &mut Vec) { + let n = y.len(); + 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); + } + } + let total = v[n * w + k]; + let vv = |nn: usize, kx: usize| v[nn * w + kx]; + let mut index = 0u64; + let mut kk = k; + for j in 0..n { + if kk == 0 { + break; + } + let rem = n - j; + let t = y[j].unsigned_abs() as usize; + // Codewords with |y_j| = t start at sum_{m < kk - t} V(rem-1, m) within their sign half. + let mut base: u64 = (0..kk - t).map(|m| vv(rem - 1, m)).sum(); + if y[j] < 0 { + base += (vv(rem - 1, kk) + vv(rem, kk)) / 2; + } + index += base; + kk -= t; + } + debug_assert!(total <= u32::MAX as u64); + if total >= 2 { + enc.uint(index as u32, total as u32); + } +} + +pub struct BandEnc<'m, 'e> { + pub m: &'m Mode, + pub enc: &'e mut RangeEncoder, + pub i: usize, + pub spread: usize, + pub remaining_bits: i32, + pvq_v: Vec, + iy: Vec, + tmp: Vec, +} + +impl BandEnc<'_, '_> { + /// `compute_theta` for a mono split (single block): measures the energy split between the + /// two halves and codes the quantised angle with the triangular pdf. + fn compute_theta_split(&mut self, xa: &[f32], xb: &[f32], b: &mut i32, lm: i32) -> (i32, i32, i32) { + let n = xa.len() as i32; + let i = self.i; + let pulse_cap = self.m.log_n[i] + lm * (1 << BITRES); + let offset = (pulse_cap >> 1) - QTHETA_OFFSET; + let qn = compute_qn(n, *b, offset, pulse_cap, false); + let tell = self.enc.tell_frac(); + let mut itheta = 0i32; + if qn != 1 { + let mid = xa.iter().map(|v| v * v).sum::().sqrt(); + let side = xb.iter().map(|v| v * v).sum::().sqrt(); + let ang = (side as f64).atan2(mid as f64) * std::f64::consts::FRAC_2_PI; + let it = ((ang * qn as f64 + 0.5).floor() as i32).clamp(0, qn) as u32; + let q = qn as u32; + let half = q >> 1; + let ft = (half + 1) * (half + 1); + let (fl, fs) = if it <= half { (it * (it + 1) >> 1, it + 1) } else { (ft - ((q + 1 - it) * (q + 2 - it) >> 1), q + 1 - it) }; + self.enc.encode(fl, fl + fs, ft); + itheta = ((it * 16384) / q) as i32; + } + let qalloc = self.enc.tell_frac() - tell; + *b -= qalloc; + let delta = if itheta == 0 { + -16384 + } else if itheta == 16384 { + 16384 + } else { + let imid = bitexact_cos(itheta); + let iside = bitexact_cos(16384 - itheta); + frac_mul16((n - 1) << 7, bitexact_log2tan(iside, imid)) + }; + (itheta, delta, qalloc) + } + + /// `quant_partition` (single block, encoder side). + pub fn quant_partition(&mut self, x: &[f32], mut b: i32, mut lm: i32) { + let n = x.len(); + let i = self.i; + let cache = self.m.cache(lm, i); + let cache_max = cache[cache[0] as usize] as i32; + if lm != -1 && b > cache_max + 12 && n > 2 { + let half = n >> 1; + let (xa, xb) = x.split_at(half); + lm -= 1; + let (itheta, delta, qalloc) = self.compute_theta_split(xa, xb, &mut b, lm); + let mut mbits = 0.max(b.min((b - delta) / 2)); + let mut sbits = b - mbits; + self.remaining_bits -= qalloc; + let mut rebalance = self.remaining_bits; + if mbits >= sbits { + self.quant_partition(xa, mbits, lm); + rebalance = mbits - (rebalance - self.remaining_bits); + if rebalance > 3 << BITRES && itheta != 0 { + sbits += rebalance - (3 << BITRES); + } + self.quant_partition(xb, sbits, lm); + } else { + self.quant_partition(xb, sbits, lm); + rebalance = sbits - (rebalance - self.remaining_bits); + if rebalance > 3 << BITRES && itheta != 16384 { + mbits += rebalance - (3 << BITRES); + } + self.quant_partition(xa, mbits, lm); + } + return; + } + let m = self.m; + let mut q = m.bits2pulses(i, lm, b); + let mut curr_bits = m.pulses2bits(i, lm, q); + self.remaining_bits -= curr_bits; + while self.remaining_bits < 0 && q > 0 { + self.remaining_bits += curr_bits; + q -= 1; + curr_bits = m.pulses2bits(i, lm, q); + self.remaining_bits -= curr_bits; + } + if q != 0 { + let k = get_pulses(q) as usize; + let mut tmp = std::mem::take(&mut self.tmp); + tmp.clear(); + tmp.extend_from_slice(x); + exp_rotation_inv(&mut tmp, n, k, self.spread); + let mut iy = std::mem::take(&mut self.iy); + iy.clear(); + iy.resize(n, 0); + pvq_search(&tmp, k, &mut iy); + encode_pulses(&iy, k, self.enc, &mut self.pvq_v); + self.iy = iy; + self.tmp = tmp; + } + } +} + +/// Codes all band shapes (`quant_all_bands`, encoder side, long blocks). `x` holds `c` channels +/// of `n` normalised MDCT coefficients each. +#[allow(clippy::too_many_arguments)] +pub fn quant_all_bands( + m: &Mode, + start: usize, + end: usize, + x: &[f32], + n_total: usize, + c: usize, + pulses: &[i32; NB_EBANDS], + spread: usize, + dual_stereo: bool, + intensity: usize, + total_bits: i32, + mut balance: i32, + enc: &mut RangeEncoder, + lm: usize, + coded_bands: usize, +) { + let mm = 1usize << lm; + debug_assert!(c == 1 || dual_stereo || intensity <= start, "stereo frames are coded as dual stereo"); + let mut ctx = BandEnc { m, enc, i: 0, spread, remaining_bits: 0, pvq_v: Vec::new(), iy: Vec::new(), tmp: Vec::new() }; + for i in start..end { + ctx.i = i; + let xo = mm * EBANDS[i] as usize; + let n = mm * EBANDS[i + 1] as usize - xo; + let tell = ctx.enc.tell_frac(); + if i != start { + balance -= tell; + } + let remaining_bits = total_bits - tell - 1; + ctx.remaining_bits = remaining_bits; + let b = if i < coded_bands { + let curr_balance = balance / 3.min(coded_bands as i32 - i as i32); + 0.max(16383.min((remaining_bits + 1).min(pulses[i] + curr_balance))) + } else { + 0 + }; + // Bands at and above `intensity` (= coded_bands) carry no bits, so they emit no symbols. + if i < coded_bands || c == 1 { + if c == 2 && dual_stereo && i < intensity { + ctx.quant_partition(&x[xo..xo + n], b / 2, lm as i32); + ctx.quant_partition(&x[n_total + xo..n_total + xo + n], b / 2, lm as i32); + } else { + ctx.quant_partition(&x[xo..xo + n], b, lm as i32); + } + } + balance += pulses[i] + tell; + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::range::dec::RangeDecoder; + + fn exp_rotation1(x: &mut [f32], len: usize, stride: usize, c: f32, s: f32) { + let ms = -s; + if len <= stride { + return; + } + for i in 0..len - stride { + let x1 = x[i]; + let x2 = x[i + stride]; + x[i + stride] = c * x2 + s * x1; + x[i] = c * x1 + ms * x2; + } + if len >= 2 * stride + 1 { + for i in (0..=len - 2 * stride - 1).rev() { + let x1 = x[i]; + let x2 = x[i + stride]; + x[i + stride] = c * x2 + s * x1; + x[i] = c * x1 + ms * x2; + } + } + } + + /// The decoder's spreading rotation (single block). + fn exp_rotation(x: &mut [f32], len: usize, k: usize, spread: usize) { + const SPREAD_FACTOR: [usize; 3] = [15, 10, 5]; + if 2 * k >= len || spread == 0 { + return; + } + let factor = SPREAD_FACTOR[spread - 1]; + let gain = len as f32 / (len + factor * k) as f32; + let theta = 0.5 * gain * gain; + let c = (0.5 * std::f32::consts::PI * theta).cos(); + let s = (0.5 * std::f32::consts::PI * (1.0 - theta)).cos(); + let mut stride2 = 0; + if len >= 8 { + stride2 = 1; + while stride2 * stride2 + stride2 < len { + stride2 += 1; + } + } + if stride2 != 0 { + exp_rotation1(x, len, stride2, s, c); + } + exp_rotation1(x, len, 1, c, s); + } + + #[test] + fn spreading_inverse_is_exact() { + for &len in &[8usize, 16, 24, 48, 96, 176] { + for k in [1usize, 2, 3, 5] { + let orig: Vec = (0..len).map(|i| ((i * 31 % 17) as f32 - 8.0) / 8.0).collect(); + let mut x = orig.clone(); + exp_rotation_inv(&mut x, len, k, SPREAD_NORMAL); + exp_rotation(&mut x, len, k, SPREAD_NORMAL); + for (a, b) in x.iter().zip(&orig) { + assert!((a - b).abs() < 1e-4, "len={len} k={k}"); + } + } + } + } + + /// PVQ decoder of RFC 6716 §4.3.4.2 (as in FilmCraft's decoder). + fn decode_pulses(n: usize, k: usize, dec: &mut RangeDecoder) -> Vec { + let w = k + 1; + let mut v = vec![0u64; (n + 1) * w]; + v[0] = 1; + for nn in 1..=n { + v[nn * w] = 1; + for kk in 1..=k { + v[nn * w + kk] = v[(nn - 1) * w + kk].saturating_add(v[nn * w + kk - 1]).saturating_add(v[(nn - 1) * w + kk - 1]); + } + } + let total = v[n * w + k]; + let mut i = if total >= 2 { dec.uint(total as u32) as u64 } else { 0 }; + let mut y = vec![0i32; n]; + let mut kk = k; + for j in 0..n { + let rem = n - j; + if kk == 0 { + break; + } + 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 { + i -= p; + } + let k0 = kk; + p -= vv(rem - 1, kk); + while p > i && kk > 0 { + kk -= 1; + p -= vv(rem - 1, kk); + } + let mag = (k0 - kk) as i32; + y[j] = if neg { -mag } else { mag }; + i -= p.min(i); + } + y + } + + fn all_codewords(n: usize, k: usize) -> Vec> { + if n == 0 { + return if k == 0 { vec![vec![]] } else { vec![] }; + } + let mut out = Vec::new(); + for t in -(k as i32)..=(k as i32) { + for mut rest in all_codewords(n - 1, k - t.unsigned_abs() as usize) { + rest.insert(0, t); + out.push(rest); + } + } + out + } + + #[test] + fn pvq_index_roundtrip_exhaustive() { + let mut v = Vec::new(); + for n in 1..6usize { + for k in 1..6usize { + for y in all_codewords(n, k) { + let mut e = RangeEncoder::new(); + encode_pulses(&y, k, &mut e, &mut v); + let buf = e.done(16); + let mut d = RangeDecoder::new(&buf); + assert_eq!(decode_pulses(n, k, &mut d), y, "n={n} k={k}"); + } + } + } + } + + #[test] + fn pvq_search_finds_good_direction() { + let x: Vec = (0..32).map(|i| ((i as f32) * 0.7).sin()).collect(); + let mut y = vec![0i32; 32]; + for k in [1usize, 4, 16, 64] { + pvq_search(&x, k, &mut y); + assert_eq!(y.iter().map(|v| v.unsigned_abs() as usize).sum::(), k); + let xy: f32 = x.iter().zip(&y).map(|(a, &b)| a * b as f32).sum(); + let yy: f32 = y.iter().map(|&b| (b * b) as f32).sum(); + let xx: f32 = x.iter().map(|a| a * a).sum(); + let corr = xy / (xx * yy).sqrt(); + assert!(corr > if k >= 16 { 0.9 } else { 0.2 }, "k={k} corr={corr}"); + } + } +} diff --git a/crates/opusenc/src/energy.rs b/crates/opusenc/src/energy.rs new file mode 100644 index 00000000..0057ae33 --- /dev/null +++ b/crates/opusenc/src/energy.rs @@ -0,0 +1,230 @@ +//! Band energy quantisation: Laplace-coded coarse energy with inter/intra prediction, fine +//! energy and the final refinement bits (RFC 6716 §4.3.2, encoder direction). +//! +//! The quantised energies are updated with exactly the decoder's arithmetic so the prediction +//! state never drifts between encoder and decoder. + +use crate::range::RangeEncoder; +use crate::rate::MAX_FINE_BITS; +use crate::tables::*; + +const LAPLACE_MINP: u32 = 1; +const LAPLACE_NMIN: u32 = 16; + +fn laplace_get_freq1(fs0: u32, decay: u32) -> u32 { + let ft = 32768 - LAPLACE_MINP * (2 * LAPLACE_NMIN) - fs0; + (ft * (16384 - decay)) >> 15 +} + +/// Encodes a Laplace-distributed integer (`fs` = probability of 0 in Q15, `decay` in Q14). +/// Returns the value actually coded (very large magnitudes are clamped to the representable +/// range). +pub fn laplace_encode(enc: &mut RangeEncoder, value: i32, mut fs: u32, decay: u32) -> i32 { + let mut fl = 0u32; + let mut coded = value; + if value != 0 { + let s: i32 = if value < 0 { -1 } else { 0 }; + let val = value.abs(); + fl = fs; + fs = laplace_get_freq1(fs, decay); + let mut i = 1; + while fs > 0 && i < val { + fs *= 2; + fl += fs + 2 * LAPLACE_MINP; + fs = (fs * decay) >> 15; + i += 1; + } + if fs == 0 { + let ndi_max = (32768 - fl + LAPLACE_MINP - 1) as i32; + let ndi_max = (ndi_max - s) >> 1; + let di = (val - i).min(ndi_max - 1); + fl = (fl as i32 + (2 * di + 1 + s) * LAPLACE_MINP as i32) as u32; + fs = LAPLACE_MINP.min(32768 - fl); + let mag = i + di; + coded = if s < 0 { -mag } else { mag }; + } else { + fs += LAPLACE_MINP; + if s == 0 { + fl += fs; + } + } + } + enc.encode_bin(fl, fl + fs, 15); + coded +} + +/// `quant_coarse_energy` for one prediction mode. `e` holds the target log2 energies (relative +/// to `E_MEANS`), `old_e` the quantised energies, updated in place. `budget` is the frame size in +/// bits. +#[allow(clippy::too_many_arguments)] +pub fn quant_coarse_energy( + start: usize, + end: usize, + e: &[f32; 2 * NB_EBANDS], + old_e: &mut [f32; 2 * NB_EBANDS], + intra: bool, + enc: &mut RangeEncoder, + c: usize, + lm: usize, + budget: i32, +) { + 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]) }; + let mut prev = [0f32; 2]; + for i in start..end { + for ch in 0..c { + let idx = i + ch * NB_EBANDS; + let oldc = old_e[idx].max(-9.0); + let x = e[idx]; + let f = x - coef * oldc - prev[ch]; + let mut qi = (f + 0.5).floor().clamp(-64.0, 64.0) as i32; + let tell = enc.tell(); + // Keep a few bits per remaining band in reserve when the budget runs out. + let bits_left = budget - tell - 3 * c as i32 * (end - i) as i32; + if i != start && bits_left < 30 { + if bits_left < 24 { + qi = qi.min(1); + } + if bits_left < 16 { + qi = qi.max(-1); + } + } + if budget - tell >= 15 { + let pi = 2 * i.min(20); + qi = laplace_encode(enc, qi, (prob[pi] as u32) << 7, (prob[pi + 1] as u32) << 6); + } else if budget - tell >= 2 { + qi = qi.clamp(-1, 1); + enc.icdf(if qi < 0 { (-2 * qi - 1) as usize } else { (2 * qi) as usize }, &SMALL_ENERGY_ICDF, 2); + } else if budget - tell >= 1 { + qi = qi.clamp(-1, 0); + enc.bit_logp(qi == -1, 1); + } else { + qi = -1; + } + let q = qi as f32; + old_e[idx] = oldc; + let tmp = coef * old_e[idx] + prev[ch] + q; + old_e[idx] = tmp; + prev[ch] = prev[ch] + q - beta * q; + } + } +} + +/// `quant_fine_energy`: raw bits refining each band's quantised energy. +pub fn quant_fine_energy( + start: usize, + end: usize, + e: &[f32; 2 * NB_EBANDS], + old_e: &mut [f32; 2 * NB_EBANDS], + fine_quant: &[i32; NB_EBANDS], + enc: &mut RangeEncoder, + c: usize, +) { + for i in start..end { + let fq = fine_quant[i]; + if fq <= 0 { + continue; + } + for ch in 0..c { + let idx = i + ch * NB_EBANDS; + let err = e[idx] - old_e[idx]; + let levels = 1i32 << fq; + let q2 = (((err + 0.5) * levels as f32).floor() as i32).clamp(0, levels - 1); + enc.bits(q2 as u32, fq as u32); + let offset = (q2 as f32 + 0.5) * (1 << (14 - fq)) as f32 / 16384.0 - 0.5; + old_e[idx] += offset; + } + } +} + +/// `quant_energy_finalise`: spends the bits left at the end of the frame on one more bit of +/// energy resolution per band, in priority order. +#[allow(clippy::too_many_arguments)] +pub fn quant_energy_finalise( + start: usize, + end: usize, + e: &[f32; 2 * NB_EBANDS], + old_e: &mut [f32; 2 * NB_EBANDS], + fine_quant: &[i32; NB_EBANDS], + fine_priority: &[i32; NB_EBANDS], + mut bits_left: i32, + enc: &mut RangeEncoder, + c: usize, +) { + for prio in 0..2 { + let mut i = start; + while i < end && bits_left >= c as i32 { + if fine_quant[i] >= MAX_FINE_BITS || fine_priority[i] != prio { + i += 1; + continue; + } + for ch in 0..c { + let idx = i + ch * NB_EBANDS; + let q2 = (e[idx] - old_e[idx] >= 0.0) as u32; + enc.bits(q2, 1); + let offset = (q2 as f32 - 0.5) * (1 << (14 - fine_quant[i] - 1)) as f32 / 16384.0; + old_e[idx] += offset; + bits_left -= 1; + } + i += 1; + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::range::dec::RangeDecoder; + + /// Laplace decoder of RFC 6716 §4.3.2.1 (as in FilmCraft's decoder). + fn laplace_decode(dec: &mut RangeDecoder, mut fs: u32, decay: u32) -> i32 { + let mut val = 0i32; + let fm = dec.decode(32768); + let mut fl = 0u32; + if fm >= fs { + val += 1; + fl = fs; + fs = laplace_get_freq1(fs, decay) + LAPLACE_MINP; + while fs > LAPLACE_MINP && fm >= fl + 2 * fs { + fs *= 2; + fl += fs; + fs = (((fs - 2 * LAPLACE_MINP) * decay) >> 15) + LAPLACE_MINP; + val += 1; + } + if fs <= LAPLACE_MINP { + let di = (fm - fl) >> 1; + val += di as i32; + fl += 2 * di * LAPLACE_MINP; + } + if fm < fl + fs { + val = -val; + } else { + fl += fs; + } + } + dec.update(fl, (fl + fs).min(32768), 32768); + val + } + + #[test] + fn laplace_roundtrip_all_models() { + for lm in 0..4 { + for intra in 0..2 { + for band in 0..21 { + let prob = &E_PROB_MODEL[lm][intra]; + let (fs, decay) = ((prob[2 * band] as u32) << 7, (prob[2 * band + 1] as u32) << 6); + let vals: Vec = (-64..=64).collect(); + let mut e = RangeEncoder::new(); + let coded: Vec = vals.iter().map(|&v| laplace_encode(&mut e, v, fs, decay)).collect(); + let size = (e.tell() as usize).div_ceil(8); + let buf = e.done(size); + let mut d = RangeDecoder::new(&buf); + for (&v, &c) in vals.iter().zip(&coded) { + assert!(c == v || (c.signum() == v.signum() && c.abs() < v.abs()), "clamp {v} -> {c}"); + assert_eq!(laplace_decode(&mut d, fs, decay), c, "lm={lm} intra={intra} band={band} v={v}"); + } + } + } + } + } +} diff --git a/crates/opusenc/src/lib.rs b/crates/opusenc/src/lib.rs new file mode 100644 index 00000000..1e7a2cec --- /dev/null +++ b/crates/opusenc/src/lib.rs @@ -0,0 +1,311 @@ +//! Clean-room pure-Rust Opus encoder (RFC 6716 CELT-only mode; RFC 7845 `OpusHead`). +//! +//! Produces fullband (20 kHz) 48 kHz Opus packets of 20 ms, mono or stereo, at a constant +//! packet size derived from the requested bitrate. Each packet is a single CELT frame (TOC +//! configuration 31, frame-count code 0). Stereo is coded as dual (independent L/R) stereo. +//! +//! Layer L0: no dependencies beyond `std`; no `unsafe`; builds for `wasm32-unknown-unknown`. +//! +//! ``` +//! use effectcraft_opusenc::OpusEncoder; +//! let mut enc = OpusEncoder::new(2, 128_000); +//! let pcm = vec![0.0f32; OpusEncoder::FRAME_SIZE * 2]; +//! let packet = enc.encode_float(&pcm); +//! assert_eq!(packet[0] >> 3, 31); // CELT-only, fullband, 20 ms +//! assert_eq!(enc.opus_head().len(), 19); +//! ``` + +// The bit-allocation code mirrors the normative integer arithmetic; keep C-like shift expressions. +#![allow(clippy::precedence, clippy::int_plus_one)] + +mod bands; +mod energy; +mod mdct; +mod range; +mod rate; +mod tables; + +use std::sync::OnceLock; + +use bands::{SPREAD_NORMAL, quant_all_bands}; +use energy::{quant_coarse_energy, quant_energy_finalise, quant_fine_energy}; +use mdct::Mdct; +use range::{BITRES, RangeEncoder}; +use rate::{Mode, compute_allocation}; +use tables::*; + +/// Samples per channel in one 20 ms frame at 48 kHz. +const N: usize = 960; +/// LM of a 20 ms frame (8 short blocks). +const LM: usize = 3; +/// Smallest CELT payload the encoder produces (bytes). +const MIN_FRAME_BYTES: usize = 16; +/// Largest CELT payload allowed in one frame (RFC 6716 §3.4, R2). +const MAX_FRAME_BYTES: usize = 1275; + +fn mode() -> &'static Mode { + static MODE: OnceLock = OnceLock::new(); + MODE.get_or_init(Mode::new) +} + +/// Opus encoder for one mono or stereo stream at 48 kHz. +pub struct OpusEncoder { + channels: usize, + /// CELT payload bytes per packet (the packet is one byte longer: the TOC). + frame_bytes: usize, + mdct: Mdct, + /// Last input sample of each channel (pre-emphasis filter memory, ×32768). + preemph_mem: [f32; 2], + /// Last `OVERLAP` pre-emphasised samples of each channel. + hist: Vec, + /// Quantised band energies, exactly as the decoder holds them. + old_band_e: [f32; 2 * NB_EBANDS], + block: Vec, + freq: Vec, +} + +impl OpusEncoder { + /// Samples per channel per packet (20 ms at 48 kHz). + pub const FRAME_SIZE: usize = N; + + /// Creates an encoder for `channels` (1 or 2) channels of 48 kHz audio at `bitrate_bps` + /// (total for all channels; packets are constant-size, clamped to 6.8–510 kb/s). + pub fn new(channels: u8, bitrate_bps: u32) -> Self { + assert!(channels == 1 || channels == 2, "Opus family-0 streams have 1 or 2 channels"); + let c = channels as usize; + let packet_bytes = (bitrate_bps as usize).div_ceil(400); + let frame_bytes = packet_bytes.saturating_sub(1).clamp(MIN_FRAME_BYTES, MAX_FRAME_BYTES); + OpusEncoder { + channels: c, + frame_bytes, + mdct: Mdct::new(N), + preemph_mem: [0.0; 2], + hist: vec![0.0; c * OVERLAP], + old_band_e: [-28.0; 2 * NB_EBANDS], + block: vec![0.0; N + OVERLAP], + freq: vec![0.0; c * N], + } + } + + /// Number of channels. + pub fn channels(&self) -> u8 { + self.channels as u8 + } + + /// Actual bitrate in bits per second (packet size × 50 packets per second). + pub fn bitrate(&self) -> u32 { + ((self.frame_bytes + 1) * 8 * 50) as u32 + } + + /// Samples at 48 kHz the decoder must discard from the start of the stream (the MDCT + /// overlap); this is the `OpusHead` pre-skip. + pub fn pre_skip(&self) -> u16 { + OVERLAP as u16 + } + + /// The 19-byte `OpusHead` identification header (RFC 7845 §5.1): version 1, channel count, + /// pre-skip, input sample rate 48000, output gain 0, channel mapping family 0. + pub fn opus_head(&self) -> Vec { + let mut h = Vec::with_capacity(19); + h.extend_from_slice(b"OpusHead"); + h.push(1); + h.push(self.channels as u8); + h.extend_from_slice(&self.pre_skip().to_le_bytes()); + h.extend_from_slice(&48_000u32.to_le_bytes()); + h.extend_from_slice(&0i16.to_le_bytes()); + h.push(0); + h + } + + /// Encodes one 20 ms frame of interleaved f32 samples in [-1, 1] (exactly + /// `FRAME_SIZE * channels` values) and returns one Opus packet (TOC byte + CELT frame). + pub fn encode_float(&mut self, pcm: &[f32]) -> Vec { + let c = self.channels; + assert_eq!(pcm.len(), N * c, "encode_float takes exactly FRAME_SIZE samples per channel"); + + // Pre-emphasis, the MDCT and the band energies. + let mut band_e = [0f32; 2 * NB_EBANDS]; + let mut log_e = [0f32; 2 * NB_EBANDS]; + let mut peak = 0f32; + for ch in 0..c { + self.block[..OVERLAP].copy_from_slice(&self.hist[ch * OVERLAP..(ch + 1) * OVERLAP]); + let mut mem = self.preemph_mem[ch]; + for j in 0..N { + let s = pcm[j * c + ch]; + let s = if s.is_finite() { s.clamp(-2.0, 2.0) * 32768.0 } else { 0.0 }; + let v = s - PREEMPH * mem; + mem = s; + self.block[OVERLAP + j] = v; + } + self.preemph_mem[ch] = mem; + self.hist[ch * OVERLAP..(ch + 1) * OVERLAP].copy_from_slice(&self.block[N..]); + peak = self.block.iter().fold(peak, |p, v| p.max(v.abs())); + let freq = &mut self.freq[ch * N..(ch + 1) * N]; + self.mdct.forward(&self.block, freq); + for i in 0..NB_EBANDS { + let lo = (EBANDS[i] as usize) << LM; + let hi = (EBANDS[i + 1] as usize) << LM; + let e = (1e-27f32 + freq[lo..hi].iter().map(|v| v * v).sum::()).sqrt(); + band_e[ch * NB_EBANDS + i] = e; + log_e[ch * NB_EBANDS + i] = e.log2() - E_MEANS[i]; + let g = 1.0 / e; + for v in freq[lo..hi].iter_mut() { + *v *= g; + } + } + } + let silence = peak < 1e-4; + let toc = (31u8 << 3) | if c == 2 { 4 } else { 0 }; + let mut packet = vec![toc]; + packet.extend_from_slice(&self.encode_celt(&log_e, silence)); + packet + } + + /// Codes one CELT frame (RFC 6716 §4.3, in bitstream order). + fn encode_celt(&mut self, log_e: &[f32; 2 * NB_EBANDS], silence: bool) -> Vec { + let m = mode(); + let c = self.channels; + let (start, end) = (0usize, NB_EBANDS); + let mut enc = RangeEncoder::new(); + if c == 1 { + for i in 0..NB_EBANDS { + self.old_band_e[i] = self.old_band_e[i].max(self.old_band_e[NB_EBANDS + i]); + } + } + if silence { + // The silence flag; the decoder treats every remaining bit as consumed. + enc.bit_logp(true, 15); + self.old_band_e = [-28.0; 2 * NB_EBANDS]; + return enc.done(2); + } + let len = self.frame_bytes; + let total_bits_raw = len as i32 * 8; + enc.bit_logp(false, 15); + // No pitch post-filter. + if enc.tell() + 16 <= total_bits_raw { + enc.bit_logp(false, 1); + } + // Long blocks only. + if enc.tell() + 3 <= total_bits_raw { + enc.bit_logp(false, 3); + } + // Coarse energy, choosing intra or inter prediction by trial (whichever is cheaper). + if enc.tell() + 3 <= total_bits_raw { + let mut enc_intra = enc.clone(); + let mut old_intra = self.old_band_e; + enc_intra.bit_logp(true, 3); + quant_coarse_energy(start, end, log_e, &mut old_intra, true, &mut enc_intra, c, LM, total_bits_raw); + enc.bit_logp(false, 3); + let mut old_inter = self.old_band_e; + quant_coarse_energy(start, end, log_e, &mut old_inter, false, &mut enc, c, LM, total_bits_raw); + if enc_intra.tell_frac() < enc.tell_frac() { + enc = enc_intra; + self.old_band_e = old_intra; + } else { + self.old_band_e = old_inter; + } + } else { + quant_coarse_energy(start, end, log_e, &mut self.old_band_e, false, &mut enc, c, LM, total_bits_raw); + } + // TF resolution: no changes. + { + let mut budget = total_bits_raw; + let mut tell = enc.tell(); + let mut logp = 4; + let tf_select_rsv = LM > 0 && tell + logp + 1 <= budget; + budget -= tf_select_rsv as i32; + for _ in start..end { + if tell + logp <= budget { + enc.bit_logp(false, logp as u32); + tell = enc.tell(); + } + logp = 5; + } + if tf_select_rsv && TF_SELECT_TABLE[LM][0] != TF_SELECT_TABLE[LM][2] { + enc.bit_logp(false, 1); + } + } + if enc.tell() + 4 <= total_bits_raw { + enc.icdf(SPREAD_NORMAL, &SPREAD_ICDF, 5); + } + // Dynamic allocation: no boosts. + let cap = m.init_caps(LM, c); + let offsets = [0i32; NB_EBANDS]; + let total_bits = total_bits_raw << BITRES; + let mut tellf = enc.tell_frac(); + for &cp in cap.iter().take(end).skip(start) { + if tellf + (6 << BITRES) < total_bits && 0 < cp { + enc.bit_logp(false, 6); + tellf = enc.tell_frac(); + } + } + let alloc_trim = 5; + if tellf + (6 << BITRES) <= total_bits { + enc.icdf(alloc_trim as usize, &TRIM_ICDF, 7); + } + let bits = (total_bits_raw << BITRES) - enc.tell_frac() - 1; + let alloc = compute_allocation(m, start, end, &offsets, &cap, alloc_trim, bits, c, LM, &mut enc); + quant_fine_energy(start, end, log_e, &mut self.old_band_e, &alloc.fine_quant, &mut enc, c); + quant_all_bands( + m, + start, + end, + &self.freq, + N, + c, + &alloc.pulses, + SPREAD_NORMAL, + alloc.dual_stereo, + alloc.intensity, + total_bits_raw << BITRES, + alloc.balance, + &mut enc, + LM, + alloc.coded_bands, + ); + let bits_left = total_bits_raw - enc.tell(); + quant_energy_finalise(start, end, log_e, &mut self.old_band_e, &alloc.fine_quant, &alloc.fine_priority, bits_left, &mut enc, c); + debug_assert!(enc.tell() <= total_bits_raw, "frame overflow: {} > {}", enc.tell(), total_bits_raw); + if c == 1 { + let (a, b) = self.old_band_e.split_at_mut(NB_EBANDS); + b.copy_from_slice(a); + } + enc.done(len) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn opus_head_layout() { + let e = OpusEncoder::new(2, 96_000); + let h = e.opus_head(); + assert_eq!(&h[..8], b"OpusHead"); + assert_eq!(h[8], 1); + assert_eq!(h[9], 2); + assert_eq!(u16::from_le_bytes([h[10], h[11]]), e.pre_skip()); + assert_eq!(u32::from_le_bytes([h[12], h[13], h[14], h[15]]), 48_000); + assert_eq!(&h[16..], &[0, 0, 0]); + } + + #[test] + fn packet_sizes_follow_bitrate() { + for (ch, rate) in [(1u8, 64_000u32), (2, 128_000), (2, 256_000)] { + let mut e = OpusEncoder::new(ch, rate); + let pcm: Vec = (0..N * ch as usize).map(|i| (i as f32 * 0.05).sin() * 0.5).collect(); + let p = e.encode_float(&pcm); + assert_eq!(p.len(), rate as usize / 400); + assert_eq!(p[0], 0xF8 | if ch == 2 { 4 } else { 0 }); + assert_eq!(e.bitrate(), rate); + } + } + + #[test] + fn silence_is_a_tiny_packet() { + let mut e = OpusEncoder::new(1, 64_000); + let p = e.encode_float(&vec![0.0; N]); + assert_eq!(p.len(), 3); + } +} diff --git a/crates/opusenc/src/mdct.rs b/crates/opusenc/src/mdct.rs new file mode 100644 index 00000000..814de037 --- /dev/null +++ b/crates/opusenc/src/mdct.rs @@ -0,0 +1,293 @@ +//! Forward MDCT with the CELT low-overlap window (RFC 6716 §4.3.7, run in the analysis +//! direction). +//! +//! A 20 ms frame is a 1920-point MDCT (960 coefficients) whose window is zero for 420 samples, +//! rises over the 120-sample overlap, is flat for 840 samples, falls over the next 120 samples +//! and is zero again for the last 420. The input block is therefore the 120 overlap samples kept +//! from the previous frame followed by the 960 new samples. The MDCT is folded into a DCT-IV, +//! which is computed with an N/2-point complex FFT. + +use crate::tables::OVERLAP; + +#[derive(Clone, Copy, Default, Debug)] +pub struct Cpx { + pub re: f32, + pub im: f32, +} + +impl Cpx { + #[inline] + fn mul(self, o: Cpx) -> Cpx { + Cpx { re: self.re * o.re - self.im * o.im, im: self.re * o.im + self.im * o.re } + } + #[inline] + fn add(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 (decimation in time, recursive). +pub struct Fft { + n: usize, + factors: Vec, + twiddle: Vec, +} + +impl Fft { + pub fn new(n: usize) -> Fft { + let mut factors = Vec::new(); + let mut m = n; + for p in [4usize, 2, 3, 5] { + while m.is_multiple_of(p) { + factors.push(p); + m /= p; + } + } + assert_eq!(m, 1, "unsupported FFT size {n}"); + let twiddle = (0..n) + .map(|k| { + let ph = -2.0 * std::f64::consts::PI * k as f64 / n as f64; + Cpx { re: ph.cos() as f32, im: ph.sin() as f32 } + }) + .collect(); + Fft { n, factors, twiddle } + } + + pub fn process(&self, data: &mut [Cpx], scratch: &mut Vec) { + scratch.clear(); + scratch.extend_from_slice(&data[..self.n]); + self.rec(&mut data[..self.n], scratch, 1, self.n, 0, 1); + } + + fn rec(&self, out: &mut [Cpx], inp: &[Cpx], in_stride: usize, n: usize, fi: usize, tw_stride: usize) { + let p = self.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 { + 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, tq) in t.iter_mut().enumerate().take(p) { + *tq = out[q * m + k].mul(self.twiddle[(q * k * tw_stride) % big]); + } + for q2 in 0..p { + let mut acc = Cpx::default(); + for (q, tq) in t.iter().enumerate().take(p) { + acc = acc.add(tq.mul(self.twiddle[((q * q2) % p) * (big / p)])); + } + out[q2 * m + k] = acc; + } + } + } +} + +/// Forward MDCT producing `n` coefficients from a block of `n + OVERLAP` samples. +pub struct Mdct { + /// Number of coefficients (half the MDCT length). + n: usize, + fft: Fft, + pre: Vec, + post: Vec, + /// Rising half of the low-overlap window (RFC 6716 §4.3.7). + pub window: [f32; OVERLAP], + scale: f32, + z: Vec, + scratch: Vec, + v: Vec, +} + +impl Mdct { + pub fn new(n: usize) -> Mdct { + let n2 = n / 2; + let pi = std::f64::consts::PI; + let pre = (0..n2) + .map(|m| { + let ph = -pi * m as f64 / n as f64; + Cpx { re: ph.cos() as f32, im: ph.sin() as f32 } + }) + .collect(); + let post = (0..n2) + .map(|p| { + let ph = -pi * (p as f64 + 0.25) / n as f64; + Cpx { re: ph.cos() as f32, im: ph.sin() as f32 } + }) + .collect(); + let mut window = [0f32; OVERLAP]; + for (i, w) in window.iter_mut().enumerate() { + let x = std::f64::consts::FRAC_PI_2 * (i as f64 + 0.5) / OVERLAP as f64; + let s = x.sin(); + *w = (std::f64::consts::FRAC_PI_2 * s * s).sin() as f32; + } + Mdct { n, fft: Fft::new(n2), pre, post, window, scale: 2.0 / n as f32, z: Vec::new(), scratch: Vec::new(), v: vec![0.0; n] } + } + + /// Transforms `block` (`n + OVERLAP` samples: the previous frame's last `OVERLAP` samples, then + /// `n` new ones) into `out[..n]`, scaled so that the RFC 6716 decoder's inverse MDCT and + /// overlap-add reproduce the input. + pub fn forward(&mut self, block: &[f32], out: &mut [f32]) { + let n = self.n; + let n2 = n / 2; + debug_assert_eq!(block.len(), n + OVERLAP); + // Place the windowed block into the (virtual) 2n-sample MDCT input u at offset `pad`. + let pad = (n - OVERLAP) / 2; + let w = &self.window; + let u = |t: usize| -> f32 { + if t < pad || t >= pad + n + OVERLAP { + return 0.0; + } + let j = t - pad; + let g = if j < OVERLAP { + w[j] + } else if j >= n { + w[OVERLAP - 1 - (j - n)] + } else { + 1.0 + }; + block[j] * g + }; + // Fold to a DCT-IV input: v = (-c_r - d, a - b_r) over the four quarters of u. + for m in 0..n2 { + self.v[m] = -u(n + n2 - 1 - m) - u(n + n2 + m); + self.v[n2 + m] = u(m) - u(n2 + n2 - 1 - m); + } + // DCT-IV via an n/2-point complex FFT. + self.z.clear(); + for m in 0..n2 { + let wv = Cpx { re: self.v[2 * m], im: self.v[n - 1 - 2 * m] }; + self.z.push(wv.mul(self.pre[m])); + } + self.fft.process(&mut self.z, &mut self.scratch); + for p in 0..n2 { + let y = self.z[p].mul(self.post[p]); + out[2 * p] = y.re * self.scale; + out[n - 1 - 2 * p] = -y.im * self.scale; + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Direct MDCT definition X[k] = sum u[t] cos(pi/N (t + 1/2 + N/2)(k + 1/2)). + fn direct(block: &[f32], n: usize, window: &[f32; OVERLAP]) -> Vec { + let pad = (n - OVERLAP) / 2; + let mut u = vec![0f64; 2 * n]; + for (j, &b) in block.iter().enumerate() { + let g = if j < OVERLAP { + window[j] + } else if j >= n { + window[OVERLAP - 1 - (j - n)] + } else { + 1.0 + }; + u[pad + j] = (b * g) as f64; + } + (0..n) + .map(|k| { + (0..2 * n).map(|t| u[t] * (std::f64::consts::PI / n as f64 * (t as f64 + 0.5 + n as f64 / 2.0) * (k as f64 + 0.5)).cos()).sum::() * 2.0 + / n as f64 + }) + .collect() + } + + #[test] + fn matches_direct_formula() { + for &n in &[240usize, 960] { + let block: Vec = (0..n + OVERLAP).map(|i| ((i * 37 % 101) as f32 - 50.0) / 7.0).collect(); + let mut m = Mdct::new(n); + let mut out = vec![0f32; n]; + m.forward(&block, &mut out); + let d = direct(&block, n, &m.window); + for k in 0..n { + assert!((out[k] as f64 - d[k]).abs() < 1e-2, "n={n} k={k}: {} vs {}", out[k], d[k]); + } + } + } + + // ---- Replica of the RFC 6716 decoder inverse MDCT (as in FilmCraft's decoder) ---- + + fn imdct_backward(n_full: usize, input: &[f32], out: &mut [f32], window: &[f32]) { + let n2 = n_full >> 1; + let n4 = n_full >> 2; + let t: Vec = (0..n2).map(|i| (2.0 * std::f64::consts::PI * (i as f64 + 0.125) / n_full as f64).cos() as f32).collect(); + let fft = Fft::new(n4); + let mut z = Vec::new(); + for i in 0..n4 { + let x1 = input[2 * i]; + let x2 = input[n2 - 1 - 2 * i]; + let yr = x2 * t[i] + x1 * t[n4 + i]; + let yi = x1 * t[i] - x2 * t[n4 + i]; + z.push(Cpx { re: yi, im: yr }); + } + fft.process(&mut z, &mut Vec::new()); + let overlap = OVERLAP; + let base = overlap >> 1; + let y = &mut out[base..base + n2]; + for i in 0..n4.div_ceil(2) { + let front = z[i]; + let back = z[n4 - 1 - i]; + let (re, im) = (front.im, front.re); + let (t0, t1) = (t[i], t[n4 + i]); + let yr0 = re * t0 + im * t1; + let yi0 = re * t1 - im * t0; + let (re, im) = (back.im, back.re); + let (t0, t1) = (t[n4 - i - 1], t[n2 - i - 1]); + let yr1 = re * t0 + im * t1; + let yi1 = re * t1 - im * t0; + y[2 * i] = yr0; + y[n2 - 1 - 2 * i] = yi0; + y[n2 - 2 - 2 * i] = yr1; + y[2 * i + 1] = yi1; + } + for i in 0..overlap / 2 { + let x1 = out[overlap - 1 - i]; + let x2 = out[i]; + let w1 = window[i]; + let w2 = window[overlap - 1 - i]; + out[i] = x2 * w2 - x1 * w1; + out[overlap - 1 - i] = x2 * w1 + x1 * w2; + } + } + + /// MDCT analysis followed by the decoder's IMDCT + TDAC overlap-add reproduces the input + /// delayed by exactly `OVERLAP` samples. + #[test] + fn perfect_reconstruction_through_decoder_imdct() { + let n = 960; + let frames = 5; + let sig: Vec = (0..n * frames).map(|i| (i as f32 * 0.013).sin() * 1000.0 + ((i * 7919) % 211) as f32 - 105.0).collect(); + let mut m = Mdct::new(n); + let mut hist = vec![0f32; OVERLAP]; + // Decoder memory: [frame output n][folded tail OVERLAP/2] laid out like decode_mem. + let mut mem = vec![0f32; n + OVERLAP]; + let mut decoded = Vec::new(); + for f in 0..frames { + let mut block = hist.clone(); + block.extend_from_slice(&sig[f * n..(f + 1) * n]); + hist.copy_from_slice(&block[n..]); + let mut coefs = vec![0f32; n]; + m.forward(&block, &mut coefs); + // Shift: the previous folded tail (at n..n+OVERLAP/2) moves to 0..OVERLAP/2. + mem.copy_within(n..n + OVERLAP / 2, 0); + imdct_backward(2 * n, &coefs, &mut mem, &m.window); + decoded.extend_from_slice(&mem[..n]); + } + let mut err = 0f64; + let mut sum = 0f64; + for t in OVERLAP..n * frames { + let d = decoded[t] as f64 - sig[t - OVERLAP] as f64; + err += d * d; + sum += (sig[t - OVERLAP] as f64).powi(2); + } + let snr = 10.0 * (sum / err).log10(); + assert!(snr > 80.0, "reconstruction SNR {snr:.1} dB"); + } +} diff --git a/crates/opusenc/src/range.rs b/crates/opusenc/src/range.rs new file mode 100644 index 00000000..daedd58a --- /dev/null +++ b/crates/opusenc/src/range.rs @@ -0,0 +1,472 @@ +//! Range encoder (RFC 6716 §5.1), the exact inverse of the §4.1 range decoder. +//! +//! Range-coded symbols are written from the front of the buffer and raw bits from the back; the +//! bit accounting ([`RangeEncoder::tell`], [`RangeEncoder::tell_frac`]) matches the decoder's, so +//! every budget decision the decoder takes can be mirrored exactly. + +const SYM_BITS: u32 = 8; +const CODE_BITS: u32 = 32; +const SYM_MAX: u32 = (1 << SYM_BITS) - 1; +const CODE_SHIFT: u32 = CODE_BITS - SYM_BITS - 1; +const CODE_TOP: u32 = 1 << (CODE_BITS - 1); +const CODE_BOT: u32 = CODE_TOP >> SYM_BITS; +/// Bit resolution of [`RangeEncoder::tell_frac`] (1/8 bit). +pub const BITRES: u32 = 3; + +/// Number of bits needed to represent `x` (0 for 0). +#[inline] +pub fn ilog(x: u32) -> i32 { + (32 - x.leading_zeros()) as i32 +} + +#[derive(Clone, Debug)] +pub struct RangeEncoder { + out: Vec, + end: Vec, + end_window: u32, + nend_bits: i32, + nbits_total: i32, + rng: u32, + low: u32, + rem: i32, + ext: u32, +} + +impl Default for RangeEncoder { + fn default() -> Self { + Self::new() + } +} + +impl RangeEncoder { + pub fn new() -> Self { + RangeEncoder { + out: Vec::new(), + end: Vec::new(), + end_window: 0, + nend_bits: 0, + nbits_total: (CODE_BITS + 1) as i32, + rng: CODE_TOP, + low: 0, + rem: -1, + ext: 0, + } + } + + /// Current range (equals the decoder's range after the same symbols). + #[cfg(test)] + pub fn rng(&self) -> u32 { + self.rng + } + + fn carry_out(&mut self, c: u32) { + if c != SYM_MAX { + let carry = c >> SYM_BITS; + if self.rem >= 0 { + self.out.push((self.rem as u32 + carry) as u8); + } + if self.ext > 0 { + let sym = (SYM_MAX + carry) & SYM_MAX; + for _ in 0..self.ext { + self.out.push(sym as u8); + } + self.ext = 0; + } + self.rem = (c & SYM_MAX) as i32; + } else { + self.ext += 1; + } + } + + #[inline] + fn normalize(&mut self) { + while self.rng <= CODE_BOT { + self.carry_out(self.low >> CODE_SHIFT); + self.low = (self.low << SYM_BITS) & (CODE_TOP - 1); + self.rng <<= SYM_BITS; + self.nbits_total += SYM_BITS as i32; + } + } + + /// Encodes the symbol with cumulative range `[fl, fh)` of total `ft`. + pub fn encode(&mut self, fl: u32, fh: u32, ft: u32) { + let r = self.rng / ft; + if fl > 0 { + self.low += self.rng - r * (ft - fl); + self.rng = r * (fh - fl); + } else { + self.rng -= r * (ft - fh); + } + self.normalize(); + } + + /// Like [`Self::encode`] with `ft = 1 << bits`. + pub fn encode_bin(&mut self, fl: u32, fh: u32, bits: u32) { + let r = self.rng >> bits; + let ft = 1u32 << bits; + if fl > 0 { + self.low += self.rng - r * (ft - fl); + self.rng = r * (fh - fl); + } else { + self.rng -= r * (ft - fh); + } + self.normalize(); + } + + /// Encodes a binary symbol whose probability of being 1 is `1/2^logp`. + pub fn bit_logp(&mut self, val: bool, logp: u32) { + let r = self.rng; + let s = r >> logp; + let r2 = r - s; + if val { + self.low += r2; + self.rng = s; + } else { + self.rng = r2; + } + self.normalize(); + } + + /// Encodes symbol `s` with an "inverse CDF" table (`ft = 1 << ftb`). + pub fn icdf(&mut self, s: usize, icdf: &[u8], ftb: u32) { + let r = self.rng >> ftb; + if s > 0 { + self.low += self.rng - r * icdf[s - 1] as u32; + self.rng = r * (icdf[s - 1] as u32 - icdf[s] as u32); + } else { + self.rng -= r * icdf[s] as u32; + } + self.normalize(); + } + + /// Encodes a uniformly distributed integer `fl` in `[0, ft)` (`ft >= 2`). + pub fn uint(&mut self, fl: u32, ft: u32) { + debug_assert!(ft > 1 && fl < ft); + let ft1 = ft - 1; + let ftb = ilog(ft1); + if ftb > 8 { + let ftb = (ftb - 8) as u32; + let top = (ft1 >> ftb) + 1; + let v = fl >> ftb; + self.encode(v, v + 1, top); + self.bits(fl & ((1 << ftb) - 1), ftb); + } else { + self.encode(fl, fl + 1, ft); + } + } + + /// Writes `bits` raw bits (0..=25) to the end of the buffer. + pub fn bits(&mut self, val: u32, bits: u32) { + if bits == 0 { + return; + } + let mut window = self.end_window; + let mut used = self.nend_bits; + if used as u32 + bits > CODE_BITS { + while used >= SYM_BITS as i32 { + self.end.push((window & SYM_MAX) as u8); + window >>= SYM_BITS; + used -= SYM_BITS as i32; + } + } + window |= val << used; + used += bits as i32; + self.end_window = window; + self.nend_bits = used; + self.nbits_total += bits as i32; + } + + /// Number of bits written so far (rounded up), identical to the decoder's `tell`. + #[inline] + pub fn tell(&self) -> i32 { + self.nbits_total - ilog(self.rng) + } + + /// Number of 1/8 bits written so far (rounded up), identical to the decoder's `tell_frac`. + pub fn tell_frac(&self) -> i32 { + let nbits = (self.nbits_total as u32) << BITRES; + let mut l = ilog(self.rng) as u32; + let mut r = self.rng >> (l - 16); + for _ in 0..BITRES { + r = (r * r) >> 15; + let b = r >> 16; + l = (l << 1) | b; + r >>= b; + } + nbits.wrapping_sub(l) as i32 + } + + /// Finalises the stream into a buffer of exactly `size` bytes (range data at the front, raw + /// bits at the back, zero padding in between; a partial raw-bit byte is merged into the byte + /// just before the full raw-bit bytes). + pub fn done(mut self, size: usize) -> Vec { + let mut l = CODE_BITS as i32 - ilog(self.rng); + let mut msk = (CODE_TOP - 1) >> l; + let mut end = (self.low.wrapping_add(msk)) & !msk; + if (end | msk) >= self.low.wrapping_add(self.rng) { + l += 1; + msk >>= 1; + end = (self.low.wrapping_add(msk)) & !msk; + } + while l > 0 { + self.carry_out(end >> CODE_SHIFT); + end = (end << SYM_BITS) & (CODE_TOP - 1); + l -= SYM_BITS as i32; + } + if self.rem >= 0 || self.ext > 0 { + self.carry_out(0); + } + let mut window = self.end_window; + let mut used = self.nend_bits; + while used >= SYM_BITS as i32 { + self.end.push((window & SYM_MAX) as u8); + window >>= SYM_BITS; + used -= SYM_BITS as i32; + } + let mut buf = vec![0u8; size]; + let full = self.end.len().min(size); + let room = size - full; + let n = self.out.len().min(room); + buf[..n].copy_from_slice(&self.out[..n]); + for (k, &b) in self.end.iter().take(full).enumerate() { + buf[size - 1 - k] = b; + } + if used > 0 && room > 0 { + buf[room - 1] |= window as u8; + } + buf + } +} + +#[cfg(test)] +pub(crate) mod dec { + //! Minimal range decoder (RFC 6716 §4.1) used to verify the encoder. + use super::*; + const CODE_EXTRA: u32 = (CODE_BITS - 2) % SYM_BITS + 1; + + pub struct RangeDecoder<'a> { + buf: &'a [u8], + offs: usize, + end_offs: usize, + end_window: u32, + nend_bits: i32, + nbits_total: i32, + rng: u32, + val: u32, + rem: u32, + } + + impl<'a> RangeDecoder<'a> { + pub fn new(buf: &'a [u8]) -> Self { + let mut d = RangeDecoder { + buf, + offs: 0, + end_offs: 0, + end_window: 0, + nend_bits: 0, + nbits_total: (CODE_BITS + 1 - ((CODE_BITS - CODE_EXTRA) / SYM_BITS) * SYM_BITS) as i32, + rng: 1 << CODE_EXTRA, + val: 0, + rem: 0, + }; + d.rem = d.read_byte(); + d.val = d.rng - 1 - (d.rem >> (SYM_BITS - CODE_EXTRA)); + d.normalize(); + d + } + fn read_byte(&mut self) -> u32 { + if self.offs < self.buf.len() { + self.offs += 1; + self.buf[self.offs - 1] as u32 + } else { + 0 + } + } + fn read_byte_from_end(&mut self) -> u32 { + if self.end_offs < self.buf.len() { + self.end_offs += 1; + self.buf[self.buf.len() - self.end_offs] as u32 + } else { + 0 + } + } + fn normalize(&mut self) { + while self.rng <= CODE_BOT { + self.nbits_total += SYM_BITS as i32; + self.rng <<= SYM_BITS; + let sym = self.rem; + self.rem = self.read_byte(); + let sym = ((sym << SYM_BITS) | self.rem) >> (SYM_BITS - CODE_EXTRA); + self.val = ((self.val << SYM_BITS).wrapping_add(SYM_MAX & !sym)) & (CODE_TOP - 1); + } + } + pub fn rng(&self) -> u32 { + self.rng + } + pub fn decode(&mut self, ft: u32) -> u32 { + let ext = self.rng / ft; + let s = self.val / ext; + ft - (s + 1).min(ft) + } + pub fn update(&mut self, fl: u32, fh: u32, ft: u32) { + let ext = self.rng / ft; + let s = ext * (ft - fh); + self.val = self.val.wrapping_sub(s); + self.rng = if fl > 0 { ext * (fh - fl) } else { self.rng - s }; + self.normalize(); + } + pub fn bit_logp(&mut self, logp: u32) -> bool { + let r = self.rng; + let d = self.val; + let s = r >> logp; + let ret = d < s; + if !ret { + self.val = d - s; + } + self.rng = if ret { s } else { r - s }; + self.normalize(); + ret + } + pub fn icdf(&mut self, icdf: &[u8], ftb: u32) -> usize { + let mut s = self.rng; + let d = self.val; + let r = s >> ftb; + let mut ret = 0usize; + let mut t; + loop { + t = s; + s = r * icdf[ret] as u32; + if d >= s { + break; + } + ret += 1; + } + self.val = d - s; + self.rng = t - s; + self.normalize(); + ret + } + pub fn uint(&mut self, ft: u32) -> u32 { + let ft1 = ft - 1; + let ftb = ilog(ft1); + if ftb > 8 { + let ftb = (ftb - 8) as u32; + let top = (ft1 >> ftb) + 1; + let s = self.decode(top); + self.update(s, s + 1, top); + (s << ftb) | self.bits(ftb) + } else { + let s = self.decode(ft); + self.update(s, s + 1, ft); + s + } + } + pub fn bits(&mut self, bits: u32) -> u32 { + let mut window = self.end_window; + let mut available = self.nend_bits; + if (available as u32) < bits { + loop { + window |= self.read_byte_from_end() << available; + available += SYM_BITS as i32; + if available > (CODE_BITS - SYM_BITS) as i32 { + break; + } + } + } + let ret = window & ((1u32 << bits) - 1); + self.end_window = window >> bits; + self.nend_bits = available - bits as i32; + self.nbits_total += bits as i32; + ret + } + pub fn tell(&self) -> i32 { + self.nbits_total - ilog(self.rng) + } + } +} + +#[cfg(test)] +mod tests { + use super::dec::RangeDecoder; + use super::*; + + struct Lcg(u32); + impl Lcg { + fn next(&mut self) -> u32 { + self.0 = self.0.wrapping_mul(1664525).wrapping_add(1013904223); + self.0 >> 8 + } + } + + #[test] + fn starts_at_one_bit() { + let e = RangeEncoder::new(); + assert_eq!(e.tell(), 1); + assert_eq!(e.tell_frac(), 8); + } + + #[test] + fn roundtrip_mixed_symbols_and_tell_matches() { + let mut rng = Lcg(7); + let icdf: [u8; 4] = [200, 120, 30, 0]; + for trial in 0..300 { + let n = 1 + (rng.next() % 300) as usize; + let mut ops = Vec::new(); + for _ in 0..n { + match rng.next() % 5 { + 0 => { + let ft = 2 + rng.next() % 1000; + ops.push((0u8, rng.next() % ft, ft)); + } + 1 => { + let logp = 1 + rng.next() % 15; + ops.push((1, rng.next().is_multiple_of(5) as u32, logp)); + } + 2 => { + let bits = 1 + rng.next() % 24; + ops.push((2, rng.next() & ((1 << bits) - 1), bits)); + } + 3 => { + let ft = 2 + rng.next() % 100_000_000; + ops.push((3, rng.next() % ft, ft)); + } + _ => ops.push((4, rng.next() % 4, 0)), + } + } + let mut e = RangeEncoder::new(); + let mut tells = Vec::new(); + for &(k, v, p) in &ops { + match k { + 0 => e.encode(v, v + 1, p), + 1 => e.bit_logp(v != 0, p), + 2 => e.bits(v, p), + 3 => e.uint(v, p), + _ => e.icdf(v as usize, &icdf, 8), + } + tells.push(e.tell()); + } + // Exactly the number of bytes `tell` says were used must suffice. + let size = (e.tell() as usize).div_ceil(8); + let final_rng = e.rng(); + let buf = e.done(size); + assert_eq!(buf.len(), size); + let mut d = RangeDecoder::new(&buf); + for (i, &(k, v, p)) in ops.iter().enumerate() { + let got = match k { + 0 => { + let s = d.decode(p); + d.update(s, s + 1, p); + s + } + 1 => d.bit_logp(p) as u32, + 2 => d.bits(p), + 3 => d.uint(p), + _ => d.icdf(&icdf, 8) as u32, + }; + assert_eq!(got, v, "trial {trial} op {i} kind {k}"); + assert_eq!(d.tell(), tells[i], "tell mismatch trial {trial} op {i}"); + } + assert_eq!(d.rng(), final_rng); + } + } +} diff --git a/crates/opusenc/src/rate.rs b/crates/opusenc/src/rate.rs new file mode 100644 index 00000000..d8665bf1 --- /dev/null +++ b/crates/opusenc/src/rate.rs @@ -0,0 +1,561 @@ +//! Pulse cache and bit allocation (RFC 6716 §4.3.3, §4.3.4.1). +//! +//! The allocation must be bit-identical to the decoder's, so this mirrors FilmCraft's +//! spec-derived `filmcraft-opus` decoder (first-party, MIT OR Apache-2.0) step by step; the only +//! difference is that the encoder *writes* the band-skip, intensity and dual-stereo symbols. + +use crate::range::{BITRES, RangeEncoder, ilog}; +use crate::tables::*; + +const MAX_PSEUDO: i32 = 40; +const LOG_MAX_PSEUDO: i32 = 6; +pub const MAX_FINE_BITS: i32 = 8; +const FINE_OFFSET: i32 = 21; +pub const QTHETA_OFFSET: i32 = 4; +pub const QTHETA_OFFSET_TWOPHASE: i32 = 16; +const ALLOC_STEPS: i32 = 6; + +/// Number of pulses represented by pseudo-pulse index `i`. +pub fn get_pulses(i: i32) -> i32 { + if i < 8 { i } else { (8 + (i & 7)) << ((i >> 3) - 1) } +} + +/// Conservative (rounded up) log2 of `val` with `frac` fractional bits. +pub fn log2_frac(mut val: u32, mut frac: u32) -> i32 { + let mut l = ilog(val); + if val & val.wrapping_sub(1) != 0 { + if l > 16 { + val = ((val - 1) >> (l - 16)) + 1; + } else { + val <<= 16 - l; + } + let mut l = (l - 1) << frac; + loop { + let b = (val >> 16) as i32; + l += b << frac; + val = (val + b as u32) >> b; + val = (val * val + 0x7FFF) >> 15; + if frac == 0 { + break; + } + frac -= 1; + } + l + (val > 0x8000) as i32 + } else { + l -= 1; + l << frac + } +} + +/// PVQ codebook size V(N, K) (saturating far above 2^32). +pub fn pvq_v(n: usize, k: usize) -> u128 { + if k == 0 { + return 1; + } + if n == 0 { + return 0; + } + let mut row = vec![0u128; k + 1]; + row[0] = 1; + for _ in 0..n { + let mut prev_diag = row[0]; + for kk in 1..=k { + let up = row[kk]; + let v = up.saturating_add(row[kk - 1]).saturating_add(prev_diag); + prev_diag = up; + row[kk] = v.min(u128::MAX >> 8); + } + row[0] = 1; + } + row[k] +} + +fn fits_in32(n: usize, k: i32) -> bool { + pvq_v(n, k as usize) <= u32::MAX as u128 +} + +/// Precomputed per-mode data. +pub struct Mode { + /// log2 of band width in 1/8 bits. + pub log_n: [i32; NB_EBANDS], + cache_index: Vec, + cache_bits: Vec, + caps: Vec, +} + +impl Mode { + pub fn new() -> Mode { + let mut log_n = [0i32; NB_EBANDS]; + for (i, l) in log_n.iter_mut().enumerate() { + *l = log2_frac((EBANDS[i + 1] - EBANDS[i]) as u32, BITRES); + } + let (cache_index, cache_bits) = compute_pulse_cache(MAX_LM as i32); + let mut m = Mode { log_n, cache_index, cache_bits, caps: Vec::new() }; + m.caps = compute_caps(&m); + m + } + + #[inline] + pub fn cache(&self, lm: i32, band: usize) -> &[u8] { + let idx = self.cache_index[((lm + 1) as usize) * NB_EBANDS + band]; + &self.cache_bits[idx.max(0) as usize..] + } + + pub fn bits2pulses(&self, band: usize, lm: i32, bits: i32) -> i32 { + let cache = self.cache(lm, band); + let mut lo = 0i32; + let mut hi = cache[0] as i32; + let bits = bits - 1; + for _ in 0..LOG_MAX_PSEUDO { + let mid = (lo + hi + 1) >> 1; + if cache[mid as usize] as i32 >= bits { + hi = mid; + } else { + lo = mid; + } + } + let lo_bits = if lo == 0 { -1 } else { cache[lo as usize] as i32 }; + if bits - lo_bits <= cache[hi as usize] as i32 - bits { lo } else { hi } + } + + pub fn pulses2bits(&self, band: usize, lm: i32, pulses: i32) -> i32 { + if pulses == 0 { 0 } else { self.cache(lm, band)[pulses as usize] as i32 + 1 } + } + + /// Per-band allocation caps for this frame (`init_caps`). + pub fn init_caps(&self, lm: usize, c: usize) -> [i32; NB_EBANDS] { + 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; + } + cap + } +} + +fn compute_pulse_cache(lm_max: i32) -> (Vec, Vec) { + let nb = NB_EBANDS; + let rows = (lm_max + 2) as usize; + let mut index = vec![-1i16; nb * rows]; + let mut entries: Vec<(usize, i32, usize)> = Vec::new(); + let mut curr = 0usize; + for i in 0..rows { + for j in 0..nb { + let n = (((EBANDS[j + 1] - EBANDS[j]) << i) >> 1) as usize; + index[i * nb + j] = -1; + 'search: for k in 0..=i { + for nn in 0..nb { + if k == i && nn >= j { + break; + } + if n == (((EBANDS[nn + 1] - EBANDS[nn]) << k) >> 1) as usize { + index[i * nb + j] = index[k * nb + nn]; + break 'search; + } + } + } + if index[i * nb + j] == -1 && n != 0 { + let mut k = 0; + while k < MAX_PSEUDO && fits_in32(n, get_pulses(k + 1)) { + k += 1; + } + index[i * nb + j] = curr as i16; + entries.push((n, k, curr)); + curr += (k + 1) as usize; + } + } + } + let mut bits = vec![0u8; curr]; + for &(n, k, off) in &entries { + let maxk = get_pulses(k) as usize; + let tmp = required_bits(n, maxk, BITRES); + for j in 1..=k { + bits[off + j as usize] = (tmp[get_pulses(j) as usize] - 1) as u8; + } + bits[off] = k as u8; + } + (index, bits) +} + +fn required_bits(n: usize, maxk: usize, frac: u32) -> Vec { + let mut out = vec![0i32; maxk + 1]; + if n == 1 { + for o in out.iter_mut().skip(1) { + *o = 1 << frac; + } + } else { + for (k, o) in out.iter_mut().enumerate().skip(1) { + *o = log2_frac(pvq_v(n, k) as u32, frac); + } + } + out +} + +fn compute_caps(m: &Mode) -> Vec { + let nb = NB_EBANDS; + let mut caps = vec![0u8; (MAX_LM + 1) * 2 * nb]; + for i in 0..=MAX_LM as i32 { + for c in 1..=2i32 { + for j in 0..nb { + let mut n0 = EBANDS[j + 1] - EBANDS[j]; + let max_bits; + if n0 << i == 1 { + max_bits = c * (1 + MAX_FINE_BITS) << BITRES; + } else { + let mut lm0 = 0; + if n0 > 2 { + n0 >>= 1; + lm0 -= 1; + } else if n0 <= 1 { + lm0 = i.min(1); + n0 <<= lm0; + } + let pcache = m.cache(lm0, j); + let mut mb = pcache[pcache[0] as usize] as i32 + 1; + let mut n = n0; + for k in 0..(i - lm0) { + mb <<= 1; + let offset = ((m.log_n[j] + ((lm0 + k) << BITRES)) >> 1) - QTHETA_OFFSET; + let num = 459 * ((2 * n - 1) * offset + mb); + let den = ((2 * n - 1) << 9) - 459; + let qb = ((num + (den >> 1)) / den).min(57); + mb += qb; + n <<= 1; + } + if c == 2 { + mb <<= 1; + let offset = ((m.log_n[j] + (i << BITRES)) >> 1) - if n == 2 { QTHETA_OFFSET_TWOPHASE } else { QTHETA_OFFSET }; + let ndof = 2 * n - 1 - (n == 2) as i32; + let p = if n == 2 { 512 } else { 487 }; + let num = p * (mb + ndof * offset); + let den = (ndof << 9) - p; + let qb = ((num + (den >> 1)) / den).min(if n == 2 { 64 } else { 61 }); + mb += qb; + } + let ndof = c * n + if c == 2 && n > 2 { 1 } else { 0 }; + let mut offset = ((m.log_n[j] + (i << BITRES)) >> 1) - FINE_OFFSET; + if n == 2 { + offset += (1 << BITRES) >> 2; + } + let num = mb + ndof * offset; + let den = (ndof - 1) << BITRES; + let qb = ((num + (den >> 1)) / den).min(MAX_FINE_BITS); + mb += c * qb << BITRES; + max_bits = mb; + } + let w = c * ((EBANDS[j + 1] - EBANDS[j]) << i); + let v = (4 * max_bits / w) - 64; + caps[((i * 2 + c - 1) as usize) * nb + j] = v.clamp(0, 255) as u8; + } + } + } + caps +} + +/// Result of the allocation process. +pub struct Allocation { + pub coded_bands: usize, + pub intensity: usize, + pub dual_stereo: bool, + pub balance: i32, + pub pulses: [i32; NB_EBANDS], + pub fine_quant: [i32; NB_EBANDS], + pub fine_priority: [i32; NB_EBANDS], +} + +/// `clt_compute_allocation`, encoder side. Stereo frames are always coded as dual (L/R) stereo +/// with no intensity bands, and no band is skipped voluntarily. +#[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, + enc: &mut RangeEncoder, +) -> Allocation { + let ci = c as i32; + let lmi = lm as i32; + let mut total = total.max(0); + let mut skip_start = start; + let skip_rsv = if total >= 1 << BITRES { 1 << BITRES } else { 0 }; + total -= skip_rsv; + let mut intensity_rsv = 0; + let mut dual_stereo_rsv = 0; + if c == 2 { + intensity_rsv = LOG2_FRAC_TABLE[end - start]; + if intensity_rsv > total { + intensity_rsv = 0; + } else { + total -= intensity_rsv; + dual_stereo_rsv = if total >= 1 << BITRES { 1 << BITRES } else { 0 }; + total -= dual_stereo_rsv; + } + } + let mut bits1 = [0i32; NB_EBANDS]; + let mut bits2 = [0i32; NB_EBANDS]; + let mut thresh = [0i32; NB_EBANDS]; + 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; + if w << lmi == 1 { + trim_offset[j] -= ci << BITRES; + } + } + let nb_alloc = ALLOC.len() as i32; + let mut lo = 1; + let mut hi = nb_alloc - 1; + loop { + let mut done = false; + let mut psum = 0; + 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; + if bitsj > 0 { + bitsj = (bitsj + trim_offset[j]).max(0); + } + bitsj += offsets[j]; + if bitsj >= thresh[j] || done { + done = true; + psum += bitsj.min(cap[j]); + } else if bitsj >= ci << BITRES { + psum += ci << BITRES; + } + } + if psum > total { + hi = mid - 1; + } else { + lo = mid + 1; + } + if lo > hi { + break; + } + } + hi = lo; + 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 }; + if bits1j > 0 { + bits1j = (bits1j + trim_offset[j]).max(0); + } + if bits2j > 0 { + bits2j = (bits2j + trim_offset[j]).max(0); + } + if lo > 0 { + bits1j += offsets[j]; + } + bits2j += offsets[j]; + if offsets[j] > 0 { + skip_start = j; + } + bits2j = (bits2j - bits1j).max(0); + 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, enc) +} + +#[allow(clippy::too_many_arguments)] +fn interp_bits2pulses( + m: &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, + lm: usize, + enc: &mut RangeEncoder, +) -> Allocation { + let ci = c as i32; + let alloc_floor = ci << BITRES; + let stereo = (c > 1) as i32; + let log_m = (lm as i32) << BITRES; + let mut lo = 0; + let mut hi = 1 << ALLOC_STEPS; + for _ in 0..ALLOC_STEPS { + let mid = (lo + hi) >> 1; + let mut psum = 0; + let mut done = false; + for j in (start..end).rev() { + let tmp = bits1[j] + (mid * bits2[j] >> ALLOC_STEPS); + if tmp >= thresh[j] || done { + done = true; + psum += tmp.min(cap[j]); + } else if tmp >= alloc_floor { + psum += alloc_floor; + } + } + if psum > total { + hi = mid; + } else { + lo = mid; + } + } + let mut bits = [0i32; NB_EBANDS]; + let mut psum = 0; + let mut done = false; + for j in (start..end).rev() { + 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 { + done = true; + } + tmp = tmp.min(cap[j]); + bits[j] = tmp; + psum += tmp; + } + let s = start; + let mut coded_bands = end; + loop { + let j = coded_bands - 1; + if j <= skip_start { + total += skip_rsv; + break; + } + let mut left = total - psum; + let span = EBANDS[coded_bands] - EBANDS[s]; + let percoeff = left / span; + left -= span * percoeff; + let rem = (left - (EBANDS[j] - EBANDS[s])).max(0); + let band_width = EBANDS[coded_bands] - EBANDS[j]; + let band_bits = bits[j] + percoeff * band_width + rem; + if band_bits >= thresh[j].max(alloc_floor + (1 << BITRES)) { + // Keep this band (and every band below it). + enc.bit_logp(true, 1); + break; + } + psum -= bits[j] + intensity_rsv; + if intensity_rsv > 0 { + intensity_rsv = LOG2_FRAC_TABLE[j - s]; + } + psum += intensity_rsv; + if band_bits >= alloc_floor { + psum += alloc_floor; + bits[j] = alloc_floor; + } else { + bits[j] = 0; + } + coded_bands -= 1; + } + let intensity = if intensity_rsv > 0 { + // No intensity stereo: every coded band is coded in full. + enc.uint((coded_bands - s) as u32, (coded_bands + 1 - s) as u32); + coded_bands + } else { + 0 + }; + if intensity <= s { + total += dual_stereo_rsv; + dual_stereo_rsv = 0; + } + let dual_stereo = if dual_stereo_rsv > 0 { + enc.bit_logp(true, 1); + true + } else { + false + }; + let mut left = total - psum; + let span = EBANDS[coded_bands] - EBANDS[s]; + let percoeff = left / span; + left -= span * percoeff; + for j in s..coded_bands { + bits[j] += percoeff * (EBANDS[j + 1] - EBANDS[j]); + } + for j in s..coded_bands { + let tmp = left.min(EBANDS[j + 1] - EBANDS[j]); + bits[j] += tmp; + left -= tmp; + } + let mut ebits = [0i32; NB_EBANDS]; + let mut fine_priority = [0i32; NB_EBANDS]; + let mut balance = 0; + for j in s..coded_bands { + let n0 = EBANDS[j + 1] - EBANDS[j]; + let n = n0 << lm; + let bit = bits[j] + balance; + let mut excess; + if n > 1 { + excess = (bit - cap[j]).max(0); + bits[j] = bit - excess; + let den = ci * n + if c == 2 && n > 2 && !dual_stereo && j < intensity { 1 } else { 0 }; + let nclogn = den * (m.log_n[j] + log_m); + let mut offset = (nclogn >> 1) - den * FINE_OFFSET; + if n == 2 { + offset += (den << BITRES) >> 2; + } + if bits[j] + offset < den * 2 << BITRES { + offset += nclogn >> 2; + } else if bits[j] + offset < den * 3 << BITRES { + offset += nclogn >> 3; + } + let mut eb = (bits[j] + offset + (den << (BITRES - 1))).max(0); + eb = (eb / den) >> BITRES; + if ci * eb > (bits[j] >> BITRES) { + eb = bits[j] >> stereo >> BITRES; + } + 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; + } else { + excess = (bit - (ci << BITRES)).max(0); + bits[j] = bit - excess; + ebits[j] = 0; + fine_priority[j] = 1; + } + 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; + fine_priority[j] = (extra_bits >= excess - balance) as i32; + excess -= extra_bits; + } + balance = excess; + } + for j in coded_bands..end { + ebits[j] = bits[j] >> stereo >> BITRES; + bits[j] = 0; + fine_priority[j] = (ebits[j] < 1) as i32; + } + Allocation { coded_bands, intensity, dual_stereo, balance, pulses: bits, fine_quant: ebits, fine_priority } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn caps_match_decoder_reference() { + let m = Mode::new(); + let first_rows: [u8; 42] = [ + 224, 224, 224, 224, 224, 224, 224, 224, 160, 160, 160, 160, 185, 185, 185, 178, 178, 168, 134, 61, 37, 224, 224, 224, 224, 224, 224, 224, 224, 240, + 240, 240, 240, 207, 207, 207, 198, 198, 183, 144, 66, 40, + ]; + assert_eq!(&m.caps[..42], &first_rows[..]); + } + + #[test] + fn pvq_counts() { + assert_eq!(pvq_v(5, 1), 10); + assert_eq!(pvq_v(1, 7), 2); + assert_eq!(pvq_v(2, 9), 36); + assert_eq!(pvq_v(3, 2), 18); + } +} diff --git a/crates/opusenc/src/tables.rs b/crates/opusenc/src/tables.rs new file mode 100644 index 00000000..00d623af --- /dev/null +++ b/crates/opusenc/src/tables.rs @@ -0,0 +1,122 @@ +//! Constant tables of the 48 kHz CELT mode (RFC 6716 §4.3, Tables 55–58 and §4.3.2). +//! +//! Transcribed from RFC 6716; copied from FilmCraft's spec-derived `filmcraft-opus` decoder +//! (first-party, MIT OR Apache-2.0) so both sides share the exact same values. + +/// Number of energy bands. +pub const NB_EBANDS: usize = 21; +/// Band edges in units of 2.5 ms MDCT bins (Table 55). +pub const EBANDS: [i32; NB_EBANDS + 1] = [0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 20, 24, 28, 34, 40, 48, 60, 78, 100]; +/// Overlap of the low-overlap window (samples at 48 kHz). +pub const OVERLAP: usize = 120; +/// Largest LM (20 ms frames). +pub const MAX_LM: usize = 3; +/// Pre-emphasis coefficient. +pub const PREEMPH: f32 = 0.850_006_1; + +/// Static allocation table (Table 57), indexed `[quality][band]`, in 1/32 bit per MDCT bin. +pub const ALLOC: [[u8; NB_EBANDS]; 11] = transpose([ + [0, 90, 110, 118, 126, 134, 144, 152, 162, 172, 200], + [0, 80, 100, 110, 119, 127, 137, 145, 155, 165, 200], + [0, 75, 90, 103, 112, 120, 130, 138, 148, 158, 200], + [0, 69, 84, 93, 104, 114, 124, 132, 142, 152, 200], + [0, 63, 78, 86, 95, 103, 113, 123, 133, 143, 200], + [0, 56, 71, 80, 89, 97, 107, 117, 127, 137, 200], + [0, 49, 65, 75, 83, 91, 101, 111, 121, 131, 200], + [0, 40, 58, 70, 78, 85, 95, 105, 115, 125, 200], + [0, 34, 51, 65, 72, 78, 88, 98, 108, 118, 198], + [0, 29, 45, 59, 66, 72, 82, 92, 102, 112, 193], + [0, 20, 39, 53, 60, 66, 76, 86, 96, 106, 188], + [0, 18, 32, 47, 54, 60, 70, 80, 90, 100, 183], + [0, 10, 26, 40, 47, 54, 64, 74, 84, 94, 178], + [0, 0, 20, 31, 39, 47, 57, 67, 77, 87, 173], + [0, 0, 12, 23, 32, 41, 51, 61, 71, 81, 168], + [0, 0, 0, 15, 25, 35, 45, 55, 65, 75, 163], + [0, 0, 0, 4, 17, 29, 39, 49, 59, 69, 158], + [0, 0, 0, 0, 12, 23, 33, 43, 53, 63, 153], + [0, 0, 0, 0, 1, 16, 26, 36, 46, 56, 148], + [0, 0, 0, 0, 0, 10, 15, 20, 30, 45, 129], + [0, 0, 0, 0, 0, 1, 1, 1, 1, 20, 104], +]); + +const fn transpose(t: [[u8; 11]; NB_EBANDS]) -> [[u8; NB_EBANDS]; 11] { + let mut out = [[0u8; NB_EBANDS]; 11]; + let mut b = 0; + while b < NB_EBANDS { + let mut q = 0; + while q < 11 { + out[q][b] = t[b][q]; + q += 1; + } + b += 1; + } + out +} + +/// Mean band energy (log2 domain) removed before coarse quantisation. +pub const E_MEANS: [f32; 25] = [ + 6.4375, 6.25, 5.75, 5.3125, 5.0625, 4.8125, 4.5, 4.375, 4.875, 4.6875, 4.5625, 4.4375, 4.875, 4.625, 4.3125, 4.5, 4.375, 4.625, 4.75, 4.4375, 3.75, 3.75, + 3.75, 3.75, 3.75, +]; + +/// Inter-frame prediction coefficient (alpha) per LM, Q15. +pub const PRED_COEF: [f32; 4] = [29440.0 / 32768.0, 26112.0 / 32768.0, 21248.0 / 32768.0, 16384.0 / 32768.0]; +/// Inter-band prediction coefficient (beta) per LM for inter frames, Q15. +pub const BETA_COEF: [f32; 4] = [30147.0 / 32768.0, 22282.0 / 32768.0, 12124.0 / 32768.0, 6554.0 / 32768.0]; +/// Beta for intra frames. +pub const BETA_INTRA: f32 = 4915.0 / 32768.0; + +/// Laplace model parameters for coarse energy: `[LM][intra][2*band + {0: prob, 1: decay}]`. +pub const E_PROB_MODEL: [[[u8; 42]; 2]; 4] = [ + [ + [ + 72, 127, 65, 129, 66, 128, 65, 128, 64, 128, 62, 128, 64, 128, 64, 128, 92, 78, 92, 79, 92, 78, 90, 79, 116, 41, 115, 40, 114, 40, 132, 26, 132, + 26, 145, 17, 161, 12, 176, 10, 177, 11, + ], + [ + 24, 179, 48, 138, 54, 135, 54, 132, 53, 134, 56, 133, 55, 132, 55, 132, 61, 114, 70, 96, 74, 88, 75, 88, 87, 74, 89, 66, 91, 67, 100, 59, 108, 50, + 120, 40, 122, 37, 97, 43, 78, 50, + ], + ], + [ + [ + 83, 78, 84, 81, 88, 75, 86, 74, 87, 71, 90, 73, 93, 74, 93, 74, 109, 40, 114, 36, 117, 34, 117, 34, 143, 17, 145, 18, 146, 19, 162, 12, 165, 10, + 178, 7, 189, 6, 190, 8, 177, 9, + ], + [ + 23, 178, 54, 115, 63, 102, 66, 98, 69, 99, 74, 89, 71, 91, 73, 91, 78, 89, 86, 80, 92, 66, 93, 64, 102, 59, 103, 60, 104, 60, 117, 52, 123, 44, + 138, 35, 133, 31, 97, 38, 77, 45, + ], + ], + [ + [ + 61, 90, 93, 60, 105, 42, 107, 41, 110, 45, 116, 38, 113, 38, 112, 38, 124, 26, 132, 27, 136, 19, 140, 20, 155, 14, 159, 16, 158, 18, 170, 13, 177, + 10, 187, 8, 192, 6, 175, 9, 159, 10, + ], + [ + 21, 178, 59, 110, 71, 86, 75, 85, 84, 83, 91, 66, 88, 73, 87, 72, 92, 75, 98, 72, 105, 58, 107, 54, 115, 52, 114, 55, 112, 56, 129, 51, 132, 40, + 150, 33, 140, 29, 98, 35, 77, 42, + ], + ], + [ + [ + 42, 121, 96, 66, 108, 43, 111, 40, 117, 44, 123, 32, 120, 36, 119, 33, 127, 33, 134, 34, 139, 21, 147, 23, 152, 20, 158, 25, 154, 26, 166, 21, 173, + 16, 184, 13, 184, 10, 150, 13, 139, 15, + ], + [ + 22, 178, 63, 114, 74, 82, 84, 83, 92, 82, 103, 62, 96, 72, 96, 67, 101, 73, 107, 72, 113, 55, 118, 52, 125, 52, 118, 52, 117, 55, 135, 49, 137, 39, + 157, 32, 145, 29, 97, 33, 77, 40, + ], + ], +]; + +/// Inverse CDF for coarse energy when only a few bits remain. +pub const SMALL_ENERGY_ICDF: [u8; 3] = [2, 1, 0]; +/// Conservative log2 in 1/8 bits of the number of coded bands (intensity reservation). +pub const LOG2_FRAC_TABLE: [i32; 24] = [0, 8, 13, 16, 19, 21, 23, 24, 26, 27, 28, 29, 30, 31, 32, 32, 33, 34, 34, 35, 36, 36, 37, 37]; +/// TF resolution changes: `[LM][4*transient + 2*tf_select + tf_flag]`. +pub const TF_SELECT_TABLE: [[i8; 8]; 4] = [[0, -1, 0, -1, 0, -1, 0, -1], [0, -1, 0, -2, 1, 0, 1, -1], [0, -2, 0, -3, 2, 0, 1, -1], [0, -2, 0, -3, 3, 0, 1, -1]]; +/// Spreading decision `{7, 2, 21, 2}/32`. +pub const SPREAD_ICDF: [u8; 4] = [25, 23, 2, 0]; +/// Allocation trim (Table 58). +pub const TRIM_ICDF: [u8; 11] = [126, 124, 119, 109, 87, 41, 19, 9, 4, 2, 0]; diff --git a/docs/architecture.md b/docs/architecture.md index 033ed6cb..c3d869c0 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -16,6 +16,7 @@ cpal or muda. Everything in L0 to L4, the egui UI and the web app also build for | L0 | `time` | `Tick` (254 016 000 000 per second), rational frame rates incl. NTSC, SMPTE and drop-frame timecode | | L0 | `geom` | Vectors, matrices, quaternions, the layer transform (anchor, position, scale, orientation, rotation) | | L0 | `color` | sRGB and linear, HSL/HSV, luminance, the 38 blend modes, label colors | +| L0 | `opusenc` | Opus (CELT) encoder: RFC 6716 CELT-only fullband 48 kHz 20 ms packets, mono/stereo, and the RFC 7845 `OpusHead`, for WebM export audio | | L1 | `raster` | Premultiplied float images, sampling, affine and projective warps, blurs, compositing (parallel with rayon) | | L1 | `keyframe` | Animated values, keyframes with temporal ease and spatial Bezier, roving, hold, velocity | | L1 | `path` | Bezier paths, path operators (trim, offset, round corners, zig zag, twist, merge…), stroking, coverage masks | diff --git a/xtask/src/main.rs b/xtask/src/main.rs index 09e15047..7e7616a5 100644 --- a/xtask/src/main.rs +++ b/xtask/src/main.rs @@ -24,6 +24,7 @@ const LAYERS: &[(&str, u8)] = &[ ("geom", 0), ("color", 0), ("testkit", 0), + ("opusenc", 0), ("raster", 1), ("keyframe", 1), ("path", 1),