Python ICA.plot_components方法代码示例

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [as 别名]
def run_ica(method):
    ica = ICA(n_components=20, method=method, random_state=0)
    t0 = time()
    ica.fit(raw, picks=picks, reject=reject)
    fit_time = time() - t0
    title = ('ICA decomposition using %s (took %.1fs)' % (method, fit_time))
    ica.plot_components(title=title)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [as 别名]
def test_plot_ica_components():
    """Test plotting of ICA solutions."""
    import matplotlib.pyplot as plt
    res = 8
    fast_test = {"res": res, "contours": 0, "sensors": False}
    raw = _get_raw()
    ica = ICA(noise_cov=read_cov(cov_fname), n_components=2,
              max_pca_components=3, n_pca_components=3)
    ica_picks = _get_picks(raw)
    with warnings.catch_warnings(record=True):
        ica.fit(raw, picks=ica_picks)
    warnings.simplefilter('always', UserWarning)
    with warnings.catch_warnings(record=True):
        for components in [0, [0], [0, 1], [0, 1] * 2, None]:
            ica.plot_components(components, image_interp='bilinear',
                                colorbar=True, **fast_test)
        plt.close('all')

        # test interactive mode (passing 'inst' arg)
        ica.plot_components([0, 1], image_interp='bilinear', inst=raw, res=16)
        fig = plt.gcf()

        # test title click
        # ----------------
        lbl = fig.axes[1].get_label()
        ica_idx = int(lbl[-3:])
        titles = [ax.title for ax in fig.axes]
        title_pos_midpoint = (titles[1].get_window_extent().extents
                              .reshape((2, 2)).mean(axis=0))
        # first click adds to exclude
        _fake_click(fig, fig.axes[1], title_pos_midpoint, xform='pix')
        assert ica_idx in ica.exclude
        # clicking again removes from exclude
        _fake_click(fig, fig.axes[1], title_pos_midpoint, xform='pix')
        assert ica_idx not in ica.exclude

        # test topo click
        # ---------------
        _fake_click(fig, fig.axes[1], (0., 0.), xform='data')

        c_fig = plt.gcf()
        labels = [ax.get_label() for ax in c_fig.axes]

        for l in ['topomap', 'image', 'erp', 'spectrum', 'variance']:
            assert_true(l in labels)

        topomap_ax = c_fig.axes[labels.index('topomap')]
        title = topomap_ax.get_title()
        assert_true(lbl == title)

    ica.info = None
    assert_raises(ValueError, ica.plot_components, 1)
    assert_raises(RuntimeError, ica.plot_components, 1, ch_type='mag')
    plt.close('all')

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [as 别名]
def test_plot_ica_components():
    """Test plotting of ICA solutions
    """
    raw = _get_raw()
    ica = ICA(noise_cov=read_cov(cov_fname), n_components=2, max_pca_components=3, n_pca_components=3)
    ica_picks = pick_types(raw.info, meg=True, eeg=False, stim=False, ecg=False, eog=False, exclude="bads")
    ica.fit(raw, picks=ica_picks)
    warnings.simplefilter("always", UserWarning)
    with warnings.catch_warnings(record=True):
        for components in [0, [0], [0, 1], [0, 1] * 7, None]:
            ica.plot_components(components)
    ica.info = None
    assert_raises(RuntimeError, ica.plot_components, 1)
    plt.close("all")

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [as 别名]
def test_plot_ica_components():
    """Test plotting of ICA solutions
    """
    raw = _get_raw()
    ica = ICA(noise_cov=read_cov(cov_fname), n_components=2,
              max_pca_components=3, n_pca_components=3)
    ica_picks = _get_picks(raw)
    ica.fit(raw, picks=ica_picks)
    warnings.simplefilter('always', UserWarning)
    with warnings.catch_warnings(record=True):
        for components in [0, [0], [0, 1], [0, 1] * 2, None]:
            ica.plot_components(components, image_interp='bilinear', res=16)
    ica.info = None
    assert_raises(RuntimeError, ica.plot_components, 1)
    plt.close('all')

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [as 别名]
def test_plot_ica_components():
    """Test plotting of ICA solutions."""
    import matplotlib.pyplot as plt
    res = 8
    fast_test = {"res": res, "contours": 0, "sensors": False}
    raw = _get_raw()
    ica = ICA(noise_cov=read_cov(cov_fname), n_components=2,
              max_pca_components=3, n_pca_components=3)
    ica_picks = _get_picks(raw)
    with warnings.catch_warnings(record=True):
        ica.fit(raw, picks=ica_picks)
    warnings.simplefilter('always', UserWarning)
    with warnings.catch_warnings(record=True):
        for components in [0, [0], [0, 1], [0, 1] * 2, None]:
            ica.plot_components(components, image_interp='bilinear',
                                colorbar=True, **fast_test)

        # test interactive mode (passing 'inst' arg)
        plt.close('all')
        ica.plot_components([0, 1], image_interp='bilinear', inst=raw, res=16)

        fig = plt.gcf()
        ax = [a for a in fig.get_children() if isinstance(a, plt.Axes)]
        lbl = ax[1].get_label()
        _fake_click(fig, ax[1], (0., 0.), xform='data')

        c_fig = plt.gcf()
        ax = [a for a in c_fig.get_children() if isinstance(a, plt.Axes)]
        labels = [a.get_label() for a in ax]

        for l in ['topomap', 'image', 'erp', 'spectrum', 'variance']:
            assert_true(l in labels)

        topomap_ax = ax[labels.index('topomap')]
        title = topomap_ax.get_title()
        assert_true(lbl == title)

    ica.info = None
    assert_raises(ValueError, ica.plot_components, 1)
    assert_raises(RuntimeError, ica.plot_components, 1, ch_type='mag')
    plt.close('all')

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [as 别名]
def test_plot_ica_components():
    """Test plotting of ICA solutions."""
    import matplotlib.pyplot as plt

    raw = _get_raw()
    ica = ICA(noise_cov=read_cov(cov_fname), n_components=2, max_pca_components=3, n_pca_components=3)
    ica_picks = _get_picks(raw)
    with warnings.catch_warnings(record=True):
        ica.fit(raw, picks=ica_picks)
    warnings.simplefilter("always", UserWarning)
    with warnings.catch_warnings(record=True):
        for components in [0, [0], [0, 1], [0, 1] * 2, None]:
            ica.plot_components(components, image_interp="bilinear", res=16, colorbar=True)

        # test interactive mode (passing 'inst' arg)
        plt.close("all")
        ica.plot_components([0, 1], image_interp="bilinear", res=16, inst=raw)

        fig = plt.gcf()
        ax = [a for a in fig.get_children() if isinstance(a, plt.Axes)]
        lbl = ax[1].get_label()
        _fake_click(fig, ax[1], (0.0, 0.0), xform="data")

        c_fig = plt.gcf()
        ax = [a for a in c_fig.get_children() if isinstance(a, plt.Axes)]
        labels = [a.get_label() for a in ax]

        for l in ["topomap", "image", "erp", "spectrum", "variance"]:
            assert_true(l in labels)

        topomap_ax = ax[labels.index("topomap")]
        title = topomap_ax.get_title()
        assert_true(lbl == title)

    ica.info = None
    assert_raises(ValueError, ica.plot_components, 1)
    assert_raises(RuntimeError, ica.plot_components, 1, ch_type="mag")
    plt.close("all")

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [as 别名]
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)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [as 别名]
    # 2) identify bad components by analyzing latent sources.

    title = "Sources related to %s artifacts (red)"

    # generate ECG epochs use detection via phase statistics

    ecg_epochs = create_ecg_epochs(raw, tmin=-0.5, tmax=0.5, picks=picks)

    ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method="ctps")
    ica.plot_scores(scores, exclude=ecg_inds, title=title % "ecg")

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

        ica.plot_sources(raw, show_picks, exclude=ecg_inds, title=title % "ecg")
        ica.plot_components(ecg_inds, title=title % "ecg", colorbar=True)

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

    # detect EOG by correlation

    eog_inds, scores = ica.find_bads_eog(raw)
    ica.plot_scores(scores, exclude=eog_inds, title=title % "eog")

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

        #        ica.plot_sources(raw, show_picks, exclude=eog_inds,
        #                         title="Sources related to EOG artifacts (red)")
        ica.plot_components(eog_inds, title="Sources related to EOG artifacts", colorbar=True)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [as 别名]
