Python ICA.find_bads_ecg方法代码示例

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]
def test_ica_labels():
    """Test ICA labels."""
    # The CTF data are uniquely well suited to testing the ICA.find_bads_
    # methods
    raw = read_raw_ctf(ctf_fname, preload=True)
    # derive reference ICA components and append them to raw
    icarf = ICA(n_components=2, random_state=0, max_iter=2, allow_ref_meg=True)
    with pytest.warns(UserWarning, match='did not converge'):
        icarf.fit(raw.copy().pick_types(meg=False, ref_meg=True))
    icacomps = icarf.get_sources(raw)
    # rename components so they are auto-detected by find_bads_ref
    icacomps.rename_channels({c: 'REF_' + c for c in icacomps.ch_names})
    # and add them to raw
    raw.add_channels([icacomps])
    # set the appropriate EEG channels to EOG and ECG
    raw.set_channel_types({'EEG057': 'eog', 'EEG058': 'eog', 'EEG059': 'ecg'})
    ica = ICA(n_components=4, random_state=0, max_iter=2, method='fastica')
    with pytest.warns(UserWarning, match='did not converge'):
        ica.fit(raw)

    ica.find_bads_eog(raw, l_freq=None, h_freq=None)
    picks = list(pick_types(raw.info, meg=False, eog=True))
    for idx, ch in enumerate(picks):
        assert '{}/{}/{}'.format('eog', idx, raw.ch_names[ch]) in ica.labels_
    assert 'eog' in ica.labels_
    for key in ('ecg', 'ref_meg', 'ecg/ECG-MAG'):
        assert key not in ica.labels_

    ica.find_bads_ecg(raw, l_freq=None, h_freq=None, method='correlation')
    picks = list(pick_types(raw.info, meg=False, ecg=True))
    for idx, ch in enumerate(picks):
        assert '{}/{}/{}'.format('ecg', idx, raw.ch_names[ch]) in ica.labels_
    for key in ('ecg', 'eog'):
        assert key in ica.labels_
    for key in ('ref_meg', 'ecg/ECG-MAG'):
        assert key not in ica.labels_

    ica.find_bads_ref(raw, l_freq=None, h_freq=None)
    picks = pick_channels_regexp(raw.ch_names, 'REF_ICA*')
    for idx, ch in enumerate(picks):
        assert '{}/{}/{}'.format('ref_meg', idx,
                                 raw.ch_names[ch]) in ica.labels_
    for key in ('ecg', 'eog', 'ref_meg'):
        assert key in ica.labels_
    assert 'ecg/ECG-MAG' not in ica.labels_

    ica.find_bads_ecg(raw, l_freq=None, h_freq=None)
    for key in ('ecg', 'eog', 'ref_meg', 'ecg/ECG-MAG'):
        assert key in ica.labels_

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]

#.........这里部分代码省略.........
        assert_true(isinstance(ica_read.info, Info))

        sources = ica.get_sources(raw)[:, :][0]
        sources2 = ica_read.get_sources(raw)[:, :][0]
        assert_array_almost_equal(sources, sources2)

        _raw1 = ica.apply(raw, exclude=[1])
        _raw2 = ica_read.apply(raw, exclude=[1])
        assert_array_almost_equal(_raw1[:, :][0], _raw2[:, :][0])

    os.remove(test_ica_fname)
    # check scrore funcs
    for name, func in get_score_funcs().items():
        if name in score_funcs_unsuited:
            continue
        scores = ica.score_sources(raw, target='EOG 061', score_func=func,
                                   start=0, stop=10)
        assert_true(ica.n_components_ == len(scores))

    # check univariate stats
    scores = ica.score_sources(raw, score_func=stats.skew)
    # check exception handling
    assert_raises(ValueError, ica.score_sources, raw,
                  target=np.arange(1))

    params = []
    params += [(None, -1, slice(2), [0, 1])]  # varicance, kurtosis idx params
    params += [(None, 'MEG 1531')]  # ECG / EOG channel params
    for idx, ch_name in product(*params):
        ica.detect_artifacts(raw, start_find=0, stop_find=50, ecg_ch=ch_name,
                             eog_ch=ch_name, skew_criterion=idx,
                             var_criterion=idx, kurt_criterion=idx)
    with warnings.catch_warnings(record=True):
        idx, scores = ica.find_bads_ecg(raw, method='ctps')
        assert_equal(len(scores), ica.n_components_)
        idx, scores = ica.find_bads_ecg(raw, method='correlation')
        assert_equal(len(scores), ica.n_components_)
        idx, scores = ica.find_bads_ecg(epochs, method='ctps')
        assert_equal(len(scores), ica.n_components_)
        assert_raises(ValueError, ica.find_bads_ecg, epochs.average(),
                      method='ctps')
        assert_raises(ValueError, ica.find_bads_ecg, raw,
                      method='crazy-coupling')

        idx, scores = ica.find_bads_eog(raw)
        assert_equal(len(scores), ica.n_components_)
        raw.info['chs'][raw.ch_names.index('EOG 061') - 1]['kind'] = 202
        idx, scores = ica.find_bads_eog(raw)
        assert_true(isinstance(scores, list))
        assert_equal(len(scores[0]), ica.n_components_)

    # check score funcs
    for name, func in get_score_funcs().items():
        if name in score_funcs_unsuited:
            continue
        scores = ica.score_sources(epochs_eog, target='EOG 061',
                                   score_func=func)
        assert_true(ica.n_components_ == len(scores))

    # check univariate stats
    scores = ica.score_sources(epochs, score_func=stats.skew)

    # check exception handling
    assert_raises(ValueError, ica.score_sources, epochs,
                  target=np.arange(1))

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]
# to definitely register this component as a bad one to be removed
# there is the ``ica.exclude`` attribute, a simple Python list
ica.exclude.extend(eog_inds)

