iBOTS Learn: Spatial Filters: PCA and ICA

Spatial Filters: PCA and ICA

Author

Dr. Ole Bialas

Download Materials

import numpy as np
from matplotlib import pyplot as plt
import mne
import mne_icalabel
from mne.preprocessing import ICA

Due to volume conduction the electrical field generated by neural activity spreads across the scalp and is picked up by all EEG sensors simultaneously. Because the intensity of the signal depends on the proximity to the source, sensors close to the neural generator will pick up more of its signal than sensors further away. Spatial filters exploit this phenomenon by separating the signal into distinct sources, each of which is a weighted combination of all sensors:

Y = W \cdot X

Where $X$ is the channels-by-time matrix of EEG recordings, $Y$ is the components-by-time matrix of extracted sources and $W$ is the components-by-channels weight matrix that indicates how strongly each channel contributes to each component. The maximum number of components is equal to the number of channels in $X$ . There are many different methods for estimating the weight matrix $W$ and in this notebook you are going to learn about two of them: principal component analysis (PCA) and independent component analysis (ICA).

PCA finds the components that account for the largest amount of variance in the data. This means that usually, with only a few components, we can explain most of our data and additional components capture only noise. Thus, PCA can be useful for noise reduction and it can also make data analysis more efficient because we can focus on a few principal components instead of all recorded channels. However, because PCA requires the weight vectors to be orthogonal, it can only find sources whose scalp topographies are orthogonal to each other. If two neural sources have spatially overlapping topographies, PCA is forced to mix them together and cannot separate them.

ICA on the other hand tries to find components that maximize the statistical independence of the reconstructed sources. Two signals are said to be statistically independent if knowing one signal does not tell you anything about the other. Unlike PCA, ICA places no orthogonality constraint on the weight vectors, so it can find sources with overlapping scalp topographies as long as their time series are statistically independent. As a tradeoff, ICA is more computationally expensive and can lead to suboptimal results if there is a mismatch between the number of extracted components and the number of independent sources in the data.

In this notebook we are going to test PCA and ICA on simulated EEG data. This provides a scenario with a known ground truth: we know exactly which brain sources contribute to the EEG signal and can compare the results of PCA and ICA against this knowledge.

Setup

Utility Functions

The PCA class wraps the implementation from sklearn to provide the same user interface as MNE’s ICA class.

from sklearn.decomposition import PCA as _PCA

class PCA:
    def __init__(self, n_components):
        self.n_components = n_components
        self._pca = _PCA(n_components=n_components)

    def fit(self, inst):
        picks = mne.pick_types(inst.info, eeg=True)
        self.info = mne.pick_info(inst.info, picks)
        if isinstance(inst, mne.Evoked):
            X = inst.get_data(picks=picks).T           # (n_times, n_channels)
        else:
            X = inst.get_data(picks=picks)             # (n_epochs, n_channels, n_times)
            X = X.transpose(0, 2, 1).reshape(-1, len(picks))
        self._pca.fit(X)
        self.components_ = self._pca.components_       # (n_components, n_channels)
        self.explained_variance_ratio_ = self._pca.explained_variance_ratio_
        return self

    def get_sources(self, epochs):
        source_info = mne.create_info(
            ch_names=[f"PCA{i:03d}" for i in range(self.n_components)],
            sfreq=epochs.info["sfreq"],
            ch_types="misc",
        )
        picks = mne.pick_types(epochs.info, eeg=True)
        X = self.components_ @ epochs.get_data(picks=picks)
        return mne.EpochsArray(X, source_info, tmin=epochs.tmin, verbose=False)

    def apply(self, epochs, include=None):
        W = self.components_ if include is None else self.components_[include]
        proj = W.T @ W
        picks = mne.pick_types(epochs.info, eeg=True)
        data = epochs.get_data().copy()
        data[:, picks, :] = proj @ data[:, picks, :]
        return mne.EpochsArray(data, epochs.info, tmin=epochs.tmin, verbose=False)

    def get_components(self):
        return self.components_.T  # mixing matrix (n_channels, n_components)

    def plot_components(self, picks=None):
        indices = list(picks if picks is not None else range(self.n_components))
        n_cols = min(5, len(indices))
        n_rows = (len(indices) + n_cols - 1) // n_cols
        fig, axes = plt.subplots(n_rows, n_cols, figsize=(3 * n_cols, 3 * n_rows), squeeze=False)
        axes = axes.flatten()
        for ax, i in zip(axes, indices):
            mne.viz.plot_topomap(self.components_[i], self.info, axes=ax, show=False)
            ax.set_title(f"PCA{i:03d}")
        for ax in axes[len(indices):]:
            ax.set_visible(False)
        fig.suptitle("PCA components")
        return fig

