Merge pull request #23 from quartiq/iir-rework

rework iir filter

CHANGELOG.md

@@ -10,6 +10,21 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 * `hbf` FIRs, symmetric FIRs, half band filters, HBF decimators and interpolators
+* `iir::Pid`, `iir:Filter` a builder for PID coefficients and the collection of standard Biquad filters
+* `iir::Biquad::{HOLD, IDENTITY, proportional}`
+* `iir::Biquad` getter/setter
+* `iir`: support for other integers (i8, i16, i128)
+* `iir::Biquad`: support for reduced DF1 state and DF2T state
+* `svf`: state variable filter
+
+### Removed
+
+* `iir::Vec5` type alias has been removed.
+* `iir_int`: integrated into `iir`.
+
+### Changed
+
+* `iir`: The biquad IIR filter API has been reworked. `IIR -> Biquad` renamed.
 
 ## [0.10.0](https://github.com/quartiq/idsp/compare/v0.9.2..v0.10.0) - 2023-07-20
 

Cargo.toml

@@ -16,16 +16,10 @@ num-complex = { version = "0.4.0", features = ["serde"], default-features = fals
 num-traits = { version = "0.2.14", features = ["libm"], default-features = false}
 
 [dev-dependencies]
-easybench = "1.0"
 rand = "0.8"
-ndarray = "0.15"
 rustfft = "6.1.0"
 # futuredsp = "0.0.6"
 # sdr = "0.7.0"
 
-[[bench]]
-name = "micro"
-harness = false
-
 [profile.release]
 debug = 1

README.md

@@ -1,49 +1,92 @@
 # Embedded DSP algorithms
 
 [![GitHub release](https://img.shields.io/github/v/release/quartiq/idsp?include_prereleases)](https://github.com/quartiq/idsp/releases)
+[![crates.io](https://img.shields.io/crates/v/idsp.svg)](https://crates.io/crates/idsp)
 [![Documentation](https://img.shields.io/badge/docs-online-success)](https://docs.rs/idsp)
 [![QUARTIQ Matrix Chat](https://img.shields.io/matrix/quartiq:matrix.org)](https://matrix.to/#/#quartiq:matrix.org)
 [![Continuous Integration](https://github.com/quartiq/idsp/actions/workflows/ci.yml/badge.svg)](https://github.com/quartiq/idsp/actions/workflows/ci.yml)
 
 This crate contains some tuned DSP algorithms for general and especially embedded use.
-Many of the algorithms are implemented on integer datatypes for reasons that become important in certain cases:
-
-* Speed: even with a hard floating point unit integer operations are faster.
-* Accuracy: single precision FP has a 24 bit mantissa, `i32` has full 32 bit.
-* No rounding errors.
-* Natural wrap around (modulo) at the integer overflow: critical for phase/frequency applications.
-* Natural definition of "full scale".
+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`
-
-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.
+## Fixed point
 
-## Two-argument arcus-tangens `atan2`
+### Cosine/Sine
 
-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.
+[`cossin()`] 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.
 
-## ComplexExt
+### Two-argument arcus-tangens
 
-An extension trait for the `num::Complex` type featuring especially a `std`-like API to the two functions above.
+[`atan2()`] 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
 
-High accuracy, zero-assumption, fully robust, forward and reciprocal PLLs with dynamically adjustable time constant and arbitrary (in the Nyquist sampling sense) capture range.
+[`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()`
 
+[`Unwrapper`], [`Accu`], [`saturating_scale()`]:
 Tools to handle, track, and unwrap phase signals or generate them.
 
-## iir_int, iir
+## Float and Fixed point
+
+## IIR/Biquad
+
+[`iir::Biquad`] are 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 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.
+
+### Comparison
+
+This is a rough feature comparison of several available `biquad` crates, with no claim for completeness, accuracy, or even fairness.
+TL;DR: `idsp` is slower but offers more features.
+
+| Feature\Crate | [`biquad-rs`](https://crates.io/crates/biquad) | [`fixed-filters`](https://crates.io/crates/fixed-filters) | `idsp::iir` |
+|---|---|---|---|
+| Floating point `f32`/`f64` | ✅ | ❌ | ✅ |
+| Fixed point `i32` | ❌ | ✅ | ✅ |
+| Parametric fixed point `i32` | ❌ | ✅ | ❌ |
+| Fixed point `i8`/`i16`/`i64`/`i128` | ❌ | ❌ | ✅ |
+| DF2T | ✅ | ❌ | ✅ |
+| Limiting/Clamping | ❌ | ✅ | ✅ |
+| Fixed point accumulator guard bits | ❌ | ❌ | ✅ |
+| Summing junction offset | ❌ | ❌ | ✅ |
+| Fixed point noise shaping | ❌ | ❌ | ✅ |
+| Configuration/state decoupling/multi-channel | ❌ | ❌ | ✅ |
+| `f32` parameter audio filter builder | ✅ | ✅ | ✅ |
+| `f64` parameter audio filter builder | ✅ | ❌ | ✅ |
+| Additional filters (I/HO) | ❌ | ❌ | ✅ |
+| `f32` PI builder | ❌ | ✅ | ✅ |
+| `f32/f64` PI²D² builder | ❌ | ❌ | ✅ |
+| PI²D² builder limits | ❌ | ❌ | ✅ |
+| Support for fixed point `a1=-2` | ❌ | ❌ | ✅ |
+
+Three crates have been compared when processing 4x1M samples (4 channels) with a biquad lowpass.
+Hardware was `thumbv7em-none-eabihf`, `cortex-m7`, code in ITCM, data in DTCM, caches enabled.
+
+| Crate | Type, features | Cycles per sample |
+|---|---|---|
+| [`biquad-rs`](https://crates.io/crates/biquad) | `f32` | 11.4 |
+| `idsp::iir` | `f32`, limits, offset | 15.5 |
+| [`fixed-filters`](https://crates.io/crates/fixed-filters) | `i32`, limits | 20.3 |
+| `idsp::iir` | `i32`, limits, offset | 23.5 |
+| `idsp::iir` | `i32`, limits, offset, noise shaping | 30.0 |
+
+## State variable filter
 
-`i32` and `f32` biquad IIR filters with robust and clean clipping and offset (anti-windup, no derivative kick, dynamically adjustable gains).
+[`svf`] is a 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.
+[`Lowpass`], [`Lockin`] are fast, infinitely cascadable, first- and second-order lowpass and the corresponding integration into a lockin amplifier algorithm.
 
 ## 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.

rustfmt.toml

@@ -0,0 +1 @@
+format_code_in_doc_comments = true

src/accu.rs

@@ -1,12 +1,14 @@
 use num_traits::ops::wrapping::WrappingAdd;
 
+/// Wrapping Accumulator
 #[derive(Copy, Clone, Default, PartialEq, Eq, Debug)]
 pub struct Accu<T> {
     state: T,
     step: T,
 }
 
 impl<T> Accu<T> {
+    /// Create a new accumulator with given initial state and step.
     pub fn new(state: T, step: T) -> Self {
         Self { state, step }
     }

src/complex.rs

@@ -4,11 +4,17 @@ use super::{atan2, cossin};
 
 /// Complex extension trait offering DSP (fast, good accuracy) functionality.
 pub trait ComplexExt<T, U> {
+    /// Unit magnitude from angle
     fn from_angle(angle: T) -> Self;
+    /// Square of magnitude
     fn abs_sqr(&self) -> U;
+    /// Log2 approximation
     fn log2(&self) -> T;
+    /// Angle
     fn arg(&self) -> T;
+    /// Staturating addition
     fn saturating_add(&self, other: Self) -> Self;
+    /// Saturating subtraction
     fn saturating_sub(&self, other: Self) -> Self;
 }
 
@@ -20,8 +26,8 @@ impl ComplexExt<i32, u32> for Complex<i32> {
     /// ```
     /// use idsp::{Complex, ComplexExt};
     /// Complex::<i32>::from_angle(0);
-    /// Complex::<i32>::from_angle(1 << 30);  // pi/2
-    /// Complex::<i32>::from_angle(-1 << 30);  // -pi/2
+    /// Complex::<i32>::from_angle(1 << 30); // pi/2
+    /// Complex::<i32>::from_angle(-1 << 30); // -pi/2
     /// ```
     fn from_angle(angle: i32) -> Self {
         let (c, s) = cossin(angle);
@@ -97,6 +103,7 @@ impl ComplexExt<i32, u32> for Complex<i32> {
 
 /// Full scale fixed point multiplication.
 pub trait MulScaled<T> {
+    /// Scaled multiplication for fixed point
     fn mul_scaled(self, other: T) -> Self;
 }
 

src/filter.rs

@@ -1,4 +1,9 @@
+/// Single inpout single output i32 filter
 pub trait Filter {
+    /// Filter configuration type.
+    ///
+    /// While the filter struct owns the state,
+    /// the configuration is decoupled to allow sharing.
     type Config;
     /// Update the filter with a new sample.
     ///
@@ -16,6 +21,9 @@ pub trait Filter {
     fn set(&mut self, x: i32);
 }
 
+/// Nyquist zero
+///
+/// Filter with a flat transfer function and a transfer function zero at Nyquist.
 #[derive(Copy, Clone, Default)]
 pub struct Nyquist(i32);
 impl Filter for Nyquist {
@@ -34,6 +42,7 @@ impl Filter for Nyquist {
     }
 }
 
+/// Repeat another filter
 #[derive(Copy, Clone)]
 pub struct Repeat<const N: usize, T>([T; N]);
 impl<const N: usize, T: Filter> Filter for Repeat<N, T> {
@@ -54,6 +63,7 @@ impl<const N: usize, T: Default + Copy> Default for Repeat<N, T> {
     }
 }
 
+/// Combine two different filters in cascade
 #[derive(Copy, Clone, Default)]
 pub struct Cascade<T, U>(T, U);
 impl<T: Filter, U: Filter> Filter for Cascade<T, U> {

src/hbf.rs

@@ -1,3 +1,7 @@
+//! Half-band filters and cascades
+//!
+//! Used to perform very efficient high-dynamic range rate changes by powers of two.
+
 use core::{
     iter::Sum,
     ops::{Add, Mul},
@@ -457,12 +461,18 @@ impl Default for HbfDecCascade {
 }
 
 impl HbfDecCascade {
+    /// Set cascade depth
+    ///
+    /// Sets the number of HBF filter stages to apply.
     #[inline]
     pub fn set_depth(&mut self, n: usize) {
         assert!(n <= 4);
         self.depth = n;
     }
 
+    /// Cascade depth
+    ///
+    /// The number of HBF filter stages to apply.
     #[inline]
     pub fn depth(&self) -> usize {
         self.depth
@@ -543,7 +553,7 @@ impl Filter for HbfDecCascade {
 #[derive(Copy, Clone, Debug)]
 pub struct HbfIntCascade {
     depth: usize,
-    pub stages: (
+    stages: (
         HbfInt<
             'static,
             f32,
@@ -586,11 +596,17 @@ impl Default for HbfIntCascade {
 }
 
 impl HbfIntCascade {
+    /// Set cascade depth
+    ///
+    /// Sets the number of HBF filter stages to apply.
     pub fn set_depth(&mut self, n: usize) {
         assert!(n <= 4);
         self.depth = n;
     }
 
+    /// Cascade depth
+    ///
+    /// The number of HBF filter stages to apply.
     pub fn depth(&self) -> usize {
         self.depth
     }