A bit of refactoring in checks, added documentation and examples of ZPKModel

This commit is contained in:
Anne de Jong 2024-08-11 13:47:24 +02:00
parent 99a8db23b8
commit c7d2cfc43f
4 changed files with 148 additions and 104 deletions

View File

@ -13,6 +13,7 @@ mod seriesbiquad;
mod zpkmodel; mod zpkmodel;
mod butter; mod butter;
pub use super::ps::FreqWeightingType;
pub use biquad::Biquad; pub use biquad::Biquad;
pub use biquadbank::BiquadBank; pub use biquadbank::BiquadBank;
pub use dummy::DummyFilter; pub use dummy::DummyFilter;

View File

@ -1,13 +1,12 @@
use std::cmp::{max, min};
use super::butter::butter_lowpass_roots; use super::butter::butter_lowpass_roots;
use super::{Biquad, SeriesBiquad, TransferFunction};
use crate::{config::*, ps::FreqWeightingType};
use itertools::{EitherOrBoth, Itertools}; use itertools::{EitherOrBoth, Itertools};
use num::{zero, Complex}; use num::{zero, Complex};
use std::cmp::{max, min};
use crate::config::*; /// Reasonable maximum order of Butterworth filters.
pub const BUTTER_MAX_ORDER: u32 = 40;
use super::{Biquad, SeriesBiquad, TransferFunction};
/// Specification of a filter for a certain type. /// Specification of a filter for a certain type.
/// ///
/// The order corresponds to the rolloff in dB/decade. order=1 means 20 /// The order corresponds to the rolloff in dB/decade. order=1 means 20
@ -48,10 +47,15 @@ pub enum FilterSpec {
/// to generate analog filters of standard type, i.e. bandpass, lowpass and /// to generate analog filters of standard type, i.e. bandpass, lowpass and
/// highpass. These can subsequentially be used to generate a digital filter. /// highpass. These can subsequentially be used to generate a digital filter.
/// ///
/// # Example ///
/// # Example: Create a digital second order Butterworth bandpass filter
/// ///
/// ```rust /// ```rust
/// use lasprs::filter::{FilterSpec, ZPKModel}; /// use lasprs::filter::{FilterSpec, ZPKModel};
/// let fs = 48000.;
///
/// let butter = ZPKModel::butter(FilterSpec::Bandpass{fl: 10., fu: 100., order: 2});
/// let mut filt = butter.bilinear(fs);
/// ///
/// ``` /// ```
/// ///
@ -230,8 +234,12 @@ impl ZPKModel {
ZPKModel { z, p, k, fwarp } ZPKModel { z, p, k, fwarp }
} }
// Set critical frequency in filter /// Set critical frequency in filter, used for bilinear transform.
fn setWarpFreq(mut self, fcrit: Flt) -> ZPKModel { ///
/// # Args
///
/// - `fcrit` - New critical frequency in \[Hz\].
pub fn setWarpFreq(mut self, fcrit: Flt) -> ZPKModel {
self.fwarp = Some(fcrit); self.fwarp = Some(fcrit);
self self
} }
@ -346,12 +354,46 @@ impl ZPKModel {
} }
} }
/// Create a Butterworth filter according to a certain specification fn check_spec(spec: FilterSpec) {
pub fn butter(spec: FilterSpec) -> ZPKModel { let check_fc = |fc| assert!(fc > 0., "Cut-off frequency should be > 0");
let check_order = |order| {
assert!(
order > 0 && order <= BUTTER_MAX_ORDER,
"Invalid filter order"
)
};
match spec {
FilterSpec::Lowpass { fc, order } => {
check_fc(fc);
check_order(order);
}
FilterSpec::Highpass { fc, order } => {
check_fc(fc);
check_order(order);
}
FilterSpec::Bandpass { fl, fu, order } => {
assert!(
fl <= fu && fl > 0.,
"Invalid cut-on and cut-off frequency specified"
);
check_order(order);
}
}
}
/// Create a Butterworth filter according to a certain specification in
/// `spec`.
///
/// # Panics
///
/// - If specified `order == 0`
/// - If order is larger than [BUTTER_MAX_ORDER].
/// - If for a bandpass filter `fl>=fu`.
/// - If `fl`, or `fu < 0`.
pub fn butter(spec: FilterSpec) -> ZPKModel {
Self::check_spec(spec);
match spec { match spec {
FilterSpec::Lowpass { fc, order } => { FilterSpec::Lowpass { fc, order } => {
assert!(fc > 0.);
assert!(order > 0);
let p = butter_lowpass_roots(fc, order as u32).collect(); let p = butter_lowpass_roots(fc, order as u32).collect();
let z = vec![]; let z = vec![];
ZPKModel { ZPKModel {
@ -365,8 +407,6 @@ impl ZPKModel {
.setWarpFreq(fc) .setWarpFreq(fc)
} }
FilterSpec::Highpass { fc, order } => { FilterSpec::Highpass { fc, order } => {
assert!(fc > 0.);
assert!(order > 0);
let p = butter_lowpass_roots(fc, order as u32).collect(); let p = butter_lowpass_roots(fc, order as u32).collect();
let z = vec![PoleOrZero::Real1(0.); order as usize]; let z = vec![PoleOrZero::Real1(0.); order as usize];
ZPKModel { ZPKModel {
@ -380,8 +420,6 @@ impl ZPKModel {
.setWarpFreq(fc) .setWarpFreq(fc)
} }
FilterSpec::Bandpass { fl, fu, order } => { FilterSpec::Bandpass { fl, fu, order } => {
assert!(fl <= fu && fl > 0.);
assert!(order > 0);
let fmid = (fl * fu).sqrt(); let fmid = (fl * fu).sqrt();
let Bw_Hz = fu - fl; let Bw_Hz = fu - fl;
let lp = Self::butter(FilterSpec::Lowpass { fc: fmid, order }); let lp = Self::butter(FilterSpec::Lowpass { fc: fmid, order });
@ -393,7 +431,8 @@ impl ZPKModel {
} }
} }
/// Apply bilinear transform to obtain series biquads from this ZPK model. /// Apply bilinear transform to obtain series biquads from this ZPK model.
/// No prewarping taken into account. /// Pre-warping is taken into account, based on settings stored in
/// [ZPKModel]. Using [ZPKModel::setWarpFreq], this can be overridden.
/// ///
/// # Args /// # Args
/// ///
@ -448,6 +487,87 @@ impl ZPKModel {
} }
SeriesBiquad::newFromBiqs(biqs) SeriesBiquad::newFromBiqs(biqs)
} }
/// Create analog filter prototype for a frequency weighting as used in
/// Sound Level Meters.
///
/// # Args
///
/// - `wt` - `[FreqWeightingType]` 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};
///
/// // Sampling frequency in Hz
/// let fs = 48000.;
///
/// let mut afilter = ZPKModel::freqWeightingFilter(FreqWeightingType::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 {
return ZPKModel {
z: vec![],
p: vec![],
k: 1.0,
fwarp: None,
};
}
let fr: Flt = 1000.;
let fL: Flt = num::Float::powf(10., 1.5);
let fH: Flt = num::Float::powf(10., 3.9);
let sq5: Flt = num::Float::powf(5., 0.5);
let fLsq = fL.powi(2);
let fHsq: Flt = fH.powi(2);
let frsq: Flt = fr.powi(2);
let fA = num::Float::powf(10., 2.45);
let D = num::Float::powf(2., 0.5);
let b = (1. / (1. - D)) * (frsq + fLsq * fHsq / frsq - D * (fLsq + fHsq));
let c = fLsq * fHsq;
let f2 = (3. - sq5) / 2. * fA;
let f3 = (3. + sq5) / 2. * fA;
let f1 = ((-b - (b.powi(2) - 4. * c).sqrt()) / 2.).sqrt();
let f4 = ((-b + (b.powi(2) - 4. * c).sqrt()) / 2.).sqrt();
let p1 = 2. * pi * f1;
let p2 = 2. * pi * f2;
let p3 = 2. * pi * f3;
let p4 = 2. * pi * f4;
let (zeros, poles) = match wt {
FreqWeightingType::Z => {
unreachable!()
}
FreqWeightingType::C => {
let zeros = vec![PoleOrZero::Real2(0., 0.)];
let poles = vec![PoleOrZero::Real2(-p1, -p1), PoleOrZero::Real2(-p4, -p4)];
(zeros, poles)
}
FreqWeightingType::A => {
let poles = vec![
PoleOrZero::Real2(-p1, -p1),
PoleOrZero::Real2(-p2, -p3),
PoleOrZero::Real2(-p4, -p4),
];
let zeros = vec![PoleOrZero::Real2(0., 0.), PoleOrZero::Real2(0., 0.)];
(zeros, poles)
}
};
ZPKModel::new(zeros, poles, 1.0).setGainAt(1000., 1.0)
}
} }
/// Enumeration describing a pole or zero, a complex conjugate pair, a single /// Enumeration describing a pole or zero, a complex conjugate pair, a single

View File

@ -6,11 +6,12 @@
//! //!
//! You will find the following stuff in this crate: //! You will find the following stuff in this crate:
//! //!
//! - Data acquisition, recording, signal generation //! - Data acquisition, recording, signal generation: [daq].
//! - Power spectra estimation, transfer function estimation tools. //! - Power spectra estimation, transfer function estimation tools: [ps].
//! - Sound Level Meter implementation. //! - Sound Level Meter implementation: [slm].
//! - Filter design tools, maybe borrowed from other crates? //! - Filter design tools. I.e. [filter::ZPKModel::butter].
//! - Tools for real time displaying of sensor data //! - Includes bilinear transforms
//! - Tools for real time displaying of sensor data: [rt].
//! //!
//! ## Note to potential users //! ## Note to potential users
//! //!

View File

@ -1,9 +1,3 @@
use crate::config::{pi, Flt};
use crate::filter::{Filter, PoleOrZero, SeriesBiquad, TransferFunction, ZPKModel};
use num::Float;
use std::default;
use std::ops::Deref;
use strum_macros::{Display, EnumMessage}; use strum_macros::{Display, EnumMessage};
/// Sound level frequency weighting type (A, C, Z) /// Sound level frequency weighting type (A, C, Z)
#[derive(Display, Debug, EnumMessage, Default, Clone)] #[derive(Display, Debug, EnumMessage, Default, Clone)]
@ -16,78 +10,6 @@ pub enum FreqWeightingType {
#[default] #[default]
Z, Z,
} }
struct FreqWeightingFilter {
// The calculated filter for this frequency weighting
filter: SeriesBiquad,
}
impl Deref for FreqWeightingFilter {
type Target = SeriesBiquad;
fn deref(&self) -> &Self::Target {
&self.filter
}
}
impl FreqWeightingFilter {
pub fn new(fs: Flt, wt: FreqWeightingType) -> FreqWeightingFilter {
if let FreqWeightingType::Z = wt {
return FreqWeightingFilter {
filter: SeriesBiquad::unit(),
};
}
let fr: Flt = 1000.;
let fL: Flt = Float::powf(10., 1.5);
let fH: Flt = Float::powf(10., 3.9);
let sq5: Flt = Float::powf(5., 0.5);
let fLsq = fL.powi(2);
let fHsq: Flt = fH.powi(2);
let frsq: Flt = fr.powi(2);
let fA = Float::powf(10., 2.45);
let D = Float::powf(2., 0.5);
let b = (1. / (1. - D)) * (frsq + fLsq * fHsq / frsq - D * (fLsq + fHsq));
let c = fLsq * fHsq;
let f2 = (3. - sq5) / 2. * fA;
let f3 = (3. + sq5) / 2. * fA;
let f1 = ((-b - (b.powi(2) - 4. * c).sqrt()) / 2.).sqrt();
let f4 = ((-b + (b.powi(2) - 4. * c).sqrt()) / 2.).sqrt();
let p1 = 2. * pi * f1;
let p2 = 2. * pi * f2;
let p3 = 2. * pi * f3;
let p4 = 2. * pi * f4;
let (zeros, poles) = match wt {
FreqWeightingType::Z => {
unreachable!()
}
FreqWeightingType::C => {
let zeros = vec![PoleOrZero::Real2(0., 0.)];
let poles = vec![PoleOrZero::Real2(-p1, -p1), PoleOrZero::Real2(-p4, -p4)];
(zeros, poles)
}
FreqWeightingType::A => {
let poles = vec![
PoleOrZero::Real2(-p1, -p1),
PoleOrZero::Real2(-p2, -p3),
PoleOrZero::Real2(-p4, -p4),
];
let zeros = vec![PoleOrZero::Real2(0., 0.), PoleOrZero::Real2(0., 0.)];
(zeros, poles)
}
};
return FreqWeightingFilter {
filter: ZPKModel::new(zeros, poles, 1.0)
.setGainAt(1000., 1.0)
.bilinear(fs),
};
}
}
#[cfg(test)] #[cfg(test)]
mod test { mod test {
use super::*; use super::*;