.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_tutorials/04_sensor_simulation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_tutorials_04_sensor_simulation.py: 04. Sensor Simulation ===================== This tutorial mainly explains the ``SensorSimulator`` class. It projects source activity through a leadfield, then adds Gaussian sensor noise with the ``alpha_SNR`` mixing rule implemented by ``SensorSimulator``. The examples cover: - fixed-orientation projection; - free-orientation EEG projection; - free-orientation MEG projection; - the effect of different ``alpha_SNR`` settings; - workflow noise-variance quantities used downstream by source estimation. .. GENERATED FROM PYTHON SOURCE LINES 21-45 Scientific motivation --------------------- Source simulation produces ground-truth dipole moments at source locations, but inverse solvers operate on sensor measurements. ``SensorSimulator`` bridges this gap by applying the forward model: ``y_clean = L x`` and then adding white sensor noise: ``y_noisy = y_clean + eta * eps``. Here ``eps`` is Gaussian noise and ``eta`` is chosen from ``alpha_SNR`` using Frobenius norms. The result is a controlled sensor-level dataset that can be passed to inverse solvers. In the current CaliBrain workflow, this sensor-level simulation also defines the three noise-variance modes used later by source estimation: - ``oracle``: use the variance of the injected sensor noise; - ``baseline``: estimate variance from the pre-stimulus sensor segment; - ``adaptive_joint_learning``: do not fix ``noise_var`` here; let the solver learn it jointly from the data. .. GENERATED FROM PYTHON SOURCE LINES 45-55 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from mne.io.constants import FIFF from calibrain import SensorSimulator, SourceSimulator RANDOM_SEED = 31 .. GENERATED FROM PYTHON SOURCE LINES 56-64 Configure source simulation --------------------------- Source amplitudes are in ``nAm``. The leadfield units determine the sensor units: - if ``L`` is in ``fT / nAm``, then sensor outputs are in ``fT``; - if ``L`` is in ``µV / nAm``, then sensor outputs are in ``µV``. .. GENERATED FROM PYTHON SOURCE LINES 64-89 .. code-block:: Python erp_config = { "tmin": -0.1, "tmax": 0.8, "stim_onset": 0.0, "sfreq": 100, "fmin": 2, "fmax": 8, "amplitude_distribution": { "median": 8.0, "sigma": 0.15, "clip": [2.0, 20.0], }, "random_erp_timing": False, "erp_min_length": 20, } source_simulator = SourceSimulator(ERP_config=erp_config) sensor_simulator = SensorSimulator() times = np.arange( erp_config["tmin"], erp_config["tmax"], 1.0 / erp_config["sfreq"], ) .. GENERATED FROM PYTHON SOURCE LINES 90-98 Build small tutorial leadfields ------------------------------- These are synthetic leadfields for demonstration only: - fixed: shape ``(n_sensors, n_sources)``; - free EEG: shape ``(n_sensors, n_sources, 3)`` in ``µV / nAm``; - free MEG: shape ``(n_sensors, n_sources, 2)`` in ``fT / nAm``. .. GENERATED FROM PYTHON SOURCE LINES 98-107 .. code-block:: Python rng = np.random.default_rng(RANDOM_SEED) n_sensors = 16 n_sources = 40 leadfield_fixed = rng.normal(scale=0.15, size=(n_sensors, n_sources)) leadfield_free_eeg = rng.normal(scale=0.03, size=(n_sensors, n_sources, 3)) leadfield_free_meg = rng.normal(scale=4.0, size=(n_sensors, n_sources, 2)) .. GENERATED FROM PYTHON SOURCE LINES 108-112 Simulate source activity ------------------------ The three source arrays match the three leadfield conventions used below. .. GENERATED FROM PYTHON SOURCE LINES 112-140 .. code-block:: Python x_fixed, active_fixed = source_simulator.simulate( n_sources=n_sources, nnz=4, orientation_type="fixed", seed=RANDOM_SEED, ) x_free_eeg, active_free_eeg = source_simulator.simulate( n_sources=n_sources, nnz=4, orientation_type="free", coil_type=FIFF.FIFFV_COIL_EEG, seed=RANDOM_SEED, ) x_free_meg, active_free_meg = source_simulator.simulate( n_sources=n_sources, nnz=4, orientation_type="free", coil_type=FIFF.FIFFV_COIL_VV_MAG_T1, seed=RANDOM_SEED, ) print("fixed source shape:", x_fixed.shape) print("free EEG source shape:", x_free_eeg.shape) print("free MEG source shape:", x_free_meg.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none fixed source shape: (40, 90) free EEG source shape: (40, 3, 90) free MEG source shape: (40, 2, 90) .. GENERATED FROM PYTHON SOURCE LINES 141-146 Fixed-orientation sensor simulation ----------------------------------- For fixed orientation, projection is simple matrix multiplication ``y_clean = L @ x``. .. GENERATED FROM PYTHON SOURCE LINES 146-166 .. code-block:: Python sensor_simulator.set_sensor_metadata( kind=FIFF.FIFFV_MEG_CH, units=FIFF.FIFF_UNIT_T, unitmult=FIFF.FIFF_UNITM_F, coil_type=FIFF.FIFFV_COIL_VV_MAG_T1, ) y_fixed_clean, y_fixed_noisy, fixed_noise, fixed_eta = sensor_simulator.simulate( x=x_fixed, L=leadfield_fixed, alpha_SNR=0.7, sensor_white_noise_std=1.0, seed=RANDOM_SEED, ) print("fixed sensor shape:", y_fixed_clean.shape) print("fixed eta:", fixed_eta) print("sensor units:", sensor_simulator.units, "unitmult:", sensor_simulator.unitmult) .. rst-class:: sphx-glr-script-out .. code-block:: none fixed sensor shape: (16, 90) fixed eta: 0.36793688114720363 sensor units: 112 (FIFF_UNIT_T) unitmult: -15 (FIFF_UNITM_F) .. GENERATED FROM PYTHON SOURCE LINES 167-172 Free-orientation EEG sensor simulation -------------------------------------- EEG free orientation uses three local components per source, so the leadfield shape is ``(n_sensors, n_sources, 3)``. .. GENERATED FROM PYTHON SOURCE LINES 172-191 .. code-block:: Python sensor_simulator.set_sensor_metadata( kind=FIFF.FIFFV_EEG_CH, units=FIFF.FIFF_UNIT_V, unitmult=FIFF.FIFF_UNITM_MU, coil_type=FIFF.FIFFV_COIL_EEG, ) y_eeg_clean, y_eeg_noisy, eeg_noise, eeg_eta = sensor_simulator.simulate( x=x_free_eeg, L=leadfield_free_eeg, alpha_SNR=0.7, sensor_white_noise_std=0.1, seed=RANDOM_SEED, ) print("free EEG sensor shape:", y_eeg_clean.shape) print("free EEG eta:", eeg_eta) .. rst-class:: sphx-glr-script-out .. code-block:: none free EEG sensor shape: (16, 90) free EEG eta: 1.1724223634943463 .. GENERATED FROM PYTHON SOURCE LINES 192-197 Free-orientation MEG sensor simulation -------------------------------------- MEG free orientation uses the reduced two-component representation, so the leadfield shape is ``(n_sensors, n_sources, 2)``. .. GENERATED FROM PYTHON SOURCE LINES 197-216 .. code-block:: Python sensor_simulator.set_sensor_metadata( kind=FIFF.FIFFV_MEG_CH, units=FIFF.FIFF_UNIT_T, unitmult=FIFF.FIFF_UNITM_F, coil_type=FIFF.FIFFV_COIL_VV_MAG_T1, ) y_meg_clean, y_meg_noisy, meg_noise, meg_eta = sensor_simulator.simulate( x=x_free_meg, L=leadfield_free_meg, alpha_SNR=0.7, sensor_white_noise_std=1.0, seed=RANDOM_SEED, ) print("free MEG sensor shape:", y_meg_clean.shape) print("free MEG eta:", meg_eta) .. rst-class:: sphx-glr-script-out .. code-block:: none free MEG sensor shape: (16, 90) free MEG eta: 13.298082191264779 .. GENERATED FROM PYTHON SOURCE LINES 217-228 Workflow noise-variance estimates --------------------------------- The workflow does not use ``alpha_SNR`` directly as a solver input. Instead, after sensor simulation it derives named noise-variance modes: - ``oracle`` uses ``var(noise)``; - ``baseline`` uses the pre-stimulus segment of ``y_noisy``; - ``adaptive_joint_learning`` passes no fixed ``noise_var`` to the solver. The baseline estimate depends on ``tmin``, ``stim_onset``, and ``sfreq``. .. GENERATED FROM PYTHON SOURCE LINES 228-242 .. code-block:: Python tmin = erp_config["tmin"] stim_onset = erp_config["stim_onset"] sfreq = erp_config["sfreq"] pre_stimulus_onset = int((stim_onset - tmin) * sfreq) y_pre = y_fixed_noisy[:, :pre_stimulus_onset] oracle_noise_var = float(np.var(fixed_noise)) baseline_noise_var = float(np.mean(np.std(y_pre, axis=1) ** 2)) print("oracle noise variance:", oracle_noise_var) print("baseline noise variance:", baseline_noise_var) print("adaptive joint learning noise variance:", None) .. rst-class:: sphx-glr-script-out .. code-block:: none oracle noise variance: 0.13790136520275284 baseline noise variance: 0.12615502818075974 adaptive joint learning noise variance: None .. GENERATED FROM PYTHON SOURCE LINES 243-251 Effect of alpha_SNR ------------------- ``alpha_SNR`` is a mixing parameter in ``[0, 1]``: - ``1.0`` means no added noise; - ``0.0`` means pure scaled noise; - intermediate values mix signal and noise. .. GENERATED FROM PYTHON SOURCE LINES 251-265 .. code-block:: Python alpha_values = [1.0, 0.7, 0.3, 0.0] alpha_outputs = [] for alpha in alpha_values: alpha_outputs.append( sensor_simulator.simulate( x=x_fixed, L=leadfield_fixed, alpha_SNR=alpha, sensor_white_noise_std=1.0, seed=RANDOM_SEED, ) ) .. GENERATED FROM PYTHON SOURCE LINES 266-272 Plot sensor traces ------------------ The first plot compares clean and noisy fixed-orientation sensor traces for one representative sensor. The second plot shows how ``alpha_SNR`` changes the observed trace. .. GENERATED FROM PYTHON SOURCE LINES 272-294 .. code-block:: Python sensor_idx = 0 fig, axes = plt.subplots(2, 1, figsize=(8, 6), sharex=True) axes[0].plot(times, y_fixed_clean[sensor_idx], label="clean") axes[0].plot(times, y_fixed_noisy[sensor_idx], label="noisy") axes[0].set( ylabel="Sensor signal (tutorial units)", title="Fixed-orientation sensor trace", ) axes[0].legend(loc="best") for alpha, (_, y_noisy_alpha, _, _) in zip(alpha_values, alpha_outputs): axes[1].plot(times, y_noisy_alpha[sensor_idx], label=f"alpha_SNR = {alpha}") axes[1].set( xlabel="Time (s)", ylabel="Sensor signal (tutorial units)", title="Effect of alpha_SNR on one sensor", ) axes[1].legend(loc="best") fig.tight_layout() .. image-sg:: /auto_tutorials/images/sphx_glr_04_sensor_simulation_001.png :alt: Fixed-orientation sensor trace, Effect of alpha_SNR on one sensor :srcset: /auto_tutorials/images/sphx_glr_04_sensor_simulation_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 295-301 Compare sensor energy across modalities --------------------------------------- The absolute magnitudes differ because the synthetic leadfields use different unit conventions. The relevant comparison here is structural, not absolute: shape compatibility and the effect of noise mixing. .. GENERATED FROM PYTHON SOURCE LINES 301-310 .. code-block:: Python energy_summary = { "fixed clean": np.linalg.norm(y_fixed_clean, ord="fro"), "fixed noisy": np.linalg.norm(y_fixed_noisy, ord="fro"), "free EEG clean": np.linalg.norm(y_eeg_clean, ord="fro"), "free MEG clean": np.linalg.norm(y_meg_clean, ord="fro"), } print(energy_summary) .. rst-class:: sphx-glr-script-out .. code-block:: none {'fixed clean': np.float64(32.89446394819959), 'fixed noisy': np.float64(35.67080778140602), 'free EEG clean': np.float64(10.481744870962862), 'free MEG clean': np.float64(1188.8813207767196)} .. GENERATED FROM PYTHON SOURCE LINES 311-322 What this stage produces ------------------------ ``SensorSimulator.simulate`` returns: - ``y_clean``: noiseless sensor data; - ``y_noisy``: sensor data after additive white noise; - ``noise``: the added noise term; - ``eta``: the scaling factor implied by ``alpha_SNR``. These outputs are the direct inputs to CaliBrain inverse solvers. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.226 seconds) .. _sphx_glr_download_auto_tutorials_04_sensor_simulation.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 04_sensor_simulation.ipynb <04_sensor_simulation.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 04_sensor_simulation.py <04_sensor_simulation.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 04_sensor_simulation.zip <04_sensor_simulation.zip>`