modewise_hilbert/
analytic.rs

1use std::f64::consts::PI;
2
3use modewise_decomposition::{Decomposition, InstantaneousModes};
4use modewise_method::MethodError;
5use ndarray::{Array3, ArrayView1, Axis};
6use rustfft::FftPlanner;
7use rustfft::num_complex::Complex64;
8
9/// How the signal is extended before the FFT, which treats it as periodic.
10#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
11pub enum Padding {
12    /// Mirror-extend to twice the length, which avoids a wrap-around jump.
13    #[default]
14    Mirror,
15    /// Use the signal as it is (matches `scipy.signal.hilbert`).
16    None,
17}
18
19/// Analytic signal `x + i·H{x}` by the FFT method: keep DC (and Nyquist for even lengths),
20/// double positive frequencies, zero negative ones.
21pub fn analytic_signal(x: ArrayView1<f64>, padding: Padding) -> Vec<Complex64> {
22    let n = x.len();
23    if n == 0 {
24        return Vec::new();
25    }
26    let (extended, offset): (Vec<f64>, usize) = match padding {
27        Padding::None => (x.to_vec(), 0),
28        // The signal followed by its reverse is periodic without a jump. Where the reversed
29        // part goes does not matter: the FFT method is unaffected by rotating the sequence.
30        Padding::Mirror => (x.iter().chain(x.iter().rev()).copied().collect(), 0),
31    };
32    let m = extended.len();
33    let mut planner = FftPlanner::new();
34    let mut buffer: Vec<Complex64> = extended.iter().map(|&v| Complex64::new(v, 0.0)).collect();
35    planner.plan_fft_forward(m).process(&mut buffer);
36    for (k, value) in buffer.iter_mut().enumerate() {
37        let gain = if k == 0 || (m % 2 == 0 && k == m / 2) {
38            1.0
39        } else if k < m.div_ceil(2) {
40            2.0
41        } else {
42            0.0
43        };
44        *value *= gain;
45    }
46    planner.plan_fft_inverse(m).process(&mut buffer);
47    let scale = 1.0 / m as f64;
48    buffer[offset..offset + n]
49        .iter()
50        .map(|z| z * scale)
51        .collect()
52}
53
54/// Instantaneous frequency in Hz from an analytic signal: centred difference of the unwrapped
55/// phase, `(φ[n+1] − φ[n−1]) · fs / 4π`; at each end, the slope of a line through 3 points (the
56/// same estimator as a 3-point Savitzky-Golay derivative). Phase steps are taken one sample at
57/// a time, so frequencies up to Nyquist are estimated without wrapping.
58pub fn instantaneous_frequency(z: &[Complex64], fs: f64) -> Vec<f64> {
59    let n = z.len();
60    let to_hz = |radians: f64| radians * fs / (2.0 * PI);
61    let step: Vec<f64> = z.windows(2).map(|w| (w[1] * w[0].conj()).arg()).collect();
62    match n {
63        1 => vec![0.0],
64        2 => vec![to_hz(step[0]); 2],
65        _ => (0..n)
66            .map(|i| {
67                let first_step = i.clamp(1, n - 2) - 1;
68                to_hz(step[first_step] + step[first_step + 1]) / 2.0
69            })
70            .collect(),
71    }
72}
73
74/// Instantaneous amplitude and frequency of every mode and channel.
75pub fn instantaneous(
76    modes: &Decomposition,
77    padding: Padding,
78) -> Result<InstantaneousModes, MethodError> {
79    let fs = modes.sampling_rate().hz();
80    let dim = modes.modes().dim();
81    let mut amplitude = Array3::zeros(dim);
82    let mut frequency = Array3::zeros(dim);
83    for (k, mode) in modes.modes().axis_iter(Axis(0)).enumerate() {
84        for (c, signal) in mode.axis_iter(Axis(0)).enumerate() {
85            let z = analytic_signal(signal, padding);
86            for (t, (value, f)) in z.iter().zip(instantaneous_frequency(&z, fs)).enumerate() {
87                amplitude[[k, c, t]] = value.norm();
88                frequency[[k, c, t]] = f;
89            }
90        }
91    }
92    InstantaneousModes::new(
93        amplitude,
94        frequency,
95        modes.centre_frequencies_hz().to_vec(),
96        modes.sampling_rate(),
97        modes.labels().to_vec(),
98    )
99}
100
101#[cfg(test)]
102mod tests {
103    use ndarray::Array1;
104
105    use super::*;
106
107    #[test]
108    fn analytic_signal_of_cosine_is_complex_exponential() {
109        for n in [64usize, 65] {
110            let k = 5.0;
111            let x = Array1::from_shape_fn(n, |t| (2.0 * PI * k * t as f64 / n as f64).cos());
112            let z = analytic_signal(x.view(), Padding::None);
113            for (t, value) in z.iter().enumerate() {
114                let expected = Complex64::from_polar(1.0, 2.0 * PI * k * t as f64 / n as f64);
115                assert!((value - expected).norm() < 1e-12, "n={n} t={t}");
116            }
117        }
118    }
119
120    #[test]
121    fn edge_cases() {
122        assert!(analytic_signal(Array1::<f64>::zeros(0).view(), Padding::Mirror).is_empty());
123        assert_eq!(instantaneous_frequency(&[], 1.0), Vec::<f64>::new());
124        assert_eq!(
125            instantaneous_frequency(&[Complex64::new(1.0, 0.0)], 1.0),
126            [0.0]
127        );
128        let quarter = [Complex64::new(1.0, 0.0), Complex64::new(0.0, 1.0)];
129        assert_eq!(
130            instantaneous_frequency(&quarter, 4.0),
131            [1.0, 1.0],
132            "a quarter turn per sample at 4 Hz is 1 Hz"
133        );
134    }
135}
modewise_hilbert/analytic.rs

modewise_hilbert/
analytic.rs
