diff --git a/examples/plot_example_dbs_data.py b/examples/plot_example_dbs_data.py index 60d82c7..8f91aae 100644 --- a/examples/plot_example_dbs_data.py +++ b/examples/plot_example_dbs_data.py @@ -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 # %% @@ -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. # %% @@ -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`. # %% @@ -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. # %% @@ -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) @@ -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, @@ -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) diff --git a/examples/plot_use_param_explorer.py b/examples/plot_use_param_explorer.py index 3205247..314e33e 100644 --- a/examples/plot_use_param_explorer.py +++ b/examples/plot_use_param_explorer.py @@ -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 -# `_ -# 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. +# `_ 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 ` -# 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 # diff --git a/examples/plot_use_parrm.py b/examples/plot_use_parrm.py index 470181f..c2eb449 100644 --- a/examples/plot_use_parrm.py +++ b/examples/plot_use_parrm.py @@ -3,13 +3,13 @@ Using PyPARRM to filter out stimulation artefacts from data =========================================================== -This example demonstrates how the PARRM algorithm :footcite:`DastinEtAl2021` -can be used to identify and remove stimulation artefacts from -electrophysiological data in the PyPARRM package. +This example demonstrates how the PARRM algorithm :footcite:`DastinEtAl2021` can be used +to identify and remove stimulation artefacts from electrophysiological data in the +PyPARRM package. """ # Author(s): -# Thomas Samuel Binns | github.com/tsbinns +# Thomas S. Binns | github.com/tsbinns # %% @@ -19,26 +19,24 @@ from pyparrm import get_example_data_paths, PARRM from pyparrm._utils._power import compute_psd -############################################################################### +######################################################################################## # Background # ---------- -# When delivering electrical stimulation to bioligical tissues for research or -# clinical purposes, it is often the case that electrophysiological recordings -# collected during this time are contaminated by stimulation artefacts. This -# contamination makes it difficult to analyse the underlying physiological or -# pathophysiological electrical activity. To this end, the Period-based -# Artefact Reconstruction and Removal Method (PARRM) was developed, enabling -# the removal of stimulation arefacts from electrophysiological recordings in a -# robust and computationally cheap manner :footcite:`DastinEtAl2021`. **N.B.** -# PARRM assumes that the artefacts are semi-regular, periodic, and linearly -# combined with the signal of interest. +# When delivering electrical stimulation to bioligical tissues for research or clinical +# purposes, it is often the case that electrophysiological recordings collected during +# this time are contaminated by stimulation artefacts. This contamination makes it +# difficult to analyse the underlying physiological or pathophysiological electrical +# activity. To this end, the Period-based Artefact Reconstruction and Removal Method +# (PARRM) was developed, enabling the removal of stimulation arefacts from +# electrophysiological recordings in a robust and computationally cheap manner +# :footcite:`DastinEtAl2021`. **N.B.** PARRM assumes that the artefacts are +# semi-regular, periodic, and linearly combined with the signal of interest. # -# To demonstrate how PARRM can be used to remove stimulation arfectacts from -# data, we will start by loading some example data. This is the same example -# data used in the MATLAB implementation of the method -# (https://github.com/neuromotion/PARRM), consisting of a single channel -# with ~100 seconds of data at a sampling frequency of 200 Hz, and containing -# stimulation artefacts with a frequency of 150 Hz. +# To demonstrate how PARRM can be used to remove stimulation arfectacts from data, we +# will start by loading some example data. This is the same example data used in the +# MATLAB implementation of the method (https://github.com/neuromotion/PARRM), consisting +# of a single channel with ~100 seconds of data at a sampling frequency of 200 Hz, and +# containing stimulation artefacts with a frequency of 150 Hz. # %% @@ -52,94 +50,85 @@ f"`data` duration: {data.shape[1] / sampling_freq :.2f} seconds" ) -############################################################################### +######################################################################################## # Finding the period of the stimulation artefacts # ----------------------------------------------- # Having loaded the example data, we can now find the period of the stimulation -# artefacts, which we require to remove them from the data. Having imported the -# PARRM object, we initialise it, providing the data, its sampling frequency, -# and the stimulation frequency. +# artefacts, which we require to remove them from the data. Having imported the PARRM +# object, we initialise it, providing the data, its sampling frequency, and the +# stimulation frequency. # -# After this, we find the period of the artefact using the :meth:`find_period -# ` method. By default, all timepoints of the -# data will be used for this, and the initial guess of the period will be taken -# as the sampling frequency divided by the artefact frequency. The settings for -# finding the period can be specified in the method call, however the default -# settings should suffice as a starting point for period estimation. +# After this, we find the period of the artefact using the +# :meth:`~pyparrm.parrm.PARRM.find_period` method. By default, all timepoints of the +# data will be used for this, and the initial guess of the period will be taken as the +# sampling frequency divided by the artefact frequency. The settings for finding the +# period can be specified in the method call, however the default settings should +# suffice as a starting point for period estimation. # -# The period is found using a grid search, with the goal of minimising the mean -# squared error between the data and the best fitting sinusoidal harmonics of -# the period found with linear regression. The process is described in detail -# in Dastin *et al.* (2021) :footcite:`DastinEtAl2021`, and is also -# demonstrated in this `video -# `_ -# from the paper's authors. +# The period is found using a grid search, with the goal of minimising the mean squared +# error between the data and the best fitting sinusoidal harmonics of the period found +# with linear regression. The process is described in detail in Dastin *et al.* (2021) +# :footcite:`DastinEtAl2021`, and is also demonstrated in this `video +# `_ from the +# paper's authors. # %% parrm = PARRM( - data=data, - sampling_freq=sampling_freq, - artefact_freq=artefact_freq, - verbose=False, # silenced to reduce pqdm output clutter + data=data, sampling_freq=sampling_freq, artefact_freq=artefact_freq, verbose=False ) parrm.find_period() print(f"Estimated artefact period: {parrm.period :.4f}") -############################################################################### +######################################################################################## # Creating the filter and removing the artefacts # ---------------------------------------------- -# Now that we have an estimate of the artefact period, we can design a filter -# to remove it from the data using the :meth:`create_filter -# ` method. When creating the filter, there -# are four key parameters: +# Now that we have an estimate of the artefact period, we can design a filter to remove +# it from the data using the :meth:`~pyparrm.parrm.PARRM.create_filter` method. When +# creating the filter, there are four key parameters: # -# 1. The size of the filter window, specified with the ``filter_half_width`` -# parameter. This should be chosen based on the timescale on which the -# artefact shape varies. If no such timescale is known, the power spectra of -# the filtered data can be inspected and the size of the filter window tuned -# until artefact-related peaks are sufficiently attenuated. +# 1. The size of the filter window, specified with the ``filter_half_width`` parameter. +# This should be chosen based on the timescale on which the artefact shape varies. If +# no such timescale is known, the power spectra of the filtered data can be inspected +# and the size of the filter window tuned until artefact-related peaks are +# sufficiently attenuated. # -# 2. The number of samples to omit from the centre of the filter window, -# specified with the ``omit_n_samples`` parameter. This parameter serves to -# control for overfitting to features of interest in the underlying data. -# For instance, if there is a physiological signal of interest known to -# occur on a particular timescale, an appropriate number of samples should -# be omitted according to this range of time. +# 2. The number of samples to omit from the centre of the filter window, specified with +# the ``omit_n_samples`` parameter. This parameter serves to control for overfitting +# to features of interest in the underlying data. For instance, if there is a +# physiological signal of interest known to occur on a particular timescale, an +# appropriate number of samples should be omitted according to this range of time. # # 3. The direction considered when building the filter, specified with the -# ``filter_direction`` parameter. 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. +# ``filter_direction`` parameter. 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, specified with the ``period_half_width`` -# parameter. The size of this window should be based on changes in the -# waveform of the artefact, which can be estimated by plotting the data on -# the timescale of the period and identifying the timescale over which -# features remain fairly constant. This parameter controls which samples are -# combined on the timescale of the period. +# 4. The period window size, specified with the ``period_half_width`` parameter. The +# size of this window should be based on changes in the waveform of the artefact, +# which can be estimated by plotting the data on the timescale of the period and +# identifying the timescale over which features remain fairly constant. This +# parameter controls which samples are combined on the timescale of the period. # # If you are unsure of what parameters are best for your data, you can use the # interactive parameter exploration tool offered in PyPARRM (see the -# :meth:`explore_filter_params ` -# method and the associated example: :doc:`plot_use_param_explorer`). +# :meth:`~pyparrm.parrm.PARRM.explore_filter_params` method and the associated example: +# :doc:`plot_use_param_explorer`). # -# Here, we specify that the filter should have a half-width of 2,000 samples, -# ignoring those 20 samples adjacent to the centre of the filter window, -# considering samples both before and after the centre of the filter window, -# and finally using a half-width of 0.01 samples in the period space. +# Here, we specify that the filter should have a half-width of 2,000 samples, ignoring +# those 20 samples adjacent to the centre of the filter window, considering samples both +# before and after the centre of the filter window, and finally using a half-width of +# 0.01 samples in the period space. # # Once the filter has been created, it can be applied using the -# :meth:`filter_data ` method, which returns -# the artefact-free data. By default, -# :meth:`filter_data ` will filter the data -# stored in the :class:`PARRM ` object, however other data -# can be provided should you wish to filter this instead. The filter itself, as -# well as a copy of the filtered data, can be accessed from the :attr:`filter -# ` and :attr:`filtered_data -# ` attributes, respectively. +# :meth:`~pyparrm.parrm.PARRM.filter_data` method, which returns the artefact-free data. +# By default, :meth:`~pyparrm.parrm.PARRM.filter_data` will filter the data stored in +# the :class:`~pyparrm.parrm.PARRM` object, however other data can be provided should +# you wish to filter this instead. The filter itself, as well as a copy of the filtered +# data, can be accessed from the :attr:`~pyparrm.parrm.PARRM.filter` and +# :attr:`~pyparrm.parrm.PARRM.filtered_data` attributes, respectively. # %% @@ -151,18 +140,18 @@ ) filtered_data = parrm.filter_data() # other data to filter can be given here -############################################################################### +######################################################################################## # Inspecting the results # ---------------------- -# Having filtered the data, we can now compare the results to the original -# data, as well as the artefact-free form of this simulated data to see how -# well the method is able to remove the underlying artefacts. +# Having filtered the data, we can now compare the results to the original data, as well +# as the artefact-free form of this simulated data to see how well the method is able to +# remove the underlying artefacts. # -# As you can see, the filtered timeseries data closely resembles the true -# artefact-free data. Furthermore, inspecting the power spectra shows just how -# well PARRM is able to attentuate the overwhelming power at the 50 Hz -# subharmonic of the stimulation artefacts, again bringing the results closely -# in line with those of the true artefact-free data. +# As you can see, the filtered timeseries data closely resembles the true artefact-free +# data. Furthermore, inspecting the power spectra shows just how well PARRM is able to +# attentuate the overwhelming power at the 50 Hz subharmonic of the stimulation +# artefacts, again bringing the results closely in line with those of the true +# artefact-free data. # %% @@ -176,12 +165,8 @@ inset_axis = axes[0].inset_axes((0.12, 0.6, 0.5, 0.35)) # main timeseries plot -axes[0].plot( - times, data[0, start:end], color="black", alpha=0.3, label="Unfiltered" -) -axes[0].plot( - times, artefact_free[0, start:end], linewidth=3, label="Artefact-free" -) +axes[0].plot(times, data[0, start:end], color="black", alpha=0.3, label="Unfiltered") +axes[0].plot(times, artefact_free[0, start:end], linewidth=3, label="Artefact-free") axes[0].plot(times, filtered_data[0, start:end], label="Filtered (PyPARRM)") axes[0].legend() axes[0].set_xlabel("Time (s)") @@ -194,23 +179,13 @@ inset_axis.patch.set_alpha(0.7) # power spectral density plot -n_freqs = sampling_freq / 2 -psd_freqs, psd_raw = compute_psd( - data[0, start:end], sampling_freq, int(n_freqs * 2) -) -_, psd_filtered = compute_psd( - filtered_data[0, start:end], sampling_freq, int(n_freqs * 2) -) -_, psd_artefact_free = compute_psd( - artefact_free[0, start:end], sampling_freq, int(n_freqs * 2) -) +n_points = int((sampling_freq / 2) * 2) +psd_freqs, psd_raw = compute_psd(data[0, start:end], sampling_freq, n_points) +_, psd_filtered = compute_psd(filtered_data[0, start:end], sampling_freq, n_points) +_, psd_artefact_free = compute_psd(artefact_free[0, start:end], sampling_freq, n_points) -axes[1].loglog( - psd_freqs, psd_raw, color="black", alpha=0.3, label="Unfiltered" -) -axes[1].loglog( - psd_freqs, psd_artefact_free, linewidth=3, label="Artefact-free" -) +axes[1].loglog(psd_freqs, psd_raw, color="black", alpha=0.3, label="Unfiltered") +axes[1].loglog(psd_freqs, psd_artefact_free, linewidth=3, label="Artefact-free") axes[1].loglog(psd_freqs, psd_filtered, label="Filtered (PyPARRM)") axes[1].legend() axes[1].set_xlabel("Log frequency (Hz)") @@ -219,16 +194,15 @@ fig.tight_layout() fig.show() -############################################################################### +######################################################################################## # Python vs. MATLAB implementation comparison # ------------------------------------------- -# Due to rounding errors, the final filtered data timeseries of PyPARRM and the -# MATLAB implementation are not perfectly identical (as shown by the failure of -# :func:`numpy.all`), however they are extremely close (as shown by the success -# of :func:`numpy.allclose`). Visual inspection of the plots further -# demonstrates this similarity. Accordingly, both implementations are suitable -# for identifying and removing stimulation artefacts from electrophysiological -# recordings. +# Due to rounding errors, the final filtered data timeseries of PyPARRM and the MATLAB +# implementation are not perfectly identical (as shown by the failure of +# :func:`numpy.all`), however they are extremely close (as shown by the success of +# :func:`numpy.allclose`). Visual inspection of the plots further demonstrates this +# similarity. Accordingly, both implementations are suitable for identifying and +# removing stimulation artefacts from electrophysiological recordings. # %% @@ -240,10 +214,7 @@ # main plot axis.plot( - times, - matlab_filtered[0, start:end], - linewidth=3, - label="Filtered (MATLAB PARRM)", + times, matlab_filtered[0, start:end], linewidth=3, label="Filtered (MATLAB PARRM)" ) axis.plot(times, filtered_data[0, start:end], label="Filtered (PyPARRM)") axis.legend() @@ -254,9 +225,7 @@ axis.set_ylim(ylim[0], ylim[1] * 3) # inset plot -inset_axis.plot( - times[:50], matlab_filtered[0, start : start + 50], linewidth=3 -) +inset_axis.plot(times[:50], matlab_filtered[0, start : start + 50], linewidth=3) inset_axis.plot(times[:50], filtered_data[0, start : start + 50]) axis.indicate_inset_zoom(inset_axis, edgecolor="black", alpha=0.4) @@ -270,9 +239,9 @@ f"{np.allclose(filtered_data, matlab_filtered)}" ) -############################################################################### +######################################################################################## # References -# ----------------------------------------------------------------------------- +# ---------- # .. footbibliography:: # %% diff --git a/pyproject.toml b/pyproject.toml index ddf5dd3..ae0352a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -3,7 +3,7 @@ build-backend = "hatchling.build" requires = ["hatchling"] [project] -authors = [{email = "t.s.binns@outlook.com", name = "Thomas Samuel Binns"}] +authors = [{email = "t.s.binns@outlook.com", name = "Thomas S. Binns"}] classifiers = [ "License :: OSI Approved :: MIT License", "Operating System :: OS Independent", diff --git a/src/pyparrm/_utils/_plotting.py b/src/pyparrm/_utils/_plotting.py index 8b2aa51..fec4529 100644 --- a/src/pyparrm/_utils/_plotting.py +++ b/src/pyparrm/_utils/_plotting.py @@ -1,8 +1,9 @@ """Tools for plotting results.""" # Author(s): -# Thomas Samuel Binns | github.com/tsbinns +# Thomas S. Binns | github.com/tsbinns +from copy import deepcopy from multiprocessing import cpu_count from matplotlib import pyplot as plt @@ -108,10 +109,10 @@ def _check_sort_init_inputs( ) -> None: """Check and sort init. inputs.""" assert parrm._period is not None, ( - "PyPARRM Internal Error: `_ParamSelection` should only be called if the " + "PyPARRM Internal Error: `_ExploreParams` should only be called if the " "period has been estimated. Please contact the PyPARRM developers." ) - self.parrm = parrm + self.parrm = deepcopy(parrm) self.parrm._verbose = False # time_range @@ -297,7 +298,7 @@ def _initialise_plot(self) -> None: layout="constrained", ) axes["lower inner"].remove() - self.figure.set_constrained_layout_pads(w_pad=0.05, h_pad=0.1) + self.figure.get_layout_engine().set(w_pad=0.05, h_pad=0.1) self.figure.suptitle("placeholder\n") # smaller title, less space self.figure.canvas.draw() # set basic layout self._update_suptitle() # set bigger title, but don't redraw! diff --git a/src/pyparrm/_utils/_power.py b/src/pyparrm/_utils/_power.py index 8e4835f..82ed79f 100644 --- a/src/pyparrm/_utils/_power.py +++ b/src/pyparrm/_utils/_power.py @@ -1,7 +1,7 @@ """Tools for computing signal power.""" # Author(s): -# Thomas Samuel Binns | github.com/tsbinns +# Thomas S. Binns | github.com/tsbinns import numpy as np from scipy.fft import fft, fftfreq # faster than numpy for real-valued inputs diff --git a/src/pyparrm/parrm.py b/src/pyparrm/parrm.py index 8ade189..9403d7b 100644 --- a/src/pyparrm/parrm.py +++ b/src/pyparrm/parrm.py @@ -1,7 +1,7 @@ """Tools for fitting PARRM filters to data.""" # Author(s): -# Thomas Samuel Binns | github.com/tsbinns +# Thomas S. Binns | github.com/tsbinns from multiprocessing import cpu_count @@ -634,9 +634,9 @@ def _fit_waves_to_data( def explore_filter_params( self, - time_range: list[int] | None = None, + time_range: list[int | float] | None = None, time_res: int | float = 0.01, - freq_range: list[int] | None = None, + freq_range: list[int | float] | None = None, freq_res: int | float = 5.0, n_jobs: int = 1, ) -> None: