JP2020043456A - 音響処理装置、音響処理方法およびプログラム - Google Patents

Abstract

【課題】音源を移動させなくてもマイクロホンアレイを校正することができる音響処理装置、音響処理方法およびプログラムを提供する。【解決手段】収音位置推定部は、Ｎ（Ｎは、３以上の整数）個の基準音源のそれぞれから収音素子までの基準音の到達時間である第１時間を前記収音素子が収音した収音信号を用いて定め、前記Ｎ個の基準音源のそれぞれから前記収音素子の位置の候補である候補位置までの前記基準音の到達時間である第２時間を取得し、前記第１時間と前記第２時間との差が小さいほど高い確率を与える確率関数を用いて、前記候補位置が収音素子の位置である可能性を示す尤度を算出し、前記尤度を最大化する候補位置を、前記収音素子の位置として定める。【選択図】図３

Description

本発明は、音響処理装置、音響処理方法およびプログラムに関する。

マイクロホンアレイは、複数のマイクロホンを備え、音源から到来した音を収音するとともに、音質の補正、雑音抑圧、音源方向の推定などの処理（以下、アレイ処理と呼ぶ）に用いられる。アレイ処理では、音源から個々のマイクロホンに到来するまでの音の伝達特性を示す伝達関数が用いられる。一般に、伝達関数は、個々のマイクロホンの位置に依存する。

しかしながら、マイクロホンの位置は常に一定とは限らない。例えば、特許文献１には、音源から個々のマイクロホンまでの伝達関数を測定用の音源を用いずに算出する音響処理装置について記載されている。この音響処理装置では、複数のマイクロホンが人型ロボットの頭部に固定されている。また、非特許文献１には、複数のマイクロホンが分散配置された細長のホースの形状を有する柔構造ロボットについて記載されている。いずれも複数のマイクロホンの位置が可変であり、マイクロホンアレイの校正において個々のマイクロホンの位置を知得することが重要となる。

特開２０１５−１５４２０７号公報

Ｙ．Ｂａｎｄｏ，Ｔ．Ｏｔｓｕｋａ，Ｋ．Ｉｔｏｙａｍａ ｅｔ ａｌ．"Ｐｏｓｔｕｒｅ ｅｓｔｉｍａｔｉｏｎ ｏｆ ｈｏｓｅ−ｓｈａｐｅｄ ｒｏｂｏｔ ｂｙ ｕｓｉｎｇ ａｃｔｉｖｅ ｍｉｃｒｏｐｈｏｎｅ ａｒｒａｙ"，Ａｄｖａｎｃｅｄ Ｒｏｂｏｔｉｃｓ，２０１５，Ｔａｙｌｏｒ ＆ Ｆｒａｎｃｉｓ，Ｖｏｌ．２９，Ｎｏ．１，ｐ．３５−４９，http://dx.doi.org/10.1080/01691864.2014.981291

マイクロホンアレイの校正では、例えば、次の措置がなされる。（ａ）各方向に設置された音源から特定の測定用の信号に基づく音を提示し、個々のマイクロホンで収音して得られる収音信号を記録したうえで解析する。（ｂ）非特許文献１に記載の柔構造ロボットでは、個々のマイクロホンの位置を予め定め、その位置にマイクロホンを設置もしくは調整する。（ｃ）特許文献１の音響処理装置では、伝達関数の算出のために動く音源を用いる。
しかしながら、（ａ）では、収音、記録の際に音源方向を与える必要がある。また、記録や解析のために多くの時間を要する。（ｂ）では、非特許文献１に記載の柔構造ロボットのように複雑な形状を有し、その配置が可変な環境では、予め定めた位置にマイクロホンを設置すること自体が困難である。（ｃ）では、音源が移動するとは限らないため、常に伝達関数を取得することができない。

本発明は上記の点に鑑みてなされたものであり、音源を移動させなくても確実にマイクロホンアレイを校正することができる音響処理装置、音響処理方法およびプログラムを提供することである。

（１）本発明は上記の課題を解決するためになされたものであり、本発明の一態様は、Ｎ（Ｎは、３以上の整数）個の基準音源のそれぞれから収音素子までの基準音の到達時間である第１時間を前記収音素子が収音した収音信号を用いて定め、前記Ｎ個の基準音源のそれぞれから前記収音素子の位置の候補である候補位置までの前記基準音の到達時間である第２時間を取得し、前記第１時間と前記第２時間との差が小さいほど高い確率を与える確率関数を用いて、前記候補位置が収音素子の位置である可能性を示す尤度を算出し、前記尤度を最大化する候補位置を、前記収音素子の位置として定める収音位置推定部を備える音響処理装置である。

（２）本発明の他の態様は、（１）の音響処理装置であって、前記収音位置推定部は、前記基準音を示す基準音信号と前記収音素子が収音した収音信号との相互相関関数に基づいて前記第１時間を算出する。

（３）本発明の他の態様は、（２）の音響処理装置であって、前記収音位置推定部は、離散時刻ごとに算出した前記相互相関関数を所定の連続関数で補間して得られる関数値を最大化する時刻に基づいて前記第１時間を算出する。

（４）本発明の他の態様は、（２）の音響処理装置であって、前記収音位置推定部は、前記相互相関関数の周波数領域における変換係数の位相成分を算出し、前記位相成分を時間領域に逆変換して得られる離散時刻ごとの逆変換関数値を所定の連続関数で補間して得られる関数値を最大化する時刻までの期間を、前記第１時間として算出する。

（５）本発明の他の態様は、（１）から（４）のいずれかの音響処理装置であって、前記基準音は、可聴帯域よりも高い周波数の成分からなる。

（６）本発明の他の態様は、（１）から（５）のいずれかの音響処理装置であって、所定の収音素子が収音した基準チャネルの収音信号と他の収音素子が収音した他チャネルの収音信号との到達時間差を要素とする到達時間差ベクトルを所定期間ごとに生成し、前記到達時間差ベクトルをクラスタごとに分類し、前記基準音源とは別個の音源に対応するクラスタを代表する到達時間差ベクトルに基づいて第１ステアリングベクトルを定める音源推定部と、前記第１ステアリングベクトルに基づく相関行列を固有値展開して固有ベクトルを算出し、音源方向に対応する応答係数を要素として有する第２ステアリングベクトルを生成し、前記第２ステアリングベクトルの各要素の二乗和を、２次以上の各次の前記固有ベクトルと前記第２ステアリングベクトルとの内積の総和で正規化した空間スペクトルを最大化する音源方向を定める音源方向推定部と、をさらに備える。

（７）本発明の他の態様は、音声処理装置における音声処理方法であって、音響処理装置における音響処理方法であって、Ｎ（Ｎは、３以上の整数）個の基準音源のそれぞれから収音素子までの基準音の到達時間である第１時間を前記収音素子が収音した収音信号を用いて定め、前記Ｎ個の基準音源のそれぞれから前記収音素子の位置の候補である候補位置までの前記基準音の到達時間である第２時間を取得し、前記第１時間と前記第２時間との差が小さいほど高い確率を与える確率関数を用いて、前記候補位置が収音素子の位置である可能性を示す尤度を算出し、前記尤度を最大化する候補位置を、前記収音素子の位置として定める収音位置推定ステップを有する音響処理方法である。

（８）本発明の他の態様は、音響処理装置のコンピュータに、Ｎ（Ｎは、３以上の整数）個の基準音源のそれぞれから収音素子までの基準音の到達時間である第１時間を前記収音素子が収音した収音信号を用いて定め、前記Ｎ個の基準音源のそれぞれから前記収音素子の位置の候補である候補位置までの前記基準音の到達時間である第２時間を取得し、前記第１時間と前記第２時間との差が小さいほど高い確率を与える確率関数を用いて、前記候補位置が収音素子の位置である可能性を示す尤度を算出し、前記尤度を最大化する候補位置を、前記収音素子の位置として定める収音位置推定手順を実行させるためのプログラムである。

上述の（１）、（７）または（８）の構成によれば、３個以上の基準音源のそれぞれから収音素子までの基準音の到達時間と、基準音源のそれぞれから候補位置までの到達時間により、それぞれの候補位置が収音素子の位置である確率を用いて尤度が算出される。そして、尤度が最大となる候補位置が収音素子の位置として定まる。そのため、３個以上の基準音源を用いて候補位置が収音素子の位置となる可能性が総合的に評価される。よって、到達時間に誤差が含まれていても、基準音源が移動しているか否かに関わらず、より確実に収音位置を定めることができる。

また、（２）の構成によれば、相互相関関数により基準音信号と収音信号との相関性を定量的に評価して、基準音の再生から収音までの第１時間を定めることができる。

また、（３）の構成によれば、離散時刻ごとに相互相関関数を定め、相互に隣接する２つの離散時刻間の時刻に対する相互相関関数を推定し、推定された相互相関関数が最大となる時刻が基準音の再生から収音までの第１時間として定められる。そのため、時間離散化による第１時間の推定誤差、ひいては収音位置の推定誤差を低減することができる。

また、（４）の構成によれば、相互相関関数の周波数領域における変換係数の振幅成分が除去されるので、相互相関関数の振幅の周波数依存性が解消される。そして、相互相関関数の振幅成分を除去して得られる時間領域の逆変換関数値が相互に隣接する２つの離散時刻間で補間される。そのため、時間離散化に加え、基準音信号の振幅の周波数特性による第１時間、ひいては収音位置の推定誤差を低減することができる。

また、（５）の特徴により、収音位置の推定に可聴帯域の音声よりも波長が短い超音波が用いられるので、第１時間の推定誤差、ひいては収音位置の推定誤差を低減することができる。また、基準音は人間に知覚されないため、人間に対する受聴環境が害されない。

また、（６）の構成によれば、Ｍチャネルの収音信号から個々の音源に対応する第１ステアリングベクトルが定まり、第１ステアリングベクトルに基づく相関行列に対するＭ個の固有ベクトルが算出される、そのうち、第１固有ベクトルは、基準音源とは別個の目的音源に対するステアリングベクトルが張られる信号部分空間の基底ベクトルとなり、第２固有ベクトルから第Ｍ固有ベクトルのいずれとも直交する。音源方向に対応する第２ステアリングベクトルと第２固有ベクトルないし第Ｍ固有ベクトルとの内積の値が最も小さくなるので、音源方向に対する値が最大となる空間スペクトルが算出される。そのため、音源数が未知であっても、位置が可変な収音素子で収音される収音信号を用いて、個々の音源の音源方向を正確に推定することができる。

本実施形態に係る音響処理システムの構成例を示す概略ブロック図である。 本実施形態に係るクラスタリング処理の例を示すフローチャートである。 本実施形態に係る収音位置推定処理の例を示すフローチャートである。 本実施形態に係る音源方向推定処理の例を示すフローチャートである。 本実施形態に係る音響処理システムの外観構成例を示す平面図である。 本実施形態に係る基準音源と収音素子の配置例を示す平面図である。 本実施形態による収音位置の推定結果の第１の例を示す表である。 本実施形態による収音位置の推定結果の第２の例を示す表である。 本実施形態による収音位置の推定結果の第３の例を示す表である。 本実施形態による収音位置の推定結果の第４の例を示す表である。 本実施形態に係る収音位置算出部により算出される尤度の分布例を示す平面図である。 本実施形態に係る音響処理システムと話者の配置例を示す平面図である。 本実施形態による音源数の推定結果の例を示す表である。 本実施形態による音源方向の推定結果の例を示す表である。

（第１の実施形態）
以下、図面を参照しながら本発明の第１の実施形態について説明する。
図１は、本実施形態に係る音響処理システム１の構成例を示す概略ブロック図である。
音響処理システム１は、音響処理装置１０、音源部２０および収音部３０を含んで構成される。音源部２０は、Ｎ個の基準音源２０−１〜２０−Ｎ（Ｎは、３以上の整数）を備える。収音部３０は、Ｍ個の収音素子３０−１〜３０−Ｍ（Ｍは、２以上の整数）を備える。

