update lock-in for integer math and PLL

dsp/src/lib.rs

@@ -5,7 +5,7 @@ use core::ops::{Add, Mul, Neg};
 
 pub type Complex<T> = (T, T);
 
-/// Round up half.
+/// Bit shift, round up half.
 ///
 /// # Arguments
 ///
@@ -20,6 +20,25 @@ pub fn shift_round(x: i32, shift: usize) -> i32 {
     (x + (1 << (shift - 1))) >> shift
 }
 
+/// Integer division, round up half.
+///
+/// # Arguments
+///
+/// `dividend` - Value to divide.
+/// `divisor` - Value that divides the dividend.
+///
+/// # Returns
+///
+/// Divided and rounded value.
+#[inline(always)]
+pub fn divide_round(dividend: i64, divisor: i64) -> i64 {
+    debug_assert!(
+        dividend as i128 + (divisor as i128 - 1) < i64::MAX as i128
+            && dividend as i128 + (divisor as i128 - 1) > i64::MIN as i128
+    );
+    (dividend + (divisor - 1)) / divisor
+}
+
 fn abs<T>(x: T) -> T
 where
     T: PartialOrd + Default + Neg<Output = T>,

dsp/src/lockin.rs

@@ -52,47 +52,36 @@
 //! the demodulation frequency. This does not require any state
 //! information and is therefore a normal function.
 
-use super::iir::{IIRState, IIR};
+use super::iir_int::{IIRState, IIR};
-use super::Complex;
+use super::pll::PLL;
-use core::f32::consts::PI;
+use super::trig::{atan2, cossin};
+use super::{divide_round, Complex};
+
+/// TODO these constants are copied from main.rs and should be
+/// shared. Additionally, we should probably store the log2 values and
+/// compute the actual values from these in main, as is done here.
+pub const SAMPLE_BUFFER_SIZE_LOG2: usize = 0;
+pub const SAMPLE_BUFFER_SIZE: usize = 1 << SAMPLE_BUFFER_SIZE_LOG2;
+
+pub const ADC_SAMPLE_TICKS_LOG2: usize = 8;
+pub const ADC_SAMPLE_TICKS: usize = 1 << ADC_SAMPLE_TICKS_LOG2;
+
+pub const ADC_BATCHES_LOG2: usize =
+    32 - SAMPLE_BUFFER_SIZE_LOG2 - ADC_SAMPLE_TICKS_LOG2;
+pub const ADC_BATCHES: usize = 1 << ADC_BATCHES_LOG2;
 
-/// 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],
+    phase_offset: u32,
+    batch_index: u32,
+    last_phase: Option<i64>,
+    last_frequency: Option<i64>,
+    pll: PLL,
+    pll_shift_frequency: u8,
+    pll_shift_phase: u8,
     iir: IIR,
     iirstate: [IIRState; 2],
 }
