Ridge detection for nonstationary multicomponent signals with time-varying wave-shape functions and its applications

We consider the adaptive non-harmonic model (ANHM) to model a non-sinusoidal oscillatory signal. Fix a small constant $\epsilon>0$ and $L\in\mathbb{N}$. Assume the signal satisfies

where $a_{\ell}(t)s_{\ell}(\phi_{\ell}(t))$ is called the $\ell$-th intrinsic model type (IMT) function that satisfies

• (C1)

  $\phi_{\ell}(t)$ is a $C^{2}$ function called the phase function of the $\ell$-th IMT function. We assume $\phi^{\prime}_{\ell}>0$;

• (C2)

  $\phi_{\ell}^{\prime}(t)$ is the instantaneous frequency (IF) of the $\ell$-th IMT function so that $|\phi^{\prime\prime}_{\ell}(t)|\leq\epsilon\phi^{\prime}_{\ell}(t)$ for all $t\in\mathbb{R}$ and $\inf_{t\in\mathbb{R}}\phi_{1}^{\prime}(t)>\Delta$ for some $\Delta>0$. When $L>1$, we assume $\phi_{j+1}^{\prime}(t)-\phi_{j}^{\prime}(t)>\Delta>0$ for all $j=1,\ldots,L-1$ and $t\in\mathbb{R}$;

• (C3)

  $a_{\ell}(t)>0$ is a $C^{1}$ function denoting the amplitude modulation (AM) of the $\ell$-th IMT function so that $|a_{\ell}^{\prime}(t)|\leq\epsilon\phi^{\prime}_{\ell}(t)$ for all $t\in\mathbb{R}$;

• (C4)

  $s_{\ell}$ is a smooth $1$-period function with mean 0 and unit $L^{2}$ norm so that $|\hat{s}_{\ell}(1)|>0$, which is called the wave-shape function (WSF) of the $\ell$-th IMT function;

• (C5)

  $\Phi(\cdot)$ is a mean $0$ stationary random process in the wide sense with finite variance;

• (C6)

  $T(t)$ is a smooth function so that its Fourier transform, $\widehat{T}$, is compactly supported in $[-\Delta,\Delta]$.

The model’s identifiability is discussed in [8]. A typical example adhering to (1) is the accelerometer signal, with $L=1$, shown in Fig. 2. Clearly, the IMT function does not oscillate sinusoidally. This model assumes one WSF for each IMT function. In essence, the $\ell$-th IMT function arises from stretching $s_{\ell}$ by $\phi_{\ell}$, subsequently scaled by $a_{\ell}$. Nevertheless, in practical instances, researchers have observed WSFs not to be fixed [2, 14, 39]. See Fig. 2 for an example. It is evident that the oscillatory pattern varies over time. Notably, both STFT and SST reveal that the strength of the fundamental component (indicated by the red arrow) weakens, while the third, sixth, and seventh harmonics (indicated by purple arrows) strengthen after 225 seconds. As a result, (1) lacks appropriateness as a model.

To properly model such signals, we consider the generalization of (1) in [21]. To motivate this generalization, assume $L=1$ to simplify the discussion. We have

where $\{\alpha_{j}\}\subset[0,\infty)$ and $\{\beta_{j}\}\subset[0,2\pi)$ are determined by Fourier coefficients of $s_{1}$. The generalization idea in [21] involves permitting time-varying characteristics for each harmonic in (2), given specific conditions. Thus, consider

where $D_{\ell}\in\mathbb{N}$ is called the harmonic order, and

1. (C7)

   $\phi_{\ell,1}$ behaves like $\phi_{\ell}$ in (C1)-(C2), $\phi_{\ell,j}\in C^{2}$, and $|j-\phi^{\prime}_{\ell,j}(t)/\phi^{\prime}_{\ell,1}(t)|\leq\epsilon^{\prime}$ for $j=2,\ldots$ and a small constant $\epsilon^{\prime}\geq 0$ for all $t\in\mathbb{R}$;

