JP2016533231A - 非定常信号の機能要素への分解 - Google Patents

Abstract

信号処理、特に、生理学的現象において代表的な非定常シグナルを機能要素に分解する方法であって、各々の上記機能要素を特徴づけるパラメータを推定する方法。本方法は、信号の時間経過において動的トレンド及び基線のトレンドを推定する手段、非定常信号を機能要素が発達するセグメントに分ける手段、及びＥＣＫを形成するための前述の機能要素に適合するＱＧＫからの重複を補正する手段、を利用する。補正は、ＥＣＫを周波数領域に変換する手段、ＱＧＫの加重及び分散のパラメータを推定する手段、パラメータ信頼性の検証、ＱＧＫの開始時間を推定する手段、各々の再帰サイクル後の非定常信号モデルを拡張する手段、好ましくないＱＧＫを除去し、かつ好ましいＱＧＫを再構成する手段、及び予め定義した階級に属するパラメータを用いて、非定常信号の部分的モデルを作成する手段を含む。

Description

発明の詳細な説明

〔技術分野〕
本発明は、全体的に、信号処理、及び、特に非定常生物医学的信号の機能要素への分解、及び、上記機能要素の各々を特徴づけるパラメータの推定に関する。

〔背景技術〕
非定常信号は、時間に依存する特性値を伴う系によって発生する。理論的立場、実験的立場及び診断的立場から上記非定常信号を理解し、分析する必要性が、生物医学的研究のほぼすべての分野において、生じている。しかし、非定常信号は、一般信号理論に従わず、また、非定常信号を扱う好ましい数学的手順は、かなり限られている。

現在の信号処理手順は、一般的に、線形定常系の理論に基づいている。従って、実データの非定常特性は、しばしば無視されるか、又は上記効果は無視できるものとみなされる。時間的不変性及び線形性に関する当該仮説によって、いくつかの応用例に関して許容できる結果が提供されているが、多くの生理現象は、極めて可変的及び非線形的であり、従って、理にかなった正確な線形近似が認められていない。一方で、膨大な数の科学的及び診断的証拠により、心電図（ＥＣＧ）及び脳波図（ＥＥＧ）のような、主要な非定常生理信号に関する、より総合的な解析は、異なる正常条件及び病的条件下における、様々な生理学的システムのメカニズム及び状態に関する、より深い理解を導き得ることが示唆される。このような背景から、ヒトの脳及び体のイメージング、コンピュータを導入した健康モニタリング、並びにＥＥＧベースの脳−コンピュータ間インターフェースといった、ますます挑戦的な多数の応用例への実現可能なアプローチとして、非定常信号処理に関する新たなアルゴリズムの開発が浮上している。同時に、コンピュータを用いた生物医学的システムの継続的な進歩によって、速さ、サイズ及び費用といった点で、洗練されたアルゴリズムの開発が、実用的かつ優れたコストパフォーマンスを示すようになってきた。

非定常信号解析は、線形定常系理論の場合のような適切な理論及び計算の枠組みによって裏付けられていない。特定のタイプの生理信号の解析においては、数学的に厳格な方法よりもむしろ発見的方法が、一般的に受け入れられている。例えば、ＥＥＧ信号解析の際、上記信号は、局所定常性概念（Ｂａｒｌｏｗ， Ｊ．Ｓ． “Ｍｅｔｈｏｄｓ ｏｆ ａｎａｌｙｓｉｓ ｏｆ ｎｏｎ−ｓｔａｔｉｏｎａｒｙ ＥＥＧｓ， ｗｉｔｈ ｅｍｐｈａｓｉｓ ｏｎ ｓｅｇｍｅｎｔａｔｉｏｎ ｔｅｃｈｎｉｑｕｅｓ： ａ ｃｏｍｐａｒａｔｉｖｅ ｒｅｖｉｅｗ”． Ｊ． Ｃｌｉｎ． Ｎｅｕｒｏｐｈｙｓｉｏｌ． １９８５， ｖｏｌ． ２， ｐｐ． ２６７−３０４ 参照）を用いて解析され得る。上記概念は、信号は本質的には非定常であるが、極少時間内において、信号パターンがごくわずかに定常でなくなることを推測している（Ｐｒｉｅｓｔｌｅｙ， Ｍ．Ｂ．： “Ｎｏｎ−Ｌｉｎｅａｒ ａｎｄ Ｎｏｎ−Ｓｔａｔｉｏｎａｒｙ Ｔｉｍｅ Ｓｅｒｉｅｓ Ａｎａｌｙｓｉｓ”． Ａｃａｄｅｍｉｃ Ｐｒｅｓｓ， Ｌｏｎｄｏｎ， １９８８ 参照）。このような準定常成分は、明確な非線形系によって発生しうる信号の種々の要素と結びついている。

これまで、非定常信号の成分の検出に関する問題点には、データ駆動型信号処理方法論及びモデル駆動型信号処理方法論を用いて対処してきた。

データ駆動型信号処理方法の利点は、複雑かつ実験的な波形の動きを詳細に記述できる点である。上記データ駆動型方法の一例は、Ｈｕａｎｇらによって構築されてきた、経験的モード分解（ＥＭＤ）であり（Ｈｕａｎｇ ｅｔ ａｌ．： “Ｔｈｅ ｅｍｐｉｒｉｃａｌ ｍｏｄｅ ｄｅｃｏｍｐｏｓｉｔｉｏｎ ａｎｄ ｔｈｅ Ｈｉｌｂｅｒｔ ｓｐｅｃｔｒｕｍ ｆｏｒ ｎｏｎｌｉｎｅａｒ ａｎｄ ｎｏｎ−ｓｔａｔｉｏｎａｒｙ ｔｉｍｅ ｓｅｒｉｅｓ ａｎａｌｙｓｉｓ”． Ｐｒｏｃ． Ｒ． Ｓｏｃ． Ｌｏｎｄｏｎ， Ｓｅｒ． Ａ， １９９８， ｖｏｌ． ４５４， ｐｐ． ９０３−９９５ 参照）、適用原理の下、内因性モード関数と呼ばれる、「振る舞いがよい」振動関数の和に非定常信号を分解する。以下の開発において、各々内因性モード関数に対する瞬時周波数を提供するＨｉｌｂｅｒｔ変換によって、ＥＭＤは、裏付けられている。これらの測定を用いて、本技術は、瞬時周波数の変化という点に関して本システムを記述する。しかし、本解決策は、多くの実験的パラメータに依存し、厳格な基準によってその選択が立証されていないため、一意的に定義されていない。ＥＭＤに関する強固な数学的基盤なしでは、この純粋に経験的及び直観的なアプローチは、多くの難解な方法論的課題を生み出す。例えば、内因性モード関数は、解析的定義なしで、反復アルゴリズムの出力として単に数字的に定義されてしまう。これらの条件下で、方法論の異なる局面についての判断は、経験論に基づく状況に応じた比較によって、行われている。

生物医学的信号処理に対するモデル駆動型方法には、関与する系に対する正確なモデルが必要である。課題の具体的な局面は、生理信号が、複数の微視的スケールの情報源によって導かれた、大域的スケールの変数である点である。物理的モデルの作成には、関連する微視的スケールの成分の厳格な特定及び合成、これら根拠に基づいた大域的モデルの作成が必要となる。しかし、関与する細胞現象及び分子現象の独特なパターンを正確に記載することが不可能であるため、上記の厳格な決定論的概念を用いた大域的スケールから微視的スケールへの変換は、一意的な解を有さない。このような不確定性を考えた場合、大域的スケールの信号に着目した細胞機構のモデル作成に関する膨大な試みは、モデルに関する生理学的詳細及び解剖学的詳細だけでなく、基礎数学的手順という点で、大きく異なる顕著に多様な発見的手順によって、裏付けられている。本モデルの矛盾の程度は知られていない。

従って、従来の方法論では成功していない、非定常時系列信号の成分を効果的に特定する必要性は未だ存在しており、それゆえ、当該モデルは、信号の特性値を分類づけるのに効果的に利用することができるだろう。

〔概要〕
本発明の目的は、既存の手順の一つ又は複数の欠点を、実質的に克服、少なくとも、改善する点にある。

本発明の第一の態様において、基本成分として準ガウシアンカーネル（ＱＧＫ）を用いて自己相似機能要素のセットへと非定常信号を分解する方法であって、個々の機能要素を特徴づけるパラメータを推定するために、以下の工程を含む方法を提供する：
非定常信号の時間経過における動的トレンド及び基線のトレンドを推定し、ガイド関数及びその適合した分割を抽出する工程；
非定常信号を要素の発達（development）の領域に分ける工程；
先行する機能要素に適合するＱＧＫからの重複を再帰的に補正し、再帰サイクルに含まれる各サイクルにおいて経験的成分カーネル（ＥＣＫ）を形成する工程、そのサブ工程は以下を含む：
周波数領域へＥＣＫを変換する工程；
ＱＧＫのパラメータの推定及びそれらの検証を行う工程；及び、
非定常信号の大域的モデル及び部分的モデルを作成する工程。

本発明の他の態様においては、前記方法を実施するための装置を提供する。

本発明に係る他の態様も開示される。

