ASCEE/lasprs
5
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

AvPowerSpectra not yet working, but going in right direction

Browse Source
This commit is contained in:
Anne de Jong 2024-07-06 19:38:49 +02:00
parent 0a847318f3
commit 565829f1fa
2 changed files with 334 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

200
src/ps/aps.rs
View File

@ -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<CrossPowerSpecra>),
/// 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<Flt>,
/// 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<Overlap>) -> Result<usize> {
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<usize> {
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<AvPowerSpectra> {
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>,
windowtype: Option<WindowType>,
overlap: Option<Overlap>,
fs_tau: Option<Flt>,
mode: Option<ApsMode>,
) -> Result<AvPowerSpectra> {
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);
}
}

166
src/ps/timebuffer.rs
View File

@ -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;
}
/// 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<VecDeque<Flt>>,
}
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<Dmat> {
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);
}
}