表面肌电信号实验手册.docx

表面肌电信号实验手册.docx

《表面肌电信号实验手册.docx》由会员分享，可在线阅读，更多相关《表面肌电信号实验手册.docx（14页珍藏版）》请在冰豆网上搜索。

表面肌电信号实验手册

实验基于sEMG时域特征特的动作识别

1、实验目的

1.了解肌电信号常用的时域分析方法；

2.利用MATLAB对肌电信号进行去噪、特征提取及动作识别；

2、实验设备

1.Wi-Fi表面肌电信号采集卡；

2.32位WindowsXP台式机（Matlab7.0软件）；

3.802.11b/g无线网卡；

3、实验内容

（1）学习信号的基本去噪方法，并用MATLAB实现；

（2）学习肌电信号常用的时域特征并利用Matlab来进行波形长度（WL）符号改变数（SSC）、过零点（ZC）、威尔逊赋值（WAMP）等特征的提取；

（3）学习神经网络信号处理方法，掌握BP神经网络的用法，将其用于肌电信号的动作识别。

学习以上三个部分，最终完成一整套肌电信号去噪、特征提取（选取一种特征）、基于特征的动作识别的MATLAB程序。

4、实验原理

（1）小波去噪

小波去噪方法是一种建立在小波变换基础上的新兴算法，基本思想是根据噪声在不同频带上的小波分解系数具有不同强度分布的特点，将各频带上的噪声对应的小系数去除，保留原始信号的小波分解系数，然后对处理后系数进行小波重构，得到纯净信号。

小波去噪的基本原理图如下

（2）特征提取

时域分析是将肌电信号看成均值为零，而方差随着信号强度的变化而变化的随机信号。

时域特征的计算复杂度低，提取比较方便。

最常用的方法有：

方差，过零点数（ZeroCrossing,ZC），Willison幅值（WillisonAmplitude,WAMP），绝对值平均值（MeanAbsoluteValue,MAV）和波形长度（Wavelength，WL）等。

在实际应用中，为了让特征可以包含更多的信息，往往选择用不同的时域特征组合形成联合特征向量。

我们主要介绍一下几种方法：

过零率（ZC）：

为波形通过零线的次数,从一定程度上反映了信号的频率特性。

为了降低零点引入的噪声，往往会引入一个阈值δ。

计算方式如下：

（1）

Willison幅值：

是由Willison提出一种对表面肌电信号的幅值变化数量进行计算的方法，经过后人的研究，对Willison幅值的阈值有了明确的范围限定，目前认为

是最合适的阈值范围。

其数学表示公式如公式（3-3）。

（2）

其中：

波形长度（WL）：

它是对某一分析窗中的波形长度的统计，波长可以体现该样本的持续时间、幅值、频率的特征。

（3）

符号改变斜率（SSC）：

为信号的的频率性能提供了一些附加信息，对于3个连续的采样点，给定阈值ω，通过下面的公式计算波峰波谷的个数。

（4）

（3）神经网络

BP神经网络又称误差反向传播（BackPropagation），它是一种多层的前向型神经网络。

在BP网络中，信号是前向传播的，而误差是反向传播的。

所谓的反向传播是指误差的调整过程是从最后的输出层依次向之前各层逐渐进行的。

标准的BP网络采用梯度下降算法，与Widrow-Hoff学习规则相似，网络权值沿着性能函数的梯度反向调整。

前向型神经网络通常具有一个或多个由sigmoid神经元构成的隐层，以及一个由线性神经元构成的输出层。

多个具有非线性传递函数的神经元层使得网络可以学习输入和输出之间的非线性关系，而线性输出层使得网络可以产生[-1，+1]之外的输出值。

由两层神经元构成的BP网络结构

（1）BP网络的训练算法

①BP算法

BP算法沿着误差函数减小最快的方向，也就是梯度的反方向改变权值和偏差，这一点与线性网络的学习算法是一致的。

BP算法的迭代计算公式可以表示为：

（1）

其中，

代表当前权值和偏差，

代表迭代产生的下一次的权值与偏差，

为当前误差函数的梯度，

代表学习速率。

②有动量的梯度下降算法

标准的梯度下降法在调整权值时，仅仅按照当前时刻的负梯度方向进行调整，并没有考虑以前各次运算步骤中的梯度方向，因此新的样本对迭代过程影响太大，可能会导致训练过程中调整方向发生震荡，导致不稳定和收敛速度慢的问题，有动量的梯度下降算法则考虑了往前时刻的贡献，其权值迭代算法为：

（2）

其中，

分别表示n时刻，n-1时刻的负梯度。

由于加入了以前时刻梯度的贡献，相当于给迭代过程添加了一个低通滤波器，使得网络忽略误差曲面上细节特征，避免了陷入局部极小点的问题。

③共轭梯度算法

尽管标准的BP算法采用梯度下降算法，权值和偏差沿误差函数下降最快的方向调整，但却并不一定是收敛最快的算法。

在改进的BP训练算法中，有一大类的算法称为共轭梯度算法。

在这一类算法中，权值和偏差沿着共轭梯度方向进行调整，通常能够获得比标准的梯度算法更快的收敛速度。

共轭梯度算法的第一次迭代都是从最陡下降的梯度方向开始。

梯度向量为：

（3）

沿着此方向进行权值和偏差的调整，公式为：

（4）

下一次搜索方向则由前两次搜索方向的共轭方向决定，表达式为：

（5）

对于系数

不同计算方法产生不同的共轭梯度算法。

a）F-R共轭梯度算法采取的系数确定方法为：

（6）

即本次迭代梯度相对于上一次迭代梯度的归一化值。

b）P-R共轭梯度算法采取的系数确定方法为：

（7）

