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:

1. Using short-term Fourier transform: mne.time_frequency.csd_fourier()

2. Using a multitaper approach: mne.time_frequency.csd_multitaper()

3. Using Morlet wavelets: mne.time_frequency.csd_morlet()

# Author: Marijn van Vliet <w.m.vanvliet@gmail.com>
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


data_path = sample.data_path()
fname_raw = data_path + '/MEG/sample/sample_audvis_raw.fif'
fname_event = data_path + '/MEG/sample/sample_audvis_raw-eve.fif'


Out:

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


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),


Out:

72 matching events found
Applying baseline correction (mode: mean)
3 projection items activated


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)


Out:

Computing cross-spectral density from epochs...
Computing CSD matrix for epoch 1
Computing CSD matrix for epoch 2
Computing CSD matrix for epoch 3
Computing CSD matrix for epoch 4
Computing CSD matrix for epoch 5
Computing CSD matrix for epoch 6
Computing CSD matrix for epoch 7
Computing CSD matrix for epoch 8
Computing CSD matrix for epoch 9
Computing CSD matrix for epoch 10
Computing CSD matrix for epoch 11
Computing CSD matrix for epoch 12
Computing CSD matrix for epoch 13
Computing CSD matrix for epoch 14
Computing CSD matrix for epoch 15
Computing CSD matrix for epoch 16
Computing CSD matrix for epoch 17
Computing CSD matrix for epoch 18
Computing CSD matrix for epoch 19
Computing CSD matrix for epoch 20
Computing CSD matrix for epoch 21
Computing CSD matrix for epoch 22
Computing CSD matrix for epoch 23
Computing CSD matrix for epoch 24
Computing CSD matrix for epoch 25
Computing CSD matrix for epoch 26
Computing CSD matrix for epoch 27
Computing CSD matrix for epoch 28
Computing CSD matrix for epoch 29
Computing CSD matrix for epoch 30
Computing CSD matrix for epoch 31
Computing CSD matrix for epoch 32
Computing CSD matrix for epoch 33
Computing CSD matrix for epoch 34
Computing CSD matrix for epoch 35
Computing CSD matrix for epoch 36
Computing CSD matrix for epoch 37
Computing CSD matrix for epoch 38
Computing CSD matrix for epoch 39
Computing CSD matrix for epoch 40
Computing CSD matrix for epoch 41
Computing CSD matrix for epoch 42
Computing CSD matrix for epoch 43
Computing CSD matrix for epoch 44
Computing CSD matrix for epoch 45
Computing CSD matrix for epoch 46
Computing CSD matrix for epoch 47
Computing CSD matrix for epoch 48
Computing CSD matrix for epoch 49
Computing CSD matrix for epoch 50
Computing CSD matrix for epoch 51
Computing CSD matrix for epoch 52
Computing CSD matrix for epoch 53
Computing CSD matrix for epoch 54
Computing CSD matrix for epoch 55
Computing CSD matrix for epoch 56
Computing CSD matrix for epoch 57
Computing CSD matrix for epoch 58
Computing CSD matrix for epoch 59
Computing CSD matrix for epoch 60
Computing CSD matrix for epoch 61
Computing CSD matrix for epoch 62
Computing CSD matrix for epoch 63
Computing CSD matrix for epoch 64
Computing CSD matrix for epoch 65
Computing CSD matrix for epoch 66
Computing CSD matrix for epoch 67
Computing CSD matrix for epoch 68
Computing CSD matrix for epoch 69
Computing CSD matrix for epoch 70
Computing CSD matrix for epoch 71
Computing CSD matrix for epoch 72
[done]
Using multitaper spectrum estimation with 7 DPSS windows
Computing cross-spectral density from epochs...
Computing CSD matrix for epoch 1
Computing CSD matrix for epoch 2
Computing CSD matrix for epoch 3
Computing CSD matrix for epoch 4
Computing CSD matrix for epoch 5
Computing CSD matrix for epoch 6
Computing CSD matrix for epoch 7
Computing CSD matrix for epoch 8
Computing CSD matrix for epoch 9
Computing CSD matrix for epoch 10
Computing CSD matrix for epoch 11
Computing CSD matrix for epoch 12
Computing CSD matrix for epoch 13
Computing CSD matrix for epoch 14
Computing CSD matrix for epoch 15
Computing CSD matrix for epoch 16
Computing CSD matrix for epoch 17
Computing CSD matrix for epoch 18
Computing CSD matrix for epoch 19
Computing CSD matrix for epoch 20
Computing CSD matrix for epoch 21
Computing CSD matrix for epoch 22
Computing CSD matrix for epoch 23
Computing CSD matrix for epoch 24
Computing CSD matrix for epoch 25
Computing CSD matrix for epoch 26
Computing CSD matrix for epoch 27
Computing CSD matrix for epoch 28
Computing CSD matrix for epoch 29
Computing CSD matrix for epoch 30
Computing CSD matrix for epoch 31
Computing CSD matrix for epoch 32
Computing CSD matrix for epoch 33
Computing CSD matrix for epoch 34
Computing CSD matrix for epoch 35
Computing CSD matrix for epoch 36
Computing CSD matrix for epoch 37
Computing CSD matrix for epoch 38
Computing CSD matrix for epoch 39
Computing CSD matrix for epoch 40
Computing CSD matrix for epoch 41
Computing CSD matrix for epoch 42
Computing CSD matrix for epoch 43
Computing CSD matrix for epoch 44
Computing CSD matrix for epoch 45
Computing CSD matrix for epoch 46
Computing CSD matrix for epoch 47
Computing CSD matrix for epoch 48
Computing CSD matrix for epoch 49
Computing CSD matrix for epoch 50
Computing CSD matrix for epoch 51
Computing CSD matrix for epoch 52
Computing CSD matrix for epoch 53
Computing CSD matrix for epoch 54
Computing CSD matrix for epoch 55
Computing CSD matrix for epoch 56
Computing CSD matrix for epoch 57
Computing CSD matrix for epoch 58
Computing CSD matrix for epoch 59
Computing CSD matrix for epoch 60
Computing CSD matrix for epoch 61
Computing CSD matrix for epoch 62
Computing CSD matrix for epoch 63
Computing CSD matrix for epoch 64
Computing CSD matrix for epoch 65
Computing CSD matrix for epoch 66
Computing CSD matrix for epoch 67
Computing CSD matrix for epoch 68
Computing CSD matrix for epoch 69
Computing CSD matrix for epoch 70
Computing CSD matrix for epoch 71
Computing CSD matrix for epoch 72
[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)


