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

Abstract

虚像音源の誤推定の発生を低減し、高精度に音源の方向を推定するために、複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する複数の入力部と、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを少なくとも二つ選択する信号選択部と、入力信号のペアを構成する少なくとも二つの入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する相対遅延時間計算部と、入力信号のペアと相対遅延時間とを用いて波動の波源の推定方向情報を周波数ごとに生成する少なくとも一つの周波数別推定方向情報生成部と、周波数別推定方向情報生成部によって生成される周波数ごとの推定方向情報を統合する統合部とを備える波源方向推定装置とする。

Description

本発明は、波源方向推定装置、波源方向推定方法、およびプログラムに関する。特に、本発明は、複数のセンサが取得した信号に基づいて波源方向を推定する波源方向推定装置、波源方向推定方法、およびプログラムに関する。

特許文献１および非特許文献１には、２つのマイクロフォンの受音信号の到達時間差から音源の方向を推定する方法が開示されている。特許文献１および非特許文献１の方法では、音波の到達時間差の確率密度関数を周波数別に求め、それらの重ね合わせにより得られた確率密度関数から到達時間差を算出し、音源方向を推定する。

特許文献２には、所定の観測点に伝搬される音と振動とを採取して、音の音源が振動源から発生した音か否かを探査する探査方法について開示されている。特許文献２の方法では、音源から伝搬する音と振動源から伝搬する表面波の振動とを同時に測定する。そして、音の音圧レベルのデータから得られた音源の方向と振動の振動レベルのデータから得られた振動源の方向とを比較し、音源からの音が音の発生を伴う振動源からの音かどうかを判定する。

国際公開第２０１８／００３１５８号 特開２０１０−２３６９４４号公報

M. Kato, Y. Senda, R. Kondo, "TDOA estimation based on phase-voting cross correlation and circular standard deviation," 25th European Signal Processing Conference (EUSIPCO), EURASIP, August 2017, p.1230-1234

特許文献１および非特許文献１の方法によれば、信号対雑音比（ＳＮＲ：Signal-to-Noise Ratio）が高い周波数帯域では、到達時間差の確率密度関数が鋭いピークを形成するため、高ＳＮＲ帯域が少なくても精度よく到達時間差を推定することができる。しかしながら、特許文献１および非特許文献１の方法においては、周波数別の到達時間差の確率密度関数を重ね合わせる際に、音源が存在しなくても、偶然位相が揃うことで重ね合わせた確率密度関数にピークが生成される。そのため、特許文献１および非特許文献１の方法には、虚像音源を誤推定してしまうという問題点があった。

特許文献２の方法によれば、音源からの音が音の発生を伴う振動源からの音か、あるいは、振動を伴わない音源からの音か否かの判定や、振動源が音を伴わない振動源であるか否かの判定を的確に行うことができる。しかしながら、特許文献２の方法では、異なるマイクロフォンの間で位相が偶然揃うことによって音源とは異なる方向の虚像音源の到達時間差を算出し、虚像音源を誤推定する可能性があるという問題点があった。

本発明の目的は、上述した課題を解決し、虚像音源の誤推定の発生を低減し、高精度に音源の方向を推定することができる波源方向推定装置を提供することにある。

本発明の一態様の波源方向推定装置は、複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する複数の入力部と、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを少なくとも二つ選択する信号選択部と、入力信号のペアを構成する少なくとも二つの入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する相対遅延時間計算部と、入力信号のペアと相対遅延時間とを用いて波動の波源の推定方向情報を周波数ごとに生成する少なくとも一つの周波数別推定方向情報生成部と、周波数別推定方向情報生成部によって生成される周波数ごとの推定方向情報を統合する統合部とを備える。

本発明の一態様の波源方向推定方法においては、情報処理装置が、複数のセンサによって取得された波動から変換される電気信号を入力信号として取得し、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを少なくとも二つ選択し、入力信号のペアを構成する少なくとも二つの入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出し、入力信号のペアと相対遅延時間とを用いて波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成し、周波数ごとの推定方向情報を統合する。

本発明の一態様のプログラムは、複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する処理と、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを少なくとも二つ選択する処理と、入力信号のペアを構成する少なくとも二つの入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する処理と、入力信号のペアと相対遅延時間とを用いて波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成する処理と、周波数ごとの推定方向情報を統合する処理とコンピュータに実行させる。

本発明によれば、虚像音源の誤推定の発生を低減し、高精度に音源の方向を推定することができる波源方向推定装置を提供することが可能になる。

本発明の第１の実施形態に係る波源方向推定装置の構成の一例を示すブロック図である。 本発明の第１の実施形態に係る波源方向推定装置における相対遅延時間計算部の処理の一例について説明するための概念図である。 本発明の第１の実施形態に係る波源方向推定装置における相対遅延時間計算部の処理の別の一例について説明するための概念図である。 本発明の第１の実施形態に係る波源方向推定装置が有する周波数別推定方向情報生成部の構成の一例を示すブロック図である。 本発明の第１の実施形態に係る波源方向推定装置が有する周波数別クロスススペクトル生成部の構成の一例を示すブロック図である。 本発明の第１の実施形態に係る波源方向推定装置に少なくとも一つのセンサを追加した構成の一例を示すブロック図である。 本発明の第１の実施形態に係る波源方向推定装置の動作の概略について説明するためのフローチャートである。 本発明の第１の実施形態に係る波源方向推定装置の周波数別推定方向情報生成部の動作について説明するためのフローチャートである。 本発明の第１の実施形態に係る波源方向推定装置の周波数別推定方向情報生成部の周波数別クロススペクトル生成部の動作について説明するためのフローチャートである。 本発明の第２の実施形態に係る波源方向推定装置の構成の一例を示すブロック図である。 本発明の各実施形態に係る波源方向推定装置を実現するハードウェア構成の一例を示すブロック図である。

以下に、本発明を実施するための形態について図面を用いて説明する。ただし、以下に述べる実施形態には、本発明を実施するために技術的に好ましい限定がされているが、発明の範囲を以下に限定するものではない。なお、以下の実施形態の説明に用いる全図においては、特に理由がない限り、同様箇所には同一符号を付す。また、以下の実施形態において、同様の構成・動作に関しては繰り返しの説明を省略する場合がある。また、図面中の矢印の向きは、一例を示すものであり、ブロック間の信号の向きを限定するものではない。

（第１の実施形態）
まず、本発明の第１の実施形態に係る波源方向推定装置について図面を参照しながら説明する。以下においては、本実施形態の波源方向推定装置が、空気や水の振動波である音波の発生源を推定する例について説明する。そのため、本実施形態の波源方向推定装置は、マイクロフォンによって電気信号に変換された振動波を検証する。なお、本実施形態の波源方向推定装置の推定対象は、音波の発生源に限定されず、振動波や電磁波などの任意の波動の発生源（波源とも呼ぶ）を推定することにも使用できる。

（構成）
図１は、本実施形態の波源方向推定装置１０の構成を表すブロック図である。波源方向推定装置１０は、入力端子１１、信号選択部１２、相対遅延時間計算部１３、周波数別推定方向情報生成部１５、および統合部１７を備える。

波源方向推定装置１０は、ｐ個の入力端子１１を備える（ｐは２以上の整数）。また、波源方向推定装置１０は、Ｒ個の周波数別推定方向情報生成部１５を備える（Ｒは１以上の整数）。図１においては、個々の入力端子１１を区別するために、符号の後にハイフンを挟んで１〜ｐの番号を付す。同様に、図１においては、個々の周波数別推定方向情報生成部１５を区別するために、符号の後にハイフンを挟んで１〜Ｒの番号を付す。

