mecomp_analysis/

chroma.rs

//! Chroma feature extraction module.
//!
//! Contains functions to compute the chromagram of a song, and
//! then from this chromagram extract the song's tone and mode
//! (minor / major).
use crate::Feature;

use super::errors::{AnalysisError, AnalysisResult};
use super::utils::{Normalize, hz_to_octs_inplace, stft};
use bitvec::vec::BitVec;
use likely_stable::{LikelyResult, likely, unlikely};
use ndarray::{Array, Array1, Array2, Axis, Order, Zip, arr2, concatenate, s};
use ndarray_stats::QuantileExt;
use noisy_float::prelude::*;

/**
 * General object holding the chroma descriptor.
 *
 * Current chroma descriptors are interval features (see
 * <https://speech.di.uoa.gr/ICMC-SMC-2014/images/VOL_2/1461.pdf>).
 *
 * Contrary to the other descriptors that can be used with streaming
 * without consequences, this one performs better if the full song is used at
 * once.
 */
#[derive(Debug, Clone)]
#[allow(clippy::module_name_repetitions)]
pub struct ChromaDesc {
    sample_rate: u32,
    n_chroma: u32,
    values_chroma: Array2<f32>,
}

impl Normalize for ChromaDesc {
    const MAX_VALUE: Feature = 1.0;
    const MIN_VALUE: Feature = 0.;
}

impl ChromaDesc {
    pub const WINDOW_SIZE: usize = 8192;
    /// The theoretical maximum value for IC1-6 is each value at (1/2)².
    /// The reason is that `extract_interval_features` computes the product of the
    /// L1-normalized chroma vector (so all of its values are <= 1) by itself.
    /// The maximum value of this is all coordinates to 1/2 (since dyads will
    /// select three values). The maximum of this is then (1/2)², so the maximum of its
    /// L2 norm is this sqrt(2 * (1/2)²) ~= 0.62. However, real-life simulations shown
    /// that 0.25 is a good ceiling value (see tests).
    pub const MAX_L2_INTERVAL: f32 = 0.25;
    /// The theoretical maximum value for IC7-10 is each value at (1/3)³.
    /// The reason is that `extract_interval_features` computes the product of the
    /// L1-normalized chroma vector (so all of its values are <= 1) by itself.
    /// The maximum value of this is all coordinates to 1/3 (since triads will
    /// select three values). The maximum of this is then (1/3)³, so the maximum of its
    /// L2 norm is this sqrt(4 * (1/3)³) ~= 0.074. However, real-life simulations shown
    /// that 0.025 is a good ceiling value (see tests).
    pub const MAX_L2_TRIAD: f32 = 0.025;
    /// We are using atan2 to keep the ratio bounded.
    pub const MAX_TRIAD_INTERVAL_RATIO: f32 = std::f32::consts::FRAC_PI_2;
    #[must_use]
    #[inline]
    pub fn new(sample_rate: u32, n_chroma: u32) -> Self {
        Self {
            sample_rate,
            n_chroma,
            values_chroma: Array2::zeros((n_chroma as usize, 0)),
        }
    }

    /**
     * Compute and store the chroma of a signal.
     *
     * Passing a full song here once instead of streaming smaller parts of the
     * song will greatly improve accuracy.
     */
    #[allow(clippy::missing_errors_doc, clippy::missing_panics_doc)]
    #[inline]
    pub fn do_(&mut self, signal: &[f32]) -> AnalysisResult<()> {
        let stft = stft(signal, Self::WINDOW_SIZE, 2205);
        let tuning = estimate_tuning(self.sample_rate, &stft, Self::WINDOW_SIZE, 0.01, 12)?;
        let chroma = chroma_stft(
            self.sample_rate,
            &stft,
            Self::WINDOW_SIZE,
            self.n_chroma,
            tuning,
        )?;
        self.values_chroma = concatenate![Axis(1), self.values_chroma, chroma];
        Ok(())
    }

