表面肌实验报告

1、武汉理工大学现代数字信号处理在前沿学科中的应用实验报告基于sEMG时域特征的动作识别学院： 信息工程学院 学号： 79 姓名： 吴志勇 班级： 电子154 实验 基于sEMG时域特征特的动作识别1、 实验目的1. 了解肌电信号常用的时域分析方法；2. 利用MATLAB对肌电信号进行去噪、特征提取及动作识别；2、 实验设备1. Wi-Fi表面肌电信号采集卡；2. 32位Windows XP台式机（Matlab 7.0软件）；3. 802.11b/g无线网卡；3、 实验内容（1）学习信号的基本去噪方法，并用MATLAB实现；（2）学习肌电信号常用的时域特征并利用Matlab来进行波形长度（WL）符

2、号改变数（SSC）、过零点（ZC）、威尔逊赋值（WAMP）等特征的提取；（3） 学习神经网络信号处理方法，掌握BP神经网络的用法，将其用于肌电信号的动作识别。学习以上三个部分，最终完成一整套肌电信号去噪、特征提取（选取一种特征）、基于特征的动作识别的MATLAB程序。4、 实验原理（1） 小波去噪小波去噪方法是一种建立在小波变换基础上的新兴算法，基本思想是根据噪声在不同频带上的小波分解系数具有不同强度分布的特点，将各频带上的噪声对应的小系数去除，保留原始信号的小波分解系数，然后对处理后系数进行小波重构，得到纯净信号。小波去噪的基本原理图如下（2） 特征提取时域分析是将肌电信号看成均值为零，而方

