WO2020003343A1

WO2020003343A1 - 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体

Info

Publication number: WO2020003343A1
Application number: PCT/JP2018/023971
Authority: WO
Inventors: 友督荒井; 玲史近藤; 裕三仙田
Original assignee: 日本電気株式会社
Priority date: 2018-06-25
Filing date: 2018-06-25
Publication date: 2020-01-02
Also published as: US20210263126A1; US11408963B2; JPWO2020003343A1; JP6933303B2

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（ｗ，ｎ）｜は以下の式１８で与えられる。ただし、式１８において、ｐは１以上Ｐ以下の整数である（Ｐは整数）。

　波源方向推定を行うときに重要な情報は位相成分であるため、式１８においては、振幅成分に適当な定数を用いている。なお、ｋの定数倍の周波数の振幅成分｜Ｕ_k（ｗ，ｎ）｜として、１．０の代わりに周波数ｋの平均クロススペクトルの絶対値｜ＳＳ₁₂（ｋ，ｎ）｜を用いてもよい。すなわち、周波数別基本クロススペクトル計算部５５１は、以下の式１９を用いて振幅成分｜Ｕ_k（ｗ，ｎ）｜を求めてもよい。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

　〔周波数寄与率計算部〕
　周波数寄与率計算部１５８には、周波数別推定方向情報計算部１５７から周波数別推定方向情報が入力される。周波数寄与率計算部１５８は、周波数別推定方向情報計算部１５７から供給された周波数別推定方向情報に基づいて、音源方向ごとの推定方向情報に対する周波数別の寄与率（周波数寄与率とも呼ぶ）を計算する。周波数寄与率計算部１５８は、算出した周波数寄与率を重み計算部１６に出力する。

　例えば、周波数寄与率計算部１５８は、以下の式３３を用いて、音源方向θがθ’のときの周波数ｋ’の周波数寄与率ｙ_θ’,k’,r（ｎ）を計算する。

　全ての音源方向θおよび全ての周波数ｋについて、式３３で表記される周波数寄与率ｙ_θ’,k’,r（ｎ）を求めると、周波数寄与率はｙ_θ,r（ｋ，ｎ）と書き換えられる。

　〔周波数統合部〕
　周波数統合部１５９には、周波数別推定方向情報計算部１５７から周波数別推定方向情報が入力される。周波数統合部１５９は、周波数別推定方向情報計算部１５７から供給された周波数別推定方向情報を統合する。すなわち、周波数統合部１５９は、周波数（周波数ビン番号）ごとに別々に算出された周波数別推定方向情報を一つに統合することによって、統合された推定方向情報を生成する。例えば、周波数統合部１５９は、個別に求めた複数の周波数別推定方向情報を合成したり、重ね合わせたりすることにより、一つの推定方向情報に統合する。周波数統合部１５９は、統合された推定方向情報を統合部１７に出力する。

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

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

　〔重み計算部〕
　重み計算部１６には、推定方向情報・周波数寄与率生成部１５－１～Ｒのそれぞれの周波数寄与率計算部１５８から周波数寄与率が入力される。重み計算部１６は、推定方向情報・周波数寄与率生成部１５－１～Ｒから供給される周波数寄与率に基づいて、マイクペア間の周波数寄与率（単に寄与率とも呼ぶ）の類似度に応じた音源方向ごとの重みを計算する。重み計算部１６は、算出した音源方向ごとの重みを統合部１７に出力する。

　ここで、３つのマイクによって方向推定する例における重み算出方法について図面を参照しながら説明する。図６は、推定方向情報・周波数寄与率生成部１５が、マイクＡ、マイクＢ、およびマイクＣに関して算出する全てのマイクペア（ＡＢ、ＢＣ、ＡＣ）の推定方向情報の一例である。図７は、音源方向θがθ’の場合のマイクペアごとの周波数寄与率を示す。重み計算部１６は、異なるマイクペアの寄与率が類似している方向θほど大きくなるような重みを計算する。

　例えば、重み計算部１６は、マイクペアＡＢの周波数寄与率ｙ_θ,AB（ｋ，ｎ）およびマイクペアＢＣの周波数寄与率ｙ_θ,BC（ｋ，ｎ）を以下の式３６に適用して２つのマイクペア間の類似度ｖ_θ，AB-BC（ｎ）を計算する。

　上記の式３６を一般化すると、重み計算部１６は、マイクペアｒ_iの周波数寄与率ｙ_θ,ｒi（ｋ，ｎ）とマイクペアｒ_jの周波数寄与率ｙ_θ,ｒｊ（ｋ，ｎ）との類似度ｖ_θ，ｒ'（ｎ）を以下の式３７に適用して算出できる。

　なお、式３７において、ｒ’は、マイクペアｒ_iとマイクペアｒ_jとの組み合わせ番号を表す。マイクペアｒ_iの周波数寄与率ｙ_θ,ｒi（ｋ，ｎ）と、マイクペアｒ_ｊの周波数寄与率ｙ_θ,ｒｊ（ｋ，ｎ）とが一致する場合、式３７で表される類似度ｖ_θ，ｒ'（ｎ）は最大値１となる。

　例えば、重み計算部１６は、以下の式３８を用いて、類似度ｖ_θ，ｒ'（ｎ）の相乗平均から音源方向θごとの重みa_θを計算する。

　また、例えば、重み計算部１６は、以下の式３９を用いて、類似度ｖ_θ，ｒ'（ｎ）の相加平均から音源方向θごとの重みa_θ（ｎ）を計算できる。

　上述のように計算される重みは、全てのマイクペアの類似度が高いほど大きくなる。

　〔統合部〕
　統合部１７には、推定方向情報・周波数寄与率生成部１５－１～Ｒから推定方向情報が入力され、重み計算部１６から重みが入力される。統合部１７は、推定方向情報・周波数寄与率生成部１５－１～Ｒから供給されるマイクペアごとの推定方向情報を全てのマイクペアについて統合し、統合した推定方向情報に重み計算部１６から供給される重みをつけた統合推定方向情報を生成する。統合部１７は、マイクペアごとの推定方向情報の総和や総乗を取ることで統合し、統合した推定方向情報に重みを掛けることによって統合推定方向情報を生成する。統合部１７は、生成した統合推定方向情報を出力する。例えば、統合部１７は、上位システム（図示しない）に統合推定方向情報を出力する。

　例えば、統合部１７は、以下の式４０を用いて、重みa_θ（ｎ）と推定方向情報Ｈ_r（θ，ｎ）の総乗との積を統合推定方向情報Ｈ（θ，ｎ）として算出する。

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

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

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

　（動作）
　次に、本実施形態の波源方向推定装置１０の動作の概略について図面を参照しながら説明する。図９は、波源方向推定装置１０の動作の概略について説明するためのフローチャートである。なお、図９のフローチャートに沿った説明においては、波源方向推定装置１０を動作の主体として説明する。

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

