modewise_emd/
ensemble.rs

1//! Noise-assisted EMD: EEMD (Wu and Huang 2009) and CEEMDAN in its improved form (Colominas,
2//! Schlotthauer and Torres 2014), each following its authors' code (eemd.m, iceemdan.m).
3
4use modewise_decomposition::{Decomposition, ModeDecomposer};
5use modewise_method::{MethodDescriptor, MethodError};
6use modewise_signals::MultichannelSeries;
7use ndarray::{Array1, ArrayView1};
8use rand::SeedableRng;
9use rand_chacha::ChaCha8Rng;
10use rand_distr::{Distribution, StandardNormal};
11use rayon::prelude::*;
12
13use crate::decomposers::validate;
14use crate::layout::{Channel, emd1, per_channel};
15use crate::sift::Stop;
16
17fn sd(x: ArrayView1<f64>) -> f64 {
18    x.std(1.0)
19}
20
21fn mean_of(parts: Vec<Array1<f64>>, n: usize) -> Array1<f64> {
22    let count = parts.len() as f64;
23    parts.into_iter().fold(Array1::zeros(n), |a, b| a + b) / count
24}
25
26/// Unit-variance white noise, one realisation per trial; channel `c` has its own stream.
27fn white_noise(seed: u64, channel: usize, trials: usize, n: usize) -> Vec<Array1<f64>> {
28    let mut rng = ChaCha8Rng::seed_from_u64(seed);
29    rng.set_stream(channel as u64);
30    (0..trials)
31        .map(|_| Array1::from_shape_simple_fn(n, || StandardNormal.sample(&mut rng)))
32        .collect()
33}
34
35fn check_noise(trials: usize, scale: f64, name: &'static str) -> Result<(), MethodError> {
36    if trials == 0 {
37        return Err(MethodError::invalid("trials", "must be at least 1"));
38    }
39    if !(scale.is_finite() && scale > 0.0) {
40        return Err(MethodError::invalid(name, "must be finite and positive"));
41    }
42    Ok(())
43}
44
45/// EEMD of one channel with the given noise realisations (already scaled): the mean over trials
46/// of each IMF index and of the residual.
47pub(crate) fn eemd_channel(
48    x: ArrayView1<f64>,
49    noises: &[Array1<f64>],
50    stop: Stop,
51    max_imfs: usize,
52) -> Channel {
53    let n = x.len();
54    let trials: Vec<Channel> = noises
55        .par_iter()
56        .map(|w| emd1((&x + w).view(), stop, Some(max_imfs)))
57        .collect();
58    let count = trials.iter().map(|(imfs, _)| imfs.len()).max().unwrap_or(0);
59    let imfs = (0..count)
60        .map(|j| {
61            mean_of(
62                trials
63                    .iter()
64                    .map(|(imfs, _)| imfs.get(j).cloned().unwrap_or_else(|| Array1::zeros(n)))
65                    .collect(),
66                n,
67            )
68        })
69        .collect();
70    let residual = mean_of(trials.into_iter().map(|(_, r)| r).collect(), n);
71    (imfs, residual)
72}
73
74/// Ensemble EMD (Wu and Huang 2009), each channel on its own: the mean of EMDs of the signal
75/// plus white noise with SD `noise_sd`·SD(signal). Every trial takes at most ⌊log₂ n⌋ − 1 IMFs
76/// unless set, as Wu and Huang's eemd.m; unlike it, a trial stops early once its residual no
77/// longer oscillates. Modes sum to the signal plus the mean noise, which shrinks with more trials.
78#[derive(Debug, Clone)]
79pub struct Eemd {
80    pub stop: Stop,
81    pub max_imfs: Option<usize>,
82    pub trials: usize,
83    pub noise_sd: f64,
84    pub seed: u64,
85}
86
87impl Eemd {
88    /// Wu and Huang's settings: 10 siftings per IMF, 100 trials, noise SD 0.2.
89    pub fn new(trials: usize, noise_sd: f64, seed: u64) -> Result<Self, MethodError> {
90        Self {
91            stop: Stop::Fixed { siftings: 10 },
92            max_imfs: None,
93            trials,
94            noise_sd,
95            seed,
96        }
97        .validated()
98    }
99
100    pub fn validated(self) -> Result<Self, MethodError> {
101        validate(self.stop, self.max_imfs)?;
102        check_noise(self.trials, self.noise_sd, "noise_sd")?;
103        Ok(self)
104    }
105}
106
107impl ModeDecomposer for Eemd {
108    fn descriptor(&self) -> MethodDescriptor {
109        MethodDescriptor {
110            name: "eemd",
111            version: env!("CARGO_PKG_VERSION"),
112            multivariate: false,
113        }
114    }
115
116    fn decompose(&self, signal: &MultichannelSeries) -> Result<Decomposition, MethodError> {
117        let n = signal.n_samples();
118        let k = self
119            .max_imfs
120            .unwrap_or_else(|| (n.ilog2() as usize).saturating_sub(1).max(1));
121        per_channel(signal, |c, x| {
122            let scale = self.noise_sd * sd(x);
123            let noises: Vec<_> = white_noise(self.seed, c, self.trials, n)
124                .into_iter()
125                .map(|w| w * scale)
126                .collect();
127            eemd_channel(x, &noises, self.stop, k)
128        })
129    }
130}
131
132/// The local mean M(y): what is left after the first IMF (y itself when it has none).
133fn local_mean(y: ArrayView1<f64>, stop: Stop) -> Array1<f64> {
134    emd1(y, stop, Some(1)).1
135}
136
137fn has_imf(y: ArrayView1<f64>, stop: Stop) -> bool {
138    !emd1(y, stop, Some(1)).0.is_empty()
139}
140
141/// CEEMDAN of one channel with unit-variance noise realisations, as iceemdan.m (SNRFlag 1).
142pub(crate) fn ceemdan_channel(
143    x: ArrayView1<f64>,
144    noises: &[Array1<f64>],
145    stop: Stop,
146    epsilon: f64,
147    max_imfs: Option<usize>,
148) -> Channel {
149    let n = x.len();
150    let scale = sd(x);
151    if scale == 0.0 || !scale.is_finite() {
152        return (Vec::new(), x.to_owned());
153    }
154    let x = &x / scale;
155    // Each realisation's EMD, with its residual as the last row, as MATLAB emd returns it.
156    let noise_modes: Vec<Vec<Array1<f64>>> = noises
157        .par_iter()
158        .map(|w| {
159            let (mut imfs, residual) = emd1(w.view(), stop, None);
160            imfs.push(residual);
161            imfs
162        })
163        .collect();
164    let first: Vec<Array1<f64>> = noise_modes
165        .par_iter()
166        .map(|e| {
167            let xi = &x + &(&e[0] * (epsilon / sd(e[0].view())));
168            local_mean(xi.view(), stop)
169        })
170        .collect();
171    let mut r = mean_of(first, n);
172    let mut modes = vec![&x - &r];
173    let mut k = 1;
174    while max_imfs.is_none_or(|m| modes.len() < m) && has_imf(r.view(), stop) {
175        let beta = epsilon * sd(r.view());
176        let means: Vec<Array1<f64>> = noise_modes
177            .par_iter()
178            .map(|e| match e.get(k) {
179                Some(noise) => local_mean((&r + &(noise * beta)).view(), stop),
180                None => local_mean(r.view(), stop),
181            })
182            .collect();
183        let next = mean_of(means, n);
184        modes.push(&r - &next);
185        r = next;
186        k += 1;
187    }
188    (modes.into_iter().map(|m| m * scale).collect(), r * scale)
189}
190
191/// Complete ensemble EMD with adaptive noise, improved form (Colominas et al. 2014), each channel
192/// on its own. Mode k is r_{k−1} − ⟨M(r_{k−1} + ε·SD(r_{k−1})·E_k(w))⟩, so modes sum to the
193/// signal exactly.
194#[derive(Debug, Clone)]
195pub struct Ceemdan {
196    pub stop: Stop,
197    pub max_imfs: Option<usize>,
198    pub trials: usize,
199    pub epsilon: f64,
200    pub seed: u64,
201}
202
203impl Ceemdan {
204    /// The authors' settings: emd.m's Rilling stop, ε = 0.2.
205    pub fn new(trials: usize, epsilon: f64, seed: u64) -> Result<Self, MethodError> {
206        Self {
207            stop: Stop::RILLING_EMD,
208            max_imfs: None,
209            trials,
210            epsilon,
211            seed,
212        }
213        .validated()
214    }
215
216    pub fn validated(self) -> Result<Self, MethodError> {
217        validate(self.stop, self.max_imfs)?;
218        check_noise(self.trials, self.epsilon, "epsilon")?;
219        Ok(self)
220    }
221}
222
223impl ModeDecomposer for Ceemdan {
224    fn descriptor(&self) -> MethodDescriptor {
225        MethodDescriptor {
226            name: "ceemdan",
227            version: env!("CARGO_PKG_VERSION"),
228            multivariate: false,
229        }
230    }
231
232    fn decompose(&self, signal: &MultichannelSeries) -> Result<Decomposition, MethodError> {
233        let n = signal.n_samples();
234        per_channel(signal, |c, x| {
235            let noises = white_noise(self.seed, c, self.trials, n);
236            ceemdan_channel(x, &noises, self.stop, self.epsilon, self.max_imfs)
237        })
238    }
239}
modewise_emd/ensemble.rs

modewise_emd/
ensemble.rs
