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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227

//! Hybrid CQT and Pseudo CQT implementations.

use crate::matrix::Matrix;

use super::transform::{CqtParams, complex_magnitude, complex_matmul, cqt};
use super::wavelet::{WaveletLengthsParams, WaveletParams, wavelet, wavelet_lengths};

/// Pseudo CQT — single FFT size, returns magnitude only (no phase).
///
/// Matches `librosa.pseudo_cqt(...)`.
/// Returns Matrix of shape (n_bins, n_frames) — real-valued magnitude.
pub fn pseudo_cqt(y: &[f64], params: &CqtParams) -> Matrix {
    let freqs = crate::convert::cqt_frequencies(
        params.n_bins,
        params.fmin,
        params.bins_per_octave,
        params.tuning,
    );

    let wl_params = WaveletLengthsParams {
        sr: params.sr,
        window: params.window,
        filter_scale: params.filter_scale,
        gamma: Some(0.0),
    };
    let (lengths, _f_cutoff) = wavelet_lengths(&freqs, &wl_params);

    // Build wavelet filters
    let wp = WaveletParams {
        sr: params.sr,
        window: params.window,
        filter_scale: params.filter_scale,
        pad_fft: true,
        norm: params.norm,
        gamma: Some(0.0),
    };
    let (filters_re, filters_im, _) = wavelet(&freqs, &wp);

    // Determine n_fft: max of filter length and 2*hop_length, rounded to power of 2
    let min_fft = (2 * params.hop_length).next_power_of_two();
    let n_fft = filters_re.cols().max(min_fft).next_power_of_two();

    // FFT the filter basis
    let n_bins_fft = n_fft / 2 + 1;
    let n_freqs = freqs.len();

    let mut planner = realfft::RealFftPlanner::<f64>::new();
    let fft = planner.plan_fft_forward(n_fft);
    let mut input_buf = fft.make_input_vec();
    let mut output_buf = fft.make_output_vec();
    let mut scratch = fft.make_scratch_vec();

    let mut basis_re_data = vec![0.0; n_freqs * n_bins_fft];
    let mut basis_im_data = vec![0.0; n_freqs * n_bins_fft];

    for k in 0..n_freqs {
        let scale = lengths[k] / n_fft as f64;
        let filter_cols = filters_re.cols();
        let pad_start = if n_fft > filter_cols {
            (n_fft - filter_cols) / 2
        } else {
            0
        };
        let copy_len = filter_cols.min(n_fft);

        // FFT real part
        input_buf.fill(0.0);
        for j in 0..copy_len {
            input_buf[pad_start + j] = filters_re.get(k, j) * scale;
        }
        fft.process_with_scratch(&mut input_buf, &mut output_buf, &mut scratch)
            .expect("FFT failed");

        let mut fft_re_re = vec![0.0; n_bins_fft];
        let mut fft_re_im = vec![0.0; n_bins_fft];
        for j in 0..n_bins_fft {
            fft_re_re[j] = output_buf[j].re;
            fft_re_im[j] = output_buf[j].im;
        }

        // FFT imag part
        input_buf.fill(0.0);
        for j in 0..copy_len {
            input_buf[pad_start + j] = filters_im.get(k, j) * scale;
        }
        fft.process_with_scratch(&mut input_buf, &mut output_buf, &mut scratch)
            .expect("FFT failed");

        for j in 0..n_bins_fft {
            let fft_im_re = output_buf[j].re;
            let fft_im_im = output_buf[j].im;
            basis_re_data[k * n_bins_fft + j] = fft_re_re[j] - fft_im_im;
            basis_im_data[k * n_bins_fft + j] = fft_re_im[j] + fft_im_re;
        }
    }

    let basis_re = Matrix::from_vec(basis_re_data, n_freqs, n_bins_fft);
    let basis_im = Matrix::from_vec(basis_im_data, n_freqs, n_bins_fft);

    // STFT with boxcar window
    let stft_params = crate::stft::StftParams {
        n_fft,
        hop_length: params.hop_length,
        win_length: n_fft,
        center: true,
        window: Some(vec![1.0; n_fft]),
    };
    let (stft_re, stft_im) = crate::stft::stft_as_real_imag(y, &stft_params);

    // response = abs(basis @ STFT)
    let (resp_re, resp_im) = complex_matmul(&basis_re, &basis_im, &stft_re, &stft_im);
    let mut result = complex_magnitude(&resp_re, &resp_im);

    // Scale by 1/sqrt(filter_length)
    if params.scale {
        for k in 0..params.n_bins {
            let s = 1.0 / lengths[k].sqrt();
            for f in 0..result.cols() {
                *result.get_mut(k, f) *= s;
            }
        }
    }

    result
}

