(57)【特許請求の範囲】
【請求項１】
時間窓処理部が、信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算するステップと、
分散平均推定部が、前記整流平滑化信号の平均および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋電信号の分散の平均の推定値を計算するステップと、
信号強度依存ノイズ分散推定部が、前記筋電信号の分散の平均の推定値および前記筋電信号の分散の分布を規定するパラメータから前記信号強度依存ノイズの分散の推定値を計算するステップとを備え、
前記筋電信号の分散の分布は、前記筋電信号の分布に対して共役な性質を有するとともに、分散の非負性を満たす分布である
筋電信号処理方法。
【請求項２】
第１の乱数発生器が、白色正規乱数を発生させるステップと、
シェイピングフィルタが、前記白色正規乱数の周波数特性を前記筋電信号と同様の周波数特性に整形するステップと、
第２の乱数発生器が、前記筋電信号の分散の平均の推定値および前記信号強度依存ノイズの分散の推定値で規定される前記筋電信号の分散の分布に従った乱数を発生させるステップと、
人工筋電信号生成部が、前記整形された乱数および前記筋電信号の分散の分布に従った乱数から人工筋電信号を生成するステップとを備えた請求項１に記載の筋電信号処理方法。
【請求項３】
前記分散平均推定部が、最大随意筋収縮時の前記整流平滑化信号、指定された筋活動度および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋活動度での筋電信号の分散の平均の推定値を計算するステップと、
前記信号強度依存ノイズ分散推定部が、前記筋電信号の分散の分布を規定するパラメータおよび前記筋活動度での筋電信号の分散の平均の推定値から前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算するステップと、
前記第２の乱数発生器が、前記筋活動度での筋電信号の分散の平均の推定値および前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される前記筋活動度での筋電信号の分散の分布に従った乱数を発生させるステップとを備えた請求項２に記載の筋電信号処理方法。
【請求項４】
前記筋電信号の分散の分布が逆ガンマ分布であり、
前記パラメータが逆ガンマ分布の形状母数である、請求項１ないし請求項３のいずれかに記載の筋電信号処理方法。
【請求項５】
信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算する時間窓処理部と、
前記整流平滑化信号の平均および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋電信号の分散の平均の推定値を計算する分散平均推定部と、
前記筋電信号の分散の平均の推定値および前記筋電信号の分散の分布を規定するパラメータから前記信号強度依存ノイズの分散の推定値を計算する信号強度依存ノイズ分散推定部とを備え、
前記筋電信号の分散の分布は、前記筋電信号の分布に対して共役な性質を有するとともに、分散の非負性を満たす分布である
筋電信号処理装置。
【請求項６】
白色正規乱数を発生させる第１の乱数発生器と、
前記白色正規乱数の周波数特性を前記筋電信号と同様の周波数特性に整形するシェイピングフィルタと、
前記筋電信号の分散の平均の推定値および前記信号強度依存ノイズの分散の推定値で規定される前記筋電信号の分散の分布に従った乱数を発生させる第２の乱数発生器と、
前記整形された乱数および前記筋電信号の分散の分布に従った乱数から人工筋電信号を生成する人工筋電信号生成部とを備えた請求項５に記載の筋電信号処理装置。
【請求項７】
前記分散平均推定部が、最大随意筋収縮時の前記整流平滑化信号、指定された筋活動度および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋活動度での筋電信号の分散の平均の推定値を計算するものであり、
前記信号強度依存ノイズ分散推定部が、前記筋電信号の分散の分布を規定するパラメータおよび前記筋活動度での筋電信号の分散の平均の推定値から前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算するものであり、
前記第２の乱数発生器が、前記筋活動度での筋電信号の分散の平均の推定値および前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される前記筋活動度での筋電信号の分散の分布に従った乱数を発生させるものである、請求項６に記載の筋電信号処理装置。
【請求項８】
前記筋電信号の分散の分布が逆ガンマ分布であり、
前記パラメータが逆ガンマ分布の形状母数である、請求項５ないし請求項７のいずれかに記載の筋電信号処理装置。
【請求項９】
信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算する手段、
前記整流平滑化信号の平均および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋電信号の分散の平均の推定値を計算する手段、および
前記筋電信号の分散の平均の推定値および前記筋電信号の分散の分布を規定するパラメータから前記信号強度依存ノイズの分散の推定値を計算する手段としてコンピュータを機能させ、
前記筋電信号の分散の分布は、前記筋電信号の分布に対して共役な性質を有するとともに、分散の非負性を満たす分布である
筋電信号処理プログラム。
【請求項１０】
白色正規乱数を発生させる手段、
前記白色正規乱数の周波数特性を前記筋電信号と同様の周波数特性に整形する手段、
前記筋電信号の分散の平均の推定値および前記信号強度依存ノイズの分散の推定値で規定される前記筋電信号の分散の分布に従った乱数を発生させる手段、および
前記整形された乱数および前記筋電信号の分散の分布に従った乱数から人工筋電信号を生成する手段としてコンピュータを機能させる請求項９に記載の筋電信号処理プログラム。
【請求項１１】
最大随意筋収縮時の前記整流平滑化信号、指定された筋活動度および前記筋電信号の整流平滑化に係るフィルタゲインの逆数から前記筋活動度での筋電信号の分散の平均の推定値を計算する手段、
前記筋電信号の分散の分布を規定するパラメータおよび前記筋活動度での筋電信号の分散の平均の推定値から前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算する手段、および
前記筋活動度での筋電信号の分散の平均の推定値および前記筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される前記筋活動度での筋電信号の分散の分布に従った乱数を発生させる手段としてコンピュータを機能させる請求項１０に記載の筋電信号処理プログラム。
【請求項１２】
前記筋電信号の分散の分布が逆ガンマ分布であり、
前記パラメータが逆ガンマ分布の形状母数である、請求項９ないし請求項１１のいずれかに記載の筋電信号処理プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、表面筋電位信号（以下、筋電信号という）の処理技術に関し、特に、整流平滑化後の筋電信号から信号強度依存ノイズが重畳された整流平滑化前の筋電信号の分散の分布を推定し、さらに当該推定した分布に基づいて人工筋電信号を生成する技術に関する。
【背景技術】
【０００２】
筋電信号は筋活動の定量評価に有用であり、従来、筋電義手の制御や、リハビリ時やトレーニング時の運動解析に応用されてきた。それらの応用においては、筋電信号から特徴量を抽出すべく、筋電信号がモデル化され用いられる（例えば、非特許文献１参照）。
【０００３】
近年、筋電信号を含む生体信号の研究が進むにつれ、信号強度依存ノイズの存在が明らかになってきた（例えば、非特許文献２参照）。すなわち、腕に思いきり力を入れると腕が震えるように、筋電信号には筋力の強さに応じた信号強度依存ノイズが重畳される。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】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.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
従来の筋電信号モデルでは筋電信号の分散は一定であると仮定しているが、実際には筋電信号には筋力に応じた信号強度依存ノイズが重畳されるため、筋電信号の分散は一定ではなく信号強度依存ノイズに応じてばらつきが生じる。すなわち、従来の筋電信号モデルは信号強度依存ノイズの存在を十分に反映できていない。
【０００６】
一方、筋電信号に重畳された信号強度依存ノイズは筋電信号を解析することで推定することができるが、計測筋電信号のサンプリング周波数は１０００〜２０００［Ｈｚ］であり信号帯域幅が広く情報量が多いため、現実には筋電計において計測筋電信号は５［Ｈｚ］程度でダウンサンプリングされて情報量を落とした整流平滑化信号に変換される。そのような筋電計から得られる整流平滑化信号から信号強度依存ノイズを推定したり、計測筋電信号を復元したりすることのニーズがあるが、従来の筋電信号モデルでは整流平滑化信号から信号強度依存ノイズを推定することも計測筋電信号を復元することも困難である。
【０００７】
さらに、特に筋電義手制御の分野において、信号強度依存ノイズを考慮して人工筋電信号を生成することへのニーズがある。筋電義手の制御には、サポートベクターマシンやニューラルネットワークが利用され、事前に識別対象動作の筋電信号のパターンを学習させる必要があるが、筋電義手の使用者の負担を考え、軽い発揮筋力での機械学習が行われる。このため、実生活において大きな筋力が発揮されると、信号強度依存ノイズの影響で識別対象動作の誤識別が起こり、筋電義手が誤作動するおそれがある。そのような誤識別や誤作動を防ぐには、実生活において発揮されると考えられるさまざまな大きさの筋力に対応した筋活動度での筋電信号を、信号強度依存ノイズを考慮した形で人工的に生成し、この人工筋電信号を用いて機械学習することが望まれる。
【０００８】
上記問題に鑑み、本発明は、整流平滑化後の筋電信号から信号強度依存ノイズが重畳された整流平滑化前の筋電信号の分散の分布を推定し、さらに当該推定した分布に基づいて人工筋電信号を生成することを課題とする。
【課題を解決するための手段】
【０００９】
本発明の一局面に従った筋電信号処理方法は、時間窓処理部が、信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算するステップと、分散平均推定部が、整流平滑化信号の平均および筋電信号の整流平滑化に係るフィルタゲインの逆数から筋電信号の分散の平均の推定値を計算するステップと、信号強度依存ノイズ分散推定部が、筋電信号の分散の平均の推定値および筋電信号の分散の分布を規定するパラメータから信号強度依存ノイズの分散の推定値を計算するステップとを備え、前記筋電信号の分散の分布は、前記筋電信号の分布に対して共役な性質を有するとともに、分散の非負性を満たす分布である。
【００１０】
同様に、筋電信号処理装置は、信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算する時間窓処理部と、整流平滑化信号の平均および筋電信号の整流平滑化に係るフィルタゲインの逆数から筋電信号の分散の平均の推定値を計算する分散平均推定部と、筋電信号の分散の平均の推定値および筋電信号の分散の分布を規定するパラメータから信号強度依存ノイズの分散の推定値を計算する信号強度依存ノイズ分散推定部とを備え、前記筋電信号の分散の分布は、前記筋電信号の分布に対して共役な性質を有するとともに、分散の非負性を満たす分布である。
【００１１】
同様に、筋電信号処理プログラムは、信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号について所定の時間長ごとに平均を計算する手段、整流平滑化信号の平均および筋電信号の整流平滑化に係るフィルタゲインの逆数から筋電信号の分散の平均の推定値を計算する手段、および筋電信号の分散の平均の推定値および筋電信号の分散の分布を規定するパラメータから信号強度依存ノイズの分散の推定値を計算する手段としてコンピュータを機能させ、前記筋電信号の分散の分布は、前記筋電信号の分布に対して共役な性質を有するとともに、分散の非負性を満たす分布である。
【００１２】
当該筋電信号処理方法、装置およびプログラムによると、信号強度依存ノイズが重畳された筋電信号を整流平滑化した整流平滑化信号から筋電信号の分散の平均の推定値および信号強度依存ノイズの分散の推定値が計算され、これら推定値を用いて、信号強度依存ノイズが重畳された整流平滑化前の筋電信号の分散の分布を推定することができる。
【００１３】
上記筋電信号処理方法は、第１の乱数発生器が、白色正規乱数を発生させるステップと、シェイピングフィルタが、白色正規乱数の周波数特性を筋電信号と同様の周波数特性に整形するステップと、第２の乱数発生器が、筋電信号の分散の平均の推定値および信号強度依存ノイズの分散の推定値で規定される筋電信号の分散の分布に従った乱数を発生させるステップと、人工筋電信号生成部が、整形された乱数および筋電信号の分散の分布に従った乱数から人工筋電信号を生成するステップとを備えていてもよい。
【００１４】
同様に、上記筋電信号処理装置は、白色正規乱数を発生させる第１の乱数発生器と、白色正規乱数の周波数特性を筋電信号と同様の周波数特性に整形するシェイピングフィルタと、筋電信号の分散の平均の推定値および信号強度依存ノイズの分散の推定値で規定される筋電信号の分散の分布に従った乱数を発生させる第２の乱数発生器と、整形された乱数および筋電信号の分散の分布に従った乱数から人工筋電信号を生成する人工筋電信号生成部とを備えていてもよい。
【００１５】
同様に、上記筋電信号処理プログラムは、白色正規乱数を発生させる手段、白色正規乱数の周波数特性を筋電信号と同様の周波数特性に整形する手段、筋電信号の分散の平均の推定値および信号強度依存ノイズの分散の推定値で規定される筋電信号の分散の分布に従った乱数を発生させる手段、および整形された乱数および筋電信号の分散の分布に従った乱数から人工筋電信号を生成する手段としてコンピュータを機能させてもよい。
【００１６】
当該筋電信号処理方法、装置およびプログラムによると、整形された白色正規乱数および筋電信号の分散の分布に従った乱数から信号強度依存ノイズを考慮した人工筋電信号を生成することができる。
【００１７】
上記筋電信号処理方法は、分散平均推定部が、最大随意筋収縮時の整流平滑化信号、指定された筋活動度および筋電信号の整流平滑化に係るフィルタゲインの逆数から当該筋活動度での筋電信号の分散の平均の推定値を計算するステップと、信号強度依存ノイズ分散推定部が、筋電信号の分散の分布を規定するパラメータおよび当該筋活動度での筋電信号の分散の平均の推定値から当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算するステップと、第２の乱数発生器が、当該筋活動度での筋電信号の分散の平均の推定値および当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される当該筋活動度での筋電信号の分散の分布に従った乱数を発生させるステップとを備えていてもよい。
【００１８】
同様に、上記筋電信号処理装置において、分散平均推定部が、最大随意筋収縮時の整流平滑化信号、指定された筋活動度および筋電信号の整流平滑化に係るフィルタゲインの逆数から当該筋活動度での筋電信号の分散の平均の推定値を計算するものであり、信号強度依存ノイズ分散推定部が、筋電信号の分散の分布を規定するパラメータおよび当該筋活動度での筋電信号の分散の平均の推定値から当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算するものであり、第２の乱数発生器が、当該筋活動度での筋電信号の分散の平均の推定値および当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される当該筋活動度での筋電信号の分散の分布に従った乱数を発生させるものであってもよい。
【００１９】
同様に、上記筋電信号処理プログラムは、最大随意筋収縮時の整流平滑化信号、指定された筋活動度および筋電信号の整流平滑化に係るフィルタゲインの逆数から当該筋活動度での筋電信号の分散の平均の推定値を計算する手段、筋電信号の分散の分布を規定するパラメータおよび当該筋活動度での筋電信号の分散の平均の推定値から当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値を計算する手段、および当該筋活動度での筋電信号の分散の平均の推定値および当該筋活動度での筋電信号に重畳される信号強度依存ノイズの分散の推定値で規定される当該筋活動度での筋電信号の分散の分布に従った乱数を発生させる手段としてコンピュータを機能させてもよい。
【００２０】
当該筋電信号処理方法、装置およびプログラムによると、任意の筋活動度を指定して、その筋活動度での人工筋電信号を生成することができる。
【００２１】
上記筋電信号処理方法、装置およびプログラムにおいて、例えば、筋電信号の分散の分布が逆ガンマ分布であってもよく、パラメータが逆ガンマ分布の形状母数であってもよい。
【発明の効果】
【００２２】
本発明によると、整流平滑化後の筋電信号から信号強度依存ノイズが重畳された整流平滑化前の筋電信号の分散の分布を推定することができる。さらに、当該推定した分布に基づいて人工筋電信号を生成することができる。
【図面の簡単な説明】
【００２３】
【図１】本発明で採用する筋電信号モデルの概略図
【図２】本発明の一実施形態に係る筋電信号処理装置を含むシステムの概略図
【図３】筋電信号処理装置のブロック図
【図４】各サンプリング周波数での筋電信号の分散の平均の推定に係る平均誤差率を示すグラフ
【図５】各時間窓での筋電信号の分散の平均の推定に係る平均誤差率を示すグラフ
【図６】各サンプリング周波数での信号強度依存ノイズの分散の推定に係る平均誤差率を示すグラフ
【図７】ある１試行における各負荷での筋電信号の分散の事後分布を示すグラフ
【図８】全試行における各筋力での筋電信号の分散の平均および信号強度依存ノイズの分散の各推定値を示すグラフ
【図９】各負荷での計測筋電信号とそれを元に生成した人工筋電信号を示すグラフ
【図１０】各負荷での計測筋電信号と人工筋電信号の平均振幅の平均絶対誤差率を示すグラフ
【図１１】各負荷での計測筋電信号と人工筋電信号とのパワースペクトル密度の相関係数を示すグラフ
【図１２】各負荷での筋電信号の分散の分布に係るパラメータの平均絶対誤差率を示すグラフ
【発明を実施するための形態】
【００２４】
以下、適宜図面を参照しながら、実施の形態を詳細に説明する。ただし、必要以上に詳細な説明は省略する場合がある。例えば、既によく知られた事項の詳細説明や実質的に同一の構成に対する重複説明を省略する場合がある。これは、以下の説明が不必要に冗長になるのを避け、当業者の理解を容易にするためである。
【００２５】
なお、発明者らは、当業者が本発明を十分に理解するために添付図面および以下の説明を提供するのであって、これらによって特許請求の範囲に記載の主題を限定することを意図するものではない。
【００２６】
≪筋電信号の分散の事後分布の推定≫
（筋電信号モデル）
図１は、本発明で採用する筋電信号モデルの概略図である。当該筋電信号モデルは、時刻ｔにおける筋電信号ｘｔを、シェイピングフィルタＨを通したホワイトガウスノイズＷおよび分散σｔ２に基づいて表現する。ここで、σｔ２は時刻ｔにおける確率変数σ２の値であり、σ２の分布は筋力ｆに応じて決まるσ￣２（本明細書ではσの平均（シグマバー）を便宜上「σ￣」と表す）および信号強度依存ノイズεによって決定される。さらに、ｘｔを整流平滑化した信号ｙｔを用いてσ２の分布を推定可能である。
【００２７】
まず、ｆとσ￣との関係は次のように表現できる。
【００２８】
【数１】
0006652252.tif 000002
【００２９】
ここで、ｋとａは定数であり、実験的に求められる。
【００３０】
次に、σ２はσ￣２とノイズεの和として表される。
【００３１】
【数２】
0006652252.tif 000003
【００３２】
ここで、ノイズεは平均０の確率変数である。
【００３３】
これにより、σ２の平均Ｅ［σ２］と分散Ｖａｒ［σ２］は、
【００３４】
【数３】
0006652252.tif 000004
【００３５】
と表される。
【００３６】
ｘｔはシェイピングフィルタＨを通したホワイトガウスノイズＷとσｔの積で表されるため、平均０、分散σｔ２の正規分布に従う。
【００３７】
【数４】
0006652252.tif 000005
【００３８】
ｙｔはｘｔをＮ次のローパスフィルタを用いて整流平滑化したものである。
【００３９】
【数５】
0006652252.tif 000006
【００４０】
ここで、ａ０，ａ１，・・・，ａＮ−１およびｂ０，ｂ１，・・・，ｂＮはフィルタ係数である。
【００４１】
（分散分布推定）
次に、ｘｔからσ２の分布を推定する方法について説明する。まず、筋電信号ｘｔが計測されたとき、分散の事後分布Ｐ（σ２｜ｘｔ）はベイズの定理を用いると次のように表現できる。
【００４２】
【数６】
0006652252.tif 000007
【００４３】
ここで、尤度Ｐ（ｘｔ｜σ２）は（５）式より、
【００４４】
【数７】
0006652252.tif 000008
【００４５】
である。また、事前分布Ｐ（σ２）にはσ２＞０を考慮して正規分布と共役の事前分布である逆ガンマ分布ＩＧ（α，β）を仮定する。
【００４６】
【数８】
0006652252.tif 000009
【００４７】
ただし、α，βは逆ガンマ分布を決定するパラメータである。このとき、（７）式の右辺分子は次のように展開できる。
【００４８】
【数９】
0006652252.tif 000010
【００４９】
ここで、κ（α，β）は定数項をまとめたものであり、
【００５０】
【数１０】
0006652252.tif 000011
【００５１】
となる。（７）式のＰ（ｘｔ）は定数、Ｐ（σ２｜ｘｔ）は確率密度関数でなければならないため、κ（α，β）とＰ（ｘｔ）がキャンセルされＰ（σ２｜ｘｔ）はＩＧ（α＋１／２，β＋ｘｔ２／２）と一致する。このとき、σ２の事後平均Ｅ［σ２｜ｘｔ］、事後分散Ｖａｒ［σ２｜ｘｔ］は、
【００５２】
【数１１】
0006652252.tif 000012
【００５３】
となる。このとき、（３）式、（４）式よりＥ［σ２｜ｘｔ］とＶａｒ［σ２｜ｘｔ］はｘｔが観測された条件下におけるσ￣２とＶａｒ［ε］の推定値にそれぞれ一致するため、αをあらかじめ設定することでσ￣２からＶａｒ［ε］を求めることができる。
【００５４】
（σ２の推定）
次に、整流平滑化後の信号ｙｔからσ￣２を推定する方法について説明する。ｙｔは（６）式で表現できるため、ｙｔの期待値は次のようになる。
【００５５】
【数１２】
0006652252.tif 000013
【００５６】
ここで、ｘｔは（５）式に従うためＥ［｜ｘｔ｜］は、
【００５７】
【数１３】
0006652252.tif 000014
【００５８】
となる。εｔが時刻ｔにおけるεの値である。また、Ｅ［ｙｔ］は時刻に依存しないとすると、（１４）式は（１５）式を用いて、
【００５９】
【数１４】
0006652252.tif 000015
【００６０】
となる。ここで、（１６）式の右辺第２項は期待値０であり、εｔ（≪σ￣）のｂｉ（＞０）による荷重平均を含むため、第１項より十分小さく、無視できる。これにより、σ￣２は（１６）式を用いて、
【００６１】
【数１５】
0006652252.tif 000016
【００６２】
と推定できる。このとき、右辺の分数係数項はフィルタゲインの逆数に相当する。
【００６３】
以上のことから、整流平滑化された筋電信号ｙｔと形状母数αを用いて（１３）式、（１７）式に基づきσ￣２とＶａｒ［ε］を推定することができる。
【００６４】
≪人工筋電信号の生成≫
整流平滑化後の信号ｙｔから整流平滑化前の信号ｘｔの分散σｔ２の事後分布が推定できれば、当該推定した分布を用いて整流平滑化前の信号ｘｔと同じ統計学的性質を持つ人工筋電信号を生成することができる。
【００６５】
（筋電信号の再現）
時刻ｔにおける人工筋電信号ｚｔは、シェイピングフィルタＨを通した正規乱数ｗ´ｔと分散σｔ２に基づいて表現することができる。
【００６６】
ｗ´ｔは筋電信号ｘｔと同様の周波数特性を持つ正規乱数であり、次数Ｍの自己回帰モデルに基づくシェイピングフィルタＨによって生成される。
【００６７】
【数１６】
0006652252.tif 000017
【００６８】
ここで、ｗｔは平均０、分散１の白色正規乱数、ｖ，ａｊ（ｊ＝１，・・・，Ｍ）はそれぞれ自己回帰モデルにおける推定誤差分散とパラメータである。また、自己回帰モデルパラメータの推定において、分散が１となるように正規化した筋電信号ｘｔを用いることによって、ｗ´ｔは平均０、分散１の正規乱数となる。
【００６９】
そして、人工筋電信号ｚｔは（３）（４）式に基づく分布により生成される乱数σｔとｗ´ｔの積で次のように定義される。すなわち、図１の筋電信号モデルにおいて筋電信号ｘｔを人工筋電信号ｚｔに置換することができる。
【００７０】
【数１７】
0006652252.tif 000018
【００７１】
以上より、（１）式のパラメータｋ，ａとσ２の分布、およびシェイピングフィルタＨのパラメータをあらかじめ設定することで、発揮筋力ｆに応じた人工筋電信号ｚｔを生成することができる。すなわち、整流平滑化後の信号ｙｔから整流平滑化前の信号ｘｔを人工筋電信号ｚｔとして復元（擬似的に再現）することができる。
【００７２】
（任意の筋活動度での人工筋電信号の生成）
さらに、任意の筋活動度での筋電信号を人工的に生成することもできる。例えば、最大随意筋収縮時の整流平滑化後の信号をｙｍａｘ、筋活動度（％ＭＶＣ）をｒｔとすると、
ｒｔ＝ｙｔ／ｙｍａｘ      （２０）
と表される。したがって、ｙｔの期待値Ｅ［ｙｔ］は、
Ｅ［ｙｔ］＝Ｅ［ｒｔ］・ｙｍａｘ      （２１）
と表され、（２１）式を（１７）式に代入すると、
【００７３】
【数１８】
0006652252.tif 000019
【００７４】
となる。すなわち、最大随意筋収縮時の整流平滑化後の信号ｙｍａｘを与え、任意の筋活動度ｒｔを指定すれば、（１３）式、（２２）式に基づいて、筋活動度ｒｔでの筋電信号ｘｔの分散の事後分布を推定することができる。そして、（１９）式に基づいて筋活動度ｒｔでの人工筋電信号ｚｔを生成することができる。
【００７５】
なお、複数チャネルの筋電信号を計測している場合、被験者のある動作中に着目するチャンネルｌおよびその動作に係る主動作筋のチャンネルｍの各筋電信号から各チャンネルの筋活動度の時間平均λｌおよびλｍを計算し、その比率λｌ／λｍを筋活動度ｒｔに乗じた値ｒｔ×λｌ／λｍを（２２）式中のＥ［ｒｔ］として使用してもよい。例えば手を開く動作について考えると、軽く手を開いても思い切り手を開いても各チャネルの筋電信号のパターンは同じようになると考えられるため、上記のように主動作筋の筋活動度の時間平均に対する着目チャンネルの筋活動度の時間平均の比率を乗じて指定された筋活動度を補正することで、動作に対応した特徴的な筋電パターンを崩さずに所望の筋活動度での人工筋電信号を生成することができる。
【００７６】
≪筋電信号処理装置≫
次に、本発明に係る筋電信号処理装置について説明する。図２は、本発明の一実施形態に係る筋電信号処理装置を含むシステムの概略図である。本システム１００は、筋電計１０と筋電信号処理装置２０とで構成される。
【００７７】
筋電計１０は、任意の個数（本実施形態では４個）の筋電センサ１１とデータ収集装置１２とを備えている。データ収集装置１２と筋電センサ１１とはパラレルケーブルまたはシリアルケーブルで接続されている。さらに、これらケーブルは、データ収集装置１２から取り外しできるようになっている。
【００７８】
筋電センサ１１は、被験者の任意の身体部位の表面筋電位信号（筋電信号）を計測するためのものである。筋電信号は筋肉の複合活動電位を波形化して表したものであり、複合活動電位とは筋肉を構成する各筋線維から発生し、容積伝導により体表面に伝搬する活動電位を合計したものをいう。複合活動電位は筋収縮の程度により変化する。例えば、筋肉が小さな力を出している場合には収縮する筋線維の数は少なく、筋線維から発生する活動電位の数も少ない。そこから力を増加させていくと収縮する筋線維の数が徐々に増加し、活動電位の数も増加する。これにより、筋電信号として観測される複合活動電位は増大する。このように、力の出し方と筋電信号とは相関関係にある。
【００７９】
表面信号は２つの表面電極間の電位差として捕捉される。具体的には、各筋電センサ１１は、２個の表面電極１０１（プラス電極とマイナス電極）を有している。これら電極１０１を対象部位の筋繊維の走行方向に沿って数ｃｍ程度の間隔で被験者の体表面に貼付する。なお、いずれか一つの筋電センサ１１は、図略のリファレンス電極をさらに有する。
【００８０】
データ収集装置１２は、図略のアンプ、フィルタ、ＡＤ変換器、ＣＰＵ、通信デバイスなどを備えており、筋電センサ１１が計測した筋電信号（筋電信号モデルにおけるｘｔ）を収集し、当該収集した信号を整流および平滑化して整流平滑化信号（筋電信号モデルにおけるｙｔ）を生成して筋電信号処理装置２０へ送信する。データ収集装置１２と筋電信号処理装置２０との間の通信には、Ｗｉ−Ｆｉ（登録商標）、Ｂｌｕｅｔｏｏｔｈ（登録商標）、ＺｉｇＢｅｅ（登録商標）、特定小電力無線などの無線通信を使用することができる。このようにデータ収集装置１２において計測筋電信号を整流平滑化してデータ量を削減することで、比較的低速の通信回線で多くのチャネルの筋電信号を筋電信号処理装置２０へ送信することができる。
【００８１】
筋電信号処理装置２０は、図１の筋電信号モデルに基づく上記の筋電信号処理を行う装置である。すなわち、筋電信号処理装置２０は、筋電計１０から整流平滑化信号ｙｔを受信して上述の筋電信号処理を行い、筋電センサ１１が計測した筋電信号ｘｔの分散の事後分布を推定し、さらに当該推定した分布に基づいて人工筋電信号を生成する。筋電信号処理装置２０が生成した人工筋電信号は、筋電信号処理装置２０または別の装置のモニター画面にリアルタイムで表示したり、筋電義手制御の機械学習用の信号として別の装置に出力したりすることができる。
【００８２】
筋電信号処理装置２０として、汎用コンピュータ（デスクトップコンピュータやノートコンピュータなど）、スマートフォン、タブレット端末などを使用することができる。これらコンピュータや端末に、上述の筋電信号処理を実現するためのプログラムまたはアプリケーションをインストールすることで、当該コンピュータや端末を筋電信号処理装置２０として利用することができる。
【００８３】
図３は、筋電信号処理装置２０のブロック図である。筋電信号処理装置２０は、時間窓処理部２１、分散平均推定部２２、信号強度依存ノイズ分散推定部２３、ホワイトガウスノイズ発生器２４、シェイピングフィルタ２５、逆ガンマ分布乱数発生器２６および人工筋電信号生成部２７を備えている。このうち、整流平滑化信号ｙｔから筋電センサ１１が計測した筋電信号ｘｔの分散の事後分布を推定するのに必要なブロックは時間窓処理部２１、分散平均推定部２２および信号強度依存ノイズ分散推定部２３の３つであり、残り４つのブロックを追加することで人工筋電信号の生成が可能となる。
【００８４】
時間窓処理部２１は、所定の時間長（時間窓）ごとに整流平滑化信号ｙｔの平均Ｅ［ｙｔ］を計算する。ここで、ｙｔは信号強度依存ノイズεが重畳された筋電信号ｘｔを整流平滑化したものであり、時間窓処理部２１は、実質的に（１４）式の演算を実行する。
【００８５】
分散平均推定部２２は、整流平滑化信号の平均Ｅ［ｙｔ］および筋電信号ｘｔの整流平滑化に係るフィルタゲインの逆数から筋電信号の分散の平均σ￣２の推定値を計算する。すなわち、分散平均推定部２２は、実質的に（１７）式の演算を実行する。
【００８６】
また、任意の筋活動度ｒｔでの人工筋電信号ｚｔを生成する場合、分散平均推定部２２は、最大随意筋収縮時の整流平滑化信号ｙｍａｘ、指定された筋活動度ｒｔおよび筋電信号の整流平滑化に係るフィルタゲインの逆数から筋活動度ｒｔでの筋電信号の分散の平均σ￣２の推定値を計算する。この場合、分散平均推定部２２は、実質的に（２２）式の演算を実行する。なお、分散平均推定部２２は、指定された筋活動度ｒｔをそのまま使用してもよいし、上述したように、主動作筋の筋活動度の時間平均に対する着目チャンネルの筋活動度の時間平均の比率を乗じて補正した値を使用してもよい。
【００８７】
最大随意筋収縮時の筋電信号を計測するには被験者に筋活動度１００（１００％ＭＶＣ）で筋力を発揮してもらう必要があるが、被験者に負担がかかるなどの理由でそれができない場合がある。そのような場合には最大随意筋収縮時の整流平滑化信号を取得できないため、分散平均推定部２２に入力する最大随意筋収縮時の整流平滑化信号ｙｍａｘとして、例えば、他の被験者の実測値などで代用した想定値を用いることができる。もちろん、被験者の最大随意筋収縮時の筋電信号が計測できる場合にはｙｍａｘとして想定値ではなく実測値を用いてもよいことは言うまでもない。
【００８８】
信号強度依存ノイズ分散推定部２３は、筋電信号の分散の平均σ￣２の推定値および筋電信号の分散の分布に係るパラメータから信号強度依存ノイズの分散Ｖａｒ［ε］の推定値を計算する。ここで、筋電信号の分散の分布として逆ガンマ分布を想定すると、当該パラメータとして逆ガンマ分布の形状母数αを採用することができる。すなわち、信号強度依存ノイズ分散推定部２３は、実質的に（１３）式の演算を実行する。
【００８９】
また、任意の筋活動度ｒｔでの人工筋電信号ｚｔを生成する場合、信号強度依存ノイズ分散推定部２３は、筋電信号の分散の分布に係るパラメータおよび筋活動度ｒｔでの筋電信号の分散の平均σ￣２の推定値から筋活動度ｒｔでの筋電信号に重畳される信号強度依存ノイズの分散Ｖａｒ［ε］の推定値を計算する。
【００９０】
ホワイトガウスノイズ発生器２４は、白色正規乱数ｗｔを発生させる。ここで、ｗｔの平均は０、分散は１である。
【００９１】
シェイピングフィルタ２５は、白色正規乱数ｗｔの周波数特性を筋電信号ｘｔと同様の周波数特性に整形して乱数ｗ´ｔを生成する。ここで、ｗ´ｔの平均は０、分散は１である。すなわち、シェイピングフィルタ２５は、実質的に（１８）式の演算を実行する。
【００９２】
逆ガンマ分布乱数発生器２６は、筋電信号の分散の平均σ￣２の推定値および信号強度依存ノイズの分散Ｖａｒ［ε］の推定値で規定される筋電信号の分散の分布に従った乱数σｔを発生させる。ここで、筋電信号の分散の分布として逆ガンマ分布を想定しているため、当該乱数σｔとして逆ガンマ分布乱数を使用する。
【００９３】
また、任意の筋活動度ｒｔでの人工筋電信号ｚｔを生成する場合、逆ガンマ分布乱数発生器２６は、筋活動度ｒｔでの筋電信号の分散の平均σ￣２の推定値および筋活動度ｒｔでの筋電信号に重畳される信号強度依存ノイズの分散Ｖａｒ［ε］の推定値で規定される筋活動度ｒｔでの筋電信号の分散の分布に従った乱数σｔを発生させる。
【００９４】
なお、逆ガンマ分布を規定するパラメータである形状母数αおよび尺度母数βは、例えば、第２種の最尤推定に基づいて、あらかじめ得られた筋電信号系列について周辺尤度を最大化することで求めることができる。
【００９５】
人工筋電信号生成部２７は、シェイピングフィルタ２５で整形後の乱数ｗ´ｔおよび筋電信号の分散の分布に従った乱数σｔから人工筋電信号ｚｔを生成する。すなわち、人工筋電信号生成部２７は、実質的に（１９）式の演算を実行する。
【００９６】
≪実証実験≫
（分散分布推定精度）
筋電信号処理装置２０による分散分布推定精度を確認するため人工的に生成した人工筋電信号を用いて実験を行った。人工筋電信号の生成は以下の手順で行うことで、分散にノイズが重畳した信号を擬似的に生成した。
（１）長さＴの離散系列｛σｔ２；ｔ＝１，・・・，Ｔ｝を逆ガンマ分布ＩＧ（α０，β０）に従う乱数で発生させる。
（２）各ｔにおいて正規分布Ｎ（０，σｔ２）に従う乱数を１サンプルずつ発生させ、これを系列｛ｘｔ｝とする。
（３）｛ｘｔ｝をサンプリング周波数Ｆｓ［Ｈｚ］で計測された筋電信号とみなす。
このとき、｛ｘｔ｝の分散の平均、および分散の分散の真値はそれぞれ以下の値となる。
【００９７】
【数１９】
0006652252.tif 000020
【００９８】
分散分布推定精度検証は、これら真値と筋電信号処理装置２０で推定したσ￣２およびＶａｒ［ε］とを比較することによって行った。
【００９９】
σ￣２の推定は、｛ｘｔ｝を整流平滑化した信号を｛ｙｔ｝とし、｛ｙｔ｝をＦｄｏｗｎ［Ｈｚ］でダウンサンプリングした信号のうち最初のＬ［ｍｓ］を用いて行った。このとき、Ｌを１０００［ｍｓ］に固定してＦｄｏｗｎを５００［Ｈｚ］、２００［Ｈｚ］、１００［Ｈｚ］、５０［Ｈｚ］、２０［Ｈｚ］、１０［Ｈｚ］と変更した場合、およびＦｄｏｗｎを１０００［Ｈｚ］に固定してＬを１０００［ｍｓ］、５００［ｍｓ］、２００［ｍｓ］、１００［ｍｓ］、５０［ｍｓ］、２０［ｍｓ］、１０［ｍｓ］と変更した場合の各条件において、真値σ￣０２とσ￣２との％誤差を算出した。比較手法として、同様の条件において｛ｘｔ｝から算出した分散の最尤推定値を用いた。
【０１００】
Ｖａｒ［ε］の推定精度検証では、Ｌを１０００［ｍｓ］に固定してＦｄｏｗｎを１０００［Ｈｚ］、５０［Ｈｚ］と変更した各条件において、Ｖａｒ［ε０］との％誤差を算出した。また、最尤推定に基づく分散の分散との比較を行った。ここで、最尤推定に基づく分散の分散とは、Ｌ［ｍｓ］の｛ｘｔ｝をＫ分割し、各分割において算出したＫ個の分散最尤推定値の分散のことをいう。なお、事前分布パラメータは経験ベイズ法を用いて決定した。
【０１０１】
各実験において、Ｔ＝１０００、Ｆｓ＝１０００とし、平滑化処理にはカットオフ周波数１［Ｈｚ］の２次のバターワースフィルタを使用した。また、％誤差は真値を２０通り（α＝１５、β＝０．５，１，１．５，・・・，１０）変更した際の平均値である。
【０１０２】
図４は、各サンプリング周波数での筋電信号の分散の平均σ￣２の推定に係る平均誤差率を示すグラフである。図５は、各時間窓での筋電信号の分散の平均σ￣２の推定に係る平均誤差率を示すグラフである。図４のグラフからはＦｄｏｗｎ≦２００で、図５のグラフからはＬ≦１００で、本発明（図１の筋電信号モデルに基づく推定）と従来技術（最尤推定）との間に有意差が認められた。図６は、各サンプリング周波数での信号強度依存ノイズの分散Ｖａｒ［ε］の推定に係る平均誤差率を示すグラフである。Bonferroni法に基づく多重比較検定を行った結果、本発明と従来技術との間に有意差が確認できた。
【０１０３】
図４のグラフを参照すると、Ｆｄｏｗｎが高い場合には本発明および従来技術ともに誤差率５［％］以下であり推定精度が高いが、Ｆｄｏｗｎが下がるにつれ従来技術の精度は悪化し、本発明は高精度を維持する。また、図５のグラフを参照すると、窓幅Ｌを短くすると本発明および従来技術ともに精度が悪化するが、Ｌ＝１０［ｍｓ］においても本発明の推定誤差は５［％］程度であり、従来技術よりも有意に小さい。このことから、本発明によると、より少ないサンプルでσ￣２を推定できると言える。
【０１０４】
また、図６のグラフを参照すると、本発明は従来技術と比較して、Ｖａｒ［ε］の推定誤差が有意に小さい。これは、従来技術の最尤法が各分割において分散一定を仮定するため、各時刻独立に重畳するノイズを十分に表現できないことが原因である。一方、本発明では各時刻独立の分散値を仮定し、それらの母集団の事後分布を推定することで、分散の分散を誤差率１０［％］以下の良好な精度で推定することができる。
【０１０５】
以上のように、本発明の一実施形態に係る筋電信号処理装置２０によると、信号強度依存ノイズが重畳された筋電信号の分散の分布を高精度に推定することができる。
【０１０６】
次に、実データに対する本発明の有効性を確認するために、筋電信号を用いて分散事後分布Ｐ（σ２｜ｘ）の推定実験を行った。筋電信号の計測では、被験者５名（平均年齢：２３．３±０．８、右利き）の上腕二頭筋上の皮膚表面にＡｇ／ＡｇＣｌ電極を貼付し、座位で上腕を下垂し手掌を上に向けて前腕を水平前方に曲げた姿勢をとらせ、肘を机につき手首に負荷をのせ姿勢を１０秒間維持させた際の筋電信号をサンプリング周波数１０００［Ｈｚ］で取得した。このとき、負荷を５００［ｇ］、１０００［ｇ］、１５００［ｇ］、２０００［ｇ］と変更し、それぞれ５試行計測した。なお、計測には日本光電工業株式会社製マルチテレメータシステムＷＥＢ−５０００（高域遮断周波数：１００［Ｈｚ］、低域遮断周波数：５．３［Ｈｚ］）を用いた。
【０１０７】
筋電信号の解析には、計測した１０秒間のうち後半の５秒間のデータを用いた。このとき、前腕および手の自重を体重比０．０２２、質量中心を肘関節からの長さ比０．３１８、上腕二頭筋の肘関節回転中心に対するモーメントアーム長を０．０３［ｍ］と仮定し、上腕二頭筋の発揮筋力を計算した。そして、負荷の増加に伴う分散事後分布Ｐ（σ２｜ｘ）の変化と、筋力の変化に伴うＥ［σ２｜ｘ］およびＶａｒ［σ２｜ｘ］の変化を求めた。
【０１０８】
なお、事前分布パラメータはあらかじめ計測した２０００［ｇ］の負荷時の筋電信号１試行分を用いて、経験ベイズ法に基づき被験者ごとに設定した。
【０１０９】
図７は、ある１試行における各負荷での筋電信号の分散の事後分布を示すグラフである。図８は、全試行における各筋力での筋電信号の分散の平均および信号強度依存ノイズの分散の各推定値を示すグラフである。図８にはSteel-Dwass法に基づく検定結果を併記している。
【０１１０】
図７のグラフを参照すると、負荷が大きくなるにつれてＰ（σ２｜ｘ）が右に移動し広がりが大きくなっている。これは、負荷の増加に伴い、事後平均と事後分散が増加傾向にあることを意味する。また、図８のグラフを参照すると、筋力と事後平均Ｅ［σ２｜ｘ］が単調増加傾向にあることから、（１）式で筋力が表現可能であることがわかる。このとき、（１）式中の指数ａを求めると、被験者順にａ＝１．１８，０．５６，０．８１，０．７４，０．６５となった。一方、非特許文献１と同様に最尤推定に基づく求めた分散σ２を使用して指数ａを求めると、ａ＝１．１７，０．５５，０．７８，０．７３，０．６５となった。このことから、本発明で算出できる事後平均Ｅ［σ２｜ｘ］は非特許文献１の方法において最尤法で求められる分散と同様の役割を果たしていることがわかる。
【０１１１】
一方、筋力に応じた事後分散Ｖａｒ［σ２｜ｘ］の増加は事後平均Ｅ［σ２｜ｘ］の増加と（１３）式により明らかであるが、特筆すべきは、本発明において事後分散Ｖａｒ［σ２｜ｘ］が事後平均Ｅ［σ２｜ｘ］の２乗に比例する点である。過去に、等尺性収縮時における筋力の分散が平均筋力のおよそ２乗に比例するという実験結果が報告されているが、本発明は、筋電信号の分散について同様の関係が成立することを示している。このことは、信号強度依存ノイズにより筋電信号の分散値に対してもノイズが重畳することを意味している。
【０１１２】
以上のように、本発明の一実施形態に係る筋電信号処理装置２０によると、筋電信号の分散分布を推定でき、さらに分散に重畳するノイズを推定することができる。
【０１１３】
（人工筋電信号の生成実験）
次に、筋電信号処理装置２０による人工筋電信号の生成実験を行った。実験では、まず健常大学生１名から筋電信号を計測し、分散分布パラメータσ￣２とＶａｒ［ε］、およびシェイピングフィルタ２５のパラメータを推定した。そして、推定したパラメータに基づき生成した人工筋電信号と計測した筋電信号との比較を行った。
【０１１４】
筋電信号の計測条件は上記と同じである。負荷を５００［ｇ］、１０００［ｇ］、１５００［ｇ］、２０００［ｇ］と変更し、それぞれ１試行計測した。パラメータ推定および比較には計測した１０秒間のうち後半の５秒間のデータを使用した。
【０１１５】
分散事後分布の推定は負荷ごとに行った。ただし、事前パラメータはあらかじめ計測した負荷２０００［ｇ］時の筋電信号を用いて第２種の最尤推定に基づき設定し、実験を通じて共通とした。また、シェイピングフィルタ２５のパラメータは負荷２０００［ｇ］時の筋電信号をBurg法によって決定し、モデル次数はベイズ情報量規準に基づきＭ＝２０とした。人工筋電信号の生成は負荷ごとに５回ずつ行い、逆ガンマ分布乱数σｔの生成にはTanizaki法を用いた。
【０１１６】
比較は、平均振幅、周波数成分、分散分布の３項目について行った。平均振幅は信号を整流平滑化した値の平均値とし、計測筋電信号における値を真値としたときの平均絶対誤差率を各負荷において算出した。周波数成分の比較では、計測筋電信号と人工筋電信号のパワースペクトル密度をそれぞれ求め、それらの間の相関関係を算出した。分散分布の比較については、筋電信号の分散の分布を逆ガンマ分布と仮定することでパラメータα、βを第２種の最尤推定によって推定し、平均振幅の場合と同様に平均絶対誤差率を算出した。
【０１１７】
図９は、各負荷での計測筋電信号とそれを元に生成した人工筋電信号を示すグラフである。図９のグラフを参照すると、負荷が大きくなるにつれて計測筋電信号と人工筋電信号の振幅がともに大きくなることがわかる。
【０１１８】
図１０は、各負荷での計測筋電信号と人工筋電信号の平均振幅の平均絶対誤差率を示すグラフである。図１０のグラフを参照すると、人工筋電信号はどの負荷においても５［％］程度の平均絶対誤差率で筋電信号の振幅を再現できていることが確認できる。
【０１１９】
図１１は、各負荷での計測筋電信号と人工筋電信号とのパワースペクトル密度の相関係数を示すグラフである。図１１のグラフを参照すると、すべての負荷において相関係数は０．８以上であり、周波数成分においても人工筋電信号が高い精度で筋電信号を再現できていることがわかる。
【０１２０】
図１２は、各負荷での筋電信号の分散の分布に係るパラメータの平均絶対誤差率を示すグラフである。図１２のグラフを参照すると、形状母数αは負荷が大きくなるにつれて平均誤差率が大きくなる傾向があるものの０．０１［％］未満と非常に小さい。一方、尺度母数βの平均誤差率はαの平均誤差率よりも大きい５−１０［％］程度である。この原因として、（１３）式および（１７）式においてαは事後分布の導出に直接用いられているが、βはそれに用いられていないため、事後分布の尺度に関して推定誤差が大きくなったことが考えられる。しかし、形状母数αは分布の本質的な形を決定するのに対して尺度母数βは分布の広がりを決定するものであるため、βが人工筋電信号の生成に及ぼす影響はαに比べて小さく問題にならないと考える。
【０１２１】
以上のように、本発明の一実施形態に係る筋電信号処理装置２０によると、計測筋電信号の振幅、周波数、分散分布を再現する人工筋電信号を生成することができる。
【０１２２】
以上のように、本発明における技術の例示として、実施の形態を説明した。そのために、添付図面および詳細な説明を提供した。
【０１２３】
したがって、添付図面および詳細な説明に記載された構成要素の中には、課題解決のために必須な構成要素だけでなく、上記技術を例示するために、課題解決のためには必須でない構成要素も含まれ得る。そのため、それらの必須ではない構成要素が添付図面や詳細な説明に記載されていることをもって、直ちに、それらの必須ではない構成要素が必須であるとの認定をするべきではない。
【０１２４】
また、上述の実施の形態は、本発明における技術を例示するためのものであるから、特許請求の範囲またはその均等の範囲において種々の変更、置き換え、付加、省略などを行うことができる。例えば、上記では筋電信号の分散の分布として逆ガンマ分布を想定したが、筋電信号の分散の分布として逆ガンマ分布以外に、筋電信号の分布に対して共役な性質を有するとともに、分散の非負性を満たす分布を適用してもよい。
【符号の説明】
【０１２５】
２０  筋電信号処理装置
２１  時間窓処理部
２２  分散平均推定部
２３  信号強度依存ノイズ分散推定部
２４  ホワイトガウスノイズ発生器（第１の乱数発生器）
２５  シェイピングフィルタ
２６  逆ガンマ分布乱数発生器（第２の乱数発生器）
２７  人工筋電信号生成部
【図１】
0006652252.tif 000021
【図２】
0006652252.tif 000022
【図３】
0006652252.tif 000023
【図４】
0006652252.tif 000024
【図５】
0006652252.tif 000025
【図６】
0006652252.tif 000026
【図７】
0006652252.tif 000027
【図８】
0006652252.tif 000028
【図９】
0006652252.tif 000029
【図１０】
0006652252.tif 000030
【図１１】
0006652252.tif 000031
【図１２】
0006652252.tif 000032