Removed unnecessary complex API with APSResult and switched to two methods: compute_last and compute_all. Improved docstrings and fixes in docstrings.

This commit is contained in:
Anne de Jong 2024-07-11 21:02:21 +02:00
parent 9db49ff8c7
commit 66ab149159
7 changed files with 66 additions and 63 deletions

View File

@ -6,7 +6,7 @@ and processing of (multi) sensor data in real time on a PC and output results.
## Documentation ## 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 ## Python bindings

View File

@ -15,11 +15,11 @@ struct Cli {
/// File name to write recording to /// File name to write recording to
filename: String, 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.)] #[arg(short, long = "duration", default_value_t = 0.)]
duration_s: Flt, 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. /// start delay will be used.
#[arg(short, long = "startdelay", default_value_t = 0.)] #[arg(short, long = "startdelay", default_value_t = 0.)]
start_delay_s: Flt, start_delay_s: Flt,

View File

@ -173,8 +173,8 @@ impl Biquad {
/// First order low pass filter (one pole in the real axis). No pre-warping /// First order low pass filter (one pole in the real axis). No pre-warping
/// correction done. /// correction done.
/// ///
/// * `fs` - Sampling frequency [Hz] /// * `fs` - Sampling frequency \[Hz\]
/// * `fc` - Cut-off frequency (-3 dB point) [Hz] /// * `fc` - Cut-off frequency (-3 dB point) \[Hz\]
pub fn firstOrderLowPass(fs: Flt, fc: Flt) -> Result<Biquad> { pub fn firstOrderLowPass(fs: Flt, fc: Flt) -> Result<Biquad> {
if fc <= 0. { if fc <= 0. {
bail!("Cuton frequency, given: should be > 0") bail!("Cuton frequency, given: should be > 0")

View File

@ -42,7 +42,7 @@ where
/// ///
/// # Args /// # Args
/// ///
/// * `freq` - The frequency in [Hz] /// * `freq` - The frequency in \[Hz\]
/// ///
/// # Returns /// # Returns
/// ///

View File

@ -22,15 +22,15 @@ impl Default for Overlap {
} }
} }
/// Result from [AvPowerspectra.compute()] /// Result from [compute method](AvPowerSpectra::compute).
pub enum ApsResult<'a> { enum ApsResult<'a> {
/// Returns all intermediate results. Useful when evolution over time should /// Returns all intermediate results. Useful when evolution over time should
/// be visible. I.e. in case of spectrograms, but also in case the /// be visible. I.e. in case of spectrograms, but also in case the
/// converging process to the average should be made visible. /// converging process to the average should be made visible.
AllIntermediateResults(Vec<CPS>), AllIntermediateResults(Vec<CPSResult>),
/// Give only last result back, the most recent value. /// Give only last result back, the most recent value.
OnlyLastResult(&'a CPS), OnlyLastResult(&'a CPSResult),
/// No new data available. Nothing given back. /// No new data available. Nothing given back.
None, None,
@ -47,7 +47,7 @@ pub enum ApsMode {
/// where new data is weighted with old data, and old data exponentially /// where new data is weighted with old data, and old data exponentially
/// backs off. This mode only makes sense when `tau >> nfft/fs` /// backs off. This mode only makes sense when `tau >> nfft/fs`
ExponentialWeighting { ExponentialWeighting {
/// Sampling frequency in [Hz], used for computing IIR filter coefficient. /// Sampling frequency in `[Hz]`, used for computing IIR filter coefficient.
fs: Flt, fs: Flt,
/// Time weighting constant, follows convention of Sound Level Meters. /// Time weighting constant, follows convention of Sound Level Meters.
/// Means the data is approximately low-pass filtered with a cut-off /// Means the data is approximately low-pass filtered with a cut-off
@ -88,7 +88,7 @@ pub struct AvPowerSpectra {
timebuf: TimeBuffer, timebuf: TimeBuffer,
// Current estimation of the power spectra // Current estimation of the power spectra
cur_est: CPS, cur_est: CPSResult,
} }
impl AvPowerSpectra { impl AvPowerSpectra {
/// The FFT Length of estimating (cross)power spectra /// The FFT Length of estimating (cross)power spectra
@ -125,7 +125,7 @@ impl AvPowerSpectra {
pub fn reset(&mut self) { pub fn reset(&mut self) {
self.N = 0; self.N = 0;
self.timebuf.reset(); 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 /// Create new averaged power spectra estimator for weighing over the full
/// amount of data supplied (no exponential spectra weighting) using /// amount of data supplied (no exponential spectra weighting) using
@ -140,7 +140,7 @@ impl AvPowerSpectra {
/// # Panics /// # Panics
/// ///
/// - When nfft is not even, or 0. /// - When nfft is not even, or 0.
/// ///
pub fn new_simple(nfft: usize) -> AvPowerSpectra { pub fn new_simple(nfft: usize) -> AvPowerSpectra {
AvPowerSpectra::build(nfft, None, None, None).unwrap() AvPowerSpectra::build(nfft, None, None, None).unwrap()
} }
@ -199,7 +199,7 @@ impl AvPowerSpectra {
overlap_keep, overlap_keep,
mode, mode,
N: 0, N: 0,
cur_est: CPS::default((0, 0, 0)), cur_est: CPSResult::default((0, 0, 0)),
timebuf: TimeBuffer::new(), timebuf: TimeBuffer::new(),
}) })
} }
@ -215,7 +215,7 @@ impl AvPowerSpectra {
// Initialize to zero // Initialize to zero
if self.cur_est.len() == 0 { if self.cur_est.len() == 0 {
assert_eq!(self.N, 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 // Update the number of blocks processed
@ -285,27 +285,18 @@ impl AvPowerSpectra {
/// ///
/// * `timedata``: New available time data. Number of columns is number of /// * `timedata``: New available time data. Number of columns is number of
/// channels, number of rows is number of frames (samples per channel). /// 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 /// # Panics
/// ///
/// If timedata.ncols() does not match number of columns in already present /// If timedata.ncols() does not match number of columns in already present
/// data. /// data.
pub fn compute<'a, 'b, T>( pub fn compute_last<'a, 'b, T>(&'a mut self, timedata: T) -> Option<&'a CPSResult>
&'a mut self,
timedata: T,
giveInterMediateResults: bool,
) -> ApsResult<'a>
where where
T: AsArray<'b, Flt, Ix2>, T: AsArray<'b, Flt, Ix2>,
{ {
// Push new data in the time buffer. // Push new data in the time buffer.
self.timebuf.push(timedata); self.timebuf.push(timedata);
// Storage for the result
let mut result = ApsResult::None;
// Flag to indicate that we have obtained one result for sure. // Flag to indicate that we have obtained one result for sure.
let mut computed_single = false; let mut computed_single = false;
@ -314,34 +305,46 @@ impl AvPowerSpectra {
// Compute cross-power spectra for current time block // Compute cross-power spectra for current time block
self.update_singleblock(&timeblock); self.update_singleblock(&timeblock);
// We ha
computed_single = true; 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 {
if computed_single && !giveInterMediateResults { return Some(&self.cur_est);
// 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);
} }
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<CPSResult>
where 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 ndarray_rand::RandomExt;
use super::CrossPowerSpecra; 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] #[test]
fn test_overlap_keep() { fn test_overlap_keep() {
@ -399,21 +402,21 @@ mod test {
let timedata_some = Dmat::random((nfft, 1), distr); let timedata_some = Dmat::random((nfft, 1), distr);
let timedata_zeros = Dmat::zeros((nfft, 1)); 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(); first_result = v.clone();
// println!("{:?}", first_result.ap(0)[0]); // println!("{:?}", first_result.ap(0)[0]);
} else { } else {
assert!(false, "Should return one value"); assert!(false, "Should return one value");
} }
let overlap_keep = AvPowerSpectra::get_overlap_keep(nfft, Overlap::NoOverlap).unwrap(); 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. // Do nothing with it.
} else { } else {
assert!(false, "Should return one value"); 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::exp(-((nfft - overlap_keep) as Flt) / (fs * tau));
// let alpha: Flt = 1.; // let alpha: Flt = 1.;
for i in 0..nfft / 2 + 1 { for i in 0..nfft / 2 + 1 {
@ -436,7 +439,7 @@ mod test {
.unwrap(); .unwrap();
let mut aps = AvPowerSpectra::build(nfft, None, None, None).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); let tf = v.tf(0, 1, None);
assert_eq!((&tf - 2.0 * Cflt::ONE).sum().abs(), 0.0); assert_eq!((&tf - 2.0 * Cflt::ONE).sum().abs(), 0.0);
} else { } else {
@ -458,7 +461,7 @@ mod test {
.unwrap(); .unwrap();
let mut aps = AvPowerSpectra::build(nfft, None, None, None).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)); let tf = v.tf(0, 1, Some(2));
assert_eq!((&tf - 2.0 * Cflt::ONE).sum().abs(), 0.0); assert_eq!((&tf - 2.0 * Cflt::ONE).sum().abs(), 0.0);
} else { } else {
@ -480,7 +483,7 @@ mod test {
None, None,
] { ] {
let mut aps = AvPowerSpectra::build(nfft, wt, None, None).unwrap(); 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); let ap = v.ap(0);
assert_abs_diff_eq!((&ap).sum().abs(), timedata_mean_square, epsilon = 1e-2); assert_abs_diff_eq!((&ap).sum().abs(), timedata_mean_square, epsilon = 1e-2);
} else { } else {

View File

@ -12,6 +12,6 @@ mod window;
use crate::config::*; use crate::config::*;
pub use aps::{ApsMode, ApsResult, AvPowerSpectra, Overlap}; pub use aps::{ApsMode, AvPowerSpectra, Overlap};
pub use ps::{CrossPowerSpecra, PowerSpectra, CPS}; pub use ps::{CrossPowerSpecra, PowerSpectra, CPSResult};
pub use window::{Window, WindowType}; pub use window::{Window, WindowType};

View File

@ -16,9 +16,9 @@ use realfft::{RealFftPlanner, RealToComplex};
/// Cross power spectra, which is a 3D array, with the following properties: /// 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 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<Cflt>; pub type CPSResult = Array3<Cflt>;
/// Extra typical methods that are of use for 3D-arrays of complex numbers, that /// Extra typical methods that are of use for 3D-arrays of complex numbers, that
/// are typically implemented as cross-power spectra. /// 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<Flt> { fn ap(&self, ch: usize) -> Array1<Flt> {
// Slice out one value for all frequencies, map to only real part, and // Slice out one value for all frequencies, map to only real part, and
// return. // return.
@ -192,7 +192,7 @@ impl PowerSpectra {
/// (nfft/2+1,timedata.ncols(), timedata.ncols()). Its content is: /// (nfft/2+1,timedata.ncols(), timedata.ncols()). Its content is:
/// [freq_index, chi, chj] = crosspower: chi*conj(chj) /// [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 where
T: AsArray<'a, Flt, Ix2>, T: AsArray<'a, Flt, Ix2>,
{ {