/// Hybrid CQT — pseudo for short filters, full CQT for long filters.
///
/// Matches `librosa.hybrid_cqt(...)`.
/// Returns (re, im) matrices of shape (n_bins, n_frames).
pub fn hybrid_cqt(y: &[f64], params: &CqtParams) -> (Matrix, Matrix) {
    let freqs = crate::convert::cqt_frequencies(
        params.n_bins,
        params.fmin,
        params.bins_per_octave,
        params.tuning,
    );

    let wl_params = WaveletLengthsParams {
        sr: params.sr,
        window: params.window,
        filter_scale: params.filter_scale,
        gamma: Some(0.0),
    };
    let (lengths, _) = wavelet_lengths(&freqs, &wl_params);

    // Find split point: bins where filter_len_pow2 < 2 * hop_length use pseudo
    let pseudo_threshold = 2 * params.hop_length;
    let mut n_full = params.n_bins;
    for k in (0..params.n_bins).rev() {
        let flen_pow2 = (lengths[k].ceil() as usize).next_power_of_two();
        if flen_pow2 < pseudo_threshold {
            n_full = k;
            break;
        }
    }

    if n_full == params.n_bins {
        // All bins use full CQT
        return cqt(y, params);
    }

    if n_full == 0 {
        // All bins use pseudo CQT
        let mag = pseudo_cqt(y, params);
        let im = Matrix::zeros(mag.rows(), mag.cols());
        return (mag, im);
    }

    // Split: lower bins (longer filters) use full CQT, upper bins use pseudo
    let full_params = CqtParams {
        sr: params.sr,
        hop_length: params.hop_length,
        fmin: params.fmin,
        n_bins: n_full,
        bins_per_octave: params.bins_per_octave,
        tuning: params.tuning,
        filter_scale: params.filter_scale,
        norm: params.norm,
        window: params.window,
        scale: params.scale,
    };
    let (full_re, full_im) = cqt(y, &full_params);

    // Pseudo CQT for upper bins
    let n_pseudo = params.n_bins - n_full;
    let pseudo_fmin = freqs[n_full];
    let pseudo_params = CqtParams {
        sr: params.sr,
        hop_length: params.hop_length,
        fmin: pseudo_fmin,
        n_bins: n_pseudo,
        bins_per_octave: params.bins_per_octave,
        tuning: 0.0, // tuning already incorporated in fmin
        filter_scale: params.filter_scale,
        norm: params.norm,
        window: params.window,
        scale: params.scale,
    };
    let pseudo_mag = pseudo_cqt(y, &pseudo_params);

    // Trim to same number of frames
    let n_frames = full_re.cols().min(pseudo_mag.cols());
    let mut out_re = vec![0.0; params.n_bins * n_frames];
    let mut out_im = vec![0.0; params.n_bins * n_frames];

    // Copy full CQT (lower bins)
    for k in 0..n_full {
        for f in 0..n_frames {
            out_re[k * n_frames + f] = full_re.get(k, f);
            out_im[k * n_frames + f] = full_im.get(k, f);
        }
    }

    // Copy pseudo CQT (upper bins)
    for k in 0..n_pseudo {
        for f in 0..n_frames {
            out_re[(n_full + k) * n_frames + f] = pseudo_mag.get(k, f);
            // im stays 0
        }
    }

    (
        Matrix::from_vec(out_re, params.n_bins, n_frames),
        Matrix::from_vec(out_im, params.n_bins, n_frames),
    )
}