JP6874766B2 - 波源方向推定装置、波源方向推定システム、波源方向推定方法および波源方向推定プログラム - Google Patents

Description

本発明は、波源方向推定装置、波源方向推定システム、波源方向推定方法および波源方向推定プログラムに関する。

上記技術分野において、非特許文献１には、２つ以上のマイクロフォンの受音信号を用いて複数の相互相関関数を求めてから音源（音波の発生源、発生場所）の方向を推定する技術が記載されている。具体的には、全てのマイクロフォンの組み合わせに対して、受音信号間の相互相関関数をそれぞれ求め、各組み合わせに対応する相互相関関数を総和して一つの相互相関関数を計算する。そして、総和により得られた相互相関関数が最大値を与える時間差を音波の到来時間差として算出することで、音波の到来方向を推定する。

２００１年、「マイクロホン・アレイズ」、第8章、シュプリンガー、ベルリン ハイデルベルグ ニューヨーク(CH.5, MICROPHONE ARRAYS, SPRINGER, BERLIN HEIDELBERG NEW YORK, 2001.) 古井、「新音響・音声工学」、近代科学社、２００６年

しかしながら、上記文献に記載の技術では、周囲のノイズレベルが高い環境の場合には、波源の方向を高精度に推定することができなかった。

本発明の目的は、上述の課題を解決する技術を提供することにある。

上記目的を達成するため、本発明に係る波源方向推定装置は、
波源で発生した信号を入力信号として取得する複数の入力信号取得手段と、
前記入力信号取得手段で取得した前記入力信号に基づいて、相互相関関数を計算する相互相関関数計算手段と、
計算した前記相互相関関数に基づいて、包絡関数を抽出する包絡関数抽出手段と、
抽出した前記包絡関数を統合して統合包絡関数を計算する統合包絡関数計算手段と、
計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する推定方向情報生成手段と、
を備えた。

上記目的を達成するため、本発明に係る波源方向推定システムは、
波源で発生した信号を入力信号として取得する複数の入力信号取得手段と、
複数の前記入力信号取得手段で取得した前記入力信号に基づいて、相互相関関数を計算する相互相関関数計算手段と、
計算した前記相互相関関数に基づいて、包絡関数を抽出する包絡関数抽出手段と、
抽出した前記包絡関数を統合して統合包絡関数を計算する統合包絡関数計算手段と、
計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する推定方向情報生成手段と、
前記推定方向情報を表示する表示手段と、
を備えた。

上記目的を達成するため、本発明に係る波源方向推定方法は、
波源で発生した信号を入力信号として取得する入力信号取得ステップと、
前記入力信号取得ステップで取得した前記入力信号に基づいて、相互相関関数を計算する相互相関関数計算ステップと、
計算した前記相互相関関数に基づいて、包絡関数を抽出する包絡関数抽出ステップと、
抽出した前記包絡関数を統合して統合包絡関数を計算する統合包絡関数計算ステップと、
計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する推定方向情報生成ステップと、
を含む。

上記目的を達成するため、本発明に係る波源方向推定プログラムは、
波源で発生した信号を入力信号として取得する入力信号取得ステップと、
前記入力信号取得ステップで取得した前記入力信号に基づいて、相互相関関数を計算する相互相関関数計算ステップと、
計算した前記相互相関関数に基づいて、包絡関数を抽出する包絡関数抽出ステップと、
抽出した前記包絡関数を統合して統合包絡関数を計算する統合包絡関数計算ステップと、
計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する推定方向情報生成ステップと、
をコンピュータに実行させる。

本発明によれば、周囲のノイズレベルが高い環境の場合であっても、波源の方向を高精度に推定することができる。

本発明の第１実施形態に係る波源方向推定装置の構成を示すブロック図である。 本発明の第２実施形態に係る波源方向推定装置の構成を示すブロック図である。 本発明の第２実施形態に係る波源方向推定装置の備える計算法選択テーブルの構成の一例を示す図である。 本発明の第２実施形態に係る波源方向推定装置のハードウェア構成を示すブロック図である。 本発明の第２実施形態に係る波源方向推定装置の処理手順を説明するフローチャートである。 本発明の第３実施形態に係る波源方向推定装置の構成を示すブロック図である。 本発明の第３実施形態に係る波源方向推定装置の処理手順を説明するフローチャートである。 本発明の第４実施形態に係る波源方向推定システムの構成を示すブロック図である。 本発明の第４実施形態に係る波源方向推定システムの表示部で可視化した画像の一例を示す図である。 本発明の第４実施形態に係る波源方向推定システムの表示部で可視化した画像の他の例を示す図である。 本発明の第４実施形態に係る波源方向推定システムにおいて低域成分を除去した相互相関関数および包絡関数の表示例を示す図である。 本発明の第４実施形態に係る波源方向推定システムにおいて一定間隔で離散化された相互相関関数の表示例を示す図である。 本発明の第４実施形態に係る波源方向推定システムにおいて離散化された包絡関数の表示例を示す図である。

以下に、図面を参照して、本発明の実施の形態について例示的に詳しく説明する。ただし、以下の実施の形態に記載されている構成要素はあくまで例示であり、本発明の技術範囲をそれらのみに限定する趣旨のものではない。なお、以下の説明中における「音声信号」とは、音声その他の音響に従って生ずる直接的の電気的変化であって、音声その他の音響を伝送するためのものをいい、音声に限定されない。

また、以下の実施形態に係る波源方向推定装置の推定対象は、空気や水の振動波である音波の発生源に限定されない。地震や地滑りなどの土や固体を媒質とする振動波の発生源にも適用できる。その場合、振動波を電気信号に変換する装置には、マイクロフォンではなく振動センサが用いられる。さらに、気体・液体・固体の振動波だけでなく、電波を用いて方向を推定する場合にも以下の実施形態に係る波源方向推定装置を適用できる。その場合、電波を電気信号に変換する装置にはアンテナが用いられる。以下の実施形態においては、波源は音源と仮定して説明する。

