1use modewise_hilbert::{Padding, analytic_signal};
2use ndarray::{Array2, ArrayView2, Axis};
3
4pub 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
28pub(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
35pub 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
55pub 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}