3、差随着信号强度的变化而变化的随机信号。时域特征的计算复杂度低，提取比较方便。最常用的方法有：方差，过零点数（Zero Crossing, ZC），Willison幅值（Willison Amplitude, WAMP），绝对值平均值 (Mean Absolute Value, MAV）和波形长度（Wave length，WL）等。在实际应用中，为了让特征可以包含更多的信息，往往选择用不同的时域特征组合形成联合特征向量。我们主要介绍一下几种方法：过零率（ZC）：为波形通过零线的次数,从一定程度上反映了信号的频率特性。为了降低零点引入的噪声，往往会引入一个阈值。计算方式如下：（1）Willison

4、幅值：是由Willison提出一种对表面肌电信号的幅值变化数量进行计算的方法，经过后人的研究，对Willison幅值的阈值有了明确的范围限定，目前认为 是最合适的阈值范围。其数学表示公式如公式（3-3）。（2）其中：波形长度（WL）：它是对某一分析窗中的波形长度的统计，波长可以体现该样本的持续时间、幅值、频率的特征。（3）符号改变斜率（SSC）：为信号的的频率性能提供了一些附加信息，对于3个连续的采样点，给定阈值，通过下面的公式计算波峰波谷的个数。（4）（3） 神经网络BP神经网络又称误差反向传播（Back Propagation），它是一种多层的前向型神经网络。在BP网络中，信号是前向传播的

5、，而误差是反向传播的。所谓的反向传播是指误差的调整过程是从最后的输出层依次向之前各层逐渐进行的。标准的BP网络采用梯度下降算法，与Widrow-Hoff学习规则相似，网络权值沿着性能函数的梯度反向调整。前向型神经网络通常具有一个或多个由sigmoid神经元构成的隐层，以及一个由线性神经元构成的输出层。多个具有非线性传递函数的神经元层使得网络可以学习输入和输出之间的非线性关系，而线性输出层使得网络可以产生-1，+1之外的输出值。由两层神经元构成的BP网络结构（1） BP网络的训练算法 BP算法BP算法沿着误差函数减小最快的方向，也就是梯度的反方向改变权值和偏差，这一点与线性网络的学习算法是一致的

6、。BP算法的迭代计算公式可以表示为： （1）其中，代表当前权值和偏差，代表迭代产生的下一次的权值与偏差，为当前误差函数的梯度，代表学习速率。 有动量的梯度下降算法标准的梯度下降法在调整权值时，仅仅按照当前时刻的负梯度方向进行调整，并没有考虑以前各次运算步骤中的梯度方向，因此新的样本对迭代过程影响太大，可能会导致训练过程中调整方向发生震荡，导致不稳定和收敛速度慢的问题，有动量的梯度下降算法则考虑了往前时刻的贡献，其权值迭代算法为： （2）其中，分别表示n时刻，n-1时刻的负梯度。由于加入了以前时刻梯度的贡献，相当于给迭代过程添加了一个低通滤波器，使得网络忽略误差曲面上细节特征，避免了陷入局部极小

7、点的问题。 共轭梯度算法尽管标准的BP算法采用梯度下降算法，权值和偏差沿误差函数下降最快的方向调整，但却并不一定是收敛最快的算法。在改进的BP训练算法中，有一大类的算法称为共轭梯度算法。在这一类算法中，权值和偏差沿着共轭梯度方向进行调整，通常能够获得比标准的梯度算法更快的收敛速度。共轭梯度算法的第一次迭代都是从最陡下降的梯度方向开始。梯度向量为： （3）沿着此方向进行权值和偏差的调整，公式为： （4）下一次搜索方向则由前两次搜索方向的共轭方向决定，表达式为： （5）对于系数不同计算方法产生不同的共轭梯度算法。a）F-R共轭梯度算法采取的系数确定方法为： （6）即本次迭代梯度相对于上一次迭代梯度

8、的归一化值。b）P-R共轭梯度算法采取的系数确定方法为： （7）即上次迭代梯度与本次迭代梯度的内积对本次梯度的归一化值。c）Scaled共轭梯度算法到目前为止，讨论过的所有共轭梯度算法都需要在每一步迭代过程中对搜索方向进行计算，这样的计算量是比较大。对此moller提出了Scaled梯度搜索算法4，在每一步迭代过程中不计算搜索方向，以减少训练过程的计算量。其基本原理是利用下面介绍的L-M算法与共轭梯度法相结合产生的。 L-M算法L-M算法其权值和阈值的更新过程为： （8）其中，e为期望输出与实际输出的误差；J为误差对权值微分的Jacobi矩阵；为标量因子。如果训练成功，误差性能函数减小，那么就

9、减小的值；反之就减小其值。5、 实验步骤（1） 数据格式转换：N=4;%通道的个数M=512;%每个通道的采样数K=10;%每种类别所取的个数P=6;%类别数ClassPermode=10;%每类动作数据的个数ClassNumber=6;%类别的个数M=512;%每个样本点数据的个数Channel=4;%通道的个数Count=ClassPermode*ClassNumber*M;sample_train=zeros(Count,Channel);%保存训练样本sample_test=zeros(Count,Channel);%保存训练样本for i=1:60 D=importdata(strc

10、at(train2Data,num2str(i),.txt); sample=D.data; sample_train(i-1)*M+1:i*M),:)=sample(:,(2:5);%sample(:,(2:5); endfor i=1:60 D=importdata(strcat(test2Data,num2str(i),.txt); sample=D.data; sample_test(i-1)*M+1:i*M),:)=sample(:,(2:5); endsave sample_train sample_trainsave sample_test sample_test程序运行后，会看

11、到文件夹中多出了sample_train和sample_test两个.mat文件。（2） 小波去噪及特征提取，选用WAMP特征的提取。load sample_test.mat;load sample_train.mat;% 参数说明 %Window=256;%分析窗口的长度M=512; %采集数据时一个data的样本数Channel=4; %采集数据的通道数Class=6; %类别数Number=10; %每个类别的个数WinLap=64; %窗口移动的间隔JudgeTime=Window/WinLap; %一个分析窗口需要移动的次数 8Count=M*Number*Class/WinLap-

12、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;%thr=0.2;% train样本小波去噪 %for i=1:Channel x=sample_train(:,i); thr,sorh,keepapp=ddencmp(den,wv,x)

13、; 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_train(:,1),r);hold on;plot(new_train(:,1),b);legend(原始信号,去噪后信号);title(train样本滤波前后信号对比);subplot(

14、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);% trin样本特征提取 %Threshold=0.006;WAMP=zeros(Channel,Count);for c=1:Count for n=1:Channel for w=

15、1:Window-1 if(abs(new_train(c-1)*WinLap+w,n)-new_train(c-1)*WinLap+w+1,n)Threshold) WAMP(n,c)=WAMP(n,c)+1; end end endendFeature_train=WAMP;%figure(3)%plot(Feature_train);% test样本小波去噪%for i=1:Channel x=sample_test(:,i); thr,sorh,keepapp=ddencmp(den,wv,x); c,l=wavedec(x,3,db1); a3=appcoef(c,l,db1,3);