　次に、波源方向推定装置１０は、複数のマイクのそれぞれに対応する入力信号から２つの入力信号を選択する（ステップＳ１２）。すなわち、波源方向推定装置１０は、複数のマイクの中から２つのマイクを選択してマイクペアとする。ここでは、波源方向推定装置１０は、複数のマイクペアを選択する。

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

　次に、波源方向推定装置１０は、選択されたマイクペアの入力信号と相対遅延時間とを用いて、そのマイクペアの入力信号間の推定方向情報と周波数寄与率とを生成する（ステップＳ１４）。すなわち、波源方向推定装置１０は、推定方向情報と周波数寄与率とをマイクペアごとに生成する。

　次に、波源方向推定装置１０は、マイクペア間の周波数寄与率の類似度に応じた音源方向ごとの重みを算出する（ステップＳ１５）。

　次に、波源方向推定装置１０は、マイクペアごとの推定方向情報に重みを付けて統合することで統合推定方向情報を算出する（ステップＳ１６）。

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

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

　以上のように、本実施形態の波源方向推定装置は、複数の入力部、信号選択部、相対遅延時間計算部、複数の推定方向情報・周波数寄与率生成部、重み計算部、および統合部を備える。

　入力部は、複数のセンサによって検出される波動に基づく電気信号を入力信号として取得する。

　信号選択部は、複数の入力信号のうち二つを組み合わせたペアを複数組選択する。

　相対遅延時間計算部は、ペアをなす二つの入力信号の供給元のセンサへの波源方向ごとの波動の到達時間差を相対遅延時間として算出する。

　推定方向情報・周波数寄与率生成部は、信号選択部によって選択された全てのペアのうち１組のペアをなす二つの入力信号を取得するとともに、取得したペアの相対遅延時間を取得する。推定方向情報・周波数寄与率生成部は、ペアをなす二つの入力信号間の周波数別相互相関関数とペアの相対遅延時間とを用いて、ペアをなす二つの入力信号間の推定方向情報と周波数別の寄与率とを生成する。

　重み計算部は、複数の推定方向情報・周波数寄与率生成部から周波数ごとの寄与率を取得し、異なるペア間の寄与率の類似度に応じて波源方向ごとの重みを計算する。

　統合部は、複数の推定方向情報・周波数寄与率生成部から推定方向情報を取得するとともに、重み計算部から重みを取得し、複数のペアごとの推定方向情報に重みをつけて統合して統合推定方向情報を生成する。

　例えば、推定方向情報・周波数寄与率生成部は、変換部、クロススペクトル計算部、平均計算部、分散計算部、周波数別クロススペクトル生成部、逆変換部、周波数別推定方向情報計算部、周波数統合部、および周波数寄与率計算部を備える。変換部は、ペアをなす二つの入力信号を周波数領域の変換信号に変換する。クロススペクトル計算部は、変換部によって変換された変換信号を用いてクロススペクトルを計算する。平均計算部は、クロススペクトル計算部によって算出されるクロススペクトルを用いて平均クロススペクトルを計算する。分散計算部は、平均計算部によって算出される平均クロススペクトルを用いて分散を計算する。周波数別クロススペクトル生成部は、平均計算部によって算出される平均クロススペクトルと、分散計算部によって算出される分散とを用いて周波数別クロススペクトルを計算する。逆変換部は、周波数別クロススペクトル生成部によって算出される周波数クロススペクトルを逆変換して周波数別相互相関関数を計算する。周波数別推定方向情報計算部は、ペアの周波数別相互相関関数と相対遅延時間とに基づいて波源方向と相関値との対応関係である周波数別推定方向情報を計算する。周波数統合部は、周波数別推定方向情報を統合して推定方向情報を計算する。周波数寄与率計算部は、周波数別推定方向情報に基づいて波源方向ごとの推定方向情報に対する周波数別の寄与率を計算する。

　例えば、重み計算部は、複数のペア間において寄与率の類似度を計算し、算出した類似度に基づいて波源方向ごとの重みを計算する。

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

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

　図１０は、本実施形態の波源方向推定装置２０の構成の一例を示すブロック図である。図１０のように、波源方向推定装置２０は、複数の入力端子２１、信号選択部２２、相対遅延時間計算部２３、複数の周波数別推定方向情報生成部２５、重み生成部２６、および統合部２７を備える。以下において、第１の実施形態の波源方向推定装置１０と同様の構成・動作については説明を省略する場合がある。

　〔入力端子〕
　複数の入力端子２１－１～ｐ（入力部とも呼ぶ）のそれぞれは、図示しないマイクに接続される（ｐは２以上の整数）。複数の入力端子２１－１～ｐのそれぞれには、異なる位置に配置されたマイクによって集音された音波（音信号とも呼ぶ）から変換された電気信号が入力信号として入力される。なお、入力端子２１－１～ｐは、第１の実施形態の入力端子２１－１～ｐと同様であるため、詳細な説明は省略する。

　〔信号選択部〕
　信号選択部２２は、入力端子２１－１～ｐに供給されるｐ個の入力信号のうちの２つの入力信号を選択する。信号選択部２２は、選択した２つの入力信号を周波数別推定方向情報生成部２５－１～Ｒに出力し、それらの入力信号の供給元であるマイクの位置情報（以下、マイク位置情報とも呼ぶ）を相対遅延時間計算部２３に出力する（Ｒは１以上の整数）。なお、信号選択部２２は、第１の実施形態の信号選択部１２と同様であるため、詳細な説明は省略する。

　〔相対遅延時間計算部〕
　相対遅延時間計算部２３には、信号選択部２２からマイク位置情報が入力される。相対遅延時間計算部２３は、予め設定された音源探索対象方向（音源方向とも呼ぶ）と、マイク位置情報とを用いて、信号選択部２２が選択した全てのマイクペアごとに相対遅延時間を計算する。相対遅延時間計算部２３は、設定された音源方向と相対遅延時間とをセットにして周波数別推定方向情報生成部２５に出力する。なお、相対遅延時間計算部２３は、第１の実施形態の相対遅延時間計算部１３と同様であるため、詳細な説明は省略する。

　〔周波数別推定方向情報生成部〕
　周波数別推定方向情報生成部２５－１～Ｒのそれぞれには、信号選択部２２によって選択された全マイクペアのうち１組のマイクペアの入力信号が入力される。また、周波数別推定方向情報生成部２５－１～Ｒのそれぞれには、入力信号の供給元のマイクペアの相対遅延時間が相対遅延時間計算部２３から入力される。周波数別推定方向情報生成部２５－１～Ｒのそれぞれは、１組のマイクペアの入力信号と相対遅延時間とを用いて、そのマイクペアごとに周波数別推定方向情報を生成する。周波数別推定方向情報生成部２５－１～Ｒは、生成したマイクペアごとの周波数別推定方向情報を重み生成部２６および統合部２７に出力する。

