{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Compute waveshape features\n\nThis example demonstrates how waveshape features can be computed with\nPyBispectra.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\nfrom numpy.random import RandomState\nfrom matplotlib import pyplot as plt\n\nfrom pybispectra import compute_fft, get_example_data_paths, WaveShape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Background\nWhen analysing signals, important information may be gleaned from a variety\nof features, including their shape. For example, in neuroscience it has been\nsuggested that non-sinusoidal waves may play important roles in physiology\nand pathology, such as waveform sharpness reflecting synchrony of synaptic\ninputs :footcite:`Sherman2016` and correlating with symptoms of Parkinson's\ndisease :footcite:`Cole2017`. Two aspects of waveshape described in recent\nliterature include: rise-decay asymmetry - how much the wave resembles a\nsawtooth pattern (also called waveform sharpness or derivative skewness); and\npeak-trough asymmetry - whether peaks (events with a positive-valued\namplitude) or troughs (events with a negative-valued amplitude) are more\ndominant in the signal (also called signal/value skewness).\n\nA common strategy for waveshape analysis involves identifying and\ncharacterising the features of waves in time-series data - see Cole *et al.*\n(2017) :footcite:`Cole2017` for an example. Naturally, it can be of interest\nto explore the shapes of signals at particular frequencies, in which case the\ntime-series data can be bandpass filtered. There is, however, a major\nlimitation to this approach: applying a bandpass filter to data can\nseriously alter non-sinusoidal signal shape, compromising any analysis of\nwaveshape before it has begun.\n\nThankfully, the bispectrum captures information about non-sinudoisal\nwaveshape, enabling spectrally-resolved analyses at a fine frequency\nresolution without the need for bandpass filtering. In particular, the\nbispectrum contains information about rise-decay asymmetry (encoded in the\nimaginary part of the bispectrum) and peak-trough asymmetry (encoded in the\nreal part of the bispectrum) :footcite:`Bartz2019`.\n\nThe bispectrum has the general form\n\n$\\textbf{B}_{kmn}(f_1,f_2)=<\\textbf{k}(f_1)\\textbf{m}(f_2)\\textbf{n}^*\n(f_2+f_1)>$ ,\n\nwhere $kmn$ is a combination of signals with Fourier coefficients\n$\\textbf{k}$, $\\textbf{m}$, and $\\textbf{n}$, respectively;\n$f_1$ and $f_2$ correspond to a lower and higher frequency,\nrespectively; and $<>$ represents the average value over epochs. When\nanalysing waveshape, we are interested in only a single signal, and as such\n$k=m=n$.\n\nFurthermore, we can normalise the bispectrum to the bicoherence,\n$\\boldsymbol{\\mathcal{B}}$, using the threenorm, $\\textbf{N}$,\n:footcite:`Shahbazi2014`\n\n$\\textbf{N}_{xxx}(f_1,f_2)=(<|\\textbf{x}(f_1)|^3><|\\textbf{x}(f_2)|^3>\n<|\\textbf{x}(f_2+f_1)|^3>)^{\\frac{1}{3}}$ ,\n\n$\\boldsymbol{\\mathcal{B}}_{xxx}(f_1,f_2)=\\Large\\frac{\\textbf{B}_{xxx}\n(f_1,f_2)}{\\textbf{N}_{xxx}(f_1,f_2)}$ ,\n\nwhere the resulting values lie in the range $[-1, 1]$.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Loading data and computing Fourier coefficients\nWe will start by loading some example data and computing the Fourier\ncoefficients using the :func:`~pybispectra.utils.compute_fft` function. This\ndata consists of sawtooth waves (information will be captured in the\nrise-decay asymmetry) and waves with a dominance of peaks or troughs\n(information will be captured in the peak-trough asymmetry), all simulated as\nbursting oscillators at 10 Hz.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# load example data\ndata_sawtooths = np.load(\n get_example_data_paths(\"sim_data_waveshape_sawtooths\")\n)\ndata_peaks_troughs = np.load(\n get_example_data_paths(\"sim_data_waveshape_peaks_troughs\")\n)\nsampling_freq = 1000 # Hz\n\n# plot timeseries data\ntimes = np.linspace(\n 0, (data_sawtooths.shape[2] / sampling_freq), data_sawtooths.shape[2]\n)\nfig, axes = plt.subplots(2, 2)\naxes[0, 0].plot(times, data_sawtooths[15, 0])\naxes[0, 1].plot(times, data_sawtooths[15, 1])\naxes[1, 0].plot(times, data_peaks_troughs[15, 0])\naxes[1, 1].plot(times, data_peaks_troughs[15, 1])\ntitles = [\n \"Ramp up sawtooth\",\n \"Ramp down sawtooth\",\n \"Peak dominance\",\n \"Trough dominance\",\n]\nfor ax, title in zip(axes.flatten(), titles):\n ax.set_title(title)\n ax.set_xlabel(\"Time (s)\")\n ax.set_ylabel(\"Amplitude (A.U.)\")\nfig.tight_layout()\n\n# add noise for numerical stability\nrandom = RandomState(44)\nsnr = 0.25\ndatasets = [data_sawtooths, data_peaks_troughs]\nfor data_idx, data in enumerate(datasets):\n datasets[data_idx] = snr * data + (1 - snr) * random.rand(*data.shape)\ndata_sawtooths = datasets[0]\ndata_peaks_troughs = datasets[1]\n\n# compute Fourier coeffs.\nfft_coeffs_sawtooths, freqs = compute_fft(\n data=data_sawtooths,\n sampling_freq=sampling_freq,\n n_points=sampling_freq,\n verbose=False,\n)\nfft_coeffs_peaks_troughs, _ = compute_fft(\n data=data_peaks_troughs,\n sampling_freq=sampling_freq,\n n_points=sampling_freq,\n verbose=False,\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting the data, we see that the sawtooth waves consist of a signal where\nthe decay steepness is greater than the rise steepness (ramp up sawtooth) and\na signal where the rise steepness is greater than the decay steepness (ramp\ndown sawtooth). Additionally, the peak and trough waves consist of a signal\nwhere peaks are most dominant, and a signal where troughs are most dominant.\nAfter loading the data, we add some noise for numerical stability.\n\n## Computing waveshape features\nTo compute waveshape, we start by initialising the\n:class:`~pybispectra.waveshape.WaveShape` class object with the Fourier\ncoefficients and the frequency information. To compute waveshape, we call the\n:meth:`~pybispectra.waveshape.WaveShape.compute` method. By default,\nwaveshape is computed for all channels and all frequency combinations,\nhowever we can also specify particular channels and combinations of interest.\n\nHere, we specify the frequency arguments\n:attr:`~pybispectra.waveshape.WaveShape.f1s` and\n:attr:`~pybispectra.waveshape.WaveShape.f2s` to compute waveshape on in the\nrange 5-35 Hz (around the frequency at which the signal features were\nsimulated). By leaving the indices argument blank, we will look at all\nchannels in the data.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# sawtooth waves\nwaveshape_sawtooths = WaveShape(\n data=fft_coeffs_sawtooths,\n freqs=freqs,\n sampling_freq=sampling_freq,\n verbose=False,\n) # initialise object\nwaveshape_sawtooths.compute(f1s=(5, 35), f2s=(5, 35)) # compute waveshape\n\n# peaks and troughs\nwaveshape_peaks_troughs = WaveShape(\n data=fft_coeffs_peaks_troughs,\n freqs=freqs,\n sampling_freq=sampling_freq,\n verbose=False,\n)\nwaveshape_peaks_troughs.compute(f1s=(5, 35), f2s=(5, 35)) # compute waveshape\n\n# return results as an array\nwaveshape_results = waveshape_sawtooths.results.get_results()\n\nprint(\n f\"Waveshape results: [{waveshape_results.shape[0]} channels x \"\n f\"{waveshape_results.shape[1]} f1s x {waveshape_results.shape[2]} f2s]\"\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that waveshape features have been computed for both channels and\nthe specified frequency combinations, averaged across our epochs. Given the\nnature of the bispectrum, entries where $f_1$ would be higher than\n$f_2$, as well as where $f_2 + f_1$ exceeds the frequency bounds\nof our data, cannot be computed. In such cases, the values corresponding to\nthose 'bad' frequency combinations are :obj:`numpy.nan`.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plotting waveshape features\nLet us now inspect the results. Information about the different waveshape\nfeatures are encoded in different aspects of the complex-valued bicoherence,\nwith peak-trough asymmetry encoded in the real part, and rise-decay asymmetry\nencoded in the imaginary part. We can therefore additionally examine the\nabsolute value of the bicoherence (i.e. the magnitude) as well as the phase\nangle to get an overall picture of the combination of peak-trough and\nrise-decay asymmetries.\n\nFor the sawtooth waves, we expect the real part of bicoherence to be ~0 and\nthe imaginary part to be non-zero at the simulated 10 Hz frequency. From the\nplots, we see that this is indeed the case. However, we also see that the\nimaginary values at the 10 Hz higher harmonics (i.e. 20 and 30 Hz) are also\nnon-zero, a product of the Fourier transform's application to non-sinusoidal\nsignals. It is also worth noting that the sign of the imaginary values varies\nfor the particular sawtooth type, with a ramp up sawtooth resulting in\npositive values, and a ramp down sawtooth resulting in negative values.\n\nInformation about the direction of the asymmetry is encoded not only in the\nsign of the bicoherence values, but also in its phase. As in Bartz *et al.*\n:footcite:`Bartz2019`, we represent phase in the range $(0, 2\\pi]$\n(travelling counter-clockwise from the positive real axis). Accordingly, a\nphase of $\\frac{1}{2}\\pi$ is seen at 10 Hz and its higher harmonics for\nthe ramp up sawtooth, with a phase of $\\frac{3}{2}\\pi$ for the ramp\ndown sawtooth. The phases and absolute values (i.e. the magnitude) therefore\ncombine information from both the real and imaginary components.\n\nIn contrast, we expect the real part of the bicoherence to be non-zero for\nsignals with peak-trough asymmetry, and the imaginary part to be ~0. Again,\nthis is what we observe. Similarly to before, the signs of the real values\nare positive for the peak-dominant signal, and negative for the\ntrough-dominant signal, which is also reflected in the phases (~0 or\n2 $\\pi$ for the peak-dominant signal, and $\\pi$ for the\ntrough-dominant signal).\n\nHere, we plotted the real and imaginary parts of the bicoherence without\ntaking the absolute value. If the particular direction of asymmetry is not of\ninterest, the absolute values can be plotted instead (by setting\n``plot_absolute=True``) to show the overall degree of asymmetry. In any case,\nthe direction of asymmetry can be inferred from the phases.\n\nFinally, note that the :class:`~matplotlib.figure.Figure` and\n:class:`~matplotlib.axes.Axes` objects can also be returned for any desired\nmanual adjustments of the plots.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"figs, axes = waveshape_sawtooths.results.plot(\n major_tick_intervals=10,\n minor_tick_intervals=2,\n cbar_range_abs=(0, 1),\n cbar_range_real=(-1, 1),\n cbar_range_imag=(-1, 1),\n cbar_range_phase=(0, 2),\n plot_absolute=False,\n show=False,\n)\ntitles = [\"Ramp up\", \"Ramp down\"]\nfor fig, title in zip(figs, titles):\n fig.suptitle(f\"{title} sawtooth\")\n fig.set_size_inches(6, 6)\n fig.show()\n\nfigs, axes = waveshape_peaks_troughs.results.plot(\n major_tick_intervals=10,\n minor_tick_intervals=2,\n cbar_range_abs=(0, 1),\n cbar_range_real=(-1, 1),\n cbar_range_imag=(-1, 1),\n cbar_range_phase=(0, 2),\n plot_absolute=False,\n show=False,\n)\ntitles = [\"Peak\", \"Trough\"]\nfor fig, title in zip(figs, titles):\n fig.suptitle(f\"{title} dominance\")\n fig.set_size_inches(6, 6)\n fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Analysing waveshape in low signal-to-noise ratio data\nDepending on the degree of signal-to-noise ratio as well as the colour of the\nnoise, the ability of the bispectrum to extract information about the\nunderlying waveshape features can vary. To alleviate this, Bartz *et al.*\n:footcite:`Bartz2019` propose utilising spatio-spectral filtering to enhance\nthe signal-to-noise ratio of the data at a frequency band of interest (which\nhas the added benefit of enabling multivariate signal analysis). Details of\nhow spatio-spectral filtering can be incorporated into waveshape analysis are\npresented in the following example: :doc:`plot_compute_waveshape_noisy_data`.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n.. footbibliography::\n\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.6"
}
},
"nbformat": 4,
"nbformat_minor": 0
}