Python mne.compute_raw_data_covariance函数代码示例

def test_cov_estimation_on_raw_segment():
    """Test estimation from raw on continuous recordings (typically empty room)
    """
    tempdir = _TempDir()
    raw = Raw(raw_fname, preload=False)
    cov = compute_raw_data_covariance(raw)
    cov_mne = read_cov(erm_cov_fname)
    assert_true(cov_mne.ch_names == cov.ch_names)
    assert_true(linalg.norm(cov.data - cov_mne.data, ord='fro')
                / linalg.norm(cov.data, ord='fro') < 1e-4)

    # test IO when computation done in Python
    cov.save(op.join(tempdir, 'test-cov.fif'))  # test saving
    cov_read = read_cov(op.join(tempdir, 'test-cov.fif'))
    assert_true(cov_read.ch_names == cov.ch_names)
    assert_true(cov_read.nfree == cov.nfree)
    assert_array_almost_equal(cov.data, cov_read.data)

    # test with a subset of channels
    picks = pick_channels(raw.ch_names, include=raw.ch_names[:5])
    cov = compute_raw_data_covariance(raw, picks=picks)
    assert_true(cov_mne.ch_names[:5] == cov.ch_names)
    assert_true(linalg.norm(cov.data - cov_mne.data[picks][:, picks],
                ord='fro') / linalg.norm(cov.data, ord='fro') < 1e-4)
    # make sure we get a warning with too short a segment
    raw_2 = raw.crop(0, 1)
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        cov = compute_raw_data_covariance(raw_2)
    assert_true(len(w) == 1)

def test_cov_estimation_on_raw_segment():
    """Estimate raw on continuous recordings (typically empty room)
    """
    raw = Raw(raw_fname)
    cov = compute_raw_data_covariance(raw)
    cov_mne = read_cov(erm_cov_fname)
    assert_true(cov_mne.ch_names == cov.ch_names)
    print (linalg.norm(cov.data - cov_mne.data, ord='fro')
            / linalg.norm(cov.data, ord='fro'))
    assert_true(linalg.norm(cov.data - cov_mne.data, ord='fro')
            / linalg.norm(cov.data, ord='fro')) < 1e-6

    # test IO when computation done in Python
    cov.save('test-cov.fif')  # test saving
    cov_read = read_cov('test-cov.fif')
    assert_true(cov_read.ch_names == cov.ch_names)
    assert_true(cov_read.nfree == cov.nfree)
    assert_true((linalg.norm(cov.data - cov_read.data, ord='fro')
            / linalg.norm(cov.data, ord='fro')) < 1e-5)

    # test with a subset of channels
    picks = pick_channels(raw.ch_names, include=raw.ch_names[:5])
    cov = compute_raw_data_covariance(raw, picks=picks)
    assert_true(cov_mne.ch_names[:5] == cov.ch_names)
    assert_true(linalg.norm(cov.data - cov_mne.data[picks][:, picks],
                ord='fro') / linalg.norm(cov.data, ord='fro')) < 1e-6

def get_epochs_and_cov(X, y, window=500):
    """return epochs from array."""
    raw_train = toMNE(X, y)
    picks = range(len(getChannelNames()))

    events = list()
    events_id = dict()
    for j, eid in enumerate(getEventNames()):
        tmp = find_events(raw_train, stim_channel=eid, verbose=False)
        tmp[:, -1] = j + 1
        events.append(tmp)
        events_id[eid] = j + 1

    # concatenate and sort events
    events = np.concatenate(events, axis=0)
    order_ev = np.argsort(events[:, 0])
    events = events[order_ev]

    epochs = Epochs(raw_train, events, events_id,
                    tmin=-(window / 500.0) + 1 / 500.0 + 0.150,
                    tmax=0.150, proj=False, picks=picks, baseline=None,
                    preload=True, add_eeg_ref=False, verbose=False)

    cov_signal = compute_raw_data_covariance(raw_train, verbose=False)
    return epochs, cov_signal