# from now on the ICA will reject this component even if no exclude
# parameter is passed, and this information will be stored to disk
# on saving

# uncomment this for reading and writing
# ica.save('my-ica.fif')
# ica = read_ica('my-ica.fif')

###############################################################################
# Exercise: find and remove ECG artifacts using ICA!
ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5)
ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
ica.plot_properties(ecg_epochs, picks=ecg_inds, psd_args={'fmax': 35.})


###############################################################################
# What if we don't have an EOG channel?
# -------------------------------------
#
# We could either:
#
# 1. make a bipolar reference from frontal EEG sensors and use as virtual EOG
#    channel. This can be tricky though as you can only hope that the frontal
#    EEG channels only reflect EOG and not brain dynamics in the prefrontal
#    cortex.
# 2. go for a semi-automated approach, using template matching.
#

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]
title = 'Sources related to %s artifacts (red)'

# generate ECG epochs use detection via phase statistics

picks = mne.pick_types(raw.info, meg=True, eeg=False, eog=True, ecg=True,
                       stim=False, exclude='bads')

# create_ecg_epochs is strange: it strips the channels of anything non M/EEG
# UNLESS picks=None
picks=None
ecg_epochs = create_ecg_epochs(raw, ch_name='ECG002', tmin=-.5, tmax=.5, picks=picks, verbose=True)

# This will work with the above, but uses MASSIVE RAM
# Not sure the ECG quality is good enough for the QRS-detector
ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps', ch_name='ECG002')
# This creates a synthetic ECG from magnetometers, probably better...
ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps', ch_name='ECG002')

fig = ica.plot_scores(scores, exclude=ecg_inds, title=title % 'ecg')
fig.savefig(img_folder + '/ica_ecg_scores.png')

show_picks = np.abs(scores).argsort()[::-1][:5]

fig = ica.plot_sources(raw, show_picks, exclude=ecg_inds, title=title % 'ecg')
fig.savefig(img_folder + '/ica_ecg_sources.png')
fig = ica.plot_components(ecg_inds, title=title % 'ecg', colorbar=True)
fig.set_size_inches(12.,8.)
fig.savefig(img_folder + '/ica_ecg_components.png')

ecg_inds = ecg_inds[:n_max_ecg]

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]
#
# License: BSD (3-clause)

import mne
from mne.io import Raw
from mne.preprocessing import ICA, create_ecg_epochs
from mne.datasets import sample

print(__doc__)

###############################################################################
# Fit ICA model using the FastICA algorithm, detect and inspect components

data_path = sample.data_path()
raw_fname = data_path + "/MEG/sample/sample_audvis_filt-0-40_raw.fif"

raw = Raw(raw_fname, preload=True)
raw.filter(1, 30, method="iir")
raw.pick_types(meg=True, eeg=False, exclude="bads", stim=True)

# longer + more epochs for more artifact exposure
events = mne.find_events(raw, stim_channel="STI 014")
epochs = mne.Epochs(raw, events, event_id=None, tmin=-0.2, tmax=0.5)

ica = ICA(n_components=0.95, method="fastica").fit(epochs)

ecg_epochs = create_ecg_epochs(raw, tmin=-0.5, tmax=0.5)
ecg_inds, scores = ica.find_bads_ecg(ecg_epochs)

