本發明涉及一種肌電信號特征提取領域，特別涉及一種融合小波包和雙譜分析的肌電信號特征提取方法。
背景技術：
表面肌電(surface electromyography,sEMG)信號是肌肉收縮時所產生的電波動，既與肌肉本身的組織生理特性有關，也與神經控制系統有關，它反映了神經肌肉的活動和功能狀態。因此，肌電信號已廣泛應用于生理醫學、康復醫學及運動醫學等領域的研究，并且成為驅動機器人、控制假肢運動以及功能性電刺激的理想控制信號。
sEMG具有非平穩特性，同時，在采集sEMG信號時易受環境影響而引入高斯噪聲。單采用小波變換、AR模型或小波包變換等傳統方法進行特征提取時，其正確識別率依然不高，導致指令控制智能硬件時出現誤動作。鑒于此，本文結合小波包的時頻分析能力和雙譜對高斯噪聲不敏感的優勢，利用小波包對sEMG的各頻帶進行分解，得到sEMG在不同頻帶內的能量分布，然后通過對小波包重構的時域信號分析，提取出雙譜特征，最后結合這兩種方法提取出特征向量。不僅能得到sEMG信號的特征信息，還能抑制噪聲，對sEMG信號識別具有很重要的意義。
技術實現要素：
本發明旨在解決以上現有技術的問題。提出了一種提高sEMG信號的識別率的融合小波包和雙譜分析的肌電信號特征提取方法。本發明的技術方案如下：
一種融合小波包和雙譜分析的肌電信號特征提取方法，其特征在于，包括以下步驟：
S1、輸入表面肌電sEMG信號，采用小波包對sEMG信號進行j層分解；S2、選取n‘個有效頻帶，求出每個頻帶的能量；S3、將步驟S2求出的頻帶能量進行歸一化處理，作為表面肌電sEMG信號的特征向量；S4、小波包重構sEMG信號的時域信號；S5、對表面肌電sEMG的時域信號進行雙譜分析，提取出雙譜特征；S6、融合小波包和雙譜分析的肌電信號特征向量。
步驟S1中采用小波包對表面肌電sEMG信號進行j層分解，小波包對sEMG進行分解的算法為：
式中為上層小波包分解結果，和為下一級分解結果，j是尺度指標，l為位置指標，n為頻率指標，k為變量，h0和h1為分解使用的多分辨率濾波器系數。
進一步的，所述步驟S2：選取n個有效頻帶，求出每個頻帶的能量具體包括步驟：假設信號的采樣頻率為2f，如果對信號進行j層小波包分解，就可以形成2j個等寬頻帶，每個區間頻寬為f/(2j)，小波包分解之后，得到j層小波包系數Cj,m,k，k＝0,1,…,2j-1,m為小波包空間位置標識；
由Parseval能量積分等式，sEMG信號x(t)在時域上的能量為：
根據該能量積分等式關聯起來可知，式子與x(t)的小波包變換系數Cj,k具有能量量綱，因此原始信號的能量分布用小波包能量譜來表示是可行的。
進一步的，將步驟S2求出的頻帶能量進行歸一化處理的公式為：
Tj,k表示頻帶能量；T′j,k表示歸一化處理的頻帶能量。
并將其當作sEMG信號的特征向量，即：
T＝[T′1,T′2,…,T′n]。
T＝[T′1,T′2,…,T′n]。
進一步的，步驟S4采用小波包重構表面肌電sEMG的時域信號，其中小波包重構算法為：
式中為上層小波包分解結果，和為下一級分解結果，j是尺度指標，l為位置指標，n為頻率指標，k為變量，h0和h1為分解使用的多分辨率濾波器系數。
進一步的，S5、對表面肌電sEMG的時域信號進行雙譜分析，提取出雙譜特征，具體包括步驟：
設一離散、平穩時間序列{x(i)}，t＝1,2,…,N，根據信號的離散傅里葉變換，雙譜定義為：
Bx(f1,f2)＝E[X(f1)X(f2)X*(f1+f2)]，
式中fi為頻率變量；X(fi)為信號的離散傅里葉變換；i＝1,2；E[]為數學期望，雙譜估計采用雙譜估計間接法，首先估計實驗數據{x(i)}的三階累積量，然后對該累積量進行二維DFT變換，這樣就得到隨機序列的雙譜估計。
進一步的，所述雙譜估計間接法算法的具體描述如下：
a、將長度為N的實驗數據{x(i)}分成K段，每段有M個數據，即N＝KM；
b、進行去均值操作，去除每段數據的均值，使將分析數據的均值為0；
c、假設{xj(i)}(i＝1,2,…,M；j＝1,2,…,K)為第j段的數據，估計每段數據的三階累積量：
表示第j段數據的三階積累量，
k1＝max{1,-m,-n}，k1表示1、-m、-n中最大值；
k2＝min{M,M-m,M-n}；k2表示M、M-m、M-n中最小值；
d、對進行統計平均，得到K組數據的累積量估計，即：
e、對三階累積量估計進行二維DFT變換，即得雙譜估計。
進一步的，所述融合小波包和雙譜分析的肌電信號特征向量，其具體步驟為：
每個樣本的雙譜值為Bj(f1,f2)，j＝1,2,…,n，n為樣本數，對每個樣本分別進行雙譜計算，取平均可得樣本雙譜為：
接下來再對每個樣本求極大值Max(Bj(f1,f2))所對應的即為特征向量，
然后與所述能量特征組合作為新的特征向量：
本發明的優點及有益效果如下：
本發明針對表面肌電信號(sEMG)的非平穩特性導致信號的正確識別率低的問題，提出一種融合小波包和雙譜分析相結合的特征提取方法。首先，采用小波包對sEMG的各頻帶進行分解，得到sEMG在不同頻帶內的能量分布；其次，通過對小波包重構的時域信號進行雙譜分析，提取出雙譜特征；最后，融合這兩種方法提取出特征向量。這種融合小波包和雙譜分析的肌電信號特征提取方法不僅能提取到高時頻分辨率的sEMG信號特征信息，還能抑制噪聲，從而提高sEMG信號的識別率。
附圖說明
圖1是本發明提供優選實施例融合小波包和雙譜分析的肌電信號特征提取方法流程圖。
具體實施方式
下面將結合本發明實施例中的附圖，對本發明實施例中的技術方案進行清楚、詳細地描述。所描述的實施例僅僅是本發明的一部分實施例。
本發明解決上述技術問題的技術方案是，
圖1所示，本發明提供了一種融合小波包和雙譜分析的肌電信號特征提取方法，其包括以下步驟：
第一步，采用小波包對sEMG進行j層分解。小波包對sEMG進行分解算法為：
式中為上層小波包分解結果，和為下一級分解結果，j是尺度指標，l為位置指標，n為頻率指標，k為變量，h0和h1為分解使用的多分辨率濾波器系數。
第二步，選取n個有效頻帶，求出每個頻帶的能量。假設信號的采樣頻率為2f，如果對信號進行j層小波包分解，就可以形成2j個等寬頻帶，每個區間頻寬為f/(2j)。小波包分解之后，得到j層小波包系數Cj,m,k，k＝0,1,…,2j-1,m為小波包空間位置標識。
由Parseval能量積分等式，sEMG信號x(t)在時域上的能量為：
根據該能量積分等式關聯起來可知，式子與x(t)的小波包變換系數Cj,k具有能量量綱，因此原始信號的能量分布用小波包能量譜來表示是可行的。
第三步，將上述的頻帶能量進行歸一化處理：
并將其當作sEMG信號的特征向量，即：
T＝[T′1,T′2,…,T′n]。
T＝[T′1,T′2,…,T′n]。
第四步，小波包重構sEMG的時域信號。小波包重構算法為：
式中為上層小波包分解結果，和為下一級分解結果，j是尺度指標，l為位置指標，n為頻率指標，k為變量，h0和h1為分解使用的多分辨率濾波器系數。
第五步，對重構的sEMG的時域信號進行雙譜分析，提取出雙譜特征。具體算法和步驟為：
設一離散、平穩時間序列{x(i)}，t＝1,2,…,N，根據信號的離散傅里葉變換，雙譜定義為：
Bx(f1,f2)＝E[X(f1)X(f2)X*(f1+f2)]，
式中fi為頻率變量；X(fi)為信號的離散傅里葉變換；i＝1,2；E[]為數學期望。雙譜估計有直接法和間接法，本文采用雙譜估計間接法，其核心思想是首先估計實驗數據{x(i)}的三階累積量，然后對該累積量進行二維DFT變換，這樣就得到隨機序列的雙譜估計。該算法的具體描述如下：
a、將長度為N的實驗數據{x(i)}分成K段，每段有M個數據，即N＝KM；
b、進行去均值操作，去除每段數據的均值，使將分析數據的均值為0；
c、假設{xj(i)}(i＝1,2,…,M；j＝1,2,…,K)為第j段的數據，估計每段數據的三階累積量：
k1＝max{1,-m,-n}，
k2＝min{M,M-m,M-n}；
d、對進行統計平均，得到K組數據的累積量估計，即：
e、對三階累積量估計進行二維DFT變換，即得雙譜估計。
第六步，融合小波包和雙譜分析的肌電信號特征向量。其具體步驟為：
每個樣本的雙譜值為Bj(f1,f2)，j＝1,2,…,n，n為樣本數，對每個樣本分別進行雙譜計算，取平均可得樣本雙譜為：
接下來再對每個樣本求極大值Max(Bj(f1,f2))所對應的即為特征向量。
然后與步驟4中的能量特征組合作為新的特征向量：
以上這些實施例應理解為僅用于說明本發明而不用于限制本發明的保護范圍。在閱讀了本發明的記載的內容之后，技術人員可以對本發明作各種改動或修改，這些等效變化和修飾同樣落入本發明權利要求所限定的范圍。