［第１実施形態］
本発明の第１実施形態としての波源方向推定装置１００について、図１を用いて説明する。波源方向推定装置１００は、取得した信号に基づいて、波源位置の方向を推定する装置である。

図１に示すように、波源方向推定装置１００は、入力信号取得部１０１と、相互相関関数計算部１０２と、包絡関数抽出部１０３と、統合包絡関数計算部１０４と、推定方向情報生成部１０５とを含む。

複数の入力信号取得部１０１は、波源で発生した信号を入力信号として取得する。相互相関関数計算部１０２は、複数の入力信号取得部１０１で取得した入力信号に基づいて、相互相関関数を計算する。包絡関数抽出部１０３は、計算した相互相関関数に基づいて、包絡関数を抽出する。統合包絡関数計算部１０４は、抽出した包絡関数を統合して統合包絡関数を計算する。推定方向情報生成部１０５は、計算した統合包絡関数に基づいて、波源の推定方向情報を生成する。

本実施形態によれば、周囲のノイズレベルが高い環境の場合であっても、波源の方向を高精度に推定することができる。

［第２実施形態］
次に本発明の第２実施形態に係る波源方向推定装置について、図２乃至図５を用いて説明する。

＜前提技術＞
上記非特許文献１および非特許文献２に記載の技術では、屋外などの周囲騒音レベルが高い環境において、遠方に存在する音源の方向を高精度に推定することが困難であった。例えば、推定対象の音源（目標音源）が、マイクロフォンから遠く離れた場所に存在する場合、目標音源から放射される音の音量は、マイクロフォンに到達する時点で大幅に小さくなる。このため、目標音源の音が周囲環境雑音に埋もれて、目標音源の方向推定精度が低下することがあった。

＜本実施形態の技術＞
図２は、本実施形態に係る波源方向推定装置２００の構成を説明するためのブロック図である。本実施形態の波源方向推定装置２００は、例えば、デジタルビデオカメラ、スマートフォン、携帯電話、ノートパソコン、パッシブソーナーなどといった装置の一部として機能する。また、不審ドローン検知、悲鳴検知、車両事故検知などの声や音に基づいて異常を検知する異常音検知装置にも搭載される。しかし、本実施形態に係る波源方向推定装置２００の適用例はこれらに限定されるものではなく、受信音から目標音源の方向推定を要求されるあらゆる波源方向推定装置に適用可能である。

波源方向推定装置２００は、相互相関関数計算部２０１と、包絡関数抽出部２０２と、統合包絡関数計算部２０３と、推定方向情報生成部２０４と、相対遅延時間計算部２０６とを備える。

入力端子２０₁〜２０_Mには、目標音源の音と集音装置であるマイクロフォン(以下、マイク)の周囲で生じる様々な雑音が混在した音とがデジタル信号(サンプル値系列)として入力される。入力端子２０₁〜２０_Mに入力された音信号を本実施形態では入力信号と呼ぶ。そして、時刻ｔにおける入力端子２０_a(ただし、ａ＝ｍ１，・・・，ｍM)の入力信号をｘａ（ｔ）と表す。

入力端子に入力される音は、集音装置であるマイクで集音される。入力端子は複数存在するので、目標音源の音を集音する場合には、端子数と同じＭ個のマイクが同時に使用される。本実施形態では、入力端子とマイクとは一対一に対応しており、ａ番目のマイクが集音した音は、ａ番目の入力端子に供給されると仮定する。したがって、ａ番目の入力端子に供給された入力信号のことを「ａ番目のマイクの入力信号」とも呼ぶ。

波源方向推定装置２００への入力音を集音するＭ個のマイクは、様々なレイアウトで配置される。最も基本的な配置は、一直線上に一列に並ぶ配置である。また、平面上に円状や格子状に配置する場合や、不規則な間隔でランダムに配置することもある。その他には、正多面体の頂点や球体上に配置することもある。入力端子に供給される入力信号は、このように配置されたマイクで集音された音であると仮定する。

波源方向推定装置２００は、目標音源の音が複数のマイクに到達する時間差を利用して音源の方向を推定する。このため、マイク間隔も重要な情報となるので、波源方向推定装置２００には、入力信号だけでなく、マイク位置情報も供給される。

入力端子がＭ個存在する場合、相互相関関数計算部２０１には、Ｍ種類の入力信号ｘａ（ｔ）が入力される。ただし、ａは入力端子番号である。相互相関関数計算部２０１は、入力端子から供給された入力信号から、適当な長さの波形を、一定の周期でずらしながら切り出す。こうして切り出した信号区間をフレーム、切り出した波形の長さをフレーム長、フレームをずらす周期をフレーム周期と呼ぶ。そして、相互相関関数計算部２０１は、切り出された信号の相互相関関数を計算する。このような切り出しは、全ての入力信号を対象に実施される。ｎ番目のフレームの時刻をｔ_nとすると、入力端子ｍ１の入力信号ｘｍ１（ｔ）と、入力端子ｍ２の入力信号ｘｍ２（ｔ）との相互相関関数ｙｍ１ｍ２（τ，ｎ）は、次の式で計算される。

ただし、τはラグサンプル、Ｔはフレーム長を表す。

相互相関関数計算部２０１は、入力信号の全ての組み合わせに対して相互相関関数を計算する。したがって、入力端子がＭ個存在する場合、相互相関関数の種類数はＭ（Ｍ−１）／２になる。例えば入力端子が８個の場合、２８種類の相互相関関数を計算する。計算量を削減する目的で、一部の組み合わせだけに対して相互相関関数を求めることも可能であるが、音源方向の推定精度が低下する。

包絡関数抽出部２０２は、相互相関関数計算部２０１から供給された相互相関関数の包絡（包絡関数）を抽出する。包絡の抽出は、全ての相互相関関数に対して独立に実施される。

