JP6769495B2 - 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置 - Google Patents

Description

本発明は、相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置に関するものである。

上記技術分野において、特許文献１には、２つの音響信号からクロススペクトルを算出し、クロススペクトルの位相成分の時間変動が少ない周波数成分を音声成分とし、時間変動の大きい部分を雑音成分として、雑音成分を抑圧する技術が開示されている。また、上記技術分野において、特許文献２には、２つの受信信号のクロススペクトルが示すコヒーレンスレベルの高い周波数帯における位相の周波数勾配から到来電波の時間差を算出して電波の到来方向を探知することが開示されている。さらに、上記技術分野において、非特許文献１および非特許文献２には、２つのマイクロフォンの受音信号を用いて、音源(音波の発生源、発生場所)の方向を推定する方法が開示されている。具体的には、２つの受音信号から、受音信号間の相互相関関数を求める。そして、相互相関関数が最大値を与える時間差を音波の到来時間差として算出することで、音波の到来方向を推定する。

特開２０１１−０３３７１７号公報 特開２００９−２８７９４２号公報

C.H.Knapp and G.C.Carter, "The generalized correlation method for estimation of time delay," IEEE Trans. Acoustics, Speech and Signal Processing, vol.24, no. 4, pp. 320-327, Aug. 1976. J.P. Ianniello, "Time delay estimation via cross-correlation in the presence of large estimation errors," IEEE Trans. Acoustics, Speech and Signal Processing, vol.30, no. 6, pp. 998-1003, Dec. 1982.

しかしながら、上記文献に記載の技術では、周囲のノイズレベルが高い環境の場合には、明瞭なピークを持つ相関関数が生成できなかった。

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

上記目的を達成するため、本発明に係る相関関数生成装置は、
波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換手段と、
複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
前記クロススペクトルの位相の分散として円周分散を計算する分散計算手段と、
前記円周分散に基づいて算出された重みによって重み付された前記クロススペクトルに基づいて、相関関数を計算する相関関数計算手段と、
を備える。

上記目的を達成するため、本発明に係る相関関数生成方法は、
波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
複数の前記入力信号取得ステップで取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
前記クロススペクトルの位相の分散として円周分散を計算する分散計算ステップと、
前記円周分散に基づいて算出された重みによって重み付された前記クロススペクトルに基づいて、相関関数を計算して生成する相関関数計算ステップと、
を含む

上記目的を達成するため、本発明に係る相関関数生成プログラムは、
波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
複数の前記入力信号取得ステップで取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
前記クロススペクトルの位相の分散として円周分散を計算する分散計算ステップと、
前記円周分散に基づいて算出された重みによって重み付された前記クロススペクトルに基づいて、相関関数を計算して生成する相関関数計算ステップと、
をコンピュータに実行させる。

上記目的を達成するため、本発明に係る波源方向推定装置は、
上記相関関数生成装置を有し、
前記相関関数生成装置により生成された相関関数に基づいて、前記波源の方向を推定する。

本発明によれば、周囲のノイズレベルが高い環境の場合であっても、明瞭なピークを持つ相関関数を生成できる。

本発明の第１実施形態に係る相関関数生成装置の構成を示すブロック図である。 本発明の第２実施形態に係る相関関数生成部を含む波源方向推定装置の構成を説明するためのブロック図である。 本発明の第２実施形態に係る波源方向推定装置の備える相関関数計算部の構成を示すブロック図である。 本発明の第２実施形態に係る波源方向推定装置の備える相関関数テーブルの一例を説明する図である。 本発明の第２実施形態に係る波源方向推定装置の備える変換関数テーブルの一例を説明する図である。 本発明の第２実施形態に係る波源方向推定装置の備える分散テーブルの一例を説明する図である。 本発明の第２実施形態に係る波源方向推定装置のハードウェア構成を説明するブロック図である。 本発明の第２実施形態に係る波源方向推定装置の処理手順を説明するフローチャートである。 本発明の第３実施形態に係る相関関数生成部を含む波源方向推定装置の構成を説明するためのブロック図である。 本発明の第４実施形態に係る相関関数生成部を含む波源方向推定装置の構成を説明するためのブロック図である。 本発明の第５実施形態に係る波源方向推定装置の備える相関関数計算部の構成を説明するためのブロック図である。 本発明の第５実施形態に係る波源方向推定装置の備える周波数別クロススペクトル生成の構成を説明するためのブロック図である。 本発明の第５実施形態に係る波源方向推定装置の備える統合相関関数計算部の構成を説明するためのブロック図である。 本発明の第５実施形態に係る波源方向推定装置により計算された周波数別クロススペクトルおよび周波数別相関関数を説明する図である。 本発明の第５実施形態に係る波源方向推定装置によりカーネル関数スペクトルが乗じられた周波数別クロススペクトルと周波数別相関関数との関係を説明する図である。 本発明の第５実施形態に係る波源方向推定装置においてカーネル関数の有無と統合相関関数との関係を説明する図である。 本発明の第５実施形態に係る波源方向推定装置においてカーネル関数スペクトルの幅の違いと統合相関関数との関係を説明する図である。 本発明の第６実施形態に係る波源方向推定装置の備える周波数別クロススペクトル生成の構成を説明するブロック図である。 本発明の第７実施形態に係る統合相関関数計算部の構成を説明するためのブロック図である。 本発明の第８実施形態に係る波源方向推定システムの構成を説明するための図である。 本発明の第８実施形態に係る波源方向推定システムの表示部で表示した画像の一例を示す図である。

以下に、本発明を実施するための形態について、図面を参照して、例示的に詳しく説明記載する。ただし、以下の実施の形態に記載されている、構成、数値、処理の流れ、機能要素などは一例に過ぎず、その変形や変更は自由であって、本発明の技術範囲を以下の記載に限定する趣旨のものではない。

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

［第１実施形態］
本発明の第１実施形態としての相関関数生成装置１００について、図１を用いて説明する。相関関数生成装置１００は、明瞭なピークを持つ相関関数を生成する装置である。

図１に示すように、相関関数生成装置１００は、複数の入力信号取得部１０１と、変換部１０２と、クロススペクトル計算部１０３と、分散計算部１０４と、相関関数計算部１０５と、を含む。

入力信号取得部１０１は、波源で発生した波を入力信号として取得する。変換部１０２は、複数の入力信号取得部１０１で取得した複数の入力信号を変換して、複数の周波数領域信号を得る。クロススペクトル計算部１０３は、複数の周波数領域信号に基づいて、クロススペクトルを計算する。分散計算部１０４は、クロススペクトルの分散を計算する。相関関数計算部１０５は、クロススペクトルと分散とに基づいて、相関関数を計算して生成する。

本実施形態によれば、周囲のノイズレベルが高い環境の場合であっても、明瞭なピークを持つ相関関数を生成できる。

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

＜本実施形態の技術＞
次に本発明の第２実施形態としての波源方向推定装置２００について、図２乃至図６を用いて説明する。

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

波源方向推定装置２００は、入力端子２０_１、入力端子２０_２、変換部２０１、クロススペクトル計算部２０２、分散計算部２０３、相関関数計算部２０４、推定方向情報生成部２０５および相対遅延時間計算部２０６を備える。なお、変換部２０１、クロススペクトル計算部２０２、分散計算部２０３および相関関数計算部２０４は、相関関数生成部２１０を構成する。

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

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

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

（変換部２０１）
変換部２０１は、入力端子２０_１と入力端子２０_２とから供給された入力信号を変換し、クロススペクトル計算部２０２へ供給する。変換は、入力信号を複数の周波数成分に分解する目的で実施される。ここでは、代表的なフーリエ変換を用いた場合について説明する。

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

