Something went wrong. Try again.
Monorepo for PackTools - tooling for Minecraft Modpacks!
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162//! Downsampling for `audioSampleRate`: a windowed-sinc polyphase filter in plain scalar `f64`.//!//! Hand-written rather than a crate: FFT resamplers pick SIMD code at runtime, so two machines of//! one architecture could write different bytes, and the per-file cache is shared between//! machines. Here the output depends only on the input and the build.//!//! The ratio is reduced to `up/down`; each output sample sits at a whole input position plus one//! of `up` phases, and every phase has a precomputed Kaiser-windowed sinc row, normalized to unit//! gain. The cutoff sits at 95% of the new Nyquist frequency. Only downsampling is supported://! upsampling never saves bytes.
const ZERO_CROSSINGS: f64 = 16.0;const KAISER_BETA: f64 = 8.6;const CUTOFF: f64 = 0.95;/// Phase tables past this size (an unusual pair of rates) fall back to computing each row.const MAX_PHASES: u64 = 4096;
fn gcd(a: u64, b: u64) -> u64 { if b == 0 { a } else { gcd(b, a % b) }}
/// The zeroth-order modified Bessel function, by its series.fn bessel_i0(x: f64) -> f64 { let mut sum = 1.0; let mut term = 1.0; let half = x / 2.0; for k in 1..64 { term *= (half / k as f64) * (half / k as f64); sum += term; if term < sum * 1e-17 { break; } } sum}
struct Filter { up: u64, down: u64, /// Taps on each side of the centre, in input samples. half: i64, /// Lowpass cutoff in cycles per input sample. fc: f64, norm: f64,}
impl Filter { fn new(from: u32, to: u32) -> Filter { let g = gcd(u64::from(from), u64::from(to)); let (up, down) = (u64::from(to) / g, u64::from(from) / g); let fc = 0.5 * (up as f64 / down as f64) * CUTOFF; let half = (ZERO_CROSSINGS / (2.0 * fc)).ceil() as i64; Filter { up, down, half, fc, norm: bessel_i0(KAISER_BETA) } }
/// The weights for an output at input position `centre + phase/up`, normalized to sum 1. fn row(&self, phase: u64) -> Vec<f64> { let offset = phase as f64 / self.up as f64; let mut row: Vec<f64> = (-self.half + 1..=self.half) .map(|k| { let x = k as f64 - offset; let window = 1.0 - (x / self.half as f64).powi(2); if window <= 0.0 { return 0.0; } let sinc = if x == 0.0 { 1.0 } else { let a = std::f64::consts::PI * 2.0 * self.fc * x; a.sin() / a }; 2.0 * self.fc * sinc * bessel_i0(KAISER_BETA * window.sqrt()) / self.norm }) .collect(); let sum: f64 = row.iter().sum(); if sum != 0.0 { row.iter_mut().for_each(|w| *w /= sum); } row }}
/// `channels` (one Vec per channel, all the same length) at `from` Hz, downsampled to `to` Hz./// Returns the input unchanged when `to` is not lower than `from`.pub fn downsample(channels: Vec<Vec<f32>>, from: u32, to: u32) -> Vec<Vec<f32>> { if to == 0 || to >= from { return channels; } let f = Filter::new(from, to); let frames = channels.first().map_or(0, Vec::len) as u64; let out_len = (frames * f.up).div_ceil(f.down); let table: Option<Vec<Vec<f64>>> = (f.up <= MAX_PHASES).then(|| (0..f.up).map(|p| f.row(p)).collect()); channels .iter() .map(|input| { (0..out_len) .map(|n| { let position = n * f.down; let centre = (position / f.up) as i64; let phase = position % f.up; let computed; let row = match &table { Some(t) => &t[phase as usize], None => { computed = f.row(phase); &computed } }; let mut acc = 0.0f64; for (j, w) in row.iter().enumerate() { let k = centre - f.half + 1 + j as i64; if k >= 0 && (k as u64) < frames { acc += f64::from(input[k as usize]) * w; } } acc as f32 }) .collect() }) .collect()}
#[cfg(test)]mod tests { use super::*;
fn sine(hz: f64, rate: u32, seconds: f64) -> Vec<f32> { (0..(f64::from(rate) * seconds) as usize) .map(|i| (2.0 * std::f64::consts::PI * hz * i as f64 / f64::from(rate)).sin() as f32 * 0.5) .collect() }
fn rms(x: &[f32]) -> f64 { (x.iter().map(|v| f64::from(*v).powi(2)).sum::<f64>() / x.len() as f64).sqrt() }
#[test] fn lengths_and_no_upsampling() { let out = downsample(vec![vec![0.0; 44_100], vec![0.0; 44_100]], 44_100, 32_000); assert_eq!((out.len(), out[0].len()), (2, 32_000)); assert_eq!(downsample(vec![vec![1.0; 10]], 22_050, 32_000)[0].len(), 10); }
#[test] fn passes_the_band_and_removes_what_would_alias() { // 1 kHz keeps its level; 20 kHz (above the new 16 kHz Nyquist) is filtered out. let low = downsample(vec![sine(1000.0, 48_000, 0.5)], 48_000, 32_000); let high = downsample(vec![sine(20_000.0, 48_000, 0.5)], 48_000, 32_000); let middle = |x: &[f32]| x[2000..x.len() - 2000].to_vec(); assert!((rms(&middle(&low[0])) - 0.5 / 2f64.sqrt()).abs() < 0.01, "{}", rms(&low[0])); assert!(rms(&middle(&high[0])) < 0.005, "{}", rms(&high[0])); }
#[test] fn deterministic_and_odd_ratios_work() { let x = sine(440.0, 44_100, 0.2); assert_eq!(downsample(vec![x.clone()], 44_100, 32_000), downsample(vec![x.clone()], 44_100, 32_000)); // 44,100 to 31,999 has a phase table too large to precompute. assert_eq!(downsample(vec![x], 44_100, 31_999)[0].len(), 6400); }}