〔入力端子〕
入力端子１１−１〜ｐ（入力部とも呼ぶ）のそれぞれは、図示しないマイクロフォン（以下、マイクとも呼ぶ）に接続される。入力端子１１−１〜ｐのそれぞれには、異なる位置に配置されたマイクによって集音された音波（音信号とも呼ぶ）から変換された電気信号が入力信号として入力される。以下において、時刻ｔにおいてｍ番目の入力端子１１−ｍに入力した入力信号をｘ_m（ｔ）と記載する（ｔ：実数、ｍ：１以上ｐ以下の整数）。

マイクは、目標音源で発生した音と、そのマイクの周囲で発生した様々な雑音とが混在した音波を集音し、集音した音波をデジタル信号（サンプル値系列とも呼ぶ）に変換する集音装置である。マイクは、目標音源からの音を集音するために、入力端子１１−１〜ｐのそれぞれに対応させて異なる位置に配置される。以下においては、ｍ番目のマイクによって集音された音波から変換された入力信号は、ｍ番目の入力端子１１−ｍに供給されるものとする。また、以下において、ｍ番目の入力端子１１−ｍに供給された入力信号のことを「ｍ番目のマイクの入力信号」とも呼ぶ。

〔信号選択部〕
信号選択部１２は、入力端子１１−１〜ｐに供給されるＰ個の入力信号のうちの２個の入力信号を選択する。信号選択部１２は、選択した２個の入力信号を周波数別推定方向情報生成部１５−１〜Ｒに出力し、それらの入力信号の供給元であるマイクの位置情報（以下、マイク位置情報とも呼ぶ）を相対遅延時間計算部１３に出力する。ここで、周波数別推定方向情報生成部１５の数Ｒは、入力信号の組み合わせの数Ｒに相当する。信号選択部１２は、２個の入力信号を選択する際に、全ての組み合わせを選択してもよいし、一部の組み合わせを選択してもよい。全ての組み合わせを選択する場合、Ｒは以下の式１で表される。

波源方向推定装置１０は、目標音源からの音が２つのマイクに到達する時間差を利用して音源の方向を推定する。マイクの間隔（以下、マイク間隔とも呼ぶ）が大きすぎると、目標音源からの音が空気や水などの媒質の影響で同一の音として観測されないために方向推定精度が低下する。また、マイク間隔が小さすぎると、２つのマイクの間の音波の到達時間差が小さくなりすぎるために方向推定精度が低下する。そのため、信号選択部１２は、式２に示すように、マイク間隔ｄが一定の範囲内に収まるマイクの入力信号を選択するのがよい（ｄ_min、ｄ_max：実数）。

例えば、信号選択部１２は、マイク間隔ｄが十分に小さい場合、マイク間隔ｄが最大の２個の入力信号を選択してもよい。また、信号選択部１２は、マイク間隔ｄが十分に小さい場合、マイク間隔ｄが大きい方から順に並べ、上位Ｒ位（ｒ＜Ｃ（ｐ，２））までの入力信号の組み合わせを選択してもよい。このように、信号選択部１２が一部の組み合わせを選択することは、方向推定精度の低下を防ぐこと以外に、計算量を低減することに繋がる。

目標音源からの音の２つのマイクへの到達時間差を求める際には、マイク位置情報も重要である。そのため、信号選択部１２は、入力信号に加えて、マイク位置情報を相対遅延時間計算部１３に出力する。

〔相対遅延時間計算部〕
相対遅延時間計算部１３には、信号選択部１２からマイク位置情報が入力される。相対遅延時間計算部１３は、マイク位置情報と音源探索対象方向とを用いて、信号選択部１２が選択した全てのマイクペア間の相対遅延時間を算出する。相対遅延時間とは、マイク間隔と音源方向とに基づいて一意に定まる音波の到達時間差のことである。例えば、音源探索対象方向は、所定の角度刻みで設定される。すなわち、相対遅延時間は、音源探索対象方向の分だけ算出される。相対遅延時間計算部１３は、算出した音源探索対象方向と相対遅延時間とをセットにして周波数別推定方向情報生成部１５に出力する。

相対遅延時間は、マイクペアの位置関係によって算出方法が異なる。以下において、マイクペアの位置関係を二つ挙げ、それらのマイクペアの位置関係ごとに相対遅延時間の算出方法を示す。

図２は、全てのマイクが同一直線上に配置される例である。図２の例では、マイクが３つの場合について説明する。ここで、音速をｃ、マイク間隔をｄ_r、音源探索対象方向（音源方向とも呼ぶ）をθとする。音源方向θは、音源１００の方向を推定するために設定される少なくとも一つの角である。このとき、音源方向θに対する相対遅延時間τ_r（θ）は、以下の式３を用いて計算できる。

マイク間隔ｄは、信号選択部１２によって選択される入力信号の組み合わせにより異なる。そのため、相対遅延時間τ_r(θ)は、組み合わせ番号ｒごとに異なる。例えば、図２のマイクペアＡＢ間の距離をｄ₁とすると、相対遅延時間τ₁（θ）は、以下の式４を用いて計算できる。

また、図２のマイクペアＡＣ間の距離をｄ₂とすると、相対遅延時間τ₂（θ）は、以下の式５を用いて計算できる。

このように、全てのマイクが同一直線上に位置する場合、ある音源に関する相対遅延時間τ_r(θ)はマイク間隔ｄに比例するが、音源方向θはどのマイクから見ても同一とみなせる。

図３は、二組のマイクペアが互いに直交する直線上に配置される例である。図３の例では、音源方向θはマイクペアによって異なる。図３のマイクペアＡＢ間の相対遅延時間τ₁（θ）は、以下の式６を用いて計算できる。

一方、図３のマイクＣＤ間の相対遅延時間τ₂（θ）は、以下の式７で計算できる。

このように、あるマイクペアを基準として、別のマイクペアの相対遅延時間τ_ｒ（θ）は、基準のマイクペアから見た音源方向θの関数として、以下の式８のように一般化できる。なお、基準とするマイクペアはどれを選んでもかまわない。

相対遅延時間計算部１３は、全ての音源探索対象方向に対して相対遅延時間を計算する。例えば、相対遅延時間計算部１３は、音源方向の探索範囲が１０度刻みで０度から９０度まで、つまり０度、１０度、２０度、・・・、９０度である場合、１０種類の相対遅延時間を計算する。そして、相対遅延時間計算部１３は、音源探索対象方向と相対遅延時間とを周波数別推定方向情報生成部１５に出力する。

〔周波数別推定方向情報生成部〕
周波数別推定方向情報生成部１５−１〜Ｒには、信号選択部１２によって選択された全マイクペアのうち１組のマイクペアの入力信号と、相対遅延時間計算部１３から供給された相対遅延時間とが入力される。周波数別推定方向情報生成部１５−１〜Ｒは、入力されたマイクペアの入力信号と相対遅延時間とを用いて、そのマイクペアの入力信号間の周波数別推定方向情報を生成する。

ここで、図４を参照しながら、周波数別推定方向情報生成部１５の詳細な構成について説明する。図４は、周波数別推定方向情報生成部１５のブロック図である。周波数別推定方向情報生成部１５は、変換部１５１、クロススペクトル計算部１５２、平均計算部１５３、分散計算部１５４、周波数別クロススペクトル生成部１５５、逆変換部１５６、および周波数別推定方向情報計算部１５７を備える。

〔変換部〕
変換部１５１には、信号選択部１２から２つの入力信号（入力信号Ａおよび入力信号Ｂ）が入力される。変換部１５１は、信号選択部１２から供給された２つの入力信号を変換信号（周波数領域信号とも呼ぶ）に変換する。変換部１５１は、入力信号を複数の周波数成分に分解するための変換を行う。例えば、変換部１５１は、フーリエ変換を用いて、入力信号を複数の周波数成分に分解する。変換部１５１は、クロススペクトル計算部１５２に変換信号を出力する。