ここで、ｊは虚数単位（−１の平方根）、ｅｘｐは指数関数を表す。また、ｋは周波数ビン番号を表し、０以上Ｋ−１以下の整数である。以下、簡略化のため、ｋのことを周波数ビン番号ではなく、単に「周波数」と呼ぶ。

（クロススペクトル計算部２０２）
クロススペクトル計算部２０２は、変換部２０１から供給される変換信号に基づいて、クロススペクトルを計算し、分散計算部２０３と相関関数計算部２０４へ伝達する。クロススペクトル計算部２０２は、変換信号Ｘ_２（ｋ,ｎ）の複素共役と変換信号Ｘ_１（ｋ,ｎ）の積を計算する。変換信号のクロススペクトルをＳ_１２（ｋ,ｎ）とすると、クロススペクトルは次のように計算される。

ここで、ｃｏｎｊ（Ｘ_２（ｋ,ｎ））は、Ｘ_２（ｋ,ｎ）の複素共役を表す。また、式（２）で計算したクロススペクトルの代わりに、振幅成分で正規化したクロススペクトルを用いてもよい。このようなクロススペクトルは以下のように計算される。

（分散計算部２０３）
分散計算部２０３は、クロススペクトル計算部２０２から供給されたクロススペクトルを用いて、クロススペクトルの分散を計算し、相関関数計算部２０４へ伝達する。本実施形態では、過去に入力されたクロススペクトルから周波数ビンごとに分散を計算する例を説明する。周波数ビン単位ではなく、複数の周波数ビンを束ねたサブバンド単位や、全周波数ビンを用いて計算してもよい。また、それらを組み合わせて分散を計算してもよい。サブバンド単位で計算すると、分散を計算する帯域数が削減できるので、演算量が低減できる。

分散計算部２０３で計算される分散とは、クロススペクトルではなく、クロススペクトルの位相の分散である。ここで、クロススペクトルの位相とは、クロススペクトルを複素数と捉えたとき、クロススペクトルの偏角のことである。また、分散とは、位相の時系列がある基準値に対してどの程度散らばっているかを示す指標、位相の分布の広がりのことである。したがって、本実施形態における分散とは、確率論における分散の定義に限定されるものではなく、標準偏差も分散の一例であるとする。分散の代表例は、確率論の定義に従ったものである。この場合、第ｎフレームの周波数ビンｋにおけるクロススペクトルＳ_１２（ｋ,ｎ）の位相をａｒｇ（Ｓ_１２（ｋ,ｎ））とすると、過去Ｌ−１フレーム分のクロススペクトルの位相ａｒｇ（Ｓ_１２（ｋ,ｎ））,ａｒｇ（Ｓ_１２（ｋ,ｎ−１））,ａｒｇ（Ｓ_１２（ｋ,ｎ−２））,・・・,ａｒｇ（Ｓ_１２（ｋ,ｎーＬ＋１））から、周波数ビンｋの分散Ｖ_１２（ｋ,ｎ）は次のように計算される。

また、標準偏差の定義に基づいて計算した場合、分散は次のように計算される。

また、あらかじめ定めた基準値との２乗誤差の総和も分散の例であるとする。この場合、分散Ｖ_１２（ｋ,ｎ）は以下のように計算される。

ここで、ａｒｇ（Ｓ_０）は基準となる位相である。他には、２乗誤差の代わりに絶対値誤差を用いてもよい。

ここで、｜ｘ｜はｘの絶対値を表す。

位相が角度情報であるという点に着目すると、方向統計学（Directional Statistics）で定義されている分散(円周分散、circular variance)や標準偏差(円周標準偏差、circular standard deviation)を用いてもよい。円周分散を用いた場合、分散Ｖ１２（ｋ,ｎ）は以下のように計算される。

ここでＲ（ｋ,ｎ）は、クロススペクトルの平均であり、方向統計学では円周平均と呼ばれる。Ｒ（ｋ,ｎ）は、次のように計算される。

ただし、式（９）のＳ_１２（ｋ,ｎ）は、式（３）の「振幅成分で正規化したクロススペクトル」に限定される。クロススペクトルの平均は、式（８）の移動平均の代わりに、次のようにリーク積分を用いてもよい。

ここで、αは実数で、０＜α＜１を満足する。

式（８）の変形型である次の分散も有効である。

また、円周標準偏差を用いた場合、分散Ｖ_１２（ｋ,ｎ）は以下のように計算される。

ここで、ｌｎ（ｘ）はｘの自然対数を表す。この他には、円周標準偏差に類する統計量として知られる角状変位（angular deviation）を用いてもよい。

以上のとおり、方向統計学を用いた場合には、円周平均に基づいて各種分散を求めることが可能である。したがって、分散の計算方法としては、上記に示した数式に限定されることなく、円周平均の多次元関数、多項式、非線形関数、およびそれらの組み合わせを用いてもよい。

以上のとおりに求めた分散に対して補正処理を追加することも、分散の精度を向上させる上では有効である。つまり、分散の代わりに、補正された分散を相関関数計算部２０４へ伝達してもよい。補正は、時間周波数平面上で周辺に存在する分散（補正対象の周辺のフレームや周波数に存在する分散）を用いて行う。代表的な例は、分散の平均値を計算する方法である。この場合、補正された分散ＶＶ_１２（ｋ,ｎ）は次のように計算される。

ここで、Ｐ_１とＰ_２とは、それぞれあらかじめ定めた定数である。この他の例としては、隣接する複数の周波数やフレームの分散にばらつきがある場合、そのフレーム・周波数に目標信号が含まれている可能性が低いので、分散を大きくする。また、隣接周波数・フレームに分散値が突出して低い値が見つかった場合、分散を低くする。

（相関関数計算部２０４）
相関関数計算部２０４は、クロススペクトル計算部２０２から供給されたクロススペクトルと、分散計算部２０３から供給された分散に基づき、相関関数を計算し、推定方向情報生成部２０５に伝達する。次に図３を参照して、相関関数計算部２０４の詳細を説明する。

図３は、本実施形態に係る波源方向推定装置２００の備える相関関数計算部２０４の構成を示すブロック図である。相関関数計算部２０４は、重み付け部２４１および逆変換部２４２を備える。

（重み付け部２４１）
重み付け部２４１は、分散計算部２０３から供給された分散に基づいて、クロススペクトル計算部２０２から供給されたクロススペクトルに重み付けを行い、重み付けされたクロススペクトルを逆変換部２４２に伝達する。

重み付け部２４１では、分散が大きい場合には、クロススペクトルに目標音源の情報が含まれている可能性が低いので、クロススペクトルの値を小さくする。代表的な方法は、あらかじめ用意した写像関数を用いて、分散から重みを計算し、その重みを乗じる方法である。このとき、重み付けされたクロススペクトルＷＳ_１２（ｋ,ｎ）は次のように計算される。

ここで、Ｇ（ｋ,ｎ）は分散に基づいて算出された重みである。Ｇ（ｋ,ｎ）は、次のように計算される。

ここで、Ｖ_１２（ｋ,ｎ）は分散計算部２０３から供給された分散である。Ｖ_１２（ｋ,ｎ）の代わりに補正された分散であるＶＶ_１２（ｋ,ｎ）を用いてもよい。また、ａとｂはそれぞれ実数であり、ａ＞０を満足する。また、次式のように、単純に閾値判定に基づいて値を切り替える方法でもよい。

ここで、Ｖ_ｔｈは閾値であり、正の実数である。この他にも、線形写像関数や高次の多項式関数、非線形関数など、他の形で表される分散の関数をゲインの計算に用いることも可能である。