def INVERSE(wdir, Subject, epoch_info, evokeds):

    # import parameters from configuration file
    from configuration import ( lambda2, method )

    # compute noise covariance from empty room data
    emptyroom_raw = mne.io.Raw(wdir + '/data/maxfilter/' + Subject + '/'+ Subject +'_empty_sss.fif')  
    noise_cov     = mne.compute_raw_data_covariance(emptyroom_raw)
    
    # compute dSPM solution
    fname_fwd     = wdir + '/data/forward/' + Subject + '/' + Subject + '_phase1_trans_sss_filt140_raw-ico5-fwd.fif'
    forward       = mne.read_forward_solution(fname_fwd, surf_ori=True)
    
    # create inverse operator
    inverse_operator = make_inverse_operator(epoch_info, forward, noise_cov, loose=0.4, depth=0.8)
    
    # Compute inverse solution
    stcs = []
    for evoked in evokeds:
        stcs.append(apply_inverse(evoked, inverse_operator, lambda2, method=method, pick_ori = None))
    
    # save a covariance picture for visual inspection
    mne.viz.plot_cov(noise_cov, epoch_info, colorbar=True, proj=True,show_svd=False,show=False)
    plt.savefig(wdir + "/plots/" + Subject + "_covmat")
    plt.close()
    
    return stcs

def test_lcmv_raw():
    """Test LCMV with raw data
    """
    raw, _, _, _, noise_cov, label, forward, _, _, _ =\
        _get_data(all_forward=False, epochs=False, data_cov=False)

    tmin, tmax = 0, 20
    start, stop = raw.time_as_index([tmin, tmax])

    # use only the left-temporal MEG channels for LCMV
    left_temporal_channels = mne.read_selection('Left-temporal')
    picks = mne.fiff.pick_types(raw.info, meg=True, exclude='bads',
                                selection=left_temporal_channels)

    data_cov = mne.compute_raw_data_covariance(raw, tmin=tmin, tmax=tmax)

    stc = lcmv_raw(raw, forward, noise_cov, data_cov, reg=0.01, label=label,
                   start=start, stop=stop, picks=picks)

    assert_array_almost_equal(np.array([tmin, tmax]),
                              np.array([stc.times[0], stc.times[-1]]),
                              decimal=2)

    # make sure we get an stc with vertices only in the lh
    vertno = [forward['src'][0]['vertno'], forward['src'][1]['vertno']]
    assert_true(len(stc.vertno[0]) == len(np.intersect1d(vertno[0],
                                                         label.vertices)))
    assert_true(len(stc.vertno[1]) == 0)

def test_maxwell_filter_additional():
    """Test processing of Maxwell filtered data"""

    # TODO: Future tests integrate with mne/io/tests/test_proc_history

    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    data_path = op.join(testing.data_path(download=False))

    file_name = 'test_move_anon'

    raw_fname = op.join(data_path, 'SSS', file_name + '_raw.fif')

    with warnings.catch_warnings(record=True):  # maxshield
        # Use 2.0 seconds of data to get stable cov. estimate
        raw = Raw(raw_fname, preload=False, proj=False,
                  allow_maxshield=True).crop(0., 2., False)

    # Get MEG channels, compute Maxwell filtered data
    raw.preload_data()
    raw.pick_types(meg=True, eeg=False)
    int_order, ext_order = 8, 3
    raw_sss = maxwell.maxwell_filter(raw, int_order=int_order,
                                     ext_order=ext_order)

    # Test io on processed data
    tempdir = _TempDir()
    test_outname = op.join(tempdir, 'test_raw_sss.fif')
    raw_sss.save(test_outname)
    raw_sss_loaded = Raw(test_outname, preload=True, proj=False,
                         allow_maxshield=True)

    # Some numerical imprecision since save uses 'single' fmt
    assert_allclose(raw_sss_loaded._data[:, :], raw_sss._data[:, :],
                    rtol=1e-6, atol=1e-20)

    # Test rank of covariance matrices for raw and SSS processed data
    cov_raw = compute_raw_data_covariance(raw)
    cov_sss = compute_raw_data_covariance(raw_sss)

    scalings = None
    cov_raw_rank = _estimate_rank_meeg_cov(cov_raw['data'], raw.info, scalings)
    cov_sss_rank = _estimate_rank_meeg_cov(cov_sss['data'], raw_sss.info,
                                           scalings)

    assert_equal(cov_raw_rank, raw.info['nchan'])
    assert_equal(cov_sss_rank, maxwell.get_num_moments(int_order, 0))