ica.plot_components(ecg_inds)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]
def test_ica_additional(method):
    """Test additional ICA functionality."""
    _skip_check_picard(method)

    tempdir = _TempDir()
    stop2 = 500
    raw = read_raw_fif(raw_fname).crop(1.5, stop).load_data()
    raw.del_proj()  # avoid warnings
    raw.set_annotations(Annotations([0.5], [0.5], ['BAD']))
    # XXX This breaks the tests :(
    # raw.info['bads'] = [raw.ch_names[1]]
    test_cov = read_cov(test_cov_name)
    events = read_events(event_name)
    picks = pick_types(raw.info, meg=True, stim=False, ecg=False,
                       eog=False, exclude='bads')[1::2]
    epochs = Epochs(raw, events, None, tmin, tmax, picks=picks,
                    baseline=(None, 0), preload=True, proj=False)
    epochs.decimate(3, verbose='error')
    assert len(epochs) == 4

    # test if n_components=None works
    ica = ICA(n_components=None, max_pca_components=None,
              n_pca_components=None, random_state=0, method=method, max_iter=1)
    with pytest.warns(UserWarning, match='did not converge'):
        ica.fit(epochs)
    # for testing eog functionality
    picks2 = np.concatenate([picks, pick_types(raw.info, False, eog=True)])
    epochs_eog = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks2,
                        baseline=(None, 0), preload=True)
    del picks2

    test_cov2 = test_cov.copy()
    ica = ICA(noise_cov=test_cov2, n_components=3, max_pca_components=4,
              n_pca_components=4, method=method)
    assert (ica.info is None)
    with pytest.warns(RuntimeWarning, match='normalize_proj'):
        ica.fit(raw, picks[:5])
    assert (isinstance(ica.info, Info))
    assert (ica.n_components_ < 5)

    ica = ICA(n_components=3, max_pca_components=4, method=method,
              n_pca_components=4, random_state=0)
    pytest.raises(RuntimeError, ica.save, '')

    ica.fit(raw, picks=[1, 2, 3, 4, 5], start=start, stop=stop2)

    # check passing a ch_name to find_bads_ecg
    with pytest.warns(RuntimeWarning, match='longer'):
        _, scores_1 = ica.find_bads_ecg(raw)
        _, scores_2 = ica.find_bads_ecg(raw, raw.ch_names[1])
    assert scores_1[0] != scores_2[0]

    # test corrmap
    ica2 = ica.copy()
    ica3 = ica.copy()
    corrmap([ica, ica2], (0, 0), threshold='auto', label='blinks', plot=True,
            ch_type="mag")
    corrmap([ica, ica2], (0, 0), threshold=2, plot=False, show=False)
    assert (ica.labels_["blinks"] == ica2.labels_["blinks"])
    assert (0 in ica.labels_["blinks"])
    # test retrieval of component maps as arrays
    components = ica.get_components()
    template = components[:, 0]
    EvokedArray(components, ica.info, tmin=0.).plot_topomap([0], time_unit='s')

    corrmap([ica, ica3], template, threshold='auto', label='blinks', plot=True,
            ch_type="mag")
    assert (ica2.labels_["blinks"] == ica3.labels_["blinks"])

    plt.close('all')

    ica_different_channels = ICA(n_components=2, random_state=0).fit(
        raw, picks=[2, 3, 4, 5])
    pytest.raises(ValueError, corrmap, [ica_different_channels, ica], (0, 0))

    # test warnings on bad filenames
    ica_badname = op.join(op.dirname(tempdir), 'test-bad-name.fif.gz')
    with pytest.warns(RuntimeWarning, match='-ica.fif'):
        ica.save(ica_badname)
    with pytest.warns(RuntimeWarning, match='-ica.fif'):
        read_ica(ica_badname)

    # test decim
    ica = ICA(n_components=3, max_pca_components=4,
              n_pca_components=4, method=method, max_iter=1)
    raw_ = raw.copy()
    for _ in range(3):
        raw_.append(raw_)
    n_samples = raw_._data.shape[1]
    with pytest.warns(UserWarning, match='did not converge'):
        ica.fit(raw, picks=picks[:5], decim=3)
    assert raw_._data.shape[1] == n_samples

    # test expl var
    ica = ICA(n_components=1.0, max_pca_components=4,
              n_pca_components=4, method=method, max_iter=1)
    with pytest.warns(UserWarning, match='did not converge'):
        ica.fit(raw, picks=None, decim=3)
    assert (ica.n_components_ == 4)
    ica_var = _ica_explained_variance(ica, raw, normalize=True)
#.........这里部分代码省略.........

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]

#.........这里部分代码省略.........

    picks_ = np.array([raw.ch_names.index(k) for k in ica.ch_names])
    if 'eeg' in ica:
        if 'ecg' in raw:
            picks_ = np.append(picks_,
                               pick_types(raw.info, meg=False, ecg=True)[0])
        else:
            logger.info('There is no ECG channel, trying to guess ECG from '
                        'magnetormeters')

    if artifact_stats is None:
        artifact_stats = dict()

    ecg_epochs = create_ecg_epochs(raw, tmin=ecg_tmin, tmax=ecg_tmax,
                                   keep_ecg=True, picks=picks_, reject=reject_)

    n_ecg_epochs_found = len(ecg_epochs.events)
    artifact_stats['ecg_n_events'] = n_ecg_epochs_found
    n_max_ecg_epochs = min(n_max_ecg_epochs, n_ecg_epochs_found)
    artifact_stats['ecg_n_used'] = n_max_ecg_epochs

    sel_ecg_epochs = np.arange(n_ecg_epochs_found)
    rng = np.random.RandomState(42)
    rng.shuffle(sel_ecg_epochs)
    ecg_ave = ecg_epochs.average()

    report.add_figs_to_section(ecg_ave.plot(), 'ECG-full', 'artifacts')
    ecg_epochs = ecg_epochs[sel_ecg_epochs[:n_max_ecg_epochs]]
    ecg_ave = ecg_epochs.average()
    report.add_figs_to_section(ecg_ave.plot(), 'ECG-used', 'artifacts')

    _put_artifact_range(artifact_stats, ecg_ave, kind='ecg')

    ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
    if len(ecg_inds) > 0:
        ecg_evoked = ecg_epochs.average()
        del ecg_epochs

        fig = ica.plot_scores(scores, exclude=ecg_inds, labels='ecg',
                              title='', show=show)

        report.add_figs_to_section(fig, 'scores ({})'.format(subject),
                                   section=comment + 'ECG',
                                   scale=img_scale)

        current_exclude = [e for e in ica.exclude]  # issue #2608 MNE
        fig = ica.plot_sources(raw, ecg_inds, exclude=ecg_inds,
                               title=title % ('components', 'ecg'), show=show)

        report.add_figs_to_section(fig, 'sources ({})'.format(subject),
                                   section=comment + 'ECG',
                                   scale=img_scale)
        ica.exclude = current_exclude

        fig = ica.plot_components(ecg_inds, ch_type=topo_ch_type,
                                  title='', colorbar=True, show=show)
        report.add_figs_to_section(fig, title % ('sources', 'ecg'),
                                   section=comment + 'ECG', scale=img_scale)
        ica.exclude = current_exclude

        ecg_inds = ecg_inds[:n_max_ecg]
        ica.exclude += ecg_inds
        fig = ica.plot_sources(ecg_evoked, exclude=ecg_inds, show=show)
        report.add_figs_to_section(fig, 'evoked sources ({})'.format(subject),
                                   section=comment + 'ECG',
                                   scale=img_scale)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]