tmin, tmax = -0.2, 0.6
baseline = None  # no baseline as high-pass is applied
reject = dict(mag=5e-12)

epochs = mne.Epochs(raw, events, event_ids, tmin, tmax,  picks=picks,
                    baseline=baseline, preload=True, reject=reject)

# Fit ICA, find and remove major artifacts
ica = ICA(n_components=0.95, random_state=0).fit(raw, decim=1, reject=reject)

# compute correlation scores, get bad indices sorted by score
eog_epochs = create_eog_epochs(raw, ch_name='MRT31-2908', reject=reject)
eog_inds, eog_scores = ica.find_bads_eog(eog_epochs, ch_name='MRT31-2908')
ica.plot_scores(eog_scores, eog_inds)  # see scores the selection is based on
ica.plot_components(eog_inds)  # view topographic sensitivity of components
ica.exclude += eog_inds[:1]  # we saw the 2nd ECG component looked too dipolar
ica.plot_overlay(eog_epochs.average())  # inspect artifact removal
ica.apply(epochs)  # clean data, default in place

evoked = [epochs[k].average() for k in event_ids]

contrast = combine_evoked(evoked, weights=[-1, 1])  # Faces - scrambled

evoked.append(contrast)

for e in evoked:
    e.plot(ylim=dict(mag=[-400, 400]))

