实验四肌电信号的肌肉疲劳估计.docx

实验四肌电信号的肌肉疲劳估计肌肉在持续的收缩过程中，会逐渐进入疲劳状态，肌肉疲劳特性的研究在康复医学、运动医学领域具有重要作用肌电信号(electromyogram,EMG)是从人体骨骼肌表面通过电极记录下来的神经肌肉活动时发放的生物电信号，它反映了神经、肌肉的功能状态，因此通过EMG研究肌肉疲劳是一个有效途径已有许多研究发现，在疲劳过程中EMG信号会出现幅度增长，功率谱朝低频方向移动等现象这些效应是由于神经传导速率的变化所引起的，会对肌电假肢的控制以及运动力量的估测等造成不利影响因此对肌肉疲劳的检测以及疲劳程度的度量显得非常必要本文通过实验采集到实验者的肌电信号，对其进行了预处理，并且定量分析估计了肌肉疲劳的过程1.EMG的采集本实验采用生理信号采集仪MP150采集肌电信号的表面电极使用一次性电极，型号为LT-301,材料为Ag/AgC"采样频率2KHz,放大倍数500倍实验者均采用坐姿，在上臂的肱三头肌，肘肌，肱二头肌，肱楼肌(肱二头肌、肱横肌分别是屈肘动作的主动肌和协同肌；肱三头肌、肘肌分别是伸肘动作的主动肌和协同肌)上分别贴上表面电极肘部动作的起始位置设置在裤缝线处，手臂尽量与水平面垂直。

动作的终止位置大约在水平位置先屈肘后伸肘实验者不间断均匀重复举重为5.5kg的哑铃，感到疲劳时记录疲劳前举重次数，然后重复举重直到肌肉无力举起为止EMG是一种非常复杂的信号，信号本身非常微弱，稳定性较差，随机性很强因此信号检测时需要注意以下相关事项：①电极位置：电极所在位置应受其他肌肉串扰的影响最小检测电极应置于肌腹的中间，尽量离其他肌肉足够远；电极对的方向应与肌梭方向平行参考电极尽可能置于肌肉最少的地方②检测电极对的距离：检测电极间隔的距离越大，拾取的信号越广越深，信号的幅值也越大,因此为了保证测量的可比性，每次测量时电极间隔的距离应固定③皮肤阻抗：人体皮肤阻抗高达10、100kQ/cm2，变化范围很大皮肤-电极界面的阻抗越大,sEMG的基线就越不稳定，运动时测得的sEMG就越不精确因此测量时应预先对皮肤进行处理，如去死皮和油脂④电磁干扰：注意日光灯、电动机等电子设备带来的电磁干扰⑤受试者：训练受试者的运动技巧以及调节能力，减少其他肌肉、具有较好的实时性1. 积分肌电值(iEMG)和均方根幅值(RMS)积分肌电值(iEMG)是指所得肌电信号经整流滤波求单位时间内曲线下面积的总和，是评价疲劳的重要手段。

