JP5593547B2 - 波動源位置演算方法及び波動源位置演算装置 - Google Patents

Description

本発明は、波動源位置演算方法及び波動源位置演算装置に関するものである。

時空間勾配法に基づく３次元音源定位センサシステムに関し、観測点の任意の時間における音圧、音圧微分及び空間的な音圧勾配を観測し、波源拘束偏微分方程式を適用することで、音源の位置を推定する手法が知られている(非特許文献)。

安藤他、"時空間勾配法に基づく３次元音源定位センサシステム"計測自動制御学会論文集 Ｖｏｌ．２９Ｎｏ．５ ５２０／５２８

しかしながら、特定の時間周波数区間に含まれる音源の情報が外乱による影響を受け、当該音源の位置の推定精度が低下する可能性があった。

本発明が解決しようとする課題は、外乱による影響を取り除き、波動源位置の演算精度を高めることができる波動源位置演算方法及び装置を提供する。

本発明は、マスキング関数を用いて、音信号の時間周波数帯域に制限をかけて、波動源の位置を演算することにより、上記課題を解決する。

本発明によれば、マスキング関数により、外乱により影響する帯域を特定し、時間周波数帯域から、当該特定された帯域を取り除き、波動源位置を演算するため、Ｓ／Ｎ比が向上し、波動源位置の演算精度を向上させることができる。

発明の実施形態に係る波動源位置推定方法を用いた波動源位置演算装置のブロック図である。 図１の音源位置演算部の演算結果であって、検出対象の方位角に対する頻度を示すグラフである。 図１の波動源位置演算装置の制御手順を示すフローチャートである。 図１の波動源位置演算装置で用いる式（２３）の閾値（γth）の特性を示すグラフである。 発明の他の実施形態に係る波動源位置推定方法を用いた波動源位置演算装置のブロック図である。 図５の波動源位置演算装置の制御手順を示すフローチャートである。 実施例及び比較例における、方位演算の実験結果を示す図である。

以下、本発明の実施形態を図面に基づいて説明する。

《第１実施形態》
図１は、本例の波動源位置推定方法を用いた波動源位置演算装置のブロック図を示す。なお、本例では、観測点となる移動体として、車両を例に説明するが、本例は車両以外の装置に適用してもよく、また必ずしも移動体である必要はない。

本例の波動源位置演算装置は、演算装置１、音センサ２及び速度センサ３を有する。演算装置１は、ＣＰＵ、ＭＰＵ等の集積回路、ＨＤＤ、メインメモリ及びキャッシュメモリ等の記録媒体等により構成される。音センサ２は、音源である他車両から発せられる音を検出するセンサである。音センサ２には、例えば、少なくとも一つ以上の一般的なマイクロホンが設けられる。速度センサ３は、本例の波動源位置演算装置を備える車両の速度を検出するセンサである。音センサ２及び速度センサ３と、演算装置１との間は、例えば車内信号を伝達するＣＡＮ（Ｃｏｎｔｒｏｌｌｅｒ Ａｒｅａ Ｎｅｔｗｏｒｋ）通信網により接続される。

次に、演算装置１の詳細な構成と本例の波動源位置演算装置の制御内容を説明する。なお、本例は、音センサ２により検出された信号に基づき、例えば非特許文献１の方位推定アルゴリズムを応用して、音源である他車両の位置を演算する。

まず音センサ２は、音源となる他車両から放射される音波を検出し、波動信号として変換部１００に送信する。なお、音センサ２は、自車両の走行音や車両以外の物体を起因とする外音も含めて検出するため、波動信号は、検出対象となる音源以外によるノイズを含んでおり、当該ノイズが、本例の演算において、外乱となる。

変換部１００は、受信した音信号を直行変換等の手法を用い、時間周波数帯域等に変換する。具体的には、受信信号の音信号波形にハン窓（Ｈａｎｎ Ｗｉｎｄｏｗ）を適用し、短時間毎に一定間隔の区間にフレーム分割し、高速離散フーリエ変換（ＦＦＴ）を用いて、フレーム毎に周波数領域の信号Ｓ（ｔ、ω）に変換する。なお、ｔは時間を、ωは周波数を示す。これにより、変換部１００は、波動信号を時間周波数の離散信号に変換する。

マスキング部１１０は、離散信号である信号Ｓ（ｔ、ω）から、外乱を含み演算に適さない領域を除去する。具体的には、後述するマスキング関数設定部により設定されるマスキング関数ｍ（ω）を、信号Ｓ（ｔ、ω）に掛け合わせて、信号Ｓ（ｔ、ω）の時間周波数の帯域に制限をかける。

位置推定演算部２００は、マスキング部１１０からの信号から、瞬時の音圧、音圧時間微分及び音圧の空間勾配を演算し、音源の位置を表す位置推定値ＥＳ（ｔ、ω）を演算する。

