diff --git a/bci_essentials/paradigm/p300_paradigm.py b/bci_essentials/paradigm/p300_paradigm.py index c793c509..d69b803b 100644 --- a/bci_essentials/paradigm/p300_paradigm.py +++ b/bci_essentials/paradigm/p300_paradigm.py @@ -15,6 +15,7 @@ def __init__( epoch_start=0, epoch_end=0.6, buffer_time=0.01, + preprocessing_window=2, ): """ Parameters @@ -34,6 +35,9 @@ def __init__( buffer_time : float, *optional* Defines the time in seconds after an epoch for which we require EEG data to ensure that all EEG is present in that epoch. - Default is `0.01`. + preprocessing_window : float, *optional* + Defines the time in seconds before and after a marker to include additional EEG data in preprocessing. + - Default is `2`. """ super().__init__(filters) @@ -54,6 +58,8 @@ def __init__( self.buffer_time = buffer_time + self.preprocessing_window = preprocessing_window + def get_eeg_start_and_end_times(self, markers, timestamps): """ Get the start and end times of the EEG data based on the markers. @@ -129,8 +135,10 @@ def process_markers(self, markers, marker_timestamps, eeg, eeg_timestamps, fsamp # Subtract the marker timestamp from the EEG timestamps so that 0 becomes the marker onset marker_eeg_timestamps = eeg_timestamps - marker_timestamp - # Create the epoch time vector - epoch_time = np.arange(self.epoch_start, self.epoch_end, 1 / fsample) + # Create the epoch time vector using the preprocessing window + epoch_time = np.arange( + -self.preprocessing_window, self.preprocessing_window, 1 / fsample + ) epoch_X = np.zeros((1, n_channels, len(epoch_time))) @@ -150,23 +158,32 @@ def process_markers(self, markers, marker_timestamps, eeg, eeg_timestamps, fsamp epoch_X[0, :, :], fsample, self.lowcut, self.highcut ) + # Trim the preprocessed epoch + trimmed_epoch_X = epoch_X[ + :, + :, + ((epoch_time >= self.epoch_start) & (epoch_time <= self.epoch_end)), + ] + # For each flash index in the marker for flash_index in flash_indices: if flash_counts[flash_index] == 0: - object_epochs[flash_index] = epoch_X + object_epochs[flash_index] = trimmed_epoch_X flash_counts[flash_index] += 1 else: object_epochs[flash_index] = np.concatenate( - (object_epochs[flash_index], epoch_X), axis=0 + (object_epochs[flash_index], trimmed_epoch_X), axis=0 ) flash_counts[flash_index] += 1 # Average all epochs for each object - object_epochs_mean = [np.zeros((n_channels, len(epoch_time)))] * num_objects + object_epochs_mean = [ + np.zeros((n_channels, trimmed_epoch_X.shape[-1])) + ] * num_objects for i in range(num_objects): object_epochs_mean[i] = np.mean(object_epochs[i], axis=0) - X = np.zeros((num_objects, n_channels, len(epoch_time))) + X = np.zeros((num_objects, n_channels, trimmed_epoch_X.shape[-1])) for i in range(num_objects): X[i, :, :] = object_epochs_mean[i] # # object_epochs_mean = np.mean(object_epochs, axis=1) diff --git a/bci_essentials/paradigm/paradigm.py b/bci_essentials/paradigm/paradigm.py index a0c20086..67232420 100644 --- a/bci_essentials/paradigm/paradigm.py +++ b/bci_essentials/paradigm/paradigm.py @@ -36,23 +36,18 @@ def __init__(self, filters=[5, 30], channel_subset=None): def _preprocess(self, eeg, fsample, lowcut, highcut, order=5): """ Preprocess EEG data with the appropriate filter type: - - If the data is continuous (i.e., shape is [channels, samples]), a - bandpass filter is used. - - - If the data is epoched (i.e., shape is [epoch, channels, samples]), - the filter type depends on the signal length relative to the filter's settling time: - - If signal length > settling time: use bandpass filter - - If signal length ≤ settling time: use lowpass filter - - The settling time is calculated by: - 1. Compute the time constant (tc) of the highpass filter: - tc = 1 / (2 * π * lowcut) - 2. Compute the settling time with the formula: - settling_time = tc * 5 * order - - In a first-order system, a rule-of-thumb is that the signal settles in approximately 5 time constants (tc * 5) - - In higher-order filters (e.g., a 5th-order filter set as the default), the settling time incrases linearly with the order of the filter (order * tc * 5) - - This is a simplification, but it provides a good approximation for the settling time of the filter. - - For more details, see the reference below: - https://www.analogictips.com/an-overview-of-filters-and-their-parameters-part-4-time-and-phase-issues/ + - If signal length > settling time: use bandpass filter + - If signal length ≤ settling time: use lowpass filter + - The settling time is calculated by: + 1. Compute the time constant (tc) of the highpass filter: + tc = 1 / (2 * π * lowcut) + 2. Compute the settling time with the formula: + settling_time = tc * 5 * order + - In a first-order system, a rule-of-thumb is that the signal settles in approximately 5 time constants (tc * 5) + - In higher-order filters (e.g., a 5th-order filter set as the default), the settling time incrases linearly with the order of the filter (order * tc * 5) + - This is a simplification, but it provides a good approximation for the settling time of the filter. + - For more details, see the reference below: + https://www.analogictips.com/an-overview-of-filters-and-their-parameters-part-4-time-and-phase-issues/ Parameters ---------- @@ -74,31 +69,21 @@ def _preprocess(self, eeg, fsample, lowcut, highcut, order=5): Preprocessed EEG. Shape is the same as `eeg`. """ + logger.debug("Preprocessing EEG data") - n_dims = len(eeg.shape) - if n_dims == 2: - logger.debug("Preprocessing continuous EEG") + # Get the length of the signal + signal_length = eeg.shape[-1] / fsample + + # Highpass filter settling time + tc = 1 / (2 * np.pi * lowcut) + settling_time = order * tc * 5 + + if signal_length > settling_time: + logger.debug("Applied bandpass filter") preprocessed_eeg = bandpass(eeg, lowcut, highcut, order, fsample) - elif n_dims == 3: - logger.debug("Preprocessing epoched EEG") - - # Get the length of the signal - signal_length = eeg.shape[2] - - # Highpass filter settling time - tc = 1 / (2 * np.pi * lowcut) - settling_time = order * tc * 5 - - if signal_length > settling_time: - logger.debug("Applied bandpass filter to epoched EEG") - preprocessed_eeg = bandpass(eeg, lowcut, highcut, order, fsample) - else: - logger.debug("Applied lowpass filter to epoched EEG") - preprocessed_eeg = lowpass(eeg, highcut, order, fsample) else: - raise ValueError( - "Preprocessing failed. EEG must be 2D (continuous) or 3D (epoched)." - ) + logger.debug("Applied lowpass filter") + preprocessed_eeg = lowpass(eeg, highcut, order, fsample) return preprocessed_eeg