Simulate Data

Run the cell below to simulate EEG data four with dipole sources using the approach described in the notebook in forward modeling.

%run simulation.py

    Reading a source space...
    [done]
    Reading a source space...
    [done]
    2 source spaces read
Adding noise to 64/64 channels (64 channels in cov)

Sphere                : origin at (0.0 0.0 0.0) mm
              radius  : 0.1 mm
Source location file  : dict()
Assuming input in millimeters
Assuming input in MRI coordinates

Positions (in meters) and orientations
2 sources
blink simulated and trace not stored
Setting up forward solutions
EEG channel type selected for re-referencing
Adding average EEG reference projection.
1 projection items deactivated
Average reference projection was added, but has not been applied yet. Use the apply_proj method to apply it.

Run the cell below to load the simulated epochs and the source activity stc.

epochs = mne.read_epochs("sim-epo.fif")
stc = mne.read_source_estimate("sim")

Reading /home/atleeri/repositories/new-learning-platform/notebooks/mne/02_temporal_and_spatial_filters/01_spatial_filters/sim-epo.fif ...
    Read a total of 1 projection items:
        Average EEG reference (1 x 64) active
    Found the data of interest:
        t =    -200.00 ...     800.00 ms
        0 CTF compensation matrices available
Not setting metadata
299 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 1)
1 projection items activated

The simulated data contains 4 sources:

An evoked response in the left auditory cortex
An evoked response in the right auditory cortex (with slightly different latency)
A steady 10 Hz (alpha) oscillation in the left visual cortex
A modulated 6 Hz (theta) oscillation in the left prefrontal cortex

In addition to these sources, the simulated data contains sensor noise and EOG (eyeblink) artifacts.

Run the cells below to plot the location and time course of the simulated sources.

fs_dir = mne.datasets.fetch_fsaverage(verbose=False)
brain = mne.viz.Brain("fsaverage", hemi="both", subjects_dir=fs_dir.parent, surf="inflated")
colors = ["blue", "red", "green", "orange"]
vertices = [(v, "lh") for v in stc.vertices[0]] + [(v, "rh") for v in stc.vertices[1]]
for (vertex, hemi), color in zip(vertices, colors):
    brain.add_foci(vertex, coords_as_verts=True, hemi=hemi, color=color)

Using pyvistaqt 3d backend.

2026-04-24 14:24:06.068 (   2.741s) [    7A2BE29A3440]vtkXOpenGLRenderWindow.:1539   ERR| The result is out of range, failed to get the converted tmp.

plt.plot(stc.times[:250], stc.data[:, :250].T)
plt.xlabel("Time [s]")
plt.ylabel("Current Dipole Moment [$A \\times m$]")
plt.legend(["Occipital alpha", "Left auditory", "Frontal theta", "Right auditory" ]);

Section 1: Selecting and Visualizing Principal Components

Background

Principal Component Analysis (PCA) decomposes the multichannel EEG signal into a set of principal components: orthogonal spatial filters ordered by the variance they explain. Each component is a weighted sum of electrode signals, and the corresponding component loading (time course) shows how strongly that spatial pattern is expressed at each moment.

Because components are constrained to be orthogonal, PCA can only separate sources that are both spatially and temporally uncorrelated. Sources with overlapping scalp distributions or correlated time courses — such as two auditory responses evoked by the same stimulus — will be merged into a single component.

Exercises

In this section, you are going to fit PCA to the simulated epoched data and select and visualize individual components using their time course and topographical distribution. This allows us to link the components to the underlying sources they capture. Here are the code snippets you’ll need:

Code	Description
`pca = PCA(n_components=k)` `pca.fit(epochs)`	Fit PCA to epoched EEG data
`pca.plot_components()`	Plot component scalp topographies
`pca.get_sources(epochs)`	Get component time courses as an Epochs object
`epochs_pca.average(picks=[i]).plot()`	Plot the average time course of component `i`

Example: Fit a PCA with 2 components to epochs and plot the components. These components capture the auditory and prefrontal sources.