ここで、非特許文献１の方位推定アルゴリズムを用いて、変換部１００及び位置推定値演算部２００による、位置推定値ＥＳ（ｔ、ω）の演算方法について、詳述する。なお、演算対象は、マスキング関数ｍ（ω）により制限がかかっていない時間周波数帯域とする。

非特許文献１に示す時空間勾配法のように、音源から観測点に到来する波面にかかる振幅成分、振幅成分の微分及び空間微分成分による、観測点における波面を物理的に表現した式を用いて、観測点における波面から、波動源の方位及び位置を推定する。観測点から音源の方位を推定するには、勾配を元にどの方向から音が到来したのか調べればよい。

当該時空間勾配法は、点音源から放射された音波を任意の点で観測したとき、その点の任意の時間の音圧をｆ（ｔ）、音圧時間微分をｆ_ｔ（ｔ）とすると、観測点における音圧の空間的な勾配ｆ_ｘ（ｔ）、ｆ_ｙ（ｔ）、ｆ_ｚ（ｔ）は、

と表現される。ｎ_ｘ、ｎ_ｙ、ｎ_ｚは、観測点から音源方向へのベクトルｎの要素であり、Ｒは観測点から音源までの距離、Ｃは音速を示す。式１乃至式３は音源が点音源であり、かつ、単一のものであれば、音圧と音圧の空間的な勾配は従属関係にある。以下では、式１乃至式３を波源拘束偏微分方程式と呼ぶ。すなわち、観測点の任意の時間における音圧、音圧時間微分および空間的な音圧勾配を計測し、波源拘束偏微分方程式を適用することで、音源の位置及び方位を演算することができる。

また、変換部１００により、入力された時間領域の音信号を適当な時間間隔の窓関数で切り出してフレーム分割し、短時間高速フーリエ変換によって周波数領域に変換した信号を、波源拘束偏微分方程式を用いて、音源の方位及び位置を推定することもできる。

信号Ｓ（ｔ、ω）から、任意の周波数における音圧をＦ（ω）、音圧時間微分をＦ_ｔ（ω）とすると、観測点における音圧の空間的な勾配音圧の空間的な勾配をＦ_ｘ（ω）、Ｆ_ｙ（ω）、Ｆ_ｚ（ω）は、

として表される。

そして、これらを方向毎独立に最小二乗法を適用することで、ベクトルｎを求めることができるが、接近車音源の水平面の方位を求めたければ、ベクトルｎの要素であるｎ_ｘ、ｎ_ｙ、ｎ_ｚのいずれか一つが求まればよい。例としてｎ_ｘの解法に着目する。（式４）を変形して、

の二乗誤差が最小となるようなｎ_ｘおよびＲを求めることを考える。その結果は、

としたとき、

で与えられる。これにより、ベクトルｎの要素であるｎ_ｘと観測点から音源までの距離Ｒが演算されるため、音源の位置が推定される。そして、上記により演算される音源の位置を示す演算値を位置推定値ＥＳ（ｔ、ω）とする。

なお、（式８）乃至（式１３）は短時間高速フーリエ変換によって得られた所定の周波数帯域について総和を計算しているため、任意の時間区間に対し唯一の位置推定結果を計算しているが、（式８）乃至（式１３）を総和ではなく周波数毎に計算することで、（式１４）乃至（式１５）の結果は、周波数毎に方位推定結果を得ることも可能である。

次に、マスキング部１１０、マスキング関数設定部１２０及び初期設定部１３０について、説明する。初期値設定部１３０は、速度センサ３からの信号に基づき、マスキング関数ｍ（ω）の初期値を設定する。外乱の一つの要因として、自車両が走行することにより発生するノイズが挙げられ、当該ノイズは、自車両との路面との摩擦音や自車両から発生する音等に影響する。そのため自車両の走行速度に応じて、音センサ２が検出するノイズを検出し、帯域フィルタを予め設定する。当該帯域フィルタは、通過帯域ω_ｃにおいてｍ（ω）［０］＝１とし、非通過帯域においてｍ（ω）［０］＝０となるようにｍ(ω)［０］を決定し、

により表される。また、ω_ｃは、速度センサ３の検出速度に応じて設定される。

マスキング関数設定部１２０は、初期値設定部１３０により設定された、（式１６）により表される帯域フィルタの関数を、マスキング関数ｍ（ω）として設定し、マスキング部１１０に送信する。

マスキング部１１０は、マスキング関数設定部１２０により設定されたマスキング関数を、変換部１００からの信号Ｓ（ｔ、ω）に掛け合わせ、当該信号の時間周波数の帯域に制限をかける。これにより、車両の速度に応じて発生するノイズを含む帯域を、信号Ｓ（ｔ、ω）の周波数帯域から除くことができる。

