(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2023086656
(43)【公開日】2023-06-22
(54)【発明の名称】成分波形抽出システム、成分波形抽出方法及びコンピュータを成分波形抽出システムとして動作させるプログラム、並びに信号解析システム、信号解析方法及びコンピュータを信号解析システムとして動作させるプログラム
(51)【国際特許分類】
A61B 5/397 20210101AFI20230615BHJP
【FI】
A61B5/397
【審査請求】未請求
【請求項の数】32
【出願形態】OL
(21)【出願番号】P 2022105826
(22)【出願日】2022-06-30
(31)【優先権主張番号】P 2021200605
(32)【優先日】2021-12-10
(33)【優先権主張国・地域又は機関】JP
(71)【出願人】
【識別番号】504174135
【氏名又は名称】国立大学法人九州工業大学
(74)【代理人】
【識別番号】100197642
【弁理士】
【氏名又は名称】南瀬 透
(74)【代理人】
【識別番号】100099508
【弁理士】
【氏名又は名称】加藤 久
(74)【代理人】
【識別番号】100182567
【弁理士】
【氏名又は名称】遠坂 啓太
(74)【代理人】
【識別番号】100219483
【弁理士】
【氏名又は名称】宇野 智也
(72)【発明者】
【氏名】永井 秀利
【テーマコード(参考)】
4C127
【Fターム(参考)】
4C127AA04
4C127FF02
4C127GG11
(57)【要約】
【課題】計測した信号波形に含まれ、ターゲットとする成分波形の取得や、この取得した成分波形の各特徴を定量的に解析することができる信号波形解析システム等の提供。
【解決手段】成分波形抽出システムは、信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得る解析部10と、係数集合の係数群に対する成分波形の適合度を算出して、適合度に基づいて成分波形を検出し、成分波形に対応する係数群を逆ウェーブレット変換により波形再現することで信号波形に含まれる成分波形を得る成分波形抽出部20とを有する。
【選択図】
図4
【特許請求の範囲】
【請求項1】
信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得る解析部と、
前記係数集合の係数群に対する成分波形の適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形に対応する係数群を逆ウェーブレット変換により波形再現することで前記信号波形に含まれる成分波形を得る成分波形抽出部とを有する、
成分波形抽出システム。
【請求項2】
信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得る解析部と、
前記係数集合の係数群に対する成分波形との適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形の特性に基づく特徴量を分析する特徴量分析部とを有する、
信号解析システム。
【請求項3】
前記適合度が、前記係数集合の係数群に成分波形の基本特性に基づいて得られるものである、請求項1に記載の成分波形抽出システム。
【請求項4】
前記基本特性が、成分波形の特性を示す基本ベクトルであり、
前記適合度は、前記基本ベクトルと前記係数集合の係数群に基づくベクトルとにより得られるものである請求項3に記載の成分波形抽出システム。
【請求項5】
前記係数集合の係数群に、成分波形に基づいて重み付けされた適合フィルタが適用された係数集合の係数群に対し、成分波形との適合度を算出する請求項4に記載の成分波形抽出システム。
【請求項6】
前記適合度と、成分波形の基本特性に基づいて重み付けされた調整フィルタを用いて、前記係数集合の係数群に包含された複数の波形と成分波形を分離し、成分波形に対応したウェーブレット係数群に調整するものである請求項5に記載の成分波形抽出システム。
【請求項7】
前記調整フィルタは、前記係数集合の係数群を、前記適合度に基づいて検出された複数の成分波形のそれぞれに対応する係数集合の係数群に分配する分配型抽出フィルタを有するものである請求項6に記載の成分波形抽出システム。
【請求項8】
前記調整フィルタは、前記係数集合の係数群から、前記適合度に基づいて検出された成分波形に対応するウェーブレット係数群を導出する導出型抽出フィルタを適用するものである請求項7に記載の成分波形抽出システム。
【請求項9】
前記成分波形抽出部により検出された成分波形の特性に基づく特徴量を解析する特徴量分析部を有し、
成分波形の強度および/または性質を示す特徴量により、前記適合度を補正する、請求項8に記載の成分波形抽出システム。
【請求項10】
前記強度を示す特徴量は、成分波形の信号の強度を示す信号強度特徴量であり、前記性質を示す特徴量は、成分波形に含まれる周波数成分の分布に基づく周波数特徴量である請求項9に記載の成分波形抽出システム。
【請求項11】
前記特徴量分析部が、前記適合度に基づいて、前記信号波形に含まれる複数の成分波形を一つのブロックとして得る請求項9または10に記載の成分波形抽出システム。
【請求項12】
前記解析部は、
前記信号波形に含まれる複数の副波形成分について、それぞれのウェーブレット係数群を算出して、それぞれのウェーブレット係数集合に基づいて係数集合の係数群を得るものであり、
前記成分波形抽出部は、
それぞれの前記係数集合の係数群に対する成分波形の適合度を算出するものである、請求項9または10に記載の成分波形抽出システム。
【請求項13】
前記解析部は、
前記信号波形に含まれる複数の副波形成分を連結させたものを1つの成分波形としてウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて係数集合の係数群を得るものである、請求項9または10に記載の成分波形抽出システム。
【請求項14】
前記適合度が、前記係数集合の係数群に成分波形の基本特性に基づいて得られるものである、請求項2に記載の信号解析システム。
【請求項15】
前記基本特性が、成分波形の特性を示す基本ベクトルであり、
前記適合度は、前記基本ベクトルと前記係数集合の係数群に基づくベクトルとにより得られるものである請求項14に記載の信号解析システム。
【請求項16】
前記係数集合の係数群に、成分波形に基づいて重み付けされた適合フィルタが適用された係数集合の係数群に対し、成分波形との適合度を算出する請求項15に記載の信号解析システム。
【請求項17】
前記適合度と、成分波形の基本特性に基づいて重み付けされた調整フィルタを用いて、前記係数集合の係数群に包含された複数の波形と成分波形を分離し、成分波形に対応したウェーブレット係数群に調整するものである請求項16に記載の信号解析システム。
【請求項18】
前記調整フィルタは、前記係数集合の係数群を、前記適合度に基づいて検出された複数の成分波形のそれぞれに対応する係数集合の係数群に分配する分配型抽出フィルタを有するものである請求項17に記載の信号解析システム。
【請求項19】
前記調整フィルタは、前記係数集合の係数群から、前記適合度に基づいて検出された成分波形に対応するウェーブレット係数群を導出する導出型抽出フィルタを適用するものである請求項18に記載の信号解析システム。
【請求項20】
前記特徴量分析部により求められる前記成分波形の強度および/または性質を示す特徴量により、前記適合度を補正する、請求項19に記載の信号解析システム。
【請求項21】
前記強度を示す特徴量は、成分波形の信号の強度を示す信号強度特徴量であり、前記性質を示す特徴量は、成分波形に含まれる周波数成分の分布に基づく周波数特徴量である請求項20に記載の信号解析システム。
【請求項22】
前記特徴量分析部が、前記適合度に基づいて、前記信号波形に含まれる複数の成分波形を一つのブロックとして得る請求項18または19に記載の信号解析システム。
【請求項23】
前記特徴量分析部は、前記ブロック内の、成分波形の変動平均と変動変化量を分析すると共に、前記変動平均および変動変化量に基づいて、ブロックの変動係数および変動度を解析するブロック変動解析部を有するものである請求項22に記載の信号解析システム。
【請求項24】
前記ブロック変動解析部は、前記ブロックの変動度に基づいて、ブロックの変動に対する、成分波形の活性度を推定するものである請求項23に記載の信号解析システム。
【請求項25】
前記特徴量分析部は、前記ブロックの変動度と前記ブロックを構成する成分波形の信号強度特徴量に基づいて、ブロックの活動係数およびブロックの活動度を解析するブロック出力特性解析部を有する請求項24に記載の信号解析システム。
【請求項26】
前記特徴量分析部は、ブロックの発生度およびブロックの特性に基づいて、ブロックの活動余力度を解析する集散度分析部を有する請求項25に記載の信号解析システム。
【請求項27】
前記ブロック出力特性解析部は、前記信号波形に含まれる複数の成分波形のうち異なる特性を示す成分波形をブロックとして得て、当該異なる特性に基づく当該ブロックの変動度と当該ブロックを構成する成分波形の信号強度特徴量に基づいて、前記活動度を解析するものである、請求項25に記載の信号解析システム。
【請求項28】
前記特徴量分析部は、前記信号波形に含まれる複数の成分波形のうち異なる特性を示す成分波形をブロックとして得て、当該異なる特性に基づくブロックの発生度およびブロックの特性に基づいて、前記活動余力度を解析するものである、請求項25に記載の信号解析システム。
【請求項29】
成分波形抽出システムにより、
信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得るステップと、
前記係数集合の係数群に対する成分波形の適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形に対応する係数群を逆ウェーブレット変換により波形再現することで前記信号波形に含まれる成分波形を得るステップとを有する、
成分波形抽出方法。
【請求項30】
信号解析システムにより、
信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得るステップと、
前記係数集合の係数群に対する成分波形との適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形の特性に基づく特徴量を分析するステップとを有する、
信号解析方法。
【請求項31】
コンピュータを、
信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得る解析部と、
前記係数集合の係数群に対する成分波形の適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形に対応する係数群を逆ウェーブレット変換により波形再現することで前記信号波形に含まれる成分波形を得る成分波形抽出部とを有する、
成分波形抽出システムとして動作させるプログラム。
【請求項32】
コンピュータを、
信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得る解析部と、
前記係数集合の係数群に対する成分波形との適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形の特性に基づく特徴量を分析する特徴量分析部とを有する、
信号解析システムとして動作させるプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、発生頻度が不規則な単発波からなる混合信号波形等を分析可能であり、そのような信号波形の成分波形を抽出する成分波形抽出システム及びその方法等、信号波形の成分波形を分析する信号解析システム及びその方法等に関する。
【背景技術】
【0002】
ある信号波形に含まれる、周期等、規則性に乏しい複数の単発波等が重なり合って形成される波形(成分波形)をターゲットとし、当該成分波形の挙動を把握しようとすることがある。このような解析対象の信号波形としては、例えば、黙声認識(いわゆる、口パク)や、心音、心電図による電気的波形、筋肉の活動等、を捉える筋電波形や、工学的には、機器や各種環境から生じる騒音、振動、または機器等の故障検知を行う異常波形などの信号波形などがあげられる。
【0003】
観測される信号波形から情報を得ようとする行為は極めて一般的である。このとき、その信号波形が、より小さな有限長の単位波形(成分波形)から構成される、あるいは含んでいる場合、成分波形の挙動を調べることで信号波形が示す情報をより高精度に捉えたいと考えることは極めて自然なことである。しかし、成分波形の形状が単純ではなく周波数成分の時間分布が不揃いであったり、成分波形に許容される変化の範囲が存在したり、成分波形以外の様々な波形が含まれていたり、成分波形の出現に周期性、規則性がなかったり、時間差のある複数の成分波形が重なり合ったりといった状況が多く存在するほど、成分波形を検出して挙動を捉えることは極めて困難になる。
【0004】
そのような成分波形の挙動を掴むために、特定のルールで計測条件を違えた多チャネルの計測データを使用し、解析モデルを用いて推定する手法も存在するが、多チャネル計測のコストや、モデルの精度や計算コストといった課題が多く、抽象的な特徴量で信号波形の特性を得るというような妥協をせざるを得ないことが一般的と言える。このような解析対象の信号波形としては、各種筋肉の活動や筋音、単純なsin波の合成としては捉えづらい心電図波形や心音、等の筋電波形や機器等で生じる駆動音や振動、故障等を示唆する異常波形などが挙げられる。
【0005】
これらの波形には、分析対象とする成分波形以外の、その他様々な波形が含まれており、分析を行うためには先ず、その中から当該成分波形を抽出する必要がある。更にまた、この成分波形は、信号波形と共に、特にその周波数は時間変動を伴うなど単純ではなく、時間毎の周波数変動を加味した挙動の把握は難しい。
【0006】
ここで、このような信号波形の解析に関し、一つの具体例として筋電の測定、解析技術に沿って説明する。一つの筋は、多数の筋繊維で構成される。ただし、各筋繊維の活動は独立ではなく、一つの神経細胞によって支配される筋繊維の集まり(以下、「運動単位」ともいう。)を単位として活動が管理されており、各運動単位が時間差を持って重なり合いながら活動することで筋活動が成立している。一つの運動単位を構成する筋繊維は同時に活動すると考えると、詳細な筋活動を捉えるには運動単位活動を捉えれば良く、各筋繊維で生じた電位変化の重畳である運動単位活動電位波形がその道具立てとなる。
【0007】
筋電計測は、各種計測方法が存在し、例えば人体の骨格を形成する筋肉の動きを測定する方法として大きく分けると針筋電と表面筋電とが存在する。針筋電は電極針を筋に刺入することで計測するものであるため、目標運動単位の活動電位をかなり直接的に計測できるが、痛みを伴うために運動中の適用は困難で、気軽に用いるのは難しい。また、得られるデータが刺入位置近辺の局所的な情報だけであり、筋全体の活動を捉えるには情報不足である点も問題である。それに対し、表面筋電は皮膚表面に電極を貼り付けるだけの非侵襲な計測で取り扱いも容易であり、広範囲の運動単位活動を捉えることが可能だが、計測されるのは複数の運動単位活動電位の重畳となる。
【0008】
詳細な筋活動を捉えるためにはこのような運動単位活動を得ることが望ましく、それによって生じる電位変化を捉えるのが近道である。また、多くの情報を内包する表面筋電から運動単位活動を得られれば非常に有用であると言える。しかし、従来技術では、前述したように多チャネルの信号やモデルを用いて推定するような手法がほとんどで、運動単位のタイプの違いや活動タイミングなどを直接的に獲得できているとは言い難い状況にあると考える。
【0009】
なお、筋活動の解析に関る技術としては、次のような文献が開示されている。非特許文献1は、筋組織、筋活動に関連する一般情報に関する。非特許文献2は、多極での計測による評価の試みであり、基本的には同種運動時の相対比較するものである。非特許文献3には、「表面筋電図からパワー、筋疲労、筋活動量の関連性に関する調査報告はない」との記述がなされている。
【先行技術文献】
【特許文献】
【0010】
【特許文献1】特開2018-047230号公報
【特許文献2】特開2002-272692号公報
【特許文献3】特開2007-202612号公報
【特許文献4】特開2016-063995号公報
【特許文献5】WO2004/023996号公報
【特許文献6】特公平05-032057号公報
【特許文献7】特開2013-244027号公報
【非特許文献】
【0011】
【非特許文献1】福永哲夫編「筋の科学辞典~構造・機能・運動~」朝倉書店(2002)
【非特許文献2】白井礼、水戸和幸「多点電極を用いた表面筋電図による動的運動時の筋疲労評価」信学技報Vol.114 No.79、MBE2014-16、p.19-22(2014)
【非特許文献3】徳安達士、松野航大、松平慎平「表面筋電位のパターンエントロピー分析による筋疲労の定量評価」、日本機械学会ロボティクス・メカトロニクス講演会2012講演論文集、2A2-B04(2012)
【非特許文献4】Hidetoshi Nagai, Computational Complexity in Continuous Estimation of Muscular Activity based on Redundant Wavelet Coefficients of Surface EMG, Proc. of Life Engineering Symposium 2015, pp. 329-334, 2015.
【非特許文献5】永井秀利「表面筋電の特徴抽出のためのウェーブレット係数集合のシフト選択法」信学技報, vol. 119, no. 391, MBE2019-78, pp. 51-56(2020)
【発明の概要】
【発明が解決しようとする課題】
【0012】
単発波が固定された周波数で完全に離散的に生じるという性質が期待できない波形は、例えば、前述した筋電に含まれる運動単位活動の波形を含めて多く存在すると考えられ、そのような波形は、その波形に含まれる単発波の信号成分が時間的、周波数的に広がりを持って存在する。
【0013】
そのため、特許文献1~3で開示されているように、連続ウェーブレット変換結果を2次元又は3次元的に表示する場合、どの部分が単発波の信号成分の一部でどの部分がノイズ的な成分かといった判断が難しい。
【0014】
計算機で信号解析を行う場合に用いられるのは、サンプリングによって得られた離散データであり、離散データに対するウェーブレット解析は離散ウェーブレット解析と言われる。離散ウェーブレット解析で一般的な多重解像度解析結果を単純に用いる場合、時間-周波数平面の分割に関する制約により、時間-周波数平面上での単発波の特定は極めて困難となる。
【0015】
特許文献4に記載されているように、冗長離散ウェーブレット解析を行えば、多重解像度解析の問題を回避できるが、連続ウェーブレット変換と類似した問題が残る。これに関し、単一チャネル信号から要素信号(例えば、単チャネルの表面筋電信号から構成要素たる運動単位活動電位)を検出、分離可能にすべく、多チャネルの信号とモデルないしは特別な計測条件とを用いる手法が存在する(特許文献5、6参照)。
【0016】
しかしながら、これらの手法は複雑な計算を必要とし、ノイズ等のモデルの条件に適合しない要素が混入した信号に対して、要素信号の推定精度が低下する。さらに深刻な問題は、推定の際に多チャネルの特徴情報を使うので、得られた要素信号に基づく分析を多チャネルの観点から行えないことである。単一チャネルの信号からの抽出、分離が実現できれば、多チャネルを活用した分析を目指す技術(例えば、特許文献7)はさらに発展できる。
【0017】
なお、前述した非特許文献に関し、非特許文献4は、冗長離散ウェーブレット解析のアルゴリズムに関する基本となる論文である。また、非特許文献5は、ウェーブレット係数抽出時刻をずらす行為の一般化を試みた論文である。現在のウェーブレット係数集合の基本となる論文だが、係数集合と成分波形の関連性については明確には論じられておらず、係数集合が成分波形抽出につながるという考えには至っていない。
【0018】
本発明は、かかる事情に鑑みてなされたもので、計測した信号波形に含まれ、ターゲットとする成分波形の取得や、この取得した成分波形の各特徴を定量的に解析することができる信号波形解析方法や、システム、プログラムを提供することを目的とする。
【課題を解決するための手段】
【0019】
本発明者は、上記課題を解決すべく鋭意研究を重ねた結果、下記の発明が上記目的に合致することを見出し、本発明に至った。すなわち、本発明は、以下の発明に係るものである。
【0020】
<1>信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得る解析部と、
前記係数集合の係数群に対する成分波形の適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形に対応する係数群を逆ウェーブレット変換により波形再現することで前記信号波形に含まれる成分波形を得る成分波形抽出部とを有する、
成分波形抽出システム。
<2>信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得る解析部と、
前記係数集合の係数群に対する成分波形との適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形の特性に基づく特徴量を分析する特徴量分析部とを有する、
信号解析システム。
<3>前記適合度が、前記係数集合の係数群に成分波形の基本特性に基づいて得られるものである、<1>に記載の成分波形抽出システム。
<4>前記基本特性が、成分波形の特性を示す基本ベクトルであり、
前記適合度は、前記基本ベクトルと前記係数集合の係数群に基づくベクトルとにより得られるものである<3>に記載の成分波形抽出システム。
<5>前記係数集合の係数群に、成分波形に基づいて重み付けされた適合フィルタが適用された係数集合の係数群に対し、成分波形との適合度を算出する<1>、<3>、および、<4>のいずれかに記載の成分波形抽出システム。
<6>前記適合度と、成分波形の基本特性に基づいて重み付けされた調整フィルタを用いて、前記係数集合の係数群に包含された複数の波形と成分波形を分離し、成分波形に対応したウェーブレット係数群に調整するものである<1>、および、<3>~<5>のいずれかに記載の成分波形抽出システム。
<7>前記調整フィルタは、前記係数集合の係数群を、前記適合度に基づいて検出された複数の成分波形のそれぞれに対応する係数集合の係数群に分配する分配型抽出フィルタを有するものである<6>に記載の成分波形抽出システム。
<8>前記調整フィルタは、前記係数集合の係数群から、前記適合度に基づいて検出された成分波形に対応するウェーブレット係数群を導出する導出型抽出フィルタを適用するものである<6>または<7>に記載の成分波形抽出システム。
<9>前記成分波形抽出部により検出された成分波形の特性に基づく特徴量を解析する特徴量分析部を有し、
成分波形の強度および/または性質を示す特徴量により、前記適合度を補正する、<1>、および、<3>~<8>のいずれかに記載の成分波形抽出システム。
<10>前記強度を示す特徴量は、成分波形の信号の強度を示す信号強度特徴量であり、前記性質を示す特徴量は、成分波形に含まれる周波数成分の分布に基づく周波数特徴量である<9>に記載の成分波形抽出システム。
<11>前記特徴量分析部が、前記適合度に基づいて、前記信号波形に含まれる複数の成分波形を一つのブロックとして得る<9>または<10>に記載の成分波形抽出システム。
<12>前記解析部は、前記信号波形に含まれる複数の副波形成分について、それぞれのウェーブレット係数群を算出して、それぞれのウェーブレット係数集合に基づいて係数集合の係数群を得るものであり、前記成分波形抽出部は、それぞれの前記係数集合の係数群に対する成分波形の適合度を算出するものである、<9>または<10>に記載の成分波形抽出システム。
<13>前記解析部は、前記信号波形に含まれる複数の副波形成分を連結させたものを1つの成分波形としてウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて係数集合の係数群を得るものである、<9>または<10>に記載の成分波形抽出システム。
<14>前記適合度が、前記係数集合の係数群に成分波形の基本特性に基づいて得られるものである、<2>に記載の信号解析システム。
<15>前記基本特性が、成分波形の特性を示す基本ベクトルであり、
前記適合度は、前記基本ベクトルと前記係数集合の係数群に基づくベクトルとにより得られるものである<14>に記載の信号解析システム。
<16>前記係数集合の係数群に、成分波形に基づいて重み付けされた適合フィルタが適用された係数集合の係数群に対し、成分波形との適合度を算出する<2>、<14>、および、<15>のいずれかに記載の信号解析システム。
<17>前記適合度と、成分波形の基本特性に基づいて重み付けされた調整フィルタを用いて、前記係数集合の係数群に包含された複数の波形と成分波形を分離し、成分波形に対応したウェーブレット係数群に調整するものである<2>、および、<14>~<16>のいずれかに記載の信号解析システム。
<18>前記調整フィルタは、前記係数集合の係数群を、前記適合度に基づいて検出された複数の成分波形のそれぞれに対応する係数集合の係数群に分配する分配型抽出フィルタを有するものである<17>に記載の信号解析システム。
<19>前記調整フィルタは、前記係数集合の係数群から、前記適合度に基づいて検出された成分波形に対応するウェーブレット係数群を導出する導出型抽出フィルタを適用するものである<17>または<18>に記載の信号解析システム。
<20>前記特徴量分析部により求められる前記成分波形の強度および/または性質を示す特徴量により、前記適合度を補正する、<2>、および、<14>~<19>のいずれかに記載の信号解析システム。
<21>前記強度を示す特徴量は、成分波形の信号の強度を示す信号強度特徴量であり、前記性質を示す特徴量は、成分波形に含まれる周波数成分の分布に基づく周波数特徴量である<20>に記載の信号解析システム。
<22>前記特徴量分析部が、前記適合度に基づいて、前記信号波形に含まれる複数の成分波形を一つのブロックとして得る<20>または<21>に記載の信号解析システム。
<23>前記特徴量分析部は、前記ブロック内の、成分波形の変動平均と変動変化量を分析すると共に、前記変動平均および変動変化量に基づいて、ブロックの変動係数および変動度を解析するブロック変動解析部を有するものである<22>に記載の信号解析システム。
<24>前記ブロック変動解析部は、前記ブロックの変動度に基づいて、ブロックの変動に対する、成分波形の活性度を推定するものである<23>に記載の信号解析システム。
<25>前記特徴量分析部は、前記ブロックの変動度と前記ブロックを構成する成分波形の信号強度特徴量に基づいて、ブロックの活動係数およびブロックの活動度を解析するブロック出力特性解析部を有する<22>~<24>のいずれかに記載の信号解析システム。
<26>前記特徴量分析部は、ブロックの発生度およびブロックの特性に基づいて、ブロックの活動余力度を解析するブロック出力変動特性解析部集散度分析部を有する<22>~<25>のいずれかに記載の信号解析システム。
<27>前記ブロック出力特性解析部は、前記信号波形に含まれる複数の成分波形のうち異なる特性を示す成分波形をブロックとして得て、当該異なる特性に基づく当該ブロックの変動度と当該ブロックを構成する成分波形の信号強度特徴量に基づいて、前記活動度を解析するものである、<25>に記載の信号解析システム。
<28>前記特徴量分析部は、前記信号波形に含まれる複数の成分波形のうち異なる特性を示す成分波形をブロックとして得て、当該異なる特性に基づくブロックの発生度およびブロックの特性に基づいて、前記活動余力度を解析するものである、<25>に記載の信号解析システム。
<29>成分波形抽出システムにより、
信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得るステップと、
前記係数集合の係数群に対する成分波形の適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形に対応する係数群を逆ウェーブレット変換により波形再現することで前記信号波形に含まれる成分波形を得るステップとを有する、
成分波形抽出方法。
<30>信号解析システムにより、
信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得るステップと、
前記係数集合の係数群に対する成分波形との適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形の特性に基づく特徴量を分析するステップとを有する、
信号解析方法。
<31>コンピュータを、
信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得る解析部と、
前記係数集合の係数群に対する成分波形の適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形に対応する係数群を逆ウェーブレット変換により波形再現することで前記信号波形に含まれる成分波形を得る成分波形抽出部とを有する、
成分波形抽出システムとして動作させるプログラム。
<32>コンピュータを、
信号波形のウェーブレット係数群を算出して、ウェーブレット係数集合に基づいて、係数集合の係数群を得る解析部と、
前記係数集合の係数群に対する成分波形との適合度を算出して、前記適合度に基づいて成分波形を検出し、前記成分波形の特性に基づく特徴量を分析する特徴量分析部とを有する、
信号解析システムとして動作させるプログラム。
【発明の効果】
【0021】
本発明によれば、計測した信号波形に含まれ、ターゲットとする成分波形の取得や、この取得した成分波形の各特徴を定量的に解析し評価することができる。
【0022】
本発明によれば、周波数が時間毎に変動する成分要素からなる波形の挙動を、単一のチャンネルによるシンプルな計測系により取得した信号波形により詳細な評価が可能となる。これは、成分波形のリアルタイムでの挙動把握にも適用できる。
【図面の簡単な説明】
【0023】
【
図1】成分波形の特徴の評価を説明するための図である。
【
図2】超球上に投射した基本ベクトルを説明するための図である。
【
図4】本発明の実施の形態にかかる成分波形抽出システムの解析フローを示す図である。
【
図5】本発明の実施の形態にかかる信号波形解析システムの解析フローを示す図である。
【
図6】本発明を用いて行った実験を説明するための図である。
【
図7】本発明を用いて行った実験を説明するための図である。
【
図8】本発明を用いて行った実験を説明するための図である。
【
図9】本発明を用いて行った実験を説明するための図である。
【
図10】本発明を用いて行った実験を説明するための図である。
【
図11】本発明を用いて行った実験を説明するための図である。
【
図12】本発明を用いて行った実験を説明するための図である。
【
図13】本発明を用いて行った実験を説明するための図である。
【
図14】本発明を用いて行った実験を説明するための図である。
【
図15】本発明を用いて行った実験を説明するための図である。
【
図16】本発明を用いて行った実験を説明するための図である。
【
図17】本発明を用いて行った実験を説明するための図である。
【
図18】本発明を用いて行った実験を説明するための図である。
【
図19】本発明を用いて行った実験を説明するための図である。
【
図20】本発明を用いて行った実験を説明するための図である。
【
図21】本発明を用いて行った実験を説明するための図である。
【
図22】本発明を用いて行った実験を説明するための図である。
【
図23】本発明を用いて行った実験を説明するための図である。
【
図24】本発明を用いて行った実験を説明するための図である。
【
図25】本発明を用いて行った実験を説明するための図である。
【
図26】本発明を用いて行った実験を説明するための図である。
【
図27】本発明を用いて行った実験を説明するための図である。
【
図28】本発明を用いて行った実験を説明するための図である。
【
図29】本発明を用いて行った実験を説明するための図である。
【
図30】本発明を用いて行った実験を説明するための図である。
【
図31】本発明を用いて行った実験を説明するための図である。
【
図32】本発明を用いて行った実験を説明するための図である。
【
図33】本発明を用いて行った実験を説明するための図である。
【
図34】本発明を用いて行った実験を説明するための図である。
【
図35】本発明を用いて行った実験を説明するための図である。
【
図36】本発明を用いて行った実験を説明するための図である。
【
図37】本発明を用いて行った実験を説明するための図である。
【
図38】本発明を用いて行った実験を説明するための図である。
【
図39】本発明を用いて行った実験を説明するための図である。
【
図40】本発明を用いて行った実験を説明するための図である。
【
図41】本発明を用いて行った実験を説明するための図である。
【
図42】本発明を用いて行った実験を説明するための図である。
【
図43】本発明を用いて行った実験を説明するための図である。
【発明を実施するための形態】
【0024】
以下に本発明の実施の形態を詳細に説明するが、以下に記載する構成要件の説明は、本発明の実施態様の一例(代表例)であり、本発明はその要旨を変更しない限り、以下の内容に限定されない。なお、本明細書において「~」という表現を用いる場合、その前後の数値を含む表現として用いる。
【0025】
[本発明のシステム、方法]
本発明においては、解析対象の信号波形から成分波形を抽出する成分波形抽出システムと、抽出した成分波形についてその特徴を定量的に分析可能な信号波形解析システムを有する。成分波形抽出システムでは、計測された信号波形に対して冗長ウェーブレット解析を行い、その信号波形の各時刻に存在する波形の特徴を求める。そして、この解析された各波形の特徴に対して、予め設定した、ターゲットとする成分波形の特性に基づき、当該波形が成分波形であるとする、適合度と共に、成分波形を抽出する。また、信号波形解析システムでは、抽出した成分波形に対して、その大きさ等の定量的な値を算出し、信号波形における、この成分波形の挙動に対する定量的な評価を可能とする。
【0026】
[解析対象]
本発明は様々な解析対象波形および成分波形についても適用できる。例えば、黙声認識(いわゆる、口パク)や、心音、心電図による電気的波形、筋肉の活動等、を捉える筋電波形、または工学的には、機器や各種環境から生じる騒音、振動波形、または機器等の故障検知を目的とした波形を対象とすることができる。本発明においては、解析対象とする波形に筋電波形、成分波形に運動単位活動電位波形を典型的な例として説明する。
【0027】
例えば、本発明を筋電波形計測に用いる場合、速筋/遅筋を形成する運動単位の活動に対して、既発明では当該運動単位の周波数分布から速筋/遅筋の時間的変動を可視化し定性的に把握するに留まるが、本発明では、ターゲットとする成分波形の抽出、及び定量的評価が可能となるため、その成分波形の対象となる、速筋/遅筋の活動の定量化ができそれぞれの筋活動の状態や疲労の程度、または筋繊維全体で見た時に、それぞれの筋活動の影響等を評価することができる。
【0028】
先ず、初めに成分波形抽出システムについて説明する。成分波形抽出システムには、各種センサ類により計測された解析対象の信号波形を解析する解析部と、解析部により求められる各種情報に基づき、解析対象の信号波形から成分波形を抽出する成分波形抽出部を有する。それぞれについて以下に説明する。
【0029】
[解析部]
解析部では、各種センサ類により計測し取得された信号波形に対して冗長ウェーブレット解析が行われ、この解析により算出された、信号波形の各時刻に存在する波形のウェーブレット係数からなる係数群(以下、「ウェーブレット係数群」という。)から、抽出しようとする成分波形に対応するウェーブレット係数からなる係数群(以下、「係数集合の係数群」という。)が求められる。
【0030】
冗長ウェーブレット解析では、使用するウェーブレットの形状に近い形状の信号変化に対して敏感に反応するため、ある波形を解析した場合、この波形を構成する様々な周波数成分が存在する時間帯毎に、この周波数成分に対応したウェーブレット係数が出力される。また、このウェーブレット係数を逆変換することによりこの周波数成分の波形が再現される。このように、ウェーブレット係数は解析対象とする波形の周波数成分と強い関係性を有するものであることから、解析部で求められる係数集合の係数群は、信号波形に存在する成分波形の特徴を示すということができる。
【0031】
係数集合の係数群を算出するために、先ずは、計測し取得された信号波形に対して冗長ウェーブレット解析が行われ、信号波形の各時刻に存在する各々の波形のウェーブレット係数群が算出される。そして、次に成分波形に対応するウェーブレット係数集合が設定され、このウェーブレット係数集合に基づいて上記の各々の波形のウェーブレット係数群から、成分波形に対応する係数集合の係数群が抽出される。
【0032】
ここで、成分波形に対応するウェーブレット係数集合について説明する。ウェーブレット係数は、上記で述べたように解析対象とする波形の周波数成分と強い関係性を有するものであり、更にその関係性は当該周波数成分が存在する時間帯とも関連することから、時間-周波数平面上で規則性を持って分布するものが存在する。ウェーブレット係数集合は、そのウェーブレット係数に対し、時間-周波数領域において、時間位置と周波数帯域との関係を規定したものである。その点を踏まえると成分波形に対応するウェーブレット係数集合(以下、「ウェーブレット係数集合」ともいう。)は、成分波形をウェーブレット係数で表した時、時間-周波数領域において、そのウェーブレット係数を時間位置と周波数帯域との関係で規定したものである。
【0033】
またウェーブレット係数集合で捉えようとする成分波形は、そのウェーブレット係数集合と同様に時間幅があり、ウェーブレット係数集合の基準時点は、成分波形が存在する時間幅から選択される時刻でもある。これにより、上記で算出された信号波形の各時刻に存在する各々の波形のウェーブレット係数を、例えば時間-周波数平面上で表した場合の各時刻において、当該時刻を基準時点とする、ウェーブレット係数集合に基づいてウェーブレット係数を抽出すれば、当該時刻に存在する成分波形に対応するウェーブレット係数群を抽出することができるようになる。つまり、ウェーブレット係数集合に対応するウェーブレット係数を全て、各周波数帯域毎に当該基準時点に基づいて調整し配置すれば、前述したように、成分波形に含まれる各周波数成分は、周波数帯域とそのウェーブレット係数の組からなる群、すなわち、係数集合の係数群として、それぞれが異なる時刻に出力されることなく、成分波形が存在する時刻に対応して出力される。
【0034】
ウェーブレット係数集合の定め方についてはどのような手法をとっても問題はないが、例えば、次のようにして定めることができる。
(1)成分波形の理想形をウェーブレット変換した結果に基づいて設定
(2)試験用の信号データサンプルのウェーブレット変換結果から抽出
(3)逆ウェーブレット変換結果が成分波形に近づくように基本的なウェーブレット係数集合から微調整
【0035】
(1)に関しては、実験結果等の各種情報から推定可能な形状、一般的に周知または公用されている形状、または類似する波形に基づき作成したもの等、各種形状に基づく波形を成分波形として用い、冗長ウェーブレット解析によりウェーブレット係数を求め、ウェーブレット係数集合を設定することができる。(2)に関しては、意図的に、ターゲットとする成分波形が出力される試験を行い、この試験で取得し信号波形を冗長ウェーブレット解析して、ウェーブレット係数集合を設定する。
【0036】
(3)に関しては、あるウェーブレット係数集合を微調整し、逆ウェーブレット変換により成分波形に近付けることにより設定する。具体的には、例えば、既発明にあるように、時間―周波数平面において、あるウェーブレット係数集合に対し、構成する各周波数帯域の領域全ての時間軸方向末端位置を左右どちらか一方に揃えたウェーブレット係数集合を設定する。
【0037】
次に、このウェーブレット係数集合に基づき、その最も低い周波数帯域の末端位置を基準とし、隣接する高周波数側の高周波数帯域の末端を調整し、更に調整した帯域と隣接する高周波数側の帯域間でも同様に順次調整してターゲットとする成分波形のウェーブレット係数集合を設定する。このウェーブレット係数集合を設定する際、調整する大きさの微調整により形成されるウェーブレット係数集合に基づきウェーブレット係数群を時間-周波数平面上にプロットした場合に、その係数群が時間軸に垂直に並び、かつ、その係数群を逆変換した結果が成分波形に近くなるように調整することが肝要である。なお、逆変換による信号波形の作成は、係数群の各係数を本来の時間位置に配置した上で、係数ごとに逆変換して得られた波形を足し合わせることで行うことができる。
【0038】
ウェーブレット係数集合のパターンの一例としては、T+S及びH-Sがある。T、Hは、ウェーブレット係数集合の1つのパターンを示すものであり、T及びHは各周波数帯域の領域の、時間軸方向末端(時間軸方向に対する大きい側)及び時間軸方向先端(時間軸方向に対する小さい側)の位置が揃ったウェーブレット係数集合を表す。また、Sは各周波数帯域の領域のシフト量を示し、例えば、T+Sの場合、そのTのパターンを基に、そのパターンの最も低い周波数帯域の時間軸方向末端位置を基準として、その末端に対して隣接する高周波数側の高周波数帯域の時間軸方向末端(時間軸方向に対する大きい側)を当該高周波数帯域の時間幅に対するSのみシフトさせ、その更に高周波数帯域においても、隣接する帯域間で同様に順次、シフトさせて形成される。
【0039】
例えば、T+Sの場合は、T+1/2、T+1/4、T+0、H-Sの場合は、H-1/2、H-1/4、H+0、または、Sがこのような値に限らず、14/32をベースとして、1/32刻みで変化させる等、抽出しようとする成分波形の特徴に応じて、シフト量Sを適宜調整して設定することができる。
【0040】
一方、解析対象信号に、“Daubechies N=2”のウェーブレットで示されるような、1つの主たるピークを持った周波数成分を含む成分波形の場合は、末端位置揃え(T)を標準設定として0以上の調整シフト量を設定いたウェーブレット係数集合が有用であり、シフト量を調整することで立ち上がりと減衰の速度をバランスよく設定し、成分波形の安定した分析が可能となる。
【0041】
図4は、本発明の実施の形態にかかる成分波形抽出システムの解析フローを示す図である。
成分波形のウェーブレット係数集合、及びウェーブレット係数群は、解析部10に記憶部(図示せず)を設定し、この記憶部に保存され、解析部10で解析される際に、当該記憶部から読み込まれ用いられるようにしてもよい。また、ウェーブレット係数群は、解析部10で求められた後、記憶部に保存されずリアルタイムで後述する成分波形抽出部20にて用いることも可能である。またウェーブレット係数集合は、予め設定しておくことも可能であり、この場合、例えば、上記の記憶部に保存しておき、解析時に記憶部から読み込まれるようにしてもよい。
【0042】
上記では解析部10で係数集合の係数群を求める実施例を示したが、解析部10に、係数集合の係数群を演算する演算部(以下、「係数群演算部」という)を設け、解析するようにしてもよい。この場合、計測し取得された信号波形に対して冗長ウェーブレット解析が行われ、信号波形の各時刻に存在する各々の波形のウェーブレット係数が算出される。次に、係数群演算部により、ウェーブレット係数集合に基づき、上記の各々の波形のウェーブレット係数群から、成分波形に対応する係数集合の係数群を抽出し求められる。これにより、解析部10での演算負荷が低減されると共に、計算処理効率を向上させることができる。また、係数群演算部を設けることで、この解析結果に基づく成分波形の挙動の可視化情報を効率的に取得することができると共に、更に、後述する信号波形解析システムと接続しその結果を用いることで、成分波形以外のその周辺の波形の状況を考慮した分析等、有用な分析を行う事が可能となる。
【0043】
信号波形を取得するセンサについては、解析対象とする信号を取得できるものであれば特に限定されるものでない。例えば、生体信号であれば、前出の表面筋電を計測する表面筋電センサの他、筋音センサ、心電図センサ、脈派センサ、血圧センサ等が、騒音や振動であれば、音圧センサ、インテンシティセンサ、加速度センサ、歪センサ等が挙げられる。
【0044】
[成分波形抽出部]
次に成分波形抽出部20について以下に説明する。成分波形抽出部20では、抽出された波形を成分波形とする確からしさを示す適合度を算出すると共に、その適合度に基づき成分波形を抽出する。適合度の算出に関しては、成分波形の基本特性に基づき行われる。
【0045】
[基本特性]
基本特性は、ウェーブレット係数集合に基づく、成分波形の形状や大きさ、信号波形中での時間的変化等、成分波形が持ちうる特徴である。具体的には、成分波形を構成するウェーブレット係数値とその周波数及び時間との関係を示す、ウェーブレット係数の周波数特性で設定され、例えば、成分波形の代表的な形状に対応する周波数特性や、成分波形が時間的あるいは何らかの外的要因等により変化することを考慮しその変化特性の過程における各種波形形状に対応する周波数特性等を設定することができる。またこれら基本特性については、解析目的や精度等に応じ適宜設定することができる。ウェーブレット係数は、先で述べたように解析対象とする波形の周波数成分と強い関係性を有するものであることから、成分波形のウェーブレット係数集合を形成するウェーブレット係数の分布に基づく成分波形を構成する各周波数成分の形状や大きさ、及び時間的変化を把握することができる。この基本特性と対比することにより解析部10で各時刻に対応して求められたウェーブレット係数集合の係数群を成分波形とみなすことの妥当性や、前記係数群が成分波形であるとした場合のウェーブレット係数値の大きさや、ウェーブレット係数の周波数分布の時間的変化等から、信号波形中での成分波形の時間的特徴を把握することができる。
【0046】
適合度の算出に用いる基本特性としては、一例として後述するように、成分波形の特徴をベクトル化したものとすることができる。このベクトルは例えば、成分波形のウェーブレット係数集合に基づき、各周波数帯域(以下、「周波数レベル」という)を次元とし、その周波数レベルに対応するウェーブレット係数をそのベクトルの成分要素として構成することができる。また、成分波形はその波形形状が信号波形中で時間的に変化するため、成分波形の時間的変化の過程における、代表的な状態を示す波形に対応するウェーブレット係数の周波数特性に基づき同様にベクトル化し、これらベクトルを各々、単位ベクトル化したものを基本ベクトルとして設定する。そして、ある時刻の波形に対して解析されたウェーブレット係数群に対しても同様にベクトル化し、基本ベクトルとの対比による成分波形として適合している度合いや、その波形が、成分波形の変化の過程のどの時点において生じたものかを分析・評価することができる。
【0047】
上記の基本ベクトルに関する具体的な例としては、後述する、速筋非疲労時、速筋疲労時、遅筋非疲労時、遅筋疲労時の波形のウェーブレット係数に基づき示すことができる。速筋非疲労時と速筋疲労時、または遅筋非疲労時と遅筋疲労時は、速筋または遅筋を成分波形とした場合、その変化の過程における代表的な状態である。速筋、遅筋という運動単位の活動の波形を、その変化する波形を含め抽出する場合は、これらのウェーブレット係数の周波数特性に基づき、各周波数レベルを次元とする4つの基本ベクトルを設定することができる。このように基本ベクトルは、成分波形の性質に応じて、1つに限らず複数個を設定することが可能である。
【0048】
[適合度及び成分波形の抽出方法]
次に、適合度及び成分波形の抽出方法について説明する。
【0049】
成分波形はその波形形状が信号波形中で時間的に変化するため、抽出の際には成分波形の形状に許容される変化の方向性(例えば、上記の筋疲労前後での変化)を加味する必要がある。成分波形の特徴をベクトル化(簡単のため、2次元で議論することとする)して見た時、その際その特徴は、ある成分波形の変化範囲とするAからBの間で変化しうるとすると、この時、サンプルデータの特徴がAからBのパス上にある時には同等の「成分波形らしさ」を持つと評価されると考えられる。
【0050】
図1は、成分波形の特徴の評価を説明するための図である。例えば、
図1における成分波形らしさはV=W>X>Y>Zと評価すべきである。ただし、パスA、Bは変化過程の一部であり、その延長線上も成分波形の妥当な変化と考えることができる場合にはV=W=Z>X>Yと評価されるべきとも言える。
【0051】
ここで更に、上記の2次元の内容を多次元空間で考え成分波形の特徴をベクトル化して表した場合、このような成分波形の波形形状の特徴はベクトルの方位に現れ、信号の強さはベクトルの長さに現れると考えることができる。
【0052】
そこで、先ずは信号の強弱を加味せずに波形形状の特徴だけを捉えるために、超球上に投射した基本ベクトルを検討する。また基本ベクトルとしては成分波形の特徴をベクトル化したものを単位ベクトル化し半径1の単位超球上に投射したものとするのが簡単である。なお、各次元の値は0以上とする。
図2は超球において原点0、超球上の点A、Bを通る平面を切り出したものとする。
【0053】
また、A、Bを基本ベクトルとし、信号波形の解析により「成分波形らしさ」の評価を行う波形の、ある特徴を示すベクトルをCとし、孤AB上の点Α、Βについて基本ベクトルA、Bとの関係を検討すると、この点Α、Βのベクトルと基本ベクトルとのなす角は∠AΑB=∠AΒBであり、超球上(ただし、各次元の値は0以上)の任意の点Cに関する∠ACBで最大の値となる。また、ある点Cを設定し点CがAB間にあることを∠ABC≦90°かつ∠BAC≦90°のときと定義すると、点Cが孤ABから遠ざかるほど点Cのベクトルと基本ベクトルA、Bとのなす角∠ACBは小さくなる。
【0054】
このように、成分波形の波形形状の特徴はベクトルの方位に現れ、また評価を行う波形とこの特徴との関係においては、基本ベクトルとの距離により検討することができる。また、この距離による検討は、評価を行う波形の特徴を示すベクトルと基本ベクトルとのなす角に基づき行うことができる。なお、
図3に示すように、点γのようにabのb側の外にある場合には、∠aγbの代わりに∠abγで評価する。
【0055】
本発明において、解析対象の信号波形のある時刻に存在する波形の特徴は、その時刻に対応するウェーブレット係数群の各周波数レベルのウェーブレット係数値によって表される。成分波形抽出部20では、各周波数レベルを多次元空間の座標軸として、成分波形のウェーブレット係数集合に基づく基本ベクトルと、解析部10により解析された、上記信号波形中の成分波形に対応する係数集合の係数群に基づくベクトルを設定し、基本ベクトルとの関係による、上記の考え方に基づいて「成分波形らしさ」を評価することができる。
【0056】
つまり、上記の考え方に基づく「成分波形らしさ」を、ある時刻の係数集合の係数群が成分波形を成しているかどうかとうかの「適合度」という評価値として算出することができる。なおこの考え方において、許容範囲(
図3のa、b)に対してどのような距離関係にあるかを調べるために、ベクトルのなす角度に基づいて評価できることを示したが、許容範囲の各端点からの距離(ユークリッド距離や開き角など、妥当なものを適宜選択すれば良い)および距離の比で評価しても良い。
【0057】
適合度が算出されると、その適合度が大きくなる時刻が成分波形の存在時刻として推定される時刻であるために、その時刻のベクトルWに対応するウェーブレット係数群を逆ウェーブレット変換することにより、成分波形が抽出される。
【0058】
成分波形は時間的に幅を持って存在するため、その成分波形の存在時刻として規定される時刻から少しずれた時刻であっても、成分波形として値の大きい適合度を持つと予想される。この点も考慮し、成分波形を捉えるために、適合度の値に限らず、適合度の変化、または適合度の値と適合度の変化との両方に基づくことで、最適な成分波形の存在時刻と形状(逆ウェーブレット変換で獲得)とを得ることができる。
【0059】
なお、上記の評価による成分波形らしさは形状的な特徴のみに基づいているが、例えば、信号の強度(大きさ)など、その他の特徴に対しても形状的な特徴を含めて評価することも可能である。この場合は、後述する、適合強度のように、その他の特徴パラメータとして適合度と掛け合わせることで評価することができ、より解析目的に沿って高精度に成分波形の抽出が可能となる。
【0060】
また、成分波形の形状が固定的であり、観測される波形の散らばりが確率的である場合、サンプリングした信号波形中の成分波形に対する「成分波形らしさ」は、上記の基本ベクトルとのユークリッド距離やマハラノビス距離などの単純な距離による評価で評価することも可能である。
【0061】
[適合度演算部]
適合度の算出に関しては、成分波形抽出部20とは別に、適合度演算部(図示せず)を設けて行うことも可能である。具体的には、適合度演算部において、成分波形のウェーブレット係数集合と、解析部10で解析された係数集合の係数群に基づき、上記と同様に、基本ベクトルを設定し、この基本ベクトルに基づき適合度を算出する。適合度演算部を設定することにより、成分波形抽出部の計算負荷の軽減を図ることができると共に、例えば、前述した、「成分波形らしさ」が形状的な特徴以外にも基づく場合、適合度を形状的な特徴によるものとして区別して取り扱うことができ、成分波形の評価を効率的に行うことができる。
【0062】
[フィルタ]
適合度の算出/評価、及び成分波形の抽出を、より効率的に行うために、適合フィルタ、調整フィルタを設定し用いることができる。適合フィルタに関しては、解析部10により解析された係数集合の係数群に対して、調整フィルタは、ある波形の適合度が算出された後の係数集合の係数群に対して、その適合度と共に用いることが可能である。以下に両フィルタについて説明する。
【0063】
[適合フィルタ(抽出フィルタ)]
成分波形の周波数成分は、信号波形の解析により得られるすべての周波数帯域に存在しているとは限らない。また、成分波形が存在する周波数帯域でも、対象とする以外の信号波形が混在する帯域やS/N比が低い帯域などが存在する可能性がある場合、成分波形に対する適合度の評価、及び抽出に対してはこれらの影響に関し考慮する必要がある。成分波形の基本ベクトルを定める際にも、必要とするのはその成分波形に含まれる周波数成分の帯域のみである。
【0064】
そこで、適合度の評価、及び成分波形の抽出を行う各周波数帯域を適正なものとするために、各周波数帯域の各々に重みを規定する適合フィルタを設定することができる。適合フィルタは周波数帯域ごとに0~1の重みを与えるものであり、成分波形の周波数成分を含まない帯域には0を、対象外の信号成分の混入が予想されるために影響力を抑制したい帯域にはその影響に応じ0~1の範囲の値(例えば、中間の値等)を設定する。適合フィルタは、適合度を計算する前に適合度評価の対象となる係数集合の係数群に掛け込む形で使用する。
【0065】
適合フィルタの適用例としては、例えば、前述したように成分波形を速筋及び遅筋を対象とした場合、比較的高周波数帯域に留まる速筋の成分波形と、低い周波数帯域まで存在する遅筋の成分波形とが混在する筋電波形信号において、速筋の成分波形を抽出するため適合度を計算する際には低周波数帯域の影響を抑制するように適合フィルタを設定することができる。また、同様に、遅筋成分のためには遅筋成分を抽出するための適合フィルタも設定することが可能である。
【0066】
成分波形の抽出や分析のために用いる係数集合の係数群としては、この適合フィルタを適用した後の値を使用することができる。この場合、例えば、筋電波形信号の上記速筋成分と遅筋成分のように、複数種類の成分波形の分析や抽出を行う場合にはそれぞれに異なる適合フィルタが設定されるため、元となる係数集合の係数群は同じであっても、それぞれの解析で使用する係数群は異なる。
【0067】
特に、速筋の成分は比較的、高周波域に留まることなどを鑑み、基本ベクトルと係数集合の係数群を比較する際に当該成分が存在する範囲外の周波数帯域のウェーブレット係数値を削減するための適合フィルタ(Мatching Filter)を設定することができる。この際、例えば、成分波形のウェーブレット係数集合の要素ごとに1または0の値で重み設定しても良いし、実数で傾斜を付けるような重みにしてもよい。また速筋対象と遅筋対象とのそれぞれで同等の処理を行うことも可能である。
【0068】
筋電波形の場合に限らず、周波数帯域に違いがある複数種類の成分波形が混在する場合(単一種類の成分波形を周波数帯域で分けて分析しようとしている場合を含む)にも、上記の速筋と遅筋の関係のように扱うことが可能である。それらの発生の同時性は、例えば各適合度の積を利用する等、後述の適合度評価値に基づいても分析できる。
【0069】
[調整フィルタ]
解析対象信号が複数種類の成分波形から構成されている時、それらを分離して抽出するには、各周波数帯域のウェーブレット係数値をそれぞれの成分波形の要素として振り分ける必要がある。そのために成分波形の種類ごとに調整フィルタを設定することができる。調整フィルタは、成分波形のウェーブレット係数集合の各周波数帯域ごとに重みを考慮して設定したものであり、分配型調整フィルタと導出型調整フィルタを有する。
【0070】
分配型調整フィルタは、複数種類の成分波形がある周波数帯域で重なりを生じている等、その周波数帯域のウェーブレット係数値がそれらの成分波形のウェーブレット係数値の合成となっている場合に、そのウェーブレット係数値を、それぞれの成分波形に振り分けるために使用される。また、導出型調整フィルタは、適合度が算出された係数集合の係数群に対して、適合フィルタにより削除された範囲等、適合度の評価に使用されなかった範囲に成分波形の周波数分布が反映されるように補完するために使用される。
【0071】
また、分配型調整フィルタは、ウェーブレット係数集合の各周波数帯域ごとに重み値が設定され、適合度と組み合わせて使用される。具体的には、各周波数帯域において、各成分波形ごとに適合度と分配型調整フィルタの値の積を求め、それらの比によってその周波数帯域のウェーブレット係数値が各成分波形に振り分けられる。そして、成分波形の分析を行う際には、これによって分配されたウェーブレット係数値が用いられる。
【0072】
適合度が算出された係数集合の係数群をWとし、Wのレベル-iのウェーブレット係数をw-iとする。また、レベル-iにおける各成分波形の分配型調整フィルタの値をDF-iとして設定する。
次に、分配対象とする複数の成分波形に対して算出された適合度とそれぞれの成分波形の分配型調整フィルタとの積を算出し、その積の比により、Wに含まれるそれぞれの成分波形に対してウェーブレット係数を振り分ける。
【0073】
例えば、筋電波形解析を例として、速筋用分配型フィルタfDF-i、遅筋用分配型フィルタの値をsDF-iとする。また、各々速筋運動単位活動電位波形としての適合度がfMR、遅筋運動単位活動電位波形としての適合度がsMRであるとき、速筋と遅筋のそれぞれの係数集合の係数群fWとsWの、レベル-iのウェーブレット係数値fw-iとsw-iは次式で与えることができる。
fw-i=w-i*(fMR*fDF-i/(fMR*fDF-i+sMR*sDF-i))
sw-i=w-i*(sMR*sDF-i/(fMR*fDF-i+sMR*sDF-i))
【0074】
以降、単に、係数集合の係数群Wと記述した場合、特に明示されていなければ、速筋運動単位活動電位波形を対象としている時はfWを、遅筋運動単位活動電位波形を対象としている場合はsWを指すこととする。
【0075】
導出型調整フィルタについて説明する。解析対象の信号波形において、例えば、成分波形の周波数帯域の一部の帯域に成分波形とは異なる信号の混入や、波形が複雑に重畳する等、解析精度に影響する信号又は状態を有することが想定される。本発明では、このような解析精度へ影響する可能性のある周波数帯域に対しては適合フィルタを用いて削除して解析することができるが、一方、このような周波数帯域にも本来、解析対象とする成分波形の要素が含まれていた可能性があり、より解析精度を高めるためには削除した周波数帯域の情報を反映するのがよい場合がある。そこで、導出型調整フィルタは、適合フィルタが適用され適合度が算出された係数集合(分配型調整フィルタの適用によって得られた係数集合も含む)の係数群の結果を利用して、適合フィルタで削除されたウェーブレット係数値を推定し補完するものである。
【0076】
導出型調整フィルタは、ウェーブレット係数集合に基づき、解析目的や確保すべき精度等に応じて、変化の過程を含め成分波形が持ちうるウェーブレット係数値の周波数分布を検討、反映したものであり、具体的には、そのウェーブレット係数値の周波数特性に基づいてその周波数帯域毎に重み付けし、成分波形が持つと推定されるウェーブレット係数値として設定したものである。導出型調整フィルタによる補完方法に関しては、導出の基準となる成分波形としての強度を評価するために、適合フィルタが適用された係数集合の係数群のウェーブレット係数値(以下、「適合フィルタ適用係数値」という。)が算出される。次に、導出型調整フィルタにも同様に適合フィルタを適用し、その適用した結果を係数集合の係数群と見做して、値の大きさを評価するためのウェーブレット係数値(以下、「導出型調整フィルタ係数値」という。)を求める。そして、適合フィルタ適用係数値と導出型調整フィルタ係数値との比と適合度を導出型調整フィルタのウェーブレット係数値にかけた値により補完される。つまり、導出型調整フィルタ係数値は、導出型調整フィルタの、適合フィルタが高い値を持つ範囲のウェーブレット係数値となることから、その値との比率で算出することで、導出型調整フィルタのウェーブレット係数値に対する信号波形中に存在しうる値の割合を求めながら、更に、その波形が実際に存在しうる程度を示す適合度を乗算することで、成分波形が本来持ちうると推定されるウェーブレット係数値が導出される。
【0077】
具体的には、解析部10により算出された係数集合の係数群におけるウェーブレット係数値をw1、...、wnとし,適合フィルタをMF1、...、MFn(ただし0.0≦MFi≦1.0)とする。この適合フィルタを適用し得られた適合度をMRとし、導出型抽出フィルタをEF1、...、EFnとする。EFiには、必要に応じて負の値を設定してもよい。
【0078】
先ず、適合フィルタを適用後のウェーブレット係数の二乗和の平方根をSwとして、次式を用いて算出する。
Sw=sqrt(Σ((wi*MFi)2))
wi*MFiは、適合フィルタ適用係数値であり、Swは適合度評価に使った部分の係数値の大きさに相当する。
また、導出型調整フィルタのウェーブレット係数値に適合フィルタを適用したものの二乗和の平方根をSEとして、次式を用いて算出する。
SE=sqrt(Σ((EFi*MFi)2))
EFi*MFiは、導出型調整フィルタ係数値であり、SEは、適合度評価に使われている範囲での導出型調整フィルタの大きさに相当する。導出型抽出フィルタは成分波形のウェーブレット係数分布を反映させて設定するので、導出型抽出フィルタの値の大きさの成分波形が存在したとした場合に、適合度評価に使われるであろう部分の係数値の大きさに相当する。
【0079】
この時,抽出する成分波形のウェーブレット係数値をew1,...,ewnとすると、ew1,...,ewnは、MR及び上記Sw、SEを用いて次式により算出される。
ewi=(wi*MFi)+sqrt(MR*Sw/SE)*EFi*(1-MFi)
本式により、例えば、適合フィルタが1、つまり適合度評価においてそのままの値を用いた部分には、係数集合の係数群の係数値がそのまま使用され、適合フィルタが0、つまり適合度評価において削除して使用しなかった部分は、成分波形の周波数特性を反映した導出型調整フィルタの値に基づいて推定される。
【0080】
なお、適合フィルタの値はその部分の係数値の信頼度のようなもので,解析への影響を抑制するため、適合フィルタにより係数集合の係数群におけるウェーブレット係数に対して適合度評価から除いたり、その値を小さくしたりすることができる。上式ewiは、その適合フィルタが適用されたそれぞれの値の比率に応じて、wiの直接の値に、成分波形の周波数特性からの推定値を重み付き加算し設計される。また、MRは,適合度評価に使った係数値をどの程度そのままに成分波形の値として信じて良いかを反映しているものとして解釈することもできるので,補充分の係数値をどの程度信じて加算するかの量の重み付けの一部に使用される。
【0081】
次に、
図4の解析フローに基づき、成分波形抽出システムにおける、成分波形を抽出する処理の流れについて説明する。先ず、解析部10により、各種センサにより計測、取得された解析対象信号波形を冗長ウェーブレット解析しウェーブレット係数群が算出され、ウェーブレット係数集合に基づき、このウェーブレット係数群から係数集合の係数群が抽出される。次に、成分波形抽出部20により、成分波形の基本特性に基づき、解析部10により抽出された係数集合の係数群に含まれる波形に対する適合度が算出される。その際に、所定の基準に基づき、成分波形の適合度とその適合度を持つ係数集合の係数群から成分波形に対応するウェーブレット係数群が抽出され、そのウェーブレット係数群におけるウェーブレット係数を逆ウェーブレット変換することにより成分波形の波形が算出、推定される。
【0082】
図4に示す適合フィルタに関しては、成分波形の適合度と成分波形に対応するウェーブレット係数群が抽出される前に、適宜、適合フィルタを適用し、係数集合の係数群から解析に不必要な周波数帯域に対してフィルターをかけて、成分波形に対応するウェーブレット係数群を効率的に抽出できるように用いることが可能である。
【0083】
調整フィルタに関しては、抽出された成分波形に対応するウェーブレット係数群に対して、成分波形の基本特性に基づき、補完するようにフィルタをかけて、抽出される波形が、より成分波形に近付くように用いることが可能である。
【0084】
以下において、上記の解析フローによる、基本ベクトルによる適合度の算出及び成分波形の抽出処理について、筋電波形を一実施例として具体的に説明する。基本ベクトルの元となる、成分波形のウェーブレット係数集合、及びウェーブレット係数群の各レベル(周波数帯域)の値はすべて正値とする。なお、成分波形の特性上、特定のレベルの値が負値であることが適切である場合には、そのレベルの正負の扱いを逆にする(例えば、ウェーブレット係数値の正負を反転させたり、半波整流時に負値のみとする等)ことで、正値と仮定した場合と同等に扱うことも可能である。
【0085】
また、ここでは成分波形の変化の方向性として特定のレベルの値が正と負の間で変化することは想定していないが、オフセット値の加算(および必要であれば定数倍)による可逆加工(ウェーブレット逆変換時には逆加工で戻した上で用いる)を施すことで本手段を適用することも可能である。
【0086】
ある時刻で取得した信号波形をサンプルとして冗長ウェーブレット解析したウェーブレット係数集合に基づき得られた係数集合の係数群を半波整流して正値のみとし、ベクトルをWとする。ただし、ここでWはゼロベクトルではないものとする。Wがゼロベクトルである場合については後述する。さらに、Wを単位ベクトル化したものをSとする。
【0087】
またこの時、低疲労時の基本ベクトルをL、高疲労時の基本ベクトルをHとすると、各ベクトルの終点は半径1の超球上に配されることになる。以降、曖昧でない限りは各基本ベクトルの終点を同じ記号で示すこととする。
【0088】
終点をS、L、Hとする単位ベクトルが成す角をθSL、θSH、θLHとする。このとき、ベクトル終点を結ぶ直線の長さは
|SL|=sqrt(2*(1-cos(θSL)))
|SH|=sqrt(2*(1-cos(θSH)))
|LH|=sqrt(2*(1-cos(θLH)))
である。各cos値は単位ベクトルの内積で得られる。
【0089】
∠LSH=φS、∠SLH=φL、∠SHL=φH、とするとき、以下の通りとなる。
cos(φS)=(1+cos(θLH)-cos(θSL)-cos(θSH))/(2*sqrt((1-cos(θSL))*(1-cos(θSH)))
cos(φL)=(1+cos(θSH)-cos(θSL)-cos(θLH))/(2*sqrt((1-cos(θSL))*(1-cos(θLH)))
cos(φH)=(1+cos(θSL)-cos(θSH)-cos(θLH))/(2*sqrt((1-cos(θSL))*(1-cos(θLH)))
なお、各分母が0になることは、後述する条件を設定することで除外されるため、問題になることはない。
【0090】
係数集合の係数群が速筋(あるいは遅筋)の運動単位活動である可能性が高ければ、Sは孤LH上か、それに近い位置に存在すると推測できる。φSは、Sが孤LH上にある時に最大となり、離れるほどに小さくなり、またこれに伴ってcos(φS)が大きくなる。これらの値を用いて運動単位活動としての適合性の高さの評価に利用することができる。なお、上記において、cos(φS)が最小の場合はcos(φS)=-sqrt((1+cos(θLH))/2)となる。
【0091】
更に孤LHの間ではなく延長線上にある場合も適合性は高いと考えることができることから、Sの適合度(Matching Rate)をMRSとし、MRS(ただし0≦MRS≦1)は次のように与えることができ、適合度が求められる。
【0092】
(1)1-cos(θSL)<ε(ただし、εは適当な小さな値)の場合
SとLとは同一と見なし、|SL|=1、|SH|=0、MRS=1.0とする。
【0093】
(2)1-cos(θSH)<ε(ただし、εは適当な小さな値)の場合
SとHとは同一と見なし、|SL|=0、|SH|=1、MRS=1.0とする。
【0094】
(3)cos(φL)≧0 かつ cos(φH)≧0の場合
SはLとHの間に存在する状況である。
mincos=-sqrt((1+cos(θLH))/2)
MRS=(1-cos(φS))/(1-mincos)
【0095】
(4)cos(φL)<0の場合
Sは弧LHのL側での外に存在する状況である。
mincos=-sqrt((1+cos(θSH))/2)
MRS=(1-cos(φL))/(1-mincos)
【0096】
(5)cos(φH)<0の場合
Sは弧LHのH側での外に存在する状況である。
mincos=-sqrt((1+cos(θSL)))/2)
MRS=(1-cos(φH))/(1-mincos)
【0097】
なお、L、Hが許容できる変化の端点として捉えることができる場合には、上記(4)、(5)の適合度は上記(3)の場合よりも低く評価する方が適切である。その場合には上記(4)、(5)でのmincosとして上記(3)のmincosと同じものを使用する。また、これまでWはゼロベクトルではないとの前提で適合度MRを定義したが、Wがゼロベクトルである場合にはMRS=0とする。
【0098】
次に、MRSが算出されると、その値が大きくなる時刻が、運動単位活動の存在時刻として推定される時刻であり、その時刻のベクトルWに対応するウェーブレット係数群を逆ウェーブレット変換し、運動単位活動電位波形が抽出される。
【0099】
上述の実施形態は本発明を実施する際の好適な形態の一例であるがこれに限るものではなく、本発明の要旨を逸脱しない範囲において種々変更、追加して実施可能である。以下にその変更・追加例を示す。
【0100】
上述の適合度評価において、成分波形の挙動が、その変化範囲において、どのような状況であるのか、またはどのような方向に変化しているのか等、解析の目的に沿った評価を行うことも可能である。これについては、成分波形の変化特性に基づき、その中で注目する特定に対応する基本ベクトルを設定し、評価対象とする波形のベクトルとの距離に基づき、その基本ベクトルに対応した適合度で評価する。例えば、上述の筋電波形解析を例に説明すると、|SL|と|SH|を用いて、Sが低疲労時と高疲労時とのいずれの側に寄っているかを評価するため、適合度MRSを低疲労側MRLSと高疲労側MRHSに、以下のように振り分けて行う。疲労に関連した何らかの評価での活用を目指すこともできる。
MRLS=MRS*|SL|2/(|SL|2+|SH|2)
MRHS=MRS*|SH|2/(|SL|2+|SH|2)
【0101】
適合度の評価に関しては、例えば、事前にテスト用の信号波形において適合度または適合強度の変化を検討する等で、閾値を設定することにより行うことも可能である。これにより、解析の目的に応じて効率的に解析を行うことができる。
【0102】
成分波形の変化特性が連続する複数の基本ベクトルV1、・・・、Vmax(ただしi=1、・・・、max-1について基本ベクトルペアViとVi+1が成す角はすべて等しい)で規定されている場合、評価対象の波形に基づくベクトルXの適合度算定に使用する基本ベクトルペア(Vi、Vi+1)を次の優先順で選択する。
【0103】
(i)上記(1)または(2)を満たす基本ベクトルペアが存在するとき、そのペアを選択する。
連続する基本ベクトル群の両端を除き、条件を満足するペアは2個存在するが、いずれを選択しても実質的な結果は同じである。
【0104】
(ii)条件(i)を満たす基本ベクトルペアが存在しない場合、上記(3)を満たすペアの内で適合度が最大となる基本ベクトルを選択する。
【0105】
(iii)条件(i)および(ii)を満たす基本ベクトルペアが存在しない場合(上記(4)または(5))、mincosが最小となる基本ベクトルペアを選択する。
【0106】
[適合強度に基づく成分波形の抽出]
適合度での評価に関しては、前述したように、信号の強弱を加味せずに波形形状の特徴のみで捉えることを考慮したが、成分波形の抽出に対し強度も加味する必要がある場合が考えられる。そこで、このような成分波形の抽出を行う際には、係数集合の係数群から、成分波形の基本特性に基づき、成分波形らしさを示す適合度とその適合度を持つウェーブレット係数群を抽出し、そのウェーブレット係数群に対応する波形の強度を示す特徴情報と適合度を乗算し、適合度を補正したものを用いて行うことができる。ここで言う強度とは、単純に成分波形の信号強度が大きいということだけではなく、成分波形が強い、あるいは単純に信号波形が存在する場合に値が大きくなるというような性質を持つ特徴量全般を含むものである。具体的には係数集合の係数群の二乗平均平方根などの信号強度特徴量や、強い成分波形ではその成分波形の各周波数成分の大きさの分布に片寄りが生じるというような性質を特徴量として示す周波数特徴量がある。なお、信号強度特徴量及び周波数特徴量の詳細については後述で説明する。
【0107】
また、上記では係数集合の係数群から適合強度を算出することが記載されているが、本システムでは最終的に成分波形として特定することが重要であるため、適合強度を算出する場合としては、その他として、適合フィルタ適用後のウェーブレット係数群や、調整フィルタ適用後のウェーブレット係数群等も必要に応じて適宜用いられる。対象の特性によっては、それら特徴情報を掛けずに適合度をそのまま用いても良い。いずれを選択しても以降の手順に違いがあるわけではないので、以降、簡単のため、適合度と上記の各々の特徴量を掛け合わせたものを特別に区別をする必要がない限りは「適合強度」という表現を用いる。ウェーブレット係数群が成分波形に近いほど、適合強度は大きくなる。
【0108】
解析対象信号に対して冗長ウェーブレット解析を行い、その結果から各時刻のウェーブレット係数群を得る。そして、得られたウェーブレット係数群に基づき、各時刻の適合強度を求め、その適合強度に基づき成分波形の存在時刻、ウェーブレット係数群を抽出し逆変換した結果を成分波形として獲得する。
【0109】
このように本発明では、計算コストの増大を避けたシンプルな処理で、正常変化の方向性を持つ成分波形の検出・抽出のためのウェーブレット係数群の適合度評価手法を提供することができる。更に本発明は計算効率を大幅に向上させた処理であるため、リアルタイムでの処理が可能である。また適合度が得られれば、その適合度に基づき、成分波形の存在時刻を推定でき、その時刻で抽出したウェーブレット係数群を逆変換することで成分波形を得ることができる。
【0110】
次に、本発明における信号波形解析システムについて説明する。
図5は、本発明の実施の形態にかかる信号波形解析システムの解析フローを示す図である。
信号波形解析システムには、成分波形抽出システムの解析部10で得られる情報に基づき、解析対象信号波形における成分波形の特徴情報を分析する特徴情報分析部30を有する。以下に特徴情報分析部30について説明する。なお、信号波形解析システムにおいて、成分波形抽出システムとは切り離して単独で作動することができるように、成分波形抽出システムの解析部10を信号波形解析システムに設けることも可能である。
【0111】
特徴情報分析部30では、成分波形抽出システムの解析部10で得られる係数集合の係数群に対し、適合度が算出され、この適合度に基づき成分波形に対応するウェーブレット係数群と適合度が抽出され、このウェーブレット係数群に基づき解析対象の信号波形における成分波形の特徴情報を分析する。またこの特徴情報の分析に関し、成分波形に対応するウェーブレット係数群が適合強度に基づき抽出される場合は、適合度に基づき係数集合の係数群から一旦、成分波形らしさを有するウェーブレット係数群とその適合度が抽出され、そのウェーブレット係数群の特徴情報(特徴情報については後述による)と適合度を乗算した値に基づき行うことができる。この特徴情報は、成分波形の強度および/または性質を示すものであり、一例として、成分波形の強度を示すものとして信号強度特徴量と、成分波形の性質を示すものとしては周波数特徴量を有する。また、それぞれの特徴量は、図示しないが強度特徴量分析部、周波数特徴量分析部により分析が行われる。
【0112】
また、特徴情報分析部30は成分波形抽出システムにも設けることが可能である。この場合、特徴情報分析部30では信号強度特徴量または周波数特徴量及び、それらの特徴量と適合度により適合強度が求められ、特徴量を加味した成分波形の抽出を行うことが可能となる。
【0113】
特徴情報分析部30に設けられる強度特徴分析部により、信号強度特徴量を算出することができる。信号強度特徴量は、抽出しようとする成分波形の強度(大きさ)を示すものであり、具体的には、適合度に基づき推定された、成分波形に対応するウェーブレット係数群の強度(WSP;Wavelet-Set Power)として表され、このウェーブレット係数群Wを構成する各周波数レベルkのウェーブレット係数をwk(kは-1~-LMAX)とするとき,WSP(W)は、次式として示される。
WSP(W) = sqrt(Σk(wk)2)
【0114】
信号強度特徴量は、適合度と乗算した値によっても成分波形の分析を行うことができる。この場合においては特に、成分波形を、周波数特性等の波形形状の特徴だけでなくその信号の強弱を加味して分析する場合に有用である。先にも述べた通り、適合度は成分波形の波形形状の特徴に基づき算出される値であるので、適合度をMRとしこの適合度に信号の強弱を表すWSPを乗算した値、つまり、MR*WSPにより信号強度特徴量として用いることができる。
【0115】
なお、MR*WSP、つまり適合度と信号強度特徴量とを乗算した値は、先に述べたように、適合強度として、信号強度も考慮して成分波形を抽出することができる。
【0116】
筋電波形信号解析においては、WSPに適合度MRを乗算したMR*WSPを運動単位活動の観点での信号強度特徴量とすることができる。この場合、MR*WSPは、時間解像度的に分離できない運動単位の活動電位の総和を反映する特徴量であり、強い筋活動を達成するために同期的に多数の運動単位が活動するほどに大きな値となる傾向を持つ。
【0117】
特徴情報分析部30に設けられる周波数特徴量分析部により、周波数特徴量を算出することができる。周波数特徴量として、適合度に基づき推定された、成分波形に対応するウェーブレット係数群Wに対してウェーブレット重心(CoB)を求める。ウェーブレット重心とは、ウェーブレット係数で集合化された係数群において、各ウェーブレット係数の、ある基準の周波数に対するモーメント量の合算を各ウェーブレット係数の合算値で割った値であり、この係数群の中心となる周波数位置を示したものである。具体的には、成分波形に対応するウェーブレット係数群Wのレベル範囲を-1、・・・、-Lmax、各レベルのウェーブレット係数をw-1、・・・、w-Lmaxとするとき、CoB(W)は、以下で与える。
(1)Σk|w-k|=0のとき、CoB(W)=0.0
(2)Σk|w-k|!=0のとき、CoB(W)=(Σk(|w-k|*2(Lmax-k))/Σk|w-k|)-1
【0118】
また、CoBは窓フーリエ変換が扱う、解析の時間幅(窓幅)よりも時間応答の短い、瞬間の周波数特性を示すものであるため、値の時間変化は著しいものになる場合がある。そこで、例えば、中~低周波数帯域での変動に比べて高周波数帯域での変動が目立ちやすく、広いダイナミックレンジでの取り扱いが必要となる等、に対応するため、取り扱う周波数帯域を対数軸で表す特性値で評価することも可能である。例えば、一例として後述するように、対数ウェーブレット重心(LCoB)を用いることができる。
【0119】
その際、例えば、取り得る値の範囲(ウェーブレット係数集合のレベル範囲で規定される)が0~1となるように正規化すると、例えば周波数変動範囲の特性が異なる複数種の信号を同じスケールで比較(指定した周波数範囲内のどのような割合位置にあるか等)する場合等に、好適である。
【0120】
そこで、レベル-1から-Lmaxのウェーブレット係数が得られている時に、その部分区間であるウェーブレット係数Ws=w-a、・・・、w-b(すなわちレベル-aから-b;1≦a<Lmax、a≦b≦Lmax)で定義され、値の範囲が0から1となるような対数ウェーブレット重心(LCoB(Ws))を、以下で与える(ただしΣkにおけるkの範囲はa≦k≦b)。
(1)Σk|w-k|=0のときLCoB(Ws)=0.0
(2)Σk|w-k|!=0のときLCoB(Ws)=log2(Σk(|w-k|*2(b-k))/Σk|w-k|)/(b-a)
なお、CoB(Ws)とLCoB(Ws)の関係はLCoB(Ws)=log2(CoB(Ws)+1)/(b-a)である。
【0121】
CoBもLCoBも、それが示すものの本質が違うわけではないため、以降、特に必要がない限りは両者を区別することなくCoBと呼称する。利用者の必要性に応じてCoBとLCoBとを読み替えればよい。
【0122】
図5に、信号波形解析システムでの適合強度による基本的な解析フローを示す。
まず、解析部10で解析対象信号を冗長ウェーブレット変換し、ウェーブレット係数群を算出する。そして、当該ウェーブレット係数群にウェーブレット係数集合を適用し、係数集合の係数群を得る。次に、成分波形抽出部20で、適合フィルタを用いて、解析部10で得られた係数集合の係数群に成分波形の基本特性に基づいて得られる適合度を適用する(
図5の丸数字1)。そして、調整フィルタを用いて、成分波形に対応するウェーブレット係数群を得る(
図5の丸数字2)。さらに、適合強度による場合、大きさを考慮した成分波形に対応するウェーブレット群を得て、成分波形の特徴情報(成分波形の特性に基づく特徴量)を分析する(
図5のケースB)。なお、適合強度によらない場合、調整フィルタを用いて得られた成分波形に対応するウェーブレット係数群により、成分波形の特徴情報を分析する(
図5のケースA)。
【0123】
[ブロックによる分析]
解析対象とする信号波形中の成分波形の挙動を分析する場合、連続して発生した複数の成分波形の時間間隔が成分波形の時間幅よりも短いために複数の成分波形が重なりながら存在したり、または例えば、振動や騒音分析のオクターブバンドによる解析のように、複数の成分波形が存在する中、あるまとまりとして解析した方がよい場合等も想定される。そこで、特徴情報分析部30では、ブロック分析部(図示せず)が設けられ、解析対象の信号波形中に複数の成分波形が含まれる場合、それらを分離せず成分波形の塊、つまり、ブロックとしてとらえ、ブロックに含まれる成分波形の平均特性として分析することも可能である。ブロック分析部では、所定の範囲内の適合度または適合強度に基づきブロックを設定する。例えば、所定の範囲として閾値を設け、適合度または適合強度が閾値以上になる波形を含むものをブロックとして設定することも可能である。
【0124】
例えば、閾値をΔAとし、適合度MRまたは適合強度がこの閾値を越えている連続区間を成分波形のブロックとして抽出する。適合度または適合強度の連続した値をs0、s1、・・・、sN、sN+1とするとs0とs1の間で閾値を越え、sNとsN+1の間で閾値を下回った場合、つまり連続で値が閾値を越えた時の、その範囲をブロックとして抽出し、そのブロック幅はN+1と規定する。なお、適合強度に基づく場合は、その適合強度は、例えばMR*WSPやMR*CoBで与えられる。
【0125】
また、成分波形が存在する時区間における適合度または適合強度の変化は、成分波形の存在時刻ではピーク波形の形状をなし、ブロック内に複数の成分波形が包含されている場合はブロック内に複数の山として現れる。そこで、このようなピーク波形として現れる適合度または適合強度の変化数によって、このブロック内に含まれる成分波形の数(以下、「成分波形推定数」という)を推定することも可能である。
【0126】
ブロック内の適合度または適合強度の値の間の変化量d0、d1、・・・、dN(ただし、di=si+1-si)を求める。ブロックの先頭と末尾では閾値を横切るため、d0>0、dN<0である。この変化量の正負反転数により、ブロックに含まれる成分波形の推定数を((反転数+1)/2)の小数点以下の端数を切り捨てて整数化した値で与えることができる。
なお、変化量の値が細かい増減変動を繰り返しながら推移し、単純に正負変化数を数えると反転数は非常に大きな数となる場合には、ある範囲の閾値を設定し、解析目的または精度を確保しながら効率的に変化数を算出することも可能である。例えば、閾値ΔZを定めて-ΔZから+ΔZまでをゼロゾーン(値0)として、この範囲を越えない限りは符号変化としては数えないこととする等も有効である。
こうした細かい変動は、ノイズの影響である可能性の他、同期して変化する成分波形の信号到達のわずかな時間差や、成分波形が生じる時間やその発生位置等の違いによる時間差等、によって生じている可能性もあり、このような状況を加味しながら解析を行うことも有用である。
【0127】
次にブロック分析部では、ブロック内の成分波形の時間的変化の大きさに基づくブロックの変動特性(以下、「ブロック変動特性」という。)、また、ブロックから出力される大きさ(以下、「ブロック出力特性」という。)及び、ブロックの出力特性に基づく成分波形集散度を分析することができる。
【0128】
ブロック変動特性は、ブロック分析部に設けられるブロック変動解析部(図示せず)により分析を行うことができる。ブロック変動特性は、ブロック内に含まれる複数の成分波形による挙動変化を平均的に表したものであり、ブロック変動度とブロック変動評価値で与えられる。ブロック変動度は、ブロック内に含まれる各成分波形による大きさや周波数特性等に基づく、ブロックの単位で表された特徴量の時間変化を示し、ブロックの平均挙動特性とブロック中の1成分波形当たりの特性変化量とに基づく値の時間変化率で表される。また、ブロック変動評価値に関しては、ある所定の時間範囲を定め、その範囲内に有するブロックの、ブロック変動度の平均値で表される。
【0129】
ブロックの平均挙動特性や、ブロックの1成分波形当たりの特性変化量は、ブロック内に含まれる成分波形の各特徴量に基づき求めることができ、また時間変化率については、ブロックの幅で除算して算出することができる。ここでブロック変動度及びブロック変動評価値の1つの実施例として筋電波形解析に用いた場合に基づいて以下に詳細に説明する。
【0130】
筋疲労が進行すると、筋繊維の信号伝達速度が低下することが知られている。伝達速度が遅くなれば、電位測定位置での電位変化の速度は緩やかになり、適合強度の値の時間変化も小さくなると推測できる。しかし、複数の運動単位活動が融合している状態では個々の運動単位での変化速度を適切に評価するのは難しいため、このような場合、ブロック変動特性を用いて評価することができる。具体的には、先ず、ブロック内での各運動単位活動電位波形の特徴量の変化量を運動単位活動電位波形の推定数、つまり、成分波形推定数で割ることで1成分波形当りの特性変化量を算出する。
【0131】
また、筋疲労が進行すると周波数成分の分布が低周波に寄ってくることも知られている。これは、ウェーブレット係数値による特性ではCoBまたはLCoBの値の低下、つまり周波数特徴量の低下となり、ブロック内に有する運動単位活動電位波形の周波数特徴量の低下となって現れる。そこで、ブロックの周波数特徴量として、各運動単位活動電位波形の周波数特徴量に適合度を乗算した値の平均値を算出する。
【0132】
それらの点から、運動単位活動がどのくらい活発か(低疲労で活動しているか)は、ベースとなる後者の値に個別の活動の強さとも言える前者の値とに基づき捉えることができると考えられるため、これをブロック当たりの運動単位活性係数とし、また、この運動単位活性係数は運動単位の疲労が進行すれば値が低下する性質を有するものとなるため、運動単位の活性度合い、すなわち、ブロックの変動特性を分析することができる。
【0133】
上述したように、筋が低疲労状態にある時の成分波形はより高い周波数までの周波数成分の増減といった周波数分布の時間変化を伴い、高疲労状態になると周波数成分分布の時間変化は低疲労時よりは低い周波数までの変化に留まる。また、筋が低疲労状態にある時には筋の信号伝達速度は高く、高疲労状態になると信号伝達速度が低下することが知られており、信号伝達速度が高ければ周波数特徴量の時間変化は大きく、速度が低ければ周波数特徴量の時間変化は緩やかになる。これを鑑み、1成分波形当りの特性変化量、つまり運動単位当りの変化量評価をsqrt((Σk(ddk
2))/CMU)として運動単位の活性度の評価に組み込む。ここで、CMUはブロックに含まれる運動単位活動電位波形の推定数であり、ddkは、ブロック内のMR*CoBサンプル(ssi、1<ssi<N、N+1;ブロック幅(BW))の値の間の特徴変化量dd0、dd1、・・・、ddN(ただし、ddi=ssi+1-ssi)である。ブロックの先頭と末尾ではブロック設定の閾値を横切るため、dd0>0、ddN<0となる。また、運動単位活性係数は分析対象とする運動単位活動電位波形に基づく評価を行うものであるため、使用する周波数特徴量については、そのブロックに含まれる運動単位活動電位波形が、分析対象の波形であること、すなわち成分波形らしさを示す適合度(MR)を周波数特徴量に乗算した値を使用している。
【0134】
従って、ブロック単位で規定される運動単位活性係数をMUActCoef(t)と表し、次式のように定義する。
(1)tはブロック外:MUActCoef(t)=0
(2)tはブロック内:MUActCoef(t)=(AVEC+sqrt((Σk(ddk
2))/CMU))*(CMU/BW)
同じブロックに属するtに関して、MUActCoef(t)は同じ値となる。
【0135】
なお、AVECはMR*CoB値の平均値であり、AVEC=sqrt(Σk(ssk)/BW=sqrt(Σk(ssk)/(N+1)(ただし、k=1,…,N)で与えられる。このように、ブロック変動特性は、ブロックの特徴量のベースとなる、ブロック内における各成分波形の特徴量平均値と、ブロックの特徴変化を示す、ブロックの1成分波形当たり特徴変化量に基づく値を、分析対象の時間幅とするブロック幅で除算することで表すことができ、ブロックの変動に関する情報を取得することができる。
【0136】
また、筋の疲労は一つの筋を構成する運動単位のすべてで均質に進行するわけではなく、低疲労の運動単位と疲労が進行した運動単位とが混じりあって活動していると考えるのが自然であろう。MUActCoefは一つのブロックに融合した運動単位活動群の平均的な活性度を評価するものであるため、疲労が進行していても低疲労のブロックと高疲労のブロックとが混在して出現する。
【0137】
それゆえ、筋疲労のような筋状態の評価は、特定の時区間内の複数のブロックの平均的な特性に基づいて行う必要がありブロック変動評価値を用いて評価することができる。具体的には、評価する所定の時間範囲を設定し、その範囲の各時刻のMUActCoef値の平均を求め、これを筋活性度評価値(Muscular Activity Evaluation Value;MActEval)と表して評価する。なお、所定の時間範囲については、解析目的や必要とする解析精度等に応じて適宜設定することが可能であり、例えば、後述の実施例では100~200[ms]を使用している。
【0138】
平均は通常の算術平均、二乗平均平方根(RMS)、調和平均等、様々な値があるが、解析目的や必要とする解析精度等に応じて、適宜選択し用いることができる。
【0139】
筋活性度評価値は、筋全体としての疲労(筋を構成する運動単位の平均的な疲労)が進行するほどに低い値となる性質を持つものである。運動単位活性係数で評価する場合には、直接的に運動単位の疲労を捉え把握することができ、所定の時間範囲を設定して筋活性度評価値で評価する場合は、例えば、筋疲労の進行で個々の運動単位活動がそれぞれの活動を補い合い、同期する必要性が増すとこの所定の時間範囲内のブロックが疎となる傾向が強まる等、の所定の時間内での総合的なブロックの挙動を捉えることが可能となる。
【0140】
筋活性度評価値の平均を求める際の考え方として、活動中の状態だけに基づいて評価する(主として運動単位の状態を基準とした評価とする)という考え方から所定の時間範囲内でブロックの部分だけを対象にする方法と、ブロックが疎になる傾向が強まるという性質まで含めて評価する(運動単位活動の同期必要性まで含めた総合的な筋活動評価とする)ブロック外の部分も値0として含める方法とが可能である。総合的な筋活動の評価を行う目的や状況に応じて、前記方法のいずれかまたは両方を組み合わせて使用することができる。また、ブロックが疎になる傾向を直接的に扱うために、所定の時間範囲の時間幅に対して、ブロックとなっている時間が占める割合を評価して組み合わせることも可能である。なお、後述の筋活動余力係数と組み合わせて筋活動の評価を行おうとしている場合には、ブロックが疎になる傾向は筋活動余力係数での評価要素としての性質の方が強いため、筋活性度評価値ではブロックの部分だけでの平均とする方が良いだろう。
【0141】
また、筋電波形信号の強度(振幅)は、筋活動強度(発揮力)との関係性が強いことは良く知られている。しかしながら振幅は、発揮力が大きくなったときに増大するだけでなく、同じ発揮力でも筋疲労が進行すると増大することも良く知られており、変化する発揮力の算定が難しい運動中に筋電波形の信号強度で筋疲労を評価することが難しくなる場合も想定される。
【0142】
その場合、運動単位活性係数(MUActCoef)は、ほぼ周波数特性のみに依存して算出されるものであるため、発揮力で変動する信号強度の影響を受けにくく、負荷(発揮力)が大きく変動する間でも、この運動単位活性係数に基づく筋活性度評価値(MActEval)の安定性は良好である。
【0143】
また、所定の時間範囲を短時間単位でも評価が行えるため、運動中の頻繁な筋活動変化への追随性も高い。こうした機能特性を持つことによりが、各種分析をリアルタイムで行うことが可能となる。
【0144】
[変形例]
ブロック変動評価値に関する変形例を以下に示す。ブロック変動評価値は、特定の時区間内のブロックの平均的な特性を示すものであるが、その特性における変化を詳細に分析しようとする場合、その変化が微細である等で、グラフ等で捉えることが難しい状況も想定される。そこで、ブロック変動評価値をvalueとして、次式により、べき乗計算を行いブロック変動評価値を数値的に拡大して、その特性をより明確に分析することができる。
Z=Y*((X(scale*(base-value))) ― zbase)
Z;ブロック変動評価値のスケール拡大値
Y;Zの値域の範囲をコントロールする実数
zbase;Zの値域の0位置をコントロールする実数
X;拡大量の基準となる正の実数
scale;拡大倍率の程度を定める非0の実数
base;拡大の中心(倍率1.0)となる値
なお、特段の事情がなければXの値は2、Yの値は1、zbaseは0として、(base-value)の値に応じたZの値の増加の程度(傾きの増加量)をscaleで設定することが好ましい。
それ以外の値については、(base-value)が取る値の範囲でのZの値の増加の程度(傾きの増加量)を考慮してXとscaleを、更にXとscaleを定めたときのZの値域を考慮してYとzbaseを設定することができる。
ここで、上記の拡大値に対する1つの実施例として筋電波形解析に用いた場合について以下に詳細に説明する。
【0145】
前述した筋活性度評価値は疲労による変化を十分敏感に捉えてはいるが、その変化は視覚的には明確に現れない場合があることに加え、投入される運動単位の疲労程度のバラつきなどによる値の変動も存在するため、そのままの値をグラフにするなどして提示しても変化による差異を視認することが難しい場合も想定される。
【0146】
筋活性度評価値は、「どの程度活発か」と「どの程度疲労しているか」の二つを同時に示していると考えることができるため、そこで本発明では疲労による変化を拡大(相対的に活性度による変化は縮小)するように、数式を定義して変換した値を用いて分析を行うことができる。また、この変換した値は、疲労による変化が表されることから筋電波形分析においては、筋疲労度変化係数(FatigueCoef)と定義して筋の疲労の解析を行うことができる。この係数は、筋疲労度が低い(活発度が高い)状態では0に近い値と取り、疲労度が高くなるほど大きな値を示す。
【0147】
差異を強調する変換式は様々考えることができ、特に制約も無く、値の本質も変わるわけではないため、この変換式によって本発明の範疇を外れることはないが、一例として、本発明では、筋活性度評価値valueに対して2(scale*(base-value))を用いることができる。筋活性度評価値は疲労が進行すると値が小さくなる傾向があるが、その変化による差の数値は小さいものであるため、差異を分かりやすくするためにべき乗で値を拡大するのがscaleの役目である。baseとscaleは、筋活性度評価値に基づく分析を行う際の感度調整としての機能を持つ。すなわち、baseは分析上の通常時(例えば疲労なしや低活動)に属する活動か否かの境界となる値であり、scaleはどの程度の小さな量の境界越えまで細かく検出したいかによって定める値である。
【0148】
baseを大きく筋疲労した時のおよその筋活性度評価値の値としてscaleを負の値とした場合には、低疲労な筋活動であるほど高い値を示すものとなるし、逆にbaseを低疲労時のおよその筋活性度評価値の値としてscaleを正の値とした場合には、筋疲労が進行するほど高い値を示すものとなる。例えば、筋疲労がどの程度進行しているかを分析する場合には、baseを電極装着時(低疲労時)に計測した値に基づいて設定(scaleは正値)し行うことができる。
【0149】
なお、後述する実施例では、低疲労時の筋活性度評価値の値(低疲労と判断する範囲での低めの値を推奨)を基準値referenceとし、base=reference/2、scale=1/(0.2*reference)としたものを筋疲労度評価係数(FatigueCoef)として行った。
【0150】
[ブロック出力特性]
次に、ブロック出力特性について説明する。ブロック出力特性は、ブロック分析部に設けられるブロック出力特性解析部(図示せず)により分析を行うことができる。ブロック出力特性は、ブロック内に含まれる成分波形の挙動により生じる大きさを示すものでありブロック出力度とブロック出力推定値で表され、これら成分波形の特徴量に基づき算出される。先ずは、ブロック出力度について説明する。
【0151】
ブロック内には複数の成分波形が重畳して存在しており、このような状況下においては、ブロックから出力される大きさは各成分波形の振幅に基づき評価されることが一般的である。この場合、各成分波形はその信号伝達速度が変化しながら変動していることが考えられるため、その信号伝達速度の変化による成分波形の周波数変化も考慮すると、大きさに対する分析精度の更なる向上が期待できる。つまり、ブロック内の各成分波形の周波数変化に伴う振幅の変動も考慮すると、より的確にブロックにより生じる大きさの評価が可能となる。
【0152】
そこで、ブロックにより生じる大きさの変化度をブロック出力度とし、ブロック内の各成分波形の大きさ(強度)を示す信号強度特徴量(WSP)とブロックの周波数特性を加味した時間的な変動特性を示す、ブロック変動度に基づき表すことができる。具体的には、ブロック内の各成分波形の、適合度(MR)を乗算した信号強度特徴(WSP)にブロック変動度を乗算し設定される。ここで、ブロック出力度の実施例として、筋電波形解析で用いる場合により詳細に説明する。
【0153】
筋電波形に基づき筋により生じ得る力(以下、「発揮力」という。)を推定する場合、所定の時間範囲の計測値の整流平均(ARV)や二乗平均平方根(RMS)を使うことが一般的である。一方、振幅の増減は筋疲労の影響も受けるため、筋疲労の進行を考慮する必要があり、振幅への依存性が低く、高い感度で筋疲労の程度を捉えることができる手法が求められる。
【0154】
そこでブロック出力特性解析部では、筋の発揮力推定のパラメータとして、ブロック変動度である運動単位活性係数と、ブロック内の成分波形である各運動単位活動波形の信号強度特徴量(WSP)とに基づく値(以下、「発揮力評価係数」という。)により、筋の発揮力の推定、評価を行うことができる。筋の発揮力は、成分波形である運動単位活動波形の存在時刻とされた瞬間のみに現れるわけではなく、運動単位活動波形が存在する時区間の間に、筋収縮の過程に依存した発揮力の変化が存在する。この変化の瞬間値を、各時刻のWSPに適合度(MR)を乗算することで近似する。具体的に、発揮力評価係数(PowerCoef(t))は、時刻tの運動単位活性係数、信号強度特徴量、適合度をそれぞれMUActCoef(t)、WSP(t)、MR(t)として、次式で表すことができる。
PowerCoef(t)=MUActCoef(t)*sqrt(MR(t)*WSP(t))
なお、前述のように、1つのブロック内のtにおいてMUActCoef(t)は同一値である。
【0155】
運動単位活性係数は筋活動の根幹を成す運動単位の観点で活動効率を捉えるものであり、信号強度特徴量(WSP)は同時に活動している運動単位(筋繊維)活動の大きさやパワー(強度)であり、発揮力評価係数は、これらの掛け合わせ、すなわち、運動単位基準での「効率×量」であり、発揮力評価の基準となる特徴量として表すことができる。
【0156】
次にブロック出力推定値について説明する。ブロック出力推定値は、ある所定の時間範囲を定め、その範囲内に有するブロックの、ブロック出力度の平均値で表される。ブロック出力推定値は、ブロック出力度がある時刻に生じるブロックの大きさ(強度)の瞬間値であることから、このブロック出力度の平均値に基づき、この所定の時間範囲における、ブロック内の成分波形の挙動による、ブロックの出力の推移等を、把握することができるものである。なお、平均値としては、通常の算術平均の他、移動平均や、二乗平均平方根(RMS)、調和平均等、様々な値があるが、解析目的や必要とする解析精度等に応じて、適宜選択し用いることができる。
【0157】
ここでブロック出力推定値について、筋電波形分析に用いた場合を例として以下に説明する。上述した発揮力評価係数は運動単位活動が励起された時の瞬間値であるが、一方、例えば負荷を連続的に加えて筋の活動状況を捉える等、ある連続した時間の中で励起された筋の活動の時間持続性の評価を必要する場合は、ブロック出力推定値を用いて分析することができる。この場合、ブロック出力推定値を、上述した発揮力評価係数に対応して発揮力推定値(Power Estimation value;PowerEst)と称して表すことができ、発揮力推定値は例えば、発揮力評価係数の移動平均によって算出することができ、発揮力評価係数は運動中の頻繁な筋活動変化への追随性も高いため、リアルタイムにかつ的確に筋の発揮力の評価が可能となる。発揮力推定値と負荷の大きさとの関係は、筋疲労状態の変化がある場合でも、筋電波形信号の振幅のRMSなどよりも遥かに良好な直線性を有する。
【0158】
一例として、速筋と遅筋による筋活動における発揮力の推定に関しては、速筋と遅筋のそれぞれで発揮力評価係数の移動平均を求めてそれぞれの発揮力を推定し、それらを合計することで筋全体の発揮力を推定することができる。
【0159】
次に成分波形集散度について説明する。解析対象信号中の成分波形の動向としては、複数の成分波形の出現がどの程度分散しているか、あるいは同期しているかといった成分波形の散らばりの程度(以下、「成分波形集散度」という)も重要な特性であると考える。
例えば、複数の成分波形の出現が近い時刻に集中していれば、それらは1つのブロックを形成し、ブロック間には空白期が出現する可能性が高い。さらに同期性が強くなるなど、ブロックを形成する成分波形の出現の仕方の画一化が進行すれば、集中によって単位時間あたりのブロック数は減少し、各ブロックの形状的な特徴のばらつきも小さなものになると予想される。
【0160】
上述したブロック出力特性では、ブロックから出力される大きさの時間変化率や、ある所定時間範囲における各ブロックの大きさの平均的な値を推定、把握できるが。一方、これらブロックが時間的な変化に伴って、上記のような、ブロック同士の融合等によるその大きさや数の変化または、所定時間範囲内における各ブロックの時間間隔等のブロックの出現分布等の変化をとらえることができれば非常に有用である。
【0161】
そこで、ブロック分析部に集散度分析部(図示せず)を設け、成分波形集散度を定量的に算出し、このような状況下において、所定時間におけるブロックによる出力の変動状況及び推移を把握し信号波形中の成分波形の動向を評価することができる。
【0162】
[成分波形集散度]
成分波形集散度(Clumping Coefficient;ClumpingCoef)に関して、上述のようにブロックの大きさ等の変動が進行すると、所定時間における以下のssrr値の変化が予測される。ただし、所定時間内のブロック(またはブロック間空白期間)の数をN、各ブロック(またはブロック間空白期間)での対象値をx1,…,xNとするとき、ave=(Σi(xi))/Nであり、ssrr=Σi(((xi-ave)/ave)2)である。
(1)ブロックの時間幅およびブロック間空白期間の時間幅に関して、所定時間内の全ブロックの時間幅の平均(ave1)に対する残差を平均値との比率とした値の二乗和(ssrr1)、および所定時間内の全ブロック間空白期間の時間幅の平均(ave2)に対する残差を平均値との比率とした値の二乗和(ssrr2)
(2)ブロックの平均値に関して、所定時間内の全ブロックでの平均(ave3)に対する残差を平均値との比率とした値の二乗和(ssrr3)
(3)ブロック内成分波形推定数に関して、所定時間内の全ブロックでの平均(ave4)に対する残差を平均値との比率とした値の二乗和(ssrr4)
(4)ブロック内のサンプル間変化量の二乗和に関して、所定時間内の全ブロックでの平均(ave5)に対する残差を平均値との比率とした値の二乗和(ssrr5)
【0163】
若しくは、成分波形集散度に関して以下のような値の変化が予測される。これらは、解析時に必要に応じて適宜選択可能である。
(1)ブロックの数(p1)、およびブロック間空白期間の数(p2)
(2)ブロックの平均値に関して、所定時間内の全ブロックでの平均に対する残差の二乗和(p3)
(3)ブロック内成分波形推定数に関して、所定時間内の全ブロックでの平均に対する残差の二乗和(p4)
(4)ブロック内のサンプル間変化量の二乗和に関して、所定時間内の全ブロックでの平均に対する残差の二乗和(p5)
【0164】
ブロックとブロック間空白期間との数や時間幅は独立した関係とは言い難いため、ssrr1とssrr2、p1とp2は、それぞれの組の両方を評価に含めることを推奨する。なお、評価に一方だけを含めるようにすることを禁ずるわけではない。p1とp2は実質的に同じことを指しているが、所定時間の幅の解析区間の境界での数え上げ状況が悪影響を与える可能性を低減するために、敢えて両方を要素として含めている。
独立性の弱い二つの特徴量の両方を独立に含めてしまうとその特性の影響力が不平等に大きくなってしまうため、ssrr1とssrr2とではsqrt(ssrr1*ssrr2)を、p1とp2とでは成分波形集散度の評価を行える状況では値が1以上であることを鑑みてsqrt((p1-1)*(p2-1))を、評価式に組み入れる上での一つの特徴量として扱うこととする。
ssrr1とssrr2の組は成分波形の集散度の違いがブロック形状の違いに現れる傾向が強い場合に向いており、p1とp2の組は集散度の違いが所定時間の幅の中のブロックの数の増減に現れる傾向が強い場合に向いている。
【0165】
ssrr3とp3、ssrr4とp4、ssrr5とp5のそれぞれの間には強い関連性があり、前者(ssrr)は変動の程度を相対量で扱うのに対し、後者(p)は変動の程度を絶対量で扱う。成分波形の集散度の違いが、平均値の違いにはあまり関係なく変動の過多に現れる傾向が強いなら、絶対量を選択する方が向いている。そうではなく、平均値が小さくなればそれに引きずられて変動の量も小さくなる傾向が強いなら、相対量を選択する方が向いている。
よって、対象とする成分波形の特性に応じて、以降においてssrr3~5で記述されているものをp3~p5に置き換えても構わない。もちろん、両者の違いに基づき別種の特徴量として扱い、両方ともに計算式の中に含めてもよい。
【0166】
所定時間の幅の解析区間を移動させつつ解析対象信号の解析を行う場合、ブロックまたはブロック間空白期間が解析区間の境界を跨ぐのが通常である。この場合に、境界に位置するブロックやブロック間空白期間から解析区間内の部分だけを切り出して計算に含めるのは不適切である。ブロックやブロック間空白期間が境界を跨ぐ場合には、その全体を含めるか除外するかを選択するようにする。この境界処理の結果として、各所定時間での計算に関与するサンプルの総数(winN)は変動する。
【0167】
これらの数値を用いて、成分波形集散度(Clumping Coef)を次の式で与える。
( (sqrt(ssrr1 * ssrr2) * ssrr3 * ssrr4 * ssrr5)(1/4) ) / winN
この式は前出の(1)~(4)のssrr(ただし、(1)については二つのssrrの相乗平均)の相乗平均を所定時間の時間幅で割ったものとなっている。なお、ssrr1~ssrr5のすべての利用が必須というわけではなく、一部だけを選択して、それらの相乗平均をwinNで割ったものを使用しても良い。ssrr1とssrr2の一方だけを使用する場合には、それらの相乗平均の代わりに選択したssrrをそのまま用いる。また,前出の(1)~(4)の一部だけを選択した場合、あるいはssrrとpとの両方を導入した結果として選択数が増えた場合は、選択したものの相乗平均を成分波形集散度の式の分子として用いることとする。
【0168】
解析対象信号中にk種の成分波形が存在し、すべての種類の成分波形の総合的な集散度を得たい場合には、以下のように扱う。
(1)すべての成分波形のブロックを重ね合わせた上で、ssrr1とssrr2、あるいはp1とp2を求める。
(2)ssrr3からssrr5のそれぞれについて、k種の成分波形のssrrをssrr1,...,ssrrkとするとき、総合的な集散度を求める際のssrrは(Πi(ssrri))(1/k)で与える。
【0169】
解析対象信号中の成分波形が離散的に出現する場合,成分波形の形状が安定し,出現が周期的であるほど、ClumpingCoefの値は小さくなる。複数の成分波形が時間的に重なりを持って出現する場合は、集中して出現するほどClumpingCoefの値は小さくなる。これにより、解析対象信号中の成分波形の出現動向の一側面を数値化して評価することができる。
成分波形集散度(ClumpingCoef)はそのまま用いても良いが、対象とする成分波形に特有の特性が存在する場合には、それと組み合わせて利用することもできる。この成分波形集散度の1つの実施例として、運動単位活動に特有の特性を加えつつ筋電波形分析に適用したものを以下に説明する。
【0170】
筋活動において負荷が大きくなると、発揮力を増大させるために多くの運動単位が同時投入される。筋を構成する運動単位の数は有限であるため、同時投入量が増えると、筋の発揮力を安定させるために、運動単位活動を時間的に広く分散させるための運動単位が不足して運動単位活動の集中が生じる。
【0171】
こうした状況は筋疲労が進行した際に多く観測されるが、それは筋疲労の直接の結果ではなく、個々の運動単位の発揮力が疲労によって低下したことで相対的に負荷が増大した結果であると言える。
【0172】
運動単位活動の集中は、適合度の推移で見ると、所定時間内に発生するブロック数の減少と、ブロックの形状特徴の画一化となって現れる。つまり、ブロック数が減少し、かつ、ブロックの各種特徴量のバラつきが失われた状況は、それらに変化を及ぼすような運動単位活動の違いを生じさせる余裕がないというような、活動限界の状態であると考えることができる。こうしたバラつき具合の大小に基づき、運動単位活動の限界に対してどの程度の余裕を筋が有しているかを評価するための特徴量を活動余力係数(Tolerance Coefficient of Muscular Activity;ToleranceCoef)と称する。
運動単位活動においては、筋活動が少ない場合にも、運動単位活動の集中と似たようなブロックの離散化が生じる。活動限界を考える上では、これらを同一視しないようにすることが必要である。そこで、低活動による離散化と比較して、高負荷や高疲労による集中での離散化の場合はブロックの幅が増大する傾向があることを鑑み、活動余力係数(ToleranceCoef)は次の式で与える。
ToleranceCoef=(ssrr2/ssrr1)*ClumpingCoef
なお、この式において筋疲労と低活動との差異に関する調整項としている(ssrr2/ssrr1)の逆数(ssrr1/ssrr2)は、筋疲労の進行による集中で、値が低下する傾向がある。そのため、成分波形集散度の評価のことだけを考えるなら、(ssrr1/ssrr2)を評価のパラメータとして導入することは可能である。
【0173】
活動余力係数の値は、ブロックのバラつきが大きい時、例えば運動単位の部分集合だけで必要な発揮力を得ることができており、励起タイミングを分散させて発揮力を安定できている時等、には大きな値となり、ブロックのバラつきが小さい時、例えば運動単位の活動を集中させなければ必要な発揮力を得られない時等、には小さな値となるであろう。
また、低疲労時のように個々の運動単位の活性度が十分に大きい状況や低活動時には、運動単位活動の集中のさせ方にも余裕があるために活動余力係数の値の変動は大きくなりやすいが、疲労が進行して運動単位の活性度が低下した時や高活動時には、集中のさせ方の余裕も低下し、活動余力係数の値の変動も小さくなりやすいという傾向もある。よって、活動余力係数に基づいて評価を行おうとする場合には、短時間での値だけではなく、もう少し長い時間幅での値の振れ幅といった動向にも目を向けることが望ましい。
【0174】
筋活動がない状態では評価できないし、筋活動がわずかである場合は評価結果が不適切となりやすいため、活動余力係数が利用できるのはある程度以上の筋活動があるときが好適である。
【0175】
成分波形集散度や、それを運動単位活動に拡張して適用した活動余力係数は、解析対象信号中での成分波形の分布を扱うものであるが、分布の違いによって生じる成分波形の活動加算の効果については含んでいない。加算効果の評価を望む場合は、成分波形集散度等に、信号強度等、評価したい加算影響を扱うための特徴量を組み合わせればよい。運動単位活動における実施例を以下に説明する。
【0176】
活動余力係数は、所定時間において現状から更なる発揮力を得るために運動単位活動を集中させる余裕がどの程度あるかを数値化したものと言うことができるが、力としての余裕を直接的に評価するものではない。力の大きさを扱うには、同所定時間における運動単位の活性度も考慮する必要がある。力としての余力の大きさを推定した値を活動余力推定値(Tolerance Estimation Value;ToleranceEst)とするとき、ToleranceEst=(ToleranceCoefR)*ave(MActCoef)で与えることができる。なお、Rは、ToleranceCoefの変化に対する重み付けであり、特にR=1/2が好適である。また、ave(MActCoef)は、ToleranceCoefを算出する所定時間でのMActCoefの平均値である。平均値としては、対象の所定時間での単純な算術平均ではなく二乗平均平方根を使っても良い。
【0177】
活動余力推定値(ToleranceEst)は、値の大小によって力としての余力の大小を示すものであるが、発揮力推定値(PowerEst)と同様に力の大きさを直接に表すものではない。余力の相対的な変化を追うだけであればそのままで問題はないが、筋活動を評価する上で関連性が強い発揮力推定値との間での比較や加算を行うには値のスケールの違いを調整する必要がある。そのため、調整項としての定数であるスケールファクター(sf)を導入した推定値を活動余力調整推定値SToleranceEst=sf*ToleranceEstとし、PowerEstとの対比の際にはこれを用いることとする。
sfとして与える値によって、発揮力と余力との対比結果の見かけは極めて大きく変化する。仮に筋発揮力の真の限界値というようなものが存在し、かつ正確に計測できるのであれば、それをsfの最適解を得るための基準とすることもできるだろう。しかし、そのような理想的な基準値が得られることはほぼ期待できず、「正解」が不明であるために、適切なsfの値を定めることは簡単ではない。sfを与えるための方法の一つとして、低疲労の状態で発揮力を弱い力から十分に大きな力まで急激にならない程度に数回往復させた際のデータに基づき、ToleranceEstにおいて外れ値を除いたおよそのピーク値peakTと、PowerEstにおいて外れ値を除いたおよそのピーク値peakPとを求め、sf=peakP/peakTを設定値とする方法を用いてもよい。得られた設定値を必要に応じて微調整してもよい。また、100%MVCでも発揮力の真の限界値とまでは言えないが、100%MVCを計測した上で、低疲労時の活動におけるPowerEst+SToleranceEstの振れ幅上限あるいは平均的な値が比較的安定した値となるように調整しても良いだろう。
【0178】
活動余力係数に基づく評価は、筋疲労との関係は強いものの、別視点での評価である。sfの設定の難しさに加え、値の時間変動(振れ幅)も大きいため、活動余力推定値においても余力の具体的な値を厳密に与えるようなものではないが、発揮力が大きいときは余力が減少し、発揮力が低下すると余力が回復してくるというような傾向や、筋疲労が進行するほど推定される余力は低下し、値の振れ幅も小さくなるというような傾向を有している。また、継続している運動の負荷の程度が被験者にとって大きいほど、活動余力推定値の観点での余力の悪化の速度は高くなる傾向があり、悪化を始めてから運動の限界に近づく過程では、値や振れ幅といった特性の時間変化が直線的であることが多い。こうした性質は、運動時の負荷量の管理や限界点の見極めなどに極めて有用であると考える。
【0179】
本発明では、成分波形の正常範囲での変化の方向性という概念の理解しやすさを鑑み、1つの実施例として非侵襲で計測できる筋電波形信号の成分波形である運動単位活動電位波形を題材として説明するが、本発明の手法の適用範囲は筋電波形信号に限られるものではない。
【0180】
加えて、波形に限らず、正常/異常の判別が必要な問題領域におけるデータ集合において、すべてのベクトル成分が0以上という条件でデータを特徴ベクトル化することが可能で、正常とされる特徴ベクトルの変動範囲に方向性が存在するような場合にも、本発明の手法を適用することができる。
【0181】
本発明での成分波形検出・抽出手法、及びそれらに基づく、成分波形の挙動分析をリアルタイムに行う手法を提供する。
解析対象とする信号波形の挙動を捉える上で、その中に含まれる成分波形の変化量や大きさ等の特性を把握し評価することは重要であり、従来技術では特定条件の付与や複雑なモデル化を構築した上で行われるものがほとんどであり、事前予測できず短時間で、かつリアルタイムに評価することはほぼ不可能であった。
【0182】
特開2016-063995号公報に記載の技術はそれを可能とすることを目指した数少ない技術の一つであるが、主に観測された波形に現れる現象を可視化することに着目し評価することを目的したものである。本発明は更に発展させ、その現象について、より定量的な評価を行うことができるように、その現象を生じさせる成分波形の各種特徴量を算出し、その特徴量に基づき成分波形の挙動を評価可能にしたものである
【0183】
本発明によれば、次のような効果を得ることができる。
・従来技術ではほぼ不可能であったような、解析対象とする信号波形中のターゲットとする成分波形の抽出、及びその成分波形の定量的な挙動分析を1つのチャンネルで、かつリアルタイムでの高精度、高感度な評価が可能になる。
・本発明を、特に筋電波形分析に用いた場合、高精度な筋疲労評価の実現により、運動中の筋疲労や発揮力のリアルタイム評価の精度を向上することができる。また、従来技術にはない高度な筋活動状態情報を提供することが可能になることで、医学、スポーツ科学の分野の、技術発展に寄与する。
・もちろん、本発明は不規則信号波形からその中に含まれ、解析のターゲットとする成分波形を抽出、及び分析することができることから、様々な解析対象波形および成分波形についても適用できる。例えば、黙声認識(いわゆる、口パク)や、心音、心電図による電気的波形、筋肉の活動等、を捉える筋電波形、または工学的には、機器や各種環境から生じる騒音、振動波形、または機器等の故障検知を目的とした波形を対象とすることができる。
【0184】
非侵襲で計測できる表面筋電から筋活動の概略を捉えるのではなく、同じ表面筋電を使いつつも筋活動を構成する運動単位活動の観点での分析を可能にする。これにより、速筋、遅筋の寄与バランスや筋疲労の進行を運動中に直接的に調査することができ、基礎医学やスポーツ医学等の発展に寄与できることに加え、筋電義手やパワーアシストなどの筋電活用機器の精緻化に役立つ。ただし、本技術の適用範囲は表面筋電に留まらない。
【実施例0185】
[実施例1]
以下、実施例により本発明を更に詳細に説明するが、本発明は、その要旨を変更しない限り以下の実施例に限定されるものではない。
まず、実施例1として、本発明の作用効果を確認するために後述するような実験を行った。
【0186】
[計測方法]
本実験における計測方法を、以下に示す。
(1)直立姿勢を取り、肘を90度に曲げて前腕を水平に保持する。
(2)手首位置にベルトをかけて加重(前腕は極力水平を保持する。加重増減時間は、それぞれおよそ3秒)した際の、上腕二頭筋の表面筋電を計測する。
【0187】
本実験における実験方法を、以下に示す。
[実験方法]
(1)上記計測方法での計測を、約20秒の脱力休憩時間を挟んで4回実施する。
(2)約1分間のダンベル昇降運動を行い、計測対象の上腕二頭筋を疲労させる。
(3)上記計測方法での計測を、約10秒の脱力休憩時間を挟んで5回実施する。
ここで、上記(1)におけるフェーズを「低疲労時」、上記(2)におけるフェーズを「筋疲労運動」、上記(3)におけるフェーズを「疲労後」と言う。つまり、低疲労時に4回、疲労後に5回、それぞれ上記計測方法により筋電波形の計測を行った。
【0188】
図6~15は、本発明を用いて行った実験を説明するための図である。
図6は、低疲労時4回、疲労後5回の各計測時の加重値の推移を示すグラフである。
図6に示すように、およそ3000ms~4000msまで加重が増加しており、以後は加重が減少している。また、各回の加重値の推移の違いは十分に小さいということができる。
【0189】
また、
図7は、各計測(低疲労時4回、疲労後5回の計測)における筋電波形の変化を示すグラフである。つまり、上記実験方法の(1)~(3)のそれぞれのタイミングで計測された筋電波形のグラフが示されている。この図から、低疲労時の4回の計測では、回ごとの筋電波形に大きな変化は確認されない。加重の増減に呼応して、振幅の増減が生じている。一方、疲労後の直後(
図7の右下)には筋電波形は大きく変化し、その後、10秒ずつの休憩を挟むに従い、次第に筋電波形の変化が収まっていることが分かる。つまり、わずか10秒ずつの休憩を挟んでも、被計測者の腕(上腕二頭筋)が次第に回復している様子が見て取れる。
【0190】
また、
図8は、各計測における表面筋電のRMS値(窓幅200ms)の変化を示すグラフである。また、
図9は、横軸に加重値、縦軸にRMS値としてプロットした結果の変化を示すグラフである。RMS値は表面筋電の振幅に依存した値で、発揮力評価のために従来技術で用いることが多い推定方法であるが、
図8,9からも分かるように、疲労前後での値の特性の大幅な変化が見られる。つまり、筋疲労の影響を無視できる範囲では安定していると言えるが、筋疲労の影響の程度を評価できない場合には、これら全ての状態の可能性を否定できないことから,全ての状態のサンプルを混ぜ合わせた分布の上で相関を考える必要があるため、RMS値から発揮力を推定する精度は極めて低くなる。
【0191】
また、
図10は、本発明を用いた解析の実施例を示す図である。なお、解析における条件を以下に示す。
・冗長ウェーブレット解析の最大深さ:レベル-11
・ウェーブレット係数集合設定:T+7/16(T+14/32)
・速筋用と遅筋用の基本ベクトルと、適合フィルタおよび分割フィルタを設定(
図15上部参照)
上記条件で、低疲労時1回目と疲労後1回目においてそれぞれ、ウェーブレット係数集合に基づいて係数集合の係数群を得て、速筋および遅筋の適合度を算出し、分配型調整フィルタ(速筋・遅筋)を適用して、成分波形(のウエーブレット係数群)を検出している。
【0192】
また、
図11は、同様に低疲労時1回目と疲労後1回目においてそれぞれ、速筋および遅筋の適合度を算出し、分配型調整フィルタ(速筋・遅筋)を適用して、速筋および遅筋のMR×LCoBをグラフ化したものである。これにより、筋疲労が起こると、高周波成分が減少する、つまり成分波形の鋭い電位変化の喪失が起こり、大きな値の頻度低下が観測されることが分かる。また、筋疲労が起こると、成分波形の電位変化速度が減少する、つまり成分波形形状の鈍化が起こり、時間変動の沈静化が観測されることが分かる。なお、筋疲労が起こると、双極誘導電極間電位の変化速度低下、つまり伝達速度の低下が起こり、個々の山の幅の増大が観測されることが分かる。
【0193】
また、
図12は、ウェーブレット係数集合の適合調整の例を示す図である。解析における条件は、上述したように以下に示すものである。
・冗長ウェーブレット解析の最大深さ:レベル-11
・ウェーブレット係数集合設定:T+7/16(T+14/32)
つまり、ベースとしたウェーブレット係数集合設定は「T+7/16(T+14/32)」であり、速筋用と遅筋用の基本ベクトルと、適合フィルタおよび分割フィルタは、ベースとした上記ウェーブレット係数集合で設定したものを使用した。
【0194】
この結果、速筋と遅筋の増減RMS比の変化を見ると(
図12の下部参照)、低疲労時はT+7/32付近で、疲労時はT+12/32付近で交差している。そのため、T+13/32を最適なウェーブレット係数集合として採用した。
【0195】
図13は、この最適なウェーブレット係数集合を用いて、
図11に示したものと同様に、低疲労時1回目と疲労後1回目においてそれぞれ、速筋および遅筋の適合度を算出し、分配型調整フィルタ(速筋・遅筋)を適用して、速筋および遅筋のMR×LCoBをグラフ化したものである。そして、さらに運動単位活性係数を算出したものである。
【0196】
また、
図14は、本実験の各計測における発揮力評価係数の変化を示したものである。各計測において、上図が速筋で、下図が遅筋の発揮力評価係数を示すグラフである。これにより、筋疲労の影響を組み込んだ評価により、表面筋電の振幅やRMS値、WRMS値(ウェーブレット係数集合の係数群の係数値の二乗平均平方根の値であり、WSPとは平均を取るための定数項で割っているかの違いだけで、値の特性としては同じである)に見られたような筋疲労前後での値の特性の大幅な変化は見られず、筋疲労による悪影響を低減できていることが分かる。
【0197】
また、
図15は、本実験の各計測における横軸を加重値、縦軸を発揮力推定値(窓幅200msの移動平均)とした結果の変化を示したものである。これにより、従来技術のRMS値での評価に比べると、加重(負荷)との相関が遥かに良好であると言える。なお、加重増大時に高め、加重減少時に低めの値を示しているのは、肘角度維持のための拮抗筋との協調による影響であると推測される。
【0198】
[実施例2-1]
次に、実施例2-1として、電気刺激時の筋音(Mechanomyogram:MMG)波形を題材としたウェーブレット係数集合設定の実施例を示す。
【0199】
まず、下記引用文献を参照して、
図16に示す仮想MMGデータを作成した。当該データは、電気刺激によって生じた筋音波形を人為的に再現したものである。
引用文献:Evaluation method for muscles, measuring mechanomyogram induced by electrical muscle Stimulation using lead zirconate titanate-based acoustic sensor”, Japanese Journal of Applied Physics 58, SLLD11 (2019)
【0200】
そして、この仮想MMGデータを「Daubechies N=2」を使って冗長ウェーブレット解析を行い、その結果から正値係数のみを抽出した。
図17が、抽出したものを3D表示させたグラフ、
図18が、抽出したものを等高線表示させたグラフである。ただし、
図18の曲線(
図18の矢印参照)は、副波形成分としてのまとまりを示すために追記したものである。
【0201】
ここで、
図17,18に示されたMMG波形が解析対象信号の成分波形であると仮定する。仮想MMG波形データを解析した結果から判断すると、MMG波形は2種類の副成分波形の複合によって形成されていると捉えることができる(
図18の2本の曲線で示されたまとまり部分)。
【0202】
冗長ウェーブレット解析において複数種類のウェーブレット係数集合を同時に出力するようにしても、計算コストはほとんど変わらない。そのため、このような成分波形抽出のためのウェーブレット係数集合設定と分析手順は、大きく分けると
(1)副成分波形のそれぞれについて独立に係数集合を設定して個別に解析を行った上で、対応する時刻(例えば
図18の点線位置)のそれぞれの適合度等で総合評価する。
(2)副成分波形ごとの係数集合を連結して1つのウェーブレット係数集合として用いることで、シンプルな性質の成分波形の場合と同じ手順で分析する。
というものになり得る。
【0203】
上記(1)は、処理の手間は増えるが、副成分波形の特徴変化があり得る場合の分析は行いやすい。それに対し、上記(2)は上記(1)と逆の性質を持ち、成分波形の変形を気にする必要がない場合には簡便な方法である。
【0204】
[実施例2-2]
続いて、MMG波形を重畳(合成)させた波形を作成した。電気刺激が与えられる時間間隔をランダムに変動させ、更に、足し合わせるMMG波形形の強度もランダムに変動させた。
【0205】
図19に示すグラフが、MMG波形を重畳させて作成した波形である。横軸の数値は、電気刺激が与えられた時刻、すなわち各MMG波形の時刻0に相当する。
【0206】
また、
図20は、
図19に示す合成波形を「Daubechies N=2」を使って冗長ウェーブレット解析を行い、その結果を等高線表示させたグラフである。ウェーブレット係数集合としてはH+0に相当する。なお、縦軸のレベル(周波数帯域)の数値の前に「Hi」と「Lo」の文字が付与されているが、これは
図18の2本の曲線の高周波数側(Hi)と低周波数側(Lo)の信号成分が含まれている周波数帯域との対応を示す印である。例えば、高周波数側(Hi)と低周波数側(Lo)のどちらにもレベル:-6が存在するが(
図18参照)、「HiL-6」は(
図20参照)、高周波数側(Hi)のレベル:-6を示している。
【0207】
更に、
図21は、
図20からウェーブレット係数が正の部分(正値係数)のみを抽出して等高線表示させたグラフである。この図から、成分波形であるMMG波形の存在時刻(MMG波形の時刻0)との関係性を見出すことは難しいと言える。
【0208】
なお、
図22は、レベル-8とレベル-7は低周波域のカーブに沿って、レベル-6から-1は高周波数域のカーブに沿ってそれぞれウェーブレット係数を抽出し、ウェーブレット係数集合を作成した例である。こちらも、正値係数のみを抽出して等高線表示させている。
この図から、取り落としているもの、取り落としそうなもの、誤検出しそうなものなどの問題も見受けられるのは確かではあるが、単純な冗長ウェーブレット解析の結果に比べると、場合によっては更なるチューニングを行う等の処理を加えることで、かなり高い抽出精度を期待できる結果になっていると言える。
【0209】
[実施例2-3]
それから、
図18の曲線に沿う形でウェーブレット係数集合を設定して解析を行った。ただし、
図22からも分かるように、低周波数帯域成分は波形の重なりが多すぎて信頼できないため、高周波数帯域成分だけで適合を調べることとし、
図23に示すような適合フィルタを設定した。
【0210】
この適合フィルタは、成分波形らしさを評価するのに高周波数帯域だけを使うことを示している。しかしながら、成分波形を得る場合には、少なくとも成分波形らしさの評価時に省いた低周波数帯域成分を推定し付加する必要がある。そのため、成分波形らしさ評価するための基本ベクトルに基づいて設定した導出型調整フィルタの例を
図24に示す。
なお、ここで設定しているウェーブレット係数集合は、成分波形の顕著な特徴要素のみを抜き出した簡素なものであるので、例えば基本ベクトルを逆変換したとしても、成分波形の形状との違いは大きい。
【0211】
そして、
図25が、成分波形(inv-WT of base wavelet-set)と、成分波形のウェーブレット係数集合を逆変換した波形(pesudo MMG wave)との対比である。tims:700msまで伸びている方の波形が、成分波形(inv-WT of base wavelet-set)である(
図27,29,30についても同じ)。
この図を見ると、成分波形の特徴をある程度残してはいるものの、差異は大きいと言える。しかし、その一方、30ms付近のピークや300ms付近のピークは正確に捉えられていると言える。
【0212】
また、導出型調整フィルタ適用の例として、ある時刻で検出された成分波形の再現例を以降に示す。
図26は
図22と同じ図であるが、
図26に示す丸印の時刻(67ms,152ms,226ms,311ms,499ms,525msの6箇所)で検出された成分波形の再現例を時刻が早い順にそれぞれ
図27(A)~27(F)に示す。
【0213】
図27(A)~27(F)の凡例に書かれた「scale」の値(各グラフの右上参照)は、それぞれの位置で合成する際に成分波形の振幅にかけた重みの値であり、成分波形のグラフはその適用結果になっている。これらの例から、各位置で合成された成分波形の振幅についても良好に捉えていることがわかる。
【0214】
さらに、より多くのウェーブレット係数を使用して規定することで、導出型調整フィルタによる成分波形の再現波形をより精細なものにすることも可能である。
図28は、
図24で示した導出型調整フィルタに、更にウェーブレット係数を追加して拡張したフィルタの例である。各係数の時間位置が示されていないので、この図だけでは情報不足であるが、負値まで含めて設定している場合の参考例として示す。
【0215】
図29に、基本ベクトルに
図28に示すフィルタを適用した結果を示す。また、
図30(A)~30(F)に、先の例での各係数集合検出結果に対して
図28のフィルタを適用した結果を示す。
これらより、導出型調整フィルタの精緻化によって成分波形の再現精度が向上していることが分かる。
【0216】
[実施例3]
更に、実施例3として、実施例1で行った計測方法および実験方法で得られた表面筋電のデータに基づいて、速筋、遅筋の発揮力推定値と活動余力調整推定値を算出した。
【0217】
図31は、各計測における速筋の活動余力推定値を示すグラフである。また、
図32は、各計測における遅筋の活動余力推定値を示すグラフである。なお、
図33は、各計測における速筋と遅筋の総合での活動余力推定値を示すグラフである。
活動余力に関する評価は、観測される運動単位活動のバラつき具合という安定性の低い特徴量に基づく評価であるため、通常は活動余力推定値の変動も大きなものになる。しかしながら、例えば、
図31の点線枠で示した部分のように、筋肉が疲労した状態で高い負荷をかけると活動のバラつき具合が変動する余裕すら乏しくなり、活動余力推定値は低変動で低い値を示すようになる。一方、休憩で疲労が回復するに従って、活動余力推定値は次第に高負荷でも高い値を示せるようになり、回復後の波形が実験開始直後の波形とかなり似たものになっていることが分かる(例えば、
図31の左上のグラフおよび右上のグラフ参照)。
【0218】
また、
図34は、各計測における速筋の発揮力推定値と活動余力調整推定値とを組み合わせたものを示すグラフである。また、
図35は、各計測における遅筋の発揮力推定値と活動余力調整推定値とを組み合わせたものを示すグラフである。なお、
図36は、各計測における速筋と遅筋の総合での発揮力推定値と活動余力調整推定値とを組み合わせたものを示すグラフである。
【0219】
上の折れ線(例えば、
図34の左上のグラフの(A))が発揮力推定値と活動余力調整推定値とを足し合わせた値を、下の折れ線(例えば、
図34の左上のグラフの(B))が発揮力推定値を示している。
【0220】
スケールファクターの調整次第でかなり変わっては来るが、発揮力推定値と活動余力調整推定値との合計は、筋疲労が進行していない状況では、変動幅は大きいものの、長期的にみると比較的安定した値を中心として変動している傾向がある。しかし、筋疲労が進行すると、負荷が大きくなるほどに活動余力調整推定値の減少が著しくなり、
図34の点線枠で示した部分のように、余力が乏しい状況で筋活動が行われている部分の量が多くなっていることが分かる。
【0221】
[実施例4]
最後に、実施例4として、以下のような実験を行った。
【0222】
[計測方法]
本実験における計測方法を、以下に示す。
(1)直立姿勢を取り、ダンベルを持って肘を伸ばして真っ直ぐにおろした状態と、肘を90度に曲げて前腕を水平にした状態との間を往復する動作を反復する。
(2)上げ下げそれぞれおよそ2秒程度で、勢いを付けずにできるだけ安定した速度で反復するようにした際の上腕二頭筋の表面筋電を計測する。
(3)100秒が経過するか昇降運動ができなくなるまで運動を継続することを1回の試行とし、休憩をはさんで3回行う。
【0223】
図37に、計測された表面筋電波形を示す。上から1回目の試行結果、2回目の試行結果、3回目の試行結果をそれぞれ示す。
図37より、2回目と3回目は途中でダンベルを持ち上げられなくなって中止していることが分かる。
【0224】
また、
図38は、各計測における速筋の活動余力推定値を示すグラフである。また、
図39は、各計測における速筋の発揮力推定値と、発揮力推定値と活動余力調整推定値との合計を示すグラフである。
【0225】
一方、
図40は、各計測における遅筋の活動余力推定値を示すグラフである。また、
図41は、各計測における遅筋の発揮力推定値と、発揮力推定値と活動余力調整推定値との合計を示すグラフである。
【0226】
なお、
図42は、各計測における速筋と遅筋の総合での活動余力推定値を示すグラフである。また、
図43は、各計測における速筋と遅筋の総合での発揮力推定値と、発揮力推定値と活動余力調整推定値との合計を示すグラフである。
【0227】
図38,40,42に示すように、速筋の活動余力推定値、遅筋の活動余力推定値、および速筋と遅筋の総合の活動余力推定値のいずれも、筋疲労が進行するにつれて値が低下していっていることが分かる。2回目や3回目の試行では、休息不足で筋疲労が十分には回復しておらず、最初から筋疲労により発揮力が低下した状態、かつ、低い活動余力推定値での運動開始となり、運動継続でさらに値が低下、つまり十分な発揮力を出すことが困難になった上に余力も失い、運動を中止していることがわかる。
また、
図39,41,42の発揮力推定値と活動余力調整推定値との合計が示すように、1回目の試行の最初の頃のように筋疲労があまり進行していない状況では、発揮力推定値が下がると活動余力調整推定値が少し大きくなって、合計の値の変動量が抑制されているのに対し、筋疲労が進行すると合計の値は大きく低下し、運動開始時と比べると非常に小さい力しか発揮できない状況になっていることが分かる。
本発明は、信号波形の解析等に利用することができ、例えば、筋力増強運動を伴うリハビリテーションにおいて、解析結果に基づき適切な負荷を適切な時間かけることができたり、医療分野やスポーツ分野等においても応用することができるため、産業上有用である。