Matlab FBCSP實作

最近在進行運動想象分類的時候想利用FBCSP(Filter Bank Common Sptial Pattern)進行特征提取，但在網上搜索一番都沒有找到實作FBCSP的matlab代碼，所以打算自己用matlab編碼，

FBCSP原理

FBCSP是基于CSP(Common Sptial Pattern)的空域提取方法通過對頻帶分片加之特征選擇演算法實作的，
下面為FBCSP的原理圖：

如上圖所示，其中Stage 1為對頻帶進行切片；
Stage 2對切邊后的各個子頻帶進行CSP濾波；
Stage 3對進行CSP濾波后的特征進行特征選擇；
Stage 4對特征選擇后的資料放到分類器里面分類，

測驗資料

用到的測驗資料為BCIcomp2002的二分類運動想象資料集，因為官方給出的資料沒有測驗集標簽，所以我把訓練集及其標簽拆分為100:40變為訓練集和測驗集，即前100個trials為訓練集，后40個trials的資料用于測驗，把分割后的資料保存到XY_data.mat中，其范式如下圖所示：

運行結果

因為代碼寫的比較創促，所以打算用參考文獻[1]的互資訊特征選擇演算法還沒有寫，打算以后補充上，分類演算法用matlab的svm工具箱，其中-c設定為2，-g設定為0.1250，也嘗試過用網格搜索法確定-c -g的引數但最后出來的效果似乎不是很好，運行結果如下所示：

代碼

檔案結構

主檔案test.m