マスキング関数の初期値によりフィルタリングされた信号Ｓ（ｔ、ω）は、位置推定演算部２００に入力され、上記の通り、位置推定値ＥＳ（ｔ、ω）が演算される。マスキング関数設定部１２０は、マスキング関数設定部１２０により入力される当該位置推定値ＥＳ（ｔ、ω）に基づいて、マスキング関数ｍ（ω）を更新し、マスキング部１１０に出力する。これにより、本例は、フィードバック演算制御を行い、演算精度を高める。

以下、本例のフィードバック演算制御方法について説明する。

マスキング関数設定部により設定されるマスキング関数をｍ（ω）［ｉ］とし、ｉは０以上の整数であって、繰り返し回数を示す。ただし、ｉ＝０は初期値を示し、初期値設定部１３０で設定された関数を用いる。

位置推定部演算部２１０において、マスキング関数ｍ（ω）［ｉ］を適用して、（式８）から（式１４）を表すと、

となる。

そして、（式１７）乃至（式２２）と、（式１４）及び（式１５）を用いて、ベクトルｎのｘ方向の要素であるｎｘ及び観測点から音源までの距離Ｒを演算し、位置推定値ＥＳ（ｔ、ω）を演算する。当該位置推定値ＥＳ（ｔ、ω）は、マスキング関数設定部１２０に入力され、マスキング関数設定部１２０は、誤差値ＥＲ（ｔ、ω）を、波源拘束偏微分方程式を用いて演算する。

誤差値ＥＲ（ｔ、ω）は、波源拘束偏微分方程式を用いて表し、位置推定値ＥＳ（ｔ、ω）を代入すると

のように表される。

次に、マスキング関数設定部１２０は、誤差値ＥＲ（ｔ、ω）から、マスキング関数ｍ（ω）［ｉ＋１］を設定する。

マスキング関数ｍ（ω）［ｉ＋１］は、

のように、ER(t、ω)の閾値γ_ｔｈに対する大きさによって表される。閾値γ_ｔｈは閾値γ_ｔｈ（ω）とし、帯域毎に車速に応じた値である。例えば、車速が所定の速度より低いときには１ｋＨｚ未満のγ_ｔｈ（ω）を閾値Ｃ_Ｌに、１ｋＨｚ以上のγ_ｔｈ（ω）を閾値Ｃ_Ｈ（＞Ｃ_Ｌ）にする。車速が当該所定の速度より高いときには２ｋＨｚ未満のγ_ｔｈ（ω）を閾値Ｃ_Ｌに、２ｋＨｚ以上のγ_ｔｈ（ω）を閾値Ｃ_Ｈ（＞Ｃ_Ｌ）にする。（式２４）より、誤差値ＥＲ（ｔ、ω）は、時間周波数帯域における評価関数であって、誤差値ＥＲ（ｔ、ω）が大きければ、外乱が含まれる可能性が高く、誤差値ＥＲ（ｔ、ω）が小さけば、外乱が含まれる可能性が低いと評価される。そして、（式２４）に示すとおり、マスキング関数ｍ（ω）［ｉ＋１］は、誤差値ＥＲ（ｔ、ω）が閾値γ_ｔｈより大きい帯域に制限をかけ、誤差値ＥＲ（ｔ、ω）が閾値γ_ｔｈより小さい帯域に制限をかける関数とする。なお、閾値γ_ｔｈは、音波の振幅レベルに応じて決定してもよい。

そして、マスキング関数設定部１２０は、（式２４）に基づきマスキング関数ｍ（ω）［ｉ］を、マスキング関数ｍ（ω）［ｉ＋１］に更新し、マスキング部１１０に出力する。マスキング部１１０は、マスキング関数ｍ（ω）［ｉ＋１］を用いて、変換部１００からの信号の時間周波数帯域に制限をかけて、位置推定値演算部２１０は、位置推定値ＥＳ（ｔ、ω）を演算する。

またマスキング関数設定部１２０は、マスキング関数ｍ（ω）［ｉ＋１］の収束判定を行う。マスキング関数設定部１２０は、前回のマスキング関数ｍ（ω）［ｉ］を保存しており、ｍ（ω）［ｉ＋１］＝ｍ（ω）［ｉ］の場合に、マスキング関数は収束した、と判断する。なお、収束判断は、例えば、予め設定された、ｉの最大値ｍａｘ（ｉ）と繰り返し回数ｉが等しくなった場合に収束はしたと判断してもよい。

そして、マスキング関数設定部１２０により、マスキング関数が収束したと判断されると、本例の制御は、上記のフィードバック演算制御から抜けだし、音源位置演算部３１０による、以下のクラスタリング工程にうつる。

