From 05cb4a9210ebb58d7e8aefdf88d0df1d88557b7f Mon Sep 17 00:00:00 2001 From: Roderick van Domburg Date: Sun, 11 Jan 2026 00:07:31 +0100 Subject: [PATCH 1/4] fix: aliasing in sinc downsampling Add rate-dependent configuration to Interpolator and use it in Sinc to apply an anti-aliasing bandwidth when downsampling. The Interpolator trait now exposes set_hz_to_hz, set_playback_hz_scale, and set_sample_hz_scale (default no-ops); Converter methods call them automatically. Fixes #174 --- CHANGELOG.md | 5 +++++ dasp_interpolate/src/lib.rs | 20 ++++++++++++++++++++ dasp_interpolate/src/sinc/mod.rs | 25 ++++++++++++++++++++----- dasp_signal/src/interpolate.rs | 14 +++++++++++++- dasp_signal/tests/interpolate.rs | 27 +++++++++++++++++++++++++++ 5 files changed, 85 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5cb66e2d..f76b93c8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,10 @@ # Unreleased +- Fixed aliasing in Sinc interpolator when downsampling by adding automatic + anti-aliasing filter cutoff adjustment. The `Interpolator` trait now includes + `set_hz_to_hz`, `set_playback_hz_scale`, and `set_sample_hz_scale` methods + (with default no-op implementations) that are called automatically by the + corresponding `Converter` methods to configure rate-dependent parameters. - Renamed `window-hanning` to `window-hann` - Made `IntoInterleavedSamples` and `IntoInterleavedSamplesIterator` stop yielding samples when the underlying signal gets exhausted. This is a breaking diff --git a/dasp_interpolate/src/lib.rs b/dasp_interpolate/src/lib.rs index cf3daf2e..97c51ae1 100644 --- a/dasp_interpolate/src/lib.rs +++ b/dasp_interpolate/src/lib.rs @@ -38,6 +38,17 @@ pub mod sinc; /// /// Implementations should keep track of the necessary data both before and after the current /// frame. +/// +/// # Rate Configuration +/// +/// Some interpolators require sample rate information to operate correctly (e.g., sinc +/// interpolation needs the rate ratio for anti-aliasing). The `set_hz_to_hz`, +/// `set_playback_hz_scale`, and `set_sample_hz_scale` methods provide alternative ways +/// to configure this - use whichever matches the information available. These methods +/// are called automatically by the corresponding `Converter` methods. +/// +/// Interpolators that don't need rate information (floor, linear) can use the default +/// no-op implementations. pub trait Interpolator { /// The type of frame over which the interpolate may operate. type Frame: Frame; @@ -53,4 +64,13 @@ pub trait Interpolator { /// /// Call this when there's a break in the continuity of the input data stream. fn reset(&mut self); + + /// Configures the interpolator from absolute sample rates. + fn set_hz_to_hz(&mut self, _source_hz: f64, _target_hz: f64) {} + + /// Configures the interpolator from playback rate scale (`source_hz / target_hz`). + fn set_playback_hz_scale(&mut self, _scale: f64) {} + + /// Configures the interpolator from sample rate scale (`target_hz / source_hz`). + fn set_sample_hz_scale(&mut self, _scale: f64) {} } diff --git a/dasp_interpolate/src/sinc/mod.rs b/dasp_interpolate/src/sinc/mod.rs index dd6d6029..5598ea9e 100644 --- a/dasp_interpolate/src/sinc/mod.rs +++ b/dasp_interpolate/src/sinc/mod.rs @@ -26,6 +26,7 @@ mod ops; pub struct Sinc { frames: ring_buffer::Fixed, idx: usize, + bandwidth: f64, } impl Sinc { @@ -49,8 +50,9 @@ impl Sinc { { assert!(frames.len() % 2 == 0); Sinc { - frames: frames, + frames, idx: 0, + bandwidth: 1.0, } } @@ -77,6 +79,7 @@ where let nl = self.idx; let nr = self.idx + 1; let depth = self.depth(); + let bandwidth = self.bandwidth; let rightmost = nl + depth; let leftmost = nr as isize - depth as isize; @@ -90,9 +93,9 @@ where (0..max_depth).fold(Self::Frame::EQUILIBRIUM, |mut v, n| { v = { - let a = PI * (phil + n as f64); + let a = PI * bandwidth * (phil + n as f64); let first = if a == 0.0 { 1.0 } else { sin(a) / a }; - let second = 0.5 + 0.5 * cos(a / depth as f64); + let second = 0.5 + 0.5 * cos(a / (depth as f64 * bandwidth)); v.zip_map(self.frames[nl - n], |vs, r_lag| { vs.add_amp( (first * second * r_lag.to_sample::()) @@ -102,9 +105,9 @@ where }) }; - let a = PI * (phir + n as f64); + let a = PI * bandwidth * (phir + n as f64); let first = if a == 0.0 { 1.0 } else { sin(a) / a }; - let second = 0.5 + 0.5 * cos(a / depth as f64); + let second = 0.5 + 0.5 * cos(a / (depth as f64 * bandwidth)); v.zip_map(self.frames[nr + n], |vs, r_lag| { vs.add_amp( (first * second * r_lag.to_sample::()) @@ -129,4 +132,16 @@ where *frame = Self::Frame::EQUILIBRIUM; } } + + fn set_hz_to_hz(&mut self, source_hz: f64, target_hz: f64) { + self.bandwidth = (target_hz / source_hz).min(1.0); + } + + fn set_playback_hz_scale(&mut self, scale: f64) { + self.bandwidth = (1.0 / scale).min(1.0); + } + + fn set_sample_hz_scale(&mut self, scale: f64) { + self.bandwidth = scale.min(1.0); + } } diff --git a/dasp_signal/src/interpolate.rs b/dasp_signal/src/interpolate.rs index fd28540d..65646b4e 100644 --- a/dasp_signal/src/interpolate.rs +++ b/dasp_signal/src/interpolate.rs @@ -31,8 +31,11 @@ where { /// Construct a new `Converter` from the source frames and the source and target sample rates /// (in Hz). + /// + /// This method calls `interpolator.set_hz_to_hz(source_hz, target_hz)` internally. #[inline] - pub fn from_hz_to_hz(source: S, interpolator: I, source_hz: f64, target_hz: f64) -> Self { + pub fn from_hz_to_hz(source: S, mut interpolator: I, source_hz: f64, target_hz: f64) -> Self { + interpolator.set_hz_to_hz(source_hz, target_hz); Self::scale_playback_hz(source, interpolator, source_hz / target_hz) } @@ -72,24 +75,33 @@ where /// Update the `source_to_target_ratio` internally given the source and target hz. /// /// This method might be useful for changing the sample rate during playback. + /// + /// This method calls `interpolator.set_hz_to_hz(source_hz, target_hz)` internally. #[inline] pub fn set_hz_to_hz(&mut self, source_hz: f64, target_hz: f64) { + self.interpolator.set_hz_to_hz(source_hz, target_hz); self.set_playback_hz_scale(source_hz / target_hz) } /// Update the `source_to_target_ratio` internally given a new **playback rate** multiplier. /// /// This method is useful for dynamically changing rates. + /// + /// This method calls `interpolator.set_playback_hz_scale(scale)` internally. #[inline] pub fn set_playback_hz_scale(&mut self, scale: f64) { + self.interpolator.set_playback_hz_scale(scale); self.source_to_target_ratio = scale; } /// Update the `source_to_target_ratio` internally given a new **sample rate** multiplier. /// /// This method is useful for dynamically changing rates. + /// + /// This method calls `interpolator.set_sample_hz_scale(scale)` internally. #[inline] pub fn set_sample_hz_scale(&mut self, scale: f64) { + self.interpolator.set_sample_hz_scale(scale); self.set_playback_hz_scale(1.0 / scale); } diff --git a/dasp_signal/tests/interpolate.rs b/dasp_signal/tests/interpolate.rs index 065ad74d..ac5fd16f 100644 --- a/dasp_signal/tests/interpolate.rs +++ b/dasp_signal/tests/interpolate.rs @@ -71,3 +71,30 @@ fn test_sinc() { None ); } + +#[test] +fn test_sinc_downsampling_antialiasing() { + const SOURCE_HZ: f64 = 2000.0; + const TARGET_HZ: f64 = 1500.0; + const SIGNAL_FREQ: f64 = 900.0; + + let source_signal = signal::rate(SOURCE_HZ).const_hz(SIGNAL_FREQ).sine(); + + let ring_buffer = ring_buffer::Fixed::from(vec![0.0; 100]); + let sinc = Sinc::new(ring_buffer); + let mut downsampled = source_signal.from_hz_to_hz(sinc, SOURCE_HZ, TARGET_HZ); + + for _ in 0..50 { + downsampled.next(); + } + + let samples: Vec = downsampled.take(1500).collect(); + + let rms = (samples.iter().map(|&s| s * s).sum::() / samples.len() as f64).sqrt(); + + assert!( + rms < 0.1, + "Expected RMS < 0.1 (well-filtered), got {:.3}. Aliasing occurred.", + rms + ); +} From 9668870b9062b51119e83cc4c68dbac9e854d740 Mon Sep 17 00:00:00 2001 From: Roderick van Domburg Date: Sun, 11 Jan 2026 13:43:42 +0100 Subject: [PATCH 2/4] fix: scale sinc kernel by bandwidth --- dasp_interpolate/src/sinc/mod.rs | 24 ++++++++++++++++++------ dasp_signal/tests/interpolate.rs | 31 ++++++++++++++++++++++--------- 2 files changed, 40 insertions(+), 15 deletions(-) diff --git a/dasp_interpolate/src/sinc/mod.rs b/dasp_interpolate/src/sinc/mod.rs index 5598ea9e..9fee5b5b 100644 --- a/dasp_interpolate/src/sinc/mod.rs +++ b/dasp_interpolate/src/sinc/mod.rs @@ -93,9 +93,15 @@ where (0..max_depth).fold(Self::Frame::EQUILIBRIUM, |mut v, n| { v = { - let a = PI * bandwidth * (phil + n as f64); - let first = if a == 0.0 { 1.0 } else { sin(a) / a }; - let second = 0.5 + 0.5 * cos(a / (depth as f64 * bandwidth)); + let t = phil + n as f64; + let a = PI * bandwidth * t; + let b = PI * t / depth as f64; + let first = if a.abs() < f64::EPSILON { + bandwidth + } else { + bandwidth * sin(a) / a + }; + let second = 0.5 + 0.5 * cos(b); v.zip_map(self.frames[nl - n], |vs, r_lag| { vs.add_amp( (first * second * r_lag.to_sample::()) @@ -105,9 +111,15 @@ where }) }; - let a = PI * bandwidth * (phir + n as f64); - let first = if a == 0.0 { 1.0 } else { sin(a) / a }; - let second = 0.5 + 0.5 * cos(a / (depth as f64 * bandwidth)); + let t = phir + n as f64; + let a = PI * bandwidth * t; + let b = PI * t / depth as f64; + let first = if a.abs() < f64::EPSILON { + bandwidth + } else { + bandwidth * sin(a) / a + }; + let second = 0.5 + 0.5 * cos(b); v.zip_map(self.frames[nr + n], |vs, r_lag| { vs.add_amp( (first * second * r_lag.to_sample::()) diff --git a/dasp_signal/tests/interpolate.rs b/dasp_signal/tests/interpolate.rs index ac5fd16f..9e2931d0 100644 --- a/dasp_signal/tests/interpolate.rs +++ b/dasp_signal/tests/interpolate.rs @@ -73,28 +73,41 @@ fn test_sinc() { } #[test] -fn test_sinc_downsampling_antialiasing() { +fn test_sinc_antialiasing() { const SOURCE_HZ: f64 = 2000.0; const TARGET_HZ: f64 = 1500.0; const SIGNAL_FREQ: f64 = 900.0; + const TAPS: usize = 100; let source_signal = signal::rate(SOURCE_HZ).const_hz(SIGNAL_FREQ).sine(); - - let ring_buffer = ring_buffer::Fixed::from(vec![0.0; 100]); + let ring_buffer = ring_buffer::Fixed::from(vec![0.0; TAPS]); let sinc = Sinc::new(ring_buffer); let mut downsampled = source_signal.from_hz_to_hz(sinc, SOURCE_HZ, TARGET_HZ); - for _ in 0..50 { + // Warm up the filter to reach steady state + for _ in 0..TAPS / 2 { downsampled.next(); } + // Measure peak amplitude in steady state let samples: Vec = downsampled.take(1500).collect(); - - let rms = (samples.iter().map(|&s| s * s).sum::() / samples.len() as f64).sqrt(); + let peak_amplitude = samples + .iter() + .map(|&s| s.abs()) + .max_by(|a, b| a.partial_cmp(b).unwrap()) + .unwrap(); + + // With Hann-windowed sinc, With 50 taps and 900 Hz being relatively close to the 750 Hz + // cutoff, we're in the early stopband region where attenuation is less than -44 dB. + // Empirical measurement shows ~41 dB for this configuration. + const EXPECTED_ATTENUATION_DB: f64 = 41.0; + let max_peak_amplitude = 10.0_f64.powf(-EXPECTED_ATTENUATION_DB / 20.0); assert!( - rms < 0.1, - "Expected RMS < 0.1 (well-filtered), got {:.3}. Aliasing occurred.", - rms + peak_amplitude < max_peak_amplitude, + "Expected ≥{:.0} dB attenuation (peak < {:.4}), got peak {:.4}", + EXPECTED_ATTENUATION_DB, + max_peak_amplitude, + peak_amplitude ); } From b94688e20dd3e1e14ab9cb13e10e63a5d80fadf5 Mon Sep 17 00:00:00 2001 From: Roderick van Domburg Date: Mon, 12 Jan 2026 22:38:16 +0100 Subject: [PATCH 3/4] refactor: replace interpolator rate methods with set_bandwidth --- CHANGELOG.md | 5 ++--- dasp_interpolate/src/lib.rs | 23 ++++------------------- dasp_interpolate/src/sinc/mod.rs | 12 ++---------- dasp_signal/src/interpolate.rs | 16 +++++----------- 4 files changed, 13 insertions(+), 43 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f76b93c8..5bda1b5a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,9 +2,8 @@ - Fixed aliasing in Sinc interpolator when downsampling by adding automatic anti-aliasing filter cutoff adjustment. The `Interpolator` trait now includes - `set_hz_to_hz`, `set_playback_hz_scale`, and `set_sample_hz_scale` methods - (with default no-op implementations) that are called automatically by the - corresponding `Converter` methods to configure rate-dependent parameters. + a `set_bandwidth` method (with default no-op implementation) that is called + automatically by `Converter`. - Renamed `window-hanning` to `window-hann` - Made `IntoInterleavedSamples` and `IntoInterleavedSamplesIterator` stop yielding samples when the underlying signal gets exhausted. This is a breaking diff --git a/dasp_interpolate/src/lib.rs b/dasp_interpolate/src/lib.rs index 97c51ae1..f691cee9 100644 --- a/dasp_interpolate/src/lib.rs +++ b/dasp_interpolate/src/lib.rs @@ -38,17 +38,6 @@ pub mod sinc; /// /// Implementations should keep track of the necessary data both before and after the current /// frame. -/// -/// # Rate Configuration -/// -/// Some interpolators require sample rate information to operate correctly (e.g., sinc -/// interpolation needs the rate ratio for anti-aliasing). The `set_hz_to_hz`, -/// `set_playback_hz_scale`, and `set_sample_hz_scale` methods provide alternative ways -/// to configure this - use whichever matches the information available. These methods -/// are called automatically by the corresponding `Converter` methods. -/// -/// Interpolators that don't need rate information (floor, linear) can use the default -/// no-op implementations. pub trait Interpolator { /// The type of frame over which the interpolate may operate. type Frame: Frame; @@ -65,12 +54,8 @@ pub trait Interpolator { /// Call this when there's a break in the continuity of the input data stream. fn reset(&mut self); - /// Configures the interpolator from absolute sample rates. - fn set_hz_to_hz(&mut self, _source_hz: f64, _target_hz: f64) {} - - /// Configures the interpolator from playback rate scale (`source_hz / target_hz`). - fn set_playback_hz_scale(&mut self, _scale: f64) {} - - /// Configures the interpolator from sample rate scale (`target_hz / source_hz`). - fn set_sample_hz_scale(&mut self, _scale: f64) {} + /// Configures filter bandwidth for anti-aliasing. No-op for non-filtering interpolators. + /// + /// For downsampling, set to `target_hz / source_hz` to prevent aliasing. + fn set_bandwidth(&mut self, _bandwidth: f64) {} } diff --git a/dasp_interpolate/src/sinc/mod.rs b/dasp_interpolate/src/sinc/mod.rs index 9fee5b5b..ae0a06a6 100644 --- a/dasp_interpolate/src/sinc/mod.rs +++ b/dasp_interpolate/src/sinc/mod.rs @@ -145,15 +145,7 @@ where } } - fn set_hz_to_hz(&mut self, source_hz: f64, target_hz: f64) { - self.bandwidth = (target_hz / source_hz).min(1.0); - } - - fn set_playback_hz_scale(&mut self, scale: f64) { - self.bandwidth = (1.0 / scale).min(1.0); - } - - fn set_sample_hz_scale(&mut self, scale: f64) { - self.bandwidth = scale.min(1.0); + fn set_bandwidth(&mut self, bandwidth: f64) { + self.bandwidth = bandwidth.clamp(f64::MIN_POSITIVE, 1.0); } } diff --git a/dasp_signal/src/interpolate.rs b/dasp_signal/src/interpolate.rs index 65646b4e..0ea020e3 100644 --- a/dasp_signal/src/interpolate.rs +++ b/dasp_signal/src/interpolate.rs @@ -32,10 +32,10 @@ where /// Construct a new `Converter` from the source frames and the source and target sample rates /// (in Hz). /// - /// This method calls `interpolator.set_hz_to_hz(source_hz, target_hz)` internally. + /// Configures the interpolator's bandwidth for anti-aliasing when downsampling. #[inline] pub fn from_hz_to_hz(source: S, mut interpolator: I, source_hz: f64, target_hz: f64) -> Self { - interpolator.set_hz_to_hz(source_hz, target_hz); + interpolator.set_bandwidth(target_hz / source_hz); Self::scale_playback_hz(source, interpolator, source_hz / target_hz) } @@ -75,33 +75,27 @@ where /// Update the `source_to_target_ratio` internally given the source and target hz. /// /// This method might be useful for changing the sample rate during playback. - /// - /// This method calls `interpolator.set_hz_to_hz(source_hz, target_hz)` internally. #[inline] pub fn set_hz_to_hz(&mut self, source_hz: f64, target_hz: f64) { - self.interpolator.set_hz_to_hz(source_hz, target_hz); + self.interpolator.set_bandwidth(target_hz / source_hz); self.set_playback_hz_scale(source_hz / target_hz) } /// Update the `source_to_target_ratio` internally given a new **playback rate** multiplier. /// /// This method is useful for dynamically changing rates. - /// - /// This method calls `interpolator.set_playback_hz_scale(scale)` internally. #[inline] pub fn set_playback_hz_scale(&mut self, scale: f64) { - self.interpolator.set_playback_hz_scale(scale); + self.interpolator.set_bandwidth(1.0 / scale); self.source_to_target_ratio = scale; } /// Update the `source_to_target_ratio` internally given a new **sample rate** multiplier. /// /// This method is useful for dynamically changing rates. - /// - /// This method calls `interpolator.set_sample_hz_scale(scale)` internally. #[inline] pub fn set_sample_hz_scale(&mut self, scale: f64) { - self.interpolator.set_sample_hz_scale(scale); + self.interpolator.set_bandwidth(scale); self.set_playback_hz_scale(1.0 / scale); } From 75110f71e46209b3a43c5f62bb10d5c2a962afde Mon Sep 17 00:00:00 2001 From: Roderick van Domburg Date: Mon, 12 Jan 2026 22:55:57 +0100 Subject: [PATCH 4/4] refactor: add windowed sinc helpers and docs --- dasp_interpolate/src/sinc/mod.rs | 96 +++++++++++++++++--------------- 1 file changed, 51 insertions(+), 45 deletions(-) diff --git a/dasp_interpolate/src/sinc/mod.rs b/dasp_interpolate/src/sinc/mod.rs index ae0a06a6..5b8b0184 100644 --- a/dasp_interpolate/src/sinc/mod.rs +++ b/dasp_interpolate/src/sinc/mod.rs @@ -1,4 +1,4 @@ -//! A sinc interpolator implementation. +//! A Hann-windowed sinc interpolator implementation. //! //! ### Required Features //! @@ -14,7 +14,7 @@ use ops::f64::{cos, sin}; mod ops; -/// Interpolator for sinc interpolation. +/// Interpolator for Hann-windowed sinc interpolation. /// /// Generally accepted as one of the better sample rate converters, although it uses significantly /// more computation. @@ -30,14 +30,15 @@ pub struct Sinc { } impl Sinc { - /// Create a new **Sinc** interpolator with the given ring buffer. + /// Create a new **Sinc** (Hann-windowed) interpolator with the given ring buffer. /// - /// The given ring buffer should have a length twice that of the desired sinc interpolation - /// `depth`. + /// The given ring buffer should have a length twice that of the desired interpolation `depth`. /// - /// The initial contents of the ring_buffer will act as padding for the interpolated signal. + /// The initial contents of the ring buffer will act as padding for the interpolated signal. /// - /// **panic!**s if the given ring buffer's length is not a multiple of `2`. + /// ### Panics + /// + /// If the given ring buffer's length is not a multiple of `2`. /// /// ### Required Features /// @@ -56,12 +57,50 @@ impl Sinc { } } - fn depth(&self) -> usize + /// Get the interpolation depth (number of taps per side). + pub fn depth(&self) -> usize where S: ring_buffer::Slice, { self.frames.len() / 2 } + + /// Compute the Hann-windowed sinc coefficient for a given time offset. + fn windowed_sinc(&self, t: f64, depth: usize, bandwidth: f64) -> f64 { + let a = PI * bandwidth * t; + let b = PI * t / depth as f64; + let sinc = if a.abs() < f64::EPSILON { + bandwidth + } else { + bandwidth * sin(a) / a + }; + let hann_window = 0.5 + 0.5 * cos(b); + sinc * hann_window + } + + /// Accumulate a single tap of the Hann-windowed sinc filter. + fn accumulate_tap( + &self, + accum: S::Element, + frame_idx: usize, + phase: f64, + depth: usize, + bandwidth: f64, + ) -> S::Element + where + S: ring_buffer::Slice, + S::Element: Frame, + ::Sample: Duplex, + { + let coeff = self.windowed_sinc(phase, depth, bandwidth); + accum.zip_map(self.frames[frame_idx], |vs, frame_sample| { + vs.add_amp( + (coeff * frame_sample.to_sample::()) + .to_sample::<::Sample>() + .to_signed_sample(), + ) + }) + } } impl Interpolator for Sinc @@ -72,7 +111,7 @@ where { type Frame = S::Element; - /// Sinc interpolation + /// Hann-windowed sinc interpolation fn interpolate(&self, x: f64) -> Self::Frame { let phil = x; let phir = 1.0 - x; @@ -91,42 +130,9 @@ where depth }; - (0..max_depth).fold(Self::Frame::EQUILIBRIUM, |mut v, n| { - v = { - let t = phil + n as f64; - let a = PI * bandwidth * t; - let b = PI * t / depth as f64; - let first = if a.abs() < f64::EPSILON { - bandwidth - } else { - bandwidth * sin(a) / a - }; - let second = 0.5 + 0.5 * cos(b); - v.zip_map(self.frames[nl - n], |vs, r_lag| { - vs.add_amp( - (first * second * r_lag.to_sample::()) - .to_sample::<::Sample>() - .to_signed_sample(), - ) - }) - }; - - let t = phir + n as f64; - let a = PI * bandwidth * t; - let b = PI * t / depth as f64; - let first = if a.abs() < f64::EPSILON { - bandwidth - } else { - bandwidth * sin(a) / a - }; - let second = 0.5 + 0.5 * cos(b); - v.zip_map(self.frames[nr + n], |vs, r_lag| { - vs.add_amp( - (first * second * r_lag.to_sample::()) - .to_sample::<::Sample>() - .to_signed_sample(), - ) - }) + (0..max_depth).fold(Self::Frame::EQUILIBRIUM, |v, n| { + let v = self.accumulate_tap(v, nl - n, phil + n as f64, depth, bandwidth); + self.accumulate_tap(v, nr + n, phir + n as f64, depth, bandwidth) }) }