实战3 - 面向康复工程的脑电信号分析和判别模型

题目来源于2020年中国研究生数学建模竞赛C题——面向康复工程的脑电信号分析和判别模型。

在本次实战的数据分析过程中，涉及以下技术内容：
（1）分类模型的训练与评价
（2）数据标准化
（3）带通滤波
（4）重采样（数据增强）
（5）模型在新数据集上预测

我先将字符矩阵的标识符、行/列的标识符存在两个矩阵里，并保存为 C_HWB_2020.mat，便于后续程序随时调用。
矩阵BiaoShiFu_all如下，

矩阵HangLie_BSF_all如下，

LJ_train_data='附件1-P300脑机接口数据S1S1_train_data.xlsx';
LJ_train_event='附件1-P300脑机接口数据S1S1_train_event.xlsx';
sheets_1=sheetnames(LJ_train_data);
sheets_2=sheetnames(LJ_train_event);

有些算法对标准化敏感，有些则不敏感。但因为后边我们要对未知字符进行预测，这些字符中有些字符在训练集中未出现过，为了使训练模型在新数据有较好的预测性能，我采用了数据标准化。而且这种标准化是将12个sheet的数据放到一块组成matrix1_all，然后对其每一列进行标准化，并保留标准化的输出物ps，用ps来对训练集（测试集）数据以及新数据进行标准化，这样便统一了标准化的标准。

matrix1_all=[];
for i=1:length(sheets_1)
    table_tem1=readtable(LJ_train_data,'Sheet',sheets_1(i));
    matrix_1=table_tem1{:,:};
    matrix1_all=[matrix1_all;matrix_1];
end

因为P300信号的频率在0.1~20Hz之间，此处我编写了带通滤波函数，滤去该频率区间外的波进行带通滤波，

for k=1:size(matrix1_all,2)
    matrix1_all(:,k)=bandfilter(matrix1_all(:,k));
