.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_compute_waveshape.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_auto_examples_plot_compute_waveshape.py: ========================== Compute waveshape features ========================== This example demonstrates how waveshape features can be computed with PyBispectra. .. GENERATED FROM PYTHON SOURCE LINES 9-13 .. code-block:: Python # Author(s): # Thomas S. Binns | github.com/tsbinns .. GENERATED FROM PYTHON SOURCE LINES 14-21 .. code-block:: Python import numpy as np from numpy.random import RandomState from matplotlib import pyplot as plt from pybispectra import compute_fft, get_example_data_paths, WaveShape .. GENERATED FROM PYTHON SOURCE LINES 22-75 Background ---------- When analysing signals, important information may be gleaned from a variety of features, including their shape. For example, in neuroscience it has been suggested that non-sinusoidal waves may play important roles in physiology and pathology, such as waveform sharpness reflecting synchrony of synaptic inputs :footcite:`Sherman2016` and correlating with symptoms of Parkinson's disease :footcite:`Cole2017`. Two aspects of waveshape described in recent literature include: rise-decay asymmetry - how much the wave resembles a sawtooth pattern (also called waveform sharpness or derivative skewness); and peak-trough asymmetry - whether peaks (events with a positive-valued amplitude) or troughs (events with a negative-valued amplitude) are more dominant in the signal (also called signal/value skewness). A common strategy for waveshape analysis involves identifying and characterising the features of waves in time-series data - see Cole *et al.* (2017) :footcite:`Cole2017` for an example. Naturally, it can be of interest to explore the shapes of signals at particular frequencies, in which case the time-series data can be bandpass filtered. There is, however, a major limitation to this approach: applying a bandpass filter to data can seriously alter non-sinusoidal signal shape, compromising any analysis of waveshape before it has begun. Thankfully, the bispectrum captures information about non-sinudoisal waveshape, enabling spectrally-resolved analyses at a fine frequency resolution without the need for bandpass filtering. In particular, the bispectrum contains information about rise-decay asymmetry (encoded in the imaginary part of the bispectrum) and peak-trough asymmetry (encoded in the real part of the bispectrum) :footcite:`Bartz2019`. The bispectrum has the general form :math:`\textbf{B}_{kmn}(f_1,f_2)=<\textbf{k}(f_1)\textbf{m}(f_2)\textbf{n}^* (f_2+f_1)>` , where :math:`kmn` is a combination of signals with Fourier coefficients :math:`\textbf{k}`, :math:`\textbf{m}`, and :math:`\textbf{n}`, respectively; :math:`f_1` and :math:`f_2` correspond to a lower and higher frequency, respectively; and :math:`<>` represents the average value over epochs. When analysing waveshape, we are interested in only a single signal, and as such :math:`k=m=n`. Furthermore, we can normalise the bispectrum to the bicoherence, :math:`\boldsymbol{\mathcal{B}}`, using the threenorm, :math:`\textbf{N}`, :footcite:`Shahbazi2014` :math:`\textbf{N}_{xxx}(f_1,f_2)=(<|\textbf{x}(f_1)|^3><|\textbf{x}(f_2)|^3> <|\textbf{x}(f_2+f_1)|^3>)^{\frac{1}{3}}` , :math:`\boldsymbol{\mathcal{B}}_{xxx}(f_1,f_2)=\Large\frac{\textbf{B}_{xxx} (f_1,f_2)}{\textbf{N}_{xxx}(f_1,f_2)}` , where the resulting values lie in the range :math:`[-1, 1]`. .. GENERATED FROM PYTHON SOURCE LINES 77-85 Loading data and computing Fourier coefficients ----------------------------------------------- We will start by loading some example data and computing the Fourier coefficients using the :func:`~pybispectra.utils.compute_fft` function. This data consists of sawtooth waves (information will be captured in the rise-decay asymmetry) and waves with a dominance of peaks or troughs (information will be captured in the peak-trough asymmetry), all simulated as bursting oscillators at 10 Hz. .. GENERATED FROM PYTHON SOURCE LINES 87-141 .. code-block:: Python # load example data data_sawtooths = np.load( get_example_data_paths("sim_data_waveshape_sawtooths") ) data_peaks_troughs = np.load( get_example_data_paths("sim_data_waveshape_peaks_troughs") ) sampling_freq = 1000 # Hz # plot timeseries data times = np.linspace( 0, (data_sawtooths.shape[2] / sampling_freq), data_sawtooths.shape[2] ) fig, axes = plt.subplots(2, 2) axes[0, 0].plot(times, data_sawtooths[15, 0]) axes[0, 1].plot(times, data_sawtooths[15, 1]) axes[1, 0].plot(times, data_peaks_troughs[15, 0]) axes[1, 1].plot(times, data_peaks_troughs[15, 1]) titles = [ "Ramp up sawtooth", "Ramp down sawtooth", "Peak dominance", "Trough dominance", ] for ax, title in zip(axes.flatten(), titles): ax.set_title(title) ax.set_xlabel("Time (s)") ax.set_ylabel("Amplitude (A.U.)") fig.tight_layout() # add noise for numerical stability random = RandomState(44) snr = 0.25 datasets = [data_sawtooths, data_peaks_troughs] for data_idx, data in enumerate(datasets): datasets[data_idx] = snr * data + (1 - snr) * random.rand(*data.shape) data_sawtooths = datasets[0] data_peaks_troughs = datasets[1] # compute Fourier coeffs. fft_coeffs_sawtooths, freqs = compute_fft( data=data_sawtooths, sampling_freq=sampling_freq, n_points=sampling_freq, verbose=False, ) fft_coeffs_peaks_troughs, _ = compute_fft( data=data_peaks_troughs, sampling_freq=sampling_freq, n_points=sampling_freq, verbose=False, ) .. image-sg:: /auto_examples/images/sphx_glr_plot_compute_waveshape_001.png :alt: Ramp up sawtooth, Ramp down sawtooth, Peak dominance, Trough dominance :srcset: /auto_examples/images/sphx_glr_plot_compute_waveshape_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 142-164 Plotting the data, we see that the sawtooth waves consist of a signal where the decay steepness is greater than the rise steepness (ramp up sawtooth) and a signal where the rise steepness is greater than the decay steepness (ramp down sawtooth). Additionally, the peak and trough waves consist of a signal where peaks are most dominant, and a signal where troughs are most dominant. After loading the data, we add some noise for numerical stability. Computing waveshape features ---------------------------- To compute waveshape, we start by initialising the :class:`~pybispectra.waveshape.WaveShape` class object with the Fourier coefficients and the frequency information. To compute waveshape, we call the :meth:`~pybispectra.waveshape.WaveShape.compute` method. By default, waveshape is computed for all channels and all frequency combinations, however we can also specify particular channels and combinations of interest. Here, we specify the frequency arguments :attr:`~pybispectra.waveshape.WaveShape.f1s` and :attr:`~pybispectra.waveshape.WaveShape.f2s` to compute waveshape on in the range 5-35 Hz (around the frequency at which the signal features were simulated). By leaving the indices argument blank, we will look at all channels in the data. .. GENERATED FROM PYTHON SOURCE LINES 166-193 .. code-block:: Python # sawtooth waves waveshape_sawtooths = WaveShape( data=fft_coeffs_sawtooths, freqs=freqs, sampling_freq=sampling_freq, verbose=False, ) # initialise object waveshape_sawtooths.compute(f1s=(5, 35), f2s=(5, 35)) # compute waveshape # peaks and troughs waveshape_peaks_troughs = WaveShape( data=fft_coeffs_peaks_troughs, freqs=freqs, sampling_freq=sampling_freq, verbose=False, ) waveshape_peaks_troughs.compute(f1s=(5, 35), f2s=(5, 35)) # compute waveshape # return results as an array waveshape_results = waveshape_sawtooths.results.get_results() print( f"Waveshape results: [{waveshape_results.shape[0]} channels x " f"{waveshape_results.shape[1]} f1s x {waveshape_results.shape[2]} f2s]" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Waveshape results: [2 channels x 31 f1s x 31 f2s] .. GENERATED FROM PYTHON SOURCE LINES 194-200 We can see that waveshape features have been computed for both channels and the specified frequency combinations, averaged across our epochs. Given the nature of the bispectrum, entries where :math:`f_1` would be higher than :math:`f_2`, as well as where :math:`f_2 + f_1` exceeds the frequency bounds of our data, cannot be computed. In such cases, the values corresponding to those 'bad' frequency combinations are :obj:`numpy.nan`. .. GENERATED FROM PYTHON SOURCE LINES 202-247 Plotting waveshape features --------------------------- Let us now inspect the results. Information about the different waveshape features are encoded in different aspects of the complex-valued bicoherence, with peak-trough asymmetry encoded in the real part, and rise-decay asymmetry encoded in the imaginary part. We can therefore additionally examine the absolute value of the bicoherence (i.e. the magnitude) as well as the phase angle to get an overall picture of the combination of peak-trough and rise-decay asymmetries. For the sawtooth waves, we expect the real part of bicoherence to be ~0 and the imaginary part to be non-zero at the simulated 10 Hz frequency. From the plots, we see that this is indeed the case. However, we also see that the imaginary values at the 10 Hz higher harmonics (i.e. 20 and 30 Hz) are also non-zero, a product of the Fourier transform's application to non-sinusoidal signals. It is also worth noting that the sign of the imaginary values varies for the particular sawtooth type, with a ramp up sawtooth resulting in positive values, and a ramp down sawtooth resulting in negative values. Information about the direction of the asymmetry is encoded not only in the sign of the bicoherence values, but also in its phase. As in Bartz *et al.* :footcite:`Bartz2019`, we represent phase in the range :math:`(0, 2\pi]` (travelling counter-clockwise from the positive real axis). Accordingly, a phase of :math:`\frac{1}{2}\pi` is seen at 10 Hz and its higher harmonics for the ramp up sawtooth, with a phase of :math:`\frac{3}{2}\pi` for the ramp down sawtooth. The phases and absolute values (i.e. the magnitude) therefore combine information from both the real and imaginary components. In contrast, we expect the real part of the bicoherence to be non-zero for signals with peak-trough asymmetry, and the imaginary part to be ~0. Again, this is what we observe. Similarly to before, the signs of the real values are positive for the peak-dominant signal, and negative for the trough-dominant signal, which is also reflected in the phases (~0 or 2 :math:`\pi` for the peak-dominant signal, and :math:`\pi` for the trough-dominant signal). Here, we plotted the real and imaginary parts of the bicoherence without taking the absolute value. If the particular direction of asymmetry is not of interest, the absolute values can be plotted instead (by setting ``plot_absolute=True``) to show the overall degree of asymmetry. In any case, the direction of asymmetry can be inferred from the phases. Finally, note that the :class:`~matplotlib.figure.Figure` and :class:`~matplotlib.axes.Axes` objects can also be returned for any desired manual adjustments of the plots. .. GENERATED FROM PYTHON SOURCE LINES 249-282 .. code-block:: Python figs, axes = waveshape_sawtooths.results.plot( major_tick_intervals=10, minor_tick_intervals=2, cbar_range_abs=(0, 1), cbar_range_real=(-1, 1), cbar_range_imag=(-1, 1), cbar_range_phase=(0, 2), plot_absolute=False, show=False, ) titles = ["Ramp up", "Ramp down"] for fig, title in zip(figs, titles): fig.suptitle(f"{title} sawtooth") fig.set_size_inches(6, 6) fig.show() figs, axes = waveshape_peaks_troughs.results.plot( major_tick_intervals=10, minor_tick_intervals=2, cbar_range_abs=(0, 1), cbar_range_real=(-1, 1), cbar_range_imag=(-1, 1), cbar_range_phase=(0, 2), plot_absolute=False, show=False, ) titles = ["Peak", "Trough"] for fig, title in zip(figs, titles): fig.suptitle(f"{title} dominance") fig.set_size_inches(6, 6) fig.show() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/images/sphx_glr_plot_compute_waveshape_002.png :alt: Ramp up sawtooth, Absolute, Real, Imaginary, Phase :srcset: /auto_examples/images/sphx_glr_plot_compute_waveshape_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_plot_compute_waveshape_003.png :alt: Ramp down sawtooth, Absolute, Real, Imaginary, Phase :srcset: /auto_examples/images/sphx_glr_plot_compute_waveshape_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_plot_compute_waveshape_004.png :alt: Peak dominance, Absolute, Real, Imaginary, Phase :srcset: /auto_examples/images/sphx_glr_plot_compute_waveshape_004.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_plot_compute_waveshape_005.png :alt: Trough dominance, Absolute, Real, Imaginary, Phase :srcset: /auto_examples/images/sphx_glr_plot_compute_waveshape_005.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 283-293 Analysing waveshape in low signal-to-noise ratio data ----------------------------------------------------- Depending on the degree of signal-to-noise ratio as well as the colour of the noise, the ability of the bispectrum to extract information about the underlying waveshape features can vary. To alleviate this, Bartz *et al.* :footcite:`Bartz2019` propose utilising spatio-spectral filtering to enhance the signal-to-noise ratio of the data at a frequency band of interest (which has the added benefit of enabling multivariate signal analysis). Details of how spatio-spectral filtering can be incorporated into waveshape analysis are presented in the following example: :doc:`plot_compute_waveshape_noisy_data`. .. GENERATED FROM PYTHON SOURCE LINES 295-298 References ----------------------------------------------------------------------------- .. footbibliography:: .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 4.657 seconds) .. _sphx_glr_download_auto_examples_plot_compute_waveshape.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_compute_waveshape.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_compute_waveshape.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_