変換部１５１には２種類の入力信号ｘ_m（ｔ）が入力される。ここで、ｍは、入力端子１１の番号である。変換部１５１は、入力端子１１から供給された入力信号から、適当な長さの波形を一定の周期でずらしながら切り出す。このように切り出された信号区間をフレーム、切り出された波形の長さをフレーム長、フレームをずらす周期をフレーム周期と呼ぶ。そして、変換部１５１は、フーリエ変換を用いて、切り出された信号を周波数領域信号に変換する。ここで、ｎをフレーム番号とし、切り出される入力信号をｘ_m（ｔ，ｎ）とする（ｔ＝０、１、・・・、Ｋ―１）。このとき、入力信号ｘ_ｍ（ｔ,ｎ）のフーリエ変換Ｘ_m（ｋ,ｎ）は、以下の式９を用いて計算できる。

なお、上記の式９において、ｊは虚数単位、ｅｘｐは指数関数を表す。また、ｋは、周波数ビン番号を表し、０以上Ｋ−１以下の整数である。以下、簡略化のため、ｋのことを周波数ビン番号ではなく、単に周波数とも呼ぶ。

〔クロススペクトル計算部〕
クロススペクトル計算部１５２には、変換部１５１から変換信号が入力される。クロススペクトル計算部１５２は、変換部１５１から供給される変換信号を用いてクロススペクトルを計算する。クロススペクトル計算部１５２は、算出したクロススペクトルを平均計算部１５３に出力する。

クロススペクトル計算部１５２は、変換信号Ｘ₂（ｋ，ｎ）の複素共役と変換信号Ｘ₁（ｋ,ｎ）との積を計算してクロススペクトルを計算する。ここで、変換信号のクロススペクトルをＳ₁₂（ｋ,ｎ）とする。このとき、クロススペクトル計算部１５２は、以下の式１０を用いてクロススペクトルを計算する。

なお、ｃｏｎｊ（Ｘ₂（ｋ，ｎ））は、Ｘ₂（ｋ,ｎ）の複素共役を表す。また、式１０ではなく、振幅成分で正規化したクロススペクトルを用いてもよい。クロススペクトル計算部１５２は、振幅成分で正規化する場合、以下の式１１を用いてクロススペクトルを計算する。

〔平均計算部〕
平均計算部１５３には、クロススペクトル計算部１５２からクロススペクトルが入力される。平均計算部１５３は、クロススペクトル計算部１５２から供給されたクロススペクトルの平均（平均クロススペクトルとも呼ぶ）を計算する。平均計算部１５３は、算出した平均クロススペクトルを分散計算部１５４および周波数別クロススペクトル生成部１５５に出力する。

ここで、平均計算部１５３が過去に入力されたクロススペクトルから周波数ビンごとに平均クロススペクトルを計算する例について説明する。なお、平均計算部１５３は、周波数ビン単位ではなく、複数の周波数ビンを束ねたサブバンド単位で平均クロススペクトルを計算してもよい。ここで、第ｎフレームの周波数ビンｋにおけるクロススペクトルをＳ₁₂（ｋ，ｎ）とする、このとき、平均計算部１５３は、以下の式１２を用いて、過去Ｌフレームから求めた平均クロススペクトルＳＳ₁₂（ｋ，ｎ）を計算する。

また、平均計算部１５３は、次のリーク積分を用いて平均クロススペクトルＳＳ₁₂（ｋ，ｎ）を計算してもよい。ただし、αは、０よりも大きく、１よりも小さい実数である。

〔分散計算部〕
分散計算部１５４には、平均計算部１５３から平均クロススペクトルが入力される。分散計算部１５４は、平均計算部１５３から供給された平均クロススペクトルを用いて分散を計算する。分散計算部１５４は、算出した分散を周波数別クロススペクトル生成部１５５に出力する。

ここで、平均クロススペクトルをＳＳ₁₂（ｋ，ｎ）とする。このとき、分散計算部１５４は、クロススペクトルの位相の分散の計算において円周分散を用いた場合、以下の式１４を用いて、分散Ｖ₁₂（ｋ，ｎ）を計算する。

なお、分散計算部１５４は、以下の式１５を用いて分散Ｖ₁₂（ｋ，ｎ）を計算してもよい。

また、分散計算部１５４は、円周標準偏差を用いる場合、以下の式１６を用いて、分散Ｖ₁₂（ｋ，ｎ）を計算する。

〔周波数別クロススペクトル生成部〕
ここで、周波数別クロススペクトル生成部１５５の構成について図面を参照しながら説明する。図５は、周波数別クロススペクトル生成部１５５の構成の一例を示すブロック図である。図５のように、周波数別クロススペクトル生成部１５５は、周波数別基本クロススペクトル計算部５５１、カーネル関数スペクトル生成部５５２、および乗算部５５３を有する。

〔周波数別基本クロススペクトル計算部〕
周波数別基本クロススペクトル計算部５５１には、平均計算部１５３から平均クロススペクトルが入力される。周波数別基本クロススペクトル計算部５５１は、平均計算部１５３から供給された平均クロススペクトルを用いて周波数別基本クロススペクトルを計算する。周波数別基本クロススペクトル計算部５５１は、算出した周波数別基本クロススペクトルを乗算部５５３に出力する。

周波数別基本クロススペクトル計算部５５１は、平均計算部１５３から供給される平均クロススペクトルＳＳ₁₂（ｋ，ｎ）を用いて、平均クロススペクトルＳＳ₁₂（ｋ，ｎ）の各周波数ｋに対応するクロススペクトル（周波数別基本クロススペクトルとも呼ぶ）を計算する。周波数別基本クロススペクトル計算部５５１は、算出する周波数別基本クロススペクトルを乗算部５５３に出力する。周波数別基本クロススペクトルは、周波数成分ごとに相関関数を計算するために算出される。周波数別基本クロススペクトル計算部５５１は、ある周波数ｋに対応する相関関数（周波数別相関関数とも呼ぶ）を後段で求めるための周波数別基本クロススペクトルを計算する。

ここで、周波数別基本クロススペクトル計算部５５１が周波数ｋの周波数別基本クロススペクトルを計算する例について詳細に説明する。周波数別基本クロススペクトル計算部５５１は、周波数ｋの平均クロススペクトルＳＳ₁₂（ｋ，ｎ）を用いて周波数別基本クロススペクトルを計算する際、位相成分と振幅成分とを予め別々に求めたのちに統合する。周波数ｋの周波数別基本クロススペクトルＵ_k（ｗ，ｎ）、その振幅成分を｜Ｕ_k（ｗ，ｎ）｜、位相成分をａｒｇ（Ｕ_k（ｗ，ｎ））とすると、以下の式１７の関係が成り立つ。ただし、式１７において、ｗは、周波数を表し、０以上Ｗ−１以下の整数である。

以下において、周波数別基本クロススペクトル計算部５５１が、周波数ｋの平均クロススペクトルＳＳ₁₂（ｋ，ｎ）を用いて、周波数別基本クロススペクトルの振幅成分｜Ｕ_k（ｗ，ｎ）｜と位相成分ａｒｇ（Ｕ_k（ｗ,ｎ））とを求める方法について説明する。

ｋの整数倍の周波数の振幅成分｜Ｕ_k（ｗ，ｎ）｜には１．０が用いられる。一方、ｋの非整数倍の周波数の位相成分はゼロにする。これを数式で表現すると、振幅成分｜Ｕ_k（ｗ，ｎ）｜は以下の式１８で与えられる。ただし、式１８において、pは1以上Ｐ以下の整数である。

波源方向推定を行うときに重要な情報は位相成分であるため、振幅成分には、式１８のように適当な定数を用いる。なお、ｋの整数倍の周波数の振幅成分｜Ｕ_k（ｗ，ｎ）｜として、１．０の代わりに｜ＳＳ₁₂（ｋ，ｎ）｜を用いてもよい。つまり、振幅成分｜Ｕ_k（ｗ，ｎ）｜を以下の式１９を用いて求めてもよい。

