pluto-sdr 0.1.1

// Author: Roman Hayn
// MIT License 2023

use std::f32::consts::PI;

/// Create a FIR Filter, which can be applied to signals (simply `Vec<f32>`).
/// ```
/// use pluto_sdr::filter::Filter;
/// // arbitrary signal
/// let signal = vec![1., 4., -3., 5., 10., -19., -4., 10., -2., -1.];
/// // coeff describe a 4 point moving average
/// let lpf_coefficients = vec![0.25, 0.25, 0.25, 0.25];
/// // create the low pass filter
/// let lpf = Filter::new(lpf_coefficients.into());
///
/// // filter the signal
/// let filtered = lpf.filter_windowed(&signal);
/// println!("{:?}", filtered);
/// ```
#[derive(PartialEq, Debug)]
pub struct Filter {
    coeff: Vec<f32>,
}

impl Filter {
    /// Create a new FIR Filter.
    /// `coeff` is of type `Vec<f32>` and contains the FIR Filters Coefficients: `[b_0, b_1, b_2, ... ]`.
    pub fn new(coeff: Vec<f32>) -> Self {
        Filter { coeff }
    }

    /// The Filter size is just the size of its coefficients
    pub fn len(&self) -> usize {
        self.coeff.len()
    }

    /// recall that windowed convolution produces a signal of length `N - M + 1`
    /// where `N` is the length of signal 1: `signal`;
    ///   and `M` is the length of signal 2: `self.coeff`
    /// The `signal` should be longer than the internal filter coefficients.
    pub fn filter_windowed(&self, signal: &Vec<f32>) -> Vec<f32> {
        signal
            .windows(self.coeff.len())
            .map(|window| window.iter().zip(&self.coeff).map(|(x, y)| x * y).sum())
            .collect()
    }

    /// Creates a moving average filter, where all coefficients are the same value of `1 / N`,
    /// where the size of the filter is `N`.
    pub fn create_moving_average(size: usize) -> Self {
        let coefficients = vec![1.0 / size as f32; size];
        Self {
            coeff: coefficients,
        }
    }

    /// Creates a Matched Filter / Square-root raised cosine filter.
    /// A Root Raised Cosine / Matched Filter has very low ISI / Inter Symbol Interference.
    /// `n_half` is half the length of the pulse (0,1,2...(n_half*2+1)).
    /// -> we want one sample in the middle of the pulse: `sinc(0)`, that is why the length of this
    /// filter is always uneven.
    /// `b` is beta / sinc roll-off and `m` is the oversampling factor.
    /// (Ported from Matlab: srrc.m, see Software Receiver Design book)
    pub fn create_srrc(n_half: usize, b: f32, m: usize) -> Self {
        // create a vector of size: N*2+1
        let mut srrc = Vec::with_capacity(n_half * 2 * m + 1);
        // offset the sinc function / rightshift
        let offset = n_half * m;
        // pre calculate a factor outside of loop
        let factor = 4.0 * b / f32::sqrt(m as f32);

        // create first half of signal, excluding middle sample
        for i in 0..(n_half * m) {
            let t = i as f32 - offset as f32;

            let cos_calc = f32::cos((1.0 + b) * PI * t / (m as f32));
            let sin_calc = f32::sin((1.0 - b) * PI * t / (m as f32)) / (4.0 * b * t / (m as f32));
            let latter_f = PI * (1.0 - 16.0 * (b * t / (m as f32)).powi(2));

            let mut val: f32 = factor * (cos_calc + sin_calc) / latter_f;
            // primitive check if NaN:
            // happens when t = (+/-) m/(4b)
            // we shall not divide by 0
            if val.is_nan() {
                val = b / f32::sqrt(2.0 * m as f32)
                    * ((1.0 + 2.0 / PI) * f32::sin(PI / (4.0 * b))
                        + (1.0 - 2.0 / PI) * f32::cos(PI / (4.0 * b)));
            }
            srrc.push(val as f32);
        }

        // calculate sinc for t=0, i.e.: t=offset, since right shifted
        srrc.push(1.0 / f32::sqrt(m as f32) * ((1.0 - b) + 4.0 * b / PI));

        // create second half of signal by reflection symmetry
        // mirrors signal at n=offset
        for i in 0..offset {
            srrc.push(srrc[offset - i - 1]);
        }

        Filter { coeff: srrc }
    }
}

#[cfg(test)]
mod test {
    use super::*;

    #[test]
    fn test_srrc_len() {
        // Size schould be 1 + 4*2*2
        let srrc = Filter::create_srrc(4, 0.5, 2);
        assert_eq!(srrc.len(), 17);
    }

    #[test]
    fn test_srrc_symmetry() {
        let srrc = Filter::create_srrc(4, 0.5, 2);

        // split srrc filter in middle and compare equality of sample i and N-i
        // e.g. 0 and len-0-1
        let len = srrc.len();
        for i in 0..len / 2 - 1 {
            assert_eq!(srrc.coeff[i], srrc.coeff[len - i - 1]);
        }
    }

    #[test]
    fn test_srrc_values() {
        let srrc = Filter::create_srrc(4, 0.5, 1);
        // these values are from octave: see Software Radio Design book:
        // -> included matlab file srrc.m
        // only really need to validate one side of the pulse (last sample is middle of srrc
        let srrc_mat = vec![-0.010105, 0.0030315, 0.042441, -0.1061, 1.1366];

        // diff/error between octave and my implementation should not exceed this
        let err = 0.0001;
        for i in 0..srrc_mat.len() {
            let err_sqrd = (srrc.coeff[i] - srrc_mat[i]).abs();
            assert!(err_sqrd < err)
        }
    }

    #[test]
    fn test_filter() {
        // N - M + 1 -> 7 - 3 + 1 = 5
        let filter = Filter::new(vec![1., 0., 1.]);
        let sig = vec![1., 2., 3., 4., 5., 6., 7.];

        let expected_result = vec![4., 6., 8., 10., 12.];
        assert_eq!(expected_result, filter.filter_windowed(&sig));
    }

    #[test]
    fn create_moving_average() {
        let vec = vec![0.2, 0.2, 0.2, 0.2, 0.2];
        let filter_ma_5 = Filter::new(vec);
        assert_eq!(filter_ma_5, Filter::create_moving_average(5));

        let vec = vec![0.25, 0.25, 0.25, 0.25];
        let filter_ma_4 = Filter::new(vec);
        assert_eq!(filter_ma_4, Filter::create_moving_average(4));
    }
}