Merge branch 'master' into feature/qspi-stream

.github/workflows/ci.yml

@@ -77,6 +77,25 @@ jobs:
           command: build
           args: --release --features semihosting
 
+  test:
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        toolchain:
+          - stable
+          - beta
+    steps:
+      - uses: actions/checkout@v2
+      - name: Install Rust ${{ matrix.toolchain }}
+        uses: actions-rs/toolchain@v1
+        with:
+          toolchain: ${{ matrix.toolchain }}
+      - name: cargo test
+        uses: actions-rs/cargo@v1
+        with:
+          command: test
+          args: --package dsp --target=x86_64-unknown-linux-gnu
+
   # Tell bors about it
   # https://github.com/rtic-rs/cortex-m-rtic/blob/8a4f9c6b8ae91bebeea0791680f89375a78bffc6/.github/workflows/build.yml#L566-L603
   ci-success:

.gitignore

@@ -1,2 +1,3 @@
 /target
+/dsp/target
 .gdb_history

Cargo.lock

@@ -178,9 +178,9 @@ dependencies = [
 
 [[package]]
 name = "cortex-m-semihosting"
-version = "0.3.5"
+version = "0.3.7"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "113ef0ecffee2b62b58f9380f4469099b30e9f9cbee2804771b4203ba1762cfa"
+checksum = "6bffa6c1454368a6aa4811ae60964c38e6996d397ff8095a8b9211b1c1f749bc"
 dependencies = [
  "cortex-m",
 ]
@@ -189,6 +189,7 @@ dependencies = [
 name = "dsp"
 version = "0.1.0"
 dependencies = [
+ "libm",
  "serde",
 ]
 
@@ -297,6 +298,12 @@ dependencies = [
  "hashbrown",
 ]
 
+[[package]]
+name = "libm"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c7d73b3f436185384286bd8098d17ec07c9a7d2388a6599f824d8502b529702a"
+
 [[package]]
 name = "log"
 version = "0.4.11"
@@ -383,9 +390,9 @@ checksum = "e2a38df5b15c8d5c7e8654189744d8e396bddc18ad48041a500ce52d6948941f"
 
 [[package]]
 name = "rtic-core"
-version = "0.3.0"
+version = "0.3.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "ab51fe832317e805f869b3d859f91aadf855c2c3da51f9b84bc645c201597158"
+checksum = "8bd58a6949de8ff797a346a28d9f13f7b8f54fa61bb5e3cb0985a4efb497a5ef"
 
 [[package]]
 name = "rtic-syntax"
@@ -424,9 +431,9 @@ checksum = "388a1df253eca08550bef6c72392cfe7c30914bf41df5269b68cbd6ff8f570a3"
 
 [[package]]
 name = "serde"
-version = "1.0.117"
+version = "1.0.118"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "b88fa983de7720629c9387e9f517353ed404164b1e482c970a90c1a4aaf7dc1a"
+checksum = "06c64263859d87aa2eb554587e2d23183398d617427327cf2b3d0ed8c69e4800"
 dependencies = [
  "serde_derive",
 ]
@@ -443,9 +450,9 @@ dependencies = [
 
 [[package]]
 name = "serde_derive"
-version = "1.0.117"
+version = "1.0.118"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "cbd1ae72adb44aab48f325a02444a5fc079349a8d804c1fc922aed3f7454c74e"
+checksum = "c84d3526699cd55261af4b941e4e725444df67aa4f9e6a3564f18030d12672df"
 dependencies = [
  "proc-macro2",
  "quote",
@@ -527,9 +534,9 @@ dependencies = [
 
 [[package]]
 name = "syn"
-version = "1.0.48"
+version = "1.0.53"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "cc371affeffc477f42a221a1e4297aedcea33d47d19b61455588bd9d8f6b19ac"
+checksum = "8833e20724c24de12bbaba5ad230ea61c3eafb05b881c7c9d3cfe8638b187e68"
 dependencies = [
  "proc-macro2",
  "quote",

Cargo.toml

@@ -18,7 +18,6 @@ exclude = [
 ]
 
 [badges]
-travis-ci = { repository = "quartiq/stabilizer", branch = "master" }
 maintenance = { status = "experimental" }
 
 [package.metadata.docs.rs]
@@ -42,6 +41,7 @@ asm-delay = "0.9.0"
 enum-iterator = "0.6.0"
 paste = "1"
 dsp = { path = "dsp" }
+ad9959 = { path = "ad9959" }
 
 [dependencies.mcp23017]
 git = "https://github.com/mrd0ll4r/mcp23017.git"
@@ -51,9 +51,6 @@ version = "0.6"
 features = ["ethernet", "proto-ipv4", "socket-tcp", "proto-ipv6"]
 default-features = false
 
-[dependencies.ad9959]
-path = "ad9959"
-
 [dependencies.stm32h7xx-hal]
 features = ["stm32h743v", "rt", "unproven", "ethernet", "quadspi"]
 git = "https://github.com/stm32-rs/stm32h7xx-hal"

dsp/Cargo.lock

@@ -4,9 +4,16 @@
 name = "dsp"
 version = "0.1.0"
 dependencies = [
+ "libm",
  "serde",
 ]
 
+[[package]]
+name = "libm"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c7d73b3f436185384286bd8098d17ec07c9a7d2388a6599f824d8502b529702a"
+
 [[package]]
 name = "proc-macro2"
 version = "1.0.24"

dsp/Cargo.toml

@@ -5,6 +5,7 @@ authors = ["Robert Jördens <rj@quartiq.de>"]
 edition = "2018"
 
 [dependencies]
+libm = "0.2.1"
 serde = { version = "1.0", features = ["derive"], default-features = false }
 
 [features]

dsp/src/lib.rs

@@ -1,4 +1,11 @@
-#![no_std]
+#![cfg_attr(not(test), no_std)]
 #![cfg_attr(feature = "nightly", feature(asm, core_intrinsics))]
 
+pub type Complex<T> = (T, T);
 pub mod iir;
+pub mod lockin;
+pub mod pll;
+pub mod unwrap;
+
+#[cfg(test)]
+mod testing;

dsp/src/lockin.rs

@@ -0,0 +1,518 @@
+//! Lock-in amplifier.
+//!
+//! Lock-in processing is performed through a combination of the
+//! following modular processing blocks: demodulation, filtering,
+//! decimation and computing the magnitude and phase from a complex
+//! signal. These processing blocks are mutually independent.
+//!
+//! # Terminology
+//!
+//! * _demodulation signal_ - A copy of the reference signal that is
+//! optionally frequency scaled and phase shifted. This is a complex
+//! signal. The demodulation signals are used to demodulate the ADC
+//! sampled signal.
+//! * _internal clock_ - A fast internal clock used to increment a
+//! counter for determining the 0-phase points of a reference signal.
+//! * _reference signal_ - A constant-frequency signal used to derive
+//! the demodulation signal.
+//! * _timestamp_ - Timestamps record the timing of the reference
+//! signal's 0-phase points. For instance, if a reference signal is
+//! provided externally, a fast internal clock increments a
+//! counter. When the external reference reaches the 0-phase point
+//! (e.g., a positive edge), the value of the counter is recorded as a
+//! timestamp. These timestamps are used to determine the frequency
+//! and phase of the reference signal.
+//!
+//! # Usage
+//!
+//! The first step is to initialize a `Lockin` instance with
+//! `Lockin::new()`. This provides the lock-in algorithms with
+//! necessary information about the demodulation and filtering steps,
+//! such as whether to demodulate with a harmonic of the reference
+//! signal and the IIR biquad filter to use. There are then 4
+//! different processing steps that can be used:
+//!
+//! * `demodulate` - Computes the phase of the demodulation signal
+//! corresponding to each ADC sample, uses this phase to compute the
+//! demodulation signal, and multiplies this demodulation signal by
+//! the ADC-sampled signal. This is a method of `Lockin` since it
+//! requires information about how to modify the reference signal for
+//! demodulation.
+//! * `filter` - Performs IIR biquad filtering of a complex
+//! signals. This is commonly performed on the signal provided by the
+//! demodulation step, but can be performed at any other point in the
+//! processing chain or omitted entirely. `filter` is a method of
+//! `Lockin` since it must hold onto the filter configuration and
+//! state.
+//! * `decimate` - This decimates a signal to reduce the load on the
+//! DAC output. It does not require any state information and is
+//! therefore a normal function.
+//! * `magnitude_phase` - Computes the magnitude and phase of the
+//! component of the ADC-sampled signal whose frequency is equal to
+//! the demodulation frequency. This does not require any state
+//! information and is therefore a normal function.
+
+use super::iir::{IIRState, IIR};
+use super::Complex;
+use core::f32::consts::PI;
+
+/// The number of ADC samples in one batch.
+pub const ADC_SAMPLE_BUFFER_SIZE: usize = 16;
+/// The number of outputs sent to the DAC for each ADC batch.
+pub const DECIMATED_BUFFER_SIZE: usize = 1;
+
+/// Treat the 2-element array as a FIFO. This allows new elements to
+/// be pushed into the array, existing elements to shift back in the
+/// array, and the last element to fall off the array.
+trait Fifo2<T> {
+    fn push(&mut self, new_element: Option<T>);
+}
+
+impl<T: Copy> Fifo2<T> for [Option<T>; 2] {
+    /// Push a new element into the array. The existing elements move
+    /// backward in the array by one location, and the current last
+    /// element is discarded.
+    ///
+    /// # Arguments
+    ///
+    /// * `new_element` - New element pushed into the front of the
+    /// array.
+    fn push(&mut self, new_element: Option<T>) {
+        // For array sizes greater than 2 it would be preferable to
+        // use a rotating index to avoid unnecessary data
+        // copying. However, this would somewhat complicate the use of
+        // iterators and for 2 elements, shifting is inexpensive.
+        self[1] = self[0];
+        self[0] = new_element;
+    }
+}
+
+/// Performs lock-in amplifier processing of a signal.
+pub struct Lockin {
+    phase_offset: f32,
+    sample_period: u32,
+    harmonic: u32,
+    timestamps: [Option<i32>; 2],
+    iir: IIR,
+    iirstate: [IIRState; 2],
+}
+
+impl Lockin {
+    /// Initialize a new `Lockin` instance.
+    ///
+    /// # Arguments
+    ///
+    /// * `phase_offset` - Phase offset (in radians) applied to the
+    /// demodulation signal.
+    /// * `sample_period` - ADC sampling period in terms of the
+    /// internal clock period.
+    /// * `harmonic` - Integer scaling factor used to adjust the
+    /// demodulation frequency. E.g., 2 would demodulate with the
+    /// first harmonic.
+    /// * `iir` - IIR biquad filter.
+    ///
+    /// # Returns
+    ///
+    /// New `Lockin` instance.
+    pub fn new(
+        phase_offset: f32,
+        sample_period: u32,
+        harmonic: u32,
+        iir: IIR,
+    ) -> Self {
+        Lockin {
+            phase_offset: phase_offset,
+            sample_period: sample_period,
+            harmonic: harmonic,
+            timestamps: [None, None],
+            iir: iir,
+            iirstate: [[0.; 5]; 2],
+        }
+    }
+
+    /// Demodulate an input signal with the complex reference signal.
+    ///
+    /// # Arguments
+    ///
+    /// * `adc_samples` - One batch of ADC samples.
+    /// * `timestamps` - Counter values corresponding to the edges of
+    /// an external reference signal. The counter is incremented by a
+    /// fast internal clock.
+    ///
+    /// # Returns
+    ///
+    /// The demodulated complex signal as a `Result`. When there are
+    /// an insufficient number of timestamps to perform processing,
+    /// `Err` is returned.
+    ///
+    /// # Assumptions
+    ///
+    /// `demodulate` expects that the timestamp counter value is equal
+    /// to 0 when the ADC samples its first input in a batch. This can
+    /// be achieved by configuring the timestamp counter to overflow
+    /// at the end of the ADC batch sampling period.
+    pub fn demodulate(
+        &mut self,
+        adc_samples: &[i16],
+        timestamps: &[u16],
+    ) -> Result<[Complex<f32>; ADC_SAMPLE_BUFFER_SIZE], &str> {
+        let sample_period = self.sample_period as i32;
+        // update old timestamps for new ADC batch
+        self.timestamps.iter_mut().for_each(|t| match *t {
+            Some(timestamp) => {
+                // Existing timestamps have aged by one ADC batch
+                // period since the last ADC batch.
+                *t = Some(
+                    timestamp - ADC_SAMPLE_BUFFER_SIZE as i32 * sample_period,
+                );
+            }
+            None => (),
+        });
+
+        // return prematurely if there aren't enough timestamps for
+        // processing
+        let old_timestamp_count =
+            self.timestamps.iter().filter(|t| t.is_some()).count();
+        if old_timestamp_count + timestamps.len() < 2 {
+            return Err("insufficient timestamps");
+        }
+
+        let mut signal = [(0., 0.); ADC_SAMPLE_BUFFER_SIZE];
+        // if we have not yet recorded any timestamps, the first
+        // reference period must be computed from the first and
+        // second timestamps in the array
+        let mut timestamp_index: usize =
+            if old_timestamp_count == 0 { 1 } else { 0 };
+
+        // compute ADC sample phases, sines/cosines and demodulate
+        signal
+            .iter_mut()
+            .zip(adc_samples.iter())
+            .enumerate()
+            .for_each(|(i, (s, sample))| {
+                let adc_sample_count = i as i32 * sample_period;
+                // index of the closest timestamp that occurred after
+                // the current ADC sample
+                let closest_timestamp_after_index: i32 = if timestamps.len() > 0
+                {
+                    // Linear search is fast because both the timestamps
+                    // and ADC sample counts are sorted. Because of this,
+                    // we only need to check timestamps that were also
+                    // greater than the last ADC sample count.
+                    while timestamp_index < timestamps.len() - 1
+                        && (timestamps[timestamp_index] as i32)
+                            < adc_sample_count
+                    {
+                        timestamp_index += 1;
+                    }
+                    timestamp_index as i32
+                } else {
+                    -1
+                };
+
+                // closest timestamp that occurred before the current
+                // ADC sample
+                let closest_timestamp_before: i32;
+                let reference_period = if closest_timestamp_after_index < 0 {
+                    closest_timestamp_before = self.timestamps[0].unwrap();
+                    closest_timestamp_before - self.timestamps[1].unwrap()
+                } else if closest_timestamp_after_index == 0 {
+                    closest_timestamp_before = self.timestamps[0].unwrap();
+                    timestamps[0] as i32 - closest_timestamp_before
+                } else {
+                    closest_timestamp_before = timestamps
+                        [(closest_timestamp_after_index - 1) as usize]
+                        as i32;
+                    timestamps[closest_timestamp_after_index as usize] as i32
+                        - closest_timestamp_before
+                };
+
+                let integer_phase: i32 = (adc_sample_count
+                    - closest_timestamp_before)
+                    * self.harmonic as i32;
+                let phase = self.phase_offset
+                    + 2. * PI * integer_phase as f32 / reference_period as f32;
+                let (sine, cosine) = libm::sincosf(phase);
+                let sample = *sample as f32;
+                s.0 = sine * sample;
+                s.1 = cosine * sample;
+            });
+
+        // record new timestamps
+        let start_index: usize = if timestamps.len() < 2 {
+            0
+        } else {
+            timestamps.len() - 2
+        };
+        timestamps[start_index..]
+            .iter()
+            .for_each(|t| self.timestamps.push(Some(*t as i32)));
+
+        Ok(signal)
+    }
+
+    /// Filter the complex signal using the supplied biquad IIR. The
+    /// signal array is modified in place.
+    ///
+    /// # Arguments
+    ///
+    /// * `signal` - Complex signal to filter.
+    pub fn filter(&mut self, signal: &mut [Complex<f32>]) {
+        signal.iter_mut().for_each(|s| {
+            s.0 = self.iir.update(&mut self.iirstate[0], s.0);
+            s.1 = self.iir.update(&mut self.iirstate[1], s.1);
+        });
+    }
+}
+
+/// Decimate the complex signal to `DECIMATED_BUFFER_SIZE`. The ratio
+/// of `ADC_SAMPLE_BUFFER_SIZE` to `DECIMATED_BUFFER_SIZE` must be a
+/// power of 2.
+///
+/// # Arguments
+///
+/// * `signal` - Complex signal to decimate.
+///
+/// # Returns
+///
+/// The decimated signal.
+pub fn decimate(
+    signal: [Complex<f32>; ADC_SAMPLE_BUFFER_SIZE],
+) -> [Complex<f32>; DECIMATED_BUFFER_SIZE] {
+    let n_k = ADC_SAMPLE_BUFFER_SIZE / DECIMATED_BUFFER_SIZE;
+    debug_assert!(
+        ADC_SAMPLE_BUFFER_SIZE == DECIMATED_BUFFER_SIZE || n_k % 2 == 0
+    );
+
+    let mut signal_decimated = [(0_f32, 0_f32); DECIMATED_BUFFER_SIZE];
+
+    signal_decimated
+        .iter_mut()
+        .zip(signal.iter().step_by(n_k))
+        .for_each(|(s_d, s)| {
+            s_d.0 = s.0;
+            s_d.1 = s.1;
+        });
+
+    signal_decimated
+}
+
+/// Compute the magnitude and phase from the complex signal. The
+/// signal array is modified in place.
+///
+/// # Arguments
+///
+/// * `signal` - Complex signal to decimate.
+pub fn magnitude_phase(signal: &mut [Complex<f32>]) {
+    signal.iter_mut().for_each(|s| {
+        let new_i = libm::sqrtf([s.0, s.1].iter().map(|i| i * i).sum());
+        let new_q = libm::atan2f(s.1, s.0);
+        s.0 = new_i;
+        s.1 = new_q;
+    });
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::testing::complex_allclose;
+
+    #[test]
+    fn array_push() {
+        let mut arr: [Option<u32>; 2] = [None, None];
+        arr.push(Some(1));
+        assert_eq!(arr, [Some(1), None]);
+        arr.push(Some(2));
+        assert_eq!(arr, [Some(2), Some(1)]);
+        arr.push(Some(10));
+        assert_eq!(arr, [Some(10), Some(2)]);
+    }
+
+    #[test]
+    fn magnitude_phase_length_1_quadrant_1() {
+        let mut signal: [Complex<f32>; 1] = [(1., 1.)];
+        magnitude_phase(&mut signal);
+        assert!(complex_allclose(
+            &signal,
+            &[(2_f32.sqrt(), PI / 4.)],
+            f32::EPSILON,
+            0.
+        ));
+
+        signal = [(3_f32.sqrt() / 2., 1. / 2.)];
+        magnitude_phase(&mut signal);
+        assert!(complex_allclose(
+            &signal,
+            &[(1., PI / 6.)],
+            f32::EPSILON,
+            0.
+        ));
+    }
+
+    #[test]
+    fn magnitude_phase_length_1_quadrant_2() {
+        let mut signal = [(-1., 1.)];
+        magnitude_phase(&mut signal);
+        assert!(complex_allclose(
+            &signal,
+            &[(2_f32.sqrt(), 3. * PI / 4.)],
+            f32::EPSILON,
+            0.
+        ));
+
+        signal = [(-1. / 2., 3_f32.sqrt() / 2.)];
+        magnitude_phase(&mut signal);
+        assert!(complex_allclose(
+            &signal,
+            &[(1_f32, 2. * PI / 3.)],
+            f32::EPSILON,
+            0.
+        ));
+    }
+
+    #[test]
+    fn magnitude_phase_length_1_quadrant_3() {
+        let mut signal = [(-1. / 2_f32.sqrt(), -1. / 2_f32.sqrt())];
+        magnitude_phase(&mut signal);
+        assert!(complex_allclose(
+            &signal,
+            &[(1_f32.sqrt(), -3. * PI / 4.)],
+            f32::EPSILON,
+            0.
+        ));
+
+        signal = [(-1. / 2., -2_f32.sqrt())];
+        magnitude_phase(&mut signal);
+        assert!(complex_allclose(
+            &signal,
+            &[((3. / 2.) as f32, -1.91063323625 as f32)],
+            f32::EPSILON,
+            0.
+        ));
+    }
+
+    #[test]
+    fn magnitude_phase_length_1_quadrant_4() {
+        let mut signal = [(1. / 2_f32.sqrt(), -1. / 2_f32.sqrt())];
+        magnitude_phase(&mut signal);
+        assert!(complex_allclose(
+            &signal,
+            &[(1_f32.sqrt(), -1. * PI / 4.)],
+            f32::EPSILON,
+            0.
+        ));
+
+        signal = [(3_f32.sqrt() / 2., -1. / 2.)];
+        magnitude_phase(&mut signal);
+        assert!(complex_allclose(
+            &signal,
+            &[(1_f32, -PI / 6.)],
+            f32::EPSILON,
+            0.
+        ));
+    }
+
+    #[test]
+    fn decimate_sample_16_decimated_1() {
+        let signal: [Complex<f32>; ADC_SAMPLE_BUFFER_SIZE] = [
+            (0.0, 1.6),
+            (0.1, 1.7),
+            (0.2, 1.8),
+            (0.3, 1.9),
+            (0.4, 2.0),
+            (0.5, 2.1),
+            (0.6, 2.2),
+            (0.7, 2.3),
+            (0.8, 2.4),
+            (0.9, 2.5),
+            (1.0, 2.6),
+            (1.1, 2.7),
+            (1.2, 2.8),
+            (1.3, 2.9),
+            (1.4, 3.0),
+            (1.5, 3.1),
+        ];
+        assert_eq!(decimate(signal), [(0.0, 1.6)]);
+    }
+
+    #[test]
+    fn lockin_demodulate_valid_0() {
+        let mut lockin = Lockin::new(
+            0.,
+            200,
+            1,
+            IIR {
+                ba: [0_f32; 5],
+                y_offset: 0.,
+                y_min: -(1 << 15) as f32,
+                y_max: (1 << 15) as f32 - 1.,
+            },
+        );
+        assert_eq!(
+            lockin.demodulate(&[0; ADC_SAMPLE_BUFFER_SIZE], &[]),
+            Err("insufficient timestamps")
+        );
+    }
+
+    #[test]
+    fn lockin_demodulate_valid_1() {
+        let mut lockin = Lockin::new(
+            0.,
+            200,
+            1,
+            IIR {
+                ba: [0_f32; 5],
+                y_offset: 0.,
+                y_min: -(1 << 15) as f32,
+                y_max: (1 << 15) as f32 - 1.,
+            },
+        );
+        assert_eq!(
+            lockin.demodulate(&[0; ADC_SAMPLE_BUFFER_SIZE], &[0],),
+            Err("insufficient timestamps")
+        );
+    }
+
+    #[test]
+    fn lockin_demodulate_valid_2() {
+        let adc_period: u32 = 200;
+        let mut lockin = Lockin::new(
+            0.,
+            adc_period,
+            1,
+            IIR {
+                ba: [0_f32; 5],
+                y_offset: 0.,
+                y_min: -(1 << 15) as f32,
+                y_max: (1 << 15) as f32 - 1.,
+            },
+        );
+        let adc_samples: [i16; ADC_SAMPLE_BUFFER_SIZE] =
+            [-8, 7, -7, 6, -6, 5, -5, 4, -4, 3, -3, 2, -2, -1, 1, 0];
+        let reference_period: u16 = 2800;
+        let initial_phase_integer: u16 = 200;
+        let timestamps: &[u16] = &[
+            initial_phase_integer,
+            initial_phase_integer + reference_period,
+        ];
+        let initial_phase: f32 =
+            -(initial_phase_integer as f32) / reference_period as f32 * 2. * PI;
+        let phase_increment: f32 =
+            adc_period as f32 / reference_period as f32 * 2. * PI;
+        let mut signal = [(0., 0.); ADC_SAMPLE_BUFFER_SIZE];
+        for (n, s) in signal.iter_mut().enumerate() {
+            let adc_phase = initial_phase + n as f32 * phase_increment;
+            let sine = adc_phase.sin();
+            let cosine = adc_phase.cos();
+            s.0 = sine * adc_samples[n] as f32;
+            s.1 = cosine * adc_samples[n] as f32;
+        }
+        let result = lockin.demodulate(&adc_samples, timestamps).unwrap();
+        assert!(
+            complex_allclose(&result, &signal, 0., 1e-5),
+            "\nsignal computed: {:?},\nsignal expected: {:?}",
+            result,
+            signal
+        );
+    }
+}

dsp/src/pll.rs

@@ -0,0 +1,97 @@
+use serde::{Deserialize, Serialize};
+
+/// Type-II, sampled phase, discrete time PLL
+///
+/// This PLL tracks the frequency and phase of an input signal with respect to the sampling clock.
+/// The transfer function is I^2,I from input phase to output phase and P,I from input phase to
+/// output frequency.
+///
+/// The PLL locks to any frequency (i.e. it locks to the alias in the first Nyquist zone) and is
+/// stable for any gain (1 <= shift <= 30). It has a single parameter that determines the loop
+/// bandwidth in octave steps. The gain can be changed freely between updates.
+///
+/// The frequency and phase settling time constants for an (any) frequency jump are `1 << shift`
+/// update cycles. The loop bandwidth is about `1/(2*pi*(1 << shift))` in units of the sample rate.
+///
+/// All math is naturally wrapping 32 bit integer. Phase and frequency are understood modulo that
+/// overflow in the first Nyquist zone. Expressing the IIR equations in other ways (e.g. single
+/// (T)-DF-{I,II} biquad/IIR) would break on overflow.
+///
+/// There are no floating point rounding errors here. But there is integer quantization/truncation
+/// error of the `shift` lowest bits leading to a phase offset for very low gains. Truncation
+/// bias is applied. Rounding is "half up". The phase truncation error can be removed very
+/// efficiently by dithering.
+///
+/// This PLL does not unwrap phase slips during lock acquisition. This can and should be
+/// implemented elsewhere by (down) scaling and then unwrapping the input phase and (up) scaling
+/// and wrapping output phase and frequency. This affects dynamic range accordingly.
+///
+/// The extension to I^3,I^2,I behavior to track chirps phase-accurately or to i64 data to
+/// increase resolution for extremely narrowband applications is obvious.
+#[derive(Copy, Clone, Default, Deserialize, Serialize)]
+pub struct PLL {
+    // last input phase
+    x: i32,
+    // filtered frequency
+    f: i32,
+    // filtered output phase
+    y: i32,
+}
+
+impl PLL {
+    /// Update the PLL with a new phase sample.
+    ///
+    /// Args:
+    /// * `input`: New input phase sample.
+    /// * `shift`: Error scaling. The frequency gain per update is `1/(1 << shift)`. The phase gain
+    ///   is always twice the frequency gain.
+    ///
+    /// Returns:
+    /// A tuple of instantaneous phase and frequency (the current phase increment).
+    pub fn update(&mut self, x: i32, shift: u8) -> (i32, i32) {
+        debug_assert!((1..=30).contains(&shift));
+        let bias = 1i32 << shift;
+        let e = x.wrapping_sub(self.f);
+        self.f = self.f.wrapping_add(
+            (bias >> 1).wrapping_add(e).wrapping_sub(self.x) >> shift,
+        );
+        self.x = x;
+        let f = self.f.wrapping_add(
+            bias.wrapping_add(e).wrapping_sub(self.y) >> (shift - 1),
+        );
+        self.y = self.y.wrapping_add(f);
+        (self.y, f)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    #[test]
+    fn mini() {
+        let mut p = PLL::default();
+        let (y, f) = p.update(0x10000, 10);
+        assert_eq!(y, 0xc2);
+        assert_eq!(f, y);
+    }
+
+    #[test]
+    fn converge() {
+        let mut p = PLL::default();
+        let f0 = 0x71f63049_i32;
+        let shift = 10;
+        let n = 31 << shift + 2;
+        let mut x = 0i32;
+        for i in 0..n {
+            x = x.wrapping_add(f0);
+            let (y, f) = p.update(x, shift);
+            if i > n / 4 {
+                assert_eq!(f.wrapping_sub(f0).abs() <= 1, true);
+            }
+            if i > n / 2 {
+                // The remaining error is removed by dithering.
+                assert_eq!(y.wrapping_sub(x).abs() < 1 << 18, true);
+            }
+        }
+    }
+}

dsp/src/testing.rs

@@ -0,0 +1,27 @@
+use super::Complex;
+
+pub fn isclose(a: f32, b: f32, rtol: f32, atol: f32) -> bool {
+    (a - b).abs() <= a.abs().max(b.abs()) * rtol + atol
+}
+
+pub fn complex_isclose(
+    a: Complex<f32>,
+    b: Complex<f32>,
+    rtol: f32,
+    atol: f32,
+) -> bool {
+    isclose(a.0, b.0, rtol, atol) && isclose(a.1, b.1, rtol, atol)
+}
+
+pub fn complex_allclose(
+    a: &[Complex<f32>],
+    b: &[Complex<f32>],
+    rtol: f32,
+    atol: f32,
+) -> bool {
+    let mut result: bool = true;
+    a.iter().zip(b.iter()).for_each(|(i, j)| {
+        result &= complex_isclose(*i, *j, rtol, atol);
+    });
+    result
+}

dsp/src/unwrap.rs

@@ -0,0 +1,83 @@
+use serde::{Deserialize, Serialize};
+
+/// Subtract `y - x` with signed overflow.
+///
+/// This is very similar to `i32::overflowing_sub(y, x)` except that the
+/// overflow indicator is not a boolean but the signum of the overflow.
+/// Additionally it's typically faster.
+///
+/// Returns:
+/// A tuple containg the (wrapped) difference `y - x` and the signum of the
+/// overflow.
+#[inline(always)]
+pub fn overflowing_sub(y: i32, x: i32) -> (i32, i8) {
+    let delta = y.wrapping_sub(x);
+    let wrap = (delta >= 0) as i8 - (y >= x) as i8;
+    (delta, wrap)
+}
+
+/// Overflow unwrapper.
+///
+/// This is unwrapping as in the phase and overflow unwrapping context, not
+/// unwrapping as in the `Result`/`Option` context.
+#[derive(Copy, Clone, Default, Deserialize, Serialize)]
+pub struct Unwrapper {
+    // last input
+    x: i32,
+    // last wraps
+    w: i32,
+}
+
+impl Unwrapper {
+    /// Unwrap a new sample from a sequence and update the unwrapper state.
+    ///
+    /// Args:
+    /// * `x`: New sample
+    ///
+    /// Returns:
+    /// A tuple containing the (wrapped) difference `x - x_old` and the
+    /// signed number of wraps accumulated by the new sample.
+    pub fn update(&mut self, x: i32) -> (i32, i32) {
+        let (dx, dw) = overflowing_sub(x, self.x);
+        self.x = x;
+        self.w = self.w.wrapping_add(dw as i32);
+        (dx, self.w)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    #[test]
+    fn mini() {
+        for (x0, x1, v) in [
+            (0i32, 0i32, 0i8),
+            (0, 1, 0),
+            (0, -1, 0),
+            (1, 0, 0),
+            (-1, 0, 0),
+            (0, 0x7fff_ffff, 0),
+            (-1, 0x7fff_ffff, -1),
+            (-2, 0x7fff_ffff, -1),
+            (-1, -0x8000_0000, 0),
+            (0, -0x8000_0000, 0),
+            (1, -0x8000_0000, 1),
+            (-0x6000_0000, 0x6000_0000, -1),
+            (0x6000_0000, -0x6000_0000, 1),
+            (-0x4000_0000, 0x3fff_ffff, 0),
+            (-0x4000_0000, 0x4000_0000, -1),
+            (-0x4000_0000, 0x4000_0001, -1),
+            (0x4000_0000, -0x3fff_ffff, 0),
+            (0x4000_0000, -0x4000_0000, 0),
+            (0x4000_0000, -0x4000_0001, 1),
+        ]
+        .iter()
+        {
+            let (dx, w) = overflowing_sub(*x1, *x0);
+            assert_eq!(*v, w, " = overflowing_sub({:#x}, {:#x})", *x0, *x1);
+            let (dx0, w0) = x1.overflowing_sub(*x0);
+            assert_eq!(w0, w != 0);
+            assert_eq!(dx, dx0);
+        }
+    }
+}

src/main.rs

@@ -60,6 +60,9 @@ const SAMPLE_FREQUENCY_KHZ: u32 = 500;
 // The desired ADC sample processing buffer size.
 const SAMPLE_BUFFER_SIZE: usize = 1;
 
+// The number of cascaded IIR biquads per channel. Select 1 or 2!
+const IIR_CASCADE_LENGTH: usize = 1;
+
 #[link_section = ".sram3.eth"]
 static mut DES_RING: ethernet::DesRing = ethernet::DesRing::new();
 
@@ -214,10 +217,11 @@ const APP: () = {
 
         pounder: Option<pounder::PounderDevices>,
 
-        #[init([[0.; 5]; 2])]
+        // Format: iir_state[ch][cascade-no][coeff]
-        iir_state: [iir::IIRState; 2],
+        #[init([[[0.; 5]; IIR_CASCADE_LENGTH]; 2])]
-        #[init([iir::IIR { ba: [1., 0., 0., 0., 0.], y_offset: 0., y_min: -SCALE - 1., y_max: SCALE }; 2])]
+        iir_state: [[iir::IIRState; IIR_CASCADE_LENGTH]; 2],
-        iir_ch: [iir::IIR; 2],
+        #[init([[iir::IIR { ba: [1., 0., 0., 0., 0.], y_offset: 0., y_min: -SCALE - 1., y_max: SCALE }; IIR_CASCADE_LENGTH]; 2])]
+        iir_ch: [[iir::IIR; IIR_CASCADE_LENGTH]; 2],
     }
 
     #[init]
@@ -809,8 +813,11 @@ const APP: () = {
         for channel in 0..adc_samples.len() {
             for sample in 0..adc_samples[0].len() {
                 let x = f32::from(adc_samples[channel][sample] as i16);
-                let y = c.resources.iir_ch[channel]
+                let mut y = x;
-                    .update(&mut c.resources.iir_state[channel], x);
+                for i in 0..c.resources.iir_state[channel].len() {
+                    y = c.resources.iir_ch[channel][i]
+                        .update(&mut c.resources.iir_state[channel][i], y);
+                }
                 // Note(unsafe): The filter limits ensure that the value is in range.
                 // The truncation introduces 1/2 LSB distortion.
                 let y = unsafe { y.to_int_unchecked::<i16>() };
@@ -887,10 +894,23 @@ const APP: () = {
                                     let state = c.resources.iir_state.lock(|iir_state|
                                         server::Status {
                                             t: time,
-                                            x0: iir_state[0][0],
+                                            x0: iir_state[0][0][0],
-                                            y0: iir_state[0][2],
+                                            y0: iir_state[0][0][2],
-                                            x1: iir_state[1][0],
+                                            x1: iir_state[1][0][0],
-                                            y1: iir_state[1][2],
+                                            y1: iir_state[1][0][2],
+                                    });
+
+                                    Ok::<server::Status, ()>(state)
+                                }),
+                                // "_b" means cascades 2nd IIR
+                                "stabilizer/iir_b/state": (|| {
+                                    let state = c.resources.iir_state.lock(|iir_state|
+                                        server::Status {
+                                            t: time,
+                                            x0: iir_state[0][IIR_CASCADE_LENGTH-1][0],
+                                            y0: iir_state[0][IIR_CASCADE_LENGTH-1][2],
+                                            x1: iir_state[1][IIR_CASCADE_LENGTH-1][0],
+                                            y1: iir_state[1][IIR_CASCADE_LENGTH-1][2],
                                     });
 
                                     Ok::<server::Status, ()>(state)
@@ -906,7 +926,7 @@ const APP: () = {
                                             return Err(());
                                         }
 
-                                        iir_ch[req.channel as usize] = req.iir;
+                                        iir_ch[req.channel as usize][0] = req.iir;
 
                                         Ok::<server::IirRequest, ()>(req)
                                     })
@@ -917,7 +937,29 @@ const APP: () = {
                                             return Err(());
                                         }
 
-                                        iir_ch[req.channel as usize] = req.iir;
+                                        iir_ch[req.channel as usize][0] = req.iir;
+
+                                        Ok::<server::IirRequest, ()>(req)
+                                    })
+                                }),
+                                "stabilizer/iir_b0/state": server::IirRequest, (|req: server::IirRequest| {
+                                    c.resources.iir_ch.lock(|iir_ch| {
+                                        if req.channel > 1 {
+                                            return Err(());
+                                        }
+
+                                        iir_ch[req.channel as usize][IIR_CASCADE_LENGTH-1] = req.iir;
+
+                                        Ok::<server::IirRequest, ()>(req)
+                                    })
+                                }),
+                                "stabilizer/iir_b1/state": server::IirRequest,(|req: server::IirRequest| {
+                                    c.resources.iir_ch.lock(|iir_ch| {
+                                        if req.channel > 1 {
+                                            return Err(());
+                                        }
+
+                                        iir_ch[req.channel as usize][IIR_CASCADE_LENGTH-1] = req.iir;
 
                                         Ok::<server::IirRequest, ()>(req)
                                     })