特許 7549879 信号波形解析システム及びそのプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-09-04

(45)【発行日】2024-09-12

(54)【発明の名称】信号波形解析システム及びそのプログラム

(51)【国際特許分類】

   A61B 5/389 20210101AFI20240905BHJP

【ＦＩ】

A61B5/389

【請求項の数】 8

(21)【出願番号】P 2020215284

(22)【出願日】2020-12-24

(65)【公開番号】P2022100970

(43)【公開日】2022-07-06

【審査請求日】2023-10-16

(73)【特許権者】

【識別番号】504174135

【氏名又は名称】国立大学法人九州工業大学

(74)【代理人】

【識別番号】100120086

【弁理士】

【氏名又は名称】▲高▼津 一也

(74)【代理人】

【識別番号】100090697

【弁理士】

【氏名又は名称】中前 富士男

(74)【代理人】

【識別番号】100176142

【弁理士】

【氏名又は名称】清井 洋平

(74)【代理人】

【氏名又は名称】来田 義弘

(72)【発明者】

【氏名】永井 秀利

【審査官】藤原 伸二

(56)【参考文献】

【文献】特開２０１９－２１１９３５（ＪＰ，Ａ）

【文献】米国特許出願公開第２０１５／０２１３１９１（ＵＳ，Ａ１）

【文献】永井秀利，表面筋電の特徴抽出のためのウェーブレット係数集合のシフト選択法，信学技報，Vol.119,No.391，日本，一般社団法人 電子情報通信学会，2020年01月25日，p.51-56

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ａ６１Ｂ ５／００－５／０５３８

Ａ６１Ｂ ５／０６－５／３９８

Ｇ０６Ｆ １７／１４

ＪＳＴＰｌｕｓ／ＪＭＥＤＰｌｕｓ／ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

検出された信号波形に含まれる成分波形を解析する信号波形解析システムにおいて、
前記信号波形に対してウェーブレット解析を行い、前記成分波形に関連するウェーブレット係数特性に基づき前記成分波形に含まれる各周波数成分の周波数帯域ごとのウェーブレット係数を算出する処理手段を備え、
前記処理手段は、算出した前記各周波数成分の周波数帯域ごとのウェーブレット係数をプロットして、所定範囲の大きさの該ウェーブレット係数の集合からなる係数集合体が点在した時間－周波数平面を表示することができ、前記表示の際に、前記ウェーブレット係数を前記時間－周波数平面上に配置する時間位置を調整して、１つの前記成分波形を構成する全ての前記ウェーブレット係数が、前記時間－周波数平面上で、１つの前記係数集合体に含まれるようにすることを特徴とする信号波形解析システム。

【請求項2】

請求項１記載の信号波形解析システムにおいて、前記各係数集合体は、前記時間－周波数平面の時間軸に垂直な直線で、大半が隙間なく縦断される形状を有することを特徴とする信号波形解析システム。

【請求項3】

請求項１又は２記載の信号波形解析システムにおいて、前記ウェーブレット解析を行う際に、前記成分波形が存在するか否かに関わらず前記信号波形に含まれる、前記成分波形に該当しない波形を閾値処理によって削減することを特徴とする信号波形解析システム。

【請求項4】

請求項３記載の信号波形解析システムにおいて、前記信号波形は筋電信号であり、脱力状態の生体を対象に検出した前記信号波形から求めた前記ウェーブレット係数の標準偏差を用いて各周波数帯域で行う前記閾値処理の閾値を定めることを特徴とする信号波形解析システム。

【請求項5】

検出された信号波形に含まれる成分波形を解析する信号波形解析システムにて用いられる信号波形解析プログラムにおいて、
前記信号波形に対してウェーブレット解析を行い、前記成分波形に関連するウェーブレット係数特性に基づき前記成分波形に含まれる各周波数成分の周波数帯域ごとのウェーブレット係数を算出し、算出した前記各周波数成分の周波数帯域ごとのウェーブレット係数をプロットして、所定範囲の大きさの該ウェーブレット係数の集合からなる係数集合体が点在した時間－周波数平面を表示する際に、前記ウェーブレット係数を前記時間－周波数平面上に配置する時間位置を調整して、１つの前記成分波形を構成する全ての前記ウェーブレット係数が、前記時間－周波数平面上で、１つの前記係数集合体に含まれるようにすることを特徴とする信号波形解析プログラム。

【請求項6】

請求項５記載の信号波形解析プログラムにおいて、前記各係数集合体は、前記時間－周波数平面の時間軸に垂直な直線で、大半が隙間なく縦断される形状を有することを特徴とする信号波形解析プログラム。

【請求項7】

請求項５又は６記載の信号波形解析プログラムにおいて、前記ウェーブレット解析を行う際に、前記成分波形が存在するか否かに関わらず前記信号波形に含まれる、前記成分波形に該当しない波形を閾値処理によって削減することを特徴とする信号波形解析プログラム。

【請求項8】

請求項７記載の信号波形解析プログラムにおいて、前記信号波形は筋電信号であり、脱力状態の生体を対象に検出した前記信号波形から求めた前記ウェーブレット係数の標準偏差を用いて各周波数帯域で行う前記閾値処理の閾値を定めることを特徴とする信号波形解析プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、発生頻度が不規則な単発波からなる混合信号波形等を分析可能な信号波形解析システム及びそのプログラムに関する。

【背景技術】

【0002】

一つの筋の活動は、筋を構成する多数の運動単位の活動の合成からなる。そのため、詳細な筋活動を捉えるためには個々の運動単位の活動を得ることが望ましく、運動単位の活動を得るには、それによって生じる電位変化を捉えるのが近道となる。電位変化を捉える方法として、電極針を筋に刺入する針筋電が挙げられる。針筋電は、電極針の刺入位置付近の運動単位活動電位を直接的に計測できる。

【0003】

しかしながら、電極針の刺入は痛みを伴い、運動中の適用は困難であることから、使用するにあたり制限が多い。更に、針筋電で得られるのは電極針の刺入位置近辺の局所的な情報であり、広範囲の情報を得るには多数の針電極の使用が必要となるが、多数の針電極の刺入は現実的ではなく、運動単位を含む筋全体の活動に関しては「木を見て森を見ず」の状態となることも問題である。

【0004】

これに対して、皮膚表面に電極を貼り付けるだけの非侵襲な計測の表面筋電は、取り扱いが容易で広範囲の情報を捉えることが可能であるが、複数の運動単位活動電位が重畳したものを計測することとなる。
ここで、個々の運動単位の活動は、神経インパルスが到達した際に単発で生じることから、運動単位活動電位変化も離散的な単発波となる。

【0005】

また、一つの筋を構成する多数の運動単位は均質に活動するわけではなく、活動速度が異なる複数のタイプの運動単位が負荷に応じて選択的に投入される。加えて、疲労によって運動単位の反応は変化するため、たとえ筋の長さや負荷に変化が無かったとしても、運動単位活動電位の単発波の形状や周期等が規則性や安定性を維持できる時間は極めて短い。

