Nonlinear mode decomposition: A noise-robust, adaptive decomposition method

The signals emanating from complex systems are usually composed of a mixture of different oscillations which, for a reliable analysis, should be separated from each other and from the inevitable background of noise. Here we introduce an adaptive decomposition tool—nonlinear mode decomposition (NMD)—which decomposes a given signal into a set of physically meaningful oscillations for any wave form, simultaneously removing the noise. NMD is based on the powerful combination of time-frequency analysis techniques—which, together with the adaptive choice of their parameters, make it extremely noise robust—and surrogate data tests used to identify interdependent oscillations and to distinguish deterministic from random activity. We illustrate the application of NMD to both simulated and real signals and demonstrate its qualitative and quantitative superiority over other approaches, such as (ensemble) empirical mode decomposition, Karhunen-Loève expansion, and independent component analysis. We point out that NMD is likely to be applicable and useful in many different areas of research, such as geophysics, finance, and the life sciences. The necessary matlab codes for running NMD are freely available for download.

Figure 1

(a) The central 50-s section of the signal $s\left(t\right)$ specified at the top, which is composed of three oscillations corrupted with white Gaussian noise of standard deviation equal to $0.5$. The oscillations are (1) near 10 Hz, with frequency modulation (chirp); (2) near 6 Hz, with amplitude modulation; (3) near 1 Hz, with both modulations and a nonsinusoidal wave form. (b),(c) Respectively, the central parts of the WFT and WT of the signal $s\left(t\right)$ shown in (a), where the ridge curves related to each oscillation are shown in each case by solid lines of the corresponding colors: 1, magenta (top line); 2, - green (second from the top); 3, brown (two bottom lines); the third oscillation is represented by two curves because of its complex wave form, while the apparent distortions in the ridge frequency profiles are due to noise; note also the logarithmic frequency scale for the WT, which is natural for it. The signal was sampled at 100 Hz for 100 s.

Figure 2

Illustration of the extraction of the $h\mathrm{th}$ harmonic ridge curve ${\omega }_p^{\left(h\right)}\left(t\right)$ based on the fundamental ridge frequency ${\omega }_p^{\left(1\right)}\left(t\right)$. At each time, starting from the expected ridge frequency $h{\omega }_p^{\left(1\right)}\left(t\right)$ of the harmonic (blue diamonds), one climbs (i.e., follows in the direction of TFR amplitude increase) to the nearest peak, which is then assigned to ${\omega }_p^{\left(h\right)}\left(t\right)$ (red circles).

Figure 3

(a) Nonlinear mode with amplitude modulation, as specified by the equation for $s\left(t\right)$ above the box. (b),(c) The corresponding WFT and WT amplitudes, respectively. (d),(e) The signal's FT, integrated over one frequency bin (vertical lines), with the parts responsible for different harmonics shown using different colors and markers; the shaded areas show the absolute values of the window functions $\widehat{g}(\omega -\xi )$ or wavelet functions ${\widehat{\psi }}^{*}({\omega }_{\psi }\xi /\omega )$ centered at the mean frequencies of the harmonics $\omega =2\pi h$ and rescaled to half of the harmonics' mean amplitudes. The signal was sampled at 50 Hz for 50 s.

Figure 4

(a) Nonlinear mode with frequency modulation, as specified by the equation for $s\left(t\right)$ above the box. (b),(c) The corresponding WFT and WT amplitudes, respectively. (d),(e) The signal's FT, integrated over one frequency bin (vertical lines), with the parts responsible for different harmonics shown using different colors and markers; the shaded areas show the absolute values of the window functions $\widehat{g}(\omega -\xi )$ or wavelet functions ${\widehat{\psi }}^{*}({\omega }_{\psi }\xi /\omega )$ centered at the mean frequencies of the harmonics $\omega =2\pi h$ and rescaled to half of the harmonics' mean amplitudes. The signal was sampled at 50 Hz for 55 s.

Figure 5

(a) The central 50-s section of the signal $s\left(t\right)$ specified by the equation at the top, which represents a single NM corrupted by Brownian noise of standard deviation equal to 4 [Brownian noise of unit deviation ${\eta }_B\left(t_n\right)\sim {\sum }_{m=1}^n{\eta }_W\left(t_{m-1}\right)\sqrt{t_m-t_{m-1}}$ is obtained as the normalized cumulative sum of the Gaussian white noise signal ${\eta }_W\left(t\right)$]. (b),(c) Central parts of the WFT and WT, respectively, of the signal $s\left(t\right)$ shown in (a). The signal was sampled at 100 Hz for 100 s.

Figure 6

The WFTs from which each true harmonic is reconstructed, with the corresponding extracted ridge curves being shown by solid red lines. The window resolution parameter $f_0$ is adjusted individually for each harmonic. After the harmonic is reconstructed, it is subtracted from the signal, so that, e.g., the first harmonic no longer appears in the WFTs of (b)–(d). Note that the color scaling in (b)–(d) covers half the amplitude range of that in (a).

Figure 7

