diff --git a/src/filter/mod.rs b/src/filter/mod.rs index 666aea4..6ca5c7e 100644 --- a/src/filter/mod.rs +++ b/src/filter/mod.rs @@ -13,6 +13,7 @@ mod seriesbiquad; mod zpkmodel; mod butter; +pub use super::ps::FreqWeightingType; pub use biquad::Biquad; pub use biquadbank::BiquadBank; pub use dummy::DummyFilter; diff --git a/src/filter/zpkmodel.rs b/src/filter/zpkmodel.rs index 84569ad..29e29fc 100644 --- a/src/filter/zpkmodel.rs +++ b/src/filter/zpkmodel.rs @@ -1,13 +1,12 @@ -use std::cmp::{max, min}; - use super::butter::butter_lowpass_roots; +use super::{Biquad, SeriesBiquad, TransferFunction}; +use crate::{config::*, ps::FreqWeightingType}; use itertools::{EitherOrBoth, Itertools}; use num::{zero, Complex}; +use std::cmp::{max, min}; -use crate::config::*; - -use super::{Biquad, SeriesBiquad, TransferFunction}; - +/// Reasonable maximum order of Butterworth filters. +pub const BUTTER_MAX_ORDER: u32 = 40; /// Specification of a filter for a certain type. /// /// The order corresponds to the rolloff in dB/decade. order=1 means 20 @@ -47,12 +46,17 @@ pub enum FilterSpec { /// Analog zero-pole-gain model for real input to real output. Provides methods /// to generate analog filters of standard type, i.e. bandpass, lowpass and /// highpass. These can subsequentially be used to generate a digital filter. -/// -/// # Example -/// +/// +/// +/// # Example: Create a digital second order Butterworth bandpass filter +/// /// ```rust /// 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); +/// /// ``` /// /// It has a transfer function that can be described as a rational function of @@ -63,7 +67,7 @@ pub enum FilterSpec { /// H(s) = k ------------ /// Π_i (s-p_i) /// ``` -/// +/// /// where `Π` denotes the product of a series, `z_i` are the zeros, and `p_i` /// are the poles. In order to have real output for a real input, the zeros and /// poles should either be real, or come in complex conjugate pairs. This is @@ -230,8 +234,12 @@ impl ZPKModel { ZPKModel { z, p, k, fwarp } } - // Set critical frequency in filter - fn setWarpFreq(mut self, fcrit: Flt) -> ZPKModel { + /// Set critical frequency in filter, used for bilinear transform. + /// + /// # Args + /// + /// - `fcrit` - New critical frequency in \[Hz\]. + pub fn setWarpFreq(mut self, fcrit: Flt) -> ZPKModel { self.fwarp = Some(fcrit); self } @@ -346,12 +354,46 @@ impl ZPKModel { } } - /// Create a Butterworth filter according to a certain specification - pub fn butter(spec: FilterSpec) -> ZPKModel { + fn check_spec(spec: FilterSpec) { + 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 { FilterSpec::Lowpass { fc, order } => { - assert!(fc > 0.); - assert!(order > 0); let p = butter_lowpass_roots(fc, order as u32).collect(); let z = vec![]; ZPKModel { @@ -365,8 +407,6 @@ impl ZPKModel { .setWarpFreq(fc) } FilterSpec::Highpass { fc, order } => { - assert!(fc > 0.); - assert!(order > 0); let p = butter_lowpass_roots(fc, order as u32).collect(); let z = vec![PoleOrZero::Real1(0.); order as usize]; ZPKModel { @@ -380,8 +420,6 @@ impl ZPKModel { .setWarpFreq(fc) } FilterSpec::Bandpass { fl, fu, order } => { - assert!(fl <= fu && fl > 0.); - assert!(order > 0); let fmid = (fl * fu).sqrt(); let Bw_Hz = fu - fl; 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. - /// 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 /// @@ -448,6 +487,87 @@ impl ZPKModel { } 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 diff --git a/src/lib.rs b/src/lib.rs index 949456a..b7f34a3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -6,11 +6,12 @@ //! //! You will find the following stuff in this crate: //! -//! - Data acquisition, recording, signal generation -//! - Power spectra estimation, transfer function estimation tools. -//! - Sound Level Meter implementation. -//! - Filter design tools, maybe borrowed from other crates? -//! - Tools for real time displaying of sensor data +//! - Data acquisition, recording, signal generation: [daq]. +//! - Power spectra estimation, transfer function estimation tools: [ps]. +//! - Sound Level Meter implementation: [slm]. +//! - Filter design tools. I.e. [filter::ZPKModel::butter]. +//! - Includes bilinear transforms +//! - Tools for real time displaying of sensor data: [rt]. //! //! ## Note to potential users //! diff --git a/src/ps/freqweighting.rs b/src/ps/freqweighting.rs index 7bdefe8..156b160 100644 --- a/src/ps/freqweighting.rs +++ b/src/ps/freqweighting.rs @@ -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}; /// Sound level frequency weighting type (A, C, Z) #[derive(Display, Debug, EnumMessage, Default, Clone)] @@ -16,78 +10,6 @@ pub enum FreqWeightingType { #[default] 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)] mod test { use super::*;