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

Moved apsmode and apssettings to their own files

Browse Source
This commit is contained in:
Anne de Jong 2024-09-28 14:22:35 +02:00
parent 08ecdf6dc4
commit 826266b8ee
11 changed files with 254 additions and 185 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/config.rs
View File

@ -26,7 +26,8 @@ cfg_if::cfg_if! {
cfg_if::cfg_if! {
if #[cfg(feature = "python-bindings")] {
pub use numpy::{IntoPyArray,PyArray, PyArray1, PyArrayDyn, PyArrayLike1, PyReadonlyArrayDyn};
pub use numpy::{IntoPyArray,PyArray, PyArray1, PyArray2, PyArray3, PyArrayDyn, PyArrayLike1,
PyArrayLike2,PyArrayLike3,PyReadonlyArrayDyn, convert::ToPyArray};
pub use pyo3::prelude::*;
pub use pyo3::exceptions::PyValueError;
pub use pyo3::{pymodule, types::PyModule, PyResult};

3
src/lib.rs
View File

@ -64,6 +64,9 @@ fn lasprs(m: &Bound<'_, PyModule>) -> PyResult<()> {
m.add_class::<slm::SLM>()?;
m.add_class::<ps::WindowType>()?;
m.add_class::<ps::Overlap>()?;
m.add_class::<ps::ApsMode>()?;
m.add_class::<ps::ApsSettings>()?;
m.add_class::<ps::AvPowerSpectra>()?;
Ok(())
}

209
src/ps/aps.rs
View File

@ -1,198 +1,25 @@
use super::timebuffer::TimeBuffer;
use super::CrossPowerSpecra;
use super::*;
use super::{timebuffer::TimeBuffer, CrossPowerSpecra};
use crate::{config::*, TransferFunction, ZPKModel};
use anyhow::{bail, Error, Result};
use derive_builder::Builder;
use freqweighting::FreqWeighting;
/// All settings used for computing averaged power spectra using Welch' method.
#[derive(Builder, Clone)]
#[builder(build_fn(validate = "Self::validate", error = "Error"))]
pub struct ApsSettings {
/// Mode of computation, see [ApsMode].
#[builder(default)]
mode: ApsMode,
/// Overlap in time segments. See [Overlap].
#[builder(default)]
overlap: Overlap,
/// Window applied to time segments. See [WindowType].
#[builder(default)]
windowType: WindowType,
/// Kind of freqency weighting. Defaults to Z
#[builder(default)]
freqWeightingType: FreqWeighting,
/// FFT Length
nfft: usize,
/// Sampling frequency
fs: Flt,
}
impl ApsSettingsBuilder {
fn validate(&self) -> Result<()> {
if !self.fs.is_some() {
bail!("Sampling frequency not given");
}
let fs = self.fs.unwrap();
if !fs.is_normal() {
bail!("Sampling frequency not a normal number")
}
if fs <= 0.0 {
bail!("Invalid sampling frequency given as parameter");
}
if self.nfft.is_none() {
bail!("nfft not specified")
};
let nfft = self.nfft.unwrap();
if nfft % 2 != 0 {
bail!("NFFT should be even")
}
if nfft == 0 {
bail!("Invalid NFFT, should be > 0.")
}
// Perform some checks on ApsMode
if let Some(ApsMode::ExponentialWeighting { tau }) = self.mode {
if tau <= 0.0 {
bail!("Invalid time weighting constant [s]. Should be > 0 if given.");
}
}
Ok(())
}
}
impl ApsSettings {
/// Returns nfft
pub fn nfft(&self) -> usize {
self.nfft
}
fn get_overlap_keep(&self) -> usize {
self.validate_get_overlap_keep().unwrap()
}
/// Returns the amount of samples to `keep` in the time buffer when
/// overlapping time segments using [TimeBuffer].
fn validate_get_overlap_keep(&self) -> Result<usize> {
let nfft = self.nfft;
let overlap_keep = match self.overlap {
Overlap::Number { N } if N >= nfft => {
bail!("Invalid overlap number of samples. Should be < nfft, which is {nfft}.")
}
// Keep 1 sample, if overlap is 1 sample etc.
Overlap::Number { N } if N < nfft => N,
// If overlap percentage is >= 100, or < 0.0 its an error
Overlap::Percentage { pct } if !(0.0..100.).contains(&pct) => {
bail!("Invalid overlap percentage. Should be >= 0. And < 100.")
}
// If overlap percentage is 0, this gives
Overlap::Percentage { pct } => ((pct * nfft as Flt) / 100.) as usize,
Overlap::NoOverlap {} => 0,
_ => unreachable!(),
};
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_keep)
}
/// Return a reasonable acoustic default with a frequency resolution around
/// ~ 10 Hz, where nfft is still an integer power of 2.
///
/// # Errors
///
/// If `fs` is something odd, i.e. < 1 kHz, or higher than 1 MHz.
///
pub fn reasonableAcousticDefault(fs: Flt, mode: ApsMode) -> Result<ApsSettings> {
if fs < 1e3 || fs > 1e6 {
bail!("Sampling frequency for reasonable acoustic data is >= 1 kHz and <= 1 MHz.");
}
let fs_div_10_rounded = (fs / 10.) as u32;
// 2^30 is about 1 million. We search for a two-power of an nfft that is
// the closest to fs/10. The frequency resolution is about fs/nfft.
let nfft = (0..30).map(|i| 2u32.pow(i) - fs_div_10_rounded).fold(
// Start wth a value that is always too large
fs as u32 * 10,
|cur, new| cur.min(new),
) as usize;
Ok(ApsSettings {
mode,
fs,
nfft,
windowType: WindowType::default(),
overlap: Overlap::default(),
freqWeightingType: FreqWeighting::default(),
})
}
/// Return sampling frequency
pub fn fs(&self) -> Flt {
self.fs
}
/// Return Nyquist frequency
pub fn fnyq(&self) -> Flt {
self.fs() / 2.
}
/// Returns a single-sided frequency array corresponding to points in Power
/// spectra computation.
pub fn getFreq(&self) -> Array1<Flt> {
let df = self.fs / self.nfft as Flt;
let K = self.nfft / 2 + 1;
Array1::linspace(0., (K - 1) as Flt * df, K)
}
}
/// 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.
#[derive(Copy, Clone, PartialEq)]
#[cfg_attr(feature = "python-bindings", pyclass)]
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 {
/// Time weighting constant, follows convention of Sound Level Meters.
/// Means the data is approximately low-pass filtered with a cut-off
/// frequency f_c of s/tau ≅ 1 → f_c = (2 * pi * tau)^-1.
tau: Flt,
},
/// Spectrogram mode. Only returns the latest estimate(s).
Spectrogram {},
}
impl Default for ApsMode {
fn default() -> Self {
ApsMode::AllAveraging {}
}
}
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)]
impl ApsMode {
#[inline]
fn __eq__(&self, other: &Self) -> bool {
self == other
}
}
/// Averaged power spectra computing engine
/// 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.
///
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pyclass)]
#[derive(Debug)]
pub struct AvPowerSpectra {
// Power spectra estimator for single block
ps: PowerSpectra,
// Settings for computing power spectra, see [ApsSettings]
settings: ApsSettings,
// The number of samples to keep in the time buffer when overlapping time
@ -454,6 +281,28 @@ impl AvPowerSpectra {
}
}
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)]
impl AvPowerSpectra {
#[new]
fn new_py(s: ApsSettings) -> AvPowerSpectra {
AvPowerSpectra::new(s)
}
#[pyo3(name = "compute")]
fn compute_py<'py>(
&mut self,
py: Python<'py>,
dat: PyArrayLike2<Flt>,
) -> Bound<'py, PyArray3<Cflt>> {
let dat = dat.as_array();
if let Some(res) = self.compute_last(dat) {
let res = res.clone();
return res.to_pyarray_bound(py);
}
panic!("No data!");
}
}
#[cfg(test)]
mod test {
use approx::assert_abs_diff_eq;

35
src/ps/apsmode.rs Normal file
View File

@ -0,0 +1,35 @@
use crate::config::*;
/// 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.
#[derive(Copy, Clone, PartialEq, Debug)]
#[cfg_attr(feature = "python-bindings", pyclass)]
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 {
/// Time weighting constant, follows convention of Sound Level Meters.
/// Means the data is approximately low-pass filtered with a cut-off
/// frequency f_c of s/tau ≅ 1 → f_c = (2 * pi * tau)^-1.
tau: Flt,
},
/// Spectrogram mode. Only returns the latest estimate(s).
Spectrogram {},
}
impl Default for ApsMode {
fn default() -> Self {
ApsMode::AllAveraging {}
}
}
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)]
impl ApsMode {
#[inline]
fn __eq__(&self, other: &Self) -> bool {
self == other
}
}

