JP2017104333A - 筋電信号処理方法、装置およびプログラム - Google Patents

筋電信号処理方法、装置およびプログラム Download PDF

Info

Publication number
JP2017104333A
JP2017104333A JP2015241474A JP2015241474A JP2017104333A JP 2017104333 A JP2017104333 A JP 2017104333A JP 2015241474 A JP2015241474 A JP 2015241474A JP 2015241474 A JP2015241474 A JP 2015241474A JP 2017104333 A JP2017104333 A JP 2017104333A
Authority
JP
Japan
Prior art keywords
signal
variance
myoelectric signal
myoelectric
distribution
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
JP2015241474A
Other languages
English (en)
Other versions
JP6652252B2 (ja
Inventor
辻 敏夫
Toshio Tsuji
敏夫 辻
雄一 栗田
Yuichi Kurita
雄一 栗田
英朗 早志
Hideaki Hayashi
英朗 早志
彬 古居
Akira Furui
彬 古居
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hiroshima University NUC
Original Assignee
Hiroshima University NUC
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hiroshima University NUC filed Critical Hiroshima University NUC
Priority to JP2015241474A priority Critical patent/JP6652252B2/ja
Publication of JP2017104333A publication Critical patent/JP2017104333A/ja
Application granted granted Critical
Publication of JP6652252B2 publication Critical patent/JP6652252B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Abstract

【課題】整流平滑化後の筋電信号から信号強度依存ノイズが重畳された整流平滑化前の筋電信号の分散の分布を推定し、さらに人工筋電信号を生成する。【解決手段】筋電信号処理方法は、信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算し、整流平滑化信号の平均および筋電信号の整流平滑化に係るフィルタゲインの逆数から筋電信号の分散の平均の推定値を計算し、筋電信号の分散の平均の推定値および筋電信号の分散の分布に係るパラメータから信号強度依存ノイズの分散の推定値を計算し、白色正規乱数を発生させ、白色正規乱数の周波数特性を筋電信号と同様の周波数特性に整形し、筋電信号の分散の平均の推定値および信号強度依存ノイズの分散の推定値で規定される筋電信号の分散の分布に従った乱数を発生させ、整形された乱数および筋電信号の分散の分布に従った乱数から人工筋電信号を生成する。【選択図】図1

Description

本発明は、表面筋電位信号(以下、筋電信号という)の処理技術に関し、特に、整流平滑化後の筋電信号から信号強度依存ノイズが重畳された整流平滑化前の筋電信号の分散の分布を推定し、さらに当該推定した分布に基づいて人工筋電信号を生成する技術に関する。
筋電信号は筋活動の定量評価に有用であり、従来、筋電義手の制御や、リハビリ時やトレーニング時の運動解析に応用されてきた。それらの応用においては、筋電信号から特徴量を抽出すべく、筋電信号がモデル化され用いられる(例えば、非特許文献1参照)。
近年、筋電信号を含む生体信号の研究が進むにつれ、信号強度依存ノイズの存在が明らかになってきた(例えば、非特許文献2参照)。すなわち、腕に思いきり力を入れると腕が震えるように、筋電信号には筋力の強さに応じた信号強度依存ノイズが重畳される。
N. Hogan and R. W. Mann, "Myoelectric signal processing: Optimal estimation applied to electromyography-part i: Derivation of the optimal myoprocessor," IEEE Trans. Biomed. Eng., vol. BME-27, no.7, pp.382-395, 1980. C. M. Harris and D. M. Wolpert, "Signal-dependent noise determines motor planning," Nature, vol. 394, no. 6695, pp. 780-784, 1998.
従来の筋電信号モデルでは筋電信号の分散は一定であると仮定しているが、実際には筋電信号には筋力に応じた信号強度依存ノイズが重畳されるため、筋電信号の分散は一定ではなく信号強度依存ノイズに応じてばらつきが生じる。すなわち、従来の筋電信号モデルは信号強度依存ノイズの存在を十分に反映できていない。
一方、筋電信号に重畳された信号強度依存ノイズは筋電信号を解析することで推定することができるが、計測筋電信号のサンプリング周波数は1000〜2000[Hz]であり信号帯域幅が広く情報量が多いため、現実には筋電計において計測筋電信号は5[Hz]程度でダウンサンプリングされて情報量を落とした整流平滑化信号に変換される。そのような筋電計から得られる整流平滑化信号から信号強度依存ノイズを推定したり、計測筋電信号を復元したりすることのニーズがあるが、従来の筋電信号モデルでは整流平滑化信号から信号強度依存ノイズを推定することも計測筋電信号を復元することも困難である。
さらに、特に筋電義手制御の分野において、信号強度依存ノイズを考慮して人工筋電信号を生成することへのニーズがある。筋電義手の制御には、サポートベクターマシンやニューラルネットワークが利用され、事前に識別対象動作の筋電信号のパターンを学習させる必要があるが、筋電義手の使用者の負担を考え、軽い発揮筋力での機械学習が行われる。このため、実生活において大きな筋力が発揮されると、信号強度依存ノイズの影響で識別対象動作の誤識別が起こり、筋電義手が誤作動するおそれがある。そのような誤識別や誤作動を防ぐには、実生活において発揮されると考えられるさまざまな大きさの筋力に対応した筋活動度での筋電信号を、信号強度依存ノイズを考慮した形で人工的に生成し、この人工筋電信号を用いて機械学習することが望まれる。
上記問題に鑑み、本発明は、整流平滑化後の筋電信号から信号強度依存ノイズが重畳された整流平滑化前の筋電信号の分散の分布を推定し、さらに当該推定した分布に基づいて人工筋電信号を生成することを課題とする。
本発明の一局面に従った筋電信号処理方法は、時間窓処理部が、信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算するステップと、分散平均推定部が、整流平滑化信号の平均および筋電信号の整流平滑化に係るフィルタゲインの逆数から筋電信号の分散の平均の推定値を計算するステップと、信号強度依存ノイズ分散推定部が、筋電信号の分散の平均の推定値および筋電信号の分散の分布に係るパラメータから信号強度依存ノイズの分散の推定値を計算するステップとを備えている。
同様に、筋電信号処理装置は、信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算する時間窓処理部と、整流平滑化信号の平均および筋電信号の整流平滑化に係るフィルタゲインの逆数から筋電信号の分散の平均の推定値を計算する分散平均推定部と、筋電信号の分散の平均の推定値および筋電信号の分散の分布に係るパラメータから信号強度依存ノイズの分散の推定値を計算する信号強度依存ノイズ分散推定部とを備えている。
同様に、筋電信号処理プログラムは、信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算する手段、整流平滑化信号の平均および筋電信号の整流平滑化に係るフィルタゲインの逆数から筋電信号の分散の平均の推定値を計算する手段、および筋電信号の分散の平均の推定値および筋電信号の分散の分布に係るパラメータから信号強度依存ノイズの分散の推定値を計算する手段としてコンピュータを機能させるものである。
当該筋電信号処理方法、装置およびプログラムによると、信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号から筋電信号の分散の平均の推定値および信号強度依存ノイズの分散の推定値が計算され、これら推定値を用いて、信号強度依存ノイズが重畳された整流平滑化前の筋電信号の分散の分布を推定することができる。
上記筋電信号処理方法は、第1の乱数発生器が、白色正規乱数を発生させるステップと、シェイピングフィルタが、白色正規乱数の周波数特性を筋電信号と同様の周波数特性に整形するステップと、第2の乱数発生器が、筋電信号の分散の平均の推定値および信号強度依存ノイズの分散の推定値で規定される筋電信号の分散の分布に従った乱数を発生させるステップと、人工筋電信号生成部が、整形された乱数および筋電信号の分散の分布に従った乱数から人工筋電信号を生成するステップとを備えていてもよい。
同様に、上記筋電信号処理装置は、白色正規乱数を発生させる第1の乱数発生器と、白色正規乱数の周波数特性を筋電信号と同様の周波数特性に整形するシェイピングフィルタと、筋電信号の分散の平均の推定値および信号強度依存ノイズの分散の推定値で規定される筋電信号の分散の分布に従った乱数を発生させる第2の乱数発生器と、整形された乱数および筋電信号の分散の分布に従った乱数から人工筋電信号を生成する人工筋電信号生成部とを備えていてもよい。
同様に、上記筋電信号処理プログラムは、白色正規乱数を発生させる手段、白色正規乱数の周波数特性を筋電信号と同様の周波数特性に整形する手段、筋電信号の分散の平均の推定値および信号強度依存ノイズの分散の推定値で規定される筋電信号の分散の分布に従った乱数を発生させる手段、および整形された乱数および筋電信号の分散の分布に従った乱数から人工筋電信号を生成する手段としてコンピュータを機能させてもよい。
当該筋電信号処理方法、装置およびプログラムによると、整形された白色正規乱数および筋電信号の分散の分布に従った乱数から信号強度依存ノイズを考慮した人工筋電信号を生成することができる。
上記筋電信号処理方法は、分散平均推定部が、最大随意筋収縮時の整流平滑化信号、指定された筋活動度および筋電信号の整流平滑化に係るフィルタゲインの逆数から当該筋活動度での筋電信号の分散の平均の推定値を計算するステップと、信号強度依存ノイズ分散推定部が、筋電信号の分散の分布に係るパラメータおよび当該筋活動度での筋電信号の分散の平均の推定値から当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算するステップと、第2の乱数発生器が、当該筋活動度での筋電信号の分散の平均の推定値および当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される当該筋活動度での筋電信号の分散の分布に従った乱数を発生させるステップとを備えていてもよい。
同様に、上記筋電信号処理装置において、分散平均推定部が、最大随意筋収縮時の整流平滑化信号、指定された筋活動度および筋電信号の整流平滑化に係るフィルタゲインの逆数から当該筋活動度での筋電信号の分散の平均の推定値を計算するものであり、信号強度依存ノイズ分散推定部が、筋電信号の分散の分布に係るパラメータおよび当該筋活動度での筋電信号の分散の平均の推定値から当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算するものであり、第2の乱数発生器が、当該筋活動度での筋電信号の分散の平均の推定値および当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される当該筋活動度での筋電信号の分散の分布に従った乱数を発生させるものであってもよい。
同様に、上記筋電信号処理プログラムは、最大随意筋収縮時の整流平滑化信号、指定された筋活動度および筋電信号の整流平滑化に係るフィルタゲインの逆数から当該筋活動度での筋電信号の分散の平均の推定値を計算する手段、筋電信号の分散の分布に係るパラメータおよび当該筋活動度での筋電信号の分散の平均の推定値から当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算する手段、および当該筋活動度での筋電信号の分散の平均の推定値および当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される当該筋活動度での筋電信号の分散の分布に従った乱数を発生させる手段としてコンピュータを機能させてもよい。
当該筋電信号処理方法、装置およびプログラムによると、任意の筋活動度を指定して、その筋活動度での人工筋電信号を生成することができる。
上記筋電信号処理方法、装置およびプログラムにおいて、例えば、筋電信号の分散の分布が逆ガンマ分布であってもよく、パラメータが逆ガンマ分布の形状母数であってもよい。
本発明によると、整流平滑化後の筋電信号から信号強度依存ノイズが重畳された整流平滑化前の筋電信号の分散の分布を推定することができる。さらに、当該推定した分布に基づいて人工筋電信号を生成することができる。
本発明で採用する筋電信号モデルの概略図 本発明の一実施形態に係る筋電信号処理装置を含むシステムの概略図 筋電信号処理装置のブロック図 各サンプリング周波数での筋電信号の分散の平均の推定に係る平均誤差率を示すグラフ 各時間窓での筋電信号の分散の平均の推定に係る平均誤差率を示すグラフ 各サンプリング周波数での信号強度依存ノイズの分散の推定に係る平均誤差率を示すグラフ ある1試行における各負荷での筋電信号の分散の事後分布を示すグラフ 全試行における各筋力での筋電信号の分散の平均および信号強度依存ノイズの分散の各推定値を示すグラフ 各負荷での計測筋電信号とそれを元に生成した人工筋電信号を示すグラフ 各負荷での計測筋電信号と人工筋電信号の平均振幅の平均絶対誤差率を示すグラフ 各負荷での計測筋電信号と人工筋電信号とのパワースペクトル密度の相関係数を示すグラフ 各負荷での筋電信号の分散の分布に係るパラメータの平均絶対誤差率を示すグラフ
以下、適宜図面を参照しながら、実施の形態を詳細に説明する。ただし、必要以上に詳細な説明は省略する場合がある。例えば、既によく知られた事項の詳細説明や実質的に同一の構成に対する重複説明を省略する場合がある。これは、以下の説明が不必要に冗長になるのを避け、当業者の理解を容易にするためである。
なお、発明者らは、当業者が本発明を十分に理解するために添付図面および以下の説明を提供するのであって、これらによって特許請求の範囲に記載の主題を限定することを意図するものではない。
≪筋電信号の分散の事後分布の推定≫
(筋電信号モデル)
図1は、本発明で採用する筋電信号モデルの概略図である。当該筋電信号モデルは、時刻tにおける筋電信号xを、シェイピングフィルタHを通したホワイトガウスノイズWおよび分散σ に基づいて表現する。ここで、σ は時刻tにおける確率変数σの値であり、σの分布は筋力fに応じて決まるσ ̄(本明細書ではσの平均(シグマバー)を便宜上「σ ̄」と表す)および信号強度依存ノイズεによって決定される。さらに、xを整流平滑化した信号yを用いてσの分布を推定可能である。
まず、fとσ ̄との関係は次のように表現できる。
Figure 2017104333
ここで、kとaは定数であり、実験的に求められる。
次に、σはσ ̄とノイズεの和として表される。
Figure 2017104333
ここで、ノイズεは平均0の確率変数である。
これにより、σの平均E[σ]と分散Var[σ]は、
Figure 2017104333
と表される。
はシェイピングフィルタHを通したホワイトガウスノイズWとσの積で表されるため、平均0、分散σ の正規分布に従う。
Figure 2017104333
はxをN次のローパスフィルタを用いて整流平滑化したものである。
Figure 2017104333
ここで、a,a,・・・,aN−1およびb,b,・・・,bはフィルタ係数である。
(分散分布推定)
次に、xからσの分布を推定する方法について説明する。まず、筋電信号xが計測されたとき、分散の事後分布P(σ|x)はベイズの定理を用いると次のように表現できる。
Figure 2017104333
ここで、尤度P(x|σ)は(5)式より、
Figure 2017104333
である。また、事前分布P(σ)にはσ>0を考慮して正規分布と共役の事前分布である逆ガンマ分布IG(α,β)を仮定する。
Figure 2017104333
ただし、α,βは逆ガンマ分布を決定するパラメータである。このとき、(7)式の右辺分子は次のように展開できる。
Figure 2017104333
ここで、κ(α,β)は定数項をまとめたものであり、
Figure 2017104333
となる。(7)式のP(x)は定数、P(σ|x)は確率密度関数でなければならないため、κ(α,β)とP(x)がキャンセルされP(σ|x)はIG(α+1/2,β+x /2)と一致する。このとき、σの事後平均E[σ|x]、事後分散Var[σ|x]は、
Figure 2017104333
となる。このとき、(3)式、(4)式よりE[σ|x]とVar[σ|x]はxが観測された条件下におけるσ ̄とVar[ε]の推定値にそれぞれ一致するため、αをあらかじめ設定することでσ ̄からVar[ε]を求めることができる。
(σの推定)
次に、整流平滑化後の信号yからσ ̄を推定する方法について説明する。yは(6)式で表現できるため、yの期待値は次のようになる。
Figure 2017104333
ここで、xは(5)式に従うためE[|x|]は、
Figure 2017104333
となる。εが時刻tにおけるεの値である。また、E[y]は時刻に依存しないとすると、(14)式は(15)式を用いて、
Figure 2017104333
となる。ここで、(16)式の右辺第2項は期待値0であり、ε(≪σ ̄)のb(>0)による荷重平均を含むため、第1項より十分小さく、無視できる。これにより、σ ̄は(16)式を用いて、
Figure 2017104333
と推定できる。このとき、右辺の分数係数項はフィルタゲインの逆数に相当する。
以上のことから、整流平滑化された筋電信号yと形状母数αを用いて(13)式、(17)式に基づきσ ̄とVar[ε]を推定することができる。
≪人工筋電信号の生成≫
整流平滑化後の信号yから整流平滑化前の信号xの分散σ の事後分布が推定できれば、当該推定した分布を用いて整流平滑化前の信号xと同じ統計学的性質を持つ人工筋電信号を生成することができる。
(筋電信号の再現)
時刻tにおける人工筋電信号zは、シェイピングフィルタHを通した正規乱数w´と分散σ に基づいて表現することができる。
w´は筋電信号xと同様の周波数特性を持つ正規乱数であり、次数Mの自己回帰モデルに基づくシェイピングフィルタHによって生成される。
Figure 2017104333
ここで、wは平均0、分散1の白色正規乱数、v,a(j=1,・・・,M)はそれぞれ自己回帰モデルにおける推定誤差分散とパラメータである。また、自己回帰モデルパラメータの推定において、分散が1となるように正規化した筋電信号xを用いることによって、w´は平均0、分散1の正規乱数となる。
そして、人工筋電信号zは(3)(4)式に基づく分布により生成される乱数σとw´の積で次のように定義される。すなわち、図1の筋電信号モデルにおいて筋電信号xを人工筋電信号zに置換することができる。
Figure 2017104333
以上より、(1)式のパラメータk,aとσの分布、およびシェイピングフィルタHのパラメータをあらかじめ設定することで、発揮筋力fに応じた人工筋電信号zを生成することができる。すなわち、整流平滑化後の信号yから整流平滑化前の信号xを人工筋電信号zとして復元(擬似的に再現)することができる。
(任意の筋活動度での人工筋電信号の生成)
さらに、任意の筋活動度での筋電信号を人工的に生成することもできる。例えば、最大随意筋収縮時の整流平滑化後の信号をymax、筋活動度(%MVC)をrとすると、
=y/ymax (20)
と表される。したがって、yの期待値E[y]は、
E[y]=E[r]・ymax (21)
と表され、(21)式を(17)式に代入すると、
Figure 2017104333
となる。すなわち、最大随意筋収縮時の整流平滑化後の信号ymaxを与え、任意の筋活動度rを指定すれば、(13)式、(22)式に基づいて、筋活動度rでの筋電信号xの分散の事後分布を推定することができる。そして、(19)式に基づいて筋活動度rでの人工筋電信号zを生成することができる。
なお、複数チャネルの筋電信号を計測している場合、被験者のある動作中に着目するチャンネルlおよびその動作に係る主動作筋のチャンネルmの各筋電信号から各チャンネルの筋活動度の時間平均λおよびλを計算し、その比率λ/λを筋活動度rに乗じた値r×λ/λを(22)式中のE[r]として使用してもよい。例えば手を開く動作について考えると、軽く手を開いても思い切り手を開いても各チャネルの筋電信号のパターンは同じようになると考えられるため、上記のように主動作筋の筋活動度の時間平均に対する着目チャンネルの筋活動度の時間平均の比率を乗じて指定された筋活動度を補正することで、動作に対応した特徴的な筋電パターンを崩さずに所望の筋活動度での人工筋電信号を生成することができる。
≪筋電信号処理装置≫
次に、本発明に係る筋電信号処理装置について説明する。図2は、本発明の一実施形態に係る筋電信号処理装置を含むシステムの概略図である。本システム100は、筋電計10と筋電信号処理装置20とで構成される。
筋電計10は、任意の個数(本実施形態では4個)の筋電センサ11とデータ収集装置12とを備えている。データ収集装置12と筋電センサ11とはパラレルケーブルまたはシリアルケーブルで接続されている。さらに、これらケーブルは、データ収集装置12から取り外しできるようになっている。
筋電センサ11は、被験者の任意の身体部位の表面筋電位信号(筋電信号)を計測するためのものである。筋電信号は筋肉の複合活動電位を波形化して表したものであり、複合活動電位とは筋肉を構成する各筋線維から発生し、容積伝導により体表面に伝搬する活動電位を合計したものをいう。複合活動電位は筋収縮の程度により変化する。例えば、筋肉が小さな力を出している場合には収縮する筋線維の数は少なく、筋線維から発生する活動電位の数も少ない。そこから力を増加させていくと収縮する筋線維の数が徐々に増加し、活動電位の数も増加する。これにより、筋電信号として観測される複合活動電位は増大する。このように、力の出し方と筋電信号とは相関関係にある。
表面信号は2つの表面電極間の電位差として捕捉される。具体的には、各筋電センサ11は、2個の表面電極101(プラス電極とマイナス電極)を有している。これら電極101を対象部位の筋繊維の走行方向に沿って数cm程度の間隔で被験者の体表面に貼付する。なお、いずれか一つの筋電センサ11は、図略のリファレンス電極をさらに有する。
データ収集装置12は、図略のアンプ、フィルタ、AD変換器、CPU、通信デバイスなどを備えており、筋電センサ11が計測した筋電信号(筋電信号モデルにおけるx)を収集し、当該収集した信号を整流および平滑化して整流平滑化信号(筋電信号モデルにおけるy)を生成して筋電信号処理装置20へ送信する。データ収集装置12と筋電信号処理装置20との間の通信には、Wi−Fi(登録商標)、Bluetooth(登録商標)、ZigBee(登録商標)、特定小電力無線などの無線通信を使用することができる。このようにデータ収集装置12において計測筋電信号を整流平滑化してデータ量を削減することで、比較的低速の通信回線で多くのチャネルの筋電信号を筋電信号処理装置20へ送信することができる。
筋電信号処理装置20は、図1の筋電信号モデルに基づく上記の筋電信号処理を行う装置である。すなわち、筋電信号処理装置20は、筋電計10から整流平滑化信号yを受信して上述の筋電信号処理を行い、筋電センサ11が計測した筋電信号xの分散の事後分布を推定し、さらに当該推定した分布に基づいて人工筋電信号を生成する。筋電信号処理装置20が生成した人工筋電信号は、筋電信号処理装置20または別の装置のモニター画面にリアルタイムで表示したり、筋電義手制御の機械学習用の信号として別の装置に出力したりすることができる。
筋電信号処理装置20として、汎用コンピュータ(デスクトップコンピュータやノートコンピュータなど)、スマートフォン、タブレット端末などを使用することができる。これらコンピュータや端末に、上述の筋電信号処理を実現するためのプログラムまたはアプリケーションをインストールすることで、当該コンピュータや端末を筋電信号処理装置20として利用することができる。
図3は、筋電信号処理装置20のブロック図である。筋電信号処理装置20は、時間窓処理部21、分散平均推定部22、信号強度依存ノイズ分散推定部23、ホワイトガウスノイズ発生器24、シェイピングフィルタ25、逆ガンマ分布乱数発生器26および人工筋電信号生成部27を備えている。このうち、整流平滑化信号yから筋電センサ11が計測した筋電信号xの分散の事後分布を推定するのに必要なブロックは時間窓処理部21、分散平均推定部22および信号強度依存ノイズ分散推定部23の3つであり、残り4つのブロックを追加することで人工筋電信号の生成が可能となる。
時間窓処理部21は、所定の時間長(時間窓)ごとに整流平滑化信号yの平均E[y]を計算する。ここで、yは信号強度依存ノイズεが重畳された筋電信号xを整流平滑化したものであり、時間窓処理部21は、実質的に(14)式の演算を実行する。
分散平均推定部22は、整流平滑化信号の平均E[y]および筋電信号xの整流平滑化に係るフィルタゲインの逆数から筋電信号の分散の平均σ ̄の推定値を計算する。すなわち、分散平均推定部22は、実質的に(17)式の演算を実行する。
また、任意の筋活動度rでの人工筋電信号zを生成する場合、分散平均推定部22は、最大随意筋収縮時の整流平滑化信号ymax、指定された筋活動度rおよび筋電信号の整流平滑化に係るフィルタゲインの逆数から筋活動度rでの筋電信号の分散の平均σ ̄の推定値を計算する。この場合、分散平均推定部22は、実質的に(22)式の演算を実行する。なお、分散平均推定部22は、指定された筋活動度rをそのまま使用してもよいし、上述したように、主動作筋の筋活動度の時間平均に対する着目チャンネルの筋活動度の時間平均の比率を乗じて補正した値を使用してもよい。
最大随意筋収縮時の筋電信号を計測するには被験者に筋活動度100(100%MVC)で筋力を発揮してもらう必要があるが、被験者に負担がかかるなどの理由でそれができない場合がある。そのような場合には最大随意筋収縮時の整流平滑化信号を取得できないため、分散平均推定部22に入力する最大随意筋収縮時の整流平滑化信号ymaxとして、例えば、他の被験者の実測値などで代用した想定値を用いることができる。もちろん、被験者の最大随意筋収縮時の筋電信号が計測できる場合にはymaxとして想定値ではなく実測値を用いてもよいことは言うまでもない。
信号強度依存ノイズ分散推定部23は、筋電信号の分散の平均σ ̄の推定値および筋電信号の分散の分布に係るパラメータから信号強度依存ノイズの分散Var[ε]の推定値を計算する。ここで、筋電信号の分散の分布として逆ガンマ分布を想定すると、当該パラメータとして逆ガンマ分布の形状母数αを採用することができる。すなわち、信号強度依存ノイズ分散推定部23は、実質的に(13)式の演算を実行する。
また、任意の筋活動度rでの人工筋電信号zを生成する場合、信号強度依存ノイズ分散推定部23は、筋電信号の分散の分布に係るパラメータおよび筋活動度rでの筋電信号の分散の平均σ ̄の推定値から筋活動度rでの筋電信号に重畳される信号強度依存ノイズの分散Var[ε]の推定値を計算する。
ホワイトガウスノイズ発生器24は、白色正規乱数wを発生させる。ここで、wの平均は0、分散は1である。
シェイピングフィルタ25は、白色正規乱数wの周波数特性を筋電信号xと同様の周波数特性に整形して乱数w´を生成する。ここで、w´の平均は0、分散は1である。すなわち、シェイピングフィルタ25は、実質的に(18)式の演算を実行する。
逆ガンマ分布乱数発生器26は、筋電信号の分散の平均σ ̄の推定値および信号強度依存ノイズの分散Var[ε]の推定値で規定される筋電信号の分散の分布に従った乱数σを発生させる。ここで、筋電信号の分散の分布として逆ガンマ分布を想定しているため、当該乱数σとして逆ガンマ分布乱数を使用する。
また、任意の筋活動度rでの人工筋電信号zを生成する場合、逆ガンマ分布乱数発生器26は、筋活動度rでの筋電信号の分散の平均σ ̄の推定値および筋活動度rでの筋電信号に重畳される信号強度依存ノイズの分散Var[ε]の推定値で規定される筋活動度rでの筋電信号の分散の分布に従った乱数σを発生させる。
なお、逆ガンマ分布を規定するパラメータである形状母数αおよび尺度母数βは、例えば、第2種の最尤推定に基づいて、あらかじめ得られた筋電信号系列について周辺尤度を最大化することで求めることができる。
人工筋電信号生成部27は、シェイピングフィルタ25で整形後の乱数w´および筋電信号の分散の分布に従った乱数σから人工筋電信号zを生成する。すなわち、人工筋電信号生成部27は、実質的に(19)式の演算を実行する。
≪実証実験≫
(分散分布推定精度)
筋電信号処理装置20による分散分布推定精度を確認するため人工的に生成した人工筋電信号を用いて実験を行った。人工筋電信号の生成は以下の手順で行うことで、分散にノイズが重畳した信号を擬似的に生成した。
(1)長さTの離散系列{σ ;t=1,・・・,T}を逆ガンマ分布IG(α,β)に従う乱数で発生させる。
(2)各tにおいて正規分布N(0,σ )に従う乱数を1サンプルずつ発生させ、これを系列{x}とする。
(3){x}をサンプリング周波数F[Hz]で計測された筋電信号とみなす。
このとき、{x}の分散の平均、および分散の分散の真値はそれぞれ以下の値となる。
Figure 2017104333
分散分布推定精度検証は、これら真値と筋電信号処理装置20で推定したσ ̄およびVar[ε]とを比較することによって行った。
σ ̄の推定は、{x}を整流平滑化した信号を{y}とし、{y}をFdown[Hz]でダウンサンプリングした信号のうち最初のL[ms]を用いて行った。このとき、Lを1000[ms]に固定してFdownを500[Hz]、200[Hz]、100[Hz]、50[Hz]、20[Hz]、10[Hz]と変更した場合、およびFdownを1000[Hz]に固定してLを1000[ms]、500[ms]、200[ms]、100[ms]、50[ms]、20[ms]、10[ms]と変更した場合の各条件において、真値σ ̄ とσ ̄との%誤差を算出した。比較手法として、同様の条件において{x}から算出した分散の最尤推定値を用いた。
Var[ε]の推定精度検証では、Lを1000[ms]に固定してFdownを1000[Hz]、50[Hz]と変更した各条件において、Var[ε]との%誤差を算出した。また、最尤推定に基づく分散の分散との比較を行った。ここで、最尤推定に基づく分散の分散とは、L[ms]の{x}をK分割し、各分割において算出したK個の分散最尤推定値の分散のことをいう。なお、事前分布パラメータは経験ベイズ法を用いて決定した。
各実験において、T=1000、F=1000とし、平滑化処理にはカットオフ周波数1[Hz]の2次のバターワースフィルタを使用した。また、%誤差は真値を20通り(α=15、β=0.5,1,1.5,・・・,10)変更した際の平均値である。
図4は、各サンプリング周波数での筋電信号の分散の平均σ ̄の推定に係る平均誤差率を示すグラフである。図5は、各時間窓での筋電信号の分散の平均σ ̄の推定に係る平均誤差率を示すグラフである。図4のグラフからはFdown≦200で、図5のグラフからはL≦100で、本発明(図1の筋電信号モデルに基づく推定)と従来技術(最尤推定)との間に有意差が認められた。図6は、各サンプリング周波数での信号強度依存ノイズの分散Var[ε]の推定に係る平均誤差率を示すグラフである。Bonferroni法に基づく多重比較検定を行った結果、本発明と従来技術との間に有意差が確認できた。
図4のグラフを参照すると、Fdownが高い場合には本発明および従来技術ともに誤差率5[%]以下であり推定精度が高いが、Fdownが下がるにつれ従来技術の精度は悪化し、本発明は高精度を維持する。また、図5のグラフを参照すると、窓幅Lを短くすると本発明および従来技術ともに精度が悪化するが、L=10[ms]においても本発明の推定誤差は5[%]程度であり、従来技術よりも有意に小さい。このことから、本発明によると、より少ないサンプルでσ ̄を推定できると言える。
また、図6のグラフを参照すると、本発明は従来技術と比較して、Var[ε]の推定誤差が有意に小さい。これは、従来技術の最尤法が各分割において分散一定を仮定するため、各時刻独立に重畳するノイズを十分に表現できないことが原因である。一方、本発明では各時刻独立の分散値を仮定し、それらの母集団の事後分布を推定することで、分散の分散を誤差率10[%]以下の良好な精度で推定することができる。
以上のように、本発明の一実施形態に係る筋電信号処理装置20によると、信号強度依存ノイズが重畳された筋電信号の分散の分布を高精度に推定することができる。
次に、実データに対する本発明の有効性を確認するために、筋電信号を用いて分散事後分布P(σ|x)の推定実験を行った。筋電信号の計測では、被験者5名(平均年齢:23.3±0.8、右利き)の上腕二頭筋上の皮膚表面にAg/AgCl電極を貼付し、座位で上腕を下垂し手掌を上に向けて前腕を水平前方に曲げた姿勢をとらせ、肘を机につき手首に負荷をのせ姿勢を10秒間維持させた際の筋電信号をサンプリング周波数1000[Hz]で取得した。このとき、負荷を500[g]、1000[g]、1500[g]、2000[g]と変更し、それぞれ5試行計測した。なお、計測には日本光電工業株式会社製マルチテレメータシステムWEB−5000(高域遮断周波数:100[Hz]、低域遮断周波数:5.3[Hz])を用いた。
筋電信号の解析には、計測した10秒間のうち後半の5秒間のデータを用いた。このとき、前腕および手の自重を体重比0.022、質量中心を肘関節からの長さ比0.318、上腕二頭筋の肘関節回転中心に対するモーメントアーム長を0.03[m]と仮定し、上腕二頭筋の発揮筋力を計算した。そして、負荷の増加に伴う分散事後分布P(σ|x)の変化と、筋力の変化に伴うE[σ|x]およびVar[σ|x]の変化を求めた。
なお、事前分布パラメータはあらかじめ計測した2000[g]の負荷時の筋電信号1試行分を用いて、経験ベイズ法に基づき被験者ごとに設定した。
図7は、ある1試行における各負荷での筋電信号の分散の事後分布を示すグラフである。図8は、全試行における各筋力での筋電信号の分散の平均および信号強度依存ノイズの分散の各推定値を示すグラフである。図8にはSteel-Dwass法に基づく検定結果を併記している。
図7のグラフを参照すると、負荷が大きくなるにつれてP(σ|x)が右に移動し広がりが大きくなっている。これは、負荷の増加に伴い、事後平均と事後分散が増加傾向にあることを意味する。また、図8のグラフを参照すると、筋力と事後平均E[σ|x]が単調増加傾向にあることから、(1)式で筋力が表現可能であることがわかる。このとき、(1)式中の指数aを求めると、被験者順にa=1.18,0.56,0.81,0.74,0.65となった。一方、非特許文献1と同様に最尤推定に基づく求めた分散σを使用して指数aを求めると、a=1.17,0.55,0.78,0.73,0.65となった。このことから、本発明で算出できる事後平均E[σ|x]は非特許文献1の方法において最尤法で求められる分散と同様の役割を果たしていることがわかる。
一方、筋力に応じた事後分散Var[σ|x]の増加は事後平均E[σ|x]の増加と(13)式により明らかであるが、特筆すべきは、本発明において事後分散Var[σ|x]が事後平均E[σ|x]の2乗に比例する点である。過去に、等尺性収縮時における筋力の分散が平均筋力のおよそ2乗に比例するという実験結果が報告されているが、本発明は、筋電信号の分散について同様の関係が成立することを示している。このことは、信号強度依存ノイズにより筋電信号の分散値に対してもノイズが重畳することを意味している。
以上のように、本発明の一実施形態に係る筋電信号処理装置20によると、筋電信号の分散分布を推定でき、さらに分散に重畳するノイズを推定することができる。
(人工筋電信号の生成実験)
次に、筋電信号処理装置20による人工筋電信号の生成実験を行った。実験では、まず健常大学生1名から筋電信号を計測し、分散分布パラメータσ ̄とVar[ε]、およびシェイピングフィルタ25のパラメータを推定した。そして、推定したパラメータに基づき生成した人工筋電信号と計測した筋電信号との比較を行った。
筋電信号の計測条件は上記と同じである。負荷を500[g]、1000[g]、1500[g]、2000[g]と変更し、それぞれ1試行計測した。パラメータ推定および比較には計測した10秒間のうち後半の5秒間のデータを使用した。
分散事後分布の推定は負荷ごとに行った。ただし、事前パラメータはあらかじめ計測した負荷2000[g]時の筋電信号を用いて第2種の最尤推定に基づき設定し、実験を通じて共通とした。また、シェイピングフィルタ25のパラメータは負荷2000[g]時の筋電信号をBurg法によって決定し、モデル次数はベイズ情報量規準に基づきM=20とした。人工筋電信号の生成は負荷ごとに5回ずつ行い、逆ガンマ分布乱数σの生成にはTanizaki法を用いた。
比較は、平均振幅、周波数成分、分散分布の3項目について行った。平均振幅は信号を整流平滑化した値の平均値とし、計測筋電信号における値を真値としたときの平均絶対誤差率を各負荷において算出した。周波数成分の比較では、計測筋電信号と人工筋電信号のパワースペクトル密度をそれぞれ求め、それらの間の相関関係を算出した。分散分布の比較については、筋電信号の分散の分布を逆ガンマ分布と仮定することでパラメータα、βを第2種の最尤推定によって推定し、平均振幅の場合と同様に平均絶対誤差率を算出した。
図9は、各負荷での計測筋電信号とそれを元に生成した人工筋電信号を示すグラフである。図9のグラフを参照すると、負荷が大きくなるにつれて計測筋電信号と人工筋電信号の振幅がともに大きくなることがわかる。
図10は、各負荷での計測筋電信号と人工筋電信号の平均振幅の平均絶対誤差率を示すグラフである。図10のグラフを参照すると、人工筋電信号はどの負荷においても5[%]程度の平均絶対誤差率で筋電信号の振幅を再現できていることが確認できる。
図11は、各負荷での計測筋電信号と人工筋電信号とのパワースペクトル密度の相関係数を示すグラフである。図11のグラフを参照すると、すべての負荷において相関係数は0.8以上であり、周波数成分においても人工筋電信号が高い精度で筋電信号を再現できていることがわかる。
図12は、各負荷での筋電信号の分散の分布に係るパラメータの平均絶対誤差率を示すグラフである。図12のグラフを参照すると、形状母数αは負荷が大きくなるにつれて平均誤差率が大きくなる傾向があるものの0.01[%]未満と非常に小さい。一方、尺度母数βの平均誤差率はαの平均誤差率よりも大きい5−10[%]程度である。この原因として、(13)式および(17)式においてαは事後分布の導出に直接用いられているが、βはそれに用いられていないため、事後分布の尺度に関して推定誤差が大きくなったことが考えられる。しかし、形状母数αは分布の本質的な形を決定するのに対して尺度母数βは分布の広がりを決定するものであるため、βが人工筋電信号の生成に及ぼす影響はαに比べて小さく問題にならないと考える。
以上のように、本発明の一実施形態に係る筋電信号処理装置20によると、計測筋電信号の振幅、周波数、分散分布を再現する人工筋電信号を生成することができる。
以上のように、本発明における技術の例示として、実施の形態を説明した。そのために、添付図面および詳細な説明を提供した。
したがって、添付図面および詳細な説明に記載された構成要素の中には、課題解決のために必須な構成要素だけでなく、上記技術を例示するために、課題解決のためには必須でない構成要素も含まれ得る。そのため、それらの必須ではない構成要素が添付図面や詳細な説明に記載されていることをもって、直ちに、それらの必須ではない構成要素が必須であるとの認定をするべきではない。
また、上述の実施の形態は、本発明における技術を例示するためのものであるから、特許請求の範囲またはその均等の範囲において種々の変更、置き換え、付加、省略などを行うことができる。例えば、上記では筋電信号の分散の分布として逆ガンマ分布を想定したが、筋電信号の分散の分布として逆ガンマ分布以外に、歪みを有する非対称の分布を適用してもよい。
20 筋電信号処理装置
21 時間窓処理部
22 分散平均推定部
23 信号強度依存ノイズ分散推定部
24 ホワイトガウスノイズ発生器(第1の乱数発生器)
25 シェイピングフィルタ
26 逆ガンマ分布乱数発生器(第2の乱数発生器)
27 人工筋電信号生成部

Claims (12)

  1. 時間窓処理部が、信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算するステップと、
    分散平均推定部が、前記整流平滑化信号の平均および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋電信号の分散の平均の推定値を計算するステップと、
    信号強度依存ノイズ分散推定部が、前記筋電信号の分散の平均の推定値および前記筋電信号の分散の分布に係るパラメータから前記信号強度依存ノイズの分散の推定値を計算するステップとを備えた筋電信号処理方法。
  2. 第1の乱数発生器が、白色正規乱数を発生させるステップと、
    シェイピングフィルタが、前記白色正規乱数の周波数特性を前記筋電信号と同様の周波数特性に整形するステップと、
    第2の乱数発生器が、前記筋電信号の分散の平均の推定値および前記信号強度依存ノイズの分散の推定値で規定される前記筋電信号の分散の分布に従った乱数を発生させるステップと、
    人工筋電信号生成部が、前記整形された乱数および前記筋電信号の分散の分布に従った乱数から人工筋電信号を生成するステップとを備えた請求項1に記載の筋電信号処理方法。
  3. 前記分散平均推定部が、最大随意筋収縮時の前記整流平滑化信号、指定された筋活動度および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋活動度での筋電信号の分散の平均の推定値を計算するステップと、
    前記信号強度依存ノイズ分散推定部が、前記筋電信号の分散の分布に係るパラメータおよび前記筋活動度での筋電信号の分散の平均の推定値から前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算するステップと、
    前記第2の乱数発生器が、前記筋活動度での筋電信号の分散の平均の推定値および前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される前記筋活動度での筋電信号の分散の分布に従った乱数を発生させるステップとを備えた請求項2に記載の筋電信号処理方法。
  4. 前記筋電信号の分散の分布が逆ガンマ分布であり、
    前記パラメータが逆ガンマ分布の形状母数である、請求項1ないし請求項3のいずれかに記載の筋電信号処理方法。
  5. 信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算する時間窓処理部と、
    前記整流平滑化信号の平均および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋電信号の分散の平均の推定値を計算する分散平均推定部と、
    前記筋電信号の分散の平均の推定値および前記筋電信号の分散の分布に係るパラメータから前記信号強度依存ノイズの分散の推定値を計算する信号強度依存ノイズ分散推定部とを備えた筋電信号処理装置。
  6. 白色正規乱数を発生させる第1の乱数発生器と、
    前記白色正規乱数の周波数特性を前記筋電信号と同様の周波数特性に整形するシェイピングフィルタと、
    前記筋電信号の分散の平均の推定値および前記信号強度依存ノイズの分散の推定値で規定される前記筋電信号の分散の分布に従った乱数を発生させる第2の乱数発生器と、
    前記整形された乱数および前記筋電信号の分散の分布に従った乱数から人工筋電信号を生成する人工筋電信号生成部とを備えた請求項5に記載の筋電信号処理装置。
  7. 前記分散平均推定部が、最大随意筋収縮時の前記整流平滑化信号、指定された筋活動度および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋活動度での筋電信号の分散の平均の推定値を計算するものであり、
    前記信号強度依存ノイズ分散推定部が、前記筋電信号の分散の分布に係るパラメータおよび前記筋活動度での筋電信号の分散の平均の推定値から前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算するものであり、
    前記第2の乱数発生器が、前記筋活動度での筋電信号の分散の平均の推定値および前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される前記筋活動度での筋電信号の分散の分布に従った乱数を発生させるものである、請求項6に記載の筋電信号処理装置。
  8. 前記筋電信号の分散の分布が逆ガンマ分布であり、
    前記パラメータが逆ガンマ分布の形状母数である、請求項5ないし請求項7のいずれかに記載の筋電信号処理装置。
  9. 信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算する手段、
    前記整流平滑化信号の平均および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋電信号の分散の平均の推定値を計算する手段、および
    前記筋電信号の分散の平均の推定値および前記筋電信号の分散の分布に係るパラメータから前記信号強度依存ノイズの分散の推定値を計算する手段としてコンピュータを機能させる筋電信号処理プログラム。
  10. 白色正規乱数を発生させる手段、
    前記白色正規乱数の周波数特性を前記筋電信号と同様の周波数特性に整形する手段、
    前記筋電信号の分散の平均の推定値および前記信号強度依存ノイズの分散の推定値で規定される前記筋電信号の分散の分布に従った乱数を発生させる手段、および
    前記整形された乱数および前記筋電信号の分散の分布に従った乱数から人工筋電信号を生成する手段としてコンピュータを機能させる請求項9に記載の筋電信号処理プログラム。
  11. 最大随意筋収縮時の前記整流平滑化信号、指定された筋活動度および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋活動度での筋電信号の分散の平均の推定値を計算する手段、
    前記筋電信号の分散の分布に係るパラメータおよび前記筋活動度での筋電信号の分散の平均の推定値から前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算する手段、および
    前記筋活動度での筋電信号の分散の平均の推定値および前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される前記筋活動度での筋電信号の分散の分布に従った乱数を発生させる手段としてコンピュータを機能させる請求項10に記載の筋電信号処理プログラム。
  12. 前記筋電信号の分散の分布が逆ガンマ分布であり、
    前記パラメータが逆ガンマ分布の形状母数である、請求項9ないし請求項11のいずれかに記載の筋電信号処理プログラム。
JP2015241474A 2015-12-10 2015-12-10 筋電信号処理方法、装置およびプログラム Active JP6652252B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2015241474A JP6652252B2 (ja) 2015-12-10 2015-12-10 筋電信号処理方法、装置およびプログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2015241474A JP6652252B2 (ja) 2015-12-10 2015-12-10 筋電信号処理方法、装置およびプログラム

Publications (2)

Publication Number Publication Date
JP2017104333A true JP2017104333A (ja) 2017-06-15
JP6652252B2 JP6652252B2 (ja) 2020-02-19

Family

ID=59058060

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2015241474A Active JP6652252B2 (ja) 2015-12-10 2015-12-10 筋電信号処理方法、装置およびプログラム

Country Status (1)

Country Link
JP (1) JP6652252B2 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109965847A (zh) * 2019-04-08 2019-07-05 清华大学 服务器和信号分析系统
JP2021129896A (ja) * 2020-02-21 2021-09-09 国立大学法人広島大学 交感神経活動推定装置、交感神経活動推定方法及びプログラム
CN114469142A (zh) * 2022-01-06 2022-05-13 中南大学 一种基于人体肌肉动力学模型和肌电信号的肌力解码方法
CN115034273A (zh) * 2021-12-27 2022-09-09 驻马店市中心医院 一种基于模式识别的肌电生物反馈设备及系统

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109965847A (zh) * 2019-04-08 2019-07-05 清华大学 服务器和信号分析系统
CN109965847B (zh) * 2019-04-08 2023-11-07 清华大学 服务器和信号分析系统
JP2021129896A (ja) * 2020-02-21 2021-09-09 国立大学法人広島大学 交感神経活動推定装置、交感神経活動推定方法及びプログラム
JP7390716B2 (ja) 2020-02-21 2023-12-04 国立大学法人広島大学 交感神経活動推定装置、交感神経活動推定方法及びプログラム
CN115034273A (zh) * 2021-12-27 2022-09-09 驻马店市中心医院 一种基于模式识别的肌电生物反馈设备及系统
CN115034273B (zh) * 2021-12-27 2023-09-01 驻马店市中心医院 一种基于模式识别的肌电生物反馈设备及系统
CN114469142A (zh) * 2022-01-06 2022-05-13 中南大学 一种基于人体肌肉动力学模型和肌电信号的肌力解码方法

Also Published As

Publication number Publication date
JP6652252B2 (ja) 2020-02-19

Similar Documents

Publication Publication Date Title
KR101666399B1 (ko) 다중 채널 표면 근전도에서의 신체 관절 운동학 정보 추출 방법, 이를 수행하기 위한 기록 매체 및 장치
CN110801226A (zh) 一种基于表面肌电信号的人体膝关节力矩测试系统方法及应用
JP5464072B2 (ja) 筋活動診断装置および方法、並びにプログラム
US11317840B2 (en) Method for real time analyzing stress using deep neural network algorithm
JP6652252B2 (ja) 筋電信号処理方法、装置およびプログラム
US8798726B2 (en) Method and apparatus for eliminating motion artifacts of bio signal using personalized bio signal pattern
Mobasser et al. Estimation of elbow-induced wrist force with EMG signals using fast orthogonal search
He et al. Single channel blind source separation on the instantaneous mixed signal of multiple dynamic sources
Triwiyanto et al. Evaluating the performance of Kalman filter on elbow joint angle prediction based on electromyography
US20160150987A1 (en) Biosignal processing apparatus and method
CN106901732B (zh) 突变状态下肌力与肌张力的测量方法和测量装置
Menegaldo et al. The influence of modeling hypothesis and experimental methodologies in the accuracy of muscle force estimation using EMG-driven models
Zhang et al. Human joint motion estimation for electromyography (EMG)-based dynamic motion control
CN113116321A (zh) 基于pso-grnn神经网络的无创连续血压测量系统
Kim et al. Missing sample recovery for wireless inertial sensor-based human movement acquisition
Bai et al. Novel time-frequency approach for muscle fatigue detection based on sEMG
KR20090001146A (ko) 가중 퍼지 소속함수 기반 신경망 알고리즘을 이용한생체신호의 이상데이터 추출 장치 및 그 방법
Johns et al. Force modelling of upper limb biomechanics using ensemble fast orthogonal search on high-density electromyography
Dong et al. Development of a fatigue-tracking system for monitoring human body movement
Yahya et al. Electromyography signal on biceps muscle in time domain analysis
Segala et al. Nonlinear smooth orthogonal decomposition of kinematic features of sawing reconstructs muscle fatigue evolution as indicated by electromyography
Mountjoy et al. Use of the fast orthogonal search method to estimate optimal joint angle for upper limb Hill-muscle models
Xiao et al. Does force myography recorded at the wrist correlate to resistance load levels during bicep curls?
JP5913493B2 (ja) 欠損生体信号推定方法
Phinyomark et al. Optimal EMG amplitude detectors for muscle-computer interface

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20181127

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20191016

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20191029

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20191219

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20200107

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20200116

R150 Certificate of patent or registration of utility model

Ref document number: 6652252

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250