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
//! Coherence test of the residual and the input, for any input.
use num_traits::Float;
use rustfft::FftNum;
use super::{Validation, ValidationError};
impl<T: Float + FftNum> Validation<T> {
/// Coherence test of the residual and the input, for any input (non-periodic: random, chirp,
/// measured operation data, ...).
///
/// The evaluated samples are split into segments of `segment_len` samples (Hann window, 50 %
/// overlap, mean removed per segment), and with the DFTs `E_l(f)`, `U_l(f)` of segment `l`
///
/// ```text
/// γ²(f) = |Σ_l E_l(f) U_l(f)*|^2 / (Σ_l |E_l(f)|^2 Σ_l |U_l(f)|^2)
/// ```
///
/// at the bins `f = 1 ..= segment_len / 2` (frequency `f / (segment_len ts)`). If the model is
/// right, the residual is noise independent of the input, and for `L` independent segments
/// `γ²(f)` follows the Beta distribution `Beta(1, L - 1)` whatever the input (it is the squared
/// projection of the Gaussian vector `E_l(f)` on the fixed direction `U_l(f)`): the bound at
/// the given `confidence` is `1 - (1 - confidence)^(1 / (L - 1))`, and the mean `1 / L`.
/// Overlapping segments are not independent; `L` is the equivalent number of independent
/// segments for the window and the overlap (Welch: `K / (1 + 2 Σ_m (1 - m/K) c_m^2)` for `K`
/// segments, `c_m` the normalized overlap of the window with itself shifted by `m` steps).
///
/// `γ²(f)` outside the bound means the residual still depends on the input at that frequency.
/// Longer segments give a finer frequency resolution but fewer segments, so a higher bound
/// (less power); at the frequencies where the input has no power the test cannot detect
/// anything (see `CoherenceTest::excited`).
pub fn coherence_test(&self, segment_len: usize, confidence: T) -> Result<CoherenceTest<T>, ValidationError> {
let spectra = self.spectra(segment_len)?;
let l = spectra.effective_segments;
let bins: Vec<usize> = (1..=segment_len / 2).collect();
Ok(CoherenceTest {
coherence: bins
.iter()
.map(|&f| {
let denom = spectra.ee[f] * spectra.uu[f];
if denom > T::zero() { spectra.eu[f].norm_sqr() / denom } else { T::zero() }
})
.collect(),
input_power: bins.iter().map(|&f| spectra.uu[f]).collect(),
bins,
segment_len,
segments: spectra.segments,
effective_segments: l,
bound: T::one() - (T::one() - confidence).powf(T::one() / (l - T::one())),
})
}
}
/// Result of `Validation::coherence_test`.
#[derive(Clone, Debug)]
pub struct CoherenceTest<T> {
/// Frequency bins `f` (frequency `f / (segment_len ts)`).
pub bins: Vec<usize>,
/// Coherence `γ²(f)` of the residual and the input at each bin.
pub coherence: Vec<T>,
/// Input power at each bin (mean `|U_l(f)|^2` over the segments, windowed), for `excited`.
pub input_power: Vec<T>,
/// Segment length \[samples\].
pub segment_len: usize,
/// Number of (overlapping) segments `K`.
pub segments: usize,
/// Equivalent number of independent segments `L`.
pub effective_segments: T,
/// Bound of `γ²(f)` at the given confidence.
pub bound: T,
}
impl<T: Float> CoherenceTest<T> {
/// Frequencies \[Hz\] of the bins for the sampling period `ts`.
pub fn frequencies(&self, ts: T) -> Vec<T> {
let df = T::one() / (T::from(self.segment_len).unwrap() * ts);
self.bins.iter().map(|&f| T::from(f).unwrap() * df).collect()
}
/// The bins where the input power is at least `relative` times its maximum (e.g. `1e-2` for
/// -20 dB): elsewhere the test has no power to detect a model error.
pub fn excited(&self, relative: T) -> Self {
let threshold = relative * self.input_power.iter().fold(T::zero(), |acc, &p| acc.max(p));
let keep: Vec<usize> = (0..self.bins.len()).filter(|&i| self.input_power[i] >= threshold).collect();
Self {
bins: keep.iter().map(|&i| self.bins[i]).collect(),
coherence: keep.iter().map(|&i| self.coherence[i]).collect(),
input_power: keep.iter().map(|&i| self.input_power[i]).collect(),
..self.clone()
}
}
/// Bins with `γ²(f) > bound`.
pub fn outside(&self) -> Vec<usize> {
self.bins.iter().zip(&self.coherence).filter(|(_, c)| **c > self.bound).map(|(&f, _)| f).collect()
}
/// Fraction of the bins with `γ²(f) > bound`; about `1 - confidence` for a right model.
pub fn fraction_outside(&self) -> T {
T::from(self.outside().len()).unwrap() / T::from(self.bins.len().max(1)).unwrap()
}
/// Mean of `γ²(f)` over the bins, to compare with `expected_coherence`.
pub fn mean_coherence(&self) -> T {
self.coherence.iter().fold(T::zero(), |acc, &c| acc + c) / T::from(self.bins.len().max(1)).unwrap()
}
/// Expectation `1 / L` of `γ²(f)` for a right model.
pub fn expected_coherence(&self) -> T {
T::one() / self.effective_segments
}
}