　ここで、図１１を参照しながら、周波数別推定方向情報生成部２５の詳細な構成について説明する。図１１は、周波数別推定方向情報生成部２５の構成の一例を示すブロック図である。周波数別推定方向情報生成部２５は、変換部２５１、クロススペクトル計算部２５２、平均計算部２５３、分散計算部２５４、周波数別クロススペクトル生成部２５５、逆変換部２５６、および周波数別推定方向情報計算部２５７を備える。すなわち、周波数別推定方向情報生成部２５は、第１の実施形態の推定方向情報・周波数寄与率生成部１５から周波数統合部１５９および周波数寄与率計算部１５８を省いた構成を有する。なお、変換部２５１、クロススペクトル計算部２５２、平均計算部２５３、分散計算部２５４、周波数別クロススペクトル生成部２５５、逆変換部２５６、および周波数別推定方向情報計算部２５７は、第１の実施形態の対応する構成と同様である。そのため、変換部２５１、クロススペクトル計算部２５２、平均計算部２５３、分散計算部２５４、周波数別クロススペクトル生成部２５５、逆変換部２５６、および周波数別推定方向情報計算部２５７についての詳細な説明は省略する。

　〔重み生成部〕
　重み生成部２６には、周波数別推定方向情報生成部２５－１～Ｒのそれぞれからマイクペアごとの周波数別推定方向情報が入力される。重み生成部２６は、周波数別推定方向情報生成部２５－１～Ｒのそれぞれから供給されるマイクペアごとの周波数別推定方向情報を用いて、推定方向情報がピークをもつ音源方向ごとの重みを計算する。重み生成部２６は、算出した音源方向ごとの重みを統合部２７に出力する。

　ここで、図１２を参照しながら、重み生成部２６の詳細な構成について説明する。図１２は、重み生成部２６の構成の一例を示すブロック図である。重み生成部２６は、複数の周波数統合部２６１－１～Ｒ、複数の波源方向候補選出部２６２－１～Ｒ、波源方向候補統合部２６３、複数の波源方向候補別周波数寄与率計算部２６４－１～Ｒ、および重み計算部２６５を有する。周波数統合部２６１－１～Ｒ、波源方向候補選出部２６２－１～Ｒ、および波源方向候補別周波数寄与率計算部２６４－１～Ｒのそれぞれは、いずれかのマイクペアに対応する。そのため、周波数統合部２６１－１～Ｒ、波源方向候補選出部２６２－１～Ｒ、および波源方向候補別周波数寄与率計算部２６４－１～Ｒのそれぞれには、マイクペアの番号（１～Ｒ）を付す。

　〔周波数統合部〕
　周波数統合部２６１－１～Ｒのそれぞれには、周波数別推定方向情報生成部２５－１～Ｒのそれぞれから周波数別推定方向情報が入力される。周波数統合部２６１－１～Ｒのそれぞれは、周波数別推定方向情報生成部２５－１～Ｒのそれぞれから供給される周波数別推定方向情報を統合し、推定方向情報を計算する。なお、周波数統合部２６１の動作は、第１の実施形態の周波数統合部１５９と同様であるため、詳細な説明は省略する。周波数統合部２６１－１～Ｒのそれぞれは、算出した推定方向情報を波源方向候補選出部２６２－１～Ｒのそれぞれに出力する。

　〔波源方向候補選出部〕
　波源方向候補選出部２６２－１～Ｒのそれぞれには、周波数統合部２６１－１～Ｒのそれぞれから推定方向情報が入力される。波源方向候補選出部２６２－１～Ｒのそれぞれは、周波数統合部２６１－１～Ｒのそれぞれから供給される推定方向情報からピークを検出し、ピークの方向を波源方向候補として選出する。波源方向候補選出部２６２－１～Ｒのそれぞれは、選出した波源方向候補を波源方向候補統合部２６３に出力する。

　例えば、波源方向候補選出部２６２は、推定方向情報が極大値を取る音源方向θを波源方向候補として選出する。また、例えば、波源方向候補選出部２６２は、音源数が３つと事前に仮定できる場合には、推定方向情報が大きい順に３つのピークを選別し、それら３つの値を取る音源方向θを波源方向候補として選出してもよい。特に、波源方向候補選出部２６２は、音源数が１つと事前に仮定できる場合には、推定方向情報が最大値となる音源方向θを波源方向候補として選出する。この場合、周波数統合部２６１から供給される推定方向情報をＨ_r（θ，ｎ）とすると、波源方向候補選出部２６２は、以下の式４１を用いて、推定方向情報Ｈ_r（θ，ｎ）が最大値を取るような推定方向情報Ｈ_r（θ，ｎ）の引数を要素とする集合を求める。

　〔波源方向候補統合部〕
　波源方向候補統合部２６３には、波源方向候補選出部２６２－１～Ｒのそれぞれから波源方向候補が入力される。波源方向候補統合部２６３は、波源方向候補選出部２６２－１～Ｒのそれぞれによって選出された波源方向候補を統合する。波源方向候補統合部２６３は、統合した波源方向候補を波源方向候補別周波数寄与率計算部２６４－１～Ｒのそれぞれに出力する。

　ここで、波源方向候補選出部２６２－１～Ｒのそれぞれによって選出されたマイクペアごとの複数の波源方向候補をベクトルφ_iとする。波源方向候補統合部２６３は、以下の式４２を用いて波源方向候補を統合し、統合した波源方向候補Φ'_i(n)を計算する。

　ここで、Φ'_i(n)の成分で重複する方向がある場合、波源方向候補統合部２６３は、重複する方向を省いて、以下の式４３を用いてマイクペアごとの波源方向候補を全てのマイクペアについて統合し、統合した波源方向候補Φ'_i(n)を計算する。ただし、式４３において、Ｉは全ての波源方向候補の数を示す。

　〔波源方向候補別周波数寄与率計算部〕
　波源方向候補別周波数寄与率計算部２６４－１～Ｒのそれぞれには、周波数別推定方向情報生成部２５－１～Ｒのそれぞれからマイクペアごとの周波数別推定方向情報が入力される。また、波源方向候補別周波数寄与率計算部２６４－１～Ｒのそれぞれには、統合された波源方向候補が波源方向候補統合部２６３から入力される。波源方向候補別周波数寄与率計算部２６４－１～Ｒのそれぞれは、マイクペアごとの周波数別推定方向情報と、統合された波源方向候補とを用いて、波源方向候補ごとの推定方向情報に対する周波数寄与率である波源方向候補別周波数寄与率を計算する。波源方向候補別周波数寄与率計算部２６４－１～Ｒのそれぞれは、算出した波源方向候補別周波数寄与率を重み計算部２６５に出力する。第１の実施形態の周波数寄与率計算部１５８は、全ての音源方向について周波数寄与率を求める。それに対し、本実施形態の波源方向候補別周波数寄与率計算部２６４は、選出された波源方向候補について周波数寄与率を求める。

