JP6531050B2

JP6531050B2 - 音源定位装置、方法、及びプログラム

Info

Publication number: JP6531050B2
Application number: JP2016032365A
Authority: JP
Inventors: 弘和亀岡; 惇鈴木
Original assignee: Nippon Telegraph and Telephone Corp
Current assignee: Nippon Telegraph and Telephone Corp
Priority date: 2016-02-23
Filing date: 2016-02-23
Publication date: 2019-06-12
Anticipated expiration: 2036-02-23
Also published as: JP2017150903A

Description

本発明は、音源定位装置、方法、及びプログラムに係り、特に、音響信号から、音源の位置を推定する音源定位装置、方法、及びプログラムに関する。

波源定位は、レーダやソナーといった幅広い応用を有している。特に、小さいアレイで、移動する波源を瞬時に定位し追跡できるようにすることは重要課題である。波源定位問題に対する従来法としては、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) アルゴリズムにより複数の音源の波源定位を行う。

＜本発明の実施の形態の原理＞
次に、音源の位置を推定する原理について説明する。

＜音源拘束偏微分方程式＞
図１に示すように、観測点の基準となる位置ベクトルを

とし、単一波源の位置ベクトルを

とする。波源の信号をg(t)、音速をc とし、単一点波源からの球面波伝播を仮定すると観測点における観測値は

と表される。ここで、

である。観測点から波源方向へ向かう単位ベクトルをn とすると、

であるため、f(r, t) の空間微分は、

となる。また、f(r, t) の時間微分は

となるので、式(1) と式(8) を式(7) に代入することでgが消去され、

のように、観測信号とその時間・空間微分のみを含む方程式を立てることができる。ただし、Ｒ = ｜r − r₀｜は観測点から波源までの距離である。この式を音源拘束式と呼ぶ（上記非特許文献１〜３）。以上のように音源拘束式は、任意の音源信号波形で成り立つ、音源の位置と空間の場の一意な関係を記述する偏微分方程式である。

＜音源拘束偏微分方程式に基づく音響信号の確率モデル化＞
図２のようなマイクロホンアレイで、観測信号の空間微分を空間差分で近似する場合を考える。図２に示す観測信号ｆの空間微分を取得するためのアレイ幾何の例では、例えばｘ方向のｆの空間微分は、（f_1,t −f_2,t）／２Ｄで近似できる。

ただし、マイクロホンアレイの配置は、観測信号の空間微分を空間差分で近似できるものであれば良く、以下の理論は図２の配置に限らない。図２のマイクロホンアレイの場合、７本のマイクロホンを用いて各時刻t_lで、基準点における信号f_0,l およびその各方向の空間差分

を得ることができる。ただし、lは離散時刻のインデックスを表す。

基準点における観測信号の時間微分を時間差分で近似することにすると、式(9) は

と表せる。ただし、n_x、n_y、n_z はそれぞれのx, y, z 方向の成分、T はサンプリング周期である。

式(10) の左辺を右辺に移項し、整理すると

が得られる。式(11) の右辺は雑音の存在や差分近似に伴う誤差により実際には必ずしも厳密に0 にはならない。そこで、式(11) の右辺を

のように誤差変数ε_x,l, ε_y,l, ε_z,l に置き換え、これらを平均が0 で互いに独立な正規確率変数（複素正規分布に従う確率変数）

と仮定する。また、観測点における観測信号を、平均が0、分散がσ² ₀の正規確率変数とする。これは、

と仮定することに相当する。

ここで、f_x,l, f_y,l, f_z,l, f_0,l を並べたベクトルとε_x,l, ε_y,l, ε_z,l, ε_0,l を並べたベクトルを

とし、f_0,0, f₁,..., f_L を連結したベクトルとε_0,0, ε₁,...,ε_L を連結したベクトルを

と表記すると、式(12) は

の形で書ける。ただし、θ= ｛R,n｝であり、A(θ) は

