Python ICA.apply方法代码示例

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import apply [as 别名]
def test_ica_ctf():
    """Test run ICA computation on ctf data with/without compensation."""
    method = 'fastica'
    raw = read_raw_ctf(ctf_fname, preload=True)
    events = make_fixed_length_events(raw, 99999)
    for comp in [0, 1]:
        raw.apply_gradient_compensation(comp)
        epochs = Epochs(raw, events, None, -0.2, 0.2, preload=True)
        evoked = epochs.average()

        # test fit
        for inst in [raw, epochs]:
            ica = ICA(n_components=2, random_state=0, max_iter=2,
                      method=method)
            with pytest.warns(UserWarning, match='did not converge'):
                ica.fit(inst)

        # test apply and get_sources
        for inst in [raw, epochs, evoked]:
            ica.apply(inst)
            ica.get_sources(inst)

    # test mixed compensation case
    raw.apply_gradient_compensation(0)
    ica = ICA(n_components=2, random_state=0, max_iter=2, method=method)
    with pytest.warns(UserWarning, match='did not converge'):
        ica.fit(raw)
    raw.apply_gradient_compensation(1)
    epochs = Epochs(raw, events, None, -0.2, 0.2, preload=True)
    evoked = epochs.average()
    for inst in [raw, epochs, evoked]:
        with pytest.raises(RuntimeError, match='Compensation grade of ICA'):
            ica.apply(inst)
        with pytest.raises(RuntimeError, match='Compensation grade of ICA'):
            ica.get_sources(inst)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import apply [as 别名]
def test_ica_eeg():
    """Test ICA on EEG."""
    method = 'fastica'
    raw_fif = read_raw_fif(fif_fname, preload=True)
    with pytest.warns(RuntimeWarning, match='events'):
        raw_eeglab = read_raw_eeglab(input_fname=eeglab_fname,
                                     montage=eeglab_montage, preload=True)
    for raw in [raw_fif, raw_eeglab]:
        events = make_fixed_length_events(raw, 99999, start=0, stop=0.3,
                                          duration=0.1)
        picks_meg = pick_types(raw.info, meg=True, eeg=False)[:2]
        picks_eeg = pick_types(raw.info, meg=False, eeg=True)[:2]
        picks_all = []
        picks_all.extend(picks_meg)
        picks_all.extend(picks_eeg)
        epochs = Epochs(raw, events, None, -0.1, 0.1, preload=True)
        evoked = epochs.average()

        for picks in [picks_meg, picks_eeg, picks_all]:
            if len(picks) == 0:
                continue
            # test fit
            for inst in [raw, epochs]:
                ica = ICA(n_components=2, random_state=0, max_iter=2,
                          method=method)
                with pytest.warns(None):
                    ica.fit(inst, picks=picks)

            # test apply and get_sources
            for inst in [raw, epochs, evoked]:
                ica.apply(inst)
                ica.get_sources(inst)

    with pytest.warns(RuntimeWarning, match='MISC channel'):
        raw = read_raw_ctf(ctf_fname2,  preload=True)
    events = make_fixed_length_events(raw, 99999, start=0, stop=0.2,
                                      duration=0.1)
    picks_meg = pick_types(raw.info, meg=True, eeg=False)[:2]
    picks_eeg = pick_types(raw.info, meg=False, eeg=True)[:2]
    picks_all = picks_meg + picks_eeg
    for comp in [0, 1]:
        raw.apply_gradient_compensation(comp)
        epochs = Epochs(raw, events, None, -0.1, 0.1, preload=True)
        evoked = epochs.average()

        for picks in [picks_meg, picks_eeg, picks_all]:
            if len(picks) == 0:
                continue
            # test fit
            for inst in [raw, epochs]:
                ica = ICA(n_components=2, random_state=0, max_iter=2,
                          method=method)
                with pytest.warns(None):
                    ica.fit(inst)

            # test apply and get_sources
            for inst in [raw, epochs, evoked]:
                ica.apply(inst)
                ica.get_sources(inst)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import apply [as 别名]