ｋを整数倍した周波数の位相成分ａｒｇ（Ｕ_k（ｗ，ｎ））には、周波数ｋの平均クロススペクトルＳＳ₁₂（ｋ，ｎ）を定数倍したものを用いる。例えば、周波数ｋ、２ｋ、３ｋ、および４ｋの位相成分には、周波数ｋの位相成分ａｒｇ（ＳＳ₁₂（ｋ，ｎ））のそれぞれを同一の倍率で整数倍したものが用いられる。すなわち、周波数ｋ、２ｋ、３ｋ、および４ｋの位相成分のそれぞれには、ａｒｇ（ＳＳ₁₂（ｋ，ｎ））、２ａｒｇ（ＳＳ_1２（ｋ，ｎ））、３ａｒｇ（ＳＳ₁₂（ｋ，ｎ））、４ａｒｇ（ＳＳ₁₂（ｋ，ｎ））が用いられる。一方、ｋの非整数倍の周波数の位相成分はゼロにする。したがって、周波数ｋに対応する周波数別基本クロススペクトルの位相成分ａｒｇ（Ｕ_k（ｗ，ｎ））は、以下の式２０を用いて計算される。なお、pは１以上Ｐ以下の整数である（Ｐ＞１）。

周波数別基本クロススペクトル計算部５５１は、式１８や式１９を用いて算出される振幅成分と、式２０を用いて算出される位相成分とを、式１７を用いて統合し、周波数ｋの周波数別基本クロススペクトルＵ_k（ｗ，ｎ）を得る。

これまで説明してきた方法では、振幅成分と位相成分を別々に求めてから周波数別基本クロススペクトルを算出した。しかし、以下の式２１に示すクロススペクトルのべき乗を用いれば、振幅成分と位相成分を求めることなく、周波数別基本クロススペクトルＵ_k（ｗ，ｎ）を求めることができる。

〔カーネル関数スペクトル生成部〕
カーネル関数スペクトル生成部５５２には、分散計算部１５４から分散が入力される。カーネル関数スペクトル生成部５５２は、分散計算部１５４から供給された分散を用いてカーネル関数スペクトルを算出する。カーネル関数スペクトルとは、カーネル関数をフーリエ変換し、その絶対値を取ったものである。なお、カーネル関数スペクトルには、カーネル関数をフーリエ変換し、その絶対値を取る代わりに二乗してもよい。また、カーネル関数スペクトルには、カーネル関数をフーリエ変換し、その絶対値の二乗としてもよい。カーネル関数スペクトル生成部５５２は、算出したカーネル関数スペクトルを乗算部５５３に出力する。

ここで、カーネル関数スペクトルをＧ（ｗ）とし、カーネル関数をｇ（τ）とする。また、カーネル関数としてガウス関数を用いる。このとき、ガウス関数は、以下の式２２で与えられる。

ただし、式２２において、ｇ_1、ｇ₂、およびｇ₃は正の実数である。ｇ₁はガウス関数の大きさを制御し、ｇ₂はガウス関数のピークの位置を制御し、ｇ₃はガウス関数の広がりを制御する。特に、ガウス関数の広がりを調整するｇ₃は、周波数別相関関数のピークの鋭さに大きな影響を与えるので重要である。すなわち、式２２は、ｇ₃が大きくなればガウス関数の広がりが大きくなることを示す。

また、以下の式２３のロジスティック分布の確率密度関数をカーネル関数として用いてもよい。ただし、式２３において、ｇ₄とｇ₅は正の実数である。

ロジスティック分布の確率密度関数は、ガウス関数と同様の形状をしているが、ガウス関数よりも裾が長い。特に、ロジスティック分布の確率密度関数の広がりを調整するｇ₅は、式２２のガウス関数におけるｇ₃の場合と同様に、周波数別相関関数のピークの鋭さに大きな影響を与えるパラメータである。また、コサイン関数や一様関数をカーネル関数に用いてもよい。

カーネル関数のパラメータのうち、カーネル関数の広がりに影響を与えるｇ₃やｇ₅は、分散計算部１５４から入力される分散を用いて決定される。ここで、これらのパラメータのことを広がり制御パラメータと呼び、ｑ（ｋ，ｎ）と表現する。したがって、カーネル関数がガウス関数の場合には、ｇ₃はｑ（ｋ，ｎ）である。分散が小さければ、周波数別相関関数のピークが鋭く、裾が狭くなるようにパラメータを変化させる。したがって、広がり制御パラメータを小さくする。

広がり制御パラメータは、予め設定した写像関数を用いて分散の値を変換することによって算出できる。例えば、分散がある閾値を上回れば広がり制御パラメータを大きな値（例えば１０）にし、分散がある閾値を下回れば小さな値（例えば０．０１）に設定する。ここで、分散をＶ₁₂（ｋ，ｎ）とし、閾値をｐ_thとする。このとき、第ｎフレームの周波数ビンｋにおける広がり制御パラメータｑ（ｋ，ｎ）は、以下の式２４を用いて計算できる。ただし、式２４において、ｑ₁およびｑ₂は、ｑ₁＞ｑ₂を満たす正の実数である。

また、広がり制御パラメータｑ（ｋ，ｎ）は、以下の式２５のように一次関数を用いて計算してもよい。ただし、式２５において、ｑ₃は０よりも大きな実数であり、ｑ₄は実数である。

ｑ_3、ｑ₄としては、例えば、および式２７で示す値を用いてもよい。

Ｌは平均計算部１５３が平均クロススペクトルを求める際に平均化したフレーム数を表す。平均クロススペクトルの誤差は平均化フレーム数Ｌに反比例するため、式２６および式２７を用いることで、平均クロススペクトルの誤差（信頼性）を考慮して広がり制御パラメータを求めることができる。

また、線形写像関数や高次の多項式関数、非線形関数などで表される分散の関数を分散の計算に用いることも可能である。また、分散をそのまま広がり制御パラメータにしてもよい。

広がり制御パラメータを求める関数は、分散だけでなく、周波数ｋの関数としてもよい。例えば、周波数ｋが大きくなるにつれて小さくなる関数を用いることができる。このような代表例としては、ｋの逆数を用いる例が挙げられる。この場合、式２５の代わりに、以下の式２８の関数を用いて広がり制御パラメータｑ（ｋ，ｎ）を計算できる。

また、式２６の代わりとしては、以下の式２９の関数を用いて広がり制御パラメータｑ（ｋ，ｎ）を計算できる。

〔乗算部〕
乗算部５５３には、周波数別基本クロススペクトル計算部５５１から周波数別基本クロススペクトルが入力され、カーネル関数スペクトル生成部５５２からカーネル関数スペクトルが入力される。乗算部５５３は、周波数別基本クロススペクトル計算部５５１から供給される周波数別基本クロススペクトルと、カーネル関数スペクトル生成部５５２から供給されるカーネル関数スペクトルとの積を計算し、周波数別クロススペクトルを算出する。乗算部５５３は、算出した周波数別クロススペクトルを逆変換部１５６に出力する。

ここで、周波数別基本クロススペクトル計算部５５１から供給される周波数別基本クロススペクトルをＵ_k（ｗ，ｎ）、カーネル関数スペクトル生成部５５２から供給されるカーネル関数スペクトルをＧ（ｗ）とする。このとき、乗算部５５３は、以下の式３０を用いて、周波数別クロススペクトルＵＭ_k（ｗ，ｎ）を計算する。

〔逆変換部〕
逆変換部１５６には、周波数別クロススペクトル生成部１５５の乗算部５５３から周波数別クロススペクトルが入力される。例えば、変換部１５１がフーリエ変換を用いる場合、逆変換部１５６は逆フーリエ変換を用いて逆変換を行う。逆変換部１５６は、周波数別クロススペクトル生成部１５５から供給された周波数別クロススペクトルの逆変換を求める。