　例えば、波源方向候補別周波数寄与率計算部２６４は、ある波源方向候補φ_iの周波数ｋの波源方向候補別周波数寄与率ｙ_φi,r（ｋ，ｎ）を以下の式４４を用いて計算する。ただし、式４４において、Ｈ_k,r（Φ_i，ｎ）は、周波数別推定方向情報を示す。

　〔重み計算部〕
　重み計算部２６５には、波源方向候補別周波数寄与率計算部２６４－１～Ｒのそれぞれから波源方向候補別周波数寄与率が入力される。重み計算部２６５は、波源方向候補別周波数寄与率計算部２６４－１～Ｒのそれぞれから供給される波源方向候補別周波数寄与率を用いて、マイクペア間の周波数寄与率の類似度に応じた波源方向候補ごとの重みを計算する。重み計算部２６５は、算出した波源方向候補ごとの重みを統合部２７に出力する。

　第１の実施形態の重み計算部１６は、全ての音源方向について重みを求める。それに対し、本実施形態の重み計算部２６５は、波源方向候補について重みを求める。重み計算部２６５が算出する波源方向候補ごとの重みをａ_Φ（ｎ）と表記する。重み計算部２６５のその他の動作は重み計算部１６と同様であるため、詳細な説明は省略する。

　〔統合部〕
　統合部２７には、周波数別推定方向情報生成部２５－１～Ｒのそれぞれから周波数別推定方向情報が入力され、重み生成部２６から波源方向候補ごとの重みが入力される。統合部２７は、周波数別推定方向情報生成部２５－１～Ｒから供給される周波数別推定方向情報を波源方向候補ごとに全てのマイクペアについて統合し、統合した波源方向候補ごとの周波数別推定方向情報に重みを付け、波源方向候補ごとの統合推定方向情報を算出する。統合部２７は、マイクペアごとの周波数別推定方向情報の総和や相乗を取ることで統合し、統合した推定方向情報に重みを掛けることによって波源方向候補ごとの統合推定方向情報を生成する。例えば、統合部２７は、以下の式４５を用いて、波源方向候補Φごとの統合推定方向情報Ｈ（Φ，ｎ）を算出する。

　ただし、式４５において、統合推定方向情報Ｈ（Φ，ｎ）は、波源方向候補Φごとに周波数別推定方向情報を統合した推定方向情報を表す。統合部２７は、周波数統合部２６１と同様に統合推定方向情報Ｈ（Φ，ｎ）を計算する。ただし、周波数統合部２６１が全ての音源方向について推定方向情報を求めるのに対し、統合部２７は波源方向候補Φについて推定方向情報を求める。

　（動作）
　次に、本実施形態の波源方向推定装置２０の動作の概略について図面を参照しながら説明する。図１３は、波源方向推定装置２０の動作の概略について説明するためのフローチャートである。なお、図１３のフローチャートに沿った説明においては、波源方向推定装置２０を動作の主体として説明する。

　図１３において、まず、波源方向推定装置２０は、複数のマイクから入力信号を取得する（ステップＳ２１）。

　次に、波源方向推定装置２０は、複数のマイクのそれぞれに対応する入力信号から２つの入力信号を選択する（ステップＳ２２）。すなわち、波源方向推定装置２０は、複数のマイクの中から２つのマイクを選択してマイクペアとする。ここでは、波源方向推定装置２０は、複数のマイクペアを選択する。

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

　次に、波源方向推定装置２０は、選択されたマイクペアの入力信号と相対遅延時間とを用いて、そのマイクペアの周波数別推定方向情報を生成する（ステップＳ２４）。すなわち、波源方向推定装置２０は、周波数別推定方向情報をマイクペアごとに生成する。

　次に、波源方向推定装置２０は、マイクペアごとの周波数別推定方向情報を統合した推定方向情報に基づいて波源方向候補を選出し、波源方向候補別に周波数寄与率を計算する（ステップＳ２５）。

　次に、波源方向推定装置２０は、波源方向候補別周波数寄与率を用いて、マイクペア間の周波数寄与率の類似度に応じた波源方向候補ごとの重みを計算する（ステップＳ２６）。

　次に、波源方向推定装置２０は、波源方向候補ごとの重みに基づいて、波源方向候補ごとに周波数別推定方向情報を統合して統合推定方向情報を算出する（ステップＳ２７）。

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

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

　以上のように、本実施形態の波源方向推定装置は、複数の入力部、信号選択部、相対遅延時間計算部、複数の周波数別推定方向情報生成部、重み生成部、および統合部を備える。

　周波数別推定方向情報生成部は、信号選択部によって選択された全てのペアのうち１組のペアをなす二つの入力信号を取得するとともに、取得したペアの相対遅延時間を取得する。周波数別推定方向情報生成部は、ペアをなす二つの入力信号を用いて算出される周波数別相互相関関数とペアの相対遅延時間とを用いて、ペアごとの周波数別推定方向情報を計算する。

　重み生成部は、複数の周波数別推定方向情報生成部からペアごとの周波数別推定方向情報を取得する。重み生成部は、取得したペアごとの周波数別推定方向情報に基づいて波源方向候補を計算し、算出した波源方向候補ごとに重みを計算する。

　統合部は、複数の周波数別推定方向情報生成部からペアごとの周波数別推定方向情報を取得するとともに、重み生成部から重みを取得する。統合部は、複数のペアごとの周波数別推定方向情報に重みをつけて統合して統合推定方向情報を生成する。

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

　例えば、重み生成部は、複数の周波数統合部、複数の波源方向候補選出部、波源方向候補統合部、複数の波源方向候補別周波数寄与率計算部、および重み計算部を有する。周波数統合部は、周波数別推定方向情報を統合してペアごとの推定方向情報を計算する。波源方向候補選出部は、ペアごとの推定方向情報からピークを検出し、ピークとなる方向をペアごとの波源方向候補として選出する。波源方向候補統合部は、複数の波源方向候補選出部からペアごとの波源方向候補を取得し、取得したペアごとの波源方向候補を統合する。波源方向候補別周波数寄与率計算部は、統合された波源方向候補と周波数別推定方向情報とに基づいて、波源方向候補ごとの推定方向情報に対する周波数別の寄与率を計算する。重み計算部は、波源方向候補ごとの推定方向情報に対する周波数別の寄与率を複数の波源方向候補別周波数寄与率計算部から取得する。重み計算部は、取得した波源方向候補ごとの推定方向情報に対する周波数別の寄与率に基づいて波源方向候補ごとの重みを計算する。

