math-dsp 0.5.18

DSP utilities: signal generation, FFT analysis, and audio analysis tools
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

use math_audio_dsp::analysis::{compute_rt60_broadband, find_db_point};

#[test]
fn test_find_db_point_flat() {
    let frequencies = vec![20.0, 100.0, 1000.0, 20000.0];
    let magnitude_db = vec![0.0, 0.0, 0.0, 0.0];

    // In a flat response, -3dB is never reached
    assert_eq!(find_db_point(&frequencies, &magnitude_db, -3.0, true), None);
    assert_eq!(
        find_db_point(&frequencies, &magnitude_db, -3.0, false),
        None
    );
}

#[test]
fn test_find_db_point_sloping_down() {
    let frequencies = vec![100.0, 200.0, 400.0, 800.0];
    let magnitude_db = vec![0.0, -2.0, -4.0, -6.0];

    // -3dB should be between 200Hz and 400Hz
    // Linear interpolation: -2.0 + t*(-4.0 - (-2.0)) = -3.0 => -2.0 - 2t = -3.0 => 2t = 1.0 => t = 0.5
    // Freq: 200 + 0.5 * (400 - 200) = 300
    let p3db = find_db_point(&frequencies, &magnitude_db, -3.0, true);
    assert!(p3db.is_some());
    assert!((p3db.unwrap() - 300.0).abs() < 1e-3);
}

#[test]
fn test_find_db_point_sloping_up() {
    let frequencies = vec![100.0, 200.0, 400.0, 800.0];
    let magnitude_db = vec![-6.0, -4.0, -2.0, 0.0];

    // -3dB should be between 200Hz and 400Hz (ascending)
    // Linear interpolation: -4.0 + t*(-2.0 - (-4.0)) = -3.0 => -4.0 + 2t = -3.0 => 2t = 1.0 => t = 0.5
    // Freq: 200 + 0.5 * (400 - 200) = 300
    let p3db = find_db_point(&frequencies, &magnitude_db, -3.0, true);
    assert!(p3db.is_some());
    assert!((p3db.unwrap() - 300.0).abs() < 1e-3);
}

#[test]
fn test_find_db_point_subwoofer() {
    let frequencies = vec![20.0, 40.0, 80.0, 120.0, 160.0, 200.0];
    let magnitude_db = vec![-10.0, -3.0, 0.0, 0.0, -3.0, -10.0];

    // Low -3dB point at 40Hz
    let low_3db = find_db_point(&frequencies, &magnitude_db, -3.0, true);
    assert!(low_3db.is_some());
    assert!((low_3db.unwrap() - 40.0).abs() < 1e-3);

    // High -3dB point at 160Hz (searching from the end)
    let high_3db = find_db_point(&frequencies, &magnitude_db, -3.0, false);
    assert!(high_3db.is_some());
    assert!((high_3db.unwrap() - 160.0).abs() < 1e-3);
}

#[test]
fn test_compute_average_response_full_band() {
    use math_audio_dsp::analysis::compute_average_response;
    let frequencies = vec![100.0, 200.0, 400.0, 800.0];
    let magnitude_db = vec![10.0, 10.0, 10.0, 10.0];

    // Flat response should average to 10.0
    let avg = compute_average_response(&frequencies, &magnitude_db, None);
    assert!((avg - 10.0).abs() < 1e-3);
}

#[test]
fn test_compute_average_response_range() {
    use math_audio_dsp::analysis::compute_average_response;
    let frequencies = vec![100.0, 200.0, 400.0, 800.0];
    let magnitude_db = vec![0.0, 10.0, 10.0, 20.0];

    // Average between 200Hz and 400Hz (where both are 10dB)
    let avg = compute_average_response(&frequencies, &magnitude_db, Some((200.0, 400.0)));
    assert!((avg - 10.0).abs() < 1e-3);
}

#[test]
fn test_compute_average_response_log_spacing() {
    use math_audio_dsp::analysis::compute_average_response;
    let frequencies = vec![10.0, 100.0, 1000.0];
    let magnitude_db = vec![0.0, 10.0, 20.0];

    // Log-weighted average should be more accurate for acoustic data
    // Here we'll just test that it returns a value within range for now
    let avg = compute_average_response(&frequencies, &magnitude_db, None);
    assert!(avg > 0.0 && avg < 20.0);
}

fn make_exponential_decay(num_samples: usize, sample_rate: f32, rt60_seconds: f32) -> Vec<f32> {
    let k = 3.0 * std::f32::consts::LN_10 / rt60_seconds;
    (0..num_samples)
        .map(|i| {
            let t = i as f32 / sample_rate;
            (-k * t).exp()
        })
        .collect()
}

fn lcg_noise(n: usize, seed: u32, amplitude: f32) -> Vec<f32> {
    let mut state = seed;
    (0..n)
        .map(|_| {
            state = state.wrapping_mul(1664525).wrapping_add(1013904223);
            let normalised = (state as f32) / (u32::MAX as f32) - 0.5;
            normalised * 2.0 * amplitude
        })
        .collect()
}

#[test]
fn test_compute_rt60_broadband_recovers_exponential_decay() {
    let sample_rate = 48_000.0_f32;
    let expected_rt60 = 0.5_f32;
    let impulse = make_exponential_decay((sample_rate * 2.0) as usize, sample_rate, expected_rt60);

    let rt60 = compute_rt60_broadband(&impulse, sample_rate);

    assert!(
        (rt60 - expected_rt60).abs() < 0.03,
        "expected RT60 near {expected_rt60:.3}s, got {rt60:.3}s"
    );
}

#[test]
fn test_compute_rt60_broadband_trims_noise_tail_before_fit() {
    let sample_rate = 48_000.0_f32;
    let expected_rt60 = 0.45_f32;
    let mut impulse =
        make_exponential_decay((sample_rate * 0.8) as usize, sample_rate, expected_rt60);
    impulse.extend(lcg_noise((sample_rate * 1.2) as usize, 0x1234_5678, 0.0005));

    let rt60 = compute_rt60_broadband(&impulse, sample_rate);

    assert!(
        (rt60 - expected_rt60).abs() < 0.06,
        "expected tail-trimmed RT60 near {expected_rt60:.3}s, got {rt60:.3}s"
    );
}

#[test]
fn test_compute_rt60_broadband_rejects_insufficient_dynamic_range() {
    let sample_rate = 48_000.0_f32;
    let impulse = make_exponential_decay((sample_rate * 0.050) as usize, sample_rate, 0.5);

    assert_eq!(compute_rt60_broadband(&impulse, sample_rate), 0.0);
}