The result of applying NMD to the signal shown in Fig. 5. (a) The reconstructed NM (black line) is compared with the true NM (gray background line); (b) the residual provided by NMD (black line) is compared with the actual background noise (gray background line).

Figure 8

The result of applying EMD (left panels) and EEMD (right panels) to the signal shown in Fig. 5. Thin red lines, where present, show the real first harmonic (in C4 for EMD and C5 for EEMD), sum of the third and fourth harmonics (in C3 for EMD and C4 for EEMD), and the seventh harmonic (in C2 for EMD and C3 for EEMD). The bottom panels show the sum of components 7 to 13. For EEMD, we have used 1000 noise realizations with standard deviations ten times smaller than that of the signal.

Figure 9

(a) The central 20-s section of the signal $s\left(t\right)$ specified by the equation at the top, representing the sum of two NMs corrupted by Gaussian white noise of standard deviation equal to 1.725; the phases of the NMs ${\varphi }_{1,2}\left(t\right)$ were obtained as ${\varphi }_{1,2}\left(t\right)={\int }_0^t{\nu }_{1,2}\left(\tau \right)d\tau$, with ${\nu }_{1,2}\left(t\right)/2\pi =(1,2)+0.01{\tilde{\eta }}_{B;1,2}\left(t\right)$, where ${\tilde{\eta }}_{B;1,2}\left(t\right)$ are two independent realizations of the unit deviation Brownian noise filtered in the range $[0.01,0.2]$ Hz. (b),(c) Central parts of the WFT and WT of the signal $s\left(t\right)$ shown in (a). The signal was sampled at 100 Hz for 100 s.

Figure 10

The result of applying NMD to the signal shown in Fig. 9. In (a) and (b) black lines indicate the two reconstructed NMs and the gray background lines show the true NMs for comparison. (c) Similarly, the black and gray lines show the residual returned by NMD and true background noise, respectively.

Figure 11

(a) An example of the blood flow signal, measured by LDF with the probe over the right wrist caput ulna (for more details, see [3, 53]); the signal is sampled at 40 Hz for 1800 s, and the panel shows its central 600 s part. (b) The WT of the signal shown in (a). Gray dotted lines partition the frequency axis into the regions within which the physiologically meaningful oscillations are located (according to [3, 52, 53, 54]; see the text). Bold colored lines show those extracted components which pass the surrogate test against noise, with thinner lines of the same color showing their higher harmonics. (c),(d) The reconstructed NMs corresponding to four upper magenta lines and two bottom brown lines on (b), respectively; the main graphs show the modes during 50 s and small insets during 600 s [as in (a)]. (e),(f) The wave forms of the oscillations shown in (c) and (d), respectively; the gray dashed line shows a pure sinusoidal wave form of amplitude equal to that of the first harmonic and is provided as a guide for the eye. The cardiac wave form (e) has four harmonics $h=1,2,3,4$, with $a_h=[1,0.52,0.37,0.16]$ and ${\phi }_h/\pi =[0,0.4,0.66,-0.97]$ [in the notation of (4)]; the wave form in (f) is characterized by $h=1,2, a_h=[1,0.08]$, and ${\phi }_h/\pi =[0,0.07]$.

Figure 12

Same as in Fig. 11, but for the blood flow measured from a different subject. The nonlinear modes shown in (c) and (d) correspond, respectively, to the four upper lines (magenta) and the bottom (green) line in (b). The cardiac wave form (e) has four harmonics $h=1,2,3,4$, with $a_h=[1,0.32,0.15,0.04]$ and ${\phi }_h/\pi =[0,0.41,0.95,-0.44]$ [in the notation of (4)]; the respiratory oscillations have only one harmonic, so that the wave form in (f) is a pure sinusoid.

Figure 13

(a) The EEG signal [measured using a ${\mathrm{BIS}}^{\mathrm{TM}}$ 4-Electrode Sensor (Covidien-Medtronic, Dublin, Ireland and Fridley, Minnesota) placed on the forehead and a purpose-built signal conditioning unit (Jožef Stefan Institute, Slovenia), as described in 59]. (b) The ECG signal (three-lead, with electrodes on shoulders and the lowest left rib, see, e.g., [60]). (c),(d) The WFTs of the EEG and ECG signals shown in (a) and (b), respectively, calculated using the default resolution parameter $f_0=1$. (e) The WFT of the EEG signal calculated using the adaptively adjusted resolution parameter $f_0=3.76$, due to which the cardiac component becomes much more visible than in (c). (f) The cardiac artifacts extracted from the EEG, with an inset showing them over the same 500 s as presented in (a),(b); gray background lines show the ECG scaled to the dimensions of the plot. (g) The wave form of the cardiac artifacts, which has four harmonics with $a_h=[1,0.33,0.24,0.11]$ and ${\phi }_h/\pi =[0,0.49,1,-0.59]$ [in the notation of (4)]; the gray dashed line represents a pure sinusoid with the amplitude of the first harmonic, for comparison. The full signals are of duration 20 min and are sampled at 80 Hz; the EEG was actually measured under anaesthesia, but the same artifacts arise under all conditions.