ENH: Add ssp_meg option (#595)

config.py

@@ -569,8 +569,8 @@ class ArbitraryContrast(TypedDict):
 
 mf_st_duration: Optional[float] = None
 """
-There are two kinds of maxfiltering: SSS (signal space separation) and tSSS
-(temporal signal space separation)
+There are two kinds of Maxwell filtering: SSS (signal space separation) and
+tSSS (temporal signal space separation)
 (see [Taulu et al., 2004](http://cds.cern.ch/record/709081/files/0401166.pdf)).
 
 If not None, apply spatiotemporal SSS (tSSS) with specified buffer
@@ -1063,6 +1063,17 @@ class ArbitraryContrast(TypedDict):
 on individual EOG epochs.
 """
 
+ssp_meg: Literal['separate', 'combined', 'auto'] = 'auto'
+"""
+Whether to compute SSP vectors for MEG channels separately (`'separate'`)
+or jointly (`'combined'`) for magnetometers and gradiomenters. When using
+Maxwell filtering, magnetometer and gradiometer signals are synthesized from
+multipole moments jointly and are no longer independent, so it can be useful to
+estimate projectors from all MEG sensors simultaneously. The default is
+`'auto'`, which will use `'combined'` when Maxwell filtering is used and
+`'separate'` otherwise.
+"""
+
 ssp_reject_ecg: Optional[
     Union[
         Dict[str, float],
@@ -2068,7 +2079,7 @@ def gen_log_kwargs(
 
 if (use_maxwell_filter and
         len(set(ch_types).intersection(('meg', 'grad', 'mag'))) == 0):
-    raise ValueError('Cannot use maxwell filter without MEG channels.')
+    raise ValueError('Cannot use Maxwell filter without MEG channels.')
 
 if (spatial_filter == 'ica' and
         ica_algorithm not in ('picard', 'fastica', 'extended_infomax')):

docs/source/changes.md

@@ -250,6 +250,10 @@ authors:
 - Patch information is now incorporated when computing surface source spaces,
   which should slightly improve the surface normals
   ({{ gh(588) }} by {{ authors.larsoner }})
+- Add [`ssp_meg`][config.ssp_meg] option for MEG SSP computation. This
+  defaults to `'auto'`, which will use `ssp_meg='combined'` for SSP computation
+  when Maxwell filtering is used.
+  ({{ gh(595) }} by {{ authors.larsoner }})
 
 ### Code health
 

docs/source/settings/preprocessing/ssp_ica.md

@@ -3,6 +3,7 @@
 ::: config.min_eog_epochs
 ::: config.n_proj_eog
 ::: config.n_proj_ecg
+::: config.ssp_meg
 ::: config.ecg_proj_from_average
 ::: config.eog_proj_from_average
 ::: config.ssp_reject_eog

scripts/preprocessing/_04b_run_ssp.py

@@ -60,6 +60,8 @@ def run_ssp(*, cfg, subject, session=None):
 
     ecg_projs = []
     ecg_epochs = create_ecg_epochs(raw)
+    if cfg.ssp_meg == 'auto':
+        cfg.ssp_meg = 'combined' if cfg.use_maxwell_filter else 'separate'
     if len(ecg_epochs) >= config.min_ecg_epochs:
         if cfg.ssp_reject_ecg == 'autoreject_global':
             reject_ecg_ = config.get_ssp_reject(
@@ -68,6 +70,7 @@ def run_ssp(*, cfg, subject, session=None):
             ecg_projs, _ = compute_proj_ecg(raw,
                                             average=cfg.ecg_proj_from_average,
                                             reject=reject_ecg_,
+                                            meg=cfg.ssp_meg,
                                             **cfg.n_proj_ecg)
         else:
             reject_ecg_ = config.get_ssp_reject(
@@ -76,6 +79,7 @@ def run_ssp(*, cfg, subject, session=None):
             ecg_projs, _ = compute_proj_ecg(raw,
                                             average=cfg.ecg_proj_from_average,
                                             reject=reject_ecg_,
+                                            meg=cfg.ssp_meg,
                                             **cfg.n_proj_ecg)
 
     if not ecg_projs:
@@ -141,6 +145,8 @@ def get_config(
         eog_proj_from_average=config.eog_proj_from_average,
         n_proj_eog=config.n_proj_eog,
         n_proj_ecg=config.n_proj_ecg,
+        ssp_meg=config.ssp_meg,
+        use_maxwell_filter=config.use_maxwell_filter,
     )
     return cfg
 

tests/configs/config_ds000248.py

@@ -37,6 +37,7 @@ def noise_cov(bp):
 spatial_filter = 'ssp'
 n_proj_eog = dict(n_mag=1, n_grad=1, n_eeg=1)
 n_proj_ecg = dict(n_mag=1, n_grad=1, n_eeg=0)
+ssp_meg = 'combined'
 ecg_proj_from_average = True
 eog_proj_from_average = False