dasp_rs/features/
harmonics.rs

1use ndarray::{Array2, ArrayView1};
2use num_complex::Complex;
3use rayon::prelude::*;
4use thiserror::Error;
5
6/// Error conditions for harmonic feature extraction and processing.
7///
8/// Enumerates specific failure modes in harmonic analysis and phase vocoding operations,
9/// providing detailed diagnostics for DSP pipeline debugging.
10#[derive(Error, Debug)]
11pub enum HarmonicsError {
12    /// Input arrays have mismatched lengths (e.g., amplitudes vs. frequency bins).
13    #[error("Length mismatch: {0} vs {1}")]
14    LengthMismatch(usize, usize),
15
16    /// Invalid input parameter (e.g., empty array, negative rate).
17    #[error("Invalid input: {0}")]
18    InvalidInput(String),
19
20    /// Numerical computation failure (e.g., division by zero, overflow).
21    #[error("Computation failed: {0}")]
22    ComputationFailed(String),
23}
24
25/// Interpolates harmonic amplitudes across frequency bins in parallel.
26///
27/// Performs linear interpolation of amplitude values at harmonic frequencies derived
28/// from a fundamental frequency grid, leveraging parallel processing for efficiency.
29///
30/// # Parameters
31/// - `x`: Amplitude spectrum as a slice of `f32` (frequency bin values).
32/// - `freqs`: Frequency bins corresponding to `x` (monotonically increasing).
33/// - `harmonics`: Harmonic multipliers (e.g., `[1.0, 2.0, 3.0]` for first three harmonics).
34///
35/// # Returns
36/// - `Ok(Array2<f32>)`: Interpolated amplitudes, shape `(n_harmonics, n_bins)`.
37/// - `Err(HarmonicsError)`: Failure due to length mismatch or invalid input.
38///
39/// # Constraints
40/// - `x.len() == freqs.len()`; enforced via error return.
41/// - `freqs` must be sorted in ascending order for correct interpolation.
42/// - Harmonic frequencies exceeding `freqs.last()` are clamped to zero.
43pub fn interp_harmonics(x: &[f32], freqs: &[f32], harmonics: &[f32]) -> Result<Array2<f32>, HarmonicsError> {
44    if x.len() != freqs.len() {
45        return Err(HarmonicsError::LengthMismatch(x.len(), freqs.len()));
46    }
47    if x.is_empty() {
48        return Err(HarmonicsError::InvalidInput("Amplitude spectrum is empty".to_string()));
49    }
50    if !freqs.windows(2).all(|w| w[0] <= w[1]) {
51        return Err(HarmonicsError::InvalidInput("Frequency bins must be sorted".to_string()));
52    }
53
54    let n_bins = freqs.len();
55    let n_harmonics = harmonics.len();
56    let mut result = Array2::zeros((n_harmonics, n_bins));
57
58    result.axis_iter_mut(ndarray::Axis(0))
59        .into_par_iter()
60        .enumerate()
61        .for_each(|(h_idx, mut row)| {
62            let h = harmonics[h_idx];
63            for (bin, &f) in freqs.iter().enumerate() {
64                let target_freq = f * h;
65                if target_freq <= *freqs.last().unwrap() {
66                    let left_idx = freqs.binary_search_by(|&x| x.partial_cmp(&target_freq).unwrap())
67                        .unwrap_or_else(|e| e.saturating_sub(1));
68                    let left_idx = left_idx.min(n_bins - 2);
69                    let right_idx = left_idx + 1;
70                    let left_freq = freqs[left_idx];
71                    let right_freq = freqs[right_idx];
72                    let alpha = if left_freq == right_freq {
73                        0.0
74                    } else {
75                        (target_freq - left_freq) / (right_freq - left_freq)
76                    };
77                    row[bin] = x[left_idx] * (1.0 - alpha) + x[right_idx] * alpha;
78                }
79            }
80        });
81
82    Ok(result)
83}
84
85/// Computes a salience map from a spectrogram via harmonic summation.
86///
87/// Sums weighted harmonic contributions across frequency bins and frames, producing
88/// a salience map for pitch detection or harmonic analysis. Parallelized over frames.
89///
90/// # Parameters
91/// - `s`: Spectrogram as `Array2<f32>`, shape `(n_bins, n_frames)`.
92/// - `freqs`: Frequency bins corresponding to spectrogram rows (monotonically increasing).
93/// - `harmonics`: Harmonic multipliers (e.g., `[1.0, 2.0]`).
94/// - `weights`: Optional harmonic weights; defaults to uniform `1.0` if `None`.
95///
96/// # Returns
97/// - `Ok(Array2<f32>)`: Salience map, shape `(n_bins, n_frames)`.
98/// - `Err(HarmonicsError)`: Failure due to dimension mismatch or invalid input.
99///
100/// # Constraints
101/// - `s.shape()[0] == freqs.len()`.
102/// - `weights.len() == harmonics.len()` if provided.
103pub fn salience(s: &Array2<f32>, freqs: &[f32], harmonics: &[f32], weights: Option<&[f32]>) -> Result<Array2<f32>, HarmonicsError> {
104    let n_bins = s.shape()[0];
105    let n_frames = s.shape()[1];
106    if n_bins != freqs.len() {
107        return Err(HarmonicsError::LengthMismatch(n_bins, freqs.len()));
108    }
109    if n_bins == 0 || n_frames == 0 {
110        return Err(HarmonicsError::InvalidInput("Spectrogram is empty".to_string()));
111    }
112    if !freqs.windows(2).all(|w| w[0] <= w[1]) {
113        return Err(HarmonicsError::InvalidInput("Frequency bins must be sorted".to_string()));
114    }
115
116    let n_harmonics = harmonics.len();
117    let default_weights = vec![1.0; n_harmonics];
118    let weights = weights.unwrap_or(&default_weights);
119    if weights.len() != n_harmonics {
120        return Err(HarmonicsError::LengthMismatch(weights.len(), n_harmonics));
121    }
122
123    let mut salience_map = Array2::zeros((n_bins, n_frames));
124    salience_map.axis_iter_mut(ndarray::Axis(1))
125        .into_par_iter()
126        .enumerate()
127        .for_each(|(frame, mut col)| {
128            for (bin, &f) in freqs.iter().enumerate() {
129                let mut total = 0.0;
130                for (h_idx, &h) in harmonics.iter().enumerate() {
131                    let harmonic_freq = f * h;
132                    if harmonic_freq <= freqs[n_bins - 1] {
133                        let left_idx = freqs.binary_search_by(|&x| x.partial_cmp(&harmonic_freq).unwrap())
134                            .unwrap_or_else(|e| e.saturating_sub(1));
135                        let left_idx = left_idx.min(n_bins - 2);
136                        let right_idx = left_idx + 1;
137                        let alpha = (harmonic_freq - freqs[left_idx]) / (freqs[right_idx] - freqs[left_idx]);
138                        let interp = s[[left_idx, frame]] * (1.0 - alpha) + s[[right_idx, frame]] * alpha;
139                        total += interp * weights[h_idx];
140                    }
141                }
142                col[bin] = total;
143            }
144        });
145
146    Ok(salience_map)
147}
148
149/// Extracts harmonic amplitudes from time-varying fundamental frequencies.
150///
151/// Interpolates amplitudes at harmonic frequencies based on frame-wise `f0` values,
152/// parallelized across frames for performance.
153///
154/// # Parameters
155/// - `x`: Amplitude spectrum as a slice of `f32`.
156/// - `f0`: Fundamental frequencies per frame.
157/// - `freqs`: Frequency bins corresponding to `x` (monotonically increasing).
158/// - `harmonics`: Harmonic multipliers.
159///
160/// # Returns
161/// - `Ok(Array2<f32>)`: Harmonic amplitudes, shape `(n_harmonics, n_frames)`.
162/// - `Err(HarmonicsError)`: Failure due to length mismatch or invalid input.
163///
164/// # Constraints
165/// - `x.len() == freqs.len()`.
166/// - `f0` must be non-empty.
167pub fn f0_harmonics(x: &[f32], f0: &[f32], freqs: &[f32], harmonics: &[f32]) -> Result<Array2<f32>, HarmonicsError> {
168    if x.len() != freqs.len() {
169        return Err(HarmonicsError::LengthMismatch(x.len(), freqs.len()));
170    }
171    if x.is_empty() || f0.is_empty() {
172        return Err(HarmonicsError::InvalidInput("Input arrays cannot be empty".to_string()));
173    }
174    if !freqs.windows(2).all(|w| w[0] <= w[1]) {
175        return Err(HarmonicsError::InvalidInput("Frequency bins must be sorted".to_string()));
176    }
177
178    let n_frames = f0.len();
179    let n_harmonics = harmonics.len();
180    let mut result = Array2::zeros((n_harmonics, n_frames));
181
182    result.axis_iter_mut(ndarray::Axis(1))
183        .into_par_iter()
184        .enumerate()
185        .for_each(|(frame, mut col)| {
186            let fund = f0[frame];
187            for (h_idx, &h) in harmonics.iter().enumerate() {
188                let target_freq = fund * h;
189                if target_freq <= freqs[freqs.len() - 1] {
190                    let left_idx = freqs.binary_search_by(|&x| x.partial_cmp(&target_freq).unwrap())
191                        .unwrap_or_else(|e| e.saturating_sub(1));
192                    let left_idx = left_idx.min(freqs.len() - 2);
193                    let right_idx = left_idx + 1;
194                    let alpha = (target_freq - freqs[left_idx]) / (freqs[right_idx] - freqs[left_idx]);
195                    col[h_idx] = x[left_idx] * (1.0 - alpha) + x[right_idx] * alpha;
196                }
197            }
198        });
199
200    Ok(result)
201}
202
203/// Performs phase vocoding for time-scale modification of a complex spectrogram.
204///
205/// Adjusts the temporal resolution of a spectrogram while preserving frequency content,
206/// implementing phase unwrapping and accumulation for coherent resynthesis.
207///
208/// # Parameters
209/// - `d`: Complex spectrogram as `Array2<Complex<f32>>`, shape `(n_bins, n_frames)`.
210/// - `rate`: Time stretching factor (>1 stretches, <1 compresses).
211/// - `hop_length`: Optional hop length between frames; defaults to `n_fft / 4`.
212/// - `n_fft`: Optional FFT size; defaults to `(n_bins - 1) * 2`.
213///
214/// # Returns
215/// - `Ok(Array2<Complex<f32>>)`: Time-scaled spectrogram with adjusted frame count.
216/// - `Err(HarmonicsError)`: Failure due to invalid rate or dimensions.
217///
218/// # Constraints
219/// - `rate > 0.0`.
220/// - `hop_length > 0` if provided.
221pub fn phase_vocoder(
222    d: &Array2<Complex<f32>>,
223    rate: f32,
224    hop_length: Option<usize>,
225    n_fft: Option<usize>,
226) -> Result<Array2<Complex<f32>>, HarmonicsError> {
227    if rate <= 0.0 {
228        return Err(HarmonicsError::InvalidInput("Rate must be positive".to_string()));
229    }
230    let n_bins = d.shape()[0];
231    let orig_frames = d.shape()[1];
232    if n_bins == 0 || orig_frames == 0 {
233        return Err(HarmonicsError::InvalidInput("Spectrogram is empty".to_string()));
234    }
235
236    let n_fft = n_fft.unwrap_or((n_bins - 1) * 2);
237    let hop = hop_length.unwrap_or(n_fft / 4);
238    if hop == 0 {
239        return Err(HarmonicsError::InvalidInput("Hop length must be positive".to_string()));
240    }
241
242    let new_frames = ((orig_frames as f32 * hop as f32) / rate / hop as f32).ceil() as usize;
243    let mut output = Array2::zeros((n_bins, new_frames));
244    let mut phase_acc = Array2::zeros((n_bins, 1));
245
246    let omega = (0..n_bins).map(|k| 2.0 * std::f32::consts::PI * k as f32 / n_fft as f32).collect::<Vec<f32>>();
247
248    for t in 0..new_frames {
249        let orig_t = (t as f32 * rate * hop as f32 / hop as f32) as usize;
250        let orig_t_next = ((t + 1) as f32 * rate * hop as f32 / hop as f32) as usize;
251
252        if orig_t >= orig_frames || orig_t_next >= orig_frames {
253            continue;
254        }
255
256        let mag = d.column(orig_t).mapv(|c| c.norm());
257        let phase = d.column(orig_t).mapv(|c| c.arg());
258        let phase_next = d.column(orig_t_next).mapv(|c| c.arg());
259        let delta_phase = phase_next - phase - &ArrayView1::from(&omega) * hop as f32;
260
261        let delta_phase_unwrapped: Vec<f32> = delta_phase.iter()
262            .map(|&dp| dp - 2.0 * std::f32::consts::PI * (dp / (2.0 * std::f32::consts::PI)).round())
263            .collect();
264
265        let phase_advance = ArrayView1::from(&delta_phase_unwrapped).mapv(|x| x / hop as f32 * (hop as f32 * rate));
266        phase_acc = phase_acc + phase_advance;
267
268        output.column_mut(t).assign(&mag.mapv(|m| Complex::from_polar(m, phase_acc[[0, 0]])));
269    }
270
271    Ok(output)
272}
273
274#[cfg(test)]
275mod tests {
276    use super::*;
277    use ndarray::array;
278
279    #[test]
280    fn test_interp_harmonics() {
281        let x = vec![0.1, 0.2, 0.3, 0.4];
282        let freqs = vec![0.0, 100.0, 200.0, 300.0];
283        let harmonics = vec![1.0, 2.0];
284        let result = interp_harmonics(&x, &freqs, &harmonics).unwrap();
285        assert_eq!(result.shape(), &[2, 4]);
286        assert_eq!(result[[0, 0]], 0.1);
287        assert_eq!(result[[1, 0]], 0.1);
288        assert_eq!(result[[0, 1]], 0.2);
289        assert_eq!(result[[1, 1]], 0.3);
290    }
291
292    #[test]
293    fn test_interp_harmonics_mismatch() {
294        let x = vec![0.1, 0.2];
295        let freqs = vec![0.0, 100.0, 200.0];
296        let harmonics = vec![1.0];
297        let result = interp_harmonics(&x, &freqs, &harmonics);
298        assert!(matches!(result, Err(HarmonicsError::LengthMismatch(2, 3))));
299    }
300
301    #[test]
302    fn test_salience() {
303        let s = array![[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9], [1.0, 1.1, 1.2]];
304        let freqs = vec![0.0, 100.0, 200.0, 300.0];
305        let harmonics = vec![1.0, 2.0];
306        let weights: Option<&[f32]> = Some(&[1.0, 0.5]);
307        let result = salience(&s, &freqs, &harmonics, weights).unwrap();
308        assert_eq!(result.shape(), &[4, 3]);
309        assert_eq!(result[[0, 0]], 0.1 + 0.1 * 0.5);
310        assert_eq!(result[[1, 0]], 0.4 + 0.35);
311    }
312
313    #[test]
314    fn test_salience_weight_mismatch() {
315        let s = array![[0.1, 0.2], [0.3, 0.4]];
316        let freqs = vec![0.0, 100.0];
317        let harmonics = vec![1.0, 2.0];
318        let weights: Option<&[f32]> = Some(&[1.0]);
319        let result = salience(&s, &freqs, &harmonics, weights);
320        assert!(matches!(result, Err(HarmonicsError::LengthMismatch(1, 2))));
321    }
322
323    #[test]
324    fn test_f0_harmonics() {
325        let x = vec![0.1, 0.2, 0.3, 0.4];
326        let f0 = vec![100.0, 150.0];
327        let freqs = vec![0.0, 100.0, 200.0, 300.0];
328        let harmonics = vec![1.0, 2.0];
329        let result = f0_harmonics(&x, &f0, &freqs, &harmonics).unwrap();
330        assert_eq!(result.shape(), &[2, 2]);
331        assert_eq!(result[[0, 0]], 0.2);
332        assert_eq!(result[[1, 0]], 0.3);
333        assert_eq!(result[[0, 1]], 0.25);
334    }
335
336    #[test]
337    fn test_f0_harmonics_empty() {
338        let x = vec![];
339        let f0 = vec![100.0];
340        let freqs = vec![];
341        let harmonics = vec![1.0];
342        let result = f0_harmonics(&x, &f0, &freqs, &harmonics);
343        assert!(matches!(result, Err(HarmonicsError::InvalidInput(_))));
344    }
345
346    #[test]
347    fn test_phase_vocoder() {
348        let d = array![
349            [Complex::new(1.0, 0.0), Complex::new(2.0, 0.0), Complex::new(3.0, 0.0), Complex::new(4.0, 0.0)],
350            [Complex::new(5.0, 0.0), Complex::new(6.0, 0.0), Complex::new(7.0, 0.0), Complex::new(8.0, 0.0)],
351        ];
352        let result = phase_vocoder(&d, 0.5, Some(1), Some(4)).unwrap();
353        assert_eq!(result.shape(), &[2, 8]);
354        assert_eq!(result[[0, 0]].norm(), 1.0);
355    }
356
357    #[test]
358    fn test_phase_vocoder_invalid_rate() {
359        let d = array![[Complex::new(1.0, 0.0)]];
360        let result = phase_vocoder(&d, 0.0, None, None);
361        assert!(matches!(result, Err(HarmonicsError::InvalidInput(_))));
362    }
363}
dasp_rs/features/harmonics.rs

dasp_rs/features/
harmonics.rs
