表面肌电信号试验手册总结

1、实验 基于sEMG时域特征特的动作识别一、实验目的, 了解肌电信号常用的时域分析方法；.利用MATLAB对肌电信号进行去噪、特征提取及动作识别；二、实验设备Wi-Fi表面肌电信号采集卡；32 位 Windows XP 台式机（Matlab 7.0 软件）；802.11b/g无线网卡；三、实验内容（1）学习信号的基本去噪方法，并用 MATLAB实现；（2）学习肌电信号常用的时域特征并利用 Matlab来进行波形长度（WL） 符号改变数（SSC）、过零点（ZC）、威尔逊赋值（WAMP）等特征的提取；（3）学习神经网络信号处理方法，掌握 BP神经网络的用法，将其用于肌 电信号的动作识别。学习以上三个

2、部分，最终完成一整套肌电信号去噪、 特征提取（选取一种特 征）、基于特征的动作识别的 MATLAB程序。四、实验原理小波去噪小波去噪方法是一种建立在小波变换基础上的新兴算法，基本思想是根据噪 声在不同频带上的小波分解系数具有不同强度分布的特点， 将各频带上的噪声对 应的小系数去除，保留原始信号的小波分解系数，然后对处理后系数进行小波重 构，得到纯净信号。小波去噪的基本原理图如下1生理信号采集与分析处理平台实验指导-武汉理工大学信息工程学院(2)特征提取时域分析是将肌电信号看成均值为零，而方差随着信号强度的变化而变化的 随机信号。时域特征的计算复杂度低，提取比较方便。最常用的方法有：方差，过零点