（逆変換部２４２）
逆変換部２４２は、重み付け部２４１から供給されたクロススペクトルの逆変換を求める。本実施形態では、変換部２０１でフーリエ変換を用いたので、逆変換には逆フーリエ変換を用いる方法について説明する。重み付け部２４１から供給されたクロススペクトルをＷＳ_１２（ｋ,ｎ）とすると、ＷＳ_１２（ｋ,ｎ）の逆変換により得られる相関関数ｓ_１２（τ,ｎ）は次のように計算される。

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

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

（推定方向情報生成部２０５）
推定方向情報生成部２０５は、相関関数計算部２０４から供給される相関関数と、相対遅延時間計算部２０６から供給される相対遅延時間と、に基づき、方向と相関値との対応関係を推定方向情報として出力する。相関関数をｓ_１２（τ,ｎ）、相対遅延時間τ(θ)とすると、推定方向情報Ｈ（θ,ｎ）は、次式で与えられる。

各方向に対して相関値が定まるので、基本的に相関値が高ければ、その方向に音源が存在する可能性が高いと判断できる。

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

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

図４Ａは、本実施形態に係る波源方向推定装置２００の備える相関関数テーブル４０１の一例を説明する図である。相関関数テーブル４０１は、入力信号４１１に関連付けて、周波数領域信号４１２、クロススペクトル４１３、分散４１４、相対遅延時間４１５および相関関数４１６を記憶する。波源方向推定装置２００は、例えば、相関関数テーブル４０１を参照して、入力信号４１１に対する相関関数４１６を算出し、推定方向情報を生成する。

図４Ｂは、本実施形態に係る波源方向推定装置２００の備える変換関数テーブル４０２の一例を説明する図である。変換関数テーブル４０２は、入力信号４１１に関連付けて変換関数４２１および逆変換関数４２２を記憶する。波源方向推定装置２００は、例えば、変換関数テーブル４０２を参照して、入力信号４１１を変換して、周波数領域信号４１２やクロススペクトル４１３などを算出する。

図４Ｃは、本実施形態に係る波源方向推定装置２００の備える分散テーブル４０３の一例を説明する図である。分散テーブル４０３は、クロススペクトル４２４に関連付けて分散手法４３１を記憶する。波源方向推定装置２００は、例えば、分散テーブル４０３を参照してクロススペクトル４２４の分散４１４を算出する。

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

ＲＡＭ５４０は、ＣＰＵ５１０が一時記憶のワークエリアとして使用するランダムアクセスメモリである。ＲＡＭ５４０には、本実施形態の実現に必要なデータを記憶する領域が確保されている。入力信号５４１は、入力端子２０_１および入力端子２０_２に入力された音などの信号である。周波数領域信号５４２は、入力信号を所定の変換方法で変換して得られた信号である。クロススペクトル５４３は、周波数領域信号５４２に基づいて計算される。分散５４４は、クロススペクトル５４３から計算されたクロススペクトル５４３の分散である。相対遅延時間５４５は、マイク間隔と音源方向とに基づいて決まる音波などの到達時間差である。相関関数５４６は、クロススペクトル５４３と分散５４４とから計算される。これらのデータは、例えば、相関関数テーブル４０１から展開されたデータである。なお、図示していないが、ＲＡＭ５４０は、変換関数テーブル４０２および分散テーブル４０３から展開されたデータを含んでいてもよい。

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

ストレージ５５０には、データベースや各種のパラメータ、あるいは本実施形態の実現に必要な以下のデータまたはプログラムが記憶されている。ストレージ５５０は、相関関数テーブル４０１、変換関数テーブル４０２および分散テーブル４０３を格納する。相関関数テーブル４０１は、図４Ａに示した、入力信号４１１と相関関数４１６などとの関係を管理するテーブルである。変換関数テーブル４０２は、図４Ｂに示した、入力信号４１１と変換関数４２１などとの関係を管理するテーブルである。分散テーブル４０３は、図４Ｃに示した、クロススペクトル４２４と分散手法４３１などとの関係を管理するテーブルである。ストレージ５５０は、さらに、変換モジュール５５１、クロススペクトル計算モジュール５５２、分散計算モジュール５５３、相関関数計算モジュール５５４、相対遅延時間計算モジュール５５５および推定方向情報生成モジュール５５６を格納する。

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

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

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

ステップＳ６０１において、波源方向推定装置２００は、入力信号を取得する。ステップＳ６０３において、波源方向推定装置２００は、入力信号を変換して周波数領域信号を取得する。ステップＳ６０５において、波源方向推定装置２００は、周波数領域信号からクロススペクトルを計算する。ステップＳ６０７において、波源方向推定装置２００は、クロススペクトルの位相の分散を計算する。ステップＳ６０９において、波源方向推定装置２００は、計算した分散に基づいて、クロススペクトルを重み付けする。ステップＳ６１１において、波源方向推定装置２００は、重み付けしたクロススペクトルを逆変換して、相関関数を計算する。ステップＳ６１３において、波源方向推定装置２００は、相対遅延時間を計算する。ステップＳ６１５において、相関関数と相対遅延時間とから推定方向情報を生成する。

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

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

クロススペクトルの位相から、２つのマイクへの音波の到達時間差、すなわち音源方向が得られる。したがって、位相の時系列（本実施形態では過去Ｌ−１フレームの位相）から得られた分散の大きさに基づいて、目標音と環境音との区別が可能になる。目標音源の方向の時間変化は小さいため、位相の分散は小さくなる。もし、目標音の発生時間が数秒程度であっても、例えば１０ミリ秒ごとに位相を求める場合には、その数秒から得られる位相の分散は小さくなる。一方で、環境雑音は方向の特定が困難なため、方向の時間変化は大きくなり、位相の分散も大きくなる。分散は周波数ごとに計算するので、目標音源が存在する周波数帯域を特定することが可能になる。

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

（平均計算部７０１）
平均計算部７０１は、クロススペクトル計算部２０２から供給されたクロススペクトルの平均を計算し、相関関数計算部２０４へ伝達する。本実施形態では、過去に入力されたクロススペクトルから周波数ビンごとに平均を計算する例を説明する。周波数ビン単位ではなく、複数の周波数ビンを束ねたサブバンド単位で計算してもよい。第ｎフレームの周波数ビンｋにおけるクロススペクトルをＳ_１２（ｋ,ｎ）とすると、過去Ｌ−１フレームから求めた平均クロススペクトルＳＳ_１２（ｋ,ｎ）は次のように計算される。

また、次のようにリーク積分を用いてもよい。

ここで、αは実数で、０＜α＜１を満足する。

＜本実施形態の効果＞
本実施形態によれば、クロススペクトルの代わりに平均クロススペクトルを相関関数計算部２０４に伝達する。このため、第２実施形態と比べると、位置の時間変動が小さな目標音源に対しては、相関関数のピーク位置、および推定方向情報の時間変動が小さくなるので、方向推定の精度が向上する。特に、自動車の走行音、航空機やドローンの飛行音、船舶の航行音(スクリュー音)など、目標物が発する音が継続的に生じる場合、クロススペクトルの平均化効果は高くなるので、相関関数のピーク位置、および推定方向情報が明確になる。

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

（分散計算部８０１）
分散計算部８０１は、平均計算部３０１から供給された平均クロススペクトルを用いて分散を計算し、相関関数計算部２０４へ伝達する。第３実施形態と比較すると、分散の計算にクロススペクトルではなく、平均クロススペクトルを用いている点が異なる。特に、方向統計学で知られている分散や標準偏差を計算する場合には、その計算過程において円周平均の計算を省略できる。平均クロススペクトルをＳＳ_１２（ｋ,ｎ）とすると、「クロススペクトルの位相の分散」の計算において円周分散を用いた場合、分散Ｖ_１２（ｋ,ｎ）は以下のように計算される。

