Tutorial

The example described below is provided in emd_example.py.

First, import all nedeed libraries:

import numpy as np
import colorednoise as cn
import matplotlib.pyplot as plt

from scope.fourier import fit_fourier
from scope.emd import emd_modes, emd_trend, emd_energy_spectrum, emd_noise_conf
from scope.utils import plot_modes, plot_signal, plot_fft_spectrum, plot_emd_spectrum

The sample signal in this example consists of an oscillatory component, an exponentially decaying trend, and a combination of white and coloured noise obeying the power law:

Input Signal

Applying EMD to Extract Intrinsic Mode Functions (IMFs)

After setting the mean of the input signal to zero, we apply EMD to obtain the set of intrinsic mode functions (IMFs):

modes = emd_modes(x, sd_thresh=1e-4)
plot_modes(t, modes)

where the sd_thresh parameter is the threshold at which the sift of each IMF stops. In our example, we obtained seven EMD modes, six of which are oscillatory IMFs, and one is a non-oscillatory residual (usually, the number of EMD modes is about \(\log_2(N)\), where \(N\) is the number of data points in the input signal).

First EMD Mode

Estimating the Empirical Trend

The empirical trend of the signal is estimated using the scope.emd.emd_trend() function. This function identifies modes with periods exceeding a fraction of the total signal duration (denoted by the cutoff parameter) and the residual, combines them into an empirical trend of the input signal, and returns a new set of modes in which all modes have periods shorter than the cutoff while the last mode represents the signal’s trend.

By default, this cutoff is set to 0.4 of the total signal length, meaning that a mode with less than 2.5 oscillation cycles is considered part of the empirical trend.

modes = emd_trend(modes, t)
trend_emd = modes[:, -1]
plot_signal(t, trend_emd, 'Trend of the signal')

For our example, the empirical trend of the signal is found to form by the last EMD mode (the residual) only:

Empirical Trend

Hence, the detrended signal is:

Detrended Signal

Fourier Analysis of the Detrended Signal

Now we estimate the parameters of the superimposed noise by applying the scope.fourier.fit_fourier() function to the detrended signal. The function returns the FFT spectrum best-fitted by a power-law model, with powers of white (if present) and coloured noise and the power-law index of coloured noise as model parameters.

For our example, the FFT spectrum shows a combination of white and coloured noise components, with the power-law index of coloured noise being 1.1±0.3 and the ratio of the white to coloured noise energies about 0.3. The scope.fourier.fit_fourier() function also estimates the confidence interval of a given value (e.g., 95%, false alarm probability = 0.05). Fourier peaks outside this confidence interval are attributed to statistically significant oscillatory processes of non-noise origin.

fit_fft = fit_fourier(x, dt, fap=0.05)
plot_fft_spectrum(fit_fft)
FFT Spectrum

Computing the EMD Energy Spectrum

The EMD energy spectrum, i.e., the relationship between the EMD modal energy vs. dominant oscillation period, is computed using the scope.emd.emd_energy_spectrum() function:

emd_sp = emd_energy_spectrum(modes, t)
cutoff_period = 0.4 * len(x) * dt  # show cutoff period
plot_emd_spectrum(emd_sp, cutoff_period)
EMD Energy Spectrum

The vertical dashed line corresponds to the cutoff period adopted in the scope.emd.emd_trend() function; all modes beyond this line are considered components of the trend.

Computing Confidence Limits for the EMD Energy Spectrum

Using the power-law index and noise energy returned by scope.fourier.fit_fourier(), we compute the confidence limits of the EMD energy spectrum using the scope.emd.emd_noise_conf() function (separately for coloured noise and, if present, white noise):

# False alarm probability
fap = 0.05

# Confidence limits for coloured noise
conf_c = emd_noise_conf(t, alpha=alpha, period_min=2*dt,
                        period_max=N*dt, num_samples=500,
                        signal_energy=fit_fft['color_energy'], fap=fap)

# Confidence limits for white noise
if fit_fft['white_energy'] > 0:  # check if there is only coloured noise model
    conf_w = emd_noise_conf(t, alpha=0, period_min=2*dt,
                            period_max=N*dt, num_samples=500,
                            signal_energy=fit_fft['white_energy'], fap=fap)

Here, the false alarm probability (fap) is set to 0.05 (95% confidence). The scope.emd.emd_noise_conf`() function generates 500 independent noise samples with the same power-law index (alpha) and energy (signal_energy) as the input. The other parameters, period_min and period_max, define the range of periods over which confidence limits are computed.

Combining the upper and lower confidence limits for white and coloured noise components:

# Upper confidence limit for the combined noises
conf_up = conf_c['up'] + conf_w['up']

# Lower confidence limit for the combined noises
conf_down = conf_c['down'] + conf_w['down']

Visualizing the EMD energy spectrum with confidence limits:

# Plot EMD spectrum
plot_emd_spectrum(emd_sp, cutoff_period, conf_period, conf_up, conf_down, conf_mean, fap)
EMD Energy Spectrum with Confidence Limits

Here: - conf_mean represents the expected mean noise energy:

conf_mean = conf_c['mean_energy'] + conf_w['mean_energy']

  • conf_period is the array of oscillation periods for which confidence limits were computed: conf_period = conf_c['period']

Identifying Statistically Significant EMD Modes

The EMD modes beyond the confidence limits are considered significant, meaning they are not likely to be caused by random noise. In our example, only one mode is found to be significant, which appears consistent with the input oscillatory component of the original signal.

Significant Mode