sonora-simd 0.1.0

SIMD abstraction layer for WebRTC Audio Processing
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

//! Scalar fallback implementations of SIMD operations.

pub(crate) fn dot_product(a: &[f32], b: &[f32]) -> f32 {
    a.iter().zip(b.iter()).map(|(x, y)| x * y).sum()
}

pub(crate) fn dual_dot_product(input: &[f32], k1: &[f32], k2: &[f32]) -> (f32, f32) {
    let mut sum1 = 0.0f32;
    let mut sum2 = 0.0f32;
    for i in 0..input.len() {
        sum1 += input[i] * k1[i];
        sum2 += input[i] * k2[i];
    }
    (sum1, sum2)
}

/// Sinc resampler convolution: dual dot product with interpolation in f64.
///
/// Matches C++ `Convolve_C` which interpolates in `double` precision.
pub(crate) fn convolve_sinc(
    input: &[f32],
    k1: &[f32],
    k2: &[f32],
    kernel_interpolation_factor: f64,
) -> f32 {
    let mut sum1 = 0.0f32;
    let mut sum2 = 0.0f32;
    for i in 0..input.len() {
        sum1 += input[i] * k1[i];
        sum2 += input[i] * k2[i];
    }
    // Interpolate in f64 matching C++ Convolve_C:
    //   static_cast<float>((1.0 - kernel_interpolation_factor) * sum1 +
    //                       kernel_interpolation_factor * sum2);
    ((1.0 - kernel_interpolation_factor) * sum1 as f64 + kernel_interpolation_factor * sum2 as f64)
        as f32
}

pub(crate) fn multiply_accumulate(acc: &mut [f32], a: &[f32], b: &[f32]) {
    for i in 0..acc.len() {
        acc[i] += a[i] * b[i];
    }
}

pub(crate) fn sum(x: &[f32]) -> f32 {
    x.iter().sum()
}

pub(crate) fn elementwise_sqrt(x: &mut [f32]) {
    for v in x.iter_mut() {
        *v = v.sqrt();
    }
}

pub(crate) fn elementwise_multiply(x: &[f32], y: &[f32], z: &mut [f32]) {
    for i in 0..z.len() {
        z[i] = x[i] * y[i];
    }
}

pub(crate) fn elementwise_accumulate(x: &[f32], z: &mut [f32]) {
    for i in 0..z.len() {
        z[i] += x[i];
    }
}

pub(crate) fn power_spectrum(re: &[f32], im: &[f32], out: &mut [f32]) {
    for i in 0..out.len() {
        out[i] = re[i] * re[i] + im[i] * im[i];
    }
}

pub(crate) fn elementwise_min(a: &[f32], b: &[f32], out: &mut [f32]) {
    for i in 0..out.len() {
        out[i] = a[i].min(b[i]);
    }
}

pub(crate) fn elementwise_max(a: &[f32], b: &[f32], out: &mut [f32]) {
    for i in 0..out.len() {
        out[i] = a[i].max(b[i]);
    }
}

pub(crate) fn complex_multiply_accumulate(
    x_re: &[f32],
    x_im: &[f32],
    h_re: &[f32],
    h_im: &[f32],
    acc_re: &mut [f32],
    acc_im: &mut [f32],
) {
    for i in 0..acc_re.len() {
        // (x_re + j*x_im) * (h_re + j*h_im) =
        //   (x_re*h_re - x_im*h_im) + j*(x_re*h_im + x_im*h_re)
        // Note: AEC3 uses conjugate multiply convention (+ for real, - for imag):
        //   real = x_re*h_re + x_im*h_im
        //   imag = x_re*h_im - x_im*h_re
        acc_re[i] += x_re[i] * h_re[i] + x_im[i] * h_im[i];
        acc_im[i] += x_re[i] * h_im[i] - x_im[i] * h_re[i];
    }
}

pub(crate) fn complex_multiply_accumulate_standard(
    x_re: &[f32],
    x_im: &[f32],
    h_re: &[f32],
    h_im: &[f32],
    acc_re: &mut [f32],
    acc_im: &mut [f32],
) {
    for i in 0..acc_re.len() {
        // Standard complex multiply: (x_re + j*x_im) * (h_re + j*h_im)
        //   real = x_re*h_re - x_im*h_im
        //   imag = x_re*h_im + x_im*h_re
        acc_re[i] += x_re[i] * h_re[i] - x_im[i] * h_im[i];
        acc_im[i] += x_re[i] * h_im[i] + x_im[i] * h_re[i];
    }
}

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

    #[test]
    fn test_dot_product_empty() {
        assert_eq!(dot_product(&[], &[]), 0.0);
    }

    #[test]
    fn test_dot_product_single() {
        assert_eq!(dot_product(&[3.0], &[4.0]), 12.0);
    }

    #[test]
    fn test_sum_empty() {
        assert_eq!(sum(&[]), 0.0);
    }
}