rosa 0.1.0

Audio feature extraction library matching librosa output
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

//! Spectral contrast matching `librosa.feature.spectral_contrast`.

use crate::dsp::db::power_to_db;
use crate::matrix::Matrix;

/// Compute spectral contrast from a magnitude spectrogram.
///
/// Matches `librosa.feature.spectral_contrast(S=S, sr=sr, freq=freq, fmin=fmin,
///   n_bands=n_bands, quantile=quantile, linear=linear)`.
///
/// Each frame is divided into octave sub-bands. For each sub-band,
/// contrast is the difference between peak and valley energy.
///
/// - `s`: magnitude spectrogram (rows = frequency bins, cols = frames)
/// - `freq`: center frequencies for each bin (length = s.rows())
/// - `fmin`: minimum frequency for first band (default: 200.0)
/// - `n_bands`: number of octave bands (default: 6)
/// - `quantile`: quantile for peaks/valleys (default: 0.02)
/// - `linear`: if true, return linear difference; if false, dB difference (default: false)
///
/// Returns a Matrix with shape (n_bands + 1, n_frames).
///
/// # Panics
///
/// Panics if `freq.len() != s.rows()`.
pub fn spectral_contrast(
    s: &Matrix,
    freq: &[f64],
    fmin: f64,
    n_bands: usize,
    quantile: f64,
    linear: bool,
) -> Matrix {
    assert_eq!(
        s.rows(),
        freq.len(),
        "freq length {} != spectrogram rows {}",
        freq.len(),
        s.rows()
    );

    let n_frames = s.cols();

    // Octave frequency boundaries: [0, fmin, 2*fmin, 4*fmin, ..., fmin*2^n_bands]
    let mut octa = vec![0.0; n_bands + 2];
    for i in 0..=n_bands {
        octa[i + 1] = fmin * 2.0_f64.powi(i as i32);
    }

    let mut valley = Matrix::zeros(n_bands + 1, n_frames);
    let mut peak = Matrix::zeros(n_bands + 1, n_frames);

    for k in 0..=n_bands {
        let f_low = octa[k];
        let f_high = octa[k + 1];

        // Find bins in current band: freq >= f_low && freq <= f_high
        let mut current_band: Vec<bool> = freq.iter().map(|&f| f >= f_low && f <= f_high).collect();

        // Find first and last index in band
        let first_idx = current_band.iter().position(|&b| b);
        let last_idx = current_band.iter().rposition(|&b| b);

        if let (Some(first), Some(_last)) = (first_idx, last_idx) {
            // Include previous bin for bands > 0
            if k > 0 && first > 0 {
                current_band[first - 1] = true;
            }

            // For last band, include all bins above the upper edge
            if k == n_bands
                && let Some(last) = last_idx
            {
                for val in &mut current_band[(last + 1)..] {
                    *val = true;
                }
            }
        }

        // Collect indices of active bins
        let active_indices: Vec<usize> = current_band
            .iter()
            .enumerate()
            .filter(|&(_, b)| *b)
            .map(|(i, _)| i)
            .collect();

        // For non-last bands, exclude the last bin (upper edge)
        let sub_indices: Vec<usize> = if k < n_bands && active_indices.len() > 1 {
            active_indices[..active_indices.len() - 1].to_vec()
        } else {
            active_indices.clone()
        };

        if sub_indices.is_empty() {
            continue;
        }

        // Number of bins used for quantile calculation
        // librosa: idx = round(quantile * sum(current_band)), max(idx, 1)
        let n_active: usize = current_band.iter().filter(|&&b| b).count();
        let idx = ((quantile * n_active as f64).round() as usize).max(1);

        // For each frame, sort the sub-band magnitudes and compute valley/peak
        for t in 0..n_frames {
            let mut band_vals: Vec<f64> = sub_indices.iter().map(|&i| s.get(i, t)).collect();
            band_vals.sort_by(|a, b| a.total_cmp(b));

            let n = band_vals.len();
            let take = idx.min(n);

            // Valley: mean of lowest `take` values
            let valley_val: f64 = band_vals[..take].iter().sum::<f64>() / take as f64;
            // Peak: mean of highest `take` values
            let peak_val: f64 = band_vals[n - take..].iter().sum::<f64>() / take as f64;

            *valley.get_mut(k, t) = valley_val;
            *peak.get_mut(k, t) = peak_val;
        }
    }

    if linear {
        peak.zip_map(&valley, |p, v| p - v)
    } else {
        let peak_db = power_to_db(peak.as_slice(), 1.0, 1e-10, Some(80.0));
        let valley_db = power_to_db(valley.as_slice(), 1.0, 1e-10, Some(80.0));
        let contrast_data: Vec<f64> = peak_db
            .iter()
            .zip(valley_db.iter())
            .map(|(&p, &v)| p - v)
            .collect();
        Matrix::from_vec(contrast_data, n_bands + 1, n_frames)
    }
}