craft-audio: Opus encoder

EffectCraft's `effectcraft-opusenc` as `opus::encoder` (and
`opus::Encoder`) behind `opus-encode`: SILK, hybrid and CELT, 48 kHz,
20 ms frames, mono or stereo. The constructors and `encode_float`
return `Result` instead of panicking. `Bandwidth` is one type for the
decoder and the encoder, and the two `Mode` enums convert into each
other. `encode` now turns on `opus-encode`.

The encoder tests (modes, ffmpeg's two Opus decoders as oracles) and
`TrackDecoder` on Opus without an `OpusHead` come with it.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
alehelandClaude Opus 5.5 authored and Dylan M. Taylor committed 2026-11-13 21:09:46 -05:00
1 parent dc200ff58c
commit bdf8cb2a19
25 files changed
+5871 -25

No files matched your search

+1 -1
View File
@@ -44,7 +44,7 @@ jobs:
# Each feature on its own, so code that needs another feature fails here.
- name: Clippy, each feature alone
run: |
for f in "" decode decode-native symphonia encode lossless-encode pcm-encode flac-encode peaks aac-decode aac-encode "aac-encode,rayon" opus-decode ac3-decode mpa-decode serde; do
for f in "" decode decode-native symphonia encode lossless-encode pcm-encode flac-encode peaks aac-decode aac-encode "aac-encode,rayon" opus-decode opus-encode ac3-decode mpa-decode serde; do
echo "::group::features: ${f:-none}"
cargo clippy -p craft-audio --all-targets --locked --no-default-features ${f:+--features "$f"} -- -D warnings
echo "::endgroup::"
+4 -3
View File
@@ -1,6 +1,6 @@
[package]
name = "craft-audio"
description = "Audio file I/O and codecs: WAV/AIFF read and write, FLAC encoder, AAC and ALAC decoders and encoders, AC-3/E-AC-3 decoder and encoder, MP3 decoder and encoder, Opus and MP1/MP2 decoders."
description = "Audio file I/O and codecs: WAV/AIFF read and write, FLAC encoder, AAC and ALAC decoders and encoders, AC-3/E-AC-3 decoder and encoder, MP3 decoder and encoder, Opus decoder and encoder, MP1/MP2 decoders."
version.workspace = true
edition.workspace = true
license.workspace = true
@@ -22,8 +22,8 @@ decode = ["decode-native", "symphonia"]
decode-native = ["aac-decode", "alac-decode", "opus-decode", "ac3-decode", "mpa-decode", "mp3-decode", "craft-media/isobmff", "craft-media/matroska", "craft-media/ogg", "dep:log"]
# The symphonia fallback for the formats with no native decoder yet.
symphonia = ["decode-native", "dep:symphonia"]
# Every encoder: the lossless file writers and the AAC, ALAC, AC-3/E-AC-3 and MP3 encoders.
encode = ["lossless-encode", "aac-encode", "alac-encode", "ac3-encode", "eac3-encode", "mp3-encode"]
# Every encoder: the lossless file writers and the AAC, ALAC, AC-3/E-AC-3, MP3 and Opus encoders.
encode = ["lossless-encode", "aac-encode", "alac-encode", "ac3-encode", "eac3-encode", "mp3-encode", "opus-encode"]
# Encode WAV/BWF/RF64, AIFF/AIFF-C and FLAC files.
lossless-encode = ["pcm-encode", "flac-encode"]
# Encode WAV/BWF/RF64 and AIFF/AIFF-C files.
@@ -34,6 +34,7 @@ flac-encode = ["craft-media/bitstream"]
aac-decode = ["craft-media/bitstream"]
aac-encode = ["craft-media/bitstream"]
opus-decode = []
opus-encode = []
ac3-decode = ["craft-media/bitstream"]
mpa-decode = ["craft-media/bitstream"]
# MPEG-1/2 Layer III (MP3) decoder (`mp3::Decoder`).
+44 -7
View File
@@ -1,7 +1,8 @@
# craft-audio
Audio file I/O and codecs for the Craft apps. Native implementations: WAV and AIFF read/write,
FLAC encoder, AAC and ALAC decoders/encoders, AC-3/E-AC-3 decoder/encoder, MP3 decoder/encoder, Opus and MP1/MP2 decoders.
FLAC encoder, AAC and ALAC decoders/encoders, AC-3/E-AC-3 decoder/encoder, MP3 decoder/encoder, Opus
decoder/encoder, MP1/MP2 decoders.
Formats without a native decoder are decoded by [symphonia](#symphonia), which is to be replaced
by native decoders.
@@ -28,14 +29,14 @@ No default features. The shared types (`AudioBuffer`, `AudioInfo`, `ChannelLayou
| `decode` | `decode`, `probe`, `Decoder`, `tracks`, `TrackDecoder`; every format in [Format support](#format-support) | `decode-native`, `symphonia` |
| `decode-native` | the same, native formats only (others are `Unsupported`); `wav_pcm`, `aiff_pcm` (a file's integer PCM bit for bit, as `RawPcm`) | the six decoders below, craft-media `isobmff`/`matroska`/`ogg`, `log` |
| `symphonia` | the fallback for formats with no native decoder; turning it on is the same as `decode` | `decode-native`, symphonia |
| `encode` | every encoder: the WAV, AIFF and FLAC file writers, and the AAC, ALAC, AC-3/E-AC-3 and MP3 encoders | `lossless-encode`, `aac-encode`, `alac-encode`, `ac3-encode`, `eac3-encode`, `mp3-encode` |
| `encode` | every encoder: the WAV, AIFF and FLAC file writers, and the AAC, ALAC, AC-3/E-AC-3, MP3 and Opus encoders | `lossless-encode`, `aac-encode`, `alac-encode`, `ac3-encode`, `eac3-encode`, `mp3-encode`, `opus-encode` |
| `lossless-encode` | WAV, AIFF and FLAC writing | `pcm-encode`, `flac-encode` |
| `pcm-encode` | `encode`, `FileEncoder` for WAV and AIFF; `wav_float_header`; `pcm_wav_file`, `pcm_wav_header`, `pcm_aiff_file`, `pcm_aiff_header` (minimal files around integer PCM bytes, for essence embedded in AAF / OMF) | none |
| `flac-encode` | `encode`, `FileEncoder` for FLAC; `flac::Encoder`, the packet encoder for callers that mux FLAC themselves | craft-media `bitstream` |
| `aac-decode`, `aac-encode` | `aac::Decoder`, `aac::Encoder`; with `aac-decode`, `aac::latm_config` / `latm_split` for LOAS / LATM frames (AAC in MPEG transport streams) | craft-media `bitstream` |
| `alac-decode` | `alac::Decoder` | none |
| `alac-encode` | `alac::Encoder` | none |
| `opus-decode` | `opus::Decoder`/`StreamDecoder` | none |
| `opus-decode`, `opus-encode` | `opus::Decoder`/`StreamDecoder`, `opus::Encoder` | none |
| `ac3-decode` | `ac3::Decoder`, `ac3::ProgrammeDecoder` | craft-media `bitstream` |
| `ac3-encode`, `eac3-encode` | `ac3::encoder::Encoder` | craft-media `bitstream` |
| `mpa-decode` | `mpa::Decoder` | craft-media `bitstream` |
@@ -264,6 +265,38 @@ Ls Rs) has no equivalent in that layout; passed unchanged, its channels go to th
are a playable `.aac` file. To write MP4 or Matroska, pass the access units to craft-media's
muxers ([MP4](../media/README.md#muxing), [Matroska](../media/README.md#muxing-1)).
### Encode Opus
```rust
use craft_audio::opus::Encoder;
let mut enc = Encoder::new(2, 128_000)?; // channels, bits per second
let frame_len = Encoder::FRAME_SIZE * 2; // 20 ms, interleaved stereo
let mut samples = buf.interleaved(); // 48 kHz audio
// The last frame is padded with silence.
let padding = (frame_len - samples.len() % frame_len) % frame_len;
samples.resize(samples.len() + padding, 0.0);
let mut packets = Vec::new();
for frame in samples.chunks_exact(frame_len) {
packets.push(enc.encode_float(frame)?);
}
let head = enc.opus_head(); // Matroska's CodecPrivate
let pre_skip = enc.pre_skip(); // samples the player drops at the start
```
The encoder takes 48 kHz only, and craft-audio has no resampler.
craft-media has no Ogg muxer, so `.opus` files cannot be written yet; Opus goes into
Matroska/WebM or MP4:
- **Matroska/WebM:** `head` is the track's `CodecPrivate`, and the pre-skip goes in
`TrackSpec::codec_delay_ns` (in nanoseconds at 48 kHz). The Matroska muxer cannot write
`DiscardPadding`, so the padding at the end stays in the file.
- **MP4:** the sample entry takes an `isobmff::OpusConfig` (`dOps`) with the same fields as the
`OpusHead`: channels, pre-skip, input rate, gain and mapping family. `TrackConfig::media_start`
set to the pre-skip skips it; an explicit `TrackConfig::edits` with the audio's length also
drops the padding.
### Waveform overviews
```rust
@@ -314,7 +347,7 @@ match craft_audio::decode(&bytes, None) {
| ALAC | native | native, packets (`alac::Encoder`; no file writer) |
| MP1, MP2 (MPEG audio Layer I/II) | native | — |
| MP3 (Layer III) | native | native, frames (`mp3::Encoder`; the frames in a row are an `.mp3` file) |
| Opus | native | — |
| Opus | native | native |
| AC-3, E-AC-3 | native | native packet encoder |
| Vorbis, ADPCM, A-law, µ-law | symphonia | — |
@@ -415,7 +448,7 @@ Each codec module's rustdoc lists the coding tools it implements.
| WAV, AIFF | Microsoft RIFF, Apple AIFF/AIFF-C | read: integer 8–32-bit, float 32/64-bit, `WAVE_FORMAT_EXTENSIBLE` masks, RF64/BW64, BWF `bext`, truncated files; write: 16/24/32-bit integer and 32-bit float, RF64 above 4 GiB, any layout, TPDF dither | round trips; hostile-input tests |
| FLAC encoder | RFC 9639 | 16, 24, 32 bits; 1–8 channels, layouts decoders would not assume for the count as a `WAVEFORMATEXTENSIBLE_CHANNEL_MASK` comment; up to 655,350 Hz; LPC, stereo decorrelation, wasted bits, STREAMINFO MD5 | bit-exact through ffmpeg, STREAMINFO MD5 checked against ffmpeg |
| AAC (`aac`) | ISO/IEC 14496-3; SBR tables from ETSI ES 201 980 Annex I | AAC-LC decode and encode (CBR, VBR, ADTS), 1–8 channels; HE-AAC (SBR) decode, see [Limits](#limits) | ffmpeg: ~135 dB SNR; PNS noise differs per sample at equal band energy. HE-AAC from AudioToolbox, 1–6 channels: same 962-sample delay as ffmpeg, ~54 dB SNR, high band level within 0.01 dB |
| Opus (`opus`) | RFC 6716, 8251, 7845 | decoder: SILK, hybrid, CELT, all frame sizes, PLC and FEC, multistream (`Decoder`) and single stream (`StreamDecoder`), 8–48 kHz output | all 36 RFC 8251 conformance cases |
| Opus (`opus`) | RFC 6716, 8251, 7845 | decoder: SILK, hybrid, CELT, all frame sizes, PLC and FEC, multistream (`Decoder`) and single stream (`StreamDecoder`), 8–48 kHz output; encoder (`Encoder`): SILK, hybrid, CELT, 48 kHz, 20 ms frames, mono or stereo | all 36 RFC 8251 conformance cases; encoder output through ffmpeg's two Opus decoders |
| AC-3/E-AC-3 (`ac3`) | decoder: ATSC A/52:2012; encoder/programme API: ETSI TS 102 366 V1.4.1 | single-syncframe and complete-programme decoding; packet encoding through discrete 7.1 (see below) | decoder comparisons and encoder output through ffmpeg, CRC/channel/timing/bit-budget regressions; synthetic quality limitations below |
| MP1/MP2 (`mpa`) | ISO/IEC 11172-3, 13818-3 | Layer I and II, MPEG-1/2/2.5, 8–48 kHz; also the Layer II rate/mode pairs the spec disallows (ffmpeg's default mono 384 kbit/s); no free format | ffmpeg-encoded streams, 73–84 dB SNR |
| MP3 decoder (`mp3`) | ISO/IEC 11172-3, 13818-3 | Layer III, MPEG-1/2/2.5, 8–48 kHz, every bit rate, CBR and VBR; mono, stereo, dual channel, joint stereo (mid/side and intensity); long, short and mixed blocks; bit reservoir; CRC read, not checked; Xing / Info / VBRI and LAME tags; no free format (read through symphonia) | ffmpeg's float decoder: 127–131 dB SNR on LAME streams of every sampling frequency (CBR, VBR, CRC, no reservoir, transients with short blocks) and on this crate's encoder's; intensity stereo and mixed blocks through LAME frames edited to use them (112–131 dB) |
@@ -583,6 +616,7 @@ the full tests with ffmpeg and the vectors required.
| MP3 decoder | written here, from ISO/IEC 11172-3 and 13818-3; shares the encoder's tables and the MP1/MP2 synthesis filter |
| MP3 encoder | FilmCraft `filmcraft-mp3enc` [4bec22b](https://github.com/storytold/filmcraft/commit/4bec22bd03bd1d995c9cf72e02a284f96d3d72cd) |
| ALAC | Apple's reference implementation [macosforge/alac@c38887c](https://github.com/macosforge/alac/tree/c38887c5c5e64a4b31108733bd79ca9b2496d987/codec), Copyright © 2011 Apple Inc., Apache-2.0; ported to Rust and modified (see the `alac` module docs and [NOTICE](../../NOTICE)) |
| Opus encoder | EffectCraft `effectcraft-opusenc` [c618429](https://github.com/storytold/effectcraft/commit/c618429252587525a2a6d4b5bccbe4ec5261e1b9) |
Where apps had overlapping implementations, the one that measured better was kept, with the
tests of both:
@@ -599,7 +633,8 @@ For the apps switching from their own code; to be removed once all three have sw
The codec modules' public items keep their names, so apps can re-export them
(`pub use craft_audio::aac::*;`). Names kept only for that are `#[doc(hidden)]`:
`ac3::Decoded` (now `ac3::Frame`), `aac::Upsample2x`, `aac::ProgramConfig::parse`,
`ac3::Decoded` (now `ac3::Frame`), `opus::encoder::OpusEncoder` (now `opus::Encoder`),
`aac::Upsample2x`, `aac::ProgramConfig::parse`,
`aac::config_layout` and `aac::sample_rate_index`. What else an app adjusts:
- `AudioInfo` has `container`, `codec`, `layout` and `damaged` instead of SoundCraft's `format`.
@@ -625,8 +660,10 @@ The codec modules' public items keep their names, so apps can re-export them
- `TrackInfo::sample_rate` and `channels` are `Option`s.
- Features are per codec. SoundCraft's `encode` (WAV, AIFF, FLAC) is now `lossless-encode`, made
of `pcm-encode` (WAV, AIFF) and `flac-encode`. `encode` turns on every encoder, including AAC,
ALAC and AC-3/E-AC-3. Writing a format whose feature is off returns `AudioError::Unsupported`.
ALAC, AC-3/E-AC-3 and Opus. Writing a format whose feature is off returns `AudioError::Unsupported`.
- AAC encoder settings it cannot use are `aac::Error::InvalidSettings` (`AudioError::Encode`).
- The Opus encoder's constructors and `encode_float` return `Result` instead of panicking on an
invalid channel count, mode/bandwidth or frame length.
- `TrackDecoder::decode` returns an error for a corrupt Opus packet, as for the other codecs;
call `conceal`.
- `FileEncoder::write` takes any per-channel slices, but refuses a block with missing channels or
+11 -4
View File
@@ -24,9 +24,9 @@
//! file bit for bit, and minimal files around such PCM, for essence embedded in other containers.
//! - [`peaks`] (feature `peaks`): min/max/RMS pyramid for waveform drawing.
//! - Codec modules for packet-level use, each with its own features: [`aac`] (`aac-decode`,
//! `aac-encode`), [`alac`] (`alac-decode`, `alac-encode`), [`opus`] (`opus-decode`), [`ac3`]
//! (`ac3-decode`), [`mpa`] (`mpa-decode`), [`mp3`] (`mp3-decode`, `mp3-encode`),
//! [`flac`] (`flac-encode`: the packet encoder).
//! `aac-encode`), [`alac`] (`alac-decode`, `alac-encode`), [`opus`] (`opus-decode`, and
//! [`opus::encoder`] with `opus-encode`), [`ac3`] (`ac3-decode`), [`mpa`] (`mpa-decode`), [`mp3`]
//! (`mp3-decode`, `mp3-encode`), [`flac`] (`flac-encode`: the packet encoder).
//!
//! The shared types below build with no features.
@@ -42,7 +42,7 @@ mod layout;
pub mod mp3;
#[cfg(feature = "mpa-decode")]
pub mod mpa;
#[cfg(feature = "opus-decode")]
#[cfg(any(feature = "opus-decode", feature = "opus-encode"))]
pub mod opus;
#[cfg(any(feature = "aac-decode", feature = "alac-decode", feature = "opus-decode", feature = "ac3-decode", feature = "mpa-decode"))]
mod pcm;
@@ -207,6 +207,13 @@ impl From<alac::Error> for AudioError {
}
}
#[cfg(feature = "opus-encode")]
impl From<opus::encoder::Error> for AudioError {
fn from(e: opus::encoder::Error) -> Self {
AudioError::Encode(e.to_string())
}
}
/// The file formats [`encode`] writes. Each needs its feature: WAV and AIFF `pcm-encode`, FLAC
/// `flac-encode`; asking for one that is not built returns [`AudioError::Unsupported`].
#[non_exhaustive]
+497
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::opus::encoder::range::{BITRES, RangeEncoder, ilog};
use crate::opus::encoder::rate::{Mode, QTHETA_OFFSET, get_pulses};
use crate::opus::encoder::tables::{EBANDS, NB_EBANDS};
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::opus::encoder::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}");
}
}
}
+230
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::opus::encoder::range::RangeEncoder;
use crate::opus::encoder::rate::MAX_FINE_BITS;
use crate::opus::encoder::tables::{BETA_COEF, BETA_INTRA, E_PROB_MODEL, NB_EBANDS, PRED_COEF, SMALL_ENERGY_ICDF};
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::opus::encoder::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}");
}
}
}
}
}
}
+294
View File
@@ -0,0 +1,294 @@
//! 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::opus::encoder::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");
}
}
+589
View File
@@ -0,0 +1,589 @@
//! Clean-room pure-Rust Opus encoder (RFC 6716 SILK, hybrid and CELT modes; RFC 7845 `OpusHead`).
//!
//! Produces 20 ms Opus packets from 48 kHz float input, mono or stereo. The coding mode and
//! audio bandwidth are chosen once per stream from the per-channel bitrate and the
//! [`Application`] (see the crate README for the table):
//!
//! - **SILK** (narrowband/mediumband/wideband): linear-predictive speech coding at an internal
//! rate of 8/12/16 kHz; packet sizes vary around the target bitrate.
//! - **Hybrid** (super-wideband/fullband): SILK wideband codes 0–8 kHz and CELT codes the bands
//! above in the same range-coded frame; constant packet size.
//! - **CELT** (fullband): MDCT transform coding; constant packet size.
//!
//! From EffectCraft's `effectcraft-opusenc` (storytold/effectcraft@c618429). Public items keep their
//! names, so EffectCraft can re-export this module.
//!
//! ```
//! # #[cfg(feature = "opus-encode")] {
//! use craft_audio::opus::encoder::{Application, Mode, Encoder};
//! let mut enc = Encoder::new(2, 128_000).unwrap();
//! let pcm = vec![0.0f32; Encoder::FRAME_SIZE * 2];
//! let packet = enc.encode_float(&pcm).unwrap();
//! assert_eq!(packet[0] >> 3, 31); // CELT-only, fullband, 20 ms
//! assert_eq!(enc.opus_head().len(), 19);
//! let voip = Encoder::with_application(1, 16_000, Application::Voip).unwrap();
//! assert_eq!(voip.mode(), Mode::Silk);
//! # }
//! ```
// 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 resample;
mod silk;
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 as CeltMode, compute_allocation};
use resample::Downsampler;
use silk::SilkEncoder;
use tables::{E_MEANS, EBANDS, NB_EBANDS, OVERLAP, PREEMPH, SPREAD_ICDF, TF_SELECT_TABLE, TRIM_ICDF};
/// 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 payload allowed in one frame (RFC 6716 §3.4, R2).
const MAX_FRAME_BYTES: usize = 1275;
/// First CELT band coded in hybrid frames (8 kHz).
const HYBRID_START_BAND: usize = 17;
fn celt_mode() -> &'static CeltMode {
static MODE: OnceLock<CeltMode> = OnceLock::new();
MODE.get_or_init(CeltMode::new)
}
/// Why an encoder could not be created or a frame could not be encoded.
#[non_exhaustive]
#[derive(Clone, Copy, Debug, PartialEq, Eq, thiserror::Error)]
pub enum Error {
/// Opus family-0 streams have 1 or 2 channels.
#[error("Opus family-0 streams have 1 or 2 channels, not {0}")]
Channels(u8),
/// The mode cannot code that bandwidth.
#[error("unsupported mode/bandwidth combination {0:?}/{1:?}")]
ModeBandwidth(Mode, Bandwidth),
/// `encode_float` takes exactly `FRAME_SIZE * channels` samples.
#[error("a frame is {expected} interleaved samples, not {got}")]
FrameLength {
/// Interleaved samples a frame takes (`FRAME_SIZE * channels`).
expected: usize,
/// Interleaved samples passed.
got: usize,
},
}
/// Intended use of the stream; selects the mode/bandwidth table.
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Hash)]
pub enum Application {
/// General audio (music, mixed content): CELT unless the bitrate is low.
#[default]
Audio,
/// Speech: SILK and hybrid over a wider bitrate range, narrower bandwidths at low rates.
Voip,
}
/// Opus coding mode (RFC 6716 §2).
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum Mode {
/// MDCT transform coding only.
Celt,
/// Linear-prediction coding only.
Silk,
/// SILK below 8 kHz plus CELT above.
Hybrid,
}
pub use crate::opus::Bandwidth;
/// Mode and bandwidth for a per-channel bitrate (bits per second per channel).
fn choose_mode(app: Application, per_channel: u32) -> (Mode, Bandwidth) {
match app {
Application::Voip => match per_channel {
0..10_000 => (Mode::Silk, Bandwidth::Narrowband),
10_000..12_000 => (Mode::Silk, Bandwidth::Mediumband),
12_000..18_000 => (Mode::Silk, Bandwidth::Wideband),
18_000..28_000 => (Mode::Hybrid, Bandwidth::SuperWideband),
28_000..40_000 => (Mode::Hybrid, Bandwidth::Fullband),
_ => (Mode::Celt, Bandwidth::Fullband),
},
Application::Audio => match per_channel {
0..16_000 => (Mode::Silk, Bandwidth::Wideband),
16_000..36_000 => (Mode::Hybrid, Bandwidth::Fullband),
_ => (Mode::Celt, Bandwidth::Fullband),
},
}
}
/// SILK side of a SILK-only or hybrid stream.
struct SilkPath {
enc: SilkEncoder,
down: Vec<Downsampler>,
/// Bits per frame aimed at (excluding the TOC byte in SILK-only mode).
target_bits: f64,
/// Bits saved (positive) or overspent by previous frames (SILK-only mode).
reservoir: f64,
}
/// The name EffectCraft's `effectcraft-opusenc` used, kept so it can re-export this module.
#[doc(hidden)]
pub type OpusEncoder = Encoder;
/// Opus encoder for one mono or stereo stream at 48 kHz.
pub struct Encoder {
channels: usize,
application: Application,
mode: Mode,
bandwidth: Bandwidth,
/// Requested bitrate (clamped), bits per second.
bitrate_bps: u32,
/// Payload bytes per packet in CELT and hybrid modes (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>,
silk: Option<SilkPath>,
}
impl Encoder {
/// Samples per channel per packet (20 ms at 48 kHz).
pub const FRAME_SIZE: usize = N;
/// Creates a general-audio encoder ([`Application::Audio`]) for `channels` (1 or 2)
/// channels of 48 kHz audio at `bitrate_bps` (total for all channels).
pub fn new(channels: u8, bitrate_bps: u32) -> Result<Self, Error> {
Self::with_application(channels, bitrate_bps, Application::Audio)
}
/// Creates an encoder for `channels` (1 or 2) channels of 48 kHz audio at `bitrate_bps`
/// (total for all channels; clamped to 6 kb/s per channel – 510 kb/s), choosing the mode
/// and bandwidth for `application`.
pub fn with_application(channels: u8, bitrate_bps: u32, application: Application) -> Result<Self, Error> {
let c = channels.max(1) as u32;
let (mode, bandwidth) = choose_mode(application, bitrate_bps.clamp(6_000 * c, 510_000) / c);
Self::build(channels, bitrate_bps, application, mode, bandwidth)
}
/// Creates an encoder with an explicit mode and bandwidth instead of the application's
/// table. Valid combinations: SILK with narrowband, mediumband or wideband; hybrid with
/// super-wideband or fullband; CELT with fullband.
///
/// Any other combination is [`Error::ModeBandwidth`]; `channels` other than 1 or 2 is
/// [`Error::Channels`].
pub fn with_mode(channels: u8, bitrate_bps: u32, mode: Mode, bandwidth: Bandwidth) -> Result<Self, Error> {
let ok = match mode {
Mode::Silk => bandwidth <= Bandwidth::Wideband,
Mode::Hybrid => bandwidth >= Bandwidth::SuperWideband,
Mode::Celt => bandwidth == Bandwidth::Fullband,
};
if !ok {
return Err(Error::ModeBandwidth(mode, bandwidth));
}
Self::build(channels, bitrate_bps, Application::Audio, mode, bandwidth)
}
fn build(channels: u8, bitrate_bps: u32, application: Application, mode: Mode, bandwidth: Bandwidth) -> Result<Self, Error> {
if channels != 1 && channels != 2 {
return Err(Error::Channels(channels));
}
let c = channels as usize;
let bitrate_bps = bitrate_bps.clamp(6_000 * c as u32, 510_000);
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);
let silk = match mode {
Mode::Celt => None,
Mode::Silk | Mode::Hybrid => {
let fs_khz = match bandwidth {
Bandwidth::Narrowband => 8,
Bandwidth::Mediumband => 12,
_ => 16,
};
let target_bits = if mode == Mode::Silk {
bitrate_bps as f64 / 50.0 - 8.0
} else {
// SILK gets most of a hybrid packet; CELT codes 8 kHz and up with the rest.
let per = (bitrate_bps / c as u32) as f64;
let share = if bandwidth == Bandwidth::SuperWideband { 0.72 } else { 0.62 };
c as f64 * (per * share).min(per - 6_000.0).max(per * 0.5) / 50.0
};
Some(SilkPath {
enc: SilkEncoder::new(fs_khz, c),
down: (0..c).map(|_| Downsampler::new(fs_khz as u32 * 1000)).collect(),
target_bits,
reservoir: 0.0,
})
}
};
Ok(Encoder {
channels: c,
application,
mode,
bandwidth,
bitrate_bps,
frame_bytes,
mdct: Mdct::new(N),
preemph_mem: [0.0; 2],
hist: vec![0.0; c * OVERLAP],
old_band_e: [0.0; 2 * NB_EBANDS],
block: vec![0.0; N + OVERLAP],
freq: vec![0.0; c * N],
silk,
})
}
/// Number of channels.
pub fn channels(&self) -> u8 {
self.channels as u8
}
/// The application the encoder was created for.
pub fn application(&self) -> Application {
self.application
}
/// Coding mode of every packet of the stream.
pub fn mode(&self) -> Mode {
self.mode
}
/// Coded audio bandwidth of every packet of the stream.
pub fn bandwidth(&self) -> Bandwidth {
self.bandwidth
}
/// Bitrate in bits per second: the constant packet size × 50 packets per second in CELT
/// and hybrid modes, the (average) target in SILK mode.
pub fn bitrate(&self) -> u32 {
match self.mode {
Mode::Silk => self.bitrate_bps,
_ => ((self.frame_bytes + 1) * 8 * 50) as u32,
}
}
/// Samples at 48 kHz the decoder must discard from the start of the stream (the `OpusHead`
/// pre-skip): the encoder's total delay. CELT: the 120-sample MDCT overlap. SILK/hybrid:
/// the decimator (83 samples, 88 for narrowband) plus the decoder's SILK delay (one
/// internal-rate sample and the resampler delay of RFC 6716 Table 54), which also rounds
/// to 120 samples.
pub fn pre_skip(&self) -> u16 {
let (fs, table54_ms) = match (self.mode, self.bandwidth) {
(Mode::Celt, _) => return OVERLAP as u16,
(_, Bandwidth::Narrowband) => (8_000, 0.538),
(_, Bandwidth::Mediumband) => (12_000, 0.692),
_ => (16_000, 0.706),
};
(resample::delay(fs) as f64 + 48_000.0 / fs as f64 + table54_ms * 48.0).round() 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
}
/// TOC byte (RFC 6716 §3.1, Table 2) of a 20 ms single-frame packet.
fn toc(&self) -> u8 {
let config: u8 = match (self.mode, self.bandwidth) {
(Mode::Silk, Bandwidth::Narrowband) => 1,
(Mode::Silk, Bandwidth::Mediumband) => 5,
(Mode::Silk, _) => 9,
(Mode::Hybrid, Bandwidth::SuperWideband) => 13,
(Mode::Hybrid, _) => 15,
(Mode::Celt, _) => 31,
};
(config << 3) | if self.channels == 2 { 4 } else { 0 }
}
/// Encodes one 20 ms frame of interleaved f32 samples in [-1, 1] (exactly
/// `FRAME_SIZE * channels` values) and returns one Opus packet, or
/// [`Error::FrameLength`] for any other length.
pub fn encode_float(&mut self, pcm: &[f32]) -> Result<Vec<u8>, Error> {
let c = self.channels;
if pcm.len() != N * c {
return Err(Error::FrameLength { expected: N * c, got: pcm.len() });
}
let mut packet = vec![self.toc()];
match self.mode {
Mode::Celt => {
let (log_e, silence) = self.celt_analysis(pcm);
let len = self.frame_bytes;
packet.extend_from_slice(&self.encode_celt(RangeEncoder::new(), len, 0, NB_EBANDS, &log_e, silence));
}
Mode::Silk => {
let x = self.silk_input(pcm);
// The SILK state exists whenever the mode uses SILK (set with the mode).
if let Some(sp) = self.silk.as_mut() {
let target = sp.target_bits + 0.25 * sp.reservoir;
let mut enc = RangeEncoder::new();
sp.enc.encode(&x, target.max(24.0), (MAX_FRAME_BYTES * 8) as i32, &mut enc);
let bytes = (enc.tell() as usize).div_ceil(8).max(2);
let lim = 4.0 * sp.target_bits;
sp.reservoir = (sp.reservoir + sp.target_bits - (bytes * 8) as f64).clamp(-lim, lim);
packet.extend_from_slice(&enc.done(bytes));
}
}
Mode::Hybrid => {
let x = self.silk_input(pcm);
let len = self.frame_bytes;
let total = (len * 8) as i32;
let celt_min = (total / 5).max(80);
let mut enc = RangeEncoder::new();
if let Some(sp) = self.silk.as_mut() {
sp.enc.encode(&x, sp.target_bits, total - celt_min, &mut enc);
}
// No redundancy (RFC 6716 §4.5.1); the flag is present when there is room.
if enc.tell() + 17 + 20 <= total {
enc.bit_logp(false, 12);
}
let (log_e, _) = self.celt_analysis(pcm);
let end = if self.bandwidth == Bandwidth::SuperWideband { 19 } else { NB_EBANDS };
packet.extend_from_slice(&self.encode_celt(enc, len, HYBRID_START_BAND, end, &log_e, false));
}
}
Ok(packet)
}
/// Decimates the input to the SILK internal rate (16-bit scale), one vector per channel.
fn silk_input(&mut self, pcm: &[f32]) -> Vec<Vec<f32>> {
let c = self.channels;
let Some(sp) = self.silk.as_mut() else { return vec![vec![]; c] };
let mut out = Vec::with_capacity(c);
let mut tmp = vec![0f32; N];
for ch in 0..c {
for (j, t) in tmp.iter_mut().enumerate() {
let s = pcm[j * c + ch];
*t = if s.is_finite() { s.clamp(-1.0, 1.0) * 32768.0 } else { 0.0 };
}
let mut y = Vec::with_capacity(N / 3);
sp.down[ch].process(&tmp, &mut y);
out.push(y);
}
out
}
/// Pre-emphasis, the MDCT, band energies (log2, relative to `E_MEANS`) and normalised
/// bands; returns the energies and whether the frame is digital silence.
fn celt_analysis(&mut self, pcm: &[f32]) -> ([f32; 2 * NB_EBANDS], bool) {
let c = self.channels;
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();
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;
}
}
}
(log_e, peak < 1e-4)
}
/// Codes one CELT frame (RFC 6716 §4.3, in bitstream order) of `len` bytes into `enc`
/// (which already holds the SILK frame in hybrid mode), coding bands `start..end`, and
/// finalises the frame.
#[allow(clippy::too_many_arguments)]
fn encode_celt(&mut self, mut enc: RangeEncoder, len: usize, start: usize, end: usize, log_e: &[f32; 2 * NB_EBANDS], silence: bool) -> Vec<u8> {
let m = celt_mode();
let c = self.channels;
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 && start == 0 {
// 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 total_bits_raw = len as i32 * 8;
if enc.tell() == 1 {
enc.bit_logp(false, 15);
}
// No pitch post-filter (only signalled in CELT-only frames).
if start == 0 && 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);
}
// The decoder forgets the energies of bands it did not code.
for ch in 0..2 {
for i in (0..start).chain(end..NB_EBANDS) {
self.old_band_e[ch * NB_EBANDS + i] = 0.0;
}
}
enc.done(len)
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn opus_head_layout() {
let e = Encoder::new(2, 96_000).unwrap();
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 = Encoder::new(ch, rate).unwrap();
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).unwrap();
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 = Encoder::new(1, 64_000).unwrap();
let p = e.encode_float(&vec![0.0; N]).unwrap();
assert_eq!(p.len(), 3);
}
#[test]
fn wrong_frame_length_is_an_error() {
let mut e = Encoder::new(2, 128_000).unwrap();
assert_eq!(e.encode_float(&vec![0.0; N]), Err(Error::FrameLength { expected: 2 * N, got: N }));
assert_eq!(e.encode_float(&[]), Err(Error::FrameLength { expected: 2 * N, got: 0 }));
assert!(e.encode_float(&vec![0.0; 2 * N]).is_ok());
}
}
+472
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
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::opus::encoder::range::{BITRES, RangeEncoder, ilog};
use crate::opus::encoder::tables::{ALLOC, EBANDS, LOG2_FRAC_TABLE, MAX_LM, NB_EBANDS};
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 {
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;
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);
}
}
+131
View File
@@ -0,0 +1,131 @@
//! 48 kHz → 8/12/16 kHz decimator feeding the SILK layer (own design).
//!
//! A linear-phase low-pass FIR (Kaiser-windowed sinc at 48 kHz) evaluated only at the retained
//! output phases. Its group delay ([`delay`]) is a whole number of 48 kHz samples, chosen so
//! that the SILK path (this filter + the decoder's internal-rate delay and resampler, RFC 6716
//! Table 54) lands within a fraction of a sample of the CELT path's 120-sample MDCT overlap.
/// Group delay of the decimator to `fs_out` in 48 kHz samples (the filter has `2 * delay + 1`
/// taps).
pub fn delay(fs_out: u32) -> usize {
if fs_out == 8_000 { 88 } else { 83 }
}
/// Zeroth-order modified Bessel function of the first kind (power series).
fn bessel_i0(x: f64) -> f64 {
let mut sum = 1.0;
let mut term = 1.0;
let q = x * x / 4.0;
for k in 1..64 {
term *= q / (k * k) as f64;
sum += term;
if term < sum * 1e-17 {
break;
}
}
sum
}
/// Windowed-sinc low-pass with cutoff `fc` (Hz at 48 kHz), unit DC gain.
fn design(fc: f64, beta: f64, delay: usize) -> Vec<f32> {
let wc = 2.0 * std::f64::consts::PI * fc / 48_000.0;
let i0b = bessel_i0(beta);
let mut h: Vec<f64> = (0..2 * delay + 1)
.map(|k| {
let m = k as f64 - delay as f64;
let sinc = if m == 0.0 { wc / std::f64::consts::PI } else { (wc * m).sin() / (std::f64::consts::PI * m) };
let r = m / delay as f64;
sinc * bessel_i0(beta * (1.0 - r * r).max(0.0).sqrt()) / i0b
})
.collect();
let s: f64 = h.iter().sum();
for v in h.iter_mut() {
*v /= s;
}
h.iter().map(|&v| v as f32).collect()
}
/// Decimates one channel of 48 kHz audio by 3, 4 or 6.
#[derive(Clone)]
pub struct Downsampler {
factor: usize,
h: Vec<f32>,
taps: usize,
/// The last `taps - 1` input samples followed by the current frame.
buf: Vec<f32>,
}
impl Downsampler {
/// `fs_out` is 8000, 12000 or 16000.
pub fn new(fs_out: u32) -> Self {
let factor = (48_000 / fs_out) as usize;
// Cutoff at 95 % of the output Nyquist frequency (pass band to ~88 %, stop band from ~110 %).
let fc = 0.95 * fs_out as f64 / 2.0;
let d = delay(fs_out);
let taps = 2 * d + 1;
Downsampler { factor, h: design(fc, 6.5, d), taps, buf: vec![0.0; taps - 1] }
}
/// Decimates `input` (a multiple of the factor long), appending to `out`. Output sample `j`
/// is the filtered input at frame position `factor * j`, delayed by [`delay`].
pub fn process(&mut self, input: &[f32], out: &mut Vec<f32>) {
debug_assert!(input.len().is_multiple_of(self.factor));
let taps = self.taps;
self.buf.truncate(taps - 1);
self.buf.extend_from_slice(input);
for j in 0..input.len() / self.factor {
let t = taps - 1 + self.factor * j;
let win = &self.buf[t + 1 - taps..=t];
// h is symmetric, so the window can be taken in forward order.
out.push(win.iter().zip(&self.h).map(|(a, b)| a * b).sum());
}
let n = self.buf.len();
self.buf.copy_within(n - (taps - 1).., 0);
self.buf.truncate(taps - 1);
}
}
#[cfg(test)]
mod tests {
use super::*;
fn tone(f: f64, n: usize) -> Vec<f32> {
(0..n).map(|i| (2.0 * std::f64::consts::PI * f * i as f64 / 48_000.0).sin() as f32).collect()
}
fn rms(x: &[f32]) -> f64 {
(x.iter().map(|&v| (v as f64).powi(2)).sum::<f64>() / x.len() as f64).sqrt()
}
#[test]
fn passes_band_and_rejects_aliases() {
for &(fs, pass, stop) in &[(16_000u32, 6_000.0, 8_600.0), (12_000, 4_500.0, 6_500.0), (8_000, 3_000.0, 4_500.0)] {
for &(f, want_pass) in &[(pass, true), (stop, false)] {
let mut d = Downsampler::new(fs);
let x = tone(f, 9600);
let mut y = Vec::new();
for c in x.chunks(960) {
d.process(c, &mut y);
}
let r = rms(&y[y.len() / 2..]) * std::f64::consts::SQRT_2;
if want_pass {
assert!((r - 1.0).abs() < 0.02, "{fs}: {f} Hz gain {r}");
} else {
assert!(r < 0.01, "{fs}: {f} Hz leaks {r}");
}
}
}
}
#[test]
fn delay_is_exact() {
let mut d = Downsampler::new(16_000);
let mut x = vec![0f32; 960];
x[300] = 1.0;
let mut y = Vec::new();
d.process(&x, &mut y);
// The impulse at 300 shows up centred at 48 kHz time 300 + 83 = 383 (output 127.67).
let peak = y.iter().enumerate().max_by(|a, b| a.1.total_cmp(b.1)).map(|(i, _)| i).unwrap_or(0);
assert!((peak as f64 - (300.0 + delay(16_000) as f64) / 3.0).abs() < 1.0, "peak at {peak}");
}
}
@@ -0,0 +1,322 @@
//! Encoder-side signal analysis for SILK (own design): windowed autocorrelation LPC, LPC → NLSF
//! conversion, open-loop pitch search and LTP codebook search.
use super::tables::{CB_NB_20, CB_WB_20, LTP_FILTER, LTP_PERIODICITY, LTP_TAPS_0, LTP_TAPS_1, LTP_TAPS_2};
/// Short-term prediction coefficients `a` (prediction `x[n] ≈ Σ a[j] x[n-1-j]`) from the
/// windowed segment `x`, with a Gaussian lag window and white-noise correction.
pub fn lpc(x: &[f32], order: usize, fs_khz: usize) -> Vec<f64> {
let n = x.len();
// Asymmetric window: sin² rise over the first half, flat, cos fall over the last 1/8.
let rise = n / 2;
let fall = n / 8;
let w = |i: usize| -> f64 {
if i < rise {
let s = (std::f64::consts::FRAC_PI_2 * (i as f64 + 0.5) / rise as f64).sin();
s * s
} else if i >= n - fall {
(std::f64::consts::FRAC_PI_2 * (i + fall - n) as f64 / fall as f64 + 0.0).cos()
} else {
1.0
}
};
let xw: Vec<f64> = x.iter().enumerate().map(|(i, &v)| v as f64 * w(i)).collect();
let mut r = vec![0f64; order + 1];
for (k, rk) in r.iter_mut().enumerate() {
*rk = xw[k..].iter().zip(&xw).map(|(a, b)| a * b).sum();
}
if r[0] <= 1e-6 {
return vec![0.0; order];
}
// White-noise correction (-40 dB) and a 50 Hz Gaussian lag window.
r[0] *= 1.0 + 1e-4;
let f0 = 2.0 * std::f64::consts::PI * 50.0 / (fs_khz as f64 * 1000.0);
for (k, rk) in r.iter_mut().enumerate().skip(1) {
*rk *= (-0.5 * (f0 * k as f64).powi(2)).exp();
}
let mut a = levinson(&r, order);
// Mild bandwidth expansion keeps the NLSFs apart.
let mut g = 1.0;
for v in a.iter_mut() {
g *= 0.996;
*v *= g;
}
a
}
/// Levinson-Durbin recursion on autocorrelation `r` (predictor convention as in [`lpc`]).
pub fn levinson(r: &[f64], order: usize) -> Vec<f64> {
let mut a = vec![0f64; order];
let mut err = r[0];
let mut tmp = vec![0f64; order];
for i in 0..order {
let mut acc = r[i + 1];
for j in 0..i {
acc -= a[j] * r[i - j];
}
let k = (acc / err).clamp(-0.9999, 0.9999);
tmp[..i].copy_from_slice(&a[..i]);
for j in 0..i {
a[j] = tmp[j] - k * tmp[i - 1 - j];
}
a[i] = k;
err *= 1.0 - k * k;
if err <= 0.0 {
break;
}
}
a
}
/// Evaluates the symmetric polynomial `p` (length `d + 1`, `d` even) as the cosine series
/// `e^{jωd/2} P(e^{jω})`.
fn cos_series(p: &[f64], w: f64) -> f64 {
let h = (p.len() - 1) / 2;
let mut s = p[h];
for k in 0..h {
s += 2.0 * p[k] * ((h - k) as f64 * w).cos();
}
s
}
fn roots(p: &[f64], out: &mut Vec<f64>) {
const GRID: usize = 1024;
let mut w0 = 0.0;
let mut f0 = cos_series(p, w0);
for g in 1..=GRID {
let w1 = std::f64::consts::PI * g as f64 / GRID as f64;
let f1 = cos_series(p, w1);
if f0 == 0.0 {
out.push(w0);
} else if f0 * f1 < 0.0 {
let (mut lo, mut hi, mut flo) = (w0, w1, f0);
for _ in 0..40 {
let mid = 0.5 * (lo + hi);
let fm = cos_series(p, mid);
if fm * flo <= 0.0 {
hi = mid;
} else {
lo = mid;
flo = fm;
}
}
out.push(0.5 * (lo + hi));
}
w0 = w1;
f0 = f1;
}
}
/// LPC → normalised line spectral frequencies in radians (ascending, `order` values), or `None`
/// when the root search fails (unstable or degenerate filter).
pub fn a2nlsf(a: &[f64]) -> Option<Vec<f64>> {
let d = a.len();
let mut c = vec![0f64; d + 2];
c[0] = 1.0;
for k in 0..d {
c[k + 1] = -a[k];
}
let mut p = vec![0f64; d + 2];
let mut q = vec![0f64; d + 2];
for k in 0..d + 2 {
p[k] = c[k] + c[d + 1 - k];
q[k] = c[k] - c[d + 1 - k];
}
// Remove the trivial roots at z = -1 (P) and z = +1 (Q).
let mut p1 = vec![0f64; d + 1];
let mut q1 = vec![0f64; d + 1];
p1[0] = p[0];
q1[0] = q[0];
for k in 1..=d {
p1[k] = p[k] - p1[k - 1];
q1[k] = q[k] + q1[k - 1];
}
let mut r = Vec::with_capacity(d);
roots(&p1, &mut r);
roots(&q1, &mut r);
if r.len() != d {
return None;
}
r.sort_by(f64::total_cmp);
Some(r)
}
/// Normalised cross-correlation of `x[n]` and `x[n - lag]` over `n` in `range`.
fn ncorr(x: &[f32], range: std::ops::Range<usize>, lag: usize) -> f64 {
let (mut xy, mut xx, mut yy) = (0f64, 0f64, 0f64);
for n in range {
let a = x[n] as f64;
let b = x[n - lag] as f64;
xy += a * b;
xx += a * a;
yy += b * b;
}
if xy <= 0.0 { 0.0 } else { xy / (xx * yy + 1e-9).sqrt() }
}
pub struct Pitch {
pub lag_index: i32,
pub contour: usize,
pub lags: [i32; 4],
pub score: f64,
}
/// Open-loop pitch search on `x` (history followed by the frame starting at `start`, `len`
/// samples, four subframes) at `fs_khz`.
pub fn pitch_search(x: &[f32], start: usize, len: usize, fs_khz: usize) -> Pitch {
let fs = fs_khz as i32;
let min_lag = 2 * fs;
let max_lag = 18 * fs - 1;
let sub = len / 4;
let nl = (max_lag - min_lag + 1) as usize;
// Whole-frame score per lag.
let full: Vec<f64> = (0..nl).map(|i| ncorr(x, start..start + len, min_lag as usize + i)).collect();
let mut best = 0usize;
for i in 0..nl {
// Slight preference for shorter lags.
if full[i] * (1.0 - 0.0004 * i as f64) > full[best] * (1.0 - 0.0004 * best as f64) {
best = i;
}
}
let mut lag = min_lag + best as i32;
// Prefer a sub-multiple when it explains the signal almost as well (octave errors).
for div in [4, 3, 2] {
let cand = (lag as f64 / div as f64).round() as i32;
let mut cbest = (0.0, 0);
for c in cand - 1..=cand + 1 {
if c >= min_lag && c <= max_lag {
let s = full[(c - min_lag) as usize];
if s > cbest.0 {
cbest = (s, c);
}
}
}
if cbest.1 != 0 && cbest.0 > 0.9 * full[(lag - min_lag) as usize] {
lag = cbest.1;
break;
}
}
let score = full[(lag - min_lag) as usize];
// Contour: per-subframe lag offsets around the chosen lag.
let cb: &[[i32; 4]] = if fs_khz == 8 { &CB_NB_20 } else { &CB_WB_20 };
let sub_score = |k: usize, l: i32| -> f64 {
let l = l.clamp(min_lag, 18 * fs);
ncorr(x, start + k * sub..start + (k + 1) * sub, l as usize)
};
let mut best = (f64::MIN, lag, 0usize);
for l0 in (lag - 2).max(min_lag)..=(lag + 2).min(max_lag) {
for (ci, offs) in cb.iter().enumerate() {
let s: f64 = (0..4).map(|k| sub_score(k, l0 + offs[k])).sum::<f64>() - 0.002 * ci as f64;
if s > best.0 {
best = (s, l0, ci);
}
}
}
let (_, l0, ci) = best;
let mut lags = [0i32; 4];
for k in 0..4 {
lags[k] = (l0 + cb[ci][k]).clamp(min_lag, 18 * fs);
}
Pitch { lag_index: l0 - min_lag, contour: ci, lags, score }
}
/// Result of the LTP codebook search.
pub struct Ltp {
pub per_index: usize,
pub ltp_index: [usize; 4],
/// Residual energy after LTP prediction / before, per subframe.
pub err: [f64; 4],
pub energy: [f64; 4],
}
fn icdf_bits(icdf: &[u8], s: usize) -> f64 {
let hi = if s == 0 { 256 } else { icdf[s - 1] as u32 };
let p = (hi - icdf[s] as u32).max(1) as f64 / 256.0;
-p.log2()
}
/// Searches the three LTP codebooks (RFC 6716 Tables 39–42) for the taps that best predict the
/// residual `r` (history + frame starting at `start`) from its past at the per-subframe `lags`.
pub fn ltp_search(r: &[f64], start: usize, sub: usize, lags: &[i32; 4]) -> Ltp {
let mut rr = [[0f64; 25]; 4];
let mut cv = [[0f64; 5]; 4];
let mut e = [0f64; 4];
for k in 0..4 {
let s0 = start + k * sub;
let lag = lags[k] as usize;
for n in s0..s0 + sub {
let t = r[n];
e[k] += t * t;
let base = n + 2 - lag;
for i in 0..5 {
let vi = r[base - i];
cv[k][i] += t * vi;
for j in 0..5 {
rr[k][i * 5 + j] += vi * r[base - j];
}
}
}
}
let books: [&[[i32; 5]]; 3] = [&LTP_TAPS_0, &LTP_TAPS_1, &LTP_TAPS_2];
let mut best = (f64::MAX, 0usize, [0usize; 4], [0f64; 4]);
for (p, book) in books.iter().enumerate() {
let mut total = icdf_bits(&LTP_PERIODICITY, p) * 1e-3 * e.iter().sum::<f64>() / 4.0;
let mut idx = [0usize; 4];
let mut errs = [0f64; 4];
for k in 0..4 {
let mut bk = (f64::MAX, 0usize, 0f64);
for (i, taps) in book.iter().enumerate() {
let b: [f64; 5] = std::array::from_fn(|j| taps[j] as f64 / 128.0);
let mut err = e[k];
for a in 0..5 {
err -= 2.0 * b[a] * cv[k][a];
for c in 0..5 {
err += b[a] * b[c] * rr[k][a * 5 + c];
}
}
let cost = err + icdf_bits(LTP_FILTER[p], i) * 1e-3 * e[k];
if cost < bk.0 {
bk = (cost, i, err.max(0.0));
}
}
total += bk.0;
idx[k] = bk.1;
errs[k] = bk.2;
}
if total < best.0 {
best = (total, p, idx, errs);
}
}
Ltp { per_index: best.1, ltp_index: best.2, err: best.3, energy: e }
}
#[cfg(test)]
mod tests {
use super::super::lpc::nlsf2a;
use super::*;
#[test]
fn nlsf_roundtrip_through_decoder_conversion() {
// A resonant 10th/16th-order predictor from a synthetic signal.
for &(order, fs) in &[(10usize, 8usize), (16, 16)] {
let mut st = [0f64; 4];
let mut seed = 1u32;
let x: Vec<f32> = (0..400)
.map(|_| {
seed = seed.wrapping_mul(1_664_525).wrapping_add(1_013_904_223);
let w = (seed >> 9) as f64 / (1u32 << 23) as f64 - 1.0;
let y = w + 1.6 * st[0] - 0.9 * st[1] + 0.3 * st[2] - 0.2 * st[3];
st = [y, st[0], st[1], st[2]];
(y * 1000.0) as f32
})
.collect();
let a = lpc(&x, order, fs);
let nlsf = a2nlsf(&a).expect("roots");
let q15: Vec<i32> = nlsf.iter().map(|w| (w / std::f64::consts::PI * 32768.0).round() as i32).collect();
let a_q12 = nlsf2a(&q15);
for j in 0..order {
assert!((a_q12[j] as f64 / 4096.0 - a[j]).abs() < 0.02, "order {order} coef {j}: {} vs {}", a_q12[j] as f64 / 4096.0, a[j]);
}
}
}
}
+370
View File
@@ -0,0 +1,370 @@
//! Fixed-point NLSF decoding, NLSF→LPC conversion and filter stabilisation
//! (RFC 6716 §4.2.7.5, with the RFC 8251 §6/§7 overflow fixes).
//!
//! This is the decoder's normative arithmetic, mirrored from FilmCraft's first-party,
//! spec-derived `filmcraft-opus` decoder (MIT OR Apache-2.0): the encoder must reconstruct the
//! quantised NLSFs and LPC coefficients bit-exactly as the decoder does.
use super::tables::{
COS_Q12, LSF_ORDER_NBMB, LSF_ORDER_WB, NLSF_CB_NBMB, NLSF_CB_WB, NLSF_MIN_SPACING_NBMB, NLSF_MIN_SPACING_WB, NLSF_PRED_NBMB, NLSF_PRED_WB,
NLSF_PSEL_NBMB, NLSF_PSEL_WB,
};
pub const MAX_LPC_ORDER: usize = 16;
#[inline]
pub fn smulwb(a: i32, b: i32) -> i32 {
((a as i64 * (b as i16) as i64) >> 16) as i32
}
#[inline]
pub fn smlawb(a: i32, b: i32, c: i32) -> i32 {
a.wrapping_add(smulwb(b, c))
}
#[inline]
pub fn smulww(a: i32, b: i32) -> i32 {
((a as i64 * b as i64) >> 16) as i32
}
#[inline]
pub fn smmul(a: i32, b: i32) -> i32 {
((a as i64 * b as i64) >> 32) as i32
}
#[inline]
pub fn rshift_round(a: i32, s: u32) -> i32 {
if s == 1 { (a >> 1) + (a & 1) } else { ((a >> (s - 1)) + 1) >> 1 }
}
#[inline]
pub fn rshift_round64(a: i64, s: u32) -> i64 {
if s == 1 { (a >> 1) + (a & 1) } else { ((a >> (s - 1)) + 1) >> 1 }
}
#[inline]
pub fn sat16(a: i32) -> i32 {
a.clamp(-32768, 32767)
}
/// `silk_INVERSE32_varQ`: approximation of `(1 << q_res) / b32`.
pub fn inverse32_varq(b32: i32, q_res: i32) -> i32 {
let b_headrm = b32.unsigned_abs().leading_zeros() as i32 - 1;
let b32_nrm = b32.wrapping_shl(b_headrm as u32);
let b32_inv = (i32::MAX >> 2) / (b32_nrm >> 16);
let mut result = b32_inv << 16;
let err_q32 = ((1i32 << 29).wrapping_sub(smulwb(b32_nrm, b32_inv))).wrapping_shl(3);
result = result.wrapping_add(((err_q32 as i64 * b32_inv as i64) >> 16) as i32);
let lshift = 61 - b_headrm - q_res;
if lshift <= 0 {
let s = (-lshift) as u32;
let lim_hi = i32::MAX >> s;
let lim_lo = i32::MIN >> s;
result.clamp(lim_lo, lim_hi) << s
} else if lshift < 32 {
result >> lshift
} else {
0
}
}
/// `silk_DIV32_varQ`: approximation of `(a32 << q_res) / b32`.
pub fn div32_varq(a32: i32, b32: i32, q_res: i32) -> i32 {
let a_headrm = a32.unsigned_abs().leading_zeros() as i32 - 1;
let a32_nrm = a32.wrapping_shl(a_headrm as u32);
let b_headrm = b32.unsigned_abs().leading_zeros() as i32 - 1;
let b32_nrm = b32.wrapping_shl(b_headrm as u32);
let b32_inv = (i32::MAX >> 2) / (b32_nrm >> 16);
let mut result = smulwb(a32_nrm, b32_inv);
let a32_nrm = a32_nrm.wrapping_sub(((smmul(b32_nrm, result) as i64) << 3) as i32);
result = smlawb(result, a32_nrm, b32_inv);
let lshift = 29 + a_headrm - b_headrm - q_res;
if lshift < 0 {
let s = (-lshift) as u32;
let lim_hi = i32::MAX >> s;
let lim_lo = i32::MIN >> s;
result.clamp(lim_lo, lim_hi) << s
} else if lshift < 32 {
result >> lshift
} else {
0
}
}
pub fn sqrt_approx(x: i32) -> i32 {
if x <= 0 {
return 0;
}
let lz = x.leading_zeros() as i32;
let frac = (x.rotate_right((24 - lz).rem_euclid(32) as u32)) & 0x7f;
let mut y = if lz & 1 != 0 { 32768 } else { 46214 };
y >>= lz >> 1;
smlawb(y, y, 213 * frac)
}
/// Decodes the stage-1/stage-2 indices into stabilised NLSFs (Q15).
pub fn nlsf_decode(indices: &[i32], wb: bool) -> [i32; MAX_LPC_ORDER] {
let order = if wb { 16 } else { 10 };
let i1 = indices[0] as usize;
let cb1: &[i32] = if wb { &NLSF_CB_WB[i1] } else { &NLSF_CB_NBMB[i1] };
let qstep = if wb { 9830 } else { 11796 };
// Residual dequantisation with backwards prediction.
let mut res_q10 = [0i32; MAX_LPC_ORDER];
let mut out_q10 = 0i32;
for i in (0..order).rev() {
let pred = if i + 1 < order {
let sel = if wb { NLSF_PSEL_WB[i1][i] } else { NLSF_PSEL_NBMB[i1][i] } as usize;
let w = if wb { NLSF_PRED_WB[sel][i] } else { NLSF_PRED_NBMB[sel][i] };
(out_q10 * w) >> 8
} else {
0
};
let mut v = indices[i + 1] << 10;
if v > 0 {
v -= 102;
} else if v < 0 {
v += 102;
}
out_q10 = smlawb(pred, v, qstep);
res_q10[i] = out_q10;
}
// Laroia weights.
let mut nlsf = [0i32; MAX_LPC_ORDER];
for k in 0..order {
let prev = if k == 0 { 0 } else { cb1[k - 1] };
let next = if k + 1 == order { 256 } else { cb1[k + 1] };
let w = (1024 / (cb1[k] - prev).max(1) + 1024 / (next - cb1[k]).max(1)).min(32767);
let w_q9 = sqrt_approx(w << 16);
let v = (cb1[k] << 7) + (res_q10[k] << 14) / w_q9.max(1);
nlsf[k] = v.clamp(0, 32767);
}
let dmin: &[i32] = if wb { &NLSF_MIN_SPACING_WB } else { &NLSF_MIN_SPACING_NBMB };
nlsf_stabilize(&mut nlsf[..order], dmin);
nlsf
}
pub fn nlsf_stabilize(nlsf: &mut [i32], dmin: &[i32]) {
let l = nlsf.len();
for _ in 0..20 {
let mut min_diff = nlsf[0] - dmin[0];
let mut idx = 0;
for i in 1..l {
let d = nlsf[i] - (nlsf[i - 1] + dmin[i]);
if d < min_diff {
min_diff = d;
idx = i;
}
}
let d = 32768 - (nlsf[l - 1] + dmin[l]);
if d < min_diff {
min_diff = d;
idx = l;
}
if min_diff >= 0 {
return;
}
if idx == 0 {
nlsf[0] = dmin[0];
} else if idx == l {
nlsf[l - 1] = 32768 - dmin[l];
} else {
let mut min_center = dmin[..idx].iter().sum::<i32>();
min_center += dmin[idx] >> 1;
let mut max_center = 32768 - dmin[idx + 1..=l].iter().sum::<i32>();
max_center -= dmin[idx] >> 1;
let c = rshift_round(nlsf[idx - 1] + nlsf[idx], 1);
let c = if min_center > max_center { c.clamp(max_center, min_center) } else { c.clamp(min_center, max_center) };
nlsf[idx - 1] = c - (dmin[idx] >> 1);
nlsf[idx] = nlsf[idx - 1] + dmin[idx];
}
}
// Fallback.
nlsf.sort_unstable();
nlsf[0] = nlsf[0].max(dmin[0]);
for i in 1..l {
nlsf[i] = nlsf[i].max(sat16(nlsf[i - 1] + dmin[i]));
}
nlsf[l - 1] = nlsf[l - 1].min(32768 - dmin[l]);
for i in (0..l - 1).rev() {
nlsf[i] = nlsf[i].min(nlsf[i + 1] - dmin[i + 1]);
}
}
fn bwexpander_32(ar: &mut [i32], mut chirp_q16: i32) {
let d = ar.len();
let chirp_minus_one = chirp_q16 - 65536;
for v in ar.iter_mut().take(d - 1) {
*v = smulww(chirp_q16, *v);
chirp_q16 += rshift_round(chirp_q16.wrapping_mul(chirp_minus_one), 16);
}
ar[d - 1] = smulww(chirp_q16, ar[d - 1]);
}
fn lpc_inverse_pred_gain(a_q12: &[i32]) -> i32 {
const QA: u32 = 24;
const A_LIMIT: i32 = 16_773_022;
let order = a_q12.len();
let mut a = [0i32; MAX_LPC_ORDER];
let mut dc = 0;
for k in 0..order {
dc += a_q12[k];
a[k] = a_q12[k] << (QA - 12);
}
if dc >= 4096 {
return 0;
}
let mut inv_gain: i32 = 1 << 30;
let mul32_frac_q31 = |x: i32, y: i32| rshift_round64(x as i64 * y as i64, 31) as i32;
for k in (1..order).rev() {
if a[k] > A_LIMIT || a[k] < -A_LIMIT {
return 0;
}
let rc_q31 = -(a[k] << (31 - QA));
let rc_mult1 = (1i32 << 30) - smmul(rc_q31, rc_q31);
inv_gain = smmul(inv_gain, rc_mult1) << 2;
if inv_gain < 107_374 {
return 0;
}
let mult2q = 32 - rc_mult1.unsigned_abs().leading_zeros() as i32;
let rc_mult2 = inverse32_varq(rc_mult1, mult2q + 30);
for n in 0..(k + 1) >> 1 {
let t1 = a[n];
let t2 = a[k - n - 1];
let v1 = rshift_round64(t1.saturating_sub(mul32_frac_q31(t2, rc_q31)) as i64 * rc_mult2 as i64, mult2q as u32);
if v1 > i32::MAX as i64 || v1 < i32::MIN as i64 {
return 0;
}
let v2 = rshift_round64(t2.saturating_sub(mul32_frac_q31(t1, rc_q31)) as i64 * rc_mult2 as i64, mult2q as u32);
if v2 > i32::MAX as i64 || v2 < i32::MIN as i64 {
return 0;
}
a[n] = v1 as i32;
a[k - n - 1] = v2 as i32;
}
}
if a[0] > A_LIMIT || a[0] < -A_LIMIT {
return 0;
}
let rc_q31 = -(a[0] << (31 - QA));
let rc_mult1 = (1i32 << 30) - smmul(rc_q31, rc_q31);
inv_gain = smmul(inv_gain, rc_mult1) << 2;
if inv_gain < 107_374 {
return 0;
}
inv_gain
}
/// Converts NLSFs (Q15) to stable Q12 LPC coefficients (`silk_NLSF2A`).
pub fn nlsf2a(nlsf: &[i32]) -> [i32; MAX_LPC_ORDER] {
let d = nlsf.len();
let ordering: &[usize] = if d == 16 { &LSF_ORDER_WB } else { &LSF_ORDER_NBMB };
let mut c_q17 = [0i32; MAX_LPC_ORDER];
for k in 0..d {
let i = (nlsf[k] >> 8) as usize;
let f = nlsf[k] & 255;
let c0 = COS_Q12[i];
let c1 = COS_Q12[(i + 1).min(128)];
c_q17[ordering[k]] = (c0 * 256 + (c1 - c0) * f + 4) >> 3;
}
let d2 = d / 2;
let find_poly = |off: usize| -> [i64; 10] {
let mut out = [0i64; 10];
out[0] = 1 << 16;
out[1] = -(c_q17[off] as i64);
for k in 1..d2 {
let ftmp = c_q17[2 * k + off] as i64;
out[k + 1] = 2 * out[k - 1] - ((ftmp * out[k] + 32768) >> 16);
for n in (2..=k).rev() {
out[n] += out[n - 2] - ((ftmp * out[n - 1] + 32768) >> 16);
}
out[1] -= ftmp;
}
out
};
let p = find_poly(0);
let q = find_poly(1);
let mut a32 = [0i32; MAX_LPC_ORDER];
for k in 0..d2 {
let ptmp = p[k + 1] + p[k];
let qtmp = q[k + 1] - q[k];
a32[k] = (-qtmp - ptmp) as i32;
a32[d - k - 1] = (qtmp - ptmp) as i32;
}
let a32 = &mut a32[..d];
// Limit the dynamic range (silk_LPC_fit with QIN = 17, QOUT = 12).
let mut a_q12 = [0i32; MAX_LPC_ORDER];
let mut i = 0;
while i < 10 {
let mut maxabs = 0i32;
let mut idx = 0;
for (k, &v) in a32.iter().enumerate() {
let av = v.saturating_abs();
if av > maxabs {
maxabs = av;
idx = k;
}
}
let maxabs = rshift_round(maxabs, 5);
if maxabs > 32767 {
let maxabs = maxabs.min(163_838);
let chirp = 65470 - ((maxabs - 32767) << 14) / ((maxabs * (idx as i32 + 1)) >> 2);
bwexpander_32(a32, chirp);
} else {
break;
}
i += 1;
}
if i == 10 {
for k in 0..d {
a_q12[k] = sat16(rshift_round(a32[k], 5));
a32[k] = a_q12[k] << 5;
}
} else {
for k in 0..d {
a_q12[k] = rshift_round(a32[k], 5);
}
}
for i in 0..16 {
if lpc_inverse_pred_gain(&a_q12[..d]) == 0 {
bwexpander_32(a32, 65536 - (2 << i));
for k in 0..d {
a_q12[k] = rshift_round(a32[k], 5);
}
} else {
break;
}
}
a_q12
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn flat_nlsf_gives_small_lpc() {
// Uniformly spaced NLSFs correspond to an (almost) flat spectrum.
let nlsf: Vec<i32> = (1..=16).map(|k| k * 32768 / 17).collect();
let a = nlsf2a(&nlsf);
assert!(a.iter().all(|v| v.abs() < 600), "{a:?}");
let nlsf: Vec<i32> = (1..=10).map(|k| k * 32768 / 11).collect();
let a = nlsf2a(&nlsf);
assert!(a[..10].iter().all(|v| v.abs() < 600), "{a:?}");
}
#[test]
fn stabilize_enforces_spacing() {
let mut v = [100, 90, 80, 30000, 30001, 30002, 31000, 32000, 32700, 32760];
nlsf_stabilize(&mut v, &NLSF_MIN_SPACING_NBMB);
assert!(v[0] >= NLSF_MIN_SPACING_NBMB[0]);
for i in 1..10 {
assert!(v[i] - v[i - 1] >= NLSF_MIN_SPACING_NBMB[i], "{v:?}");
}
assert!(32768 - v[9] >= NLSF_MIN_SPACING_NBMB[10]);
}
#[test]
fn inverse_approximations() {
for &b in &[65536, 100_000, 1 << 20, 81920, 1_686_110_208] {
let inv = inverse32_varq(b, 47) as f64;
let exact = 2f64.powi(47) / b as f64;
assert!((inv - exact).abs() / exact < 1e-3, "{b}");
}
let q = div32_varq(65536, 131072, 16) as f64;
assert!((q - 32768.0).abs() < 2.0);
}
}
+714
View File
@@ -0,0 +1,714 @@
//! SILK layer encoder (RFC 6716 §4.2, encoder direction), 20 ms frames at 8, 12 or 16 kHz.
//!
//! Each frame is analysed (VAD, open-loop pitch, windowed-autocorrelation LPC, LTP codebook
//! search), its NLSFs are quantised, and the excitation is chosen sample by sample by running
//! the decoder's exact fixed-point reconstruction (LTP and LPC synthesis with the dequantised
//! Q12/Q14 coefficients, quantisation offsets and the pseudo-random sign sequence), so the
//! encoder always holds the decoder's state and there is no drift. A noise-feedback loop
//! shapes the coding noise with a bandwidth-expanded LPC envelope, and a per-frame search over
//! a gain offset hits the bit budget. Stereo input is coded as mid/side with the decoder's
//! stereo prediction (§4.2.8). No LBRR, no NLSF interpolation.
pub(crate) mod analysis;
pub(crate) mod lpc;
mod nlsf;
mod pulses;
#[allow(dead_code)]
pub(crate) mod tables;
use lpc::{MAX_LPC_ORDER, div32_varq, inverse32_varq, nlsf2a, rshift_round, sat16, smlawb, smulwb, smulww};
use tables::{
DELTA_GAIN, GAIN_MSB, LTP_FILTER, LTP_PERIODICITY, LTP_SCALE, LTP_TAPS_0, LTP_TAPS_1, LTP_TAPS_2, MID_ONLY, NLSF_EXT, NLSF_INTERP, NLSF_SEL_NBMB,
NLSF_SEL_WB, NLSF_STAGE1, NLSF_STAGE2, PITCH_CONTOUR, PITCH_HIGH, PITCH_LOW, STEREO_PRED_JOINT, STEREO_WEIGHTS_Q13, TYPE_ACTIVE, TYPE_INACTIVE,
UNIFORM3, UNIFORM4, UNIFORM5, UNIFORM8,
};
use crate::opus::encoder::range::RangeEncoder;
const MAX_FRAME: usize = 320;
const LTP_ORDER: usize = 5;
/// Frame RMS (16-bit scale) below which a frame is coded as inactive (VAD off).
const VAD_RMS: f32 = 30.0;
/// Normalised pitch correlation above which a frame is a voicing candidate.
const VOICING_THRESHOLD: f64 = 0.45;
/// Minimum LTP prediction gain (energy ratio) for a frame to be coded as voiced.
const MIN_LTP_GAIN: f64 = 1.3;
/// Rate penalty (squared pulse units per unit magnitude) of the excitation quantiser.
const RD_LAMBDA: f64 = 0.15;
/// Extra gain offset (log2) for inactive frames.
const INACTIVE_OFFSET: f64 = 1.0;
/// Noise feedback: the coding noise is shaped by `1 / A(z/γ)` (A the short-term predictor).
/// γ = 0 leaves it white (closed-loop DPCM, best waveform SNR at high rates), γ = 1 shapes it
/// like the spectral envelope (open-loop residual coding, needed at low rates). γ grows with the
/// gain offset, i.e. as the quantisation gets coarser.
fn nf_gamma(offset: f64) -> f64 {
(0.75 + 0.05 * offset).clamp(0.75, 0.9)
}
/// `silk_log2lin` (decoder arithmetic).
fn log2lin(in_log_q7: i32) -> i32 {
if in_log_q7 < 0 {
return 0;
}
if in_log_q7 >= 3967 {
return i32::MAX;
}
let out = 1i32 << (in_log_q7 >> 7);
let frac = in_log_q7 & 0x7F;
if in_log_q7 < 2048 {
out + ((out * smlawb(frac, frac * (128 - frac), -174)) >> 7)
} else {
out + (out >> 7) * smlawb(frac, frac * (128 - frac), -174)
}
}
fn gain_q16(index: i32) -> i32 {
log2lin((smulwb(0x1D1C71, index) + 2090).min(3967))
}
/// The decoder's LPC analysis filter on 16-bit samples.
fn lpc_analysis_filter(out: &mut [i16], input: &[i16], b: &[i32], d: usize) {
let len = out.len();
for ix in d..len {
let mut acc: i32 = 0;
for j in 0..d {
acc = acc.wrapping_add((input[ix - 1 - j] as i32).wrapping_mul(b[j]));
}
let v = ((input[ix] as i32) << 12).wrapping_sub(acc);
out[ix] = sat16(rshift_round(v, 12)) as i16;
}
for v in out.iter_mut().take(d.min(len)) {
*v = 0;
}
}
/// The part of the decoder's channel state that the reconstruction depends on.
#[derive(Clone)]
struct DecState {
prev_gain_q16: i32,
s_lpc_q14: [i32; MAX_LPC_ORDER],
out_buf: [i16; MAX_FRAME],
last_gain_index: i32,
}
impl DecState {
fn new() -> Self {
DecState { prev_gain_q16: 65536, s_lpc_q14: [0; MAX_LPC_ORDER], out_buf: [0; MAX_FRAME], last_gain_index: 10 }
}
/// The decoder's reset of the side channel when it is coded again after mid-only frames.
fn reset_side(&mut self) {
self.out_buf = [0; MAX_FRAME];
self.s_lpc_q14 = [0; MAX_LPC_ORDER];
self.last_gain_index = 10;
}
}
/// Coded parameters of one frame (everything except gains and pulses).
#[derive(Clone, Default)]
struct Params {
vad: bool,
signal_type: usize,
quant_offset_type: usize,
nlsf_ix: [i32; 17],
a_q12: [i32; MAX_LPC_ORDER],
lag_index: i32,
contour: usize,
per_index: usize,
ltp_index: [usize; 4],
pitch_l: [i32; 4],
/// RMS (16-bit scale) of the open-loop prediction residual per subframe.
res_rms: [f32; 4],
}
/// Gains, seed and pulses chosen by the quantiser.
struct Coded {
gain_ix: [i32; 4],
gains_q16: [i32; 4],
seed: usize,
pulses: Vec<i32>,
}
#[derive(Clone)]
struct Channel {
fs_khz: usize,
len: usize,
sub: usize,
order: usize,
dec: DecState,
/// The last two frames of input (internal rate, 16-bit scale).
hist: Vec<f32>,
/// Recent output errors (decoded − input), most recent first, for noise feedback.
nf: [f64; MAX_LPC_ORDER],
}
impl Channel {
fn new(fs_khz: usize) -> Self {
let len = 20 * fs_khz;
Channel {
fs_khz,
len,
sub: 5 * fs_khz,
order: if fs_khz == 16 { 16 } else { 10 },
dec: DecState::new(),
hist: vec![0.0; 2 * len],
nf: [0.0; MAX_LPC_ORDER],
}
}
fn push_history(&mut self, x: &[f32]) {
let l = self.len;
self.hist.copy_within(l.., 0);
self.hist[l..].copy_from_slice(x);
}
/// Frame analysis: VAD, signal type, pitch, LPC/NLSF quantisation and LTP taps.
fn analyze(&self, x: &[f32]) -> Params {
let l = self.len;
let sub = self.sub;
let order = self.order;
let wb = order == 16;
let mut buf = self.hist.clone();
buf.extend_from_slice(x);
let start = 2 * l;
let rms = (x.iter().map(|v| v * v).sum::<f32>() / l as f32).sqrt();
let vad = rms > VAD_RMS;
let pitch = if vad { Some(analysis::pitch_search(&buf, start, l, self.fs_khz)) } else { None };
let mut voiced = pitch.as_ref().is_some_and(|p| p.score > VOICING_THRESHOLD);
let a = analysis::lpc(&buf[start - l / 4..], order, self.fs_khz);
let target: Vec<f64> = match analysis::a2nlsf(&a) {
Some(w) => w.iter().map(|v| v / std::f64::consts::PI * 32768.0).collect(),
None => (1..=order).map(|k| k as f64 * 32768.0 / (order + 1) as f64).collect(),
};
let quant = |signal_type: usize| {
let (ix, nlsf) = nlsf::quantize(&target, signal_type, wb);
(ix, nlsf2a(&nlsf[..order]))
};
let mut p = Params { vad, ..Default::default() };
let (mut ix, mut a_q12) = quant(if voiced { 2 } else { 1 });
// Open-loop residual with the quantised predictor.
let mut r = vec![0f64; buf.len()];
for n in order..buf.len() {
let mut acc = buf[n] as f64;
for j in 0..order {
acc -= a_q12[j] as f64 / 4096.0 * buf[n - 1 - j] as f64;
}
r[n] = acc;
}
if voiced && let Some(pitch) = &pitch {
let ltp = analysis::ltp_search(&r, start, sub, &pitch.lags);
let gain = ltp.energy.iter().sum::<f64>() / ltp.err.iter().sum::<f64>().max(1e-9);
if gain >= MIN_LTP_GAIN {
p.lag_index = pitch.lag_index;
p.contour = pitch.contour;
p.pitch_l = pitch.lags;
p.per_index = ltp.per_index;
p.ltp_index = ltp.ltp_index;
} else {
voiced = false;
(ix, a_q12) = quant(1);
}
}
// Gains follow the short-term residual level (the long-term predictor then only
// reduces the number of pulses), which keeps the noise-to-signal ratio steady.
for k in 0..4 {
let e: f64 = r[start + k * sub..start + (k + 1) * sub].iter().map(|v| v * v).sum();
p.res_rms[k] = (e / sub as f64).sqrt() as f32;
}
p.signal_type = if !vad {
0
} else if voiced {
2
} else {
1
};
p.quant_offset_type = 0;
p.nlsf_ix = ix;
p.a_q12 = a_q12;
p
}
/// Chooses gains and pulses for target `x` by running the decoder's reconstruction, and
/// advances the mirrored decoder state. `offset` scales the gains (log2), `floor` gives
/// minimum gains (Q16) per subframe. Returns the coded values and the decoded frame.
fn quantize(&mut self, p: &Params, x: &[f32], offset: f64, floor: Option<&[i32; 4]>, seed: usize) -> (Coded, Vec<i16>) {
let l = self.len;
let sub = self.sub;
let order = self.order;
let st = &mut self.dec;
let nf = &mut self.nf;
let off = if p.signal_type == 0 { offset + INACTIVE_OFFSET } else { offset };
let gamma = nf_gamma(off);
let mut nf_c = [0f64; MAX_LPC_ORDER];
let mut gpow = 1.0;
for j in 0..order {
gpow *= gamma;
nf_c[j] = p.a_q12[j] as f64 / 4096.0 * gpow;
}
// Gains (independent coding of the first subframe, deltas for the rest).
let mut gain_ix = [0i32; 4];
let mut gains_q16 = [0i32; 4];
let mut prev = st.last_gain_index;
for k in 0..4 {
let mut g = p.res_rms[k] as f64 * off.exp2();
if let Some(f) = floor {
g = g.max(f[k] as f64 / 65536.0);
}
let want = (((g * 65536.0).max(1.0).log2() * 128.0 - 2090.0) / (1_907_825.0 / 65536.0)).round().clamp(0.0, 63.0) as i32;
if k == 0 {
gain_ix[0] = want;
prev = want.max(prev - 16).clamp(0, 63);
} else {
let mut best = (i32::MAX, 0i32, prev);
for ind in 0..41 {
let tmp = ind - 4;
let thr = 2 * 36 - 64 + prev;
let d = if tmp > thr { prev + (tmp << 1) - thr } else { prev + tmp }.clamp(0, 63);
let e = (d - want).abs();
if e < best.0 {
best = (e, ind, d);
}
}
gain_ix[k] = best.1;
prev = best.2;
}
gains_q16[k] = gain_q16(prev);
}
st.last_gain_index = prev;
// Analysis-by-synthesis excitation quantisation through the decoder's exact synthesis.
let voiced = p.signal_type == 2;
let offset_q10 = [[100, 240], [32, 100]][p.signal_type >> 1][p.quant_offset_type];
let off16 = offset_q10 << 4;
let mut seed_v = seed as i32;
let mut pulses = vec![0i32; l];
let mut xq = vec![0i16; l];
let mut s_lpc = [0i32; MAX_FRAME / 4 + MAX_LPC_ORDER];
s_lpc[..MAX_LPC_ORDER].copy_from_slice(&st.s_lpc_q14);
let mut s_ltp = [0i16; MAX_FRAME];
let mut s_ltp_q15 = [0i32; 2 * MAX_FRAME];
let ltp_mem = l;
let mut ltp_buf_idx = ltp_mem;
let a_q12 = &p.a_q12;
let ltp_scale_q14 = 15565;
for k in 0..4 {
let mut b_q14 = [0i32; LTP_ORDER];
if voiced {
let ix = p.ltp_index[k];
let taps = match p.per_index {
0 => LTP_TAPS_0[ix],
1 => LTP_TAPS_1[ix],
_ => LTP_TAPS_2[ix],
};
for i in 0..LTP_ORDER {
b_q14[i] = taps[i] << 7;
}
}
let gain_q10 = gains_q16[k] >> 6;
let mut inv_gain_q31 = inverse32_varq(gains_q16[k], 47);
let gain_adj_q16 = if gains_q16[k] != st.prev_gain_q16 {
let g = div32_varq(st.prev_gain_q16, gains_q16[k], 16);
for v in s_lpc[..MAX_LPC_ORDER].iter_mut() {
*v = smulww(g, *v);
}
g
} else {
1 << 16
};
st.prev_gain_q16 = gains_q16[k];
let lag = p.pitch_l[k];
if voiced {
let lag_u = lag as usize;
if k == 0 {
let start_idx = (ltp_mem as isize - lag as isize - order as isize - (LTP_ORDER / 2) as isize).max(0) as usize;
lpc_analysis_filter(&mut s_ltp[start_idx..ltp_mem], &st.out_buf[start_idx..ltp_mem], a_q12, order);
inv_gain_q31 = smulwb(inv_gain_q31, ltp_scale_q14) << 2;
for i in 0..(lag_u + LTP_ORDER / 2).min(ltp_mem) {
s_ltp_q15[ltp_buf_idx - i - 1] = smulwb(inv_gain_q31, s_ltp[ltp_mem - i - 1] as i32);
}
} else if gain_adj_q16 != 1 << 16 {
for i in 0..(lag_u + LTP_ORDER / 2).min(ltp_buf_idx) {
let q = ltp_buf_idx - i - 1;
s_ltp_q15[q] = smulww(gain_adj_q16, s_ltp_q15[q]);
}
}
}
let inv_g = 16_777_216.0 / gain_q10.max(1) as f64;
for i in 0..sub {
let n = k * sub + i;
seed_v = seed_v.wrapping_mul(196_314_165).wrapping_add(907_633_515);
let ltp_q14 = if voiced {
let base = ltp_buf_idx as isize - lag as isize + (LTP_ORDER / 2) as isize;
let g = |o: isize| -> i32 { if base - o >= 0 { s_ltp_q15[(base - o) as usize] } else { 0 } };
let mut pred = 2i32;
for o in 0..LTP_ORDER {
pred = smlawb(pred, g(o as isize), b_q14[o]);
}
pred << 1
} else {
0
};
let mut lpc_pred = (order >> 1) as i32;
for j in 0..order {
lpc_pred = smlawb(lpc_pred, s_lpc[MAX_LPC_ORDER + i - 1 - j], a_q12[j]);
}
let lpc_q14 = lpc_pred.saturating_mul(16);
// Target excitation (Q14) for this sample.
let fb: f64 = (0..order).map(|j| nf_c[j] * nf[j]).sum();
let tgt = (x[n] as f64 + fb).clamp(-32768.0, 32767.0) * inv_g - lpc_q14 as f64 - ltp_q14 as f64;
let sgn = if seed_v < 0 { -1.0 } else { 1.0 };
let u = sgn * tgt;
let pf = (u - off16 as f64) / 16384.0;
let p0 = if pf > 0.0 { (pf + 0.078_125).round() } else { (pf - 0.078_125).round() } as i32;
let level = |q: i32| -> i32 {
let mut e = q << 14;
if e > 0 {
e -= 80 << 4;
} else if e < 0 {
e += 80 << 4;
}
e + off16
};
let mut best = (f64::MAX, 0i32);
for q in [p0 - 1, p0, p0 + 1, 0] {
let q = q.clamp(-pulses::MAX_PULSE, pulses::MAX_PULSE);
let d = (level(q) as f64 - u) / 16384.0;
let cost = d * d + RD_LAMBDA * q.abs() as f64;
if cost < best.0 {
best = (cost, q);
}
}
let q = best.1;
pulses[n] = q;
let mut e = level(q);
if seed_v < 0 {
e = -e;
}
seed_v = seed_v.wrapping_add(q);
let res = e.wrapping_add(ltp_q14);
if voiced {
s_ltp_q15[ltp_buf_idx] = res << 1;
ltp_buf_idx += 1;
}
let v = res.saturating_add(lpc_q14);
s_lpc[MAX_LPC_ORDER + i] = v;
xq[n] = sat16(rshift_round(smulww(v, gain_q10), 8)) as i16;
nf.copy_within(0..MAX_LPC_ORDER - 1, 1);
nf[0] = xq[n] as f64 - x[n] as f64;
}
s_lpc.copy_within(sub..sub + MAX_LPC_ORDER, 0);
}
st.s_lpc_q14.copy_from_slice(&s_lpc[..MAX_LPC_ORDER]);
st.out_buf[..l].copy_from_slice(&xq);
(Coded { gain_ix, gains_q16, seed, pulses }, xq)
}
/// Writes the frame's indices and pulses (`decode_indices` + `decode_pulses` order, coded
/// independently).
fn write(&self, enc: &mut RangeEncoder, p: &Params, c: &Coded) {
let ix = p.signal_type * 2 + p.quant_offset_type;
if p.vad {
enc.icdf(ix - 2, &TYPE_ACTIVE, 8);
} else {
enc.icdf(ix, &TYPE_INACTIVE, 8);
}
enc.icdf((c.gain_ix[0] >> 3) as usize, GAIN_MSB[p.signal_type], 8);
enc.icdf((c.gain_ix[0] & 7) as usize, &UNIFORM8, 8);
for k in 1..4 {
enc.icdf(c.gain_ix[k] as usize, &DELTA_GAIN, 8);
}
let wb = self.order == 16;
let i1 = p.nlsf_ix[0] as usize;
enc.icdf(i1, NLSF_STAGE1[(p.signal_type >> 1) + if wb { 2 } else { 0 }], 8);
for i in 0..self.order {
let cb = if wb { 8 + NLSF_SEL_WB[i1][i] as usize } else { NLSF_SEL_NBMB[i1][i] as usize };
let v = p.nlsf_ix[i + 1] + 4;
if v <= 0 {
enc.icdf(0, NLSF_STAGE2[cb], 8);
enc.icdf((-v) as usize, &NLSF_EXT, 8);
} else if v >= 8 {
enc.icdf(8, NLSF_STAGE2[cb], 8);
enc.icdf((v - 8) as usize, &NLSF_EXT, 8);
} else {
enc.icdf(v as usize, NLSF_STAGE2[cb], 8);
}
}
// No NLSF interpolation.
enc.icdf(4, &NLSF_INTERP, 8);
if p.signal_type == 2 {
let half = (self.fs_khz >> 1) as i32;
enc.icdf((p.lag_index / half) as usize, &PITCH_HIGH, 8);
let low = match self.fs_khz {
8 => PITCH_LOW[0],
12 => PITCH_LOW[1],
_ => PITCH_LOW[2],
};
enc.icdf((p.lag_index % half) as usize, low, 8);
enc.icdf(p.contour, if self.fs_khz == 8 { PITCH_CONTOUR[1] } else { PITCH_CONTOUR[3] }, 8);
enc.icdf(p.per_index, &LTP_PERIODICITY, 8);
for k in 0..4 {
enc.icdf(p.ltp_index[k], LTP_FILTER[p.per_index], 8);
}
// LTP scaling index 0 (15565 in Q14), matching `quantize`.
enc.icdf(0, &LTP_SCALE, 8);
}
enc.icdf(c.seed, &UNIFORM4, 8);
pulses::encode(enc, &c.pulses, p.signal_type, p.quant_offset_type);
}
}
/// Quantises a stereo prediction weight (Q13) to the decoder's grid (§4.2.7.1): returns the
/// (UNIFORM3, UNIFORM5, joint) indices and the reconstructed value.
fn quant_pred(w: f64) -> ([usize; 3], i32) {
let mut best = (f64::MAX, [0usize; 3], 0i32);
for i in 0..15 {
let low = STEREO_WEIGHTS_Q13[i];
let step = smulwb(STEREO_WEIGHTS_Q13[i + 1] - low, 6554);
for j in 0..5 {
let v = low + step * (2 * j as i32 + 1);
let e = (v as f64 - w).abs();
if e < best.0 {
best = (e, [i % 3, j, i / 3], v);
}
}
}
(best.1, best.2)
}
/// Per-frame preparation shared by every rate-loop trial.
struct Prep {
x: [Vec<f32>; 2],
params: [Params; 2],
mid_only: bool,
pred_ix: [[usize; 3]; 2],
pred_q13: [i32; 2],
}
/// SILK encoder for one mono or stereo (mid/side) stream.
#[derive(Clone)]
pub struct SilkEncoder {
fs_khz: usize,
channels: usize,
ch: [Channel; 2],
pred_prev_q13: [i32; 2],
prev_mid_only: bool,
/// Last mid sample of the previous frame (for the stereo predictor's low-pass).
prev_mid: f32,
frame: usize,
/// Gain offset chosen for the previous frame (search starting point).
pub last_offset: f64,
}
impl SilkEncoder {
pub fn new(fs_khz: usize, channels: usize) -> Self {
SilkEncoder {
fs_khz,
channels,
ch: [Channel::new(fs_khz), Channel::new(fs_khz)],
pred_prev_q13: [0; 2],
prev_mid_only: false,
prev_mid: 0.0,
frame: 0,
last_offset: 1.0,
}
}
/// Samples per channel of one 20 ms frame at the internal rate.
pub fn frame_len(&self) -> usize {
20 * self.fs_khz
}
fn prepare(&self, input: &[Vec<f32>]) -> Prep {
let l = self.frame_len();
if self.channels == 1 {
let x = input[0].clone();
let p = self.ch[0].analyze(&x);
return Prep { x: [x, vec![0.0; l]], params: [p, Params::default()], mid_only: true, pred_ix: [[0; 3]; 2], pred_q13: [0; 2] };
}
let mid: Vec<f32> = (0..l).map(|i| 0.5 * (input[0][i] + input[1][i])).collect();
let side: Vec<f32> = (0..l).map(|i| 0.5 * (input[0][i] - input[1][i])).collect();
let lp = |j: usize| -> f64 {
let a = if j == 0 { self.prev_mid } else { mid[j - 1] };
let c = if j + 1 < l { mid[j + 1] } else { mid[l - 1] };
(a as f64 + 2.0 * mid[j] as f64 + c as f64) / 4.0
};
// Least-squares prediction of the side from the low-passed mid and the mid.
let (mut r00, mut r01, mut r11, mut c0, mut c1, mut ss, mut mm) = (0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64);
for j in 0..l {
let (a, b, s) = (lp(j), mid[j] as f64, side[j] as f64);
r00 += a * a;
r01 += a * b;
r11 += b * b;
c0 += a * s;
c1 += b * s;
ss += s * s;
mm += b * b;
}
// Ridge regularisation: the low-passed mid and the mid are nearly collinear.
let lam = 0.02 * 0.5 * (r00 + r11) + 1e-9;
let (a00, a11) = (r00 + lam, r11 + lam);
let det = a00 * a11 - r01 * r01;
let (p0, p1) = ((c0 * a11 - c1 * r01) / det, (c1 * a00 - c0 * r01) / det);
let lim = 13732.0;
let quant = |p0: f64, p1: f64| {
let (ix0, q0) = quant_pred(((p0 + p1) * 8192.0).clamp(-lim, lim));
let (ix1, q1) = quant_pred((p1 * 8192.0).clamp(-lim, lim));
([ix0, ix1], [q0 - q1, q1])
};
// Side residual with the decoder's interpolated predictor.
let fs = self.fs_khz as i32;
let denom_q16 = (1 << 16) / (8 * fs);
let interp = 8 * self.fs_khz;
let residual = |pred_q13: [i32; 2]| -> Vec<f32> {
let d0 = rshift_round((pred_q13[0] - self.pred_prev_q13[0]) as i16 as i32 * denom_q16 as i16 as i32, 16);
let d1 = rshift_round((pred_q13[1] - self.pred_prev_q13[1]) as i16 as i32 * denom_q16 as i16 as i32, 16);
let (mut pp0, mut pp1) = (self.pred_prev_q13[0] + d0, self.pred_prev_q13[1] + d1);
(0..l)
.map(|j| {
if j + 1 < interp {
pp0 += d0;
pp1 += d1;
} else {
pp0 = pred_q13[0];
pp1 = pred_q13[1];
}
(side[j] as f64 - (pp0 as f64 * lp(j) + pp1 as f64 * mid[j] as f64) / 8192.0) as f32
})
.collect()
};
let energy = |v: &[f32]| -> f64 { v.iter().map(|&x| (x as f64).powi(2)).sum() };
let (mut pred_ix, mut pred_q13) = quant(p0, p1);
let mut res = residual(pred_q13);
let mut res_e = energy(&res);
// Never predict worse than not predicting at all.
if res_e > ss {
(pred_ix, pred_q13) = quant(0.0, 0.0);
res = residual(pred_q13);
res_e = energy(&res);
}
let res_rms = (res_e / l as f64).sqrt();
let mid_rms = (mm / l as f64).sqrt();
let mid_only = res_rms < 0.01 * mid_rms + 1.0;
let pm = self.ch[0].analyze(&mid);
let ps = if mid_only { Params::default() } else { self.ch[1].analyze(&res) };
Prep { x: [mid, res], params: [pm, ps], mid_only, pred_ix, pred_q13 }
}
/// Codes one frame into `enc` with gain offset `offset`.
fn code(&mut self, prep: &Prep, offset: f64, enc: &mut RangeEncoder) {
let stereo = self.channels == 2;
let side_vad = stereo && !prep.mid_only && prep.params[1].vad;
enc.bit_logp(prep.params[0].vad, 1);
enc.bit_logp(false, 1);
if stereo {
enc.bit_logp(side_vad, 1);
enc.bit_logp(false, 1);
let [a, b] = prep.pred_ix;
enc.icdf(5 * a[2] + b[2], &STEREO_PRED_JOINT, 8);
enc.icdf(a[0], &UNIFORM3, 8);
enc.icdf(a[1], &UNIFORM5, 8);
enc.icdf(b[0], &UNIFORM3, 8);
enc.icdf(b[1], &UNIFORM5, 8);
if !side_vad {
enc.icdf(prep.mid_only as usize, &MID_ONLY, 8);
}
}
let seed = self.frame & 3;
let (cm, _) = self.ch[0].quantize(&prep.params[0], &prep.x[0], offset, None, seed);
self.ch[0].write(enc, &prep.params[0], &cm);
if stereo {
if !prep.mid_only {
if self.prev_mid_only {
self.ch[1].dec.reset_side();
}
let (cs, _) = self.ch[1].quantize(&prep.params[1], &prep.x[1], offset, Some(&cm.gains_q16), seed);
self.ch[1].write(enc, &prep.params[1], &cs);
}
self.pred_prev_q13 = prep.pred_q13;
self.prev_mid_only = prep.mid_only;
}
}
/// Encodes one 20 ms frame (`input`: one internal-rate, 16-bit-scale vector per channel)
/// into `enc`, aiming at `target_bits` total and never exceeding `max_bits`.
pub fn encode(&mut self, input: &[Vec<f32>], target_bits: f64, max_bits: i32, enc: &mut RangeEncoder) {
let prep = self.prepare(input);
let start_bits = enc.tell();
let trial = |s: &Self, o: f64| -> i32 {
let mut t = s.clone();
let mut e = enc.clone();
t.code(&prep, o, &mut e);
e.tell() - start_bits
};
// Bisection on the gain offset: the smallest offset (finest quantisation) within budget.
let target = target_bits.min(max_bits as f64);
// Bits are not strictly monotonic in the offset (very coarse quantisation feeds large
// errors back), so walk an integer grid from the previous frame's offset to find a
// bracket [infeasible, feasible], then bisect inside it.
const MIN_O: f64 = -3.0;
const MAX_O: f64 = 6.0;
let mut best_any = (i32::MAX, MAX_O);
let mut probe = |o: f64| -> bool {
let b = trial(self, o);
if b < best_any.0 {
best_any = (b, o);
}
b as f64 <= target
};
let g0 = self.last_offset.round().clamp(MIN_O, MAX_O);
let mut bracket = None;
if probe(g0) {
let mut hi = g0;
loop {
if hi - 1.0 < MIN_O {
bracket = Some((None, hi));
break;
}
if probe(hi - 1.0) {
hi -= 1.0;
} else {
bracket = Some((Some(hi - 1.0), hi));
break;
}
}
} else {
let mut lo = g0;
while lo + 1.0 <= MAX_O {
if probe(lo + 1.0) {
bracket = Some((Some(lo), lo + 1.0));
break;
}
lo += 1.0;
}
}
let mut o = match bracket {
Some((Some(mut lo), mut hi)) => {
for _ in 0..4 {
let m = 0.5 * (lo + hi);
if probe(m) {
hi = m;
} else {
lo = m;
}
}
hi
}
Some((None, hi)) => hi,
None => best_any.1,
};
// Hard limit (hybrid packets must leave room for CELT).
while trial(self, o) > max_bits && o < 12.0 {
o += 1.0;
}
self.last_offset = o;
self.code(&prep, o, enc);
for ch in 0..self.channels {
self.ch[ch].push_history(&prep.x[ch]);
}
if self.channels == 2 {
self.prev_mid = prep.x[0][self.frame_len() - 1];
}
self.frame += 1;
}
}
+131
View File
@@ -0,0 +1,131 @@
//! NLSF quantisation (encoder side of RFC 6716 §4.2.7.5): stage-1 vector codebook preselection
//! and a backwards-predictive stage-2 scalar quantiser that reproduces the decoder's residual
//! dequantisation exactly. The search strategy is our own.
use super::lpc::{MAX_LPC_ORDER, nlsf_decode, smlawb, sqrt_approx};
use super::tables::{
NLSF_CB_NBMB, NLSF_CB_WB, NLSF_EXT, NLSF_PRED_NBMB, NLSF_PRED_WB, NLSF_PSEL_NBMB, NLSF_PSEL_WB, NLSF_SEL_NBMB, NLSF_SEL_WB, NLSF_STAGE1,
NLSF_STAGE2,
};
/// Number of stage-1 candidates refined with stage 2.
const SURVIVORS: usize = 8;
/// Rate weight (weighted squared Q15 error per bit).
const LAMBDA: f64 = 60.0;
fn icdf_bits(icdf: &[u8], s: usize) -> f64 {
let hi = if s == 0 { 256 } else { icdf[s - 1] as u32 };
let p = (hi - icdf[s] as u32).max(1) as f64 / 256.0;
-p.log2()
}
/// Bits of a stage-2 index `q` (-10..=10) with stage-2 table `tab`.
fn stage2_bits(tab: &[u8], q: i32) -> f64 {
let v = q + 4;
if v <= 0 {
icdf_bits(tab, 0) + icdf_bits(&NLSF_EXT, (-v) as usize)
} else if v >= 8 {
icdf_bits(tab, 8) + icdf_bits(&NLSF_EXT, (v - 8) as usize)
} else {
icdf_bits(tab, v as usize)
}
}
/// Quantises target NLSFs (`target`, Q15 as f64, ascending) and returns the 1 + order indices
/// in the decoder's layout (`[stage1, stage2...]`) together with the reconstructed NLSFs.
pub fn quantize(target: &[f64], signal_type: usize, wb: bool) -> ([i32; 17], [i32; MAX_LPC_ORDER]) {
let order = target.len();
// Laroia-style weights.
let w: Vec<f64> = (0..order)
.map(|k| {
let prev = if k == 0 { 0.0 } else { target[k - 1] };
let next = if k + 1 == order { 32768.0 } else { target[k + 1] };
1.0 / (target[k] - prev).max(50.0) + 1.0 / (next - target[k]).max(50.0)
})
.collect();
let s1_tab = NLSF_STAGE1[(signal_type >> 1) + if wb { 2 } else { 0 }];
let n_cb = 32;
let mut cands: Vec<(f64, usize)> = (0..n_cb)
.map(|i1| {
let cb: &[i32] = if wb { &NLSF_CB_WB[i1] } else { &NLSF_CB_NBMB[i1] };
let d: f64 = (0..order).map(|k| w[k] * ((cb[k] << 7) as f64 - target[k]).powi(2)).sum();
(d + LAMBDA * icdf_bits(s1_tab, i1), i1)
})
.collect();
cands.sort_by(|a, b| a.0.total_cmp(&b.0));
let qstep = if wb { 9830 } else { 11796 };
let mut best = (f64::MAX, [0i32; 17], [0i32; MAX_LPC_ORDER]);
for &(_, i1) in cands.iter().take(SURVIVORS) {
let cb1: &[i32] = if wb { &NLSF_CB_WB[i1] } else { &NLSF_CB_NBMB[i1] };
let mut ix = [0i32; 17];
ix[0] = i1 as i32;
let mut bits = icdf_bits(s1_tab, i1);
// Exact decoder weights per coefficient.
let w_q9: Vec<i32> = (0..order)
.map(|k| {
let prev = if k == 0 { 0 } else { cb1[k - 1] };
let next = if k + 1 == order { 256 } else { cb1[k + 1] };
let ww = (1024 / (cb1[k] - prev).max(1) + 1024 / (next - cb1[k]).max(1)).min(32767);
sqrt_approx(ww << 16).max(1)
})
.collect();
let mut out_q10 = 0i32;
for i in (0..order).rev() {
let pred = if i + 1 < order {
let sel = if wb { NLSF_PSEL_WB[i1][i] } else { NLSF_PSEL_NBMB[i1][i] } as usize;
let wp = if wb { NLSF_PRED_WB[sel][i] } else { NLSF_PRED_NBMB[sel][i] };
(out_q10 * wp) >> 8
} else {
0
};
let cb = if wb { 8 + NLSF_SEL_WB[i1][i] as usize } else { NLSF_SEL_NBMB[i1][i] as usize };
let tab = NLSF_STAGE2[cb];
let want_q10 = (target[i] - (cb1[i] << 7) as f64) * w_q9[i] as f64 / 16384.0;
let step = qstep as f64 * 1024.0 / 65536.0;
let q0 = ((want_q10 - pred as f64) / step).round() as i32;
let mut bq = (f64::MAX, 0i32, 0i32);
for q in (q0 - 1).max(-10)..=(q0 + 1).min(10) {
let mut v = q << 10;
if v > 0 {
v -= 102;
} else if v < 0 {
v += 102;
}
let o = smlawb(pred, v, qstep);
let nl = (cb1[i] << 7) as f64 + ((o << 14) / w_q9[i]) as f64;
let cost = w[i] * (nl - target[i]).powi(2) + LAMBDA * stage2_bits(tab, q);
if cost < bq.0 {
bq = (cost, q, o);
}
}
ix[i + 1] = bq.1;
out_q10 = bq.2;
bits += stage2_bits(tab, bq.1);
}
let nlsf = nlsf_decode(&ix, wb);
let d: f64 = (0..order).map(|k| w[k] * (nlsf[k] as f64 - target[k]).powi(2)).sum();
let cost = d + LAMBDA * bits;
if cost < best.0 {
best = (cost, ix, nlsf);
}
}
(best.1, best.2)
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn quantised_nlsfs_track_the_target() {
for wb in [false, true] {
let order = if wb { 16 } else { 10 };
// A smooth, plausible NLSF vector.
let target: Vec<f64> = (0..order).map(|k| 32768.0 * (k as f64 + 0.8 + 0.3 * (k as f64).sin()) / (order as f64 + 1.0)).collect();
let (ix, nlsf) = quantize(&target, 1, wb);
assert!(ix[1..=order].iter().all(|v| (-10..=10).contains(v)));
let max_err = (0..order).map(|k| (nlsf[k] as f64 - target[k]).abs()).fold(0.0, f64::max);
assert!(max_err < 1200.0, "wb={wb} max err {max_err}: {:?} vs {target:?}", &nlsf[..order]);
}
}
}
@@ -0,0 +1,113 @@
//! Excitation coding (encoder side of RFC 6716 §4.2.7.8): rate level, per-block pulse counts
//! with LSB escapes, the recursive shell code, LSBs and signs, written in exactly the order the
//! decoder reads them.
use super::tables::{EXC_LSB, EXC_SIGN, PULSE_COUNT, RATE_LEVEL, SHELL_2, SHELL_4, SHELL_8, SHELL_16};
use crate::opus::encoder::range::RangeEncoder;
/// Largest pulse magnitude the quantiser produces (keeps every 16-sample block codable).
pub const MAX_PULSE: i32 = 500;
pub fn icdf_bits(icdf: &[u8], s: usize) -> f32 {
let hi = if s == 0 { 256 } else { icdf[s - 1] as u32 };
let p = (hi - icdf[s] as u32).max(1) as f32 / 256.0;
-p.log2()
}
fn count_table(nls: usize) -> &'static [u8] {
if nls == 10 { PULSE_COUNT[10] } else { PULSE_COUNT[9] }
}
fn shell(enc: &mut RangeEncoder, v: &[i32], count: i32) {
let n = v.len();
if n == 1 || count == 0 {
return;
}
let tab = match n {
16 => SHELL_16[(count - 1) as usize],
8 => SHELL_8[(count - 1) as usize],
4 => SHELL_4[(count - 1) as usize],
_ => SHELL_2[(count - 1) as usize],
};
let (a, b) = v.split_at(n / 2);
let left: i32 = a.iter().sum();
enc.icdf(left as usize, tab, 8);
shell(enc, a, left);
shell(enc, b, count - left);
}
/// Codes the frame's pulses (`len` a multiple of 16).
pub fn encode(enc: &mut RangeEncoder, pulses: &[i32], signal_type: usize, quant_offset_type: usize) {
let nblk = pulses.len() / 16;
let mut nls = vec![0usize; nblk];
let mut sums = vec![0i32; nblk];
let mut hi = vec![0i32; pulses.len()];
for b in 0..nblk {
let blk = &pulses[b * 16..(b + 1) * 16];
let mut s = 0usize;
loop {
let sum: i32 = blk.iter().map(|p| p.abs() >> s).sum();
if sum <= 16 || s == 10 {
sums[b] = sum.min(16);
break;
}
s += 1;
}
nls[b] = s;
for (h, p) in hi[b * 16..(b + 1) * 16].iter_mut().zip(blk) {
*h = p.abs() >> s;
}
}
// Rate level with the cheapest pulse-count symbols.
let rl_tab = RATE_LEVEL[signal_type >> 1];
let mut best = (f32::MAX, 0usize);
for rl in 0..9 {
let mut bits = icdf_bits(rl_tab, rl);
for b in 0..nblk {
bits += if nls[b] == 0 { icdf_bits(PULSE_COUNT[rl], sums[b] as usize) } else { icdf_bits(PULSE_COUNT[rl], 17) };
}
if bits < best.0 {
best = (bits, rl);
}
}
let rl = best.1;
enc.icdf(rl, rl_tab, 8);
for b in 0..nblk {
if nls[b] == 0 {
enc.icdf(sums[b] as usize, PULSE_COUNT[rl], 8);
} else {
enc.icdf(17, PULSE_COUNT[rl], 8);
for j in 1..nls[b] {
enc.icdf(17, count_table(j), 8);
}
enc.icdf(sums[b] as usize, count_table(nls[b]), 8);
}
}
for b in 0..nblk {
if sums[b] > 0 {
shell(enc, &hi[b * 16..(b + 1) * 16], sums[b]);
}
}
for b in 0..nblk {
if nls[b] > 0 {
for p in &pulses[b * 16..(b + 1) * 16] {
let a = p.abs();
for j in (0..nls[b]).rev() {
enc.icdf(((a >> j) & 1) as usize, &EXC_LSB, 8);
}
}
}
}
let base = 7 * (quant_offset_type + (signal_type << 1));
for b in 0..nblk {
let p = sums[b] | ((nls[b] as i32) << 5);
if p > 0 {
let tab = EXC_SIGN[base + ((p & 0x1F) as usize).min(6)];
for &q in &pulses[b * 16..(b + 1) * 16] {
if q != 0 {
enc.icdf((q > 0) as usize, tab, 8);
}
}
}
}
}
@@ -0,0 +1,543 @@
//! SILK probability tables and codebooks, transcribed from RFC 6716 §4.2 (Tables 4–53).
//!
//! PDFs are written as in the RFC (`{...}/256`) and converted at compile time to the inverse-CDF
//! form consumed by the range coder.
//!
//! Copied from FilmCraft's first-party, spec-derived `filmcraft-opus` decoder (MIT OR
//! Apache-2.0) so encoder and decoder share exactly the same values.
/// Converts a PDF (summing to 256) into an inverse CDF (`256 - cdf[i+1]`).
pub const fn icdf<const N: usize>(pdf: &[u8; N]) -> [u8; N] {
let mut out = [0u8; N];
let mut acc: u32 = 0;
let mut i = 0;
while i < N {
acc += pdf[i] as u32;
out[i] = (256 - acc) as u8;
i += 1;
}
out
}
pub const STEREO_PRED_JOINT: [u8; 25] = icdf(&[7, 2, 1, 1, 1, 10, 24, 8, 1, 1, 3, 23, 92, 23, 3, 1, 1, 8, 24, 10, 1, 1, 1, 2, 7]);
pub const UNIFORM3: [u8; 3] = icdf(&[85, 86, 85]);
pub const UNIFORM5: [u8; 5] = icdf(&[51, 51, 52, 51, 51]);
pub const STEREO_WEIGHTS_Q13: [i32; 16] = [-13732, -10050, -8266, -7526, -6500, -5000, -2950, -820, 820, 2950, 5000, 6500, 7526, 8266, 10050, 13732];
pub const MID_ONLY: [u8; 2] = icdf(&[192, 64]);
pub const TYPE_INACTIVE: [u8; 2] = icdf(&[26, 230]);
pub const TYPE_ACTIVE: [u8; 4] = icdf(&[24, 74, 148, 10]);
pub const GAIN_MSB: [&[u8]; 3] =
[&icdf::<8>(&[32, 112, 68, 29, 12, 1, 1, 1]), &icdf::<8>(&[2, 17, 45, 60, 62, 47, 19, 4]), &icdf::<8>(&[1, 3, 26, 71, 94, 50, 9, 2])];
pub const UNIFORM8: [u8; 8] = icdf(&[32, 32, 32, 32, 32, 32, 32, 32]);
pub const DELTA_GAIN: [u8; 41] =
icdf(&[6, 5, 11, 31, 132, 21, 8, 4, 3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]);
pub const NLSF_STAGE1: [&[u8]; 4] = [
&icdf::<32>(&[44, 34, 30, 19, 21, 12, 11, 3, 3, 2, 16, 2, 2, 1, 5, 2, 1, 3, 3, 1, 1, 2, 2, 2, 3, 1, 9, 9, 2, 7, 2, 1]),
&icdf::<32>(&[1, 10, 1, 8, 3, 8, 8, 14, 13, 14, 1, 14, 12, 13, 11, 11, 12, 11, 10, 10, 11, 8, 9, 8, 7, 8, 1, 1, 6, 1, 6, 5]),
&icdf::<32>(&[31, 21, 3, 17, 1, 8, 17, 4, 1, 18, 16, 4, 2, 3, 1, 10, 1, 3, 16, 11, 16, 2, 2, 3, 2, 11, 1, 4, 9, 8, 7, 3]),
&icdf::<32>(&[1, 4, 16, 5, 18, 11, 5, 14, 15, 1, 3, 12, 13, 14, 14, 6, 14, 12, 2, 6, 1, 12, 12, 11, 10, 3, 10, 5, 1, 1, 1, 3]),
];
pub const NLSF_STAGE2: [&[u8]; 16] = [
&icdf::<9>(&[1, 1, 1, 15, 224, 11, 1, 1, 1]),
&icdf::<9>(&[1, 1, 2, 34, 183, 32, 1, 1, 1]),
&icdf::<9>(&[1, 1, 4, 42, 149, 55, 2, 1, 1]),
&icdf::<9>(&[1, 1, 8, 52, 123, 61, 8, 1, 1]),
&icdf::<9>(&[1, 3, 16, 53, 101, 74, 6, 1, 1]),
&icdf::<9>(&[1, 3, 17, 55, 90, 73, 15, 1, 1]),
&icdf::<9>(&[1, 7, 24, 53, 74, 67, 26, 3, 1]),
&icdf::<9>(&[1, 1, 18, 63, 78, 58, 30, 6, 1]),
&icdf::<9>(&[1, 1, 1, 9, 232, 9, 1, 1, 1]),
&icdf::<9>(&[1, 1, 2, 28, 186, 35, 1, 1, 1]),
&icdf::<9>(&[1, 1, 3, 42, 152, 53, 2, 1, 1]),
&icdf::<9>(&[1, 1, 10, 49, 126, 65, 2, 1, 1]),
&icdf::<9>(&[1, 4, 19, 48, 100, 77, 5, 1, 1]),
&icdf::<9>(&[1, 1, 14, 54, 100, 72, 12, 1, 1]),
&icdf::<9>(&[1, 1, 15, 61, 87, 61, 25, 4, 1]),
&icdf::<9>(&[1, 7, 21, 50, 77, 81, 17, 1, 1]),
];
pub const NLSF_SEL_NBMB: [[u8; 10]; 32] = [
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 3, 1, 2, 2, 1, 2, 1, 1, 1],
[2, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 2, 2, 2, 2, 1, 2, 1, 1, 1],
[2, 3, 3, 3, 3, 2, 2, 2, 2, 2],
[0, 5, 3, 3, 2, 2, 2, 2, 1, 1],
[0, 2, 2, 2, 2, 2, 2, 2, 2, 1],
[2, 3, 6, 4, 4, 4, 5, 4, 5, 5],
[2, 4, 5, 5, 4, 5, 4, 6, 4, 4],
[2, 4, 4, 7, 4, 5, 4, 5, 5, 4],
[4, 3, 3, 3, 2, 3, 2, 2, 2, 2],
[1, 5, 5, 6, 4, 5, 4, 5, 5, 5],
[2, 7, 4, 6, 5, 5, 5, 5, 5, 5],
[2, 7, 5, 5, 5, 5, 5, 6, 5, 4],
[3, 3, 5, 4, 4, 5, 4, 5, 4, 4],
[2, 3, 3, 5, 5, 4, 4, 4, 4, 4],
[2, 4, 4, 6, 4, 5, 4, 5, 5, 5],
[2, 5, 4, 6, 5, 5, 5, 4, 5, 4],
[2, 7, 4, 5, 4, 5, 4, 5, 5, 5],
[2, 5, 4, 6, 7, 6, 5, 6, 5, 4],
[3, 6, 7, 4, 6, 5, 5, 6, 4, 5],
[2, 7, 6, 4, 4, 4, 5, 4, 5, 5],
[4, 5, 5, 4, 6, 6, 5, 6, 5, 4],
[2, 5, 5, 6, 5, 6, 4, 6, 4, 4],
[4, 5, 5, 5, 3, 7, 4, 5, 5, 4],
[2, 3, 4, 5, 5, 6, 4, 5, 5, 4],
[2, 3, 2, 3, 3, 4, 2, 3, 3, 3],
[1, 1, 2, 2, 2, 2, 2, 3, 2, 2],
[4, 5, 5, 6, 6, 6, 5, 6, 4, 5],
[3, 5, 5, 4, 4, 4, 4, 3, 3, 2],
[2, 5, 3, 7, 5, 5, 4, 4, 5, 4],
[4, 4, 5, 4, 5, 6, 5, 6, 5, 4],
];
pub const NLSF_SEL_WB: [[u8; 16]; 32] = [
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[2, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 1, 1, 1, 0, 3],
[2, 5, 5, 3, 7, 4, 4, 5, 2, 5, 4, 5, 5, 4, 3, 3],
[0, 2, 1, 2, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1],
[0, 6, 5, 4, 6, 4, 7, 5, 4, 4, 4, 5, 5, 4, 4, 3],
[0, 3, 5, 5, 4, 3, 3, 5, 3, 3, 3, 3, 3, 3, 2, 4],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 2, 6, 3, 7, 2, 5, 3, 4, 5, 5, 4, 3, 3, 2, 3],
[0, 6, 2, 6, 6, 4, 5, 4, 6, 5, 4, 4, 5, 3, 3, 3],
[2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[2, 2, 3, 4, 5, 3, 3, 3, 3, 3, 3, 3, 2, 2, 1, 3],
[2, 2, 3, 3, 4, 3, 3, 3, 3, 3, 3, 3, 3, 2, 1, 3],
[3, 4, 4, 4, 6, 4, 4, 5, 3, 5, 4, 4, 5, 4, 3, 4],
[0, 6, 4, 5, 4, 7, 5, 2, 6, 5, 7, 4, 4, 3, 5, 3],
[0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0],
[1, 6, 5, 7, 5, 4, 5, 3, 4, 5, 4, 4, 4, 3, 3, 4],
[1, 3, 3, 4, 4, 3, 3, 5, 2, 3, 3, 5, 5, 5, 3, 4],
[2, 3, 3, 2, 2, 2, 3, 2, 1, 2, 1, 2, 1, 1, 1, 4],
[0, 2, 3, 5, 3, 3, 2, 2, 2, 1, 1, 0, 0, 0, 0, 0],
[3, 4, 3, 5, 3, 3, 2, 2, 1, 1, 1, 1, 1, 2, 2, 4],
[2, 6, 3, 7, 7, 4, 5, 4, 5, 3, 5, 3, 3, 2, 3, 3],
[2, 3, 5, 6, 6, 3, 5, 3, 4, 4, 3, 3, 3, 3, 2, 4],
[1, 3, 3, 4, 4, 4, 4, 3, 5, 5, 5, 3, 1, 1, 1, 1],
[2, 5, 3, 6, 6, 4, 7, 4, 4, 5, 3, 4, 4, 3, 3, 3],
[0, 6, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 6, 6, 3, 5, 2, 5, 5, 3, 4, 4, 7, 7, 4, 4, 4],
[3, 3, 7, 3, 5, 4, 3, 3, 3, 2, 2, 3, 3, 3, 2, 3],
[0, 0, 1, 0, 0, 0, 2, 1, 2, 1, 1, 2, 2, 2, 1, 1],
[0, 3, 2, 5, 3, 3, 2, 3, 2, 1, 0, 0, 1, 0, 0, 1],
[3, 5, 5, 4, 7, 5, 3, 3, 2, 3, 2, 2, 1, 0, 1, 0],
[2, 3, 5, 3, 4, 3, 3, 3, 2, 1, 2, 6, 4, 0, 0, 0],
];
pub const NLSF_EXT: [u8; 7] = icdf(&[156, 60, 24, 9, 4, 2, 1]);
pub const NLSF_PRED_NBMB: [[i32; 9]; 2] = [[179, 138, 140, 148, 151, 149, 153, 151, 163], [116, 67, 82, 59, 92, 72, 100, 89, 92]];
pub const NLSF_PRED_WB: [[i32; 15]; 2] = [
[175, 148, 160, 176, 178, 173, 174, 164, 177, 174, 196, 182, 198, 192, 182],
[68, 62, 66, 60, 72, 117, 85, 90, 118, 136, 151, 142, 160, 142, 155],
];
pub const NLSF_PSEL_NBMB: [[u8; 9]; 32] = [
[0, 1, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 1, 1, 0, 0, 0, 0, 1, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 1, 1, 0, 0, 0, 1, 0],
[0, 1, 1, 0, 0, 1, 1, 0, 0],
[0, 0, 1, 1, 0, 1, 0, 1, 1],
[0, 0, 1, 1, 0, 0, 1, 1, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 1, 1, 1, 1, 1, 0],
[0, 1, 0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 0],
[1, 0, 1, 1, 0, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 0, 1, 0],
[0, 0, 1, 1, 0, 1, 0, 1, 0],
[0, 0, 1, 1, 1, 0, 1, 1, 1],
[0, 1, 1, 0, 0, 1, 1, 1, 0],
[0, 0, 0, 1, 1, 1, 0, 1, 0],
[0, 1, 1, 0, 0, 1, 0, 1, 0],
[0, 1, 1, 0, 0, 0, 1, 1, 0],
[0, 0, 0, 0, 0, 1, 1, 1, 1],
[0, 0, 1, 1, 0, 0, 0, 1, 1],
[0, 0, 0, 1, 0, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 1, 1, 0, 1, 0],
[1, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 1, 0, 1, 0, 1],
[1, 0, 1, 1, 0, 1, 1, 1, 1],
];
pub const NLSF_PSEL_WB: [[u8; 15]; 32] = [
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
[0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0],
[0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1],
[0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0],
[0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0],
[0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1],
[0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
[0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0],
[0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0],
[0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1],
[0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0],
];
pub const NLSF_CB_NBMB: [[i32; 10]; 32] = [
[12, 35, 60, 83, 108, 132, 157, 180, 206, 228],
[15, 32, 55, 77, 101, 125, 151, 175, 201, 225],
[19, 42, 66, 89, 114, 137, 162, 184, 209, 230],
[12, 25, 50, 72, 97, 120, 147, 172, 200, 223],
[26, 44, 69, 90, 114, 135, 159, 180, 205, 225],
[13, 22, 53, 80, 106, 130, 156, 180, 205, 228],
[15, 25, 44, 64, 90, 115, 142, 168, 196, 222],
[19, 24, 62, 82, 100, 120, 145, 168, 190, 214],
[22, 31, 50, 79, 103, 120, 151, 170, 203, 227],
[21, 29, 45, 65, 106, 124, 150, 171, 196, 224],
[30, 49, 75, 97, 121, 142, 165, 186, 209, 229],
[19, 25, 52, 70, 93, 116, 143, 166, 192, 219],
[26, 34, 62, 75, 97, 118, 145, 167, 194, 217],
[25, 33, 56, 70, 91, 113, 143, 165, 196, 223],
[21, 34, 51, 72, 97, 117, 145, 171, 196, 222],
[20, 29, 50, 67, 90, 117, 144, 168, 197, 221],
[22, 31, 48, 66, 95, 117, 146, 168, 196, 222],
[24, 33, 51, 77, 116, 134, 158, 180, 200, 224],
[21, 28, 70, 87, 106, 124, 149, 170, 194, 217],
[26, 33, 53, 64, 83, 117, 152, 173, 204, 225],
[27, 34, 65, 95, 108, 129, 155, 174, 210, 225],
[20, 26, 72, 99, 113, 131, 154, 176, 200, 219],
[34, 43, 61, 78, 93, 114, 155, 177, 205, 229],
[23, 29, 54, 97, 124, 138, 163, 179, 209, 229],
[30, 38, 56, 89, 118, 129, 158, 178, 200, 231],
[21, 29, 49, 63, 85, 111, 142, 163, 193, 222],
[27, 48, 77, 103, 133, 158, 179, 196, 215, 232],
[29, 47, 74, 99, 124, 151, 176, 198, 220, 237],
[33, 42, 61, 76, 93, 121, 155, 174, 207, 225],
[29, 53, 87, 112, 136, 154, 170, 188, 208, 227],
[24, 30, 52, 84, 131, 150, 166, 186, 203, 229],
[37, 48, 64, 84, 104, 118, 156, 177, 201, 230],
];
pub const NLSF_CB_WB: [[i32; 16]; 32] = [
[7, 23, 38, 54, 69, 85, 100, 116, 131, 147, 162, 178, 193, 208, 223, 239],
[13, 25, 41, 55, 69, 83, 98, 112, 127, 142, 157, 171, 187, 203, 220, 236],
[15, 21, 34, 51, 61, 78, 92, 106, 126, 136, 152, 167, 185, 205, 225, 240],
[10, 21, 36, 50, 63, 79, 95, 110, 126, 141, 157, 173, 189, 205, 221, 237],
[17, 20, 37, 51, 59, 78, 89, 107, 123, 134, 150, 164, 184, 205, 224, 240],
[10, 15, 32, 51, 67, 81, 96, 112, 129, 142, 158, 173, 189, 204, 220, 236],
[8, 21, 37, 51, 65, 79, 98, 113, 126, 138, 155, 168, 179, 192, 209, 218],
[12, 15, 34, 55, 63, 78, 87, 108, 118, 131, 148, 167, 185, 203, 219, 236],
[16, 19, 32, 36, 56, 79, 91, 108, 118, 136, 154, 171, 186, 204, 220, 237],
[11, 28, 43, 58, 74, 89, 105, 120, 135, 150, 165, 180, 196, 211, 226, 241],
[6, 16, 33, 46, 60, 75, 92, 107, 123, 137, 156, 169, 185, 199, 214, 225],
[11, 19, 30, 44, 57, 74, 89, 105, 121, 135, 152, 169, 186, 202, 218, 234],
[12, 19, 29, 46, 57, 71, 88, 100, 120, 132, 148, 165, 182, 199, 216, 233],
[17, 23, 35, 46, 56, 77, 92, 106, 123, 134, 152, 167, 185, 204, 222, 237],
[14, 17, 45, 53, 63, 75, 89, 107, 115, 132, 151, 171, 188, 206, 221, 240],
[9, 16, 29, 40, 56, 71, 88, 103, 119, 137, 154, 171, 189, 205, 222, 237],
[16, 19, 36, 48, 57, 76, 87, 105, 118, 132, 150, 167, 185, 202, 218, 236],
[12, 17, 29, 54, 71, 81, 94, 104, 126, 136, 149, 164, 182, 201, 221, 237],
[15, 28, 47, 62, 79, 97, 115, 129, 142, 155, 168, 180, 194, 208, 223, 238],
[8, 14, 30, 45, 62, 78, 94, 111, 127, 143, 159, 175, 192, 207, 223, 239],
[17, 30, 49, 62, 79, 92, 107, 119, 132, 145, 160, 174, 190, 204, 220, 235],
[14, 19, 36, 45, 61, 76, 91, 108, 121, 138, 154, 172, 189, 205, 222, 238],
[12, 18, 31, 45, 60, 76, 91, 107, 123, 138, 154, 171, 187, 204, 221, 236],
[13, 17, 31, 43, 53, 70, 83, 103, 114, 131, 149, 167, 185, 203, 220, 237],
[17, 22, 35, 42, 58, 78, 93, 110, 125, 139, 155, 170, 188, 206, 224, 240],
[8, 15, 34, 50, 67, 83, 99, 115, 131, 146, 162, 178, 193, 209, 224, 239],
[13, 16, 41, 66, 73, 86, 95, 111, 128, 137, 150, 163, 183, 206, 225, 241],
[17, 25, 37, 52, 63, 75, 92, 102, 119, 132, 144, 160, 175, 191, 212, 231],
[19, 31, 49, 65, 83, 100, 117, 133, 147, 161, 174, 187, 200, 213, 227, 242],
[18, 31, 52, 68, 88, 103, 117, 126, 138, 149, 163, 177, 192, 207, 223, 239],
[16, 29, 47, 61, 76, 90, 106, 119, 133, 147, 161, 176, 193, 209, 224, 240],
[15, 21, 35, 50, 61, 73, 86, 97, 110, 119, 129, 141, 175, 198, 218, 237],
];
pub const NLSF_MIN_SPACING_NBMB: [i32; 11] = [250, 3, 6, 3, 3, 3, 4, 3, 3, 3, 461];
pub const NLSF_MIN_SPACING_WB: [i32; 17] = [100, 3, 40, 3, 3, 3, 5, 14, 14, 10, 11, 3, 8, 9, 7, 3, 347];
pub const NLSF_INTERP: [u8; 5] = icdf(&[13, 22, 29, 11, 181]);
pub const LSF_ORDER_NBMB: [usize; 10] = [0, 9, 6, 3, 4, 5, 8, 1, 2, 7];
pub const LSF_ORDER_WB: [usize; 16] = [0, 15, 8, 7, 4, 11, 12, 3, 2, 13, 10, 5, 6, 9, 14, 1];
pub const COS_Q12: [i32; 129] = [
4096, 4095, 4091, 4085, 4076, 4065, 4052, 4036, 4017, 3997, 3973, 3948, 3920, 3889, 3857, 3822, 3784, 3745, 3703, 3659, 3613, 3564, 3513, 3461,
3406, 3349, 3290, 3229, 3166, 3102, 3035, 2967, 2896, 2824, 2751, 2676, 2599, 2520, 2440, 2359, 2276, 2191, 2106, 2019, 1931, 1842, 1751, 1660,
1568, 1474, 1380, 1285, 1189, 1093, 995, 897, 799, 700, 601, 501, 401, 301, 201, 101, 0, -101, -201, -301, -401, -501, -601, -700, -799, -897,
-995, -1093, -1189, -1285, -1380, -1474, -1568, -1660, -1751, -1842, -1931, -2019, -2106, -2191, -2276, -2359, -2440, -2520, -2599, -2676, -2751,
-2824, -2896, -2967, -3035, -3102, -3166, -3229, -3290, -3349, -3406, -3461, -3513, -3564, -3613, -3659, -3703, -3745, -3784, -3822, -3857,
-3889, -3920, -3948, -3973, -3997, -4017, -4036, -4052, -4065, -4076, -4085, -4091, -4095, -4096,
];
pub const PITCH_HIGH: [u8; 32] = icdf(&[3, 3, 6, 11, 21, 30, 32, 19, 11, 10, 12, 13, 13, 12, 11, 9, 8, 7, 6, 4, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1]);
pub const PITCH_LOW: [&[u8]; 3] =
[&icdf::<4>(&[64, 64, 64, 64]), &icdf::<6>(&[43, 42, 43, 43, 42, 43]), &icdf::<8>(&[32, 32, 32, 32, 32, 32, 32, 32])];
pub const PITCH_DELTA: [u8; 21] = icdf(&[46, 2, 2, 3, 4, 6, 10, 15, 26, 38, 30, 22, 15, 10, 7, 6, 4, 4, 2, 2, 2]);
pub const PITCH_CONTOUR: [&[u8]; 4] = [
&icdf::<3>(&[143, 50, 63]),
&icdf::<11>(&[68, 12, 21, 17, 19, 22, 30, 24, 17, 16, 10]),
&icdf::<12>(&[91, 46, 39, 19, 14, 12, 8, 7, 6, 5, 5, 4]),
&icdf::<34>(&[33, 22, 18, 16, 15, 14, 14, 13, 13, 10, 9, 9, 8, 6, 6, 6, 5, 4, 4, 4, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1]),
];
pub const CB_NB_10: [[i32; 2]; 3] = [[0, 0], [1, 0], [0, 1]];
pub const CB_NB_20: [[i32; 4]; 11] = [
[0, 0, 0, 0],
[2, 1, 0, -1],
[-1, 0, 1, 2],
[-1, 0, 0, 1],
[-1, 0, 0, 0],
[0, 0, 0, 1],
[0, 0, 1, 1],
[1, 1, 0, 0],
[1, 0, 0, 0],
[0, 0, 0, -1],
[1, 0, 0, -1],
];
pub const CB_WB_10: [[i32; 2]; 12] = [[0, 0], [0, 1], [1, 0], [-1, 1], [1, -1], [-1, 2], [2, -1], [-2, 2], [2, -2], [-2, 3], [3, -2], [-3, 3]];
pub const CB_WB_20: [[i32; 4]; 34] = [
[0, 0, 0, 0],
[0, 0, 1, 1],
[1, 1, 0, 0],
[-1, 0, 0, 0],
[0, 0, 0, 1],
[1, 0, 0, 0],
[-1, 0, 0, 1],
[0, 0, 0, -1],
[-1, 0, 1, 2],
[1, 0, 0, -1],
[-2, -1, 1, 2],
[2, 1, 0, -1],
[-2, 0, 0, 2],
[-2, 0, 1, 3],
[2, 1, -1, -2],
[-3, -1, 1, 3],
[2, 0, 0, -2],
[3, 1, 0, -2],
[-3, -1, 2, 4],
[-4, -1, 1, 4],
[3, 1, -1, -3],
[-4, -1, 2, 5],
[4, 2, -1, -3],
[4, 1, -1, -4],
[-5, -1, 2, 6],
[5, 2, -1, -4],
[-6, -2, 2, 6],
[-5, -2, 2, 5],
[6, 2, -1, -5],
[-7, -2, 3, 8],
[6, 2, -2, -6],
[5, 2, -2, -5],
[8, 3, -2, -7],
[-9, -3, 3, 9],
];
pub const LTP_PERIODICITY: [u8; 3] = icdf(&[77, 80, 99]);
pub const LTP_FILTER: [&[u8]; 3] = [
&icdf::<8>(&[185, 15, 13, 13, 9, 9, 6, 6]),
&icdf::<16>(&[57, 34, 21, 20, 15, 13, 12, 13, 10, 10, 9, 10, 9, 8, 7, 8]),
&icdf::<32>(&[15, 16, 14, 12, 12, 12, 11, 11, 11, 10, 9, 9, 9, 9, 8, 8, 8, 8, 7, 7, 6, 6, 5, 4, 5, 4, 4, 4, 3, 4, 3, 2]),
];
pub const LTP_TAPS_0: [[i32; 5]; 8] = [
[4, 6, 24, 7, 5],
[0, 0, 2, 0, 0],
[12, 28, 41, 13, -4],
[-9, 15, 42, 25, 14],
[1, -2, 62, 41, -9],
[-10, 37, 65, -4, 3],
[-6, 4, 66, 7, -8],
[16, 14, 38, -3, 33],
];
pub const LTP_TAPS_1: [[i32; 5]; 16] = [
[13, 22, 39, 23, 12],
[-1, 36, 64, 27, -6],
[-7, 10, 55, 43, 17],
[1, 1, 8, 1, 1],
[6, -11, 74, 53, -9],
[-12, 55, 76, -12, 8],
[-3, 3, 93, 27, -4],
[26, 39, 59, 3, -8],
[2, 0, 77, 11, 9],
[-8, 22, 44, -6, 7],
[40, 9, 26, 3, 9],
[-7, 20, 101, -7, 4],
[3, -8, 42, 26, 0],
[-15, 33, 68, 2, 23],
[-2, 55, 46, -2, 15],
[3, -1, 21, 16, 41],
];
pub const LTP_TAPS_2: [[i32; 5]; 32] = [
[-6, 27, 61, 39, 5],
[-11, 42, 88, 4, 1],
[-2, 60, 65, 6, -4],
[-1, -5, 73, 56, 1],
[-9, 19, 94, 29, -9],
[0, 12, 99, 6, 4],
[8, -19, 102, 46, -13],
[3, 2, 13, 3, 2],
[9, -21, 84, 72, -18],
[-11, 46, 104, -22, 8],
[18, 38, 48, 23, 0],
[-16, 70, 83, -21, 11],
[5, -11, 117, 22, -8],
[-6, 23, 117, -12, 3],
[3, -8, 95, 28, 4],
[-10, 15, 77, 60, -15],
[-1, 4, 124, 2, -4],
[3, 38, 84, 24, -25],
[2, 13, 42, 13, 31],
[21, -4, 56, 46, -1],
[-1, 35, 79, -13, 19],
[-7, 65, 88, -9, -14],
[20, 4, 81, 49, -29],
[20, 0, 75, 3, -17],
[5, -9, 44, 92, -8],
[1, -3, 22, 69, 31],
[-6, 95, 41, -12, 5],
[39, 67, 16, -4, 1],
[0, -6, 120, 55, -36],
[-13, 44, 122, 4, -24],
[81, 5, 11, 3, 7],
[2, 0, 9, 10, 88],
];
pub const LTP_SCALE: [u8; 3] = icdf(&[128, 64, 64]);
pub const UNIFORM4: [u8; 4] = icdf(&[64, 64, 64, 64]);
pub const RATE_LEVEL: [&[u8]; 2] = [&icdf::<9>(&[15, 51, 12, 46, 45, 13, 33, 27, 14]), &icdf::<9>(&[33, 30, 36, 17, 34, 49, 18, 21, 18])];
pub const PULSE_COUNT: [&[u8]; 11] = [
&icdf::<18>(&[131, 74, 25, 8, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]),
&icdf::<18>(&[58, 93, 60, 23, 7, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]),
&icdf::<18>(&[43, 51, 46, 33, 24, 16, 11, 8, 6, 3, 3, 3, 2, 1, 1, 2, 1, 2]),
&icdf::<18>(&[17, 52, 71, 57, 31, 12, 5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]),
&icdf::<18>(&[6, 21, 41, 53, 49, 35, 21, 11, 6, 3, 2, 2, 1, 1, 1, 1, 1, 1]),
&icdf::<18>(&[7, 14, 22, 28, 29, 28, 25, 20, 17, 13, 11, 9, 7, 5, 4, 4, 3, 10]),
&icdf::<18>(&[2, 5, 14, 29, 42, 46, 41, 31, 19, 11, 6, 3, 2, 1, 1, 1, 1, 1]),
&icdf::<18>(&[1, 2, 4, 10, 19, 29, 35, 37, 34, 28, 20, 14, 8, 5, 4, 2, 2, 2]),
&icdf::<18>(&[1, 2, 2, 5, 9, 14, 20, 24, 27, 28, 26, 23, 20, 15, 11, 8, 6, 15]),
&icdf::<18>(&[1, 1, 1, 6, 27, 58, 56, 39, 25, 14, 10, 6, 3, 3, 2, 1, 1, 2]),
&icdf::<18>(&[2, 1, 6, 27, 58, 56, 39, 25, 14, 10, 6, 3, 3, 2, 1, 1, 2, 0]),
];
pub const SHELL_16: [&[u8]; 16] = [
&icdf::<2>(&[126, 130]),
&icdf::<3>(&[56, 142, 58]),
&icdf::<4>(&[25, 101, 104, 26]),
&icdf::<5>(&[12, 60, 108, 64, 12]),
&icdf::<6>(&[7, 35, 84, 87, 37, 6]),
&icdf::<7>(&[4, 20, 59, 86, 63, 21, 3]),
&icdf::<8>(&[3, 12, 38, 72, 75, 42, 12, 2]),
&icdf::<9>(&[2, 8, 25, 54, 73, 59, 27, 7, 1]),
&icdf::<10>(&[2, 5, 17, 39, 63, 65, 42, 18, 4, 1]),
&icdf::<11>(&[1, 4, 12, 28, 49, 63, 54, 30, 11, 3, 1]),
&icdf::<12>(&[1, 4, 8, 20, 37, 55, 57, 41, 22, 8, 2, 1]),
&icdf::<13>(&[1, 3, 7, 15, 28, 44, 53, 48, 33, 16, 6, 1, 1]),
&icdf::<14>(&[1, 2, 6, 12, 21, 35, 47, 48, 40, 25, 12, 5, 1, 1]),
&icdf::<15>(&[1, 1, 4, 10, 17, 27, 37, 47, 43, 33, 21, 9, 4, 1, 1]),
&icdf::<16>(&[1, 1, 1, 8, 14, 22, 33, 40, 43, 38, 28, 16, 8, 1, 1, 1]),
&icdf::<17>(&[1, 1, 1, 1, 13, 18, 27, 36, 41, 41, 34, 24, 14, 1, 1, 1, 1]),
];
pub const SHELL_8: [&[u8]; 16] = [
&icdf::<2>(&[127, 129]),
&icdf::<3>(&[53, 149, 54]),
&icdf::<4>(&[22, 105, 106, 23]),
&icdf::<5>(&[11, 61, 111, 63, 10]),
&icdf::<6>(&[6, 35, 86, 88, 36, 5]),
&icdf::<7>(&[4, 20, 59, 87, 62, 21, 3]),
&icdf::<8>(&[3, 13, 40, 71, 73, 41, 13, 2]),
&icdf::<9>(&[3, 9, 27, 53, 70, 56, 28, 9, 1]),
&icdf::<10>(&[3, 8, 19, 37, 57, 61, 44, 20, 6, 1]),
&icdf::<11>(&[3, 7, 15, 28, 44, 54, 49, 33, 17, 5, 1]),
&icdf::<12>(&[1, 7, 13, 22, 34, 46, 48, 38, 28, 14, 4, 1]),
&icdf::<13>(&[1, 1, 11, 22, 27, 35, 42, 47, 33, 25, 10, 1, 1]),
&icdf::<14>(&[1, 1, 6, 14, 26, 37, 43, 43, 37, 26, 14, 6, 1, 1]),
&icdf::<15>(&[1, 1, 4, 10, 20, 31, 40, 42, 40, 31, 20, 10, 4, 1, 1]),
&icdf::<16>(&[1, 1, 3, 8, 16, 26, 35, 38, 38, 35, 26, 16, 8, 3, 1, 1]),
&icdf::<17>(&[1, 1, 2, 6, 12, 21, 30, 36, 38, 36, 30, 21, 12, 6, 2, 1, 1]),
];
pub const SHELL_4: [&[u8]; 16] = [
&icdf::<2>(&[127, 129]),
&icdf::<3>(&[49, 157, 50]),
&icdf::<4>(&[20, 107, 109, 20]),
&icdf::<5>(&[11, 60, 113, 62, 10]),
&icdf::<6>(&[7, 36, 84, 87, 36, 6]),
&icdf::<7>(&[6, 24, 57, 82, 60, 23, 4]),
&icdf::<8>(&[5, 18, 39, 64, 68, 42, 16, 4]),
&icdf::<9>(&[6, 14, 29, 47, 61, 52, 30, 14, 3]),
&icdf::<10>(&[1, 15, 23, 35, 51, 50, 40, 30, 10, 1]),
&icdf::<11>(&[1, 1, 21, 32, 42, 52, 46, 41, 18, 1, 1]),
&icdf::<12>(&[1, 6, 16, 27, 36, 42, 42, 36, 27, 16, 6, 1]),
&icdf::<13>(&[1, 5, 12, 21, 31, 38, 40, 38, 31, 21, 12, 5, 1]),
&icdf::<14>(&[1, 3, 9, 17, 26, 34, 38, 38, 34, 26, 17, 9, 3, 1]),
&icdf::<15>(&[1, 3, 7, 14, 22, 29, 34, 36, 34, 29, 22, 14, 7, 3, 1]),
&icdf::<16>(&[1, 2, 5, 11, 18, 25, 31, 35, 35, 31, 25, 18, 11, 5, 2, 1]),
&icdf::<17>(&[1, 1, 4, 9, 15, 21, 28, 32, 34, 32, 28, 21, 15, 9, 4, 1, 1]),
];
pub const SHELL_2: [&[u8]; 16] = [
&icdf::<2>(&[128, 128]),
&icdf::<3>(&[42, 172, 42]),
&icdf::<4>(&[21, 107, 107, 21]),
&icdf::<5>(&[12, 60, 112, 61, 11]),
&icdf::<6>(&[8, 34, 86, 86, 35, 7]),
&icdf::<7>(&[8, 23, 55, 90, 55, 20, 5]),
&icdf::<8>(&[5, 15, 38, 72, 72, 36, 15, 3]),
&icdf::<9>(&[6, 12, 27, 52, 77, 47, 20, 10, 5]),
&icdf::<10>(&[6, 19, 28, 35, 40, 40, 35, 28, 19, 6]),
&icdf::<11>(&[4, 14, 22, 31, 37, 40, 37, 31, 22, 14, 4]),
&icdf::<12>(&[3, 10, 18, 26, 33, 38, 38, 33, 26, 18, 10, 3]),
&icdf::<13>(&[2, 8, 13, 21, 29, 36, 38, 36, 29, 21, 13, 8, 2]),
&icdf::<14>(&[1, 5, 10, 17, 25, 32, 38, 38, 32, 25, 17, 10, 5, 1]),
&icdf::<15>(&[1, 4, 7, 13, 21, 29, 35, 36, 35, 29, 21, 13, 7, 4, 1]),
&icdf::<16>(&[1, 2, 5, 10, 17, 25, 32, 36, 36, 32, 25, 17, 10, 5, 2, 1]),
&icdf::<17>(&[1, 2, 4, 7, 13, 21, 28, 34, 36, 34, 28, 21, 13, 7, 4, 2, 1]),
];
pub const EXC_LSB: [u8; 2] = icdf(&[136, 120]);
pub const EXC_SIGN: [&[u8]; 42] = [
&icdf::<2>(&[2, 254]),
&icdf::<2>(&[207, 49]),
&icdf::<2>(&[189, 67]),
&icdf::<2>(&[179, 77]),
&icdf::<2>(&[174, 82]),
&icdf::<2>(&[163, 93]),
&icdf::<2>(&[157, 99]),
&icdf::<2>(&[58, 198]),
&icdf::<2>(&[245, 11]),
&icdf::<2>(&[238, 18]),
&icdf::<2>(&[232, 24]),
&icdf::<2>(&[225, 31]),
&icdf::<2>(&[220, 36]),
&icdf::<2>(&[211, 45]),
&icdf::<2>(&[1, 255]),
&icdf::<2>(&[210, 46]),
&icdf::<2>(&[190, 66]),
&icdf::<2>(&[178, 78]),
&icdf::<2>(&[169, 87]),
&icdf::<2>(&[162, 94]),
&icdf::<2>(&[152, 104]),
&icdf::<2>(&[48, 208]),
&icdf::<2>(&[242, 14]),
&icdf::<2>(&[235, 21]),
&icdf::<2>(&[224, 32]),
&icdf::<2>(&[214, 42]),
&icdf::<2>(&[205, 51]),
&icdf::<2>(&[190, 66]),
&icdf::<2>(&[1, 255]),
&icdf::<2>(&[162, 94]),
&icdf::<2>(&[152, 104]),
&icdf::<2>(&[147, 109]),
&icdf::<2>(&[144, 112]),
&icdf::<2>(&[141, 115]),
&icdf::<2>(&[138, 118]),
&icdf::<2>(&[8, 248]),
&icdf::<2>(&[203, 53]),
&icdf::<2>(&[187, 69]),
&icdf::<2>(&[176, 80]),
&icdf::<2>(&[168, 88]),
&icdf::<2>(&[161, 95]),
&icdf::<2>(&[154, 102]),
];
pub const LBRR_FLAGS: [&[u8]; 2] = [&icdf::<3>(&[53, 53, 150]), &icdf::<7>(&[41, 20, 29, 41, 15, 28, 82])];
+123
View File
@@ -0,0 +1,123 @@
//! 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];
+41 -7
View File
@@ -1,20 +1,31 @@
//! Clean-room pure-Rust Opus (RFC 6716, updated by RFC 8251; RFC 7845 `OpusHead`).
//!
//! Decoding (feature `opus-decode`): [`Decoder`] (multistream, surround), [`StreamDecoder`] (one
//! elementary stream), [`OpusHead`] and the packet parser [`packet`]. From FilmCraft's
//! `filmcraft-opus` (storytold/filmcraft@8fda2f2).
//! - Decoding (feature `opus-decode`): [`Decoder`] (multistream, surround), [`StreamDecoder`] (one
//! elementary stream), [`OpusHead`] and the packet parser [`packet`]. From FilmCraft's
//! `filmcraft-opus` (storytold/filmcraft@8fda2f2).
//! - Encoding (feature `opus-encode`): [`Encoder`] (in [`encoder`]) (SILK, hybrid and CELT modes, 48 kHz,
//! 20 ms, mono or stereo). From EffectCraft's `effectcraft-opusenc` (storytold/effectcraft@c618429).
//!
//! Public items keep their names, so FilmCraft can re-export this module.
//! Public items keep their names, so FilmCraft and EffectCraft can re-export this module. The two
//! halves come from separate crates and do not share code yet (each has its own range coder, MDCT
//! and band tables). [`Bandwidth`] is one type for both (`encoder::Bandwidth` is the same type);
//! the two `Mode` enums keep their apps' variant names and convert into each other.
// Index loops and long argument lists follow the specification's pseudo-code (FilmCraft allows
// these lints workspace-wide).
// Index loops and long argument lists follow the specification's pseudo-code (FilmCraft and
// EffectCraft allow these lints workspace-wide).
#![allow(clippy::needless_range_loop, clippy::too_many_arguments)]
#[cfg(feature = "opus-decode")]
mod decode;
#[cfg(feature = "opus-encode")]
pub mod encoder;
#[cfg(feature = "opus-decode")]
pub use decode::{Decoder, Error, Mode, OpusHead, Packet, Result, SAMPLE_RATES, StreamDecoder, Toc, packet};
#[cfg(feature = "opus-encode")]
pub use encoder::Encoder;
/// Coded audio bandwidth (RFC 6716 §2, Table 1).
/// Coded audio bandwidth (RFC 6716 §2, Table 1), for the decoder and the encoder alike.
#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub enum Bandwidth {
/// 4 kHz (8 kHz sampling).
@@ -29,9 +40,32 @@ pub enum Bandwidth {
Fullband,
}
#[cfg(all(feature = "opus-decode", feature = "opus-encode"))]
impl From<Mode> for encoder::Mode {
fn from(m: Mode) -> Self {
match m {
Mode::SilkOnly => encoder::Mode::Silk,
Mode::Hybrid => encoder::Mode::Hybrid,
Mode::CeltOnly => encoder::Mode::Celt,
}
}
}
#[cfg(all(feature = "opus-decode", feature = "opus-encode"))]
impl From<encoder::Mode> for Mode {
fn from(m: encoder::Mode) -> Self {
match m {
encoder::Mode::Silk => Mode::SilkOnly,
encoder::Mode::Hybrid => Mode::Hybrid,
encoder::Mode::Celt => Mode::CeltOnly,
}
}
}
/// Channel order for Vorbis-ordered surround (mapping family 1, RFC 7845 §5.1.1.2): output channel
/// → decoded channel, giving the WAV/SMPTE order (L R C LFE back… side…) the rest of craft-audio
/// uses. From FilmCraft's `filmcraft-codecs`.
#[cfg(feature = "opus-decode")]
pub fn wav_order(channels: usize) -> Option<&'static [usize]> {
Some(match channels {
3 => &[0, 2, 1],
+3 -2
View File
@@ -1,5 +1,5 @@
//! Tests of the codecs ported from FilmCraft, built as one test binary: each file is a module, as
//! in their original crates.
//! Tests of the codecs ported from FilmCraft and EffectCraft, built as one test binary: each file
//! is a module, as in their original crates.
// Index loops follow the specifications, as in FilmCraft (which allows this lint workspace-wide).
#![allow(clippy::needless_range_loop)]
@@ -16,5 +16,6 @@ mod mp3_decoder_oracle;
mod mp3_encoder_oracle;
mod mpa_oracle;
mod opus;
mod opusenc;
mod routing_oracle;
mod track_decoder;
+286
View File
@@ -0,0 +1,286 @@
//! Shared helpers for the SILK/hybrid round-trip tests: synthetic signals, encoding, decoding
//! with FilmCraft's spec-derived decoder, delay-searched SNR and an Ogg Opus writer for ffmpeg.
#![allow(dead_code)]
use craft_audio::opus::encoder::OpusEncoder;
pub const SR: f64 = 48_000.0;
pub const FRAME: usize = OpusEncoder::FRAME_SIZE;
/// Two-pole resonator state.
struct Reson {
y1: f64,
y2: f64,
}
impl Reson {
fn run(&mut self, x: f64, f: f64, bw: f64) -> f64 {
let r = (-std::f64::consts::PI * bw / SR).exp();
let a1 = 2.0 * r * (2.0 * std::f64::consts::PI * f / SR).cos();
let a2 = -r * r;
let g = 1.0 - r;
let y = g * x + a1 * self.y1 + a2 * self.y2;
self.y2 = self.y1;
self.y1 = y;
y
}
}
/// Speech-like test signal: a glottal pulse train (100–200 Hz, gliding pitch) through three
/// formant resonators whose frequencies change per syllable, with syllable envelopes, short
/// pauses and a little aspiration noise.
pub fn speech(len: usize, seed: u32) -> Vec<f32> {
let vowels: [[f64; 3]; 5] =
[[730.0, 1090.0, 2440.0], [270.0, 2290.0, 3010.0], [530.0, 1840.0, 2480.0], [570.0, 840.0, 2410.0], [300.0, 870.0, 2240.0]];
let syl = (0.24 * SR) as usize;
let mut out = Vec::with_capacity(len);
let mut res = [Reson { y1: 0.0, y2: 0.0 }, Reson { y1: 0.0, y2: 0.0 }, Reson { y1: 0.0, y2: 0.0 }];
let mut phase = 0.0f64;
let mut s = seed;
let mut prev_glottal = 0.0;
for n in 0..len {
let t = n as f64 / SR;
let k = n / syl;
let pos = (n % syl) as f64 / syl as f64;
// Every fourth syllable is a pause.
let pause = k % 4 == 3;
let env = if pause { 0.0 } else { (std::f64::consts::PI * pos).sin().powf(0.7) };
let f0 = 100.0 + 50.0 * (1.0 + (2.0 * std::f64::consts::PI * 0.9 * t).sin()) + 15.0 * ((k * 7) % 5) as f64 - 30.0 * pos;
phase += f0 / SR;
if phase >= 1.0 {
phase -= 1.0;
}
// Rosenberg-like glottal flow; its derivative excites the vocal tract.
let g = if phase < 0.4 {
0.5 * (1.0 - (std::f64::consts::PI * phase / 0.4).cos())
} else if phase < 0.56 {
(std::f64::consts::FRAC_PI_2 * (phase - 0.4) / 0.16).cos()
} else {
0.0
};
let dg = (g - prev_glottal) * 40.0;
prev_glottal = g;
s = s.wrapping_mul(1_664_525).wrapping_add(1_013_904_223);
let noise = ((s >> 8) as f64 / (1u32 << 24) as f64 - 0.5) * 0.02;
let v = vowels[(k * 3 + 1) % 5];
let w = vowels[(k * 3 + 4) % 5];
let mix = pos;
let mut x = dg + noise;
let mut y = 0.0;
for (i, r) in res.iter_mut().enumerate() {
let f = v[i] * (1.0 - mix) + w[i] * mix;
y += r.run(x, f, 60.0 + 40.0 * i as f64) * [1.0, 0.6, 0.3][i];
}
x = y * env * 1.2;
out.push(x as f32);
}
let peak = out.iter().fold(0f32, |m, v| m.max(v.abs()));
let g = 0.5 / peak.max(1e-9);
out.iter().map(|v| v * g).collect()
}
/// Music-like test signal: three-note chords of harmonic tones (six partials each), changing
/// every half second with soft attacks and releases.
pub fn music(len: usize) -> Vec<f32> {
let chords: [[f64; 3]; 4] = [[261.63, 329.63, 392.0], [220.0, 277.18, 329.63], [174.61, 220.0, 261.63], [196.0, 246.94, 293.66]];
let seg = (0.5 * SR) as usize;
(0..len)
.map(|n| {
let t = n as f64 / SR;
let k = n / seg;
let pos = (n % seg) as f64 / seg as f64;
let env = (pos * 20.0).min(1.0) * ((1.0 - pos) * 10.0).min(1.0);
let mut x = 0.0;
for (j, &f) in chords[k % 4].iter().enumerate() {
for h in 1..=6 {
x += (2.0 * std::f64::consts::PI * f * h as f64 * t + j as f64).sin() / h as f64;
}
}
(0.12 * env * x) as f32
})
.collect()
}
pub fn interleave(planar: &[Vec<f32>]) -> 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;
}
}
pcm
}
/// Encodes planar input (padding the tail so the encoder delay is flushed).
pub fn encode(mut enc: OpusEncoder, planar: &[Vec<f32>]) -> (OpusEncoder, Vec<Vec<u8>>) {
let c = planar.len();
let len = planar[0].len();
let total = len + enc.pre_skip() as usize + 16;
let frames = total.div_ceil(FRAME);
let mut pcm = interleave(planar);
pcm.resize(frames * FRAME * c, 0.0);
let packets = pcm.chunks(FRAME * c).map(|f| enc.encode_float(f).unwrap()).collect();
(enc, packets)
}
/// Decodes with filmcraft-opus, drops the pre-skip and returns planar channels (full length).
pub fn decode(head: &[u8], packets: &[Vec<u8>]) -> Vec<Vec<f32>> {
let mut dec = craft_audio::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);
}
}
out
}
/// Waveform SNR (dB) of `dec[i + shift]` against `orig[i]`, excluding the first and last 20 ms.
pub fn snr_at(orig: &[f32], dec: &[f32], shift: isize) -> f64 {
let (mut s, mut e) = (0f64, 0f64);
for i in FRAME..orig.len() - FRAME {
let j = i as isize + shift;
let d = if j >= 0 && (j as usize) < dec.len() { dec[j as usize] as f64 } else { 0.0 };
s += (orig[i] as f64).powi(2);
e += (orig[i] as f64 - d).powi(2);
}
10.0 * (s / e.max(1e-30)).log10()
}
/// Best SNR over residual delays of ±`range` samples: (SNR dB, delay).
pub fn snr_search(orig: &[f32], dec: &[f32], range: isize) -> (f64, isize) {
(-range..=range).map(|d| (snr_at(orig, dec, d), d)).fold((f64::MIN, 0), |a, b| if b.0 > a.0 { b } else { a })
}
/// Segmental SNR (dB): mean over 20 ms segments holding signal (each clamped to [-10, 40]).
pub fn seg_snr(orig: &[f32], dec: &[f32], shift: isize) -> f64 {
let mut acc = 0f64;
let mut n = 0;
let mut i = FRAME;
while i + FRAME <= orig.len() - FRAME {
let (mut s, mut e) = (0f64, 0f64);
for k in i..i + FRAME {
let j = k as isize + shift;
let d = if j >= 0 && (j as usize) < dec.len() { dec[j as usize] as f64 } else { 0.0 };
s += (orig[k] as f64).powi(2);
e += (orig[k] as f64 - d).powi(2);
}
if s / FRAME as f64 > 1e-5 {
acc += (10.0 * (s / e.max(1e-30)).log10()).clamp(-10.0, 40.0);
n += 1;
}
i += FRAME;
}
acc / n.max(1) as f64
}
/// Average bitrate (bits per second) of a packet sequence of 20 ms packets.
pub fn bitrate(packets: &[Vec<u8>]) -> f64 {
packets.iter().map(|p| p.len() * 8).sum::<usize>() as f64 / (packets.len() as f64 * 0.02)
}
/// Low-pass filtered copy (FFT brick wall) used to compare band-limited modes fairly.
pub fn lowpass(x: &[f32], fc: f64) -> Vec<f32> {
let n = x.len().next_power_of_two();
let mut re: Vec<f64> = x.iter().map(|&v| v as f64).collect();
re.resize(n, 0.0);
let mut im = vec![0f64; n];
fft(&mut re, &mut im, false);
for k in 0..n {
let f = (k.min(n - k)) as f64 * SR / n as f64;
if f > fc {
re[k] = 0.0;
im[k] = 0.0;
}
}
fft(&mut re, &mut im, true);
re[..x.len()].iter().map(|&v| (v / n as f64) as f32).collect()
}
pub fn fft(re: &mut [f64], im: &mut [f64], inverse: bool) {
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 = if inverse { 2.0 } else { -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;
}
}
// ---- 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());
out.extend_from_slice(&seq.to_le_bytes());
out.extend_from_slice(&0u32.to_le_bytes());
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());
}
/// An Ogg Opus file holding `packets` for `len` samples of audio.
pub 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
}
@@ -0,0 +1,81 @@
//! ffmpeg as a second, external oracle: every mode/bandwidth written as Ogg Opus must decode
//! without errors with both ffmpeg's native `opus` decoder and `libopus`, honour the pre-skip
//! and end granule, and reproduce the input. Skipped when ffmpeg is not installed.
use super::common;
use std::process::Command;
use common::*;
use craft_audio::opus::encoder::{Bandwidth, Mode, OpusEncoder};
/// ffmpeg, found and skipped by craft-media's oracle rules (EffectCraft looked in Homebrew's
/// directory, then `PATH`).
fn ffmpeg() -> Option<std::path::PathBuf> {
crate::testkit::ffmpeg_or_skip("opus encoder ffmpeg oracle")
}
/// Decodes an Ogg Opus file with ffmpeg's decoder `codec`; returns interleaved f32 output.
fn ffmpeg_decode(exe: &std::path::Path, codec: &str, file: &[u8], channels: usize, tag: &str) -> Option<Vec<f32>> {
let has = Command::new(exe).args(["-hide_banner", "-decoders"]).output().ok()?;
if !String::from_utf8_lossy(&has.stdout).lines().any(|l| l.split_whitespace().nth(1) == Some(codec)) {
eprintln!("ffmpeg has no {codec} decoder; skipping it");
return None;
}
let path = std::env::temp_dir().join(format!("effectcraft-opusenc-{}-{tag}.opus", std::process::id()));
std::fs::write(&path, file).unwrap();
let out = Command::new(exe)
.args(["-v", "error", "-c:a", codec, "-i"])
.arg(&path)
.args(["-f", "f32le", "-ac", &channels.to_string(), "-ar", "48000", "-"])
.output()
.expect("run ffmpeg");
let _ = std::fs::remove_file(&path);
let stderr = String::from_utf8_lossy(&out.stderr);
assert!(out.status.success(), "ffmpeg ({codec}) failed: {stderr}");
// ffmpeg's native decoder also makes its raw-PCM muxer print a timestamp complaint for
// narrowband SILK streams written by libopus itself; that line is about ffmpeg's own
// timestamping, not the bitstream, so it is the one message tolerated.
let errors: Vec<&str> = stderr.lines().filter(|l| !l.trim().is_empty() && !l.contains("non monotonically increasing dts to muxer")).collect();
assert!(errors.is_empty(), "ffmpeg ({codec}) reported errors: {errors:?}");
Some(out.stdout.as_chunks::<4>().0.iter().map(|b| f32::from_le_bytes([b[0], b[1], b[2], b[3]])).collect())
}
#[test]
fn ffmpeg_decodes_every_mode() {
let Some(exe) = ffmpeg() else {
eprintln!("ffmpeg not installed; skipping the oracle test");
return;
};
let len = 48_000;
let sp = speech(len + 30_000, 3);
let left = sp[..len].to_vec();
let right: Vec<f32> = sp[30_000..].iter().map(|v| v * 0.7).collect();
let cases = [
(1usize, Mode::Silk, Bandwidth::Narrowband, 8_000u32, 3_400.0, 8.0),
(1, Mode::Silk, Bandwidth::Mediumband, 11_000, 5_000.0, 10.0),
(1, Mode::Silk, Bandwidth::Wideband, 16_000, 7_000.0, 12.0),
(2, Mode::Silk, Bandwidth::Wideband, 32_000, 7_000.0, 12.0),
(1, Mode::Hybrid, Bandwidth::SuperWideband, 24_000, 24_000.0, 10.0),
(1, Mode::Hybrid, Bandwidth::Fullband, 32_000, 24_000.0, 10.0),
(2, Mode::Hybrid, Bandwidth::Fullband, 64_000, 24_000.0, 10.0),
(2, Mode::Celt, Bandwidth::Fullband, 96_000, 24_000.0, 10.0),
];
for (c, mode, bw, rate, fc, min_snr) in cases {
let planar: Vec<Vec<f32>> = if c == 1 { vec![left.clone()] } else { vec![left.clone(), right.clone()] };
let (enc, packets) = encode(OpusEncoder::with_mode(c as u8, rate, mode, bw).unwrap(), &planar);
let file = ogg_opus(&enc, &packets, len);
for codec in ["opus", "libopus"] {
let Some(dec) = ffmpeg_decode(&exe, codec, &file, c, &format!("{mode:?}{bw:?}{c}{codec}")) else {
continue;
};
assert_eq!(dec.len(), len * c, "{codec}: pre-skip and end granule must be honoured");
for (ch, x) in planar.iter().enumerate() {
let d: Vec<f32> = (0..len).map(|i| dec[i * c + ch]).collect();
let (snr, delay) = snr_search(&lowpass(x, fc), &d, 6);
eprintln!("ffmpeg {codec:7} {mode:?}/{bw:?} {c} ch @ {} kb/s, channel {ch}: SNR {snr:.1} dB (residual delay {delay})", rate / 1000);
assert!(snr > min_snr, "{codec} {mode:?}/{bw:?}: SNR {snr:.1} dB");
}
}
}
}
+5
View File
@@ -0,0 +1,5 @@
//! EffectCraft's `effectcraft-opusenc` integration tests.
mod common;
mod ffmpeg;
mod modes;
+286
View File
@@ -0,0 +1,286 @@
//! SILK, hybrid and CELT round trips through FilmCraft's spec-derived Opus decoder: every
//! mode/bandwidth decodes, waveform quality on speech-like and music-like signals, bitrate
//! tracking, stereo, the application tables and delay alignment.
use super::common;
use common::*;
use craft_audio::opus::encoder::{Application, Bandwidth, Error, Mode, OpusEncoder};
struct Run {
snr: f64,
seg: f64,
delay: isize,
kbps: f64,
packets: Vec<Vec<u8>>,
}
/// Encodes `planar`, decodes, and measures each channel against `reference` (searching a
/// residual delay of up to ±6 samples).
fn run(enc: OpusEncoder, planar: &[Vec<f32>], reference: &[Vec<f32>]) -> Vec<Run> {
let (enc, packets) = encode(enc, planar);
let dec = decode(&enc.opus_head(), &packets);
let kbps = bitrate(&packets) / 1000.0;
reference
.iter()
.zip(&dec)
.map(|(r, d)| {
assert!(d.len() >= r.len(), "decoded {} < {}", d.len(), r.len());
assert!(d.iter().all(|v| v.is_finite()));
let (snr, delay) = snr_search(r, d, 6);
Run { snr, seg: seg_snr(r, d, delay), delay, kbps, packets: packets.clone() }
})
.collect()
}
fn mono(mode: Mode, bw: Bandwidth, rate: u32, x: &[f32], fc: f64) -> Run {
let reference = if fc < 20_000.0 { lowpass(x, fc) } else { x.to_vec() };
run(OpusEncoder::with_mode(1, rate, mode, bw).unwrap(), &[x.to_vec()], &[reference]).remove(0)
}
/// Upper edge of the band a SILK bandwidth keeps (for band-limited references).
fn band_edge(bw: Bandwidth) -> f64 {
match bw {
Bandwidth::Narrowband => 3_400.0,
Bandwidth::Mediumband => 5_000.0,
Bandwidth::Wideband => 7_000.0,
_ => 24_000.0,
}
}
const LEN: usize = 96_000;
#[test]
fn every_mode_and_bandwidth_decodes() {
let sp = speech(LEN + 40_000, 11);
let left = sp[..LEN].to_vec();
let right: Vec<f32> = sp[40_000..].iter().map(|v| v * 0.6).collect();
for (mode, bw, config, rate) in [
(Mode::Silk, Bandwidth::Narrowband, 1u8, 8_000u32),
(Mode::Silk, Bandwidth::Mediumband, 5, 11_000),
(Mode::Silk, Bandwidth::Wideband, 9, 16_000),
(Mode::Hybrid, Bandwidth::SuperWideband, 13, 24_000),
(Mode::Hybrid, Bandwidth::Fullband, 15, 32_000),
(Mode::Celt, Bandwidth::Fullband, 31, 64_000),
] {
for c in [1usize, 2] {
let planar = if c == 1 { vec![left.clone()] } else { vec![left.clone(), right.clone()] };
let reference: Vec<Vec<f32>> = planar.iter().map(|x| lowpass(x, band_edge(bw))).collect();
let enc = OpusEncoder::with_mode(c as u8, rate * c as u32, mode, bw).unwrap();
assert_eq!((enc.mode(), enc.bandwidth()), (mode, bw));
let res = run(enc, &planar, &reference);
for p in &res[0].packets {
assert_eq!(p[0], (config << 3) | if c == 2 { 4 } else { 0 }, "TOC of {mode:?}/{bw:?}");
assert!(p.len() >= 3 && p.len() <= 1276);
}
for (ch, r) in res.iter().enumerate() {
eprintln!(
"{mode:?}/{bw:?} {c} ch @ {} kb/s: {:.1} kb/s, channel {ch} SNR {:.1} dB, segSNR {:.1} dB, residual delay {}",
rate * c as u32 / 1000,
r.kbps,
r.snr,
r.seg,
r.delay
);
assert!(r.snr > 8.0, "{mode:?}/{bw:?} {c} ch channel {ch}: SNR {:.1} dB", r.snr);
assert!(r.delay.abs() <= 1, "{mode:?}/{bw:?}: residual delay {} (pre-skip mismatch)", r.delay);
}
}
}
}
#[test]
fn silk_speech_quality() {
let x = speech(LEN, 7);
for (bw, rate, min_snr) in [
(Bandwidth::Narrowband, 8_000u32, 10.0),
(Bandwidth::Mediumband, 11_000, 14.0),
(Bandwidth::Wideband, 12_000, 14.0),
(Bandwidth::Wideband, 16_000, 18.0),
(Bandwidth::Wideband, 24_000, 24.0),
] {
let r = mono(Mode::Silk, bw, rate, &x, band_edge(bw));
eprintln!("SILK {bw:?} speech @ {} kb/s: {:.1} kb/s, SNR {:.1} dB, segSNR {:.1} dB", rate / 1000, r.kbps, r.snr, r.seg);
assert!(r.snr > min_snr && r.seg > min_snr - 3.0, "SILK {bw:?} @ {rate}: SNR {:.1} / seg {:.1}", r.snr, r.seg);
}
}
#[test]
fn silk_beats_celt_on_speech_at_low_rate() {
let x = speech(LEN, 5);
// Both measured against the full-band input.
let silk = mono(Mode::Silk, Bandwidth::Wideband, 12_000, &x, 24_000.0);
let celt = mono(Mode::Celt, Bandwidth::Fullband, 12_000, &x, 24_000.0);
eprintln!(
"12 kb/s mono speech: SILK WB {:.1} dB (seg {:.1}) at {:.1} kb/s vs CELT {:.1} dB (seg {:.1}) at {:.1} kb/s",
silk.snr, silk.seg, silk.kbps, celt.snr, celt.seg, celt.kbps
);
assert!(silk.snr > celt.snr + 3.0, "SILK {:.1} dB vs CELT {:.1} dB", silk.snr, celt.snr);
assert!(silk.seg > celt.seg + 3.0, "SILK seg {:.1} dB vs CELT seg {:.1} dB", silk.seg, celt.seg);
}
#[test]
fn hybrid_quality() {
let sp = speech(LEN, 9);
let mu = music(LEN);
for (name, x, bw, rate, min_snr) in [
("speech", &sp, Bandwidth::SuperWideband, 24_000u32, 18.0),
("speech", &sp, Bandwidth::Fullband, 32_000, 20.0),
("music", &mu, Bandwidth::SuperWideband, 24_000, 9.0),
("music", &mu, Bandwidth::Fullband, 32_000, 10.0),
] {
let r = mono(Mode::Hybrid, bw, rate, x, 24_000.0);
eprintln!("hybrid {bw:?} {name} @ {} kb/s: {:.1} kb/s, SNR {:.1} dB, segSNR {:.1} dB", rate / 1000, r.kbps, r.snr, r.seg);
assert!(r.snr > min_snr, "hybrid {bw:?} {name}: SNR {:.1} dB", r.snr);
// Constant packet size.
assert!(r.packets.iter().all(|p| p.len() == rate as usize / 400));
}
}
#[test]
fn silk_music_quality() {
let x = music(LEN);
for (rate, min_snr) in [(12_000u32, 5.0), (16_000, 8.0), (24_000, 13.0)] {
let r = mono(Mode::Silk, Bandwidth::Wideband, rate, &x, 7_000.0);
eprintln!("SILK WB music @ {} kb/s: {:.1} kb/s, SNR {:.1} dB", rate / 1000, r.kbps, r.snr);
assert!(r.snr > min_snr, "SILK music @ {rate}: SNR {:.1} dB", r.snr);
}
}
#[test]
fn bitrate_tracking() {
let sp = speech(LEN, 13);
let mu = music(LEN);
for (name, x) in [("speech", &sp), ("music", &mu)] {
for (mode, bw, rate) in [
(Mode::Silk, Bandwidth::Narrowband, 6_000u32),
(Mode::Silk, Bandwidth::Narrowband, 9_000),
(Mode::Silk, Bandwidth::Mediumband, 11_000),
(Mode::Silk, Bandwidth::Wideband, 12_000),
(Mode::Silk, Bandwidth::Wideband, 20_000),
(Mode::Silk, Bandwidth::Wideband, 32_000),
(Mode::Hybrid, Bandwidth::SuperWideband, 20_000),
(Mode::Hybrid, Bandwidth::Fullband, 36_000),
] {
let r = mono(mode, bw, rate, x, band_edge(bw));
let ratio = r.kbps * 1000.0 / rate as f64;
eprintln!(
"{name} {mode:?}/{bw:?} target {} kb/s: actual {:.2} kb/s ({:+.0} %), SNR {:.1} dB",
rate / 1000,
r.kbps,
(ratio - 1.0) * 100.0,
r.snr
);
assert!((0.75..=1.25).contains(&ratio), "{name} {mode:?}/{bw:?} @ {rate}: {:.2} kb/s", r.kbps);
}
}
}
#[test]
fn stereo_silk_and_hybrid() {
let sp = speech(LEN + 40_000, 21);
let a = sp[..LEN].to_vec();
let b: Vec<f32> = sp[40_000..].iter().map(|v| v * 0.7).collect();
let mu = music(LEN);
let cases: Vec<(&str, Vec<f32>, Vec<f32>)> = vec![
("dual mono", a.clone(), a.clone()),
("panned", a.iter().map(|v| v * 0.9).collect(), a.iter().map(|v| v * 0.35).collect()),
("independent", a.clone(), b),
("speech+music", a.iter().zip(&mu).map(|(x, y)| 0.8 * x + 0.3 * y).collect(), a.iter().zip(&mu).map(|(x, y)| 0.3 * x + 0.8 * y).collect()),
];
for (name, l, r) in &cases {
for (mode, bw, rate, min_snr) in [(Mode::Silk, Bandwidth::Wideband, 32_000u32, 9.0), (Mode::Hybrid, Bandwidth::Fullband, 64_000, 9.0)] {
let fc = band_edge(bw);
let res = run(OpusEncoder::with_mode(2, rate, mode, bw).unwrap(), &[l.clone(), r.clone()], &[lowpass(l, fc), lowpass(r, fc)]);
eprintln!(
"stereo {name} {mode:?}/{bw:?} @ {} kb/s: {:.1} kb/s, L {:.1} dB, R {:.1} dB",
rate / 1000,
res[0].kbps,
res[0].snr,
res[1].snr
);
assert!(res[0].snr > min_snr && res[1].snr > min_snr, "stereo {name} {mode:?}: L {:.1} R {:.1}", res[0].snr, res[1].snr);
}
}
}
#[test]
fn application_tables_and_pre_skip() {
let cases = [
(Application::Voip, 1u8, 8_000u32, Mode::Silk, Bandwidth::Narrowband),
(Application::Voip, 1, 11_000, Mode::Silk, Bandwidth::Mediumband),
(Application::Voip, 1, 16_000, Mode::Silk, Bandwidth::Wideband),
(Application::Voip, 1, 24_000, Mode::Hybrid, Bandwidth::SuperWideband),
(Application::Voip, 1, 32_000, Mode::Hybrid, Bandwidth::Fullband),
(Application::Voip, 1, 48_000, Mode::Celt, Bandwidth::Fullband),
(Application::Voip, 2, 32_000, Mode::Silk, Bandwidth::Wideband),
(Application::Audio, 1, 12_000, Mode::Silk, Bandwidth::Wideband),
(Application::Audio, 1, 24_000, Mode::Hybrid, Bandwidth::Fullband),
(Application::Audio, 2, 48_000, Mode::Hybrid, Bandwidth::Fullband),
(Application::Audio, 1, 64_000, Mode::Celt, Bandwidth::Fullband),
(Application::Audio, 2, 128_000, Mode::Celt, Bandwidth::Fullband),
];
for (app, c, rate, mode, bw) in cases {
let e = OpusEncoder::with_application(c, rate, app).unwrap();
assert_eq!((e.mode(), e.bandwidth(), e.application()), (mode, bw, app), "{app:?} {c} ch @ {rate}");
assert_eq!(e.pre_skip(), 120, "every mode is aligned to a 120-sample delay");
}
assert_eq!(OpusEncoder::new(2, 128_000).unwrap().mode(), Mode::Celt);
assert_eq!(OpusEncoder::new(1, 12_000).unwrap().application(), Application::Audio);
}
/// EffectCraft's encoder panicked here; craft-audio returns an error (never crash).
#[test]
fn invalid_mode_bandwidth_is_an_error() {
assert_eq!(OpusEncoder::with_mode(1, 16_000, Mode::Silk, Bandwidth::Fullband).err(), Some(Error::ModeBandwidth(Mode::Silk, Bandwidth::Fullband)));
assert_eq!(OpusEncoder::new(3, 128_000).err(), Some(Error::Channels(3)));
assert_eq!(OpusEncoder::new(0, 128_000).err(), Some(Error::Channels(0)));
}
#[test]
fn extreme_inputs_stay_decodable() {
let len = 24_000;
let mut s = 1u32;
let noise: Vec<f32> = (0..len)
.map(|_| {
s = s.wrapping_mul(1_664_525).wrapping_add(1_013_904_223);
(s >> 8) as f32 / (1u32 << 23) as f32 - 1.0
})
.collect();
let square: Vec<f32> = (0..len).map(|i| if (i / 60) % 2 == 0 { 0.999 } else { -0.999 }).collect();
let clicks: Vec<f32> = (0..len).map(|i| if i % 4_801 == 0 { 1.0 } else { 0.0 }).collect();
let dc = vec![0.7f32; len];
let nan: Vec<f32> = (0..len).map(|i| if i % 1000 == 0 { f32::NAN } else { 0.1 * (i as f32 * 0.01).sin() }).collect();
for (name, x) in [("noise", &noise), ("square", &square), ("clicks", &clicks), ("dc", &dc), ("nan", &nan)] {
for (mode, bw, rate) in
[(Mode::Silk, Bandwidth::Narrowband, 6_000u32), (Mode::Silk, Bandwidth::Wideband, 40_000), (Mode::Hybrid, Bandwidth::Fullband, 24_000)]
{
for c in [1usize, 2] {
let planar: Vec<Vec<f32>> = if c == 1 { vec![x.clone()] } else { vec![x.clone(), x.iter().map(|v| -0.5 * v).collect()] };
let (enc, packets) = encode(OpusEncoder::with_mode(c as u8, rate * c as u32, mode, bw).unwrap(), &planar);
assert!(packets.iter().all(|p| p.len() >= 3 && p.len() <= 1276), "{name} {mode:?}: packet sizes");
let dec = decode(&enc.opus_head(), &packets);
let peak = dec.iter().flatten().fold(0f32, |m, v| m.max(v.abs()));
let snr = if name == "nan" { f64::NAN } else { snr_search(&planar[0], &dec[0], 2).0 };
eprintln!("{name} {mode:?}/{bw:?} {c} ch: peak {peak:.2}, SNR {snr:.1} dB, {:.1} kb/s", bitrate(&packets) / 1000.0);
assert!(dec.iter().all(|ch| ch.iter().all(|v| v.is_finite() && v.abs() <= 4.0)), "{name} {mode:?}/{bw:?} {c} ch: bad output");
}
}
}
}
#[test]
fn silk_and_hybrid_silence() {
for (mode, bw, rate) in [(Mode::Silk, Bandwidth::Wideband, 16_000u32), (Mode::Hybrid, Bandwidth::Fullband, 32_000)] {
let x = vec![0f32; 48_000];
let (enc, packets) = encode(OpusEncoder::with_mode(1, rate, mode, bw).unwrap(), std::slice::from_ref(&x));
let dec = decode(&enc.opus_head(), &packets);
let peak = dec[0].iter().fold(0f32, |m, v| m.max(v.abs()));
let kbps = bitrate(&packets) / 1000.0;
eprintln!("{mode:?} silence: {kbps:.2} kb/s, decoded peak {peak:.2e}");
assert!(peak < 1e-3, "{mode:?} silence decodes to peak {peak}");
if mode == Mode::Silk {
assert!(kbps < 4.0, "SILK silence costs {kbps:.2} kb/s");
}
}
}
+19 -1
View File
@@ -1,7 +1,8 @@
//! `TrackDecoder` on streams an app's own demuxer hands over without codec configuration:
//! ADTS-framed AAC (MPEG-TS).
//! ADTS-framed AAC (MPEG-TS) and Opus without an `OpusHead` (RTP).
use craft_audio::aac::{Encoder as AacEncoder, EncoderConfig};
use craft_audio::opus::Encoder as OpusEncoder;
use craft_audio::{CodecId, CodecParams, TrackDecoder};
fn tone(n: usize) -> Vec<f32> {
@@ -32,3 +33,20 @@ fn adts_aac_without_extradata() {
let raw = &units[0][7..];
assert!(TrackDecoder::new(&CodecParams::new(CodecId::Aac), &[raw]).is_err());
}
#[test]
fn opus_without_opus_head() {
let mut enc = OpusEncoder::new(1, 64_000).unwrap();
let pcm = tone(OpusEncoder::FRAME_SIZE * 10);
let packets: Vec<Vec<u8>> = pcm.chunks(OpusEncoder::FRAME_SIZE).map(|f| enc.encode_float(f).unwrap()).collect();
let params = CodecParams::new(CodecId::Opus).with_channels(1).with_sample_rate(48_000);
let mut dec = TrackDecoder::new(&params, &[]).unwrap();
assert_eq!((dec.sample_rate(), dec.channels()), (48_000, 1));
for p in &packets {
assert_eq!(dec.decode(p).unwrap().frames(), OpusEncoder::FRAME_SIZE);
}
// More than two channels need the mapping an OpusHead carries.
assert!(TrackDecoder::new(&CodecParams::new(CodecId::Opus).with_channels(6), &[]).is_err());
}