fix: scale sinc kernel by bandwidth

Author: roderickvd

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::<f64>())
@@ -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::<f64>())

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<f64> = downsampled.take(1500).collect();
-
-    let rms = (samples.iter().map(|&s| s * s).sum::<f64>() / 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
     );
 }