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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159

//! Linear Predictive Coding using Burg's method.
//!
//! Matches `librosa.lpc(y, order=order)`.
//!
//! ## Numerical conditioning
//!
//! For pure sinusoidal signals, order 2 captures essentially all signal energy
//! (a\[1\] ≈ -2cos(ω), a\[2\] ≈ 1). After order 2, the denominator `den` drops
//! by ~6 orders of magnitude, making higher orders numerically sensitive to
//! floating-point accumulation differences. Librosa's numba JIT uses SIMD and
//! potentially FMA instructions that produce slightly different accumulation
//! than Rust's sequential operations. A pure-Python (non-numba) reimplementation
//! also diverges from librosa by ~7e-4, confirming this is a compilation-level
//! numerical difference, not an algorithmic bug. Well-conditioned signals (e.g.
//! noise) match at 1e-10; ill-conditioned pure tones diverge ~3e-4 (order 4)
//! to ~9e-3 (order 8).

/// Compute LP (Linear Prediction) coefficients using Burg's method.
///
/// Returns AR coefficients of length `order + 1` where `a[0] = 1.0`.
///
/// Matches `librosa.lpc(y, order=order)`.
///
/// # Panics
/// - If `order < 1`
/// - If `y.len() <= order`
pub fn lpc(y: &[f64], order: usize) -> Vec<f64> {
    let n = y.len();
    assert!(order >= 1, "lpc: order must be >= 1, got {order}");
    assert!(
        n > order,
        "lpc: signal length ({n}) must be > order ({order})"
    );

    // Epsilon for numerical stability, matching librosa's `util.tiny(den)` which
    // returns `np.finfo(float64).tiny`. Librosa divides by `den[0] + epsilon`
    // (line 1083 in audio.py) and comments that the `den <= 0` early-exit check
    // "can be removed if we keep the epsilon bias". We follow the same approach:
    // no early break on `den <= 0`, rely on epsilon to prevent division by zero.
    let epsilon: f64 = f64::MIN_POSITIVE; // ~2.2e-308

    // Librosa initializes BOTH arrays with [0] = 1.0 via `ar_coeffs_prev = ar_coeffs.copy()`.
    // The Levinson-Durbin recursion at j = i+1 requires `ar_coeffs_prev[0] = 1.0` to copy
    // the reflection coefficient into the new position. In librosa, the reference swap
    // `ar_coeffs_prev, ar_coeffs = ar_coeffs, ar_coeffs_prev` alternates which physical
    // array is "prev", but since both started with [0] = 1 and the j loop never touches
    // index 0, the invariant holds. Here we only init one array but reset ar_coeffs[0] = 1.0
    // inside the loop body (after the swap) to achieve the same effect.
    let mut ar_coeffs = vec![0.0; order + 1];
    let mut ar_coeffs_prev = vec![0.0; order + 1];
    ar_coeffs[0] = 1.0;

    // Initialize forward and backward prediction errors
    let mut forward_error: Vec<f64> = y[1..].to_vec();
    let mut backward_error: Vec<f64> = y[..n - 1].to_vec();

    // Initial denominator (DEN_M from eqn 16 of Marple)
    let mut den: f64 = forward_error.iter().map(|x| x * x).sum::<f64>()
        + backward_error.iter().map(|x| x * x).sum::<f64>();

    for i in 0..order {
        // Reflection coefficient (Marple eq. 15)
        // Use epsilon bias instead of breaking on den <= 0 (matches librosa)
        let reflect_numer: f64 = forward_error
            .iter()
            .zip(backward_error.iter())
            .map(|(f, b)| b * f)
            .sum::<f64>();
        let reflect = -2.0 * reflect_numer / (den + epsilon);

        // Levinson-Durbin update (eq. 5).
        // The swap mirrors librosa's `ar_coeffs_prev, ar_coeffs = ar_coeffs, ar_coeffs_prev`
        // which swaps array references in numba (not values). We swap Vec handles, same effect.
        std::mem::swap(&mut ar_coeffs, &mut ar_coeffs_prev);
        for j in 1..=i + 1 {
            // When j == i+1: ar_coeffs_prev[0] must be 1.0 for the reflection coefficient
            // to be copied into position. See the initialization comment above.
            ar_coeffs[j] = ar_coeffs_prev[j] + reflect * ar_coeffs_prev[i + 1 - j];
        }
        // Reset [0] = 1.0 so that when this array becomes ar_coeffs_prev in the next
        // iteration's swap, ar_coeffs_prev[0] = 1.0 is guaranteed.
        ar_coeffs[0] = 1.0;

        // Update prediction errors (eqs. 13-14).
        // Librosa uses `fwd_pred_error_tmp = fwd_pred_error` (a reference, not copy) then
        // creates new arrays via `fwd + reflect_coeff * bwd` (numpy `+` allocates). We
        // achieve the same by computing both into fresh arrays from the originals.
        // Note: in librosa, `reflect_coeff` has shape (1,) — not (order+1,) — so it
        // broadcasts correctly across the error arrays. This is because the caller sets
        // `shape[0] = 1` before creating reflect_coeff (line 1031 of audio.py).
        let fe_len = forward_error.len();
        let mut new_forward = vec![0.0; fe_len];
        let mut new_backward = vec![0.0; fe_len];
        for k in 0..fe_len {
            new_forward[k] = forward_error[k] + reflect * backward_error[k];
            new_backward[k] = backward_error[k] + reflect * forward_error[k];
        }

        // Update denominator (eq. 17).
        // Uses the UPDATED (but not yet shifted) errors, matching librosa lines 1121-1122.
        // For ill-conditioned signals (e.g. pure tone after order 2), den can go negative
        // here due to floating-point rounding. We don't break — the epsilon bias on the
        // next iteration's division handles it gracefully.
        let q = 1.0 - reflect * reflect;
        den = q * den
            - new_backward[fe_len - 1] * new_backward[fe_len - 1]
            - new_forward[0] * new_forward[0];

        // Shift: drop first element of forward, last of backward
        if fe_len > 1 {
            forward_error = new_forward[1..].to_vec();
            backward_error = new_backward[..fe_len - 1].to_vec();
        } else {
            break;
        }
    }

    // Ensure a[0] = 1.0 in the result
    ar_coeffs[0] = 1.0;
    ar_coeffs
}

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

    #[test]
    fn test_lpc_order_1() {
        // Simple constant signal: LPC should give a[0]=1, a[1] ≈ -1
        let y: Vec<f64> = vec![1.0; 100];
        let coeffs = lpc(&y, 1);
        assert_eq!(coeffs.len(), 2);
        assert!((coeffs[0] - 1.0).abs() < 1e-15);
        // For constant signal, reflection coefficient is -1
        assert!((coeffs[1] - (-1.0)).abs() < 1e-12);
    }

    #[test]
    fn test_lpc_returns_correct_length() {
        let y: Vec<f64> = (0..100).map(|i| (i as f64 * 0.1).sin()).collect();
        let coeffs = lpc(&y, 4);
        assert_eq!(coeffs.len(), 5);
        assert!((coeffs[0] - 1.0).abs() < 1e-15);
    }

    #[test]
    #[should_panic]
    fn test_lpc_panics_if_order_zero() {
        let y = vec![1.0; 10];
        lpc(&y, 0);
    }

    #[test]
    #[should_panic]
    fn test_lpc_panics_if_signal_too_short() {
        let y = vec![1.0; 3];
        lpc(&y, 5);
    }
}