From 6787d0a7f17f2879be184e157b0c45471dcff9f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=B6rdens?= Date: Fri, 12 Jan 2024 17:46:27 +0100 Subject: [PATCH] docs --- README.md | 45 ++++++++++++---- src/iir/biquad.rs | 128 +++++++++++++++++++++++++++++----------------- src/iir/pid.rs | 4 +- src/num.rs | 6 +-- 4 files changed, 121 insertions(+), 62 deletions(-) diff --git a/README.md b/README.md index fc61986..1988c93 100644 --- a/README.md +++ b/README.md @@ -11,38 +11,63 @@ Many of the algorithms are implemented on integer (fixed point) datatypes. One comprehensive user for these algorithms is [Stabilizer](https://github.com/quartiq/stabilizer). -## Cosine/Sine `cossin` +## Fixed point + +### Cosine/Sine [`cossin()`] This uses a small (128 element or 512 byte) LUT, smart octant (un)mapping, linear interpolation and comprehensive analysis of corner cases to achieve a very clean signal (4e-6 RMS error, 9e-6 max error, 108 dB SNR typ), low spurs, and no bias with about 40 cortex-m instruction per call. It computes both cosine and sine (i.e. the complex signal) at once given a phase input. -## Two-argument arcus-tangens `atan2` +### Two-argument arcus-tangens [`atan2()`] This returns a phase given a complex signal (a pair of in-phase/`x`/cosine and quadrature/`y`/sine). The RMS phase error is less than 5e-6 rad, max error is less than 1.2e-5 rad, i.e. 20.5 bit RMS, 19.1 bit max accuracy. The bias is minimal. -## PLL, RPLL +## [`PLL`], [`RPLL`] High accuracy, zero-assumption, fully robust, forward and reciprocal PLLs with dynamically adjustable time constant and arbitrary (in the Nyquist sampling sense) capture range, and noise shaping. -## Unwrapper, Accu, saturating_scale +## [`Unwrapper`], [`Accu`], [`saturating_scale()`] Tools to handle, track, and unwrap phase signals or generate them. -## IIR/Biquad +## Float and Fixed point + +## [`iir`]/[`iir::Biquad`] Fixed point (`i8`, `i16`, `i32`, `i64`) and floating point (`f32`, `f64`) biquad IIR filters. -Robust and clean clipping and offset (anti-windup, no derivative kick, dynamically adjustable gains) suitable for PID controller applications. -Direct Form 1 and Direct Form 2 Transposed supported. +Robust and clean clipping and offset (anti-windup, no derivative kick, dynamically adjustable gains and gain limits) suitable for PID controller applications. +Three kinds of filter actions: Direct Form 1, Direct Form 2 Transposed, and Direct Form 1 with noise shaping supported. Coefficient sharing for multiple channels. -## State variable filter +Compared to [`biquad-rs`](https://crates.io/crates/biquad) this crates adds several additional important features: + +* fixed point implementations (`i32`, `i64`, etc.) in addition to `f32`/`f64` floating point +* additional [`iir::Filter`] builders (I/HO) +* decoupled [`iir::Biquad`] configuration and flat `[T; N]` state +* [`iir::Pid`] builder +* DF1 noise shaping for fixed point +* proper output limiting/clamping before feedback ("anti-windup") +* summing junction offset (for PID controller applications) + +Compared to [`fixed-filters`](https://crates.io/crates/fixed-filters) this crate: + +* Supports unified floating point and fixed point API +* decoupled [`iir::Biquad`] configuration and flat `[T; N]` state +* [`iir::Pid`] builder +* additional [`iir::Filter`] builders (I/HO) +* noise shaping for fixed point +* summing junction offset (for PID controller applications) + +## [`svf`] State variable filter Simple IIR state variable filter simultaneously providing highpass, lowpass, bandpass, and notch filtering of a signal. -## Lowpass, Lockin +## [`Lowpass`], [`Lockin`] Fast, infinitely cascadable, first- and second-order lowpass and the corresponding integration into a lockin amplifier algorithm. -## FIR filters +## FIR filters: [`hbf::HbfDec`], [`hbf::HbfInt`], [`hbf::HbfDecCascade`], [`hbf::HbfIntCascade`] Fast `f32` symmetric FIR filters, optimized half-band filters, half-band filter decimators and integators and cascades. +These are used in [`stabilizer-stream`](https://github.com/quartiq/stabilizer-stream) for online PSD calculation on log +frequency scale for arbitrarily large amounts of data. diff --git a/src/iir/biquad.rs b/src/iir/biquad.rs index 4b64447..95de077 100644 --- a/src/iir/biquad.rs +++ b/src/iir/biquad.rs @@ -1,51 +1,85 @@ use num_traits::{AsPrimitive, Float}; use serde::{Deserialize, Serialize}; -use crate::FilterNum; +use crate::Coefficient; -/// Filter architecture +/// Biquad IIR filter +/// +/// A biquadratic IIR filter supports up to two zeros and two poles in the transfer function. +/// It can be used to implement a wide range of responses to input signals. +/// +/// The Biquad performs the following operation to compute a new output sample `y0` from a new +/// input sample `x0` given its configuration and previous samples: +/// +/// `y0 = clamp(b0*x0 + b1*x1 + b2*x2 - a1*y1 - a2*y2 + u, min, max)` +/// +/// This implementation here saves storage and improves caching opportunities by decoupling +/// filter configuration (coefficients, limits and offset) from filter state +/// and thus supports both (a) sharing a single filter between multiple states ("channels") and (b) +/// rapid switching of filters (tuning, transfer) for a given state without copying either +/// state of configuration. +/// +/// # Filter architecture /// /// Direct Form 1 (DF1) and Direct Form 2 transposed (DF2T) are the only IIR filter /// structures with an (effective bin the case of TDF2) single summing junction /// this allows clamping of the output before feedback. /// -/// DF1 allows atomic coefficient change because only x/y are pipelined. +/// DF1 allows atomic coefficient change because only inputs and outputs are pipelined. /// The summing junctuion pipelining of TDF2 would require incremental /// coefficient changes and is thus less amenable to online tuning. /// -/// DF2T needs less state storage (2 instead of 4). This is in addition to the coefficient storage -/// (5 plus 2 limits plus 1 offset) -/// This implementation already saves storage by decoupling coefficients/limits and offset from state -/// and thus supports both (a) sharing a single filter between multiple states ("channels") and (b) -/// rapid switching of filters (tuning, transfer) for a given state without copying. +/// DF2T needs less state storage (2 instead of 4). This is in addition to the coefficient +/// storage (5 plus 2 limits plus 1 offset) /// /// DF2T is less efficient and accurate for fixed-point architectures as quantization /// happens at each intermediate summing junction in addition to the output quantization. This is /// especially true for common `i64 + i32 * i32 -> i64` MACC architectures. +/// One could use wide state storage for fixed point DF2T but that would negate the storage +/// and processing advantages. +/// +/// # Coefficients +/// +/// `ba: [T; 5] = [b0, b1, b2, a1, a2]` is the coefficients type. +/// To represent the IIR coefficients, this contains the feed-forward +/// coefficients `b0, b1, b2` followed by the feed-back coefficients +/// `a1, a2`, all five normalized such that `a0 = 1`. +/// +/// The summing junction of the filter also receives an offset `u`. +/// +/// The filter applies clamping such that `min <= y <= max`. +/// +/// See [`crate::iir::Filter`] and [`crate::iir::Pid`] for ways to generate coefficients. +/// +/// # Fixed point /// -/// # Coefficients and state +/// Coefficient scaling (see [`Coefficient`]) is fixed and optimized such that -2 is exactly +/// representable. This is tailored to low-passes, PID, II etc, where the integration rule is +/// [1, -2, 1]. /// -/// `[T; 5]` is the coefficients type. +/// There are two guard bits in the accumulator before clamping/limiting. +/// While this isn't enough to cover the worst case accumulator, it does catch many real world +/// overflow cases. +/// +/// # State /// /// To represent the IIR state (input and output memory) during [`Biquad::update()`] -/// this contains the two previous inputs and output `[x1, x2, y1, y2]` +/// the DF1 state contains the two previous inputs and output `[x1, x2, y1, y2]` /// concatenated. Lower indices correspond to more recent samples. -/// To represent the IIR coefficients, this contains the feed-forward -/// coefficients `[b0, b1, b2]` followd by the negated feed-back coefficients -/// `[a1, a2]`, all five normalized such that `a0 = 1`. -/// Note that between filter [`Biquad::update()`] the `xy` state contains -/// `[x0, x1, y0, y1]`. /// -/// The IIR coefficients can be mapped to other transfer function -/// representations, for example as described in +/// In the DF2T case the state contains `[b1*x1 + b2*x2 - a1*y1 - a2*y2, b2*x1 - a2*y1]` /// -/// # IIR filter as PID controller +/// In the DF1 case with first order noise shaping, the state contains `[x1, x2, y1, y2, e1]` +/// where `e0` is the accumulated quantization error. /// -/// Contains the coeeficients `ba`, the summing junction offset `u`, and the -/// output limits `min` and `max`. Data is represented in floating-point -/// for all internal signals, input and output. +/// # PID controller /// -/// This implementation achieves several important properties: +/// The IIR coefficients can be mapped to other transfer function +/// representations, for example PID controllers as described in +/// and +/// . +/// +/// Using a Biquad as a template for a PID controller achieves several important properties: /// /// * Its transfer function is universal in the sense that any biquadratic /// transfer function can be implemented (high-passes, gain limits, second @@ -66,16 +100,6 @@ use crate::FilterNum; /// coefficients/offset sets. /// * Cascading multiple IIR filters allows stable and robust /// implementation of transfer functions beyond bequadratic terms. -/// -/// See also . -/// -/// -/// Offset and limiting disabled to suit lowpass applications. -/// Coefficient scaling fixed and optimized such that -2 is representable. -/// Tailored to low-passes, PID, II etc, where the integration rule is [1, -2, 1]. -/// Since the relevant coefficients `a1` and `a2` are negated, we also negate the -/// stored `y1` and `y2` in the state. -/// Note that `xy` contains the negative `y1` and `y2`, such that `-a1` #[derive(Copy, Clone, Debug, Deserialize, Serialize, PartialEq, PartialOrd)] pub struct Biquad { ba: [T; 5], @@ -84,7 +108,7 @@ pub struct Biquad { max: T, } -impl Default for Biquad { +impl Default for Biquad { fn default() -> Self { Self { ba: [T::ZERO; 5], @@ -95,7 +119,7 @@ impl Default for Biquad { } } -impl From<[T; 5]> for Biquad { +impl From<[T; 5]> for Biquad { fn from(ba: [T; 5]) -> Self { Self { ba, @@ -106,7 +130,7 @@ impl From<[T; 5]> for Biquad { impl From<&[C; 6]> for Biquad where - T: FilterNum + AsPrimitive, + T: Coefficient + AsPrimitive, C: Float + AsPrimitive, { fn from(ba: &[C; 6]) -> Self { @@ -124,7 +148,7 @@ where impl From<&Biquad> for [C; 6] where - T: FilterNum + AsPrimitive, + T: Coefficient + AsPrimitive, C: 'static + Copy, { fn from(value: &Biquad) -> Self { @@ -140,7 +164,7 @@ where } } -impl Biquad { +impl Biquad { /// A "hold" filter that ingests input and maintains output /// /// ``` @@ -193,7 +217,7 @@ impl Biquad { /// `y0 = clamp(b0*x0 + b1*x1 + b2*x2 - a1*y1 - a2*y2 + u, min, max)`. /// /// ``` - /// # use idsp::FilterNum; + /// # use idsp::Coefficient; /// # use idsp::iir::*; /// assert_eq!(Biquad::::IDENTITY.ba()[0], i32::ONE); /// assert_eq!(Biquad::::HOLD.ba()[3], -i32::ONE); @@ -207,7 +231,7 @@ impl Biquad { /// See [`Biquad::ba()`]. /// /// ``` - /// # use idsp::FilterNum; + /// # use idsp::Coefficient; /// # use idsp::iir::*; /// let mut i = Biquad::default(); /// i.ba_mut()[0] = i32::ONE; @@ -321,7 +345,7 @@ impl Biquad { /// Compute input-referred (`x`) offset. /// /// ``` - /// # use idsp::FilterNum; + /// # use idsp::Coefficient; /// # use idsp::iir::*; /// let mut i = Biquad::proportional(3); /// i.set_u(3); @@ -349,7 +373,7 @@ impl Biquad { /// ``` /// /// ``` - /// # use idsp::FilterNum; + /// # use idsp::Coefficient; /// # use idsp::iir::*; /// let mut i = Biquad::proportional(-i32::ONE); /// i.set_input_offset(1); @@ -384,18 +408,28 @@ impl Biquad { /// /// ## `N=5` Direct Form 1 with first order noise shaping /// + /// ``` + /// # use idsp::iir::*; + /// let mut xy = [1, 2, 3, 4, 5]; + /// let x0 = 6; + /// let y0 = Biquad::IDENTITY.update(&mut xy, x0); + /// assert_eq!(y0, x0); + /// assert_eq!(xy, [x0, 1, y0, 3, 5]); + /// ``` + /// /// `xy` contains: /// * On entry: `[x1, x2, y1, y2, e1]` /// * On exit: `[x0, x1, y0, y1, e0]` /// - /// Note: This isn't useful for floating point. + /// Note: This is only useful for fixed point filters. /// /// ## `N=2` Direct Form 2 transposed /// - /// Note: Don't use this for fixed point. Quantization happens at each state store operation. - /// Ideally the state would be `T::ACCU` but then for fixed point it would use equal amount - /// of storage compared to DF1 for little to no gain in performance and none in functionality. - /// There are also no guard bits. + /// Note: This is only useful for floating point filters. + /// Don't use this for fixed point: Quantization happens at each state store operation. + /// Ideally the state would be `[T::ACCU; 2]` but then for fixed point it would use equal amount + /// of storage compared to DF1 for no gain in performance and loss in functionality. + /// There are also no guard bits here. /// /// `xy` contains: /// * On entry: `[b1*x1 + b2*x2 - a1*y1 - a2*y2, b2*x1 - a2*y1]` diff --git a/src/iir/pid.rs b/src/iir/pid.rs index c0159ce..388ec09 100644 --- a/src/iir/pid.rs +++ b/src/iir/pid.rs @@ -1,7 +1,7 @@ use num_traits::{AsPrimitive, Float}; use serde::{Deserialize, Serialize}; -use crate::FilterNum; +use crate::Coefficient; /// PID controller builder /// @@ -156,7 +156,7 @@ impl Pid { /// /// # Panic /// Will panic in debug mode on fixed point coefficient overflow. - pub fn build>(&self) -> Result<[C; 5], PidError> + pub fn build>(&self) -> Result<[C; 5], PidError> where T: AsPrimitive, { diff --git a/src/num.rs b/src/num.rs index 0c94963..4ab24d5 100644 --- a/src/num.rs +++ b/src/num.rs @@ -1,7 +1,7 @@ use num_traits::{AsPrimitive, Float, Num}; /// Helper trait unifying fixed point and floating point coefficients/samples -pub trait FilterNum: 'static + Copy + Num + AsPrimitive { +pub trait Coefficient: 'static + Copy + Num + AsPrimitive { /// Multiplicative identity const ONE: Self; /// Negative multiplicative identity, equal to `-Self::ONE`. @@ -41,7 +41,7 @@ pub trait FilterNum: 'static + Copy + Num + AsPrimitive { macro_rules! impl_float { ($T:ty) => { - impl FilterNum for $T { + impl Coefficient for $T { const ONE: Self = 1.0; const NEG_ONE: Self = -1.0; const ZERO: Self = 0.0; @@ -82,7 +82,7 @@ impl_float!(f64); macro_rules! impl_int { ($T:ty, $U:ty, $A:ty, $Q:literal) => { - impl FilterNum for $T { + impl Coefficient for $T { const ONE: Self = 1 << $Q; const NEG_ONE: Self = -1 << $Q; const ZERO: Self = 0;