Skip to content

Commit

Permalink
Update formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
@tsbinns
tsbinns committed Sep 18, 2024
1 parent 80a0af7 commit 159502b
Show file tree
Hide file tree
Showing 7 changed files with 217 additions and 273 deletions.
123 changes: 52 additions & 71 deletions examples/plot_example_dbs_data.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,15 +3,15 @@
Real data example: stimulation artefacts in human ECoG data
===========================================================
This example demonstrates the utility of PARRM :footcite:`DastinEtAl2021` on a
genuine recording of human brain activity from local field potentials (LFPs) at
the site of deep brain stimulation in the subthalamic nucleus, and from
electrocorticography (ECoG) at the cortex.
This example demonstrates the utility of PARRM :footcite:`DastinEtAl2021` on a genuine
recording of human brain activity from local field potentials (LFPs) at the site of deep
brain stimulation in the subthalamic nucleus, and from electrocorticography (ECoG) at
the cortex.
"""

# Author(s):
# Timon Merk | github.com/timonmerk
# Thomas Samuel Binns | github.com/tsbinns
# Timon Merk | github.com/timonmerk
# Thomas S. Binns | github.com/tsbinns

# %%

Expand All @@ -21,24 +21,23 @@
from pyparrm import get_example_data_paths, PARRM
from pyparrm._utils._power import compute_psd

###############################################################################
########################################################################################
# Background
# ----------
# Deep brain stimulation (DBS) is an established treatment for several
# disorders such as Parkinson's disease, epilepsy, and dystonia. In subjects
# undergoing DBS, brain activity can also be recorded. Unfortunately, the
# quality of signals recorded at not only the site of stimulation, but often
# from distal regions too, can be detrimentally affected by stimulation-related
# artefacts.
# Deep brain stimulation (DBS) is an established treatment for several disorders such as
# Parkinson's disease, epilepsy, and dystonia. In subjects undergoing DBS, brain
# activity can also be recorded. Unfortunately, the quality of signals recorded at not
# only the site of stimulation, but often from distal regions too, can be detrimentally
# affected by stimulation-related artefacts.
#
# In this example from a Parkinson's disease patient undergoing DBS of the
# subthalamic nucleus - a common target for stimulation - we show how PARRM can
# be used to remove stimulation artefacts from both LFP recordings of the
# subthalamic nucleus and ECoG recordings of cortical activity.
# In this example from a Parkinson's disease patient undergoing DBS of the subthalamic
# nucleus - a common target for stimulation - we show how PARRM can be used to remove
# stimulation artefacts from both LFP recordings of the subthalamic nucleus and ECoG
# recordings of cortical activity.
#
# We start by loading some ECoG and LFP data collected during rest, each
# consisting of an individual channel spanning 60 seconds, with a sampling
# frequency of 1,000 Hz. DBS was delivered at a frequency of 130 Hz.
# We start by loading some ECoG and LFP data collected during rest, each consisting of
# an individual channel spanning 60 seconds, with a sampling frequency of 1,000 Hz. DBS
# was delivered at a frequency of 130 Hz.

# %%

Expand All @@ -47,27 +46,25 @@
sampling_freq = 1000 # Hz
artefact_freq = 130 # Hz

###############################################################################
########################################################################################
# Finding the artefact period and filtering the data
# --------------------------------------------------
# When handling data from multiple sites, the decision must be made whether to
# estimate the period of the stimulation artefacts for this data separately, or
# together. Assuming that the stimulation artefacts have a similar effect on
# the subthalamic and cortical signals, we can take advantage of the more
# efficient approach of estimating the period for all data simultaneously.
# However, if differences in the period between recording sites are of a
# particular concern, we can also estimate the periods separately. Below, we
# estimate the periods on: the ECoG and LFP data together; the ECoG data alone;
# and the LFP data alone. In this case, the estimate of the periods are
# When handling data from multiple sites, the decision must be made whether to estimate
# the period of the stimulation artefacts for this data separately, or together.
# Assuming that the stimulation artefacts have a similar effect on the subthalamic and
# cortical signals, we can take advantage of the more efficient approach of estimating
# the period for all data simultaneously. However, if differences in the period between
# recording sites are of a particular concern, we can also estimate the periods
# separately. Below, we estimate the periods on: the ECoG and LFP data together; the
# ECoG data alone; and the LFP data alone. In this case, the estimate of the periods are
# identical to 5 decimal places. Nevertheless, if you can afford the increased
# computational cost, you may see improved performance for estimating the
# period independently for recordings from different sites.
# computational cost, you may see improved performance for estimating the period
# independently for recordings from different sites.
#
# Having estimated the period, we proceed to create a filter and apply it to
# the data. The following examples give a detailed explanation of the process
# for finding the artefact period and filtering the data, as well as for
# finding the optimal filter parameters for your data: :doc:`plot_use_parrm`;
# :doc:`plot_use_param_explorer`.
# Having estimated the period, we proceed to create a filter and apply it to the data.
# The following examples give a detailed explanation of the process for finding the
# artefact period and filtering the data, as well as for finding the optimal filter
# parameters for your data: :doc:`plot_use_parrm`; :doc:`plot_use_param_explorer`.

# %%

Expand Down Expand Up @@ -106,26 +103,24 @@
f"Artefact period (only LFP): {parrm_lfp.period :.7f}"
)

###############################################################################
########################################################################################
# Inspecting the results
# ----------------------
# Having filtered the data, we can now compare the results to the original,
# artefact-contaminated signals. The plots below show the entire 60 second
# timeseries of the data, as well as the power spectral densities across this
# period.
# artefact-contaminated signals. The plots below show the entire 60 second timeseries of
# the data, as well as the power spectral densities across this period.
#
# For both the ECoG and LFP data, the reduction in power at the 130 Hz
# stimulation frequency and its harmonics (260 Hz and 390 Hz) is readily
# apparent. Whilst some spikes in power at these frequencies remain, the
# activity is much less broadband, and as such the information in the
# neighbouring frequencies is less compromised. Finally, it can be seen that
# PARRM performs well for both the ECoG data from the cortex - located distally
# from the site of stimulation - as well as for the LFP data - collected at the
# site of DBS.
# For both the ECoG and LFP data, the reduction in power at the 130 Hz stimulation
# frequency and its harmonics (260 Hz and 390 Hz) is readily apparent. Whilst some
# spikes in power at these frequencies remain, the activity is much less broadband, and
# as such the information in the neighbouring frequencies is less compromised. Finally,
# it can be seen that PARRM performs well for both the ECoG data from the cortex -
# located distally from the site of stimulation - as well as for the LFP data -
# collected at the site of DBS.
#
# Altogether, it is clear that PARRM is a powerful tool for removing periodic
# stimulation artefacts from electrophysiological data in a
# computationally-efficient manner.
# stimulation artefacts from electrophysiological data in a computationally-efficient
# manner.

# %%

Expand All @@ -138,15 +133,9 @@
inset_end_idx = int(30.5 * sampling_freq)

n_freqs = sampling_freq / 2
psd_freqs, psd_unfiltered = compute_psd(
data_ecog_lfp, sampling_freq, int(n_freqs * 2)
)
_, psd_filtered_ecog = compute_psd(
filtered_ecog, sampling_freq, int(n_freqs * 2)
)
_, psd_filtered_lfp = compute_psd(
filtered_lfp, sampling_freq, int(n_freqs * 2)
)
psd_freqs, psd_unfiltered = compute_psd(data_ecog_lfp, sampling_freq, int(n_freqs * 2))
_, psd_filtered_ecog = compute_psd(filtered_ecog, sampling_freq, int(n_freqs * 2))
_, psd_filtered_lfp = compute_psd(filtered_lfp, sampling_freq, int(n_freqs * 2))

fig = plt.figure(figsize=(12, 8), layout="constrained")
subfigs = fig.subfigures(2, 1, hspace=0.07)
Expand All @@ -169,10 +158,7 @@
label="Unfiltered",
)
axes[0].plot(
times,
filtered_data[start_idx : end_idx + 1],
color="#1f77b4",
label="Filtered",
times, filtered_data[start_idx : end_idx + 1], color="#1f77b4", label="Filtered"
)
inset_axis_time.plot(
inset_times,
Expand All @@ -181,18 +167,13 @@
alpha=0.7,
)
inset_axis_time.plot(
inset_times,
filtered_data[inset_start_idx : inset_end_idx + 1],
color="#1f77b4",
inset_times, filtered_data[inset_start_idx : inset_end_idx + 1], color="#1f77b4"
)
inset_axis_time.set_xlim(30, 30.5)
inset_data = filtered_data[inset_start_idx : inset_end_idx + 1]
inset_ylim_pad = (np.max(inset_data) - np.min(inset_data)) * 0.1
inset_ylim = np.array(
(
np.min(inset_data) - inset_ylim_pad,
np.max(inset_data) + inset_ylim_pad,
)
(np.min(inset_data) - inset_ylim_pad, np.max(inset_data) + inset_ylim_pad)
)
inset_axis_time.set_ylim(inset_ylim)
axes[0].indicate_inset_zoom(inset_axis_time, edgecolor="black", alpha=0.4)
Expand Down
119 changes: 56 additions & 63 deletions examples/plot_use_param_explorer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,95 +3,88 @@
Exploring the best filter parameters for your data
==================================================
This example demonstrates how the PyPARRM interactive parameter explorer can be
used to find the best filter settings for your data.
This example demonstrates how the PyPARRM interactive parameter explorer can be used to
find the best filter settings for your data.
"""