def test_cov_estimation_on_raw_segment():
    """Estimate raw on continuous recordings (typically empty room)
    """
    raw = Raw(raw_fname)
    cov = mne.compute_raw_data_covariance(raw)
    cov_mne = mne.Covariance(erm_cov_fname)
    assert cov_mne.ch_names == cov.ch_names
    print (linalg.norm(cov.data - cov_mne.data, ord='fro')
            / linalg.norm(cov.data, ord='fro'))
    assert (linalg.norm(cov.data - cov_mne.data, ord='fro')
            / linalg.norm(cov.data, ord='fro')) < 1e-6

    def _fit_xdawn(self, raw, epochs=None):
        # Epochs can be computed from raw info, but might have been computed already
        if epochs is None:
            epochs = self._epochs(raw, self._events())

        sig_cov = mne.compute_raw_data_covariance(raw, picks=self._picks(raw), verbose=False).data
        evoked_cov = np.cov(epochs.average().data)
        evals, evecs = eigh(evoked_cov, sig_cov)
        # sort eigenvectors by their corresponding eigenvalues and normalize them
        evecs = evecs[:, np.argsort(evals)[::-1]]
        self._V = evecs / np.sqrt(np.sum(evecs ** 2, axis=0))
        self._A = inv(self._V.T)

def comp_noise(sub_id):
    """ This function compute and save the noise covariance matrix from
    the empty room recordings
    """

    fname = "sub_%d_empty_room.fif" % sub_id
    outname = "sub_%d_empty_room-cov.fif" % sub_id
    raw = mne.fiff.Raw(fname, preload=True)

    picks = mne.fiff.pick_types(raw.info, meg=True, eeg=True, stim=False,
                                eog=True, exclude='bads')
    cov = mne.compute_raw_data_covariance(raw, picks=picks, reject=None)

    cov.save(outname)

def localize_activity(subj, tmax=np.Inf, clean_only=True, reg=0):
    import find_good_segments as fgs

    data_path = env.data + '/MEG_data/fifs/'
    fwd_path = env.data + '/MEG_data/analysis/rest/'
    fwd_fname = fwd_path + subj + '_rest_LP100_HP0.6_CP3_DS300_raw-5-fwd.fif'

    # preloading makes computing the covariance a lot faster
    raw = mne.fiff.Raw(data_path + subj + '_rest_LP100_HP0.6_CP3_DS300_raw.fif', preload=True)

    picks = mne.fiff.pick_channels_regexp(raw.info['ch_names'], 'M..-*')

    # we don't need to choose picks because we only work with MEG channels, and
    # all channels are good
    fwd = mne.read_forward_solution(fwd_fname)

    # Load real data as templates
    if clean_only:
        start, end, num_chans = fgs.find_good_segments(subj, threshold=3500e-15)
        if start == 0:
            start = start + 3

        if end - start > tmax:
            start = end - tmax

        cov = mne.compute_raw_data_covariance(raw, tmin=start, tmax=end, picks=picks)
    else:
        cov = mne.compute_raw_data_covariance(raw, picks=picks)

    weights = calculate_weights(fwd, cov, reg=reg)

    data, times = raw[picks, raw.time_as_index(start):raw.time_as_index(end)]
    print 'Multiplying data by beamformer weights...'
    sol = np.dot(weights, data)
    src = mne.SourceEstimate(sol, [fwd['src'][0]['vertno'], fwd['src'][1]['vertno']], times[0], times[1] - times[0])
    return src

def test_xdawn_fit():
    """Test Xdawn fit."""
    # get data
    raw, events, picks = _get_data()
    epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    preload=True, baseline=None, verbose=False)
    # =========== Basic Fit test =================
    # test base xdawn
    xd = Xdawn(n_components=2, correct_overlap='auto',
               signal_cov=None, reg=None)
    xd.fit(epochs)
    # with this parameters, the overlapp correction must be False
    assert_equal(xd.correct_overlap, False)
    # no overlapp correction should give averaged evoked
    evoked = epochs['cond2'].average()
    assert_array_equal(evoked.data, xd.evokeds_['cond2'].data)

    # ========== with signal cov provided ====================
    # provide covariance object
    signal_cov = compute_raw_data_covariance(raw, picks=picks)
    xd = Xdawn(n_components=2, correct_overlap=False,
               signal_cov=signal_cov, reg=None)
    xd.fit(epochs)
    # provide ndarray
    signal_cov = np.eye(len(picks))
    xd = Xdawn(n_components=2, correct_overlap=False,
               signal_cov=signal_cov, reg=None)
    xd.fit(epochs)
    # provide ndarray of bad shape
    signal_cov = np.eye(len(picks) - 1)
    xd = Xdawn(n_components=2, correct_overlap=False,
               signal_cov=signal_cov, reg=None)
    assert_raises(ValueError, xd.fit, epochs)
    # provide another type
    signal_cov = 42
    xd = Xdawn(n_components=2, correct_overlap=False,
               signal_cov=signal_cov, reg=None)
    assert_raises(ValueError, xd.fit, epochs)
    # fit with baseline correction and ovverlapp correction should throw an
    # error
    epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    preload=True, baseline=(None, 0), verbose=False)

    xd = Xdawn(n_components=2, correct_overlap=True)
    assert_raises(ValueError, xd.fit, epochs)

