# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
import numpy as np
from .._fiff.constants import FIFF
from .._fiff.pick import pick_info
from ..baseline import _log_rescale, rescale
from ..epochs import Epochs
from ..event import make_fixed_length_events
from ..evoked import EvokedArray
from ..fixes import _safe_svd
from ..label import BiHemiLabel, Label
from ..parallel import parallel_func
from ..source_estimate import _make_stc
from ..time_frequency.multitaper import (
_compute_mt_params,
_mt_spectra,
_psd_from_mt,
_psd_from_mt_adaptive,
)
from ..time_frequency.tfr import cwt, morlet
from ..utils import ProgressBar, _check_option, _pl, _validate_type, logger, verbose
from .inverse import (
INVERSE_METHODS,
_assemble_kernel,
_check_or_prepare,
_check_ori,
_pick_channels_inverse_operator,
_subject_from_inverse,
combine_xyz,
)
def _restrict_K_to_lbls(labels, K, noise_norm, vertno, pick_ori):
"""Use labels to choose desired sources in the kernel."""
verts_to_use = [[], []]
# create mask for K by compiling original vertices from vertno in labels
for ii in range(len(labels)):
lab = labels[ii]
# handle BiHemi labels; ok so long as no overlap w/ single hemi labels
if lab.hemi == "both":
l_verts = np.intersect1d(vertno[0], lab.lh.vertices)
r_verts = np.intersect1d(vertno[1], lab.rh.vertices) # output sorted
verts_to_use[0] += list(l_verts)
verts_to_use[1] += list(r_verts)
else:
hidx = 0 if lab.hemi == "lh" else 1
verts = np.intersect1d(vertno[hidx], lab.vertices)
verts_to_use[hidx] += list(verts)
# check that we don't have overlapping vertices in our labels
for ii in range(2):
if len(np.unique(verts_to_use[ii])) != len(verts_to_use[ii]):
raise RuntimeError(
"Labels cannot have overlapping vertices. "
"Please select labels with unique vertices "
"and try again."
)
# turn original vertex numbers from vertno into indices for K
K_mask = np.searchsorted(vertno[0], verts_to_use[0])
r_kmask = np.searchsorted(vertno[1], verts_to_use[1]) + len(vertno[0])
K_mask = np.hstack((K_mask, r_kmask))
# record which original vertices are at each index in out_K
hemis = ("lh", "rh")
ki_keys = [
(hemis[hi], verts_to_use[hi][ii])
for hi in range(2)
for ii in range(len(verts_to_use[hi]))
]
ki_vals = list(range(len(K_mask)))
k_idxs = dict(zip(ki_keys, ki_vals))
# mask K, handling the orientation issue
len_allverts = len(vertno[0]) + len(vertno[1])
if len(K) == len_allverts:
assert pick_ori == "normal"
out_K = K[K_mask]
else:
# here, K = [x0, y0, z0, x1, y1, z1 ...]
# we need to drop x, y and z of unused vertices
assert not pick_ori == "normal", pick_ori
assert len(K) == 3 * len_allverts, (len(K), len_allverts)
out_len = len(K_mask) * 3
out_K = K[0:out_len] # get the correct-shaped array
for di in range(3):
K_pick = K[di::3]
out_K[di::3] = K_pick[K_mask] # set correct values for out
out_vertno = verts_to_use
if noise_norm is not None:
out_nn = noise_norm[K_mask]
else:
out_nn = None
return out_K, out_nn, out_vertno, k_idxs
def _prepare_source_params(
inst,
inverse_operator,
label=None,
lambda2=1.0 / 9.0,
method="dSPM",
nave=1,
pca=True,
pick_ori="normal",
prepared=False,
method_params=None,
use_cps=True,
):
"""Prepare inverse operator and params for spectral / TFR analysis."""
inv = _check_or_prepare(
inverse_operator, nave, lambda2, method, method_params, prepared
)
#
# Pick the correct channels from the data
#
sel = _pick_channels_inverse_operator(inst.ch_names, inv)
logger.info("Picked %d channels from the data", len(sel))
logger.info("Computing inverse...")
#
# Simple matrix multiplication followed by combination of the
# three current components
#
# This does all the data transformations to compute the weights for the
# eigenleads
#
# K shape: (3 x n_sources, n_channels) or (n_sources, n_channels)
# noise_norm shape: (n_sources, 1)
# vertno: [lh_verts, rh_verts]
k_idxs = None
if not isinstance(label, Label | BiHemiLabel):
whole_K, whole_noise_norm, whole_vertno, _ = _assemble_kernel(
inv, None, method, pick_ori, use_cps=use_cps
)
if isinstance(label, list):
K, noise_norm, vertno, k_idxs = _restrict_K_to_lbls(
label, whole_K, whole_noise_norm, whole_vertno, pick_ori
)
else:
assert not label
K, noise_norm, vertno = whole_K, whole_noise_norm, whole_vertno
elif isinstance(label, Label | BiHemiLabel):
K, noise_norm, vertno, _ = _assemble_kernel(
inv, label, method, pick_ori, use_cps=use_cps
)
if pca:
U, s, Vh = _safe_svd(K, full_matrices=False)
rank = np.sum(s > 1e-8 * s[0])
K = s[:rank] * U[:, :rank]
Vh = Vh[:rank]
logger.info("Reducing data rank %d -> %d", len(s), rank)
else:
Vh = None
is_free_ori = inverse_operator["source_ori"] == FIFF.FIFFV_MNE_FREE_ORI
return K, sel, Vh, vertno, is_free_ori, noise_norm, k_idxs
@verbose
def source_band_induced_power(
epochs,
inverse_operator,
bands,
label=None,
lambda2=1.0 / 9.0,
method="dSPM",
nave=1,
n_cycles=5,
df=1,
use_fft=False,
decim=1,
baseline=None,
baseline_mode="logratio",
pca=True,
n_jobs=None,
prepared=False,
method_params=None,
use_cps=True,
*,
verbose=None,
):
"""Compute source space induced power in given frequency bands.
Parameters
----------
epochs : instance of Epochs
The epochs.
inverse_operator : instance of InverseOperator
The inverse operator.
bands : dict
Example : bands = dict(alpha=[8, 9]).
label : Label | list of Label
Restricts the source estimates to a given label or list of labels. If
labels are provided in a list, power will be averaged over vertices.
lambda2 : float
The regularization parameter of the minimum norm.
method : "MNE" | "dSPM" | "sLORETA" | "eLORETA"
Use minimum norm, dSPM (default), sLORETA, or eLORETA.
nave : int
The number of averages used to scale the noise covariance matrix.
n_cycles : float | array of float
Number of cycles. Fixed number or one per frequency.
df : float
Delta frequency within bands.
use_fft : bool
Do convolutions in time or frequency domain with FFT.
decim : int
Temporal decimation factor.
baseline : None (default) or tuple, shape (2,)
The time interval to apply baseline correction. If None do not apply
it. If baseline is (a, b) the interval is between "a (s)" and "b (s)".
If a is None the beginning of the data is used and if b is None then b
is set to the end of the interval. If baseline is equal to (None, None)
all the time interval is used.
baseline_mode : 'mean' | 'ratio' | 'logratio' | 'percent' | 'zscore' | 'zlogratio'
Perform baseline correction by
- subtracting the mean of baseline values ('mean')
- dividing by the mean of baseline values ('ratio')
- dividing by the mean of baseline values and taking the log
('logratio')
- subtracting the mean of baseline values followed by dividing by
the mean of baseline values ('percent')
- subtracting the mean of baseline values and dividing by the
standard deviation of baseline values ('zscore')
- dividing by the mean of baseline values, taking the log, and
dividing by the standard deviation of log baseline values
('zlogratio')
pca : bool
If True, the true dimension of data is estimated before running
the time-frequency transforms. It reduces the computation times
e.g. with a dataset that was maxfiltered (true dim is 64).
%(n_jobs)s
prepared : bool
If True, do not call :func:`prepare_inverse_operator`.
method_params : dict | None
Additional options for eLORETA. See Notes of :func:`apply_inverse`.
.. versionadded:: 0.16
%(use_cps_restricted)s
.. versionadded:: 0.20
%(verbose)s
Returns
-------
stcs : dict of SourceEstimate (or VolSourceEstimate)
The estimated source space induced power estimates in shape
(n_vertices, n_frequencies, n_samples) if label=None or label=label.
For lists of one or more labels, the induced power estimate has shape
(n_labels, n_frequencies, n_samples).
""" # noqa: E501
_check_option("method", method, INVERSE_METHODS)
freqs = np.concatenate(
[np.arange(band[0], band[1] + df / 2.0, df) for _, band in bands.items()]
)
powers, _, vertno = _source_induced_power(
epochs,
inverse_operator,
freqs,
label=label,
lambda2=lambda2,
method=method,
nave=nave,
n_cycles=n_cycles,
decim=decim,
use_fft=use_fft,
pca=pca,
n_jobs=n_jobs,
with_plv=False,
prepared=prepared,
method_params=method_params,
use_cps=use_cps,
)
Fs = epochs.info["sfreq"] # sampling in Hz
stcs = dict()
subject = _subject_from_inverse(inverse_operator)
_log_rescale(baseline, baseline_mode) # for early failure
for name, band in bands.items():
idx = [k for k, f in enumerate(freqs) if band[0] <= f <= band[1]]
# average power in band + mean over epochs
power = np.mean(powers[:, idx, :], axis=1)
# Run baseline correction
power = rescale(
power,
epochs.times[::decim],
baseline,
baseline_mode,
copy=False,
verbose=False,
)
tmin = epochs.times[0]
tstep = float(decim) / Fs
stc = _make_stc(
power,
vertices=vertno,
tmin=tmin,
tstep=tstep,
subject=subject,
src_type=inverse_operator["src"].kind,
)
stcs[name] = stc
logger.info("[done]")
return stcs
def _prepare_tfr(data, decim, pick_ori, Ws, K, source_ori):
"""Prepare TFR source localization."""
n_times = data[:, :, ::decim].shape[2]
n_freqs = len(Ws)
n_sources = K.shape[0]
is_free_ori = False
if source_ori == FIFF.FIFFV_MNE_FREE_ORI and pick_ori is None:
is_free_ori = True
n_sources //= 3
shape = (n_sources, n_freqs, n_times)
return shape, is_free_ori
@verbose
def _compute_pow_plv(
data,
K,
sel,
Ws,
source_ori,
use_fft,
Vh,
with_power,
with_plv,
pick_ori,
decim,
noise_norm=None,
verbose=None,
):
"""Aux function for induced power and PLV."""
shape, is_free_ori = _prepare_tfr(data, decim, pick_ori, Ws, K, source_ori)
power = np.zeros(shape, dtype=np.float64) # power or raw TFR
# phase lock
plv = np.zeros(shape, dtype=np.complex128) if with_plv else None
for epoch in data:
epoch = epoch[sel] # keep only selected channels
if Vh is not None:
epoch = np.dot(Vh, epoch) # reducing data rank
power_e, plv_e = _single_epoch_tfr(
data=epoch,
is_free_ori=is_free_ori,
K=K,
Ws=Ws,
use_fft=use_fft,
decim=decim,
shape=shape,
with_plv=with_plv,
with_power=with_power,
)
power += power_e
if with_plv:
plv += plv_e
if noise_norm is not None:
power *= noise_norm[:, :, np.newaxis] ** 2
return power, plv
def _single_epoch_tfr(
data, is_free_ori, K, Ws, use_fft, decim, shape, with_plv, with_power
):
"""Compute single trial TFRs, either ITC, power or raw TFR."""
tfr_e = np.zeros(shape, dtype=np.float64) # power or raw TFR
# phase lock
plv_e = np.zeros(shape, dtype=np.complex128) if with_plv else None
n_sources, _, n_times = shape
for f, w in enumerate(Ws):
tfr_ = cwt(data, [w], use_fft=use_fft, decim=decim)
tfr_ = np.asfortranarray(tfr_.reshape(len(data), -1))
# phase lock and power at freq f
if with_plv:
plv_f = np.zeros((n_sources, n_times), dtype=np.complex128)
tfr_f = np.zeros((n_sources, n_times), dtype=np.float64)
for k, t in enumerate([np.real(tfr_), np.imag(tfr_)]):
sol = np.dot(K, t)
sol_pick_normal = sol
if is_free_ori:
sol_pick_normal = sol[2::3]
if with_plv:
if k == 0: # real
plv_f += sol_pick_normal
else: # imag
plv_f += 1j * sol_pick_normal
if is_free_ori:
logger.debug("combining the current components...")
sol = combine_xyz(sol, square=with_power)
elif with_power:
sol *= sol
tfr_f += sol
del sol
tfr_e[:, f, :] += tfr_f
del tfr_f
if with_plv:
plv_f /= np.abs(plv_f)
plv_e[:, f, :] += plv_f
del plv_f
return tfr_e, plv_e
def _get_label_power(power, labels, vertno, k_idxs):
"""Average power across vertices in labels."""
(_, ps1, ps2) = power.shape
# construct out array with correct shape
out_power = np.zeros(shape=(len(labels), ps1, ps2))
# for each label, compile list of vertices we want
for li in np.arange(len(labels)):
lab = labels[li]
hemis = ("lh", "rh")
all_vnums = [[], []]
if lab.hemi == "both":
all_vnums[0] = np.intersect1d(lab.lh.vertices, vertno[0])
all_vnums[1] = np.intersect1d(lab.rh.vertices, vertno[1])
else:
assert lab.hemi == "lh" or lab.hemi == "rh"
h_id = 0 if lab.hemi == "lh" else 1
all_vnums[h_id] = np.intersect1d(vertno[h_id], lab.vertices)
verts = [(hemis[hi], vn) for hi in range(2) for vn in all_vnums[hi]]
# restrict power to relevant vertices in label
lab_mask = np.array([False] * len(power))
for vert in verts:
lab_mask[k_idxs[vert]] = True # k_idxs[vert] gives power row index
lab_power = power[lab_mask] # only pass through rows we want
assert lab_power.shape == (len(verts), ps1, ps2)
# set correct out values for label
out_power[li, :, :] = np.mean(lab_power, axis=0)
assert out_power.shape == (len(labels), ps1, ps2)
return out_power
@verbose
def _source_induced_power(
epochs,
inverse_operator,
freqs,
label=None,
lambda2=1.0 / 9.0,
method="dSPM",
nave=1,
n_cycles=5,
decim=1,
use_fft=False,
pca=True,
pick_ori="normal",
n_jobs=None,
with_plv=True,
zero_mean=False,
prepared=False,
method_params=None,
use_cps=True,
verbose=None,
):
"""Aux function for source induced power."""
if label:
_validate_type(
label,
types=(Label, BiHemiLabel, list, tuple, None),
type_name=("Label or BiHemiLabel", "list of labels", "None"),
)
if isinstance(label, list | tuple):
for item in label:
_validate_type(
item,
types=(Label, BiHemiLabel),
type_name=("Label or BiHemiLabel"),
)
if len(label) > 1 and with_plv:
raise RuntimeError(
"Phase-locking value cannot be calculated "
"when averaging induced power within "
"labels. Please set `with_plv` to False, pass a "
"single `label=label`, or set `label=None`."
)
epochs_data = epochs.get_data(copy=False)
K, sel, Vh, vertno, is_free_ori, noise_norm, k_id = _prepare_source_params(
inst=epochs,
inverse_operator=inverse_operator,
label=label,
lambda2=lambda2,
method=method,
nave=nave,
pca=pca,
pick_ori=pick_ori,
prepared=prepared,
method_params=method_params,
use_cps=use_cps,
)
inv = inverse_operator
parallel, my_compute_source_tfrs, n_jobs = parallel_func(
_compute_pow_plv, n_jobs, max_jobs=len(epochs_data)
)
Fs = epochs.info["sfreq"] # sampling in Hz
logger.info("Computing source power ...")
Ws = morlet(Fs, freqs, n_cycles=n_cycles, zero_mean=zero_mean)
out = parallel(
my_compute_source_tfrs(
data=data,
K=K,
sel=sel,
Ws=Ws,
source_ori=inv["source_ori"],
use_fft=use_fft,
Vh=Vh,
with_plv=with_plv,
with_power=True,
pick_ori=pick_ori,
decim=decim,
noise_norm=noise_norm,
)
for data in np.array_split(epochs_data, n_jobs)
)
power = sum(o[0] for o in out) # power shape: (n_verts, n_freqs, n_samps)
power /= len(epochs_data) # average power over epochs
if isinstance(label, Label | BiHemiLabel):
logger.info(
f"Outputting power for {len(power)} vertices in label {label.name}."
)
elif isinstance(label, list):
power = _get_label_power(power, label, vertno, k_id)
logger.info(
"Averaging induced power across vertices within labels "
f"for {len(label)} label{_pl(label)}."
)
else:
assert not label
if with_plv:
plv = sum(o[1] for o in out)
plv = np.abs(plv)
plv /= len(epochs_data) # average power over epochs
else:
plv = None
return power, plv, vertno
@verbose
def source_induced_power(
epochs,
inverse_operator,
freqs,
label=None,
lambda2=1.0 / 9.0,
method="dSPM",
nave=1,
n_cycles=5,
decim=1,
use_fft=False,
pick_ori=None,
baseline=None,
baseline_mode="logratio",
pca=True,
n_jobs=None,
*,
return_plv=True,
zero_mean=False,
prepared=False,
method_params=None,
use_cps=True,
verbose=None,
):
"""Compute induced power and phase lock.
Computation can optionally be restricted in a label.
Parameters
----------
epochs : instance of Epochs
The epochs.
inverse_operator : instance of InverseOperator
The inverse operator.
freqs : array
Array of frequencies of interest.
label : Label | list of Label
Restricts the source estimates to a given label or list of labels. If
labels are provided in a list, power will be averaged over vertices within each
label.
lambda2 : float
The regularization parameter of the minimum norm.
method : "MNE" | "dSPM" | "sLORETA" | "eLORETA"
Use minimum norm, dSPM (default), sLORETA, or eLORETA.
nave : int
The number of averages used to scale the noise covariance matrix.
n_cycles : float | array of float
Number of cycles. Fixed number or one per frequency.
decim : int
Temporal decimation factor.
use_fft : bool
Do convolutions in time or frequency domain with FFT.
pick_ori : None | "normal"
If "normal", rather than pooling the orientations by taking the norm,
only the radial component is kept. This is only implemented
when working with loose orientations.
baseline : None (default) or tuple of length 2
The time interval to apply baseline correction.
If None do not apply it. If baseline is (a, b)
the interval is between "a (s)" and "b (s)".
If a is None the beginning of the data is used
and if b is None then b is set to the end of the interval.
If baseline is equal to (None, None) all the time
interval is used.
baseline_mode : 'mean' | 'ratio' | 'logratio' | 'percent' | 'zscore' | 'zlogratio'
Perform baseline correction by
- subtracting the mean of baseline values ('mean')
- dividing by the mean of baseline values ('ratio')
- dividing by the mean of baseline values and taking the log
('logratio')
- subtracting the mean of baseline values followed by dividing by
the mean of baseline values ('percent')
- subtracting the mean of baseline values and dividing by the
standard deviation of baseline values ('zscore')
- dividing by the mean of baseline values, taking the log, and
dividing by the standard deviation of log baseline values
('zlogratio')
pca : bool
If True, the true dimension of data is estimated before running
the time-frequency transforms. It reduces the computation times
e.g. with a dataset that was maxfiltered (true dim is 64).
%(n_jobs)s
return_plv : bool
If True, return the phase-locking value array. Else, only return power.
.. versionadded:: 1.6
zero_mean : bool
Make sure the wavelets are zero mean.
prepared : bool
If True, do not call :func:`prepare_inverse_operator`.
method_params : dict | None
Additional options for eLORETA. See Notes of :func:`apply_inverse`.
%(use_cps_restricted)s
.. versionadded:: 0.20
%(verbose)s
Returns
-------
power : array
The induced power array with shape (n_sources, n_freqs, n_samples) if
label=None or label=label. For lists of one or more labels, the induced
power estimate has shape (n_labels, n_frequencies, n_samples).
plv : array
The phase-locking value array with shape (n_sources, n_freqs,
n_samples). Only returned if ``return_plv=True``.
""" # noqa: E501
_check_option("method", method, INVERSE_METHODS)
_check_ori(pick_ori, inverse_operator["source_ori"], inverse_operator["src"])
power, plv, vertno = _source_induced_power(
epochs,
inverse_operator,
freqs,
label=label,
lambda2=lambda2,
method=method,
nave=nave,
n_cycles=n_cycles,
decim=decim,
use_fft=use_fft,
pick_ori=pick_ori,
pca=pca,
n_jobs=n_jobs,
with_plv=return_plv,
method_params=method_params,
zero_mean=zero_mean,
prepared=prepared,
use_cps=use_cps,
)
# Run baseline correction
power = rescale(power, epochs.times[::decim], baseline, baseline_mode, copy=False)
outs = (power, plv) if return_plv else power
return outs
@verbose
def compute_source_psd(
raw,
inverse_operator,
lambda2=1.0 / 9.0,
method="dSPM",
tmin=0.0,
tmax=None,
fmin=0.0,
fmax=200.0,
n_fft=2048,
overlap=0.5,
pick_ori=None,
label=None,
nave=1,
pca=True,
prepared=False,
method_params=None,
inv_split=None,
bandwidth="hann",
adaptive=False,
low_bias=False,
n_jobs=None,
return_sensor=False,
dB=False,
*,
verbose=None,
):
"""Compute source power spectral density (PSD).
Parameters
----------
raw : instance of Raw
The raw data.
inverse_operator : instance of InverseOperator
The inverse operator.
lambda2 : float
The regularization parameter.
method : "MNE" | "dSPM" | "sLORETA"
Use minimum norm, dSPM (default), sLORETA, or eLORETA.
tmin : float
The beginning of the time interval of interest (in seconds).
Use 0. for the beginning of the file.
tmax : float | None
The end of the time interval of interest (in seconds). If None
stop at the end of the file.
fmin : float
The lower frequency of interest.
fmax : float
The upper frequency of interest.
n_fft : int
Window size for the FFT. Should be a power of 2.
overlap : float
The overlap fraction between windows. Should be between 0 and 1.
0 means no overlap.
pick_ori : None | "normal"
If "normal", rather than pooling the orientations by taking the norm,
only the radial component is kept. This is only implemented
when working with loose orientations.
label : Label
Restricts the source estimates to a given label.
nave : int
The number of averages used to scale the noise covariance matrix.
pca : bool
If True, the true dimension of data is estimated before running
the time-frequency transforms. It reduces the computation times
e.g. with a dataset that was maxfiltered (true dim is 64).
prepared : bool
If True, do not call :func:`prepare_inverse_operator`.
method_params : dict | None
Additional options for eLORETA. See Notes of :func:`apply_inverse`.
.. versionadded:: 0.16
inv_split : int or None
Split inverse operator into inv_split parts in order to save memory.
.. versionadded:: 0.17
bandwidth : float | str
The bandwidth of the multi taper windowing function in Hz.
Can also be a string (e.g., 'hann') to use a single window.
For backward compatibility, the default is 'hann'.
.. versionadded:: 0.17
adaptive : bool
Use adaptive weights to combine the tapered spectra into PSD
(slow, use n_jobs >> 1 to speed up computation).
.. versionadded:: 0.17
low_bias : bool
Only use tapers with more than 90%% spectral concentration within
bandwidth.
.. versionadded:: 0.17
%(n_jobs)s
It is only used if adaptive=True.
.. versionadded:: 0.17
return_sensor : bool
If True, return the sensor PSDs as an EvokedArray.
.. versionadded:: 0.17
dB : bool
If True (default False), return output it decibels.
.. versionadded:: 0.17
%(verbose)s
Returns
-------
stc_psd : instance of SourceEstimate | VolSourceEstimate
The PSD of each of the sources.
sensor_psd : instance of EvokedArray
The PSD of each sensor. Only returned if ``return_sensor`` is True.
See Also
--------
compute_source_psd_epochs
Notes
-----
Each window is multiplied by a window before processing, so
using a non-zero overlap is recommended.
This function is different from :func:`compute_source_psd_epochs` in that:
1. ``bandwidth='hann'`` by default, skipping multitaper estimation
2. For convenience it wraps
:func:`mne.make_fixed_length_events` and :class:`mne.Epochs`.
Otherwise the two should produce identical results.
"""
tmin = 0.0 if tmin is None else float(tmin)
overlap = float(overlap)
if not 0 <= overlap < 1:
raise ValueError(f"Overlap must be at least 0 and less than 1, got {overlap}")
n_fft = int(n_fft)
duration = ((1.0 - overlap) * n_fft) / raw.info["sfreq"]
events = make_fixed_length_events(raw, 1, tmin, tmax, duration)
epochs = Epochs(raw, events, 1, 0, (n_fft - 1) / raw.info["sfreq"], baseline=None)
out = compute_source_psd_epochs(
epochs,
inverse_operator,
lambda2,
method,
fmin,
fmax,
pick_ori,
label,
nave,
pca,
inv_split,
bandwidth,
adaptive,
low_bias,
True,
n_jobs,
prepared,
method_params,
return_sensor=True,
)
source_data = 0.0
sensor_data = 0.0
count = 0
for stc, evoked in out:
source_data += stc.data
sensor_data += evoked.data
count += 1
assert count > 0 # should be guaranteed by make_fixed_length_events
sensor_data /= count
source_data /= count
if dB:
np.log10(sensor_data, out=sensor_data)
sensor_data *= 10.0
np.log10(source_data, out=source_data)
source_data *= 10.0
evoked.data = sensor_data
evoked.nave = count
stc.data = source_data
out = stc
if return_sensor:
out = (out, evoked)
return out
def _compute_source_psd_epochs(
epochs,
inverse_operator,
lambda2=1.0 / 9.0,
method="dSPM",
fmin=0.0,
fmax=200.0,
pick_ori=None,
label=None,
nave=1,
pca=True,
inv_split=None,
bandwidth=4.0,
adaptive=False,
low_bias=True,
n_jobs=None,
prepared=False,
method_params=None,
return_sensor=False,
use_cps=True,
):
"""Generate compute_source_psd_epochs."""
logger.info(f"Considering frequencies {fmin} ... {fmax} Hz")
if label:
# TODO: add multi-label support
# since `_prepare_source_params` can handle a list of labels now,
# multi-label support should be within reach for psd calc as well
_validate_type(
label,
types=(Label, BiHemiLabel, None),
type_name=("Label or BiHemiLabel", "None"),
)
K, sel, Vh, vertno, is_free_ori, noise_norm, _ = _prepare_source_params(
inst=epochs,
inverse_operator=inverse_operator,
label=label,
lambda2=lambda2,
method=method,
nave=nave,
pca=pca,
pick_ori=pick_ori,
prepared=prepared,
method_params=method_params,
use_cps=use_cps,
)
# Simplify code with a tiny (rel. to other computations) penalty for eye
# mult
Vh = np.eye(K.shape[1]) if Vh is None else Vh
# split the inverse operator
if inv_split is not None:
K_split = np.array_split(K, inv_split)
else:
K_split = [K]
# compute DPSS windows
n_times = len(epochs.times)
sfreq = epochs.info["sfreq"]
dpss, eigvals, adaptive = _compute_mt_params(
n_times, sfreq, bandwidth, low_bias, adaptive, verbose=False
)
n_tapers = len(dpss)
try:
n_epochs = len(epochs)
except RuntimeError:
n_epochs = len(epochs.events)
extra = f"on at most {n_epochs} epochs"
else:
extra = f"on {n_epochs} epochs"
if isinstance(bandwidth, str):
bandwidth = f"{bandwidth} windowing"
else:
bandwidth = f"{n_tapers} tapers with bandwidth {bandwidth:0.1f} Hz"
logger.info(f"Using {bandwidth} {extra}")
if adaptive:
parallel, my_psd_from_mt_adaptive, n_jobs = parallel_func(
_psd_from_mt_adaptive, n_jobs
)
else:
weights = np.sqrt(eigvals)[np.newaxis, :, np.newaxis]
subject = _subject_from_inverse(inverse_operator)
iter_epochs = ProgressBar(epochs, max_value=n_epochs)
evoked_info = pick_info(epochs.info, sel, verbose=False)
for k, e in enumerate(iter_epochs):
data = np.dot(Vh, e[sel]) # reducing data rank
# compute tapered spectra in sensor space
x_mt, freqs = _mt_spectra(data, dpss, sfreq)
if k == 0:
freq_mask = (freqs >= fmin) & (freqs <= fmax)
fstep = np.mean(np.diff(freqs))
with evoked_info._unlock():
evoked_info["sfreq"] = 1.0 / fstep
freqs = freqs[freq_mask]
# sensor space PSD
x_mt_sensor = np.empty(
(len(sel), x_mt.shape[1], x_mt.shape[2]), dtype=x_mt.dtype
)
for i in range(n_tapers):
x_mt_sensor[:, i, :] = np.dot(Vh.T, x_mt[:, i, :])
if adaptive:
out = parallel(
my_psd_from_mt_adaptive(x, eigvals, freq_mask)
for x in np.array_split(x_mt_sensor, min(n_jobs, len(x_mt_sensor)))
)
sensor_psd = np.concatenate(out)
else:
x_mt_sensor = x_mt_sensor[:, :, freq_mask]
sensor_psd = _psd_from_mt(x_mt_sensor, weights)
# allocate space for output
psd = np.empty((K.shape[0], np.sum(freq_mask)))
# Optionally, we split the inverse operator into parts to save memory.
# Without splitting the tapered spectra in source space have size
# (n_vertices x n_tapers x n_times / 2)
pos = 0
for K_part in K_split:
# allocate space for tapered spectra in source space
x_mt_src = np.empty(
(K_part.shape[0], x_mt.shape[1], x_mt.shape[2]), dtype=x_mt.dtype
)
# apply inverse to each taper (faster than equiv einsum)
for i in range(n_tapers):
x_mt_src[:, i, :] = np.dot(K_part, x_mt[:, i, :])
# compute the psd
if adaptive:
out = parallel(
my_psd_from_mt_adaptive(x, eigvals, freq_mask)
for x in np.array_split(x_mt_src, min(n_jobs, len(x_mt_src)))
)
this_psd = np.concatenate(out)
else:
x_mt_src = x_mt_src[:, :, freq_mask]
this_psd = _psd_from_mt(x_mt_src, weights)
psd[pos : pos + K_part.shape[0], :] = this_psd
pos += K_part.shape[0]
# combine orientations
if is_free_ori and pick_ori is None:
psd = combine_xyz(psd, square=False)
if noise_norm is not None:
psd *= noise_norm**2
out = _make_stc(
psd,
tmin=freqs[0],
tstep=fstep,
vertices=vertno,
subject=subject,
src_type=inverse_operator["src"].kind,
)
if return_sensor:
comment = f"Epoch {k} PSD"
out = (
out,
EvokedArray(sensor_psd, evoked_info.copy(), freqs[0], comment, nave),
)
# we return a generator object for "stream processing"
yield out
iter_epochs.update(n_epochs) # in case some were skipped
@verbose
def compute_source_psd_epochs(
epochs,
inverse_operator,
lambda2=1.0 / 9.0,
method="dSPM",
fmin=0.0,
fmax=200.0,
pick_ori=None,
label=None,
nave=1,
pca=True,
inv_split=None,
bandwidth=4.0,
adaptive=False,
low_bias=True,
return_generator=False,
n_jobs=None,
prepared=False,
method_params=None,
return_sensor=False,
use_cps=True,
verbose=None,
):
"""Compute source power spectral density (PSD) from Epochs.
This uses the multi-taper method to compute the PSD for each epoch.
Parameters
----------
epochs : instance of Epochs
The raw data.
inverse_operator : instance of InverseOperator
The inverse operator.
lambda2 : float
The regularization parameter.
method : "MNE" | "dSPM" | "sLORETA" | "eLORETA"
Use minimum norm, dSPM (default), sLORETA, or eLORETA.
fmin : float
The lower frequency of interest.
fmax : float
The upper frequency of interest.
pick_ori : None | "normal"
If "normal", rather than pooling the orientations by taking the norm,
only the radial component is kept. This is only implemented
when working with loose orientations.
label : Label
Restricts the source estimates to a given label.
nave : int
The number of averages used to scale the noise covariance matrix.
pca : bool
If True, the true dimension of data is estimated before running
the time-frequency transforms. It reduces the computation times
e.g. with a dataset that was maxfiltered (true dim is 64).
inv_split : int or None
Split inverse operator into inv_split parts in order to save memory.
bandwidth : float | str
The bandwidth of the multi taper windowing function in Hz.
Can also be a string (e.g., 'hann') to use a single window.
adaptive : bool
Use adaptive weights to combine the tapered spectra into PSD
(slow, use n_jobs >> 1 to speed up computation).
low_bias : bool
Only use tapers with more than 90%% spectral concentration within
bandwidth.
return_generator : bool
Return a generator object instead of a list. This allows iterating
over the stcs without having to keep them all in memory.
%(n_jobs)s
It is only used if adaptive=True.
prepared : bool
If True, do not call :func:`prepare_inverse_operator`.
method_params : dict | None
Additional options for eLORETA. See Notes of :func:`apply_inverse`.
.. versionadded:: 0.16
return_sensor : bool
If True, also return the sensor PSD for each epoch as an EvokedArray.
.. versionadded:: 0.17
%(use_cps_restricted)s
.. versionadded:: 0.20
%(verbose)s
Returns
-------
out : list (or generator object)
A list (or generator) for the source space PSD (and optionally the
sensor PSD) for each epoch.
See Also
--------
compute_source_psd
"""
# use an auxiliary function so we can either return a generator or a list
stcs_gen = _compute_source_psd_epochs(
epochs,
inverse_operator,
lambda2=lambda2,
method=method,
fmin=fmin,
fmax=fmax,
pick_ori=pick_ori,
label=label,
nave=nave,
pca=pca,
inv_split=inv_split,
bandwidth=bandwidth,
adaptive=adaptive,
low_bias=low_bias,
n_jobs=n_jobs,
prepared=prepared,
method_params=method_params,
return_sensor=return_sensor,
use_cps=use_cps,
)
if return_generator:
# return generator object
return stcs_gen
else:
# return a list
stcs = list()
for stc in stcs_gen:
stcs.append(stc)
return stcs
Datasets
Models
