共空间模式 CSP 和滤波器组共空间模式 FBCSP 的 MATLAB 实现

共有四个函数

cspfunc：针对两类脑电和标签提取csp特征并获取滤波器

cspfilt：使用获取的csp滤波器滤波

fbcspfunc：针对两类脑电和标签提取fbcsp特征并获取滤波器（fbcsp 就是分频段组合的 csp，可以是csps？）

fbcspfilt：使用获取的fbcsp滤波器滤波

需要注意的是，必须先进行 func 计算获取滤波器参数，才能进行 filt

cspfunc 函数（主要的函数，借鉴上边博主的demo，大家如果觉得这篇有用也请一定去支持一下他）：

function [cspfeature cspw] = cspfunc(eeg,label,m)
% eeg:分段后的三维脑电（已滤波），采样点 * 通道 * 实验次数
% label：脑电的标签，列向量，值只有 1 和 2
% m：csp 参数

% cspfeature:提取到的fbcsp 特征，主要用于训练，如果是测试，需要用 fbcspfilt 函数滤波
% cspw:fbcsp 滤波器

%check and initializations
EEG_Channels = size(eeg,2); % 通道数目
EEG_Trials = size(eeg,3); % 实验次数
classLabels = unique(label);% Return non-repeating values
EEG_Classes = length(classLabels); % 标签类个数

covMatrix = cell(EEG_Classes,1); % 协方差矩阵
% Computing the normalized covariance matrices for each trial
trialCov = zeros(EEG_Channels,EEG_Channels,EEG_Trials);
for i = 1:EEG_Trials
    E = eeg(:,:,i)';
    EE = E*E';
    trialCov(:,:,i) = EE./trace(EE);  % 计算协方差矩阵
end
clear E;
clear EE;
% 计算每一类样本数据的空间协方差之和
for i = 1:EEG_Classes
    covMatrix{i} = mean(trialCov(:,:,label == classLabels(i)),3);
end
% 计算两类数据的空间协方差之和
covTotal = covMatrix{1} + covMatrix{2};
% 计算特征向量和特征矩阵
[Uc,Dt] = eig(covTotal);
% 特征值要降序排列
eigenvalues = diag(Dt);
[eigenvalues,egIndex] = sort(eigenvalues, 'descend');% 降序
Ut = Uc(:,egIndex);
% 矩阵白化
P = diag(sqrt(1./eigenvalues))*Ut';
% 矩阵P作用求公共特征向量transformedCov1 
transformedCov1 = P*covMatrix{1}*P';
%计算公共特征向量transformedCov1的特征向量和特征矩阵
[U1,D1] = eig(transformedCov1);
eigenvalues = diag(D1);
[eigenvalues,egIndex] = sort(eigenvalues, 'descend');% 降序排列
U1 = U1(:, egIndex);
% 计算投影矩阵W
CSPMatrix = U1' * P;
% 计算特征矩阵
FilterPairs = m;       % CSP特征选择参数m    CSP特征为2*m个
features_train = zeros(EEG_Trials, 2*FilterPairs+1);

Filter = CSPMatrix([1:FilterPairs (end-FilterPairs+1):end],:);
cspw = Filter;

%extracting the CSP features from each trial
for t=1:EEG_Trials    
    % projecting the data onto the CSP filters    
    projectedTrial_train = Filter * eeg(:,:,t)';    
    
    % generating the features as the log variance of the projected signals
    variances_train = var(projectedTrial_train,0,2);  
    
    for f=1:length(variances_train)
        features_train(t,f) = log(variances_train(f));
        % features_train(t,f) = log(variances_train(f)/sum(variances_train));   %修改后对应公式
    end
    
end
cspfeature = features_train(:,1:2*m);

end

cspfilt函数：

function cspfiltedeeg = cspfilt(eeg,cspw)
%使用获取到的 csp 滤波器滤波，获取方式见 cspfunc

Filter = cspw;
m = size(cspw,1)/2;

