JP2009106574A - 呼吸音解析装置及びプログラム - Google Patents

Abstract

【課題】ラ音の検出を精度良く行うことができる呼吸音解析装置を提供する。
【解決手段】呼吸音データをヒルベルト変換して解析信号を算出し、解析信号から包絡線、瞬時周波数及び瞬時帯域幅を算出し、包絡線からラ音の候補となるピークを抽出し、抽出したピークに対応して瞬時周波数又は瞬時帯域幅のうち少なくとも１つを算出し、算出した瞬時周波数又は瞬時帯域幅のうち少なくとも１つに基づいてラ音を検出する。
【選択図】図４

Description

本発明は、呼吸音解析装置及びプログラムに関する。

呼吸音を聴取し、診断することは聴診として古くから用いられてきた診断方法である。呼吸音の中の異常音は副雑音と呼ばれ、様々な種類の音が知られている。副雑音の中でもラ音は重要で、その自動的な検出はこれまでにも試みられている（例えば、特許文献１参照）。
特表２００１−５０５０８５号公報

呼吸信号にはさまざまな、ノイズ信号が含まれる。例えば、音声や咳、心音、マイクを擦る摩擦音など非常に多岐に渡る。このノイズをいかに誤検出せず、かつ目的とするラ音を精度よく検出するかが、大きな課題であった。

ラ音の中でも断続性ラ音（クラックル）の検出は難しく、特にコースクラックル（水泡音）の検出精度が低い問題があった。

本発明は、以上のような問題に鑑みてなされたものであり、ラ音の検出を精度良く行うことができる呼吸音解析装置を提供することを目的としている。

本発明の呼吸音解析装置は、被検者の呼吸音を測定し呼吸音データを出力する呼吸音測定部と、前記呼吸音測定部により測定された呼吸音データに基づいて呼吸音の解析を行う呼吸音解析部と、を有し、前記呼吸音解析部は、前記呼吸音データをヒルベルト変換して解析信号を算出し、解析信号から包絡線、瞬時周波数及び瞬時帯域幅を算出し、包絡線からラ音の候補となるピークを抽出し、抽出したピークに対応して瞬時周波数又は瞬時帯域幅のうち少なくとも１つを算出し、算出した瞬時周波数又は瞬時帯域幅のうち少なくとも１つに基づいてラ音を検出することを特徴としている。

本発明のプログラムは、入力された呼吸音データをヒルベルト変換して解析信号を算出するステップと、解析信号から包絡線、瞬時周波数及び瞬時帯域幅を算出するステップと、包絡線からラ音の候補となるピークを抽出するステップと、抽出したピークに対応して瞬時周波数又は瞬時帯域幅のうち少なくとも１つを算出するステップと、算出した瞬時周波数又は瞬時帯域幅のうち少なくとも１つに基づいてラ音を検出するステップと、をコンピュータに実行させることを特徴としている。

本発明によれば、瞬時周波数、瞬時帯域幅を用いることにより、ラ音の検出を精度良く行うことができる。

以下、図面に基づいて本発明の実施形態について説明するが、一例であり、本実施形態に限定するものではない。

（第１の実施形態）
［装置構成］
図１は、本実施形態に係る呼吸音解析装置の構成図である。呼吸音解析装置１は、呼吸音測定部１０及び呼吸音解析部３０から構成され、通信媒体Ｌを介して互いに接続されている。この通信媒体Ｌは、有線であっても無線であってもよい。無線により接続すれば、呼吸音測定部１０に用いられているセンサに不要な振動を与える機会を減らすことができるとともに、被検者の行動の制約を少なくすることができる。

呼吸音測定部１０は、入力される呼吸音を呼吸音信号に変換するコンデンサマイク等の密着型のマイクロフォン１１、マイクロフォン１１により変換された呼吸音信号を増幅する増幅器１２、増幅器１２により増幅された呼吸音信号を平滑化処理する平滑回路１３、平滑回路１３により処理された呼吸音信号を呼吸音データに変換するＡ／Ｄ変換器１４、及びＡ／Ｄ変換器１４により変換された呼吸音データを呼吸音解析部３０に送信するＩ／Ｆ１５から構成されている。ここで、Ａ／Ｄ変換器１４は、例えばサンプリング周波数Ｆｓ＝１１．０２５ｋＨｚで呼吸音信号から呼吸音データを生成する。

