modewise_fc/
measures.rs

1use modewise_hilbert::{Padding, analytic_signal};
2use ndarray::{Array2, ArrayView2, Axis};
3
4/// Pearson correlation between the rows (channels) of a channel × sample array.
5/// A channel with zero variance gets NaN in its row and column.
6pub fn correlation(x: ArrayView2<f64>) -> Array2<f64> {
7    let n = x.ncols() as f64;
8    let mut z = x.to_owned();
9    let mut flat = vec![false; x.nrows()];
10    for (mut row, flat) in z.axis_iter_mut(Axis(0)).zip(&mut flat) {
11        let mean = row.sum() / n;
12        row.mapv_inplace(|v| v - mean);
13        let norm = row.dot(&row).sqrt();
14        if norm > 0.0 {
15            row /= norm;
16        } else {
17            *flat = true;
18        }
19    }
20    let mut r = z.dot(&z.t()).mapv(|v| v.clamp(-1.0, 1.0));
21    for i in 0..r.nrows() {
22        r[[i, i]] = 1.0;
23    }
24    mask_rows(&mut r, &flat);
25    r
26}
27
28/// Fisher z = atanh(r). The diagonal is undefined (r = 1) and becomes NaN.
29pub(crate) fn fisher_z(r: Array2<f64>) -> Array2<f64> {
30    let mut z = r.mapv(f64::atanh);
31    z.diag_mut().fill(f64::NAN);
32    z
33}
34
35/// Unit phase vectors `e^{iφ}` of each row's analytic signal, as real and imaginary parts,
36/// and which rows are all zero (no phase).
37pub fn unit_phases(x: ArrayView2<f64>, padding: Padding) -> (Array2<f64>, Array2<f64>, Vec<bool>) {
38    let (rows, n) = x.dim();
39    let (mut re, mut im) = (Array2::zeros((rows, n)), Array2::zeros((rows, n)));
40    let mut silent = vec![false; rows];
41    for (c, row) in x.axis_iter(Axis(0)).enumerate() {
42        let z = analytic_signal(row, padding);
43        silent[c] = z.iter().all(|v| v.norm() == 0.0);
44        for (t, v) in z.iter().enumerate() {
45            let a = v.norm();
46            if a > 0.0 {
47                re[[c, t]] = v.re / a;
48                im[[c, t]] = v.im / a;
49            }
50        }
51    }
52    (re, im, silent)
53}
54
55/// Phase locking value `|mean_t e^{i(φ_a − φ_b)}|` for every pair of rows of unit phase vectors.
56pub fn phase_locking(re: &ArrayView2<f64>, im: &ArrayView2<f64>, silent: &[bool]) -> Array2<f64> {
57    let n = re.ncols() as f64;
58    let real = re.dot(&re.t()) + im.dot(&im.t());
59    let imaginary = im.dot(&re.t()) - re.dot(&im.t());
60    let mut plv = Array2::from_shape_fn(real.dim(), |(i, j)| {
61        (real[[i, j]].hypot(imaginary[[i, j]]) / n).min(1.0)
62    });
63    mask_rows(&mut plv, silent);
64    plv
65}
66
67fn mask_rows(m: &mut Array2<f64>, masked: &[bool]) {
68    for (i, &bad) in masked.iter().enumerate() {
69        if bad {
70            m.row_mut(i).fill(f64::NAN);
71            m.column_mut(i).fill(f64::NAN);
72        }
73    }
74}