前言
信号功率谱密度（Power Spectral Density）计算
- 基于 FFT 计算功率谱密度
- 基于 scipy.signal.welch 计算功率谱密度
- 基于 mne.time_frequency.psd_array_multitaper 计算功率谱密度
特定频带绝对功率（Absolute Power）、相对功率(Relative Power)计算
References

前言

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

EEG 信号频带功率计算相关推荐

基于动态空间滤波的受损EEG信号的鲁棒性学习
文章来源于微信公众号(茗创科技),欢迎有兴趣的朋友搜索关注. 导读利用实验室环境外记录的脑电建立机器学习模型需要对噪声数据和随机丢失通道具有鲁棒性的方法.这一需求在处理稀疏脑电蒙太奇(1-6个通道) ...
多模态自编码器从EEG信号预测fNIRS静息态
导读本研究介绍了一种深度学习架构,用于评估40名癫痫患者的多模态脑电图(EEG)和功能性近红外光谱(fNIRS)记录.长短期记忆网络和卷积神经网络集成在一个多模态序列到序列的自编码器中.训练后的神经 ...
运动想象| EEG信号、共空间模式算法（CSP）
摘要作为一种特殊的人机交互模式,脑-机接口(brain-computer interface, BCI)技术成为了当前信息交互的研究热点.其中脑电信号(electroencephalography, ...
使用CNN-LSTM混合深度学习分类基于MUSE采集的运动想象EEG信号
点击上面"脑机接口社区"关注我们更多技术干货第一时间送达脑机接口(BrainComputer Interfaces)技术是将人脑与外部设备建立起直接的通路,在智能助残.人机工程 ...
信号归一化功率_信号的频谱频谱密度功率谱密度能量谱密度
这是我在知乎上的一个回答,鉴于很多朋友对这几个概念不是很清楚,就在公众中发一下. 这几个概念,对于刚学信号系统的同学甚至对于很多信号处理的老手来说,都是分不清楚的,下面我们就一一解释这几个概念. ...
基于眼球追踪和脑电波EEG信号的学习者注意力量化
原文地址一.引言这篇文章是非侵入式课堂追踪系统(NiCATS 2021IEEE)的工作后续,本研究的重点是了解如何使用眼球追踪数据(例如平均注视次数)和脑电图(EEG)信号数据来测量学生的注意力, ...
手把手教你使用 1D 卷积和 LSTM 混合模型做 EEG 信号识别
本文是由CSDN用户[Memory逆光]授权分享.主要介绍了使用 1D 卷积和 LSTM 混合模型做 EEG 信号识别.感谢Memory逆光! 内容包括:1. 数据集(1.1 数据集下载.1.2 数据 ...
信号归一化功率_UE低发射功率余量分析
1.功率余量基础原理 1.1 功率余量PH 云南 (图1) 1.2 功率余量报告PHR PHR,全称是Power Headroom Report,中文为功率余量报告,即UE向网侧报告功率余量的过程.这 ...
能量信号和功率信号的分别
首先要明确一点,这两种信号概念是建立在无穷大的时间积分的基础上的. 一.能量与功率判断一个信号是能量信号还是功率信号,首先需要计算其能量和功率. 不想看公式的可以直接跳到加粗的结论部分. 对于信号f ...

EEG 信号频带功率计算

目录

前言