ホーム > 特許ランキング > 国立大学法人京都大学
知財求人 - 知財ポータルサイト「IP Force」
▶ 国立大学法人京都大学の特許一覧 ▶ 株式会社モノプロダイムの特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2022074196
(43)【公開日】2022-05-18
(54)【発明の名称】呼吸曲線抽出方法と異常呼吸事象抽出方法
(51)【国際特許分類】
A61B 5/11 20060101AFI20220511BHJP
A61B 5/113 20060101ALI20220511BHJP
【ＦＩ】
A61B5/11 100
A61B5/113
【審査請求】未請求
【請求項の数】9
【出願形態】ＯＬ
(21)【出願番号】P 2020184035
(22)【出願日】2020-11-03
(71)【出願人】
【識別番号】504132272
【氏名又は名称】国立大学法人京都大学
(71)【出願人】
【識別番号】518060723
【氏名又は名称】株式会社モノプロダイム
(74)【代理人】
【識別番号】100183575
【弁理士】
【氏名又は名称】老田 政憲
(72)【発明者】
【氏名】堀口 由貴男
(72)【発明者】
【氏名】高橋 裕人
(72)【発明者】
【氏名】中西 弘明
(72)【発明者】
【氏名】椹木 哲夫
(72)【発明者】
【氏名】村瀬 亨
【テーマコード（参考）】
4C038
【Ｆターム（参考）】
4C038SS09
4C038ST04
4C038VA15
4C038VB01
4C038VB33
(57)【要約】
【課題】正確に呼吸曲線を抽出できる呼吸曲線抽出方法と無呼吸事象検出方法を提供すること。
【解決手段】２次元の人体の振動波形を測定し、上記測定した第１データを周波数領域データへ変換する変換工程と、上記変換工程で変換された第２データを、モルフォロジ成分分析をし、第２データを得るＧＭＣＡ工程と、上記ＧＭＣＡ工程で得られた第３データを、離散信号データへ変換する逆変換工程と、を含む呼吸曲線抽出法を用いる。
【選択図】 図１
【解決手段】２次元の人体の振動波形を測定し、上記測定した第１データを周波数領域データへ変換する変換工程と、上記変換工程で変換された第２データを、モルフォロジ成分分析をし、第２データを得るＧＭＣＡ工程と、上記ＧＭＣＡ工程で得られた第３データを、離散信号データへ変換する逆変換工程と、を含む呼吸曲線抽出法を用いる。
【選択図】 図１
【特許請求の範囲】
【請求項1】
２次元の人体の振動波形を測定し、前記測定した第１データを周波数領域データへ変換する変換工程と、
前記変換工程で変換された第２データを、モルフォロジ成分分析で、信号表現のスパース性を拠り所にした一般化モルフォロジ成分分析をし、第２データを得るＧＭＣＡ工程と、
前記ＧＭＣＡ工程で得られた第３データを、離散信号データへ変換する逆変換工程と、を含む呼吸曲線抽出法。
前記変換工程で変換された第２データを、モルフォロジ成分分析で、信号表現のスパース性を拠り所にした一般化モルフォロジ成分分析をし、第２データを得るＧＭＣＡ工程と、
前記ＧＭＣＡ工程で得られた第３データを、離散信号データへ変換する逆変換工程と、を含む呼吸曲線抽出法。
【請求項2】
前記一般化モルフォロジ成分分析は、
前記振動波形の複数の源信号に重みをつけて、足し合わせたものが、前記振動波形と一致するように，かつ，各前記源信号ができる最小の信号コードである基底で表現できるように，反復計算で前記源信号と前記重みを求める請求項１記載の呼吸曲線抽出法。
前記振動波形の複数の源信号に重みをつけて、足し合わせたものが、前記振動波形と一致するように，かつ，各前記源信号ができる最小の信号コードである基底で表現できるように，反復計算で前記源信号と前記重みを求める請求項１記載の呼吸曲線抽出法。
【請求項3】
前記一般化モルフォロジ成分分析の反復計算において、初期値として、呼吸周波数帯域に相当する係数を１，それ以外の係数を０とする請求項２記載の呼吸曲線抽出法。
【請求項4】
前記変換工程では、離散コサイン変換を使用する請求項１～３のいずれか１項に記載の呼吸曲線抽出法。
【請求項5】
呼吸曲線を複数のエポック分割する分割工程と、
前記複数のエポックを複数のクラスターに分類する分類工程と、
異常呼吸区間のクラスターを抽出する異常呼吸区間を抽出工程と、
を含む異常呼吸事象検出方法。
前記複数のエポックを複数のクラスターに分類する分類工程と、
異常呼吸区間のクラスターを抽出する異常呼吸区間を抽出工程と、
を含む異常呼吸事象検出方法。
【請求項6】
前記分割工程では、特異スペクトル変換を用いる請求項５記載の異常呼吸事象検出方法。
【請求項7】
前記分類工程では、階層的クラスタリング方法を用いる請求項５又は６記載の異常呼吸事象検出方法。
【請求項8】
前記呼吸曲線として、請求項１～３のいずれか１項で抽出した呼吸曲線を用いる請求項５～７のいずれか１項に記載の異常呼吸事象検出方法。
【請求項9】
前記異常呼吸区間は、無呼吸事象区間である請求項５～８のいずれか１項に記載の異常呼吸事象検出方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、呼吸曲線抽出方法と異常呼吸事象抽出方法に関するものである。呼吸曲線を抽出し、特に無呼吸現象や、呼吸系の異常現象を抽出する。
【背景技術】
【0002】
従来、呼吸曲線を抽出する方法があった。特許文献１では、呼吸周波数の標準偏差の逆数により、求めていた。
【先行技術文献】
【特許文献】
【0003】
【特許文献1】特許第５６７００７１号公報
【発明の概要】
【発明が解決しようとする課題】
【0004】
しかし、特許文献１の方法では、正確に呼吸曲線を抽出できなかった。
【0005】
よって、本願の課題は、正確に呼吸曲線を抽出できる呼吸曲線抽出方法と異常呼吸事象抽出方法を提供することである。
【課題を解決するための手段】
【0006】
上記課題を解決するために、２次元の人体の振動波形を測定し、上記測定した第１データを周波数領域データへ変換する変換工程と、上記変換工程で変換された第２データを、モルフォロジ成分分析で、信号表現のスパース性を拠り所にした一般化モルフォロジ成分分析（Ｇｅｎｅｒａｌｉｚｅｄ Ｍｏｒｐｈｏｌｏｇｉｃａｌ Ｃｏｍｐｏｎｅｎｔ Ａｎａｌｙｓｉｓ； ＧＭＣＡ）をし、第２データを得るＧＭＣＡ工程と、上記ＧＭＣＡ工程で得られた第３データを、離散信号データへ変換する逆変換工程と、を含む呼吸曲線抽出法を用いる。
【0007】
また、呼吸曲線を複数のエポック分割する分割工程と、上記複数のエポックを複数のクラスターに分類する分類工程と、異常呼吸区間のクラスターを抽出する異常呼吸区間を抽出工程と、を含む異常呼吸事象検出方法を用いる。
【発明の効果】
【0008】
本発明の呼吸曲線抽出方法と異常呼吸事象検出方法によれば、正確に呼吸曲線を抽出でき、また、異常呼吸事象を抽出できる。
【図面の簡単な説明】
【0009】
【図1】実施の形態の呼吸運動抽出法の概略を示す図
【図2】実施の形態の圧力測定を説明する平面図
【図3】実施の形態のＢＳＳ処理で検討した初期値設定時の相関係数の平均値を示す図
【図4】実施の形態のＢＳＳ処理の各初期値条件における相関係数の分布を示す図
【図5】ｎｉｓｈｉｄａ，ＩＣＡ，ＧＭＣＡのそれぞれの手法で得られた呼吸曲線とＲＩＰｓｕｍの相関係数の平均値を示す図
【図6】各手法における相関係数の分布を示す図
【図7】初期条件Ｉｎｉｔｉａｌ１０のＧＭＣＡを用いた呼吸曲線抽出結果において，被検者ごとに相関係数の平均値を示す図
【図8】特異スペクトル変換の概略を示す図
【図9】ＰＳＧ検査で計測されたＲＩＰ信号からＲＩＰｓｕｍを算出し，ＳＳＴを用いて呼吸曲線の分割を行った結果の例である図
【図10】部分呼吸曲線の結合の階層構造を表すデンドログラムである図
【図11】（ａ）～（ｄ）呼吸曲線を分割しクラスタリングを行った結果を示す図
【発明を実施するための形態】
【0010】
（実施の形態１）一般化モルフォロジ成分分析を用いた呼吸曲線抽出
＜呼吸曲線抽出法の概要＞
最初に、人の圧力信号（観測信号）を解析し、呼吸曲線を得ることを説明する。
＜呼吸運動抽出法＞
図１に，実施の形態の呼吸運動抽出法の概略を示す。
＜呼吸曲線抽出法の概要＞
最初に、人の圧力信号（観測信号）を解析し、呼吸曲線を得ることを説明する。
＜呼吸運動抽出法＞
図１に，実施の形態の呼吸運動抽出法の概略を示す。
【0011】
１-１：圧力信号を測定
測定点のエポック平均圧力値を差し引くことで，各測定点の圧力時系列からバイアスを取り除く。バイアスが取り除かれた圧力時系列で観測信号行列Ｘを作成する。計測された体圧分布時系列は，１睡眠エポック（３０秒間）ごとに以下のＧＭＣＡ（以下で説明）で処理される。
測定点のエポック平均圧力値を差し引くことで，各測定点の圧力時系列からバイアスを取り除く。バイアスが取り除かれた圧力時系列で観測信号行列Ｘを作成する。計測された体圧分布時系列は，１睡眠エポック（３０秒間）ごとに以下のＧＭＣＡ（以下で説明）で処理される。
【0012】
１－２：データ変換
ＸにＤＣＴをする。ＸにＤＣＴを施して得られた係数行列をΘＸとなる。
１－３：ＢＳＳ処理（ＧＭＣＡ）
ΘＸに対してＧＭＣＡを適用すると，混合行列と原信号行列の推定値~Ａと~ΘＳが求まる。
ＸにＤＣＴをする。ＸにＤＣＴを施して得られた係数行列をΘＸとなる。
１－３：ＢＳＳ処理（ＧＭＣＡ）
ΘＸに対してＧＭＣＡを適用すると，混合行列と原信号行列の推定値~Ａと~ΘＳが求まる。
【0013】
１－４：逆データ変換
その後，~ΘＳに、逆ＤＣＴを施す。源信号行列の推定値˜Ｓが得られる。~Ｓにはｎ個の源信号推定値~ｓｊが含まれる。~Ａの列ベクトル~ａｊは，対応する源信号~ｓｊが各観測点の信号にどれだけ含まれているかを表している。以下では，~ａｊの要素~ａｊｉを重み係数とよぶ。源信号の中で最も呼吸運動に近いと判断された源信号が，ＧＭＣＡにより推定された呼吸曲線である。この呼吸曲線に対する重み係数ベクトルは，面状圧力センサ上の呼吸運動の空間分布を表す。
その後，~ΘＳに、逆ＤＣＴを施す。源信号行列の推定値˜Ｓが得られる。~Ｓにはｎ個の源信号推定値~ｓｊが含まれる。~Ａの列ベクトル~ａｊは，対応する源信号~ｓｊが各観測点の信号にどれだけ含まれているかを表している。以下では，~ａｊの要素~ａｊｉを重み係数とよぶ。源信号の中で最も呼吸運動に近いと判断された源信号が，ＧＭＣＡにより推定された呼吸曲線である。この呼吸曲線に対する重み係数ベクトルは，面状圧力センサ上の呼吸運動の空間分布を表す。
【0014】
【0015】
＜１：圧力信号を測定、体圧分布時系列が含む情報＞
面状圧力センサは多数の測定点で圧力を測定し，体重による２次元の圧力分布を計測する装置である。この装置で計測される体圧分布時系列には，呼吸努力や心拍によって引き起こされる周期的な圧力変化が含まれている。
面状圧力センサは多数の測定点で圧力を測定し，体重による２次元の圧力分布を計測する装置である。この装置で計測される体圧分布時系列には，呼吸努力や心拍によって引き起こされる周期的な圧力変化が含まれている。
【0016】
ある１点の圧力センサが測定する信号σは、以下のようにモデル化できる。
σ（ｔ）＝ＳＨ（ｔ）＋ＳＲ（ｔ）＋Ｍ（ｔ）＋Ｂ（ｔ）＋Ｎ（ｔ）
ここで，ＳＨとＳＲは、それぞれ呼吸と心拍による信号，Ｍは、寝返りなどの体動による信号，Ｂは体重がかかることによるバイアス，Ｎはノイズを表す。
面状圧力センサで、計測された体圧分布時系列σ（ｔ）から呼吸運動や心拍運動の情報を適切に分離抽出することができれば，被測定者の睡眠状態の推定に役立つ。
σ（ｔ）＝ＳＨ（ｔ）＋ＳＲ（ｔ）＋Ｍ（ｔ）＋Ｂ（ｔ）＋Ｎ（ｔ）
ここで，ＳＨとＳＲは、それぞれ呼吸と心拍による信号，Ｍは、寝返りなどの体動による信号，Ｂは体重がかかることによるバイアス，Ｎはノイズを表す。
面状圧力センサで、計測された体圧分布時系列σ（ｔ）から呼吸運動や心拍運動の情報を適切に分離抽出することができれば，被測定者の睡眠状態の推定に役立つ。
【0017】
＜１－２データ変換＞
離散コサイン変換（Ｄｉｓｃｒｅｔｅ Ｃｏｓｉｎｅ Ｔｒａｎｓｆｏｒｍ； ＤＣＴ）は離散信号を周波数領域へ変換する手法の一つである。Ｎサンプルの離散信号Ｘ０，．．．，ＸＮ－１に対して，Ｎ個の係数Ｃ０，．．．，ＣＮ－１に変換されるとき，ＤＣＴは次のように定義される。
・・・数１
・・・数２
離散コサイン変換（Ｄｉｓｃｒｅｔｅ Ｃｏｓｉｎｅ Ｔｒａｎｓｆｏｒｍ； ＤＣＴ）は離散信号を周波数領域へ変換する手法の一つである。Ｎサンプルの離散信号Ｘ０，．．．，ＸＮ－１に対して，Ｎ個の係数Ｃ０，．．．，ＣＮ－１に変換されるとき，ＤＣＴは次のように定義される。
【数1】
【数2】
【0018】
離散コサイン変換は、離散信号を周波数領域へ変換する方法の一つであり、他の変換方法でもよい。
後ででてくる逆離散コサイン変換も同様に、周波数領域を離散信号へ変換する方法の一つであり、他の変換方法でもよい。
以下で実施の形態を説明するが、発明の１つの例示であり、これに発明は限定されない。
＜１－３：ＢＳＳ処理（ＧＭＣＡ）＞：一般化モルフォロジ成分分析
ｎ個の信号源から出力した信号をｍ個のセンサで計測して観測信号を得たとき，観測信号から未知の元の信号，すなわち源信号を分離，推定する問題をブラインド信号源分離（ｂｌｉｎｄ ｓｏｕｒｃｅ ｓｅｐａｒａｔｉｏｎ； ＢＳＳ）という。
後ででてくる逆離散コサイン変換も同様に、周波数領域を離散信号へ変換する方法の一つであり、他の変換方法でもよい。
以下で実施の形態を説明するが、発明の１つの例示であり、これに発明は限定されない。
＜１－３：ＢＳＳ処理（ＧＭＣＡ）＞：一般化モルフォロジ成分分析
ｎ個の信号源から出力した信号をｍ個のセンサで計測して観測信号を得たとき，観測信号から未知の元の信号，すなわち源信号を分離，推定する問題をブラインド信号源分離（ｂｌｉｎｄ ｓｏｕｒｃｅ ｓｅｐａｒａｔｉｏｎ； ＢＳＳ）という。
【0019】
ｍ個の観測点でそれぞれ測定された長さＮの観測信号をｘｉ（ｉ＝１，．．．，ｍ）とする。
観測信号を並べた行列をＸ：＝［ｘＴ１，．．．，ｘＴｍ］Ｔとする。
一方，ｎ個ある源信号をｓｊ（ｊ＝１，．．．，ｎ）とする。
源信号を並べた行列をＳ：＝［ｓＴ１，．．．，ｓＴｎ］Ｔとする。ｍ×ｎの混合行列をＡとすると，観測信号行列Ｘと源信号行列Ｓの関係は、Ｘ＝ＡＳと表される。Ｘは既知，Ａ，Ｓは未知であり一般に解くことはできないため，ＢＳＳでは何らかの仮定に基づいて信号源分離を実現する。ＢＳＳの手法の例としては主成分分析，独立成分分析，非負行列因子分解などがあげられる。
観測信号を並べた行列をＸ：＝［ｘＴ１，．．．，ｘＴｍ］Ｔとする。
一方，ｎ個ある源信号をｓｊ（ｊ＝１，．．．，ｎ）とする。
源信号を並べた行列をＳ：＝［ｓＴ１，．．．，ｓＴｎ］Ｔとする。ｍ×ｎの混合行列をＡとすると，観測信号行列Ｘと源信号行列Ｓの関係は、Ｘ＝ＡＳと表される。Ｘは既知，Ａ，Ｓは未知であり一般に解くことはできないため，ＢＳＳでは何らかの仮定に基づいて信号源分離を実現する。ＢＳＳの手法の例としては主成分分析，独立成分分析，非負行列因子分解などがあげられる。
【0020】
冗長な信号表現基底の辞書を与えた場合に信号や画像がそれらの一部の基底を用いて表現できることを仮定し，信号分解とノイズ除去を行う手法としてモルフォロジ成分分析（Ｍｏｒｐｈｏｌｏｇｉｃａｌ Ｃｏｍｐｏｎｅｎｔ Ａｎａｌｙｓｉｓ；ＭＣＡ）（Ｊ．－Ｌ Ｓｔａｒｃｋ，Ｙ．Ａ，Ｊ．Ｂｏｂｉｎ，Ｍｉｃｈａｅｌ Ｅｌａｄ，ａｎｄ Ｄ． Ｄｏｎｏｈｏ．Ｍｏｒｐｈｏｌｏｇｉｃａｌ ｃｏｍｐｏｎｅｎｔ ａｎａｌｙｓｉｓ．Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ ＳＰＩＥ －Ｔｈｅ Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｓｏｃｉｅｔｙ ｆｏｒ Ｏｐｔｉｃａｌ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｖｏｌ．５９１４，，０８ ２００５．）が提案された。ＭＣＡは，ある信号ｙが辞書Ｄ中の特定の基底Φｋでスパースに表現（ほとんどの基底に対する係数がゼロになる）され，Ｄ個の形態成分の線形結合で表せると仮定する。ただし，辞書Ｄは正規直交基底ΦｋをＤ＝［ΦＴ１，．．．，ΦＴＤ］Ｔとなるように結合させたものである。このとき，ｙは、数３となる。
・・・数３
【数3】
【0021】
ただし，αｋはｙをΦｋに対応づける係数ベクトルである。
係数ベクトルαはスパースであることから次の数４の最小化問題が得られる。
・・・数４
∥．．．∥０は、ｌ０ノルムであり，０でない成分の個数を表す。この最小化問題を解く方法としてＧＭＣＡが提案された（一般化ＭＣＡ）。ここで信号ｙに固有のスパース分解が存在し，ｙ＝α*Ｄがその解であるとき，α*＝ΔＤ（ｙ）とする。
係数ベクトルαはスパースであることから次の数４の最小化問題が得られる。
【数4】
∥．．．∥０は、ｌ０ノルムであり，０でない成分の個数を表す。この最小化問題を解く方法としてＧＭＣＡが提案された（一般化ＭＣＡ）。ここで信号ｙに固有のスパース分解が存在し，ｙ＝α*Ｄがその解であるとき，α*＝ΔＤ（ｙ）とする。
【0022】
ＧＭＣＡは，源信号が辞書Ｄ中の特定の基底Φｋでスパースに表現され，Ｄ個の形態成分の線形結合で表せると仮定する。ＧＭＣＡは信号表現のスパース性を拠り所にしてＢＳＳを試みる。
各源信号ｓｊが辞書Ｄ中の特定の基底ΦｋによりＤ個の形態成分の線形結合でスパースに表現されると仮定すると，源信号ｓｊは、数５となる。
・・・数５
各源信号ｓｊが辞書Ｄ中の特定の基底ΦｋによりＤ個の形態成分の線形結合でスパースに表現されると仮定すると，源信号ｓｊは、数５となる。
【数5】
【0023】
ただし，αｊｋはｓｊをΦｋに対応づける係数ベクトルである。
源信号ｓｊが辞書Ｄを用いてスパースに表されることから，観測信号行列Ｘの分解に関して次の数６の最適化問題が定義できる。
・・・数６
源信号ｓｊが辞書Ｄを用いてスパースに表されることから，観測信号行列Ｘの分解に関して次の数６の最適化問題が定義できる。
【数6】
【0024】
ｌ０ノルム（∥・∥０）をｌ１ノルムに緩和しラグランジュ未定乗数λｊを導入することで，数６の問題は以下の数７に近似できる。
・・・数７
【数7】
【0025】
ただし，∥・∥２
２はフロベニウスノルムであり，∥Ｐ∥２
２＝Ｔｒａｃｅ（ＰＴＰ）である。
ＧＭＣＡでは，反復縮小しきい値アルゴリズムで数７の最適化問題を解き，観測信号行列Ｘから源信号行列と混合行列の推定値~Ｓと~Ａを求める。すなわち，λｊを大きな値から徐々に小さくすることで，スパース性についての制約を緩和しながら最適化計算を繰り返し，式７の解を求める。
ＧＭＣＡでは，反復縮小しきい値アルゴリズムで数７の最適化問題を解き，観測信号行列Ｘから源信号行列と混合行列の推定値~Ｓと~Ａを求める。すなわち，λｊを大きな値から徐々に小さくすることで，スパース性についての制約を緩和しながら最適化計算を繰り返し，式７の解を求める。
【0026】
式７を解くアルゴリズムでは源信号行列Ｓを反復計算により推定するとき，各計算時にＤによる変換を行う必要があり計算コストが多くかかる。ＧＭＣＡの高速アルゴリズムは観測信号ｘｉに固有のスパース分解が存在すると仮定し，係数空間で信号源分離を行う。観測信号行列ｘｉのスパース分解の係数行列をΘＸ：＝［ΔＤ（ｘ１）Ｔ，．．．，ΔＤ（ｘｍ）Ｔ］Ｔ，源信号行列ｓｉのスパース分解の係数行列をΘＳ：＝［ΔＤ（ｓ１）Ｔ，．．．，ΔＤ（ｓｎ）Ｔ］Ｔと定義する。高速アルゴリズムは数７を次の最適化問題で近似する。
・・・数８
【数8】
【0027】
辞書Ｄが正規直交基底でない限り数７と数８は等価ではない。辞書Ｄが冗長である場合、厳密な数学的証明は与えられないが，実験により良い結果が得られる。
数８でＡを固定するとΘＳは次数により求められる。
・・・数９
数８でＡを固定するとΘＳは次数により求められる。
【数9】
【0028】
ただし，Ａ†は疑似逆行列である。また，Ｓλはソフトしきい値作用素（ｓｏｆｔ－ｔｈｒｅｓｈｏｌｄｉｎｇ）で，ｙ［ｔ］を入力としたとき，次のように作用する。
・・・数１０
数８でΘＳを固定するとＡは最小二乗推定となり，次数１１により求められる。
【数10】
数８でΘＳを固定するとＡは最小二乗推定となり，次数１１により求められる。
【0029】
【数11】
数８を解く反復縮小しきい値アルゴリズムは次のようになる。
１．観測信号ｘｉに離散コサイン変換を適用し，ΘＸを計算する。
２．初期しきい値λ０，しきい値の最小値λｍｉｎ，混合行列の初期値Ａｉｎｉを設定する。
３．１ 数９によりΘＳを推定する。
３．２ 数１１によりＡを推定する。
３．３ しきい値λを小さくする。
３．４ λ＞λｍｉｎの場合，３．１に戻る。
【0030】
４．しきい値λ＝λｍｉｎとなった場合，計算を終了する。
しきい値の最小値λｍｉｎをどれほどまで小さくするかは，最終的に求まる源信号のスパース性を左右する。その値が大きい場合には源信号の推定値はより少数の基底で表され，小さい場合にはより多数の基底で表される。しきい値λが大きいとき，源信号は最も重要な係数から推定され，しきい値λが小さくなるにつれて，混合行列と源信号行列の積と観測信号行列の誤差が小さくなる。この反復縮小しきい値アルゴリズムは，データの最も重要な特徴を最初に推定し，徐々に細部を推定することでロバスト性を確保している。
しきい値の最小値λｍｉｎをどれほどまで小さくするかは，最終的に求まる源信号のスパース性を左右する。その値が大きい場合には源信号の推定値はより少数の基底で表され，小さい場合にはより多数の基底で表される。しきい値λが大きいとき，源信号は最も重要な係数から推定され，しきい値λが小さくなるにつれて，混合行列と源信号行列の積と観測信号行列の誤差が小さくなる。この反復縮小しきい値アルゴリズムは，データの最も重要な特徴を最初に推定し，徐々に細部を推定することでロバスト性を確保している。
【0031】
１－３－Ａ 一般化モルフォロジ成分分析の初期値設定
数１のように体圧分布時系列には主として呼吸周波数帯域の圧力変化，心拍周波数帯域の圧力変化，寝返りなどの体動による圧力変化が含まれている。この事前知識をＧＭＣＡに与えることで呼吸運動抽出精度が向上する。
与える事前知識として混合行列Ａｉｎｉに呼吸周波数帯域の情報がえられる。
ＧＭＣＡに与える初期値として以下の４条件を検討する。
数１のように体圧分布時系列には主として呼吸周波数帯域の圧力変化，心拍周波数帯域の圧力変化，寝返りなどの体動による圧力変化が含まれている。この事前知識をＧＭＣＡに与えることで呼吸運動抽出精度が向上する。
与える事前知識として混合行列Ａｉｎｉに呼吸周波数帯域の情報がえられる。
ＧＭＣＡに与える初期値として以下の４条件を検討する。
【0032】
１．Ｒａｎｄｏｍ：Ａｉｎｉに標準正規分布に従うランダムな値を与える。
２．ＰＣＡｂａｓｅｄ：観測信号行列の係数行列の分散共分散行列Σ＝ΘＸΘＴ Ｘの固有ベクトルｖｊ（ｊ＝１，２，．．，ｎ）を主成分分析（ＰＣＡ）により求める。Ａｉｎｉはｎ個の固有ベクトルｖｊ（ｊ＝１，２，．．，ｍ）により与える。
２．ＰＣＡｂａｓｅｄ：観測信号行列の係数行列の分散共分散行列Σ＝ΘＸΘＴ Ｘの固有ベクトルｖｊ（ｊ＝１，２，．．，ｎ）を主成分分析（ＰＣＡ）により求める。Ａｉｎｉはｎ個の固有ベクトルｖｊ（ｊ＝１，２，．．，ｍ）により与える。
【0033】
３．Ｉｎｉｔｉａｌ１：一つ目の源信号の係数ベクトルΘｓ１について，呼吸周波数帯域に相当する係数を１，それ以外の係数を０とする。重み係数ベクトルについては，一つ目の源信号の係数ベクトルΘｓ１に対応する重み係数ベクトルａ１をａ１＝ΘＸΘＴ
ｓ１として与える。残りの重み係数ベクトルａｊ（ｊ＝１，２，．．．，ｎ）には標準正規分布に従うランダムな値を与える。
【0034】
４．Ｉｎｉｔｉａｌ１０：Ｉｎｉｔｉａｌ１と同様に，一つ目の源信号の係数ベクトルΘｓ１について呼吸周波数帯域に相当する係数を１，それ以外の係数を０とする。
さらに，二つ目の源信号の係数ベクトルΘｓ２については心拍周波数帯域に相当する係数を１，それ以外の係数を０とする。
さらに，二つ目の源信号の係数ベクトルΘｓ２については心拍周波数帯域に相当する係数を１，それ以外の係数を０とする。
【0035】
三つ目の源信号の係数ベクトルΘｓ３については呼吸周波数帯域よりも低周波の周波数帯域に相当する係数を１，それ以外の係数を０とする。これは，寝返りなどの体動成分に相当する。
四つ目の源信号の係数ベクトルΘｓ４については呼吸周波数帯域と心拍周波数帯域の間の周波数帯域に相当する係数を１，それ以外の係数を０とする。
四つ目の源信号の係数ベクトルΘｓ４については呼吸周波数帯域と心拍周波数帯域の間の周波数帯域に相当する係数を１，それ以外の係数を０とする。
【0036】
残りの源信号の係数ベクトルは心拍周波数帯域よりも高周波の周波数帯域を分割し，分割した周波数帯域に相当する係数を１，それ以外の係数を０とする。
以上の設定により定義された源信号行列の係数行列ΘＳの疑似逆行列Θ†ｓと観測信号行列の係数ベクトルΘＸを用いて混合行列の初期値をＡｉｎｉ＝ΘＸΘ†Ｓと与える。
以上の設定により定義された源信号行列の係数行列ΘＳの疑似逆行列Θ†ｓと観測信号行列の係数ベクトルΘＸを用いて混合行列の初期値をＡｉｎｉ＝ΘＸΘ†Ｓと与える。
【0037】
＜１－４：逆データ変換＞
ＤＣＴの逆変換を逆離散コサイン変換（Ｉｎｖｅｒｓｅ ＤＣＴ； ＩＤＣＴ）といい，次のように定義される。
・・・数１２
ＤＣＴ により得られる係数は離散信号をさまざまな周波数のコサイン関数の和で表現したときのそれぞれの周波数のコサイン関数の重みを表す。
ＤＣＴの逆変換を逆離散コサイン変換（Ｉｎｖｅｒｓｅ ＤＣＴ； ＩＤＣＴ）といい，次のように定義される。
【数12】
ＤＣＴ により得られる係数は離散信号をさまざまな周波数のコサイン関数の和で表現したときのそれぞれの周波数のコサイン関数の重みを表す。
【0038】
＜実施例＞
一般化モルフォロジ成分分析を用いた呼吸曲線の抽出
１：データセット
体圧分布時系列の計測には，住友理工株式会社製ＳｍａｒｔＲｕｂｂｅｒ（登録商標）のセンサ１１（以降，ＳＲのセンサ１１）を用いた。使用されたＳＲのセンサ１１は，寸法が１０５ｃｍ×９０ｃｍのマット状の圧力センサで，３２×２５＝８００点の計測点をもつ。各測定点の計測周期は２０Ｈｚである。圧力値はおよそ１００～２０００の範囲の整数で出力される。解析対象としたデータセットは，太田睡眠科学センターでＰＳＧ検査と並行して計測された。
一般化モルフォロジ成分分析を用いた呼吸曲線の抽出
１：データセット
体圧分布時系列の計測には，住友理工株式会社製ＳｍａｒｔＲｕｂｂｅｒ（登録商標）のセンサ１１（以降，ＳＲのセンサ１１）を用いた。使用されたＳＲのセンサ１１は，寸法が１０５ｃｍ×９０ｃｍのマット状の圧力センサで，３２×２５＝８００点の計測点をもつ。各測定点の計測周期は２０Ｈｚである。圧力値はおよそ１００～２０００の範囲の整数で出力される。解析対象としたデータセットは，太田睡眠科学センターでＰＳＧ検査と並行して計測された。
【0039】
ＰＳＧ検査とは、睡眠時無呼吸症候群（ＳＡＳ）を確定するための、ポリソムノグラフィー（Ｐｏｌｙｓｏｍｎｏｇｒａｐｈｙ； ＰＳＧ）という検査である。これは医療機関に泊りがけで行う検査法で、頭や顔、体の必要な部位にテープで電極を貼りつけ、実際に一晩眠りながら脳波や呼吸、眼球、筋肉の動きなどを記録し、睡眠の状態について調べるものである。この結果を解析し、正常な睡眠の経過と比較することで、睡眠時無呼吸症候群（ＳＡＳ）の診断を行う。
【0040】
計測に際し，ＳＲのセンサ１１は検査室で被検者（人１０）が寝るベッドのマットレス上に敷かれた。ＰＳＧ検査では被検者の胸部（ｃｈｅｓｔ）と腹部（ａｂｄｏｍｅｎ）に呼吸インダクタンス・プレチスモグラフ（ＲＩＰ）ベルトが取り付けられ呼吸運動が計測された。胸部と腹部のＲＩＰベルトの出力電位の和であるＲＩＰｓｕｍ信号ｒｓｕｍを，体圧分布時系列の処理結果と照合する参考信号とした。呼吸運動抽出結果の評価指標には，ＲＩＰｓｕｍとの相関係数を用いた。提案手法が出力する源信号推定値~ｓｊは複数あるため，±０．５秒の時間シフトを許して相関係数を計算し，相関係数の絶対値が最大となる源信号を呼吸曲線，その時の相関係数の絶対値を評価指標とした。
【0041】
２：呼吸曲線の抽出精度
２３名の被検者（人１０）の一晩の睡眠時体圧分布時系列に対してＧＭＣＡを用いた呼吸曲線抽出法を適用し呼吸曲線の抽出精度を評価した。源信号数の設定はｎ＝１０である。合計２２１５８個のエポックが解析対象となり，各エポックにおいて睡眠時体圧分布時系列から呼吸曲線を抽出し相関係数を算出した。
２３名の被検者（人１０）の一晩の睡眠時体圧分布時系列に対してＧＭＣＡを用いた呼吸曲線抽出法を適用し呼吸曲線の抽出精度を評価した。源信号数の設定はｎ＝１０である。合計２２１５８個のエポックが解析対象となり，各エポックにおいて睡眠時体圧分布時系列から呼吸曲線を抽出し相関係数を算出した。
【0042】
図３は，１－３―Ａで検討した初期値設定時の相関係数の平均値を示す。相関係数の平均値はＩｎｉｔｉａｌ１０，Ｉｎｉｔｉａｌ１，Ｒａｎｄｏｍ，ＰＣＡｂａｓｅｄの順番に高い結果となった。Ｉｎｉｔｉａｌ１０，Ｉｎｉｔｉａｌ１の相関係数の平均値が高くなったことから，体圧分布時系列は呼吸周波数帯域の信号を含むという事前知識に基づいた初期値を用いることで呼吸運動の抽出精度が向上することがわかる。
【0043】
図４は，各初期値条件における相関係数の分布を示す。Ｉｎｉｔｉａｌ１０，Ｉｎｉｔｉａｌ１，Ｒａｎｄｏｍの初期値条件において相関係数０．１と０．８５付近に２つのピークが存在する。Ｉｎｉｔｉａｌ１０の分布は他の条件の分布と比較して相関係数０．１のピーク近傍の分布が小さく，相関係数０．８５のピークの近傍の分布が大きい。３つの初期値条件に関して２つのピークの位置に大きな違いはない。すなわち，Ｉｎｉｔｉａｌ１０の他の初期値条件に対する呼吸運動の抽出精度の向上はＲＩＰｓｕｍ により近い呼吸曲線を抽出しているのではなく，他の初期値条件では呼吸曲線を抽出できていないエポックにおいて呼吸曲線の抽出を可能にしていることを示している。
【0044】
次に，同じ２３名の睡眠時体圧分布時系列に対して，西田らの手法（西田佳史，武田正資，森武俊，溝博，佐藤知正．圧力センサによる睡眠中の呼吸・体位の無侵襲・無拘束な計測．日本ロボット学会誌，Ｖｏｌ．１６，Ｎｏ．５，ｐｐ．７０５－７１１，ｊｕｌ１９９８．：複数のセンサからの信号の位相差から呼吸曲線を計算する方法），独立成分分析（ＩＣＡ），ＧＭＣＡを用いた呼吸曲線抽出法を適用し呼吸曲線抽出精度を評価する。ＧＭＣＡを用いて呼吸曲線を抽出する際の初期値にはＩｎｉｔｉａｌ１０を用いた。ＩＣＡでは，データ行列に対して基底数５でＩＣＡを実行することで得られる源信号を抽出結果とした。
【0045】
図５は，西田らの方法，ＩＣＡ，ＧＭＣＡのそれぞれの手法で得られた呼吸曲線とＲＩＰｓｕｍの相関係数の平均値を示す。ＧＭＣＡを用いた呼吸曲線抽出法は他の手法と比較して精度良く呼吸曲線を抽出することができることが，このグラフから分かる。
図６は，各手法における相関係数の分布を示す。ＧＭＣＡはＩＣＡと比較して，ピークの位置の相関係数が高く，よりＲＩＰｓｕｍに近い呼吸曲線を抽出することができていることが分かる。
図６は，各手法における相関係数の分布を示す。ＧＭＣＡはＩＣＡと比較して，ピークの位置の相関係数が高く，よりＲＩＰｓｕｍに近い呼吸曲線を抽出することができていることが分かる。
【0046】
図７は，初期条件Ｉｎｉｔｉａｌ１０のＧＭＣＡを用いた呼吸曲線抽出結果において，被検者ごとに相関係数の平均値を示す。１２番の被検者については呼吸曲線の抽出精度が低い結果となった。この被検者は２３名の被検者の中でも検査中の睡眠時間が短く，何度も離床を繰り返していた。このような計測不可能な時間が含まれているため抽出精度が低くなっている可能性がある。しかし，安定して体圧分布時系列を計測できている時間においても抽出精度はあまり高くない。このような被検者によっては抽出精度が悪くなってしまう原因の調査と対策は今後の課題である。
【0047】
（実施の形態２）呼吸曲線から無呼吸事象検出
上記で、呼吸曲線が抽出された。ここでは、その呼吸曲線から無呼吸事象を検出する。なお、呼吸曲線は別の方法で得たものを用いてもよい。
＜２－１：分割＞
無呼吸事象を検出するためにはさまざまな計測方法で計測される呼吸曲線を正常呼吸区間と無呼吸区間に判別する必要がある。呼吸曲線には被検者の体格や呼吸方式，呼吸数の違いによる被検者間の差が存在する。同一被検者に対する検査において被検者の状態の変化によって呼吸曲線は変化する。検査機器と被検者の関係は睡眠中の被検者の動きにより容易に変化し，その変化により呼吸曲線は変化する。呼吸曲線は正常呼吸と無呼吸の変化だけでなく，このような他の要因による変化を多分に含む。このため，すべての呼吸曲線の変化をモデル化することで適切に呼吸曲線を分割することは困難である。そこで実施の形態では，ノンパラメトリックな時系列解析手法であるＳＳＴを用いて呼吸曲線を分割する。
上記で、呼吸曲線が抽出された。ここでは、その呼吸曲線から無呼吸事象を検出する。なお、呼吸曲線は別の方法で得たものを用いてもよい。
＜２－１：分割＞
無呼吸事象を検出するためにはさまざまな計測方法で計測される呼吸曲線を正常呼吸区間と無呼吸区間に判別する必要がある。呼吸曲線には被検者の体格や呼吸方式，呼吸数の違いによる被検者間の差が存在する。同一被検者に対する検査において被検者の状態の変化によって呼吸曲線は変化する。検査機器と被検者の関係は睡眠中の被検者の動きにより容易に変化し，その変化により呼吸曲線は変化する。呼吸曲線は正常呼吸と無呼吸の変化だけでなく，このような他の要因による変化を多分に含む。このため，すべての呼吸曲線の変化をモデル化することで適切に呼吸曲線を分割することは困難である。そこで実施の形態では，ノンパラメトリックな時系列解析手法であるＳＳＴを用いて呼吸曲線を分割する。
【0048】
（ＳＳＴ）
特異スペクトル変換（Ｓｉｎｇｕｌａｒ Ｓｐｅｃｔｒｕｍ Ｔｒａｎｓｆｏｒｍａｔｉｏｎ； ＳＳＴ）（井出剛，杉山将．異常検知と変化検知．講談社，２０１５．）は特異値分解により時系列データの特徴パターンを抽出し，それに基づいて変化度を求める手法である。
図８に概略を示す。
時系列データｙ＝｛ｙ１，ｙ２，．．．，ｙｔ，．．．｝を考える．時刻ｔの直前のｗ個のデータをまとめた部分時系列ｘ（ｔ）をｘ（ｔ）＝［ｙ（ｔ-ｗ），ｙ（ｔ-ｗ＋１），．．．，ｙ（ｔ-１）］Ｔとする。時刻ｔの過去側と未来側において，部分時系列を１時刻ずつずらしてｎ個並べ，ｗ行ｎ列の２つのデータ行列Ｘ（ｔ），Ｚ（ｔ）をＸ（ｔ）＝［ｘ（ｔ-ｎ＋１），．．．，ｘ（ｔ-１），ｘ（ｔ）］
Ｚ（ｔ）＝［ｘ（ｔ-ｎ＋１＋Ｌ），．．．，ｘ（ｔ-１＋Ｌ），ｘ（ｔ＋Ｌ）］のように構成する。
特異スペクトル変換（Ｓｉｎｇｕｌａｒ Ｓｐｅｃｔｒｕｍ Ｔｒａｎｓｆｏｒｍａｔｉｏｎ； ＳＳＴ）（井出剛，杉山将．異常検知と変化検知．講談社，２０１５．）は特異値分解により時系列データの特徴パターンを抽出し，それに基づいて変化度を求める手法である。
図８に概略を示す。
時系列データｙ＝｛ｙ１，ｙ２，．．．，ｙｔ，．．．｝を考える．時刻ｔの直前のｗ個のデータをまとめた部分時系列ｘ（ｔ）をｘ（ｔ）＝［ｙ（ｔ-ｗ），ｙ（ｔ-ｗ＋１），．．．，ｙ（ｔ-１）］Ｔとする。時刻ｔの過去側と未来側において，部分時系列を１時刻ずつずらしてｎ個並べ，ｗ行ｎ列の２つのデータ行列Ｘ（ｔ），Ｚ（ｔ）をＸ（ｔ）＝［ｘ（ｔ-ｎ＋１），．．．，ｘ（ｔ-１），ｘ（ｔ）］
Ｚ（ｔ）＝［ｘ（ｔ-ｎ＋１＋Ｌ），．．．，ｘ（ｔ-１＋Ｌ），ｘ（ｔ＋Ｌ）］のように構成する。
【0049】
Ｘ（ｔ）とＺ（ｔ）をそれぞれ履歴行列，テスト行列とよぶ。Ｌは、履歴行列とテスト行列の相互位置を定める非負整数である。履歴行列は時刻ｔ-ｗ-ｎ＋１からｔのデータから構成され，時刻ｔに対して過去の情報が含まれる。テスト行列は時刻ｔ-ｗ-ｎ＋１＋Ｌからｔ＋Ｌのデータから構成され，時刻ｔに対して未来の情報が含まれる。
【0050】
次に、履歴行列Ｘ（ｔ）とテスト行列Ｚ（ｔ）に特異値分解を行い，特異値の最も大きい左特異ベクトルをそれぞれα（ｔ），β（ｔ）とする。α（ｔ），β（ｔ）はそれぞれ過去と未来の特徴的なパターンを表現している。時系列データｙの振る舞いが時刻ｔの前後で変化しない場合，α（ｔ）とβ（ｔ）は一致する。逆に，ｙの振る舞いが時刻ｔの前後で変化したとき，α（ｔ）とβ（ｔ）は一致しない。α（ｔ）とβ（ｔ）の類似度を評価することでｙの振る舞いの変化を数値化することができる。変化度ａ（ｔ）は，ａ（ｔ）＝１-α（ｔ）Ｔβ（ｔ）と定義される。変化度ａ（ｔ）はα（ｔ）とβ（ｔ）が一致している時に最小となり，α（ｔ）とβ（ｔ）が直交している時に最大となる。
【0051】
呼吸曲線に対してＳＳＴを適用し、変化度ａ（ｔ）を算出し，変化度ａ（ｔ）が極大値をとる点を変化点とする。変化点間の呼吸曲線を部分呼吸曲線と定義する。部分呼吸曲線には少なくとも一呼吸波形は含まれることが望ましい。睡眠中の一呼吸にかかる時間は３～６秒であるので，実施の形態において時刻ｔの前後６秒間において変化度ａ（ｔ）が最大となるとき，時刻ｔを変化点と定義する。
【0052】
図９は，ＰＳＧ検査で計測されたＲＩＰ信号からＲＩＰｓｕｍを算出し，ＳＳＴを用いて呼吸曲線の分割を行った結果の例である。ＲＩＰ信号のサンプリングレートは１６［Ｈｚ］であり，ＳＳＴのパラメータはｗ＝ｎ＝２０，Ｌ＝１０である。上段はＲＩＰｓｕｍであり，下段はＳＳＴにより算出した変化度ａ（ｔ）である。赤線は変化点を示す。呼吸曲線を分割した結果，ＲＩＰｓｕｍの振幅が低下している時刻８６［ｓ］から１０５［ｓ］付近の区間を部分呼吸曲線として取り出すことができている。実施の形態ではＳＳＴのパラメータは同様の設定ｗ＝ｎ＝２０，Ｌ＝１０を採用する。
【0053】
＜２－２：クラスタリング＞部分呼吸曲線の分類
本節では，部分呼吸曲線の特徴量を定義し，クラスタリングにより部分呼吸曲線を自動的に分類する手法について説明する。
（階層的クラスタリング）
クラスタリングとは，外的基準を用いずにデータを類似するものをまとめたクラスターに分け，自動的に分類を行う方法である。階層的クラスタリングはクラスタリングの手法の一つであり，データがそれぞれ個々のクラスターの状態から類似のクラスターを順に結合し，クラスターの階層構造を生成する手法である。クラスターの階層構造を樹形図で表現したものをデンドログラムという。
階層的クラスタリングの流れは次のようになる。
本節では，部分呼吸曲線の特徴量を定義し，クラスタリングにより部分呼吸曲線を自動的に分類する手法について説明する。
（階層的クラスタリング）
クラスタリングとは，外的基準を用いずにデータを類似するものをまとめたクラスターに分け，自動的に分類を行う方法である。階層的クラスタリングはクラスタリングの手法の一つであり，データがそれぞれ個々のクラスターの状態から類似のクラスターを順に結合し，クラスターの階層構造を生成する手法である。クラスターの階層構造を樹形図で表現したものをデンドログラムという。
階層的クラスタリングの流れは次のようになる。
【0054】
１．個々のデータをクラスターとして，全ての組み合わせのクラスター間の距離を求める。
２．最も距離の近いクラスター対を結合する。
３．クラスター間の距離を再計算する。
４．クラスター数が２以上のとき，２．に戻る。
２．最も距離の近いクラスター対を結合する。
３．クラスター間の距離を再計算する。
４．クラスター数が２以上のとき，２．に戻る。
【0055】
Ｗａｒｄ法はクラスター間の距離を再計算する手法の一つであり，クラスター内の分散が最も小さくなるようにクラスターを順次結合していく方法である。個々のデータ点間の距離はユークリッド距離を用いる。Ｎ個のデータの集合Ｘ＝｛ｘ１，．．．，ｘＮ｝があり，ＸがＣ個のクラスターＧ１，．．．，ＧＣに分割されているとする。クラスターＧｉの重心Ｍ（Ｇｉ）はクラスターＧｉの要素数｜Ｇｉ｜を用いて，Ｍ（Ｇｉ）＝１｜Ｇｉ｜Σｘｎ∈Ｇｉｘｎとなる。クラスターＧｉに関して重心と各個体との距離の差の２乗和Ｅ（ＧｉはＥ（Ｇｉ）＝Σｘｎ∈Ｇｉ∥ｘｎ-Ｍ（Ｇｉ）∥２となる。ここでクラスター間の距離としてΔＥをΔＥ（Ｇｉ，Ｇｊ）＝Ｅ（Ｇｉ∪Ｇｊ）-Ｅ（Ｇｉ）-Ｅ（Ｇｊ）とおき，ΔＥが最小となるクラスター対を結合する。クラスター内分散の総和Ｅは、Ｅ＝ΣＣｉ＝１Ｅ（Ｇｉ）となり，Ｗａｒｄ法はＥが最小となるクラスター構造を生成する。
【0056】
呼吸曲線は、呼吸運動による身体の動きを計測したものである。睡眠中の計測機器と被検者間の関係は被検者の寝返りなどにより容易に変化する。また，呼吸運動による身体の動きは被検者による差や，呼吸状態，睡眠状態により変化する。このように，呼吸曲線はさまざまな要因により変化するため，正確にモデリングすることは困難である。ＰＳＧ検査における無呼吸低呼吸の判定基準も呼吸曲線のモデルを網羅しているわけではない。そこで，データ同士を教師なしにクラスターにまとめる手法であるクラスタリングを用いて部分呼吸曲線を分類する。
【0057】
実施の形態では，部分呼吸曲線から自己回帰係数と部分呼吸曲線振幅比を用いて部分呼吸曲線の特徴を抽出し，部分呼吸曲線の特徴量をデータとして階層的クラスタリングを行い，部分呼吸曲線を分類する。部分呼吸曲線を自動的に分類することで正常呼吸区間や無呼吸区間に分類することを目的とする。また，被検者間や呼吸方式，センシングの状態により分類されることが好ましい。
【0058】
＜実施例＞
無呼吸区間の抽出
１：データセット
太田睡眠科学センターで実施されたＰＳＧ検査で計測された胸部と腹部のＲＩＰベルトの出力値の和ＲＩＰｓｕｍを算出し，ＲＩＰｓｕｍ信号を呼吸曲線として解析を行った。２３名の被検者のデータからランダムに６エポックの区間を５０区間抽出し，計３００エポックを解析の対象とした。
ＰＳＧ検査のデータは、３０秒を１睡眠エポックとして区切られている。解析対象のＰＳＧ検査のデータには，専門の検査技師により，１０秒以上の無呼吸区間が発生したエポックに閉塞性無呼吸などの呼吸異常のラベルがつけられている。呼吸異常のラベルが付けられていないエポックを正常呼吸エポック，呼吸異常のラベルが付けられているエポックを無呼吸低呼吸エポックとよぶことにする。また，１０秒以上の無呼吸区間が発生することを無呼吸イベントとよぶことにする。
無呼吸区間の抽出
１：データセット
太田睡眠科学センターで実施されたＰＳＧ検査で計測された胸部と腹部のＲＩＰベルトの出力値の和ＲＩＰｓｕｍを算出し，ＲＩＰｓｕｍ信号を呼吸曲線として解析を行った。２３名の被検者のデータからランダムに６エポックの区間を５０区間抽出し，計３００エポックを解析の対象とした。
ＰＳＧ検査のデータは、３０秒を１睡眠エポックとして区切られている。解析対象のＰＳＧ検査のデータには，専門の検査技師により，１０秒以上の無呼吸区間が発生したエポックに閉塞性無呼吸などの呼吸異常のラベルがつけられている。呼吸異常のラベルが付けられていないエポックを正常呼吸エポック，呼吸異常のラベルが付けられているエポックを無呼吸低呼吸エポックとよぶことにする。また，１０秒以上の無呼吸区間が発生することを無呼吸イベントとよぶことにする。
【0059】
２：部分呼吸曲線のクラスター構造
まず，５０個の長さ６エポックのＲＩＰｓｕｍ信号を２－１のＳＳＴの手法を用いて部分呼吸曲線に分割した。分割された部分呼吸曲線の数は６２５個となった。次に，２－２の階層的クラスタリングの手法を用いて部分呼吸曲線のクラスタリングを行った。図１０は，部分呼吸曲線の結合の階層構造を表すデンドログラムである。縦軸はクラスターを結合した時のクラスター内の分散の総和を表す。今回の解析では，クラスター数が１３となる段階でクラスターの結合を停止し，部分呼吸曲線のラベルとする。つまり，部分呼吸曲線は１３個のクラスターに分類される。赤色の横線はクラスターの結合を停止した段階を示す。赤線より下で結合された部分呼吸曲線同士は同じクラスターに属し，赤線より上で結合した部分呼吸曲線同士は異なるクラスターに属する。図１０で左に位置するクラスターから順にクラスター１，クラスター２，．．．，クラスター１３とする。
まず，５０個の長さ６エポックのＲＩＰｓｕｍ信号を２－１のＳＳＴの手法を用いて部分呼吸曲線に分割した。分割された部分呼吸曲線の数は６２５個となった。次に，２－２の階層的クラスタリングの手法を用いて部分呼吸曲線のクラスタリングを行った。図１０は，部分呼吸曲線の結合の階層構造を表すデンドログラムである。縦軸はクラスターを結合した時のクラスター内の分散の総和を表す。今回の解析では，クラスター数が１３となる段階でクラスターの結合を停止し，部分呼吸曲線のラベルとする。つまり，部分呼吸曲線は１３個のクラスターに分類される。赤色の横線はクラスターの結合を停止した段階を示す。赤線より下で結合された部分呼吸曲線同士は同じクラスターに属し，赤線より上で結合した部分呼吸曲線同士は異なるクラスターに属する。図１０で左に位置するクラスターから順にクラスター１，クラスター２，．．．，クラスター１３とする。
【0060】
３：無呼吸区間に対応するクラスターの特定
無呼吸区間に対応するクラスターの特定ＰＳＧ検査で付けられたラベルを用いて無呼吸区間に対応するクラスターを特定する。正常呼吸区間に対応するクラスターに属する部分呼吸曲線を正常部分呼吸曲線，無呼吸区間に対応するクラスターに属する部分呼吸曲線を無呼吸部分呼吸曲線とする。正常部分呼吸曲線の次に無呼吸部分呼吸曲線が続き，その無呼吸区間の長さが１０秒以上となる時，正常部分呼吸曲線と無呼吸部分呼吸曲線の間の変化点を含むエポックは、無呼吸低呼吸エポックとなる。無呼吸低呼吸エポックには無呼吸部分呼吸曲線の開始点が含まれる。
無呼吸区間に対応するクラスターの特定ＰＳＧ検査で付けられたラベルを用いて無呼吸区間に対応するクラスターを特定する。正常呼吸区間に対応するクラスターに属する部分呼吸曲線を正常部分呼吸曲線，無呼吸区間に対応するクラスターに属する部分呼吸曲線を無呼吸部分呼吸曲線とする。正常部分呼吸曲線の次に無呼吸部分呼吸曲線が続き，その無呼吸区間の長さが１０秒以上となる時，正常部分呼吸曲線と無呼吸部分呼吸曲線の間の変化点を含むエポックは、無呼吸低呼吸エポックとなる。無呼吸低呼吸エポックには無呼吸部分呼吸曲線の開始点が含まれる。
【0061】
呼吸曲線に占める正常呼吸区間の割合は呼吸曲線に占める無呼吸区間の割合よりも大きい。このため，無呼吸部分呼吸曲線は正常呼吸区間に続き，その開始点が無呼吸イベントとなる場合が多い。つまり，あるクラスターが無呼吸区間に対応するならば，そのクラスターに属する部分呼吸曲線の開始点が正常呼吸エポックに含まれる割合は低く，無呼吸低呼吸エポックに含まれる割合が高くなる。各クラスターについてそのクラスターに属する部分呼吸曲線の開始点が，正常呼吸エポックに含まれる数と無呼吸低呼吸エポックに含まれる数を集計することで無呼吸区間に対応するクラスターを検討する。
【0062】
表１は，クラスターごとに，正常呼吸，無呼吸低呼吸エポックにおける部分呼吸曲線の開始点の数を示す。クラスターの部分呼吸曲線数に対する無呼吸低呼吸エポックに含まれる部分呼吸曲線の開始点の数の比が大きいクラスターは，大きい順にクラスター７，３，２である。これらのクラスターは無呼吸区間に対応するクラスターである可能性が高い。しかし，いずれのクラスターについても正常呼吸エポックに含まれる部分呼吸曲線の開始点の数の方が、無呼吸低呼吸エポックに含まれる部分呼吸曲線の開始点の数よりも大きい。
【表1】
【0063】
図１１は、呼吸曲線を分割しクラスタリングを行った結果を示す。縦線は変化点であり，縦線で区切られた呼吸曲線が部分呼吸曲線である。部分呼吸曲線の左上に示した数字はその部分呼吸曲線が属するクラスター番号を示す。ｔ＝３０，６０，．．．，１５０［ｓ］の縦線は、エポックの境界を示す。ＰＳＧ検査で検査技師に付けられたラベルをエポックの区間の下部に示す。ＯＳＡ，ＯＳＨはそのエポックにそれぞれ閉塞性無呼吸，閉塞性低呼吸のラベルが付けられたことを示す。
【0064】
図１１（ａ）ではｔ＝６０～９０［ｓ］のエポックにＯＳＨのラベルが付けられている。ｔ＝８５～１００［ｓ］の部分呼吸曲線はクラスター３に属し，振幅が小さく低呼吸区間に対応している。
図１１（ｂ）ではｔ＝３０～６０［ｓ］，ｔ＝６０～９０［ｓ］とｔ＝１５０～１８０［ｓ］のエポックに閉塞性無呼吸のラベルが付けられている。呼吸曲線の振幅から考えると，ｔ＝３５～１００［ｓ］の区間は６０秒付近の振幅の増加を除き無呼吸区間とみなすことができる。この区間はクラスター３と７に属する部分呼吸曲線に対応している。最後のエポックの無呼吸区間は最終エポックの終了時刻の５秒前から始まっている。
図１１（ｂ）ではｔ＝３０～６０［ｓ］，ｔ＝６０～９０［ｓ］とｔ＝１５０～１８０［ｓ］のエポックに閉塞性無呼吸のラベルが付けられている。呼吸曲線の振幅から考えると，ｔ＝３５～１００［ｓ］の区間は６０秒付近の振幅の増加を除き無呼吸区間とみなすことができる。この区間はクラスター３と７に属する部分呼吸曲線に対応している。最後のエポックの無呼吸区間は最終エポックの終了時刻の５秒前から始まっている。
【0065】
実施の形態では呼吸曲線を６エポックごとに分離して解析を行ったため，最後のエポックの判定に用いられた無呼吸区間のデータが解析対象に含まれていない。この次のエポックを含めて解析を行うことで，ｔ＝３０～６０［ｓ］の無呼吸イベントと同様に，クラスター３に属する部分呼吸曲線の存在を検出できるはずである。
【0066】
図１１（ｃ）は、ｔ＝６０～９０［ｓ］とｔ＝１２０～１５０［ｓ］のエポックに閉塞性低呼吸のラベルが付けられている。呼吸曲線の振幅の変化は、図１１（ａ），図１１（ｂ）よりも緩やかであるが，クラスター７に属するｔ＝７０～１００［ｓ］の二つの部分呼吸曲線が無呼吸区間に対応している。振幅変化の違いと同様に無呼吸区間に対応するクラスターが変化している。
【0067】
図１１（ｄ）において，ｔ＝１３５～１５５［ｓ］の二つの部分呼吸曲線はクラスター２に属する。ｔ＝１２０～１５０［ｓ］のエポックの閉塞性低呼吸のラベルが付けられた無呼吸区間とクラスター２に属する二つの部分呼吸曲線の区間には対応がある。クラスター２の部分呼吸曲線についても，無呼吸区間との一致を確認することができた。
【0068】
図１１と以上の議論から，クラスター２，３，７に属する部分呼吸曲線は無呼吸区間に対応していることがわかる。しかし，表１によると，クラスター２，３，７に属する部分呼吸曲線の内，半分以上は正常呼吸エポックを開始点としている。図１１から，クラスター２，３，７に属する部分呼吸曲線が実際に無呼吸区間に対応しているにもかかわらず，その開始点が正常呼吸エポックとなる場合が二つある。
【0069】
一つ目は，正常呼吸区間に続いてクラスター２，３，７に属する部分呼吸曲線が表れ，短時間で正常呼吸に戻る場合である。図１１ｂのｔ＝１４０～１５０［ｓ］の部分呼吸曲線が一つの例である。この部分呼吸曲線は前後の部分呼吸曲線に対して振幅が小さく無呼吸区間である。しかし，区間の長さが１０秒と短く無呼吸イベントの開始点として判定されない。その結果，表１では正常呼吸エポックを開始点としている部分呼吸曲線として集計される。
【0070】
二つ目は，複数の部分呼吸曲線の区間で構成された無呼吸区間がエポックをまたいで存在し，エポックの境界以降に部分呼吸曲線の開始点が存在する場合である。これは図１１ｃのｔ＝１２８［ｓ］以降の区間に見られる。この区間はクラスター７に属する部分呼吸曲線三つで構成されており，呼吸曲線の振幅から見ても無呼吸区間に対応している。無呼吸区間がエポックの境界を超えて存在する場合には，前のエポックしか無呼吸低呼吸エポックにならない。ｔ＝１２８～１５８［ｓ］の部分呼吸曲線の開始点は無呼吸低呼吸エポックに含まれる。しかし，ｔ＝１５８～１７０［ｓ］とｔ＝１７０～１８０［ｓ］の部分呼吸曲線の開始点は正常呼吸エポックとして集計される。
【0071】
表１を元に無呼吸区間に対応する可能性の高いクラスターの部分呼吸曲線区間が正常呼吸区間か無呼吸区間のいずれに対応するかを集計する。表１の正常呼吸エポックに開始する部分呼吸曲線の内，前述した部分呼吸曲線が無呼吸区間に対応する２つのパターンを無呼吸区間として集計する。つまり，正常呼吸エポックを開始点とするクラスター２，３，７に属する部分呼吸曲線について，続く部分呼吸曲線がクラスター２，３，７以外に属し部分呼吸曲線の長さが短い場合と，無呼吸低呼吸エポックを開始点とするクラスター２，３，７に含まれる部分呼吸曲線から続く場合とを、解析結果を目視で確認し、集計する。この条件に当てはまらない正常呼吸エポックに開始する部分呼吸曲線は正常呼吸区間として集計する。また，無呼吸低呼吸エポックに開始する部分呼吸曲線は無呼吸区間として集計する。
【0072】
表２は，集計結果として，クラスター２，３，７の部分呼吸曲線と正常呼吸区間と無呼吸区間の対応を示す。クラスター２，３，７のすべてについて，部分呼吸曲線の６６％以上が無呼吸区間に相当している。表２においてクラスター２の正常呼吸区間とされた部分呼吸曲線の中には，無呼吸低呼吸エポックの開始点の少し前から始まりほとんどの時区間が無呼吸低呼吸エポックに含まれるものが６個存在した。これを考慮に入れるとクラスター２の部分呼吸曲線について無呼吸区間に相当する割合は８０％に達する。以上の議論から，クラスター２，３，７は無呼吸区間に対応する部分呼吸曲線の集合であると結論づける。
【表2】
【0073】
（全体として）
実施の形態は、それぞれ組み合わせることができる。
各実施の形態はそれぞれ組み合わせることができる。
上記は無呼吸現象を抽出したが、ＣＯＶＩＤ－１９、ＳＡＲＳ（重症急性呼吸器症候群）などの呼吸系の異常も抽出できる。
実施の形態は、それぞれ組み合わせることができる。
各実施の形態はそれぞれ組み合わせることができる。
上記は無呼吸現象を抽出したが、ＣＯＶＩＤ－１９、ＳＡＲＳ（重症急性呼吸器症候群）などの呼吸系の異常も抽出できる。
【産業上の利用可能性】
【0074】
本発明の目覚まし方法およびそれを用いた機器は、時計、目覚まし時計、携帯機器、携帯電話で利用できる。
【符号の説明】
【0075】
１０ 人
１１ センサ
１２ 制御部
１１ センサ
１２ 制御部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
