Add code for using the USRP N210 receiver

[varenius]

We need code to actually measure stuff, by getting data from the USRP N210 receiver. Some test-code was written but not put on github recently by @hlinander and @symmsaur , I paste it here so we don't forget it:

use std::env::set_var;

use anyhow::{Context, Result};
use num_complex::Complex;
use tap::Pipe;
use uhd::{self, StreamCommand, StreamCommandType, StreamTime, TuneRequest, Usrp};
use rustfft::{FftPlanner, FftDirection, algorithm::Radix4};
use rustfft::Fft;
use plotters::prelude::*;

const CHANNEL: usize = 0;
const NUM_SAMPLES: usize = 1e6 as usize;

pub fn main() -> Result<()> {
    set_var("RUST_LOG", "DEBUG");
    env_logger::init();

    log::info!("Starting receive test");

    // let mut usrp = Usrp::find("")
    //     .context("Failed to open device list")?
    //     .drain(..)
    //     .next()
    //     .context("Failed to find a valid USRP to attach to")?
    //     .pipe(|addr| Usrp::open(&addr))
    //     .context("Failed to find properly open the USRP")?;

    let mut usrp = Usrp::open("addr=192.168.5.31").unwrap();

    usrp.set_rx_sample_rate(1e6, CHANNEL)?;
    usrp.set_rx_antenna("TX/RX", CHANNEL)?;
    usrp.set_rx_frequency(&TuneRequest::with_frequency(1.42044057517862e9), CHANNEL)?;

    let mut receiver = usrp
        .get_rx_stream(&uhd::StreamArgs::<Complex<i16>>::new("sc16"))
        .unwrap();

    let mut buffer = uhd::alloc_boxed_slice::<Complex<i16>, NUM_SAMPLES>();

    receiver.send_command(&StreamCommand {
        command_type: StreamCommandType::CountAndDone(buffer.len() as u64),
        time: StreamTime::Now,
    })?;
    let status = receiver.receive_simple(buffer.as_mut())?;

    log::info!("{:?}", status);
    log::info!("{:?}", &buffer[..16]);

    let mut planner = FftPlanner::<f32>::new();

    let fft_samples: usize = 4096;
    let mut fft_buffer: Vec<num_complex::Complex<f32>> = buffer[0..fft_samples].iter().copied()
        .map(|x| {
            Complex::<f32>::new(x.re as f32, x.im as f32)
        }).collect();
    let fft = Radix4::new(fft_samples, FftDirection::Forward);
    fft.process(&mut fft_buffer);

    let backend = BitMapBackend::new("plot.png", (800, 600)).into_drawing_area();

    //let freqs = (0..fft_samples).map(|x| NUM_SAMPLES)

    let chart = ChartBuilder::on(&backend)
                    .build_cartesian_2d(0..fft_samples, -1.0...1.0).unwrap();

    chart.draw_series(
        LineSeries::new(
            (fft_buffer.iter().enumerate().map(|(idx, x)| (idx, x.re))),
            &BLACK)
    ).unwrap();
    //buffer[0..4096]

    Ok(())
}

[varenius]

Will need to set gain as well. The USRP c++ API has a method to set gain without specified gain element name at https://github.com/EttusResearch/uhd/blob/master/host/include/uhd/usrp/multi_usrp.hpp#L1073 this is not available in the uhd crate; instead we need to specify a gain name at https://github.com/samcrow/uhd-rust/blob/master/uhd/src/usrp.rs#L675 However, the cpp function above just uses the gain name ALL_GAINS and then calls the other function witht this. According to https://github.com/EttusResearch/uhd/blob/master/host/lib/usrp/multi_usrp.cpp#L49 the ALL_GAINS variable is set to empty string, so possibly it would work just calling the uhd crate set_rx_gain with name = "".

[varenius]

We also probably want to use the function https://github.com/samcrow/uhd-rust/blob/1ec3f512568b9d0180b1c59107cdbf52506b068c/uhd/src/usrp.rs#L643 to get rid of DC spike in center (has been a problem in the past with the old system)

[varenius]

Good documentation on I/Q sampling and SDR applications: https://pysdr.org/content/sampling.html#quadrature-sampling

[hlinander]

Thanks! Started preparing a PR for this but since the UHD crate requires system libraries I didn't want to add this straight into the current dependencies. Was planning to start with a separate workspace but then I got side-tracked.

[varenius]

Tried to evolve the code a little. Mixed results, mainly because I don't know what I'm doing ;). But, what I tried to do is

• Set gain for USRP
• Loop through all samples and do fft for chunks of FFT_POINTS each.
• Stack the resulting power spectrum = abs(fft)
• Plot the result

The code below compiles, but does not produce a plot. It writes out the 256-len vector, which I put into python to make a plot, which looks very wrong :). Anyway, I just paste it here for now:

