embedded_dsp/
statistics.rs

1//! Statistics functions (mean, variance, standard deviation, RMS, power, min, max, absmax, absmin, entropy, Kullback-Leibler, LogSumExp).
2
3#[allow(unused_imports)]
4use crate::math::FloatMath;
5use crate::types::*;
6
7// --- Mean ---
8
9pub fn mean_f32(src: &[f32], result: &mut f32) -> Status {
10    if src.is_empty() {
11        return Status::LengthError;
12    }
13    let mut sum = 0.0f32;
14    for &val in src {
15        sum += val;
16    }
17    *result = sum / (src.len() as f32);
18    Status::Success
19}
20
21pub fn mean_f64(src: &[f64], result: &mut f64) -> Status {
22    if src.is_empty() {
23        return Status::LengthError;
24    }
25    let mut sum = 0.0f64;
26    for &val in src {
27        sum += val;
28    }
29    *result = sum / (src.len() as f64);
30    Status::Success
31}
32
33pub fn mean_q31(src: &[q31], result: &mut q31) -> Status {
34    if src.is_empty() {
35        return Status::LengthError;
36    }
37    let mut sum: q63 = 0;
38    for &val in src {
39        sum += val.to_bits() as q63;
40    }
41    *result = q31::from_bits((sum / (src.len() as q63)) as i32);
42    Status::Success
43}
44
45pub fn mean_q15(src: &[q15], result: &mut q15) -> Status {
46    if src.is_empty() {
47        return Status::LengthError;
48    }
49    let mut sum: i32 = 0;
50    for &val in src {
51        sum += val.to_bits() as i32;
52    }
53    *result = q15::from_bits((sum / (src.len() as i32)) as i16);
54    Status::Success
55}
56
57pub fn mean_q7(src: &[q7], result: &mut q7) -> Status {
58    if src.is_empty() {
59        return Status::LengthError;
60    }
61    let mut sum: i32 = 0;
62    for &val in src {
63        sum += val.to_bits() as i32;
64    }
65    *result = q7::from_bits((sum / (src.len() as i32)) as i8);
66    Status::Success
67}
68
69// --- Variance ---
70
71pub fn var_f32(src: &[f32], result: &mut f32) -> Status {
72    if src.len() <= 1 {
73        return Status::LengthError;
74    }
75    let mut mean = 0.0f32;
76    mean_f32(src, &mut mean);
77    let mut sum_sq = 0.0f32;
78    for &val in src {
79        let diff = val - mean;
80        sum_sq += diff * diff;
81    }
82    *result = sum_sq / ((src.len() - 1) as f32);
83    Status::Success
84}
85
86pub fn var_f64(src: &[f64], result: &mut f64) -> Status {
87    if src.len() <= 1 {
88        return Status::LengthError;
89    }
90    let mut mean = 0.0f64;
91    mean_f64(src, &mut mean);
92    let mut sum_sq = 0.0f64;
93    for &val in src {
94        let diff = val - mean;
95        sum_sq += diff * diff;
96    }
97    *result = sum_sq / ((src.len() - 1) as f64);
98    Status::Success
99}
100
101pub fn var_q31(src: &[q31], result: &mut q31) -> Status {
102    if src.len() <= 1 {
103        return Status::LengthError;
104    }
105    let mut m = q31::ZERO;
106    mean_q31(src, &mut m);
107    let mut sum_sq: u64 = 0;
108    for &val in src {
109        let diff = (val.to_bits() as i64) - (m.to_bits() as i64);
110        sum_sq += ((diff * diff) >> 31) as u64;
111    }
112    *result = q31::from_bits(
113        ((sum_sq / (src.len() - 1) as u64) as i64).clamp(0, i32::MAX as i64) as i32,
114    );
115    Status::Success
116}
117
118pub fn var_q15(src: &[q15], result: &mut q15) -> Status {
119    if src.len() <= 1 {
120        return Status::LengthError;
121    }
122    let mut m = q15::ZERO;
123    mean_q15(src, &mut m);
124    let mut sum_sq: u32 = 0;
125    for &val in src {
126        let diff = (val.to_bits() as i32) - (m.to_bits() as i32);
127        sum_sq += ((diff * diff) >> 15) as u32;
128    }
129    *result = q15::from_bits(
130        ((sum_sq / (src.len() - 1) as u32) as i32).clamp(0, i16::MAX as i32) as i16,
131    );
132    Status::Success
133}
134
135pub fn var_q7(src: &[q7], result: &mut q7) -> Status {
136    if src.len() <= 1 {
137        return Status::LengthError;
138    }
139    let mut m = q7::ZERO;
140    mean_q7(src, &mut m);
141    let mut sum_sq: u32 = 0;
142    for &val in src {
143        let diff = (val.to_bits() as i32) - (m.to_bits() as i32);
144        sum_sq += ((diff * diff) >> 7) as u32;
145    }
146    *result =
147        q7::from_bits(((sum_sq / (src.len() - 1) as u32) as i32).clamp(0, i8::MAX as i32) as i8);
148    Status::Success
149}
150
151// --- Standard Deviation ---
152
153pub fn std_f32(src: &[f32], result: &mut f32) -> Status {
154    let mut v = 0.0f32;
155    let status = var_f32(src, &mut v);
156    if status == Status::Success {
157        *result = v.sqrt();
158    }
159    status
160}
161
162pub fn std_f64(src: &[f64], result: &mut f64) -> Status {
163    let mut v = 0.0f64;
164    let status = var_f64(src, &mut v);
165    if status == Status::Success {
166        *result = v.sqrt();
167    }
168    status
169}
170
171pub fn std_q31(src: &[q31], result: &mut q31) -> Status {
172    let mut v = q31::ZERO;
173    let status = var_q31(src, &mut v);
174    if status == Status::Success {
175        let _ = crate::fast_math::sqrt_q31(v, result);
176    }
177    status
178}
179
180pub fn std_q15(src: &[q15], result: &mut q15) -> Status {
181    let mut v = q15::ZERO;
182    let status = var_q15(src, &mut v);
183    if status == Status::Success {
184        let _ = crate::fast_math::sqrt_q15(v, result);
185    }
186    status
187}
188
189pub fn std_q7(src: &[q7], result: &mut q7) -> Status {
190    let mut v = q7::ZERO;
191    let status = var_q7(src, &mut v);
192    if status == Status::Success {
193        let n = (v.to_bits().max(0) as u32) << 7;
194        *result = q7::from_bits(crate::math::isqrt_u32(n).min(i8::MAX as u32) as i8);
195    }
196    status
197}
198
199// --- RMS ---
200
201pub fn rms_f32(src: &[f32], result: &mut f32) -> Status {
202    if src.is_empty() {
203        return Status::LengthError;
204    }
205    let mut sum_sq = 0.0f32;
206    for &val in src {
207        sum_sq += val * val;
208    }
209    *result = (sum_sq / (src.len() as f32)).sqrt();
210    Status::Success
211}
212
213pub fn rms_q31(src: &[q31], result: &mut q31) -> Status {
214    if src.is_empty() {
215        return Status::LengthError;
216    }
217    let mut sum_sq: u64 = 0;
218    for &val in src {
219        let v = val.to_bits() as i64;
220        sum_sq += ((v * v) >> 31) as u64;
221    }
222    let mean_sq = q31::from_bits((sum_sq / (src.len() as u64)).min(i32::MAX as u64) as i32);
223    let _ = crate::fast_math::sqrt_q31(mean_sq, result);
224    Status::Success
225}
226
227pub fn rms_q15(src: &[q15], result: &mut q15) -> Status {
228    if src.is_empty() {
229        return Status::LengthError;
230    }
231    let mut sum_sq: u32 = 0;
232    for &val in src {
233        let v = val.to_bits() as i32;
234        sum_sq += ((v * v) >> 15) as u32;
235    }
236    let mean_sq = q15::from_bits((sum_sq / (src.len() as u32)).min(i16::MAX as u32) as i16);
237    let _ = crate::fast_math::sqrt_q15(mean_sq, result);
238    Status::Success
239}
240
241// --- Power ---
242
243pub fn power_f32(src: &[f32], result: &mut f32) -> Status {
244    if src.is_empty() {
245        return Status::LengthError;
246    }
247    let mut sum_sq = 0.0f32;
248    for &val in src {
249        sum_sq += val * val;
250    }
251    *result = sum_sq;
252    Status::Success
253}
254
255pub fn power_q31(src: &[q31], result: &mut q63) -> Status {
256    if src.is_empty() {
257        return Status::LengthError;
258    }
259    let mut sum_sq: q63 = 0;
260    for &val in src {
261        let v = val.to_bits() as i64;
262        sum_sq += (v * v) >> 14;
263    }
264    *result = sum_sq;
265    Status::Success
266}
267
268pub fn power_q15(src: &[q15], result: &mut q63) -> Status {
269    if src.is_empty() {
270        return Status::LengthError;
271    }
272    let mut sum_sq: q63 = 0;
273    for &val in src {
274        let v = val.to_bits() as i32;
275        sum_sq += (v * v) as q63;
276    }
277    *result = sum_sq;
278    Status::Success
279}
280
281pub fn power_q7(src: &[q7], result: &mut q31) -> Status {
282    if src.is_empty() {
283        return Status::LengthError;
284    }
285    let mut sum_sq = q31::ZERO;
286    for &val in src {
287        let v = val.to_bits() as i32;
288        sum_sq += q31::from_bits(v * v);
289    }
290    *result = sum_sq;
291    Status::Success
292}
293
294// --- Min & Max ---
295
296pub fn min_f32(src: &[f32], result: &mut f32, index: &mut usize) -> Status {
297    if src.is_empty() {
298        return Status::LengthError;
299    }
300    let mut min_val = src[0];
301    let mut min_idx = 0;
302    for (i, &val) in src.iter().enumerate().skip(1) {
303        if val < min_val {
304            min_val = val;
305            min_idx = i;
306        }
307    }
308    *result = min_val;
309    *index = min_idx;
310    Status::Success
311}
312
313pub fn max_f32(src: &[f32], result: &mut f32, index: &mut usize) -> Status {
314    if src.is_empty() {
315        return Status::LengthError;
316    }
317    let mut max_val = src[0];
318    let mut max_idx = 0;
319    for (i, &val) in src.iter().enumerate().skip(1) {
320        if val > max_val {
321            max_val = val;
322            max_idx = i;
323        }
324    }
325    *result = max_val;
326    *index = max_idx;
327    Status::Success
328}
329
330pub fn min_q31(src: &[q31], result: &mut q31, index: &mut usize) -> Status {
331    if src.is_empty() {
332        return Status::LengthError;
333    }
334    let mut min_val = src[0];
335    let mut min_idx = 0;
336    for (i, &val) in src.iter().enumerate().skip(1) {
337        if val < min_val {
338            min_val = val;
339            min_idx = i;
340        }
341    }
342    *result = min_val;
343    *index = min_idx;
344    Status::Success
345}
346
347pub fn max_q31(src: &[q31], result: &mut q31, index: &mut usize) -> Status {
348    if src.is_empty() {
349        return Status::LengthError;
350    }
351    let mut max_val = src[0];
352    let mut max_idx = 0;
353    for (i, &val) in src.iter().enumerate().skip(1) {
354        if val > max_val {
355            max_val = val;
356            max_idx = i;
357        }
358    }
359    *result = max_val;
360    *index = max_idx;
361    Status::Success
362}
363
364pub fn min_q15(src: &[q15], result: &mut q15, index: &mut usize) -> Status {
365    if src.is_empty() {
366        return Status::LengthError;
367    }
368    let mut min_val = src[0];
369    let mut min_idx = 0;
370    for (i, &val) in src.iter().enumerate().skip(1) {
371        if val < min_val {
372            min_val = val;
373            min_idx = i;
374        }
375    }
376    *result = min_val;
377    *index = min_idx;
378    Status::Success
379}
380
381pub fn max_q15(src: &[q15], result: &mut q15, index: &mut usize) -> Status {
382    if src.is_empty() {
383        return Status::LengthError;
384    }
385    let mut max_val = src[0];
386    let mut max_idx = 0;
387    for (i, &val) in src.iter().enumerate().skip(1) {
388        if val > max_val {
389            max_val = val;
390            max_idx = i;
391        }
392    }
393    *result = max_val;
394    *index = max_idx;
395    Status::Success
396}
397
398pub fn min_q7(src: &[q7], result: &mut q7, index: &mut usize) -> Status {
399    if src.is_empty() {
400        return Status::LengthError;
401    }
402    let mut min_val = src[0];
403    let mut min_idx = 0;
404    for (i, &val) in src.iter().enumerate().skip(1) {
405        if val < min_val {
406            min_val = val;
407            min_idx = i;
408        }
409    }
410    *result = min_val;
411    *index = min_idx;
412    Status::Success
413}
414
415pub fn max_q7(src: &[q7], result: &mut q7, index: &mut usize) -> Status {
416    if src.is_empty() {
417        return Status::LengthError;
418    }
419    let mut max_val = src[0];
420    let mut max_idx = 0;
421    for (i, &val) in src.iter().enumerate().skip(1) {
422        if val > max_val {
423            max_val = val;
424            max_idx = i;
425        }
426    }
427    *result = max_val;
428    *index = max_idx;
429    Status::Success
430}
431
432// --- Absmax & Absmin ---
433
434pub fn absmax_f32(src: &[f32], result: &mut f32, index: &mut usize) -> Status {
435    if src.is_empty() {
436        return Status::LengthError;
437    }
438    let mut max_val = src[0].abs();
439    let mut max_idx = 0;
440    for (i, &val) in src.iter().enumerate().skip(1) {
441        let abs_val = val.abs();
442        if abs_val > max_val {
443            max_val = abs_val;
444            max_idx = i;
445        }
446    }
447    *result = max_val;
448    *index = max_idx;
449    Status::Success
450}
451
452pub fn absmin_f32(src: &[f32], result: &mut f32, index: &mut usize) -> Status {
453    if src.is_empty() {
454        return Status::LengthError;
455    }
456    let mut min_val = src[0].abs();
457    let mut min_idx = 0;
458    for (i, &val) in src.iter().enumerate().skip(1) {
459        let abs_val = val.abs();
460        if abs_val < min_val {
461            min_val = abs_val;
462            min_idx = i;
463        }
464    }
465    *result = min_val;
466    *index = min_idx;
467    Status::Success
468}
469
470// --- Entropy, KL Divergence, LogSumExp ---
471
472pub fn entropy_f32(src: &[f32]) -> f32 {
473    let mut ent = 0.0f32;
474    for &p in src {
475        if p > 0.0 {
476            ent -= p * p.ln();
477        }
478    }
479    ent
480}
481
482pub fn kullback_leibler_f32(p: &[f32], q: &[f32]) -> f32 {
483    let len = p.len().min(q.len());
484    let mut kl = 0.0f32;
485    for i in 0..len {
486        if p[i] > 0.0 && q[i] > 0.0 {
487            kl += p[i] * (p[i] / q[i]).ln();
488        }
489    }
490    kl
491}
492
493pub fn logsumexp_f32(src: &[f32]) -> f32 {
494    if src.is_empty() {
495        return 0.0;
496    }
497    let mut max_v = src[0];
498    for &v in src.iter().skip(1) {
499        if v > max_v {
500            max_v = v;
501        }
502    }
503    let mut sum_exp = 0.0f32;
504    for &v in src {
505        sum_exp += (v - max_v).exp();
506    }
507    max_v + sum_exp.ln()
508}
embedded_dsp/statistics.rs

embedded_dsp/
statistics.rs
