Added loose gain value to SeriesBiquad. Implemented pink noise. Bup to 0.6.5

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "lasprs"
-version = "0.6.4"
+version = "0.6.5"
 edition = "2021"
 authors = ["J.A. de Jong <j.a.dejong@ascee.nl>"]
 description = "Library for Acoustic Signal Processing (Rust edition, with optional Python bindings via pyo3)"

src/filter/mod.rs

@@ -24,13 +24,14 @@ pub use dummy::DummyFilter;
 pub use seriesbiquad::SeriesBiquad;
 pub use zpkmodel::{PoleOrZero, ZPKModel, FilterSpec};
 
-/// Implementations of this trait are able to DSP-filter input data.
+/// Implementations of this trait are able to DSP-filter input data. The filter
+/// trait is implemented by, for example, [Biquad], [SeriesBiquad], and
+/// [BiquadBank].
 pub trait Filter: Send + Debug {
-    //! The filter trait is implemented by, for example, [Biquad], [SeriesBiquad], and [BiquadBank].
-
     /// Filter input to generate output. A vector of output floats is generated with the same
     /// length as input.
     fn filter(&mut self, input: &[Flt]) -> Vd;
+
     /// Reset the filter state(s). In essence, this makes sure that all memory of the past is
     /// forgotten.
     fn reset(&mut self);

src/filter/seriesbiquad.rs

@@ -13,6 +13,8 @@ use anyhow::{bail, Result};
 /// See (tests)
 ///
 pub struct SeriesBiquad {
+    // A loose gain value
+    gain: Flt,
     biqs: Vec<Biquad>,
 }
 
@@ -23,7 +25,7 @@ impl SeriesBiquad {
     fn __repr__(&self) -> String {
         format!("{self:?}")
     }
-    
+
     /// Create new series filter set. See [SeriesBiquad::new()]
     ///
     #[new]
@@ -36,7 +38,7 @@ impl SeriesBiquad {
     pub fn unit_py() -> SeriesBiquad {
         SeriesBiquad::unit()
     }
-    
+
     /// See: [SeriesBiquad::filter]
     #[pyo3(name = "filter")]
     pub fn filter_py<'py>(
@@ -59,7 +61,12 @@ impl SeriesBiquad {
     ///
     pub fn newFromBiqs(biqs: Vec<Biquad>) -> SeriesBiquad {
         assert!(biqs.len() > 0);
-        SeriesBiquad { biqs }
+        SeriesBiquad { biqs, gain: 1.0 }
+    }
+
+    /// Set an overall gain value to amplify / attenuate overall output
+    pub fn setGain(&mut self, g: Flt) {
+        self.gain = g;
     }
 
     /// Return reference to internally stored biquads
@@ -94,7 +101,7 @@ impl SeriesBiquad {
             bail!("No filter coefficients given!");
         }
 
-        Ok(SeriesBiquad { biqs })
+        Ok(SeriesBiquad { biqs, gain: 1.0 })
     }
 
     /// Unit impulse response series biquad. Input = output
@@ -107,6 +114,14 @@ impl SeriesBiquad {
     fn clone_dyn(&self) -> Box<dyn Filter> {
         Box::new(self.clone())
     }
+
+    /// Filter input signal in place
+    pub fn filter_inout(&mut self, inout: &mut [Flt]) {
+        for biq in self.biqs.iter_mut() {
+            biq.filter_inout(inout);
+        }
+        inout.iter_mut().for_each(|io| *io *= self.gain);
+    }
 }
 
 impl Filter for SeriesBiquad {
@@ -115,9 +130,7 @@ impl Filter for SeriesBiquad {
     //!
     fn filter(&mut self, input: &[Flt]) -> Vd {
         let mut inout = input.to_vec();
-        for biq in self.biqs.iter_mut() {
-            biq.filter_inout(&mut inout);
-        }
+        self.filter_inout(&mut inout);
         inout
     }
     fn reset(&mut self) {
@@ -130,7 +143,7 @@ impl Filter for SeriesBiquad {
 impl<'a, T: AsArray<'a, Flt>> TransferFunction<'a, T> for SeriesBiquad {
     fn tf(&self, fs: Flt, freq: T) -> Ccol {
         let freq = freq.into();
-        let mut res = self.biqs.first().unwrap().tf(fs, freq);
+        let mut res = self.biqs.first().unwrap().tf(fs, freq) * self.gain;
         for biq in self.biqs.iter().skip(1) {
             res = &res * biq.tf(fs, freq);
         }

src/siggen/mod.rs

@@ -21,6 +21,7 @@
 //!
 //! ```
 mod siggen;
+mod noise;
 mod siggenchannel;
 mod source;
 mod siggencmd;

src/siggen/noise.rs

@@ -0,0 +1,130 @@
+//! (Colored) noise Source implementations. Uses bilinear transform for a fixe and
+//! pre-calculated set of analogue poles and zeros.
+//!
+use super::source::SourceImpl;
+use crate::{config::*, PoleOrZero, SeriesBiquad, TransferFunction, ZPKModel};
+use rand::{rngs::SmallRng, Rng, SeedableRng};
+use rand_distr::StandardNormal;
+
+/// The number of poles in the filter that turns white noise into pink noise.
+const PINKNOISE_ANALOG_ORDER: usize = 10;
+
+/// White noise source. Can be colored by applying a color filter to the source
+#[derive(Clone, Debug)]
+pub struct WhiteNoise {
+    // SmallRng is a cheap random number generator
+    rng: SmallRng,
+}
+impl WhiteNoise {
+    pub fn new() -> Self {
+        WhiteNoise {
+            rng: SmallRng::from_entropy(),
+        }
+    }
+}
+impl SourceImpl for WhiteNoise {
+    fn genSignal_unscaled(&mut self, sig: &mut dyn ExactSizeIterator<Item = &mut Flt>) {
+        sig.for_each(|s| {
+            *s = self.rng.sample(StandardNormal);
+        });
+    }
+    fn reset(&mut self, _fs: Flt) {}
+    fn clone_dyn(&self) -> Box<dyn SourceImpl> {
+        Box::new(self.clone())
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct ColoredNoise {
+    // White noise generator
+    wn: WhiteNoise,
+    tmp: Vec<Flt>,
+    analogue_blueprint: ZPKModel,
+    filter: SeriesBiquad,
+}
+impl ColoredNoise {
+    /// Generate a colored noise signal source that outputs pink noise (-3 dB /
+    /// octave ) from 20 Hz to 20 kHz.
+    pub fn newPinkNoise() -> Self {
+        let twopi = 2. * pi;
+        let fl = 10.;
+        let fu = 20e3;
+        // Exponentially spread poles
+        let fpoles_real: Vec<Flt> = (0..PINKNOISE_ANALOG_ORDER)
+            .map(|i| fl * (fu / fl).powf((i as Flt) / (PINKNOISE_ANALOG_ORDER as Flt - 1.)))
+            .collect();
+
+        // Zeros as geometric means inbetween
+        let fzeros_real: Vec<Flt> = (0..PINKNOISE_ANALOG_ORDER - 1)
+            .map(|i| (fpoles_real[i] * fpoles_real[i + 1]).sqrt())
+            .collect();
+
+        let poles: Vec<PoleOrZero> = fpoles_real
+            .into_iter()
+            .map(|f| PoleOrZero::Real1(-twopi * f))
+            .collect();
+        let zeros: Vec<PoleOrZero> = fzeros_real
+            .into_iter()
+            .map(|f| PoleOrZero::Real1(-twopi * f))
+            .collect();
+        let gain = 1.;
+
+        let analogue_blueprint = ZPKModel::new(zeros, poles, gain);
+        let filter = analogue_blueprint.bilinear(480000.);
+        Self {
+            wn: WhiteNoise::new(),
+            tmp: vec![],
+            analogue_blueprint,
+            filter,
+        }
+    }
+}
+
+impl SourceImpl for ColoredNoise {
+    fn genSignal_unscaled(&mut self, sig: &mut dyn ExactSizeIterator<Item = &mut Flt>) {
+        // Make temporary array of same size as output iterator
+        self.tmp.resize(sig.len(), 0.);
+
+        // Generate white noise signal in place
+        self.wn.genSignal_unscaled(&mut self.tmp.iter_mut());
+
+        // Then, filter it with noise shaping filter
+        self.filter.filter_inout(&mut self.tmp);
+
+        // And copy over to output iterator
+        sig.zip(&self.tmp).for_each(|(sig, src)| {
+            *sig = *src;
+        });
+    }
+    fn reset(&mut self, fs: Flt) {
+        // Pre-computed analog poles
+        self.filter = self.analogue_blueprint.bilinear(fs);
+        let fnyq = fs / 2.;
+
+        // We set the gain such that the signal power is the same as that for
+        // the white noise that enters it. We do this by calculating the sum of
+        // the power in each bin up to the Nyquist frequency, and calculating a
+        // gain that is the inverse of this value.
+
+        // Choose a number of points. As long as its fine enough over the whole
+        // range (no peaks / dips are missing) At 48 kHz this means we sample
+        // the tf with a resolution of 10 Hz. Fine enough for slowly varying
+        // filter gains.
+        let N = 2400;
+        let freq = Array1::linspace(20., fnyq, 2000);
+        let tf_power_sum = self
+            .filter
+            .tf(fs, freq.view())
+            .mapv(|t| t.abs().powi(2))
+            .sum();
+        // The above is sum of squares. The white noise signal is unit
+        // variance, meaning also unit linear level. So to compute gain, we
+        // have to take the square root of this ratio.
+        let gain = ((N as Flt) / tf_power_sum).sqrt();
+        // dbg!(gain);
+        self.filter.setGain(gain);
+    }
+    fn clone_dyn(&self) -> Box<dyn SourceImpl> {
+        Box::new(self.clone())
+    }
+}

src/siggen/siggen.rs

@@ -21,6 +21,9 @@ use std::slice::IterMut;
 #[derive(Clone, Debug)]
 #[cfg_attr(feature = "python-bindings", pyclass)]
 pub struct Siggen {
+    // Sampling frequency in Hz
+    fs: Option<Flt>,
+
     // The source dynamic signal. Noise, a sine wave, sweep, etc
     source: Source,
 
@@ -71,6 +74,7 @@ impl Siggen {
     /// - `source` - Source function
     pub fn new(nchannels: usize, source: Source) -> Siggen {
         Siggen {
+            fs: None,
             source,
             channels: vec![SiggenChannelConfig::new(); nchannels],
             source_buf: vec![],
@@ -141,6 +145,8 @@ impl Siggen {
         match msg {
             SiggenCommand::ChangeSource { src } => {
                 self.source = src;
+
+                self.reset(self.fs.unwrap_or(48e3));
                 Ok(())
             }
             SiggenCommand::ResetSiggen { fs } => {
@@ -240,6 +246,7 @@ impl Siggen {
     /// * fs: (New) Sampling frequency \[Hz\]
     ///
     pub fn reset(&mut self, fs: Flt) {
+        self.fs = Some(fs);
         self.source.reset(fs);
         self.channels.iter_mut().for_each(|x| x.reset(fs))
     }

src/siggen/source.rs

@@ -1,5 +1,6 @@
 //! All sources for a signal generator. Sine waves, sweeps, noise, etc.
 use super::sweep::{SweepParams, SweepType};
+use super::noise::{ColoredNoise, WhiteNoise};
 use crate::config::*;
 use std::fmt::Debug;
 use std::ops::{Deref, DerefMut};
@@ -9,8 +10,6 @@ const twopi: Flt = 2.0 * pi;
 use crate::config::*;
 use anyhow::{bail, Result};
 use rand::prelude::*;
-use rand::rngs::ThreadRng;
-use rand_distr::StandardNormal;
 
 /// Signal source for a signal generator. A signal source is capable of creating
 /// new signal data.
@@ -42,7 +41,13 @@ impl Source {
     /// Create a white noise signal source
     pub fn newWhiteNoise() -> Source {
         Source {
-            src: Box::new(WhiteNoise { rng: SmallRng::from_entropy()}),
+            src: Box::new(WhiteNoise::new()),
+        }
+    }
+    /// Create a pink noise signal source
+    pub fn newPinkNoise() -> Source {
+        Source {
+            src: Box::new(ColoredNoise::newPinkNoise()),
         }
     }
 
@@ -89,6 +94,11 @@ impl Source {
         Self::newWhiteNoise()
     }
     #[staticmethod]
+    #[pyo3(name = "newPinkNoise")]
+    fn newPinkNoise_py() -> Source {
+        Self::newPinkNoise()
+    }
+    #[staticmethod]
     #[pyo3(name = "newSweep")]
     fn newSweep_py(
         fs: Flt,
@@ -119,23 +129,6 @@ impl SourceImpl for Silence {
         Box::new(self.clone())
     }
 }
-/// White noise source. Can be colored by applying a color filter to the source
-#[derive(Clone, Debug)]
-struct WhiteNoise {
-    // SmallRng is a cheap random number generator
-    rng: SmallRng
-}
-impl SourceImpl for WhiteNoise {
-    fn genSignal_unscaled(&mut self, sig: &mut dyn ExactSizeIterator<Item = &mut Flt>) {
-        sig.for_each(|s| {
-            *s = self.rng.sample(StandardNormal);
-        });
-    }
-    fn reset(&mut self, _fs: Flt) {}
-    fn clone_dyn(&self) -> Box<dyn SourceImpl> {
-        Box::new(self.clone())
-    }
-}
 
 /// Sine wave, with configurable frequency
 #[derive(Clone, Debug)]