16、 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_test(:,i)=s4;endfigure(2)subplot(4,1,1);plot(sample_test(:,1),r);hold on;plot(new_test(:,1),b);legend(原始信号,去噪后信号);title(test样本滤波前后信号对比);subplot(4,1,2);plot(sample_test(:,2),r);hold on;plot(new_t

17、est(:,2),b);subplot(4,1,3);plot(sample_test(:,3),r);hold on;plot(new_test(:,3),b);subplot(4,1,4);plot(sample_test(:,4),r);hold on;plot(new_test(:,4),b);% test样本特征提取 %Threshold=0.006;WAMP=zeros(Channel,Count);for c=1:Count for n=1:Channel for w=1:Window-1 if(abs(new_test(c-1)*WinLap+w,n)-new_test(c-1

18、)*WinLap+w+1,n)Threshold) WAMP(n,c)=WAMP(n,c)+1; end end endendFeature_test=WAMP;save feature Feature_train Feature_test;train样本和test样本滤波前后信号对比。（3） 模式识别：clear all; close all;load feature.matGap=80;K=77;Count=477;%共有477组特征值P=6;%6类动作Feature=Feature_test;%统一接口test=Feature_train;%统一接口net1=newff(minmax(F

19、eature),7,6,tansig,purelin,trainlm);%net = newff ( A, B, C ,trainFun)% % 参数：% % A：一个n2的矩阵，第i行元素为输入信号xi的最小值和最大值；% B：一个k维行向量，其元素为网络中各层节点数；% C：一个k维字符串行向量，每一分量为对应层神经元的激活函数；常用的有：% 线性函数purelin，对数S形转移函数logsig 双曲正切S形函数 tansig% trainFun ：为学习规则采用的训练算法。也可以根据需要修改%设置训练的输出目标矩阵ze1=zeros(1,Gap)+1;ze0=zeros(1,Gap);z

20、e1_6=zeros(1,K)+1;ze0_6=zeros(1,K);t=ze1 ze0 ze0 ze0 ze0 ze0_6 ze0 ze1 ze0 ze0 ze0 ze0_6 ze0 ze0 ze1 ze0 ze0 ze0_6 ze0 ze0 ze0 ze1 ze0 ze0_6 ze0 ze0 ze0 ze0 ze1 ze0_6 ze0 ze0 ze0 ze0 ze0 ze1_6;%神经网络进行训练tic ; net1.trainParam.show=200; net1.trainParam.epochs=10000; net1.trainParam.goal=0.1; net1.train

21、Param.lr=0.01; net1,tr=train(net1,Feature,t); Y=sim(net1,Feature); YS=sim(net1,test); toc; % 一些重要的网络配置参数如下：% net.trainparam.goal ：神经网络训练的目标误差% net.trainparam.show ： 显示中间结果的周期% net.trainparam.epochs ：最大迭代次数% net.trainParam.lr ： 学习率% 语法： net, tr = train( net, X, Y )% 参数：% X：网络实际输入% Y：网络应有输出% tr：训练跟踪信息

22、% Y1：网络实际输出% 语法：Y=sim(net,X)% net：网络% X：输入给网络的N矩阵，其中K为网络输入个数，N为数据样本数% Y：输出矩阵QN，其中Q为网络输出个数%判断输出结果为第几类 maxY=max(Y,1); for i=1:Count for j=1:P if(maxY(i)=Y(j,i) class(i)=j;end end end maxYS=max(YS,1); for i=1:Count for j=1:P if(maxYS(i)=YS(j,i) class_test(i)=j;end end endT=Gap;%80wrong=zeros(1,Count);for k=1:Gap if (class_test(k)=1) wrong(k)=wrong(k)+class_test(k);end if (class_test(k+T)=2) wrong(k+T)=wrong(k+T)+class_test(k+T);end if (class_test(k+2*T)=3) wrong(k+2*T)=wrong(k+2*T)+class_test(k+2*T);end if (class_test(