Notice how the left and right auditory sources are captured in a single component because their time course is highly correlated.

pca = PCA(n_components=2)
pca.fit(epochs)
pca.plot_components();

Exercise: Fit a PCA with 4 components to epochs and plot the components. Which sources are the additional components capturing?

Solution

pca = PCA(n_components=4)
pca.fit(epochs)
pca.plot_components();

Exercise: Fit a PCA with 10 components to epochs and plot the components. What activity is captured by the additional components?

Solution

pca = PCA(n_components=10)
pca.fit(epochs)
pca.plot_components();

Example: Use pca.get_sources() to get the time course of the components (from the last PCA with 10 components) and plot the average time course of the first component.

epochs_pca = pca.get_sources(epochs)
epochs_pca.average(picks=[0]).plot();

Exercise: Use pca.get_sources() to get the time course of the components and plot the average time course of the second component.

Solution

epochs_pca = pca.get_sources(epochs)
epochs_pca.average(picks=[1]).plot();

Exercise: Use pca.get_sources() to get the time course of the components and plot the average time course of the last component.

Solution

epochs_pca = pca.get_sources(epochs)
epochs_pca.average(picks=[-1]).plot();

Section 2: Applying PCA Components as Spatial Filters

Background

PCA has two practical applications in EEG analysis. First, the explained variance ratio quantifies how much of the total signal variance each component accounts for. Plotting these values helps determine how many components are needed to capture the dominant signal, with a steep initial drop-off indicating a clear boundary between signal and noise.

Second, PCA can be used as a spatial filter. By projecting the data into component space and then projecting back into sensor space, we obtain a version of the data that retains only the selected components — effectively denoising the data in a data-driven way.

Exercises

In this section you are going to fit PCA to the epoched data and visualize the amount of variance in the data captured by the components. You are also going to use the PCA components as spatial filters by projecting them into sensor-space. The resulting data has the same number of channels as the original but only contains information from the selected components.

Code	Description
`pca.explained_variance_ratio_`	Fraction of total variance per component
`np.cumsum(pca.explained_variance_ratio_)`	Cumulative variance explained
`pca.apply(epochs, include=[i, j])`	Reconstruct epochs retaining only the selected components
`epochs.average().plot_joint()`	Plot the average ERP with scalp topomaps

Example: Compute a PCA with n_components=5 and plot the .explained_variance_ratio_ for every component.

pca = PCA(n_components=5)
pca.fit(epochs)

plt.plot(pca.explained_variance_ratio_)
plt.xlabel("Component Number")
plt.ylabel("Ratio of Variance Explained");

Exercise: Compute a PCA with n_components=15 and plot the .explained_variance_ratio_ for every component.

Solution

pca = PCA(n_components=15)
pca.fit(epochs)

plt.plot(pca.explained_variance_ratio_)
plt.xlabel("Component Number")
plt.ylabel("Ratio of Variance Explained")

Exercise: Plot the cumulative sum (np.cumsum) of the .explained_variance_ratio_. How much of the variance in the data is accounted for by the PCA components?

Solution

plt.plot(np.cumsum(pca.explained_variance_ratio_))

Exercise: Compute a PCA with n_components=15 on the average ERP (epochs.average()) and plot the .explained_variance_ratio_ as well as its cumulative sum.

Why does the PCA require fewer components to account for the ERP compared to the unaveraged epochs?

Solution

pca = PCA(n_components=15)
pca.fit(epochs.average())

plt.plot(np.cumsum(pca.explained_variance_ratio_))
plt.xlabel("Component Number")
plt.ylabel("Ratio of Variance Explained")

PCA can be used as a spatial filter by transforming the data using the PCA components and then transforming it back into sensor space. The resulting data has the same number of channels as the original but only contains information from the selected component(s).

Example: Use pca.apply() to reconstruct the epochs keeping only the first component. Then, plot the average ERP from the transformed epochs.

epochs_transformed = pca.apply(epochs, include=[0])
epochs_transformed.average().plot_joint();

Projections have already been applied. Setting proj attribute to True.

Exercise: Use pca.apply() to reconstruct the epochs keeping only the second component. Then, plot the average ERP from the transformed epochs.

Solution

epochs_transformed = pca.apply(epochs, include=[1])
epochs_transformed.average().plot_joint();

Projections have already been applied. Setting proj attribute to True.

Exercise: Use pca.apply() to reconstruct the epochs keeping only the first two components. Then, plot the average ERP from the transformed epochs.