呼吸音解析部３０は、呼吸音測定部１０からの呼吸音データを受信するＩ／Ｆ３１、Ｉ／Ｆ３１で受信された呼吸音データをプログラムに従って処理するＣＰＵ３２、ＣＰＵ３２での処理に必要なプログラムやデータ等を記憶するＲＯＭ３３、ＣＰＵ３２での処理に必要なプログラムやデータ等を一時的に記憶するＲＡＭ３４、ＣＰＵ３２での処理結果等をハードディスク、ＤＶＤ−Ｒ、ＣＤ−Ｒ等に保存する外部記憶装置３５、各種データを入力する入力部３６、及びＣＰＵ３２での処理結果等を表示する表示部３７から構成されている。

呼吸音解析部３０は、専用の情報処理装置で構成されていてもよいし、汎用のパーソナルコンピュータで構成されていてもよい。パーソナルコンピュータであれば、持ち運びが容易にできる携帯情報端末（ＰＤＡ）であることが好ましい。

本実施形態においては、呼吸音測定部１０にＡ／Ｄ変換器１４を設けたが、呼吸音解析部３０にＡ／Ｄ変換器を設けてもよい。

［呼吸音測定処理］
図２は、本実施形態に係る呼吸音測定処理のフロー図である。この呼吸音測定処理フローは、ＲＯＭ３３内の呼吸音測定プログラムに基づいて、ＣＰＵ３２により実行されるフローである。予め入力部３６により被検者を特定するＩＤ等は入力されているものとする。

まず、ＣＰＵ３２は、入力部３６から呼吸音測定の開始が指示されたか否かを判断する（ステップＳ１０）。呼吸音測定の開始が指示されたと判断すると（ステップＳ１０；Ｙｅｓ）、ＣＰＵ３２は、呼吸音測定部１０からの呼吸音データを外部記憶装置３５に保存する動作を開始させる（ステップＳ４１）。このとき、被検者のＩＤ及び時刻に対応付けて呼吸音データは保存される。呼吸音測定の開始が指示されていないと判断すると（ステップＳ１０；Ｎｏ）、ステップＳ１０に戻り呼吸音測定の開始が指示されるまで待機する。

次に、ＣＰＵ３２は、呼吸音測定の終了が指示されたか否かを判断する（ステップＳ１２）。呼吸音測定の終了が指示されたと判断すると（ステップＳ１２；Ｙｅｓ）、ＣＰＵ３２は、呼吸音測定部１０からの呼吸音デーを外部記憶装置３５に保存する動作を終了させる（ステップＳ１３）。呼吸音測定の終了が指示されていないと判断すると（ステップＳ１２；Ｎｏ）、ステップＳ１２に戻り呼吸音測定の終了が指示されるまで待機する。

［解析パラメータ］
まず、本発明の解析に用いるパラメータについて説明する。

解析信号ｚ（ｔ）は複素数の信号であり、その実部は原信号と等しく、虚部は原信号のヒルベルト変換である。原信号をｘ（ｔ）とし、ｘ（ｔ）の解析信号をｚ（ｔ）とすると、

と表される。但し、

であり、Ａ（ｔ）、ψ（ｔ）は、それぞれｘ（ｔ）の包絡線及び瞬時位相を示す。

瞬時周波数ｆ（ｔ）は、

で表される。

瞬時帯域幅ｂ（ｔ）は、

で表される。

［ラ音検出の考え方］
ラ音の注目すべき特徴は時間的局在性である。本実施形態では、エネルギー集中指数ＥＣＩ（Ｅｎｅｒｇｙ Ｃｏｎｃｅｎｔｒａｔｉｏｎ Ｉｎｄｅｘ）という指標を用いてラ音を検出する。

