Tested the cubic interpolation now properly, also against alternative maths

src/dsp/helpers.rs

@@ -901,14 +901,14 @@ fn fclampc<F: Flt>(x: F, mi: f64, mx: f64) -> F {
 ///```
 /// use hexodsp::dsp::helpers::cubic_interpolate;
 ///
-/// let buf [f32; 9] = [1.0, 0.2, 0.4, 0.5, 0.7, 0.9, 1.0, 0.3, 0.3];
+/// let buf : [f32; 9] = [1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2];
 /// let pos = 3.3_f32;
 ///
-/// let i = pos.floor();
+/// let i = pos.floor() as usize;
 /// let f = pos.fract();
 ///
 /// let res = cubic_interpolate(&buf[..], buf.len(), i, f);
-/// assert_eq!((res - 0.4).abs() < 0.2);
+/// assert!((res - 0.67).abs() < 0.2_f32);
 ///```
 #[inline]
 pub fn cubic_interpolate<F: Flt>(data: &[F], len: usize, index: usize, fract: F) -> F {
@@ -931,7 +931,31 @@ pub fn cubic_interpolate<F: Flt>(data: &[F], len: usize, index: usize, fract: F)
     let a = w + v + (x2 - x0) * f(0.5);
     let b_neg = w + a;
 
-    (((a * fract) - b_neg) * fract + c) * fract + x0
+    let res = (((a * fract) - b_neg) * fract + c) * fract + x0;
+
+    // let rr2 =
+    //     x0 + f::<F>(0.5) * fract * (
+    //         x1 - xm1 + fract * (
+    //             f::<F>(4.0) * x1
+    //             + f::<F>(2.0) * xm1
+    //             - f::<F>(5.0) * x0
+    //             - x2
+    //             + fract * (f::<F>(3.0) * (x0 - x1) - xm1 + x2)));
+
+    // eprintln!(
+    //     "index={} fract={:6.4} xm1={:6.4} x0={:6.4} x1={:6.4} x2={:6.4} = {:6.4} <> {:6.4}",
+    //     index, fract.to_f64().unwrap(), xm1.to_f64().unwrap(), x0.to_f64().unwrap(), x1.to_f64().unwrap(), x2.to_f64().unwrap(),
+    //     res.to_f64().unwrap(),
+    //     rr2.to_f64().unwrap()
+    // );
+
+    // eprintln!(
+    //     "index={} fract={:6.4} xm1={:6.4} x0={:6.4} x1={:6.4} x2={:6.4} = {:6.4}",
+    //     index, fract.to_f64().unwrap(), xm1.to_f64().unwrap(), x0.to_f64().unwrap(), x1.to_f64().unwrap(), x2.to_f64().unwrap(),
+    //     res.to_f64().unwrap(),
+    // );
+
+    res
 }
 
 #[derive(Debug, Clone, Default)]
@@ -1067,6 +1091,10 @@ impl<F: Flt> DelayBuffer<F> {
         let offs = s_offs.floor().to_usize().unwrap_or(0) % len;
         let fract = s_offs.fract();
 
+        // (offs + 1) offset for compensating that self.wr points to the next
+        // unwritten position.
+        // Additional (offs + 1 + 1) offset for cubic_interpolate, which
+        // interpolates into the past through the delay buffer.
         let i = (self.wr + len) - (offs + 2);
         let res = cubic_interpolate(data, len, i, f::<F>(1.0) - fract);
         //        eprintln!(

tests/delay_buffer.rs

@@ -110,6 +110,51 @@ fn check_cubic_interpolate() {
     use crate::helpers::cubic_interpolate;
     let data = [1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0];
 
+    let mut samples_out = vec![];
+    let mut pos = 0.0_f32;
+    let pos_inc = 0.5_f32;
+    for _ in 0..30 {
+        let i = pos.floor() as usize;
+        let f = pos.fract();
+        samples_out.push(cubic_interpolate(&data[..], data.len(), i, f));
+        pos += pos_inc;
+    }
+    assert_vec_feq!(
+        samples_out,
+        vec![
+            1.0,
+            1.01875,
+            0.9,
+            0.85,
+            0.8,
+            0.75,
+            0.7,
+            0.65,
+            0.6,
+            0.55,
+            0.5,
+            0.45,
+            0.4,
+            0.35000002,
+            0.3,
+            0.25,
+            0.2,
+            0.15,
+            0.1,
+            -0.018750004,
+            0.0,
+            0.49999997,
+            1.0,
+            1.01875,
+            0.9,
+            0.85,
+            0.8,
+            0.75,
+            0.7,
+            0.65
+        ]
+    );
+
     let mut samples_out = vec![];
     let mut pos = 0.0_f32;
     let pos_inc = 0.1_f32;
@@ -148,23 +193,7 @@ fn check_delaybuffer_cubic_interpolation() {
 
     let mut samples_out = vec![];
     let mut pos = 0.0;
-    let pos_inc = 0.5;
+    let pos_inc = 0.1;
-    for _ in 0..20 {
-        samples_out.push(buf.cubic_interpolate_at_s(pos));
-        pos += pos_inc;
-    }
-
-    assert_vec_feq!(
-        samples_out,
-        vec![
-            1.0, 0.95, 0.9, 0.85, 0.8, 0.75, 0.7, 0.65, 0.6, 0.55, 0.5, 0.45, 0.4, 0.35000002, 0.3,
-            0.25, 0.2, 0.15, 0.1, 0.05
-        ]
-    );
-
-    let mut samples_out = vec![];
-    let mut pos = 0.0;
-    let pos_inc = 0.2;
     for _ in 0..30 {
         samples_out.push(buf.cubic_interpolate_at_s(pos));
         pos += pos_inc;
@@ -173,10 +202,54 @@ fn check_delaybuffer_cubic_interpolation() {
     assert_vec_feq!(
         samples_out,
         vec![
-            1.0, 0.98, 0.96, 0.94, 0.91999996, 0.9, 0.88, 0.85999995, 0.84, 0.82, 0.8, 0.78, 0.76,
+            1.0, 1.03455, 1.0504, 1.05085, 1.0392, 1.01875, 0.99279994, 0.9646499, 0.9375999,
-            0.73999995, 0.71999997, 0.6999999, 0.67999995, 0.65999997, 0.6399999, 0.61999995,
+            0.91494995, 0.9, 0.89, 0.87999994, 0.86999995, 0.85999995, 0.84999996, 0.84, 0.83,
-            0.59999996, 0.58, 0.56, 0.54, 0.52000004, 0.50000006, 0.48000008, 0.4600001,
+            0.82, 0.80999994, 0.8, 0.79, 0.78000003, 0.77000004, 0.76, 0.75, 0.74, 0.73, 0.72,
-            0.44000012, 0.42000014
+            0.71000004
+        ]
+    );
+
+    let mut samples_out = vec![];
+    let mut pos = 0.0;
+    let pos_inc = 0.5;
+    for _ in 0..30 {
+        samples_out.push(buf.cubic_interpolate_at_s(pos));
+        pos += pos_inc;
+    }
+
+    assert_vec_feq!(
+        samples_out,
+        vec![
+            1.0,
+            1.01875,
+            0.9,
+            0.85,
+            0.8,
+            0.75,
+            0.7,
+            0.65,
+            0.6,
+            0.55,
+            0.5,
+            0.45,
+            0.4,
+            0.35000002,
+            0.3,
+            0.25,
+            0.2,
+            0.15,
+            0.1,
+            0.043750003,
+            0.0,
+            -0.00625,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0
         ]
     );
 }