1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
|
"""
.. _ex-eeg-csd:
=====================================================
Transform EEG data using current source density (CSD)
=====================================================
This script shows an example of how to use CSD
:footcite:`PerrinEtAl1987,PerrinEtAl1989,Cohen2014,KayserTenke2015`.
CSD takes the spatial Laplacian of the sensor signal (derivative in both
x and y). It does what a planar gradiometer does in MEG. Computing these
spatial derivatives reduces point spread. CSD transformed data have a sharper
or more distinct topography, reducing the negative impact of volume conduction.
"""
# Authors: Alex Rockhill <aprockhill@mailbox.org>
#
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
# %%
# sphinx_gallery_thumbnail_number = 6
import matplotlib.pyplot as plt
import numpy as np
import mne
from mne.datasets import sample
print(__doc__)
data_path = sample.data_path()
# %%
# Load sample subject data
meg_path = data_path / "MEG" / "sample"
raw = mne.io.read_raw_fif(meg_path / "sample_audvis_raw.fif")
raw = raw.pick(picks=["eeg", "eog", "ecg", "stim"], exclude="bads").load_data()
events = mne.find_events(raw)
raw.set_eeg_reference(projection=True).apply_proj()
# %%
# Plot the raw data and CSD-transformed raw data:
raw_csd = mne.preprocessing.compute_current_source_density(raw)
raw.plot()
raw_csd.plot()
# %%
# Also look at the power spectral densities:
raw.compute_psd().plot(picks="data", exclude="bads", amplitude=False)
raw_csd.compute_psd().plot(picks="data", exclude="bads", amplitude=False)
# %%
# CSD can also be computed on Evoked (averaged) data.
# Here we epoch and average the data so we can demonstrate that.
event_id = {
"auditory/left": 1,
"auditory/right": 2,
"visual/left": 3,
"visual/right": 4,
"smiley": 5,
"button": 32,
}
epochs = mne.Epochs(raw, events, event_id=event_id, tmin=-0.2, tmax=0.5, preload=True)
evoked = epochs["auditory"].average()
# %%
# First let's look at how CSD affects scalp topography:
times = np.array([-0.1, 0.0, 0.05, 0.1, 0.15])
evoked_csd = mne.preprocessing.compute_current_source_density(evoked)
evoked.plot_joint(title="Average Reference", show=False)
evoked_csd.plot_joint(title="Current Source Density")
# %%
# CSD has parameters ``stiffness`` and ``lambda2`` affecting smoothing and
# spline flexibility, respectively. Let's see how they affect the solution:
fig, ax = plt.subplots(4, 4, layout="constrained")
fig.set_size_inches(10, 10)
for i, lambda2 in enumerate([0, 1e-7, 1e-5, 1e-3]):
for j, m in enumerate([5, 4, 3, 2]):
this_evoked_csd = mne.preprocessing.compute_current_source_density(
evoked, stiffness=m, lambda2=lambda2
)
this_evoked_csd.plot_topomap(
0.1, axes=ax[i, j], contours=4, time_unit="s", colorbar=False, show=False
)
ax[i, j].set_title(f"stiffness={m}\nλ²={lambda2}")
# %%
# References
# ----------
# .. footbibliography::
|