@@ -102,31 +91,36 @@ impl Lockin {
     ///
     /// # 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.
+    /// * `phase_offset` - Phase offset applied to the demodulation
+    /// signal.
     /// * `iir` - IIR biquad filter.
+    /// * `pll_shift_frequency` - See PLL::update().
+    /// * `pll_shift_phase` - See PLL::update().
     ///
     /// # Returns
     ///
     /// New `Lockin` instance.
     pub fn new(
-        phase_offset: f32,
-        sample_period: u32,
         harmonic: u32,
+        phase_offset: u32,
         iir: IIR,
+        pll_shift_frequency: u8,
+        pll_shift_phase: u8,
     ) -> Self {
         Lockin {
-            phase_offset: phase_offset,
+            harmonic,
-            sample_period: sample_period,
+            phase_offset,
-            harmonic: harmonic,
+            batch_index: 0,
-            timestamps: [None, None],
+            last_phase: None,
-            iir: iir,
+            last_frequency: None,
-            iirstate: [[0.; 5]; 2],
+            pll: PLL::default(),
+            pll_shift_frequency,
+            pll_shift_phase,
+            iir,
+            iirstate: [[0; 5]; 2],
         }
     }
 
@@ -135,120 +129,88 @@ impl Lockin {
     /// # Arguments
     ///
     /// * `adc_samples` - One batch of ADC samples.
-    /// * `timestamps` - Counter values corresponding to the edges of
+    /// * `timestamp` - Counter value corresponding to the edges of an
-    /// an external reference signal. The counter is incremented by a
+    /// external reference signal. The counter is incremented by a
-    /// fast internal clock.
+    /// fast internal clock. Each ADC sample batch can contain 0 or 1
+    /// timestamps.
     ///
     /// # 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],
+        timestamp: Option<u32>,
-    ) -> Result<[Complex<f32>; ADC_SAMPLE_BUFFER_SIZE], &str> {
+    ) -> Result<[Complex<i32>; SAMPLE_BUFFER_SIZE], &str> {
-        let sample_period = self.sample_period as i32;
+        let frequency: i64;
-        // update old timestamps for new ADC batch
+        let phase: i64;
-        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
+        match timestamp {
-        // processing
+            Some(t) => {
-        let old_timestamp_count =
+                let res = self.pll.update(
-            self.timestamps.iter().filter(|t| t.is_some()).count();
+                    t as i32,
-        if old_timestamp_count + timestamps.len() < 2 {
+                    self.pll_shift_frequency,
-            return Err("insufficient timestamps");
+                    self.pll_shift_phase,
+                );
+                phase = res.0 as u32 as i64;
+                frequency = res.1 as u32 as i64;
+                self.last_phase = Some(phase);
+                self.last_frequency = Some(frequency);
+            }
+            None => match self.last_phase {
+                Some(t) => {
+                    phase = t;
+                    frequency = self.last_frequency.unwrap();
+                }
+                None => {
+                    self.batch_index += 1;
+                    if self.batch_index == ADC_BATCHES as u32 {
+                        self.batch_index = 0;
+                    }
+                    return Err("insufficient timestamps");
+                }
+            },
         }
 
-        let mut signal = [(0., 0.); ADC_SAMPLE_BUFFER_SIZE];
+        let demodulation_frequency = divide_round(
-        // if we have not yet recorded any timestamps, the first
+            1 << (64 - SAMPLE_BUFFER_SIZE_LOG2 - ADC_BATCHES_LOG2),
-        // reference period must be computed from the first and
+            frequency,
-        // second timestamps in the array
+        ) as u32;
-        let mut timestamp_index: usize =
+        let demodulation_initial_phase = divide_round(
-            if old_timestamp_count == 0 { 1 } else { 0 };
+            (((self.batch_index as i64) << (32 - ADC_BATCHES_LOG2)) - phase)
+                << 32,
+            frequency,
+        ) as u32;
 
-        // compute ADC sample phases, sines/cosines and demodulate
+        let mut demodulation_signal = [(0_i32, 0_i32); SAMPLE_BUFFER_SIZE];
-        signal
+
+        demodulation_signal
             .iter_mut()
             .zip(adc_samples.iter())
             .enumerate()
             .for_each(|(i, (s, sample))| {
-                let adc_sample_count = i as i32 * sample_period;
+                let sample_phase = (self.harmonic.wrapping_mul(
-                // index of the closest timestamp that occurred after
+                    (demodulation_frequency.wrapping_mul(i as u32))
-                // the current ADC sample
+                        .wrapping_add(demodulation_initial_phase),
-                let closest_timestamp_after_index: i32 = if timestamps.len() > 0
+                ))
-                {
+                .wrapping_add(self.phase_offset);
-                    // Linear search is fast because both the timestamps
+                let (cos, sin) = cossin(sample_phase as i32);
-                    // and ADC sample counts are sorted. Because of this,
+                // cos/sin take up 32 bits and sample takes up 16
-                    // we only need to check timestamps that were also
+                // bits. Make this fit into a 32 bit result.
-                    // greater than the last ADC sample count.
+                s.0 = ((*sample as i64 * cos as i64) >> 16) as i32;
-                    while timestamp_index < timestamps.len() - 1
+                s.1 = ((*sample as i64 * sin as i64) >> 16) as i32;
-                        && (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
+        if self.batch_index < ADC_BATCHES as u32 - 1 {
-        let start_index: usize = if timestamps.len() < 2 {
+            self.batch_index += 1;
-            0
         } else {
-            timestamps.len() - 2
+            self.batch_index = 0;
-        };
+            self.last_phase = Some(self.last_phase.unwrap() - (1 << 32));
-        timestamps[start_index..]
+        }
-            .iter()
-            .for_each(|t| self.timestamps.push(Some(*t as i32)));
 
-        Ok(signal)
+        Ok(demodulation_signal)
     }
 
     /// Filter the complex signal using the supplied biquad IIR. The
@@ -257,7 +219,7 @@ impl Lockin {
     /// # Arguments
     ///
     /// * `signal` - Complex signal to filter.
-    pub fn filter(&mut self, signal: &mut [Complex<f32>]) {
+    pub fn filter(&mut self, signal: &mut [Complex<i32>]) {
         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);
@@ -266,8 +228,8 @@ impl Lockin {
 }
 
 /// Decimate the complex signal to `DECIMATED_BUFFER_SIZE`. The ratio
-/// of `ADC_SAMPLE_BUFFER_SIZE` to `DECIMATED_BUFFER_SIZE` must be a
+/// of `SAMPLE_BUFFER_SIZE` to `DECIMATED_BUFFER_SIZE` must be a power
-/// power of 2.
+/// of 2.
 ///
 /// # Arguments
 ///
@@ -277,14 +239,12 @@ impl Lockin {
 ///
 /// The decimated signal.
 pub fn decimate(
-    signal: [Complex<f32>; ADC_SAMPLE_BUFFER_SIZE],
+    signal: [Complex<i32>; SAMPLE_BUFFER_SIZE],
-) -> [Complex<f32>; DECIMATED_BUFFER_SIZE] {
+) -> [Complex<i32>; DECIMATED_BUFFER_SIZE] {
-    let n_k = ADC_SAMPLE_BUFFER_SIZE / DECIMATED_BUFFER_SIZE;
+    let n_k = SAMPLE_BUFFER_SIZE / DECIMATED_BUFFER_SIZE;
-    debug_assert!(
+    debug_assert!(SAMPLE_BUFFER_SIZE == DECIMATED_BUFFER_SIZE || n_k % 2 == 0);
-        ADC_SAMPLE_BUFFER_SIZE == DECIMATED_BUFFER_SIZE || n_k % 2 == 0
-    );
 
-    let mut signal_decimated = [(0_f32, 0_f32); DECIMATED_BUFFER_SIZE];
+    let mut signal_decimated = [(0_i32, 0_i32); DECIMATED_BUFFER_SIZE];
 
     signal_decimated
         .iter_mut()
@@ -302,11 +262,13 @@ pub fn decimate(
 ///
 /// # Arguments
 ///
-/// * `signal` - Complex signal to decimate.
+/// * `signal` - Complex signal for which the magnitude and phase
-pub fn magnitude_phase(signal: &mut [Complex<f32>]) {
+/// should be computed. TODO currently, we compute the square of the
+/// magnitude. This should be changed to be the actual magnitude.
+pub fn magnitude_phase(signal: &mut [Complex<i32>]) {
     signal.iter_mut().for_each(|s| {
-        let new_i = libm::sqrtf([s.0, s.1].iter().map(|i| i * i).sum());
+        let new_i = [s.0, s.1].iter().map(|i| i * i).sum();
-        let new_q = libm::atan2f(s.1, s.0);
+        let new_q = atan2(s.1, s.0);
         s.0 = new_i;
         s.1 = new_q;
     });
@@ -315,204 +277,115 @@ pub fn magnitude_phase(signal: &mut [Complex<f32>]) {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::testing::complex_allclose;
 
     #[test]
-    fn array_push() {
+    /// Ensure that the demodulation signals are within some tolerance
-        let mut arr: [Option<u32>; 2] = [None, None];
+    /// band of the target value given the phase and frequency values
-        arr.push(Some(1));
+    /// provided by the PLL.
-        assert_eq!(arr, [Some(1), None]);
+    fn demodulate() {
-        arr.push(Some(2));
+        const PLL_SHIFT_FREQUENCY: u8 = 4;
-        assert_eq!(arr, [Some(2), Some(1)]);
+        const PLL_SHIFT_PHASE: u8 = 3;
-        arr.push(Some(10));
+        const HARMONIC: u32 = 1;
-        assert_eq!(arr, [Some(10), Some(2)]);
+        const PHASE_OFFSET: u32 = 0;
-    }
-
-    #[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.,
+            HARMONIC,
-            200,
+            PHASE_OFFSET,
-            1,
+            IIR { ba: [0; 5] },
-            IIR {
+            PLL_SHIFT_FREQUENCY,
-                ba: [0_f32; 5],
+            PLL_SHIFT_PHASE,
-                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]
+        // Duplicate the PLL outside demodulate so that we don't test
-    fn lockin_demodulate_valid_1() {
+        // its behavior.
-        let mut lockin = Lockin::new(
+        let mut tracking_pll = PLL::default();
-            0.,
+        let mut tracking_phase: i32 = 0;
-            200,
+        let mut tracking_frequency: i32 = 0;
-            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]
+        const REFERENCE_FREQUENCY: usize = 10_000;
-    fn lockin_demodulate_valid_2() {
+        let mut reference_edge: usize = REFERENCE_FREQUENCY;
-        let adc_period: u32 = 200;
+
-        let mut lockin = Lockin::new(
+        // Ensure that we receive at most 1 timestamp per batch.
-            0.,
+        debug_assert!(
-            adc_period,
+            REFERENCE_FREQUENCY >= SAMPLE_BUFFER_SIZE * ADC_SAMPLE_TICKS
-            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];
+        for batch in 0..100_000 {
-        let reference_period: u16 = 2800;
+            let tick: usize = batch * ADC_SAMPLE_TICKS * SAMPLE_BUFFER_SIZE;
-        let initial_phase_integer: u16 = 200;
+            let timestamp: Option<u32>;
-        let timestamps: &[u16] = &[
+
-            initial_phase_integer,
+            // When the reference edge occurred during the current
-            initial_phase_integer + reference_period,
+            // batch acquisition, register the timestamp and update
-        ];
+            // the tracking PLL.
-        let initial_phase: f32 =
+            if reference_edge >= tick
-            -(initial_phase_integer as f32) / reference_period as f32 * 2. * PI;
+                && reference_edge < tick + ADC_SAMPLE_TICKS * SAMPLE_BUFFER_SIZE
-        let phase_increment: f32 =
+            {
-            adc_period as f32 / reference_period as f32 * 2. * PI;
+                timestamp = Some(reference_edge as u32);
-        let mut signal = [(0., 0.); ADC_SAMPLE_BUFFER_SIZE];
+
-        for (n, s) in signal.iter_mut().enumerate() {
+                let tracking_update = tracking_pll.update(
-            let adc_phase = initial_phase + n as f32 * phase_increment;
+                    reference_edge as i32,
-            let sine = adc_phase.sin();
+                    PLL_SHIFT_FREQUENCY,
-            let cosine = adc_phase.cos();
+                    PLL_SHIFT_PHASE,
-            s.0 = sine * adc_samples[n] as f32;
+                );
-            s.1 = cosine * adc_samples[n] as f32;
+                tracking_phase = tracking_update.0;
+                tracking_frequency = tracking_update.1;
+
+                reference_edge += REFERENCE_FREQUENCY;
+            } else {
+                timestamp = None;
+            }
+
+            let timestamp_before_batch = if tracking_phase > tick as i32 {
+                // There can be at most 1 reference edge per batch, so
+                // this will necessarily place the timestamp prior to
+                // the current batch.
+                tracking_phase - tracking_frequency
+            } else {
+                tracking_phase
+            };
+
+            let initial_phase = (((tick as f64
+                - timestamp_before_batch as f64)
+                / tracking_frequency as f64
+                * (1_i64 << 32) as f64)
+                .round()
+                % u32::MAX as f64) as u32;
+            let frequency = ((ADC_SAMPLE_TICKS as f64
+                / tracking_frequency as f64
+                * (1_i64 << 32) as f64)
+                .round()
+                % u32::MAX as f64) as u32;
+
+            match lockin.demodulate(&[i16::MAX; SAMPLE_BUFFER_SIZE], timestamp)
+            {
+                Ok(v) => {
+                    println!("batch      : {}", batch);
+                    for sample in 0..SAMPLE_BUFFER_SIZE {
+                        const TOL: i32 = 50_000;
+                        let cos = v[sample].0;
+                        let sin = v[sample].1;
+
+                        let (mut target_cos, mut target_sin) = cossin(
+                            HARMONIC
+                                .wrapping_mul(
+                                    (frequency.wrapping_mul(sample as u32))
+                                        .wrapping_add(initial_phase),
+                                )
+                                .wrapping_add(PHASE_OFFSET)
+                                as i32,
+                        );
+                        target_cos /= 2;
+                        target_sin /= 2;
+
+                        println!("sample     : {}", sample);
+                        println!("tol        : {}", TOL);
+                        println!("cos, target: {}, {}", cos, target_cos);
+                        println!("sin, target: {}, {}", sin, target_sin);
+                        assert!((cos - target_cos).abs() < TOL);
+                        assert!((sin - target_sin).abs() < TOL);
+                    }
+                }
+                Err(_) => {}
+            }
         }
-        let result = lockin.demodulate(&adc_samples, timestamps).unwrap();
-        assert!(
-            complex_allclose(&result, &signal, 0., 1e-5),
-            "\nsignal computed: {:?},\nsignal expected: {:?}",
-            result,
-            signal
-        );
     }
 }