end
%标准化
[~,ps]=mapstd(matrix1_all'); 
matrix1_all=[];

滤波函数如下，

function y=bandfilter(x)
filorder = 20;  %滤波器阶数
cutf1 =0.001;  %滤波频率1
cutf2 = 23;  %滤波频率2
samplerate=250;  %采样频率
d = designfilt('bandpassfir','FilterOrder',filorder, ...
         'CutoffFrequency1',cutf1,'CutoffFrequency2',cutf2, ...
         'SampleRate',samplerate);
y = filtfilt(d,x);

end

P300信号位于刺激发生后的600ms内，因为采样频率是250Hz，所以P300信号也就是位于150个采样点内。注意每个训练样本是150×20的矩阵，需要将其拉成一行，

warning('off');
Xdata =[];
Ydata=[];
for i=1:length(sheets_1)
    
    table_tem1=readtable(LJ_train_data,'Sheet',sheets_1(i));
    table_tem2=readtable(LJ_train_event,'Sheet',sheets_2(i));
    matrix_1=table_tem1{:,:};
    matrix_2=table_tem2{:,:};
    %
    %进行带通滤波
    for k=1:size(matrix_1,2)
        matrix_1(:,k)=bandfilter(matrix_1(:,k));
    end
    %
    %标准化
    bz=mapstd('apply',matrix_1',ps);
    matrix_1=bz';
    
    for j=1:length(sel_index)
        
        start_index=matrix_2(sel_index(j),2);
        sel_data=matrix_1(start_index:(start_index+150),:);
        %拉成一行
        Xdata = [Xdata;reshape(sel_data,1,[])];
        hangxie=HangLie_BSF_all(BiaoShiFu_all==matrix_2(1,1),:);
        Ydata=[Ydata;categorical(ismember(matrix_2(sel_index(j),1),hangxie))];
        
    end
    
end

因为正样本（包含P300信号）只有120个，而负样本有600个，正负样本不均衡，需要对正样本进行重采样，使其数据和负样本相当。此处我采用随机抽取两个正样本数据求平均值的方式，来增加正样本数量。这样既使得正负样本比例达到平衡，又增加了新的正样本数据，相对于通过重复采样来增加正样本的方法来说， 这种方法对后续识别新的字符更有利。

d_Ydata=double(Ydata);
hang_2=find(d_Ydata==2);
Xdata_chong =Xdata;
Ydata_chong =Ydata;
rng(1);
for k=1:480
    r = randi([1,length(hang_2)],2,1);
    Xdata_chong =[Xdata_chong;mean(Xdata(hang_2(r),:))];
    Ydata_chong =[Ydata_chong;Ydata(hang_2(r(1)),:)];
end

直接参考help文档里关于划分训练集和测试集的代码，

%合并特征数据与目标数据
Table_data_S1=array2table(Xdata_chong);
Table_data_S1.Ydata=Ydata_chong;

%划分训练集和测试集
rng(1); % For reproducibility of the data partition
c = cvpartition(Ydata_chong,'Holdout',0.2);
trainingIdx = training(c); % Training set indices
S1Train = Table_data_S1(trainingIdx,:);
testIdx = test(c); % Test set indices
S1Test = Table_data_S1(testIdx,:);

通过分类学习Classification Learner app，对S1Train进行分类，比对下各种分类器的效果，见下图，



发现线性模型、SVM、KNN和集成模型精确度超多80%，接下来我们在S1数据上使用fitcauto函数对这四种分别进行训练，来详细比对各学习器的性能，
线性模型训练结果如下，

%训练分类模型
options = struct('UseParallel',true);
Mdl = fitcauto(S1Train,'Ydata','HyperparameterOptimizationOptions',options,"Learners",'linear');

%评估分类模型性能
testAccuracy = 1 - loss(Mdl,S1Test,'Ydata')
confusionchart(S1Test.Ydata,predict(Mdl,S1Test))

SVM训练结果如下，

%训练分类模型
options = struct('UseParallel',true);
Mdl = fitcauto(S1Train,'Ydata','HyperparameterOptimizationOptions',options,"Learners",'svm');

%评估分类模型性能
testAccuracy = 1 - loss(Mdl,S1Test,'Ydata')
confusionchart(S1Test.Ydata,predict(Mdl,S1Test))

KNN训练结果如下，

%训练分类模型
options = struct('UseParallel',true);
Mdl = fitcauto(S1Train,'Ydata','HyperparameterOptimizationOptions',options,"Learners",'knn');

%评估分类模型性能
testAccuracy = 1 - loss(Mdl,S1Test,'Ydata')
confusionchart(S1Test.Ydata,predict(Mdl,S1Test))

集成模型如下，

%训练分类模型
options = struct('UseParallel',true);
Mdl = fitcauto(S1Train,'Ydata','HyperparameterOptimizationOptions',options,"Learners",'knn');

%评估分类模型性能
testAccuracy = 1 - loss(Mdl,S1Test,'Ydata')
confusionchart(S1Test.Ydata,predict(Mdl,S1Test))

对比发现，SVM性能最好。

使用训练的SVM模型分别对S1~S5的新数据进行预测，

LJ_test_data='附件1-P300脑机接口数据S1S1_test_data.xlsx';
LJ_test_event='附件1-P300脑机接口数据S1S1_test_event.xlsx';
%取出每个文件里包含sheet的名称
sheets_1=sheetnames(LJ_test_data);
sheets_2=sheetnames(LJ_test_event);

%预测新字符
warning('off');
varname=Table_data_S1.Properties.VariableNames;
cell_predict={};%存储每个被试的预测结果
for i=1:length(sheets_1)
    %取出每个被试的data和event数据
    Xdata =[];
    Y_BiaoQian=[];
    table_tem1=readtable(LJ_test_data,'Sheet',sheets_1(i));
    table_tem2=readtable(LJ_test_event,'Sheet',sheets_2(i));
    matrix_1=table_tem1{:,:};
    matrix_2=table_tem2{:,:};
    %
    %进行带通滤波
    for k=1:size(matrix_1,2)
        matrix_1(:,k)=bandfilter(matrix_1(:,k));
    end
    %
    %按训练数据进行标准化
    bz=mapstd('apply',matrix_1',ps);
    matrix_1=bz';
    
    %取出每次闪烁对应的数据矩阵，并拉直成一行
    for j=1:length(sel_index)
        
        start_index=matrix_2(sel_index(j),2);
        sel_data=matrix_1(start_index:(start_index+150),:);
        Xdata = [Xdata;reshape(sel_data,1,[])];
        %Y_BiaoQian存储每次闪烁对应的行/列的标识符
        Y_BiaoQian=[Y_BiaoQian;matrix_2(sel_index(j),1)];
        
    end
    %求出预测结果为"true"对应的行/列的标识符
    Table_predict=array2table(Xdata);
    Table_predict.Properties.VariableNames=varname(1:end-1);
    Y_predict=predict(Mdl,Table_predict);
    cell_predict{i}=Y_BiaoQian(Y_predict==categorical("true"));
end

将每个被试对待识别字符的识别结果汇总（被试S2和S3只有前9个字符的数据），分别取行标识符里频率最高和列标识符里频率最高的标识符，由两者来确定待识别字符。
第一个待识别字符的行/列的标识符汇总频率表如下，最高频率的是3和7，对应字符“M”。

第二个待识别字符的行/列的标识符汇总频率表如下，最高频率的是1和12，对应字符“F”。

第三个待识别字符的行/列的标识符汇总频率表如下，最高频率的是6和7，对应字符“5”。

第四个待识别字符的行/列的标识符汇总频率表如下，最高频率的是5和10，对应字符“2”。

第五个待识别字符的行/列的标识符汇总频率表如下，最高频率的是2和9，对应字符“I”。

第六个待识别字符的行/列的标识符汇总频率表如下，最高频率的是4和8，对应字符“T”。

第七个待识别字符的行/列的标识符汇总频率表如下，最高频率的是2和11，对应字符“K”。

第八个待识别字符的行/列的标识符汇总频率表如下，最高频率的是4和12，对应字符“X”。

第九个待识别字符的行/列的标识符汇总频率表如下，最高频率的是1和7，对应字符“A”。

第十个待识别字符的行/列的标识符汇总频率表如下，最高频率的是6和12，对应字符“0”。

问题二：对于每个被试，可使用NCA算法求出151×20个特征的重要度，在使用reshape重新排成151×20的矩阵，按行求和，即得出每个通道的重要度得分。
问题三：使用一半数据作为有标签数据，另一半作为无标签数据，然后首先使用有标签数据训练一个模型SVMModel，使用如下函数去预测无标签数据，

[label,score] = predict(SVMModel,X)

每个无标签数据都有一个得分score（score为两列，正值所在的那一列即为预测的类别），我们把得分最高的那个无标签数据赋上类别标签，加入到训练集中，重新训练模型SVMModel，反复使用这种方法，就可以使用少量的有标签数据，逐渐给无标签数据打上标签。
问题四：最后一问是上述三问所得方法的实践题，不再赘述。

（1）训练的模型较好，但放到新数据上预测时变得很差，大概率是没使用训练集数据的标准化方法来标准化新数据。否则，那就是新数据和训练数据有结构性的重要差异，训练得到的经验无法推广到差异性太大的新数据上。
（2）滤波对提升模型的性能有一定的功效。
（3）数据增强时，生成新数据比重复采样的方法要好一些。