    /**
     * Get the song's interval features.
     *
     * Return the 6 pitch class set categories, as well as the major, minor,
     * diminished and augmented triads.
     *
     * See this paper <https://speech.di.uoa.gr/ICMC-SMC-2014/images/VOL_2/1461.pdf>
     * for more information ("Timbre-invariant Audio Features for Style Analysis of Classical
     * Music").
     */
    #[inline]
    pub fn get_value(&mut self) -> Vec<Feature> {
        let mut raw_features = chroma_interval_features(&self.values_chroma);
        let (mut interval_class, mut interval_class_mode) =
            raw_features.view_mut().split_at(Axis(0), 6);
        // Compute those two norms separately because the values for the IC1-6 and IC7-10 don't
        // have the same range.
        let l2_norm_interval_class = interval_class.dot(&interval_class).sqrt();
        let l2_norm_interval_class_mode = interval_class_mode.dot(&interval_class_mode).sqrt();
        if l2_norm_interval_class > 0. {
            interval_class /= l2_norm_interval_class;
        }
        if l2_norm_interval_class_mode > 0. {
            interval_class_mode /= l2_norm_interval_class_mode;
        }
        let mut features = raw_features.mapv_into_any(|x| self.normalize(x)).to_vec();

        let normalized_l2_norm_interval_class =
            (2. * l2_norm_interval_class / Self::MAX_L2_INTERVAL - 1.).min(1.);
        features.push(normalized_l2_norm_interval_class);
        let normalized_l2_norm_interval_class_mode =
            (2. * l2_norm_interval_class_mode / Self::MAX_L2_TRIAD - 1.).min(1.);
        features.push(normalized_l2_norm_interval_class_mode);
        let angle = (20. * l2_norm_interval_class_mode).atan2(l2_norm_interval_class + 1e-12_f32);
        let normalized_ratio = 2. * angle / Self::MAX_TRIAD_INTERVAL_RATIO - 1.;
        features.push(normalized_ratio);
        features
    }
}

