mne_compute_proj_eog.py

#!/usr/bin/env python
r"""Compute SSP/PCA projections for EOG artifacts.

You can do for example:

$ mne compute_proj_eog -i sample_audvis_raw.fif \
                       --l-freq 1 --h-freq 35 \
                       --rej-grad 3000 --rej-mag 4000 --rej-eeg 100

or

$ mne compute_proj_eog -i sample_audvis_raw.fif \
                       --l-freq 1 --h-freq 35 \
                       --rej-grad 3000 --rej-mag 4000 --rej-eeg 100 \
                       --proj sample_audvis_ecg-proj.fif

to exclude ECG artifacts from projection computation.
"""
from __future__ import print_function

# Authors : Alexandre Gramfort, Ph.D.
#           Martin Luessi, Ph.D.

from mne.externals.six import string_types
import os
import sys
import mne


def run():
    """Run command."""
    from mne.commands.utils import get_optparser

    parser = get_optparser(__file__)

    parser.add_option("-i", "--in", dest="raw_in",
                      help="Input raw FIF file", metavar="FILE")
    parser.add_option("--tmin", dest="tmin", type="float",
                      help="Time before event in seconds", default=-0.2)
    parser.add_option("--tmax", dest="tmax", type="float",
                      help="Time after event in seconds", default=0.2)
    parser.add_option("-g", "--n-grad", dest="n_grad", type="int",
                      help="Number of SSP vectors for gradiometers",
                      default=2)
    parser.add_option("-m", "--n-mag", dest="n_mag", type="int",
                      help="Number of SSP vectors for magnetometers",
                      default=2)
    parser.add_option("-e", "--n-eeg", dest="n_eeg", type="int",
                      help="Number of SSP vectors for EEG", default=2)
    parser.add_option("--l-freq", dest="l_freq", type="float",
                      help="Filter low cut-off frequency in Hz",
                      default=1)
    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)
    parser.add_option("-a", "--average", dest="average", action="store_true",
                      help="Compute SSP after averaging",
                      default=False)
    parser.add_option("--proj", dest="proj",
                      help="Use SSP projections from a fif file.",
                      default=None)
    parser.add_option("--filtersize", dest="filter_length", type="int",
                      help="Number of taps to use for filtering",
                      default=2048)
    parser.add_option("-j", "--n-jobs", dest="n_jobs", type="int",
                      help="Number of jobs to run in parallel", default=1)
    parser.add_option("--rej-grad", dest="rej_grad", type="float",
                      help="Gradiometers rejection parameter in fT/cm (peak "
                      "to peak amplitude)", default=2000)
    parser.add_option("--rej-mag", dest="rej_mag", type="float",
                      help="Magnetometers rejection parameter in fT (peak to "
                      "peak amplitude)", default=3000)
    parser.add_option("--rej-eeg", dest="rej_eeg", type="float",
                      help="EEG rejection parameter in uV (peak to peak "
                      "amplitude)", default=50)
    parser.add_option("--rej-eog", dest="rej_eog", type="float",
                      help="EOG rejection parameter in uV (peak to peak "
                      "amplitude)", default=1e9)
    parser.add_option("--avg-ref", dest="avg_ref", action="store_true",
                      help="Add EEG average reference proj",
                      default=False)
    parser.add_option("--no-proj", dest="no_proj", action="store_true",
                      help="Exclude the SSP projectors currently in the "
                      "fiff file",  default=False)
    parser.add_option("--bad", dest="bad_fname",
                      help="Text file containing bad channels list "
                      "(one per line)", default=None)
    parser.add_option("--event-id", dest="event_id", type="int",
                      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.)
    parser.add_option("-c", "--channel", dest="ch_name", type="string",
                      help="Custom EOG channel(s), comma separated",
                      default=None)

    options, args = parser.parse_args()

    raw_in = options.raw_in

    if raw_in is None:
        parser.print_help()
        sys.exit(1)

    tmin = options.tmin
    tmax = options.tmax
    n_grad = options.n_grad
    n_mag = options.n_mag
    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
    n_jobs = options.n_jobs
    reject = dict(grad=1e-13 * float(options.rej_grad),
                  mag=1e-15 * float(options.rej_mag),
                  eeg=1e-6 * float(options.rej_eeg),
                  eog=1e-6 * float(options.rej_eog))
    avg_ref = options.avg_ref
    no_proj = options.no_proj
    bad_fname = options.bad_fname
    event_id = options.event_id
    proj_fname = options.proj
    raw_event_fname = options.raw_event_fname
    tstart = options.tstart
    ch_name = options.ch_name

    if bad_fname is not None:
        with open(bad_fname, 'r') as fid:
            bads = [w.rstrip() for w in fid.readlines()]
        print('Bad channels read : %s' % bads)
    else:
        bads = []

    if raw_in.endswith('_raw.fif') or raw_in.endswith('-raw.fif'):
        prefix = raw_in[:-8]
    else:
        prefix = raw_in[:-4]

    eog_event_fname = prefix + '_eog-eve.fif'

    if average:
        eog_proj_fname = prefix + '_eog_avg-proj.fif'
    else:
        eog_proj_fname = prefix + '_eog-proj.fif'

    raw = mne.io.read_raw_fif(raw_in, preload=preload)

    if raw_event_fname is not None:
        raw_event = mne.io.read_raw_fif(raw_event_fname)
    else:
        raw_event = raw

    flat = None  # XXX : not exposed to the user
    projs, events = mne.preprocessing.compute_proj_eog(
        raw=raw, raw_event=raw_event, tmin=tmin, tmax=tmax, n_grad=n_grad,
        n_mag=n_mag, n_eeg=n_eeg, l_freq=l_freq, h_freq=h_freq,
        average=average, filter_length=filter_length,
        n_jobs=n_jobs, reject=reject, flat=flat, bads=bads,
        avg_ref=avg_ref, no_proj=no_proj, event_id=event_id,
        eog_l_freq=eog_l_freq, eog_h_freq=eog_h_freq,
        tstart=tstart, ch_name=ch_name, copy=False)

    raw.close()

    if raw_event_fname is not None:
        raw_event.close()

    if proj_fname is not None:
        print('Including SSP projections from : %s' % proj_fname)
        # append the eog projs, so they are last in the list
        projs = mne.read_proj(proj_fname) + projs

    if isinstance(preload, string_types) and os.path.exists(preload):
        os.remove(preload)

    print("Writing EOG projections in %s" % eog_proj_fname)
    mne.write_proj(eog_proj_fname, projs)

    print("Writing EOG events in %s" % eog_event_fname)
    mne.write_events(eog_event_fname, events)

is_main = (__name__ == '__main__')
if is_main:
    run()