// An highlighted block
load XY_data.mat
train_data=x_train_new(128*4:end,:,:);  %訓練集(拋棄前4s)
train_label=y_train_new;                %訓練集標簽
test_data=x_test_new(128*4:end,:,:);    %測驗集(拋棄前4s）
test_label=y_test_new;                  %測驗集標簽
sampleRate=128; %定義采樣率 
[result,proj,classNum]=FBCSP(train_data,train_label,sampleRate,2); %對訓練集資料進行FBCSP
% result=MIfeature_select(result,k); %對所得訓練集資料進行互資訊特征選擇
model=libsvmtrain(train_label,result,'-c 2 -g 0.1250'); %訓練模型
fea=FBCSPOnline(test_data,proj,classNum,sampleRate,2);  %對測驗集資料進行FBCSP
%fea=MIfeature_select    %對測驗集資料進行MIfeature_select
libsvmpredict(test_label,fea,model);    

FBCSP.m函式

%author:mao date:2020-10-8
%programme:the Filter Bank Common Sptial Pattern algorithm
%input:train_data 3維 EEG資料，其中，第一維是采樣點，第二維是通道數量，第三維度是trials大小
%      train_label train_data對應的標簽
%      sampleRate 采樣率
%      m CSP的m引數
%output:features_train 融合后各子頻帶后的特征陣列
%       projMAll 由各子頻帶計算所得的投影矩陣
%       classNum 待分類的類別數量
function [features_train,projM_All,classNum]=FBCSP(train_data,train_label,sampleRate,m)

[q,p,k]=size(train_data);%獲取總的trial次數

%% acquire and combine feature of different frequency bands
features_train=[];      %宣告訓練集csp特征融合陣列
filter_data=zeros(size(train_data));
freq=[8 12 16 20 24 28 32]; %設定子頻帶頻率
classNum=max(train_label);  %獲取類別數量
projM_All=zeros(p,p,max(train_label)*(size(freq,2)-1)); %申請個頻帶的拼接投影矩陣空間
for i=1:size(freq,2)
    lower=freq(i);  %獲取低頻
    if lower==freq(size(freq,2))
        break;
    end
    higher=freq(i+1);%獲取高頻
    %對各子頻帶進行濾波
    filter_tmp=[];
    for j=1:k   %對每個trial進行回圈濾波,filter()函式可以濾波3維資料？
        filter_tmp=filter_param(train_data(:,:,j),lower,higher,sampleRate,4);
        filter_data(:,:,j)=filter_tmp;
    end
    % 計算csp濾波器，用csp濾波器進行特征提取
    projM=cspProjMatrix(filter_data,train_label); %要回圈保存投影矩陣用于在線CSP濾波
    projM_All(:,:,1+(i-1)*classNum:i*classNum)=projM;   %存盤當前頻帶投影矩陣
    feature=[];  %宣告本子頻帶特征矩陣
    for b=1:k    %回圈提取特征
        feature(b,:)=cspFeature(projM,filter_data(:,:,b),m); %第三個引數m不要超過通道數的一半，不然會出現重復特征
    end
    tmp_data=feature;
    features_train=[features_train,tmp_data]; %拼接個自頻帶特征矩陣
end

FBCSPOnline.m函式

%author:mao date:2020-10-8
%programme:the Filter Bank Common Sptial Pattern algorithm
%input:train_data 3維 EEG資料，其中，第一維是采樣點，第二維是通道數量，第三維度是trials大小
%                 2維 EEG資料，其中，第一位是采樣點，第二維是通道數量
%      projMAll 由訓練集計算所得個子頻帶CSP投影矩陣
%      classNum 待分類的類別數
%      sampleRate 采樣率
%      m CSP的m引數
%output:features 融合后各子頻帶后的特征陣列
function features_train=FBCSPOnline(train_data,projMAll,classNum,sampleRate,m)
if ndims(train_data)==3 %輸入EEG資料為3維
    %% acquire and combine feature of different frequency bands
    [q,p,k]=size(train_data);%獲取總的trial次數
    filter_data=zeros(size(train_data));
    features_train=[];      %宣告訓練集csp特征融合陣列
    freq=[8 12 16 20 24 28 32]; %設定子頻帶頻率
    for i=1:size(freq,2)
        lower=freq(i);  %獲取低頻
        if lower==freq(size(freq,2))
            break;
        end
        higher=freq(i+1);%獲取高頻
        %對各子頻帶進行濾波
        filter_tmp=[];
        for j=1:k   %對每個trial進行回圈濾波,matlab中的filter()函式可以濾波3維資料？
            filter_tmp=filter_param(train_data(:,:,j),lower,higher,sampleRate,4);
            filter_data(:,:,j)=filter_tmp;
        end
        feature=[];  %宣告本子頻帶特征矩陣
        for b=1:k    %回圈提取特征
            feature(b,:)=cspFeature(projMAll(:,:,1+(i-1)*classNum:i*classNum),filter_data(:,:,b),m); %第三個引數m不要超過通道數的一半，不然會出現重復特征
        end
        tmp_data=feature;
        features_train=[features_train,tmp_data]; %拼接個自頻帶特征矩陣
    end
else                  
    %% 輸入EEG資料為2維
    features_train=[];      %宣告訓練集csp特征融合陣列
    freq=[8 12 16 20 24 28 32]; %設定子頻帶頻率
    for i=1:size(freq,2)
        lower=freq(i);  %獲取低頻
        if lower==freq(size(freq,2))
            break;
        end
        higher=freq(i+1);%獲取高頻
        %對各子頻帶進行濾波
        filter_data=filter_param(train_data,lower,higher,sampleRate,4);
        feature=[];  %宣告本子頻帶特征矩陣
        feature=cspFeature(projMAll(:,:,1+(i-1)*classNum:i*classNum),filter_data,m); %第三個引數m不要超過通道數的一半，不然會出現重復特征
        tmp_data=feature';
        features_train=[features_train,tmp_data]; %拼接個自頻帶特征矩陣
    end
end

filter_param.m濾波函式

%% 對資料進行濾波
%輸入：data 待濾波EEG資料
%   low        高通濾波引數設定
%   high        低通濾波引數設定
%   sampleRate          采樣率
%   filterorder  butterworth濾波器階數
%回傳：filterdata       濾波后EEG資料
function filterdata=filter_param(data,low,high,sampleRate,filterorder)
%% 設定濾波引數
 filtercutoff = [low*2/sampleRate high*2/sampleRate]; 
[filterParamB, filterParamA] = butter(filterorder,filtercutoff);
 filterdata= filter( filterParamB, filterParamA, data);

cspProjMatrix.m (計算CSP濾波器）

function projM=cspProjMatrix(x,y)
%基于共空間模式演算法計算出一個投影矩陣
%輸入引數：
%        x:3維 EEG資料，其中，第一維是采樣點，第二維是通道數量，第三維度是trials大小
%        y: 一維列向量標簽 范圍是從1到分類數量，長度與x的第三維保持一致
%注意：這里y標簽只能從1開始，往后延，不能用-1 1這種標簽格式

trialNo=length(y); %獲取標簽長度
classNo=max(y);    %獲取標簽類別數量
channelNo=length(x(1,:,1)); %獲取通道數量

for k=1:classNo    %對每一類進行訓練
    N_a=sum(y==k); %number of trials for class k，當前類的trials數量
    N_b=trialNo-N_a;
    R_a=zeros(channelNo,channelNo); %申請[通道數量*通道數量] 方陣大小的空間
    R_b=zeros(channelNo,channelNo);
    for i=1:trialNo
        R=x(:,:,i)'*x(:,:,i); 
        %R=cov(x(:,:,i)); 
        R=R/trace(R);
        if y(i)==k   %當前類
            R_a=R_a+R;
        else         %其他類
            R_b=R_b+R;
        end
    end
    R_a=R_a/N_a;
    R_b=R_b/N_b;
    [V,D]=svd(R_a+R_b);  %矩陣奇異值分解
    W=D^(-0.5)*V';       %P白化矩陣，P矩陣回傳為W
    S_a=W*R_a*W';
    [V,D]=svd(S_a);
    projM(:,:,k)=W'*V;  %投影矩陣， 最后投影矩陣的大小為 [通道數量 通道數量 類別數量] 其中第三維度為每個類的濾波器
end 

cspFeature.m (利用CSP濾波器進行過濾）

function feature=cspFeature(projM,x,m)
%%%%%%通過投影矩陣進行特征提取
% 輸入引數：
%      projM: csp投影矩陣
%          x: 一個時間視窗的2維EEG資料，其中，第一維是采樣點；第二維是通道
%          m: 投影資料矩陣的第一列和最后一列的個數，
% 輸出引數：
%    feature: 從列向量中提取到的特征


classNo=length(projM(1,1,:));  %獲取類別數量
channelNo=size(x,2);           %獲取通道數量
feature=[];                    %宣告特征矩陣
for k=1:classNo                %classNo為類數量
    Z=x*projM(:,:,k); %projected data matrix
    for j=1:m
        feature=[feature; var(Z(:,j)); var(Z(:,channelNo-j+1))];  %var(A) 算矩陣A沒列方差，此時默認是除N-1
        %variances of the first and last m columns(第1和最后m列的方差)
    end
end
feature=log(feature/sum(feature)); 

參考文獻

[1] Ang K K , Chin Z Y , Zhang H , et al. Filter Bank Common Spatial Pattern (FBCSP) in Brain-Computer Interface[C]// IEEE International Joint Conference on Neural Networks. IEEE, 2008.