【0006】

従って、表面筋電で検出された重畳した信号波形から個々の運動単位活動電位を分離抽出するのは困難である。筋活動に限らず、生体活動は必ずしも規則性を持たない単発活動の集合体であることが多く、観測される生体信号も単発活動に依存した信号の重畳である可能性が高いと考えられる。不規則な単発波形（例えば異常信号）が重畳した信号波形として観測され得る分野は、生体活動を分析する分野に限定されず、単発波形の時間位置（発生時刻）等の抽出や分析は様々な分野で求められる。

【0007】

このような信号を解析するためには局所的な特性に敏感である必要があるため、時間解像度が高いウェーブレット解析が適している。ウェーブレット解析を用いれば、単発波が固定された周波数で完全に離散的に生じるようなものを容易に検出できる。

【先行技術文献】

【特許文献】

【0008】

【文献】特開２０１８－０４７２３０号公報

【文献】特開２００２－２７２６９２号公報

【文献】特開２００７－２０２６１２号公報

【文献】特開２０１６－０６３９９５号公報

【文献】ＷＯ２００４／０２３９９６号公報

【文献】特公平０５－０３２０５７号公報

【文献】特開２０１３－２４４０２７号公報

【発明の概要】

【発明が解決しようとする課題】

【0009】

しかしながら、単発波が固定された周波数で完全に離散的に生じるという性質が期待できない波形は、筋電に含まれる運動単位活動の波形を含めて多く存在すると考えられる。そのような波形は、単発波の信号成分が時間的、周波数的に広がりを持って存在する。そのため、特許文献１、２、３で開示されているように、連続ウェーブレット変換結果を２次元又は３次元的に表示する場合、どの部分が単発波の信号成分の一部でどの部分がノイズ的な成分かといった判断が難しい。

【0010】

計算機で信号解析を行う場合に用いられるのは、サンプリングによって得られた離散データであり、離散データに対するウェーブレット解析は離散ウェーブレット解析と言われる。離散ウェーブレット解析で一般的な多重解像度解析結果を単純に用いる場合、時間－周波数平面の分割に関する制約により、時間－周波数平面上での単発波の特定は極めて困難となる。特許文献４に記載されているように、冗長離散ウェーブレット解析を行えば、多重解像度解析の問題を回避できるが、連続ウェーブレット変換と類似した問題が残る。
これに関し、単一チャネル信号から要素信号（例えば、単チャネルの表面筋電信号から構成要素たる運動単位活動電位）を検出、分離可能にすべく、多チャネルの信号とモデルないしは特別な計測条件とを用いる手法が存在する（特許文献５、６参照）。

【0011】

しかしながら、これらの手法は複雑な計算を必要とし、ノイズ等のモデルの条件に適合しない要素が混入した信号に対して、要素信号の推定精度が低下する。さらに深刻な問題は、推定の際に多チャネルの特徴情報を使うので、得られた要素信号に基づく分析を多チャネルの観点から行えないことである。単一チャネル信号からの抽出、分離が実現できれば、多チャネルを活用した分析を目指す技術（例えば、特許文献７）はさらに発展できる。

【0012】

本発明は、かかる事情に鑑みてなされたもので、成分波形の周波数成分が時間変動を伴う等単純でない信号波形であっても、成分波形が存在する時刻に対応して、その成分波形に含まれる各周波数成分を定量的に抽出可能な信号波形解析システム及びそのプログラムを提供することを目的とする。

【課題を解決するための手段】

【0013】

前記目的に沿う第１の発明に係る信号波形解析システムは、検出された信号波形に含まれる成分波形を解析する信号波形解析システムにおいて、前記信号波形に対してウェーブレット解析を行い、解析対象の前記成分波形に関連するウェーブレット係数特性に基づき前記成分波形に含まれる各周波数成分の周波数帯域ごとのウェーブレット係数を算出する処理手段を備え、前記処理手段は、前記各周波数成分の周波数帯域ごとのウェーブレット係数を、前記信号波形に前記成分波形が存在する時刻に対応して出力することを特徴としたものである。

【0014】

前記目的に沿う第２の発明に係る信号波形解析プログラムは、検出された信号波形に含まれる成分波形を解析する信号波形解析システムにて用いられる信号波形解析プログラムにおいて、前記信号波形に対してウェーブレット解析を行い、解析対象の前記成分波形に関連するウェーブレット係数特性に基づき前記成分波形に含まれる各周波数成分の周波数帯域ごとのウェーブレット係数を算出するにあたり、前記各周波数成分の周波数帯域ごとのウェーブレット係数を、前記信号波形に前記成分波形が存在する時刻に対応して出力することを特徴としたものである。

【0015】

なお、ウェーブレット解析の手法としては、計算機での処理に向き、特に、リアルタイムの解析を必要とする場合は、冗長離散ウェーブレット解析が好ましい。以降の本発明の実施の形態においては、冗長離散ウェーブレット解析を中心に説明するが、必要な時間－周波数領域のウェーブレット係数を得られればよく、第１の発明に係る信号波形解析システム及び第２の発明に係る信号波形解析プログラムに採用できるウェーブレット解析は冗長離散ウェーブレット解析に限定されない。

【0016】

第１の発明に係る信号波形解析システム及び第２の発明に係る信号波形解析プログラムが取り扱う、信号波形の成分波形は、発生頻度が不規則な単発波又はそれに類する波であることが多く、一般的に、複数の周波数成分によって構成されており、各周波数成分が存在する時間帯は周波数帯域によって異なると考えるのが自然である。成分波形の存在時間は短く、成分波形に含まれる周期の数も極めて少ないと考えられるため、検出対象に周期性を求めるフーリエ解析は、短時間フーリエ解析の手法を含めて適さない。よって、単一の信号波形の成分波形を解析するには、フーリエ解析の利用に比べ、ウェーブレット解析の利用が遥かに望ましい。

【0017】

ウェーブレット解析では、使用するウェーブレットの形状に近い形状の信号変化に対して敏感に反応することから、成分波形の周波数成分が存在する時間帯が周波数ごとに異なる場合、周波数成分に反応して大きなウェーブレット係数が出力される時間帯が周波数成分ごとに異なる。そのため、成分波形の発生頻度が低く、成分波形が解析対象の信号中に孤立して存在する場合であれば、ウェーブレット解析から得られる値を、例えば時間－周波数平面上で表すと時間的及び周波数的に広がりを持った大きなウェーブレット係数値の塊として成分波形が存在することとなる。

【0018】

しかしながら、多数の成分波形が時間的に近接して解析対象信号に含まれる場合、ウェーブレット解析結果に基づく信号強度分布を時間－周波数平面上に単純に表示すると、例えば、それぞれ異なる成分波形に属する２つのウェーブレット係数が、時間－周波数平面上で同一時間位置の高周波域と低周波域とに現れる。その結果、時間－周波数平面上のどの領域にどのような信号強度の塊が存在しているのかが不明瞭となり、成分波形の検出は極めて困難になる。

