.. 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_03_ensemblesift.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_03_ensemblesift.py: Ensemble sifting ================ This tutorial introduces the ensemble sift. This is a noise-assisted approach to sifting which overcomes some of the problems in the original sift. .. GENERATED FROM PYTHON SOURCE LINES 10-12 Lets make a simulated signal to get started. Here we generate two bursts of oscillations and some dynamic noise. .. GENERATED FROM PYTHON SOURCE LINES 12-43 .. code-block:: Python import emd import numpy as np import matplotlib.pyplot as plt sample_rate = 512 seconds = 2 time_vect = np.linspace(0, seconds, seconds*sample_rate) # Create an amplitude modulated burst am = -np.cos(2*np.pi*1*time_vect) * 2 am[am < 0] = 0 burst = am*np.sin(2*np.pi*42*time_vect) # Create some noise with increasing amplitude am = np.linspace(0, 1, sample_rate*seconds)**2 + .1 np.random.seed(42) n = am * np.random.randn(sample_rate*seconds,) # Signal is burst + noise x = burst + n # Quick summary figure plt.figure() plt.subplot(311) plt.plot(burst) plt.subplot(312) plt.plot(n) plt.subplot(313) plt.plot(x) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_001.png :alt: emd tutorial 01 sift 03 ensemblesift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 44-46 The standard sifting methods do not perform well on this signal. Here we run ``emd.sift.sift`` and plot the result. .. GENERATED FROM PYTHON SOURCE LINES 46-51 .. code-block:: Python imf_opts = {'sd_thresh': 0.05} imf = emd.sift.sift(burst+n, imf_opts=imf_opts) emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_002.png :alt: emd tutorial 01 sift 03 ensemblesift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 52-68 We can see that the 42Hz bursting activity is split across two or three IMFs. The first burst is split between IMFs 1 and 2 whereas the second burst is split between IMFs 2 and 3. As these events are all the same frequency, we want them to be in the the IMF to make subsequent spectrum analysis easier. The problem in this case is the dynamic noise. The peaks and troughs of the first burst are only slightly distorted by the additive noise (whose variance is low in the first half of the signal). Therefore most of the first burst appears in the first IMF with a small segment in the second IMF. In contrast, the high variance noise under the second burst is large enough to consitute a feature in itself. The second burst is completely pushed out of IMF1. The combined problems of noise sensitivity and transient signals can lead to poor sift results. .. GENERATED FROM PYTHON SOURCE LINES 70-72 Noise-assisted sifting ^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 74-86 One solution to these issues is to try and normalise the amount of noise through a signal by adding a small amount of white noise to the signal before sifting. We do this many times, creating an ensemble of sift processes each with a separate white noise added. The final set of IMFs is taken to be the average across the whole ensemble. This might seem like an odd solution. It relies on the effect of the additive white noise cancelling out across the whole ensemble, whilst the true signal (which is present across the whole ensemble) should survive the averaging. Here we try this out by extracting out first IMF five times. Once without noise and four times with noise. .. GENERATED FROM PYTHON SOURCE LINES 86-97 .. code-block:: Python imf = np.zeros((x.shape[0], 5)) # Standard sift imf[:, 0] = emd.sift.get_next_imf(x, **imf_opts)[0][:, 0] # Additive noise sifts noise_variance = .25 for ii in range(4): imf[:, ii+1] = emd.sift.get_next_imf(x + np.random.randn(x.shape[0],)*noise_variance, **imf_opts)[0][:, 0] .. GENERATED FROM PYTHON SOURCE LINES 98-99 Next we plot our IMF from the standard sift and the four noise-added ensemble sifts .. GENERATED FROM PYTHON SOURCE LINES 99-110 .. code-block:: Python plt.figure() for ii in range(5): plt.subplot(5, 1, ii+1) plt.plot(imf[:, ii]) if ii == 0: plt.ylabel('Normal') else: plt.ylabel('Ens. {0}'.format(ii)) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_003.png :alt: emd tutorial 01 sift 03 ensemblesift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 111-118 We can see that the first burst appears clearly in the first IMF for the standard sift as before. This burst is suppressed in each of the ensembles. Some of it remains but crucially different parts of the signal are attenuated in each of the four ensembles. The average of these four noise-added signals shows a strong suppression of the first burst. .. GENERATED FROM PYTHON SOURCE LINES 118-127 .. code-block:: Python plt.figure() plt.subplot(211) plt.plot(imf[:, 0]) plt.ylabel('Normal') plt.subplot(212) plt.plot(imf[:, 1:].mean(axis=1)) plt.ylabel('Ens. Avg') .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_004.png :alt: emd tutorial 01 sift 03 ensemblesift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(44.222222222222214, 0.5, 'Ens. Avg') .. GENERATED FROM PYTHON SOURCE LINES 128-132 Some of the first burst still remains. We can increase the attenuation by increasing the variance of the added noise. Here, we run the ensemble of four again but we increase the ``noise_variance`` from `0.25` to `1`. .. GENERATED FROM PYTHON SOURCE LINES 132-143 .. code-block:: Python imf = np.zeros((x.shape[0], 5)) # Standard sift imf[:, 0] = emd.sift.get_next_imf(x, **imf_opts)[0][:, 0] # Additive noise sifts noise_variance = 1 for ii in range(4): imf[:, ii+1] = emd.sift.get_next_imf(x + np.random.randn(x.shape[0],)*noise_variance, **imf_opts)[0][:, 0] .. GENERATED FROM PYTHON SOURCE LINES 144-146 If we plot the average of this new ensemble, we can see that the first burst is now completely removed from the first IMF. .. GENERATED FROM PYTHON SOURCE LINES 146-153 .. code-block:: Python plt.figure() plt.subplot(211) plt.plot(imf[:, 0]) plt.subplot(212) plt.plot(imf[:, 1:].mean(axis=1)) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_005.png :alt: emd tutorial 01 sift 03 ensemblesift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 154-156 Ensemble Sifting ^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 158-173 This process is implemented for a whole sifting run in ``emd.sift.ensemble_sift``. This function works much like ``emd.sift.sift`` with a few extra options for controlling the noise ensembles. * `nensembles` defines the number of parallel sifts to compute * `ensemble_noise` defines the noise variance relative to the standard deviation of the input signal * `nprocesses` allow for parallel processing of the ensembles to speed up computation Next, we call ``emd.sift.ensemble_sift`` to run through a whole sift of our signal. Note that it is ofen a good idea to limit the total number of IMFs in an ensemble_sift. If we allow the sift to compute all possible IMFs then there is a chance that some processes in the ensemble might complete with one more or one less IMF than the others - this will break the averaging so here, we only compute the first five IMFs. .. GENERATED FROM PYTHON SOURCE LINES 175-176 First, we compute the sift with `2` ensembles. .. GENERATED FROM PYTHON SOURCE LINES 176-180 .. code-block:: Python imf = emd.sift.ensemble_sift(burst+n, max_imfs=5, nensembles=2, nprocesses=6, ensemble_noise=1, imf_opts=imf_opts) emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_006.png :alt: emd tutorial 01 sift 03 ensemblesift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 181-185 This does a reasonable job. Certainly better than the standard sift we ran at the start of the tutorial. There is still some mixing of signals visible though. To help with this we next compute the ensemble sift with four separate noise added sifts. .. GENERATED FROM PYTHON SOURCE LINES 185-189 .. code-block:: Python imf = emd.sift.ensemble_sift(burst+n, max_imfs=5, nensembles=4, nprocesses=6, ensemble_noise=1, imf_opts=imf_opts) emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_007.png :alt: emd tutorial 01 sift 03 ensemblesift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 190-192 This is better again but still could be improved. Finally, we run a sift with an ensemble of 24 separate sifts. .. GENERATED FROM PYTHON SOURCE LINES 192-196 .. code-block:: Python imf = emd.sift.ensemble_sift(burst+n, max_imfs=5, nensembles=24, nprocesses=6, ensemble_noise=1, imf_opts=imf_opts) emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_008.png :alt: emd tutorial 01 sift 03 ensemblesift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_03_ensemblesift_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 197-200 This does a good job. Both bursts are clearly recovered with smooth amplitude modulations and the noise in the first IMF shows the smooth increase in variance that we defined at the start. .. GENERATED FROM PYTHON SOURCE LINES 202-204 Further Reading & References ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 206-215 Wu, Z., & Huang, N. E. (2009). Ensemble Empirical Mode Decomposition: A Noise-Assisted Data Analysis Method. Advances in Adaptive Data Analysis, 1(1), 1–41. https://doi.org/10.1142/s1793536909000047 Wu, Z., & Huang, N. E. (2004). A study of the characteristics of white noise using the empirical mode decomposition method. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 460(2046), 1597–1611. https://doi.org/10.1098/rspa.2003.1221 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.201 seconds) .. _sphx_glr_download_emd_tutorials_01_sifting_emd_tutorial_01_sift_03_ensemblesift.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_03_ensemblesift.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: emd_tutorial_01_sift_03_ensemblesift.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_