From 302c2cbfd38cad02730ed776a8396c6807c2d0d9 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Thu, 15 Aug 2024 13:10:40 +0200 Subject: [PATCH] Updated test code to also pass on f32 compilation. Improved doc tests. --- Cargo.toml | 5 ++++ src/config.rs | 2 ++ src/filter/biquad.rs | 45 ++++++++++++++++++++---------- src/filter/butter.rs | 4 ++- src/ps/aps.rs | 3 ++ src/rt/rtaps.rs | 18 ++++++++---- src/siggen.rs | 66 +++++++++++++++++++++----------------------- 7 files changed, 88 insertions(+), 55 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 46a9690..b36a913 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -86,6 +86,8 @@ realfft = "3.3.0" # Fast Mutex parking_lot = "0.12.3" + +# Builder code derive_builder = "0.20.0" [dev-dependencies] @@ -93,6 +95,9 @@ ndarray-rand = "0.14.0" [features] default = ["f64", "cpal-api", "record"] +# Use this to test if everything works well in f32 +# default = ["f32", "cpal-api", "record"] + # Use this for debugging extensions # default = ["f64", "python-bindings", "record", "cpal-api"] diff --git a/src/config.rs b/src/config.rs index fb4735b..749bd22 100644 --- a/src/config.rs +++ b/src/config.rs @@ -10,6 +10,8 @@ cfg_if::cfg_if! { pub type Flt = f64; /// Ratio between circumference and diameter of a circle pub const pi: Flt = std::f64::consts::PI; + + } else if #[cfg(feature="f32")] { /// Floating-point value, compile time option to make it either f32, or f64 diff --git a/src/filter/biquad.rs b/src/filter/biquad.rs index 5ef06e7..c7816de 100644 --- a/src/filter/biquad.rs +++ b/src/filter/biquad.rs @@ -221,9 +221,11 @@ impl Biquad { /// /// The analog filter is defined as: /// + /// ```math /// b0 + b1*s + b2*s^2 /// H(s) = -------------------- /// a0 + a1*s + a2*s^2 + /// ``` /// /// # Args /// @@ -258,8 +260,6 @@ impl Biquad { let Ksq = K.powi(2); // let a0fac = a2a * Ksq + a1a * K + a0a; - println!("Ksq = {Ksq}"); - println!("a0fac = {a0fac}"); // Coefficient b0 let b0 = (b2a * Ksq + b1a * K + b0a) / a0fac; // Coefficient b1 @@ -358,8 +358,6 @@ impl Biquad { } else { [1., 0., 0.] }; - println!("b = {b:?}"); - println!("a = {a:?}"); Biquad::bilinear(fs, &b, &a, fwarp) } } @@ -418,14 +416,16 @@ mod test { let fs = 1e5; let fc = 10.; let b = Biquad::firstOrderMovingAverage(fs, fc).unwrap(); + println!("b = {b:#?}"); let mut freq = Dcol::from_elem(5, 0.); freq[1] = fc; freq[2] = fs / 2.; let tf = b.tf(fs, freq.view()); - // println!("{:?}", tf); - assert_abs_diff_eq!(tf[0].re, 1., epsilon = 1e-6); - assert_abs_diff_eq!(tf[0].im, 0.); - assert_abs_diff_eq!(tf[1].abs(), 1. / Flt::sqrt(2.), epsilon = 1e-6); + + let epsilon = Flt::EPSILON * 150.; + assert_abs_diff_eq!(tf[0].re, 1., epsilon = 10. * epsilon); + assert_abs_diff_eq!(tf[0].im, 0., epsilon = epsilon); + assert_abs_diff_eq!(tf[1].abs(), 1. / Flt::sqrt(2.), epsilon = 1e-4); } #[test] fn test_bilinear() { @@ -436,11 +436,21 @@ mod test { let b2 = Biquad::bilinear_zpk(fs, None, Some(PoleOrZero::Real1(-omgc)), Some(omgc), None); println!("b1 = {b1:?}"); println!("b2 = {b2:?}"); - assert_abs_diff_eq!((b1.tf(fs, &[0.])[0] - Cflt::ONE).abs(), 0., epsilon = 1e-9); - assert_abs_diff_eq!((b2.tf(fs, &[0.])[0] - Cflt::ONE).abs(), 0., epsilon = 1e-9); + let epsilon = Flt::EPSILON * 10.; + println!("Epsilon = {epsilon}"); + assert_abs_diff_eq!( + (b1.tf(fs, &[0.])[0] - Cflt::ONE).abs(), + 0., + epsilon = epsilon + ); + assert_abs_diff_eq!( + (b2.tf(fs, &[0.])[0] - Cflt::ONE).abs(), + 0., + epsilon = epsilon + ); // assert_eq!(b1, b2); - assert_abs_diff_eq!((b1.tf(fs, &[fs / 2.])[0]).abs(), 0., epsilon = 1e-9); - assert_abs_diff_eq!((b2.tf(fs, &[fs / 2.])[0]).abs(), 0., epsilon = 1e-9); + assert_abs_diff_eq!((b1.tf(fs, &[fs / 2.])[0]).abs(), 0., epsilon = epsilon); + assert_abs_diff_eq!((b2.tf(fs, &[fs / 2.])[0]).abs(), 0., epsilon = epsilon); } #[test] fn test_firstOrderHighPass() { @@ -449,9 +459,14 @@ mod test { let b3 = Biquad::firstOrderHighPass(fs, fc).unwrap(); println!("b3 = {b3:?}"); - assert_abs_diff_eq!((b3.tf(fs, &[0.])[0]).abs(), 0., epsilon = 1e-9); - assert_abs_diff_eq!((b3.tf(fs, &[(fs-fs/1e9) / 2.])[0]).abs(), 1., epsilon = 1e-9); - assert_abs_diff_eq!((b3.tf(fs, &[fc])[0]).abs(), (0.5).sqrt(), epsilon = 1e-9); + let epsilon = Flt::EPSILON * 10.; + assert_abs_diff_eq!((b3.tf(fs, &[0.])[0]).abs(), 0., epsilon = epsilon); + assert_abs_diff_eq!( + (b3.tf(fs, &[(fs - fs / 1e9) / 2.])[0]).abs(), + 1., + epsilon = epsilon + ); + assert_abs_diff_eq!((b3.tf(fs, &[fc])[0]).abs(), (0.5).sqrt(), epsilon = epsilon); // let freq = &[0., 10.,100.,1000., 2000.]; // println!("{:?}", b3.tf(fs, freq)); } diff --git a/src/filter/butter.rs b/src/filter/butter.rs index 91e71e3..896c54b 100644 --- a/src/filter/butter.rs +++ b/src/filter/butter.rs @@ -3,11 +3,13 @@ //! is a fine source to understand the theory presented here. //! A Butterworth lowpass filter has the form //! +//! ```math //! 1 //! |H(s)|^2 = ------------------ //! 1+ (omega/omega_c)^(2n) +//! ``` //! -//! where n is the order of the filter. +//! where `n`` is the order of the filter. use approx::abs_diff_eq; diff --git a/src/ps/aps.rs b/src/ps/aps.rs index 7a98bd7..3f2bdd8 100644 --- a/src/ps/aps.rs +++ b/src/ps/aps.rs @@ -30,6 +30,9 @@ pub struct ApsSettings { impl ApsSettingsBuilder { fn validate(&self) -> Result<()> { + if !self.fs.is_some() { + bail!("Sampling frequency not given"); + } let fs = self.fs.unwrap(); if !fs.is_normal() { diff --git a/src/rt/rtaps.rs b/src/rt/rtaps.rs index cc08e9c..f1983e2 100644 --- a/src/rt/rtaps.rs +++ b/src/rt/rtaps.rs @@ -52,7 +52,6 @@ impl RtAps { let mut meta: Option> = None; if let Some(msg) = rx.recv_timeout(std::time::Duration::from_millis(10)).ok() { - match msg { InStreamMsg::StreamStarted(new_meta) => { aps.reset(); @@ -96,7 +95,6 @@ impl RtAps { } // Move last_cps into mutex. if let Some(last_cps) = last_cps.take() { - *last_result_lock = Some(RtApsComm::NewResult(last_cps)); } } @@ -134,15 +132,25 @@ impl Drop for RtAps { mod test { use std::time::Duration; + use anyhow::{anyhow, bail, Result}; + use super::*; use crate::{daq::StreamMgr, ps::ApsSettingsBuilder}; #[test] - fn test_rtaps1() -> Result<()> { + fn test_rtaps1() -> Result<(), anyhow::Error> { { - let settings = ApsSettingsBuilder::default().nfft(2048).build().unwrap(); let mut smgr = StreamMgr::new(); - let rtaps = RtAps::new(&mut smgr, settings); smgr.startDefaultInputStream()?; + let meta = smgr + .getStreamMetaData(crate::daq::StreamType::Input) + .ok_or_else(|| anyhow!("Stream is not running"))?; + + let settings = ApsSettingsBuilder::default() + .nfft(2048) + .fs(meta.samplerate) + .build() + .unwrap(); + let rtaps = RtAps::new(&mut smgr, settings); thread::sleep(Duration::from_secs(2)); drop(rtaps); } diff --git a/src/siggen.rs b/src/siggen.rs index f6c4d81..324d043 100644 --- a/src/siggen.rs +++ b/src/siggen.rs @@ -22,7 +22,7 @@ //! ``` use super::config::*; use super::filter::Filter; -use dasp_sample::{ FromSample, Sample }; +use dasp_sample::{FromSample, Sample}; use rayon::prelude::*; use std::fmt::Debug; use std::iter::ExactSizeIterator; @@ -154,7 +154,7 @@ pub struct Siggen { // Output buffers (for filtered source signal) chout_buf: Vec>, } -#[cfg(feature="python-bindings")] +#[cfg(feature = "python-bindings")] #[cfg_attr(feature = "python-bindings", pymethods)] impl Siggen { #[pyo3(name = "newWhiteNoise")] @@ -217,12 +217,9 @@ impl Siggen { /// Set the DC offset for all channels pub fn setDCOffset(&mut self, dc: &[Flt]) { - self.channels - .iter_mut() - .zip(dc) - .for_each(|(ch, dc)| { - ch.DCOffset = *dc; - }); + self.channels.iter_mut().zip(dc).for_each(|(ch, dc)| { + ch.DCOffset = *dc; + }); } /// Create a sine wave signal generator @@ -244,7 +241,9 @@ impl Siggen { /// Creates *interleaved* output signal pub fn genSignal(&mut self, out: &mut [T]) - where T: Sample + FromSample + Debug, Flt: Sample + where + T: Sample + FromSample + Debug, + Flt: Sample, { let nch = self.nchannels(); let nsamples: usize = out.len() / nch; @@ -252,17 +251,20 @@ impl Siggen { // Create source signal self.source_buf.resize(nsamples, 0.0); - self.source.genSignal_unscaled(&mut self.source_buf.iter_mut()); + self.source + .genSignal_unscaled(&mut self.source_buf.iter_mut()); // println!("Source signal: {:?}", self.source_buf); // Write output while casted to the correct type // Iterate over each channel, and counter self.chout_buf.resize(nch, vec![]); - for (channelno, (channel, chout)) in self.channels + for (channelno, (channel, chout)) in self + .channels .iter_mut() .zip(self.chout_buf.iter_mut()) - .enumerate() { + .enumerate() + { chout.resize(nsamples, 0.0); // Create output signal, overwrite chout @@ -298,12 +300,9 @@ impl Siggen { /// as number of channels in signal generator. pub fn setMute(&mut self, mute: &[bool]) { assert!(mute.len() == self.nchannels()); - self.channels - .iter_mut() - .zip(mute) - .for_each(|(s, m)| { - s.setMute(*m); - }); + self.channels.iter_mut().zip(mute).for_each(|(s, m)| { + s.setMute(*m); + }); } } @@ -386,7 +385,10 @@ impl SiggenChannelConfig { #[cfg(test)] mod test { + use approx::assert_abs_diff_eq; + use super::*; + use crate::Flt; #[test] fn test_whitenoise() { @@ -408,14 +410,15 @@ mod test { siggen.reset(10.0); siggen.setAllMute(false); siggen.genSignal(&mut s1); + siggen.reset(10.0); siggen.genSignal(&mut s2); let absdiff = s1 .iter() .zip(s2.iter()) - .map(|(s1, s2)| { Flt::abs(*s1 - *s2) }) + .map(|(s1, s2)| Flt::abs(*s1 - *s2)) .sum::(); - assert!(absdiff < 1e-10); + assert_abs_diff_eq!(absdiff, 0., epsilon = Flt::EPSILON * 100.); } #[test] @@ -438,23 +441,18 @@ mod test { siggen.genSignal(&mut signal[Nframes / 2..]); // Mean square of the signal - let ms1 = - signal - .iter() - .step_by(2) - .map(|s1| { *s1 * *s1 }) - .sum::() / (Nframes as Flt); + let ms1 = signal.iter().step_by(2).map(|s1| *s1 * *s1).sum::() / (Nframes as Flt); println!("ms1: {}", ms1); - let ms2 = - signal - .iter() - .skip(1) - .step_by(2) - .map(|s1| { *s1 * *s1 }) - .sum::() / (Nframes as Flt); + let ms2 = signal + .iter() + .skip(1) + .step_by(2) + .map(|s1| *s1 * *s1) + .sum::() + / (Nframes as Flt); - assert!(Flt::abs(ms1 - 0.5) < 1e-12); + assert_abs_diff_eq!(Flt::abs(ms1 - 0.5) , 0., epsilon= Flt::EPSILON * 1e3); assert_eq!(ms2, 0.0); }