Biguad, biquadbank, dummy filters, etc. A bit on all topics

This commit is contained in:
Anne de Jong 2024-07-03 20:01:12 +02:00
parent 7315939cbd
commit 7e9cf734d0
9 changed files with 188 additions and 140 deletions

View File

@ -10,7 +10,7 @@ use num::Complex;
/// This implementation only allows for normalized coefficients (a_0 = 1). It /// This implementation only allows for normalized coefficients (a_0 = 1). It
/// performs the following relation of output to input: /// performs the following relation of output to input:
/// ///
/// ``` /// ```math
/// y[n] = - a_1 * y[n-1] - a_2 * y[n-2] /// y[n] = - a_1 * y[n-1] - a_2 * y[n-2]
/// + b_0 * x[n] + b_1 * x[n-1] + b_2 * x[n-2] /// + b_0 * x[n] + b_1 * x[n-1] + b_2 * x[n-2]
/// ``` /// ```
@ -21,7 +21,7 @@ use num::Complex;
/// ///
/// The transfer function is: /// The transfer function is:
/// ///
/// ``` /// ```math
/// b_0 + b_1 z^-1 + b_2 * z^-2 /// b_0 + b_1 z^-1 + b_2 * z^-2
/// H[z] = ----------------------------- /// H[z] = -----------------------------
/// 1 + a_1 z^-1 + a_2 * z^-2 /// 1 + a_1 z^-1 + a_2 * z^-2
@ -62,11 +62,32 @@ impl Biquad {
#[pyo3(name = "firstOrderHighPass")] #[pyo3(name = "firstOrderHighPass")]
#[staticmethod] #[staticmethod]
/// See: [Biquad::firstOrderHighPass()] /// See: [Biquad::firstOrderHighPass()]
pub fn firstOrderHighPass_py(fs: Flt, cuton_Hz: Flt) -> PyResult<Biquad> { pub fn firstOrderHighPass_py(fs: Flt, fc: Flt) -> PyResult<Biquad> {
Ok(Biquad::firstOrderHighPass(fs, cuton_Hz)?) Ok(Biquad::firstOrderHighPass(fs, fc)?)
} }
#[pyo3(name = "filter")]
/// See: [Biquad::firstOrderLowPass()]
#[pyo3(name = "firstOrderLowPass")]
#[staticmethod]
pub fn firstOrderLowPass_py(fs: Flt, fc: Flt) -> PyResult<Biquad> {
Ok(Biquad::firstOrderLowPass(fs, fc)?)
}
/// See: [Biquad::tf()]
#[pyo3(name = "tf")]
pub fn tf_py<'py>(
&self,
py: Python<'py>,
fs: Flt,
freq: PyArrayLike1<Flt>,
) -> PyResult<PyArr1Cflt<'py>> {
let freq = freq.as_array();
let res = PyArray1::from_array_bound(py, &self.tf(fs, freq));
Ok(res)
}
/// See: [Biquad::filter()] /// See: [Biquad::filter()]
#[pyo3(name = "filter")]
pub fn filter_py<'py>( pub fn filter_py<'py>(
&mut self, &mut self,
py: Python<'py>, py: Python<'py>,
@ -151,22 +172,21 @@ impl Biquad {
/// First order low pass filter (one pole in the real axis). No pre-warping /// First order low pass filter (one pole in the real axis). No pre-warping
/// correction done. /// correction done.
///
/// * `fs` - Sampling frequency [Hz]
/// * `fc` - Cut-off frequency (-3 dB point) [Hz]
pub fn firstOrderLowPass(fs: Flt, fc: Flt) -> Result<Biquad> { pub fn firstOrderLowPass(fs: Flt, fc: Flt) -> Result<Biquad> {
match fc { if fc <= 0. {
x if fc <= 0. => bail!("Cuton frequency should be > 0"), bail!("Cuton frequency, given: should be > 0")
x if fc >= fs / 2. => bail!("Cuton frequency should be smaller than Nyquist frequency"),
_ => (),
} }
let w0: Flt = 2. * pi * fc / fs; if fc >= fs / 2. {
let cw = Flt::cos(w0); bail!("Cuton frequency should be smaller than Nyquist frequency")
let b0: Flt = 2. * pi * fc * (cw + 1.) / (2. * pi * fc * cw + 2. * pi * fc - cw + 1.); }
let b1: Flt = 2. * pi * fc * (cw + 1.) / (2. * pi * fc * cw + 2. * pi * fc - cw + 1.); let b0 = pi*fc/(pi*fc+fs);
let b2: Flt = 0.; let b1 = b0;
let a1: Flt = (2. * pi * fc * cw + 2. * pi * fc + cw - 1.) let a1 = (pi*fc-fs)/(pi*fc+fs);
/ (2. * pi * fc * cw + 2. * pi * fc - cw + 1.);
let a2: Flt = 0.;
Ok(Biquad::fromCoefs(b0, b1, b2, a1, a2)) Ok(Biquad::fromCoefs(b0, b1, 0., a1, 0.))
} }
/// Filter input signal, output by overwriting input slice. /// Filter input signal, output by overwriting input slice.
@ -204,8 +224,9 @@ impl Filter for Biquad {
Box::new(*self) Box::new(*self)
} }
} }
impl TransferFunction for Biquad { impl<'a, T: AsArray<'a, Flt>> TransferFunction<'a, T> for Biquad {
fn tf(&self, fs: Flt, freq: VdView) -> Ccol { fn tf(&self, fs: Flt, freq: T) -> Ccol {
let freq = freq.into();
let res = freq.mapv(|f| { let res = freq.mapv(|f| {
let z = Complex::exp(I * 2. * pi * f / fs); let z = Complex::exp(I * 2. * pi * f / fs);
let num = self.b0 + self.b1 / z + self.b2 / z / z; let num = self.b0 + self.b1 / z + self.b2 / z / z;
@ -232,15 +253,16 @@ mod test {
#[test] #[test]
fn test_firstOrderLowpass() { fn test_firstOrderLowpass() {
let fs = 10.; let fs = 1e5;
let fc = 1.; let fc = 10.;
let b = Biquad::firstOrderLowPass(fs, fc).unwrap(); let b = Biquad::firstOrderLowPass(fs, fc).unwrap();
let mut freq = Dcol::from_elem((3), 0.); let mut freq = Dcol::from_elem(5, 0.);
freq[1] = fc; freq[1] = fc;
freq[2] = fs/2.; freq[2] = fs / 2.;
let tf = b.tf(fs, freq.view()); let tf = b.tf(fs, freq.view());
assert_abs_diff_eq!(tf[0].re,1.); // println!("{:?}", tf);
assert_abs_diff_eq!(tf[0].im,0.); assert_abs_diff_eq!(tf[0].re, 1., epsilon = 1e-6);
assert_abs_diff_eq!(tf[1].abs(),1./Flt::sqrt(2.)); assert_abs_diff_eq!(tf[0].im, 0.);
assert_abs_diff_eq!(tf[1].abs(), 1. / Flt::sqrt(2.), epsilon = 1e-6);
} }
} }

View File

@ -66,8 +66,8 @@ impl BiquadBank {
&mut self, &mut self,
py: Python<'py>, py: Python<'py>,
input: PyArrayLike1<Flt>, input: PyArrayLike1<Flt>,
) -> PyResult<&'py PyArray1<Flt>> { ) -> PyResult<PyArr1Flt<'py>> {
Ok(self.filter(input.as_slice()?).into_pyarray(py)) Ok(self.filter(input.as_slice()?).into_pyarray_bound(py))
} }
#[pyo3(name = "reset")] #[pyo3(name = "reset")]
/// See: [BiquadBank::reset()] /// See: [BiquadBank::reset()]

24
src/filter/dummy.rs Normal file
View File

@ -0,0 +1,24 @@
use itertools::Itertools;
use super::*;
/// A Dummy fillter just does 'nothing', its input equals its output
#[derive(Clone, Copy, Debug)]
pub struct DummyFilter;
impl Filter for DummyFilter {
fn filter(&mut self, input: &[Flt]) -> Vec<Flt> {
// Just returns an allocated copy
input.to_vec()
}
fn reset(&mut self) { }
fn clone_dyn(&self) -> Box<dyn Filter> {
Box::new(*self)
}
}
impl<'a, T: AsArray<'a, Flt>> TransferFunction<'a, T> for DummyFilter {
fn tf(&self, _fs: Flt, freq: T) -> Ccol {
let freq = freq.into();
Ccol::ones(freq.len())
}
}

View File

@ -7,6 +7,7 @@ pub use super::config::*;
mod biquad; mod biquad;
mod biquadbank; mod biquadbank;
mod dummy;
mod seriesbiquad; mod seriesbiquad;
pub use biquad::Biquad; pub use biquad::Biquad;
@ -31,11 +32,14 @@ pub trait Filter: Send {
/// Implementations are able to generate transfer functions of itself /// Implementations are able to generate transfer functions of itself
pub trait TransferFunction: Send { pub trait TransferFunction<'a, T>: Send
where
T: AsArray<'a, Flt>,
{
/// Compute frequency response (i.e. transfer function from input to output) /// Compute frequency response (i.e. transfer function from input to output)
/// ///
/// Args /// Args
fn tf(&self, fs: Flt, freq: VdView) -> Ccol; fn tf(&self, fs: Flt, freq: T) -> Ccol;
} }
impl Clone for Box<dyn Filter> { impl Clone for Box<dyn Filter> {

View File

@ -22,11 +22,9 @@ pub struct SeriesBiquad {
#[cfg(feature = "python-bindings")] #[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)] #[cfg_attr(feature = "python-bindings", pymethods)]
impl SeriesBiquad { impl SeriesBiquad {
#[new]
/// Create new series filter set. See [SeriesBiquad::new()] /// Create new series filter set. See [SeriesBiquad::new()]
/// ///
#[new]
pub fn new_py<'py>(coefs: PyReadonlyArrayDyn<Flt>) -> PyResult<Self> { pub fn new_py<'py>(coefs: PyReadonlyArrayDyn<Flt>) -> PyResult<Self> {
Ok(SeriesBiquad::new(coefs.as_slice()?)?) Ok(SeriesBiquad::new(coefs.as_slice()?)?)
} }

View File

@ -1 +1,5 @@
//! Averaged power spectra module //! Averaged power spectra module. 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,

View File

@ -1,3 +1,5 @@
//! Power spectra estimator, that uses a Windowed FFT to estimate cross-power
//! spectra. Window functions are documented in the `window` module.
use crate::config::*; use crate::config::*;
use ndarray::parallel::prelude::*; use ndarray::parallel::prelude::*;
use num::pow::Pow; use num::pow::Pow;
@ -14,7 +16,10 @@ use std::mem::MaybeUninit;
use realfft::{RealFftPlanner, RealToComplex}; use realfft::{RealFftPlanner, RealToComplex};
/// Singlesided cross-Power spectra computation engine. /// Singlesided cross-Power spectra computation engine.
struct PowerSpectra { ///
/// Computes the signal(s) auto power and cross-power spectrum in each frequency
/// bin.
pub struct PowerSpectra {
// Window used in estimator // Window used in estimator
pub window: Window, pub window: Window,
// The window power, is corrected for in power spectra estimants // The window power, is corrected for in power spectra estimants
@ -100,6 +105,11 @@ impl PowerSpectra {
/// Compute cross power spectra from input time data. First axis is /// Compute cross power spectra from input time data. First axis is
/// frequency, second axis is channel i, third axis is channel j. /// frequency, second axis is channel i, third axis is channel j.
///
/// # Argument
///
/// * `tdata` - Input time data. This is a 2D array, where the first axis is
/// time and the second axis is the channel number.
pub fn compute<'a, T>(&mut self, tdata: T) -> Array3<Cflt> pub fn compute<'a, T>(&mut self, tdata: T) -> Array3<Cflt>
where where
T: Into<ArrayView<'a, Flt, Ix2>>, T: Into<ArrayView<'a, Flt, Ix2>>,
@ -130,9 +140,10 @@ impl PowerSpectra {
Zip::from(&mut out) Zip::from(&mut out)
.and(chi) .and(chi)
.and(chj) .and(chj)
.for_each(|out, chi, chjc| .for_each(|out, chi, chjc|{
// Loop over frequency components // Loop over frequency components
*out = 0.5 * chi * chjc *out = 0.5 * chi * chjc;
}
); );
// The DC component has no 0.5 correction, as it only // The DC component has no 0.5 correction, as it only
@ -275,12 +286,14 @@ mod test {
#[test] #[test]
fn test_parseval_with_window() { fn test_parseval_with_window() {
const nfft: usize = 48000; // A sufficiently high value is required here, to show that it works.
let rect = Window::new(WindowType::Hann, nfft); const nfft: usize = 2usize.pow(20);
let mut ps = PowerSpectra::newFromWindow(rect); let window = Window::new(WindowType::Hann, nfft);
// let window = Window::new(WindowType::Rect, nfft);
let mut ps = PowerSpectra::newFromWindow(window);
// Start with a time signal // Start with a time signal
let t: Dmat = Dmat::random((nfft, 1), StandardNormal); let t: Dmat = 2.*Dmat::random((nfft, 1), StandardNormal);
let tavg = t.sum()/(nfft as Flt); let tavg = t.sum()/(nfft as Flt);
let t_dc_power = tavg.powi(2); let t_dc_power = tavg.powi(2);
@ -295,7 +308,9 @@ mod test {
let fpower = power.sum().abs(); let fpower = power.sum().abs();
assert_ulps_eq!(t_dc_power, power[(0,0,0)].abs(), epsilon = Flt::EPSILON * (nfft as Flt).powi(2)); assert_ulps_eq!(t_dc_power, power[(0,0,0)].abs(), epsilon = Flt::EPSILON * (nfft as Flt).powi(2));
assert_ulps_eq!(signal_pwr, fpower, epsilon = Flt::EPSILON * (nfft as Flt).powi(2));
// This one fails when nfft is too short.
assert_ulps_eq!(signal_pwr, fpower, epsilon = 1e-2);
} }

View File

@ -1,56 +1,70 @@
//! Window functions designed for Welch' method. Implementations are given for common window //! Window functions designed for Welch' method. Implementations are given for
//! functions, as well as optimal 'jump' values `R`, that result in a certain overlap. //! the 5 classical window functions:
//!
//! * Rect - rectangular
//! * Hann - Von Hann window (sometimes wrongly called "Hanning")
//! * Bartlett
//! * Hamming
//! * Blackman
//!
#![allow(non_snake_case)] #![allow(non_snake_case)]
use crate::config::*; use crate::config::*;
#[macro_use] use strum_macros::Display;
use strum_macros::{Display};
fn linspace(nfft: usize) -> Dcol {
Dcol::linspace(0., nfft as Flt, nfft)
}
/// Von Hann window, often misnamed as the 'Hanning' window. /// Von Hann window, often misnamed as the 'Hanning' window.
fn hann(nfft: usize) -> (Dcol, usize) { fn hann(nfft: usize) -> Dcol {
let nfftF = nfft as Flt; let nfftF = nfft as Flt;
( Dcol::from_iter((0..nfft).map(|i| {
// The Window // 0.5*(1-cos(2*pi*i/(n-1)))
linspace(nfft).mapv(|i| (pi * i / (nfftF+1.)).sin().powi(2)), 0.5 * (1. - Flt::cos(2. * pi * i as Flt / (nfftF - 1.)))
// The hopp size }))
(nfft) / 2,
)
} }
fn rect(nfft: usize) -> (Dcol, usize) { /// Rectangular window
(Dcol::ones(nfft), nfft) fn rect(nfft: usize) -> Dcol {
Dcol::ones(nfft)
} }
fn blackman(N: usize) -> (Dcol, usize) { // Blackman window
let a0 = 7938. / 18608.; fn blackman(N: usize) -> Dcol {
let a1 = 9240. / 18608.; // Exact a0 coefficient, approximate value is 0.42
let a2 = 1430. / 18608.; const a0: Flt = 7938. / 18608.;
// Exact a1 coefficient, approximate value is 0.50
const a1: Flt = 9240. / 18608.;
// Exact a2 coefficient, approximate value is 0.08
const a2: Flt = 1430. / 18608.;
let Nf = N as Flt; let Nf = N as Flt;
let lin = linspace(N);
( Dcol::from_iter((0..N).map(|i| {
a0 - a1 * ((2. * pi / Nf) * lin.clone()).mapv(|x| x.cos()) let iF = i as Flt;
+ a2 * ((4. * pi / Nf) * lin).mapv(|x| x.cos()), a0 - a1 * Flt::cos(2. * pi * iF / (Nf - 1.)) + a2 * Flt::cos(4. * pi * iF / (Nf - 1.))
// The hop size }))
N / 3,
)
} }
fn bartlett(nfft: usize) -> (Dcol, usize) { fn bartlett(N: usize) -> Dcol {
let Nf = nfft as Flt; let Nf = N as Flt;
( Dcol::from_iter((0..N).map(|i| {
(1. - (2. * (linspace(nfft) - (Nf - 1.) / 2.) / Nf)).mapv(|x| x.abs()), let iF = i as Flt;
// The hop size if i <= (N - 1) / 2 {
nfft / 2, 2. * iF / (Nf - 1.)
) } else {
2. - 2. * iF / (Nf - 1.)
}
}))
} }
fn hamming(nfft: usize) -> (Dcol, usize) { fn hamming(N: usize) -> Dcol {
let alpha = 25.0 / 46.0; let Nf = N as Flt;
let beta = (1. - alpha) / 2.; // Approx 0.54
let Nf = nfft as Flt; const a0: Flt = 25.0 / 46.0;
( // Approx 0.46
alpha + 2. * beta * ((2. * pi / (Nf - 0.)) * linspace(nfft)).mapv(|x| x.sin()), const a1: Flt = 1. - a0;
// The hop size
(nfft) / 2, Dcol::from_iter((0..N).map(|i|
// Alessio et al.
a0 - a1 * Flt::cos(2. * pi * i as Flt / (Nf - 1.))
// end of map
)
// end from iter
) )
} }
@ -90,13 +104,14 @@ impl Window {
if nfft % 2 != 0 { if nfft % 2 != 0 {
panic!("Requires even nfft"); panic!("Requires even nfft");
} }
let (win, R) = match w { let win = match w {
WindowType::Hann => hann(nfft), WindowType::Hann => hann(nfft),
WindowType::Hamming => hamming(nfft), WindowType::Hamming => hamming(nfft),
WindowType::Rect => rect(nfft), WindowType::Rect => rect(nfft),
WindowType::Bartlett => bartlett(nfft), WindowType::Bartlett => bartlett(nfft),
WindowType::Blackman => blackman(nfft), WindowType::Blackman => blackman(nfft),
}; };
let R = nfft/2;
Window { w, win, R } Window { w, win, R }
} }
/// Convenience function that returns the length of the window. /// Convenience function that returns the length of the window.
@ -107,62 +122,20 @@ impl Window {
// //
#[cfg(test)] #[cfg(test)]
mod test { mod test {
use approx::assert_abs_diff_eq;
use super::*; use super::*;
#[test]
fn test_linspace() {
assert!(linspace(2)[0] == 0.);
// println!("{:?}", linspace(3));
assert!(linspace(3)[1] == 1.);
assert!(linspace(4).len() == 4);
}
#[test] #[test]
fn test_cola_hann() { fn test_windows(){
let nfft = 66; let Hann = hann(11);
let hann = Window::new(WindowType::Hann, nfft); let Hamming = hamming(11);
println!("{:?}", hann.win); let Bartlett = bartlett(11);
let Blackmann = bartlett(11);
let mut hanntot = Dcol::zeros(hann.len() * 4); // let h = hann(11);
assert_eq!(Hann[5] , 1.);
assert!(2 * hann.R == nfft); assert_eq!(Hamming[5] , 1.);
hanntot.slice_mut(s![0..nfft]).assign(&hann.win); assert_eq!(Bartlett[5] , 1.);
hanntot assert_eq!(Blackmann[5] , 1.);
.slice_mut(s![hann.R..nfft + hann.R])
.scaled_add(1.0, &hann.win);
hanntot
.slice_mut(s![nfft..2 * nfft])
.scaled_add(1.0, &hann.win);
hanntot
.slice_mut(s![nfft + hann.R..2 * nfft + hann.R])
.scaled_add(1.0, &hann.win);
hanntot
.slice_mut(s![2 * nfft..3 * nfft])
.scaled_add(1.0, &hann.win);
hanntot
.slice_mut(s![2 * nfft + hann.R..3 * nfft + hann.R])
.scaled_add(1.0, &hann.win);
println!("{:?}", hanntot);
}
#[test]
fn tets_cola_hamming() {
let nfft = 25;
let ham = Window::new(WindowType::Hamming, nfft);
let mut hamtot = Dcol::zeros(ham.len() * 3);
assert!(2 * ham.R == nfft);
hamtot.slice_mut(s![0..nfft]).scaled_add(1.0, &ham.win);
// println!("{:?}", hamtot);
hamtot
.slice_mut(s![ham.R..nfft + ham.R])
.scaled_add(1.0, &ham.win);
hamtot
.slice_mut(s![nfft..2 * nfft])
.scaled_add(1.0, &ham.win);
// println!("{:?}", hamtot);
// hantot.slice_mut(s![1+2*han1.R..nfft+1+2*han1.R]).scaled_add(1.0, &han2.win);
} }
} }

View File

@ -0,0 +1,8 @@
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 {
}