def preprocess_ICA_fif_to_ts(fif_file, ECG_ch_name, EoG_ch_name, l_freq, h_freq, down_sfreq, variance, is_sensor_space, data_type):
    import os
    import numpy as np

    import mne
    from mne.io import Raw
    from mne.preprocessing import ICA, read_ica
    from mne.preprocessing import create_ecg_epochs, create_eog_epochs
    from mne.report import Report

    from nipype.utils.filemanip import split_filename as split_f

    report = Report()

    subj_path, basename, ext = split_f(fif_file)
    (data_path, sbj_name) = os.path.split(subj_path)
    print data_path

    # Read raw
    # If None the compensation in the data is not modified.
    # If set to n, e.g. 3, apply gradient compensation of grade n as for
    # CTF systems (compensation=3)
    raw = Raw(fif_file, preload=True)

    # select sensors
    select_sensors = mne.pick_types(raw.info, meg=True, ref_meg=False,
                                    exclude='bads')
    picks_meeg = mne.pick_types(raw.info, meg=True, eeg=True, exclude='bads')

    # save electrode locations
    sens_loc = [raw.info['chs'][i]['loc'][:3] for i in select_sensors]
    sens_loc = np.array(sens_loc)

    channel_coords_file = os.path.abspath("correct_channel_coords.txt")
    print '*** ' + channel_coords_file + '***'
    np.savetxt(channel_coords_file, sens_loc, fmt='%s')

    # save electrode names
    sens_names = np.array([raw.ch_names[pos] for pos in select_sensors],dtype = "str")

    # AP 21032016 
#    channel_names_file = os.path.join(data_path, "correct_channel_names.txt") 
    channel_names_file = os.path.abspath("correct_channel_names.txt")
    np.savetxt(channel_names_file,sens_names , fmt = '%s')
 
    ### filtering + downsampling
    raw.filter(l_freq=l_freq, h_freq=h_freq, picks=picks_meeg,
               method='iir', n_jobs=8)
#    raw.filter(l_freq = l_freq, h_freq = h_freq, picks = picks_meeg,
#               method='iir')
#    raw.resample(sfreq=down_sfreq, npad=0)

    ### 1) Fit ICA model using the FastICA algorithm
    # Other available choices are `infomax` or `extended-infomax`
    # We pass a float value between 0 and 1 to select n_components based on the
    # percentage of variance explained by the PCA components.
    ICA_title = 'Sources related to %s artifacts (red)'
    is_show = False # visualization
    reject = dict(mag=4e-12, grad=4000e-13)

    # check if we have an ICA, if yes, we load it
    ica_filename = os.path.join(subj_path,basename + "-ica.fif")  
    if os.path.exists(ica_filename) is False:
        ica = ICA(n_components=variance, method='fastica', max_iter=500) # , max_iter=500
        ica.fit(raw, picks=select_sensors, reject=reject) # decim = 3, 

        has_ICA = False
    else:
        has_ICA = True
        print ica_filename + '   exists!!!'
        ica = read_ica(ica_filename)
        ica.exclude = []

    # 2) identify bad components by analyzing latent sources.
    # generate ECG epochs use detection via phase statistics

    # if we just have exclude channels we jump these steps
