diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..f17d192 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,12 @@ +# Changelog + +## v0.3.0 +- `FrequencySpectrum::min()` and `FrequencySpectrum::max()` + now return tuples/pairs (issue #6) +- `FrequencySpectrum` now has convenient methods to get + the value of a desired frequency from the underlying vector + of FFT results (issue #8) +- `FrequencySpectrum` now has a field + getter for `frequency_resolution` + (issue #5) + +For all issues see: https://github.com/phip1611/spectrum-analyzer/milestone/2 diff --git a/Cargo.toml b/Cargo.toml index 3b09bc2..5a0d3b2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,7 +4,7 @@ description = """ A simple and fast `no_std` library to get the frequency spectrum of a digital signal (e.g. audio) using FFT. It follows the KISS principle and consists of simple building blocks/optional features. """ -version = "0.2.2" +version = "0.3.0" authors = ["Philipp Schuster "] edition = "2018" keywords = ["fft", "spectrum", "frequencies", "audio", "dsp"] @@ -18,6 +18,7 @@ documentation = "https://docs.rs/spectrum-analyzer" [dependencies] rustfft = "5.0.1" +float-cmp = "0.8.0" # approx. compare floats [dev-dependencies] minimp3 = "0.5.1" diff --git a/README.md b/README.md index 1322712..b016dde 100644 --- a/README.md +++ b/README.md @@ -1,17 +1,21 @@ # Rust: library for frequency spectrum analysis using FFT A simple and fast `no_std` library to get the frequency spectrum of a digital signal (e.g. audio) using FFT. -It follows the KISS principle and consists of simple building blocks/optional features. +It follows the KISS principle and consists of simple building blocks/optional features. In short, this is +a convenient wrapper around the great [rustfft](https://crates.io/crates/rustfft) +library. **I'm not an expert on digital signal processing. Code contributions are highly welcome! :)** ## How to use +Most tips and comments are located inside the code, so please check out the repository on +Github! Anyway, the most basic usage looks like this: ```rust use spectrum_analyzer::{samples_fft_to_spectrum, FrequencyLimit}; use spectrum_analyzer::windows::hann_window; fn main() { // This lib also works in `no_std` environments! - let samples: &[f32] = get_samples(); // TODO implement + let samples: &[f32] = get_samples(); // TODO you need to implement the samples source // apply hann window for smoothing; length must be a power of 2 for the FFT let hann_window = hann_window(&samples[0..4096]); // calc spectrum @@ -34,30 +38,10 @@ fn main() { } ``` -### How to scale values -```rust -// e.g. like this -fn get_scale_to_one_fn_factory() -> SpectrumTotalScaleFunctionFactory{ - Box::new( - move |min: f32, max: f32, average: f32, median: f32| { - Box::new( - move |x| x/max - ) - } - ) -} -fn main() { - // ... - let spectrum_hann_window = samples_fft_to_spectrum( - &hann_window, - 44100, - FrequencyLimit::All, - None, - // optional total scaling at the end; see doc comments - Some(get_scale_to_one_fn_factory()), - ); -} -``` +## Scaling the frequency values/amplitudes +As already mentioned, there are lots of comments in the code. Short story is: +Type `ComplexSpectrumScalingFunction` can do anything whereas `BasicSpectrumScalingFunction` +is easier to write, especially for Rust beginners. ## Performance *Measurements taken on i7-8650U @ 3 Ghz (Single-Core) with optimized build* diff --git a/examples/mp3-samples.rs b/examples/mp3-samples.rs index dc25cb1..f10cedf 100644 --- a/examples/mp3-samples.rs +++ b/examples/mp3-samples.rs @@ -27,9 +27,7 @@ use audio_visualizer::spectrum::staticc::plotters_png_file::spectrum_static_plot use audio_visualizer::test_support::TEST_OUT_DIR; use minimp3::{Decoder as Mp3Decoder, Error as Mp3Error, Frame as Mp3Frame}; use spectrum_analyzer::windows::{blackman_harris_4term, hamming_window, hann_window}; -use spectrum_analyzer::{ - samples_fft_to_spectrum, FrequencyLimit, SpectrumTotalScaleFunctionFactory, -}; +use spectrum_analyzer::{samples_fft_to_spectrum, FrequencyLimit, scaling}; use std::fs::File; use std::time::Instant; @@ -124,7 +122,7 @@ fn to_spectrum_and_plot(samples: &[f32], filename: &str, frequency_limit: Freque 44100, frequency_limit, None, - Some(get_scale_to_one_fn_factory()), + Some(scaling::complex::scale_to_zero_to_one()), ); println!( "[Measurement]: FFT to Spectrum with no window with {} samples took: {}µs", @@ -137,7 +135,7 @@ fn to_spectrum_and_plot(samples: &[f32], filename: &str, frequency_limit: Freque 44100, frequency_limit, None, - Some(get_scale_to_one_fn_factory()), + Some(scaling::complex::scale_to_zero_to_one()), ); println!( "[Measurement]: FFT to Spectrum with Hamming window with {} samples took: {}µs", @@ -150,7 +148,7 @@ fn to_spectrum_and_plot(samples: &[f32], filename: &str, frequency_limit: Freque 44100, frequency_limit, None, - Some(get_scale_to_one_fn_factory()), + Some(scaling::complex::scale_to_zero_to_one()), ); println!( "[Measurement]: FFT to Spectrum with Hann window with {} samples took: {}µs", @@ -163,7 +161,7 @@ fn to_spectrum_and_plot(samples: &[f32], filename: &str, frequency_limit: Freque 44100, frequency_limit, None, - Some(get_scale_to_one_fn_factory()), + Some(scaling::complex::scale_to_zero_to_one()), ); println!("[Measurement]: FFT to Spectrum with Blackmann Harris 4-term window with {} samples took: {}µs", samples.len(), now.elapsed().as_micros()); let now = Instant::now(); @@ -172,7 +170,7 @@ fn to_spectrum_and_plot(samples: &[f32], filename: &str, frequency_limit: Freque 44100, frequency_limit, None, - Some(get_scale_to_one_fn_factory()), + Some(scaling::complex::scale_to_zero_to_one()), ); println!("[Measurement]: FFT to Spectrum with Blackmann Harris 7-term window with {} samples took: {}µs", samples.len(), now.elapsed().as_micros()); @@ -216,10 +214,6 @@ fn to_spectrum_and_plot(samples: &[f32], filename: &str, frequency_limit: Freque ); } -fn get_scale_to_one_fn_factory() -> SpectrumTotalScaleFunctionFactory { - Box::new(move |_min: f32, max: f32, _average: f32, _median: f32| Box::new(move |x| x / max)) -} - fn read_mp3(file: &str) -> Vec { let samples = read_mp3_to_mono(file); let samples = samples.into_iter().map(|x| x as f32).collect::>(); diff --git a/src/lib.rs b/src/lib.rs index 789ec26..747361f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -24,6 +24,8 @@ SOFTWARE. //! A simple and fast `no_std` library to get the frequency spectrum of a digital signal //! (e.g. audio) using FFT. It follows the KISS principle and consists of simple building //! blocks/optional features. +//! +//! In short, this is a convenient wrapper around the great `rustfft` library. #![no_std] @@ -43,15 +45,25 @@ use rustfft::{Fft, FftDirection}; pub use crate::frequency::{Frequency, FrequencyValue}; pub use crate::limit::FrequencyLimit; -pub use crate::spectrum::{FrequencySpectrum, SpectrumTotalScaleFunctionFactory}; +pub use crate::spectrum::{FrequencySpectrum, ComplexSpectrumScalingFunction}; use core::convert::identity; mod frequency; mod limit; mod spectrum; +pub mod scaling; +pub mod windows; #[cfg(test)] mod tests; -pub mod windows; + +/// Definition of a simple function that gets applied on each frequency magnitude +/// in the spectrum. This is easier to write, especially for Rust beginners. +/// Everything that can be achieved with this, can also be achieved with parameter +/// `total_scaling_fn`. +/// +/// The scaling only affects the value/amplitude of the frequency +/// but not the frequency itself. +pub type SimpleSpectrumScalingFunction<'a> = &'a dyn Fn(f32) -> f32; /// Takes an array of samples (length must be a power of 2), /// e.g. 2048, applies an FFT (using library `rustfft`) on it @@ -63,11 +75,13 @@ pub mod windows; /// e.g. `44100/(16384/2) == 5.383Hz`, i.e. more samples => better accuracy /// * `sampling_rate` sampling_rate, e.g. `44100 [Hz]` /// * `frequency_limit` Frequency limit. See [`FrequencyLimit´] -/// * `per_element_scaling_fn` Optional per element scaling function, e.g. `20 * log(x)`. -/// To see where this equation comes from, check out -/// this paper: -/// https://www.sjsu.edu/people/burford.furman/docs/me120/FFT_tutorial_NI.pdf -/// * `total_scaling_fn` See [`crate::spectrum::SpectrumTotalScaleFunctionFactory`]. +/// * `per_element_scaling_fn` See [`crate::SimpleSpectrumScalingFunction`] for details. +/// This is easier to write, especially for Rust beginners. Everything +/// that can be achieved with this, can also be achieved with +/// parameter `total_scaling_fn`. +/// See [`crate::scaling`] for example implementations. +/// * `total_scaling_fn` See [`crate::spectrum::SpectrumTotalScaleFunctionFactory`] for details. +/// See [`crate::scaling`] for example implementations. /// /// ## Returns value /// New object of type [`FrequencySpectrum`]. @@ -75,8 +89,8 @@ pub fn samples_fft_to_spectrum( samples: &[f32], sampling_rate: u32, frequency_limit: FrequencyLimit, - per_element_scaling_fn: Option<&dyn Fn(f32) -> f32>, - total_scaling_fn: Option, + per_element_scaling_fn: Option, + total_scaling_fn: Option, ) -> FrequencySpectrum { // With FFT we transform an array of time-domain waveform samples // into an array of frequency-domain spectrum samples @@ -150,13 +164,19 @@ fn fft_result_to_frequency_to_magnitude_map( sampling_rate: u32, frequency_limit: FrequencyLimit, per_element_scaling_fn: Option<&dyn Fn(f32) -> f32>, - total_scaling_fn: Option, + total_scaling_fn: Option, ) -> FrequencySpectrum { let maybe_min = frequency_limit.maybe_min(); let maybe_max = frequency_limit.maybe_max(); let samples_len = fft_result.len(); + // see documentation of fft_calc_frequency_resolution for better explanation + let frequency_resolution = fft_calc_frequency_resolution( + sampling_rate, + samples_len as u32, + ); + // collect frequency => frequency value in Vector of Pairs/Tuples let frequency_vec = fft_result .into_iter() @@ -167,7 +187,9 @@ fn fft_result_to_frequency_to_magnitude_map( // calc index => corresponding frequency .map(|(fft_index, complex)| { ( - fft_index_to_corresponding_frequency(fft_index, samples_len as u32, sampling_rate), + // corresponding frequency of each index of FFT result + // see documentation of fft_calc_frequency_resolution for better explanation + fft_index as f32 * frequency_resolution, complex, ) }) @@ -202,41 +224,52 @@ fn fft_result_to_frequency_to_magnitude_map( .collect::>(); // create spectrum object - let fs = FrequencySpectrum::new(frequency_vec); + let spectrum = FrequencySpectrum::new( + frequency_vec, + frequency_resolution, + ); + // optionally scale if let Some(total_scaling_fn) = total_scaling_fn { - fs.apply_total_scaling_fn(total_scaling_fn) + spectrum.apply_complex_scaling_fn(total_scaling_fn) } - fs + + spectrum } -/// Calculate what index in the FFT result corresponds to what frequency. +/// Calculate the frequency resolution of the FFT. It is determined by the sampling rate +/// in Hertz and N, the number of samples given into the FFT. With the frequency resolution, +/// we can determine the corresponding frequency of each index in the FFT result buffer. /// /// ## Parameters -/// * `fft_index` Index in FFT result buffer. If `samples.len() == 2048` then this is in `{0, 1, ..., 1023}` /// * `samples_len` Number of samples put into the FFT /// * `sampling_rate` sampling_rate, e.g. `44100 [Hz]` /// /// ## Return value +/// Frequency resolution in Hertz. +/// +/// ## More info +/// * https://www.researchgate.net/post/How-can-I-define-the-frequency-resolution-in-FFT-And-what-is-the-difference-on-interpreting-the-results-between-high-and-low-frequency-resolution +/// * https://stackoverflow.com/questions/4364823/ #[inline(always)] -fn fft_index_to_corresponding_frequency( - fft_index: usize, - samples_len: u32, +fn fft_calc_frequency_resolution( sampling_rate: u32, + samples_len: u32, ) -> f32 { // Explanation for the algorithm: // https://stackoverflow.com/questions/4364823/ // samples : [0], [1], [2], [3], ... , ..., [2047] => 2048 samples for example // FFT Result : [0], [1], [2], [3], ... , ..., [2047] - // Relevant part of FFT Result: [0], [1], [2], [3], ... , [1023] + // Relevant part of FFT Result: [0], [1], [2], [3], ... , [1023] => first N/2 results important // ^ ^ // Frequency : 0Hz, .................... Sampling Rate/2 - // 0Hz is also called (e.g. 22050Hz @ 44100H sampling rate) + // 0Hz is also called (e.g. 22050Hz for 44100H sampling rate) // "DC Component" - // frequency step/resolution is for example: 1/1024 * 44100 - // 1024: relevant FFT result, 2048 samples, 44100 sample rate + // frequency step/resolution is for example: 1/2048 * 44100 + // 2048 samples, 44100 sample rate - fft_index as f32 / samples_len as f32 * sampling_rate as f32 + // equal to: 1.0 / samples_len as f32 * sampling_rate as f32 + sampling_rate as f32 / samples_len as f32 } diff --git a/src/scaling.rs b/src/scaling.rs new file mode 100644 index 0000000..76cbafe --- /dev/null +++ b/src/scaling.rs @@ -0,0 +1,67 @@ +/* +MIT License + +Copyright (c) 2021 Philipp Schuster + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ +//! This module contains convenient public transform functions that you can use +//! as parameters in [`crate::samples_fft_to_spectrum`]. + +/// Practical implementations for [`crate::SimpleSpectrumScalingFunction`]. +pub mod basic { + /// Calculates the base 10 logarithm of each frequency magnitude and + /// multiplies it with 20. This scaling is quite common, you can + /// find more information for example here: + /// https://www.sjsu.edu/people/burford.furman/docs/me120/FFT_tutorial_NI.pdf + /// + /// ## Usage + /// ```rust + ///use spectrum_analyzer::{samples_fft_to_spectrum, scaling, FrequencyLimit}; + ///let window = [0.0, 0.1, 0.2, 0.3]; // add real data here + ///let spectrum = samples_fft_to_spectrum( + /// &window, + /// 44100, + /// FrequencyLimit::All, + /// Some(&scaling::basic::scale_20_times_log10), + /// None, + /// ); + /// ``` + pub fn scale_20_times_log10(frequency_magnitude: f32) -> f32 { + 20.0 * frequency_magnitude.log10() + } +} + +/// Practical implementations for [`crate::spectrum::ComplexSpectrumScalingFunction`]. +pub mod complex { + use alloc::boxed::Box; + use crate::ComplexSpectrumScalingFunction; + + /// Returns a function factory that generates a function that scales + /// each frequency value/amplitude in the spectrum to interval `[0.0; 1.0]`. + pub fn scale_to_zero_to_one() -> ComplexSpectrumScalingFunction { + Box::new( + move |_min: f32, max: f32, _average: f32, _median: f32| { + Box::new( + move |x| x / max + ) + } + ) + } +} diff --git a/src/spectrum.rs b/src/spectrum.rs index 37e3076..35273c3 100644 --- a/src/spectrum.rs +++ b/src/spectrum.rs @@ -31,12 +31,23 @@ use core::cell::{Cell, Ref, RefCell}; /// Describes the type for a function factory that generates a function that can scale/normalize /// the data inside [`FrequencySpectrum`]. -/// This can be used to subtract `min` value from all values for example , if `min` -/// is `> 0`. The signature is the following: +/// +/// **Complex** means it's not basic. It has nothing to do with complex numbers. +/// +/// This can archive exactly the same as [`crate::SimpleSpectrumScalingFunction`] +/// but is capable of doing more complex logic, i.e. you have access to min, max, +/// or median! +/// +/// This can be used for example to subtract `min` value from all values, +/// if `min` is `> 0`. The signature is the following: /// `(min: f32, max: f32, average: f32, median: f32) -> fn(f32) -> f32` /// i.e. you provide a function which generates a function that gets -/// applied to each element. -pub type SpectrumTotalScaleFunctionFactory = +/// applied to each element. The input arguments are automatically calculated +/// by [`FrequencySpectrum`]. +/// +/// The scaling only affects the value/amplitude of the frequency +/// but not the frequency itself. +pub type ComplexSpectrumScalingFunction = Box Box f32>>; /// Convenient wrapper around the processed FFT result which describes each frequency and @@ -46,29 +57,57 @@ pub type SpectrumTotalScaleFunctionFactory = #[derive(Debug)] pub struct FrequencySpectrum { /// Raw data. Vector is sorted from lowest - /// frequency to highest. + /// frequency to highest and data is normalized/scaled + /// according to all applied scaling functions. data: RefCell>, - /// Average value of frequency value/magnitude/amplitude. + /// Frequency resolution of the examined samples in Hertz, + /// i.e the frequency steps between elements in the vector + /// inside field [`data`]. + frequency_resolution: f32, + /// Average value of frequency value/magnitude/amplitude + /// corresponding to data in [`FrequencySpectrum::data`]. average: Cell, - /// Median value of frequency value/magnitude/amplitude. + /// Median value of frequency value/magnitude/amplitude + /// corresponding to data in [`FrequencySpectrum::data`]. median: Cell, - /// Minimum value of frequency value/magnitude/amplitude. - min: Cell, - /// Maximum value of frequency value/magnitude/amplitude. - max: Cell, + /// Pair of (frequency, frequency value/magnitude/amplitude) where + /// frequency value is **minimal** inside the spectrum. + /// Corresponding to data in [`FrequencySpectrum::data`]. + min: Cell<(Frequency, FrequencyValue)>, + /// Pair of (frequency, frequency value/magnitude/amplitude) where + /// frequency value is **maximum** inside the spectrum. + /// Corresponding to data in [`FrequencySpectrum::data`]. + max: Cell<(Frequency, FrequencyValue)>, } impl FrequencySpectrum { /// Creates a new object. Calculates several metrics on top of /// the passed vector. + /// + /// ## Parameters + /// * `data` Vector with all ([`Frequency`], [`FrequencyValue`])-tuples + /// * `frequency_resolution` Resolution in Hertz. This equals to + /// `data[1].0 - data[0].0`. #[inline(always)] - pub fn new(data: Vec<(Frequency, FrequencyValue)>) -> Self { + pub fn new( + data: Vec<(Frequency, FrequencyValue)>, + frequency_resolution: f32, + ) -> Self { let obj = Self { data: RefCell::new(data), + frequency_resolution, + + // default/placeholder values average: Cell::new(FrequencyValue::from(-1.0)), median: Cell::new(FrequencyValue::from(-1.0)), - min: Cell::new(FrequencyValue::from(-1.0)), - max: Cell::new(FrequencyValue::from(-1.0)), + min: Cell::new(( + Frequency::from(-1.0), + FrequencyValue::from(-1.0), + )), + max: Cell::new(( + Frequency::from(-1.0), + FrequencyValue::from(-1.0), + )), }; // IMPORTANT!! obj.calc_statistics(); @@ -81,11 +120,11 @@ impl FrequencySpectrum { /// ## Parameters /// * `total_scaling_fn` See [`crate::spectrum::SpectrumTotalScaleFunctionFactory`]. #[inline(always)] - pub fn apply_total_scaling_fn(&self, total_scaling_fn: SpectrumTotalScaleFunctionFactory) { + pub fn apply_complex_scaling_fn(&self, total_scaling_fn: ComplexSpectrumScalingFunction) { let scale_fn = (total_scaling_fn)( // into() => FrequencyValue => f32 - self.min.get().val(), - self.max.get().val(), + self.min.get().1.val(), + self.max.get().1.val(), self.average.get().val(), self.median.get().val(), ); @@ -101,42 +140,203 @@ impl FrequencySpectrum { self.calc_statistics(); } - /// Getter for lazily evaluated `average`. + /// Getter for [`FrequencySpectrum::average`]. #[inline(always)] pub fn average(&self) -> FrequencyValue { self.average.get() } - /// Getter for lazily evaluated `median`. + /// Getter for [`FrequencySpectrum::median`]. #[inline(always)] pub fn median(&self) -> FrequencyValue { self.median.get() } - /// Getter for lazily evaluated `max`. + /// Getter for [`FrequencySpectrum::max`]. #[inline(always)] - pub fn max(&self) -> FrequencyValue { + pub fn max(&self) -> (Frequency, FrequencyValue) { self.max.get() } - /// Getter for lazily evaluated `min`. + /// Getter for [`FrequencySpectrum::min`]. #[inline(always)] - pub fn min(&self) -> FrequencyValue { + pub fn min(&self) -> (Frequency, FrequencyValue) { self.min.get() } - /// [`max()`] - [`min()`] + /// Returns [`FrequencySpectrum::max().1`] - [`FrequencySpectrum::min().1`], + /// i.e. the range of the frequency values (not the frequencies itself, + /// but their amplitude/value). #[inline(always)] pub fn range(&self) -> FrequencyValue { - self.max() - self.min() + self.max().1 - self.min().1 } - /// Getter for `data`. + /// Getter for [`FrequencySpectrum::data`]. #[inline(always)] pub fn data(&self) -> Ref> { self.data.borrow() } + /// Getter for [`FrequencySpectrum::frequency_resolution`]. + #[inline(always)] + pub fn frequency_resolution(&self) -> f32 { + self.frequency_resolution + } + + /// Returns the value of the given frequency from the spectrum either exactly or approximated. + /// If `search_fr` is not exactly given in the spectrum, i.e. due to the + /// [`self::frequency_resolution`], this function takes the two closest + /// neighbors/points (A, B), put a linear function through them and calculates + /// the point C in the middle. This is done by using + /// [`calculate_point_between_points`]. + /// + /// ## Panics + /// If parameter `search_fr` (frequency) is below the lowest or the maximum + /// frequency, this function panics! This is because the user provide + /// the min/max frequency when the spectrum is created and knows about it. + /// This is similar to an intended "out of bounds"-access. + /// + /// ## Parameters + /// - `search_fr` The frequency of that you want the amplitude/value in the spectrum. + /// + /// ## Return + /// Either exact value of approximated value, determined by [`self::frequency_resolution`]. + #[inline(always)] + pub fn freq_val_exact(&self, search_fr: f32) -> FrequencyValue { + let data = self.data.borrow(); + + // lowest frequency in the spectrum + // TODO use minFrequency() and maxFrequency() + let (min_fr, min_fr_val) = data[0]; + // highest frequency in the spectrum + let (max_fr, max_fr_val) = data[data.len() - 1]; + + // https://docs.rs/float-cmp/0.8.0/float_cmp/ + let equals_min_fr = float_cmp::approx_eq!(f32, min_fr.val(), search_fr, ulps = 3); + let equals_max_fr = float_cmp::approx_eq!(f32, max_fr.val(), search_fr, ulps = 3); + + // Fast return if possible + if equals_min_fr { + return min_fr_val; + } + if equals_max_fr { + return max_fr_val; + } + // bounds check + if search_fr < min_fr.val() || search_fr > max_fr.val() { + panic!("Frequency {}Hz is out of bounds [{}; {}]!", search_fr, min_fr.val(), max_fr.val()); + } + + // We search for Point C (x=search_fr, y=???) between Point A and Point B iteratively. + // Point B is always the successor of A. + + for two_points in data.iter().as_slice().windows(2) { + let point_a = two_points[0]; + let point_b = two_points[1]; + let point_a_x = point_a.0.val(); + let point_a_y = point_a.1; + let point_b_x = point_b.0.val(); + let point_b_y = point_b.1.val(); + + // check if we are in the correct window; we are in the correct window + // iff point_a_x <= search_fr <= point_b_x + if search_fr > point_b_x { + continue; + } + + return if float_cmp::approx_eq!(f32, point_a_x, search_fr, ulps = 3) { + // directly return if possible + point_a_y + } else { + calculate_y_coord_between_points( + (point_a_x, point_a_y.val()), + (point_b_x, point_b_y), + search_fr, + ).into() + } + } + + panic!("Here be dragons"); + } + + /// Returns the frequency closest to parameter `search_fr` in the spectrum. For example + /// if the spectrum looks like this: + /// ```text + /// Vector: [0] [1] [2] [3] + /// Frequency 100 Hz 200 Hz 300 Hz 400 Hz + /// Fr Value 0.0 1.0 0.5 0.1 + /// ``` + /// then `get_frequency_value_closest(320)` will return `(300.0, 0.5)`. + /// + /// ## Panics + /// If parameter `search_fre` (frequency) is below the lowest or the maximum + /// frequency, this function panics! + /// + /// ## Parameters + /// - `search_fr` The frequency of that you want the amplitude/value in the spectrum. + /// + /// ## Return + /// Closest matching point in spectrum, determined by [`self::frequency_resolution`]. + #[inline(always)] + pub fn freq_val_closest(&self, search_fr: f32) -> (Frequency, FrequencyValue) { + let data = self.data.borrow(); + + // lowest frequency in the spectrum + // TODO use minFrequency() and maxFrequency() + let (min_fr, min_fr_val) = data[0]; + // highest frequency in the spectrum + let (max_fr, max_fr_val) = data[data.len() - 1]; + + // https://docs.rs/float-cmp/0.8.0/float_cmp/ + let equals_min_fr = float_cmp::approx_eq!(f32, min_fr.val(), search_fr, ulps = 3); + let equals_max_fr = float_cmp::approx_eq!(f32, max_fr.val(), search_fr, ulps = 3); + + // Fast return if possible + if equals_min_fr { + return (min_fr, min_fr_val); + } + if equals_max_fr { + return (max_fr, max_fr_val); + } + + // bounds check + if search_fr < min_fr.val() || search_fr > max_fr.val() { + panic!("Frequency {}Hz is out of bounds [{}; {}]!", search_fr, min_fr.val(), max_fr.val()); + } + + for two_points in data.iter().as_slice().windows(2) { + let point_a = two_points[0]; + let point_b = two_points[1]; + let point_a_x = point_a.0; + let point_a_y = point_a.1; + let point_b_x = point_b.0; + let point_b_y = point_b.1; + + // check if we are in the correct window; we are in the correct window + // iff point_a_x <= search_fr <= point_b_x + if search_fr > point_b_x.val() { + continue; + } + + return if float_cmp::approx_eq!(f32, point_a_x.val(), search_fr, ulps = 3) { + // directly return if possible + (point_a_x, point_a_y) + } else { + // absolute difference + let delta_to_a = search_fr - point_a_x.val(); + // let delta_to_b = point_b_x.val() - search_fr; + if delta_to_a / self.frequency_resolution < 0.5 { + (point_a_x, point_a_y) + } else { + (point_b_x, point_b_y) + } + } + } + + panic!("Here be dragons"); + } + /// Returns a `BTreeMap`. The key is of type u32. /// (`f32` is not `Ord`, hence we can't use it as key.) You can optionally specify a /// scale function, e.g. multiply all frequencies with 1000 for better @@ -176,35 +376,35 @@ impl FrequencySpectrum { /// Calculates min, max, median and average of the frequency values/magnitudes/amplitudes. #[inline(always)] fn calc_statistics(&self) { - let data = self.data.borrow(); - // first: order all by frequency value in ascending order - let mut vals = data - .iter() - // map to only value - .map(|(_fr, val)| val) - // f64 to prevent overflow - // .map(|v| v as f64) - .collect::>(); - vals.sort(); + let mut data_sorted = self.data.borrow().clone(); + data_sorted.sort_by(|(_l_fr, l_fr_val), (_r_fr, r_fr_val)| { + // compare by frequency value, from min to max + l_fr_val.cmp(r_fr_val) + }); // sum - let sum: f32 = vals + let sum: f32 = data_sorted .iter() - .map(|fr_val| fr_val.val()) + .map(|fr_val| fr_val.1.val()) .fold(0.0, |a, b| a + b); - let avg = sum / vals.len() as f32; + let avg = sum / data_sorted.len() as f32; let average: FrequencyValue = avg.into(); let median = { - // we assume that vals.length() is always even, because + // we assume that data_sorted.length() is always even, because // it must be a power of 2 (for FFT) - let a = *vals[vals.len() / 2 - 1]; - let b = *vals[vals.len() / 2]; + let a = data_sorted[data_sorted.len() / 2 - 1].1; + let b = data_sorted[data_sorted.len() / 2].1; (a + b) / 2.0.into() }; - let min = *vals[0]; - let max = *vals[vals.len() - 1]; + // because we sorted the vector a few lines above + // by the value, the following lines are correct + let min = data_sorted[0]; + let max = data_sorted[data_sorted.len() - 1]; + + // check that I get the comparison right (and not from max to min) + debug_assert!(min.1 <= max.1); self.min.replace(min); self.max.replace(max); @@ -228,12 +428,68 @@ impl FrequencySpectrum { } }*/ +/// Calculates the y coordinate of Point C between two given points A and B +/// if the x-coordinate of C is known. It does that by putting a linear function +/// through the two given points. +/// +/// ## Parameters +/// - `(x1, y1)` x and y of point A +/// - `(x2, y2)` x and y of point B +/// - `x_coord` x coordinate of searched point C +/// +/// ## Return Value +/// y coordinate of searched point C +#[inline(always)] +fn calculate_y_coord_between_points((x1, y1): (f32, f32), (x2, y2): (f32, f32), x_coord: f32) -> f32 { + // e.g. Points (100, 1.0) and (200, 0.0) + // y=f(x)=-0.01x + c + // 1.0 = f(100) = -0.01x + c + // c = 1.0 + 0.01*100 = 2.0 + // y=f(180)=-0.01*180 + 2.0 + + + // gradient, anstieg + let slope = (y2 - y1)/(x2 - x1); + // calculate c in y=f(x)=slope * x + c + let c = y1 - slope * x1; + + slope * x_coord + c +} + #[cfg(test)] mod tests { use super::*; #[test] - fn test_spectrum() { + fn test_calculate_point_between_points() { + assert_eq!( + // expected y coordinate + 0.5, + calculate_y_coord_between_points( + (100.0, 1.0), + (200.0, 0.0), + 150.0, + ), + "Must calculate middle point between points by laying a linear function through the two points" + ); + assert!( + // https://docs.rs/float-cmp/0.8.0/float_cmp/ + float_cmp::approx_eq!( + f32, + 0.2, + calculate_y_coord_between_points( + (100.0, 1.0), + (200.0, 0.0), + 180.0, + ), + ulps = 3 + ), + "Must calculate arbitrary point between points by laying a linear function through the two points" + ); + } + + #[test] + fn test_spectrum_basic() { let spectrum = vec![ (0.0_f32, 5.0_f32), (50.0, 50.0), @@ -249,57 +505,219 @@ mod tests { .into_iter() .map(|(fr, val)| (fr.into(), val.into())) .collect::>(); - let spectrum = FrequencySpectrum::new(spectrum); - - assert_eq!( - (0.0.into(), 5.0.into()), - spectrum.data()[0], - "Vector must be ordered" + let spectrum = FrequencySpectrum::new( + spectrum, + 50.0, ); - assert_eq!( - (50.0.into(), 50.0.into()), - spectrum.data()[1], - "Vector must be ordered" - ); - assert_eq!( - (100.0.into(), 100.0.into()), - spectrum.data()[2], - "Vector must be ordered" - ); - assert_eq!( - (150.0.into(), 150.0.into()), - spectrum.data()[3], - "Vector must be ordered" - ); - assert_eq!( - (200.0.into(), 100.0.into()), - spectrum.data()[4], - "Vector must be ordered" - ); - assert_eq!( - (250.0.into(), 20.0.into()), - spectrum.data()[5], - "Vector must be ordered" + + // test inner vector is ordered + { + assert_eq!( + (0.0.into(), 5.0.into()), + spectrum.data()[0], + "Vector must be ordered" + ); + assert_eq!( + (50.0.into(), 50.0.into()), + spectrum.data()[1], + "Vector must be ordered" + ); + assert_eq!( + (100.0.into(), 100.0.into()), + spectrum.data()[2], + "Vector must be ordered" + ); + assert_eq!( + (150.0.into(), 150.0.into()), + spectrum.data()[3], + "Vector must be ordered" + ); + assert_eq!( + (200.0.into(), 100.0.into()), + spectrum.data()[4], + "Vector must be ordered" + ); + assert_eq!( + (250.0.into(), 20.0.into()), + spectrum.data()[5], + "Vector must be ordered" + ); + assert_eq!( + (300.0.into(), 0.0.into()), + spectrum.data()[6], + "Vector must be ordered" + ); + assert_eq!( + (450.0.into(), 200.0.into()), + spectrum.data()[7], + "Vector must be ordered" + ); + } + + // test getters + { + assert_eq!((300.0.into(), 0.0.into()), spectrum.min(), "min() must work"); + assert_eq!((450.0.into(), 200.0.into()), spectrum.max(), "max() must work"); + assert_eq!(200.0 - 0.0, spectrum.range().val(), "range() must work"); + assert_eq!(78.125, spectrum.average().val(), "average() must work"); + assert_eq!( + (50 + 100) as f32 / 2.0, + spectrum.median().val(), + "median() must work" + ); + assert_eq!( + 50.0, + spectrum.frequency_resolution(), + "frequency resolution must be returned" + ); + } + + // test get frequency exact + { + assert_eq!( + 5.0, + spectrum.freq_val_exact(0.0).val(), + ); + assert_eq!( + 50.0, + spectrum.freq_val_exact(50.0).val(), + ); + assert_eq!( + 150.0, + spectrum.freq_val_exact(150.0).val(), + ); + assert_eq!( + 100.0, + spectrum.freq_val_exact(200.0).val(), + ); + assert_eq!( + 20.0, + spectrum.freq_val_exact(250.0).val(), + ); + assert_eq!( + 0.0, + spectrum.freq_val_exact(300.0).val(), + ); + assert_eq!( + 100.0, + spectrum.freq_val_exact(375.0).val(), + ); + assert_eq!( + 200.0, + spectrum.freq_val_exact(450.0).val(), + ); + } + + // test get frequency closest + { + assert_eq!( + (0.0.into(), 5.0.into()), + spectrum.freq_val_closest(0.0), + ); + assert_eq!( + (50.0.into(), 50.0.into()), + spectrum.freq_val_closest(50.0), + ); + assert_eq!( + (450.0.into(), 200.0.into()), + spectrum.freq_val_closest(450.0), + ); + assert_eq!( + (450.0.into(), 200.0.into()), + spectrum.freq_val_closest(448.0), + ); + assert_eq!( + (450.0.into(), 200.0.into()), + spectrum.freq_val_closest(400.0), + ); + assert_eq!( + (50.0.into(), 50.0.into()), + spectrum.freq_val_closest(47.3), + ); + assert_eq!( + (50.0.into(), 50.0.into()), + spectrum.freq_val_closest(51.3), + ); + } + } + + #[test] + #[should_panic] + fn test_spectrum_get_frequency_value_exact_panic_below_min() { + let spectrum_vector = vec![ + (0.0_f32, 5.0_f32), + (450.0, 200.0), + ]; + + let spectrum = spectrum_vector + .into_iter() + .map(|(fr, val)| (fr.into(), val.into())) + .collect::>(); + let spectrum = FrequencySpectrum::new( + spectrum, + 50.0, ); - assert_eq!( - (300.0.into(), 0.0.into()), - spectrum.data()[6], - "Vector must be ordered" + + spectrum.freq_val_exact(-1.0).val(); + } + + #[test] + #[should_panic] + fn test_spectrum_get_frequency_value_exact_panic_below_max() { + let spectrum_vector = vec![ + (0.0_f32, 5.0_f32), + (450.0, 200.0), + ]; + + let spectrum = spectrum_vector + .into_iter() + .map(|(fr, val)| (fr.into(), val.into())) + .collect::>(); + let spectrum = FrequencySpectrum::new( + spectrum, + 50.0, ); - assert_eq!( - (450.0.into(), 200.0.into()), - spectrum.data()[7], - "Vector must be ordered" + + spectrum.freq_val_exact(451.0).val(); + } + + #[test] + #[should_panic] + fn test_spectrum_get_frequency_value_closest_panic_below_min() { + let spectrum_vector = vec![ + (0.0_f32, 5.0_f32), + (450.0, 200.0), + ]; + + let spectrum = spectrum_vector + .into_iter() + .map(|(fr, val)| (fr.into(), val.into())) + .collect::>(); + let spectrum = FrequencySpectrum::new( + spectrum, + 50.0, ); - assert_eq!(0.0, spectrum.min().val(), "min() must work"); - assert_eq!(200.0, spectrum.max().val(), "max() must work"); - assert_eq!(200.0 - 0.0, spectrum.range().val(), "range() must work"); - assert_eq!(78.125, spectrum.average().val(), "average() must work"); - assert_eq!( - (50 + 100) as f32 / 2.0, - spectrum.median().val(), - "median() must work" + spectrum.freq_val_closest(-1.0); + } + + #[test] + #[should_panic] + fn test_spectrum_get_frequency_value_closest_panic_below_max() { + let spectrum_vector = vec![ + (0.0_f32, 5.0_f32), + (450.0, 200.0), + ]; + + let spectrum = spectrum_vector + .into_iter() + .map(|(fr, val)| (fr.into(), val.into())) + .collect::>(); + let spectrum = FrequencySpectrum::new( + spectrum, + 50.0, ); + + spectrum.freq_val_closest(451.0); } } diff --git a/src/tests/mod.rs b/src/tests/mod.rs index 2251db4..67f9ad6 100644 --- a/src/tests/mod.rs +++ b/src/tests/mod.rs @@ -23,12 +23,11 @@ SOFTWARE. */ use crate::tests::sine::sine_wave_audio_data_multiple; use crate::windows::{blackman_harris_4term, blackman_harris_7term, hamming_window, hann_window}; -use crate::{samples_fft_to_spectrum, FrequencyLimit, SpectrumTotalScaleFunctionFactory}; +use crate::{samples_fft_to_spectrum, FrequencyLimit, ComplexSpectrumScalingFunction}; use alloc::boxed::Box; use alloc::vec::Vec; use audio_visualizer::spectrum::staticc::plotters_png_file::spectrum_static_plotters_png_visualize; use audio_visualizer::waveform::staticc::plotters_png_file::waveform_static_plotters_png_visualize; -use audio_visualizer::waveform::staticc::png_file::waveform_static_png_visualize; use audio_visualizer::Channels; /// Directory with test samples (e.g. mp3) can be found here. @@ -60,11 +59,11 @@ fn test_spectrum_and_visualize_sine_waves_50_1000_3777hz() { .collect::>(); // FFT frequency accuracy is: sample_rate / (N / 2) - // 44100/(16384/2) = 5.383Hz + // 44100/(4096/2) = 21.5Hz // get a window that we want to analyze - // 1/44100 * 16384 => 0.3715s - let window = &sine_audio[0..2048]; + // 1/44100 * 4096 => 0.0928s + let window = &sine_audio[0..4096]; let no_window = &window[..]; let hamming_window = hamming_window(no_window); @@ -118,6 +117,18 @@ fn test_spectrum_and_visualize_sine_waves_50_1000_3777hz() { false, ); + // test getters match spectrum + // we use Hann windowed spectrum because the accuracy is much better than + // with no window! + assert!(spectrum_hann_window.freq_val_exact(50.0).val() > 0.85); + assert!(spectrum_hann_window.freq_val_closest(50.0).1.val() > 0.85); + assert!(spectrum_hann_window.freq_val_exact(1000.0).val() > 0.85); + assert!(spectrum_hann_window.freq_val_closest(1000.0).1.val() > 0.85); + assert!(spectrum_hann_window.freq_val_exact(3777.0).val() > 0.85); + assert!(spectrum_hann_window.freq_val_closest(3777.0).1.val() > 0.85); + assert!(spectrum_hann_window.freq_val_exact(500.0).val() < 0.00001); + assert!(spectrum_hann_window.freq_val_closest(500.0).1.val() < 0.00001); + /*for (fr, vol) in spectrum.iter() { // you will experience inaccuracies here // TODO add further smoothing / noise reduction @@ -127,6 +138,6 @@ fn test_spectrum_and_visualize_sine_waves_50_1000_3777hz() { }*/ } -fn get_scale_to_one_fn_factory() -> SpectrumTotalScaleFunctionFactory { +fn get_scale_to_one_fn_factory() -> ComplexSpectrumScalingFunction { Box::new(move |_min: f32, max: f32, _average: f32, _median: f32| Box::new(move |x| x / max)) } diff --git a/src/tests/sine.rs b/src/tests/sine.rs index 7d5aee9..47a4277 100644 --- a/src/tests/sine.rs +++ b/src/tests/sine.rs @@ -67,7 +67,7 @@ pub fn sine_wave_audio_data_multiple( .collect:: f32>>>(); // How many samples to generate with each sine wave function - let sample_count = (sampling_rate as f32 * (duration_ms as f32 / 1000_f32)) as usize; + let sample_count = (sampling_rate as f32 * (duration_ms as f32 / 1000.0)) as usize; // Calculate the final sine wave let mut sine_wave = Vec::with_capacity(sample_count); diff --git a/src/windows.rs b/src/windows.rs index 3a25359..cfb4845 100644 --- a/src/windows.rs +++ b/src/windows.rs @@ -24,7 +24,7 @@ SOFTWARE. //! Several window functions which you can apply before doing the FFT. //! For more information: //! - https://en.wikipedia.org/wiki/Window_function -//! - https://www.youtube.com/watch?v=dCeHOf4cJE0 +//! - https://www.youtube.com/watch?v=dCeHOf4cJE0 (FFT and windowing by Texas Instruments) use alloc::vec::Vec; use core::f32::consts::PI;