【0019】

一方、成分波形を構成する周波数成分と当該周波数成分が存在する時間帯とを考えると、一つの成分波形を構成する周波数成分として強い関係を持つウェーブレット係数（このウェーブレット係数の群を、以下、「ウェーブレット係数群」とも言う）は、時間－周波数平面上で規則性を持って分布すると考えられる。なお、当該ウェーブレット係数群の各ウェーブレット係数間には、ウェーブレット解析における計算上の直接的な関係や時間帯の包含関係が存在するとは限らないし、そうした関係が存在する必要もない。

【0020】

ここで、ウェーブレット係数群について、時間－周波数領域において、各ウェーブレット係数の時間位置と周波数帯域との関係を規定したものを、ウェーブレット係数特性としてウェーブレット係数集合と定義し、このウェーブレット係数集合の時間－周波数領域における基準時間位置を基準時点と定義する。ウェーブレット係数集合で捉えようとする成分波形は、ウェーブレット係数集合と同様に時間幅があり、ウェーブレット係数集合の基準時点は、成分波形が存在する時刻でもある。

【0021】

これにより、単一チャネルの解析対象信号であっても、この解析対象信号をウェーブレット解析した結果を、例えば時間－周波数平面上で表した場合の各時刻において、当該時刻を基準時点とするウェーブレット係数集合に基づいてウェーブレット係数群を抽出すれば、当該時刻に存在する成分波形を抽出及び分析することができる。つまり、ウェーブレット係数集合に対応するウェーブレット係数を全て、各周波数帯域毎に当該基準時点に基づいて調整し配置すれば、前述したように，成分波形に含まれる各周波数成分は，周波数帯域とそのウェーブレット係数の組からなる群として，それぞれが異なる時刻に出力されることなく，成分波形が存在する時刻に対応して出力される。また、この結果を時間-周波数平面上に表示すると、上記時刻に対応して、時間軸に垂直な方向に値の大きなウェーブレット係数の並びが出現することになり、解析対象信号中の成分波形の分布の視認や分析を容易に行うことができる。

【0022】

解析対象信号を計算機に取り込む際のサンプリング周波数は、通常、解析対象信号に含まれる周波数成分の最大周波数値に基づいて定めるが、本発明に係る成分波形の分析では、それに加えて、成分波形がどれだけ時間的に近接して含まれているか、それらをどれだけの時間解像度で抽出したいかを考慮して、サンプリング周波数を調整する。分離したい時間幅に対し、それを周期とする周波数の２倍以上のサンプリング周波数とすることが好適である。

【発明の効果】

【0023】

第１の発明に係る信号波形解析システム及び第２の発明に係る信号波形解析プログラムは、信号波形に対してウェーブレット解析を行い、解析対象の成分波形に関連するウェーブレット係数特性に基づき成分波形に含まれる各周波数成分の周波数帯域ごとのウェーブレット係数を算出する処理手段を備え、処理手段は、この各周波数成分の周波数帯域ごとのウェーブレット係数を、信号波形に成分波形が存在する時刻に対応して出力するようにするので、成分波形の周波数成分が時間変動を伴う等単純でない信号波形であっても、成分波形が存在する時刻に対応して、その成分波形に含まれる各周波数成分を定量的に抽出することができる。また、例えば、時間－周波数平面上で表すと、解析対象信号での成分波形の時間的挙動や個々の成分波形を構成する周波数成分及び大きさを推定、分析可能となる。なお、周波数成分の時間変動の特性に関する情報は、ウェーブレット係数の配置調整設定に包含されていることは言うまでもない。

【0024】

また、第１の発明及び第２の発明により、従来手法に見られる、多チャンネルでの計測や、又はモデルや特別な計測条件等との併用は特に必要なく、単一のチャンネルから計測される信号により、詳細な成分波形の分布、分析が可能となる。

【図面の簡単な説明】

【0025】

【図1】本発明の一実施の形態に係る信号波形解析システムの説明図である。

【図2】時間－周波数平面の説明図である。

【図3】（Ａ）、（Ｂ）はそれぞれ、ウェーブレット解析の説明図である。

【図4】（Ａ）、（Ｂ）はそれぞれ、ウェーブレット解析の説明図である。

【図5】ウェーブレット解析の説明図である。

【図6】（Ａ）、（Ｂ）、（Ｃ）はそれぞれ、ウェーブレット係数集合の選択例である。

【図7】（Ａ）、（Ｂ）、（Ｃ）はそれぞれ、ウェーブレット係数集合の選択例である。

【図8】（Ａ）、（Ｂ）はそれぞれ、加重が約５０％ＭＶＣの状態の表面筋電信号を基に従来手法を適用して得た時間－周波数平面である。

【図9】（Ａ）、（Ｂ）はそれぞれ、加重が約９０％ＭＶＣの状態の表面筋電信号を基に従来手法を適用して得た時間－周波数平面である。

【図10】（Ａ）、（Ｂ）はそれぞれ、図８（Ａ）、（Ｂ）の時間－周波数平面を等高線化した時間－周波数平面である。

【図11】（Ａ）、（Ｂ）はそれぞれ、図９（Ａ）、（Ｂ）の時間－周波数平面を等高線化した時間－周波数平面である。

【図12】（Ａ）、（Ｂ）はそれぞれ、加重が約５０％ＭＶＣの状態の表面筋電信号を基にＴ＋１／２のウェーブレット係数集合を用いて得た等高線化された時間－周波数平面である。

【図13】（Ａ）、（Ｂ）はそれぞれ、加重が約９０％ＭＶＣの状態の表面筋電信号を基にＴ＋１／２のウェーブレット係数集合を用いて得た等高線化された時間－周波数平面である。

【図14】図１３（Ａ）の一部拡大図である。

【図15】（Ａ）～（Ｆ）はそれぞれ、時間－周波数平面上の所定の時刻のウェーブレット係数を逆変換して得た波形である。

【発明を実施するための形態】

【0026】

続いて、添付した図面を参照しつつ、本発明を具体化した実施の形態につき説明し、本発明の理解に供する。
図１に示すように、本発明の一実施の形態に係る信号波形解析システム１０は、人体に取り付けられて、表面筋電信号（波形信号の一例）を検出するセンサ部１１と、センサ部１１から出力されるアナログの表面筋電信号を増幅するアンプ１２と、増幅されたアナログの表面筋電信号をデジタル化するＡ／Ｄ変換器１３と、デジタル化された表面筋電信号（「デジタル化された表面筋電信号」を、以下、単に「表面筋電信号」とも言う）を対象に冗長離散ウェーブレット解析等を行って表面筋電信号に含まれる成分波形を解析する処理手段１４とを備えている。

【0027】

