Theory
========

emd_period_energy
**************

As mentioned in the example section, the total energy and (dominant) period of \
each EMD mode are required for constructing an EMD energy spectrum. The total \
modal energy is estimated by summing up squares of instantaneous amplitudes of \
each EMD mode. The dominant period of each EMD mode is estimated by best-fitting \
the global wavelet spectrum `Torrence et al. (1998) <https://psl.noaa.gov/people/gilbert.p.compo/Torrence_compo1998.pdf>`_ \
of the mode with a Gaussian + Parabolic function, performed in the :func:`scope.emd.emd_period_energy` \
function. The position and standard deviation of the Gaussian peak are used for \
the dominant EMD modal period and the uncertainty of this estimation.

fit_fourier
**************

In the :func:`scope.fourier.fit_fourier` function, we fit the FFT spectrum by a \
power-law model in log-log scale to extract the power-law index and energy of \
the noise component of the signal. Firstly, we must note that, at each Fourier \
frequency, the Fourier power :math:`I(f_{j})` follows a chi-squared distribution \
with 2 degrees of freedom, denoted as:

.. math::
   I(f_{j}) = P(f_{j}) \chi_{2}^{2}/2

where :math:`P(f_{j})` is the true power spectrum, and :math:`\chi_{2}^{2}` is \
a random variable distributed :math:`chi^{2}` with 2 degrees of freedom. Since \
the least squares method assumes that the input data set is Gaussian-distributed, \
we cannot directly apply this method to best-fit the FFT power spectrum. Instead, \
we should consider the mean of the :math:`\chi_{2}^{2}/2` term. In log scale, \
:math:`\left\langle \mathrm{log}(\chi^{2}_{2}/2) \right\rangle = -0.25068` \
(`Vaughan (2005) <https://doi.org/10.1051/0004-6361:20041453>`_). This term \
corresponds to the bias that will be introduced to the fitting if one directly \
implements the least squares method. Hence, we shall include this term in the \
model function such that the least squares fitting will not be 'biased'. We \
also note that the value of this bias term is independent of the choice of \
normalization of the FFT power spectrum.

The power law model we used in the :func:`scope.fourier.fit_fourier` function \
is a superposition of white and coloured noise components, given by:

.. math::
   P(f) = P_{c}(f) + P_{w}(f) = Z_{c} f^{-\alpha} + Z_{w},

where :math:`Z_{c}` and :math:`Z_{w}` are the proportionality constants of \
coloured and white noises, respectively, and  :math:`\alpha` is the power law \
index of coloured noise. 

After obtaining the proportionality constants from the debiased least squares \
fit, we can estimate the energy of each noise type, :math:`E_{c/w}` using:

.. math::
   E_{c/w} \propto nf \times Z_{c/w},

where :math:`nf` is the number of Fourier frequencies, which does not include \
0 Hz and the Nyquist frequency. And estimate the confidence limit for a given \
false alarm probability (fap) as :math:`-\ln\left(1-(1-\mathrm{fap})^{1/nf}\right)\times P(f_{j})`.

See `fft_fit_example.py <https://github.com/Warwick-Solar/scope/blob/main/examples/fft_fit_example.py>`_ \
for an example use of the :func:`scope.fourier.fit_fourier` function.

emd_noise_conf
**************

`Kolotkov et al. (2016) <https://doi.org/10.1051/0004-6361/201628306>`_ showed \
that the dyadic property of EMD (the center frequencies of consecutive IMFs \
tend to have a ratio close to 2) results in the following relationship between \
modal energy and modal period of a noise signal:

.. math::
   E_{m}P_{m}^{1-\alpha} = \text{const,}

where the parameter :math:`\alpha` is the power-law index used for characterizing \
the colour of noise in the Fourier analysis.

See the following image for a visual representation of the relationship:

.. image:: ../_static/mc_emd_spectra.png

See `energy_period_relation.py <https://github.com/Warwick-Solar/scope/blob/main/examples/energy_period_relation.py>`_ demonstrating the relationship between :math:`E_{m}` and :math:`P_{m}` for various types of noise (values of :math:`\alpha`). Here, we set the false alarm probability to 0.05.

It was shown that the energy of each EMD mode, :math:`E_m` has a chi-squared \
distribution with :math:`k` degrees of freedom (DoF). In contrast to the Fourier \
power, for which the number of DoF is 2 for all Fourier harmonics, the number of \
DoF :math:`k` of EMD modal energy :math:`E_m` is usually :math:`>2` and varies with \
the mode number (hence, with the period). Thus, for given noise parameters (i.e., /
the energies of the white and coloured components + the power-law index of the \
coloured component) determined with the :func:`scope.fourier.fit_fourier` function, \
the :func:`scope.emd.emd_noise_conf` function first estimates the number of DoF :math:`k` \
over the entire range of EMD modal periods and then estimates the confidence limits using \
the percent-point function of the chi-squared distribution with :math:`k` DoF. 

The :func:`scope.emd.emd_noise_conf` function generates 500 (by default) independent /
noise samples with the same power-law index and energy as the input and performs the /
EMD analysis on them. It extracts the dominant period and modal energy for each IMF /
by calling the :func:`scope.emd.emd_period_energy` function. The \
:func:`scope.emd.emd_noise_fit` function fits the chi-squared distribution to /
the histogram of modal energy :math:`E_m` for each mode number to extract the /
mean energy and :math:`k`. We obtain the mean period, mean energy, and number of /
DoF :math:`k` for each mode number. The empirically established relationships between /
both mean energy vs mean period and :math:`k` vs mean period are best-fitted with /
power-law functions (linear functions in log-log scale). These best-fit functions /
are then used to construct the confidence limits over the whole range of EMD modal /
periods. As, for :math:`E_m`, the number of DoF :math:`k>2`, resulting in a /
non-monotonic distribution of :math:`E_m`, we get two confidence limits (upper /
and lower) in the EMD energy spectrum. In practice, the modes above the upper /
confidence limit are of greater interest as more energetic. This analysis does not /
apply to the very first IMF (the EMD mode with the shortest timescale) as its /
energy does not obey the chi-squared distribution.