.. 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_05_iterated_masksift.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_05_iterated_masksift.py: Iterated Masked sifting ======================== Masking signals are a powerful method for reducing the mode mixing problem when running EMD on noisy or transient signals. However, the standard methods for selecting masks often perform poorly in real data and selecting masks by hand can be cumbersome and prone to error. This tutorial introduces iterated masking EMD (itEMD), a recent sifting method that automates mask signal frequencies based on the dynamics present in the data. We will show how it automatically dentifies oscillations and minimises mode mixing without prespecification of masking signals. The method is published in the following paper: Marco S. Fabus, Andrew J. Quinn, Catherine E. Warnaby, and Mark W. Woolrich (2021). Automatic decomposition of electroptysiological data into distinct nonsinusoidal oscillatory modes. Journal of Neurophysiology 2021 126:5, 1670-1684. https://doi.org/10.1152/jn.00315.2021 .. GENERATED FROM PYTHON SOURCE LINES 24-26 Issues with Mask Selection ^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 28-29 Let's start by simulating a signal with a 2Hz nonsinusoidal oscillation and 30Hz beta bursts. .. GENERATED FROM PYTHON SOURCE LINES 29-62 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt import emd np.random.seed(0) sample_rate = 256 seconds = 3 num_samples = sample_rate*seconds time_vect = np.linspace(0, seconds, num_samples) # Create an amplitude modulation freq = 2 am = np.sin(2*np.pi*freq*time_vect) am[am < 0] = 0 # Non-sinusoidal intermittend 2Hz signal x1 = np.sin(np.sin(np.sin(np.sin(np.sin(np.sin(np.sin(2*np.pi*freq*time_vect))))))) / 0.5 # Bursting 30Hz oscillation x2 = am * 0.5 * np.cos(2*np.pi*30*time_vect) # Sum them together x = x1 + x2 # Quick summary figure plt.figure(figsize=(8, 4)) plt.plot(time_vect, x) plt.xlim(0, 3) # sphinx_gallery_thumbnail_number = 6 .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_001.png :alt: emd tutorial 01 sift 05 iterated masksift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (0.0, 3.0) .. GENERATED FROM PYTHON SOURCE LINES 63-64 Let's run the default mask sift, where mask frequency is determined from zero crossings of the signal. .. GENERATED FROM PYTHON SOURCE LINES 64-68 .. code-block:: Python imf = emd.sift.mask_sift(x, max_imfs=4) emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_002.png :alt: emd tutorial 01 sift 05 iterated masksift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 69-74 We can see that the result is not great. Bursting beta is in IMF-1, but the non-sinusoidal oscillation is split between IMF-2 and IMF-4. This is because the waveform is highly non-sinusoidal and spans a lot of frequencies, confusing the default mask sift. Let's try the same thing, but now with custom mask frequencies. .. GENERATED FROM PYTHON SOURCE LINES 74-79 .. code-block:: Python mask = np.array([60, 30, 24, 2, 1, 0.5])/sample_rate imf = emd.sift.mask_sift(x, max_imfs=5, mask_freqs=mask) emd.plotting.plot_imfs(imf, time_vect=time_vect) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_003.png :alt: emd tutorial 01 sift 05 iterated masksift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 80-87 Better! Most of the 2Hz oscillation is now in IMF-4 and most of bursting beta in IMF-2. However, it's still not perfect - beta is slightly split between IMF-1 and IMF-2 and we have an unnecessary 'empty' mode IMF-3. Finding a mask that balances all of these issues can be an arduous manual process. How could we do this automatically? Let's come back to our original sift and look at the mean instantaneous frequencies of the modes weighted by instantaneous amplitude. .. GENERATED FROM PYTHON SOURCE LINES 87-92 .. code-block:: Python imf = emd.sift.mask_sift(x, max_imfs=4) IP,IF,IA = emd.spectra.frequency_transform(imf, sample_rate, 'nht') print(np.average(IF, 0, weights=IA**2)) .. rst-class:: sphx-glr-script-out .. code-block:: none [27.04427481 6.20693296 1.63994241 1.9635715 ] .. GENERATED FROM PYTHON SOURCE LINES 93-95 Automated Mask Selection ^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 97-100 Even though the modes are split, frequencies of the first two modes are close to the ground truths of 30Hz and 2Hz that we know. What happens if we take these frequencies as the mask for masked sift? .. GENERATED FROM PYTHON SOURCE LINES 100-110 .. code-block:: Python # Compute the default mask sift imf = emd.sift.mask_sift(x, max_imfs=4) # Take the mean IF as new mask, compute mask sift again and plot IP,IF,IA = emd.spectra.frequency_transform(imf, sample_rate, 'nht') mask_1 = np.average(IF, 0, weights=IA**2) / sample_rate imf = emd.sift.mask_sift(x, max_imfs=4, mask_freqs=mask_1) emd.plotting.plot_imfs(imf, time_vect=time_vect) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_004.png :alt: emd tutorial 01 sift 05 iterated masksift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 111-115 We have now removed the empty mode and the bursting beta in is entirely IMF-1. It's certainly an improvement on the default sift. Let's repeat this process ten times and track what happens to average frequencies of the first two IMFs. .. GENERATED FROM PYTHON SOURCE LINES 115-141 .. code-block:: Python # Compute the default mask sift imf, mask_0 = emd.sift.mask_sift(x, max_imfs=4, ret_mask_freq=True) # Take the mean IF as new mask, compute mask sift again and plot IP,IF,IA = emd.spectra.frequency_transform(imf, sample_rate, 'nht') mask = np.average(IF, 0, weights=IA**2) / sample_rate # Save all masks: mask_all = np.zeros((3, 12)) mask_all[:, 0] = mask_0[:3] mask_all[:, 1] = mask[:3] for i in range(10): imf = emd.sift.mask_sift(x, max_imfs=4, mask_freqs=mask) IP,IF,IA = emd.spectra.frequency_transform(imf, sample_rate, 'nht') mask = np.average(IF, 0, weights=IA**2) / sample_rate mask_all[:, i+2] = mask[:3] plt.figure() for i in range(2): plt.plot(mask_all[i, :]*sample_rate) plt.xlabel('Iteration #') plt.ylabel('Mask frequency [Hz]') plt.show() .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_005.png :alt: emd tutorial 01 sift 05 iterated masksift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 142-147 The iteration process rapidly converged on the ground truth frequencies, 2Hz and 30Hz! This is the essence of iterated masking EMD (itEMD). By finding the equilibrium between masks and IMF frequencies, we automatically extract oscillations of interest in a data-driven way. Let's apply itEMD directly and show the result. .. GENERATED FROM PYTHON SOURCE LINES 147-151 .. code-block:: Python imf = emd.sift.iterated_mask_sift(x, sample_rate=sample_rate, max_imfs=3) emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_006.png :alt: emd tutorial 01 sift 05 iterated masksift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 152-158 And that's it! We have identified both the bursting beta and a non-sinusoidal oscillation automatically and without 'empty' modes between them. This decomposition has a clear Hilbert-Huang transform in the amplitude modulations of the 30Hz signal and within-cycle instantaneous frequency variability of the 2Hz signal are clearly visible. .. GENERATED FROM PYTHON SOURCE LINES 158-168 .. code-block:: Python IP, IF, IA = emd.spectra.frequency_transform(imf, sample_rate, 'nht') f, hht = emd.spectra.hilberthuang(IF, IA, (1, 50, 49, 'log'), sum_time=False) plt.figure(figsize=(8,4)) ax = plt.subplot(111) emd.plotting.plot_hilberthuang(hht, time_vect, f, log_y=True, cmap='ocean_r', ax=ax, time_lims=(0.5, 2)) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_007.png :alt: Hilbert-Huang Transform :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 169-171 Oscillations in Noise ^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 173-176 Lets take a look at how itEMD performs on a couple of noisy signals. Here we look at a dynamic 15Hz oscillation embedded in white noise. The itEMD is able to isolate the oscillation into a single component. .. GENERATED FROM PYTHON SOURCE LINES 176-189 .. code-block:: Python # Define and simulate a simple signal peak_freq = 15 sample_rate = 256 seconds = 3 noise_std = .4 x = emd.simulate.ar_oscillator(peak_freq, sample_rate, seconds, noise_std=noise_std, random_seed=42, r=.96)[:, 0] x = x*1e-4 imf = emd.sift.iterated_mask_sift(x, sample_rate=sample_rate, max_imfs=3) emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_008.png :alt: emd tutorial 01 sift 05 iterated masksift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 190-193 And, finally, we look at the how iterated mask sifting handles the signal from the Holospectrum tutorial. This case contains several oscillations with complex amplitude modulations. These are nicely separated into two IMFs. .. GENERATED FROM PYTHON SOURCE LINES 193-210 .. code-block:: Python seconds = 4 sample_rate = 200 t = np.linspace(0, seconds, seconds*sample_rate) # First we create a slow 4.25Hz oscillation with a 0.5Hz amplitude modulation slow = np.sin(2*np.pi*5*t) * (.5+(np.cos(2*np.pi*.5*t)/2)) # Second, we create a faster 37Hz oscillation that is amplitude modulated by the first. fast = .5*np.sin(2*np.pi*37*t) * (slow+(.5+(np.cos(2*np.pi*.5*t)/2))) # We create our signal by summing the oscillation and adding some noise x = slow+fast + np.random.randn(*t.shape)*.05 imf = emd.sift.iterated_mask_sift(x, sample_rate=sample_rate, max_imfs=5) emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_009.png :alt: emd tutorial 01 sift 05 iterated masksift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_05_iterated_masksift_009.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 211-225 Depending on the nature of the data, some itEMD options may need to be modified. For instance, the max_iter argument can be increased if itEMD reaches the maximum number of iterations, or iter_th can be modified to make the convergence criterion more stringent. To isolate very transient and infrequent oscillations, it may be a good idea to try different instantaneous amplitude weighting for the iteration process (keyword argument w_method). Finally, if the data is non-stationary such that different masks might be appropriate for different segments of the signal then the time-series should first be segmented and itEMD applied to quasi-stationary segments as it assumes one equilibrium can be found. For more information, see the documentation for emd and Fabus et al (2021) Journal of Neurophysiology.. .. GENERATED FROM PYTHON SOURCE LINES 227-229 Further Reading & References ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 231-235 Marco S. Fabus, Andrew J. Quinn, Catherine E. Warnaby, and Mark W. Woolrich (2021). Automatic decomposition of electroptysiological data into distinct nonsinusoidal oscillatory modes. Journal of Neurophysiology 2021 126:5, 1670-1684. https://doi.org/10.1152/jn.00315.2021 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.602 seconds) .. _sphx_glr_download_emd_tutorials_01_sifting_emd_tutorial_01_sift_05_iterated_masksift.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_05_iterated_masksift.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: emd_tutorial_01_sift_05_iterated_masksift.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_