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:

  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>
# 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.

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:

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...
    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]
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...
    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:

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...
    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')
  • short-term Fourier transform, 15.8-20.0 Hz.
  • adaptive multitapers, 15.8-20.0 Hz.
  • Morlet wavelet transform, 16.0-20.0 Hz.

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

Estimated memory usage: 195 MB

Gallery generated by Sphinx-Gallery