From 5a4dc0e2aab21d49baf3cf7aa4d119a7d957ec31 Mon Sep 17 00:00:00 2001 From: Weird Constructor Date: Wed, 20 Jul 2022 03:57:33 +0200 Subject: [PATCH] Tested the cubic interpolation now properly, also against alternative maths --- src/dsp/helpers.rs | 48 ++++++++++++++---- tests/delay_buffer.rs | 115 ++++++++++++++++++++++++++++++++++-------- 2 files changed, 132 insertions(+), 31 deletions(-) diff --git a/src/dsp/helpers.rs b/src/dsp/helpers.rs index b8526fd..5df3b93 100644 --- a/src/dsp/helpers.rs +++ b/src/dsp/helpers.rs @@ -901,14 +901,14 @@ fn fclampc(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(data: &[F], len: usize, index: usize, fract: F) -> F { @@ -931,7 +931,31 @@ pub fn cubic_interpolate(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::(0.5) * fract * ( + // x1 - xm1 + fract * ( + // f::(4.0) * x1 + // + f::(2.0) * xm1 + // - f::(5.0) * x0 + // - x2 + // + fract * (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,14 +1091,18 @@ impl DelayBuffer { 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::(1.0) - fract); -// eprintln!( -// "cubic at={} ({:6.4}) res={:6.4}", -// i % len, -// s_offs.to_f64().unwrap(), -// res.to_f64().unwrap() -// ); + // eprintln!( + // "cubic at={} ({:6.4}) res={:6.4}", + // i % len, + // s_offs.to_f64().unwrap(), + // res.to_f64().unwrap() + // ); res } diff --git a/tests/delay_buffer.rs b/tests/delay_buffer.rs index c305ada..d27a1d5 100644 --- a/tests/delay_buffer.rs +++ b/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; - 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; + let pos_inc = 0.1; 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, - 0.73999995, 0.71999997, 0.6999999, 0.67999995, 0.65999997, 0.6399999, 0.61999995, - 0.59999996, 0.58, 0.56, 0.54, 0.52000004, 0.50000006, 0.48000008, 0.4600001, - 0.44000012, 0.42000014 + 1.0, 1.03455, 1.0504, 1.05085, 1.0392, 1.01875, 0.99279994, 0.9646499, 0.9375999, + 0.91494995, 0.9, 0.89, 0.87999994, 0.86999995, 0.85999995, 0.84999996, 0.84, 0.83, + 0.82, 0.80999994, 0.8, 0.79, 0.78000003, 0.77000004, 0.76, 0.75, 0.74, 0.73, 0.72, + 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 ] ); }