// Functions below are Rust versions of python notebooks by AudioLabs Erlang
// (<https://www.audiolabs-erlangen.de/resources/MIR/FMP/C0/C0.html>)
#[allow(
    clippy::missing_errors_doc,
    clippy::missing_panics_doc,
    clippy::module_name_repetitions
)]
#[must_use]
#[inline]
pub fn chroma_interval_features(chroma: &Array2<f32>) -> Array1<f32> {
    let chroma = normalize_feature_sequence(&(chroma * 15.).exp());
    let templates = arr2(&[
        [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
        [1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
        [0, 0, 1, 0, 0, 0, 0, 1, 1, 0],
        [0, 0, 0, 1, 0, 0, 1, 0, 0, 1],
        [0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
        [0, 0, 0, 0, 0, 0, 1, 1, 0, 0],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    ]);
    let interval_feature_matrix = extract_interval_features(&chroma, &templates);
    interval_feature_matrix.mean_axis(Axis(1)).unwrap()
}

#[must_use]
#[inline]
pub fn extract_interval_features(chroma: &Array2<f32>, templates: &Array2<i32>) -> Array2<f32> {
    let n_templates = templates.shape()[1];
    let n_chroma = chroma.shape()[0]; // should be 12
    let n_cols = chroma.shape()[1];

    let chroma_t = chroma.t();
    let mut f_intervals: Array2<f32> = Array::zeros((n_templates, n_cols));

    // precompute active indices for all templates
    let active_indices: Vec<Vec<usize>> = templates
        .axis_iter(Axis(1))
        .map(|t| {
            t.iter()
                .enumerate()
                .filter_map(|(i, &v)| (v == 1).then_some(i))
                .collect()
        })
        .collect();

    // For each column in chroma, compute the product of values at shifted active indices
    for (col_idx, col) in chroma_t.rows().into_iter().enumerate() {
        for (tmpl_idx, indices) in active_indices.iter().enumerate() {
            let sum = (0..n_chroma)
                .map(|shift| {
                    indices
                        .iter()
                        .map(|&idx| col[(idx + shift) % n_chroma])
                        .product::<f32>()
                })
                .sum();
            f_intervals[(tmpl_idx, col_idx)] = sum;
        }
    }

    f_intervals
}

#[inline]
pub fn normalize_feature_sequence(feature: &Array2<f32>) -> Array2<f32> {
    let mut normalized_sequence = feature.to_owned();
    for mut column in normalized_sequence.columns_mut() {
        let sum: f32 = column.iter().copied().map(f32::abs).sum();
        if likely(sum >= 0.0001) {
            column /= sum;
        }
    }

    normalized_sequence
}

// All the functions below are more than heavily inspired from
// librosa"s code: https://github.com/librosa/librosa/blob/main/librosa/feature/spectral.py#L1165
// chroma(22050, n_fft=5, n_chroma=12)
//
// Could be precomputed, but it takes very little time to compute it
// on the fly compared to the rest of the functions, and we'd lose the
// possibility to tweak parameters.
#[allow(
    clippy::missing_errors_doc,
    clippy::missing_panics_doc,
    clippy::module_name_repetitions,
    clippy::missing_inline_in_public_items
)]
pub fn chroma_filter(
    sample_rate: u32,
    n_fft: usize,
    n_chroma: u32,
    tuning: f32,
) -> AnalysisResult<Array2<f32>> {
    let ctroct = 5.0;
    let octwidth = 2.;
    #[allow(clippy::cast_precision_loss)]
    let n_chroma2 = (n_chroma >> 1) as f32;
    #[allow(clippy::cast_precision_loss)]
    let n_chroma_float = n_chroma as f32;

    #[allow(clippy::cast_precision_loss)]
    let frequencies = Array::linspace(0., sample_rate as f32, n_fft + 1);

    let mut freq_bins = frequencies;
    hz_to_octs_inplace(&mut freq_bins, tuning, n_chroma);
    freq_bins *= n_chroma_float;
    freq_bins[0] = (1.5).mul_add(-n_chroma_float, freq_bins[1]);
    let mut binwidth_bins = Array::ones(freq_bins.raw_dim());
    binwidth_bins
        .slice_mut(s![0..freq_bins.len() - 1])
        .assign(&(&freq_bins.slice(s![1..]) - &freq_bins.slice(s![..-1])).mapv(|x| x.max(1.)));

    let mut d: Array2<f32> = Array::zeros((n_chroma as usize, (freq_bins).len()));
    for (idx, mut row) in d.rows_mut().into_iter().enumerate() {
        #[allow(clippy::cast_precision_loss)]
        row.fill(idx as f32);
    }

    d.zip_mut_with(&freq_bins, |d_elem, &fb| {
        let x = -*d_elem + fb;
        let x = n_chroma_float.mul_add(10., x + n_chroma2);
        *d_elem = x % n_chroma_float - n_chroma2;
    });
    d.zip_mut_with(&binwidth_bins, |d_elem, &bb| {
        let x = *d_elem / bb;
        *d_elem = (-2. * x * x).exp();
    });

    let mut wts = d;
    // Normalize by computing the l2-norm over the columns
    for mut col in wts.columns_mut() {
        let sum = col.pow2().sum().sqrt();
        if sum >= f32::MIN_POSITIVE {
            col /= sum;
        }
    }

    // Apply Gaussian tuning curve
    freq_bins = (-0.5 * ((freq_bins / n_chroma_float - ctroct) / octwidth).powi(2)).exp();

    wts *= &freq_bins;

    // np.roll(), np bro
    let mut b = Array2::zeros(wts.dim());
    b.slice_mut(s![-3.., ..]).assign(&wts.slice(s![..3, ..]));
    b.slice_mut(s![..-3, ..]).assign(&wts.slice(s![3.., ..]));

    wts = b;
    let non_aliased = 1 + n_fft / 2;
    Ok(wts.slice_move(s![.., ..non_aliased]))
}

#[allow(clippy::missing_errors_doc, clippy::missing_panics_doc)]
#[allow(clippy::missing_inline_in_public_items)]
pub fn pip_track(
    sample_rate: u32,
    spectrum: &Array2<f32>,
    n_fft: usize,
) -> AnalysisResult<(Vec<f32>, Vec<f32>)> {
    #[allow(clippy::cast_precision_loss)]
    let sample_rate_float = sample_rate as f32;
    let fmin = 150.0;
    let fmax = 4000.0.min(sample_rate_float / 2.0);
    let threshold = 0.1;

    let fft_freqs = Array::linspace(0., sample_rate_float / 2., 1 + n_fft / 2);

    let length = spectrum.shape()[1];

    let freq_mask = fft_freqs
        .iter()
        .map(|&f| (fmin <= f) && (f < fmax))
        .collect::<BitVec>();

    let ref_value = spectrum.map_axis(Axis(0), |x| {
        let first = *x.first().expect("empty spectrum axis");
        let max = x.fold(first, |acc, &elem| acc.max(elem));
        threshold * max
    });

    // compute number of taken columns and beginning / end indices
    let freq_mask_len = freq_mask.len();
    let (taken_columns, beginning, end) = freq_mask.iter().enumerate().fold(
        (0, freq_mask_len, 0),
        |(taken, beginning, end), (i, b)| {
            b.then(|| (taken + 1, beginning.min(i), end.max(i + 1)))
                .unwrap_or((taken, beginning, end))
        },
    );

    // Validate that a valid frequency range was found
    if beginning >= end {
        return Err(AnalysisError::AnalysisError(String::from(
            "in chroma: no valid frequency range found",
        )));
    }
    // There will be at most taken_columns * length elements in pitches / mags
    let mut pitches = Vec::with_capacity(taken_columns * length);
    let mut mags = Vec::with_capacity(taken_columns * length);

    let zipped = Zip::indexed(spectrum.slice(s![beginning..end - 3, ..]))
        .and(spectrum.slice(s![beginning + 1..end - 2, ..]))
        .and(spectrum.slice(s![beginning + 2..end - 1, ..]));

    // No need to handle the last column, since freq_mask[length - 1] is
    // always going to be `false` for 22.5kHz
    zipped.for_each(|(i, j), &before_elem, &elem, &after_elem| {
        if elem > ref_value[j] && after_elem <= elem && before_elem < elem {
            let avg = 0.5 * (after_elem - before_elem);
            let mut shift = (2.0).mul_add(elem, -after_elem - before_elem);
            if shift.abs() < f32::MIN_POSITIVE {
                shift += 1.;
            }
            shift = avg / shift;
            #[allow(clippy::cast_precision_loss)]
            pitches.push(((i + beginning + 1) as f32 + shift) * sample_rate_float / n_fft as f32);
            mags.push((0.5 * avg).mul_add(shift, elem));
        }
    });

    Ok((pitches, mags))
}

// Only use this with strictly positive `frequencies`.
#[allow(clippy::missing_errors_doc, clippy::missing_panics_doc)]
#[inline]
pub fn pitch_tuning(
    mut frequencies: Array1<f32>,
    resolution: f32,
    bins_per_octave: u32,
) -> AnalysisResult<f32> {
    if unlikely(frequencies.is_empty()) {
        return Ok(0.0);
    }
    hz_to_octs_inplace(&mut frequencies, 0.0, 12);
    #[allow(clippy::cast_precision_loss)]
    frequencies.mapv_inplace(|x| (bins_per_octave as f32 * x).fract());

    // Wrap values from [0,1) to [-0.5, 0.5), then shift back up to [0,1)
    frequencies.mapv_inplace(|x| if x >= 0.5 { x - 0.5 } else { x + 0.5 });

    #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
    let indexes = (frequencies / resolution).mapv(|x| x as usize);
    #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
    let mut counts: Array1<usize> = Array::zeros(resolution.recip() as usize);
    for &idx in &indexes {
        counts[idx] += 1;
    }
    let max_index = counts
        .argmax()
        .map_err_unlikely(|e| AnalysisError::AnalysisError(format!("in chroma: {e}")))?;

    // Return the bin with the most reoccurring frequency.
    #[allow(clippy::cast_precision_loss)]
    Ok((100. * resolution).mul_add(max_index as f32, -50.) / 100.)
}

#[allow(clippy::missing_errors_doc, clippy::missing_panics_doc)]
#[inline]
pub fn estimate_tuning(
    sample_rate: u32,
    spectrum: &Array2<f32>,
    n_fft: usize,
    resolution: f32,
    bins_per_octave: u32,
) -> AnalysisResult<f32> {
    let (pitch, mag) = pip_track(sample_rate, spectrum, n_fft)?;

    if unlikely(pitch.is_empty()) {
        return Ok(0.);
    }

    let (filtered_pitch, filtered_mag): (Vec<N32>, Vec<N32>) = pitch
        .iter()
        .zip(&mag)
        .filter(|&(&p, _)| p > 0.)
        .map(|(x, y)| (n32(*x), n32(*y)))
        .unzip();

    let mut mag_copy = filtered_mag.clone();
    let mid = mag_copy.len() / 2;
    let threshold = *mag_copy
        .select_nth_unstable_by(mid, |a, b| a.partial_cmp(b).unwrap())
        .1;

    let pitch = filtered_pitch
        .iter()
        .zip(&filtered_mag)
        .filter_map(
            |(&p, &m)| {
                if m >= threshold { Some(p.into()) } else { None }
            },
        )
        .collect::<Array1<f32>>();
    pitch_tuning(pitch, resolution, bins_per_octave)
}

#[allow(
    clippy::missing_errors_doc,
    clippy::missing_panics_doc,
    clippy::module_name_repetitions
)]
#[inline]
pub fn chroma_stft(
    sample_rate: u32,
    spectrum: &Array2<f32>, // shape: (window_length / 2 + 1, signal.len().div_ceil(hop_length))
    n_fft: usize,
    n_chroma: u32,
    tuning: f32,
) -> AnalysisResult<Array2<f32>> {
    let mut raw_chroma = chroma_filter(sample_rate, n_fft, n_chroma, tuning)?;

    raw_chroma = raw_chroma.dot(&spectrum.pow2());

    // We want to maximize cache locality, and are iterating over columns,
    // so let's make sure our array is in column-major order.
    raw_chroma = raw_chroma
        .to_shape((raw_chroma.dim(), Order::ColumnMajor))
        .map_err_unlikely(|_| {
            AnalysisError::AnalysisError(String::from("in chroma: failed to reorder array"))
        })?
        .to_owned();

    Zip::from(raw_chroma.columns_mut()).for_each(|mut row| {
        let sum = row.sum(); // we know that our values are positive, so no need to use abs
        if sum >= f32::MIN_POSITIVE {
            row /= sum;
        }
    });

    Ok(raw_chroma)
}

