From bdf8cb2a19e9fd0994fa3e8bc382cee0a82ccb7e Mon Sep 17 00:00:00 2001 From: alehel Date: Fri, 30 Oct 2026 12:31:22 +0100 Subject: [PATCH] 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 --- .github/workflows/audio.yml | 2 +- crates/audio/Cargo.toml | 7 +- crates/audio/README.md | 51 +- crates/audio/src/lib.rs | 15 +- crates/audio/src/opus/encoder/bands.rs | 497 ++++++++++++ crates/audio/src/opus/encoder/energy.rs | 230 ++++++ crates/audio/src/opus/encoder/mdct.rs | 294 ++++++++ crates/audio/src/opus/encoder/mod.rs | 589 +++++++++++++++ crates/audio/src/opus/encoder/range.rs | 472 ++++++++++++ crates/audio/src/opus/encoder/rate.rs | 561 ++++++++++++++ crates/audio/src/opus/encoder/resample.rs | 131 ++++ .../audio/src/opus/encoder/silk/analysis.rs | 322 ++++++++ crates/audio/src/opus/encoder/silk/lpc.rs | 370 +++++++++ crates/audio/src/opus/encoder/silk/mod.rs | 714 ++++++++++++++++++ crates/audio/src/opus/encoder/silk/nlsf.rs | 131 ++++ crates/audio/src/opus/encoder/silk/pulses.rs | 113 +++ crates/audio/src/opus/encoder/silk/tables.rs | 543 +++++++++++++ crates/audio/src/opus/encoder/tables.rs | 123 +++ crates/audio/src/opus/mod.rs | 48 +- crates/audio/tests/codecs/main.rs | 5 +- crates/audio/tests/codecs/opusenc/common.rs | 286 +++++++ crates/audio/tests/codecs/opusenc/ffmpeg.rs | 81 ++ crates/audio/tests/codecs/opusenc/mod.rs | 5 + crates/audio/tests/codecs/opusenc/modes.rs | 286 +++++++ crates/audio/tests/codecs/track_decoder.rs | 20 +- 25 files changed, 5871 insertions(+), 25 deletions(-) create mode 100644 crates/audio/src/opus/encoder/bands.rs create mode 100644 crates/audio/src/opus/encoder/energy.rs create mode 100644 crates/audio/src/opus/encoder/mdct.rs create mode 100644 crates/audio/src/opus/encoder/mod.rs create mode 100644 crates/audio/src/opus/encoder/range.rs create mode 100644 crates/audio/src/opus/encoder/rate.rs create mode 100644 crates/audio/src/opus/encoder/resample.rs create mode 100644 crates/audio/src/opus/encoder/silk/analysis.rs create mode 100644 crates/audio/src/opus/encoder/silk/lpc.rs create mode 100644 crates/audio/src/opus/encoder/silk/mod.rs create mode 100644 crates/audio/src/opus/encoder/silk/nlsf.rs create mode 100644 crates/audio/src/opus/encoder/silk/pulses.rs create mode 100644 crates/audio/src/opus/encoder/silk/tables.rs create mode 100644 crates/audio/src/opus/encoder/tables.rs create mode 100644 crates/audio/tests/codecs/opusenc/common.rs create mode 100644 crates/audio/tests/codecs/opusenc/ffmpeg.rs create mode 100644 crates/audio/tests/codecs/opusenc/mod.rs create mode 100644 crates/audio/tests/codecs/opusenc/modes.rs diff --git a/.github/workflows/audio.yml b/.github/workflows/audio.yml index 0104f848e..0517efe6f 100644 --- a/.github/workflows/audio.yml +++ b/.github/workflows/audio.yml @@ -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::" diff --git a/crates/audio/Cargo.toml b/crates/audio/Cargo.toml index 822c8a8ad..9b6200198 100644 --- a/crates/audio/Cargo.toml +++ b/crates/audio/Cargo.toml @@ -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`). diff --git a/crates/audio/README.md b/crates/audio/README.md index abd3228f4..57f17054a 100644 --- a/crates/audio/README.md +++ b/crates/audio/README.md @@ -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 diff --git a/crates/audio/src/lib.rs b/crates/audio/src/lib.rs index 3ed684245..f8441a51d 100644 --- a/crates/audio/src/lib.rs +++ b/crates/audio/src/lib.rs @@ -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 for AudioError { } } +#[cfg(feature = "opus-encode")] +impl From 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] diff --git a/crates/audio/src/opus/encoder/bands.rs b/crates/audio/src/opus/encoder/bands.rs new file mode 100644 index 000000000..604aab595 --- /dev/null +++ b/crates/audio/src/opus/encoder/bands.rs @@ -0,0 +1,497 @@ +//! Band shape coding: PVQ search and codeword indexing, spreading, recursive band splitting +//! (RFC 6716 §4.3.4, encoder direction). +//! +//! The encoder codes 20 ms long-block frames only (no transients, no TF changes), and stereo +//! as dual (independent L/R) stereo, so a band is always coded by `quant_partition` with a single +//! block. Every bit-budget computation mirrors the decoder so both sides take the same +//! split/pulse decisions. + +use crate::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 = x.iter().map(|v| v.abs()).collect(); + let sum: f32 = abs.iter().sum(); + if !sum.is_finite() || sum <= 1e-15 { + y[0] = k as i32; + return; + } + let mut placed = 0usize; + let mut xy = 0f32; + let mut yy = 0f32; + if k > n / 2 { + // Pre-project onto the pyramid (rounding down) and place the rest greedily. + let r = (k as f32 - 1.0) / sum; + for j in 0..n { + let p = (abs[j] * r).floor() as i32; + y[j] = p; + placed += p as usize; + xy += abs[j] * p as f32; + yy += (p * p) as f32; + } + } + while placed < k { + let mut best = 0usize; + let mut best_num = -1f32; + let mut best_den = 1f32; + for j in 0..n { + let num = xy + abs[j]; + let num = num * num; + let den = yy + 2.0 * y[j] as f32 + 1.0; + if num * best_den > best_num * den { + best_num = num; + best_den = den; + best = j; + } + } + xy += abs[best]; + yy += 2.0 * y[best] as f32 + 1.0; + y[best] += 1; + placed += 1; + } + for j in 0..n { + if x[j] < 0.0 { + y[j] = -y[j]; + } + } +} + +/// Encodes the PVQ codeword `y` (`sum |y| = k`) as the index the RFC 6716 §4.3.4.2 decoder maps +/// back to `y`. +pub fn encode_pulses(y: &[i32], k: usize, enc: &mut RangeEncoder, v: &mut Vec) { + let n = y.len(); + let w = k + 1; + v.clear(); + v.resize((n + 1) * w, 0); + v[0] = 1; + for nn in 1..=n { + v[nn * w] = 1; + for kk in 1..=k { + let a = v[(nn - 1) * w + kk]; + let b = v[nn * w + kk - 1]; + let c = v[(nn - 1) * w + kk - 1]; + v[nn * w + kk] = a.saturating_add(b).saturating_add(c); + } + } + let total = v[n * w + k]; + let vv = |nn: usize, kx: usize| v[nn * w + kx]; + let mut index = 0u64; + let mut kk = k; + for j in 0..n { + if kk == 0 { + break; + } + let rem = n - j; + let t = y[j].unsigned_abs() as usize; + // Codewords with |y_j| = t start at sum_{m < kk - t} V(rem-1, m) within their sign half. + let mut base: u64 = (0..kk - t).map(|m| vv(rem - 1, m)).sum(); + if y[j] < 0 { + base += (vv(rem - 1, kk) + vv(rem, kk)) / 2; + } + index += base; + kk -= t; + } + debug_assert!(total <= u32::MAX as u64); + if total >= 2 { + enc.uint(index as u32, total as u32); + } +} + +pub struct BandEnc<'m, 'e> { + pub m: &'m Mode, + pub enc: &'e mut RangeEncoder, + pub i: usize, + pub spread: usize, + pub remaining_bits: i32, + pvq_v: Vec, + iy: Vec, + tmp: Vec, +} + +impl BandEnc<'_, '_> { + /// `compute_theta` for a mono split (single block): measures the energy split between the + /// two halves and codes the quantised angle with the triangular pdf. + fn compute_theta_split(&mut self, xa: &[f32], xb: &[f32], b: &mut i32, lm: i32) -> (i32, i32, i32) { + let n = xa.len() as i32; + let i = self.i; + let pulse_cap = self.m.log_n[i] + lm * (1 << BITRES); + let offset = (pulse_cap >> 1) - QTHETA_OFFSET; + let qn = compute_qn(n, *b, offset, pulse_cap, false); + let tell = self.enc.tell_frac(); + let mut itheta = 0i32; + if qn != 1 { + let mid = xa.iter().map(|v| v * v).sum::().sqrt(); + let side = xb.iter().map(|v| v * v).sum::().sqrt(); + let ang = (side as f64).atan2(mid as f64) * std::f64::consts::FRAC_2_PI; + let it = ((ang * qn as f64 + 0.5).floor() as i32).clamp(0, qn) as u32; + let q = qn as u32; + let half = q >> 1; + let ft = (half + 1) * (half + 1); + let (fl, fs) = if it <= half { (it * (it + 1) >> 1, it + 1) } else { (ft - ((q + 1 - it) * (q + 2 - it) >> 1), q + 1 - it) }; + self.enc.encode(fl, fl + fs, ft); + itheta = ((it * 16384) / q) as i32; + } + let qalloc = self.enc.tell_frac() - tell; + *b -= qalloc; + let delta = if itheta == 0 { + -16384 + } else if itheta == 16384 { + 16384 + } else { + let imid = bitexact_cos(itheta); + let iside = bitexact_cos(16384 - itheta); + frac_mul16((n - 1) << 7, bitexact_log2tan(iside, imid)) + }; + (itheta, delta, qalloc) + } + + /// `quant_partition` (single block, encoder side). + pub fn quant_partition(&mut self, x: &[f32], mut b: i32, mut lm: i32) { + let n = x.len(); + let i = self.i; + let cache = self.m.cache(lm, i); + let cache_max = cache[cache[0] as usize] as i32; + if lm != -1 && b > cache_max + 12 && n > 2 { + let half = n >> 1; + let (xa, xb) = x.split_at(half); + lm -= 1; + let (itheta, delta, qalloc) = self.compute_theta_split(xa, xb, &mut b, lm); + let mut mbits = 0.max(b.min((b - delta) / 2)); + let mut sbits = b - mbits; + self.remaining_bits -= qalloc; + let mut rebalance = self.remaining_bits; + if mbits >= sbits { + self.quant_partition(xa, mbits, lm); + rebalance = mbits - (rebalance - self.remaining_bits); + if rebalance > 3 << BITRES && itheta != 0 { + sbits += rebalance - (3 << BITRES); + } + self.quant_partition(xb, sbits, lm); + } else { + self.quant_partition(xb, sbits, lm); + rebalance = sbits - (rebalance - self.remaining_bits); + if rebalance > 3 << BITRES && itheta != 16384 { + mbits += rebalance - (3 << BITRES); + } + self.quant_partition(xa, mbits, lm); + } + return; + } + let m = self.m; + let mut q = m.bits2pulses(i, lm, b); + let mut curr_bits = m.pulses2bits(i, lm, q); + self.remaining_bits -= curr_bits; + while self.remaining_bits < 0 && q > 0 { + self.remaining_bits += curr_bits; + q -= 1; + curr_bits = m.pulses2bits(i, lm, q); + self.remaining_bits -= curr_bits; + } + if q != 0 { + let k = get_pulses(q) as usize; + let mut tmp = std::mem::take(&mut self.tmp); + tmp.clear(); + tmp.extend_from_slice(x); + exp_rotation_inv(&mut tmp, n, k, self.spread); + let mut iy = std::mem::take(&mut self.iy); + iy.clear(); + iy.resize(n, 0); + pvq_search(&tmp, k, &mut iy); + encode_pulses(&iy, k, self.enc, &mut self.pvq_v); + self.iy = iy; + self.tmp = tmp; + } + } +} + +/// Codes all band shapes (`quant_all_bands`, encoder side, long blocks). `x` holds `c` channels +/// of `n` normalised MDCT coefficients each. +#[allow(clippy::too_many_arguments)] +pub fn quant_all_bands( + m: &Mode, + start: usize, + end: usize, + x: &[f32], + n_total: usize, + c: usize, + pulses: &[i32; NB_EBANDS], + spread: usize, + dual_stereo: bool, + intensity: usize, + total_bits: i32, + mut balance: i32, + enc: &mut RangeEncoder, + lm: usize, + coded_bands: usize, +) { + let mm = 1usize << lm; + debug_assert!(c == 1 || dual_stereo || intensity <= start, "stereo frames are coded as dual stereo"); + let mut ctx = BandEnc { m, enc, i: 0, spread, remaining_bits: 0, pvq_v: Vec::new(), iy: Vec::new(), tmp: Vec::new() }; + for i in start..end { + ctx.i = i; + let xo = mm * EBANDS[i] as usize; + let n = mm * EBANDS[i + 1] as usize - xo; + let tell = ctx.enc.tell_frac(); + if i != start { + balance -= tell; + } + let remaining_bits = total_bits - tell - 1; + ctx.remaining_bits = remaining_bits; + let b = if i < coded_bands { + let curr_balance = balance / 3.min(coded_bands as i32 - i as i32); + 0.max(16383.min((remaining_bits + 1).min(pulses[i] + curr_balance))) + } else { + 0 + }; + // Bands at and above `intensity` (= coded_bands) carry no bits, so they emit no symbols. + if i < coded_bands || c == 1 { + if c == 2 && dual_stereo && i < intensity { + ctx.quant_partition(&x[xo..xo + n], b / 2, lm as i32); + ctx.quant_partition(&x[n_total + xo..n_total + xo + n], b / 2, lm as i32); + } else { + ctx.quant_partition(&x[xo..xo + n], b, lm as i32); + } + } + balance += pulses[i] + tell; + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::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 = (0..len).map(|i| ((i * 31 % 17) as f32 - 8.0) / 8.0).collect(); + let mut x = orig.clone(); + exp_rotation_inv(&mut x, len, k, SPREAD_NORMAL); + exp_rotation(&mut x, len, k, SPREAD_NORMAL); + for (a, b) in x.iter().zip(&orig) { + assert!((a - b).abs() < 1e-4, "len={len} k={k}"); + } + } + } + } + + /// PVQ decoder of RFC 6716 §4.3.4.2 (as in FilmCraft's decoder). + fn decode_pulses(n: usize, k: usize, dec: &mut RangeDecoder) -> Vec { + let w = k + 1; + let mut v = vec![0u64; (n + 1) * w]; + v[0] = 1; + for nn in 1..=n { + v[nn * w] = 1; + for kk in 1..=k { + v[nn * w + kk] = v[(nn - 1) * w + kk].saturating_add(v[nn * w + kk - 1]).saturating_add(v[(nn - 1) * w + kk - 1]); + } + } + let total = v[n * w + k]; + let mut i = if total >= 2 { dec.uint(total as u32) as u64 } else { 0 }; + let mut y = vec![0i32; n]; + let mut kk = k; + for j in 0..n { + let rem = n - j; + if kk == 0 { + break; + } + let vv = |nn: usize, kx: usize| v[nn * w + kx]; + let mut p = (vv(rem - 1, kk) + vv(rem, kk)) / 2; + let neg = i >= p; + if neg { + i -= p; + } + let k0 = kk; + p -= vv(rem - 1, kk); + while p > i && kk > 0 { + kk -= 1; + p -= vv(rem - 1, kk); + } + let mag = (k0 - kk) as i32; + y[j] = if neg { -mag } else { mag }; + i -= p.min(i); + } + y + } + + fn all_codewords(n: usize, k: usize) -> Vec> { + if n == 0 { + return if k == 0 { vec![vec![]] } else { vec![] }; + } + let mut out = Vec::new(); + for t in -(k as i32)..=(k as i32) { + for mut rest in all_codewords(n - 1, k - t.unsigned_abs() as usize) { + rest.insert(0, t); + out.push(rest); + } + } + out + } + + #[test] + fn pvq_index_roundtrip_exhaustive() { + let mut v = Vec::new(); + for n in 1..6usize { + for k in 1..6usize { + for y in all_codewords(n, k) { + let mut e = RangeEncoder::new(); + encode_pulses(&y, k, &mut e, &mut v); + let buf = e.done(16); + let mut d = RangeDecoder::new(&buf); + assert_eq!(decode_pulses(n, k, &mut d), y, "n={n} k={k}"); + } + } + } + } + + #[test] + fn pvq_search_finds_good_direction() { + let x: Vec = (0..32).map(|i| ((i as f32) * 0.7).sin()).collect(); + let mut y = vec![0i32; 32]; + for k in [1usize, 4, 16, 64] { + pvq_search(&x, k, &mut y); + assert_eq!(y.iter().map(|v| v.unsigned_abs() as usize).sum::(), k); + let xy: f32 = x.iter().zip(&y).map(|(a, &b)| a * b as f32).sum(); + let yy: f32 = y.iter().map(|&b| (b * b) as f32).sum(); + let xx: f32 = x.iter().map(|a| a * a).sum(); + let corr = xy / (xx * yy).sqrt(); + assert!(corr > if k >= 16 { 0.9 } else { 0.2 }, "k={k} corr={corr}"); + } + } +} diff --git a/crates/audio/src/opus/encoder/energy.rs b/crates/audio/src/opus/encoder/energy.rs new file mode 100644 index 000000000..8d4f357d0 --- /dev/null +++ b/crates/audio/src/opus/encoder/energy.rs @@ -0,0 +1,230 @@ +//! Band energy quantisation: Laplace-coded coarse energy with inter/intra prediction, fine +//! energy and the final refinement bits (RFC 6716 §4.3.2, encoder direction). +//! +//! The quantised energies are updated with exactly the decoder's arithmetic so the prediction +//! state never drifts between encoder and decoder. + +use crate::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 = (-64..=64).collect(); + let mut e = RangeEncoder::new(); + let coded: Vec = vals.iter().map(|&v| laplace_encode(&mut e, v, fs, decay)).collect(); + let size = (e.tell() as usize).div_ceil(8); + let buf = e.done(size); + let mut d = RangeDecoder::new(&buf); + for (&v, &c) in vals.iter().zip(&coded) { + assert!(c == v || (c.signum() == v.signum() && c.abs() < v.abs()), "clamp {v} -> {c}"); + assert_eq!(laplace_decode(&mut d, fs, decay), c, "lm={lm} intra={intra} band={band} v={v}"); + } + } + } + } + } +} diff --git a/crates/audio/src/opus/encoder/mdct.rs b/crates/audio/src/opus/encoder/mdct.rs new file mode 100644 index 000000000..6a42cce96 --- /dev/null +++ b/crates/audio/src/opus/encoder/mdct.rs @@ -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, + twiddle: Vec, +} + +impl Fft { + pub fn new(n: usize) -> Fft { + let mut factors = Vec::new(); + let mut m = n; + for p in [4usize, 2, 3, 5] { + while m.is_multiple_of(p) { + factors.push(p); + m /= p; + } + } + assert_eq!(m, 1, "unsupported FFT size {n}"); + let twiddle = (0..n) + .map(|k| { + let ph = -2.0 * std::f64::consts::PI * k as f64 / n as f64; + Cpx { re: ph.cos() as f32, im: ph.sin() as f32 } + }) + .collect(); + Fft { n, factors, twiddle } + } + + pub fn process(&self, data: &mut [Cpx], scratch: &mut Vec) { + scratch.clear(); + scratch.extend_from_slice(&data[..self.n]); + self.rec(&mut data[..self.n], scratch, 1, self.n, 0, 1); + } + + fn rec(&self, out: &mut [Cpx], inp: &[Cpx], in_stride: usize, n: usize, fi: usize, tw_stride: usize) { + let p = self.factors[fi]; + let m = n / p; + if m == 1 { + for (i, o) in out.iter_mut().enumerate().take(p) { + *o = inp[i * in_stride]; + } + } else { + for q in 0..p { + self.rec(&mut out[q * m..(q + 1) * m], &inp[q * in_stride..], in_stride * p, m, fi + 1, tw_stride * p); + } + } + let big = self.n; + let mut t = [Cpx::default(); 5]; + for k in 0..m { + for (q, tq) in t.iter_mut().enumerate().take(p) { + *tq = out[q * m + k].mul(self.twiddle[(q * k * tw_stride) % big]); + } + for q2 in 0..p { + let mut acc = Cpx::default(); + for (q, tq) in t.iter().enumerate().take(p) { + acc = acc.add(tq.mul(self.twiddle[((q * q2) % p) * (big / p)])); + } + out[q2 * m + k] = acc; + } + } + } +} + +/// Forward MDCT producing `n` coefficients from a block of `n + OVERLAP` samples. +pub struct Mdct { + /// Number of coefficients (half the MDCT length). + n: usize, + fft: Fft, + pre: Vec, + post: Vec, + /// Rising half of the low-overlap window (RFC 6716 §4.3.7). + pub window: [f32; OVERLAP], + scale: f32, + z: Vec, + scratch: Vec, + v: Vec, +} + +impl Mdct { + pub fn new(n: usize) -> Mdct { + let n2 = n / 2; + let pi = std::f64::consts::PI; + let pre = (0..n2) + .map(|m| { + let ph = -pi * m as f64 / n as f64; + Cpx { re: ph.cos() as f32, im: ph.sin() as f32 } + }) + .collect(); + let post = (0..n2) + .map(|p| { + let ph = -pi * (p as f64 + 0.25) / n as f64; + Cpx { re: ph.cos() as f32, im: ph.sin() as f32 } + }) + .collect(); + let mut window = [0f32; OVERLAP]; + for (i, w) in window.iter_mut().enumerate() { + let x = std::f64::consts::FRAC_PI_2 * (i as f64 + 0.5) / OVERLAP as f64; + let s = x.sin(); + *w = (std::f64::consts::FRAC_PI_2 * s * s).sin() as f32; + } + Mdct { n, fft: Fft::new(n2), pre, post, window, scale: 2.0 / n as f32, z: Vec::new(), scratch: Vec::new(), v: vec![0.0; n] } + } + + /// Transforms `block` (`n + OVERLAP` samples: the previous frame's last `OVERLAP` samples, then + /// `n` new ones) into `out[..n]`, scaled so that the RFC 6716 decoder's inverse MDCT and + /// overlap-add reproduce the input. + pub fn forward(&mut self, block: &[f32], out: &mut [f32]) { + let n = self.n; + let n2 = n / 2; + debug_assert_eq!(block.len(), n + OVERLAP); + // Place the windowed block into the (virtual) 2n-sample MDCT input u at offset `pad`. + let pad = (n - OVERLAP) / 2; + let w = &self.window; + let u = |t: usize| -> f32 { + if t < pad || t >= pad + n + OVERLAP { + return 0.0; + } + let j = t - pad; + let g = if j < OVERLAP { + w[j] + } else if j >= n { + w[OVERLAP - 1 - (j - n)] + } else { + 1.0 + }; + block[j] * g + }; + // Fold to a DCT-IV input: v = (-c_r - d, a - b_r) over the four quarters of u. + for m in 0..n2 { + self.v[m] = -u(n + n2 - 1 - m) - u(n + n2 + m); + self.v[n2 + m] = u(m) - u(n2 + n2 - 1 - m); + } + // DCT-IV via an n/2-point complex FFT. + self.z.clear(); + for m in 0..n2 { + let wv = Cpx { re: self.v[2 * m], im: self.v[n - 1 - 2 * m] }; + self.z.push(wv.mul(self.pre[m])); + } + self.fft.process(&mut self.z, &mut self.scratch); + for p in 0..n2 { + let y = self.z[p].mul(self.post[p]); + out[2 * p] = y.re * self.scale; + out[n - 1 - 2 * p] = -y.im * self.scale; + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Direct MDCT definition X[k] = sum u[t] cos(pi/N (t + 1/2 + N/2)(k + 1/2)). + fn direct(block: &[f32], n: usize, window: &[f32; OVERLAP]) -> Vec { + let pad = (n - OVERLAP) / 2; + let mut u = vec![0f64; 2 * n]; + for (j, &b) in block.iter().enumerate() { + let g = if j < OVERLAP { + window[j] + } else if j >= n { + window[OVERLAP - 1 - (j - n)] + } else { + 1.0 + }; + u[pad + j] = (b * g) as f64; + } + (0..n) + .map(|k| { + (0..2 * n).map(|t| u[t] * (std::f64::consts::PI / n as f64 * (t as f64 + 0.5 + n as f64 / 2.0) * (k as f64 + 0.5)).cos()).sum::() + * 2.0 + / n as f64 + }) + .collect() + } + + #[test] + fn matches_direct_formula() { + for &n in &[240usize, 960] { + let block: Vec = (0..n + OVERLAP).map(|i| ((i * 37 % 101) as f32 - 50.0) / 7.0).collect(); + let mut m = Mdct::new(n); + let mut out = vec![0f32; n]; + m.forward(&block, &mut out); + let d = direct(&block, n, &m.window); + for k in 0..n { + assert!((out[k] as f64 - d[k]).abs() < 1e-2, "n={n} k={k}: {} vs {}", out[k], d[k]); + } + } + } + + // ---- Replica of the RFC 6716 decoder inverse MDCT (as in FilmCraft's decoder) ---- + + fn imdct_backward(n_full: usize, input: &[f32], out: &mut [f32], window: &[f32]) { + let n2 = n_full >> 1; + let n4 = n_full >> 2; + let t: Vec = (0..n2).map(|i| (2.0 * std::f64::consts::PI * (i as f64 + 0.125) / n_full as f64).cos() as f32).collect(); + let fft = Fft::new(n4); + let mut z = Vec::new(); + for i in 0..n4 { + let x1 = input[2 * i]; + let x2 = input[n2 - 1 - 2 * i]; + let yr = x2 * t[i] + x1 * t[n4 + i]; + let yi = x1 * t[i] - x2 * t[n4 + i]; + z.push(Cpx { re: yi, im: yr }); + } + fft.process(&mut z, &mut Vec::new()); + let overlap = OVERLAP; + let base = overlap >> 1; + let y = &mut out[base..base + n2]; + for i in 0..n4.div_ceil(2) { + let front = z[i]; + let back = z[n4 - 1 - i]; + let (re, im) = (front.im, front.re); + let (t0, t1) = (t[i], t[n4 + i]); + let yr0 = re * t0 + im * t1; + let yi0 = re * t1 - im * t0; + let (re, im) = (back.im, back.re); + let (t0, t1) = (t[n4 - i - 1], t[n2 - i - 1]); + let yr1 = re * t0 + im * t1; + let yi1 = re * t1 - im * t0; + y[2 * i] = yr0; + y[n2 - 1 - 2 * i] = yi0; + y[n2 - 2 - 2 * i] = yr1; + y[2 * i + 1] = yi1; + } + for i in 0..overlap / 2 { + let x1 = out[overlap - 1 - i]; + let x2 = out[i]; + let w1 = window[i]; + let w2 = window[overlap - 1 - i]; + out[i] = x2 * w2 - x1 * w1; + out[overlap - 1 - i] = x2 * w1 + x1 * w2; + } + } + + /// MDCT analysis followed by the decoder's IMDCT + TDAC overlap-add reproduces the input + /// delayed by exactly `OVERLAP` samples. + #[test] + fn perfect_reconstruction_through_decoder_imdct() { + let n = 960; + let frames = 5; + let sig: Vec = (0..n * frames).map(|i| (i as f32 * 0.013).sin() * 1000.0 + ((i * 7919) % 211) as f32 - 105.0).collect(); + let mut m = Mdct::new(n); + let mut hist = vec![0f32; OVERLAP]; + // Decoder memory: [frame output n][folded tail OVERLAP/2] laid out like decode_mem. + let mut mem = vec![0f32; n + OVERLAP]; + let mut decoded = Vec::new(); + for f in 0..frames { + let mut block = hist.clone(); + block.extend_from_slice(&sig[f * n..(f + 1) * n]); + hist.copy_from_slice(&block[n..]); + let mut coefs = vec![0f32; n]; + m.forward(&block, &mut coefs); + // Shift: the previous folded tail (at n..n+OVERLAP/2) moves to 0..OVERLAP/2. + mem.copy_within(n..n + OVERLAP / 2, 0); + imdct_backward(2 * n, &coefs, &mut mem, &m.window); + decoded.extend_from_slice(&mem[..n]); + } + let mut err = 0f64; + let mut sum = 0f64; + for t in OVERLAP..n * frames { + let d = decoded[t] as f64 - sig[t - OVERLAP] as f64; + err += d * d; + sum += (sig[t - OVERLAP] as f64).powi(2); + } + let snr = 10.0 * (sum / err).log10(); + assert!(snr > 80.0, "reconstruction SNR {snr:.1} dB"); + } +} diff --git a/crates/audio/src/opus/encoder/mod.rs b/crates/audio/src/opus/encoder/mod.rs new file mode 100644 index 000000000..da29b3f62 --- /dev/null +++ b/crates/audio/src/opus/encoder/mod.rs @@ -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 = 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, + /// 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, + /// Quantised band energies, exactly as the decoder holds them. + old_band_e: [f32; 2 * NB_EBANDS], + block: Vec, + freq: Vec, + silk: Option, +} + +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::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 { + 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 { + 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 { + 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 { + 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, 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> { + 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::()).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 { + 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 = (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()); + } +} diff --git a/crates/audio/src/opus/encoder/range.rs b/crates/audio/src/opus/encoder/range.rs new file mode 100644 index 000000000..daedd58ae --- /dev/null +++ b/crates/audio/src/opus/encoder/range.rs @@ -0,0 +1,472 @@ +//! Range encoder (RFC 6716 §5.1), the exact inverse of the §4.1 range decoder. +//! +//! Range-coded symbols are written from the front of the buffer and raw bits from the back; the +//! bit accounting ([`RangeEncoder::tell`], [`RangeEncoder::tell_frac`]) matches the decoder's, so +//! every budget decision the decoder takes can be mirrored exactly. + +const SYM_BITS: u32 = 8; +const CODE_BITS: u32 = 32; +const SYM_MAX: u32 = (1 << SYM_BITS) - 1; +const CODE_SHIFT: u32 = CODE_BITS - SYM_BITS - 1; +const CODE_TOP: u32 = 1 << (CODE_BITS - 1); +const CODE_BOT: u32 = CODE_TOP >> SYM_BITS; +/// Bit resolution of [`RangeEncoder::tell_frac`] (1/8 bit). +pub const BITRES: u32 = 3; + +/// Number of bits needed to represent `x` (0 for 0). +#[inline] +pub fn ilog(x: u32) -> i32 { + (32 - x.leading_zeros()) as i32 +} + +#[derive(Clone, Debug)] +pub struct RangeEncoder { + out: Vec, + end: Vec, + end_window: u32, + nend_bits: i32, + nbits_total: i32, + rng: u32, + low: u32, + rem: i32, + ext: u32, +} + +impl Default for RangeEncoder { + fn default() -> Self { + Self::new() + } +} + +impl RangeEncoder { + pub fn new() -> Self { + RangeEncoder { + out: Vec::new(), + end: Vec::new(), + end_window: 0, + nend_bits: 0, + nbits_total: (CODE_BITS + 1) as i32, + rng: CODE_TOP, + low: 0, + rem: -1, + ext: 0, + } + } + + /// Current range (equals the decoder's range after the same symbols). + #[cfg(test)] + pub fn rng(&self) -> u32 { + self.rng + } + + fn carry_out(&mut self, c: u32) { + if c != SYM_MAX { + let carry = c >> SYM_BITS; + if self.rem >= 0 { + self.out.push((self.rem as u32 + carry) as u8); + } + if self.ext > 0 { + let sym = (SYM_MAX + carry) & SYM_MAX; + for _ in 0..self.ext { + self.out.push(sym as u8); + } + self.ext = 0; + } + self.rem = (c & SYM_MAX) as i32; + } else { + self.ext += 1; + } + } + + #[inline] + fn normalize(&mut self) { + while self.rng <= CODE_BOT { + self.carry_out(self.low >> CODE_SHIFT); + self.low = (self.low << SYM_BITS) & (CODE_TOP - 1); + self.rng <<= SYM_BITS; + self.nbits_total += SYM_BITS as i32; + } + } + + /// Encodes the symbol with cumulative range `[fl, fh)` of total `ft`. + pub fn encode(&mut self, fl: u32, fh: u32, ft: u32) { + let r = self.rng / ft; + if fl > 0 { + self.low += self.rng - r * (ft - fl); + self.rng = r * (fh - fl); + } else { + self.rng -= r * (ft - fh); + } + self.normalize(); + } + + /// Like [`Self::encode`] with `ft = 1 << bits`. + pub fn encode_bin(&mut self, fl: u32, fh: u32, bits: u32) { + let r = self.rng >> bits; + let ft = 1u32 << bits; + if fl > 0 { + self.low += self.rng - r * (ft - fl); + self.rng = r * (fh - fl); + } else { + self.rng -= r * (ft - fh); + } + self.normalize(); + } + + /// Encodes a binary symbol whose probability of being 1 is `1/2^logp`. + pub fn bit_logp(&mut self, val: bool, logp: u32) { + let r = self.rng; + let s = r >> logp; + let r2 = r - s; + if val { + self.low += r2; + self.rng = s; + } else { + self.rng = r2; + } + self.normalize(); + } + + /// Encodes symbol `s` with an "inverse CDF" table (`ft = 1 << ftb`). + pub fn icdf(&mut self, s: usize, icdf: &[u8], ftb: u32) { + let r = self.rng >> ftb; + if s > 0 { + self.low += self.rng - r * icdf[s - 1] as u32; + self.rng = r * (icdf[s - 1] as u32 - icdf[s] as u32); + } else { + self.rng -= r * icdf[s] as u32; + } + self.normalize(); + } + + /// Encodes a uniformly distributed integer `fl` in `[0, ft)` (`ft >= 2`). + pub fn uint(&mut self, fl: u32, ft: u32) { + debug_assert!(ft > 1 && fl < ft); + let ft1 = ft - 1; + let ftb = ilog(ft1); + if ftb > 8 { + let ftb = (ftb - 8) as u32; + let top = (ft1 >> ftb) + 1; + let v = fl >> ftb; + self.encode(v, v + 1, top); + self.bits(fl & ((1 << ftb) - 1), ftb); + } else { + self.encode(fl, fl + 1, ft); + } + } + + /// Writes `bits` raw bits (0..=25) to the end of the buffer. + pub fn bits(&mut self, val: u32, bits: u32) { + if bits == 0 { + return; + } + let mut window = self.end_window; + let mut used = self.nend_bits; + if used as u32 + bits > CODE_BITS { + while used >= SYM_BITS as i32 { + self.end.push((window & SYM_MAX) as u8); + window >>= SYM_BITS; + used -= SYM_BITS as i32; + } + } + window |= val << used; + used += bits as i32; + self.end_window = window; + self.nend_bits = used; + self.nbits_total += bits as i32; + } + + /// Number of bits written so far (rounded up), identical to the decoder's `tell`. + #[inline] + pub fn tell(&self) -> i32 { + self.nbits_total - ilog(self.rng) + } + + /// Number of 1/8 bits written so far (rounded up), identical to the decoder's `tell_frac`. + pub fn tell_frac(&self) -> i32 { + let nbits = (self.nbits_total as u32) << BITRES; + let mut l = ilog(self.rng) as u32; + let mut r = self.rng >> (l - 16); + for _ in 0..BITRES { + r = (r * r) >> 15; + let b = r >> 16; + l = (l << 1) | b; + r >>= b; + } + nbits.wrapping_sub(l) as i32 + } + + /// Finalises the stream into a buffer of exactly `size` bytes (range data at the front, raw + /// bits at the back, zero padding in between; a partial raw-bit byte is merged into the byte + /// just before the full raw-bit bytes). + pub fn done(mut self, size: usize) -> Vec { + let mut l = CODE_BITS as i32 - ilog(self.rng); + let mut msk = (CODE_TOP - 1) >> l; + let mut end = (self.low.wrapping_add(msk)) & !msk; + if (end | msk) >= self.low.wrapping_add(self.rng) { + l += 1; + msk >>= 1; + end = (self.low.wrapping_add(msk)) & !msk; + } + while l > 0 { + self.carry_out(end >> CODE_SHIFT); + end = (end << SYM_BITS) & (CODE_TOP - 1); + l -= SYM_BITS as i32; + } + if self.rem >= 0 || self.ext > 0 { + self.carry_out(0); + } + let mut window = self.end_window; + let mut used = self.nend_bits; + while used >= SYM_BITS as i32 { + self.end.push((window & SYM_MAX) as u8); + window >>= SYM_BITS; + used -= SYM_BITS as i32; + } + let mut buf = vec![0u8; size]; + let full = self.end.len().min(size); + let room = size - full; + let n = self.out.len().min(room); + buf[..n].copy_from_slice(&self.out[..n]); + for (k, &b) in self.end.iter().take(full).enumerate() { + buf[size - 1 - k] = b; + } + if used > 0 && room > 0 { + buf[room - 1] |= window as u8; + } + buf + } +} + +#[cfg(test)] +pub(crate) mod dec { + //! Minimal range decoder (RFC 6716 §4.1) used to verify the encoder. + use super::*; + const CODE_EXTRA: u32 = (CODE_BITS - 2) % SYM_BITS + 1; + + pub struct RangeDecoder<'a> { + buf: &'a [u8], + offs: usize, + end_offs: usize, + end_window: u32, + nend_bits: i32, + nbits_total: i32, + rng: u32, + val: u32, + rem: u32, + } + + impl<'a> RangeDecoder<'a> { + pub fn new(buf: &'a [u8]) -> Self { + let mut d = RangeDecoder { + buf, + offs: 0, + end_offs: 0, + end_window: 0, + nend_bits: 0, + nbits_total: (CODE_BITS + 1 - ((CODE_BITS - CODE_EXTRA) / SYM_BITS) * SYM_BITS) as i32, + rng: 1 << CODE_EXTRA, + val: 0, + rem: 0, + }; + d.rem = d.read_byte(); + d.val = d.rng - 1 - (d.rem >> (SYM_BITS - CODE_EXTRA)); + d.normalize(); + d + } + fn read_byte(&mut self) -> u32 { + if self.offs < self.buf.len() { + self.offs += 1; + self.buf[self.offs - 1] as u32 + } else { + 0 + } + } + fn read_byte_from_end(&mut self) -> u32 { + if self.end_offs < self.buf.len() { + self.end_offs += 1; + self.buf[self.buf.len() - self.end_offs] as u32 + } else { + 0 + } + } + fn normalize(&mut self) { + while self.rng <= CODE_BOT { + self.nbits_total += SYM_BITS as i32; + self.rng <<= SYM_BITS; + let sym = self.rem; + self.rem = self.read_byte(); + let sym = ((sym << SYM_BITS) | self.rem) >> (SYM_BITS - CODE_EXTRA); + self.val = ((self.val << SYM_BITS).wrapping_add(SYM_MAX & !sym)) & (CODE_TOP - 1); + } + } + pub fn rng(&self) -> u32 { + self.rng + } + pub fn decode(&mut self, ft: u32) -> u32 { + let ext = self.rng / ft; + let s = self.val / ext; + ft - (s + 1).min(ft) + } + pub fn update(&mut self, fl: u32, fh: u32, ft: u32) { + let ext = self.rng / ft; + let s = ext * (ft - fh); + self.val = self.val.wrapping_sub(s); + self.rng = if fl > 0 { ext * (fh - fl) } else { self.rng - s }; + self.normalize(); + } + pub fn bit_logp(&mut self, logp: u32) -> bool { + let r = self.rng; + let d = self.val; + let s = r >> logp; + let ret = d < s; + if !ret { + self.val = d - s; + } + self.rng = if ret { s } else { r - s }; + self.normalize(); + ret + } + pub fn icdf(&mut self, icdf: &[u8], ftb: u32) -> usize { + let mut s = self.rng; + let d = self.val; + let r = s >> ftb; + let mut ret = 0usize; + let mut t; + loop { + t = s; + s = r * icdf[ret] as u32; + if d >= s { + break; + } + ret += 1; + } + self.val = d - s; + self.rng = t - s; + self.normalize(); + ret + } + pub fn uint(&mut self, ft: u32) -> u32 { + let ft1 = ft - 1; + let ftb = ilog(ft1); + if ftb > 8 { + let ftb = (ftb - 8) as u32; + let top = (ft1 >> ftb) + 1; + let s = self.decode(top); + self.update(s, s + 1, top); + (s << ftb) | self.bits(ftb) + } else { + let s = self.decode(ft); + self.update(s, s + 1, ft); + s + } + } + pub fn bits(&mut self, bits: u32) -> u32 { + let mut window = self.end_window; + let mut available = self.nend_bits; + if (available as u32) < bits { + loop { + window |= self.read_byte_from_end() << available; + available += SYM_BITS as i32; + if available > (CODE_BITS - SYM_BITS) as i32 { + break; + } + } + } + let ret = window & ((1u32 << bits) - 1); + self.end_window = window >> bits; + self.nend_bits = available - bits as i32; + self.nbits_total += bits as i32; + ret + } + pub fn tell(&self) -> i32 { + self.nbits_total - ilog(self.rng) + } + } +} + +#[cfg(test)] +mod tests { + use super::dec::RangeDecoder; + use super::*; + + struct Lcg(u32); + impl Lcg { + fn next(&mut self) -> u32 { + self.0 = self.0.wrapping_mul(1664525).wrapping_add(1013904223); + self.0 >> 8 + } + } + + #[test] + fn starts_at_one_bit() { + let e = RangeEncoder::new(); + assert_eq!(e.tell(), 1); + assert_eq!(e.tell_frac(), 8); + } + + #[test] + fn roundtrip_mixed_symbols_and_tell_matches() { + let mut rng = Lcg(7); + let icdf: [u8; 4] = [200, 120, 30, 0]; + for trial in 0..300 { + let n = 1 + (rng.next() % 300) as usize; + let mut ops = Vec::new(); + for _ in 0..n { + match rng.next() % 5 { + 0 => { + let ft = 2 + rng.next() % 1000; + ops.push((0u8, rng.next() % ft, ft)); + } + 1 => { + let logp = 1 + rng.next() % 15; + ops.push((1, rng.next().is_multiple_of(5) as u32, logp)); + } + 2 => { + let bits = 1 + rng.next() % 24; + ops.push((2, rng.next() & ((1 << bits) - 1), bits)); + } + 3 => { + let ft = 2 + rng.next() % 100_000_000; + ops.push((3, rng.next() % ft, ft)); + } + _ => ops.push((4, rng.next() % 4, 0)), + } + } + let mut e = RangeEncoder::new(); + let mut tells = Vec::new(); + for &(k, v, p) in &ops { + match k { + 0 => e.encode(v, v + 1, p), + 1 => e.bit_logp(v != 0, p), + 2 => e.bits(v, p), + 3 => e.uint(v, p), + _ => e.icdf(v as usize, &icdf, 8), + } + tells.push(e.tell()); + } + // Exactly the number of bytes `tell` says were used must suffice. + let size = (e.tell() as usize).div_ceil(8); + let final_rng = e.rng(); + let buf = e.done(size); + assert_eq!(buf.len(), size); + let mut d = RangeDecoder::new(&buf); + for (i, &(k, v, p)) in ops.iter().enumerate() { + let got = match k { + 0 => { + let s = d.decode(p); + d.update(s, s + 1, p); + s + } + 1 => d.bit_logp(p) as u32, + 2 => d.bits(p), + 3 => d.uint(p), + _ => d.icdf(&icdf, 8) as u32, + }; + assert_eq!(got, v, "trial {trial} op {i} kind {k}"); + assert_eq!(d.tell(), tells[i], "tell mismatch trial {trial} op {i}"); + } + assert_eq!(d.rng(), final_rng); + } + } +} diff --git a/crates/audio/src/opus/encoder/rate.rs b/crates/audio/src/opus/encoder/rate.rs new file mode 100644 index 000000000..c360f9da7 --- /dev/null +++ b/crates/audio/src/opus/encoder/rate.rs @@ -0,0 +1,561 @@ +//! Pulse cache and bit allocation (RFC 6716 §4.3.3, §4.3.4.1). +//! +//! The allocation must be bit-identical to the decoder's, so this mirrors FilmCraft's +//! spec-derived `filmcraft-opus` decoder (first-party, MIT OR Apache-2.0) step by step; the only +//! difference is that the encoder *writes* the band-skip, intensity and dual-stereo symbols. + +use crate::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, + cache_bits: Vec, + caps: Vec, +} + +impl Mode { + pub fn new() -> Mode { + let mut log_n = [0i32; NB_EBANDS]; + for (i, l) in log_n.iter_mut().enumerate() { + *l = log2_frac((EBANDS[i + 1] - EBANDS[i]) as u32, BITRES); + } + let (cache_index, cache_bits) = compute_pulse_cache(MAX_LM as i32); + let mut m = Mode { log_n, cache_index, cache_bits, caps: Vec::new() }; + m.caps = compute_caps(&m); + m + } + + #[inline] + pub fn cache(&self, lm: i32, band: usize) -> &[u8] { + let idx = self.cache_index[((lm + 1) as usize) * NB_EBANDS + band]; + &self.cache_bits[idx.max(0) as usize..] + } + + pub fn bits2pulses(&self, band: usize, lm: i32, bits: i32) -> i32 { + let cache = self.cache(lm, band); + let mut lo = 0i32; + let mut hi = cache[0] as i32; + let bits = bits - 1; + for _ in 0..LOG_MAX_PSEUDO { + let mid = (lo + hi + 1) >> 1; + if cache[mid as usize] as i32 >= bits { + hi = mid; + } else { + lo = mid; + } + } + let lo_bits = if lo == 0 { -1 } else { cache[lo as usize] as i32 }; + if bits - lo_bits <= cache[hi as usize] as i32 - bits { lo } else { hi } + } + + pub fn pulses2bits(&self, band: usize, lm: i32, pulses: i32) -> i32 { + if pulses == 0 { 0 } else { self.cache(lm, band)[pulses as usize] as i32 + 1 } + } + + /// Per-band allocation caps for this frame (`init_caps`). + pub fn init_caps(&self, lm: usize, c: usize) -> [i32; NB_EBANDS] { + let mut cap = [0i32; NB_EBANDS]; + for (i, cp) in cap.iter_mut().enumerate() { + let n = (EBANDS[i + 1] - EBANDS[i]) << lm; + *cp = (self.caps[NB_EBANDS * (2 * lm + c - 1) + i] as i32 + 64) * c as i32 * n >> 2; + } + cap + } +} + +fn compute_pulse_cache(lm_max: i32) -> (Vec, Vec) { + let nb = NB_EBANDS; + let rows = (lm_max + 2) as usize; + let mut index = vec![-1i16; nb * rows]; + let mut entries: Vec<(usize, i32, usize)> = Vec::new(); + let mut curr = 0usize; + for i in 0..rows { + for j in 0..nb { + let n = (((EBANDS[j + 1] - EBANDS[j]) << i) >> 1) as usize; + index[i * nb + j] = -1; + 'search: for k in 0..=i { + for nn in 0..nb { + if k == i && nn >= j { + break; + } + if n == (((EBANDS[nn + 1] - EBANDS[nn]) << k) >> 1) as usize { + index[i * nb + j] = index[k * nb + nn]; + break 'search; + } + } + } + if index[i * nb + j] == -1 && n != 0 { + let mut k = 0; + while k < MAX_PSEUDO && fits_in32(n, get_pulses(k + 1)) { + k += 1; + } + index[i * nb + j] = curr as i16; + entries.push((n, k, curr)); + curr += (k + 1) as usize; + } + } + } + let mut bits = vec![0u8; curr]; + for &(n, k, off) in &entries { + let maxk = get_pulses(k) as usize; + let tmp = required_bits(n, maxk, BITRES); + for j in 1..=k { + bits[off + j as usize] = (tmp[get_pulses(j) as usize] - 1) as u8; + } + bits[off] = k as u8; + } + (index, bits) +} + +fn required_bits(n: usize, maxk: usize, frac: u32) -> Vec { + let mut out = vec![0i32; maxk + 1]; + if n == 1 { + for o in out.iter_mut().skip(1) { + *o = 1 << frac; + } + } else { + for (k, o) in out.iter_mut().enumerate().skip(1) { + *o = log2_frac(pvq_v(n, k) as u32, frac); + } + } + out +} + +fn compute_caps(m: &Mode) -> Vec { + let nb = NB_EBANDS; + let mut caps = vec![0u8; (MAX_LM + 1) * 2 * nb]; + for i in 0..=MAX_LM as i32 { + for c in 1..=2i32 { + for j in 0..nb { + let mut n0 = EBANDS[j + 1] - EBANDS[j]; + + let max_bits = if n0 << i == 1 { + 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); + } +} diff --git a/crates/audio/src/opus/encoder/resample.rs b/crates/audio/src/opus/encoder/resample.rs new file mode 100644 index 000000000..9045cb86d --- /dev/null +++ b/crates/audio/src/opus/encoder/resample.rs @@ -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 { + let wc = 2.0 * std::f64::consts::PI * fc / 48_000.0; + let i0b = bessel_i0(beta); + let mut h: Vec = (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, + taps: usize, + /// The last `taps - 1` input samples followed by the current frame. + buf: Vec, +} + +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) { + 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 { + (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::() / 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}"); + } +} diff --git a/crates/audio/src/opus/encoder/silk/analysis.rs b/crates/audio/src/opus/encoder/silk/analysis.rs new file mode 100644 index 000000000..439331076 --- /dev/null +++ b/crates/audio/src/opus/encoder/silk/analysis.rs @@ -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 { + 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 = 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 { + 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) { + 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> { + 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, 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 = (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::() - 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] = [<P_TAPS_0, <P_TAPS_1, <P_TAPS_2]; + let mut best = (f64::MAX, 0usize, [0usize; 4], [0f64; 4]); + for (p, book) in books.iter().enumerate() { + let mut total = icdf_bits(<P_PERIODICITY, p) * 1e-3 * e.iter().sum::() / 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 = (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 = 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]); + } + } + } +} diff --git a/crates/audio/src/opus/encoder/silk/lpc.rs b/crates/audio/src/opus/encoder/silk/lpc.rs new file mode 100644 index 000000000..e4e9ae624 --- /dev/null +++ b/crates/audio/src/opus/encoder/silk/lpc.rs @@ -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::(); + min_center += dmin[idx] >> 1; + let mut max_center = 32768 - dmin[idx + 1..=l].iter().sum::(); + 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 = (1..=16).map(|k| k * 32768 / 17).collect(); + let a = nlsf2a(&nlsf); + assert!(a.iter().all(|v| v.abs() < 600), "{a:?}"); + let nlsf: Vec = (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); + } +} diff --git a/crates/audio/src/opus/encoder/silk/mod.rs b/crates/audio/src/opus/encoder/silk/mod.rs new file mode 100644 index 000000000..d82e537e0 --- /dev/null +++ b/crates/audio/src/opus/encoder/silk/mod.rs @@ -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, +} + +#[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, + /// 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::() / 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 = 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::() / ltp.err.iter().sum::().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) { + 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, <P_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, <P_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; 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]) -> 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 = (0..l).map(|i| 0.5 * (input[0][i] + input[1][i])).collect(); + let side: Vec = (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 { + 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], 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; + } +} diff --git a/crates/audio/src/opus/encoder/silk/nlsf.rs b/crates/audio/src/opus/encoder/silk/nlsf.rs new file mode 100644 index 000000000..59fca68cb --- /dev/null +++ b/crates/audio/src/opus/encoder/silk/nlsf.rs @@ -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 = (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 = (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 = (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]); + } + } +} diff --git a/crates/audio/src/opus/encoder/silk/pulses.rs b/crates/audio/src/opus/encoder/silk/pulses.rs new file mode 100644 index 000000000..7ec4fa6bf --- /dev/null +++ b/crates/audio/src/opus/encoder/silk/pulses.rs @@ -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); + } + } + } + } +} diff --git a/crates/audio/src/opus/encoder/silk/tables.rs b/crates/audio/src/opus/encoder/silk/tables.rs new file mode 100644 index 000000000..a4789a414 --- /dev/null +++ b/crates/audio/src/opus/encoder/silk/tables.rs @@ -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(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])]; diff --git a/crates/audio/src/opus/encoder/tables.rs b/crates/audio/src/opus/encoder/tables.rs new file mode 100644 index 000000000..1985202cc --- /dev/null +++ b/crates/audio/src/opus/encoder/tables.rs @@ -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]; diff --git a/crates/audio/src/opus/mod.rs b/crates/audio/src/opus/mod.rs index cf2accd33..b67e1ef3d 100644 --- a/crates/audio/src/opus/mod.rs +++ b/crates/audio/src/opus/mod.rs @@ -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 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 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], diff --git a/crates/audio/tests/codecs/main.rs b/crates/audio/tests/codecs/main.rs index 5db6af4b0..279246dc6 100644 --- a/crates/audio/tests/codecs/main.rs +++ b/crates/audio/tests/codecs/main.rs @@ -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; diff --git a/crates/audio/tests/codecs/opusenc/common.rs b/crates/audio/tests/codecs/opusenc/common.rs new file mode 100644 index 000000000..571c9826c --- /dev/null +++ b/crates/audio/tests/codecs/opusenc/common.rs @@ -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 { + 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 { + 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]) -> Vec { + let c = planar.len(); + let len = planar[0].len(); + let mut pcm = vec![0f32; len * c]; + for (ch, p) in planar.iter().enumerate() { + for (i, &v) in p.iter().enumerate() { + pcm[i * c + ch] = v; + } + } + pcm +} + +/// Encodes planar input (padding the tail so the encoder delay is flushed). +pub fn encode(mut enc: OpusEncoder, planar: &[Vec]) -> (OpusEncoder, Vec>) { + 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]) -> Vec> { + let mut dec = craft_audio::opus::Decoder::new(head).expect("OpusHead"); + dec.set_trim_pre_skip(true); + let mut out: Vec> = vec![Vec::new(); dec.channels()]; + for p in packets { + let planes = dec.decode(Some(p)).expect("decode"); + for (o, p) in out.iter_mut().zip(planes) { + o.extend_from_slice(&p); + } + } + 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]) -> f64 { + packets.iter().map(|p| p.len() * 8).sum::() 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 { + let n = x.len().next_power_of_two(); + let mut re: Vec = 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, 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], len: usize) -> Vec { + let mut out = Vec::new(); + ogg_page(&mut out, &enc.opus_head(), 0x02, 0, 0); + let vendor = b"effectcraft-opusenc"; + let mut tags = b"OpusTags".to_vec(); + tags.extend_from_slice(&(vendor.len() as u32).to_le_bytes()); + tags.extend_from_slice(vendor); + tags.extend_from_slice(&0u32.to_le_bytes()); + ogg_page(&mut out, &tags, 0, 0, 1); + let pre_skip = enc.pre_skip() as u64; + for (i, p) in packets.iter().enumerate() { + let last = i + 1 == packets.len(); + let granule = if last { pre_skip + len as u64 } else { ((i + 1) * FRAME) as u64 }; + ogg_page(&mut out, p, if last { 0x04 } else { 0 }, granule, i as u32 + 2); + } + out +} diff --git a/crates/audio/tests/codecs/opusenc/ffmpeg.rs b/crates/audio/tests/codecs/opusenc/ffmpeg.rs new file mode 100644 index 000000000..b428582c1 --- /dev/null +++ b/crates/audio/tests/codecs/opusenc/ffmpeg.rs @@ -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 { + 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> { + 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 = 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> = 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 = (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"); + } + } + } +} diff --git a/crates/audio/tests/codecs/opusenc/mod.rs b/crates/audio/tests/codecs/opusenc/mod.rs new file mode 100644 index 000000000..f2a886463 --- /dev/null +++ b/crates/audio/tests/codecs/opusenc/mod.rs @@ -0,0 +1,5 @@ +//! EffectCraft's `effectcraft-opusenc` integration tests. + +mod common; +mod ffmpeg; +mod modes; diff --git a/crates/audio/tests/codecs/opusenc/modes.rs b/crates/audio/tests/codecs/opusenc/modes.rs new file mode 100644 index 000000000..7060173d4 --- /dev/null +++ b/crates/audio/tests/codecs/opusenc/modes.rs @@ -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>, +} + +/// Encodes `planar`, decodes, and measures each channel against `reference` (searching a +/// residual delay of up to ±6 samples). +fn run(enc: OpusEncoder, planar: &[Vec], reference: &[Vec]) -> Vec { + 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 = 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> = 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 = sp[40_000..].iter().map(|v| v * 0.7).collect(); + let mu = music(LEN); + let cases: Vec<(&str, Vec, Vec)> = 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 = (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 = (0..len).map(|i| if (i / 60) % 2 == 0 { 0.999 } else { -0.999 }).collect(); + let clicks: Vec = (0..len).map(|i| if i % 4_801 == 0 { 1.0 } else { 0.0 }).collect(); + let dc = vec![0.7f32; len]; + let nan: Vec = (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> = 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"); + } + } +} diff --git a/crates/audio/tests/codecs/track_decoder.rs b/crates/audio/tests/codecs/track_decoder.rs index 0607db4ac..fee8f1f97 100644 --- a/crates/audio/tests/codecs/track_decoder.rs +++ b/crates/audio/tests/codecs/track_decoder.rs @@ -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 { @@ -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> = 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(¶ms, &[]).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()); +}