2. (C8)

   $a_{\ell,1}(t)$ behaves like $a_{\ell}$ in (C3), $a_{\ell,j}\in C^{1}$ and $|a^{\prime}_{\ell,l}(t)|\leq\epsilon^{\prime}\phi^{\prime}_{\ell,1}(t)$ for all $t\in\mathbb{R}$.

Conditions (C7) and (C8) capture the time-varying WSF, assuming slow variation. This model is a simplified version of that in [21] to facilitate discussion. For $j\geq 1$, $a_{\ell,j}(t)\cos(2\pi\phi_{\ell,j}(t))$ is called the $j$-th harmonic of the $\ell$-th IMT function, and $a_{\ell,1}(t)\cos(2\pi\phi_{\ell,1}(t))$ is the fundamental component. When $D_{\ell}=1$, the $\ell$-th IMT function oscillates sinusoidally. In this work, we assume the knowledge of $L$ and apply [31] to estimate $D_{\ell}$ to avoid the distraction. Although there are tools when $D_{\ell}=1$ for all $\ell$ [37, 32, 19, 31], estimating $L$ in general is challenging. The slowly varying IF condition (C2) can be extended to a fast varying one [18], which requires more technical details. To stay focused, we will discuss this generalization in future work.

Due to the non-stationary nature of ANHM, conventional time series methods may be limited. A common strategy to analyze such signals involves employing TF analysis algorithms [15]. The core notion behind TF analysis is a “divide-and-conquer” approach. Locally assuming a “stable”, “time-invariant”, or “stationary” structure, we estimate this structure accurately in patches. By combining local estimates, we construct a TFR, which enables solving signal processing tasks like AM and IF estimation of each IMT function, signal decomposition into essential components (IMT functions and their harmonics), denoising, and dynamic feature extraction. This paper concentrates on STFT and SST, and we refer readers to [15] for other TF analysis algorithms.

Below, we explain the overall idea. Mathematically, for a function $f$ satisfying mild conditions (e.g., a tempered distribution), the STFT is defined as

where $t\in\mathbb{R}$ is the time, $\xi\in\mathbb{R}$ is the frequency, and $h$ is a chosen window that is smooth and decays sufficiently fast (e.g., the Gaussian function). The spectrogram is $|V^{(h)}_{f}(t,\xi)|^{2}$, yielding a function on $\mathbb{R}^{2}$. When $L>1$, the essential support of $\hat{h}$ should reside in $[-\Delta/2,\Delta/2]$ or the spectral interference of different IMT functions might happen. Generally, any TF analysis method’s output is a function on $\mathbb{R}\times\mathbb{R}^{+}$, known as the TFR. We call $\mathbb{R}\times\mathbb{R}^{+}$ the TF domain. See Fig. 1(b) for an example, where the input signal follows ANHM with $L=1$, where the fundamental component is weaker compared to its second harmonic and diminishes over time, as evident in the spectrogram. Nonetheless, STFT encounters the uncertainty principle [30]; that is, the TFR appears ’blurred’ due to the chosen window. To address this, SST [12], a reassignment algorithm [3] variant, and its generalizations like second-order SST [27] were introduced. They utilize the phase information within the complex-valued TFR obtained through STFT to sharpen the STFT-based TFR. Let $S_{f}^{(h,[1])}:\mathbb{R}\times\mathbb{R}^{+}\to\mathbb{C}$ and $S_{f}^{(h,[2])}:\mathbb{R}\times\mathbb{R}^{+}\to\mathbb{C}$ denote the original SST [12] and second-order SST [27] TFRs of function $f$, respectively.

To proceed, the conventional procedure entails extracting ridges connected to harmonic IFs. In contrast to STFT, a sharper TFR, as determined by SST, can enhance ridge estimation accuracy [24]. It has been theoretically established that when a signal adheres to ANHM with sinusoidal WSFs, the ridges closely approximate IFs of harmonics of all IMT functions [13, 12]. Once ridges are accurately extracted, IMTs can be reconstructed through a reconstruction formula. Specifically, suppose the extracted ridge $c(t)$ approximates IF of $a_{\ell,j}(t)\cos(2\pi\phi_{\ell,j}(t))$ in (3) for some $\ell$ and $j$. Then, if $\phi_{\ell,j}^{\prime}(t)$ does not intersect with and is away from IF of any other harmonics, the following formula holds:

where $q=1,2$ and user-defined bandwidth $\Delta_{q}>0$. This formula yields estimates of $a_{\ell,j}(t)$ through absolute value extraction and $\phi_{\ell,j}(t)$ through phase unwrapping [1]. With these estimates, further signal processing tasks become feasible. Regrettably, the condition that $\phi_{\ell,j}^{\prime}(t)$ must not intersect with or approach the IF of other IMT function’s harmonics, as required by (5), is excessively stringent and practically infeasible. Thus, an alternative approach becomes necessary. Various mode retrieval algorithms for STFT exist, such as ridge-based skeleton [6], integration approach [19], basin attractor utilization [22, 23, 19], penalization with detected ridges [7], and others. However, these methods face similar limitations posed by non-sinusoidal WSFs.

See Section A for a discussion of the numerical implementation of STFT and SST. Below, we denote the discretized TFR of $f$, either determined by STFT or SST, as $\mathbf{R}\in\mathbb{C}^{N\times M}$, where $N$ is the number of the sampling points of the signal with the sampling period $\Delta_{t}$, and $M$ is the number of the frequency domain ticks, each bin spaced by $\frac{1}{2M\Delta_{t}}$.

The key step in the signal processing process outlined in Section 2.2 involves an accurate RD algorithm. RD algorithms can be broadly categorized into three groups. One involves detecting ridges individually and gathering all ridges through an iterative peeling scheme [6, 8, 10]. The second is the multiridge scheme, which detects multiple ridges concurrently [7]. The third is the preprocessing scheme, where the TF domain is segmented to include one ridge per portion or the TFR is enhanced before applying either peeling or multiridge algorithms [19]. Given that our proposed algorithm adheres to the peeling scheme and the focus is handling WSFs, our focus remains on reviewing algorithms within this category.

In the peeling scheme, once we obtain a ridge $c_{1}:[N]\rightarrow[M]$, where $[N]:=\{1,\cdots,N\}$, from $\mathbf{R}_{1}:=\mathbf{R}$, we set iteratively the following masked TFR, $\ell=2,\ldots,L+1$:

where $B_{\ell-1,n}:=[c_{\ell-1}(n)-\eta_{-}(n),c_{\ell-1}(n)+\eta_{+}(n)]$, $c_{\ell-1}$ is the ridge extracted from $\mathbf{R}_{\ell-1}$, and $\eta_{-}(n)\geq 0$ and $\eta_{+}(n)\geq 0$ is the bandwidth selected by the user. Note that we assume the knowledge of $L$ in this work. Practically, the selection of $\eta_{-}(n)$ and $\eta_{+}(n)$ depends on the TFR and profoundly impacts RD performance. Several methods exist for masking the TFR. For instance, one can opt for a constant frequency bandwidth (CFB) approach, where $\eta_{-}(n)$ and $\eta_{+}(n)$ remain constant for all $n$ [10]. Alternatively, the bandwidth can adapt to noise levels, known as the varying frequency bandwidth (VFB) [10]. Another option involves the modulation-based (MB) approach, which relies on estimating spectral spreading due to chirp [10].

Numerous algorithms exist for fitting a ridge to the (masked) TFR $\mathbf{R}_{\ell}$ in each iteration. These algorithms share the common goal of approximating a curve within the TFR to capture maximum energy while meeting specific conditions. A frequently used method is to solve the following path optimization problem [6, 8]:

where $\Delta c\in[M]^{N-1}:=\{(m_{1},\cdots,m_{N-1}):m_{k}\in[M],\,1\leq k\leq N-1\}$ so that $\Delta c(j):=c(j+1)-c(j)$ for $j=1,\ldots,N-1$, which is the numerical differentiation of the curve $c$, $\lambda_{1}>0$ is the penalty term constraining the regularity of the fit curve $c^{*}$, and $\widetilde{\mathbf{R}}(\ell,q)=\log\frac{|\mathbf{R}(\ell,q)|}{\sum_{i=1}^{N}\sum_{j=1}^{M}|\mathbf{R}(i,j)|}$ is a normalization of the matrix $\mathbf{R}$. In practice, various approaches to solving the optimization problem (9) exist, including the simulated annealing algorithm [6]. However, we recommend the more efficient penalized forward-backward greedy algorithm, referred to as the single curve extraction algorithm (Single-RD) [6, 8]. Note that we could further enhance regularity by considering second-order numerical differentiation in (9), or even incorporating available noise information [6]. For simplicity, we omit these aspects from the current discussion.

Several other algorithms proposed in [17, 10, 19] could be viewed as solving variations of (9). For example, the frequency modulation (FM) approach proposed in [10] solves $c^{*}=\operatorname*{argmax}_{c:[N]\rightarrow[M]}\sum_{j=1}^{N}\left|\mathbf{R}(j,c(j))\right|^{2}$ such that $|\Delta c(j)-\frac{K}{N^{2}}\hat{q}(j,c(j))|\leq C$ for some $C\geq 0$, where $\hat{q}$ is the chirp rate estimate [10, eq. (6)]. This optimization problem could be solved by, for example, the recurring principle [19, eq. (11)], leading to the FM-RD algorithm. In [19], to handle substantial noise, the concept of a relevant ridge portion (RRP) is introduced, giving rise to the RRP-RD algorithm. RRPs correspond to the desired ridges and are employed to segment the TF domain into specific structures known as basins of attraction [23]. The spline least-square approximation is applied to RRPs for RD.

As an example, consider the simplest scenario with only one IMT function ($L=1$ in (3)) that oscillates sinusoidally ($a_{1,j}=0$ for all $j\geq 2$). In this case, the ridge $c^{*}$ determined by any of these RD algorithms provides an accurate estimate of the associated fundamental component’s IF [12], followed by other signal processing processes. Nonetheless, this simplest scenario is rarely encountered in practice.

To illustrate the challenges posed by existing RD algorithms, consider scenarios closer to real-world situations. First, when $L>1$ in (3) and all IMT functions oscillate sinusoidally, we can apply any of the aforementioned RD algorithms to estimate the IFs associated with all IMT functions, assuming the IF separation condition in (C2) holds. Second, when $L=1$ in (3) but the WSF significantly deviates from being sinusoidal, we can apply the above RD algorithms based on (C7) to estimate all ridges linked to the harmonics. The ridge associated with the lowest IF can then be assumed as the fundamental component’s IF. In these more realistic scenarios, successfully applying existing RD algorithms becomes more challenging. An example in Fig. 1(f) demonstrates this complexity, where the strength of the fundamental component is weak in the middle while the strength of the second harmonic dominates. Such shifts in strengths can easily confuse existing RD algorithms.

In real-world scenarios, more intricate situations often challenge the aforementioned RD algorithms. When $L>1$ in (3) and the WSFs deviate significantly from sinusoidal patterns, the TFR can become complex with multiple curves representing different harmonics of distinct IMT functions. This complexity involves interference among harmonics from different IMT functions and can lead to mode mixing. For instance, as shown in Fig. 3, two IMT functions exhibit non-sinusoidal WSFs. Here, the second harmonic of the first IMT function (red arrows) intersects with the fundamental component of the second IMT function (blue arrows), making it challenging to accurately determine the IF of the fundamental component for each IMT function. Given that numerous signal processing tasks, such as SAMD [11], rely on accurately acquiring the phase of each IMT function and considering the prevalence of such signals, particularly in biomedical applications, there is a compelling need for a new RD algorithm.

Given an input signal $f_{1}=f$ that satisfies (3) with $L\geq 1$, let $\mathbf{R}_{1}\in\mathbb{C}^{N\times M}$ represent the discretized TFR of $f_{1}$, obtained from STFT or SST. Choose $K\geq 1$ as the intended number of harmonics to extract simultaneously. Using the notations from (9), we fit $K$ ridges for the harmonics of an IMT function by solving the optimization problem:

where $e_{k}\in\mathbb{R}^{K}$, $k=1,\ldots,K$, is a unit vector with $e_{k}(k)=1$, $\lambda_{k}\geq 0$ is the smoothness penalty as in (9), and $\mu_{k}\geq 0$ stands for the similarity penalty that captures the time-varying wave-shape condition stated in (C7). Note that $\mu_{k}>0$ forces the estimated IF of the $k$-th harmonic, $e_{k}^{\top}\mathbf{c}$ to satisfy $e_{k}^{\top}\mathbf{c}(\ell)\approx ke_{1}^{\top}\mathbf{c}(\ell)$ at time $\ell\Delta_{t}$, and hence the geometric structure is respected. The output $\mathbf{c}^{*}\in[M]^{K\times N}$ represents the IFs of the first $K$ harmonics of an IMT function of the input signal $f$.

We simplify (10) by introducing a set

where $\beta\geq 0$ is the chosen “bandwidth” that echos the similarity penalty $\mu_{k}$ in (10), so that (10) becomes

Solving (12) directly can be inefficient. We recommend a segmentwise approach. This method leverages the slowly-varying IF’s uniform continuity by approximating the curve with a finite sequence of “thin rectangles” $\{[t_{q-1},\,t_{q}]\times[L_{q},U_{q}]\}_{q=1}^{Q}$. First, estimate the fundamental IF $\hat{c}$ by running the algorithm on $[t_{0},\,t_{1}]$ over the whole frequency domain. Then, for each segment $[t_{q},\,t_{q+1}]$, run the algorithm over a frequency band $[-B_{q},B_{q}]+\hat{c}(t_{q})$, where $B_{q}>0$ is chosen properly to ensure the box $[t_{q},\,t_{q+1}]\times[\hat{c}(t_{q})-B_{q},\hat{c}(t_{q})+B_{q}]$ covers the ridge. This process determine the fundamental IF over $[t_{q},t_{q+1}]$. Note that a smaller $B_{q}$ reduces the search space and improves efficiency, but it must not be too small to avoid missing the ridge. We suggest using $B_{q}=1$ and $t_{q}-t_{q-1}=1$ for all $q$, based on the slowly-varying IF assumption to ensure the ridge is covered. This section of the algorithm is referred to as MHRD and is summarized in Algorithm 2.

After identifying one IMT function’s harmonic ridges using MHRD, we employ the SAMD algorithm [11] to extract the corresponding IMT function ($\mathbf{x}^{[1]}_{1}$ in Algorithm 1) related to the detected ridges with its fundamental component’s IF. In short, SAMD is an optimization-based fitting algorithm that uses the estimated fundamental component of an IMT function ($\mathbf{x}_{\ell,1}^{[i]}$ in Algorithm 1 using (5) with $e_{1}^{\top}\mathbf{c}_{(\ell)}$) as inputs. The loss function accounts for the time-varying WSF, and its minimizer provides an estimate of the IMT function. See [11] for details.

Subtracting this extracted IMT function from the input signal $f_{1}$ results in a new signal denoted as $f_{2}$, containing $(L-1)$ IMT functions. The discretized TFR of $f_{2}$, determined by STFT or SST, is represented as $\mathbf{R}_{2}\in\mathbb{C}^{N\times M}$. Note that ridges related to the extracted IMT function do not exist in $\mathbf{R}_{2}$. This step thus serves as a replacement for (8).

Having obtained the initial IMT function and $\mathbf{R}_{2}$, we can proceed to iterate Step 1 (MHRD) on $f_{2}$ and $\mathbf{R}_{2}$ to extract the harmonic ridges of one of the remaining IMT functions. Then, Step 2 is once again executed to extract the associated IMT function using the SAMD algorithm. This iterative process is repeated for a total of $L$ cycles, resulting in the complete set of harmonic ridges for all IMT functions. These initial ridge estimates, denoted as $c_{1}^{*},\ldots,c_{L}^{*}$, are ordered based on the fundamental components’ IFs, from low to high.