features_train = zeros(size(eeg,3), 2*m+1);

%extracting the CSP features from each trial
for t=1:size(eeg,3)    
    % projecting the data onto the CSP filters    
    projectedTrial_train = Filter * eeg(:,:,t)';    
    
    % generating the features as the log variance of the projected signals
    variances_train = var(projectedTrial_train,0,2);  
    
    for f=1:length(variances_train)
        features_train(t,f) = log(variances_train(f));
        % features_train(t,f) = log(variances_train(f)/sum(variances_train));   %修改后对应公式
    end
    
end
cspfiltedeeg = features_train(:,1:2*m);

end

fbcspfunc 函数：

function [fbcspfeature,fbcspw] = fbcspfunc(eeg,label,m,bands,fs,filterorder)
% 调用函数 cspfunc，针对 bands 中的每个频带进行 csp，然后拼接特征
% eeg:分段后的三维脑电（已滤波），采样点 * 通道 * 实验次数
% label：脑电的标签，列向量，值只有 1 和 2
% m：csp 参数
% bands:不同的频带，如[1 4;4 7;8 13]
% fs:采样率
% filterorder:滤波器阶数

% fbcspfeature:提取到的fbcsp 特征，主要用于训练，如果是测试，需要用 fbcspfilt 函数滤波
% fbcspw:fbcsp 滤波器

fbcspfeature = [];
fbcspw=[];

for i = 1:size(bands,1)
    udfilter=designfilt('bandpassiir','FilterOrder',filterorder, ...
    'HalfPowerFrequency1',bands(i,1),'HalfPowerFrequency2',bands(i,2), ...
    'SampleRate',fs);
    filtedeeg = zeros(size(eeg));
    for trial = 1:size(eeg,3)
        for ch = 1:size(eeg,2)
            filtedeeg(:,ch,trial)=filter(udfilter,eeg(:,ch,trial));
        end
    end

    [cspfeature cspw] = cspfunc(filtedeeg,label,m);

    fbcspfeature = cat(2,fbcspfeature,cspfeature);

    fbcspw = cat(3,fbcspw,cspw);

end

end

fbcspfilt 函数：

function fbcspfiltedeeg = fbcspfilt(eeg,fbcspw,bands,fs,filterorder)
%使用获取到的 csp 滤波器滤波，获取方式见 fbcspfunc


Filter = fbcspw;
m = size(fbcspw,1)/2;

fbcspfiltedeeg = [];
for i=1:size(fbcspw,3)
    udfilter=designfilt('bandpassiir','FilterOrder',filterorder, ...
    'HalfPowerFrequency1',bands(i,1),'HalfPowerFrequency2',bands(i,2), ...
    'SampleRate',fs);
    filtedeeg = zeros(size(eeg));
    for trial = 1:size(eeg,3)
        for ch = 1:size(eeg,2)
            filtedeeg(:,ch,trial)=filter(udfilter,eeg(:,ch,trial));
        end
    end

    features_train = zeros(size(filtedeeg,3), 2*m+1);
    %extracting the CSP features from each trial
    for t=1:size(filtedeeg,3)
        % projecting the data onto the CSP filters
        projectedTrial_train = Filter(:,:,i) * filtedeeg(:,:,t)';

        % generating the features as the log variance of the projected signals
        variances_train = var(projectedTrial_train,0,2);

        for f=1:length(variances_train)
            features_train(t,f) = log(variances_train(f));
            % features_train(t,f) = log(variances_train(f)/sum(variances_train));   %修改后对应公式
        end

    end
    cspfiltedeeg = features_train(:,1:2*m);
    fbcspfiltedeeg = cat(2,fbcspfiltedeeg,cspfiltedeeg);
end
end

按照格式使用函数就行

如：

[fbcspfeature fbcspw]= fbcspfunc(eeg,label,2,[1 4;4 7;8 13],512,4);
fbcspfiltedeeg = fbcspfilt(eeg,fbcspw,[1 4;4 7;8 13],512,4);