【EEG分类】脑电信号提取PSD功率谱密度特征

【EEG 分类】脑电信号提取PSD 功率谱密度特征

1. 序⾔

脑电信号是⼀种⾮平稳的随机信号，⼀般⽽⾔随机信号的持续时间是⽆限长的，因此随机信号的总能量是⽆限的，⽽随机过程的任意⼀个样本函数都不满⾜绝对可积条件，所以其傅⾥叶变换不存在。

麻醉系统不过，尽管随机信号的总能量是⽆限的，但其平均功率却是有限的，因此，要对随机信号的频域进⾏分析，应从功率谱出发进⾏研究才有意义。正因如此，在研究中经常使⽤功率谱密度（Power spectral density, PSD）来分析脑电信号的频域特性。

本⽂简单的演⽰了对脑电信号提取功率谱密度特征然后分类的基本流程。希望对那些尚未⼊门、⾯对 BCI 任务不知所措的新⼿能有⼀点启发。

2. 功率谱密度理论基础简述

功率谱密度描是对随机变量均⽅值的量度，是单位频率的平均功率量纲。对功率谱在频域上积分就可以得到信号的平均功率。功率谱密度是⼀个以频率为⾃变量的映射，反映了在频率成分上信号有多少功率。

我们假定⼀个随机过程，并定义⼀个截断阈值，随机过程的截断过程就可以定义为

则该随机过程的能量可定义为

对能量函数求导，就可以获得平均功率。

根据 Parseval 定理（即能量从时域⾓度和频域⾓度来看都是相等的）可得：

这⾥是经过傅⾥叶变换后的形式。由于随机过程被限定在了⼀个有限的时间区间之间，所以对随机过程的傅⾥叶变换不再受限。另外我们还需要注意到，是⼀个随机变量，因此为了得到最终总体的平均功率，还需要求取随机变量的期望值。

由此，通过求取时的极限，就可以得到原始随机过程的平均功率。

