phip1611 / spectrum-analyzer

An easy to use and fast `no_std` library (with `alloc`) to get the frequency spectrum of a digital signal (e.g. audio) using FFT.
MIT License
125 stars 20 forks source link

This crate says it's no_std but actually isn't (switch to no_std compatible FFT crate) #15

Closed Michael-F-Bryan closed 3 years ago

Michael-F-Bryan commented 3 years ago

This crate's main dependency, rustfft, unconditionally requires the standard library and therefore (despite the #![no_std] in lib.rs) the spectrum_analyzer crate actually isn't no_std.

The easiest way to check this is by finding a target that doesn't have a standard library (I used thumbv7em-none-eabihf - a common STM32 chi[) and running cargo check.

$ rustup target add thumbv7em-none-eabihf
info: downloading component 'rust-std' for 'thumbv7em-none-eabihf'
info: installing component 'rust-std' for 'thumbv7em-none-eabihf'

$ cargo check --target thumbv7em-none-eabihf
    Checking num-traits v0.2.14
    Checking log v0.4.14
error[E0463]: can't find crate for `std`
  --> /home/michael/.cargo/registry/src/github.com-1ecc6299db9ec823/num-traits-0.2.14/src/lib.rs:21:1
   |
21 | extern crate std;
   | ^^^^^^^^^^^^^^^^^ can't find crate
   |
   = note: the `thumbv7em-none-eabihf` target may not be installed

error: aborting due to previous error

The error in that message points to num-traits, but that's not an issue because num-traits has a std feature you can disable.

The problem is that the following crates aren't no_std and don't provide an enabled-by-default std feature for controlling the std feature of its dependencies:

phip1611 commented 3 years ago

Oh dang! Thanks for pointing this out @Michael-F-Bryan ! But I guess there is an FFT implementation that is really no_std out there. I will look for one soon and then evaluate performance differences and eventually replace the fft library or make this feature-dependent.

I'm confident I can work on this in the next 7 days.

phip1611 commented 3 years ago

@Michael-F-Bryan quick solution will be to use the complex functions in https://docs.rs/microfft/. I compared the results and it looks promising. Also, microfft internally compares it's results in tests against rustfft. With that, this crate should be finally really no_std.

Long-time idea: Add Cargo features and enable switching the FFT implementation (including STD-variants)

Michael-F-Bryan commented 3 years ago

Thanks for looking into this @phip1611.

I ended up adapting an implementation from Rosetta Code to our constraints (no allocations at runtime) so we don't have an urgent need for it.

I like the idea of using cargo features for enabling/disabling different implementations depending on their trade-offs. I could see you having something like this in your lib.rs:

// lib.rs

// Expose nice wrappers around both implementations
#[cfg(feature = "rustfft")]
mod rust;
#[cfg(feature = "microfft")]
mod micro;

// then provide a sane default they can use which will prefer the
// most performant implementation depending on the enabled 
// features.
cfg_if::cfg_if! {
  if #[cfg(feature = "rustfft")] {
    pub use rust::samples_fft_to_spectrum;
  } else if #[cfg(feature = "microfft")] {
    pub use micro::samples_fft_to_spectrum;
  }
}
phip1611 commented 3 years ago

@Michael-F-Bryan Unfortunately it seems like this is not trivial because "microfft" has a constellation of dependencies which results in this bug:

https://gitlab.com/ra_kete/microfft-rs/-/merge_requests/11

This problem is also in my crate if I use microfft as dependency, i.e. it's transitive. And just adding resolver = "2" results in other problems.

TL;DR: the code works but it doesn't build because of a Cargo bug. 🤔

UPDATE Damn, it works finally :D

phip1611 commented 3 years ago

next release will compile and work on no_std, but it still needs an allocator

brainstorm commented 5 months ago

I ended up adapting an implementation from Rosetta Code to our constraints (no allocations at runtime) so we don't have an urgent need for it.

@Michael-F-Bryan Could you please share that implementation, please? I'm looking for a fully no_std and no alloc crate and so far this mention of yours is the closest I've found after browsing github/crates.io/lib.rs :)

Michael-F-Bryan commented 5 months ago

@brainstorm it's been a while, but I think I was referring to this code: https://github.com/hotg-ai/spectrogram-computer/blob/master/src/fft.rs

brainstorm commented 5 months ago

@brainstorm it's been a while, but I think I was referring to this code: https://github.com/hotg-ai/spectrogram-computer/blob/master/src/fft.rs

404 :-S

Michael-F-Bryan commented 5 months ago

Sorry, looks like that repo is private. It's not exactly zero-allocation, but you should be able to swap Vec out for something like arrayvec::ArrayVec pretty easily.

