Added octave band filter design code. First steps into SLM. Still needs proper testing

Cargo.toml

@@ -16,7 +16,7 @@ crate-type = ["cdylib", "rlib",]
 
 [dependencies]
 # Error handling
-anyhow = "1.0.75"
+anyhow = "1.0.86"
 
 # Numerics
 # Optional future feature for ndarray: blas
@@ -90,6 +90,12 @@ parking_lot = "0.12.3"
 # Builder code
 derive_builder = "0.20.0"
 
+# Stack-allocated vectors
+smallvec = "1.13.2"
+
+# Compile time constant floating point operations
+softfloat = "1.0.0"
+
 [dev-dependencies]
 ndarray-rand = "0.14.0"
 

src/config.rs

@@ -36,7 +36,7 @@ if #[cfg(feature = "python-bindings")] {
 }
 }
 pub use ndarray::prelude::*;
-pub use ndarray::{Array1, Array2, ArrayView1};
+pub use ndarray::{Array1, Array2, ArrayView1, ArrayViewMut1};
 
 pub use ndarray::Zip;
 use num::complex::Complex;

src/filter/biquad.rs

@@ -30,7 +30,11 @@ use num::Complex;
 ///  And the frequency response can be found by filling in in above equation z =
 ///  exp(i*omega/fs), where fs is the sampling frequency and omega is the radian
 ///  frequency at which the transfer function is evaluated.
+///  
+/// ## Implementation details
 ///
+/// The implementaion is so-called "Direct-form 2", see
+/// [https://en.wikipedia.org/wiki/Digital_biquad_filter].
 pub struct Biquad {
     // State parameters
     w1: Flt,
@@ -124,6 +128,29 @@ impl Biquad {
         }
     }
 
+    /// Re-initialize state. *This is an advanced function. You should know what
+    /// you are doing!*. If not, please use any other function like
+    /// [Biquad::reset].
+    pub fn setNextOutputX0(&mut self, out: Flt) {
+        let (b0, b1, b2, a1, a2) = (self.b0, self.b1, self.b2, self.a1, self.a2);
+        let w = out / (b1 + b2 - b0 * (a1 + a2));
+        self.w1 = w;
+        self.w2 = w;
+    }
+
+    /// Change the gain value such that it matches `val` at frequency `freq`.
+    /// Does not change the phase at the given frequency.
+    pub fn setGainAt(mut self, freq: Flt, required_gain: Flt) -> Biquad {
+        assert!(required_gain > 0.);
+        let freq = [freq];
+        let cur_gain_at_freq = self.tf(-1.0, &freq)[0].abs();
+        let gain_fac = required_gain / cur_gain_at_freq;
+        self.b0 *= gain_fac;
+        self.b1 *= gain_fac;
+        self.b2 *= gain_fac;
+        self
+    }
+
     /// Construct a Biquad with 0 initial state from coefficients given as
     /// arguments.
     ///
@@ -202,17 +229,22 @@ impl Biquad {
         Ok(Biquad::fromCoefs(b0, b1, 0., a1, 0.))
     }
 
+    #[inline]
+    /// Filter single sample, outputs by overwriting input sample.
+    pub fn filter_inout_single(&mut self, sample: &mut Flt) {
+        let w0 = *sample - self.a1 * self.w1 - self.a2 * self.w2;
+        let yn = self.b0 * w0 + self.b1 * self.w1 + self.b2 * self.w2;
+        self.w2 = self.w1;
+        self.w1 = w0;
+        *sample = yn;
+    }
+
     /// Filter input signal, output by overwriting input slice.
     #[inline]
     pub fn filter_inout(&mut self, inout: &mut [Flt]) {
         for sample in inout.iter_mut() {
-            let w0 = *sample - self.a1 * self.w1 - self.a2 * self.w2;
+            self.filter_inout_single(sample);
-            let yn = self.b0 * w0 + self.b1 * self.w1 + self.b2 * self.w2;
-            self.w2 = self.w1;
-            self.w1 = w0;
-            *sample = yn;
         }
-        // println!("{:?}", inout);
     }
 
     /// Create new biquad using bilinear transform. Optionally pre-warps the
@@ -470,4 +502,21 @@ mod test {
         // let freq = &[0., 10.,100.,1000., 2000.];
         // println!("{:?}", b3.tf(fs, freq));
     }
+
+    #[test]
+    fn test_setOutput1() {
+        let mut f = Biquad::firstOrderHighPass(10., 1.).unwrap();
+        f.setNextOutputX0(1.0);
+        let mut sample = 0.;
+        f.filter_inout_single(&mut sample);
+        assert_abs_diff_eq!(sample, 1.0);
+    }
+    #[test]
+    fn test_setOutput2() {
+        let mut f = Biquad::bilinear_zpk(1.0, None, Some(PoleOrZero::Real1(-1.)), None, None);
+        f.setNextOutputX0(4.2);
+        let mut sample = 0.;
+        f.filter_inout_single(&mut sample);
+        assert_abs_diff_eq!(sample, 4.2, epsilon = 1e-6);
+    }
 }