包絡とは、相互相関関数を波形と捉えた場合に、その波形から高周波数成分が除去されたものである。つまり、関数値の急激な時間変化が除去された、緩やかな形状を有する関数を「相関関数の包絡」と定義する。したがって、定数倍や単純な線形変換では、包絡を得ることはできない。

最も単純な包絡抽出方法としては、相互相関関数の絶対値をラグサンプル軸方向に移動平均する方法が挙げられる。しかし移動平均では、相関関数のピークと非ピークとの差が大幅に小さくなるので、相関関数のピーク位置に基づいて方向推定を行う本実施形態の方法では有効性は低くなる。したがって、ピークと非ピークとの差は小さくならずに、関数の形状変化が緩やかになることが、包絡の抽出では望ましい。

以上の条件を考慮して、本実施形態では、包絡の抽出にヒルベルト変換を用いた例を説明する。まず、相互相関関数ｙｍ１ｍ２（τ，ｎ）のフーリエ変換を計算する。求めたフーリエ変換をＹｍ１ｍ２（ｋ，ｎ）（ただし、ｋは周波数ビン番号、ｋ＝０，１，・・・，Ｋ−１）とすると、以下の写像変換を行う。

そして、写像関数Ｚｍ１ｍ２（ｋ，ｎ）の逆フーリエ変換の複素数の絶対値を取って得られた関数ｚｍ１ｍ２（τ，ｎ）が、相互相関関数ｙｍ１ｍ２（τ，ｎ）の包絡関数となる。

他には、非特許文献２に記載されているケプストラム変換や、線型予測分析を用いた方法でも包絡を抽出できる。これらの方法は、音声信号処理の分野において、音声信号のパワースペクトルからスペクトルの包絡成分を抽出する際に用いられている。ヒルベルト変換と比較すると、計算量は多くなるものの、包絡を高精度に抽出することが可能である。

相対遅延時間計算部２０６は、入力されたマイク位置情報と音源探索対象方向とから、２つのマイクペア間の相対遅延時間を求める。相対遅延時間とは、マイク間隔と音源方向とに基づいて一意に定まる音波の到達時間差のことである。音速をｃと仮定し、ある２つのマイクの間隔をｄ、音源方向、つまり音の到来方向をθとした場合、音波の到達時間差である相対遅延時間τは次の式で計算される。

相対遅延時間計算部２０６では、全ての方向とマイクペアに対して、上記のτを計算する。

音源探索方向が１軸である場合、例えば垂直方向を無視して水平方向だけを対象とする場合は、方向パラメータは一次元である。一方、２軸である場合、つまり水平方向と垂直方向との両方の角度を用いて音源方向を推定する場合、方向パラメータは二次元となる。本実施形態では、方向パラメータが２次元の場合について説明する。また、水平方向の角度を方位角、垂直方向の角度を仰角と呼ぶこととする。

音源方向推定では、探索範囲にある全ての方向に対して、包絡関数を統合したものである統合包絡関数を計算する。そこで、相対遅延時間計算部２０６では、各マイクペアに対して、全ての方向に対応する相対遅延時間を計算する。例えば、水平方向の探索範囲が１０度刻みで０度から９０度まで、つまり０度、１０度、２０度、・・・９０度であり、垂直方向の探索範囲が１０度刻みで０度から６０度まで、つまり０度、１０度、２０度、・・・６０度である場合を考える。この場合、相対遅延時間計算部２０６は、ある特定のマイクペアに対して７０種類の相対遅延時間を計算する。これを全てのマイクペアに対して計算するので、マイク数がＭ個の場合、相対遅延時間は合計で７０・Ｍ（Ｍ−１）／２種類計算される。

相対遅延時間計算部２０６に入力されるマイク位置情報は、相対遅延時間の計算に用いるマイク間隔を得るために必要となる。したがって、マイク位置情報は、全てのマイクペアに対するマイク間隔(相対的な位置関係)の場合もあれば、マイク位置座標(絶対的な位置関係)の場合もある。本実施形態では、マイクが３次元空間上に配置されていると仮定し、マイク位置情報が３次元座標である場合について説明する。

相対遅延時間の計算方法は次の通りである。まず、三次元空間におけるマイク位置座標ｐ_ｍ ^Ｔ、音源方向を表す単位ベクトルである信号伝搬ベクトルｕ(θ、Φ)をそれぞれ次のように定義する。

ただし、θは方位角、Φは仰角を表す。また、Ａ^Ｔは、行列Ａの転置行列を表す。このとき、マイクｍ１とマイクｍ２との組み合わせに対する相対遅延時間τ_m1m2（θ、Φ）は、次のように与えられる。

ただし、ｃは音速である。式（５）で得られる相対遅延時間は実数なので、統合包絡関数計算部２０３には整数化した相対遅延時間τ_m1m2（θ、Φ）が供給される。

統合包絡関数計算部２０３は、相対遅延時間計算部２０６から供給された相対遅延時間を用いて、包絡関数抽出部２０２から供給された包絡関数を統合して、一つの統合包絡関数を算出する。統合包絡関数は、複数の包絡関数を混合したり、重ね合わせることにより生成される。相対遅延時間計算部２０６から供給された、方位角θと仰角Φとに対する相対遅延時間τ_m1m2（θ、Φ）に基づいて、全てのマイクペアに対する包絡関数ｈ_m1m2（τ、n）を重ね合わせて、方位角θと仰角Φとに対する統合包絡関数求める。包絡関数抽出部２０２から供給された包絡関数を ｚ_m1m2（τ、n）とすると、統合包絡関数 Ｈ_n（τ、n）は次式で計算される。

すなわち、和による重ね合わせ、つまり各包絡関数の総和を統合包絡関数とする。また、総和では無く総乗を用いてもよい。その場合、統合包絡関数は次のように計算される。

この他にも、包絡関数を個別に定数倍してから総和を求める重み付け和や、包絡関数の指数や対数を求めてから総和や総乗を計算する方法も同様に用いることが可能である。

