255 lines
9.1 KiB
Rust
255 lines
9.1 KiB
Rust
#![allow(dead_code)]
|
|
|
|
use std::path::PathBuf;
|
|
|
|
/// Plain SNR (dB) of `test` against `reference` over the common length.
|
|
pub fn snr_db(reference: &[f32], test: &[f32]) -> f64 {
|
|
let n = reference.len().min(test.len());
|
|
let mut s = 0f64;
|
|
let mut e = 0f64;
|
|
for i in 0..n {
|
|
s += (reference[i] as f64).powi(2);
|
|
e += (reference[i] as f64 - test[i] as f64).powi(2);
|
|
}
|
|
if e == 0.0 { 200.0 } else { 10.0 * (s.max(1e-20) / e).log10() }
|
|
}
|
|
|
|
const NBANDS: usize = 21;
|
|
const BANDS: [usize; NBANDS + 1] = [0, 2, 4, 6, 8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 68, 80, 96, 120, 156, 200];
|
|
const WIN: usize = 480;
|
|
const STEP: usize = 120;
|
|
const NFREQ: usize = 240;
|
|
/// Calibration: white noise at 48 dB SNR maps to quality 0 (RFC 6716 §6.1).
|
|
const CAL: f64 = 0.144;
|
|
|
|
fn spectra(x: &[f32], channels: usize, ch: usize, rate: u32) -> Vec<[f64; NFREQ]> {
|
|
// Work at 48 kHz bin spacing (100 Hz per bin); lower rates use proportionally fewer bins.
|
|
let scale = 48000 / rate as usize;
|
|
let win = WIN / scale;
|
|
let step = STEP / scale;
|
|
let n = x.len() / channels;
|
|
let nframes = if n >= win { (n - win) / step + 1 } else { 0 };
|
|
let w: Vec<f64> = (0..win).map(|i| 0.5 - 0.5 * (2.0 * std::f64::consts::PI * (i as f64 + 0.5) / win as f64).cos()).collect();
|
|
let nb = NFREQ / scale;
|
|
let mut out = Vec::with_capacity(nframes);
|
|
let mut buf = vec![(0f64, 0f64); win];
|
|
for f in 0..nframes {
|
|
for i in 0..win {
|
|
buf[i] = (x[(f * step + i) * channels + ch] as f64 * w[i], 0.0);
|
|
}
|
|
let spec = fft(&buf);
|
|
let mut ps = [0f64; NFREQ];
|
|
for (k, p) in ps.iter_mut().enumerate().take(nb) {
|
|
*p = spec[k].0 * spec[k].0 + spec[k].1 * spec[k].1;
|
|
}
|
|
out.push(ps);
|
|
}
|
|
out
|
|
}
|
|
|
|
/// Recursive mixed-radix DFT (sizes with factors 2, 3, 5).
|
|
fn fft(x: &[(f64, f64)]) -> Vec<(f64, f64)> {
|
|
let n = x.len();
|
|
if n == 1 {
|
|
return x.to_vec();
|
|
}
|
|
let p = [2usize, 3, 5].into_iter().find(|p| n.is_multiple_of(*p)).unwrap_or(n);
|
|
let m = n / p;
|
|
let subs: Vec<Vec<(f64, f64)>> = (0..p).map(|r| fft(&x.iter().skip(r).step_by(p).copied().collect::<Vec<_>>())).collect();
|
|
let mut out = vec![(0f64, 0f64); n];
|
|
for (k, o) in out.iter_mut().enumerate() {
|
|
let mut acc = (0f64, 0f64);
|
|
for (r, sub) in subs.iter().enumerate() {
|
|
let v = sub[k % m];
|
|
let ang = -2.0 * std::f64::consts::PI * (r * k) as f64 / n as f64;
|
|
let (c, s) = (ang.cos(), ang.sin());
|
|
acc.0 += v.0 * c - v.1 * s;
|
|
acc.1 += v.0 * s + v.1 * c;
|
|
}
|
|
*o = acc;
|
|
}
|
|
out
|
|
}
|
|
|
|
/// Opus quality metric in the spirit of `opus_compare` (RFC 6716 §6.1): a per-band,
|
|
/// masking-weighted spectral distortion between reference and test, mapped so that 100 means
|
|
/// identical output and 0 is the pass threshold (white noise at 48 dB SNR). Phase-insensitive
|
|
/// (compares power spectra), as the RFC requires tolerance to resampler phase differences.
|
|
pub fn opus_quality(reference: &[f32], test: &[f32], channels: usize, rate: u32) -> f64 {
|
|
let n = reference.len().min(test.len());
|
|
let (reference, test) = (&reference[..n], &test[..n]);
|
|
let scale = 48000 / rate as usize;
|
|
let mut total_err = 0f64;
|
|
let mut count = 0usize;
|
|
for ch in 0..channels {
|
|
let xr = spectra(reference, channels, ch, rate);
|
|
let xt = spectra(test, channels, ch, rate);
|
|
let nf = xr.len();
|
|
// Band energies of the reference with temporal and frequency spreading (masking).
|
|
let mut mask = vec![[0f64; NBANDS]; nf];
|
|
for f in 0..nf {
|
|
for b in 0..NBANDS {
|
|
let (lo, hi) = (BANDS[b] / scale, (BANDS[b + 1] / scale).max(BANDS[b] / scale + 1));
|
|
let e: f64 = xr[f][lo..hi.min(NFREQ / scale)].iter().sum::<f64>() / (hi - lo) as f64;
|
|
mask[f][b] = e;
|
|
}
|
|
}
|
|
for f in 0..nf {
|
|
for b in 1..NBANDS {
|
|
mask[f][b] += 0.1 * mask[f][b - 1];
|
|
}
|
|
for b in (0..NBANDS - 1).rev() {
|
|
mask[f][b] += 0.03 * mask[f][b + 1];
|
|
}
|
|
}
|
|
for f in 1..nf {
|
|
for b in 0..NBANDS {
|
|
mask[f][b] += 0.5 * mask[f - 1][b];
|
|
}
|
|
}
|
|
for f in (0..nf.saturating_sub(1)).rev() {
|
|
for b in 0..NBANDS {
|
|
mask[f][b] += 0.1 * mask[f + 1][b];
|
|
}
|
|
}
|
|
for f in 0..nf {
|
|
let mut ef = 0f64;
|
|
let mut nb = 0;
|
|
for b in 0..NBANDS {
|
|
let (lo, hi) = (BANDS[b] / scale, BANDS[b + 1] / scale);
|
|
if hi <= lo {
|
|
continue;
|
|
}
|
|
let m = mask[f][b] + 1e-9 * (WIN * WIN) as f64 / (scale * scale) as f64;
|
|
let mut eb = 0f64;
|
|
for k in lo..hi {
|
|
let d = xt[f][k].sqrt() - xr[f][k].sqrt();
|
|
eb += d * d;
|
|
}
|
|
eb /= (hi - lo) as f64 * m;
|
|
ef += eb;
|
|
nb += 1;
|
|
}
|
|
total_err += (ef / nb.max(1) as f64).powi(2);
|
|
count += 1;
|
|
}
|
|
}
|
|
let err = (total_err / count.max(1) as f64).sqrt() / CAL;
|
|
100.0 * (1.0 - 0.5 * (1.0 + err).ln() / 1.13f64.ln())
|
|
}
|
|
|
|
/// ffmpeg (see `filmcraft_testkit::oracle`); `None` after printing `SKIPPED`.
|
|
pub fn ffmpeg() -> Option<PathBuf> {
|
|
filmcraft_testkit::ffmpeg_or_skip("opus oracle")
|
|
}
|
|
|
|
pub fn fixtures() -> PathBuf {
|
|
filmcraft_testkit::fixtures_dir("opus/oracle")
|
|
}
|
|
|
|
/// Minimal Ogg demuxer: returns (packets, last granule position).
|
|
pub fn ogg_packets(data: &[u8]) -> (Vec<Vec<u8>>, i64) {
|
|
let mut packets = Vec::new();
|
|
let mut cur = Vec::new();
|
|
let mut pos = 0;
|
|
let mut granule = 0i64;
|
|
while pos + 27 <= data.len() && &data[pos..pos + 4] == b"OggS" {
|
|
let g = i64::from_le_bytes(data[pos + 6..pos + 14].try_into().unwrap());
|
|
let nseg = data[pos + 26] as usize;
|
|
let table = &data[pos + 27..pos + 27 + nseg];
|
|
let mut body = pos + 27 + nseg;
|
|
for &l in table {
|
|
cur.extend_from_slice(&data[body..body + l as usize]);
|
|
body += l as usize;
|
|
if l < 255 {
|
|
packets.push(std::mem::take(&mut cur));
|
|
}
|
|
}
|
|
if g >= 0 {
|
|
granule = g;
|
|
}
|
|
pos = body;
|
|
}
|
|
(packets, granule)
|
|
}
|
|
|
|
pub struct Rng(pub u64);
|
|
impl Rng {
|
|
pub fn next(&mut self) -> f32 {
|
|
self.0 ^= self.0 << 13;
|
|
self.0 ^= self.0 >> 7;
|
|
self.0 ^= self.0 << 17;
|
|
((self.0 >> 40) as f32 / (1u64 << 24) as f32) * 2.0 - 1.0
|
|
}
|
|
}
|
|
|
|
/// Logarithmic sine sweep 40 Hz -> 20 kHz at 48 kHz.
|
|
pub fn sweep(n: usize, amp: f32) -> Vec<f32> {
|
|
let (f0, f1) = (40f64, 20000f64);
|
|
let t_total = n as f64 / 48000.0;
|
|
let k = (f1 / f0).ln();
|
|
(0..n)
|
|
.map(|i| {
|
|
let t = i as f64 / 48000.0;
|
|
((2.0 * std::f64::consts::PI * f0 * t_total / k * ((t / t_total * k).exp() - 1.0)).sin() as f32) * amp
|
|
})
|
|
.collect()
|
|
}
|
|
|
|
/// Music-like: harmonic chords with vibrato changing every 0.5 s, decaying noise bursts.
|
|
pub fn music(n: usize, seed: u64) -> Vec<f32> {
|
|
let mut r = Rng(seed);
|
|
let notes = [220.0f32, 277.2, 329.6, 440.0, 349.2, 523.3, 392.0, 659.3];
|
|
(0..n)
|
|
.map(|i| {
|
|
let t = i as f32 / 48000.0;
|
|
let chord = (t * 2.0) as usize;
|
|
let mut s = 0.0;
|
|
for k in 0..3 {
|
|
let f = notes[(chord * 3 + k * 2) % notes.len()] * (1.0 + 0.003 * (2.0 * std::f32::consts::PI * 5.0 * t).sin());
|
|
for h in 1..6 {
|
|
s += (2.0 * std::f32::consts::PI * f * h as f32 * t).sin() * 0.08 / h as f32;
|
|
}
|
|
}
|
|
let env = (-(t * 2.0).fract() * 3.0).exp();
|
|
s * (0.5 + 0.5 * env) + 0.02 * r.next() * env
|
|
})
|
|
.collect()
|
|
}
|
|
|
|
/// Speech-like: glottal pulses through formant resonators with a syllabic envelope.
|
|
pub fn speech(n: usize, seed: u64) -> Vec<f32> {
|
|
let mut r = Rng(seed);
|
|
let fs = 48000.0f32;
|
|
let formants = [(700.0f32, 130.0f32), (1220.0, 70.0), (2600.0, 160.0)];
|
|
let mut states = [[0f32; 2]; 3];
|
|
let mut phase = 0f32;
|
|
let mut out = Vec::with_capacity(n);
|
|
for i in 0..n {
|
|
let t = i as f32 / fs;
|
|
let f0 = 120.0 + 30.0 * (2.0 * std::f32::consts::PI * 0.7 * t).sin();
|
|
phase += f0 / fs;
|
|
let pulse = if phase >= 1.0 {
|
|
phase -= 1.0;
|
|
1.0
|
|
} else {
|
|
0.0
|
|
};
|
|
let syll = (2.0 * std::f32::consts::PI * 3.0 * t).sin().max(0.0);
|
|
let voiced = if (t * 1.3).fract() < 0.8 { syll } else { 0.0 };
|
|
let fric = if (t * 1.3).fract() > 0.85 { 0.05 * r.next() } else { 0.0 };
|
|
let mut s = 0f32;
|
|
for (j, &(f, bw)) in formants.iter().enumerate() {
|
|
let rr = (-std::f32::consts::PI * bw / fs).exp();
|
|
let a1 = 2.0 * rr * (2.0 * std::f32::consts::PI * f / fs).cos();
|
|
let a2 = -rr * rr;
|
|
let y = pulse * voiced + a1 * states[j][0] + a2 * states[j][1];
|
|
states[j][1] = states[j][0];
|
|
states[j][0] = y;
|
|
s += y;
|
|
}
|
|
out.push(s * 0.05 + fric);
|
|
}
|
|
out
|
|
}
|