で与えられる。式(13), (15) より、εは平均が0、分散共分散行列が

の複素正規分布

に従う。

（A(θ) は正則）であるので、f は

と表され、式(27) より、

が言える。従って、観測信号およびその空間差分が与えられた下での最尤音源位置^θは

により得られる。

＜複数音源の定位アルゴリズム＞
以上のf の確率モデル化により、音源が複数個存在する場合、および雑音が存在する場合の観測信号の確率分布を導くことができる。音源インデックスをk とし、音源k に由来する観測信号の成分、音源位置パラメータをそれぞれf^(k)、θ^(k) とする。また、f^(k) の分散をσ^(k) ₀ ² とする。式(32) より、

となる。また、加法雑音をｖとし、観測信号を

とする。f⁽¹⁾,・・・,f^(K),ｖが互いに独立であれば、観測信号y は

に従う。ただし、Γはｖの分散共分散行列である。以上より、複数の音源と雑音が存在する場合の各音源の最尤音源位置

は、観測信号ｙが与えられた下で

を解くことにより得られる。

y を不完全データ、

を完全データと見なすことで、以上の最尤推定問題に対しExpectation-Maximization (EM) アルゴリズムを適用することができる。完全データ対数尤度log p(x｜θ) は

で与えられるので、y が与えられた下でのlog p(x｜θ) のx に関する条件付き期待値（Q 関数）は

で与えられる。ただし、

はx に関係する項のみについての等号を意味する。この関数が増大するようにθを更新するステップ（M ステップ）と、更新したθをθ´に代入し、

と

を計算するステップ（E ステップ）を繰り返すことでp(y｜θ) を局所的に最大にするθを求めることができる。

完全データx と不完全データy の関係は

と書けるので、

はそれぞれ

で与えられる。以上より、以下の初期設定、Ｅステップ、Ｍステップからなるアルゴリズムを得る。

（初期ステップ）
θを初期設定する。

（E ステップ）
θをθ´に代入し、式(41) により

を計算する。

（Ｍステップ）
下式によりθを更新する。

＜Mステップ更新式＞
Ｍステップでは、

ができるだけ大きくなるように

を更新する。

を最大にする

の同時最適解を解析的に求めることは難しいが、座標勾配法によりそれぞれの変数に関して

が最大となるように反復更新することで

を局所最大化することができる（EM アルゴリズムでは、M ステップで補助関数が単調に増大することが保証されていれば収束性は保証される）。以下に、M ステップの更新方法を２例示す。

＜例１＞
＜n^(k)の更新式＞
n^(k)は単位ベクトルなので、

の下で

ができるだけ小さくなるようにn^(k) を更新する。この制約つき最適化問題は、例えば、

のようなラグランジアンを用いてLagrange 未定乗数法で解くことができる。Ａ(θ^(k)) は

のようにn^(k)に依存する項とそうでない項に分解できるので、Ｌ(n^(k)) のn^(k) に関する偏微分を0 と置くことにより、

を得る。ただし、E_i,jは、i 行j 列目の要素のみが1 で残りは0 であるような4×4 行列である。

より、あとは

となるようにγ^(k)を二分法などで探索し、式(52)に代入すれば良い。

＜音源距離R^(k)の更新式＞

とする。上記と同様、Ａ(θ^(k)) は

のようにρ^(k) に依存する項とそうでない項に分解することができるので、Q(θ,θ´) のρ^(k)に関する偏微分を0 と置くことにより、

を得る。

＜σ^(k) ₀ ² の更新式＞
上記と同様、Σ^(k)-1を

のようにσ^(k) _0,m ² に依存する項とそうでない項に分解できるので、Q(θ,θ´) のσ^(k) ₀ ² に関する偏微分を0 と置くことにより、

を得る。

＜雑音分散共分散行列Γの更新式＞
雑音の分散共分散行列を