そして、推定方向情報生成部２０４は、算出された統合包絡関数を推定方向情報として出力する。各方向(方位角と仰角との組み合わせ)に対して関数値が定まるので、基本的に関数値が高ければ、その方向に音源が存在する可能性が高いと判断できる。

このような推定方向情報は、様々な形で利用される。例えば、関数が複数のピークを有する場合には、各ピークを到来方向とする音源が複数存在すると考えられる。したがって、各音源の方向を同時に推定できるだけでなく、音源数の推定にも用いることが可能である。

また、関数のピークと非ピークとの差分に基づき、音源の存在可能性を判定することも可能である。もし、ピークと非ピークとの差分が大きければ音源の存在可能性が高いと判定できる。同時に、推定方向の信頼性も高いと判断できる。もし、音源数が１つと事前に仮定できる場合には、関数値が最大となる方向を推定方向情報として出力してもよい。この場合、推定方向情報は、統合包絡関数では無く、方向そのものとなる。

図３は、本実施形態に係る波源方向推定装置２００の備える計算法選択テーブルの構成の一例を示す図である。計算法選択テーブル３０１は、計算対象３１１に関連付けて計算法３１２を記憶する。例えば、包絡関数抽出法が計算対象の場合、ヒルベルト変換やケプストラム変換、線型予測分析などの計算法があり、波源方向推定装置２００は、計算法選択テーブル３０１を参照して計算法を選択してもよい。

図４は、本実施形態に係る波源方向推定装置２００のハードウェア構成を示すブロック図である。

ＣＰＵ(Central Processing Unit)４１０は演算制御用のプロセッサであり、プログラムを実行することで図２の波源方向推定装置２００の機能構成部を実現する。ＲＯＭ(Read Only Memory)４２０は、初期データおよびプログラムなどの固定データおよびプログラムを記憶する。また、通信制御部４３０は、ネットワークを介してその他の装置などと通信する。なお、ＣＰＵ４１０は１つに限定されず、複数のＣＰＵであっても、あるいは画像処理用のＧＰＵ(Graphics Processing Unit)を含んでもよい。また、通信制御部４３０は、ＣＰＵ４１０とは独立したＣＰＵを有して、ＲＡＭ(Random Access Memory)４４０の領域に送受信データを書き込みあるいは読み出しするのが望ましい。また、ＲＡＭ４４０とストレージ４５０との間でデータを転送するＤＭＡＣ(Direct Memory Access Controller)を設けるのが望ましい（図示なし）。さらに、入出力インタフェース４６０は、ＣＰＵ４１０とは独立したＣＰＵを有して、ＲＡＭ４４０の領域に入出力データを書き込みあるいは読み出しするのが望ましい。したがって、ＣＰＵ４１０は、ＲＡＭ４４０にデータが受信あるいは転送されたことを認識してデータを処理する。また、ＣＰＵ４１０は、処理結果をＲＡＭ４４０に準備し、後の送信あるいは転送は通信制御部４３０やＤＭＡＣ、あるいは入出力インタフェース４６０に任せる。

ＲＡＭ４４０は、ＣＰＵ４１０が一時記憶のワークエリアとして使用するランダムアクセスメモリである。ＲＡＭ４４０には、本実施形態の実現に必要なデータを記憶する領域が確保されている。入力信号データ４４１は、マイクなどの集音装置で集音した音信号データや入力信号取得装置などに入力され、取得された信号データなどである。

相互相関関数データ４４２は、相互相関関数計算部２０１で計算した相互相関関数に関するデータである。包絡関数データ４４３は、包絡関数抽出部２０２で抽出した包絡関数に関するデータである。統合包絡関数データ４４４は、統合包絡関数計算部２０３で計算した統合包絡関数に関するデータである。推定方向情報データ４４５は、推定方向情報生成部２０４で生成した推定方向情報に関するデータである。計算法４４６は、例えば、包絡関数の抽出法に関するデータであり、図３に示した計算法選択テーブル３０１から展開されたデータである。

入出力データ４４７は、入出力インタフェース４６０を介して入出力されるデータである。送受信データ４４８は、通信制御部４３０を介して送受信されるデータである。また、ＲＡＭ４４０は、各種アプリケーションモジュールを実行するためのアプリケーション実行領域４４９を有する。

ストレージ４５０には、データベースや各種のパラメータ、あるいは本実施形態の実現に必要な以下のデータまたはプログラムが記憶されている。ストレージ４５０は、計算法選択テーブル３０１を格納する。計算法選択テーブル３０１は、図３に示した、計算対象と計算法との関係を管理するテーブルである。

ストレージ４５０は、さらに、入力信号取得モジュール４５１と、相互相関関数計算モジュール４５２と、包絡関数抽出モジュール４５３と、統合包絡関数計算モジュール４５４とを格納する。また、ストレージ４５０は、推定方向情報生成モジュール４５５と、相対遅延時間計算モジュール４５６とを格納する。

入力信号取得モジュール４５１は、入力信号を取得するモジュールである。相互相関関数計算モジュール４５２は、入力信号に基づいて相互相関関数を計算するモジュールである。包絡関数抽出モジュール４５３は、相互相関関数に基づいて包絡関数を抽出するモジュールである。統合包絡関数計算モジュール４５４は、包絡関数を統合して統合包絡関数を計算するモジュールである。推定方向情報生成モジュール４５５は、統合包絡関数に基づいて波源の推定方向情報を生成するモジュールである。相対遅延時間計算モジュール４５６は、相対遅延時間を計算するモジュールである。これらのモジュール４５１〜４５６は、ＣＰＵ４１０によりＲＡＭ４４０のアプリケーション実行領域４４９に読み出され、実行される。制御プログラム４５７は、波源方向推定装置２００の全体を制御するためのプログラムである。