S (f )f S (f )f X (t )t 0X X (t )t 0X (t )=t 0{X (t )0

∣t ∣≤t o

∣t ∣>t o E =X (t )dt =X (t )dt X t 0∫−t 02∫−∞t 02

P =X (t )dt

X t 02t o 1∫−∞t 02P =X (f )df

X t 02t o 1∫−∞∣∣t 02∣∣2X (f )t 0X (t )t 0X [−t ,t ]00P X t 0=E P =E X (f )df高压储罐

P X t 0[X t 0]2t o 1∫−∞[∣∣t 02∣∣2]t →0∞P X

将式中被积函数单独提取出来，定义为 :

这样⼀来，平均功率可以表⽰为。通过这种定义⽅式，函数可以表征每⼀个最⼩极限单位的频率分量所拥有的功率⼤⼩，因此我们把称为功率谱密度。

3. Matlab 中 PSD 函数的使⽤

功率谱密度的估计⽅法有很多。总体来讲可以分为两⼤类：传统的⾮参数⽅法，和现代的参数⽅法。

本节不对理论知识做详细的叙述，感兴趣的可以深⼊查阅⽂献，这⾥只介绍⼀下有哪些⽅法，以及他们在 matlab 当中的使⽤。

3.1 传统⾮参数⽅法估计 PSD 最简单的⽅法是周期图法，先对信号做 FFT 变换，然后求平⽅，periodogram 函数实现了这个功能。不过周期图法估计的⽅差随采样点数N 的增加⽽增加，不是很建议使⽤。

=E X (f )df =df P X t →∞0lim

2t 01∫−∞[∣t 0∣]∫−∞t →∞0lim 2t 0E X (f )[∣t 0∣]S (f )S (f )=t →∞0lim

2t 0E X (f )[∣t 0∣]

P X S (f )df ∫−∞∞

S (f )S (f )

另⼀种⾃相关⽅法，基于维纳⾟钦定律：信号的功率谱估计等于该信号⾃相关函数的离散DTFT，不过我没有在 matlab ⾥到对应的函数，如果有知道的朋友请告诉我⼀下。

最常⽤的函数是 pwelch 函数，利⽤ welch ⽅法来求 PSD，这也是最推荐使⽤的。

3.2 参数⽅法估计 PSD

包括 pconv、pburg、pyulear 等⼏个⽅法。

这些⽅法我没⽤过，所以也不敢随便乱说。

4. 实验⽰例

给出从 EEG 信号中提取功率谱特征并分类的简单范例。

4.1 实验数据

本⽂选⽤的实验数据为的，使⽤的数据为 sp1s_aa_1000Hz.mat。

这个数据集中，受试者坐在⼀张椅⼦上，⼿臂放在桌⼦上，⼿指放在电脑键盘的标准打字位置。被试需要⽤⾷指和⼩指依照⾃⼰选择的顺序按相应的键。实验的⽬标是预测按键前130毫秒⼿指运动的⽅向（左 OR 右）。

在 matlab 中导⼊数据。

%%导⼊数据

%1000 Hz 记录了500 ms

load('sp1s_aa_1000Hz.mat');

%采样率1000 Hz

srate =1000;

[frames, channels, epochs]=size(x_train);

rightwards =sum(y_train);

leftwards =length(y_train)- rightwards;

fprintf('⼀共有 %d 个训练样本，其中往左运动有 %d 个，往右运动有 %d 个\n',...

length(y_train), leftwards, rightwards);

⼀共有 316 个训练样本，其中往左运动有 159 个，往右运动有 157 个

4.2 提取特征

我们使⽤ welch 法来提取功率谱密度，利⽤ pwelch 函数计算功率谱，使⽤ bandpower 函数可以提取特定频段的功率信息，所以分别提δθαβ

取、、、节律的功率。最后取各通道平均功率的前12个点（根据 f 来看，前 12 个点基本覆盖了 0到 40Hz 的频带）

%%提取 PSD 特征

function[power_features]=ExtractPowerSpectralFeature(eeg_data, srate)

%从 EEG 信号中提取功率谱特征

% Parameters:

% eeg_data:[channels, frames]的 EEG 信号数据

% srate: int，采样率

% Returns:

% eeg_segments:[1, n_features] vector

%%计算各个节律频带的信号功率

[pxx, f]=pwelch(eeg_data,[],[],[], srate);

宋嘉宝

power_delta =bandpower(pxx, f,[0.5,4],'psd');

power_theta =bandpower(pxx, f,[4,8],'psd');

power_alpha =bandpower(pxx, f,[8,14],'psd');

湖北高职学校

power_beta =bandpower(pxx, f,[14,30],'psd');

%求 pxx 在通道维度上的平均值

mean_pxx =mean(pxx,2);

%简单地堆叠起来构成特征（可以⽤更⾼级地⽅法，⽐如考虑通道之间的相关性的⽅法构成特征向量）

power_features =[

power_delta, power_theta,...

power_alpha, power_beta,...

mean_pxx(1:12)';

];

end

然后对每个样本都提取特征，构造⼀个⼆维矩阵 X_train_features。

X_train_features =[];

for i =1:epochs

%取出数据华北油田

eeg_data =squeeze(x_train(:,:, i));

feature =ExtractPowerSpectralFeature(eeg_data, srate);

X_train_features =[X_train_features; feature];

end

%原始的 y_train 是⾏向量，展开成列向量

y_train =y_train(:);

4.3 分类

使⽤ SVM 进⾏简单的分类任务，由于只是简单演⽰，所以不划分训练集、交叉验证集。

%由于只是简单演⽰，所以不划分训练集、交叉验证集

model =fitcsvm(X_train_features, y_train,...

'Standardize',true,'KernelFunction','RBF','KernelScale','auto','verbose',1);

耐热钢焊接y_pred = model.predict(X_train_features);

acc =sum(y_pred == y_train)/length(y_pred);

fprintf('Train Accuracy: %.2f%%\n', acc *100);

结果如下：

|========================================================================================================================== =========|

|========================================================================================================================== =========|

| 0 |active| 316 | 9.968454e-01 | 2.000000e+00 | 1.000000e+00 | 0 | 0.000000e+00 | 0.000000e+00 |

| 350 |active| 316 | 5.175246e-05 | 9.741516e-04 | 5.129944e-04 | 312 | -1.850933e+02 | 5.967449e-16 |

由于 DeltaGradient，收敛时退出活动集。

Train Accuracy: 94.62%

Extra Update (2021.11)

Update On 2021/11/04。

之前写这篇⽂章⼏乎是随⼿写的，⽆论是⽂本还是代码都写得很粗糙。但不曾想这篇⽂章这么受欢迎，以⾄于之前粗糙的地⽅使得不断有朋友来提出疑问。⽐如这样。

我在这⾥再加⼀些解释。

在本⽂的 demo ⾥，我调⽤了 pwelch 函数。

[pxx, f] = pwelch(eeg_data, [], [], [], srate);

这个函数的具体调⽤格式应该是：

[pxx,f] = pwelch(x,window,noverlap,f,fs)

本文发布于:2024-09-21 14:38:35，感谢您对本站的认可！

本文链接：https://www.17tex.com/xueshu/400490.html

上一篇：双边噪声功率谱密度公式

下一篇：离散信号的功率谱密度

标签：功率信号函数提取密度特征平均功率过程

留言与评论（共有 0 条评论）