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()
fname_raw = data_path + '/MEG/sample/sample_audvis_raw.fif'
fname_event = data_path + '/MEG/sample/sample_audvis_raw-eve.fif'
raw = mne.io.read_raw_fif(fname_raw)
events = mne.read_events(fname_event)
```

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
Ready.
Current compensation grade : 0
```

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

Out:

```
Not setting metadata
Not setting metadata
72 matching events found
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)
```

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()
plt.suptitle('adaptive multitapers')
csd_wav.mean().plot()
plt.suptitle('Morlet wavelet transform')
```

**Total running time of the script:** ( 0 minutes 28.545 seconds)

**Estimated memory usage:** 223 MB