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

embedded_dsp/
statistics.rs