センサ部１１は、表面筋電信号の成分波形を捉えるのに十分な時間解像度が得られるようなサンプリング周波数で計測するように設定されている。なお、センサ部１１は、計測対象者の皮膚に対して接触するタイプであってもよいし、接触しないタイプであってもよい。表面筋電信号は、従来、サンプリング周波数を１～３ｋＨｚとするのに対し、本実施の形態では、成分波形抽出の時間解像度を高めることを主な目的として、サンプリング周波数を従来の値より高い値（例えば、１０～２０ｋＨｚかそれ以上）にしている。
表面筋電信号の検出時にデジタルフィルタ等によりノイズ低減処理を行ってもよいが、成分波形の抽出には成分波形に含まれる最も高い周波数成分が分解能に影響を及ぼすため、ノイズ低減処理を行う場合は、表面筋電信号の波形に残す高周波成分にできるだけ余裕を残すように行うのが好ましい。

【0028】

処理手段１４は、例えば、モニタ付きのコンピュータによって構成可能であり、ウェーブレット係数集合を得るための信号波形解析プログラムがインストールされている。処理手段１４は、Ａ／Ｄ変換器１３から入力される表面筋電信号に対し、冗長離散ウェーブレット解析（採用するウェーブレットには、一般的な多重解像度解析で用いられる種類を選ぶことが可能）、あるいは、それに準じて特徴抽出に必要な時間－周波数領域のウェーブレット係数を取得可能なウェーブレット解析を行って、成分波形に含まれる各周波数成分の周波数帯域に対応するウェーブレット係数を求めることができる。また、当該ウェーブレット係数をプロットし、時間－周波数平面上にグラフとして表示することも可能である（図２参照）。このような構成を有する処理手段１４により、波形信号からサンプリングされた値を基にウェーブレット解析を行って、各周波数帯域に対応するウェーブレット係数が求められ、また解析結果は必要に応じて時間－周波数平面上に表示できることとなる。

【0029】

時間－周波数平面は、図２に示すように、縦軸を周波数軸とし、横軸を時間軸とする平面であって、平面上の各単位領域にウェーブレット係数が入力されている。なお、図２は、各周波数帯域における、ウェーブレット係数の単位領域がわかりやすいように多重解像度解析の結果を示しているが、冗長ウェーブレット解析では１サンプル時刻ずつずらしたすべての単位領域でのウェーブレット係数を求めている。時間－周波数平面では、所定範囲の大きさ（例えば、絶対値が所定値以上）のウェーブレット係数の集合からなる係数集合体（塊）が点在している。図２に示された時間－周波数平面は、色の濃淡でウェーブレット係数の大きさを表現しており、色が濃い単位領域ほどウェーブレット係数が大きいことを意味する。

【0030】

以下に、先ずは波形信号からサンプリングされるデータ（以下、「サンプリング値」とも言う）に基づき時間－周波数平面の単位領域にプロットされるウェーブレット係数を求める処理について説明する。なお、以下では、Ｄａｕｂｅｃｈｉｅｓ’Ｎ＝２のウェーブレットを採用し、周波数帯域の数を３つとした場合を説明する。
ウェーブレット係数の算出には、図３（Ａ）に示すように、サンプリング値が入力されるテーブル１６と、それぞれ時間－周波数平面の異なる周波数帯域に対応するテーブル１７、１８、１９とが利用される。各テーブル１６、１７、１８、１９は、それぞれ数値を格納可能な複数のバッファを有している。
なお、各テーブルの横方向に示される、０、－1、－2等の番号はサンプル値が入力される時間（以降、時間位置ともいう。）を示す。マイナスの符号は過去に入力されたデータであることを意味し、例えば、「－１」のデータは、後述する流れに沿って０番目に入力されているデータの前に入力されたものである。またテーブル番号として示される１６、１７、１８は周波数帯域を示し、番号の値が大きくなるほど、低周波数であることを示す。

【0031】

Ａ／Ｄ変換器１３から新たにサンプリング値が処理手段１４に与えられると、当該サンプリング値は、図３（Ａ）に示すように、先ずテーブル１６の０番目のバッファに入力され、次のサンプル値が０番目に入力される毎に、その前にバッファに入力されていたサンプル値が左へ順次移動することで次のバッファに入力される。そして、図３（Ｂ）に示すように、テーブル１６の０番目～－３番目の各バッファにサンプル値が入力されると、この４つの値（幅４で定義されるＤａｕｂｅｃｈｉｅｓ’Ｎ＝２の場合）からウェーブレット係数と、そのウェーブレット係数で示される信号成分を取り除いた残りの低周波信号成分の値（以下、単に「残存信号成分」とも言う）が算出される。なお、テーブル１６の－１番目～－３番目のバッファにはそれぞれ、１個前のサンプリング値、２個前のサンプリング値及び３個前のサンプリング値が入力されている。
そして、算出されたウェーブレット係数は、テーブル１７の－３番目のバッファに入力され、当該ウェーブレット係数と同時に算出された残存信号成分の値がテーブル１６の－３番目のバッファに上書きされる。

【0032】

上記のフローで解析が進むと、次に、図４（Ａ）に示すように、テーブル１６の－３番目のバッファ、－５番目のバッファ、－７番目のバッファ及び－９番目のバッファにそれぞれ入力されている４つの値からウェーブレット係数と残存信号成分が算出され、算出されたウェーブレット係数がテーブル１８の－９番目のバッファに入力され、同時に算出された残存信号成分の値がテーブル１６の－９番目のバッファに上書きされる。なお、上書き前の値については、上記フローと同様に解析が進められているので、テーブル１６の－９番目のバッファには、テーブル１７のバッファに入力される値として６個前に算出されたウェーブレット係数と同時に算出された残存信号成分の値が入力されていたことになる。

【0033】

同様にして、図４（Ｂ）に示すように、テーブル１６の－９番目のバッファ、－１３番目のバッファ、－１７番目のバッファ及び－２１番目のバッファにそれぞれ入力されている４つの値からウェーブレット係数が算出され、算出されたウェーブレット係数がテーブル１９の－２１番目のバッファに入力され、そのウェーブレット係数と同時に算出された残存信号成分の値がテーブル１６の－２１番目のバッファに上書きされる。なお、テーブル１６の－２１番目のバッファには、上書き前は、テーブル１８のバッファに入力される値として１２個前に算出されたウェーブレット係数と同時に算出された残存信号成分の値が入力されていたことになる。

【0034】

そして、図５に示すように、テーブル１７の－２１番目のバッファ、テーブル１８の－２１番目のバッファ及びテーブル１９の－２１番目のバッファにそれぞれ入力されているウェーブレット係数が、時間－周波数平面の単位領域にプロットされるウェーブレット係数として出力され、テーブル１６、１７、１８、１９それぞれの－２１番目のバッファに入力されていた数値が消去される。

【0035】

テーブル１７、１８、１９に関して、前述したように、テーブル１７が時間－周波数平面の最も高い周波数帯域に対応し、テーブル１８が時間－周波数平面の２番目に高い周波数帯域に対応し、テーブル１９が時間－周波数平面の最も低い周波数帯域に対応する。テーブル１７、１８、１９の－２１番目のバッファに出力されたウェーブレット係数は、時間－周波数平面の周波数帯域において、それぞれに対応する周波数単位の領域の、時間位置が等しい位置にプロット（入力）される。