168
src/ps/apssettings.rs Normal file
View File

@ -0,0 +1,168 @@
use super::*;
use crate::config::*;
use anyhow::{bail, Error, Result};
use derive_builder::Builder;
/// All settings used for computing averaged power spectra using Welch' method.
#[derive(Builder, Clone, Debug)]
#[cfg_attr(feature = "python-bindings", pyclass)]
#[builder(build_fn(validate = "Self::validate", error = "Error"))]
pub struct ApsSettings {
/// Mode of computation, see [ApsMode].
#[builder(default)]
pub mode: ApsMode,
/// Overlap in time segments. See [Overlap].
#[builder(default)]
pub overlap: Overlap,
/// Window applied to time segments. See [WindowType].
#[builder(default)]
pub windowType: WindowType,
/// Kind of freqency weighting. Defaults to Z
#[builder(default)]
pub freqWeightingType: FreqWeighting,
/// FFT Length
pub nfft: usize,
/// Sampling frequency
pub fs: Flt,
}
impl ApsSettingsBuilder {
fn validate(&self) -> Result<()> {
if !self.fs.is_some() {
bail!("Sampling frequency not given");
}
let fs = self.fs.unwrap();
if !fs.is_normal() {
bail!("Sampling frequency not a normal number")
}
if fs <= 0.0 {
bail!("Invalid sampling frequency given as parameter");
}
if self.nfft.is_none() {
bail!("nfft not specified")
};
let nfft = self.nfft.unwrap();
if nfft % 2 != 0 {
bail!("NFFT should be even")
}
if nfft == 0 {
bail!("Invalid NFFT, should be > 0.")
}
// Perform some checks on ApsMode
if let Some(ApsMode::ExponentialWeighting { tau }) = self.mode {
if tau <= 0.0 {
bail!("Invalid time weighting constant [s]. Should be > 0 if given.");
}
}
Ok(())
}
}
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)]
impl ApsSettings {
#[new]
fn new(
mode: ApsMode,
overlap: Overlap,
windowType: WindowType,
freqWeightingType: FreqWeighting,
nfft: usize,
fs: Flt,
) -> ApsSettings {
ApsSettings {
mode,
overlap,
windowType,
freqWeightingType,
nfft,
fs,
}
}
}
impl ApsSettings {
/// Returns the amount of samples to keep in overlapping blocks of power
/// spectra.
pub fn get_overlap_keep(&self) -> usize {
self.validate_get_overlap_keep().unwrap()
}
/// Returns the amount of samples to `keep` in the time buffer when
/// overlapping time segments using [TimeBuffer].
fn validate_get_overlap_keep(&self) -> Result<usize> {
let nfft = self.nfft;
let overlap_keep = match self.overlap {
Overlap::Number { N } if N >= nfft => {
bail!("Invalid overlap number of samples. Should be < nfft, which is {nfft}.")
}
// Keep 1 sample, if overlap is 1 sample etc.
Overlap::Number { N } if N < nfft => N,
// If overlap percentage is >= 100, or < 0.0 its an error
Overlap::Percentage { pct } if !(0.0..100.).contains(&pct) => {
bail!("Invalid overlap percentage. Should be >= 0. And < 100.")
}
// If overlap percentage is 0, this gives
Overlap::Percentage { pct } => ((pct * nfft as Flt) / 100.) as usize,
Overlap::NoOverlap {} => 0,
_ => unreachable!(),
};
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_keep)
}
/// Return a reasonable acoustic default with a frequency resolution around
/// ~ 10 Hz, where nfft is still an integer power of 2.
///
/// # Errors
///
/// If `fs` is something odd, i.e. < 1 kHz, or higher than 1 MHz.
///
pub fn reasonableAcousticDefault(fs: Flt, mode: ApsMode) -> Result<ApsSettings> {
if fs < 1e3 || fs > 1e6 {
bail!("Sampling frequency for reasonable acoustic data is >= 1 kHz and <= 1 MHz.");
}
let fs_div_10_rounded = (fs / 10.) as u32;
// 2^30 is about 1 million. We search for a two-power of an nfft that is
// the closest to fs/10. The frequency resolution is about fs/nfft.
let nfft = (0..30).map(|i| 2u32.pow(i) - fs_div_10_rounded).fold(
// Start wth a value that is always too large
fs as u32 * 10,
|cur, new| cur.min(new),
) as usize;
Ok(ApsSettings {
mode,
fs,
nfft,
windowType: WindowType::default(),
overlap: Overlap::default(),
freqWeightingType: FreqWeighting::default(),
})
}
/// Return sampling frequency
pub fn fs(&self) -> Flt {
self.fs
}
/// Return Nyquist frequency
pub fn fnyq(&self) -> Flt {
self.fs / 2.
}
/// Returns a single-sided frequency array corresponding to points in Power
/// spectra computation.
pub fn getFreq(&self) -> Array1<Flt> {
let df = self.fs / self.nfft as Flt;
let K = self.nfft / 2 + 1;
Array1::linspace(0., (K - 1) as Flt * df, K)
}
}

