- MNE version: MNE Python 3.11.7
- operating system: Windows 10
My final analysis had these ugly distorced waveforms:
Did anyone get something similar? What might be going wrong?
I’m pasting the whole code i’m using below:
# %% [markdown]
#
#
# # Análise fNIRS
#
# Script para análise dos experimentos realizados no escopo do projeto Sonhos e regulação emocional.
#
#
# %%
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
from IPython.display import clear_output
# %%
from itertools import compress
import matplotlib.pyplot as plt
import numpy as np
import mne
# pasta onde estão armazenados os dados coletados nos experimentos
fnirs_data_folder = "C:/Users/CogLab/Documents/0_Exp_Brasil_Pasta/DATA/fNIRS Data/NIRSTAR Data"
# caminho completo do conjunto de dados a ser analisado: pasta do arquivo no sistema + nome da pasta do participante
fnirs_cw_amplitude_dir = fnirs_data_folder + "/2024-06-17_013_eu_fNIRS_isolado_Sem_cruz"
raw_intensity = mne.io.read_raw_nirx(fnirs_cw_amplitude_dir, verbose=True)
raw_intensity.load_data()
# %% [markdown]
# ## Providing more meaningful annotation information
#
# First, we attribute more meaningful names to the trigger codes which are
# stored as annotations. Second, we include information about the duration of
# each stimulus, which was 5 seconds for all conditions in this experiment.
# Third, we remove the trigger code 15, which signaled the start and end
# of the experiment and is not relevant to our analysis.
#
#
# %%
print(raw_intensity.annotations.description)
# %% [markdown]
# O parâmetro para a função set_durations é duração da imagem na tela - ou seja, a duração do estÃmulo.
# %%
#escolhemos 9 segundos como duração do estÃmulo (1s de estÃmulo antes da instrução, 2s de instrução, 6s de estÃmulo após desaparcimento da instrução)
raw_intensity.annotations.set_durations(9)
raw_intensity.annotations.rename(
{"2.0": "Negativo/Visualizar", "4.0": "Negativo/Reavaliar", "8.0": "Neutro"}
)
unwanted = np.nonzero(raw_intensity.annotations.description == "3.0")
raw_intensity.annotations.delete(unwanted)
unwanted = np.nonzero(raw_intensity.annotations.description == "6.0")
raw_intensity.annotations.delete(unwanted)
unwanted = np.nonzero(raw_intensity.annotations.description == "BAD_SATURATED")
raw_intensity.annotations.delete(unwanted)
# %%
print(raw_intensity.annotations.description)
# %% [markdown]
# ## Viewing location of sensors over brain surface
#
# Here we validate that the location of sources-detector pairs and channels
# are in the expected locations. Source-detector pairs are shown as lines
# between the optodes, channels (the mid point of source-detector pairs) are
# optionally shown as orange dots. Source are optionally shown as red dots and
# detectors as black.
#
#
# %%
subjects_dir = mne.datasets.sample.data_path() / "subjects"
brain = mne.viz.Brain(
"fsaverage", subjects_dir=subjects_dir, background="w", cortex="0.5"
)
brain.add_sensors(
raw_intensity.info,
trans="fsaverage",
fnirs=["channels", "pairs", "sources", "detectors"],
)
brain.show_view(azimuth=20, elevation=60, distance=400)
# %% [markdown]
# ## Selecting channels appropriate for detecting neural responses
#
# First we remove channels that are too close together (short channels) to
# detect a neural response (less than 1 cm distance between optodes).
# These short channels can be seen in the figure above.
# To achieve this we pick all the channels that are not considered to be short.
#
#
# %%
picks = mne.pick_types(raw_intensity.info, meg=False, fnirs=True)
dists = mne.preprocessing.nirs.source_detector_distances(
raw_intensity.info, picks=picks
)
raw_intensity.pick(picks[dists > 0.01])
raw_intensity.plot(
n_channels=len(raw_intensity.ch_names), duration=500, show_scrollbars=False
)
# %% [markdown]
# ## Converting from raw intensity to optical density
#
# The raw intensity values are then converted to optical density.
#
#
# %%
raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)
raw_od.plot(n_channels=len(raw_od.ch_names), duration=500, show_scrollbars=False)
# %% [markdown]
# ## Evaluating the quality of the data
#
# At this stage we can quantify the quality of the coupling
# between the scalp and the optodes using the scalp coupling index. This
# method looks for the presence of a prominent synchronous signal in the
# frequency range of cardiac signals across both photodetected signals.
#
# In this example the data is clean and the coupling is good for all
# channels, so we will not mark any channels as bad based on the scalp
# coupling index.
#
#
# %%
sci = mne.preprocessing.nirs.scalp_coupling_index(raw_od)
fig, ax = plt.subplots(layout="constrained")
ax.hist(sci)
ax.set(xlabel="Scalp Coupling Index", ylabel="Count", xlim=[0, 1])
# %% [markdown]
# In this example we will mark all channels with a SCI less than 0.5 as bad
# (this dataset is quite clean, so no channels are marked as bad).
#
#
# %%
raw_od.info["bads"] = list(compress(raw_od.ch_names, sci < 0.5))
# %% [markdown]
# Now we're gonna plot the bad channels in two images: one is the signals, and the other is the channels represented on the head.
# %%
raw_od
raw_od.copy().pick(range(20)).plot(duration=300, n_channels=len(raw_od.ch_names), clipping=None)
plt.rcParams["figure.figsize"] = (6, 6) # (w, h)
raw_od.plot_sensors()
# %% [markdown]
# At this stage it is appropriate to inspect your data
# (for instructions on how to use the interactive data visualisation tool
# see `tut-visualize-raw`)
# to ensure that channels with poor scalp coupling have been removed.
# If your data contains lots of artifacts you may decide to apply
# artifact reduction techniques as described in `ex-fnirs-artifacts`.
#
#
# %% [markdown]
# ## Converting from optical density to haemoglobin
#
# Next we convert the optical density data to haemoglobin concentration using
# the modified Beer-Lambert law.
#
#
# %%
raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od, ppf=0.1)
raw_haemo.plot(n_channels=len(raw_haemo.ch_names), duration=500, show_scrollbars=False)
# %% [markdown]
# ## Removing heart rate from signal
#
# The haemodynamic response has frequency content predominantly below 0.5 Hz.
# An increase in activity around 1 Hz can be seen in the data that is due to
# the person's heart beat and is unwanted. So we use a low pass filter to
# remove this. A high pass filter is also included to remove slow drifts
# in the data.
#
#
# %%
raw_haemo_unfiltered = raw_haemo.copy()
raw_haemo.filter(0.05, 0.7, h_trans_bandwidth=0.2, l_trans_bandwidth=0.02)
for when, _raw in dict(Before=raw_haemo_unfiltered, After=raw_haemo).items():
fig = _raw.compute_psd().plot(
average=True, amplitude=False, picks="data", exclude="bads"
)
fig.suptitle(f"{when} filtering", weight="bold", size="x-large")
# %% [markdown]
# ## Extract epochs
#
# Now that the signal has been converted to relative haemoglobin concentration,
# and the unwanted heart rate component has been removed, we can extract epochs
# related to each of the experimental conditions.
#
# First we extract the events of interest and visualise them to ensure they are
# correct.
#
#
# %%
events, event_dict = mne.events_from_annotations(raw_haemo)
fig = mne.viz.plot_events(events, event_id=event_dict, sfreq=raw_haemo.info["sfreq"])
# %% [markdown]
# Next we define the range of our epochs, the rejection criteria,
# baseline correction, and extract the epochs. We visualise the log of which
# epochs were dropped.
#
#
# %%
reject_criteria = dict(hbo=80e-6) #O valor original era reject criteria = dict(hbo=80e-6)
#-2: 2 segundos antes de aparecer a imagem
#9: 9 segundos após o evento
# 25: 9 segundos de imagem, 4 de tempo estimado do SAM, 12 de ITI = 25
tmin, tmax = -2, 9
epochs = mne.Epochs(
raw_haemo,
events,
event_id=event_dict,
tmin=tmin,
tmax=tmax,
reject=reject_criteria,
reject_by_annotation=True,
proj=True,
#baseline=(10, 40),
baseline=(-2, -1),
preload=True,
detrend=None,
verbose=True,
)
epochs.plot_drop_log()
# %% [markdown]
# ## View consistency of responses across trials
#
# Now we can view the haemodynamic response for our condition.
# We visualise the response for both the oxy- and deoxyhaemoglobin, and
# observe the expected peak in HbO at around 6 seconds consistently across
# trials, and the consistent dip in HbR that is slightly delayed relative to
# the HbO peak.
#
#
# %%
epochs["Negativo/Visualizar"].plot_image(
combine="mean",
vmin=-30,
vmax=30,
ts_args=dict(ylim=dict(hbo=[-15, 15], hbr=[-15, 15])),
)
# %%
epochs["Negativo/Reavaliar"].plot_image(
combine="mean",
vmin=-30,
vmax=30,
ts_args=dict(ylim=dict(hbo=[-15, 15], hbr=[-15, 15])),
)
# %% [markdown]
# We can also view the epoched data for the control condition and observe
# that it does not show the expected morphology.
#
#
# %%
epochs["Neutro"].plot_image(
combine="mean",
vmin=-30,
vmax=30,
ts_args=dict(ylim=dict(hbo=[-15, 15], hbr=[-15, 15])),
)
# %% [markdown]
# ## View consistency of responses across channels
#
# Similarly we can view how consistent the response is across the optode
# pairs that we selected. All the channels in this data are located over the
# prefrontal cortex, and all channels show a similar pattern in the data.
#
#
# %%
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(20, 8), layout="constrained")
clims = dict(hbo=[-30, 30], hbr=[-30, 30])
# Plotagem das condições nos subplots
epochs["Negativo/Visualizar"].average().plot_image(axes=axes[0, :], clim=clims)
epochs["Negativo/Reavaliar"].average().plot_image(axes=axes[1, :], clim=clims)
epochs["Neutro"].average().plot_image(axes=axes[2, :], clim=clims)
# Definição dos tÃtulos para cada gráfico
conditions = ["Negativo/Visualizar", "Negativo/Reavaliar", "Neutro"]
for row, condition in enumerate(conditions):
for ax in axes[row, :]:
old_title = ax.get_title()
title = "{}: {}".format(condition, old_title)
ax.set_title(title)
# Exibição do gráfico
clear_output()
display(fig)
# %% [markdown]
# ## Plot standard fNIRS response image
#
# Next we generate the most common visualisation of fNIRS data: plotting
# both the HbO and HbR on the same figure to illustrate the relation between
# the two signals.
#
#
# %%
import matplotlib.pyplot as plt
import mne
evoked_dict = {
"Neg_Vis/HbO": epochs["Negativo/Visualizar"].average(picks="hbo"),
"Neg_Vis/HbR": epochs["Negativo/Visualizar"].average(picks="hbr"),
"Neg_Rea/HbO": epochs["Negativo/Reavaliar"].average(picks="hbo"),
"Neg_Rea/HbR": epochs["Negativo/Reavaliar"].average(picks="hbr"),
"Neu/HbO": epochs["Neutro"].average(picks="hbo"),
"Neu/HbR": epochs["Neutro"].average(picks="hbr"),
}
# Rename channels until the encoding of frequency in ch_name is fixed
for condition in evoked_dict:
evoked_dict[condition].rename_channels(lambda x: x[:-4])
# Define cores para cada condição
color_dict = {
"Neg_Vis/HbO": "#8B0000", # Vermelho escuro
"Neg_Vis/HbR": "#FF6347", # Vermelho claro
"Neg_Rea/HbO": "#006400", # Verde escuro
"Neg_Rea/HbR": "#98FB98", # Verde claro
"Neu/HbO": "#00008B", # Azul escuro
"Neu/HbR": "#87CEFA" # Azul claro
}
# Plot compare evokeds
fig, ax = plt.subplots(figsize=(10, 6))
mne.viz.plot_compare_evokeds(
evoked_dict, combine="mean", ci=0.95, colors=color_dict, axes=ax, show=False
)
# Ajuste a legenda para ter 2 linhas e 3 colunas na parte superior da figura
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, labels, ncol=3, loc='upper center', bbox_to_anchor=(0.5, 1.15), fontsize='small')
plt.show()
# %% [markdown]
# ## View topographic representation of activity (Oxyhemoglobin)
#
# Next we view how the topographic activity changes throughout the response (Oxyhemoglobin only)
#
# %%
# Define times within the valid range of your data (Os números abaixo são, respectivamente, o tempo mÃnimo e o tempo máximo a ser plotado, além do intervalo para plotar mais uma cabeça)
times = np.arange(-2, 9, 2.0)
# Now, use these corrected times in your plotting function
topomap_args = dict(extrapolate="local")
epochs["Negativo/Visualizar"].average(picks="hbo").plot_joint(
times=times, topomap_args=topomap_args, title="Negativo/Visualizar_Hbo"
)
epochs["Negativo/Reavaliar"].average(picks="hbo").plot_joint(
times=times, topomap_args=topomap_args, title="Negativo/Reavaliar_Hbo"
)
epochs["Neutro"].average(picks="hbo").plot_joint(
times=times, topomap_args=topomap_args, title="Neutro_Hbo"
)
# %%
print(raw_intensity.info['chs'])
# %% [markdown]
# ## Compare tapping of left and right hands
#
# Finally we generate topo maps for the left and right conditions to view
# the location of activity. First we visualise the HbO activity.
#
#
# %%
times = np.arange(-2.0, 9.0, 1.0)
epochs["Negativo/Visualizar"].average(picks="hbo").plot_topomap(times=times, **topomap_args)
epochs["Negativo/Reavaliar"].average(picks="hbo").plot_topomap(times=times, **topomap_args)
epochs["Neutro"].average(picks="hbo").plot_topomap(times=times, **topomap_args)
# %% [markdown]
# And we also view the HbR activity for the two conditions.
#
#
# %%
epochs["Negativo/Visualizar"].average(picks="hbr").plot_topomap(times=times, **topomap_args)
epochs["Negativo/Reavaliar"].average(picks="hbr").plot_topomap(times=times, **topomap_args)
epochs["Neutro"].average(picks="hbr").plot_topomap(times=times, **topomap_args)
# %% [markdown]
# And we can plot the comparison at a single time point for two conditions.
#
#
# %%
fig, axes = plt.subplots(
nrows=2,
ncols=4,
figsize=(9, 5),
gridspec_kw=dict(width_ratios=[1, 1, 1, 0.1]),
layout="constrained",
)
vlim = (-8, 8)
ts = 6.0
evoked_left = epochs["Negativo/Visualizar"].average()
evoked_right = epochs["Negativo/Reavaliar"].average()
evoked_right = epochs["Neutro"].average()
evoked_left.plot_topomap(
ch_type="hbo", times=ts, axes=axes[0, 0], vlim=vlim, colorbar=False, **topomap_args
)
evoked_left.plot_topomap(
ch_type="hbr", times=ts, axes=axes[1, 0], vlim=vlim, colorbar=False, **topomap_args
)
evoked_right.plot_topomap(
ch_type="hbo", times=ts, axes=axes[0, 1], vlim=vlim, colorbar=False, **topomap_args
)
evoked_right.plot_topomap(
ch_type="hbr", times=ts, axes=axes[1, 1], vlim=vlim, colorbar=False, **topomap_args
)
evoked_diff = mne.combine_evoked([evoked_left, evoked_right], weights=[1, -1])
evoked_diff.plot_topomap(
ch_type="hbo", times=ts, axes=axes[0, 2:], vlim=vlim, colorbar=True, **topomap_args
)
evoked_diff.plot_topomap(
ch_type="hbr", times=ts, axes=axes[1, 2:], vlim=vlim, colorbar=True, **topomap_args
)
for column, condition in enumerate(["Negativo_Visualizar", "Negativo_Reavaliar", "Visuazlizar-Reavaliar"]):
for row, chroma in enumerate(["HbO", "HbR"]):
axes[row, column].set_title("{}: {}".format(chroma, condition))
# %% [markdown]
# Lastly, we can also look at the individual waveforms to see what is
# driving the topographic plot above.
#
#
# %%
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 4), layout="constrained")
mne.viz.plot_evoked_topo(
epochs["Visualizar"].average(picks="hbo"), color="b", axes=axes, legend=False
)
mne.viz.plot_evoked_topo(
epochs["Reavaliar"].average(picks="hbo"), color="r", axes=axes, legend=False
)
# Tidy the legend:
leg_lines = [line for line in axes.lines if line.get_c() == "b"][:1]
leg_lines.append([line for line in axes.lines if line.get_c() == "r"][0])
fig.legend(leg_lines, ["Visualizar", "Reavaliar"], loc="lower right")