(57)【特許請求の範囲】
【請求項１】
マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置であって、
複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出する空間差分算出部と、
前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号を入力として、各周波数の観測時間周波数成分を出力すると共に、前記空間差分算出部によって前記複数の方向の各々に対して算出された前記観測信号の差分を入力として、前記複数の方向の各々に対して、各周波数の観測時間周波数成分を出力する時間周波数展開部と、
前記時間周波数展開部により出力された、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分に基づいて、音源拘束偏微分方程式の周波数領域表現を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定する音源位置推定部と、
を含む音源定位装置。
【請求項２】
前記確率密度値は、以下の式で表わされる請求項１記載の音源定位装置。
【数１】
0006529451.tif 000079
0006529451.tif 000080
0006529451.tif 000081
ただし、ｙｍが、前記基準のマイクロホンの周波数ｍの観測時間周波数成分、及び前記複数の方向の各々に対する周波数ｍの観測時間周波数成分を表し、Γｍは、前記加法雑音の周波数ｍの分散共分散行列であり、θ(k)は、音源ｋの位置を表し、σx,m(k)2、σy,m(k)2、σz,m(k)2は、音源ｋからの音源信号の方向ｘ、ｙ、ｚの差分における周波数ｍの成分エネルギーを表し、σ0,m(k)2は、音源ｋからの音源信号の前記基準のマイクロホンにおける周波数ｍの成分エネルギーを表し、Ｒは、音源までの距離を表し、ｃは、音速を表し、Ｌは、周波数のインデッックスを規定するための定数であり、Ｔは、サンプリング周期を表し、ｎx、ｎy、ｎzは、音源へ向かう単位ベクトルの方向ｘ、ｙ、ｚの成分を表す。
【請求項３】
前記音源位置推定部は、ＥＭ（Expectation-Maximization）アルゴリズムにより、前記確率密度値が大きくなるように、複数の音源ｋの各々までの距離Ｒ(k)、前記複数の音源ｋの各々へ向かう単位ベクトルｎ(k)、前記複数の音源ｋの各々の前記基準のマイクロホンにおける周波数ｍの成分エネルギーσ0,m(k)2、及び前記加法雑音の周波数ｍの分散共分散行列Γｍを繰り返し更新することにより、前記複数の音源の各々の位置を推定する請求項２記載の音源定位装置。
【請求項４】
マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置における音源定位方法であって、
空間差分算出部が、複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出し、
時間周波数展開部が、前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号を入力として、各周波数の観測時間周波数成分を出力すると共に、前記空間差分算出部によって前記複数の方向の各々に対して算出された前記観測信号の差分を入力として、前記複数の方向の各々に対して、各周波数の観測時間周波数成分を出力し、
音源位置推定部が、前記時間周波数展開部により出力された、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分に基づいて、音源拘束偏微分方程式の周波数領域表現を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定する
音源定位方法。
【請求項５】
前記確率密度値は、以下の式で表わされる請求項４記載の音源定位方法。
【数２】
0006529451.tif 000082
0006529451.tif 000083
0006529451.tif 000084
ただし、ｙｍが、前記基準のマイクロホンの周波数ｍの観測時間周波数成分、及び前記複数の方向の各々に対する周波数ｍの観測時間周波数成分を表し、Γｍは、前記加法雑音の周波数ｍの分散共分散行列であり、θ(k)は、音源ｋの位置を表し、σx,m(k)2、σy,m(k)2、σz,m(k)2は、音源ｋからの音源信号の方向ｘ、ｙ、ｚの差分における周波数ｍの成分エネルギーを表し、σ0,m(k)2は、音源ｋからの音源信号の前記基準のマイクロホンにおける周波数ｍの成分エネルギーを表し、Ｒは、音源までの距離を表し、ｃは、音速を表し、Ｌは、周波数のインデッックスを規定するための定数であり、Ｔは、サンプリング周期を表し、ｎx、ｎy、ｎzは、音源へ向かう単位ベクトルの方向ｘ、ｙ、ｚの成分を表す。
【請求項６】
前記音源位置推定部が推定することでは、ＥＭ（Expectation-Maximization）アルゴリズムにより、前記確率密度値が大きくなるように、複数の音源ｋの各々までの距離Ｒ(k)、前記複数の音源ｋの各々へ向かう単位ベクトルｎ(k)、前記複数の音源ｋの各々の前記基準のマイクロホンにおける周波数ｍの成分エネルギーσ0,m(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.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、上記の波源拘束偏微分方程式に基づく手法は単一波源に対して成立する方程式をベースとしているため、複数の音源を同時に定位することはできない。また、雑音が存在する場合など、観測音響信号が偏微分方程式から逸脱する場合に脆弱であるという欠点を有している。
【０００６】
本発明は、上記事情を鑑みてなされたものであり、雑音が存在する場合であっても、複数の音源を同時に定位することができる音源定位装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００７】
上記の目的を達成するために本発明に係る音源定位装置は、マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置であって、複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出する空間差分算出部と、前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号を入力として、各周波数の観測時間周波数成分を出力すると共に、前記空間差分算出部によって前記複数の方向の各々に対して算出された前記観測信号の差分を入力として、前記複数の方向の各々に対して、各周波数の観測時間周波数成分を出力する時間周波数展開部と、前記時間周波数展開部により出力された、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分に基づいて、音源拘束偏微分方程式の周波数領域表現を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定する音源位置推定部と、を含んで構成されている。
【０００８】
本発明に係る音源定位方法は、マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置における音源定位方法であって、空間差分算出部が、複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出し、時間周波数展開部が、前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号を入力として、各周波数の観測時間周波数成分を出力すると共に、前記空間差分算出部によって前記複数の方向の各々に対して算出された前記観測信号の差分を入力として、前記複数の方向の各々に対して、各周波数の観測時間周波数成分を出力し、音源位置推定部が、前記時間周波数展開部により出力された、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分に基づいて、音源拘束偏微分方程式の周波数領域表現を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定する。
【０００９】
本発明に係るプログラムは、上記の音源定位装置の各部としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００１０】
以上説明したように、本発明の音源定位装置、方法、及びプログラムによれば、音源拘束偏微分方程式の周波数領域表現を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分の確率密度値を大きくするように、前記複数の音源の各々の位置を推定することにより、雑音が存在する場合であっても、複数の音源を同時に定位することができる、という効果が得られる。
【図面の簡単な説明】
【００１１】
【図１】点音源から観測点ｒへ到来する球面波を示す図である。
【図２】マイクロホンアレイの配置の一例を示す図である。
【図３】本発明の実施の形態に係る音源定位装置の構成を示す概略図である。
【図４】本発明の実施の形態に係る音源定位装置における音源定位処理ルーチンの内容を示すフローチャートである。
【図５】フレーム幅がL = 16 の場合の実験結果を示す図である。
【図６】フレーム幅がL = 32 の場合の実験結果を示す図である。
【図７】フレーム幅がL = 64 の場合の実験結果を示す図である。
【発明を実施するための形態】
【００１２】
以下、図面を参照して本発明の実施の形態を詳細に説明する。本発明で提案する技術は、音響信号から波源位置を推定することを目的とした信号処理技術である。
【００１３】
＜本発明の実施の形態の概要＞
本発明の実施の形態は、上述した従来手法の利点を併せ持つ、小領域・瞬時観測による複数音源の波源定位を可能にする技術である。
【００１４】
本発明の実施の形態では、音源拘束偏微分方程式の周波数領域表現をベースにした音響信号およびその空間差分の確率分布を構築することで複数の音源と加法雑音が存在する場合の混合信号およびその空間差分の確率分布を記述し、Expectation-Maximization (EM) アルゴリズムにより複数の音源の波源定位を行う。
【００１５】
＜本発明の実施の形態の原理＞
次に、音源の位置を推定する原理について説明する。
【００１６】
＜音源拘束偏微分方程式＞
図１に示すように、観測点の基準となる位置ベクトルを
【数１】
0006529451.tif 000002
とし、単一波源の位置ベクトルを
【数２】
0006529451.tif 000003
とする。波源の信号をg(t)、音速をc とし、単一点波源からの球面波伝播を仮定すると観測点における観測値は
【数３】
0006529451.tif 000004
と表される。ここで、
【数４】
0006529451.tif 000005
である。観測点から波源方向へ向かう単位ベクトルをn とすると、
【数５】
0006529451.tif 000006
であるため、f(r, t) の空間微分は、
【数６】
0006529451.tif 000007
となる。また、f(r, t) の時間微分は
【数７】
0006529451.tif 000008
となるので、式(1) と式(8) を式(7) に代入することでgが消去され、
【数８】
0006529451.tif 000009
のように、観測信号とその時間・空間微分のみを含む方程式を立てることができる。ただし、Ｒ = ｜r − r0｜は観測点から波源までの距離である。この式を音源拘束式と呼ぶ（上記非特許文献１〜３）。以上のように音源拘束式は、任意の音源信号波形で成り立つ、音源の位置と空間の場の一意な関係を記述する偏微分方程式である。
【００１７】
＜音源拘束偏微分方程式に基づく音響信号の確率モデル化＞
図２のようなマイクロホンアレイで、観測信号の空間微分を空間差分で近似する場合を考える。図２に示す観測信号ｆの空間微分を取得するためのアレイ幾何の例では、例えばｘ方向のｆの空間微分は、（f1,t −f2,t）／２Ｄで近似できる。
【００１８】
ただし、マイクロホンアレイの配置は、観測信号の空間微分を空間差分で近似できるものであれば良く、以下の理論は図２の配置に限らない。図２のマイクロホンアレイの場合、７本のマイクロホンを用いて各時刻tlで、基準点における信号f0,l およびその各方向の空間差分
【数９】
0006529451.tif 000010
0006529451.tif 000011
0006529451.tif 000012
を得ることができる。ただし、lは離散時刻のインデックスを表す。
【００１９】
基準点における観測信号の時間微分を時間差分で近似することにすると、式(9) は
【数１０】
0006529451.tif 000013
と表せる。ただし、nx、ny、nz はそれぞれ のx, y, z 方向の成分、T はサンプリング周期である。
【００２０】
式(10) の左辺を右辺に移項すると
【数１１】
0006529451.tif 000014
が得られる。ここで、f0,l, fx,l, fy,l, fz,l を窓関数で窓掛けして取得された信号とする。切り出し区間の両端点の影響を無視できるものとすると、式(11) は周波数領域で
【数１２】
0006529451.tif 000015
と表される。ただし、F0,m, Fx,m, Fy,m, Fz,m はf0,m, fx,m, fy,m, fz,m の離散Fourier 変換であり、m は離散周波数インデックスである。
【００２１】
式(12) の右辺は雑音の存在や差分近似に伴う誤差により実際には必ずしも厳密に0 にはならない。
【００２２】
そこで、式(11) の右辺を
【数１３】
0006529451.tif 000016
のように誤差変数εx,m, εy,m, εz,m に置き換え、これらを平均が0 で互いに独立な正規確率変数（複素正規分布に従う確率変数）
【数１４】
0006529451.tif 000017
と仮定する。また、観測点における観測信号の各周波数成分を平均が0、分散がσ20,m の正規確率変数とする。これは、
【数１５】
0006529451.tif 000018
と仮定することに相当する。
【００２３】
ここで、Fx,m, Fy,m, Fz,m, F0,m を並べたベクトルとεx,m, εy,m, εz,m, ε0,m を並べたベクトルを
【数１６】
0006529451.tif 000019
0006529451.tif 000020
とし、f0,..., fL/2 を連結したベクトルとε0,...,εL/2を連結したベクトルを
【数１７】
0006529451.tif 000021
と表記すると、式(13) は
【数１８】
0006529451.tif 000022
の形で書ける。ただし、θ= ｛R,n｝ であり、A(θ) は
【数１９】
0006529451.tif 000023
で与えられる。式(14), (16) より、εは平均が0、分散共分散行列が
【数２０】
0006529451.tif 000024
の複素正規分布
【数２１】
0006529451.tif 000025
に従う。
【００２４】
【数２２】
0006529451.tif 000026
（A(θ) は正則）であるので、f は
【数２３】
0006529451.tif 000027
と表され、式(24) より、
【数２４】
0006529451.tif 000028
が言える。従って、観測信号およびその空間差分が与えられた下での最尤音源位置^θは
【数２５】
0006529451.tif 000029
により得られる。
【００２５】
＜複数音源の定位アルゴリズム＞
以上のf の確率モデル化により、音源が複数個存在する場合、および雑音が存在する場合の観測信号の確率分布を導くことができる。音源インデックスをk とし、音源k に由来する観測信号の成分、音源位置パラメータをそれぞれf(k)、θ(k) とする。また、f(k) の周波数m の成分エネルギーをσ(k)0,m2 とする。式(29) より、
【数２６】
0006529451.tif 000030
となる。また、加法雑音をｖとし、観測信号を
【数２７】
0006529451.tif 000031
とする。f(1),・・・,f(K),ｖが互いに独立であれば、観測信号y は
【数２８】
0006529451.tif 000032
に従う。ただし、Γはｖの分散共分散行列である。以上より、複数の音源と雑音が存在する場合の各音源の最尤音源位置
【数２９】
0006529451.tif 000033
は、観測信号
【数３０】
0006529451.tif 000034
が与えられた下で
【数３１】
0006529451.tif 000035
を解くことにより得られる。
【００２６】
y を不完全データ、
【数３２】
0006529451.tif 000036
を完全データと見なすことで、以上の最尤推定問題に対しExpectation-Maximization (EM) アルゴリズムを適用することができる。完全データ対数尤度log p(x｜θ) は
【数３３】
0006529451.tif 000037
0006529451.tif 000038
で与えられるので、y が与えられた下でのlog p(x｜θ) のx に関する条件付き期待値（Q 関数）は
【数３４】
0006529451.tif 000039
で与えられる。ただし、
【数３５】
0006529451.tif 000040
はx に関係する項のみについての等号を意味する。この関数が増大するようにθを更新するステップ（M ステップ）と、更新したθをθ´に代入し、
【数３６】
0006529451.tif 000041
と
【数３７】
0006529451.tif 000042
を計算するステップ（E ステップ）を繰り返すことでp(y｜θ) を局所的に最大にするθを求めることができる。
【００２７】
完全データx と不完全データy の関係は
【数３８】
0006529451.tif 000043
と書けるので、
【数３９】
0006529451.tif 000044
はそれぞれ
【数４０】
0006529451.tif 000045
で与えられる。以上より、以下の初期設定、Ｅステップ、Ｍステップからなるアルゴリズムを得る。
【００２８】
（初期ステップ）
θを初期設定する。
【００２９】
（E ステップ）
θをθ´に代入し、式(39) により
【数４１】
0006529451.tif 000046
を計算する。
【００３０】
（Ｍステップ）
下式によりθを更新する。
【００３１】
【数４２】
0006529451.tif 000047
【００３２】
＜Mステップ更新式＞
A(θ(k)) はC1(θ(k)),...,CL(θ(k)) を対角に並べたブロック対角行列なので、
【数４３】
0006529451.tif 000048
と書かれる。ただし、
【数４４】
0006529451.tif 000049
である。また、Φ(k)m はΦ(k)の4×4のブロック対角成分である。
【００３３】
M ステップではQ(θ,θ´) ができるだけ大きくなるようにn(k), R(k),σ(k)0,m2,Γ を更新する。
【００３４】
Q(θ,θ´) を最大にするn(k), R(k),σ(k)0,m2,Γの同時最適解を解析的に求めることは難しいが、座標勾配法によりそれぞれの変数に関してQ(θ,θ´) が最大となるように反復更新することでp(y｜θ) を局所最大化することができる（EM アルゴリズムでは、M ステップで補助関数が単調に増大することが保証されていれば収束性は保証される）。以下に、M ステップの更新方法を1 例示す。
【００３５】
＜n(k)の更新式＞
n(k)は単位ベクトルなので、
【数４５】
0006529451.tif 000050
の下で
【数４６】
0006529451.tif 000051
ができるだけ小さくなるようにn(k) を更新する。この制約つき最適化問題は、
【数４７】
0006529451.tif 000052
のようなラグランジアンを用いてLagrange 未定乗数法で解くことができる。i 行j 列目の要素のみが1 で残りは0 であるような4×4 行列をEi,jとすると、Cm(θ(k)) は
【数４８】
0006529451.tif 000053
のようにn(k)に依存する項とそうでない項に分解できるので、Ｌ(n(k)) のn(k) に関する偏微分
【数４９】
0006529451.tif 000054
を0 と置くことにより、
【数５０】
0006529451.tif 000055
を得る。ただし、[・]i,j は行列のi 行j 列目の成分を表す。
【００３６】
【数５１】
0006529451.tif 000056
より、あとはn(k)x2＋n(k)y2＋n(k)z2= 1 となるようにγ(k)を二分法などで探索し、式(52)〜(54) に代入すれば良い。
【００３７】
＜音源距離R(k)の更新式＞
【数５２】
0006529451.tif 000057
とする。上記と同様、Cm(θ(k)) は
【数５３】
0006529451.tif 000058
のようにρ(k) に依存する項とそうでない項に分解することができるので、Q(θ,θ´) のρ(k)に関する偏微分
【数５４】
0006529451.tif 000059
を0 と置くことにより、
【数５５】
0006529451.tif 000060
を得る。行列要素ごとの表記にすると
【数５６】
0006529451.tif 000061
となる。ただし、
【数５７】
0006529451.tif 000062
である。式(60) の分子における
【数５８】
0006529451.tif 000063
はFast Fourier Transform (FFT) を用いて効率的に計算することができる。
【００３８】
＜σ(k)0,m2 の更新式＞
上記と同様、Σ(k)-1mを
【数５９】
0006529451.tif 000064
のようにσ(k)0,m2 に依存する項とそうでない項に分解できるので、Q(θ,θ´) のσ(k)0,m2 に関する偏微分を0 と置くことにより、
【数６０】
0006529451.tif 000065
を得る。
【００３９】
＜雑音分散共分散行列Γの更新式＞
雑音の分散共分散行列を
【数６１】
0006529451.tif 000066
のように、正規化分散共分散行列モデルWm と周波数mの成分エネルギーν2mの積で表し、ν2mを推定すべき変数とする。後述するがWm は空間無相関モデルや拡散音場モデルなどから導かれる定数行列である。Q(θ,θ´) のν2 に関する偏微分を0 と置くことにより、
【数６２】
0006529451.tif 000067
を得る。
【００４０】
＜正規化分散共分散行列モデルW の設定方法＞
ここでは雑音の空間相関行列から正規化分散共分散行列モデルWm の設定例を述べる。図2 のような7 本のマイクロホンの配置を想定する。ここで，fi,0,...,fi,L-1のFourier 変換をFi,0,...,Fi,L-1とする。~fm = (F0,m,...,F6,m)Tおよびfm = (Fx,m, Fy,m, Fz,m, F0,m)T の関係は
【数６３】
0006529451.tif 000068
と書かれることから、~fm= (F0,m,...,F6,m)T の分散共分散行列をΨm とすると、fmの分散共分散行列はBΨmBT となる。従って、例えば空間的に無相関で等しいパワーの雑音を仮定する場合、Ψm は単位行列となるため、fmの分散共分散行列Wm を
【数６４】
0006529451.tif 000069
と置けば良い。
【００４１】
ある区域内で、エネルギー密度が一様でかつすべての方向に対するエネルギーの流れが等しい確率であるとみなせる分布をしている音場を拡散音場といい、残響環境の音場を良く近似的に表すことが知られている。拡散音場においては、２点間の空間相関係数が距離d にのみ依存し、
【数６５】
0006529451.tif 000070
で与えられる。従って、拡散性雑音を仮定する場合、図2 のようなアレイ幾何の例では、~fm =(F0,m,...,F6,m)Tの分散共分散行列Ψm は
【数６６】
0006529451.tif 000071
となる。これを用いて、fmの分散共分散行列Wm をBΨmBT と置けば良い。
【００４２】
＜システム構成＞
次に、マイクロホンアレイにより入力された音響信号から、複数の音源の位置を推定する音源定位装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【００４３】
図３に示すように、本発明の実施の形態に係る音源定位装置１００は、ＣＰＵと、ＲＡＭと、音源定位処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【００４４】
図３に示すように、音源定位装置１００は、入力部１０と、演算部２０と、出力部９０とを備えている。
【００４５】
入力部１０は、上記図２に示すようなマイクロホンアレイの各マイクロホンから出力された、複数の音源からの音源信号が混じっている音響信号（以後、観測信号）の時系列データを受け付ける。
【００４６】
演算部２０は、空間差分算出部２２と、時間周波数展開部２４と、音源位置推定部２５と、を含んで構成されている。
【００４７】
空間差分算出部２２は、マイクロホンアレイの各マイクロホンから出力された観測信号から、各時刻tlで、基準点のマイクロホンにおける観測信号f0,lを取得すると共に、以下の式に従って、各方向ｘ、ｙ、ｚの空間差分ｆx,l，ｆy,l，ｆz,lを算出する。
【００４８】
【数６７】
0006529451.tif 000072
0006529451.tif 000073
0006529451.tif 000074
【００４９】
時間周波数展開部２４は、空間差分算出部２２により得られた、基準点のマイクロホンにおける各時刻ｔlの観測信号f0,lから、各周波数ｍの観測時間周波数成分Ｆ0,mを計算する。また、時間周波数展開部２４は、空間差分算出部２２により得られた、各時刻ｔlの各方向ｘ、ｙ、ｚの空間差分ｆx,l，ｆy,l，ｆz,lから、各周波数ｍの観測時間周波数成分Ｆx,m，Ｆy,m，Ｆz,mを計算する。本実施の形態においては、短時間フーリエ変換やウェーブレット変換などの時間周波数展開を行う。
【００５０】
音源位置推定部２５は、時間周波数展開部２４において取得した各周波数ｍの観測時間周波数成分Ｆx,m，Ｆy,m，Ｆz,m，Ｆ0,mからなる観測周波数成分ｙに基づいて、ＥＭアルゴリズムを用いて、音源拘束偏微分方程式の周波数領域表現を用いて定められた、複数の音源と加法雑音が存在する場合における、複数の音源の各々の位置を条件とした、観測周波数成分ｙの確率分布Ｐ（ｙ｜θ）を大きくするように、複数の音源の各々の位置を推定する。
【００５１】
音源位置推定部２５は、期待値算出部２６と、変数更新部２８と、収束判定部３０とを備えている。
【００５２】
期待値算出部２６は、時間周波数展開部２４において取得した各周波数ｍの観測時間周波数成分Ｆx,m，Ｆy,m，Ｆz,m，Ｆ0,mからなる観測周波数成分ｙと、初期設定された、又は前回更新された各音源ｋの音源位置θ(k)と、初期設定された、又は前回更新された各音源ｋの音源信号の成分エネルギーσ(k)0,m2に基づいて、上記式(39) により
【数６８】
0006529451.tif 000075
を計算する。
【００５３】
変数更新部２８は、期待値算出部２６によって算出された
【数６９】
0006529451.tif 000076
に基づいて、Q(θ,θ´) ができるだけ大きくなるように、上記式（５２）〜式（５４）、式（６０）、式（６４）、式（６６）、式（６７）に従って、各音源ｋの方向ベクトルn(k),音源距離R(k),成分エネルギーσ(k)0,m2,雑音共分散行列Γ を更新する。
【００５４】
収束判定部３０は、予め定められた収束判定条件を満たすまで、期待値算出部２６及び変数更新部２８による各処理を繰り返させる。収束判定条件としては、例えば、予め定められた繰り返し回数に到達することである。
【００５５】
収束判定条件を満たしたときに、最終的に得られた各音源ｋの方向ベクトルn(k)、音源距離R(k)を、各音源ｋの位置の推定結果として、出力部９０により出力する。
【００５６】
＜音源定位装置の作用＞
次に、本実施の形態に係る音源定位装置１００の作用について説明する。
【００５７】
入力部１０において、マイクロホンアレイの各マイクロホンから出力された観測信号の時系列データを受け付けると、音源定位装置１００は、図４に示す音源定位処理ルーチンを実行する。
【００５８】
まず、ステップＳ１２０では、マイクロホンアレイの各マイクロホンから入力された観測信号の時系列データから、各時刻tlで、基準点のマイクロホンにおける観測信号f0,lを取得すると共に、各方向ｘ、ｙ、ｚの空間差分ｆx,l，ｆy,l，ｆz,lを算出する。
【００５９】
ステップＳ１２１では、上記ステップＳ１２０で得られた基準点のマイクロホンにおける各時刻ｔlの観測信号f0,lから、各周波数ｍの観測時間周波数成分Ｆ0,mを計算する。また、各時刻ｔlの各方向ｘ、ｙ、ｚの空間差分ｆx,l，ｆy,l，ｆz,lから、各周波数ｍの観測時間周波数成分Ｆx,m，Ｆy,m，Ｆz,mを計算する。
【００６０】
ステップＳ１２２では、各音源ｋの音源位置θ(k)と各音源ｋの音源信号の成分エネルギーσ(k)0,m2とに初期値を設定する。
【００６１】
そして、ステップＳ１２３では、上記ステップＳ１２１で取得した各周波数ｍの観測時間周波数成分Ｆx,m，Ｆy,m，Ｆz,m，Ｆ0,mからなる観測周波数成分ｙと、上記ステップＳ１２２で初期設定された、又は後述するステップＳ１２４で前回更新された各音源ｋの音源位置θ(k)及び音源信号の成分エネルギーσ(k)0,m2とに基づいて、上記式(39) により
【数７０】
0006529451.tif 000077
を計算する。
【００６２】
ステップＳ１２４では、上記ステップＳ１２３で算出された
【数７１】
0006529451.tif 000078
に基づいて、Q(θ,θ´) ができるだけ大きくなるように、上記式（５２）〜式（５４）、式（６０）、式（６４）、式（６６）、式（６７）に従って、各音源ｋの方向ベクトルn(k),音源距離R(k),成分エネルギーσ(k)0,m2,雑音共分散行列Γ を更新する。
【００６３】
ステップＳ１２５において、予め定められた収束判定条件を満たしたか否かを判定し、収束判定条件を満たしていない場合には、上記ステップＳ１２３へ戻る。一方、収束判定条件を満たした場合には、ステップＳ１２６へ進む。
【００６４】
ステップＳ１２６では、上記ステップＳ１２４で最終的に得られた各音源ｋの方向ベクトルn(k),音源距離R(k)を、各音源ｋの位置の推定結果として、出力部９０により出力して、音源定位処理ルーチンを終了する。
【００６５】
＜実験＞
以下の条件で残響環境下の音源定位実験を行った。
【００６６】
（実験条件）
音源数: 1(+ 拡散雑音)
音源位置: マイク中心＋ [-1.73; 1.0; 0.0], マイク中心＋[2.0; 2.0; 0.0], マイク中心＋[0.0; -2.0; 0.0]
部屋サイズ: [6.0; 10.0; 8.0] (中心にマイクを配置)
壁面の反射係数: 0.01, 0.5, 0.8（残響の影響の大きさに相当）
マイク間隔: 0.01, 0.1 [m]
フレーム幅: 16, 32, 64 [点]
実験フレーム数: 2x10(ファイル) (無音区間は含まれていない)
【００６７】
マイク位置が３か所あるので、1 条件当たり実質60 回のデータとなる。
【００６８】
評価指標として、誤差の二乗和平方根(rad)を用いた。
【００６９】
図５〜７に実験結果を示す。単一音源のみの存在を仮定した尤度関数（式(28)）を最大にする方法（従来法に相当し，図中の”OneSrc" はこの方法を意味する。）と比べ、高精度な定位が行えていることが分かった。
【００７０】
以上説明したように、本実施の形態に係る音源定位装置によれば、音源拘束偏微分方程式の周波数領域表現を用いて定められた、複数の音源と加法雑音が存在する場合における、複数の音源の各々の位置を条件とした、基準のマイクロホンの各周波数の観測時間周波数成分、及び複数の方向の各々に対する各周波数の観測時間周波数成分の確率分布を大きくするように、複数の音源の各々の位置を推定することにより、雑音が存在する場合であっても、複数の音源を同時に定位することができる。
【００７１】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【００７２】
例えば、上述の音源定位装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【００７３】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【００７４】
１０   入力部
２０   演算部
２２   空間差分算出部
２４   時間周波数展開部
２５   音源位置推定部
２６   期待値算出部
２８   変数更新部
３０   収束判定部
９０   出力部
１００ 音源定位装置
【図１】
0006529451.tif 000085
【図２】
0006529451.tif 000086
【図３】
0006529451.tif 000087
【図４】
0006529451.tif 000088
【図５】
0006529451.tif 000089
【図６】
0006529451.tif 000090
【図７】
0006529451.tif 000091