salsa_dev@brage:~/uhd_test/src$ cat main.rs 
use std::env::set_var;
use anyhow::{Result};
use uhd::{self, StreamCommand, StreamCommandType, StreamTime, TuneRequest, Usrp};
use rustfft::{Fft, FftDirection, num_complex::Complex, algorithm::Radix4};
use plotters::prelude::*;

const CHANNEL: usize = 0; // USRP input channel
const TINT: usize = 10; // integration time, seconds
const SAMP_RATE: f64 = 5e6 ; // sample rate, Hz
const NUM_SAMPLES: usize = TINT * SAMP_RATE as usize; // total number of samples expected
const FFT_POINTS: usize = 256; // Number of points in FFT, setting spectral resolution
const NUM_STACK: usize = NUM_SAMPLES % FFT_POINTS;

pub fn main() -> Result<()> {
    set_var("RUST_LOG", "DEBUG");
    env_logger::init();

    log::info!("Starting receive test");

    let mut usrp = Usrp::open("addr=192.168.5.31").unwrap();

    usrp.set_rx_sample_rate(SAMP_RATE, CHANNEL)?;
    usrp.set_rx_gain(40.0, CHANNEL, "")?; // empty string should mean all gains
    usrp.set_rx_antenna("TX/RX", CHANNEL)?;
    usrp.set_rx_frequency(&TuneRequest::with_frequency(1.42044057517862e9), CHANNEL)?;
    usrp.set_rx_dc_offset_enabled(true, CHANNEL)?;

    //let fr = usrp.get_rx_frequency_range(CHANNEL).unwrap();
    //print!("{:?}",fr);
    //
    //let gains = usrp.get_rx_gain_names(CHANNEL).unwrap();
    //print!("{:?}",gains);
    //
    //let gain = usrp.get_rx_gain(CHANNEL, "").unwrap();
    //print!("{:?}",gain);

    let mut receiver = usrp
        .get_rx_stream(&uhd::StreamArgs::<Complex<i16>>::new("sc16"))
        .unwrap();

    let mut buffer = uhd::alloc_boxed_slice::<Complex<i16>, NUM_SAMPLES>();

    receiver.send_command(&StreamCommand {
        command_type: StreamCommandType::CountAndDone(buffer.len() as u64),
        time: StreamTime::Now,
    })?;
    let status = receiver.receive_simple(buffer.as_mut())?;

    log::info!("{:?}", status);
    log::info!("{:?}", &buffer[..16]);

    // setup fft
    let fft = Radix4::new(FFT_POINTS, FftDirection::Forward);
    // array to store power spectrum (abs of FFT result)
    let mut fft_abs: [f64; FFT_POINTS] = [0.0; FFT_POINTS];
    // Loop through the samples, taking FFF_POINTS each time
    for n in 0..NUM_STACK {
        let mut fft_buffer: Vec<Complex<f64>> = buffer[n*FFT_POINTS..(n+1)*FFT_POINTS].iter().copied()
            .map(|x| {
                Complex::<f64>::new(x.re as f64, x.im as f64)
            }).collect();
        // Do the FFT
        fft.process(&mut fft_buffer);
        // Add absolute values to stacked spectrum
        for i in 0..FFT_POINTS {
            fft_abs[i] = fft_abs[i] + fft_buffer[i].norm();
        }
    }
    // Normalise spectrum by number of stackings
    for i in 0..FFT_POINTS {
        fft_abs[i] = fft_abs[i] / ( NUM_STACK as f64);
    }

    log::info!("{:?}", &fft_abs);

    let backend = BitMapBackend::new("plot.png", (800, 600)).into_drawing_area();

    //let freqs = (0..fft_samples).map(|x| NUM_SAMPLES)

    let mut chart = ChartBuilder::on(&backend)
                    .build_cartesian_2d(0..FFT_POINTS, -1.0..1.0).unwrap();

    chart.draw_series(LineSeries::new(fft_abs.iter().copied().enumerate(),RED)).unwrap();

    Ok(())
}

[varenius]

Tried measuring on GNSS satellites, much stronger than HI, for testing. Realised that the FFT output is flipped; now the edges were in the middle. Seems empirically fixed by doing

        // Add absolute values to stacked spectrum
        // Seems the pos/neg halves of spectrum are flipped, so reflip them
        for i in 0..FFT_POINTS/2 {
            fft_abs[i+FFT_POINTS/2] = fft_abs[i+FFT_POINTS/2] + fft_buffer[i].norm();
            fft_abs[i] = fft_abs[i] + fft_buffer[i+FFT_POINTS/2].norm();
        }

Took a spectrum with same bandwidth, int time and frequency towards a GNSS satellite to compare this script (left panel) with current software (right panel):