音源位置演算部３１０は、位置推定値ＥＳ（ｔ、ω）から、周波数毎に位置推定値を統計的に演算し、検出対象となる全方位を任意の幅毎にグループ分けをして、それぞれの領域をビンによるヒストグラム化して、音源の位置を演算する。演算結果の例を図２に示す。図２は、検出対象の方位角に対する頻度を示すグラフである。本例では、検出対象となる方位を−９０°から＋９０°の間にとり、６つの領域にグループ分けをする。そして、位置推定値ＥＳ（ｔ、ω）に基づき、一定の時間内において、フレーム毎に方位を演算し、当該６つの領域に割当て、割当てられた頻度をカウントする。図２の頻度とは、カウントされた頻度に相当する。そして、音源位置演算部３１０は、頻度の高い領域内の位置推定値から、音源の位置を演算する。例えば、図２において、高い頻度である３０〜６０°を方位とし、３０〜６０°の方位に属する位置推定値から距離Ｒを演算することにより、音源位置が演算される。

次に、図３を用いて、本例の制御手順をフローチャートで説明する。図３は、本例の波動源位置演算装置の制御手順を示すフローチャートであり、本例の波動源位置演算方法を示す。

まずステップＳ１にて、音センサ２は車外環境の音波を検出し、速度センサ３は自車両の速度を検出する。ステップＳ２にて、変換部１００は、ステップＳ１にて検出された音波の波動信号を離散信号Ｓ（ｔ、ω）に変換する。次にステップＳ３にて、初期値設定部１３０は、ステップＳ１にて検出された車速に基づき、マスキング関数の初期値ｍ（ω）［０］を設定する。

ステップＳ４にて、マスキング部１１０は、設定されたマスキング関数ｍ（ω）［０］を用いて、離散信号の時間周波数帯域Ｓ（ｔ、ω）に制限をかける。ステップＳ５にて、位置推定値演算部２００は、制限のかかった離散信号から、位置推定値ＥＳ（ｔ、ω）［ｉ］を演算する。ここで、ｉは、ステップＳ４からステップ８までの繰り返し回数を示し、ステップＳ３にて初期値を設定した直後の繰り返し回数をｉ＝０とする。

次にステップＳ６にて、マスキング関数設定部１２０は、誤差値ＥＲ（ｔ、ω）［ｉ］を演算する。ステップＳ７にて、マスキング関数設定部１２０は、演算された誤差値ＥＲ（ｔ、ω）に基づき、マスキング関数ｍ（ω）［ｉ＋１］を設定する。そして、ステップＳ８にて、マスキング設定部１２０は、マスキング関数ｍ（ω）［ｉ＋１］とマスキング関数ｍ（ω）［ｉ］を比較し、マスキング関数が収束したか否かを判定する。

マスキング関数が収束していない場合、ステップＳ４に戻り、マスキング部１１０は、ステップＳ７にて設定されたマスキング関数ｍ（ω）［ｉ＋１］を用いて、離散信号の時間周波数帯域Ｓ（ｔ、ω）に制限をかけ、ステップＳ５にて、位置推定値演算部２００は、制限のかかった離散信号から、位置推定値ＥＳ（ｔ、ω）［ｉ＋１］を演算する。

一方、マスキング関数が収束した場合、音源位置演算部３１０は、ステップＳ７にて演算された位置推定値ＥＳ（ｔ、ω）［ｉ］から、音源の位置を演算し、本例の処理を終了する。

上記のとおり、本例の波動源位置演算方法又は波動源位置演算装置は、離散信号を用い、波動信号の振幅成分、波動信号の振幅の微分成分、または、波動信号の空間的な勾配の少なくとも一つのパラメータを含む式により、音源の位置を示す位置推定値を演算する。そして、当該式を用いて表される評価関数に位置推定値を代入して、評価関数の値に応じて、マスキング関数を設定し、マスキング関数により、波動信号の時間周波数帯域に制限をかける。これにより、本例は、外乱を含む周波数帯域を効率的に特定し、音源の位置を演算する過程において時間周波数帯域に制限をかけるため、Ｓ／Ｎ比が向上し、音源の位置を演算する精度を高めることができる。また、音源の位置を演算する過程において、外乱を含む帯域が取り除かれ、音源位置の精度を低下させる帯域における演算が行われないため、演算負荷を軽減させつつ、音源の位置を演算する精度を高めることができる。

また本例は、式（２３）で示す評価関数を用いて、外乱を含む帯域を特定し、フィードバックをかけて、時間周波数の帯域を制限し、音源の位置を演算する。これにより、外乱を含む周波数帯域を効率的に特定し、音源の位置を演算する過程において時間周波数帯域に制限をかけるため、Ｓ／Ｎ比が向上し、音源の位置を演算する精度を高めることができる。

また本例は、当該評価関数の値を、自車両の速度に基づく閾値と比較することにより評価する。これにより、自車両の速度の変化に応じて、外乱を含む周波数帯域を効率的に特定することができるため、音源の位置を演算する精度を高めることができる。