式（２３）の変形型である次の分散も有効である。

また、円周標準偏差を用いた場合、分散Ｖ_１２（ｋ,ｎ）は以下のように計算される。

＜本実施形態の効果＞
本実施形態によれば、分散計算部は、平均クロススペクトル計算部で計算された平均クロススペクトルを利用して、分散を計算する。このため、分散計算部において平均クロススペクトルの計算が不要になる。したがって、第３実施形態よりも少ない計算量で分散を計算することができる。

［第５実施形態］
次に本発明の第５実施形態に係る波源方向推定装置について、図９乃至図１２Ｄを用いて説明する。本実施形態に係る波源方向推定装置は、上記第２実施形態の波源方向推定装置と比べると、相関関数計算部２０４に代えて、相関関数計算部９０４を有する点で異なる。図９は、本実施形態に係る相関関数計算部９０４の構成を説明するためのブロック図である。本実施形態に係る相関関数計算部９０４は、上記第２実施形態の相関関数計算部２０４と比べると、重み付け部２４１と逆変換部２４２に代えて、周波数別クロススペクトル生成部９４１_１，９４１_２,・・・，９４１_Ｋと統合相関関数計算部９４２とを有する点で異なる。その他の構成および動作は、第２実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。

図９は、相関関数計算部９０４のブロック図である。相関関数計算部９０４は、周波数別クロススペクトル生成部９４１_１，９４１_２,・・・，９４１_Ｋと統合相関関数計算部９４２とを備える。周波数別クロススペクトル生成部９４１_１，９４１_２,・・・，９４１_Ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ_１２（ｋ,ｎ）と分散計算部２０３から供給される分散ｐ（ｋ,ｎ）とを用いて、Ｓ_１２（ｋ,ｎ）の各周波数ｋに対応するクロススペクトルを計算し、周波数別クロススペクトルとして統合相関関数計算部９４２に伝達する。周波数別クロススペクトルの生成は、周波数成分ごとに相関関数を計算するために行われる。つまり、ある周波数ｋに対応する相関関数（周波数別相関関数と呼ぶ）を後段で求めるために、周波数別クロススペクトルを計算する。

相関関数計算部９０４は、統合相関関数計算部９４２および周波数別クロススペクトル生成部９４１_１，９４１_２,・・・，９４１_Ｋから供給される周波数別クロススペクトルに基づき統合相関関数を計算し、推定方向情報生成部２０５へ伝達する。

次に、図１０を参照して、ある周波数ｋの周波数別クロススペクトルを計算する周波数別クロススペクトル生成部９４１_ｋについて詳細に説明する。その後、図１１を参照して、統合相関関数計算部９４２の詳細について説明する。

＜周波数別クロススペクトル生成部９４１_ｋ＞
図１０は、周波数別クロススペクトル生成部９４１_ｋのブロック図である。周波数別クロススペクトル生成部９４１_ｋは、周波数別クロススペクトル生成部１０１１、カーネル関数スペクトル生成部１０１２、乗算部１０１３を有する。周波数別クロススペクトル生成部９４１_ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ_１２（ｋ,ｎ）と分散計算部２０３から供給される分散を用いて、Ｓ_１２（ｋ,ｎ）の周波数ｋに対応するクロススペクトルを計算し、周波数別クロススペクトルを統合相関関数計算部９４２へ伝達する。

（周波数別クロススペクトル生成部９４１）
周波数別クロススペクトル生成部９４１_１，９４１_２,・・・，９４１_Ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ_１２（ｋ,ｎ）を用いて、Ｓ_１２（ｋ,ｎ）の各周波数ｋに対応するクロススペクトルを計算し、周波数別クロススペクトルとして乗算部１０１３に伝達する。周波数別クロススペクトルは、周波数成分ごとに相関関数を計算するために行われる。つまり、ある周波数ｋに対応する相関関数（周波数別相関関数と呼ぶ）を後段で求めるために、周波数別クロススペクトルを計算する。

次に、ある周波数ｋの周波数別クロススペクトルを計算する周波数別クロススペクトル生成部９４１_ｋについて詳細に説明する。周波数別クロススペクトル生成部９４１_ｋでは、周波数ｋのクロススペクトルＳ_１２（ｋ,ｎ）に基づいて周波数別クロススペクトルを計算する際、位相成分と振幅成分をあらかじめ別々に求めたのちに統合する。周波数ｋの周波数別クロススペクトルＵｋ（ｗ,ｎ）をとし、その振幅成分を｜Ｕ_ｋ（ｗ,ｎ）｜、位相成分をａｒｇ（Ｕ_ｋ（ｗ,ｎ））とすると、次のような関係が成立する。

ここで、ｗは周波数を表し、０以上Ｗ−１以下の整数である。以下、周波数ｋのクロススペクトルＳ_１２（ｋ,ｎ）から、周波数別クロススペクトルの振幅成分|Ｕ_ｋ（ｗ,ｎ）|と位相成分ａｒｇ（Ｕ_ｋ（ｗ,ｎ））とを求める方法について説明する。

振幅成分|Ｕ_ｋ（ｗ,ｎ）|は、ｋを整数倍した周波数には、１．０を用いる。一方、周波数ｋの非定数倍した周波数の位相成分はゼロにする。これを数式で表現すると、振幅成分|Ｕ_ｋ（ｗ,ｎ）|は以下のように与えられる。

ここで、ｐは１以上Ｐ以下の整数である。波源方向推定を行うときに重要な情報は位相成分であるため、振幅成分はこのように適当な定数を用いる。その他には、１．０の代わりに|Ｓ_１２（ｋ,ｎ）|を用いてもよい。つまり、振幅成分|Ｕ_ｋ（ｗ,ｎ）|を次式のように求めてもよい。

位相成分ａｒｇ（Ｕ_ｋ（ｗ,ｎ））は、ｋを整数倍した周波数には、周波数ｋのクロススペクトルＳ_１２（ｋ,ｎ）を定数倍したものを用いる。例えば、周波数ｋ、２ｋ、３ｋ、４ｋの位相成分には、周波数ｋの位相成分ａｒｇ（Ｓ_１２（ｋ,ｎ））をそれぞれ同一の倍率で整数倍したもの、つまりａｒｇ（Ｓ_１２（ｋ,ｎ））、２ａｒｇ（Ｓ_１２（ｋ,ｎ））、３ａｒｇ（Ｓ_１２（ｋ,ｎ））、４ａｒｇ（Ｓ_１２（ｋ,ｎ））を用いる。一方、周波数ｋの非定数倍した周波数の位相成分はゼロにする。したがって、周波数ｋに対応する周波数別クロススペクトルの位相成分ａｒｇ（Ｕ_ｋ（ｗ,ｎ））は、以下のように計算される。

ここで、ｐは１以上Ｐ以下の整数である。また、Ｐは１よりも大きい整数である。

以上の方法で求めた振幅成分と位相成分とを、先に記載した式（２６）を用いて統合し、周波数ｋの周波数別クロススペクトルＵ_ｋ（ｗ,ｎ）を得る。

これまで説明してきた方法では、振幅成分と位相成分とを別々に求めてから周波数別スペクトルを得ていた。しかし、次に示す数式に示すようにクロススペクトルのべき乗を用いれば、振幅成分と位相成分とを求めることなく周波数別スペクトルＵ_ｋ（ｗ,ｎ）を求めることが可能である。

