人体疲劳程度检测，生理信号处理

写在前面：

参加比赛做了检测人体人体疲劳程度的手环，题主现在已经毕业。检测的生理指标有肌电、脉搏、血氧含量。对时域和频域进行分析，通过支持向量机分级模型得出疲劳等级。

疲劳分级

本系统根据斯坦福嗜睡量表分为四个疲劳等级：清醒(0级)、轻度疲劳(1级)、中度疲劳(2级)、重度疲劳(3级)

表1 斯坦福嗜睡量表(SSS)

评分	疲劳程度	疲劳等级
1	感觉非常振奋、充满生机与活力	0级
2	身体机能处于较高水平，不过不在峰值，可以集中注意力	0级
3	很清醒，不过身体及思维都比较放松，能及时反应但不够灵敏	1级
4	有些倦意、松懈	1级
5	充满倦意，不再想保持清醒，非常松懈	2级
6	开始打瞌睡，头晕眼花，不再与睡意作斗争，只想躺下来休息	3级
7	睡眠初期，开始出现梦境	3级
X	沉睡	3级

检测指标

根据相关研究，疲劳状态通过选取2个或2个以上检测指标均认为有效，且现有研究中使用频率最高的检测指标主要为肌电、脉搏和血氧饱和度信号，考虑到数据采集的便捷性以及疲劳程度的相关性，本文选取了这3个生理指标作为驾驶疲劳的检测标准。

（1）肌电(EMG)特征参数提取

研究发现EMG信号随着人体疲劳程度的加深，其幅值会增大，频率会降低，因此

通过分析EMG信号来检测驾驶员的疲劳状态。通过对处理后的数据进行时域分析和频域分析，得到肌电信号的特征参数：积分肌电值(IEMG)、均方根值(RMS)、平均功率频率(MPF)、中位频率(MF)。

IEMG是指在一定时间内肌肉中参与活动的运动单位放电总量，通常其幅值随疲劳程度的加深而增加。RMS反映了运动单位数量的募集程度及其放电频率的同步化程度。其计算公式分别为：

IEMG=∫t1t2∣X(t)∣dtI E M G=\int_{t_{1}}^{t_{2}}|X(t)| d t IEMG=∫t1t2∣X(t)∣dt

RMS=1T∑t=1Tx2(t)R M S=\sqrt{\frac{1}{T} \sum_{t=1}^{T} x^{2}(t)} RMS=T1t=1∑Tx2(t)
其中，t1为采集肌电信号开始时刻，t2为采集肌电信号结束时刻，*X(t)*为肌电信号大小。

将表面肌电(EMG)信号进行快速傅里叶变换( FFT)，获得EMG信号的功率谱，反映 EMG信号在不同频率范围内的强度。MF和MPF为常用的肌肉疲劳的评价指标。其计算公式为：

MPF=∫0∞f⋅PSD(f)df∫0∞PSD(f)dfM P F=\frac{\int_{0}^{\infty} f \cdot P S D(f) d f}{\int_{0}^{\infty} P S D(f) d f} MPF=∫0∞PSD(f)df∫0∞f⋅PSD(f)df

∫0MFPSD(f)df=∫MFf0PSD(f)df=12∫0f0PSD(f)df\int_{0}^{M F} P S D(f) d f=\int_{M F}^{f_{0}} P S D(f) d f=\frac{1}{2} \int_{0}^{f_{0}} P S D(f) d f ∫0MFPSD(f)df=∫MFf0PSD(f)df=21∫0f0PSD(f)df
其中，PSD(f)为功率谱密度，MF为中位频率，f为肌电信号频率，f0为频率上限，即采样频率的一半，等于500 Hz，MPF与MF下降表明被试者出现了肌肉疲劳。

（2）脉搏特征参数提取

通过对处理后的脉搏信号数据进行时域分析和频域分析，得到脉搏信号的特征参数：

RR间期均值(MEAN)、RR间期标准差(SDNN)、低频功率(LF)、高频功率(HF)。在脉搏信号中，首先运用小波变换检测出各主波的波峰，求出两个相邻主波波峰之间的时间差，记作RRi(i=1,2,3……)，则RRi的均值即为RR间期均值、RRi的标准差即为RR间期标准差，并运用快速傅里叶变换(FFT)将时域指标变为频域指标，得到LF和HF。其计算公式为：
MEAN=∑i=1NRRiNM E A N=\sum_{i=1}^{N} \frac{R R_{i}}{N} MEAN=i=1∑NNRRi