〔図面の簡単な説明〕
本発明の一つ又は複数の実施形態は、以下の図面を参照して記載されるであろう：
図１Ａは、頭皮の表面に取り付けた二つのセンサーの間で測定される電位の典型的な波形を示す；
図１Ｂは、ＥＥＧ信号が示す極めて非定常的な特徴を例示する；
図２Ａは、聴覚オッドボールパラダイムにおいて検出される典型的なＥＲＰ信号を示す；
図２Ｂは、図２Ａに示した各セグメントを、そのセグメントに適合するパラメータを有する準ガウシアンカーネルを用いて置換した典型的な結果を示す；
図２Ｃは、準ガウシアンカーネル成分の和として得られたモデルを、比較のための図２Ａに示した元のＥＲＰ信号とともに示す；
図３は、生理信号を機能要素に適応的に分解し、その後、各機能要素を特徴づけるパラメータを決定する、波形分析方法の概略フロー図を示す；
図４Ａは、図３に記載の方法中で応用した、動的トレンド及び基線のトレンドの推定の概略フロー図を示す；
図４Ｂは、典型的な生理信号、マルチパスの移動窓による平均化を離散的な生理信号に適用して得られた基線のトレンド、及びトレンドを除去した基線の信号を示す；
図５Ａ及び図５Ｂは、ゼロ交差点から結果的に得られた分割点、並びに大域的及び局所的な最小値での分割点を各々、示す；
図５Ｃは、ゼロ交差点（実線の縦線）と最小値（点線の縦線）とにおいて、それぞれ、ＥＥＧ信号を分割した例を示す；
図６Ａは、高い時間分解能を持つ分割点を提供する初期信号及びガイド関数を示す；
図６Ｂ及び図６Ｃは、図６Ａに記載のガイド関数及び初期信号を、縦線で示した分割点とともに再度示す；
図７Ａは、図６Ａに記載の初期信号、及び低い時間分解能による信号分割を提供するガイド関数を示す；
図７Ｂ及び図７Ｃは、図７Ａに記載のガイド関数及び初期関数を、縦線で示す分割点とともに再度示す；
図８Ａは、準ガウシアンカーネルの構成に用いたガウシアン成分を示す；
図８Ｂは、準ガウシアンカーネルの一例を示す；
図９Ａは、有限成分の和への挿入箇所の分解を例示する；
図９Ｂは、三つの三角形基底関数の和へのハット関数の分解を示す；
図９Ｃは、同様の基底関数アルゴリズム中で使用した、三角形基底関数を示す；
図９Ｄは、三角形基底関数の複素スペクトルの実部と虚部とを示す；
図１０Ａは、パラメータσの推定に用いた振幅テンプレートを示す；
図１０Ｂは、パラメータσの推定過程における、正規化した実験振幅スペクトル、並びに異なる位置の振幅テンプレート、Ｔ１、Ｔ２及びＴ３を例示する；
図１１Ａは、階級Ｅに属する正規化振幅スペクトルに、テンプレートを重ね合わせた結果を示す；
図１１Ｂは、階級Ｇに属する正規化振幅スペクトルに、テンプレートを重ね合わせた結果を示す；
図１２は、図中央のパネルに記載の、ＥＥＧ信号の異なるセグメントの正規化振幅スペクトルに重ね合わせた振幅テンプレートのいくつかの結果を例示する；
図１３Ａは、低い時間分解能を用いたＥＥＧ分割の例を示す；
図１３Ｂは、図１３Ａに記載の１０セグメントの振幅スペクトルを示す；
図１３Ｃは、振幅と周波数との両方を正規化した後、図１３Ｂの振幅スペクトルを重ね合わせたものを示す；
図１４Ａは、パラメータτの推定値が得られることを基準にした場合の横座標周波数スケール上の点ａ及びｂを示す；
図１４Ｂは、横座標スケール上の点ｂが、複素スペクトルの虚部の最小値によって定義された場合を例示する；
図１５Ａは、話している人の呼吸の波形例を示す；
図１５Ｂは、図１５Ａの信号において特定した５４個のＱＧＫの波形を示す；
図１５Ｃは、特定した５４個のＱＧＫの和として合成されたモデルによって、図１５Ａの信号を再現した場合の、正確性を示す。

〔発明の詳細な説明〕
工程及び／又は特徴に係る添付の図面の一つ又は複数において、参照がなされ、同一の参照番号が含まれるが、逆の意図が示されていない限り、これらの工程及び／又は特徴は、本説明の目的のため、同一の機能又は作用を有している。

非定常信号を記載するため、重要な工程がモデルを構築している。この課題に対する従来アプローチは、「物理学に基づいた」モデルによって裏付けられている：その目的は、発生源の明白な物理モデルという点で、信号を記述することである（Ｃａｎｄｙ Ｊ．Ｖ． “Ｍｏｄｅｌ−ｂａｓｅｄ ｓｉｇｎａｌ ｐｒｏｃｅｓｓｉｎｇ”． Ｊｏｈｎ Ｗｉｌｅｙ ＆ Ｓｏｎｓ， Ｎｅｗ Ｊｅｒｓｅｙ， ２００６ 参照）。しかし、非定常信号は、通常、複数の非線形的発生源から発生し、これらの相関性は未だ明確でない。

これらの物理学に基づいたアプローチとは異なり、本明細書に開示の方法は、カオスに基づいたモデリング（ＣＢＭ）アプローチによって裏付けられる。これらは、発生源の振る舞いを定義する決定論的要素及び確率論的要素の両方を組み合わせている。このアプローチの一般的なフレームワークは、生物医学的信号の発生の基礎を成す一般的なメカニズムに関する文献によって紹介されている。

生物医学的信号解析に関する一般的な側面は、生理学的プロセスは、巨視的スケールで振る舞う複数の初期発生源によって生成される大域的スケール変数である、という点である。この観点において、生理学的プロセスは、分子的要素及び細胞的要素の非常に複雑な系の外部に対するパフォーマンス、つまり、決定論的要素及び確率論的要素で定義されるパフォーマンスを反映している。必要とされる機能要素数は非常に多く、細胞機構の詳細は非常に複雑であるため、該当する細胞全体のモデルは実現不可能である、という事実により、大きな困難が生じている。しかし、複数の微視的スケール変数によって生成される集団効果に関する、数多くの理論的研究及び計算機科学的研究が、大域的スケールにおいて重要である細胞機構に関する一般的側面に、そうでないものから、光を投じた。よって、局所的な微視的スケールのボリュームにおける、多くの初期発生源の同期的活性化は、大域的スケールにおいて測定される信号における検出可能な変化を生成するための必要条件である、というのが一般的な物理学的概念である。これらの処理に関する主要な側面は、脳波図（ＥＥＧ）、心電図（ＥＣＧ）及び筋電図（ＥＭＧ）のような、臨床的に重要な生物医学的信号に関して、幅広く研究がなされてきた。明確な生理学的要素の指標として、ピークを持つ大域的スケールの波形を示唆する一般的なフレームワークは、ＥＥＧの生成の基礎を成すメカニズムの参照によって示され得る。

図１Ａ中に示すように、ＥＥＧは頭皮の表面に取り付けたセンサーを用いて測定され得る。電位は、頭皮下の脳の皮質領域内のニューロンの活性によって引き起こされる。単一ニューロンの活性は、十分高くないため、頭皮の表面に取り付けたセンサー間で測定可能な電位を引き起こすことができない。むしろ、接近して存在する皮質ニューロンの大きな集団が同期して作用すると、活性の和によって、細胞外双極子が形成され、その結果、測定電位が得られる。

頭皮の表面で記録されるこれらの集合電位の変動及びダイナミクスは、タイムラグなしで皮質の神経処理装置の活性に従う。この観点において、ＥＥＧ信号は、ヒトの脳における情報処理の基本メカニズムをリアルタイムに研究するための、最も重要な客観的指標の一つとしての機能を果たしている。図１Ｂ中に示すＥＥＧ信号の一例は、脳の電気的活性は、典型的に極めて不規則かつ非定常的であり、様々な活性パターンを生成することを示す。

ＥＥＧ信号という観点で情報処理を理解する鍵は、ＥＥＧ信号中の変化及び認知的課題中の特定刺激との時間関係を検出することによって提供される。ＥＥＧ信号のこのような特定の断片は、事象関連電位（ＥＲＰｓ）と呼ばれている。図２Ａは、聴覚オッドボールパラダイムにおいて誘発される典型的なＥＲＰを示す。オッドボールパラダイムは、２つのタイプの音程、すなわち、バックグラウンド刺激及びターゲット刺激を用いてＥＥＧ信号の特徴に影響を与える試験である。図２Ａ中に示すＥＥＧ信号において、ｔ＝０において、ターゲット刺激を与える。この状況の後に、連続した正負の波形の揺れが生じる。重要な実験的及び臨床的証拠を元に、刺激を与えてから特定の時間間隔において観察される波形のピークは、様々な皮質ニューロンネットワークに関連する機能エンティティとして見なされている。聴覚オッドボールパラダイムにおいて、このようなピーク電位は、各々Ｎ１００、Ｐ２００、Ｎ２００及びＰ３００と示されている後期潜時(late latency)ＥＲＰとして知られている。数字は、刺激後の負のピーク（Ｎ）又は正のピーク（Ｐ）の出現のおおよその時間（単位はミリ秒）を表している。Ｐ３００の特定の側面は、反復ＥＲＰの平均を検討する際、後期潜時要素は、巨大な一つのピークとして現れる点である。しかし、単一のＥＲＰ波形、すなわち、元のＥＥＧ信号を検討する際、後期潜時ＥＲＰであるＰ３００は、通常、Ｐ３ａ及びＰ３ｂと示されている２つの補要素からなる複雑な電位を示し得る。

