説明

呼吸音解析装置及びプログラム

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

【発明の詳細な説明】
【技術分野】
【0001】
本発明は、呼吸音解析装置及びプログラムに関する。
【背景技術】
【0002】
呼吸音を聴取し、診断することは聴診として古くから用いられてきた診断方法である。呼吸音の中の異常音は副雑音と呼ばれ、様々な種類の音が知られている。副雑音の中でもラ音は重要で、その自動的な検出はこれまでにも試みられている(例えば、特許文献1参照)。
【特許文献1】特表2001−505085号公報
【発明の開示】
【発明が解決しようとする課題】
【0003】
呼吸信号にはさまざまな、ノイズ信号が含まれる。例えば、音声や咳、心音、マイクを擦る摩擦音など非常に多岐に渡る。このノイズをいかに誤検出せず、かつ目的とするラ音を精度よく検出するかが、大きな課題であった。
【0004】
ラ音の中でも断続性ラ音(クラックル)の検出は難しく、特にコースクラックル(水泡音)の検出精度が低い問題があった。
【0005】
本発明は、以上のような問題に鑑みてなされたものであり、ラ音の検出を精度良く行うことができる呼吸音解析装置を提供することを目的としている。
【課題を解決するための手段】
【0006】
本発明の呼吸音解析装置は、被検者の呼吸音を測定し呼吸音データを出力する呼吸音測定部と、前記呼吸音測定部により測定された呼吸音データに基づいて呼吸音の解析を行う呼吸音解析部と、を有し、前記呼吸音解析部は、前記呼吸音データをヒルベルト変換して解析信号を算出し、解析信号から包絡線、瞬時周波数及び瞬時帯域幅を算出し、包絡線からラ音の候補となるピークを抽出し、抽出したピークに対応して瞬時周波数又は瞬時帯域幅のうち少なくとも1つを算出し、算出した瞬時周波数又は瞬時帯域幅のうち少なくとも1つに基づいてラ音を検出することを特徴としている。
【0007】
本発明のプログラムは、入力された呼吸音データをヒルベルト変換して解析信号を算出するステップと、解析信号から包絡線、瞬時周波数及び瞬時帯域幅を算出するステップと、包絡線からラ音の候補となるピークを抽出するステップと、抽出したピークに対応して瞬時周波数又は瞬時帯域幅のうち少なくとも1つを算出するステップと、算出した瞬時周波数又は瞬時帯域幅のうち少なくとも1つに基づいてラ音を検出するステップと、をコンピュータに実行させることを特徴としている。
【発明の効果】
【0008】
本発明によれば、瞬時周波数、瞬時帯域幅を用いることにより、ラ音の検出を精度良く行うことができる。
【発明を実施するための最良の形態】
【0009】
以下、図面に基づいて本発明の実施形態について説明するが、一例であり、本実施形態に限定するものではない。
【0010】
(第1の実施形態)
[装置構成]
図1は、本実施形態に係る呼吸音解析装置の構成図である。呼吸音解析装置1は、呼吸音測定部10及び呼吸音解析部30から構成され、通信媒体Lを介して互いに接続されている。この通信媒体Lは、有線であっても無線であってもよい。無線により接続すれば、呼吸音測定部10に用いられているセンサに不要な振動を与える機会を減らすことができるとともに、被検者の行動の制約を少なくすることができる。
【0011】
呼吸音測定部10は、入力される呼吸音を呼吸音信号に変換するコンデンサマイク等の密着型のマイクロフォン11、マイクロフォン11により変換された呼吸音信号を増幅する増幅器12、増幅器12により増幅された呼吸音信号を平滑化処理する平滑回路13、平滑回路13により処理された呼吸音信号を呼吸音データに変換するA/D変換器14、及びA/D変換器14により変換された呼吸音データを呼吸音解析部30に送信するI/F15から構成されている。ここで、A/D変換器14は、例えばサンプリング周波数Fs=11.025kHzで呼吸音信号から呼吸音データを生成する。
【0012】
呼吸音解析部30は、呼吸音測定部10からの呼吸音データを受信するI/F31、I/F31で受信された呼吸音データをプログラムに従って処理するCPU32、CPU32での処理に必要なプログラムやデータ等を記憶するROM33、CPU32での処理に必要なプログラムやデータ等を一時的に記憶するRAM34、CPU32での処理結果等をハードディスク、DVD−R、CD−R等に保存する外部記憶装置35、各種データを入力する入力部36、及びCPU32での処理結果等を表示する表示部37から構成されている。
【0013】
呼吸音解析部30は、専用の情報処理装置で構成されていてもよいし、汎用のパーソナルコンピュータで構成されていてもよい。パーソナルコンピュータであれば、持ち運びが容易にできる携帯情報端末(PDA)であることが好ましい。
【0014】
本実施形態においては、呼吸音測定部10にA/D変換器14を設けたが、呼吸音解析部30にA/D変換器を設けてもよい。
【0015】
[呼吸音測定処理]
図2は、本実施形態に係る呼吸音測定処理のフロー図である。この呼吸音測定処理フローは、ROM33内の呼吸音測定プログラムに基づいて、CPU32により実行されるフローである。予め入力部36により被検者を特定するID等は入力されているものとする。
【0016】
まず、CPU32は、入力部36から呼吸音測定の開始が指示されたか否かを判断する(ステップS10)。呼吸音測定の開始が指示されたと判断すると(ステップS10;Yes)、CPU32は、呼吸音測定部10からの呼吸音データを外部記憶装置35に保存する動作を開始させる(ステップS41)。このとき、被検者のID及び時刻に対応付けて呼吸音データは保存される。呼吸音測定の開始が指示されていないと判断すると(ステップS10;No)、ステップS10に戻り呼吸音測定の開始が指示されるまで待機する。
【0017】
次に、CPU32は、呼吸音測定の終了が指示されたか否かを判断する(ステップS12)。呼吸音測定の終了が指示されたと判断すると(ステップS12;Yes)、CPU32は、呼吸音測定部10からの呼吸音デーを外部記憶装置35に保存する動作を終了させる(ステップS13)。呼吸音測定の終了が指示されていないと判断すると(ステップS12;No)、ステップS12に戻り呼吸音測定の終了が指示されるまで待機する。
【0018】
[解析パラメータ]
まず、本発明の解析に用いるパラメータについて説明する。
【0019】
解析信号z(t)は複素数の信号であり、その実部は原信号と等しく、虚部は原信号のヒルベルト変換である。原信号をx(t)とし、x(t)の解析信号をz(t)とすると、
【0020】
【数1】

