fix: aliasing in sinc downsampling

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

dasp_interpolate/src/lib.rs

@@ -23,7 +23,7 @@
 //! To enable all of the above features in a `no_std` context, enable the **all-no-std** feature.
 
 #![cfg_attr(not(feature = "std"), no_std)]
 #![cfg_attr(not(feature = "std"), feature(core_intrinsics))]
 
 use dasp_frame::Frame;
 
@@ -53,4 +53,9 @@
     ///
     /// 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) {}
 }

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<S> {
     frames: ring_buffer::Fixed<S>,
     idx: usize,
+    bandwidth: f64,
 }
 
 impl<S> Sinc<S> {
-    /// 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<S> Sinc<S> {
     {
         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,
+        <S::Element as Frame>::Sample: Duplex<f64>,
+    {
+        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::<f64>())
+                    .to_sample::<<S::Element as Frame>::Sample>()
+                    .to_signed_sample(),
+            )
+        })
+    }
 }
 
 impl<S> Interpolator for Sinc<S>
@@ -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::<f64>())
-                            .to_sample::<<Self::Frame as Frame>::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::<f64>())
-                        .to_sample::<<Self::Frame as Frame>::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);
+    }
 }

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);
     }
 

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<f64> = 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
+    );
+}