【0036】

ここで、テーブル１７の各バッファに入力されているウェーブレット係数に対応する時間幅は、図５に示すように、テーブル１６のバッファ２個分の時間幅に等しく（テーブル１７のハッチングしている領域の時間幅であり）、テーブル１８の各バッファに入力されているウェーブレット係数に対応する時間幅は、テーブル１６のバッファ４個分の時間幅に等しく（テーブル１８のハッチングしている領域の時間幅であり）、テーブル１９の各バッファに入力されているウェーブレット係数に対応する時間幅は、テーブル１６のバッファ８個分の時間幅に等しい（テーブル１９のハッチングしている領域の時間幅である）。

【0037】

その後、テーブル１６、１７、１８、１９の各バッファに入力されている数値がそれぞれ１つ左のバッファに移動（例えば、テーブル１６の－２０番目のバッファに入力されている数値は、テーブル１６の－２１番目のバッファに移動）し、テーブル１６の０番目のバッファは空となって新たなサンプリング値が入力可能な状態となる。なお、テーブルをリングバッファとして実装しておけばテーブル先端位置の変更だけで良く、直接的なデータの移動は不要である。

【0038】

周波数帯域ごとに使用するウェーブレットを変更する場合には、ある帯域での算出において参照するテーブル１６のインデックス位置と入力するテーブルのインデックス位置が変更される。例えばテーブル１８のウェーブレット係数算出時に幅６で定義される種類のウェーブレットを用いる場合，テーブル１６の－３番目のバッファ、－５番目のバッファ、－７番目のバッファ、－９番目のバッファ、－１１番目のバッファ、及び－１３番目のバッファにそれぞれ入力されている６つの値（このとき図３とは異なり、－１０番目から－１３番目もテーブル１７算出時の残存信号成分）から算出したウェーブレット係数をテーブル１８の－１３番目のバッファに入力する。
処理手段１４は、上述した処理を繰り返すことによって、係数集合体が点在した時間－周波数平面を得る。

【0039】

ここで、前述したように成分波形は、一般的に、複数の周波数成分によって構成されており、また各周波数成分が存在する時間帯も周波数帯域によって異なり様々な形状が存在する。このような成分波形の分布の推定及び分析を正確に行うためには、それらの成分波形の形状に応じた解析が必要である。
そこで、処理手段１４では、先ず、分析しようとする成分波形の、時間－周波数領域におけるウェーブレット係数特性（以下、「ウェーブレット係数集合」ともいう）を解析により把握し、そのウェーブレット係数集合に基づいて解析対象とする信号波形に含まれる成分波形を分析する。つまり、成分波形の時間－周波数特性に応じて予め、把握される、上記ウェーブレット係数集合を、例えば、上述のようなテーブルとして設定、解析し、信号波形に成分波形が存在する時刻に対応して、その成分波形に含まれる周波数成分とそのウェーブレット係数を出力することができる。

【0040】

なお、上記の分析しようとする成分波形について、処理手段１４により解析するために、例えば、信号波形を解析する前に当該波形の形状等を初期条件として入力する等の設定が必要になるが、その方法に限定されるものではない。例えば、記憶手段等を設け、代表的な成分波形を予め格納しておくことや、また、既にウェーブレット係数集合が既知、又は推定可能な成分波形の場合については、そのウェーブレット係数集合を直接、入力することや、後述するような各種ウェーブレット係数集合を予め記憶手段に格納しておき、分析を行う際に、成分波形に近いウェーブレット係数集合を選択し、各周波数帯域におけるウェーブレット係数の時間位置を調整して成分波形の分析を行う等、様々な方法が想定される。

【0041】

図５に示す例では上記に示す、１つのウェーブレット係数集合の例を示しており、このウェーブレット係数集合を用いた場合は、処理手段１４は、時間－周波数平面にプロットするウェーブレット係数として、テーブル１７、１８、１９で示される各周波数帯域のそれぞれで同一の時間位置に存在する－２１番目のバッファのウェーブレット係数を選出されるように調整される。このウェーブレット係数集合は後述するＨ＋０に相当する。以下にこのようなウェーブレット係数集合について具体的に説明する。

【0042】

例えば、ウェーブレット係数集合としては、上述の図５に示す以外にも、図６（Ａ）、（Ｂ）、（Ｃ）、図7（Ａ）、（Ｂ）、（Ｃ）に示すパターンなどが想定される。テーブル１６、１７、１８、１９は、これらのウェーブレット係数集合の各領域２０のウェーブレット係数が計算により出揃うまで保持できるように設定される。ウェーブレット係数集合の各領域２０のウェーブレット係数を保持するバッファは、各領域の時間軸方向先端（左端）の時間位置に位置するバッファである。ウェーブレット係数集合の各領域２０のウェーブレット係数が出揃った時、それらウェーブレット係数は当該ウェーブレット係数集合の基準時点のウェーブレット係数の群として出力される。

【0043】

なお、ウェーブレット係数を選出する領域２０を各周波数帯域で異なるようにする方法（ウェーブレット係数集合を得る方法）の基本形として、「永井秀利、表面筋電の特徴抽出のためのウェーブレット係数集合のシフト選択法、信学技報、ｖｏｌ．１１９、ｎｏ．３９１、ＭＢＥ２０１９－７８、ｐｐ．５１－５６（２０２０）」（以下、「文献Ｐ」と言う）に記載されたシフト選択法が挙げられる。

【0044】

シフト選択法は、従来用いられているような多重解像度解析での直接的な計算上の関係を持つようなパターン（Ｈ＋０のパターン）に反して，計算上の関係を無視した配置で得たウェーブレット係数のＲＭＳ等の時刻ごとの信号特徴量で信号波形の全体的な特徴を調べることに有効性があることを確認したものである。更に、当該手法では、連続的な配置操作により、上記の全体的な特徴も連続的に変化する性質を持つことも確認した。

【0045】

この結果を踏まえ、本発明は更に検討を進め、このようなシフト選択法での、単純に無根拠に配置を変更して信号波形の全体特性を比較分析するのではなく、当該信号波形に含まれる個々の成分波形に着目できるようにしたものである。本発明では対象とする成分波形の形状に応じた、シフト選択法の範囲に留まらないウェーブレット係数集合を検討、設定すると共に、更に当該係数集合に応じて各周波数帯域のウェーブレット係数の時間位置を調整して配置し、各時刻のウェーブレット係数値の分布を調べることで、信号波形中の成分波形の挙動を分析できること、つまり、信号波形の細部まで分析できることを見出したものである。

【0046】

