Repository: mne-tools/mne-python

Issue:
Cortical signal suppression with magnetometers and EEG
### Describe the new feature or enhancement

The [cortical signal suppression](https://pubmed.ncbi.nlm.nih.gov/30604048/) technique takes advantage of the differential sensitivities to subcortical areas between channel types to filter out cortical signals. In the paper as well as in the already existing MNE-Python implementation, this is done with gradiometers and magnetometers. Going through the paper, however, I do not yet see any theoretical reason why it must necessarily be these two channel types. I would expect the same could be done with any two channel arrays with this differential sensitivity, and some amount of linear dependence, e.g. magnetometers and EEG. In the preliminary tests I have done, this does in fact seem to be the case; redoing the CSS example with mag and eeg channel types also produces a good result, though one must also reduce n_proj from 6 to 1 or 2.

Perhaps the functionality could be rewritten to offer the user more flexibility in applying this method? As it stands now, only systems with both magnetometers and gradiometers can take advantage of this very useful technique, which is a relatively limited amount of labs/datasets.

### Describe your proposed implementation

Rewrite the function such that the user can specify which two channel arrays to use, also maybe specifying some default behaviours, i.e. if left unspecified, checking which channel types are available and choosing appropriate ones, also changing the n_proj defaults accordingly (lower for mag+eeg, higher for grad+mag).

Some optional things:
1) Adding an option to return both the cortically filtered and subcortically filtered signal, in the event that someone actually wants to suppress the subcortical activity in the signal instead.
2) Expanding the example to show the added functionality.

### Describe possible alternatives


### Additional context

_No response_

Discussion:
larsoner: While in principle it's more general than `mag_picks` and `grad_picks`, I'm -0.5 because of the code churn we'll cause people and because the existing param names match the reference publication.