上記様式で名づけられる異なる正負の変動は、重要な生理学的システムに密接に関連する一連の活性化処理の客観的指標として現れる。特定の構造機能機構を有する機能系によって生成される要素の指標としてのピークを持つ波形の概念は、科学的及び臨床的に非常に重要なほぼすべての生理学的プロセスへと従来応用されている。広く受け入れられたアプローチを考えると、種々の異なる手動又はコンピュータを導入した測定技術（ピーク感知手段）は、生理信号をピークの振幅及びピークの潜時へと変換する。これら広く実験的な技術の致命的な限界は、離れた時間点におけるピークの振幅及び潜時への信号の変換によって、重要な機能的情報及び臨床的情報を含んだ生物学的波形のダイナミクスを特定できない点である。計算的観点において、これらすべての技術の欠点は、要素の時間的重複を分解できない点である。重複を修正しなければ、振幅及び潜時のパラメータ測定は、非常に不確定なエラーの影響を受け得る。

図３は、単なる波形のピークではなく、特定の波形を持つ（正又は負の）振れとして要素を捉える波形解析の方法１００の概念フロー図を示す。より具体的には、方法１００は、非定常信号ｖ（ｔ）を機能要素に適応的に分解し、その後、時間−周波数解析及びパラメータ推定手順を適用し、動的モデル及び関連パラメータという形で、各々の機能要素を特定する。非定常信号ｖ（ｔ）は、ＥＥＧ信号、心電図（ＥＣＧ）信号、筋電図（ＥＭＧ）信号、呼吸由来の信号等であり得る。

コンピュータシステム内で実行される、一つ又は複数のアプリケーションプログラムのようなソフトウェアを用いて、方法１００を実施し得る。代わりに、一つ又は複数の集積回路のような方法１００の機能又は副次機能を果たす、専用ハードウェアにおいて、方法１００を実行してもよい。このような専用ハードウェアは、デジタル信号処理装置、又は一つもしくは複数のマイクロプロセッサ及び関連メモリを含み得る。

方法１００は、工程１０１から始まり、アナログ信号である非定常信号ｖ（ｔ）を調整し、標本値ｖ_ｍ＝ｖ（ｔ_ｍ）からなる離散入力数列とする。より具体的には、非定常信号ｖ（ｔ）を標準電圧の範囲内で増幅し、一定のサンプリング間隔ｔ_ｍ＝ｍΔ（ｍは整数値であり、Δはサンプリング間隔である）において抽出する。この手順により、標本セット

が形成される。非定常信号ｖ（ｔ）は、帯域制限信号として見なすことができ、サンプリング間隔Δは、抽出された非定常信号ｖ（ｔ）の正確な再現に必要な最小のサンプリングレートである、ナイキストレートを超えるように選択する。

詳細は以下に記載するが、その後、方法１００は、非定常信号ｖ（ｔ）の離散的表現である、離散入力数列

上で作用し、入力数列

をセグメントに分割し、各々のセグメントに時間−周波数解析を適用し、セグメントに対応するパラメータと動的モデルとを、各々のセグメントにおいて特定する。パラメータ数列は以上のように形成され、パラメータ数列は、使用者が定義した基準に応じた各々のセグメントを分類し、非定常信号ｖ（ｔ）又は定義されたセグメントを合成するために使用される。

よって、非定常信号ｖ（ｔ）、及び離散対応値

は、典型的に高い周波数のノイズ及び低い周波数のバックグラウンド活性を含んでいる。これらの不適切な活性は、通常、ローパス及びハイパスフィルタを用いて、更なる解析から取り除かれる。しかし、重要な機能要素及び不適切な活性の両方の、非常に不規則及び非定常的な性質が原因で、時間不変的周波数の特徴を持つフィルタの適用は、機能要素のダイナミクスを、著しくかつ予期せず歪ませ、機能的に重要な情報が失われ得る。

この課題に対する一般的な解決策は、適応フィルタの使用である。このような適応フィルタは、信号のダイナミクスにおいてトレンドを変えるために動作を適応させなくてはならない。最も一般的にトレンド推定に用いられるデジタルフィルタは、移動平均フィルタ又はボックスカーフィルタであり、標本各々を、問題になっている標本の平均及び任意の数の連続する前後の点に、置換する有限インパルス応答を用いたデジタルフィルタである。しかし、ボックスカーフィルタの周波数応答は、最適ではなく、単純な移動平均は、トレンド軌道において著しい歪みの原因となる。このような歪みを防ぐため、工程１０２では、入力数列ｖから、高い周波数のノイズ及び低い周波数のバックグラウンド活性を省くため、マルチパス移動平均（ＭＰＭＡ）フィルタを用いる。より具体的には、２つのＭＰＭＡフィルタは、入力数列ｖにおける動的及び基線トレンド推定のため、入力数列ｖに適用される。

ＭＰＭＡは、反復処理であり、入力数列に移動窓平均を適用し、一回目の反復の入力数列は入力数列ｖであり、それに続く次の反復の入力数列は先行の反復の結果となる。従って、ｋ回目の反復において、各々、入力数列は

であり、出力数列は

であり、式中、標本数はＭ＋１、ｋは１からＫの値を取る反復数であり、

である。

ｍ個目の標本に対する出力数列

は、以下のように与えられる：

式中、Ｐは、窓の広さを定義し、当該窓の広さはＮ_Ｗ＝２Ｐ＋１である。

Ｐの連続したその前の点及びＰの連続したその後の点によって

が定義されることは、方程式（１）より明らかである。しかし、数列ｖの定義によると、問題となっている処理の第一の標本は、

である。結果として、ｍ＋ｊ＜０である値

は、計算から省かれる。これは、ｎが−１から−Ｐの値を取るｖ_ｎ＝０のゼロ値、及び拡張数列ｖからなるゼロパディングとしてみなされ得る。

従来の移動窓平均によって、結果的に得られた数列のスペクトル中の高周波数の歪みが生成される一方で、工程１０２において適用されるＭＰＭＡは、スペクトルから徐々にこのようなずれを除去する。ｍがＰからＭ−Ｐの値を取る残差

は、各反復に従って小さくなる。ｋ回反復後のトレンド推定における改善は、反復間の標本値変化を測定することで評価する。正確性は、以下の平均二乗誤差で表す：

反復は、所望の改善が達成されると終了する。ＥＥＧ、ＥＲＰ、ＥＭＧ、及びＥＣＧ等の、典型的な電気生理信号に適用する場合は、以下の条件を満たす場合、ｋ回反復後に、反復は終了する：

式中、Ｍ_ｖは出力信号ｖ^ｋの平均値であり、Ｔ_ｒは、無次元の閾値である。これら適用において、０．０５が妥当な閾値Ｔ_ｒである場合、反復は、典型的に５反復後に終了する。

図４Ａは、工程１０２中に適用する動的及び基線のトレンド推定の概要フロー図を示す。異なる２つのＭＰＭＡは、フィルタリングされた数列を形成するため、ＭＰＭＡ１ユニット及びＭＰＭＡ２ユニット中で入力数列ｖに適用される：動的トレンドは

であり、基線のトレンドは

である。ＭＰＭＡ１及びＭＰＭＡ２は、異なる窓パラメータＮ_Ｗ１及びＮ_Ｗ２を有し、

の条件を満たさなくてはならない。入力数列ｖ、並びに変換数列ｄ及びｂは、結合され、２つの出力数列ｇ及びｕが生成される。

数列

は、

と共に、基線の再トレンド化信号、すなわち、基線のトレンドｂを差し引いた後の初期数列ｖを意味する。

図４Ｂは、基線のトレンドが除去された生理学的信号ｕの計算を例示する。実線は、Δ＝０．００４ｓで抽出した典型的なＥＥＧ信号ｖの２０１点での一部分である。点線は、３１点の窓にＭＰＭＡ２を適用し、５回反復を用いた結果得られた基線ｂである。ＥＥＧ信号ｖから基線ｂを差し引いた後の、トレンドを除去した数列ｕは破線で示されている。

を満たす数列

は、ガイド関数と呼ばれている。ガイド関数は、信号の時間経過から、不適切な高周波数及び低周波数の要素を取り除いたものを用いて、得られる。これらの純化処理は、重要な機能要素の探索のガイドとしての数列ｇの役割を裏付ける。

図３を再度見ると、工程１０２中のトレンド除去に続いて、工程１０３が続き、ガイド関数ｇに基づいた分割が実施される。工程１０３中に実施される分割の目的は、機能要素に関連した信号のセグメントを定義する分割点を推定することである。分割点は、図５Ａに図示されたゼロ交差（ＺＣ）、並びに図５Ｂに図示されたガイド関数のモジュールの大域的な最小値及び局所的な最小値として定義される。図５Ｃは、縦の実線で示されたゼロ交差点を有するガイド関数ｇ及び縦の破線で示された最小値の例を示す。

より具体的には、点ｔ_ｍは、以下の条件では、ゼロ交差点として限定される：

点ｔ_ｍは、以下の条件では、

の大域的な最小値又は局所的な最小値である：

分割点の数列

はこのように形成され、セグメントの長さ

は、以下を満たす分割点

によって定義される：

一般的に、ガイド関数ｇは、トレンドを除去した信号ｕの近似として考えられ得る。近似の正確性は、ＭＰＭＡ１の異なるパラメータを用いて修正することができる。従って、分割手順の時間分解能は、ＭＰＭＡ１のパラメータの適切な選択によって規定される。

