.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "emd_tutorials/03_cycle_ananlysis/emd_tutorial_03_cycle_02_statistics.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_03_cycle_ananlysis_emd_tutorial_03_cycle_02_statistics.py: Cycle statistics and comparisons ================================ Here we will use the 'cycle' submodule of EMD to identify and analyse individual cycles of an oscillatory signal .. GENERATED FROM PYTHON SOURCE LINES 9-12 Simulating a noisy signal ^^^^^^^^^^^^^^^^^^^^^^^^^ Firstly we will import emd and simulate a signal .. GENERATED FROM PYTHON SOURCE LINES 12-34 .. code-block:: Python import emd import numpy as np import matplotlib.pyplot as plt from scipy import ndimage # Define and simulate a simple signal peak_freq = 15 sample_rate = 256 seconds = 10 noise_std = .4 x = emd.simulate.ar_oscillator(peak_freq, sample_rate, seconds, noise_std=noise_std, random_seed=42, r=.96)[:, 0] t = np.linspace(0, seconds, seconds*sample_rate) # Plot the first 5 seconds of data plt.figure(figsize=(10, 2)) plt.plot(t[:sample_rate*4], x[:sample_rate*4], 'k') # sphinx_gallery_thumbnail_number = 5 .. image-sg:: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_001.png :alt: emd tutorial 03 cycle 02 statistics :srcset: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 35-41 Extract IMFs & find cycles ^^^^^^^^^^^^^^^^^^^^^^^^^^ We next run a mask sift with the default parameters to isolate the 15Hz oscillation. There is only one clear oscillatory signal in this simulation. This is extracted in IMF-3 whilst the remaining IMFs contain low-amplitude noise. .. GENERATED FROM PYTHON SOURCE LINES 41-49 .. code-block:: Python # Run a mask sift imf = emd.sift.mask_sift(x, max_imfs=5) # Visualise the IMFs emd.plotting.plot_imfs(imf[:sample_rate*4, :]) .. image-sg:: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_002.png :alt: emd tutorial 03 cycle 02 statistics :srcset: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 50-53 Next we locate the cycle indices from the instantaneous phase of our IMFs. We do this twice, once to identify all cycles and a second time to identify only 'good' cycles based on the cycle validation check from the previous tutorial. .. GENERATED FROM PYTHON SOURCE LINES 53-61 .. code-block:: Python # Extract frequency information IP, IF, IA = emd.spectra.frequency_transform(imf, sample_rate, 'nht') # Extract cycle locations all_cycles = emd.cycles.get_cycle_vector(IP, return_good=False) good_cycles = emd.cycles.get_cycle_vector(IP, return_good=True) .. GENERATED FROM PYTHON SOURCE LINES 62-76 We can customise the parts of the signal in which we look for cycles by defining a mask. This is a binary vector indicating which samples in a time-series should be included in the cycle detection. This could be useful for several reasons, we can mask our sections of signal with artefacts, limit cycle detection to a specific period during a task or limit cycle detection to periods where there is a high amplitude oscillation. Here we will apply a low amplitude threshold to identify good cycles which have amplitude values strictly above the 33th percentile of amplitude values in the dataset - excluding the lowest amplitude cycles. Note that the whole cycle must be in the valid part of the mask to be included, a cycle will be excluded if a single sample within it is masked out. .. GENERATED FROM PYTHON SOURCE LINES 76-82 .. code-block:: Python thresh = np.percentile(IA[:, 2], 33) mask = IA[:, 2] > thresh mask_cycles = emd.cycles.get_cycle_vector(IP, return_good=True, mask=mask) .. GENERATED FROM PYTHON SOURCE LINES 83-105 We can compute a variety of metric from our cycles using the ``emd.cycles.get_cycle_stat`` function. This is a simple helper function which takes in a set of cycle timings (the output from ``emd.cycles.get_cycle_vector``) and any time-series of interest (such as instaneous amplitude or frequency). The function then computes a metric from the time-series within each cycle. The computed metric is defined by the ``func`` argument, this can be any function which takes a vector input and returns a single-number. Often we will use se the numpy built-in functions to compute simple metrics (such as ``np.max`` or ``np.mean``) but we can use a custom user-defined function as well. Finally we can define whether to return the result in ``full`` or ``compressed`` format. The full form returns a vector of the same length as the input vector in which the indices for each cycle contains the its cycle-stat whilst, the compressed form returns a vector containing single values for each cycle in turn. For instance, the following example computes the maximum instantaneous amplitude for all detected cycles in IMF-3 and returns the result in the full-vector format. .. GENERATED FROM PYTHON SOURCE LINES 105-115 .. code-block:: Python cycle_amp = emd.cycles.get_cycle_stat(all_cycles[:, 2], IA[:, 2], out='samples', func=np.max) # Make a summary figure plt.figure(figsize=(10, 4)) plt.plot(t[:sample_rate*4], imf[:sample_rate*4, 2], 'k') plt.plot(t[:sample_rate*4], IA[:sample_rate*4, 2], 'b') plt.plot(t[:sample_rate*4], cycle_amp[:sample_rate*4], 'r') plt.legend(['IMF-3', 'Instantaneous Amplitude', 'Cycle-max Amplitude']) .. image-sg:: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_003.png :alt: emd tutorial 03 cycle 02 statistics :srcset: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 116-123 We can see that the original IMF in black and its instantaneous amplitude in blue. The red line is then the full-format output containing the cycle maximum amplitude. This nicely corresponds to the peak amplitude for each cycle as seen in blue. The next section computes the average instantaneous frequency within each cycle, again returning the result in full format. .. GENERATED FROM PYTHON SOURCE LINES 123-132 .. code-block:: Python cycle_freq = emd.cycles.get_cycle_stat(all_cycles[:, 2], IF[:, 2], out='samples', func=np.mean) # Make a summary figure plt.figure(figsize=(10, 4)) plt.plot(t[:sample_rate*4], IF[:sample_rate*4, 2], 'b') plt.plot(t[:sample_rate*4], cycle_freq[:sample_rate*4], 'r') plt.legend(['Instantaneous Frequency', 'Cycle-mean frequency']) .. image-sg:: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_004.png :alt: emd tutorial 03 cycle 02 statistics :srcset: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 133-138 We can get a nice visualisation of cycle-average frequency by overlaying the full stat vector onto the Hilbert-Huang transform. This is similar to the plot above but now we can see the signal amplitude values in the colour-scale of the HHT (hotter colours show higher amplitudes). Here we plot the cycle-average frequency for cycles above our amplitude thresholdover the HHT .. GENERATED FROM PYTHON SOURCE LINES 138-166 .. code-block:: Python # Compute cycle freq using amplitude masked-cycle indices cycle_freq = emd.cycles.get_cycle_stat(mask_cycles[:, 2], IF[:, 2], out='samples', func=np.mean) # Carrier frequency histogram definition freq_range = (3, 25, 64) # Compute the 2d Hilbert-Huang transform (power over time x carrier frequency) f, hht = emd.spectra.hilberthuang(IF, IA, freq_range, mode='amplitude', sum_time=False) # Add a little smoothing to help visualisation shht = ndimage.gaussian_filter(hht, 1) # Make a summary plot plt.figure(figsize=(10, 7)) plt.subplots_adjust(hspace=.3) plt.subplot(211) plt.plot(t[:sample_rate*4], imf[:sample_rate*4, 2], 'k') plt.plot((0, 4), (thresh, thresh), 'k:') plt.xlim(0, 4) plt.title('IMF') plt.subplot(212) plt.pcolormesh(t[:sample_rate*4], f, shht[:, :sample_rate*4], cmap='hot_r', vmin=0) plt.plot(t[:sample_rate*4], cycle_freq[:sample_rate*4], 'k') plt.title('Hilbert-Huang Transform') plt.xlabel('Time (seconds)') plt.ylabel('Frequency (Hz)') .. image-sg:: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_005.png :alt: IMF, Hilbert-Huang Transform :srcset: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(92.09722222222221, 0.5, 'Frequency (Hz)') .. GENERATED FROM PYTHON SOURCE LINES 167-180 Compressed cycle stats ^^^^^^^^^^^^^^^^^^^^^^^^^^ The full-format output is useful for visualisation and validation, but often we only want to deal with a single number summarising each cycle. The compressed format provides this simplified output. Note that the first value of the compressed format contains the average for missing cycles in the analysis (where the value in the cycles vector equals zero) we will discard this for the following analyses as we are focusing on the properties of well formed oscillatory cycles. For a first example, we compute the average frequency and amplitude of all cycles. We then make a scatter plot to explore any relationship between amplitude and frequency. .. GENERATED FROM PYTHON SOURCE LINES 180-198 .. code-block:: Python # Compute cycle average frequency for all cycles and masked cycles all_cycle_freq = emd.cycles.get_cycle_stat(all_cycles[:, 2], IF[:, 2], func=np.mean) mask_cycle_freq = emd.cycles.get_cycle_stat(mask_cycles[:, 2], IF[:, 2], func=np.mean) # Compute cycle frequency range for all cycles and for masked cycles all_cycle_amp = emd.cycles.get_cycle_stat(all_cycles[:, 2], IA[:, 2], func=np.mean) mask_cycle_amp = emd.cycles.get_cycle_stat(mask_cycles[:, 2], IA[:, 2], func=np.mean) # Make a summary figures plt.figure() plt.plot(all_cycle_freq, all_cycle_amp, 'o') plt.plot(mask_cycle_freq, mask_cycle_amp, 'o') plt.xlabel('Cycle average frequency (Hz)') plt.ylabel('Cycle average amplitude') plt.plot((9, 22), (thresh, thresh), 'k:') plt.legend(['All-cycles', 'Masked-cycles', 'Amp thresh']) .. image-sg:: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_006.png :alt: emd tutorial 03 cycle 02 statistics :srcset: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 199-215 We see that high amplitude cycles are closely clustered around 15Hz - the peak frequency of our simulated oscillation. Lower amplitude cycles are noisier and have a wider frequency distribution. The rejected bad-cycles tend to have low amplitudes and come from a wide frequency distribution. A small number of cycles pass the amplitude threshold but are rejected by the cycle quality checks. These cycles may have phase distortions or other artefacts which have lead to ``emd.cycles.get_cycle_vector`` to remove them from the set of `good` cycles. We can include more complex user-defined functions to generate cycle stats. Here we compute a range of cycle stats in compressed format (discarding the first value in the output). We compute the cycle average frequency and cycle-max amplitude for all cycles and again for only the good cycles. We can then make a scatter plot to explore any relationship between amplitude and frequency. .. GENERATED FROM PYTHON SOURCE LINES 218-229 We can include more complicated metrics in user-specified functions. Here we compute the Degree of Non-linearity (DoN; https://doi.org/10.1371/journal.pone.0168108) of each cycle as an indication of the extent to which a cycle contains non-sinudoisal content. Note that the original DoN uses the zero-crossing frequency rather than mean frequency as a normalising factor. These factors are highly correlated so, for simplicity, we use the mean here. Here we compute the degree of non-linearity for all cycles and good cycles separately and plot the results as a function of cycle average frequency .. GENERATED FROM PYTHON SOURCE LINES 229-254 .. code-block:: Python # Compute cycle average frequency for all cycles and masked cycles all_cycle_freq = emd.cycles.get_cycle_stat(all_cycles[:, 2], IF[:, 2], func=np.mean) mask_cycle_freq = emd.cycles.get_cycle_stat(mask_cycles[:, 2], IF[:, 2], func=np.mean) # Define a simple function to compute the range of a set of values def degree_nonlinearity(x): return np.std((x - x.mean()) / x.mean()) # Compute cycle freuquency range for all cycles and for masked cycles all_cycle_freq_don = emd.cycles.get_cycle_stat(all_cycles[:, 2], IF[:, 2], func=degree_nonlinearity) cycle_freq_don = emd.cycles.get_cycle_stat(mask_cycles[:, 2], IF[:, 2], func=degree_nonlinearity) # Make a summary figures plt.figure() plt.plot(all_cycle_freq, all_cycle_freq_don, 'o') plt.plot(mask_cycle_freq, cycle_freq_don, 'o') plt.xlabel('Cycle average frequency (Hz)') plt.ylabel('Cycle IF don (Hz)') plt.legend(['All-cycles', 'Masked-cycles']) .. image-sg:: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_007.png :alt: emd tutorial 03 cycle 02 statistics :srcset: /emd_tutorials/03_cycle_ananlysis/images/sphx_glr_emd_tutorial_03_cycle_02_statistics_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 255-260 The majority of cycles with very high degree of non-linearity in this simulation have been rejected by either the amplitude threshold or the cycle quality checks. The surviving cycles (in orange) are tightly clustered around 15Hz peak frequency with a relatively low degree of non-linearity. We have not defined any non-linearity in this simulation. .. GENERATED FROM PYTHON SOURCE LINES 262-264 Further Reading & References ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 266-270 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.274 seconds) .. _sphx_glr_download_emd_tutorials_03_cycle_ananlysis_emd_tutorial_03_cycle_02_statistics.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_03_cycle_02_statistics.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: emd_tutorial_03_cycle_02_statistics.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_