(plot in python, still didn't get the rust plot to work) Similar, but not quite as smooth. I wonder if the difference is due to FFT windowing or similar, which could bring down the noise?

[varenius]

One thing wrong was that the power spectrum is abs(x)**2, but I had abs(x). Corrected via:

    // Normalise spectrum by number of stackings
    // and do **2 to get power spectrum
    for i in 0..FFT_POINTS {
        fft_abs[i] = fft_abs[i]*fft_abs[i] / ( NUM_STACK as f64);
    }

With this, the GNSS signal looks very similar comparing new (left) with old (right):

However, looking carefully, there is more noise on the new plot. This is much more clear if we try to do the same comparison but for an HI spectrum:

This needs to be understood and fixed. Ideas for now are:

• We have an IIR filter in the gnuradio processing at https://github.com/varenius/salsa/blob/main/Control_program/receiver.py#L62. This is not used for rust (yet).
• The gnuradio flow has a BlackmanHarris window (https://github.com/varenius/salsa/blob/main/Control_program/receiver.py#L58) for the FFT. Not sure if the current use of rustfft has any (default) window?
• These two panels both show 128 points. However, the old SALSA software does the FFT on 4096 points, stacks these, and then averages down afterwards, while the current rustfft is doing the FFT at 128 pts directly. This may be the first thing to check!

[varenius]

AHA! It was the FFTsize, here is with left (rustfft) doing 4096 points:

Now it's just to make rust filter RFI and average like the old software, and we should be done with the proof-of-concept, and can integrate it in the full software!

[hlinander]

Nice! :partying_face:

[varenius]

Finally got a plot from rust! (Spent 2h wondering why the plot was all white, before I realised it was because all points were above the plotted range...):

Still no median filtering or averaging, will have to think more about that.

[varenius]

For median filtering we can use https://github.com/regexident/median with window size 21 on 4096 sampled FFT. Then simple average via for loop gives, finally, apples to apples:

The only problem now is that the FFT step (doing one FFT per 4096 elements of all the sampled data) takes almost as long as the integration time. Using "planner" instead of directly chosing algorithm improves this a little (likely due to enabled accelerations) but not quite there yet. In the old code there is an IIR averaging filter which decimates by factor 10 (introduced by student many years ago, see https://github.com/varenius/salsa/blob/main/Control_program/receiver.py#L62). I'm not 100% sure how this works, but maybe something similar could be implemented to get real-time performance?

EDIT: Maybe one of https://crates.io/keywords/iir ?

Current code snippet for filtering and averaging

    let mut ymax : f64 = 0.0;
    let mut ymin : f64 = 0.0;
    // Normalise spectrum by number of stackings,
    // do **2 to get power spectrum, and median filter
    // also lot max/min for plotting
    let mut filter = Filter::new(21);
    for i in 0..FFT_POINTS {
        fft_abs[i] = fft_abs[i]*fft_abs[i] / ( NUM_STACK as f64);
        fft_abs[i] = filter.consume(fft_abs[i]);

        if fft_abs[i] > ymax {
            ymax = fft_abs[i];
        }
        if fft_abs[i] < ymin {
            ymin = fft_abs[i];
        }
    }

    // Average spectrum to save data
    let mut fft_avg: [f64; AVG_POINTS] = [0.0; AVG_POINTS];
    for i in 0..AVG_POINTS {
        let mut avg = 0.0;
            for j in NUM_AVG*i..NUM_AVG*(i+1) {
                avg = avg+fft_abs[j];
            }
        fft_avg[i] = avg/(NUM_AVG as f64);
    }

[varenius]

Thinking a bit more: Maybe the higher processing speed for SALSA is also related to the fact that the FFT works on the stream, rather than the recorded buffer. In current rust implementation, we first record all (e.g. 10 seconds) data, THEN process with the FFT. It should be possible to FFT as soon as each chunk (4096 elements) arrive from the USRP. In this way, we can do the FFT while sampling, which would mean little or no extra waiting time after obs has finished?

[varenius]

I'm more and more convinced we want to do FFTs of incoming samples in paralell to the sampling happening. Not sure how to do this (neither in general, nor in rust specifically). But I note there are different streaming modes available for the USRP, see https://github.com/samcrow/uhd-rust/blob/1ec3f512568b9d0180b1c59107cdbf52506b068c/uhd/src/stream/mod.rs#L150 which is a wapper for https://github.com/EttusResearch/uhd/blob/master/host/include/uhd/usrp/usrp.h#L62.

Edit: futuresdr seems to have solved this, maybe we could borrow from them? See e.g. https://github.com/FutureSDR/FutureSDR/blob/main/src/blocks/fft.rs and examples in the repo.

[varenius]

I just remembered (new to rust) that we can compile both with "cargo run" and "cargo build --release", the latter with more optimisation. I did a test where I compiled uhd_test/src/main.rs in both ways and observed + fft for 30 seconds. With "run", the fft of the 30 sec recording takes 21 seconds. With "release" the FFT takes less than 1 second! So, in principle we don't need to process this stream-wise as long as we compile with "--release".

[varenius]

Re-arranged the code a bit for simple frequency-switched measurements, works nicely! Comparison:

Tentatively the result with rust is mayhaps a little better/clearer. Possibly due to different median filtering schemes. Still, good stuff! Next step is to integrate this into the full code, to get the plot in the web interface etc.

[varenius]

Closing this in favour of #34.