図６Ａ及び図７Ａは、ＭＰＭＡ１の異なるパラメータを用いて、トレンドを除去したあるＥＥＧ信号又はトレンドを除去した同一のＥＥＧ信号ｕ（灰色線）に適合した２つのガイド関数（破線）を示す。より具体的には、図６Ｂ及び図７Ｂ中に別々に示されたガイド関数ｇ１及びガイド関数ｇ２は、各々三点の窓、及び七点の窓を用いて計算されている。両ケースとも、２反復を用いて計算されている。図６Ａに記載の結果は、トレンドを除去した数列ｕ及びガイド関数ｇ１の間のずれが少ないことを示す。よって、ガイド関数ｇ１は、トレンドを除去した信号を詳細に再現し、結果として、トレンドを除去した数列ｕの実質的に厳密な再現として見なすことができる。縦線で示す分割点は１１のゼロ交差点及び８の最小値を含む。

逆に、図７Ａに記載の結果、及び図７Ｂとの図６Ｂの比較は、より広い窓を用いた同じＭＰＭＡ手順によって、分割点から最小値を取り除き、結果として、複雑なピークを持つ波形を巨大な一つの波形によって置換することを示す。いくつかの場合、処理の重要な点は、機能的に不適切な活性を取り除き、ＥＥＧ、ＥＭＧ、ＥＣＧ等の、大域的スケールの生理学的信号の時間領域において、有力な動的トレンドを強調する点である。しかし、不規則的な高周波数の信号の変動は、物理学的に重要な発生源から生成され得る。例えば、ＥＥＧの場合、３０Ｈｚを超える周波数を示す活性は、ガンマリズムとして知られている。よって、装置の特徴に依存して、高分解能の分割及び低分解能の分割の両方が、これらの異なる組み合わせと同様に、複素非定常信号の要素構成の検討のための強力な計算道具を提供する。

工程１０３中で適用される信号変換は、時間領域において実施されるが、周波数領域の概念とも関連し得る。図６Ａ及び図７Ａから明らかなように、ＭＰＭＡ１中で使用される狭い窓は、より高い周波数の帯域幅においてローパスフィルタリングを実施し、一方で、ＭＰＭＡ２中で使用される広い窓は、比較的低い周波数の帯域幅においてローパスフィルタリングを実施する。しかし、当初の信号は、非定常プロセスであるため、固定された「カットオフ」周波数とともに従来の周波数帯に用いることは不明瞭である。

再度、図３を参照すると、工程１０３中の分割に続いて、方法１００は、セグメント再帰として、工程１０４から１０８を実行し、当該再帰の連続したサイクルは、定義された信号セグメントの連続を扱う。各々のセグメントに対して、要素モデルが特定され、当該要素モデルのパラメータは、対応する単一信号断片に適合している。このセグメント中、信号に適合した要素モデル各々は、以下のセグメントへと拡張されるため、工程１０４から１０８の反復工程の目的は、時間的要素の重複を解消することでもある。要素の時間的重複の解消は、一般的な方法では裏付けられていない信号処理の多くの分野において緊急課題である。この課題への如何なるアプローチも、要素の発達において、予想できる構造及び予想パターンを記述するため、及び過去の軌道から要素の今後の値を予測するため、要素モデル作成に関与しなくてはならない点に、主要な難しさがある。

「物理学に基づいた」アプローチ（Ｃａｎｄｙ Ｊ．Ｖ． “Ｍｏｄｅｌ−ｂａｓｅｄ ｓｉｇｎａｌ ｐｒｏｃｅｓｓｉｎｇ”． Ｊｏｈｎ Ｗｉｌｅｙ ＆ Ｓｏｎｓ， Ｎｅｗ Ｊｅｒｓｅｙ， ２００６ 参照）の欠点は、発生源の明確な物理学モデルという点において、通常、信号に係る十分な情報を得ることができないという事実に生じている。物理学に基づいたアプローチに対して、本明細書に開示の方法は、決定論的カオスという形で、非定常信号の主要要素を定義するカオスに基づいたモデリング（ＣＢＭ）を用いた一意的な方法論を紹介している。このカオスに基づいたモデリングアプローチは、正規化効果の実験論的結果（Ｍｅｌｋｏｎｉａｎ Ｄ， Ｂｌｕｍｅｎｔｈａｌ Ｔ， Ｇｏｒｄｏｎ Ｅ． “Ｎｕｍｅｒｉｃａｌ Ｆｏｕｒｉｅｒ ｔｒａｎｓｆｏｒｍ ｓｐｅｃｔｒｏｓｃｏｐｙ ｏｆ ＥＭＧ ｈａｌｆ−ｗａｖｅｓ： ｆｒａｇｍｅｎｔａｒｙ−ｄｅｃｏｍｐｏｓｉｔｉｏｎ−ｂａｓｅｄ ａｐｐｒｏａｃｈ ｔｏ ｎｏｎｓｔａｔｉｏｎａｒｙ ｓｉｇｎａｌ ａｎａｌｙｓｉｓ”． Ｂｉｏｌ Ｃｙｂｅｒｎ．， １９９９， ｖｏｌ． ８１， ｐｐ． ４５７−４６７ 参照）によって、奨励されている。この研究結果は、決定論的かつ不規則な要素の干渉が原因となり、微視的スケールからの発生源の活性が、初期発生源の物理学的性質よりむしろ、大域的スケールにおいて確率論的不変性を反映した統計学的手順へと、変換されることを示す。よって、大域的スケールの生理学的要素は、複数の初期発生源の累積統計的集合として、発達する。

この発見的基準の普遍的かつ大域的スケールの要素モデルを数式化するため、本方法は、ランダム変数の集合の漸近挙動を定義する手段として、中心極限定理を引用する（Ｇｎｅｄｅｎｋｏ Ｂ， Ｋｏｌｍｏｇｏｒｏｖ Ａ． “Ｌｉｍｉｔ ｄｉｓｔｒｉｂｕｔｉｏｎｓ ｆｏｒ ｓｕｍｓ ｏｆ ｉｎｄｅｐｅｎｄｅｎｔ ｒａｎｄｏｍ ｖａｒｉａｂｌｅｓ”． Ａｄｄｉｓｏｎ−Ｗｅｓｌｅｙ Ｐｕｂｌｉｓｈｉｎｇ Ｃｏｍｐａｎｙ， Ｒｅａｄｉｎｇ， Ｍａｓｓ．， １９５４ 参照）。平均

及び分散

とともに、一定数のランダム変数

の和と仮定すると、指数フーリエ変換

の漸近挙動は、

のようなガウス関数に対して極めて一般的な条件下において、収束し（ｓｅｅ Ｐａｐｏｕｌｉｓ Ａ． “Ｔｈｅ Ｆｏｕｒｉｅｒ ｉｎｔｅｇｒａｌ ａｎｄ ｉｔｓ ａｐｐｌｉｃａｔｉｏｎｓ”． ＭｃＧｒａｗ−Ｈｉｌｌ Ｂｏｏｋ Ｃｏ， ＮｅｗＹｏｒｋ， １９６２ 参照）

式中、Ｆ（ｉω）は、指数フーリエ変換（複素変数ｉωの関数、

）、

及び

である。

よって、単一発生源の影響は、大域的スケールの変数として意味を持つ概略パラメータσ及びηへと変換される。

物理学用語では、

は、多数の微視的スケールの発生源から生成される大域的スケールの時間プロセスの振幅スペクトルとして、見なされ得る。

この理論的予測の実験的立証は、聴覚の事象関連電位の研究によって、提供されてきた（Ｍｅｌｋｏｎｉａｎ Ｄ， Ｇｏｒｄｏｎ Ｅ， Ｒｅｎｎｉｅ Ｃ， Ｂａｈｒａｍａｌｉ Ｈ． Ｄｙｎａｍｉｃ ｓｐｅｃｔｒａｌ ａｎａｌｙｓｉｓ ｏｆ ｅｖｅｎｔ−ｒｅｌａｔｅｄ ｐｏｔｅｎｔｉａｌｓ． Ｅｌｅｃｔｒｏｅｎｃｅｐｈａｌｏｇｒａｐｈｙ ａｎｄ Ｃｌｉｎｉｃａｌ Ｎｅｕｒｏｐｈｙｓｉｏｌｏｇｙ， １９９８， ｖｏｌ．１０８， ｐｐ． ２５１−２５９ 参照）。動的スペクトル解析は、ガウス関数が、聴覚の事象関連電位の主要な後期潜時要素の振幅スペクトルのプロファイルに関する普遍的な記述を提供することを示した。

方程式（７）の周波数領域の複素スペクトルの時間領域対応値は、無限時間スケールの間隔に関連した該当する文献において定義される。方程式（８）と同様に、時間領域の解決策は以下の形のガウス関数である傾向がある：

逆に、開示の方法の信号処理のルーティーンは、生理学的に実現可能なシステム中において、一過性と見なされ得る因果時間関数に対応する。

無限時間スケールから半無限時間スケールへと時間領域解析を移動させるため、開示の方法は、ｔ≧０において、２つのガウス関数の和として定義される、原因関数ｗ（ｔ）を導入し、

式中、

であり、Ｕ（ｔ）は単位ステップ関数であり、εは無次元定数である。

方程式（１１）及び（１２）と方程式（７）を比較すると、方程式（７）からのパラメータη（平均）は、定数εを用いて、パラメータσ（分散）の比として、方程式（１１）及び（１２）中に記載されることがわかる。定数εの値は、非定常信号のタイプに依存する。従って、開示の方法を用いた、ヒト及び動物のＥＥＧに関する解析及びコンピュータシミュレーションは、これらの応用例に対する普遍定数として、ε＝２であることを示唆する。この推定値は、定数εのデフォルト値だと考えられる。