#    if len(ica.exclude)==0:
    n_max_ecg = 3
    n_max_eog = 2

    # check if ECG_ch_name is in the raw channels
    if ECG_ch_name in raw.info['ch_names']:
        ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5,
                                       picks=select_sensors,
                                       ch_name=ECG_ch_name)
    # if not  a synthetic ECG channel is created from cross channel average
    else:
        ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5,
                                       picks=select_sensors)

    # ICA for ECG artifact
    # threshold=0.25 come default
    ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
    print scores
    print '\n len ecg_inds *** ' + str(len(ecg_inds)) + '***\n'
    if len(ecg_inds) > 0:
        ecg_evoked = ecg_epochs.average()

        fig1 = ica.plot_scores(scores, exclude=ecg_inds,
#.........这里部分代码省略.........

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]
def runICA(raw,saveRoot,name):

    saveRoot = saveRoot    
    icaList = [] 
    ica = []
    n_max_ecg = 3   # max number of ecg components 
#    n_max_eog_1 = 2 # max number of vert eog comps
#    n_max_eog_2 = 2 # max number of horiz eog comps          
    ecg_source_idx, ecg_scores, ecg_exclude = [], [], []
    eog_source_idx, eog_scores, eog_exclude = [], [], []
    #horiz = 1       # will later be modified to horiz = 0 if no horizontal EOG components are identified                   
    ica = ICA(n_components=0.90,n_pca_components=64,max_pca_components=100,noise_cov=None)
        
    ica.fit(raw)
    #*************
    eog_picks = mne.pick_types(raw.info, meg=False, eeg=False, stim=False, eog=True, ecg=False, emg=False)[0]
    ecg_picks = mne.pick_types(raw.info, meg=False, eeg=False, stim=False, ecg=True, eog=False, emg=False)[0]
    ica_picks = mne.pick_types(raw.info, meg=True, eeg=False, eog=False, ecg=False,
                   stim=False, exclude='bads')
    ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=ica_picks)
    ecg_evoked = ecg_epochs.average()
    eog_evoked = create_eog_epochs(raw, tmin=-.5, tmax=.5, picks=ica_picks).average()

    ecg_source_idx, ecg_scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
    eog_source_idx, eog_scores = ica.find_bads_eog(raw,ch_name=raw.ch_names[eog_picks].encode('ascii', 'ignore'))
       
    # defining a title-frame for later use
    title = 'Sources related to %s artifacts (red)'

    # extracting number of ica-components and plotting their topographies
    source_idx = range(0, ica.n_components_)
    ica_plot = ica.plot_components(source_idx, ch_type="mag")                                           

    # select ICA sources and reconstruct MEG signals, compute clean ERFs
    # Add detected artefact sources to exclusion list
    # We now add the eog artefacts to the ica.exclusion list
    if not ecg_source_idx:
        print("No ECG components above threshold were identified for subject " + name +
        " - selecting the component with the highest score under threshold")
        ecg_exclude = [np.absolute(ecg_scores).argmax()]
        ecg_source_idx=[np.absolute(ecg_scores).argmax()]
    elif ecg_source_idx:
        ecg_exclude += ecg_source_idx[:n_max_ecg]
    ica.exclude += ecg_exclude

    if not eog_source_idx:
        if np.absolute(eog_scores).any>0.3:
            eog_exclude=[np.absolute(eog_scores).argmax()]
            eog_source_idx=[np.absolute(eog_scores).argmax()]
            print("No EOG components above threshold were identified " + name +
            " - selecting the component with the highest score under threshold above 0.3")
        elif not np.absolute(eog_scores).any>0.3:
            eog_exclude=[]
            print("No EOG components above threshold were identified" + name)
    elif eog_source_idx:
         eog_exclude += eog_source_idx

    ica.exclude += eog_exclude

    print('########## saving')
    if len(eog_exclude) == 0:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg_none' + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 1:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 2:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 3:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '_' + map(str, ecg_exclude)[2] + '.pdf', format = 'pdf')
    elif len(eog_exclude) == 1:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg_none' + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 1:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 2:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 3:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '_' + map(str, ecg_exclude)[2] + '.pdf', format = 'pdf')
    elif len(eog_exclude) == 2:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] + '_' + map(str, eog_exclude)[1] +
            '-ecg_none' + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 1:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] + '_' + map(str, eog_exclude)[1] +
            '-ecg' + map(str, ecg_exclude)[0] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 2:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] + '_' + map(str, eog_exclude)[1] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 3:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] + '_' + map(str, eog_exclude)[1] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '_' + map(str, ecg_exclude)[2] + '.pdf', format = 'pdf')
    
    # plot the scores for the different components highlighting in red that/those related to ECG
    scores_plots_ecg=ica.plot_scores(ecg_scores, exclude=ecg_source_idx, title=title % 'ecg')
    scores_plots_ecg.savefig(saveRoot + name + '_ecg_scores.pdf', format = 'pdf')
    scores_plots_eog=ica.plot_scores(eog_scores, exclude=eog_source_idx, title=title % 'eog')
#.........这里部分代码省略.........

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]
def runICA(raw,saveRoot,name):

    saveRoot = saveRoot    
    icaList = [] 
    ica = []
    n_max_ecg = 3   # max number of ecg components 