音響処理装置１０は、Ｎチャネルの基準音信号を生成し、生成した基準音信号をそれぞれのチャネルに対応する基準音源２０−１〜２０−Ｎに出力し、基準音を発生させる。
音響処理装置１０には、収音素子３０−１〜３０−ＭからＭチャネルの音響信号が入力される。音響処理装置１０は、出力した基準音信号と入力された音響信号に基づいて基準音源２０−１〜２０−Ｎのそれぞれから放射された時刻から収音素子３０−１〜３０−Ｍのそれぞれまで到達する基準音の到達時間を第１時間として算出する。他方、音響処理装置１０は、各基準音源について候補位置までの基準音の到達時間を第２時間として算出する。候補位置は、収音素子の位置の候補である。音響処理装置１０は、第１時間と第２時間との時間差が小さいほど高い確率を与える確率関数を用いて、候補位置が収音素子の位置である可能性を示す尤度を算出する。そして、音響処理装置１０は、尤度を最大化する候補位置を、収音素子の位置として定める。
なお、本実施形態では、「Ａを最大化するＢ」とは、ある条件のもとで、できるだけ大きいＡを与えるＢを意味し、Ａが絶対的に最大になることを意味するとは限らない。例えば、Ｂを求める過程において、より小さいＡを与えるＢが求まることや、異なる条件のもとでＡを最大化するＢとして異なるＢが求まることもありうる。

音源部２０は、基準音を発する音源である。基準音は、収音素子３０−１〜３０−Ｍのそれぞれの位置を定めるために用いる音である。基準音源２０−１〜２０−Ｎは、音響処理装置１０から入力される基準音信号に基づいて基準音を発する。基準音は、必ずしも可聴帯域（例えば、２０Ｈｚ−２０ｋＨｚ）の成分を含んでいなくてもよい。基準音は、可聴帯域よりも周波数が高い成分からなる超音波であってもよい。基準音源２０−１〜２０−Ｎは、例えば、超音波エミッタ、圧電アクチュエータなどの電気音響変換器である。基準音源２０−１〜２０−Ｎの位置は、それぞれ固定されていてもよい。

収音部３０は、収音素子３０−１〜３０−Ｍを含んで構成される。収音素子３０−１〜３０−Ｍは、それぞれ自部に到来する音を収音する。収音素子３０−１〜３０−Ｍは、例えば、ＭＥＭＳ（Ｍｉｃｒｏ−ｅｌｅｃｔｒｏ−ｍｅｃｈａｎｉｃａｌ Ｓｙｓｔｅｍｓ）マイクロホン、コンデンサマイクロホン、などのいずれでもよい。基準音源２０−１〜２０−Ｎが超音波を発する場合には、収音素子３０−１〜３０−Ｍは、可聴帯域の音波とより周波数が高い周波数帯域の超音波の両者を収音できればよい。収音素子３０−１〜３０−Ｍは、それぞれ収音した音波の波形を示す各チャネルの音響信号として収音信号を生成し、生成したＭチャネルの収音信号を音響処理装置１０に出力する。収音素子３０−１〜３０−Ｍの位置は、それぞれ可変である。収音部３０の構成例については、後述する。

次に、音響処理装置１０の構成例について説明する。音響処理装置１０は、入出力部１１と、収音位置推定部１２と、音源推定部１３と、音源方向推定部１４と、を含んで構成される。

入出力部１１には、収音位置推定部１２の基準音源制御部１２１からＮチャネルの基準音信号が入力される。入出力部１１は、入力されたＮチャネルの基準音信号を、それぞれのチャネルに対応する基準音源２０−１〜２０−Ｎに出力する。
入出力部１１には、Ｍチャネルの収音信号が収音素子３０−１〜３０−Ｍから入力される。入出力部１１は、入力されたＭチャネルの収音信号を収音位置推定部１２と音源推定部１３にそれぞれ出力する。入出力部１１は、例えば、入出力インタフェースである。

収音位置推定部１２は、基準音源制御部１２１と、到達時間算出部１２２と、候補位置情報記憶部１２３と、収音位置算出部１２４と、を含んで構成される。
基準音源制御部１２１は、基準音を再生させる基準音源（以下、再生音源と呼ぶ）を選択し、その再生タイミングを制御する。基準音として出力する音響信号は、特性が既知の音響信号であればよい。かかる音響信号として、例えば、スイープ正弦波のパルス（ｐｕｌｓｅ ｏｆ ｓｗｅｐｔ ｓｉｎｕｓｏｉｄ）、ガウシアンパルス、Ｍ系列などのいずれが用いられてもよい。基準音源制御部１２１は、例えば、所定の再生周期ごとに再生音源を巡回的に切り替える。基準音源制御部１２１は、切り替えにより選択される再生音源に対応する第ｎ（ｎは、１からＮまでのいずれかの整数）チャネルの基準音信号として所定の波形を示す音響信号を割り当て、その他の基準音源に対するチャネルには音響信号を割り当てない。基準音源制御部１２１は、Ｎチャネルの基準音信号を、入出力部１１を経由して基準音源２０−１〜２０−Ｎに出力するとともに到達時間算出部１２２に出力する。よって、再生音源である基準音源２０−ｎは再生音源として基準音を放射し、再生音源以外の基準音源は基準音を放射しない。

到達時間算出部１２２は、収音素子３０−１〜３０−Ｍから入力されるＭチャネルの収音信号と、基準音源制御部１２１から入力されるＮチャネルの基準音信号のうち第ｎチャネルの基準音信号を用いて、再生音源である基準音源２０−ｎが発する基準音の収音素子３０−１〜３０−Ｍのそれぞれに伝搬するまでの到達時間（ＴｏＡ：Ｔｉｍｅ ｏｆ Ａｒｒｉｖａｌ）τ_ｍｎ（ｍは、収音素子３０−ｍを示すインデックス、即ち、１からＭまでのいずれかの整数）を算出する。到達時間算出部１２２は、算出した到達時間τ_ｍｎを収音位置算出部１２４に出力する。

到達時間算出部１２２は、例えば、第ｍチャネルの収音信号ｘ_ｍ（ｔ）と第ｎチャネルの基準音信号ｅ_ｎ（ｔ）との相互相関関数を算出し、基準音の放射が開始される時刻から、相互相関関数が極大となる時刻ｔまでの期間を到達時間τ_ｍｎとして定める。
より具体的には、到達時間算出部１２２は、式（１）に示す関係を用いて到達時間τ_{ｍｎ［ｊ］}を定める。

但し、（１）に示す関係は、再生音源が再生周期ごとに巡回的に選択され、各再生周期に１つの観測期間が設定される場合を前提としている。式（１）において、ｊは、第ｊ再生周期を示すインデックスである。ａｒｇ ｍａｘ_{ｔｎ［ｊ］≦ｔ≦ｔｎ［ｊ］＋ｔｒ}（…）は、時刻ｔ_ｎ［ｊ］から時刻ｔ_{ｎ［ｊ］＋ｔｒ}までの観測期間のうち、…を最大とする時刻ｔを示す。ｔ_ｎ［ｊ］は、第ｊ再生周期の開始時の時刻を示す。基準音としてパルス（上記のスイープ正弦波のパルス、ガウシアンパルスが該当）が用いられる場合には、到達時間算出部１２２は、時刻ｔ_ｎ［ｊ］として、パルスの波形の立ち上がり時刻（オンセット時刻）を用いてもよい。また、到達時間算出部１２２は、この時刻ｔ_ｎ［ｊ］を基準時刻０としてもよい。ｔ_ｒは、最大到達時間を示す。最大到達時間ｔ_ｒは、測定可能とする到達時間の最大値に相当する。従って、最大到達時間ｔ_ｒを、測定する空間の大きさに応じて予め到達時間算出部１２２に設定しておけばよい。例えば、収音素子３０−１〜３０−Ｍがいずれも１個の第１の部材に設置され、基準音源２０−１〜２０−Ｎがいずれも１個の第２の部材に設置され、第２の部材が第１の部材に空間的に接している場合には、最大到達時間ｔ_ｒは、アレイ長と配置距離の総和に対応する音波の到達時間よりも大きければよい。アレイ長は、第１の部材において収音素子３０−１〜３０−Ｍが配置されている領域の大きさの特徴量である。配置距離は、第２の部材において基準音源２０−１〜２０−Ｎが配置されている領域の大きさの特徴量である。
なお、以下の説明では、基準音源２０−ｎから収音素子３０−ｍに基準音が到着する時点までの到達時間を第１時間と呼ぶ。到達時間算出部１２２は、算出した第１時間を示す第１時間情報を収音位置算出部１２４に出力する。

なお、式（１）に示す関係を用いて算出した到達時間τ_{ｍｎ［ｊ］}の精度は、収音信号ｘ_ｍ（ｔ）と基準音信号ｅ_ｎ（ｔ）のサンプリング周期よりも高くすることができない。サンプリング周期は、サンプリング周波数の逆数に相当し、相互相関はサンプリング周期ごとに離散化された時刻ｔごとに算出される。現実の到達時間が、サンプリング周期の整数倍とならなければ、算出される到達時間との誤差が生じる。
そこで、到達時間算出部１２２は、収音信号ｘ_ｍ（ｔ）と基準音信号ｅ_ｎ（ｔ）の相互相関を複数の離散化された時刻間で補間し、補間により得られる相互相関関数が極大となる時刻ｔを到達時間τ_{ｍｎ［ｊ］}として定める。

より具体的には、到達時間算出部１２２は、式（２）に示すようにＧＣＣ−ＰＨＡＴ（Ｇｅｎｅｒａｌｉｚｅｄ Ｃｒｏｓｓ−Ｃｏｒｒｅｌａｔｉｏｎ Ｐｈａｓｅ Ｔｒａｎｓｆｏｒｍ；一般化相互相関関数位相変換）に基づいて到達時間τ_{ｍｎ［ｊ］}を定めてもよい。

式（２）において、Ｅ_ｎ（ω）、Ｘ_ｍ（ω）は、それぞれ基準音信号ｅ_ｎ（ｔ）、収音信号ｘ_ｍ（ｔ）に対して短時間フーリエ変換（ＳＴＦＴ：Ｓｈｏｒｔ Ｔｉｍｅ Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ）を行って得られる周波数領域の変換係数を示す。…^*は、…の複素共役を示す。Ｆ^−１［…］は、…に対する短時間フーリエ変換の逆変換を示す。Ｔ［…］は、…に対する時間領域における補間を示す。補間において、後述の包絡関数が利用可能である。即ち、式（２）は、基準音信号ｅ_ｎ（ｔ）と収音信号ｘ_ｍ（ｔ）の相互相関関数の周波数領域の変換係数を、その絶対値で除算し、その変換係数の位相成分を算出するステップと、算出した位相成分に対して短時間フーリエ逆変換を行って時間領域における逆変換信号を算出するステップと、算出された逆変換信号を構成する離散化時刻ごとの逆変換信号値を補間して得られる補間値の最大値を到達時間τ_{ｍｎ［ｊ］}として定めるステップを示す。算出された逆変換信号は、相互相関関数の位相の変動として表れる振幅の変動成分を示す。位相成分は、絶対値が１に正規化された複素数であり、ｅ^ｊφ（φは、位相）と表される。従って、基準音の振幅の周波数依存性による影響が排除される。絶対値で除算して振幅による影響を除去して位相成分を取得することは、白色化（ｗｈｉｔｅｎｉｎｇ）とも呼ばれる。なお、逆変換信号値の補間において、後述の包絡関数を用いることが可能である。

また、到達時間算出部１２２は、式（３）に示すように離散化された時刻ごとの相互相関関数の絶対値を複数の時刻間で補間して得られる補間値が極大となる時刻ｔを到達時間τ_{ｍｎ［ｊ］}として定めてもよい。