パラメータσ＝１及び定数ε＝２を満たす関数ｇ_Ａ（ｔ）及び−ｇ_Ｂ（ｔ）を、図８Ａ中に各々、破線及び点線で示す。これら数式の両端は、

である。

図８Ａ中、実線で示す、関数ｇ（ｔ）は、無限範囲の奇関数である。関数ｇ（ｔ）のフーリエ変換は以下であり：

式中、

である。

ｇ（ｔ）は奇関数であるため、

及び

である。

方程式（１７）からの最後の三角積分は、以下の表記が用いられるｗ（ｔ）のフーリエ正弦変換である。

方程式（１０）によってｗ（ｔ）を置換した後、この積分は以下の形式のように、解析的に記載される（Ｅｒｄｅｌｙｉ， Ａ． （Ｅｄ）． “Ｔａｂｌｅｓ ｏｆ ｉｎｔｅｇｒａｌ ｔｒａｎｓｆｏｒｍｓ“， Ｖｏｌ． １． ＭｃＧｒａｗ−Ｈｉｌｌ， Ｎｅｗ Ｙｏｒｋ， １９５４， ｐ． ７８ 参照）：

この式は、複素スペクトル

の虚数部分として見なすことができ、式中、

は、振幅スペクトルであり、かつ、

は、位相スペクトルである。

従って、方程式（２１）及び（２２）の周波数関数の組み合わせは、ｔ≧０において、ｗ（ｔ）の時間領域信号を補足的に定義する周波数領域に対応するものと見なすことができる。この関数を図８Ｂに示す。ｗ（ｔ）が最大に達する近似値は、

である。最大値は、

である。

物理的には、ｗ（ｔ）は、複数の微視的スケールの発生源の集合の活性化によって、ｔ＝０において誘導される一過的なプロセスと見なすことができる。別々の発生源の振る舞いは、確率論的要素及び決定論的要素の両方に規定され、決定論的カオスとして特徴づけることができる。このような振る舞いの形態の主要な特徴は、決定論的カオスを生成している系の非線形性である。この観点において、方程式（１１）で記載されるプロセスは、以下の非線形微分方程式システムで規定される動的システムのインパルス応答と見なすことができる：

式中、

は、

の導関数を意味する。

合成非定常プロセスのモデルは、以下の形式で定義される：

式中、κ_ｉ及びτ_ｉは、加重及び開始時間を意味する自由なパラメータである。

準ガウシアンカーネル（ＱＧＫ）と呼ばれる、要素モデルは以下のように定義される：

上記式は、３つのパラメータのセット、すなわち、

によって特徴づけられている。

一般的な観点では、ｖ_Ｎ（ｔ）の推定は、パラメータ

がランダム変数であるという事実に従っている。非定常性の特定の側面は、新たな発生源が、基礎を成すシステムに結合した際、時刻τ_ｉにおいてｖ（ｔ）の変化を生み出すシステムの構造的要素であるということである。決定論的カオスの重要な側面は、方程式（１０）〜（１２）のガウス依存性(Gaussian dependences)は、決定論的プロセスであるが、これら時間発展は、基礎を成す初期ランダム事象の漸近挙動として発展する、という事実によって、明らかにされている。微視的スケールにおける該当する変換の確率論的側面は、マルコフ過程、特に不均一な出生死滅過程という形で計算されてきた（Ｍｅｌｋｏｎｉａｎ Ｄ， Ｂｌｕｍｅｎｔｈａｌ Ｔ． Ｄｙｎａｍｉｃｓ ｏｆ ｔｈｅ ｅｙｅｂｌｉｎｋ ＥＭＧ ａｔ ｇｌｏｂａｌ ａｎｄ ｍｉｃｒｏｓｃｏｐｉｃ ｓｃａｌｅｓ． Ｔｈｅ Ｏｐｅｎ Ｃｙｂｅｒｎｅｔｉｃｓ ａｎｄ Ｓｙｓｔｅｍｉｃｓ Ｊｏｕｒｎａｌ，２０１３， ｖｏｌ． ７， ｐｐ． １１−２２ 参照）。

方法１００は、工程１０４から１０８を実行し、方程式（２４）の各々ＱＧＫ

に対するパラメータ

を特定する。特定アルゴリズムは再帰として構成され、当該再帰のガイド原則は方程式（２４）中の「ｎ」と「ｎ−１」との間の以下のような関係から導かれる：

関数ｖ_ｎ（ｔ）は、概念上はモデルであるため、この関数ｖ_ｎ（ｔ）は、トレンドを除去した数列ｕ（ｔ）の近似、すなわち、ｔ＞０において

であると見なすことができる。ＱＧＫ

の近似への寄与は、ｔ＝τ_ｎにおける

の初期値から始まる。つまり、以下のようである。

τ_１＝０とすると、

ｎ＞１のとき、

である。

方程式（２８）は、前項の影響がない場合において、再帰の最初のサイクルに対応している。方程式（２９）は、関数

によって前項の影響を受ける再帰ｎ回目に対応している。この関数は、再帰の前のサイクルにおいて得られる、主要なパラメータ

（ｉ＝１，・・，ｎ−１）によって解析的に定義される。

の役割は、時間要素の修正である。この観点において、方程式（２９）の右項は、解析的ＱＧＫ γ_ｎ（ｔ）の重複を修正した準実験的対応値を表す。

ＱＧＫ γ_ｎ（ｔ）のパラメータは、複素非線形関係によって、時間領域において相互に関連がある。周波数領域の利点は、周波数領域関係が、連続オーダーで実施される別々の手順によって、これらの関係を分離し、パラメータκ_ｎ、σ_ｎ及びτ_ｎの連続を推定する方法を提供する点である。対応する周波数領域の操作の重要な局面は、κ_ｎ及びσ_ｎの推定値が、解析下で、時間関数の時間シフトに対して不変である振幅スペクトルのプロファイルによって、定義される、という点である。これは、パラメータτ_ｎの暫定的な推定値（κ_ｎ及びσ_ｎの推定段階中）として、分割点

の値を用いる方法を提供する。従って、方程式（２８）及び（２９）は、以下の形式中、時間から周波数への変換の目的のために、提示される：

標準的な手順によって、時間から周波数への領域変換を裏付けるため、方程式（３１）の波形要素を

だけ移動させ、単一の方程式として記載する：

間隔は、

であり、

を満たす。この方程式は、経験的成分カーネル（ＥＣＫ）と名付けられる。

工程１０５中、ＥＣＫ ｙ_ｎ（ｔ）は、有限余弦フーリエ変換及び有限正弦フーリエ変換の数値推定を用いて、時間領域から周波数領域へと変換される。より具体的には、Ｍｅｌｋｏｎｉａｎ （２０１０）中（Ｍｅｌｋｏｎｉａｎ Ｄ． “Ｓｉｍｉｌａｒ ｂａｓｉｓ ｆｕｎｃｔｉｏｎ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｎｕｍｅｒｉｃａｌ ｅｓｔｉｍａｔｉｏｎ ｏｆ Ｆｏｕｒｉｅｒ ｉｎｔｅｇｒａｌｓ”． Ｎｕｍｅｒ Ａｌｇｏｒ．， ２０１０， ｖｏｌ． ５４， ｐｐ． ７３−１００ 参照）に開示の、相似基底関数(similar basis function)（ＳＢＦ）アルゴリズムを用いており、その内容は参照によって組み込まれる。

ＳＢＦアルゴリズムは、多項式挿入を用いた三角積分の数値推定における、最大精度を提供する、ファイロン様式の原型である。高速フーリエ変換（ＦＦＴ）及びその改良によって実施される、従来のデジタルスペクトル解析に反して、ＳＢＦアルゴリズムは、時間領域及び周波数領域の両方において、計算点の厳密なグリッドを扱う必要がない。同じことは、工程１０５中、時間から周波数への変換に従う信号セグメントの長さが可変である本ケースにおいても、特に有利である。周波数領域において、計算点の柔軟な選択の主な利点は、対数周波数スケールを用いた周波数領域依存性を表すことができる点である。

ＳＢＦアルゴリズムは、普遍的な手順で実現されているため、添え字「ｎ」（再帰における回数）は、アルゴリズムの以下の記載において省略する。

ＳＢＦアルゴリズムは、以下の有限余弦及び正弦フーリエ変換の推定に対応している：

式中、ω＝２πｆは角速度、及びｆは周波数である。これら方程式中に使用される表記

は、ｔ≧λにおいて、ｙ（ｔ）＝０と仮定する。これらは、ｙ（ｔ）の複素スペクトルの実数部分及び複素部分であり、

を満たし、式中、

である。

対応する振幅スペクトルは以下である：

デジタル形式へのｙ（ｔ）の変換は、

を満たす有限数の節点

中の、標本値

によってｙ（ｔ）を特定する。節点は均等に分けられる必要はない。よって、計算は、データ点のセット

に対応している。これらデータ点を用いて、ｙ（ｔ）は、区分線形多項式ｈ（ｔ）によって標本に挿入される。挿入は、以下の挿入条件を満たさなければいけない：

式中、

である。

の条件下で、ＳＢＦアルゴリズムは、相似基底関数：

（式中、ａ_ｉは加重係数であり、

