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

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

/// Build the DCT-II ortho matrix of shape `[n_mfcc, n_mels]`.
///
/// `D[k, n] = c[k] * cos(pi * k * (2n+1) / (2 * n_mels))`
/// where `c[0] = sqrt(1/n_mels)`, `c[k>0] = sqrt(2/n_mels)`.
pub fn dct_matrix(n_mfcc: usize, n_mels: usize) -> Matrix {
    let scale0 = (1.0 / n_mels as f64).sqrt();
    let scale = (2.0 / n_mels as f64).sqrt();
    Matrix::from_fn(n_mfcc, n_mels, |k, n| {
        let s = if k == 0 { scale0 } else { scale };
        let angle = std::f64::consts::PI * k as f64 * (2 * n + 1) as f64 / (2.0 * n_mels as f64);
        s * angle.cos()
    })
}

/// Parameters for MFCC computation.
pub struct MfccParams {
    pub n_mfcc: usize,
    /// Reference power for power_to_db. None = use max(mel). Default: Some(1.0).
    /// - Some(1.0): matches librosa.feature.mfcc(y=y) default
    /// - None: matches librosa.power_to_db(mel, ref=np.max) + mfcc(S=...)
    pub ref_power: Option<f64>,
    pub amin: f64,
    pub top_db: Option<f64>,
}

impl Default for MfccParams {
    fn default() -> Self {
        Self {
            n_mfcc: 20,
            ref_power: Some(1.0),
            amin: 1e-10,
            top_db: Some(80.0),
        }
    }
}

/// Compute MFCCs from a power mel spectrogram.
///
/// Input:  `mel [n_mels, T]` — linear-scale power mel spectrogram
/// Output: `[n_mfcc, T]`
///
/// Steps: `power_to_db(mel)` → DCT-II ortho along axis=0 → keep first `n_mfcc` rows.
///
/// When `ref_power` is None, uses max(mel) as reference (matches
/// `librosa.power_to_db(mel, ref=np.max)` + `librosa.feature.mfcc(S=...)`).
pub fn mfcc(mel: &Matrix, params: &MfccParams) -> Matrix {
    let n_mels = mel.rows();
    let ref_val = params.ref_power.unwrap_or_else(|| {
        mel.as_slice()
            .iter()
            .cloned()
            .fold(f64::NEG_INFINITY, f64::max)
    });
    let log_mel_flat = power_to_db(mel.as_slice(), ref_val, params.amin, params.top_db);
    let log_mel = Matrix::from_vec(log_mel_flat, n_mels, mel.cols());

    let d = dct_matrix(params.n_mfcc, n_mels);
    d.matmul(&log_mel)
}

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

    #[test]
    fn test_dct_matrix_orthonormal() {
        // D @ D^T should be identity (up to floating point)
        let d = dct_matrix(20, 128);
        let dt = d.transpose();
        let eye = d.matmul(&dt);
        for k in 0..20 {
            for l in 0..20 {
                let expected = if k == l { 1.0 } else { 0.0 };
                let diff = (eye.get(k, l) - expected).abs();
                assert!(
                    diff < 1e-12,
                    "D @ D^T [{k},{l}] = {}, expected {expected}",
                    eye.get(k, l)
                );
            }
        }
    }

    #[test]
    fn test_mfcc_shape() {
        let mel = Matrix::from_fn(128, 50, |r, c| (r + c) as f64 * 0.01 + 1e-6);
        let result = mfcc(&mel, &MfccParams::default());
        assert_eq!(result.shape(), (20, 50));
    }
}