In practice, it is not guaranteed which IMT function will be extracted in Steps 1 and 2. It is possible that the IF and phase of its associated fundamental component is affected by the harmonics of other IMT functions with lower IF. Thus, with the initial estimate $c_{1}^{*},\ldots,c_{L}^{*}$, we repeat Steps 1 and 2 on $f_{1}$ and $\mathbf{R}_{1}$ by replacing $\mathbf{S}^{(K)}(\beta)$ in (12) by

and estimating the phase via (5). This ridge estimate is less influenced by harmonics from other IMT functions due to the IF separation condition (C2). We then repeat Steps 1 and 2 on $f_{2}$ and $\mathbf{R}_{2}$ with $\mathbf{S}^{(K)}$ modified in the same way as (13) by considering $c_{2}^{*}$ . We continue this process iteratively with $c_{3}^{*},c_{4}^{*},\ldots$ and so on until ridges of all $L$ IMT functions are extracted, denoted as $\bar{c}_{1}^{*},\ldots,\bar{c}_{L}^{*}$.

While it is valid to conclude the process after these iterations, it is possible to enhance the results through an additional iteration of the $L$ cycles with $\bar{c}_{1}^{*},\ldots,\bar{c}_{L}^{*}$. This reiteration step has the potential to yield improved outcomes, and we recommend considering this option for $I\geq 1$ times. It is worth noting that the output includes the decomposition of all IMT functions from the input signal $f$.

Solving the optimization problem (10) can be time-consuming, especially when $K$ is large and $\beta$ is wide. However, if the phase of the fundamental component of the targeted IMT function is known or can be well estimated in advance, the algorithm can be accelerated. Let $\widetilde{\mathbf{R}}$, $\lambda_{k}$ and $\mu_{k}$ be the same as (10), and $c_{1}:[N]\rightarrow[M]$ be the known IF of the fundamental component; that is, the extra auxiliary information. In this case, we only need to estimate the higher-order harmonics, and hence (10) is reduced to

where $\mathcal{L}_{c_{1}}^{(k)}(\mathbf{f}):=\sum_{j=1}^{N}\left|\tilde{\mathbf{R}}\left(j,\mathbf{f}(j)\right)\right|-\lambda_{k}\sum_{j=2}^{N}|\mathbf{f}(j)-\mathbf{f}(j-1)|^{2}-\mu_{k}\sum_{\ell=1}^{N}|\mathbf{f}(j)-kc_{1}(j)|^{2}$ for $k\in\{2,\cdots,K\}$ and $\mathbf{f}:[N]\rightarrow[M]$. The following proposition shows that the optimization problem (14) can be reduced to several single-valued curve optimization problems, where the proof can be found in Section C.

This proposition allows us to estimate harmonics’ IFs by solving a modified single curve extraction program with respect to $c_{1}$ for each $2\leq k\leq K$; that is,

instead of solving (14) directly. This reduces the complexity of the original problem. This auxiliary information is often available in many applications. For example, when analyzing photoplethysmogram (PPG) and electrocardiogram (ECG) signals, we can use the cardiac phase from the ECG as an accurate phase estimate for the cardiac component in the PPG signal. This allows us to speed up SAMD-MHRD when analyzing the cardiac component in the PPG signal.

To optimize the performance of SAMD-MHRD for a given $K$, selecting suitable penalties $\lambda_{1},\ldots,\lambda_{K}$ is crucial. To address this, we propose a simple grid search approach with the following rationale. Assuming that the algorithm accurately detects the ridges corresponding to the true IFs, we can consider masking the extracted curves on the TFR, similar to (8), by appropriately choosing $\eta_{-}$ and $\eta_{+}$. To simplify the discussion, we apply the CFB approach [10] below, while VFB and MB can also be considered. The resulting TFR should exhibit reduced signal components. This insight prompts us to employ the $\alpha$-Rényi entropy [4, 33], where $\alpha>0$, as a metric to quantify and compare the energy distribution along the frequency axis of both the masked and original TFRs.

Consider the parameter set