式（３）において、ｅｎｖ（…）は、…の包絡関数（Ｅｎｖｅｌｏｐｅ Ｆｕｎｃｔｉｏｎ）を示す。包絡関数は、時刻ｔの連続関数であって、複数の離散化された時刻のそれぞれに対する関数値として相互相関関数の絶対値｜（ｅ_ｎ＊ｘ_ｍ）（ｔ）｜を補間するために用いられる。包絡関数として、例えば、時刻ｔに対するｐ次関数（ｐは、２以上の整数、例えば、８）を用いることができる。到達時間算出部１２２は、包絡関数として、例えば、時刻ｔ_ｎ［ｊ］を基準時とし、（ｔ−ｔ_ｎ［ｊ］）のｐ’乗（ｐ’は、０からｐまでの整数）をそれぞれ所定の係数ｃ_ｐで乗じて得られる乗算値の総和を算出する。

基準音として、超音波が用いられる場合には、到達時間算出部１２２は、相互相関関数を算出する前に、第ｍチャネルの収音信号ｘ_ｍ（ｔ）に対して高域通過フィルタリングを行って、可聴帯域よりも周波数が高い成分である高域成分を抽出してもよい。基準音信号ｅ_ｎ（ｔ）との相互相関関数を算出する際、到達時間算出部１２２は、抽出された高域成分を用いればよい。これにより、可聴帯域の音波が提示されている状況でも、式（１）〜（３）のいずれかに示す相互相関関数もしくはその位相成分への可聴帯域の成分による影響が排除される。

候補位置情報記憶部１２３には、予め候補位置情報を記憶させておく。候補位置情報は、収音素子３０−１〜３０−Ｍの位置の候補である候補位置を示す情報である。候補位置と情報として、例えば、グリッドマップ（Ｇｒｉｄ Ｍａｐ）が利用可能である。グリッドマップとは、二次元平面上に離散化して配置された複数の候補位置を示す情報である。グリッドマップでは、複数の候補位置が格子点上に配置される。格子点は、一般に互いに直交する２つの座標軸方向のそれぞれに等間隔に分布する。また、音響処理システム１が設置されるロボットの自律移動に用いられるＳＬＡＭ（Ｓｉｍｕｌｔａｎｅｏｕｓ Ｌｏｃａｌｉｚａｔｉｏｎ ａｎｄ Ｍａｐｐｉｎｇ；同時自己位置推定地図生成）方式を実行して生成されたグリッドマップが適用されてもよい。ＳＬＡＭは、レーダセンサ等の検出部を用いて、その時点における障害物までの方向と距離を検出し、候補位置ごとに障害物の存否もしくは存在確率を候補位置情報として定める手法である。

収音位置算出部１２４は、候補位置情報記憶部１２３に記憶された候補位置情報を参照し、基準音源２０−ｎから各候補位置に基準音が伝搬するまでの到達時間を算出する。以下の説明では、基準音源２０−ｎから各候補位置に基準音が伝搬する到達時間を第２時間と呼ぶ。なお、基準音源２０−ｎの位置が固定されている場合には、さらに基準音源２０−ｎごとの第２時間を含み、これらが候補位置ごとに対応付けて候補位置情報が構成されてもよい。収音位置算出部１２４は、逐次に第２時間を算出することなく、候補位置情報を参照して注目する候補位置に対応する第２時間を取得することができる。

収音位置算出部１２４は、到達時間算出部１２２から入力される第1時間情報が示す基準音源２０−ｎと収音素子３０−ｍの組ごとの第1時間と、基準音源２０−ｎと収音素子３０−ｍの候補位置の組ごとの第２時間に対する尤度を最大化する収音素子３０−ｍの位置を収音位置ξ_Ｍｍとして定める。収音位置算出部１２４は、この尤度を、基準音源２０−ｎの位置がその候補位置である可能性の度合いを示す指標として算出する。収音位置算出部１２４は、例えば、式（４）に示す関係を満たす収音位置ξ_Ｍｍを定める。

式（４）において、σは、予め定められた標準偏差を示す。τ_Ｇｎは、基準音源２０−ｎからいずれかの候補位置ξ_Ｇまで基準音が伝搬する第２時間を示す。つまり、基準音源２０−ｎならびに収音素子３０−ｍに対する第２時間τ_Ｇｎから第１時間τ_{ｍｎ［ｊ］}の差分の確率分布が、その差分に対するガウス関数（正規分布）となることが仮定されている。よって、式（４）は、収音位置算出部１２４が、収音素子３０−ｍの候補位置ξ_Ｇの尤度として、第２時間τ_Ｇｎから第１時間τ_{ｍｎ［ｊ］}の差分の確率分布の基準音源２０−ｎ間、ならびに観測対象とする再生周期ｊ間の積を算出するステップと、尤度を最大とする候補位置ξ_Ｇを収音素子３０−ｍの収音位置ξ_Ｍｍとして定めるステップを有することを示す。
収音位置算出部１２４は、収音素子３０−ｍそれぞれの収音位置ξ_Ｍｍを示す収音位置情報を音源推定部１３と音源方向推定部１４に出力する。

なお、式（４）において標準偏差σとして、例えば、基準音源２０−ｎの位置の測定誤差を音速ｃで除算して得られる誤差時間が用いられてもよい。式（４）は、第２時間τ_Ｇｎから第１時間τ_{ｍｎ［ｊ］}の差分の基準音源２０−ｎ間の二乗和が大きいほど、尤度が低くなることを示す。これにより、第２時間τ_Ｇｎから第１時間τ_{ｍｎ［ｊ］}の差分の大きさが基準音源２０−ｎ間で総合して評価でき、評価された差分が大きいほど低くなるように尤度が与えられる。

なお、第２時間τ_Ｇｎから第１時間τ_{ｍｎ［ｊ］}の差分の確率分布を示す確率密度関数は、ガウス関数に限られない。第２時間τ_Ｇｎから第１時間τ_{ｍｎ［ｊ］}の差分の大きさが基準音源２０−ｎ間で総合して評価でき、評価された差分が大きいほど低くなるように尤度が算出される関数であればよい。より具体的には、第２時間τ_Ｇｎから第１時間τ_{ｍｎ［ｊ］}の差分の絶対値が０のときに関数値が最大となり、その差分の絶対値が大きいほど０に漸近するように関数値が減少する連続関数であればよい。例えば、ガウス関数に代えて、ロジスティック分布、双曲線正割分布、ラプラス分布などが用いられてもよい。また、尤度を算出する際、収音位置算出部１２４は、必ずしも複数の再生周期ｊ間で確率分布の積を算出しなくてもよく、１つの再生周期で基準音源２０−ｎごとに算出した確率分布の基準音源２０−ｎの積を尤度として採用してもよい。

音源推定部１３は、到達時間差ベクトル生成部１３１と、音源特定部１３２と、を含んで構成される。
到達時間差ベクトル生成部１３１は、収音素子３０−１〜３０−Ｍから入出力部１１を経由して入力されるＭチャネルの収音信号を用いて、所定期間ごとに到達時間差ベクトルを生成する。

より具体的には、到達時間差ベクトル生成部１３１は、到達時間差（ＴＤＯＡ：Ｔｉｍｅ Ｄｉｆｆｅｒｅｎｃｅ ｏｆ Ａｒｒｉｖａｌ）をチャネルｍごとに算出する。ＴＤＯＡは、所定の収音素子（以下の例では、収音素子３０−１）に基準音源２０−１〜２０−Ｎとは別個の音源からの音波が到来する時刻と他の収音素子３０−ｍにその音源からの音波が到来する時刻との時間差である。基準音源２０−１〜２０−Ｎとは別個の音源は、アレイ処理の目的とする音源、例えば、発話音声、楽音、動作音などアレイ処理の目的とする音源である。
到達時間差ベクトル生成部１３１は、例えば、式（５）に示すようにＧＣＣ−ＰＨＡＴに基づいて第ｍチャネルにおける第ｆフレームのＴＤＯＡ τ_ｍ（ｆ）を算出することができる。

式（５）において、ω_Ｌ、ω_Ｈは、それぞれＴＤＯＡの算出に係る周波数帯域の最低周波数（例えば、２０〜１００Ｈｚ）、最高周波数（例えば、４ｋＨｚ〜２０ｋＨｚ）を示す。基準音として超音波が用いられる場合には、ω_Ｈは、基準音の周波数帯域の下限よりも低い周波数であってもよい。Ｘ_ｍ（ω，ｆ）は、フレームｆにおける第ｍチャネルの収音信号ｘ_ｍ（ｔ，ｆ）の周波数領域の変換係数を示す。即ち、式（５）は、フレームｆにおける第１チャネルの収音信号ｘ_１（ｔ，ｆ）と第ｍチャネルの収音信号ｘ_ｍ（ｔ，ｆ）との相互相関関数の周波数領域の変換係数Ｘ_１（ω，ｆ）Ｘ^＊ _ｍ（ω，ｆ）を、その変換係数の絶対値で正規化するステップと、正規化して得られる変換係数の位相成分に対して逆フーリエ変換するステップと、逆フーリエ変換により得られる白色化された時間領域の相互相関関数の値が最大となる時刻τをＴＤＯＡ τ_ｍ（ｆ）として定めるステップを含む。但し、τは、−Ｄ_ｍ／ｃからＤ_ｍ／ｃまでの範囲内の値に制限される。Ｄ_ｍ、ｃは、それぞれアレイ長、音速を示す。この制限により、ＴＤＯＡ τ_ｍ（ｆ）は、収音素子３０−１〜３０−ｍの実現可能な分布のもとで実現可能とする値の範囲内に制限される。
そして、到達時間差ベクトル生成部１３１は、第１チャネルを除く、第２〜第Ｍチャネルのそれぞれについて算出したＴＤＯＡ τ_ｍ（ｆ）を要素とするベクトル［τ_２（ｆ），…，τ_Ｍ（ｆ）］^Ｔを到達時間差ベクトルとして定める。よって、各１個の到達時間差ベクトルは、Ｍ−１次元のベクトルとなる。到達時間差ベクトル生成部１３１は、生成した到達時間差ベクトルを音源特定部１３２に出力する。

なお、式（５）において、各チャネルの変換係数Ｘ_１（ω，ｆ）、Ｘ_ｍ（ω，ｆ）に代えて、その時点までのＦフレーム間（つまり、第ｆ−Ｆ＋１フレームから第ｆフレーム）の時間平均値＜Ｘ_１（ω，ｆ）＞、＜Ｘ_ｍ（ω，ｆ）＞が用いられてもよい。
また、音源から到達する音のレベルが低い場合には、収音信号ｘ_ｍ（ｔ，ｆ）に含まれるノイズが相対的に顕著になる。そこで、到達時間差ベクトル生成部１３１は、収音部３０で収音される収音信号全体の強度の尺度として収音レベルＥ（ｆ）を算出し、収音レベルＥ（ｆ）が所定の収音レベルの閾値Ｔ_Ｅ以上であるフレームｆにおいて到達時間差ベクトルを生成すると判定し、収音レベルＥ（ｆ）か閾値Ｔ_Ｅ未満となるフレームｆにおいて到達時間差ベクトルを生成しないと判定してもよい。
到達時間差ベクトル生成部１３１は、例えば、式（６）に示すように、各チャネルの収音信号ｘ_ｍ（ｔ，ｆ）から収音レベルＥ（ｆ）を算出することができる。

式（６）において、Ｘ（ω，ｆ）は、各チャネルｍに対する変換係数Ｘ_ｍ（ω，ｆ）を要素とするＭ次元のベクトル［Ｘ_１（ω，ｆ），…，Ｘ_Ｍ（ω，ｆ）］^Ｔである。つまり、各チャネルｍに対する変換係数Ｘ_ｍ（ω，ｆ）の絶対値のチャネル間の二乗和の周波数ω_Ｌからω_Ｈまでの帯域間の平均値が収音レベルＥ（ｆ）として算出される。よって、収音レベルＥ（ｆ）が高いフレームからノイズの影響が少ない到達時間差ベクトルが取得される。

