點(diǎn)擊上面"腦機(jī)接口社區(qū)"關(guān)注我們

更多技術(shù)干貨第一時(shí)間送達(dá)

本文是由CSDN用戶[frostime]授權(quán)分享。主要介紹了腦電信號(hào)提取PSD功率譜密度特征，包括：功率譜密度理論基礎(chǔ)、matlab中PSD函數(shù)的使用介紹以及實(shí)驗(yàn)示例。感謝 frostime！1. 序言

腦電信號(hào)是一種非平穩(wěn)的隨機(jī)信號(hào)，一般而言隨機(jī)信號(hào)的持續(xù)時(shí)間是無限長的，因此隨機(jī)信號(hào)的總能量是無限的，而隨機(jī)過程的任意一個(gè)樣本函數(shù)都不滿足絕對可積條件，所以其傅里葉變換不存在。

不過，盡管隨機(jī)信號(hào)的總能量是無限的，但其平均功率卻是有限的，因此，要對隨機(jī)信號(hào)的頻域進(jìn)行分析，應(yīng)從功率譜出發(fā)進(jìn)行研究才有意義。正因如此，在研究中經(jīng)常使用功率譜密度(Power spectral density, PSD)來分析腦電信號(hào)的頻域特性。

本文簡單的演示了對腦電信號(hào)提取功率譜密度特征然后分類的基本流程。希望對那些尚未入門、面對 BCI 任務(wù)不知所措的新手能有一點(diǎn)啟發(fā)。

2. 功率譜密度理論基礎(chǔ)簡述

功率譜密度描是對隨機(jī)變量均方值的量度，是單位頻率的平均功率量綱。對功率譜在頻域上積分就可以得到信號(hào)的平均功率。

功率譜密度 是一個(gè)以頻率 為自變量的映射， 反映了在頻率成分 上信號(hào)有多少功率。

我們假定一個(gè)隨機(jī)過程 ，并定義一個(gè)截?cái)嚅撝?，隨機(jī)過程 的截?cái)噙^程 就可以定義為

則該隨機(jī)過程的能量可定義為

對能量函數(shù)求導(dǎo)，就可以獲得平均功率。

根據(jù) Parseval 定理(即能量從時(shí)域角度和頻域角度來看都是相等的)可得：

這里 是 經(jīng)過傅里葉變換后的形式。由于隨機(jī)過程 被限定在了一個(gè)有限的時(shí)間區(qū)間 之間，所以對隨機(jī)過程的傅里葉變換不再受限。另外我們還需要注意到， 是一個(gè)隨機(jī)變量，因此為了得到最終總體的平均功率，還需要求取隨機(jī)變量的期望值。

由此，通過求取 時(shí)的極限，就可以得到原始隨機(jī)過程的平均功率 。

將式中被積函數(shù)單獨(dú)提取出來，定義為 :

這樣一來，平均功率 可以表示為 。通過這種定義方式，函數(shù) 可以表征每一個(gè)最小極限單位的頻率分量所擁有的功率大小，因此我們把 稱為功率譜密度。

3. Matlab 中?PSD 函數(shù)的使用

功率譜密度的估計(jì)方法有很多。總體來講可以分為兩大類：傳統(tǒng)的非參數(shù)方法，和現(xiàn)代的參數(shù)方法。

在這里插入圖片描述

本節(jié)不對理論知識(shí)做詳細(xì)的敘述，感興趣的可以深入查閱文獻(xiàn)，這里只介紹一下有哪些方法，以及他們在 matlab 當(dāng)中的使用。

3.1 傳統(tǒng)非參數(shù)方法估計(jì) PSD

最簡單的方法是周期圖法，先對信號(hào)做 FFT 變換，然后求平方，periodogram 函數(shù)實(shí)現(xiàn)了這個(gè)功能。不過周期圖法估計(jì)的方差隨采樣點(diǎn)數(shù)N的增加而增加，不是很建議使用。

另一種自相關(guān)方法，基于維納辛欽定律：信號(hào)的功率譜估計(jì)等于該信號(hào)自相關(guān)函數(shù)的離散DTFT，不過我沒有在 matlab 里找到對應(yīng)的函數(shù)，如果有知道的朋友請告訴我一下。

最常用的函數(shù)是 pwelch 函數(shù)，利用 welch 方法來求 PSD，這也是最推薦使用的。

3.2 參數(shù)方法估計(jì) PSD

包括 pconv、pburg、pyulear 等幾個(gè)方法。

這些方法我沒用過，所以也不敢隨便亂說。

4. 實(shí)驗(yàn)示例

給出從 EEG 信號(hào)中提取功率譜特征并分類的簡單范例。

4.1 實(shí)驗(yàn)數(shù)據(jù)