(カーネル関数スペクトル生成部１０１２)
カーネル関数スペクトル生成部１０１２は、分散計算部２０３から供給された分散に基づいてカーネル関数スペクトルを求め、乗算部１０１３へカーネル関数スペクトルを出力する。カーネル関数スペクトルとは、カーネル関数をフーリエ変換し、その絶対値を取ったものである。絶対値を取る代わりに二乗してもよい。また、絶対値の二乗としてもよい。本実施形態では、カーネル関数スペクトルをＧ（ｗ）、カーネル関数をｇ（τ）と表記して説明する。カーネル関数としては、ガウス関数を用いる。この時、ガウス関数は、次のような数式で与えられる。

ここで、ｇ_１,ｇ_２,ｇ_３は正の実数である。ｇ_１はガウス関数の大きさ、ｇ_２はガウス関数のピークの位置、ｇ_３はガウス関数の広がりを制御する。特に、ガウス関数の広がりを調整するｇ_３は、周波数別相関関数のピークの鋭さに大きな影響を与えるので重要である。式（３１）から分かるように、ｇ_３が大きくなればガウス関数の広がりは大きくなる。

他には、以下に示すロジスティック関数を用いてもよい。

ここで、ｇ_１とｇ_２とは正の実数である。ロジスティック関数は、ガウス関数と同様の形状をしているが、ガウス関数よりも裾が長いという性質を持つ。特に、ロジスティック関数の広がりを調整するｇ_５は、ガウス関数におけるｇ_３の場合と同様、周波数別相関関数のピークの鋭さに大きな影響を与える重要なパラメータである。この他にも、コサイン関数や一様関数を用いてもよい。

カーネル関数のパラメータのうち、カーネル関数の広がりに影響を与えるｇ_３やｇ_５は、分散計算部２０３から供給された分散に基づいて決定する。本実施形態では、これらのパラメータのことを広がり制御パラメータと呼び、ｑ（ｋ,ｎ）で表現する。したがって、カーネル関数がガウス関数の場合には、ｇ_３＝ｑ（ｋ,ｎ）である。分散が小さければ、周波数別相関関数のピークが鋭く、裾が狭くなるようにパラメータを変化させる。したがって、広がり制御パラメータを小さくする。

広がり制御パラメータの基本的な算出方法は、あらかじめ設定した写像関数を用いて、分散の値を広がり制御パラメータに変換する方法である。例えば、分散がある閾値を上回れば広がり制御パラメータを大きな値（例えば１０）に、閾値を下回れば小さな値（例えば０．０１）に設定する。この場合、分散をＶ_１２（ｋ,ｎ）、閾値をｐ_ｔｈとすると、第ｎフレームの周波数ビンｋにおける広がり制御パラメータｑ（ｋ,ｎ）は次のように計算される。

ここで、ｑ_１とｑ_２とは正の実数であり、ｑ_１＞ｑ_２を満足する。また、次のように一次関数を用いて計算してもよい。

ここで、ｑ_３とｑ_４とはそれぞれ実数であり、ｑ_３＞０を満足する。この他にも、線形写像関数や高次の多項式関数、非線形関数など、他の形で表される分散の関数を分散の計算に用いることも可能である。また、分散をそのまま広がり制御パラメータにしてもよい（ｑ（ｋ,ｎ）＝Ｖ_１２（ｋ,ｎ）とすると同じ)。

広がり制御パラメータを求める関数は、分散だけでなく、周波数ｋの関数としてもよい。例えば、周波数が高くなるにつれて小さくなる関数とする。このような代表例としてｋの逆数を用いる例を説明する。この場合、式（３３）の代わりに、次の関数を用いて広がり制御パラメータｑ（ｋ,ｎ）を計算する。

また、式（３４）の代わりとしては、次の関数を用いて広がり制御パラメータｑ（ｋ,ｎ）を計算する。

（乗算部１０１３）
乗算部１０１３は、周波数別クロススペクトル生成部１０１１から供給される周波数別クロススペクトルと、カーネル関数スペクトル生成部１０１２から供給されるカーネル関数スペクトルと、の積を計算し、計算によって得られた新たな周波数クロススペクトルを統合相関関数計算部９４２に伝達する。周波数別クロススペクトル生成部１０１１から供給される周波数別スペクトルをＵ_ｋ（ｗ,ｎ）、カーネル関数スペクトル生成部１０１２から供給されるカーネル関数スペクトルをＧ（ｗ）とすると、計算により得られる新たな周波数別クロススペクトルＵＭ_ｋ（ｗ,ｎ）は次のように計算される。

＜統合相関関数計算部＞
図１１は、統合相関関数計算部９４２のブロック図である。統合相関関数計算部９４２は、逆変換部１１２１_１,１１２１_２,・・・,１１２１_Ｋと、統合部１１２２を備える。

逆変換部１１２１_１,１１２１_２,・・・,１１２１_Ｋは、周波数別クロススペクトル生成部９４１_１,９４１_２,・・・,９４１_Ｋから供給される周波数別クロススペクトルの逆変換を行い、周波数別相関関数として統合部１１２２にそれぞれ伝達する。本実施形態では、変換部２０１でフーリエ変換を用いたので、逆変換には逆フーリエ変換を用いる方法について説明する。周波数別クロススペクトル生成部９４１_ｋから供給された周波数別クロススペクトルをＵＭ_ｋ（ｗ,ｎ）とすると、ＵＭ_ｋ（ｗ,ｎ）の逆変換により得られる周波数別相関関数ｕ_ｋ（τ,ｎ）は次のように計算される。

統合部１１２２は、逆変換部１１２１_１,１１２１_２,・・・,１１２１_Ｋから供給される周波数別相関関数を統合し、統合相関関数として推定方向情報生成部２０５へ伝達する。個別に求めた複数の周波数別相関関数を混合したり、重ね合わせたりすることにより、一つの相関関数を求める。統合方法に単純な和を用いた場合、統合部１１２２は、周波数別相関関数の総和を計算する。統合相関関数をｕ（τ,ｎ）とすると、ｕ（τ,ｎ）は次のように計算される。

また、総和ではなく総乗を用いてもよい。この場合、ｕ（τ,ｎ）は次のように計算される。

目標音の存在する周波数や、目的音のパワーが大きい周波数があらかじめ判明している場合には、その周波数に対応する周波数別相関関数だけ用いて、統合された相関関数を求めてもよい。また、重み付けという形で、統合における周波数別相関関数の影響度を制御してもよい。例えば目標音の存在する周波数の集合をΩとすると、周波数の選択によりｕ（τ,ｎ）を求める場合は、次のように計算される。

また、重み付けを用いる場合には、ｕ（τ,ｎ）は次のように計算される。

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

＜周波数別クロススペクトルの説明＞
以上の方法により周波数別クロススペクトルを計算すると、周波数別クロススペクトルを逆変換して得られる周波数別相関関数のピークが鋭くなり、相関関数のピーク位置が明確になる。相関関数のピーク位置に基づいて波源方向推定を行う本発明の場合には、ピークが鋭くなれば音源方向推定の精度が向上する。さらに、Ｐの値が大きくなるにつれて、すなわちｋを整数倍した周波数の成分が増えるにつれて、相関関数のピークは鋭くなる。図１２Ａに、この様子を示す。ここで、図中のＱは３よりも大きい整数である。Ｐ＝１の場合、つまり位相成分が一つしかない場合、これを逆変換して得られる相関関数は、このようにピーク位置が不明確な相関関数になる。Ｐが大きくなると、図１２Ａに示す通り、相関関数のピークは鋭くなる。