音源特定部１３２は、到達時間差ベクトル生成部１３１から入力される到達時間差ベクトルを収集する。音源特定部１３２は、十分に多数（少なくともＭより多い。典型的には、１００〜１０００個）の到達時間差ベクトルを分類（クラスタリング）する。音源特定部１３２は、クラスタリング手法として、予めクラスタ数を指定せずに実行可能な手法を用いる。クラスタリングにより得られるクラスタが音源に対応付けるため、音源数に相当するクラスタ数Ｃが未知である場合にも適用できるためである。

音源特定部１３２は、クラスタリング手法として、例えば、Ａｆｆｉｎｉｔｙ Ｐｒｏｐａｇａｔｉｏｎ（ＡＰ法；親和性伝搬法）を用いる。ＡＰ法は、到達時間差ベクトルｖ_ｉに対して初期値に依存せず一意にクラスタを定めることができる。ノイズや残響が多い環境や、同時発話が生じる場合など、音響環境が劣悪な場合においても、ノイズ、残響なども個々の音源とするクラスタが期待される。
ＡＰ法は、各クラスタを代表するｅｘａｍｐｌａｒ（代表ベクトル）を再帰的に求める手法である。ＡＰ法は、図２に示すように次のステップを有するクラスタリング手法である。

図２は、ＡＰ法に基づくクラスタリング手法の一例を示すフローチャートである。
（ステップＳ１０２）音源特定部１３２は、クラスタリング対象とする多数の到達時間差ベクトルのうち、各２つの到達時間ベクトルからなるデータペアごとに、２つの到達時間差ベクトルｖ_ｉ，ｖ_ｊ間の類似度ｓ（ｖ_ｉ，ｖ_ｊ）を算出する。類似度ｓ（ｖ_ｉ，ｖ_ｊ）として、到達時間差ベクトルｖ_ｉ，ｖ_ｊ間の負のユークリッド距離が適用可能である。負のユークリッド距離とは、通常のユークリッド距離に対して正負を反転して算出される。
なお、図２に示す例では、１つのデータペアを構成する２つの到達時間差ベクトルｖ_ｉ，ｖ_ｊが共通（つまり、ｉ＝ｊ）となることがある。その場合、類似度ｓ（ｖ_ｉ，ｖ_ｊ）は、最大値０となる。その後、ステップＳ１０４の処理に進む。

（ステップＳ１０４）音源特定部１３２は、データペアごとにｒｅｓｐｏｎｓｉｂｉｌｉｔｙ（適切性）ｒ（ｖ_ｉ，ｖ_ｊ）とａｖａｉｌａｂｉｌｉｔｙ（可用性）ａ（ｖ_ｉ，ｖ_ｊ）の初期値を、それぞれ０と設定する。
ｒｅｓｐｏｎｓｉｂｉｌｉｔｙ ｒ（ｖ_ｉ，ｖ_ｊ）は、到達時間差ベクトルｖ_ｊが到達時間差ベクトルｖ_ｉのｅｘａｍｐｌａｒとして適切な度合いを示す数値である。
ａｖａｉｌａｂｉｌｉｔｙ ａ（ｖ_ｉ，ｖ_ｊ）は、到達時間差ベクトルｖ_ｉが到達時間差ベクトルｖ_ｊのｅｘａｍｐｌａｒとして選択することの適切な度合いを示す数値である。その後、ステップＳ１０６の処理に進む。

（ステップＳ１０６）音源特定部１３２は、データペアごとに式（７）、（８）に示す関係を用いて、ｒｅｓｐｏｎｓｉｂｉｌｉｔｙ ｒ（ｖ_ｉ，ｖ_ｊ）とａｖａｉｌａｂｉｌｉｔｙ ａ（ｖ_ｉ，ｖ_ｊ）を算出する。

式（７）、（８）において、λは、ダンピングファクタを示す。ダンピングファクタλは、算出されるｒｅｓｐｏｎｓｉｂｉｌｉｔｙ ｒ（ｖ_ｉ，ｖ_ｊ）やａｖａｉｌａｂｉｌｉｔｙ ａ（ｖ_ｉ，ｖ_ｊ）の振動を防止もしくは緩和するためのパラメータである。振動は、ステップＳ１０６の処理の繰り返しにより発生することがある。ダンピングファクタλは、０より大きく、１より小さい所定の実数である。一般に、λの値が大きいほど振動が緩和されるが、ｒｅｓｐｏｎｓｉｂｉｌｉｔｙ ｒ（ｖ_ｉ，ｖ_ｊ）やａｖａｉｌａｂｉｌｉｔｙ ａ（ｖ_ｉ，ｖ_ｊ）の収束に要する計算回数が多くなる。
ρ（ｖ_ｉ，ｖ_ｊ）、α（ｖ_ｉ，ｖ_ｊ）は、それぞれｒｅｓｐｏｎｓｉｂｉｌｉｔｙ ｒ（ｖ_ｉ，ｖ_ｊ）、ａｖａｉｌａｂｉｌｉｔｙ ａ（ｖ_ｉ，ｖ_ｊ）の伝搬値を示す。

音源特定部１３２は、式（９）、（１０）に示す関係を用いて伝搬値ρ（ｖ_ｉ，ｖ_ｊ）、α（ｖ_ｉ，ｖ_ｊ）を算出することができる。その後、ステップＳ１０８の処理に進む。

（ステップＳ１０８）音源特定部１３２は、ｒｅｓｐｏｎｓｉｂｉｌｉｔｙ ｒ（ｖ_ｉ，ｖ_ｊ）とａｖａｉｌａｂｉｌｉｔｙ ａ（ｖ_ｉ，ｖ_ｊ）が収束したか否かを判定する。音源特定部１３２は、例えば、今回のｒｅｓｐｏｎｓｉｂｉｌｉｔｙ ｒ（ｖ_ｉ，ｖ_ｊ）と前回のｒｅｓｐｏｎｓｉｂｉｌｉｔｙ ｒ（ｖ_ｉ，ｖ_ｊ）との差分のデータペア間の二乗和を第１判定パラメータとして算出し、今回のａｖａｉｌａｂｉｌｉｔｙ ａ（ｖ_ｉ，ｖ_ｊ）と前回のａｖａｉｌａｂｉｌｉｔｙ ａ（ｖ_ｉ，ｖ_ｊ）との差分のデータペア間の二乗和を第２判定パラメータとして算出する。音源特定部１３２は、算出した第１判定パラメータが予め設定された第１判定パラメータの閾値よりも小さく、かつ、算出した第２判定パラメータが予め設定された第２判定パラメータの閾値よりも小さいとき、ｒｅｓｐｏｎｓｉｂｉｌｉｔｙ ｒ（ｖ_ｉ，ｖ_ｊ）とａｖａｉｌａｂｉｌｉｔｙ ａ（ｖ_ｉ，ｖ_ｊ）が収束したと判定する。音源特定部１３２は、算出した第１判定パラメータが予め設定された第１判定パラメータの閾値以上、または、算出した第２判定パラメータが予め設定された第２判定パラメータの閾値以上であるとき、ｒｅｓｐｏｎｓｉｂｉｌｉｔｙ ｒ（ｖ_ｉ，ｖ_ｊ）とａｖａｉｌａｂｉｌｉｔｙ ａ（ｖ_ｉ，ｖ_ｊ）が収束していないと判定する。音源特定部１３２は、収束していないと判定するとき、ステップＳ１０６の処理を繰り返し、収束したと判定するとき、ステップＳ１１０の処理に進む。

（ステップＳ１１０）音源特定部１３２は、到達時間差ベクトルｖ_ｉのｅｘａｍｐｌａｒとして、ｒｅｓｐｏｎｓｉｂｉｌｉｔｙ ｒ（ｖ_ｉ，ｖ_ｊ）とａｖａｉｌａｂｉｌｉｔｙ ａ（ｖ_ｉ，ｖ_ｊ）との和が最大となる到達時間差ベクトルｖ_ｊを選択する。音源特定部１３２は、ある１つのｅｘａｍｐｌａｒを共通とする到達時間差ベクトルｖ_ｉ同士が、そのｅｘａｍｐｌａｒを代表とする各１つのクラスタに属すると判定する。このようにして、音源特定部１３２は、クラスタごとにｅｘａｍｐｌａｒとそのクラスタに属する到達時間差ベクトルｖ_ｉを特定することができる。音源特定部１３２は、互いに異なるｅｘａｍｐｌａｒの数をクラスタ数Ｃとして計数し、各クラスタに属する到達時間差ベクトルｖ_ｉの数をクラスタサイズとして算出することができる。その後、図２に示す処理を終了する。

なお、上記ではクラスタリング手法がＡＰ法である場合を例にしたが、これには限られない。音源特定部１３２は、到達時間差ベクトルのクラスタリングにおいて、予めクラスタ数を定めずに実現できるクラスタリング手法を利用することができる。かかるクラスタリング手法として、例えば、階層型クラスタリングを利用することができる。階層型クラスタリングとは、１個もしくは複数の到達時間差ベクトルからなる各２つのクラスタ同士の距離を評価し、最も距離が小さい２つのクラスタ同士を１つのクラスタに統合する処理を順次実行する手法である。階層型クラスタリングにおけるクラスタ間の距離の評価方法として、ウォード法、群平均法、最短距離法、最長距離法など、いずれの手法も利用することができる。

音源特定部１３２は、生成したＣ個のクラスタのうち、クラスタサイズが所定のクラスタサイズの閾値以上となるクラスタを採用し、その閾値よりも少ないクラスタを棄却してもよい。
より具体的には、音源特定部１３２は、生成したＣ個のクラスタをクラスタサイズＮ_［ｉ］（１≦ｉ≦Ｃ）の昇順にソートし、クラスタサイズＮ_［ｉ］がその閾値Ｔ_Ｎ（例えば、８０〜２００）以上となるＣ’個のクラスタｉを採用し、閾値Ｔ_Ｎ未満となるクラスタｉを棄却する。そして、音源特定部１３２は、より順位が低い隣接クラスタｉ−１のクラスタサイズＮ_{［ｉ−１］}に対する注目クラスタｉのクラスタサイズＮ_［ｉ］との比を隣接クラスタサイズ比Ｎ’_［ｉ］として算出する。そして、音源特定部１３２は、隣接クラスタサイズ比Ｎ’_［ｉ］が所定の隣接クラスタサイズ比の閾値Ｔ_Ｒ（Ｔ_Ｒは、１よりも有意に大きい実数、例えば、１．３〜１．７）よりも高くなるクラスタｉ’を特定する。そして、音源特定部１３２は、第ｉ’クラスタ〜第Ｃクラスタを採用し、第１クラスタ〜第ｉ’−１クラスタを棄却する。音源特定部１３２は、採用されたクラスタの個数を音源数ＮＳＳ（Ｎｕｍｂｅｒ ｏｆ Ｓｏｕｎｄ Ｓｏｕｒｃｅ）をＣ’−ｉ’＋１個として推定することができる。これにより、ノイズ源、残響、同時発話などにより形成される小規模なクラスタが除去される。

音源特定部１３２は、採用したクラスタｉごとにｅｘａｍｐｌａｒとなる到達時間差ベクトルｖ_ｉに基づいてステアリングベクトルＡ_［ｉ］（ω）を生成する。ステアリングベクトルは、式（１１）に示すように第ｍチャネルと第１チャネルとの間のＴＤＯＡ τ’_［ｉ］ｍによる位相差に基づく応答係数ｅ^{ｊωτ’［ｉ］ｍ}を第ｍ次元の要素とするＭ次元のベクトルである。但し、第１次元の要素を１とする。

