OUT-2: Opus (CELT) encoder (effectcraft-opusenc)

New L0 crate effectcraft-opusenc: clean-room RFC 6716 CELT-only encoder
(fullband 48 kHz, 20 ms packets, mono/stereo as dual L/R stereo, constant
packet size from the bitrate) with the RFC 7845 OpusHead and pre-skip.
Forward low-overlap MDCT, inter/intra coarse energy, fine/final energy,
the decoder's exact allocation, recursive band splitting, PVQ search and
codeword indexing, silence frames.

crates/export depends on it; filmcraft-opus (dev) is the round-trip oracle
in tests/opus_roundtrip.rs, with ffmpeg reading an in-test Ogg Opus file
as a second oracle when installed.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Brandon Thomas 2026-10-02 21:41:29 -04:00
parent e911cc5c93
commit b936ed0055
15 changed files with 2911 additions and 0 deletions

6
Cargo.lock generated
View File

@ -1218,6 +1218,7 @@ version = "0.1.0"
dependencies = [
"effectcraft-color",
"effectcraft-media",
"effectcraft-opusenc",
"effectcraft-project",
"effectcraft-raster",
"effectcraft-render",
@ -1225,6 +1226,7 @@ dependencies = [
"filmcraft-aac",
"filmcraft-h264enc",
"filmcraft-isobmff",
"filmcraft-opus",
"filmcraft-prores",
"gif",
"image",
@ -1338,6 +1340,10 @@ dependencies = [
"thiserror 2.0.21",
]
[[package]]
name = "effectcraft-opusenc"
version = "0.1.0"
[[package]]
name = "effectcraft-path"
version = "0.1.0"

View File

@ -32,6 +32,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"
@ -66,6 +67,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"] }
@ -129,6 +132,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

View File

@ -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 }

View File

@ -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<Vec<u8>>) {
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<u8>], len: usize) -> Vec<Vec<f32>> {
let mut dec = filmcraft_opus::Decoder::new(head).expect("OpusHead");
dec.set_trim_pre_skip(true);
let mut out: Vec<Vec<f32>> = 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<f32>], bitrate: u32) -> (Vec<Vec<u8>>, Vec<Vec<f32>>) {
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<f32> {
(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<f32> {
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<f32> {
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<f64> {
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<f64> = (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<f32> = x.iter().zip(&d[0]).map(|(a, b)| a - b).collect();
let ee = band_energies(&err);
let octave_snr: Vec<f64> = 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<f32> = 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<u8>, 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<u8>], len: usize) -> Vec<u8> {
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<Vec<f32>> =
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<f32> = 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<f32> = (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})");
}
}
}

12
crates/opusenc/Cargo.toml Normal file
View File

@ -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]

33
crates/opusenc/README.md Normal file
View File

@ -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.

497
crates/opusenc/src/bands.rs Normal file
View File

@ -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<f32> = 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<u64>) {
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<u64>,
iy: Vec<i32>,
tmp: Vec<f32>,
}
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::<f32>().sqrt();
let side = xb.iter().map(|v| v * v).sum::<f32>().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<f32> = (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<i32> {
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<Vec<i32>> {
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<f32> = (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::<usize>(), 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}");
}
}
}

View File

@ -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<i32> = (-64..=64).collect();
let mut e = RangeEncoder::new();
let coded: Vec<i32> = 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}");
}
}
}
}
}
}

311
crates/opusenc/src/lib.rs Normal file
View File

@ -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<Mode> = 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<f32>,
/// Quantised band energies, exactly as the decoder holds them.
old_band_e: [f32; 2 * NB_EBANDS],
block: Vec<f32>,
freq: Vec<f32>,
}
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<u8> {
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<u8> {
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::<f32>()).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<u8> {
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<f32> = (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);
}
}

293
crates/opusenc/src/mdct.rs Normal file
View File

@ -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<usize>,
twiddle: Vec<Cpx>,
}
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<Cpx>) {
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<Cpx>,
post: Vec<Cpx>,
/// Rising half of the low-overlap window (RFC 6716 §4.3.7).
pub window: [f32; OVERLAP],
scale: f32,
z: Vec<Cpx>,
scratch: Vec<Cpx>,
v: Vec<f32>,
}
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<f64> {
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::<f64>() * 2.0
/ n as f64
})
.collect()
}
#[test]
fn matches_direct_formula() {
for &n in &[240usize, 960] {
let block: Vec<f32> = (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<f32> = (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<f32> = (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");
}
}

472
crates/opusenc/src/range.rs Normal file
View File

@ -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<u8>,
end: Vec<u8>,
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<u8> {
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);
}
}
}

561
crates/opusenc/src/rate.rs Normal file
View File

@ -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<i16>,
cache_bits: Vec<u8>,
caps: Vec<u8>,
}
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<i16>, Vec<u8>) {
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<i32> {
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<u8> {
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);
}
}

View File

@ -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];

View File

@ -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 |

View File

@ -24,6 +24,7 @@ const LAYERS: &[(&str, u8)] = &[
("geom", 0),
("color", 0),
("testkit", 0),
("opusenc", 0),
("raster", 1),
("keyframe", 1),
("path", 1),