where $\delta_{\lambda}>0$ is small, $K$ satisfies $1-(K-1)\delta_{\lambda}>0$, and $\Lambda\subset(0,\infty)$ is a finite set. Note that while $\mathcal{S}_{\delta_{\lambda},\Lambda}$ looks like a high-dimensional space, it is effectively a one-dimensional space. In this context, the ordering $\lambda_{k}<\lambda_{\ell}$ for $k>\ell$ adheres to the slowly varying condition $|\phi_{j}^{\prime\prime}|\leq\epsilon j\phi_{1}^{\prime}$ for $j\geq 1$. For the parameter $\beta$ required to satisfy the WSF condition in $\mathbf{S}^{(K)}(\beta)$ (11), we consider a finite subset $M\subset(0,2^{-1})$.

Take the TFR $\mathbf{R}$. For each $\lambda_{1}\in\Lambda$ and $\beta\in M$, extract the curves by SAMD-MHRD with $(\lambda_{1},\cdots,\lambda_{K})\in S_{\delta_{\lambda},\Lambda}$ and $\beta$, and denote the result to be $\mathbf{c}^{*}(\lambda_{1},\beta)$. Then, “mask” $\mathbf{R}$ by the ridges $\mathbf{c}^{*}(\lambda_{1},\beta)$ following (8), and denote the masked TFR as $\mathbf{R}_{(\lambda_{1},\beta),\Delta}$, where the masking size $\Delta>0$ respects the one used in the reconstruction formula (5). Define a sequence $(q^{(\lambda_{1},\beta)}(1),\cdots,q^{(\lambda_{1},\beta)}(N))$ by $q^{(\lambda_{1},\beta)}(\ell):=R_{\alpha}(|\mathbf{R}_{(\lambda_{1},\beta),\eta}(\ell,\cdot)|)$, $\ell\in[N]$. Finally, we choose

as the optimal $(\lambda_{1},\beta)$ in $\Lambda\times M$. In practice, $\Lambda$ is chosen to range from $10^{-1}$ to $10^{1}$ and $M$ is chosen to range from $2^{-6}$ to $2^{-1}$, both with a uniform grid in the log scale.

In this first example, we compare different RD algorithms on signals with complicated WSF dynamics. Denote the smoothed standard Brownian motion $X_{1}(t)=W\ast K_{B}(t)$, where $W$ is the standard Brownian motion, $\ast$ indicates convolution, and $K_{B}$ is the Gaussian function with the standard deviation $B=20$ seconds. Over $[0,50]$, set

and

where $X_{2}$ is an independent and identical copy of $X_{1}$, $U_{\ell}$ are independent white Gaussian process with mean 0 and variance $0.1$, $B^{\prime}=5$ seconds, and

Note that the term $U_{\ell}\ast K_{B^{\prime}}(t)$ models the time-varying WSF. Consider a non-stationary noise $\Phi$, where for $n=1,\ldots,5,000$, $\Phi(n)$ follows an ARMA$(1,1)$ process with Student’s t4 innovation process, and for $n=5,000,\ldots,10,000$, $\Phi(n)$ i.i.d. follows a Student’s t5 distribution. We assume $\Phi$ is independent of $X_{1}$ and $X_{2}$ The simulated signal is defined over $[0,50]$ with the sampling rate $200$Hz so that

where $n=1,\ldots,10,000$, $s_{1}(t):=A_{1}(t)x_{1}(t)$, $\sigma\geq 0$ is determined by the desired signal-to-noise ratio (SNR), defined as $20\log_{10}\frac{\texttt{STD}(s_{1})}{\texttt{STD}(\Phi)}$ with STD standing for the standard deviation, $x_{1}(t)=D_{1}\cos(2\pi\phi_{1}(t))+(u_{1}+u_{2})\cos(2\pi\phi_{2}(t))+u_{1}\cos(2\pi\phi_{3}(t))$, $D_{1}\in(0,1]$ models the intensity of the fundamental component, and $u_{1},u_{2}$ are independent random variables uniformly distributed on $[0,1)$ that are independent of $X(t)$ and $\Phi(n)$. See Figure 8 in Section B for a realization of $Y_{1}(t)$.

