Components Along the Ridges of a Synchrosqueezed Transform

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.

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)
#> [1] 1000    3
sqrt(colMeans(Modes^2))
#> [1] 0.46084922 0.17525752 0.04381498