音源特定部１３２は、ステアリングベクトルＡ_［ｉ］（ω）に基づいて音源活性期間（ＳＳＡＰ：Ｓｏｕｎｄ Ｓｏｕｒｃｅ Ａｃｔｉｖｉｔｙ Ｐｅｒｉｏｄ）を定めてもよい。ＳＳＡＰは、有効な音源が提示されている期間である。音源特定部１３２は、所定の長さの観測区間であるウィンドウｗごとに、収音信号ｘ_ｍ（ｔ）と音源ｉに係るステアリングベクトルＡ_{［ｉ，ｗ］}（ω）を用いて各フレームｆについてＭＵＳＩＣ（Ｍｕｌｔｉｐｌｅ Ｓｉｇｎａｌ Ｃｌａｓｓｉｆｉｃａｔｉｏｎ）空間スペクトルＰ_{［ｉ，ｗ］}（ｆ）を算出する。但し、１つのウィンドウは、１個以上のフレームからなる。１つのウィンドウが複数のフレームからなる場合には、隣接するウィンドウに共通のフレームが含まれてもよい。音源特定部１３２は、ＭＵＳＩＣ空間スペクトルＰ_{［ｉ，ｗ］}（ｆ）（単に、空間スペクトルとも呼ばれる）を式（１２）に示す関係を用いて算出することができる。

式（１２）において、［…］^Ｈは、ベクトル又は行列…の共役転置を示す。ｅ_ｍ（ω）は、相関行列Ｒ（ω）の第ｍ固有ベクトルを示す。相関行列Ｒ（ω）は、第ｍ行第ｎ列の要素として、第ｍチャネルの収音信号ｘ_ｍ（ｔ）と第ｎチャネルの収音信号ｘ_ｎ（ｔ）の相互相関係数の周波数領域の変換係数を有する。即ち、音源特定部１３２は、収音信号ｘ_ｍ（ｔ）の周波数領域の変換係数Ｘ_ｍ（ω）を要素とする収音ベクトルＸ（ω）＝［Ｘ_１（ω），…，Ｘ_Ｍ（ω）］^Ｔから式（１３）に示す関係を用いて相関行列Ｒ（ω）を算出することができる。

式（１３）において、Ｅ［…］は、…の期待値を示す。
そして、音源特定部１３２は、その時点までのＦ_ｐ個のフレームの期間内のＭＵＳＩＣ空間スペクトルＰ_{［ｉ，ｗ］}（ｆ）のヒストグラムを、各音源ｉについて生成する。音源特定部１３２は、フレームごとのＭＵＳＩＣ空間スペクトルＰが、所定数（例えば、１０〜２０）の階級のうちいずれの階級の値であるかを判定し、ＭＵＳＩＣ空間スペクトルＰがそれぞれの階級内の値をとるフレームの数を度数として示すデータをヒストグラムとして定める。所定数の階級は、ＭＵＳＩＣ空間スペクトルＰとして取りうる値域を細分化した値の区間である。音源特定部１３２は、ＭＵＳＩＣ空間スペクトルＰが０となるフレームを計数対象のフレームから除外してもよい。

音源特定部１３２は、生成したヒストグラムを構成するフレームごとの階級についてクラスタリングを行って、閾値Ｔ_ＰよりもＭＵＳＩＣ空間スペクトルＰの値が大きい活性領域と、閾値Ｔ_ＰよりもＭＵＳＩＣ空間スペクトルＰの値が小さい不活性領域とに分類する。音源特定部１３２は、クラスタリングにおいて、生成されるクラスタ数をパラメータとして２に設定可能とする手法、例えば、ｋ−ｍｅａｎｓ法を用いる。音源特定部１３２は、活性領域に対する階級の最小値を閾値Ｔ_Ｐとして定める。音源特定部１３２は、活性領域に分類されたフレームをＳＳＡＰと判定し、それ以外のフレームをＳＳＡＰではないと判定することができる。なお、音源特定部１３２は、複数のＭＵＳＩＣ空間スペクトルＰを活性領域と不活性領域に分類できる手法としてクラスタリング以外の手法、例えば、サポートベクターマシンを用いてもよい。
そして、音源特定部１３２は、音源ｉと音源ｉ’のそれぞれについて、各フレームがＳＳＡＰであるか否かを示すＳＳＡＰパターンを形成する。

次に、音源特定部１３２は、音源ｉと音源ｉ’（但し、ｉ≠ｉ’）の組ごとにＳＳＡＰパターンの類似性を評価する。音源特定部１３２は、類似性の尺度として、例えば、Ｆ値（Ｆ−ｍｅａｓｕｒｅ）を用いることができる。Ｆ値ＦＭは、式（１４）に示すように適合率（Ｐｒｅｓｉｃｉｏｎ Ｒａｔｅ）ＰＲと再現率（Ｒｅｃａｌｌ Ｒａｔｅ）ＲＲとの調和平均として算出される。

再現率ＲＲは、音源ｉに係るＳＳＡＰのフレーム数に対する音源ｉ、ｉ’の両者についてＳＳＡＰと判定されたフレームの数の割合として定義される。適合率ＰＲは、音源ｉ’に係るＳＳＡＰのフレーム数に対する音源ｉ、ｉ’の両者についてＳＳＡＰと判定されたフレームの数の割合として定義される。Ｆ値ＦＭは０以上１以下の値をとりうる。Ｆ値ＦＭが０とは、音源ｉに係るＳＳＡＰパターンと音源ｉ’に係るＳＳＡＰパターンとが異なることを意味する。Ｆ値ＦＭが１とは、音源ｉに係るＳＳＡＰパターンと音源ｉ’に係るＳＳＡＰパターンが同様であることを示す。そこで、音源特定部１３２は、Ｆ値ＦＭが所定のＦ値の閾値Ｔ_ＦＭ（例えば、０．４〜０．６）よりも大きい場合、音源ｉに係るＳＳＡＰパターンと音源ｉ’に係るＳＳＡＰパターンが類似し、Ｆ値ＦＭが所定のＦ値の閾値Ｔ_ＦＭ以下であるとき、音源ｉに係るＳＳＡＰパターンと音源ｉ’に係るＳＳＡＰパターンが類似しないと判定することができる。

なお、類似性の尺度として、Ｆ値以外にも相関係数などが利用可能である。音源特定部１３２は、例えば、音源ごとに所定の評価区間内の各フレームについてＳＳＡＰで否かにより、１か０の値を要素として含むＳＳＡＰパターンベクトルを生成する。そして、音源特定部１３２は、ＳＳＡＰパターンベクトル間の内積を相関係数として算出する。音源特定部１３２は、算出した相関係数の値が所定の閾値よりも大きいか否かにより、対応するＳＳＡＰパターンが類似するか否かを判定することができる。

そして、音源特定部１３２は、音源ｉと音源ｉ’とでＳＳＡＰパターンが類似すると判定するとき、音源ｉと音源ｉ’のいずれか一方の音源とその一方の音源に係るステアリングベクトルを棄却し、残った音源とその残った音源に係るステアリングベクトルを採用する。音源特定部１３２は、例えば、音源ｉと音源ｉ’のそれぞれに対応するクラスタのクラスタサイズを比較し、クラスタサイズが大きい方のクラスタに対応する音源を採用すればよい。音源特定部１３２は、既に算出した音源数から棄却した音源の数を差し引く。
音源特定部１３２は、採用した音源ｉごとのステアリングベクトルＡ_［ｉ］（ω）を音源方向推定部１４に出力する。

一般に各１つの音源は１または複数の周波数成分を含む音を提示するが、複数の周波数成分の音が同時に提示されるとき、周波数成分間でＳＳＡＰパターンも互いに類似するはずである。他方、到達時間差ベクトルのクラスタリングを周波数ごとに実行すると、複数の周波数成分のそれぞれについて異なるクラスタが得られるため、共通の音源から発されたにも関わらず、それぞれ別個の音源として検出される可能性がある。上記のようにＳＳＡＰパターンを評価して、互いに類似する複数のＳＳＡＰパターンにそれぞれ対応する音源を共通の音源と判定し、いずれか一方の音源を棄却することで、過不足なく音源を検出することができる。

音源方向推定部１４は、収音位置算出部１２４から入力される収音位置ξ_Ｍｍと音源特定部１３２から入力されるステアリングベクトルＡ_［ｉ］（ω）に基づいて、音源ｉの方向θ_［ｉ］を推定する。
より具体的には、音源方向推定部１４は、音源ｉに係るステアリングベクトルＡ_［ｉ］（ω）の共役転置Ａ^Ｈ _［ｉ］（ω）にステアリングベクトルＡ_［ｉ］（ω）を乗算して相関行列Ｒ_[ｉ]（ω）を算出する。音源方向推定部１４は、相関行列Ｒ_[ｉ]（ω）に対して固有値展開を行って、Ｍ個の固有値と、それぞれに対応する固有ベクトルを算出する。音源方向推定部１４は、固有ベクトルｖ_[ｉ]ｍ（ω）の順序ｍは、対応する固有値λ_ｍの絶対値の降順に定める。

音源方向推定部１４は、水平面内の方向θのそれぞれについて収音位置ξ_Ｍｍに対するステアリングベクトルＡ_ξ（ω，θ）を生成する。音源方向推定部１４は、ステアリングベクトルＡ_ξ（ω，θ）の第ｍ次元の要素値である応答係数として、収音位置ξ_Ｍｍ、ξ_Ｍ１間のＴＤＯＡ τ（ξ_Ｍｍ，ξ_Ｍ１，θ）に基づく位相成分ｅ^{ｊω（τ（ξＭｍ，ξＭ１，θ））}を定める。但し、ＴＤＯＡを算出する際、収音素子３０−１〜３０−Ｍに方向θから平面波が到来することを仮定する。この仮定では、収音位置ξ_Ｍｍ、ξ_Ｍ１間のＴＤＯＡ τ（ξ_Ｍｍ，ξ_Ｍ１，θ）は、｛（ξ_Ｍｍ，ｘ−ξ_Ｍ１，ｘ）ｃｏｓθ＋（ξ_Ｍｍ，ｙ−ξ_Ｍ１，ｙ）ｓｉｎθ｝／ｃとなる。ここで、ξ_Ｍｍ，ｘ、ξ_Ｍｍ，ｙは、それぞれ収音位置ξ_Ｍｍのξ_ｘ座標、ξ_ｙ座標を示す。ξ_Ｍ１，ｘ、ξ_Ｍ１，ｙは、それぞれ収音位置ξ_Ｍ１のξ_ｘ座標、ξ_ｙ座標を示す。ξ_ｘ座標、ξ_ｙ座標は、２次元直交座標系を構成する互いに直交するξ_ｘ、ξ_ｙ方向の座標値である。また、方向θは、ξ_ｘ方向を基準方向（０°）としてなす右回りの角度がθとなる方向である。
そして、音源方向推定部１４は、ステアリングベクトルＡ_ξ（ω，θ）と第２固有ベクトルｖ_[ｉ]２（ω）から第Ｍ固有ベクトルｖ_[ｉ]Ｍ（ω）を用いて、式（１５）に示す関係を用いて音源ｉの音源方向θ_[ｉ]を定めることができる。

式（１５）の右辺のカッコ内で表される部分は、ＭＵＳＩＣ空間スペクトルを示す。式（１５）は、ステアリングベクトルＡ_ξ（ω，θ）の各要素の二乗和を、ステアリングベクトルＡ_ξ（ω，θ）と第ｍ固有ベクトルｖ_[ｉ]ｍ（ω）との内積の絶対値の第２固有ベクトルｖ_[ｉ]２（ω）から第Ｍ固有ベクトルｖ_[ｉ]Ｍ（ω）までの総和で正規化するステップと、正規化して得られる値をさらに周波数ω間で累積してＭＵＳＩＣ空間スペクトルを算出するステップと、ＭＵＳＩＣ空間スペクトルが最大となるθを音源ｉの方向θ_［ｉ］として定めるステップを示す。