　例えば、重み計算部は、複数のペア間の波源方向候補別の周波数寄与率の類似度を計算し、算出した類似度に基づいて波源方向候補ごとの重みを求める。

　以上で説明したように、本実施形態の波源方向推定装置は、マイクペアごとに波源方向候補を求め、波源方向候補ごとの重みに基づいて統合推定方向情報を求める。そのため、虚像音源を誤推定することなく高精度に音源の方向を推定することができるという効果に加えて、計算量を低減できるという効果を有する。

　（第３の実施形態）
　次に、本発明の第３の実施形態に係る波源方向推定装置について図面を参照しながら説明する。本実施形態の波源方向推定装置は、第１および第２の実施形態の波源推定装置を上位概念化した構成である。なお、本実施形態においては、波源を音源に限定せず、任意の波動の波源の方向を推定対象とする。

　図１４は、本実施形態の波源方向推定装置３０の構成の一例を示すブロック図である。図１４のように、波源方向推定装置３０は、入力端子３１、信号選択部３２、相対遅延時間計算部３３、および統合推定方向情報算出部３５を備える。以下において、第１の実施形態の波源方向推定装置１０や、第２の実施形態の波源方向推定装置２０と同様の構成・動作については説明を省略する場合がある。

　〔入力端子〕
　複数の入力端子３１－１～ｐ（入力部とも呼ぶ）のそれぞれは、図示しないマイクに接続される（ｐは２以上の整数）。複数の入力端子３１－１～ｐのそれぞれには、異なる位置に配置されたセンサによって集音された波動（波動信号とも呼ぶ）から変換された電気信号が入力信号として入力される。なお、入力端子３１－１～ｐは、第１および第２の実施形態の対応する構成と同様であるため、詳細な説明は省略する。

　〔信号選択部〕
　信号選択部３２は、入力端子３１－１～ｐに供給されるｐ個の入力信号のうちの２つの入力信号（入力信号のペアとも呼ぶ）を複数組選択する。信号選択部３２は、選択した複数組の入力信号のペアを統合推定方向情報算出部３５に出力し、それらの入力信号の供給元であるセンサの位置情報（以下、センサ位置情報とも呼ぶ）を相対遅延時間計算部３３に出力する。なお、信号選択部３２は、第１および第２の実施形態の対応する構成と同様であるため、詳細な説明は省略する。

　〔相対遅延時間計算部〕
　相対遅延時間計算部３３には、信号選択部３２からセンサ位置情報が入力される。相対遅延時間計算部３３は、センサ位置情報と音源探索対象方向（音源方向とも呼ぶ）とを用いて、信号選択部３２が選択した全ての入力信号のペア間の相対遅延時間を計算する。相対遅延時間計算部３３は、設定された音源方向と相対遅延時間とをセットにして統合推定方向情報算出部３５に出力する。なお、相対遅延時間計算部３３は、第１および第２の実施形態の対応する構成と同様であるため、詳細な説明は省略する。

　〔統合推定方向情報算出部〕
　統合推定方向情報算出部３５には、複数組の入力信号のペアが信号選択部３２から入力される。また、統合推定方向情報算出部３５には、信号選択部３２によって選択された全ての入力信号のペアの相対遅延時間が相対遅延時間計算部３３から入力される。

　統合推定方向情報算出部３５は、複数組の入力信号のペアと、それらの入力信号のペアの相対遅延時間とを用いて、複数の入力信号のペアに関する周波数ごとの推定方向情報（周波数別推定方向情報とも呼ぶ）を生成する。統合推定方向情報算出部３５は、生成した入力信号のペアごとの周波数別推定方向情報に基づいて、波源方向ごとの推定方向情報に対する周波数寄与率を算出する。統合推定方向情報算出部３５は、算出した周波数寄与率に基づいて、入力信号のペア間の寄与率の類似度に応じた波源方向ごとの重みを計算する。統合推定方向情報算出部３５は、入力信号のペアごとの推定方向情報に重みをつけて統合することによって統合推定方向情報を生成する。統合推定方向情報算出部３５は、生成した統合推定方向情報を出力する。

　第１の実施形態の波源方向推定装置１０に対応させると、統合推定方向情報算出部３５は、推定方向情報・周波数寄与率生成部１５、重み計算部１６、および統合部１７を含む構成である。また、第２の実施形態の波源方向推定装置２０に対応させると、統合推定方向情報算出部３５は、周波数別推定方向情報生成部２５、重み生成部２６、および統合部２７を含む構成である。

　以上のように、本実施形態の波源方向推定装置は、複数の入力部、信号選択部、相対遅延時間計算部、および統合推定方向情報算出部を備える。

　信号選択部は、複数の入力信号のうち二つを組み合わせたペアを複数組選択する。相対遅延時間計算部は、ペアをなす二つの入力信号の供給元のセンサへの波源方向ごとの波動の到達時間差を相対遅延時間として算出する。

　統合推定方向情報算出部は、ペアをなす入力信号とペアの相対遅延時間とを用いて周波数別推定方向情報をペアごとに生成する。統合推定方向情報算出部は生成したペアごとの周波数別推定方向情報を用いて波源方向ごとの推定方向情報に対する周波数別の寄与率を算出する。統合推定方向情報算出部は、算出した寄与率に基づいて異なるペア間の寄与率の類似度に応じた波源方向ごとの重みを計算し、全てのペアの推定方向情報に重みをつけて統合することによって統合推定方向情報を生成する。

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

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

　本実施形態の波源方向推定装置によれば、虚像音源を誤推定することなく、高精度に音源の方向を推定することができる。

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

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

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

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

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

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

　以上のように、本実施形態の波源方向推定装置は、複数の入力部、信号選択部、相対遅延時間計算部、統合推定方向情報算出部、および波源方向算出部を備える。例えば、波源方向算出部は、一定時間ごとに統合推定方向情報が最大になる時刻に対応する方向を波源方向として算出する。本実施形態の波源方向推定装置によれば、虚像音源を誤推定することなく、高精度に音源の方向を推定することができる。

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

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

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

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

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

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

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

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

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

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

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

　以上、実施形態を参照して本発明を説明してきたが、本発明は上記実施形態に限定されるものではない。本発明の構成や詳細には、本発明のスコープ内で当業者が理解し得る様々な変更をすることができる。

