[med-svn] [python-mne] 187/353: ENH : expose qrs_detector params in mne_compute_ecg_proj.py (add --tstart option)
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:55 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to tag 0.4
in repository python-mne.
commit fc767b8957a6cfc27180a4006d0b561276777169
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Thu May 24 19:10:22 2012 +0200
ENH : expose qrs_detector params in mne_compute_ecg_proj.py (add --tstart option)
---
bin/mne_compute_proj_ecg.py | 14 +++++++++++++-
bin/mne_compute_proj_eog.py | 14 +++++++++++++-
mne/artifacts/ecg.py | 36 +++++++++++++++++++++++++++-------
mne/artifacts/eog.py | 10 +++++++---
mne/preprocessing/ssp.py | 47 ++++++++++++++++++++++++++++++++++++++-------
5 files changed, 102 insertions(+), 19 deletions(-)
diff --git a/bin/mne_compute_proj_ecg.py b/bin/mne_compute_proj_ecg.py
index acdc104..b05cc0d 100755
--- a/bin/mne_compute_proj_ecg.py
+++ b/bin/mne_compute_proj_ecg.py
@@ -42,6 +42,12 @@ if __name__ == '__main__':
parser.add_option("--h-freq", dest="h_freq", type="float",
help="Filter high cut-off frequency in Hz",
default=100)
+ parser.add_option("--ecg-l-freq", dest="ecg_l_freq", type="float",
+ help="Filter low cut-off frequency in Hz used for ECG event detection",
+ default=5)
+ parser.add_option("--ecg-h-freq", dest="ecg_h_freq", type="float",
+ help="Filter high cut-off frequency in Hz used for ECG event detection",
+ default=35)
parser.add_option("-p", "--preload", dest="preload",
help="Temporary file used during computation (to save memory)",
default=True)
@@ -85,6 +91,8 @@ if __name__ == '__main__':
help="ID to use for events", default=999)
parser.add_option("--event-raw", dest="raw_event_fname",
help="raw file to use for event detection", default=None)
+ parser.add_option("--tstart", dest="tstart", type="float",
+ help="Start artifact detection after tstart seconds", default=0.)
options, args = parser.parse_args()
@@ -101,6 +109,8 @@ if __name__ == '__main__':
n_eeg = options.n_eeg
l_freq = options.l_freq
h_freq = options.h_freq
+ ecg_l_freq = options.ecg_l_freq
+ ecg_h_freq = options.ecg_h_freq
average = options.average
preload = options.preload
filter_length = options.filter_length
@@ -116,6 +126,7 @@ if __name__ == '__main__':
event_id = options.event_id
proj_fname = options.proj
raw_event_fname = options.raw_event_fname
+ tstart = options.tstart
if bad_fname is not None:
bads = [w.rstrip().split()[0] for w in open(bad_fname).readlines()]
@@ -146,7 +157,8 @@ if __name__ == '__main__':
tmin, tmax, n_grad, n_mag, n_eeg,
l_freq, h_freq, average, filter_length,
n_jobs, ch_name, reject,
- bads, avg_ref, no_proj, event_id)
+ bads, avg_ref, no_proj, event_id,
+ ecg_l_freq, ecg_h_freq, tstart)
raw.close()
diff --git a/bin/mne_compute_proj_eog.py b/bin/mne_compute_proj_eog.py
index 08c4968..770348d 100755
--- a/bin/mne_compute_proj_eog.py
+++ b/bin/mne_compute_proj_eog.py
@@ -48,6 +48,12 @@ if __name__ == '__main__':
parser.add_option("--h-freq", dest="h_freq", type="float",
help="Filter high cut-off frequency in Hz",
default=35)
+ parser.add_option("--eog-l-freq", dest="eog_l_freq", type="float",
+ help="Filter low cut-off frequency in Hz used for EOG event detection",
+ default=1)
+ parser.add_option("--eog-h-freq", dest="eog_h_freq", type="float",
+ help="Filter high cut-off frequency in Hz used for EOG event detection",
+ default=10)
parser.add_option("-p", "--preload", dest="preload",
help="Temporary file used during computation (to save memory)",
default=True)
@@ -88,6 +94,8 @@ if __name__ == '__main__':
help="ID to use for events", default=998)
parser.add_option("--event-raw", dest="raw_event_fname",
help="raw file to use for event detection", default=None)
+ parser.add_option("--tstart", dest="tstart", type="float",
+ help="Start artifact detection after tstart seconds", default=0.)
options, args = parser.parse_args()
@@ -104,6 +112,8 @@ if __name__ == '__main__':
n_eeg = options.n_eeg
l_freq = options.l_freq
h_freq = options.h_freq
+ eog_l_freq = options.eog_l_freq
+ eog_h_freq = options.eog_h_freq
average = options.average
preload = options.preload
filter_length = options.filter_length
@@ -118,6 +128,7 @@ if __name__ == '__main__':
event_id = options.event_id
proj_fname = options.proj
raw_event_fname = options.raw_event_fname
+ tstart = options.tstart
if bad_fname is not None:
bads = [w.rstrip().split()[0] for w in open(bad_fname).readlines()]
@@ -148,7 +159,8 @@ if __name__ == '__main__':
tmin, tmax, n_grad, n_mag, n_eeg,
l_freq, h_freq, average, filter_length,
n_jobs, reject, bads,
- avg_ref, no_proj, event_id)
+ avg_ref, no_proj, event_id,
+ eog_l_freq, eog_h_freq, tstart)
raw.close()
diff --git a/mne/artifacts/ecg.py b/mne/artifacts/ecg.py
index 18b7c43..456097a 100644
--- a/mne/artifacts/ecg.py
+++ b/mne/artifacts/ecg.py
@@ -5,7 +5,7 @@ from ..filter import band_pass_filter
def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3,
- low_pass=5, high_pass=35):
+ l_freq=5, h_freq=35, tstart=0):
"""Detect QRS component in ECG channels.
QRS is the main wave on the heart beat.
@@ -22,10 +22,12 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3,
number of std from mean to include for detection
n_thresh: int
max number of crossings
- low_pass: float
+ l_freq: float
Low pass frequency
- high_pass: float
+ h_freq: float
High pass frequency
+ tstart: float
+ Start detection after tstart seconds.
Returns
-------
@@ -34,12 +36,16 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3,
"""
win_size = round((60.0 * sfreq) / 120.0)
- filtecg = band_pass_filter(ecg, sfreq, low_pass, high_pass)
- n_points = len(filtecg)
+ filtecg = band_pass_filter(ecg, sfreq, l_freq, h_freq)
absecg = np.abs(filtecg)
init = int(sfreq)
+ n_samples_start = int(init * tstart)
+ absecg = absecg[n_samples_start:]
+
+ n_points = len(absecg)
+
maxpt = np.empty(3)
maxpt[0] = np.max(absecg[:init])
maxpt[1] = np.max(absecg[init:init * 2])
@@ -73,10 +79,13 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3,
a = np.array(numcross)[b]
clean_events = time[b[a < n_thresh]]
+ clean_events += n_samples_start
+
return clean_events
-def find_ecg_events(raw, event_id=999, ch_name=None):
+def find_ecg_events(raw, event_id=999, ch_name=None, tstart=0.0,
+ l_freq=5, h_freq=35, qrs_threshold=0.6):
"""Find ECG peaks
Parameters
@@ -89,6 +98,15 @@ def find_ecg_events(raw, event_id=999, ch_name=None):
The name of the channel to use for ECG peak detection.
The argument is mandatory if the dataset contains no ECG
channels.
+ tstart: float
+ Start detection after tstart seconds. Useful when beginning
+ of run is noisy.
+ l_freq: float
+ Low pass frequency
+ h_freq: float
+ High pass frequency
+ qrs_threshold: float
+ Between 0 and 1. qrs detection threshold
Returns
-------
@@ -122,7 +140,11 @@ def find_ecg_events(raw, event_id=999, ch_name=None):
ecg, times = raw[ch_ECG, :]
# detecting QRS and generating event file
- ecg_events = qrs_detector(info['sfreq'], ecg.ravel())
+ ecg_events = qrs_detector(info['sfreq'], ecg.ravel(), tstart=tstart,
+ thresh_value=qrs_threshold, l_freq=l_freq,
+ h_freq=h_freq)
+ import ipdb; ipdb.set_trace()
+
n_events = len(ecg_events)
average_pulse = n_events * 60.0 / (times[-1] - times[0])
print ("Number of ECG events detected : %d (average pulse %d / min.)"
diff --git a/mne/artifacts/eog.py b/mne/artifacts/eog.py
index 15ff8da..c9d7431 100644
--- a/mne/artifacts/eog.py
+++ b/mne/artifacts/eog.py
@@ -5,7 +5,7 @@ from .. import fiff
from ..filter import band_pass_filter
-def find_eog_events(raw, event_id=998):
+def find_eog_events(raw, event_id=998, l_freq=1, h_freq=10):
"""Locate EOG artifacts
Parameters
@@ -14,6 +14,10 @@ def find_eog_events(raw, event_id=998):
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
Returns
-------
@@ -49,8 +53,8 @@ def find_eog_events(raw, event_id=998):
indexmax = np.argmax(temp)
- # easy to detect peaks with this filtering.
- filteog = band_pass_filter(eog[indexmax], sampRate, 1, 10)
+ # easier to detect peaks with filtering.
+ filteog = band_pass_filter(eog[indexmax], sampRate, l_freq, h_freq)
# detecting eog blinks and generating event file
diff --git a/mne/preprocessing/ssp.py b/mne/preprocessing/ssp.py
index 711df3f..60713af 100644
--- a/mne/preprocessing/ssp.py
+++ b/mne/preprocessing/ssp.py
@@ -14,7 +14,8 @@ from ..artifacts import find_ecg_events, find_eog_events
def _compute_exg_proj(mode, raw, raw_event, tmin, tmax,
n_grad, n_mag, n_eeg, l_freq, h_freq,
average, filter_length, n_jobs, ch_name,
- reject, bads, avg_ref, no_proj, event_id):
+ reject, bads, avg_ref, no_proj, event_id,
+ exg_l_freq, exg_h_freq, tstart):
"""Compute SSP/PCA projections for ECG or EOG artifacts
Note: raw has to be constructed with preload=True (or string)
@@ -79,6 +80,15 @@ def _compute_exg_proj(mode, raw, raw_event, tmin, tmax,
event_id: int
ID to use for events
+ exg_l_freq: float
+ Low pass frequency applied for filtering EXG channel
+
+ exg_h_freq: float
+ High pass frequency applied for filtering EXG channel
+
+ tstart: float
+ Start artifact detection after tstart seconds.
+
Returns
-------
proj : list
@@ -108,10 +118,12 @@ def _compute_exg_proj(mode, raw, raw_event, tmin, tmax,
if mode == 'ECG':
print 'Running ECG SSP computation'
events, _, _ = find_ecg_events(raw_event, ch_name=ch_name,
- event_id=event_id)
+ event_id=event_id, l_freq=exg_l_freq,
+ h_freq=exg_h_freq, tstart=tstart)
elif mode == 'EOG':
print 'Running EOG SSP computation'
- events = find_eog_events(raw_event, event_id=event_id)
+ events = find_eog_events(raw_event, event_id=event_id,
+ l_freq=exg_l_freq, h_freq=exg_h_freq)
else:
ValueError("mode must be 'ECG' or 'EOG'")
@@ -162,7 +174,8 @@ def compute_proj_ecg(raw, raw_event=None, tmin=-0.2, tmax=0.4,
average=False, filter_length=2048, n_jobs=1, ch_name=None,
reject=dict(grad=2000e-13, mag=3000e-15, eeg=50e-6,
eog=250e-6), bads=[], avg_ref=False, no_proj=False,
- event_id=999):
+ event_id=999, ecg_l_freq=5, ecg_h_freq=35,
+ tstart=0.):
"""Compute SSP/PCA projections for ECG artifacts
Note: raw has to be constructed with preload=True (or string)
@@ -224,6 +237,15 @@ def compute_proj_ecg(raw, raw_event=None, tmin=-0.2, tmax=0.4,
event_id: int
ID to use for events
+ ecg_l_freq: float
+ Low pass frequency applied for filtering ECG channel
+
+ ecg_h_freq: float
+ High pass frequency applied for filtering ECG channel
+
+ tstart: float
+ Start artifact detection after tstart seconds.
+
Returns
-------
proj : list
@@ -236,7 +258,8 @@ def compute_proj_ecg(raw, raw_event=None, tmin=-0.2, tmax=0.4,
projs, ecg_events = _compute_exg_proj('ECG', raw, raw_event, tmin, tmax,
n_grad, n_mag, n_eeg, l_freq, h_freq,
average, filter_length, n_jobs, ch_name,
- reject, bads, avg_ref, no_proj, event_id)
+ reject, bads, avg_ref, no_proj, event_id,
+ ecg_l_freq, ecg_h_freq, tstart)
return projs, ecg_events
@@ -246,7 +269,7 @@ def compute_proj_eog(raw, raw_event=None, tmin=-0.2, tmax=0.2,
average=False, filter_length=2048, n_jobs=1,
reject=dict(grad=2000e-13, mag=3000e-15, eeg=500e-6,
eog=np.inf), bads=[], avg_ref=False, no_proj=False,
- event_id=998):
+ event_id=998, eog_l_freq=1, eog_h_freq=10, tstart=0.):
"""Compute SSP/PCA projections for EOG artifacts
Note: raw has to be constructed with preload=True (or string)
@@ -308,6 +331,15 @@ def compute_proj_eog(raw, raw_event=None, tmin=-0.2, tmax=0.2,
event_id: int
ID to use for events
+ eog_l_freq: float
+ Low pass frequency applied for filtering E0G channel
+
+ eog_h_freq: float
+ High pass frequency applied for filtering E0G channel
+
+ tstart: float
+ Start artifact detection after tstart seconds.
+
Returns
-------
proj : list
@@ -320,6 +352,7 @@ def compute_proj_eog(raw, raw_event=None, tmin=-0.2, tmax=0.2,
projs, eog_events = _compute_exg_proj('EOG', raw, raw_event, tmin, tmax,
n_grad, n_mag, n_eeg, l_freq, h_freq,
average, filter_length, n_jobs, None,
- reject, bads, avg_ref, no_proj, event_id)
+ reject, bads, avg_ref, no_proj, event_id,
+ eog_l_freq, eog_h_freq, tstart)
return projs, eog_events
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git
More information about the debian-med-commit
mailing list