5
src/ps/fft.rs
View File

@ -14,6 +14,11 @@ pub struct FFT {
// nfft stored as float, this is how it is required most often
nfftF: Flt,
}
impl Debug for FFT {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.debug_struct("Forward FFT engine, lenfth: {self.nfftF}").finish()
}
}
impl FFT {
/// Create new FFT from given nfft

3
src/ps/freqweighting.rs
View File

@ -25,6 +25,9 @@ impl FreqWeighting {
fn all() -> Vec<Self> {
Self::iter().collect()
}
fn __str__(&self) -> String {
format!("{self}-weighting")
}
#[staticmethod]
#[pyo3(name = "default")]
fn default_py() -> Self {

6
src/ps/mod.rs
View File

@ -10,12 +10,16 @@ mod ps;
mod timebuffer;
mod window;
mod freqweighting;
mod apssettings;
mod overlap;
mod apsmode;
use crate::config::*;
pub use freqweighting::FreqWeighting;
pub use overlap::Overlap;
pub use aps::{ApsSettings, ApsSettingsBuilder,ApsMode, AvPowerSpectra};
pub use apssettings::{ApsSettings, ApsSettingsBuilder};
pub use apsmode::ApsMode;
pub use aps::AvPowerSpectra;
pub use ps::{CrossPowerSpecra, PowerSpectra, CPSResult};
pub use window::{Window, WindowType};

1
src/ps/ps.rs
View File

@ -85,6 +85,7 @@ impl CrossPowerSpecra for CPSResult {
/// example the computations of spectrograms, or Welch' method of spectral
/// estimation.
///
#[derive(Debug)]
pub struct PowerSpectra {
/// Window used in estimator. The actual Window in here is normalized with
/// the square root of the Window power. This safes one division when

2
src/ps/timebuffer.rs
View File

@ -8,7 +8,7 @@ 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.
#[derive(Default)]
#[derive(Default, Debug)]
pub struct TimeBuffer {
data: Vec<VecDeque<Flt>>,
}

2
src/ps/window.rs
View File

@ -108,7 +108,7 @@ impl WindowType {
/// Window (taper) computed from specified window type.
#[derive(Clone)]
#[derive(Clone, Debug)]
pub struct Window {
/// The enum from which it is generated
pub w: WindowType,