modewise_emd/
decomposers.rs

1use modewise_decomposition::{Decomposition, ModeDecomposer};
2use modewise_method::{MethodDescriptor, MethodError};
3use modewise_signals::MultichannelSeries;
4use ndarray::{Array3, Axis};
5use tracing::warn;
6
7use crate::directions::directions;
8use crate::layout::{emd1, per_channel};
9use crate::sift::{Stop, decompose};
10use modewise_fourier::centroid_hz;
11
12pub(crate) fn validate(stop: Stop, max_imfs: Option<usize>) -> Result<(), MethodError> {
13    match stop {
14        Stop::Rilling { sd, sd2, tol } => {
15            if ![sd, sd2, tol].iter().all(|v| v.is_finite() && *v > 0.0) {
16                return Err(MethodError::invalid(
17                    "stop",
18                    "rilling thresholds must be finite and positive",
19                ));
20            }
21        }
22        Stop::SNumber { count: 0 } | Stop::Fixed { siftings: 0 } => {
23            return Err(MethodError::invalid("stop", "the count must be at least 1"));
24        }
25        _ => {}
26    }
27    if max_imfs == Some(0) {
28        return Err(MethodError::invalid("max_imfs", "must be at least 1"));
29    }
30    Ok(())
31}
32
33/// Empirical mode decomposition (Huang et al. 1998) of each channel on its own.
34///
35/// Channels can give different numbers of IMFs. Modes are aligned by IMF index from the top:
36/// the last mode is every channel's IMF 1 (highest frequency), mode 0 is the residual, and a
37/// channel with fewer IMFs has zero modes (NaN frequency) just above its residual.
38#[derive(Debug, Clone)]
39pub struct Emd {
40    stop: Stop,
41    max_imfs: Option<usize>,
42}
43
44impl Emd {
45    pub fn new(stop: Stop, max_imfs: Option<usize>) -> Result<Self, MethodError> {
46        validate(stop, max_imfs)?;
47        Ok(Self { stop, max_imfs })
48    }
49}
50
51impl ModeDecomposer for Emd {
52    fn descriptor(&self) -> MethodDescriptor {
53        MethodDescriptor {
54            name: "emd",
55            version: env!("CARGO_PKG_VERSION"),
56            multivariate: false,
57        }
58    }
59
60    fn decompose(&self, signal: &MultichannelSeries) -> Result<Decomposition, MethodError> {
61        per_channel(signal, |_, x| emd1(x, self.stop, self.max_imfs))
62    }
63}
64
65/// Multivariate EMD (Rehman and Mandic 2010): envelopes from projections of all channels along
66/// `num_directions` Hammersley directions, so every channel gets the same IMFs.
67#[derive(Debug, Clone)]
68pub struct Memd {
69    stop: Stop,
70    max_imfs: Option<usize>,
71    num_directions: usize,
72    memd_flat_rule: bool,
73}
74
75impl Memd {
76    pub fn new(
77        stop: Stop,
78        max_imfs: Option<usize>,
79        num_directions: usize,
80    ) -> Result<Self, MethodError> {
81        validate(stop, max_imfs)?;
82        if num_directions < 6 {
83            return Err(MethodError::invalid(
84                "num_directions",
85                "must be at least 6, as memd.m requires",
86            ));
87        }
88        Ok(Self {
89            stop,
90            max_imfs,
91            num_directions,
92            memd_flat_rule: false,
93        })
94    }
95}
96
97impl Memd {
98    /// Reproduces memd.m exactly, including its rule that a projection with every value below
99    /// 1e-5 has no extrema. For comparisons with memd.m results only.
100    #[doc(hidden)]
101    pub fn with_memd_flat_rule(mut self) -> Self {
102        self.memd_flat_rule = true;
103        self
104    }
105}
106
107impl ModeDecomposer for Memd {
108    fn descriptor(&self) -> MethodDescriptor {
109        MethodDescriptor {
110            name: "memd",
111            version: env!("CARGO_PKG_VERSION"),
112            multivariate: true,
113        }
114    }
115
116    fn decompose(&self, signal: &MultichannelSeries) -> Result<Decomposition, MethodError> {
117        let (n_channels, n) = signal.data().dim();
118        if n_channels < 3 {
119            return Err(MethodError::invalid(
120                "signal",
121                "MEMD needs at least 3 channels, as memd.m; use emd for fewer",
122            ));
123        }
124        if self.num_directions < 2 * n_channels {
125            warn!(
126                directions = self.num_directions,
127                channels = n_channels,
128                "fewer directions than twice the channels; envelopes may be poorly sampled"
129            );
130        }
131        let data = signal.data();
132        let x = data.t();
133        let (imfs, residual) = decompose(
134            x,
135            &directions(self.num_directions, n_channels),
136            self.stop,
137            self.max_imfs,
138            self.memd_flat_rule,
139        );
140        let fs = signal.sampling_rate();
141        let mut modes = Array3::zeros((imfs.len() + 1, n_channels, n));
142        modes.index_axis_mut(Axis(0), 0).assign(&residual.t());
143        for (j, imf) in imfs.iter().enumerate() {
144            modes
145                .index_axis_mut(Axis(0), imfs.len() - j)
146                .assign(&imf.t());
147        }
148        let centres = modes
149            .axis_iter(Axis(0))
150            .map(|m| centroid_hz(m, fs.hz()))
151            .collect();
152        Decomposition::new(modes, centres, fs, signal.labels().to_vec(), true)
153    }
154}
modewise_emd/decomposers.rs

modewise_emd/
decomposers.rs
