写在前面 ：

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

疲劳分级

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

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

检测指标

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

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

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

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

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

​
I E M G = ∫ t 1 t 2 ∣ X ( t ) ∣ d t I E M G=\int_{t_{1}}^{t_{2}}|X(t)| d t IEMG=∫t1​t2​​∣X(t)∣dt

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

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

M P F = ∫ 0 ∞ f ⋅ P S D ( f ) d f ∫ 0 ∞ P S D ( f ) d f M 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​
​
∫ 0 M F P S D ( f ) d f = ∫ M F f 0 P S D ( f ) d f = 1 2 ∫ 0 f 0 P S D ( f ) d f \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 ∫0MF​PSD(f)df=∫MFf0​​PSD(f)df=21​∫0f0​​PSD(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。其计算公式为：
M E A N = ∑ i = 1 N R R i N M E A N=\sum_{i=1}^{N} \frac{R R_{i}}{N} MEAN=i=1∑N​NRRi​​

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

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

当人体处于疲劳状态时，SpO2含量有所下降[9]。通过对处理后的实验数据进行时域分析，得到SpO2的特征参数：血氧标准差(RSD)。根据所得结果的波动程度，便可准确的反映出被试者的疲劳状态。其计算公式为：
R S D = 1 N − 1 ∑ i = 1 N ( R i − R A V G ) 2 R S D=\sqrt{\frac{1}{N-1} \sum_{i=1}^{N}\left(R_{i}-R A V G\right)^{2}} RSD=N−11​i=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~=1
    x=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:N
    kk=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)));
end

if i~=1
    xl=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:l
    namelist(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*df
         pp=pp+(px(N)+px(N+1))/2*df;
         N=N+1;
    end
    n_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:l
    namelist(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 实验设备简介