(57)【特許請求の範囲】
【請求項１】
マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置であって、
複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出する空間差分算出部と、
前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号と、前記複数の方向の各々に対して算出された前記観測信号の差分とに基づいて、音源拘束偏微分方程式を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各時刻の観測信号、及び前記複数の方向の各々に対する各時刻の観測信号の差分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定する音源位置推定部と、
を含み、
前記確率密度値は、以下の式で表わされる音源定位装置。
【数１】
0006531050.tif 000082
0006531050.tif 000083
0006531050.tif 000084
ただし、ｆ(k)が、音源ｋからの観測される音源信号を表し、Γは、前記加法雑音の分散共分散行列であり、θ(k)は、音源ｋの位置を表し、σx(k)2、σy(k)2、σz(k)2は、音源ｋからの観測される音源信号の方向ｘ、ｙ、ｚの差分における分散を表し、σ0(k)2は、音源ｋからの観測される音源信号の前記基準のマイクロホンにおける分散を表し、Ｒは、音源までの距離を表し、ｃは、音速を表し、Ｔは、サンプリング周期を表し、ｎx、ｎy、ｎzは、音源へ向かう単位ベクトルの方向ｘ、ｙ、ｚの成分を表す。
【請求項２】
前記音源位置推定部は、ＥＭ（Expectation-Maximization）アルゴリズムにより、前記確率密度値が大きくなるように、複数の音源ｋの各々までの距離Ｒ(k)、前記複数の音源ｋの各々へ向かう単位ベクトルｎ(k)、前記複数の音源ｋの各々の前記基準のマイクロホンにおける分散0(k)2、及び前記加法雑音の分散共分散行列Γを繰り返し更新することにより、前記複数の音源の各々の位置を推定する請求項１記載の音源定位装置。
【請求項３】
マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置における音源定位方法であって、
空間差分算出部が、複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出し、
音源位置推定部が、前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号と、前記複数の方向の各々に対して算出された前記観測信号の差分とに基づいて、音源拘束偏微分方程式を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各時刻の観測信号、及び前記複数の方向の各々に対する各時刻の観測信号の差分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定し、
前記確率密度値は、以下の式で表わされる音源定位方法。
【数２】
0006531050.tif 000085
0006531050.tif 000086
0006531050.tif 000087
ただし、ｆ(k)が、音源ｋからの観測される音源信号を表し、Γは、前記加法雑音の分散共分散行列であり、θ(k)は、音源ｋの位置を表し、σx(k)2、σy(k)2、σz(k)2は、音源ｋからの観測される音源信号の方向ｘ、ｙ、ｚの差分における分散を表し、σ0(k)2は、音源ｋからの観測される音源信号の前記基準のマイクロホンにおける分散を表し、Ｒは、音源までの距離を表し、ｃは、音速を表し、Ｔは、サンプリング周期を表し、ｎx、ｎy、ｎzは、音源へ向かう単位ベクトルの方向ｘ、ｙ、ｚの成分を表す。
【請求項４】
前記音源位置推定部が推定することでは、ＥＭ（Expectation-Maximization）アルゴリズムにより、前記確率密度値が大きくなるように、複数の音源ｋの各々までの距離Ｒ(k)、前記複数の音源ｋの各々へ向かう単位ベクトルｎ(k)、前記複数の音源ｋの各々の前記基準のマイクロホンにおける分散0(k)2、及び前記加法雑音の分散共分散行列Γを繰り返し更新することにより、前記複数の音源の各々の位置を推定する請求項３記載の音源定位方法。
【請求項５】
請求項１又は２に記載の音源定位装置の各部としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音源定位装置、方法、及びプログラムに係り、特に、音響信号から、音源の位置を推定する音源定位装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
波源定位は、レーダやソナーといった幅広い応用を有している。特に、小さいアレイで、移動する波源を瞬時に定位し追跡できるようにすることは重要課題である。波源定位問題に対する従来法としては、Multiple Signal Classication (MUSIC) 法、Generalized Cross-Correlation methods with Phase Transform (GCC-PHAT) 法、波源拘束偏微分方程式に基づく手法（非特許文献１〜３）などがある。
【０００３】
MUSIC 法やGCC-PHAT 法は、音源に対し平面波を仮定し各音源のセンサ間での到来時間差を定位の手がかりとするため、一般にアレイサイズは大きい方が有利となる。また、いずれもセンサアレイの受信信号間の自己相関関数や相互相関関数といった、統計量に基づく手法であるため、音源を高い精度で定位するためには観測時間幅を十分長く取る必要がある。このため、これらの手法は小さいアレイサイズと瞬時的な観測のみによる波源定位には必ずしも向いていない。一方、波源拘束偏微分方程式に基づく手法は、各時刻ごとに成立する音響信号の時空間偏微分方程式を元に音源定位を行うもので、理論的には瞬時の小領域観測のみで波源定位を行うことが可能である。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】藤田悠哉, 小野順貴, 安藤繁, "有限時間窓と離散フーリエ変換の利用を可能にする音源定位の高速厳密解法とその実験" 日本音響学会2006 年秋季研究発表会講演論文集, 3-1-3, pp. 483-484, Sep. 2006.
【非特許文献２】S. Ando, N. Ono, T. Nara, "Direct algebraic method for sound source localization with nest resolution both in time and frequency," in Proc. ICSV14, Jul. 2007.
【非特許文献３】小山翔一, 栗原徹, 安藤繁, "偏微分方程式の空間荷重積分による瞬時音源定位," 日本音響学会2008 年秋季研究発表会講演論文集, 2-8-20, pp. 679-682, Sep. 2008.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、上記の波源拘束偏微分方程式に基づく手法は単一波源に対して成立する方程式をベースとしているため、複数の音源を同時に定位することはできない。また、雑音が存在する場合など、観測音響信号が偏微分方程式から逸脱する場合に脆弱であるという欠点を有している。
【０００６】
本発明は、上記事情を鑑みてなされたものであり、雑音が存在する場合であっても、複数の音源を同時に定位することができる音源定位装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００７】
上記の目的を達成するために本発明に係る音源定位装置は、マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置であって、複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出する空間差分算出部と、前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号と、前記複数の方向の各々に対して算出された前記観測信号の差分とに基づいて、音源拘束偏微分方程式を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各時刻の観測信号、及び前記複数の方向の各々に対する各時刻の観測信号の差分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定する音源位置推定部と、を含んで構成されている。
【０００８】
本発明に係る音源定位方法は、マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置における音源定位方法であって、空間差分算出部が、複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出し、音源位置推定部が、前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号と、前記複数の方向の各々に対して算出された前記観測信号の差分とに基づいて、音源拘束偏微分方程式を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各時刻の観測信号、及び前記複数の方向の各々に対する各時刻の観測信号の差分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定する。
【０００９】
本発明に係るプログラムは、上記の音源定位装置の各部としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００１０】
以上説明したように、本発明の音源定位装置、方法、及びプログラムによれば、音源拘束偏微分方程式を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各時刻の観測信号、及び前記複数の方向の各々に対する各時刻の観測信号の差分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定することにより、雑音が存在する場合であっても、複数の音源を同時に定位することができる、という効果が得られる。
【図面の簡単な説明】
【００１１】
【図１】点音源から観測点ｒへ到来する球面波を示す図である。
【図２】マイクロホンアレイの配置の一例を示す図である。
【図３】本発明の実施の形態に係る音源定位装置の構成を示す概略図である。
【図４】本発明の実施の形態に係る音源定位装置における音源定位処理ルーチンの内容を示すフローチャートである。
【図５】マイクロホンアレイ基準点と音源位置を示す図である。
【図６】一音源、雑音なし、壁の反射係数0.01 の場合の定位結果を示す図である
【図７】一音源、雑音あり、壁の反射係数0.01 の場合の定位結果を示す図である。
【図８】一音源、雑音あり、壁の反射係数0.25 の場合の定位結果を示す図である。
【図９】一音源、雑音あり、壁の反射係数0.5 の場合の定位結果を示す図である。
【図１０】マイクロホンアレイ基準点と音源位置を示す図である。
【図１１】二音源、雑音あり、壁の反射係数0.01 の場合の定位結果を示す図である
【発明を実施するための形態】
【００１２】
以下、図面を参照して本発明の実施の形態を詳細に説明する。本発明で提案する技術は、音響信号から波源位置を推定することを目的とした信号処理技術である。
【００１３】
＜本発明の実施の形態の概要＞
本発明の実施の形態は、上述した従来手法の利点を併せ持つ、小領域・瞬時観測による複数音源の波源定位を可能にする技術である。
【００１４】
本発明の実施の形態では、音源拘束偏微分方程式の時間領域表現をベースにした音響信号の確率分布を構築することで複数の音源と加法雑音が存在する場合の混合信号の確率分布を記述し、Expectation-Maximization (EM) アルゴリズムにより複数の音源の波源定位を行う。
【００１５】
＜本発明の実施の形態の原理＞
次に、音源の位置を推定する原理について説明する。
【００１６】
＜音源拘束偏微分方程式＞
図１に示すように、観測点の基準となる位置ベクトルを
【数１】
0006531050.tif 000002
とし、単一波源の位置ベクトルを
【数２】
0006531050.tif 000003
とする。波源の信号をg(t)、音速をc とし、単一点波源からの球面波伝播を仮定すると観測点における観測値は
【数３】
0006531050.tif 000004
と表される。ここで、
【数４】
0006531050.tif 000005
である。観測点から波源方向へ向かう単位ベクトルをn とすると、
【数５】
0006531050.tif 000006
であるため、f(r, t) の空間微分は、
【数６】
0006531050.tif 000007
となる。また、f(r, t) の時間微分は
【数７】
0006531050.tif 000008
となるので、式(1) と式(8) を式(7) に代入することでgが消去され、
【数８】
0006531050.tif 000009
のように、観測信号とその時間・空間微分のみを含む方程式を立てることができる。ただし、Ｒ = ｜r − r0｜は観測点から波源までの距離である。この式を音源拘束式と呼ぶ（上記非特許文献１〜３）。以上のように音源拘束式は、任意の音源信号波形で成り立つ、音源の位置と空間の場の一意な関係を記述する偏微分方程式である。
【００１７】
＜音源拘束偏微分方程式に基づく音響信号の確率モデル化＞
図２のようなマイクロホンアレイで、観測信号の空間微分を空間差分で近似する場合を考える。図２に示す観測信号ｆの空間微分を取得するためのアレイ幾何の例では、例えばｘ方向のｆの空間微分は、（f1,t −f2,t）／２Ｄで近似できる。
【００１８】
ただし、マイクロホンアレイの配置は、観測信号の空間微分を空間差分で近似できるものであれば良く、以下の理論は図２の配置に限らない。図２のマイクロホンアレイの場合、７本のマイクロホンを用いて各時刻tlで、基準点における信号f0,l およびその各方向の空間差分
【数９】
0006531050.tif 000010
0006531050.tif 000011
0006531050.tif 000012
を得ることができる。ただし、lは離散時刻のインデックスを表す。
【００１９】
基準点における観測信号の時間微分を時間差分で近似することにすると、式(9) は
【数１０】
0006531050.tif 000013
と表せる。ただし、nx、ny、nz はそれぞれ のx, y, z 方向の成分、T はサンプリング周期である。
【００２０】
式(10) の左辺を右辺に移項し、整理すると
【数１１】
0006531050.tif 000014
が得られる。式(11) の右辺は雑音の存在や差分近似に伴う誤差により実際には必ずしも厳密に0 にはならない。そこで、式(11) の右辺を
【数１２】
0006531050.tif 000015
のように誤差変数εx,l, εy,l, εz,l に置き換え、これらを平均が0 で互いに独立な正規確率変数（複素正規分布に従う確率変数）
【数１３】
0006531050.tif 000016
と仮定する。また、観測点における観測信号を、平均が0、分散がσ20の正規確率変数とする。これは、
【数１４】
0006531050.tif 000017
と仮定することに相当する。
【００２１】
ここで、fx,l, fy,l, fz,l, f0,l を並べたベクトルとεx,l, εy,l, εz,l, ε0,l を並べたベクトルを
【数１５】
0006531050.tif 000018
とし、f0,0, f1,..., fL を連結したベクトルとε0,0, ε1,...,εL を連結したベクトルを
【数１６】
0006531050.tif 000019
と表記すると、式(12) は
【数１７】
0006531050.tif 000020
の形で書ける。ただし、θ= ｛R,n｝ であり、A(θ) は
【数１８】
0006531050.tif 000021
で与えられる。式(13), (15) より、εは平均が0、分散共分散行列が
【数１９】
0006531050.tif 000022
の複素正規分布
【数２０】
0006531050.tif 000023
に従う。
【００２２】
【数２１】
0006531050.tif 000024
（A(θ) は正則）であるので、f は
【数２２】
0006531050.tif 000025
と表され、式(27) より、
【数２３】
0006531050.tif 000026
が言える。従って、観測信号およびその空間差分が与えられた下での最尤音源位置^θは
【数２４】
0006531050.tif 000027
0006531050.tif 000028
により得られる。
【００２３】
＜複数音源の定位アルゴリズム＞
以上のf の確率モデル化により、音源が複数個存在する場合、および雑音が存在する場合の観測信号の確率分布を導くことができる。音源インデックスをk とし、音源k に由来する観測信号の成分、音源位置パラメータをそれぞれf(k)、θ(k) とする。また、f(k) の分散をσ(k)02 とする。式(32) より、
【数２５】
0006531050.tif 000029
となる。また、加法雑音をｖとし、観測信号を
【数２６】
0006531050.tif 000030
とする。f(1),・・・,f(K),ｖが互いに独立であれば、観測信号y は
【数２７】
0006531050.tif 000031
に従う。ただし、Γはｖの分散共分散行列である。以上より、複数の音源と雑音が存在する場合の各音源の最尤音源位置
【数２８】
0006531050.tif 000032
は、観測信号ｙが与えられた下で
【数２９】
0006531050.tif 000033
を解くことにより得られる。
【００２４】
y を不完全データ、
【数３０】
0006531050.tif 000034
を完全データと見なすことで、以上の最尤推定問題に対しExpectation-Maximization (EM) アルゴリズムを適用することができる。完全データ対数尤度log p(x｜θ) は
【数３１】
0006531050.tif 000035
で与えられるので、y が与えられた下でのlog p(x｜θ) のx に関する条件付き期待値（Q 関数）は
【数３２】
0006531050.tif 000036
で与えられる。ただし、
【数３３】
0006531050.tif 000037
はx に関係する項のみについての等号を意味する。この関数が増大するようにθを更新するステップ（M ステップ）と、更新したθをθ´に代入し、
【数３４】
0006531050.tif 000038
と
【数３５】
0006531050.tif 000039
を計算するステップ（E ステップ）を繰り返すことでp(y｜θ) を局所的に最大にするθを求めることができる。
【００２５】
完全データx と不完全データy の関係は
【数３６】
0006531050.tif 000040
と書けるので、
【数３７】
0006531050.tif 000041
はそれぞれ
【数３８】
0006531050.tif 000042
で与えられる。以上より、以下の初期設定、Ｅステップ、Ｍステップからなるアルゴリズムを得る。
【００２６】
（初期ステップ）
θを初期設定する。
【００２７】
（E ステップ）
θをθ´に代入し、式(41) により
【数３９】
0006531050.tif 000043
を計算する。
【００２８】
（Ｍステップ）
下式によりθを更新する。
【００２９】
【数４０】
0006531050.tif 000044
【００３０】
＜Mステップ更新式＞
Ｍステップでは、
【数４１】
0006531050.tif 000045
ができるだけ大きくなるように
【数４２】
0006531050.tif 000046
を更新する。
【００３１】
【数４３】
0006531050.tif 000047
を最大にする
【数４４】
0006531050.tif 000048
の同時最適解を解析的に求めることは難しいが、座標勾配法によりそれぞれの変数に関して
【数４５】
0006531050.tif 000049
が最大となるように反復更新することで
【数４６】
0006531050.tif 000050
を局所最大化することができる（EM アルゴリズムでは、M ステップで補助関数が単調に増大することが保証されていれば収束性は保証される）。以下に、M ステップの更新方法を２例示す。
【００３２】
＜例１＞
＜n(k)の更新式＞
n(k)は単位ベクトルなので、
【数４７】
0006531050.tif 000051
の下で
【数４８】
0006531050.tif 000052
ができるだけ小さくなるようにn(k) を更新する。この制約つき最適化問題は、例えば、
【数４９】
0006531050.tif 000053
のようなラグランジアンを用いてLagrange 未定乗数法で解くことができる。Ａ(θ(k)) は
【数５０】
0006531050.tif 000054
【数５１】
0006531050.tif 000055
のようにn(k)に依存する項とそうでない項に分解できるので、Ｌ(n(k)) のn(k) に関する偏微分を0 と置くことにより、
【数５２】
0006531050.tif 000056
を得る。ただし、Ei,jは、i 行j 列目の要素のみが1 で残りは0 であるような4×4 行列である。
【００３３】
【数５３】
0006531050.tif 000057
より、あとは
【数５４】
0006531050.tif 000058
となるようにγ(k)を二分法などで探索し、式(52)に代入すれば良い。
【００３４】
＜音源距離R(k)の更新式＞
【数５５】
0006531050.tif 000059
とする。上記と同様、Ａ(θ(k)) は
【数５６】
0006531050.tif 000060
のようにρ(k) に依存する項とそうでない項に分解することができるので、Q(θ,θ´) のρ(k)に関する偏微分を0 と置くことにより、
【数５７】
0006531050.tif 000061
を得る。
【００３５】
＜σ(k)02 の更新式＞
上記と同様、Σ(k)-1を
【数５８】
0006531050.tif 000062
のようにσ(k)0,m2 に依存する項とそうでない項に分解できるので、Q(θ,θ´) のσ(k)02 に関する偏微分を0 と置くことにより、
【数５９】
0006531050.tif 000063
を得る。
【００３６】
＜雑音分散共分散行列Γの更新式＞
雑音の分散共分散行列を
【数６０】
0006531050.tif 000064
のように、正規化分散共分散行列モデルW と雑音のエネルギーν2の積で表し、ν2mを変数とする。W は空間無相関モデルや拡散音場モデルなどから導かれる定数行列である。Q(θ,θ´) のν2 に関する偏微分を0 と置くことにより、
【数６１】
0006531050.tif 000065
を得る。
【００３７】
＜例２＞
＜音源位置ベクトルr(k)の更新式＞
この例では例１と到来方向の更新の仕方のみが異なる。この例では
【数６２】
0006531050.tif 000066
を変数とする。この場合はノルムの制約は不要なので、制約なし最適化問題として、Q(θ，θ´) を最大にするr(k) を求めれば良い。A(θ(k))は
【数６３】
0006531050.tif 000067
のようにr(k) に依存する項とそうでない項に分解できるので，Q(θ，θ´)のr(k) に関する偏微分を0と置くことにより、
【数６４】
0006531050.tif 000068
を得る。音源距離R(k)の更新式、σ(k)02 の更新式、雑音分散共分散行列Γの更新式は例１と同様である。
【００３８】
＜逆行列計算＞
式(41) より、Ｅステップでは
【数６５】
0006531050.tif 000069
の逆行列計算が必要である。ここでは点音源が一つと雑音源が一つの場合と、点音源が二つの場合にこの逆行列計算が効率的に行えることを示す。
【００３９】
＜１音源と１雑音源の場合＞
１音源と１雑音源の場合、(ＨΛＨT)-1 は
【数６６】
0006531050.tif 000070
と書ける。ここで、Woodbury の公式
【数６７】
0006531050.tif 000071
を用いると、式（６１）は
【数６８】
0006531050.tif 000072
と書ける。Γはブロック対角行列、Ｖ1は帯行列（ブロック三重対角行列）なので、Γ-1＋Ｖ1 は帯行列となり、(Γ-1＋Ｖ1)Ｖ1 の計算はCholesky 分解を用いて効率的に行うことができる。
【００４０】
＜２音源の場合＞
２音源の場合，(ＨΛＨT)-1 は
【数６９】
0006531050.tif 000073
と書ける。上記同様、Woodbury の公式を用いると、式(65) は
【数７０】
0006531050.tif 000074
と書ける。Ｖ1 とＶ2 はいずれも帯行列なので、Ｖ2＋Ｖ1も帯行列となり、(Ｖ2＋Ｖ1)‐1Ｖ1の計算はCholesky 分解を用いて効率的に行うことができる。
【００４１】
＜システム構成＞
次に、マイクロホンアレイにより入力された音響信号から、複数の音源の位置を推定する音源定位装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【００４２】
図３に示すように、本発明の第１の実施の形態に係る音源定位装置１００は、ＣＰＵと、ＲＡＭと、音源定位処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【００４３】
図３に示すように、音源定位装置１００は、入力部１０と、演算部２０と、出力部９０とを備えている。
【００４４】
入力部１０は、上記図２に示すようなマイクロホンアレイの各マイクロホンから出力された、複数の音源からの音源信号が混じっている音響信号（以後、観測信号）の時系列データを受け付ける。
【００４５】
演算部２０は、空間差分算出部２２と、音源位置推定部２５と、を含んで構成されている。
【００４６】
空間差分算出部２２は、マイクロホンアレイの各マイクロホンから出力された観測信号から、各時刻tlで、基準点のマイクロホンにおける観測信号f0,lを取得すると共に、以下の式に従って、各方向ｘ、ｙ、ｚの空間差分ｆx,l，ｆy,l，ｆz,lを算出する。
【００４７】
【数７１】
0006531050.tif 000075
0006531050.tif 000076
0006531050.tif 000077
【００４８】
音源位置推定部２５は、空間差分算出部２２により得られた、基準点のマイクロホンにおける各時刻ｔlの観測信号f0,lと、各時刻ｔlの各方向ｘ、ｙ、ｚの空間差分ｆx,l，ｆy,l，ｆz,lとからなる観測信号ｙに基づいて、ＥＭアルゴリズムを用いて、音源拘束偏微分方程式を用いて定められた、複数の音源と加法雑音が存在する場合における、複数の音源の各々の位置を条件とした、観測信号ｙの確率分布Ｐ（ｙ｜θ）を大きくするように、複数の音源の各々の位置を推定する。
【００４９】
音源位置推定部２５は、期待値算出部２６と、変数更新部２８と、収束判定部３０とを備えている。
【００５０】
期待値算出部２６は、空間差分算出部２２により得られた、基準点のマイクロホンにおける各時刻ｔlの観測信号f0,lと、各時刻ｔlの各方向ｘ、ｙ、ｚの空間差分ｆx,l，ｆy,l，ｆz,lとからなる観測信号ｙと、初期設定された、又は前回更新された各音源ｋの音源位置θ(k)と、初期設定された、又は前回更新された各音源ｋの音源信号の分散σ(k)02に基づいて、上記式(41) により
【数７２】
0006531050.tif 000078
を計算する。
【００５１】
変数更新部２８は、期待値算出部２６によって算出された
【数７３】
0006531050.tif 000079
に基づいて、Q(θ,θ´) ができるだけ大きくなるように、上記式（５２）、式（５４）、式（５６）〜式（５８）に従って、各音源ｋの方向ベクトルn(k),音源距離R(k),分散σ(k)02,雑音共分散行列Γ を更新する。なお、上述した例２のように、上記式（６０）、式（５４）、式（５６）〜式（５８）に従って、各音源ｋの位置ベクトルｒ(k),音源距離R(k),分散σ(k)02,雑音共分散行列Γ を更新する。
【００５２】
収束判定部３０は、予め定められた収束判定条件を満たすまで、期待値算出部２６及び変数更新部２８による各処理を繰り返させる。収束判定条件としては、例えば、予め定められた繰り返し回数に到達することである。
【００５３】
収束判定条件を満たしたときに、最終的に得られた各音源ｋの方向ベクトルn(k)、音源距離R(k)を、各音源ｋの位置の推定結果として、出力部９０により出力する。
【００５４】
＜音源定位装置の作用＞
次に、本実施の形態に係る音源定位装置１００の作用について説明する。
【００５５】
入力部１０において、マイクロホンアレイの各マイクロホンから出力された観測信号の時系列データを受け付けると、音源定位装置１００は、図４に示す音源定位処理ルーチンを実行する。
【００５６】
まず、ステップＳ１２０では、マイクロホンアレイの各マイクロホンから入力された観測信号の時系列データから、各時刻tlで、基準点のマイクロホンにおける観測信号f0,lを取得すると共に、各方向ｘ、ｙ、ｚの空間差分ｆx,l，ｆy,l，ｆz,lを算出する。
【００５７】
ステップＳ１２２では、各音源ｋの音源位置θ(k)と各音源ｋの音源信号の分散σ(k)02とに初期値を設定する。
【００５８】
そして、ステップＳ１２３では、上記ステップＳ１２０で取得した各時刻ｔlの観測信号f0,l、及び各時刻ｔlの各方向ｘ、ｙ、ｚの空間差分ｆx,l，ｆy,l，ｆz,lからなる観測信号ｙと、上記ステップＳ１２２で初期設定された、又は後述するステップＳ１２４で前回更新された各音源ｋの音源位置θ(k)及び分散σ(k)02とに基づいて、上記式(41) により
【数７４】
0006531050.tif 000080
を計算する。
【００５９】
ステップＳ１２４では、上記ステップＳ１２３で算出された
【数７５】
0006531050.tif 000081
に基づいて、Q(θ,θ´) ができるだけ大きくなるように、上記式（５２）、式（５４）、式（５６）〜式（５８）に従って、に従って、各音源ｋの方向ベクトルn(k),音源距離R(k),分散σ(k)02,雑音共分散行列Γ を更新する。
【００６０】
ステップＳ１２５において、予め定められた収束判定条件を満たしたか否かを判定し、収束判定条件を満たしていない場合には、上記ステップＳ１２３へ戻る。一方、収束判定条件を満たした場合には、ステップＳ１２６へ進む。
【００６１】
ステップＳ１２６では、上記ステップＳ１２４で最終的に得られた各音源ｋの方向ベクトルn(k),音源距離R(k)を、各音源ｋの位置の推定結果として、出力部９０により出力して、音源定位処理ルーチンを終了する。
【００６２】
＜実験＞
図５のように単一音源とマイクロホンアレイを配置し，以下の条件で雑音・残響環境下の音源定位実験を行った。
【００６３】
（実験条件）
音源数: 1
部屋の反響を考慮した反射係数: 0.01, 0.5, 0.8
観測時間長: 64 点(4ms)
マイク間隔: 1cm
【００６４】
図６〜９に、各条件の定位結果を示す。
【００６５】
図６は、一音源、雑音なしの場合で、壁の反射係数が0.01 の場合の定位結果を示す。各点は異なる初期値ごとに得られる推定音源方向を表す。180 度の違いは正解と見なしてよいため、どの初期値からも真の音源方向を正しく推定できていることが分かる。
【００６６】
図７は、一音源、雑音あり、壁の反射係数0.01 の場合の定位結果を示す。各点は異なる初期値ごとに得られる推定音源方向を表す。
【００６７】
図８は、一音源、雑音あり、壁の反射係数0.25 の場合の定位結果を示す。各点は異なる初期値ごとに得られる推定音源方向を表す。
【００６８】
図９は、一音源、雑音あり、壁の反射係数0.5 の場合の定位結果を示す。各点は異なる初期値ごとに得られる推定音源方向を表す。
【００６９】
また、図１０のように二音源とマイクロホンアレイを配置し、以下の条件で雑音・残響環境下の音源定位実験を行った。
【００７０】
（実験条件）
音源数: 2
部屋の反響を考慮した反射係数：0.01
観測時間長: 64 点(4ms)
マイク間隔: 1cm
【００７１】
図１１に定位結果を示す。図１１は、二音源、雑音あり、壁の反射係数0.01 の場合の定位結果を示す。各点は異なる初期値ごとに得られる推定音源方向を表す。180 度の違いは正解と見なしてよいため、どの初期値からも真の音源方向を正しく推定できていることが分かる。
【００７２】
以上説明したように、本実施の形態に係る音源定位装置によれば、音源拘束偏微分方程式を用いて定められた、複数の音源と加法雑音が存在する場合における、複数の音源の各々の位置を条件とした、基準のマイクロホンの各時刻の観測信号、及び複数の方向の各々に対する各時刻の観測信号の確率分布を大きくするように、複数の音源の各々の位置を推定することにより、雑音が存在する場合であっても、複数の音源を同時に定位することができる。
【００７３】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【００７４】
例えば、上述の音源定位装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【００７５】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【００７６】
１０   入力部
２０   演算部
２２   空間差分算出部
２５   音源位置推定部
２６   期待値算出部
２８   変数更新部
３０   収束判定部
９０   出力部
１００ 音源定位装置
【図１】
0006531050.tif 000088
【図２】
0006531050.tif 000089
【図３】
0006531050.tif 000090
【図４】
0006531050.tif 000091
【図５】
0006531050.tif 000092
【図６】
0006531050.tif 000093
【図７】
0006531050.tif 000094
【図８】
0006531050.tif 000095
【図９】
0006531050.tif 000096
【図１０】
0006531050.tif 000097
【図１１】
0006531050.tif 000098