modewise_emd/
sift.rs

1//! The sifting engine shared by EMD and MEMD, following memd.m: envelopes of projections along
2//! direction vectors, averaged; one channel with direction [1] is plain EMD.
3
4use ndarray::{Array1, Array2, ArrayView2, Axis};
5
6use crate::boundary::mirror;
7use crate::extrema::{local_peaks, zero_crossings};
8use crate::spline::interpolate;
9
10/// When sifting one mode stops.
11#[derive(Debug, Clone, Copy, PartialEq)]
12pub enum Stop {
13    /// Rilling, Flandrin and Gonçalvès (2003): with σ = |envelope mean| / mode amplitude, stop
14    /// once σ < `sd` on all but a fraction `tol` of samples and σ < `sd2` everywhere.
15    Rilling { sd: f64, sd2: f64, tol: f64 },
16    /// Huang et al. (2003) S-number: stop after `count` consecutive siftings in which the
17    /// extrema and zero crossings differ by at most one (memd.m `fix_h`).
18    SNumber { count: usize },
19    /// A fixed number of siftings per mode (Wu and Huang 2009).
20    Fixed { siftings: usize },
21}
22
23impl Stop {
24    /// emd.m's default.
25    pub const RILLING_EMD: Self = Self::Rilling {
26        sd: 0.05,
27        sd2: 0.5,
28        tol: 0.05,
29    };
30    /// memd.m's default.
31    pub const RILLING_MEMD: Self = Self::Rilling {
32        sd: 0.075,
33        sd2: 0.75,
34        tol: 0.075,
35    };
36}
37
38/// Limits from the reference codes.
39pub(crate) const MAX_SIFTINGS: usize = 1000;
40
41pub(crate) struct Envelope {
42    mean: Array2<f64>,
43    amplitude: Array1<f64>,
44    extrema: Vec<usize>,
45    zeros: Vec<usize>,
46}
47
48fn project(m: ArrayView2<f64>, direction: &[f64]) -> Vec<f64> {
49    m.rows()
50        .into_iter()
51        .map(|row| row.iter().zip(direction).map(|(a, b)| a * b).sum())
52        .collect()
53}
54
55/// Mean envelope over all directions with enough extrema, or `None` when an envelope cannot be
56/// built (the reference codes then stop sifting with a zero mean).
57fn envelope(m: ArrayView2<f64>, directions: &[Vec<f64>], flat_rule: bool) -> Option<Envelope> {
58    let (n, dims) = m.dim();
59    let mut mean = Array2::zeros((n, dims));
60    let mut amplitude = Array1::zeros(n);
61    let (mut extrema, mut zeros) = (
62        Vec::with_capacity(directions.len()),
63        Vec::with_capacity(directions.len()),
64    );
65    let mut used = 0usize;
66    for direction in directions {
67        let y = project(m, direction);
68        let (indmin, indmax) = local_peaks(&y, flat_rule);
69        extrema.push(indmin.len() + indmax.len());
70        zeros.push(zero_crossings(&y));
71        if indmin.len() + indmax.len() < 3 {
72            continue;
73        }
74        let (lower, upper) = mirror(&indmin, &indmax, &y).ok()?;
75        let lo = interpolate(
76            &lower.positions,
77            m.select(Axis(0), &lower.samples).view(),
78            n,
79        );
80        let hi = interpolate(
81            &upper.positions,
82            m.select(Axis(0), &upper.samples).view(),
83            n,
84        );
85        let half = (&hi - &lo) / 2.0;
86        amplitude += &half.map_axis(Axis(1), |r| r.dot(&r).sqrt());
87        mean += &((&hi + &lo) / 2.0);
88        used += 1;
89    }
90    if used == 0 {
91        return Some(Envelope {
92            mean: Array2::zeros((n, dims)),
93            amplitude: Array1::zeros(n),
94            extrema: vec![0; directions.len()],
95            zeros,
96        });
97    }
98    mean /= used as f64;
99    amplitude /= used as f64;
100    Some(Envelope {
101        mean,
102        amplitude,
103        extrema,
104        zeros,
105    })
106}
107
108/// Whether to stop sifting `m`, and the envelope mean to subtract if not.
109struct Check {
110    stop: bool,
111    mean: Array2<f64>,
112}
113
114struct StopState {
115    rule: Stop,
116    counter: usize,
117    siftings: usize,
118}
119
120impl StopState {
121    fn new(rule: Stop) -> Self {
122        Self {
123            rule,
124            counter: 0,
125            siftings: 0,
126        }
127    }
128
129    fn check(&mut self, m: ArrayView2<f64>, directions: &[Vec<f64>], flat_rule: bool) -> Check {
130        let Some(env) = envelope(m, directions, flat_rule) else {
131            return Check {
132                stop: true,
133                mean: Array2::zeros(m.dim()),
134            };
135        };
136        let stop = match self.rule {
137            Stop::Rilling { sd, sd2, tol } => {
138                let mut sx = env.mean.map_axis(Axis(1), |r| r.dot(&r).sqrt());
139                // As the reference: scale only when every sample has an amplitude.
140                if env.amplitude.iter().all(|&a| a != 0.0) {
141                    sx /= &env.amplitude;
142                }
143                let n = sx.len() as f64;
144                let above = sx.iter().filter(|&&v| v > sd).count() as f64 / n;
145                let wild = sx.iter().any(|&v| v > sd2);
146                !((above > tol || wild) && env.extrema.iter().any(|&e| e > 2))
147            }
148            Stop::SNumber { count } => {
149                let unsettled = env
150                    .extrema
151                    .iter()
152                    .zip(&env.zeros)
153                    .all(|(&e, &z)| e.abs_diff(z) > 1);
154                if unsettled {
155                    self.counter = 0;
156                    false
157                } else {
158                    self.counter += 1;
159                    self.counter >= count
160                }
161            }
162            Stop::Fixed { siftings } => self.siftings >= siftings,
163        };
164        Check {
165            stop,
166            mean: env.mean,
167        }
168    }
169}
170
171/// True when every projection has fewer than three extrema: what is left is the residual.
172fn is_residual(r: ArrayView2<f64>, directions: &[Vec<f64>], flat_rule: bool) -> bool {
173    directions.iter().all(|d| {
174        let (indmin, indmax) = local_peaks(&project(r, d), flat_rule);
175        indmin.len() + indmax.len() < 3
176    })
177}
178
179/// Intrinsic mode functions (samples × dims each), highest frequency first, and the residual.
180pub(crate) fn decompose(
181    x: ArrayView2<f64>,
182    directions: &[Vec<f64>],
183    rule: Stop,
184    max_imfs: Option<usize>,
185    flat_rule: bool,
186) -> (Vec<Array2<f64>>, Array2<f64>) {
187    let scale = x.iter().fold(0.0f64, |a, v| a.max(v.abs()));
188    let mut r = x.to_owned();
189    let mut imfs = Vec::new();
190    while max_imfs.is_none_or(|k| imfs.len() < k) && !is_residual(r.view(), directions, flat_rule) {
191        let mut m = r.clone();
192        let mut state = StopState::new(rule);
193        let mut check = state.check(m.view(), directions, flat_rule);
194        // A mode this small would only add spurious extrema from rounding.
195        if m.iter().fold(0.0f64, |a, v| a.max(v.abs())) < 1e-10 * scale {
196            break;
197        }
198        while !check.stop && state.siftings < MAX_SIFTINGS {
199            m -= &check.mean;
200            state.siftings += 1;
201            check = state.check(m.view(), directions, flat_rule);
202        }
203        r -= &m;
204        imfs.push(m);
205    }
206    (imfs, r)
207}
208
209#[cfg(test)]
210mod tests {
211    use ndarray::Array2;
212
213    use super::*;
214
215    fn column(f: impl Fn(f64) -> f64) -> Array2<f64> {
216        Array2::from_shape_fn((200, 1), |(t, _)| f(t as f64))
217    }
218
219    #[test]
220    fn s_number_counts_only_siftings_with_extrema_and_zero_crossings_within_one() {
221        let up = vec![vec![1.0]];
222        // Riding on an offset: many extrema, no zero crossings, so not settled.
223        let riding = column(|t| (0.3 * t).sin() + 3.0);
224        let mut state = StopState::new(Stop::SNumber { count: 1 });
225        assert!(!state.check(riding.view(), &up, false).stop);
226        assert_eq!(state.counter, 0);
227        // Centred: extrema and zero crossings agree.
228        let centred = column(|t| (0.3 * t).sin());
229        assert!(state.check(centred.view(), &up, false).stop);
230        let mut two = StopState::new(Stop::SNumber { count: 2 });
231        assert!(!two.check(centred.view(), &up, false).stop);
232        assert!(two.check(centred.view(), &up, false).stop);
233    }
234}
modewise_emd/sift.rs

modewise_emd/
sift.rs