は、相似基底関数（ＳＢＦ）である）の和として、間隔［０，λ］間に形成される。ＳＢＦは、相似関係によって定義される：

この単純にパラメータ化された時間スケーリングは、「三角形基底関数」（ＴＢＦ）と呼ばれる基本的な有限要素、

より、ＳＢＦのファミリーを生成する。ＴＢＦは、図９Ｃ中に図示されている単位直角三角形である。

有限要素を用いて、挿入式ｈ（ｔ）は、以下の形式で表される：

式中、

は、ハット関数である。ハット関数

は、以下のように定義される：

上記式は、メッシュ

における式である。図９Ａは、方程式（４１）中に記載の様式で、ハット関数

の連続（標本１、２、３、４）に続く、「０」標本へと取り付けられたＴＢＦの和として、挿入式ｈ（ｔ）の分解を例示する。

方程式（３７）及び（４１）を比較すると、方程式（４１）のハット関数は、ＴＢＦによって方程式（３７）中で置換されていることがわかる。この結合において、有用だとわかる幾何学法則は、ハット関数（三角形ａｂｃ）を３つのＴＢＦ（三角形ｏｈｃ、三角形ｆｏｄ及び三角形ｇｏａ）の和に分解することを示す、図９Ｂによって、図示されている。方程式（３７）中のすべてのハット関数にこの分解を適用すると、挿入係数の推定のため、以下の方程式となる：

式中、

である。

数値解という観点において、方程式（４２）の構築の主な利点は、ＴＢＦからの余弦フーリエ積分及び正弦フーリエ積分が、以下のような解析的形式で表される点である：

図９Ｄは、Ｒｅ（方程式（４５））及びＩｍ（方程式（４６））と図示した、これら方程式を示す。

フーリエ変換の理論（相似定理）によると、時間領域における横座標の圧縮は、横座標の拡張に加えて周波数領域における縦座標の縮小に対応する。これらの処理は、ＳＢＦ及びＴＢＦのフーリエ積分間の以下の相似関係を明らかにする：

これらの方程式を用いて、方程式（３８）の補間式(interpolant)の余弦積分及び正弦積分は、以下の解析方程式で表される：

ｈ（ｔ）はｙ（ｔ）への補間式であるため、これらのフーリエ変換は、方程式（３３）及び（３４）の有限フーリエ積分の推定値として見なされる。結果として、

この観点において、方程式（４９）及び（５０）中の関数Ｒｃ及びＲｓの、その解析式である（４５）及び（４６）による置換は、以下の計算方程式を導く：

方程式（４３）からの荷重係数推定後、方程式（５２）及び（５３）は、データ点へと直接適用される。

ＳＢＦアルゴリズムの利点は、角周波数の任意のセットにおいて、周波数特性を計算できる点である。以下のパラメータ特定手順は、対数周波数スケールを用いる。従って、周波数特性の標本

は、複数点

において推定され、
式中、１を超えるＣは、対数横座標スケールにおいてサンプリングレートを定義する定数であり、ω_Ａは、初期周波数である。

１０に対して、１０、２０、５０及び１００標本に対応するＣの典型値は、各々１．２５８９、１．１２２、１．０４７１、及び１．０２３３である。

振幅スペクトルは、同じ点において以下のように計算される：

ｎ番目の再帰サイクルと仮定すると、対応する振幅スペクトル

は、κ_ｎ（加重）パラメータ及びσ_ｎ（分散）パラメータの特定の役割を果たす。この目的のため、方法１００（図３）は、次に工程１０６に進む。

角周波数ω_Ａが、

のように選択される場合、加重は、

のように定義される。

この観点で決定的に重要なものは、実験的振幅スペクトルのモデルとしての方程式（２１）のガウシアン振幅スペクトルの妥当性試験である。当該試験は、方程式（５５）から計算された正規化振幅スペクトルに適用され、方程式（５６）からパラメータκ_ｎの推定値を用いる：

方程式（２１）のモデルと、このスペクトルとの比較は、テンプレートとして、ガウシアンスペクトルを用いて、対数周波数スケールにおいて、実施される。後者は、以下のように解析的に定義され、

式中、νは、実数の無次元変数である。この関数を図１０Ａに示す。数値計算のため、テンプレートは、複数点の標本値によって定義され、

式中、１を超えるＣは、方程式（５４）中と同じ定数であり、Ｍ＋１は、テンプレートの標本数である。

よって、テンプレートは、

である、セグメント

に拡張する。境界点

は、矢印Ａ及びＢによって図１０Ａに図示されている。

テンプレートと比較される実験的な正規化振幅スペクトルの標本は、対数角周波数スケール上の以下の点に対応している：

周波数範囲（

）及びサンプリングレートは、Ｉパラメータ及びＣパラメータによって定義される。これらパラメータの適切な選択は、デジタル形式において、振幅スペクトルの十分正確な表示方法を提供する。振幅スペクトルの周波数の範囲は、テンプレートが定義した周波数の範囲を超えなくてはいけない。従って、Ｉ＞Ｍである。

テンプレート適合は、横座標スケールを通して、ガウシアン振幅テンプレートを動かすことによって、実行される。これは、テンプレート標本が、変化しない間、横座標上の位置が、定数Ｃによる、周波数の乗法によって、連動して移動する段階的処理である。よって、第一の移動は、点

の数列への方程式（５７）の複数点の数列の変化によって実行される。ｋ番目の移動の後の周波数点の数列は、

である。図１０Ｂは、実験的な正規化振幅スペクトル

（黒い実線）、並びに点Ｔ１、Ｔ２及びＴ３の異なる移動におけるガウシアン振幅テンプレートのセットを例示する。

テンプレートスペクトル及び実験的な振幅スペクトルの間のずれは、平均二乗誤差（ＭＳＥ）を用いて、各々のテンプレートの位置において評価され：

式中、ｋは位置番号である（ｋ＝１、・・・、Ｉ−Ｍ＋１）。

角周波数の点では、テンプレートの初期位置は、実験的スペクトルの周波数範囲

を含んでおり、式中、

である。ｋ回の移動後、テンプレートに含まれる周波数範囲は、

であり、式中

である。

最良のテンプレート適合は、ＭＳＥの最小値によって示される。対応する移動回数をＫとすると、パラメータ

は、対数角周波数スケールにおいて、ＭＳＥが最小であるテンプレートの位置を定義する。この適合が、方程式（５５）の実験的な振幅スペクトルの適切なモデルとして、方程式（２１）の振幅スペクトルを実証するかどうかを決定するため、これら関数間の残差を、標本によって定義される無次元誤差関数

として見なし、

式中、ｊは、０からＭの値を取り、

である。

ここでは、

の両方は、等しい初期値

を持つ無次元関数である。従って、

と仮定する。周波数の増加に伴う誤差関数の変化は、周波数の範囲

を推定するため、閾値ε_Ｔを用いて評価され、上記範囲においては、

である。

境界角周波数ω_Ｊにおける

の最終値は、

として示される。この推定値は、以下の図式に従って１から５の整数値を取るカテゴリー化パラメータξを用いて、テンプレート適合の質を確認するために用いられる：

階級、すなわち、ξの値によって、モデリングの結果は以下のように評価される：Ｅ−優れている(excellent)、Ｇ−よい(good)、Ｓ−満足できる(satisfactory)、Ｗ−不十分(weak)、Ｕ−不成功(unsuccessful)。

指令に依存して、階級Ｗ又は階級Ｕの場合、加重は無視される、すなわち、κ_ｎ＝０である。

図１１Ａ及び図１１Ｂは、閾値

を用いてテンプレート適合の２つの実施例の誤差関数を図示する。対応する境界周波数は、図１１Ａ及び図１１Ｂ中、縦の点線によって図示されている。これらの線は、図１１Ａ及び図１１Ｂ中の振幅スペクトルと、各々、０．０６５及び０．３０１の値において交差する。その結果、対応するモデルは、各々、階級Ｅ及びＧに属する。

最良な適合は、テンプレート標本が、物理エンティティ、つまり、対数横座標スケール上の

の対応サンプルにおける角周波数の値、と関連し得るテンプレート標本における、テンプレートのいくつかの固定位置を定義する。この観点において、テンプレート標本は、以下の形式の角周波数の関数として、表され、

式中、下付き文字「Ｆ」は、テンプレートの最終位置を示し、σ_ｎは、ｎ番目のサイクルの再帰における方程式（２１）より得られたσパラメータを示す。

任意のω_ｍにおいて、所要のパラメータσ_ｎは、以下のように計算される

計算の安定した正確性のため、角周波数ω_ｍは、

が０．４から０．８の値を取る周波数の範囲から選択される。

図１２は、テンプレート適合の典型的な結果を例示する。中央のパネルは、縦の点線で示された分割点を有するＥＥＧの一部を示す。分割点間で抽出された６個の連続したＥＥＧ断片の振幅スペクトルは、上側のパネル及び下側のパネルにおいて、実線で表される。対応するテンプレートは、灰色の線で記載される。テンプレート適合の例は各々、以下の情報によって、補足される：モデルの階級（Ｅ又はＢ）、角周波数

及び値

である。

図１３Ａ中に同じ信号を再描写するため、分割点の間の信号断片の振幅スペクトルが、図１３Ｂ中に示される。図１３Ｃは、振幅と周波数との両方を正規化した後のこれらスペクトル特性の重ね合わせを示す。これらの点は、ガウシアン振幅テンプレート（黒い線）との適合が、低い周波数の範囲又は中間の周波数の範囲において優れているという点において、典型的である。コンピュータシミュレーションは、高い周波数の偏差は、機能要素と混同した不適切な活性によって生成されることを示してきた。

