extractSSTModes: Components Along the Ridges of a Synchrosqueezed Transform

In xemd: eXtended Empirical Mode Decomposition

Components Along the Ridges of a Synchrosqueezed Transform

Description

Extracts ridges of a synchrosqueezed transform one after another by dynamic programming and reconstructs a component from the coefficients in a band around each ridge.

Usage

extractSSTModes(transform, frequency, modeCount, penalty = 1000,
                exclusionWidth = 4, reconstructionWidth, morletCenter)

Arguments

transform	Numeric or complex matrix of finite values, samples by frequency cells, such as the transform of transformSST.
frequency	Increasing positive frequencies of the cells, one per column of transform, such as the frequency of transformSST.
modeCount	Whole number \ge 1 of ridges. No default; curve_ext_multi of the Synchrosqueezing Toolbox requires the value.
penalty	Dimensionless penalty \lambda \ge 0 of a jump between cells; default 1000, the value that curve_ext_multi assigns.
exclusionWidth	Whole number b \ge 0 of cells on each side of a ridge that are removed before the next extraction; default 4, the value that curve_ext_multi assigns.
reconstructionWidth	Whole number w \ge 0 of cells on each side of a ridge in its reconstruction band. No default; a value not below the number of cells gives the complete inversion, as synsq_cwt_iw does without a band.
morletCenter	Center \mu > 0 of the Morlet wavelet that produced transform, in normalized frequency. No default.

Details

With T_{n,j} the coefficients, f_j the frequencies and z_j = \log_2 f_j, the cost of a cell is E_{n,j} = -\log(|T_{n,j}|^2 / Q + 10^{-8}), Q the sum of |T|^2 over the matrix. The accumulated cost is F_{1,j} = E_{1,j} and

F_{n,j} = E_{n,j} + \min_k \left[F_{n-1,k} + \lambda (z_j - z_k)^2\right]

over every cell k; the ridge ends at the first cell of least F_{N,j}, and each earlier cell is the first whose residual in that recurrence is below 10^{-8} in modulus. After each ridge C_n, the cells from C_n - b to C_n + b, clipped to the cells, take the value \sqrt{\epsilon}, \epsilon the machine precision, and the next ridge is extracted with the energy recomputed. Component k sums the original coefficients of its band,

c_k[n] = \frac{1}{C} \mathrm{Re} \sum_{j = C_{k,n} - w}^{C_{k,n} + w} T_{n,j},

over the cells of the band, with C = C_\psi / (\sqrt{2\pi} \, 2 \ln 2) and C_\psi = \int_0^\infty \hat\psi(\xi) / \xi \, d\xi for the Morlet wavelet of center morletCenter. The bands may overlap, the components need not add up to the signal, and no residue is returned. The components are defined by the selected ridge bands, without imposing the IMF conditions of EMD.

Value

Double matrix, samples by components, one column per ridge.

Conditions

A transform of zero energy stops the function with an error of class xemdZeroEnergy. A sample where no cell satisfies the recurrence within the tolerance stops it with an error.

References

Daubechies, I., Lu, J. and Wu, H.-T. (2009). Synchrosqueezed wavelet transforms: a tool for empirical mode decomposition. arXiv:0912.2437.

Brevdo, E. and Thakur, G. Synchrosqueezing Toolbox, revision 7e9fec0c6c9ed478dafac4479c0658d24fc6d8ef, functions curve_ext (original code by Jianfeng Lu), curve_ext_multi and synsq_cwt_iw, https://github.com/ebrevdo/synchrosqueezing.

The extraction of the ridges is adapted from curve_ext (original author Jianfeng Lu, maintainer Eugene Brevdo) and curve_ext_multi of the Synchrosqueezing Toolbox, as ported in SynchWave 1.1.2, revision 9ac80b2c3e60b417d0ee5567a66e8e64faa1a979 (port by Dongik Jang), and the inversion from SynchWave 1.1.2 (MATLAB original by Eugene Brevdo; R port by Dongik Jang, Hee-Seok Oh and Donghoh Kim), LGPL-2.0-or-later, which xemd distributes under GPL-3 as section 3 of the GNU Library General Public License version 2 permits; see the NOTICE file of the installed package.

See Also

transformSST, decomposeEWT

Examples

SamplingInterval <- 1 / 1000
Time <- seq_len(1000) * SamplingInterval
Signal <- cos(2 * pi * 2 * Time) + 0.25 * cos(2 * pi * 24 * Time) +
    0.0625 * cos(2 * pi * 288 * Time)
SST <- transformSST(Signal, samplingInterval = SamplingInterval,
                    voiceCount = 16, morletCenter = 2 * pi,
                    threshold = sqrt(.Machine$double.eps))
Modes <- extractSSTModes(SST$transform, frequency = SST$frequency,
                         modeCount = 3, reconstructionWidth = 0,
                         morletCenter = 2 * pi)
dim(Modes)
sqrt(colMeans(Modes^2))

xemd documentation built on Nov. 9, 2026, 1:09 a.m.