Note
Go to the end to download the full example code.
Compute source power spectral density (PSD) in a label#
Returns an STC file containing the PSD (in dB) of each of the sources within a label.
# Authors: Alexandre Gramfort <alexandre.gramfort@inria.fr>
#
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
import matplotlib.pyplot as plt
import mne
from mne import io
from mne.datasets import sample
from mne.minimum_norm import compute_source_psd, read_inverse_operator
print(__doc__)
Set parameters
data_path = sample.data_path()
meg_path = data_path / "MEG" / "sample"
raw_fname = meg_path / "sample_audvis_raw.fif"
fname_inv = meg_path / "sample_audvis-meg-oct-6-meg-inv.fif"
fname_label = meg_path / "labels" / "Aud-lh.label"
# Setup for reading the raw data
raw = io.read_raw_fif(raw_fname, verbose=False)
events = mne.find_events(raw, stim_channel="STI 014")
inverse_operator = read_inverse_operator(fname_inv)
raw.info["bads"] = ["MEG 2443", "EEG 053"]
# picks MEG gradiometers
picks = mne.pick_types(
raw.info, meg=True, eeg=False, eog=True, stim=False, exclude="bads"
)
tmin, tmax = 0, 120 # use the first 120s of data
fmin, fmax = 4, 100 # look at frequencies between 4 and 100Hz
n_fft = 2048 # the FFT size (n_fft). Ideally a power of 2
label = mne.read_label(fname_label)
stc = compute_source_psd(
raw,
inverse_operator,
lambda2=1.0 / 9.0,
method="dSPM",
tmin=tmin,
tmax=tmax,
fmin=fmin,
fmax=fmax,
pick_ori="normal",
n_fft=n_fft,
label=label,
dB=True,
)
stc.save("psd_dSPM", overwrite=True)
320 events found on stim channel STI 014
Event IDs: [ 1 2 3 4 5 32]
Reading inverse operator decomposition from /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif...
Reading inverse operator info...
[done]
Reading inverse operator decomposition...
[done]
305 x 305 full covariance (kind = 1) found.
Read a total of 4 projection items:
PCA-v1 (1 x 102) active
PCA-v2 (1 x 102) active
PCA-v3 (1 x 102) active
Average EEG reference (1 x 60) active
Noise covariance matrix read.
22494 x 22494 diagonal covariance (kind = 2) found.
Source covariance matrix read.
22494 x 22494 diagonal covariance (kind = 6) found.
Orientation priors read.
22494 x 22494 diagonal covariance (kind = 5) found.
Depth priors read.
Did not find the desired covariance matrix (kind = 3)
Reading a source space...
Computing patch statistics...
Patch information added...
Distance information added...
[done]
Reading a source space...
Computing patch statistics...
Patch information added...
Distance information added...
[done]
2 source spaces read
Read a total of 4 projection items:
PCA-v1 (1 x 102) active
PCA-v2 (1 x 102) active
PCA-v3 (1 x 102) active
Average EEG reference (1 x 60) active
Source spaces transformed to the inverse solution coordinate frame
Not setting metadata
70 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 3)
3 projection items activated
Considering frequencies 4 ... 100 Hz
Preparing the inverse operator for use...
Scaled noise and source covariance from nave = 1 to nave = 1
Created the regularized inverter
Created an SSP operator (subspace dimension = 3)
Created the whitener using a noise covariance matrix with rank 302 (3 small eigenvalues omitted)
Computing noise-normalization factors (dSPM)...
[done]
Picked 305 channels from the data
Computing inverse...
Eigenleads need to be weighted ...
Reducing data rank 33 -> 33
Using hann windowing on at most 70 epochs
0%| | : 0/70 [00:00<?, ?it/s]
1%|▏ | : 1/70 [00:00<00:07, 9.03it/s]
3%|▎ | : 2/70 [00:00<00:04, 13.79it/s]
4%|▍ | : 3/70 [00:00<00:03, 18.01it/s]
6%|▌ | : 4/70 [00:00<00:03, 21.44it/s]
7%|▋ | : 5/70 [00:00<00:02, 24.14it/s]
9%|▊ | : 6/70 [00:00<00:02, 25.57it/s]
10%|█ | : 7/70 [00:00<00:02, 27.07it/s]
11%|█▏ | : 8/70 [00:00<00:02, 28.15it/s]
13%|█▎ | : 9/70 [00:00<00:02, 29.05it/s]
14%|█▍ | : 10/70 [00:00<00:02, 30.00it/s]
16%|█▌ | : 11/70 [00:00<00:01, 31.23it/s]
17%|█▋ | : 12/70 [00:00<00:01, 30.84it/s]
19%|█▊ | : 13/70 [00:00<00:01, 31.85it/s]
20%|██ | : 14/70 [00:00<00:01, 32.92it/s]
21%|██▏ | : 15/70 [00:00<00:01, 33.89it/s]
23%|██▎ | : 16/70 [00:00<00:01, 34.58it/s]
24%|██▍ | : 17/70 [00:00<00:01, 35.04it/s]
26%|██▌ | : 18/70 [00:00<00:01, 35.10it/s]
27%|██▋ | : 19/70 [00:00<00:01, 35.86it/s]
29%|██▊ | : 20/70 [00:00<00:01, 36.62it/s]
30%|███ | : 21/70 [00:00<00:01, 37.31it/s]
31%|███▏ | : 22/70 [00:00<00:01, 37.89it/s]
33%|███▎ | : 23/70 [00:00<00:01, 37.72it/s]
34%|███▍ | : 24/70 [00:00<00:01, 38.19it/s]
36%|███▌ | : 25/70 [00:00<00:01, 38.55it/s]
37%|███▋ | : 26/70 [00:00<00:01, 38.91it/s]
39%|███▊ | : 27/70 [00:00<00:01, 39.34it/s]
40%|████ | : 28/70 [00:00<00:01, 39.61it/s]
41%|████▏ | : 29/70 [00:00<00:01, 39.19it/s]
43%|████▎ | : 30/70 [00:00<00:01, 39.45it/s]
44%|████▍ | : 31/70 [00:00<00:00, 39.87it/s]
46%|████▌ | : 32/70 [00:00<00:00, 40.31it/s]
47%|████▋ | : 33/70 [00:00<00:00, 40.65it/s]
49%|████▊ | : 34/70 [00:00<00:00, 40.43it/s]
50%|█████ | : 35/70 [00:00<00:00, 40.76it/s]
51%|█████▏ | : 36/70 [00:00<00:00, 41.12it/s]
53%|█████▎ | : 37/70 [00:00<00:00, 41.55it/s]
54%|█████▍ | : 38/70 [00:00<00:00, 41.80it/s]
56%|█████▌ | : 39/70 [00:01<00:00, 41.96it/s]
57%|█████▋ | : 40/70 [00:01<00:00, 41.53it/s]
59%|█████▊ | : 41/70 [00:01<00:00, 41.83it/s]
60%|██████ | : 42/70 [00:01<00:00, 42.06it/s]
61%|██████▏ | : 43/70 [00:01<00:00, 42.24it/s]
63%|██████▎ | : 44/70 [00:01<00:00, 42.32it/s]
64%|██████▍ | : 45/70 [00:01<00:00, 42.39it/s]
66%|██████▌ | : 46/70 [00:01<00:00, 41.98it/s]
67%|██████▋ | : 47/70 [00:01<00:00, 42.19it/s]
69%|██████▊ | : 48/70 [00:01<00:00, 42.36it/s]
70%|███████ | : 49/70 [00:01<00:00, 42.53it/s]
71%|███████▏ | : 50/70 [00:01<00:00, 42.63it/s]
73%|███████▎ | : 51/70 [00:01<00:00, 42.08it/s]
74%|███████▍ | : 52/70 [00:01<00:00, 42.31it/s]
76%|███████▌ | : 53/70 [00:01<00:00, 42.32it/s]
77%|███████▋ | : 54/70 [00:01<00:00, 42.52it/s]
79%|███████▊ | : 55/70 [00:01<00:00, 42.71it/s]
80%|████████ | : 56/70 [00:01<00:00, 42.85it/s]
81%|████████▏ | : 57/70 [00:01<00:00, 42.53it/s]
83%|████████▎ | : 58/70 [00:01<00:00, 42.71it/s]
84%|████████▍ | : 59/70 [00:01<00:00, 42.83it/s]
86%|████████▌ | : 60/70 [00:01<00:00, 42.83it/s]
87%|████████▋ | : 61/70 [00:01<00:00, 42.71it/s]
89%|████████▊ | : 62/70 [00:01<00:00, 42.53it/s]
90%|█████████ | : 63/70 [00:01<00:00, 41.99it/s]
91%|█████████▏| : 64/70 [00:01<00:00, 42.02it/s]
93%|█████████▎| : 65/70 [00:01<00:00, 42.11it/s]
94%|█████████▍| : 66/70 [00:01<00:00, 42.19it/s]
96%|█████████▌| : 67/70 [00:01<00:00, 42.21it/s]
97%|█████████▋| : 68/70 [00:01<00:00, 41.68it/s]
99%|█████████▊| : 69/70 [00:01<00:00, 41.87it/s]
100%|██████████| : 70/70 [00:01<00:00, 41.94it/s]
100%|██████████| : 70/70 [00:01<00:00, 39.88it/s]
Writing STC to disk...
[done]
View PSD of sources in label
plt.plot(stc.times, stc.data.T)
plt.xlabel("Frequency (Hz)")
plt.ylabel("PSD (dB)")
plt.title("Source Power Spectrum (PSD)")
plt.show()
Total running time of the script: (0 minutes 6.676 seconds)
Estimated memory usage: 9 MB