#[cfg(test)]
mod test {
    use super::*;
    use crate::{
        SAMPLE_RATE,
        decoder::{Decoder as _, MecompDecoder as Decoder},
        utils::stft,
    };
    use ndarray::{Array2, arr1, arr2};
    use ndarray_npy::ReadNpyExt as _;
    use std::{fs::File, path::Path};

    #[test]
    fn test_chroma_interval_features() {
        let file = File::open("data/chroma.npy").unwrap();
        let chroma = Array2::<f64>::read_npy(file).unwrap().mapv(|x| x as f32);
        let features = chroma_interval_features(&chroma);
        let expected_features = arr1(&[
            0.038_602_84,
            0.021_852_81,
            0.042_243_79,
            0.063_852_78,
            0.073_111_48,
            0.025_125_66,
            0.003_198_99,
            0.003_113_08,
            0.001_074_33,
            0.002_418_61,
        ]);
        for (expected, actual) in expected_features.iter().zip(&features) {
            // original test wanted 000_000_01
            assert!(
                0.000_000_1 > (expected - actual.abs()),
                "{expected} !~= {actual}"
            );
        }
    }

    #[test]
    fn test_extract_interval_features() {
        let file = File::open("data/chroma-interval.npy").unwrap();
        let chroma = Array2::<f64>::read_npy(file).unwrap().mapv(|x| x as f32);
        let templates = arr2(&[
            [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
            [1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 1, 0, 0, 0, 0, 1, 1, 0],
            [0, 0, 0, 1, 0, 0, 1, 0, 0, 1],
            [0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
            [0, 0, 0, 0, 0, 0, 1, 1, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        ]);

        let file = File::open("data/interval-feature-matrix.npy").unwrap();
        let expected_interval_features = Array2::<f64>::read_npy(file).unwrap().mapv(|x| x as f32);

        let interval_features = extract_interval_features(&chroma, &templates);
        for (expected, actual) in expected_interval_features
            .iter()
            .zip(interval_features.iter())
        {
            assert!(
                0.000_000_1 > (expected - actual).abs(),
                "{expected} !~= {actual}"
            );
        }
    }

    #[test]
    fn test_normalize_feature_sequence() {
        let array = arr2(&[[0.1, 0.3, 0.4], [1.1, 0.53, 1.01]]);
        let expected_array = arr2(&[
            [0.083_333_33, 0.361_445_78, 0.283_687_94],
            [0.916_666_67, 0.638_554_22, 0.716_312_06],
        ]);

        let normalized_array = normalize_feature_sequence(&array);

        assert!(!array.is_empty() && !expected_array.is_empty());

        for (expected, actual) in normalized_array.iter().zip(expected_array.iter()) {
            assert!(
                0.000_000_1 > (expected - actual).abs(),
                "{expected} !~= {actual}"
            );
        }
    }

    #[test]
    fn test_chroma_desc() {
        let song = Decoder::new()
            .unwrap()
            .decode(Path::new("data/s16_mono_22_5kHz.flac"))
            .unwrap();
        let mut chroma_desc = ChromaDesc::new(SAMPLE_RATE, 12);
        chroma_desc.do_(&song.samples).unwrap();
        let expected_values = [
            -0.342_925_13,
            -0.628_034_23,
            -0.280_950_96,
            0.086_864_59,
            0.244_460_82,
            -0.572_325_7,
            0.232_920_65,
            0.199_811_46,
            -0.585_944_06,
            -0.067_842_96,
        ];
        for (expected, actual) in expected_values.iter().zip(chroma_desc.get_value().iter()) {
            // original test wanted absolute error < 0.0000001
            let relative_error = (expected - actual).abs() / expected.abs();
            assert!(
                relative_error < 0.01,
                "relative error: {relative_error}, expected: {expected}, actual: {actual}"
            );
        }
    }

    #[test]
    fn test_chroma_stft_decode() {
        let signal = Decoder::new()
            .unwrap()
            .decode(Path::new("data/s16_mono_22_5kHz.flac"))
            .unwrap()
            .samples;
        let stft = stft(&signal, 8192, 2205);

        let file = File::open("data/chroma.npy").unwrap();
        let expected_chroma = Array2::<f64>::read_npy(file).unwrap().mapv(|x| x as f32);

        let chroma = chroma_stft(22050, &stft, 8192, 12, -0.049_999_999_999_999_99).unwrap();

        assert!(!chroma.is_empty() && !expected_chroma.is_empty());

        for (expected, actual) in expected_chroma.iter().zip(chroma.iter()) {
            // original test wanted absolute error < 0.0000001
            let relative_error = (expected - actual).abs() / expected.abs();
            assert!(
                relative_error < 0.01,
                "relative error: {relative_error}, expected: {expected}, actual: {actual}"
            );
        }
    }

    #[test]
    fn test_estimate_tuning() {
        let file = File::open("data/spectrum-chroma.npy").unwrap();
        let arr = Array2::<f64>::read_npy(file).unwrap().mapv(|x| x as f32);

        let tuning = estimate_tuning(22050, &arr, 2048, 0.01, 12).unwrap();
        assert!(
            0.000_001 > (-0.099_999_999_999_999_98 - tuning).abs(),
            "{tuning} !~= -0.09999999999999998"
        );
    }

    #[test]
    fn test_chroma_estimate_tuning_empty_fix() {
        assert!(0. == estimate_tuning(22050, &Array2::zeros((8192, 1)), 8192, 0.01, 12).unwrap());
    }

    #[test]
    fn test_estimate_tuning_decode() {
        let signal = Decoder::new()
            .unwrap()
            .decode(Path::new("data/s16_mono_22_5kHz.flac"))
            .unwrap()
            .samples;
        let stft = stft(&signal, 8192, 2205);

        let tuning = estimate_tuning(22050, &stft, 8192, 0.01, 12).unwrap();
        assert!(
            0.000_001 > (-0.049_999_999_999_999_99 - tuning).abs(),
            "{tuning} !~= -0.04999999999999999"
        );
    }

    #[test]
    fn test_pitch_tuning() {
        let file = File::open("data/pitch-tuning.npy").unwrap();
        let pitch = Array1::<f64>::read_npy(file).unwrap();
        #[allow(clippy::cast_possible_truncation)]
        let pitch = pitch.mapv(|x| x as f32);

        let tuned = pitch_tuning(pitch, 0.05, 12).unwrap();
        assert!(f32::EPSILON > (tuned + 0.1).abs(), "{tuned} != -0.1");
    }

    #[test]
    fn test_pitch_tuning_no_frequencies() {
        let frequencies = arr1(&[]);
        let tuned = pitch_tuning(frequencies, 0.05, 12).unwrap();
        assert!(f32::EPSILON > tuned.abs(), "{tuned} != 0");
    }

    #[test]
    fn test_pip_track() {
        let file = File::open("data/spectrum-chroma.npy").unwrap();
        let spectrum = Array2::<f64>::read_npy(file).unwrap().mapv(|x| x as f32);

        let mags_file = File::open("data/spectrum-chroma-mags.npy").unwrap();
        let expected_mags = Array1::<f64>::read_npy(mags_file)
            .unwrap()
            .mapv(|x| x as f32);

        let pitches_file = File::open("data/spectrum-chroma-pitches.npy").unwrap();
        let expected_pitches = Array1::<f64>::read_npy(pitches_file)
            .unwrap()
            .mapv(|x| x as f32);

        let (mut pitches, mut mags) = pip_track(22050, &spectrum, 2048).unwrap();
        pitches.sort_by(|a, b| a.partial_cmp(b).unwrap());
        mags.sort_by(|a, b| a.partial_cmp(b).unwrap());

        for (expected_pitches, actual_pitches) in expected_pitches.iter().zip(pitches.iter()) {
            // original test wanted 000_000_01
            assert!(
                0.001 > (expected_pitches - actual_pitches).abs(),
                "{expected_pitches} !~= {actual_pitches}"
            );
        }
        for (expected_mags, actual_mags) in expected_mags.iter().zip(mags.iter()) {
            // original test wanted 000_000_01
            assert!(
                0.001 > (expected_mags - actual_mags).abs(),
                "{expected_mags} !~= {actual_mags}"
            );
        }
    }

    #[test]
    fn test_chroma_filter() {
        let file = File::open("data/chroma-filter.npy").unwrap();
        let expected_filter = Array2::<f64>::read_npy(file).unwrap().mapv(|x| x as f32);

        let filter = chroma_filter(22050, 2048, 12, -0.1).unwrap();

        assert!(filter.iter().all(|&x| x > 0.));

        for (expected, actual) in expected_filter.iter().zip(filter.iter()) {
            // original test wanted 0.000_000_001
            assert!(
                0.000_1 > (expected - actual).abs(),
                "{expected} !~= {actual}"
            );
        }
    }

    #[rstest::rstest]
    // High 6 should be a major triad, 7 minor, 8 diminished and 9 augmented.
    #[case::major_triad("data/chroma/Cmaj.ogg", 6)]
    #[case::major_triad("data/chroma/Dmaj.ogg", 6)]
    #[case::minor_triad("data/chroma/Cmin.ogg", 7)]
    #[case::diminished_triad("data/chroma/Cdim.ogg", 8)]
    #[case::augmented_triad("data/chroma/Caug.ogg", 9)]
    fn test_end_result_triads(
        #[case] path: &str,
        #[case] expected_dominant_chroma_feature_index: usize,
    ) {
        let song = Decoder::new().unwrap().decode(Path::new(path)).unwrap();
        let mut chroma_desc = ChromaDesc::new(SAMPLE_RATE, 12);
        chroma_desc.do_(&song.samples).unwrap();
        let chroma_values = chroma_desc.get_value();

        let mut indices: Vec<usize> = (0..chroma_values.len()).collect();
        indices.sort_by(|&i, &j| chroma_values[j].partial_cmp(&chroma_values[i]).unwrap());
        assert!(indices[0] == expected_dominant_chroma_feature_index);
        for (i, v) in chroma_values.into_iter().enumerate() {
            if i >= 6 && i <= 10 {
                if i == expected_dominant_chroma_feature_index {
                    assert!(v > 0.8);
                } else {
                    assert!(v < 0.0);
                }
            }
        }
    }

    #[test]
    fn test_end_l2_norm_dyad() {
        let song = Decoder::new()
            .unwrap()
            .decode(Path::new("data/chroma/dyad_tritone_IC6.ogg"))
            .unwrap();
        let mut chroma_desc = ChromaDesc::new(SAMPLE_RATE, 12);
        chroma_desc.do_(&song.samples).unwrap();
        let chroma_values = chroma_desc.get_value();
        assert!(chroma_values[10] > 0.9);
    }

    #[test]
    fn test_end_l2_norm_mode() {
        let song = Decoder::new()
            .unwrap()
            .decode(Path::new("data/chroma/Cmaj_triads.ogg"))
            .unwrap();
        let mut chroma_desc = ChromaDesc::new(SAMPLE_RATE, 12);
        chroma_desc.do_(&song.samples).unwrap();
        let chroma_values = chroma_desc.get_value();
        assert!(chroma_values[11] > 0.9);
    }

    #[test]
    fn test_end_l2_norm_ratio() {
        let song = Decoder::new()
            .unwrap()
            .decode(Path::new("data/chroma/triad_aug_maximize_ratio.ogg"))
            .unwrap();
        let mut chroma_desc = ChromaDesc::new(SAMPLE_RATE, 12);
        chroma_desc.do_(&song.samples).unwrap();
        let chroma_values = chroma_desc.get_value();
        assert!(chroma_values[12] > 0.7);
    }

    #[rstest::rstest]
    // Test all 12 intervals.
    #[case::minor_second("data/chroma/minor_second.ogg", 0)]
    #[case::major_second("data/chroma/major_second.ogg", 1)]
    #[case::minor_third("data/chroma/minor_third.ogg", 2)]
    #[case::major_third("data/chroma/major_third.ogg", 3)]
    #[case::perfect_fourth("data/chroma/perfect_fourth.ogg", 4)]
    #[case::tritone("data/chroma/tritone.ogg", 5)]
    #[case::perfect_fifth("data/chroma/perfect_fifth.ogg", 4)]
    #[case::minor_sixth("data/chroma/minor_sixth.ogg", 3)]
    #[case::major_sixth("data/chroma/major_sixth.ogg", 2)]
    #[case::minor_seventh("data/chroma/minor_seventh.ogg", 1)]
    #[case::major_seventh("data/chroma/major_seventh.ogg", 0)]
    fn test_end_result_intervals(
        #[case] path: &str,
        #[case] expected_dominant_chroma_feature_index: usize,
    ) {
        let song = Decoder::new().unwrap().decode(Path::new(path)).unwrap();
        let mut chroma_desc = ChromaDesc::new(SAMPLE_RATE, 12);
        chroma_desc.do_(&song.samples).unwrap();
        let chroma_values = chroma_desc.get_value();

        let mut indices: Vec<usize> = (0..chroma_values.len()).collect();
        indices.sort_by(|&i, &j| chroma_values[j].partial_cmp(&chroma_values[i]).unwrap());
        assert_eq!(indices[0], expected_dominant_chroma_feature_index);
        for (i, v) in chroma_values.into_iter().enumerate() {
            if i < 6 {
                if i == expected_dominant_chroma_feature_index {
                    assert!(v > 0.9);
                } else {
                    assert!(v < 0.0);
                }
            }
        }
    }
}