matlab emd功率谱密度,【脑电信 分类】脑电信 提取PSD功率谱密度特征

脑电信 是一种非平稳的随机信 ，一般而言随机信 的持续时间是无限长的，因此随机信 的总能量是无限的，而随机过程的任意一个样本函数都不满足绝对可积条件，所以其傅里叶变换不存在。

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

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

功率谱密度描是对随机变量均方值的量度，是单位频率的平均功率量纲。对功率谱在频域上积分就可以得到信 的平均功率。

功率谱密度 是一个以频率 为自变量的映射， 反映了在频率成分 上信 有多少功率。

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

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

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

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

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

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

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

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

3. Matlab 中 PSD 函数的使用

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

实验使用的数据

这个数据集中，受试者坐在一张椅子上，手臂放在桌子上，手指放在电脑键盘的标准打字位置。被试需要用食指和小指依照自己选择的顺序按相应的键。实验的目标是预测按键前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);

结果如下：

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

|   Iteration  | Set  |   Set Size   |  Feasibility  |     Delta     |      KKT      |  Number of   |   Objective   |   Constraint  |

|              |      |              |      Gap      |    Gradient   |   Violation   |  Supp. Vec.  |               |   Violation   |

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

|            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%

相关资源：Yalefree雅乐简谱打谱软件_打谱软件-WindowsServer工具类资源…