表面肌電信號試驗手冊總結

1、實驗 基于sEMG時域特征特的動作識別一、實驗目的, 了解肌電信號常用的時域分析方法；.利用MATLAB對肌電信號進行去噪、特征提取及動作識別；二、實驗設備Wi-Fi表面肌電信號采集卡；32 位 Windows XP 臺式機（Matlab 7.0 軟件）；802.11b/g無線網卡；三、實驗內容（1）學習信號的基本去噪方法，并用 MATLAB實現；（2）學習肌電信號常用的時域特征并利用 Matlab來進行波形長度（WL） 符號改變數（SSC）、過零點（ZC）、威爾遜賦值（WAMP）等特征的提?。唬?）學習神經網絡信號處理方法，掌握 BP神經網絡的用法，將其用于肌 電信號的動作識別。學習以上三個

2、部分，最終完成一整套肌電信號去噪、 特征提?。ㄟx取一種特 征）、基于特征的動作識別的 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 ):它是對某一分析窗中的波形長度的統(tǒng)計，波長可以體現該樣本的持續(xù)時間、幅值、頻率的特征。1 N 1WL=? x(i 1)-x(i)(3)N 1符號改變斜率(SSC：為信號的的頻率性能提供了一些附加信息，對于 3 個連續(xù)的采樣點，給定閾值 ，通過下面的公式計算波峰波谷的個數。(4)為-為二 X -X 1 一，i

5、=1, ,N(3)神經網絡BP神經網絡又稱誤差反向傳播(Back Propagation),它是一種多層的前向 型神經網絡。在BP網絡中，信號是前向傳播的，而誤差是反向傳播的。所謂的 反向傳播是指誤差的調整過程是從最后的輸出層依次向之前各層逐漸進行的。標準的BP網絡采用梯度下降算法，與 Widrow-Hoff學習規(guī)則相似，網絡權值沿著 性能函數的梯度反向調整。生理信號采集與分析處理平臺實驗指導-武漢理工大學信息工程學院前向型神經網絡通常具有一個或多個由 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)整，并沒有考慮以前各次運算步驟中的梯度方向， 因此新的樣本對迭代過程影響 太大，可能會導致訓練過程中調整方向發(fā)生震蕩， 導致不穩(wěn)定和收斂速度慢的問 題，有動量的梯度下降算法則考慮了往前時刻的貢獻，其權值迭代算法為：wj (n 1) = wj (n)(1 - )D(n)= D(n -1)其中，D(n), D(n -1)分別表示n時刻，n-1時刻的負梯度。由于加入了以前時刻梯度的貢獻，相當于給迭代過程添加了一個低通濾波器， 使得網絡忽略誤差生理信號采集與分析處理平臺實驗指導-武漢理工大學信息工程學院曲面上細節(jié)特征，避免了陷入局部極小點的問題共腕梯度算法盡管標準的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維行向量，其元素為網絡中各層節(jié)點數；生理信號采集與分析處理平臺實驗指導-武漢理工大學信息工程學院%

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