diff --git a/src/filter/mod.rs b/src/filter/mod.rs index f41b2cc..3903ca1 100644 --- a/src/filter/mod.rs +++ b/src/filter/mod.rs @@ -14,7 +14,7 @@ mod zpkmodel; mod butter; mod octave; -pub use super::ps::FreqWeightingType; +pub use super::ps::FreqWeighting; pub use biquad::Biquad; pub use biquadbank::BiquadBank; pub use octave::{StandardFilterDescriptor, G, FREQ_REF}; diff --git a/src/filter/octave.rs b/src/filter/octave.rs index 4f1f781..8f45dd7 100644 --- a/src/filter/octave.rs +++ b/src/filter/octave.rs @@ -1,6 +1,5 @@ use crate::{Flt, ZPKModel}; use anyhow::{anyhow, bail, Result}; -use clap::Error; use num::{traits::float, Float}; use rayon::iter::Filter; use softfloat::F64; @@ -67,11 +66,11 @@ pub const FREQ_REF: Flt = 1000.; /// /// ```rust /// use lasprs::filter::*; -/// fn main() -> anyhow::Result<()> { -/// let desc = StandardFilterDescriptor::Octave("16")?; -/// let filter = desc.genFilter().bilinear(48e3); -/// Ok(()) -/// } +/// # fn main() -> anyhow::Result<()> { +/// let desc = StandardFilterDescriptor::Octave("16")?; +/// let filter = desc.genFilter().bilinear(48e3); +/// # Ok(()) +/// # } /// ``` /// /// ## Create a one-third octave band bandpass filter that has the frequency of 42 in its pass-band @@ -79,11 +78,27 @@ pub const FREQ_REF: Flt = 1000.; /// ```rust /// /// use lasprs::filter::*; -/// fn main() -> anyhow::Result<()> { -/// let desc = StandardFilterDescriptor::filterForFreq(3, 42.)?; -/// let filter = desc.genFilter().bilinear(48e3); -/// Ok(()) -/// } +/// # fn main() -> anyhow::Result<()> { +/// let desc = StandardFilterDescriptor::filterForFreq(3, 42.)?; +/// let filter = desc.genFilter().bilinear(48e3); +/// # Ok(()) +/// # } +/// ``` +/// +/// ## Create a set of octave band filters +/// +/// ```rust +/// use lasprs::filter::*; +/// # fn main() -> anyhow::Result<()> { +/// // Generate custom set.. +/// let desc = StandardFilterDescriptor::genOctaveFilterSet("16", "16k").unwrap(); +/// // Or, when lazy: just generate the full set +/// let desc = StandardFilterDescriptor::fullOctaveFilterSet(); +/// let filters: Vec = desc.iter() +/// .map(|desc| desc.genFilter().bilinear(48e3)) +/// .collect(); +/// # Ok(()) +/// # } /// ``` #[derive(PartialEq, Clone, Debug)] pub struct StandardFilterDescriptor { @@ -222,38 +237,40 @@ impl StandardFilterDescriptor { } /// Creates a set of octave filters. - pub fn genOctaveFilterSet(low_f: Option, high_f: Option) -> Result> + pub fn genOctaveFilterSet(low_f: T, high_f: T) -> Result> where - T: TryInto, + T: TryInto, { - let xmin = if let Some(low_f) = low_f { - low_f.try_into()?.x - } else { - -OCTAVE_NAMES_OFFSET as i32 - }; - let xmax = if let Some(high_f) = high_f { - high_f.try_into()?.x - } else { - (OCTAVE_NOMINAL_MIDBAND_NAMES.len() - 1) as i32 - }; + let xmin = low_f.try_into()?.x; + let xmax = high_f.try_into()?.x; Ok((xmin..=xmax).map(|x| Self::Octave(x).unwrap()).collect()) } + /// Generate a full third octave bandpass filter set + pub fn fullThirdOctaveFilterSet() -> Vec { + Self::genThirdOctaveFilterSet( + THIRDOCTAVE_NOMINAL_MIDBAND_NAMES[0], + *(THIRDOCTAVE_NOMINAL_MIDBAND_NAMES.last().unwrap()), + ) + .unwrap() + } + + /// Generate a full octave bandpass filter set + pub fn fullOctaveFilterSet() -> Vec { + Self::genOctaveFilterSet( + OCTAVE_NOMINAL_MIDBAND_NAMES[0], + *(OCTAVE_NOMINAL_MIDBAND_NAMES.last().unwrap()), + ) + .unwrap() + } + /// Creates a set of one-third octave bandpass filters. - pub fn genThirdOctaveFilterSet(low_f: Option, high_f: Option) -> Result> + pub fn genThirdOctaveFilterSet(low_f: T, high_f: T) -> Result> where - T: TryInto, + T: TryInto, { - let xmin = if let Some(low_f) = low_f { - low_f.try_into()?.x - } else { - -THIRDOCTAVE_NAMES_OFFSET as i32 - }; - let xmax = if let Some(high_f) = high_f { - high_f.try_into()?.x - } else { - (THIRDOCTAVE_NOMINAL_MIDBAND_NAMES.len() - 1) as i32 - }; + let xmin = low_f.try_into()?.x; + let xmax = high_f.try_into()?.x; Ok((xmin..=xmax) .map(|x| Self::ThirdOctave(x).unwrap()) .collect()) diff --git a/src/filter/seriesbiquad.rs b/src/filter/seriesbiquad.rs index b899f22..fce0149 100644 --- a/src/filter/seriesbiquad.rs +++ b/src/filter/seriesbiquad.rs @@ -103,6 +103,7 @@ impl SeriesBiquad { SeriesBiquad::new(filter_coefs).unwrap() } + #[allow(dead_code)] fn clone_dyn(&self) -> Box { Box::new(self.clone()) } diff --git a/src/filter/zpkmodel.rs b/src/filter/zpkmodel.rs index 0c146bc..1518436 100644 --- a/src/filter/zpkmodel.rs +++ b/src/filter/zpkmodel.rs @@ -1,6 +1,6 @@ use super::butter::butter_lowpass_roots; use super::{Biquad, SeriesBiquad, TransferFunction}; -use crate::{config::*, ps::FreqWeightingType}; +use crate::{config::*, ps::FreqWeighting}; use itertools::{EitherOrBoth, Itertools}; use num::{zero, Complex}; use std::cmp::{max, min}; @@ -73,7 +73,7 @@ pub enum FilterSpec { /// poles should either be real, or come in complex conjugate pairs. This is /// enforced by the way the poles and zero's are internally stored. /// -#[derive(Clone, Debug, Default)] +#[derive(Clone, Debug)] #[cfg_attr(feature = "python-bindings", pyclass)] pub struct ZPKModel { // List of zeros @@ -87,6 +87,12 @@ pub struct ZPKModel { // transform to create digital filter of this analogue one. fwarp: Option, } +impl Default for ZPKModel { + fn default() -> Self { + ZPKModel{z:vec![], p:vec![], k: 1.0, fwarp: None} + } + +} impl<'a, T: AsArray<'a, Flt>> TransferFunction<'a, T> for ZPKModel { fn tf(&self, _fs: Flt, freq: T) -> Ccol { let freq = freq.into(); @@ -494,26 +500,26 @@ impl ZPKModel { /// /// # Args /// - /// - `wt` - `[FreqWeightingType]` to use. i.e. A-weighting. + /// - `wt` - `[FreqWeighting]` to use. i.e. A-weighting. /// /// # Examples /// /// ## Get part of pulse response of digital A-filter at 48 kHz /// /// ``` - /// use lasprs::filter::{ZPKModel, FreqWeightingType, Filter}; + /// use lasprs::filter::{ZPKModel, FreqWeighting, Filter}; /// /// // Sampling frequency in Hz /// let fs = 48000.; /// - /// let mut afilter = ZPKModel::freqWeightingFilter(FreqWeightingType::A).bilinear(fs); + /// let mut afilter = ZPKModel::freqWeightingFilter(FreqWeighting::A).bilinear(fs); /// let mut data = [0.; 1000]; /// data[0] = 1.0; /// let out = afilter.filter(&data); /// ``` /// - pub fn freqWeightingFilter(wt: FreqWeightingType) -> ZPKModel { - if let FreqWeightingType::Z = wt { + pub fn freqWeightingFilter(wt: FreqWeighting) -> ZPKModel { + if let FreqWeighting::Z = wt { return ZPKModel { z: vec![], p: vec![], @@ -548,15 +554,15 @@ impl ZPKModel { let p4 = 2. * pi * f4; let (zeros, poles) = match wt { - FreqWeightingType::Z => { + FreqWeighting::Z => { unreachable!() } - FreqWeightingType::C => { + FreqWeighting::C => { let zeros = vec![PoleOrZero::Real2(0., 0.)]; let poles = vec![PoleOrZero::Real2(-p1, -p1), PoleOrZero::Real2(-p4, -p4)]; (zeros, poles) } - FreqWeightingType::A => { + FreqWeighting::A => { let poles = vec![ PoleOrZero::Real2(-p1, -p1), PoleOrZero::Real2(-p2, -p3), diff --git a/src/ps/aps.rs b/src/ps/aps.rs index 3f2bdd8..0840976 100644 --- a/src/ps/aps.rs +++ b/src/ps/aps.rs @@ -4,7 +4,7 @@ use super::*; use crate::config::*; use anyhow::{bail, Error, Result}; use derive_builder::Builder; -use freqweighting::FreqWeightingType; +use freqweighting::FreqWeighting; /// All settings used for computing averaged power spectra using Welch' method. #[derive(Builder, Clone)] @@ -21,7 +21,7 @@ pub struct ApsSettings { windowType: WindowType, /// Kind of freqency weighting. Defaults to Z #[builder(default)] - freqWeightingType: FreqWeightingType, + freqWeightingType: FreqWeighting, /// FFT Length nfft: usize, /// Sampling frequency @@ -71,7 +71,7 @@ impl ApsSettings { self.validate_get_overlap_keep().unwrap() } /// Unpack all, returns parts in tuple - pub fn get(self) -> (ApsMode, Overlap, WindowType, FreqWeightingType, usize, Flt) { + pub fn get(self) -> (ApsMode, Overlap, WindowType, FreqWeighting, usize, Flt) { ( self.mode, self.overlap, @@ -134,7 +134,7 @@ impl ApsSettings { nfft, windowType: WindowType::default(), overlap: Overlap::default(), - freqWeightingType: FreqWeightingType::default(), + freqWeightingType: FreqWeighting::default(), }) } diff --git a/src/ps/fft.rs b/src/ps/fft.rs index 4bda503..e4585b7 100644 --- a/src/ps/fft.rs +++ b/src/ps/fft.rs @@ -11,11 +11,13 @@ pub struct FFT { timescratch: Vec, // rounded down nfft/2 half_nfft_rounded: usize, + // nfft stored as float, this is how it is required most often nfftF: Flt, } impl FFT { /// Create new FFT from given nfft + #[allow(dead_code)] pub fn newFromNFFT(nfft: usize) -> FFT { let mut planner = RealFftPlanner::::new(); let fft = planner.plan_fft_forward(nfft); diff --git a/src/ps/freqweighting.rs b/src/ps/freqweighting.rs index 156b160..2b57b37 100644 --- a/src/ps/freqweighting.rs +++ b/src/ps/freqweighting.rs @@ -1,7 +1,7 @@ use strum_macros::{Display, EnumMessage}; /// Sound level frequency weighting type (A, C, Z) #[derive(Display, Debug, EnumMessage, Default, Clone)] -pub enum FreqWeightingType { +pub enum FreqWeighting { /// A-weighting A, /// C-weighting @@ -15,9 +15,9 @@ mod test { use super::*; #[test] fn test() { - let a = FreqWeightingType::A; - let c = FreqWeightingType::C; - let z = FreqWeightingType::Z; + let a = FreqWeighting::A; + let c = FreqWeighting::C; + let z = FreqWeighting::Z; println!("A-weighting: {a}"); println!("C-weighting: {c}"); println!("Z-weighting: {z}"); diff --git a/src/ps/mod.rs b/src/ps/mod.rs index c91add0..e1e87bc 100644 --- a/src/ps/mod.rs +++ b/src/ps/mod.rs @@ -13,7 +13,7 @@ mod freqweighting; use crate::config::*; -pub use freqweighting::FreqWeightingType; +pub use freqweighting::FreqWeighting; pub use aps::{ApsSettings, ApsSettingsBuilder,ApsMode, AvPowerSpectra, Overlap}; pub use ps::{CrossPowerSpecra, PowerSpectra, CPSResult}; pub use window::{Window, WindowType}; diff --git a/src/slm/mod.rs b/src/slm/mod.rs index fd2370d..58f39a0 100644 --- a/src/slm/mod.rs +++ b/src/slm/mod.rs @@ -2,10 +2,70 @@ //! //! Provides structs and helpers (SLMBuilder) for creating configurated Sound //! Level Meters. -//! +//! +//! # Usage examples +//! +//! ## Simple over-all level SLM +//! +//! This example creates a simple SLM that processes at 48kHz, uses fast time +//! weighting and A frequency weighting: +//! +//! ``` +//! # use anyhow::Result; +//! use lasprs::slm::*; +//! use lasprs::filter::StandardFilterDescriptor; +//! use ndarray::Array1; +//! # fn main() -> Result<()> { +//! +//! // Generate an overall filter (no filter at all) +//! let desc = StandardFilterDescriptor::Overall().unwrap(); +//! +//! // Generate settings +//! let settings = SLMSettingsBuilder::default() +//! .fs(48e3) +//! .freqWeighting(FreqWeighting::A) +//! .timeWeighting(TimeWeighting::Fast) +//! .filterDescriptors(&[desc]).build().unwrap(); +//! +//! let mut slm = SLM::new(settings); +//! // Generate some data. Yes, this is not the most spectacular set +//! let mut data = Array1::zeros(48000); +//! data[0] = 1.; +//! +//! // Now apply some data. This is a kind of the SLM-s impulse response +//! let res = slm.run(&data.as_slice().unwrap()).unwrap(); +//! +//! // Only one channel of result data +//! assert_eq!(res.len(), 1); +//! +//! let res = &res[0]; +//! println!("Data is: {res:#?}"); +//! +//! // Get the equivalent level: +//! let Leq = slm.Leq()[0]; +//! +//! # Ok::<() , anyhow::Error>(()) +//! # } +//! ``` +//! +//! ## One-third octave band, plus overall +//! +//! ``` +//! use lasprs::filter::StandardFilterDescriptor; +//! // Generate the default set of one-third octave band filters +//! let mut desc = StandardFilterDescriptor::fullThirdOctaveFilterSet(); +//! desc.push(StandardFilterDescriptor::Overall().unwrap()); +//! +//! // Rest of code is the same as in previous example +//! +//! ``` +//! mod settings; mod tw; mod slm; pub use slm::SLM; -pub use settings::SLMSettings; -pub use tw::TimeWeightingType; \ No newline at end of file +pub use settings::{SLMSettings, SLMSettingsBuilder}; +pub use tw::TimeWeighting; +pub use crate::ps::FreqWeighting; + +const SLM_MAX_CHANNELS: usize = 64; \ No newline at end of file diff --git a/src/slm/settings.rs b/src/slm/settings.rs index 64aa6f3..481b5cd 100644 --- a/src/slm/settings.rs +++ b/src/slm/settings.rs @@ -1,19 +1,46 @@ +use super::{TimeWeighting, SLM_MAX_CHANNELS}; +use crate::{filter::StandardFilterDescriptor, Flt, FreqWeighting}; +use anyhow::Result; +use clap::builder; use derive_builder::Builder; -use smallvec::SmallVec; -use crate::{Flt, FreqWeightingType, filter::StandardFilterDescriptor}; -use super::TimeWeightingType; - - +use smallvec::{smallvec, SmallVec}; +/// Settings used to create a Sound Level Meter. #[derive(Builder, Clone)] +#[builder(setter(into))] pub struct SLMSettings { + /// Sampling frequency in \[Hz\] pub fs: Flt, + /// Reference level, in units of measured quantity. For sound pressure in + /// air, this is typically 20 μPa. This is also the default value. + #[builder(default = "2e-5")] pub Lref: Flt, - pub freqWeighting: FreqWeightingType, - pub timeWeighting: TimeWeightingType, - pub filterDescriptors: SmallVec<[StandardFilterDescriptor; 64]>, - + /// Frequency weightin A/C/Z applied to data. Defaults to [FreqWeighting::default()]. + #[builder(default)] + pub freqWeighting: FreqWeighting, + /// For time-dependent output, the time weighthing applied (Fast / Slow) + pub timeWeighting: TimeWeighting, + /// Descriptors for the filters, maximum of 64, which is a reasonable amount + /// and - if all used - might already choke a computer. + pub filterDescriptors: Vec, } -impl SLMSettings { - -} \ No newline at end of file +#[cfg(test)] +mod test { + use super::*; + use anyhow::Result; + + #[test] + fn test_slmsettings1() -> Result<()> { + let desc = StandardFilterDescriptor::genFilterSetInRange(1, 100., 5e3, true).unwrap(); + + let _ = SLMSettingsBuilder::default() + .fs(1e3) + .freqWeighting(FreqWeighting::A) + .timeWeighting(TimeWeighting::Slow) + .filterDescriptors(desc) + .build() + .unwrap(); + + Ok(()) + } +} diff --git a/src/slm/slm.rs b/src/slm/slm.rs index 6045053..08c61e5 100644 --- a/src/slm/slm.rs +++ b/src/slm/slm.rs @@ -5,23 +5,30 @@ use rayon::iter::{IntoParallelRefIterator, ParallelIterator}; use rayon::prelude::*; use smallvec::SmallVec; -use super::settings::SLMSettings; +use super::{settings::SLMSettings, SLM_MAX_CHANNELS}; use crate::{config::*, filter::Filter}; -use crate::{Biquad, Dcol, Flt, FreqWeightingType, PoleOrZero, SeriesBiquad, ZPKModel}; +use crate::{Biquad, Dcol, Flt, FreqWeighting, PoleOrZero, SeriesBiquad, ZPKModel}; +#[derive(Debug, Clone)] struct SLMChannel { + // Statistics to store stat: SLMStat, + // The bandpass filter for this channel bp: SeriesBiquad, + // The rectifier filter rect_lowpass_up: Biquad, + // The rectifier filter for decreasing values. Only for asymmetric time + // weighting (case of impulse weighting: [TimeWeighting::Impulse]) rect_lowpass_down: Option, } /// Sound Level Meter +#[derive(Debug, Clone)] pub struct SLM { // Number of samples processed after last run() is called. N: usize, Lrefsq: Flt, prefilter: SeriesBiquad, - channels: SmallVec<[SLMChannel; 64]>, + channels: SmallVec<[SLMChannel; SLM_MAX_CHANNELS]>, } impl SLM { @@ -179,7 +186,7 @@ impl SLM { } } -#[derive(Clone, Default)] +#[derive(Debug, Clone, Default)] /// Quantities defined as powers, i.e. square of amplitude struct SLMStat { // Max signal power @@ -192,3 +199,45 @@ struct SLMStat { // Last obtained signal power, after last time run() is called. Pt_last: Flt, } + +#[cfg(test)] +mod test { + use crate::{ + siggen::Siggen, + slm::{SLMSettingsBuilder, TimeWeighting}, + Flt, FreqWeighting, StandardFilterDescriptor, + }; + + use super::SLM; + + #[test] + fn test_slm1() { + const fs: Flt = 48e3; + const N: usize = (fs/8.) as usize; + + let desc = StandardFilterDescriptor::Overall().unwrap(); + + let settings = SLMSettingsBuilder::default() + .fs(fs) + .timeWeighting(TimeWeighting::Fast) + .freqWeighting(FreqWeighting::Z) + .filterDescriptors(&[desc]) + .build() + .unwrap(); + + let mut siggen = Siggen::newSine(1, 1000.); + siggen.setAllMute(false); + siggen.reset(fs); + let mut data = vec![0.; N]; + siggen.genSignal(&mut data); + // println!("{:#?}", data); + + let mut slm = SLM::new(settings); + // println!("{slm:#?}"); + let res = slm.run(&data).unwrap(); + let res = &res[0]; + println!("{slm:#?}"); + println!("{:#?}", &res[res.len()- 100..]); + + } +} diff --git a/src/slm/tw.rs b/src/slm/tw.rs index 6540ae6..5391b9b 100644 --- a/src/slm/tw.rs +++ b/src/slm/tw.rs @@ -1,14 +1,16 @@ use crate::Flt; +/// Time weighting to use in level detection of Sound Level Meter. #[derive(Clone, Copy)] -pub enum TimeWeightingType { - /// Slow time weighting +pub enum TimeWeighting { + /// Slow time weighting ~ 1 s Slow, - /// Fast time weighting + /// Fast time weighting ~ 1/8 s Fast, - /// Impulse time weighting + /// Impulse time weighting ~ 30 ms Impulse, /// A custom symmetric time weighting CustomSymmetric { + /// Custom time constant [s] t: Flt, }, /// A custom symmetric time weighting @@ -19,14 +21,23 @@ pub enum TimeWeightingType { tdown: Flt, }, } -impl TimeWeightingType { +impl Default for TimeWeighting { + fn default() -> Self { + TimeWeighting::Fast + } +} +impl TimeWeighting { /// get the analog poles of the single pole lowpass filter required for - /// getting the 'rectified' level (detector phase of SLM). + /// getting the 'rectified' level (detection phase of SLM). pub fn getLowpassPoles(&self) -> (Flt, Option) { - use TimeWeightingType::*; + use TimeWeighting::*; match self { Slow => (-1.0, None), - Fast => (-1. / 8., None), + Fast => + // Time constant is 1/8 s, pole is at -8 rad/s + { + (-8., None) + } Impulse => { // For the impulse time weighting, some source says ~ 2.9 dB/s // drop for the decay @@ -41,7 +52,7 @@ impl TimeWeightingType { // with 10*log10(exp(-1.0/tau)) = -10/ln(10)/tau ≅ -4.34/tau // dB/s where ln denotes the natural logarithm. So suppose we // have 1.5 s, we indeed get a decay rate of 2.9 dB/s - (-35e-3, Some(-1.5)) + (-1. / 35e-3, Some(-1. / 1.5)) } CustomSymmetric { t } => { assert!(*t > 0.); diff --git a/tools/watchdoc.sh b/tools/watchdoc.sh index fac447f..a08ed14 100755 --- a/tools/watchdoc.sh +++ b/tools/watchdoc.sh @@ -6,6 +6,6 @@ # $ cargo install cargo-watch cargo-docserver` # ``` # -cargo watch -s "clear && cargo rustdoc -p lasprs --lib && cargo docserve" +cargo watch -s "clear && cargo doc --no-deps -p lasprs --lib && cargo docserve" # Then open: ${browser} http://localhost:4000