Solution

epochs_transformed = pca.apply(epochs, include=[0, 1])
epochs_transformed.average().plot_joint();

Projections have already been applied. Setting proj attribute to True.

Exercise: How many components do you have to keep to make the transformed ERP look practically indistinguishable from the original plotted below?

epochs.average().plot_joint();

Projections have already been applied. Setting proj attribute to True.

Solution

epochs_transformed = pca.apply(epochs, include=[0, 1, 2, 3])
epochs_transformed.average().plot_joint();

Projections have already been applied. Setting proj attribute to True.

Section 3: Independent Component Analysis (ICA)

Background

Independent Component Analysis (ICA) is a blind source separation technique that decomposes the multichannel EEG signal into independent components — spatial filters optimized for statistical independence rather than orthogonality. Unlike PCA, ICA can separate sources with overlapping scalp distributions as long as their time courses are statistically independent.

ICA is particularly useful in EEG for two purposes: identifying and characterizing individual neural sources, and removing artifacts such as eye blinks and muscle activity that are captured in dedicated components.

The number of components should be chosen based on the expected number of sources in the data — using too few prevents clean source separation, while too many increases the risk of converging to a suboptimal solution.

Exercises

In this section you are going to fit ICA to the simulated epochs to observe how it performs compared to PCA. You are also going to apply ICA as a spatial filter to remove artifacts like eyeblinks and see how the ICLabel algorithm can be applied to automatically label components.

Code	Description
`ica = ICA(n_components=k)` `ica.fit(epochs)`	Fit ICA to epoched EEG data
`ica.plot_components()`	Plot component scalp topographies
`ica.get_sources(epochs)`	Get component time courses as an Epochs object
`ica.apply(epochs.copy(), exclude=[i])`	Remove component `i` and return cleaned epochs
`mne.combine_evoked([erp1, erp2], weights=[-1, 1])`	Compute the difference between two ERPs

Example: Fit an ICA on epochs with 3 components and plot the components. Note that the ICA isn’t separating the sources perfectly — the last component captures the prefrontal source but also includes activity from the occipital source.

ica = ICA(n_components=3, random_state=97)
ica.fit(epochs)
ica.plot_components();

Fitting ICA to data using 64 channels (please be patient, this may take a while)
    Applying projection operator with 1 vector (pre-whitener computation)
    Applying projection operator with 1 vector (pre-whitener application)
Selecting by number: 3 components
    Applying projection operator with 1 vector (pre-whitener application)
Fitting ICA took 0.6s.

Exercise: Fit an ICA to epochs with 7 components and plot the components. Which additional sources are you capturing?

Solution

ica = ICA(n_components=7, random_state=97)
ica.fit(epochs)
ica.plot_components();

Fitting ICA to data using 64 channels (please be patient, this may take a while)
    Applying projection operator with 1 vector (pre-whitener computation)
    Applying projection operator with 1 vector (pre-whitener application)
Selecting by number: 7 components
    Applying projection operator with 1 vector (pre-whitener application)
Fitting ICA took 0.5s.

Exercise: Fit an ICA on epochs with 6 components and plot the components. Can you see a separation of the left and right auditory sources?

Solution

ica = ICA(n_components=6, random_state=97)
ica.fit(epochs)
ica.plot_components();

Fitting ICA to data using 64 channels (please be patient, this may take a while)
    Applying projection operator with 1 vector (pre-whitener computation)
    Applying projection operator with 1 vector (pre-whitener application)
Selecting by number: 6 components
    Applying projection operator with 1 vector (pre-whitener application)
Fitting ICA took 0.6s.

Exercise: Use ica.get_sources() to get the time course of the components and plot the average time course of the first two components (works exactly the same as with pca.get_sources()).

Solution

epochs_ica = ica.get_sources(epochs)
epochs_ica.average(picks=[0, 1]).plot();

    Applying projection operator with 1 vector (pre-whitener application)

Component 4 represents EOG (eyeblink) artifacts. These are ubiquitous in EEG recordings and can be recognized by the focal frontal topography. We can also plot the time course of the component. If it represents EOG artifacts, it should be flat for most of the signal with large spikes in epochs where the (simulated) participant blinked.

epochs_ica = ica.get_sources(epochs)
epochs_ica[152:160].plot(picks=[4], scalings=8);

    Applying projection operator with 1 vector (pre-whitener application)
