diff --git a/CHANGELOG.md b/CHANGELOG.md index 5cb66e2d..5bda1b5a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,9 @@ # Unreleased +- Fixed aliasing in Sinc interpolator when downsampling by adding automatic + anti-aliasing filter cutoff adjustment. The `Interpolator` trait now includes + 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 cf3daf2e..f691cee9 100644 --- a/dasp_interpolate/src/lib.rs +++ b/dasp_interpolate/src/lib.rs @@ -53,4 +53,9 @@ pub trait Interpolator { /// /// Call this when there's a break in the continuity of the input data stream. fn reset(&mut self); + + /// 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 dd6d6029..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. @@ -26,17 +26,19 @@ mod ops; pub struct Sinc { frames: ring_buffer::Fixed, idx: usize, + bandwidth: f64, } 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 /// @@ -49,17 +51,56 @@ impl Sinc { { assert!(frames.len() % 2 == 0); Sinc { - frames: frames, + frames, idx: 0, + bandwidth: 1.0, } } - 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 @@ -70,13 +111,14 @@ 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; 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; @@ -88,30 +130,9 @@ where depth }; - (0..max_depth).fold(Self::Frame::EQUILIBRIUM, |mut v, n| { - v = { - let a = PI * (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); - 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 a = PI * (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); - 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) }) } @@ -129,4 +150,8 @@ where *frame = Self::Frame::EQUILIBRIUM; } } + + 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 fd28540d..0ea020e3 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). + /// + /// Configures the interpolator's bandwidth for anti-aliasing when downsampling. #[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_bandwidth(target_hz / source_hz); Self::scale_playback_hz(source, interpolator, source_hz / target_hz) } @@ -74,6 +77,7 @@ where /// This method might be useful for changing the sample rate during playback. #[inline] pub fn set_hz_to_hz(&mut self, source_hz: f64, target_hz: f64) { + self.interpolator.set_bandwidth(target_hz / source_hz); self.set_playback_hz_scale(source_hz / target_hz) } @@ -82,6 +86,7 @@ where /// This method is useful for dynamically changing rates. #[inline] pub fn set_playback_hz_scale(&mut self, scale: f64) { + self.interpolator.set_bandwidth(1.0 / scale); self.source_to_target_ratio = scale; } @@ -90,6 +95,7 @@ where /// This method is useful for dynamically changing rates. #[inline] pub fn set_sample_hz_scale(&mut self, scale: f64) { + self.interpolator.set_bandwidth(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..9e2931d0 100644 --- a/dasp_signal/tests/interpolate.rs +++ b/dasp_signal/tests/interpolate.rs @@ -71,3 +71,43 @@ fn test_sinc() { None ); } + +#[test] +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; TAPS]); + let sinc = Sinc::new(ring_buffer); + let mut downsampled = source_signal.from_hz_to_hz(sinc, SOURCE_HZ, TARGET_HZ); + + // 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 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!( + peak_amplitude < max_peak_amplitude, + "Expected ≥{:.0} dB attenuation (peak < {:.4}), got peak {:.4}", + EXPECTED_ATTENUATION_DB, + max_peak_amplitude, + peak_amplitude + ); +}