#    n_max_eog_1 = 2 # max number of vert eog comps
#    n_max_eog_2 = 2 # max number of horiz eog comps          
    ecg_source_idx, ecg_scores, ecg_exclude = [], [], []
    eog_source_idx, eog_scores, eog_exclude = [], [], []
    #horiz = 1       # will later be modified to horiz = 0 if no horizontal EOG components are identified                   
    ica = ICA(n_components=0.90,n_pca_components=64,max_pca_components=100,noise_cov=None)
    
    fit_picks = mne.pick_types(raw.info, meg=True, eeg=True, eog=False, ecg=False,
                   stim=False, exclude='bads')    
    ica.fit(raw, picks=fit_picks)
    #ica.fit(raw)
    #*************
    eog_picks = mne.pick_types(raw.info, meg=False, eeg=False, stim=False, eog=True, ecg=False, emg=False)[0]
    ecg_picks = mne.pick_types(raw.info, meg=False, eeg=False, stim=False, ecg=True, eog=False, emg=False)[0]
    ica_picks = mne.pick_types(raw.info, meg=True, eeg=True, eog=False, ecg=False,
                   stim=False, exclude='bads')
    ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=ica_picks)
    ecg_evoked = ecg_epochs.average()
    eog_evoked = create_eog_epochs(raw, tmin=-.5, tmax=.5, picks=ica_picks).average()

    ecg_source_idx, ecg_scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
    eog_source_idx, eog_scores = ica.find_bads_eog(raw,ch_name=raw.ch_names[eog_picks].encode('ascii', 'ignore'))
    #eog_source_idx_2, eog_scores_2 = ica.find_bads_eog(raw,ch_name='EOG002')
    #if not eog_source_idx_2:
    #    horiz = 0
    
    #show_picks = np.abs(scores).argsort()[::-1][:5]
    #ica.plot_sources(raw, show_picks, exclude=ecg_inds, title=title % 'ecg')
    
        
    # defining a title-frame for later use
    title = 'Sources related to %s artifacts (red)'
  

    # extracting number of ica-components and plotting their topographies
    source_idx = range(0, ica.n_components_)
    #ica_plot = ica.plot_components(source_idx, ch_type="mag")
    ica_plot = ica.plot_components(source_idx)
                                          
    #ica_plot = ica.plot_components(source_idx)

    # select ICA sources and reconstruct MEG signals, compute clean ERFs
    # Add detected artefact sources to exclusion list
    # We now add the eog artefacts to the ica.exclusion list
    if not ecg_source_idx:
        print("No ECG components above threshold were identified for subject " + name +
        " - selecting the component with the highest score under threshold")
        ecg_exclude = [np.absolute(ecg_scores).argmax()]
        ecg_source_idx=[np.absolute(ecg_scores).argmax()]
    elif ecg_source_idx:
        ecg_exclude += ecg_source_idx[:n_max_ecg]
    ica.exclude += ecg_exclude

    if not eog_source_idx:
        if np.absolute(eog_scores).any>0.3:
            eog_exclude=[np.absolute(eog_scores).argmax()]
            eog_source_idx=[np.absolute(eog_scores).argmax()]
            print("No EOG components above threshold were identified " + name +
            " - selecting the component with the highest score under threshold above 0.3")
        elif not np.absolute(eog_scores).any>0.3:
            eog_exclude=[]
            print("No EOG components above threshold were identified" + name)
    elif eog_source_idx:
         eog_exclude += eog_source_idx

    ica.exclude += eog_exclude

    print('########## saving')
    if len(eog_exclude) == 0:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg_none' + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 1:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 2:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 3:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '_' + map(str, ecg_exclude)[2] + '.pdf', format = 'pdf')
    elif len(eog_exclude) == 1:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg_none' + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 1:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 2:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 3:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '_' + map(str, ecg_exclude)[2] + '.pdf', format = 'pdf')
    elif len(eog_exclude) == 2:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] + '_' + map(str, eog_exclude)[1] +
            '-ecg_none' + '.pdf', format = 'pdf')