Using qt as 2D backend.

Example: Apply an ICA filter to a copy of epochs and exclude the EOG component to obtain epochs_clean. Then, compute the difference ERP between the original and cleaned data and plot it to see what was removed.

epochs_clean = ica.apply(epochs.copy(), exclude=[4])
diff = mne.combine_evoked([epochs.average(), epochs_clean.average()], weights=(-1, 1))
diff.plot();

Applying ICA to Epochs instance
    Applying projection operator with 1 vector (pre-whitener application)
    Transforming to ICA space (6 components)
    Zeroing out 1 ICA component
    Projecting back using 64 PCA components

Exercise: Apply an ICA filter to a copy of epochs and exclude component 5 to obtain epochs_clean. Then, compute the difference ERP between the original and cleaned data and plot it to see what was removed.

Solution

epochs_clean = ica.apply(epochs.copy(), exclude=[5])
diff = mne.combine_evoked([epochs.average(), epochs_clean.average()], weights=(-1, 1))
diff.plot();

Applying ICA to Epochs instance
    Applying projection operator with 1 vector (pre-whitener application)
    Transforming to ICA space (6 components)
    Zeroing out 1 ICA component
    Projecting back using 64 PCA components

Demo: The ICLabel algorithm (Pion-Tonachini et al., 2019) is an automated procedure to label ICA components as brain signals or different kinds of artifacts. It uses a classifier trained on a crowdsourced database of manually labeled components. The project also provides a tutorial website that allows you to practice labeling components.

The code below applies ICLabel to the ica fitted on epochs and prints the labels stored in ica.labels_.

mne_icalabel.label_components(epochs, ica, method="iclabel")
ica.labels_

{'brain': [np.int64(0), np.int64(1), np.int64(2), np.int64(3)],
 'muscle': [],
 'eog': [np.int64(4)],
 'ecg': [],
 'line_noise': [],
 'ch_noise': [],
 'other': [np.int64(5)]}

Section 4: Combining Components Across Recordings

Background

A practical challenge in ICA-based analysis is that components are not guaranteed to be consistent across recordings: the order, polarity, and even the number of components may differ between participants or sessions. This makes it difficult to directly compare or average ICA results across recordings.

The corrmap() function addresses this by using a template component from one recording and searching for the most spatially correlated component in each of the other recordings. It handles sign flips automatically and stores the matched component indices in the .labels_ attribute of each ICA object, making it straightforward to apply consistent artifact removal or source selection across datasets.

Exercises

In this section you are going to fit ICA to two separate halves of the simulated data and use corrmap() to identify the same neural sources across both recordings.

Code	Description
`ica1 = ICA(n_components=k)` `ica1.fit(epochs1)`	Fit ICA to the first recording
`mne.preprocessing.corrmap([ica1, ica2], template=(0, i), label="name")`	Find components matching template component `i` of `ica1` across all ICA objects
`ica1.labels_`	Dictionary mapping label names to matched component indices

Run the cell below to split epochs into two equally long segments epochs1 and epochs2. Assume these are separate recordings we want to analyze with ICA.

n = len(epochs) // 2
epochs1 = epochs[:n]
epochs2 = epochs[n:]

Example: Fit ica1 with 5 components to epochs1 and plot the components.

ica1 = ICA(n_components=5, random_state=101)
ica1.fit(epochs1)
ica1.plot_components();

Fitting ICA to data using 64 channels (please be patient, this may take a while)
    Applying projection operator with 1 vector (pre-whitener computation)
    Applying projection operator with 1 vector (pre-whitener application)
Selecting by number: 5 components
    Applying projection operator with 1 vector (pre-whitener application)
Fitting ICA took 0.5s.

Exercise: Fit ica2 with 5 components to epochs2 and plot the components. Notice how both the order and polarity of the components may differ.

Solution

ica2 = ICA(n_components=5, random_state=101)
ica2.fit(epochs2)
ica2.plot_components();

Fitting ICA to data using 64 channels (please be patient, this may take a while)
    Applying projection operator with 1 vector (pre-whitener computation)
    Applying projection operator with 1 vector (pre-whitener application)
Selecting by number: 5 components
    Applying projection operator with 1 vector (pre-whitener application)
Fitting ICA took 0.2s.

Example: Use the corrmap() function to find matching components across ica1 and ica2. The template argument determines the template used for searching. The first number refers to the ICA and the second number to the component - (0, 0) refers to the first component in ica1. The indices of the matching components are stored in the .labels_ of the ICA objects.