即上次迭代梯度与本次迭代梯度的内积对本次梯度的归一化值。

c）Scaled共轭梯度算法

到目前为止，讨论过的所有共轭梯度算法都需要在每一步迭代过程中对搜索方向进行计算，这样的计算量是比较大。

对此moller提出了Scaled梯度搜索算法[4]，在每一步迭代过程中不计算搜索方向，以减少训练过程的计算量。

其基本原理是利用下面介绍的L-M算法与共轭梯度法相结合产生的。

④L-M算法

L-M算法其权值和阈值的更新过程为：

（8）

其中，e为期望输出与实际输出的误差；J为误差对权值微分的Jacobi矩阵；

为标量因子。

如果训练成功，误差性能函数减小，那么就减小

的值；反之就减小其值。

5、实验步骤

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;%窗口移动的间隔

JudgeTime=Window/WinLap;%一个分析窗口需要移动的次数8

Count=M*Number*Class/WinLap-Window/WinLap+1;%所有数据需要分析的次数477

ClassCount=M*Number/WinLap-Window/WinLap+1;%一类数据需要的分析次数77

GapCount=M*Number/WinLap;%训练样本两类动作之间的间隔80

ClassOne=1;

ClassTwo=2;

ClassThree=3;

ClassFour=4;

ClassFive=5;

ClassSix=6;

%thr=0.2;

%%%%%train样本小波去噪%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fori=1:

Channel

x=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;

end

figure

（1）

subplot（4,1,1）;plot（sample_train（:

1）,'r'）;holdon;plot（new_train（:

1）,'b'）;legend（'原始信号','去噪后信号'）;title（'train样本滤波前后信号对比'）;

subplot（4,1,2）;plot（sample_train（:

2）,'r'）;holdon;plot（new_train（:

2）,'b'）;

subplot（4,1,3）;plot（sample_train（:

3）,'r'）;holdon;plot（new_train（:

3）,'b'）;

subplot（4,1,4）;plot（sample_train（:

4）,'r'）;holdon;plot（new_train（:

4）,'b'）;

train样本去噪前后的信号对比如下图所示

仿照对train样本进行小波去噪的方法，写出利用MATLAB对test样本进行小波去噪的程序。

（2）特征提取

clearall;

closeall;

load'sample_train.mat';

%%%%%参数说明%%%%%%

Window=256;%分析窗口的长度

M=512;%采集数据时一个data的样本数

Channel=4;%采集数据的通道数

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=sample_train;%接口

WL=zeros（Channel,Count）;

forc=1:

Count

forn=1:

Channel

forw=1:

Window-1WL（n,c）=WL（n,c）+abs（sample（（c-1）*WinLap+w,n）-sample（（c-1）*WinLap+w+1,n））;

end

end

end

Feature_train=WL;

仿照求WL（波长）特征的程序，任选其他一种特征，写出MATLAB程序。

注意，train样本和test样本都需要做特征提取。

将最终求得的train和test样本的特征存入feature.mat。

（3）动作识别

clearall;

closeall;

load'feature.mat'

Gap=80;

K=77;

Count=477;%共有477组特征值

P=6;%6类动作

Feature=Feature_test;%统一接口

test=Feature_train;%统一接口

net=newff（minmax（Feature）,[5,6],{'tansig','purelin'},'trainlm'）;

%net=newff（A,B,{C},‘trainFun’）

%

%参数：

%

%A：

一个n×2的矩阵，第i行元素为输入信号xi的最小值和最大值；

%B：

一个k维行向量，其元素为网络中各层节点数；

%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=[ze1ze0ze0ze0ze0ze0_6

ze0ze1ze0ze0ze0ze0_6

ze0ze0ze1ze0ze0ze0_6

ze0ze0ze0ze1ze0ze0_6

ze0ze0ze0ze0ze1ze0_6

ze0ze0ze0ze0ze0ze1_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.show：

显示中间结果的周期

%net.trainparam.epochs　：

最大迭代次数

%net.trainParam.lr：

学习率

%语法：

[net,tr]=train（net,X,Y）

%参数：

%X：

网络实际输入

%Y：

网络应有输出

%tr：

训练跟踪信息

%Y1：

网络实际输出

%语法：

Y=sim（net,X）

%net：

网络

%X：

输入给网络的Ｋ×N矩阵，其中K为网络输入个数，N为数据样本数

%Y：

输出矩阵Q×N，其中Q为网络输出个数

%判断输出结果为第几类

maxY=max（Y,[],1）;

fori=1:

Count

forj=1:

P

if（maxY（i）==Y（j,i））class（i）=j;end

end

end

maxYS=max（YS,[],1）;

fori=1:

Count

forj=1:

P

if（maxYS（i）==YS（j,i））class_test（i）=j;end

end

end

T=Gap;%80

wrong=zeros（1,Count）;

fork=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（k+3*T）~=4）

wrong（k+3*T）=wrong（k+3*T）+class_test（k+3*T）;end

if（class_test（k+4*T）~=5）

wrong（k+4*T）=wrong（k+4*T）+class_test（k+4*T）;end

end

fork=1:

K

if（class_test（k+5*T）~=6）

wrong（k+5*T）=wrong（k+5*T）+class_test（k+5*T）;end

end

wrongall=sum（wrong）;

figure

（1）;

i=1:

Count;

plot（i,class_test（i）,'ob'）;

holdon;

%stem（i,class（i）,'b'）;

plot（i,wrong（i）,'or'）;legend（'正确的类别','错误的类别'）;title（'分类结果'）;

accuracy=1-wrongall/Count%分类准确率

最终的识别结果如下图所示。

改变BP网络隐含层节点数以及初始化各层之间的权值和阈值，可以提高最终的识别率。