docs

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<T>`] 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<T>`] 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.

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 <https://arxiv.org/abs/1508.06319>
+/// 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
+/// <https://hackmd.io/IACbwcOTSt6Adj3_F9bKuw> and
+/// <https://arxiv.org/abs/1508.06319>.
+///
+/// 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 <https://hackmd.io/IACbwcOTSt6Adj3_F9bKuw>.
-///
-///
-/// 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<T> {
     ba: [T; 5],
@@ -84,7 +108,7 @@ pub struct Biquad<T> {
     max: T,
 }
 
-impl<T: FilterNum> Default for Biquad<T> {
+impl<T: Coefficient> Default for Biquad<T> {
     fn default() -> Self {
         Self {
             ba: [T::ZERO; 5],
@@ -95,7 +119,7 @@ impl<T: FilterNum> Default for Biquad<T> {
     }
 }
 
-impl<T: FilterNum> From<[T; 5]> for Biquad<T> {
+impl<T: Coefficient> From<[T; 5]> for Biquad<T> {
     fn from(ba: [T; 5]) -> Self {
         Self {
             ba,
@@ -106,7 +130,7 @@ impl<T: FilterNum> From<[T; 5]> for Biquad<T> {
 
 impl<T, C> From<&[C; 6]> for Biquad<T>
 where
-    T: FilterNum + AsPrimitive<C>,
+    T: Coefficient + AsPrimitive<C>,
     C: Float + AsPrimitive<T>,
 {
     fn from(ba: &[C; 6]) -> Self {
@@ -124,7 +148,7 @@ where
 
 impl<T, C> From<&Biquad<T>> for [C; 6]
 where
-    T: FilterNum + AsPrimitive<C>,
+    T: Coefficient + AsPrimitive<C>,
     C: 'static + Copy,
 {
     fn from(value: &Biquad<T>) -> Self {
@@ -140,7 +164,7 @@ where
     }
 }
 
-impl<T: FilterNum> Biquad<T> {
+impl<T: Coefficient> Biquad<T> {
     /// A "hold" filter that ingests input and maintains output
     ///
     /// ```
@@ -193,7 +217,7 @@ impl<T: FilterNum> Biquad<T> {
     /// `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::<i32>::IDENTITY.ba()[0], i32::ONE);
     /// assert_eq!(Biquad::<i32>::HOLD.ba()[3], -i32::ONE);
@@ -207,7 +231,7 @@ impl<T: FilterNum> Biquad<T> {
     /// 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<T: FilterNum> Biquad<T> {
     /// 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<T: FilterNum> Biquad<T> {
     /// ```
     ///
     /// ```
-    /// # 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<T: FilterNum> Biquad<T> {
     ///
     /// ## `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]`

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<T: Float> Pid<T> {
     ///
     /// # Panic
     /// Will panic in debug mode on fixed point coefficient overflow.
-    pub fn build<C: FilterNum + AsPrimitive<T>>(&self) -> Result<[C; 5], PidError>
+    pub fn build<C: Coefficient + AsPrimitive<T>>(&self) -> Result<[C; 5], PidError>
     where
         T: AsPrimitive<C>,
     {

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<Self::ACCU> {
+pub trait Coefficient: 'static + Copy + Num + AsPrimitive<Self::ACCU> {
     /// Multiplicative identity
     const ONE: Self;
     /// Negative multiplicative identity, equal to `-Self::ONE`.
@@ -41,7 +41,7 @@ pub trait FilterNum: 'static + Copy + Num + AsPrimitive<Self::ACCU> {
 
 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;