【0021】
と表される。但し、
【0022】
【数2】

【0023】
【数3】

【0024】
であり、A(t)、ψ(t)は、それぞれx(t)の包絡線及び瞬時位相を示す。
【0025】
瞬時周波数f(t)は、
【0026】
【数4】

【0027】
で表される。
【0028】
瞬時帯域幅b(t)は、
【0029】
【数5】

【0030】
で表される。
【0031】
[ラ音検出の考え方]
ラ音の注目すべき特徴は時間的局在性である。本実施形態では、エネルギー集中指数ECI(Energy Concentration Index)という指標を用いてラ音を検出する。
【0032】
エネルギー集中指数ECIは、包絡線の最大値(極大値)の近傍における原波形のエネルギーE1と、包絡線の最大値を中心としてラ音1つ分程度の時間幅における原波形のエネルギーE2との比(E1/E2)で表される。エネルギー集中指数ECIは、波形のエネルギーがどのくらい最大値付近に集中しているかの指標となり、ラ音の指標とできる。
【0033】
本アルゴリズムでは、包絡線の最大値付近における瞬時周波数の平均値の逆数をTとして(Tは、そのパルス波形の振動の周期の推定値となる)、E1の値はTの2倍の時間幅に含まれるエネルギー(データの2乗和)、E2の値はTの6倍の時間幅に含まれるエネルギーとした。
【0034】
瞬時周波数は、ラ音部分ではラ音の(主要)周波数を安定的に示す。信号のない部分では小刻みにばらつく。瞬時帯域幅は、ラ音のピークにおいてほぼゼロになる。ラ音部分で下に凸のお椀型の変化を示す。信号のない部分では大きくばらつく。
【0035】
本実施形態では、ラ音の候補点において、エネルギー集中指数ECIがある閾値以上(波形の時間的局在性が大きい)、瞬時周波数の変動が少ない、瞬時帯域幅が小さい、及び候補点の近傍における瞬時帯域幅の積分値(データの総和)(以下、BSと称する)がある閾値以上(瞬時帯域幅の変化がお椀型であるか否かの判断)の場合、ラ音であると判定する。
【0036】
[ラ音検出における呼吸音解析処理]
図3は、第1の実施形態に係る呼吸音解析処理のフロー図である。
【0037】
本呼吸音解析処理は、取得した呼吸音データを解析し、ラ音の検出を行う処理である。ROM33内の呼吸音解析プログラム(本発明のプログラム)に基づいて、CPU32により実行される。尚、予め入力部36により被検者を特定するID等は入力されているものとする。
【0038】
まず、CPU32は、入力部36から呼吸音解析の指示が入力されたか否かを判断する(ステップS20)。呼吸音解析の指示が入力されたと判断すると(ステップS20;Yes)、CPU32は、被検者IDに対応する呼吸音データを外部記憶装置35から読み出し、RAM34にロードする(ステップS21)。呼吸音解析の指示が入力されていないと判断すると(ステップS20;No)、ステップS20に戻り呼吸音解析の指示が入力されるまで待機する。
【0039】
次に、CPU32は、RAM34にロードされた呼吸音データ(サンプリング周波数Fs=11.025kHz)を131072点毎に区分する(ステップS22)。
【0040】
次に、CPU32は、各区間で心音等の雑音を除去するために、カットオフ周波数50Hzのハイパスフィルタをかける(ステップS23)。
【0041】
次に、CPU32は、データをフーリエ変換し、直流成分を1倍、正の周波数成分を2倍、負の周波数成分をゼロとして、逆フーリエ変換を行うことにより解析信号z(t)を得る(ステップS24)。フーリエ変換を用いる方法のほかにヒルベルト変換器(ディジタルフィルタ)を用いる方法もある。
【0042】
次に、CPU32は、包絡線A(t)、瞬時周波数f(t)、瞬時帯域幅b(t)を前掲の式により算出する。微分は差分に置き換えた。それぞれの結果を図4に示す。図4(a)は原信号x(t)、図4(b)は包絡線a(t)、図4(c)はエネルギー集中指数ECI、図4(d)は瞬時周波数f(t)、図4(e)は瞬時帯域幅b(t)である。さらにf(t)の時間微分(差分)をとって2乗した量v(t)を算出する(ステップS25)。
【0043】
次に、CPU32は、重みを3角形型、幅を55点としてA(t)の移動平均As(t)を算出する。同様にb(t)及びv(t)についても、幅を35点としてそれぞれの移動平均bs(t)及びvs(t)を算出する(ステップS26)。
【0044】
次に、CPU32は、移動平均As(t)が極大となる位置を探してラ音の位置の候補とする(ステップS27)。
【0045】
次に、CPU32は、それぞれの候補位置において、エネルギー集中指数ECIを算出する。候補位置を中心として、32点のデータ幅におけるf(t)の平均値FAを求める。ECI=E1/E2とし、E1を求める際のデータ数N1=2×Fs/FA、E2を求める際のデータ数N2=6Fs/FAとする。原信号x(t)において、候補位置を中心としてN1点の幅におけるx(t)の2乗和をE1とし、候補位置を中心としてN2点の幅におけるx(t)の2乗和をE2としてエネルギー集中指数ECIを算出する(ステップS28)。
【0046】
次に、CPU32は、候補位置において、瞬時帯域幅の積分値BSを計算する。瞬時帯域幅b(t)において、候補位置を中心として、442点の幅におけるb(t)の総和を求め、BSとする(ステップS29)。
【0047】
ラ音の候補位置において、エネルギー集中指数ECIが0.65以上、vs(t)が3.04×1010以下、瞬時帯域幅の移動平均bs(t)が154以下、候補点の近傍における瞬時帯域幅の積分値BSが12以上、である場合、その位置の波形をラ音と決定する(ステップS30)。図4においては、図4(a)の原信号x(t)上に3角形の印で示した。
【0048】
以上のように、本実施形態によれば、瞬時周波数、瞬時帯域幅を用いることにより、ラ音の検出を精度良く行うことができる。このため、ラ音のない部分(ノイズ部分)を判別することができ、入力信号の振幅変化にロバストなアルゴリズムが実現できる。また、瞬時周波数を用いて、ラ音の(主要)周波数などを容易に求めることが可能である。
【0049】
本実施形態では、瞬時帯域幅の変化がお椀型であるか否かの判断も組み込まれており、精度を高める上で好ましい。
【0050】
(第2の実施形態)
装置構成及び呼吸音測定処理については、第1の実施形態と同様であり、説明を省略する。
【0051】
[ラ音検出の考え方]
第1の実施形態ではエネルギー集中指数ECIを用いてラ音の検出を行ったが、第2の実施形態においては、後述するクレストインデックス又は包絡線のピーク値まわりの2次モーメントを用いてラ音の検出を行う。
【0052】
クレストインデックスCI(波高指数)は、波高率と同じであり、包絡線の最大値を現波形の実効値で除した値と定義する。波形がパルスに近いほど大きな値を示す。
【0053】
2次モーメントは、確率における分散と同じであり、分布の広がりの程度を示す。ラ音の包絡線では、相対的に低い値を示す。
【0054】
本実施形態では、ラ音の候補点において、クレストインデックスCIがある閾値以上又は2次モーメントの値がある閾値以下(波形の時間的局在性が大きい)、瞬時周波数の変動が少ない、及び瞬時帯域幅が小さい場合、ラ音であると判定する。
【0055】
[ラ音検出における呼吸音解析処理]
図5は、第2の実施形態に係る呼吸音解析処理のフロー図である。一例としてクレストインデックスCIを用いてラ音の検出する場合について説明するが、2次モーメントを用いる場合も基本的に同様である。
【0056】
本呼吸音解析処理は、取得した呼吸音データを解析し、ラ音の検出を行う処理である。ROM33内の呼吸音解析プログラムに基づいて、CPU32により実行される。尚、予め入力部36により被検者を特定するID等は入力されているものとする。
【0057】
まず、CPU32は、入力部36から呼吸音解析の指示が入力されたか否かを判断する(ステップS40)。呼吸音解析の指示が入力されたと判断すると(ステップS40;Yes)、CPU32は、被検者IDに対応する呼吸音データを外部記憶装置35から読み出し、RAM34にロードする(ステップS41)。呼吸音解析の指示が入力されていないと判断すると(ステップS40;No)、ステップS40に戻り呼吸音解析の指示が入力されるまで待機する。
【0058】
次に、CPU32は、RAM34にロードされた呼吸音データ(サンプリング周波数Fs=11.025kHz)を131072点毎に区分する(ステップS42)。
【0059】
次に、CPU32は、各区間で心音等の雑音を除去するために、カットオフ周波数50Hzのハイパスフィルタをかける(ステップS43)。
【0060】
次に、CPU32は、データをフーリエ変換し、直流成分を1倍、正の周波数成分を2倍、負の周波数成分をゼロとして、逆フーリエ変換を行うことにより解析信号z(t)を得る(ステップS44)。フーリエ変換を用いる方法のほかにヒルベルト変換器(ディジタルフィルタ)を用いる方法もある。
【0061】
次に、CPU32は、包絡線A(t)、瞬時周波数f(t)、瞬時帯域幅b(t)を前掲の式により算出する。微分は差分に置き換えた。さらにf(t)の時間微分(差分)をとって2乗した量v(t)を算出する(ステップS45)。
【0062】
次に、CPU32は、重みを3角形型、幅を55点としてA(t)の移動平均As(t)を算出する。同様にb(t)及びv(t)についても、幅を35点としてそれぞれの移動平均bs(t)及びvs(t)を算出する(ステップS46)。
【0063】
次に、CPU32は、移動平均As(t)が極大となる位置を探してラ音の位置の候補とする。但し、移動平均As(t)の平均値を求めて、その3.2倍以下のピークは候補から除く(ステップS47)。
【0064】
次に、CPU32は、それぞれのピークのクレストインデックスCIを算出する。クレストインデックスCIを算出するデータ長は、波形のピーク位置(ラ音の候補)を中心とする1024点の幅とする(ステップS48)。
【0065】
ラ音の候補位置において、クレストインデックスCIが2.6以上、瞬時帯域幅の移動平均bs(t)が154以下、及びvs(t)が3.04×1010以下である場合、その位置の波形をラ音と決定する(ステップS49)。
【0066】
以上のように、本実施形態によれば、瞬時周波数、瞬時帯域幅を用いることにより、ラ音の検出を精度良く行うことができる。このため、ラ音のない部分(ノイズ部分)を判別することができ、入力信号の振幅変化にロバストなアルゴリズムが実現できる。また、瞬時周波数を用いて、ラ音の(主要)周波数などを容易に求めることが可能である。
【0067】
クレストインデックスもしくは包絡線の2次モーメントを計算する範囲(データ長)は、固定長でも良いが、波形のピークおよびその近傍における瞬時周波数をもとに、データ長をそれぞれのピークごとに決定することによって(そのピークの瞬時周波数の逆数の定数倍など)、検出能力を改善することができる。また瞬時周波数を用いてラ音を選択的に検出する、ノイズによる誤検出を減らすなどの処理を行うことができる。
【0068】
本発明における瞬時周波数の波形の粗さに関わる指標とは、例えば、信号のある区間の偏差等のものをいう。偏差が少ないほど瞬時周波数の変動が少なくラ音と判定されやすくなる。
【0069】
本発明における瞬時帯域幅の波形概形とは、ローフィルターを用いた信号の低周波成分でもよいし、区間の移動平均でもよいし、包絡線的なものでもよい。
【0070】
瞬時周波数が振幅の多い波形の場合には、包絡線にするラ音でない部分にも下に凸の領域が現れる場合もあるが、この場合には、瞬時周波数の波形は粗くなるので、ラ音としては検出されない。本発明の検出方法は単独で用いても良いし、他の手法と組み合わせて用いることも好ましい。
【0071】
ファインクラックル(高音性断続性ラ音)の検出については、適当なカットオフ周波数のハイパスフィルタを通して波形整形を行った後、本呼吸音解析処理を行うことによって検出精度を向上させることができる。
【0072】
本発明に用いる呼吸音信号は、例えば電子聴診器により採取することができる。採取装置は電子聴診器のように、耳で聞く機能が必ずしも必要ではなく、体表面から体内の音信号を採取できるものであればよい。
【図面の簡単な説明】
【0073】
【図1】本実施形態に係る呼吸音解析装置の構成図である。
【図2】本実施形態に係る呼吸音測定処理のフロー図である。
【図3】第1の実施形態に係る呼吸音解析処理のフロー図である。
【図4】第1の実施形態に係る各波形を示す概念図である。
【図5】第2の実施形態に係る呼吸音解析処理のフロー図である。
【符号の説明】
【0074】
10 呼吸音測定部
30 呼吸音解析部
32 CPU
33 ROM

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

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図4】
image rotate

【図5】
image rotate


【公開番号】特開2009−106574(P2009−106574A)
【公開日】平成21年5月21日(2009.5.21)
【国際特許分類】
【出願番号】特願2007−282734(P2007−282734)
【出願日】平成19年10月31日(2007.10.31)
【出願人】(303000420)コニカミノルタエムジー株式会社 (2,950)
【出願人】(592243405)
【出願人】(507360003)