呼吸インピーダンスを推定するための方法及び装置
呼吸器系の音響インピーダンスは、被験者の気道に生成される振動から推定できる。インピーダンスは、結果として生じる流量の振動と圧力の振動との間の周波数依存的な関係を記述する。インピーダンスが吸気から呼気へ変化するとき、インピーダンスは高い時間解像度で推定されなければならない。生理的目的のために充分に短い時間間隔でインピーダンスを確実に推定する方法が提供される。不確定性原理の単純なバージョンが、離散的時間及び周波数に対して導出された。この原理に従って、最適時間―周波数解像度を与える離散的時間―周波数変換が開発された。変換は正規直交であり、離散的時間―周波数領域の分散の分析を可能にする。ノイズが流量及び圧力両方に存在するという仮定の下、インピーダンスは、時間―周波数領域の二変量の最小二乗分析に従う。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、強制的圧力振動から呼吸インピーダンスを推定するための方法と装置とに関する。
【背景技術】
【０００２】
人間の呼吸器系の呼吸又は音響インピーダンスは、気道、肺及び胸壁の抵抗、追従及び慣性に関する情報を得るために測定できる。この情報は、慢性閉塞性肺疾患（ＣＯＰＤ）、喘息及び気管支炎のような様々な呼吸器疾患の性質及び発病度を診断する際に有効である。
【０００３】
ページ363-374のEuropean Respiratory Journal Vol.29 No.2のＤｅｌｌａｃａ等による「Expiratory Flow Limitation Detected by Forced Oscillation and Negative Expiratory Pressure」（参考文献１）に説明されているような強制的振動技術（ＦＯＴ）では、人が正常に呼吸している間、音響波が呼吸器系に導かれ、反応が呼吸インピーダンスを決定するために測定される。
【発明の概要】
【発明が解決しようとする課題】
【０００４】
呼吸インピーダンスは、流量及び圧力に関して音響波から生じている振動における周波数依存的な関係を記述する。呼吸インピーダンスが（幾つかの疾患及び他の病状におけるような）吸気から呼気まで変化する所で、呼吸インピーダンスは、精細な時間解像度で推定されなければならない。しかしながら、従来技術では、生理的目的のために充分に短い（すなわち、吸入又は呼気の期間より短い）時間間隔の呼吸インピーダンスを推定するための技術の信頼性に対して、要件がほとんど与えられていなかった。
【０００５】
従って、短い時間間隔の呼吸インピーダンスを確実に推定するための方法と装置とのためのニーズがある。
【課題を解決するための手段】
【０００６】
本発明の第１の態様によると、被験者の気道と連通する患者インターフェイス装置内に圧力波を生成するステップと、流量及び圧力を表わすそれぞれの時系列を作るため被験者の気道と患者インターフェイス装置とを含む空気制御システム内のガスの流量及び圧力を決定するステップと、それぞれの時系列を時間−周波数領域へ変換して変換された時系列を作るステップと、それぞれの変換された時系列から流量及び圧力のパワーを推定するステップと、変換された時系列に基づいて流量及び圧力のそれぞれのクロススペクトルを推定するステップと、推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定するステップとを有する呼吸インピーダンスを推定する方法が提供される。
【０００７】
本発明の第２の態様によると、患者インターフェイス装置と、被験者の気道のガスの圧力、流量又はボリュームの振動を生成するための励振源と、患者インターフェイス装置及び被験者の気道により規定される空気圧回路のガスの流量及び圧力を決定し、流量及び圧力を表しているそれぞれの時系列の値を出力するための手段と、プロセッサとを有し、前記プロセッサは、それぞれの時系列を時間―周波数領域へ変換し、それぞれ変換された時系列から、流量及び圧力のパワーを推定し、それぞれ変換された時系列に基づいて、流量及び圧力のそれぞれのクロススペクトルを推定し、推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定する、呼吸インピーダンスを推定するための装置が提供される。
【０００８】
本発明の第３の態様によると、適切なプロセッサ又はコンピュータによる実行の際、前記プロセッサ又はコンピュータが上述の方法を実施するコンピュータ可読媒体に組み込まれたコンピュータ可読のコードを具備する、コンピュータプログラムが提供される。
【０００９】
呼吸インピーダンスの推定は、気道の妨害を評価するか、又は疾患の発病度を推定するために推定を使用する診断ツールによりなされる。従って、診断ツールは、また、呼吸インピーダンスに影響を及ぼすべき治療の（薬理学的又は他の）効果を評価するためにも用いられる。
【００１０】
よって、第４の本発明の態様によると、上述されたような呼吸インピーダンスを測定するステップと、測定された呼吸インピーダンスを基礎として生理的状態を診断するステップとを有する生理的状態を診断する方法が提供される。
【００１１】
呼吸インピーダンスの推定は、また、又は代わりに、病状の治療に使用される機械、例えば気道閉塞を打ち消すために用いられる非侵襲性のベンチレータの設定を適応させるため、又は、生理的状態に対する治療法（例えばどの特定の薬物又は装置を使用するか、薬物投与量等）を決定する際に用いられる。
【００１２】
よって、本発明の第５の態様によると、上述されたように呼吸インピーダンスを測定するステップと、測定された呼吸インピーダンスに基づいて生理的状態の治療を決定し及び／又は処置するステップとを有する、生理的状態の治療方法が提供される。
【００１３】
本発明の好ましい実施例が、以下の図面を参照して、単なる例示として、ここで詳細に説明されるだろう。
【図面の簡単な説明】
【００１４】
【図１】図１は、本発明による強制的圧力振動から呼吸インピーダンスを推定するための装置のブロック図である。
【図２】図２は、本発明による呼吸インピーダンスを推定するために用いられる線形モデルの模式図である。
【図３】図３Ａ及び図３Ｂは、それぞれ離散的時間領域の巡回性推定及び不確定度を例示する。
【図４】図４Ａ及び図４Ｂは、Δｔ２＋Ｎ２Δｆ２が最小である時系列を例示する。
【図５】図５は、時間―周波数領域の二変量の最小二乗法を例示する。
【図６】図６は、複素平面内の信頼領域を例示する。
【図７】図７は、５つの異なる周波数での強制的圧力振動に対する信頼限界を持つ時間及び周波数の関数として、健常な患者の呼吸インピーダンスを例示する。
【図８】図８は、５つの異なる周波数での強制的圧力振動に対する信頼限界を持つ時間及び周波数の関数として、ＣＯＰＤを持つ患者の呼吸インピーダンスを例示する。
【図９】図９は、本発明による呼吸インピーダンスを推定する方法のステップを図示する流れ図である。
【図１０】図１０は、本発明による装置により実施される処理ステップをより詳細に図示する流れ図である。
【図１１】図１１は、時系列のイベントの数Ｎの関数として、行列の最小及び最大の特異値を示す。
【図１２】図１２は、Ｎ＝１６に対するＣＨＣの固有ベクトルの重心と時間Δｔでの不確定度との間の関係を例示する。
【図１３】図１３は、Ｎ＝１６に対する周波数Δｆでの不確定度と時間Δｔでの不確定度との間の関係を例示する。
【図１４】図１４は、ｆｎに対するνｍｉｎ，ｎをプロットしているグラフである。
【図１５】図１５は、Ｎ＝１６に対するＮΔｆ対Δｔのプロットにおける不確定度を例示する。
【図１６】図１６は、複素平面において推定値Ｚｔ，ｍに対する１００（１−α）％信頼限界の推定を例示する。
【図１７】図１７は、健常な被験者に対する時間及び周波数の関数として二乗コヒーレンスを例示する。
【図１８】図１８は、ＣＯＰＤを持つ被験者に対する時間及び周波数の関数として二乗コヒーレンスを例示する。
【発明を実施するための形態】
【００１５】
上述のように、呼吸器系の音響インピーダンスは、マウスピース又はフェースマスクの強制的圧力振動（強制的振動技術、ＦＯＴ）から推定されることができる。この「呼吸インピーダンス」は気道、肺及び胸壁の抵抗、追従及び慣性に依存するので、呼吸インピーダンスは、様々な呼吸器疾患の発病度及び性質に対する洞察を提供する。周波数依存的なインピーダンスは、しばしば時間的に変動する。慢性閉塞性肺疾患において、気道抵抗は、呼気の間に通常増大する。睡眠無呼吸又はいびきにおいて、抵抗は、吸気の間に増大する。これは、短い時間間隔（呼吸の期間より短い）におけるインピーダンスの信頼性が高い推定を要する。
【００１６】
呼吸インピーダンスは、周波数ｆの関数として、流量ｑと圧力ｐとの間の線形関係を記述する複素数値の伝達関数として規定される。
ｐ（ｆ）＝ｑ（ｆ）ｚ（ｆ）（１）
ここで、インピーダンスはｚ（ｆ）により示される。この式は、ｚ（ｆ）が一定で、時間が無限である場合だけ成り立つ。システムが一時的に安定の場合、ｚ（ｆ）は短い時間間隔で依然推定される。しかしながら、これは、信号が時間及び周波数領域両方において鋭くは特定できないと言う不確定性原理により制限される。加えて、短い時間間隔の使用は、推定を偶然のイベントに益々依存させる。問題は、これが、呼吸サイクルの一部で安定なだけであるｚ（ｆ）のような伝達関数の推定にどのくらい影響を及ぼすかである。
【００１７】
以下では、一時的に安定なプロセスの伝達関数が、離散的時間―周波数領域の線形回帰により推定される。不確定性関係は、離散的時間及び周波数に対して導出され、最適時間―周波数解像度を得るために用いられる。使用されるインピーダンス推定器の統計特性の分析は、呼吸メカニクスの周波数依存的なパラメータを得て、その推定された値及び信頼限界両方が時間的にフォローできる。
【００１８】
発明の詳細な説明の後に含まれるテーブル（表１）は、本発明の以下の説明で使用されるシンボル及び略語一覧を提供する。
【００１９】
方法
発明の詳細な説明のこのセクションは、本発明による方法及び装置の数学的基礎を説明しておく。本発明のより数学的でない説明は、以下の「説明」セクションにおいて提供される。
【００２０】
本発明による装置２の模式図が図１に示される。簡潔には、使用されるＦＯＴ装置２は、被験者の気道を含む空気制御システムを作るために被験者の気道に接続される患者インターフェイス装置４から成る。患者インターフェイス装置４は、マスク又はマウスピースのような被験者の気道との気送接続を供給するために適切な任意の装置である。患者インターフェイス装置４は、比較的低い周波数、例えば８Ｈｚから２４Ｈｚまでの周波数を持つ圧力波、流量変化又はボリューム変化を生成する励振源６に動作的に結合される。本発明の例示的な実施例では、患者インターフェイス装置４は拡声器である。
【００２１】
図示の実施例では、患者は、患者インターフェイス装置４と励振源６との間に位置されるニューモタコヘッド８に近いメッシュ有線の抵抗９を通じて空気を呼吸する。強制的波は、患者の気道及び肺において部分的に反射され、ＦＯＴ装置２において測定されるエアフロー及び圧力の振動を導く。ニューモタコヘッド８を通るエアフローは、差圧トランスジューサ１０を使用して測定され、圧力はマウスピース４の近くの第２の圧力トランスジューサ１２を使用して測定される。ＡＤコンバータ１２は、流量及び圧力のそれぞれの時系列｛ｑｔ｝｛ｐｔ｝を作り、分析のためプロセッサ又はコンピュータ１６に供給され、ここで、整数ｔは離散的時間指標である。コンピュータ又はプロセッサ１６は、ＡＤコンバータ１４及びアンプ１８を介して拡声器６により生成される音響波の周波数を制御するための信号も供給できる。装置２の特定の実施例の更なる詳細は、「測定装置」と名付けられた以下のセクションで提供される。
【００２２】
本発明は、また、他のガス流送出装置が、励振源６に加えて被験者の気道に結合できることを考察する。例えば、圧力サポートシステム又はベンチレータは、例えば、励振源６が被験者の気道を振動させている間、ガスの流れを供給するために患者インターフェイス装置４に結合できる。また、圧力サポートシステム又はベンチレータ内で圧力及び流量を測定することを含む、流量及び圧力測定が、空気制御システムに沿って任意の位置でなされ得ることに留意されるべきである。その後、測定された圧力及び流量は、任意の従来の技術を用いて患者の気道のような空気制御システムのガスの流量及び圧力を決定するために使用できる。
【００２３】
本発明は、また、励振源６が独立型装置よりもむしろ圧力サポートシステム又はベンチレータの部品であることも考察する。例えば、圧力サポートシステム又はベンチレータは、患者へ送られるガスの流量／圧力を制御するためのバルブを含む。本発明は、低周波圧力波、流量変化又はボリューム変化をもたらすために斯様なバルブを振動させることを考察する。もちろん、低周波圧力波、流量変化又はボリューム変化を生成できるピストンのような他の任意の装置又はシステムも、励振源６として使用できる。
【００２４】
圧力及び／又は流量センサを使用して圧力及び／又は流量を測定するよりはむしろ、本発明は、また、流量又は圧力が、圧力サポートシステム又はベンチレータによって設定又は制御できることも考察する。この場合、設定された圧力又は設定された流量は、差圧トランスジューサ１０又は第２の圧力トランスデューサ１２の出力としてよりもむしろ現目的のための圧力又は流量として使用される。
【００２５】
２つのサンプルの間の間隔（秒）は、サンプルレート（Ｈｚ）の逆数である。このセクションでは、時間依存及び周波数依存インピーダンスが、比較的単純な線形代数学を使用してこれらの時系列から導かれる。線形代数学をよく知らない読者は、主要な結果が言葉で要約されている以下の「説明」セクションを参照されたい。
【００２６】
単純な線形モデル偶然のイベントを説明するために、二変量の時系列｛ｑｔ、ｐｔ｝が二変量の確率プロセス｛Ｑｔ、Ｐｔ｝の実現であると仮定する。決定論的な変数は小文字で書かれ、ランダム（確率）変数（ＲＶ）は大文字で書かれる。時間が無限に長い（ｔ＝．．．，−１，０，１．．．）と更に仮定する。ここで図２を参照すると、励振源入力ｘｔは、インパルス反応シーケンス｛ｈｑｘ，ｌ｝及び｛ｈｐｘ、ｌ｝を持つ線形時間不変のフィルタ（２及び３）を通ってフィルタリングされる。これは、以下の線形モデルにより予測されるように、流量及び圧力を与える。
（２）
【００２７】
これらの「真」の値がそれぞれの分散σＱ２及びσＰ２を持つゼロ平均ホワイトノイズの２つの独立源により妨げられると仮定する。ノイズのこれらの源は、「エラー」項Ｑｅ，ｔ及びＰｅ，ｔを与えるために、線形時間不変のフィルタ（１及び４）も通過する。完全なモデルは、このとき、以下のように書かれる。
（３）
ここで、平均値μＱ及びμＰは定数である。結果として生じるＱｔ及びＰｔが正常に分配されるので、仮定された二変量のプロセス｛Ｑｔ、Ｐｔ｝は、完全に変動しない（静止している）。
【００２８】
循環性仮定
説明された静止プロセスが無限に長い一方、静止がＮ個のイベントの有限シーケンスに対してのみ仮定できている過渡現象を説明することが発明の目的である。これは、最後のイベントが最初のイベントに先行するという点で、時間が循環していると仮定することにより解決できる。結果として生じる時系列は、このとき、時計のような円に沿って定期的にＮイベントを配置することにより表される。円から離れずに、同一方向で円を永久に追求できる。結果として生じる無限の時系列は、Ｎで周期的である。すなわち、任意の整数ｊに対して、以下の通りである。
ｘｔ＝ｘｔ＋ｊＮ （４）
【００２９】
ｕ＝ｔ−ｌ＋ｊＮと置き換えると、式２の第１のコンボリューションは、以下のように書きかえられる。
【００３０】
仮定された周期的キャラクタｘｔ（式４）に対して、ｘｕ−ｊＮ＝ｘｕである。上記の式の無限和は、「長さＮに周期化されたｈｑｘ，ｔ−ｕとして、以下の式のように示されて規定できる。
（５）
結果として、
（６）
【００３１】
ｈｑｘ，ｔ−ｕｏもＮにおいて周期的であるので、（インデックスはゼロからＮ−１の範囲に常にある（これは便利である）ように）ｔ−ｕが負である場合、インデックスｔ−ｕはＮの分増大できる。シーケンスが列ベクトルにより表されるとき、式６は以下のように書ける（Ｎ＝４である場合）。
（７）
【００３２】
この式のＮｘＮ行列は、巡回行列式（端で循環しながら、前の行が右にシフトするバージョンの行を持つ正方行列）である。第１の行がｈｑｘＨにより示されるとき、肩文字Ｈは、エルミート転置を表わし、巡回行列はｃｉｒｃ｛ｈｑｘＨ｝とも書かれる。
、
及び
とする。このとき、ｑｐ＝Ｃｑｘｘである。同様に、ｐＰ＝Ｃｐｘｘなので、式３は以下のようになる。
（８）
【００３３】
離散的周波数領域
周波数依存的なインピーダンスを導くために、式８のモデルは、離散的周波数領域へ変換される必要がある。時間―周波数分析へのステップアップとして、これは、ここで短く要約される。所与の時間系列
から離散的周波数領域への変換は、時系列を離散的周波数
（ｎ＝０，．．．，Ｎ−１））を持つ一組の調和振動に区切る。斯様な振動は、以下のようなＮ―ベクトルにより記述される。
（９）
ここで、ω≡ｅｘｐ（２πｉｎｔ／Ｎ）及びｉ２＝−１である。オイラーの定理により、
ｅｘｐ（２πｉｎｔ／Ｎ）＝ｃｏｓ（２πｎｔ／Ｎ）＋ｉｓｉｎ（２πｎｔ／Ｎ） （１０）
よって、ベクトルｆｎは、実際に、周波数ｆｎを持つ実軸及び虚軸に沿った複合振動を説明する。容易に以下の式が示される。
（１１）
【００３４】
結果として、二乗ノルム
はＮに等しい一方、異なる調和周波数の２つの振動は直交する。
【００３５】
ｘの離散的周波数領域への変換は、以下の式により規定される「正規直交離散的フーリエ変換」（ＯＤＦＴ）行列により実施される。ここで、ｎｔ＝０，．．．，Ｎ−１である。
（１２）
この変換は、Ｎ―ベクトルＦｘを作り、これらの成分は、内積
である。式１１から、Ｆの行（及び列）が正規直交であるので、行列はユニタリであることがわかる。
ＦＦＨ＝ＦＨＦ＝Ｉ （１３）
ここでＩは、ＮｘＮの単位行列である。これは、Ｎ振動からｘの合成を可能にする。
（１４）
【００３６】
従って、ｘは、周波数ｆｎを持つＮの調和振動の重み付けられた和として書かれる。
【００３７】
ベクトルｆｎは巡回行列の固有ベクトルであるので、以下のように書くことができる。Ｃｑｘｆｎ＝ｆｎＨｑｘ，ｎ （１５）
【００３８】
複素数Ｈｑｘ，ｎは、以下の式により規定されるように、周波数ｆｎに対してｘからｑまでの伝達関数である。
（１６）
【００３９】
Ｈｐｘ、ｎが同様に規定されるとき、周波数ｆｎに対するインピーダンスは、以下の式である。
（１７）
【００４０】
固有値が対角行列
に入れられ、列として対応する固有ベクトルをＮｘＮ行列ＦＨへ入れられるとき、式１５は以下の式になる。
（１８）
【００４１】
巡回行列は、正規行列である（これらは、これらのエルミート転置行列と代えられる）。従って、式１８は、ＮｘＮ正規行列がＮ個の正規直交固有ベクトル（この場合、Ｎ個の正規直交列ベクトルｆｎ）を持つというスペクトル定理と一致している。
【００４２】
離散的時間及び周波数に対する不確定性原理
の成分は、時間領域において鋭く局所的であるが、周波数領域においては不確かである。これとは逆に、フーリエ変換Ｆｘの成分は、周波数領域において鋭く局所的であるが、時間領域においては不確かである。両方の領域で小さな不確定度を持つ信号を得るために、不確定さの尺度が、規定されなければならない。
【００４３】
時間が巡回性であると仮定すると、フォーブズ及びアロンゾのアプローチが、最も適当なようである（参照文献３を参照）。ここで図３Ａを参照すると、時間領域で、二乗成分｜ｘｔ｜２は、各イベントの重みに対応する図３Ａの各質点のサイズを持ち、原点で中央に置かれる、複素平面内の半径ｒｔを持つ円上の規則的な間隔の質点として置かれる。時間は、２つのイベントの間の円弧長により表される。円の周囲は全体の期間Ｎ（この具体例ではＮ＝１２）であるので、半径は
に等しい。時刻ｔでの所与のイベントは、複素平面のポイントｒｔωｔに位置される。フォーブズ及びアロンゾ（参照文献３）は、不確定度の幾つかの尺度を提案した。ここで、わずかに異なる尺度が使われ、これは、時間ｔでの各イベントと基準時間ｔＲでのイベントとの間の二乗距離の加重平均として式１９に規定される。
（１９）
【００４４】
不確定度のこの尺度は、二次形式の比率として容易に表されることができるという利点を持つ。一般性の損失なしに、ｔＲ＝０（図３Ａに示されるように）と仮定できる。このとき、二乗された不確定度Δｔ２は、以下の式である。
（２０）
ここで、対角行列Ω≡ｄｉａｇ｛１，ω１，．．．，ωＮ−１｝である。下記の式を規定することも可能である。
（２１）
ここで、Ａ≡Ω−Ｉである。物理的に、Δｔ２は、基準ポイントＲ（時間ｔＲに対応する）を通る（円の平面と直角をなす）軸の周りの質点のセットの第２の慣性モーメントである。図３Ｂに図示されるように、この慣性モーメントは、重心に直接関係する。幾何学的に、Δｔは、Ｒから、重心を通る垂直線が円と交差する２つのポイントのうちの一方のポイントまでの距離に等しい。それで、Δｔは、重心からＲまでの水平距離により、完全に決定される。重心がＲと一致する場合、ｘの重さは時間ｔＲに完全に集中し、Δｔは、ゼロである。重心が円の反対のポイントにある場合、重さはそのポイントに完全に集中し、Δｔは最大であり、２ｒｔに等しい。よって、図３Ｂにおいて、重心Ｂを通る垂直線は、ポイントＣで円と交差する。不確定度Δｔは距離
と等しい。
【００４５】
周波数領域において、Ｆｘの成分は、周波数ｆｎ≡ｎ／Ｎと関連している。同じ成分が周波数ｎ／Ｎ＋ｊに対して得られ、ここで、ｊは任意の整数である（式９、式１０を参照）。従って、Ｆｘの成分は、単位ユニット期間を持って、周波数領域において周期的である。このように、不確定度の匹敵する尺度は、基準周波数ｆＲに関係して、周波数領域で規定できる。Ｆｘの成分の二乗の大きさは、複素平面の円上の質点として、再び配置される。ここで、半径は、ｒｆ≡１／（２π）に等しいので、円周はユニティ、１に等しい。ｆＲ＝０の場合、下記の式である。
（２２）
【００４６】
対角行列Ωは、「時間シフトオペレータ」Ｔ≡ｃｉｒｃ｛０，．．．，０，１｝に関係する（２２）。Ｔによる左からのｘの乗算は、巡回的にｘの成分を一つ下方へシフトさせる。式１８から容易にＦＴ−１＝ΩＦ（「シフト定理」）であることがわかるので、以下の式となる。
【００４７】
Ｆがユニタリであるので、以下の式となる。
とすると、以下の式となる。
（２３）
【００４８】
ここで、以下の行列式を考える。
【００４９】
時間及び周波数の全体の不確定度は、以下の比率の式により表わされる。
（２４）
【００５０】
この比率（レイリー比率）は、Ｃの最小二乗特異値から最大二乗特異値までの範囲の値をとることができる。これは、離散的時間及び周波数に対する対応する不確定性原理を、以下の式のように作る。
（２５）
【００５１】
従って、全体の不確定度Δｔ２＋Ｎ２Δｆ２が少なくともσｍｉｎ２に等しくなければならないので、Δｔ及びΔｆ両方が非常に小さいことは不可能である。全体の不確定度は、ｘがσｍｉｎ２に対応するＣＨＣの固有ベクトルである場合、最小である。ＣＨＣの固有ベクトルは、離散的マチュー関数として考えられる。ＣＨＣの興味深い特性は、Ｆと代われることであるので、ＣＨＣとＦとは、同じ固有ベクトルを共有する。ｖｍｉｎに対応するＦの固有値はユニタリに等しいので、Ｆｖｍｉｎ＝ｖｍｉｎである。図４Ａに示されるように、これは、ｘがｖｍｉｎと等しい場合、Δｔ＝ＮΔｆを意味する。図４は、合計Δｔ２＋Ｎ２Δｆ２が最小である時系列を示す。図４Ａでは、固有ベクトルは、基準時間ｔＲ＝０及び基準周波数ｆＲ＝０（上のグラフ）に対して導出される。これは、フーリエ行列（単位固有値を持つ）の固有ベクトルでもあるので、時系列は離散的周波数領域（下のグラフ）への変換後、不変のままである。これらのグラフにおいて、時間軸及び周波数軸の両方とも、視覚的理由のために巡回的にシフトされる。図４Ｂは、基準周波数ｆＲ＝１／４である以外、図４Ａに対応する。図４Ｂでは、黒丸が実数値を表し、白丸が虚数値を表し、十字が絶対値を示す。
【００５２】
不確定性原理のより詳細な分析のため下記のセクション「不確定性原理の導出」を参照されたい。
【００５３】
離散的時間―周波数領域
ＯＤＦＴは、Ｎ個のイベントの時系列をＮ個の異なる周波数を持つ一組の振動へ区切られる。時間―周波数分析を実施する簡易なやり方は、短く続くシーケンスのＭ（Ｍ＜Ｎ）個のイベントに対してＯＤＦＴを繰り返し計算することである。従って、系列が周波数ｆｍ≡ｍ／Ｍ（ｍ＝０，．．．，Ｍ−１）を持つＭ個の振動に区切られ、これは離散的「短時間フーリエ変換」になる。各短く続く振動がｖｍｉｎの成分で「窓」をつけられる（ウィンドウ化される）場合、時間及び周波数についての結果的な不確定度は最小である。
【００５４】
Ｍ個のイベントのシーケンスに対して、斯様なウィンドウ化された振動は、以下の式のようにＭベクトルにより記述される。
（２６）
【００５５】
例が、図４Ｂに示される。ｖｍｉｎ（ｍ）は、周波数領域において実際にｖｍｉｎの巡回的にシフトされたバージョンであることに留意されたい。図４に例示されて、下記のセクション「不確定性原理の導出」で証明されているように（式Ｂ１９を参照）、これは、（ｆｍが基準周波数として採用される場合）固有ベクトルと同じ時間及び周波数の不確定度を持つ。ｖｍｉｎ（ｍ）の成分をｖｍｉｎ，ｔ（ｍ）により示すとする。短く続く振動を記述するＮベクトルは、このとき、ｖｍｉｎ（ｍ）の中間にゼロを挿入することにより得られる（よって、ベクトルの単一モードの構造が完全なままである）。これは、以下の式を与える。
（２７）
ここで、ＫはＭ／２以下の最も小さな整数である。この振動は、時間ｔ＝０及び周波数ｆｍ＝ｍ／Ｍを中央に置く。一般に、以下の式
（２８）
の振動は、時間ｔ及び周波数ｆｍを中央に置く。
【００５６】
最小二乗法という意味では、ｘとｈｍ，ｔとの間の最も適合する線形関係は、ｈｍ，ｔによりスパンされた（複素）一次元の部分空間上へのｘの直交射影からわかる。ｈｍ，ｔが単位ベクトル（式２７を参照）として規定されることを考慮すると、この射影は、以下の式となる。
（２９）
【００５７】
所与の周波数ｆｍに対して、ＮｘＮ行列Ｍｍを以下の式のように規定する。
（３０）
【００５８】
これは、共振周波数ｆｍを持つ帯域フィルタを行う巡回行列である。Ｍｍによる左からのｘの乗算は、波形符号又は「波符号」により示されるＮベクトルを与える。
（３１）
【００５９】
の成分は、内積
である。
【００６０】
十分な時間―周波数変換（ＴＦＴ）は、ｍ＝０，．．．，Ｍ−１に対する全ての周波数ｆｍを含む。ＭＮｘＮ行列のＭを以下の式のように規定する。
（３２）
【００６１】
離散的時間―周波数領域（時間―周波数領域）への変換は、（所与のＭに固有な）Ｍによる左からのｘの乗算により実施される。これは、ＭＮ―ベクトル
を与える。バック変換は、再びＭＨによる乗算を意味し、
を与える。下記の「時間―周波数変換の導出」というセクションにおいて、Ｍの列は正規直交であることが示されるので、以下の式となる。
ＭＨＭ＝Ｉ（３３）
【００６２】
ＭＨＭｘ＝Ｉｘ＝ｘなので、よって、バック変換は元の時系列を回復する。これは、以下の式に従って、ｘの時間―周波数合成を可能にする。
（３４）
【００６３】
このように、時系列はＭＮ個の短く続く振動の加重加算として書かれ、各々が時間ｔ及び周波数ｆｍに中心が置かれる（式１４を参照）。所与のＭに対して、これらの振動は、式２５の不確定性原理に従って、最適時間―周波数解像度を与える。下記の「不確定性原理の導出」セクションに示されるように、振動は、Ｍ’＝１０・Δｔゼロ以外の要素で丸められたガウス分布により近似できる。
【００６４】
励振源への入力
励振源入力ｘとして、一組の調和振動を使用することが便利である。データを分析するために、このとき、これらの振動の周波数がＴＦＴ共振周波数ｆｍと一致するように、ＴＦＴのスケールＭが選択できる。最初に、周波数ｆｎ及び振幅ａｍを持つ単一の正弦波入力を考える。オイラーの定理（式１０）により、実数値の調和振動は、複素数値とその複素共役との和である。よって、励振源入力は、下記の式のように書ける。
（３５）
ここで、星印は複素共役を示す。式９及び式１０から、ｆｎ＊＝ｆＮ−ｎがわかる。ＮがＭの整数倍である場合、周波数ｆｎ≡ｎ／Ｍは、ＴＦＴ共振周波数ｆｍ≡ｍ／Ｍとマッチできる。その後、ｆｎ＝ｆｍであるように、ｍが選択される場合、以下の式となる。
（３６）
【００６５】
Ｍが十分に大きく、対応するΔｆが、共振周波数ｆｎ＝ｆｍを持つＴＦＴフィルタにより周波数ｆＮ−ｎで振動を抑制するのに十分小さい場合、後の近似値は有効である。「相補的周波数」ｆＮ−ｎ（ここでは、ｆＭ−ｍに等しい）での成分も式３６でＭＭ−ｍをＭｍに加算することにより含まれる場合、時間―周波数領域の励振源入力の十分な表現が得られる。
【００６６】
式８のモデルは、Ｍにより両方の式を左から乗算することにより、時間―周波数領域に変換される。対象の周波数（例えば、式３５により記述された強制的振動の周波数）に注目すると、Ｍｍにより式８の上の部分を乗算することは、以下の式となる。
（３７）
ｍがゼロでない場合、ＭｍμＱがゼロであることが容易に示される。よって、
（３８）
【００６７】
巡回行列全てが代えられる
であることがわかる。ｆｎ＝ｆｍであるようにｎ及びｍが選ばれる場合、このとき、式３６及び式１５を用いて、以下の式となる。
（３９）
【００６８】
ｆｎ＝ｆｍに対するＨｑｘ、ｎがｂＱ、ｍにより示される（対応するＨｐｘ、ｎがｂＰ、ｍにより示される）とき、式８のモデルは、以下の式に低減される。
（４０）
【００６９】
ｆｎ＝ｆｍに対するインピーダンスは、以下の式となる。
ｚｍ≡ｂＰ、ｍ／ｂＱ、ｍ （４１）
【００７０】
これらの振動が時間―周波数領域においてよく離隔されるほど、使用されるＴＦＴのΔｆが十分に小さいならば、他の振動のセットは、式３５の励振源入力に加えられる。
【００７１】
時間―周波数領域の分散分析
ＴＦＴは、時間―周波数領域の分散の分析を可能にする（「時間―周波数領域の分散分析」というこのセクションを参照）。例えば、流れのノイズの二乗ノルム（式８を参照）は以下の式に区切られる。
（４２）
【００７２】
のランダム変数は、ノイズＱｅの「ＴＦＴサンプルパワースペクトル」と呼ばれている。転送がノイズ１からＱｅ，ｔ（図２）までのフィルタ１の伝達関数が、ｆｍについてのＴＦＴフィルタの通過帯域で相対的に平坦であると仮定する。このとき、
は、あたかもゼロ平均及び分散σＱ、ｍ２Ｉを持つホワイトノイズから導出されるかのように、ほぼ分散される。結果として、
は、等価自由度ηを持つ近似のカイ二乗分布に従う。（相補的周波数ｆＭ−ｍの成分と同様に）循環仮定から独立している
のそれらの成分だけが含まれているηの導出のため、以下のセクション「時間―周波数領域の分散分析」を参照されたい。
【００７３】
時間―周波数領域の二変量の最小二乗法
各時間ｔで、図２のモデルにより記述されている、しばらくの間変動しないプロセスが始まると、ここで仮定される。当該プロセスは、Ｎ個のイベントを含み、循環的である（すなわち、Ｎで周期的）とみなされる。時間―周波数領域では、モデルは（式４０から）以下の式となる。
（４３）
ここで、ベクトルは時間ｔで始まる「短く続く」Ｎ個のベクトルである。定数ｂＱ、ｍ及びｂＰ、ｍは、これらがプロセスが始まる時間に特有であるとみなされるので、同様にインデックスｔを得た。従って、モデルへの入力と同様に定数は、（２つのプロセスが時間的に重複する場合であっても）異なる時間ｔで異なる。
【００７４】
ｂＱ、ｔ、ｍ及びｂＰ、ｔ、ｍの最も適切な推定値は、エラーベクトルの二乗ノルムを最小化することにより、
及び
から導出できる。これは、最小二乗推定値
及び
を生じる。推定された関係は、以下の式である。
（４４）
ここで、推定された「予測された」ベクトル
は、推定された「エラー」
に正規直交している（同様に、
は
と正規直交している）。図５では、ベクトルは矢印として示され、全ての変数は時間インデックスｔ及び周波数インデックスｍ（周波数ｆｍ＝ｍ／Ｍに対応する）を参照している。特に、図５では、励振源入力
は一定の変数である一方で、流量
及び圧力
はランダム変数である。ベクトル
及び
は、
によりスパンされた複素一次元の部分空間（ラインｌ）への
及び
の射影である。インピーダンス推定値は、
から分かる。エラーベクトル
及び
は
に垂直である。各ベクトルは、長さがノルムに等しい矢印により表される。ベクトル間の角度は、ｃｏｓ２φ及びｃｏｓ２θが、
とそれぞれ
及び
との間の二乗コヒーレンスに等しいように描かれる。
【００７５】
と
との間の線形関係の程度は、以下の式のような「ＴＦＴは二乗サンプルコヒーレンス」の表現になる。
（４５）
【００７６】
ランダム変数
は、ｘからＱまでの「ＴＦＴサンプルクロススペクトル」と呼ばれている。幾何学的に、ＫＱｘ，ｔ，ｍ２はｃｏｓ２φと等しく、ここで、φは図５の
から
までの角度である。同様に、ＫＰｘ，ｔ，ｍ２はｃｏｓ２θと等しい。最小二乗推定値は、標準的やり方で導出される。これは、以下の式を与える。
（４６）
【００７７】
従って、
は、ｘからＱまでのサンプルクロススペクトルとｘのサンプルパワースペクトルとの比率である。
が同様のやり方で導出されるとき、インピーダンスは以下の式により推定される。
（４７）
【００７８】
「真の」インピーダンスｚｔ、ｍは、実数部分ｒｔ、ｍ（「抵抗」）及び虚数部分ｘｔ、ｍ（「リアクタンス」）を持つ。対応する推定値は、以下の式、
の実数部分及び虚数部分である。
（４８）
【００７９】
図６Ａに示されるように、ｂＱ、ｔ、ｍに対する１００（１−α）％信頼領域は、
を中心にした複素平面の円により区切られている。信頼領域は、以下の式により与えられる。
（４９）
ここで、
であり、ＦαはＦ分布の上の１００（１−α）％である（以下のセクション「信頼限界の導出」を参照）。励振源入力ｘｔは、無次元変数であり、
は（図５から
なので）流量のユニットを持つ。循環仮定に依存する値は、ここで放棄される（サンプルの有効数は、ＮＳ＝Ｎ−Ｍ＋１である）。原点が円内にある場合（このとき、
は１００・α％信頼レベルでゼロと著しくは異ならない）、関係は拒絶される。よって、ＡＱ、ｔ、ｍ＜１の場合、有意性が得られることが、式４９から分かる。
に対する信頼領域は比較可能な円により区切られ、
と
との間の関係は、ＡＰ、ｔ、ｍ＜１の場合、重要であると考えられる。
【００８０】
に対する信頼領域の控え目な推定値は、
及び
が非相関であるという仮定の下、以下の「信頼限界の導出」セクションで導出される。当該領域は、以下の式により記述される。
（５０）
【００８１】
図６Ｂを参照すると、限界円は推定値
の周りに同心円ではない。
の実数部分及び虚数部分が、
及び
により表わされる。垂直矢印は、実数部分に対する上位の信頼限界を示し、水平矢印は、インピーダンスの虚数部分に対する下位の信頼限界を示す。
の所与の実現のための「真の」ｚｔ，ｍの分布は、明らかに斜めに対称形である。図６Ｂでは、歪度は角度βにより決定される。歪度は、図５のβと角度φとの間の関係から明らかなように、完全に流量のノイズに起因する。
（５１）
【００８２】
流量のノイズがゼロである場合、このとき、φ＝０、ＫＱｘ、ｔ、ｍ２＝１及びＡＱ、ｔ、ｍ２＝０であるので、β＝０であり、信頼領域が
について対称形である。
【００８３】
下記の「信頼限界の導出」セクションは、流量のノイズ又は圧力のノイズがゼロであると仮定される場合に起こることの更なる分析を与える。それぞれの推定値は、このとき、以下の式に等しい。
（５２）
【００８４】
これは、時間―周波数領域における流量から圧力まで（又はその逆）の線形回帰となる。変数
及び
は、それぞれのサンプルパワー及びクロススペクトルである。式５２の推定値と
との主要な違いは、推定されたインピーダンスの大きさにある。大きさは、以下の式を通じて関係付けられる。
（５３）
【００８５】
ノイズが流量及び圧力に存在する場合、インピーダンスの大きさは、
が
の代わりに使用される場合過小評価される（
が使用される場合過大評価される）。
【００８６】
設定
装置２の使用中、流量及び圧力が、８００Ｈｚのサンプルレートで記録された。ＴＦＴは、Δｔ＝５０で実施された（これは、５０／８００＝０．０６２５ｓに達する）。これは、Ｍ＝４πΔｔ２＝３１，４１６のＴＦＴフィルタの幅を必要とする。しかしながら、ＴＦＴフィルタの係数は、Ｍ’＝１０・Δｔ＝５００で丸められたガウス分布により近似された（以下の「不確定性原理の導出」セクションを参照されたい）。周波数の関連する不確定度は、Δｆ＝１／（４πΔｔ）＝０．００１６である（又は、０．００１６ｘ８００＝１．２７Ｈｚ）。ＴＦＴは、８、１２、１６、２０、２４Ｈｚの課されたＦＯＴ周波数に対応して、ｆｍ＝０．０１、０．０１５、０．０２、０．０２５、０．０３に対して導出された。サンプルパワー及びクロススペクトルは、ＮＳ＝５００に対して（すなわち、５００／８００＝０．６２５ｓの時間間隔にわたって）計算された。循環仮定から独立しているデータだけが使われた。これは、η＝６．７２を持つサンプルスペクトルの近似のカイ二乗分布を導く。インピーダンスに対する信頼限界は、９０％の信頼レベル（Ｆα＝４．０５）に対して導出された。時間―周波数スペクトルは、最初の時系列と時間的に位置合わせされた。
【００８７】
例
図７は、通常の被験者に対する時間及び周波数の関数として、推定されたインピーダンス
の実数部分及び虚数部分を示す。青いラインは実数部分
を表し、赤いラインは虚数部分
を表す。グレイのラインは、図６Ｂに図示されたように、信頼限界を表す。８Ｈｚの周波数帯域の右側の軸は、インピーダンスの成分のスケールを示す。８Ｈｚで、実数部分
は正であり、呼吸サイクルの間で変動する一方、虚数部分
はほぼゼロに等しく、比較的一定のままである。信頼限界は、時間の経過で、
及び
の変化の統計的有意性について直観的な印象を与える。より高い周波数で、
は、
が徐々により正になっていくと共に、減少する。被験者が要求により吸い込むとき、信号は激しく妨げられ、これは信頼領域の大きい増加が伴われた。関連するコヒーレントなスペクトルは、下記のセクション「信頼限界の導出」に示される。
【００８８】
図８は、呼気の間、
の大きな負のスイングを持つＣＯＰＤを持つ患者に対するインピーダンスを示す（これは、明らかに、信頼限界からみて吸気値と著しく異なる）。斯様な負のスイングは、以前にＣＯＰＤにおいて記述されていて、呼気の間、気道の準最大の虚脱のため、流量制限に関係する。おそらく患者自身の呼吸の高周波成分のため、強制的振動を妨げるので、信頼領域は、吸気から呼気への転換点で最大であり、またこの逆も成立する。干渉性スペクトルのために下記のセクション「信頼限界の導出」を参照されたい。
【００８９】
説明
図９は、本発明による方法のステップを要約するフローチャートである。この例示的実施例では、圧力波は、被験者の気道と連通する患者インターフェイス装置４としてマウスピースに接続された拡声器を使用して生成される（ステップ１０１）。マウスピース４を通るガスの流れ及び圧力の振動が、測定のそれぞれの時系列を与えるために測定される（ステップ１０３）。定常状態において、呼吸インピーダンスは、静かに呼吸している被験者のマウスピース４の強制的圧力振動から推定できる。呼吸器系の機械的特性が吸気から呼気にしばしば変化するので、当該方法は一時的な「安定性」の状況の下でインピーダンスを推定する。主要な課題は、インピーダンスが周波数依存的な量であるということであり、そして、高い時間的解像度が必然的に低い周波解像度に至るということであり、これによって、推定を歪める。データが離散的時系列（マウスピース４の流量及び圧力のシーケンシャル測定）から成るので、これは短い時系列の最適時間―周波数分析を求める。
【００９０】
上記のセクション「方法」にて説明されたように、本発明は、時間が循環している（下記参照）とみなされる場合、短い時系列が変動しないという概念に基づく。説明されたように、時間及び周波数領域の全体の不確定度に下位限界を配置する、不確定性原理の新しいバージョンが、斯様な時系列に対して導出された。ステップ１０５では、この原理は、各時系列の相違を特定の時間及び周波数と関係している成分に区切るために用いられる（「離散的時間―周波数領域」又は「時間―周波数領域」への変換を通じて）。その後、時間―周波数領域の最小二乗法分析は、最適時間―周波数解像度を持つインピーダンスのバイアスされていない推定値を与える（ステップ１０７）。最後に、インピーダンスに対する信頼限界が、時間及び周波数の関数として作られる。主要なステップが、以下に要約され説明される。
【００９１】
線形モデル
呼吸インピーダンスの推定値は、単純な線形モデルに基づく（式３、図２）。装置２のマウスピース４の振動は、拡声器６により誘発される。流量及び圧力の結果として生じる変化は、患者の呼吸器系（その系の音響効果）の機械的特性に依存する。モデルにおいて、これは、拡声器入力を修正する２つの線形時間不変のフィルタ（図２の２及び３）により表される。流量及び圧力の結果として生じる変動は、流れ及び圧力両方のノイズの２つの独立源により妨げられる。呼吸インピーダンスは、２つの線形フィルタ（２及び３）により決定される。
【００９２】
ノイズの可能性がある源は、以下の通りである。１）強制的振動の周波数範囲内の患者の呼吸の成分。これは、おそらく最も重要な課題である（参考文献２を参照）。２）患者による吸い込み、咳嗽、又は他の動き。３）心臓性振動（特に気道抵抗が低いとき）。４）呼吸インピーダンス自体の時間依存性の変動。５）ランダムな測定エラー。
【００９３】
呼吸器系の非線形インタラクションが無視できないという指標があったが、線形モデルは、以前、同様の実験（参考文献２を参照）の間、流量と圧力との間の周波数依存的な関係を説明するために用いられた。しかしながら、ＦＯＴ装置２の流量と圧力との間の関係の線形記述は、圧力差が平均絶対圧と比較して小さいという事実によりサポートされる。当該関係を記述するためにインピーダンスを使用することは、これらの状況の下で期待できる定常波のインタラクションの簡略化である。しかしながら、主要な制限は、呼吸器系が安定と考えられる時間にある。
【００９４】
循環仮定
偶然のイベントを説明するために、流量及び圧力が所与の確率分布を持つランダム変数（ＲＶ）であると仮定された。よって、流量及び圧力のＮ個の対の測定値の短いシーケンスは、有限二変量の確率過程（確率プロセス）からのサンプルとして考えられる（Ｎ個の経時的に並べられた対のＲＶのセット）。斯様なプロセスにおいて、流量と圧力との間の周波数依存的な関係は、二変量の確率プロセスが（２次のオーダーで）安定していると仮定される場合、安定なインピーダンスによってのみ記述できる。これは、同時に起こる値と連続した値との間の共分散だけでなく、流量及び圧力の期待値が一定であることを意味する。しかしながら、連続した値の間の共分散は、開始及び終了効果のため、（これらの共分散がゼロに等しくない限り）有限シーケンスに対して一定でありえない。
【００９５】
しかしながら、この課題は、最後のイベントが最初のイベントに先行するという意味で、時間が循環すると仮定することにより解決される。このようにして、連続した値間の共分散は一定である一方、シーケンスは依然有限数のＲＶから成る。これは、シーケンスの離散的及び有限の性質の直接的な結果である。全体の時間（Ｎ＝１２個のイベントの経時的シーケンス）を見ることが唯一可能である時計を考える。あなたが１時（時刻ｔ１）及び２時（時刻ｔ２）を見る場合、ｔ１からｔ２まで１時間が経過したこともあり得るし、１３時間が経過した又は時刻ｔ２でのイベントがｔ１の１１時間前に起こったこともあり得る。他の知識なしで、ある時間と、ある時間に対して１２の整数倍数をプラス又はマイナスした時間との間を識別することは可能ではない。この不確定性は、「エイリアシング」に直接関係する（これは、斯様な時系列から連続信号を再構成しようとする場合に起こる）。しかし、測定されたシーケンスの最後の部分の偶然の変化によって、最初の部分に続かない変数両方の変動が生じる。これによってインピーダンスの推定におけるバイアスが生じる。しかしながら、これは、適用された時間―周波数分析（下記を参照）において除外される。
【００９６】
不確定性原理
「時間の不確定性（度）」は、一連の成分が時間的にどのくらい拡がるかを表す時系列の特性である。離散的、有限及び巡回の時系列に対して、これは、図３Ａに示されるような図を与える。時系列の成分は、時間の円の周りに規則的に配置される質点により表される（各質点は、成分の二乗の絶対値に等しい）。不確定度Δｔは、基準ポイント（ここでは、ｔ＝０のポイント、式１９）に関して定められ、重心に密接に関係する（式Ｂ３）。時系列が一つだけゼロ以外の成分（時刻ｔ＝０で）を持つ場合、このとき、重心は、ｔ＝０に位置され、Δｔはゼロである。全ての成分が等しい重さを持つ場合、このとき、重心は中心に位置され、Δｔは高い（時間領域の最大の拡がりと関連している）。
【００９７】
同様に、「周波数の不確定性（度）」Δｆは、周波数領域での拡がりを表す（式２２）。離散的フーリエ変換を通じて、Ｎ個の成分を持つ各時系列は、所与の振幅及び周波数（離散的周波数領域の表現）を持つＮ個の調和振動の和として書ける。離散的時系列に対して、これらの周波数も、循環的である。周波数領域の拡がりは、図３Ａのものと比較できるような円の周りにＮ個の質点として二乗振幅を配置することにより視覚化できる（各質点は特定の周波数に対応する）。
【００９８】
連続的な時間及び周波数に対して、以下の式の不確定性関係により表される時間及び周波数の不確定度のトレードオフがあることが示された。
（５４）
【００９９】
これは、（特定の条件の下）連続的で無限の時間に対して成り立つ。Δｔが小さい場合、Δｆは、積が少なくとも１／４πに等しくなければならないので大きくなければならず、またこの逆も成立する。この不等式は、粒子の位置及び運動量に対するハイゼンベルクの有名な不確定性原理に直接関係する。しかしながら、離散的時系列に対して、積ΔｔΔｆの下の限界は、ゼロである（下記の「不確定性原理の導出」を参照）。しかし、Δｔ及びΔｆの連合値の制限がある。上記の「方法」セクション（式２５）において、離散的で、有限及び循環の時系列に対して、以下の式が示される。
（５５）
ここで、σｍｉｎはＮの関数である。よってΔｔ及びΔｆ両方が非常に小さくなることができない。和Δｔ２＋Ｎ２Δｆ２は、離散的フーリエ変換が最初の時系列（図４Ａ）に正確に等しい時系列に対して最小である。Ｎが大きくなるにつれて、時間及び周波数の最小の連合の不確定度を持つこの時系列は、時間のガウス関数に近づく。この時系列は、量子力学における２つの引きつけられている粒子の最小の不確定性状態（いわゆる「量子調和発振器」）に類似している。フォーブズ等（参照文献４を参照）は、斯様な発振器に対する波動式から、同様の時系列を導出した。しかしながら、今の場合、この最小の不確定性状態は、離散的及び有限時系列に対するΔｔ及びΔｆの定義から簡単に導出される。下記の「不確定性原理の導出」において、式５５の不確定性原理が大きなＮに対するガボールの不確定性原理と、どのくらい関係するかが示される（図１４を参照）。
【０１００】
離散的時間―周波数領域
提示された時間―周波数変換（ＴＦＴ）は、離散的及び有限時系列を、Ｍ＜ＮであるＭ個の成分の短く続く振動のセットへ区切る（時間―周波数領域への変換）。各振動は、式５５の原理に従う両方の領域の最小の不確定度を持って、特定の時間及び周波数の周りに集まる。ＴＦＴは、ウィンドウ化された「短期フーリエ変換」又は「実行しているフーリエ変換」となる。結果が第１の振動が始まる時間から独立しているように、短く続く振動は、時間的に互いに最大限重なり合う。時系列の始まりの前に始まり、時系列の終わりの後に終わる振動だけが、循環仮定に依存している（これらは、バイアスを回避するために、分析において除外される）。以下の「時間―周波数変換の導出」セクションにおいて、変換が「正規直交である」（式３３）ことが示され、これは、離散的時間―周波数領域への変換後時系列の全ての情報が保存されることを意味する。何も得られていないし、何も失われていない。逆変換は、元の時系列を回復する。
【０１０１】
正確には、離散的時間―周波数領域と呼ぶことはできず、しかし、選ばれた値Ｍに対する時間―周波数領域だけについてである。この整数は各短期振動の継続時間を決定するが、これらの振動が集まる異なる周波数の数に等しい（ＴＦＴ共振周波数）。以下の「不確定性原理の導出」で示されるように、式５５の不確定性原理による最適のΔｔ及びΔｆはＭに依存する（式Ｂ８においてＮをＭにより置換する、図１１も参照）。もっともなことだが、Ｍが大きいと、対応する最適なΔｔが大きくなり、対応するΔｆが小さくなる。現在の分析では、拡声器により誘発される異なる振動を識別可能であるのにΔｆが十分小さいように、Ｍは選択される。
【０１０２】
拡声器への入力
例示的な実施例では、我々は拡声器６への入力として５つの調和振動のセットを使用した（それぞれ８、１２、１６、２０及び２４Ｈｚの周波数）。分析のために、Ｍの値（実際の有効値、式Ｂ６参照）は、不確定度が時間ユニットでΔｔ＝０．０６２５ｓ及びΔｆ＝１．２７Ｈｚで表されるように、選択された。ＴＦＴ係数は、これら５つの周波数に対してのみ計算され、これは５つの共振周波数を持つ帯域通過フィルタリングになる。１つの従来技術のアプローチは、拡声器入力のために単一の周波数を使用する。これは、時間―周波数分析のより少ない要件をもたらす（これらは、矩形の窓を持つ帯域通過フィルタを使用し、本発明のものより所与のΔｔに対してより大きいΔｆを持つ）。単一の周波数に対して、周波数解像度は、患者の呼吸のより高い高調波を除去するために依然有効であるが、あまり重要でない。他方では、インピーダンスの周波数依存は、呼吸メカニクス上に付加的情報を供給できる。代わりに、広帯域のノイズが適用でき、これは、検討された周波数で比較的小さなパワーを持って、入力信号のパワーが全ての周波数上に徐々に分散されるという不利な点を持つ（拡声器の同じ全体出力で、より低いＳＮ比）。更に他のアプローチでは、調和周波数での振動の間の相互作用を防止するために、換気された患者に非調和シヌソイドでできたベンチレータ波形が適用された。しかしながら、斯様なアプローチは、自由に呼吸する被験者に適用できない。調和周波数での呼吸器系の幾つかの干渉は、これは広帯域の刺激の後、系の周波数反応で観察されなかったが、実際に本発明で生じてもよい。
【０１０３】
時間―周波数領域の分散分析
ＴＦＴの正規直交特性は、確率過程の分散を時間―周波数依存の成分に区切ることを可能にし、これは時間―周波数パワースペクトルを与える（下記の「時間―周波数領域の分散の分析」を参照）。ＴＦＴサンプルパワーは、総サンプル分散に対して所与の時間及び周波数についての短く続く振動の寄与を表す。よって、ＴＦＴパワースペクトルの分析は、時間―周波数領域の分散の分析になる（式４２）。モデル（式４３）において、一緒に二変量の巡回静止確率過程を形成する、流量及び圧力のノイズの独立源があると仮定される。循環仮定に依存する値が放棄されるとき、時間―周波数領域の多くのＮＳ＝５００の連続した値が、このプロセスからのサンプルとしてみなされる（８００Ｈｚのサンプルレートで０．６２５ｓのエピソード）。各時点（1秒につき８００回）で、このタイプの新規な確率過程が以前の分散と必ずしも等しいわけではない分散で開始すると仮定された。よって、測定は５００の値の最大に重なり合うエピソードに再分割され、それぞれは異なる静止確率過程からのサンプルとしてみなされた。ノイズの源は「白色（ホワイト）である」とみなされなかったが、パワーはＴＦＴフィルタの各周波数バンドで一定であると仮定された。このとき、５００の値の各エピソードの間の平均パワーは、等価の自由度が導出された近似のカイ二乗分布に従う。
【０１０４】
時間―周波数領域の二変量の最小二乗法
ＮＳ＝５００の連続した値の各エピソードに対して、拡声器入力と流量（ｂＱ、ｔ、ｍ）との間の線形関係及び入力と圧力（ｂＰ、ｔ、ｍ）との間の関係の係数は、時間―周波数領域の単純な線形回帰により推定された（式４４、図５）。これは、時刻ｔ及び周波数ｆｍに関連したバイアスがされていない推定値に結果としてなる（それぞれ、
及び
式４６を参照）。各推定値は、時間―周波数領域のＮＳ値上の時間平均化されたパワー及びクロススペクトルに基づく。その後、インピーダンスは、比率
により推定された。
【０１０５】
Ｎの選択は、系が静止しているとみなされる時間に依存する。インピーダンスが各吸気及び呼気の中でほぼ安定であると仮定すると、呼吸のために、〜０．５ｓのエピソードがもっともらしいようである。示された例において、これは、よく機能している（図７及び図８）。しかしながら、ほとんどの従来技術において、スペクトルは、より長いエピソード、通常１０ｓより長いエピソードにわたって平均化されている。幾つかの平均化されたスペクトル値は、１６ｓの全体のエピソードにわたって、〜０．６５ｓの重なっていないセグメントから導出された。しかしながら、重なっていないセグメントの使用は、ここで使用される最大に重なり合うセグメント（すなわち、ＴＦＴで使用されるＭ値のセグメント）と比較して、同じ長さのエピソードに対して非常に少ない自由度を与える。大きな課題は、１６ｓのエピソードの使用が、通常ＣＯＰＤのような疾患においてより発音される吸気と呼気との間のあり得る生理的違いを平均化するということである。その場合、１６ｓにわたる流量と圧力との間の低い干渉性は、ノイズの反映だけでなく、（疾患の一部である）インピーダンスの呼吸内可変性の反映でもある。想定された静止の期間は、あまり短く選択されるべきではない。図８の例では、信頼区間は吸気からの呼気への転換点で比較的広く、逆も成り立つ。これは、患者の自身の呼吸の高周波成分に起因する。これらの変化は０．５ｓより明らかに短く、その結果、これらは推定のノイズとして現れる。従って、想定された静止の期間は、これらの変化をノイズとして「解釈する」のに十分長くなければならない。
【０１０６】
インピーダンス推定値
が２つの正常に分散されたＲＶ（
及び
）の比率として導出されるので、これは、期待値を持たないコーシー分布に従う。それで、厳密に言うと、
は、期待値と本当の値との間の差という意味において、バイアスを持たない。まだ、
及び
がバイアスされていない正常に分散されたＲＶであるので、
は、近似により期待値ｚｔ、ｍ＝ｂＰ，ｔ，ｍ／ｂＱ，ｔ，ｍで通常分散される。Ｄａｒｏｃｚｙ及びＨａｎｔｏｓ（参照文献２）は、拡声器入力から流量及び圧力それぞれまでの回帰係数の比率としてインピーダンス推定値を導出する類似のアプローチに従っている。患者及び装置の機械的特性の単純なモデルに基づいて、彼らは、ノイズが流量又は圧力の何れかだけで生じると仮定される場合、系統的誤差が導入されると主張した。斯様なバイアスに対する実験証拠は見つかった。しかしながら、他者は、拡声器が圧力波及び主要なノイズを生成し、（圧力のノイズがゼロである一方）患者自身の呼吸の高周波成分が、純粋な流量源であると仮定した。ＴＦＴが用いられた場合、推定値は式５２で
であるだろう。他方では、他の文献は、式５２の
と類似した推定値を使用した。
及び
は、圧力又は流量のノイズがゼロであるという仮定の下、最小二乗推定値であるので、これは、現実と一致していない場合バイアスを導き、このことは、数人の著者によりすでに指摘されていた（参照文献２を参照）。流量と圧力との間の干渉性が高くて、さもなければ推定値を放棄するよう勧められる場合、このタイプのエラーは最小であると議論された。これは、干渉性が統一の傾向がある場合、
及び
両方が
（よって、本当のｚｔ，ｍ）に近づくことが分かる式５３に従っている。式５０及び５１も、流量のノイズが、インピーダンスに対する信頼領域に影響を及ぼすことを示す（特に、歪度パラメータβ、図６Ｂを参照）。１つの以前の文献では、帯域通過フィルタリングされた圧力は、時間の関数（上記参照）として、インピーダンスを得るために、流れにより割られた。結果として生じるインピーダンスは、流量及び圧力のサンプル分散の等価な寄与を持つ「全体の最小二乗」推定値に、数値的に等しい。これらの値が、しばしば、（入力から流量までの干渉性が入力から圧力までの干渉性と等しい場合）二変量の最小二乗法で得られたものと同じ範囲にあることが期待できる。しかしながら、この「推定値」は、ランダム誤差に非常に影響されやすい平均値から導出されない。
【０１０７】
よってこれまで、推定されたインピーダンス上の信頼限界は、相対的に時間が長い間隔でのみ導出されてきた。本発明の実施例の時間―周波数に依存する信頼限界は、短く続くシーケンスが静止確率過程からのサンプルであるという仮説の有効性の継続的印象を提供する。これらは時間の経過でインピーダンスの変化の重要性を試験し、（図７の吸い込む間のような）信頼できない推定を自動的に拒絶するための根拠を提供することを可能にし、これはリアルタイムに監視する際の実際的な利点である。
【０１０８】
図１０は、更に詳細に本発明による装置２（及び特にコンピュータ１６により実施される処理ステップ）により実施される方法を例示するフローチャートである。図９のように、ステップ１２１で、圧力波は被験者により用いられているマウスピース４に接続される拡声器６を使用して生成され、マウスピース４を通る空気の流量及び圧力の振動が測定されて、測定のそれぞれの時系列を与えるためにデジタル化される（ステップ１２３）。このように、流量及び圧力は、離散的時間ｔ、ｑｔ及びｐｔの関数として与えられる。
【０１０９】
その後、ステップ１２５で、時系列は、離散的時間ｔ及び異なる周波数ｆｍ＝ｍ／Ｍの関数として、流量
及び圧力
を与えるために離散的時間―周波数領域へ各々変換される。ここで、ｍは周波数指標であり、Ｍは時間―周波数フィルタの有効幅である。周波数は、課されたＦＯＴ周波数の周波数に合うように、選択される。
【０１１０】
ステップ１２５は、流量時系列ｑｔに対する以下の式を評価するステップを有する。
（５６）
ここで、ｉ２＝−１であり、加重因子ｗｌは時間のガウス関数
により近似され、Δｔは時間における選ばれた不確定度である。フィルタは、Ｋ１＝５・Δｔで丸められる。加重因子ｗｌがガウス関数により近似されるとここで説明される一方、三角窓、区分的線形近似又は多項式関数のような他の窓又は加重因子が本発明により考察されることも理解されるべきである。
【０１１１】
時間及び周波数の関数として圧力
は、以下の式
（５６ａ）
から同様に導出される。
【０１１２】
その後、ステップ１２７において、パワー及びクロススペクトルが時間及び周波数の関数として導出される。流量のパワーは、以下の式により与えられる。
（５７）
ここで、ＮＳは時間―周波数領域のサンプルの数であり、ＮＳが奇数の場合Ｋ２＝１／２（ＮＳ−１）であり、｜｜は絶対値を表わす。
【０１１３】
ＰＰ，ｔ，ｍにより表わされる圧力のパワーは、以下の式から同様に導出される。
（５７ａ）
【０１１４】
時間―周波数領域の拡声器入力
から流量
までのクロススペクトルは、以下の式により与えられる。
（５８）
ここで、＊は複素共役を示す。
【０１１５】
拡声器入力から圧力ＣｘＰ，ｔ，ｍまでのクロススペクトルは、同様に
及び
から以下の式により導出される。
（５８ａ）
【０１１６】
ステップ１２９では、パワー及びクロススペクトルは、拡声器入力から流量及び圧力それぞれへの伝達関数を決定するために用いられる。
【０１１７】
特に、拡声器入力から流量への伝達関数は、クロススペクトルとパワーとの比により以下の式により与えられる。
（５９）
拡声器入力から圧力への伝達関数ＢｘＰ，ｔ，ｍは同様に以下の式により導出される。
（５９ａ）
【０１１８】
よって、ステップ１３１で、呼吸器系のインピーダンスは、実数成分ＲｒＳｔ，ｍ及び虚数成分ＸｒＳｔ，ｍを持つ伝達関数ブロックの比率から決定できる。
（６０）
【０１１９】
更にまた、呼吸インピーダンス上の信頼限界は、以下の通りに導出できる（ステップ１３３）。
【０１２０】
パワー及びクロススペクトルの自由度ηの等価数は、以下の式により与えられる。
（６１）
ここでＰはＮｘＮ行列であり、Ｈはエルミート転置行列を表わし、ｔｒ｛｝は行列のトレースを表す。行列Ｐは、以下の式のように規定される。
Ｐ＝Ｓ（Ｍｍ＋ＭＭ−ｍ）／ＮＳ （６２）
ここで、Ｓは主対角線上にＮＳ個の１（及び残りはゼロ）を持つ対角行列Ｓである。すなわち、
Ｓ＝ｄｉａｇ｛０，．．．，０，１，．．．，１，０，．．．，０｝ （６３）
【０１２１】
（シフトされる）時間―周波数変換行列Ｍｍは巡回行列であり、第１の行のエントリは、ｌ＝−Ｋ３，．．．，Ｋ３及びＫ３＝１／２（Ｎ−１）に対してｗｌｅ２πｉｍｌ／Ｍであり、ここで、Ｎは奇数である。Ｆαは、自由度２、η−２でのＦ分布の上位１００（１−α）％ポイントとしてラベル付けられる。
【０１２２】
と
との間の二乗コヒーレンスは、以下の式により与えられる。
（６４）
と
との間の二乗コヒーレンスＫｘＰ，ｔ，ｍ２は同様に以下の式から導出される。
（６４ａ）
【０１２３】
流量に関係する変数Ａｑ，ｔ，ｍ２は、以下の式として規定される。
（６５）
Ａｐ，ｔ，ｍ２で示される圧力に関係する類似の変数は、同様に以下の式で導出される。
（６５ａ）
【０１２４】
従って、呼吸インピーダンスの実数部分及び虚数部分に対する１００（１−α）％信頼は、それぞれ、Ｒｒｓ・ｃ１±ｃ２及びＸｒｓ・ｃ１±ｃ２により与えられ、ｃ１＝１＋Ａｑ，ｔ，ｍ・ｃ３、ｃ２＝｜Ｚｒｓｔ，ｍ｜・ｃ３、及び
である。
【０１２５】
所与の時間及び周波数に対する呼吸インピーダンスの推定値は、Ａｑ，ｔ，ｍ≧第１の閾値、又はＡｐ，ｔ，ｍ≧第２の閾値の何れかである場合、重要でないとして拒絶される。ある具体例では、第１の閾値及び第２の閾値は１に設定されるので、所与の時間及び周波数に対する呼吸インピーダンスの推定値は、Ａｑ，ｔ，ｍ≧１、又はＡｐ，ｔ，ｍ≧１の何れかである場合、重要でないとして拒絶される。
【０１２６】
最初と最後の効果を取扱うために、アルゴリズムは、幾つかの巡回データバッファを含む。
【０１２７】
測定装置
測定エアフローは、ニューモタコヘッド８及び内蔵型圧力トランスジューサ１０（例えばJaeger Masterscreen pneumotach type BF/IEC601-1, Hoechberg, Germany）で測定できる。マウスピース４での圧力は、差圧トランスジューサ１２（例えばHans Rudolph Pneumotach amplifier 1 series 1110, Shawnee, KS）で周囲空気を基準にして測定できる。圧力振動は、アンプ１８（例えばHarman Kardon HK970 amplifier, Washington DC）により増幅されるＡＤ変換カード１４（例えばNational Instruments PCI-6221, Dallas,TX）を通じたパーソナルコンピュータ１６（例えばHewlett Packard Compac dc 7600, PaloAlto, CA）からのアナログ出力信号により制御される拡声器６（例えばJaeger Masterscreen IOS, Hoechberg, Germany）によりＦＯＴ装置２内で生成できる。被験者は、ニューモタコヘッド８に接続されたメッシュワイヤ抵抗９を通じて呼吸できる。アナログ入力信号は、８００Ｈｚのサンプルレートで同じＡＤ変換カード１４を通じてデジタルシーケンスに変換でき、パーソナルコンピュータ１６に保存できる。出力信号の生成及びデータの分析は、適切なコンピュータソフトウェアを実行しているコンピュータ又はプロセッサ１６により実施できる。
【０１２８】
測定は、多くの呼吸サイクル、例えば９０秒の静かな呼吸をカバーする期間の間で、装置２を使用してなされる。
【０１２９】
不確定性原理の導出
時間及び周波数の不確定度（性）の幾何学的な解釈
式２０の規定によると、不確定度Δｔは、時系列ｘ≡｛ｘｔ｝の「重心」に直接関する。時系列のイベントは、値｜ｘｔ｜２を持つ質点として表わされ、原点を中心に、複素平面内の半径ｒｔを持つ円の周辺部に定間隔に配置される。各イベントは、複素平面内のポイントｒｔωｔで起こり、ここで、
及び
である。加重平均値
は、以下の式のように規定できる。
（Ｂ１）
【０１３０】
重心は
に位置される。行列式では、以下の式となる。
（Ｂ２）
ここで、Ω≡ｄｉａｇ｛１，ω，．．．，ωＮー１｝である。Ωがユニタリ行列であるという事実を使用して、二乗不確定度Δｔ２は、このとき、以下のように低減できる。
（Ｂ３）
【０１３１】
これは、図１２の基準点Ｒを通る軸（円の平面と直角をなす）についての第２の慣性モーメントに対する「平行軸定理」に従う。図１２は、時間の不確定度ΔｔとＮ＝１６に対するＣＨＣの固有ベクトルの重心との間の関係を例示する。固有ベクトルは、原点を中心として、半径ｒｔを持つ複素平面内の時間円に沿ってマップされる。各固有ベクトルの重心は、黒いドットにより表される。Ｃの対応する特異値は（基準点Ｒに近い）右側で重心に対して最小であり、図の左へ行くと徐々に増大する。各固有ベクトルに対する不確定度Δｔは、Ｒから重心を通る垂直線と円との交差点の１つまでの傾斜距離である。
【０１３２】
図１３の表現では、式１１から
であることが分かり、ここで、
は図１２のＡからＲまでの距離である。ピタゴラスの定理を使用して、
から、
であることが分かる。
【０１３３】
周波数領域において、不確定度Δｆは同様のやり方で幾何学的に解釈できる。
【０１３４】
行列Ｃの特異値
式２５の不確定性原理は、Ｃの二乗最小特異値及び二乗最大特異値であるσｍｉｎ２及びσｍａｘ２により決定される。特異値σは、
ｄｅｔ（ＣＨＣ−σ２Ｉ）＝０（Ｂ４）
から導出できる。ここで、ｄｅｔ（）は行列の決定要素を表わす。Ｎ＝２に対して、σｍｉｎ２＝２（２−√２）／π２及びσｍａｘ２＝２（２＋√２）／π２であることが容易に分かる。より高いＮに対して、Ｃの特異値（ＣＨＣの固有値）を導出するための様々な戦略が開発されてきた。図１１は、Ｎの関数として、σｍｉｎ及びσｍａｘを示す。より大きなＮ（例えば、Ｎ＞１５）に対して、σｍｉｎは
に近づき、σｍａｘはＮ√２／πに近づくことが分かる。時間円の半径に関して、このことは、
及び
を意味する。
【０１３５】
ｖを特異値σに対応するＣＨＣの単位固有ベクトルとする。正規方程式（ＣＨＣ−σ２Ｉ）ｖ＝０は、以下のように書ける。
（Ｂ５）
【０１３６】
第１の項は、実軸に沿った基準点Ｒからのｖの成分の出発に依存する（式Ｂ３を参照）。第２の項は、（離散的及び巡回の時間に対する）ｖの成分の２次の導関数とみることができる。式Ｂ５は、量子調和振動に対するシュレディンガー方程式に相当する。解は、離散的正規直交マシュー（Mathieu）関数と考えられる。σｍｉｎに対応する固有ベクトルｖｍｉｎは、時間の単一モードの関数である（図４Ａのような）。大きなＮに対して、ｖｍｉｎはベクトルｖ’ｍｉｎにより近似でき、その成分ν’ｍｉｎ，ｔは以下の式のような時間のガウス関数である。
（Ｂ６）
ここで、Ｎは、ｔ≧（１／２）Ｎの場合、ｔから減算されなければならない。ｔ＞５・Δｔの場合、成分ν’ｍｉｎ，ｔはゼロに近いので、ガウス関数は、最大値付近の１０・Δｔ個の最も大きい値だけを使用して、残りをゼロに設定することにより丸められる。結果として近似された固有ベクトルｖ’及び差ベクトルｄ≡ｖ−ｖ’をコールして、ノルム
は、近似値に作られる誤差の尺度として採ることができる。適切な数学的ソフトウェアを使用して、Ｎ＞７５の場合、相対的な誤差
が０．００５未満であることが分かる。
【０１３７】
ＣＨＣがＦと代わるので、これらの行列は、同じ固有ベクトルを共有する。Ｆ４＝Ｉであるので、Ｆの固有値は、１、−１、ｉ及び−ｉである。結果として、各ユニット固有ベクトルｖに対して、時間の不確定度は、周波数の不確定度に直接関係する。
（Ｂ７）
【０１３８】
これは、各固有ベクトルｖに対して、全体の不確定度が時間及び周波数にわたって均一に分散されることを意味する。
（Ｂ８）
よって、各固有ベクトルｖに対して、Δｔ＝ＮΔｆ＝σ／√２である。式Ｂ６によるｖｍｉｎのガウス近似は、以下のように書き直せる。
従って、全体の不確定度σｍｉｎ２は、このガウス関数の分散に等しい。式Ｂ８は、また、Ｎ＞１５に対して、最小及び最大の不確定度に対応するΔｔの値がそれぞれ、
及び
であることも意味する。後者は、２ｒｔが時間円内の最大可能なΔｔであるので、直観的に合理的である（図１３を参照）。図１３は、Ｎ＝１６に対する時間の不確定度Δｔと周波数の不確定度Δｆとの間の関係を例示する。全ての可能性があるＮ―ベクトルｘに対するアクセス可能な領域は、半径σｍｉｎ及びσｍａｘ（Ｃの最小及び最大特異値）を持つ２つの円により囲まれている。図１３において、ｒｔは時間円の半径であり、黒点はＣＨＣの固有ベクトルに属し、ｖｍｉｎ及びｖｍａｘはσｍｉｎ及びσｍａｘに対応する固有ベクトルであり、ｅ０及びｅ８はｔ＝０及びｔ＝８に対応する正準ベクトルであり、ｆ０及びｆ８は周波数ｆｎ＝０及びｆｎ＝８／１６＝１／２を持つ調和振動である。
【０１３９】
上述されたように、図１２は、Ｎ＝１６に対する時間円の全ての固有ベクトルｖに対する重心を示す。ＣＨＣは正規行列であるので、スペクトル定理に従ってＮ個の直交固有ベクトルを持つ。式Ｂ３からわかるように、対応するΔｔは、Ｒから複素平面の重心までの水平距離に依存する。Δｔは、図１２から、Ｒから重心を通る垂直線と円との交点までの弦の長さとして、読み取れる。Ｒに最も近い重心は、σｍｉｎと一致する。更に離れたところにある重心は、σ（及び
）の増大する値に対応する。比較的大きなＮに対して、ｖｍｉｎの重心とＲとの間の距離が式Ｂ３に従って１／４に近づき、そのΔｔ，ｍｉｎを見つけることは
に近づく。同じ固有値σ２＝Ｎ２／π２を持つ２つの直交固有ベクトルがあるという意味で、ＣＨＣの多数の固有値がある。Ｎが４により割り切れる場合のみ、斯様な多数の固有値が発生するだろう。これら２つの直交固有ベクトルの重心は両方とも虚数軸上に位置する。これらは、同じΔｔを持ち、
に等しい。対応するΔｆは
に等しい。これらの２つの固有ベクトルは、時間及び周波数領域において最大限拡がった（重心が時間及び周波数円両方の中心に位置する場合、得られるだろう値と同一である）。σのより高い値に対応する固有ベクトルは、（Ｒに対する）より大きな「不確定度」と関係するが、時間的により小さな拡がりを持つ。最大の総不確定度σｍａｘが固有ベクトルｖｍａｘに対して到達され、Δｔが２ｒｔにほぼ等しくなるまで、ベクトルの重みは時間及び周波数両方の円の反対側に実際にますます集中される。図１２の明らかな対称性に留意されたい（重心は、虚数軸に反映される）。
【０１４０】
図１３はΔｔの関数としてＮΔｆを示す。式２５の不確定性原理のため、アクセス可能な領域は、半径σｍｉｎ及びσｍａｘで原点の周りに中心がある２つの円により境界づけられる。これは、これらの円の間の全ての値が可能であるというわけではないが、円の外側の全ての値が不可能であることを意味する。ＣＨＣの固有ベクトルに帰属する座標は、全て同一ライン上に位置する（式Ｂ７のため）。幾つかの極端な場合も示される。１つは、正規ベクトル
である。この「重み」
は、ｔ＝０に完全に集結されているので、Δｔはゼロである。周波数領域への変換は
（Ｂ９）
を与える。ここで、ｆｎは式９として規定され、符号「＊」は複素共役を表わし、１は１だけを含むＮベクトルである。よって、ｅ０の重みは周波数領域の全ての周波数上に均一に分散され、
又は
である。他の極端な場合は、ｅ８である。（正規ベクトルｅｔを｛０，．．．０，１，０，．．．，０｝として規定し、ここで、１はゼロから始めてｔ番目の位置に表わされる。）ｅ８の重みは、複素平面内のポイント（−ｒｔ，０）に完全に集中されているのでΔｔは最大であり、２ｒｔに等しい。周波数領域において、
（Ｂ１０）
【０１４１】
列ベクトルｆ８＊は周波数ｆｎ＝８／１６＝１／２を持つ調和振動を記述する。その二乗成分も全ての周波数上に均一に分散されるので、
である。時間領域から始めると、他の極端な場合は、その重みが時間領域内に均等に分散され、ｆｎ＝０で周波数領域において鋭く集中する（
及びΔｆ＝０）振動ｆ０と、その重みが時間領域内に均等に分散され、ｆｎ＝１／２で周波数領域において鋭く集中する（
及びＮΔｆ＝２ｒｔ）ｆ８とである。Δｔ及びＮΔｆに対する可能性がある値が範囲［０，２ｒｔ］に限定されるので、全ての可能性がある組合せ（Δｔ，ＮΔｆ）がｖｍｉｎ，ｅ０，ｆ８，ｖｍａｘ，ｅ８及びｆ０に対する座標により囲まれる図１３の領域に限定されることが期待できる。これらの極端なベクトルは、通常の規則の実例でもある：（Ｆによる前置乗算を通じて）時系列ベクトルｘのフーリエ変換をすることは、図１３のプロット線の（Δｔ，ＮΔｆ）座標を交換する。別の表現では：ｘ及びＦｘに対する座標が同一のラインに反映される。これは、式２１及び式２３の規定を使用して、容易に検証される。
【０１４２】
ガボールの不確定性原理との比較
１９４６年のガボールにより説明された不確定性原理（下記の参照文献５を参照）は、連続時間及び周波数に対する積ΔｔΔｆに対する下位限界を決定する。
（Ｂ１１）
【０１４３】
この不等式は、式２５による離散的時間と周波数に対する不確定性原理にどのように関係するだろうか。式でＢ１１では、時間が無限のライン（又は、無限の半径を持つ円）上にマップされているが、Δｔ２は、また、ｔ＝０であるポイントに対して垂直な軸の周りの第２の慣性モーメントと規定される。式Ｂ１１の導出において、平均時間がゼロである（重心がｔ＝０に位置される）と仮定される。同じことは、周波数領域の拡がりに対しても満たされる（周波数が離散的でも周期的でもなく、連続的で無限であるという事実に関係なく）。
【０１４４】
式Ｂ１１と式２５との不等式の間の主要な違いは、ガボールの不確定性原理によると、ΔｔΔｆの下位限界ｔがΔｔの任意のゼロでない値を持つ（式Ｂ６と同じ形式の）ｔのガウス関数に対して達成される（参照文献５を参照）という事実にある。これは、一組の線形の独立（無限の）ベクトルである。しかしながら、式２５による総不確定度Δｔ２＋ＮΔｆ２の下位限界は、一つのΔｔの値（
に等しい）に対応するＣＨＣのたった一つの固有空間（ｖｍｉｎにわたる一次元複素部分空間）に対して達成される。他方では、この固有空間のベクトルは、ガボールの不確定性原理の下位限界を達成する（Ｎ→∞の限界状況において）。ＣＨＣの各固有ベクトルに対して、式Ｂ８から以下の式が成り立つ。
（Ｂ１２）
【０１４５】
ｖｍｉｎ（又は対応する固有空間の他の任意のベクトル）に対して、σは、Ｎが大きくなるにつれて、
に近づく（図１１）。このとき数の近似は、Ｎ→∞だと、以下の式を示す。
（Ｂ１３）
【０１４６】
しかしながら、有限のＮに対して、１／（４π）はΔｔΔｆに対する絶対の下位限界ではない。ｅ０及びｆ０に対して積はゼロであり、小さな乱数の追加は、わずかに異なるベクトルに対してゼロに近いことを示す。
【０１４７】
（Ｍイベントからなる、１＜Ｍ≦Ｎ）より短い時系列に対して得られたｖｍｉｎベクトルから導出されるＮ個のベクトルを考える場合、２つの不確定性原理の関係はより明らかになる。ｖｍｉｎが８ｘ８行列ＣＨＣから導出される図１４の例を参照されたい。Ｆｖｍｉｎの成分は、周波数ｆｎ（黒いドット）の関数として示される。これらの周波数は、１／８の倍数である。（ベクトルの単一モードの構造が円形の時間ベース上に完全なままであるように）８ベクトルｖｍｉｎの２つの最も小さな成分間に８個のゼロのブロックを挿入し、これは結果的に１６−ベクトルとなる。
（Ｂ１４）
【０１４８】
これは、Δｔの僅かな増大（０．７５８３から０．７９００へ）とほとんど重要でないΔｆの増大（０．０９４８から０．０９４９へ）を与える。図１４の白丸は、８個のゼロのブロックが時間領域の固有ベクトルの２つの最も小さな値の間に挿入されるとき周波数領域の補間された値を表し、実線は、時間領域の加えられたゼロの数が無限に向かうにつれて補間された値を表す。結果として生じる１６−ベクトルは、１／１６の倍数である周波数（白丸）で、周波数領域の補間を導く。これは、（時間領域の「ゼロ詰め」を通じた）周波数領域の良く知られた補間の形式である。
【０１４９】
このやり方で、所与のＮに対して、Ｎ−１の多くのベクトルは、ＭｘＭ行列ＣＨＣから得られるより小さなＭ次元ｖｍｉｎベクトルにゼロを加えることにより導出できる。自明の例は、Ｍ＝１に対して「１―ベクトル」ｖｍｉｎ＝１にゼロを加えることにより導出されるｅ０である。図１５は、Ｎ＝１６に対するＮΔｆ対Δｔのプロット線の対応する不確定度を示す。これらの値は、Δｔ＝０、
及び
のラインにより区切られるゾーンＩに現れる。Ｍ＝１に対する値はｅ０と一致し、Ｍ＝２に対する値は矢印により示され、より高いＭに対する値は（Ｎ＝１６に対するｖｍｉｎに一致する）識別ラインに位置されるＭ＝Ｎ＝１６に対する値に徐々に近づく。比較的大きなＭに対して、Δｔ及びΔｆがゼロの付加によりほとんど変化しないのに対し関連するσは
に近づくので、以下の式となる。
（Ｂ１５）
【０１５０】
結果として、比較的大きなＭに対する値は、図１５の双曲線ｘ・ｙ＝Ｎ／（４π）より僅かに下にプロットされる。時間が進みＮが大きくなるにつれ、値のトレイルはこの双曲線に沿って成長し、新しい値が識別ラインに現れる。ΔｆがΔｔに対してプロットされる場合、関係のある双曲線はｘ・ｙ＝１／（４π）により与えられ、ガボールの不確定性原理の下位限界で異なるガウス分布に対応する全ての可能性がある値を含む。
【０１５１】
正式の数学的証明なしで、Ｍ＜Ｎに対するｖｍｉｎから導出されるＮ−ベクトルは、図１５の（Δｔ，ＮΔｆ）平面のアクセス可能な領域までの更なる限界を持つ。これらのベクトルがＭがゼロでないイベントのシーケンスに対して最小の可能性がある総不確定度を持つので、図１５の黒いドットによりマークされる領域の外側の値が発生しないと期待できる。これは、時間領域の重みの変更なしで、周波数領域のベクトルを巡回的にシフトさせる（周波数円の反対側の部分）、Ω８とゾーンＩのベクトルを前置乗算することにより導出される、ゾーンＩＩの値に対しても成り立つ。よって、これらのベクトルのΔｆは、Δｔの変化なしで最大化される。これとは逆に、ゾーンＶＩのポイントは、時間領域のベクトルを巡回的にシフトさせる（これによりΔｆの変化なしでΔｔを最大化する）、Ｔ８と前置乗算された、ゾーンＩのベクトルから得られる。ゾーンＶのベクトルは、Ｔ８と前置乗算を通じてゾーンＩＩのベクトルから生じる。ゾーンＶＩＩＩのベクトルは、フーリエ変換をする（Ｆによる前置乗算は識別ラインの反映を意味する）ことにより、ゾーンＩのベクトルから生じる。Ω８によるこれらのベクトルの前置乗算は、ゾーンＩＩＩのベクトルを与え、Ｔ８による前置乗算は、ゾーンＶＩＩのベクトルを与え、Ω８による後者のベクトルの前置乗算は、ゾーンＩＶのベクトルを与える。（一つの時間周波数ゾーンから他の時間周波数ゾーンへ行かせる多くのやり方があることは明らかであり、例えば、ゾーンＶＩのベクトルのフーリエ変換を採用することはそれらをゾーンＩＩＩへ移動させる）。結論として、Δｔ、ＮΔｆの可能性がある組合せが、図１５の黒いドット（及び白丸）によりアウトラインが引かれたスペード状の領域に限定されることが期待できる。
【０１５２】
ゼロではない平均時間及び周波数に対する不確定性原理
式１９の基準時間ｔＲが、ゼロではないことが可能な場合、そのとき、以下の式が成り立つ。
（Ｂ１６）
【０１５３】
同様に、ゼロではない基準周波数
に対して、以下の式が成り立つ。
（Ｂ１７）
【０１５４】
及び
とする。
【０１５５】
ＣがｔＲ＝０及びｆＲ＝０に関して依然規定されるとき、ＣＲＨＣＲ及びＣＨＣは同様に以下の式になる。
（Ｂ１８）
【０１５６】
これは、任意の整数ｋ、ｎに対して有効であるＴｋΩｎ＝ω−ｎｋΩｎＴｋという関係を使用して容易に示される。それで、（σ２，ｖ）がＣＨＣの固有値−固有ベクトル対である場合、このとき、以下の式が成り立つ。
（Ｂ１９）
これは、
が固有値σ２を持つＣＲＨＣＲの固有ベクトルであることを意味する。時間領域のｔＲ及び周波数領域のｆＲに対する
の不確定度は、両方の領域のゼロに対するｖの不確定度に等しいことが成り立つ。これは、思いがけないことではない。（
による）周波数領域の巡回シフト及び（
による）時間領域の巡回シフトの後でも、
の成分の二乗絶対値（「重み」）のｖのものと同一である。
【０１５７】
時間―周波数変換の導出
時間―周波数変換行列Ｍの列は、正規直交である。
ＭＨＭ＝Ｉ（Ｃ１）
【０１５８】
これは、以下のように証明できる。積が巡回行列の和として以下のように書ける。
（Ｃ２）
【０１５９】
Ｔ−１によるｈｍＨの後置乗算は、その成分が右へ巡回的にシフトされることを意味するので、巡回Ｍｍは以下のように書ける。
（Ｃ３）
【０１６０】
巡回は正規行列なので、これらは、エルミート転置行列と置換される。
（Ｃ４）
ここで、
はｈｍの「複素自己相関」である。
【０１６１】
その成分は以下の通りである。
（Ｃ５）
【０１６２】
Ｎ―ベクトルｈｍの規定（式２６及び式２７）によるとこれはΩｍＮ／Ｍｈ０としても書ける。ＴｋΩｎ＝ω−ｎｋΩｎＴｋを使用して、ゼロ以外の成分は、下記の式に等しいことが成り立つ。
（Ｃ６）
ここで、
である。これらの成分がｍにわたって加算されるとき、下記の式となる。
（Ｃ７）
この式は、｛ｅ２πｉｍｕ／Ｍ｝がｍ＝０，．．．，Ｍ−１に対して等比数列である事実から成り立つ。ｈｍが単位ベクトルであるとみなされるので、
であることに留意されたい。結果として、式Ｃ２の和は、以下の式に書き直せる。
（Ｃ８）
これは、式Ｃ１の証明を完成する。
【０１６３】
正規直交性特性（式Ｃ１）は、ｈ０が、不確定度Δｔ及びＭ’＝１０・Δｔでゼロ以外の成分を持つ、丸められ正規化されたガウス分布により近似される場合（式Ｂ４を使用して）、ｈ０の選択に依存しないので（
と仮定する）、依然成り立つ。
【０１６４】
時間―周波数合成
正規直交性特性の直接的な推論は、所与の時系列ベクトルｘが短く後続する振動の加重和として書けるということである（式３４）。式２８の規定を使用して、以下の式が成り立つ。
（Ｃ９）
【０１６５】
時間―周波数領域の分散の分析
ＴＦＴの正規直交性特性の他の直接的な推論は、Ｎ―ベクトルｘの「エネルギー」（二乗されたノルム）が変換後保存される。
（Ｄ１）
【０１６６】
これは、時間―周波数領域の「分散分析」（ＡＮＯＶＡ）を可能にする。上記は、また、ｘのエネルギーがＭ周波数依存の成分に区切られることを意味する。
（Ｄ２）
【０１６７】
ランダムなＮ―ベクトルＸが、期待値Ｅ｛Ｘ｝＝０及び分散ｖａｒ｛Ｘ｝＝σ２Ｉを持つ多変量の正規分布を持つと仮定する。このとき、二次形式
は、Ｎ自由度のカイ二乗分布に従う。式Ｄ２によると、
は、Ｍ周波数依存の成分
に区切られる。二次形式
は、周波数ｆｍに対するＸの「ＴＦＴサンプルパワースペクトル」と呼ばれる。ホワイトノイズに対して、このＲＶの期待値は、以下の通りである。
ｔｒ｛・｝は行列のトレースを表わす。式Ｃ４を使用すると、ｔｒ｛ＭｍＨＭｍ｝＝ＮｈｍＨｈｍ＝Ｎが成り立つので、以下の式が成り立つ。
（Ｄ３）
【０１６８】
しかしながら、上記規定によるサンプルパワースペクトルは、循環仮定のためバイアスに影響される。
の成分は、Ｍ個のゼロ以外の要素を持つ線形時間不変のフィルタとの循環畳みこみによりＸから導出される（式２６−式３１）。このフィルタの構成は、
の最初のＭ−Ｋ−１個及び最後のＫ個の成分が循環仮定に依存するようになされ（これらの整数はゼロより大きいとする）、ここで、ＫはＭ／２以下の最大の整数である。このように、循環仮定によるバイアスがないパワー推定値は、以下の式で規定される「選択行列」Ｓで得られる。
（Ｄ４）
【０１６９】
Ｓによる
の前置乗算は、循環仮定とは独立している
の成分を選択し、これは、バイアスされていない二次形式
を与える。相補的周波数ｆＭ−ｍも含まれるとき、「全」バイアスがないサンプルパワーは、以下の式で規定される。
（Ｄ５）
【０１７０】
ｔ―ｆ領域の「サンプル」の数は、ＮＳ＝Ｎ−Ｍ＋１である。
【０１７１】
サンプルパワー（σ２により割られた）は、等価の自由度（ＥＤＯＦ）を持つ近似のカイ二乗分布に従う。
（Ｄ６）
【０１７２】
及び
に対して、
とすると、ゼロ平均ホワイトノイズＸに対して以下の式が成り立つ。
（Ｄ７）
これは以下の式を与える。
（Ｄ８）
【０１７３】
ここでは、ηは数学的コンピュータソフトウェアを使用して、式Ｄ８によって数値的に計算される。
【０１７４】
信頼限界の導出
及び
の平均及び分散
流量
のノイズがゼロ平均ホワイトノイズであるとみなされるので、式４３を使用して、以下の式が成り立つ。
（Ｅ１）
【０１７５】
の分散はσＱ，ｍ２Ｉである。積
が決定論的であるので、以下の式が成り立つ。
（Ｅ２）
【０１７６】
の平均及び分散は、式４６から以下の式が成り立つ。
（Ｅ３）
（Ｅ４）
【０１７７】
結果として、
は、ｂＱ，ｔ，ｍのバイアスされていない推定値である。
に対する平均及び分散は、類似の態様で、以下の式となる。
（Ｅ５）
【０１７８】
信頼限界
及び
の信頼限界は、最小二乗分析の標準的アプローチからわかる。式４３及び式４４の組み合わせは、以下の式を与える。
（Ｅ６）
【０１７９】
と
との間の直交性のため（図５を参照）、以下の式が成り立つ。
（Ｅ７）
【０１８０】
ランダム変数
は、
の形式を持ち、よって、（１つの自由度で）カイ二乗変数として（正確に）分散される。相補的周波数ｆＭ−ｍ＝１−ｍ／Ｍでの
も推定に含まれるとき、結果として生じるカイ二乗変数は２つの自由度を持つ。ノイズが全体のバイアスされていないパワースペクトルにより推定されるとき（式Ｄ５）、式Ｅ７は、ノイズの全パワーをＥＤＯＦη＝２＋（η−２）を持つ直交成分に分解させる。その後、式Ｅ７の右辺の２つの直交成分の比率は、（２，ηー２）程度の自由度でＦ分布に従う。この分布の上位１００（１−α）％のポイントをＦαにより示す。このとき、
（Ｅ８）
【０１８１】
ＴＦＴサンプルコヒーレンス（式４５）の規定を使用し、ＫＱ，ｔ，ｍ２＞０とすると、これは以下の式に書ける。
（Ｅ９）
【０１８２】
ＡＱ，ｔ，ｍ２によりこの不等の右辺を書き直すと、
（Ｅ１０）
【０１８３】
よって、
に対する１００（１−α）％信頼領域は、中心
及び半径
を持つ複素平面内の円により記述される。
に対する信頼領域は、同様の態様で導出される。
（Ｅ１１）
【０１８４】
に対する信頼領域に到達するために、
により式１０の両辺を乗算し、
を置換すると、以下の式を与える。
（Ｅ１２）
【０１８５】
複素変数の実数部分及び虚数部分を使用して、これは以下の式に再編成できる。
（Ｅ１３）
【０１８６】
次に、
により式Ｅ１１の両辺を割る。
により
を置換し（式４７）、Ｕにより
を置換すると、
（Ｅ１４）
【０１８７】
に対する信頼領域が、ここでツーステップアプローチで導出される。式Ｅ１３によると、Ｕは、１００（１−α）％の機会で、中心
及び半径
を持つ円内にある（図１６Ａ）。
及び
が相関しないと仮定する。このとき、所与のＵの値に対して、ｚｔ，ｍは、１００（１−α）％の機会で、中心Ｕ及び半径｜Ｕ｜ＡＰ，ｔ，ｍを持つ第２の円内にある（式Ｅ１４）。｜Ｕ｜に対する２つの極値、図１６ＡのＵｍｉｎ、Ｕｍａｘの大きさを考える。｜ｚｔ，ｍ｜に対する対応する極値は、図１６Ｂのｚｍｉｎ及びｚｍａｘの大きさである。式Ｅ１３及び式Ｅ１４から、
（Ｅ１５）
【０１８８】
に対する１００（１−α）％信頼領域の控え目な限界が、ここでｚｍｉｎ及びｚｍａｘを含む円により与えられ、その中心が
及び原点を通るライン上にある。結果として、対応する信頼領域は、以下の式により記述される。
（Ｅ１６）
【０１８９】
この領域は、中心
及び半径
を持つ図１６Ｂ及び図１６Ｃの円により区切られる。この円は、
の周りに同心円ではない。所与の
のｚｔ，ｍに対する可能な値の分布は、
に対して非対称である（図１６Ｄ）。この非対称性は、図１６Ｃのｔａｎβにより表される。式Ｅ１６から、以下の式が成り立つ。
（Ｅ１７）
【０１９０】
特別なケース１：流れにノイズがない。流れ内のノイズがゼロである、
である場合、そのとき、
及び
である（式４３及び式４６）。ｚｔ，ｍの対応する推定値を
により示す。
と仮定すると
である。
【０１９１】
この場合、
なので、
（Ｅ１８）
【０１９２】
これは、「通常の最小二乗推定値」である。このとき、ＫＱｘ，ｔ，ｍ２＝１及びＡＱ，ｔ，ｍ＝０である。ＡＱ，ｔ，ｍ＝０であるが、信頼限界は、式Ｅ１６により依然記述される。ここで、半径
を持つ信頼円は、
を中心に持つ。非対称パラメータｔａｎβはゼロに等しい。
【０１９３】
特別なケース２：圧力のノイズがない。圧力のノイズがゼロである、
である場合、そのとき、
及び
である。ｚｔ，ｍの対応する推定値が
により示されるとき、
（Ｅ１９）
【０１９４】
これは、「データ最小二乗法推定値」である。ここで、ＫＰｘ，ｔ，ｍ２＝１及びＡＰ，ｔ，ｍ＝０である。式Ｅ１６から、信頼円は、中心
及び半径
を持つ。非対称パラメータｔａｎβはＡＱ，ｔ，ｍに等しい。２つの推定値
及び
は、流量と圧力との間の二乗コヒーレンスを通じて関係がある。
（Ｅ２０）
【０１９５】
例：時間及び周波数の関数としての二乗コヒーレンス。図１７は、図７と同じ記録のための二乗コヒーレンスを示す。この例では、拡声器入力と流れとの間のコヒーレンスは、入力と圧力との間のコヒーレンスより概して低かった。図から明らかなように、コヒーレンスは、時間及び周波数両方の関数として変化する。最も小さなコヒーレンスは、吸気から呼気への転換点及びその逆の転換点で、最も低い周波数で得られる。飲み込むことは、特に入力と流れとの間のコヒーレンスへの深い影響を持つ。図１８は、ＣＯＰＤを持つ患者に対して図８に示される推定されたインピーダンスに対する時間及び周波数の関数としての二乗コヒーレンスを示す。呼吸フェーズ間の転換点のコヒーレンスの低下は、図１７の通常の例より非常に多く発音される。図１７及び図１８両方において、太線は、異なる共振周波数（８〜２４Ｈｚ）に対する時間の関数として、拡声器入力と流れとの間の二乗コヒーレンスＫＱｘ，ｔ，ｍ２を表し、細い線は、拡声器入力と圧力との間の二乗コヒーレンスＫＰｘ，ｔ，ｍ２を表す。「流量」は、エアフローの低周波成分を表し、「ｉｎｓｐ」は吸気を表し、「ｅｘｐ」は呼気を表す。
【０１９６】
結論として、最適時間―周波数解像度を持つ被験者の呼吸器系の伝達関数ブロックを推定するための方法が提示される。
【０１９７】
上述のように、呼吸インピーダンスの推定は、気道の妨害を評価するか、又は疾患の発病度を推定するための推定を使用する診断ツールによりなされる。従って、診断ツールは、呼吸インピーダンスに影響を及ぼすべき（薬理学的な又は他の）治療の効果を評価するためにも用いられる。
【０１９８】
例えば、呼吸インピーダンスの推定は、（ｉ）慢性閉塞性肺疾患（ＣＯＰＤ）の呼気の流量限界；（ｉｉ）発病度自体が時間の経過により吸入された気道拡張器の効果（例えば交感神経作用薬又は副交感神経作用薬）を評価するために用いられる（これは、調査設定に、更に、患者が標準強制的呼吸操作を実施できない臨床設定に適用できる）ＣＯＰＤ又は喘息の気道閉塞の当該発病度；（ｉｉｉ）喘息の気道閉塞又は睡眠の間のＣＯＰＤ；又は、（ｉｖ）睡眠無呼吸―呼吸低下症候群の疑いがある患者の睡眠中の上位の気道閉塞を検出するために使用できる。
【０１９９】
呼吸インピーダンスの推定は、また又は代わりに、内科的疾患の治療で使用される機械の設定を適応させるために用いられる。例えば、呼吸インピーダンスの推定が、ＣＯＰＤの非侵襲性ベンチレーションの設定の調整において用いられる。呼吸インピーダンスは、連続的に気道インピーダンスの重症度に関する情報を提供するためにベンチレータへの入力として用いられる（信頼値により示される信頼できない値は廃棄される）。この情報は、また、呼気の流量限界がちょうど克服される二相性陽気道圧のレベルを調整するためにも用いられる。他の例として、呼吸インピーダンスの推定は、援助として閉塞性睡眠時無呼吸症候群を持つ患者の持続気道陽圧のレベルをガイドするための補助として用いられる。
【０２００】
本発明は、図面及び前述の説明で例示され詳述されてきたが、斯様な図例及び説明は、図示的及び例示的であって、限定的でないとみなされるべきであり、本発明は、開示された実施例に限定されない。
【０２０１】
開示された実施例のバリエーションは、図面、明細書及び添付の請求の範囲の検討から、請求された本発明を実施する当業者により理解され遂行できる。請求項において、用語「有する」は他の要素又はステップを除外しないし、不定冠詞「ａ」又は「ａｎ」は複数を除外しない。単一のプロセッサ又は他のユニットが、請求項に引用される幾つかの部材の機能を成し遂げてもよい。特定の手段が相互に異なる従属請求項において引用されているという単なる事実は、これらの手段の組合せが効果的に使用できないことを示さない。コンピュータプログラムは、他のハードウェアと共に又は一部として光記憶媒体又はソリッドステート媒体のような適切な媒体に格納／配布されてもよいが、インターネット、他の有線又は無線通信システムを介してのような他の形式で配信されてもよい。請求項内の参照符号は、範囲を限定するものとして解釈されるべきではない。
【０２０２】
参照文献
１ DellacらによるExpiratory flow limitation detected by forced oscillation and negative expiratory pressure, European Respiratory Journal Vol.29 No.2, pages363-3742
２ Darｏczy B,Hantos Z.によるAn improved forced oscillation estimation of respiratory impedance. Int J Biomed Comput 13:221-235, 1982
３ ForbesGW,AlonsoMA.によるConsistent analogs of the Fourier uncertainty relation. Am J Physics 69:340-347, 2001
４ Forbes GW, Alonso MA, SiegmanAE.によるUncertainty relations and minimal uncertainty states for the discrete Fourier transform and the Fourier series. J Phys A: MathGen 36:7027-7047, 2003
５ GaborD.によるTheory of communications. J Inst Electr Eng93:429-457, 1946
【０２０３】
表１
シンボル及び省略形
時間及び周波数
ｆｍ、ｆｎ 離散的周波数ｍ／Ｍ（又はｎ／Ｎ）
ｆＲ Δｆに対する基準周波数
ｆｎ 周波数ｆｎを持つ調和振動を記述するＮ―ベクトル
ｍ ＴＦＴに対する周波数指標
Ｍ（Ｍ’） ＴＦＴフィルタの幅（又は切り詰めたガウス分布を使用した有効幅）
ｎ ＯＤＦＴに対する周波数指標
Ｎ 時系列数のイベントの数
ＮＳ 確率過程からのｔ―ｆ領域のサンプルの数
ｒｔ、ｒｆ 時間円（又は周波数円）の半径
ｔ 離散的時間指標
ｔＲ Δｔに対する基準時間
Ａ，Ｂ，Ｃ 時間及び周波数の不確定度に関係する行列
Ｆ ＯＤＦＴ行列
Ｍ(Ｍｍ) ＴＦＴ行列（又は周波数ｆｍに対する部分的なＴＦＴ行列）
Ｔ タイムシフトオペレータ；ｃｉｒｃ｛０，．．．，０，１｝
Δｔ、Δｆ 時間又は周波数の不確定度
ω 複素指数；ｅｘｐ（２πｉ／Ｎ）
σｍｉｎ、σｍａｘ Ｃの最小及び最大の単一の値
Ω 主対角線上にωｔを持つ対角行列；ｄｉａｇ｛ω０，．．．，ωＮ−１｝
変数の名前
ｐ 圧力
ｑ 流量
ｘ 拡声器入力（又は、未知の変数）
ｚ 呼吸インピーダンス
変数のタイプ 例
離散的時間の関数としての決定論的な変数 ｑｔ
離散的時間の関数としてのランダム変数 Ｑｔ
決定論的なＮ―ベクトル（小文字） ｑ
ランダムなＮ―ベクトル（大文字） Ｑ
ｔ―ｆ領域のベクトル（周波数ｆｍが中心）
ｔ―ｆ領域のベクトル（周波数ｆｍが中心で時刻ｔから開始）
ｔ―ｆ領域の線形回帰により「予測される」決定論的なベクトル
ｔ―ｆ領域の線形回帰に対するランダムな「エラー」
ｔ―ｆ領域のエラーの推定値
リニアフィルタ
ｈｑｘ ｘからｑまでのフィルタの期分けされたインパルス応答を持つベクトル
Ｈｑｘ，ｎ 周波数ｆｎに対するｘからｑまでの伝達関数
Ｃｑｘ 拡声器入力と流量との間のフィルタを記述する巡回行列
Λｘｑ 主対角線上にＨｑｘ，ｎを持つ対角行列
統計値
ＡＱ，ｔ，ｍ
に対する信頼領域への流量のノイズの寄与を決定するｔ―ｆ領域の変数
Ｆα （２、η−２）自由度でのＦ検定のための上位１００（１−α）％の信頼限界
ＫＱｘ，ｔ，ｍ２ ｔ―ｆ領域の流量と拡声器入力との間の二乗サンプルコヒーレンス
時間及び周波数の関数としてインピーダンスｚｔ，ｍの二変量の最小推定値
流量のノイズがゼロであるという仮定の下のｚｔ，ｍの推定値
α タイプＩＩエラー
β 信頼領域の中心から
の偏差を決定する角度
η 等価自由度
μＱ 平均の流量
σＱ２ 流量のノイズの分散
通常の数学的シンボル
ｉ 虚数（ｉ２＝−１）
Ｉ ＮｘＮの恒等行列
≡ 規定により等しいこと
｜ｘｔ｜ 複素数値ｘｔの大きさ
Ｎ―ベクトルｘの２−ノルム
ｘＨ ベクトルｘのエルミート転置
ｘ＊ ｘの複素共役
｛・｝ シーケンス（通常、列ベクトルとして表される）
ｄｉａｇ｛・｝ 主対角線上にシーケンスを持つ対角行列
ｃｉｒｃ｛・｝ 第１の行にシーケンスを巡回行列
省略形
ＦＯＴ 強制的振動技術
ｔ―ｆ領域 時間―周波数領域
ＴＦＴ 時間―周波数変換
ＯＤＦＴ 正規直交離散フーリエ変換
ＲＶ ランダム変数
【技術分野】
【０００１】
本発明は、強制的圧力振動から呼吸インピーダンスを推定するための方法と装置とに関する。
【背景技術】
【０００２】
人間の呼吸器系の呼吸又は音響インピーダンスは、気道、肺及び胸壁の抵抗、追従及び慣性に関する情報を得るために測定できる。この情報は、慢性閉塞性肺疾患（ＣＯＰＤ）、喘息及び気管支炎のような様々な呼吸器疾患の性質及び発病度を診断する際に有効である。
【０００３】
ページ363-374のEuropean Respiratory Journal Vol.29 No.2のＤｅｌｌａｃａ等による「Expiratory Flow Limitation Detected by Forced Oscillation and Negative Expiratory Pressure」（参考文献１）に説明されているような強制的振動技術（ＦＯＴ）では、人が正常に呼吸している間、音響波が呼吸器系に導かれ、反応が呼吸インピーダンスを決定するために測定される。
【発明の概要】
【発明が解決しようとする課題】
【０００４】
呼吸インピーダンスは、流量及び圧力に関して音響波から生じている振動における周波数依存的な関係を記述する。呼吸インピーダンスが（幾つかの疾患及び他の病状におけるような）吸気から呼気まで変化する所で、呼吸インピーダンスは、精細な時間解像度で推定されなければならない。しかしながら、従来技術では、生理的目的のために充分に短い（すなわち、吸入又は呼気の期間より短い）時間間隔の呼吸インピーダンスを推定するための技術の信頼性に対して、要件がほとんど与えられていなかった。
【０００５】
従って、短い時間間隔の呼吸インピーダンスを確実に推定するための方法と装置とのためのニーズがある。
【課題を解決するための手段】
【０００６】
本発明の第１の態様によると、被験者の気道と連通する患者インターフェイス装置内に圧力波を生成するステップと、流量及び圧力を表わすそれぞれの時系列を作るため被験者の気道と患者インターフェイス装置とを含む空気制御システム内のガスの流量及び圧力を決定するステップと、それぞれの時系列を時間−周波数領域へ変換して変換された時系列を作るステップと、それぞれの変換された時系列から流量及び圧力のパワーを推定するステップと、変換された時系列に基づいて流量及び圧力のそれぞれのクロススペクトルを推定するステップと、推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定するステップとを有する呼吸インピーダンスを推定する方法が提供される。
【０００７】
本発明の第２の態様によると、患者インターフェイス装置と、被験者の気道のガスの圧力、流量又はボリュームの振動を生成するための励振源と、患者インターフェイス装置及び被験者の気道により規定される空気圧回路のガスの流量及び圧力を決定し、流量及び圧力を表しているそれぞれの時系列の値を出力するための手段と、プロセッサとを有し、前記プロセッサは、それぞれの時系列を時間―周波数領域へ変換し、それぞれ変換された時系列から、流量及び圧力のパワーを推定し、それぞれ変換された時系列に基づいて、流量及び圧力のそれぞれのクロススペクトルを推定し、推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定する、呼吸インピーダンスを推定するための装置が提供される。
【０００８】
本発明の第３の態様によると、適切なプロセッサ又はコンピュータによる実行の際、前記プロセッサ又はコンピュータが上述の方法を実施するコンピュータ可読媒体に組み込まれたコンピュータ可読のコードを具備する、コンピュータプログラムが提供される。
【０００９】
呼吸インピーダンスの推定は、気道の妨害を評価するか、又は疾患の発病度を推定するために推定を使用する診断ツールによりなされる。従って、診断ツールは、また、呼吸インピーダンスに影響を及ぼすべき治療の（薬理学的又は他の）効果を評価するためにも用いられる。
【００１０】
よって、第４の本発明の態様によると、上述されたような呼吸インピーダンスを測定するステップと、測定された呼吸インピーダンスを基礎として生理的状態を診断するステップとを有する生理的状態を診断する方法が提供される。
【００１１】
呼吸インピーダンスの推定は、また、又は代わりに、病状の治療に使用される機械、例えば気道閉塞を打ち消すために用いられる非侵襲性のベンチレータの設定を適応させるため、又は、生理的状態に対する治療法（例えばどの特定の薬物又は装置を使用するか、薬物投与量等）を決定する際に用いられる。
【００１２】
よって、本発明の第５の態様によると、上述されたように呼吸インピーダンスを測定するステップと、測定された呼吸インピーダンスに基づいて生理的状態の治療を決定し及び／又は処置するステップとを有する、生理的状態の治療方法が提供される。
【００１３】
本発明の好ましい実施例が、以下の図面を参照して、単なる例示として、ここで詳細に説明されるだろう。
【図面の簡単な説明】
【００１４】
【図１】図１は、本発明による強制的圧力振動から呼吸インピーダンスを推定するための装置のブロック図である。
【図２】図２は、本発明による呼吸インピーダンスを推定するために用いられる線形モデルの模式図である。
【図３】図３Ａ及び図３Ｂは、それぞれ離散的時間領域の巡回性推定及び不確定度を例示する。
【図４】図４Ａ及び図４Ｂは、Δｔ２＋Ｎ２Δｆ２が最小である時系列を例示する。
【図５】図５は、時間―周波数領域の二変量の最小二乗法を例示する。
【図６】図６は、複素平面内の信頼領域を例示する。
【図７】図７は、５つの異なる周波数での強制的圧力振動に対する信頼限界を持つ時間及び周波数の関数として、健常な患者の呼吸インピーダンスを例示する。
【図８】図８は、５つの異なる周波数での強制的圧力振動に対する信頼限界を持つ時間及び周波数の関数として、ＣＯＰＤを持つ患者の呼吸インピーダンスを例示する。
【図９】図９は、本発明による呼吸インピーダンスを推定する方法のステップを図示する流れ図である。
【図１０】図１０は、本発明による装置により実施される処理ステップをより詳細に図示する流れ図である。
【図１１】図１１は、時系列のイベントの数Ｎの関数として、行列の最小及び最大の特異値を示す。
【図１２】図１２は、Ｎ＝１６に対するＣＨＣの固有ベクトルの重心と時間Δｔでの不確定度との間の関係を例示する。
【図１３】図１３は、Ｎ＝１６に対する周波数Δｆでの不確定度と時間Δｔでの不確定度との間の関係を例示する。
【図１４】図１４は、ｆｎに対するνｍｉｎ，ｎをプロットしているグラフである。
【図１５】図１５は、Ｎ＝１６に対するＮΔｆ対Δｔのプロットにおける不確定度を例示する。
【図１６】図１６は、複素平面において推定値Ｚｔ，ｍに対する１００（１−α）％信頼限界の推定を例示する。
【図１７】図１７は、健常な被験者に対する時間及び周波数の関数として二乗コヒーレンスを例示する。
【図１８】図１８は、ＣＯＰＤを持つ被験者に対する時間及び周波数の関数として二乗コヒーレンスを例示する。
【発明を実施するための形態】
【００１５】
上述のように、呼吸器系の音響インピーダンスは、マウスピース又はフェースマスクの強制的圧力振動（強制的振動技術、ＦＯＴ）から推定されることができる。この「呼吸インピーダンス」は気道、肺及び胸壁の抵抗、追従及び慣性に依存するので、呼吸インピーダンスは、様々な呼吸器疾患の発病度及び性質に対する洞察を提供する。周波数依存的なインピーダンスは、しばしば時間的に変動する。慢性閉塞性肺疾患において、気道抵抗は、呼気の間に通常増大する。睡眠無呼吸又はいびきにおいて、抵抗は、吸気の間に増大する。これは、短い時間間隔（呼吸の期間より短い）におけるインピーダンスの信頼性が高い推定を要する。
【００１６】
呼吸インピーダンスは、周波数ｆの関数として、流量ｑと圧力ｐとの間の線形関係を記述する複素数値の伝達関数として規定される。
ｐ（ｆ）＝ｑ（ｆ）ｚ（ｆ）（１）
ここで、インピーダンスはｚ（ｆ）により示される。この式は、ｚ（ｆ）が一定で、時間が無限である場合だけ成り立つ。システムが一時的に安定の場合、ｚ（ｆ）は短い時間間隔で依然推定される。しかしながら、これは、信号が時間及び周波数領域両方において鋭くは特定できないと言う不確定性原理により制限される。加えて、短い時間間隔の使用は、推定を偶然のイベントに益々依存させる。問題は、これが、呼吸サイクルの一部で安定なだけであるｚ（ｆ）のような伝達関数の推定にどのくらい影響を及ぼすかである。
【００１７】
以下では、一時的に安定なプロセスの伝達関数が、離散的時間―周波数領域の線形回帰により推定される。不確定性関係は、離散的時間及び周波数に対して導出され、最適時間―周波数解像度を得るために用いられる。使用されるインピーダンス推定器の統計特性の分析は、呼吸メカニクスの周波数依存的なパラメータを得て、その推定された値及び信頼限界両方が時間的にフォローできる。
【００１８】
発明の詳細な説明の後に含まれるテーブル（表１）は、本発明の以下の説明で使用されるシンボル及び略語一覧を提供する。
【００１９】
方法
発明の詳細な説明のこのセクションは、本発明による方法及び装置の数学的基礎を説明しておく。本発明のより数学的でない説明は、以下の「説明」セクションにおいて提供される。
【００２０】
本発明による装置２の模式図が図１に示される。簡潔には、使用されるＦＯＴ装置２は、被験者の気道を含む空気制御システムを作るために被験者の気道に接続される患者インターフェイス装置４から成る。患者インターフェイス装置４は、マスク又はマウスピースのような被験者の気道との気送接続を供給するために適切な任意の装置である。患者インターフェイス装置４は、比較的低い周波数、例えば８Ｈｚから２４Ｈｚまでの周波数を持つ圧力波、流量変化又はボリューム変化を生成する励振源６に動作的に結合される。本発明の例示的な実施例では、患者インターフェイス装置４は拡声器である。
【００２１】
図示の実施例では、患者は、患者インターフェイス装置４と励振源６との間に位置されるニューモタコヘッド８に近いメッシュ有線の抵抗９を通じて空気を呼吸する。強制的波は、患者の気道及び肺において部分的に反射され、ＦＯＴ装置２において測定されるエアフロー及び圧力の振動を導く。ニューモタコヘッド８を通るエアフローは、差圧トランスジューサ１０を使用して測定され、圧力はマウスピース４の近くの第２の圧力トランスジューサ１２を使用して測定される。ＡＤコンバータ１２は、流量及び圧力のそれぞれの時系列｛ｑｔ｝｛ｐｔ｝を作り、分析のためプロセッサ又はコンピュータ１６に供給され、ここで、整数ｔは離散的時間指標である。コンピュータ又はプロセッサ１６は、ＡＤコンバータ１４及びアンプ１８を介して拡声器６により生成される音響波の周波数を制御するための信号も供給できる。装置２の特定の実施例の更なる詳細は、「測定装置」と名付けられた以下のセクションで提供される。
【００２２】
本発明は、また、他のガス流送出装置が、励振源６に加えて被験者の気道に結合できることを考察する。例えば、圧力サポートシステム又はベンチレータは、例えば、励振源６が被験者の気道を振動させている間、ガスの流れを供給するために患者インターフェイス装置４に結合できる。また、圧力サポートシステム又はベンチレータ内で圧力及び流量を測定することを含む、流量及び圧力測定が、空気制御システムに沿って任意の位置でなされ得ることに留意されるべきである。その後、測定された圧力及び流量は、任意の従来の技術を用いて患者の気道のような空気制御システムのガスの流量及び圧力を決定するために使用できる。
【００２３】
本発明は、また、励振源６が独立型装置よりもむしろ圧力サポートシステム又はベンチレータの部品であることも考察する。例えば、圧力サポートシステム又はベンチレータは、患者へ送られるガスの流量／圧力を制御するためのバルブを含む。本発明は、低周波圧力波、流量変化又はボリューム変化をもたらすために斯様なバルブを振動させることを考察する。もちろん、低周波圧力波、流量変化又はボリューム変化を生成できるピストンのような他の任意の装置又はシステムも、励振源６として使用できる。
【００２４】
圧力及び／又は流量センサを使用して圧力及び／又は流量を測定するよりはむしろ、本発明は、また、流量又は圧力が、圧力サポートシステム又はベンチレータによって設定又は制御できることも考察する。この場合、設定された圧力又は設定された流量は、差圧トランスジューサ１０又は第２の圧力トランスデューサ１２の出力としてよりもむしろ現目的のための圧力又は流量として使用される。
【００２５】
２つのサンプルの間の間隔（秒）は、サンプルレート（Ｈｚ）の逆数である。このセクションでは、時間依存及び周波数依存インピーダンスが、比較的単純な線形代数学を使用してこれらの時系列から導かれる。線形代数学をよく知らない読者は、主要な結果が言葉で要約されている以下の「説明」セクションを参照されたい。
【００２６】
単純な線形モデル偶然のイベントを説明するために、二変量の時系列｛ｑｔ、ｐｔ｝が二変量の確率プロセス｛Ｑｔ、Ｐｔ｝の実現であると仮定する。決定論的な変数は小文字で書かれ、ランダム（確率）変数（ＲＶ）は大文字で書かれる。時間が無限に長い（ｔ＝．．．，−１，０，１．．．）と更に仮定する。ここで図２を参照すると、励振源入力ｘｔは、インパルス反応シーケンス｛ｈｑｘ，ｌ｝及び｛ｈｐｘ、ｌ｝を持つ線形時間不変のフィルタ（２及び３）を通ってフィルタリングされる。これは、以下の線形モデルにより予測されるように、流量及び圧力を与える。
（２）
【００２７】
これらの「真」の値がそれぞれの分散σＱ２及びσＰ２を持つゼロ平均ホワイトノイズの２つの独立源により妨げられると仮定する。ノイズのこれらの源は、「エラー」項Ｑｅ，ｔ及びＰｅ，ｔを与えるために、線形時間不変のフィルタ（１及び４）も通過する。完全なモデルは、このとき、以下のように書かれる。
（３）
ここで、平均値μＱ及びμＰは定数である。結果として生じるＱｔ及びＰｔが正常に分配されるので、仮定された二変量のプロセス｛Ｑｔ、Ｐｔ｝は、完全に変動しない（静止している）。
【００２８】
循環性仮定
説明された静止プロセスが無限に長い一方、静止がＮ個のイベントの有限シーケンスに対してのみ仮定できている過渡現象を説明することが発明の目的である。これは、最後のイベントが最初のイベントに先行するという点で、時間が循環していると仮定することにより解決できる。結果として生じる時系列は、このとき、時計のような円に沿って定期的にＮイベントを配置することにより表される。円から離れずに、同一方向で円を永久に追求できる。結果として生じる無限の時系列は、Ｎで周期的である。すなわち、任意の整数ｊに対して、以下の通りである。
ｘｔ＝ｘｔ＋ｊＮ （４）
【００２９】
ｕ＝ｔ−ｌ＋ｊＮと置き換えると、式２の第１のコンボリューションは、以下のように書きかえられる。
【００３０】
仮定された周期的キャラクタｘｔ（式４）に対して、ｘｕ−ｊＮ＝ｘｕである。上記の式の無限和は、「長さＮに周期化されたｈｑｘ，ｔ−ｕとして、以下の式のように示されて規定できる。
（５）
結果として、
（６）
【００３１】
ｈｑｘ，ｔ−ｕｏもＮにおいて周期的であるので、（インデックスはゼロからＮ−１の範囲に常にある（これは便利である）ように）ｔ−ｕが負である場合、インデックスｔ−ｕはＮの分増大できる。シーケンスが列ベクトルにより表されるとき、式６は以下のように書ける（Ｎ＝４である場合）。
（７）
【００３２】
この式のＮｘＮ行列は、巡回行列式（端で循環しながら、前の行が右にシフトするバージョンの行を持つ正方行列）である。第１の行がｈｑｘＨにより示されるとき、肩文字Ｈは、エルミート転置を表わし、巡回行列はｃｉｒｃ｛ｈｑｘＨ｝とも書かれる。
、
及び
とする。このとき、ｑｐ＝Ｃｑｘｘである。同様に、ｐＰ＝Ｃｐｘｘなので、式３は以下のようになる。
（８）
【００３３】
離散的周波数領域
周波数依存的なインピーダンスを導くために、式８のモデルは、離散的周波数領域へ変換される必要がある。時間―周波数分析へのステップアップとして、これは、ここで短く要約される。所与の時間系列
から離散的周波数領域への変換は、時系列を離散的周波数
（ｎ＝０，．．．，Ｎ−１））を持つ一組の調和振動に区切る。斯様な振動は、以下のようなＮ―ベクトルにより記述される。
（９）
ここで、ω≡ｅｘｐ（２πｉｎｔ／Ｎ）及びｉ２＝−１である。オイラーの定理により、
ｅｘｐ（２πｉｎｔ／Ｎ）＝ｃｏｓ（２πｎｔ／Ｎ）＋ｉｓｉｎ（２πｎｔ／Ｎ） （１０）
よって、ベクトルｆｎは、実際に、周波数ｆｎを持つ実軸及び虚軸に沿った複合振動を説明する。容易に以下の式が示される。
（１１）
【００３４】
結果として、二乗ノルム
はＮに等しい一方、異なる調和周波数の２つの振動は直交する。
【００３５】
ｘの離散的周波数領域への変換は、以下の式により規定される「正規直交離散的フーリエ変換」（ＯＤＦＴ）行列により実施される。ここで、ｎｔ＝０，．．．，Ｎ−１である。
（１２）
この変換は、Ｎ―ベクトルＦｘを作り、これらの成分は、内積
である。式１１から、Ｆの行（及び列）が正規直交であるので、行列はユニタリであることがわかる。
ＦＦＨ＝ＦＨＦ＝Ｉ （１３）
ここでＩは、ＮｘＮの単位行列である。これは、Ｎ振動からｘの合成を可能にする。
（１４）
【００３６】
従って、ｘは、周波数ｆｎを持つＮの調和振動の重み付けられた和として書かれる。
【００３７】
ベクトルｆｎは巡回行列の固有ベクトルであるので、以下のように書くことができる。Ｃｑｘｆｎ＝ｆｎＨｑｘ，ｎ （１５）
【００３８】
複素数Ｈｑｘ，ｎは、以下の式により規定されるように、周波数ｆｎに対してｘからｑまでの伝達関数である。
（１６）
【００３９】
Ｈｐｘ、ｎが同様に規定されるとき、周波数ｆｎに対するインピーダンスは、以下の式である。
（１７）
【００４０】
固有値が対角行列
に入れられ、列として対応する固有ベクトルをＮｘＮ行列ＦＨへ入れられるとき、式１５は以下の式になる。
（１８）
【００４１】
巡回行列は、正規行列である（これらは、これらのエルミート転置行列と代えられる）。従って、式１８は、ＮｘＮ正規行列がＮ個の正規直交固有ベクトル（この場合、Ｎ個の正規直交列ベクトルｆｎ）を持つというスペクトル定理と一致している。
【００４２】
離散的時間及び周波数に対する不確定性原理
の成分は、時間領域において鋭く局所的であるが、周波数領域においては不確かである。これとは逆に、フーリエ変換Ｆｘの成分は、周波数領域において鋭く局所的であるが、時間領域においては不確かである。両方の領域で小さな不確定度を持つ信号を得るために、不確定さの尺度が、規定されなければならない。
【００４３】
時間が巡回性であると仮定すると、フォーブズ及びアロンゾのアプローチが、最も適当なようである（参照文献３を参照）。ここで図３Ａを参照すると、時間領域で、二乗成分｜ｘｔ｜２は、各イベントの重みに対応する図３Ａの各質点のサイズを持ち、原点で中央に置かれる、複素平面内の半径ｒｔを持つ円上の規則的な間隔の質点として置かれる。時間は、２つのイベントの間の円弧長により表される。円の周囲は全体の期間Ｎ（この具体例ではＮ＝１２）であるので、半径は
に等しい。時刻ｔでの所与のイベントは、複素平面のポイントｒｔωｔに位置される。フォーブズ及びアロンゾ（参照文献３）は、不確定度の幾つかの尺度を提案した。ここで、わずかに異なる尺度が使われ、これは、時間ｔでの各イベントと基準時間ｔＲでのイベントとの間の二乗距離の加重平均として式１９に規定される。
（１９）
【００４４】
不確定度のこの尺度は、二次形式の比率として容易に表されることができるという利点を持つ。一般性の損失なしに、ｔＲ＝０（図３Ａに示されるように）と仮定できる。このとき、二乗された不確定度Δｔ２は、以下の式である。
（２０）
ここで、対角行列Ω≡ｄｉａｇ｛１，ω１，．．．，ωＮ−１｝である。下記の式を規定することも可能である。
（２１）
ここで、Ａ≡Ω−Ｉである。物理的に、Δｔ２は、基準ポイントＲ（時間ｔＲに対応する）を通る（円の平面と直角をなす）軸の周りの質点のセットの第２の慣性モーメントである。図３Ｂに図示されるように、この慣性モーメントは、重心に直接関係する。幾何学的に、Δｔは、Ｒから、重心を通る垂直線が円と交差する２つのポイントのうちの一方のポイントまでの距離に等しい。それで、Δｔは、重心からＲまでの水平距離により、完全に決定される。重心がＲと一致する場合、ｘの重さは時間ｔＲに完全に集中し、Δｔは、ゼロである。重心が円の反対のポイントにある場合、重さはそのポイントに完全に集中し、Δｔは最大であり、２ｒｔに等しい。よって、図３Ｂにおいて、重心Ｂを通る垂直線は、ポイントＣで円と交差する。不確定度Δｔは距離
と等しい。
【００４５】
周波数領域において、Ｆｘの成分は、周波数ｆｎ≡ｎ／Ｎと関連している。同じ成分が周波数ｎ／Ｎ＋ｊに対して得られ、ここで、ｊは任意の整数である（式９、式１０を参照）。従って、Ｆｘの成分は、単位ユニット期間を持って、周波数領域において周期的である。このように、不確定度の匹敵する尺度は、基準周波数ｆＲに関係して、周波数領域で規定できる。Ｆｘの成分の二乗の大きさは、複素平面の円上の質点として、再び配置される。ここで、半径は、ｒｆ≡１／（２π）に等しいので、円周はユニティ、１に等しい。ｆＲ＝０の場合、下記の式である。
（２２）
【００４６】
対角行列Ωは、「時間シフトオペレータ」Ｔ≡ｃｉｒｃ｛０，．．．，０，１｝に関係する（２２）。Ｔによる左からのｘの乗算は、巡回的にｘの成分を一つ下方へシフトさせる。式１８から容易にＦＴ−１＝ΩＦ（「シフト定理」）であることがわかるので、以下の式となる。
【００４７】
Ｆがユニタリであるので、以下の式となる。
とすると、以下の式となる。
（２３）
【００４８】
ここで、以下の行列式を考える。
【００４９】
時間及び周波数の全体の不確定度は、以下の比率の式により表わされる。
（２４）
【００５０】
この比率（レイリー比率）は、Ｃの最小二乗特異値から最大二乗特異値までの範囲の値をとることができる。これは、離散的時間及び周波数に対する対応する不確定性原理を、以下の式のように作る。
（２５）
【００５１】
従って、全体の不確定度Δｔ２＋Ｎ２Δｆ２が少なくともσｍｉｎ２に等しくなければならないので、Δｔ及びΔｆ両方が非常に小さいことは不可能である。全体の不確定度は、ｘがσｍｉｎ２に対応するＣＨＣの固有ベクトルである場合、最小である。ＣＨＣの固有ベクトルは、離散的マチュー関数として考えられる。ＣＨＣの興味深い特性は、Ｆと代われることであるので、ＣＨＣとＦとは、同じ固有ベクトルを共有する。ｖｍｉｎに対応するＦの固有値はユニタリに等しいので、Ｆｖｍｉｎ＝ｖｍｉｎである。図４Ａに示されるように、これは、ｘがｖｍｉｎと等しい場合、Δｔ＝ＮΔｆを意味する。図４は、合計Δｔ２＋Ｎ２Δｆ２が最小である時系列を示す。図４Ａでは、固有ベクトルは、基準時間ｔＲ＝０及び基準周波数ｆＲ＝０（上のグラフ）に対して導出される。これは、フーリエ行列（単位固有値を持つ）の固有ベクトルでもあるので、時系列は離散的周波数領域（下のグラフ）への変換後、不変のままである。これらのグラフにおいて、時間軸及び周波数軸の両方とも、視覚的理由のために巡回的にシフトされる。図４Ｂは、基準周波数ｆＲ＝１／４である以外、図４Ａに対応する。図４Ｂでは、黒丸が実数値を表し、白丸が虚数値を表し、十字が絶対値を示す。
【００５２】
不確定性原理のより詳細な分析のため下記のセクション「不確定性原理の導出」を参照されたい。
【００５３】
離散的時間―周波数領域
ＯＤＦＴは、Ｎ個のイベントの時系列をＮ個の異なる周波数を持つ一組の振動へ区切られる。時間―周波数分析を実施する簡易なやり方は、短く続くシーケンスのＭ（Ｍ＜Ｎ）個のイベントに対してＯＤＦＴを繰り返し計算することである。従って、系列が周波数ｆｍ≡ｍ／Ｍ（ｍ＝０，．．．，Ｍ−１）を持つＭ個の振動に区切られ、これは離散的「短時間フーリエ変換」になる。各短く続く振動がｖｍｉｎの成分で「窓」をつけられる（ウィンドウ化される）場合、時間及び周波数についての結果的な不確定度は最小である。
【００５４】
Ｍ個のイベントのシーケンスに対して、斯様なウィンドウ化された振動は、以下の式のようにＭベクトルにより記述される。
（２６）
【００５５】
例が、図４Ｂに示される。ｖｍｉｎ（ｍ）は、周波数領域において実際にｖｍｉｎの巡回的にシフトされたバージョンであることに留意されたい。図４に例示されて、下記のセクション「不確定性原理の導出」で証明されているように（式Ｂ１９を参照）、これは、（ｆｍが基準周波数として採用される場合）固有ベクトルと同じ時間及び周波数の不確定度を持つ。ｖｍｉｎ（ｍ）の成分をｖｍｉｎ，ｔ（ｍ）により示すとする。短く続く振動を記述するＮベクトルは、このとき、ｖｍｉｎ（ｍ）の中間にゼロを挿入することにより得られる（よって、ベクトルの単一モードの構造が完全なままである）。これは、以下の式を与える。
（２７）
ここで、ＫはＭ／２以下の最も小さな整数である。この振動は、時間ｔ＝０及び周波数ｆｍ＝ｍ／Ｍを中央に置く。一般に、以下の式
（２８）
の振動は、時間ｔ及び周波数ｆｍを中央に置く。
【００５６】
最小二乗法という意味では、ｘとｈｍ，ｔとの間の最も適合する線形関係は、ｈｍ，ｔによりスパンされた（複素）一次元の部分空間上へのｘの直交射影からわかる。ｈｍ，ｔが単位ベクトル（式２７を参照）として規定されることを考慮すると、この射影は、以下の式となる。
（２９）
【００５７】
所与の周波数ｆｍに対して、ＮｘＮ行列Ｍｍを以下の式のように規定する。
（３０）
【００５８】
これは、共振周波数ｆｍを持つ帯域フィルタを行う巡回行列である。Ｍｍによる左からのｘの乗算は、波形符号又は「波符号」により示されるＮベクトルを与える。
（３１）
【００５９】
の成分は、内積
である。
【００６０】
十分な時間―周波数変換（ＴＦＴ）は、ｍ＝０，．．．，Ｍ−１に対する全ての周波数ｆｍを含む。ＭＮｘＮ行列のＭを以下の式のように規定する。
（３２）
【００６１】
離散的時間―周波数領域（時間―周波数領域）への変換は、（所与のＭに固有な）Ｍによる左からのｘの乗算により実施される。これは、ＭＮ―ベクトル
を与える。バック変換は、再びＭＨによる乗算を意味し、
を与える。下記の「時間―周波数変換の導出」というセクションにおいて、Ｍの列は正規直交であることが示されるので、以下の式となる。
ＭＨＭ＝Ｉ（３３）
【００６２】
ＭＨＭｘ＝Ｉｘ＝ｘなので、よって、バック変換は元の時系列を回復する。これは、以下の式に従って、ｘの時間―周波数合成を可能にする。
（３４）
【００６３】
このように、時系列はＭＮ個の短く続く振動の加重加算として書かれ、各々が時間ｔ及び周波数ｆｍに中心が置かれる（式１４を参照）。所与のＭに対して、これらの振動は、式２５の不確定性原理に従って、最適時間―周波数解像度を与える。下記の「不確定性原理の導出」セクションに示されるように、振動は、Ｍ’＝１０・Δｔゼロ以外の要素で丸められたガウス分布により近似できる。
【００６４】
励振源への入力
励振源入力ｘとして、一組の調和振動を使用することが便利である。データを分析するために、このとき、これらの振動の周波数がＴＦＴ共振周波数ｆｍと一致するように、ＴＦＴのスケールＭが選択できる。最初に、周波数ｆｎ及び振幅ａｍを持つ単一の正弦波入力を考える。オイラーの定理（式１０）により、実数値の調和振動は、複素数値とその複素共役との和である。よって、励振源入力は、下記の式のように書ける。
（３５）
ここで、星印は複素共役を示す。式９及び式１０から、ｆｎ＊＝ｆＮ−ｎがわかる。ＮがＭの整数倍である場合、周波数ｆｎ≡ｎ／Ｍは、ＴＦＴ共振周波数ｆｍ≡ｍ／Ｍとマッチできる。その後、ｆｎ＝ｆｍであるように、ｍが選択される場合、以下の式となる。
（３６）
【００６５】
Ｍが十分に大きく、対応するΔｆが、共振周波数ｆｎ＝ｆｍを持つＴＦＴフィルタにより周波数ｆＮ−ｎで振動を抑制するのに十分小さい場合、後の近似値は有効である。「相補的周波数」ｆＮ−ｎ（ここでは、ｆＭ−ｍに等しい）での成分も式３６でＭＭ−ｍをＭｍに加算することにより含まれる場合、時間―周波数領域の励振源入力の十分な表現が得られる。
【００６６】
式８のモデルは、Ｍにより両方の式を左から乗算することにより、時間―周波数領域に変換される。対象の周波数（例えば、式３５により記述された強制的振動の周波数）に注目すると、Ｍｍにより式８の上の部分を乗算することは、以下の式となる。
（３７）
ｍがゼロでない場合、ＭｍμＱがゼロであることが容易に示される。よって、
（３８）
【００６７】
巡回行列全てが代えられる
であることがわかる。ｆｎ＝ｆｍであるようにｎ及びｍが選ばれる場合、このとき、式３６及び式１５を用いて、以下の式となる。
（３９）
【００６８】
ｆｎ＝ｆｍに対するＨｑｘ、ｎがｂＱ、ｍにより示される（対応するＨｐｘ、ｎがｂＰ、ｍにより示される）とき、式８のモデルは、以下の式に低減される。
（４０）
【００６９】
ｆｎ＝ｆｍに対するインピーダンスは、以下の式となる。
ｚｍ≡ｂＰ、ｍ／ｂＱ、ｍ （４１）
【００７０】
これらの振動が時間―周波数領域においてよく離隔されるほど、使用されるＴＦＴのΔｆが十分に小さいならば、他の振動のセットは、式３５の励振源入力に加えられる。
【００７１】
時間―周波数領域の分散分析
ＴＦＴは、時間―周波数領域の分散の分析を可能にする（「時間―周波数領域の分散分析」というこのセクションを参照）。例えば、流れのノイズの二乗ノルム（式８を参照）は以下の式に区切られる。
（４２）
【００７２】
のランダム変数は、ノイズＱｅの「ＴＦＴサンプルパワースペクトル」と呼ばれている。転送がノイズ１からＱｅ，ｔ（図２）までのフィルタ１の伝達関数が、ｆｍについてのＴＦＴフィルタの通過帯域で相対的に平坦であると仮定する。このとき、
は、あたかもゼロ平均及び分散σＱ、ｍ２Ｉを持つホワイトノイズから導出されるかのように、ほぼ分散される。結果として、
は、等価自由度ηを持つ近似のカイ二乗分布に従う。（相補的周波数ｆＭ−ｍの成分と同様に）循環仮定から独立している
のそれらの成分だけが含まれているηの導出のため、以下のセクション「時間―周波数領域の分散分析」を参照されたい。
【００７３】
時間―周波数領域の二変量の最小二乗法
各時間ｔで、図２のモデルにより記述されている、しばらくの間変動しないプロセスが始まると、ここで仮定される。当該プロセスは、Ｎ個のイベントを含み、循環的である（すなわち、Ｎで周期的）とみなされる。時間―周波数領域では、モデルは（式４０から）以下の式となる。
（４３）
ここで、ベクトルは時間ｔで始まる「短く続く」Ｎ個のベクトルである。定数ｂＱ、ｍ及びｂＰ、ｍは、これらがプロセスが始まる時間に特有であるとみなされるので、同様にインデックスｔを得た。従って、モデルへの入力と同様に定数は、（２つのプロセスが時間的に重複する場合であっても）異なる時間ｔで異なる。
【００７４】
ｂＱ、ｔ、ｍ及びｂＰ、ｔ、ｍの最も適切な推定値は、エラーベクトルの二乗ノルムを最小化することにより、
及び
から導出できる。これは、最小二乗推定値
及び
を生じる。推定された関係は、以下の式である。
（４４）
ここで、推定された「予測された」ベクトル
は、推定された「エラー」
に正規直交している（同様に、
は
と正規直交している）。図５では、ベクトルは矢印として示され、全ての変数は時間インデックスｔ及び周波数インデックスｍ（周波数ｆｍ＝ｍ／Ｍに対応する）を参照している。特に、図５では、励振源入力
は一定の変数である一方で、流量
及び圧力
はランダム変数である。ベクトル
及び
は、
によりスパンされた複素一次元の部分空間（ラインｌ）への
及び
の射影である。インピーダンス推定値は、
から分かる。エラーベクトル
及び
は
に垂直である。各ベクトルは、長さがノルムに等しい矢印により表される。ベクトル間の角度は、ｃｏｓ２φ及びｃｏｓ２θが、
とそれぞれ
及び
との間の二乗コヒーレンスに等しいように描かれる。
【００７５】
と
との間の線形関係の程度は、以下の式のような「ＴＦＴは二乗サンプルコヒーレンス」の表現になる。
（４５）
【００７６】
ランダム変数
は、ｘからＱまでの「ＴＦＴサンプルクロススペクトル」と呼ばれている。幾何学的に、ＫＱｘ，ｔ，ｍ２はｃｏｓ２φと等しく、ここで、φは図５の
から
までの角度である。同様に、ＫＰｘ，ｔ，ｍ２はｃｏｓ２θと等しい。最小二乗推定値は、標準的やり方で導出される。これは、以下の式を与える。
（４６）
【００７７】
従って、
は、ｘからＱまでのサンプルクロススペクトルとｘのサンプルパワースペクトルとの比率である。
が同様のやり方で導出されるとき、インピーダンスは以下の式により推定される。
（４７）
【００７８】
「真の」インピーダンスｚｔ、ｍは、実数部分ｒｔ、ｍ（「抵抗」）及び虚数部分ｘｔ、ｍ（「リアクタンス」）を持つ。対応する推定値は、以下の式、
の実数部分及び虚数部分である。
（４８）
【００７９】
図６Ａに示されるように、ｂＱ、ｔ、ｍに対する１００（１−α）％信頼領域は、
を中心にした複素平面の円により区切られている。信頼領域は、以下の式により与えられる。
（４９）
ここで、
であり、ＦαはＦ分布の上の１００（１−α）％である（以下のセクション「信頼限界の導出」を参照）。励振源入力ｘｔは、無次元変数であり、
は（図５から
なので）流量のユニットを持つ。循環仮定に依存する値は、ここで放棄される（サンプルの有効数は、ＮＳ＝Ｎ−Ｍ＋１である）。原点が円内にある場合（このとき、
は１００・α％信頼レベルでゼロと著しくは異ならない）、関係は拒絶される。よって、ＡＱ、ｔ、ｍ＜１の場合、有意性が得られることが、式４９から分かる。
に対する信頼領域は比較可能な円により区切られ、
と
との間の関係は、ＡＰ、ｔ、ｍ＜１の場合、重要であると考えられる。
【００８０】
に対する信頼領域の控え目な推定値は、
及び
が非相関であるという仮定の下、以下の「信頼限界の導出」セクションで導出される。当該領域は、以下の式により記述される。
（５０）
【００８１】
図６Ｂを参照すると、限界円は推定値
の周りに同心円ではない。
の実数部分及び虚数部分が、
及び
により表わされる。垂直矢印は、実数部分に対する上位の信頼限界を示し、水平矢印は、インピーダンスの虚数部分に対する下位の信頼限界を示す。
の所与の実現のための「真の」ｚｔ，ｍの分布は、明らかに斜めに対称形である。図６Ｂでは、歪度は角度βにより決定される。歪度は、図５のβと角度φとの間の関係から明らかなように、完全に流量のノイズに起因する。
（５１）
【００８２】
流量のノイズがゼロである場合、このとき、φ＝０、ＫＱｘ、ｔ、ｍ２＝１及びＡＱ、ｔ、ｍ２＝０であるので、β＝０であり、信頼領域が
について対称形である。
【００８３】
下記の「信頼限界の導出」セクションは、流量のノイズ又は圧力のノイズがゼロであると仮定される場合に起こることの更なる分析を与える。それぞれの推定値は、このとき、以下の式に等しい。
（５２）
【００８４】
これは、時間―周波数領域における流量から圧力まで（又はその逆）の線形回帰となる。変数
及び
は、それぞれのサンプルパワー及びクロススペクトルである。式５２の推定値と
との主要な違いは、推定されたインピーダンスの大きさにある。大きさは、以下の式を通じて関係付けられる。
（５３）
【００８５】
ノイズが流量及び圧力に存在する場合、インピーダンスの大きさは、
が
の代わりに使用される場合過小評価される（
が使用される場合過大評価される）。
【００８６】
設定
装置２の使用中、流量及び圧力が、８００Ｈｚのサンプルレートで記録された。ＴＦＴは、Δｔ＝５０で実施された（これは、５０／８００＝０．０６２５ｓに達する）。これは、Ｍ＝４πΔｔ２＝３１，４１６のＴＦＴフィルタの幅を必要とする。しかしながら、ＴＦＴフィルタの係数は、Ｍ’＝１０・Δｔ＝５００で丸められたガウス分布により近似された（以下の「不確定性原理の導出」セクションを参照されたい）。周波数の関連する不確定度は、Δｆ＝１／（４πΔｔ）＝０．００１６である（又は、０．００１６ｘ８００＝１．２７Ｈｚ）。ＴＦＴは、８、１２、１６、２０、２４Ｈｚの課されたＦＯＴ周波数に対応して、ｆｍ＝０．０１、０．０１５、０．０２、０．０２５、０．０３に対して導出された。サンプルパワー及びクロススペクトルは、ＮＳ＝５００に対して（すなわち、５００／８００＝０．６２５ｓの時間間隔にわたって）計算された。循環仮定から独立しているデータだけが使われた。これは、η＝６．７２を持つサンプルスペクトルの近似のカイ二乗分布を導く。インピーダンスに対する信頼限界は、９０％の信頼レベル（Ｆα＝４．０５）に対して導出された。時間―周波数スペクトルは、最初の時系列と時間的に位置合わせされた。
【００８７】
例
図７は、通常の被験者に対する時間及び周波数の関数として、推定されたインピーダンス
の実数部分及び虚数部分を示す。青いラインは実数部分
を表し、赤いラインは虚数部分
を表す。グレイのラインは、図６Ｂに図示されたように、信頼限界を表す。８Ｈｚの周波数帯域の右側の軸は、インピーダンスの成分のスケールを示す。８Ｈｚで、実数部分
は正であり、呼吸サイクルの間で変動する一方、虚数部分
はほぼゼロに等しく、比較的一定のままである。信頼限界は、時間の経過で、
及び
の変化の統計的有意性について直観的な印象を与える。より高い周波数で、
は、
が徐々により正になっていくと共に、減少する。被験者が要求により吸い込むとき、信号は激しく妨げられ、これは信頼領域の大きい増加が伴われた。関連するコヒーレントなスペクトルは、下記のセクション「信頼限界の導出」に示される。
【００８８】
図８は、呼気の間、
の大きな負のスイングを持つＣＯＰＤを持つ患者に対するインピーダンスを示す（これは、明らかに、信頼限界からみて吸気値と著しく異なる）。斯様な負のスイングは、以前にＣＯＰＤにおいて記述されていて、呼気の間、気道の準最大の虚脱のため、流量制限に関係する。おそらく患者自身の呼吸の高周波成分のため、強制的振動を妨げるので、信頼領域は、吸気から呼気への転換点で最大であり、またこの逆も成立する。干渉性スペクトルのために下記のセクション「信頼限界の導出」を参照されたい。
【００８９】
説明
図９は、本発明による方法のステップを要約するフローチャートである。この例示的実施例では、圧力波は、被験者の気道と連通する患者インターフェイス装置４としてマウスピースに接続された拡声器を使用して生成される（ステップ１０１）。マウスピース４を通るガスの流れ及び圧力の振動が、測定のそれぞれの時系列を与えるために測定される（ステップ１０３）。定常状態において、呼吸インピーダンスは、静かに呼吸している被験者のマウスピース４の強制的圧力振動から推定できる。呼吸器系の機械的特性が吸気から呼気にしばしば変化するので、当該方法は一時的な「安定性」の状況の下でインピーダンスを推定する。主要な課題は、インピーダンスが周波数依存的な量であるということであり、そして、高い時間的解像度が必然的に低い周波解像度に至るということであり、これによって、推定を歪める。データが離散的時系列（マウスピース４の流量及び圧力のシーケンシャル測定）から成るので、これは短い時系列の最適時間―周波数分析を求める。
【００９０】
上記のセクション「方法」にて説明されたように、本発明は、時間が循環している（下記参照）とみなされる場合、短い時系列が変動しないという概念に基づく。説明されたように、時間及び周波数領域の全体の不確定度に下位限界を配置する、不確定性原理の新しいバージョンが、斯様な時系列に対して導出された。ステップ１０５では、この原理は、各時系列の相違を特定の時間及び周波数と関係している成分に区切るために用いられる（「離散的時間―周波数領域」又は「時間―周波数領域」への変換を通じて）。その後、時間―周波数領域の最小二乗法分析は、最適時間―周波数解像度を持つインピーダンスのバイアスされていない推定値を与える（ステップ１０７）。最後に、インピーダンスに対する信頼限界が、時間及び周波数の関数として作られる。主要なステップが、以下に要約され説明される。
【００９１】
線形モデル
呼吸インピーダンスの推定値は、単純な線形モデルに基づく（式３、図２）。装置２のマウスピース４の振動は、拡声器６により誘発される。流量及び圧力の結果として生じる変化は、患者の呼吸器系（その系の音響効果）の機械的特性に依存する。モデルにおいて、これは、拡声器入力を修正する２つの線形時間不変のフィルタ（図２の２及び３）により表される。流量及び圧力の結果として生じる変動は、流れ及び圧力両方のノイズの２つの独立源により妨げられる。呼吸インピーダンスは、２つの線形フィルタ（２及び３）により決定される。
【００９２】
ノイズの可能性がある源は、以下の通りである。１）強制的振動の周波数範囲内の患者の呼吸の成分。これは、おそらく最も重要な課題である（参考文献２を参照）。２）患者による吸い込み、咳嗽、又は他の動き。３）心臓性振動（特に気道抵抗が低いとき）。４）呼吸インピーダンス自体の時間依存性の変動。５）ランダムな測定エラー。
【００９３】
呼吸器系の非線形インタラクションが無視できないという指標があったが、線形モデルは、以前、同様の実験（参考文献２を参照）の間、流量と圧力との間の周波数依存的な関係を説明するために用いられた。しかしながら、ＦＯＴ装置２の流量と圧力との間の関係の線形記述は、圧力差が平均絶対圧と比較して小さいという事実によりサポートされる。当該関係を記述するためにインピーダンスを使用することは、これらの状況の下で期待できる定常波のインタラクションの簡略化である。しかしながら、主要な制限は、呼吸器系が安定と考えられる時間にある。
【００９４】
循環仮定
偶然のイベントを説明するために、流量及び圧力が所与の確率分布を持つランダム変数（ＲＶ）であると仮定された。よって、流量及び圧力のＮ個の対の測定値の短いシーケンスは、有限二変量の確率過程（確率プロセス）からのサンプルとして考えられる（Ｎ個の経時的に並べられた対のＲＶのセット）。斯様なプロセスにおいて、流量と圧力との間の周波数依存的な関係は、二変量の確率プロセスが（２次のオーダーで）安定していると仮定される場合、安定なインピーダンスによってのみ記述できる。これは、同時に起こる値と連続した値との間の共分散だけでなく、流量及び圧力の期待値が一定であることを意味する。しかしながら、連続した値の間の共分散は、開始及び終了効果のため、（これらの共分散がゼロに等しくない限り）有限シーケンスに対して一定でありえない。
【００９５】
しかしながら、この課題は、最後のイベントが最初のイベントに先行するという意味で、時間が循環すると仮定することにより解決される。このようにして、連続した値間の共分散は一定である一方、シーケンスは依然有限数のＲＶから成る。これは、シーケンスの離散的及び有限の性質の直接的な結果である。全体の時間（Ｎ＝１２個のイベントの経時的シーケンス）を見ることが唯一可能である時計を考える。あなたが１時（時刻ｔ１）及び２時（時刻ｔ２）を見る場合、ｔ１からｔ２まで１時間が経過したこともあり得るし、１３時間が経過した又は時刻ｔ２でのイベントがｔ１の１１時間前に起こったこともあり得る。他の知識なしで、ある時間と、ある時間に対して１２の整数倍数をプラス又はマイナスした時間との間を識別することは可能ではない。この不確定性は、「エイリアシング」に直接関係する（これは、斯様な時系列から連続信号を再構成しようとする場合に起こる）。しかし、測定されたシーケンスの最後の部分の偶然の変化によって、最初の部分に続かない変数両方の変動が生じる。これによってインピーダンスの推定におけるバイアスが生じる。しかしながら、これは、適用された時間―周波数分析（下記を参照）において除外される。
【００９６】
不確定性原理
「時間の不確定性（度）」は、一連の成分が時間的にどのくらい拡がるかを表す時系列の特性である。離散的、有限及び巡回の時系列に対して、これは、図３Ａに示されるような図を与える。時系列の成分は、時間の円の周りに規則的に配置される質点により表される（各質点は、成分の二乗の絶対値に等しい）。不確定度Δｔは、基準ポイント（ここでは、ｔ＝０のポイント、式１９）に関して定められ、重心に密接に関係する（式Ｂ３）。時系列が一つだけゼロ以外の成分（時刻ｔ＝０で）を持つ場合、このとき、重心は、ｔ＝０に位置され、Δｔはゼロである。全ての成分が等しい重さを持つ場合、このとき、重心は中心に位置され、Δｔは高い（時間領域の最大の拡がりと関連している）。
【００９７】
同様に、「周波数の不確定性（度）」Δｆは、周波数領域での拡がりを表す（式２２）。離散的フーリエ変換を通じて、Ｎ個の成分を持つ各時系列は、所与の振幅及び周波数（離散的周波数領域の表現）を持つＮ個の調和振動の和として書ける。離散的時系列に対して、これらの周波数も、循環的である。周波数領域の拡がりは、図３Ａのものと比較できるような円の周りにＮ個の質点として二乗振幅を配置することにより視覚化できる（各質点は特定の周波数に対応する）。
【００９８】
連続的な時間及び周波数に対して、以下の式の不確定性関係により表される時間及び周波数の不確定度のトレードオフがあることが示された。
（５４）
【００９９】
これは、（特定の条件の下）連続的で無限の時間に対して成り立つ。Δｔが小さい場合、Δｆは、積が少なくとも１／４πに等しくなければならないので大きくなければならず、またこの逆も成立する。この不等式は、粒子の位置及び運動量に対するハイゼンベルクの有名な不確定性原理に直接関係する。しかしながら、離散的時系列に対して、積ΔｔΔｆの下の限界は、ゼロである（下記の「不確定性原理の導出」を参照）。しかし、Δｔ及びΔｆの連合値の制限がある。上記の「方法」セクション（式２５）において、離散的で、有限及び循環の時系列に対して、以下の式が示される。
（５５）
ここで、σｍｉｎはＮの関数である。よってΔｔ及びΔｆ両方が非常に小さくなることができない。和Δｔ２＋Ｎ２Δｆ２は、離散的フーリエ変換が最初の時系列（図４Ａ）に正確に等しい時系列に対して最小である。Ｎが大きくなるにつれて、時間及び周波数の最小の連合の不確定度を持つこの時系列は、時間のガウス関数に近づく。この時系列は、量子力学における２つの引きつけられている粒子の最小の不確定性状態（いわゆる「量子調和発振器」）に類似している。フォーブズ等（参照文献４を参照）は、斯様な発振器に対する波動式から、同様の時系列を導出した。しかしながら、今の場合、この最小の不確定性状態は、離散的及び有限時系列に対するΔｔ及びΔｆの定義から簡単に導出される。下記の「不確定性原理の導出」において、式５５の不確定性原理が大きなＮに対するガボールの不確定性原理と、どのくらい関係するかが示される（図１４を参照）。
【０１００】
離散的時間―周波数領域
提示された時間―周波数変換（ＴＦＴ）は、離散的及び有限時系列を、Ｍ＜ＮであるＭ個の成分の短く続く振動のセットへ区切る（時間―周波数領域への変換）。各振動は、式５５の原理に従う両方の領域の最小の不確定度を持って、特定の時間及び周波数の周りに集まる。ＴＦＴは、ウィンドウ化された「短期フーリエ変換」又は「実行しているフーリエ変換」となる。結果が第１の振動が始まる時間から独立しているように、短く続く振動は、時間的に互いに最大限重なり合う。時系列の始まりの前に始まり、時系列の終わりの後に終わる振動だけが、循環仮定に依存している（これらは、バイアスを回避するために、分析において除外される）。以下の「時間―周波数変換の導出」セクションにおいて、変換が「正規直交である」（式３３）ことが示され、これは、離散的時間―周波数領域への変換後時系列の全ての情報が保存されることを意味する。何も得られていないし、何も失われていない。逆変換は、元の時系列を回復する。
【０１０１】
正確には、離散的時間―周波数領域と呼ぶことはできず、しかし、選ばれた値Ｍに対する時間―周波数領域だけについてである。この整数は各短期振動の継続時間を決定するが、これらの振動が集まる異なる周波数の数に等しい（ＴＦＴ共振周波数）。以下の「不確定性原理の導出」で示されるように、式５５の不確定性原理による最適のΔｔ及びΔｆはＭに依存する（式Ｂ８においてＮをＭにより置換する、図１１も参照）。もっともなことだが、Ｍが大きいと、対応する最適なΔｔが大きくなり、対応するΔｆが小さくなる。現在の分析では、拡声器により誘発される異なる振動を識別可能であるのにΔｆが十分小さいように、Ｍは選択される。
【０１０２】
拡声器への入力
例示的な実施例では、我々は拡声器６への入力として５つの調和振動のセットを使用した（それぞれ８、１２、１６、２０及び２４Ｈｚの周波数）。分析のために、Ｍの値（実際の有効値、式Ｂ６参照）は、不確定度が時間ユニットでΔｔ＝０．０６２５ｓ及びΔｆ＝１．２７Ｈｚで表されるように、選択された。ＴＦＴ係数は、これら５つの周波数に対してのみ計算され、これは５つの共振周波数を持つ帯域通過フィルタリングになる。１つの従来技術のアプローチは、拡声器入力のために単一の周波数を使用する。これは、時間―周波数分析のより少ない要件をもたらす（これらは、矩形の窓を持つ帯域通過フィルタを使用し、本発明のものより所与のΔｔに対してより大きいΔｆを持つ）。単一の周波数に対して、周波数解像度は、患者の呼吸のより高い高調波を除去するために依然有効であるが、あまり重要でない。他方では、インピーダンスの周波数依存は、呼吸メカニクス上に付加的情報を供給できる。代わりに、広帯域のノイズが適用でき、これは、検討された周波数で比較的小さなパワーを持って、入力信号のパワーが全ての周波数上に徐々に分散されるという不利な点を持つ（拡声器の同じ全体出力で、より低いＳＮ比）。更に他のアプローチでは、調和周波数での振動の間の相互作用を防止するために、換気された患者に非調和シヌソイドでできたベンチレータ波形が適用された。しかしながら、斯様なアプローチは、自由に呼吸する被験者に適用できない。調和周波数での呼吸器系の幾つかの干渉は、これは広帯域の刺激の後、系の周波数反応で観察されなかったが、実際に本発明で生じてもよい。
【０１０３】
時間―周波数領域の分散分析
ＴＦＴの正規直交特性は、確率過程の分散を時間―周波数依存の成分に区切ることを可能にし、これは時間―周波数パワースペクトルを与える（下記の「時間―周波数領域の分散の分析」を参照）。ＴＦＴサンプルパワーは、総サンプル分散に対して所与の時間及び周波数についての短く続く振動の寄与を表す。よって、ＴＦＴパワースペクトルの分析は、時間―周波数領域の分散の分析になる（式４２）。モデル（式４３）において、一緒に二変量の巡回静止確率過程を形成する、流量及び圧力のノイズの独立源があると仮定される。循環仮定に依存する値が放棄されるとき、時間―周波数領域の多くのＮＳ＝５００の連続した値が、このプロセスからのサンプルとしてみなされる（８００Ｈｚのサンプルレートで０．６２５ｓのエピソード）。各時点（1秒につき８００回）で、このタイプの新規な確率過程が以前の分散と必ずしも等しいわけではない分散で開始すると仮定された。よって、測定は５００の値の最大に重なり合うエピソードに再分割され、それぞれは異なる静止確率過程からのサンプルとしてみなされた。ノイズの源は「白色（ホワイト）である」とみなされなかったが、パワーはＴＦＴフィルタの各周波数バンドで一定であると仮定された。このとき、５００の値の各エピソードの間の平均パワーは、等価の自由度が導出された近似のカイ二乗分布に従う。
【０１０４】
時間―周波数領域の二変量の最小二乗法
ＮＳ＝５００の連続した値の各エピソードに対して、拡声器入力と流量（ｂＱ、ｔ、ｍ）との間の線形関係及び入力と圧力（ｂＰ、ｔ、ｍ）との間の関係の係数は、時間―周波数領域の単純な線形回帰により推定された（式４４、図５）。これは、時刻ｔ及び周波数ｆｍに関連したバイアスがされていない推定値に結果としてなる（それぞれ、
及び
式４６を参照）。各推定値は、時間―周波数領域のＮＳ値上の時間平均化されたパワー及びクロススペクトルに基づく。その後、インピーダンスは、比率
により推定された。
【０１０５】
Ｎの選択は、系が静止しているとみなされる時間に依存する。インピーダンスが各吸気及び呼気の中でほぼ安定であると仮定すると、呼吸のために、〜０．５ｓのエピソードがもっともらしいようである。示された例において、これは、よく機能している（図７及び図８）。しかしながら、ほとんどの従来技術において、スペクトルは、より長いエピソード、通常１０ｓより長いエピソードにわたって平均化されている。幾つかの平均化されたスペクトル値は、１６ｓの全体のエピソードにわたって、〜０．６５ｓの重なっていないセグメントから導出された。しかしながら、重なっていないセグメントの使用は、ここで使用される最大に重なり合うセグメント（すなわち、ＴＦＴで使用されるＭ値のセグメント）と比較して、同じ長さのエピソードに対して非常に少ない自由度を与える。大きな課題は、１６ｓのエピソードの使用が、通常ＣＯＰＤのような疾患においてより発音される吸気と呼気との間のあり得る生理的違いを平均化するということである。その場合、１６ｓにわたる流量と圧力との間の低い干渉性は、ノイズの反映だけでなく、（疾患の一部である）インピーダンスの呼吸内可変性の反映でもある。想定された静止の期間は、あまり短く選択されるべきではない。図８の例では、信頼区間は吸気からの呼気への転換点で比較的広く、逆も成り立つ。これは、患者の自身の呼吸の高周波成分に起因する。これらの変化は０．５ｓより明らかに短く、その結果、これらは推定のノイズとして現れる。従って、想定された静止の期間は、これらの変化をノイズとして「解釈する」のに十分長くなければならない。
【０１０６】
インピーダンス推定値
が２つの正常に分散されたＲＶ（
及び
）の比率として導出されるので、これは、期待値を持たないコーシー分布に従う。それで、厳密に言うと、
は、期待値と本当の値との間の差という意味において、バイアスを持たない。まだ、
及び
がバイアスされていない正常に分散されたＲＶであるので、
は、近似により期待値ｚｔ、ｍ＝ｂＰ，ｔ，ｍ／ｂＱ，ｔ，ｍで通常分散される。Ｄａｒｏｃｚｙ及びＨａｎｔｏｓ（参照文献２）は、拡声器入力から流量及び圧力それぞれまでの回帰係数の比率としてインピーダンス推定値を導出する類似のアプローチに従っている。患者及び装置の機械的特性の単純なモデルに基づいて、彼らは、ノイズが流量又は圧力の何れかだけで生じると仮定される場合、系統的誤差が導入されると主張した。斯様なバイアスに対する実験証拠は見つかった。しかしながら、他者は、拡声器が圧力波及び主要なノイズを生成し、（圧力のノイズがゼロである一方）患者自身の呼吸の高周波成分が、純粋な流量源であると仮定した。ＴＦＴが用いられた場合、推定値は式５２で
であるだろう。他方では、他の文献は、式５２の
と類似した推定値を使用した。
及び
は、圧力又は流量のノイズがゼロであるという仮定の下、最小二乗推定値であるので、これは、現実と一致していない場合バイアスを導き、このことは、数人の著者によりすでに指摘されていた（参照文献２を参照）。流量と圧力との間の干渉性が高くて、さもなければ推定値を放棄するよう勧められる場合、このタイプのエラーは最小であると議論された。これは、干渉性が統一の傾向がある場合、
及び
両方が
（よって、本当のｚｔ，ｍ）に近づくことが分かる式５３に従っている。式５０及び５１も、流量のノイズが、インピーダンスに対する信頼領域に影響を及ぼすことを示す（特に、歪度パラメータβ、図６Ｂを参照）。１つの以前の文献では、帯域通過フィルタリングされた圧力は、時間の関数（上記参照）として、インピーダンスを得るために、流れにより割られた。結果として生じるインピーダンスは、流量及び圧力のサンプル分散の等価な寄与を持つ「全体の最小二乗」推定値に、数値的に等しい。これらの値が、しばしば、（入力から流量までの干渉性が入力から圧力までの干渉性と等しい場合）二変量の最小二乗法で得られたものと同じ範囲にあることが期待できる。しかしながら、この「推定値」は、ランダム誤差に非常に影響されやすい平均値から導出されない。
【０１０７】
よってこれまで、推定されたインピーダンス上の信頼限界は、相対的に時間が長い間隔でのみ導出されてきた。本発明の実施例の時間―周波数に依存する信頼限界は、短く続くシーケンスが静止確率過程からのサンプルであるという仮説の有効性の継続的印象を提供する。これらは時間の経過でインピーダンスの変化の重要性を試験し、（図７の吸い込む間のような）信頼できない推定を自動的に拒絶するための根拠を提供することを可能にし、これはリアルタイムに監視する際の実際的な利点である。
【０１０８】
図１０は、更に詳細に本発明による装置２（及び特にコンピュータ１６により実施される処理ステップ）により実施される方法を例示するフローチャートである。図９のように、ステップ１２１で、圧力波は被験者により用いられているマウスピース４に接続される拡声器６を使用して生成され、マウスピース４を通る空気の流量及び圧力の振動が測定されて、測定のそれぞれの時系列を与えるためにデジタル化される（ステップ１２３）。このように、流量及び圧力は、離散的時間ｔ、ｑｔ及びｐｔの関数として与えられる。
【０１０９】
その後、ステップ１２５で、時系列は、離散的時間ｔ及び異なる周波数ｆｍ＝ｍ／Ｍの関数として、流量
及び圧力
を与えるために離散的時間―周波数領域へ各々変換される。ここで、ｍは周波数指標であり、Ｍは時間―周波数フィルタの有効幅である。周波数は、課されたＦＯＴ周波数の周波数に合うように、選択される。
【０１１０】
ステップ１２５は、流量時系列ｑｔに対する以下の式を評価するステップを有する。
（５６）
ここで、ｉ２＝−１であり、加重因子ｗｌは時間のガウス関数
により近似され、Δｔは時間における選ばれた不確定度である。フィルタは、Ｋ１＝５・Δｔで丸められる。加重因子ｗｌがガウス関数により近似されるとここで説明される一方、三角窓、区分的線形近似又は多項式関数のような他の窓又は加重因子が本発明により考察されることも理解されるべきである。
【０１１１】
時間及び周波数の関数として圧力
は、以下の式
（５６ａ）
から同様に導出される。
【０１１２】
その後、ステップ１２７において、パワー及びクロススペクトルが時間及び周波数の関数として導出される。流量のパワーは、以下の式により与えられる。
（５７）
ここで、ＮＳは時間―周波数領域のサンプルの数であり、ＮＳが奇数の場合Ｋ２＝１／２（ＮＳ−１）であり、｜｜は絶対値を表わす。
【０１１３】
ＰＰ，ｔ，ｍにより表わされる圧力のパワーは、以下の式から同様に導出される。
（５７ａ）
【０１１４】
時間―周波数領域の拡声器入力
から流量
までのクロススペクトルは、以下の式により与えられる。
（５８）
ここで、＊は複素共役を示す。
【０１１５】
拡声器入力から圧力ＣｘＰ，ｔ，ｍまでのクロススペクトルは、同様に
及び
から以下の式により導出される。
（５８ａ）
【０１１６】
ステップ１２９では、パワー及びクロススペクトルは、拡声器入力から流量及び圧力それぞれへの伝達関数を決定するために用いられる。
【０１１７】
特に、拡声器入力から流量への伝達関数は、クロススペクトルとパワーとの比により以下の式により与えられる。
（５９）
拡声器入力から圧力への伝達関数ＢｘＰ，ｔ，ｍは同様に以下の式により導出される。
（５９ａ）
【０１１８】
よって、ステップ１３１で、呼吸器系のインピーダンスは、実数成分ＲｒＳｔ，ｍ及び虚数成分ＸｒＳｔ，ｍを持つ伝達関数ブロックの比率から決定できる。
（６０）
【０１１９】
更にまた、呼吸インピーダンス上の信頼限界は、以下の通りに導出できる（ステップ１３３）。
【０１２０】
パワー及びクロススペクトルの自由度ηの等価数は、以下の式により与えられる。
（６１）
ここでＰはＮｘＮ行列であり、Ｈはエルミート転置行列を表わし、ｔｒ｛｝は行列のトレースを表す。行列Ｐは、以下の式のように規定される。
Ｐ＝Ｓ（Ｍｍ＋ＭＭ−ｍ）／ＮＳ （６２）
ここで、Ｓは主対角線上にＮＳ個の１（及び残りはゼロ）を持つ対角行列Ｓである。すなわち、
Ｓ＝ｄｉａｇ｛０，．．．，０，１，．．．，１，０，．．．，０｝ （６３）
【０１２１】
（シフトされる）時間―周波数変換行列Ｍｍは巡回行列であり、第１の行のエントリは、ｌ＝−Ｋ３，．．．，Ｋ３及びＫ３＝１／２（Ｎ−１）に対してｗｌｅ２πｉｍｌ／Ｍであり、ここで、Ｎは奇数である。Ｆαは、自由度２、η−２でのＦ分布の上位１００（１−α）％ポイントとしてラベル付けられる。
【０１２２】
と
との間の二乗コヒーレンスは、以下の式により与えられる。
（６４）
と
との間の二乗コヒーレンスＫｘＰ，ｔ，ｍ２は同様に以下の式から導出される。
（６４ａ）
【０１２３】
流量に関係する変数Ａｑ，ｔ，ｍ２は、以下の式として規定される。
（６５）
Ａｐ，ｔ，ｍ２で示される圧力に関係する類似の変数は、同様に以下の式で導出される。
（６５ａ）
【０１２４】
従って、呼吸インピーダンスの実数部分及び虚数部分に対する１００（１−α）％信頼は、それぞれ、Ｒｒｓ・ｃ１±ｃ２及びＸｒｓ・ｃ１±ｃ２により与えられ、ｃ１＝１＋Ａｑ，ｔ，ｍ・ｃ３、ｃ２＝｜Ｚｒｓｔ，ｍ｜・ｃ３、及び
である。
【０１２５】
所与の時間及び周波数に対する呼吸インピーダンスの推定値は、Ａｑ，ｔ，ｍ≧第１の閾値、又はＡｐ，ｔ，ｍ≧第２の閾値の何れかである場合、重要でないとして拒絶される。ある具体例では、第１の閾値及び第２の閾値は１に設定されるので、所与の時間及び周波数に対する呼吸インピーダンスの推定値は、Ａｑ，ｔ，ｍ≧１、又はＡｐ，ｔ，ｍ≧１の何れかである場合、重要でないとして拒絶される。
【０１２６】
最初と最後の効果を取扱うために、アルゴリズムは、幾つかの巡回データバッファを含む。
【０１２７】
測定装置
測定エアフローは、ニューモタコヘッド８及び内蔵型圧力トランスジューサ１０（例えばJaeger Masterscreen pneumotach type BF/IEC601-1, Hoechberg, Germany）で測定できる。マウスピース４での圧力は、差圧トランスジューサ１２（例えばHans Rudolph Pneumotach amplifier 1 series 1110, Shawnee, KS）で周囲空気を基準にして測定できる。圧力振動は、アンプ１８（例えばHarman Kardon HK970 amplifier, Washington DC）により増幅されるＡＤ変換カード１４（例えばNational Instruments PCI-6221, Dallas,TX）を通じたパーソナルコンピュータ１６（例えばHewlett Packard Compac dc 7600, PaloAlto, CA）からのアナログ出力信号により制御される拡声器６（例えばJaeger Masterscreen IOS, Hoechberg, Germany）によりＦＯＴ装置２内で生成できる。被験者は、ニューモタコヘッド８に接続されたメッシュワイヤ抵抗９を通じて呼吸できる。アナログ入力信号は、８００Ｈｚのサンプルレートで同じＡＤ変換カード１４を通じてデジタルシーケンスに変換でき、パーソナルコンピュータ１６に保存できる。出力信号の生成及びデータの分析は、適切なコンピュータソフトウェアを実行しているコンピュータ又はプロセッサ１６により実施できる。
【０１２８】
測定は、多くの呼吸サイクル、例えば９０秒の静かな呼吸をカバーする期間の間で、装置２を使用してなされる。
【０１２９】
不確定性原理の導出
時間及び周波数の不確定度（性）の幾何学的な解釈
式２０の規定によると、不確定度Δｔは、時系列ｘ≡｛ｘｔ｝の「重心」に直接関する。時系列のイベントは、値｜ｘｔ｜２を持つ質点として表わされ、原点を中心に、複素平面内の半径ｒｔを持つ円の周辺部に定間隔に配置される。各イベントは、複素平面内のポイントｒｔωｔで起こり、ここで、
及び
である。加重平均値
は、以下の式のように規定できる。
（Ｂ１）
【０１３０】
重心は
に位置される。行列式では、以下の式となる。
（Ｂ２）
ここで、Ω≡ｄｉａｇ｛１，ω，．．．，ωＮー１｝である。Ωがユニタリ行列であるという事実を使用して、二乗不確定度Δｔ２は、このとき、以下のように低減できる。
（Ｂ３）
【０１３１】
これは、図１２の基準点Ｒを通る軸（円の平面と直角をなす）についての第２の慣性モーメントに対する「平行軸定理」に従う。図１２は、時間の不確定度ΔｔとＮ＝１６に対するＣＨＣの固有ベクトルの重心との間の関係を例示する。固有ベクトルは、原点を中心として、半径ｒｔを持つ複素平面内の時間円に沿ってマップされる。各固有ベクトルの重心は、黒いドットにより表される。Ｃの対応する特異値は（基準点Ｒに近い）右側で重心に対して最小であり、図の左へ行くと徐々に増大する。各固有ベクトルに対する不確定度Δｔは、Ｒから重心を通る垂直線と円との交差点の１つまでの傾斜距離である。
【０１３２】
図１３の表現では、式１１から
であることが分かり、ここで、
は図１２のＡからＲまでの距離である。ピタゴラスの定理を使用して、
から、
であることが分かる。
【０１３３】
周波数領域において、不確定度Δｆは同様のやり方で幾何学的に解釈できる。
【０１３４】
行列Ｃの特異値
式２５の不確定性原理は、Ｃの二乗最小特異値及び二乗最大特異値であるσｍｉｎ２及びσｍａｘ２により決定される。特異値σは、
ｄｅｔ（ＣＨＣ−σ２Ｉ）＝０（Ｂ４）
から導出できる。ここで、ｄｅｔ（）は行列の決定要素を表わす。Ｎ＝２に対して、σｍｉｎ２＝２（２−√２）／π２及びσｍａｘ２＝２（２＋√２）／π２であることが容易に分かる。より高いＮに対して、Ｃの特異値（ＣＨＣの固有値）を導出するための様々な戦略が開発されてきた。図１１は、Ｎの関数として、σｍｉｎ及びσｍａｘを示す。より大きなＮ（例えば、Ｎ＞１５）に対して、σｍｉｎは
に近づき、σｍａｘはＮ√２／πに近づくことが分かる。時間円の半径に関して、このことは、
及び
を意味する。
【０１３５】
ｖを特異値σに対応するＣＨＣの単位固有ベクトルとする。正規方程式（ＣＨＣ−σ２Ｉ）ｖ＝０は、以下のように書ける。
（Ｂ５）
【０１３６】
第１の項は、実軸に沿った基準点Ｒからのｖの成分の出発に依存する（式Ｂ３を参照）。第２の項は、（離散的及び巡回の時間に対する）ｖの成分の２次の導関数とみることができる。式Ｂ５は、量子調和振動に対するシュレディンガー方程式に相当する。解は、離散的正規直交マシュー（Mathieu）関数と考えられる。σｍｉｎに対応する固有ベクトルｖｍｉｎは、時間の単一モードの関数である（図４Ａのような）。大きなＮに対して、ｖｍｉｎはベクトルｖ’ｍｉｎにより近似でき、その成分ν’ｍｉｎ，ｔは以下の式のような時間のガウス関数である。
（Ｂ６）
ここで、Ｎは、ｔ≧（１／２）Ｎの場合、ｔから減算されなければならない。ｔ＞５・Δｔの場合、成分ν’ｍｉｎ，ｔはゼロに近いので、ガウス関数は、最大値付近の１０・Δｔ個の最も大きい値だけを使用して、残りをゼロに設定することにより丸められる。結果として近似された固有ベクトルｖ’及び差ベクトルｄ≡ｖ−ｖ’をコールして、ノルム
は、近似値に作られる誤差の尺度として採ることができる。適切な数学的ソフトウェアを使用して、Ｎ＞７５の場合、相対的な誤差
が０．００５未満であることが分かる。
【０１３７】
ＣＨＣがＦと代わるので、これらの行列は、同じ固有ベクトルを共有する。Ｆ４＝Ｉであるので、Ｆの固有値は、１、−１、ｉ及び−ｉである。結果として、各ユニット固有ベクトルｖに対して、時間の不確定度は、周波数の不確定度に直接関係する。
（Ｂ７）
【０１３８】
これは、各固有ベクトルｖに対して、全体の不確定度が時間及び周波数にわたって均一に分散されることを意味する。
（Ｂ８）
よって、各固有ベクトルｖに対して、Δｔ＝ＮΔｆ＝σ／√２である。式Ｂ６によるｖｍｉｎのガウス近似は、以下のように書き直せる。
従って、全体の不確定度σｍｉｎ２は、このガウス関数の分散に等しい。式Ｂ８は、また、Ｎ＞１５に対して、最小及び最大の不確定度に対応するΔｔの値がそれぞれ、
及び
であることも意味する。後者は、２ｒｔが時間円内の最大可能なΔｔであるので、直観的に合理的である（図１３を参照）。図１３は、Ｎ＝１６に対する時間の不確定度Δｔと周波数の不確定度Δｆとの間の関係を例示する。全ての可能性があるＮ―ベクトルｘに対するアクセス可能な領域は、半径σｍｉｎ及びσｍａｘ（Ｃの最小及び最大特異値）を持つ２つの円により囲まれている。図１３において、ｒｔは時間円の半径であり、黒点はＣＨＣの固有ベクトルに属し、ｖｍｉｎ及びｖｍａｘはσｍｉｎ及びσｍａｘに対応する固有ベクトルであり、ｅ０及びｅ８はｔ＝０及びｔ＝８に対応する正準ベクトルであり、ｆ０及びｆ８は周波数ｆｎ＝０及びｆｎ＝８／１６＝１／２を持つ調和振動である。
【０１３９】
上述されたように、図１２は、Ｎ＝１６に対する時間円の全ての固有ベクトルｖに対する重心を示す。ＣＨＣは正規行列であるので、スペクトル定理に従ってＮ個の直交固有ベクトルを持つ。式Ｂ３からわかるように、対応するΔｔは、Ｒから複素平面の重心までの水平距離に依存する。Δｔは、図１２から、Ｒから重心を通る垂直線と円との交点までの弦の長さとして、読み取れる。Ｒに最も近い重心は、σｍｉｎと一致する。更に離れたところにある重心は、σ（及び
）の増大する値に対応する。比較的大きなＮに対して、ｖｍｉｎの重心とＲとの間の距離が式Ｂ３に従って１／４に近づき、そのΔｔ，ｍｉｎを見つけることは
に近づく。同じ固有値σ２＝Ｎ２／π２を持つ２つの直交固有ベクトルがあるという意味で、ＣＨＣの多数の固有値がある。Ｎが４により割り切れる場合のみ、斯様な多数の固有値が発生するだろう。これら２つの直交固有ベクトルの重心は両方とも虚数軸上に位置する。これらは、同じΔｔを持ち、
に等しい。対応するΔｆは
に等しい。これらの２つの固有ベクトルは、時間及び周波数領域において最大限拡がった（重心が時間及び周波数円両方の中心に位置する場合、得られるだろう値と同一である）。σのより高い値に対応する固有ベクトルは、（Ｒに対する）より大きな「不確定度」と関係するが、時間的により小さな拡がりを持つ。最大の総不確定度σｍａｘが固有ベクトルｖｍａｘに対して到達され、Δｔが２ｒｔにほぼ等しくなるまで、ベクトルの重みは時間及び周波数両方の円の反対側に実際にますます集中される。図１２の明らかな対称性に留意されたい（重心は、虚数軸に反映される）。
【０１４０】
図１３はΔｔの関数としてＮΔｆを示す。式２５の不確定性原理のため、アクセス可能な領域は、半径σｍｉｎ及びσｍａｘで原点の周りに中心がある２つの円により境界づけられる。これは、これらの円の間の全ての値が可能であるというわけではないが、円の外側の全ての値が不可能であることを意味する。ＣＨＣの固有ベクトルに帰属する座標は、全て同一ライン上に位置する（式Ｂ７のため）。幾つかの極端な場合も示される。１つは、正規ベクトル
である。この「重み」
は、ｔ＝０に完全に集結されているので、Δｔはゼロである。周波数領域への変換は
（Ｂ９）
を与える。ここで、ｆｎは式９として規定され、符号「＊」は複素共役を表わし、１は１だけを含むＮベクトルである。よって、ｅ０の重みは周波数領域の全ての周波数上に均一に分散され、
又は
である。他の極端な場合は、ｅ８である。（正規ベクトルｅｔを｛０，．．．０，１，０，．．．，０｝として規定し、ここで、１はゼロから始めてｔ番目の位置に表わされる。）ｅ８の重みは、複素平面内のポイント（−ｒｔ，０）に完全に集中されているのでΔｔは最大であり、２ｒｔに等しい。周波数領域において、
（Ｂ１０）
【０１４１】
列ベクトルｆ８＊は周波数ｆｎ＝８／１６＝１／２を持つ調和振動を記述する。その二乗成分も全ての周波数上に均一に分散されるので、
である。時間領域から始めると、他の極端な場合は、その重みが時間領域内に均等に分散され、ｆｎ＝０で周波数領域において鋭く集中する（
及びΔｆ＝０）振動ｆ０と、その重みが時間領域内に均等に分散され、ｆｎ＝１／２で周波数領域において鋭く集中する（
及びＮΔｆ＝２ｒｔ）ｆ８とである。Δｔ及びＮΔｆに対する可能性がある値が範囲［０，２ｒｔ］に限定されるので、全ての可能性がある組合せ（Δｔ，ＮΔｆ）がｖｍｉｎ，ｅ０，ｆ８，ｖｍａｘ，ｅ８及びｆ０に対する座標により囲まれる図１３の領域に限定されることが期待できる。これらの極端なベクトルは、通常の規則の実例でもある：（Ｆによる前置乗算を通じて）時系列ベクトルｘのフーリエ変換をすることは、図１３のプロット線の（Δｔ，ＮΔｆ）座標を交換する。別の表現では：ｘ及びＦｘに対する座標が同一のラインに反映される。これは、式２１及び式２３の規定を使用して、容易に検証される。
【０１４２】
ガボールの不確定性原理との比較
１９４６年のガボールにより説明された不確定性原理（下記の参照文献５を参照）は、連続時間及び周波数に対する積ΔｔΔｆに対する下位限界を決定する。
（Ｂ１１）
【０１４３】
この不等式は、式２５による離散的時間と周波数に対する不確定性原理にどのように関係するだろうか。式でＢ１１では、時間が無限のライン（又は、無限の半径を持つ円）上にマップされているが、Δｔ２は、また、ｔ＝０であるポイントに対して垂直な軸の周りの第２の慣性モーメントと規定される。式Ｂ１１の導出において、平均時間がゼロである（重心がｔ＝０に位置される）と仮定される。同じことは、周波数領域の拡がりに対しても満たされる（周波数が離散的でも周期的でもなく、連続的で無限であるという事実に関係なく）。
【０１４４】
式Ｂ１１と式２５との不等式の間の主要な違いは、ガボールの不確定性原理によると、ΔｔΔｆの下位限界ｔがΔｔの任意のゼロでない値を持つ（式Ｂ６と同じ形式の）ｔのガウス関数に対して達成される（参照文献５を参照）という事実にある。これは、一組の線形の独立（無限の）ベクトルである。しかしながら、式２５による総不確定度Δｔ２＋ＮΔｆ２の下位限界は、一つのΔｔの値（
に等しい）に対応するＣＨＣのたった一つの固有空間（ｖｍｉｎにわたる一次元複素部分空間）に対して達成される。他方では、この固有空間のベクトルは、ガボールの不確定性原理の下位限界を達成する（Ｎ→∞の限界状況において）。ＣＨＣの各固有ベクトルに対して、式Ｂ８から以下の式が成り立つ。
（Ｂ１２）
【０１４５】
ｖｍｉｎ（又は対応する固有空間の他の任意のベクトル）に対して、σは、Ｎが大きくなるにつれて、
に近づく（図１１）。このとき数の近似は、Ｎ→∞だと、以下の式を示す。
（Ｂ１３）
【０１４６】
しかしながら、有限のＮに対して、１／（４π）はΔｔΔｆに対する絶対の下位限界ではない。ｅ０及びｆ０に対して積はゼロであり、小さな乱数の追加は、わずかに異なるベクトルに対してゼロに近いことを示す。
【０１４７】
（Ｍイベントからなる、１＜Ｍ≦Ｎ）より短い時系列に対して得られたｖｍｉｎベクトルから導出されるＮ個のベクトルを考える場合、２つの不確定性原理の関係はより明らかになる。ｖｍｉｎが８ｘ８行列ＣＨＣから導出される図１４の例を参照されたい。Ｆｖｍｉｎの成分は、周波数ｆｎ（黒いドット）の関数として示される。これらの周波数は、１／８の倍数である。（ベクトルの単一モードの構造が円形の時間ベース上に完全なままであるように）８ベクトルｖｍｉｎの２つの最も小さな成分間に８個のゼロのブロックを挿入し、これは結果的に１６−ベクトルとなる。
（Ｂ１４）
【０１４８】
これは、Δｔの僅かな増大（０．７５８３から０．７９００へ）とほとんど重要でないΔｆの増大（０．０９４８から０．０９４９へ）を与える。図１４の白丸は、８個のゼロのブロックが時間領域の固有ベクトルの２つの最も小さな値の間に挿入されるとき周波数領域の補間された値を表し、実線は、時間領域の加えられたゼロの数が無限に向かうにつれて補間された値を表す。結果として生じる１６−ベクトルは、１／１６の倍数である周波数（白丸）で、周波数領域の補間を導く。これは、（時間領域の「ゼロ詰め」を通じた）周波数領域の良く知られた補間の形式である。
【０１４９】
このやり方で、所与のＮに対して、Ｎ−１の多くのベクトルは、ＭｘＭ行列ＣＨＣから得られるより小さなＭ次元ｖｍｉｎベクトルにゼロを加えることにより導出できる。自明の例は、Ｍ＝１に対して「１―ベクトル」ｖｍｉｎ＝１にゼロを加えることにより導出されるｅ０である。図１５は、Ｎ＝１６に対するＮΔｆ対Δｔのプロット線の対応する不確定度を示す。これらの値は、Δｔ＝０、
及び
のラインにより区切られるゾーンＩに現れる。Ｍ＝１に対する値はｅ０と一致し、Ｍ＝２に対する値は矢印により示され、より高いＭに対する値は（Ｎ＝１６に対するｖｍｉｎに一致する）識別ラインに位置されるＭ＝Ｎ＝１６に対する値に徐々に近づく。比較的大きなＭに対して、Δｔ及びΔｆがゼロの付加によりほとんど変化しないのに対し関連するσは
に近づくので、以下の式となる。
（Ｂ１５）
【０１５０】
結果として、比較的大きなＭに対する値は、図１５の双曲線ｘ・ｙ＝Ｎ／（４π）より僅かに下にプロットされる。時間が進みＮが大きくなるにつれ、値のトレイルはこの双曲線に沿って成長し、新しい値が識別ラインに現れる。ΔｆがΔｔに対してプロットされる場合、関係のある双曲線はｘ・ｙ＝１／（４π）により与えられ、ガボールの不確定性原理の下位限界で異なるガウス分布に対応する全ての可能性がある値を含む。
【０１５１】
正式の数学的証明なしで、Ｍ＜Ｎに対するｖｍｉｎから導出されるＮ−ベクトルは、図１５の（Δｔ，ＮΔｆ）平面のアクセス可能な領域までの更なる限界を持つ。これらのベクトルがＭがゼロでないイベントのシーケンスに対して最小の可能性がある総不確定度を持つので、図１５の黒いドットによりマークされる領域の外側の値が発生しないと期待できる。これは、時間領域の重みの変更なしで、周波数領域のベクトルを巡回的にシフトさせる（周波数円の反対側の部分）、Ω８とゾーンＩのベクトルを前置乗算することにより導出される、ゾーンＩＩの値に対しても成り立つ。よって、これらのベクトルのΔｆは、Δｔの変化なしで最大化される。これとは逆に、ゾーンＶＩのポイントは、時間領域のベクトルを巡回的にシフトさせる（これによりΔｆの変化なしでΔｔを最大化する）、Ｔ８と前置乗算された、ゾーンＩのベクトルから得られる。ゾーンＶのベクトルは、Ｔ８と前置乗算を通じてゾーンＩＩのベクトルから生じる。ゾーンＶＩＩＩのベクトルは、フーリエ変換をする（Ｆによる前置乗算は識別ラインの反映を意味する）ことにより、ゾーンＩのベクトルから生じる。Ω８によるこれらのベクトルの前置乗算は、ゾーンＩＩＩのベクトルを与え、Ｔ８による前置乗算は、ゾーンＶＩＩのベクトルを与え、Ω８による後者のベクトルの前置乗算は、ゾーンＩＶのベクトルを与える。（一つの時間周波数ゾーンから他の時間周波数ゾーンへ行かせる多くのやり方があることは明らかであり、例えば、ゾーンＶＩのベクトルのフーリエ変換を採用することはそれらをゾーンＩＩＩへ移動させる）。結論として、Δｔ、ＮΔｆの可能性がある組合せが、図１５の黒いドット（及び白丸）によりアウトラインが引かれたスペード状の領域に限定されることが期待できる。
【０１５２】
ゼロではない平均時間及び周波数に対する不確定性原理
式１９の基準時間ｔＲが、ゼロではないことが可能な場合、そのとき、以下の式が成り立つ。
（Ｂ１６）
【０１５３】
同様に、ゼロではない基準周波数
に対して、以下の式が成り立つ。
（Ｂ１７）
【０１５４】
及び
とする。
【０１５５】
ＣがｔＲ＝０及びｆＲ＝０に関して依然規定されるとき、ＣＲＨＣＲ及びＣＨＣは同様に以下の式になる。
（Ｂ１８）
【０１５６】
これは、任意の整数ｋ、ｎに対して有効であるＴｋΩｎ＝ω−ｎｋΩｎＴｋという関係を使用して容易に示される。それで、（σ２，ｖ）がＣＨＣの固有値−固有ベクトル対である場合、このとき、以下の式が成り立つ。
（Ｂ１９）
これは、
が固有値σ２を持つＣＲＨＣＲの固有ベクトルであることを意味する。時間領域のｔＲ及び周波数領域のｆＲに対する
の不確定度は、両方の領域のゼロに対するｖの不確定度に等しいことが成り立つ。これは、思いがけないことではない。（
による）周波数領域の巡回シフト及び（
による）時間領域の巡回シフトの後でも、
の成分の二乗絶対値（「重み」）のｖのものと同一である。
【０１５７】
時間―周波数変換の導出
時間―周波数変換行列Ｍの列は、正規直交である。
ＭＨＭ＝Ｉ（Ｃ１）
【０１５８】
これは、以下のように証明できる。積が巡回行列の和として以下のように書ける。
（Ｃ２）
【０１５９】
Ｔ−１によるｈｍＨの後置乗算は、その成分が右へ巡回的にシフトされることを意味するので、巡回Ｍｍは以下のように書ける。
（Ｃ３）
【０１６０】
巡回は正規行列なので、これらは、エルミート転置行列と置換される。
（Ｃ４）
ここで、
はｈｍの「複素自己相関」である。
【０１６１】
その成分は以下の通りである。
（Ｃ５）
【０１６２】
Ｎ―ベクトルｈｍの規定（式２６及び式２７）によるとこれはΩｍＮ／Ｍｈ０としても書ける。ＴｋΩｎ＝ω−ｎｋΩｎＴｋを使用して、ゼロ以外の成分は、下記の式に等しいことが成り立つ。
（Ｃ６）
ここで、
である。これらの成分がｍにわたって加算されるとき、下記の式となる。
（Ｃ７）
この式は、｛ｅ２πｉｍｕ／Ｍ｝がｍ＝０，．．．，Ｍ−１に対して等比数列である事実から成り立つ。ｈｍが単位ベクトルであるとみなされるので、
であることに留意されたい。結果として、式Ｃ２の和は、以下の式に書き直せる。
（Ｃ８）
これは、式Ｃ１の証明を完成する。
【０１６３】
正規直交性特性（式Ｃ１）は、ｈ０が、不確定度Δｔ及びＭ’＝１０・Δｔでゼロ以外の成分を持つ、丸められ正規化されたガウス分布により近似される場合（式Ｂ４を使用して）、ｈ０の選択に依存しないので（
と仮定する）、依然成り立つ。
【０１６４】
時間―周波数合成
正規直交性特性の直接的な推論は、所与の時系列ベクトルｘが短く後続する振動の加重和として書けるということである（式３４）。式２８の規定を使用して、以下の式が成り立つ。
（Ｃ９）
【０１６５】
時間―周波数領域の分散の分析
ＴＦＴの正規直交性特性の他の直接的な推論は、Ｎ―ベクトルｘの「エネルギー」（二乗されたノルム）が変換後保存される。
（Ｄ１）
【０１６６】
これは、時間―周波数領域の「分散分析」（ＡＮＯＶＡ）を可能にする。上記は、また、ｘのエネルギーがＭ周波数依存の成分に区切られることを意味する。
（Ｄ２）
【０１６７】
ランダムなＮ―ベクトルＸが、期待値Ｅ｛Ｘ｝＝０及び分散ｖａｒ｛Ｘ｝＝σ２Ｉを持つ多変量の正規分布を持つと仮定する。このとき、二次形式
は、Ｎ自由度のカイ二乗分布に従う。式Ｄ２によると、
は、Ｍ周波数依存の成分
に区切られる。二次形式
は、周波数ｆｍに対するＸの「ＴＦＴサンプルパワースペクトル」と呼ばれる。ホワイトノイズに対して、このＲＶの期待値は、以下の通りである。
ｔｒ｛・｝は行列のトレースを表わす。式Ｃ４を使用すると、ｔｒ｛ＭｍＨＭｍ｝＝ＮｈｍＨｈｍ＝Ｎが成り立つので、以下の式が成り立つ。
（Ｄ３）
【０１６８】
しかしながら、上記規定によるサンプルパワースペクトルは、循環仮定のためバイアスに影響される。
の成分は、Ｍ個のゼロ以外の要素を持つ線形時間不変のフィルタとの循環畳みこみによりＸから導出される（式２６−式３１）。このフィルタの構成は、
の最初のＭ−Ｋ−１個及び最後のＫ個の成分が循環仮定に依存するようになされ（これらの整数はゼロより大きいとする）、ここで、ＫはＭ／２以下の最大の整数である。このように、循環仮定によるバイアスがないパワー推定値は、以下の式で規定される「選択行列」Ｓで得られる。
（Ｄ４）
【０１６９】
Ｓによる
の前置乗算は、循環仮定とは独立している
の成分を選択し、これは、バイアスされていない二次形式
を与える。相補的周波数ｆＭ−ｍも含まれるとき、「全」バイアスがないサンプルパワーは、以下の式で規定される。
（Ｄ５）
【０１７０】
ｔ―ｆ領域の「サンプル」の数は、ＮＳ＝Ｎ−Ｍ＋１である。
【０１７１】
サンプルパワー（σ２により割られた）は、等価の自由度（ＥＤＯＦ）を持つ近似のカイ二乗分布に従う。
（Ｄ６）
【０１７２】
及び
に対して、
とすると、ゼロ平均ホワイトノイズＸに対して以下の式が成り立つ。
（Ｄ７）
これは以下の式を与える。
（Ｄ８）
【０１７３】
ここでは、ηは数学的コンピュータソフトウェアを使用して、式Ｄ８によって数値的に計算される。
【０１７４】
信頼限界の導出
及び
の平均及び分散
流量
のノイズがゼロ平均ホワイトノイズであるとみなされるので、式４３を使用して、以下の式が成り立つ。
（Ｅ１）
【０１７５】
の分散はσＱ，ｍ２Ｉである。積
が決定論的であるので、以下の式が成り立つ。
（Ｅ２）
【０１７６】
の平均及び分散は、式４６から以下の式が成り立つ。
（Ｅ３）
（Ｅ４）
【０１７７】
結果として、
は、ｂＱ，ｔ，ｍのバイアスされていない推定値である。
に対する平均及び分散は、類似の態様で、以下の式となる。
（Ｅ５）
【０１７８】
信頼限界
及び
の信頼限界は、最小二乗分析の標準的アプローチからわかる。式４３及び式４４の組み合わせは、以下の式を与える。
（Ｅ６）
【０１７９】
と
との間の直交性のため（図５を参照）、以下の式が成り立つ。
（Ｅ７）
【０１８０】
ランダム変数
は、
の形式を持ち、よって、（１つの自由度で）カイ二乗変数として（正確に）分散される。相補的周波数ｆＭ−ｍ＝１−ｍ／Ｍでの
も推定に含まれるとき、結果として生じるカイ二乗変数は２つの自由度を持つ。ノイズが全体のバイアスされていないパワースペクトルにより推定されるとき（式Ｄ５）、式Ｅ７は、ノイズの全パワーをＥＤＯＦη＝２＋（η−２）を持つ直交成分に分解させる。その後、式Ｅ７の右辺の２つの直交成分の比率は、（２，ηー２）程度の自由度でＦ分布に従う。この分布の上位１００（１−α）％のポイントをＦαにより示す。このとき、
（Ｅ８）
【０１８１】
ＴＦＴサンプルコヒーレンス（式４５）の規定を使用し、ＫＱ，ｔ，ｍ２＞０とすると、これは以下の式に書ける。
（Ｅ９）
【０１８２】
ＡＱ，ｔ，ｍ２によりこの不等の右辺を書き直すと、
（Ｅ１０）
【０１８３】
よって、
に対する１００（１−α）％信頼領域は、中心
及び半径
を持つ複素平面内の円により記述される。
に対する信頼領域は、同様の態様で導出される。
（Ｅ１１）
【０１８４】
に対する信頼領域に到達するために、
により式１０の両辺を乗算し、
を置換すると、以下の式を与える。
（Ｅ１２）
【０１８５】
複素変数の実数部分及び虚数部分を使用して、これは以下の式に再編成できる。
（Ｅ１３）
【０１８６】
次に、
により式Ｅ１１の両辺を割る。
により
を置換し（式４７）、Ｕにより
を置換すると、
（Ｅ１４）
【０１８７】
に対する信頼領域が、ここでツーステップアプローチで導出される。式Ｅ１３によると、Ｕは、１００（１−α）％の機会で、中心
及び半径
を持つ円内にある（図１６Ａ）。
及び
が相関しないと仮定する。このとき、所与のＵの値に対して、ｚｔ，ｍは、１００（１−α）％の機会で、中心Ｕ及び半径｜Ｕ｜ＡＰ，ｔ，ｍを持つ第２の円内にある（式Ｅ１４）。｜Ｕ｜に対する２つの極値、図１６ＡのＵｍｉｎ、Ｕｍａｘの大きさを考える。｜ｚｔ，ｍ｜に対する対応する極値は、図１６Ｂのｚｍｉｎ及びｚｍａｘの大きさである。式Ｅ１３及び式Ｅ１４から、
（Ｅ１５）
【０１８８】
に対する１００（１−α）％信頼領域の控え目な限界が、ここでｚｍｉｎ及びｚｍａｘを含む円により与えられ、その中心が
及び原点を通るライン上にある。結果として、対応する信頼領域は、以下の式により記述される。
（Ｅ１６）
【０１８９】
この領域は、中心
及び半径
を持つ図１６Ｂ及び図１６Ｃの円により区切られる。この円は、
の周りに同心円ではない。所与の
のｚｔ，ｍに対する可能な値の分布は、
に対して非対称である（図１６Ｄ）。この非対称性は、図１６Ｃのｔａｎβにより表される。式Ｅ１６から、以下の式が成り立つ。
（Ｅ１７）
【０１９０】
特別なケース１：流れにノイズがない。流れ内のノイズがゼロである、
である場合、そのとき、
及び
である（式４３及び式４６）。ｚｔ，ｍの対応する推定値を
により示す。
と仮定すると
である。
【０１９１】
この場合、
なので、
（Ｅ１８）
【０１９２】
これは、「通常の最小二乗推定値」である。このとき、ＫＱｘ，ｔ，ｍ２＝１及びＡＱ，ｔ，ｍ＝０である。ＡＱ，ｔ，ｍ＝０であるが、信頼限界は、式Ｅ１６により依然記述される。ここで、半径
を持つ信頼円は、
を中心に持つ。非対称パラメータｔａｎβはゼロに等しい。
【０１９３】
特別なケース２：圧力のノイズがない。圧力のノイズがゼロである、
である場合、そのとき、
及び
である。ｚｔ，ｍの対応する推定値が
により示されるとき、
（Ｅ１９）
【０１９４】
これは、「データ最小二乗法推定値」である。ここで、ＫＰｘ，ｔ，ｍ２＝１及びＡＰ，ｔ，ｍ＝０である。式Ｅ１６から、信頼円は、中心
及び半径
を持つ。非対称パラメータｔａｎβはＡＱ，ｔ，ｍに等しい。２つの推定値
及び
は、流量と圧力との間の二乗コヒーレンスを通じて関係がある。
（Ｅ２０）
【０１９５】
例：時間及び周波数の関数としての二乗コヒーレンス。図１７は、図７と同じ記録のための二乗コヒーレンスを示す。この例では、拡声器入力と流れとの間のコヒーレンスは、入力と圧力との間のコヒーレンスより概して低かった。図から明らかなように、コヒーレンスは、時間及び周波数両方の関数として変化する。最も小さなコヒーレンスは、吸気から呼気への転換点及びその逆の転換点で、最も低い周波数で得られる。飲み込むことは、特に入力と流れとの間のコヒーレンスへの深い影響を持つ。図１８は、ＣＯＰＤを持つ患者に対して図８に示される推定されたインピーダンスに対する時間及び周波数の関数としての二乗コヒーレンスを示す。呼吸フェーズ間の転換点のコヒーレンスの低下は、図１７の通常の例より非常に多く発音される。図１７及び図１８両方において、太線は、異なる共振周波数（８〜２４Ｈｚ）に対する時間の関数として、拡声器入力と流れとの間の二乗コヒーレンスＫＱｘ，ｔ，ｍ２を表し、細い線は、拡声器入力と圧力との間の二乗コヒーレンスＫＰｘ，ｔ，ｍ２を表す。「流量」は、エアフローの低周波成分を表し、「ｉｎｓｐ」は吸気を表し、「ｅｘｐ」は呼気を表す。
【０１９６】
結論として、最適時間―周波数解像度を持つ被験者の呼吸器系の伝達関数ブロックを推定するための方法が提示される。
【０１９７】
上述のように、呼吸インピーダンスの推定は、気道の妨害を評価するか、又は疾患の発病度を推定するための推定を使用する診断ツールによりなされる。従って、診断ツールは、呼吸インピーダンスに影響を及ぼすべき（薬理学的な又は他の）治療の効果を評価するためにも用いられる。
【０１９８】
例えば、呼吸インピーダンスの推定は、（ｉ）慢性閉塞性肺疾患（ＣＯＰＤ）の呼気の流量限界；（ｉｉ）発病度自体が時間の経過により吸入された気道拡張器の効果（例えば交感神経作用薬又は副交感神経作用薬）を評価するために用いられる（これは、調査設定に、更に、患者が標準強制的呼吸操作を実施できない臨床設定に適用できる）ＣＯＰＤ又は喘息の気道閉塞の当該発病度；（ｉｉｉ）喘息の気道閉塞又は睡眠の間のＣＯＰＤ；又は、（ｉｖ）睡眠無呼吸―呼吸低下症候群の疑いがある患者の睡眠中の上位の気道閉塞を検出するために使用できる。
【０１９９】
呼吸インピーダンスの推定は、また又は代わりに、内科的疾患の治療で使用される機械の設定を適応させるために用いられる。例えば、呼吸インピーダンスの推定が、ＣＯＰＤの非侵襲性ベンチレーションの設定の調整において用いられる。呼吸インピーダンスは、連続的に気道インピーダンスの重症度に関する情報を提供するためにベンチレータへの入力として用いられる（信頼値により示される信頼できない値は廃棄される）。この情報は、また、呼気の流量限界がちょうど克服される二相性陽気道圧のレベルを調整するためにも用いられる。他の例として、呼吸インピーダンスの推定は、援助として閉塞性睡眠時無呼吸症候群を持つ患者の持続気道陽圧のレベルをガイドするための補助として用いられる。
【０２００】
本発明は、図面及び前述の説明で例示され詳述されてきたが、斯様な図例及び説明は、図示的及び例示的であって、限定的でないとみなされるべきであり、本発明は、開示された実施例に限定されない。
【０２０１】
開示された実施例のバリエーションは、図面、明細書及び添付の請求の範囲の検討から、請求された本発明を実施する当業者により理解され遂行できる。請求項において、用語「有する」は他の要素又はステップを除外しないし、不定冠詞「ａ」又は「ａｎ」は複数を除外しない。単一のプロセッサ又は他のユニットが、請求項に引用される幾つかの部材の機能を成し遂げてもよい。特定の手段が相互に異なる従属請求項において引用されているという単なる事実は、これらの手段の組合せが効果的に使用できないことを示さない。コンピュータプログラムは、他のハードウェアと共に又は一部として光記憶媒体又はソリッドステート媒体のような適切な媒体に格納／配布されてもよいが、インターネット、他の有線又は無線通信システムを介してのような他の形式で配信されてもよい。請求項内の参照符号は、範囲を限定するものとして解釈されるべきではない。
【０２０２】
参照文献
１ DellacらによるExpiratory flow limitation detected by forced oscillation and negative expiratory pressure, European Respiratory Journal Vol.29 No.2, pages363-3742
２ Darｏczy B,Hantos Z.によるAn improved forced oscillation estimation of respiratory impedance. Int J Biomed Comput 13:221-235, 1982
３ ForbesGW,AlonsoMA.によるConsistent analogs of the Fourier uncertainty relation. Am J Physics 69:340-347, 2001
４ Forbes GW, Alonso MA, SiegmanAE.によるUncertainty relations and minimal uncertainty states for the discrete Fourier transform and the Fourier series. J Phys A: MathGen 36:7027-7047, 2003
５ GaborD.によるTheory of communications. J Inst Electr Eng93:429-457, 1946
【０２０３】
表１
シンボル及び省略形
時間及び周波数
ｆｍ、ｆｎ 離散的周波数ｍ／Ｍ（又はｎ／Ｎ）
ｆＲ Δｆに対する基準周波数
ｆｎ 周波数ｆｎを持つ調和振動を記述するＮ―ベクトル
ｍ ＴＦＴに対する周波数指標
Ｍ（Ｍ’） ＴＦＴフィルタの幅（又は切り詰めたガウス分布を使用した有効幅）
ｎ ＯＤＦＴに対する周波数指標
Ｎ 時系列数のイベントの数
ＮＳ 確率過程からのｔ―ｆ領域のサンプルの数
ｒｔ、ｒｆ 時間円（又は周波数円）の半径
ｔ 離散的時間指標
ｔＲ Δｔに対する基準時間
Ａ，Ｂ，Ｃ 時間及び周波数の不確定度に関係する行列
Ｆ ＯＤＦＴ行列
Ｍ(Ｍｍ) ＴＦＴ行列（又は周波数ｆｍに対する部分的なＴＦＴ行列）
Ｔ タイムシフトオペレータ；ｃｉｒｃ｛０，．．．，０，１｝
Δｔ、Δｆ 時間又は周波数の不確定度
ω 複素指数；ｅｘｐ（２πｉ／Ｎ）
σｍｉｎ、σｍａｘ Ｃの最小及び最大の単一の値
Ω 主対角線上にωｔを持つ対角行列；ｄｉａｇ｛ω０，．．．，ωＮ−１｝
変数の名前
ｐ 圧力
ｑ 流量
ｘ 拡声器入力（又は、未知の変数）
ｚ 呼吸インピーダンス
変数のタイプ 例
離散的時間の関数としての決定論的な変数 ｑｔ
離散的時間の関数としてのランダム変数 Ｑｔ
決定論的なＮ―ベクトル（小文字） ｑ
ランダムなＮ―ベクトル（大文字） Ｑ
ｔ―ｆ領域のベクトル（周波数ｆｍが中心）
ｔ―ｆ領域のベクトル（周波数ｆｍが中心で時刻ｔから開始）
ｔ―ｆ領域の線形回帰により「予測される」決定論的なベクトル
ｔ―ｆ領域の線形回帰に対するランダムな「エラー」
ｔ―ｆ領域のエラーの推定値
リニアフィルタ
ｈｑｘ ｘからｑまでのフィルタの期分けされたインパルス応答を持つベクトル
Ｈｑｘ，ｎ 周波数ｆｎに対するｘからｑまでの伝達関数
Ｃｑｘ 拡声器入力と流量との間のフィルタを記述する巡回行列
Λｘｑ 主対角線上にＨｑｘ，ｎを持つ対角行列
統計値
ＡＱ，ｔ，ｍ
に対する信頼領域への流量のノイズの寄与を決定するｔ―ｆ領域の変数
Ｆα （２、η−２）自由度でのＦ検定のための上位１００（１−α）％の信頼限界
ＫＱｘ，ｔ，ｍ２ ｔ―ｆ領域の流量と拡声器入力との間の二乗サンプルコヒーレンス
時間及び周波数の関数としてインピーダンスｚｔ，ｍの二変量の最小推定値
流量のノイズがゼロであるという仮定の下のｚｔ，ｍの推定値
α タイプＩＩエラー
β 信頼領域の中心から
の偏差を決定する角度
η 等価自由度
μＱ 平均の流量
σＱ２ 流量のノイズの分散
通常の数学的シンボル
ｉ 虚数（ｉ２＝−１）
Ｉ ＮｘＮの恒等行列
≡ 規定により等しいこと
｜ｘｔ｜ 複素数値ｘｔの大きさ
Ｎ―ベクトルｘの２−ノルム
ｘＨ ベクトルｘのエルミート転置
ｘ＊ ｘの複素共役
｛・｝ シーケンス（通常、列ベクトルとして表される）
ｄｉａｇ｛・｝ 主対角線上にシーケンスを持つ対角行列
ｃｉｒｃ｛・｝ 第１の行にシーケンスを巡回行列
省略形
ＦＯＴ 強制的振動技術
ｔ―ｆ領域 時間―周波数領域
ＴＦＴ 時間―周波数変換
ＯＤＦＴ 正規直交離散フーリエ変換
ＲＶ ランダム変数
【特許請求の範囲】
【請求項１】
患者インターフェイス装置及び被験者の気道を含む空気制御システムを作るため前記患者インターフェイス装置を被験者の気道に結合させるステップと、前記患者インターフェイス装置に動作的に結合される励振源を介して被験者の気道のガスの圧力、流量又はボリュームを振動するステップと、前記空気制御システムのガスの流量及び圧力を決定し、流量及び圧力を表すそれぞれの時系列を作るステップと、それぞれの時系列を時間―周波数領域に変換して変換時系列を作るステップと、変換された時系列から圧力及び流量のパワーを推定するステップと、変換された時系列に基づいて流量及び圧力のそれぞれのクロススペクトルを推定するステップと、推定されたパワー及び推定されたクロススペクトルから被験者の呼吸インピーダンスを推定するステップとを有する、呼吸インピーダンスを推定する方法。
【請求項２】
推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定するステップは、それぞれのパワー及びクロススペクトルから、圧力波から流量及び圧力へのそれぞれの伝達関数を決定するステップと、流量及び圧力伝達関数から呼吸インピーダンスを推定するステップとを有する、請求項１に記載の方法。
【請求項３】
推定されたインピーダンスの信頼限界を決定するステップを更に有する、請求項１又は２に記載の方法。
【請求項４】
流量及び圧力に対する測定の時系列は、ｑｔ及びｐｔでそれぞれ示され、それぞれの時系列を時間―周波数領域に変換するステップは、離散的時間ｔ及び異なる周波数ｆｍ＝ｍ／Ｍの関数として流量
及び圧力
を与えるために、流量の時系列に対して
と圧力の時系列に対して
とを評価するステップを有し、ここで、ｍは周波数指標であり、Ｍは時間―周波数フィルタの実効幅であり、ｉ２＝−１であり、加重ファクタｗｌは、時間の窓関数
により近似され、Δｔは時間の不確定度である、請求項１、２又は３に記載の方法。
【請求項５】
窓関数がガウス関数、三角関数、区分的線形近似関数又は多項式関数である、請求項４に記載の方法。
【請求項６】
流量及び圧力のパワーを推定するステップは、流量に対する
及び圧力に対する
を評価するステップを有し、ＮＳは時間―周波数領域のサンプルの数であり、Ｋ２＝１／２（ＮＳ−１）であり、ＮＳは奇数であり、｜・｜は絶対値を表わす、請求項４に記載の方法。
【請求項７】
生成された圧力波を持つ流量及び圧力のクロススペクトルを推定するステップが、流量に対する
と、圧力に対する
とを評価するステップを有し、＊は複素共役を示し、
は圧力波を表す、請求項６に記載の方法。
【請求項８】
推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定するステップは、
を評価することにより圧力波から流量への伝達関数と、
を評価することにより圧力波から圧力への伝達関数とを決定するステップを有し、ここで、呼吸インピーダンスは、実数成分ＲｒＳｔ，ｍ及び虚数成分ＸｒＳｔ，ｍを持つ
を評価することにより推定される、請求項７に記載の方法。
【請求項９】
を使用して圧力波と流量との間の二乗コヒーレンスを決定するステップと、
を使用して圧力波と圧力との間の二乗コヒーレンスを決定するステップと、
パワー及びクロススペクトルの自由度ηの数を決定するステップと、
から流量に関係する変数Ａｑ，ｔ，ｍ２を決定するステップと、
から圧力に関係する変数Ａｐ，ｔ，ｍ２を決定するステップと、Ｒｒｓ・ｃ１±ｃ２及びＸｒｓ・ｃ１±ｃ２それぞれから推定されたインピーダンスの実数部分及び虚数部分に対する１００（１−α）％の信頼限界を決定するステップであって、ここでｃ１＝１＋Ａｑ，ｔ，ｍ・ｃ３、ｃ２＝｜ＺｒＳｔ，ｍ｜・ｃ３、及び
であるステップとを含む、推定されたインピーダンス上の信頼限界を決定するステップを更に有する、請求項８に記載の方法。
【請求項１０】
Ａｑ，ｔ，ｍ≧第１の閾値又はＡｐ，ｔ，ｍ≧第２の閾値の何れかである場合、所与の時間及び周波数に対する呼吸インピーダンスの推定値を拒絶するステップを更に有する、請求項９に記載の方法。
【請求項１１】
第１の閾値が１であり、第２の閾値が１である、請求項１０に記載の方法。
【請求項１２】
変換された時系列に基づいて流量及び圧力のそれぞれのクロススペクトルを推定するステップが、入力としてガスの振動圧力、流量又はボリュームと出力として流量及び圧力とを持つモデルを使用して、最小二乗基準に基づいてなされる、請求項１に記載の方法。
【請求項１３】
前記空気制御システムのガスの流量を決定するステップが、前記空気制御システムに動作的に結合された流量センサを使用して流量を測定することにより達成され、前記空気制御システムのガスの圧力を決定するステップが、前記空気制御システムに動作的に結合された圧力センサを使用して達成される、請求項１に記載の方法。
【請求項１４】
患者インターフェイス装置と、被験者の気道のガスの圧力、流量又はボリュームの振動を生成するための励振源と、患者インターフェイス装置及び被験者の気道により規定される空気圧回路のガスの流量及び圧力を決定し、流量又は圧力を表しているそれぞれの時系列の値を出力するための手段と、プロセッサとを有し、前記プロセッサは、それぞれの時系列を時間―周波数領域へ変換し、それぞれ変換された時系列に基づいて、流量及び圧力のパワーを推定し、それぞれ変換された時系列に基づいて、流量及び圧力のそれぞれのクロススペクトルを推定し、推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定する、呼吸インピーダンスを推定するための装置。
【請求項１５】
前記空気圧回路のガスの流量及び圧力を決定するための手段が、前記空気制御システムに動作的に結合された流量センサと、前記空気制御システムに動作的に結合された圧力センサとを有する、請求項１４に記載の装置。
【請求項１６】
変換された時系列に基づいて流量及び圧力のそれぞれのクロススペクトルを推定することが、入力としてガスの振動する圧力、流量又はボリュームと、出力として流量及び圧力とを持つモデルを使用して、最小二乗基準に基づいてなされる、請求項１４に記載の装置。
【請求項１７】
前記プロセッサが、流量及び圧力それぞれのパワー及びクロススペクトルから、圧力波から流量及び圧力へのそれぞれの伝達関数を決定し、流量及び圧力の伝達関数から呼吸インピーダンスを推定することにより、推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定する、請求項１４に記載の装置。
【請求項１８】
前記プロセッサは、推定された呼吸インピーダンスに対する信頼限界を決定する、請求項１４又は１５に記載の装置。
【請求項１９】
適切なプロセッサ又はコンピュータによる実行の際、前記プロセッサ又はコンピュータが請求項１乃至１３の何れか一項に記載の方法を実施するコンピュータ可読媒体に組み込まれたコンピュータ可読のコードを具備する、コンピュータプログラム。
【請求項１】
患者インターフェイス装置及び被験者の気道を含む空気制御システムを作るため前記患者インターフェイス装置を被験者の気道に結合させるステップと、前記患者インターフェイス装置に動作的に結合される励振源を介して被験者の気道のガスの圧力、流量又はボリュームを振動するステップと、前記空気制御システムのガスの流量及び圧力を決定し、流量及び圧力を表すそれぞれの時系列を作るステップと、それぞれの時系列を時間―周波数領域に変換して変換時系列を作るステップと、変換された時系列から圧力及び流量のパワーを推定するステップと、変換された時系列に基づいて流量及び圧力のそれぞれのクロススペクトルを推定するステップと、推定されたパワー及び推定されたクロススペクトルから被験者の呼吸インピーダンスを推定するステップとを有する、呼吸インピーダンスを推定する方法。
【請求項２】
推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定するステップは、それぞれのパワー及びクロススペクトルから、圧力波から流量及び圧力へのそれぞれの伝達関数を決定するステップと、流量及び圧力伝達関数から呼吸インピーダンスを推定するステップとを有する、請求項１に記載の方法。
【請求項３】
推定されたインピーダンスの信頼限界を決定するステップを更に有する、請求項１又は２に記載の方法。
【請求項４】
流量及び圧力に対する測定の時系列は、ｑｔ及びｐｔでそれぞれ示され、それぞれの時系列を時間―周波数領域に変換するステップは、離散的時間ｔ及び異なる周波数ｆｍ＝ｍ／Ｍの関数として流量
及び圧力
を与えるために、流量の時系列に対して
と圧力の時系列に対して
とを評価するステップを有し、ここで、ｍは周波数指標であり、Ｍは時間―周波数フィルタの実効幅であり、ｉ２＝−１であり、加重ファクタｗｌは、時間の窓関数
により近似され、Δｔは時間の不確定度である、請求項１、２又は３に記載の方法。
【請求項５】
窓関数がガウス関数、三角関数、区分的線形近似関数又は多項式関数である、請求項４に記載の方法。
【請求項６】
流量及び圧力のパワーを推定するステップは、流量に対する
及び圧力に対する
を評価するステップを有し、ＮＳは時間―周波数領域のサンプルの数であり、Ｋ２＝１／２（ＮＳ−１）であり、ＮＳは奇数であり、｜・｜は絶対値を表わす、請求項４に記載の方法。
【請求項７】
生成された圧力波を持つ流量及び圧力のクロススペクトルを推定するステップが、流量に対する
と、圧力に対する
とを評価するステップを有し、＊は複素共役を示し、
は圧力波を表す、請求項６に記載の方法。
【請求項８】
推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定するステップは、
を評価することにより圧力波から流量への伝達関数と、
を評価することにより圧力波から圧力への伝達関数とを決定するステップを有し、ここで、呼吸インピーダンスは、実数成分ＲｒＳｔ，ｍ及び虚数成分ＸｒＳｔ，ｍを持つ
を評価することにより推定される、請求項７に記載の方法。
【請求項９】
を使用して圧力波と流量との間の二乗コヒーレンスを決定するステップと、
を使用して圧力波と圧力との間の二乗コヒーレンスを決定するステップと、
パワー及びクロススペクトルの自由度ηの数を決定するステップと、
から流量に関係する変数Ａｑ，ｔ，ｍ２を決定するステップと、
から圧力に関係する変数Ａｐ，ｔ，ｍ２を決定するステップと、Ｒｒｓ・ｃ１±ｃ２及びＸｒｓ・ｃ１±ｃ２それぞれから推定されたインピーダンスの実数部分及び虚数部分に対する１００（１−α）％の信頼限界を決定するステップであって、ここでｃ１＝１＋Ａｑ，ｔ，ｍ・ｃ３、ｃ２＝｜ＺｒＳｔ，ｍ｜・ｃ３、及び
であるステップとを含む、推定されたインピーダンス上の信頼限界を決定するステップを更に有する、請求項８に記載の方法。
【請求項１０】
Ａｑ，ｔ，ｍ≧第１の閾値又はＡｐ，ｔ，ｍ≧第２の閾値の何れかである場合、所与の時間及び周波数に対する呼吸インピーダンスの推定値を拒絶するステップを更に有する、請求項９に記載の方法。
【請求項１１】
第１の閾値が１であり、第２の閾値が１である、請求項１０に記載の方法。
【請求項１２】
変換された時系列に基づいて流量及び圧力のそれぞれのクロススペクトルを推定するステップが、入力としてガスの振動圧力、流量又はボリュームと出力として流量及び圧力とを持つモデルを使用して、最小二乗基準に基づいてなされる、請求項１に記載の方法。
【請求項１３】
前記空気制御システムのガスの流量を決定するステップが、前記空気制御システムに動作的に結合された流量センサを使用して流量を測定することにより達成され、前記空気制御システムのガスの圧力を決定するステップが、前記空気制御システムに動作的に結合された圧力センサを使用して達成される、請求項１に記載の方法。
【請求項１４】
患者インターフェイス装置と、被験者の気道のガスの圧力、流量又はボリュームの振動を生成するための励振源と、患者インターフェイス装置及び被験者の気道により規定される空気圧回路のガスの流量及び圧力を決定し、流量又は圧力を表しているそれぞれの時系列の値を出力するための手段と、プロセッサとを有し、前記プロセッサは、それぞれの時系列を時間―周波数領域へ変換し、それぞれ変換された時系列に基づいて、流量及び圧力のパワーを推定し、それぞれ変換された時系列に基づいて、流量及び圧力のそれぞれのクロススペクトルを推定し、推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定する、呼吸インピーダンスを推定するための装置。
【請求項１５】
前記空気圧回路のガスの流量及び圧力を決定するための手段が、前記空気制御システムに動作的に結合された流量センサと、前記空気制御システムに動作的に結合された圧力センサとを有する、請求項１４に記載の装置。
【請求項１６】
変換された時系列に基づいて流量及び圧力のそれぞれのクロススペクトルを推定することが、入力としてガスの振動する圧力、流量又はボリュームと、出力として流量及び圧力とを持つモデルを使用して、最小二乗基準に基づいてなされる、請求項１４に記載の装置。
【請求項１７】
前記プロセッサが、流量及び圧力それぞれのパワー及びクロススペクトルから、圧力波から流量及び圧力へのそれぞれの伝達関数を決定し、流量及び圧力の伝達関数から呼吸インピーダンスを推定することにより、推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定する、請求項１４に記載の装置。
【請求項１８】
前記プロセッサは、推定された呼吸インピーダンスに対する信頼限界を決定する、請求項１４又は１５に記載の装置。
【請求項１９】
適切なプロセッサ又はコンピュータによる実行の際、前記プロセッサ又はコンピュータが請求項１乃至１３の何れか一項に記載の方法を実施するコンピュータ可読媒体に組み込まれたコンピュータ可読のコードを具備する、コンピュータプログラム。
【図１】
【図２】
【図３Ａ】
【図３Ｂ】
【図４Ａ−４Ｂ】
【図５】
【図６Ａ】
【図６Ｂ】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６Ａ−１６Ｄ】
【図１７】
【図１８】
【図２】
【図３Ａ】
【図３Ｂ】
【図４Ａ−４Ｂ】
【図５】
【図６Ａ】
【図６Ｂ】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６Ａ−１６Ｄ】
【図１７】
【図１８】
【公表番号】特表２０１３−５１２７１８（Ｐ２０１３−５１２７１８Ａ）
【公表日】平成２５年４月１８日（２０１３．４．１８）
【国際特許分類】
【出願番号】特願２０１２−５４１６０７（Ｐ２０１２−５４１６０７）
【出願日】平成２２年１１月２４日（２０１０．１１．２４）
【国際出願番号】ＰＣＴ／ＩＢ２０１０／０５５３９２
【国際公開番号】ＷＯ２０１１／０６７６９８
【国際公開日】平成２３年６月９日（２０１１．６．９）
【出願人】（５９００００２４８）コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ (12,071)
【Ｆターム（参考）】
【公表日】平成２５年４月１８日（２０１３．４．１８）
【国際特許分類】
【出願日】平成２２年１１月２４日（２０１０．１１．２４）
【国際出願番号】ＰＣＴ／ＩＢ２０１０／０５５３９２
【国際公開番号】ＷＯ２０１１／０６７６９８
【国際公開日】平成２３年６月９日（２０１１．６．９）
【出願人】（５９００００２４８）コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ (12,071)
【Ｆターム（参考）】
[ Back to top ]