Tutorial
===================

The example described below is provided in `emd_example.py <https://github.com/Warwick-Solar/scope/blob/main/examples/emd_example.py>`_.

First, import all  nedeed libraries:

.. code-block:: python

    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:

.. image:: ../_static/input_signal.png
   :align: center
   :alt: 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):

.. code-block:: python

   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 \
:math:`\log_2(N)`, where :math:`N` is the number of data points in the input signal).

.. image:: ../_static/1st_EMD.png
   :align: center
   :alt: First EMD Mode

Estimating the Empirical Trend
------------------------------

The empirical trend of the signal is estimated using the :func:`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.

.. code-block:: python

   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:

.. image:: ../_static/trend_signal.png
   :align: center
   :alt: Empirical Trend

Hence, the detrended signal is:

.. image:: ../_static/detrended_signal.png
   :align: center
   :alt: Detrended Signal

Fourier Analysis of the Detrended Signal
----------------------------------------

Now we estimate the parameters of the superimposed noise by applying the :func:`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 :func:`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.

.. code-block:: python

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

.. image:: ../_static/FFT_spectrum.png
   :align: center
   :alt: 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 :func:`scope.emd.emd_energy_spectrum` function:

.. code-block:: python

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

.. image:: ../_static/emd_spectrum.png
   :align: center
   :alt: EMD Energy Spectrum

The vertical dashed line corresponds to the cutoff period adopted in the :func:`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 :func:`scope.fourier.fit_fourier`, \
we compute the confidence limits of the EMD energy spectrum using the \
:func:`scope.emd.emd_noise_conf` function (separately for coloured noise and, if present, white noise):

.. code-block:: python

   # 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 :func:`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:

.. code-block:: python

   # 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:

.. code-block:: python

   # Plot EMD spectrum
   plot_emd_spectrum(emd_sp, cutoff_period, conf_period, conf_up, conf_down, conf_mean, fap)

.. image:: ../_static/emd_spectrum_with_conf.png
   :align: center
   :alt: 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.

.. image:: ../_static/significant_mode.png
   :align: center
   :alt: Significant Mode