Aluode/PerceptionLabPortable
0
1# Authors: The MNE-Python contributors.2# License: BSD-3-Clause3# Copyright the MNE-Python contributors.4 5import numpy as np6 7from .._fiff.meas_info import create_info8from .._fiff.pick import _picks_to_idx, pick_channels, pick_types9from ..annotations import _annotations_starts_stops10from ..epochs import BaseEpochs, Epochs11from ..evoked import Evoked12from ..filter import filter_data13from ..io import BaseRaw, RawArray14from ..utils import int_like, logger, sum_squared, verbose, warn15 16 17@verbose18def qrs_detector(19 sfreq,20 ecg,21 thresh_value=0.6,22 levels=2.5,23 n_thresh=3,24 l_freq=5,25 h_freq=35,26 tstart=0,27 filter_length="10s",28 verbose=None,29):30 """Detect QRS component in ECG channels.31 32 QRS is the main wave on the heart beat.33 34 Parameters35 ----------36 sfreq : float37 Sampling rate38 ecg : array39 ECG signal40 thresh_value : float | str41 qrs detection threshold. Can also be "auto" for automatic42 selection of threshold.43 levels : float44 number of std from mean to include for detection45 n_thresh : int46 max number of crossings47 l_freq : float48 Low pass frequency49 h_freq : float50 High pass frequency51 %(tstart_ecg)s52 %(filter_length_ecg)s53 %(verbose)s54 55 Returns56 -------57 events : array58 Indices of ECG peaks.59 """60 win_size = int(round((60.0 * sfreq) / 120.0))61 62 filtecg = filter_data(63 ecg,64 sfreq,65 l_freq,66 h_freq,67 None,68 filter_length,69 0.5,70 0.5,71 phase="zero-double",72 fir_window="hann",73 fir_design="firwin2",74 )75 76 ecg_abs = np.abs(filtecg)77 init = int(sfreq)78 79 n_samples_start = int(sfreq * tstart)80 ecg_abs = ecg_abs[n_samples_start:]81 82 n_points = len(ecg_abs)83 84 maxpt = np.empty(3)85 maxpt[0] = np.max(ecg_abs[:init])86 maxpt[1] = np.max(ecg_abs[init : init * 2])87 maxpt[2] = np.max(ecg_abs[init * 2 : init * 3])88 89 init_max = np.mean(maxpt)90 91 if thresh_value == "auto":92 thresh_runs = np.arange(0.3, 1.1, 0.05)93 elif isinstance(thresh_value, str):94 raise ValueError('threshold value must be "auto" or a float')95 else:96 thresh_runs = [thresh_value]97 98 # Try a few thresholds (or just one)99 clean_events = list()100 for thresh_value in thresh_runs:101 thresh1 = init_max * thresh_value102 numcross = list()103 time = list()104 rms = list()105 ii = 0106 while ii < (n_points - win_size):107 window = ecg_abs[ii : ii + win_size]108 if window[0] > thresh1:109 max_time = np.argmax(window)110 time.append(ii + max_time)111 nx = np.sum(112 np.diff(((window > thresh1).astype(np.int64) == 1).astype(int))113 )114 numcross.append(nx)115 rms.append(np.sqrt(sum_squared(window) / window.size))116 ii += win_size117 else:118 ii += 1119 120 if len(rms) == 0:121 rms.append(0.0)122 time.append(0.0)123 time = np.array(time)124 rms_mean = np.mean(rms)125 rms_std = np.std(rms)126 rms_thresh = rms_mean + (rms_std * levels)127 b = np.where(rms < rms_thresh)[0]128 a = np.array(numcross)[b]129 ce = time[b[a < n_thresh]]130 131 ce += n_samples_start132 if ce.size > 0: # We actually found an event133 clean_events.append(ce)134 135 if clean_events:136 # pick the best threshold; first get effective heart rates137 rates = np.array(138 [60.0 * len(cev) / (len(ecg) / float(sfreq)) for cev in clean_events]139 )140 141 # now find heart rates that seem reasonable (infant through adult142 # athlete)143 idx = np.where(np.logical_and(rates <= 160.0, rates >= 40.0))[0]144 if idx.size > 0:145 ideal_rate = np.median(rates[idx]) # get close to the median146 else:147 ideal_rate = 80.0 # get close to a reasonable default148 149 idx = np.argmin(np.abs(rates - ideal_rate))150 clean_events = clean_events[idx]151 else:152 clean_events = np.array([])153 154 return clean_events155 156 157@verbose158def find_ecg_events(159 raw,160 event_id=999,161 ch_name=None,162 tstart=0.0,163 l_freq=5,164 h_freq=35,165 qrs_threshold="auto",166 filter_length="10s",167 return_ecg=False,168 reject_by_annotation=True,169 verbose=None,170):171 """Find ECG events by localizing the R wave peaks.172 173 Parameters174 ----------175 raw : instance of Raw176 The raw data.177 %(event_id_ecg)s178 %(ch_name_ecg)s179 %(tstart_ecg)s180 %(l_freq_ecg_filter)s181 qrs_threshold : float | str182 Between 0 and 1. qrs detection threshold. Can also be "auto" to183 automatically choose the threshold that generates a reasonable184 number of heartbeats (40-160 beats / min).185 %(filter_length_ecg)s186 return_ecg : bool187 Return the ECG data. This is especially useful if no ECG channel188 is present in the input data, so one will be synthesized (only works if MEG189 channels are present in the data). Defaults to ``False``.190 %(reject_by_annotation_all)s191 192 .. versionadded:: 0.18193 %(verbose)s194 195 Returns196 -------197 ecg_events : array198 The events corresponding to the peaks of the R waves.199 ch_ecg : int | None200 Index of channel used.201 average_pulse : float202 The estimated average pulse. If no ECG events could be found, this will203 be zero.204 ecg : array | None205 The ECG data of the synthesized ECG channel, if any. This will only206 be returned if ``return_ecg=True`` was passed.207 208 See Also209 --------210 create_ecg_epochs211 compute_proj_ecg212 """213 skip_by_annotation = ("edge", "bad") if reject_by_annotation else ()214 del reject_by_annotation215 idx_ecg = _get_ecg_channel_index(ch_name, raw)216 if idx_ecg is not None:217 logger.info(f"Using channel {raw.ch_names[idx_ecg]} to identify heart beats.")218 ecg = raw.get_data(picks=idx_ecg)219 else:220 ecg, _ = _make_ecg(raw, start=None, stop=None)221 assert ecg.ndim == 2 and ecg.shape[0] == 1222 ecg = ecg[0]223 # Deal with filtering the same way we do in raw, i.e. filter each good224 # segment225 onsets, ends = _annotations_starts_stops(226 raw, skip_by_annotation, "reject_by_annotation", invert=True227 )228 ecgs = list()229 max_idx = (ends - onsets).argmax()230 for si, (start, stop) in enumerate(zip(onsets, ends)):231 # Only output filter params once (for info level), and only warn232 # once about the length criterion (longest segment is too short)233 use_verbose = verbose if si == max_idx else "error"234 ecgs.append(235 filter_data(236 ecg[start:stop],237 raw.info["sfreq"],238 l_freq,239 h_freq,240 [0],241 filter_length,242 0.5,243 0.5,244 1,245 "fir",246 None,247 copy=False,248 phase="zero-double",249 fir_window="hann",250 fir_design="firwin2",251 verbose=use_verbose,252 )253 )254 ecg = np.concatenate(ecgs)255 256 # detecting QRS and generating events. Since not user-controlled, don't257 # output filter params here (hardcode verbose=False)258 ecg_events = qrs_detector(259 raw.info["sfreq"],260 ecg,261 tstart=tstart,262 thresh_value=qrs_threshold,263 l_freq=None,264 h_freq=None,265 verbose=False,266 )267 268 # map ECG events back to original times269 remap = np.empty(len(ecg), int)270 offset = 0271 for start, stop in zip(onsets, ends):272 this_len = stop - start273 assert this_len >= 0274 remap[offset : offset + this_len] = np.arange(start, stop)275 offset += this_len276 assert offset == len(ecg)277 278 if ecg_events.size > 0:279 ecg_events = remap[ecg_events]280 else:281 ecg_events = np.array([])282 283 n_events = len(ecg_events)284 duration_sec = len(ecg) / raw.info["sfreq"] - tstart285 duration_min = duration_sec / 60.0286 average_pulse = n_events / duration_min287 logger.info(288 f"Number of ECG events detected : {n_events} "289 f"(average pulse {average_pulse} / min.)"290 )291 292 ecg_events = np.array(293 [294 ecg_events + raw.first_samp,295 np.zeros(n_events, int),296 event_id * np.ones(n_events, int),297 ]298 ).T299 300 out = (ecg_events, idx_ecg, average_pulse)301 ecg = ecg[np.newaxis] # backward compat output 2D302 if return_ecg:303 out += (ecg,)304 return out305 306 307def _get_ecg_channel_index(ch_name, inst):308 """Get ECG channel index, if no channel found returns None."""309 if ch_name is None:310 ecg_idx = pick_types(311 inst.info,312 meg=False,313 eeg=False,314 stim=False,315 eog=False,316 ecg=True,317 emg=False,318 ref_meg=False,319 exclude="bads",320 )321 else:322 if ch_name not in inst.ch_names:323 raise ValueError(f"{ch_name} not in channel list ({inst.ch_names})")324 ecg_idx = pick_channels(inst.ch_names, include=[ch_name])325 326 if len(ecg_idx) == 0:327 return None328 329 if len(ecg_idx) > 1:330 warn(331 f"More than one ECG channel found. Using only {inst.ch_names[ecg_idx[0]]}."332 )333 334 return ecg_idx[0]335 336 337@verbose338def create_ecg_epochs(339 raw,340 ch_name=None,341 event_id=999,342 picks=None,343 tmin=-0.5,344 tmax=0.5,345 l_freq=8,346 h_freq=16,347 reject=None,348 flat=None,349 baseline=None,350 preload=True,351 keep_ecg=False,352 reject_by_annotation=True,353 decim=1,354 verbose=None,355):356 """Conveniently generate epochs around ECG artifact events.357 358 %(create_ecg_epochs)s359 360 .. note:: Filtering is only applied to the ECG channel while finding361 events. The resulting ``ecg_epochs`` will have no filtering362 applied (i.e., have the same filter properties as the input363 ``raw`` instance).364 365 Parameters366 ----------367 raw : instance of Raw368 The raw data.369 %(ch_name_ecg)s370 %(event_id_ecg)s371 %(picks_all)s372 tmin : float373 Start time before event.374 tmax : float375 End time after event.376 %(l_freq_ecg_filter)s377 %(reject_epochs)s378 %(flat)s379 %(baseline_epochs)s380 preload : bool381 Preload epochs or not (default True). Must be True if382 keep_ecg is True.383 keep_ecg : bool384 When ECG is synthetically created (after picking), should it be added385 to the epochs? Must be False when synthetic channel is not used.386 Defaults to False.387 %(reject_by_annotation_epochs)s388 389 .. versionadded:: 0.14.0390 %(decim)s391 392 .. versionadded:: 0.21.0393 %(verbose)s394 395 Returns396 -------397 ecg_epochs : instance of Epochs398 Data epoched around ECG R wave peaks.399 400 See Also401 --------402 find_ecg_events403 compute_proj_ecg404 405 Notes406 -----407 If you already have a list of R-peak times, or want to compute R-peaks408 outside MNE-Python using a different algorithm, the recommended approach is409 to call the :class:`~mne.Epochs` constructor directly, with your R-peaks410 formatted as an :term:`events` array (here we also demonstrate the relevant411 default values)::412 413 mne.Epochs(raw, r_peak_events_array, tmin=-0.5, tmax=0.5,414 baseline=None, preload=True, proj=False) # doctest: +SKIP415 """416 has_ecg = "ecg" in raw or ch_name is not None417 if keep_ecg and (has_ecg or not preload):418 raise ValueError(419 "keep_ecg can be True only if the ECG channel is "420 "created synthetically and preload=True."421 )422 423 events, _, _, ecg = find_ecg_events(424 raw,425 ch_name=ch_name,426 event_id=event_id,427 l_freq=l_freq,428 h_freq=h_freq,429 return_ecg=True,430 reject_by_annotation=reject_by_annotation,431 )432 433 picks = _picks_to_idx(raw.info, picks, "all", exclude=())434 435 # create epochs around ECG events and baseline (important)436 ecg_epochs = Epochs(437 raw,438 events=events,439 event_id=event_id,440 tmin=tmin,441 tmax=tmax,442 proj=False,443 flat=flat,444 picks=picks,445 reject=reject,446 baseline=baseline,447 reject_by_annotation=reject_by_annotation,448 preload=preload,449 decim=decim,450 )451 452 if keep_ecg:453 # We know we have created a synthetic channel and epochs are preloaded454 ecg_raw = RawArray(455 ecg,456 create_info(457 ch_names=["ECG-SYN"], sfreq=raw.info["sfreq"], ch_types=["ecg"]458 ),459 first_samp=raw.first_samp,460 )461 with ecg_raw.info._unlock():462 ignore = ["ch_names", "chs", "nchan", "bads"]463 for k, v in raw.info.items():464 if k not in ignore:465 ecg_raw.info[k] = v466 syn_epochs = Epochs(467 ecg_raw,468 events=ecg_epochs.events,469 event_id=event_id,470 tmin=tmin,471 tmax=tmax,472 proj=False,473 picks=[0],474 baseline=baseline,475 decim=decim,476 preload=True,477 )478 ecg_epochs = ecg_epochs.add_channels([syn_epochs])479 480 return ecg_epochs481 482 483@verbose484def _make_ecg(inst, start, stop, reject_by_annotation=False, verbose=None):485 """Create ECG signal from cross channel average."""486 if not any(c in inst for c in ["mag", "grad"]):487 raise ValueError(488 "Generating an artificial ECG channel can only be done for MEG data."489 )490 for ch in ["mag", "grad"]:491 if ch in inst:492 break493 logger.info(494 "Reconstructing ECG signal from {}".format(495 {"mag": "Magnetometers", "grad": "Gradiometers"}[ch]496 )497 )498 picks = pick_types(inst.info, meg=ch, eeg=False, ref_meg=False)499 500 # Handle start/stop501 msg = (502 "integer arguments for the start and stop parameters are "503 "not supported for Epochs and Evoked objects. Please "504 "consider using float arguments specifying start and stop "505 "time in seconds."506 )507 begin_param_name = "tmin"508 if isinstance(start, int_like):509 if isinstance(inst, BaseRaw):510 # Raw has start param, can just use int511 begin_param_name = "start"512 else:513 raise ValueError(msg)514 515 end_param_name = "tmax"516 if isinstance(start, int_like):517 if isinstance(inst, BaseRaw):518 # Raw has stop param, can just use int519 end_param_name = "stop"520 else:521 raise ValueError(msg)522 523 kwargs = {begin_param_name: start, end_param_name: stop}524 525 if isinstance(inst, BaseRaw):526 reject_by_annotation = "omit" if reject_by_annotation else None527 ecg, times = inst.get_data(528 picks,529 return_times=True,530 **kwargs,531 reject_by_annotation=reject_by_annotation,532 )533 elif isinstance(inst, BaseEpochs):534 ecg = np.hstack(inst.copy().get_data(picks, **kwargs))535 times = inst.times536 elif isinstance(inst, Evoked):537 ecg = inst.get_data(picks, **kwargs)538 times = inst.times539 return ecg.mean(0, keepdims=True), times540