方法１００（図３）に再び目を向けると、工程１０６に続いて、方法１００は、開始時間τ_ｎが推定される工程１０７へと続く。この手順は、方程式（５５）を用いた振幅スペクトルの正規化と同様の方法において正規化した後、正規化振幅スペクトル

及び複素スペクトルの対応する虚部を用いている：

式中、下付き文字ｎは、再帰サイクル数である。

図１４Ａ及び図１４Ｂは、ＡＳ（正規化振幅スペクトル

）及びＩｍ（複素スペクトル

の正規化虚部）と示すこれら関数の典型的な形を例示する。これら関数の一般的関係は、

の形式中、位相スペクトル

を含む。

方程式（２２）の位相スペクトルと同様、

は、以下のように、角周波数への線形従属として表される。

方程式（２２）及び方程式（６９）の位相スペクトル間の差は、シフト係数

によって、説明される。

物理的に、この移動は、要素開始時間τ_ｎ及び分割点

の値の間の差を示す。この観点で、

である。

先行の再帰サイクルにおける、ＳＢＦアルゴリズムの手順は、解析下において、ＥＣＫの未知の開始時間の置換として、分割点を取る。従って、分割点は、数値積分が開始する時刻としての機能を果たす。κ_ｎ及びσ_ｎパラメータの元である振幅スペクトルは、対応する時間関数の時間移動に対して不変であるため、この条件は、当初定義したパラメータ値において、ほとんど効果を有さない。

τ_ｎ推定は、

に従った方程式（７１）を参照する。

この方程式を用いるため、未知のパラメータβ_ｎを推定すべきである。方程式（６８）及び方程式（６９）を参照すると、パラメータβ_ｎは、関数

の間の以下の関係によって、完全に定義される：

この方程式からパラメータβ_ｎを見出すため、角周波数の単一点に対して

の推定値を有することは、十分である。しかし、計算誤差を減らすため、２つの推定値は、異なる角周波数

において、計算される。

と仮定すると、これら周波数の選択は、ω_１ｎからω_２ｎの間隔において、

は、ωの増加と共に減少すると推定される。

図１４Ａ及び図１４Ｂ中、横座標スケール上、文字ａで示される縦の点線は、

を満たす点として、角周波数ω_１ｎを定義する。方程式（７３）への参照は、この点において

であることを示す。方程式（７２）へのこの値の挿入は、開始時間の第一の推定値を提供する：

方程式（２１）及び方程式（２２）の形式中、周波数領域モデルに従うと、

の期待される形式は、図１４Ａ中に図示されるように、極小値に続くゼロ交差によって特徴づけられる。このような場合を考えると、角周波数ω_２ｎは、

が横座標スケールと交差する、すなわち、

を満たす点として定義される。

図１４Ａ中、ω_２ｎは、横座標スケール上のｂから始まる縦の点線によって示される。位相スペクトル

の値は、

を定義する。

方程式（７２）は、以下のように、位相移動の第二の推定値を提供する。

いくつかの場合、ＥＣＫの不適切な要素は、

の期待される周波数領域のプロファイルを歪める。図１４Ｂ中のプロットＩｍは、横座標スケールを交差する代わりに、最小に到達する

の特徴的な例である。対応する角周波数は、横座標スケール上のｂから始まる縦の点線によって示される。開始時間の第二の推定値を得るため、この周波数は、方程式（７５）中、ω_２ｎを置換する。

τ_ｎ［１］及びτ_ｎ［２］を仮定すると、開始時間の最終的な選択は、σ＝σ_ｎを満たす方程式（１９）から

の間のずれの、異なるτ_ｎの値に対する評価が元になっている。当該ずれは、以下の平均二乗誤差によって計算され、

式中、標本

は、τ_ｎの推定値であるτ_ｎ［１］又はτ_ｎ［２］を用いて計算され、Ｎは、

のように、選択される整数である。

ＭＳＥ［１］及びＭＳＥ［２］が、τ_ｎ［１］及びτ_ｎ［２］と関連があると仮定すると、ＭＳＥ［１］≦ＭＳＥ［２］の場合、τ_ｎ＝τ_ｎ［１］である。それ以外は、τ_ｎ＝τ_ｎ［２］である。

よって、工程１０６及び工程１０７は、検討中のセグメントに対するパラメータ

を推定する。工程１０８は、方法１００がその後、工程１０４から工程１０８が、次の分割のパラメータ

を推定する工程１０４へと戻る前に、方程式（２６）を用いて、以下の部分モデルを推定する。

図２Ｂは、対応するＥＣＫに適合するパラメータ

を有するＱＧＫと共に、図２Ａ（Ｎ１００、ｐ２００、Ｎ２００、Ｐ３ａ及びＰ３ｂ）中に示す各々機能要素の置換の典型的な結果を図示する。セグメントは、図２Ａ中の実験的な対応結果として、同じ記号Ｎ１００、Ｐ２００、Ｎ２００、Ｐ３ａ及びＰ３ｂによって、明示されている。図２Ｃは、比較として図２Ａ中に示す当初のＥＲＰ信号と共に、要素ＱＧＫの和、すなわち、大域的モデルを図示する。信号分解及び特定要素のパラメータは、発明の方法の一般的な手順から完全に導かれたものであり、以降で適用される現象に適合させるための如何なる修正も行っていないことから、これは、満足できる程度の適合である。

この結論の指示における更なる証拠は、異なる発生源によって生成される非定常信号を用いて、得られている。この点を例示するため、図１５Ａは、Δ＝０．０５６ｓにおいて抽出された、話している人の呼吸の波形の６０秒セグメントを示す。プロットは、３点又は４７点の窓フィルタに対して５つのパスを用いて、呼吸センサーからトレンドを除去した信号を表す。ｙ軸ユニットは、生の呼吸センサーの離散的な出力である。

図１５Ｂは、特定した５４個のＱＧＫの重ね合わせを描く。すべてのモデル要素の和として合成された呼吸波形の大域的モデルは、図１５Ｃ中、実線で示される。同じ図において、点線は、図１５Ａから再現された当初の呼吸波形である。

ＥＥＧと同様に、呼吸の波形は、高度な非定常信号である。しかし、電気の発生源（皮質ニューロン）によって生成されるＥＥＧと異なり、呼吸器の動きは、胸壁の動きを規定する数多くの生体機械的要素に依存する。よって、呼吸器の動きの基礎となる物理学的プロセスは、脳波の電気力学とは、著しく異なっている。それにもかかわらず、本明細書に開示の方法は、両方の適用に対し、等しく効果的な要素の特定ツールを提供する。実験的波形への顕著に正確なＱＧＫの適合は、異なる生物種及び異なる実験条件下において記録される、様々な電気生理学的信号（ＥＥＧ、ＥＭＧ、ＥＣＧ）及び呼吸の波形にとって典型的である。これらの結果は、発明の方法論のカオスに基づくモデリング（ＣＢＭ）のアプローチを強く裏付ける。

すべてのセグメントが処理され、パラメータξ、κ、σ、τが工程１０３において形成された各々のセグメントに対して推定された後、方法１００は、モデリングの結果が、指標としてパラメータξを用いて検討される、工程１０９へと続く。各々のセグメントに対するパラメータξの値は、定義された閾値ξ_Ｔと比較される。ξ_ｎ≧ξ_Ｔの場合、対応する要素のモデルは、更なる処理段階から除かれる。ξ_Ｔの典型値は、４及び５である。

Ｎ_Ｉの不適当な要素を除いた後、残っている要素は、Ｎ_Ｓ−１が分割点の数である、ｎ＝１からＭ＝Ｎ_Ｓ−Ｎ_Ｉの連続的かつ順番に、番号が付け替えられる。対応するパラメータのセットは、以下のように記述され、

式中、

は、対応するパラメータのＭ値を各々含むベクトルである。

ベクトルτの要素が、定義されたモデル要素の発達が開始する時間点の有限集合を定義する一方で、ベクトルκ及びベクトルσの要素は、モデル要素の波形を定義する。この観点において、ベクトル

は、データ点のセット

を提供し、これは時系列として、パラメータκの値の数列を表す。さらに、ベクトル

は、パラメータσの値の時系列を提供する。

時系列は、信号処理の多くの方法が適用される数値データの標準形式である。よって、信号を標準パラメータのセットに分解する開示の方法１００は、データ正規化、データ変換、主成分分析、自己相関、自己回帰モデル、データクラスタリング、予測等、といった、非定常信号処理の、従来の時系列解析の方法へと適用する方法を提供する。

要素モデル（信号解析）のセットへの非定常信号の分解の結果を用いて、方法１００は、異なる問題を解決するための方法を提供し、これが信号合成である。これは、あらかじめ定義された階級に属する要素の複合体の形成に関与する。

合成のための、要素の収集は、工程１１０によって実行されるパラメータ階級分け手順から開始する。開始時間は、変化しないままで、階級分けは、各々特定された要素モデルに対して、パラメータξ、κ及びσへと適用される。より具体的には、パラメータの部分集合は、以下の条件に従って、モデルパラメータＰ_Ｍの全体集合から定義され：

式中、

は、対応するパラメータの境界（最初と最後との）値を定義する定数である。

方程式（７８）より、パラメータの部分集合は、以下のように示され、

式中、オーバーハット（＾）は、選択された要素を意味し、Ｌは、要素番号である。明らかに、Ｌ≦Ｍを満たす。

