oximedia-audio-analysis 0.1.7

Audio analysis tools for media including loudness, spectrum, and speech detection
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

//! Total harmonic distortion (THD) measurement.

use oxifft::Complex;

/// Compute total harmonic distortion (THD).
///
/// THD is the ratio of the sum of harmonic powers to the fundamental power.
///
/// # Arguments
/// * `samples` - Audio samples
/// * `sample_rate` - Sample rate in Hz
///
/// # Returns
/// THD value (0-1, lower is better)
#[must_use]
pub fn total_harmonic_distortion(samples: &[f32], _sample_rate: f32) -> f32 {
    if samples.is_empty() {
        return 0.0;
    }

    let fft_size = samples.len().next_power_of_two();

    // Prepare FFT input
    let buffer: Vec<Complex<f64>> = samples
        .iter()
        .map(|&s| Complex::new(f64::from(s), 0.0))
        .chain(std::iter::repeat(Complex::new(0.0, 0.0)))
        .take(fft_size)
        .collect();

    let buffer = oxifft::fft(&buffer);

    // Compute magnitude spectrum
    let magnitude: Vec<f32> = buffer[..fft_size / 2]
        .iter()
        .map(|c| c.norm() as f32)
        .collect();

    // Find fundamental frequency peak
    let fundamental_bin = magnitude
        .iter()
        .enumerate()
        .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
        .map_or(0, |(i, _)| i);

    if fundamental_bin == 0 {
        return 0.0;
    }

    let fundamental_power = magnitude[fundamental_bin].powi(2);

    // Sum harmonic powers (2f, 3f, 4f, 5f)
    let mut harmonic_power = 0.0;
    for harmonic in 2..=5 {
        let harmonic_bin = fundamental_bin * harmonic;
        if harmonic_bin < magnitude.len() {
            harmonic_power += magnitude[harmonic_bin].powi(2);
        }
    }

    if fundamental_power > 0.0 {
        (harmonic_power / fundamental_power).sqrt()
    } else {
        0.0
    }
}

/// Compute THD+N (THD plus noise).
#[must_use]
pub fn thd_plus_noise(samples: &[f32], _sample_rate: f32) -> f32 {
    if samples.is_empty() {
        return 0.0;
    }

    let fft_size = samples.len().next_power_of_two();

    let buffer: Vec<Complex<f64>> = samples
        .iter()
        .map(|&s| Complex::new(f64::from(s), 0.0))
        .chain(std::iter::repeat(Complex::new(0.0, 0.0)))
        .take(fft_size)
        .collect();

    let buffer = oxifft::fft(&buffer);

    let magnitude: Vec<f32> = buffer[..fft_size / 2]
        .iter()
        .map(|c| c.norm() as f32)
        .collect();

    let fundamental_bin = magnitude
        .iter()
        .enumerate()
        .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
        .map_or(0, |(i, _)| i);

    if fundamental_bin == 0 {
        return 0.0;
    }

    let fundamental_power = magnitude[fundamental_bin].powi(2);
    let total_power: f32 = magnitude.iter().map(|&m| m.powi(2)).sum();

    let noise_plus_harmonics = total_power - fundamental_power;

    if fundamental_power > 0.0 {
        (noise_plus_harmonics / fundamental_power).sqrt()
    } else {
        0.0
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_thd_clean_signal() {
        // Pure sine wave should have very low THD
        let sample_rate = 44100.0;
        let samples: Vec<f32> = (0..4096)
            .map(|i| (2.0 * std::f32::consts::PI * 440.0 * i as f32 / sample_rate).sin())
            .collect();

        let thd = total_harmonic_distortion(&samples, sample_rate);
        assert!(thd < 0.1);
    }

    #[test]
    fn test_thd_distorted_signal() {
        // Clipped sine wave should have higher THD
        let sample_rate = 44100.0;
        let samples: Vec<f32> = (0..4096)
            .map(|i| {
                let x = (2.0 * std::f32::consts::PI * 440.0 * i as f32 / sample_rate).sin();
                x.clamp(-0.5, 0.5) // Clipping
            })
            .collect();

        let thd = total_harmonic_distortion(&samples, sample_rate);
        assert!(thd > 0.05);
    }
}