.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "emd_tutorials/02_spectrum_analysis/emd_tutorial_02_spectrum_03_crossfrequency.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_02_spectrum_analysis_emd_tutorial_02_spectrum_03_crossfrequency.py: Cross-Frequency Coupling ======================== The spectrum tools in EMD can be used to explore cross-frequency interactions in oscillatory signals. The simplest approach is to look at how the phase, frequency or amplitude of a high-frequency signal interacts with the phase of a low frequency signal. This can be extended into 2 or 3 dimensions by exploring the relationship between low-frequency phase and the Hilbert-Huang Transform and Holospectrum. This tutorial shows some examples of these analyses with three signals with varying phase-amplitude coupling profiles. .. GENERATED FROM PYTHON SOURCE LINES 15-19 Simulating a signal with amplitude modulations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ First of all, we import EMD alongside numpy and matplotlib. We will also use scipy's ndimage module to smooth our results for visualisation later. .. GENERATED FROM PYTHON SOURCE LINES 19-26 .. code-block:: Python # sphinx_gallery_thumbnail_number = 8 import matplotlib.pyplot as plt import numpy as np import emd .. GENERATED FROM PYTHON SOURCE LINES 27-33 Next, we create three simulated signal to analyse. Each signal will be composed of two oscillations - a low-frequency signal at 5Hz and a high-frequency signal at 37Hz. The three variants of this signal differ in how the amplitude of the high frequency signal varies with the phase of the low frequency signal. We vary the width of the high-frequency burst across the three signals. .. GENERATED FROM PYTHON SOURCE LINES 33-55 .. code-block:: Python seconds = 60 sample_rate = 200 t = np.linspace(0, seconds, seconds*sample_rate) # First we create a slow 5Hz oscillation and a fast 37Hz oscillation slow = np.sin(2*np.pi*5*t - np.pi/2) fast = np.cos(2*np.pi*37*t) # Next we create three different amplitude modulation signals for the fast # oscillation. One sinusoidal, one wide modulation and one narrow modulation. # These cases differ by the duration of the high frequency burst. fast_am = 0.5*slow + 0.5 fast_am_narrow = fast_am**3 fast_am_wide = 1 - (0.5*-slow + 0.5)**3 # We create our signal by summing the oscillation and adding some noise x = slow + fast * fast_am + np.random.randn(*t.shape) * .05 x_wide = slow + fast * fast_am_wide + np.random.randn(*t.shape) * .05 x_narrow = slow + fast * fast_am_narrow + np.random.randn(*t.shape) * .05 .. GENERATED FROM PYTHON SOURCE LINES 56-57 Let's take a look at the three signals. .. GENERATED FROM PYTHON SOURCE LINES 57-66 .. code-block:: Python plt.figure() plt.subplot(311) plt.plot(t[:seconds*3], x_narrow[:seconds*3]) plt.subplot(312) plt.plot(t[:seconds*3], x[:seconds*3]) plt.subplot(313) plt.plot(t[:seconds*3], x_wide[:seconds*3]) .. image-sg:: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_001.png :alt: emd tutorial 02 spectrum 03 crossfrequency :srcset: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 67-71 The narrow amplitude modulated signal is plotted on the top, the sinusoidally modulated case in the middle and the wide amplitude modulation on the bottom. All three high-frequency signals peak at the same point in the low-frequency cycle. .. GENERATED FROM PYTHON SOURCE LINES 73-76 Next we run a mask sift on these signals to create IMFs for each case before running the frequency transforms to get the instantaneous phase, frequency and phase. .. GENERATED FROM PYTHON SOURCE LINES 76-94 .. code-block:: Python # Define a mask sift config config = emd.sift.get_config('mask_sift') config['max_imfs'] = 7 config['mask_freqs'] = 50/sample_rate config['mask_amp_mode'] = 'ratio_sig' config['imf_opts/sd_thresh'] = 0.05 config['extrema_opts/method'] = 'rilling' imf = emd.sift.mask_sift(x, **config) IP, IF, IA = emd.spectra.frequency_transform(imf, sample_rate, 'hilbert') imf_wide = emd.sift.mask_sift(x_wide, **config) IP_wide, IF_wide, IA_wide = emd.spectra.frequency_transform(imf_wide, sample_rate, 'hilbert') imf_narrow = emd.sift.mask_sift(x_narrow, **config) IP_narrow, IF_narrow, IA_narrow = emd.spectra.frequency_transform(imf_narrow, sample_rate, 'hilbert') .. GENERATED FROM PYTHON SOURCE LINES 95-98 We plot up the IMFs for the sinusoidally modulated case and see that the fast signal is isolated into the first IMF and the low-frequency signal is in the fourth. The rest of the IMFs largely contain noise .. GENERATED FROM PYTHON SOURCE LINES 98-101 .. code-block:: Python emd.plotting.plot_imfs(imf[:sample_rate*2, :]) .. image-sg:: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_002.png :alt: emd tutorial 02 spectrum 03 crossfrequency :srcset: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 102-104 1D Phase-amplitude coupling with instantaneous amplitude ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 106-109 Let's zoom in on the fast signals containing the simulated amplitude modulations. Here, we plot the first IMF for each of the three signals with different modulation widths. .. GENERATED FROM PYTHON SOURCE LINES 109-126 .. code-block:: Python plt.figure() plt.subplot(411) plt.plot(t[:seconds*3], imf_narrow[:seconds*3, 0]) plt.plot(t[:seconds*3], IA_narrow[:seconds*3, 0]) plt.subplot(412) plt.plot(t[:seconds*3], imf[:seconds*3, 0]) plt.plot(t[:seconds*3], IA[:seconds*3, 0]) plt.subplot(413) plt.plot(t[:seconds*3], imf_wide[:seconds*3, 0]) plt.plot(t[:seconds*3], IA_wide[:seconds*3, 0]) plt.subplot(414) plt.plot(t[:seconds*3], IA_narrow[:seconds*3, 0]) plt.plot(t[:seconds*3], IA[:seconds*3, 0]) plt.plot(t[:seconds*3], IA_wide[:seconds*3, 0]) plt.legend(['Narrow AM', 'Sinusoidal AM', 'Wide AM']) .. image-sg:: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_003.png :alt: emd tutorial 02 spectrum 03 crossfrequency :srcset: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 127-131 The narrow modulation is on top-row, the sinusoidal modulation in the second and the wide modulation is in the third. All three amplitude modulations on top of each other in the bottom row. By eye, the different modulation widths are clear but it is perhaps less obvious how we can quantify this effect. .. GENERATED FROM PYTHON SOURCE LINES 133-142 One simple approach is to explore how the instantaneous amplitude of each case varies with the phase of the low-frequency signal component. We can compute this with ``emd.cycles.bin_by_phase``. The inputs take a phase signal which is segmented into a set of time-bins in which the second input is averaged. (A very similar alternative analysis could be run using ``emd.cycles.phase_align``). We compute the average high-frequency IA as a function of low-frequency phase for each of our three examples and plot the result: .. GENERATED FROM PYTHON SOURCE LINES 142-155 .. code-block:: Python ia_by_phase, ia_by_phase_var, phase_bins = emd.cycles.bin_by_phase(IP[:, 4], IA[:, 0], nbins=24) ia_by_phase_wide, ia_by_phase_wide_var, _ = emd.cycles.bin_by_phase(IP[:, 4], IA_wide[:, 0], nbins=24) ia_by_phase_narrow, ia_by_phase_narrow_var, _ = emd.cycles.bin_by_phase(IP[:, 4], IA_narrow[:, 0], nbins=24) plt.figure() plt.errorbar(phase_bins, ia_by_phase, yerr=ia_by_phase_var) plt.errorbar(phase_bins, ia_by_phase_wide, yerr=ia_by_phase_wide_var) plt.errorbar(phase_bins, ia_by_phase_narrow, yerr=ia_by_phase_narrow_var) plt.legend(['Fast Amp', 'Fast Amp Wide', ' Fast Amp Narrow']) plt.xlabel('Slow Phase') plt.ylabel('IMF-1 IA') .. image-sg:: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_004.png :alt: emd tutorial 02 spectrum 03 crossfrequency :srcset: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(42.722222222222214, 0.5, 'IMF-1 IA') .. GENERATED FROM PYTHON SOURCE LINES 156-161 Now, we see the three amplitude modulation profiles directly linked to theta phase - there is a clear peak in high-frequency amplitude at one point in theta phase confirming the presence of phase-amplitude coupling in this signal. We can also see the different in modulation width in the three signals. .. GENERATED FROM PYTHON SOURCE LINES 164-166 2D Phase-amplitude coupling with the HHT ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 168-174 We can run a two-dimensional equivalent to the above analysis by exploring how a whole Hilbert-Huang transform varies across low-frequency phase. This gives us a bit more information than our first analysis. Specifically we can also see which frequency the amplitude modulated signal is peaking at. The first step is to compute the Hilbert-Huang Transform (HHT) for our three signals. .. GENERATED FROM PYTHON SOURCE LINES 174-182 .. code-block:: Python freq_edges, freq_centres = emd.spectra.define_hist_bins(10, 75, 75, 'log') f, hht = emd.spectra.hilberthuang(IF, IA, freq_edges, mode='amplitude', sum_time=False) f, hht_wide = emd.spectra.hilberthuang(IF_wide, IA_wide, freq_edges, mode='amplitude', sum_time=False) f, hht_narrow = emd.spectra.hilberthuang(IF_narrow, IA_narrow, freq_edges, mode='amplitude', sum_time=False) .. GENERATED FROM PYTHON SOURCE LINES 183-187 A quick summary figure shows us the HHT representation of our amplitude modulated signals. The recurring high-frequency bursts are visible in all three cases with the duration of each burst differing between the narrow, sinusoidal and widde cases. .. GENERATED FROM PYTHON SOURCE LINES 187-198 .. code-block:: Python plt.figure() plt.subplot(411) plt.plot(t[:seconds*3], x[:seconds*3]) plt.subplot(412) plt.pcolormesh(t[:seconds*3], f, hht_narrow[:, :seconds*3], cmap='hot_r', shading='auto') plt.subplot(413) plt.pcolormesh(t[:seconds*3], f, hht[:, :seconds*3], cmap='hot_r', shading='auto') plt.subplot(414) plt.pcolormesh(t[:seconds*3], f, hht_wide[:, :seconds*3], cmap='hot_r', shading='auto') .. image-sg:: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_005.png :alt: emd tutorial 02 spectrum 03 crossfrequency :srcset: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 199-203 We can use ``emd.cycles.bin_by_phase`` to explore how high dimensional signals vary with a phase time-course - as long as the first dimension of the signal to be averaged matches the length of the phase signal. Here, we bin the HHT by low-frequency phase for each of our three signals. .. GENERATED FROM PYTHON SOURCE LINES 203-208 .. code-block:: Python hht_by_phase, _, phase_centres = emd.cycles.bin_by_phase(IP[:, 4], hht.T) hht_by_phase_wide, _, _ = emd.cycles.bin_by_phase(IP[:, 4], hht_wide.T) hht_by_phase_narrow, _, _ = emd.cycles.bin_by_phase(IP[:, 4], hht_narrow.T) .. GENERATED FROM PYTHON SOURCE LINES 209-210 Let's make a quick summary figure of the phase-resolved HHTs. .. GENERATED FROM PYTHON SOURCE LINES 210-226 .. code-block:: Python plt.figure() plt.subplot(131) plt.pcolormesh(phase_centres, freq_centres, hht_by_phase_narrow.T, vmax=0.25, cmap='hot_r', shading='auto') plt.xlabel('Phase (rads') plt.ylabel('Frequency (Hz)') plt.title('Narrow AM') plt.subplot(132) plt.pcolormesh(phase_centres, freq_centres, hht_by_phase.T, vmax=0.25, cmap='hot_r', shading='auto') plt.xlabel('Phase (rads') plt.title('Sinusoidal AM') plt.subplot(133) plt.pcolormesh(phase_centres, freq_centres, hht_by_phase_wide.T, vmax=0.25, cmap='hot_r', shading='auto') plt.xlabel('Phase (rads') plt.title('Wide AM') .. image-sg:: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_006.png :alt: Narrow AM, Sinusoidal AM, Wide AM :srcset: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Wide AM') .. GENERATED FROM PYTHON SOURCE LINES 227-236 We can see that all three cases peak around the same point in phase (around pi/2) at at the same frequency (37Hz) but clearly differ in their modulation width. Though this is visible here, we still don't directly quantify this feature with the Hilbert-Huang Transform. We can quantify this difference by looking at the frequency content of the high-frequency signal amplitude envelopes. It appears like the signal with wide modulation has a lower frequency amplitude modulation How can we measure this? .. GENERATED FROM PYTHON SOURCE LINES 238-240 Quantifying the frequency of amplitude modulations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 242-252 The Holospectrum provides a method for direct quantification of the frequency and timing of amplitude modulations in oscillatory signal. The concept and method for Holospectrum estimation are covered in detail in the previous tutorial - here we will only cover its application in cross-frequency coupling. First, we compute the second-layer sift of our amplitude modulated signals. We apply a mask-sift at the second level to the instantaneous amplitude of our high frequency signals. The frequency transform of the resulting second-layer IMFs is then computed using the Hilbert transform. .. GENERATED FROM PYTHON SOURCE LINES 252-268 .. code-block:: Python # Reduce the mask frequencies for the second layer sift - otherwise we use the # same config as the first layer. config['mask_freqs'] = [10/sample_rate/ii for ii in range(1, 10)] imf2 = emd.sift.mask_sift_second_layer(IA[:, 0], config['mask_freqs'], sift_args=config) IP2, IF2, IA2 = emd.spectra.frequency_transform(imf2, sample_rate, 'hilbert') imf2_wide = emd.sift.mask_sift_second_layer(IA_wide[:, 0], config['mask_freqs'], sift_args=config) IP2_wide, IF2_wide, IA2_wide = emd.spectra.frequency_transform(imf2_wide, sample_rate, 'hilbert') imf2_narrow = emd.sift.mask_sift_second_layer(IA_narrow[:, 0], config['mask_freqs'], sift_args=config) IP2_narrow, IF2_narrow, IA2_narrow = emd.spectra.frequency_transform(imf2_narrow, sample_rate, 'hilbert') .. GENERATED FROM PYTHON SOURCE LINES 269-272 Next we compute the Holospectrum itself. We define a two sets of frequency bins for the first layer oscillations and second layer amplitude modulations before computing the Holospectrum using ``emd.spectra.holospectrum``. .. GENERATED FROM PYTHON SOURCE LINES 272-283 .. code-block:: Python freq_edges, freq_centres = emd.spectra.define_hist_bins(1, 75, 75, 'log') freq_edges_low, freq_centres_low = emd.spectra.define_hist_bins(2, 10, 20, 'linear') fcarrier, fam, holo = emd.spectra.holospectrum(IF[:, 0, None], IF2, IA2, freq_edges, freq_edges_low) fcarrier, fam, holo_narrow = emd.spectra.holospectrum(IF_narrow[:, 0, None], IF2_narrow, IA2_narrow, freq_edges, freq_edges_low) fcarrier, fam, holo_wide = emd.spectra.holospectrum(IF_wide[:, 0, None], IF2_wide, IA2_wide, freq_edges, freq_edges_low) # Holo is [second-layer frequencies x first-layer frequencies] print(holo.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none (75, 20) .. GENERATED FROM PYTHON SOURCE LINES 284-290 The output holospectrum is a 2d matrix of first-layer frequencies by second-layer frequencies (we sum across time dimension by default). We next make a summary plot with this 2d Holospectrum for each case in the top row and the average holospectrum across first-layer frequencies (ie expressing the energy in signals with given amplitude modulation frequencies across all first-layer frequencies.) .. GENERATED FROM PYTHON SOURCE LINES 290-321 .. code-block:: Python plt.figure() plt.subplot(231) plt.pcolormesh(freq_centres_low, freq_centres, holo_narrow, cmap='hot_r', shading='auto') plt.ylabel('Frequency (Hz)') plt.xticks(np.arange(3, 10, 2)) plt.title('Narrow AM') plt.subplot(234) plt.plot(freq_centres_low, holo_narrow[62, :]) plt.xlabel('AM Frequency (Hz)') plt.xticks(np.arange(3, 10, 2)) plt.ylabel('Amplitude') plt.subplot(232) plt.pcolormesh(freq_centres_low, freq_centres, holo, cmap='hot_r', shading='auto') plt.xticks(np.arange(3, 10, 2)) plt.title('Sinusoidal AM') plt.subplot(235) plt.plot(freq_centres_low, holo[62, :]) plt.xticks(np.arange(3, 10, 2)) plt.xlabel('AM Frequency (Hz)') plt.subplot(233) plt.pcolormesh(freq_centres_low, freq_centres, holo_wide, cmap='hot_r', shading='auto') plt.xticks(np.arange(3, 10, 2)) plt.title('Wide AM') plt.subplot(236) plt.plot(freq_centres_low, holo_wide[62, :]) plt.xticks(np.arange(3, 10, 2)) plt.xlabel('AM Frequency (Hz)') .. image-sg:: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_007.png :alt: Narrow AM, Sinusoidal AM, Wide AM :srcset: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 23.52222222222222, 'AM Frequency (Hz)') .. GENERATED FROM PYTHON SOURCE LINES 322-342 The holospectra on the top row show the distribution of energy across frequency and amplitude modulation frequency within each signal. The bottom row sums the holospectra across the y-axis to summarise just the distribution of energy across apmlitude modulation frequencies. The sinusoidal signal has a clear peak with a frequency of 37Hz with amplitude modulations of 5Hz - exactly as we would expect from this simulation. The other signals have a similar peak but slightly skewed to higher or lower amplitude modulation frequencies. The signal with narrow modulation has higher amplitude modulation frequencies - reflecting the faster/sharper amplitude modulation profile of the signal. In contrast, the signal with wide amplitude modulations skews towards slower amplitude modulations - reflectings its slower, flatter amplitude modulation profile. The holospectrum provides a convenient summary of the amplitude modulations in a signal, but doesn't explicitly link them to the phase of a lower frequency signal. To complete a full cross-frequency coupling analysis with the Holospectrum we need to show not only that our high frequency signal has amplitude modulations but that those amplitude modulations are specifically linked to a n observed low-frequency signal. .. GENERATED FROM PYTHON SOURCE LINES 345-347 3D phase-amplitude coupling with the Holospectrum ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 349-356 We can link the Holospectrum to low frequency phase in exactly the same way as we analysed the instantaneous amplitude and Hilbert-Huang Transforms. First, we have to recompute the holospectra whilst preserving the time dimension in the output. By default, the holospectrum is summed over time before being returned - we can ask for the full 3D holospectrum to be returned by setting ``squash_time=False`` in the ``emd.spectra.holospectrum`` call. .. GENERATED FROM PYTHON SOURCE LINES 356-371 .. code-block:: Python # Define a new set of second-layer frequencies - slightly wider than the last one. freq_edges_low, freq_centres_low = emd.spectra.define_hist_bins(.5, 15, 32, 'linear') fcarrier, fam, holot = emd.spectra.holospectrum(IF[:, 0, None], IF2, IA2, freq_edges, freq_edges_low, sum_time=False) fcarrier, fam, holot_wide = emd.spectra.holospectrum(IF_wide[:, 0, None], IF2_wide, IA2_wide, freq_edges, freq_edges_low, sum_time=False) fcarrier, fam, holot_narrow = emd.spectra.holospectrum(IF_narrow[:, 0, None], IF2_narrow, IA2_narrow, freq_edges, freq_edges_low, sum_time=False) # holot is [second-layer frequencies x first-layer frequencies x time samples] print(holot.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none (75, 32, 12000) .. GENERATED FROM PYTHON SOURCE LINES 372-375 Now that we have a time-varying holospectrum estimate - we can compute its average in different low-frequency bins by using ``emd.cycles.bin_by_phase`` as for our previous analyses. .. GENERATED FROM PYTHON SOURCE LINES 375-383 .. code-block:: Python holo_by_phase, _, phase_bins = emd.cycles.bin_by_phase(IP[:, 4], np.swapaxes(holot, 0, 2)) holo_by_phase_wide, _, _ = emd.cycles.bin_by_phase(IP[:, 4], np.swapaxes(holot_wide, 0, 2)) holo_by_phase_narrow, _, _ = emd.cycles.bin_by_phase(IP[:, 4], np.swapaxes(holot_narrow, 0, 2)) # holo_by_phase is [phase bins x second-layer frequencies x first-layer frequencies] print(holo_by_phase.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none (24, 32, 75) .. GENERATED FROM PYTHON SOURCE LINES 384-395 ``holo_by_phase`` now contains a summary of the energy in each signal separated by first-layer frequencies, second layer-frequencies and by the phase of our low frequency signal. If the amplitude modulations we observed in our signals are explicitly linked to the 5Hz phase - we should be able to see it here. We can summarise ``holo_by_phase`` a few ways. To start we we plot the high-frequency instantaneous amplitude on the top-row for reference. The second row contains the first-layer frequencies in the Holospectrum as a function of 5Hz phase, finally the third row contains the second-layer amplitude modulation frequencies across 5Hz phase. .. GENERATED FROM PYTHON SOURCE LINES 395-430 .. code-block:: Python plt.figure(figsize=(8, 8)) plt.subplot(331) plt.errorbar(phase_bins, ia_by_phase_narrow, yerr=ia_by_phase_narrow_var) plt.ylim(0, 1) plt.ylabel('High-frequency\ninstantaneous amplitudee') plt.title('Narrow AM') plt.subplot(332) plt.errorbar(phase_bins, ia_by_phase, yerr=ia_by_phase_var) plt.ylim(0, 1) plt.title('Sinusoidal AM') plt.subplot(333) plt.errorbar(phase_bins, ia_by_phase_wide, yerr=ia_by_phase_wide_var) plt.ylim(0, 1) plt.title('Wide AM') plt.subplot(334) plt.pcolormesh(phase_centres, freq_centres, holo_by_phase_narrow.sum(axis=1).T, cmap='hot_r', shading='auto') plt.ylabel('Frequency (Hz)') plt.subplot(335) plt.pcolormesh(phase_centres, freq_centres, holo_by_phase.sum(axis=1).T, cmap='hot_r', shading='auto') plt.subplot(336) plt.pcolormesh(phase_centres, freq_centres, holo_by_phase_wide.sum(axis=1).T, cmap='hot_r', shading='auto') plt.subplot(337) plt.pcolormesh(phase_centres, freq_centres_low, holo_by_phase_narrow.sum(axis=2).T, cmap='hot_r', shading='auto') plt.ylabel('Amplitude Modulation\nFrequency (Hz)') plt.xlabel('Slow phase') plt.subplot(338) plt.pcolormesh(phase_centres, freq_centres_low, holo_by_phase.sum(axis=2).T, cmap='hot_r', shading='auto') plt.xlabel('Slow phase') plt.subplot(339) plt.pcolormesh(phase_centres, freq_centres_low, holo_by_phase_wide.sum(axis=2).T, cmap='hot_r', shading='auto') plt.xlabel('Slow phase') .. image-sg:: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_008.png :alt: Narrow AM, Sinusoidal AM, Wide AM :srcset: /emd_tutorials/02_spectrum_analysis/images/sphx_glr_emd_tutorial_02_spectrum_03_crossfrequency_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 58.7222222222222, 'Slow phase') .. GENERATED FROM PYTHON SOURCE LINES 431-450 The second row here is very similar to the HHT-by-phase plot in our previous section and confirms that there is a peak in 37Hz power around pi/2 in the low-frequency phase. Again, the width of the modulation is visible as stretching of the peak in the x-axis. The third row gives us some new information about how amplitude modulation frequency of the 37Hz frequency signal changes across the phase of the 5Hz signal. This is a flat profile at 5Hz for the sinsoidal signal but changes for the narrow and wide case. Critically the narrow amplitude modulation has an increase in amplitude modulation frequency around pi/2 - the point in the low-frequency signal where the high-frequency signal peaks. In contrast, the wide modulation signal has a lower amplitude modulation frequency at the same point. This reflects the fast-sharp peak in the narrow modualtion and the flat peak in the wide modulation. All three signals are clearly linked to the low-frequency phase. Interestingly, the holospectrum is able to quantify the non-linear differences in amplitude modulation frequency driven by the diferences in modulation width. .. GENERATED FROM PYTHON SOURCE LINES 452-454 Further Reading & References ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 456-468 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 Huang, N. E., Hu, K., Yang, A. C. C., Chang, H.-C., Jia, D., Liang, W.-K., … Wu, Z. (2016). On Holo-Hilbert spectral analysis: a full informational spectral representation for nonlinear and non-stationary data. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374(2065), 20150206. https://doi.org/10.1098/rsta.2015.0206 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.129 seconds) .. _sphx_glr_download_emd_tutorials_02_spectrum_analysis_emd_tutorial_02_spectrum_03_crossfrequency.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_02_spectrum_03_crossfrequency.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: emd_tutorial_02_spectrum_03_crossfrequency.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_