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

Author: roderickvd

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

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

dasp_interpolate/src/sinc/mod.rs

@@ -26,6 +26,7 @@ mod ops;
 pub struct Sinc<S> {
     frames: ring_buffer::Fixed<S>,
     idx: usize,
+    bandwidth: f64,
 }
 
 impl<S> Sinc<S> {
@@ -49,8 +50,9 @@ impl<S> Sinc<S> {
     {
         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::<f64>())
@@ -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::<f64>())
@@ -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);
+    }
 }

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

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<f64> = downsampled.take(1500).collect();
+
+    let rms = (samples.iter().map(|&s| s * s).sum::<f64>() / samples.len() as f64).sqrt();
+
+    assert!(
+        rms < 0.1,
+        "Expected RMS < 0.1 (well-filtered), got {:.3}. Aliasing occurred.",
+        rms
+    );
+}