エネルギー集中指数ＥＣＩは、包絡線の最大値（極大値）の近傍における原波形のエネルギーＥ１と、包絡線の最大値を中心としてラ音１つ分程度の時間幅における原波形のエネルギーＥ２との比（Ｅ１／Ｅ２）で表される。エネルギー集中指数ＥＣＩは、波形のエネルギーがどのくらい最大値付近に集中しているかの指標となり、ラ音の指標とできる。

本アルゴリズムでは、包絡線の最大値付近における瞬時周波数の平均値の逆数をＴとして（Ｔは、そのパルス波形の振動の周期の推定値となる）、Ｅ１の値はＴの２倍の時間幅に含まれるエネルギー（データの２乗和）、Ｅ２の値はＴの６倍の時間幅に含まれるエネルギーとした。

瞬時周波数は、ラ音部分ではラ音の（主要）周波数を安定的に示す。信号のない部分では小刻みにばらつく。瞬時帯域幅は、ラ音のピークにおいてほぼゼロになる。ラ音部分で下に凸のお椀型の変化を示す。信号のない部分では大きくばらつく。

本実施形態では、ラ音の候補点において、エネルギー集中指数ＥＣＩがある閾値以上（波形の時間的局在性が大きい）、瞬時周波数の変動が少ない、瞬時帯域幅が小さい、及び候補点の近傍における瞬時帯域幅の積分値（データの総和）（以下、ＢＳと称する）がある閾値以上（瞬時帯域幅の変化がお椀型であるか否かの判断）の場合、ラ音であると判定する。

［ラ音検出における呼吸音解析処理］
図３は、第１の実施形態に係る呼吸音解析処理のフロー図である。

本呼吸音解析処理は、取得した呼吸音データを解析し、ラ音の検出を行う処理である。ＲＯＭ３３内の呼吸音解析プログラム（本発明のプログラム）に基づいて、ＣＰＵ３２により実行される。尚、予め入力部３６により被検者を特定するＩＤ等は入力されているものとする。

まず、ＣＰＵ３２は、入力部３６から呼吸音解析の指示が入力されたか否かを判断する（ステップＳ２０）。呼吸音解析の指示が入力されたと判断すると（ステップＳ２０；Ｙｅｓ）、ＣＰＵ３２は、被検者ＩＤに対応する呼吸音データを外部記憶装置３５から読み出し、ＲＡＭ３４にロードする（ステップＳ２１）。呼吸音解析の指示が入力されていないと判断すると（ステップＳ２０；Ｎｏ）、ステップＳ２０に戻り呼吸音解析の指示が入力されるまで待機する。

次に、ＣＰＵ３２は、ＲＡＭ３４にロードされた呼吸音データ（サンプリング周波数Ｆｓ＝１１．０２５ｋＨｚ）を１３１０７２点毎に区分する（ステップＳ２２）。

次に、ＣＰＵ３２は、各区間で心音等の雑音を除去するために、カットオフ周波数５０Ｈｚのハイパスフィルタをかける（ステップＳ２３）。

次に、ＣＰＵ３２は、データをフーリエ変換し、直流成分を１倍、正の周波数成分を２倍、負の周波数成分をゼロとして、逆フーリエ変換を行うことにより解析信号ｚ（ｔ）を得る（ステップＳ２４）。フーリエ変換を用いる方法のほかにヒルベルト変換器（ディジタルフィルタ）を用いる方法もある。

次に、ＣＰＵ３２は、包絡線Ａ（ｔ）、瞬時周波数ｆ（ｔ）、瞬時帯域幅ｂ（ｔ）を前掲の式により算出する。微分は差分に置き換えた。それぞれの結果を図４に示す。図４（ａ）は原信号ｘ（ｔ）、図４（ｂ）は包絡線ａ（ｔ）、図４（ｃ）はエネルギー集中指数ＥＣＩ、図４（ｄ）は瞬時周波数ｆ（ｔ）、図４（ｅ）は瞬時帯域幅ｂ（ｔ）である。さらにｆ（ｔ）の時間微分（差分）をとって２乗した量ｖ（ｔ）を算出する（ステップＳ２５）。