def test_ica_full_data_recovery(method):
    """Test recovery of full data when no source is rejected."""
    # Most basic recovery
    _skip_check_picard(method)
    raw = read_raw_fif(raw_fname).crop(0.5, stop).load_data()
    events = read_events(event_name)
    picks = pick_types(raw.info, meg=True, stim=False, ecg=False,
                       eog=False, exclude='bads')[:10]
    with pytest.warns(RuntimeWarning, match='projection'):
        epochs = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks,
                        baseline=(None, 0), preload=True)
    evoked = epochs.average()
    n_channels = 5
    data = raw._data[:n_channels].copy()
    data_epochs = epochs.get_data()
    data_evoked = evoked.data
    raw.set_annotations(Annotations([0.5], [0.5], ['BAD']))
    methods = [method]
    for method in methods:
        stuff = [(2, n_channels, True), (2, n_channels // 2, False)]
        for n_components, n_pca_components, ok in stuff:
            ica = ICA(n_components=n_components, random_state=0,
                      max_pca_components=n_pca_components,
                      n_pca_components=n_pca_components,
                      method=method, max_iter=1)
            with pytest.warns(UserWarning, match=None):  # sometimes warns
                ica.fit(raw, picks=list(range(n_channels)))
            raw2 = ica.apply(raw.copy(), exclude=[])
            if ok:
                assert_allclose(data[:n_channels], raw2._data[:n_channels],
                                rtol=1e-10, atol=1e-15)
            else:
                diff = np.abs(data[:n_channels] - raw2._data[:n_channels])
                assert (np.max(diff) > 1e-14)

            ica = ICA(n_components=n_components, method=method,
                      max_pca_components=n_pca_components,
                      n_pca_components=n_pca_components, random_state=0)
            with pytest.warns(None):  # sometimes warns
                ica.fit(epochs, picks=list(range(n_channels)))
            epochs2 = ica.apply(epochs.copy(), exclude=[])
            data2 = epochs2.get_data()[:, :n_channels]
            if ok:
                assert_allclose(data_epochs[:, :n_channels], data2,
                                rtol=1e-10, atol=1e-15)
            else:
                diff = np.abs(data_epochs[:, :n_channels] - data2)
                assert (np.max(diff) > 1e-14)

            evoked2 = ica.apply(evoked.copy(), exclude=[])
            data2 = evoked2.data[:n_channels]
            if ok:
                assert_allclose(data_evoked[:n_channels], data2,
                                rtol=1e-10, atol=1e-15)
            else:
                diff = np.abs(evoked.data[:n_channels] - data2)
                assert (np.max(diff) > 1e-14)
    pytest.raises(ValueError, ICA, method='pizza-decomposision')

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import apply [as 别名]
def test_ica_full_data_recovery():
    """Test recovery of full data when no source is rejected"""
    # Most basic recovery
    raw = Raw(raw_fname).crop(0.5, stop, False)
    raw.load_data()
    events = read_events(event_name)
    picks = pick_types(raw.info, meg=True, stim=False, ecg=False,
                       eog=False, exclude='bads')[:10]
    with warnings.catch_warnings(record=True):  # bad proj
        epochs = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks,
                        baseline=(None, 0), preload=True)
    evoked = epochs.average()
    n_channels = 5
    data = raw._data[:n_channels].copy()
    data_epochs = epochs.get_data()
    data_evoked = evoked.data
    for method in ['fastica']:
        stuff = [(2, n_channels, True), (2, n_channels // 2, False)]
        for n_components, n_pca_components, ok in stuff:
            ica = ICA(n_components=n_components,
                      max_pca_components=n_pca_components,
                      n_pca_components=n_pca_components,
                      method=method, max_iter=1)
            with warnings.catch_warnings(record=True):
                ica.fit(raw, picks=list(range(n_channels)))
            raw2 = ica.apply(raw, exclude=[], copy=True)
            if ok:
                assert_allclose(data[:n_channels], raw2._data[:n_channels],
                                rtol=1e-10, atol=1e-15)
            else:
                diff = np.abs(data[:n_channels] - raw2._data[:n_channels])
                assert_true(np.max(diff) > 1e-14)

            ica = ICA(n_components=n_components,
                      max_pca_components=n_pca_components,
                      n_pca_components=n_pca_components)
            with warnings.catch_warnings(record=True):
                ica.fit(epochs, picks=list(range(n_channels)))
            epochs2 = ica.apply(epochs, exclude=[], copy=True)
            data2 = epochs2.get_data()[:, :n_channels]
            if ok:
                assert_allclose(data_epochs[:, :n_channels], data2,
                                rtol=1e-10, atol=1e-15)
            else:
                diff = np.abs(data_epochs[:, :n_channels] - data2)
                assert_true(np.max(diff) > 1e-14)

            evoked2 = ica.apply(evoked, exclude=[], copy=True)
            data2 = evoked2.data[:n_channels]
            if ok:
                assert_allclose(data_evoked[:n_channels], data2,
                                rtol=1e-10, atol=1e-15)
            else:
                diff = np.abs(evoked.data[:n_channels] - data2)
                assert_true(np.max(diff) > 1e-14)
    assert_raises(ValueError, ICA, method='pizza-decomposision')

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import apply [as 别名]
def test_ica_rank_reduction():
    """Test recovery of full data when no source is rejected"""
    # Most basic recovery
    raw = Raw(raw_fname).crop(0.5, stop, False)
    raw.load_data()
    picks = pick_types(raw.info, meg=True, stim=False, ecg=False,
                       eog=False, exclude='bads')[:10]
    n_components = 5
    max_pca_components = len(picks)
    for n_pca_components in [6, 10]:
        with warnings.catch_warnings(record=True):  # non-convergence
            warnings.simplefilter('always')
            ica = ICA(n_components=n_components,
                      max_pca_components=max_pca_components,
                      n_pca_components=n_pca_components,
                      method='fastica', max_iter=1).fit(raw, picks=picks)

        rank_before = raw.estimate_rank(picks=picks)
        assert_equal(rank_before, len(picks))
        raw_clean = ica.apply(raw, copy=True)
        rank_after = raw_clean.estimate_rank(picks=picks)
        # interaction between ICA rejection and PCA components difficult
        # to preduct. Rank_after often seems to be 1 higher then
        # n_pca_components
        assert_true(n_components < n_pca_components <= rank_after <=
                    rank_before)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import apply [as 别名]
def test_ica_rank_reduction(method):
    """Test recovery ICA rank reduction."""
    _skip_check_picard(method)
    # Most basic recovery
    raw = read_raw_fif(raw_fname).crop(0.5, stop).load_data()
    picks = pick_types(raw.info, meg=True, stim=False, ecg=False,
                       eog=False, exclude='bads')[:10]
    n_components = 5
    max_pca_components = len(picks)
    for n_pca_components in [6, 10]:
        with pytest.warns(UserWarning, match='did not converge'):
            ica = ICA(n_components=n_components,
                      max_pca_components=max_pca_components,
                      n_pca_components=n_pca_components,
                      method=method, max_iter=1).fit(raw, picks=picks)

        rank_before = _compute_rank_int(raw.copy().pick(picks), proj=False)
        assert_equal(rank_before, len(picks))
        raw_clean = ica.apply(raw.copy())
        rank_after = _compute_rank_int(raw_clean.copy().pick(picks),
                                       proj=False)
        # interaction between ICA rejection and PCA components difficult
        # to preduct. Rank_after often seems to be 1 higher then
        # n_pca_components
        assert (n_components < n_pca_components <= rank_after <=
                rank_before)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import apply [as 别名]
def test_ica_twice(method):
    """Test running ICA twice."""
    _skip_check_picard(method)
    raw = read_raw_fif(raw_fname).crop(1.5, stop).load_data()
    picks = pick_types(raw.info, meg='grad', exclude='bads')
    n_components = 0.9
    max_pca_components = None
    n_pca_components = 1.1
    ica1 = ICA(n_components=n_components, method=method,
               max_pca_components=max_pca_components,
               n_pca_components=n_pca_components, random_state=0)

    ica1.fit(raw, picks=picks, decim=3)
    raw_new = ica1.apply(raw, n_pca_components=n_pca_components)
    ica2 = ICA(n_components=n_components, method=method,
               max_pca_components=max_pca_components,
               n_pca_components=1.0, random_state=0)
    ica2.fit(raw_new, picks=picks, decim=3)
    assert_equal(ica1.n_components_, ica2.n_components_)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import apply [as 别名]
def test_ica_twice():
    """Test running ICA twice."""
    raw = read_raw_fif(raw_fname).crop(1.5, stop).load_data()
    picks = pick_types(raw.info, meg='grad', exclude='bads')
    n_components = 0.9
    max_pca_components = None
    n_pca_components = 1.1
    with warnings.catch_warnings(record=True):
        ica1 = ICA(n_components=n_components,
                   max_pca_components=max_pca_components,
                   n_pca_components=n_pca_components, random_state=0)

        ica1.fit(raw, picks=picks, decim=3)
        raw_new = ica1.apply(raw, n_pca_components=n_pca_components)
        ica2 = ICA(n_components=n_components,
                   max_pca_components=max_pca_components,
                   n_pca_components=1.0, random_state=0)
        ica2.fit(raw_new, picks=picks, decim=3)
        assert_equal(ica1.n_components_, ica2.n_components_)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import apply [as 别名]
def test_ica_twice():
    """Test running ICA twice"""
    raw = io.Raw(raw_fname, preload=True).crop(0, stop, False).crop(1.5)
    picks = pick_types(raw.info, meg="grad", exclude="bads")
    n_components = 0.9
    max_pca_components = None
    n_pca_components = 1.1
    with warnings.catch_warnings(record=True):
        ica1 = ICA(
            n_components=n_components,
            max_pca_components=max_pca_components,
            n_pca_components=n_pca_components,
            random_state=0,
        )

        ica1.fit(raw, picks=picks, decim=3)
        raw_new = ica1.apply(raw, n_pca_components=n_pca_components)
        ica2 = ICA(
            n_components=n_components, max_pca_components=max_pca_components, n_pca_components=1.0, random_state=0
        )
        ica2.fit(raw_new, picks=picks, decim=3)
        assert_equal(ica1.n_components_, ica2.n_components_)

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

###############################################################################
# Note that nothing is yet removed from the raw data. To remove the effects of
# the rejected components,
# :meth:`the apply method <mne.preprocessing.ICA.apply>` must be called.
# Here we apply it on the copy of the first ten seconds, so that the rest of
# this tutorial still works as intended.
raw_copy = raw.copy().crop(0, 10)
ica.apply(raw_copy)
raw_copy.plot()  # check the result

###############################################################################
# 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

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

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

    ###########################################################################
    # 3) Assess component selection and unmixing quality

    # estimate average artifact
    ecg_evoked = ecg_epochs.average()
    # plot ECG sources + selection
    ica.plot_sources(ecg_evoked, exclude=ecg_inds)
    # plot ECG cleaning
    ica.plot_overlay(ecg_evoked, exclude=ecg_inds)

    eog_evoked = create_eog_epochs(raw, tmin=-0.5, tmax=0.5, picks=picks).average()
    # plot EOG sources + selection
    # ica.plot_sources(eog_evoked, exclude=eog_inds)
    ica.plot_overlay(eog_evoked, exclude=eog_inds)  # plot EOG cleaning

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

    ##########################################################################
    # Apply the solution to Raw, Epochs or Evoked like this:
    raw_ica = ica.apply(raw, copy=False)
    raw_ica.save(data_path + "tone_task-%s-tsss-mc-autobad-ica_raw.fif" % (condition), overwrite=True)

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

    Parameters
    ----------
    subject : string
        the subject id to be loaded
    """
    raw = mne.io.Raw(data_folder + "%s_bp-raw.fif" % subject, preload=True)
    raw.set_montage = montage
    raw.apply_proj()
    # raw.resample(512, n_jobs=2)

    # ICA Part
    ica = ICA(n_components=None,
              max_pca_components=40,
              method='fastica',
              max_iter=256)

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

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

    # maximum number of components to reject
    n_max_eog = 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="Ext4",
    #                                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(data_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(data_folder + "pics/%s_ecg_sources.png" % subject)
    #     fig = ica.plot_components(ecg_inds, title=title % ('ecg', subject),
    #                               colorbar=True)
    #     fig.savefig(data_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(data_folder + "pics/%s_ecg_sources_ave.png" % subject)
    #
    # # plot ECG cleaning
    # ica.plot_overlay(ecg_evoked, exclude=ecg_inds)
    # fig.savefig(data_folder + "pics/%s_ecg_sources_clean_ave.png" % subject)

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

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

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

    del eog_epochs

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

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import apply [as 别名]
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
raw_ica.save(save_folder + "%s_%s_filtered_ica_mc_tsss-raw.fif" % (
    subject, condition),
             overwrite=True)

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

#.........这里部分代码省略.........
        fig9 = ica.plot_sources(eog_evoked, exclude=eog_inds, show=is_show)  # plot EOG sources + selection
        fig10 = ica.plot_overlay(eog_evoked, exclude=eog_inds, show=is_show)  # plot EOG cleaning

        fig = [fig9, fig10]
        report.add_figs_to_section(fig, captions=['Time-locked EOG sources',
                                                  'EOG overlay'], section = 'ICA - EOG')

    fig11 = ica.plot_overlay(raw, show=is_show)
    report.add_figs_to_section(fig11, captions=['Signal'], section = 'Signal quality') 
   
    ### plot all topographies and time seris of the ICA components
    n_ica_components = ica.mixing_matrix_.shape[1]
    
    n_topo = 10;
    n_fig  = n_ica_components/n_topo;
    n_plot = n_ica_components%n_topo;

    print '*************** n_fig = ' + str(n_fig) + ' n_plot = ' + str(n_plot) + '********************'
    fig = []
    t_start = 0
    t_stop = None # 60 if we want to take the fist 60s
    for n in range(0,n_fig):
        fig_tmp = ica.plot_components(range(n_topo*n,n_topo*(n+1)),title='ICA components', show=is_show)    
        fig.append(fig_tmp)
        fig_tmp = ica.plot_sources(raw, range(n_topo*n,n_topo*(n+1)), 
                                    start = t_start, stop  = t_stop, 
                                    title='ICA components')     
        fig.append(fig_tmp)
    
#    if n_plot > 0:
#        fig_tmp = ica.plot_components(range(n_fig*n_topo,n_ica_components), title='ICA components', show=is_show)    
#        fig.append(fig_tmp)
#        fig_tmp = ica.plot_sources(raw, range(n_fig*n_topo,n_ica_components), 
#                                        start = t_start, stop  = t_stop,
#                                        title='ICA components')     
#        fig.append(fig_tmp)   
#   
#    for n in range(0, len(fig)):
#        report.add_figs_to_section(fig[n], captions=['TOPO'], section = 'ICA Topo Maps')   
     
    if n_plot > 5:
        n_fig_l  = n_plot//5
                
        print '*************************** ' + str(n_fig_l) + ' *********************************'
        for n in range(0,n_fig_l):
            print range(n_fig*n_topo+5*n, n_fig*n_topo+5*(n+1))
            fig_tmp = ica.plot_components(range(n_fig*n_topo+5*n, n_fig*n_topo+5*(n+1)),title='ICA components')    
            fig.append(fig_tmp)
            fig_tmp = ica.plot_sources(raw, range(n_fig*n_topo+5*n, n_fig*n_topo+5*(n+1)), 
                                    start = t_start, stop  = t_stop, 
                                    title='ICA components')     
            fig.append(fig_tmp)
        
        print range(n_fig*n_topo+5*(n+1),n_ica_components)
        fig_tmp = ica.plot_components(range(n_fig*n_topo+5*(n+1),n_ica_components), title='ICA components')    
        fig.append(fig_tmp)
        fig_tmp = ica.plot_sources(raw, range(n_fig*n_topo+5*(n+1),n_ica_components), 
                                        start = t_start, stop  = t_stop, 
                                        title='ICA components')     
        fig.append(fig_tmp)   
        
    for n in range(0, len(fig)):
        report.add_figs_to_section(fig[n], captions=['TOPO'], section = 'ICA Topo Maps')       
    
    report_filename = os.path.join(subj_path,basename + "-report.html")
    print '******* ' + report_filename
    report.save(report_filename, open_browser=False, overwrite=True)
        
        
    # 3) apply ICA to raw data and save solution and report
    # check the amplitudes do not change
#    raw_ica_file = os.path.abspath(basename[:i_raw] + 'ica-raw.fif')
    raw_ica_file = os.path.join(subj_path, basename + '-preproc-raw.fif')
    raw_ica = ica.apply(raw)

    raw_ica.resample(sfreq=down_sfreq, npad=0)

    raw_ica.save(raw_ica_file, overwrite=True)

    # save ICA solution
    print ica_filename
    if has_ICA is False:
        ica.save(ica_filename)

    # 4) save data
    data_noIca, times = raw[select_sensors, :]
    data, times = raw_ica[select_sensors, :]

    print data.shape
    print raw.info['sfreq']

    ts_file = os.path.abspath(basename + "_ica.npy")
    np.save(ts_file, data)
    print '***** TS FILE ' + ts_file + '*****'

    if is_sensor_space:
        return ts_file, channel_coords_file, channel_names_file, raw.info['sfreq']
    else:
#        return raw_ica, channel_coords_file, channel_names_file, raw.info['sfreq']
        return raw_ica_file, channel_coords_file, channel_names_file, raw.info['sfreq']

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import apply [as 别名]
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).fit(raw, decim=6, 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
epochs_cln = ica.apply(epochs, copy=True)  # clean data, default in place

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

contrast = evoked[1] - evoked[0]

evoked.append(contrast)

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

plt.show()

# estimate noise covarariance
noise_cov = mne.compute_covariance(epochs_cln, tmax=0)