EEG 信号频带功率计算

前言

       EEG 信号与大脑的活动和状态密切相关，在脑机接口等领域中有着广泛的应用，由此也衍生出了大量的 EEG 信号分析和特征提取方法，不同的子领域间侧重点可能略有差异，但一些基础的方法是相通的。例如，EEG 信号的功率计算就是一个十分基础且应用极其广泛的分析方法。本文以Raphael Vallet 的文章为基础，记录了一些我认为在 EEG 信号功率计算中需要注意的地方。（Raphael Vallet 写得十分清晰，可能的话，推荐直接看原文。）

信号功率谱密度（Power Spectral Density）计算

       功率谱密度（Power Spectral Density） 表征了信号功率在频域的分布状况。下图为一示例：

横轴一般单位为 Hz，纵轴在不同的领域中往往使用不同的单位，常用的有dB/Hz、W/Hz 等，在 EEG 信号分析中通常使用 μV²/Hz。本文不会涉及信号频域分析的数学推导和理论细节（感兴趣的朋友可以查阅教科书或其他资料），而是直接介绍3种基于 python 的计算 EEG 信号功率谱密度的方式。

基于 FFT 计算功率谱密度

       第一种方式根据功率谱密度的定义直接进行计算。首先对时序信号进行快速傅里叶变换（Fast Fourier Transformation，FFT），得到各频率对应的信号幅度和相位，通常对 FFT 结果的幅度取平方以表示功率谱密度，以下为一个示例

import numpy.fft as fft
import matplotlib.pyplot as plt
# EEG_data 维度为 1 * seq_len
complex_array = fft.fft(EEG_data)
power = np.abs(complex_array) ** 2
Fs = sample_rate  # 采样频率
T = 1 / Fs  # 采样周期
L = len(EEG_data)  # 信号长度
t = np.array([i * T for i in range(L)])
freqs = fft.fftfreq(t.size, t[1] - t[0])  # 得到分解波的频率序列
plt.plot(freqs[freqs > 0], pows[freqs > 0], color='orangered', label='Frequency')
plt.show()

得到的功率谱密度图如下（服务器上字体有些问题，请见谅）：

       直接 FFT 计算得到的功率谱密度 “has a good frequency resolution but way too much variance.”——Raphael Vallet

基于 scipy.signal.welch 计算功率谱密度

       第二种方法利用了 scipy.signal.welch，在直接 FFT 方法上做了一些改进，具体可以参考官方文档。以下为一个示例：

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# 示例 EEG 数据
data = np.loadtxt('data.txt')
# 定义采样率和时间向量
sf = 100
time = np.arange(data.size) / sf
# 定义窗口长度为4秒
win = 4 * sf
freqs, psd = signal.welch(data, sf, nperseg=win)
# 画功率谱密度图
sns.set(font_scale=1.2, style='white')
plt.figure(figsize=(8, 4))
plt.plot(freqs, psd, color='k', lw=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (V^2 / Hz)')
plt.ylim([0, psd.max() * 1.1])
plt.title("Welch's periodogram")
plt.xlim([0, freqs.max()])
sns.despine()

得到的功率谱密度图如下：

       利用 welch 计算得到的功率谱密度 “has a low variance, at the cost of a lower frequency resolution.”——Raphael Vallet

基于 mne.time_frequency.psd_array_multitaper 计算功率谱密度

       第三种方法是 multitaper 方法，最早由 David J. Thompson 于1982开发，用以克服经典频谱估算技术的一些限制，而 mne 则是一个专门用来处理 EEG, EMG, ECG 等生理信号的 Python 库。以下为一个示例：

import matplotlib.pyplot as plt
from mne.time_frequency import psd_array_multitaper
psd, freqs = psd_array_multitaper(data, sf, adaptive=True,normalization='full', verbose=0)
plt.plot(freqs, psd)
plt.show()

得到的功率谱密度图如下：

       利用 multitaper 计算得到的功率谱密度 “has the advantages of the two previous methods: high frequency resolution and low variance.”——Raphael Vallet
       三种计算方式的对比图如下所示：

特定频带绝对功率（Absolute Power）、相对功率(Relative Power)计算

       在 EEG 信号分析中，常常需要在功率谱密度的基础上计算特定频带的功率，而这也是我最初犯迷糊的地方——理论上，特定频带的功率只需要对这段频率范围的功率谱密度积分即可，然而上述三种方法得到的功率谱密度都是一系列离散的值，应该如何计算呢？
       最开始我用的方式是直接求和，如下所示：

import numpy as np
from mne.time_frequency import psd_array_multitaper
# 计算功率谱密度
psd, freqs = psd_array_multitaper(data, sf, adaptive=True,normalization='full', verbose=0)
# 指定频带范围
band = [lower_freq, upper_freq]
idx_band = np.logical_and(freqs >= band[0], freqs <= band[1])
band_power = sum(psd[idx_band])

       显然这种计算方式得到的结果和真实的频带功率相差甚远。那应该如何计算呢？可以使用 composite Simpson’s rule 进行估算，具体示例如下，使用了 scipy 提供的 simps 方法：

import numpy as np
from scipy.integrate import simps
from mne.time_frequency import psd_array_multitaper
# 计算功率谱密度
psd, freqs = psd_array_multitaper(data, sf, adaptive=True,normalization='full', verbose=0)
# 计算频率分辨率
freq_res = freqs[1] - freqs[0]
# 指定频带范围
band = [lower_freq, upper_freq]
idx_band = np.logical_and(freqs >= band[0], freqs <= band[1])
# psd 为一系列离散值，无法直接积分，因此采用 simps 做抛物线近似
band_power = simps(psd[idx_band], dx=freq_res)

       采用这种方法估算的频带功率更为准确。有了特定频带绝对功率的估算方式后，频带的相对功率计算便一目了然了band_power_relative = band_power_absolute / total_power_absolute。
       有了上述基础，类似于 spectral edge frequency 95（SEF95，the frequency such that 95% of the power in the spectrum lies below this value.） 的特征计算就是水到渠成的事情了。

References

1、https://raphaelvallat.com/bandpower.html
2、https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html
3、https://mathworld.wolfram.com/SimpsonsRule.html