また本例は、マスキング関数の初期値に、自車両の速度を反映させる。これにより、自車両の速度の変化に応じて、外乱を含む周波数帯域を効率的に特定することができるため、音源の位置を演算する精度を高めることができる。

なお本例は、音源の波動信号からの音圧、音圧時間微分及び音圧の空間勾配をパラメータとして用いて、波源拘束偏微分方程式を表すが、音圧、音圧時間微分及び音圧の空間勾配の代わりに、粒子速度や媒質のインピーダンス等を用いてもよい。

例えば、音圧を用いて壁が剛で均一な断面をもつ長い管の内部の音場（すなわち１次元の波動）を音圧で表現する（例えば文献 橘秀樹、“サウンドインテンシティ”オーム社 ｐ２６−３７を参照（文献２））場合、進行波と後退波の音圧を夫々、f、 gの関数であらわすと任意の点の音圧は、１次元の波動方程式の一般解で表され、

となる。この現象を粒子速度で表現すると、

と表現できる。ここで、ρ_０は平均密度、ｃは音速であり、ｆとｇは、空間的及び時間的境界条件に依存する関数であって、それぞれｘの正及び負方向への速度ｃで伝播する波を示す。更に、音圧と粒子速度の間には、媒質の音響インピーダンスを用いて、

の関係がある。すなわち、粒子速度及び媒質の音響インピーダンスが分かれば、音圧を演算することができる（例えば文献 Ｐ．Ａ．Ｎｅｌｓｏｎ＆Ｓ．Ｊ．Ｅｌｌｉｏｔｔ、“Ａｃｔｉｖｅ Ｃｏｎｔｒｏｌ ｏｆ Ｓｏｕｎｄ”、ｐ１６、ｐ１７を参照（文献３））。また、（式２５）、（式２６）及び（式２７）は、３次元に拡張できる（文献２及び３を参照）。そのため、音圧、音圧時間微分及び音圧の空間勾配は、粒子速度や音響インピーダンスを用いて演算することができ、（式１）〜（式３）を導出することができる。なお、粒子速度は、周知の方法（例えば、ｈｔｔｐ：／／ｗｗｗ．ｍｉｃｒｏｆｌｏｗｎ．ｃｏｍ／ または ｈｔｔｐ：／／ｔｏｙｏ．ｃｏ．ｊｐ／pｃｂ／ｏｎｋｙｏｒｕｓｈｉ／ 要約すると、空気分子の振動速度を直接計測する音響粒子速度センサ）を用いることにより検出される。

また本例は、音圧時間微分を、周知の方法により、直接、観測してもよい（例えば文献 Ｎ．Ｏｎｏ、Ａ．Ｓａｉｔｏ、ａｎｄＳ．Ａｎｄｏ、“Ｂｉｏ−ｍｉｍｃｒｙ Ｓｏｕｎｄ Ｓｏｕｒｃｅ Ｌｏｃａｌｉｚａｔｉｏｎ ｗｉｔｈ Ｇｉｍｂａｌ Ｄｉａｐｈａｒａｇｍ” ＩＥＥＪ Ｔｒａｎｓ．ＳＭ、Ｂｏｌ．１２３ Ｎｏ．３、２００３を参照）。

また本例の波動拘束偏微分方程式は、他の方程式により表現され、例えば波動方程式から展開される音圧（波動の振幅）または粒子速度を用いた音場表現方式、または、その組み合わせを用いることで、表現される。

また本例は、波源拘束偏微分方程式を、波動信号の振幅成分、波動信号の振幅の微分成分、及び、波動信号の空間的な勾配のパラメータを用いて表したが、波動信号の振幅成分、波動信号の振幅の微分成分、及び、波動信号の空間的な勾配のパラメータのうち、少なくとも一つのパラメータにより、波源光源微分方程式を表してもよい。

また本例は、上記の波源拘束偏微分方程式の代わりに、音源と観測点が所望の関係の方程式として表現されていればよく、例えば、自動車などのトンネルや高速道路の外壁が存在する領域では、ダクトをモデルとした虚像法による式を用いてもよい（文献２を参照）。

また、本例は、位置推定値演算部２１０又は音源位置演算部３１０において、観測点から音源までの距離を、音センサ２の信号に基づき演算するが、カメラや電波送受信機等を用いて測定してもよい。

本例は、初期値設定部１３０にて、速度センサ３からの検出速度に応じて、マスキング関数ｍ（ω）の初期値を設定するが、当該初期値は予め設定されている一定の値として、マスキング関数ｍ（ω）の初期値を、所定の周波数帯域で信号を通過させ、それ以外の周波数帯域では信号に制限をかけるように設定してもよい。また、本例は、位置推定値演算部で演算されるパラメータに応じて、マスキング関数ｍ（ω）の初期値を設定してもよい。また、本例において、速度センサ３を省略してもよい。

