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
import matplotlib.pyplot as plt
import mne
from mne import io
from mne.datasets import sample
from mne.minimum_norm import read_inverse_operator, compute_source_psd
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
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 | : 1/70 [00:00<00:07, 9.42it/s]
3%|2 | : 2/70 [00:00<00:04, 14.49it/s]
4%|4 | : 3/70 [00:00<00:03, 18.71it/s]
6%|5 | : 4/70 [00:00<00:03, 21.87it/s]
7%|7 | : 5/70 [00:00<00:02, 24.51it/s]
9%|8 | : 6/70 [00:00<00:02, 26.73it/s]
10%|# | : 7/70 [00:00<00:02, 28.43it/s]
11%|#1 | : 8/70 [00:00<00:02, 29.91it/s]
13%|#2 | : 9/70 [00:00<00:01, 31.13it/s]
14%|#4 | : 10/70 [00:00<00:01, 32.28it/s]
16%|#5 | : 11/70 [00:00<00:01, 33.42it/s]
17%|#7 | : 12/70 [00:00<00:01, 34.10it/s]
19%|#8 | : 13/70 [00:00<00:01, 34.71it/s]
20%|## | : 14/70 [00:00<00:01, 35.35it/s]
21%|##1 | : 15/70 [00:00<00:01, 35.93it/s]
23%|##2 | : 16/70 [00:00<00:01, 36.61it/s]
24%|##4 | : 17/70 [00:00<00:01, 37.05it/s]
26%|##5 | : 18/70 [00:00<00:01, 37.47it/s]
27%|##7 | : 19/70 [00:00<00:01, 37.75it/s]
29%|##8 | : 20/70 [00:00<00:01, 38.11it/s]
30%|### | : 21/70 [00:00<00:01, 38.46it/s]
31%|###1 | : 22/70 [00:00<00:01, 38.69it/s]
33%|###2 | : 23/70 [00:00<00:01, 39.06it/s]
34%|###4 | : 24/70 [00:00<00:01, 39.32it/s]
36%|###5 | : 25/70 [00:00<00:01, 39.48it/s]
37%|###7 | : 26/70 [00:00<00:01, 39.57it/s]
39%|###8 | : 27/70 [00:00<00:01, 39.74it/s]
40%|#### | : 28/70 [00:00<00:01, 39.67it/s]
41%|####1 | : 29/70 [00:00<00:01, 39.47it/s]
43%|####2 | : 30/70 [00:00<00:01, 39.34it/s]
44%|####4 | : 31/70 [00:00<00:00, 39.22it/s]
46%|####5 | : 32/70 [00:00<00:00, 39.08it/s]
47%|####7 | : 33/70 [00:00<00:00, 38.96it/s]
49%|####8 | : 34/70 [00:00<00:00, 38.88it/s]
50%|##### | : 35/70 [00:00<00:00, 38.80it/s]
51%|#####1 | : 36/70 [00:00<00:00, 38.68it/s]
53%|#####2 | : 37/70 [00:00<00:00, 38.60it/s]
54%|#####4 | : 38/70 [00:01<00:00, 38.59it/s]
56%|#####5 | : 39/70 [00:01<00:00, 38.54it/s]
57%|#####7 | : 40/70 [00:01<00:00, 38.53it/s]
59%|#####8 | : 41/70 [00:01<00:00, 38.53it/s]
60%|###### | : 42/70 [00:01<00:00, 38.52it/s]
61%|######1 | : 43/70 [00:01<00:00, 38.45it/s]
63%|######2 | : 44/70 [00:01<00:00, 38.41it/s]
64%|######4 | : 45/70 [00:01<00:00, 38.38it/s]
66%|######5 | : 46/70 [00:01<00:00, 38.30it/s]
67%|######7 | : 47/70 [00:01<00:00, 38.22it/s]
69%|######8 | : 48/70 [00:01<00:00, 38.13it/s]
70%|####### | : 49/70 [00:01<00:00, 38.05it/s]
71%|#######1 | : 50/70 [00:01<00:00, 38.02it/s]
73%|#######2 | : 51/70 [00:01<00:00, 38.01it/s]
74%|#######4 | : 52/70 [00:01<00:00, 37.99it/s]
76%|#######5 | : 53/70 [00:01<00:00, 37.94it/s]
77%|#######7 | : 54/70 [00:01<00:00, 38.07it/s]
79%|#######8 | : 55/70 [00:01<00:00, 38.37it/s]
80%|######## | : 56/70 [00:01<00:00, 38.38it/s]
81%|########1 | : 57/70 [00:01<00:00, 38.18it/s]
83%|########2 | : 58/70 [00:01<00:00, 38.38it/s]
84%|########4 | : 59/70 [00:01<00:00, 38.61it/s]
86%|########5 | : 60/70 [00:01<00:00, 38.75it/s]
87%|########7 | : 61/70 [00:01<00:00, 38.90it/s]
89%|########8 | : 62/70 [00:01<00:00, 39.03it/s]
90%|######### | : 63/70 [00:01<00:00, 39.23it/s]
91%|#########1| : 64/70 [00:01<00:00, 39.42it/s]
93%|#########2| : 65/70 [00:01<00:00, 39.55it/s]
94%|#########4| : 66/70 [00:01<00:00, 39.54it/s]
96%|#########5| : 67/70 [00:01<00:00, 39.68it/s]
97%|#########7| : 68/70 [00:01<00:00, 39.78it/s]
99%|#########8| : 69/70 [00:01<00:00, 39.89it/s]
100%|##########| : 70/70 [00:01<00:00, 39.99it/s]
100%|##########| : 70/70 [00:01<00:00, 38.19it/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 7.485 seconds)
Estimated memory usage: 10 MB