入出力インタフェース４６０は、入出力機器との入出力データをインタフェースする。入出力インタフェース４６０には、表示部４６１、操作部４６２、が接続される。また、入出力インタフェース４６０には、さらに、記憶媒体４６４が接続されてもよい。さらに、音声出力部であるスピーカ４６３や、音声入力部であるマイク、あるいは、ＧＰＳ位置判定部が接続されてもよい。なお、図４に示したＲＡＭ４４０やストレージ４５０には、波源方向推定装置２００が有する汎用の機能や他の実現可能な機能に関するプログラムやデータは図示されていない。

図５は、本実施形態に係る波源方向推定装置２００の処理手順を示すフローチャートである。このフローチャートは、図４のＣＰＵ４１０がＲＡＭ４４０を使用して実行し、図２の波源方向推定装置２００の機能構成部を実現する。

ステップＳ５０１において、波源方向推定装置２００は、入力信号を取得する。ステップＳ５０３において、波源方向推定装置２００は、取得した入力信号に基づいて相互相関関数を計算する。ステップＳ５０５において、波源方向推定装置２００は、計算した相互相関関数に基づいて、包絡関数を抽出する。ステップＳ５０７において、波源方向推定装置２００は、相対遅延時間を計算する。ステップＳ５０９において、波源方向推定装置２００は、計算した相対遅延時間を用いて統合包絡関数を計算する。ステップＳ５１１において、計算した統合包絡関数に基づいて、方向推定情報を生成する。

本実施形態によれば、入力信号に含まれる目標音の到来方向、すなわち目標物体が存在する方向を推定する。環境雑音レベルが高い環境において、目標物が発する音を手掛かりに、目標物が存在する方向を推定する場合に有効である。環境雑音の例としては、繁華街や街頭、街道沿い、人や自動車が多く集まる場所が挙げられる。また、目標物の例としては、人間や動物、自動車、航空機、船舶、水上バイク、ドローン（小型無人機）が挙げられる。

例えば、屋外のテーマパークや展示会場などに接近する不審な自動車・船舶・ドローンなどを検知し、それらの方向を推定することで、不審者や不審物を効率的に取り締ることが可能である。また、音源方向推定を複数箇所で実施することで、目標音源の位置を特定できる。これにより、環境雑音レベルが高い環境でも、銃声や自動車の衝突音の発生箇所などを正確に特定することが可能となる。

［第３実施形態］
次に本発明の第３実施形態に係る波源方向推定装置について、図６および図７を用いて説明する。図６は、本実施形態に係る波源方向推定装置６００の構成を説明するためのブロック図である。本実施形態に係る波源方向推定装置６００は、上記第２実施形態と比べると、変換部６０１、雑音抑圧部６０２、クロススペクトル計算部６０３および包絡関数抽出部６０４を有する点で異なる。その他の構成および動作は、第２実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。

変換部６０１は、入力端子２０₁〜２０_Mから供給された入力信号を変換し、雑音抑圧部６０２へ供給する。変換は、入力信号を複数の周波数成分に分解する目的で実施される。ここでは、代表的なフーリエ変換を用いた場合について説明する。

第２実施形態の相互相関関数計算部２０１と同様に、変換部６０１にはＭ種類の入力信号ｘｍ（ｔ）が入力される。ただし、ｍは入力端子番号である。変換部６０１は、入力端子から供給された入力信号から、適当な長さの波形を一定の周期でずらしながら切り出す。そして、フーリエ変換を用いて切り出された信号を周波数領域信号に変換する。ｎをフレーム番号とし、切り出した入力信号ｘｍ（ｔ，ｎ）（ｔ＝０，１，・・・，Ｔ−１）とすると、ｘｍ（ｔ，ｎ）のフーリエ変換Ｘｍ（ｋ，ｎ）は次のように計算される。

ただし、ｊは虚数単位、ｅｘｐは指数関数を表す。

雑音抑圧部６０２は、変換部６０１から供給された変換信号に基づいて、入力信号に含まれる雑音成分を抑圧し、雑音抑圧信号としてクロススペクトル計算部６０３へ供給する。

雑音成分が低周波数帯域に存在していることが事前に判明している場合には、低周波数の信号成分をゼロにする方法が有効である。入力信号のサンプリングレートをＦｓとすると、Ｘｍ（ｋ，ｎ）の周波数Ｂヘルツ以下の信号成分が抑圧された雑音抑圧信号Ｙｍ（ｋ，ｎ）は次の式で計算される。

もし、特定の周波数帯域に雑音が存在する場合には、その周波数帯域の信号成分を同様の方法で抑圧すればよい。

この他には、入力信号に目的音源の音が含まれていない時間帯に、入力信号の平均パワースペクトルを求め、その平均パワースペクトルに基づいて雑音成分を抑圧する方法がある。このような雑音抑圧方法には、スペクトル減算法、ウィーナーフィルタ法、ＭＭＳＥ
ＳＴＳＡ法（Minimum Mean-Square-Error Short-Time Spectral Amplitude Estimator）などが知られている。マイク設置場所に定常的に存在する環境雑音を抑圧する場合に有効である。

クロススペクトル計算部６０３は、雑音抑圧部６０２から供給された雑音抑圧信号に基づいて、クロススペクトルを計算し、包絡関数抽出部６０４へ伝達する。クロススペクトル計算部６０３は、全てのマイクペアに対して２つの雑音抑圧信号の積を計算する。入力端子ｍ１の雑音抑圧信号Ｙｍ１（ｋ，ｎ）と、入力端子ｍ２の雑音抑圧信号Ｙｍ２（ｋ，ｎ）とのクロススペクトルＰｍ１ｍ２（ｋ，ｎ）は、次の式で計算される。

ただし、Ｘバーは、Ｘの複素共役を表す。クロススペクトル計算部６０３は、入力信号の全ての組み合わせに対してクロススペクトルを計算する。したがって、入力端子がＭ個存在する場合、クロススペクトルの種類数はＭ（Ｍ−１）／２になる。