次に、ＣＰＵ３２は、重みを３角形型、幅を５５点としてＡ（ｔ）の移動平均Ａｓ（ｔ）を算出する。同様にｂ（ｔ）及びｖ（ｔ）についても、幅を３５点としてそれぞれの移動平均ｂｓ（ｔ）及びｖｓ（ｔ）を算出する（ステップＳ２６）。

次に、ＣＰＵ３２は、移動平均Ａｓ（ｔ）が極大となる位置を探してラ音の位置の候補とする（ステップＳ２７）。

次に、ＣＰＵ３２は、それぞれの候補位置において、エネルギー集中指数ＥＣＩを算出する。候補位置を中心として、３２点のデータ幅におけるｆ（ｔ）の平均値ＦＡを求める。ＥＣＩ＝Ｅ１／Ｅ２とし、Ｅ１を求める際のデータ数Ｎ１＝２×Ｆｓ／ＦＡ、Ｅ２を求める際のデータ数Ｎ２＝６Ｆｓ／ＦＡとする。原信号ｘ（ｔ）において、候補位置を中心としてＮ１点の幅におけるｘ（ｔ）の２乗和をＥ１とし、候補位置を中心としてＮ２点の幅におけるｘ（ｔ）の２乗和をＥ２としてエネルギー集中指数ＥＣＩを算出する（ステップＳ２８）。

次に、ＣＰＵ３２は、候補位置において、瞬時帯域幅の積分値ＢＳを計算する。瞬時帯域幅ｂ（ｔ）において、候補位置を中心として、４４２点の幅におけるｂ（ｔ）の総和を求め、ＢＳとする（ステップＳ２９）。

ラ音の候補位置において、エネルギー集中指数ＥＣＩが０．６５以上、ｖｓ（ｔ）が３．０４×１０¹⁰以下、瞬時帯域幅の移動平均ｂｓ（ｔ）が１５４以下、候補点の近傍における瞬時帯域幅の積分値ＢＳが１２以上、である場合、その位置の波形をラ音と決定する（ステップＳ３０）。図４においては、図４（ａ）の原信号ｘ（ｔ）上に３角形の印で示した。

以上のように、本実施形態によれば、瞬時周波数、瞬時帯域幅を用いることにより、ラ音の検出を精度良く行うことができる。このため、ラ音のない部分（ノイズ部分）を判別することができ、入力信号の振幅変化にロバストなアルゴリズムが実現できる。また、瞬時周波数を用いて、ラ音の（主要）周波数などを容易に求めることが可能である。

本実施形態では、瞬時帯域幅の変化がお椀型であるか否かの判断も組み込まれており、精度を高める上で好ましい。

（第２の実施形態）
装置構成及び呼吸音測定処理については、第１の実施形態と同様であり、説明を省略する。

［ラ音検出の考え方］
第１の実施形態ではエネルギー集中指数ＥＣＩを用いてラ音の検出を行ったが、第２の実施形態においては、後述するクレストインデックス又は包絡線のピーク値まわりの２次モーメントを用いてラ音の検出を行う。

クレストインデックスＣＩ（波高指数）は、波高率と同じであり、包絡線の最大値を現波形の実効値で除した値と定義する。波形がパルスに近いほど大きな値を示す。

２次モーメントは、確率における分散と同じであり、分布の広がりの程度を示す。ラ音の包絡線では、相対的に低い値を示す。

本実施形態では、ラ音の候補点において、クレストインデックスＣＩがある閾値以上又は２次モーメントの値がある閾値以下（波形の時間的局在性が大きい）、瞬時周波数の変動が少ない、及び瞬時帯域幅が小さい場合、ラ音であると判定する。

［ラ音検出における呼吸音解析処理］
図５は、第２の実施形態に係る呼吸音解析処理のフロー図である。一例としてクレストインデックスＣＩを用いてラ音の検出する場合について説明するが、２次モーメントを用いる場合も基本的に同様である。

本呼吸音解析処理は、取得した呼吸音データを解析し、ラ音の検出を行う処理である。ＲＯＭ３３内の呼吸音解析プログラムに基づいて、ＣＰＵ３２により実行される。尚、予め入力部３６により被検者を特定するＩＤ等は入力されているものとする。

