(57)【要約】  周囲のノイズレベルが高い環境の場合であっても、明瞭なピークを持つ相関関数を生成する。相関関数生成装置であって、波源で発生した波を入力信号として取得する複数の入力信号取得部と、入力信号取得部で取得した複数の入力信号を複数の周波数領域信号に変換する変換部と、周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算部と、クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算部と、周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算部と、を備えた。
2018003158.tif 000002
【特許請求の範囲】
【請求項１】
波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
前記入力信号取得手段で取得した複数の前記入力信号を複数の周波数領域信号に変換する変換手段と、
前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算手段と、
前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算手段と、
を備えた相関関数生成装置。
【請求項２】
前記統合相関関数計算手段は、
前記周波数別クロススペクトルの逆変換により周波数別相関関数を生成する周波数別相関関数生成手段と、
前記周波数別相関関数を統合して一つの統合相関関数を生成する統合相関関数生成手段と、
を有する請求項１に記載の相関関数生成装置。
【請求項３】
前記統合相関関数計算手段は、
前記周波数別クロススペクトルを統合して統合クロススペクトルを生成する統合クロススペクトル生成手段と、
前記統合クロススペクトルの逆変換により統合相関関数を生成する統合相関関数生成手段と、
を有する請求項１に記載の相関関数生成装置。
【請求項４】
前記周波数別クロススペクトル計算手段は、
前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段を有し、
前記周波数別基本クロススペクトルを前記周波数別クロススペクトルとして求める請求項１乃至３のいずれか１項に記載の相関関数生成装置。
【請求項５】
前記周波数別クロススペクトル計算手段は、
前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
カーネル関数スペクトルを記憶するカーネル関数記憶手段と、
前記周波数別基本クロススペクトルと前記カーネル関数スペクトルとを乗じて、前記周波数別クロススペクトルを求める乗算手段と、
を有する請求項１乃至３のいずれか１項に記載の相関関数生成装置。
【請求項６】
前記周波数別クロススペクトル計算手段は、
前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
前記クロススペクトルに基づいて、カーネル関数スペクトルを計算するカーネル関数スペクトル計算手段と、
前記周波数別基本クロススペクトルと前記カーネル関数スペクトルとを乗じて、前記周波数別クロススペクトルを求める乗算手段と、
を有する請求項１乃至３のいずれか１項に記載の相関関数生成装置。
【請求項７】
波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
を含む相関関数生成方法。
【請求項８】
波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
をコンピュータに実行させる相関関数生成プログラム。
【請求項９】
請求項１乃至６のいずれか１項に記載の相関関数生成装置と、
統合相関関数に基づいて、波源の推定方向情報を生成する推定方向情報生成手段と、
を備えた波源方向推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置に関する。
【背景技術】
【０００２】
上記技術分野において、非特許文献１および非特許文献２には、２つのマイクロフォンの受音信号を用いて、音源(音波の発生源、発生場所)の方向を推定する方法が記載されている。具体的には、２つの受音信号から、受音信号間の相互相関関数を求める。そして、相互相関関数が最大値を与える時間差を音波の到来時間差として算出することで、音波の到来方向を推定する技術が開示されている。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】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.
【非特許文献３】Y. Ephraim and D. Malah, “Speech enhancement using a minimum mean-square error short-time spectral amplitude estimator,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-32, no. 6, pp. 1109-1121, Dec. 1984.
【非特許文献４】R. Martin, "Spectral subtraction based on minimum statistics," Proc. of EUSPICO-94, pp.1182-1185, Sept. 1994
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、上記文献に記載の技術では、周囲のノイズレベルが高い環境の場合には、明瞭なピークを持つ相関関数が生成できなかった。また、波源の方向を高精度に推定することができなかった。【０００５】
本発明の目的は、上述の課題を解決する技術を提供することにある。
【課題を解決するための手段】
【０００６】
上記目的を達成するため、本発明に係る相関関数生成装置は、
波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
前記入力信号取得手段で取得した複数の前記入力信号を複数の周波数領域信号に変換する変換手段と、
前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算手段と、
前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算手段と、
を備えた。【０００７】
上記目的を達成するため、本発明に係る相関関数生成方法は、
波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
を含む。【０００８】
上記目的を達成するため、本発明に係る相関関数生成プログラムは、
波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
をコンピュータに実行させる。【０００９】
上記目的を達成するため、本発明に係る波源方向推定装置は、
上記相関関数生成装置と、
統合相関関数に基づいて、波源の推定方向情報を生成する推定方向情報生成手段と、
を備えた。
【発明の効果】
【００１０】
本発明によれば、周囲のノイズレベルが高い環境の場合であっても、明瞭なピークを持つ相関関数を生成できる。また、波源の方向を高精度に推定することができる。
【図面の簡単な説明】
【００１１】
【図１】本発明の第１実施形態に係る情報処理装置の構成を示すブロック図である。
【図２Ａ】本発明の第２実施形態に係る波源方向推定装置の構成を示すブロック図である。
【図２Ｂ】本発明の第２実施形態に係る波源方向推定装置の備える周波数別クロススペクトル計算部の構成を示すブロック図である。
【図２Ｃ】本発明の第２実施形態に係る波源方向推定装置の備える統合相関関数計算部の構成を示すブロック図である。
【図３Ａ】本発明の第２実施形態に係る波源方向推定装置により得られた周波数別相関関数の一例について説明する図である。
【図３Ｂ】本発明の第２実施形態に係る波源方向推定装置により得られた周波数別相関関数を統合した統合相関関数の一例について説明する図である。
【図４】本発明の第２実施形態に係る波源方向推定装置の備える統合相関関数テーブルの構成の一例を示す図である。
【図５】本発明の第２実施形態に係る波源方向推定装置のハードウェア構成を示すブロック図である。
【図６】本発明の第２実施形態に係る波源方向推定装置の処理手順を説明するフローチャートである。
【図７】本発明の第３実施形態に係る波源方向推定装置の備える統合相関関数生成部の構成を示すブロック図である。
【図８Ａ】本発明の第４実施形態に係る波源方向推定装置の構成を示すブロック図である。
【図８Ｂ】本発明の第４実施形態に係る波源方向推定装置の備える周波数別クロススペクトル計算部の構成を示すブロック図である。
【図９】本発明の第４実施形態に係る波源方向推定装置の周波数別クロススペクトル計算部における、カーネル関数スペクトルを乗じた周波数別クロススペクトルと、周波数別相関関数との関係を示す図である。
【図１０】本発明の第４実施形態に係る波源方向推定装置の周波数別クロススペクトル計算部における、カーネル関数による周波数別相関関数の大きさ制御の効果について説明する図である。
【図１１】本発明の第４実施形態に係る波源方向推定装置の周波数別クロススペクトル計算部における、カーネル関数スペクトル幅の違いと、統合相関関数との関係を示す図である。
【図１２】本発明の第５実施形態に係る波源方向推定装置の備える周波数別クロススペクトル計算部の構成を説明するための図である。
【図１３Ａ】本発明の第６実施形態に係る波源方向推定装置の構成を説明するためのブロック図である。
【図１３Ｂ】本発明の第６実施形態に係る波源方向推定装置の備える周波数別クロススペクトルの構成を説明するためのブロック図である。
【図１４】本発明の第７実施形態に係る波源方向推定システムの構成を説明するための図である。
【図１５】本発明の第７実施形態に係る波源方向推定システムの表示部に表示される画像の一例を示す図である。
【発明を実施するための形態】
【００１２】
以下に、本発明を実施するための形態について、図面を参照して、例示的に詳しく説明記載する。ただし、以下の実施の形態に記載されている、構成、数値、処理の流れ、機能要素などは一例に過ぎず、その変形や変更は自由であって、本発明の技術範囲を以下の記載に限定する趣旨のものではない。
【００１３】
また、以下の実施形態に係る波源方向推定装置の推定対象は、空気や水の振動波である音波の発生源に限定されない。地震や地滑りなどの土や固体を媒質とする振動波の発生源にも適用できる。その場合、振動波を電気信号に変換する装置には、マイクロフォンではなく振動センサが用いられる。さらに、気体・液体・固体の振動波だけでなく、電波を用いて方向を推定する場合にも以下の実施形態に係る波源方向推定装置を適用できる。その場合、電波を電気信号に変換する装置にはアンテナが用いられる。以下の実施形態においては、波源は音源と仮定して説明する。
【００１４】
［第１実施形態］
本発明の第１実施形態としての相関関数生成装置１００について、図１を用いて説明する。相関関数生成装置１００は、入力信号に基づいて相関関数を生成する装置である。
【００１５】
図１に示すように、相関関数生成装置１００は、入力信号取得部１０１と、変換部１０２と、クロススペクトル計算部１０３と、周波数別クロススペクトル計算部１０４と、統合相関関数計算部１０５とを含む。
【００１６】
複数の入力信号取得部１０１は、波源で発生した波を入力信号として取得する。変換部１０２は、入力信号取得手段で取得した複数の前記入力信号を複数の周波数領域信号に変換する。クロススペクトル計算部１０３は、周波数領域信号に基づいて、クロススペクトルを計算する。周波数別クロススペクトル計算部１０４は、クロススペクトルに基づいて、周波数別クロススペクトルを計算する。統合相関関数計算部１０５は、周波数別クロススペクトルに基づいて、統合相関関数を計算する。
【００１７】
本実施形態によれば、周囲のノイズレベルが高い環境の場合であっても、明瞭なピークを持つ相関関数を生成できる。また、波源の方向を高精度に推定することができる。
【００１８】
［第２実施形態］
次に本発明の第２実施形態に係る波源方向推定装置について、図２Ａ乃至図６を用いて説明する。
【００１９】
＜前提技術＞
上記非特許文献１および非特許文献２に記載の技術では、屋外などの周囲騒音レベルが高い環境において、遠方に存在する音源の方向を高精度に推定することが困難であった。例えば、推定対象の音源（目標音源）が、マイクロフォンから遠く離れた場所に存在する場合、目標音源から放射される音の音量は、マイクロフォンに到達する時点で大幅に小さくなる。このため、目標音源の音が周囲環境雑音に埋もれて、明瞭なピークを持つ相関関数を生成できなかった。このため、目標音源の方向推定精度が低下することがあった。
【００２０】
＜本実施形態の技術＞図２Ａは、本実施形態に係る波源方向推定装置の構成を示すブロック図である。図２Ｂは、本実施形態に係る波源方向推定装置の備える統合相関関数計算部の構成を示すブロック図である。
【００２１】
本実施形態に係る波源方向推定装置２００は、例えば、デジタルビデオカメラ、スマートフォン、携帯電話、ノートパソコン、パッシブソーナーなどといった装置の一部として機能する。また、不審ドローン検知、悲鳴検知、車両事故検知などの声や音に基づいて異常を検知する異常音検知装置にも搭載される。しかし、本実施形態に係る波源方向推定装置２００の適用例はこれらに限定されるものではなく、受信音から目標音源の方向推定を要求されるあらゆる波源方向推定装置に適用可能である。
【００２２】
波源方向推定装置２００は、入力端子２０１と、入力端子２０２と、変換部２０１と、クロススペクトル計算部２０２と、周波数別クロススペクトル計算部２０３１〜２０３ｋと、を備える。波源方向推定装置２００は、さらに、統合相関関数計算部２０４と、推定方向情報生成部２０５と、相対遅延時間計算部２０６とを備える。
【００２３】
入力端子２０１と入力端子２０２とには、目標音源の音と、集音装置であるマイクロフォン(以下、マイク)の周囲で生じる様々な雑音が混在した音とがデジタル信号(サンプル値系列)として入力される。入力端子２０１と入力端子２０２とに入力された音信号を本実施形態では入力信号と呼ぶ。そして、時刻ｔにおける入力端子２０１の入力信号をｘ１（ｔ）と、入力端子２０２の入力信号をｘ２（ｔ）と表す。
【００２４】
入力端子に入力される音は、集音装置であるマイクで集音される。入力端子は複数存在するので、目標音源の音を集音する場合には、端子数と同じ２個のマイクが同時に使用される。本実施形態では、入力端子とマイクとは一対一に対応しており、ｍ番目のマイクが集音した音は、ｍ番目の入力端子に供給されると仮定する。したがって、ｍ番目の入力端子に供給された入力信号のことを「ｍ番目のマイクの入力信号」とも呼ぶ。
【００２５】
波源方向推定装置２００は、目標音源の音が２つのマイクに到達する時間差を利用して音源の方向を推定する。このため、マイク間隔も重要な情報となるので、波源方向推定装置２００には、入力信号だけでなく、マイク位置情報も供給される。
【００２６】
変換部２０１は、入力端子２０１と入力端子２０２とから供給された入力信号を変換し、クロススペクトル計算部２０２へ供給する。変換は、入力信号を複数の周波数成分に分解する目的で実施される。ここでは、代表的なフーリエ変換を用いた場合について説明する。
【００２７】
変換部２０１には２種類の入力信号ｘｍ（ｔ）が入力される。ここで、ｍは入力端子番号である。変換部２０１は、入力端子から供給された入力信号から、適当な長さの波形を一定の周期でずらしながら切り出す。こうして切り出した信号区間をフレーム、切り出した波形の長さをフレーム長、フレームをずらす周期をフレーム周期と呼ぶ。そして、フーリエ変換を用いて切り出された信号を周波数領域信号に変換する。ｎをフレーム番号とし、切り出した入力信号ｘｍ（ｔ，ｎ）（ｔ＝０，１，・・・，Ｋ−１）とすると、ｘｍ（ｔ，ｎ）のフーリエ変換Ｘｍ（ｋ，ｎ）は次のように計算される。
【００２８】
【数１】
2018003158.tif 000003
ここで、ｊは虚数単位（−１の平方根）、ｅｘｐは指数関数を表す。また、ｋは周波数ビン番号を表し、０以上Ｋ−１以下の整数である。以下、簡略化のため、ｋのことを周波数ビン番号ではなく、単に「周波数」と呼ぶ。
【００２９】
クロススペクトル計算部２０２は、変換部２０１から供給される変換信号に基づいて、クロススペクトルを計算し、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋへ伝達する。クロススペクトル計算部２０２は、変換信号Ｘ２（ｋ，ｎ）の複素共役と変換信号Ｘ１（ｋ，ｎ）との積を計算する。変換信号のクロススペクトルをＳ１２（ｋ，ｎ）とすると、クロススペクトルは次のように計算される。
【００３０】
【数２】
2018003158.tif 000004
ここで、ｃｏｎｊ（Ｘ２（ｋ，ｎ））は、Ｘ２（ｋ，ｎ）の複素共役を表す。
【００３１】
＜周波数別クロススペクトル計算部＞
周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ１２（ｋ，ｎ）を用いて、Ｓ１２（ｋ，ｎ）の各周波数ｋに対応するクロススペクトルを計算し、周波数別クロススペクトルとして統合相関関数計算部２０４に伝達する。周波数別クロススペクトルの計算は、周波数成分ごとに相関関数を計算するために行われる。つまり、ある周波数ｋに対応する相関関数(周波数別相関関数と呼ぶ)を後段で求めるために、周波数別クロススペクトルを計算する。
【００３２】
次に、ある周波数kの周波数別クロススペクトルを計算する周波数別クロススペクトル計算部２０３ｋについて詳細に説明する。図２Ｂは、周波数別クロススペクトル計算部２０３ｋのブロック図である。周波数別クロススペクトル計算部２０３ｋは、周波数別基本クロススペクトル計算部２０３１ｋを有する。周波数別クロススペクトル計算部２０３ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ１２（ｋ，ｎ）を用いて、周波数別基本クロススペクトルを計算し、周波数別クロススペクトルとして統合相関関数計算部２０４へ伝達する。周波数別基本クロススペクトル計算部２０３１ｋでは、周波数ｋのクロススペクトルＳ１２（ｋ，ｎ）に基づいて、周波数別基本クロススペクトルを計算する際、位相成分と振幅成分とをあらかじめ別々に求めたのちに統合する。周波数ｋの周波数別基本クロススペクトルＵｋ（ｗ，ｎ）、その振幅成分を｜Ｕｋ（ｗ，ｎ）｜、位相成分ａｒｇ（Ｕｋ（ｗ，ｎ））とすると、次のような関係が成立する。
【００３３】
【数３】
2018003158.tif 000005
ここで、ｗは周波数を表し、０以上Ｗ−１以下の整数である。以下、周波数ｋのクロススペクトルＳ１２（ｋ，ｎ）から、周波数別基本クロススペクトルの振幅成分｜Ｕｋ（ｗ，ｎ）｜と位相成分ａｒｇ（Ｕｋ（ｗ，ｎ））とを求める方法について説明する。
【００３４】
振幅成分｜Ｕｋ（ｗ，ｎ）｜は、ｋを整数倍した周波数には、１．０を用いる。一方、周波数kの非定数倍した周波数の位相成分はゼロにする。これを数式で表現すると、振幅成分｜Ｕｋ（ｗ，ｎ）｜は以下のように与えられる。
【００３５】
【数４】
2018003158.tif 000006
ここで、ｐは１以上Ｐ以下の整数である。波源方向推定を行うときに重要な情報は位相成分であるため、振幅成分はこのように適当な定数を用いる。その他には、１．０の代わりに｜Ｓ１２（ｋ，ｎ）｜を用いてもよい。つまり、振幅成分｜Ｕｋ（ｗ，ｎ）｜を次式のように求めてもよい。
【００３６】
【数５】
2018003158.tif 000007
位相成分ａｒｇ（Ｕｋ（ｗ，ｎ））は、ｋを整数倍した周波数には、周波数ｋのクロススペクトルＳ１２（ｋ，ｎ）を定数倍したものを用いる。例えば、周波数ｋ、２ｋ、３ｋ、４ｋの位相成分には、周波数ｋの位相成分ａｒｇ（Ｓ１２（ｋ，ｎ））をそれぞれ同一の倍率で整数倍したもの、つまりａｒｇ（Ｓ１２（ｋ，ｎ））、２ａｒｇ（Ｓ１２（ｋ，ｎ））、３ａｒｇ（Ｓ１２（ｋ，ｎ））、４ａｒｇ（Ｓ１２（ｋ，ｎ））を用いる。一方、周波数ｋの非定数倍した周波数の位相成分はゼロにする。したがって、周波数ｋに対応する周波数別基本クロススペクトルの位相成分ａｒｇ（Ｕｋ（ｗ，ｎ））は、以下のように計算される。
【００３７】
【数６】
2018003158.tif 000008
ここで、ｐは１以上Ｐ以下の整数である。また、Ｐは１よりも大きい整数である。
【００３８】
以上の方法で求めた振幅成分と位相成分とを、上述の式（３）を用いて統合し、周波数ｋの周波数別基本クロススペクトルＵｋ（ｗ，ｎ）を得る。
【００３９】
これまで説明してきた方法では、振幅成分と位相成分とを別々に求めてから周波数別スペクトルを得ていた。しかし、次に示す数式に示すようにクロススペクトルのべき乗を用いれば、振幅成分と位相成分とを求めることなく周波数別スペクトルＵｋ（ｗ，ｎ）を求めることが可能である。
【００４０】
【数７】
2018003158.tif 000009
統合相関関数計算部２０４は、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋから供給される周波数別クロススペクトルに基づき統合相関関数を計算し、推定方向情報生成部２０５へ伝達する。
【００４１】
＜統合相関関数計算部＞
図２Ｃは、本実施形態に係る波源方向推定装置２００の備える統合相関関数計算２０４部の構成を示すブロック図である。統合相関関数計算部２０４は、周波数別相関関数生成部２４１１，２４１２，・・・，２４１Ｋと、統合部２４２を備える。
【００４２】
周波数別相関関数生成部２４１１，２４１２，・・・，２４１Ｋは、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋから供給される周波数別クロススペクトルの逆変換を行い、周波数別相関関数として統合部２４２にそれぞれ伝達する。本実施形態では、変換部２０１でフーリエ変換を用いたので、逆変換には逆フーリエ変換を用いる方法について説明する。周波数別クロススペクトル計算部２０３ｋから供給された周波数別クロススペクトルをＵｋ（ｗ，ｎ）とすると、Ｕｋ（ｗ，ｎ）の逆変換により得られる周波数別相関関数ｕｋ（τ，ｎ）は次のように計算される。
【００４３】
【数８】
2018003158.tif 000010
統合部２４２は、周波数別相関関数生成部２４１１，２４１２，・・・，２４１Ｋから供給される周波数別相関関数を統合し、統合相関関数として推定方向情報生成部２０５へ伝達する。個別に求めた複数の周波数別相関関数を混合したり、重ね合わせたりすることにより、一つの相関関数を求める。統合方法に単純な和を用いた場合、統合部２４２は、周波数別相関関数の総和を計算する。統合相関関数をｕ（τ，ｎ）とすると、ｕ（τ，ｎ）は次のように計算される。
【００４４】
【数９】
2018003158.tif 000011
また、総和では無く総乗を用いてもよい。この場合、ｕ（τ，ｎ）は次のように計算される。
【００４５】
【数１０】
2018003158.tif 000012
目標音の存在する周波数や、目的音のパワーが大きい周波数があらかじめ判明している場合には、その周波数に対応する周波数別相関関数だけ用いて、統合された相関関数を求めてもよい。また、重み付けという形で、統合における周波数別相関関数の影響度を制御してもよい。例えば目標音の存在する周波数の集合をΩとすると、周波数の選択によりｕ（τ，ｎ）を求める場合は、次のように計算される。
【００４６】
【数１１】
2018003158.tif 000013
また、重み付けを用いる場合には、ｕ（τ，ｎ）は次のように計算される。
【００４７】
【数１２】
2018003158.tif 000014
ここで、ａとｂとは実数で、ａ＞ｂ＞０を満足する。このように、目的音が存在する周波数の周波数別相関関数を重点的に用いて統合すると、雑音などの非目的音の影響が小さい相関関数を生成できるため、方向推定精度が向上する。
【００４８】
相対遅延時間計算部２０６は、入力されたマイク位置情報と音源探索対象方向とから、マイクペア間の相対遅延時間を求め、音源探索対象方向とセットで推定方向情報生成部２０５に伝達する。相対遅延時間とは、マイク間隔と音源方向とに基づいて一意に定まる音波の到達時間差のことである。音速をｃと仮定し、２つのマイクの間隔をｄ、音源方向、つまり音の到来方向をθとした場合、音源方向θに対する相対遅延時間τ(θ)は次の式で計算される。
【００４９】
【数１３】
2018003158.tif 000015
相対遅延時間は、全ての音源探索対象方向に対して計算される。例えば、方向の探索範囲が１０度刻みで０度から９０度まで、つまり０度、１０度、２０度、・・・、９０度である場合、１０種類の相対遅延時間を計算する。そして、探索対象の方向と相対遅延時間をペアで推定方向情報生成部２０５へ供給する。
【００５０】
推定方向情報生成部２０５は、統合相関関数計算部２０４から供給される統合相関関数と、相対遅延時間計算部２０６から供給される相対遅延時間に基づき、方向と相関値との対応関係を推定方向情報として出力する。相関関数をｕ（τ，ｎ）、相対遅延時間τ（θ）をとすると、推定方向情報Ｈ（θ，ｎ）は、次式で与えられる。
【００５１】
【数１４】
2018003158.tif 000016
各方向に対して相関値が定まるので、基本的に相関値が高ければ、その方向に音源が存在する可能性が高いと判断できる。
【００５２】
このような推定方向情報は、様々な形で利用される。例えば、関数が複数のピークを有する場合には、各ピークを到来方向とする音源が複数存在すると考えられる。したがって、各音源の方向を同時に推定できるだけでなく、音源数の推定にも用いることが可能である。
【００５３】
また、相関関数のピークと非ピークとの差分に基づき、音源の存在可能性を判定することも可能である。もし、ピークと非ピークとの差分が大きければ音源の存在可能性が高いと判定できる。同時に、推定方向の信頼性も高いと判断できる。もし、音源数が１つと事前に仮定できる場合には、相関値が最大となる方向を推定方向情報として出力してもよい。この場合、推定方向情報は、方向と相関値との対応関係ではなく、方向そのものとなる。
【００５４】
＜周波数別クロススペクトルの説明＞
以上の方法により周波数別クロススペクトルを計算すると、周波数別クロススペクトルを逆変換して得られる周波数別相関関数のピークが鋭くなり、相関関数のピーク位置が明確になる。相関関数のピーク位置に基づいて波源方向推定を行う本実施形態の場合には、ピークが鋭くなれば音源方向推定の精度が向上する。さらに、Ｐの値が大きくなるにつれて、すなわちｋを整数倍した周波数の成分が増えるにつれて、相関関数のピークは鋭くなる。図３Ａに、この様子を示す。ここで、図中のＱは３よりも大きい整数である。Ｐ＝１の場合、つまり位相成分が一つしかない場合、これを逆変換して得られる相関関数は、このようにピーク位置が不明確な相関関数になる。Ｐが大きくなると、図３Ａに示す通り、相関関数のピークは鋭くなる。
【００５５】
従って本発明では、周波数別クロススペクトルは、ある周波数kのクロススペクトルを基に「周波数kを整数倍した周波数ｐｋの位相成分に、周波数kの位相成分ａｒｇ（Ｓ１２（ｋ，ｎ））をｐ倍したものを割り当てたもの」として定義される。ここで、ｐは１以上の整数である。つまり、周波数別クロススペクトルは、その位相成分ａｒｇ（Ｕｋ（ｗ，ｎ））が少なくとも以下の式を満足するものとして定義される。
【００５６】
【数１５】
2018003158.tif 000017
加えてｐは、ｐ＝１と２、ｐ＝１と３、ｐ＝２と３など、２つ以上に限定される。ｐが１のみの場合には、周波数kの成分だけを抜き出して周波数別クロススペクトルが生成されることになるが、方向推定精度は従来技術と同等となり、方向推定の高精度を達成できない。なお、図３Ａで説明した通り、ｐ＝１，２，３，４・・・とｐの数が増えれば、つまり、周波数ｐｋの位相成分への割り当てが増えれば、周波数別相関関数のピークは鋭くなるため、方向推定の精度は向上する。
【００５７】
周波数別スペクトルを計算する際にPを大きくする効果について、図３Ｂを用いて説明する。図３Ｂは、異なるＰ＝１とＰ＝５との場合に対する、２つの周波数別相関関数と、それらを統合した統合相関関数の例を示した図である。図３Ｂの下段に示す通り、周波数別相関関数のピークは周期的に出現し、そのピーク間隔は周波数ｋに反比例する。周波数ｋが高くなれば、隣接する２つの周波数別相関関数のピークが接近し、相関関数の重なりからピークは目立たなくなる。Ｐ＝１の場合、周波数別相関関数のピークは鋭くないため、統合して得られた統合関数のピークも滑らかである。このため、推定方向は不明確となる。一方、Ｐ＝５の場合、周波数別相関関数、およびその統合相関関数のピークは鋭くなる。このため、推定方向は明確になり、推定精度の向上につながる。
【００５８】
図４は、本実施形態に係る波源方向推定装置２００の備える統合相関関数テーブル４０１の構成の一例を示す図である。統合相関関数テーブル４０１は、入力信号４１１に関連付けて周波数領域信号４１２と、クロススペクトル４１３と、周波数別クロススペクトル４１４と、統合相関関数４１５とを記憶する。波源方向推定装置２００は、入力信号を取得するたびに統合相関関数を計算してもよいし、入力信号に対応する統合相関関数をあらかじめ求めておき、統合相関関数テーブル４０１を参照して統合相関関数を算出してもよい。
【００５９】
図５は、本実施形態に係る波源方向推定装置２００のハードウェア構成を示すブロック図である。
【００６０】
ＣＰＵ（Central Processing Unit）５１０は演算制御用のプロセッサであり、プログラムを実行することで図２Ａの波源方向推定装置２００の機能構成部を実現する。ＲＯＭ（Read Only Memory）５２０は、初期データおよびプログラムなどの固定データおよびプログラムを記憶する。また、通信制御部５３０は、ネットワークを介してその他の装置などと通信する。なお、ＣＰＵ５１０は１つに限定されず、複数のＣＰＵであっても、あるいは画像処理用のＧＰＵ（Graphics Processing Unit）を含んでもよい。また、通信制御部５３０は、ＣＰＵ５１０とは独立したＣＰＵを有して、ＲＡＭ（Random Access Memory）５４０の領域に送受信データを書き込みあるいは読み出しするのが望ましい。また、ＲＡＭ５４０とストレージ５５０との間でデータを転送するＤＭＡＣ（Direct Memory Access Controller）を設けるのが望ましい（図示なし）。さらに、入出力インタフェース５６０は、ＣＰＵ５１０とは独立したＣＰＵを有して、ＲＡＭ５４０の領域に入出力データを書き込みあるいは読み出しするのが望ましい。したがって、ＣＰＵ５１０は、ＲＡＭ５４０にデータが受信あるいは転送されたことを認識してデータを処理する。また、ＣＰＵ５１０は、処理結果をＲＡＭ５４０に準備し、後の送信あるいは転送は通信制御部５３０やＤＭＡＣ、あるいは入出力インタフェース５６０に任せる。
【００６１】
ＲＡＭ５４０は、ＣＰＵ５１０が一時記憶のワークエリアとして使用するランダムアクセスメモリである。ＲＡＭ５４０には、本実施形態の実現に必要なデータを記憶する領域が確保されている。入力信号５４１は、マイクなどの集音装置で集音した音信号データや入力信号取得装置などに入力され、取得された信号データなどである。
【００６２】
周波数領域信号５４２は、入力信号５４１を変換部２０１で変換した信号である。クロススペクトル５４３は、クロススペクトル計算部２０２で計算したスペクトルである。周波数別クロススペクトル５４４は、周波数別クロススペクトル計算部２０３ｋで計算したスペクトルである。統合相関関数５４５は、統合相関関数計算部２０４で計算した関数である。
【００６３】
入出力データ５４６は、入出力インタフェース５６０を介して入出力されるデータである。送受信データ５４７は、ネットワークインタフェース５３０を介して送受信されるデータである。また、ＲＡＭ５４０は、各種アプリケーションモジュールを実行するためのアプリケーション実行領域５４８を有する。
【００６４】
ストレージ５５０には、データベースや各種のパラメータ、あるいは本実施形態の実現に必要な以下のデータまたはプログラムが記憶されている。ストレージ５５０は、統合相関関数テーブル４０１を格納する。統合相関関数テーブル４０１は、図４に示した、入力信号と統合相関関数との関係を管理するテーブルである。
【００６５】
ストレージ５５０は、さらに、変換モジュール５５１と、クロススペクトル計算モジュール５５２と、周波数別クロススペクトル計算モジュール５５３と、統合相関関数計算モジュール５５４とを格納する。また、ストレージ５５０は、推定方向情報生成モジュール５５５と、相対遅延時間計算モジュール５５６とを格納する。
【００６６】
変換モジュール５５１は、入力信号を周波数領域信号に変換するモジュールである。クロススペクトル計算モジュール５５２は、周波数領域信号に基づいてクロススペクトルを計算するモジュールである。周波数別クロススペクトル計算モジュール５５３は、クロススペクトルを用いて、周波数別クロススペクトルを計算するモジュールである。統合相関関数計算モジュール５５４は、周波数別クロススペクトルに基づいて、統合相関関数を計算するモジュールである。
【００６７】
推定方向情報生成モジュール５５５は、統合包絡関数に基づいて波源の推定方向情報を生成するモジュールである。相対遅延時間計算モジュール５５６は、相対遅延時間を計算するモジュールである。これらのモジュール５５１〜５５６は、ＣＰＵ５１０によりＲＡＭ５４０のアプリケーション実行領域５４８に読み出され、実行される。制御プログラム５５７は、波源方向推定装置２００の全体を制御するためのプログラムである。
【００６８】
入出力インタフェース５６０は、入出力機器との入出力データをインタフェースする。入出力インタフェース５６０には、表示部５６１、操作部５６２、が接続される。また、入出力インタフェース５６０には、さらに、記憶媒体５６４が接続されてもよい。さらに、音声出力部であるスピーカ５６３や、音声入力部であるマイク、あるいは、ＧＰＳ位置判定部が接続されてもよい。なお、図５に示したＲＡＭ５４０やストレージ５５０には、波源方向推定装置２００が有する汎用の機能や他の実現可能な機能に関するプログラムやデータは図示されていない。
【００６９】
図６本実施形態に係る波源方向推定装置２００の処理手順を示すフローチャートである。このフローチャートは、図５のＣＰＵ５１０がＲＡＭ５４０を使用して実行し、図２の波源方向推定装置２００の機能構成部を実現する。
【００７０】
ステップＳ６０１において、波源方向推定装置２００は、入力信号を取得する。ステップＳ６０３において、波源方向推定装置２００の変換部２０１は、入力端子２０１および入力端子２０２から供給された入力信号を変換する。変換部２０１は、変換して得られた周波数領域信号をクロススペクトル計算部２０２へ供給する。ステップＳ６０４において、クロススペクトル計算部２０２は、供給された変換信号に基づいて、クロススペクトルを計算する。クロススペクトル計算部２０２は、計算したクロススペクトルを周波数別クロススペクトル計算部２０３１，２０３ｋ，・・・２０３Ｋへ伝達する。
【００７１】
ステップＳ６０７において、周波数別クロススペクトル計算部２０３１，２０３ｋ，・・・２０３Ｋは、クロススペクトルの各周波数ｋに対応するクロススペクトルを計算する。すなわち、周波数別クロススペクトル計算部２０３１，２０３ｋ，・・・２０３Ｋは、周波数別クロススペクトルを計算する。そして、周波数別クロススペクトル計算部２０３１，２０３ｋ，・・・２０３Ｋは、周波数別クロススペクトルを統合相関関数計算部２０４へ伝達する。
【００７２】
ステップＳ６０９において、周波数別相関関数生成部２４１１，２４１２，・・・，２４１Ｋは、周波数別クロススペクトルを逆変換して、周波数別相関関数を計算する。ステップＳ６１１において、統合部２４２は、周波数別相関関数を統合して、統合相関関数を計算する。
【００７３】
ステップＳ６１３において、相対遅延時間計算部２０６は、マイク位置情報と音源探索対象方向とから、マイクペア間の相対遅延時間を計算する。ステップＳ６１５において、推定方向情報生成部２０５は、統合相関関数と相対遅延時間とから推定方向情報を生成する。
【００７４】
本実施形態によれば、入力信号に含まれる目標音の到来方向、すなわち目標物体が存在する方向を推定する。環境雑音レベルが高い環境において、目標物が発する音を手掛かりに、目標物が存在する方向を推定する場合に有効である。環境雑音の例としては、繁華街や街頭、街道沿い、人や自動車が多く集まる場所が挙げられる。また、目標物の例としては、人間や動物、自動車、航空機、船舶、水上バイク、ドローン（小型無人機）が挙げられる。
【００７５】
例えば、屋外のテーマパークや展示会場などに接近する不審な自動車・船舶・ドローンなどを検知し、それらの方向を推定することで、不審者や不審物を効率的に取り締ることが可能である。また、音源方向推定を複数箇所で実施することで、目標音源の位置を特定できる。これにより、環境雑音レベルが高い環境でも、悲鳴や銃声、自動車の衝突音の発生箇所などを正確に特定することが可能となる。
【００７６】
［第３実施形態］
次に本発明の第３実施形態に係る波源方向推定装置について、図７を用いて説明する。図７は、本実施形態に係る波源方向推定装置の備える統合相関関数生成部７０４の構成を説明するためのブロック図である。本実施形態に係る波源方向推定装置の備える統合相関関数生成部７０４は、上記第２実施形態の統合相関関数生成部２０４と比べると、周波数別相関関数生成部２４１１，２４１２，・・・，２４１Ｋと統合部２４２とに代えて、統合部７４１と統合相関関数生成部７４２とを有する点で異なる。その他の構成および動作は、第２実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
【００７７】
統合部７４１は、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋから供給される周波数別クロススペクトルを統合し、統合クロススペクトルとして統合相関関数生成部７４２へ伝達する。個別に求めた複数の周波数別クロススペクトルを混合したり、重ね合わせたりすることにより、一つの統合クロススペクトルを求める。統合には、第２実施形態の統合部２４２と同様に総和や総乗が用いられる。統合に総和を用いた場合、統合クロススペクトルＵ（ｋ，ｎ）は次のように計算される。
【００７８】
【数１６】
2018003158.tif 000018
また、総乗を用いた場合、統合クロススペクトルＵ（ｋ，ｎ）は次のように計算される。
【００７９】
【数１７】
2018003158.tif 000019
第２実施形態の統合部２４２と同様に、目標音源の存在する周波数や、目的音源のパワーが大きい周波数があらかじめ判明している場合には、統合クロススペクトルＵ（ｋ，ｎ）を生成する際に補正してもよい。第２実施形態と同様に、周波数の選択や重み付けという形で影響度を制御する。例えば目標音の存在する周波数の集合をΩとすると、帯域の選択により統合クロススペクトルＵ（ｋ，ｎ）を求める場合は、次のように計算する。
【００８０】
【数１８】
2018003158.tif 000020
また、重み付けを用いる場合には、Ｕ（ｋ，ｎ）は次のように計算される。
【００８１】
【数１９】
2018003158.tif 000021
ここで、ａとｂとは実数で、ａ＞ｂ＞０を満足する。このように、目的音が存在する周波数の周波数別相関関数を重点的に用いて統合すると、雑音などの非目的音の影響が小さい相関関数を生成できるため、方向推定精度が向上する。
【００８２】
統合相関関数生成部７４２は、統合部７４１から供給される統合クロススペクトルの逆変換を行い、統合相関関数として推定方向情報生成部２０５に伝達する。本実施形態でも、逆変換には逆フーリエ変換を用いる方法について説明する。統合部７４１から供給された統合クロススペクトルをＵ（ｋ，ｎ）とすると、Ｕ（ｋ，ｎ）の逆変換により得られる統合相関関数ｕ（τ，ｎ）は次のように計算される。
【００８３】
【数２０】
2018003158.tif 000022
本実施形態によれば、周波数別クロススペクトルを統合してから逆変換を行い、統合相関関数を得る。このため、周波数別クロススペクトルごとに逆変換を行っていた第２実施形態と比べると、逆変換の回数が少なくなる。したがって、第２実施形態よりも少ない計算量で統合相関関数を求めることができる。
【００８４】
［第４実施形態］
次に本発明の第４実施形態に係る波源方向推定装置について、図８Ａ乃至図１１を用いて説明する。図８Ａは、本実施形態に係る波源方向推定装置８００の構成を説明するためのブロック図である。本実施形態に係る波源方向推定装置８００は、上記第２実施形態と比べると、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋに代えて、周波数別クロススペクトル計算部８０３１，８０３２，・・・，８０３Ｋを有する点で異なる。その他の構成および動作は、第1実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
【００８５】
図８Ｂは、周波数別クロススペクトル計算部８０３ｋのブロック図である。周波数別クロススペクトル計算部８０３ｋは、周波数別基本クロススペクトル計算部２０３１ｋ、カーネル関数スペクトル記憶部８３１、乗算部８３２を有する。周波数別基本クロススペクトル計算部２０３１ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ１２（ｋ，ｎ）を用いて、Ｓ１２（ｋ，ｎ）の周波数ｋに対応するクロススペクトルを計算し、周波数別基本クロススペクトルとして乗算部８３２へ伝達する。周波数別基本クロススペクトル計算部２０３１ｋの動作は、その出力先を除いて第２実施形態の周波数別基本クロススペクトル計算部２０３１ｋと同様であるため、詳しい説明は省略する。
【００８６】
カーネル関数スペクトル記憶部８３１は、カーネル関数スペクトルを記憶しており、乗算部８３２へカーネル関数スペクトルを出力する。カーネル関数スペクトルとは、カーネル関数をフーリエ変換し、その絶対値を取ったものである。絶対値を取る代わりに二乗してもよい。カーネル関数としては、ガウス関数を用いる。ガウス関数は、あらかじめ与えられた３つの実数ｇ１，ｇ２，ｇ３を用いて、次のような数式で与えられる。
【００８７】
【数２１】
2018003158.tif 000023
ｇ１はガウス関数の大きさ、ｇ２はガウス関数のピークの位置、ｇ３はガウス関数の広がりを制御する。特に、ガウス関数の広がりを調整するｇ３は、周波数別相関関数のピークの鋭さに大きな影響を与えるので重要である。式（２１）から分かるように、ｇ３が大きくなればガウス関数の広がりは大きくなる。
【００８８】
他には、以下に示すロジスティック関数を用いてもよい。
【００８９】
【数２２】
2018003158.tif 000024
ここで、ｇ１とｇ２とは実数である。ロジスティック関数は、ガウス関数と同様の形状をしているが、ガウス関数よりも裾が長いという性質を持つ。特に、ロジスティック関数の広がりを調整するｇ５は、ガウス関数におけるｇ３の場合と同様、周波数別相関関数のピークの鋭さに大きな影響を与える重要なパラメータである。この他にも、コサイン関数や一様関数を用いてもよい。
【００９０】
カーネル関数で用いられるパラメータｇ１〜ｇ５は、定数ではなく、周波数ｋに応じて異なる値を用いてもよい。つまり、ｇ１（ｋ）〜ｇ５（ｋ）のように周波数kの関数としてもよい。例えば、ｇ３を周波数ｋの関数ｇ３（ｋ）とし、周波数が高くなるにつれて小さな値を持つ関数とする。このような代表例としてｋの逆数を関数ｇ３（ｋ）関数とする場合、ｇ３（ｋ）は次のように与えられる。
【００９１】
【数２３】
2018003158.tif 000025
ここで、Ｇ３は実数である。このとき、カーネル関数Ｇ（ｋ）は、周波数ｋが高くなるにつれて、ピークが鋭く裾が狭い関数になる。
【００９２】
乗算部８３２は、周波数別基本クロススペクトル計算部２０３１ｋから供給される周波数別基本クロススペクトルと、カーネル関数スペクトル記憶部８３１から供給されるカーネル関数スペクトルの積を計算し、周波数別クロススペクトルとして統合相関関数計算部２０４に伝達する。周波数別基本クロススペクトル計算部２０３１ｋから供給される周波数別基本スペクトルをＵｋ（ｗ，ｎ）、カーネル関数スペクトル記憶部８３１から供給されるカーネル関数スペクトルをＧ（ｗ）とすると、周波数別クロススペクトルＵＭｋ（ｗ，ｎ）は次のように計算される。
【００９３】
【数２４】
2018003158.tif 000026
このように、周波数別基本クロススペクトルにカーネル関数スペクトルを乗じることで、統合相関関数計算部２０４にある周波数別相関関数生成部２４１ｋで得られる周波数別相関関数の大きさを変更できる。図９は、カーネル関数スペクトルが乗じられた周波数別クロススペクトルと、周波数別相関関数の関係を示している。比較のため、カーネル関数スペクトルを乗じる前の周波数別クロススペクトルも掲載する。図９の左側図に示すように、カーネル関数スペクトルを乗じていなければ、高い周波数まで成分が存在するため、周波数別相関関数のピークは鋭くなる。一方で、図９の中央図、右側図のようにカーネル関数スペクトルを乗じると、高い周波数の成分が減衰するので、周波数別相関関数のピークの鋭さは小さくなる。つまり、カーネル関数スペクトルのピークが鋭くなる(カーネル関数スペクトルの裾が狭くなる)につれて、周波数別相関関数のピークの鋭さは小さくなる。また、図９の右側図のように、周波数別相関関数の裾が大きく広がると、隣接する山の裾が重なり、谷が浅い周波数別相関関数が得られる。
【００９４】
ここで、カーネル関数とカーネル関数スペクトルの形状の関係について補足する。フーリエ変換の性質により、形状の関係は逆になる。カーネル関数のピークが鋭く、裾が狭くなるにつれて、カーネル関数スペクトルのピークは平坦に近づき、裾が広くなる。ガウス関数の広がりを調整するｇ３との関係も含めて説明すると、ｇ３が大きくなれば、ガウス関数の広がりは大きくなるが、そのスペクトルの広がりは小さくなる。
【００９５】
カーネル関数による周波数別相関関数の大きさ制御の効果について、図１０において説明する。図１０は、カーネル関数の有無と統合相関関数の関係を示した図である。図１０の（ａ）に示すようなカーネル関数無しの場合、周波数別相関関数ｕ１（τ，ｎ）〜ｕ３（τ，ｎ）のピーク位置は近いが、ｕ１（τ，ｎ）〜ｕ３（τ，ｎ）の幅が細いため、統合時に大きなピークを形成できない。このため、ピークの位置が明確にならない。一方、図１０の（ｂ）に示すようなカーネル関数有りの場合、周波数別相関関数の幅は太くなっているため、統合によりｕ１（τ，ｎ）〜ｕ３（τ，ｎ）は大きなピークを形成できる。このため、（ａ）のカーネル関数無しの場合に比べて、ピークの位置が明確になる。
【００９６】
さらに、カーネル関数による周波数別相関関数の大きさ制御の別の効果について、図１１で説明する。図１１は、カーネル関数スペクトルの幅の違いと、統合相関関数との関係を示した図である。図９の右側図に示すように、幅の広いカーネル関数スペクトルを用いると、相関関数の周期性により、谷が浅い周波数別相関関数が形成される。したがって、図１１の（ｃ）に示す通り、谷が浅い周波数別相関関数を統合すると、谷が浅い、つまりピークの目立たない統合相関関数が生成される。一方、図９の中央図に示すように、幅の狭いカーネル関数スペクトルを用いた場合には、図９の右図よりも谷が深い周波数別相関関数が形成される。したがって、図１１の（ｄ）に示す通り、ピークが明確な統合相関関数が生成される。
【００９７】
本実施形態では、カーネル関数のフーリエ変換により得られたカーネル関数スペクトルと、周波数別基本クロススペクトルとの積を計算しているが、フーリエ変換の性質により時間領域でも実現できる。周波数別クロススペクトル計算部８０３Ｋでは無く、統合相関関数計算部２０４の中にある周波数別相関関数生成部２４１ｋの後段にカーネル関数を畳み込む「畳み込み演算部」を設け、周波数別相関関数生成部２４１ｋから供給される周波数別相関関数にカーネル関数を畳み込んでもよい。ただし、畳み込み演算は計算量が多いため、本実施形態のように周波数領域で積を計算するほうが効率的である。
【００９８】
本実施形態によれば、周波数別基本クロススペクトルにカーネル関数スペクトルを乗じて周波数別クロススペクトルを生成する。このため、逆変換により得られる周波数別相関関数の幅が太くなり、統合相関関数のピークが明確になる。特に、個別の周波数別相関関数のピーク位置が近いにも関わらず、それぞれ鋭いピークを持つ場合、補正を行うことにより統合相関関数のピークの明確化効果が高くなる。
【００９９】
［第５実施形態］
次に本発明の第５実施形態に係る波源方向推定装置について、図１２を用いて説明する。図１２は、本実施形態に係る波源方向推定装置の備える周波数別クロススペクトル計算部の構成を説明するための図である。本実施形態に係る波源方向推定装置の備える周波数別クロススペクトル計算部１２０３ｋは、上記第４実施形態と比べると、カーネル関数スペクトル記憶部８３１に代えて、カーネル関数スペクトル生成部１２３１を有する点で異なる。その他の構成および動作は、第４実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
【０１００】
カーネル関数スペクトル生成部１２３１は、クロススペクトル計算部２０２から供給されるクロススペクトルに用いてカーネル関数スペクトルを生成し、生成したカーネル関数スペクトルを乗算部８３２に伝達する。カーネル関数スペクトル生成部１２３１は、供給されたクロススペクトルを分析することで、入力信号に目標音が存在する可能性を求め、存在可能性を反映した形状を有するカーネル関数スペクトルを生成する。基本的に、存在可能性が低ければ、幅が狭い、広がりが小さいカーネル関数スペクトルを生成する。これにより、周波数別相関関数のピークは低くなるため、統合相関関数に誤ったピークが出現する可能性を低減できる。
【０１０１】
目標音の存在可能性を求める方法として、入力信号のＳＮＲ（目的信号対雑音比：Signal-to-Noise Ratio）を推定する方法について説明する。まず、供給されたクロススペクトルの絶対値を算出する。一般的には、変換部２０１で得られたフーリエ変換を二乗したものを入力信号パワースペクトルと呼ぶが、本実施形態ではクロススペクトルの絶対値を入力信号パワースペクトルとして扱う。次に、入力信号パワースペクトルに基づき、入力信号に含まれる雑音成分(非目的音成分)のパワースペクトルを推定する。入力信号パワースペクトルをＰＸ（ｋ，ｎ）とすると、ＰＸ（ｋ，ｎ）は以下のように計算される。
【０１０２】
【数２５】
2018003158.tif 000027
次に、入力信号パワースペクトルに基づいて、雑音成分のパワースペクトルを推定する。ここでは、非特許文献３に記載の方法を用いる。推定雑音パワースペクトルを、入力信号パワースペクトルが供給され始めた推定初期のパワースペクトルを平均化したものとする。この場合、推定が開始された直後には目的音が含まれないという条件を満たす必要がある。推定雑音パワースペクトルをＰＮ（ｋ，ｎ）とすると、ＰＮ（ｋ，ｎ）は以下のように計算される。
【０１０３】
【数２６】
2018003158.tif 000028
ここで、Ｎ０はあらかじめ定めた整数である。
【０１０４】
他の方法としては、推定雑音パワースペクトルを、入力信号パワースペクトルの最小値（最小統計量）から求める方法が非特許文献４に開示されている。この方法では、一定時間内における入力信号パワースペクトルの最小値を周波数ごとに保持し、その最小値から雑音成分を推定する。入力信号パワースペクトルの最小値は、雑音のパワースペクトルとスペクトル形状と似ているため、雑音パワースペクトルの推定値として用いることができる。
【０１０５】
推定雑音パワースペクトルが得られたあと、入力信号パワースペクトルとの比をとってＳＮＲの推定値とする。入力信号パワースペクトルをＰＸ（ｋ，ｎ）、推定雑音パワースペクトルをＰＮ（ｋ，ｎ）とすると、推定ＳＮＲγ（ｋ，ｎ）は、次のように計算される。
【０１０６】
【数２７】
2018003158.tif 000029
そして、この推定ＳＮＲであるγ（ｋ，ｎ）をそのまま目的音の存在可能性ｑ（ｋ，ｎ）に用いる。
【０１０７】
このようにして得られた推定ＳＮＲは、非特許文献３において、推定a-posteriori  ＳＮＲと呼ばれている。推定ＳＮＲには、推定a-posteriori  ＳＮＲの代わりに、非特許文献３に記載の方法で得られる推定a-priori  ＳＮＲを用いてもよい。a-priori  ＳＮＲの推定には、雑音成分の抑圧を行ってからＳＮＲを推定するので、計算量は多くなるもののa-posteriori  ＳＮＲよりも高い推定精度を達成できる。
【０１０８】
入力信号パワースペクトルと推定雑音パワースペクトルを用いて目標音の存在可能性を計算する方法は、推定ＳＮＲのような両者の比に限定されない。例えば、比の代わりに、両者の差分を用いてもよい。また、単純な大小関係を用いてもよい。
【０１０９】
クロススペクトルの分析により目標音が存在する可能性を求める方法としては、パワースペクトルを用いる方法に限定されない。他の代表的な例は、クロススペクトルの位相成分を分析する方法が挙げられる。位相成分の分析方法として、クロススペクトルの群遅延（位相成分を周波数方向に微分したもの）を用いた方法について説明する。はじめに、クロススペクトルの群遅延を求める。群遅延をｇｄ（ｋ，ｎ）とすると、クロススペクトルＳ１２（ｋ，ｎ）の群遅延は次のように算出できる。
【０１１０】
【数２８】
2018003158.tif 000030
そして、ｇｄ（ｋ，ｎ）の平均値を計算し、その平均値からの乖離度を存在可能性とする。例えば、ガウス関数を用いて目標音の存在可能性を計算する場合、存在可能性ｑ（ｋ，ｎ）は次のように計算できる。
【０１１１】
【数２９】
2018003158.tif 000031
ここで、ｑ０は正の実数である。また、ｇｄ（ｋ，ｎ）バーは、ｇｄ（ｋ，ｎ）を周波数方向に平均化したものである。平均化には様々な方法があるが、例えば、以下のような相加平均を用いることができる。
【０１１２】
【数３０】
2018003158.tif 000032
式（２９）を参照すると、ｇｄ（ｋ，ｎ）がｇｄ（ｋ，ｎ）バーに近ければ、ｑ（ｋ，ｎ）は１に近づく。一方、ｇｄ（ｋ，ｎ）がｇｄ（ｋ，ｎ）バーから遠ざかるにつれて、ｑ（ｋ，ｎ）は０に近づく。
【０１１３】
次に得られた存在可能性を用いて、カーネル関数スペクトルを生成する。ここでは、カーネル関数スペクトルの元となるカーネル関数のパラメータを制御する例について説明する。また、カーネル関数としては、ガウス関数を用いる例を説明する。目的音の存在可能性が高いときは、ｇ３を小さくする。これにより、存在可能性が高くなるほど、ｇ（τ）の幅は狭くなり、ｇ（τ）ピークが強調された形状に近づく。目的音の存在可能性からｇ３を求めるには、存在可能性の逆数を変数とする一次関数を用いる。このとき、存在可能性をｑ（ｋ，ｎ）とすると、ｇ３は次のように計算される。
【０１１４】
【数３１】
2018003158.tif 000033
ここで、ａ１とｂ１とは実数であり、ａ１＞０．０、ｂ１＞０．０を満足する。目的音の存在可能性ｑ（ｋ，ｎ）からｇ３を求める関数は、１次関数に限定されない。シグモイド関数や高次の多項式関数、非線形関数など、他の形で表される関数を1次関数の代わりに用いることも可能である。
【０１１５】
カーネル関数として、ロジスティック関数を用いる場合には、ｇ３と同様の方法でｇ５を計算すればよい。その結果、目的音の存在可能性が高ければ、ｇ５は小さくなるので、カーネル関数ｇ（τ）の幅は狭く、ピークが強調された形状に近づく。
【０１１６】
このように存在可能性からパラメータを生成した後、カーネル関数、およびカーネル関数スペクトルを生成する。
【０１１７】
本実施形態によれば、目的音の存在可能性を求め、その可能性に基づいてカーネル関数を算出する。可能性が高い場合には、カーネル関数スペクトルの幅は広がり、形状は平坦なものに近づく。逆に可能性が低い場合には、カーネル関数スペクトルの幅は狭くなる。これにより、目的音が存在する周波数の周波数別相関関数のピークが高くなり、目的音が存在しない周波数の周波数別相関関数のピークが低くなる。以上から、第４実施形態よりも統合相関関数のピークが強調され、目的音の方向推定精度が向上する。特に、非目的音の周波数別相関関数のピークが低くなるので、統合相関関数に誤ったピークが出現する可能性を低減できる。
【０１１８】
［第６実施形態］
次に本発明の第６実施形態に係る波源方向推定装置について、図１３Ａおよび図１３Ｂを用いて説明する。図１３Ａは、本実施形態に係る波源方向推定装置１３００の構成を説明するためのブロック図である。本実施形態に係る波源方向推定装置１３００は、上記第３実施形態と比べると、統合相関関数計算部２０４に代えて、統合相関関数計算部１３０４を有する点で異なる。その他の構成および動作は、第３実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
【０１１９】
図１３Ｂは、本発明の第６実施形態に係る波源方向推定装置の備える周波数別クロススペクトルの構成を説明するためのブロック図である。本実施形態に係る周波数別クロススペクトル計算部２０３ｋは、上記第３実施形態と比べると、統合部７４１に代えて、統合クロススペクトル生成部１３４１を有する点で異なる。その他の構成および動作は、第３実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
【０１２０】
統合クロススペクトル生成部１３４１は、クロススペクトル計算部２０２から供給されるクロススペクトルに基づき、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋから供給される周波数別クロススペクトルを統合し、統合クロススペクトルとして統合相関関数生成部７４２へ伝達する。第３実施形態では、目標音の存在する周波数や、目的音のパワーが大きい周波数があらかじめ判明している場合について説明した。本実施形態では、供給されたクロススペクトルを分析することで、入力信号に目標音が存在する可能性を求め、存在可能性に基づいて統合する。
【０１２１】
まず、供給されたクロススペクトルに基づき、目標音の存在可能性を求める。存在可能性の計算には、第５実施形態で説明した方法を同様に用いることができる。次に、求めた存在可能性を用いて、周波数クロススペクトルの統合を行う。まず、目標音の存在可能性をｑ（ｋ，ｎ）とすると、目標音の存在する可能性が高い周波数の集合Ωをｑ（ｋ，ｎ）に基づいて求める。ある周波数kに対するｑ（ｋ，ｎ）があらかじめ定めた閾値θｑを超えれば、その周波数を集合Ωの要素とする。これを数式で表現すると、次のようになる。
【０１２２】
【数３２】
2018003158.tif 000034
集合Ωが定まれば、第３実施形態で説明した方法を用いればよい。具体的には、式（１７）や式（１８）に示した計算式で求めればよい。
【０１２３】
また、存在可能性ｑ（ｋ，ｎ）を用いて重みを計算し、その重みを用いて重み付け和による統合を行ってもよい。重み関数をη（ｑ（ｋ，ｎ））とすると、統合クロススペクトルＵ（ｋ，ｎ）を次のように計算する。
【０１２４】
【数３３】
2018003158.tif 000035
ただし、重み関数η（ｑ（ｗ，ｎ））は、大きなｑ（ｗ，ｎ）に対して大きな値を取る単調増加関数とする。
【０１２５】
本実施形態によれば、クロススペクトルに基づき目標音の存在可能性を求めてから、その存在可能性を用いて統合クロススペクトルを計算する。このため、あらかじめ目標音の存在可能性が不明の状態でも、統合クロススペクトル生成時の帯域選択や重み付けが適切に行えるようになり、高い推定精度を達成できる。
【０１２６】
［第７実施形態］
次に本発明の第７実施形態に係る波源方向推定システムについて、図１４および図１５を用いて説明する。図１４は、本実施形態に係る波源方向推定システム１４００の構成を説明するための図である。本実施形態に係る波源方向推定システム１４００は、上記第２実施形態に係る波源方向推定装置２００を用いている。したがって、第２実施形態と同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
【０１２７】
本実施形態に係る波源方向推定システム１４００は、マイク１４０１、マイク１４０２、ＡＤ変換部１４０１および表示部１４０２を有する。なお、本実施形態では、波源方向推定装置２００の代わりに波源方向推定装置８００、または波源方向推定装置１３００を用いることが可能である。また、波源が音源であるという仮定で説明するため、マイクを用いた例を説明するが、音源以外の場合には、その波源から放射される波動を受信して電気信号に変換できる各種センサが、マイクの代わりに用いられる。
【０１２８】
マイク１４０１およびマイク１４０２は、推定対象である目標物体から生じる音を含めた装置周辺の音を電気信号に変換し、ＡＤ変換部１４０１へ伝達する。音が伝わる媒質が空気媒質である場合、音は空気の振動としてマイクに到達する。マイクは、到達した空気の振動を電気信号に変換する。
【０１２９】
ＡＤ変換部１４０１は、マイク１４０１、およびマイク１４０２から供給された音の電気信号をデジタル信号に変換し、入力端子２０１、および入力端子２０２に伝達する。
【０１３０】
表示部１４０２は、波源方向推定装置２００から供給された推定方向情報を画像などの可視化データに変換し、ディスプレイなどの表示装置に表示する。最も基本的な可視化方法は、ある時刻における相関関数を２次元グラフで表示する方法である。その際、横軸に方向、縦軸に相関値を表示する。ある時刻だけでなく、相関関数の時間変化を３次元で表示する方法も有効である。時間変化を表示することにより、目標音源の出現の明確化、目標音源の移動パタン、目標音源の移動方向の予測などが可能になる。３次元ではなく、２次元平面に投影する方法も有効である。３次元だと表示されたときに裏側が見づらいという問題がある。上から投影した平面上に表示すれば、死角が無くなり一覧性が向上する。相関値は、色の濃淡ではなく、等高線で表現してもよい。
【０１３１】
図１５は、本実施形態に係る波源方向推定システム１４００の表示部１４０２で表示した画像の一例を示す図であり、波源方向推定装置２００から供給された推定方向情報から得られたものである。これは、本実施形態の効果を確認する目的で取得した。例の作成には、街頭環境において、時刻２０秒〜２５秒に方位３０度で悲鳴が発生する状況の音を用いた。数センチ間隔で設置した二つのマイクを用いて集音した。
【０１３２】
図１５は、色が黒いほど相関値が高いことを表している。方位角の範囲は０から１８０度である。縦軸は時刻を表している。図１５を参照すると、時刻２０秒〜２５秒あたりで方位３０度の相関値が高くなっていることが分かる。このことから、悲鳴が時刻２０秒〜２５秒に発生していること、悲鳴の発生方向が約３０度であることが分かる。
【０１３３】
本実施形態によれば、推定方向情報を画像などの可視化データとして表示するので、ユーザが波源の方向推定情報を視覚的に把握することができる。
【０１３４】
［他の実施形態］
以上、実施形態を参照して本願発明を説明したが、本願発明は上記実施形態に限定されるものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。また、それぞれの実施形態に含まれる別々の特徴を如何様に組み合わせたシステムまたは装置も、本発明の範疇に含まれる。
【０１３５】
また、本発明は、複数の機器から構成されるシステムに適用されてもよいし、単体の装置に適用されてもよい。さらに、本発明は、実施形態の機能を実現する情報処理プログラムが、システムあるいは装置に直接あるいは遠隔から供給される場合にも適用可能である。したがって、本発明の機能をコンピュータで実現するために、コンピュータにインストールされるプログラム、あるいはそのプログラムを格納した媒体、そのプログラムをダウンロードさせるＷＷＷ(World Wide Web)サーバも、本発明の範疇に含まれる。特に、少なくとも、上述した実施形態に含まれる処理ステップをコンピュータに実行させるプログラムを格納した非一時的コンピュータ可読媒体（non-transitory computer readable medium）は本発明の範疇に含まれる。
【０１３６】
［実施形態の他の表現］
上記の実施形態の一部または全部は、以下の付記のようにも記載されうるが、以下には限られない。（付記１）
波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
前記入力信号取得手段で取得した複数の前記入力信号を複数の周波数領域信号に変換する変換手段と、
前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算手段と、
前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算手段と、
を備えた相関関数生成装置。（付記２）
前記統合相関関数計算手段は、
前記周波数別クロススペクトルの逆変換により周波数別相関関数を生成する周波数別相関関数生成手段と、
前記周波数別相関関数を統合して一つの統合相関関数を生成する統合相関関数生成手段と、
を有する付記１に記載の相関関数生成装置。（付記３）
前記統合相関関数計算手段は、
前記周波数別クロススペクトルを統合して統合クロススペクトルを生成する統合クロススペクトル生成手段と、
前記統合クロススペクトルの逆変換により統合相関関数を生成する統合相関関数生成手段と、
を有する付記１に記載の相関関数生成装置。（付記４）
前記周波数別クロススペクトル計算手段は、
前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段を有し、
前記周波数別基本クロススペクトルを前記周波数別クロススペクトルとして求める付記１乃至３のいずれか１項に記載の相関関数生成装置。（付記５）
前記周波数別クロススペクトル計算手段は、
前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
カーネル関数スペクトルを記憶するカーネル関数記憶手段と、
前記周波数別基本クロススペクトルと前記カーネル関数スペクトルとを乗じて、前記周波数別クロススペクトルを求める乗算手段と、
を有する付記１乃至３のいずれか１項に記載の相関関数生成装置。（付記６）
前記周波数別クロススペクトル計算手段は、
前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
前記クロススペクトルに基づいて、カーネル関数スペクトルを計算するカーネル関数スペクトル計算手段と、
前記周波数別基本クロススペクトルと前記カーネル関数スペクトルとを乗じて、前記周波数別クロススペクトルを求める乗算手段と、
を有する付記１乃至３のいずれか１項に記載の相関関数生成装置。（付記７）
波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
を含む相関関数生成方法。（付記８）
波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
をコンピュータに実行させる相関関数生成プログラム。（付記９）
付記１乃至６のいずれか１項に記載の相関関数生成装置と、
統合相関関数に基づいて、波源の推定方向情報を生成する推定方向情報生成手段と、
を備えた波源方向推定装置。
【０１３７】
この出願は、２０１６年６月２９日に出願された日本出願特願２０１６−１２８４８６を基礎とする優先権を主張し、その開示の全てをここに取り込む。
【図１】
2018003158.tif 000036
【図２Ａ】
2018003158.tif 000037
【図２Ｂ】
2018003158.tif 000038
【図２Ｃ】
2018003158.tif 000039
【図３Ａ】
2018003158.tif 000040
【図３Ｂ】
2018003158.tif 000041
【図４】
2018003158.tif 000042
【図５】
2018003158.tif 000043
【図６】
2018003158.tif 000044
【図７】
2018003158.tif 000045
【図８Ａ】
2018003158.tif 000046
【図８Ｂ】
2018003158.tif 000047
【図９】
2018003158.tif 000048
【図１０】
2018003158.tif 000049
【図１１】
2018003158.tif 000050
【図１２】
2018003158.tif 000051
【図１３Ａ】
2018003158.tif 000052
【図１３Ｂ】
2018003158.tif 000053
【図１４】
2018003158.tif 000054
【図１５】
2018003158.tif 000055
【国際調査報告】
2018003158.tif 000056
2018003158.tif 000057
2018003158.tif 000058
2018003158.tif 000059