def test_lcmv_raw():
    """Test LCMV with raw data
    """
    forward = mne.read_forward_solution(fname_fwd)
    label = mne.read_label(fname_label)
    noise_cov = mne.read_cov(fname_cov)
    raw = mne.fiff.Raw(fname_raw, preload=False)

    tmin, tmax = 0, 20
    # Setup for reading the raw data
    raw.info['bads'] = ['MEG 2443', 'EEG 053']  # 2 bads channels

    # Set up pick list: EEG + MEG - bad channels (modify to your needs)
    left_temporal_channels = mne.read_selection('Left-temporal')
    picks = mne.fiff.pick_types(raw.info, meg=True, eeg=False, stim=True,
                                eog=True, exclude='bads',
                                selection=left_temporal_channels)

    noise_cov = mne.read_cov(fname_cov)
    noise_cov = mne.cov.regularize(noise_cov, raw.info,
                                   mag=0.05, grad=0.05, eeg=0.1, proj=True)

    start, stop = raw.time_as_index([tmin, tmax])

    # use only the left-temporal MEG channels for LCMV
    picks = mne.fiff.pick_types(raw.info, meg=True, exclude='bads',
                                selection=left_temporal_channels)

    data_cov = mne.compute_raw_data_covariance(raw, tmin=tmin, tmax=tmax)

    stc = lcmv_raw(raw, forward, noise_cov, data_cov, reg=0.01, label=label,
                   start=start, stop=stop, picks=picks)

    assert_array_almost_equal(np.array([tmin, tmax]),
                              np.array([stc.times[0], stc.times[-1]]),
                              decimal=2)

    # make sure we get an stc with vertices only in the lh
    vertno = [forward['src'][0]['vertno'], forward['src'][1]['vertno']]
    assert_true(len(stc.vertno[0]) == len(np.intersect1d(vertno[0],
                                                         label.vertices)))
    assert_true(len(stc.vertno[1]) == 0)

def create_noise_covariance_matrix(fname_empty_room, fname_out=None, verbose=None):
    """Creates the noise covariance matrix from an empty room file"""

    print ">>>> estimate noise covariance matrix from empty room file..."

    # read in data
    raw_empty = mne.fiff.Raw(fname_empty_room,
                             verbose=verbose)
    # filter data

    # pick only MEG channels
    picks = mne.fiff.pick_types(raw_empty.info,
                                meg=True,
                                exclude='bads')

    # calculate noise-covariance matrix
    noise_cov_mat = mne.compute_raw_data_covariance(raw_empty,
                                                    picks=picks,
                                                    verbose=verbose)
    # write noise-covariance matrix to disk
    if fname_out is not None:
        mne.write_cov(fname_out, noise_cov_mat)

    return noise_cov_mat

Estimate covariance matrix from a raw FIF file
==============================================

"""
# Author: Alexandre Gramfort <[email protected]>
#
# License: BSD (3-clause)

print __doc__

import mne
from mne import fiff
from mne.datasets import sample

data_path = sample.data_path('.')
fname = data_path + '/MEG/sample/sample_audvis_raw.fif'

raw = fiff.Raw(fname)

# Set up pick list: MEG + STI 014 - bad channels
cov = mne.compute_raw_data_covariance(raw, reject=dict(eeg=40e-6, eog=150e-6))
print cov

###############################################################################
# Show covariance
import pylab as pl
pl.figure()
pl.imshow(cov.data, interpolation="nearest")
pl.title('Full covariance matrix')
pl.show()

# License: BSD (3-clause)

print __doc__

import mne
from mne import fiff
from mne.datasets import sample

data_path = sample.data_path()
fname = data_path + '/MEG/sample/sample_audvis_raw.fif'

raw = fiff.Raw(fname)

include = []  # or stim channels ['STI 014']
raw.info['bads'] += ['EEG 053']  # bads + 1 more

# pick EEG channels
picks = fiff.pick_types(raw.info, meg=True, eeg=True, stim=False, eog=True,
                                            include=include, exclude='bads')
# setup rejection
reject = dict(eeg=80e-6, eog=150e-6)

# Compute the covariance from the raw data
cov = mne.compute_raw_data_covariance(raw, picks=picks, reject=reject)
print cov

###############################################################################
# Show covariance
mne.viz.plot_cov(cov, raw.info, colorbar=True, proj=True)
# try setting proj to False to see the effect