modewise_hilbert/
spectrum.rs

1use modewise_decomposition::InstantaneousModes;
2use modewise_method::MethodError;
3use modewise_spectral::{FrequencyGrid, TimeFrequencyMap};
4use ndarray::Array3;
5
6/// What each sample contributes to its frequency bin.
7#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
8pub enum Weight {
9    /// Amplitude squared (the `emd` package default).
10    #[default]
11    Power,
12    Amplitude,
13}
14
15pub struct HilbertSpectrum {
16    pub map: TimeFrequencyMap,
17    /// Samples whose instantaneous frequency fell outside the grid (or was not finite).
18    pub dropped: usize,
19    pub total: usize,
20}
21
22/// Hilbert spectrum, channel × frequency × sample: each sample's weight is added to the bin
23/// holding its instantaneous frequency, summed over modes (Huang et al. 1998).
24pub fn hilbert_spectrum(
25    modes: &InstantaneousModes,
26    grid: &FrequencyGrid,
27    weight: Weight,
28) -> Result<HilbertSpectrum, MethodError> {
29    let (n_modes, n_channels, n_samples) = modes.amplitude().dim();
30    let mut power = Array3::zeros((n_channels, grid.n_bins(), n_samples));
31    let amplitude = modes.amplitude();
32    let frequency = modes.frequency_hz();
33    let mut dropped = 0;
34    for k in 0..n_modes {
35        for c in 0..n_channels {
36            for t in 0..n_samples {
37                match grid.bin(frequency[[k, c, t]]) {
38                    Some(bin) => {
39                        let a = amplitude[[k, c, t]];
40                        power[[c, bin, t]] += match weight {
41                            Weight::Power => a * a,
42                            Weight::Amplitude => a,
43                        };
44                    }
45                    None => dropped += 1,
46                }
47            }
48        }
49    }
50    let fs = modes.sampling_rate().hz();
51    let times = (0..n_samples).map(|t| t as f64 / fs).collect();
52    let map = TimeFrequencyMap::new(power, grid.centres(), times, modes.labels().to_vec())?
53        .with_sampling_rate(modes.sampling_rate())?;
54    Ok(HilbertSpectrum {
55        map,
56        dropped,
57        total: n_modes * n_channels * n_samples,
58    })
59}
modewise_hilbert/spectrum.rs

modewise_hilbert/
spectrum.rs