plt.show()

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [as 别名]
ica.fit(raw, picks=picks, decim=3, reject=dict(mag=4e-12, grad=4000e-13))

# create EOG epochs to improve detection by correlation
picks = mne.pick_types(raw.info, meg=True, eog=True)
eog_epochs = create_eog_epochs(raw, picks=picks)

eog_inds, scores = ica.find_bads_eog(eog_epochs)  # inds sorted!

ica.plot_scores(scores, exclude=eog_inds)  # inspect metrics used

show_picks = np.abs(scores).argsort()[::-1][:5]  # indices of top five scores

# detected artifacts drawn in red (via exclude)
ica.plot_sources(raw, show_picks, exclude=eog_inds, start=0., stop=3.0)
ica.plot_components(eog_inds, colorbar=False)  # show component sensitivites

ica.exclude += eog_inds[:1]  # mark first for exclusion

###############################################################################
# 3) check detection and visualize artifact rejection

# estimate average artifact
eog_evoked = eog_epochs.average()
ica.plot_sources(eog_evoked)  # latent EOG sources + selction
ica.plot_overlay(eog_evoked)  # overlay raw and clean EOG artifacts

# check the amplitudes do not change
ica.plot_overlay(raw)  # EOG artifacts remain

###############################################################################

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [as 别名]
def ica_method(raw, picks, plot='n', save ='n'):
    
    
    ###############################################################################
    # 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 = ICA(n_components=0.95, method='fastica')
    
    picks = mne.pick_types(raw.info, meg=False, eeg=True, eog=False,
                           stim=False, exclude='bads')
    
    ica.fit(raw, picks=picks, decim=3)
    
    # maximum number of components to reject
    n_max_eog =  1  # here we don't expect horizontal EOG components
    
    ###############################################################################
    # 2) identify bad components by analyzing latent sources.
    
    
    # detect EOG by correlation
    
    eog_inds, scores = ica.find_bads_eog(raw, threshold=2.5)    
    show_picks = np.abs(scores).argsort()[::-1][:5]
    eog_inds = eog_inds[:n_max_eog]
    ica.exclude += eog_inds
    
    ###############################################################################
    # 3) Assess component selection and unmixing quality
    eog_evoked = create_eog_epochs(raw, tmin=-.5, tmax=.5, picks=picks).average()
    
    if plot=='y':
        
        title = 'Sources related to %s artifacts (red)'
        ica.plot_scores(scores, exclude=eog_inds, title=title % 'eog', labels='eog')
        if save=='y':
            pylab.savefig('2.png')

        ica.plot_sources(raw, show_picks, exclude=eog_inds, title=title % 'eog')
        if save=='y':
            pylab.savefig('3.png')

        ica.plot_components(eog_inds, title=title % 'eog', colorbar=True)
        if save=='y':

            pylab.savefig('4.png')

        ica.plot_overlay(raw)  # EOG artifacts remain
        if save=='y':

            pylab.savefig('5.png')

        ica.plot_sources(eog_evoked, exclude=eog_inds)  # plot EOG sources + selection
        if save=='y':

            pylab.savefig('6.png')

        ica.plot_overlay(eog_evoked, exclude=eog_inds)  # plot EOG cleaning
        if save=='y':
            pylab.savefig('7.png')


    ica.apply(raw, exclude=eog_inds)    
    eeg_only_after=raw.pick_types(meg=False, eeg=True)    
    

    return eeg_only_after
    

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [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 plot_components [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 plot_components [as 别名]
 # %%
# eog_inds = [10]
# ica.exclude += eog_inds

fig = ica.plot_scores(scores, exclude=eog_inds,
                      title=title % ('eog', subject))
fig.savefig(save_folder + "pics/%s_%s_eog_scores.png" % (subject,
                                                         condition))
                                                         

fig = ica.plot_sources(eog_average, exclude=None)
fig.savefig(save_folder + "pics/%s_%s_eog_source.png" % (subject,
                                                         condition))

fig = ica.plot_components(ica.exclude, title=title % ('eog', subject),
                          colorbar=True)
fig.savefig(save_folder + "pics/%s_%s_eog_component.png" % (subject,
                                                            condition))
fig = ica.plot_overlay(eog_average, exclude=None, show=False)                                                                
fig.savefig(save_folder + "pics/%s_%s_eog_excluded.png" % (subject,
                                                            condition))


# del eog_epochs, eog_average

##########################################################################
# Apply the solution to Raw, Epochs or Evoked like this:
raw_ica = ica.apply(raw)
ica.save(save_folder + "%s_%s-ica.fif" % (subject, condition))  # save ICA
# componenets
# Save raw with ICA removed

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import plot_components [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")