ここで、周波数別クロススペクトル生成部１５５から供給される周波数別クロススペクトルをＵＭ_k（ｗ，ｎ）とする。このとき、逆変換部１５６は、以下の式３１を用いて、ＵＭ_k（ｗ，ｎ）を逆変換して周波数別相互相関関数ｕ_k（τ，ｎ）を計算する。

〔周波数別推定方向情報計算部〕
周波数別推定方向情報計算部１５７には、逆変換部１５６から周波数別相互相関関数が入力され、相対遅延時間計算部１３から相対遅延時間が入力される。周波数別推定方向情報計算部１５７は、逆変換部１５６から供給された周波数別相互相関関数と、相対遅延時間計算部１３から供給された相対遅延時間とを用いて、方向と相関値との対応関係を周波数別推定方向情報として求める。周波数別推定方向情報計算部１５７は、求めた周波数別推定方向情報を統合部１７に出力する。

ここで、周波数別相互相関関数をｕ_k（τ，ｎ）とし、相対遅延時間をτ_r（θ）とする。このとき、周波数別推定方向情報計算部１５７は、以下の式３２を用いて、周波数別推定方向情報Ｈ_k,r（θ，ｎ）を算出する。

式３２を用いれば、各方向θに対して相関値が定まるので、相関値が高い方向に音源が存在する可能性が高いと判断できる。

〔統合部〕
統合部１７には、周波数別推定方向情報生成部１５−１〜Ｒから周波数別推定方向情報が入力される。統合部１７は、周波数別推定方向情報生成部１５−１〜Ｒから供給される周波数別推定方向情報を統合して統合推定方向情報を算出する。統合部１７は、個別に求めた複数の周波数別推定方向情報を混合したり、重ね合わせたりすることによって一つの推定方向情報を求める。統合部１７は、算出した統合推定方向情報を出力する。例えば、統合部１７は、上位システム（図示しない）に統合推定方向情報を出力する。

例えば、統合部１７は、まず、入力信号の組み合わせ数（Ｒ個）分の周波数別推定方向情報Ｈ_k,r（θ，ｎ）を統合することによって周波数別統合推定方向情報Ｈ_k（θ，ｎ）を算出する。そして、統合部１７は、算出した周波数別統合推定方向情報を全周波数について統合することによって統合推定方向情報Ｈ（θ，ｎ）を算出する。

例えば、統合部１７は、周波数別推定方向情報Ｈ_k,r（θ，ｎ）の総乗を計算することによって周波数別統合推定方向情報Ｈ_k（θ，ｎ）を算出する。このとき、統合部１７は、以下の式３３を用いて周波数別統合推定方向情報Ｈ_k（θ，ｎ）を算出する。

また、例えば、統合部１７は、周波数別推定方向情報Ｈ_k,r（θ，ｎ）の総和を計算することによって周波数別統合推定方向情報Ｈ_k（θ，ｎ）を算出してもよい。このとき、統合部１７は、以下の式３４を用いて周波数別統合推定方向情報Ｈ_k（θ，ｎ）を算出する。

統合部１７は、統合推定方向情報Ｈ（θ，ｎ）の算出においては、周波数別統合推定方向情報Ｈ_k（θ，ｎ）の周波数ｋについての総和または総乗を計算する。

例えば、統合部１７は、以下の式３５を用いて、周波数別統合推定方向情報Ｈ_k（θ，ｎ）の周波数ｋについての総和を統合推定方向情報Ｈ（θ，ｎ）として計算する。

また、例えば、統合部１７は、以下の式３６を用いて、周波数別統合推定方向情報Ｈ_k（θ，ｎ）の周波数ｋについての総乗を統合推定方向情報Ｈ（θ，ｎ）として計算する。

目標音の存在する周波数や、目的音のパワーが大きい周波数が予め判明している場合、統合部１７は、その周波数に対応する周波数別統合推定方向情報だけを用いて、統合推定方向情報を求めてもよい。また、統合部１７は、重み付けという形で、統合における周波数別統合推定方向情報の影響度を制御してもよい。例えば、目標音の存在する周波数の集合をΩとすると、統合部１７は、以下の式３７を用いて、周波数の選択により統合推定方向情報Ｈ（θ，ｎ）を求めることができる。

また、重み付けを用いる場合には、統合部１７は、以下の式３８を用いて統合推定方向情報Ｈ（θ，ｎ）を算出できる。ただし、式３８において、ａおよびｂは、ａ＞ｂ＞０を満たす実数である。

このように、目的音が存在する周波数の周波数別統合推定方向情報を重点的に用いて統合すると、雑音などの非目的音の影響が小さい相関関数を生成できるため、方向推定精度が向上する。

また、統合部１７は、別の算出方法を用いて統合推定方向情報Ｈ（θ，ｎ）を算出してもよい。例えば、統合部１７は、まず、全周波数について周波数別推定方向情報Ｈ_k,r（θ，ｎ）を統合した入力信号組み合わせ別統合推定方向情報Ｈ_r（θ，ｎ）を算出する。そして、統合部１７は、入力信号の全組み合わせについて入力信号組み合わせ別統合推定方向情報を統合した統合推定方向情報Ｈ（θ，ｎ）を算出してもよい。

以上が、本実施形態の波源方向推定装置１０の構成についての説明である。

なお、図６のように、波源方向推定装置１０にマイクなどのセンサ１１０を少なくとも一つ追加した構成も本実施形態の範囲に含まれる。それぞれのセンサ１１０は、インターネットやイントラネットなどネットワークやケーブルを介して、波源方向推定装置１０のいずれかの入力端子１１に接続される。

例えば、センサ１１０は、音波を検出する場合、マイクロフォンによって実現される。例えば、センサ１１０は、振動波を検出する場合、振動センサによって実現される。例えば、センサ１１０は、電磁波を検出する場合、アンテナによって実現される。なお、センサ１１０は、検知対象の波動を電気信号に変換できるものであれば、その形態に限定を加えない。

（動作）
次に、本実施形態の波源方向推定装置１０の動作について図面を参照しながら説明する。

〔波源方向推定〕
まず、図７のフローチャートを用いて、波源方向推定装置１０の動作の概略について説明する。なお、図７のフローチャートに沿った説明においては、波源方向推定装置１０を動作の主体として説明する。

図７において、まず、波源方向推定装置１０は、複数のマイクから電気信号（入力信号とも呼ぶ）を入力する（ステップＳ１１１）。

次に、波源方向推定装置１０は、複数のマイクに対応する入力信号から２つの入力信号を選択する（ステップＳ１１２）。

次に、波源方向推定装置１０は、選択した２つの入力信号の供給元である２つのマイクの間隔（マイク間隔とも呼ぶ）と、設定された音源探索対象方向とに基づいて相対遅延時間を算出する（ステップＳ１１３）。

次に、波源方向推定装置１０は、選択された２つの入力信号と相対遅延時間とを用いて、周波数ごとに推定方向情報（周波数別推定方向情報とも呼ぶ）を生成する（ステップＳ１１４）。

次に、波源方向推定装置１０は、周波数ごとの推定方向情報を統合して統合推定方向情報を算出する（ステップＳ１１５）。

そして、波源方向推定装置１０は、統合推定方向情報を出力する（ステップＳ１１６）。

以上が、波源方向推定装置１０の動作の概略についての説明である。

〔周波数別推定方向情報生成〕
次に、図８のフローチャートを用いて、波源方向推定装置１０の周波数別推定方向情報生成部１５の動作について説明する。図８のフローチャートの処理は、図７のフローチャートのステップＳ１１４を細分化したものである。なお、図８のフローチャートに沿った説明においては、周波数別推定方向情報生成部１５を動作の主体として説明する。

図８において、まず、周波数別推定方向情報生成部１５は、信号選択部１２によって選択された２つの入力信号と、それらの入力信号の相対遅延時間とを入力する（ステップＳ１２１）。

次に、周波数別推定方向情報生成部１５は、２つの入力信号を周波数領域信号（変換信号とも呼ぶ）に変換する（ステップＳ１２２）。