# Author(s):
# Thomas Samuel Binns | github.com/tsbinns
# Thomas S. Binns | github.com/tsbinns

###############################################################################
########################################################################################
# Background
# ----------
# Once you have an estimate of the artefact period, you need to design a filter
# to remove the artefacts from your data. There are four key filter parameters:
# Once you have an estimate of the artefact period, you need to design a filter to
# remove the artefacts from your data. There are four key filter parameters:
#
# 1. The size of the filter window. This should be chosen based on the
# timescale of which the artefact shape varies.
# 1. The size of the filter window. This should be chosen based on the timescale of
# which the artefact shape varies.
#
# 2. The number of samples to omit from the centre of the filter window. This
# parameter serves to control for overfitting to features of interest in the
# underlying data.
# 2. The number of samples to omit from the centre of the filter window. This parameter
# serves to control for overfitting to features of interest in the underlying data.
#
# 3. The direction considered when building the filter. This can be used to
# control whether the filter window takes only previous samples, future
# samples, or both previous and future samples into account, based on their
# position relative to the centre of the filter window.
# 3. The direction considered when building the filter. This can be used to control
# whether the filter window takes only previous samples, future samples, or both
# previous and future samples into account, based on their position relative to the
# centre of the filter window.
#
# 4. The period window size. The size of this window should be based on changes
# in the waveform of the artefact. This parameter controls which samples are
# combined on the timescale of the period.
# 4. The period window size. The size of this window should be based on changes in the
# waveform of the artefact. This parameter controls which samples are combined on the
# timescale of the period.
#
# These parameters are described in detail in Dastin *et al.* (2021)
# :footcite:`DastinEtAl2021` and the supplementary `video
# <https://ars.els-cdn.com/content/image/1-s2.0-S2667237521000102-mmc2.mp4>`_
# from the paper's authors. Crucially, visualising the effects of these
# parameters can help greatly for identifying the best settings to use when
# creating a filter for your data.
# <https://ars.els-cdn.com/content/image/1-s2.0-S2667237521000102-mmc2.mp4>`_ from the
# paper's authors. Crucially, visualising the effects of these parameters can help
# greatly for identifying the best settings to use when creating a filter for your data.

