表面肌實驗報告

1、武漢理工大學(xué)現(xiàn)代數(shù)字信號處理在前沿學(xué)科中的應(yīng)用實驗報告基于sEMG時域特征的動作識別學(xué)院： 信息工程學(xué)院 學(xué)號： 1049731503279 姓名： 吳志勇 班級： 電子154 實驗 基于sEMG時域特征特的動作識別1、 實驗?zāi)康?. 了解肌電信號常用的時域分析方法；2. 利用MATLAB對肌電信號進行去噪、特征提取及動作識別；2、 實驗設(shè)備1. Wi-Fi表面肌電信號采集卡；2. 32位Windows XP臺式機（Matlab 7.0軟件）；3. 802.11b/g無線網(wǎng)卡；3、 實驗內(nèi)容（1）學(xué)習(xí)信號的基本去噪方法，并用MATLAB實現(xiàn)；（2）學(xué)習(xí)肌電信號常用的時域特征并利用Matlab來

2、進行波形長度（WL）符號改變數(shù)（SSC）、過零點（ZC）、威爾遜賦值（WAMP）等特征的提取；（3） 學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)信號處理方法，掌握BP神經(jīng)網(wǎng)絡(luò)的用法，將其用于肌電信號的動作識別。學(xué)習(xí)以上三個部分，最終完成一整套肌電信號去噪、特征提取（選取一種特征）、基于特征的動作識別的MATLAB程序。4、 實驗原理（1） 小波去噪小波去噪方法是一種建立在小波變換基礎(chǔ)上的新興算法，基本思想是根據(jù)噪聲在不同頻帶上的小波分解系數(shù)具有不同強度分布的特點，將各頻帶上的噪聲對應(yīng)的小系數(shù)去除，保留原始信號的小波分解系數(shù)，然后對處理后系數(shù)進行小波重構(gòu)，得到純凈信號。小波去噪的基本原理圖如下（2） 特征提取時域分析是將肌電

