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 vf = v as f64 / 2147483648.0;
171        let s = vf.sqrt();
172        *result = (s * 2147483647.0).clamp(0.0, 2147483647.0) as q31;
173    }
174    status
175}
176
177pub fn std_q15(src: &[q15], result: &mut q15) -> Status {
178    let mut v = 0;
179    let status = var_q15(src, &mut v);
180    if status == Status::Success {
181        let vf = v as f32 / 32768.0;
182        let s = vf.sqrt();
183        *result = (s * 32767.0).clamp(0.0, 32767.0) as q15;
184    }
185    status
186}
187
188pub fn std_q7(src: &[q7], result: &mut q7) -> Status {
189    let mut v = 0;
190    let status = var_q7(src, &mut v);
191    if status == Status::Success {
192        let vf = v as f32 / 128.0;
193        let s = vf.sqrt();
194        *result = (s * 127.0).clamp(0.0, 127.0) as q7;
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 as i64;
220        sum_sq += ((v * v) >> 31) as u64;
221    }
222    let mean_sq = sum_sq / (src.len() as u64);
223    let vf = mean_sq as f64 / 2147483648.0;
224    let r = vf.sqrt();
225    *result = (r * 2147483647.0).clamp(0.0, 2147483647.0) as q31;
226    Status::Success
227}
228
229pub fn rms_q15(src: &[q15], result: &mut q15) -> Status {
230    if src.is_empty() {
231        return Status::LengthError;
232    }
233    let mut sum_sq: u32 = 0;
234    for &val in src {
235        let v = val as i32;
236        sum_sq += ((v * v) >> 15) as u32;
237    }
238    let mean_sq = sum_sq / (src.len() as u32);
239    let vf = mean_sq as f32 / 32768.0;
240    let r = vf.sqrt();
241    *result = (r * 32767.0).clamp(0.0, 32767.0) as q15;
242    Status::Success
243}
244
245// --- Power ---
246
247pub fn power_f32(src: &[f32], result: &mut f32) -> Status {
248    if src.is_empty() {
249        return Status::LengthError;
250    }
251    let mut sum_sq = 0.0f32;
252    for &val in src {
253        sum_sq += val * val;
254    }
255    *result = sum_sq;
256    Status::Success
257}
258
259pub fn power_q31(src: &[q31], result: &mut q63) -> Status {
260    if src.is_empty() {
261        return Status::LengthError;
262    }
263    let mut sum_sq: q63 = 0;
264    for &val in src {
265        let v = val as i64;
266        sum_sq += (v * v) >> 14;
267    }
268    *result = sum_sq;
269    Status::Success
270}
271
272pub fn power_q15(src: &[q15], result: &mut q63) -> Status {
273    if src.is_empty() {
274        return Status::LengthError;
275    }
276    let mut sum_sq: q63 = 0;
277    for &val in src {
278        let v = val as i32;
279        sum_sq += (v * v) as q63;
280    }
281    *result = sum_sq;
282    Status::Success
283}
284
285pub fn power_q7(src: &[q7], result: &mut q31) -> Status {
286    if src.is_empty() {
287        return Status::LengthError;
288    }
289    let mut sum_sq: q31 = 0;
290    for &val in src {
291        let v = val as i32;
292        sum_sq += v * v;
293    }
294    *result = sum_sq;
295    Status::Success
296}
297
298// --- Min & Max ---
299
300pub fn min_f32(src: &[f32], result: &mut f32, index: &mut usize) -> Status {
301    if src.is_empty() {
302        return Status::LengthError;
303    }
304    let mut min_val = src[0];
305    let mut min_idx = 0;
306    for (i, &val) in src.iter().enumerate().skip(1) {
307        if val < min_val {
308            min_val = val;
309            min_idx = i;
310        }
311    }
312    *result = min_val;
313    *index = min_idx;
314    Status::Success
315}
316
317pub fn max_f32(src: &[f32], result: &mut f32, index: &mut usize) -> Status {
318    if src.is_empty() {
319        return Status::LengthError;
320    }
321    let mut max_val = src[0];
322    let mut max_idx = 0;
323    for (i, &val) in src.iter().enumerate().skip(1) {
324        if val > max_val {
325            max_val = val;
326            max_idx = i;
327        }
328    }
329    *result = max_val;
330    *index = max_idx;
331    Status::Success
332}
333
334pub fn min_q31(src: &[q31], result: &mut q31, index: &mut usize) -> Status {
335    if src.is_empty() {
336        return Status::LengthError;
337    }
338    let mut min_val = src[0];
339    let mut min_idx = 0;
340    for (i, &val) in src.iter().enumerate().skip(1) {
341        if val < min_val {
342            min_val = val;
343            min_idx = i;
344        }
345    }
346    *result = min_val;
347    *index = min_idx;
348    Status::Success
349}
350
351pub fn max_q31(src: &[q31], result: &mut q31, index: &mut usize) -> Status {
352    if src.is_empty() {
353        return Status::LengthError;
354    }
355    let mut max_val = src[0];
356    let mut max_idx = 0;
357    for (i, &val) in src.iter().enumerate().skip(1) {
358        if val > max_val {
359            max_val = val;
360            max_idx = i;
361        }
362    }
363    *result = max_val;
364    *index = max_idx;
365    Status::Success
366}
367
368pub fn min_q15(src: &[q15], result: &mut q15, index: &mut usize) -> Status {
369    if src.is_empty() {
370        return Status::LengthError;
371    }
372    let mut min_val = src[0];
373    let mut min_idx = 0;
374    for (i, &val) in src.iter().enumerate().skip(1) {
375        if val < min_val {
376            min_val = val;
377            min_idx = i;
378        }
379    }
380    *result = min_val;
381    *index = min_idx;
382    Status::Success
383}
384
385pub fn max_q15(src: &[q15], result: &mut q15, index: &mut usize) -> Status {
386    if src.is_empty() {
387        return Status::LengthError;
388    }
389    let mut max_val = src[0];
390    let mut max_idx = 0;
391    for (i, &val) in src.iter().enumerate().skip(1) {
392        if val > max_val {
393            max_val = val;
394            max_idx = i;
395        }
396    }
397    *result = max_val;
398    *index = max_idx;
399    Status::Success
400}
401
402pub fn min_q7(src: &[q7], result: &mut q7, index: &mut usize) -> Status {
403    if src.is_empty() {
404        return Status::LengthError;
405    }
406    let mut min_val = src[0];
407    let mut min_idx = 0;
408    for (i, &val) in src.iter().enumerate().skip(1) {
409        if val < min_val {
410            min_val = val;
411            min_idx = i;
412        }
413    }
414    *result = min_val;
415    *index = min_idx;
416    Status::Success
417}
418
419pub fn max_q7(src: &[q7], result: &mut q7, index: &mut usize) -> Status {
420    if src.is_empty() {
421        return Status::LengthError;
422    }
423    let mut max_val = src[0];
424    let mut max_idx = 0;
425    for (i, &val) in src.iter().enumerate().skip(1) {
426        if val > max_val {
427            max_val = val;
428            max_idx = i;
429        }
430    }
431    *result = max_val;
432    *index = max_idx;
433    Status::Success
434}
435
436// --- Absmax & Absmin ---
437
438pub fn absmax_f32(src: &[f32], result: &mut f32, index: &mut usize) -> Status {
439    if src.is_empty() {
440        return Status::LengthError;
441    }
442    let mut max_val = src[0].abs();
443    let mut max_idx = 0;
444    for (i, &val) in src.iter().enumerate().skip(1) {
445        let abs_val = val.abs();
446        if abs_val > max_val {
447            max_val = abs_val;
448            max_idx = i;
449        }
450    }
451    *result = max_val;
452    *index = max_idx;
453    Status::Success
454}
455
456pub fn absmin_f32(src: &[f32], result: &mut f32, index: &mut usize) -> Status {
457    if src.is_empty() {
458        return Status::LengthError;
459    }
460    let mut min_val = src[0].abs();
461    let mut min_idx = 0;
462    for (i, &val) in src.iter().enumerate().skip(1) {
463        let abs_val = val.abs();
464        if abs_val < min_val {
465            min_val = abs_val;
466            min_idx = i;
467        }
468    }
469    *result = min_val;
470    *index = min_idx;
471    Status::Success
472}
473
474// --- Entropy, KL Divergence, LogSumExp ---
475
476pub fn entropy_f32(src: &[f32]) -> f32 {
477    let mut ent = 0.0f32;
478    for &p in src {
479        if p > 0.0 {
480            ent -= p * p.ln();
481        }
482    }
483    ent
484}
485
486pub fn kullback_leibler_f32(p: &[f32], q: &[f32]) -> f32 {
487    let len = p.len().min(q.len());
488    let mut kl = 0.0f32;
489    for i in 0..len {
490        if p[i] > 0.0 && q[i] > 0.0 {
491            kl += p[i] * (p[i] / q[i]).ln();
492        }
493    }
494    kl
495}
496
497pub fn logsumexp_f32(src: &[f32]) -> f32 {
498    if src.is_empty() {
499        return 0.0;
500    }
501    let mut max_v = src[0];
502    for &v in src.iter().skip(1) {
503        if v > max_v {
504            max_v = v;
505        }
506    }
507    let mut sum_exp = 0.0f32;
508    for &v in src {
509        sum_exp += (v - max_v).exp();
510    }
511    max_v + sum_exp.ln()
512}
embedded_dsp/statistics.rs

embedded_dsp/
statistics.rs