Then, we realize independently $Y_{1}(t)$ over $[0,50]$ for $100$ times. For a fair comparison with existing algorithms, we utilize the provided code from the authors [10, 19], following their recommended guidelines.³³3The code for FM-RD is obtained through private communication with the authors. The code for RRP-RD is available at https://github.com/Nils-Laurent/RRP-RD. However, since these algorithms are not tailored to handle WSFs, we need to make minor adjustments and assume the number of harmonics is known. In the case of FM-RD, we incorporate the peeling step (8) as a post-processing procedure. Essentially, we run FM-RD iteratively and apply TFR masking three times. Among the three extracted curves, we identify the one with the lowest frequency value as the result for the fundamental component’s ridge. Regarding RRP-RD, we concurrently detect three relevant ridge portions on the TF domain and extract the corresponding ridges. From these, we choose the ridge associated with the lowest frequency value as the outcome. If RRP-RD is unsuccessful in identifying three ridge portions, we rerun it to detect two RRPs. If this attempt also fails, we proceed with a single RRP. Finally, we assign the ridge with the lowest frequency as the final outcome.

Let $c_{1}^{(\square)}$ be the detected ridge by the algorithm $\square$, which could be S, FM, RRP or MH, standing for the Single-RD, FM-RD, RRP-RD or SAMD-MHRD. To compare different algorithms, we consider the following metric:

We run the Wilcoxon signed-rank test on $100$ realizations to test the null hypothesis that the median of $\delta^{(\square)}-\delta^{(\textsf{MH})}$ is zero, with a significance level of $0.05$ and Bonferroni correction.

The results are depicted in Fig. 4, illustrating the distributions of ${\delta^{(\square)}}$ for a fixed $D_{1}$ (rows: $D_{1}=0.1$, $D_{1}=0.2$, $D_{1}=0.5$) and a constant noise level (columns: noise-free, SNR = 5 dB, SNR = 0 dB) with asterisks indicating instances where $p<0.05/27$. The outcomes reveal that in the absence of noise, all algorithms yield comparable results. However, when encountering scenarios with a weak fundamental component and noise interference, FM-RD and RRP-RD exhibit limitations compared to SAMD-MHRD. More numerical results for RRP-RD with different number of ridge portions can be found in Section B.

Following the same notations used in Section 4.1, we consider the case when $L=2$. Denote $X_{3}$ and $X_{4}$ as two independent copies of the smoothed standard Brownian motion and are independent of $X_{1}$, $X_{2}$, $\Phi$, $u_{1}$ and $u_{2}$. Over $[0,50]$, set two random processes

and

where $\phi_{0}^{*}(t)=2.33t+0.2t^{2}$. The second simulated signal is defined over $[0,50]$ with the sampling rate $200$Hz so that

where $s_{2}(t):=A_{2}(t)x_{2}(t)$, $x_{2}(t)=D_{2}\cos(2\pi\phi_{1}^{*}(t))+(u_{1}^{*}+u_{2}^{*})\cos(2\pi\phi_{2}^{*}(t))+(u_{1}^{*})\cos(2\pi\phi_{3}^{*}(t))\big{]}$, $D_{2}\in(0,1]$ models the intensity of the fundamental component, $u_{1}^{*}$ and $u_{2}^{*}$ are the independent copies of $u_{1}$ and $\sigma\geq 0$ is again determined by the desired SNR. Then we realize 100 copies of $Y_{2}$ for each selected SNR.

See Fig. 5 for the RD results and a decomposition example of SAMD-MHRD. Both IMT functions are successfully recovered. We then run SAMD-MHRD and evaluate the reconstruction error of the $j$-th IMT function as

where the superscript $q$ indicates the $q$-th iteration. See Fig. 6 for the comparison of $\delta^{[1]}_{j}$ and $\delta^{[3]}_{j}$ under different SNRs. As SNR decreases, errors increase. Additionally, the second IMT exhibits significantly larger recovery error than the first IMT, as tested with the Wilcoxon test at a significance level of $0.005$, which is reasonable due to its dependence on the harmonics of the first IMT. Iteration improvement is statistically significant except for the second IMT at high noise levels.