dealing with life's imperfections (Part 6/12) 🕊️

At one of my old apartments, our (old, at the time) router sat under an (old, at the time) microwave. I remember this arrangement fondly, filled with nostalga – plus wifi was always solid in the kitchen, which did matter at the time. I remember many a dinner watching compiler output as I streamed Netflix on my Nexus 7. Those were nice days. Only hiccup was, from time to time, I’d get hit with a spinning loading bar right when the show was getting good. Eventually we realized it was whenever anyone flipped on the microwave – no doubt our delicious kraft mac cup’s boiling water was jamming the AP’s ability to hear anything.

Turns out the world is pretty noisy. It also turns out the RF world is possibly even more noisy – in no small part due to the lack of physical stimulus when something is noisy in RF. While an ear-piercing whine coming from your fridge will bother you until you replace it – there are very few things in the RF spectrum that cause you the same amount of physical pain. Humans tend to be pretty bad at fixing stuff that’s “out of sight, out of mind”. The nice folks over at the FCC have done an incredible amount of work keeping the bands orderly (emissions testing, specific bands for ISM equipment like microwaves), but nevertheless, things break, manufacturing swaps an IC out and causes a new unintentional emission, or QA doesn’t catch a bad unit. It tends to go largely unnoticed.

That is, except when its your microwave jamming your router. In that case, you do what any normal, well-adjusted and safety conscious person does – we used enough tin foil between the router and microwave to stop enough of the magnetron’s output from jamming the router. Fixed. This isn’t even a one-off, at a different (and, much more well maintained) apartment, my fancy LED-powered track lights would completely jam the 144MHz (2m) ham radio band when turned on (and I have the waterfalls and videos to prove it). The solution there was to (you guessed it), just leave the lights off if I wanted to play radios; fixing the problem for real was just plain too involved. Doubly so as a renter. Ham radio by candle light.

And here we are – full circle. We had our chance to reap, and now comes the time to sow.

When a signal moves through freespace from the transmitter to a receiver, many unspeakable – and unknowable – horrors have interacted with that signal in a number of (usually) unpleasant ways. We have noises produced by human sources (like other users of the same radio service, or a microwave), noises produced by non-human sources (calm down alien people – things like lightning or, say, the sun), or just plain physical interactions (attenuation of part of the signal’s frequency but not others, phase shifts across the bandwidth, movement between the transmitter and receiver, clock differences in crystals, heat in the SDR adding RF noise, grounding issues, etc).

Let’s attempt to correct some of those issues before we demodulate any data. Attempting to demodulate uncorrected bursts has never worked for me (yet?) – making this a particularly important bug-fixing (and bug-creating) section of the codebase. Problems in signal cleanup prior to demodulation can trivially cascade throughout the stack, rendering perfectly fine bursts undecodable.

Frequency Offset Calculations

Let’s start off with the easiest one to reason about – two radios, both tuned to the exact same frequency, will still be slightly off in frequency from each other. This means if you transmit exactly on frequency according to one radio, the receiving SDR may see it a few kHz off center. Let’s go ahead and account for that drift.

There are two types of frequency offset calculations we’ll do. First is computing the coarse frequency offset, to correct errors on the order of kHz, and the second is computing the fine frequency offset, to adjust the frequency on the order of 10s to 100s of Hz.

The first core concept here is that we need to remember that a consistent change in phase over time signifies a frequency offset.

This means that if we measure the phase difference between samples that should notionally be the same, and average the phase difference between them, we can determine how much we’d expect the frequency to be off.