また本例は、マスキング関数設定部１２０及び初期設定部１３０において、マスキング関数の初期値に通過帯域ω_ｃを予め設定し、マスキング関数を更新するが、帯域毎の離散信号Ｓ（ｔ、ω）の大きさ（例えば、振幅の絶対値等）に応じて、

としてｍ(ω)を設定してもよい。なおＳ_ｔｈは、予め設定される閾値である、なお、離散信号Ｓ（ｔ、ω）の大きさは、離散信号の波面の瞬時振幅を用いてもよいし、当該瞬時振幅の微分値を用いてもよい。

なお、式２３において、分母｜Ｆ_ｘ（τ、ω）｜^２の値が０に近づく場合（すなわちｘ軸に対して直行する方向から到来する波源）、誤差ＥＲ(t、 ω)は、理想状態にあっても発散する。そのため、分母｜Ｆ_ｘ（τ、ω）｜^２の値が０に近づく場合、
以下の式を導入する。

ここでεは任意の微小な値である。これにより、（式２３）の発散を防ぐことができる。

またｘ軸に対して直行する方向から到来する波源に対する対応方法として、（式２４）のγthをｎ_ｘに依存する関数としてもよい。図４は、演算部３１０の演算結果における、ｘ軸の角度に対する閾値（γth）の特性を示すグラフである。例えば、ｘ軸に対して直行する方向を０°とした場合、図４に示すように、０°付近に限りγthの値を、０°付近以外の角度領域の値より大きくする。そして、マスキング関数設定部１２０は、（式２３）の計算結果が発散する場合、図４に示す特性を（式２３）の閾値（γth）に設定し、マスキング関数ｍ（ω）を演算する。これにより（式２９）を用いなくても、ｘ軸に対して直行する方向の音源の位置を演算することができる。

また本例は、観測点を音源とし、音センサ２により音波を検出するが、観測対象は音源に限らず波動源であればよく、音センサ２は当該波動源からの波を検出してもよい。そして、上記の本例の波動源位置演算装置及び波動源位置演算方法は、観測点を音源に限らず、波動源に適用し、波動源の位置を演算することができる。

なお、本例の変換部１００が本発明の「変換手段」に相当し、位置推定値演算部２００が「位置推定値演算手段」に、マスキング関数設定部１２０が「マスキング関数設定手段」に、マスクキング部１１０が「マスキング手段」に、音源位置演算部３１０が「波動源位置演算手段」に相当する。

《第２実施形態》
発明の他の実施形態に係る波動源位置演算装置及び波動源位置演算方法を、図５を用いて説明する。本例は上述した第１実施形態に対して、演算装置１における制御内容の一部が異なる。これ以外の構成で上述した第１実施形態と同じ構成は、その記載を適宜、援用する。図５は、本発明の波動源位置演算装置のブロック図を示す。

図５に示すように、速度センサ３及び初期値設定部１３０を設けていない点を除き、基本的な構成は第１実施形態の波動源位置演算装置と同様であるが、制御手順が異なる。

マスキング関数設定部１２０は、位置推定値演算部２００により演算された位置推定値ＥＳ（ｔ、ω）に基づいて、（式１７）〜（式２３）を用いて、誤差値ＥＲ（ｔ、ω）を算出し、（式２４）に基づき、マスキング関数ｍ（ω）を演算する。

マスキング部１１０は、マスキング関数ｍ（ω）を用い、位置推定値演算部２００から入力される位置推定値ＥＳ（ｔ、ω）の周波数の帯域に制限をかける。そして、音源位置演算部３１０は、マスキング関数ｍ（ω）を通過した帯域の位置推定値ＥＳ（ｔ、ω）から、上述した第１実施形態と同様に、統計的に音源の方位角と観測点から音源までの距離を演算し、音源の位置を演算する。

次に、本例の制御手順を、図６を用いて、説明する。図６は本例の波動源位置演算装置の制御手順であり、本例の波動源位置演算方法のフローチャートを示す。

まずステップＳ１１にて、音センサ２は車外環境の音波を検出する。ステップＳ１２にて、変換部１００は、ステップＳ１１にて検出された音波の波動信号を離散信号Ｓ（ｔ、ω）に変換する。ステップＳ１３にて、位置推定値演算部２００は、離散信号から、位置推定値ＥＳ（ｔ、ω）を演算する。次にステップＳ１４にて、マスキング関数設定部１２０は、予め記録媒体等に記録されている理論値と、ステップＳ１３にて演算された位置推定値とを比較し、理論値との誤差値ＥＲ（ｔ、ω）を演算する。