src/filter/mod.rs

@@ -12,10 +12,12 @@ mod dummy;
 mod seriesbiquad;
 mod zpkmodel;
 mod butter;
+mod octave;
 
 pub use super::ps::FreqWeightingType;
 pub use biquad::Biquad;
 pub use biquadbank::BiquadBank;
+pub use octave::{StandardFilterDescriptor, G, FREQ_REF};
 pub use dummy::DummyFilter;
 pub use seriesbiquad::SeriesBiquad;
 pub use zpkmodel::{PoleOrZero, ZPKModel, FilterSpec};

src/filter/octave.rs

@@ -0,0 +1,542 @@
+use crate::{Flt, ZPKModel};
+use anyhow::{anyhow, bail, Result};
+use clap::Error;
+use num::{traits::float, Float};
+use rayon::iter::Filter;
+use softfloat::F64;
+use std::{borrow::Cow, cmp::Ordering};
+
+/// Names of standard octave filters
+const OCTAVE_NOMINAL_MIDBAND_NAMES: [&str; 12] = [
+    "8", "16", "31.5", "63", "125", "250", "500", "1k", "2k", "4k", "8k", "16k",
+];
+const OCTAVE_NAMES_OFFSET: i32 = 7;
+
+const MIN_MIDBAND_FREQ: Flt = 8.;
+const MAX_MIDBAND_FREQ: Flt = 20e3;
+
+const THIRDOCTAVE_NOMINAL_MIDBAND_NAMES: [&'static str; 33] = [
+    "12.5", "16", "20", "25", "31.5", "40", "50", "63", "80", "100", "125", "160", "200", "250",
+    "315", "400", "500", "630", "800", "1k", "1.25k", "1.6k", "2k", "2.5k", "3.15k", "4k", "5k",
+    "6.3k", "8k", "10k", "12.5k", "16k", "20k",
+];
+const THIRDOCTAVE_NAMES_OFFSET: i32 = 19;
+
+/// Return the num x-value for a certain 'name', like '16', or '1k'
+fn nominal_octave_designator(name: &str) -> Result<i32> {
+    debug_assert!(OCTAVE_NOMINAL_MIDBAND_NAMES[OCTAVE_NAMES_OFFSET as usize] == "1k");
+    Ok(OCTAVE_NOMINAL_MIDBAND_NAMES
+        .iter()
+        .position(|i| *i == name)
+        .ok_or(anyhow!(
+            "Cannot find name in list of OCTAVE_NOMINAL_MIDBAND_NAMES"
+        ))? as i32
+        - OCTAVE_NAMES_OFFSET)
+}
+
+fn nominal_thirdoctave_designator(name: &str) -> Result<i32> {
+    debug_assert!(THIRDOCTAVE_NOMINAL_MIDBAND_NAMES[THIRDOCTAVE_NAMES_OFFSET as usize] == "1k");
+    Ok(THIRDOCTAVE_NOMINAL_MIDBAND_NAMES
+        .iter()
+        .position(|i| *i == name)
+        .ok_or(anyhow!(
+            "Cannot find name in list of THIRDOCTAVE_NOMINAL_MIDBAND_NAMES"
+        ))? as i32
+        - THIRDOCTAVE_NAMES_OFFSET)
+}
+// Raise a^b. In const-mode.
+const fn powf(a: Flt, b: Flt) -> Flt {
+    let a = a as f64;
+    let b = b as f64;
+    let a = softfloat::F64::from_native_f64(a);
+    let b = softfloat::F64::from_native_f64(b);
+    softfloat::F64::exp(b.mul(a.ln())).to_native_f64() as Flt
+}
+
+/// Octave ratio. We use G_10, which is 10^(3/10) ≅ 1.995
+pub const G: Flt = powf(10., 0.3);
+/// Reference freuqency, 1kHz
+pub const FREQ_REF: Flt = 1000.;
+
+/// Standard filter descriptor. Used to generate bandpass filters that are
+/// compliant with IEC 61260 (1995).
+///
+/// # Examples
+///
+/// ## Create a 16 Hz octave band digital filter running at 48kHz.
+///
+/// ```rust
+/// use lasprs::filter::*;
+/// 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
+///
+/// ```rust
+///
+/// use lasprs::filter::*;
+/// fn main() -> anyhow::Result<()> {
+///     let desc = StandardFilterDescriptor::filterForFreq(3, 42.)?;
+///     let filter = desc.genFilter().bilinear(48e3);
+///     Ok(())
+/// }
+/// ```
+#[derive(PartialEq, Clone, Debug)]
+pub struct StandardFilterDescriptor {
+    /// b and x. Bandwidth and offset w.r.t. reference frequency.
+    ///
+    /// Band width fraction of an octave. 1 means full octave of bandwidth. 3
+    /// means 1/3th octave, 6 means 1/6th, and so on.
+    ///
+    /// If bx is None, it means we do not filter at all (an overall channel)
+    b: u32,
+    x: i32,
+}
+
+impl StandardFilterDescriptor {
+    /// Create analog filter specification from descriptor
+    pub fn genFilter(&self) -> ZPKModel {
+        let order = 5;
+        if let Some((fl, fu)) = self.fl_fh() {
+            ZPKModel::butter(crate::FilterSpec::Bandpass { fl, fu, order })
+        } else {
+            ZPKModel::default()
+        }
+    }
+    // Check whether a certain midband frequency of created
+    // StandardFilterDescriptor is in the allowed range. This is a helper
+    // function that is used to check wheter created StandardFilterDescriptors
+    // are valid.
+    fn check_fmid_in_range(&self) -> Result<()> {
+        if let Some(fm) = self.fm() {
+            if fm < MIN_MIDBAND_FREQ {
+                bail!(
+                    "Invalid x. Computed filter center frequency is {} Hz, which is too low. Lowest allowed is {} Hz",
+                    fm, MIN_MIDBAND_FREQ
+                )
+            } else if fm > 20e3 {
+                bail!(
+                    "Invalid x. Computed filter center frequency is {} Hz, which is too high. Highest allowed is {} Hz",
+                    fm, MAX_MIDBAND_FREQ
+                )
+            }
+        }
+
+        Ok(())
+    }
+    /// Create new standard filter descriptor `b` from given relative bandwidth
+    /// and band designator `x`. If not sure what `x` and `b` are, see
+    /// documentation on [StandardFilterDescriptor::genFilterSetByDesignator].
+    pub fn build(b: u32, x: i32) -> Result<StandardFilterDescriptor> {
+        let desc = StandardFilterDescriptor { b, x };
+        desc.check_fmid_in_range()?;
+        match b {
+            0 => Ok(desc),
+            1 => Ok(desc),
+            3 => Ok(desc),
+            6 => Ok(desc),
+            12 => Ok(desc),
+            24 => Ok(desc),
+            _ => bail!(
+                "Bandwidth {} is invalid. Please choose a value from 0, 1, 3, 6, 12 or 24",
+                b
+            ),
+        }
+    }
+    /// Generate filter descriptor. Practically applies no filtering at all.
+    pub fn Overall() -> Result<StandardFilterDescriptor> {
+        Ok(StandardFilterDescriptor { b: 0, x: 0 })
+    }
+    /// Generate filter descriptor for octave band.
+    ///
+    /// # Args
+    ///
+    /// - `band_descr` - band designator. Can be '1k', or 0.
+    pub fn Octave<T>(band_descr: T) -> Result<StandardFilterDescriptor>
+    where
+        T: TryInto<OctaveBandDescriptor, Error = anyhow::Error>,
+    {
+        let x = band_descr.try_into()?.x;
+        Ok(StandardFilterDescriptor { b: 1, x })
+    }
+    /// Generate filter descriptor for one-third octave band.
+    ///
+    /// # Args
+    ///
+    /// - `x` - band designator
+    pub fn ThirdOctave<T>(band_descr: T) -> Result<StandardFilterDescriptor>
+    where
+        T: TryInto<ThirdOctaveBandDescriptor, Error = anyhow::Error>,
+    {
+        let x = band_descr.try_into()?.x;
+        Ok(StandardFilterDescriptor { b: 3, x })
+    }
+
+    /// Searches for a filter with `1/b` relative bandwidth w.r.t one octave
+    /// that has frequency `f` in its pass-band.
+    ///
+    pub fn filterForFreq(b: u32, f: Flt) -> Result<StandardFilterDescriptor> {
+        if f < MIN_MIDBAND_FREQ || f > MAX_MIDBAND_FREQ {
+            bail!("Invalid frequency. Please use search frequency between 8 Hz and 20 kHz")
+        }
+        match b {
+            0 => Self::Overall(),
+            1 | 3 | 6 | 12 | 24 => {
+                let mut desc = StandardFilterDescriptor { b, x: 0 };
+
+                let f_in_range = |desc: &StandardFilterDescriptor| {
+                    let (fl, fh) = desc.fl_fh().unwrap();
+                    // println!("fl: {fl}, fh: {fh}");
+                    if f < fl {
+                        Ordering::Less
+                    } else if f > fh {
+                        Ordering::Greater
+                    } else {
+                        Ordering::Equal
+                    }
+                };
+
+                loop {
+                    let fm = desc.fm().unwrap();
+                    // eprintln!("Fmid: {fm:.2e}");
+                    // eprintln!("desc: {desc:#?}");
+                    let ord = f_in_range(&desc);
+                    // Bands for midband frequencies are a bit wider here
+                    if fm < MIN_MIDBAND_FREQ - 3. || fm > MAX_MIDBAND_FREQ * 1.1 {
+                        bail!("Frequency not in range");
+                    }
+                    match ord {
+                        Ordering::Equal => break,
+                        Ordering::Less => desc = StandardFilterDescriptor { b, x: desc.x - 1 },
+                        Ordering::Greater => desc = StandardFilterDescriptor { b, x: desc.x + 1 },
+                    }
+                }
+                Ok(desc)
+            }
+            _ => Self::build(b, 0),
+        }
+    }
+
+    /// Creates a set of octave filters.
+    pub fn genOctaveFilterSet<T>(low_f: Option<T>, high_f: Option<T>) -> Result<Vec<Self>>
+    where
+        T: TryInto<OctaveBandDescriptor, Error = Error>,
+    {
+        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
+        };
+        Ok((xmin..=xmax).map(|x| Self::Octave(x).unwrap()).collect())
+    }
+
+    /// Creates a set of one-third octave bandpass filters.
+    pub fn genThirdOctaveFilterSet<T>(low_f: Option<T>, high_f: Option<T>) -> Result<Vec<Self>>
+    where
+        T: TryInto<ThirdOctaveBandDescriptor, Error = Error>,
+    {
+        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
+        };
+        Ok((xmin..=xmax)
+            .map(|x| Self::ThirdOctave(x).unwrap())
+            .collect())
+    }
+
+    /// Generate a filter set using designators
+    ///
+    /// # Args
+    ///
+    /// - `b` - Inverse of the relative bandwidth w.r.t. one octave. `b=0` means
+    ///   overall, `b=1` is one octave, `b=3`` is one-third, etc.
+    /// - `xmin` - Band designator of lowest band. Midband frequency can be computed as [FREQ_REF]*[G]^(`xmin/b`)
+    /// - `xmax` - Band designator of lowest band. Midband frequency can be computed as [FREQ_REF]*[G]^(`xmax/b`)
+    /// - `include_overall` - If `true`, adds an overall filter (a no-op) as the last designator in the list
+    pub fn genFilterSetByDesignator(
+        b: u32,
+        xmin: i32,
+        xmax: i32,
+        include_overall: bool,
+    ) -> Result<Vec<Self>> {
+        if xmin > xmax {
+            bail!("xmin should be <= xmax");
+        }
+        let cap = (xmax - xmin) as usize + if include_overall { 1 } else { 0 };
+        let mut res = Vec::with_capacity(cap);
+
+        for x in xmin..=xmax {
+            res.push(StandardFilterDescriptor::build(b, x)?);
+        }
+
+        if include_overall {
+            res.push(StandardFilterDescriptor::Overall()?)
+        }
+
+        Ok(res)
+    }
+
+    /// Convenience function for creating a filter bank. Creates a set of
+    /// standard filters with relative bandwidth `b`, that has `fl` in the
+    /// lowest bandpass filter and `fu` in the highest.
+    ///
+    /// # Other args
+    ///
+    /// - `include_overall` - If `true`, adds an overall filter (a no-op) as the
+    pub fn genFilterSetInRange(
+        b: u32,
+        fl: Flt,
+        fu: Flt,
+        include_overall: bool,
+    ) -> Result<Vec<Self>> {
+        let xmin = StandardFilterDescriptor::filterForFreq(b, fl)?.x;
+        let xmax = StandardFilterDescriptor::filterForFreq(b, fu)?.x;
+        StandardFilterDescriptor::genFilterSetByDesignator(b, xmin, xmax, include_overall)
+    }
+
+    /// Returns the midband frequency in \[Hz\]
+    pub fn fm(&self) -> Option<Flt> {
+        if self.b == 0 {
+            None
+        } else {
+            let b = self.b as Flt;
+            let x = self.x as Flt;
+            Some(FREQ_REF * G.powf(x / b))
+        }
+    }
+
+    /// Cuton frequency and cut-off frequency, in \[Hz\].
+    /// Returns none if it does not apply, for [FilterDescriptor::Overall].
+    pub fn fl_fh(&self) -> Option<(Flt, Flt)> {
+        match self.b {
+            0 => None,
+            b => {
+                let fm = self.fm().unwrap();
+                let b = b as Flt;
+                let fl = fm * G.powf(-1. / (2. * b));
+                let fu = fm * G.powf(1. / (2. * b));
+                Some((fl, fu))
+            }
+        }
+    }
+    /// Give a common name to filter, specifically the filters are named after
+    /// the midband frequency.
+    pub fn name(&self) -> Cow<'static, str> {
+        let x = self.x;
+        match self.b {
+            0 => Cow::Borrowed("Overall"),
+            1 => OctaveBandDescriptor { x }.name(),
+            3 => ThirdOctaveBandDescriptor { x }.name(),
+            6 => {
+                if x % 2 == 0 {
+                    ThirdOctaveBandDescriptor { x: x / 2 }.name()
+                } else {
+                    Default::default()
+                }
+            }
+            12 => {
+                if x % 2 == 0 {
+                    StandardFilterDescriptor {
+                        b: self.b / 2,
+                        x: self.x / 2,
+                    }
+                    .name()
+                } else {
+                    Default::default()
+                }
+            }
+            _ => unreachable!(),
+        }
+    }
+}
+
+/// A valid descriptor for a standard one-third octave band
+pub struct ThirdOctaveBandDescriptor {
+    x: i32,
+}
+/// A valid descriptor for a standard octave band
+pub struct OctaveBandDescriptor {
+    x: i32,
+}
+
+impl TryFrom<i32> for OctaveBandDescriptor {
+    type Error = anyhow::Error;
+    fn try_from(x: i32) -> Result<Self, Self::Error> {
+        if x + OCTAVE_NAMES_OFFSET < 0
+            || x + OCTAVE_NAMES_OFFSET >= OCTAVE_NOMINAL_MIDBAND_NAMES.len() as i32
+        {
+            bail!(
+                "Invalid filter designator x. Should be >= -{OCTAVE_NAMES_OFFSET} and < {}",
+                OCTAVE_NOMINAL_MIDBAND_NAMES.len() as i32 - OCTAVE_NAMES_OFFSET
+            );
+        }
+        Ok(Self { x })
+    }
+}
+impl TryFrom<&str> for OctaveBandDescriptor {
+    type Error = anyhow::Error;
+    fn try_from(name: &str) -> Result<Self, Self::Error> {
+        Ok(Self {
+            x: nominal_octave_designator(name)?,
+        })
+    }
+}
+impl TryFrom<i32> for ThirdOctaveBandDescriptor {
+    type Error = anyhow::Error;
+    fn try_from(x: i32) -> Result<Self, Self::Error> {
+        if x + THIRDOCTAVE_NAMES_OFFSET < 0
+            || x + THIRDOCTAVE_NAMES_OFFSET >= THIRDOCTAVE_NOMINAL_MIDBAND_NAMES.len() as i32
+        {
+            bail!(
+                "Invalid filter designator x. Should be >= -{THIRDOCTAVE_NAMES_OFFSET} and < {}",
+                THIRDOCTAVE_NOMINAL_MIDBAND_NAMES.len() as i32 - THIRDOCTAVE_NAMES_OFFSET
+            );
+        }
+        Ok(Self { x })
+    }
+}
+impl TryFrom<&str> for ThirdOctaveBandDescriptor {
+    type Error = anyhow::Error;
+    fn try_from(name: &str) -> Result<Self, Self::Error> {
+        Ok(Self {
+            x: nominal_thirdoctave_designator(name)?,
+        })
+    }
+}
+
+trait BandDescriptor {
+    fn name(&self) -> Cow<'static, str>;
+}
+impl BandDescriptor for OctaveBandDescriptor {
+    fn name(&self) -> Cow<'static, str> {
+        Cow::Borrowed(
+            OCTAVE_NOMINAL_MIDBAND_NAMES
+                .get((self.x + OCTAVE_NAMES_OFFSET) as usize)
+                .map(|s| *s)
+                .unwrap_or_default(),
+        )
+    }
+}
+impl BandDescriptor for ThirdOctaveBandDescriptor {
+    fn name(&self) -> Cow<'static, str> {
+        Cow::Borrowed(
+            THIRDOCTAVE_NOMINAL_MIDBAND_NAMES
+                .get((self.x + THIRDOCTAVE_NAMES_OFFSET) as usize)
+                .map(|s| *s)
+                .unwrap_or_default(),
+        )
+    }
+}
+
+// impl Into<OctaveBandDescriptor> for &str {
+//     fn into(self) -> OctaveBandDescriptor {
+//         OctaveBandDescriptor{x: self}
+//     }
+// }
+#[cfg(test)]
+mod test {
+    use super::*;
+    use approx::assert_abs_diff_eq;
+
+    #[test]
+    fn test_finder() {
+        // assert_eq!(
+        //     StandardFilterDescriptor::filterForFreq(0, 1000.).unwrap(),
+        //     StandardFilterDescriptor::Overall().unwrap()
+        // );
+        // assert_eq!(
+        //     StandardFilterDescriptor::filterForFreq(1, 1e3).unwrap(),
+        //     StandardFilterDescriptor::Octave(0).unwrap()
+        // );
+        assert_eq!(
+            StandardFilterDescriptor::filterForFreq(1, 8.).unwrap(),
+            StandardFilterDescriptor::Octave(-OCTAVE_NAMES_OFFSET).unwrap()
+        );
+        // assert_eq!(
+        //     StandardFilterDescriptor::filterForFreq(3, 1000.).unwrap(),
+        //     StandardFilterDescriptor::ThirdOctave(0).unwrap()
+        // );
+        // assert_eq!(
+        //     StandardFilterDescriptor::filterForFreq(3, 12.).unwrap(),
+        //     StandardFilterDescriptor::ThirdOctave(-THIRDOCTAVE_NAMES_OFFSET).unwrap()
+        // );
+    }
+
+    #[test]
+    fn test_builders() {
+        assert_eq!(
+            StandardFilterDescriptor::Octave("8").unwrap(),
+            StandardFilterDescriptor::Octave(-OCTAVE_NAMES_OFFSET).unwrap()
+        );
+        assert_eq!(
+            StandardFilterDescriptor::Octave("2k").unwrap(),
+            StandardFilterDescriptor::Octave(1).unwrap()
+        );
+        assert_eq!(
+            StandardFilterDescriptor::ThirdOctave("12.5").unwrap(),
+            StandardFilterDescriptor::ThirdOctave(-THIRDOCTAVE_NAMES_OFFSET).unwrap()
+        );
+        assert_eq!(
+            StandardFilterDescriptor::ThirdOctave("2k").unwrap(),
+            StandardFilterDescriptor::ThirdOctave(3).unwrap()
+        );
+    }
+    #[test]
+    #[should_panic]
+    fn out_range_octave1() {
+        StandardFilterDescriptor::Octave("4").unwrap();
+    }
+    #[test]
+    #[should_panic]
+    fn out_range_octave2() {
+        StandardFilterDescriptor::Octave("7").unwrap();
+    }
+
+    #[test]
+    fn test_name_and_approx() {
+        assert_eq!(
+            StandardFilterDescriptor::filterForFreq(1, 16e3).unwrap(),
+            StandardFilterDescriptor::Octave("16k").unwrap(),
+        );
+    }
+
+    #[test]
+    fn test_names() {
+        assert_eq!(StandardFilterDescriptor::Octave(1).unwrap().name(), "2k");
+        assert_eq!(
+            StandardFilterDescriptor::ThirdOctave(1).unwrap().name(),
+            "1.25k"
+        );
+    }
+
+    #[test]
+    fn test_octave() {
+        assert_eq!(nominal_octave_designator("1k").unwrap(), 0);
+        assert_eq!(nominal_octave_designator("2k").unwrap(), 1);
+    }
+    #[test]
+    fn test_thirdoctave() {
+        assert_eq!(nominal_thirdoctave_designator("1k").unwrap(), 0);
+        assert_eq!(nominal_thirdoctave_designator("2k").unwrap(), 3);
+    }
+
+    #[test]
+    fn test_G() {
+        assert_abs_diff_eq!(G, (10 as Flt).powf(3. / 10.));
+    }
+}

