kai
/
pounder_test
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

move lock-in code to main.rs

master
Matt Huszagh 2021-01-12 10:45:34 -08:00
parent 891aad3f17
commit e14aa8b613
3 changed files with 110 additions and 413 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
dsp/src/lib.rs
View File

@ -116,7 +116,6 @@ where
pub mod iir;
pub mod iir_int;
pub mod lockin;
pub mod pll;
pub mod trig;
pub mod unwrap;

392
dsp/src/lockin.rs
View File

@ -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<i64>,
last_frequency: Option<i64>,
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<u32>,
) -> Result<[Complex<i32>; 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<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);
});
}
}
/// 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<i32>; SAMPLE_BUFFER_SIZE],
) -> [Complex<i32>; 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<i32>]) {
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<u32>;
// 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(_) => {}
}
}
}
}

130
src/main.rs
View File

@ -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<hal::stm32::I2C2>,
@ -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::<i16>() };
// 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);
}