次に、周波数別推定方向情報生成部１５は、変換信号を用いてクロススペクトルを算出する（ステップＳ１２３）。

次に、周波数別推定方向情報生成部１５は、クロススペクトルを用いて平均クロススペクトルを算出する（ステップＳ１２４）。

次に、周波数別推定方向情報生成部１５は、平均クロススペクトルを用いて分散を算出する（ステップＳ１２５）。

次に、周波数別推定方向情報生成部１５は、平均クロススペクトルおよび分散を用いて周波数別クロススペクトルを算出する（ステップＳ１２６）。

次に、周波数別推定方向情報生成部１５は、周波数別クロススペクトルを用いて周波数別相互相関関数を算出する（ステップＳ１２７）。

次に、周波数別推定方向情報生成部１５は、周波数別相互相関関数と相対遅延時間とを用いて周波数別推定方向情報を算出する（ステップＳ１２８）。

そして、周波数別推定方向情報生成部１５は、周波数別推定方向情報を統合部１７に出力する（ステップＳ１２９）。

以上が、周波数別推定方向情報生成部１５の動作についての説明である。

〔周波数別クロススペクトル生成〕
次に、図９のフローチャートを用いて、波源方向推定装置１０の周波数別推定方向情報生成部１５に含まれる周波数別クロススペクトル生成部１５５の動作について説明する。図９のフローチャートの処理は、図８のフローチャートのステップＳ１２５を細分化したものである。なお、図９のフローチャートに沿った説明においては、周波数別クロススペクトル生成部１５５を動作の主体として説明する。

図９において、まず、周波数別クロススペクトル生成部１５５は、平均計算部１５３から平均クロススペクトルを入力し、分散計算部１５４から分散を入力する（ステップＳ１３１）。

次に、周波数別クロススペクトル生成部１５５は、平均クロススペクトルを用いて周波数別基本クロススペクトルを算出する（ステップＳ１３２）。

また、周波数別クロススペクトル生成部１５５は、分散を用いてカーネル関数スペクトルを算出する（ステップＳ１３３）。なお、ステップＳ１３２の処理と、ステップＳ１３３の処理とは、並行して行ってもよいし、逐次的に行ってもよい。

次に、周波数別クロススペクトル生成部１５５は、周波数別基本クロススペクトルとカーネル関数スペクトルとの積を計算し、周波数別クロススペクトルを算出する（ステップＳ１３４）。

そして、周波数別クロススペクトル生成部１５５は、算出した周波数別クロススペクトルを逆変換部１５６に出力する（ステップＳ１３５）。

以上が、周波数別クロススペクトル生成部１５５の動作についての説明である。

以上のように、本実施形態の波源方向推定装置は、複数の入力部、信号選択部、相対遅延時間計算部、少なくとも一つの周波数別推定方向情報生成部、および統合部を備える。複数の入力部は、複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する。信号選択部は、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを少なくとも二つ選択する。相対遅延時間計算部は、入力信号のペアを構成する少なくとも二つの入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する。少なくとも一つの周波数別推定方向情報生成部は、入力信号のペアと相対遅延時間とを用いて波動の波源の推定方向情報を周波数ごとに生成する。統合部は、周波数別推定方向情報生成部によって生成される周波数ごとの推定方向情報を統合する。

例えば、信号選択部は、センサの間隔に基づいて、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを選択する。

例えば、相対遅延時間計算部は、一つの入力信号のペアの供給元であるセンサのペアに対する波源探索方向を基準として、信号選択部によって選択された全ての入力信号のペアの相対遅延時間を、基準である波源探索方向の関数として計算する。

例えば、周波数別推定方向情報生成部は、変換部、クロススペクトル計算部、平均計算部、分散計算部、周波数別クロススペクトル生成部、逆変換部、および推定方向情報計算部を有する。変換部は、ペアをなす少なくとも２つの入力信号を周波数領域の変換信号に変換する。クロススペクトル計算部は、変換手段によって変換された変換信号を用いてクロススペクトルを計算する。平均計算部は、クロススペクトル計算部によって算出されるクロススペクトルを用いて平均クロススペクトルを計算する。分散計算部は、平均計算部によって算出される平均クロススペクトルを用いて分散を計算する。周波数別クロススペクトル生成部は、平均計算部によって算出される平均クロススペクトルと、分散計算部によって算出される分散とを用いて周波数別クロススペクトルを計算する。逆変換部は、周波数別クロススペクトル生成部によって算出される周波数別クロススペクトルを逆変換して周波数別相互相関関数を計算する。推定方向情報計算部は、逆変換部によって算出される周波数別相互相関関数と相対遅延時間とを用いて周波数別推定周波数ごとに推定方向情報を計算する。

例えば、周波数別クロススペクトル生成部は、周波数別基本クロススペクトル計算部、カーネル関数スペクトル生成部、および乗算部を含む。周波数別基本クロススペクトル計算部は、平均計算部から平均クロススペクトルを取得し、取得した平均クロススペクトルを用いて周波数別基本クロススペクトルを計算する。カーネル関数スペクトル生成部は、分散計算部から分散を取得し、取得した分散を用いてカーネル関数スペクトルを計算する。乗算部は、周波数別基本クロススペクトル計算部によって算出される周波数別基本クロススペクトルと、カーネル関数スペクトル生成部によって算出されるカーネル関数スペクトルとの積を計算して周波数別クロススペクトルを計算する。

例えば、統合部は、複数の周波数ごとの推定方向情報を複数の入力信号のペアについて統合した周波数別統合推定方向情報を算出する。そして、統合部は、算出した周波数別統合推定方向情報を全ての周波数について統合して統合推定方向情報を算出する。

例えば、統合部は、複数の周波数ごとの推定方向情報を全ての周波数について統合した入力信号組み合わせ別統合推定方向情報を算出する。統合部は、算出した入力信号組み合わせ別統合推定方向情報を全ての入力信号の組み合わせについて統合して統合推定方向情報を算出する。

例えば、波源方向推定装置は、複数の入力部のそれぞれに対応して配置されるセンサを備える。

本実施形態の波源方向推定装置は、マイクペア間の相互相関関数から推定方向情報を求め、複数のマイクペア間で推定方向情報を統合する。その結果、本実施形態の波源方向推定装置によれば、マイクペア間で位相が偶然揃うことに起因して発生する、音源方向以外の方向における推定方向情報の偽ピークが小さくなり、虚像音源の誤推定の発生を低減し、高精度に音源の方向を推定することができる。

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

本実施形態の波源方向推定装置によって推定される統合推定方向情報は、様々な形で利用できる。例えば、統合推定方向情報が複数のピークを有する場合には、各ピークを到来方向とする複数の音源が存在すると推定される。したがって、統合推定方向情報を用いれば、各音源の方向を同時に推定できるだけでなく、音源の数を推定できる。

（第２の実施形態）
次に、本発明の第２の実施形態に係る波源方向推定装置について図面を参照しながら説明する。本実施形態の波源方向推定装置は、第１の実施形態の波源方向推定装置に波源方向算出部を追加した構成を有する。

図１０は、本実施形態の波源方向推定装置２０の構成を表すブロック図である。波源方向推定装置２０は、入力端子２１、信号選択部２２、相対遅延時間計算部２３、周波数別推定方向情報生成部２５、統合部２７、および波源方向算出部２８を備える。なお、入力端子２１、信号選択部２２、相対遅延時間計算部２３、周波数別推定方向情報生成部２５、および統合部２７は、第１の実施形態の波源方向推定装置１０の対応する構成と同様であるため、詳細な説明を省略する。

〔波源方向算出部〕
波源方向算出部２８には、統合部２７から統合推定方向情報が入力される。波源方向算出部２８は、統合推定方向情報を用いて、波源方向を計算する。波源方向算出部２８は、計算した波源方向を出力する。