Out:

Computing cross-spectral density from epochs...
Computing CSD matrix for epoch 1
Computing CSD matrix for epoch 2
Computing CSD matrix for epoch 3
Computing CSD matrix for epoch 4
Computing CSD matrix for epoch 5
Computing CSD matrix for epoch 6
Computing CSD matrix for epoch 7
Computing CSD matrix for epoch 8
Computing CSD matrix for epoch 9
Computing CSD matrix for epoch 10
Computing CSD matrix for epoch 11
Computing CSD matrix for epoch 12
Computing CSD matrix for epoch 13
Computing CSD matrix for epoch 14
Computing CSD matrix for epoch 15
Computing CSD matrix for epoch 16
Computing CSD matrix for epoch 17
Computing CSD matrix for epoch 18
Computing CSD matrix for epoch 19
Computing CSD matrix for epoch 20
Computing CSD matrix for epoch 21
Computing CSD matrix for epoch 22
Computing CSD matrix for epoch 23
Computing CSD matrix for epoch 24
Computing CSD matrix for epoch 25
Computing CSD matrix for epoch 26
Computing CSD matrix for epoch 27
Computing CSD matrix for epoch 28
Computing CSD matrix for epoch 29
Computing CSD matrix for epoch 30
Computing CSD matrix for epoch 31
Computing CSD matrix for epoch 32
Computing CSD matrix for epoch 33
Computing CSD matrix for epoch 34
Computing CSD matrix for epoch 35
Computing CSD matrix for epoch 36
Computing CSD matrix for epoch 37
Computing CSD matrix for epoch 38
Computing CSD matrix for epoch 39
Computing CSD matrix for epoch 40
Computing CSD matrix for epoch 41
Computing CSD matrix for epoch 42
Computing CSD matrix for epoch 43
Computing CSD matrix for epoch 44
Computing CSD matrix for epoch 45
Computing CSD matrix for epoch 46
Computing CSD matrix for epoch 47
Computing CSD matrix for epoch 48
Computing CSD matrix for epoch 49
Computing CSD matrix for epoch 50
Computing CSD matrix for epoch 51
Computing CSD matrix for epoch 52
Computing CSD matrix for epoch 53
Computing CSD matrix for epoch 54
Computing CSD matrix for epoch 55
Computing CSD matrix for epoch 56
Computing CSD matrix for epoch 57
Computing CSD matrix for epoch 58
Computing CSD matrix for epoch 59
Computing CSD matrix for epoch 60
Computing CSD matrix for epoch 61
Computing CSD matrix for epoch 62
Computing CSD matrix for epoch 63
Computing CSD matrix for epoch 64
Computing CSD matrix for epoch 65
Computing CSD matrix for epoch 66
Computing CSD matrix for epoch 67
Computing CSD matrix for epoch 68
Computing CSD matrix for epoch 69
Computing CSD matrix for epoch 70
Computing CSD matrix for epoch 71
Computing CSD matrix for epoch 72
[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()