SDNN=1N∑i=1N(RRi−MEAN)2S D N N=\sqrt{\frac{1}{N} \sum_{i=1}^{N}\left(R R_{i}-M E A N\right)^{2}} SDNN=N1i=1∑N(RRi−MEAN)2
其中，N表示正常脉搏搏动的总次数，RRi表示第i个相邻主波波峰之间的时间差，MEAN表示RR间期均值，SDNN表示RR间期标准差。

（3）血氧饱和度(SpO2)特征参数提取

当人体处于疲劳状态时，SpO2含量有所下降[9]。通过对处理后的实验数据进行时域分析，得到SpO2的特征参数：血氧标准差(RSD)。根据所得结果的波动程度，便可准确的反映出被试者的疲劳状态。其计算公式为：
RSD=1N−1∑i=1N(Ri−RAVG)2R S D=\sqrt{\frac{1}{N-1} \sum_{i=1}^{N}\left(R_{i}-R A V G\right)^{2}} RSD=N−11i=1∑N(Ri−RAVG)2
其中，RAVG表示血氧饱和度平均值，N表示采集样本中血氧数据总个数，Ri表示第i个样本的血氧饱和度值。

程序运行-matlab

程序采用MATLAB运行，运行时候，需 更改文件路径，滤波器用于滤除采集的生理信号的噪声。脉搏、肌电信号在采集过程中均为低频微弱信号，在实验中常会混入各种噪声，包括基线漂移，50Hz工频干扰以及周围环境的噪声干扰等，故采用巴特沃斯滤波器处理信号中的噪声。血氧信号在采集过程中，由于其频率、幅度变化比较稳定，故一般不需对其进行数据预处理

原数据文件：百度网盘提取码：mmb7

巴特沃斯滤波器