クロススペクトルは相互相関関数のフーリエ変換であるので、クロススペクトルも広義の相互相関関数と言える。したがって、ここでクロススペクトルを逆フーリエ変換により相互相関関数に変換すれば、包絡関数抽出部６０４の代わりに、第２実施形態に記載した包絡関数抽出部２０２を用いることが可能である。本実施形態では、クロススペクトルの状態で包絡を抽出する例について説明する。

包絡関数抽出部６０４は、クロススペクトル計算部６０３から供給されたクロススペクトルに基づいて、相互相関関数の包絡を抽出し、包絡関数として統合包絡関数計算部２０３へ伝達する。包絡の抽出は、全てのクロススペクトルに対して独立に実施される。第２実施形態と同様に、本実施形態でもヒルベルト変換を用いた例を説明する。はじめに、以下の写像変換を行う。

そして、写像関数Ｚｍ１ｍ２（ｋ，ｎ）の逆フーリエ変換の複素数の絶対値を取って得られた関数ｚｍ１ｍ２（τ，ｎ）が、包絡関数となる。本実施形態では、相互相関関数のフーリエ変換であるクロススペクトルに対してヒルベルト変換を適用している。このため、第２実施形態のようにヒルベルト変換を行う前のフーリエ変換が不要である。他には、第２実施形態と同様にケプストラム変換や、線型予測分析を用いた方法でも包絡を抽出できる。これらの方法は、ヒルベルト変換と比較して計算量が多くなるものの、包絡を高精度に抽出することが可能である。

以上の構成により、波源方向推定装置６００は、入力信号に含まれる雑音成分の抑圧により目標音を強調してから、目標音の到来方向、すなわち目標物体が存在する方向を推定する。

図７は、本実施形態に係る波源方向推定装置６００の処理手順を示すフローチャートである。このフローチャートは、図４のＣＰＵ４１０がＲＡＭ４４０を使用して実行し、図６の波源方向推定装置６００の機能構成部を実現する。なお、図５と同様の処理には同じステップ番号を付して説明を省略する。

ステップＳ７０１において、波源方向推定装置６００の変換部６０１は、入力信号を変換して、変換した入力信号（変換信号）を雑音抑圧部６０２へ供給する。ステップＳ７０３において、波源方向推定装置６００の雑音抑圧部６０２は、供給された変換信号に基づいて、入力信号に含まれる雑音成分を抑圧して、雑音抑圧信号を生成する。ステップＳ７０５において、波源方向推定装置６００のクロススペクトル計算部６０３は、雑音抑圧信号に基づいて、クロススペクトルを計算し、包絡関数抽出部６０４へ供給する。ステップＳ７０７において、波源方向推定装置６００の包絡関数抽出部６０４は、供給されたクロススペクトルに基づいて、相互相関関数の包絡を抽出し、包絡関数として統合包絡関数計算部２０３へ供給する。

本実施形態によれば、特定の周波数帯域に高いレベルのパワーを有する雑音成分が存在する環境において、目標音源の方向を推定する場合に有効である。例えば、自動車走行雑音のパワーは低域に集中していることから、低域を除去することで、高速道路や街道沿いなどの自動車の走行雑音レベルが高い環境で、高精度に音源方向推定を行うことが可能である。

また、目標音のレベルが環境雑音よりも高い帯域が事前に判明している場合にも有効である。例えば、バイクやヘリコプター、ドローンは人間の声や自動車走行音と比べると、高周波数帯域のパワーが大きい。このため、バイクやヘリコプターなどの方向を推定したい場合にも、低域成分を除去し、目標音が支配的な高域成分を強調することで、音源方向推定の推定精度が向上する。

［第４実施形態］
次に本発明の第４実施形態に係る波源方向推定システムについて、図８乃至図１２を用いて説明する。図８は、本実施形態に係る波源方向推定システム８００の構成を説明するための図である。本実施形態に係る波源方向推定システム８００は、上記第２実施形態に係る波源方向推定装置または第３実施形態に係る波源方向推定装置を用いている。したがって、第２実施形態および第３実施形態と同じ構成および動作については同じ符号を付してその詳しい説明を省略する。

本実施形態に係る波源方向推定システム８００は、マイク８０₁〜８０_M、ＡＤ変換部８０１および表示部８０２を有する。なお、本実施形態では、波源方向推定装置２００の代わりに波源方向推定装置６００を用いることが可能である。また、波源が音源であるという仮定で説明するため、マイクを用いた例を説明するが、音源以外の場合には、その波源から放射される波動を受信して電気信号に変換できる各種センサが、マイクの代わりに用いられる。

マイク８０₁〜８０_Mは、推定対象である目標物体から生じる音を含めた装置周辺の音を電気信号に変換し、ＡＤ変換部８０１へ伝達する。音が伝わる媒質が空気媒質である場合、音は空気の振動としてマイクに到達する。マイクは、到達した空気の振動を電気信号に変換する。

ＡＤ変換部８０１は、マイク８０₁〜８０_Mから供給された音の電気信号をデジタル信号に変換し、入力端子２０₁〜２０_Mに伝達する。

表示部８０２は、波源方向推定装置２００（６００）から供給された推定方向情報を画像などの可視化データに変換し、ディスプレイなどの表示装置に表示する。最も基本的な可視化方法は、球体などの３次元の立体に投影する方法である。その際、相関値をヒートマップと呼ばれる色の濃淡などで表現する。立体では無く、円などの２次元平面に投影する方法も有効である。球体だと表示されたときに裏側が見づらいという問題がある。正面や上から投影した平面上に表示すれば、死角が無くなり一覧性が向上する。相関値は、色の濃淡ではなく、等高線で表現してもよい。相関値の時間変化を把握しやすくする上では、横軸を方向、縦軸を時間にした二次元平面にプロットした方法も有効である。その際、仰角と方位角は別々にプロットしたほうがよい。

