import numpy as np from .peak_finder import peak_finder from .. import pick_types, pick_channels from ..utils import logger, verbose from ..filter import band_pass_filter from ..epochs import Epochs @verbose def find_eog_events(raw, event_id=998, l_freq=1, h_freq=10, filter_length='10s', ch_name=None, tstart=0, verbose=None): """Locate EOG artifacts Parameters ---------- raw : instance of Raw The raw data. event_id : int The index to assign to found events. low_pass : float Low pass frequency. high_pass : float High pass frequency. filter_length : str | int | None Number of taps to use for filtering. ch_name: str | None If not None, use specified channel(s) for EOG tstart : float Start detection after tstart seconds. verbose : bool, str, int, or None If not None, override default verbose level (see mne.verbose). Returns ------- eog_events : array Events. """ # Getting EOG Channel eog_inds = _get_eog_channel_index(ch_name, raw) logger.info('EOG channel index for this subject is: %s' % eog_inds) eog, _ = raw[eog_inds, :] eog_events = _find_eog_events(eog, event_id=event_id, l_freq=l_freq, h_freq=h_freq, sampling_rate=raw.info['sfreq'], first_samp=raw.first_samp, filter_length=filter_length, tstart=tstart) return eog_events def _find_eog_events(eog, event_id, l_freq, h_freq, sampling_rate, first_samp, filter_length='10s', tstart=0.): """Helper function""" logger.info('Filtering the data to remove DC offset to help ' 'distinguish blinks from saccades') # filtering to remove dc offset so that we know which is blink and saccades fmax = np.minimum(45, sampling_rate / 2.0 - 0.75) # protect Nyquist filteog = np.array([band_pass_filter(x, sampling_rate, 2, fmax, filter_length=filter_length) for x in eog]) temp = np.sqrt(np.sum(filteog ** 2, axis=1)) indexmax = np.argmax(temp) # easier to detect peaks with filtering. filteog = band_pass_filter(eog[indexmax], sampling_rate, l_freq, h_freq, filter_length=filter_length) # detecting eog blinks and generating event file logger.info('Now detecting blinks and generating corresponding events') temp = filteog - np.mean(filteog) n_samples_start = int(sampling_rate * tstart) if np.abs(np.max(temp)) > np.abs(np.min(temp)): eog_events, _ = peak_finder(filteog[n_samples_start:], extrema=1) else: eog_events, _ = peak_finder(filteog[n_samples_start:], extrema=-1) eog_events += n_samples_start n_events = len(eog_events) logger.info("Number of EOG events detected : %d" % n_events) eog_events = np.c_[eog_events + first_samp, np.zeros(n_events), event_id * np.ones(n_events)] return eog_events def _get_eog_channel_index(ch_name, inst): if isinstance(ch_name, str): # Check if multiple EOG Channels if ',' in ch_name: ch_name = ch_name.split(',') else: ch_name = [ch_name] eog_inds = pick_channels(inst.ch_names, include=ch_name) if len(eog_inds) == 0: raise ValueError('%s not in channel list' % ch_name) else: logger.info('Using channel %s as EOG channel%s' % ( " and ".join(ch_name), '' if len(eog_inds) < 2 else 's')) elif ch_name is None: eog_inds = pick_types(inst.info, meg=False, eeg=False, stim=False, eog=True, ecg=False, emg=False, ref_meg=False, exclude='bads') if len(eog_inds) == 0: logger.info('No EOG channels found') logger.info('Trying with EEG 061 and EEG 062') eog_inds = pick_channels(inst.ch_names, include=['EEG 061', 'EEG 062']) if len(eog_inds) != 2: raise RuntimeError('EEG 61 or EEG 62 channel not found !!') else: raise ValueError('Could not find EOG channel.') return eog_inds @verbose def create_eog_epochs(raw, ch_name=None, event_id=998, picks=None, tmin=-0.5, tmax=0.5, l_freq=1, h_freq=10, reject=None, flat=None, baseline=None, verbose=None): """Conveniently generate epochs around EOG artifact events Parameters ---------- raw : instance of Raw The raw data ch_name : str The name of the channel to use for ECG peak detection. The argument is mandatory if the dataset contains no ECG channels. event_id : int The index to assign to found events picks : array-like of int | None (default) Indices of channels to include (if None, all channels are used). tmin : float Start time before event. tmax : float End time after event. l_freq : float Low pass frequency. h_freq : float High pass frequency. reject : dict | None Rejection parameters based on peak to peak amplitude. Valid keys are 'grad' | 'mag' | 'eeg' | 'eog' | 'ecg'. If reject is None then no rejection is done. You should use such parameters to reject big measurement artifacts and not ECG for example flat : dict | None Rejection parameters based on flatness of signal Valid keys are 'grad' | 'mag' | 'eeg' | 'eog' | 'ecg' If flat is None then no rejection is done. verbose : bool, str, int, or None If not None, override default verbose level (see mne.verbose). baseline : tuple or list of length 2, or None The time interval to apply rescaling / baseline correction. If None do not apply it. If baseline is (a, b) the interval is between "a (s)" and "b (s)". If a is None the beginning of the data is used and if b is None then b is set to the end of the interval. If baseline is equal ot (None, None) all the time interval is used. If None, no correction is applied. Returns ------- ecg_epochs : instance of Epochs Data epoched around ECG r-peaks. """ events = find_eog_events(raw, ch_name=ch_name, event_id=event_id, l_freq=l_freq, h_freq=h_freq) # create epochs around EOG events eog_epochs = Epochs(raw, events=events, event_id=event_id, tmin=tmin, tmax=tmax, proj=False, reject=reject, flat=flat, picks=picks, baseline=baseline, preload=True) return eog_epochs