したがって、周波数別クロススペクトル生成部１０１１で得られる周波数別クロススペクトルは、ある周波数ｋのクロススペクトルを基に「周波数ｋを整数倍した周波数ｐｋの位相成分に、周波数ｋの位相成分ａｒｇ（Ｓ_１２（ｋ,ｎ））をｐ倍したものを割り当てたもの」として定義される。ここでｐは１以上の整数である。つまり、周波数別クロススペクトルは、その位相成分ａｒｇ（Ｕｋ（ｗ,ｎ））が少なくとも以下の式を満足するものとして定義される。

加えてｐは、ｐ＝１と２、ｐ＝１と３、ｐ＝２と３など、２つ以上に限定される。ｐが１のみの場合には、周波数ｋの成分だけを抜き出して周波数別クロススペクトルが生成されることになるが、方向推定精度は従来技術と同等となり、方向推定の高精度を達成できない。なお、図１２Ａで説明した通り、ｐ＝１，２，３，４,・・・とｐの数が増えれば、つまり周波数ｐｋの位相成分への割り当てが増えれば、周波数別相関関数のピークは鋭くなるため、方向推定の精度は向上する。

＜周波数別クロススペクトルにカーネル関数スペクトルを乗じる効果の説明＞
周波数別クロススペクトルにカーネル関数スペクトルを乗じることで、周波数別相関関数の形状を変更できる。図１２Ｂは、カーネル関数スペクトルが乗じられた周波数別クロススペクトルと、周波数別相関関数の関係を示している。比較のため、カーネル関数スペクトルを乗じる前の周波数別クロススペクトルも掲載する。図１２Ｂの左側図に示すように、カーネル関数スペクトルを乗じていなければ、高い周波数まで成分が存在するため、周波数別相関関数のピークは鋭くなる。一方で、図１２Ｂの中央図、右側図のようにカーネル関数スペクトルを乗じると、高い周波数の成分が減衰するので、周波数別相関関数のピークの鋭さは小さくなる。つまり、カーネル関数スペクトルのピークが鋭くなる（カーネル関数スペクトルの裾が狭くなる）につれて、周波数別相関関数のピークの鋭さは小さくなる。また、図１２Ｂの右側図のように、周波数別相関関数の裾が大きく広がると、隣接する山の裾が重なり、谷が浅い周波数別相関関数が得られる。

ここで、カーネル関数とカーネル関数スペクトルとの形状の関係について補足する。フーリエ変換の性質により、形状の関係は逆になる。カーネル関数のピークが鋭く、裾が狭くなるにつれて、カーネル関数スペクトルのピークは平坦に近づき、裾が広くなる。ガウス関数の広がりを調整するｇ_３との関係も含めて説明すると、ｇ_３が大きくなれば、ガウス関数の広がりは大きくなるが、そのスペクトルの広がりは小さくなる。

カーネル関数による周波数別相関関数の大きさ制御の効果について、図１２Ｃにおいて説明する。図１２Ｃは、カーネル関数の有無と統合相関関数との関係を示した図である。図１２Ｃに示すようなカーネル関数無しの場合（１２１０）、周波数別相関関数ｕ_１（τ,ｎ）〜ｕ_３（τ,ｎ）のピーク位置は近いが、ｕ_１（τ,ｎ）〜ｕ_３（τ,ｎ）の幅が細いため、統合時に大きなピークを形成できない。このため、ピークの位置が明確にならない。一方、図１２Ｃに示すようなカーネル関数有りの場合（１２２０）、周波数別相関関数の幅は太くなっているため、統合によりｕ_１（τ,ｎ）〜ｕ_３（τ,ｎ）は大きなピークを形成できる。このため、カーネル関数無しの場合（１２１０）に比べて、ピークの位置が明確になる。

さらに、カーネル関数による周波数別相関関数の大きさ制御の別の効果について、図１２Ｄを参照して説明する。図１２Ｄは、カーネル関数スペクトルの幅の違いと、統合相関関数と、の関係を示した図である。図１２Ｂの右側図に示すように、幅の広いカーネル関数スペクトルを用いると、相関関数の周期性により、谷が浅い周波数別相関関数が形成される。したがって、図１２Ｄの１２３０に示す通り、谷が浅い周波数別相関関数を統合すると、谷が浅い、つまりピークの目立たない統合相関関数が生成される。一方、図１２Ｂの中央図に示すように、幅の狭いカーネル関数スペクトルを用いた場合には、図１２Ｂの右側図よりも谷が深い周波数別相関関数が形成される。したがって、図１２Ｄの１２４０に示す通り、ピークが明確な統合相関関数が生成される。

本実施形態では、カーネル関数のフーリエ変換により得られたカーネル関数スペクトルと、周波数別クロススペクトルと、の積を計算しているが、フーリエ変換の性質により時間領域でも実現できる。周波数別クロススペクトル計算部ではなく、統合相関関数計算部９４２の中にある逆変換部１１２１_ｋの後段にカーネル関数を畳み込む「畳み込み演算部」を設け、逆変換部１１２１_ｋから供給される周波数別相関関数にカーネル関数を畳み込んでもよい。ただし、畳み込み演算は計算量が多いため、本実施形態のように周波数領域で積を計算するほうが効率的である。

＜本実施形態の効果＞
本実施形態によれば、周波数別クロススペクトルを逆変換して得られる周波数別相関関数のピークが鋭くなり、相関関数のピーク位置が明確になる。相関関数のピーク位置に基づいて波源方向推定を行う本実施形態の場合には、音源方向推定の精度が向上する。また、周波数別クロススペクトルにカーネル関数スペクトルを乗じることで、周波数別相関関数の形状を変更できる。その際、カーネル関数の広さは、分散の大きさに応じて変更するので、目標音源が存在する周波数の周波数別相関関数のピークを強調できる。したがって、周波数ごとの相関関数を求めず、また周波数別相関関数のピーク強調を行わない第２実施形態と比較すると、相関関数のピークが強調できるので、方向推定精度が向上する。

［第６実施形態］
次に本発明の第６実施形態に係る波源方向推定装置について、図１３を用いて説明する。本実施形態に係る波源方向推定装置は、上記第５実施形態の波源方向推定装置と比べると、周波数別クロススペクトル生成部９４１_ｋに代えて、周波数別クロススペクトル生成部１３４１_ｋを有する点で異なる。図１３は、本実施形態に係る周波数別クロススペクトル生成部１３４１_ｋの構成を説明するためのブロック図である。本実施形態に係る周波数別クロススペクトル生成部１３４１_ｋは、上記第５実施形態の周波数別クロススペクトル生成部９４１_ｋと比べると、カーネル関数スペクトル生成部１０１２に代えて、カーネル関数スペクトル記憶部１３１４とカーネル関数スペクトル選択部１３１２とを有する点で異なる。その他の構成および動作は、第５実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。

図１３は、周波数別クロススペクトル生成部１３４１_ｋのブロック図である。周波数別クロススペクトル生成部１３４１_ｋは、周波数別クロススペクトル生成部１０１１、カーネル関数スペクトル選択部１３１２、カーネル関数スペクトル記憶部１３１４、乗算部１０１３を有する。周波数別クロススペクトル生成部１０１１_ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ_１２（ｋ,ｎ）と分散計算部２０３から供給される分散を用いて、Ｓ_１２（ｋ,ｎ）の周波数ｋに対応するクロススペクトルを計算し、周波数別クロススペクトルを統合相関関数計算部９４２へ伝達する。