のように、正規化分散共分散行列モデルW と雑音のエネルギーν²の積で表し、ν² _mを変数とする。W は空間無相関モデルや拡散音場モデルなどから導かれる定数行列である。Q(θ,θ´) のν² に関する偏微分を0 と置くことにより、

を得る。

＜例２＞
＜音源位置ベクトルr(k)の更新式＞
この例では例１と到来方向の更新の仕方のみが異なる。この例では

を変数とする。この場合はノルムの制約は不要なので、制約なし最適化問題として、Q(θ，θ´) を最大にするr(k) を求めれば良い。A(θ(k))は

のようにr(k) に依存する項とそうでない項に分解できるので，Q(θ，θ´)のr(k) に関する偏微分を0と置くことにより、

を得る。音源距離R^(k)の更新式、σ^(k) ₀ ² の更新式、雑音分散共分散行列Γの更新式は例１と同様である。

＜逆行列計算＞
式(41) より、Ｅステップでは

の逆行列計算が必要である。ここでは点音源が一つと雑音源が一つの場合と、点音源が二つの場合にこの逆行列計算が効率的に行えることを示す。

＜１音源と１雑音源の場合＞
１音源と１雑音源の場合、(ＨΛＨ^T)^-1 は

と書ける。ここで、Woodbury の公式

を用いると、式（６１）は

と書ける。Γはブロック対角行列、Ｖ₁は帯行列（ブロック三重対角行列）なので、Γ^-1＋Ｖ₁ は帯行列となり、(Γ^-1＋Ｖ₁)Ｖ₁ の計算はCholesky 分解を用いて効率的に行うことができる。

＜２音源の場合＞
２音源の場合，(ＨΛＨ^T)^-1 は

と書ける。上記同様、Woodbury の公式を用いると、式(65) は

と書ける。Ｖ₁ とＶ₂ はいずれも帯行列なので、Ｖ₂＋Ｖ₁も帯行列となり、(Ｖ₂＋Ｖ₁)‐1Ｖ₁の計算はCholesky 分解を用いて効率的に行うことができる。

＜システム構成＞
次に、マイクロホンアレイにより入力された音響信号から、複数の音源の位置を推定する音源定位装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。

図３に示すように、本発明の第１の実施の形態に係る音源定位装置１００は、ＣＰＵと、ＲＡＭと、音源定位処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。

図３に示すように、音源定位装置１００は、入力部１０と、演算部２０と、出力部９０とを備えている。

入力部１０は、上記図２に示すようなマイクロホンアレイの各マイクロホンから出力された、複数の音源からの音源信号が混じっている音響信号（以後、観測信号）の時系列データを受け付ける。

演算部２０は、空間差分算出部２２と、音源位置推定部２５と、を含んで構成されている。

空間差分算出部２２は、マイクロホンアレイの各マイクロホンから出力された観測信号から、各時刻t_lで、基準点のマイクロホンにおける観測信号f_0,lを取得すると共に、以下の式に従って、各方向ｘ、ｙ、ｚの空間差分ｆ_x,l，ｆ_y,l，ｆ_z,lを算出する。

音源位置推定部２５は、空間差分算出部２２により得られた、基準点のマイクロホンにおける各時刻ｔ_lの観測信号f_0,lと、各時刻ｔ_lの各方向ｘ、ｙ、ｚの空間差分ｆ_x,l，ｆ_y,l，ｆ_z,lとからなる観測信号ｙに基づいて、ＥＭアルゴリズムを用いて、音源拘束偏微分方程式を用いて定められた、複数の音源と加法雑音が存在する場合における、複数の音源の各々の位置を条件とした、観測信号ｙの確率分布Ｐ（ｙ｜θ）を大きくするように、複数の音源の各々の位置を推定する。

音源位置推定部２５は、期待値算出部２６と、変数更新部２８と、収束判定部３０とを備えている。

