Note
Go to the end to download the full example code
2. Cycle-by-cycle algorithm¶
Demonstration of the cycle-by-cycle approach.
In the last tutorial notebook, I described the conventional approach for analyzing time-varying properties of neural oscillations, and in this notebook, we will go over our alternative approach. The fundamental goal of this approach is to characterize neural oscillations directly in the time domain. However, this is not straightforward because it attempts to extract the properties of the oscillatory component, despite the large amount of noise. Specifically, there are two very difficult problems:
What are the features of the oscillation? How do they vary over time?
During what times is the oscillation present?
The cycle-by-cycle approach deploys a few strategies to approach these questions. As its name indicates, this algorithm segments the signal into individual cycles and then analyzes their features separately from one another. Normally, some preprocessing is recommended to aid in localizing peaks and troughs (eliminating high frequency power that mostly do not comprise the oscillator of interest). Additionally, a burst detection approach is applied to define the segments of the signal to be analyzed for their oscillatory properties.
During this process (as with all data analyses), it is important to be aware if the data is being processed appropriately. As signal processing is complicated, it is very beneficial to visualize the measured features along with the raw data to assure they make sense.
0. Preprocess signal¶
A crucial part of the cycle-by-cycle approach is the ability to localize the peaks and troughs of the oscillation. Therefore, some preprocessing of the signal is often useful in order to make these extrema more apparent, i.e. isolate the oscillation component and minimize the nonoscillatory components. One effective way of doing this is by applying a lowpass filter. The choice of cutoff frequency is very important. The cutoff frequency should not be low enough in order to remove high frequency “noise” that interferes with extrema localization but not so low that it deforms the shape of the oscillation of interest. In order to assess this, the user should plot the filtered signal in comparison to the original signal.
import numpy as np
import pandas as pd
from neurodsp.filt import filter_signal
from neurodsp.plts import plot_time_series
from neurodsp.sim import sim_combined
from bycycle import Bycycle
from bycycle.cyclepoints import find_extrema, find_zerox
from bycycle.cyclepoints.zerox import find_flank_zerox
from bycycle.plts import plot_cyclepoints_array
from bycycle.utils.download import load_bycycle_data
pd.options.display.max_columns = 10
# Load data
sig = load_bycycle_data('ca1.npy', folder='data')
fs = 1250
# Filter settings
f_theta = (4, 10)
f_lowpass = 30
n_seconds_filter = .1
# Lowpass filter
sig_low = filter_signal(sig, fs, 'lowpass', f_lowpass,
n_seconds=n_seconds_filter, remove_edges=False)
# Plot signal
times = np.arange(0, len(sig)/fs, 1/fs)
xlim = (2, 5)
tidx = np.logical_and(times >= xlim[0], times < xlim[1])
plot_time_series(times[tidx], [sig[tidx], sig_low[tidx]], colors=['k', 'k'], alpha=[.5, 1], lw=2)
1. Localize peaks and troughs¶
In order to characterize the oscillation, it is useful to know the precise times of peaks and troughs. For one, this will allow us to compute the periods and rise-decay symmetries of the individual cycles. To do this, the signal is first narrow-bandpass filtered in order to estimate “zero-crossings.” Then, in between these zerocrossings, the absolute maxima and minima are found and labeled as the peaks and troughs, respectively.
# Narrowband filter signal
n_seconds_theta = .75
sig_narrow = filter_signal(sig, fs, 'bandpass', f_theta,
n_seconds=n_seconds_theta, remove_edges=False)
# Find rising and falling zerocrossings (narrowband)
rise_xs = find_flank_zerox(sig_narrow, 'rise')
decay_xs = find_flank_zerox(sig_narrow, 'decay')
# Find peaks and troughs (this function also does the above)
peaks, troughs = find_extrema(sig_low, fs, f_theta,
filter_kwargs={'n_seconds':n_seconds_theta})
plot_cyclepoints_array(sig_low, fs, peaks=peaks, troughs=troughs, xlim=(12, 15))
Note the filter characteristics used in the process of finding peaks and troughs¶
# Plot frequency response of bandpass filter
sig_filt = filter_signal(sig, fs, 'bandpass', (4, 10), n_seconds=.75, plot_properties=True)
2. Localize rise and decay midpoints¶
In addition to localizing the peaks and troughs of a cycle, we also want to get more information about the rise and decay periods. For instance, these flanks may have deflections if the peaks or troughs are particularly sharp. In order to gauge a dimension of this, we localize midpoints for each of the rise and decay segments. These midpoints are defined as the times at which the voltage crosses halfway between the adjacent peak and trough voltages. If this threshold is crossed multiple times, then the median time is chosen as the flank midpoint. This is not perfect; however, this is rare, and most of these cycles should be removed by burst detection.
Note: Plotting midpoints and extrema may also be performed using the dataframe output from
compute_features() with the plot_cyclepoints() function.
rises, decays = find_zerox(sig_low, peaks, troughs)
plot_cyclepoints_array(sig_low, fs, xlim=(13, 14), peaks=peaks, troughs=troughs,
rises=rises, decays=decays)
3. Compute features of each cycle¶
After these 4 points of each cycle are localized, we compute some simple statistics for each
cycle. The main cycle-by-cycle object, Bycycle, has a df_features dataframe
attribute containing cycle features and sample locations of cyclepoints in the signal. Each entry
or row in either dataframe is a cycle and each column is a property of that cycle (see table below).
The four main features are:
amplitude (volt_amp) - average voltage change of the rise and decay
period (period) - time between consecutive troughs (or peaks, if default is changed)
rise-decay symmetry (time_rdsym) - fraction of the period in the rise period
peak-trough symmetry (time_ptsym) - fraction of the period in the peak period
Note that a warning appears here because no burst detection parameters are provided. This is addressed in section #4.
bm = Bycycle()
bm.fit(sig, fs, f_theta)
/Users/ryanhammonds/projects/bycycle/bycycle/objs/fit.py:31: UserWarning:
No burst detection thresholds are provided. This is not recommended. Please
inspect your data and choose appropriate parameters for 'thresholds'.
Default burst detection parameters are likely not well suited for your
desired application.
warnings.warn("""
bm.df_features
# Dataframe columns are directly accessible via attributes
bm.volt_amp
array([2.3585, 3.2485, 2.3595, 1.298 , 1.8865, 2.5945, 2.4725, 2.262 ,
2.486 , 1.9365, 2.834 , 2.263 , 2.724 , 2.6 , 2.0375, 2.9975,
2.628 , 3.254 , 3.2795, 2.5465, 2.2205, 2.7695, 2.3775, 1.9505,
1.741 , 2.0585, 2.9705, 2.25 , 2.2135, 1.7865, 1.6535, 2.5215,
2.1585, 2.915 , 2.751 , 3.5615, 2.4685, 2.572 , 2.4155, 2.1445,
2.1505, 2.2405, 2.5145, 2.5305, 1.718 , 2.2285, 2.5425, 2.709 ,
2.4675, 2.844 , 1.691 , 2.289 , 2.353 , 2.1625, 2.021 , 2.3845,
2.1105, 1.6185, 2.5705, 3.103 , 2.8055, 3.046 , 2.9755, 3.0415,
2.598 , 2.933 , 2.7675, 2.3675, 2.398 , 2.3805, 2.311 , 2.4265,
2.5105, 2.94 , 2.24 , 2.0995, 2.1035, 2.253 , 2.0635, 2.028 ,
2.4255, 2.5195, 2.3495, 1.5075, 1.8275, 1.8425, 2.151 , 2.7045,
1.9765, 2.604 , 2.3405, 2.78 , 2.3375, 2.2315, 1.67 , 1.8525,
2.8815, 1.75 , 2.5065, 2.825 , 2.1985, 3.097 , 2.7375, 3.2905,
3.1635, 2.8975, 2.8155, 2.996 , 2.9965, 2.7015, 3.04 , 2.305 ,
2.244 , 2.5645, 2.295 , 2.0615, 2.345 , 2.455 , 2.5535, 3.3205,
2.273 , 2.7745, 2.834 , 1.9065, 2.1815, 1.665 , 2.5985, 2.4485,
2.383 , 2.634 , 2.195 , 2.012 , 2.144 , 2.464 , 2.01 , 1.8145,
2.1355, 2.7455, 2.9545, 2.1015, 3.263 , 2.4565, 1.9715, 3.6035,
2.8145, 2.985 , 2.1055, 2.685 , 2.389 , 2.452 , 2.086 , 2.2215,
3.3515, 2.1385, 2.077 , 2.327 , 2.416 , 1.4395, 2.3335, 2.528 ,
1.8405, 2.6875, 2.7815, 2.2745, 1.758 , 3.3585, 3.2585, 2.0315,
2.1295, 2.112 , 2.922 , 3.0285, 2.693 , 2.0835, 2.9245, 2.425 ,
2.396 , 2.471 , 2.2825, 1.992 , 2.3655, 2.8855, 2.7055, 2.524 ,
2.078 , 2.096 , 1.869 , 2.507 , 2.4005, 1.934 , 2.7115, 2.2135,
2.148 , 2.941 , 2.308 , 2.448 , 2.4435, 2.191 , 2.0785, 2.275 ,
2.6175, 2.652 , 2.626 , 1.1475, 2.5735, 2.4535, 2.5035, 1.854 ,
1.7325, 1.975 , 2.03 , 2.3775, 2.178 , 3.593 , 3.3015, 1.8365,
2.4115, 2.6385, 2.219 , 2.809 , 1.9085, 2.8425, 2.4405, 1.769 ,
2.0855, 2.629 , 2.198 , 2.6595, 2.448 , 2.3225, 2.5055, 2.3505,
2.2105, 2.2345, 2.106 , 2.1225, 2.58 , 3.141 , 1.869 , 1.4585,
2.3895, 2.794 , 2.583 , 2.4005, 3.0155, 2.729 , 2.6875, 3.053 ,
3.532 , 3.4045, 3.291 , 2.8395, 2.6235, 2.3375, 2.365 , 3.0125,
2.559 , 2.535 , 2.5665, 2.1485, 2.425 , 3.0695, 2.4675, 2.575 ,
2.451 , 2.3355, 3.143 , 2.936 , 3.2855, 2.9005, 2.731 , 2.9885,
2.8745, 3.2725, 3.6055, 2.6855, 3.295 , 2.685 , 2.654 , 1.973 ,
3.424 , 2.143 , 3.111 , 3.16 , 2.274 , 2.483 , 2.351 , 2.3425,
2.07 , 2.444 , 2.5625, 1.7155, 2.9495, 2.6025, 2.276 , 3.2865,
2.463 , 4.9315, 2.2645, 2.7765, 2.106 , 3.0895, 2.264 , 3.247 ,
2.0285, 2.3405, 2.4605, 2.297 , 2.743 , 2.707 , 2.7485, 2.4155,
2.6865, 2.96 , 2.918 , 2.3305, 2.698 , 2.568 , 2.8555, 2.297 ,
2.866 , 2.954 , 2.595 , 2.71 , 2.4545, 2.2445, 2.45 , 2.9615,
2.294 , 2.2655, 1.9505, 1.69 , 2.833 , 2.3775, 1.7635, 2.594 ,
3.058 , 2.45 , 2.158 , 2.251 , 3.2435, 2.6345, 3.519 , 2.923 ,
3.1255, 2.2635, 2.166 , 2.6145, 2.6945, 2.3325, 2.859 , 2.565 ,
2.5275, 2.6815, 2.632 , 2.5675, 2.3695, 2.365 , 2.702 , 3.2545,
2.64 , 2.725 , 2.7845, 3.1965, 2.4835, 2.6015, 2.974 , 2.8825,
2.792 , 2.573 , 2.2875, 3.072 , 2.2 , 2.1365, 2.432 , 2.951 ,
2.2995, 2.3255, 2.901 , 3.004 , 2.3845, 2.817 , 2.5265, 2.759 ,
2.053 , 2.8505, 2.5425, 2.056 , 2.1595, 2.9675, 3.1405, 3.518 ,
2.9615, 2.6845, 2.5545, 2.0325, 2.548 , 2.4265, 2.236 , 2.4585,
2.103 , 2.632 , 2.055 , 1.8525, 2.645 , 2.356 , 3.1205, 2.4585,
2.8445, 2.715 , 2.7055, 3.0595, 3.057 , 3.026 , 2.6045, 2.6185,
2.68 , 2.351 , 2.774 , 3.204 , 2.0665, 2.379 , 2.3305, 1.9235,
2.688 , 2.02 , 2.7655, 2.333 , 2.2695, 2.642 , 2.578 , 1.6665,
1.662 , 2.585 , 2.6235, 2.052 , 2.647 , 1.9835, 2.3125, 2.2435,
2.8685, 2.8355, 2.479 , 2.686 , 2.6615, 2.343 , 2.839 , 2.929 ,
2.815 , 3.166 , 2.876 , 2.511 , 1.806 , 2.3045, 2.4035, 2.2395,
2.047 , 2.3925, 3.1375, 2.519 , 2.7135, 2.927 , 2.26 , 1.511 ,
2.2075, 2.331 , 2.8535, 1.728 ])
4. Determine parts of signal in oscillatory burst¶
Note above that the signal is segmented into cycles and the dataframe provides properties for each segment of the signal. However, if no oscillation is apparent in the signal at a given time, the properties for these “cycles” are meaningless. Therefore, it is useful to have a binary indicator for each cycle that indicates whether the cycle being analyzed is truly part of an oscillatory burst or not. Recently, significant interest has emerged in detecting bursts in signals and analyzing their properties (see e.g. Feingold et al., PNAS, 2015). Nearly all efforts toward burst detection relies on amplitude thresholds, but this can be disadvantageous because these algorithms will behave very differently on signals where oscillations are common versus rare.
In our approach, we employ an alternative technique for burst detection. There are 3 thresholds that need to be met in order for a cycle to be classified as part of an oscillatory burst.
amplitude consistency - consecutive rises and decays should be comparable in magnitude.
The amplitude consistency of a cycle is equal to the maximum relative difference between rises and decay amplitudes across all pairs of adjacent rises and decays that include one of the flanks in the cycle (3 pairs)
e.g. if a rise is 10mV and a decay is 7mV, then its amplitude consistency is 0.7.
period consistency - consecutive cycles should be comparable in duration
The period consistency is equal to the maximu relative difference between all pairs of adjacent periods that include the cycle of interest (2 pairs: current + previous cycles and current + next cycles)
e.g. if the previous, current, and next cycles have periods 60ms, 100ms, and 120ms, respectively, then the period consistency is min(60/100, 100/120) = 0.6.
monotonicity - the rise and decay flanks of the cycle should be mostly monotonic
The monotonicity is the fraction of samples that the instantaneous derivative (numpy.diff) is consistent with the direction of the flank.
e.g. if in the rise, the instantaneous derivative is 90% positive, and in the decay, the instantaneous derivative is 80% negative, then the monotonicity of the cycle would be 0.85 ((0.9+0.8)/2)
Below, we load a simulated signal and then define 3 sets of thresholds ranging from liberal to conservative.
Load a simulated signal and apply a lowpass filter¶
# Simulate a signal
n_seconds = 10
fs = 1000 # Sampling rate
f_alpha = (8, 12)
components = {'sim_bursty_oscillation': {'freq': 10, 'enter_burst': .1, 'leave_burst': .1},
'sim_powerlaw': {'f_range': (2, None)}}
sig = sim_combined(n_seconds, fs, components=components, component_variances=(5, 1))
# Apply a lowpass filter to remove high frequency power that interferes with extrema localization
sig = filter_signal(sig, fs, 'lowpass', 30, n_seconds=.2, remove_edges=False)
Visualizing burst detection settings¶
Below, we visualize how the burst detector determined which cycles were part of an oscillatory burst. The top plot shows a markup of the time series. The portions of the signal in red were determined to be parts of bursts. Signals in black were not part of bursts. Magenta and cyan dots denote detected peaks and troughs, respectively. Highlights indicate cycles marked as not part of a burst because they did not meet certain thresholds:
red highlight: amplitude consistency threshold violation
yellow highlight: period consistency threshold violation
green highlight: monotonicity threshold violation The plots below show the relevant features for each cycle as well as the threshold (dotted lines), where we can see the highlights appear if the features went below the threshold.
Note there is an optional “band amplitude fraction” threshold. This is currently unused (set to 0) , but is present in case users want to add an amplitude threshold to this algorithm.
Burst detection settings: too liberal¶
The following burst detection thresholds (defined in burst_kwargs) are too low, so some portions of the signal that do not have much apparent oscillatory burst are still labeled as if they do.
thresholds = {
'amp_fraction': 0,
'amp_consistency': .2,
'period_consistency': .45,
'monotonicity': .7,
'min_n_cycles': 3
}
bm = Bycycle(thresholds=thresholds)
bm.fit(sig, fs, f_alpha)
bm.plot(figsize=(16, 3))
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
Burst detection settings: too conservative¶
These new burst detection thresholds seem to be set too high (too strict) as the algorithm is not able to detect the bursts that are present.
thresholds = {
'amp_fraction': 0,
'amp_consistency': .75,
'period_consistency': .7,
'monotonicity': .9,
'min_n_cycles': 3
}
bm = Bycycle(thresholds=thresholds)
bm.fit(sig, fs, f_alpha)
bm.plot(figsize=(16, 3))
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
More appropriate burst detection thresholds¶
The conservative thresholds were then lowered, and we can see now that the algorithms correctly identifies parts of the 3 bursting periods. Therefore, for a signal with this level of noise, we expect these parameters to be pretty good.
Notice that adding a small amplitude fraction threshold (e.g. 0.3) helps remove some false positives that may occur, like that around 1.5 seconds. Consisitency features for non-bursting cycles on the edges of bursts may be recomputed to reduce false negatives when entering or exiting a burst.
thresholds = {
'amp_fraction': .3,
'amp_consistency': .4,
'period_consistency': .5,
'monotonicity': .8,
'min_n_cycles': 3
}
bm = Bycycle(thresholds=thresholds)
bm.fit(sig, fs, f_alpha)
bm.plot(figsize=(16, 3))
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
/Users/ryanhammonds/projects/bycycle/.env/lib/python3.11/site-packages/neurodsp/plts/style.py:69: MatplotlibDeprecationWarning: MarkerStyle(None) is deprecated since 3.6; support will be removed two minor releases later. Use MarkerStyle('') to construct an empty MarkerStyle.
line.set(**{style : next(values)})
Total running time of the script: ( 0 minutes 2.435 seconds)