.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "emd_tutorials/01_sifting/emd_tutorial_01_sift_02_siftindetail.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_emd_tutorials_01_sifting_emd_tutorial_01_sift_02_siftindetail.py: The sift in detail ================== Here, we will run through the different steps of the sift and get to know some of the lower-level functions which are used by the core sift functions. There are four levels of functions which are used in the sift. We will take a look at each of these steps in turn using a simulated time-series. .. GENERATED FROM PYTHON SOURCE LINES 13-14 Lets make a simulated signal to get started. .. GENERATED FROM PYTHON SOURCE LINES 14-39 .. code-block:: Python import emd import numpy as np import matplotlib.pyplot as plt sample_rate = 1000 seconds = 10 num_samples = sample_rate*seconds time_vect = np.linspace(0, seconds, num_samples) freq = 5 # Change extent of deformation from sinusoidal shape [-1 to 1] nonlinearity_deg = .25 # Change left-right skew of deformation [-pi to pi] nonlinearity_phi = -np.pi/4 # Create a non-linear oscillation x = emd.simulate.abreu2010(freq, nonlinearity_deg, nonlinearity_phi, sample_rate, seconds) x += np.cos(2 * np.pi * 1 * time_vect) # Add a simple 1Hz sinusoid x -= np.sin(2 * np.pi * 2.2e-1 * time_vect) # Add a very slow cycle # sphinx_gallery_thumbnail_number = 7 .. GENERATED FROM PYTHON SOURCE LINES 40-42 Sifting ^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 44-49 The top-level of options configure the sift itself. These options vary between the type of sift that is being performed and options don't generalise between different variants of the sift. Here we will run a standard sift on our simulation. .. GENERATED FROM PYTHON SOURCE LINES 49-59 .. code-block:: Python # Get the default configuration for a sift config = emd.sift.get_config('sift') # Adjust the threshold for accepting an IMF config['imf_opts/sd_thresh'] = 0.05 imf = emd.sift.sift(x) emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_001.png :alt: emd tutorial 01 sift 02 siftindetail :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 60-71 .. code-block:: Python # Get the default configuration for a sift config = emd.sift.get_config('sift') # Adjust the threshold for accepting an IMF config['imf_opts/sd_thresh'] = 0.05 config['extrema_opts/method'] = 'rilling' imf = emd.sift.sift(x) emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_002.png :alt: emd tutorial 01 sift 02 siftindetail :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 74-82 Internally the ``sift`` function calls a set of lower level functions to extract the IMFs. These functions are call in a hierarchy when you run ``sift`` it will call ``get_next_imf`` behind the scenes. Similarly, ``get_next_imf`` makes use of ``interp_envelope`` and so on. - ``get_next_imf`` extracts the fastest IMF from an input signal - ``interp_envelope`` find the interpolated envelope of a signal. - ``get_padded_extrema`` identifies the location and magnitude of signal extrema with optional padding. .. GENERATED FROM PYTHON SOURCE LINES 85-86 We will run through each of these functions in now, giving some examples of their use and options. .. GENERATED FROM PYTHON SOURCE LINES 88-90 IMF Extraction ^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 92-107 After the top-level sift function, the next layer is IMF extraction as implemented in ``emd.sift.get_next_imf``. This uses the envelope interpolation and extrema detection to carry out the sifting iterations on a time-series to return a single intrinsic mode function. This is the main function used when implementing novel types of sift. For instance, the ensemble sift uses this ``emd.sift.get_next_imf`` to extract IMFs from many repetitions of the signal with small amounts of noise added. Similarly the mask sift calls ``emd.sift.get_next_imf`` after adding a mask signal to the data. Here we call ``get_next_imf`` repeatedly on a signal and its residuals to implement a very simple sift. We extract the first IMF, subtract it from the data and then extract the second and third IMFs. We then plot the original signal, the IMFs and the residual. .. GENERATED FROM PYTHON SOURCE LINES 107-129 .. code-block:: Python # Extract the options for get_next_imf - these can be customised here at this point. imf_opts = config['imf_opts'] # Extract first IMF from the signal imf1, continue_sift = emd.sift.get_next_imf(x[:, None], **imf_opts) # Extract second IMF from the signal with the first IMF removed imf2, continue_sift = emd.sift.get_next_imf(x[:, None]-imf1, **imf_opts) # Extract third IMF from the signal with the first and second IMFs removed imf3, continue_sift = emd.sift.get_next_imf(x[:, None]-imf1-imf2, **imf_opts) # The residual is the signal component left after removing the IMFs residual = x[:, None] - imf1 - imf2 - imf3 # Contactenate our IMFs into one array imfs_manual = np.c_[imf1, imf2, imf3, residual] # Visualise emd.plotting.plot_imfs(imfs_manual) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_003.png :alt: emd tutorial 01 sift 02 siftindetail :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 130-131 These IMFs should be identical to the IMFs obtained using ``emd.sift.sift`` above. .. GENERATED FROM PYTHON SOURCE LINES 134-136 Envelope interpolation ^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 138-143 A large part of IMF exraction is the computation of an upper and lower envelope of the signal. This is done through interpolation using ``emd.sift.interp_envelope`` and the options in the ``envelope`` section of the config. .. GENERATED FROM PYTHON SOURCE LINES 143-162 .. code-block:: Python # Extract envelope options env_opts = config['envelope_opts'] # Compute upper and lower envelopes upper_env = emd.sift.interp_envelope(x, mode='upper', **env_opts) lower_env = emd.sift.interp_envelope(x, mode='lower', **env_opts) # Compute average envelope avg_env = (upper_env+lower_env) / 2 # Visualise plt.figure(figsize=(12, 6)) plt.plot(x, 'k') plt.plot(upper_env, 'r') plt.plot(lower_env, 'b') plt.plot(avg_env, 'g') plt.legend(['Signal', 'Upper Envelope', 'Lower Envelope', 'Avg. Envelope']) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_004.png :alt: emd tutorial 01 sift 02 siftindetail :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 163-166 Subtracting the upper and lower envelopes from the signal removes slow dynamics from the signal. Next, we subtract the average envelope from our signal. .. GENERATED FROM PYTHON SOURCE LINES 166-176 .. code-block:: Python # Plot the signal with the average of the upper and lower envelopes subtracted alongside IMF1 plt.figure(figsize=(12, 6)) plt.subplot(211) plt.plot(x-avg_env, 'k') plt.title('Signal-Average Envelope') plt.subplot(212) plt.plot(imf1, 'b') plt.title('IMF1') .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_005.png :alt: Signal-Average Envelope, IMF1 :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'IMF1') .. GENERATED FROM PYTHON SOURCE LINES 177-186 Note that the simple subtraction is very similar to IMF1 as extracted above. In real data, several iterations of envelope computation and subtraction may be required to identify a well-formed IMF. In this case there is a small amplitude error in the IMF at the very start. This is due to uncertainty in the envelope interpolation at the edges. This can sometimes be reduced by changing the interpolation and extrema padding options but is hard to completely overcome. It is often sensible to treat the first and last couple of cycles in an IMF with caution. .. GENERATED FROM PYTHON SOURCE LINES 188-190 Extrema detection and padding ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 192-196 Finally, the lowest-level functions involve extrema detection and padding as implemented in the ``emd.sift.get_padded_extrema`` function. This is a simple function which identifies extrema using ``scipy.signal``. Here we identify peaks and troughs without any padding applied. .. GENERATED FROM PYTHON SOURCE LINES 196-207 .. code-block:: Python max_locs, max_mag = emd.sift.get_padded_extrema(x, pad_width=0, mode='peaks') min_locs, min_mag = emd.sift.get_padded_extrema(x, pad_width=0, mode='troughs') plt.figure(figsize=(12, 3)) plt.plot(x, 'k') plt.plot(max_locs, max_mag, 'or') plt.plot(min_locs, min_mag, 'ob') plt.legend(['Signal', 'Maxima', 'Minima']) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_006.png :alt: emd tutorial 01 sift 02 siftindetail :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 208-211 Extrema padding is used to stablise the envelope at the edges of the time-series. The ``emd.sift.get_padded_extrema`` function identifies and pads extrema in a time-series. This calls the ``emd.sift.find_extrema`` internally. .. GENERATED FROM PYTHON SOURCE LINES 211-236 .. code-block:: Python plt.figure(figsize=(12, 7)) # Method 1 - np.pad max_locs, max_mag = emd.sift.get_padded_extrema(x, pad_width=2, mode='peaks', method='numpypad') min_locs, min_mag = emd.sift.get_padded_extrema(x, pad_width=2, mode='troughs', method='numpypad') plt.subplot(211) plt.plot(x, 'k') plt.plot(max_locs, max_mag, 'or') plt.plot(min_locs, min_mag, 'ob') plt.legend(['Signal', 'Maxima', 'Minima']) plt.title('default numpy.pad method') # Method 2 - Rillings method max_locs, max_mag = emd.sift.get_padded_extrema(x, pad_width=4, mode='peaks', method='rilling') min_locs, min_mag = emd.sift.get_padded_extrema(x, pad_width=4, mode='troughs', method='rilling') plt.subplot(212) plt.plot(x, 'k') plt.plot(max_locs, max_mag, 'or') plt.plot(min_locs, min_mag, 'ob') plt.legend(['Signal', 'Maxima', 'Minima']) plt.title('Rillings Method') .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_007.png :alt: default numpy.pad method, Rillings Method :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Rillings Method') .. GENERATED FROM PYTHON SOURCE LINES 237-240 Extrema padding is used to stablise the envelope at the edges of the time-series. The ``emd.sift.get_padded_extrema`` function identifies and pads extrema in a time-series. This calls the ``emd.sift.find_extrema`` internally. .. GENERATED FROM PYTHON SOURCE LINES 240-263 .. code-block:: Python plt.figure(figsize=(12, 7)) env, extr = emd.sift.interp_envelope(x, ret_extrema=True, extrema_opts={'method': 'numpypad'}, trim=False) plt.subplot(211) plt.plot(x, 'k') plt.plot(extr[0], extr[1], 'or') plt.plot(extr[2], extr[3], 'ob') plt.plot(np.arange(extr[0][0], extr[0][-1]), env[1]) plt.plot(np.arange(extr[2][0], extr[2][-1]), env[0]) plt.legend(['Signal', 'Maxima', 'Minima', 'Upper Env', 'Lower Env']) plt.title('np.pad Method') env, extr = emd.sift.interp_envelope(x, ret_extrema=True, extrema_opts={'method': 'rilling'}, trim=False) plt.subplot(212) plt.plot(x, 'k') plt.plot(extr[0], extr[1], 'or') plt.plot(extr[2], extr[3], 'ob') plt.plot(np.arange(extr[0][0], extr[0][-1]), env[1]) plt.plot(np.arange(extr[2][0], extr[2][-1]), env[0]) plt.legend(['Signal', 'Maxima', 'Minima', 'Upper Env', 'Lower Env']) plt.title('Rillings Method') .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_008.png :alt: np.pad Method, Rillings Method :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Rillings Method') .. GENERATED FROM PYTHON SOURCE LINES 267-280 The extrema detection and padding arguments are specified in the config dict under the extrema, mag_pad and loc_pad keywords. These are passed directly into ``emd.sift.get_padded_extrema`` when running the sift. The default extrema padding uses Rilling's method, but we can switch this to use a built in numpy function ``np.pad``. The ``mag_pad`` and ``loc_pad`` dictionaries are passed into np.pad to define the padding in the y-axis (extrema magnitude) and x-axis (extrema time-point) respectively. Note that ``np.pad`` takes a mode as a positional argument - this must be included as a keyword argument here. Lets try customising the extrema padding. First we get the 'extrema' options from a nested config then try changing a couple of options .. GENERATED FROM PYTHON SOURCE LINES 280-323 .. code-block:: Python ext_opts = config['extrema_opts'] ext_opts['method'] = 'numpypad' # The default options max_locs, max_mag = emd.sift.get_padded_extrema(x, mode='peaks', **ext_opts) min_locs, min_mag = emd.sift.get_padded_extrema(x, mode='troughs', **ext_opts) plt.figure(figsize=(12, 12)) plt.subplot(311) plt.plot(x, 'k') plt.plot(max_locs, max_mag, 'or') plt.plot(min_locs, min_mag, 'ob') plt.legend(['Signal', 'Maxima', 'Minima']) plt.title('Default') # Increase the pad width to 5 extrema ext_opts['pad_width'] = 5 max_locs, max_mag = emd.sift.get_padded_extrema(x, mode='peaks', **ext_opts) min_locs, min_mag = emd.sift.get_padded_extrema(x, mode='troughs', **ext_opts) plt.subplot(312) plt.plot(x, 'k') plt.plot(max_locs, max_mag, 'or') plt.plot(min_locs, min_mag, 'ob') plt.legend(['Signal', 'Maxima', 'Minima']) plt.title('Increased pad width') # Change the y-axis padding to 'reflect' rather than 'median' (this option is # for illustration and not recommended for actual sifting....) ext_opts['mag_pad_opts']['mode'] = 'reflect' del ext_opts['mag_pad_opts']['stat_length'] max_locs, max_mag = emd.sift.get_padded_extrema(x, mode='peaks', **ext_opts) min_locs, min_mag = emd.sift.get_padded_extrema(x, mode='troughs', **ext_opts) plt.subplot(313) plt.plot(x, 'k') plt.plot(max_locs, max_mag, 'or') plt.plot(min_locs, min_mag, 'ob') plt.legend(['Signal', 'Maxima', 'Minima']) plt.title('Reflected extrema and increased pad width') .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_009.png :alt: Default, Increased pad width, Reflected extrema and increased pad width :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_02_siftindetail_009.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Reflected extrema and increased pad width') .. GENERATED FROM PYTHON SOURCE LINES 324-326 Further Reading & References ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 328-333 Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H. H., Zheng, Q., … Liu, H. H. (1998). The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 454(1971), 903–995. https://doi.org/10.1098/rspa.1998.0193 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.156 seconds) .. _sphx_glr_download_emd_tutorials_01_sifting_emd_tutorial_01_sift_02_siftindetail.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: emd_tutorial_01_sift_02_siftindetail.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: emd_tutorial_01_sift_02_siftindetail.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_