mne.preprocessing.corrmap([ica1, ica2], template=(0, 0), label="auditory/left")
ica1.labels_, ica2.labels_

Median correlation with constructed map: 1.000
Displaying selected ICs per subject.
At least 1 IC detected for each subject.

({'auditory/left': [np.int64(0)]}, {'auditory/left': [np.int64(0)]})

Exercise: Use the corrmap() function to find matching components from the right auditory components (component 2 in ica1) and label them as "auditory/right".

Solution

mne.preprocessing.corrmap([ica1, ica2], template=(0,1), label="auditory/right")
ica1.labels_, ica2.labels_

Median correlation with constructed map: 0.998
Displaying selected ICs per subject.
At least 1 IC detected for each subject.

({'auditory/left': [np.int64(0)], 'auditory/right': [np.int64(1)]},
 {'auditory/left': [np.int64(0)], 'auditory/right': [np.int64(2)]})

Section 5: Bonus: Applying ICA to Real EEG Data

Now that you have built some intuition for how PCA and ICA work it’s time to apply these methods to real data. In this section you are going to combine what you have learned so far by applying PCA and ICA to remove eyeblink artifacts from real EEG data.

Here’s a reminder of the essential code snippets:

Code	Description
`pca = PCA(n_components=k)` `pca.fit(epochs)`	Fit PCA to epoched EEG data
`np.cumsum(pca.explained_variance_ratio_)`	Cumulative variance explained
`ica = ICA(n_components=k)` `ica.fit(epochs)`	Fit ICA to epoched EEG data
`ica.plot_components()`	Plot component scalp topographies
`mne_icalabel.label_components(epochs, ica, method="iclabel")`	Classify components using ICLabel and store results in `ica.labels_`
`ica.apply(epochs.copy(), include=ica.labels_["brain"])`	Retain only brain components and return cleaned epochs
`mne.combine_evoked([erp1, erp2], weights=[-1, 1])`	Compute the difference between two ERPs

Run the cell below to load and epoch the MNE sample data and plot the average auditory ERP.

data_path = mne.datasets.sample.data_path()
meg_path = data_path / "MEG" / "sample"
raw_fname = meg_path / "sample_audvis_raw.fif"
event_fname = meg_path / "sample_audvis_raw-eve.fif"
raw = mne.io.read_raw_fif(raw_fname, preload=True, verbose=False)
raw.filter(1, 40, verbose=False)
event_id = {
    "auditory/left": 1,
    "auditory/right": 2,
}
epochs = mne.Epochs(
    raw,
    events= mne.read_events(event_fname),
    event_id = event_id,
    tmin=-0.2,
    tmax=0.5,
    picks = mne.pick_types(raw.info, eeg=True),
    decim=4,
    baseline=None,
    preload=True,
    verbose=False
)
epochs["auditory/right"].average().plot();

Exercise: Fit PCA to epochs with the maximum number of components (59) and plot the cumulative sum of the explained variance ratio. How does this compare to the simulated data?

Solution

pca = PCA(n_components=59)
pca.fit(epochs)
plt.plot(np.cumsum(pca.explained_variance_ratio_));

Exercise: Fit an ICA to epochs and choose the number of components so that you capture most of the variance in the data (use the plot from the previous exercise to inform your decision).

Solution

ica = ICA(n_components=30, random_state=101)
ica.fit(epochs)
ica.plot_components();

Fitting ICA to data using 59 channels (please be patient, this may take a while)
Selecting by number: 30 components
Fitting ICA took 1.8s.

Exercise: Identify components that correspond to eyeblinks and use ica.apply() to apply a spatial filter that removes these components.

Bonus: Instead of manually selecting components, use mne_icalabel.label_components and only include components labeled as "brain".

Solution

mne_icalabel.label_components(epochs, ica, method="iclabel")
epochs_clean = ica.apply(epochs.copy(), include=ica.labels_["brain"]);

Applying ICA to Epochs instance
    Transforming to ICA space (30 components)
    Zeroing out 12 ICA components
    Projecting back using 59 PCA components

Exercise: Compute the difference ERP between the original and cleaned epochs and plot it to visualize the activity that was removed.

Solution

diff = mne.combine_evoked([epochs.average(), epochs_clean.average()], weights=(-1, 1))
diff.plot();