JP5868866B2 - 呼吸インピーダンスを推定するための装置の作動方法及び当該装置 - Google Patents
呼吸インピーダンスを推定するための装置の作動方法及び当該装置 Download PDFInfo
- Publication number
- JP5868866B2 JP5868866B2 JP2012541607A JP2012541607A JP5868866B2 JP 5868866 B2 JP5868866 B2 JP 5868866B2 JP 2012541607 A JP2012541607 A JP 2012541607A JP 2012541607 A JP2012541607 A JP 2012541607A JP 5868866 B2 JP5868866 B2 JP 5868866B2
- Authority
- JP
- Japan
- Prior art keywords
- pressure
- time
- flow
- frequency
- processor
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000000241 respiratory effect Effects 0.000 title claims description 65
- 238000000034 method Methods 0.000 title claims description 47
- 238000001228 spectrum Methods 0.000 claims description 45
- 230000005284 excitation Effects 0.000 claims description 21
- 238000012546 transfer Methods 0.000 claims description 18
- 238000005259 measurement Methods 0.000 claims description 7
- 230000001131 transforming effect Effects 0.000 claims 1
- 239000013598 vector Substances 0.000 description 77
- 230000006870 function Effects 0.000 description 47
- 239000011159 matrix material Substances 0.000 description 47
- 230000000875 corresponding effect Effects 0.000 description 29
- 125000004122 cyclic group Chemical group 0.000 description 24
- 230000010355 oscillation Effects 0.000 description 22
- 238000004458 analytical method Methods 0.000 description 21
- 238000009795 derivation Methods 0.000 description 19
- 230000001419 dependent effect Effects 0.000 description 16
- 230000005484 gravity Effects 0.000 description 15
- 230000029058 respiratory gaseous exchange Effects 0.000 description 15
- 208000006545 Chronic Obstructive Pulmonary Disease Diseases 0.000 description 12
- 238000013459 approach Methods 0.000 description 12
- 230000009466 transformation Effects 0.000 description 12
- 210000002345 respiratory system Anatomy 0.000 description 10
- 230000008569 process Effects 0.000 description 9
- 238000005309 stochastic process Methods 0.000 description 9
- 239000003570 air Substances 0.000 description 7
- 238000006243 chemical reaction Methods 0.000 description 7
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 7
- 230000000737 periodic effect Effects 0.000 description 7
- 238000010586 diagram Methods 0.000 description 6
- 201000010099 disease Diseases 0.000 description 6
- 230000000694 effects Effects 0.000 description 6
- 208000000884 Airway Obstruction Diseases 0.000 description 5
- 238000000540 analysis of variance Methods 0.000 description 5
- 230000008859 change Effects 0.000 description 5
- 238000012417 linear regression Methods 0.000 description 5
- 230000004962 physiological condition Effects 0.000 description 5
- 238000011282 treatment Methods 0.000 description 5
- 238000004891 communication Methods 0.000 description 4
- 230000000295 complement effect Effects 0.000 description 4
- 230000004044 response Effects 0.000 description 4
- 201000002859 sleep apnea Diseases 0.000 description 4
- 208000006673 asthma Diseases 0.000 description 3
- 230000015572 biosynthetic process Effects 0.000 description 3
- 230000003993 interaction Effects 0.000 description 3
- 210000004072 lung Anatomy 0.000 description 3
- 238000005192 partition Methods 0.000 description 3
- 230000002829 reductive effect Effects 0.000 description 3
- 238000003786 synthesis reaction Methods 0.000 description 3
- 241000238413 Octopus Species 0.000 description 2
- 241000529895 Stercorarius Species 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000001427 coherent effect Effects 0.000 description 2
- 238000004590 computer program Methods 0.000 description 2
- 230000001276 controlling effect Effects 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 229940079593 drug Drugs 0.000 description 2
- 239000003814 drug Substances 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 230000000144 pharmacologic effect Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 208000023504 respiratory system disease Diseases 0.000 description 2
- 230000002123 temporal effect Effects 0.000 description 2
- 210000000779 thoracic wall Anatomy 0.000 description 2
- 230000036962 time dependent Effects 0.000 description 2
- 230000001052 transient effect Effects 0.000 description 2
- 206010011224 Cough Diseases 0.000 description 1
- 238000001134 F-test Methods 0.000 description 1
- 102000005717 Myeloma Proteins Human genes 0.000 description 1
- 108010045503 Myeloma Proteins Proteins 0.000 description 1
- 206010041235 Snoring Diseases 0.000 description 1
- 239000000150 Sympathomimetic Substances 0.000 description 1
- 206010067775 Upper airway obstruction Diseases 0.000 description 1
- 239000012080 ambient air Substances 0.000 description 1
- 230000002051 biphasic effect Effects 0.000 description 1
- 238000009530 blood pressure measurement Methods 0.000 description 1
- 206010006451 bronchitis Diseases 0.000 description 1
- 239000000872 buffer Substances 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000000747 cardiac effect Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 208000035475 disorder Diseases 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000007274 generation of a signal involved in cell-cell signaling Effects 0.000 description 1
- 230000003434 inspiratory effect Effects 0.000 description 1
- 230000002045 lasting effect Effects 0.000 description 1
- 230000000670 limiting effect Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 208000001797 obstructive sleep apnea Diseases 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 239000000734 parasympathomimetic agent Substances 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 230000005610 quantum mechanics Effects 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 238000012883 sequential measurement Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 230000000638 stimulation Effects 0.000 description 1
- 230000002459 sustained effect Effects 0.000 description 1
- 230000009747 swallowing Effects 0.000 description 1
- 229940127230 sympathomimetic drug Drugs 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 238000009423 ventilation Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/08—Detecting, measuring or recording devices for evaluating the respiratory organs
- A61B5/085—Measuring impedance of respiratory organs or lung elasticity
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/7253—Details of waveform analysis characterised by using transforms
- A61B5/7257—Details of waveform analysis characterised by using transforms using Fourier transforms
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Pulmonology (AREA)
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Surgery (AREA)
- General Health & Medical Sciences (AREA)
- Biophysics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Physiology (AREA)
- Animal Behavior & Ethology (AREA)
- Pathology (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Mathematical Physics (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Psychiatry (AREA)
- Signal Processing (AREA)
- Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
- Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)
Description
p(f)=q(f)z(f)(1)
ここで、インピーダンスはz(f)により示される。この式は、z(f)が一定で、時間が無限である場合だけ成り立つ。システムが一時的に安定の場合、z(f)は短い時間間隔で依然推定される。しかしながら、これは、信号が時間及び周波数領域両方において鋭くは特定できないと言う不確定性原理により制限される。加えて、短い時間間隔の使用は、推定を偶然のイベントに益々依存させる。問題は、これが、呼吸サイクルの一部で安定なだけであるz(f)のような伝達関数の推定にどのくらい影響を及ぼすかである。
発明の詳細な説明のこのセクションは、本発明による方法及び装置の数学的基礎を説明しておく。本発明のより数学的でない説明は、以下の「説明」セクションにおいて提供される。
(2)
(3)
ここで、平均値μQ及びμPは定数である。結果として生じるQt及びPtが正常に分配されるので、仮定された二変量のプロセス{Qt、Pt}は、完全に変動しない(静止している)。
説明された静止プロセスが無限に長い一方、静止がN個のイベントの有限シーケンスに対してのみ仮定できている過渡現象を説明することが発明の目的である。これは、最後のイベントが最初のイベントに先行するという点で、時間が循環していると仮定することにより解決できる。結果として生じる時系列は、このとき、時計のような円に沿って定期的にNイベントを配置することにより表される。円から離れずに、同一方向で円を永久に追求できる。結果として生じる無限の時系列は、Nで周期的である。すなわち、任意の整数jに対して、以下の通りである。
xt=xt+jN (4)
(5)
結果として、
(6)
(7)
、
及び
とする。このとき、qp=Cqxxである。同様に、pP=Cpxxなので、式3は以下のようになる。
(8)
周波数依存的なインピーダンスを導くために、式8のモデルは、離散的周波数領域へ変換される必要がある。時間―周波数分析へのステップアップとして、これは、ここで短く要約される。所与の時間系列
から離散的周波数領域への変換は、時系列を離散的周波数
(n=0,...,N−1))を持つ一組の調和振動に区切る。斯様な振動は、以下のようなN―ベクトルにより記述される。
(9)
ここで、ω≡exp(2πint/N)及びi2=−1である。オイラーの定理により、
exp(2πint/N)=cos(2πnt/N)+isin(2πnt/N) (10)
よって、ベクトルfnは、実際に、周波数fnを持つ実軸及び虚軸に沿った複合振動を説明する。容易に以下の式が示される。
(11)
はNに等しい一方、異なる調和周波数の2つの振動は直交する。
(12)
この変換は、N―ベクトルFxを作り、これらの成分は、内積
である。式11から、Fの行(及び列)が正規直交であるので、行列はユニタリであることがわかる。
FFH=FHF=I (13)
ここでIは、NxNの単位行列である。これは、N振動からxの合成を可能にする。
(14)
(16)
(17)
に入れられ、列として対応する固有ベクトルをNxN行列FHへ入れられるとき、式15は以下の式になる。
(18)
の成分は、時間領域において鋭く局所的であるが、周波数領域においては不確かである。これとは逆に、フーリエ変換Fxの成分は、周波数領域において鋭く局所的であるが、時間領域においては不確かである。両方の領域で小さな不確定度を持つ信号を得るために、不確定さの尺度が、規定されなければならない。
に等しい。時刻tでの所与のイベントは、複素平面のポイントrtωtに位置される。フォーブズ及びアロンゾ(参照文献3)は、不確定度の幾つかの尺度を提案した。ここで、わずかに異なる尺度が使われ、これは、時間tでの各イベントと基準時間tRでのイベントとの間の二乗距離の加重平均として式19に規定される。
(19)
(20)
ここで、対角行列Ω≡diag{1,ω1,...,ωN−1}である。下記の式を規定することも可能である。
(21)
ここで、A≡Ω−Iである。物理的に、Δt 2は、基準ポイントR(時間tRに対応する)を通る(円の平面と直角をなす)軸の周りの質点のセットの第2の慣性モーメントである。図3Bに図示されるように、この慣性モーメントは、重心に直接関係する。幾何学的に、Δtは、Rから、重心を通る垂直線が円と交差する2つのポイントのうちの一方のポイントまでの距離に等しい。それで、Δtは、重心からRまでの水平距離により、完全に決定される。重心がRと一致する場合、xの重さは時間tRに完全に集中し、Δtは、ゼロである。重心が円の反対のポイントにある場合、重さはそのポイントに完全に集中し、Δtは最大であり、2rtに等しい。よって、図3Bにおいて、重心Bを通る垂直線は、ポイントCで円と交差する。不確定度Δtは距離
と等しい。
(22)
とすると、以下の式となる。
(23)
(24)
(25)
ODFTは、N個のイベントの時系列をN個の異なる周波数を持つ一組の振動へ区切られる。時間―周波数分析を実施する簡易なやり方は、短く続くシーケンスのM(M<N)個のイベントに対してODFTを繰り返し計算することである。従って、系列が周波数fm≡m/M(m=0,...,M−1)を持つM個の振動に区切られ、これは離散的「短時間フーリエ変換」になる。各短く続く振動がvminの成分で「窓」をつけられる(ウィンドウ化される)場合、時間及び周波数についての結果的な不確定度は最小である。
(26)
(27)
ここで、KはM/2以下の最も小さな整数である。この振動は、時間t=0及び周波数fm=m/Mを中央に置く。一般に、以下の式
(28)
の振動は、時間t及び周波数fmを中央に置く。
(29)
(30)
(31)
である。
(32)
を与える。バック変換は、再びMHによる乗算を意味し、
を与える。下記の「時間―周波数変換の導出」というセクションにおいて、Mの列は正規直交であることが示されるので、以下の式となる。
MHM=I(33)
(34)
励振源入力xとして、一組の調和振動を使用することが便利である。データを分析するために、このとき、これらの振動の周波数がTFT共振周波数fmと一致するように、TFTのスケールMが選択できる。最初に、周波数fn及び振幅amを持つ単一の正弦波入力を考える。オイラーの定理(式10)により、実数値の調和振動は、複素数値とその複素共役との和である。よって、励振源入力は、下記の式のように書ける。
(35)
ここで、星印は複素共役を示す。式9及び式10から、fn *=fN−nがわかる。NがMの整数倍である場合、周波数fn≡n/Mは、TFT共振周波数fm≡m/Mとマッチできる。その後、fn=fmであるように、mが選択される場合、以下の式となる。
(36)
(37)
mがゼロでない場合、MmμQがゼロであることが容易に示される。よって、
(38)
であることがわかる。fn=fmであるようにn及びmが選ばれる場合、このとき、式36及び式15を用いて、以下の式となる。
(39)
(40)
zm≡bP、m/bQ、m (41)
TFTは、時間―周波数領域の分散の分析を可能にする(「時間―周波数領域の分散分析」というこのセクションを参照)。例えば、流れのノイズの二乗ノルム(式8を参照)は以下の式に区切られる。
(42)
は、あたかもゼロ平均及び分散σQ、m 2Iを持つホワイトノイズから導出されるかのように、ほぼ分散される。結果として、
は、等価自由度ηを持つ近似のカイ二乗分布に従う。(相補的周波数fM−mの成分と同様に)循環仮定から独立している
のそれらの成分だけが含まれているηの導出のため、以下のセクション「時間―周波数領域の分散分析」を参照されたい。
各時間tで、図2のモデルにより記述されている、しばらくの間変動しないプロセスが始まると、ここで仮定される。当該プロセスは、N個のイベントを含み、循環的である(すなわち、Nで周期的)とみなされる。時間―周波数領域では、モデルは(式40から)以下の式となる。
(43)
ここで、ベクトルは時間tで始まる「短く続く」N個のベクトルである。定数bQ、m及びbP、mは、これらがプロセスが始まる時間に特有であるとみなされるので、同様にインデックスtを得た。従って、モデルへの入力と同様に定数は、(2つのプロセスが時間的に重複する場合であっても)異なる時間tで異なる。
及び
から導出できる。これは、最小二乗推定値
及び
を生じる。推定された関係は、以下の式である。
(44)
ここで、推定された「予測された」ベクトル
は、推定された「エラー」
に正規直交している(同様に、
は
と正規直交している)。図5では、ベクトルは矢印として示され、全ての変数は時間インデックスt及び周波数インデックスm(周波数fm=m/Mに対応する)を参照している。特に、図5では、励振源入力
は一定の変数である一方で、流量
及び圧力
はランダム変数である。ベクトル
及び
は、
によりスパンされた複素一次元の部分空間(ラインl)への
及び
の射影である。インピーダンス推定値は、
から分かる。エラーベクトル
及び
は
に垂直である。各ベクトルは、長さがノルムに等しい矢印により表される。ベクトル間の角度は、cos2φ及びcos2θが、
とそれぞれ
及び
との間の二乗コヒーレンスに等しいように描かれる。
との間の線形関係の程度は、以下の式のような「TFTは二乗サンプルコヒーレンス」の表現になる。
(45)
は、xからQまでの「TFTサンプルクロススペクトル」と呼ばれている。幾何学的に、KQx,t,m 2はcos2φと等しく、ここで、φは図5の
から
までの角度である。同様に、KPx,t,m 2はcos2θと等しい。最小二乗推定値は、標準的やり方で導出される。これは、以下の式を与える。
(46)
は、xからQまでのサンプルクロススペクトルとxのサンプルパワースペクトルとの比率である。
が同様のやり方で導出されるとき、インピーダンスは以下の式により推定される。
(47)
の実数部分及び虚数部分である。
(48)
を中心にした複素平面の円により区切られている。信頼領域は、以下の式により与えられる。
(49)
ここで、
であり、FαはF分布の上の100(1−α)%である(以下のセクション「信頼限界の導出」を参照)。励振源入力xtは、無次元変数であり、
は(図5から
なので)流量のユニットを持つ。循環仮定に依存する値は、ここで放棄される(サンプルの有効数は、NS=N−M+1である)。原点が円内にある場合(このとき、
は100・α%信頼レベルでゼロと著しくは異ならない)、関係は拒絶される。よって、AQ、t、m<1の場合、有意性が得られることが、式49から分かる。
に対する信頼領域は比較可能な円により区切られ、
と
との間の関係は、AP、t、m<1の場合、重要であると考えられる。
及び
が非相関であるという仮定の下、以下の「信頼限界の導出」セクションで導出される。当該領域は、以下の式により記述される。
(50)
の周りに同心円ではない。
の実数部分及び虚数部分が、
及び
により表わされる。垂直矢印は、実数部分に対する上位の信頼限界を示し、水平矢印は、インピーダンスの虚数部分に対する下位の信頼限界を示す。
の所与の実現のための「真の」zt,mの分布は、明らかに斜めに対称形である。図6Bでは、歪度は角度βにより決定される。歪度は、図5のβと角度φとの間の関係から明らかなように、完全に流量のノイズに起因する。
(51)
について対称形である。
(52)
及び
は、それぞれのサンプルパワー及びクロススペクトルである。式52の推定値と
との主要な違いは、推定されたインピーダンスの大きさにある。大きさは、以下の式を通じて関係付けられる。
(53)
が
の代わりに使用される場合過小評価される(
が使用される場合過大評価される)。
装置2の使用中、流量及び圧力が、800Hzのサンプルレートで記録された。TFTは、Δt=50で実施された(これは、50/800=0.0625sに達する)。これは、M=4πΔt 2=31,416のTFTフィルタの幅を必要とする。しかしながら、TFTフィルタの係数は、M’=10・Δt=500で丸められたガウス分布により近似された(以下の「不確定性原理の導出」セクションを参照されたい)。周波数の関連する不確定度は、Δf=1/(4πΔt)=0.0016である(又は、0.0016x800=1.27Hz)。TFTは、8、12、16、20、24Hzの課されたFOT周波数に対応して、fm=0.01、0.015、0.02、0.025、0.03に対して導出された。サンプルパワー及びクロススペクトルは、NS=500に対して(すなわち、500/800=0.625sの時間間隔にわたって)計算された。循環仮定から独立しているデータだけが使われた。これは、η=6.72を持つサンプルスペクトルの近似のカイ二乗分布を導く。インピーダンスに対する信頼限界は、90%の信頼レベル(Fα=4.05)に対して導出された。時間―周波数スペクトルは、最初の時系列と時間的に位置合わせされた。
図7は、通常の被験者に対する時間及び周波数の関数として、推定されたインピーダンス
の実数部分及び虚数部分を示す。青いラインは実数部分
を表し、赤いラインは虚数部分
を表す。グレイのラインは、図6Bに図示されたように、信頼限界を表す。8Hzの周波数帯域の右側の軸は、インピーダンスの成分のスケールを示す。8Hzで、実数部分
は正であり、呼吸サイクルの間で変動する一方、虚数部分
はほぼゼロに等しく、比較的一定のままである。信頼限界は、時間の経過で、
及び
の変化の統計的有意性について直観的な印象を与える。より高い周波数で、
は、
が徐々により正になっていくと共に、減少する。被験者が要求により吸い込むとき、信号は激しく妨げられ、これは信頼領域の大きい増加が伴われた。関連するコヒーレントなスペクトルは、下記のセクション「信頼限界の導出」に示される。
の大きな負のスイングを持つCOPDを持つ患者に対するインピーダンスを示す(これは、明らかに、信頼限界からみて吸気値と著しく異なる)。斯様な負のスイングは、以前にCOPDにおいて記述されていて、呼気の間、気道の準最大の虚脱のため、流量制限に関係する。おそらく患者自身の呼吸の高周波成分のため、強制的振動を妨げるので、信頼領域は、吸気から呼気への転換点で最大であり、またこの逆も成立する。干渉性スペクトルのために下記のセクション「信頼限界の導出」を参照されたい。
図9は、本発明による方法のステップを要約するフローチャートである。この例示的実施例では、圧力波は、被験者の気道と連通する患者インターフェイス装置4としてマウスピースに接続された拡声器を使用して生成される(ステップ101)。マウスピース4を通るガスの流れ及び圧力の振動が、測定のそれぞれの時系列を与えるために測定される(ステップ103)。定常状態において、呼吸インピーダンスは、静かに呼吸している被験者のマウスピース4の強制的圧力振動から推定できる。呼吸器系の機械的特性が吸気から呼気にしばしば変化するので、当該方法は一時的な「安定性」の状況の下でインピーダンスを推定する。主要な課題は、インピーダンスが周波数依存的な量であるということであり、そして、高い時間的解像度が必然的に低い周波解像度に至るということであり、これによって、推定を歪める。データが離散的時系列(マウスピース4の流量及び圧力のシーケンシャル測定)から成るので、これは短い時系列の最適時間―周波数分析を求める。
呼吸インピーダンスの推定値は、単純な線形モデルに基づく(式3、図2)。装置2のマウスピース4の振動は、拡声器6により誘発される。流量及び圧力の結果として生じる変化は、患者の呼吸器系(その系の音響効果)の機械的特性に依存する。モデルにおいて、これは、拡声器入力を修正する2つの線形時間不変のフィルタ(図2の2及び3)により表される。流量及び圧力の結果として生じる変動は、流れ及び圧力両方のノイズの2つの独立源により妨げられる。呼吸インピーダンスは、2つの線形フィルタ(2及び3)により決定される。
偶然のイベントを説明するために、流量及び圧力が所与の確率分布を持つランダム変数(RV)であると仮定された。よって、流量及び圧力のN個の対の測定値の短いシーケンスは、有限二変量の確率過程(確率プロセス)からのサンプルとして考えられる(N個の経時的に並べられた対のRVのセット)。斯様なプロセスにおいて、流量と圧力との間の周波数依存的な関係は、二変量の確率プロセスが(2次のオーダーで)安定していると仮定される場合、安定なインピーダンスによってのみ記述できる。これは、同時に起こる値と連続した値との間の共分散だけでなく、流量及び圧力の期待値が一定であることを意味する。しかしながら、連続した値の間の共分散は、開始及び終了効果のため、(これらの共分散がゼロに等しくない限り)有限シーケンスに対して一定でありえない。
「時間の不確定性(度)」は、一連の成分が時間的にどのくらい拡がるかを表す時系列の特性である。離散的、有限及び巡回の時系列に対して、これは、図3Aに示されるような図を与える。時系列の成分は、時間の円の周りに規則的に配置される質点により表される(各質点は、成分の二乗の絶対値に等しい)。不確定度Δtは、基準ポイント(ここでは、t=0のポイント、式19)に関して定められ、重心に密接に関係する(式B3)。時系列が一つだけゼロ以外の成分(時刻t=0で)を持つ場合、このとき、重心は、t=0に位置され、Δtはゼロである。全ての成分が等しい重さを持つ場合、このとき、重心は中心に位置され、Δtは高い(時間領域の最大の拡がりと関連している)。
(54)
(55)
ここで、σminはNの関数である。よってΔt及びΔf両方が非常に小さくなることができない。和Δt 2+N2Δf 2は、離散的フーリエ変換が最初の時系列(図4A)に正確に等しい時系列に対して最小である。Nが大きくなるにつれて、時間及び周波数の最小の連合の不確定度を持つこの時系列は、時間のガウス関数に近づく。この時系列は、量子力学における2つの引きつけられている粒子の最小の不確定性状態(いわゆる「量子調和発振器」)に類似している。フォーブズ等(参照文献4を参照)は、斯様な発振器に対する波動式から、同様の時系列を導出した。しかしながら、今の場合、この最小の不確定性状態は、離散的及び有限時系列に対するΔt及びΔfの定義から簡単に導出される。下記の「不確定性原理の導出」において、式55の不確定性原理が大きなNに対するガボールの不確定性原理と、どのくらい関係するかが示される(図14を参照)。
提示された時間―周波数変換(TFT)は、離散的及び有限時系列を、M<NであるM個の成分の短く続く振動のセットへ区切る(時間―周波数領域への変換)。各振動は、式55の原理に従う両方の領域の最小の不確定度を持って、特定の時間及び周波数の周りに集まる。TFTは、ウィンドウ化された「短期フーリエ変換」又は「実行しているフーリエ変換」となる。結果が第1の振動が始まる時間から独立しているように、短く続く振動は、時間的に互いに最大限重なり合う。時系列の始まりの前に始まり、時系列の終わりの後に終わる振動だけが、循環仮定に依存している(これらは、バイアスを回避するために、分析において除外される)。以下の「時間―周波数変換の導出」セクションにおいて、変換が「正規直交である」(式33)ことが示され、これは、離散的時間―周波数領域への変換後時系列の全ての情報が保存されることを意味する。何も得られていないし、何も失われていない。逆変換は、元の時系列を回復する。
例示的な実施例では、我々は拡声器6への入力として5つの調和振動のセットを使用した(それぞれ8、12、16、20及び24Hzの周波数)。分析のために、Mの値(実際の有効値、式B6参照)は、不確定度が時間ユニットでΔt=0.0625s及びΔf=1.27Hzで表されるように、選択された。TFT係数は、これら5つの周波数に対してのみ計算され、これは5つの共振周波数を持つ帯域通過フィルタリングになる。1つの従来技術のアプローチは、拡声器入力のために単一の周波数を使用する。これは、時間―周波数分析のより少ない要件をもたらす(これらは、矩形の窓を持つ帯域通過フィルタを使用し、本発明のものより所与のΔtに対してより大きいΔfを持つ)。単一の周波数に対して、周波数解像度は、患者の呼吸のより高い高調波を除去するために依然有効であるが、あまり重要でない。他方では、インピーダンスの周波数依存は、呼吸メカニクス上に付加的情報を供給できる。代わりに、広帯域のノイズが適用でき、これは、検討された周波数で比較的小さなパワーを持って、入力信号のパワーが全ての周波数上に徐々に分散されるという不利な点を持つ(拡声器の同じ全体出力で、より低いSN比)。更に他のアプローチでは、調和周波数での振動の間の相互作用を防止するために、換気された患者に非調和シヌソイドでできたベンチレータ波形が適用された。しかしながら、斯様なアプローチは、自由に呼吸する被験者に適用できない。調和周波数での呼吸器系の幾つかの干渉は、これは広帯域の刺激の後、系の周波数反応で観察されなかったが、実際に本発明で生じてもよい。
TFTの正規直交特性は、確率過程の分散を時間―周波数依存の成分に区切ることを可能にし、これは時間―周波数パワースペクトルを与える(下記の「時間―周波数領域の分散の分析」を参照)。TFTサンプルパワーは、総サンプル分散に対して所与の時間及び周波数についての短く続く振動の寄与を表す。よって、TFTパワースペクトルの分析は、時間―周波数領域の分散の分析になる(式42)。モデル(式43)において、一緒に二変量の巡回静止確率過程を形成する、流量及び圧力のノイズの独立源があると仮定される。循環仮定に依存する値が放棄されるとき、時間―周波数領域の多くのNS=500の連続した値が、このプロセスからのサンプルとしてみなされる(800Hzのサンプルレートで0.625sのエピソード)。各時点(1秒につき800回)で、このタイプの新規な確率過程が以前の分散と必ずしも等しいわけではない分散で開始すると仮定された。よって、測定は500の値の最大に重なり合うエピソードに再分割され、それぞれは異なる静止確率過程からのサンプルとしてみなされた。ノイズの源は「白色(ホワイト)である」とみなされなかったが、パワーはTFTフィルタの各周波数バンドで一定であると仮定された。このとき、500の値の各エピソードの間の平均パワーは、等価の自由度が導出された近似のカイ二乗分布に従う。
NS=500の連続した値の各エピソードに対して、拡声器入力と流量(bQ、t、m)との間の線形関係及び入力と圧力(bP、t、m)との間の関係の係数は、時間―周波数領域の単純な線形回帰により推定された(式44、図5)。これは、時刻t及び周波数fmに関連したバイアスがされていない推定値に結果としてなる(それぞれ、
及び
式46を参照)。各推定値は、時間―周波数領域のNS値上の時間平均化されたパワー及びクロススペクトルに基づく。その後、インピーダンスは、比率
により推定された。
が2つの正常に分散されたRV(
及び
)の比率として導出されるので、これは、期待値を持たないコーシー分布に従う。それで、厳密に言うと、
は、期待値と本当の値との間の差という意味において、バイアスを持たない。まだ、
及び
がバイアスされていない正常に分散されたRVであるので、
は、近似により期待値zt、m=bP,t,m/bQ,t,mで通常分散される。Daroczy及びHantos(参照文献2)は、拡声器入力から流量及び圧力それぞれまでの回帰係数の比率としてインピーダンス推定値を導出する類似のアプローチに従っている。患者及び装置の機械的特性の単純なモデルに基づいて、彼らは、ノイズが流量又は圧力の何れかだけで生じると仮定される場合、系統的誤差が導入されると主張した。斯様なバイアスに対する実験証拠は見つかった。しかしながら、他者は、拡声器が圧力波及び主要なノイズを生成し、(圧力のノイズがゼロである一方)患者自身の呼吸の高周波成分が、純粋な流量源であると仮定した。TFTが用いられた場合、推定値は式52で
であるだろう。他方では、他の文献は、式52の
と類似した推定値を使用した。
及び
は、圧力又は流量のノイズがゼロであるという仮定の下、最小二乗推定値であるので、これは、現実と一致していない場合バイアスを導き、このことは、数人の著者によりすでに指摘されていた(参照文献2を参照)。流量と圧力との間の干渉性が高くて、さもなければ推定値を放棄するよう勧められる場合、このタイプのエラーは最小であると議論された。これは、干渉性が統一の傾向がある場合、
及び
両方が
(よって、本当のzt,m)に近づくことが分かる式53に従っている。式50及び51も、流量のノイズが、インピーダンスに対する信頼領域に影響を及ぼすことを示す(特に、歪度パラメータβ、図6Bを参照)。1つの以前の文献では、帯域通過フィルタリングされた圧力は、時間の関数(上記参照)として、インピーダンスを得るために、流れにより割られた。結果として生じるインピーダンスは、流量及び圧力のサンプル分散の等価な寄与を持つ「全体の最小二乗」推定値に、数値的に等しい。これらの値が、しばしば、(入力から流量までの干渉性が入力から圧力までの干渉性と等しい場合)二変量の最小二乗法で得られたものと同じ範囲にあることが期待できる。しかしながら、この「推定値」は、ランダム誤差に非常に影響されやすい平均値から導出されない。
及び圧力
を与えるために離散的時間―周波数領域へ各々変換される。ここで、mは周波数指標であり、Mは時間―周波数フィルタの有効幅である。周波数は、課されたFOT周波数の周波数に合うように、選択される。
(56)
ここで、i2=−1であり、加重因子wlは時間のガウス関数
により近似され、Δtは時間における選ばれた不確定度である。フィルタは、K1=5・Δtで丸められる。加重因子wlがガウス関数により近似されるとここで説明される一方、三角窓、区分的線形近似又は多項式関数のような他の窓又は加重因子が本発明により考察されることも理解されるべきである。
は、以下の式
(56a)
から同様に導出される。
(57)
ここで、NSは時間―周波数領域のサンプルの数であり、NSが奇数の場合K2=1/2(NS−1)であり、||は絶対値を表わす。
(57a)
から流量
までのクロススペクトルは、以下の式により与えられる。
(58)
ここで、*は複素共役を示す。
及び
から以下の式により導出される。
(58a)
(59)
拡声器入力から圧力への伝達関数BxP,t,mは同様に以下の式により導出される。
(59a)
(60)
(61)
ここでPはNxN行列であり、Hはエルミート転置行列を表わし、tr{}は行列のトレースを表す。行列Pは、以下の式のように規定される。
P=S(Mm+MM−m)/NS (62)
ここで、Sは主対角線上にNS個の1(及び残りはゼロ)を持つ対角行列Sである。すなわち、
S=diag{0,...,0,1,...,1,0,...,0} (63)
との間の二乗コヒーレンスは、以下の式により与えられる。
(64)
と
との間の二乗コヒーレンスKxP,t,m 2は同様に以下の式から導出される。
(64a)
(65)
Ap,t,m 2で示される圧力に関係する類似の変数は、同様に以下の式で導出される。
(65a)
である。
測定エアフローは、ニューモタコヘッド8及び内蔵型圧力トランスジューサ10(例えばJaeger Masterscreen pneumotach type BF/IEC601-1, Hoechberg, Germany)で測定できる。マウスピース4での圧力は、差圧トランスジューサ12(例えばHans Rudolph Pneumotach amplifier 1 series 1110, Shawnee, KS)で周囲空気を基準にして測定できる。圧力振動は、アンプ18(例えばHarman Kardon HK970 amplifier, Washington DC)により増幅されるAD変換カード14(例えばNational Instruments PCI-6221, Dallas,TX)を通じたパーソナルコンピュータ16(例えばHewlett Packard Compac dc 7600, PaloAlto, CA)からのアナログ出力信号により制御される拡声器6(例えばJaeger Masterscreen IOS, Hoechberg, Germany)によりFOT装置2内で生成できる。被験者は、ニューモタコヘッド8に接続されたメッシュワイヤ抵抗9を通じて呼吸できる。アナログ入力信号は、800Hzのサンプルレートで同じAD変換カード14を通じてデジタルシーケンスに変換でき、パーソナルコンピュータ16に保存できる。出力信号の生成及びデータの分析は、適切なコンピュータソフトウェアを実行しているコンピュータ又はプロセッサ16により実施できる。
時間及び周波数の不確定度(性)の幾何学的な解釈
式20の規定によると、不確定度Δtは、時系列x≡{xt}の「重心」に直接関する。時系列のイベントは、値|xt|2を持つ質点として表わされ、原点を中心に、複素平面内の半径rtを持つ円の周辺部に定間隔に配置される。各イベントは、複素平面内のポイントrtωtで起こり、ここで、
及び
である。加重平均値
は、以下の式のように規定できる。
(B1)
に位置される。行列式では、以下の式となる。
(B2)
ここで、Ω≡diag{1,ω,...,ωNー1}である。Ωがユニタリ行列であるという事実を使用して、二乗不確定度Δt 2は、このとき、以下のように低減できる。
(B3)
であることが分かり、ここで、
は図12のAからRまでの距離である。ピタゴラスの定理を使用して、
から、
であることが分かる。
式25の不確定性原理は、Cの二乗最小特異値及び二乗最大特異値であるσmin 2及びσmax 2により決定される。特異値σは、
det(CHC−σ2I)=0(B4)
から導出できる。ここで、det()は行列の決定要素を表わす。N=2に対して、σmin 2=2(2−√2)/π2及びσmax 2=2(2+√2)/π2であることが容易に分かる。より高いNに対して、Cの特異値(CHCの固有値)を導出するための様々な戦略が開発されてきた。図11は、Nの関数として、σmin及びσmaxを示す。より大きなN(例えば、N>15)に対して、σminは
に近づき、σmaxはN√2/πに近づくことが分かる。時間円の半径に関して、このことは、
及び
を意味する。
(B5)
(B6)
ここで、Nは、t≧(1/2)Nの場合、tから減算されなければならない。t>5・Δtの場合、成分ν’min,tはゼロに近いので、ガウス関数は、最大値付近の10・Δt個の最も大きい値だけを使用して、残りをゼロに設定することにより丸められる。結果として近似された固有ベクトルv’及び差ベクトルd≡v−v’をコールして、ノルム
は、近似値に作られる誤差の尺度として採ることができる。適切な数学的ソフトウェアを使用して、N>75の場合、相対的な誤差
が0.005未満であることが分かる。
(B7)
(B8)
よって、各固有ベクトルvに対して、Δt=NΔf=σ/√2である。式B6によるvminのガウス近似は、以下のように書き直せる。
従って、全体の不確定度σmin 2は、このガウス関数の分散に等しい。式B8は、また、N>15に対して、最小及び最大の不確定度に対応するΔtの値がそれぞれ、
及び
であることも意味する。後者は、2rtが時間円内の最大可能なΔtであるので、直観的に合理的である(図13を参照)。図13は、N=16に対する時間の不確定度Δtと周波数の不確定度Δfとの間の関係を例示する。全ての可能性があるN―ベクトルxに対するアクセス可能な領域は、半径σmin及びσmax(Cの最小及び最大特異値)を持つ2つの円により囲まれている。図13において、rtは時間円の半径であり、黒点はCHCの固有ベクトルに属し、vmin及びvmaxはσmin及びσmaxに対応する固有ベクトルであり、e0及びe8はt=0及びt=8に対応する正準ベクトルであり、f0及びf8は周波数fn=0及びfn=8/16=1/2を持つ調和振動である。
)の増大する値に対応する。比較的大きなNに対して、vminの重心とRとの間の距離が式B3に従って1/4に近づき、そのΔt,minを見つけることは
に近づく。同じ固有値σ2=N2/π2を持つ2つの直交固有ベクトルがあるという意味で、CHCの多数の固有値がある。Nが4により割り切れる場合のみ、斯様な多数の固有値が発生するだろう。これら2つの直交固有ベクトルの重心は両方とも虚数軸上に位置する。これらは、同じΔtを持ち、
に等しい。対応するΔfは
に等しい。これらの2つの固有ベクトルは、時間及び周波数領域において最大限拡がった(重心が時間及び周波数円両方の中心に位置する場合、得られるだろう値と同一である)。σのより高い値に対応する固有ベクトルは、(Rに対する)より大きな「不確定度」と関係するが、時間的により小さな拡がりを持つ。最大の総不確定度σmaxが固有ベクトルvmaxに対して到達され、Δtが2rtにほぼ等しくなるまで、ベクトルの重みは時間及び周波数両方の円の反対側に実際にますます集中される。図12の明らかな対称性に留意されたい(重心は、虚数軸に反映される)。
である。この「重み」
は、t=0に完全に集結されているので、Δtはゼロである。周波数領域への変換は
(B9)
を与える。ここで、fnは式9として規定され、符号「*」は複素共役を表わし、1は1だけを含むNベクトルである。よって、e0の重みは周波数領域の全ての周波数上に均一に分散され、
又は
である。他の極端な場合は、e8である。(正規ベクトルetを{0,...0,1,0,...,0}として規定し、ここで、1はゼロから始めてt番目の位置に表わされる。)e8の重みは、複素平面内のポイント(−rt,0)に完全に集中されているのでΔtは最大であり、2rtに等しい。周波数領域において、
(B10)
である。時間領域から始めると、他の極端な場合は、その重みが時間領域内に均等に分散され、fn=0で周波数領域において鋭く集中する(
及びΔf=0)振動f0と、その重みが時間領域内に均等に分散され、fn=1/2で周波数領域において鋭く集中する(
及びNΔf=2rt)f8とである。Δt及びNΔfに対する可能性がある値が範囲[0,2rt]に限定されるので、全ての可能性がある組合せ(Δt,NΔf)がvmin,e0,f8,vmax,e8及びf0に対する座標により囲まれる図13の領域に限定されることが期待できる。これらの極端なベクトルは、通常の規則の実例でもある:(Fによる前置乗算を通じて)時系列ベクトルxのフーリエ変換をすることは、図13のプロット線の(Δt,NΔf)座標を交換する。別の表現では:x及びFxに対する座標が同一のラインに反映される。これは、式21及び式23の規定を使用して、容易に検証される。
1946年のガボールにより説明された不確定性原理(下記の参照文献5を参照)は、連続時間及び周波数に対する積ΔtΔfに対する下位限界を決定する。
(B11)
に等しい)に対応するCHCのたった一つの固有空間(vminにわたる一次元複素部分空間)に対して達成される。他方では、この固有空間のベクトルは、ガボールの不確定性原理の下位限界を達成する(N→∞の限界状況において)。CHCの各固有ベクトルに対して、式B8から以下の式が成り立つ。
(B12)
に近づく(図11)。このとき数の近似は、N→∞だと、以下の式を示す。
(B13)
(B14)
及び
のラインにより区切られるゾーンIに現れる。M=1に対する値はe0と一致し、M=2に対する値は矢印により示され、より高いMに対する値は(N=16に対するvminに一致する)識別ラインに位置されるM=N=16に対する値に徐々に近づく。比較的大きなMに対して、Δt及びΔfがゼロの付加によりほとんど変化しないのに対し関連するσは
に近づくので、以下の式となる。
(B15)
式19の基準時間tRが、ゼロではないことが可能な場合、そのとき、以下の式が成り立つ。
(B16)
に対して、以下の式が成り立つ。
(B17)
とする。
(B18)
(B19)
これは、
が固有値σ2を持つCR HCRの固有ベクトルであることを意味する。時間領域のtR及び周波数領域のfRに対する
の不確定度は、両方の領域のゼロに対するvの不確定度に等しいことが成り立つ。これは、思いがけないことではない。(
による)周波数領域の巡回シフト及び(
による)時間領域の巡回シフトの後でも、
の成分の二乗絶対値(「重み」)のvのものと同一である。
時間―周波数変換行列Mの列は、正規直交である。
MHM=I(C1)
(C2)
(C3)
(C4)
ここで、
はhmの「複素自己相関」である。
(C5)
(C6)
ここで、
である。これらの成分がmにわたって加算されるとき、下記の式となる。
(C7)
この式は、{e2πimu/M}がm=0,...,M−1に対して等比数列である事実から成り立つ。hmが単位ベクトルであるとみなされるので、
であることに留意されたい。結果として、式C2の和は、以下の式に書き直せる。
(C8)
これは、式C1の証明を完成する。
と仮定する)、依然成り立つ。
正規直交性特性の直接的な推論は、所与の時系列ベクトルxが短く後続する振動の加重和として書けるということである(式34)。式28の規定を使用して、以下の式が成り立つ。
(C9)
TFTの正規直交性特性の他の直接的な推論は、N―ベクトルxの「エネルギー」(二乗されたノルム)が変換後保存される。
(D1)
(D2)
は、N自由度のカイ二乗分布に従う。式D2によると、
は、M周波数依存の成分
に区切られる。二次形式
は、周波数fmに対するXの「TFTサンプルパワースペクトル」と呼ばれる。ホワイトノイズに対して、このRVの期待値は、以下の通りである。
tr{・}は行列のトレースを表わす。式C4を使用すると、tr{Mm HMm}=Nhm Hhm=Nが成り立つので、以下の式が成り立つ。
(D3)
の成分は、M個のゼロ以外の要素を持つ線形時間不変のフィルタとの循環畳みこみによりXから導出される(式26−式31)。このフィルタの構成は、
の最初のM−K−1個及び最後のK個の成分が循環仮定に依存するようになされ(これらの整数はゼロより大きいとする)、ここで、KはM/2以下の最大の整数である。このように、循環仮定によるバイアスがないパワー推定値は、以下の式で規定される「選択行列」Sで得られる。
(D4)
の前置乗算は、循環仮定とは独立している
の成分を選択し、これは、バイアスされていない二次形式
を与える。相補的周波数fM−mも含まれるとき、「全」バイアスがないサンプルパワーは、以下の式で規定される。
(D5)
(D6)
に対して、
とすると、ゼロ平均ホワイトノイズXに対して以下の式が成り立つ。
(D7)
これは以下の式を与える。
(D8)
及び
の平均及び分散
流量
のノイズがゼロ平均ホワイトノイズであるとみなされるので、式43を使用して、以下の式が成り立つ。
(E1)
が決定論的であるので、以下の式が成り立つ。
(E2)
(E3)
(E4)
は、bQ,t,mのバイアスされていない推定値である。
に対する平均及び分散は、類似の態様で、以下の式となる。
(E5)
及び
の信頼限界は、最小二乗分析の標準的アプローチからわかる。式43及び式44の組み合わせは、以下の式を与える。
(E6)
との間の直交性のため(図5を参照)、以下の式が成り立つ。
(E7)
は、
の形式を持ち、よって、(1つの自由度で)カイ二乗変数として(正確に)分散される。相補的周波数fM−m=1−m/Mでの
も推定に含まれるとき、結果として生じるカイ二乗変数は2つの自由度を持つ。ノイズが全体のバイアスされていないパワースペクトルにより推定されるとき(式D5)、式E7は、ノイズの全パワーをEDOFη=2+(η−2)を持つ直交成分に分解させる。その後、式E7の右辺の2つの直交成分の比率は、(2,ηー2)程度の自由度でF分布に従う。この分布の上位100(1−α)%のポイントをFαにより示す。このとき、
(E8)
(E9)
(E10)
に対する100(1−α)%信頼領域は、中心
及び半径
を持つ複素平面内の円により記述される。
に対する信頼領域は、同様の態様で導出される。
(E11)
により式10の両辺を乗算し、
を置換すると、以下の式を与える。
(E12)
(E13)
により式E11の両辺を割る。
により
を置換し(式47)、Uにより
を置換すると、
(E14)
及び半径
を持つ円内にある(図16A)。
及び
が相関しないと仮定する。このとき、所与のUの値に対して、zt,mは、100(1−α)%の機会で、中心U及び半径|U|AP,t,mを持つ第2の円内にある(式E14)。|U|に対する2つの極値、図16AのUmin、Umaxの大きさを考える。|zt,m|に対する対応する極値は、図16Bのzmin及びzmaxの大きさである。式E13及び式E14から、
(E15)
及び原点を通るライン上にある。結果として、対応する信頼領域は、以下の式により記述される。
(E16)
及び半径
を持つ図16B及び図16Cの円により区切られる。この円は、
の周りに同心円ではない。所与の
のzt,mに対する可能な値の分布は、
に対して非対称である(図16D)。この非対称性は、図16Cのtanβにより表される。式E16から、以下の式が成り立つ。
(E17)
である場合、そのとき、
及び
である(式43及び式46)。zt,mの対応する推定値を
により示す。
と仮定すると
である。
なので、
(E18)
を持つ信頼円は、
を中心に持つ。非対称パラメータtanβはゼロに等しい。
である場合、そのとき、
及び
である。zt,mの対応する推定値が
により示されるとき、
(E19)
及び半径
を持つ。非対称パラメータtanβはAQ,t,mに等しい。2つの推定値
及び
は、流量と圧力との間の二乗コヒーレンスを通じて関係がある。
(E20)
1 DellacらによるExpiratory flow limitation detected by forced oscillation and negative expiratory pressure, European Respiratory Journal Vol.29 No.2, pages363-3742
2 Daroczy B,Hantos Z.によるAn improved forced oscillation estimation of respiratory impedance. Int J Biomed Comput 13:221-235, 1982
3 ForbesGW,AlonsoMA.によるConsistent analogs of the Fourier uncertainty relation. Am J Physics 69:340-347, 2001
4 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
5 GaborD.によるTheory of communications. J Inst Electr Eng93:429-457, 1946
シンボル及び省略形
時間及び周波数
fm、fn 離散的周波数m/M(又はn/N)
fR Δfに対する基準周波数
fn 周波数fnを持つ調和振動を記述するN―ベクトル
m TFTに対する周波数指標
M(M’) TFTフィルタの幅(又は切り詰めたガウス分布を使用した有効幅)
n ODFTに対する周波数指標
N 時系列数のイベントの数
NS 確率過程からのt―f領域のサンプルの数
rt、rf 時間円(又は周波数円)の半径
t 離散的時間指標
tR Δtに対する基準時間
A,B,C 時間及び周波数の不確定度に関係する行列
F ODFT行列
M(Mm) TFT行列(又は周波数fmに対する部分的なTFT行列)
T タイムシフトオペレータ;circ{0,...,0,1}
Δt、Δf 時間又は周波数の不確定度
ω 複素指数;exp(2πi/N)
σmin、σmax Cの最小及び最大の単一の値
Ω 主対角線上にωtを持つ対角行列;diag{ω0,...,ωN−1}
変数の名前
p 圧力
q 流量
x 拡声器入力(又は、未知の変数)
z 呼吸インピーダンス
変数のタイプ 例
離散的時間の関数としての決定論的な変数 qt
離散的時間の関数としてのランダム変数 Qt
決定論的なN―ベクトル(小文字) q
ランダムなN―ベクトル(大文字) Q
t―f領域のベクトル(周波数fmが中心)
t―f領域のベクトル(周波数fmが中心で時刻tから開始)
t―f領域の線形回帰により「予測される」決定論的なベクトル
t―f領域の線形回帰に対するランダムな「エラー」
t―f領域のエラーの推定値
リニアフィルタ
hqx xからqまでのフィルタの期分けされたインパルス応答を持つベクトル
Hqx,n 周波数fnに対するxからqまでの伝達関数
Cqx 拡声器入力と流量との間のフィルタを記述する巡回行列
Λxq 主対角線上にHqx,nを持つ対角行列
統計値
AQ,t,m
に対する信頼領域への流量のノイズの寄与を決定するt―f領域の変数
Fα (2、η−2)自由度でのF検定のための上位100(1−α)%の信頼限界
KQx,t,m 2 t―f領域の流量と拡声器入力との間の二乗サンプルコヒーレンス
時間及び周波数の関数としてインピーダンスzt,mの二変量の最小推定値
流量のノイズがゼロであるという仮定の下のzt,mの推定値
α タイプIIエラー
β 信頼領域の中心から
の偏差を決定する角度
η 等価自由度
μQ 平均の流量
σQ 2 流量のノイズの分散
通常の数学的シンボル
i 虚数(i2=−1)
I NxNの恒等行列
≡ 規定により等しいこと
|xt| 複素数値xtの大きさ
N―ベクトルxの2−ノルム
xH ベクトルxのエルミート転置
x* xの複素共役
{・} シーケンス(通常、列ベクトルとして表される)
diag{・} 主対角線上にシーケンスを持つ対角行列
circ{・} 第1の行にシーケンスを巡回行列
省略形
FOT 強制的振動技術
t―f領域 時間―周波数領域
TFT 時間―周波数変換
ODFT 正規直交離散フーリエ変換
RV ランダム変数
Claims (18)
- 患者インターフェイス装置と、被験者の気道のガスの圧力、流量又はボリュームの振動を生成するための励振源と、患者インターフェイス装置及び被験者の気道により規定される空気圧回路のガスの流量及び圧力を決定し、流量又は圧力を表しているそれぞれの時系列の値を出力するための手段と、プロセッサとを有する呼吸インピーダンスを推定するための装置の作動方法であって、前記励振源により被験者の気道のガスの圧力、流量又はボリュームを振動させるステップと、前記手段により、前記患者インターフェイス装置及び被験者の気道を含む空気制御システムのガスの流量及び圧力を決定し、流量及び圧力を表すそれぞれの時系列を作るステップと、前記プロセッサにより、それぞれの時系列を時間―周波数領域に変換して変換時系列を作るステップと、前記プロセッサにより、時間及び周波数の関数として、変換された時系列から圧力及び流量のパワーを推定するステップと、前記プロセッサにより、時間及び周波数の関数として、変換された時系列に基づいて流量及び圧力のそれぞれのクロススペクトルを推定するステップと、前記プロセッサにより、推定されたパワー及び推定されたクロススペクトルから被験者の呼吸インピーダンスを推定するステップと、前記プロセッサにより、推定されたインピーダンスの信頼限界を決定するステップとを有する、装置の作動方法。
- 推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを前記プロセッサにより推定するステップは、それぞれのパワー及びクロススペクトルから、圧力波から流量及び圧力へのそれぞれの伝達関数を決定するステップと、流量及び圧力伝達関数から呼吸インピーダンスを推定するステップとを有する、請求項1に記載の装置の作動方法。
- 流量及び圧力に対する測定の時系列は、qt及びptでそれぞれ示され、前記プロセッサにより、それぞれの時系列を時間―周波数領域に変換するステップは、離散的時間t及び異なる周波数fm=m/Mの関数として流量
及び圧力
を与えるために、流量の時系列に対して
と圧力の時系列に対して
とを評価するステップを有し、ここで、mは周波数指標であり、Mは時間―周波数フィルタの実効幅であり、i2=−1であり、加重ファクタwlは、時間の窓関数
により近似され、Δtは時間の不確定度である、請求項1又は2に記載の装置の作動方法。 - 窓関数がガウス関数、三角関数、区分的線形近似関数又は多項式関数である、請求項3に記載の装置の作動方法。
- 流量及び圧力のパワーを前記プロセッサにより推定するステップは、流量に対する
及び圧力に対する
を評価するステップを有し、NSは時間―周波数領域のサンプルの数であり、K2=1/2(NS−1)であり、NSは奇数であり、|・|は絶対値を表わす、請求項3に記載の装置の作動方法。 - 生成された圧力波を持つ流量及び圧力のクロススペクトルを前記プロセッサにより推定するステップが、流量に対する
と、圧力に対する
とを評価するステップを有し、*は複素共役を示し、
は圧力波を表す、請求項5に記載の装置の作動方法。 - 推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを前記プロセッサにより推定するステップは、
を評価することにより圧力波から流量への伝達関数と、
を評価することにより圧力波から圧力への伝達関数とを決定するステップを有し、ここで、呼吸インピーダンスは、実数成分RrSt,m及び虚数成分XrSt,mを持つ
を評価することにより推定される、請求項6に記載の装置の作動方法。 -
を使用して圧力波と流量との間の二乗コヒーレンスを前記プロセッサにより決定するステップと、
を使用して圧力波と圧力との間の二乗コヒーレンスを前記プロセッサにより決定するステップと、
パワー及びクロススペクトルの自由度ηの数を前記プロセッサにより決定するステップと、
から流量に関係する変数Aq,t,m 2を決定するステップと、
から圧力に関係する変数Ap,t,m 2を前記プロセッサにより決定するステップと、Rrs・c1±c2及びXrs・c1±c2それぞれから推定されたインピーダンスの実数部分及び虚数部分に対する100(1−α)%の信頼限界を前記プロセッサにより決定するステップであって、ここでc1=1+Aq,t,m・c3、c2=|ZrSt,m|・c3、及び
であるステップとを含む、推定されたインピーダンス上の信頼限界を前記プロセッサにより決定するステップを更に有する、請求項7に記載の装置の作動方法。 - Aq,t,m≧第1の閾値又はAp,t,m≧第2の閾値の何れかである場合、前記プロセッサにより、所与の時間及び周波数に対する呼吸インピーダンスの推定値を拒絶するステップを更に有する、請求項8に記載の装置の作動方法。
- 第1の閾値が1であり、第2の閾値が1である、請求項9に記載の装置の作動方法。
- 選択された信頼レベル、流量の推定されたパワー及びクロススペクトルの自由度の数、並びに前記励振源と前記流量との間の二乗コヒーレンスを使用して、前記プロセッサにより、推定された呼吸インピーダンスに対する時間―周波数領域における流量に関係する変数を決定するステップと、
選択された信頼レベル、圧力の推定されたパワー及びクロススペクトルの自由度の数、並びに前記励振源と前記圧力との間の二乗コヒーレンスを使用して、前記プロセッサにより、推定された呼吸インピーダンスに対する時間―周波数領域における圧力に関係する変数を決定するステップと、
決定された流量に関係する変数が第1の閾値以上であるか又は決定された圧力に関係する変数が第2の閾値以上である場合、前記プロセッサにより、推定された呼吸インピーダンスを拒絶するステップとを更に有する、請求項1に記載の装置の作動方法。 - 変換された時系列に基づいて流量及び圧力のそれぞれのクロススペクトルを前記プロセッサにより推定するステップが、入力としてガスの振動圧力、流量又はボリュームと出力として流量及び圧力とを持つモデルを使用して、最小二乗基準に基づいてなされる、請求項1に記載の装置の作動方法。
- 前記空気制御システムのガスの流量を前記プロセッサにより決定するステップが、前記空気制御システムに動作的に結合された流量センサを使用して流量を測定することにより達成され、前記空気制御システムのガスの圧力を前記プロセッサにより決定するステップが、前記空気制御システムに動作的に結合された圧力センサを使用して達成される、請求項1に記載の装置の作動方法。
- 患者インターフェイス装置と、被験者の気道のガスの圧力、流量又はボリュームの振動を生成するための励振源と、患者インターフェイス装置及び被験者の気道により規定される空気圧回路のガスの流量及び圧力を決定し、流量又は圧力を表しているそれぞれの時系列の値を出力するための手段と、プロセッサとを有し、前記プロセッサは、それぞれの時系列を時間―周波数領域へ変換し、それぞれ変換された時系列に基づいて、時間及び周波数の関数として流量及び圧力のパワーを推定し、それぞれ変換された時系列に基づいて、時間及び周波数の関数として流量及び圧力のそれぞれのクロススペクトルを推定し、推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定し、推定された呼吸インピーダンスに対する信頼限界を決定する、呼吸インピーダンスを推定するための装置。
- 前記空気圧回路のガスの流量及び圧力を決定するための手段が、前記空気制御システムに動作的に結合された流量センサと、前記空気制御システムに動作的に結合された圧力センサとを有する、請求項14に記載の装置。
- 変換された時系列に基づいて流量及び圧力のそれぞれのクロススペクトルを推定することが、入力としてガスの振動する圧力、流量又はボリュームと、出力として流量及び圧力とを持つモデルを使用して、最小二乗基準に基づいてなされる、請求項14に記載の装置。
- 前記プロセッサが、流量及び圧力それぞれのパワー及びクロススペクトルから、圧力波から流量及び圧力へのそれぞれの伝達関数を決定し、流量及び圧力の伝達関数から呼吸インピーダンスを推定することにより、推定されたパワー及びクロススペクトルから被験者の呼吸インピーダンスを推定する、請求項14に記載の装置。
- 前記プロセッサが、
選択された信頼レベル、流量の推定されたパワー及びクロススペクトルの自由度の数、並びに前記励振源と前記流量との間の二乗コヒーレンスを使用して、推定された呼吸インピーダンスに対する時間―周波数領域における流量に関係する変数を決定し、
選択された信頼レベル、圧力の推定されたパワー及びクロススペクトルの自由度の数、並びに前記励振源と前記圧力との間の二乗コヒーレンスを使用して、推定された呼吸インピーダンスに対する時間―周波数領域における圧力に関係する変数を決定し、
決定された流量に関係する変数が第1の閾値以上であるか又は決定された圧力に関係する変数が第2の閾値以上である場合、推定された呼吸インピーダンスを拒絶する、請求項14に記載の装置。
Applications Claiming Priority (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
EP09177892 | 2009-12-03 | ||
EP09177892.8 | 2009-12-03 | ||
PCT/IB2010/055392 WO2011067698A2 (en) | 2009-12-03 | 2010-11-24 | Method and apparatus for estimating respiratory impedance |
Publications (3)
Publication Number | Publication Date |
---|---|
JP2013512718A JP2013512718A (ja) | 2013-04-18 |
JP2013512718A5 JP2013512718A5 (ja) | 2014-01-16 |
JP5868866B2 true JP5868866B2 (ja) | 2016-02-24 |
Family
ID=43877324
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2012541607A Active JP5868866B2 (ja) | 2009-12-03 | 2010-11-24 | 呼吸インピーダンスを推定するための装置の作動方法及び当該装置 |
Country Status (7)
Country | Link |
---|---|
US (1) | US9649050B2 (ja) |
EP (1) | EP2506765B1 (ja) |
JP (1) | JP5868866B2 (ja) |
CN (1) | CN102639056B (ja) |
BR (1) | BR112012013000A2 (ja) |
RU (1) | RU2606107C2 (ja) |
WO (1) | WO2011067698A2 (ja) |
Families Citing this family (15)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP6600460B2 (ja) * | 2011-12-30 | 2019-10-30 | コーニンクレッカ フィリップス エヌ ヴェ | 呼吸治療装置と一体化される気道インピーダンス測定 |
JP6031896B2 (ja) * | 2012-08-29 | 2016-11-24 | アイシン精機株式会社 | 呼吸信号推定装置及び呼吸信号推定方法 |
EP2934316B1 (en) * | 2012-12-19 | 2022-09-07 | Koninklijke Philips N.V. | Detection of respiratory disorders |
US10064583B2 (en) * | 2013-08-07 | 2018-09-04 | Covidien Lp | Detection of expiratory airflow limitation in ventilated patient |
US11175421B2 (en) | 2014-01-10 | 2021-11-16 | Cgg Services Sas | Device and method for mitigating cycle-skipping in full waveform inversion |
EP2945084A1 (en) * | 2014-05-12 | 2015-11-18 | Electrosalus Biyomedikal Sanayi ve Ticaret Anonim Sirketi | Auscultation data acquisition, communication and evaluation system incorporating mobile facilities |
WO2016004004A1 (en) * | 2014-07-01 | 2016-01-07 | Kosmo Technologies, Inc. | Methods and devices for positioning of a mandible of a subject for determining an optimal airway opening |
JP6859330B2 (ja) * | 2015-09-28 | 2021-04-14 | コーニンクレッカ フィリップス エヌ ヴェKoninklijke Philips N.V. | 漏れが存在する患者回路のコンプライアンスを推定する方法及びシステム |
ES2926626T3 (es) | 2016-02-03 | 2022-10-27 | Cognita Labs Llc | Prueba de función pulmonar basada en la técnica de oscilación forzada |
WO2019018938A1 (en) | 2017-07-26 | 2019-01-31 | Thorasys Thoracic Medical Systems Inc. | METHOD AND SYSTEM FOR ACQUIRING OSCILLOMETRY MEASUREMENTS |
US11433202B2 (en) * | 2017-12-18 | 2022-09-06 | Koninklijke Philips N.V. | Interactive guidance related to a subject's expiratory flow limitation results |
WO2020093176A1 (en) * | 2018-11-09 | 2020-05-14 | Thorasys Thoracic Medical Systems Inc. | Modular oscillometry device with dynamic calibration |
GB2583117B (en) * | 2019-04-17 | 2021-06-30 | Sonocent Ltd | Processing and visualising audio signals |
US11872344B2 (en) * | 2019-09-30 | 2024-01-16 | Koninklijke Philips N.V. | Detecting and treating COPD-OSA overlap syndrome |
CN116458872B (zh) * | 2023-06-13 | 2023-09-05 | 汶上县人民医院 | 一种呼吸数据的分析方法及系统 |
Family Cites Families (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
SU1635959A1 (ru) * | 1987-07-09 | 1991-03-23 | Петрозаводский государственный университет им.О.В.Куусинена | Способ определени нарушений бронхиальной проходимости |
US5318038A (en) * | 1993-01-19 | 1994-06-07 | Trustees Of Boston University | Infant respiratory impedance measuring apparatus and methods using forced oscillations |
US5555880A (en) * | 1994-01-31 | 1996-09-17 | Southwest Research Institute | High frequency oscillatory ventilator and respiratory measurement system |
US5885225A (en) * | 1996-01-23 | 1999-03-23 | Boys Town National Research Hospital | System and method for the measurement of evoked otoacoustic emissions |
US5795304A (en) * | 1996-03-27 | 1998-08-18 | Drexel University | System and method for analyzing electrogastrophic signal |
US6142952A (en) * | 1997-10-29 | 2000-11-07 | The Board Of Regents, The University Of Texas System | Method and apparatus for detection and diagnosis of airway obstruction degree |
US6066101A (en) | 1998-04-20 | 2000-05-23 | University Of Maryland | Airflow perturbation device and method for measuring respiratory resistance |
FR2791248B1 (fr) * | 1999-03-24 | 2001-08-24 | Georges Kehyayan | Dispositif d'analyse de bruits auscultatoires, en particulier de bruits respiratoires |
ITMI20021273A1 (it) * | 2002-06-11 | 2003-12-11 | Milano Politecnico | Sistema e metodo per la rilevazione automatica della limitazione del flusso espiratorio |
EP1596704B1 (en) | 2003-01-30 | 2019-08-07 | Compumedics Medical Innovation Pty Ltd | Algorithm for automatic positive air pressure titration |
NZ551074A (en) * | 2004-05-04 | 2010-08-27 | Univ Dalhousie | Method of assessment of airway variability in airway hyperresponsiveness |
JP5198162B2 (ja) | 2008-03-10 | 2013-05-15 | チェスト株式会社 | 呼吸インピーダンス測定装置及びその測定方法 |
-
2010
- 2010-11-24 WO PCT/IB2010/055392 patent/WO2011067698A2/en active Application Filing
- 2010-11-24 EP EP10800788.1A patent/EP2506765B1/en active Active
- 2010-11-24 JP JP2012541607A patent/JP5868866B2/ja active Active
- 2010-11-24 RU RU2012127567A patent/RU2606107C2/ru active
- 2010-11-24 CN CN201080054374.1A patent/CN102639056B/zh active Active
- 2010-11-24 BR BR112012013000A patent/BR112012013000A2/pt not_active Application Discontinuation
- 2010-11-24 US US13/511,207 patent/US9649050B2/en active Active
Also Published As
Publication number | Publication date |
---|---|
WO2011067698A2 (en) | 2011-06-09 |
EP2506765A2 (en) | 2012-10-10 |
CN102639056A (zh) | 2012-08-15 |
WO2011067698A3 (en) | 2011-09-01 |
RU2012127567A (ru) | 2014-01-20 |
CN102639056B (zh) | 2015-08-05 |
RU2606107C2 (ru) | 2017-01-10 |
BR112012013000A2 (pt) | 2018-08-28 |
JP2013512718A (ja) | 2013-04-18 |
US20120289852A1 (en) | 2012-11-15 |
US9649050B2 (en) | 2017-05-16 |
EP2506765B1 (en) | 2021-01-06 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP5868866B2 (ja) | 呼吸インピーダンスを推定するための装置の作動方法及び当該装置 | |
Charlton et al. | An assessment of algorithms to estimate respiratory rate from the electrocardiogram and photoplethysmogram | |
US20210298666A1 (en) | Methods and apparatus for monitoring chronic disease | |
KR102091167B1 (ko) | 수면 스테이지를 결정하는 시스템 및 방법 | |
EP3232912B1 (en) | System and method for cardiorespiratory sleep stage classification | |
CN113710151A (zh) | 用于检测呼吸障碍的方法和设备 | |
Sankar et al. | Performance study of various adaptive filter algorithms for noise cancellation in respiratory signals | |
US20080082018A1 (en) | Systems and methods for respiratory event detection | |
Moussavi | Fundamentals of respiratory sounds and analysis | |
Nguyen-Ky et al. | Measuring and reflecting depth of anesthesia using wavelet and power spectral density | |
Maes et al. | A fan-based, low-frequent, forced oscillation technique apparatus | |
US20220095952A1 (en) | Processor and method for determining a respiratory signal | |
JP2018000224A (ja) | 呼吸検出装置、呼吸検出方法および呼吸検出用プログラム | |
Yadollahi et al. | Respiratory Flow | |
JP2022507035A (ja) | 呼吸不全の検出および/または予測のための気道内圧の継続的管理のための方法および装置 | |
Bijaoui et al. | Mechanical output impedance of the lung determined from cardiogenic oscillations | |
WO2018002268A1 (en) | Improved methods for determining respiratory properties and system therefor | |
WO2021188786A1 (en) | Systems and methods for using a sound-producing breathing device to perform breathing exercises, and/or determine pulmonary function | |
Muthusamy et al. | An overview of respiratory airflow estimation techniques: Acoustic vs non-acoustic | |
JP7457102B2 (ja) | 呼吸数を処理するシステム | |
JP2021515620A (ja) | 高周波胸壁オシレータ | |
Reynolds et al. | Unrestrained acoustic plethysmograph for measuring specific airway resistance in mice | |
JP6784388B2 (ja) | 呼吸診断装置及び呼吸診断プログラム | |
Ionescu et al. | The Respiratory Impedance | |
Scott | The determination of respiratory impedance by the oscillatory airflow technique |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20131120 |
|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20131120 |
|
RD02 | Notification of acceptance of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7422 Effective date: 20131120 |
|
A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20141030 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20141208 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20150305 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20150909 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20151102 |
|
TRDD | Decision of grant or rejection written | ||
A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20151208 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20160106 |
|
R150 | Certificate of patent or registration of utility model |
Ref document number: 5868866 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |