"音データを検出する音検出部と、検出された音データに基づいて、方向及び時間ごとのMUSICスペクトルを演算するMUSICスペクトル演算部と、各方向の音源の有無による状態、状態の遷移を記述する状態遷移モデル、及び音源有りの状態及び音源無しの状態に観測されるMUSICスペクトルを記述する観測モデルを定め、MUSICスペクトルの時系列データに基づいて、前記観測モデル及び前記状態遷移モデルのモデルパラメータの事後分布を推定する、モデルパラメータ推定部と、推定されたモデルパラメータの事後分布に基づいて、方向及び時間ごとの音源存在事後確率のパーティクルをサンプルすることによって、音源定位を行う音源定位部と、を備えた音源定位システム。前記観測モデルとしてガウス混合モデルを使用する請求項1に記載の音源定位システム。音データを検出するステップと、検出された音データに基づいて、方向及び時間ごとのMUSICスペクトルを演算するステップと、各方向の音源の有無による状態、状態の遷移を記述する状態遷移モデル、及び音源有りの状態及び音源無しの状態に観測されるMUSICスペクトルを記述する観測モデルを定め、MUSICスペクトルの時系列データに基づいて、前記観測モデル及び前記状態遷移モデルのモデルパラメータの事後分布を推定するステップと、推定されたモデルパラメータの事後分布に基づいて、方向及び時間ごとの音源存在事後確率のパーティクルをサンプルすることによって、音源定位を行うステップと、を含む音源定位方法。前記観測モデルとしてガウス混合モデルを使用する請求項3に記載の音源定位方法。音源定位を行うステップは、P個のパーティクルをサンプルするステップと、各パーティクルについて重みを計算するステップと、各パーティクルの重みを正規化するステップと、各パーティクルの重みを使用してリサンプリングを行うステップと、を含む請求項3または4に記載の音源定位方法。したがって、動的に音環境が変化する状況において、MUSICスペクトルの閾値を自動的に定めることのできる音源定位システム及び音源定位方法に対するニーズがある。本発明の第1の態様による音源定位システムは、音データを検出する音検出部と、検出された音データに基づいて、方向及び時間ごとのMUSICスペクトルを演算するMUSICスペクトル演算部と、を備えている。本システムは、各方向の音源の有無による状態、状態の遷移を記述する状態遷移モデル、及び音源有りの状態及び音源無しの状態に観測されるMUSICスペクトルを記述する観測モデルを定め、MUSICスペクトルの時系列データに基づいて、前記観測モデル及び前記状態遷移モデルのモデルパラメータの事後分布を推定する、モデルパラメータ推定部と、推定されたモデルパラメータの事後分布に基づいて、方向及び時間ごとの音源存在事後確率のパーティクルをサンプルすることによって、音源定位を行う音源定位部と、をさらに備えている。本態様の音源定位システムによれば、観測モデル及び状態遷移モデルのモデルパラメータの事後分布を推定し、推定されたモデルパラメータの事後分布に基づいて音源定位を行うので、動的に音環境が変化する状況においても、手動で閾値を設定する必要なく、ロバストに音源定位を行うことができる。本発明の第1の態様の第1の実施形態による音源定位システムは、前記観測モデルとしてガウス混合モデルを使用する。本実施形態によれば、ガウス分布を使用することで解析的な計算が可能となる。本発明の第2の態様による音源定位システムは、音データを検出するステップと、検出された音データに基づいて、方向及び時間ごとのMUSICスペクトルを演算するステップと、を含む。本方法は、各方向の音源の有無による状態、状態の遷移を記述する状態遷移モデル、及び音源有りの状態及び音源無しの状態に観測されるMUSICスペクトルを記述する観測モデルを定め、MUSICスペクトルの時系列データに基づいて、前記観測モデル及び前記状態遷移モデルのモデルパラメータの事後分布を推定するステップと、推定されたモデルパラメータの事後分布に基づいて、方向及び時間ごとの音源存在事後確率のパーティクルをサンプルすることによって、音源定位を行うステップと、をさらに含む。本態様の音源定位方法によれば、観測モデル及び状態遷移モデルのモデルパラメータの事後分布を推定し、推定されたモデルパラメータの事後分布に基づいて音源定位を行うので、動的に音環境が変化する状況においても、手動で閾値を設定する必要なく、ロバストに音源定位を行うことができる。本発明の第2の態様の第1の実施形態による音源定位方法は、前記観測モデルとしてガウス混合モデルを使用する。本実施形態によれば、ガウス分布を使用することで解析的な計算が可能となる。本発明の第2の態様の第2の実施形態による音源定位方法において、音源定位を行うステップは、P個のパーティクルをサンプルするステップと、各パーティクルについて重みを計算するステップと、各パーティクルの重みを正規化するステップと、各パーティクルの重みを使用してリサンプリングを行うステップと、を含む。本実施形態によれば、推定されたモデルパラメータの分布に基づいてパーティクルをサンプルすることによって、方向及び時間ごとの音源事後確率のパーティクルを簡単な手順で求めることができる。図1は、本発明の一実施形態による音源定位システム100の構成を示す図である。音源定位システム100は、音検出部101と、MUSICスペクトル演算部103と、モデルパラメータ推定部105と、音源定位部107とを含む。音検出部101は、M個のマイクロフォンからなるマイクロフォンアレイであってもよい。図2は、M個のマイクロフォン1011からなるマイクロフォンアレイ101の構成を示す図である。図2において、M=8である。一例として、8個のマイクロフォン1011は水平面上に配置されており、音源定位システム100は、音源が、水平面上のどの方向に存在するかを定める。一例として、方向の解像度は5°であり、360°/5°=72のいずれかの方向を定める。たとえば、音検出部であるマイクロフォンアレイによって、Mチャネルの音響信号が得られる。また、周波数ビンごとに水平面上のD方向D=72からの伝達関数が与えられているとする。音源定位システム100は、N個の音源方向を定める。ただし、同時に定位可能な最大音源数Nmaxは、マイクロフォンの数未満である。N≦Nmax<MMUSICスペクトル演算部103におけるMUSICスペクトルの算出方法について説明する。より詳細は、文献ROSchmidtMultipleEmitterLocationandSignalParameterEstimationIEEETransonAntennasandPropagationvol34no3pp2762801986PDanesandJBonnalInformationTheoreticDetectionofBroadbandSourcesinaCoherentBeamspaceMUSICSchemeinProcofIROS20102011pp19761981に記載されている。MUSIC法は、時間周波数領域において適用される。具体的には、サンプリング周波数16000Hzで、窓長12pt、シフト幅160ptの短時間フーリエ変換を行っている。をMチャネル音響信号の時間フレームτ、周波数ビンωにおける入力Mチャネル音響信号の複素振幅ベクトルとする。各周波数ビンω、ΔT[sec]間隔の時刻tに対して、入力信号の自己相関行列の計算、の固有値分解、固有ベクトルと伝達関数を用いたMUSICスペクトルの計算を行う。上記の乃至について以下に説明する。1入力信号の自己相関行列の計算入力信号の自己相関行列は、時間ΔTで観測したサンプル値の相関として計算する。ただし、はエルミート転置、は時刻tにおける時間フレームを表す。入力ベクトルのM個の要素は各チャネルに対応する。2固有値分解を次のように固有値分解する。ここで、は固有ベクトル、は固有値からなる対角行列である。は、とのM個の固有ベクトルで表せ、となる。固有値は降順に並べられているものとする。入力信号にN個の音源が含まれる場合、固有値からまでが、音源のエネルギに対応する大きな値を有する。それに対し、残りの固有値からまではマイクロフォンの観測ノイズなどに対応する小さな値を有する。ここで、重要な点は、からのノイズに対応する固有ベクトルは、音源方向に対応する伝達関数ベクトルと直交するという点であるROSchmidtMultipleEmitterLocationandSignalParameterEstimationIEEETransonAntennasandPropagationvol34no3pp2762801986。固有ベクトルと伝達関数を用いたMUSICスペクトルの計算MUSICスペクトルは以下の式にしたがって計算する。ただし、は方向d、周波数ビンωに対応するM次元の伝達関数ベクトルである。これらの伝達関数は、マイクロフォンアレイを使用して事前に測定したものである。ここで、観測されうる最大の音源数はNmaxである。そのため、からまでの固有ベクトルは、音源方向dに対応する伝達関数と直交する。したがって、式の分母は音源方向のdに対して0となる。すなわち、式のMUSICスペクトルは発散する。ただし、実際には、壁からの反射音などの影響でMUSICスペクトルは発散せず鋭いピークとして観測されることが多い。つぎに、周波数ビンごとのMUSICスペクトルを以下の式にしたがって合算する。ここで、は、周波数ビンωにおける最大固有値である。本実施形態では、音声信号を対象とするため、とした。つぎに、モデルパラメータ推定部105の機能について説明する。モデルパラメータ推定部105は、変分ベイズ隠れマルコフモデルを使用する。状態ベクトルとしてD次元の2値ベクトルを使用する。各次元のベクトル値は、その方向に音源が存在するか否かを示す。また、MUSICスペクトルがガウス分布にしたがう観測値であるとみなし、音源有りの場合及び音源無しの場合に対応するガウス分布からなるガウス混合分布を観測モデルとする。ガウス分布を使用する理由は、複数の周波数ビンの値を加算して対数とったMUSICスペクトルが近似的にガウス分布とみなせることと、ガウス分布を使用することで解析的な計算が可能となることである。図3は、対数スケールのMUSICスペクトルの分布を示す図である。図3の横軸は、対数スケールのMUSICスペクトルを示す。対数スケールのMUSICスペクトルは、以下の式によって求める。図3の縦軸は、観測数を表す。図3に点線で示す音源が存在しないオフのときのガウス分布は、狭くかつ値の小さなMUSICスペクトルの領域に形成され、図3に実線で示す音源が存在するオンときのガウス分布は、広くかつ値の大きなMUSICスペクトルの領域を覆っている。モデルパラメータ推定部105で使用される観測モデルは以下の式で表せる。ただし、はx=yのときそうでなければ、を表す。また、は、平均μ、精度λの正規分布の確率密度関数を表す。パラメータμとλには、正規ガンマ分布を使用する。ただし,Nは平均m精度Lの正規分布であり、以下の式で表せる。Gは形状a尺度bのガンマ分布であり、以下の式で表せる。正規分布のβ、ガンマ分布のaは事前分布の影響の大きさを表す。本実施形態では、事前情報よりも学習時に観測されるデータを重視するため、に設定する。m0は平均パラメータμの事前情報で与えられた平均値を示す。本実施形態では、およそ25程度が適している。あるいは、VB-HMMの学習に用いる観測値のサンプル平均でも問題ない。b0は精度パラメータλの事前情報で与えられたばらけ具合を示す。実験的に500に設定している。あるいは、VB-HMMの学習に用いる観測値のサンプル分散程度でも問題ない。図4は、VB-HMMの確率変数間の条件付き独立性を示すグラフィカルモデルを示す図である。VB-HMMは、状態遷移確率のパラメータθkや、観測確率のパラメータμ、λが数値ではなく、確率変数として扱われる点が通常のHMMと異なる。モデルパラメータ推定部105は、これらのパラメータの確率分布を学習する。モデルパラメータ推定部105で使用される状態遷移モデルは基本的に、各方向ビンdについて、前状態で音源がない場合と音源がある場合から、次状態で音源が出現する、継続する、消滅するといった遷移を考える。本実施形態においては、移動する音源についても考慮するために、以下の表1に示すように前状態の組み合わせから成る四つの場合を考える。すなわち、前時刻の同方向ビンに音源が存在するかどうかと前時刻の隣接方向ビンのいずれかに音源が存在するかによって分類する。たとえば、θ1は前時刻に当該方向d及び隣接ビンd±1に音源が存在しない状態から音源が出現する確率、θ2は、前時刻に方向dに音源が存在しないが、隣接ビンd±1には音源が存在したため、その音源が方向dに移動してきてとなる確率を表す。状態遷移確率は、以下の式によって表せる。ここで、は表1に従って、方向ビンdの周りの前状態の値によって条件kに合致するときにその他の場合は0を返す条件識別関数である。初期状態としては、音源は存在しない、すなわちすべてのdに対してとする。状態遷移パラメータであるには、式の共役事前分布としてベータ分布を用いる。ただしはパラメータc、dを持つβ分布の確率密度関数である。モデルパラメータ推定部105におけるVB-HMMの学習は、事後分布を以下のように因数分解可能な分布に近似して推定する。1Tは、時刻1からTまでの確率変数の集合を表す。一般的なVB-HMMの推論については、文献MJBealVariationalAlgorithmsforApproximateBayesianInferencePhDdissertationGatsbyComputationalNeuroscienceUnitUniversityColledgeLondon2003に記載されている。はそれぞれのkに対し、式に示すパラメータを持つベータ分布となり、は、式で表されるパラメータを持つ正規ガウス分布となる。ただし、変数stdjはstd0のとき、std01、また、std1のときstd11となる変数である。式、に用いられる正規分布の十分統計量はと定義する。また、は式の分布による期待値演算子である。各時刻の状態変数と状態遷移の期待値は次のように計算する。ただし、αとβはそれぞれ前向き、後ろ向き再帰式により計算される。ここで、である。遷移、観測確率の幾何平均は以下の式で表せる。ここで、は、以下の式で定義されるディガンマ関数である。式及びは、ともに添え字j、kを動かしたとき総和が1になるように正規化されている。は、状態遷移の条件kに関する前向き確率である。図5は、モデルパラメータ推定部105による、モデルパラメータの分布を推定する手順を示す流れ図である。図5のステップS1010において、モデルパラメータ推定部105は、初期値を定める。初期値は、式及び式の値に対して、たとえば、以下手順によって定めることができる。式の左辺は、時刻t、方向ビンdにおいて音源が存在しない場合はstd01かつstd10、音源が存在する場合はstd00かつstd11となる2値変数の期待値である。観測値xtdが適当な閾値を超える場合は、たとえば、0810802のように設定する。08の代わりに1にしてもほぼ同様に動作する。式の左辺も、xtdが閾値を超えるか否かに従って算出する。この値は、std101の2通りと、k1~4のいずれかでfk1となる4通りの組み合わせで、計8通りの組み合わせが存在する。xtdの閾値処理と、前時刻のxt1dの閾値処理、前時刻近傍方向のxt1d±1の閾値処理結果から、表1を参照して、該当するkの値を求める。たとえば、前時刻xt1dは閾値以下だが、隣接ビンのxt1d1では閾値を超えた場合はk2にあたる。xtdが閾値を超えた場合、08、その他の7通りの組み合わせでは、7のように設定する。図5のステップS1020において、モデルパラメータ推定部105は、式及び式を使用して、遷移、観測確率の幾何平均を求める。図5のステップS1030において、モデルパラメータ推定部105は、ステップS1020で求めた遷移、観測確率の幾何平均、式及び式を使用して、αとβを計算する。図5のステップS1040において、モデルパラメータ推定部105は、ステップS1030で求めたαとβ、式及び式を使用して、各時刻の状態変数と状態遷移の期待値を求める。図5のステップS1050において、モデルパラメータ推定部105は、ステップS1040で求めた各時刻の状態変数と状態遷移の期待値、式から式を使用して、モデルパラメータの事後分布を推定する。図5のステップS1060において、モデルパラメータ推定部105は、収束判定を行う。具体的に、式で算出されるパラメータβmabの値の変動がなくなることによって収束判定を行う。収束していないと判断されれば、ステップS1020に戻る。収束していると判断されれば、処理を終了する。つぎに、音源定位部107の機能について説明する。音源定位部107は、モデルパラメータ推定部105が推定した、モデルパラメータの事後分布に基づいて、パーティクルフィルタを使用して、複数音源の存在事後確率計算を行う。パーティクルフィルタの推定対象はMUSICスペクトルの時系列データが与えられたときの、各方向ビンにおける音源存在事後確率である。この分布を、P個のパーティクルを用いて以下のように近似計算する。ただしはパーティクルpの重み、は状態ベクトルの値である。図6は、音源定位部107が、各方向ビンにおける音源存在事後確率を表す、P個のパーティクルを求める手順を示す流れ図である。図6のステップS2010において、音源定位部107は、P個のパーティクルをサンプルする。Pは以下のように定める。Pが大きいほど式の近似は正確になるが、Pの大きさに比例した計算時間がかかる。そこで、Pを決める一般的な手続きとしては、実用的な近似を得るほど十分大きなPを与え、このようなPの処理時間が大きすぎる場合は、Pの大きさを適宜縮小する。本実施形態においては、近似結果が収束し、理速度も十分高速な値として、P=500とした。P個のパーティクルのサンプルは、以下の式で表せる分布を使用して行う。ただし、xtdが極大値をとるdのとき、C1でその他の場合は、C0となる。上記の分布の重みには、以下の式で表せるマハラノビス距離を使用する。時刻tにおいて,式で計算される分布qは、計D個のビンに対してオンであるか、またはオフであるかの確率を与える。サンプル方法としては、各dについてaC0の場合、j0とする。つまりbC1の場合、j01それぞれの場合の分布qの確率を参照する。たとえば、の場合、0~1の区間から一様乱数を生成し、その値が08以下の場合は、、08を超えた場合はのように値を決定する。図6のステップS2020において、音源定位部107は、以下の式にしたがって、各パーティクルについて重みwpを算出する。式及び式の状態遷移、観測確率は、モデルパラメータ推定部105で使用された式及び式の事後分布で積分消去することで計算することができる。分布の共役性を用いるとこの積分計算は、次のように解析的に求まる。ただし、は、平均m精度l自由度nのStudentt分布である。さらに、最大の音源数をNmaxに抑えるため、状態ベクトルに存在する音源数がNmaxを超える場合には観測確率はとする。図6のステップS2030において、音源定位部107は、各パーティクルの重みwpがとなるよう正規化する。図6のステップS2040において、処理を終了するかどうか判断する。たとえば、スイッチの状態によって処理を終了するかどうか判断してもよい。処理を終了しないと判断した場合には、ステップS2050に進む。そうでなければ処理を終了する。図6のステップS2050において、音源定位部107は、リサンプリングを行う。リサンプリングは、パーティクルpが持つ値を、そのパーティクルの重みに比例した確率で複製することによって行う。方法の一例としては、以下の処理をp’1~Pまで繰り返す。a0~1区間から一様乱数を生成rp’とする。bp1~Pまでi.rp’←rp’wpii.rp’0になったら、として,pのループを抜ける。iii.wp’←1Pcaに戻る。つぎに、評価実験について説明する。評価実験は、本実施形態による音源定位システム及び固定閾値を使用する従来の音源定位システムの性能を比較した。モデルパラメータ推定部105による、オフラインでのVB-HMMでの学習は、一人の話者がマイクロフォンの周囲を発話しながら動く音響信号で行った。図7は、オンラインの音源定位実験に使用した音源の配置を示す図である。二人の話者301及び303が、マイクロフォンアレイ101の周囲を移動しながら発話し、固定されたスピーカ201から音楽が再生されている。オフライン及びオンラインで用いられた信号の長さはともに20secである。パラメータの設定は以下のとおりである。Nmax3α011β01a01b0500また、パーティクル数はP=500とした。実験で使用した室内の残響時間はRT20840msecであった。図8は、従来の音源定位システムによるオンライン音源定位の結果を示す図である。図の横軸は秒単位の時間を示し、図の縦軸は度単位の方向を示す。従来の音源定位システムの閾値はPthres232527に設定されている。図8、、は、それぞれ、閾値が232527の場合の音源定位の結果を示す。図8、、においては、閾値を超えるビンを、音源が存在するとして黒く表示している。図8、、において、固定されたスピーカ及び移動する話者が黒で表示されている。しかし、図8、の実線の囲みで示すように、閾値を低く設定した場合は音源の誤検出が頻発する。図9は、本実施形態の音源定位システムによるオンライン音源定位の結果を示す図である。図の横軸は秒単位の時間を示し、図の縦軸は度単位の方向を示す。本実施形態の音源定位システムの初期値はm0232527に設定されている。図9、、は、それぞれ、初期値が232527の場合の音源定位の結果を示す。図9、、においては、事後分布の音源存在確率が095以上のビンを、音源が存在するとして黒く表示している。図9、、において、固定されたスピーカ及び移動する話者が黒で表示されている。図8、の実線の囲みの位置に対応する、図9、の点線の囲みの位置において音源の誤検出は少ない。このことは、本実施形態の音源定位システムにおいて、学習の初期値にかかわらず音源の誤検出が少ないことを意味する。また、本実施形態の音源定位システムにおいて音源存在確率の閾値を095から100まで動かして結果を検証し、これらの閾値に対してもロバストに同様の結果を示すことを確認した。この結果から、本実施形態の音源定位システムのモデルパラメータ推定部105によるオフライン学習及び音源定位部107によるオンライン定位の枠組みが、自動的に音源定位に適したパラメータに収束することが確認できた。さらに、今回の実験条件から、本実施形態の音源定位方法は、1音源のみで学習した場合であっても、複数音源に対して安定した音源定位結果をもたらすことが実証された。"