なお、本発明ではシフト選択法が持つ連続性に基づき、あるウェーブレット係数集合に対しても、そのウェーブレット係数集合を基本形としてシフト選択法での配置操作手法を利用し、各周波数帯域の領域を少しずらしながら比較、分析することも可能である。そのようにすることで例えば、周波数帯域間のスケーリング関係により、例えば、同種形状で波長が異なるような成分波形は、ウェーブレット係数集合の周波数軸上での分布が高周波帯域寄りか低周波帯域寄りかという違いで現れる。

【0047】

図６（Ａ）は、ウェーブレット係数集合を構成する領域２０の時間位置が各周波数帯域で１／２シフトしている例であり、図６（Ｂ）は同時間位置が各周波数帯域で１／４シフトしている例であり、図６（Ｃ）はシフトが無い例である。以下、各領域２０のシフト量をＳとして、ウェーブレット係数の選出される領域２０の時間位置が各周波数帯域で異なる本実施の形態を、シフト量を含むＴ＋Ｓとして表現する。図６（Ａ）、（Ｂ）、（Ｃ）に示す例は、それぞれＴ＋１／２、Ｔ＋１／４、Ｔ（Ｔ＋０）として表される。

【0048】

なお、Ｔは、ウェーブレット係数集合の１つのパターンを示すものであり、各周波数帯域の領域２０の時間軸方向末端（右側）の位置が揃ったウェーブレット係数集合を表す。またＴ＋Ｓは、そのＴのパターンを基に、そのパターンの最も低い周波数帯域の末端を基準とし、この末端に対して隣接する高周波数側の高周波数帯域の末端（右端）をＳのみシフトさせ、その更に高周波数帯域も、隣接する帯域間で同様に順次、シフトさせて形成できる。

【0049】

また、図７（Ａ）、（Ｂ）に示すウェーブレット係数集合は、各周波数帯域の領域２０の時間軸方向先端（左端）を揃えたものを基準したＨ に対して、各周波数帯域の時間位置をシフトさせたものを示している。図７（Ａ）、（Ｂ）、（Ｃ）に示す例は、それぞれＨ－１／２、Ｈ－１／４、Ｈ（Ｈ＋０）であり、Ｈのパターンについては、従来技術で示したものである（図７（Ｃ）参照）。

【0050】

本発明による成分波形の分析においては、分析が必要なターゲットとする成分波形に応じたウェーブレット係数集合を用いて行う。具体的には例えば、高周波数域で強いピークを持ちつつ、周波数が減少しつつ減衰、消失するような性質の成分波形については、シフト量はその減衰の速度に関連することを加味し、Ｈ－１．０やＨ－０．５等のＨ－Ｓの形態に基づいて、成分波形の挙動を明確に把握することが有効である。

【0051】

一方、解析対象信号に、Ｄａｕｂｅｃｈｉｅｓ‘ Ｎ＝２のウェーブレットで示されるような、一つの主たるピークを持った周波数成分を含む成分波形の場合は、末端位置揃え（Ｔ）を標準設定として０以上のシフト量を設定したウェーブレット係数集合が有用であり、シフト量を調整する事で立ち上がりと減衰の速度をバランスよく設定し、成分波形の安定した分析が可能となる。

【0052】

またこのように予め、ターゲットとする成分波形の特性が推定可能な場合は、前述にも少し触れたように、記憶手段を設け、上記のＴやＨのパターンを格納しておき、例えば、成分波形の解析によるウェーブレット係数集合との併用、又は成分波形を予め、解析することなく、後述するようにこれらのパターンに基づいて各周波数帯域の時間位置を適宜、細かく調整して分析を行うこと等も可能である。

【0053】

ウェーブレット係数には正の値及び負の値が存在するところ、本実施の形態では、選出したウェーブレット係数の逆変換波形が成分波形の形状に近いほど各ウェーブレット係数が正の大きな値に揃うようにウェーブレット係数集合を構成する領域を定めている。逆に、成分波形の形状に近いほど負の大きな値を取りやすいのであれば、負の値に揃うように領域を定める。これは、プロットしたウェーブレット係数集合体が明瞭になるようにするためである。なお、正又は負に揃えるために、一部の周波数領域だけでウェーブレット係数の値の正負を反転させる処理を行ってもよい。

【0054】

また、成分波形の時間幅に変動があるような場合は、スケール性を考慮して、ウェーブレットの形状が単独で成分波形に近いウェーブレットの種類を選ぶようにするのが好ましい。例えば、表面筋電とその構成要素である運動単位活動電位波形とを対象にする場合、サポートが狭く運動単位活動電位波形に形状が比較的近いＤａｕｂｅｃｈｉｅｓ’Ｎ＝２がウェーブレットとして推奨される。仮に成分波形が周波数帯域ごとに異なる波形特徴を有することが判明している場合、周波数帯域により利用するウェーブレットの種類を変えてもよい。

【0055】

ウェーブレット係数集合における各周波数帯域をシフトさせながら、選出するウェーブレット係数の時間位置を周波数帯域ごとに調整するにあたり、係数集合体を構成するウェーブレット係数の逆変換波形を成分波形全体の形状に近付くようにしなくてもよい。例えば、成分波形の後半の形状のみに近付けるような選択をし、成分波形の後半を検出した上で成分波形全体を推測するような処理を行うことができる。但し、この場合、対象外とした成分波形の部分から得られるウェーブレット係数が抽出精度を低下させる可能性があることに注意が必要である。

【0056】

また、本実施の形態では、成分波形を明瞭化するため、ウェーブレット解析時に閾値処理（定常波低減処理又はウェーブレット縮退処理）を行って定常波を削減する。これは、成分波形がほとんど含まれない信号波形を定常波と扱い、解析対象の信号波形から定常波の信号成分を削減するものであり、定常波をノイズと見なしたノイズ低減処理に相当する。本実施の形態では、定常波を削減する閾値処理に、ｓｏｆｔ ｔｈｒｅｓｈｏｌｄｉｎｇを採用しているが、その他を採用してもよい。なお、閾値処理は、ウェーブレット解析時ではなく、ウェーブレット解析後にテーブルの所定のバッファから出力されたウェーブレット係数に対して行うようにしてもよい。

【0057】

閾値処理の閾値を定める手法は数多く存在する。本実施の形態のように信号波形が表面筋電信号の場合、被験者（即ち、生体）が脱力状態の信号を定常波としてウェーブレット解析し、得られた信号成分のほぼ全てをノイズであるとみなして削減するように閾値を定めることが多い。この際、全ての周波数帯域で閾値を同一にする方法と、周波数帯域ごとに閾値を個別に定める方法とに大きく分けられるが、表面筋電信号を対象とする場合、低周波帯域と高周波帯域とで信号強度（ウェーブレット係数の大きさ）にかなりの違いがある（低周波帯域の信号強度が高周波帯域の信号強度に比べてかなり大きい）ため、全ての周波数帯域で閾値を同一にすると、低周波ノイズが多く残ってしまったり、高周波成分を削り過ぎて情報が失われたりする。

【0058】

