(57)【要約】  生体音解析装置は、生体音の経時的な変化を示す第１生体音情報を取得する取得手段（１００）と、第１生体音情報に対して、生体音に含まれる雑音を示す第１雑音情報を強調させる第１の処理を行うことで、第２生体音情報を生成する処理手段（２３０）と、第２生体音情報の互いに隣接する期間における相関関係を示す相関情報（Ｐ（ｔ））を算出する算出手段（２４０）と、相関情報に基づいて、生体音に含まれる連続的な雑音を示す第２雑音情報を出力する出力手段（２５０）とを備える。このようにすれば、取得した生体音に含まれる連続的な雑音を好適に解析することが可能である。
2017130417.tif 000002
【特許請求の範囲】
【請求項１】
生体音の経時的な変化を示す第１生体音情報を取得する取得手段と、
前記第１生体音情報に対して、前記生体音に含まれる雑音を示す第１雑音情報を強調させる第１の処理を行うことで、第２生体音情報を生成する処理手段と、
前記第２生体音情報の互いに隣接する期間における相関関係を示す相関情報を算出する算出手段と、
前記相関情報に基づいて、前記生体音に含まれる連続的な雑音を示す第２雑音情報を出力する出力手段と
を備えることを特徴とする生体音解析装置。
【請求項２】
前記処理手段は、前記第１の処理として、前記第１生体音情報に対して逆フーリエ変換処理を行うことで算出されるケプストラム情報から、所定次数以下の成分を低減させる処理を行うことを特徴とする請求項１に記載の生体音解析装置。
【請求項３】
前記処理手段は、前記第１の処理を行う前に、所定期間内の複数の前記ケプストラム情報を平均して得られる特性に対応する成分を、前記複数のケプストラム情報の各々から低減させる第２の処理を行うことを特徴とする請求項２に記載の生体音解析装置。
【請求項４】
前記処理手段は、前記第１の処理を行った後に、前記第２生体音情報の所定基準値からの乖離量に応じて定まる情報量を所定区間ごとに算出し、前記情報量が第１閾値以下の成分を低減させる第３の処理を行うことで、第２生体音情報を生成することを特徴とする請求項３に記載の生体音解析装置。
【請求項５】
前記出力手段は、所定周期における前記相関情報が第２閾値を所定時間以上連続して超える期間の割合を算出し、前記期間の割合が所定割合以上である場合に、前記第２雑音情報を出力することを特徴とする請求項１から４のいずれか一項に記載の生体音解析装置。
【請求項６】
前記出力手段は、前記第２雑音情報のうち、所定の周波数以下の成分の割合が第３閾値以上であるものを、類鼾音を示す第３雑音情報として出力することを特徴とする請求項１から５のいずれか一項に記載の生体音解析装置。
【請求項７】
生体音の経時的な変化を示す第１生体音情報を取得する取得工程と、
前記第１生体音情報に対して、前記生体音に含まれる雑音を示す第１雑音情報を強調させる第１の処理を行うことで、第２生体音情報を生成する処理工程と、
前記第２生体音情報の互いに隣接する期間における相関関係を示す相関情報を算出する算出工程と、
前記相関情報に基づいて、前記生体音に含まれる連続的な雑音を示す第２雑音情報を出力する出力工程と
を備えることを特徴とする生体音解析方法。
【請求項８】
生体音の経時的な変化を示す第１生体音情報を取得する取得工程と、
前記第１生体音情報に対して、前記生体音に含まれる雑音を示す第１雑音情報を強調させる第１の処理を行うことで、第２生体音情報を生成する処理工程と、
前記第２生体音情報の互いに隣接する期間における相関関係を示す相関情報を算出する算出工程と、
前記相関情報に基づいて、前記生体音に含まれる連続的な雑音を示す第２雑音情報を出力する出力工程と
をコンピュータに実行させることを特徴とするコンピュータプログラム。
【請求項９】
請求項８に記載のコンピュータプログラムが記録されていることを特徴とする記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、例えば連続性ラ音等の雑音を含む生体音を解析する生体音解析装置及び生体音解析方法、並びにコンピュータプログラム及び記録媒体の技術分野に関する。
【背景技術】
【０００２】
この種の装置として、電子聴診器等によって検出される生体の呼吸音について、そこに含まれる副雑音（即ち、正常な呼吸音とは異なる音）を検出しようとするものが知られている。例えば特許文献１では、ＦＦＴ（Fast Fourier Transform：高速フーリエ変換）処理によって、振幅スペクトル、位相スペクトル、パワースペクトルを算出し、パワースペクトルの局所分散値が閾値を超えているか否かによって正常呼吸音と連続性ラ音とを分別するという技術が開示されている。また、特許文献２では、自己相関関数上のピーク位置を検出した後に、誤検出と判断できるピークを削除して、連続性ラ音を検出するという技術が開示されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００４−３５７７５８号公報
【特許文献２】特開２０１４−２８１１１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
上述した特許文献１及び２に記載されている技術では、ＦＦＴ処理を行ってからスペクトルのピークを利用して連続性ラ音を検出している。しかしながら、連続性ラ音は、個人差や測定環境の影響等によって、ピークの現れ方が一様ではない。このため、ピークのみから連続性ラ音を検出しようとすると、検出精度が低下してしまう。よって、特許文献１及び２を含む従来の技術は、生体音に含まれる連続性ラ音を正確に検出することができないという技術的問題点を有している。【０００５】
本発明が解決しようとする課題には、上記のようなものが一例として挙げられる。本発明は、生体音に含まれる雑音を好適に解析可能な生体音解析装置及び生体音解析方法、並びにコンピュータプログラム及び記録媒体を提供することを課題とする。
【課題を解決するための手段】
【０００６】
上記課題を解決するための生体音解析装置は、生体音の経時的な変化を示す第１生体音情報を取得する取得手段と、前記第１生体音情報に対して、前記生体音に含まれる雑音を示す第１雑音情報を強調させる第１の処理を行うことで、第２生体音情報を生成する処理手段と、前記第２生体音情報の互いに隣接する期間における相関関係を示す相関情報を算出する算出手段と、前記相関情報に基づいて、前記生体音に含まれる連続的な雑音を示す第２雑音情報を出力する出力手段とを備える。【０００７】
上記課題を解決するための生体音解析方法は、生体音の経時的な変化を示す第１生体音情報を取得する取得工程と、前記第１生体音情報に対して、前記生体音に含まれる雑音を示す第１雑音情報を強調させる第１の処理を行うことで、第２生体音情報を生成する処理工程と、前記第２生体音情報の互いに隣接する期間における相関関係を示す相関情報を算出する算出工程と、前記相関情報に基づいて、前記生体音に含まれる連続的な雑音を示す第２雑音情報を出力する出力工程とを備える。【０００８】
上記課題を解決するためのコンピュータプログラムは、生体音の経時的な変化を示す第１生体音情報を取得する取得工程と、前記第１生体音情報に対して、前記生体音に含まれる雑音を示す第１雑音情報を強調させる第１の処理を行うことで、第２生体音情報を生成する処理工程と、前記第２生体音情報の互いに隣接する期間における相関関係を示す相関情報を算出する算出工程と、前記相関情報に基づいて、前記生体音に含まれる連続的な雑音を示す第２雑音情報を出力する出力工程とをコンピュータに実行させる。【０００９】
上記課題を解決するための記録媒体は、上述したコンピュータプログラムが記録されている。
【図面の簡単な説明】
【００１０】
【図１】本実施例に係る生体音解析装置の構成を示すブロック図である。
【図２】本実施例に係る生体音解析装置の動作の流れを示すフローチャートである。
【図３】呼吸音の時間軸波形の一例を示す波形図である。
【図４】呼吸音のパワースペクトルの一例を示すパワースペクトログラムである。
【図５】パワースペクトルのケプストラムへの変換方法を示すプロセス図である。
【図６】ＣＭＮ処理を示すプロセス図である。
【図７】ＣＭＮ処理後のパワースペクトルの一例を示すパワースペクトログラムである。
【図８】リフタリング処理を示すプロセス図である。
【図９】リフタリング処理後のパワースペクトルの一例を示すパワースペクトログラムである。
【図１０】連続性ラ音が発生している区間における隣接フレーム間での相関係数の算出方法を示すプロセス図である。
【図１１】連続性ラ音が発生していない区間における隣接フレーム間での相関係数の算出方法を示すプロセス図である。
【図１２】連続性ラ音を含む呼吸音の相関係数の一例を示すグラフである。
【図１３】連続性ラ音を含まない呼吸音の相関係数の一例を示すグラフである。
【図１４】連続性ラ音傾向の算出方法を示す概念図である。
【図１５】ＫＬ情報量の数値化概念を示す３次元グラフ（その１）である。
【図１６】ＫＬ情報量の数値化概念を示す３次元グラフ（その２）である。
【図１７】ＫＬ情報量の計算手法を示す概念図である。
【図１８】ＫＬ情報量を用いた強調処理前のパワースペクトルの一例を示すパワースペクトログラムである。
【図１９】ＫＬ情報量を用いた強調処理後のパワースペクトルの一例を示すパワースペクトログラムである。
【図２０】類鼾音を含む生体音のパワースペクトルの一例を示すパワースペクトログラムである。
【図２１】笛声音を含む生体音のパワースペクトルの一例を示すパワースペクトログラムである。
【図２２】類鼾音及び笛声音の周波数傾向を示すグラフである。
【図２３】図２０に示すパワースペクトログラムの類鼾音傾向を算出する工程を示すプロセス図である。
【図２４】図２１に示すパワースペクトログラムの類鼾音傾向を算出する工程を示すプロセス図である。
【発明を実施するための形態】
【００１１】
＜１＞
本実施形態に係る生体音解析装置は、生体音の経時的な変化を示す第１生体音情報を取得する取得手段と、前記第１生体音情報に対して、前記生体音に含まれる雑音を示す第１雑音情報を強調させる第１の処理を行うことで、第２生体音情報を生成する処理手段と、前記第２生体音情報の互いに隣接する期間における相関関係を示す相関情報を算出する算出手段と、前記相関情報に基づいて、前記生体音に含まれる連続的な雑音を示す第２雑音情報を出力する出力手段とを備える。
【００１２】
本実施形態に係る生体音解析装置によれば、その動作時には、先ず取得手段によって、第１生体音情報が取得される。なお、「生体音」とは、生体が発する音であり、典型的には呼吸音である。また、「第１生体音情報」とは、生体音の経時的な変化を示す情報であり、例えば生体音を示す時間軸波形として取得される。
【００１３】
第１生体音情報が取得されると、処理手段によって、第１生体音情報に対する第１の処理が行われる。第１の処理は、生体音に含まれる雑音を示す第１雑音情報を強調させる処理である。なお「雑音」とは、正常な生体音には含まれない異常音（例えば、呼吸音に含まれる副雑音等）である。また第１雑音情報は、例えば雑音の有無を示す情報や、雑音の種別、雑音の強度等を示す情報である。
【００１４】
第１の処理としては、例えばＣＭＮ（Cepstral Mean Normalization）処理や、リフタリング処理（具体的には、ケプストラムの低次ケフレンシー成分をカットする処理）等の既存の処理を用いることができる。第１の処理を行うことで、第１生体音情報と比べて第１雑音情報が強調された、第２生体音情報が生成される。
【００１５】
第２生体音情報が生成されると、算出手段によって、第２生体音情報の相関情報が算出される。ここでの「相関情報」とは、第２生体音情報の互いに隣接する期間における相関関係を示す情報である。相関情報は、例えば第２生体音情報を複数のフレームに分割して、隣接するフレーム間でのケプストラム情報の相関係数として算出することができる。この場合、隣接するフレーム間の相関関係が強いほど、相関情報（相関係数）は大きな値として算出される。
【００１６】
相関情報が算出されると、出力手段によって、生体音に含まれる連続的な雑音（例えば、連続性ラ音等）を示す第２雑音情報が出力される。本実施形態では特に、第２雑音情報は、上述した相関情報に基づいて出力される。具体的には、相関情報が大きい（即ち、隣接する期間での相関関係が強い）ものが、連続的な雑音と判定され、その判定結果に応じて第２雑音情報が出力される。
【００１７】
第２雑音情報で示される連続的な雑音は、その特性上、スペクトログラム上に時間的に連続して現れるものであるため、他の音（例えば、断続性ラ音等）と比べると、隣接する期間の相関関係も強くなると考えられる。よって、相関情報に基づいて連続的な雑音を判定すれば、時間的な連続性を考慮して、より正確な第２雑音情報を出力することが可能となる。
【００１８】
＜２＞
本実施形態に係る生体音解析装置の一態様では、前記処理手段は、前記第１の処理として、前記第１生体音情報に対してフーリエ変換処理を行うことでパワースペクトル情報を算出し、パワースペクトル情報の対数をとって逆フーリエ変換処理を行うことで算出されるケプストラム情報から、所定次数以下の成分を低減させる処理を行う。
【００１９】
この態様によれば、第１の処理として、いわゆるリフタリング処理が行われる。リフタリング処理では、ケプストラム情報から所定次数以下の成分が低減されることで、雑音の特徴部分を際立たせることができる。よって、第１雑音情報が強調された第２生体音情報を好適に生成することができる。
【００２０】
＜３＞
上述したケプストラム情報の所定次数以下の成分を低減させる態様では、前記処理手段は、前記第１の処理を行う前に、所定期間内の複数の前記ケプストラム情報を平均して得られる特性に対応する成分を、前記複数のケプストラム情報の各々から低減させる第２の処理を行ってもよい。
【００２１】
この場合、第１の処理（リフタリング処理）の前に、第２の処理が行われる。第２の処理は、いわゆるＣＭＮ処理であり、ケプストラム情報から、センサや環境等の影響で定常的に畳み込まれている特性を低減することができる。この特性は、所定期間（例えば、１呼吸周期）内の複数のケプストラム情報を平均して得られる。
【００２２】
第１の処理の前に第２の処理を行えば、雑音情報以外の音（即ち、不要な成分）を低減することができるため、より正確な第２雑音情報を出力することが可能となる。
【００２３】
＜４＞
上述したケプストラム情報から平均した特性に対応する成分を低減する態様では、前記処理手段は、前記第１の処理を行った後に、前記第２生体音情報の所定基準値からの乖離量に応じて定まる情報量を所定区間ごとに算出し、前記情報量が第１閾値以下の成分を低減させる第３の処理を行うことで、第２生体音情報を生成してもよい。
【００２４】
この場合、第１の処理（リフタリング処理）の後に、第３の処理が行われることで、第２生体情報が生成される。第３の処理では、まず第２生体音情報の所定基準値からの乖離量に応じて定まる情報量が所定区間ごとに算出される。ここでの「情報量」とは、第２生体音情報における特徴部分を数値化したものであり、その一例としてＫＬ（カルバック・ライブラー）情報量が挙げられる。情報量が算出されると、情報量が第１閾値以下の成分が低減される。なお、「第１閾値」は、特徴部分ではないと判断できる程度に情報量が小さいことを判定するための閾値である。
【００２５】
上述した第３の処理によれば、情報量が小さい（特徴的ではない）成分が低減されるため、雑音を示す第１雑音情報がより強調されることになる。
【００２６】
＜５＞
本実施形態に係る生体音解析装置の他の態様では、前記出力手段は、所定周期における前記相関情報が第２閾値を所定時間以上連続して超える期間の割合を算出し、前記期間の割合が所定割合以上である場合に、前記第２雑音情報を出力する。
【００２７】
この態様によれば、出力手段では、まず所定周期（例えば、１呼吸周期）における相関情報が第２閾値を所定時間以上連続して超える期間の割合が算出される。即ち、相関情報が比較的高い状態のまま維持されていた期間の割合が算出される。なお、「第２閾値」は、相関関数が雑音と判断できる程度に高いことを判定するための閾値である。また「所定時間」とは、瞬間的な相関情報の変動を除外するために設定される時間である。
【００２８】
出力手段では更に、算出された期間の割合が所定割合以上であるか否かが判定され、所定割合上である場合に第２雑音情報が出力される（即ち、生体音が連続的な雑音を含んでいると判断される）。言い換えれば、期間の割合が所定割合未満である場合には、第２雑音情報は出力されない（即ち、生体音が連続的な雑音を含んでいないと判断される）。なお、「所定割合」は、期間の割合が雑音を検出できるまでに高いことを判定するための閾値として設定される値である。
【００２９】
上述したように、相関関数が高い状態を維持した期間の割合を利用すれば、一時的に相関関数が高くなるようなケースを除外できるため、より正確な第２雑音情報を出力することが可能となる。
【００３０】
＜６＞
本実施形態に係る生体音解析装置の他の態様では、前記出力手段は、前記第２雑音情報のうち、所定の周波数以下の成分の割合が第３閾値以上であるものを、類鼾音を示す第３雑音情報として出力する。
【００３１】
この態様によれば、出力手段では、第２雑音情報のうち、所定の周波数以下の成分の割合が第３閾値以上となるものが存在するか否かが判定される。なお、「所定の周波数」とは、類鼾音に特有の周波数成分を判別するための閾値であり、例えば１２０Ｈｚ〜１３０Ｈｚ程度の値として設定される。また、「第３閾値」は、所定の周波数以下の成分の割合が、類鼾音と判断できる程度に高いか否かを判定するために設定される閾値である。
【００３２】
類鼾音は、連続性ラ音の中でも比較的低い周波数帯成分を含むものとして知られている。よって、上述した判定によれば、連続性ラ音が類鼾音か否かを好適に判定することができる。なお、類鼾音を示す第３雑音情報は、第２雑音情報と共に出力されてもよいし、第２雑音情報に代えて出力されてもよい。
【００３３】
＜７＞
本実施形態に係る生体音解析方法は、生体音の経時的な変化を示す生体音情報を取得する第１取得工程と、前記生体音情報に含まれ得る外部騒音の特性を示す特性情報を取得する第２取得工程と、前記生体音情報及び前記特性情報に基づいて、前記生体音に含まれる雑音を示す雑音情報を出力する出力工程とを備える。
【００３４】
本実施形態に係る生体音解析方法によれば、上述した本実施形態に係る生体音解析装置と同様に、正確な雑音情報を出力することができる。
【００３５】
なお、本実施形態に係る生体音解析方法においても、上述した本実施形態に係る生体音解析装置における各種態様と同様の各種態様を採ることが可能である。
【００３６】
＜８＞
本実施形態に係るコンピュータプログラムは、生体音の経時的な変化を示す生体音情報を取得する第１取得工程と、前記生体音情報に含まれ得る外部騒音の特性を示す特性情報を取得する第２取得工程と、前記生体音情報及び前記特性情報に基づいて、前記生体音に含まれる雑音を示す雑音情報を出力する出力工程とをコンピュータに実行させる。
【００３７】
本実施形態に係るコンピュータプログラムによれば、上述した本実施形態に係る生体音解析方法と同様の処理をコンピュータに実行させることができるため、正確な雑音情報を出力することができる。
【００３８】
なお、本実施形態に係るコンピュータプログラムにおいても、上述した本実施形態に係る生体音解析装置における各種態様と同様の各種態様を採ることが可能である。
【００３９】
＜９＞
本実施形態に係る記録媒体は、上述したコンピュータプログラムが記録されている。
【００４０】
本実施形態に係る記録媒体によれば、上述したコンピュータプログラムをコンピュータにより実行させることにより、正確な雑音情報を出力することができる。
【００４１】
本実施形態に係る生体音解析装置及び生体音解析方法、並びにコンピュータプログラム及び記録媒体の作用及び他の利得については、以下に示す実施例において、より詳細に説明する。
【実施例】
【００４２】
以下では、生体音解析装置及び生体音解析方法、並びにコンピュータプログラム及び記録媒体の実施例について、図面を参照しながら詳細に説明する。なお、以下では、呼吸音の解析を行う生体音解析装置を例に挙げて説明する。【００４３】
＜第１実施例＞
第１実施例に係る生体音解析装置について、図１から図１４を参照して説明する。【００４４】
＜装置構成＞
先ず、本実施例に係る生体音解析装置の構成について、図１を参照して説明する。ここに図１は、本実施例に係る生体音解析装置の構成を示すブロック図である。【００４５】
図１において、本実施例に係る生体音解析装置は、呼吸音取得部１００と、処理部２００と、結果表示部３００とを備えて構成されている。【００４６】
呼吸音取得部１００は、生体の呼吸音を呼吸音信号として取得可能に構成されたセンサである。呼吸音取得部１００は、例えばＥＣＭ（Electret Condenser Microphone）やピエゾを利用したマイク、振動センサ等で構成されている。また、呼吸音取得部１００は、生体の呼吸音を呼吸音信号として取得可能に構成されたセンサだけでなく、センサからの呼吸音信号を取得するものを含んでいてもよい。呼吸音取得部１００で取得された呼吸音信号は、処理部２００における時間周波数解析部２１０に出力される構成となっている。なお、呼吸音取得部１００は、「取得手段」の一具体例である。【００４７】
処理部２００は、複数の演算回路やメモリ等を含んで構成されている。処理部２００は、時間周波数解析部２１０と、ケプストラム解析部２２０と、雑音強調処理部２３０と、相関情報算出部２４０と、連続性ラ音判定部２５０を備えて構成されている。【００４８】
時間周波数解析部２１０は、呼吸音取得部１００で取得された呼吸音情報に対して時間周波数解析処理を実行する。具体的には、時間周波数解析部２１０は、呼吸音情報に対してＦＦＴ処理等を実行可能に構成されている。時間周波数解析部２１０の解析結果は、ケプストラム解析部２２０に出力される構成となっている。【００４９】
ケプストラム解析部２２０は、時間周波数解析された呼吸音情報（パワースペクトル）を、ケプストラムに変換する。具体的には、ケプストラム解析部２２０は、呼吸音情報に対して逆フーリエ変換を実行可能に構成されている。ケプストラム解析部２２０の解析結果は、雑音強調処理部２３０に出力される構成となっている。【００５０】
雑音強調処理部２３０は、ケプストラム化された呼吸音情報に対し、含まれている雑音を強調させる雑音強調処理を行う。雑音強調処理は、例えばＣＭＮ処理やリフタリング処理等である。この雑音強調処理については、後の動作説明において詳しく述べる。雑音強調処理部２３０の処理結果は、相関情報算出部２４０に出力される構成となっている。雑音強調処理部２３０は、「処理手段」の一具体例である。【００５１】
相関情報算出部２４０は、雑音強調処理が実行された呼吸音情報から、相関情報を算出する。相関情報は、呼吸音情報の互いに隣接するフレーム間の相関関係を示す情報であり、フレーム間でケプストラムを比較することで算出できる。相関情報算出部２４０で算出された相関情報は、連続性ラ音判定部２４０に出力される構成となっている。相関情報算出部２４０は、「算出手段」の一具体例である。【００５２】
連続性ラ音判定部２５０は、相関情報に基づいて、呼吸音に連続性ラ音が含まれているか否かを判定する。具体的には、連続性ラ音判定部２５０は、相関情報を用いて連続性ラ音らしさを示す連続性ラ音傾向を算出し、その連続性ラ音傾向から、呼吸音に連続性ラ音が含まれているか否かを判定する。連続性ラ音判定部２５０における判定結果は、結果表示部３００に出力される構成となっている。連続性ラ音判定部２５０は、「出力手段」の一具体例である。【００５３】
以上のように、処理部２００は、呼吸音取得部１００で取得された呼吸音情報に基づいて、生体音に連続性ラ音が含まれているか否かを判定できる。また、処理部２００は、生体音に連続性ラ音が含まれているか否かだけでなく、連続性ラ音の強度等を出力可能に構成されてもよい。【００５４】
結果表示部３００は、例えば液晶モニタ等のディスプレイとして構成されており、処理部２００から出力される各種情報を画像データとして表示する。【００５５】
＜動作説明＞
次に、本実施例に係る生体音解析装置の動作について、図２を参照して説明する。ここに図２は、本実施例に係る生体音解析装置の動作の流れを示すフローチャートである。【００５６】
図２において、本実施例に係る生体音解析装置の動作時には、先ず呼吸音取得部１００において、生体の呼吸音を示す呼吸音信号が取得される（ステップＳ１０１）。呼吸音信号が取得されると、時間周波数解析部２１０において、呼吸音信号が所定の時間フレーム単位に区切られ（ステップＳ１０２）、ＦＦＴ処理による時間周波数解析が行われる（ステップＳ１０３）。【００５７】
以下では、呼吸音信号の時間周波数解析について、図３及び図４を参照して詳細に説明する。ここに図３は、呼吸音の時間軸波形の一例を示す波形図である。また図４は、呼吸音のパワースペクトルの一例を示すパワースペクトログラムである。【００５８】
図３に示すように、呼吸音情報は時間軸波形として取得される。なお、時間軸波形は、「第１生体音情報」の一具体例である。この時間軸波形を時間周波数解析してパワースペクトルにする際には、まず時間軸波形の時刻ｎの値ｘ（ｎ）に対して、長さＮのフレーム単位で、Ｎポイント短時間フーリエ変換及び対数変換を施す。これにより、時間周波数解析波形ＰＬｏｇ［ｎ，ω］が算出される。なお、フレームを切り出す際には、長さＮの窓関数ｗ（ｍ）を用いる。【００５９】
より具体的には、Ｎポイント短時間フーリエ変換は、下記数式（１）を用いて行われる。【００６０】
【数１】
2017130417.tif 000003
また、対数変換は、下記数式（２）を用いて行われる。【００６１】
【数２】
2017130417.tif 000004
図４に示すように、時間周波数解析波形ＰＬｏｇ［ｎ，ω］は、時間を横軸、周波数を縦軸とするスペクトログラムとして取得される。なお、スペクトログラムの濃淡は、パワーを示しており、明るいほどパワーが大きい。【００６２】
図２に戻り、呼吸音情報の時間周波数解析によって得られた時間周波数解析波形ＰＬｏｇ［ｎ，ω］（即ち、パワースペクトル）は、ケプストラム解析部２２０においてケプストラムに変換される（ステップＳ１０４）。【００６３】
以下では、パワースペクトルのケプストラムへの変換方法について、図５を参照して詳細に説明する。ここに図５は、パワースペクトルのケプストラムへの変換方法を示すプロセス図である。【００６４】
図５に示すように、パワースペクトルをケプストラムへの変換する際には、各時刻ｎのパワースペクトル（ログスケール）を逆フーリエ変換することで、ケプストラムＣ［ｎ，ｌ］を計算すればよい。具体的には、下記数式（３）を計算すればよい。【００６５】
【数３】
2017130417.tif 000005
なお、ｌは離散化されたケフレンシ―軸での次数であり、ｌ＝０，・・・，Ｎである。【００６６】
このようにして得られるケプストラムを利用すれば、後述する雑音強調処理、及び相関係数の算出を好適に行うことができる。【００６７】
ちなみに、ケプストラムをパワースペクトルへ逆変換する場合には、下記数式（４）を計算すればよい。【００６８】
【数４】
2017130417.tif 000006
再び図２に戻り、ケプストラムが取得されると、雑音強調処理部２３０において、呼吸音に含まれる雑音を強調するための雑音強調処理が実行される（ステップＳ１０５）。本実施例では、雑音強調処理としてＣＭＮ処理及びリフタリング処理が実行される。【００６９】
以下では、ＣＭＮ処理及びリフタリング処理について、図６から図９を参照して詳細に説明する。ここに図６は、ＣＭＮ処理を示すプロセス図であり、図７は、ＣＭＮ処理後のパワースペクトルの一例を示すパワースペクトログラムである。また図８は、リフタリング処理を示すプロセス図であり、図９は、リフタリング処理後のパワースペクトルの一例を示すパワースペクトログラムである。【００７０】
図６に示すように、雑音強調処理部２３０では、ケプストラムに対して、まずＣＭＮ処理が実行される。ＣＭＮ処理は、「第２の処理」の一具体例である。【００７１】
ＣＭＮ処理では、センサや環境等の定常的に畳み込まれている特性が、ケプストラムからキャンセル（即ち、除去）される。この特性Ｃｍｅａｎ［ｌ］は、下記数式（５）を用いて算出できる。【００７２】
【数５】
2017130417.tif 000007
算出した特性Ｃｍｅａｎ［ｌ］は、下記数式（６）のように、ケプストラムＣ［ｎ，ｌ］から差し引かれる。【００７３】
【数６】
2017130417.tif 000008
このようにして得られたＣ’［ｎ，ｌ］は、定常的に畳み込まれている特性の除去により、雑音が強調されたものとなる。【００７４】
図７に示すように、ＣＭＮ処理後のケプストラムＣ’［ｎ，ｌ］をパワースペクトルに逆変換すると、図４で示したＣＭＮ処理前のパワースペクトルと比べて、雑音部分が強調されていることが分かる。【００７５】
図８に示すように、ＣＭＮ処理後には、リフタリング処理が実行される。リフタリング処理は「第１の処理」の一具体例である。【００７６】
リフタリング処理では、ケプストラムから所定次数以下のケフレンシー成分がカットされる。この所定次数は、連続性ラ音を好適に強調できるような値に調整すればよい。なお、図８に示す例では、ケフレンシー０．１８ｍｓ以下の成分がカットされている。【００７７】
上述したリフタリング処理を実現するためには、適切なカット次数ｌｉｆｔを決定した上で、下記数式（７）を計算すればよい。【００７８】
【数７】
2017130417.tif 000009
リフタリング処理によれば、周波数特性の平坦な成分がカットされるため、比較的急峻な成分が強調されることになる。【００７９】
図９に示すように、リフタリング処理後のケプストラムＣ”［ｎ，ｌ］をパワースペクトルに逆変換すると、図７で示したリフタリング処理前のパワースペクトルと比べても、更に雑音部分が強調されていることが分かる。【００８０】
以上のように、ＣＭＮ処理及びリフタリング処理を実行すれば、呼吸音情報に含まれている雑音を強調することができる。なお、ＣＭＮ処理及びリフタリング処理によって強調される雑音情報は、「第１雑音情報」の一具体例であり、ＣＭＮ処理及びリフタリング処理後の呼吸音情報は、「第２生体音情報」の一具体例である。【００８１】
再び図２に戻り、強調処理が終了すると、全ての時間フレームに対して強調処理が終了したか否かが判定される（ステップＳ１０６）。全ての時間フレームに対して強調処理が終了していない場合には（ステップＳ１０６：ＮＯ）、他の時間フレームについて、ステップＳ１０３以降の処理が繰り返し実行される。一方、全ての時間フレームに対して強調処理が終了している場合には（ステップＳ１０６：ＹＥＳ）、相関情報算出部２４０において、相関情報が算出される（ステップＳ１０７）。【００８２】
以下では、相関情報の算出方法について、図１０から図１３を参照して詳細に説明する。ここに図１０は、連続性ラ音が発生している区間における隣接フレーム間での相関係数の算出方法を示すプロセス図であり、図１１は、連続性ラ音が発生していない区間における隣接フレーム間での相関係数の算出方法を示すプロセス図である。また図１２は、連続性ラ音を含む呼吸音の相関係数の一例を示すグラフであり、図１３は、連続性ラ音を含まない呼吸音の相関係数の一例を示すグラフである。【００８３】
図１０及び図１１に示すように、相関情報は、時間軸上で互いに隣接するフレーム（例えば、ｎフレームとｎ＋１フレーム）の相関関係を示す相関係数として算出される。【００８４】
ここで特に、スペクトル周波数軸上では、ピークの数が極めて少ない（例えば、１０２４の周波数分割数に対して、ピーク数は３である）。また、連続性ラ音の周波数は時間的に変化し得る特徴がある。隣接フレーム間でスペクトルの相関をとった場合、少数のピーク周波数の位置がわずかでも異なれば、連続的な変化であるのにも関わらず相関係数は低い値となる。これに対して、ケプストラムに変換すると、調波構造に関する情報がケフレンシー軸上全体にわたって現れてくるので、ピーク周波数が変化することによる影響が低減され、隣接フレーム間での相関の高低を的確に数値化することができる。【００８５】
具体的には、図１０に示す呼吸音の第nフレームに対応する時刻においては連続性ラ音が含まれているが、この場合には相関係数が高い値として算出される。一方で、図１１に示す呼吸音の第ｍフレームに対応する時刻においては連続性ラ音が含まれていないため、相関係数が低い値として算出される。【００８６】
図１２に示すように、連続性ラ音である類鼾音を含む呼吸音の相関係数には、連続して高い値となる部分が存在している。一方、図１３に示す連続性ラ音を含まない呼吸音の相関係数には、図１２のように連続して高い値となる部分は存在していない。以上の結果からも、相関係数が連続性ラ音を判定するパラメータとして有効であることが分かる。【００８７】
再び図２に戻り、相関情報が算出されると、全ての時間フレームに対して相関情報が算出されたか否かが判定される（ステップＳ１０８）。全ての時間フレームに対して相関情報が算出されていない場合には（ステップＳ１０８：ＮＯ）、他の時間フレームについて、ステップＳ１０７の処理が実行される。一方、全ての時間フレームに対して相関情報が算出された場合には（ステップＳ１０８：ＹＥＳ）、連続性ラ音判定部２５０において、連続性ラ音傾向Ｃ１が算出される（ステップＳ１０７）。【００８８】
以下では、連続性ラ音傾向Ｃ１の算出方法について、図１４を参照して詳細に説明する。ここに図１４は、連続性ラ音傾向の算出方法を示す概念図である。【００８９】
図１４に示すように、連続性ラ音傾向Ｃ１は、相関係数Ｐ（ｔ）が所定の閾値Ｐｔｈｒｅを所定の時間幅以上連続して超えている期間の、１呼吸周期に対する割合として算出される。なお、閾値Ｐｔｈｒｅは、「第２閾値」の一具体例であり、所定の時間幅は「所定時間」の一具体例である。【００９０】
より具体的には、連続性ラ音傾向Ｃ１を算出する際には、相関係数Ｐ（ｔ）が所定の閾値Ｐｔｈｒｅを超えている場合に“１”、所定の閾値Ｐｔｈｒｅ以下である場合に“０”となる関数ｇ（ｔ）が用いられる。連続性ラ音傾向Ｃ１は、関数ｇ（ｔ）を用いて、下記数式（８）により算出される。【００９１】
【数８】
2017130417.tif 000010
算出された連続性ラ音傾向Ｃ１は、所定閾値Ｃｔｈｒｅより大きいか否かが判定される（ステップＳ１１０）。なお、閾値Ｃｔｈｒｅは、連続性ラ音傾向Ｃ１が十分に高いか否かを判定するための閾値であり、実際の連続性ラ音のデータ等に基づいて予め設定されている。【００９２】
連続性ラ音傾向Ｃ１が所定の閾値Ｃｔｈｒｅより大きいと判定された場合（ステップＳ１１０：ＹＥＳ）、呼吸音には連続性ラ音が含まれると判定される（ステップＳ１１１）。一方で、連続性ラ音傾向Ｃ１が所定の閾値Ｃｔｈｒｅより大きくないと判定された場合（ステップＳ１１０：ＮＯ）、呼吸音には連続性ラ音が含まれないと判定される（ステップＳ１１２）。【００９３】
以上説明した一連の処理による判定結果は、結果表示部３００に出力される。これにより、結果表示部３００では、生体音に連続性ラ音が含まれるか否かを示す情報が表示されることになる。なお、生体音に連続性ラ音が含まれるか否かを示す情報は、「第２雑音情報」の一具体例である。【００９４】
＜第１実施例の効果＞
次に、これまで説明してきた第１実施例に係る生体音解析装置によって得られる技術的効果について説明する。【００９５】
第１実施例に係る生体音解析装置によれば、呼吸音情報に含まれる雑音を強調させた後、相関係数を用いて連続性ラ音傾向の有無が判定される。相関関数は、既に説明したように、時間的連続性を示すパラメータである。このため、相関関数を利用すれば、ノイズに埋もれたピーク成分や、非定常で周波数が連続的に変化しているピーク成分を検出することができる。この結果、環境等に起因するノイズの影響を好適に低減して、正確に連続性ラ音を判別することが可能となる。【００９６】
＜第２実施例＞
次に、第２実施例に係る生体音解析装置について、図１５から図１９を参照して説明する。なお、第２実施例は、既に説明した第１実施例と多くの部分で共通しており、一部の構成や動作のみが異なるものである。このため、以下では第１実施例と異なる部分についてのみ詳細に説明し、他の重複する部分については、適宜説明を省略するものとする。【００９７】
＜第１実施例との違い＞
第２実施例に係る生体音解析装置は、雑音強調処理の部分で第１実施例と異なっている。具体的には、第２実施例は、雑音強調処理として、ＣＭＮ処理及びリフタリング処理に加えて、ＫＬ情報量を用いた処理を実行可能に構成されている。【００９８】
＜ＫＬ情報量を用いた強調処理＞
以下では、ＫＬ情報量を用いた強調処理について、図１５から図１７を参照して具体的に説明する。ここに図１５及び図１６は夫々、情報量の数値化概念を示す３次元グラフであり、図１７は、情報量の計算手法を示す概念図である。【００９９】
ＫＬ情報量は、観測値Ｐと基準値Ｑ（例えば、理論値、モデル値、予測値等）とを用いて算出されるパラメータであり、基準値Ｑに対して特徴のある観測値Ｐが現れると、ＫＬ情報量は大きな値として算出される。ＫＬ情報量ＤＫＬは、下記数式（９）を用いて算出することができる。【０１００】
【数９】
2017130417.tif 000011
本実施例では、スペクトログラム上のパワー分布を、確率分布とみなしてＫＬ情報ＤＫＬを算出している。ＫＬ情報量ＤＫＬの算出には、ＣＭＮ処理とリフタリング処理が行われたスペクトルパワーＰｏｗｅｒ［ｎ，ω］が用いられる（ｎ：離散化された時刻インデックス、ω：離散化された周波数インデックス）。なお、時間周波数領域の情報を持った他の物理量（例えば、振幅スペクトル）を用いて計算することも可能である。【０１０１】
図１７に示すように、局所的なＫＬ情報量を算出するため、所定の時間方向幅ｗと、所定の周波数方向幅ｈとが事前に設定されている。そして、その枠を図の矢印に沿うようにスペクトログラム上順次走査していくことで、各時間周波数領域上のあるポイント［ｎ，ω］まわりの情報が数値化される。【０１０２】
まずは、［ｎ，ω］まわりの点［ｉ，ｊ］（ｉ＝ｎ−ｗ／２，・・・，ｎ＋ｗ／２、ｊ＝ω−ｈ／２，・・・，ω＋ｈ／２）で下記数式（１０）及び（１１）を夫々計算し、基準値Ｑ［ｉ，ｊ］及び観測値Ｐ［ｉ，ｊ］を求める。【０１０３】
【数１０】
2017130417.tif 000012
【０１０４】
【数１１】
2017130417.tif 000013
そして、算出したＱ［ｉ，ｊ］及びＰ［ｉ，ｊ］を用いて、下記数式（１２）を計算し、ＫＬ情報量ＤＫＬを算出する。これを、全てのｎ，ωについて計算すればよい。【０１０５】
【数１２】
2017130417.tif 000014
算出されたＫＬ情報量ＤＫＬは、所定の閾値Ｄｔｈｒｅと比較され、ＫＬ情報量ＤＫＬが閾値Ｄｔｈｒｅ以下となった成分がカットされる。これにより、ＫＬ情報量ＤＫＬが所定の閾値Ｄｔｈｒｅよりも大きくなる成分が強調されることになる。なお、閾値Ｄｔｈｒｅは、「第１閾値」の一具体例であり、ＫＬ情報量ＤＫＬを利用した強調処理は、「第３の処理」の一具体例である。【０１０６】
＜第２実施例の効果＞
次に、第２実施例に係る生体音解析装置によって得られる技術的効果について、図１８及び図１９を参照して詳細に説明する。ここに図１８は、ＫＬ情報量を用いた強調処理前のパワースペクトルの一例を示すパワースペクトログラムであり、図１９は、ＫＬ情報量を用いた強調処理後のパワースペクトルの一例を示すパワースペクトログラムである。【０１０７】
図１８において、第１実施例でも説明したように、ＣＭＮ処理及びリフタリング処理を行えば、呼吸音情報に含まれる雑音を好適に強調させることができる。しかしながら、図に示す例のように、ＣＭＮ処理及びリフタリング処理行っただけでは、検出すべき成分（即ち、連続性ラ音成分）だけでなく、不要な成分が残ってしまうことがある。このように不要な成分が検出されていると、相関係数を算出する上で悪影響が生じてしまうと考えられる。【０１０８】
これに対し、上述したＫＬ情報量ＤＫＬを用いた強調処理を行えば、不要な成分を更に低減することが可能である。即ち、ＫＬ情報量ＤＫＬを用いた強調処理を行えば、ＣＭＮ処理及びリフタリング処理だけを行う場合と比較して、連続性ラ音に対する検出感度を向上させることができる。【０１０９】
図１９において、ＫＬ情報量ＤＫＬを用いた強調処理後のパワースペクトルでは、図１８で残されていた不要な成分が除去されている。このように、検出すべき成分だけを残すことができれば、強調処理の結果から、より適切な相関係数を算出することができる。【０１１０】
従って、第２実施例に係る生体音解析装置よれば、呼吸音に含まれる連続性ラ音をより正確に判別することが可能である。【０１１１】
＜第３実施例＞
次に、第３実施例に係る生体音解析装置について、図２０から図２４を参照して説明する。なお、第３実施例は、既に説明した第１及び第２実施例と多くの部分で共通しており、一部の構成や動作のみが異なるものである。このため、以下では第１実施例と異なる部分についてのみ詳細に説明し、他の重複する部分については、適宜説明を省略するものとする。【０１１２】
＜第１及び第２実施例との違い＞
第１及び第２実施例に係る生体音解析装置は、呼吸音に連続性ラ音が含まれているか否かを最終的な出力としている。これに対して、第３実施例に係る生体音解析装置は、呼吸音に含まれている連続性ラ音が類鼾音であるか否かを出力可能に構成されている。【０１１３】
＜類鼾音及び笛声音の特性＞
以下では、連続性ラ音である類鼾音及び笛声音の特性について、図２０から図２２を参照して詳細に説明する。ここに図２０は、類鼾音を含む生体音のパワースペクトルの一例を示すパワースペクトログラムであり、図２１は、笛声音を含む生体音のパワースペクトルの一例を示すパワースペクトログラムである。また図２２は、類鼾音及び笛声音の周波数傾向を示すグラフである。【０１１４】
図２０及び図２１を比べると分かるように、類鼾音と笛声音とでは、発生する周波数帯域が明確に異なっている。特に、１２０Ｈｚ〜１３０Ｈｚを閾値として見ると、類鼾音は閾値以下の部分でも多く発生しているのに対し、笛声音は閾値以下の部分には殆ど発生していない。なお、図２０及び図２１に示すスペクトログラムはあくまで一例であるが、本願発明者の研究によれば、その他の例においても同様の傾向があることが判明している。【０１１５】
図２２に示すように、１２０Ｈｚ〜１３０Ｈｚの閾値以下には類鼾音だけが含まれ、笛声音は含まれない。この特性を利用すれば、連続性ラ音が類鼾音であるのか、それとも笛声音であるのかを判別することができる。類鼾音を適切に判別するためには、以下で説明する類鼾音傾向Ｒを用いることが好ましい。【０１１６】
＜類鼾音傾向の算出＞
類鼾音傾向Ｒの算出方法について、図２３及び図２４を参照して詳細に説明する。ここに図２３は、図２０に示すパワースペクトログラムの類鼾音傾向を算出する工程を示すプロセス図である。また図２４は、図２１に示すパワースペクトログラムの類鼾音傾向を算出する工程を示すプロセス図である。【０１１７】
図２３及び図２４に示すように、類鼾音傾向Ｒを算出する際には、まず隣接するフレーム間の相関係数Ｐ（ｔ）が所定の閾値Ｐｔｈｒｅ以上となっている区間が、連続性ラ音発生区間として検出される。なお、連続性ラ音発生区間は複数の区間として検出されてもよい。【０１１８】
そして、連続性ラ音発生区間における１２０Ｈｚ以下のパワーの割合が類鼾音傾向Ｒとして算出される。算出された類鼾音傾向Ｒは、所定の閾値Ｒｔｈｒｅと比較され、類鼾音傾向Ｒが所定の閾値Ｒｔｈｒｅ以上である場合には、連続性ラ音が類鼾音であると判定される。一方、類鼾音傾向Ｒが所定の閾値Ｒｔｈｒｅ未満である場合には、連続性ラ音が類鼾音でない（例えば、笛声音である）と判定される。【０１１９】
＜第３実施例の効果＞
最後に、第３実施例に係る生体音解析装置によって得られる技術的効果について説明する。【０１２０】
第３実施例に係る生体音解析装置によれば、類鼾音の周波数特性を利用して、連続性ラ音が類鼾音であるか否かを判別することができる。また、本実施例では特に、類鼾音であるか否かを判別するために用いる類鼾音傾向Ｒが、相関係数Ｐ（ｔ）を用いて算出される。このため、時間的連続性が考慮されることになり、より正確な類鼾音傾向Ｒが算出される。【０１２１】
以上のように、第３実施例に係る生体音解析装置によれば、呼吸音に含まれる連続性ラ音について、より詳細で正確な解析が行える。【０１２２】
本発明は、上述した実施形態に限られるものではなく、特許請求の範囲及び明細書全体から読み取れる発明の要旨或いは思想に反しない範囲で適宜変更可能であり、そのような変更を伴う生体音解析装置及び生体音解析方法、並びにコンピュータプログラム及び記録媒体もまた本発明の技術的範囲に含まれるものである。
【符号の説明】
【０１２３】
１００  呼吸音取得部
２００  処理部
２１０  時間周波数解析部
２２０  ケプストラム解析部
２３０  雑音強調処理部
２４０  相関情報算出部
２５０  連続性ラ音判定部
３００  結果表示部
Ｐ（t）  相関係数
Ｃ１  連続性ラ音傾向
ＤＫＬ  ＫＬ情報量
Ｒ  類鼾音傾向
【図１】
2017130417.tif 000015
【図２】
2017130417.tif 000016
【図３】
2017130417.tif 000017
【図４】
2017130417.tif 000018
【図５】
2017130417.tif 000019
【図６】
2017130417.tif 000020
【図７】
2017130417.tif 000021
【図８】
2017130417.tif 000022
【図９】
2017130417.tif 000023
【図１０】
2017130417.tif 000024
【図１１】
2017130417.tif 000025
【図１２】
2017130417.tif 000026
【図１３】
2017130417.tif 000027
【図１４】
2017130417.tif 000028
【図１５】
2017130417.tif 000029
【図１６】
2017130417.tif 000030
【図１７】
2017130417.tif 000031
【図１８】
2017130417.tif 000032
【図１９】
2017130417.tif 000033
【図２０】
2017130417.tif 000034
【図２１】
2017130417.tif 000035
【図２２】
2017130417.tif 000036
【図２３】
2017130417.tif 000037
【図２４】
2017130417.tif 000038
【国際調査報告】
2017130417.tif 000039
2017130417.tif 000040
2017130417.tif 000041
2017130417.tif 000042