まず、ＣＰＵ３２は、入力部３６から呼吸音解析の指示が入力されたか否かを判断する（ステップＳ４０）。呼吸音解析の指示が入力されたと判断すると（ステップＳ４０；Ｙｅｓ）、ＣＰＵ３２は、被検者ＩＤに対応する呼吸音データを外部記憶装置３５から読み出し、ＲＡＭ３４にロードする（ステップＳ４１）。呼吸音解析の指示が入力されていないと判断すると（ステップＳ４０；Ｎｏ）、ステップＳ４０に戻り呼吸音解析の指示が入力されるまで待機する。

次に、ＣＰＵ３２は、ＲＡＭ３４にロードされた呼吸音データ（サンプリング周波数Ｆｓ＝１１．０２５ｋＨｚ）を１３１０７２点毎に区分する（ステップＳ４２）。

次に、ＣＰＵ３２は、各区間で心音等の雑音を除去するために、カットオフ周波数５０Ｈｚのハイパスフィルタをかける（ステップＳ４３）。

次に、ＣＰＵ３２は、データをフーリエ変換し、直流成分を１倍、正の周波数成分を２倍、負の周波数成分をゼロとして、逆フーリエ変換を行うことにより解析信号ｚ（ｔ）を得る（ステップＳ４４）。フーリエ変換を用いる方法のほかにヒルベルト変換器（ディジタルフィルタ）を用いる方法もある。

次に、ＣＰＵ３２は、包絡線Ａ（ｔ）、瞬時周波数ｆ（ｔ）、瞬時帯域幅ｂ（ｔ）を前掲の式により算出する。微分は差分に置き換えた。さらにｆ（ｔ）の時間微分（差分）をとって２乗した量ｖ（ｔ）を算出する（ステップＳ４５）。

次に、ＣＰＵ３２は、重みを３角形型、幅を５５点としてＡ（ｔ）の移動平均Ａｓ（ｔ）を算出する。同様にｂ（ｔ）及びｖ（ｔ）についても、幅を３５点としてそれぞれの移動平均ｂｓ（ｔ）及びｖｓ（ｔ）を算出する（ステップＳ４６）。

次に、ＣＰＵ３２は、移動平均Ａｓ（ｔ）が極大となる位置を探してラ音の位置の候補とする。但し、移動平均Ａｓ（ｔ）の平均値を求めて、その３．２倍以下のピークは候補から除く（ステップＳ４７）。

次に、ＣＰＵ３２は、それぞれのピークのクレストインデックスＣＩを算出する。クレストインデックスＣＩを算出するデータ長は、波形のピーク位置（ラ音の候補）を中心とする１０２４点の幅とする（ステップＳ４８）。

ラ音の候補位置において、クレストインデックスＣＩが２．６以上、瞬時帯域幅の移動平均ｂｓ（ｔ）が１５４以下、及びｖｓ（ｔ）が３．０４×１０¹⁰以下である場合、その位置の波形をラ音と決定する（ステップＳ４９）。

以上のように、本実施形態によれば、瞬時周波数、瞬時帯域幅を用いることにより、ラ音の検出を精度良く行うことができる。このため、ラ音のない部分（ノイズ部分）を判別することができ、入力信号の振幅変化にロバストなアルゴリズムが実現できる。また、瞬時周波数を用いて、ラ音の（主要）周波数などを容易に求めることが可能である。

クレストインデックスもしくは包絡線の２次モーメントを計算する範囲（データ長）は、固定長でも良いが、波形のピークおよびその近傍における瞬時周波数をもとに、データ長をそれぞれのピークごとに決定することによって（そのピークの瞬時周波数の逆数の定数倍など）、検出能力を改善することができる。また瞬時周波数を用いてラ音を選択的に検出する、ノイズによる誤検出を減らすなどの処理を行うことができる。

本発明における瞬時周波数の波形の粗さに関わる指標とは、例えば、信号のある区間の偏差等のものをいう。偏差が少ないほど瞬時周波数の変動が少なくラ音と判定されやすくなる。

本発明における瞬時帯域幅の波形概形とは、ローフィルターを用いた信号の低周波成分でもよいし、区間の移動平均でもよいし、包絡線的なものでもよい。