一方、各周波数帯域のウェーブレット係数値に基づいて周波数帯域ごとに独立して閾値を定めると、時間－周波数平面上に乱雑な信号成分残存が発生し易く、分析の大きな妨げとなる。そこで、本実施の形態では、パラメータで効果を調整可能な以下に記す手順で周波数帯域ごとの閾値を定めるようにしている。

【0059】

（１）脱力状態の生体を対象に検出された表面筋電信号に対し、ウェーブレット解析を行ってウェーブレット係数を求める。これにより各周波数帯域において得られたウェーブレット係数の個数をｎとする。
（２）得られたウェーブレット係数には例外的に極度に大きなノイズやノイズ以外の信号成分が含まれている可能性があるため、パラメータとして外れ値割合基準ＯＢを与え、周波数帯域ごとのウェーブレット係数について、大きい順に上位１００×（ＯＢ／ｌｏｇ（ｎ））％を外れ値として除外する。ＯＢの適切な値は計測環境等に依存する。一例として、本実施の形態では、ＯＢ＝０．２５とした。

【0060】

（３）周波数帯域ｉに対し、平均を０として各周波数帯域のウェーブレット係数の標準偏差ＳＤｉを求める。ここで、全周波数帯域のＳＤｉの最大値をｍａｘ＿ＳＤとする。
（４）閾値重みＴＷをパラメータとして与え、例えば、以下の式１、式２のいずれかを周波数帯域ｉの閾値とする。

【0061】

ＳＤｉ×ＴＷ×（１＋ｌｏｇ（ｍａｘ＿ＳＤ／ＳＤｉ）） ・・・式１
ＳＤｉ×ＴＷ×（１＋√（ｍａｘ＿ＳＤ／ＳＤｉ）） ・・・式２

【0062】

式１、式２のいずれが適切か、ＴＷの適当な値は取り扱うデータの種類によって異なる。経験的には式１のほうが式２に比べて良好な結果を示すことが多く、本実施の形態では、一例として、ＴＷ＝２．０とした。
ここで、式１、式２から、各周波数帯域で行う閾値処理の閾値は、ウェーブレット係数の標準偏差を用いて定められることとなる。

【0063】

Ｘ軸をウェーブレット係数集合の基準時刻、Ｙ軸をウェーブレット係数の周波数帯域（複数のウェーブレット係数が採用されている周波数帯域については、想定されるウェーブレット係数集合の特性に応じて仮想的に周波数帯域を分割しそれぞれに割り当てることも可能）、Ｚ軸をウェーブレット係数値（複数のウェーブレット係数が採用されている周波数帯域については、平均化処理や加算処理を行うことも可能）として、使用したウェーブレットとの適合度の観点で評価した信号強度が時間－周波数平面上で分布する様を、Ｚ軸が正値の範囲（逆変換波形が成分波形に近いほど負に大きなウェーブレット係数となるように設定した場合は、負値の範囲）で３次元プロットして視覚化することにより、Ｙ軸に平行な方向（時間軸に垂直な方向）に近い強信号領域のまとまりとして、成分波形の存在を観測することができる。なお、成分波形の観測が容易にできるように、表示するＺ軸の範囲を調整してもよいことは言うまでもない。

【0064】

解析対象の信号波形に含まれる成分波形の数が多い場合、時間－周波数平面上で値をプロットすると、低い周波数帯域で、成分波形の周波数成分とウェーブレット係数で表される強信号領域の融合が生じ易くなる。サンプリング周波数が十分であれば、高い周波数帯域で融合が生じることは少ない。このような状況では、成分波形に対して最適と考えられるウェーブレット係数集合に基づいて、時間軸方向に各領域２０を更に少しシフトさせながら調整したウェーブレット係数集合を用いて３次元プロットすることで、個々の成分波形の分布が判別し易くなることがある。融合した成分波形の周波数特性に違いがある場合には、特に効果的である。

【0065】

ウェーブレット係数集合をウェーブレット係数値が正の領域（逆変換波形が成分波形に近いほど負に大きなウェーブレット係数となるように設定した場合は、負の領域）で時間－周波数平面上に３次元プロットすることによって、個々の運動単位活動は、例えば、時間軸に略垂直な涙滴状の分布等、その運動単位活動の波形が信号波形に存在する時刻に対応して出力、表示され、観測可能となる。多くの運動単位活動が存在する場合には涙滴状分布の融合も生じるが、該運動単位活動の、より効率的な分析を行う場合は、上述したように少しずつシフトさせたり、シフト量を小さくするようにしたりすることで（例えば、Ｔ＋１／２、即ち、Ｔ＋０．５を基準設定とする場合、Ｔ＋０．３５やＴ＋０．５５等）、融合した個々の分布を分析する助けとなる。

【0066】

成分波形の分布に基づいて切り出したウェーブレット係数集合をウェーブレット逆変換することによって、信号波形中に実際に存在する成分波形の推定形状を得ることができる。単チャネルの表面筋電信号から運動単位活動を分離、抽出する場合、荒く刻むのであればシフト選択法で言うところのＴ＋１／２を基準設定の採用が好適であり、より細かく刻むのであればＴ＋６／１６ないしはＴ＋７／１６の採用が好適である。

【0067】

例えば、サンプリング周波数が２０（ｋＨｚ）の場合では、運動単位活動の涙滴状分布がピークとなる周波数帯域が大きく分けて２種類存在する。一方のピークは速筋の運動単位の活動であり、他方のピークは遅筋の運動単位の活動であると推測される。該当の時刻のウェーブレット係数集合から分布の性質に基づいて分離、抽出を行い、ウェーブレット逆変換で波形を求めることによって、速筋と遅筋それぞれの運動単位活動電位波形を推測することができる。

【0068】

各運動単位活動の涙滴状分布は時間軸方向に山形に推移する（図１３（Ａ）、（Ｂ）参照）。その幅は、観測される信号が変化する時間幅であり、信号伝達速度との関係が強いと考えることができる。筋疲労の前後で比較すると、筋疲労後には低周波数域の値が増大して時間幅が広い涙滴状分布が増えているようにも見受けられ、筋疲労時に生じる筋繊維伝導速度の低下を反映している可能性が高い。よって、涙滴状分布の分析は、筋疲労の分析に新たな観点を与えられるのが明らかである。

【0069】

このように、不規則で密な成分波形で構成されており、成分波形の解析困難な信号波形として表面筋電信号に係る信号解析を例に述べてきたが、本発明は表面筋電信号に限定されるものではない。本発明の重要な点は、対象とする成分波形の出現に非常に強く反応するという点である。加えて、解析自体はリアルタイムに実施できることから、非常に低頻度で一過性の単発波形の出現監視のような用途にも使うことができる。

【0070】