この手法において、第１固有ベクトルｖ_[ｉ]Ｍ（ω）は、音源ｉに対するステアリングベクトルが張られる信号部分空間の基底ベクトルとなり、第２固有ベクトルｖ_[ｉ]２（ω）から第Ｍ固有ベクトルｖ_[ｉ]Ｍ（ω）のいずれとも直交する。方向θが音源ｉの方向θ_［ｉ］となるとき、ステアリングベクトルＡ_ξ（ω，θ）と第ｍ固有ベクトルｖ_[ｉ]ｍ（ω）（ｍは、２以上）との内積の値が最も小さくなるので、ＭＵＳＩＣ空間スペクトルは、方向θが音源ｉの方向θ_［ｉ］となるとき最大となることが期待される。これにより、推定された音源ごとに、位置が可変な収音素子３０−ｍで収音される収音信号ｘ_ｍ（ｔ）を用いて音源方向が推定される。

なお、音源方向推定部１４は、音源方向θ_[ｉ]に代えて、音源ｉに対する音源位置ξ_[ｉ]を定めてもよい。但し、音源方向推定部１４は、候補位置ξ_Ｇのそれぞれについて収音位置ξ_Ｍｍに対するステアリングベクトルＡ_ξ（ω，ξ_Ｇ）を生成する。音源方向推定部１４は、ステアリングベクトルＡ_ξ（ω，ξ_Ｇ）の第ｍ次元の要素値として収音位置ξ_Ｍｍ、ξ_Ｍ１間のＴＤＯＡに基づく位相成分ｅ^{ｊω（τ（ξＭｍ，ξＧ）−τ（ξＭ１，ξＧ））}に定める際、音源位置ξ_ｉから収音位置ξ_Ｍｍに球面波が到来することを仮定してもよい。即ち、音源方向推定部１４は、音源位置ξ_Ｇから収音素子３０−ｍの収音位置ξ_Ｍｍまでの距離ｄ_ｍを音速ｃで除算して到来時刻τ（ξ_Ｍｍ，ξ_Ｇ）を算出する。また、第ｍ次元の要素値を定める際、音源位置ξ_Ｇから収音位置ξ_Ｍｍへの距離減衰をさらに考慮してもよい。即ち、音源方向推定部１４は、ステアリングベクトルＡ_ξ（ω，ξ_Ｇ）の第ｍ次元の要素値を、上記の位相成分ｅ^{ｊω（τ（ξＭｍ，ξＧ）−τ（ξＭ１，ξＧ））}にさらに距離ｄ_１に対する距離ｄ_ｍの比ｄ_ｍ／ｄ_１で除算して得られる値に定める。
そして、音源方向推定部１４は、ステアリングベクトルＡ_ξ（ω，ξ_Ｇ）と第２固有ベクトルｖ_[ｉ]２（ω）から第Ｍ固有ベクトルｖ_[ｉ]Ｍ（ω）を用いて、式（１６）に示す関係を用いて音源ｉの音源位置ξ_[ｉ]を定めることができる。

式（１６）の右辺のカッコ内で表される部分は、ＭＵＳＩＣ空間スペクトルを示す。式（１６）は、ステアリングベクトルＡ_ξ（ω，ξ_Ｇ）の各要素の二乗和を、ステアリングベクトルＡ_ξ（ω，ξ_Ｇ）と第ｍ固有ベクトルｖ_[ｉ]ｍ（ω）との内積の絶対値の第２固有ベクトルｖ_[ｉ]２（ω）から第Ｍ固有ベクトルｖ_[ｉ]Ｍ（ω）までの総和で正規化するステップと、正規化して得られる値をさらに周波数ω間で累積してＭＵＳＩＣ空間スペクトルを算出するステップと、ＭＵＳＩＣ空間スペクトルが最大となるξ_Ｇを音源ｉの位置ξ_［ｉ］として定めるステップを示す。
この手法により、候補位置ξ_Ｇに仮想的に設置された音源からの音波の拡散による到来方向や音量の収音素子３０−１〜３０−ｍ間の差異を考慮してステアリングベクトルＡ_ξ（ω，ξ_Ｇ）が算出される。そのため、音源が収音素子３０−１〜３０−ｍに比較的近接している環境下であっても、その音源位置を正確に特定することができる。

（収音位置推定処理）
次に、本実施形態に係る収音位置推定処理について説明する。図３は、本実施形態に係る収音位置推定処理の例を示すフローチャートである。
（ステップＳ１１２）到達時間算出部１２２は、到達時間算出部１２２は、第ｍチャネルの収音信号ｘ_ｍ（ｔ）と第ｎチャネルの基準音信号ｅ_ｎ（ｔ）との相互相関関数を算出する。その後、ステップＳ１１４の処理に進む。
（ステップＳ１１４）到達時間算出部１２２は、式（１）〜（３）のいずれかに示す関係に従って、相互相関関数もしくはその時間領域の位相成分が極大となる時刻ｔを算出する。到達時間算出部１２２は、基準音源２０−ｎから基準音が放射される時刻から算出した時刻までの期間を到達時間τ_{ｍｎ［ｊ］}（第１時間）として定める。その後、ステップＳ１１６の処理に進む。

（ステップＳ１１６）収音位置算出部１２４は、第１時間である到達時間τ_{ｍｎ［ｊ］}と候補位置情報が示す候補位置ξ_Ｇのそれぞれに対応する第２時間である到達時間τ_Ｇに基づいて、その候補位置ξ_Ｇが収音位置ξ_Ｍｍとなる可能性を示す尤度を式（４）に基づいて算出する。その後、ステップＳ１１８の処理に進む。
（ステップＳ１１８）収音位置算出部１２４は、算出した尤度を最大化する候補位置ξ_Ｇを収音位置ξ_Ｍｍとして定める。その後、図３に示す処理を終了する。

（音源方向推定処理）
次に、本実施形態に係る音源方向推定処理について説明する。図４は、本実施形態に係る音源方向推定処理の例を示すフローチャートである。
（ステップＳ１２２）収音位置算出部１２４は、図３に示す処理を実行して収音位置ξ_Ｍｍを推定する。その後、ステップＳ１２４の処理に進む。
（ステップＳ１２４）到達時間差ベクトル生成部１３１は、Ｍチャネルの収音信号ｘ_１〜ｘ_Ｍを用いて、所定の期間ごとに第ｍチャネルと第１チャネルとのＴＤＯＡを要素とする到達時間差ベクトルを生成する。その後、ステップＳ１２６の処理に進む。

（ステップＳ１２６）音源特定部１３２は、到達時間差ベクトルに対してクラスタリング処理を行い、クラスタとクラスタごとの代表ベクトル（ｅｘａｍｐｌａｒ）を定める。その後、ステップＳ１２８の処理に進む。
（ステップＳ１２８）音源特定部１３２は、定めたクラスタの全部または一部のそれぞれを音源に対応するクラスタとして定める。音源特定部１３２は、定めたクラスタの代表ベクトルの要素に基づいて応答係数を算出し、算出した応答係数を要素として有するステアリングベクトルを式（１１）に従って音源に対応するステアリングベクトルとして生成する。その後、ステップＳ１３０の処理に進む。

（ステップＳ１３０）音源方向推定部１４は、ステップＳ１２８で生成されたステアリングベクトルに基づいて相関行列を算出し、算出した相関行列に対して固有値展開を行って、第１固有ベクトルから第Ｍ固有ベクトルを算出する。音源方向推定部１４は、音源方向に対応する別個のステアリングベクトルを生成し、別個に生成したステアリングベクトルと第１固有ベクトルから第Ｍ固有ベクトルを用いて、式（１５）に示す空間スペクトルを算出する。その後、ステップＳ１３２の処理に進む。
（ステップＳ１３２）音源方向推定部１４は、算出した空間スペクトルを最大化する音源方向を、ステップＳ１２８で生成されたステアリングベクトルに対応する音源の音源方向として特定する。その後、図４に示す処理を終了する。

（外観構成）
次に、本実施形態に係る音響処理システム１の外観構成例について説明する。図５は、本実施形態に係る音響処理システム１の外観構成例を示す平面図である。
図５は、柔構造ロボットに設置されている音響処理システム１を例示する。
音響処理装置１０は、円盤状の筐体内に格納されている。基準音源の個数、収音素子の個数は、それぞれ３、５である。基準音源２０−１〜２０−３は、それぞれ筐体表面の外周部に設置されている。収音素子３０−１〜３０−５は、それぞれ柔構造ロボットの管部（ホース）の長手方向にほぼ等間隔に設置されている。管部の長さ、筐体の直径が、それぞれアレイ長、配置距離に相当する。柔構造ロボットの一端は、音響処理装置１０を格納する筐体に接続されている。柔構造ロボットは、管部の形状を変えながら平面上を爬行可能としている。そのため、収音素子３０−１〜３０−５全体の位置も、収音素子３０−１〜３０−５の相互間の相対的な位置関係も可変である。

（評価実験）
上記の実施形態で実行される音響処理方法の有効性を検証するため、出願人は、次に説明する評価実験を行った。評価実験を行った実験室の内面は、グラスウォールで覆われ、残響時間の尺度であるＲＴ２０が０．４〜０．５［ｓ］である。縦、横の寸法は、それぞれ４［ｍ］、７．５［ｍ］である。その環境における、平均Ｓ／Ｎ比は、１０−１５［ｄＢ］である。
評価実験は、図５に示す外観構成を有する音響処理システム１を用い、図６に示す基準音源２０−１〜２０−３と収音素子３０−１〜３０−５の配置のもとで行われた。基準音源２０−１〜２０−３は、原点Ｏを中心とする半径ｄ_Ｅの円周上に等角度（１２０°）間隔で配置された。半径ｄ_Ｅは、０．０４［ｍ］である。収音素子３０−１、３０−２、３０−３は、原点からの距離を特性長ｄ_Ｍとし、それぞれ２次元直交座標系のξ_ｘ方向（横方向）とは逆方向に、原点からξ_ｙ方向（縦方向）に、原点からξ_ｙ方向とは逆方向に離れた位置に配置された。収音素子３０−４、３０−５は、それぞれ座標［ｄ_Ｍ，ｄ_Ｍ／２］、座標［ｄ_Ｍ，−ｄ_Ｍ／２］の位置に配置された。但し、収音部３０の特性長ｄ_Ｍを、０．０４［ｍ］、０．０８［ｍ］、０．１２［ｍ］の３通りとした。

基準音信号、収音信号のサンプリング周波数をそれぞれ１９２［ｋＨｚ］とし、サンプルごとの信号値のビットサイズを１６ビットとした。但し、音源推定部１３、音源方向推定部１４は、収音信号のサンプリング周波数を１６［ｋＨｚ］にダウンサンプリングして、それぞれの処理を実行した。また、フレーム長、シフト長を、それぞれ５１２［サンプル］、１６０［サンプル］とした。

評価実験は、第１の実験と、第２の実験とを含む。第１の実験では、収音位置推定部１２により推定された収音位置の精度を評価した。
ここで、基準音の周波数帯域を、２０−４０［ｋＨｚ］、３０−５０［ｋＨｚ］、４０−６０［ｋＨｚ］の３通りとし、基準音信号としてスイープ正弦波のパルスを用いた。基準音信号の信号長ｔ_Ｅを、６４［サンプル］、１２８［サンプル］、２５６［サンプル］、５１２［サンプル］、１０２４［サンプル］、２０４８［サンプル］の６通りとした。また、最大到達時間ｔ_ｒを、４０９６［サンプル］とした。

候補位置の集合であるグリッドマップとして、縦、横の大きさが、それぞれ１［ｍ］の２次元のマップを用いた。但し、グリッドマップの中心を原点Ｏにおいた。
収音位置の算出に用いる到達時間τ_{ｍｎ［ｊ］}を定める際、式（１）〜（３）の３通りの手法を用い、それぞれについて３０回試行した。但し、到達時間τ_{ｍｎ［ｊ］}を定めるための基準音の再生周期ｊは、各１周期のみとした。そして、推定された収音位置ξ_Ｍｍのｘ座標ξ_Ｍｍ，ｘ、ｙ座標ξ_Ｍｍ，ｙの少なくともいずれか一方が−０．３［ｍ］未満もしくは０．３［ｍ］より大きい試行を無効な試行として判定し、ξ_Ｍｍ，ｘ、ξ_Ｍｍ，ｙがそれぞれ−０．３［ｍ］以上０．３［ｍ］以下となる試行を有効な試行として判定する。そして、無効な試行の回数の全試行回数に対する割合を、失敗率（ＦＲ：Ｆａｉｌｕｒｅ Ｒａｔｅ）として定めた。

