audio_lab/dsp/
frequency.rs

1use crate::signal::FreqSignal;
2use crate::signal::data::RealData;
3use ndarray::{Array1, Zip};
4
5#[derive(Debug, thiserror::Error)]
6pub enum FrequencyError {
7    #[error("group delay requires at least two frequency bins")]
8    TooFewFrequencyBins,
9}
10
11pub fn group_delay(signal: &FreqSignal) -> Result<RealData, FrequencyError> {
12    // The frequency axis is an `XAxis`, which is always strictly increasing,
13    // so only the length needs checking.
14    if signal.num_freq_bins() < 2 {
15        return Err(FrequencyError::TooFewFrequencyBins);
16    }
17
18    // Weight h[n] by n in place: FFT(n·h[n]) / FFT(h[n]) gives the group delay.
19    let mut weighted_time = signal.to_time();
20    let ramp = Array1::from_iter((0..weighted_time.num_time_steps()).map(|n| n as f64));
21    for mut channel in weighted_time.channel_iter_mut() {
22        channel *= &ramp;
23    }
24    let weighted_freq = weighted_time.into_freq();
25
26    // Re(N / D) = Re(N · conj(D)) / |D|², which avoids a sqrt and a full
27    // complex division per bin. |D|² < 1e-30 is the same cutoff as |D| < 1e-15.
28    let group_delay = Zip::from(weighted_freq.freq_data())
29        .and(signal.freq_data())
30        .map_collect(|numerator, denominator| {
31            let denominator_power = denominator.norm_sqr();
32            if denominator_power < 1e-30 {
33                0.0
34            } else {
35                (numerator * denominator.conj()).re / denominator_power
36            }
37        });
38
39    Ok(RealData::from_axis(signal.freq_axis().clone(), group_delay)
40        .expect("freq bins originate from FreqSignal"))
41}
42
43#[cfg(test)]
44mod tests {
45    use crate::signal::TimeSignal;
46    use approx::assert_abs_diff_eq;
47    use ndarray::arr2;
48
49    use super::*;
50
51    #[test]
52    fn group_delay_matches_direct_complex_division() {
53        let time = TimeSignal::new(
54            arr2(&[
55                [0.1, 0.7, -0.4, 0.2, 0.05, -0.1, 0.0, 0.02, 0.3],
56                [0.0, 0.0, 1.0, 0.5, 0.25, 0.0, 0.0, 0.0, 0.0],
57            ]),
58            9.0,
59        )
60        .unwrap();
61        let freq = time.to_freq();
62
63        let mut weighted = time.time_data().to_owned();
64        for mut channel in weighted.outer_iter_mut() {
65            for (n, sample) in channel.iter_mut().enumerate() {
66                *sample *= n as f64;
67            }
68        }
69        let weighted_freq = TimeSignal::new(weighted, 9.0).unwrap().into_freq();
70
71        let delay = group_delay(&freq).unwrap();
72        for ch in 0..2 {
73            for (bin, actual) in delay.channel(ch).iter().enumerate() {
74                let expected = (weighted_freq.channel(ch)[bin] / freq.channel(ch)[bin]).re;
75                assert_abs_diff_eq!(*actual, expected, epsilon = 1e-10);
76            }
77        }
78    }
79
80    #[test]
81    fn group_delay_matches_integer_sample_delay() {
82        let delayed = TimeSignal::new(arr2(&[[0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]]), 8.0)
83            .unwrap()
84            .into_freq();
85
86        let delay = group_delay(&delayed).unwrap();
87
88        for value in delay
89            .channel(0)
90            .iter()
91            .skip(1)
92            .take(delay.num_data_points() - 2)
93        {
94            assert_abs_diff_eq!(*value, 2.0, epsilon = 1e-10);
95        }
96    }
97}
audio_lab/dsp/frequency.rs

audio_lab/dsp/
frequency.rs
