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

//! STFT, ISTFT, mel filterbank, and mel spectrogram.
mod convert;
mod mel;
mod transform;

pub use convert::{fft_frequencies, hz_to_mel, mel_frequencies, mel_to_hz};
pub use mel::{MelFilterParams, MelNorm, MelSpectrogramParams, mel_filter, mel_spectrogram};
pub use transform::{IstftParams, StftParams, istft, spectrogram, stft, stft_as_real_imag};

use crate::matrix::Matrix;

/// Separate a complex spectrogram into magnitude and phase.
///
/// Matches `librosa.magphase`.
///
/// # Parameters
/// - `re`, `im`: Real and imaginary parts of a complex spectrogram (same shape)
/// - `power`: Exponent for magnitude (1.0 = amplitude, 2.0 = power)
///
/// # Returns
/// `(magnitude, phase)` where:
/// - `magnitude[i,j] = (re[i,j]² + im[i,j]²).sqrt().powf(power)`
/// - `phase[i,j] = atan2(im[i,j], re[i,j])`
///
/// For zero magnitudes, phase is 0.0.
pub fn magphase(re: &Matrix, im: &Matrix, power: f64) -> (Matrix, Matrix) {
    assert_eq!(re.shape(), im.shape());
    let (rows, cols) = re.shape();
    let mut mag = Matrix::zeros(rows, cols);
    let mut phase = Matrix::zeros(rows, cols);
    for r in 0..rows {
        for c in 0..cols {
            let re_val = re.get(r, c);
            let im_val = im.get(r, c);
            let m = (re_val * re_val + im_val * im_val).sqrt();
            *mag.get_mut(r, c) = if power == 1.0 { m } else { m.powf(power) };
            *phase.get_mut(r, c) = if m == 0.0 { 0.0 } else { im_val.atan2(re_val) };
        }
    }
    (mag, phase)
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::f64::consts::PI;

    #[test]
    fn test_magphase_known_values() {
        // 2x2 matrix with known complex values
        let re = Matrix::from_vec(vec![3.0, 0.0, -1.0, 1.0], 2, 2);
        let im = Matrix::from_vec(vec![4.0, 1.0, 0.0, -1.0], 2, 2);

        let (mag, phase) = magphase(&re, &im, 1.0);

        // mag = sqrt(re² + im²)
        assert!((mag.get(0, 0) - 5.0).abs() < 1e-12);
        assert!((mag.get(0, 1) - 1.0).abs() < 1e-12);
        assert!((mag.get(1, 0) - 1.0).abs() < 1e-12);
        assert!((mag.get(1, 1) - 2.0_f64.sqrt()).abs() < 1e-12);

        // phase = atan2(im, re)
        assert!((phase.get(0, 0) - (4.0_f64).atan2(3.0)).abs() < 1e-12);
        assert!((phase.get(0, 1) - PI / 2.0).abs() < 1e-12); // atan2(1, 0)
        assert!((phase.get(1, 0) - PI).abs() < 1e-12); // atan2(0, -1)
    }

    #[test]
    fn test_magphase_power() {
        let re = Matrix::from_vec(vec![3.0, 0.0], 1, 2);
        let im = Matrix::from_vec(vec![4.0, 0.0], 1, 2);

        let (mag, _phase) = magphase(&re, &im, 2.0);
        assert!((mag.get(0, 0) - 25.0).abs() < 1e-12); // 5^2
        assert!((mag.get(0, 1) - 0.0).abs() < 1e-12);
    }

    #[test]
    fn test_magphase_zero_phase() {
        let re = Matrix::from_vec(vec![0.0], 1, 1);
        let im = Matrix::from_vec(vec![0.0], 1, 1);
        let (_mag, phase) = magphase(&re, &im, 1.0);
        assert_eq!(phase.get(0, 0), 0.0);
    }

    #[test]
    fn test_magphase_roundtrip() {
        let re = Matrix::from_vec(vec![1.0, -2.0, 3.0, -4.0, 5.0, -6.0], 2, 3);
        let im = Matrix::from_vec(vec![7.0, -8.0, 9.0, -10.0, 11.0, -12.0], 2, 3);

        let (mag, phase) = magphase(&re, &im, 1.0);
        let (rows, cols) = re.shape();

        for r in 0..rows {
            for c in 0..cols {
                let m = mag.get(r, c);
                let p = phase.get(r, c);
                let re_recovered = m * p.cos();
                let im_recovered = m * p.sin();
                assert!(
                    (re_recovered - re.get(r, c)).abs() < 1e-10,
                    "re roundtrip failed at ({r},{c})"
                );
                assert!(
                    (im_recovered - im.get(r, c)).abs() < 1e-10,
                    "im roundtrip failed at ({r},{c})"
                );
            }
        }
    }
}