この場合、解析時のテーブルを大きくして保持時間を長くしておけば、単発波形発生の一定時間前からの分析も可能である。また、本発明における成分波形の分析能力は、成分波形に周期性や反復性がある信号波形において、その異常検出にも有効である。例えば、フーリエ解析で想定する単純なｓｉｎ波の合成では表現することが難しい心電図の波形の分析に適用することもできる。Ｒ波のピークを基準時点としたウェーブレット係数集合を特徴ベクトルとしてｃｏｓ類似度で比較する等による波形変動の分析を行うことも可能である。生体情報以外でも、例えば稼働中の機械装置において反復的に動作しているソレノイドアクチュエータの駆動時衝撃音がある信号波形中の成分波形として存在するような場合は、その駆動状態を分析することにも適用可能である。汚れ等により幾分かの駆動速度の低下が生じれば、それによる動作タイミングのずれや駆動時衝撃音の高周波成分の低下が生じる可能性が高く、その検出に本発明が有効に機能すると考える。

【実施例】

【0071】

次に、本発明の作用効果を確認するために行った実験について説明する。
実験では、直立姿勢で立った状態で肘を９０度屈曲して前腕を水平にした状態で手首に加重した際の上腕二頭筋の単チャネルの表面筋電信号を計測し、当該表面筋電信号に含まれる運動単位活動の可視化を行った。実験では、サンプリング周波数が２０（ｋＨｚ）であり、使用したウェーブレット解析法はＤａｕｂｅｃｈｉｅｓ‘ Ｎ＝２の冗長離散ウェーブレット解析であった。

【0072】

手首への加重が約５０％ＭＶＣ（最大随意筋力）の状態及び約９０％ＭＶＣの状態について、従来手法（図７（Ｃ）に示すＨ（Ｈ＋０）のウェーブレット係数集合）を適用して得た時間－周波数平面を、図８～図１１に示す。
図８（Ａ）、（Ｂ）はそれぞれ、手首への加重が約５０％ＭＶＣにおける筋疲労前後の時間－周波数平面であり、図９（Ａ）、（Ｂ）はそれぞれ、手首への加重が約９０％ＭＶＣにおける筋疲労前後の時間－周波数平面である。図８（Ａ）、（Ｂ）、図９（Ａ）、（Ｂ）の時間－周波数平面は全て、各単位領域のウェーブレット係数の大きさが色の濃淡で表されている。図８（Ａ）、（Ｂ）、図９（Ａ）、（Ｂ）の時間－周波数平面からは、信号特徴の変化があることは確認できるものの、細部はつぶれて不明瞭であり、筋疲労前後の違いも判別困難であった。

【0073】

図１０（Ａ）、（Ｂ）はそれぞれ、図８（Ａ）、（Ｂ）の時間－周波数平面でウェーブレット係数の値が正の範囲を等高線で表した時間－周波数平面であり、図１１（Ａ）、（Ｂ）はそれぞれ、図９（Ａ）、（Ｂ）の時間－周波数平面でウェーブレット係数の値が正の範囲を等高線で表した時間－周波数平面である。図１０（Ａ）、（Ｂ）、図１１（Ａ）、（Ｂ）の時間－周波数平面から、筋疲労前後で信号特徴に変化があることは確認できるが、時間－周波数平面上に存在する係数集合体（塊）は大きく傾き、うねっているものや、島のように個別に存在するもの等、係数集合体間の繋がりにも規則性が乏しく、運動単位活動を分析できなかった。
よって、実験で計測された表面筋電信号を基に成分波形の運動単位活動を分析するために、Ｈのパターンは不適切であることが分かる。

【0074】

次に、同じ表面筋電信号に対して、Ｔ＋１／２のパターンを適用して得た時間－周波数平面でウェーブレット係数の値が正の範囲を等高線により表したものを図１２（Ａ）、（Ｂ）、図１３（Ａ）、（Ｂ）に示す。図１２（Ａ）、（Ｂ）、図１３（Ａ）、（Ｂ）の時間－周波数平面から、等高線によって表される各係数集合体が、涙滴のように垂れ下がった塊として表現され、各涙滴状の塊が各運動単位活動に相当することが確認できた。つまり、成分波形である各運動単位活動に関して、その各運動単位活動に含まれる各周波数成分の周波数帯域ごとのウェーブレット係数が、当該各運動単位活動の波形が存在する時刻に対して出力されていることが確認できる。

【0075】

また、このように、成分波形に応じたウェーブレット係数集合を用い調整することによって、１つの成分波形を構成する全てのウェーブレット係数が、時間－周波数平面上の複数の係数集合体に分かれるのではなく、１つの係数集合体に含まれるようにすることができることが確認できる。
従来の研究から、速筋の運動単位活動は高周波であり、遅筋の運動単位活動は低周波であることが確認されていることから、高周波数域に存在する係数集合体は速筋の運動単位活動に対応し、低周波数域に存在する係数集合体は遅筋の運動単位活動に対応すると考えられる。涙滴状の塊の等高線で表現された高さ（紙面に対して垂直な方向の高さ）は運動単位活動のエネルギーの大きさに対応している。

【0076】

図１４は図１３（Ａ）の一部を拡大したものである。図１５（Ａ）、（Ｂ）、（Ｃ）は、それぞれ、個々の涙滴状の塊を考慮せずに、図１４に示す時刻（１）、（２）、（３）に存在するウェーブレット係数の全てを、本来の時間－周波数関係に戻して逆変換した波形である。図１５（Ｄ）、（Ｅ）、（Ｆ）は、それぞれ、図１４に示す時刻（１）の高周波数域（縦軸の－５～－７の範囲）、時刻（２）の高～中周波数域（縦軸の－５～－８の範囲）、時刻（３）の高周波数域（縦軸の－５～－７の範囲）に存在するウェーブレット係数（涙滴状の塊において時間軸に垂直な直線で隙間なく縦断される部分のウェーブレット係数）を本来の時間－周波数関係に戻して逆変換した波形である。

【0077】

図１５（Ｄ）、（Ｅ）、（Ｆ）に示す波形は、神経インパルスの到達で運動単位が励起された時の電位変化に類似した波形であることが確認できた。これに対し、図１５（Ｂ）では運動単位励起時の電位変化に近い波形は見られるが、図１５（Ａ）、（Ｃ）では運動単位励起時の電位変化からかけ離れた波形となっている。

【0078】

Ｔ＋１／２のウェーブレット係数集合で得た運動単位活動電位の推定波形（図１５（Ｄ）、（Ｅ）、（Ｆ）に示す波形）は、ピークからの減衰時の負へのオーバーシュートが大きすぎるようにも見えるが、採用するウェーブレット係数集合のシフト量を調整（Ｔ＋１／２からＴ＋６／１６やＴ＋７／１６に）することで改善できるものと考えられる。

【0079】

以上、本発明の実施の形態を説明したが、本発明は、上記した形態に限定されるものでなく、要旨を逸脱しない条件の変更等は全て本発明の適用範囲である。
例えば、閾値処理は省略してもよい。
また、運動単位群の活動の時間－周波数平面上での分布を可視化する機能は省略可能である。

【符号の説明】

【0080】

１０：信号波形解析システム、１１：センサ部、１２：アンプ、１３：Ａ／Ｄ変換器、１４：処理手段、１６、１７、１８、１９：テーブル、２０：領域

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】