Aluode/PerceptionLabPortable
0
1# Authors: The MNE-Python contributors.2# License: BSD-3-Clause3# Copyright the MNE-Python contributors.4 5import collections.abc as abc6from functools import partial7 8import numpy as np9 10from .._fiff.meas_info import Info, create_info11from .._fiff.pick import _picks_to_idx12from ..filter import filter_data13from ..utils import (14 _validate_type,15 fill_doc,16 logger,17)18from ._covs_ged import _ssd_estimate19from ._mod_ged import _get_spectral_ratio, _ssd_mod20from .base import _GEDTransformer21 22 23@fill_doc24class SSD(_GEDTransformer):25 """26 Signal decomposition using the Spatio-Spectral Decomposition (SSD).27 28 SSD seeks to maximize the power at a frequency band of interest while29 simultaneously minimizing it at the flanking (surrounding) frequency bins30 (considered noise). It extremizes the covariance matrices associated with31 signal and noise :footcite:`NikulinEtAl2011`.32 33 SSD can either be used as a dimensionality reduction method or a34 ‘denoised’ low rank factorization method :footcite:`HaufeEtAl2014b`.35 36 Parameters37 ----------38 %(info_not_none)s Must match the input data.39 filt_params_signal : dict40 Filtering for the frequencies of interest.41 filt_params_noise : dict42 Filtering for the frequencies of non-interest.43 reg : float | str | None (default)44 Which covariance estimator to use.45 If not None (same as 'empirical'), allow regularization for covariance46 estimation. If float, shrinkage is used (0 <= shrinkage <= 1). For str47 options, reg will be passed to method :func:`mne.compute_covariance`.48 n_components : int | None (default None)49 The number of components to extract from the signal.50 If None, the number of components equal to the rank of the data are51 returned (see ``rank``).52 picks : array of int | None (default None)53 The indices of good channels.54 sort_by_spectral_ratio : bool (default True)55 If set to True, the components are sorted according to the spectral56 ratio.57 See Eq. (24) in :footcite:`NikulinEtAl2011`.58 return_filtered : bool (default False)59 If return_filtered is True, data is bandpassed and projected onto the60 SSD components.61 n_fft : int (default None)62 If sort_by_spectral_ratio is set to True, then the SSD sources will be63 sorted according to their spectral ratio which is calculated based on64 :func:`mne.time_frequency.psd_array_welch`. The n_fft parameter sets the65 length of FFT used. The default (None) will use 1 second of data.66 See :func:`mne.time_frequency.psd_array_welch` for more information.67 cov_method_params : dict | None (default None)68 As in :class:`mne.decoding.SPoC`69 The default is None.70 restr_type : "restricting" | "whitening" | "ssd" | None71 Restricting transformation for covariance matrices before performing72 generalized eigendecomposition.73 If "restricting" only restriction to the principal subspace of signal_cov74 will be performed.75 If "whitening", covariance matrices will be additionally rescaled according76 to the whitening for the signal_cov.77 If "ssd", simplified version of "whitening" is performed.78 If None, no restriction will be applied. Defaults to "ssd".79 80 .. versionadded:: 1.1181 rank : None | dict | ‘info’ | ‘full’82 As in :class:`mne.decoding.SPoC`83 This controls the rank computation that can be read from the84 measurement info or estimated from the data, which determines the85 maximum possible number of components.86 See Notes of :func:`mne.compute_rank` for details.87 We recommend to use 'full' when working with epoched data.88 89 Attributes90 ----------91 filters_ : array, shape (``n_channels or less``, n_channels)92 The spatial filters to be multiplied with the signal.93 patterns_ : array, shape (``n_channels or less``, n_channels)94 The patterns for reconstructing the signal from the filtered data.95 96 References97 ----------98 .. footbibliography::99 """100 101 def __init__(102 self,103 info,104 filt_params_signal,105 filt_params_noise,106 reg=None,107 n_components=None,108 picks=None,109 sort_by_spectral_ratio=True,110 return_filtered=False,111 n_fft=None,112 cov_method_params=None,113 *,114 restr_type="whitening",115 rank=None,116 ):117 """Initialize instance."""118 self.info = info119 self.filt_params_signal = filt_params_signal120 self.filt_params_noise = filt_params_noise121 self.reg = reg122 self.n_components = n_components123 self.picks = picks124 self.sort_by_spectral_ratio = sort_by_spectral_ratio125 self.return_filtered = return_filtered126 self.n_fft = n_fft127 self.cov_method_params = cov_method_params128 self.restr_type = restr_type129 self.rank = rank130 131 cov_callable = partial(132 _ssd_estimate,133 reg=reg,134 cov_method_params=cov_method_params,135 info=info,136 picks=picks,137 n_fft=n_fft,138 filt_params_signal=filt_params_signal,139 filt_params_noise=filt_params_noise,140 rank=rank,141 sort_by_spectral_ratio=sort_by_spectral_ratio,142 )143 super().__init__(144 n_components=n_components,145 cov_callable=cov_callable,146 mod_ged_callable=_ssd_mod,147 restr_type=restr_type,148 )149 150 def _validate_params(self, X):151 if isinstance(self.info, float): # special case, mostly for testing152 self.sfreq_ = self.info153 else:154 _validate_type(self.info, Info, "info")155 self.sfreq_ = self.info["sfreq"]156 dicts = {"signal": self.filt_params_signal, "noise": self.filt_params_noise}157 for param, dd in [("l", 0), ("h", 0), ("l", 1), ("h", 1)]:158 key = ("signal", "noise")[dd]159 if param + "_freq" not in dicts[key]:160 raise ValueError(161 f"{param + '_freq'} must be defined in filter parameters for {key}"162 )163 val = dicts[key][param + "_freq"]164 if not isinstance(val, int | float):165 _validate_type(val, ("numeric",), f"{key} {param}_freq")166 # check freq bands167 if (168 self.filt_params_noise["l_freq"] > self.filt_params_signal["l_freq"]169 or self.filt_params_signal["h_freq"] > self.filt_params_noise["h_freq"]170 ):171 raise ValueError(172 "Wrongly specified frequency bands!\n"173 "The signal band-pass must be within the noise "174 "band-pass!"175 )176 self.freqs_signal_ = (177 self.filt_params_signal["l_freq"],178 self.filt_params_signal["h_freq"],179 )180 self.freqs_noise_ = (181 self.filt_params_noise["l_freq"],182 self.filt_params_noise["h_freq"],183 )184 _validate_type(self.sort_by_spectral_ratio, (bool,), "sort_by_spectral_ratio")185 _validate_type(self.n_fft, ("numeric", None), "n_fft")186 self.n_fft_ = min(187 int(self.n_fft if self.n_fft is not None else self.sfreq_),188 X.shape[-1],189 )190 _validate_type(self.return_filtered, (bool,), "return_filtered")191 if isinstance(self.info, Info):192 ch_types = self.info.get_channel_types(picks=self.picks, unique=True)193 if len(ch_types) > 1:194 raise ValueError(195 "At this point SSD only supports fitting "196 f"single channel types. Your info has {len(ch_types)} types."197 )198 _validate_type(self.cov_method_params, (abc.Mapping, None), "cov_method_params")199 200 def _check_X(self, X, *, y=None, fit=False):201 """Check input data."""202 X = self._check_data(X, y=y, fit=fit, atleast_3d=False)203 n_chan = X.shape[-2]204 if isinstance(self.info, Info) and n_chan != self.info["nchan"]:205 raise ValueError(206 "Info must match the input data."207 f"Found {n_chan} channels but expected {self.info['nchan']}."208 )209 return X210 211 def fit(self, X, y=None):212 """Estimate the SSD decomposition on raw or epoched data.213 214 Parameters215 ----------216 X : array, shape ([n_epochs, ]n_channels, n_times)217 The input data from which to estimate the SSD. Either 2D array218 obtained from continuous data or 3D array obtained from epoched219 data.220 y : None221 Ignored; exists for compatibility with scikit-learn pipelines.222 223 Returns224 -------225 self : instance of SSD226 Returns the modified instance.227 """228 X = self._check_X(X, y=y, fit=True)229 self._validate_params(X)230 if isinstance(self.info, Info):231 info = self.info232 else:233 info = create_info(X.shape[-2], self.sfreq_, ch_types="eeg")234 self.picks_ = _picks_to_idx(info, self.picks, none="data", exclude="bads")235 236 super().fit(X, y)237 238 logger.info("Done.")239 return self240 241 def transform(self, X):242 """Estimate epochs sources given the SSD filters.243 244 Parameters245 ----------246 X : array, shape ([n_epochs, ]n_channels, n_times)247 The input data from which to estimate the SSD. Either 2D array248 obtained from continuous data or 3D array obtained from epoched249 data.250 251 Returns252 -------253 X_ssd : array, shape ([n_epochs, ]n_components, n_times)254 The processed data.255 """256 X = self._check_X(X)257 # For the case where n_epochs dimension is absent.258 if X.ndim == 2:259 X = np.expand_dims(X, axis=0)260 X_aux = X[..., self.picks_, :]261 if self.return_filtered:262 X_aux = filter_data(X_aux, self.sfreq_, **self.filt_params_signal)263 X_ssd = super().transform(X_aux).squeeze()264 265 return X_ssd266 267 def fit_transform(self, X, y=None, **fit_params):268 """Fit SSD to data, then transform it.269 270 Fits transformer to ``X`` and ``y`` with optional parameters ``fit_params``, and271 returns a transformed version of ``X``.272 273 Parameters274 ----------275 X : array, shape ([n_epochs, ]n_channels, n_times)276 The input data from which to estimate the SSD. Either 2D array obtained from277 continuous data or 3D array obtained from epoched data.278 y : None279 Ignored; exists for compatibility with scikit-learn pipelines.280 **fit_params : dict281 Additional fitting parameters passed to the :meth:`mne.decoding.SSD.fit`282 method. Not used for this class.283 284 Returns285 -------286 X_ssd : array, shape ([n_epochs, ]n_components, n_times)287 The processed data.288 """289 # use parent TransformerMixin method but with custom docstring290 return super().fit_transform(X, y=y, **fit_params)291 292 def get_spectral_ratio(self, ssd_sources):293 """Get the spectal signal-to-noise ratio for each spatial filter.294 295 Spectral ratio measure for best n_components selection296 See :footcite:`NikulinEtAl2011`, Eq. (24).297 298 Parameters299 ----------300 ssd_sources : array301 Data projected to SSD space.302 303 Returns304 -------305 spec_ratio : array, shape (n_channels)306 Array with the sprectal ratio value for each component.307 sorter_spec : array, shape (n_channels)308 Array of indices for sorting spec_ratio.309 310 References311 ----------312 .. footbibliography::313 """314 spec_ratio, sorter_spec = _get_spectral_ratio(315 ssd_sources=ssd_sources,316 sfreq=self.sfreq_,317 n_fft=self.n_fft_,318 freqs_signal=self.freqs_signal_,319 freqs_noise=self.freqs_noise_,320 )321 return spec_ratio, sorter_spec322 323 def inverse_transform(self):324 """Not implemented yet."""325 raise NotImplementedError("inverse_transform is not yet available.")326 327 def apply(self, X):328 """Remove selected components from the signal.329 330 This procedure will reconstruct M/EEG signals from which the dynamics331 described by the excluded components is subtracted332 (denoised by low-rank factorization).333 See :footcite:`HaufeEtAl2014b` for more information.334 335 .. note:: Unlike in other classes with an apply method,336 only NumPy arrays are supported (not instances of MNE objects).337 338 Parameters339 ----------340 X : array, shape ([n_epochs, ]n_channels, n_times)341 The input data from which to estimate the SSD. Either 2D array342 obtained from continuous data or 3D array obtained from epoched343 data.344 345 Returns346 -------347 X : array, shape ([n_epochs, ]n_channels, n_times)348 The processed data.349 """350 X_ssd = self.transform(X)351 pick_patterns = self.patterns_[: self.n_components].T352 X = pick_patterns @ X_ssd353 return X354