瞬時周波数が振幅の多い波形の場合には、包絡線にするラ音でない部分にも下に凸の領域が現れる場合もあるが、この場合には、瞬時周波数の波形は粗くなるので、ラ音としては検出されない。本発明の検出方法は単独で用いても良いし、他の手法と組み合わせて用いることも好ましい。

ファインクラックル（高音性断続性ラ音）の検出については、適当なカットオフ周波数のハイパスフィルタを通して波形整形を行った後、本呼吸音解析処理を行うことによって検出精度を向上させることができる。

本発明に用いる呼吸音信号は、例えば電子聴診器により採取することができる。採取装置は電子聴診器のように、耳で聞く機能が必ずしも必要ではなく、体表面から体内の音信号を採取できるものであればよい。

本実施形態に係る呼吸音解析装置の構成図である。 本実施形態に係る呼吸音測定処理のフロー図である。 第１の実施形態に係る呼吸音解析処理のフロー図である。 第１の実施形態に係る各波形を示す概念図である。 第２の実施形態に係る呼吸音解析処理のフロー図である。

符号の説明

１０ 呼吸音測定部
３０ 呼吸音解析部
３２ ＣＰＵ
３３ ＲＯＭ

Claims (8)

1. 被検者の呼吸音を測定し呼吸音データを出力する呼吸音測定部と、
   前記呼吸音測定部により測定された呼吸音データに基づいて呼吸音の解析を行う呼吸音解析部と、
   を有し、
   前記呼吸音解析部は、
   前記呼吸音データをヒルベルト変換して解析信号を算出し、
   解析信号から包絡線、瞬時周波数及び瞬時帯域幅を算出し、
   包絡線からラ音の候補となるピークを抽出し、
   抽出したピークに対応して瞬時周波数又は瞬時帯域幅のうち少なくとも１つを算出し、
   算出した瞬時周波数又は瞬時帯域幅のうち少なくとも１つに基づいてラ音を検出することを特徴とする呼吸音解析装置。
2. 前記呼吸音解析部は、
   抽出したピークに対応してエネルギー集中指数、クレストインデックス又は２次モーメントのうち少なくとも１つを算出し、
   算出したエネルギー集中指数、クレストインデックス又は２次モーメントのうち少なくとも１つを含めてラ音を検出することを特徴とする請求項１に記載の呼吸音解析装置。
3. 前記瞬時周波数における波形の粗さに係る指標に基づいてラ音を検出することを特徴とする請求項１又は２に記載の呼吸音解析装置。
4. 前記瞬時帯域幅の波形概形が下に凸であることを示す指標に基づいてラ音を検出することを特徴とする請求項１〜３の何れか１項に記載の呼吸音解析装置。
5. 前記呼吸音解析部は、前記ラ音の候補として抽出したピークに対応する瞬時帯域幅の積分値を算出し、算出された瞬時帯域幅の積分値を含めてラ音を検出することを特徴とする請求項１〜４の何れか１項に記載の呼吸音解析装置。
6. 前記呼吸音解析部は、前記包絡線の移動平均を算出し、算出した移動平均のピークをラ音の候補となるピークとして抽出することを特徴とする請求項１〜５の何れか１項に記載の呼吸音解析装置。
7. 前記エネルギー集中指数、前記クレストインデックス又は前記２次モーメントを算出する際のデータ長は、抽出したピークにおける瞬時周波数を基に決定することを特徴とする請求項１〜６の何れか１項に記載の呼吸音解析装置。
8. 入力された呼吸音データをヒルベルト変換して解析信号を算出するステップと、
   解析信号から包絡線、瞬時周波数及び瞬時帯域幅を算出するステップと、
   包絡線からラ音の候補となるピークを抽出するステップと、
   抽出したピークに対応して瞬時周波数又は瞬時帯域幅のうち少なくとも１つを算出するステップと、
   算出した瞬時周波数又は瞬時帯域幅のうち少なくとも１つに基づいてラ音を検出するステップと、
   をコンピュータに実行させることを特徴とするプログラム。

Images