Note
Click here to download the full example code
Compute a cross-spectral density (CSD) matrix#
A cross-spectral density (CSD) matrix is similar to a covariance matrix, but in the time-frequency domain. It is the first step towards computing sensor-to-sensor coherence or a DICS beamformer.
This script demonstrates the three methods that MNE-Python provides to compute the CSD:
Using short-term Fourier transform:
mne.time_frequency.csd_fourier()Using a multitaper approach:
mne.time_frequency.csd_multitaper()Using Morlet wavelets:
mne.time_frequency.csd_morlet()
# Author: Marijn van Vliet <w.m.vanvliet@gmail.com>
# License: BSD-3-Clause
from matplotlib import pyplot as plt
import mne
from mne.datasets import sample
from mne.time_frequency import csd_fourier, csd_multitaper, csd_morlet
print(__doc__)
In the following example, the computation of the CSD matrices can be
performed using multiple cores. Set n_jobs to a value >1 to select the
number of cores to use.
n_jobs = 1
Loading the sample dataset.
data_path = sample.data_path()
meg_path = data_path / 'MEG' / 'sample'
fname_raw = meg_path / 'sample_audvis_raw.fif'
fname_event = meg_path / 'sample_audvis_raw-eve.fif'
raw = mne.io.read_raw_fif(fname_raw)
events = mne.read_events(fname_event)
Opening raw data file /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...
Read a total of 3 projection items:
PCA-v1 (1 x 102) idle
PCA-v2 (1 x 102) idle
PCA-v3 (1 x 102) idle
Range : 25800 ... 192599 = 42.956 ... 320.670 secs
Ready.
By default, CSD matrices are computed using all MEG/EEG channels. When interpreting a CSD matrix with mixed sensor types, be aware that the measurement units, and thus the scalings, differ across sensors. In this example, for speed and clarity, we select a single channel type: gradiometers.
picks = mne.pick_types(raw.info, meg='grad')
# Make some epochs, based on events with trigger code 1
epochs = mne.Epochs(raw, events, event_id=1, tmin=-0.2, tmax=1,
picks=picks, baseline=(None, 0),
reject=dict(grad=4000e-13), preload=True)
Not setting metadata
72 matching events found
Setting baseline interval to [-0.19979521315838786, 0.0] sec
Applying baseline correction (mode: mean)
3 projection items activated
Loading data for 72 events and 722 original time points ...
0 bad epochs dropped
Computing CSD matrices using short-term Fourier transform and (adaptive) multitapers is straightforward:
csd_fft = csd_fourier(epochs, fmin=15, fmax=20, n_jobs=n_jobs)
csd_mt = csd_multitaper(epochs, fmin=15, fmax=20, adaptive=True, n_jobs=n_jobs)
Removing projector <Projection | PCA-v1, active : True, n_channels : 102>
Removing projector <Projection | PCA-v2, active : True, n_channels : 102>
Removing projector <Projection | PCA-v3, active : True, n_channels : 102>
Computing cross-spectral density from epochs...
0%| | CSD epoch blocks : 0/72 [00:00<?, ?it/s]
1%|1 | CSD epoch blocks : 1/72 [00:00<00:02, 26.58it/s]
4%|4 | CSD epoch blocks : 3/72 [00:00<00:01, 45.22it/s]
7%|6 | CSD epoch blocks : 5/72 [00:00<00:01, 52.26it/s]
10%|9 | CSD epoch blocks : 7/72 [00:00<00:01, 55.85it/s]
12%|#2 | CSD epoch blocks : 9/72 [00:00<00:01, 58.27it/s]
15%|#5 | CSD epoch blocks : 11/72 [00:00<00:01, 59.95it/s]
18%|#8 | CSD epoch blocks : 13/72 [00:00<00:00, 61.24it/s]
21%|## | CSD epoch blocks : 15/72 [00:00<00:00, 62.18it/s]
24%|##3 | CSD epoch blocks : 17/72 [00:00<00:00, 62.81it/s]
26%|##6 | CSD epoch blocks : 19/72 [00:00<00:00, 63.23it/s]
29%|##9 | CSD epoch blocks : 21/72 [00:00<00:00, 63.80it/s]
32%|###1 | CSD epoch blocks : 23/72 [00:00<00:00, 64.35it/s]
35%|###4 | CSD epoch blocks : 25/72 [00:00<00:00, 64.73it/s]
38%|###7 | CSD epoch blocks : 27/72 [00:00<00:00, 65.13it/s]
40%|#### | CSD epoch blocks : 29/72 [00:00<00:00, 65.34it/s]
43%|####3 | CSD epoch blocks : 31/72 [00:00<00:00, 65.39it/s]
46%|####5 | CSD epoch blocks : 33/72 [00:00<00:00, 65.54it/s]
49%|####8 | CSD epoch blocks : 35/72 [00:00<00:00, 65.73it/s]
51%|#####1 | CSD epoch blocks : 37/72 [00:00<00:00, 65.84it/s]
54%|#####4 | CSD epoch blocks : 39/72 [00:00<00:00, 65.90it/s]
57%|#####6 | CSD epoch blocks : 41/72 [00:00<00:00, 65.93it/s]
60%|#####9 | CSD epoch blocks : 43/72 [00:00<00:00, 65.98it/s]
62%|######2 | CSD epoch blocks : 45/72 [00:00<00:00, 66.09it/s]
65%|######5 | CSD epoch blocks : 47/72 [00:00<00:00, 66.26it/s]
68%|######8 | CSD epoch blocks : 49/72 [00:00<00:00, 66.24it/s]
71%|####### | CSD epoch blocks : 51/72 [00:00<00:00, 66.33it/s]
74%|#######3 | CSD epoch blocks : 53/72 [00:00<00:00, 66.33it/s]
76%|#######6 | CSD epoch blocks : 55/72 [00:00<00:00, 66.36it/s]
79%|#######9 | CSD epoch blocks : 57/72 [00:00<00:00, 66.41it/s]
82%|########1 | CSD epoch blocks : 59/72 [00:00<00:00, 66.54it/s]
85%|########4 | CSD epoch blocks : 61/72 [00:00<00:00, 66.61it/s]
88%|########7 | CSD epoch blocks : 63/72 [00:00<00:00, 66.69it/s]
90%|######### | CSD epoch blocks : 65/72 [00:00<00:00, 66.78it/s]
93%|#########3| CSD epoch blocks : 67/72 [00:01<00:00, 66.83it/s]
96%|#########5| CSD epoch blocks : 69/72 [00:01<00:00, 66.94it/s]
99%|#########8| CSD epoch blocks : 71/72 [00:01<00:00, 66.94it/s]
100%|##########| CSD epoch blocks : 72/72 [00:01<00:00, 65.96it/s]
[done]
Removing projector <Projection | PCA-v1, active : True, n_channels : 102>
Removing projector <Projection | PCA-v2, active : True, n_channels : 102>
Removing projector <Projection | PCA-v3, active : True, n_channels : 102>
Using multitaper spectrum estimation with 7 DPSS windows
Computing cross-spectral density from epochs...
0%| | CSD epoch blocks : 0/72 [00:00<?, ?it/s]
1%|1 | CSD epoch blocks : 1/72 [00:00<00:15, 4.61it/s]
3%|2 | CSD epoch blocks : 2/72 [00:00<00:11, 5.86it/s]
4%|4 | CSD epoch blocks : 3/72 [00:00<00:10, 6.47it/s]
6%|5 | CSD epoch blocks : 4/72 [00:00<00:09, 6.82it/s]
7%|6 | CSD epoch blocks : 5/72 [00:00<00:09, 7.06it/s]
8%|8 | CSD epoch blocks : 6/72 [00:00<00:09, 7.22it/s]
10%|9 | CSD epoch blocks : 7/72 [00:00<00:08, 7.35it/s]
11%|#1 | CSD epoch blocks : 8/72 [00:01<00:08, 7.45it/s]
12%|#2 | CSD epoch blocks : 9/72 [00:01<00:08, 7.52it/s]
14%|#3 | CSD epoch blocks : 10/72 [00:01<00:08, 7.58it/s]
15%|#5 | CSD epoch blocks : 11/72 [00:01<00:07, 7.63it/s]
17%|#6 | CSD epoch blocks : 12/72 [00:01<00:07, 7.67it/s]
18%|#8 | CSD epoch blocks : 13/72 [00:01<00:07, 7.70it/s]
19%|#9 | CSD epoch blocks : 14/72 [00:01<00:07, 7.73it/s]
21%|## | CSD epoch blocks : 15/72 [00:01<00:07, 7.76it/s]
22%|##2 | CSD epoch blocks : 16/72 [00:02<00:07, 7.79it/s]
24%|##3 | CSD epoch blocks : 17/72 [00:02<00:07, 7.81it/s]
25%|##5 | CSD epoch blocks : 18/72 [00:02<00:06, 7.83it/s]
26%|##6 | CSD epoch blocks : 19/72 [00:02<00:06, 7.84it/s]
28%|##7 | CSD epoch blocks : 20/72 [00:02<00:06, 7.86it/s]
29%|##9 | CSD epoch blocks : 21/72 [00:02<00:06, 7.87it/s]
31%|### | CSD epoch blocks : 22/72 [00:02<00:06, 7.88it/s]
32%|###1 | CSD epoch blocks : 23/72 [00:02<00:06, 7.89it/s]
33%|###3 | CSD epoch blocks : 24/72 [00:03<00:06, 7.90it/s]
35%|###4 | CSD epoch blocks : 25/72 [00:03<00:05, 7.92it/s]
36%|###6 | CSD epoch blocks : 26/72 [00:03<00:05, 7.92it/s]
38%|###7 | CSD epoch blocks : 27/72 [00:03<00:05, 7.93it/s]
39%|###8 | CSD epoch blocks : 28/72 [00:03<00:05, 7.93it/s]
40%|#### | CSD epoch blocks : 29/72 [00:03<00:05, 7.94it/s]
42%|####1 | CSD epoch blocks : 30/72 [00:03<00:05, 7.94it/s]
43%|####3 | CSD epoch blocks : 31/72 [00:03<00:05, 7.94it/s]
44%|####4 | CSD epoch blocks : 32/72 [00:04<00:05, 7.95it/s]
46%|####5 | CSD epoch blocks : 33/72 [00:04<00:04, 7.95it/s]
47%|####7 | CSD epoch blocks : 34/72 [00:04<00:04, 7.96it/s]
49%|####8 | CSD epoch blocks : 35/72 [00:04<00:04, 7.96it/s]
50%|##### | CSD epoch blocks : 36/72 [00:04<00:04, 7.96it/s]
51%|#####1 | CSD epoch blocks : 37/72 [00:04<00:04, 7.96it/s]
53%|#####2 | CSD epoch blocks : 38/72 [00:04<00:04, 7.96it/s]
54%|#####4 | CSD epoch blocks : 39/72 [00:04<00:04, 7.96it/s]
56%|#####5 | CSD epoch blocks : 40/72 [00:05<00:04, 7.96it/s]
57%|#####6 | CSD epoch blocks : 41/72 [00:05<00:03, 7.96it/s]
58%|#####8 | CSD epoch blocks : 42/72 [00:05<00:03, 7.96it/s]
60%|#####9 | CSD epoch blocks : 43/72 [00:05<00:03, 7.97it/s]
61%|######1 | CSD epoch blocks : 44/72 [00:05<00:03, 7.97it/s]
62%|######2 | CSD epoch blocks : 45/72 [00:05<00:03, 7.93it/s]
64%|######3 | CSD epoch blocks : 46/72 [00:05<00:03, 7.79it/s]
65%|######5 | CSD epoch blocks : 47/72 [00:06<00:03, 7.71it/s]
67%|######6 | CSD epoch blocks : 48/72 [00:06<00:03, 7.72it/s]
68%|######8 | CSD epoch blocks : 49/72 [00:06<00:02, 7.74it/s]
69%|######9 | CSD epoch blocks : 50/72 [00:06<00:02, 7.75it/s]
71%|####### | CSD epoch blocks : 51/72 [00:06<00:02, 7.77it/s]
72%|#######2 | CSD epoch blocks : 52/72 [00:06<00:02, 7.78it/s]
74%|#######3 | CSD epoch blocks : 53/72 [00:06<00:02, 7.79it/s]
75%|#######5 | CSD epoch blocks : 54/72 [00:06<00:02, 7.80it/s]
76%|#######6 | CSD epoch blocks : 55/72 [00:07<00:02, 7.81it/s]
78%|#######7 | CSD epoch blocks : 56/72 [00:07<00:02, 7.82it/s]
79%|#######9 | CSD epoch blocks : 57/72 [00:07<00:01, 7.83it/s]
81%|######## | CSD epoch blocks : 58/72 [00:07<00:01, 7.84it/s]
82%|########1 | CSD epoch blocks : 59/72 [00:07<00:01, 7.85it/s]
83%|########3 | CSD epoch blocks : 60/72 [00:07<00:01, 7.86it/s]
85%|########4 | CSD epoch blocks : 61/72 [00:07<00:01, 7.86it/s]
86%|########6 | CSD epoch blocks : 62/72 [00:07<00:01, 7.87it/s]
88%|########7 | CSD epoch blocks : 63/72 [00:08<00:01, 7.87it/s]
89%|########8 | CSD epoch blocks : 64/72 [00:08<00:01, 7.88it/s]
90%|######### | CSD epoch blocks : 65/72 [00:08<00:00, 7.89it/s]
92%|#########1| CSD epoch blocks : 66/72 [00:08<00:00, 7.89it/s]
93%|#########3| CSD epoch blocks : 67/72 [00:08<00:00, 7.90it/s]
94%|#########4| CSD epoch blocks : 68/72 [00:08<00:00, 7.90it/s]
96%|#########5| CSD epoch blocks : 69/72 [00:08<00:00, 7.91it/s]
97%|#########7| CSD epoch blocks : 70/72 [00:08<00:00, 7.90it/s]
99%|#########8| CSD epoch blocks : 71/72 [00:09<00:00, 7.90it/s]
100%|##########| CSD epoch blocks : 72/72 [00:09<00:00, 7.91it/s]
100%|##########| CSD epoch blocks : 72/72 [00:09<00:00, 7.85it/s]
[done]
When computing the CSD with Morlet wavelets, you specify the exact frequencies at which to compute it. For each frequency, a corresponding wavelet will be constructed and convolved with the signal, resulting in a time-frequency decomposition.
The CSD is constructed by computing the correlation between the
time-frequency representations between all sensor-to-sensor pairs. The
time-frequency decomposition originally has the same sampling rate as the
signal, in our case ~600Hz. This means the decomposition is over-specified in
time and we may not need to use all samples during our CSD computation, just
enough to get a reliable correlation statistic. By specifying decim=10,
we use every 10th sample, which will greatly speed up the computation and
will have a minimal effect on the CSD.
frequencies = [16, 17, 18, 19, 20]
csd_wav = csd_morlet(epochs, frequencies, decim=10, n_jobs=n_jobs)
Removing projector <Projection | PCA-v1, active : True, n_channels : 102>
Removing projector <Projection | PCA-v2, active : True, n_channels : 102>
Removing projector <Projection | PCA-v3, active : True, n_channels : 102>
Computing cross-spectral density from epochs...
0%| | CSD epoch blocks : 0/72 [00:00<?, ?it/s]
1%|1 | CSD epoch blocks : 1/72 [00:00<00:07, 10.11it/s]
3%|2 | CSD epoch blocks : 2/72 [00:00<00:05, 12.12it/s]
4%|4 | CSD epoch blocks : 3/72 [00:00<00:05, 12.95it/s]
6%|5 | CSD epoch blocks : 4/72 [00:00<00:05, 13.41it/s]
7%|6 | CSD epoch blocks : 5/72 [00:00<00:04, 13.70it/s]
8%|8 | CSD epoch blocks : 6/72 [00:00<00:04, 13.88it/s]
10%|9 | CSD epoch blocks : 7/72 [00:00<00:04, 14.01it/s]
11%|#1 | CSD epoch blocks : 8/72 [00:00<00:04, 14.10it/s]
12%|#2 | CSD epoch blocks : 9/72 [00:00<00:04, 14.20it/s]
14%|#3 | CSD epoch blocks : 10/72 [00:00<00:04, 14.26it/s]
15%|#5 | CSD epoch blocks : 11/72 [00:00<00:04, 14.29it/s]
17%|#6 | CSD epoch blocks : 12/72 [00:00<00:04, 14.36it/s]
18%|#8 | CSD epoch blocks : 13/72 [00:00<00:04, 14.41it/s]
19%|#9 | CSD epoch blocks : 14/72 [00:00<00:04, 14.46it/s]
21%|## | CSD epoch blocks : 15/72 [00:01<00:03, 14.51it/s]
22%|##2 | CSD epoch blocks : 16/72 [00:01<00:03, 14.55it/s]
24%|##3 | CSD epoch blocks : 17/72 [00:01<00:03, 14.59it/s]
25%|##5 | CSD epoch blocks : 18/72 [00:01<00:03, 14.61it/s]
26%|##6 | CSD epoch blocks : 19/72 [00:01<00:03, 14.65it/s]
28%|##7 | CSD epoch blocks : 20/72 [00:01<00:03, 14.68it/s]
29%|##9 | CSD epoch blocks : 21/72 [00:01<00:03, 14.71it/s]
31%|### | CSD epoch blocks : 22/72 [00:01<00:03, 14.72it/s]
32%|###1 | CSD epoch blocks : 23/72 [00:01<00:03, 14.72it/s]
33%|###3 | CSD epoch blocks : 24/72 [00:01<00:03, 14.75it/s]
35%|###4 | CSD epoch blocks : 25/72 [00:01<00:03, 14.76it/s]
36%|###6 | CSD epoch blocks : 26/72 [00:01<00:03, 14.78it/s]
38%|###7 | CSD epoch blocks : 27/72 [00:01<00:03, 14.79it/s]
39%|###8 | CSD epoch blocks : 28/72 [00:01<00:02, 14.81it/s]
40%|#### | CSD epoch blocks : 29/72 [00:01<00:02, 14.84it/s]
42%|####1 | CSD epoch blocks : 30/72 [00:02<00:02, 14.85it/s]
43%|####3 | CSD epoch blocks : 31/72 [00:02<00:02, 14.87it/s]
44%|####4 | CSD epoch blocks : 32/72 [00:02<00:02, 14.91it/s]
46%|####5 | CSD epoch blocks : 33/72 [00:02<00:02, 14.93it/s]
47%|####7 | CSD epoch blocks : 34/72 [00:02<00:02, 14.95it/s]
49%|####8 | CSD epoch blocks : 35/72 [00:02<00:02, 14.97it/s]
50%|##### | CSD epoch blocks : 36/72 [00:02<00:02, 14.98it/s]
51%|#####1 | CSD epoch blocks : 37/72 [00:02<00:02, 15.00it/s]
53%|#####2 | CSD epoch blocks : 38/72 [00:02<00:02, 15.01it/s]
54%|#####4 | CSD epoch blocks : 39/72 [00:02<00:02, 15.03it/s]
56%|#####5 | CSD epoch blocks : 40/72 [00:02<00:02, 15.05it/s]
57%|#####6 | CSD epoch blocks : 41/72 [00:02<00:02, 15.06it/s]
58%|#####8 | CSD epoch blocks : 42/72 [00:02<00:01, 15.08it/s]
60%|#####9 | CSD epoch blocks : 43/72 [00:02<00:01, 15.09it/s]
61%|######1 | CSD epoch blocks : 44/72 [00:02<00:01, 15.10it/s]
62%|######2 | CSD epoch blocks : 45/72 [00:03<00:01, 15.10it/s]
64%|######3 | CSD epoch blocks : 46/72 [00:03<00:01, 15.10it/s]
65%|######5 | CSD epoch blocks : 47/72 [00:03<00:01, 15.10it/s]
67%|######6 | CSD epoch blocks : 48/72 [00:03<00:01, 15.10it/s]
68%|######8 | CSD epoch blocks : 49/72 [00:03<00:01, 15.10it/s]
69%|######9 | CSD epoch blocks : 50/72 [00:03<00:01, 15.11it/s]
71%|####### | CSD epoch blocks : 51/72 [00:03<00:01, 15.10it/s]
72%|#######2 | CSD epoch blocks : 52/72 [00:03<00:01, 15.11it/s]
74%|#######3 | CSD epoch blocks : 53/72 [00:03<00:01, 15.11it/s]
75%|#######5 | CSD epoch blocks : 54/72 [00:03<00:01, 15.11it/s]
76%|#######6 | CSD epoch blocks : 55/72 [00:03<00:01, 15.12it/s]
78%|#######7 | CSD epoch blocks : 56/72 [00:03<00:01, 15.09it/s]
79%|#######9 | CSD epoch blocks : 57/72 [00:03<00:00, 15.10it/s]
81%|######## | CSD epoch blocks : 58/72 [00:03<00:00, 15.10it/s]
82%|########1 | CSD epoch blocks : 59/72 [00:03<00:00, 15.11it/s]
83%|########3 | CSD epoch blocks : 60/72 [00:04<00:00, 15.11it/s]
85%|########4 | CSD epoch blocks : 61/72 [00:04<00:00, 15.09it/s]
86%|########6 | CSD epoch blocks : 62/72 [00:04<00:00, 15.09it/s]
88%|########7 | CSD epoch blocks : 63/72 [00:04<00:00, 15.07it/s]
89%|########8 | CSD epoch blocks : 64/72 [00:04<00:00, 15.04it/s]
90%|######### | CSD epoch blocks : 65/72 [00:04<00:00, 15.03it/s]
92%|#########1| CSD epoch blocks : 66/72 [00:04<00:00, 15.03it/s]
93%|#########3| CSD epoch blocks : 67/72 [00:04<00:00, 15.04it/s]
94%|#########4| CSD epoch blocks : 68/72 [00:04<00:00, 15.04it/s]
96%|#########5| CSD epoch blocks : 69/72 [00:04<00:00, 15.04it/s]
97%|#########7| CSD epoch blocks : 70/72 [00:04<00:00, 15.04it/s]
99%|#########8| CSD epoch blocks : 71/72 [00:04<00:00, 15.04it/s]
100%|##########| CSD epoch blocks : 72/72 [00:04<00:00, 15.04it/s]
100%|##########| CSD epoch blocks : 72/72 [00:04<00:00, 14.93it/s]
[done]
The resulting mne.time_frequency.CrossSpectralDensity objects have a
plotting function we can use to compare the results of the different methods.
We’re plotting the mean CSD across frequencies.
csd_fft.mean().plot()
plt.suptitle('short-term Fourier transform')
csd_mt.mean().plot()
plt.suptitle('adaptive multitapers')
csd_wav.mean().plot()
plt.suptitle('Morlet wavelet transform')
Total running time of the script: ( 0 minutes 20.751 seconds)
Estimated memory usage: 196 MB