3、数(Zero Crossing, ZC),Willison幅值(Willison Amplitude, WAMP),绝对值平均值 (Mean Absolute Value, MAV)和波形长度 (Wave length, WL)等。在实际应用中，为了让特征可以包含更多的信息，往 往选择用不同的时域特征组合形成联合特征向量。我们主要介绍一下几种方法：过零率(ZC):为波形通过零线的次数，从一定程度上反映了信号的频率特 性。为了降低零点引入的噪声，往往会引入一个阈值6。计算方式如下：sgn(-Xk + 人+)(“ 一 乂女书之讨(1)Willison幅值：是由Willison提出一种对表面肌电信号

4、的幅值变化数量进行 计算的方法，经过后人的研究，对 Willison幅值的阈值有了明确的范围限定，目 前认为50100叩 是最合适的阈值范围。其数学表示公式如公式(3-3)。NWAMP = f x -xi +i其中:1 f (x) = .0if x阈值 otherwise波形长度(WL ):它是对某一分析窗中的波形长度的统计，波长可以体现该样本的持续时间、幅值、频率的特征。1 N 1WL=? x(i 1)-x(i)(3)N 1符号改变斜率(SSC：为信号的的频率性能提供了一些附加信息，对于 3 个连续的采样点，给定阈值 ，通过下面的公式计算波峰波谷的个数。(4)为-为二 X -X 1 一，i

5、=1, ,N(3)神经网络BP神经网络又称误差反向传播(Back Propagation),它是一种多层的前向 型神经网络。在BP网络中，信号是前向传播的，而误差是反向传播的。所谓的 反向传播是指误差的调整过程是从最后的输出层依次向之前各层逐渐进行的。标准的BP网络采用梯度下降算法，与 Widrow-Hoff学习规则相似，网络权值沿着 性能函数的梯度反向调整。生理信号采集与分析处理平台实验指导-武汉理工大学信息工程学院前向型神经网络通常具有一个或多个由 sigmoid神经元构成的隐层，以及一 个由线性神经元构成的输出层。多个具有非线性传递函数的神经元层使得网络可以学习输入和输出之间的非线性关系

6、，而线性输出层使得网络可以产生-1, +1之外的输出值。输入隐层a1 =tan sig(IW 1, p -b1)输出层a2 = purelin (IW 2,1a1，b2)由两层神经元构成的BP网络结构BP网络的训练算法BP算法BP算法沿着误差函数减小最快的方向，也就是梯度的反方向改变权值和偏 差，这一点与线性网络的学习算法是一致的。BP算法的迭代计算公式可以表示为：Xj = Xk-akgk(1)其中，Xk代表当前权值和偏差,Xk书代表迭代产生的下一次的权值与偏差, gk为当前误差函数的梯度，ak代表学习速率。 有动量的梯度下降算法标准的梯度下降法在调整权值时，仅仅按照当前时刻的负梯度方向进行调

7、(2)整，并没有考虑以前各次运算步骤中的梯度方向， 因此新的样本对迭代过程影响 太大，可能会导致训练过程中调整方向发生震荡， 导致不稳定和收敛速度慢的问 题，有动量的梯度下降算法则考虑了往前时刻的贡献，其权值迭代算法为：wj (n 1) = wj (n)(1 - )D(n)= D(n -1)其中，D(n), D(n -1)分别表示n时刻，n-1时刻的负梯度。由于加入了以前时刻梯度的贡献，相当于给迭代过程添加了一个低通滤波器， 使得网络忽略误差生理信号采集与分析处理平台实验指导-武汉理工大学信息工程学院曲面上细节特征，避免了陷入局部极小点的问题共腕梯度算法尽管标准的BP算法采用梯度下降算法，权值

8、和偏差沿误差函数下降最快的 方向调整，但却并不一定是收敛最快的算法。在改进的 BP训练算法中，有一大 类的算法称为共腕梯度算法。在这一类算法中，权值和偏差沿着共腕梯度方向进 行调整，通常能够获得比标准的梯度算法更快的收敛速度。共腕梯度算法的第一次迭代都是从最陡下降的梯度方向开始。梯度向量为：p0 - -g0沿着此方向进行权值和偏差的调整，公式为：Xk .1 =Xk akgk下一次搜索方向则由前两次搜索方向的共腕方向决定，表达式为：Pk = -gk ，Pk对于系数久不同计算方法产生不同的共腕梯度算法。F-R共腕梯度算法采取的系数确定方法为：Tgkgk-k gkg一即本次迭代梯度相对于上一次迭代梯

9、度的归一化值。P-R共腕梯度算法采取的系数确定方法为：b _应-k gk9即上次迭代梯度与本次迭代梯度的内积对本次梯度的归一化值。Scaled共腕梯度算法到目前为止，讨论过的所有共腕梯度算法都需要在每一步迭代过程中对搜索 方向进行计算，这样的计算量是比较大。对此 moller提出了 Scaled梯度搜索算 法4,在每一步迭代过程中不计算搜索方向，以减少训练过程的计算量。其基本 原理是利用下面介绍的L-M算法与共腕梯度法相结合产生的。L-M算法L-M算法其权值和阈值的更新过程为：xg =Xk -JTJ +NI，Te(8)生理信号采集与分析处理平台实验指导-武汉理工大学信息工程学院其中，e为期望输

10、出与实际输出的误差；J为误差对权值微分的Jacobi矩阵； N为标量因子。如果训练成功，误差性能函数减小，那么就减小 N的值；反之就 减小其值。五、实验步骤1,认真阅读本实验的原理与方法。2.利用MATLAB实现对肌电信号的去噪、特征提取及动作识别。(1)小波去噪load sample_test.mat;load sample_train.mat;%参数说明 % %Window=256;%分析窗口的长度M=512;%采集数据时一个data的样本数Channel=4; %采集数据的通道数Class=6; %类别数Number=10; %每个类别的个数WinLap=64; %窗口移动的间隔Judg

11、eTime=Window/WinLap; %一个分析窗口需要移动的次数8Count=M*Number*Class/WinLap-Window/WinLap+1; % 所有数据需要分析的 次数 477ClassCount=M*Number/WinLap-Window/WinLap+1;%一类数据需要的分析次数 77GapCount=M*Number/WinLap;% 训练样本两类动作之间的间隔80ClassOne=1;ClassTwo=2;ClassThree=3;ClassFour=4;生理信号采集与分析处理平台实验指导-武汉理工大学信息工程学院ClassFive=5;ClassSix=6;%

12、thr=0.2;% train 样 本 小 波 去 噪 % %for i=1:Channelx=sample_train(:,i);thr,sorh,keepapp=ddencmp(den,wv,x);c,l=wavedec(x,3,db1);a3=appcoef(c,l,db1,3);d3=detcoef(c,l,3);d2=detcoef(c,l,2);d1=detcoef(c,l,1);s4=wdencmp(gbl,c,l,db1,3,thr,sorh,keepapp);new_train(:,i)=s4;endfigure(1)subplot(4,1,1);plot(sample_tr

13、ain(:,1),r);hold on;plot(new_train(:,1),b);legend( 原始信号,去噪后信号);title(train样本滤波前后信号对比);subplot(4,1,2);plot(sample_train(:,2),r);hold on;plot(new_train(:,2),b);subplot(4,1,3);plot(sample_train(:,3),r);hold on;plot(new_train(:,3),b);subplot(4,1,4);plot(sample_train(:,4),r);hold on;plot(new_train(:,4),b

14、);train样又去噪前后的信号对比如下图而入生理信号采集与分析处理平台实验指导-武汉理工大学信息工程学院train样本滤波前后信号对比 TOC o 1-5 h z 1，卜s一二原始信号j00.5111.52去噪后信号I.521pII01i1a1BQ0.511.52Z533.5X 10*仿照对train样本进行小波去噪的方法，写出利用MATLAB对test样本进行 小波去噪的程序。(2)特征提取clear all;close all;load sample_train.mat;%参数说明%Window=256;%分析窗口的长度M=512;%采集数据时一个data的样本数Channel=4; %

15、采集数据的通道数Class=6; %类别数Number=10; %每个类别的个数WinLap=64; %窗口移动的间隔JudgeTime=Window/WinLap; %一个分析窗口需要移动的次数Count=M*Number*Class/WinLap-Window/WinLap+1;% 所有数据需要分析的次数CountClass=M*Number/WinLap-Window/WinLap+1; % 一类数据需要的分析次数CountGap=M*Number/WinLap;%训练样本两类动作之间的间隔生理信号采集与分析处理平台实验指导-武汉理工大学信息工程学院% 求波长 %sample=sampl

16、e_train; %H 口 WL=zeros(Channel,Count);for c=1:Countfor n=1:Channelfor w=1:Window-1WL(n,c)=WL(n,c)+abs(sample(c-1)*WinLap+w,n)-sample(c-1)*WinLap+w+1,n); end endendFeature_train=WL;仿照求WL (波长)特征的程序，任选其他一种特征，写出 MATLAB程序。 注意，train样本和test样本都需要做特征提取。将最终求得的 train和test样本 的特征存入feature.mat(3)动作识别clear all;clo

17、se all; load feature.mat Gap=80;K=77;Count=477;%共有477组特征值P=6;%6类动作Feature=Feature_test;噬一接 口test=Feature_train;%E一接 口trainFunnet=newff(minmax(Feature),5,6,tansig,purelin,trainlm);%net = newff ( A, B, C,%参数：% A: 一个nx 2的矩阵，第i行元素为输入信号xi的最小值和最大值;% B: 一个k维行向量，其元素为网络中各层节点数；生理信号采集与分析处理平台实验指导-武汉理工大学信息工程学院%

18、C: 一个k维字符串行向量，每一分量为对应层神经元的激活函数； 常用的有： % 线性函数purelin,对数S形转移函数logsig双曲正切S形函数tansig % trainFun :为学习规则采用的训练算法。也可以根据需要修改%设置训练的输出目标矩阵ze1=zeros(1,Gap)+1;ze0=zeros(1,Gap);ze1_6=zeros(1,K)+1;ze0_6=zeros(1,K);t=ze1 ze0 ze0 ze0 ze0 ze0_6ze0 ze1 ze0 ze0 ze0 ze0_6ze0 ze0 ze1 ze0 ze0 ze0_6ze0 ze0 ze0 ze1 ze0 ze0_

19、6ze0 ze0 ze0 ze0 ze1 ze0_6ze0 ze0 ze0 ze0 ze0 ze1_6;%神经网络进行训练tic ;net.trainParam.show=200;net.trainParam.epochs=10000;net.trainParam.goal=0.01;net.trainParam.lr=0.01;net,tr=train(net,Feature,t);Y=sim(net,Feature);YS=sim(net,test);toc;% 一些重要的网络配置参数如下：%net.trainparam.goal :神经网络训练的目标误差%net.trainparam.s

20、how: 显示中间结果的周期%net.trainparam.epochs:最大迭代次数%net.trainParam.lr :学习率生理信号采集与分析处理平台实验指导-武汉理工大学信息工程学院% 语法：net, tr = train( net, X, Y )%参数：%X：网络实际输入%Y:网络应有输出%tr：训练跟踪信息%Y1:网络实际输出% 语法：Y=sim(net,X)N为数据样本%net:网络%X:输入给网络的KX N矩阵，其中K为网络输入个数,数% Y:输出矩阵QXN,其中Q为网络输出个数%判断输出结果为第几类maxY=max(Y,1);for i=1:Countfor j=1:Pif(maxY(i)=Y(j,i) class(i)=j;end endendmaxYS=max(YS,1);for i=1:Countfor j=1:Pif(maxYS(i)=YS(j,i) class_test(i)=j;end endendT=Gap;%80wrong=zeros(1,Count);for k=1:Gapif (class_test(k)=1)wrong(k)=wrong(k)+class_te