/// Guess at the center frequency offset between our rx
/// channel and the transmitter's tx channel.
///
/// The provided length (`len`) is the length of every
/// element within `seq` -- the number of IQ samples
/// we're comparing in each chunk. The rate (`rate`) is
/// the sampling rate, in Hz, and the sequence (`seq`)
/// are blocks of repeating IQ.
///
/// In the case of the Schmidl-Cox preamble, the `len`
/// would be something like the order (16), and the
/// `seq` would contain `order * count` entries
/// (of `len` each), at the mode's `rate` (2.5MHz/2.5e6)
pub fn frequency_offset(
  len: usize,
  rate: f64,
  seq: &[&[IQ]],
) -> f64 {
  // Total number of IQ samples (also, therefore, time)
  // that `seq` spans.
  let nsamples = (seq.len() - 1) * len;

  // For each pair of `seq` slices, multiply each IQ sample
  // by the complex conjugate of the repeated IQ sample,
  // "rewinding" it to determine if there is any "drift".
  //
  // If there is no frequency shift at all, doing this will
  // result in `1+0j` -- a phase of 0.

  let mut phase = 0.0f64;
  for [last, now] in seq.array_windows::<2>() {
    for (iq1, iq2) in last.iter().zip(now.iter()) {
      phase += (*iq1 * iq2.conj()).phase();
    }
  }

  // since we added up the phase offsets, we may now divide
  // it (to compute the average), and adjust the phase by
  // the sample rate to get the estimated difference in
  // frequency (in hz) based on the systemic change in phase.
  //
  // remember: changing phase continuously over time is
  // just a change in frequency.

  phase /= nsamples as f64;
  phase /= len as f64;
  let step = 1.0 / rate;
  phase / (TAU * step)
}

We can use this technique anywhere IQ samples repeat in the time-domain; within pigeon we use this technique three different ways to correct systemic frequency offsets.

Coarse Frequency Offset

If we look back to the construction of our Preamble, we’ll be using the smallest repeating sequence first – the bins of 8 IQ samples to compute the CFO. Using the approach above, what this means is that we’ll be iterating over each sample (except the last 8), and compare the phase to the next, closest, repeating sample.

At the end of this, we have the accumulated phase error over len(iq) - 8 samples. We can then use the estimated frequency offset to frequency shift all the IQ samples, and fix any remaining errors by doing a fine frequency offset calculation.

Fine Frequency Offset

Just as above, we’ll be using the Preamble’s structure to compare repeating samples’ phase to do our FFO calculations, but this time, we’ll be doing it between each of the two iFFT generated IQ blocks – 64 samples apart, not 8.

Coarse frequency alignment will check for phase differences between samples very close to each other since there’s an assumed large frequency difference, which will represent as relatively large phase changes. Fine frequency alignment assumes very small frequency changes, so we want to look at changes over a much longer stretch of time.

Shifting Frequency

Given a frequency offset, the next (very logical) question is how do we shift our IQ values in frequency? Doing this in the frequency domain is conceptually possible (and maybe even possible for some specific uses), but not with OFDM, unfortunately. With OFDM, alignment is critical because taking an FFT that overlaps multiple symbols will cause corruption. As a result, we can shift the frequency in the time domain by mixing our IQ stream with a carrier wave representing our desired shift by doing complex multiplication between the IQ stream and our generated carrier wave.

/// shift time-series IQ values (`dst`) in frequency by
/// mixing our signal with a carrier wave at `freq` Hz,
/// at some sampling rate (`rate_f`).
fn shift(
  rate_f: f64,
  freq: f64,
  dst: &mut [IQ]
) {
    let rate = rate_f as usize;
    for (idx, sample) in dst.iter_mut().enumerate() {
        // for each time step (`ts`) within a second (there
        // are `rate` number of samples per second), generate
        // a carrier wave at `freq` hz.
        let ts = (idx % rate) as f64;

        // Turn the time step (`ts`) into a value between
        // 0 and 1 (0 being the exact first sample at
        // the start of a second and 1 being the exact last
        // sample of a second), multiply by the desired
        // frequency (`freq`), and convert from 0-1 to
        // radians (steps around the circle) as an input
        // to sin/cosine.
        let step = (ts / rate_f) * TAU * freq;

        // generate our carrier wave as an IQ sample,
        // and multiply (mix) the signal in `dst` by
        // it, shifting it in frequency.
        *sample = *sample * iq!(
          math::cos(step) as f32,
          math::sin(step) as f32
        );
    }
}

Old radio hands may recognize this as a frequency mixer (or just “mixer”). We have our input signal, some carrier wave, and mixing them causes a shift in frequency by the cairrer – nothing too crazy here!

Channel Estimation

In addition to changes in frequency, the “channel” (free-space path that the signal(s!) took from the transmitter to the receiver) also alters the signal, usually in frequency specific ways – the signal is altered differently at different places in the frequency domain! This can look like the phase slowly shifting from one end of the frequency domain to the other, or even frequency specific attenuation.