図９Ａは、本実施形態に係る波源方向推定システム８００の表示部で可視化した画像の一例を示す図であり、波源方向推定装置２００（６００）から供給された推定方向情報から得られたものである。図９Ｂは、本実施形態に係る波源方向推定システム８００の表示部で可視化した画像の他の例を示す図であり、包絡関数ではなく、相互相関関数を統合して得た推定方向情報を画像化したものである。

これらは、本実施形態の特徴である包絡抽出の効果を確認する目的で、同一の入力信号から同一のフレーム時刻に取得した。例の作成には、正面方向から接近するドローン(小型無人機)の飛翔音を用いた。多面体の各頂点に配置された複数のマイクを用いて飛翔音を集音した。

図９Ａは、色が黒いほど相関値が高いことを表し、図９Ｂでは色が白いほど相関値が高くなっている。方位角の範囲は０から１８０度である。また、仰角の範囲は０から９０度で、中心に近いほど仰角が大きくなる。したがって、半球体上にプロットした３次元データを上から平面上に投影した画像となっている。図９Ａと図９Ｂとを比較すると、図９Ａの図は、目標音源の方向が明確であるのに対して、図９Ｂの図は、同等の白さの箇所が複数あるため、目標音源の存在とその方向が不明確であることが分かる。

図１０は、本実施形態に係る波源方向推定システムにおいて低域成分を除去した相互相関関数および包絡関数の表示例を示す図である。本実施形態の特徴である包絡の抽出は、目標音が高い周波数帯域に強いパワーを有する場合に有効である。このことは、第３実施形態のように低域成分の除去により、高域成分を強調した音の場合でも同様である。これを図１０を用いて説明する。

図１０に示す通り、高周波数成分のみが存在する音の相互相関関数は、低域成分の欠落が原因で、相関関数のピーク幅が狭くなる。もし包絡抽出を行わず、ピーク幅が狭い状態で相関関数を重ね合わせると、雑音などの影響によりピーク位置が僅かにずれるだけで、重ね合わせによるピークの強調が困難になる。場合によっては、複数のピークが出現し、目標音源のピークが強調されなくなる。

ピーク位置のずれは、雑音だけでなく音速のゆらぎにも影響を受ける。式（３）に示す通り、重ね合わせの基準となる相対遅延時間は、音速の影響を受けるためである。音速は、空気などの媒質の温度に応じて変化する。このため、音源からマイクまでの空気の温度が一様ではない実環境では、音速に揺らぎが生じ、相対遅延時間のずれがもたらされる。

また、ピーク位置のずれは、推定範囲の解像度、つまり方位角と仰角との角度間隔にも影響を受ける。第２実施形態では１０度刻みに探索する方法を説明したが、これを２０度、３０度と大きくすると、式（５）から明らかなように相対遅延時間の時間間隔が粗くなる。

図１１は、本実施形態に係る波源方向推定システムにおいて一定間隔で離散化された相互相関関数の表示例を示す図である。推定範囲の解像度低下がピーク位置のずれをもたらす例について、図１１を用いて説明する。破線は、離散化前の相互相関関数を表す。このように離散化の間隔が大きいと、相関関数のピーク位置は、本来とは大きく異なる位置に一致することがある。このずれは離散化の間隔が小さくなるにつれて発生しにくくなるが、相関関数の重ね合わせに要する計算量が大幅に増大する。

図１２は、本実施形態に係る波源方向推定システムにおいて離散化された包絡関数の表示例を示す図である。図１２に示すように、包絡関数を離散化するのであれば、図１１と同じ離散化の間隔であってもピーク位置の大幅なずれを回避できる。

本実施形態によれば、推定方向情報を画像などの可視化データとして表示するので、ユーザが波源の方向推定情報を視覚的に把握することができる。

［他の実施形態］
以上、実施形態を参照して本願発明を説明したが、本願発明は上記実施形態に限定されるものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。また、それぞれの実施形態に含まれる別々の特徴を如何様に組み合わせたシステムまたは装置も、本発明の範疇に含まれる。

また、本発明は、複数の機器から構成されるシステムに適用されてもよいし、単体の装置に適用されてもよい。さらに、本発明は、実施形態の機能を実現する情報処理プログラムが、システムあるいは装置に直接あるいは遠隔から供給される場合にも適用可能である。したがって、本発明の機能をコンピュータで実現するために、コンピュータにインストールされるプログラム、あるいはそのプログラムを格納した媒体、そのプログラムをダウンロードさせるＷＷＷ(World Wide Web)サーバも、本発明の範疇に含まれる。特に、少なくとも、上述した実施形態に含まれる処理ステップをコンピュータに実行させるプログラムを格納した非一時的コンピュータ可読媒体（non-transitory computer readable medium）は本発明の範疇に含まれる。

