From b66fa8163c45a6684a70d3e0caad5cd2b31461d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Mon, 17 Jun 2024 14:07:42 +0200 Subject: [PATCH 01/16] cic/sweptsine wip --- src/cic.rs | 121 +++++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 4 ++ src/sweptsine.rs | 29 ++++++++++++ 3 files changed, 154 insertions(+) create mode 100644 src/cic.rs create mode 100644 src/sweptsine.rs diff --git a/src/cic.rs b/src/cic.rs new file mode 100644 index 0000000..717bf6b --- /dev/null +++ b/src/cic.rs @@ -0,0 +1,121 @@ +use core::iter::{repeat, Repeat, Take}; + +use num_traits::{Num, WrappingAdd}; + +/// Combs stage +#[derive(Clone, PartialEq, PartialOrd, Debug)] +pub struct Comb([T; N]); + +impl Default for Comb { + fn default() -> Self { + Self([T::zero(); N]) + } +} + +impl Comb { + /// Ingest a new sample into the filter and return its current output. + pub fn update(&mut self, x: T) -> T { + self.0.iter_mut().fold(x, |x, c| { + let y = x - *c; + *c = x; + y + }) + } +} + +/// Integrators stage +#[derive(Clone, PartialEq, PartialOrd, Debug)] +pub struct Integrator([T; N]); + +impl Default for Integrator { + fn default() -> Self { + Self([T::zero(); N]) + } +} + +impl Integrator { + /// Ingest a new sample into the filter and return its current output. + pub fn update(&mut self, x: T) -> T { + self.0.iter_mut().fold(x, |x, i| { + *i = x.wrapping_add(i); + *i + }) + } +} + +/// Cascaded integrator comb interpolatpr +/// +/// Order `N` where `N = 3` is cubic. +#[derive(Clone, Debug)] +pub struct CicInterpolator { + rate: usize, + combs: Comb, + up: Take>, + integrators: Integrator, +} + +impl CicInterpolator { + /// Return the filter gain + pub fn gain(&self) -> usize { + self.rate.pow(N as _) + } + + /// Create a new zero-initialized filter with the given rate change. + pub fn new(rate: usize) -> Self { + Self { + rate, + combs: Default::default(), + up: repeat(T::zero()).take(rate), + integrators: Default::default(), + } + } + + /// Optionally ingest a new low-rate sample and + /// retrieve the next output. + pub fn update(&mut self, x: Option) -> T { + if let Some(x) = x { + self.up = repeat(self.combs.update(x)).take(self.rate); + } + self.integrators.update(self.up.next().unwrap()) + } +} + +/// Cascaded integrator comb decimator. +/// +/// Order `N` where `N = 3` is cubic. +#[derive(Clone, Debug)] +pub struct CicDecimator { + rate: usize, + integrators: Integrator, + index: usize, + combs: Comb, +} + +impl CicDecimator { + /// Return the filter gain + pub fn gain(&self) -> usize { + self.rate.pow(N as _) + } + + /// Create a new zero-initialized filter with the given rate change. + pub fn new(rate: usize) -> Self { + Self { + rate, + combs: Default::default(), + index: 0, + integrators: Default::default(), + } + } + + /// Ingest a new high-rate sample and optionally retrieve next output. + pub fn update(&mut self, x: T) -> Option { + let x = self.integrators.update(x); + if self.index == 0 { + self.index = self.rate; + Some(self.combs.update(x)) + } else { + self.index -= 1; + None + } + } +} diff --git a/src/lib.rs b/src/lib.rs index 6fa71cd..ce97191 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -32,6 +32,10 @@ pub use num::*; mod dsm; pub mod svf; pub use dsm::*; +mod sweptsine; +pub use sweptsine::*; +mod cic; +pub use cic::*; #[cfg(test)] pub mod testing; diff --git a/src/sweptsine.rs b/src/sweptsine.rs new file mode 100644 index 0000000..96a36f9 --- /dev/null +++ b/src/sweptsine.rs @@ -0,0 +1,29 @@ +/// Exponential sweep generator +pub struct SweepIter { + /// Rate of increase + pub a: u32, + /// Current + pub f: u64, +} + +impl Iterator for SweepIter { + type Item = u64; + + fn next(&mut self) -> Option { + const BIAS: u64 = 1 << (32 - 1); + self.f += (self.a as u64) * ((self.f + BIAS) >> 32); + Some(self.f) + } +} + +pub fn sweep(k: u32, f1: u32, f2: u32) -> impl Iterator { + let it = SweepIter { + a: 1, + f: (f1 as u64) << 32, + }; + it.take_while(move |f| (f >> 32) as u32 <= f2) + .scan(0u64, |p, f| { + *p = p.wrapping_add(f); + Some((*p >> 32) as _) + }) +} From 8230e85692d8e4e691cd28621018296fd3eb66f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Mon, 17 Jun 2024 22:47:59 +0200 Subject: [PATCH 02/16] simplify cic --- src/cic.rs | 81 +++++++++++++++++++++++------------------------------- 1 file changed, 35 insertions(+), 46 deletions(-) diff --git a/src/cic.rs b/src/cic.rs index 717bf6b..2f73b51 100644 --- a/src/cic.rs +++ b/src/cic.rs @@ -2,44 +2,28 @@ use core::iter::{repeat, Repeat, Take}; use num_traits::{Num, WrappingAdd}; -/// Combs stage -#[derive(Clone, PartialEq, PartialOrd, Debug)] -pub struct Comb([T; N]); +/// Comb stage +#[derive(Clone, Copy, PartialEq, PartialOrd, Debug, Default)] +pub struct Comb(T); -impl Default for Comb { - fn default() -> Self { - Self([T::zero(); N]) - } -} - -impl Comb { +impl Comb { /// Ingest a new sample into the filter and return its current output. pub fn update(&mut self, x: T) -> T { - self.0.iter_mut().fold(x, |x, c| { - let y = x - *c; - *c = x; - y - }) + let y = x - self.0; + self.0 = x; + y } } /// Integrators stage -#[derive(Clone, PartialEq, PartialOrd, Debug)] -pub struct Integrator([T; N]); +#[derive(Clone, Copy, PartialEq, PartialOrd, Debug, Default)] +pub struct Integrator(T); -impl Default for Integrator { - fn default() -> Self { - Self([T::zero(); N]) - } -} - -impl Integrator { +impl Integrator { /// Ingest a new sample into the filter and return its current output. pub fn update(&mut self, x: T) -> T { - self.0.iter_mut().fold(x, |x, i| { - *i = x.wrapping_add(i); - *i - }) + self.0 = self.0.wrapping_add(&x); + self.0 } } @@ -47,16 +31,16 @@ impl Integrator { /// /// Order `N` where `N = 3` is cubic. #[derive(Clone, Debug)] -pub struct CicInterpolator { +pub struct CicInterpolator { rate: usize, - combs: Comb, + combs: [Comb; N], up: Take>, - integrators: Integrator, + integrators: [Integrator; N], } -impl CicInterpolator { +impl CicInterpolator { /// Return the filter gain - pub fn gain(&self) -> usize { + pub const fn gain(&self) -> usize { self.rate.pow(N as _) } @@ -64,9 +48,9 @@ impl CicInterpolator { pub fn new(rate: usize) -> Self { Self { rate, - combs: Default::default(), + combs: [Comb::default(); N], up: repeat(T::zero()).take(rate), - integrators: Default::default(), + integrators: [Integrator::default(); N], } } @@ -74,9 +58,12 @@ impl CicInterpolator { /// retrieve the next output. pub fn update(&mut self, x: Option) -> T { if let Some(x) = x { - self.up = repeat(self.combs.update(x)).take(self.rate); + let x = self.combs.iter_mut().fold(x, |x, c| c.update(x)); + self.up = repeat(x).take(self.rate); } - self.integrators.update(self.up.next().unwrap()) + self.integrators + .iter_mut() + .fold(self.up.next().unwrap(), |x, i| i.update(x)) } } @@ -84,14 +71,14 @@ impl CicInterpolator { /// /// Order `N` where `N = 3` is cubic. #[derive(Clone, Debug)] -pub struct CicDecimator { +pub struct CicDecimator { rate: usize, - integrators: Integrator, + integrators: [Integrator; N], index: usize, - combs: Comb, + combs: [Comb; N], } -impl CicDecimator { +impl CicDecimator { /// Return the filter gain pub fn gain(&self) -> usize { self.rate.pow(N as _) @@ -99,20 +86,22 @@ impl CicDecimator { /// Create a new zero-initialized filter with the given rate change. pub fn new(rate: usize) -> Self { + debug_assert!(rate > 0); Self { rate, - combs: Default::default(), + combs: [Comb::default(); N], index: 0, - integrators: Default::default(), + integrators: [Integrator::default(); N], } } /// Ingest a new high-rate sample and optionally retrieve next output. pub fn update(&mut self, x: T) -> Option { - let x = self.integrators.update(x); + let x = self.integrators.iter_mut().fold(x, |x, i| i.update(x)); if self.index == 0 { - self.index = self.rate; - Some(self.combs.update(x)) + self.index = self.rate - 1; + let x = self.combs.iter_mut().fold(x, |x, c| c.update(x)); + Some(x) } else { self.index -= 1; None From 0529d5307a1f55caf3f36d46dbb8efa25ac7f465 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Tue, 18 Jun 2024 10:21:05 +0200 Subject: [PATCH 03/16] cic tweaks --- src/cic.rs | 105 +++++++++++++++++++++++++---------------------------- 1 file changed, 50 insertions(+), 55 deletions(-) diff --git a/src/cic.rs b/src/cic.rs index 2f73b51..adbccec 100644 --- a/src/cic.rs +++ b/src/cic.rs @@ -1,56 +1,29 @@ -use core::iter::{repeat, Repeat, Take}; +use core::ops::AddAssign; use num_traits::{Num, WrappingAdd}; -/// Comb stage -#[derive(Clone, Copy, PartialEq, PartialOrd, Debug, Default)] -pub struct Comb(T); - -impl Comb { - /// Ingest a new sample into the filter and return its current output. - pub fn update(&mut self, x: T) -> T { - let y = x - self.0; - self.0 = x; - y - } -} - -/// Integrators stage -#[derive(Clone, Copy, PartialEq, PartialOrd, Debug, Default)] -pub struct Integrator(T); - -impl Integrator { - /// Ingest a new sample into the filter and return its current output. - pub fn update(&mut self, x: T) -> T { - self.0 = self.0.wrapping_add(&x); - self.0 - } -} - /// Cascaded integrator comb interpolatpr /// /// Order `N` where `N = 3` is cubic. #[derive(Clone, Debug)] pub struct CicInterpolator { rate: usize, - combs: [Comb; N], - up: Take>, - integrators: [Integrator; N], + combs: [T; N], + up: T, + index: usize, + integrators: [T; N], } -impl CicInterpolator { - /// Return the filter gain - pub const fn gain(&self) -> usize { - self.rate.pow(N as _) - } - +impl CicInterpolator { /// Create a new zero-initialized filter with the given rate change. pub fn new(rate: usize) -> Self { + debug_assert!(rate > 0); Self { rate, - combs: [Comb::default(); N], - up: repeat(T::zero()).take(rate), - integrators: [Integrator::default(); N], + combs: [T::zero(); N], + up: T::zero(), + index: 0, + integrators: [T::zero(); N], } } @@ -58,12 +31,26 @@ impl CicInterpolator) -> T { if let Some(x) = x { - let x = self.combs.iter_mut().fold(x, |x, c| c.update(x)); - self.up = repeat(x).take(self.rate); + let x = self.combs.iter_mut().fold(x, |x, c| { + let y = x - *c; + *c = x; + y + }); + debug_assert_eq!(self.index, 0); + self.index = self.rate - 1; + self.up = x; + } else { + self.index -= 1; } - self.integrators - .iter_mut() - .fold(self.up.next().unwrap(), |x, i| i.update(x)) + self.integrators.iter_mut().fold(self.up, |x, i| { + *i += x; + *i + }) + } + + /// Return the filter gain + pub const fn gain(&self) -> usize { + self.rate.pow(N as _) } } @@ -73,38 +60,46 @@ impl CicInterpolator { rate: usize, - integrators: [Integrator; N], + integrators: [T; N], index: usize, - combs: [Comb; N], + combs: [T; N], } -impl CicDecimator { - /// Return the filter gain - pub fn gain(&self) -> usize { - self.rate.pow(N as _) - } - +impl CicDecimator { /// Create a new zero-initialized filter with the given rate change. pub fn new(rate: usize) -> Self { debug_assert!(rate > 0); Self { rate, - combs: [Comb::default(); N], + combs: [T::zero(); N], index: 0, - integrators: [Integrator::default(); N], + integrators: [T::zero(); N], } } /// Ingest a new high-rate sample and optionally retrieve next output. pub fn update(&mut self, x: T) -> Option { - let x = self.integrators.iter_mut().fold(x, |x, i| i.update(x)); + let x = self.integrators.iter_mut().fold(x, |x, i| { + // Overflow is OK if bitwidth is sufficient (input * gain) + *i = i.wrapping_add(&x); + *i + }); if self.index == 0 { self.index = self.rate - 1; - let x = self.combs.iter_mut().fold(x, |x, c| c.update(x)); + let x = self.combs.iter_mut().fold(x, |x, c| { + let y = x - *c; + *c = x; + y + }); Some(x) } else { self.index -= 1; None } } + + /// Return the filter gain + pub fn gain(&self) -> usize { + self.rate.pow(N as _) + } } From 5c9afb5f16623a31a5995fcb06e4978afc0e6e42 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Tue, 18 Jun 2024 10:29:03 +0200 Subject: [PATCH 04/16] cic docs --- src/cic.rs | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/cic.rs b/src/cic.rs index adbccec..2cdbe80 100644 --- a/src/cic.rs +++ b/src/cic.rs @@ -7,10 +7,16 @@ use num_traits::{Num, WrappingAdd}; /// Order `N` where `N = 3` is cubic. #[derive(Clone, Debug)] pub struct CicInterpolator { + /// Rate change (> 0) rate: usize, + /// Comb/differentiator stages combs: [T; N], - up: T, + /// Zero order hold in the middle combined with the upsampler + /// This is equivalent to a single comb + (delta-)upsampler + integrator + zoh: T, + /// Rate change index (count down) index: usize, + /// Integrator stages integrators: [T; N], } @@ -21,7 +27,7 @@ impl CicInterpolator { Self { rate, combs: [T::zero(); N], - up: T::zero(), + zoh: T::zero(), index: 0, integrators: [T::zero(); N], } @@ -38,11 +44,11 @@ impl CicInterpolator { }); debug_assert_eq!(self.index, 0); self.index = self.rate - 1; - self.up = x; + self.zoh = x; } else { self.index -= 1; } - self.integrators.iter_mut().fold(self.up, |x, i| { + self.integrators.iter_mut().fold(self.zoh, |x, i| { *i += x; *i }) @@ -59,9 +65,13 @@ impl CicInterpolator { /// Order `N` where `N = 3` is cubic. #[derive(Clone, Debug)] pub struct CicDecimator { + /// Rate change (> 0) rate: usize, + /// Integration stages integrators: [T; N], + /// Rate change count down index: usize, + /// Comb/differentiator stages combs: [T; N], } From f171d04efbbcbdbc9cde8e27116326a466897978 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Tue, 18 Jun 2024 16:07:54 +0200 Subject: [PATCH 05/16] cic docs --- src/cic.rs | 173 +++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 128 insertions(+), 45 deletions(-) diff --git a/src/cic.rs b/src/cic.rs index 2cdbe80..5790ec7 100644 --- a/src/cic.rs +++ b/src/cic.rs @@ -1,17 +1,19 @@ use core::ops::AddAssign; -use num_traits::{Num, WrappingAdd}; +use num_traits::{AsPrimitive, Num, WrappingAdd}; -/// Cascaded integrator comb interpolatpr +/// Cascaded integrator comb structure /// /// Order `N` where `N = 3` is cubic. #[derive(Clone, Debug)] -pub struct CicInterpolator { +pub struct Cic { /// Rate change (> 0) rate: usize, /// Comb/differentiator stages combs: [T; N], - /// Zero order hold in the middle combined with the upsampler + /// Zero order hold behind comb sections. + /// Interpolator: In the middle combined with the upsampler + /// Decimator: After combs to support `get()` /// This is equivalent to a single comb + (delta-)upsampler + integrator zoh: T, /// Rate change index (count down) @@ -20,7 +22,11 @@ pub struct CicInterpolator { integrators: [T; N], } -impl CicInterpolator { +impl Cic +where + T: Num + AddAssign + WrappingAdd + Copy + 'static, + usize: AsPrimitive, +{ /// Create a new zero-initialized filter with the given rate change. pub fn new(rate: usize) -> Self { debug_assert!(rate > 0); @@ -33,17 +39,44 @@ impl CicInterpolator { } } + pub fn rate(&self) -> usize { + self.rate + } + + pub fn set_rate(&mut self, rate: usize) { + debug_assert!(rate > 0); + self.rate = rate; + } + + pub fn clear(&mut self) { + self.zoh = T::zero(); + self.combs = [T::zero(); N]; + self.integrators = [T::zero(); N]; + } + + pub fn settle_interpolate(&mut self, x: T) { + self.clear(); + let g = self.gain(); + let i = self.integrators.last_mut().unwrap_or(&mut self.zoh); + *i = x * g; + } + + pub fn settle_decimate(&mut self, x: T) { + self.clear(); + self.zoh = x * self.gain(); + } + /// Optionally ingest a new low-rate sample and /// retrieve the next output. - pub fn update(&mut self, x: Option) -> T { + pub fn interpolate(&mut self, x: Option) -> T { if let Some(x) = x { + debug_assert_eq!(self.index, 0); + self.index = self.rate - 1; let x = self.combs.iter_mut().fold(x, |x, c| { let y = x - *c; *c = x; y }); - debug_assert_eq!(self.index, 0); - self.index = self.rate - 1; self.zoh = x; } else { self.index -= 1; @@ -54,41 +87,8 @@ impl CicInterpolator { }) } - /// Return the filter gain - pub const fn gain(&self) -> usize { - self.rate.pow(N as _) - } -} - -/// Cascaded integrator comb decimator. -/// -/// Order `N` where `N = 3` is cubic. -#[derive(Clone, Debug)] -pub struct CicDecimator { - /// Rate change (> 0) - rate: usize, - /// Integration stages - integrators: [T; N], - /// Rate change count down - index: usize, - /// Comb/differentiator stages - combs: [T; N], -} - -impl CicDecimator { - /// Create a new zero-initialized filter with the given rate change. - pub fn new(rate: usize) -> Self { - debug_assert!(rate > 0); - Self { - rate, - combs: [T::zero(); N], - index: 0, - integrators: [T::zero(); N], - } - } - /// Ingest a new high-rate sample and optionally retrieve next output. - pub fn update(&mut self, x: T) -> Option { + pub fn decimate(&mut self, x: T) -> Option { let x = self.integrators.iter_mut().fold(x, |x, i| { // Overflow is OK if bitwidth is sufficient (input * gain) *i = i.wrapping_add(&x); @@ -101,15 +101,98 @@ impl CicDecimator { *c = x; y }); - Some(x) + self.zoh = x; + Some(self.zoh) } else { self.index -= 1; None } } + pub fn tick(&self) -> bool { + self.index == 0 + } + + pub fn get_interpolate(&self) -> T { + *self.integrators.last().unwrap_or(&self.zoh) + } + + pub fn get_decimate(&self) -> T { + self.zoh + } + /// Return the filter gain - pub fn gain(&self) -> usize { - self.rate.pow(N as _) + pub fn gain(&self) -> T { + self.rate.pow(N as _).as_() + } + + /// Return the right shift amount (the log2 of gain()) + /// + /// Panics if the gain is not a power of two. + pub fn log2_gain(&self) -> u32 { + let s = (31 - self.rate.leading_zeros()) * N as u32; + debug_assert_eq!(1 << s, self.rate.pow(N as _) as _); + s + } + + /// Return the impulse response length + pub fn response_length(&self) -> usize { + N * (self.rate - 1) + } +} + +#[cfg(test)] +mod test { + use super::*; + + #[test] + fn new() { + let _ = Cic::::new(1); + } + + #[test] + fn identity() { + let mut int = Cic::::new(1); + let mut dec = Cic::::new(1); + for x in 0..100 { + assert_eq!(int.interpolate(Some(x)), x); + assert_eq!(dec.decimate(x), Some(x)); + assert_eq!(x, int.get_interpolate()); + assert_eq!(x, dec.get_decimate()); + } + } + + #[test] + fn response_length_gain() { + let mut int = Cic::::new(33); + let x = 99; + for i in 0..2 * int.response_length() { + let y = int.interpolate(if int.tick() { Some(x) } else { None }); + assert_eq!(y, int.get_interpolate()); + if i < int.response_length() { + assert!(y < x * int.gain()); + } else { + assert_eq!(y, x * int.gain()); + } + } + } + + #[test] + fn settle() { + let x = 99; + let mut int = Cic::::new(33); + int.settle_interpolate(x); + let mut dec = Cic::::new(33); + dec.settle_decimate(x); + for _ in 0..100 { + let y = int.interpolate(if int.tick() { Some(x) } else { None }); + assert_eq!(y, x*int.gain()); + assert_eq!(y, int.get_interpolate()); + assert_eq!(dec.get_decimate(), x*dec.gain()); + if let Some(y) = dec.decimate(x) { + assert_eq!(y, dec.get_decimate()); + } + } + } } From 192a5e8a5324f621b68a102d203bfaecfd107938 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Tue, 18 Jun 2024 21:02:48 +0200 Subject: [PATCH 06/16] cic tests --- Cargo.toml | 2 + src/cic.rs | 207 ++++++++++++++++++++++++++++++----------------- src/sweptsine.rs | 2 +- 3 files changed, 138 insertions(+), 73 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 65678b6..0f35cbe 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -16,6 +16,8 @@ num-complex = { version = "0.4.0", features = ["serde"], default-features = fals num-traits = { version = "0.2.14", features = ["libm"], default-features = false} [dev-dependencies] +quickcheck = "1.0.3" +quickcheck_macros = "1.0.0" rand = "0.8" rustfft = "6.1.0" # futuredsp = "0.0.6" diff --git a/src/cic.rs b/src/cic.rs index 5790ec7..244a659 100644 --- a/src/cic.rs +++ b/src/cic.rs @@ -1,77 +1,109 @@ use core::ops::AddAssign; -use num_traits::{AsPrimitive, Num, WrappingAdd}; +use num_traits::{AsPrimitive, Num, Pow, WrappingAdd, WrappingSub}; /// Cascaded integrator comb structure /// /// Order `N` where `N = 3` is cubic. #[derive(Clone, Debug)] pub struct Cic { - /// Rate change (> 0) - rate: usize, - /// Comb/differentiator stages - combs: [T; N], + /// Rate change (fast/slow - 1) + /// Interpolator: output/input - 1 + /// Decimator: input/output - 1 + rate: u32, + /// Up/downsampler state (count down) + index: u32, /// Zero order hold behind comb sections. /// Interpolator: In the middle combined with the upsampler - /// Decimator: After combs to support `get()` + /// Decimator: After combs to support `get_decimate()` /// This is equivalent to a single comb + (delta-)upsampler + integrator zoh: T, - /// Rate change index (count down) - index: usize, - /// Integrator stages + /// Comb/differentiator state + combs: [T; N], + /// Integrator state integrators: [T; N], } impl Cic where - T: Num + AddAssign + WrappingAdd + Copy + 'static, - usize: AsPrimitive, + T: Num + AddAssign + WrappingAdd + WrappingSub + Pow + Copy + 'static, + u32: AsPrimitive, { /// Create a new zero-initialized filter with the given rate change. - pub fn new(rate: usize) -> Self { - debug_assert!(rate > 0); + pub fn new(rate: u32) -> Self { Self { rate, - combs: [T::zero(); N], - zoh: T::zero(), index: 0, + zoh: T::zero(), + combs: [T::zero(); N], integrators: [T::zero(); N], } } - pub fn rate(&self) -> usize { + /// Filter order + /// + /// * 0: zero order hold + /// * 1: linear + /// * 2: quadratic + /// * 3: cubic interpolation/decimation + /// + /// etc. + pub const fn order(&self) -> usize { + N + } + + /// Rate change + /// + /// `fast/slow - 1` + pub const fn rate(&self) -> u32 { self.rate } - pub fn set_rate(&mut self, rate: usize) { - debug_assert!(rate > 0); + /// Set the rate change + /// + /// `fast/slow - 1` + pub fn set_rate(&mut self, rate: u32) { self.rate = rate; } + /// Zero-initialize the filter state pub fn clear(&mut self) { + self.index = 0; self.zoh = T::zero(); self.combs = [T::zero(); N]; self.integrators = [T::zero(); N]; } + /// Establish a settled filter state pub fn settle_interpolate(&mut self, x: T) { self.clear(); + *self.combs.first_mut().unwrap_or(&mut self.zoh) = x; let g = self.gain(); - let i = self.integrators.last_mut().unwrap_or(&mut self.zoh); - *i = x * g; + if let Some(i) = self.integrators.last_mut() { + *i = x * g; + } } + /// Establish a settled filter state pub fn settle_decimate(&mut self, x: T) { - self.clear(); + for i in self.integrators.iter_mut() { + *i = x; + // x *= self.rate; + } + for c in self.combs.iter_mut() { + *c = x; + } self.zoh = x * self.gain(); } /// Optionally ingest a new low-rate sample and /// retrieve the next output. + /// + /// A new sample must be supplied at the correct time (see [`Cic::tick()`]) pub fn interpolate(&mut self, x: Option) -> T { if let Some(x) = x { debug_assert_eq!(self.index, 0); - self.index = self.rate - 1; + self.index = self.rate; let x = self.combs.iter_mut().fold(x, |x, c| { let y = x - *c; *c = x; @@ -95,9 +127,9 @@ where *i }); if self.index == 0 { - self.index = self.rate - 1; + self.index = self.rate; let x = self.combs.iter_mut().fold(x, |x, c| { - let y = x - *c; + let y = x.wrapping_sub(c); *c = x; y }); @@ -109,90 +141,121 @@ where } } - pub fn tick(&self) -> bool { + /// Accepts/provides new slow-rate sample + /// + /// Interpolator: accepts new input sample + /// Decimator: returns new output sample + pub const fn tick(&self) -> bool { self.index == 0 } + /// Current interpolator output pub fn get_interpolate(&self) -> T { *self.integrators.last().unwrap_or(&self.zoh) } + /// Current decimator output pub fn get_decimate(&self) -> T { self.zoh } - /// Return the filter gain + /// Filter gain pub fn gain(&self) -> T { - self.rate.pow(N as _).as_() + (self.rate.as_() + T::one()).pow(N) } - /// Return the right shift amount (the log2 of gain()) + /// Right shift amount /// - /// Panics if the gain is not a power of two. - pub fn log2_gain(&self) -> u32 { - let s = (31 - self.rate.leading_zeros()) * N as u32; - debug_assert_eq!(1 << s, self.rate.pow(N as _) as _); - s + /// `log2(gain())` if gain is a power of two, + /// otherwise an upper bound. + pub const fn gain_log2(&self) -> u32 { + (u32::BITS - self.rate.leading_zeros()) * N as u32 } - /// Return the impulse response length - pub fn response_length(&self) -> usize { - N * (self.rate - 1) + /// Impulse response length + pub const fn response_length(&self) -> usize { + self.rate as usize * N } } #[cfg(test)] mod test { + use core::cmp::Ordering; + use super::*; + use quickcheck_macros::quickcheck; - #[test] - fn new() { - let _ = Cic::::new(1); + #[quickcheck] + fn new(rate: u32) { + let _ = Cic::::new(rate); } - #[test] - fn identity() { - let mut int = Cic::::new(1); - let mut dec = Cic::::new(1); - for x in 0..100 { - assert_eq!(int.interpolate(Some(x)), x); - assert_eq!(dec.decimate(x), Some(x)); - assert_eq!(x, int.get_interpolate()); + #[quickcheck] + fn identity_dec(x: Vec) { + let mut dec = Cic::<_, 3>::new(0); + for x in x { + assert_eq!(x, dec.decimate(x).unwrap()); assert_eq!(x, dec.get_decimate()); } } - #[test] - fn response_length_gain() { - let mut int = Cic::::new(33); - let x = 99; - for i in 0..2 * int.response_length() { - let y = int.interpolate(if int.tick() { Some(x) } else { None }); - assert_eq!(y, int.get_interpolate()); - if i < int.response_length() { - assert!(y < x * int.gain()); - } else { - assert_eq!(y, x * int.gain()); + #[quickcheck] + fn identity_int(x: Vec) { + const N: usize = 3; + let mut int = Cic::<_, N>::new(0); + for x in x { + assert_eq!(x >> N, int.interpolate(Some(x >> N))); + assert_eq!(x >> N, int.get_interpolate()); + } + } + + #[quickcheck] + fn response_length_gain_settle(x: Vec, rate: u32) { + let mut int = Cic::<_, 3>::new(rate); + let shift = int.gain_log2(); + if shift >= 32 { + return; + } + assert!(int.gain() <= 1 << shift); + for x in x { + while !int.tick() { + int.interpolate(None); + } + let y_last = int.get_interpolate(); + let y_want = x as i64 * int.gain(); + for i in 0..2 * int.response_length() { + let y = int.interpolate(if int.tick() { Some(x as i64) } else { None }); + assert_eq!(y, int.get_interpolate()); + if i < int.response_length() { + match y_want.cmp(&y_last) { + Ordering::Greater => assert!((y_last..y_want).contains(&y)), + Ordering::Less => assert!((y_want..y_last).contains(&(y - 1))), + Ordering::Equal => assert_eq!(y_want, y), + } + } else { + assert_eq!(y, y_want); + } } } } - #[test] - fn settle() { - let x = 99; - let mut int = Cic::::new(33); - int.settle_interpolate(x); - let mut dec = Cic::::new(33); - dec.settle_decimate(x); + #[quickcheck] + fn settle(rate: u32, x: i32) { + let mut int = Cic::::new(rate); + if int.gain_log2() >= 32 { + return; + } + int.settle_interpolate(x as _); + //let mut dec = Cic::::new(rate); + //dec.settle_decimate(x as _); for _ in 0..100 { - let y = int.interpolate(if int.tick() { Some(x) } else { None }); - assert_eq!(y, x*int.gain()); + let y = int.interpolate(if int.tick() { Some(x as _) } else { None }); + assert_eq!(y, x as i64 * int.gain()); assert_eq!(y, int.get_interpolate()); - assert_eq!(dec.get_decimate(), x*dec.gain()); - if let Some(y) = dec.decimate(x) { - assert_eq!(y, dec.get_decimate()); - } + //assert_eq!(dec.get_decimate(), x*dec.gain()); + //if let Some(y) = dec.decimate(x) { + // assert_eq!(y, dec.get_decimate()); + //} } - } } diff --git a/src/sweptsine.rs b/src/sweptsine.rs index 96a36f9..cce1d4a 100644 --- a/src/sweptsine.rs +++ b/src/sweptsine.rs @@ -16,7 +16,7 @@ impl Iterator for SweepIter { } } -pub fn sweep(k: u32, f1: u32, f2: u32) -> impl Iterator { +pub fn sweep(_k: u32, f1: u32, f2: u32) -> impl Iterator { let it = SweepIter { a: 1, f: (f1 as u64) << 32, From ca634a8469fb4b3710f9d3a3526e237330739456 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Wed, 19 Jun 2024 00:36:15 +0200 Subject: [PATCH 07/16] cic tests --- src/cic.rs | 95 ++++++++++++++++++++++++++---------------------------- 1 file changed, 46 insertions(+), 49 deletions(-) diff --git a/src/cic.rs b/src/cic.rs index 244a659..c3ccd38 100644 --- a/src/cic.rs +++ b/src/cic.rs @@ -74,6 +74,42 @@ where self.integrators = [T::zero(); N]; } + /// Accepts/provides new slow-rate sample + /// + /// Interpolator: accepts new input sample + /// Decimator: returns new output sample + pub const fn tick(&self) -> bool { + self.index == 0 + } + + /// Current interpolator output + pub fn get_interpolate(&self) -> T { + *self.integrators.last().unwrap_or(&self.zoh) + } + + /// Current decimator output + pub fn get_decimate(&self) -> T { + self.zoh + } + + /// Filter gain + pub fn gain(&self) -> T { + (self.rate.as_() + T::one()).pow(N) + } + + /// Right shift amount + /// + /// `log2(gain())` if gain is a power of two, + /// otherwise an upper bound. + pub const fn gain_log2(&self) -> u32 { + (u32::BITS - self.rate.leading_zeros()) * N as u32 + } + + /// Impulse response length + pub const fn response_length(&self) -> usize { + self.rate as usize * N + } + /// Establish a settled filter state pub fn settle_interpolate(&mut self, x: T) { self.clear(); @@ -85,15 +121,12 @@ where } /// Establish a settled filter state + /// + /// Unimplemented! pub fn settle_decimate(&mut self, x: T) { - for i in self.integrators.iter_mut() { - *i = x; - // x *= self.rate; - } - for c in self.combs.iter_mut() { - *c = x; - } + self.clear(); self.zoh = x * self.gain(); + unimplemented!(); } /// Optionally ingest a new low-rate sample and @@ -140,42 +173,6 @@ where None } } - - /// Accepts/provides new slow-rate sample - /// - /// Interpolator: accepts new input sample - /// Decimator: returns new output sample - pub const fn tick(&self) -> bool { - self.index == 0 - } - - /// Current interpolator output - pub fn get_interpolate(&self) -> T { - *self.integrators.last().unwrap_or(&self.zoh) - } - - /// Current decimator output - pub fn get_decimate(&self) -> T { - self.zoh - } - - /// Filter gain - pub fn gain(&self) -> T { - (self.rate.as_() + T::one()).pow(N) - } - - /// Right shift amount - /// - /// `log2(gain())` if gain is a power of two, - /// otherwise an upper bound. - pub const fn gain_log2(&self) -> u32 { - (u32::BITS - self.rate.leading_zeros()) * N as u32 - } - - /// Impulse response length - pub const fn response_length(&self) -> usize { - self.rate as usize * N - } } #[cfg(test)] @@ -246,16 +243,16 @@ mod test { return; } int.settle_interpolate(x as _); - //let mut dec = Cic::::new(rate); - //dec.settle_decimate(x as _); + // let mut dec = Cic::::new(rate); + // dec.settle_decimate(x as _); for _ in 0..100 { let y = int.interpolate(if int.tick() { Some(x as _) } else { None }); assert_eq!(y, x as i64 * int.gain()); assert_eq!(y, int.get_interpolate()); - //assert_eq!(dec.get_decimate(), x*dec.gain()); - //if let Some(y) = dec.decimate(x) { - // assert_eq!(y, dec.get_decimate()); - //} + // assert_eq!(dec.get_decimate(), x as i64 * dec.gain()); + // if let Some(y) = dec.decimate(x as _) { + // assert_eq!(y, x as i64 * dec.gain()); + // } } } } From c865db4e909bf9c82c1069a2352c675c25351726 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Thu, 20 Jun 2024 21:40:58 +0200 Subject: [PATCH 08/16] sweep wip --- Cargo.toml | 3 +++ src/cic.rs | 1 - src/sweptsine.rs | 52 ++++++++++++++++++++++++++++++++++++------------ 3 files changed, 42 insertions(+), 14 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 0f35cbe..d841170 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -23,5 +23,8 @@ rustfft = "6.1.0" # futuredsp = "0.0.6" # sdr = "0.7.0" +[features] +std = [] + [profile.release] debug = 1 diff --git a/src/cic.rs b/src/cic.rs index c3ccd38..94bbbd5 100644 --- a/src/cic.rs +++ b/src/cic.rs @@ -16,7 +16,6 @@ pub struct Cic { /// Zero order hold behind comb sections. /// Interpolator: In the middle combined with the upsampler /// Decimator: After combs to support `get_decimate()` - /// This is equivalent to a single comb + (delta-)upsampler + integrator zoh: T, /// Comb/differentiator state combs: [T; N], diff --git a/src/sweptsine.rs b/src/sweptsine.rs index cce1d4a..4341bc0 100644 --- a/src/sweptsine.rs +++ b/src/sweptsine.rs @@ -1,29 +1,55 @@ +use crate::cossin; + /// Exponential sweep generator pub struct SweepIter { /// Rate of increase - pub a: u32, + pub a: i32, /// Current - pub f: u64, + pub f: i64, } impl Iterator for SweepIter { - type Item = u64; + type Item = i64; fn next(&mut self) -> Option { - const BIAS: u64 = 1 << (32 - 1); - self.f += (self.a as u64) * ((self.f + BIAS) >> 32); + self.f = self + .f + .wrapping_add((self.a as i64).wrapping_mul(self.f >> 32)); Some(self.f) } } -pub fn sweep(_k: u32, f1: u32, f2: u32) -> impl Iterator { - let it = SweepIter { - a: 1, - f: (f1 as u64) << 32, - }; - it.take_while(move |f| (f >> 32) as u32 <= f2) - .scan(0u64, |p, f| { +impl SweepIter { + /// Log swept frequency + pub fn sweep(a: i32, f1: i32, f2: i32) -> impl Iterator { + Self { + a, + f: (f1 as i64) << 32, + } + .map(|f| (f >> 32) as i32) + .take_while(move |f| *f <= f2) + } + + /// Log swept sine + pub fn swept_sine(a: i32, p0: i32, f1: i32, f2: i32) -> impl Iterator { + Self::sweep(a, f1, f2).scan(p0, |p, f| { *p = p.wrapping_add(f); - Some((*p >> 32) as _) + Some(cossin::cossin(*p).1) }) + } +} + +#[cfg(test)] +mod test { + use super::*; + // use quickcheck_macros::quickcheck; + + #[test] + fn new() { + let r = 0x3654; + let f1 = r; + let f2: i32 = i32::MAX >> 1; + let it = SweepIter::swept_sine(r, 0x0, f1, f2); + println!("{}", it.count()); + } } From 2e7445d4167d40bc66710134ff01b2ec5684f42b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Mon, 28 Oct 2024 22:32:06 +0100 Subject: [PATCH 09/16] sweptsine: restructure, good math --- src/sweptsine.rs | 119 +++++++++++++++++++++++++++++++++++++---------- 1 file changed, 95 insertions(+), 24 deletions(-) diff --git a/src/sweptsine.rs b/src/sweptsine.rs index 4341bc0..730ecaf 100644 --- a/src/sweptsine.rs +++ b/src/sweptsine.rs @@ -1,16 +1,19 @@ -use crate::cossin; +use crate::{cossin::cossin, Complex}; +use num_traits::Float; -/// Exponential sweep generator -pub struct SweepIter { +/// Exponential sweep +#[derive(Clone, Debug, PartialEq, PartialOrd)] +pub struct Sweep { /// Rate of increase pub a: i32, /// Current pub f: i64, } -impl Iterator for SweepIter { +impl Iterator for Sweep { type Item = i64; + #[inline] fn next(&mut self) -> Option { self.f = self .f @@ -19,22 +22,83 @@ impl Iterator for SweepIter { } } -impl SweepIter { - /// Log swept frequency - pub fn sweep(a: i32, f1: i32, f2: i32) -> impl Iterator { - Self { - a, - f: (f1 as i64) << 32, +impl Sweep { + /// Create a new exponential sweep + #[inline] + pub fn new(a: i32, f: i64) -> Self { + Self { a, f } + } +} + +/// Sync Sweep parameter error +#[derive(Clone, Copy, Debug, PartialEq, PartialOrd)] +pub enum Error { + /// Start out of bounds + Start, + /// End out of bounds + End, + /// Length invalid + Length, +} + +/// Exponential synchronized sweep +pub struct SyncExpSweep { + sweep: Sweep, + f_end: i32, +} + +impl SyncExpSweep { + /// Create new sweep + pub fn new(f_start: i32, octaves: u32, samples_octave: u32) -> Result { + if f_start.checked_ilog2().ok_or(Error::Start)? + octaves >= 31 { + return Err(Error::End); + } + if samples_octave == 0 { + return Err(Error::Length); } - .map(|f| (f >> 32) as i32) - .take_while(move |f| *f <= f2) + let a = Float::round( + (Float::exp2((samples_octave as f64).recip()) - 1.0) * (1i64 << 32) as f64, + ) as i32; + Ok(Self { + sweep: Sweep::new(a, (f_start as i64) << 32), + f_end: f_start << octaves, + }) } +} + +impl Iterator for SyncExpSweep { + type Item = i64; - /// Log swept sine - pub fn swept_sine(a: i32, p0: i32, f1: i32, f2: i32) -> impl Iterator { - Self::sweep(a, f1, f2).scan(p0, |p, f| { - *p = p.wrapping_add(f); - Some(cossin::cossin(*p).1) + #[inline] + fn next(&mut self) -> Option { + self.sweep + .next() + .filter(|f| ((f >> 32) as i32) < self.f_end) + } +} + +/// Accumulating oscillator +pub struct AccuOsc { + f: T, + p: i64, +} + +impl From for AccuOsc { + #[inline] + fn from(value: T) -> Self { + Self { f: value, p: 0 } + } +} + +impl> Iterator for AccuOsc { + type Item = Complex; + + #[inline] + fn next(&mut self) -> Option { + self.f.next().map(|f| { + self.p = self.p.wrapping_add(f); + let (re, im) = cossin((self.p >> 32) as _); + Complex { re, im } }) } } @@ -42,14 +106,21 @@ impl SweepIter { #[cfg(test)] mod test { use super::*; - // use quickcheck_macros::quickcheck; + use crate::ComplexExt as _; #[test] - fn new() { - let r = 0x3654; - let f1 = r; - let f2: i32 = i32::MAX >> 1; - let it = SweepIter::swept_sine(r, 0x0, f1, f2); - println!("{}", it.count()); + fn len() { + let f_start = 0x7fff; + let octaves = 16; + let samples_octave = 1 << 16; + let it: Vec<_> = + AccuOsc::from(SyncExpSweep::new(f_start, octaves, samples_octave).unwrap()).collect(); + let f1 = + it[it.len() - 1].arg().wrapping_sub(it[it.len() - 2].arg()) as f32 / f_start as f32; + println!("octaves {}", f1 / (1 << octaves) as f32 - 1.0); + println!( + "length {}", + it.len() as f32 / (octaves * samples_octave) as f32 - 1.0 + ); } } From d4b7c520245e5895b0b960ae65c79fb0351b4626 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Tue, 29 Oct 2024 17:03:22 +0000 Subject: [PATCH 10/16] sweptsine: traits, names, api --- src/accu.rs | 2 +- src/sweptsine.rs | 23 ++++++++++++++--------- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/src/accu.rs b/src/accu.rs index 6bf2d91..7f5502d 100644 --- a/src/accu.rs +++ b/src/accu.rs @@ -1,7 +1,7 @@ use num_traits::ops::wrapping::WrappingAdd; /// Wrapping Accumulator -#[derive(Copy, Clone, Default, PartialEq, Eq, Debug)] +#[derive(Copy, Clone, Default, PartialEq, PartialOrd, Debug)] pub struct Accu { state: T, step: T, diff --git a/src/sweptsine.rs b/src/sweptsine.rs index 730ecaf..4f79fe1 100644 --- a/src/sweptsine.rs +++ b/src/sweptsine.rs @@ -32,7 +32,7 @@ impl Sweep { /// Sync Sweep parameter error #[derive(Clone, Copy, Debug, PartialEq, PartialOrd)] -pub enum Error { +pub enum SweepError { /// Start out of bounds Start, /// End out of bounds @@ -42,6 +42,7 @@ pub enum Error { } /// Exponential synchronized sweep +#[derive(Clone, Debug, PartialEq, PartialOrd)] pub struct SyncExpSweep { sweep: Sweep, f_end: i32, @@ -49,16 +50,19 @@ pub struct SyncExpSweep { impl SyncExpSweep { /// Create new sweep - pub fn new(f_start: i32, octaves: u32, samples_octave: u32) -> Result { - if f_start.checked_ilog2().ok_or(Error::Start)? + octaves >= 31 { - return Err(Error::End); + /// + /// * f_start: initial phase increment + /// * octavaes: number of octaves + /// * samples: number of samples per octave + pub fn new(f_start: i32, octaves: u32, samples: u32) -> Result { + if f_start.checked_ilog2().ok_or(SweepError::Start)? + octaves >= 31 { + return Err(SweepError::End); } - if samples_octave == 0 { - return Err(Error::Length); + if samples == 0 { + return Err(SweepError::Length); } - let a = Float::round( - (Float::exp2((samples_octave as f64).recip()) - 1.0) * (1i64 << 32) as f64, - ) as i32; + let a = Float::round((Float::exp2((samples as f64).recip()) - 1.0) * (1i64 << 32) as f64) + as i32; Ok(Self { sweep: Sweep::new(a, (f_start as i64) << 32), f_end: f_start << octaves, @@ -78,6 +82,7 @@ impl Iterator for SyncExpSweep { } /// Accumulating oscillator +#[derive(Clone, Debug, PartialEq, PartialOrd)] pub struct AccuOsc { f: T, p: i64, From 0d1057670e5aad20f93271173e786e4cecccbe33 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Wed, 30 Oct 2024 01:15:59 +0000 Subject: [PATCH 11/16] sweptsine: use exp_m1 despite not base 2 --- src/sweptsine.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/sweptsine.rs b/src/sweptsine.rs index 4f79fe1..07813e5 100644 --- a/src/sweptsine.rs +++ b/src/sweptsine.rs @@ -1,5 +1,5 @@ use crate::{cossin::cossin, Complex}; -use num_traits::Float; +use num_traits::{Float, FloatConst}; /// Exponential sweep #[derive(Clone, Debug, PartialEq, PartialOrd)] @@ -61,8 +61,8 @@ impl SyncExpSweep { if samples == 0 { return Err(SweepError::Length); } - let a = Float::round((Float::exp2((samples as f64).recip()) - 1.0) * (1i64 << 32) as f64) - as i32; + let a = + Float::round(Float::exp_m1(f64::LN_2() / samples as f64) * (1i64 << 32) as f64) as i32; Ok(Self { sweep: Sweep::new(a, (f_start as i64) << 32), f_end: f_start << octaves, From 2aadc6028112d6cc776322138f5b3adc09c65715 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Wed, 30 Oct 2024 15:42:52 +0000 Subject: [PATCH 12/16] swepsine: fix numerics --- src/sweptsine.rs | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/sweptsine.rs b/src/sweptsine.rs index 07813e5..62e7ae3 100644 --- a/src/sweptsine.rs +++ b/src/sweptsine.rs @@ -58,11 +58,12 @@ impl SyncExpSweep { if f_start.checked_ilog2().ok_or(SweepError::Start)? + octaves >= 31 { return Err(SweepError::End); } - if samples == 0 { + let f0 = f_start as f64 / (1i64 << 32) as f64; + let k = Float::round(f0 * samples as f64 / f64::LN_2()); + if k < 1.0 { return Err(SweepError::Length); } - let a = - Float::round(Float::exp_m1(f64::LN_2() / samples as f64) * (1i64 << 32) as f64) as i32; + let a = Float::round(Float::exp_m1(f0 / k) * (1i64 << 32) as f64) as i32; Ok(Self { sweep: Sweep::new(a, (f_start as i64) << 32), f_end: f_start << octaves, From 2988deff7ccf09cc640b9fb1e376ee6906ec3eb6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Thu, 14 Nov 2024 12:50:13 +0000 Subject: [PATCH 13/16] cic: tweak --- src/cic.rs | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/cic.rs b/src/cic.rs index 94bbbd5..032731b 100644 --- a/src/cic.rs +++ b/src/cic.rs @@ -131,7 +131,7 @@ where /// Optionally ingest a new low-rate sample and /// retrieve the next output. /// - /// A new sample must be supplied at the correct time (see [`Cic::tick()`]) + /// A new sample must be supplied at the correct time (when [`Cic::tick()`] is true) pub fn interpolate(&mut self, x: Option) -> T { if let Some(x) = x { debug_assert_eq!(self.index, 0); @@ -158,7 +158,10 @@ where *i = i.wrapping_add(&x); *i }); - if self.index == 0 { + if let Some(index) = self.index.checked_sub(1) { + self.index = index; + None + } else { self.index = self.rate; let x = self.combs.iter_mut().fold(x, |x, c| { let y = x.wrapping_sub(c); @@ -167,9 +170,6 @@ where }); self.zoh = x; Some(self.zoh) - } else { - self.index -= 1; - None } } } From 86c347b8b55666dde03612162566e0b75cc371c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Mon, 23 Dec 2024 10:28:03 +0000 Subject: [PATCH 14/16] cic: simplify clear, clippy lints --- src/cic.rs | 5 +---- src/hbf.rs | 8 ++++---- 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/src/cic.rs b/src/cic.rs index 032731b..f622862 100644 --- a/src/cic.rs +++ b/src/cic.rs @@ -67,10 +67,7 @@ where /// Zero-initialize the filter state pub fn clear(&mut self) { - self.index = 0; - self.zoh = T::zero(); - self.combs = [T::zero(); N]; - self.integrators = [T::zero(); N]; + *self = Self::new(self.rate); } /// Accepts/provides new slow-rate sample diff --git a/src/hbf.rs b/src/hbf.rs index ad404b5..d5fba19 100644 --- a/src/hbf.rs +++ b/src/hbf.rs @@ -185,8 +185,8 @@ macro_rules! impl_half_i { } impl_half_i!(i8 i16 i32 i64 i128); -impl<'a, T: Copy + Zero + Add + Mul + Sum + Half, const M: usize, const N: usize> Filter - for HbfDec<'a, T, M, N> +impl + Sum + Half, const M: usize, const N: usize> Filter + for HbfDec<'_, T, M, N> { type Item = T; @@ -259,8 +259,8 @@ impl<'a, T: Copy + Zero + Add + Mul + Sum, const M: usize, const N: } } -impl<'a, T: Copy + Zero + Add + Mul + Sum, const M: usize, const N: usize> Filter - for HbfInt<'a, T, M, N> +impl + Sum, const M: usize, const N: usize> Filter + for HbfInt<'_, T, M, N> { type Item = T; From ee9b18185e37cda19b8ca9c88eaa5f09e2cd08ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Mon, 23 Dec 2024 10:35:08 +0000 Subject: [PATCH 15/16] testing: doc and simplify --- src/testing.rs | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/src/testing.rs b/src/testing.rs index 7660896..88bcb00 100644 --- a/src/testing.rs +++ b/src/testing.rs @@ -1,5 +1,7 @@ +//! Tools to test and benchmark algorithms #![allow(dead_code)] use super::Complex; +use num_traits::Float; /// Maximum acceptable error between a computed and actual value given fixed and relative /// tolerances. @@ -13,24 +15,28 @@ use super::Complex; /// /// # Returns /// Maximum acceptable error. -pub fn max_error(a: f64, b: f64, rtol: f64, atol: f64) -> f64 { +pub fn max_error(a: T, b: T, rtol: T, atol: T) -> T { rtol * a.abs().max(b.abs()) + atol } -pub fn isclose(a: f64, b: f64, rtol: f64, atol: f64) -> bool { - (a - b).abs() <= a.abs().max(b.abs()) * rtol + atol +/// Return whether two numbers are within absolute plus relative tolerance +pub fn isclose(a: T, b: T, rtol: T, atol: T) -> bool { + (a - b).abs() <= max_error(a, b, rtol, atol) } -pub fn isclosef(a: f32, b: f32, rtol: f32, atol: f32) -> bool { - (a - b).abs() <= a.abs().max(b.abs()) * rtol + atol +/// Return whether all values are close +pub fn allclose(a: &[T], b: &[T], rtol: T, atol: T) -> bool { + a.iter().zip(b).all(|(a, b)| isclose(*a, *b, rtol, atol)) } -pub fn complex_isclose(a: Complex, b: Complex, rtol: f32, atol: f32) -> bool { - isclosef(a.re, b.re, rtol, atol) && isclosef(a.im, b.im, rtol, atol) +/// Return whether both real and imaginary component are close +pub fn complex_isclose(a: Complex, b: Complex, rtol: T, atol: T) -> bool { + isclose(a.re, b.re, rtol, atol) && isclose(a.im, b.im, rtol, atol) } -pub fn complex_allclose(a: &[Complex], b: &[Complex], rtol: f32, atol: f32) -> bool { +/// Return whether all complex values are close +pub fn complex_allclose(a: &[Complex], b: &[Complex], rtol: T, atol: T) -> bool { a.iter() .zip(b) - .all(|(&i, &j)| complex_isclose(i, j, rtol, atol)) + .all(|(a, b)| complex_isclose(*a, *b, rtol, atol)) } From 7bcc9bb0e894f7b62d9632db675499f9bcb18582 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Mon, 23 Dec 2024 10:55:21 +0000 Subject: [PATCH 16/16] sweptsine: docs, style, comment --- src/cic.rs | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/cic.rs b/src/cic.rs index f622862..73e40a0 100644 --- a/src/cic.rs +++ b/src/cic.rs @@ -14,8 +14,8 @@ pub struct Cic { /// Up/downsampler state (count down) index: u32, /// Zero order hold behind comb sections. - /// Interpolator: In the middle combined with the upsampler - /// Decimator: After combs to support `get_decimate()` + /// Interpolator: Combined with the upsampler + /// Decimator: To support `get_decimate()` zoh: T, /// Comb/differentiator state combs: [T; N], @@ -80,7 +80,7 @@ where /// Current interpolator output pub fn get_interpolate(&self) -> T { - *self.integrators.last().unwrap_or(&self.zoh) + self.integrators.last().copied().unwrap_or(self.zoh) } /// Current decimator output @@ -109,7 +109,11 @@ where /// Establish a settled filter state pub fn settle_interpolate(&mut self, x: T) { self.clear(); - *self.combs.first_mut().unwrap_or(&mut self.zoh) = x; + if let Some(c) = self.combs.first_mut() { + *c = x; + } else { + self.zoh = x; + } let g = self.gain(); if let Some(i) = self.integrators.last_mut() { *i = x * g; @@ -143,6 +147,7 @@ where self.index -= 1; } self.integrators.iter_mut().fold(self.zoh, |x, i| { + // Overflow is not OK *i += x; *i })