期待値算出部２６は、空間差分算出部２２により得られた、基準点のマイクロホンにおける各時刻ｔ_lの観測信号f_0,lと、各時刻ｔ_lの各方向ｘ、ｙ、ｚの空間差分ｆ_x,l，ｆ_y,l，ｆ_z,lとからなる観測信号ｙと、初期設定された、又は前回更新された各音源ｋの音源位置θ^(k)と、初期設定された、又は前回更新された各音源ｋの音源信号の分散σ^(k) ₀ ²に基づいて、上記式(41) により

を計算する。

変数更新部２８は、期待値算出部２６によって算出された

に基づいて、Q(θ,θ´) ができるだけ大きくなるように、上記式（５２）、式（５４）、式（５６）〜式（５８）に従って、各音源ｋの方向ベクトルn^(k),音源距離R^(k),分散σ^(k) ₀ ²,雑音共分散行列Γ を更新する。なお、上述した例２のように、上記式（６０）、式（５４）、式（５６）〜式（５８）に従って、各音源ｋの位置ベクトルｒ^(k),音源距離R^(k),分散σ^(k) ₀ ²,雑音共分散行列Γ を更新する。

収束判定部３０は、予め定められた収束判定条件を満たすまで、期待値算出部２６及び変数更新部２８による各処理を繰り返させる。収束判定条件としては、例えば、予め定められた繰り返し回数に到達することである。

収束判定条件を満たしたときに、最終的に得られた各音源ｋの方向ベクトルn^(k)、音源距離R^(k)を、各音源ｋの位置の推定結果として、出力部９０により出力する。

＜音源定位装置の作用＞
次に、本実施の形態に係る音源定位装置１００の作用について説明する。

入力部１０において、マイクロホンアレイの各マイクロホンから出力された観測信号の時系列データを受け付けると、音源定位装置１００は、図４に示す音源定位処理ルーチンを実行する。

まず、ステップＳ１２０では、マイクロホンアレイの各マイクロホンから入力された観測信号の時系列データから、各時刻t_lで、基準点のマイクロホンにおける観測信号f_0,lを取得すると共に、各方向ｘ、ｙ、ｚの空間差分ｆ_x,l，ｆ_y,l，ｆ_z,lを算出する。

ステップＳ１２２では、各音源ｋの音源位置θ^(k)と各音源ｋの音源信号の分散σ^(k) ₀ ²とに初期値を設定する。

そして、ステップＳ１２３では、上記ステップＳ１２０で取得した各時刻ｔ_lの観測信号f_0,l、及び各時刻ｔ_lの各方向ｘ、ｙ、ｚの空間差分ｆ_x,l，ｆ_y,l，ｆ_z,lからなる観測信号ｙと、上記ステップＳ１２２で初期設定された、又は後述するステップＳ１２４で前回更新された各音源ｋの音源位置θ^(k)及び分散σ^(k) ₀ ²とに基づいて、上記式(41) により

を計算する。

ステップＳ１２４では、上記ステップＳ１２３で算出された

に基づいて、Q(θ,θ´) ができるだけ大きくなるように、上記式（５２）、式（５４）、式（５６）〜式（５８）に従って、に従って、各音源ｋの方向ベクトルn^(k),音源距離R^(k),分散σ^(k) ₀ ²,雑音共分散行列Γ を更新する。

ステップＳ１２５において、予め定められた収束判定条件を満たしたか否かを判定し、収束判定条件を満たしていない場合には、上記ステップＳ１２３へ戻る。一方、収束判定条件を満たした場合には、ステップＳ１２６へ進む。

ステップＳ１２６では、上記ステップＳ１２４で最終的に得られた各音源ｋの方向ベクトル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 度の違いは正解と見なしてよいため、どの初期値からも真の音源方向を正しく推定できていることが分かる。

以上説明したように、本実施の形態に係る音源定位装置によれば、音源拘束偏微分方程式を用いて定められた、複数の音源と加法雑音が存在する場合における、複数の音源の各々の位置を条件とした、基準のマイクロホンの各時刻の観測信号、及び複数の方向の各々に対する各時刻の観測信号の確率分布を大きくするように、複数の音源の各々の位置を推定することにより、雑音が存在する場合であっても、複数の音源を同時に定位することができる。

なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

例えば、上述の音源定位装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。

また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。

１０入力部
２０演算部
２２空間差分算出部
２５音源位置推定部
２６期待値算出部
２８変数更新部
３０収束判定部
９０出力部
１００音源定位装置

Claims

マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置であって、
複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出する空間差分算出部と、
前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号と、前記複数の方向の各々に対して算出された前記観測信号の差分とに基づいて、音源拘束偏微分方程式を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各時刻の観測信号、及び前記複数の方向の各々に対する各時刻の観測信号の差分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定する音源位置推定部と、
を含み、
前記確率密度値は、以下の式で表わされる音源定位装置。

ただし、ｆ(k)が、音源ｋからの観測される音源信号を表し、Γは、前記加法雑音の分散共分散行列であり、θ ^(k) は、音源ｋの位置を表し、σ _x ^(k)2 、σ _y ^(k)2 、σ _z ^(k)2 は、音源ｋからの観測される音源信号の方向ｘ、ｙ、ｚの差分における分散を表し、σ ₀ ^(k)2 は、音源ｋからの観測される音源信号の前記基準のマイクロホンにおける分散を表し、Ｒは、音源までの距離を表し、ｃは、音速を表し、Ｔは、サンプリング周期を表し、ｎ _x 、ｎ _y 、ｎ _z は、音源へ向かう単位ベクトルの方向ｘ、ｙ、ｚの成分を表す。
前記音源位置推定部は、ＥＭ（Expectation-Maximization）アルゴリズムにより、前記確率密度値が大きくなるように、複数の音源ｋの各々までの距離Ｒ^(k)、前記複数の音源ｋの各々へ向かう単位ベクトルｎ^(k)、前記複数の音源ｋの各々の前記基準のマイクロホンにおける分散₀ ^(k)2、及び前記加法雑音の分散共分散行列Γを繰り返し更新することにより、前記複数の音源の各々の位置を推定する請求項１記載の音源定位装置。
マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置における音源定位方法であって、
空間差分算出部が、複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出し、
音源位置推定部が、前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号と、前記複数の方向の各々に対して算出された前記観測信号の差分とに基づいて、音源拘束偏微分方程式を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各時刻の観測信号、及び前記複数の方向の各々に対する各時刻の観測信号の差分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定し、
前記確率密度値は、以下の式で表わされる音源定位方法。

ただし、ｆ(k)が、音源ｋからの観測される音源信号を表し、Γは、前記加法雑音の分散共分散行列であり、θ ^(k) は、音源ｋの位置を表し、σ _x ^(k)2 、σ _y ^(k)2 、σ _z ^(k)2 は、音源ｋからの観測される音源信号の方向ｘ、ｙ、ｚの差分における分散を表し、σ ₀ ^(k)2 は、音源ｋからの観測される音源信号の前記基準のマイクロホンにおける分散を表し、Ｒは、音源までの距離を表し、ｃは、音速を表し、Ｔは、サンプリング周期を表し、ｎ _x 、ｎ _y 、ｎ _z は、音源へ向かう単位ベクトルの方向ｘ、ｙ、ｚの成分を表す。
前記音源位置推定部が推定することでは、ＥＭ（Expectation-Maximization）アルゴリズムにより、前記確率密度値が大きくなるように、複数の音源ｋの各々までの距離Ｒ^(k)、前記複数の音源ｋの各々へ向かう単位ベクトルｎ^(k)、前記複数の音源ｋの各々の前記基準のマイクロホンにおける分散₀ ^(k)2、及び前記加法雑音の分散共分散行列Γを繰り返し更新することにより、前記複数の音源の各々の位置を推定する請求項３記載の音源定位方法。
請求項１又は２に記載の音源定位装置の各部としてコンピュータを機能させるためのプログラム。