　上記の実施形態の一部又は全部は、以下の付記のようにも記載されうるが、以下には限られない。
（付記１）
　複数のセンサによって検出される波動に基づく電気信号を入力信号として取得する複数の入力手段と、
　複数の前記入力信号のうち二つを組み合わせたペアを複数組選択する信号選択手段と、
　前記ペアをなす二つの前記入力信号の供給元の前記センサへの波源方向ごとの波動の到達時間差を相対遅延時間として算出する相対遅延時間計算手段と、
　前記ペアをなす前記入力信号と前記ペアの相対遅延時間とを用いて周波数別推定方向情報を前記ペアごとに生成し、生成した前記ペアごとの前記周波数別推定方向情報を用いて波源方向ごとの推定方向情報に対する周波数別の寄与率を算出し、算出した前記寄与率に基づいて異なる前記ペア間の前記寄与率の類似度に応じた波源方向ごとの重みを計算し、全ての前記ペアの前記推定方向情報に前記重みをつけて統合することによって統合推定方向情報を生成する統合推定方向情報算出手段とを備える波源方向推定装置。
（付記２）
　前記信号選択手段は、
　前記センサの間隔に基づいて、複数の前記入力信号から二つの前記入力信号を組み合わせた前記ペアを選択する付記１に記載の波源方向推定装置。
（付記３）
　前記相対遅延時間計算手段は、
　一つの前記ペアをなす二つの前記入力信号の供給元である前記センサにおける波源探索方向を基準として、前記信号選択手段によって選択された全ての前記ペアの前記相対遅延時間を、基準である前記波源探索方向の関数として計算する付記１または２に記載の波源方向推定装置。
（付記４）
　前記統合推定方向情報算出手段は、
　前記信号選択手段によって選択された全ての前記ペアのうち１組の前記ペアをなす二つの前記入力信号を取得するとともに、取得した前記ペアの前記相対遅延時間を取得し、前記ペアをなす二つの前記入力信号間の周波数別相互相関関数と前記ペアの前記相対遅延時間とを用いて、前記ペアをなす二つの前記入力信号間の前記推定方向情報と周波数別の前記寄与率とを生成する複数の推定方向情報・周波数寄与率生成手段と、
　複数の前記推定方向情報・周波数寄与率生成手段から周波数ごとの前記寄与率を取得し、異なる前記ペア間の前記寄与率の類似度に応じて波源方向ごとの前記重みを計算する重み計算手段と、
　複数の前記推定方向情報・周波数寄与率生成手段から前記推定方向情報を取得するとともに、前記重み計算手段から前記重みを取得し、複数の前記ペアごとの前記推定方向情報に前記重みをつけて統合して前記統合推定方向情報を生成する統合手段とを有する付記１乃至３のいずれか一項に記載の波源方向推定装置。
（付記５）
　前記推定方向情報・周波数寄与率生成手段は、
　前記ペアをなす二つの前記入力信号を周波数領域の変換信号に変換する変換手段と、
　前記変換手段によって変換された前記変換信号を用いてクロススペクトルを計算するクロススペクトル計算手段と、
　前記クロススペクトル計算手段によって算出される前記クロススペクトルを用いて平均クロススペクトルを計算する平均計算手段と、
　前記平均計算手段によって算出される前記平均クロススペクトルを用いて分散を計算する分散計算手段と、
　前記平均計算手段によって算出される前記平均クロススペクトルと、前記分散計算手段によって算出される分散とを用いて周波数別クロススペクトルを計算する周波数別クロススペクトル生成手段と、
　前記周波数別クロススペクトル生成手段によって算出される前記周波数別クロススペクトルを逆変換して周波数別相互相関関数を計算する逆変換手段と、
　前記ペアの前記周波数別相互相関関数と前記相対遅延時間とに基づいて波源方向と相関値との対応関係である周波数別推定方向情報を計算する周波数別推定方向情報計算手段と、
　前記周波数別推定方向情報を統合して前記推定方向情報を計算する周波数統合手段と、
　前記周波数別推定方向情報に基づいて波源方向ごとの前記推定方向情報に対する周波数別の前記寄与率を計算する周波数寄与率計算手段とを有する付記４に記載の波源方向推定装置。
（付記６）
　前記重み計算手段は、
　複数の前記ペア間において前記寄与率の類似度を計算し、算出した前記類似度に基づいて波源方向ごとの重みを計算する付記４または５に記載の波源方向推定装置。
（付記７）
　前記統合推定方向情報算出手段は、
　前記信号選択手段によって選択された全ての前記ペアのうち１組の前記ペアをなす二つの前記入力信号を取得するとともに、取得した前記ペアの前記相対遅延時間を取得し、前記ペアをなす二つの前記入力信号を用いて算出される周波数別相互相関関数と前記ペアの前記相対遅延時間とを用いて、前記ペアごとの周波数別推定方向情報を計算する複数の周波数別推定方向情報生成手段と、
　複数の前記周波数別推定方向情報生成手段から前記ペアごとの前記周波数別推定方向情報を取得し、取得した前記ペアごとの前記周波数別推定方向情報に基づいて波源方向候補を計算し、算出した前記波源方向候補ごとに前記重みを計算する重み生成手段と、
　複数の前記周波数別推定方向情報生成手段から前記ペアごとの前記周波数別推定方向情報を取得するとともに、前記重み生成手段から前記重みを取得し、複数の前記ペアごとの前記周波数別推定方向情報に前記重みをつけて統合して前記統合推定方向情報を生成する統合手段とを有する付記１乃至３のいずれか一項に記載の波源方向推定装置。
（付記８）
　前記周波数別推定方向情報生成手段は、
　前記ペアをなす二つの前記入力信号を周波数領域の変換信号に変換する変換手段と、
　前記変換手段によって変換された前記変換信号を用いてクロススペクトルを計算するクロススペクトル計算手段と、
　前記クロススペクトル計算手段によって算出される前記クロススペクトルを用いて平均クロススペクトルを計算する平均計算手段と、
　前記平均計算手段によって算出される前記平均クロススペクトルを用いて分散を計算する分散計算手段と、
　前記平均計算手段によって算出される前記平均クロススペクトルと、前記分散計算手段によって算出される分散とを用いて周波数別クロススペクトルを計算する周波数別クロススペクトル生成手段と、
　前記周波数別クロススペクトル生成手段によって算出される前記周波数別クロススペクトルを逆変換して周波数別相互相関関数を計算する逆変換手段と、
　前記ペアの前記周波数別相互相関関数と前記相対遅延時間とに基づいて波源方向と相関値との対応関係である周波数別推定方向情報を計算する周波数別推定方向情報計算手段とを有する付記７に記載の波源方向推定装置。
（付記９）
　前記重み生成手段は、
　前記周波数別推定方向情報を統合して前記ペアごとの前記推定方向情報を計算する複数の周波数統合手段と、
　前記ペアごとの前記推定方向情報からピークを検出し、ピークとなる方向を前記ペアごとの波源方向候補として選出する複数の波源方向候補選出手段と、
　複数の前記波源方向候補選出手段から前記ペアごとの前記波源方向候補を取得し、取得した前記ペアごとの前記波源方向候補を統合する波源方向候補統合手段と、
　統合された前記波源方向候補と前記周波数別推定方向情報とに基づいて、前記波源方向候補ごとの前記推定方向情報に対する周波数別の前記寄与率を計算する複数の波源方向候補別周波数寄与率計算手段と、
　前記波源方向候補ごとの前記推定方向情報に対する周波数別の前記寄与率を複数の前記波源方向候補別周波数寄与率計算手段から取得し、取得した前記波源方向候補ごとの前記推定方向情報に対する周波数別の前記寄与率に基づいて前記波源方向候補ごとの前記重みを計算する重み計算手段とを有する付記８に記載の波源方向推定装置。
（付記１０）
　前記重み計算手段は、
　複数の前記ペア間の前記波源方向候補ごとの周波数別の前記寄与率の類似度を計算し、算出した前記類似度に基づいて前記波源方向候補ごとの前記重みを求める付記９に記載の波源方向推定装置。
（付記１１）
　前記統合推定方向情報算出手段によって算出される前記統合推定方向情報に基づいて波源方向を算出する波源方向算出手段を備える付記１乃至１０のいずれか一項に記載の波源方向推定装置。
（付記１２）
　前記波源方向算出手段は、
　一定時間ごとに前記統合推定方向情報が最大になる時刻に対応する方向を波源方向として算出する付記１１に記載の波源方向推定装置。
（付記１３）
　複数の前記入力手段のそれぞれに対応して配置される前記センサを備える付記１乃至１２のいずれか一項に記載の波源方向推定装置。
（付記１４）
　情報処理装置が、
　複数のセンサによって検出される波動に基づく電気信号を入力信号として取得し、
　複数の前記入力信号のうち二つを組み合わせたペアを複数組選択し、
　前記ペアをなす二つの前記入力信号の供給元の前記センサへの波源方向ごとの波動の到達時間差を相対遅延時間として算出し、
　前記ペアをなす前記入力信号と前記ペアの相対遅延時間とを用いて周波数別推定方向情報を前記ペアごとに生成し、
　生成した前記ペアごとの前記周波数別推定方向情報を用いて波源方向ごとの推定方向情報に対する周波数別の寄与率を算出し、
　算出した前記寄与率に基づいて異なる前記ペア間の前記寄与率の類似度に応じた波源方向ごとの重みを計算し、
　全ての前記ペアの前記推定方向情報に前記重みをつけて統合することによって統合推定方向情報を生成する波源方向推定方法。
（付記１５）
　複数のセンサによって検出される波動に基づく電気信号を入力信号として取得する処理と、
　複数の前記入力信号のうち二つを組み合わせたペアを複数組選択する処理と、
　前記ペアをなす二つの前記入力信号の供給元の前記センサへの波源方向ごとの波動の到達時間差を相対遅延時間として算出する処理と、
　前記ペアをなす前記入力信号と前記ペアの相対遅延時間とを用いて周波数別推定方向情報を前記ペアごとに生成する処理と、
　生成した前記ペアごとの前記周波数別推定方向情報を用いて波源方向ごとの推定方向情報に対する周波数別の寄与率を算出する処理と、
　算出した前記寄与率に基づいて異なる前記ペア間の前記寄与率の類似度に応じた波源方向ごとの重みを計算する処理と、
　全ての前記ペアの前記推定方向情報に前記重みをつけて統合することによって統合推定方向情報を生成する処理とをコンピュータに実行させるプログラムを記録させたプログラム記録媒体。

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