src/filter/seriesbiquad.rs

@@ -62,6 +62,11 @@ impl SeriesBiquad {
         SeriesBiquad { biqs }
     }
 
+    /// Return reference to internally stored biquads
+    pub fn getBiquads(&self) -> &Vec<Biquad> {
+        &self.biqs
+    }
+
     /// Create a new series biquad, having an arbitrary number of biquads.
     ///
     /// # Arguments

src/filter/zpkmodel.rs

@@ -176,9 +176,10 @@ impl ZPKModel {
     /// - `poles` - list like struct of poles. Can be a `Vec<ZeroOrPole>` or an
     ///   `&[ZeroOrPole]`.
     /// - `k` - linear gain.
-    pub fn new<T>(zeros: T, poles: T, k: Flt) -> ZPKModel
+    pub fn new<T, U>(zeros: T, poles: U, k: Flt) -> ZPKModel
     where
         T: Into<Vec<PoleOrZero>>,
+        U: Into<Vec<PoleOrZero>>,
     {
         let z = zeros.into();
         let p = poles.into();

src/lib.rs

@@ -41,6 +41,7 @@ pub mod ps;
 pub mod siggen;
 use filter::*;
 pub mod rt;
+pub mod slm;
 
 /// A Python module implemented in Rust.
 #[cfg(feature = "python-bindings")]

src/slm/mod.rs

@@ -1,15 +1,11 @@
 //! Sound Level Meter (SLM) module.
-//! 
+//!
 //! Provides structs and helpers (SLMBuilder) for creating configurated Sound
 //! Level Meters.
-//! 
+//!
-
+mod settings;
-
+mod tw;
-/// Sound Level Meter
+mod slm;
-struct SLM {
+pub use slm::SLM;
-
+pub use settings::SLMSettings;
-}
+pub use tw::TimeWeightingType;
-
-impl SLM {
-
-}

src/slm/settings.rs

@@ -0,0 +1,19 @@
+use derive_builder::Builder;
+use smallvec::SmallVec;
+use crate::{Flt, FreqWeightingType, filter::StandardFilterDescriptor};
+use super::TimeWeightingType;
+
+
+#[derive(Builder, Clone)]
+pub struct SLMSettings {
+    pub fs: Flt,
+    pub Lref: Flt,
+    pub freqWeighting: FreqWeightingType,
+    pub timeWeighting: TimeWeightingType,
+    pub filterDescriptors: SmallVec<[StandardFilterDescriptor; 64]>,
+
+}
+
+impl  SLMSettings {
+    
+}

src/slm/slm.rs

@@ -0,0 +1,194 @@
+use derive_builder::Builder;
+use itertools::Itertools;
+use ndarray::ArrayView1;
+use rayon::iter::{IntoParallelRefIterator, ParallelIterator};
+use rayon::prelude::*;
+use smallvec::SmallVec;
+
+use super::settings::SLMSettings;
+use crate::{config::*, filter::Filter};
+use crate::{Biquad, Dcol, Flt, FreqWeightingType, PoleOrZero, SeriesBiquad, ZPKModel};
+struct SLMChannel {
+    stat: SLMStat,
+    bp: SeriesBiquad,
+    rect_lowpass_up: Biquad,
+    rect_lowpass_down: Option<Biquad>,
+}
+
+/// Sound Level Meter
+pub struct SLM {
+    // Number of samples processed after last run() is called.
+    N: usize,
+    Lrefsq: Flt,
+    prefilter: SeriesBiquad,
+    channels: SmallVec<[SLMChannel; 64]>,
+}
+
+impl SLM {
+    // Create simple first order lowpass filter with unit D.C. gain and given
+    // real pole.
+    fn lpfilter_from_pole(fs: Flt, p: PoleOrZero) -> Biquad {
+        Biquad::bilinear_zpk(fs, None, Some(p), Some(1.0), None).setGainAt(0., 1.)
+    }
+    /// Create new Sound Level Meter from given settings
+    pub fn new(settings: SLMSettings) -> Self {
+        let fs = settings.fs;
+        let prefilter = ZPKModel::freqWeightingFilter(settings.freqWeighting).bilinear(fs);
+        let channels = settings
+            .filterDescriptors
+            .iter()
+            .map(|descriptor| {
+                // Generate bandpass filter
+                let bp = descriptor.genFilter().bilinear(fs);
+                // Initalize statistics with defaults
+                let stat = SLMStat::default();
+
+                // Generate rectifier filter for upwards
+                let poles = settings.timeWeighting.getLowpassPoles();
+
+                let rect_lowpass_up = Self::lpfilter_from_pole(fs, PoleOrZero::Real1(poles.0));
+
+                let rect_lowpass_down = if let Some(p) = poles.1 {
+                    Some(Self::lpfilter_from_pole(fs, PoleOrZero::Real1(p)))
+                } else {
+                    None
+                };
+                SLMChannel {
+                    stat,
+                    bp,
+                    rect_lowpass_up,
+                    rect_lowpass_down,
+                }
+            })
+            .collect();
+        SLM {
+            prefilter,
+            channels,
+            Lrefsq: settings.Lref.powi(2),
+            N: 0,
+        }
+    }
+    /// Push new time data through sound level meter. Returns L(t) data for each
+    /// channel.
+    ///
+    /// # Args
+    ///
+    /// - `td`: Time data
+    pub fn run(&mut self, td: &[Flt]) -> Option<Vec<Vec<Flt>>> {
+        if td.len() == 0 {
+            return None;
+        }
+        let prefiltered = self.prefilter.filter(td);
+        let Lt_iter = self.channels.par_iter_mut().map(|ch| {
+            let mut tmp = ch.bp.filter(&prefiltered);
+            let mut N = self.N;
+
+            // Filtered squared
+            let mut filtered_squared = {
+                let mut tmp_view = ArrayViewMut1::from(&mut tmp);
+                tmp_view.mapv_inplace(|a| a * a);
+                tmp_view
+            };
+
+            // Update Lpk, Leq
+            filtered_squared.for_each(|sample_pwr| {
+                let new_pk = sample_pwr.abs();
+                if new_pk > ch.stat.Ppk {
+                    ch.stat.Ppk = new_pk;
+                }
+                // Update equivalent level
+                ch.stat.Peq = (ch.stat.Peq * N as Flt + sample_pwr) / (N as Flt + 1.);
+                N += 1;
+            });
+
+            // Run filtered_squared signal throug rectifier
+            if let Some(rectifier_down) = &mut ch.rect_lowpass_down {
+                filtered_squared.mapv_inplace(|sample_sq| {
+                    let mut fup = sample_sq;
+                    let mut fdown = sample_sq;
+
+                    // Filter in up-filter
+                    let rectifier_up = &mut ch.rect_lowpass_up;
+                    rectifier_up.filter_inout_single(&mut fup);
+                    // Filter in down-filter
+                    rectifier_down.filter_inout_single(&mut fdown);
+
+                    // Check who wins
+                    if fup > fdown {
+                        rectifier_down.setNextOutputX0(fup);
+                        fup
+                    } else {
+                        rectifier_up.setNextOutputX0(fdown);
+                        fdown
+                    }
+                });
+            } else {
+                // Filter in place
+                let rectifier = &mut ch.rect_lowpass_up;
+                rectifier.filter_inout(filtered_squared.as_slice_mut().unwrap());
+            }
+
+            // Update max signal power gotten so far
+            let rectified = &mut filtered_squared;
+            rectified.for_each(|val| {
+                if *val > ch.stat.Pmax {
+                    ch.stat.Pmax = *val;
+                }
+            });
+            // Update last signal power coming from SLM
+            ch.stat.Pt_last = *filtered_squared.last().unwrap();
+            tmp
+        });
+        let Lt: Vec<_> = Lt_iter.collect();
+        self.N += td.len();
+        Some(Lt)
+    }
+
+    /// Number of channels in SLM
+    pub fn nch(&self) -> usize {
+        self.channels.len()
+    }
+
+    fn levels_from<T>(&self, stat_returner: T) -> Dcol
+    where
+        T: Fn(&SLMChannel) -> Flt,
+    {
+        Dcol::from_iter(
+            self.channels
+                .iter()
+                .map(|ch| 20. * Flt::log10(stat_returner(ch) / self.Lrefsq)),
+        )
+    }
+
+    /// Get max levels for each channel
+    pub fn Lmax(&self) -> Dcol {
+        self.levels_from(|ch| ch.stat.Pmax)
+    }
+    /// Get peak levels for each channel
+    pub fn Lpk(&self) -> Dcol {
+        self.levels_from(|ch| ch.stat.Ppk)
+    }
+
+    /// Get equivalent levels for each channel
+    pub fn Leq(&self) -> Dcol {
+        self.levels_from(|ch| ch.stat.Peq)
+    }
+    /// Get last value of level vs time
+    pub fn Ltlast(&self) -> Dcol {
+        self.levels_from(|ch| ch.stat.Pt_last)
+    }
+}
+
+#[derive(Clone, Default)]
+/// Quantities defined as powers, i.e. square of amplitude
+struct SLMStat {
+    // Max signal power
+    Pmax: Flt,
+    // Peak signal power
+    Ppk: Flt,
+    // Equivalent signal power
+    Peq: Flt,
+
+    // Last obtained signal power, after last time run() is called.
+    Pt_last: Flt,
+}

src/slm/tw.rs

@@ -0,0 +1,57 @@
+use crate::Flt;
+#[derive(Clone, Copy)]
+pub enum TimeWeightingType {
+    /// Slow time weighting
+    Slow,
+    /// Fast time weighting
+    Fast,
+    /// Impulse time weighting
+    Impulse,
+    /// A custom symmetric time weighting
+    CustomSymmetric {
+        t: Flt,
+    },
+    /// A custom symmetric time weighting
+    CustomAsymmetric {
+        /// Time weighting when level is increasing
+        tup: Flt,
+        /// Time weighting when level is decreasing
+        tdown: Flt,
+    },
+}
+impl TimeWeightingType {
+    /// get the analog poles of the single pole lowpass filter required for
+    /// getting the 'rectified' level (detector phase of SLM).
+    pub fn getLowpassPoles(&self) -> (Flt, Option<Flt>) {
+        use TimeWeightingType::*;
+        match self {
+            Slow => (-1.0, None),
+            Fast => (-1. / 8., None),
+            Impulse => {
+                // For the impulse time weighting, some source says ~ 2.9 dB/s
+                // drop for the decay
+                // [https://www.nti-audio.com/en/support/know-how/fast-slow-impulse-time-weighting-what-do-they-mean].
+                //
+                // Other source
+                // [https://support.dewesoft.com/en/support/solutions/articles/14000139949-exponential-averaging-fast-f-slow-s-impulse-i-]
+                // say a time constant of 1.5 s. Are they compatible?
+
+                // Compute decay rate in dB/s from the filter time constant. An
+                // initial value drops as exp(-t/tau). So in 1 s the level drops
+                // 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))
+            }
+            CustomSymmetric { t } => {
+                assert!(*t > 0.);
+                (-*t, None)
+            }
+            CustomAsymmetric { tup, tdown } => {
+                assert!(*tup > 0.);
+                assert!(*tdown > 0.);
+                (-*tup, Some(-*tdown))
+            }
+        }
+    }
+}