表面肌电信号实验手册分析

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

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

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

小波去噪的基本原理图如下
（2） 特征提取
时域分析是将肌电信号看成均值为零，而方差随着信号强度的变化而变化的随机信号。

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

最常用的方法有：方差，过零点数（Zero Crossing, ZC ），Willison 幅值（Willison Amplitude, WAMP ），绝对值平均值 (Mean Absolute Value, MA V ）和波形长度（Wave length ，WL ）等。

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

我们主要介绍一下几种方法：
过零率（ZC ）：为波形通过零线的次数,从一定程度上反映了信号的频率特性。

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

计算方式如下：
)(),sgn(11δ≥-+-++k k k k x x x x
（1）
Willison 幅值：是由Willison 提出一种对表面肌电信号的幅值变化数量进行计算的方法，经过后人的研究，对Willison 幅值的阈值有了明确的范围限定，目前认为V μ100~50 是最合适的阈值范围。

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

∑=+-=N
t i i x x f WAMP 1
1
（2）
其中：⎩⎨
⎧>=otherwise x if x f 阈值
01
)(
波形长度（WL ）：它是对某一分析窗中的波形长度的统计，波长可以体现该样本的持续时间、幅值、频率的特征。

∑-=-+=1
1
)
()1(1N i i x i x N WL
（3）
符号改变斜率（SSC ）：为信号的的频率性能提供了一些附加信息，对于3个连续的采样点，给定阈值ω，通过下面的公式计算波峰波谷的个数。

()()()N i x x x x i i i i ,,1,11 =≥-⨯-+-ω
（4）
（3） 神经网络
BP 神经网络又称误差反向传播（Back Propagation ），它是一种多层的前向型神经网络。

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

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

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

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

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

+
1,
1IW 1
b 1
p
+
1
,2IW 2
b 1
a
隐层)(tan 111,11b p IW sig a +=输出层
)
(211,22b a IW purelin a +=y
a =3输入
由两层神经元构成的BP 网络结构
（1） BP 网络的训练算法
① BP 算法
BP 算法沿着误差函数减小最快的方向，也就是梯度的反方向改变权值和偏差，这一点与线性网络的学习算法是一致的。

BP 算法的迭代计算公式可以表示为：
k k k k g a x x -=+1 （1）
其中，k x 代表当前权值和偏差，1+k x 代表迭代产生的下一次的权值与偏差，
k g 为当前误差函数的梯度，k a 代表学习速率。

② 有动量的梯度下降算法
标准的梯度下降法在调整权值时，仅仅按照当前时刻的负梯度方向进行调整，并没有考虑以前各次运算步骤中的梯度方向，因此新的样本对迭代过程影响太大，可能会导致训练过程中调整方向发生震荡，导致不稳定和收敛速度慢的问题，有动量的梯度下降算法则考虑了往前时刻的贡献，其权值迭代算法为：
)1()()1[()()1(-+-+=+n D n D n w n w ij ij ααη （2）
其中，)1(),(-n D n D 分别表示n 时刻，n-1时刻的负梯度。

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

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

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

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

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

梯度向量为：
00g p -= （3）
沿着此方向进行权值和偏差的调整，公式为：
k k k k g a x x +=+1 （4）
下一次搜索方向则由前两次搜索方向的共轭方向决定，表达式为：
1-+-=k k k k p g p β （5）
对于系数k β不同计算方法产生不同的共轭梯度算法。

a ）F-R 共轭梯度算法采取的系数确定方法为：
1
1--=k T k k T
k k g g g
g β （6）
即本次迭代梯度相对于上一次迭代梯度的归一化值。

b ）P-R 共轭梯度算法采取的系数确定方法为：
1
11---∆=k T k k T k k g g g
g β （7）
即上次迭代梯度与本次迭代梯度的内积对本次梯度的归一化值。

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

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

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

④ L-M 算法
L-M 算法其权值和阈值的更新过程为：
e J I J J x x T T k k 11][-++-=μ （8）
其中，e为期望输出与实际输出的误差；J为误差对权值微分的Jacobi矩阵；μ为标量因子。

如果训练成功，误差性能函数减小，那么就减小μ的值；反之就减小其值。

五、实验步骤
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样本小波去噪%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%
for i=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');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');
train样本去噪前后的信号对比如下图所示
仿照对train样本进行小波去噪的方法，写出利用MATLAB对test样本进行小波去噪的程序。

（2）特征提取
clear all;
close all;
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);
for c=1:Count
for n=1:Channel
for w=1:Window-1
WL(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）动作识别
clear all;
close all;
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=[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 ;
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);
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
end
T=Gap;%80
wrong=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(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
for k=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');
hold on;
%stem(i,class(i),'b');
plot(i,wrong(i),'or'); legend('正确的类别','错误的类别');title('分类结果'); accuracy=1-wrongall/Count%分类准确率
最终的识别结果如下图所示。

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