図７、図８、図９は、それぞれ式（１）、（２）、（３）に示す手法を用いて算出した到達時間τ_{ｍｎ［ｊ］}に基づいて推定された収音位置ξ_Ｍｍの推定結果を示す。但し、いずれも特性長ｄ_Mを０．１２［ｍ］とした。
各行は基準音の信号長ｔ_Ｅを示し、各列は基準音の周波数帯域を示す。各条件での推定結果は、誤差の平均値、標準偏差、失敗率を含む。誤差として、推定された収音位置ξ_Ｍｍと現実の収音素子３０−ｍの位置との間のユークリッド距離を用いた。平均値、標準偏差、失敗率は、それぞれ±記号の左隣、±記号の右隣、カッコ内に示されている。平均値、標準偏差、失敗率は、それぞれ各条件での複数の試行ならびに収音素子間で算出される。

図７−図９に示す推定結果を周波数帯域間で比較すると、３０−５０［ｋＨｚ］の帯域に係る推定結果が最も良好となる傾向が認められる。図９の第６行に示す例では、２０−４０［ｋＨｚ］、３０−５０［ｋＨｚ］、４０−６０［ｋＨｚ］の帯域のそれぞれについて、平均値は１７．０、１１．９、１３．３［ｍｍ］となった。
図７−図９に示す推定結果を信号長ｔ_Ｅ間で比較すると、信号長ｔ_Ｅ間が長いほど推定結果が良好になる傾向が認められる。図８の第３列に示す例では、６４、１２８、２５６、５１２、１０２４、２０４８［サンプル］の信号長ｔ_Ｅのそれぞれについて、平均値は、１５．０、１５．５、１４．６、１４．９、１４．９、１３．７［ｍｍ］となった。
図７に示す推定結果では、第６行第３列に示すように、信号長ｔ_Ｅが１０２４［サンプル］、周波数帯域が３０−５０［ｋＨｚ］となる基準音に係る推定結果が最良と判断される。図８に示す推定結果では、第７行第３列に示すように、信号長ｔ_Ｅが２０４８［サンプル］、周波数帯域が３０−５０［ｋＨｚ］となる基準音に係る推定結果が最良と判断される。図９に示す推定結果では、第７行第３列に示すように、信号長ｔ_Ｅが２０４８［サンプル］、周波数帯域が３０−５０［ｋＨｚ］となる基準音に係る推定結果が最良と判断される。

図７−図９間で推定結果を比較すると、図８、図９に示す推定結果は、図７に示す推定結果よりも良好となる傾向がある。このことは、式（１）に示す手法よりも、式（２）、（３）に示す手法を用いた方が収音位置ξ_Ｍｍの推定精度が向上することを示す。
なお、図７−図９に示す推定結果のうち、失敗率が有意に０％よりも高くなるケースは、信号長ｔ_Ｅが６４［サンプル］、周波数帯域が４０−６０［ｋＨｚ］である場合だけであって、その他のケースでは、いずれも失敗率は０％となった。

図１０は、収音位置の推定結果の到達時間τ_{ｍｎ［ｊ］}の算出手法と特性長ｄ_Mの依存性の例を示す。
図１０は、各行に到達時間τ_{ｍｎ［ｊ］}の算出手法を示し、各列に特性長ｄ_Mを示す。到達時間τ_{ｍｎ［ｊ］}の算出手法間で推定結果を比較すると、式（３）に示す手法に基づく推定結果が、他の手法に基づく推定結果よりも良好である。例えば、特性長ｄ_Mが０．１２［m］である場合、式（１）、（２）、（３）のそれぞれについて平均値は、１５．９［ｍｍ］、１３．７［ｍｍ］、１１．４［ｍｍ］となった。また、特性長ｄ_Mが大きくなるほど、誤差も増加する傾向がある。例えば、式（３）に手法に基づく手法では、特性長ｄ_Mが０．０４［ｍ］、０．０８［ｍ］、０．１２［m］のそれぞれについて平均値は、７．２［ｍｍ］、１０．４［ｍｍ］、１１．４［ｍｍ］となった。

なお、図６に示す円Ｃ１、Ｃ２、Ｃ３は、特性長ｄ_Mが０．１２［ｍ］であって式（３）を用いて算出された到達時間τ_１３、τ_２３、τ_３３にそれぞれ対応する距離ｃ・τ_１３、ｃ・τ_２３、ｃ・τ_３３をそれぞれ半径とし、収音素子３０−１、３０−２、３０−３の収音位置を中心とする円である。２次元平面においては、少なくとも２個の到達時間を知得できれば、その到達時間に基づく円の交点に相当する位置が収音位置として推定されるはずである。しかしながら、算出される到達時間には誤差が含まれるため、２つの円、例えば、円Ｃ１、Ｃ２の交点の位置は、現実の収音位置、例えば、収音素子３０−３の位置から離れてしまう。また、各２つの円の交点は、必ずしも収音素子３０−３の近傍に存在するとは限らない。図６に示す例でも、むしろ他の収音素子に近接する交点や、収音素子が配置されている領域外に所在する交点も存在する。上述のように、収音位置算出部１２４は、式（４）に示す関係を用いて収音素子３０−ｍが存在する可能性の度合いを示す尤度を算出することにより、算出した尤度に基づいて確からしい収音位置ξ_Ｍｍを定めることができる。

図１１は、半径ｄ_Ｅが０．０４［ｍ］である場合に、図６に示す基準音源２０−１〜２０−３ならびに収音素子３０−１〜３０−３のもとで得られた到達時間τ_１３、τ_２３、τ_３３を用いて算出される尤度の分布を示す。濃淡は尤度を示し、明るい部分ほど尤度が高く、暗い部分ほど尤度が低いことを示す。この例では、収音素子３０−３の近傍において尤度が最も高くなる。このことは、式（４）に示す関係を用いて収音位置ξ_Ｍｍを正確に定めることができることを裏付ける。

第２の実験では、音源方向推定部１４により式（１５）に示す関係を用いて推定された音源方向の精度を評価した。但し、図１２に示すように基準音源２０−１〜２０−３からなる音源部２０の中心と、収音素子３０−１〜３０−５からなる収音部３０の中心を原点に配置した。また、６名の話者Ｓｐ０１−Ｓｐ０６を音源として参加させ、それぞれ原点Ｏの方向に向けて配置させた。話者Ｓｐ０１、Ｓｐ０３、Ｓｐ０４、Ｓｐ０６の位置を、矩形の領域の各頂点に配置し、話者Ｓｐ０２、Ｓｐ０５の位置を、原点Ｏからそれぞれξ_ｙ方向、ξ_ｙ方向の逆方向に０．７［ｍ］離れた位置に配置した。矩形の領域は、ξ_ｘ方向の座標ξ_ｘが−０．７［ｍ］以上０．７［ｍ］以下であって、ξ_ｙ方向の座標ξ_ｙが−０．７［ｍ］以上０．８［ｍ］以下となる領域である。また、話者Ｓｐ０１−Ｓｐ０６の高さを、いずれも原点Ｏから０．４［ｍ］高い位置とした。

まず、音源数Ｌ、すなわち６名の話者のうち現実に発話する話者の数として、２、３、４、５、６の５通りについて、音源推定部１３が推定した音源数を評価した。発話音声として、ＡＴＲ（Ａｄｖａｎｃｅｄ Ｔｅｌｅｃｏｍｍｕｎｉｃａｔｉｏｎ Ｒｅｓｅａｒｃｈ；国際電気通信基礎技術研究所）デジタル音声データベースの２１６単語セットのうち最初の１０単語の音声を用いた。当該データベースは、音声学的にバランスが取れた日本語の複数の単語からなる。各話者位置について、６個の単語セットを用いた。６個の単語セットは全６０単語からなり、６個のうち３個の単語セットは、女性の発話であり、残りの３個の単語セットは男性の発話である。そして、各話者の１０単語の発話を音源数の推定ならびに式（１１）に示すステアリングベクトルの生成に用いた。

図１３は、音源推定部１３により推定された音源数を推定結果として示す。図１３の各行は、話者数、つまり現実の音源数を示し、各列は、特性長ｄ_Ｍを示す。各行各列の推定結果は、ＳＳＡＰパターンに基づく共通と判定される音源の棄却前の音源数と棄却後の音源数をその順序で示す。各条件ともに音源数が概ね話者数と一致し、音源数が適切に推定されることが示される。特性長ｄ_Ｍが０．０４［ｍ］であって、話者数が２、３である場合には、推定された音源数である２．４、３．１と差が生じるが、その差は０．４以下となる。また、ＳＳＡＰパターンによる音源の棄却により、推定される音源数２、３が発話者数にそれぞれ一致することは、その処理により余分に検出された音源が棄却され、より正確に音源数が推定されることを示す。
また、特性長ｄ_Ｍが０．１２［ｍ］である場合には、音源数が正確に推定される。このことは、柔構造ロボットに設置され、配置が可変なマイクロホンアレイでも、その大きさを拡張すれば音源数をより正確に推定できることを示す。

図１４は、音源方向推定部１４により式（１５）に示す関係を用いて推定された音源方向を推定結果として示す。図１４は、各行に到達時間τ_{ｍｎ［ｊ］}の算出手法を示し、各列に特性長ｄ_Mを示す。各列各行に推定誤差の平均値、標準偏差を、それぞれ±記号の左隣、右隣に示す。推定誤差は、現実の音源方向と推定された音源方向θ_［ｉ］との差である。特性長ｄ_M間で推定結果を比較すると、特性長ｄ_Mが大きくなるほど、推定誤差が小さくなる傾向が認められる。また、式（３）の手法に基づく到達時間τ_{ｍｎ［ｊ］}に基づく音源方向が最も推定結果が良好となる。例えば、特性長ｄ_Mが０．１２［m］であり、式（３）の手法を用いる場合には、推定誤差の平均値、標準偏差は、それぞれ５．６５°、３．７７°となった。図１４に示す推定結果も、マイクロホンアレイの大きさを拡張すれば音源方向をより正確に推定できることを示す。

以上に説明したように、本実施形態に係る音響処理装置１０は、収音位置推定部１２を備える。収音位置推定部１２は、基準音源２０−１〜２０−Ｎ（Ｎは３以上）のそれぞれから収音素子３０−ｍまでの基準音の到達時間τ_{ｍｎ［ｊ］}である第１時間を収音素子３０−ｍが収音した収音信号を用いて定める。そして、収音位置推定部１２は、基準音源２０−１〜２０−Ｎのそれぞれから収音素子３０−ｍの位置の候補である候補位置ξ_Ｇまでの基準音の到達時間τ_Ｇである第２時間を算出する。その後、収音位置推定部１２は、第１時間と第２時間との差が小さいほど高い確率を与える確率関数を用いて、候補位置が収音素子の位置である可能性を示す尤度を算出し、当該尤度を最大化する候補位置ξ_Ｇを、収音素子３０−ｍの位置ξ_Ｍｍとして定める。
この構成によれば、３個以上の基準音源２０−１〜２０−Ｎから収音素子３０−ｍまでの基準音の到達時間τ_{ｍｎ［ｊ］}と基準音源２０−１〜２０−Ｎから候補位置ξ_Ｇまでの到達時間τ_Ｇにより候補位置ξ_Ｇが収音素子３０−ｍの位置である確率を用いて尤度が算出される。そして、尤度が最大となる候補位置ξ_Ｇが収音素子３０−ｍの位置ξ_Ｍｍとして定まる。そのため、３個以上の基準音源２０−１〜２０−Ｎを用いて候補位置ξ_Ｇが収音素子３０−ｍの位置ξ_Ｍｍとなる可能性が総合的に評価されるので、到達時間τ_{ｍｎ［ｊ］}に誤差が含まれていても、基準音源２０−１〜２０−Ｎが移動しているか否かに関わらず、より確実に収音素子３０−ｍの位置ξ_Ｍｍを定めることができる。

