使用mne Python库获取EEG信号的平均频带功率

Question

我正在寻找一种方法来从 edf 文件中获取 EEG 信号中通道的伽马波段平均频率，但我无法弄清楚如何做到这一点。我在线检查了各种来源，发现为了做到这一点，我需要使用韦尔奇方法从信号中获取PSD（功率谱密度），但我无法找到一种方法来做到这一点使用多国图书馆。到目前为止我所能完成的所有工作都附在下面。我将非常感谢任何形式的帮助。

import mne
file = "H S1 EC.edf"
data = mne.io.read_raw_edf(file)
raw_data = data.get_data()
info = data.info
channels = data.ch_names

Answer 1

要获得信号的功率谱，您可以：

import mne

file = "H S1 EC.edf"
data = mne.io.read_raw_edf(file, preload=True)
psds, freqs = mne.time_frequency.psd_welch(data)

这为您提供了所有频率的频谱功率，因此要获得平均伽玛功率，您必须对与您的频率相对应的 psds 数组切片进行平均感兴趣的频段（与psds中的功率值相对应的频率可以在freqs数组中找到）。

然而，您似乎对平均伽马功率不感兴趣，而是对伽马频率感兴趣，这要棘手得多。 EEG 中的伽马功率通常较弱，并且该频率范围内的信号可能主要由肌肉伪影（以及可能的微扫视伪影）主导 - 因此，我怀疑您是否能够识别可靠的伽马峰值（不是伪影）。我建议您看一下 Joerg Hipp 和 Markus Siegel 于 2013 年发表的论文“Dissociatingneural gamma-bandactivity from cranial and eyemuscleactivity in EEG”，他们对预处理和分析数据以获得良好的估计给出了一些很好的建议神经元伽玛能量。

Answer 2

要在 MNE 中使用 Welch 方法获取 PSD，可以使用 psd_welch() 函数。下面是一个帮助您入门的示例代码片段：

import mne

file = "H S1 EC.edf"
data = mne.io.read_raw_edf(file)

# Define the frequency range of interest (gamma band)
freq_range = (30, 100)

# Obtain the PSD using Welch's method
psds, freqs = mne.time_frequency.psd_welch(data, fmin=freq_range[0], fmax=freq_range[1], n_fft=1024)

# Find the index of the frequency bin with the highest power
gamma_index = psds.mean(axis=0).argmax()

# Obtain the frequency corresponding to the highest power bin
gamma_freq = freqs[gamma_index]

# Print the average gamma frequency for each channel
for i, channel in enumerate(data.ch_names):
    print("Channel {}: {:.2f} Hz".format(channel, psds[i,gamma_index]))

在此代码中，我们首先将感兴趣的频率范围定义为伽马频带 (30-100 Hz)。然后，我们使用 psd_welch() 函数使用 Welch 方法（窗口大小为 1024 个样本）来获取该频率范围的 PSD。生成的 PSD 存储在 psds 变量中，相应的频率存储在 freqs 变量中。

然后我们找到每个通道具有最高功率的频率仓的索引，并获得相应的频率。最后，我们打印出每个通道的平均伽玛频率。

请注意，psd_welch() 函数将 PSD 作为 2D 数组返回，其中每行对应一个通道，每列对应一个频率仓。为了获得每个通道的平均伽玛频率，我们首先取频率轴（轴=1）的平均值以获得每个通道的平均功率，然后找到每个通道的最高功率箱的索引。

或者，要使用 Welch 方法获取 EEG 信号的 PSD 并计算伽马带的平均频率，您可以使用 YASA 库。下面是一些示例代码，展示了如何执行此操作：

import mne
import yasa

# Load EEG data
file = "H S1 EC.edf"
data = mne.io.read_raw_edf(file)
raw_data = data.get_data()
sf = data.info['sfreq']  # Sampling frequency

# Compute the power spectral density using Welch's method
psd, freqs = yasa.psd_welch(raw_data, sf=sf, fmax=100)

# Calculate the average frequency of the gamma band (30-100 Hz)
gamma_freqs = freqs[(freqs >= 30) & (freqs <= 100)]
gamma_psd = psd[:, (freqs >= 30) & (freqs <= 100)]
avg_gamma_freq = (gamma_freqs * gamma_psd).sum() / gamma_psd.sum()

print("Average gamma frequency:", avg_gamma_freq)

此代码使用 MNE 加载 EEG 数据，使用 YASA 的 psd_welch() 函数计算 PSD，然后通过取以下值的加权平均值来计算伽玛波段的平均频率：伽马范围内的频率段。请注意，psd_welch() 函数负责对 EEG 信号进行加窗和重叠，这是获得准确的 PSD 估计的重要步骤。

此外，您可以参考我的github 代码 如何计算脑电图信号的平均频带功率