Here's the source code:

use alloc::vec::Vec;
use core::f64::consts::PI;
use num_complex::Complex;

pub fn fft(samples: &[i16], _sample_rate: usize) -> Vec<i16> {
    let samples: Vec<_> = samples
        .iter()
        .map(|sample| normalize(*sample))
        .map(|sample| Complex::new(sample, 0.0))
        .collect();

    let mut buffer_a = Vec::new();
    let mut buffer_b = Vec::new();

    let frequency_domain = raw_fft(&samples, &mut buffer_a, &mut buffer_b);

    frequency_domain
        .iter()
        .copied()
        .map(|freq| denormalize(freq.re))
        .collect()
}

fn normalize(n: i16) -> f64 {
    (i16::max_value() as f64 - n as f64) / i16::max_value() as f64
}

fn denormalize(n: f64) -> i16 {
    (n * i16::max_value() as f64).round() as i16
}

/// A Fast Fourier Transform implementation copied from [Rosetta Code][rc],
/// modified to avoid allocating memory.
///
/// [rc]: https://rosettacode.org/wiki/Fast_Fourier_transform#Rust
pub fn raw_fft<'a>(
    input: &[Complex<f64>],
    buf_a: &'a mut Vec<Complex<f64>>,
    buf_b: &'a mut Vec<Complex<f64>>,
) -> &'a mut [Complex<f64>] {
    // round n (length) up to a power of 2:
    let n = input.len().next_power_of_two();

    // clear the previous buffer
    buf_a.clear();
    buf_a.reserve(n);
    // copy the input into a buffer
    buf_a.extend_from_slice(input);
    // pad out the rest with zeroes
    buf_a.extend(core::iter::repeat(Complex::default()).take(n - input.len()));

    // then copy it across to our second buffer
    buf_b.clear();
    buf_b.extend_from_slice(&buf_a);

    fft_recursive(buf_a, buf_b, n, 1);

    for element in buf_a.iter_mut() {
        *element /= n as f64;
    }

    buf_a
}

fn fft_recursive(
    buf_a: &mut [Complex<f64>],
    buf_b: &mut [Complex<f64>],
    n: usize,
    step: usize,
) {
    if step >= n {
        return;
    }

    fft_recursive(buf_b, buf_a, n, step * 2);
    fft_recursive(&mut buf_b[step..], &mut buf_a[step..], n, step * 2);
    let (left, right) = buf_a.split_at_mut(n / 2);

    for i in (0..n).step_by(step * 2) {
        let t = (Complex::new(0.0, -PI) * (i as f64) / (n as f64)).exp()
            * buf_b[i + step];
        left[i / 2] = buf_b[i] + t;
        right[i / 2] = buf_b[i] - t;
    }
}

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

    // Testing data from http://www.sccon.ca/sccon/fft/fft3.htm

    #[test]
    fn impulse() {
        let mut input = [Complex::new(0.0, 0.0); 8];
        input[0] = Complex::new(1.000, 0.000);
        let should_be = [Complex::new(0.125, 0.0); 8];
        let mut a = Vec::new();
        let mut b = Vec::new();

        let got = raw_fft(&input, &mut a, &mut b);

        assert_eq!(got, should_be);
    }

    #[test]
    fn shifted_impulse() {
        let mut input = [Complex::new(0.0, 0.0); 8];
        input[1] = Complex::new(1.000, 0.000);
        let should_be = [
            Complex::new(0.125, 0.000),
            Complex::new(0.088, -0.088),
            Complex::new(0.000, -0.125),
            Complex::new(-0.088, -0.088),
            Complex::new(-0.125, 0.000),
            Complex::new(-0.088, 0.088),
            Complex::new(0.000, 0.125),
            Complex::new(0.088, 0.088),
        ];
        let mut a = Vec::new();
        let mut b = Vec::new();

        let got = raw_fft(&input, &mut a, &mut b);

        println!("{:?}", got);
        assert!(
            got.iter().all(|c| c.norm() == 0.125),
            "Each result should have the same magnitude"
        );
        for (i, (got, should_be)) in got.iter().zip(&should_be).enumerate() {
            let error = got - should_be;
            assert!(
                error.norm() < 0.001,
                "{:?} != {:?} at index {}",
                got,
                should_be,
                i
            );
        }
    }

    #[test]
    fn all_zeroes() {
        let input = [Complex::new(0.0, 0.0); 8];
        let should_be = [Complex::new(0.0, 0.0); 8];
        let mut a = Vec::new();
        let mut b = Vec::new();

        let got = raw_fft(&input, &mut a, &mut b);

        assert_eq!(got, should_be);
    }
}