diff --git a/dsp/src/lib.rs b/dsp/src/lib.rs index 1161943..3f090b5 100644 --- a/dsp/src/lib.rs +++ b/dsp/src/lib.rs @@ -116,7 +116,6 @@ where pub mod iir; pub mod iir_int; -pub mod lockin; pub mod pll; pub mod trig; pub mod unwrap; diff --git a/dsp/src/lockin.rs b/dsp/src/lockin.rs deleted file mode 100644 index 4ab6f02..0000000 --- a/dsp/src/lockin.rs +++ /dev/null @@ -1,392 +0,0 @@ -//! 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_int::{IIRState, IIR}; -use super::pll::PLL; -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; - -pub const DECIMATED_BUFFER_SIZE: usize = 1; - -/// Performs lock-in amplifier processing of a signal. -pub struct Lockin { - harmonic: u32, - phase_offset: u32, - batch_index: u32, - last_phase: Option, - last_frequency: Option, - pll: PLL, - pll_shift_frequency: u8, - pll_shift_phase: u8, - iir: IIR, - iirstate: [IIRState; 2], -} - -impl Lockin { - /// Initialize a new `Lockin` instance. - /// - /// # Arguments - /// - /// * `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( - harmonic: u32, - phase_offset: u32, - iir: IIR, - pll_shift_frequency: u8, - pll_shift_phase: u8, - ) -> Self { - Lockin { - harmonic, - phase_offset, - batch_index: 0, - last_phase: None, - last_frequency: None, - pll: PLL::default(), - pll_shift_frequency, - pll_shift_phase, - iir, - iirstate: [[0; 5]; 2], - } - } - - /// Demodulate an input signal with the complex reference signal. - /// - /// # Arguments - /// - /// * `adc_samples` - One batch of ADC samples. - /// * `timestamp` - Counter value corresponding to the edges of an - /// external reference signal. The counter is incremented by a - /// 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. - pub fn demodulate( - &mut self, - adc_samples: &[i16], - timestamp: Option, - ) -> Result<[Complex; SAMPLE_BUFFER_SIZE], &str> { - let frequency: i64; - let phase: i64; - - match timestamp { - Some(t) => { - let res = self.pll.update( - t as i32, - self.pll_shift_frequency, - 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 => { - if self.batch_index < ADC_BATCHES as u32 - 1 { - self.batch_index += 1; - } else { - self.batch_index = 0; - } - return Err("insufficient timestamps"); - } - }, - } - - let demodulation_frequency = divide_round( - 1 << (64 - SAMPLE_BUFFER_SIZE_LOG2 - ADC_BATCHES_LOG2), - frequency, - ) as u32; - let demodulation_initial_phase = divide_round( - (((self.batch_index as i64) << (32 - ADC_BATCHES_LOG2)) - phase) - << 32, - frequency, - ) as u32; - - let mut demodulation_signal = [(0_i32, 0_i32); SAMPLE_BUFFER_SIZE]; - - demodulation_signal - .iter_mut() - .zip(adc_samples.iter()) - .enumerate() - .for_each(|(i, (s, sample))| { - let sample_phase = (self.harmonic.wrapping_mul( - (demodulation_frequency.wrapping_mul(i as u32)) - .wrapping_add(demodulation_initial_phase), - )) - .wrapping_add(self.phase_offset); - let (cos, sin) = cossin(sample_phase as i32); - // cos/sin take up 32 bits and sample takes up 16 - // bits. Make this fit into a 32 bit result. - s.0 = ((*sample as i64 * cos as i64) >> 16) as i32; - s.1 = ((*sample as i64 * sin as i64) >> 16) as i32; - }); - - if self.batch_index < ADC_BATCHES as u32 - 1 { - self.batch_index += 1; - } else { - self.batch_index = 0; - self.last_phase = Some(self.last_phase.unwrap() - (1 << 32)); - } - - Ok(demodulation_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]) { - 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 `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; SAMPLE_BUFFER_SIZE], -) -> [Complex; DECIMATED_BUFFER_SIZE] { - let n_k = SAMPLE_BUFFER_SIZE / DECIMATED_BUFFER_SIZE; - debug_assert!(SAMPLE_BUFFER_SIZE == DECIMATED_BUFFER_SIZE || n_k % 2 == 0); - - let mut signal_decimated = [(0_i32, 0_i32); 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 for which the magnitude and phase -/// 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]) { - signal.iter_mut().for_each(|s| { - let new_i = [s.0, s.1].iter().map(|i| i * i).sum(); - let new_q = atan2(s.1, s.0); - s.0 = new_i; - s.1 = new_q; - }); -} - -#[cfg(test)] -mod tests { - use super::*; - - #[test] - /// Ensure that the demodulation signals are within some tolerance - /// band of the target value given the phase and frequency values - /// provided by the PLL. - fn demodulate() { - const PLL_SHIFT_FREQUENCY: u8 = 4; - const PLL_SHIFT_PHASE: u8 = 3; - const HARMONIC: u32 = 1; - const PHASE_OFFSET: u32 = 0; - let mut lockin = Lockin::new( - HARMONIC, - PHASE_OFFSET, - IIR { ba: [0; 5] }, - PLL_SHIFT_FREQUENCY, - PLL_SHIFT_PHASE, - ); - - // Duplicate the PLL outside demodulate so that we don't test - // its behavior. - let mut tracking_pll = PLL::default(); - let mut tracking_phase: i32 = 0; - let mut tracking_frequency: i32 = 0; - - const REFERENCE_FREQUENCY: usize = 10_000; - let mut reference_edge: usize = REFERENCE_FREQUENCY; - - // Ensure that we receive at most 1 timestamp per batch. - debug_assert!( - REFERENCE_FREQUENCY >= SAMPLE_BUFFER_SIZE * ADC_SAMPLE_TICKS - ); - - for batch in 0..100_000 { - let tick: usize = batch * ADC_SAMPLE_TICKS * SAMPLE_BUFFER_SIZE; - let timestamp: Option; - - // When the reference edge occurred during the current - // batch acquisition, register the timestamp and update - // the tracking PLL. - if reference_edge >= tick - && reference_edge < tick + ADC_SAMPLE_TICKS * SAMPLE_BUFFER_SIZE - { - timestamp = Some(reference_edge as u32); - - let tracking_update = tracking_pll.update( - reference_edge as i32, - PLL_SHIFT_FREQUENCY, - PLL_SHIFT_PHASE, - ); - 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(_) => {} - } - } - } -} diff --git a/src/main.rs b/src/main.rs index c5a9698..538fb2a 100644 --- a/src/main.rs +++ b/src/main.rs @@ -60,14 +60,34 @@ use heapless::{consts::*, String}; // The number of ticks in the ADC sampling timer. The timer runs at 100MHz, so the step size is // equal to 10ns per tick. // Currently, the sample rate is equal to: Fsample = 100/256 MHz = 390.625 KHz -const ADC_SAMPLE_TICKS: u16 = 256; +const ADC_SAMPLE_TICKS_LOG2: u16 = 8; +const ADC_SAMPLE_TICKS: u16 = 1 << ADC_SAMPLE_TICKS_LOG2; // The desired ADC sample processing buffer size. -const SAMPLE_BUFFER_SIZE: usize = 8; +const SAMPLE_BUFFER_SIZE_LOG2: usize = 3; +const SAMPLE_BUFFER_SIZE: usize = 1 << SAMPLE_BUFFER_SIZE_LOG2; + +// The number of ADC batches in one timer overflow period. +pub const ADC_BATCHES_LOG2: usize = + 32 - SAMPLE_BUFFER_SIZE_LOG2 - ADC_SAMPLE_TICKS_LOG2 as usize; +pub const ADC_BATCHES: usize = 1 << ADC_BATCHES_LOG2; // The number of cascaded IIR biquads per channel. Select 1 or 2! const IIR_CASCADE_LENGTH: usize = 1; +// TODO should these be global consts? + +// Frequency scaling factor for lock-in harmonic demodulation. +const HARMONIC: u32 = 1; + +// Phase offset applied to the lock-in demodulation signal. +const PHASE_OFFSET: u32 = 0; + +// The PLL locks to an external reference signal. See `PLL::update()` +// for a description of shift_frequency and shift_phase. +const PLL_SHIFT_FREQUENCY: u8 = 4; +const PLL_SHIFT_PHASE: u8 = 3; + #[link_section = ".sram3.eth"] static mut DES_RING: ethernet::DesRing = ethernet::DesRing::new(); @@ -84,7 +104,11 @@ mod timers; use adc::{Adc0Input, Adc1Input}; use dac::{Dac0Output, Dac1Output}; -use dsp::iir; +use dsp::{ + divide_round, iir, iir_int, + pll::PLL, + trig::{atan2, cossin}, +}; use pounder::DdsOutput; #[cfg(not(feature = "semihosting"))] @@ -205,6 +229,12 @@ const APP: () = { adcs: (Adc0Input, Adc1Input), dacs: (Dac0Output, Dac1Output), input_stamper: digital_input_stamper::InputStamper, + pll: PLL, + batch_index: u32, + reference_phase: i64, + reference_frequency: i64, + iir_lockin: iir_int::IIR, + iir_state_lockin: [iir_int::IIRState; 2], eeprom_i2c: hal::i2c::I2c, @@ -927,6 +957,13 @@ const APP: () = { #[cfg(not(feature = "pounder_v1_1"))] let pounder_stamper = None; + let pll = PLL::default(); + let batch_index = 0; + let reference_phase = 0; + let reference_frequency = 0; + let iir_lockin = iir_int::IIR { ba: [0; 5] }; + let iir_state_lockin = [[0; 5]; 2]; + // Start sampling ADCs. sampling_timer.start(); timestamp_timer.start(); @@ -942,6 +979,13 @@ const APP: () = { pounder: pounder_devices, pounder_stamper, + pll, + batch_index, + reference_phase, + reference_frequency, + iir_lockin, + iir_state_lockin, + eeprom_i2c, net_interface: network_interface, eth_mac, @@ -949,7 +993,7 @@ const APP: () = { } } - #[task(binds=DMA1_STR4, resources=[pounder_stamper, adcs, dacs, iir_state, iir_ch, dds_output, input_stamper], priority=2)] + #[task(binds=DMA1_STR4, resources=[pounder_stamper, adcs, dacs, iir_state, iir_ch, dds_output, input_stamper, pll, iir_lockin, iir_state_lockin, batch_index, reference_phase, reference_frequency], priority=2)] fn process(c: process::Context) { if let Some(stamper) = c.resources.pounder_stamper { let pounder_timestamps = stamper.acquire_buffer(); @@ -965,24 +1009,71 @@ const APP: () = { c.resources.dacs.1.acquire_buffer(), ]; - let _timestamp = c.resources.input_stamper.latest_timestamp(); + let reference_phase = c.resources.reference_phase; + let reference_frequency = c.resources.reference_frequency; - 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 mut y = 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::() }; - // Convert to DAC code - dac_samples[channel][sample] = y as u16 ^ 0x8000; - } + if let Some(t) = c.resources.input_stamper.latest_timestamp() { + let res = c.resources.pll.update( + t as i32, + PLL_SHIFT_FREQUENCY, + PLL_SHIFT_PHASE, + ); + *reference_phase = res.0 as u32 as i64; + *reference_frequency = res.1 as u32 as i64; } + let demodulation_frequency = divide_round( + 1 << (64 - SAMPLE_BUFFER_SIZE_LOG2 - ADC_BATCHES_LOG2), + *reference_frequency, + ) as u32; + let batch_index = c.resources.batch_index; + let demodulation_initial_phase = divide_round( + (((*batch_index as i64) << (32 - ADC_BATCHES_LOG2)) + - *reference_phase) + << 32, + *reference_frequency, + ) as u32; + + if *batch_index < ADC_BATCHES as u32 - 1 { + *batch_index += 1; + } else { + *batch_index = 0; + *reference_phase -= 1 << 32; + } + + let [dac0, dac1] = dac_samples; + let iir_lockin = c.resources.iir_lockin; + let iir_state_lockin = c.resources.iir_state_lockin; + + dac0.iter_mut().zip(dac1.iter_mut()).enumerate().for_each( + |(i, (d0, d1))| { + let sample_phase = (HARMONIC.wrapping_mul( + (demodulation_frequency.wrapping_mul(i as u32)) + .wrapping_add(demodulation_initial_phase), + )) + .wrapping_add(PHASE_OFFSET); + let (cos, sin) = cossin(sample_phase as i32); + + let mut signal = (0_i32, 0_i32); + + signal.0 = ((adc_samples[0][i] as i16 as i64 * cos as i64) + >> 16) as i32; + signal.1 = ((adc_samples[0][i] as i16 as i64 * sin as i64) + >> 16) as i32; + + signal.0 = + iir_lockin.update(&mut iir_state_lockin[0], signal.0); + signal.1 = + iir_lockin.update(&mut iir_state_lockin[1], signal.0); + + let magnitude = signal.0 * signal.0 + signal.1 * signal.1; + let phase = atan2(signal.0, signal.0); + + *d0 = (magnitude >> 16) as i16 as u16; + *d1 = (phase >> 16) as i16 as u16; + }, + ); + if let Some(dds_output) = c.resources.dds_output { let builder = dds_output.builder().update_channels( &[pounder::Channel::Out0.into()], @@ -994,7 +1085,6 @@ const APP: () = { builder.write_profile(); } - let [dac0, dac1] = dac_samples; c.resources.dacs.0.release_buffer(dac0); c.resources.dacs.1.release_buffer(dac1); }