Updated test code to also pass on f32 compilation. Improved doc tests.

This commit is contained in:
Anne de Jong 2024-08-15 13:10:40 +02:00
parent 02bb55a721
commit 302c2cbfd3
7 changed files with 88 additions and 55 deletions

View File

@ -86,6 +86,8 @@ realfft = "3.3.0"
# Fast Mutex # Fast Mutex
parking_lot = "0.12.3" parking_lot = "0.12.3"
# Builder code
derive_builder = "0.20.0" derive_builder = "0.20.0"
[dev-dependencies] [dev-dependencies]
@ -93,6 +95,9 @@ ndarray-rand = "0.14.0"
[features] [features]
default = ["f64", "cpal-api", "record"] 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 # Use this for debugging extensions
# default = ["f64", "python-bindings", "record", "cpal-api"] # default = ["f64", "python-bindings", "record", "cpal-api"]

View File

@ -10,6 +10,8 @@ cfg_if::cfg_if! {
pub type Flt = f64; pub type Flt = f64;
/// Ratio between circumference and diameter of a circle /// Ratio between circumference and diameter of a circle
pub const pi: Flt = std::f64::consts::PI; pub const pi: Flt = std::f64::consts::PI;
} }
else if #[cfg(feature="f32")] { else if #[cfg(feature="f32")] {
/// Floating-point value, compile time option to make it either f32, or f64 /// Floating-point value, compile time option to make it either f32, or f64

View File