カーネル関数スペクトル記憶部１３１４は、カーネル関数スペクトルを記憶しており、カーネル関数スペクトル選択部１３１２へカーネル関数スペクトルを伝達する。カーネル関数スペクトル記憶部１３１４は、異なる分散に対応する複数のカーネル関数スペクトルを記憶する。例えば、０．１から０．１刻みで１０種類の分散（０．１,０．２,．．．,１．０）に対応するカーネル関数スペクトルを記憶しておく。すなわち、０．１から１．０までの全て分散値に対してカーネル関数スペクトルをあらかじめ計算しておき、それらをカーネル関数スペクトル記憶部１３１４に記憶しておく。この場合、１０種類のカーネル関数スペクトルが記憶される。本実施形態では、分散が０．１、０．２、・・・１．０のカーネル関数スペクトルをＧ_０．１（ｗ）、Ｇ_０．２（ｗ）、・・・Ｇ_１．０（ｗ）と表記する。

第５実施形態で説明した通り、ピークが明確な統合相関関数を生成する上では、様々な幅を有するカーネル関数を用意することが望ましい。したがって、多くの分散値に対応するカーネル関数スペクトルを用意する、例えば、０．００１から０．００１刻みで１０，０００種類の分散（０．００１,０．００２,０．００３,・・・）に対応するカーネル関数スペクトルを記憶することが望ましい。しかし、記憶容量が増大する問題が生じる。したがって、記憶容量の増大を抑えつつ、ピークが明確な統合相関関数を生成するためには、０．０５から０．０５刻み、もしくは０．１刻みで約５〜２０種類の分散に対応するカーネル関数スペクトルを用意する必要がある。なお、刻み幅は必ずしも均等である必要はない。統合相関関数のピークを強調する上では、０に近い分散値では細かく、０から遠ざかるにつれて刻み幅を広くしたほうがよい。また、分散値がある閾値を超えた場合には、十分に大きな値を有する１種類の分散値に対応するカーネル関数スペクトルを用意すればよい。

カーネル関数スペクトル選択部１３１２は、分散計算部２０３から供給された分散に基づいて、カーネル関数スペクトル記憶部１３１４に記憶されている複数のカーネル関数スペクトルの中から、一つのカーネル関数スペクトルを選択し、乗算部１０１３へ伝達する。カーネル関数スペクトル選択部１３１２は、分散計算部２０３から供給された分散値に最も近い分散値で生成されたカーネル関数スペクトルを選択する。例えば、供給された分散値が０から０．１５未満の値であれば、カーネル関数スペクトルＧ_０．１（ｗ）を選択する。同様に、供給された分散値が０．１５以上から０．２５未満の値であれば、カーネル関数スペクトルＧ_０．２（ｗ）を選択する。

＜本実施形態の効果＞
本実施形態によれば、あらかじめカーネル関数を計算しておくので、方向推定処理中にはカーネル関数の計算が不要になる。したがって、第５実施形態よりも少ない計算量で方向推定を行うことができる。

［第７実施形態］
次に本発明の第７実施形態に係る波源方向推定装置について、図１４を用いて説明する。本実施形態に係る波源方向推定装置は、上記第５実施形態の波源方向推定装置と比べると、統合相関関数計算部９４２に代えて、統合相関関数計算部１４４２を有する点で異なる。図１４は、本実施形態に係る統合相関関数計算部１４４２の構成を説明するためのブロック図である。本実施形態に係る統合相関関数計算部１４４２は、上記第５実施形態の統合相関関数計算部９４２と比べると、逆変換部１１２１_１，１１２１_２,・・・，１１２１_Ｋと統合部１１２２とに代えて、統合部１４２１と逆変換部１４２２とを有する点で異なる。その他の構成および動作は、第５実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。

図１４は、統合相関関数計算部１４４２のブロック図である。統合相関関数計算部１４４２は、統合部１４２１と逆変換部１４２２とを備える。

統合部１４２１は、周波数別クロススペクトル生成部９４１_１，９４１_２,・・・，９４１_Ｋから供給される周波数別クロススペクトルを統合し、統合クロススペクトルとして逆変換部１４２２へ伝達する。個別に求めた複数の周波数別クロススペクトルを混合したり、重ね合わせたりすることにより、一つの統合クロススペクトルを求める。統合には、第５実施形態の統合部１１２２と同様に総和や総乗が用いられる。統合に総和を用いた場合、周波数別クロススペクトル生成部９４１_ｋから供給された周波数別クロススペクトルをＵＭ_ｋ（０,ｎ）,ＵＭ_ｋ（１,ｎ）,・・・,ＵＭ_ｋ（Ｗ・１,ｎ）とすると、統合クロススペクトルＵ（ｋ,ｎ）は次のように計算される。

また、総乗を用いた場合、統合クロススペクトルＵ（ｋ,ｎ）は次のように計算される。

第５実施形態の統合部１１２２と同様に、目標音源の存在する周波数や、目的音源のパワーが大きい周波数があらかじめ判明している場合には、統合クロススペクトルＵ（ｋ,ｎ）を生成する際に補正してもよい。この場合、第５実施形態と同様に、周波数の選択や重み付けという形で影響度を制御する。例えば目標音の存在する周波数の集合をΩとすると、帯域の選択により統合クロススペクトルＵ（ｋ,ｎ）を求める場合は、次のように計算する。

また、重み付けを用いる場合には、Ｕ（ｋ,ｎ）は次のように計算される。

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

逆変換部１４２２は、統合部１４２１から供給される統合クロススペクトルの逆変換を行い、統合相関関数として推定方向情報生成部２０５に伝達する。本実施形態でも、逆変換には逆フーリエ変換を用いる方法について説明する。統合部１４２１から供給された統合クロススペクトルをＵ（ｋ,ｎ）とすると、Ｕ（ｋ,ｎ）の逆変換により得られる統合相関関数ｕ（τ,ｎ）は次のように計算される。

＜本実施形態の効果＞
本実施形態によれば、周波数別クロススペクトルを統合してから逆変換を行い、統合相関関数を得る。このため、周波数別クロススペクトルごとに逆変換を行っていた第５実施形態と比べると、逆変換の回数が少なくなる。したがって、第５実施形態よりも少ない計算量で統合相関関数を求めることができる。

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

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

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

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

表示部１５０２は、波源方向推定装置２００から供給された推定方向情報を画像などの可視化データに変換し、ディスプレイなどの表示装置に表示する。基本的な可視化方法は、ある時刻における相関関数を２次元グラフで表示する方法である。その際、横軸に方向、縦軸に相関値を表示する。ある時刻だけでなく、相関関数の時間変化を３次元で表示する方法も有効である。時間変化を表示することにより、目標音源の出現の明確化、目標音源の移動パタン、目標音源の移動方向の予測などが可能になる。３次元ではなく、２次元平面に投影する方法も有効である。３次元だと表示されたときに裏側が見づらいという問題がある。上から投影した平面上に表示すれば、死角が無くなり一覧性が向上する。相関値は、色の濃淡ではなく、等高線で表現してもよい。

図１５Ｂは、本実施形態に係る波源方向推定システム１５００の表示部１５０２で表示した画像１５１０の一例を示す図であり、波源方向推定装置２００から供給された推定方向情報から得られたものである。これは、本実施形態の効果を確認する目的で取得した。例の作成には、屋外環境において、方位４０度にドローンの飛翔音が発生する状況の音を用いた。数センチ間隔で設置した二つのマイクを用いて集音した。

図１５Ｂは、色が黒いほど相関値が高いことを表している。方位角の範囲は０から１８０度である。横軸は時刻を表している。図１５Ｂを参照すると、方位４０度の相関値が高くなっていることが分かる。このことから、ドローンの到来方向が約４０度であることが分かる。

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

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

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