###############################################################################
########################################################################################
# Exploring the filter parameters
# -------------------------------
# Once you have an an estimate of the artefact period (see the following
# example for a detailed explanation of how to do this: :doc:`plot_use_parrm`),
# you can now visualise the effects of the filter parameters to find those that
# best remove the stimulation artefacts from the data. To do this, we call the
# :meth:`explore_filter_params <pyparrm.parrm.PARRM.explore_filter_params>`
# method.
# Once you have an an estimate of the artefact period (see the following example for a
# detailed explanation of how to do this: :doc:`plot_use_parrm`), you can now visualise
# the effects of the filter parameters to find those that best remove the stimulation
# artefacts from the data. To do this, we call the
# :meth:`~pyparrm.parrm.PARRM.explore_filter_params` method.
#
# The explorer consists of four plots, and a set of controls:
#
# * The textboxes and buttons in the bottom left of the explorer allow you to
# manipulate the filter settings. The default settings of the explorer
# reflect the settings that would be used when computing the filter if you
# did not specify any inputs. You can change: the period half-width
# (``period_half_width``); the filter half-width (``filter_half_width``);
# the number of omitted samples (``omit_n_samples``); and the filter
# direction (``filter_direction``). The ranges of possible values that can be
# * The textboxes and buttons in the bottom left of the explorer allow you to manipulate
# the filter settings. The default settings of the explorer reflect the settings that
# would be used when computing the filter if you did not specify any inputs. You can
# change: the period half-width (``period_half_width``); the filter half-width
# (``filter_half_width``); the number of omitted samples (``omit_n_samples``); and the
# filter direction (``filter_direction``). The ranges of possible values that can be
# entered in the textboxes are shown in brackets.
#
# * The top left plot shows the data samples in the period space according to
# the current period half-width. Just below this is an overview plot of all
# the samples in the period space. The red area shows the region currently
# being viewed in the top left plot. The period space can be moved through
# using the left and right arrow keys. With these plots, you can find a
# suitable period half-width, such that the waveform of the artefact remains
# fairly constant (i.e. a straight line could easily be drawn through the
# plotted samples in period space). Furthermore, as you change the filter
# half-width, these plots will reflect the change in the number of available
# samples.
# * The top left plot shows the data samples in the period space according to the
# current period half-width. Just below this is an overview plot of all the samples in
# the period space. The red area shows the region currently being viewed in the top
# left plot. The period space can be moved through using the left and right arrow
# keys. With these plots, you can find a suitable period half-width, such that the
# waveform of the artefact remains fairly constant (i.e. a straight line could easily
# be drawn through the plotted samples in period space). Furthermore, as you change
# the filter half-width, these plots will reflect the change in the number of
# available samples.
#
# * Most importantly, the two plots on the right-hand side show the unfiltered
# and filtered data according to the current filter settings in the time
# domain (top plot) and in the frequency domain (bottom plot). Both plots are
# updated in real-time according to changes in the filter settings. The
# power spectral density is particularly useful for identifying which filters
# reduce the signal's power at the artefact frequency (as well as at the
# harmonics and sub-harmonics). In some cases, the filter settings may not
# produce a valid filter, in which case a warning on the time domain plot
# will be printed.
# * Most importantly, the two plots on the right-hand side show the unfiltered and
# filtered data according to the current filter settings in the time domain (top plot)
# and in the frequency domain (bottom plot). Both plots are updated in real-time
# according to changes in the filter settings. The power spectral density is
# particularly useful for identifying which filters reduce the signal's power at the
# artefact frequency (as well as at the harmonics and sub-harmonics). In some cases,
# the filter settings may not produce a valid filter, in which case a warning on the
# time domain plot will be printed.
#
# If computational cost is a concern, a limited time range of the data and its
# resolution can be specified, as well as the frequency range and resolution of
# the power spectral density. If multiple channels are present in your data,
# these can be navigated between using the up and down arrow keys. Finally, the
# title of the explorer gives you an overview of the current filter settings,
# as well as which channel's data is currently being viewed.
# resolution can be specified, as well as the frequency range and resolution of the
# power spectral density. If multiple channels are present in your data, these can be
# navigated between using the up and down arrow keys. Finally, the title of the explorer
# gives you an overview of the current filter settings, as well as which channel's data
# is currently being viewed.
#
# The method can be called as ``parrm.explore_filter_params()``. The image
# below shows the different aspects of the parameter explorer window.
# The image below shows the different aspects of the parameter explorer window.

###############################################################################
########################################################################################
# .. figure:: ../_static/param_explorer.png
# :alt: parameter explorer window overview
#
Expand Down
Loading

0 comments on commit 159502b

Please sign in to comment.