波源方向算出部２８における波源方向の計算方法を以下で詳細に説明する。統合部２７から入力された統合推定方向情報は、ピークが大きいほど信頼性（音源の存在可能性）が高い。そのため、例えば音源数が１つと事前に仮定できる場合には、波源方向算出部２８は、統合推定方向情報が最大となる方向を推定方向として出力する。ここで、統合部２７から入力された統合推定方向情報をＨ（θ，ｎ）とする。波源方向算出部２８は、以下の式３９を用いて、統合推定方向情報Ｈ（θ，ｎ）が最大値を取るような統合推定方向情報Ｈ（θ，ｎ）の引数を要素とする集合を波源方向Θとして算出できる。なお、式３９において、θは全ての波源方向または波源方向候補を表す。

また、統合推定方向情報のピークが閾値を超える場合には、波源方向算出部２８は、閾値を超えるピークを有する方向を音源とみなし、閾値を超える方向を推定方向として出力することもできる。

また、本実施形態の波源方向推定装置は、一定時間Ｔごとに統合推定方向情報が最大になる時刻に対応する方向を音源方向として推定することもできる。ただし、音源の方向は、一定時間Ｔの間で変化しない、あるいは変化の大きさが無視できるほど小さいと仮定する。このように仮定することによって、波源方向の推定精度を向上させることができる。

以上のように、本実施形態の波源方向推定装置は、統合手段によって算出される統合推定方向情報に基づいて波動の波源方向を算出する波源方向算出手段を備える。例えば、波源方向算出手段は、一定時間ごとに統合推定方向情報が最大になる時刻に対応する方向を波源方向として算出する。本実施形態の波源方向推定装置によれば、虚像音源を誤推定することなく、高精度に音源の方向を推定することができる。

（ハードウェア）
ここで、各実施形態に係る波源方向推定装置の処理を実行するハードウェア構成について、図１１の情報処理装置９０を一例として挙げて説明する。なお、図１１の情報処理装置９０は、各実施形態の波源方向推定装置の処理を実行するための構成例であって、本発明の範囲を限定するものではない。

図１１のように、情報処理装置９０は、プロセッサ９１、主記憶装置９２、補助記憶装置９３、入出力インターフェース９５および通信インターフェース９６を備える。図１１においては、インターフェースをＩ／Ｆ（Interface）と略して表記する。プロセッサ９１、主記憶装置９２、補助記憶装置９３、入出力インターフェース９５および通信インターフェース９６は、バス９９を介して互いにデータ通信可能に接続される。また、プロセッサ９１、主記憶装置９２、補助記憶装置９３および入出力インターフェース９５は、通信インターフェース９６を介して、インターネットやイントラネットなどのネットワークに接続される。

プロセッサ９１は、補助記憶装置９３等に格納されたプログラムを主記憶装置９２に展開し、展開されたプログラムを実行する。本実施形態においては、情報処理装置９０にインストールされたソフトウェアプログラムを用いる構成とすればよい。プロセッサ９１は、本実施形態に係る波源方向推定装置による処理を実行する。

主記憶装置９２は、プログラムが展開される領域を有する。主記憶装置９２は、例えばＤＲＡＭ（Dynamic Random Access Memory）などの揮発性メモリとすればよい。また、ＭＲＡＭ（Magnetoresistive Random Access Memory）などの不揮発性メモリを主記憶装置９２として構成・追加してもよい。

補助記憶装置９３は、種々のデータを記憶する。補助記憶装置９３は、ハードディスクやフラッシュメモリなどのローカルディスクによって構成される。なお、種々のデータを主記憶装置９２に記憶させる構成とし、補助記憶装置９３を省略することも可能である。

入出力インターフェース９５は、情報処理装置９０と周辺機器とを接続するためのインターフェースである。通信インターフェース９６は、規格や仕様に基づいて、インターネットやイントラネットなどのネットワークを通じて、外部のシステムや装置に接続するためのインターフェースである。入出力インターフェース９５および通信インターフェース９６は、外部機器と接続するインターフェースとして共通化してもよい。

情報処理装置９０には、必要に応じて、キーボードやマウス、タッチパネルなどの入力機器を接続するように構成してもよい。それらの入力機器は、情報や設定の入力に使用される。なお、タッチパネルを入力機器として用いる場合は、表示機器の表示画面が入力機器のインターフェースを兼ねる構成とすればよい。プロセッサ９１と入力機器との間のデータ通信は、入出力インターフェース９５に仲介させればよい。

また、情報処理装置９０には、情報を表示するための表示機器を備え付けてもよい。表示機器を備え付ける場合、情報処理装置９０には、表示機器の表示を制御するための表示制御装置（図示しない）が備えられていることが好ましい。表示機器は、入出力インターフェース９５を介して情報処理装置９０に接続すればよい。

また、情報処理装置９０には、必要に応じて、ディスクドライブを備え付けてもよい。ディスクドライブは、バス９９に接続される。ディスクドライブは、プロセッサ９１と図示しない記録媒体（プログラム記録媒体）との間で、記録媒体からのデータ・プログラムの読み出し、情報処理装置９０の処理結果の記録媒体への書き込みなどを仲介する。記録媒体は、例えば、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）などの光学記録媒体で実現できる。また、記録媒体は、ＵＳＢ（Universal Serial Bus）メモリやＳＤ（Secure Digital）カードなどの半導体記録媒体や、フレキシブルディスクなどの磁気記録媒体、その他の記録媒体によって実現してもよい。

以上が、各実施形態に係る波源方向推定装置を可能とするためのハードウェア構成の一例である。なお、図１１のハードウェア構成は、各実施形態に係る波源方向推定装置の演算処理を実行するためのハードウェア構成の一例であって、本発明の範囲を限定するものではない。また、各実施形態に係る波源方向推定装置に関する処理をコンピュータに実行させるプログラムも本発明の範囲に含まれる。さらに、各実施形態に係るプログラムを記録したプログラム記録媒体も本発明の範囲に含まれる。

各実施形態の波源方向推定装置の構成要素は、任意に組み合わせることができる。また、各実施形態の波源方向推定装置の構成要素は、ソフトウェアによって実現してもよいし、回路によって実現してもよい。

