(57)【要約】【課題】音源数が未知の場合であっても、複数の音源を同時に定位することができるようにする。【解決手段】音源位置推定部２５が、観測データＹが与えられたときの、複数の音源ｋの位置を表す変数Ｎ及び変数ρ、各周波数に対する各時刻において支配的となる音源のインデックスを示すインジケータを表す変数Ｚ、各音源ｋが支配的になる確率πkを定めるための変数Ｖ、複数の方向の各々に対する各周波数の観測時間周波数成分の分散を表す変数λ、並びに各音源ｋ及び雑音の各周波数の時間周波数成分のパワーを表す変数ζ及び変数γを含む未知変数Θの事後分布ｐ（Θ｜Ｙ）と変関数ｑ（Θ）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき目的関数を最小化するように、変数Ｎ、変数ρ、変数Ｚ、変数Ｖ、変数λ、変数ζ、変数γの各々の分布を表す各パラメータを推定し、各音源ｋの位置を推定する。【選択図】図３
2018142917.tif 000002
【特許請求の範囲】
【請求項１】
マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置であって、
複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出する空間差分算出部と、
前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号を入力として、各周波数の観測時間周波数成分を出力すると共に、前記空間差分算出部によって前記複数の方向の各々に対して算出された前記観測信号の差分を入力として、前記複数の方向の各々に対して、各周波数の観測時間周波数成分を出力する時間周波数展開部と、
前記時間周波数展開部により出力された、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分に基づいて、
前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分からなる観測データＹが与えられたときの、複数の音源ｋの位置を表す変数Ｎ及び変数ρ、各周波数に対する各時刻において支配的となる音源のインデックスを示すインジケータを表す変数Ｚ、各音源ｋが支配的になる確率πkを定めるための変数Ｖ、前記複数の方向の各々に対する各周波数の観測時間周波数成分の分散を表す変数λ、並びに各音源ｋ及び雑音の各周波数の時間周波数成分のパワーを表す変数ζ及び変数γを含む未知変数Θの事後分布ｐ（Θ｜Ｙ）と変関数ｑ（Θ）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき前記目的関数を最小化するように、前記変数Ｎ、前記変数ρ、前記変数Ｚ、前記変数Ｖ、前記変数λ、前記変数ζ、前記変数γの各々の分布を表す各パラメータを推定するパラメータ推定部と、
推定された各音源ｋの位置を表す変数Ｎの分布を表すパラメータ、及び変数ρの分布を表すパラメータに基づいて、各音源ｋの位置を推定する音源位置推定部と、
を含む音源定位装置。
【請求項２】
前記確率πkは、ｋが大きいほど指数的に小さくなるように定められる請求項１記載の音源定位装置。
【請求項３】
前記変関数ｑ（Θ）を、ｑ（Ｎ）ｑ（ρ）ｑ（Ｚ）ｑ（Ｖ）ｑ（λ）ｑ（ζ）ｑ（γ）で近似し、
音源数の打ち切り数Ｋ*を予め定め、ｚω,lを、周波数ｍに対する各時刻ｌにおいて支配的となる音源のインデックスを示すインジケータとし、
Ｋ*＋１以上となるｋ’についてｑ（ｚω,l＝ｋ’）＝０とした請求項１又は２記載の音源定位装置。
【請求項４】
マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置における音源定位方法であって、
空間差分算出部が、複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出し、
時間周波数展開部が、前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号を入力として、各周波数の観測時間周波数成分を出力すると共に、前記空間差分算出部によって前記複数の方向の各々に対して算出された前記観測信号の差分を入力として、前記複数の方向の各々に対して、各周波数の観測時間周波数成分を出力し、
パラメータ推定部が、前記時間周波数展開部により出力された、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分に基づいて、
前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分からなる観測データＹが与えられたときの、複数の音源ｋの位置を表す変数Ｎ及び変数ρ、各周波数に対する各時刻において支配的となる音源のインデックスを示すインジケータを表す変数Ｚ、各音源ｋが支配的になる確率πkを定めるための変数Ｖ、前記複数の方向の各々に対する各周波数の観測時間周波数成分の分散を表す変数λ、並びに各音源ｋ及び雑音の各周波数の時間周波数成分のパワーを表す変数ζ及び変数γを含む未知変数Θの事後分布ｐ（Θ｜Ｙ）と変関数ｑ（Θ）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき前記目的関数を最小化するように、前記変数Ｎ、前記変数ρ、前記変数Ｚ、前記変数Ｖ、前記変数λ、前記変数ζ、前記変数γの各々の分布を表す各パラメータを推定し、
音源位置推定部が、推定された各音源ｋの位置を表す変数Ｎの分布を表すパラメータ、及び変数ρの分布を表すパラメータに基づいて、各音源ｋの位置を推定する
音源定位方法。
【請求項５】
前記確率πkは、ｋが大きいほど指数的に小さくなるように定められる請求項４記載の音源定位方法。
【請求項６】
前記変関数ｑ（Θ）を、ｑ（Ｎ）ｑ（ρ）ｑ（Ｚ）ｑ（Ｖ）ｑ（λ）ｑ（ζ）ｑ（γ）で近似し、
音源数の打ち切り数Ｋ*を予め定め、ｚω,lを、周波数ｍに対する各時刻ｌにおいて支配的となる音源のインデックスを示すインジケータとし、
Ｋ*＋１以上となるｋ’についてｑ（ｚω,l＝ｋ’）＝０とした請求項４又は５記載の音源定位方法。
【請求項７】
請求項１〜請求項３の何れか１項に記載の音源定位装置の各部としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音源定位装置、方法、及びプログラムに係り、特に、音響信号から、音源の位置を推定する音源定位装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
波源定位は、レーダやソナーといった幅広い応用を有している。特に、小さいアレイで、移動する波源を瞬時に定位し追跡できるようにすることは重要課題である。波源定位問題に対する従来法としては、Multiple Signal Classication (MUSIC) 法、Generalized Cross-Correlation methods with Phase Transform (GCC-PHAT) 法、波源拘束偏微分方程式に基づく手法などがある。
【０００３】
MUSIC 法やGCC-PHAT 法は、音源に対し平面波を仮定し各音源のセンサ間での到来時間差を定位の手がかりとするため、一般にアレイサイズは大きい方が有利となる。また、いずれもセンサアレイの受信信号間の自己相関関数や相互相関関数といった、統計量に基づく手法であるため、音源を高い精度で定位するためには観測時間幅を十分長く取る必要がある。このため、これらの手法は小さいアレイサイズと瞬時的な観測のみによる波源定位には必ずしも向いていない。一方、波源拘束偏微分方程式に基づく手法は、各時刻ごとに成立する音響信号の時空間偏微分方程式を元に音源定位を行うもので、理論的には瞬時の小領域観測のみで波源定位を行うことが可能である。
【０００４】
ただし、この手法は単一波源に対して成立する方程式をベースとしているため、雑音や複数の点音源が存在する場合などのように観測音響信号が偏微分方程式から逸脱する場合に脆弱であるという欠点を有している。本枠組においてランダム雑音や複数の点音源が存在する場合への拡張を可能にするため、波源拘束偏微分方程式を基にしたアレイ観測信号の確率モデルとそれに基づく波源定位アルゴリズムが非特許文献１で提案されている。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】鈴木惇, 亀岡弘和,"波源拘束差分方程式に基づく音響信号の確率モデル化と複数音源定位アルゴリズム",日本音響学会講演論文集, pp.615-618, 2016年3月
【発明の概要】
【発明が解決しようとする課題】
【０００６】
上記非特許文献１の手法では音源数を仮定する必要があるが、実環境では音源数が未知の場合が多い。従って、音源数を仮定せずともモデルの複雑度を適応させながら実際の音源数に合わせて音源定位を行えるようにすることが望ましい。
【０００７】
本発明は、上記事情を鑑みてなされたものであり、音源数が未知の場合であっても、複数の音源を同時に定位することができる音源定位装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記の目的を達成するために本発明に係る音源定位装置は、マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置であって、複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出する空間差分算出部と、前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号を入力として、各周波数の観測時間周波数成分を出力すると共に、前記空間差分算出部によって前記複数の方向の各々に対して算出された前記観測信号の差分を入力として、前記複数の方向の各々に対して、各周波数の観測時間周波数成分を出力する時間周波数展開部と、前記時間周波数展開部により出力された、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分に基づいて、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分からなる観測データＹが与えられたときの、複数の音源ｋの位置を表す変数Ｎ及び変数ρ、各周波数に対する各時刻において支配的となる音源のインデックスを示すインジケータを表す変数Ｚ、
各音源ｋが支配的になる確率πkを定めるための変数Ｖ、前記複数の方向の各々に対する各周波数の観測時間周波数成分の分散を表す変数λ、並びに各音源ｋ及び雑音の各周波数の時間周波数成分のパワーを表す変数ζ及び変数γを含む未知変数Θの事後分布ｐ（Θ｜Ｙ）と変関数ｑ（Θ）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき前記目的関数を最小化するように、前記変数Ｎ、前記変数ρ、前記変数Ｚ、前記変数Ｖ、変数λ、変数ζ、前記変数γの各々の分布を表す各パラメータを推定するパラメータ推定部と、推定された各音源ｋの位置を表す変数Ｎの分布を表すパラメータ、及び変数ρの分布を表すパラメータに基づいて、各音源ｋの位置を推定する音源位置推定部と、を含んで構成されている。
【０００９】
本発明に係る音源定位方法は、マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置における音源定位方法であって、空間差分算出部が、複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出し、時間周波数展開部が、前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号を入力として、各周波数の観測時間周波数成分を出力すると共に、前記空間差分算出部によって前記複数の方向の各々に対して算出された前記観測信号の差分を入力として、前記複数の方向の各々に対して、各周波数の観測時間周波数成分を出力し、パラメータ推定部が、前記時間周波数展開部により出力された、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分に基づいて、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分からなる観測データＹが与えられたときの、複数の音源ｋの位置を表す変数Ｎ及び変数ρ、各周波数に対する各時刻において支配的となる音源のインデックスを示すインジケータを表す変数Ｚ、各音源ｋが支配的になる確率πkを定めるための変数Ｖ、前記複数の方向の各々に対する各周波数の観測時間周波数成分の分散を表す変数λ、並びに各音源ｋ及び雑音の各周波数の時間周波数成分のパワーを表す変数ζ及び変数γを含む未知変数Θの事後分布ｐ（Θ｜Ｙ）と変関数ｑ（Θ）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき前記目的関数を最小化するように、前記変数Ｎ、前記変数ρ、前記変数Ｚ、前記変数Ｖ、変数λ、変数ζ、前記変数γの各々の分布を表す各パラメータを推定し、音源位置推定部が、推定された各音源ｋの位置を表す変数Ｎの分布を表すパラメータ、及び変数ρの分布を表すパラメータに基づいて、各音源ｋの位置を推定する。
【００１０】
本発明に係るプログラムは、上記の音源定位装置の各部としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００１１】
以上説明したように、本発明の音源定位装置、方法、及びプログラムによれば、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分からなる観測データＹが与えられたときの、複数の音源ｋの位置を表す変数Ｎ及び変数ρ、各周波数に対する各時刻において支配的となる音源のインデックスを示すインジケータを表す変数Ｚ、各音源ｋが支配的になる確率πkを定めるための変数Ｖ、前記複数の方向の各々に対する各周波数の観測時間周波数成分の分散を表す変数λ、並びに各音源ｋ及び雑音の各周波数の時間周波数成分のパワーを表す変数ζ及び変数γを含む未知変数Θの事後分布ｐ（Θ｜Ｙ）と変関数ｑ（Θ）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき前記目的関数を最小化するように、前記変数Ｎ、前記変数ρ、前記変数Ｚ、前記変数Ｖ、変数λ、変数ζ、前記変数γの各々の分布を表す各パラメータを推定し、各音源ｋの位置を推定することにより、音源数が未知の場合であっても、複数の音源を同時に定位することができる、という効果が得られる。
【図面の簡単な説明】
【００１２】
【図１】点音源から観測点ｒへ到来する球面波を示す図である。
【図２】マイクロホンアレイの配置の一例を示す図である。
【図３】本発明の実施の形態に係る音源定位装置の構成を示す概略図である。
【図４】本発明の実施の形態に係る音源定位装置における音源定位処理ルーチンの内容を示すフローチャートである。
【図５】実験における音源位置とマイクロホン位置を示す図である。
【図６】残響時間0.5 s、観測長2779 ms の条件での実験結果を示す図である。
【図７】残響時間0.5 s、観測長1665 ms の条件での実験結果を示す図である。
【図８】残響時間0.2 s、観測長2779 ms の条件での実験結果を示す図である。
【図９】残響時間0.2 s、観測長1665 ms の条件での実験結果を示す図である。
【発明を実施するための形態】
【００１３】
以下、図面を参照して本発明の実施の形態を詳細に説明する。本発明で提案する技術は、音響信号から波源位置を推定することを目的とした信号処理技術である。
【００１４】
＜本発明の実施の形態の概要＞
本発明の実施の形態は、上述した従来手法の利点を併せ持つ、小領域・瞬時観測による複数音源の波源定位を可能にする技術である。
【００１５】
本発明の実施の形態では、波源拘束偏微分方程式の時間周波数領域表現をベースにした音響信号の確率分布と、雑音を含む全音源のスパース性の仮定（複数の音源が混在する音響信号の時間周波数表現において、各時間周波数点で高々一つの音源のみが支配的であるという仮定）に基づき、変分推論アルゴリズムにより各時間周波数点でどの音源が支配的らしいかを推定しながら各音源の波源定位を行う。
【００１６】
さらに、本発明の実施の形態では、Dirichlet過程混合モデルをヒントにしたアイディアにより、音源数を仮定せずともモデルの複雑度を適応させながら実際の音源数に合わせて音源定位を行うことができる。
【００１７】
＜本発明の実施の形態の原理＞
次に、音源の位置を推定する原理について説明する。
【００１８】
＜音源拘束偏微分方程式＞
図１に示すように、観測点の基準となる位置ベクトルを
2018142917.tif 000003
とし、単一波源の位置ベクトルを
2018142917.tif 000004
とする。波源の信号をg(t)、音速をc とし、単一点波源からの球面波伝播を仮定すると観測点における観測値は
【００１９】
【数１】
2018142917.tif 000005
【００２０】
と表される。ここで、
【００２１】
【数２】
2018142917.tif 000006
【００２２】
である。観測点から波源方向へ向かう単位ベクトルを
2018142917.tif 000007
とすると、
【００２３】
【数３】
2018142917.tif 000008
2018142917.tif 000009
【００２４】
であるため、
2018142917.tif 000010
の空間微分は、
【００２５】
【数４】
2018142917.tif 000011
2018142917.tif 000012
2018142917.tif 000013
【００２６】
となる。また、
2018142917.tif 000014
の時間微分は
【００２７】
【数５】
2018142917.tif 000015
【００２８】
となるので、式(1) と式(8) を式(7) に代入することでg が消去され、
【００２９】
【数６】
2018142917.tif 000016
【００３０】
のように、観測信号とその時間・空間微分のみを含む方程式を立てることができる。ただし、
2018142917.tif 000017
は観測点から波源までの距離である。この式を音源拘束式と呼ぶ。以上のように音源拘束式は、任意の音源信号波形で成り立つ、音源の位置と空間の場の一意な関係を記述する偏微分方程式である。
【００３１】
＜音源拘束式に基づく音響信号の確率モデル化＞
観測信号の空間微分を空間差分で近似するため、以下では図２に示すマイクロホンアレイを仮定する。ただし、マイクロホンの配置は観測信号の空間微分を空間差分で近似できるものであれば良く、以下の理論は図２の配置に限らない。図２のマイクロホンアレイの場合、7 本のマイクロホンを用いて各時刻tj で、基準点における信号f0,jおよびその各方向の空間差分
2018142917.tif 000018
2018142917.tif 000019
2018142917.tif 000020
を得ることができる。ただし、j は標本時刻のインデックスを表す。
【００３２】
このとき式(9) は
【００３３】
【数７】
2018142917.tif 000021
【００３４】
と表せる。ただし、i = x,y,z で、nx、 ny、 nzはそれぞれn のx、 y、 z 方向の成分である。式(10)の左辺を右辺に移項すると
【００３５】
【数８】
2018142917.tif 000022
【００３６】
が得られる。ここで、f0,j、fx,j、fy,j、fz,jを窓関数で窓掛けして取得された信号とする。切り出し区間の両端点の影響を無視できるものとすると、式(11) は周波数領域で
【００３７】
【数９】
2018142917.tif 000023
【００３８】
と表される。ただし、F0,m、 Fx,m、 Fy,m、 Fz,mはf0,j、fx,j、fy,j、fz,jの離散Fourier 変換であり、mは離散周波数インデックスである。
【００３９】
式(12) の右辺は差分近似に伴う誤差により実際には必ずしも厳密に0 にはならない。そこで、式(11) の右辺を
【００４０】
【数１０】
2018142917.tif 000024
【００４１】
のように誤差変数
2018142917.tif 000025
に置き換え、これらを平均が0 で互いに独立な正規確率変数（複素正規分布に従う確率変数）
【００４２】
【数１１】
2018142917.tif 000026
【００４３】
と仮定する。また、観測点における観測信号の各周波数成分を平均が0、分散がσ20,m の正規確率変数とする。これは、
【００４４】
【数１２】
2018142917.tif 000027
2018142917.tif 000028
【００４５】
と仮定することに相当する。
【００４６】
ここで、
2018142917.tif 000029
を並べたベクトルと
2018142917.tif 000030
を並べたベクトルを
【００４７】
【数１３】
2018142917.tif 000031
2018142917.tif 000032
【００４８】
とする。式(13) は
【００４９】
【数１４】
2018142917.tif 000033
【００５０】
の形で書ける。ただし
2018142917.tif 000034
であり、
2018142917.tif 000035
は
【００５１】
【数１５】
2018142917.tif 000036
【００５２】
で与えられる。式(14)、(16) より、
2018142917.tif 000037
は平均が
2018142917.tif 000038
、分散共分散行列
2018142917.tif 000039
が
【００５３】
【数１６】
2018142917.tif 000040
【００５４】
の複素正規分布
【００５５】
【数１７】
2018142917.tif 000041
【００５６】
に従う。
2018142917.tif 000042
であるので、
2018142917.tif 000043
は
【００５７】
【数１８】
2018142917.tif 000044
【００５８】
と表される。よって式(22) より、
2018142917.tif 000045
の確率密度関数
【００５９】
【数１９】
2018142917.tif 000046
【００６０】
を得る。以上より、観測スペクトルおよびその空間差分
2018142917.tif 000047
が与えられた下で、単一音源を定位する問題は、
【００６１】
【数２０】
2018142917.tif 000048
2018142917.tif 000049
【００６２】
を解く最尤推定問題に帰着する。
【００６３】
＜音源のスパース性を活用した複数音源の定位アルゴリズム＞
音声信号や楽音など実世界の音響信号の多くは時間周波数成分がスパースである。従って、複数の音源が同時に混在する場合でも、各時間周波数点では高々一つの音源のみが支配的であると仮定できる場合が多い。この音源の時間周波数成分のスパース性の仮定と以上の
2018142917.tif 000050
の確率モデル化に基づき、音源が複数個存在する場合、および雑音が存在する場合の観測信号の確率分布を導くことができる。
【００６４】
信号の切り出しフレームの時刻のインデックスをl、音源インデックスをk = 0,...,K とし、k = 0は雑音、k ≠0は点音源に対応するものとする。また、点音源kの位置を
2018142917.tif 000051
とする。ここで、雑音を含む全音源の時間周波数成分のスパース性を仮定し、周波数m、時刻l においてzm,l 番目の音源のみが非零のパワーをもち、それ以外の音源のパワーを0 とする。このとき、所与のzm,lの下での観測信号の時間周波数成分とその空間差分
2018142917.tif 000052
（以後、観測信号）の条件付き確率密度関数は
【００６５】
【数２１】
2018142917.tif 000053
【００６６】
で与えられる。ここで、
【００６７】
【数２２】
2018142917.tif 000054
2018142917.tif 000055
【００６８】
である。また、
2018142917.tif 000056
は雑音の時間周波数成分の分散共分散行列で、周波数にのみ依存する正規化分散共分散行列モデル
2018142917.tif 000057
と時刻にも依存する雑音のパワー
2018142917.tif 000058
の積
【００６９】
【数２３】
2018142917.tif 000059
【００７０】
で表されるものとする。拡散性雑音を仮定する場合の
2018142917.tif 000060
の設定方法については後述する。
2018142917.tif 000061
はすべての未知パラメータ
2018142917.tif 000062
を表す。
2018142917.tif 000063
の事前確率を
2018142917.tif 000064
とすると、観測信号
2018142917.tif 000065
の確率密度関数（
2018142917.tif 000066
の尤度関数）は
【００７１】
【数２４】
2018142917.tif 000067
【００７２】
と書ける。以上より、複数の音源と雑音が存在する場合の各音源の位置
2018142917.tif 000068
を推定する問題は、観測信号
2018142917.tif 000069
が与えられた下で
【００７３】
【数２５】
2018142917.tif 000070
【００７４】
を解く最尤推定問題に帰着する。この最適化問題の大域解は解析的に解くことはできないが、
2018142917.tif 000071
を潜在変数としたExpectation-Maximization (EM) アルゴリズムにより停留点を探索することができる。
【００７５】
＜ノンパラメトリックベイズモデリング＞
実環境では音源数が未知の場合が多い。上記の定式化では音源数K を既知と仮定したが、音源数を仮定せずともモデルの複雑度を適応させながら実際の音源数に合わせて音源定位を行えるようにすることが望ましい。そこで、本実施の形態では上記の生成モデルをDirichlet 過程混合モデルに拡張する。k = 0 を雑音、
2018142917.tif 000072
を点音源のインデックスとし、
2018142917.tif 000073
を加算無限次元の離散分布
2018142917.tif 000074
に従って生成される確率変数とする。
【００７６】
【数２６】
2018142917.tif 000075
【００７７】
ここで、点(m, l)において雑音が支配的である（zm,l＝０となる）確率π0を、超パラメータ
2018142917.tif 000076
のベータ分布に従って生成される変数
【００７８】
【数２７】
2018142917.tif 000077
2018142917.tif 000078
【００７９】
とし、点音源k ≠ 0が支配的である（zm,l＝ｋとなる）確率πkを、棒折り過程に従って決まる変数
【００８０】
【数２８】
2018142917.tif 000079
2018142917.tif 000080
【００８１】
とする。以上のプロセスで生成される
2018142917.tif 000081
の期待値は、k が大きいほど指数的に小さくなる傾向を持つため、大きいk に対応した音源ほどアクティブになる確率が低くなる。よって、観測信号からパラメータを推論する際、必要最小限の音源インデックス数の混合モデルで観測信号を説明しようとする効果がもたらされる。以上の生成モデルにおいて、全未知変数Θ は以下となる。
【００８２】
【数２９】
2018142917.tif 000082
2018142917.tif 000083
2018142917.tif 000084
2018142917.tif 000085
2018142917.tif 000086
2018142917.tif 000087
2018142917.tif 000088
2018142917.tif 000089
【００８３】
＜変分推論アルゴリズム＞
以上の生成モデル化により、
【００８４】
【数３０】
2018142917.tif 000090
2018142917.tif 000091
2018142917.tif 000092
【００８５】
と書ける。また、その他の事前分布p(V)、p(N)、p(ρ)、p(λ)、p(ζ) 及びp(γ)を
【００８６】
【数３１】
2018142917.tif 000093
2018142917.tif 000094
2018142917.tif 000095
2018142917.tif 000096
2018142917.tif 000097
【００８７】
と仮定すると、観測信号
2018142917.tif 000098
と未知変数Θの同時分布
【００８８】
【数３２】
2018142917.tif 000099
【００８９】
を具体的に記述することができる。ただし、
2018142917.tif 000100
はvon Mise-Fisher分布、
2018142917.tif 000101
は実正規分布、
2018142917.tif 000102
はガンマ分布を表す。Θの事後分布
2018142917.tif 000103
を解析的に得ることは難しいが、変分推論法に基づき近似分布を反復計算により得ることができる。変分推論は、
2018142917.tif 000104
を満たす非負の変関数
2018142917.tif 000105
を、事後分布
2018142917.tif 000106
との間のKullback-Leibler(KL)ダイバージェンス
【００９０】
【数３３】
2018142917.tif 000107
【００９１】
が小さくなるように求める方法であり、上記式（５５）を目的関数とする。
2018142917.tif 000108
を
【００９２】
【数３４】
2018142917.tif 000109
【００９３】
のように近似できると仮定すると、
2018142917.tif 000110
について反復的に式(55)の目的関数を最小化することで
2018142917.tif 000111
の近似分布を得ることができる。また、
2018142917.tif 000112
に関して、以下の打ち切り近似
【００９４】
【数３５】
2018142917.tif 000113
【００９５】
を行う。この近似はモデルの複雑度（音源数）を固定したことを意味するのではなくｑの関数空間をある範囲に限定したことを意味する。従って、
2018142917.tif 000114
をできるだけ良く近似したければ
2018142917.tif 000115
を大きくして
2018142917.tif 000116
が取りうる範囲を広くすれば良い。
【００９６】
式(55) の各q についての変分を0 と置くことで、以下を得る。これらを変分事後分布更新式という。
【００９７】
【数３６】
2018142917.tif 000117
2018142917.tif 000118
2018142917.tif 000119
2018142917.tif 000120
2018142917.tif 000121
2018142917.tif 000122
2018142917.tif 000123
【００９８】
ただし
2018142917.tif 000124
は
2018142917.tif 000125
の
2018142917.tif 000126
に関する期待値を意味し、X が連続変数の場合は
2018142917.tif 000127
、離散変数の場合は
2018142917.tif 000128
で与えられる。また、
2018142917.tif 000129
は、
2018142917.tif 000130
の中のX 以外のすべての要素の集合を表す。導出は次節で述べるが、それぞれの変分事後分布更新式は以下のような形となる。
【００９９】
【数３７】
2018142917.tif 000131
2018142917.tif 000132
2018142917.tif 000133
2018142917.tif 000134
2018142917.tif 000135
2018142917.tif 000136
2018142917.tif 000137
【０１００】
＜変分更新式の導出＞
＜結合分布＞
上記で立てた生成モデルによりlog p(Θ，Ｙ) は具体的に以下のように書ける。
【０１０１】
【数３８】
2018142917.tif 000138
2018142917.tif 000139
【０１０２】
＜音源方向N の変分事後分布更新式＞
2018142917.tif 000140
の中でN に関係する項は
【０１０３】
【数３９】
2018142917.tif 000141
【０１０４】
であり、
【０１０５】
【数４０】
2018142917.tif 000142
【０１０６】
となるので、N の変分事後分布更新式は
【０１０７】
【数４１】
2018142917.tif 000143
2018142917.tif 000144
【０１０８】
となる。よって、期待値
2018142917.tif 000145
は以下となる。
【０１０９】
【数４２】
2018142917.tif 000146
【０１１０】
＜音源距離（の逆数）ρの変分事後分布更新式＞
2018142917.tif 000147
の中でρに関係する項は
【０１１１】
【数４３】
2018142917.tif 000148
【０１１２】
であり、
【０１１３】
【数４４】
2018142917.tif 000149
【０１１４】
となるので、ρの変分事後分布更新式は、
【０１１５】
【数４５】
2018142917.tif 000150
2018142917.tif 000151
2018142917.tif 000152
【０１１６】
となる。よって、期待値
2018142917.tif 000153
は以下となる。
【０１１７】
【数４６】
2018142917.tif 000154
2018142917.tif 000155
【０１１８】
＜アクティブ音源インデックスZ の変分事後分布更新式＞
2018142917.tif 000156
の中でZ に関係する項は
【０１１９】
【数４７】
2018142917.tif 000157
【０１２０】
であり、
【０１２１】
【数４８】
2018142917.tif 000158
【０１２２】
となる。ただし、
2018142917.tif 000159
は、ｋ＝０のとき、
【０１２３】
【数４９】
2018142917.tif 000160
【０１２４】
ｋ≠０のとき、
【０１２５】
【数５０】
2018142917.tif 000161
【０１２６】
である。従って、Z の変分事後分布更新式は、
【０１２７】
【数５１】
2018142917.tif 000162
2018142917.tif 000163
【０１２８】
となる。よって、期待値
2018142917.tif 000164
は以下となる。
【０１２９】
【数５２】
2018142917.tif 000165
【０１３０】
＜棒折り比V の変分事後分布更新式＞
2018142917.tif 000166
の中でV に関係する項は
【０１３１】
【数５３】
2018142917.tif 000167
2018142917.tif 000168
【０１３２】
であり、
【０１３３】
【数５４】
2018142917.tif 000169
【０１３４】
となるので、V の変分事後分布更新式は、
【０１３５】
【数５５】
2018142917.tif 000170
2018142917.tif 000171
2018142917.tif 000172
【０１３６】
となる。よって、期待値
2018142917.tif 000173
は以下となる。
【０１３７】
【数５６】
2018142917.tif 000174
2018142917.tif 000175
【０１３８】
ただし、
2018142917.tif 000176
はディガンマ関数である。
【０１３９】
＜誤差変数
2018142917.tif 000177
の精度（分散の逆数）λの変分事後分布更新式＞
2018142917.tif 000178
の中でλ に関係する項は
【０１４０】
【数５７】
2018142917.tif 000179
【０１４１】
であり、
【０１４２】
【数５８】
2018142917.tif 000180
【０１４３】
となるので、λ の変分事後分布更新式は、
【０１４４】
【数５９】
2018142917.tif 000181
2018142917.tif 000182
2018142917.tif 000183
【０１４５】
となる。よって、期待値
2018142917.tif 000184
は以下となる．
【０１４６】
【数６０】
2018142917.tif 000185
2018142917.tif 000186
【０１４７】
＜音源パワー（の逆数）ζの変分事後分布更新式＞
2018142917.tif 000187
の中でζに関係する項は
【０１４８】
【数６１】
2018142917.tif 000188
【０１４９】
であり、
【０１５０】
【数６２】
2018142917.tif 000189
【０１５１】
となるので、ζの変分事後分布更新式は、
【０１５２】
【数６３】
2018142917.tif 000190
2018142917.tif 000191
2018142917.tif 000192
【０１５３】
となる。よって、期待値
【０１５４】
【数６４】
2018142917.tif 000193
2018142917.tif 000194
【０１５５】
＜雑音パワー（の逆数）γの変分事後分布更新式＞
2018142917.tif 000195
の中でγに関係する項は
【０１５６】
【数６５】
2018142917.tif 000196
【０１５７】
であり、
【０１５８】
【数６６】
2018142917.tif 000197
【０１５９】
となるので、γの変分事後分布更新式は、
【０１６０】
【数６７】
2018142917.tif 000198
2018142917.tif 000199
2018142917.tif 000200
【０１６１】
となる。よって、期待値
2018142917.tif 000201
は以下となる。
【０１６２】
【数６８】
2018142917.tif 000202
2018142917.tif 000203
【０１６３】
＜雑音分散共分散行列W＞
ここでは、上記図2 のような7 本のマイクロホンの配置を想定した場合の雑音共分散行列Wmの設定例を述べる。ここで、
2018142917.tif 000204
のFourier 変換を
2018142917.tif 000205
とする。
【０１６４】
2018142917.tif 000206
および
2018142917.tif 000207
の関係は
【０１６５】
【数６９】
2018142917.tif 000208
【０１６６】
と書けるため、
2018142917.tif 000209
の分散共分散行列を
2018142917.tif 000210
とすると、
2018142917.tif 000211
の分散共分散行列は
2018142917.tif 000212
となる。従って、例えば空間的に無相関で等しいパワーの雑音を仮定する場合、
2018142917.tif 000213
は単位行列となるため、
2018142917.tif 000214
の分散共分散行列
2018142917.tif 000215
を
【０１６７】
【数７０】
2018142917.tif 000216
【０１６８】
と置けば良い。
【０１６９】
ある区域内で、エネルギー密度が一様でかつすべての方向に対するエネルギーの流れが等しい確率であるとみなせる分布をしている音場を拡散音場といい、残響環境の音場を良く近似的に表すことが知られている。拡散音場においては、2 点間の空間相関係数が距離d にのみ依存し、
【０１７０】
【数７１】
2018142917.tif 000217
【０１７１】
で与えられる。従って、拡散性雑音を仮定する場合、図2 のようなアレイ幾何の例では、
2018142917.tif 000218
の分散共分散行列
2018142917.tif 000219
は
【０１７２】
【数７２】
2018142917.tif 000220
【０１７３】
となる。これを用いて、
2018142917.tif 000221
の分散共分散行列
2018142917.tif 000222
を
2018142917.tif 000223
と置けば良い。このとき、
2018142917.tif 000224
は、
【０１７４】
【数７３】
2018142917.tif 000225
【０１７５】
のような対角行列となる。
【０１７６】
＜変分推論アルゴリズムのまとめ＞
各変数の変分事後分布更新式は
【０１７７】
【数７４】
2018142917.tif 000226
2018142917.tif 000227
2018142917.tif 000228
2018142917.tif 000229
2018142917.tif 000230
2018142917.tif 000231
2018142917.tif 000232
【０１７８】
で与えられ、各分布のパラメータの更新式は以下の通りである。
【０１７９】
【数７５】
2018142917.tif 000233
2018142917.tif 000234
2018142917.tif 000235
2018142917.tif 000236
2018142917.tif 000237
2018142917.tif 000238
2018142917.tif 000239
2018142917.tif 000240
2018142917.tif 000241
2018142917.tif 000242
2018142917.tif 000243
2018142917.tif 000244
2018142917.tif 000245
【０１８０】
ここで
2018142917.tif 000246
は、k = 0 のとき、
【０１８１】
【数７６】
2018142917.tif 000247
【０１８２】
ｋ≠０のとき、
【０１８３】
【数７７】
2018142917.tif 000248
2018142917.tif 000249
【０１８４】
である。また、更新式中に出てくる期待値は以下のように与えられる。
【０１８５】
【数７８】
2018142917.tif 000250
2018142917.tif 000251
2018142917.tif 000252
2018142917.tif 000253
2018142917.tif 000254
2018142917.tif 000255
2018142917.tif 000256
2018142917.tif 000257
2018142917.tif 000258
2018142917.tif 000259
2018142917.tif 000260
2018142917.tif 000261
【０１８６】
＜システム構成＞
次に、マイクロホンアレイにより入力された音響信号から、複数の音源の位置を推定する音源定位装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【０１８７】
図３に示すように、本発明の実施の形態に係る音源定位装置１００は、ＣＰＵと、ＲＡＭと、音源定位処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【０１８８】
図３に示すように、音源定位装置１００は、入力部１０と、演算部２０と、出力部９０とを備えている。
【０１８９】
入力部１０は、上記図２に示すようなマイクロホンアレイの各マイクロホンから出力された、複数の音源からの音源信号が混じっている音響信号（以後、観測信号）の時系列データを受け付ける。
【０１９０】
演算部２０は、空間差分算出部２２と、時間周波数展開部２４と、音源位置推定部２５と、を含んで構成されている。
【０１９１】
空間差分算出部２２は、マイクロホンアレイの各マイクロホンから出力された観測信号から、各時刻tjで、基準点のマイクロホンにおける観測信号f0,jを取得すると共に、以下の式に従って、各方向ｘ、ｙ、ｚの空間差分ｆx,j，ｆy,j，ｆz,jを算出する。
【０１９２】
【数７９】
2018142917.tif 000262
2018142917.tif 000263
2018142917.tif 000264
【０１９３】
時間周波数展開部２４は、空間差分算出部２２により得られた、基準点のマイクロホンにおける各時刻ｔｊの観測信号f0,jから、各周波数ｍの観測時間周波数成分Ｆ0,mを計算する。また、時間周波数展開部２４は、空間差分算出部２２により得られた、各時刻ｔjの各方向ｘ、ｙ、ｚの空間差分ｆx,j，ｆy,j，ｆz,jから、各周波数ｍの観測時間周波数成分Ｆx,m，Ｆy,m，Ｆz,mを計算する。本実施の形態においては、短時間フーリエ変換やウェーブレット変換などの時間周波数展開を行う。
【０１９４】
音源位置推定部２５は、時間周波数展開部２４において取得した各周波数ｍの観測時間周波数成分Ｆx,m，Ｆy,m，Ｆz,m，Ｆ0,mからなる観測周波数成分ｙに基づいて、各周波数ｍの観測時間周波数成分Ｆx,m，Ｆy,m，Ｆz,m，Ｆ0,mからなる観測データＹが与えられたときの、複数の音源ｋの位置を表す変数Ｎ及び変数ρ、各周波数に対する各時刻において支配的となる音源のインデックスを示すインジケータを表す変数Ｚ、各音源ｋが支配的になる確率πkを定めるための変数Ｖ、複数の方向の各々に対する各周波数の観測時間周波数成分の分散を表す変数λ、並びに各音源ｋ及び雑音の各周波数の時間周波数成分のパワーを表す変数ζ及び変数γを含む未知変数Θの事後分布ｐ（Θ｜Ｙ）と変関数ｑ（Θ）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき目的関数を最小化するように、変数Ｎ、変数ρ、変数Ｚ、変数Ｖ、変数λ、変数ζ、変数γの各々の分布を表す各パラメータを推定し、各音源ｋの位置を推定する。
【０１９５】
具体的には、音源位置推定部２５は、変数更新部２８と、収束判定部３０とを備えている。
【０１９６】
変数更新部２８は、まず、変数Ｎ、変数ρ、変数Ｚ、変数Ｖ、変数λ、変数ζ、変数γの各々の分布を表す各パラメータ
2018142917.tif 000265
(以後、変分パラメータと称する)の初期値を設定する。
【０１９７】
また、変数更新部２８は、観測データＹと、変分パラメータ
2018142917.tif 000266
に基づいて、上記式（１３２）〜式（１４６）に従って、変分パラメータ
2018142917.tif 000267
を更新する。
【０１９８】
収束判定部３０は、予め定められた収束判定条件を満足するか否かを判定し、収束判定条件を満足していない場合には、変数更新部２８の処理を繰り返す。収束判定部３０は、収束判定条件を満たしたときに、最終的に得られた
2018142917.tif 000268
に基づいて、各音源ｋの方向ベクトルn(k)、音源距離R(k)を、各音源ｋの位置の推定結果として、出力部９０により出力する。
【０１９９】
収束判定条件としては、反復計算回数が予め定めた回数に達したことを用いればよい。なお、一回のパラメータ更新によるパラメータの変化率がほぼ1になったと見なせたことを、収束判定条件として用いてもよい。
【０２００】
＜音源定位装置の作用＞
次に、本実施の形態に係る音源定位装置１００の作用について説明する。
【０２０１】
入力部１０において、マイクロホンアレイの各マイクロホンから出力された観測信号の時系列データを受け付けると、音源定位装置１００は、図４に示す音源定位処理ルーチンを実行する。
【０２０２】
まず、ステップＳ１２０では、マイクロホンアレイの各マイクロホンから入力された観測信号の時系列データから、各時刻tjで、基準点のマイクロホンにおける観測信号f0,jを取得すると共に、各方向ｘ、ｙ、ｚの空間差分ｆx,j，ｆy,j，ｆz,jを算出する。
【０２０３】
ステップＳ１２１では、上記ステップＳ１２０で得られた基準点のマイクロホンにおける各時刻ｔjの観測信号f0,jから、各周波数ｍの観測時間周波数成分Ｆ0,mを計算する。また、各時刻ｔjの各方向ｘ、ｙ、ｚの空間差分ｆx,j，ｆy,j，ｆz,jから、各周波数ｍの観測時間周波数成分Ｆx,m，Ｆy,m，Ｆz,mを計算する。
【０２０４】
ステップＳ１２２では、変分パラメータ
2018142917.tif 000269
の初期値を設定する。
【０２０５】
ステップＳ１２４では、上記ステップＳ１２１によって計算された観測データＹと、初期値、又は前回更新された
2018142917.tif 000270
に基づいて、上記式（１３２）〜式（１４６）に従って、変分パラメータ
2018142917.tif 000271
を更新する。
【０２０６】
ステップＳ１２５において、予め定められた収束判定条件を満たしたか否かを判定し、収束判定条件を満たしていない場合には、上記ステップＳ１２４へ戻る。一方、収束判定条件を満たした場合には、ステップＳ１２６へ進む。
【０２０７】
ステップＳ１２６では、上記ステップＳ１２４で最終的に得られた
2018142917.tif 000272
に基づいて、各音源ｋの方向ベクトルn(k),音源距離R(k)を、各音源ｋの位置の推定結果として、出力部９０により出力して、音源定位処理ルーチンを終了する。
【０２０８】
＜実験＞
提案手法の性能を検証するため、残響環境下の複数音源定位の数値シミュレーションを行った。今回、x 方向、y 方向の2 次元モデルを使用した。部屋サイズは6 m×10 m×4 m とし、中心に0.03m間隔で7 つのマイクロホンを上記図2 のように配置した。壁面の反射係数は0:7308 及び0:4566（Sabineの残響公式による残響時間がそれぞれ0.5 s 及び0.2 s に相当）とした。マイクロホンのサンプリング周波数は32 kHz、短時間Fourier 変換のフレーム長は64 ms（オーバーラップは32 ms）、観測信号全体の長さは2779 ms 及び1665 ms とした。音源数は3 個、音源の位置は部屋の中心からそれぞれ(1, 0, 0) m、(-0.5, 0.87, 0) m、(-0.5,-0.87, 0) m とした（図５参照）。図５では、バツ印が音源位置を示し、マル印がマイクロホンアレイの中心位置を示している。音源信号は話速バリエーション型音声データベース（SRV-DB）を利用した。提案法として、本実施の形態の手法である変分推論に基づく手法（VBEM）、正しい音源数（3 個）を仮定したEM アルゴリズムによる手法（EM1）、誤った音源数（6 個）を仮定したEM アルゴリズムによる手法（EM2）を評価し、従来法であるMUSIC 法との比較を行った。EM アルゴリズムにおいては定常雑音を仮定し
2018142917.tif 000273
とした。どの手法においても、検出閾値をもとに音源の有無及び方向を推定し、そのうち真の音源方向との角度の誤差が±τ以内のものがあれば正解、なければ脱落誤り、検出された音源方向のうちどの真の音源方向にも属さないものを誤挿入として、F 尺度を算出した。各τについて、検出閾値を変化させたときの最も高いF 尺度のプロットを図６〜図９に示す。図６では、残響時間0.5 s、観測長2779 ms の条件での各手法の定位精度を示している。図７では、残響時間0.5 s、観測長1665 ms の条件での各手法の定位精度を示している。図８では、残響時間0.2 s、観測長2779 ms の条件での各手法の定位精度を示している。図９では、残響時間0.2 s、観測長1665 ms の条件での各手法の定位精度を示している。多くの場合、提案法である変分推論に基づく手法が他の手法より高精度な定位を達成していることが確認できた。
【０２０９】
以上説明したように、本実施の形態に係る音源定位装置によれば、基準のマイクロホンの各周波数の観測時間周波数成分、及び複数の方向の各々に対する各周波数の観測時間周波数成分からなる観測データＹが与えられたときの、複数の音源ｋの位置を表す変数Ｎ及び変数ρ、各周波数に対する各時刻において支配的となる音源のインデックスを示すインジケータを表す変数Ｚ、各音源ｋが支配的になる確率πkを定めるための変数Ｖ、複数の方向の各々に対する各周波数の観測時間周波数成分の分散を表す変数λ、並びに各音源ｋ及び雑音の各周波数の時間周波数成分のパワーを表す変数ζ及び変数γを含む未知変数Θの事後分布ｐ（Θ｜Ｙ）と変関数ｑ（Θ）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき目的関数を最小化するように、変数Ｎ、変数ρ、変数Ｚ、変数Ｖ、変数λ、変数ζ、変数γの各々の分布を表す各パラメータを推定し、各音源ｋの位置を推定することにより、音源数が未知の場合であっても、複数の音源を同時に定位することができる。
【０２１０】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０２１１】
例えば、上述の音源定位装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０２１２】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【０２１３】
１０   入力部
２０   演算部
２２   空間差分算出部
２４   時間周波数展開部
２５   音源位置推定部
２８   変数更新部
３０   収束判定部
９０   出力部
１００ 音源定位装置
【図１】
2018142917.tif 000274
【図２】
2018142917.tif 000275
【図３】
2018142917.tif 000276
【図４】
2018142917.tif 000277
【図５】
2018142917.tif 000278
【図６】
2018142917.tif 000279
【図７】
2018142917.tif 000280
【図８】
2018142917.tif 000281
【図９】
2018142917.tif 000282