modewise_mvmd/
vmd.rs

1use modewise_decomposition::{Decomposition, ModeDecomposer};
2use modewise_method::{MethodDescriptor, MethodError};
3use modewise_signals::MultichannelSeries;
4use ndarray::{Array2, Array3, Axis};
5use rayon::prelude::*;
6
7use crate::{Mvmd, MvmdParams};
8
9/// Variational Mode Decomposition of each channel on its own: MVMD with one channel, which is
10/// VMD (Rehman and Aftab 2019, section III; same grid, mirror and update as vmdpy). Mode k is
11/// the k-th lowest centre frequency of its channel, so modes are not aligned across channels.
12#[derive(Debug, Clone)]
13pub struct Vmd {
14    inner: Mvmd,
15}
16
17impl Vmd {
18    pub fn new(params: MvmdParams) -> Result<Self, MethodError> {
19        Ok(Self {
20            inner: Mvmd::new(params)?,
21        })
22    }
23}
24
25impl ModeDecomposer for Vmd {
26    fn descriptor(&self) -> MethodDescriptor {
27        MethodDescriptor {
28            name: "vmd",
29            version: env!("CARGO_PKG_VERSION"),
30            multivariate: false,
31        }
32    }
33
34    fn decompose(&self, signal: &MultichannelSeries) -> Result<Decomposition, MethodError> {
35        let fs = signal.sampling_rate();
36        let rows: Vec<_> = signal
37            .data()
38            .axis_iter(Axis(0))
39            .map(|r| r.to_owned())
40            .collect();
41        let channels = rows
42            .into_par_iter()
43            .map(|row| {
44                let single =
45                    MultichannelSeries::new(row.insert_axis(Axis(0)), fs, vec!["x".into()])
46                        .map_err(|e| MethodError::ShapeMismatch(e.to_string()))?;
47                self.inner.decompose(&single)
48            })
49            .collect::<Result<Vec<_>, _>>()?;
50        let k = self.inner.params().num_modes;
51        let mut modes = Array3::zeros((k, signal.n_channels(), signal.n_samples()));
52        let mut centres = Array2::zeros((k, signal.n_channels()));
53        for (c, d) in channels.iter().enumerate() {
54            modes
55                .index_axis_mut(Axis(1), c)
56                .assign(&d.modes().index_axis(Axis(1), 0));
57            centres
58                .column_mut(c)
59                .assign(&ndarray::ArrayView1::from(d.centre_frequencies_hz()));
60        }
61        Decomposition::per_channel(modes, centres, fs, signal.labels().to_vec())
62    }
63}