Claims

　複数のセンサによって検出される波動に基づく電気信号を入力信号として取得する複数の入力手段と、
　複数の前記入力信号のうち二つを組み合わせたペアを複数組選択する信号選択手段と、
　前記ペアをなす二つの前記入力信号の供給元の前記センサへの波源方向ごとの波動の到達時間差を相対遅延時間として算出する相対遅延時間計算手段と、
　前記ペアをなす前記入力信号と前記ペアの相対遅延時間とを用いて周波数別推定方向情報を前記ペアごとに生成し、生成した前記ペアごとの前記周波数別推定方向情報を用いて波源方向ごとの推定方向情報に対する周波数別の寄与率を算出し、算出した前記寄与率に基づいて異なる前記ペア間の前記寄与率の類似度に応じた波源方向ごとの重みを計算し、全ての前記ペアの前記推定方向情報に前記重みをつけて統合することによって統合推定方向情報を生成する統合推定方向情報算出手段とを備える波源方向推定装置。
　前記信号選択手段は、
　前記センサの間隔に基づいて、複数の前記入力信号から二つの前記入力信号を組み合わせた前記ペアを選択する請求項１に記載の波源方向推定装置。
　前記相対遅延時間計算手段は、
　一つの前記ペアをなす二つの前記入力信号の供給元である前記センサにおける波源探索方向を基準として、前記信号選択手段によって選択された全ての前記ペアの前記相対遅延時間を、基準である前記波源探索方向の関数として計算する請求項１または２に記載の波源方向推定装置。
　前記統合推定方向情報算出手段は、
　前記信号選択手段によって選択された全ての前記ペアのうち１組の前記ペアをなす二つの前記入力信号を取得するとともに、取得した前記ペアの前記相対遅延時間を取得し、前記ペアをなす二つの前記入力信号間の周波数別相互相関関数と前記ペアの前記相対遅延時間とを用いて、前記ペアをなす二つの前記入力信号間の前記推定方向情報と周波数別の前記寄与率とを生成する複数の推定方向情報・周波数寄与率生成手段と、
　複数の前記推定方向情報・周波数寄与率生成手段から周波数ごとの前記寄与率を取得し、異なる前記ペア間の前記寄与率の類似度に応じて波源方向ごとの前記重みを計算する重み計算手段と、
　複数の前記推定方向情報・周波数寄与率生成手段から前記推定方向情報を取得するとともに、前記重み計算手段から前記重みを取得し、複数の前記ペアごとの前記推定方向情報に前記重みをつけて統合して前記統合推定方向情報を生成する統合手段とを有する請求項１乃至３のいずれか一項に記載の波源方向推定装置。
　前記推定方向情報・周波数寄与率生成手段は、
　前記ペアをなす二つの前記入力信号を周波数領域の変換信号に変換する変換手段と、
　前記変換手段によって変換された前記変換信号を用いてクロススペクトルを計算するクロススペクトル計算手段と、
　前記クロススペクトル計算手段によって算出される前記クロススペクトルを用いて平均クロススペクトルを計算する平均計算手段と、
　前記平均計算手段によって算出される前記平均クロススペクトルを用いて分散を計算する分散計算手段と、
　前記平均計算手段によって算出される前記平均クロススペクトルと、前記分散計算手段によって算出される分散とを用いて周波数別クロススペクトルを計算する周波数別クロススペクトル生成手段と、
　前記周波数別クロススペクトル生成手段によって算出される前記周波数別クロススペクトルを逆変換して周波数別相互相関関数を計算する逆変換手段と、
　前記ペアの前記周波数別相互相関関数と前記相対遅延時間とに基づいて波源方向と相関値との対応関係である周波数別推定方向情報を計算する周波数別推定方向情報計算手段と、
　前記周波数別推定方向情報を統合して前記推定方向情報を計算する周波数統合手段と、
　前記周波数別推定方向情報に基づいて波源方向ごとの前記推定方向情報に対する周波数別の前記寄与率を計算する周波数寄与率計算手段とを有する請求項４に記載の波源方向推定装置。
　前記重み計算手段は、
　複数の前記ペア間において前記寄与率の類似度を計算し、算出した前記類似度に基づいて波源方向ごとの重みを計算する請求項４または５に記載の波源方向推定装置。
　前記統合推定方向情報算出手段は、
　前記信号選択手段によって選択された全ての前記ペアのうち１組の前記ペアをなす二つの前記入力信号を取得するとともに、取得した前記ペアの前記相対遅延時間を取得し、前記ペアをなす二つの前記入力信号を用いて算出される周波数別相互相関関数と前記ペアの前記相対遅延時間とを用いて、前記ペアごとの周波数別推定方向情報を計算する複数の周波数別推定方向情報生成手段と、
　複数の前記周波数別推定方向情報生成手段から前記ペアごとの前記周波数別推定方向情報を取得し、取得した前記ペアごとの前記周波数別推定方向情報に基づいて波源方向候補を計算し、算出した前記波源方向候補ごとに前記重みを計算する重み生成手段と、
　複数の前記周波数別推定方向情報生成手段から前記ペアごとの前記周波数別推定方向情報を取得するとともに、前記重み生成手段から前記重みを取得し、複数の前記ペアごとの前記周波数別推定方向情報に前記重みをつけて統合して前記統合推定方向情報を生成する統合手段とを有する請求項１乃至３のいずれか一項に記載の波源方向推定装置。
　前記周波数別推定方向情報生成手段は、
　前記ペアをなす二つの前記入力信号を周波数領域の変換信号に変換する変換手段と、
　前記変換手段によって変換された前記変換信号を用いてクロススペクトルを計算するクロススペクトル計算手段と、
　前記クロススペクトル計算手段によって算出される前記クロススペクトルを用いて平均クロススペクトルを計算する平均計算手段と、
　前記平均計算手段によって算出される前記平均クロススペクトルを用いて分散を計算する分散計算手段と、
　前記平均計算手段によって算出される前記平均クロススペクトルと、前記分散計算手段によって算出される分散とを用いて周波数別クロススペクトルを計算する周波数別クロススペクトル生成手段と、
　前記周波数別クロススペクトル生成手段によって算出される前記周波数別クロススペクトルを逆変換して周波数別相互相関関数を計算する逆変換手段と、
　前記ペアの前記周波数別相互相関関数と前記相対遅延時間とに基づいて波源方向と相関値との対応関係である周波数別推定方向情報を計算する周波数別推定方向情報計算手段とを有する請求項７に記載の波源方向推定装置。
　前記重み生成手段は、
　前記周波数別推定方向情報を統合して前記ペアごとの前記推定方向情報を計算する複数の周波数統合手段と、
　前記ペアごとの前記推定方向情報からピークを検出し、ピークとなる方向を前記ペアごとの波源方向候補として選出する複数の波源方向候補選出手段と、
　複数の前記波源方向候補選出手段から前記ペアごとの前記波源方向候補を取得し、取得した前記ペアごとの前記波源方向候補を統合する波源方向候補統合手段と、
　統合された前記波源方向候補と前記周波数別推定方向情報とに基づいて、前記波源方向候補ごとの前記推定方向情報に対する周波数別の前記寄与率を計算する複数の波源方向候補別周波数寄与率計算手段と、
　前記波源方向候補ごとの前記推定方向情報に対する周波数別の前記寄与率を複数の前記波源方向候補別周波数寄与率計算手段から取得し、取得した前記波源方向候補ごとの前記推定方向情報に対する周波数別の前記寄与率に基づいて前記波源方向候補ごとの前記重みを計算する重み計算手段とを有する請求項８に記載の波源方向推定装置。
　前記重み計算手段は、
　複数の前記ペア間の前記波源方向候補ごとの周波数別の前記寄与率の類似度を計算し、算出した前記類似度に基づいて前記波源方向候補ごとの前記重みを求める請求項９に記載の波源方向推定装置。
　前記統合推定方向情報算出手段によって算出される前記統合推定方向情報に基づいて波源方向を算出する波源方向算出手段を備える請求項１乃至１０のいずれか一項に記載の波源方向推定装置。
　前記波源方向算出手段は、
　一定時間ごとに前記統合推定方向情報が最大になる時刻に対応する方向を波源方向として算出する請求項１１に記載の波源方向推定装置。
　複数の前記入力手段のそれぞれに対応して配置される前記センサを備える請求項１乃至１２のいずれか一項に記載の波源方向推定装置。
　情報処理装置が、
　複数のセンサによって検出される波動に基づく電気信号を入力信号として取得し、
　複数の前記入力信号のうち二つを組み合わせたペアを複数組選択し、
　前記ペアをなす二つの前記入力信号の供給元の前記センサへの波源方向ごとの波動の到達時間差を相対遅延時間として算出し、
　前記ペアをなす前記入力信号と前記ペアの相対遅延時間とを用いて周波数別推定方向情報を前記ペアごとに生成し、
　生成した前記ペアごとの前記周波数別推定方向情報を用いて波源方向ごとの推定方向情報に対する周波数別の寄与率を算出し、
　算出した前記寄与率に基づいて異なる前記ペア間の前記寄与率の類似度に応じた波源方向ごとの重みを計算し、
　全ての前記ペアの前記推定方向情報に前記重みをつけて統合することによって統合推定方向情報を生成する波源方向推定方法。
　複数のセンサによって検出される波動に基づく電気信号を入力信号として取得する処理と、
　複数の前記入力信号のうち二つを組み合わせたペアを複数組選択する処理と、
　前記ペアをなす二つの前記入力信号の供給元の前記センサへの波源方向ごとの波動の到達時間差を相対遅延時間として算出する処理と、
　前記ペアをなす前記入力信号と前記ペアの相対遅延時間とを用いて周波数別推定方向情報を前記ペアごとに生成する処理と、
　生成した前記ペアごとの前記周波数別推定方向情報を用いて波源方向ごとの推定方向情報に対する周波数別の寄与率を算出する処理と、
　算出した前記寄与率に基づいて異なる前記ペア間の前記寄与率の類似度に応じた波源方向ごとの重みを計算する処理と、
　全ての前記ペアの前記推定方向情報に前記重みをつけて統合することによって統合推定方向情報を生成する処理とをコンピュータに実行させるプログラムを記録させたプログラム記録媒体。