以上、実施形態を参照して本発明を説明してきたが、本発明は上記実施形態に限定されるものではない。本発明の構成や詳細には、本発明のスコープ内で当業者が理解し得る様々な変更をすることができる。
上記の実施形態の一部又は全部は、以下の付記のようにも記載されうるが、以下には限られない。
（付記１）
複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する複数の入力手段と、
複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択する信号選択手段と、
前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する相対遅延時間計算手段と、
前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに生成する少なくとも一つの周波数別推定方向情報生成手段と、
前記周波数別推定方向情報生成手段によって生成される周波数ごとの前記推定方向情報を統合する統合手段とを備える波源方向推定装置。
（付記２）
前記信号選択手段は、
前記センサの間隔に基づいて、複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを選択する付記１に記載の波源方向推定装置。
（付記３）
前記相対遅延時間計算手段は、
一つの前記入力信号のペアの供給元である前記センサのペアに対する前記波源探索方向を基準として、前記信号選択手段によって選択された全ての入力信号のペアの相対遅延時間を、基準である前記波源探索方向の関数として計算する付記１または２に記載の波源方向推定装置。
（付記４）
前記周波数別推定方向情報生成手段は、
ペアをなす少なくとも２つの前記入力信号を周波数領域の変換信号に変換する変換手段と、
前記変換手段によって変換された前記変換信号を用いてクロススペクトルを計算するクロススペクトル計算手段と、
前記クロススペクトル計算手段によって算出される前記クロススペクトルを用いて平均クロススペクトルを計算する平均計算手段と、
前記平均計算手段によって算出される平均クロススペクトルを用いて分散を計算する分散計算手段と、
前記平均計算手段によって算出される前記平均クロススペクトルと、前記分散計算手段によって算出される分散とを用いて周波数別クロススペクトルを計算する周波数別クロススペクトル生成手段と、
前記周波数別クロススペクトル生成手段によって算出される前記周波数別クロススペクトルを逆変換して周波数別相互相関関数を計算する逆変換手段と、
前記逆変換手段によって算出される前記周波数別相互相関関数と前記相対遅延時間とを用いて周波数ごとに前記推定方向情報を計算する周波数別推定方向情報計算手段とを有する付記１乃至３のいずれか一項に記載の波源方向推定装置。
（付記５）
前記周波数別クロススペクトル生成手段は、
前記平均計算手段から前記平均クロススペクトルを取得し、取得した前記平均クロススペクトルを用いて周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
前記分散計算手段から前記分散を取得し、取得した前記分散を用いてカーネル関数スペクトルを計算するカーネル関数スペクトル生成手段と、
前記周波数別基本クロススペクトル計算手段によって算出される前記周波数別基本クロススペクトルと、前記カーネル関数スペクトル生成手段によって算出される前記カーネル関数スペクトルとの積を計算して前記周波数別クロススペクトルを計算する乗算手段とを含む付記４に記載の波源方向推定装置。
（付記６）
前記統合手段は、
複数の周波数ごとの前記推定方向情報を複数の前記入力信号のペアについて統合した周波数別統合推定方向情報を算出し、算出した前記周波数別統合推定方向情報を全ての周波数について統合して統合推定方向情報を算出する付記１乃至５のいずれか一項に記載の波源方向推定装置。
（付記７）
前記統合手段は、
複数の周波数ごとの前記推定方向情報を全ての周波数について統合した入力信号組み合わせ別統合推定方向情報を算出し、算出した前記入力信号組み合わせ別統合推定方向情報を全ての入力信号の組み合わせについて統合して統合推定方向情報を算出する付記１乃至５のいずれか一項に記載の波源方向推定装置。
（付記８）
前記統合手段によって算出される前記統合推定方向情報に基づいて前記波動の波源方向を算出する波源方向算出手段を備える付記１乃至７のいずれか一項に記載の波源方向推定装置。
（付記９）
前記波源方向算出手段は、
一定時間ごとに前記統合推定方向情報が最大になる時刻に対応する方向を波源方向として算出する付記８に記載の波源方向推定装置。
（付記１０）
複数の前記入力手段のそれぞれに対応して配置される前記センサを備える付記１乃至９のいずれか一項に記載の波源方向推定装置。
（付記１１）
情報処理装置が、
複数のセンサによって取得された波動から変換される電気信号を入力信号として取得し、
複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択し、
前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出し、
前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成し、
周波数ごとの前記推定方向情報を統合する波源方向推定方法。
（付記１２）
複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する処理と、
複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択する処理と、
前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する処理と、
前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成する処理と、
周波数ごとの前記推定方向情報を統合する処理とコンピュータに実行させるプログラムを記録させたプログラム記録媒体。

１０、２０ 波源方向推定装置
１１、２１ 入力端子
１２、２２ 信号選択部
１３、２３ 相対遅延時間計算部
１５、２５ 周波数別推定方向情報生成部
１７、２７ 統合部
２８ 波源方向算出部
１５１ 変換部
１５２ クロススペクトル計算部
１５３ 平均計算部
１５４ 分散計算部
１５５ 周波数別クロススペクトル生成部
１５６ 逆変換部
１５７ 周波数別推定方向情報計算部
５５１ 周波数別基本クロススペクトル計算部
５５２ カーネル関数スペクトル生成部
５５３ 乗算部

Claims (12)

1. 複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する複数の入力手段と、
   複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択する信号選択手段と、
   前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する相対遅延時間計算手段と、
   前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに生成する少なくとも一つの周波数別推定方向情報生成手段と、
   前記周波数別推定方向情報生成手段によって生成される周波数ごとの前記推定方向情報を統合する統合手段とを備える波源方向推定装置。
2. 前記信号選択手段は、
   前記センサの間隔に基づいて、複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを選択する請求項１に記載の波源方向推定装置。
3. 前記相対遅延時間計算手段は、
   一つの前記入力信号のペアの供給元である前記センサのペアに対する前記波源探索方向を基準として、前記信号選択手段によって選択された全ての入力信号のペアの相対遅延時間を、基準である前記波源探索方向の関数として計算する請求項１または２に記載の波源方向推定装置。
4. 前記周波数別推定方向情報生成手段は、
   ペアをなす少なくとも２つの前記入力信号を周波数領域の変換信号に変換する変換手段と、
   前記変換手段によって変換された前記変換信号を用いてクロススペクトルを計算するクロススペクトル計算手段と、
   前記クロススペクトル計算手段によって算出される前記クロススペクトルを用いて平均クロススペクトルを計算する平均計算手段と、
   前記平均計算手段によって算出される平均クロススペクトルを用いて分散を計算する分散計算手段と、
   前記平均計算手段によって算出される前記平均クロススペクトルと、前記分散計算手段によって算出される分散とを用いて周波数別クロススペクトルを計算する周波数別クロススペクトル生成手段と、
   前記周波数別クロススペクトル生成手段によって算出される前記周波数別クロススペクトルを逆変換して周波数別相互相関関数を計算する逆変換手段と、
   前記逆変換手段によって算出される前記周波数別相互相関関数と前記相対遅延時間とを用いて周波数ごとに前記推定方向情報を計算する周波数別推定方向情報計算手段とを有する請求項１乃至３のいずれか一項に記載の波源方向推定装置。
5. 前記周波数別クロススペクトル生成手段は、
   前記平均計算手段から前記平均クロススペクトルを取得し、取得した前記平均クロススペクトルを用いて周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
   前記分散計算手段から前記分散を取得し、取得した前記分散を用いてカーネル関数スペクトルを計算するカーネル関数スペクトル生成手段と、
   前記周波数別基本クロススペクトル計算手段によって算出される前記周波数別基本クロススペクトルと、前記カーネル関数スペクトル生成手段によって算出される前記カーネル関数スペクトルとの積を計算して前記周波数別クロススペクトルを計算する乗算手段とを含む請求項４に記載の波源方向推定装置。
6. 前記統合手段は、
   複数の周波数ごとの前記推定方向情報を複数の前記入力信号のペアについて統合した周波数別統合推定方向情報を算出し、算出した前記周波数別統合推定方向情報を全ての周波数について統合して統合推定方向情報を算出する請求項１乃至５のいずれか一項に記載の波源方向推定装置。
7. 前記統合手段は、
   複数の周波数ごとの前記推定方向情報を全ての周波数について統合した入力信号組み合わせ別統合推定方向情報を算出し、算出した前記入力信号組み合わせ別統合推定方向情報を全ての入力信号の組み合わせについて統合して統合推定方向情報を算出する請求項１乃至５のいずれか一項に記載の波源方向推定装置。
8. 前記統合手段によって算出される前記統合推定方向情報に基づいて前記波動の波源方向を算出する波源方向算出手段を備える請求項７に記載の波源方向推定装置。
9. 前記波源方向算出手段は、
   一定時間ごとに前記統合推定方向情報が最大になる時刻に対応する方向を波源方向として算出する請求項８に記載の波源方向推定装置。
10. 複数の前記入力手段のそれぞれに対応して配置される前記センサを備える請求項１乃至９のいずれか一項に記載の波源方向推定装置。
11. 情報処理装置が、
    複数のセンサによって取得された波動から変換される電気信号を入力信号として取得し、
    複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択し、
    前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出し、
    前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成し、
    周波数ごとの前記推定方向情報を統合する波源方向推定方法。
12. 複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する処理と、
    複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択する処理と、
    前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する処理と、
    前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成する処理と、
    周波数ごとの前記推定方向情報を統合する処理とコンピュータに実行させるプログラムを記録させたプログラム記録媒体。

Images