本文選用的實(shí)驗(yàn)數(shù)據(jù)為BCI Competition Ⅱ的任務(wù)四，使用的數(shù)據(jù)為 sp1s_aa_1000Hz.mat。

實(shí)驗(yàn)使用的數(shù)據(jù)

這個(gè)數(shù)據(jù)集中，受試者坐在一張椅子上，手臂放在桌子上，手指放在電腦鍵盤的標(biāo)準(zhǔn)打字位置。被試需要用食指和小指依照自己選擇的順序按相應(yīng)的鍵。實(shí)驗(yàn)的目標(biāo)是預(yù)測按鍵前130毫秒手指運(yùn)動(dòng)的方向(左 OR 右)。

在 matlab 中導(dǎo)入數(shù)據(jù)。

%% 導(dǎo)入數(shù)據(jù)
% 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 個(gè)訓(xùn)練樣本，其中往左運(yùn)動(dòng)有 %d 個(gè)，往右運(yùn)動(dòng)有 %d 個(gè)\n',...
length(y_train), leftwards, rightwards);
一共有?316?個(gè)訓(xùn)練樣本，其中往左運(yùn)動(dòng)有?159?個(gè)，往右運(yùn)動(dòng)有?157?個(gè)

4.2 提取特征

我們使用 welch 法來提取功率譜密度，利用 pwelch 函數(shù)計(jì)算功率譜，使用 bandpower 函數(shù)可以提取特定頻段的功率信息，所以分別提取 、、、節(jié)律的功率。最后取各通道平均功率的前12個(gè)點(diǎn)(根據(jù) f 來看，前 12 個(gè)點(diǎn)基本覆蓋了 0到 40Hz 的頻帶)

%% 提取 PSD 特征
function [power_features] = ExtractPowerSpectralFeature(eeg_data, srate)
% 從 EEG 信號(hào)中提取功率譜特征
% Parameters:
% eeg_data: [channels, frames] 的 EEG 信號(hào)數(shù)據(jù)
% srate: int, 采樣率
% Returns:
% eeg_segments: [1, n_features] vector

%% 計(jì)算各個(gè)節(jié)律頻帶的信號(hào)功率
[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);
% 簡單地堆疊起來構(gòu)成特征(可以用更高級(jí)地方法，比如考慮通道之間的相關(guān)性的方法構(gòu)成特征向量)
power_features = [
power_delta, power_theta, ...
power_alpha, power_beta, ...
mean_pxx(1:12)';
];

end

然后對每個(gè)樣本都提取特征，構(gòu)造一個(gè)二維矩陣 X_train_features。

X_train_features = [];
for i = 1:epochs
% 取出數(shù)據(jù)
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 進(jìn)行簡單的分類任務(wù)，由于只是簡單演示，所以不劃分訓(xùn)練集、交叉驗(yàn)證集。

% 由于只是簡單演示，所以不劃分訓(xùn)練集、交叉驗(yàn)證集
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);

結(jié)果如下：

|===================================================================================================================================|
|???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，收斂時(shí)退出活動(dòng)集。
Train?Accuracy:?94.62%
作者博客：https://blog.csdn.net/frostime/article/details/106967703

文章來源于網(wǎng)絡(luò)，僅用于學(xué)術(shù)交流，不用于商業(yè)行為

若有侵權(quán)及疑問，請后臺(tái)留言，管理員即時(shí)刪侵！

更多閱讀

EEG偽影類型詳解和過濾工具的匯總(一)

你真的了解腦機(jī)接口技術(shù)嗎？

清華張鈸院士專刊文章：邁向第三代人工智能(全文收錄)

腦機(jī)接口拼寫器是否真的安全？華中科技大學(xué)研究團(tuán)隊(duì)對此做了相關(guān)研究

腦機(jī)接口和卷積神經(jīng)網(wǎng)絡(luò)的初學(xué)指南(一)

腦電數(shù)據(jù)處理分析教程匯總(eeglab, mne-python)

P300腦機(jī)接口及數(shù)據(jù)集處理

快速入門腦機(jī)接口：BCI基礎(chǔ)(一)

如何快速找到腦機(jī)接口社區(qū)的歷史文章？

腦機(jī)接口BCI學(xué)習(xí)交流QQ群：515148456

微信群請掃碼添加，Rose拉你進(jìn)群

(請務(wù)必填寫備注，eg. 姓名+單位+專業(yè)/領(lǐng)域/行業(yè))

長按關(guān)注我們

歡迎點(diǎn)個(gè)在看鼓勵(lì)一下

總結(jié)

以上是生活随笔為你收集整理的tensorflow提取mel谱特征_【脑电信号分类】脑电信号提取PSD功率谱密度特征的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。