3、信號看成均值為零，而方差隨著信號強度的變化而變化的隨機信號。時域特征的計算復(fù)雜度低，提取比較方便。最常用的方法有：方差，過零點數(shù)（Zero Crossing, ZC），Willison幅值（Willison Amplitude, WAMP），絕對值平均值 (Mean Absolute Value, MAV）和波形長度（Wave length，WL）等。在實際應(yīng)用中，為了讓特征可以包含更多的信息，往往選擇用不同的時域特征組合形成聯(lián)合特征向量。我們主要介紹一下幾種方法：過零率（ZC）：為波形通過零線的次數(shù),從一定程度上反映了信號的頻率特性。為了降低零點引入的噪聲，往往會引入一個閾值。計算方式如下：

4、（1）Willison幅值：是由Willison提出一種對表面肌電信號的幅值變化數(shù)量進行計算的方法，經(jīng)過后人的研究，對Willison幅值的閾值有了明確的范圍限定，目前認為 是最合適的閾值范圍。其數(shù)學(xué)表示公式如公式（3-3）。（2）其中：波形長度（WL）：它是對某一分析窗中的波形長度的統(tǒng)計，波長可以體現(xiàn)該樣本的持續(xù)時間、幅值、頻率的特征。（3）符號改變斜率（SSC）：為信號的的頻率性能提供了一些附加信息，對于3個連續(xù)的采樣點，給定閾值，通過下面的公式計算波峰波谷的個數(shù)。（4）（3） 神經(jīng)網(wǎng)絡(luò)BP神經(jīng)網(wǎng)絡(luò)又稱誤差反向傳播（Back Propagation），它是一種多層的前向型神經(jīng)網(wǎng)絡(luò)。在BP網(wǎng)

5、絡(luò)中，信號是前向傳播的，而誤差是反向傳播的。所謂的反向傳播是指誤差的調(diào)整過程是從最后的輸出層依次向之前各層逐漸進行的。標準的BP網(wǎng)絡(luò)采用梯度下降算法，與Widrow-Hoff學(xué)習(xí)規(guī)則相似，網(wǎng)絡(luò)權(quán)值沿著性能函數(shù)的梯度反向調(diào)整。前向型神經(jīng)網(wǎng)絡(luò)通常具有一個或多個由sigmoid神經(jīng)元構(gòu)成的隱層，以及一個由線性神經(jīng)元構(gòu)成的輸出層。多個具有非線性傳遞函數(shù)的神經(jīng)元層使得網(wǎng)絡(luò)可以學(xué)習(xí)輸入和輸出之間的非線性關(guān)系，而線性輸出層使得網(wǎng)絡(luò)可以產(chǎn)生-1，+1之外的輸出值。由兩層神經(jīng)元構(gòu)成的BP網(wǎng)絡(luò)結(jié)構(gòu)（1） BP網(wǎng)絡(luò)的訓(xùn)練算法 BP算法BP算法沿著誤差函數(shù)減小最快的方向，也就是梯度的反方向改變權(quán)值和偏差，這一點與線性

6、網(wǎng)絡(luò)的學(xué)習(xí)算法是一致的。BP算法的迭代計算公式可以表示為： （1）其中，代表當前權(quán)值和偏差，代表迭代產(chǎn)生的下一次的權(quán)值與偏差，為當前誤差函數(shù)的梯度，代表學(xué)習(xí)速率。 有動量的梯度下降算法標準的梯度下降法在調(diào)整權(quán)值時，僅僅按照當前時刻的負梯度方向進行調(diào)整，并沒有考慮以前各次運算步驟中的梯度方向，因此新的樣本對迭代過程影響太大，可能會導(dǎo)致訓(xùn)練過程中調(diào)整方向發(fā)生震蕩，導(dǎo)致不穩(wěn)定和收斂速度慢的問題，有動量的梯度下降算法則考慮了往前時刻的貢獻，其權(quán)值迭代算法為： （2）其中，分別表示n時刻，n-1時刻的負梯度。由于加入了以前時刻梯度的貢獻，相當于給迭代過程添加了一個低通濾波器，使得網(wǎng)絡(luò)忽略誤差曲面上細節(jié)特

7、征，避免了陷入局部極小點的問題。 共軛梯度算法盡管標準的BP算法采用梯度下降算法，權(quán)值和偏差沿誤差函數(shù)下降最快的方向調(diào)整，但卻并不一定是收斂最快的算法。在改進的BP訓(xùn)練算法中，有一大類的算法稱為共軛梯度算法。在這一類算法中，權(quán)值和偏差沿著共軛梯度方向進行調(diào)整，通常能夠獲得比標準的梯度算法更快的收斂速度。共軛梯度算法的第一次迭代都是從最陡下降的梯度方向開始。梯度向量為： （3）沿著此方向進行權(quán)值和偏差的調(diào)整，公式為： （4）下一次搜索方向則由前兩次搜索方向的共軛方向決定，表達式為： （5）對于系數(shù)不同計算方法產(chǎn)生不同的共軛梯度算法。a）F-R共軛梯度算法采取的系數(shù)確定方法為： （6）即本次迭代梯

8、度相對于上一次迭代梯度的歸一化值。b）P-R共軛梯度算法采取的系數(shù)確定方法為： （7）即上次迭代梯度與本次迭代梯度的內(nèi)積對本次梯度的歸一化值。c）Scaled共軛梯度算法到目前為止，討論過的所有共軛梯度算法都需要在每一步迭代過程中對搜索方向進行計算，這樣的計算量是比較大。對此moller提出了Scaled梯度搜索算法4，在每一步迭代過程中不計算搜索方向，以減少訓(xùn)練過程的計算量。其基本原理是利用下面介紹的L-M算法與共軛梯度法相結(jié)合產(chǎn)生的。 L-M算法L-M算法其權(quán)值和閾值的更新過程為： （8）其中，e為期望輸出與實際輸出的誤差；J為誤差對權(quán)值微分的Jacobi矩陣；為標量因子。如果訓(xùn)練成功，誤

9、差性能函數(shù)減小，那么就減小的值；反之就減小其值。5、 實驗步驟（1） 數(shù)據(jù)格式轉(zhuǎn)換：N=4;%通道的個數(shù)M=512;%每個通道的采樣數(shù)K=10;%每種類別所取的個數(shù)P=6;%類別數(shù)ClassPermode=10;%每類動作數(shù)據(jù)的個數(shù)ClassNumber=6;%類別的個數(shù)M=512;%每個樣本點數(shù)據(jù)的個數(shù)Channel=4;%通道的個數(shù)Count=ClassPermode*ClassNumber*M;sample_train=zeros(Count,Channel);%保存訓(xùn)練樣本sample_test=zeros(Count,Channel);%保存訓(xùn)練樣本for i=1:60 D=impo

10、rtdata(strcat('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_tra

11、in sample_trainsave sample_test sample_test程序運行后，會看到文件夾中多出了sample_train和sample_test兩個.mat文件。（2） 小波去噪及特征提取，選用WAMP特征的提取。load 'sample_test.mat'load 'sample_train.mat'% 參數(shù)說明 %Window=256;%分析窗口的長度M=512; %采集數(shù)據(jù)時一個data的樣本數(shù)Channel=4; %采集數(shù)據(jù)的通道數(shù)Class=6; %類別數(shù)Number=10; %每個類別的個數(shù)WinLap=64; %窗口移動的間隔

12、JudgeTime=Window/WinLap; %一個分析窗口需要移動的次數(shù) 8Count=M*Number*Class/WinLap-Window/WinLap+1; %所有數(shù)據(jù)需要分析的次數(shù) 477ClassCount=M*Number/WinLap-Window/WinLap+1; %一類數(shù)據(jù)需要的分析次數(shù) 77GapCount=M*Number/WinLap; %訓(xùn)練樣本兩類動作之間的間隔 80ClassOne=1;ClassTwo=2;ClassThree=3;ClassFour=4;ClassFive=5;ClassSix=6;%thr=0.2;% train樣本小波去噪 %fo

13、r 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;endfi

14、gure(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)

15、,'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=1:Window-1 if(abs(new_train(c-1)*WinLap+w,n)-new_train(c-1)*W

16、inLap+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); d3=detcoef(c,l,3

17、); 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樣本濾波前后信號對比');subpl

18、ot(4,1,2);plot(sample_test(:,2),'r');hold on;plot(new_test(:,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=

19、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)*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 '

20、;feature.mat'Gap=80;K=77;Count=477;%共有477組特征值P=6;%6類動作Feature=Feature_test;%統(tǒng)一接口test=Feature_train;%統(tǒng)一接口net1=newff(minmax(Feature),7,6,'tansig','purelin','trainlm');%net = newff ( A, B, C ,trainFun)% % 參數(shù)：% % A：一個n×2的矩陣，第i行元素為輸入信號xi的最小值和最大值；% B：一個k維行向量，其元素為網(wǎng)絡(luò)中各層節(jié)點數(shù)；%

21、 C：一個k維字符串行向量，每一分量為對應(yīng)層神經(jīng)元的激活函數(shù)；常用的有：% 線性函數(shù)purelin，對數(shù)S形轉(zhuǎn)移函數(shù)logsig 雙曲正切S形函數(shù) tansig% trainFun ：為學(xué)習(xí)規(guī)則采用的訓(xùn)練算法。也可以根據(jù)需要修改%設(shè)置訓(xùn)練的輸出目標矩陣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 ze

22、0_6 ze0 ze0 ze0 ze0 ze1 ze0_6 ze0 ze0 ze0 ze0 ze0 ze1_6;%神經(jīng)網(wǎng)絡(luò)進行訓(xùn)練tic ; net1.trainParam.show=200; net1.trainParam.epochs=10000; net1.trainParam.goal=0.1; net1.trainParam.lr=0.01; net1,tr=train(net1,Feature,t); Y=sim(net1,Feature); YS=sim(net1,test); toc; % 一些重要的網(wǎng)絡(luò)配置參數(shù)如下：% net.trainparam.goal ：神經(jīng)網(wǎng)絡(luò)訓(xùn)練的

23、目標誤差% net.trainparam.show ： 顯示中間結(jié)果的周期% net.trainparam.epochs ：最大迭代次數(shù)% net.trainParam.lr ： 學(xué)習(xí)率% 語法： net, tr = train( net, X, Y )% 參數(shù)：% X：網(wǎng)絡(luò)實際輸入% Y：網(wǎng)絡(luò)應(yīng)有輸出% tr：訓(xùn)練跟蹤信息% Y1：網(wǎng)絡(luò)實際輸出% 語法：Y=sim(net,X)% net：網(wǎng)絡(luò)% X：輸入給網(wǎng)絡(luò)的×N矩陣，其中K為網(wǎng)絡(luò)輸入個數(shù)，N為數(shù)據(jù)樣本數(shù)% Y：輸出矩陣Q×N，其中Q為網(wǎng)絡(luò)輸出個數(shù)%判斷輸出結(jié)果為第幾類 maxY=max(Y,1); for i=1:C

24、ount 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(k+3*T)