基于分类任务的信号(EEG)处理
更多技术干货第一时间送达
本文由网友Jon_Snow_Stark授权分享
对于很多从没做过信号处理以及分类的小伙伴来说，面对信号数据真的是不知道该如何下手。既然大家来看这篇博文，我相信，大家与我一样，其中的痛楚无需多言。下面我就与大家分享一下，这段时间我对这一问题的感悟（新手上路，如有错误，请大家批评指正）。
首先我们需要对这一问题有一个宏观的认识。我们采集到的原始脑电信号包含了很多的噪声与干扰，这些在分类中都会影响分类性能，那么我们就需要先对信号去噪去干扰。经过这一步的处理后得到的就是较为纯净的信号了，接下来就要对信号提取特征了，之后利用分类器对信号的特征进行学习，构建模型，然后进行分类。
那么该如何提取特征呢？在我处理的信号中，我需要分delta(1-3hz)、theta(4_7hz)、alpha (8-13hz)、beta(14-30hz)四个频带，因此我利用滤波器，先分别得到这四个频带的信号，然后再在这四个频带上提取特征。提取特征后我们利用特征矩阵进行分类任务。
下面我们来对基于脑电信号的分类任务进行分析：
1、我们首先采集到的是脑电的原始信号，例如有BDF格式、CNT格式等等，我们利用matlab处理不能直接读取这些格式的文件，那么我们就需要对这些格式的文件进行处理，将其保存为.mat格式的文件，这样matlab就能对其直接进行运算了。处理这些格式的文件可以使用eeglab工具包，这是一个图形化界面的处理工具包，当然也可以调用相关函数进行处理。关于eeglab工具包的详细信息就不在这里阐述了，大家可以参考eeglab工具包的官网进行下载和了解。
也可以查看《收藏学习！EEGLAB处理脑电教程17篇汇总》、《收藏学习！EEGLAB 视频教程 汇总》等关于EEGLAB方面的教程。
2、对原始脑电信号的格式文件处理后得到相关的.mat格式文件，此时便可以开始对我们的信号进行处理了。在这里假设我们采集了32通道的脑电信号，采样率为1000Hz，采集了10分钟，转化后的.mat格式文件的shape为32*600000，那么行向量就表示一个通道所采集的脑电信号，每一个通道有10min*60s*1000Hz个样本点。我们可以通过eeglab去除一些不好的样本点，也可以利用一些算法进行去噪去干扰，例如去除工频干扰、利用滤波器滤波、去除眼电等，可以参考EEG预处理思维导图。
图片来源于网络(来源CSDN用户:yjinyyzyq)
3、提取特征时可以按通道提取，即将32个通道分开进行提取，在这里我们提取频域信号功率谱密度。因为我需要分析delta、theta、alpha、bete四个频带，因此我先对当前处理的通道信号提取四个频带的信号。然后在每个频带上求其对应的功率谱密度。如果一个标签对应1s的数据，那么我们需要将信号按1s进行分段，即每1000个采样点求一个PSD值；然后将四个频带分别求出的特征值进行排列。此时一个通道（shape为1*600000）的数据就会变为shape为4*600的特征数据集，对32个通道的数据提取特征后，如果排列成3为则shape为32*4*600，如果排列成二维这共有32*4=128行，则shape为128*600，此时其对应的标签的shape应该为1*600。
这里给出一个简要的代码：
for i =1:channals %channals为通道数 delta = butter_bandpass_filter(data(:,i), 1, 3, frequency,3); theta = butter_bandpass_filter(data(:,i), 4, 7, frequency,3); alpha = butter_bandpass_filter(data(:,i), 8, 13, frequency,3); beta = butter_bandpass_filter(data(:,i), 14, 30, frequency,3);
for j =1:sample % sample为采样点数 psd_delta = [psd_delta;compute_psd(delta(((j-1)*frequency*time_win+1):j*frequency*time_win, 1))]; ... end.........end功率谱估计
功率谱估计是频域分析的主要手段。它的意义在于把幅度随时间变化的脑电波变换为脑电功率随频率变化的谱图，从而可直观地观察到脑电节律的分布与变换情况。
经典功率谱估计采用的是传统傅里叶变换分析方法(又称线性谱估计)，主要包括直接法(又称周期图法)和间接法(又称自相关法)两种。
周期图法直接对观测数据进行快速傅里叶变换，得到功率谱。
自相关法先估计自相关函数，再计算功率谱。
由于周期图法直接使用FFT快速计算，所以应用比较广泛。
经典谱估计优点是计算效率高，缺点是频率分辨率低，常用于频率分辨率要求不高的场合。
本篇文章主要介绍直接法(周期图法)的原理以及实现
直接法又称周期图法，这种功率谱也具有周期性，常称为周期图。
周期图是信号功率谱的一个有偏估值，而且当信号序列的长度增大到无穷时，估值的方差不趋于零。
[说明：随着所取的信号序列长度的不同，所得到的周期图也不同，这种现象称为随机起伏。由于随机起伏大，使用周期图不能得到比较稳定的估值]
功率谱直接法案例说明
clear;%采样频率Fs=1000;n=0:1/Fs:1;%生成含有噪声的信号序列xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
%设置矩形窗window=boxcar(length(xn));
nfft=1024;%直接法[Pxx,f]=periodogram(xn,window,nfft,Fs);plot(f,10*log10(Pxx));xlabel('频率/Hz');ylabel('功率/dB')4、此时就可以将处理后的数据和标签输入到分类器中进行训练模型了。例如使用SVM进行分类，我们需要将特征矩阵转换为样本*维度的shape形式，那么就对特征矩阵和标枪向量同时进行转置后输入到svmtrain函数中，再利用predict函数进行预测。