Adding sensor space noise

import mne
import matplotlib.pyplot as plt

from mne.datasets import sample

from meegsim.location import select_random
from meegsim.simulate import SourceSimulator
from meegsim.waveform import narrowband_oscillation

First, we load the head model and associated source space:

# Paths
data_path = sample.data_path() / "MEG" / "sample"
fwd_path = data_path / "sample_audvis-meg-eeg-oct-6-fwd.fif"
raw_path = data_path / "sample_audvis_raw.fif"

# Load the prerequisites: fwd, src, and info
fwd = mne.read_forward_solution(fwd_path)
fwd = mne.convert_forward_solution(fwd, force_fixed=True)
raw = mne.io.read_raw(raw_path)
src = fwd["src"]
info = raw.info

# Pick EEG channels only
eeg_idx = mne.pick_types(info, eeg=True)
info_eeg = mne.pick_info(info, eeg_idx)
fwd_eeg = fwd.pick_channels(info_eeg.ch_names)

# Simulation parameters
sfreq = 250
duration = 60
seed = 123

sim = SourceSimulator(src)

# Select some vertices randomly
sim.add_point_sources(
    location=select_random,
    waveform=narrowband_oscillation,
    location_params=dict(n=3),
    waveform_params=dict(fmin=8, fmax=12),
    snr=2.5,
    snr_params=dict(fmin=8, fmax=12),
    names=["s1", "s2", "s3"],
)

sim.add_noise_sources(location=select_random, location_params=dict(n=10))

sc = sim.simulate(sfreq, duration, fwd=fwd, random_state=seed)
raw = sc.to_raw(fwd, info)

noise_levels = [0.05, 0.25, 0.5]
n_levels = len(noise_levels)
fig, axes = plt.subplots(ncols=n_levels, figsize=(3 * n_levels, 3))

for i_level, noise_level in enumerate(noise_levels):
    raw = sc.to_raw(fwd, info, sensor_noise_level=noise_level)

    spec = raw.compute_psd(fmax=60, n_fft=sfreq, n_overlap=sfreq // 2, n_per_seg=sfreq)
    spec.plot(axes=axes[i_level], amplitude=False)

    axes[i_level].set_title(f"{noise_level=}")
    axes[i_level].set_xlabel("Frequency (Hz)")

fig.tight_layout()

Reading forward solution from /home/docs/mne_data/MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif...
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    2 source spaces read
    Desired named matrix (kind = 3523 (FIFF_MNE_FORWARD_SOLUTION_GRAD)) not available
    Read MEG forward solution (7498 sources, 306 channels, free orientations)
    Desired named matrix (kind = 3523 (FIFF_MNE_FORWARD_SOLUTION_GRAD)) not available
    Read EEG forward solution (7498 sources, 60 channels, free orientations)
    Forward solutions combined: MEG, EEG
    Source spaces transformed to the forward solution coordinate frame
    Average patch normals will be employed in the rotation to the local surface coordinates....
    Converting to surface-based source orientations...
    [done]
Opening raw data file /home/docs/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.
/home/docs/checkouts/readthedocs.org/user_builds/meegsim/envs/latest/lib/python3.12/site-packages/meegsim/simulate.py:160: UserWarning: Ignoring the provided value of local SNR since global adjustment is enabled. To enable the local adjustment, set snr_mode to 'local' when initializing the SourceSimulator.
  warnings.warn(
Projecting source estimate to sensor space...
[done]
Creating RawArray with float64 data, n_channels=59, n_times=15000
    Range : 0 ... 14999 =      0.000 ...    59.996 secs
Ready.
Projecting source estimate to sensor space...
[done]
Creating RawArray with float64 data, n_channels=59, n_times=15000
    Range : 0 ... 14999 =      0.000 ...    59.996 secs
Ready.
Effective window size : 1.000 (s)
Plotting power spectral density (dB=True).
Projecting source estimate to sensor space...
[done]
Creating RawArray with float64 data, n_channels=59, n_times=15000
    Range : 0 ... 14999 =      0.000 ...    59.996 secs
Ready.
Effective window size : 1.000 (s)
Plotting power spectral density (dB=True).
Projecting source estimate to sensor space...
[done]
Creating RawArray with float64 data, n_channels=59, n_times=15000
    Range : 0 ... 14999 =      0.000 ...    59.996 secs
Ready.
Effective window size : 1.000 (s)
Plotting power spectral density (dB=True).
/home/docs/checkouts/readthedocs.org/user_builds/meegsim/checkouts/latest/examples/building_blocks/05_plot_sensor_space_noise.py:72: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  fig.tight_layout()