From fbee4be0f622d780ebc47e492e47fee8e2f2bbe8 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Wed, 30 Dec 2020 14:19:47 -0500 Subject: [PATCH 01/16] update code docs --- neurodsp/rhythm/swm.py | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/neurodsp/rhythm/swm.py b/neurodsp/rhythm/swm.py index 54bc1dd7..880880ee 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 @@ -21,18 +21,18 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, win_len : float Window length, in seconds. win_spacing : float - Minimum window spacing, in seconds. - max_iterations : int + Minimum spacing between start positions of successive windows, in seconds. + max_iterations : int, optional, default: 500 Maximum number of iterations of potential changes in window placement. - temperature : float - Temperature parameter. Controls probability of accepting a new window. + temperature : float, optional, default: 1 + Controls the probability of accepting a new window. window_starts_custom : 1d array, optional Custom pre-set locations of initial windows. Returns ------- avg_window : 1d array - The average waveform in 'sig' in the frequency 'f_range' triggered on 'trigger'. + The average pattern discovered in the input signal. window_starts : 1d array Indices at which each window begins for the final set of windows. costs : 1d array @@ -44,14 +44,15 @@ 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. Examples -------- @@ -87,8 +88,7 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, for iter_num in range(1, max_iterations): - # Pick a random window position to randomly replace with a - # new window to improve cross-window similarity + # Pick a random window position to replace, to improve cross-window similarity window_idx_replace = random_window_idx[iter_num] # Find a new allowed position for the window @@ -103,7 +103,7 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, # Calculate the acceptance probability p_accept = np.exp(-delta_cost / float(temperature)) - # Accept update to J with a certain probability + # Accept update, based on the calculated acceptance probability if np.random.rand() < p_accept: # Update costs & windows @@ -158,6 +158,8 @@ def _find_new_window_idx(window_starts, spacing_n_samps, n_samp, tries_limit=100 new_samp = np.random.randint(n_samp) dists = np.abs(window_starts - new_samp) + # TODO note: I think this should be less than. + # See comments in Issue #196 on this function if np.min(dists) > spacing_n_samps: break From 08d8a9b85256538d230cc98956e84a2eb4039f9b Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Wed, 30 Dec 2020 16:14:40 -0500 Subject: [PATCH 02/16] update & add more text to SWM tut --- .../rhythm/plot_SlidingWindowMatching.py | 103 +++++++++++++++--- 1 file changed, 85 insertions(+), 18 deletions(-) diff --git a/tutorials/rhythm/plot_SlidingWindowMatching.py b/tutorials/rhythm/plot_SlidingWindowMatching.py index 9c07c479..581a602c 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,10 +11,21 @@ # 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 `_ # ################################################################################################### @@ -39,6 +50,8 @@ # Load neural signal # ------------------ # +# First, we will load a segment of ECoG data, as an example time series. +# ################################################################################################### @@ -49,23 +62,58 @@ 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. # ################################################################################################### @@ -74,13 +122,21 @@ win_len = .055 win_spacing = .2 +################################################################################################### + # 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) +avg_window, window_starts, costs = sliding_window_matching(sig, fs, win_len, win_spacing) ################################################################################################### +# Examine the Results +# ~~~~~~~~~~~~~~~~~~~ # -# You can plot the resulting pattern with :func:`~.plot_swm_pattern`. +# 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. +# +# In order to visualize the resulting pattern, we can use +# :func:`~.plot_swm_pattern`. # ################################################################################################### @@ -90,10 +146,21 @@ ################################################################################################### # -# 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. # From d705ad93b73187e8aaa76f250d8d6c15ae28d2d9 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Thu, 31 Dec 2020 12:17:58 -0500 Subject: [PATCH 03/16] update word noise -> ap & misc --- tutorials/filt/plot_Filtering.py | 345 ++++++++++++++++++++++++ tutorials/sim/plot_SimulateAperiodic.py | 8 +- 2 files changed, 351 insertions(+), 2 deletions(-) create mode 100644 tutorials/filt/plot_Filtering.py diff --git a/tutorials/filt/plot_Filtering.py b/tutorials/filt/plot_Filtering.py new file mode 100644 index 00000000..a61f5d49 --- /dev/null +++ b/tutorials/filt/plot_Filtering.py @@ -0,0 +1,345 @@ +""" +Filtering +========= + +Apply digital filters to neural signals, including highpass, lowpass, bandpass & bandstop filters. + +This tutorial primarily covers the ``neurodsp.filt`` module. +""" + +################################################################################################### +# Filtering with NeuroDSP +# ----------------------- +# +# The :func:`~.filter_signal` function is the main function for filtering using NeuroDSP. +# +# Sections +# ~~~~~~~~ +# +# This tutorial contains the following sections: +# +# 1. Bandpass filter: extract a single oscillation from a signal +# 2. Highpass, lowpass, and bandstop filters: remove power in unwanted frequency ranges +# 3. Time-frequency resolution trade off: changing the filter length +# 4. Infinite-impulse-response (IIR) filter option +# 5. Beta bandpass filter on a neural signal +# + +################################################################################################### + +# sphinx_gallery_thumbnail_number = 12 + +# Import filter function +from neurodsp.filt import filter_signal + +# Import simulation code for creating test data +from neurodsp.sim import sim_combined +from neurodsp.utils import set_random_seed, create_times + +# Import utilities for loading and plotting data +from neurodsp.utils.download import load_ndsp_data +from neurodsp.plts.time_series import plot_time_series + +################################################################################################### + +# Set the random seed, for consistency simulating data +set_random_seed(0) + +################################################################################################### + +# General setting for simulations +fs = 1000 +n_seconds = 5 + +# Set the default aperiodic exponent +exp = -1 + +# Generate a times vector, for plotting +times = create_times(n_seconds, fs) + +################################################################################################### +# 1. Bandpass filter +# ------------------ +# +# Extract signal within a specific frequency range (e.g. theta, 4-8 Hz). +# + +################################################################################################### + +# Set the frequency in our simulated signal +freq = 6 + +# Set up simulation for a signal with aperiodic activity and an oscillation +components = {'sim_powerlaw' : {'exponent' : exp}, + 'sim_oscillation' : {'freq' : 6}} +variances = [0.1, 1] + +# Simulate our signal +sig = sim_combined(n_seconds, fs, components, variances) + +################################################################################################### + +# Define a frequency range to filter the data +f_range = (4, 8) + +# Bandpass filter the data, across the band of interest +sig_filt = filter_signal(sig, fs, 'bandpass', f_range) + +################################################################################################### + +# Plot filtered signal +plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered']) + +################################################################################################### +# +# Notice that the edges of the filtered signal are clipped (no red). +# +# Edge artifact removal is done by default in NeuroDSP filtering, because +# the signal samples at the edges only experienced part of the filter. +# +# To bypass this feature, set `remove_edges=False`, but at your own risk! +# + +################################################################################################### +# 2. Highpass, lowpass, and bandstop filters +# ------------------------------------------ +# +# 2a. Highpass filter +# ~~~~~~~~~~~~~~~~~~~ +# +# Remove low frequency drift from the data. +# + +################################################################################################### + +# Settings for the rhythmic components in the data +freq1 = 3 +freq2 = 0.5 + +# Set up simulation for a signal with aperiodic activity, an oscillation, and low frequency drift +components = {'sim_powerlaw' : {'exponent' : exp}, + 'sim_oscillation' : [{'freq' : freq1}, {'freq' : freq2}]} +variances = [0.1, 1, 1] + +# Generate a signal including low-frequency activity +sig = sim_combined(n_seconds, fs, components, variances) + +################################################################################################### + +# Filter the data with a highpass filter +f_range = (2, None) +sig_filt = filter_signal(sig, fs, 'highpass', f_range) + +################################################################################################### + +# Plot filtered signal +plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered']) + +################################################################################################### +# 2b. Lowpass filter +# ~~~~~~~~~~~~~~~~~~ +# +# Remove high frequency activity from the data. +# + +################################################################################################### + +# Filter the data +f_range = (None, 20) +sig_filt = filter_signal(sig, fs, 'lowpass', f_range) + +################################################################################################### + +# Plot filtered signal +plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered']) + +################################################################################################### +# 2c. Bandstop filter +# ~~~~~~~~~~~~~~~~~~~ +# +# Next let's try a bandstop filter, to remove 60Hz line noise from data. +# +# Notice that it is necessary to set a non-default filter length because +# a filter of length 3 cycles of a 58Hz oscillation would not attenuate +# the 60Hz oscillation much (try this yourself!). +# + +################################################################################################### + +# Generate a signal, with a low frequency oscillation and 60 Hz line noise +components = {'sim_oscillation' : [{'freq' : 6}, {'freq' : 60}]} +variances = [1, 0.2] +sig = sim_combined(n_seconds, fs, components, variances) + +################################################################################################### + +# Filter the data +f_range = (58, 62) +sig_filt = filter_signal(sig, fs, 'bandstop', f_range, n_seconds=0.5) + +################################################################################################### + +# Plot filtered signal +plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered']) + +################################################################################################### +# +# You might sometimes see a user warning that warns about the level of attenuation. +# +# You will see this warning whenever the filter you construct has a frequency response that +# does not hit a certain level of attenuation in the stopband. By default, the warning appears +# if the level of attenuation does not go below 20dB. +# +# You can check filter properties by plotting the frequency response when you apply a filter. +# + +################################################################################################### + +# Apply a short filter +# In this case, we won't achieve our desired attenuation +sig_filt = filter_signal(sig, fs, 'bandstop', f_range, + n_seconds=0.25, plot_properties=True) + +###################################################################################################v + +# This user warning disappears if we elongate the filter +sig_filt = filter_signal(sig, fs, 'bandstop', f_range, + n_seconds=1, plot_properties=True) + +################################################################################################### +# 3. Time-frequency resolution trade off +# -------------------------------------- +# +# With longer filter kernels, we get improved frequency resolution, but worse time resolution. +# +# Two bandpass filters (one long and one short) +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Notice that the short filter preserves the start of the oscillation better than the +# long filter (i.e. better time resolution). +# +# Notice that the long filter correctly removed the 1Hz oscillation, but the short filter +# did not (i.e. better frequency resolution). +# + +################################################################################################### + +# Reset simulation settings for this example +fs = 100 +n_seconds = 3 +times = create_times(n_seconds, fs) + +# Generate a signal with aperiodic activity, an oscillation and, and low frequency drift +components = {'sim_powerlaw' : {'exponent' : exp}, + 'sim_oscillation' : [{'freq' : 6}, {'freq' : 1}]} +variances = [0.1, 1, 1] +sig = sim_combined(n_seconds, fs, components, variances) + +# Set the first second to 0 +sig[:fs] = 0 + +################################################################################################### + +# Define the frequency band of interest +f_range = (4, 8) + +# Filter the data +sig_filt_short = filter_signal(sig, fs, 'bandpass', f_range, n_seconds=.1) +sig_filt_long = filter_signal(sig, fs, 'bandpass', f_range, n_seconds=1) + +################################################################################################### + +# Plot filtered signal +plot_time_series(times, [sig, sig_filt_short, sig_filt_long], + ['Raw', 'Short Filter', 'Long Filter']) + +################################################################################################### + +# Visualize the kernels and frequency responses +print('Short filter') +sig_filt_short = filter_signal(sig, fs, 'bandpass', f_range, n_seconds=.1, + plot_properties=True) +print('\n\nLong filter') +sig_filt_long = filter_signal(sig, fs, 'bandpass', f_range, n_seconds=1, + plot_properties=True) + +################################################################################################### +# 4. Infinite impulse response (IIR) filter option +# ------------------------------------------------ +# +# So far, the filters that we've been using are finite impulse response (FIR) filters. +# +# These filters are nice because we have good control over their properties, +# by manipulating the time-frequency resolution trade off through the filter length. +# +# However, sometimes we may not be as concerned with the precise filter properties, +# and so there is a faster option: IIR filters. +# +# We often use these filters for removing 60 Hz line noise. +# +# Here we apply a 3rd order Butterworth filter to remove 60Hz noise. +# +# Notice that some edge artifacts remain. +# + +################################################################################################### + +# Reset simulation settings +n_seconds = 1 +fs = 1000 +times = create_times(n_seconds, fs) + +# Generate a signal, with a low frequency oscillation and 60 Hz line noise +components = {'sim_oscillation' : [{'freq' : 6}, {'freq' : 60}]} +variances = [1, 0.2] +sig = sim_combined(n_seconds, fs, components, variances) + +################################################################################################### + +# Filter the data +f_range = (58, 62) +sig_filt = filter_signal(sig, fs, 'bandstop', f_range, + filter_type='iir', butterworth_order=3) + +################################################################################################### + +# Plot filtered signal +plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered']) + +################################################################################################### +# 5. Beta bandpass filter on neural signal +# ---------------------------------------- +# +# Next, we'll load an example time series of real data, and filter that. +# + +################################################################################################### + +# Download, if needed, and load example data file +sig = load_ndsp_data('sample_data_1.npy', folder='data') + +# Set sampling rate, and create a times vector for plotting +fs = 1000 +times = create_times(len(sig)/fs, fs) + +################################################################################################### + +# Filter the data +f_range = (13, 30) +sig_filt, kernel = filter_signal(sig, fs, 'bandpass', f_range, n_cycles=3, + plot_properties=True, return_filter=True) + +################################################################################################### + +# Plot filtered signal +plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered'], xlim=[2, 5]) + +################################################################################################### +# +# Notice that in the filtered time series, the resulting oscillation appears to be +# more sinusoidal than the original signal really is. +# +# If you are interested in this problem, and how to deal with it, you should check out +# `bycycle `_, +# which is a tool for time-domain analyses of waveform shape. +# diff --git a/tutorials/sim/plot_SimulateAperiodic.py b/tutorials/sim/plot_SimulateAperiodic.py index 0b73cb4c..c950e567 100644 --- a/tutorials/sim/plot_SimulateAperiodic.py +++ b/tutorials/sim/plot_SimulateAperiodic.py @@ -56,7 +56,11 @@ # Set the exponent for brown noise, which is -2 exponent = -2 +<<<<<<< HEAD # Simulate brown noise powerlaw activity +======= +# Simulate powerlaw activity +>>>>>>> update word noise -> ap & misc br_noise = sim_powerlaw(n_seconds, fs, exponent) ################################################################################################### @@ -82,12 +86,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)) From 0e15ad63f827f6cc523510ee1581f484c4e02065 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Thu, 31 Dec 2020 13:15:10 -0500 Subject: [PATCH 04/16] Update + add descriptions --- tutorials/rhythm/plot_LaggedCoherence.py | 209 ++++++++++++++++++----- 1 file changed, 170 insertions(+), 39 deletions(-) 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. +# From 9e16c7d154e40748372a1a42da7e3205cb3722f4 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Thu, 31 Dec 2020 13:20:32 -0500 Subject: [PATCH 05/16] update docs & code org --- neurodsp/rhythm/lc.py | 32 ++++++++++++++++++++++---------- 1 file changed, 22 insertions(+), 10 deletions(-) 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 From e8755883f326094955c1c56c5aaac2f9448a73c3 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Thu, 31 Dec 2020 13:48:45 -0500 Subject: [PATCH 06/16] Revert "update word noise -> ap & misc" This reverts commit 856a1fd2edde7ef442f5904c552c39c2baf5e9b4. --- tutorials/filt/plot_Filtering.py | 26 +++++++++++-------------- tutorials/sim/plot_SimulateAperiodic.py | 4 ---- 2 files changed, 11 insertions(+), 19 deletions(-) diff --git a/tutorials/filt/plot_Filtering.py b/tutorials/filt/plot_Filtering.py index a61f5d49..f04f0640 100644 --- a/tutorials/filt/plot_Filtering.py +++ b/tutorials/filt/plot_Filtering.py @@ -51,9 +51,6 @@ fs = 1000 n_seconds = 5 -# Set the default aperiodic exponent -exp = -1 - # Generate a times vector, for plotting times = create_times(n_seconds, fs) @@ -69,8 +66,8 @@ # Set the frequency in our simulated signal freq = 6 -# Set up simulation for a signal with aperiodic activity and an oscillation -components = {'sim_powerlaw' : {'exponent' : exp}, +# Set up simulation for a signal with an oscillation + noise +components = {'sim_powerlaw' : {'exponent' : 0}, 'sim_oscillation' : {'freq' : 6}} variances = [0.1, 1] @@ -116,8 +113,8 @@ freq1 = 3 freq2 = 0.5 -# Set up simulation for a signal with aperiodic activity, an oscillation, and low frequency drift -components = {'sim_powerlaw' : {'exponent' : exp}, +# Set up simulation for a signal with an oscillation + noise + low frequency activity +components = {'sim_powerlaw' : {'exponent' : 0}, 'sim_oscillation' : [{'freq' : freq1}, {'freq' : freq2}]} variances = [0.1, 1, 1] @@ -126,7 +123,7 @@ ################################################################################################### -# Filter the data with a highpass filter +# Filter the data f_range = (2, None) sig_filt = filter_signal(sig, fs, 'highpass', f_range) @@ -157,7 +154,7 @@ # 2c. Bandstop filter # ~~~~~~~~~~~~~~~~~~~ # -# Next let's try a bandstop filter, to remove 60Hz line noise from data. +# Next let's try a bandstop filter, for the example use case of removing 60Hz noise from data. # # Notice that it is necessary to set a non-default filter length because # a filter of length 3 cycles of a 58Hz oscillation would not attenuate @@ -210,7 +207,8 @@ # 3. Time-frequency resolution trade off # -------------------------------------- # -# With longer filter kernels, we get improved frequency resolution, but worse time resolution. +# With longer filter kernels, we get improved frequency resolution, +# but worse time resolution. # # Two bandpass filters (one long and one short) # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -229,8 +227,8 @@ n_seconds = 3 times = create_times(n_seconds, fs) -# Generate a signal with aperiodic activity, an oscillation and, and low frequency drift -components = {'sim_powerlaw' : {'exponent' : exp}, +# Generate a signal with an oscillation, noise, and a low frequency oscillation +components = {'sim_powerlaw' : {'exponent' : 0}, 'sim_oscillation' : [{'freq' : 6}, {'freq' : 1}]} variances = [0.1, 1, 1] sig = sim_combined(n_seconds, fs, components, variances) @@ -275,7 +273,7 @@ # However, sometimes we may not be as concerned with the precise filter properties, # and so there is a faster option: IIR filters. # -# We often use these filters for removing 60 Hz line noise. +# We often use these filters when removing 60 Hz line noise. # # Here we apply a 3rd order Butterworth filter to remove 60Hz noise. # @@ -310,8 +308,6 @@ # 5. Beta bandpass filter on neural signal # ---------------------------------------- # -# Next, we'll load an example time series of real data, and filter that. -# ################################################################################################### diff --git a/tutorials/sim/plot_SimulateAperiodic.py b/tutorials/sim/plot_SimulateAperiodic.py index c950e567..f008e41d 100644 --- a/tutorials/sim/plot_SimulateAperiodic.py +++ b/tutorials/sim/plot_SimulateAperiodic.py @@ -56,11 +56,7 @@ # Set the exponent for brown noise, which is -2 exponent = -2 -<<<<<<< HEAD -# Simulate brown noise powerlaw activity -======= # Simulate powerlaw activity ->>>>>>> update word noise -> ap & misc br_noise = sim_powerlaw(n_seconds, fs, exponent) ################################################################################################### From e9ce7072e36dfe0f40ea80a65b3013488d60af2e Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Thu, 7 Jan 2021 17:13:54 -0800 Subject: [PATCH 07/16] svm algorithm update --- neurodsp/rhythm/swm.py | 35 ++++++++++++++++++++++++----------- 1 file changed, 24 insertions(+), 11 deletions(-) diff --git a/neurodsp/rhythm/swm.py b/neurodsp/rhythm/swm.py index 880880ee..5646215e 100644 --- a/neurodsp/rhythm/swm.py +++ b/neurodsp/rhythm/swm.py @@ -93,11 +93,15 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, # Find a new allowed position for the window window_starts_temp = np.copy(window_starts) + + current_window_idx = window_starts_temp[window_idx_replace] + window_starts_temp[window_idx_replace] = _find_new_window_idx( - window_starts, spacing_n_samps, len(sig) - win_n_samps) + len(sig) - win_n_samps, spacing_n_samps, current_window_idx) # 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] # Calculate the acceptance probability @@ -130,7 +134,6 @@ def _compute_cost(sig, window_starts, win_n_samps): # Get all windows and z-score them n_windows = len(window_starts) windows = np.zeros((n_windows, win_n_samps)) - 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) @@ -149,19 +152,29 @@ def _compute_cost(sig, window_starts, win_n_samps): return cost -def _find_new_window_idx(window_starts, spacing_n_samps, n_samp, tries_limit=1000): +def _find_new_window_idx(max_start, spacing_n_samps, current_window_idx, tries_limit=1000): """Find a new sample for the starting window.""" - for n_try in range(tries_limit): + for _ in range(tries_limit): - # 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) + # Find adjacent window starts + prev_idx = current_window_idx - (2*spacing_n_samps) + next_idx = current_window_idx + (2*spacing_n_samps) + + # Create allowed random sample positions/indices + new_samps = np.arange(prev_idx + spacing_n_samps, next_idx - spacing_n_samps) + new_samps = new_samps[np.where((new_samps >= 0) & (new_samps <= max_start))[0]] + + if len(new_samps) <= 1: + # All new samples are too close to adjacent window starts + continue + else: + # Drop current window index + new_samps = np.delete(new_samps, np.where(new_samps == current_window_idx)[0][0]) - # TODO note: I think this should be less than. - # See comments in Issue #196 on this function - if np.min(dists) > spacing_n_samps: - break + # Randomly select allowed sample position + new_samp = np.random.choice(new_samps) + break else: raise RuntimeError('SWM algorithm has difficulty finding a new window. \ From c6174be9ed4655302ca50bc38f8c499f7e0b516d Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Tue, 26 Jan 2021 13:00:41 -0500 Subject: [PATCH 08/16] add more doc descriptions --- neurodsp/rhythm/swm.py | 45 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 40 insertions(+), 5 deletions(-) diff --git a/neurodsp/rhythm/swm.py b/neurodsp/rhythm/swm.py index 5646215e..fa75ec75 100644 --- a/neurodsp/rhythm/swm.py +++ b/neurodsp/rhythm/swm.py @@ -19,9 +19,9 @@ 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 spacing between start positions of successive windows, in seconds. + Minimum spacing between windows, in seconds. This is G in the original paper. max_iterations : int, optional, default: 500 Maximum number of iterations of potential changes in window placement. temperature : float, optional, default: 1 @@ -53,6 +53,10 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, - 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 version, as compared to the original implementation + in [2], which has more available options. + - This version may also be much smaller than the original, as this implementation does not + currently include the Markov Chain Monte Carlo sampling to calculate distances. Examples -------- @@ -101,7 +105,6 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, # 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] # Calculate the acceptance probability @@ -129,7 +132,22 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, 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, which is proportional to the difference between pairs of 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. + + Returns + ------- + cost: float + The calculated cost of the current set of windows. + """ # Get all windows and z-score them n_windows = len(window_starts) @@ -153,7 +171,24 @@ def _compute_cost(sig, window_starts, win_n_samps): def _find_new_window_idx(max_start, spacing_n_samps, current_window_idx, tries_limit=1000): - """Find a new sample for the starting window.""" + """Find a new sample for the starting window. + + Parameters + ---------- + max_start : int + The largest possible start index for the window. + spacing_n_samps : int + The minimum spacing required between windows, in samples. + current_window_idx : int + The current index of the window. + tries_limit : int, optional, default: False + The maximum number of iterations to attempt to find a new window. + + Returns + ------- + new_samp : int + The index for the new window. + """ for _ in range(tries_limit): From 983122db2f6a35e522599148529f7a690d4f64a9 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Tue, 26 Jan 2021 14:41:42 -0800 Subject: [PATCH 09/16] find new windows uses new positions --- neurodsp/rhythm/swm.py | 67 +++++++++++++++++++++++------------------- 1 file changed, 37 insertions(+), 30 deletions(-) diff --git a/neurodsp/rhythm/swm.py b/neurodsp/rhythm/swm.py index fa75ec75..06b85e87 100644 --- a/neurodsp/rhythm/swm.py +++ b/neurodsp/rhythm/swm.py @@ -83,8 +83,6 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, 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) # Calculate initial cost costs = np.zeros(max_iterations) @@ -92,16 +90,9 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, for iter_num in range(1, max_iterations): - # Pick a random window position to replace, to improve cross-window similarity - window_idx_replace = random_window_idx[iter_num] - - # Find a new allowed position for the window - window_starts_temp = np.copy(window_starts) - - current_window_idx = window_starts_temp[window_idx_replace] - - window_starts_temp[window_idx_replace] = _find_new_window_idx( - len(sig) - win_n_samps, spacing_n_samps, current_window_idx) + # Find a new, random window start + window_starts_temp = _find_new_window_idx(np.copy(window_starts), len(sig) - win_n_samps, + spacing_n_samps) # Calculate the cost & the change in the cost function cost_temp = _compute_cost(sig, window_starts_temp, win_n_samps) @@ -170,49 +161,65 @@ def _compute_cost(sig, window_starts, win_n_samps): return cost -def _find_new_window_idx(max_start, spacing_n_samps, current_window_idx, tries_limit=1000): +def _find_new_window_idx(windows, max_start, spacing_n_samps, tries_limit=1000): """Find a new sample for the starting window. Parameters ---------- + windows : 1d array + Indices at which each window begins. max_start : int The largest possible start index for the window. spacing_n_samps : int The minimum spacing required between windows, in samples. - current_window_idx : int - The current index of the window. tries_limit : int, optional, default: False The maximum number of iterations to attempt to find a new window. Returns ------- - new_samp : int - The index for the new window. + windows : 1d array + Indices, including the updated window, at which each window begins. """ + new_win = None + for _ in range(tries_limit): + # Randomly select current window + current_idx = np.random.choice(range(len(windows))) + current_win = windows[current_idx] + # Find adjacent window starts - prev_idx = current_window_idx - (2*spacing_n_samps) - next_idx = current_window_idx + (2*spacing_n_samps) + if current_idx != 0: + lower_bound = windows[current_idx - 1] + spacing_n_samps + else: + lower_bound = 0 - # Create allowed random sample positions/indices - new_samps = np.arange(prev_idx + spacing_n_samps, next_idx - spacing_n_samps) - new_samps = new_samps[np.where((new_samps >= 0) & (new_samps <= max_start))[0]] + if current_idx != len(windows) - 1: + upper_bound = windows[current_idx + 1] - spacing_n_samps + else: + upper_bound = max_start - if len(new_samps) <= 1: - # All new samples are too close to adjacent window starts + # Create allowed random positions/indices + new_wins = np.arange(lower_bound, upper_bound + 1) + new_wins = new_wins[np.where((new_wins >= 0) & (new_wins <= max_start))[0]] + + # Drop the current window index so it can't be randomly sampled + new_wins = np.delete(new_wins, np.where(new_wins == current_win)[0][0]) + + if len(new_wins) == 0: + # All new samples are too close to adjacent window starts, try another random position continue - else: - # Drop current window index - new_samps = np.delete(new_samps, np.where(new_samps == current_window_idx)[0][0]) # Randomly select allowed sample position - new_samp = np.random.choice(new_samps) + new_win= np.random.choice(new_wins) + break - else: + if new_win is None: raise RuntimeError('SWM algorithm has difficulty finding a new window. \ Try increasing the spacing parameter.') - return new_samp + windows[current_idx] = new_win + + return windows From f8b899dd61af89a1dee20d551f693df65cbd88a0 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Mon, 22 Feb 2021 12:33:40 -0800 Subject: [PATCH 10/16] typo: smaller -> slower --- neurodsp/rhythm/swm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neurodsp/rhythm/swm.py b/neurodsp/rhythm/swm.py index 06b85e87..b6f51224 100644 --- a/neurodsp/rhythm/swm.py +++ b/neurodsp/rhythm/swm.py @@ -55,7 +55,7 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, so that the algorithm does not converge on a low frequency motif. - This implementation is a minimal version, as compared to the original implementation in [2], which has more available options. - - This version may also be much smaller than the original, as this implementation does not + - This version may also be much slower than the original, as this implementation does not currently include the Markov Chain Monte Carlo sampling to calculate distances. Examples From 0f3f6f282f76770bd1dc2593ba3e15f5303ea7f3 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Wed, 2 Jun 2021 20:19:52 -0700 Subject: [PATCH 11/16] swm updated --- neurodsp/rhythm/swm.py | 225 +++++++++++++----------------- neurodsp/tests/rhythm/test_swm.py | 8 +- 2 files changed, 101 insertions(+), 132 deletions(-) diff --git a/neurodsp/rhythm/swm.py b/neurodsp/rhythm/swm.py index b6f51224..dcc49a05 100644 --- a/neurodsp/rhythm/swm.py +++ b/neurodsp/rhythm/swm.py @@ -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 @@ -22,21 +22,20 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, Window length, in seconds. This is L in the original paper. win_spacing : float Minimum spacing between windows, in seconds. This is G in the original paper. - max_iterations : int, optional, default: 500 + max_iterations : int, optional, default: 100 Maximum number of iterations of potential changes in window placement. - temperature : float, optional, default: 1 - Controls the 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 pattern discovered in the input signal. + 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 ---------- @@ -53,173 +52,143 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=500, - 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 version, as compared to the original implementation - in [2], which has more available options. - - This version may also be much slower than the original, as this implementation does not - currently include the Markov Chain Monte Carlo sampling to calculate distances. + - 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, + ... max_iterations=100, 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) + windows = np.array([sig[start:start+win_len] for start in window_starts]) - # Calculate initial cost - costs = np.zeros(max_iterations) - costs[0] = _compute_cost(sig, window_starts, win_n_samps) + # New window bounds + lower_bounds, upper_bounds = _compute_bounds(window_starts, win_spacing, 0, len(sig) - win_len) - for iter_num in range(1, max_iterations): + # Remove low variance windows to speed up runtime + if var_thresh != None: - # Find a new, random window start - window_starts_temp = _find_new_window_idx(np.copy(window_starts), len(sig) - win_n_samps, - spacing_n_samps) + thresh = np.array([np.var(sig[loc:loc+win_len]) > var_thresh for loc in window_starts]) - # 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] + windows = windows[thresh] + window_starts = window_starts[thresh] + lower_bounds = lower_bounds[thresh] + upper_bounds = upper_bounds[thresh] - # Calculate the acceptance probability - p_accept = np.exp(-delta_cost / float(temperature)) + # Modified SWM procedure + window_idxs = np.arange(len(windows)).astype(int) - # Accept update, based on the calculated acceptance probability - if np.random.rand() < p_accept: + corrs, variance = _compute_cost(sig, window_starts, win_len) + mae = np.mean(np.abs(windows - windows.mean(axis=0))) - # Update costs & windows - costs[iter_num] = cost_temp - window_starts = window_starts_temp + for _ in range(max_iterations): - else: + # Randomly shuffle order of windows + np.random.shuffle(window_idxs) - # Update costs - costs[iter_num] = costs[iter_num - 1] + for win_idx in window_idxs: - # 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) + # 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)) - return avg_window, window_starts, costs + # Accept new window if correlation increases and variance decreases + _corrs, _variance = _compute_cost(sig, _window_starts, win_len) + if _corrs[win_idx].sum() > corrs[win_idx].sum() and _variance < variance: -def _compute_cost(sig, window_starts, win_n_samps): - """Compute the cost, which is proportional to the difference between pairs of windows. + 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) - 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. + # Phase optimization + _window_starts = window_starts.copy() - Returns - ------- - cost: float - The calculated cost of the current set of windows. - """ + for shift in np.arange(-int(win_len/2), int(win_len/2)): - # Get all windows and z-score them - n_windows = len(window_starts) - windows = np.zeros((n_windows, win_n_samps)) - 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) + _starts = _window_starts + shift - # 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) + # Skip windows shifts that are out-of-bounds + if (_starts[0] < 0) or (_starts[-1] > len(sig) - win_len): + continue - # Calculate cost function, which is the average difference, roughly - cost = np.sum(dists) / float(2 * (n_windows - 1)) + _windows = np.array([sig[start:start+win_len] for start in _starts]) - return cost + _mae = np.mean(np.abs(_windows - _windows.mean(axis=0))) + if _mae < mae: + window_starts = _starts.copy() + windows = _windows.copy() + mae = _mae -def _find_new_window_idx(windows, max_start, spacing_n_samps, tries_limit=1000): - """Find a new sample for the starting window. + 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, start=None, end=None): + """Compute the cost, as corrleation coefficients and variance across windows. Parameters ---------- - windows : 1d array - Indices at which each window begins. - max_start : int - The largest possible start index for the window. - spacing_n_samps : int - The minimum spacing required between windows, in samples. - tries_limit : int, optional, default: False - The maximum number of iterations to attempt to find a new window. + 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. Returns ------- - windows : 1d array - Indices, including the updated window, at which each window begins. + corrs: 2d array + Window correlation matrix. + variance: float + Sum of the variance across windows. """ - new_win = None - - for _ in range(tries_limit): - - # Randomly select current window - current_idx = np.random.choice(range(len(windows))) - current_win = windows[current_idx] - - # Find adjacent window starts - if current_idx != 0: - lower_bound = windows[current_idx - 1] + spacing_n_samps - else: - lower_bound = 0 - - if current_idx != len(windows) - 1: - upper_bound = windows[current_idx + 1] - spacing_n_samps - else: - upper_bound = max_start + windows = np.array([sig[start:start+win_n_samps] for start in window_starts]) - # Create allowed random positions/indices - new_wins = np.arange(lower_bound, upper_bound + 1) - new_wins = new_wins[np.where((new_wins >= 0) & (new_wins <= max_start))[0]] + corrs = np.corrcoef(windows) - # Drop the current window index so it can't be randomly sampled - new_wins = np.delete(new_wins, np.where(new_wins == current_win)[0][0]) + variance = windows.var(axis=1).sum() - if len(new_wins) == 0: - # All new samples are too close to adjacent window starts, try another random position - continue + return corrs, variance - # Randomly select allowed sample position - new_win= np.random.choice(new_wins) - break +def _compute_bounds(window_starts, win_spacing, start=None, end=None): - if new_win is None: - raise RuntimeError('SWM algorithm has difficulty finding a new window. \ - Try increasing the spacing parameter.') + lower_bounds = window_starts[:-1] + win_spacing + lower_bounds = np.insert(lower_bounds, 0, start) - windows[current_idx] = new_win + upper_bounds = window_starts[1:] - win_spacing + upper_bounds = np.insert(upper_bounds, len(upper_bounds), end) - return windows + 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) From 1ece0960c8c00f4024844d09608470a8ee14e23d Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Thu, 3 Jun 2021 10:03:14 -0400 Subject: [PATCH 12/16] drop old filtering tut --- tutorials/filt/plot_Filtering.py | 341 ------------------------------- 1 file changed, 341 deletions(-) delete mode 100644 tutorials/filt/plot_Filtering.py diff --git a/tutorials/filt/plot_Filtering.py b/tutorials/filt/plot_Filtering.py deleted file mode 100644 index f04f0640..00000000 --- a/tutorials/filt/plot_Filtering.py +++ /dev/null @@ -1,341 +0,0 @@ -""" -Filtering -========= - -Apply digital filters to neural signals, including highpass, lowpass, bandpass & bandstop filters. - -This tutorial primarily covers the ``neurodsp.filt`` module. -""" - -################################################################################################### -# Filtering with NeuroDSP -# ----------------------- -# -# The :func:`~.filter_signal` function is the main function for filtering using NeuroDSP. -# -# Sections -# ~~~~~~~~ -# -# This tutorial contains the following sections: -# -# 1. Bandpass filter: extract a single oscillation from a signal -# 2. Highpass, lowpass, and bandstop filters: remove power in unwanted frequency ranges -# 3. Time-frequency resolution trade off: changing the filter length -# 4. Infinite-impulse-response (IIR) filter option -# 5. Beta bandpass filter on a neural signal -# - -################################################################################################### - -# sphinx_gallery_thumbnail_number = 12 - -# Import filter function -from neurodsp.filt import filter_signal - -# Import simulation code for creating test data -from neurodsp.sim import sim_combined -from neurodsp.utils import set_random_seed, create_times - -# Import utilities for loading and plotting data -from neurodsp.utils.download import load_ndsp_data -from neurodsp.plts.time_series import plot_time_series - -################################################################################################### - -# Set the random seed, for consistency simulating data -set_random_seed(0) - -################################################################################################### - -# General setting for simulations -fs = 1000 -n_seconds = 5 - -# Generate a times vector, for plotting -times = create_times(n_seconds, fs) - -################################################################################################### -# 1. Bandpass filter -# ------------------ -# -# Extract signal within a specific frequency range (e.g. theta, 4-8 Hz). -# - -################################################################################################### - -# Set the frequency in our simulated signal -freq = 6 - -# Set up simulation for a signal with an oscillation + noise -components = {'sim_powerlaw' : {'exponent' : 0}, - 'sim_oscillation' : {'freq' : 6}} -variances = [0.1, 1] - -# Simulate our signal -sig = sim_combined(n_seconds, fs, components, variances) - -################################################################################################### - -# Define a frequency range to filter the data -f_range = (4, 8) - -# Bandpass filter the data, across the band of interest -sig_filt = filter_signal(sig, fs, 'bandpass', f_range) - -################################################################################################### - -# Plot filtered signal -plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered']) - -################################################################################################### -# -# Notice that the edges of the filtered signal are clipped (no red). -# -# Edge artifact removal is done by default in NeuroDSP filtering, because -# the signal samples at the edges only experienced part of the filter. -# -# To bypass this feature, set `remove_edges=False`, but at your own risk! -# - -################################################################################################### -# 2. Highpass, lowpass, and bandstop filters -# ------------------------------------------ -# -# 2a. Highpass filter -# ~~~~~~~~~~~~~~~~~~~ -# -# Remove low frequency drift from the data. -# - -################################################################################################### - -# Settings for the rhythmic components in the data -freq1 = 3 -freq2 = 0.5 - -# Set up simulation for a signal with an oscillation + noise + low frequency activity -components = {'sim_powerlaw' : {'exponent' : 0}, - 'sim_oscillation' : [{'freq' : freq1}, {'freq' : freq2}]} -variances = [0.1, 1, 1] - -# Generate a signal including low-frequency activity -sig = sim_combined(n_seconds, fs, components, variances) - -################################################################################################### - -# Filter the data -f_range = (2, None) -sig_filt = filter_signal(sig, fs, 'highpass', f_range) - -################################################################################################### - -# Plot filtered signal -plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered']) - -################################################################################################### -# 2b. Lowpass filter -# ~~~~~~~~~~~~~~~~~~ -# -# Remove high frequency activity from the data. -# - -################################################################################################### - -# Filter the data -f_range = (None, 20) -sig_filt = filter_signal(sig, fs, 'lowpass', f_range) - -################################################################################################### - -# Plot filtered signal -plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered']) - -################################################################################################### -# 2c. Bandstop filter -# ~~~~~~~~~~~~~~~~~~~ -# -# Next let's try a bandstop filter, for the example use case of removing 60Hz noise from data. -# -# Notice that it is necessary to set a non-default filter length because -# a filter of length 3 cycles of a 58Hz oscillation would not attenuate -# the 60Hz oscillation much (try this yourself!). -# - -################################################################################################### - -# Generate a signal, with a low frequency oscillation and 60 Hz line noise -components = {'sim_oscillation' : [{'freq' : 6}, {'freq' : 60}]} -variances = [1, 0.2] -sig = sim_combined(n_seconds, fs, components, variances) - -################################################################################################### - -# Filter the data -f_range = (58, 62) -sig_filt = filter_signal(sig, fs, 'bandstop', f_range, n_seconds=0.5) - -################################################################################################### - -# Plot filtered signal -plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered']) - -################################################################################################### -# -# You might sometimes see a user warning that warns about the level of attenuation. -# -# You will see this warning whenever the filter you construct has a frequency response that -# does not hit a certain level of attenuation in the stopband. By default, the warning appears -# if the level of attenuation does not go below 20dB. -# -# You can check filter properties by plotting the frequency response when you apply a filter. -# - -################################################################################################### - -# Apply a short filter -# In this case, we won't achieve our desired attenuation -sig_filt = filter_signal(sig, fs, 'bandstop', f_range, - n_seconds=0.25, plot_properties=True) - -###################################################################################################v - -# This user warning disappears if we elongate the filter -sig_filt = filter_signal(sig, fs, 'bandstop', f_range, - n_seconds=1, plot_properties=True) - -################################################################################################### -# 3. Time-frequency resolution trade off -# -------------------------------------- -# -# With longer filter kernels, we get improved frequency resolution, -# but worse time resolution. -# -# Two bandpass filters (one long and one short) -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# Notice that the short filter preserves the start of the oscillation better than the -# long filter (i.e. better time resolution). -# -# Notice that the long filter correctly removed the 1Hz oscillation, but the short filter -# did not (i.e. better frequency resolution). -# - -################################################################################################### - -# Reset simulation settings for this example -fs = 100 -n_seconds = 3 -times = create_times(n_seconds, fs) - -# Generate a signal with an oscillation, noise, and a low frequency oscillation -components = {'sim_powerlaw' : {'exponent' : 0}, - 'sim_oscillation' : [{'freq' : 6}, {'freq' : 1}]} -variances = [0.1, 1, 1] -sig = sim_combined(n_seconds, fs, components, variances) - -# Set the first second to 0 -sig[:fs] = 0 - -################################################################################################### - -# Define the frequency band of interest -f_range = (4, 8) - -# Filter the data -sig_filt_short = filter_signal(sig, fs, 'bandpass', f_range, n_seconds=.1) -sig_filt_long = filter_signal(sig, fs, 'bandpass', f_range, n_seconds=1) - -################################################################################################### - -# Plot filtered signal -plot_time_series(times, [sig, sig_filt_short, sig_filt_long], - ['Raw', 'Short Filter', 'Long Filter']) - -################################################################################################### - -# Visualize the kernels and frequency responses -print('Short filter') -sig_filt_short = filter_signal(sig, fs, 'bandpass', f_range, n_seconds=.1, - plot_properties=True) -print('\n\nLong filter') -sig_filt_long = filter_signal(sig, fs, 'bandpass', f_range, n_seconds=1, - plot_properties=True) - -################################################################################################### -# 4. Infinite impulse response (IIR) filter option -# ------------------------------------------------ -# -# So far, the filters that we've been using are finite impulse response (FIR) filters. -# -# These filters are nice because we have good control over their properties, -# by manipulating the time-frequency resolution trade off through the filter length. -# -# However, sometimes we may not be as concerned with the precise filter properties, -# and so there is a faster option: IIR filters. -# -# We often use these filters when removing 60 Hz line noise. -# -# Here we apply a 3rd order Butterworth filter to remove 60Hz noise. -# -# Notice that some edge artifacts remain. -# - -################################################################################################### - -# Reset simulation settings -n_seconds = 1 -fs = 1000 -times = create_times(n_seconds, fs) - -# Generate a signal, with a low frequency oscillation and 60 Hz line noise -components = {'sim_oscillation' : [{'freq' : 6}, {'freq' : 60}]} -variances = [1, 0.2] -sig = sim_combined(n_seconds, fs, components, variances) - -################################################################################################### - -# Filter the data -f_range = (58, 62) -sig_filt = filter_signal(sig, fs, 'bandstop', f_range, - filter_type='iir', butterworth_order=3) - -################################################################################################### - -# Plot filtered signal -plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered']) - -################################################################################################### -# 5. Beta bandpass filter on neural signal -# ---------------------------------------- -# - -################################################################################################### - -# Download, if needed, and load example data file -sig = load_ndsp_data('sample_data_1.npy', folder='data') - -# Set sampling rate, and create a times vector for plotting -fs = 1000 -times = create_times(len(sig)/fs, fs) - -################################################################################################### - -# Filter the data -f_range = (13, 30) -sig_filt, kernel = filter_signal(sig, fs, 'bandpass', f_range, n_cycles=3, - plot_properties=True, return_filter=True) - -################################################################################################### - -# Plot filtered signal -plot_time_series(times, [sig, sig_filt], ['Raw', 'Filtered'], xlim=[2, 5]) - -################################################################################################### -# -# Notice that in the filtered time series, the resulting oscillation appears to be -# more sinusoidal than the original signal really is. -# -# If you are interested in this problem, and how to deal with it, you should check out -# `bycycle `_, -# which is a tool for time-domain analyses of waveform shape. -# From 8b0985bdfae64b2360c4eefa29819d2245220b47 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Thu, 3 Jun 2021 10:13:39 -0400 Subject: [PATCH 13/16] fix doctest for plotting SWM --- neurodsp/plts/rhythm.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) 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) """ From ac8d8f4f017b7d77cd48ef76f14ccd26db8c5264 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Thu, 3 Jun 2021 10:13:51 -0400 Subject: [PATCH 14/16] fix tutorial for SWM update --- tutorials/rhythm/plot_SlidingWindowMatching.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/tutorials/rhythm/plot_SlidingWindowMatching.py b/tutorials/rhythm/plot_SlidingWindowMatching.py index 581a602c..8214fc64 100644 --- a/tutorials/rhythm/plot_SlidingWindowMatching.py +++ b/tutorials/rhythm/plot_SlidingWindowMatching.py @@ -32,6 +32,8 @@ # sphinx_gallery_thumbnail_number = 2 +import numpy as np + # Import the sliding window matching function from neurodsp.rhythm import sliding_window_matching @@ -125,7 +127,7 @@ ################################################################################################### # Apply the sliding window matching algorithm to the time series -avg_window, window_starts, costs = sliding_window_matching(sig, fs, win_len, win_spacing) +windows, window_starts = sliding_window_matching(sig, fs, win_len, win_spacing) ################################################################################################### # Examine the Results @@ -141,6 +143,9 @@ ################################################################################################### +# Compute the average window +avg_window = np.mean(windows, 0) + # Plot the discovered pattern plot_swm_pattern(avg_window) From 10143f6f1fdb67b1d564839165cc9a3e63e2adef Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Thu, 3 Jun 2021 10:35:16 -0400 Subject: [PATCH 15/16] lints --- neurodsp/rhythm/swm.py | 52 +++++++++++++------ .../rhythm/plot_SlidingWindowMatching.py | 1 + 2 files changed, 36 insertions(+), 17 deletions(-) diff --git a/neurodsp/rhythm/swm.py b/neurodsp/rhythm/swm.py index dcc49a05..14cb9ece 100644 --- a/neurodsp/rhythm/swm.py +++ b/neurodsp/rhythm/swm.py @@ -9,7 +9,7 @@ @multidim() def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=100, - window_starts_custom=None, var_thresh=None): + window_starts_custom=None, var_thresh=None): """Find recurring patterns in a time series using the sliding window matching algorithm. Parameters @@ -69,11 +69,11 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=100, Search for reoccuring patterns using sliding window matching in a simulated beta signal: >>> from neurodsp.sim import sim_combined - >>> components={'sim_bursty_oscillation': {'freq': 20, 'phase':'min'}, - ... 'sim_powerlaw': {'f_range': (2, None)}} + >>> 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, - ... max_iterations=100, var_thresh=.5) + >>> 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 @@ -86,15 +86,15 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=100, else: window_starts = window_starts_custom - windows = np.array([sig[start:start+win_len] for start in window_starts]) + windows = np.array([sig[start:start + win_len] for start in window_starts]) - # New window bounds + # 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 != None: + if var_thresh is not None: - thresh = np.array([np.var(sig[loc:loc+win_len]) > var_thresh for loc in window_starts]) + 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] @@ -117,7 +117,7 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=100, # 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)) + upper_bounds[win_idx] + 1)) # Accept new window if correlation increases and variance decreases _corrs, _variance = _compute_cost(sig, _window_starts, win_len) @@ -127,7 +127,8 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=100, 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) + lower_bounds, upper_bounds = _compute_bounds(\ + window_starts, win_spacing, 0, len(sig) - win_len) # Phase optimization _window_starts = window_starts.copy() @@ -140,7 +141,7 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=100, if (_starts[0] < 0) or (_starts[-1] > len(sig) - win_len): continue - _windows = np.array([sig[start:start+win_len] for start in _starts]) + _windows = np.array([sig[start:start + win_len] for start in _starts]) _mae = np.mean(np.abs(_windows - _windows.mean(axis=0))) @@ -149,12 +150,13 @@ def sliding_window_matching(sig, fs, win_len, win_spacing, max_iterations=100, windows = _windows.copy() mae = _mae - lower_bounds, upper_bounds = _compute_bounds(window_starts, win_spacing, 0, len(sig) - win_len) + 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, start=None, end=None): +def _compute_cost(sig, window_starts, win_n_samps): """Compute the cost, as corrleation coefficients and variance across windows. Parameters @@ -174,7 +176,7 @@ def _compute_cost(sig, window_starts, win_n_samps, start=None, end=None): Sum of the variance across windows. """ - windows = np.array([sig[start:start+win_n_samps] for start in window_starts]) + windows = np.array([sig[start:start + win_n_samps] for start in window_starts]) corrs = np.corrcoef(windows) @@ -183,12 +185,28 @@ def _compute_cost(sig, window_starts, win_n_samps, start=None, end=None): return corrs, variance -def _compute_bounds(window_starts, win_spacing, start=None, end=None): +def _compute_bounds(window_starts, win_spacing, start, end): + """Compute bounds on a new window. + + 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. + + Returns + ------- + lower_bounds, upper_bounds : 1d array + Computed upper and lower bounds for the window position. + """ lower_bounds = window_starts[:-1] + win_spacing lower_bounds = np.insert(lower_bounds, 0, start) - upper_bounds = window_starts[1:] - win_spacing + upper_bounds = window_starts[1:] - win_spacing upper_bounds = np.insert(upper_bounds, len(upper_bounds), end) return lower_bounds, upper_bounds diff --git a/tutorials/rhythm/plot_SlidingWindowMatching.py b/tutorials/rhythm/plot_SlidingWindowMatching.py index 8214fc64..2786a5ed 100644 --- a/tutorials/rhythm/plot_SlidingWindowMatching.py +++ b/tutorials/rhythm/plot_SlidingWindowMatching.py @@ -84,6 +84,7 @@ # # Sliding window matching can be applied with the # :func:`~.sliding_window_matching` function. +# ################################################################################################### # Data Preprocessing From 5a9b8b6b8f15818ded0481c2575308d6c3d5980a Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Thu, 3 Jun 2021 11:47:31 -0700 Subject: [PATCH 16/16] update swm tutorial --- tutorials/rhythm/plot_SlidingWindowMatching.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tutorials/rhythm/plot_SlidingWindowMatching.py b/tutorials/rhythm/plot_SlidingWindowMatching.py index 2786a5ed..8343385a 100644 --- a/tutorials/rhythm/plot_SlidingWindowMatching.py +++ b/tutorials/rhythm/plot_SlidingWindowMatching.py @@ -42,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 ################################################################################################### @@ -59,6 +60,7 @@ # 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 @@ -123,12 +125,12 @@ # 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 -windows, window_starts = sliding_window_matching(sig, fs, win_len, win_spacing) +windows, window_starts = sliding_window_matching(sig, fs, win_len, win_spacing, var_thresh=.5) ################################################################################################### # Examine the Results