diff --git a/crates/dsp/src/offline.rs b/crates/dsp/src/offline.rs index 8bbc80a6..6884c4ec 100644 --- a/crates/dsp/src/offline.rs +++ b/crates/dsp/src/offline.rs @@ -216,16 +216,95 @@ fn resample_ratio(ch: &[Vec], ratio: f64, out_len: usize) -> Vec> ch.iter().map(|c| resample_channel(c, ratio, out_len, &table)).collect() } +/// Most kernel weights a [`Polyphase`] bank may hold (16 MiB of f32). +const MAX_POLYPHASE_WEIGHTS: usize = 1 << 22; + +/// The windowed-sinc kernel of [`resample_channel`] precomputed for every phase of a rational +/// rate change `from` → `to` (reduced to `step` / `phases`). Output sample `n` reads input time +/// `n * step / phases`: whole part `i`, phase `p`. Its taps are input samples `i + first ..`, one +/// weight each, so every output is one dot product over contiguous memory instead of a table +/// lookup and interpolation per tap. The weights are the ones `resample_channel` computes, taps +/// beyond the kernel's half width weigh 0, and taps are summed in the same order. +struct Polyphase { + phases: u64, + step: u64, + /// Offset of the first tap from `i`. + first: i64, + taps: usize, + /// `phases * taps` weights, phase-major. + weights: Vec, +} + +impl Polyphase { + /// `None` when the bank would be too large (rates with no small common ratio); the caller + /// then falls back to [`resample_channel`]. + fn new(from: u32, to: u32, table: &SincTable) -> Option { + let g = gcd(from, to); + let (step, phases) = (u64::from(from / g), u64::from(to / g)); + let hw = table.half_width; + let first = (-hw).floor() as i64; + let last = hw.ceil() as i64 + 1; + let taps = usize::try_from(last - first + 1).ok()?; + if usize::try_from(phases).ok()?.checked_mul(taps)? > MAX_POLYPHASE_WEIGHTS { + return None; + } + let mut weights = Vec::with_capacity(taps * phases as usize); + for p in 0..phases { + let frac = p as f64 / phases as f64; + for k in first..first + taps as i64 { + let d = frac - k as f64; + weights.push(if d.abs() <= hw { table.at(d) } else { 0.0 }); + } + } + Some(Polyphase { phases, step, first, taps, weights }) + } + + fn run(&self, x: &[f32], out_len: usize) -> Vec { + let len = x.len() as i64; + let mut out = Vec::with_capacity(out_len); + for n in 0..out_len as u64 { + let Some(pos) = n.checked_mul(self.step) else { break }; + let (i, p) = ((pos / self.phases) as i64, (pos % self.phases) as usize); + let w = self.weights.get(p * self.taps..(p + 1) * self.taps).unwrap_or_default(); + let start = i + self.first; + let acc = match usize::try_from(start).ok().and_then(|s| x.get(s..s + self.taps)) { + // Inside the signal: a plain dot product. + Some(xs) => xs.iter().zip(w).fold(0.0f32, |a, (v, k)| a + v * k), + // At the edges, taps outside the signal are skipped, as `resample_channel` does. + None => (0..self.taps as i64) + .filter(|k| (0..len).contains(&(start + k))) + .fold(0.0f32, |a, k| a + x.get((start + k) as usize).copied().unwrap_or(0.0) * w.get(k as usize).copied().unwrap_or(0.0)), + }; + out.push(acc); + } + out.resize(out_len, 0.0); + out + } +} + +fn gcd(mut a: u32, mut b: u32) -> u32 { + while b != 0 { + (a, b) = (b, a % b); + } + a.max(1) +} + /// High-quality windowed-sinc sample-rate conversion from `from` Hz to `to` Hz. pub fn resample(ch: &[Vec], from: u32, to: u32) -> Vec> { if from == 0 || to == 0 || from == to { return ch.to_vec(); } let ratio = f64::from(to) / f64::from(from); + let cutoff = if ratio < 1.0 { ratio * 0.94 } else { 0.97 }; + let table = SincTable::new(cutoff); + let bank = Polyphase::new(from, to, &table); ch.iter() .map(|c| { let out_len = ((c.len() as f64 * ratio).round() as usize).min(MAX_OFFLINE_SAMPLES); - resample_ratio(std::slice::from_ref(c), ratio, out_len).into_iter().next().unwrap_or_default() + match &bank { + Some(b) => b.run(c, out_len), + None => resample_channel(c, ratio, out_len, &table), + } }) .collect() } @@ -652,6 +731,36 @@ pub fn peak_db(ch: &[Vec]) -> f32 { mod tests { use super::*; + /// The polyphase bank gives what the per-tap kernel gives (it only moves where the weights + /// are computed), at common rate pairs, both directions, and at the signal's edges. + #[test] + fn polyphase_resample_matches_the_per_tap_kernel() { + let x: Vec = (0..3_000).map(|i| (i as f32 * 0.071).sin() * 0.6 + (i as f32 * 0.53).sin() * 0.3).collect(); + let pairs = [(44_100, 48_000), (48_000, 44_100), (48_000, 96_000), (96_000, 48_000), (22_050, 48_000), (8_000, 44_100), (48_000, 22_050)]; + for (from, to) in pairs { + let ratio = f64::from(to) / f64::from(from); + let out_len = ((x.len() as f64 * ratio).round() as usize).min(MAX_OFFLINE_SAMPLES); + let cutoff = if ratio < 1.0 { ratio * 0.94 } else { 0.97 }; + let table = SincTable::new(cutoff); + let reference = resample_channel(&x, ratio, out_len, &table); + let bank = Polyphase::new(from, to, &table).expect("small common ratio"); + let fast = bank.run(&x, out_len); + assert_eq!(fast.len(), reference.len(), "{from}->{to}"); + let worst = fast.iter().zip(&reference).map(|(a, b)| (a - b).abs()).fold(0.0f32, f32::max); + assert!(worst < 1e-5, "{from}->{to}: worst difference {worst}"); + assert_eq!(resample(std::slice::from_ref(&x), from, to)[0], fast, "resample uses the bank"); + } + // Signals shorter than the kernel are all edge. + let short = [0.5f32, -0.25, 0.75]; + let table = SincTable::new(0.97); + let reference = resample_channel(&short, 48_000.0 / 44_100.0, 3, &table); + let fast = Polyphase::new(44_100, 48_000, &table).expect("bank").run(&short, 3); + assert!(fast.iter().zip(&reference).all(|(a, b)| (a - b).abs() < 1e-6), "{fast:?} {reference:?}"); + // Rates with no small common ratio fall back to the per-tap kernel. + assert!(Polyphase::new(44_100, 96_001, &table).is_none()); + assert_eq!(resample(&[short.to_vec()], 44_100, 96_001)[0].len(), 7); + } + #[test] fn fade_shapes_endpoints() { for s in FadeShape::ALL {