diff --git a/src/ps/aps.rs b/src/ps/aps.rs index 72dd5e2..b7a0c9c 100644 --- a/src/ps/aps.rs +++ b/src/ps/aps.rs @@ -7,7 +7,10 @@ use anyhow::{bail, Result}; /// Can be provided as a percentage of the block size, or as a number of /// samples. pub enum Overlap { + /// Overlap specified as a percentage of the total FFT length. Value should + /// be 0<=pct<100. Percentage(Flt), + /// Number of samples to overlap Number(usize), } impl Default for Overlap { @@ -16,42 +19,136 @@ impl Default for Overlap { } } +/// Result from [AvPowerspectra.compute()] +pub enum ApsResult { + /// 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), + + /// Give only last result back, the most recent value. + OnlyLastResult(CrossPowerSpecra), + + /// No new data available. Nothing given back. + None, +} + +/// The 'mode' used in computing averaged power spectra. When providing data in +/// blocks to the [AvPowerSpectra] the resulting 'current estimate' responds +/// differently, depending on the model. +pub enum ApsMode { + /// Averaged over all data provided. New averages can be created by calling + /// `AvPowerSpectra::reset()` + AllAveraging, + /// In this mode, the `AvPowerSpectra` works a bit like a sound level meter, + /// 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] + fs: Flt, + /// Time weighting constant, its inverse is its approximate -3 dB point + /// for forgetting old history. + tau: Flt, + }, + /// Spectrogram mode. Only returns the latest estimates. + Spectrogram, +} +impl Default for ApsMode { + fn default() -> Self { + ApsMode::AllAveraging + } +} + /// Averaged power spectra computing engine -struct AvPowerSpectra { +/// Used to compute power spectra estimations on +/// long datasets, where nfft << length of data. This way, the variance of a +/// single periodogram is suppressed with increasing number of averages. +/// +/// For more information, see the book on numerical recipes. +/// +pub struct AvPowerSpectra { + // Power spectra estimator for single block ps: PowerSpectra, - overlap_skip: usize, - // Time constant for time-weighted power spectra. If None, it averages out - // over all data. - fs_tau: Option, - /// The number of frames already used in the average + + mode: ApsMode, + + // The number of samples to keep in the time buffer when overlapping time + // blocks + overlap_keep: usize, + + /// The number of blocks of length [self.nfft()] already used in the average N: usize, + /// Storage for sample data. + timebuf: TimeBuffer, + // Current estimation of the power spectra - cur_est: Cmat + cur_est: Cmat, } impl AvPowerSpectra { - fn get_overlap_skip(nfft: usize, overlap: Option) -> Result { - let overlap = overlap.unwrap_or_default(); + /// The FFT Length of estimating (cross)power spectra + pub fn nfft(&self) -> usize { + self.ps.nfft() + } + /// Returns the amount of samples to `keep` in the time buffer when + /// overlapping time segments using [TimeBuffer]. + fn get_overlap_keep(nfft: usize, overlap: Overlap) -> Result { + let overlap_keep = match overlap { + Overlap::Number(i) if i >= nfft => { + bail!("Invalid overlap number of samples. Should be < nfft, which is {nfft}.") + } + // Keep 1 sample, if overlap is 1 sample etc. + Overlap::Number(i) if i < nfft => i, - let overlap_skip = match overlap { - Overlap::Number(i) if i >= nfft => bail!("Invalid overlap number of samples"), - Overlap::Number(i) if i < nfft => nfft - i, - Overlap::Percentage(p) if p >= 100. => { + // If overlap percentage is >= 100, or < 0.0 its an error + Overlap::Percentage(p) if p >= 100. || p < 0.0 => { bail!("Invalid overlap percentage. Should be >= 0. And < 100.") } - Overlap::Percentage(p) => nfft - ((p * nfft as Flt) / 100.) as usize, + // If overlap percentage is 0, this gives + Overlap::Percentage(p) => ((p * nfft as Flt) / 100.) as usize, _ => unreachable!(), }; - if overlap_skip == 0 || overlap_skip > nfft { + if overlap_keep >= nfft { bail!("Computed overlap results in invalid number of overlap samples. Please make sure the FFT length is large enough, when high overlap percentages are required."); } - Ok(overlap_skip) + Ok(overlap_keep) } + + /// Create new averaged power spectra estimator for weighing over the full + /// amount of data supplied (no exponential spectra weighting) using + /// sensible defaults (Hann window, 50% overlap). + /// + /// # Args + /// + /// * `nfft` - FFT Length + pub fn new_simple(nfft: usize) -> Result { + return AvPowerSpectra::new(nfft, None, None, None); + } + + /// Create power spectra estimator which weighs either all data + /// (`fs_tau=None`), or uses exponential weighting. Other parameters can be + /// provided, like the overlap (Some(Overlap)). Otherwise: all parameters + /// have sensible defaults. + /// + /// # Exponential weighting + /// + /// This can be used to `follow` a power spectrum as a function of time. + /// New data is weighted more heavily. Note that this way of looking at + /// spectra is not 'exact', and therefore should not be used for + /// spectrograms. + /// + /// # Args + /// + /// - `nfft` - The discrete Fourier Transform length used in the estimation. + /// - `windowtype` - Window Type. The window type to use Hann, etc. + /// - `overlap` - Amount of overlap in + /// - `mode` - The mode in which the + /// pub fn new( nfft: usize, - wt: Option, + windowtype: Option, overlap: Option, - fs_tau: Option, + mode: Option, ) -> Result { if nfft % 2 != 0 { bail!("NFFT should be even") @@ -59,22 +156,73 @@ impl AvPowerSpectra { if nfft == 0 { bail!("Invalid NFFT") } + let windowtype = windowtype.unwrap_or_default(); + let window = Window::new(windowtype, nfft); + let mode = mode.unwrap_or_default(); + let overlap = overlap.unwrap_or_default(); - if let Some(x) = fs_tau { - if x <= 0.0 { - bail!("Invalid time weighting constant [s]. Should be > 0 if given.") + // Perform some checks on ApsMode + if let ApsMode::ExponentialWeighting { fs, tau } = mode { + if fs <= 0.0 { + bail!("Invalid sampling frequency given as parameter"); + } + if tau <= 0.0 { + bail!("Invalid time weighting constant [s]. Should be > 0 if given."); } } - let window = Window::new(wt.unwrap_or_default(), nfft); let ps = PowerSpectra::newFromWindow(window); - let overlap_skip = Self::get_overlap_skip(nfft, overlap)?; + let overlap_keep = Self::get_overlap_keep(nfft, overlap)?; + Ok(AvPowerSpectra { ps, - overlap_skip, - fs_tau, + overlap_keep, + mode, N: 0, - cur_est: Cmat::default((0,0)) + cur_est: Cmat::default((0, 0)), + timebuf: TimeBuffer::new(), }) } + /// Computes average (cross)power spectra, and returns only the most recent + /// estimate, if any can be given back. Only gives back a result when enough + /// data is available. + /// + /// # Args + /// + /// * `timedata``: New available time data. Number of columns is number of + /// channels, number of rows is number of frames (samples per channel). + /// * `giveInterMediateResults` + /// + /// # Panics + /// + /// If timedata.ncols() does not match number of columns in already present + /// data. + fn compute<'a, T>(&mut self, timedata: T, giveInterMediateResults: bool) -> ApsResult + where + T: AsArray<'a, Flt, Ix2>, + { + // Push new data in the time buffer. + self.timebuf.push(timedata); + + // 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 + let Cpsnew = self.ps.compute(&timeblock); + // Ok, mode-determined + match self.mode { + ApsMode::AllAveraging => {} + ApsMode::ExponentialWeighting { fs, tau } => {} + ApsMode::Spectrogram => {} + } + // if self.cur_est.len() > 0 + } + ApsResult::None + } + /// See [AvPowerSpectra](compute()) + fn compute_last<'a, T>(&mut self, timedata: T) -> ApsResult + where + T: AsArray<'a, Flt, Ix2>, + { + return self.compute(timedata, false); + } } diff --git a/src/ps/timebuffer.rs b/src/ps/timebuffer.rs index 90e4f1b..d27b104 100644 --- a/src/ps/timebuffer.rs +++ b/src/ps/timebuffer.rs @@ -1,8 +1,162 @@ +//! A time buffer is used as intermediate storage for samples, to spare up +//! enough values to do 'something' with. Typical application: getting enough +//! samples to perform an FFT. The use case is computing average power spectra. use super::*; -/// A time buffer is used as intermediate storage for samples, to spare up -/// enough values to do 'something' with. Typical application: getting enough -/// samples to perform an FFT. The use case is computing average power spectra. -/// -struct TimeBuffer { +use crate::config::*; +use std::collections::VecDeque; -} \ No newline at end of file +/// TimeBuffer, storage to add blocks of data in a ring buffer, that can be +/// extracted by blocks of other size. Also, we can keep samples in a buffer to +/// create, for example, overlapping windows of time data. +pub struct TimeBuffer { + data: Vec>, +} +impl TimeBuffer { + /// Allocate new empty TimeBuffer. + pub fn new() -> TimeBuffer { + TimeBuffer { data: Vec::new() } + } + /// Reset the buffer. Clears all data in it. + pub fn reset(&mut self) { + *self = Self::new(); + } + /// Push new data to the buffer. Keep + /// + /// # Panics + /// + /// When the number of columns does not match with the first time this + /// method is called. + pub fn push<'a, T>(&mut self, new_data: T) + where + T: AsArray<'a, Flt, Ix2>, + { + let new_data = new_data.into(); + let nch = new_data.shape()[1]; + if self.data.len() == 0 { + self.data = vec![VecDeque::new(); nch]; + } + if self.data.len() != nch { + panic!("BUG: Number of channels does not match with time buffer"); + } + for (new_dat, buf) in new_data.columns().into_iter().zip(self.data.iter_mut()) { + let slice = new_dat.as_slice().unwrap(); + buf.extend(slice); + } + } + /// Return the number of samples that is currently stored + pub fn nsamples(&self) -> usize { + if let Some(q) = self.data.get(0) { + return q.len(); + }; + 0 + } + /// Pop samples from the queue, only returns Some(Dmat) when there are + /// enough samples to return. Never returns less samples than + /// `nsamples_requested`. + pub fn pop(&mut self, nsamples_requested: usize, nsamples_keep: usize) -> Option { + if self.data.len() == 0 { + return None; + } + if nsamples_keep > nsamples_requested { + panic!("BUG: Cannot keep more samples than requested to return"); + } + + debug_assert!(self.data.len() > 0); + let c1 = unsafe { self.data.get_unchecked(0) }; + + let nsamples_available = c1.len(); + debug_assert!(nsamples_available == self.nsamples()); + if nsamples_available < nsamples_requested { + println!("Less available than requested"); + return None; + } + // Number of channels + let nchannels = self.data.len(); + + // Create output array, .f() means fortran contiguous array + let mut res = Dmat::zeros((nsamples_requested, nchannels).f()); + { + for (mut col, dat) in res.columns_mut().into_iter().zip(&mut self.data) { + // This expect should never happen, as above it is called with .f(). + let col_slice = col + .as_slice_mut() + .expect("Data is not contiguous on the sample axis!"); + + // Split data of current channel in two slices of the underlying data + let (dat_slice1, dat_slice2) = dat.as_slices(); + + if dat_slice1.len() >= nsamples_requested { + // Only copy from the first slice + col_slice.copy_from_slice(&dat_slice1[..nsamples_requested]); + } else { + let slice1len = dat_slice1.len(); + // Copy from first slice + col_slice[..slice1len].copy_from_slice(&dat_slice1); + // Copy rest from second slice + col_slice[slice1len..nsamples_requested] + .copy_from_slice(&dat_slice2[..nsamples_requested - slice1len]); + } + + // From documentation: + // Removes the specified range from the deque in bulk, returning all removed + // elements as an iterator. If the iterator is dropped before being fully + // consumed, it drops the remaining removed elements. + dat.drain(0..nsamples_requested - nsamples_keep); + } + } + let c1 = unsafe { self.data.get_unchecked(0) }; + let nsamples_available = c1.len(); + debug_assert!(self + .data + .iter() + .map(|ch| ch.len()) // Iterator over all channels, transform to the amount of samples in each channel. + .all(|v| v == nsamples_available)); + + Some(res) + } +} + +#[cfg(test)] +mod test { + use crate::{Dcol, Dmat}; + use ndarray::s; + + use super::TimeBuffer; + + #[test] + fn test_timebuffer1() { + let mut t1 = Dmat::zeros((10, 1)); + t1[[9, 0]] = 1.0; + let mut tb = TimeBuffer::new(); + + tb.push(&t1); + let out = tb.pop(10, 2).unwrap(); + assert_eq!(out.len(), 10); + let out = tb.pop(3, 0); + assert!(out.is_none()); + let out = tb.pop(2, 0).unwrap(); + // println!("{:?}", out); + assert_eq!(out.len(), 2); + assert_eq!(tb.nsamples(), 0); + assert_eq!(out[[1, 0]], 1.0); + + tb.reset(); + assert_eq!(tb.nsamples(), 0); + } + + #[test] + fn test_timebuffer2() { + let mut tb = TimeBuffer::new(); + let tlin = Dcol::linspace(0., 100., 101); + // println!("{:?}", tlin); + let mut t2 = Dmat::zeros((101, 0)); + t2.push_column(tlin.view()).unwrap(); + t2.push_column(tlin.view()).unwrap(); + tb.push(&t2); + assert_eq!(tb.nsamples(), 101); + + let tres = tb.pop(50, 49).unwrap(); + assert_eq!(tres.shape(), [50, 2]); + assert_eq!(t2.slice(s![..50, ..2]), tres); + } +}