#.........这里部分代码省略.........

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]
def compute_ica(subject):
    """Function will compute ICA on raw and apply the ICA.

    params:
    subject : str
        the subject id to be loaded
    """
    raw = Raw(save_folder + "%s_filtered_data_mc_raw_tsss.fif" % subject,
              preload=True)

    # ICA Part
    ica = ICA(n_components=0.95, method='fastica', max_iter=256)

    picks = mne.pick_types(raw.info, meg=True, eeg=True,
                           stim=False, exclude='bads')

    ica.fit(raw, picks=picks, decim=decim, reject=reject)

    # maximum number of components to reject
    n_max_ecg, n_max_eog = 3, 1

    ##########################################################################
    # 2) identify bad components by analyzing latent sources.
    title = 'Sources related to %s artifacts (red) for sub: %s'

    # generate ECG epochs use detection via phase statistics
    ecg_epochs = create_ecg_epochs(raw, ch_name="ECG002",
                                   tmin=-.5, tmax=.5, picks=picks)
    n_ecg_epochs_found = len(ecg_epochs.events)
    sel_ecg_epochs = np.arange(0, n_ecg_epochs_found, 10)
    ecg_epochs = ecg_epochs[sel_ecg_epochs]

    ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
    fig = ica.plot_scores(scores, exclude=ecg_inds,
                          title=title % ('ecg', subject))
    fig.savefig(save_folder + "pics/%s_ecg_scores.png" % subject)

    if ecg_inds:
        show_picks = np.abs(scores).argsort()[::-1][:5]

        fig = ica.plot_sources(raw, show_picks, exclude=ecg_inds,
                               title=title % ('ecg', subject), show=False)
        fig.savefig(save_folder + "pics/%s_ecg_sources.png" % subject)
        fig = ica.plot_components(ecg_inds, title=title % ('ecg', subject),
                                  colorbar=True)
        fig.savefig(save_folder + "pics/%s_ecg_component.png" % subject)

        ecg_inds = ecg_inds[:n_max_ecg]
        ica.exclude += ecg_inds

    # estimate average artifact
    ecg_evoked = ecg_epochs.average()
    del ecg_epochs

    # plot ECG sources + selection
    fig = ica.plot_sources(ecg_evoked, exclude=ecg_inds)
    fig.savefig(save_folder + "pics/%s_ecg_sources_ave.png" % subject)

    # plot ECG cleaning
    ica.plot_overlay(ecg_evoked, exclude=ecg_inds)
    fig.savefig(save_folder + "pics/%s_ecg_sources_clean_ave.png" % subject)

    # DETECT EOG BY CORRELATION
    # HORIZONTAL EOG
    eog_epochs = create_eog_epochs(raw, ch_name="EOG001")
    eog_inds, scores = ica.find_bads_eog(raw)
    fig = ica.plot_scores(scores, exclude=eog_inds,
                          title=title % ('eog', subject))
    fig.savefig(save_folder + "pics/%s_eog_scores.png" % subject)

    fig = ica.plot_components(eog_inds, title=title % ('eog', subject),
                              colorbar=True)
    fig.savefig(save_folder + "pics/%s_eog_component.png" % subject)

    eog_inds = eog_inds[:n_max_eog]
    ica.exclude += eog_inds

    del eog_epochs

    ##########################################################################
    # Apply the solution to Raw, Epochs or Evoked like this:
    raw_ica = ica.apply(raw, copy=False)
    ica.save(save_folder + "%s-ica.fif" % subject)  # save ICA componenets
    # Save raw with ICA removed
    raw_ica.save(save_folder + "%s_filtered_ica_mc_raw_tsss.fif" % subject,
                 overwrite=True)
    plt.close("all")

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]

