點擊上面"腦機接口社區"關注我們
更多技術干貨第一時間送達
腦電信號是一種非平穩的隨機信號，一般而言隨機信號的持續時間是無限長的，因此隨機信號的總能量是無限的，而隨機過程的任意一個樣本函數都不滿足絕對可積條件，所以其傅里葉變換不存在。
不過，盡管隨機信號的總能量是無限的，但其平均功率卻是有限的，因此，要對隨機信號的頻域進行分析，應從功率譜出發進行研究才有意義。正因如此，在研究中經常使用功率譜密度(Power spectral density, PSD)來分析腦電信號的頻域特性。
本文簡單的演示了對腦電信號提取功率譜密度特征然后分類的基本流程。希望對那些尚未入門、面對 BCI 任務不知所措的新手能有一點啟發。
2. 功率譜密度理論基礎簡述功率譜密度描是對隨機變量均方值的量度，是單位頻率的平均功率量綱。對功率譜在頻域上積分就可以得到信號的平均功率。
功率譜密度 是一個以頻率 為自變量的映射， 反映了在頻率成分 上信號有多少功率。
我們假定一個隨機過程 ，并定義一個截斷閾值 ，隨機過程 的截斷過程 就可以定義為
則該隨機過程的能量可定義為
對能量函數求導，就可以獲得平均功率。
根據 Parseval 定理(即能量從時域角度和頻域角度來看都是相等的)可得：
這里 是 經過傅里葉變換后的形式。由于隨機過程 被限定在了一個有限的時間區間 之間，所以對隨機過程的傅里葉變換不再受限。另外我們還需要注意到， 是一個隨機變量，因此為了得到最終總體的平均功率，還需要求取隨機變量的期望值。
由此，通過求取 時的極限，就可以得到原始隨機過程的平均功率 。
將式中被積函數單獨提取出來，定義為 :
這樣一來，平均功率 可以表示為 。通過這種定義方式，函數 可以表征每一個最小極限單位的頻率分量所擁有的功率大小，因此我們把 稱為功率譜密度。
3. Matlab 中?PSD 函數的使用功率譜密度的估計方法有很多。總體來講可以分為兩大類：傳統的非參數方法，和現代的參數方法。
本節不對理論知識做詳細的敘述，感興趣的可以深入查閱文獻，這里只介紹一下有哪些方法，以及他們在 matlab 當中的使用。
3.1 傳統非參數方法估計 PSD
最簡單的方法是周期圖法，先對信號做 FFT 變換，然后求平方，periodogram 函數實現了這個功能。不過周期圖法估計的方差隨采樣點數N的增加而增加，不是很建議使用。
另一種自相關方法，基于維納辛欽定律：信號的功率譜估計等于該信號自相關函數的離散DTFT，不過我沒有在 matlab 里找到對應的函數，如果有知道的朋友請告訴我一下。
最常用的函數是 pwelch 函數，利用 welch 方法來求 PSD，這也是最推薦使用的。
3.2 參數方法估計 PSD
包括 pconv、pburg、pyulear 等幾個方法。
這些方法我沒用過，所以也不敢隨便亂說。
4. 實驗示例給出從 EEG 信號中提取功率譜特征并分類的簡單范例。
4.1 實驗數據
本文選用的實驗數據為BCI Competition Ⅱ的任務四，使用的數據為 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);
結果如下：
|===================================================================================================================================|
|???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%
文章來源于網絡，僅用于學術交流，不用于商業行為
若有侵權及疑問，請后臺留言，管理員即時刪侵！
更多閱讀
EEG偽影類型詳解和過濾工具的匯總(一)
你真的了解腦機接口技術嗎？
清華張鈸院士專刊文章：邁向第三代人工智能(全文收錄)
腦機接口拼寫器是否真的安全？華中科技大學研究團隊對此做了相關研究
腦機接口和卷積神經網絡的初學指南(一)
腦電數據處理分析教程匯總(eeglab, mne-python)
P300腦機接口及數據集處理
快速入門腦機接口：BCI基礎(一)
如何快速找到腦機接口社區的歷史文章？
腦機接口BCI學習交流QQ群：515148456
微信群請掃碼添加，Rose拉你進群
(請務必填寫備注，eg. 姓名+單位+專業/領域/行業)
長按關注我們
歡迎點個在看鼓勵一下
)