最後に、方法１００は、工程１１０からのパラメータの部分集合

が、以下の形式において非定常信号の部分モデルの合成のため使用される、工程１１０において終了する。

階級分け及び合成は、複数の目的を有する手順であり、異なる目的を追求できる。例えば、極端に小さい振幅、又は極端に大きい振幅を持つ要素の除外は、信号の不適当な要素を取り除くことができる。κ及びσの特定値の選択によって、振動要素等を強調できる。

モデルパラメータ

のフルセット及び／又は、選択されたパラメータのセット

は、更なる処理の前に、コンピュータシステムの記憶媒体に蓄積されるか、外部コンピュータシステムへと伝達され得る。更なる処理は、最適デコンボリューション、発見的クラスタリングアルゴリズム、再帰的なベイズ推定等、といった異なったモデルに基づいた信号処理方法を含み得る。

前記の事項は、本発明のいくつかの実施形態のみを記述しており、発明の目的及び方針から外れない限りは、改良及び／又は変更を行ってもよく、実施形態は、例示されるにすぎず、これに制限されるものではない。

図１Ａは、頭皮の表面に取り付けた二つのセンサーの間で測定される電位の典型的な波形を示す。 図１Ｂは、ＥＥＧ信号が示す極めて非定常的な特徴を例示する。 図２Ａは、聴覚オッドボールパラダイムにおいて検出される典型的なＥＲＰ信号を示す。 図２Ｂは、図２Ａに示した各セグメントを、そのセグメントに適合するパラメータを有する準ガウシアンカーネルを用いて置換した典型的な結果を示す。 図２Ｃは、準ガウシアンカーネル成分の和として得られたモデルを、比較のための図２Ａに示した元のＥＲＰ信号とともに示す。 図３は、生理信号を機能要素に適応的に分解し、その後、各機能要素を特徴づけるパラメータを決定する、波形分析方法の概略フロー図を示す。 図４Ａは、図３に記載の方法中で応用した、動的トレンド及び基線のトレンドの推定の概略フロー図を示す。 図４Ｂは、典型的な生理信号、マルチパスの移動窓による平均化を離散的な生理信号に適用して得られた基線のトレンド、及びトレンドを除去した基線の信号を示す。 図５Ａ及び図５Ｂは、ゼロ交差点から結果的に得られた分割点、並びに大域的及び局所的な最小値での分割点を各々、示す。 図５Ａ及び図５Ｂは、ゼロ交差点から結果的に得られた分割点、並びに大域的及び局所的な最小値での分割点を各々、示す。 図５Ｃは、ゼロ交差点（実線の縦線）と最小値（点線の縦線）とにおいて、それぞれ、ＥＥＧ信号を分割した例を示す。 図６Ａは、高い時間分解能を持つ分割点を提供する初期信号及びガイド関数を示す；図６Ｂ及び図６Ｃは、図６Ａに記載のガイド関数及び初期信号を、縦線で示した分割点とともに再度示す。 図７Ａは、図６Ａに記載の初期信号、及び低い時間分解能による信号分割を提供するガイド関数を示す；図７Ｂ及び図７Ｃは、図７Ａに記載のガイド関数及び初期関数を、縦線で示す分割点とともに再度示す。 図８Ａは、準ガウシアンカーネルの構成に用いたガウシアン成分を示す。 図８Ｂは、準ガウシアンカーネルの一例を示す。 図９Ａは、有限成分の和への挿入箇所の分解を例示する；図９Ｂは、三つの三角形基底関数の和へのハット関数の分解を示す；図９Ｃは、同様の基底関数アルゴリズム中で使用した、三角形基底関数を示す；図９Ｄは、三角形基底関数の複素スペクトルの実部と虚部とを示す。 図１０Ａは、パラメータσの推定に用いた振幅テンプレートを示す。 図１０Ｂは、パラメータσの推定過程における、正規化した実験振幅スペクトル、並びに異なる位置の振幅テンプレート、Ｔ１、Ｔ２及びＴ３を例示する。 図１１Ａは、階級Ｅに属する正規化振幅スペクトルに、テンプレートを重ね合わせた結果を示す。 図１１Ｂは、階級Ｇに属する正規化振幅スペクトルに、テンプレートを重ね合わせた結果を示す。 図１２は、図中央のパネルに記載の、ＥＥＧ信号の異なるセグメントの正規化振幅スペクトルに重ね合わせた振幅テンプレートのいくつかの結果を例示する。 図１３Ａは、低い時間分解能を用いたＥＥＧ分割の例を示す。 図１３Ｂは、図１３Ａに記載の１０セグメントの振幅スペクトルを示す。 図１３Ｃは、振幅と周波数との両方を正規化した後、図１３Ｂの振幅スペクトルを重ね合わせたものを示す。 図１４Ａは、パラメータτの推定値が得られることを基準にした場合の横座標周波数スケール上の点ａ及びｂを示す。 図１４Ｂは、横座標スケール上の点ｂが、複素スペクトルの虚部の最小値によって定義された場合を例示する。 図１５Ａは、話している人の呼吸の波形例を示す。 図１５Ｂは、図１５Ａの信号において特定した５４個のＱＧＫの波形を示す。 図１５Ｃは、特定した５４個のＱＧＫの和として合成されたモデルによって、図１５Ａの信号を再現した場合の、正確性を示す。

Claims (7)

1. コンピュータによって実行される、生理学的現象を示す非定常信号を、基本成分としてＱＧＫ（準ガウシアンカーネル）を用いて自己相似機能要素のセットへと分解する方法であって、各々の機能要素を特徴づけるパラメータを推定するために、以下の工程を含む方法：
   信号の時間経過における動的トレンド及び基線トレンドを推定するため、マルチパス移動平均フィルタを適用し、トレンドを用いて、機能要素が発達する信号断片を示すガイド関数、並びに分割点として、ゼロ交差点、並びにガイド関数のモジュールの大域的及び局所的な最小値を用いるガイド関数の適合した分割を抽出する工程；
   上記非定常信号を要素の発達の領域に分けて、上記非定常信号を、上記ガイド関数の分割点に対応する点において分割する工程；
   先行する機能要素に適合するＱＧＫからの重複を補正し、ＥＣＫ（経験的成分カーネル）を形成する工程であって、再帰的なサブ工程は以下を含む工程：
   相似基底関数（ＳＢＦ）アルゴリズムを用いて、周波数領域へＥＣＫを変換する工程；
   ガウシアンプロファイルを用いて、ＥＣＫの正規化振幅スペクトルと周波数領域のテンプレートとの間の最善の適合を見出し、評価することによって、κ（加重）及びσ（分散）パラメータを推定し、かつ、Ｅ、Ｇ、Ｓ、Ｗ及びＵの分類階級を用いたパラメータ信頼性を検証する工程；
   位相スペクトルの定義された値に対応する特性点を見出すことによって、ＥＣＫの複素スペクトルの虚部から、特定のκ及びσパラメータと共に、ＱＧＫの開始時間τを推定する工程；
   新たに特定したＱＧＫを加えることによって、再帰サイクル中の先行するサイクルにおいて特定されたＱＧＫの和から構成された、非定常信号モデルを拡張する工程；及び
   特定結果の評価のため、分類パラメータξを用いて、不都合なＱＧＫを除去した後、特定された要素を再構成する工程；
   ＱＧＫの集合の選択によって、非定常信号の部分モデルを作成する工程であって、上記ＱＧＫの集合のパラメータは予め定義した階級に属する工程。
2. 以下に従って分類パラメータξは、１から５の整数値を取る、請求項１に記載の方法：
   

   

   式中、Ｅ−優れている、Ｇ−よい、Ｓ−満足できる、Ｗ−不十分、及びＵ−不成功である。
3. 非定常信号にローパス無限インパルス応答フィルタを適用することによって、動的トレンドを推定する、請求項１に記載の方法。
4. 上記ガイド関数の適応的な分割は、分割点として、ゼロ交差点のみを用いる、請求項１又は３に記載の方法。
5. 各々のＥＣＫに適合する、各々のＱＧＫの初期値は、当該各々のＥＣＫの初期値における分割点と一致している、請求項１に記載の方法。
6. 推定工程は、ＥＣＫの複素スペクトルの位相及びガウシアン複素スペクトルの位相の間において、最善の適合を見出すことによって、各ＥＣＫに適合する各ＱＧＫの修正された初期値を推定する、サブ工程をさらに含む、請求項５に記載の方法。
7. 非定常信号を機能要素に分解するための装置であって、各々の機能要素を特徴づけるパラメータを推定するために以下の手段を含む装置：
   信号の時間経過における動的トレンド及び基線のトレンドを推定し、ガイド関数を抽出及び分割する手段；
   非定常信号を機能要素が発達するセグメントに分ける手段；
   ＥＣＫを形成するため、前述の機能要素に適合するＱＧＫからの重複を補正する手段であって、以下を含む手段：
   周波数領域へ上記ＥＣＫを変換する手段；
   上記ＱＧＫの加重及び分散のパラメータを推定し、パラメータ信頼性を検証する手段；
   上記ＱＧＫの開始時間を推定する手段；
   再帰サイクルに含まれる各サイクル後に、非定常信号のモデルを拡張する手段；
   好ましくないＱＧＫの除去及び残りのＱＧＫの再構成を行う手段；及び
   予め定義した階級に属するパラメータ集合によって非定常信号の部分モデルを作成する手段。

Images