Patch:
diff --git a/doc/changes/devel/12725.newfeature.rst b/doc/changes/devel/12725.newfeature.rst
new file mode 100644
index 000000000..3775b851a
--- /dev/null
+++ b/doc/changes/devel/12725.newfeature.rst
@@ -0,0 +1,2 @@
+Add default ``spatial_colors="auto"`` to :func:`mne.viz.plot_evoked_white` and
+:meth:`mne.Evoked.plot_white` to enable spatial colors by default, by `Eric Larson`_.
diff --git a/mne/evoked.py b/mne/evoked.py
index 024a6d14f..37e4cedbf 100644
--- a/mne/evoked.py
+++ b/mne/evoked.py
@@ -745,6 +745,8 @@ class Evoked(
         time_unit="s",
         sphere=None,
         axes=None,
+        *,
+        spatial_colors="auto",
         verbose=None,
     ):
         return plot_evoked_white(
@@ -755,6 +757,7 @@ class Evoked(
             time_unit=time_unit,
             sphere=sphere,
             axes=axes,
+            spatial_colors=spatial_colors,
             verbose=verbose,
         )
 
diff --git a/mne/preprocessing/_css.py b/mne/preprocessing/_css.py
index f3f4ba09f..ebd9a77ee 100644
--- a/mne/preprocessing/_css.py
+++ b/mne/preprocessing/_css.py
@@ -43,12 +43,12 @@ def cortical_signal_suppression(
         gradiometer, and EEG channels.
     %(picks_good_data)s
     mag_picks : array-like of int
-        Array of the magnetometer channel indices that will be used to find
-        the reference data. If None (default), all magnetometers will
+        Array of the first set of channel indices that will be used to find
+        the common temporal subspace. If None (default), all magnetometers will
         be used.
     grad_picks : array-like of int
-        Array of the gradiometer channel indices that will be used to find
-        the reference data. If None (default), all gradiometers will
+        Array of the second set of channel indices that will be used to find
+        the common temporal subspace. If None (default), all gradiometers will
         be used.
     n_proj : int
         The number of projection vectors.
@@ -57,15 +57,19 @@ def cortical_signal_suppression(
     Returns
     -------
     evoked_subcortical : instance of Evoked
-        The evoked object with cortical contributions to the EEG data
-        suppressed.
+        The evoked object with contributions from the ``mag_picks`` and ``grad_picks``
+        channels removed from the ``picks`` channels.
 
     Notes
     -----
-    This method removes the common signal subspace between the magnetometer
-    data and the gradiometer data from the EEG data. This is done by a temporal
-    projection using ``n_proj`` number of projection vectors. For reference,
-    see :footcite:`Samuelsson2019`.
+    This method removes the common signal subspace between two sets of
+    channels (``mag_picks`` and ``grad_picks``) from a set of channels
+    (``picks``) via a temporal projection using ``n_proj`` number of
+    projection vectors. In the reference publication :footcite:`Samuelsson2019`,
+    the joint subspace between magnetometers and gradiometers is used to
+    suppress the cortical signal in the EEG data. In principle, other
+    combinations of sensor types (or channels) could be used to suppress
+    signals from other sources.
 
     References
     ----------
diff --git a/mne/utils/docs.py b/mne/utils/docs.py
index 3cb0042ed..0f1236617 100644
--- a/mne/utils/docs.py
+++ b/mne/utils/docs.py
@@ -4131,6 +4131,17 @@ smooth : float in [0, 1)
     The smoothing factor to be applied. Default 0 is no smoothing.
 """
 
+docdict["spatial_colors"] = """\
+spatial_colors : bool | 'auto'
+    If True, the lines are color coded by mapping physical sensor
+    coordinates into color values. Spatially similar channels will have
+    similar colors. Bad channels will be dotted. If False, the good
+    channels are plotted black and bad channels red. If ``'auto'``, uses
+    True if channel locations are present, and False if channel locations
+    are missing or if the data contains only a single channel. Defaults to
+    ``'auto'``.
+"""
+
 docdict["spatial_colors_psd"] = """\
 spatial_colors : bool
     Whether to color spectrum lines by channel location. Ignored if
diff --git a/mne/viz/evoked.py b/mne/viz/evoked.py
index 3707c0e42..8bad99ed2 100644
--- a/mne/viz/evoked.py
+++ b/mne/viz/evoked.py
@@ -1053,14 +1053,7 @@ def plot_evoked(
            Plot GFP for EEG instead of RMS. Label RMS traces correctly as such.
     window_title : str | None
         The title to put at the top of the figure.
-    spatial_colors : bool | 'auto'
-        If True, the lines are color coded by mapping physical sensor
-        coordinates into color values. Spatially similar channels will have
-        similar colors. Bad channels will be dotted. If False, the good
-        channels are plotted black and bad channels red. If ``'auto'``, uses
-        True if channel locations are present, and False if channel locations
-        are missing or if the data contains only a single channel. Defaults to
-        ``'auto'``.
+    %(spatial_colors)s
     zorder : str | callable
         Which channels to put in the front or back. Only matters if
         ``spatial_colors`` is used.
@@ -1494,6 +1487,8 @@ def plot_evoked_white(
     time_unit="s",
     sphere=None,
     axes=None,
+    *,
+    spatial_colors="auto",
     verbose=None,
 ):
     """Plot whitened evoked response.
@@ -1522,6 +1517,9 @@ def plot_evoked_white(
         List of axes to plot into.
 
         .. versionadded:: 0.21.0
+    %(spatial_colors)s
+
+        .. versionadded:: 1.8.0
     %(verbose)s
 
     Returns
@@ -1644,32 +1642,38 @@ def plot_evoked_white(
         raise RuntimeError("Wrong axes inputs")
 
     titles_ = _handle_default("titles")
-    if has_sss:
-        titles_["meg"] = "MEG (combined)"
-
     colors = [plt.cm.Set1(i) for i in np.linspace(0, 0.5, len(noise_cov))]
     ch_colors = _handle_default("color", None)
     iter_gfp = zip(evokeds_white, noise_cov, rank_list, colors)
 
-    # the first is by law the best noise cov, on the left we plot that one.
-    if not has_sss:
-        evokeds_white[0].plot(
-            unit=False,
-            axes=axes_evoked,
-            hline=[-1.96, 1.96],
-            show=False,
-            time_unit=time_unit,
-            spatial_colors=False,
-        )
-    else:
-        for (ch_type, picks), ax in zip(picks_list, axes_evoked):
-            ax.plot(times, evokeds_white[0].data[picks].T, color="k", lw=0.5)
-            for hline in [-1.96, 1.96]:
-                ax.axhline(hline, color="red", linestyle="--", lw=2)
-            ax.set(
-                title="%s (%d channel%s)"
-                % (titles_[ch_type], len(picks), _pl(len(picks)))
-            )
+    # The first is by law the best noise cov, on the left we plot that one.
+    # When we have data in SSS / MEG-combined mode, we have to do some info
+    # hacks to get it to plot all channels in the same axes, namely setting
+    # the channel unit (most important) and coil type (for consistency) of
+    # all MEG channels to be the same.
+    meg_idx = sss_title = None
+    if has_sss:
+        titles_["meg"] = "MEG (combined)"
+        meg_idx = [
+            pi for pi, (ch_type, _) in enumerate(picks_list) if ch_type == "meg"
+        ][0]
+        # Hack the MEG channels to all be the same type so they get plotted together
+        picks = picks_list[meg_idx][1]
+        for key in ("coil_type", "unit"):  # update both
+            use = evokeds_white[0].info["chs"][picks[0]][key]
+            for pick in picks:
+                evokeds_white[0].info["chs"][pick][key] = use
+        sss_title = f"{titles_['meg']} ({len(picks)} channel{_pl(picks)})"
+    evokeds_white[0].plot(
+        unit=False,
+        axes=axes_evoked,
+        hline=[-1.96, 1.96],
+        show=False,
+        time_unit=time_unit,
+        spatial_colors=spatial_colors,
+    )
+    if has_sss:
+        axes_evoked[meg_idx].set(title=sss_title)
 
     # Now plot the GFP for all covs if indicated.
     for evoked_white, noise_cov, rank_, color in iter_gfp:
diff --git a/tutorials/evoked/40_whitened.py b/tutorials/evoked/40_whitened.py
index 84c989541..a254111ce 100644
--- a/tutorials/evoked/40_whitened.py
+++ b/tutorials/evoked/40_whitened.py
@@ -5,8 +5,7 @@
 Plotting whitened data
 ======================
 
-This tutorial demonstrates how to plot :term:`whitened <whitening>`
-evoked data.
+This tutorial demonstrates how to plot :term:`whitened <whitening>` evoked data.
 
 Data are whitened for many processes, including dipole fitting, source
 localization and some decoding algorithms. Viewing whitened data thus gives
diff --git a/tutorials/forward/90_compute_covariance.py b/tutorials/forward/90_compute_covariance.py
index 637d02ddb..da0d40d01 100644
--- a/tutorials/forward/90_compute_covariance.py
+++ b/tutorials/forward/90_compute_covariance.py
@@ -21,7 +21,7 @@ import mne
 from mne.datasets import sample
 
 # %%
-# Source estimation method such as MNE require a noise estimations from the
+# Source estimation methods such as MNE require noise estimates from the
 # recordings. In this tutorial we cover the basics of noise covariance and
 # construct a noise covariance matrix that can be used when computing the
 # inverse solution. For more information, see :ref:`minimum_norm_estimates`.
diff --git a/tutorials/inverse/30_mne_dspm_loreta.py b/tutorials/inverse/30_mne_dspm_loreta.py
index 581029126..ce44292dc 100644
--- a/tutorials/inverse/30_mne_dspm_loreta.py
+++ b/tutorials/inverse/30_mne_dspm_loreta.py
@@ -32,7 +32,7 @@ events = mne.find_events(raw, stim_channel="STI 014")
 event_id = dict(aud_l=1)  # event trigger and conditions
 tmin = -0.2  # start of each epoch (200ms before the trigger)
 tmax = 0.5  # end of each epoch (500ms after the trigger)
-raw.info["bads"] = ["MEG 2443", "EEG 053"]
+raw.info["bads"] = ["MEG 2443", "EEG 053"]  # mark known bad channels
 baseline = (None, 0)  # means from the first instant to t = 0
 reject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)
 
diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py
index ed9a14fc5..c76a5c5f5 100644
--- a/tutorials/inverse/80_brainstorm_phantom_elekta.py
+++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py
@@ -46,7 +46,7 @@ raw = read_raw_fif(raw_fname)
 
 events = find_events(raw, "STI201")
 raw.plot(events=events)
-raw.info["bads"] = ["MEG1933", "MEG2421"]
+raw.info["bads"] = ["MEG1933", "MEG2421"]  # known bad channels
 
 # %%
 # The data has strong line frequency (60 Hz and harmonics) and cHPI coil


Test patch:
diff --git a/mne/preprocessing/tests/test_css.py b/mne/preprocessing/tests/test_css.py
index d68909083..f9293b249 100644
--- a/mne/preprocessing/tests/test_css.py
+++ b/mne/preprocessing/tests/test_css.py
@@ -34,11 +34,13 @@ def test_cortical_signal_suppression():
     ave.data[eeg_ind][1, :] = np.sin(2 * np.pi * 239 * ave.times) * np.mean(
         np.abs(ave.data[eeg_ind][1, :])
     )
-    ave_f = cortical_signal_suppression(ave)
-    cort_power = np.sum(np.abs(ave.data[eeg_ind][0, :]))
-    deep_power = np.sum(np.abs(ave.data[eeg_ind][1, :]))
-    cort_power_f = np.sum(np.abs(ave_f.data[eeg_ind][0, :]))
-    deep_power_f = np.sum(np.abs(ave_f.data[eeg_ind][1, :]))
-    rel_SNR_gain = (deep_power_f / deep_power) / (cort_power_f / cort_power)
-    assert rel_SNR_gain > 0
-    assert ave_f.data.shape == ave.data.shape
+    # include test for gh-12373, that you can use MAG+EEG if you want
+    for mag_picks, ind in ((None, eeg_ind), ("eeg", mag_ind)):
+        ave_f = cortical_signal_suppression(ave, mag_picks=mag_picks)
+        assert ave_f.data.shape == ave.data.shape
+        cort_power = np.linalg.norm(ave.data[ind][0, :])
+        deep_power = np.linalg.norm(ave.data[ind][1, :])
+        cort_power_f = np.linalg.norm(ave_f.data[ind][0, :])
+        deep_power_f = np.linalg.norm(ave_f.data[ind][1, :])
+        rel_SNR_gain = (deep_power_f / deep_power) / (cort_power_f / cort_power)
+        assert rel_SNR_gain > 3, mag_picks
diff --git a/mne/viz/tests/test_evoked.py b/mne/viz/tests/test_evoked.py
index d4053315f..0d062fb67 100644
--- a/mne/viz/tests/test_evoked.py
+++ b/mne/viz/tests/test_evoked.py
@@ -376,7 +376,7 @@ def test_plot_white():
         evoked.plot_white(cov, rank={"mag": 10})
     evoked.plot_white(cov, rank={"mag": 1, "grad": 8, "eeg": 2}, time_unit="s")
     fig = evoked.plot_white(cov, rank={"mag": 1}, time_unit="s")  # test rank
-    evoked.plot_white(cov, rank={"grad": 8}, time_unit="s", axes=fig.axes)
+    evoked.plot_white(cov, rank={"grad": 8}, time_unit="s", axes=fig.axes[:4])
     with pytest.raises(ValueError, match=r"must have shape \(4,\), got \(2,"):
         evoked.plot_white(cov, axes=fig.axes[:2])
     with pytest.raises(ValueError, match="When not using SSS"):
@@ -385,18 +385,19 @@ def test_plot_white():
     plt.close("all")
 
     fig = plot_evoked_white(evoked, [cov, cov])
-    assert len(fig.axes) == 3 * 2
-    axes = np.array(fig.axes).reshape(3, 2)
+    assert len(fig.axes) == 3 * 3
+    axes = np.array(fig.axes[:6]).reshape(3, 2)
     plot_evoked_white(evoked, [cov, cov], axes=axes)
     with pytest.raises(ValueError, match=r"have shape \(3, 2\), got"):
         plot_evoked_white(evoked, [cov, cov], axes=axes[:, :1])
 
     # Hack to test plotting of maxfiltered data
-    evoked_sss = _get_epochs(picks="meg").average()
+    evoked_sss = _get_epochs(picks=("meg", "eeg")).average()
+    evoked_sss.set_eeg_reference(projection=True).apply_proj()
     sss = dict(sss_info=dict(in_order=80, components=np.arange(80)))
     with evoked_sss.info._unlock():
         evoked_sss.info["proc_history"] = [dict(max_info=sss)]
-    evoked_sss.plot_white(cov, rank={"meg": 64})
+    evoked_sss.plot_white([cov, cov], rank={"meg": 64})
     with pytest.raises(ValueError, match="When using SSS"):
         evoked_sss.plot_white(cov, rank={"grad": 201})
     evoked_sss.plot_white(cov, time_unit="s")