Thankfully, we have known values in our signal – the pilot subcarriers. We can process each subcarrier, and correct phase and amplitude imbalances in the frequency domain by interpolating those offsets. The interpolate function takes the form of a function which, technically, does both interpolation (filling in values between two known values) as well as extrapolation (filling in values “past” a known value based on prior data).

I’ve opted to use linear interpolation (since it’s simple to implement and I have no evidence I should be using anything else, and I get to say I’m “lerping”). At each step, we’ll go through our array at each index (we’ll call this index into the buffer x), with two known values (the “leftmost” known bin, which we’ll call x0, and its value in the buffer, which we’ll call y0, and the next known bin to its “right”, which we’ll call x1, and its value in the buffer, y1). We’re solving for y, which is what the interpolated value at our current offset (x) should be.

/// Interpolate (or extrapolate) values for `dst`,
/// given the known values within the `dst` buffer
/// at the indexes provided in `known`.
pub fn interpolate(dst: &mut [IQ], known: &[usize]) {
    let mut known = known.iter();
    let mut left = *known.next().unwrap();
    let mut right = *known.next().unwrap();

    // here we're going to iterate through our values,
    // and keep track of, relative to our current
    // index, where the first known value to our
    // "left" is, and where the first known value
    // to our "right" is.
    //
    // If there are no known values to the "left"
    // or "right" of our index, they will be the
    // be the closest 2 "leftmost" or "rightmost"
    // known values, for extrapolation rather than
    // interpolation.

    for idx in 0..dst.len() {
        if idx > right
            && let Some(n) = known.next()
        {
            // We're past the "right" side,
            // so let's shift forward.
            left = right;
            right = *n;
        }

        if idx == left || idx == right {
            // if we're at a known value we don't want
            // to change anything.
            continue;
        }

        let x0 = left as f32;
        let y0 = dst[left];
        let x1 = right as f32;
        let y1 = dst[right];
        let x = idx as f32;

        ..
    }
}

Now, to actually do the interpolation, we can use see how far along we are between x0 and x1, and scale our y value evenly at each step, according to where we want to start and end.

Easy peasy. Now, what we’re going to want to do is lerp (linear interpolate), independently, the difference in phase between the known pilot phases, as well as the magnitude between the known pilot magnitudes.

Let’s start with magnitude, since it’s a bit easier. Here we can fill in the blanks in the equation above with our abs operator, which will return the absolute value (magnitude) of our interpolated complex number correction (the y0 and y1 value).

Which I wound up implementing (unsupprisingly) as:

let mag = y0.abs()
    + (x - x0)
    * (y1.abs() - y0.abs())
    / (x1 - x0);

And, now, let’s do the same, but with phase. Only trick here is that we’re interpolating values between -PI and PI, rather than from 0 to 1. This poses a bit of a pickle, since we need to be sure that we’re going in the “right” direction. This can only really work right if we have something that’s radian-aware that can give us the shortest distance between two radian values.

/// Return the minimum distance in modulo-space between
/// two Radian values.
fn delta(left: f32, right: f32) -> f32 {
    let clockwise = right - left;
    let counterclockwise = TAU - clockwise.abs();
    if clockwise.abs() < counterclockwise.abs() {
        clockwise
    } else {
        if clockwise.is_negative() {
            counterclockwise
        } else {
            -counterclockwise
        }
    }
}

I figure just showing code is a bit easier to understand than trying to write all the math out – it just kinda confuses things. Now, let’s go ahead and use our new delta function to interpolate our new phase value, just like we did for mag:

This translates pretty much exactly like you’d expect:

let phase = {
    let y0 = y0.phase();
    let y1 = y1.phase();
    let y_d = -y1.delta(y0);
    *y0 + (x - x0) * y_d / (x1 - x0)
};

Now, since we have our phase and mag, we can create our new interpolated complex number, since those are “just” polar coordinates, and assign that back into our buffer:

dst[idx] = IQ::polar(phase, mag);

This leaves us with complex numbers which have been corrected for the phase and magnitude distortions based on our pilot subcarriers.

Let’s talk about forward error correction →