均方根振幅(RMS)用来描述一段时间内肌电的平均变化特征，指此段时间内所有振幅的均方根值，但不能反映肌电信号的细节变化疲劳时肌电信号的振幅增高，势必引起RMS的增加可以通过比较不同时期的RMS,确定疲劳发生的时间和疲劳的程度大多数研究结果表明从初始态到疲劳态，EMG的时域值在运动至疲劳过程中的总趋势是上升，反映了参加工作的运动单位的数量其计算公式如下：t+TiEMG=(3.1.1)t(3.1.2)2. 过零点数(ZC)肌电信号源于中枢神经所发送的电脉冲，肌电信号的强度还与电脉冲的频率有关，所以，过零点数有时也作为肌电信号的一个特征，该特征值定义如下：N-1(3.1.3)ZC=Zsgn(-x(i)x(，+l)i-o3. 其中，sgn(x)=

其计算方法如下式子所示：JRMMPF=ocJp(砒0MF®]oojp(n#=Jp(n#=Hp(n#0MF20(3.2.1)(3.2.2)其中P（f）为肌电功率谱3. 基于matlab的肌电信号处理结果及分析（matlab程序见附录）1. ，matlab编程读取肌电信号，如下图3.1所示：2. 利用窗函数法设计的FIR带通滤波器对原肌电信号进行处理，得到图3.2,如下所示：（图3.2滤波后的肌电信号）为了清楚的对比疲劳前后信号的幅度变化，将其按照疲劳时间分为疲劳前后两部分，得到图3.3所示幅值变换情况由图3.3可以看到，疲劳前，幅度变化明显小于疲劳之后图3.3疲劳前后肌电信号对比）2004006008001000频率/Hz疲劳前平均功率谱图4O3O2O1OD-10-2040-2002004006008001000频率/Hz3U疲劳后平均功率谱图ooO21COP瑜格林（图3.4平均功率谱图）matlab编程计算每个实验时域和频域四个参数，可以得到下表(B表示疲劳前,A表示疲劳后)姓名参数BABABABABAIEMGRMSMPFMFZCZJH0.34940.67470.51831.0653193.2039210.765080.421586.4802775305WL0.10300.72720.62151.1192182.8710202.412579.009981.4253304365MDR0.34980.66160.54311.0259199.5212206.341886.921184.1783298315LZD0.20900.37700.29320.6199205.7864197.810395.583889.7699342329HZS0.41970.32840.61620.6622190.9405227.382973.0457100.0210249231从表中我们可以看出，HZS的积分肌电值是减小的，这是不符合理论的；平均功率频率只有LZD的是减小的，其他实验者的按实际都应该是减小，但这里是增加的；中位频率也只有MDR和LZD两组符合实际原理。

究其原因，一方面，实验条件的有限性；另一方面实验者本身疲劳时间的反应可能已经疲劳或者还没疲劳，导致实验计算时疲劳前后分段有误3. 参考文献[1] 王奎，刘建红，宋刚.sEMG技术在评价运动性疲劳方而的方法及应用[J].安徽体育科技,2004,23(3)：49-51王健，刘如海.肌肉疲劳的表面肌电信号特征研究与展望[J].中国体育科技，2003,39(2)：4-7附录clear;clc;七读入数据fs=2000;fidin=fopen('*);%用£open命令打开位于工作区间的文件fidout=fopen('vr.txt',，沪);%fopen命令打开在工作区I'司创建的文件vr.zxt用于存放将来提取的信号数据while~feof(fidin)%判断是否为文件末尾tline=fgetl(fidin);%从文件读行数据fprintf(fidout,'%s\n\n',tline);%把数据写入上面所创建的文件中vr.txt,并且换行为下次写入做准备endfclose(fidout);out=importdata(1vr.txtr;figure(1);plot(out);title(1原肌电信v1);xlabel(1时in]/s1);ylabel(1度/V);霍窗函数法FIR滤波器处理N=201;冬窗长window=boxcar(N);舍矩形窗b=firl(200,0.01,'high1,window);sf=filter(b,lzout);b=firl(200,0.5rwindow);sf=filter(b,lzsf);figure(2);plot(sf);title(1滤波后信号】);xlabel(1时间/s1);ylabel('幅度/V1);nl=ceil(28/50*length(sf))-1;xl=sf(1:nl,1);x2=sf(nl+1:length(sf),1);figure(3);subplot(121)/plot(xl);title('疲劳前1);xlabel(1时间/sf);ylabel(*幅度/V);axis([0150000-45]);subplot(122);plot(x2);title(，疲劳后，);xlabel('时间/s1);ylabel('幅度/Vf);号积分肌电值iemgl=sum(abs(xl))/length(xl);iemg2=sum(abs(x2))/length(x2);号均方根值rmsl=sqrt(sum(xl.A2)/length(xl));rms2=sqrt(sum(x2.人2)/length(x2));号平均功率频率MPFl和中位频率MFILl=length(xl);cxl=xcorr(xl,'unbiased1);告计算xlG勺自相关函数cxkl=fft(oxi,LI);pxl=abs(cxkl)勇求功率谱密度pxxl=10*logl0(pxl);fl=(0:Ll-l)*fs/Ll;figure(4);subplot(3,2,1)zplot(f1(1:Ll/2),pxxl(1:Ll/2))xlabel(，频率/Hz1);ylabel(1功率谱/dB1);title(-平均功率谱图;gridon曹做功率谱图dfl=fs/Ll;pl=(sum(pxl(1:Ll/2-1))+sum(pxl(1:Ll/2)))/2.*df1;pfl=(sum(pxl(l:Ll/2-l).*[l:Ll/2-l]1.*df1)+sum(pxl(l:Ll/2).*[1:L1/2]1.*dfl))/2*dfl;MPFl=pfl/pl;%求平均功率频率Nl=l;ppl=0;whileabs(ppl-pl/2)>(pxl(Nl)+pxl(Nl+1))/2*df1ppl=ppl+(pxl(Nl)+pxl(Nl+1))/2*df1;N1=N1+1;endn_l=(Nl+Nl+l)/2;MFl=dfl*n_l;%求中值频率L2=length(x2);cx2=xcorr(x2,1unbiased1);cxk2=fft(cx2rL2);px2=abs(cxk2);*求功率谱密度pxx2=10*log!0(px2);f2=(0:L2-l)*fs/L2;subplot(3,2,2)zplot(f2(l:L2/2),pxx2(l:L2/2))xlabel(，频率/Hz1);ylabel(1功率谱/dB1);title(•平均功率谱图gridon号做功率谱图df2=fs/L2;p2=(sum(px2(1:L2/2-1))+sum(px2(1:L2/2)))/2.*df2;pf2=(sum(px2(1:L2/2-l).*[l:L2/2-l]1.*df2)+sum(px2(1:L2/2).*[1:L2/2]1.*df2))/2*df2;MPF2=pf2/p2;%求平均功率频率N2=l;pp2=0;whileabs(pp2-p2/2)>(px2(N2)+px2(N2+1))/2*df2pp2=pp2+(px2(N2)+px2(N2+1))/2*df2;N2=N2+1;endn_2=(N2+N2+l)/2;MF2=df2*n_2;%求中值频率IEMG=[iemgl,iemg2]RMS=[rmslzrms2]MPF=[MPF1,MPF2]。