ステップＳ１５にて、マスキング関数設定部１２０は、演算された誤差値ＥＲ（ｔ、ω）に基づき、マスキング関数ｍ（ω）を設定する。ステップＳ１６にてマスキング部１１０は、設定されたマスキング関数ｍ（ω）を用いて、位置推定値の時間周波数帯域ＥＳ（ｔ、ω）に制限をかける。ステップＳ１７にて、制限のかかった位置推定値ＥＳ（ｔ、ω）から、音源の位置を演算し、本例の処理を終了する。

上記のとおり、本例は、波動信号の振幅の微分成分、または、前記波動信号の空間的な勾配の少なくとも一つをパラメータとして有する方程式に含まれるパラメータを用いて、誤差値ＥＲ（ｔ、ω）の評価関数に基づき、外乱を含む帯域を特定し、位置推定値ＥＳ（ｔ、ω）の周波数の帯域に制限をかけて、音源の位置を演算する。

これにより、本例は、外乱を含む周波数帯域を効率的に特定し、音源の位置を演算する過程において時間周波数帯域に制限をかけるため、Ｓ／Ｎ比が向上し、音源の位置を演算する精度を高めることができる。また、音源の位置を演算する過程において、外乱を含む帯域が取り除かれ、音源位置の精度を低下させる帯域における演算が行われないため、演算負荷を軽減させつつ、音源の位置を演算する精度を高めることができる。

また本例は、マスキング関数設定部１２０において、時間周波数の帯域に応じて重みづけをする関数を用いてもよい。例えば、１ｋＨｚ、１．５ｋＨｚ、２ｋＨｚの帯域で、マスキング関数ｍ_１（ω）を、ｍ_１（１ｋＨｚ）＝２、ｍ_１（１．５ｋＨｚ）＝０．８、ｍ_１（２ｋＨｚ）＝１と設定し、誤差値ＥＲ（ｔ、ω）に基づくマスキング関数に掛け合わせて、新たなマスキング関数とする。これにより、例えば、重み付けが無い状態で、演算部３１０による演算結果が５０°の方向の頻度値Ｆ＝１となる場合、本例では頻度値Ｆ＝１×２＋１×０．８＋１×１＝３．８となるように、重み付けを行うことができる。

《第３実施形態》
発明の他の実施形態に係る波動源位置演算装置及び波動源位置演算方法を説明する。本例は上述した第１又は２実施形態に対して、マスキング関数設定部１２０において、評価関数である誤差値の演算方法及びマスキング関数が異なる。

これ以外の構成で上述した第１又は２実施形態と同じ構成は、その記載を適宜、援用する。

本例のマスキング関数設定部１２０において、マスキング関数ｍ（ω）は、複数の短時間フレームに渡った統計量を反映させることもできる。例えば、

のように、複数の時間フレームを示すτを導入することで時間周波数領域に対して最適な帯域を通過させるフィルタを設計することができる。

これにより本例は、複数の時間フレームを統計的に、マスキング関数ｍ（ω）に反映させることができるため、外乱を含む周波数帯域を効率的に特定することができ、Ｓ／Ｎ比を向上させ、音源の位置を演算する精度を高めることができる。

以下、本発明をさらに具体化した実施例を説明する。

本例の波動源位置演算装置を車両に搭載し、試験路を走行し、他車両の音源位置を演算し、演算結果の評価を行った。自車両には、２本のマイクを設けた。測定対象となる他車両は、２台準備し、１台目は後方から自車両を追い抜き、２台目は自車両の前方から接近し、自車両の後方に抜けていくように走行する。そして、（式３０）から（式３５）と（式３７）及び（式２９）を用いて、マスキング関数ｍ（ω）を計算し、当該マスキング関数ｍ（ω）を用いて、音波の離散信号に制限をかけて、再び位置推定値ＥＳ（ｔ、ω）を演算し、演算された方位をプロットした。

また比較例は、車両に搭載されるマイクの位置、測定対象となる他車両及び他車両の走行条件を同様にするが、マスキング関数ｍ（ω）の演算及びマスキング関数ｍ（ω）に基づく信号の時間周波数帯域のフィルタリングは行わず、本例の帯域制限をしない信号から、光源拘束偏微分方程式を用いて、位置推定値を演算し、演算された方位をプロットした。

また実施例及び比較例において、録音条件はサンプリング周波数４８ｋＨｚ、量子化ビット数１６ｂｉｔであり、録音された音響信号には、懸念された流体雑音はあまりみられず、聴感上も、自車エンジン音に加えて他車走行音が取得できていることを確認している。また、分析条件は高い時間分解能を得るために、フレーム長を１２８点（２．６７ｍｓ）、フレームシフトを６４点（１３３ｍｓ）とした。

それぞれのプロットした結果を図７に示す。図７は、時間に対する測定音源の方位を示す。なお、縦軸の方位は、角度を、−１．０〜＋１．０に正規化した値（Ｎｘ）とする。横軸は、時間を示す。図７において、○印は、実施例の値を示し、×印は比較例の値を示す。また△印及び□印は、車両の位置座標をＧＰＳで測定し、自車両の位置を基準とする他車両の方位に変換したものを示す。
≪考察≫

