From 1fbc6741cb7f89f5e1c7e79603009e1f7216dcc3 Mon Sep 17 00:00:00 2001 From: Weird Constructor Date: Tue, 19 Jul 2022 20:47:46 +0200 Subject: [PATCH] fixing interpolation --- src/dsp/helpers.rs | 66 +++++++++++++++++++++++++++++-------------- src/dsp/node_sampl.rs | 44 ++--------------------------- tests/node_allp.rs | 32 ++++++++++----------- tests/node_comb.rs | 43 +++++++++------------------- tests/node_delay.rs | 36 +++++++++++------------ tests/node_sampl.rs | 12 ++++---- 6 files changed, 101 insertions(+), 132 deletions(-) diff --git a/src/dsp/helpers.rs b/src/dsp/helpers.rs index 3fa7c27..4491216 100644 --- a/src/dsp/helpers.rs +++ b/src/dsp/helpers.rs @@ -892,6 +892,48 @@ fn fclampc(x: F, mi: f64, mx: f64) -> F { x.max(f(mi)).min(f(mx)) } + +/// Hermite / Cubic interpolation of a buffer full of samples at the given _index_. +/// _len_ is the buffer length to consider and wrap the index into. And _fract_ is the +/// fractional part of the index. +/// +/// Commonly used like this: +/// +///``` +/// 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 pos = 3.3_f32; +/// +/// let i = pos.floor(); +/// let f = pos.fract(); +/// +/// let res = cubic_interpolate(&buf[..], buf.len(), i, f); +/// assert_eq!((res - 0.4).abs() < 0.2); +///``` +#[inline] +pub fn cubic_interpolate(data: &[F], len: usize, index: usize, fract: F) -> F { + // Hermite interpolation, take from + // https://github.com/eric-wood/delay/blob/main/src/delay.rs#L52 + // + // Thanks go to Eric Wood! + // + // For the interpolation code: + // MIT License, Copyright (c) 2021 Eric Wood + let xm1 = data[(index - 1) % len]; + let x0 = data[index % len]; + let x1 = data[(index + 1) % len]; + let x2 = data[(index + 2) % len]; + + let c = (x1 - xm1) * f(0.5); + let v = x0 - x1; + let w = c + v; + let a = w + v + (x2 - x0) * f(0.5); + let b_neg = w + a; + + (((a * fract) - b_neg) * fract + c) * fract + x0 +} + #[derive(Debug, Clone, Default)] pub struct DelayBuffer { data: Vec, @@ -1006,7 +1048,8 @@ impl DelayBuffer { /// Fetch a sample from the delay buffer at the given offset. /// - /// * `s_offs` - Sample offset in samples. + /// * `s_offs` - Sample offset in samples into the past of the [DelayBuffer] + /// from the current write (or the "now") position. #[inline] pub fn cubic_interpolate_at_s(&self, s_offs: F) -> F { let data = &self.data[..]; @@ -1016,26 +1059,7 @@ impl DelayBuffer { let i = (self.wr + len) - offs; - // Hermite interpolation, take from - // https://github.com/eric-wood/delay/blob/main/src/delay.rs#L52 - // - // Thanks go to Eric Wood! - // - // For the interpolation code: - // MIT License, Copyright (c) 2021 Eric Wood - let xm1 = data[(i + 1) % len]; - let x0 = data[i % len]; - let x1 = data[(i - 1) % len]; - let x2 = data[(i - 2) % len]; - - let c = (x1 - xm1) * f(0.5); - let v = x0 - x1; - let w = c + v; - let a = w + v + (x2 - x0) * f(0.5); - let b_neg = w + a; - - let fract = fract as F; - (((a * fract) - b_neg) * fract + c) * fract + x0 + cubic_interpolate(data, len, i, f::(1.0) - fract) } #[inline] diff --git a/src/dsp/node_sampl.rs b/src/dsp/node_sampl.rs index 176d25c..e5d5708 100644 --- a/src/dsp/node_sampl.rs +++ b/src/dsp/node_sampl.rs @@ -2,7 +2,7 @@ // This file is a part of HexoDSP. Released under GPL-3.0-or-later. // See README.md and COPYING for details. -use super::helpers::Trigger; +use super::helpers::{Trigger, cubic_interpolate}; use crate::dsp::{at, denorm, denorm_offs, inp, out}; //, inp, denorm, denorm_v, inp_dir, at}; use crate::dsp::{DspNode, LedPhaseVals, NodeContext, NodeId, ProcBuf, SAtom}; use crate::nodes::{NodeAudioContext, NodeExecContext}; @@ -154,26 +154,7 @@ impl Sampl { let f = self.phase.fract(); self.phase = j as f64 + f + sr_factor * speed; - // Hermite interpolation, take from - // https://github.com/eric-wood/delay/blob/main/src/delay.rs#L52 - // - // Thanks go to Eric Wood! - // - // For the interpolation code: - // MIT License, Copyright (c) 2021 Eric Wood - let xm1 = sample_data[(i + 1) % sd_len]; - let x0 = sample_data[i % sd_len]; - let x1 = sample_data[(i - 1) % sd_len]; - let x2 = sample_data[(i - 2) % sd_len]; - - let c = (x1 - xm1) * 0.5; - let v = x0 - x1; - let w = c + v; - let a = w + v + (x2 - x0) * 0.5; - let b_neg = w + a; - - let f = (1.0 - f) as f32; - (((a * f) - b_neg) * f + c) * f + x0 + cubic_interpolate(&sample_data[..], sd_len, i, (1.0 - f) as f32) } #[allow(clippy::many_single_char_names)] @@ -188,26 +169,7 @@ impl Sampl { let f = self.phase.fract(); self.phase = (i % sd_len) as f64 + f + sr_factor * speed; - // Hermite interpolation, take from - // https://github.com/eric-wood/delay/blob/main/src/delay.rs#L52 - // - // Thanks go to Eric Wood! - // - // For the interpolation code: - // MIT License, Copyright (c) 2021 Eric Wood - let xm1 = sample_data[(i - 1) % sd_len]; - let x0 = sample_data[i % sd_len]; - let x1 = sample_data[(i + 1) % sd_len]; - let x2 = sample_data[(i + 2) % sd_len]; - - let c = (x1 - xm1) * 0.5; - let v = x0 - x1; - let w = c + v; - let a = w + v + (x2 - x0) * 0.5; - let b_neg = w + a; - - let f = f as f32; - (((a * f) - b_neg) * f + c) * f + x0 + cubic_interpolate(&sample_data[..], sd_len, i, f as f32) } #[allow(clippy::float_cmp)] diff --git a/tests/node_allp.rs b/tests/node_allp.rs index 0680bce..a516ff2 100644 --- a/tests/node_allp.rs +++ b/tests/node_allp.rs @@ -38,35 +38,35 @@ fn check_node_allp() { // starts with original signal * -0.7 let mut v = vec![0.7; (2.0 * 44.1_f32).ceil() as usize]; // silence for 1ms, which is the internal delay of the allpass - v.append(&mut vec![0.0; (1.0 * 44.1_f32).floor() as usize - 2]); + v.append(&mut vec![0.0; (1.0 * 44.1_f32).floor() as usize - 3]); // allpass feedback of the original signal for 2ms: // XXX: the smearing before and after the allpass is due to the // cubic interpolation! - v.append(&mut vec![-0.03748519, 0.37841395, 0.5260659]); + v.append(&mut vec![-0.01606, 0.13158, 0.54748]); v.append(&mut vec![0.51; (2.0 * 44.1_f32).ceil() as usize - 3]); // 1ms allpass silence like before: - v.append(&mut vec![0.54748523, 0.13158606, -0.016065884]); - v.append(&mut vec![0.0; (1.0 * 44.1_f32).floor() as usize - 5]); + v.append(&mut vec![0.5260659, 0.37841395, -0.03748519]); + v.append(&mut vec![0.0; (1.0 * 44.1_f32).floor() as usize - 6]); // 2ms the previous 1.0 * 0.7 fed back into the filter, // including even more smearing due to cubic interpolation: v.append(&mut vec![ - -0.0019286226, - 0.04086761, - -0.1813516, - -0.35157663, - -0.36315754, - -0.35664573, + -0.00035427228, + 0.006157537, + -0.005423375, + -0.1756484, + -0.39786762, + -0.3550714, ]); v.append(&mut vec![-0.357; (2.0 * 44.1_f32).floor() as usize - 5]); v.append(&mut vec![ - -0.3550714, - -0.39786762, - -0.1756484, - -0.005423375, - 0.006157537, - -0.00035427228, + -0.35664573, + -0.36315754, + -0.35157663, + -0.1813516, + 0.04086761, + -0.0019286226, ]); v.append(&mut vec![0.0; 10]); diff --git a/tests/node_comb.rs b/tests/node_comb.rs index b855609..f86a550 100644 --- a/tests/node_comb.rs +++ b/tests/node_comb.rs @@ -36,45 +36,28 @@ fn check_node_comb_1() { (0, 216), (11, 221), (22, 216), - (3370, 206), - (3381, 248), - (3391, 191), - (6740, 185), - (6751, 207), - (6761, 195), - (10131, 215), - (10142, 210), - (10153, 213), - (10164, 201), - (20338, 187), - (20349, 184) + (3381, 184), + (3391, 189), + (3402, 195), + (3413, 213), + (3424, 198), + (3435, 203), + (6815, 188), + (13587, 196), + (13598, 210), + (13609, 207), + (13620, 193) ] ); pset_n_wait(&mut matrix, &mut node_exec, comb_1, "time", 0.030); let fft = run_and_get_avg_fft4096_now(&mut node_exec, 180); - assert_eq!( - fft, - vec![(1001, 206), (2993, 196), (3004, 219), (3994, 197), (6998, 211), (8000, 201)] - ); + assert_eq!(fft, vec![(1001, 186), (3015, 186), (7031, 191)]); pset_n_wait(&mut matrix, &mut node_exec, comb_1, "g", 0.999); let fft = run_and_get_avg_fft4096_now(&mut node_exec, 1000); - assert_eq!( - fft, - vec![ - (0, 2003), - (11, 1015), - (991, 1078), - (1001, 1837), - (2003, 1059), - (2993, 1420), - (3004, 1775), - (3994, 1297), - (4005, 1485) - ] - ); + assert_eq!(fft, vec![(0, 2008), (11, 1017)]); } #[test] diff --git a/tests/node_delay.rs b/tests/node_delay.rs index d137a60..8e5977c 100644 --- a/tests/node_delay.rs +++ b/tests/node_delay.rs @@ -68,16 +68,16 @@ fn check_node_delay_1() { 0.0, 0.0, // delayed burst of sine for 100ms: - 0.047408286, - -0.17181452, - 0.2669317, - -0.22377986, - 0.000059626997, - 0.24652793, - -0.30384338, - 0.2087649, - -0.070256576, - 0.000003647874, + 0.039102618, + -0.16390327, + 0.27611724, + -0.2608055, + 0.060164057, + 0.20197779, + -0.28871512, + 0.21515398, + -0.081471935, + 0.0023831273, // silence afterwards: 0.0, 0.0, @@ -119,8 +119,8 @@ fn check_node_delay_2() { vec![ // 10ms smoothing time for "inp" 0.001133, // 30ms delaytime just mixing the 0.5: - 0.5, 0.5, 0.5, // the delayed smoothing ramp (10ms): - 0.951113, // the delay + input signal: + 0.5, 0.5, 0.5, // the delayed smoothing ramp (10ms): + 0.9513626, // the delay + input signal: 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 ] ); @@ -172,15 +172,15 @@ fn check_node_delay_time_mod() { let fft = run_and_get_fft4096_now(&mut node_exec, 110); // Expect a sine sweep over a // range of low frequencies: - assert_eq!(fft[0], (86, 111)); - assert_eq!(fft[5], (237, 114)); - assert_eq!(fft[10], (517, 110)); + assert_eq!(fft[0], (86, 112)); + assert_eq!(fft[5], (237, 112)); + assert_eq!(fft[10], (517, 111)); // Sweep upwards: run_for_ms(&mut node_exec, 300.0); let fft = run_and_get_fft4096_now(&mut node_exec, 122); - assert_eq!(fft[0], (2498, 122)); - assert_eq!(fft[7], (2681, 122)); + assert_eq!(fft[0], (2509, 123)); + assert_eq!(fft[7], (2821, 123)); // Sweep at mostly highest point: run_for_ms(&mut node_exec, 700.0); @@ -274,7 +274,7 @@ fn check_node_delay_fb() { let idxs_big = collect_signal_changes(&res.0[..], 50); // We expect the signal to be delayed by 20ms: - assert_eq!(idxs_big, vec![(221, 106), (442, 53)]); + assert_eq!(idxs_big, vec![(220, 106), (440, 53)]); } #[test] diff --git a/tests/node_sampl.rs b/tests/node_sampl.rs index 02ae956..66391aa 100644 --- a/tests/node_sampl.rs +++ b/tests/node_sampl.rs @@ -153,7 +153,7 @@ fn check_node_sampl_reverse() { matrix.set_param(dir_p, SAtom::setting(1)); let (rms, min, max) = run_and_get_l_rms_mimax(&mut node_exec, 50.0); - assert_rmsmima!((rms, min, max), (0.50059, -0.9997, 0.9997)); + assert_rmsmima!((rms, min, max), (0.5003373, -0.9997, 0.9997)); let fft = run_and_get_fft4096(&mut node_exec, 800, 20.0); assert_eq!(fft[0], (441, 1023)); @@ -164,11 +164,11 @@ fn check_node_sampl_reverse() { matrix.set_param(freq_p, SAtom::param(-0.1)); let fft = run_and_get_fft4096(&mut node_exec, 800, 20.0); - assert_eq!(fft[0], (215, 880)); + assert_eq!(fft[0], (215, 881)); matrix.set_param(freq_p, SAtom::param(-0.2)); let fft = run_and_get_fft4096(&mut node_exec, 800, 20.0); - assert_eq!(fft[0], (108, 986)); + assert_eq!(fft[0], (108, 987)); matrix.set_param(freq_p, SAtom::param(-0.4)); let fft = run_and_get_fft4096(&mut node_exec, 800, 20.0); @@ -176,7 +176,7 @@ fn check_node_sampl_reverse() { matrix.set_param(freq_p, SAtom::param(-0.5)); let fft = run_and_get_fft4096(&mut node_exec, 800, 20.0); - assert_eq!(fft[0], (11, 999)); + assert_eq!(fft[0], (11, 1000)); matrix.set_param(freq_p, SAtom::param(0.2)); let fft = run_and_get_fft4096(&mut node_exec, 800, 20.0); @@ -617,8 +617,8 @@ fn check_node_sampl_rev_2() { res.0, 5000, vec![ - 0.9999773, 0.886596, 0.77321476, 0.6598335, 0.5464522, 0.43307102, 0.31968975, - 0.20630851, 0.09292727 + 0.0, 0.88664144, 0.7732602, 0.6598789, 0.54649764, 0.4331164, 0.31973514, 0.20635389, + 0.09297263 ] );