基于脑电的特征频段预处理
本文对一段脑电信号进行预处理，包括工频干扰消除、基线漂移消除，对预处理后的脑电信号进行频谱分析，分别提取theta、delta、alpha、beta、gamma、piper节律信息，并分析各特定频带信号的时域、能量等特征。
软件环境及数据背景：
MATLAB2016a，data.mat。
本数据采样频率为500Hz，采集电极部位：Fp1-F3、F3-C3、C3-P3、P3-O1、O1-T5、T5-T3、T3-F7、F7-Fp1、Fz-F4、F4-C4 、C4-P4、P4-O2、O2-T6 、T6-T4、T4-F8、F8-Fz。
16通道的EEG数据通过脑电监测系统(Nicolet Monitor，U.S.A)进行记录，系统采样率为500Hz，灵敏度为400-800uV/cm，监测系统对记录到的脑电数据自动进行0.15Hz-50Hz的预滤波处理。
data.mat中有两个被试：zcy和zly。其中，zcy_0和zly_0是给药前脑电，zcy_1和zly_1是给药1小时脑电，zcy_2和zly_2是给药2小时脑电，zcy_3和zly_3是给药3小时脑电。
预处理：
在脑电分析之前，一般会对信号去除基线漂移，然后去除工频干扰，再根据所需要的信号的范围进行带通滤波，由于采集系统记录数据时已经进行了0.15Hz-50Hz的滤波处理，包含了人类脑电活动频率，故不再做带通滤波。
去除基线漂移
关于脑电基线漂移的官方定义没有找到，我的理解是指实验系统噪声随时间的缓慢变化。网上浏览了些方法，大同小异，主要有“detrend函数”、“拟合基线”，“低通滤波”和“小波变换”等方法。但是MATLAB自带detrend函数只能用于消除线性趋势，不适用于消除非线性趋势，很明显这段脑电信号将不会是一直具有线性趋势的；低通滤波没有尝试，因为脑电的δ频带为0.5–3HZ；小波还没有做了解，就先采用了拟合基线的方法，将其封装为movBaseline去除基线漂移函数，代码如下：
function [ y,t,f,Originaled,Foriginaled,movOriginaledBline ] = movBaseline( Original,channel,fs,start,finish )
%movBaseline 函数用于提取原始信号的某个通道一段采样点数据，返回提取后的数据及其单边功率谱，并对其做去除基线处理，返回时域信号
% Original 原始信号
% channel 原始信号采集通道
% fs 原始信号采样频率
% start 原始信号开始的采样点
% finish 原始信号结束的采样点
% t 时间矢量
% f 单边频率矢量
% Originaled 原始信号截取后的一段时域信号
% Foriginaled Originaled的单边功率谱
% y 拟合基线
% movOriginaledBline Originaled去除基线后的时域信号
T=1/fs; %采样周期
L=finish-start; %信号长度
t=(0:L-1)*T; %时间矢量
Originaled=Original(channel,start+1:finish); %获取通道channel中start:finish采样点数据
Fori=fft(Originaled); %原始信号做傅里叶变换
P2FOriginaled=Fori/L.*conj(Fori)/L; %双边幅频谱
Foriginaled=P2FOriginaled(1:L/2+1); %单边幅频谱
Foriginaled(2:end-1)=2*Foriginaled(2:end-1); %P1(1)是直流分量
f=fs*(0:(L/2))/L; %单边频率矢量
[p,~,mu] = polyfit(t,Originaled,4); %获得4次多项式系数向量p
y=polyval(p,t,[],mu); %拟合基线
movOriginaledBline=Originaled-y; %去除基线后信号
end代码中的L-1和start+1均是为了保证矢量长度是finish-start，比如1000:2000是1001个点而不是1000；关于函数的输入输出参数已经做了较详细备注，不再解释。挑几点说明一下，可以知道zcy_0数据采样点共有244606个，采样频率是500HZ，故采样大概持续了489秒约8分钟，在MATLAB绘制时波形时，全部显示出来基本就看不清了，故设置采样点起始和终止参数start、finish以截取我们想要查看区间，建议不要超过2000个采样点即4秒。
fft得到是对称的双边谱，双边谱与单边谱本质上是一回事，只是表现形式不同，双边谱在负频率上的值与正频率上的值相等，等于单边谱值的一半，注意这时的量纲并不是功率。
多项式去基线方法，首先使用polyfit函数根据输入信号的特征得到一个拟合曲线的系数向量，拟合原理为最小二乘法，再使用polyval函数得到拟合基线信号，将输入信号减去这个拟合基线信号即可，效果如下：
可以看出拟合基线的0秒和2秒的幅值分别约为-1和4，对应于去除基线后的信号与原始信号的差值。
去除工频干扰
在EEG信号记录过程中，存在两种伪迹，分别为：生理伪迹和非生理伪迹。其中，非生理伪迹主要包括50HZ的工频干扰等，50HZ是因为我们国家的一般用电的交流电工频是50HZ。去除工频干扰使用50HZ陷波器进行滤波，这里陷波器是基于零极点分布原理设计的，参考文献《脑电信号中工频干扰去除的综合研究_王兵》
这部分重点在于得到陷波器幅频响应H
代码如下：
function [ H,W,t,f,movfundFrecy,FmovfundFrecy ] = movFundFrecy( inSignal,f0,fs,k )
%movFundFrecy IIR陷波器去除工频干扰
% inSignal 输入信号
% f0 工频
% fs 采样频率
% K调节参数
% H 复频响应矢量
% W 频率向量
% t 时间矢量
% f 单侧频率矢量
% movfundFrecy 去工频后时域信号
% FmovfundFrecy 去工频后信号的幅频谱
T=1/fs; %采样周期
L=length(inSignal); %信号长度
t=(0:L-1)*T; %时间矢量
f=fs*(0:(L/2))/L; %单侧频率矢量
w0=2*pi*f0/fs; %角频率
B=[1 -2*cos(w0) 1]; %系统传递函数分子
A=[1 -2*k*cos(w0) k*k]; %系统传递函数分母
[H,W]=freqz(B,A,512,fs); %H为幅频响应矢量，W为频率矢量
movfundFrecy=filter(B,A,inSignal); %去工频
FmovfundFrecy=fft(movfundFrecy); %去基线去工频后信号傅里叶变换
P2FmovfundFrecy=FmovfundFrecy/L.*conj(FmovfundFrecy)/L; %双侧幅频谱
FmovfundFrecy=P2FmovfundFrecy(1:L/2+1); %单侧幅频谱
FmovfundFrecy(2:end-1)=2*FmovfundFrecy(2:end-1); %P1(1)是直流分量
end备注应该很清楚了，没什么需要讲的，freqz是滤波器频率响应函数，效果如下：
从幅频谱可以看出来，原始信号在50HZ附近就没有幅值了，看不出来多明显效果，但是总归是一个必要步骤。
求当前通道各节律功率
这部分主要使用到两个函数，pwelch函数用来获得当前输入信号的功率谱，返回功率谱和频率矢量，bandpower函数根据提供的谱密度和所需的频率范围节点，获得该频率段内的平均功率。代码如下：
function [ pxx,fpow,powerFeatures ] = powerSort( inSignal,fs )
%powerSort 求功率谱密度以及各个节律频带的信号功率
% inSignal 输入信号
% fs 采样频率
% pxx 功率谱密度
% fpow 频率向量
% powerFeatures 各节律频带的信号功率组成的数组
%使用 welch 法来提取功率谱密度
[pxx, fpow] = pwelch(inSignal, [], [], [], fs); %对去基线去工频的信号求功率谱密度
%计算各个节律频带的信号平均功率
power_delta = bandpower(pxx, fpow, [0.5, 3], 'psd');
power_theta = bandpower(pxx, fpow, [4, 7], 'psd');
power_alpha = bandpower(pxx, fpow, [8, 13], 'psd');
power_beta = bandpower(pxx, fpow, [14, 30], 'psd');
power_gamma = bandpower(pxx, fpow, [31, 60], 'psd');%从功率谱可以看出50HZ以后就基本没有幅度了
%各节律平均功率数组
powerFeatures=[power_delta,power_theta,power_alpha,power_beta,power_gamma];
end其中δ、θ、α、β和γ频带分别为0.5–3HZ、4–7HZ、8–13HZ、18–30HZ以及>30HZ，具体是根据人在不同状态下的表现来划分的，至于中间衔接的地方空着1HZ我也不清楚为什么，知道的根据需要修改。zcy_0通道1的展现效果如下：
可以看出主要能量集中在δ和θ频段，即0~7HZ。
被测zcy各时期、通道和频段的功率特征显示
为了直观感受一个被测的全部特征信息，这里对用药前和用药后3个小时4个时期，每个时期16个通道，每个通道5段节律进行了遍历显示，代码如下：
load data; %导入数据，包含两组被测，zcy和
zcy=struct('zcy0',zcy_0,'zcy1',zcy_1,'zcy2',zcy_2,'zcy3',zcy_3); %将四个时期数据组成结构体
names=fieldnames(zcy); %获得被测zcy结构体键
allPowFeat={zeros(16,5) zeros(16,5) zeros(16,5) zeros(16,5) }; %若使用三维数组绘图时需要进一步处理
% 此处若全部绘制每个时期每个通道的相关图片，工作量很大，以下仅供手动部分信号比较观察
% for channel=1:2 %手动选择通道
% figure(channel);
% sinPicEEGpro(zcy_0,channel, 10000,11000); %手动选择输入信号和开始结束端点
% end
% 求取用药前后每个时期16个通道5个节律的功率并保存
for i=1:length(names) %遍历用药前后四个时期
key=names(i); %获得当前键
key=key{1};
for channel=1:16 %遍历十六个通道
[powerFeatures] = sinPowEEGpro(zcy.(key),channel,0,length(zcy.(key))); %获取某通道全部节律平均功率
allPowFeat{i}(channel,:)=powerFeatures; %被测zcy全部节律平均功率保存
end
end
% 同一时期16个通道5个节律柱状图
for times=1:4
subplot(2,2,times);bar3(allPowFeat{times});
switch times %设置当前图标题
case 1
title('用药前');
case 2
title('用药后一小时');
case 3
title('用药后两小时');
case 4
title('用药后三小时');
end
set(gca,'xticklabel',{'delta','theta','alpha','beta','gamma'}); %设置各节律名称
end其中第二段被注释掉的部分主要是用于观察某段特定采样点的以上提到的相关波形绘制，sinPicEEGpro函数是对前面提到的几个函数的实例调用，处理某个通道某段数据，其他地方后面也都做了详细注释，没什么需要特别讲的，运行效果如下：
被遮挡部分可以用三维旋转来观察，可以发现随着用药时间推移低频段δ频带有显著增加；用药后的一小时内α频带提升迅速，但在过后的二三小时有基本恢复；β频带和γ频带保持同增同降，在用药一小时内将至最低，在接下来的第二小时出现高幅增长，第三小时开始有恢复用药前趋势；θ频带则一直保持平稳没有太大变化。
根据不同频带的生理特征，可以推测，用药前被测处于较紧张状态，用药后一小时处于情绪较平稳状态，第二个小时处于兴奋状态，第三个小时高频段下降，低频段上升明显，推测处于熟睡状态。这里纯属个人推测，因为我也不知道这个用药用的是什么药，被测是否正常成年人，不同频带的生理特征也只是根据书面的理解。
最后，有两个问题，一、最好还是把时间段细分看看，比如五分钟十分钟，我这里采用的是完整信号，除了用药前8分钟，后面三段信号都约有一个小时；二、不同通道的活跃度即幅值也差别也很大，这是因为不同通道对应与人脑的不同区域，所以至少要分区域讨论，生理这块我不了解，其实我也蛮想知道多通道后续一般怎么处理，比如说分区或是加权等等，没有仔细去找到相关文献，有兴趣的欢迎留言或者私信我。
或扫码搜索关注微信公众号“燕山徐一”，后台回复“脑电预处理”获取代码。
