From 5e414e8eb4b4658615b0b77d28a70ee55a32d3a8 Mon Sep 17 00:00:00 2001 From: Rob Morgan Date: Fri, 7 Aug 2026 09:02:48 +0800 Subject: [PATCH] fix: band-limit batch resampling (ROADMAP Stage 17) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit resample_sinc now scales its kernel cutoff when downsampling — stopband at the OUTPUT Nyquist with the same margin policy as the streaming kernel (stopband edge at the fold, not the -6 dB point; a cutoff exactly at the fold left near-Nyquist content in the transition band at ~-10 dB). Tap span widens to keep the zero-crossing count under the dilated kernel, and the Kaiser window moved to a per-call lookup table so bessel_i0 leaves the per-tap loop. AudioBuffer::resample switches from unfiltered cubic to the band-limited sinc. Measured: 18 kHz tone downsampled 2:1 folded at essentially full level before (1.9 dB rejection) and measures 89.8 dB down after; a 23 kHz tone through 48 -> 44.1 kHz conversion no longer images into the audible band. Honest scope note: pitch_shift itself measured 69.3 dB SFDR before AND after on in-band content - its pipeline only downsamples on pitch-up, where foldable content maps above Nyquist anyway - so the audible beneficiaries are the public resample APIs (review finding D5's practical weight was there, not in the shifter). New gates: downsample alias rejection >= 60 dB (resample.rs), 48->44.1 ultrasonic-fold bound (types.rs). Full suites, clippy, docs, desktop green. Co-Authored-By: Claude Fable 5 --- ROADMAP.md | 13 ++++++ src/core/resample.rs | 108 +++++++++++++++++++++++++++++++++++-------- src/core/types.rs | 38 ++++++++++++++- 3 files changed, 139 insertions(+), 20 deletions(-) diff --git a/ROADMAP.md b/ROADMAP.md index 9432a37..15f05d0 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -353,6 +353,19 @@ dependency) is complete — the fixed PV is available to audition. Automation: auto +> **Status (2026-08-07): implementation landed** (branch +> `fix/batch-resampler-antialiasing`): `resample_sinc` cutoff-scales when +> downsampling with the streaming path's stopband margin, and the Kaiser +> window moved to a per-call lookup table (Bessel out of the tap loop); +> `AudioBuffer::resample` switched from unfiltered cubic to the +> band-limited sinc. Measured: 2:1 downsample alias rejection 1.9 → +> 89.8 dB; 48→44.1 ultrasonic fold gated. Honest scope note: pitch_shift +> itself measured 69.3 dB SFDR before AND after on in-band content — its +> pipeline only downsamples on pitch-up, where foldable content maps out +> of band anyway, so the audible beneficiaries are the public resample +> APIs, not the shifter. Remaining: none (gates in CI via the standard +> suite); the stage closes on merge. + ### Why `pitch_shift()` downsamples through `resample_sinc`, whose cutoff never diff --git a/src/core/resample.rs b/src/core/resample.rs index 07f5342..7d0c9ad 100644 --- a/src/core/resample.rs +++ b/src/core/resample.rs @@ -86,17 +86,54 @@ pub fn resample_sinc(input: &[f32], output_len: usize, lobes: usize) -> Vec return vec![]; } let lobes = lobes.max(1); - if input.len() < 2 * lobes { - return resample_cubic(input, output_len); - } let ratio = (input.len() - 1) as f64 / (output_len.max(1) - 1).max(1) as f64; + // Anti-aliasing (ROADMAP Stage 17): when downsampling, scale the + // kernel cutoff so the stopband lands at the OUTPUT Nyquist instead + // of the input's, and widen the tap span to keep the same number of + // zero crossings under the dilated kernel. Without this the kernel + // passed the full input band and downsampling folded everything + // above the output Nyquist back into the audible range — the + // pitch-shift-up path aliased on bright material. + // Same margin policy as the streaming kernel (`cutoff_for_step`): + // shrink the passband a further [`STREAM_SINC_CUTOFF_SCALE`] (ramped + // in just past unity) so the finite kernel's STOPBAND — not its + // -6 dB point — lands at the fold frequency; a cutoff exactly at the + // fold leaves near-Nyquist content in the transition band at ~-10 dB. + let cutoff = if ratio > 1.0 { + let t = ((ratio - 1.0) / (STREAM_SINC_CUTOFF_RAMP_END - 1.0)).min(1.0); + let margin = 1.0 - (1.0 - STREAM_SINC_CUTOFF_SCALE) * t; + margin / ratio + } else { + 1.0 + }; + let half_span = (lobes as f64 / cutoff).ceil() as isize; + if input.len() < 2 * half_span as usize { + return resample_cubic(input, output_len); + } let mut output = Vec::with_capacity(output_len); - // Pre-compute Kaiser window for the sinc kernel. - // Beta = 6.0 gives ~60 dB stopband attenuation, good for audio. + // Kaiser window (beta = 6.0, ~60 dB stopband), sampled into a lookup + // table once per call so `bessel_i0` stays out of the per-tap loop. let beta = 6.0f64; let bessel_beta = bessel_i0(beta); + const WINDOW_TABLE: usize = 2_048; + let window_table: Vec = (0..=WINDOW_TABLE) + .map(|k| { + let t = k as f64 / WINDOW_TABLE as f64; + bessel_i0(beta * (1.0 - t * t).max(0.0).sqrt()) / bessel_beta + }) + .collect(); + let window_at = |t: f64| -> f64 { + let t = t.abs(); + if t >= 1.0 { + return 0.0; + } + let pos = t * WINDOW_TABLE as f64; + let k = pos as usize; + let frac = pos - k as f64; + window_table[k] * (1.0 - frac) + window_table[k + 1] * frac + }; for i in 0..output_len { let pos = i as f64 * ratio; @@ -106,37 +143,31 @@ pub fn resample_sinc(input: &[f32], output_len: usize, lobes: usize) -> Vec let mut sample = 0.0f64; let mut weight_sum = 0.0f64; - // Convolve with the windowed sinc kernel - let start = -(lobes as isize) + 1; - let end = lobes as isize + 1; + // Convolve with the (cutoff-scaled) windowed sinc kernel. + let start = -half_span + 1; + let end = half_span + 1; for j in start..end { let idx = center + j; if idx < 0 || idx >= input.len() as isize { continue; } - let x = frac - j as f64; + let x = (frac - j as f64) * cutoff; let sinc_val = if x.abs() < 1e-10 { 1.0 } else { let pi_x = std::f64::consts::PI * x; pi_x.sin() / pi_x }; - - // Kaiser window - let t = (j as f64 - frac) / lobes as f64; - let window = if t.abs() <= 1.0 { - bessel_i0(beta * (1.0 - t * t).max(0.0).sqrt()) / bessel_beta - } else { - 0.0 - }; + let window = window_at((j as f64 - frac) * cutoff / lobes as f64); let w = sinc_val * window; sample += input[idx as usize] as f64 * w; weight_sum += w; } - // Normalize to preserve DC gain + // Normalize to preserve DC gain (also absorbs the cutoff's + // kernel-gain factor). if weight_sum.abs() > 1e-10 { sample /= weight_sum; } @@ -1152,6 +1183,47 @@ mod tests { )); } + /// ROADMAP Stage 17: downsampling must band-limit. An 18 kHz tone + /// downsampled 2:1 folds to 4050 Hz at the output rate; before the + /// cutoff scaling it arrived at essentially full level (1.9 dB + /// rejection), after it measures ~90 dB down. + #[test] + fn test_resample_sinc_antialiases_downsampling() { + let sr = 44_100.0f64; + let n = 44_100usize; + let goertzel = |seg: &[f32], rate: f64, freq: f64| -> f64 { + let w = 2.0 * std::f64::consts::PI * freq / rate; + let coeff = 2.0 * w.cos(); + let (mut s1, mut s2) = (0.0f64, 0.0f64); + for &x in seg { + let s0 = x as f64 + coeff * s1 - s2; + s2 = s1; + s1 = s0; + } + ((s1 * s1 + s2 * s2 - coeff * s1 * s2).max(0.0)).sqrt() / (seg.len() as f64 / 2.0) + }; + let tone: Vec = (0..n) + .map(|i| (2.0 * std::f64::consts::PI * 18_000.0 * i as f64 / sr).sin() as f32) + .collect(); + let down = resample_sinc_default(&tone, n / 2); + let alias = goertzel(&down[2_000..down.len() - 2_000], sr / 2.0, 4_050.0); + let reference: Vec = (0..n) + .map(|i| (2.0 * std::f64::consts::PI * 5_000.0 * i as f64 / sr).sin() as f32) + .collect(); + let ref_down = resample_sinc_default(&reference, n / 2); + let passband = goertzel(&ref_down[2_000..ref_down.len() - 2_000], sr / 2.0, 5_000.0); + let rejection_db = 20.0 * (passband / alias.max(1e-12)).log10(); + assert!( + rejection_db > 60.0, + "downsample alias rejection {rejection_db:.1} dB (measured ~90 dB \ + with cutoff scaling, 1.9 dB without)" + ); + assert!( + passband > 0.9, + "passband level dropped through the downsample: {passband:.3}" + ); + } + #[test] fn test_bessel_i0_known_values() { // I0(0) = 1.0 diff --git a/src/core/types.rs b/src/core/types.rs index 2d44f34..019d035 100644 --- a/src/core/types.rs +++ b/src/core/types.rs @@ -624,11 +624,15 @@ impl AudioBuffer { }; } - // Resample each channel independently using cubic interpolation + // Resample each channel independently with the band-limited + // windowed-sinc converter (ROADMAP Stage 17): cubic had no + // anti-aliasing at all, so 48 -> 44.1 kHz conversions folded + // ultrasonic content into the audible band. let mut output = Vec::with_capacity(target_frames * nc); for ch in 0..nc { let channel_data: Vec = self.data.iter().skip(ch).step_by(nc).copied().collect(); - let resampled = crate::core::resample::resample_cubic(&channel_data, target_frames); + let resampled = + crate::core::resample::resample_sinc_default(&channel_data, target_frames); // Interleave: store in scratch, will interleave below if ch == 0 { output.resize(target_frames * nc, 0.0); @@ -2569,6 +2573,36 @@ mod tests { // --- resample tests --- + /// ROADMAP Stage 17: 48 -> 44.1 kHz conversion must reject content + /// above the target Nyquist instead of folding it (cubic folded it at + /// full level). A 23 kHz tone at 48 kHz folds to ~21.1 kHz at 44.1. + #[test] + fn test_resample_rejects_above_target_nyquist() { + let sr_in = 48_000.0f64; + let n = 48_000usize; + let data: Vec = (0..n) + .map(|i| (2.0 * std::f64::consts::PI * 23_000.0 * i as f64 / sr_in).sin() as f32) + .collect(); + let buf = AudioBuffer::new(data, 48_000, Channels::Mono); + let out = buf.resample(44_100); + let seg = &out.data[2_000..out.data.len() - 2_000]; + let fold_hz = 44_100.0 - 23_000.0; // 21.1 kHz image + let w = 2.0 * std::f64::consts::PI * fold_hz / 44_100.0; + let coeff = 2.0 * w.cos(); + let (mut s1, mut s2) = (0.0f64, 0.0f64); + for &x in seg { + let s0 = x as f64 + coeff * s1 - s2; + s2 = s1; + s1 = s0; + } + let alias = + ((s1 * s1 + s2 * s2 - coeff * s1 * s2).max(0.0)).sqrt() / (seg.len() as f64 / 2.0); + assert!( + alias < 0.03, + "48->44.1 folds ultrasonic content into the audible band: image level {alias:.4}" + ); + } + #[test] fn test_resample_same_rate() { let buf = AudioBuffer::from_mono(vec![1.0, 2.0, 3.0], 44100);