［実施形態の他の表現］
上記の実施形態の一部または全部は、以下の付記のようにも記載されうるが、以下には限られない。
（付記１）
波源で発生した信号を入力信号として取得する複数の入力信号取得手段と、
前記入力信号取得手段で取得した前記入力信号に基づいて、相互相関関数を計算する相互相関関数計算手段と、
計算した前記相互相関関数に基づいて、包絡関数を抽出する包絡関数抽出手段と、
抽出した前記包絡関数を統合して統合包絡関数を計算する統合包絡関数計算手段と、
計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する推定方向情報生成手段と、
を備えた波源方向推定装置。
（付記２）
前記相互相関関数計算手段は、前記入力信号を変換して周波数領域信号を得る変換手段を備え、
前記周波数領域信号に基づいて、相互相関関数を計算する付記１に記載の波源方向推定装置。
（付記３）
前記相互相関関数計算手段は、前記周波数領域信号に基づいて、前記入力信号に含まれる雑音を抑圧して強調信号を得る雑音抑圧手段をさらに備え、
前記強調信号に基づいて、相互相関関数を計算する付記２に記載の波源方向推定装置。
（付記４）
前記雑音抑圧手段は、所定の周波数帯域の雑音を抑圧する付記３に記載の波源方向推定装置。
（付記５）
前記入力信号取得手段の位置情報に基づいて、前記入力信号の相対遅延時間を計算する相対遅延時間計算手段をさらに備え、
前記統合包絡関数計算手段は、前記相対遅延時間を用いて、統合包絡関数を計算する付記１乃至４のいずれか１項に記載の波源方向推定装置。
（付記６）
波源で発生した信号を入力信号として取得する複数の入力信号取得手段と、
複数の前記入力信号取得手段で取得した前記入力信号に基づいて、相互相関関数を計算する相互相関関数計算手段と、
計算した前記相互相関関数に基づいて、包絡関数を抽出する包絡関数抽出手段と、
抽出した前記包絡関数を統合して統合包絡関数を計算する統合包絡関数計算手段と、
計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する推定方向情報生成手段と、
前記推定方向情報を表示する表示手段と、
を備えた波源方向推定システム。
（付記７）
波源で発生した信号を入力信号として取得する入力信号取得ステップと、
前記入力信号取得ステップで取得した前記入力信号に基づいて、相互相関関数を計算する相互相関関数計算ステップと、
計算した前記相互相関関数に基づいて、包絡関数を抽出する包絡関数抽出ステップと、
抽出した前記包絡関数を統合して統合包絡関数を計算する統合包絡関数計算ステップと、
計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する推定方向情報生成ステップと、
を含む波源方向推定方法。
（付記８）
波源で発生した信号を入力信号として取得する入力信号取得ステップと、
前記入力信号取得ステップで取得した前記入力信号に基づいて、相互相関関数を計算する相互相関関数計算ステップと、
計算した前記相互相関関数に基づいて、包絡関数を抽出する包絡関数抽出ステップと、
抽出した前記包絡関数を統合して統合包絡関数を計算する統合包絡関数計算ステップと、
計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する推定方向情報生成ステップと、
をコンピュータに実行させる波源方向推定プログラム。

この出願は、２０１６年６月１５日に出願された日本出願特願２０１６−１１８４５５を基礎とする優先権を主張し、その開示の全てをここに取り込む。

Claims (7)

1. 波源で発生した信号に対応する複数の入力信号に基づいて、相互相関関数を計算する相互相関関数計算手段と、
   計算した前記相互相関関数に基づいて、包絡関数を抽出する包絡関数抽出手段と、
   前記入力信号に対応する複数のマイクロフォン、複数の振動センサ又は複数のアンテナの位置情報に基づいて、前記入力信号の相対遅延時間を計算する相対遅延時間計算手段と、
   計算した前記相対遅延時間を用いて、抽出した前記包絡関数の総和若しくは総乗又は抽出した前記包絡関数を用いて求めた値の総和若しくは総乗を計算することにより、抽出した前記包絡関数を統合して統合包絡関数を計算する統合包絡関数計算手段と、
   計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する推定方向情報生成手段と、
   を備えた波源方向推定装置。
2. 前記相互相関関数計算手段は、前記入力信号を変換して周波数領域信号を得る変換手段を備え、
   前記周波数領域信号に基づいて、前記相互相関関数を計算する請求項１に記載の波源方向推定装置。
3. 前記相互相関関数計算手段は、前記周波数領域信号に基づいて、前記入力信号に含まれる雑音を抑圧して強調信号を得る雑音抑圧手段をさらに備え、
   前記強調信号に基づいて、前記相互相関関数を計算する請求項２に記載の波源方向推定装置。
4. 前記雑音抑圧手段は、所定の周波数帯域の雑音を抑圧する請求項３に記載の波源方向推定装置。
5. 波源で発生した信号に対応する複数の入力信号に基づいて、相互相関関数を計算する相互相関関数計算手段と、
   計算した前記相互相関関数に基づいて、包絡関数を抽出する包絡関数抽出手段と、
   前記入力信号に対応する複数のマイクロフォン、複数の振動センサ又は複数のアンテナの位置情報に基づいて、前記入力信号の相対遅延時間を計算する相対遅延時間計算手段と、
   計算した前記相対遅延時間を用いて、抽出した前記包絡関数の総和若しくは総乗又は抽出した前記包絡関数を用いて求めた値の総和若しくは総乗を計算することにより、抽出した前記包絡関数を統合して統合包絡関数を計算する統合包絡関数計算手段と、
   計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する推定方向情報生成手段と、
   前記推定方向情報を表示する表示手段と、
   を備えた波源方向推定システム。
6. 波源で発生した信号に対応する複数の入力信号に基づいて、相互相関関数を計算し、
   計算した前記相互相関関数に基づいて、包絡関数を抽出し、
   前記入力信号に対応する複数のマイクロフォン、複数の振動センサ又は複数のアンテナの位置情報に基づいて、前記入力信号の相対遅延時間を計算し、
   計算した前記相対遅延時間を用いて、抽出した前記包絡関数の総和若しくは総乗又は抽出した前記包絡関数を用いて求めた値の総和若しくは総乗を計算することにより、抽出した前記包絡関数を統合して統合包絡関数を計算し、
   計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する
   波源方向推定方法。
7. 波源で発生した信号に対応する複数の入力信号に基づいて、相互相関関数を計算する相互相関関数計算ステップと、
   計算した前記相互相関関数に基づいて、包絡関数を抽出する包絡関数抽出ステップと、
   前記入力信号に対応する複数のマイクロフォン、複数の振動センサ又は複数のアンテナの位置情報に基づいて、前記入力信号の相対遅延時間を計算する相対遅延時間計算ステップと、
   計算した前記相対遅延時間を用いて、抽出した前記包絡関数の総和若しくは総乗又は抽出した前記包絡関数を用いて求めた値の総和若しくは総乗を計算することにより、抽出した前記包絡関数を統合して統合包絡関数を計算する統合包絡関数計算ステップと、
   計算した前記統合包絡関数に基づいて、前記波源の推定方向情報を生成する推定方向情報生成ステップと、
   をコンピュータに実行させる波源方向推定プログラム。

Images