図７に示すとおり、矢印Ａで示される、実施例の検出範囲の方が、矢印Ｂで示される、比較例の検出範囲より広くなっている。これにより、実施例は、ノイズを多く含むと判断された領域を取り除いて音源の方位を演算するため、Ｓ／Ｎ比が向上し、音源を追従する角度範囲が広くなり、他車両の距離が遠い範囲まで正しい演算が得られていることが確認できる。

以上の結果より、マスキング関数を演算して外乱を含む帯域を特定して、特定された帯域を除いた信号に基づいて、音源の位置を演算することにより、検出精度を高めつつ、検出対象範囲を広げることができる。

１…演算装置
２…音センサ
３…速度センサ
１００…変換部
１１０…マスキング部
１２０…マスキング関数設定部
１３０…初期値設定部
２００…位置推定値演算部
３１０…演算部

Claims (7)

1. 観測点において、波動源が放射する波動信号を検出し、当該波動源の位置を演算する波動源位置演算方法であって、
   前記波動信号を時間周波数の離散信号に変換する工程と、
   前記離散信号を用い、前記波動信号の振幅成分、前記波動信号の振幅の微分成分、または、前記波動信号の空間的な勾配の少なくとも一つのパラメータを含む式により、前記時間周波数毎に、前記波動源の位置を表す位置推定値を演算する位置推定値演算工程と、
   前記式を用いて表される評価関数に前記位置推定値（ＥＳ［ｉ］（ただし、ｉは０以上の整数））を代入することで、外乱を含む可能性を誤差値として演算し、前記時間周波数の帯域を制限するマスキング関数を設定するマスキング関数設定工程と、
   前記マスキング関数により前記時間周波数の帯域に制限をかけるマスキング工程と、
   前記位置推定値から波動源の位置を演算する波動源位置演算工程とを含み、
   前記マスキング関数設定工程は、前記誤差値が所定の値未満である場合は通過帯域を示し、かつ、前記誤差値が前記所定の値より大きい場合は非通過帯域を示す関数を、前記マスキング関数として演算する
   ことを特徴とする波動源位置演算方法。
2. 前記位置推定値演算工程は、前記マスキング工程により時間周波数の帯域に制限をかけた前記離散信号を用いて、前記位置推定値を演算することを特徴とする
   請求項１記載の波動源位置演算方法。
3. 前記マスキング工程は、前記位置推定値演算工程により演算された位置推定値の時間周波数の帯域に制限をかけ、
   前記波動源位置演算工程は、前記マスキング工程により、時間周波数の帯域に制限をかけた前記位置推定値から前記波動源の位置を演算することを特徴とする
   請求項１又は２記載の波動源位置演算方法。
4. 前記位置推定値演算工程は、前記マスキング関数（ｍ［ｉ＋１］）に基づき、時間周波数の帯域に制限をかけた前記離散信号を用いて、前記位置推定値（ＥＳ［ｉ＋１］）を演算することを特徴とする請求項１記載の波動源位置演算方法。
5. 前記マスキング関数の初期値は、
   前記波動信号の振幅成分、前記波動信号の振幅の微分成分に応じて設定される値、
   前記観測点の移動速度に応じて設定される値、または、予め設定されている値であることを特徴とする請求項１〜４のいずれか一項に記載の波動源位置演算方法。
6. 前記観測点を含む移動体の速度を検出する工程とをさらに含み、
   前記マスキング関数設定工程は、前記検出速度に応じて前記所定の値を設定することを特徴とする請求項１記載の波動源位置演算方法。
7. 観測点において、波動源が放射する波動信号を検出し、当該波動源の位置を演算する波動源位置演算装置であって、
   前記波動信号を時間周波数の離散信号に変換する変換手段と、
   前記離散信号を用い、前記波動信号の振幅成分、前記波動信号の振幅の微分成分、または、前記波動信号の空間的な勾配の少なくとも一つのパラメータを含む式により、前記時間周波数毎に、前記波動源の位置を表す位置推定値を演算する位置推定値演算手段と、
   前記式を用いて表される評価関数に前記位置推定値を代入することで、外乱を含む可能性を誤差値として演算し、前記時間周波数の帯域を制限するマスキング関数を設定するマスキング関数設定手段と、
   前記マスキング関数により前記時間周波数の帯域に制限をかけるマスキング手段と、
   前記位置推定値から波動源の位置を演算する波動源位置演算手段とを含み、
   前記マスキング関数設定手段は、
   前記誤差値が所定の値未満である場合は通過帯域を示し、かつ、前記誤差値が前記所定の値より大きい場合は非通過帯域を示す関数を、前記マスキング関数として演算する
   ことを特徴とする波動源位置演算装置。

Images