@ -221,9 +221,11 @@ impl Biquad {
/// ///
/// The analog filter is defined as: /// The analog filter is defined as:
/// ///
/// ```math
/// b0 + b1*s + b2*s^2 /// b0 + b1*s + b2*s^2
/// H(s) = -------------------- /// H(s) = --------------------
/// a0 + a1*s + a2*s^2 /// a0 + a1*s + a2*s^2
/// ```
/// ///
/// # Args /// # Args
/// ///
@ -258,8 +260,6 @@ impl Biquad {
let Ksq = K.powi(2); let Ksq = K.powi(2);
// //
let a0fac = a2a * Ksq + a1a * K + a0a; let a0fac = a2a * Ksq + a1a * K + a0a;
println!("Ksq = {Ksq}");
println!("a0fac = {a0fac}");
// Coefficient b0 // Coefficient b0
let b0 = (b2a * Ksq + b1a * K + b0a) / a0fac; let b0 = (b2a * Ksq + b1a * K + b0a) / a0fac;
// Coefficient b1 // Coefficient b1
@ -358,8 +358,6 @@ impl Biquad {
} else { } else {
[1., 0., 0.] [1., 0., 0.]
}; };
println!("b = {b:?}");
println!("a = {a:?}");
Biquad::bilinear(fs, &b, &a, fwarp) Biquad::bilinear(fs, &b, &a, fwarp)
} }
} }
@ -418,14 +416,16 @@ mod test {
let fs = 1e5; let fs = 1e5;
let fc = 10.; let fc = 10.;
let b = Biquad::firstOrderMovingAverage(fs, fc).unwrap(); let b = Biquad::firstOrderMovingAverage(fs, fc).unwrap();
println!("b = {b:#?}");
let mut freq = Dcol::from_elem(5, 0.); let mut freq = Dcol::from_elem(5, 0.);
freq[1] = fc; freq[1] = fc;
freq[2] = fs / 2.; freq[2] = fs / 2.;
let tf = b.tf(fs, freq.view()); let tf = b.tf(fs, freq.view());
// println!("{:?}", tf);
assert_abs_diff_eq!(tf[0].re, 1., epsilon = 1e-6); let epsilon = Flt::EPSILON * 150.;
assert_abs_diff_eq!(tf[0].im, 0.); assert_abs_diff_eq!(tf[0].re, 1., epsilon = 10. * epsilon);
assert_abs_diff_eq!(tf[1].abs(), 1. / Flt::sqrt(2.), epsilon = 1e-6); 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] #[test]
fn test_bilinear() { fn test_bilinear() {
@ -436,11 +436,21 @@ mod test {
let b2 = Biquad::bilinear_zpk(fs, None, Some(PoleOrZero::Real1(-omgc)), Some(omgc), None); let b2 = Biquad::bilinear_zpk(fs, None, Some(PoleOrZero::Real1(-omgc)), Some(omgc), None);
println!("b1 = {b1:?}"); println!("b1 = {b1:?}");
println!("b2 = {b2:?}"); println!("b2 = {b2:?}");
assert_abs_diff_eq!((b1.tf(fs, &[0.])[0] - Cflt::ONE).abs(), 0., epsilon = 1e-9); let epsilon = Flt::EPSILON * 10.;
assert_abs_diff_eq!((b2.tf(fs, &[0.])[0] - Cflt::ONE).abs(), 0., epsilon = 1e-9); 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_eq!(b1, b2);
assert_abs_diff_eq!((b1.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 = 1e-9); assert_abs_diff_eq!((b2.tf(fs, &[fs / 2.])[0]).abs(), 0., epsilon = epsilon);
} }
#[test] #[test]
fn test_firstOrderHighPass() { fn test_firstOrderHighPass() {
@ -449,9 +459,14 @@ mod test {
let b3 = Biquad::firstOrderHighPass(fs, fc).unwrap(); let b3 = Biquad::firstOrderHighPass(fs, fc).unwrap();
println!("b3 = {b3:?}"); println!("b3 = {b3:?}");
assert_abs_diff_eq!((b3.tf(fs, &[0.])[0]).abs(), 0., epsilon = 1e-9); let epsilon = Flt::EPSILON * 10.;
assert_abs_diff_eq!((b3.tf(fs, &[(fs-fs/1e9) / 2.])[0]).abs(), 1., epsilon = 1e-9); assert_abs_diff_eq!((b3.tf(fs, &[0.])[0]).abs(), 0., epsilon = epsilon);
assert_abs_diff_eq!((b3.tf(fs, &[fc])[0]).abs(), (0.5).sqrt(), epsilon = 1e-9); 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.]; // let freq = &[0., 10.,100.,1000., 2000.];
// println!("{:?}", b3.tf(fs, freq)); // println!("{:?}", b3.tf(fs, freq));
} }

View File

@ -3,11 +3,13 @@
//! is a fine source to understand the theory presented here. //! is a fine source to understand the theory presented here.
//! A Butterworth lowpass filter has the form //! A Butterworth lowpass filter has the form
//! //!
//! ```math
//! 1 //! 1
//! |H(s)|^2 = ------------------ //! |H(s)|^2 = ------------------
//! 1+ (omega/omega_c)^(2n) //! 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; use approx::abs_diff_eq;

View File

@ -30,6 +30,9 @@ pub struct ApsSettings {
impl ApsSettingsBuilder { impl ApsSettingsBuilder {
fn validate(&self) -> Result<()> { fn validate(&self) -> Result<()> {
if !self.fs.is_some() {
bail!("Sampling frequency not given");
}
let fs = self.fs.unwrap(); let fs = self.fs.unwrap();
if !fs.is_normal() { if !fs.is_normal() {

View File

@ -52,7 +52,6 @@ impl RtAps {
let mut meta: Option<Arc<StreamMetaData>> = None; let mut meta: Option<Arc<StreamMetaData>> = None;
if let Some(msg) = rx.recv_timeout(std::time::Duration::from_millis(10)).ok() { if let Some(msg) = rx.recv_timeout(std::time::Duration::from_millis(10)).ok() {
match msg { match msg {
InStreamMsg::StreamStarted(new_meta) => { InStreamMsg::StreamStarted(new_meta) => {
aps.reset(); aps.reset();
@ -96,7 +95,6 @@ impl RtAps {
} }
// Move last_cps into mutex. // Move last_cps into mutex.
if let Some(last_cps) = last_cps.take() { if let Some(last_cps) = last_cps.take() {
*last_result_lock = Some(RtApsComm::NewResult(last_cps)); *last_result_lock = Some(RtApsComm::NewResult(last_cps));
} }
} }
@ -134,15 +132,25 @@ impl Drop for RtAps {
mod test { mod test {
use std::time::Duration; use std::time::Duration;
use anyhow::{anyhow, bail, Result};
use super::*; use super::*;
use crate::{daq::StreamMgr, ps::ApsSettingsBuilder}; use crate::{daq::StreamMgr, ps::ApsSettingsBuilder};
#[test] #[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 mut smgr = StreamMgr::new();
let rtaps = RtAps::new(&mut smgr, settings);
smgr.startDefaultInputStream()?; 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)); thread::sleep(Duration::from_secs(2));
drop(rtaps); drop(rtaps);
} }

View File

@ -217,10 +217,7 @@ impl Siggen {
/// Set the DC offset for all channels /// Set the DC offset for all channels
pub fn setDCOffset(&mut self, dc: &[Flt]) { pub fn setDCOffset(&mut self, dc: &[Flt]) {
self.channels self.channels.iter_mut().zip(dc).for_each(|(ch, dc)| {
.iter_mut()
.zip(dc)
.for_each(|(ch, dc)| {
ch.DCOffset = *dc; ch.DCOffset = *dc;
}); });
} }
@ -244,7 +241,9 @@ impl Siggen {
/// Creates *interleaved* output signal /// Creates *interleaved* output signal
pub fn genSignal<T>(&mut self, out: &mut [T]) pub fn genSignal<T>(&mut self, out: &mut [T])
where T: Sample + FromSample<Flt> + Debug, Flt: Sample where
T: Sample + FromSample<Flt> + Debug,
Flt: Sample,
{ {
let nch = self.nchannels(); let nch = self.nchannels();
let nsamples: usize = out.len() / nch; let nsamples: usize = out.len() / nch;
@ -252,17 +251,20 @@ impl Siggen {
// Create source signal // Create source signal
self.source_buf.resize(nsamples, 0.0); 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); // println!("Source signal: {:?}", self.source_buf);
// Write output while casted to the correct type // Write output while casted to the correct type
// Iterate over each channel, and counter // Iterate over each channel, and counter
self.chout_buf.resize(nch, vec![]); self.chout_buf.resize(nch, vec![]);
for (channelno, (channel, chout)) in self.channels for (channelno, (channel, chout)) in self
.channels
.iter_mut() .iter_mut()
.zip(self.chout_buf.iter_mut()) .zip(self.chout_buf.iter_mut())
.enumerate() { .enumerate()
{
chout.resize(nsamples, 0.0); chout.resize(nsamples, 0.0);
// Create output signal, overwrite chout // Create output signal, overwrite chout
@ -298,10 +300,7 @@ impl Siggen {
/// as number of channels in signal generator. /// as number of channels in signal generator.
pub fn setMute(&mut self, mute: &[bool]) { pub fn setMute(&mut self, mute: &[bool]) {
assert!(mute.len() == self.nchannels()); assert!(mute.len() == self.nchannels());
self.channels self.channels.iter_mut().zip(mute).for_each(|(s, m)| {
.iter_mut()
.zip(mute)
.for_each(|(s, m)| {
s.setMute(*m); s.setMute(*m);
}); });
} }
@ -386,7 +385,10 @@ impl SiggenChannelConfig {
#[cfg(test)] #[cfg(test)]
mod test { mod test {
use approx::assert_abs_diff_eq;
use super::*; use super::*;
use crate::Flt;
#[test] #[test]
fn test_whitenoise() { fn test_whitenoise() {
@ -408,14 +410,15 @@ mod test {
siggen.reset(10.0); siggen.reset(10.0);
siggen.setAllMute(false); siggen.setAllMute(false);
siggen.genSignal(&mut s1); siggen.genSignal(&mut s1);
siggen.reset(10.0);
siggen.genSignal(&mut s2); siggen.genSignal(&mut s2);
let absdiff = s1 let absdiff = s1
.iter() .iter()
.zip(s2.iter()) .zip(s2.iter())
.map(|(s1, s2)| { Flt::abs(*s1 - *s2) }) .map(|(s1, s2)| Flt::abs(*s1 - *s2))
.sum::<Flt>(); .sum::<Flt>();
assert!(absdiff < 1e-10); assert_abs_diff_eq!(absdiff, 0., epsilon = Flt::EPSILON * 100.);
} }
#[test] #[test]
@ -438,23 +441,18 @@ mod test {
siggen.genSignal(&mut signal[Nframes / 2..]); siggen.genSignal(&mut signal[Nframes / 2..]);
// Mean square of the signal // Mean square of the signal
let ms1 = let ms1 = signal.iter().step_by(2).map(|s1| *s1 * *s1).sum::<Flt>() / (Nframes as Flt);
signal
.iter()
.step_by(2)
.map(|s1| { *s1 * *s1 })
.sum::<Flt>() / (Nframes as Flt);
println!("ms1: {}", ms1); println!("ms1: {}", ms1);
let ms2 = let ms2 = signal
signal
.iter() .iter()
.skip(1) .skip(1)
.step_by(2) .step_by(2)
.map(|s1| { *s1 * *s1 }) .map(|s1| *s1 * *s1)
.sum::<Flt>() / (Nframes as Flt); .sum::<Flt>()
/ (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); assert_eq!(ms2, 0.0);
} }