comnum_math/stats.rs
1//! Second- and fourth-order statistics of complex sequences.
2//!
3//! - [`Moments`]: the moments $\mathbb{E}[r^{2}]$, $\mathbb{E}|r|^{2}$, $\mathbb{E}[r^{4}]$,
4//! $\mathbb{E}[r^{2}|r|^{2}]$ and $\mathbb{E}|r|^{4}$, of a sample (the means) or of a
5//! discrete distribution (the points and their probabilities): the cumulants of the
6//! blind estimators, the M2M4 estimator and the radius of the constant modulus
7//! algorithm all start from them;
8//! - [`Correlation`]: the cross-correlation $\sum_{k} a_{k} b_{k}^{*}$ of two sequences
9//! and their energies, with the normalised correlation
10//!
11//! $$
12//! \rho = \frac{\left|\sum_{k} a_{k} b_{k}^{*}\right|}{\sqrt{\sum_{k} |a_{k}|^{2}
13//! \sum_{k} |b_{k}|^{2}}} \in [0, 1],
14//! $$
15//!
16//! the quality of a frame synchronisation, the effective SNR of a sequence against
17//! another and the test of the singularity of a 2 × 2 equaliser.
18//!
19//! # Examples
20//! ```
21//! use comnum_math::C64;
22//! use comnum_math::stats::{Correlation, Moments};
23//! let a = [C64::new(1.0, 0.0), C64::new(0.0, 1.0)];
24//! let m = Moments::sample(&a);
25//! assert_eq!((m.m21, m.m42), (1.0, 1.0));
26//! assert_eq!(m.m20, C64::new(0.0, 0.0)); // 1 and -1
27//! let b: Vec<C64> = a.iter().map(|v| v * C64::new(0.0, 2.0)).collect();
28//! assert!((Correlation::between(&a, &b).rho() - 1.0).abs() < 1e-15);
29//! ```
30
31use crate::C64;
32
33/// The moments of order two and four of a complex random variable $r$ (see the module).
34#[derive(Clone, Copy, Debug, PartialEq)]
35pub struct Moments {
36 /// $\mathbb{E}[r^{2}]$.
37 pub m20: C64,
38 /// $\mathbb{E}|r|^{2}$.
39 pub m21: f64,
40 /// $\mathbb{E}[r^{4}]$.
41 pub m40: C64,
42 /// $\mathbb{E}[r^{2}|r|^{2}]$.
43 pub m41: C64,
44 /// $\mathbb{E}|r|^{4}$.
45 pub m42: f64,
46}
47
48impl Moments {
49 /// The sample moments of `x`: the sums divided by the number of samples (`NaN`
50 /// without a sample).
51 pub fn sample(x: &[C64]) -> Moments {
52 let n = x.len() as f64;
53 let s = Moments::sums(x.iter().map(|&v| (v, 1.0)));
54 Moments { m20: s.m20 / n, m21: s.m21 / n, m40: s.m40 / n, m41: s.m41 / n, m42: s.m42 / n }
55 }
56
57 /// The moments of the distribution of the points `x` with the probabilities `p` (the
58 /// pairs of the shorter of the two).
59 pub fn of_distribution(x: &[C64], p: &[f64]) -> Moments {
60 Moments::sums(x.iter().copied().zip(p.iter().copied()))
61 }
62
63 /// The weighted sums $\sum_{k} w_{k} r_{k}^{2}$, … (a weight of 1 is exact).
64 fn sums(rw: impl Iterator<Item = (C64, f64)>) -> Moments {
65 let zero = C64::new(0.0, 0.0);
66 let mut m = Moments { m20: zero, m21: 0.0, m40: zero, m41: zero, m42: 0.0 };
67 for (x, w) in rw {
68 let a = x.norm_sqr();
69 let x2 = x * x;
70 m.m20 += x2 * w;
71 m.m21 += a * w;
72 m.m40 += x2 * x2 * w;
73 m.m41 += x2 * a * w;
74 m.m42 += a * a * w;
75 }
76 m
77 }
78}
79
80/// The cross-correlation of two sequences and their energies (see the module).
81#[derive(Clone, Copy, Debug, PartialEq)]
82pub struct Correlation {
83 /// $\sum_{k} a_{k} b_{k}^{*}$.
84 pub cross: C64,
85 /// $\sum_{k} |a_{k}|^{2}$.
86 pub energy_a: f64,
87 /// $\sum_{k} |b_{k}|^{2}$.
88 pub energy_b: f64,
89}
90
91impl Correlation {
92 /// The correlation of `a` with `b` (the pairs of the shorter of the two).
93 pub fn between(a: &[C64], b: &[C64]) -> Correlation {
94 let mut c = Correlation { cross: C64::new(0.0, 0.0), energy_a: 0.0, energy_b: 0.0 };
95 for (u, v) in a.iter().zip(b) {
96 c.cross += u * v.conj();
97 c.energy_a += u.norm_sqr();
98 c.energy_b += v.norm_sqr();
99 }
100 c
101 }
102
103 /// The normalised correlation $\rho = |c| / \sqrt{E_{a} E_{b}}$, 0 when a sequence
104 /// has no energy.
105 pub fn rho(&self) -> f64 {
106 if self.energy_a > 0.0 && self.energy_b > 0.0 {
107 self.cross.norm() / (self.energy_a * self.energy_b).sqrt()
108 } else {
109 0.0
110 }
111 }
112
113 /// Its square $\rho^{2} = |c|^{2} / (E_{a} E_{b})$ (computed as such, not as the square
114 /// of [`Correlation::rho`]), 0 when a sequence has no energy.
115 pub fn rho_squared(&self) -> f64 {
116 if self.energy_a > 0.0 && self.energy_b > 0.0 {
117 self.cross.norm_sqr() / (self.energy_a * self.energy_b)
118 } else {
119 0.0
120 }
121 }
122}
123
124#[cfg(test)]
125mod tests {
126 use super::*;
127
128 #[test]
129 fn the_moments_of_a_sample_are_those_of_its_uniform_distribution() {
130 let x: Vec<C64> = (0..7).map(|k| C64::new(f64::from(k) - 3.0, 0.5 * f64::from(k))).collect();
131 let s = Moments::sample(&x);
132 let d = Moments::of_distribution(&x, &[1.0 / 7.0; 7]);
133 assert!((s.m21 - d.m21).abs() < 1e-14 && (s.m42 - d.m42).abs() < 1e-12);
134 assert!((s.m40 - d.m40).norm() < 1e-12 && (s.m41 - d.m41).norm() < 1e-12 && (s.m20 - d.m20).norm() < 1e-14);
135 let e = |k: i32| x.iter().map(|v| v.norm_sqr().powi(k)).sum::<f64>() / 7.0;
136 assert_eq!((s.m21, s.m42), (e(1), e(2)));
137 }
138
139 #[test]
140 fn the_correlation_is_one_for_proportional_sequences_and_zero_without_energy() {
141 let a = [C64::new(1.0, 2.0), C64::new(-0.5, 0.1), C64::new(0.0, -1.0)];
142 let b: Vec<C64> = a.iter().map(|v| v * C64::new(-1.5, 0.7)).collect();
143 let c = Correlation::between(&a, &b);
144 assert!((c.rho() - 1.0).abs() < 1e-15 && (c.rho_squared() - 1.0).abs() < 1e-15);
145 let z = Correlation::between(&a, &[C64::new(0.0, 0.0); 3]);
146 assert_eq!((z.rho(), z.rho_squared()), (0.0, 0.0));
147 }
148}