［実施形態の他の表現］
上記の実施形態の一部又は全部は、以下の付記のようにも記載されうるが、以下には限られない。
（付記１）
波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換手段と、
複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
前記クロススペクトルの分散を計算する分散計算手段と、
前記クロススペクトルと前記分散とに基づいて、相関関数を計算する相関関数計算手段と、
を備える相関関数生成装置。
（付記２）
前記分散計算手段は、前記クロススペクトルの位相に基づいて分散を計算する付記１に記載の相関関数生成装置。
（付記３）
前記分散計算手段は、前記クロススペクトルの円周分散を計算する付記１に記載の相関関数生成装置。
（付記４）
前記分散計算手段は、前記入力信号取得手段で過去に取得された入力信号から変換された周波数領域信号に基づいて計算されたクロススペクトルに基づいて分散を計算する付記１に記載の相関関数生成装置。
（付記５）
前記相関関数計算手段は、
前記分散に基づいて前記クロススペクトルを重み付けする重み付け手段を備え、
重み付されたクロススペクトルに基づいて、相関関数を計算する、付記１乃至４のいずれか１項に記載の相関関数生成装置。
（付記６）
前記クロススペクトルを平均して平均クロススペクトルを計算する平均計算手段をさらに備え、
前記分散計算手段は、前記平均クロススペクトルから分散を計算する付記１乃至５のいずれか１項に記載の相関関数生成装置。
（付記７）
前記相関関数計算手段は、
前記分散と前記クロススペクトルとに基づいて、複数の第１周波数別クロススペクトルを生成する周波数別クロススペクトル生成手段と、
複数の前記第１周波数別クロススペクトルを統合して一つの相関関数を計算する統合相関関数計算手段と、
を備える付記１乃至６のいずれか１項に記載の相関関数生成装置。
（付記８）
前記統合相関関数計算手段は、
複数の前記第１周波数別クロススペクトルの逆変換を計算し、複数の周波数別相関関数を得る第１逆変換手段と、
前記第１逆変換手段により計算された複数の前記周波数別相関関数を統合して一つの相関関数を計算する第１統合手段と、
を備える付記７に記載の相関関数生成装置。
（付記９）
前記周波数別クロススペクトル生成手段は、
前記分散に基づいて、カーネル関数スペクトルを生成するカーネル関数スペクトル生成手段と、
複数の前記第１周波数別クロススペクトルに、前記カーネル関数スペクトルを乗じた複数の第２周波数別クロススペクトルを得る乗算手段と、
を備える付記７または８に記載の相関関数生成装置。
（付記１０）
前記周波数別クロススペクトル生成手段は、
複数のカーネル関数スペクトルを記憶するカーネル関数スペクトル記憶手段と、
前記分散に基づいて、前記カーネル関数スペクトル記憶手段に記憶された複数の前記カーネル関数スペクトルから一つのカーネル関数スペクトルを選択するカーネル関数スペクトル選択手段と、
をさらに備え、
前記乗算手段は、複数の前記第１周波数別クロススペクトルに、選択された前記一つのカーネル関数スペクトルを乗じた複数の第２周波数別クロススペクトルを得る、付記９に記載の相関関数生成装置。
（付記１１）
前記統合相関関数計算手段は、
複数の前記第１周波数別クロススペクトルを統合して、一つの統合クロススペクトルを得る第２統合手段と、
前記一つの統合クロススペクトルの逆変換を計算して一つの相関関数を得る第２逆変換手段と、
を備える付記７乃至１０のいずれか１項に記載の相関関数生成装置。
（付記１２）
波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
前記クロススペクトルの分散を計算する分散計算ステップと、
前記クロススペクトルと前記分散とに基づいて、相関関数を計算して生成する相関関数計算ステップと、
を含む相関関数生成方法。
（付記１３）
波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
前記クロススペクトルの分散を計算する分散計算ステップと、
前記クロススペクトルと前記分散とに基づいて、相関関数を計算して生成する相関関数計算ステップと、
をコンピュータに実行させる相関関数生成プログラム。
（付記１４）
付記１乃至１１のいずれか１項に記載の相関関数生成装置を有し、
前記相関関数生成装置により生成された相関関数に基づいて、前記波源の方向を推定する波源方向推定装置。

Claims (10)

1. 波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
   複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換手段と、
   複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
   前記クロススペクトルの位相の分散として円周分散を計算する分散計算手段と、
   前記円周分散に基づいて算出された重みによって重み付された前記クロススペクトルに基づいて、相関関数を計算する相関関数計算手段と、
   を備える相関関数生成装置。
2. 前記クロススペクトルを平均して平均クロススペクトルを計算する平均計算手段をさらに備え、
   前記分散計算手段は、前記平均クロススペクトルから前記円周分散を計算する請求項１に記載の相関関数生成装置。
3. 前記相関関数計算手段は、
   前記円周分散と前記クロススペクトルとに基づいて、複数の周波数別クロススペクトルを生成する複数の周波数別クロススペクトル生成手段と、
   前記複数の周波数別クロススペクトルを統合して一つの相関関数を計算する統合相関関数計算手段と、
   を備える請求項１または２に記載の相関関数生成装置。
4. 前記周波数別クロススペクトル生成手段は、
   前記円周分散に基づいて、カーネル関数スペクトルを生成するカーネル関数スペクトル生成手段と、
   周波数別の前記クロススペクトルに、前記カーネル関数スペクトルを乗じて前記周波数別クロススペクトルを得る乗算手段と、
   を備える請求項３に記載の相関関数生成装置。
5. 前記周波数別クロススペクトル生成手段は、
   複数のカーネル関数スペクトルを記憶するカーネル関数スペクトル記憶手段と、
   前記円周分散に基づいて、前記カーネル関数スペクトル記憶手段に記憶された複数の前記カーネル関数スペクトルから一つのカーネル関数スペクトルを選択するカーネル関数スペクトル選択手段と、
   周波数別の前記クロススペクトルに、選択された前記一つのカーネル関数スペクトルを乗じて前記周波数別クロススペクトルを得る乗算手段と、
   を備える請求項３に記載の相関関数生成装置。
6. 前記統合相関関数計算手段は、
   前記複数の周波数別クロススペクトルの逆変換を計算し、複数の周波数別相関関数を得る複数の第１逆変換手段と、
   前記複数の第１逆変換手段により計算された複数の前記周波数別相関関数を統合して一つの相関関数を計算する第１統合手段と、
   を備える請求項３乃至５のいずれか１項に記載の相関関数生成装置。
7. 前記統合相関関数計算手段は、
   前記複数の周波数別クロススペクトルを統合して、一つの統合クロススペクトルを得る第２統合手段と、
   前記一つの統合クロススペクトルの逆変換を計算して一つの相関関数を得る第２逆変換手段と、
   を備える請求項３乃至５のいずれか１項に記載の相関関数生成装置。
8. 波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
   複数の前記入力信号取得ステップで取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
   複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
   前記クロススペクトルの位相の分散として円周分散を計算する分散計算ステップと、
   前記円周分散に基づいて算出された重みによって重み付された前記クロススペクトルに基づいて、相関関数を計算して生成する相関関数計算ステップと、
   を含む相関関数生成方法。
9. 波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
   複数の前記入力信号取得ステップで取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
   複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
   前記クロススペクトルの位相の分散として円周分散を計算する分散計算ステップと、
   前記円周分散に基づいて算出された重みによって重み付された前記クロススペクトルに基づいて、相関関数を計算して生成する相関関数計算ステップと、
   をコンピュータに実行させる相関関数生成プログラム。
10. 請求項１乃至７のいずれか１項に記載の相関関数生成装置を有し、
    前記相関関数生成装置により生成された相関関数に基づいて、前記波源の方向を推定する波源方向推定装置。