%带通滤波
%x-输入信号
%fs-输入信号采样频率
%M-滤波器半阶数
%fl-下截止频率
%fh-上截止频率
%xl-返回的滤波结果,其长度为原序列长度减少2M,去掉了前M和后M个不正确的点
function xl=beltfilter(x,fs,M,fl,fh)
N=length(x)-2*M;
[i j]=size(x);
if i~=1x=x';
end
k=1:M;
w=0.5+0.5*cos(pi*k/M);wl=2*pi*fl/fs;
wh=2*pi*fh/fs;
h(1)=(wh-wl)/pi;
h(2:M+1)=(sin(wh*k)-sin(wl*k))./(pi*k).*w;for k=1:Nkk=k-1+M;xl(k)=x(kk+1)*h(1)+sum(h(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
endif i~=1xl=xl';
end

积分肌电值(IEMG)、均方根值(RMS)、平均功率频率(MPF)、中位频率(MF)

clc;
fs=1000;
%%%==================批量读入文件=======================
path = 'C:\Users\hp\Desktop\EMG\数据\';
namelist = dir([path,'*.txt']);
l = length(namelist);
P = cell(1,1);%定义一个细胞数组，用于存放所有txt文件for i = 1:lnamelist(i).name;%获得该路径下的文件名，如1.txt是相对路径filename{i} = [path,namelist(i).name];%通过字符串拼接获得的就是绝对路径了s = load(filename{i});x=s(:,1);%%=================批量处理数据======================x=x/1000;figure(1);subplot(311),plot(x);title("原信号波形");y=fft(x);    %对信号进行快速Fourier变换mag=abs(y);     %求得Fourier变换后的振幅subplot(3,1,2) ,plot(mag);axis([0 500 0 1.1*max(mag)]);title("傅里叶变换");grid on;x=beltfilter(x,fs,80,0,250);L=length(x);f=(0:L-1)*fs/L;y=fft(x);    %对信号进行快速Fourier变换mag=abs(y);     %求得Fourier变换后的振幅subplot(3,1,3) ,plot(f,mag);axis([0 500 0 1.1*max(mag)]);title("带通滤波器");grid on;iemg=sum(abs(x))/length(x);  %求积分肌电值rms=sqrt(sum(x.^2)/length(x)); %求均方根值L=length(x);f=(0:L-1)*fs/L;cx=xcorr(x,'unbiased');cxk=fft(cx,L);px=abs(cxk);%求功率谱密度pxx=10*log10(px);f=(0:L-1)*fs/L;grid on  %做功率谱图df=fs/L;p=(sum(px(1:L/2-1))+sum(px(1:L/2)))/2.*df;pf=(sum(px(1:L/2-1).*[1:L/2-1]'.*df)+sum(px(1:L/2).*[1:L/2]'.*df))/2*df;MPF=pf/p; %求平均功率频率N=1;pp=0;while abs(pp-p/2)>(px(N)+px(N+1))/2*dfpp=pp+(px(N)+px(N+1))/2*df;N=N+1;endn_1=(N+N+1)/2;MF=df*n_1; %求中值频率fprintf("\n积分肌电值      %.4f\n",iemg);fprintf("均分根值        %.4f\n",rms);fprintf("平均功率频率    %.4f\n",MPF);fprintf("中值频率        %.4f\n\n",MF);%%=================记录数据=====================result1(i)=iemg;result2(i)=rms;result3(i)=MPF;result4(i)=MF;
end%%%==================数据文件的写入=======================
fid = fopen('C:\Users\hp\Desktop\EMG\数据\结果\iemg.txt','wt');
fprintf(fid,'%g\n',result1);
fclose(fid);fid = fopen('C:\Users\hp\Desktop\EMG\数据\结果\rms.txt','wt');
fprintf(fid,'%g\n',result2);
fclose(fid);fid = fopen('C:\Users\hp\Desktop\EMG\数据\结果\MPF.txt','wt');
fprintf(fid,'%g\n',result3);
fclose(fid);fid = fopen('C:\Users\hp\Desktop\EMG\数据\结果\MF.txt','wt');
fprintf(fid,'%g\n',result4);
fclose(fid);%%======================绘图===========================
figure(2)
n=1:l;
subplot(411)
plot(n,result1)
axis([0 1.1*l 0 1.1*max(result1)]);
%suptitle('肌电的4种指标')
title("肌电积分值(IEMG)");
xlabel("时间(min)");subplot(412)
plot(n,result2)
axis([0 1.1*l 0 1.1*max(result2)]);
title("均方根值(RMS)");
xlabel("时间(min)");subplot(413)
plot(n,result3)
axis([0 1.1*l 0 1.1*max(result3)]);
title("平均功率频度(MPF)");
xlabel("时间(min)");subplot(414)
plot(n,result4)
axis([0 1.1*l 0 1.1*max(result4)]);
title("中位频率(MF)");
xlabel("时间(min)");

血氧均值、血氧均方差（RSD）

clc;
%%%==================批量读入文件=======================
path = 'C:\Users\hp\Desktop\SPO2\数据\';
namelist = dir([path,'*.txt']);
l = length(namelist);
P = cell(1,1);%定义一个细胞数组，用于存放所有txt文件for i = 1:lnamelist(i).name;%获得该路径下的文件名，如1.txt是相对路径filename{i} = [path,namelist(i).name];%通过字符串拼接获得的就是绝对路径了s = load(filename{i});x=s(:,1);figure(1)plot(x);title("SPO2数据");axis([0 100 0 1.1*max(x)]);%%=================批量处理数据======================L=length(x);RAVG=sum(x)/L; %%血氧均值RSD=sqrt(  sum((x-RAVG).^2)/L  );  %%血氧均方差fprintf("\n  血氧均方差:  %.4f\n\n",RSD);result(i)=RSD;RAVG1(i)=RAVG;
end%%%==================数据文件的写入=======================
fid = fopen('C:\Users\hp\Desktop\SPO2\数据\结果\result.txt','wt');
fprintf(fid,'%g\n',result);
fclose(fid);%%======================绘图===========================
figure(2)
n=1:l;
subplot(211)
plot(n,RAVG1)
title("血氧平均值");
grid on;% set(gca,'xminortick','on');  %坐标刻度显示
% set(gca,'ticklength',[0.02 0.025]);
% set(gca,'tickdir','out');xlabel("时间(min)");
ylabel("含量(%)");subplot(212)
plot(n,result)
title("血氧均方差")
axis([0 1.1*l 0 1.1*max(result)]);
xlabel("时间(min)");
grid on;

实验设备

表2 实验设备简介

设备名称	数量	功能简介
SEN0240肌电传感器	1个	采集驾驶员的肌电信号
MAX30102脉搏血氧传感器	1个	采集驾驶员的脉搏信号和血氧信号
STM32F103单片机	2个	设备的控制和数据处理单元
HC-05蓝牙模块	2个	实现手环系统和车外警示系统的连接与通信
RGB七彩LED灯	1个	依据不同颜色，显示出驾驶员当前所处的疲劳状态
车外警示灯	1个	给予道路上其他驾驶员警示作用
振动模块	1个	通过振动，提醒驾驶员应注意休息
语音模块	1个	通过语音提示，提醒驾驶员应注意休息