また、収音位置推定部１２は、基準音を示す基準音信号と収音素子３０−ｍが収音した収音信号ｘ_ｍ（ｔ）との相互相関関数に基づいて第１時間を算出する。
この構成によれば、相互相関関数により基準音信号と収音信号ｘ_ｍ（ｔ）との相関性を定量的に評価して、基準音の再生から収音までの第１時間を定めることができる。

また、収音位置推定部１２は、離散時刻ごとに算出した相互相関関数を所定の連続関数で補間して得られる関数値を最大化する時刻に基づいて第１時間を算出する。
この構成によれば、離散時刻ごとに相互相関関数を定め、相互に隣接する２つの離散時刻間の時刻に対する相互相関関数を推定し、推定された相互相関関数が最大となる時刻が基準音の再生から収音までの第１時間として定められる。そのため、時間離散化による第１時間の推定誤差、ひいては収音位置ξ_Ｍｍの推定誤差を低減することができる。

また、音響処理装置１０は、相互相関関数の周波数領域における変換係数の位相成分を算出する。そして、位相成分を時間領域に逆変換して得られる離散時刻ごとの逆変換関数値を所定の連続関数で補間して得られる関数値を最大化する時刻までの期間を第１時間として算出する。
この構成によれば、相互相関関数の周波数領域における変換係数の振幅成分が除去されるので、相互相関関数の振幅の周波数依存性が解消される。そして、相互相関関数の振幅成分を除去して得られる時間領域の逆変換関数値が相互に隣接する２つの離散時刻間で補間される。そのため、時間離散化に加え、基準音信号の振幅の周波数特性による第１時間、ひいては収音位置ξ_Ｍｍの推定誤差を低減することができる。

また、基準音は可聴帯域よりも高い成分の成分からなる超音波である。
そのため、収音位置ξ_Ｍｍの推定に可聴帯域の音声よりも波長が短い超音波が用いられるので、第１時間の推定誤差、ひいては収音位置ξ_Ｍｍの推定誤差を低減することができる。また、基準音は人間に知覚されないため、人間に対する受聴環境が害されない。

また、収音位置推定部１２は、音源推定部１３と、音源方向推定部１４とを備える。
音源推定部１３は、所定の収音素子が収音した基準チャネルの収音信号と他の収音素子が収音した他チャネルの収音信号との到達時間差を要素とする到達時間差ベクトルを所定期間ごとに生成する。そして、音源推定部１３は、到達時間差ベクトルをクラスタごとに分類し、基準音源とは別個の音源に対応するクラスタを代表する到達時間差ベクトルに基づいて第１ステアリングベクトルを定める。
音源方向推定部１４は、第１ステアリングベクトルに基づく相関行列を固有値展開して固有ベクトルを算出し、音源方向に対応する応答係数を要素として有する第２ステアリングベクトルを生成し、第２ステアリングベクトルの各要素の二乗和を、２次以上の各次の前記固有ベクトルと前記第２ステアリングベクトルとの内積の総和で正規化した空間スペクトルを最大化する音源方向を定める。
この構成によれば、Ｍチャネルの収音信号ｘ_ｍ（ｔ）から個々の音源に対応する第１ステアリングベクトルが定まり、第１ステアリングベクトルに基づく相関行列に対するＭ個の固有ベクトルが算出される、そのうち、第１固有ベクトルは、基準音源とは別個の目的音源に対するステアリングベクトルが張られる信号部分空間の基底ベクトルとなり、第２固有ベクトルから第Ｍ固有ベクトルのいずれとも直交する。音源方向θに対応する第２ステアリングベクトルと第ｍ固有ベクトルとの内積の値が最も小さくなるので、音源方向に対する値が最大となるように空間スペクトルが算出される。そのため、音源数が未知であっても、位置が可変な収音素子３０−ｍで収音される収音信号ｘ_ｍを用いて、個々の音源の音源方向を正確に推定することができる。

以上、図面を参照してこの発明の実施形態について詳しく説明してきたが、具体的な構成は上述のものに限られることはなく、この発明の要旨を逸脱しない範囲内において様々な設計変更等をすることが可能である。

例えば、上述の音響処理システム１において、音響処理装置１０と音源部２０が別体である場合を例にしたが、これには限られない。音響処理装置１０は、音源部２０と一体に構成されてもよい。
また、収音素子３０−１〜３０−Ｍは、必ずしも柔構造ロボットに設置されていなくてもよい。収音素子３０−１〜３０−Ｍは、１個もしくは複数の他の物体、例えば、車両模型、内視鏡装置、ドローンなどの移動体、通信機器、計測器、その他、収音を主用途としない機器に設置されてもよい。
入出力部１１と、基準音源２０−１〜２０−Ｎならびに収音素子３０−１〜３０−Ｍとの間の各種の信号の入出力は、有線でもよいし、無線でもよい。

上述の実施形態に係る音響処理装置１０は、２次元空間における収音素子３０−１〜３０−Ｍもしくはその他の音源の方向もしくは位置の推定を実行する場合を例にしたが、これには限られない。音響処理装置１０は、３次元空間における収音素子３０−１〜３０−Ｍもしくはその他の音源の方向もしくは位置の推定を実行する場合に応用されてもよい。
音響処理装置１０は、少なくとも収音位置推定部１２を備えていればよい。音響処理装置１０は、収音位置推定部１２が推定した収音位置を用いて、収音素子３０−１〜３０−Ｍで取得される収音信号に対するその他のアレイ処理を実行するアレイ処理部を備えてもよい。その場合、音源推定部１３と音源方向推定部１４が省略されてもよい。アレイ処理部は、例えば、音源方向推定部１４とは別個の方式（例えば、ＭＵＳＩＣ法）を用いた音源方向推定、雑音抑圧、音源分離などのいずれかの処理、もしくはそれらの処理の任意の組み合わせを実行してもよい。

なお、上述の実施形態における音響処理装置１０の一部、例えば、収音位置推定部１２、音源推定部１３および音源方向推定部１４の全部または一部をコンピュータで実現するようにしてもよい。その場合、この制御機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現してもよい。なお、ここでいう「コンピュータシステム」とは、音響処理装置１０に内蔵されたコンピュータシステムであって、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）などの１個以上のプロセッサの他、ＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）や周辺機器等のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＣＤ−ＲＯＭ（Ｃｏｍｐａｃｔ Ｃｏｍｐａｃｔ Ｄｉｓｋ）等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含んでもよい。また上記プログラムは、前述した機能の一部を実現するためのものであってもよく、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであってもよい。

また、上述した実施形態及び変形例における音響処理装置１０の一部、例えば、または全部を、ＬＳＩ（Ｌａｒｇｅ Ｓｃａｌｅ Ｉｎｔｅｇｒａｔｉｏｎ）等の集積回路として実現してもよい。音響処理装置１０の各機能ブロックは個別にプロセッサ化してもよいし、一部、または全部を集積してプロセッサ化してもよい。また、集積回路化の手法はＬＳＩに限らず専用回路、または汎用プロセッサで実現してもよい。また、半導体技術の進歩によりＬＳＩに代替する集積回路化の技術が出現した場合、当該技術による集積回路を用いてもよい。

１…音響処理システム、１０…音響処理装置、１１…入出力部、１２…収音位置推定部、１３…音源推定部、１４…音源方向推定部、２０…音源部、２０−１〜２０−Ｎ…基準音源、３０…収音部、３０−１〜３０−Ｍ…収音素子、１２１…基準音源制御部、１２２…到達時間算出部、１２３…候補位置情報記憶部、１２４…収音位置算出部、１３１…到達時間差ベクトル生成部、１３２…音源特定部

Claims (8)

1. Ｎ（Ｎは、３以上の整数）個の基準音源のそれぞれから収音素子までの基準音の到達時間である第１時間を前記収音素子が収音した収音信号を用いて定め、
   前記Ｎ個の基準音源のそれぞれから前記収音素子の位置の候補である候補位置までの前記基準音の到達時間である第２時間を取得し、
   前記第１時間と前記第２時間との差が小さいほど高い確率を与える確率関数を用いて、前記候補位置が収音素子の位置である可能性を示す尤度を算出し、前記尤度を最大化する候補位置を、前記収音素子の位置として定める収音位置推定部
   を備える音響処理装置。
2. 前記収音位置推定部は、
   前記基準音を示す基準音信号と前記収音素子が収音した収音信号との相互相関関数に基づいて前記第１時間を算出する
   請求項１に記載の音響処理装置。
3. 前記収音位置推定部は、
   離散時刻ごとに算出した前記相互相関関数を所定の連続関数で補間して得られる関数値を最大化する時刻に基づいて前記第１時間を算出する
   請求項２に記載の音響処理装置。
4. 前記収音位置推定部は、
   前記相互相関関数の周波数領域における変換係数の位相成分を算出し、前記位相成分を時間領域に逆変換して得られる離散時刻ごとの逆変換関数値を所定の連続関数で補間して得られる関数値を最大化する時刻までの期間を、
   前記第１時間として算出する
   請求項２に記載の音響処理装置。
5. 前記基準音は、可聴帯域よりも高い周波数の成分からなる
   請求項１から請求項４のいずれか一項に記載の音響処理装置。
6. 所定の収音素子が収音した基準チャネルの収音信号と他の収音素子が収音した他チャネルの収音信号との到達時間差を要素とする到達時間差ベクトルを所定期間ごとに生成し、
   前記到達時間差ベクトルをクラスタごとに分類し、前記基準音源とは別個の音源に対応するクラスタを代表する到達時間差ベクトルに基づいて第１ステアリングベクトルを定める音源推定部と、
   前記第１ステアリングベクトルに基づく相関行列を固有値展開して固有ベクトルを算出し、
   音源方向に対応する応答係数を要素として有する第２ステアリングベクトルを生成し、
   前記第２ステアリングベクトルの各要素の二乗和を、２次以上の各次の前記固有ベクトルと前記第２ステアリングベクトルとの内積の総和で正規化した空間スペクトルを最大化する音源方向を定める音源方向推定部と、
   をさらに備える請求項１から請求項５のいずれか一項に記載の音響処理装置。
7. 音響処理装置における音響処理方法であって、
   Ｎ（Ｎは、３以上の整数）個の基準音源のそれぞれから収音素子までの基準音の到達時間である第１時間を前記収音素子が収音した収音信号を用いて定め、
   前記Ｎ個の基準音源のそれぞれから前記収音素子の位置の候補である候補位置までの前記基準音の到達時間である第２時間を取得し、
   前記第１時間と前記第２時間との差が小さいほど高い確率を与える確率関数を用いて、前記候補位置が収音素子の位置である可能性を示す尤度を算出し、前記尤度を最大化する候補位置を、前記収音素子の位置として定める収音位置推定ステップ
   を有する音響処理方法。
8. 音響処理装置のコンピュータに、
   Ｎ（Ｎは、３以上の整数）個の基準音源のそれぞれから収音素子までの基準音の到達時間である第１時間を前記収音素子が収音した収音信号を用いて定め、
   前記Ｎ個の基準音源のそれぞれから前記収音素子の位置の候補である候補位置までの前記基準音の到達時間である第２時間を取得し、
   前記第１時間と前記第２時間との差が小さいほど高い確率を与える確率関数を用いて、前記候補位置が収音素子の位置である可能性を示す尤度を算出し、前記尤度を最大化する候補位置を、前記収音素子の位置として定める収音位置推定手順
   を実行させるためのプログラム。