From 66ab149159436127a960e89874f4189b11694713 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Thu, 11 Jul 2024 21:02:21 +0200 Subject: [PATCH] Removed unnecessary complex API with APSResult and switched to two methods: compute_last and compute_all. Improved docstrings and fixes in docstrings. --- README.md | 2 +- src/bin/lasp_record.rs | 4 +- src/filter/biquad.rs | 4 +- src/filter/mod.rs | 2 +- src/ps/aps.rs | 105 +++++++++++++++++++++-------------------- src/ps/mod.rs | 4 +- src/ps/ps.rs | 8 ++-- 7 files changed, 66 insertions(+), 63 deletions(-) diff --git a/README.md b/README.md index 5fe1656..ff8714b 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ and processing of (multi) sensor data in real time on a PC and output results. ## Documentation -Documentation is provided at [doc.rs](https://docs.rs/lasprs/0.2.1/lasprs). +Documentation is provided at [doc.rs](https://docs.rs/lasprs/latest/lasprs). ## Python bindings diff --git a/src/bin/lasp_record.rs b/src/bin/lasp_record.rs index dae556a..7234cc2 100644 --- a/src/bin/lasp_record.rs +++ b/src/bin/lasp_record.rs @@ -15,11 +15,11 @@ struct Cli { /// File name to write recording to filename: String, - /// Recording duration in [s]. Rounds down to whole seconds. If not specified, records until user presses a key + /// Recording duration in \[s\]. Rounds down to whole seconds. If not specified, records until user presses a key #[arg(short, long = "duration", default_value_t = 0.)] duration_s: Flt, - /// Start delay in [s]. Rounds down to whole seconds. If not specified, no + /// Start delay in \[s\]. Rounds down to whole seconds. If not specified, no /// start delay will be used. #[arg(short, long = "startdelay", default_value_t = 0.)] start_delay_s: Flt, diff --git a/src/filter/biquad.rs b/src/filter/biquad.rs index c37415c..a86c396 100644 --- a/src/filter/biquad.rs +++ b/src/filter/biquad.rs @@ -173,8 +173,8 @@ impl Biquad { /// First order low pass filter (one pole in the real axis). No pre-warping /// correction done. /// - /// * `fs` - Sampling frequency [Hz] - /// * `fc` - Cut-off frequency (-3 dB point) [Hz] + /// * `fs` - Sampling frequency \[Hz\] + /// * `fc` - Cut-off frequency (-3 dB point) \[Hz\] pub fn firstOrderLowPass(fs: Flt, fc: Flt) -> Result { if fc <= 0. { bail!("Cuton frequency, given: should be > 0") diff --git a/src/filter/mod.rs b/src/filter/mod.rs index 6ddd075..da9e5e1 100644 --- a/src/filter/mod.rs +++ b/src/filter/mod.rs @@ -42,7 +42,7 @@ where /// /// # Args /// - /// * `freq` - The frequency in [Hz] + /// * `freq` - The frequency in \[Hz\] /// /// # Returns /// diff --git a/src/ps/aps.rs b/src/ps/aps.rs index b612ead..4ce6fb4 100644 --- a/src/ps/aps.rs +++ b/src/ps/aps.rs @@ -22,15 +22,15 @@ impl Default for Overlap { } } -/// Result from [AvPowerspectra.compute()] -pub enum ApsResult<'a> { +/// Result from [compute method](AvPowerSpectra::compute). +enum ApsResult<'a> { /// Returns all intermediate results. Useful when evolution over time should /// be visible. I.e. in case of spectrograms, but also in case the /// converging process to the average should be made visible. - AllIntermediateResults(Vec), + AllIntermediateResults(Vec), /// Give only last result back, the most recent value. - OnlyLastResult(&'a CPS), + OnlyLastResult(&'a CPSResult), /// No new data available. Nothing given back. None, @@ -47,7 +47,7 @@ pub enum ApsMode { /// where new data is weighted with old data, and old data exponentially /// backs off. This mode only makes sense when `tau >> nfft/fs` ExponentialWeighting { - /// Sampling frequency in [Hz], used for computing IIR filter coefficient. + /// Sampling frequency in `[Hz]`, used for computing IIR filter coefficient. fs: Flt, /// Time weighting constant, follows convention of Sound Level Meters. /// Means the data is approximately low-pass filtered with a cut-off @@ -88,7 +88,7 @@ pub struct AvPowerSpectra { timebuf: TimeBuffer, // Current estimation of the power spectra - cur_est: CPS, + cur_est: CPSResult, } impl AvPowerSpectra { /// The FFT Length of estimating (cross)power spectra @@ -125,7 +125,7 @@ impl AvPowerSpectra { pub fn reset(&mut self) { self.N = 0; self.timebuf.reset(); - self.cur_est = CPS::zeros((0,0,0)); + self.cur_est = CPSResult::zeros((0, 0, 0)); } /// Create new averaged power spectra estimator for weighing over the full /// amount of data supplied (no exponential spectra weighting) using @@ -140,7 +140,7 @@ impl AvPowerSpectra { /// # Panics /// /// - When nfft is not even, or 0. - /// + /// pub fn new_simple(nfft: usize) -> AvPowerSpectra { AvPowerSpectra::build(nfft, None, None, None).unwrap() } @@ -199,7 +199,7 @@ impl AvPowerSpectra { overlap_keep, mode, N: 0, - cur_est: CPS::default((0, 0, 0)), + cur_est: CPSResult::default((0, 0, 0)), timebuf: TimeBuffer::new(), }) } @@ -215,7 +215,7 @@ impl AvPowerSpectra { // Initialize to zero if self.cur_est.len() == 0 { assert_eq!(self.N, 0); - self.cur_est = CPS::zeros(Cpsnew.raw_dim().f()); + self.cur_est = CPSResult::zeros(Cpsnew.raw_dim().f()); } // Update the number of blocks processed @@ -285,27 +285,18 @@ impl AvPowerSpectra { /// /// * `timedata``: New available time data. Number of columns is number of /// channels, number of rows is number of frames (samples per channel). - /// * `giveInterMediateResults` - If true: returns a vector of all - /// intermediate results. Useful when plotting spectra over time. /// /// # Panics /// /// If timedata.ncols() does not match number of columns in already present /// data. - pub fn compute<'a, 'b, T>( - &'a mut self, - timedata: T, - giveInterMediateResults: bool, - ) -> ApsResult<'a> + pub fn compute_last<'a, 'b, T>(&'a mut self, timedata: T) -> Option<&'a CPSResult> where T: AsArray<'b, Flt, Ix2>, { // Push new data in the time buffer. self.timebuf.push(timedata); - // Storage for the result - let mut result = ApsResult::None; - // Flag to indicate that we have obtained one result for sure. let mut computed_single = false; @@ -314,34 +305,46 @@ impl AvPowerSpectra { // Compute cross-power spectra for current time block self.update_singleblock(&timeblock); - // We ha computed_single = true; - if giveInterMediateResults { - // Initialize with empty vector - if let ApsResult::None = result { - result = ApsResult::AllIntermediateResults(Vec::new()) - } - } - - // So - if let ApsResult::AllIntermediateResults(v) = &mut result { - v.push(self.cur_est.clone()); - } } - - if computed_single && !giveInterMediateResults { - // We have computed it once, but we are not interested in - // intermediate results. So we return a reference to the last data. - return ApsResult::OnlyLastResult(&self.cur_est); + if computed_single { + return Some(&self.cur_est); } - result + None } - /// See [AvPowerSpectra](compute()) - pub fn compute_last<'a, T>(&mut self, timedata: T) -> ApsResult + + /// Computes average (cross)power spectra, and returns all intermediate + /// estimates that can be calculated. This is useful when plotting spectra + /// as a function of time, and intermediate results need also be plotted. + /// + /// # Args + /// + /// * `timedata``: New available time data. Number of columns is number of + /// channels, number of rows is number of frames (samples per channel). + /// + /// # Panics + /// + /// If timedata.ncols() does not match number of columns in already present + /// data. + pub fn compute_all<'a, 'b, T>(&'a mut self, timedata: T) -> Vec where - T: AsArray<'a, Flt, Ix2>, + T: AsArray<'b, Flt, Ix2>, { - return self.compute(timedata, false); + // Push new data in the time buffer. + self.timebuf.push(timedata); + + // Storage for the result + let mut result = Vec::new(); + + // Iterate over all blocks that can come, + while let Some(timeblock) = self.timebuf.pop(self.nfft(), self.overlap_keep) { + // Compute cross-power spectra for current time block + self.update_singleblock(&timeblock); + + result.push(self.cur_est.clone()); + } + + result } } @@ -352,9 +355,9 @@ mod test { use ndarray_rand::RandomExt; use super::CrossPowerSpecra; - use crate::{config::*, ps::ApsResult}; + use crate::config::*; - use super::{ApsMode, AvPowerSpectra, Overlap, WindowType, CPS}; + use super::{ApsMode, AvPowerSpectra, CPSResult, Overlap, WindowType}; #[test] fn test_overlap_keep() { @@ -399,21 +402,21 @@ mod test { let timedata_some = Dmat::random((nfft, 1), distr); let timedata_zeros = Dmat::zeros((nfft, 1)); - let mut first_result = CPS::zeros((0, 0, 0)); + let mut first_result = CPSResult::zeros((0, 0, 0)); - if let ApsResult::OnlyLastResult(v) = aps.compute_last(timedata_some.view()) { + if let Some(v) = aps.compute_last(timedata_some.view()) { first_result = v.clone(); // println!("{:?}", first_result.ap(0)[0]); } else { assert!(false, "Should return one value"); } let overlap_keep = AvPowerSpectra::get_overlap_keep(nfft, Overlap::NoOverlap).unwrap(); - if let ApsResult::OnlyLastResult(_) = aps.compute_last(&timedata_zeros) { + if let Some(_) = aps.compute_last(&timedata_zeros) { // Do nothing with it. } else { assert!(false, "Should return one value"); } - if let ApsResult::OnlyLastResult(v) = aps.compute_last(&timedata_zeros) { + if let Some(v) = aps.compute_last(&timedata_zeros) { let alpha = Flt::exp(-((nfft - overlap_keep) as Flt) / (fs * tau)); // let alpha: Flt = 1.; for i in 0..nfft / 2 + 1 { @@ -436,7 +439,7 @@ mod test { .unwrap(); let mut aps = AvPowerSpectra::build(nfft, None, None, None).unwrap(); - if let ApsResult::OnlyLastResult(v) = aps.compute_last(&timedata) { + if let Some(v) = aps.compute_last(&timedata) { let tf = v.tf(0, 1, None); assert_eq!((&tf - 2.0 * Cflt::ONE).sum().abs(), 0.0); } else { @@ -458,7 +461,7 @@ mod test { .unwrap(); let mut aps = AvPowerSpectra::build(nfft, None, None, None).unwrap(); - if let ApsResult::OnlyLastResult(v) = aps.compute_last(&timedata) { + if let Some(v) = aps.compute_last(&timedata) { let tf = v.tf(0, 1, Some(2)); assert_eq!((&tf - 2.0 * Cflt::ONE).sum().abs(), 0.0); } else { @@ -480,7 +483,7 @@ mod test { None, ] { let mut aps = AvPowerSpectra::build(nfft, wt, None, None).unwrap(); - if let ApsResult::OnlyLastResult(v) = aps.compute_last(&timedata) { + if let Some(v) = aps.compute_last(&timedata) { let ap = v.ap(0); assert_abs_diff_eq!((&ap).sum().abs(), timedata_mean_square, epsilon = 1e-2); } else { diff --git a/src/ps/mod.rs b/src/ps/mod.rs index 32ffdfd..7cefbff 100644 --- a/src/ps/mod.rs +++ b/src/ps/mod.rs @@ -12,6 +12,6 @@ mod window; use crate::config::*; -pub use aps::{ApsMode, ApsResult, AvPowerSpectra, Overlap}; -pub use ps::{CrossPowerSpecra, PowerSpectra, CPS}; +pub use aps::{ApsMode, AvPowerSpectra, Overlap}; +pub use ps::{CrossPowerSpecra, PowerSpectra, CPSResult}; pub use window::{Window, WindowType}; diff --git a/src/ps/ps.rs b/src/ps/ps.rs index 3843dd3..46bbde5 100644 --- a/src/ps/ps.rs +++ b/src/ps/ps.rs @@ -16,9 +16,9 @@ use realfft::{RealFftPlanner, RealToComplex}; /// Cross power spectra, which is a 3D array, with the following properties: /// /// - The first index is the frequency index, starting at DC, ending at nfft/2. -/// - The second, and third index result in [i,j] = C_ij = p_i * conj(p_j) +/// - The second, and third index result in `[i,j]` = C_ij = p_i * conj(p_j) /// -pub type CPS = Array3; +pub type CPSResult = Array3; /// Extra typical methods that are of use for 3D-arrays of complex numbers, that /// are typically implemented as cross-power spectra. @@ -47,7 +47,7 @@ pub trait CrossPowerSpecra { } -impl CrossPowerSpecra for CPS { +impl CrossPowerSpecra for CPSResult { fn ap(&self, ch: usize) -> Array1 { // Slice out one value for all frequencies, map to only real part, and // return. @@ -192,7 +192,7 @@ impl PowerSpectra { /// (nfft/2+1,timedata.ncols(), timedata.ncols()). Its content is: /// [freq_index, chi, chj] = crosspower: chi*conj(chj) /// - pub fn compute<'a, T>(&mut self, tdata: T) -> CPS + pub fn compute<'a, T>(&mut self, tdata: T) -> CPSResult where T: AsArray<'a, Flt, Ix2>, {