Simulate an oscillator on the cortex.
(fwd, idx)
| 55 | |
| 56 | |
| 57 | def _simulate_data(fwd, idx): # Somewhere on the frontal lobe by default |
| 58 | """Simulate an oscillator on the cortex.""" |
| 59 | pytest.importorskip("nibabel") |
| 60 | source_vertno = fwd["src"][0]["vertno"][idx] |
| 61 | |
| 62 | sfreq = 50.0 # Hz. |
| 63 | times = np.arange(10 * sfreq) / sfreq # 10 seconds of data |
| 64 | signal = np.sin(20 * 2 * np.pi * times) # 20 Hz oscillator |
| 65 | signal[: len(times) // 2] *= 2 # Make signal louder at the beginning |
| 66 | signal *= 1e-9 # Scale to be in the ballpark of MEG data |
| 67 | |
| 68 | # Construct a SourceEstimate object that describes the signal at the |
| 69 | # cortical level. |
| 70 | stc = mne.SourceEstimate( |
| 71 | signal[np.newaxis, :], |
| 72 | vertices=[[source_vertno], []], |
| 73 | tmin=0, |
| 74 | tstep=1 / sfreq, |
| 75 | subject="sample", |
| 76 | ) |
| 77 | |
| 78 | # Create an info object that holds information about the sensors |
| 79 | info = mne.create_info(fwd["info"]["ch_names"], sfreq, ch_types="grad") |
| 80 | with info._unlock(): |
| 81 | info.update(fwd["info"]) # Merge in sensor position information |
| 82 | # heavily decimate sensors to make it much faster |
| 83 | info = mne.pick_info(info, np.arange(info["nchan"])[::5]) |
| 84 | fwd = mne.pick_channels_forward(fwd, info["ch_names"]) |
| 85 | |
| 86 | # Run the simulated signal through the forward model, obtaining |
| 87 | # simulated sensor data. |
| 88 | raw = mne.apply_forward_raw(fwd, stc, info) |
| 89 | |
| 90 | # Add a little noise |
| 91 | random = np.random.RandomState(42) |
| 92 | noise = random.randn(*raw._data.shape) * 1e-14 |
| 93 | raw._data += noise |
| 94 | |
| 95 | # Define a single epoch (weird baseline but shouldn't matter) |
| 96 | epochs = mne.Epochs( |
| 97 | raw, |
| 98 | [[0, 0, 1]], |
| 99 | event_id=1, |
| 100 | tmin=0, |
| 101 | tmax=raw.times[-1], |
| 102 | baseline=(0.0, 0.0), |
| 103 | preload=True, |
| 104 | ) |
| 105 | evoked = epochs.average() |
| 106 | |
| 107 | # Compute the cross-spectral density matrix |
| 108 | csd = csd_morlet(epochs, frequencies=[10, 20], n_cycles=[5, 10], decim=5) |
| 109 | |
| 110 | labels = mne.read_labels_from_annot("sample", hemi="lh", subjects_dir=subjects_dir) |
| 111 | label = [label for label in labels if np.isin(source_vertno, label.vertices)] |
| 112 | assert len(label) == 1 |
| 113 | label = label[0] |
| 114 | vertices = np.intersect1d(label.vertices, fwd["src"][0]["vertno"]) |
no test coverage detected