diff --git a/neurodsp/plts/rhythm.py b/neurodsp/plts/rhythm.py index d7ae774e..35cb59d3 100644 --- a/neurodsp/plts/rhythm.py +++ b/neurodsp/plts/rhythm.py @@ -24,6 +24,7 @@ def plot_swm_pattern(pattern, ax=None, **kwargs): -------- Plot the average pattern from a sliding window matching analysis: + >>> import numpy as np >>> from neurodsp.sim import sim_combined >>> from neurodsp.rhythm import sliding_window_matching >>> sig = sim_combined(n_seconds=10, fs=500, @@ -31,7 +32,8 @@ def plot_swm_pattern(pattern, ax=None, **kwargs): ... 'sim_bursty_oscillation': {'freq': 20, ... 'enter_burst': .25, ... 'leave_burst': .25}}) - >>> avg_window, _, _ = sliding_window_matching(sig, fs=500, win_len=0.05, win_spacing=0.5) + >>> windows, _ = sliding_window_matching(sig, fs=500, win_len=0.05, win_spacing=0.5) + >>> avg_window = np.mean(windows) >>> plot_swm_pattern(avg_window) """ diff --git a/neurodsp/rhythm/lc.py b/neurodsp/rhythm/lc.py index 13ed968a..fe7e3b2f 100644 --- a/neurodsp/rhythm/lc.py +++ b/neurodsp/rhythm/lc.py @@ -23,13 +23,14 @@ def compute_lagged_coherence(sig, fs, freqs, n_cycles=3, return_spectrum=False): fs : float Sampling rate, in Hz. freqs : 1d array or list of float - If array, frequency values to estimate with morlet wavelets. + The frequency values at which to estimate lagged coherence. + If array, defines the frequency values to use. If list, define the frequency range, as [freq_start, freq_stop, freq_step]. The `freq_step` is optional, and defaults to 1. Range is inclusive of `freq_stop` value. n_cycles : float or list of float, default: 3 - Number of cycles of each frequency to use to compute lagged coherence. - If a single value, the same number of cycles is used for each frequency value. - If a list or list_like, then should be a n_cycles corresponding to each frequency. + The number of cycles to use to compute lagged coherence, for each frequency. + If a single value, the same number of cycles is used for each frequency. + If a list or list_like, there should be a value corresponding to each frequency. return_spectrum : bool, optional, default: False If True, return the lagged coherence for all frequency values. Otherwise, only the mean lagged coherence value across the frequency range is returned. @@ -87,7 +88,7 @@ def compute_lagged_coherence(sig, fs, freqs, n_cycles=3, return_spectrum=False): def lagged_coherence_1freq(sig, fs, freq, n_cycles): - """Compute the lagged coherence of a frequency using the hanning-taper FFT method. + """Compute the lagged coherence at a particular frequency. Parameters ---------- @@ -98,12 +99,17 @@ def lagged_coherence_1freq(sig, fs, freq, n_cycles): freq : float The frequency at which to estimate lagged coherence. n_cycles : float - Number of cycles at the examined frequency to use to compute lagged coherence. + The number of cycles of the given frequency to use to compute lagged coherence. Returns ------- - float + lc : float The computed lagged coherence value. + + Notes + ----- + - Lagged coherence is computed using hanning-tapered FFTs. + - The returned lagged coherence value is bound between 0 and 1. """ # Determine number of samples to be used in each window to compute lagged coherence @@ -113,20 +119,26 @@ def lagged_coherence_1freq(sig, fs, freq, n_cycles): chunks = split_signal(sig, n_samps) n_chunks = len(chunks) - # For each chunk, calculate the Fourier coefficients at the frequency of interest + # Create the window to apply to each chunk hann_window = hann(n_samps) + + # Create the frequency vector, finding the frequency value of interest fft_freqs = np.fft.fftfreq(n_samps, 1 / float(fs)) fft_freqs_idx = np.argmin(np.abs(fft_freqs - freq)) + # Calculate the Fourier coefficients across chunks for the frequency of interest fft_coefs = np.zeros(n_chunks, dtype=complex) for ind, chunk in enumerate(chunks): fourier_coef = np.fft.fft(chunk * hann_window) fft_coefs[ind] = fourier_coef[fft_freqs_idx] - # Compute the lagged coherence value + # Compute lagged coherence across data segments lcs_num = 0 for ind in range(n_chunks - 1): lcs_num += fft_coefs[ind] * np.conj(fft_coefs[ind + 1]) lcs_denom = np.sqrt(np.sum(np.abs(fft_coefs[:-1])**2) * np.sum(np.abs(fft_coefs[1:])**2)) - return np.abs(lcs_num / lcs_denom) + # Normalize the lagged coherence value + lc_val = np.abs(lcs_num / lcs_denom) + + return lc_val diff --git a/neurodsp/rhythm/swm.py b/neurodsp/rhythm/swm.py index 54bc1dd7..14cb9ece 100644 --- a/neurodsp/rhythm/swm.py +++ b/neurodsp/rhythm/swm.py @@ -1,4 +1,4 @@ -"""The sliding window matching algorithm for identifying rhythmic components of a neural signal.""" +"""The sliding window matching algorithm for identifying recurring patterns in a neural signal.""" import numpy as np @@ -8,8 +8,8 @@ ################################################################################################### @multidim() -def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, - temperature=1, window_starts_custom=None): +def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=100, + window_starts_custom=None, var_thresh=None): """Find recurring patterns in a time series using the sliding window matching algorithm. Parameters @@ -19,24 +19,23 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, fs : float Sampling rate, in Hz. win_len : float - Window length, in seconds. + Window length, in seconds. This is L in the original paper. win_spacing : float - Minimum window spacing, in seconds. - max_iterations : int + Minimum spacing between windows, in seconds. This is G in the original paper. + max_iterations : int, optional, default: 100 Maximum number of iterations of potential changes in window placement. - temperature : float - Temperature parameter. Controls probability of accepting a new window. - window_starts_custom : 1d array, optional + window_starts_custom : 1d array, optional, default: None Custom pre-set locations of initial windows. + var_thresh: float, opational, default: None + Removes initial windows with variance below a set threshold. This speeds up + runtime proportional to the number of low variance windows in the data. Returns ------- - avg_window : 1d array - The average waveform in 'sig' in the frequency 'f_range' triggered on 'trigger'. + windows : 2d array + Putative patterns discovered in the input signal. window_starts : 1d array Indices at which each window begins for the final set of windows. - costs : 1d array - Cost function value at each iteration. References ---------- @@ -44,125 +43,170 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, van der Eerden, J. (2017). Discovering recurring patterns in electrophysiological recordings. Journal of Neuroscience Methods, 275, 66-79. DOI: 10.1016/j.jneumeth.2016.11.001 - Matlab Code: https://github.com/bartgips/SWM + .. [2] Matlab Code implementation: https://github.com/bartgips/SWM Notes ----- - - Apply a highpass filter if looking at high frequency activity, so that it does - not converge on a low frequency motif. - - Parameters `win_len` and `win_spacing` should be chosen to be about the size of the - motif of interest, and the N derived should be about the number of occurrences. + - The `win_len` parameter should be chosen to be about the size of the motif of interest. + The larger this window size, the more likely the pattern to reflect slower patterns. + - The `win_spacing` parameter also determines the number of windows that are used. + - If looking at high frequency activity, you may want to apply a highpass filter, + so that the algorithm does not converge on a low frequency motif. + - This implementation is a minimal, modified version, as compared to the original + implementation in [2], which has more available options. + - This version has the following changes to speed up convergence: + + 1. Each iteration is similar to an epoch, randomly moving all windows in + random order. The original implementation randomly selects windows and + does not guarantee even resampling. + 2. New window acceptance is determined via increased correlation coefficients + and reduced ivariance across windows. + 3. Phase optimization / realignment to escape local minima. + Examples -------- Search for reoccuring patterns using sliding window matching in a simulated beta signal: >>> from neurodsp.sim import sim_combined - >>> sig = sim_combined(n_seconds=10, fs=500, - ... components={'sim_powerlaw': {'f_range': (2, None)}, - ... 'sim_bursty_oscillation': {'freq': 20, - ... 'enter_burst': .25, - ... 'leave_burst': .25}}) - >>> avg_window, window_starts, costs = sliding_window_matching(sig, fs=500, win_len=0.05, - ... win_spacing=0.20) + >>> components = {'sim_bursty_oscillation' : {'freq' : 20, 'phase' : 'min'}, + ... 'sim_powerlaw' : {'f_range' : (2, None)}} + >>> sig = sim_combined(10, fs=500, components=components, component_variances=(1, .05)) + >>> windows, starts = sliding_window_matching(sig, fs=500, win_len=0.05, + ... win_spacing=0.05, var_thresh=.5) """ # Compute window length and spacing in samples - win_n_samps = int(win_len * fs) - spacing_n_samps = int(win_spacing * fs) + win_len = int(win_len * fs) + win_spacing = int(win_spacing * fs) # Initialize window positions if window_starts_custom is None: - window_starts = np.arange(0, len(sig) - win_n_samps, 2 * spacing_n_samps) + window_starts = np.arange(0, len(sig) - win_len, win_spacing).astype(int) else: window_starts = window_starts_custom - n_windows = len(window_starts) - # Randomly sample windows with replacement - random_window_idx = np.random.choice(range(n_windows), size=max_iterations) + windows = np.array([sig[start:start + win_len] for start in window_starts]) + + # Compute new window bounds + lower_bounds, upper_bounds = _compute_bounds(window_starts, win_spacing, 0, len(sig) - win_len) + + # Remove low variance windows to speed up runtime + if var_thresh is not None: + + thresh = np.array([np.var(sig[loc:loc + win_len]) > var_thresh for loc in window_starts]) + + windows = windows[thresh] + window_starts = window_starts[thresh] + lower_bounds = lower_bounds[thresh] + upper_bounds = upper_bounds[thresh] + + # Modified SWM procedure + window_idxs = np.arange(len(windows)).astype(int) + + corrs, variance = _compute_cost(sig, window_starts, win_len) + mae = np.mean(np.abs(windows - windows.mean(axis=0))) - # Calculate initial cost - costs = np.zeros(max_iterations) - costs[0] = _compute_cost(sig, window_starts, win_n_samps) + for _ in range(max_iterations): - for iter_num in range(1, max_iterations): + # Randomly shuffle order of windows + np.random.shuffle(window_idxs) - # Pick a random window position to randomly replace with a - # new window to improve cross-window similarity - window_idx_replace = random_window_idx[iter_num] + for win_idx in window_idxs: - # Find a new allowed position for the window - window_starts_temp = np.copy(window_starts) - window_starts_temp[window_idx_replace] = _find_new_window_idx( - window_starts, spacing_n_samps, len(sig) - win_n_samps) + # Find a new, random window start + _window_starts = window_starts.copy() + _window_starts[win_idx] = np.random.choice(np.arange(lower_bounds[win_idx], + upper_bounds[win_idx] + 1)) - # Calculate the cost & the change in the cost function - cost_temp = _compute_cost(sig, window_starts_temp, win_n_samps) - delta_cost = cost_temp - costs[iter_num - 1] + # Accept new window if correlation increases and variance decreases + _corrs, _variance = _compute_cost(sig, _window_starts, win_len) - # Calculate the acceptance probability - p_accept = np.exp(-delta_cost / float(temperature)) + if _corrs[win_idx].sum() > corrs[win_idx].sum() and _variance < variance: - # Accept update to J with a certain probability - if np.random.rand() < p_accept: + corrs = _corrs.copy() + variance = _variance + window_starts = _window_starts.copy() + lower_bounds, upper_bounds = _compute_bounds(\ + window_starts, win_spacing, 0, len(sig) - win_len) - # Update costs & windows - costs[iter_num] = cost_temp - window_starts = window_starts_temp + # Phase optimization + _window_starts = window_starts.copy() - else: + for shift in np.arange(-int(win_len/2), int(win_len/2)): - # Update costs - costs[iter_num] = costs[iter_num - 1] + _starts = _window_starts + shift - # Calculate average window - avg_window = np.zeros(win_n_samps) - for w_ind in range(n_windows): - avg_window = avg_window + sig[window_starts[w_ind]:window_starts[w_ind] + win_n_samps] - avg_window = avg_window / float(n_windows) + # Skip windows shifts that are out-of-bounds + if (_starts[0] < 0) or (_starts[-1] > len(sig) - win_len): + continue - return avg_window, window_starts, costs + _windows = np.array([sig[start:start + win_len] for start in _starts]) + + _mae = np.mean(np.abs(_windows - _windows.mean(axis=0))) + + if _mae < mae: + window_starts = _starts.copy() + windows = _windows.copy() + mae = _mae + + lower_bounds, upper_bounds = _compute_bounds(\ + window_starts, win_spacing, 0, len(sig) - win_len) + + return windows, window_starts def _compute_cost(sig, window_starts, win_n_samps): - """Compute the cost, which is proportional to the difference between pairs of windows.""" + """Compute the cost, as corrleation coefficients and variance across windows. + + Parameters + ---------- + sig : 1d array + Time series. + window_starts : list of int + The list of window start definitions. + win_n_samps : int + The length of each window, in samples. - # Get all windows and z-score them - n_windows = len(window_starts) - windows = np.zeros((n_windows, win_n_samps)) + Returns + ------- + corrs: 2d array + Window correlation matrix. + variance: float + Sum of the variance across windows. + """ - for ind, window in enumerate(window_starts): - temp = sig[window:window_starts[ind] + win_n_samps] - windows[ind] = (temp - np.mean(temp)) / np.std(temp) + windows = np.array([sig[start:start + win_n_samps] for start in window_starts]) - # Calculate distances for all pairs of windows - dists = [] - for ind1 in range(n_windows): - for ind2 in range(ind1 + 1, n_windows): - window_diff = windows[ind1] - windows[ind2] - dist_temp = np.sum(window_diff**2) / float(win_n_samps) - dists.append(dist_temp) + corrs = np.corrcoef(windows) - # Calculate cost function, which is the average difference, roughly - cost = np.sum(dists) / float(2 * (n_windows - 1)) + variance = windows.var(axis=1).sum() - return cost + return corrs, variance -def _find_new_window_idx(window_starts, spacing_n_samps, n_samp, tries_limit=1000): - """Find a new sample for the starting window.""" +def _compute_bounds(window_starts, win_spacing, start, end): + """Compute bounds on a new window. - for n_try in range(tries_limit): + Parameters + ---------- + window_starts : list of int + The list of window start definitions. + win_spacing : float + Minimum spacing between windows, in seconds. + start, end : int + Start and end edges for the possible window. - # Generate a random sample & check how close it is to other window starts - new_samp = np.random.randint(n_samp) - dists = np.abs(window_starts - new_samp) + Returns + ------- + lower_bounds, upper_bounds : 1d array + Computed upper and lower bounds for the window position. + """ - if np.min(dists) > spacing_n_samps: - break + lower_bounds = window_starts[:-1] + win_spacing + lower_bounds = np.insert(lower_bounds, 0, start) - else: - raise RuntimeError('SWM algorithm has difficulty finding a new window. \ - Try increasing the spacing parameter.') + upper_bounds = window_starts[1:] - win_spacing + upper_bounds = np.insert(upper_bounds, len(upper_bounds), end) - return new_samp + return lower_bounds, upper_bounds diff --git a/neurodsp/tests/rhythm/test_swm.py b/neurodsp/tests/rhythm/test_swm.py index 53e3260c..d156c7c9 100644 --- a/neurodsp/tests/rhythm/test_swm.py +++ b/neurodsp/tests/rhythm/test_swm.py @@ -11,12 +11,12 @@ def test_sliding_window_matching(tsig): win_len, win_spacing = 1, 0.5 - pattern, starts, costs = sliding_window_matching(tsig, FS, win_len, win_spacing) - assert pattern.shape[-1] == int(FS * win_len) + windows, starts = sliding_window_matching(tsig, FS, win_len, win_spacing, var_thresh=0.1) + assert windows.shape[-1] == int(FS * win_len) def test_sliding_window_matching_2d(tsig2d): win_len, win_spacing = 1, 0.5 - pattern, starts, costs = sliding_window_matching(tsig2d, FS, win_len, win_spacing) - assert pattern.shape[-1] == int(FS * win_len) + windows, starts = sliding_window_matching(tsig2d, FS, win_len, win_spacing, var_thresh=0.1) + assert windows.shape[-1] == int(FS * win_len) diff --git a/tutorials/rhythm/plot_LaggedCoherence.py b/tutorials/rhythm/plot_LaggedCoherence.py index a82061b5..6a4a0782 100644 --- a/tutorials/rhythm/plot_LaggedCoherence.py +++ b/tutorials/rhythm/plot_LaggedCoherence.py @@ -2,7 +2,7 @@ Lagged Coherence ================ -Compute lagged coherence on neural signals. +Calculate the rhythmicity of neural signal with the lagged coherence algorithm. This tutorial primarily covers the :func:`~.compute_lagged_coherence` function. """ @@ -13,7 +13,13 @@ # # Lagged coherence is a measure to quantify the rhythmicity of neural signals. # -# For more details on the lagged coherence measure see Fransen et al., 2015, Neuroimage. +# Lagged coherence works by quantifying phase consistency between non-overlapping data fragments, +# calculated with Fourier coefficients. The consistency of the phase differences across +# epochs indexes the rhythmicity of the signal. Lagged coherence can be calculated for a +# particular frequency, and/or across a range of frequencies of interest. +# +# The lagged coherence algorithm is described in +# `Fransen et al, 2015 `_ # ################################################################################################### @@ -25,11 +31,11 @@ # Import the lagged coherence function from neurodsp.rhythm import compute_lagged_coherence -# Import simulation code for creating test data +# Import functions for simulating data from neurodsp.sim import sim_powerlaw, sim_combined from neurodsp.utils import set_random_seed, create_times -# Import utilities for loading and plotting data +# Import utilities for loading and plotting from neurodsp.utils.download import load_ndsp_data from neurodsp.plts.time_series import plot_time_series from neurodsp.plts.rhythm import plot_lagged_coherence @@ -43,27 +49,52 @@ # Simulate an Example Signal # -------------------------- # -# First, we'll simulate an example signal that starts with a burst of alpha activity, followed -# by a period of only aperiodic activity. +# First, let's start by creating some simulated data, to which we can apply lagged coherence. +# + +################################################################################################### +# Simulation Settings +# ~~~~~~~~~~~~~~~~~~~ +# +# We'll start with an example signal that starts with a burst of alpha activity, +# followed by a period of only aperiodic activity. # ################################################################################################### -# Set time and sampling rate -n_seconds_burst = 1 -n_seconds_noise = 2 +# Set the sampling rate fs = 1000 +# Set time for each segment of the simulated signal, in seconds +t_osc = 1 # oscillation +t_ap = 2 # aperiodic + +# Set the frequency of the oscillation +freq = 10 + +# Set the exponent value for the aperiodic activity +exp = -1 + # Create a times vector -times = create_times(n_seconds_burst + n_seconds_noise, fs) +times = create_times(t_osc + t_ap, fs) + +################################################################################################### +# Create the simulated data +# ~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Next, we will create our simulated data, by concatenating simulated segments of data +# with and without a rhythm. +# + +################################################################################################### # Simulate a signal component with an oscillation -components = {'sim_powerlaw' : {'exponent' : 0}, - 'sim_oscillation' : {'freq' : 10}} -s1 = sim_combined(n_seconds_burst, fs, components, [0.1, 1]) +components = {'sim_oscillation' : {'freq' : freq}, + 'sim_powerlaw' : {'exponent' : exp, 'f_range' : (1, None)}} +s1 = sim_combined(t_osc, fs, components, [1, 0.5]) -# Simulate a signal component with just noise -s2 = sim_powerlaw(n_seconds_noise, fs, 0, variance=0.1) +# Simulate a signal component with only aperiodic activity +s2 = sim_powerlaw(t_ap, fs, exp, variance=0.5) # Join signals together to approximate a 'burst' sig = np.append(s1, s2) @@ -74,41 +105,98 @@ plot_time_series(times, sig) ################################################################################################### -# Compute lagged coherence for an alpha oscillation -# ------------------------------------------------- +# Compute lagged coherence on simulated data +# ------------------------------------------ # # We can compute lagged coherence with the :func:`~.compute_lagged_coherence` function. # +################################################################################################### +# Data Preprocessing & Algorithm Settings +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The lagged coherence calculates the FFT across segments of the input data, and examines +# phase properties of frequencies of interest. As a spectral method that relies of the Fourier +# transform, the data does not have to be filtered prior to applying lagged coherence. +# +# You do have to specify which frequencies to examine the lagged coherence for, which is +# provided to the function as the `freqs` input. This input, which can specify either a list +# or a range of frequencies, indicates which frequencies to analyze with lagged coherence. +# +# An optional setting is the number of cycles to use at each frequency. This parameter +# controls the segment size used for each frequency. +# + ################################################################################################### # Set the frequency range to compute lagged coherence across f_range = (8, 12) +################################################################################################### +# Apply Lagged Coherence +# ~~~~~~~~~~~~~~~~~~~~~~ +# +# Next, we can apply the lagged coherence algorithm to get the average lagged coherence +# across the alpha range. +# + +################################################################################################### + # Compute lagged coherence lag_coh_alpha = compute_lagged_coherence(sig, fs, f_range) +################################################################################################### +# Lagged coherence result +# ~~~~~~~~~~~~~~~~~~~~~~~ +# +# The resulting lagged coherence value is bound between 0 and 1, with higher values +# indicating greater rhythmicity in the signal. +# + +################################################################################################### + # Check the resulting value print('Lagged coherence = ', lag_coh_alpha) +################################################################################################### +# +# The calculated lagged coherence value in the alpha range is high, meaning our measured +# lagged coherence value is indicating a large amount of rhythmicity across the analyzed +# frequencies. This is expected, since our simulated signal contains an alpha oscillation. +# + ################################################################################################### # Compute lagged coherence across the frequency spectrum # ------------------------------------------------------ # -# Notice that lagged coherence peaks around 10Hz (the frequency of our oscillation). +# What we calculated above was the average lagged coherence across a frequency range of +# interest, specifically the alpha range of (8, 12). # -# However, the peak is not very specific to that frequency. +# Instead of looking at the average lagged coherence across a range, we can also +# calculated the lagged coherence for each frequency, and then look at the +# the distribution of lagged coherence values. # +# To do so, we can set the `return_spectrum` parameter to be True, to indicated the function +# to return the full spectrum of lagged coherence values, as opposed to the average +# across the given range. +# + +################################################################################################### + +# Set the frequency range to compute the spectrum of LC values across +lc_range = (5, 40) ################################################################################################### # Calculate lagged coherence across a frequency range -lag_coh_by_f, freqs = compute_lagged_coherence(sig, fs, (5, 40), - return_spectrum=True) +lag_coh_by_f, freqs = compute_lagged_coherence(sig, fs, lc_range, return_spectrum=True) ################################################################################################### # -# You can plot the lagged coherence results with :func:`~.plot_lagged_coherence`. +# Our outputs for lagged coherence are now a vector of lagged coherence values, as well +# as a vector of the frequencies that they correspond to. +# +# To visualize this result, we can use the :func:`~.plot_lagged_coherence` function. # ################################################################################################### @@ -117,34 +205,57 @@ plot_lagged_coherence(freqs, lag_coh_by_f) ################################################################################################### -# Compute lagged coherence for time segments with and without burst -# ----------------------------------------------------------------- # -# Note that lagged coherence is greater when analyzing a neural signal that has a burst in -# the frequency range of interest, compared to a signal that does not have an oscillation. +# In these results, the lagged coherence peaks around 10 Hz. This is expected, as that is the +# frequency of the oscillation that we simulated, and the signal contains no other rhythms. +# +# Notice, however, that the peak around 10 Hz is not very specific to that frequency (it's +# quite broad). This reflects the frequency resolution of the measure. +# +# The frequency resolution is controlled by the `n_cycles` parameter. You can explore how +# the spectrum of lagged coherence values varies as you chance the `n_cycles` input. +# + +################################################################################################### +# Compute lagged coherence for segments with and without bursts +# ------------------------------------------------------------- +# +# Another factor we may want to keep in mind is that we are analyzing a signal with a +# bursty oscillation, and averaging the lagged coherence across the entire time range. +# +# To examine how the measure varies when the oscillation is and is not present, we can +# restrict our analysis to particular segments of the signal. +# +# In the following, we will apply lagged coherence separately to the bursting and +# non-bursting segments of the data. # ################################################################################################### -# Calculate coherence for data with the burst - the 1st second of data -lag_coh_burst = compute_lagged_coherence(sig[0:fs], fs, f_range) +# Calculate coherence for data segment with the oscillation present +lc_burst = compute_lagged_coherence(sig[0:t_osc*fs], fs, f_range) + +# Calculate coherence for data segment without the alpha present burst +lc_noburst = compute_lagged_coherence(sig[t_osc*fs:2*fs*t_ap], fs, f_range) -# Calculate coherence for data without the burst - the 2nd second of data -lag_coh_noburst = compute_lagged_coherence(sig[fs:2*fs], fs, f_range) +################################################################################################### -print('Lagged coherence, bursting = ', lag_coh_burst) -print('Lagged coherence, not bursting = ', lag_coh_noburst) +print('Lagged coherence, bursting = ', lc_burst) +print('Lagged coherence, not bursting = ', lc_noburst) ################################################################################################### # Compute lagged coherence of an example neural signal # ---------------------------------------------------- # -# Next, let's apply lagged coherence to some real data. +# Finally, let's apply the lagged coherence algorithm to some real data. +# +# First we can load and plot a segment of real data, which in this case is +# a segment of ECoG data with a beta oscillation. # ################################################################################################### -# Download, if needed, and load example data files +# Download, if needed, and load example data file sig = load_ndsp_data('sample_data_1.npy', folder='data') sig_filt_true = load_ndsp_data('sample_data_1_filt.npy', folder='data') @@ -157,13 +268,33 @@ # Plot example signal plot_time_series(times, sig) +################################################################################################### +# +# Now let's apply lagged coherence to the loaded data. In this data, we suspect rhythmicity +# in the beta range, so that is the range that we will examine with lagged coherence. +# + ################################################################################################### # Set the frequency range to compute lagged coherence across -f_range = (13, 30) +beta_range = (13, 30) -# Compute lagged coherence -lag_coh_beta = compute_lagged_coherence(sig, fs, f_range) +# Compute lagged coherence across the beta range in the real data +lc_betas, freqs_beta = compute_lagged_coherence(sig, fs, beta_range, return_spectrum=True) + +################################################################################################### + +# Plot the distribution of lagged coherence values in the real data +plot_lagged_coherence(freqs_beta, lc_betas) -# Check lagged coherence result -print('Lagged coherence = ', lag_coh_beta) +################################################################################################### +# Concluding Notes +# ~~~~~~~~~~~~~~~~ +# +# In the above, we can see a pattern of lagged coherence across the examined range, that +# is consistent with beta rhythmicity. +# +# To further explore this, we might want to examine the robustness by trying different +# values for `n_cycles`, and/or comparing the this result to peaks in the power spectrum, +# and/or with burst detection results. +# diff --git a/tutorials/rhythm/plot_SlidingWindowMatching.py b/tutorials/rhythm/plot_SlidingWindowMatching.py index 9c07c479..8343385a 100644 --- a/tutorials/rhythm/plot_SlidingWindowMatching.py +++ b/tutorials/rhythm/plot_SlidingWindowMatching.py @@ -2,7 +2,7 @@ Sliding Window Matching ======================= -Find recurrent patterns in a neural signal using Sliding Window Matching. +Find recurring patterns in neural signals using Sliding Window Matching. This tutorial primarily covers the :func:`~.sliding_window_matching` function. """ @@ -11,16 +11,29 @@ # Overview # -------- # -# This notebook shows how to use sliding window matching (SWM) for identifying recurring -# patterns in a neural signal, like the shape of an oscillatory waveform. +# Non-periodic or non-sinusoidal properties can be difficult to assess in frequency domain +# methods. To try and address this, the sliding window matching (SWM) algorithm has been +# proposed for detecting and measuring recurring, but unknown, patterns in time series data. +# Patterns of interest may be transient events, and/or the waveform shape of neural oscillations. # -# For more details on Sliding Window Matching see Gips et al., 2017, J Neuro Methods. +# In this example, we will explore applying the SWM algorithm to some LFP data. +# +# The SWM approach tries to find recurring patterns (or motifs) in the data, using sliding +# windows. An iterative process samples window randomly, and compares each to the average +# window. The goal is to find a selection of windows that look maximally like the average +# window, at which point the occurrences of the window have been detected, and the average +# window pattern can be examined. +# +# The sliding window matching algorithm is described in +# `Gips et al, 2017 `_ # ################################################################################################### # sphinx_gallery_thumbnail_number = 2 +import numpy as np + # Import the sliding window matching function from neurodsp.rhythm import sliding_window_matching @@ -29,6 +42,7 @@ from neurodsp.plts.rhythm import plot_swm_pattern from neurodsp.plts.time_series import plot_time_series from neurodsp.utils import set_random_seed, create_times +from neurodsp.utils.norm import normalize_sig ################################################################################################### @@ -39,61 +53,122 @@ # Load neural signal # ------------------ # +# First, we will load a segment of ECoG data, as an example time series. +# ################################################################################################### # Download, if needed, and load example data files sig = load_ndsp_data('sample_data_1.npy', folder='data') +sig = normalize_sig(sig, mean=0, variance=1) # Set sampling rate, and create a times vector for plotting fs = 1000 times = create_times(len(sig)/fs, fs) +################################################################################################### +# +# Next, we can visualize this data segment. As we can see this segment of data has +# some prominent bursts of oscillations, in this case, in the beta frequency. +# + ################################################################################################### # Plot example signal plot_time_series(times, sig) ################################################################################################### -# Apply sliding window matching to neural signal -# ---------------------------------------------- +# Apply sliding window matching +# ----------------------------- # -# The loaded neural signal has a beta oscillation, that we can attempt to analyze -# with the sliding window matching approach. -# -# We will define the window length to be about 1 cycle, which should roughly extract -# the waveform shape of the neural oscillation. +# The beta oscillation in our data segment looks like it might have some non-sinusoidal +# properties. We can investigate this with sliding window matching. # # Sliding window matching can be applied with the # :func:`~.sliding_window_matching` function. # +################################################################################################### +# Data Preprocessing +# ~~~~~~~~~~~~~~~~~~ +# +# Typically, the input signal does not have to be filtered into a band of interest to use SWM. +# +# If the goal is to characterize non-sinusoidal rhythms, you typically won't want to +# apply a filter that will smooth out the features of interest. +# +# However, if the goal is to characterize higher frequency activity, it can be useful to +# apply a highpass filter, so that the method does not converge on a lower frequency motif. +# +# In our case, the beta rhythm of interest is the most prominent, low frequency, feature of the +# data, so we won't apply a filter. +# + +################################################################################################### +# Algorithm Settings +# ~~~~~~~~~~~~~~~~~~ +# +# The SWM algorithm has some algorithm specific settings that need to be applied, including: +# +# - `win_len` : the length of the window, defined in seconds +# - `win_spacing` : the minimum distance between windows, also defined in seconds +# +# The length of the window influences the patterns that are extracted from the data. +# Typically, you want to set the window length to match the expected timescale of the +# patterns under study. +# +# For our purposes, we will define the window length to be about 1 cycle of a beta oscillation, +# which should help the algorithm to find the waveform shape of the neural oscillation. +# + ################################################################################################### # Define window length & minimum window spacing, both in seconds win_len = .055 -win_spacing = .2 +win_spacing = .055 + +################################################################################################### # Apply the sliding window matching algorithm to the time series -avg_window, window_starts, costs = sliding_window_matching(sig, fs, win_len, win_spacing, - max_iterations=500) +windows, window_starts = sliding_window_matching(sig, fs, win_len, win_spacing, var_thresh=.5) ################################################################################################### +# Examine the Results +# ~~~~~~~~~~~~~~~~~~~ +# +# What we got back from the SWM function are the calculate average window, the list +# of indices in the data of the windows, and the calculated costs for each iteration of +# the algorithm run. # -# You can plot the resulting pattern with :func:`~.plot_swm_pattern`. +# In order to visualize the resulting pattern, we can use +# :func:`~.plot_swm_pattern`. # ################################################################################################### +# Compute the average window +avg_window = np.mean(windows, 0) + # Plot the discovered pattern plot_swm_pattern(avg_window) ################################################################################################### # -# Notice that the beta cycles have sharper troughs than peaks, and the average window is -# a beta cycle with a sharp trough. +# In the above average pattern, that looks to capture a beta rhythm, we can notice some +# waveform shape of the extracted rhythm. +# + +################################################################################################### +# Concluding Notes +# ~~~~~~~~~~~~~~~~ +# +# One thing to keep in mind is that the SWM algorithm includes a random element of sampling +# and comparing the windows - meaning it is not deterministic. Because of this, results +# can change with different random seeds. # -# One thing to explore is how these results change by changing the random seed. +# To explore this, go back and change the random seed, and see how the output changes. # -# Using more data and increasing the number of iterations helps the robustness of the algorithm. +# You can also set the number of iterations that the algorithm sweeps through. Increasing +# the number of iterations, and using longer data segments, can help improve the robustness +# of the algorithm results. # diff --git a/tutorials/sim/plot_SimulateAperiodic.py b/tutorials/sim/plot_SimulateAperiodic.py index 0b73cb4c..f008e41d 100644 --- a/tutorials/sim/plot_SimulateAperiodic.py +++ b/tutorials/sim/plot_SimulateAperiodic.py @@ -56,7 +56,7 @@ # Set the exponent for brown noise, which is -2 exponent = -2 -# Simulate brown noise powerlaw activity +# Simulate powerlaw activity br_noise = sim_powerlaw(n_seconds, fs, exponent) ################################################################################################### @@ -82,12 +82,12 @@ # be applied to the simulated data before being returned. # # Here we will apply a high-pass filter. We can see that the resulting signal has much less -# low-frequency activity than the first one. +# low-frequency drift than the first one. # ################################################################################################### -# Simulate filtered brown noise with a 1 Hz highpass filter +# Simulate highpass-filtered brown noise with a 1Hz highpass filter f_hipass_brown = 1 brown_filt = sim_powerlaw(n_seconds, fs, exponent, f_range=(f_hipass_brown, None))