#.........这里部分代码省略.........

        _raw1 = ica.apply(raw, exclude=[1])
        _raw2 = ica_read.apply(raw, exclude=[1])
        assert_array_almost_equal(_raw1[:, :][0], _raw2[:, :][0])

    os.remove(test_ica_fname)
    # check scrore funcs
    for name, func in score_funcs.items():
        if name in score_funcs_unsuited:
            continue
        scores = ica.score_sources(raw, target="EOG 061", score_func=func, start=0, stop=10)
        assert_true(ica.n_components_ == len(scores))

    # check univariate stats
    scores = ica.score_sources(raw, score_func=stats.skew)
    # check exception handling
    assert_raises(ValueError, ica.score_sources, raw, target=np.arange(1))

    params = []
    params += [(None, -1, slice(2), [0, 1])]  # varicance, kurtosis idx params
    params += [(None, "MEG 1531")]  # ECG / EOG channel params
    for idx, ch_name in product(*params):
        ica.detect_artifacts(
            raw,
            start_find=0,
            stop_find=50,
            ecg_ch=ch_name,
            eog_ch=ch_name,
            skew_criterion=idx,
            var_criterion=idx,
            kurt_criterion=idx,
        )

    idx, scores = ica.find_bads_ecg(raw, method="ctps")
    assert_equal(len(scores), ica.n_components_)
    idx, scores = ica.find_bads_ecg(raw, method="correlation")
    assert_equal(len(scores), ica.n_components_)
    idx, scores = ica.find_bads_ecg(epochs, method="ctps")
    assert_equal(len(scores), ica.n_components_)
    assert_raises(ValueError, ica.find_bads_ecg, epochs.average(), method="ctps")
    assert_raises(ValueError, ica.find_bads_ecg, raw, method="crazy-coupling")

    idx, scores = ica.find_bads_eog(raw)
    assert_equal(len(scores), ica.n_components_)
    raw.info["chs"][raw.ch_names.index("EOG 061") - 1]["kind"] = 202
    idx, scores = ica.find_bads_eog(raw)
    assert_true(isinstance(scores, list))
    assert_equal(len(scores[0]), ica.n_components_)

    # check score funcs
    for name, func in score_funcs.items():
        if name in score_funcs_unsuited:
            continue
        scores = ica.score_sources(epochs_eog, target="EOG 061", score_func=func)
        assert_true(ica.n_components_ == len(scores))

    # check univariate stats
    scores = ica.score_sources(epochs, score_func=stats.skew)

    # check exception handling
    assert_raises(ValueError, ica.score_sources, epochs, target=np.arange(1))

    # ecg functionality
    ecg_scores = ica.score_sources(raw, target="MEG 1531", score_func="pearsonr")

    with warnings.catch_warnings(record=True):  # filter attenuation warning

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_ecg [as 别名]
def preprocess_ICA_fif_to_ts(fif_file, ECG_ch_name, EoG_ch_name, l_freq, h_freq):
    # ------------------------ Import stuff ------------------------ #
    import os
    import mne
    import sys
    from mne.io import Raw
    from mne.preprocessing import ICA
    from mne.preprocessing import create_ecg_epochs, create_eog_epochs
    from nipype.utils.filemanip import split_filename as split_f
    from reportGen import generateReport
    import pickle

    subj_path, basename, ext = split_f(fif_file)
    # -------------------- Delete later ------------------- #
    subj_name = subj_path[-5:]
    results_dir = subj_path[:-6]
    # results_dir += '2016'
    subj_path = results_dir + '/' + subj_name
    if not os.path.exists(subj_path):
        try:
            os.makedirs(subj_path)
        except OSError:
            sys.stderr.write('ica_preproc: problem creating directory: ' + subj_path)
    ########################################################
    # Read raw
    #   If None the compensation in the data is not modified. If set to n, e.g. 3, apply
    #   gradient compensation of grade n as for CTF systems (compensation=3)
    print(fif_file)
    print(EoG_ch_name)
    #  ----------------------------- end Import stuff ----------------- #
    # EoG_ch_name = "EOG061, EOG062"

    # ------------- Load raw ------------- #
    raw = Raw(fif_file, preload=True)
    # select sensors
    select_sensors = mne.pick_types(raw.info, meg=True, ref_meg=False, exclude='bads')
    picks_meeg = mne.pick_types(raw.info, meg=True, eeg=True, exclude='bads')

    # filtering
    raw.filter(l_freq=l_freq, h_freq=h_freq, picks=picks_meeg, method='iir', n_jobs=1)

    # if ECG_ch_name == 'EMG063':
    if ECG_ch_name in raw.info['ch_names']:
        raw.set_channel_types({ECG_ch_name: 'ecg'})  # Without this files with ECG_ch_name = 'EMG063' fail
        # ECG_ch_name = 'ECG063'
    if EoG_ch_name == 'EMG065,EMG066,EMG067,EMG068':   # Because ica.find_bads_eog... can process max 2 EoG channels
        EoG_ch_name = 'EMG065,EMG067'                 # it won't fail if I specify 4 channels, but it'll use only first
                                                      # EMG065 and EMG066 are for vertical eye movements and
                                                      # EMG067 and EMG068 are for horizontal

    # print rnk
    rnk = 'N/A'
    # 1) Fit ICA model using the FastICA algorithm
    # Other available choices are `infomax` or `extended-infomax`
    # We pass a float value between 0 and 1 to select n_components based on the
    # percentage of variance explained by the PCA components.
    reject = dict(mag=10e-12, grad=10000e-13)
    flat = dict(mag=0.1e-12, grad=1e-13)
    # check if we have an ICA, if yes, we load it
    ica_filename = os.path.join(subj_path, basename + "-ica.fif")
    raw_ica_filename = os.path.join(subj_path, basename + "_ica_raw.fif")
    info_filename = os.path.join(subj_path, basename + "_info.pickle")
    # if os.path.exists(ica_filename) == False:
    ica = ICA(n_components=0.99, method='fastica')  # , max_iter=500
    ica.fit(raw, picks=select_sensors, reject=reject, flat=flat)  # decim = 3,
    # has_ICA = False
    # else:
    #     has_ICA = True
    #     ica = read_ica(ica_filename)
    #     ica.exclude = []
    # ica.labels_ = dict() # to avoid bug; Otherwise it'll throw an exception

    ica_sources_filename = subj_path + '/' + basename + '_ica_timecourse.fif'

    # if not os.path.isfile(ica_sources_filename):
    icaSrc = ica.get_sources(raw, add_channels=None, start=None, stop=None)
    icaSrc.save(ica_sources_filename, picks=None, tmin=0, tmax=None, buffer_size_sec=10,
                drop_small_buffer=False, proj=False, fmt='single', overwrite=True, split_size='2GB', verbose=None)
    icaSrc = None
    # if has_ICA == False:
    # ica.save(ica_filename)
    # return
    # 2) identify bad components by analyzing latent sources.
    # generate ECG epochs use detection via phase statistics

    # check if ECG_ch_name is in the raw channels
    # import ipdb; ipdb.set_trace()
    if ECG_ch_name in raw.info['ch_names']:
        ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=select_sensors, ch_name=ECG_ch_name)
    # if not  a synthetic ECG channel is created from cross channel average
    else:
        ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=select_sensors)
    # ICA for ECG artifact
    # threshold=0.25 come defualt
    ecg_inds, ecg_scores = ica.find_bads_ecg(ecg_epochs, method='ctps', threshold=0.25)
    # if len(ecg_inds) > 0:
    ecg_evoked = ecg_epochs.average()
    ecg_epochs = None    # ecg_epochs use too much memory
    n_max_ecg = 3
    ecg_inds = ecg_inds[:n_max_ecg]
#.........这里部分代码省略.........