.. 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_04_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_04_masksift.py: Masked sifting ============== This tutorial introduces some of the issues that standard EMD algorithms can have with intermitent signals and shows how the Masked sift can resolve them. .. GENERATED FROM PYTHON SOURCE LINES 10-11 Lets make a simulated signal to get started. .. GENERATED FROM PYTHON SOURCE LINES 11-49 .. code-block:: Python import emd import numpy as np import matplotlib.pyplot as plt seconds = 5 sample_rate = 1024 time_vect = np.linspace(0, seconds, seconds*sample_rate) # Create an amplitude modulation am = np.sin(2*np.pi*time_vect) am[am < 0] = 0 # Create a 25Hz signal and introduce the amplitude modulation xx = am*np.sin(2*np.pi*25*time_vect) # Create a non-modulated 6Hz signal yy = .5*np.sin(2*np.pi*6*time_vect) # Sum the 25Hz and 6Hz components together xy = xx+yy # Make a quick summary plot plt.figure(figsize=(6,9)) plt.subplots_adjust(hspace=0.3) plt.subplot(311) plt.plot(xy) plt.title('Full Signal') plt.subplot(312) plt.plot(yy) plt.title('Slow Signal') plt.subplot(313) plt.plot(xx) plt.title('Transient Fast Signal') plt.xlabel('Samples') # sphinx_gallery_thumbnail_number = 2 .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_001.png :alt: Full Signal, Slow Signal, Transient Fast Signal :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 69.7222222222222, 'Samples') .. GENERATED FROM PYTHON SOURCE LINES 50-56 This signal doesn't contain any noise and only has two frequency components so should be straightforward to sift. Unfortunately, as the 25Hz signal component disappears completely for parts of the signal the EMD doesn't quite do what we'd want it to. Here we run a default sift and plot the IMFs. .. GENERATED FROM PYTHON SOURCE LINES 56-60 .. code-block:: Python imf = emd.sift.sift(xy, max_imfs=3) emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_002.png :alt: emd tutorial 01 sift 04 masksift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 61-79 The signals are well separated when both oscillations are present. However in time periods where the fast 25Hz signal disappears the slower signal jumps up to become part of the fast component. We'd prefer the separation into narrow band components as seen in the simulations above... This happens as EMD is a locally adaptive algorithm - the peaks and troughs in the signal define the time-scales that are analysed for a given part of the signal. So, the first IMF will always find the fastest peaks for every part of the signal even if the definition of 'fast' might be different in different segments. The masked sift is a potential solution to this problem. This is a simple trick which effectively puts a lower bound on the frequency content that can enter a particular IMF. We will add a known masking signal to our time-series before running ``emd.sift.get_next_imf``. Any signals which are lower in frequency than this mask should then be ignored by the sift in favour of this known signal. Finally, we can remove the known mask to recover our IMF. .. GENERATED FROM PYTHON SOURCE LINES 82-84 Adding a single mask ^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 86-87 Here we make a 30Hz mask and plot it next to a segment of our time-series. .. GENERATED FROM PYTHON SOURCE LINES 87-102 .. code-block:: Python mask = 0.1*np.sin(2*np.pi*30*time_vect) plt.figure(figsize=(10, 6)) plt.subplot(121) plt.plot(xy) plt.plot(mask) plt.legend(['Signal', 'Mask']) plt.xlim(0, 1000) plt.subplot(122) plt.plot(xy + mask, color='k') plt.xlim(0, 1000) plt.legend(['Signal + Mask']) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_003.png :alt: emd tutorial 01 sift 04 masksift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 103-107 We see that the masking signal is close in frequency to the fast burst but much faster than the 6Hz signal. Next we identify our next IMF on the raw signal with and without the mask .. GENERATED FROM PYTHON SOURCE LINES 107-111 .. code-block:: Python imf_raw, _ = emd.sift.get_next_imf(xy) imf_mask, _ = emd.sift.get_next_imf(xy+mask) .. GENERATED FROM PYTHON SOURCE LINES 112-117 The normal IMF in the top panel has the problem we saw earlier, the slow signal is leaking into the fast IMF. The masked IMF successfully suppresses this slow signal, replacing it with the mask frequency. Finally, subtracting the mask removes everything but the 25Hz oscillation which now correctly disappears between bursts. .. GENERATED FROM PYTHON SOURCE LINES 117-135 .. code-block:: Python plt.figure(figsize=(6,9)) plt.subplots_adjust(hspace=0.3) plt.subplot(311) plt.plot(imf_raw) plt.xlim(0, 1000) plt.title('Normal IMF') plt.gca().set_xticklabels([]) plt.subplot(312) plt.plot(imf_mask) plt.xlim(0, 1000) plt.title('Masked IMF') plt.gca().set_xticklabels([]) plt.subplot(313) plt.plot(imf_mask - mask[:, np.newaxis]) plt.xlim(0, 1000) plt.title('Masked IMF with Mask removed') .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_004.png :alt: Normal IMF, Masked IMF, Masked IMF with Mask removed :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Masked IMF with Mask removed') .. GENERATED FROM PYTHON SOURCE LINES 136-137 This effect is more obvious if we look at the whole time-courses without zooming in .. GENERATED FROM PYTHON SOURCE LINES 137-152 .. code-block:: Python plt.figure(figsize=(6, 9)) plt.subplots_adjust(hspace=0.3) plt.subplot(311) plt.plot(imf_raw) plt.title('Normal IMF') plt.gca().set_xticklabels([]) plt.subplot(312) plt.plot(imf_mask) plt.title('Masked IMF') plt.gca().set_xticklabels([]) plt.subplot(313) plt.plot(imf_mask - mask[:, np.newaxis]) plt.title('Masked IMF with Mask removed') .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_005.png :alt: Normal IMF, Masked IMF, Masked IMF with Mask removed :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Masked IMF with Mask removed') .. GENERATED FROM PYTHON SOURCE LINES 153-158 This masking process is implemented in ``emd.sift.get_next_imf_mask`` which works much like ``emd.sift.get_next_imf`` with a couple of extra options for adding masks. We can specify the frequency and amplitude of the mask to be applied whilst isolating our IMF. .. GENERATED FROM PYTHON SOURCE LINES 160-162 Choosing mask frequency ^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 164-171 It is important that the mask frequency is approximately equal to the signal component we want to isolate. If we use a mask of too high or too low frequency then the procedure will not work. Here we apply a range of masks to illustrate this effect. We apply a high frequency mask, a low frequency mask and a 'just-right' mask close to the 30Hz signal. .. GENERATED FROM PYTHON SOURCE LINES 171-197 .. code-block:: Python # Masks should be specified in normalised frequencies between 0 and .5 where 0.5 is half the sampling rate high_mask_freq = 150/sample_rate imf_high_mask, _ = emd.sift.get_next_imf_mask(xy, high_mask_freq, 2) med_mask_freq = 30/sample_rate imf_med_mask, _ = emd.sift.get_next_imf_mask(xy, med_mask_freq, 2) low_mask_freq = 2/sample_rate imf_low_mask, _ = emd.sift.get_next_imf_mask(xy, low_mask_freq, 2) plt.figure(figsize=(6, 9)) plt.subplots_adjust(hspace=0.3) plt.subplot(311) plt.plot(imf_high_mask) plt.title('150Hz mask - too high') plt.ylim(-1, 1) plt.subplot(312) plt.plot(imf_med_mask) plt.ylim(-1, 1) plt.title('30Hz mask - just right') plt.subplot(313) plt.plot(imf_low_mask) plt.ylim(-1, 1) plt.title('2Hz mask - too low') .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_006.png :alt: 150Hz mask - too high, 30Hz mask - just right, 2Hz mask - too low :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, '2Hz mask - too low') .. GENERATED FROM PYTHON SOURCE LINES 198-202 We can see that the high frequency mask has suppressed both signal components whist the low frequency mask has allowed them both through. In contrast, when the mask frequency is closest to the fastest dynamics present in the signal - it is able to successfully isolate them. .. GENERATED FROM PYTHON SOURCE LINES 204-219 Choosing a mask frequency may not be so simple in real data applications where multiple signal components may be present alongside measurement noise. We recommend a rule-of-thumb based on simulations in Fosso & Molinas 2017 (http://arxiv.org/abs/1709.05547). A given mask with frequency f Hz, can be expected to remove frequencies below 0.7*f Hz. As an example, a mask with frequency 50Hz will keep oscillatory content of 40Hz but remove content of 30Hz. The flexibility provided by masked sifting overcomes several issue with conventional sifting but this flexibility can also lead to problems. If you do use a highly tuned set of masks for an analysis, we recommend careful checking of the IMFs before proceeding to the main analysis. You should ensure that the IMFs are well separated - the components should have low correlations/orthogonality scores - and that the instantaneous frequency content of the IMFs do not strongly overlap. .. GENERATED FROM PYTHON SOURCE LINES 221-223 Running a full mask sift ^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 225-233 Until now, we've only been separating two signals with a single pass of the sift. ``emd.sift.mask_sift`` uses ``emd.sift.get_next_imf_mask`` internally to run a whole set of sifts using the masking method. Each IMF is isolated with a separate mask which decreases in frequency for each successive IMF. Here we run a ``mask_sift`` using mask frequencies starting at 30Hz. This will reduce by one half for each successive IMF - the second mask will be 15Hz, the third is 7.5Hz and so on. .. GENERATED FROM PYTHON SOURCE LINES 233-238 .. code-block:: Python imf, mask_freqs = emd.sift.mask_sift(xy, mask_freqs=30/sample_rate, ret_mask_freq=True, max_imfs=4) print(mask_freqs * sample_rate) .. rst-class:: sphx-glr-script-out .. code-block:: none [30. 15. 7.5 3.75] .. GENERATED FROM PYTHON SOURCE LINES 239-245 We can see that this sift nicely separates the two components. The first IMF contains the 25Hz bursting signal which returns to a flat line between events. The second IMF contains very low amplitude noise. This is as the mask frequency of 15Hz for the second mask is still too high to isolate the oscillation of 6Hz - so IMF 2 is essentially flat. The third IMF with a mask frequency of 7.5Hz is about right to isolate the 6Hz signal. .. GENERATED FROM PYTHON SOURCE LINES 245-248 .. code-block:: Python emd.plotting.plot_imfs(imf) .. image-sg:: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_007.png :alt: emd tutorial 01 sift 04 masksift :srcset: /emd_tutorials/01_sifting/images/sphx_glr_emd_tutorial_01_sift_04_masksift_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 249-251 Further Reading & References ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 253-262 Deering R, Kaiser JF (2005) The use of a masking signal to improve empirical mode decomposition. IEEE International Conference on Acoustics, Speech, and Signal Processing. vol. 4. p. iv/485–iv/488. https://doi.org/10.1109/ICASSP.2005.1416051. Andrew J. Quinn, Vítor Lopes-dos-Santos, Norden Huang, Wei-Kuang Liang, Chi-Hung Juan, Jia-Rong Yeh, Anna C. Nobre, David Dupret, and Mark W. Woolrich (2001) Within-cycle instantaneous frequency profiles report oscillatory waveform dynamics Journal of Neurophysiology 126:4, 1190-1208 https://doi.org/10.1152/jn.00201.2021 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.583 seconds) .. _sphx_glr_download_emd_tutorials_01_sifting_emd_tutorial_01_sift_04_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_04_masksift.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: emd_tutorial_01_sift_04_masksift.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_