JP2017151221A

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

Info

Publication number: JP2017151221A
Application number: JP2016032366A
Authority: JP
Inventors: 弘和亀岡; Hirokazu Kameoka; 惇鈴木; Jun Suzuki
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: 2017-08-31
Anticipated expiration: 2036-02-23
Also published as: JP6488245B2

Abstract

【課題】雑音が存在する場合であっても、複数の音源を同時に定位することができるようにする。
【解決手段】空間差分算出部２２が、複数の方向に対し、観測信号の差分を算出し、時間周波数展開部２４が、基準のマイクロホンの観測信号を入力として、各周波数の観測時間周波数成分を出力すると共に、複数の方向に対して算出された観測信号の差分を入力として、複数の方向に対して、各周波数の観測時間周波数成分を出力する。音源位置推定部２５が、音源拘束偏微分方程式の周波数領域表現を用いて定められた、複数の音源と加法雑音が存在する場合における、複数の音源の位置を条件とした、基準のマイクロホンの各周波数の観測時間周波数成分、及び複数の方向に対する各周波数の観測時間周波数成分の確率密度値を大きくし、かつ、周波数領域表現における各時刻の各周波数成分において、高々１つの音源のみが支配的になるように、複数の音源の各々の位置を推定する。
【選択図】図３

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.

しかしながら、上記の波源拘束偏微分方程式に基づく手法は単一波源に対して成立する方程式をベースとしているため、複数の音源を同時に定位することはできない。また、雑音が存在する場合など、観測音響信号が偏微分方程式から逸脱する場合に脆弱であるという欠点を有している。

本発明は、上記事情を鑑みてなされたものであり、雑音が存在する場合であっても、複数の音源を同時に定位することができる音源定位装置、方法、及びプログラムを提供することを目的とする。

上記の目的を達成するために本発明に係る音源定位装置は、マイクロホンアレイにより入力された複数の音源からの音源信号が混合された観測信号から、前記複数の音源の各々の位置を推定する音源定位装置であって、複数の方向の各々に対し、前記マイクロホンアレイのうち、前記方向に並んだマイクロホンのペアにより入力された前記観測信号の差分を算出する空間差分算出部と、前記マイクロホンアレイのうち、基準のマイクロホンにより入力された前記観測信号を入力として、各周波数の観測時間周波数成分を出力すると共に、前記空間差分算出部によって前記複数の方向の各々に対して算出された前記観測信号の差分を入力として、前記複数の方向の各々に対して、各周波数の観測時間周波数成分を出力する時間周波数展開部と、前記時間周波数展開部により出力された、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分に基づいて、音源拘束偏微分方程式の周波数領域表現を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分の確率密度値を大きくし、かつ、周波数領域表現における各時刻の各周波数成分において、高々１つの音源のみが支配的になるように、前記複数の音源の各々の位置を推定する音源位置推定部と、を含んで構成されている。

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

本発明に係るプログラムは、上記の音源定位装置の各部としてコンピュータを機能させるためのプログラムである。

以上説明したように、本発明の音源定位装置、方法、及びプログラムによれば、音源拘束偏微分方程式の周波数領域表現を用いて定められた、前記複数の音源と加法雑音が存在する場合における、前記複数の音源の各々の位置を条件とした、前記基準のマイクロホンの各周波数の観測時間周波数成分、及び前記複数の方向の各々に対する各周波数の観測時間周波数成分の確率密度値を大きくし、かつ、周波数領域表現における各時刻の各周波数成分において、高々１つの音源のみが支配的になるように、前記複数の音源の各々の位置を推定することにより、雑音が存在する場合であっても、複数の音源を同時に定位することができる、という効果が得られる。

点音源から観測点ｒへ到来する球面波を示す図である。マイクロホンアレイの配置の一例を示す図である。本発明の実施の形態に係る音源定位装置の構成を示す概略図である。本発明の実施の形態に係る音源定位装置における音源定位処理ルーチンの内容を示すフローチャートである。フレーム幅がL = 16 の場合の実験結果を示す図である。フレーム幅がL = 32 の場合の実験結果を示す図である。フレーム幅がL = 64 の場合の実験結果を示す図である。

以下、図面を参照して本発明の実施の形態を詳細に説明する。本発明で提案する技術は、音響信号から波源位置を推定することを目的とした信号処理技術である。

＜本発明の実施の形態の概要＞
本発明の実施の形態は、上述した従来手法の利点を併せ持つ、小領域・瞬時観測による複数音源の波源定位を可能にする技術である。

本発明の実施の形態では、音源拘束偏微分方程式の周波数領域表現をベースにした音響信号の確率分布と、雑音を含む全音源のスパース性の仮定（複数の音源が混在する音響信号の時間周波数表現において、各時間周波数点で高々一つの音源のみが支配的であるという仮定）に基づき、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) の左辺を右辺に移項すると

が得られる。ここで、f_0,l, f_x,l, f_y,l, f_z,l を窓関数で窓掛けして取得された信号とする。切り出し区間の両端点の影響を無視できるものとすると、式(11) は周波数領域で

と表される。ただし、F_0,m, F_x,m, F_y,m, F_z,m はf_0,m, f_x,m, f_y,m, f_z,m の離散Fourier 変換であり、m は離散周波数インデックスである。

式(12) の右辺は雑音の存在や差分近似に伴う誤差により実際には必ずしも厳密に0 にはならない。

そこで、式(11) の右辺を

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

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

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

ここで、F_x,m, F_y,m, F_z,m, F_0,m を並べたベクトルとε_x,m, ε_y,m, ε_z,m, ε_0,m を並べたベクトルを

とし、f₀,..., f_L/2 を連結したベクトルとε₀,...,ε_L/2を連結したベクトルを

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

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

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

の複素正規分布

に従う。

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

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

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

により得られる。

＜複数音源の定位アルゴリズム＞
音声信号や楽音など実世界の音響信号の多くは時間周波数成分がスパースである。従って、複数の音源が存在する場合でも、各時間周波数点では高々一つの音源のみが支配的であると仮定できる場合が多い。この音源の時間周波数成分のスパース性の仮定と以上のｆの確率モデル化に基づき、音源が複数個存在する場合、および雑音が存在する場合の観測信号の確率分布を導くことができる。

信号の切り出しフレームの時刻のインデックスをｎ = 0,...,N-1、音源インデックスをk = 1,...,K + 1、音源k の音源位置パラメータをθ^(k) = ｛R^(k)，n^(k)｝、観測基準点における時間周波数成分のパワー

をσ_0,m,n ^(k)2 とする。また、k = K +1 は雑音に対応するものとする。ここで、雑音を含む全音源の時間周波数成分のスパース性を仮定し、周波数m，時刻n においてｚ_m,n番目の音源のみが支配的であるとき、それ以外の音源のパワーを0 とする。このとき、所与のｚ_m,n の下での観測信号の時間周波数成分とその空間差分

（以後単に観測信号と呼ぶ）の条件付き確率密度関数は

で与えられる。ただし、Σ_m,n ^(k)は

である。また、Γ_m,n は雑音の時間周波数成分の分散共分散行列で、周波数にのみ依存する正規化分散共分散行列モデルＷ_m と時刻にも依存する雑音のパワーν² _m,nの積

で表されるものとする。Ｗ_m の設定方法については後述する。

は、各音源ｋの音源位置θ^(k)と各音源ｋの音源信号の成分エネルギーσ^(k) _0,m ²とを含むすべての未知パラメータを表す。ｚ_m,n の事前確率をＰ(z_m,n= k) = α^(k) _m,n とすると、観測信号

の確率密度関数（

の尤度関数）は、以下の式となる。

以上より、複数の音源と雑音が存在する場合の各音源の最尤音源位置

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

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

y を不完全データと見なし、ｙと

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

は

で与えられるので、y が与えられた下での

のＺに関する条件付き期待値（Q 関数）は

で与えられる。この関数が増大するように

を更新するステップ（M ステップ）と、更新した
′
に代入し、

を計算するステップ（E ステップ）とを繰り返すことで

を局所的に最大にする

を求めることができる。以上より、本発明の実施の形態で説明する複数音源定位アルゴリズムは、以下の初期設定、Ｅステップ、Ｍステップからなる。

（初期ステップ）

を初期設定する。

（E ステップ）

に代入し、式(40) によりη^(k) _m,n を計算する。

（Ｍステップ）

が増大するように

を更新し、Ｅステップに戻る。

M ステップでは

を最大にするn^(k), R^(k),σ^(k) _0,m,n ²,Γ_m,n，α^(k) _m,nの同時最適解を解析的に求めることは難しいが、座標勾配法によりそれぞれの変数に関して

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

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

＜混合比α^(k) _m,nの更新式＞
ここまでは混合比α^(k) _m,n を時刻ｎと周波数ｍに依存する変数と見なしたが、すべてのn において等しい場合、すべてのm において等しい場合、すべてのm, n において等しい場合など、さまざまなバリエーションが考えられる。そこで、ここではα^(k) _m,n を以下の(a)〜(d)の場合の更新式を導く。

(a) そのまま変数として扱う場合
(b) 時刻n に依らない変数に制約する場合（α^(k) _m,n = π^(k) _m ）
(c) 周波数m に依らない変数に制約する場合（α^(k) _m,n =ρ^(k) _n ）
(d) 時刻n にも周波数m にも依らない変数に制約する場合（α^(k) _m,n = ｗ^(k)）

の形に制約する場合

の中でα^(k) _m,n に依存する項は

である。ただし、＝^ξはξに関係する項のみについての等号を表すものとする。いずれの場合でもα^(k) _m,n は

を満たす必要がある。従って、α^(k) _m,nの更新式は(a)〜(d) のケースでは

のようなラグランジアンを用いたLagrange 未定乗数法により得ることができる。それぞれのラグランジアンをα^(k) _m,n、π^(k) _m、ρ^(k) _n に関して偏微分し、0 と置くことにより、

を得る。

＜音源方向n^(k)の更新式＞
n^(k)は単位ベクトルなので、

の制約下で

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

のようなラグランジアンを用いてLagrange 未定乗数法で解くことができる。i 行j 列目の要素のみが1 で残りは0 であるような4×4 行列をE_i,jとすると、C_m(θ^(k)) は

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

を0 と置くことにより、

を得る。あとはn^(k) _x ²＋n^(k) _y ²＋n^(k) _z ²= 1 となるようにγ^(k)を二分法などで探索し、式(55)〜(57) に代入すれば良い。

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

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

のようにＳ^(k) に依存する項とそうでない項に分解することができるので、

のＳ^(k)に関する偏微分

を0 と置くことにより、

を得る。行列要素ごとの表記にすると

となる。ただし、

である。式(62) の分子における

はFast Fourier Transform (FFT) を用いて効率的に計算することができる。

＜パワーσ^(k) _0,m,n ² の更新式＞
上記と同様、Σ^(k)-1 _m,nを

のようにσ^(k) _0,m,n ² に依存する項とそうでない項に分解できるので、

のσ^(k) _0,m,n ² に関する偏微分を0 と置くことにより、

を得る。

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

のように、正規化分散共分散行列モデルW_m と周波数ｍ、時刻ｎにおけるエネルギーν² _m,nの積で表し、ν² _m,nを推定すべき変数とする。後述するがW_m は空間無相関モデルや拡散音場モデルなどから導かれる定数行列である。

のν² _m,n に関する偏微分

を0 と置くことにより、

を得る。

また、雑音エネルギーが時刻に依存しない場合を考える。この場合、雑音エネルギーはν² _mで表され、雑音の分散共分散行列は

となる。上記と同様、

のν² _mに関する偏微分

を0 と置くことにより、

を得る。

＜正規化分散共分散行列モデルW の設定方法＞
ここでは雑音の空間相関行列から正規化分散共分散行列モデルW_m の設定例を述べる。図2 のような7 本のマイクロホンの配置を想定する。ここで，f_i,0,...,f_i,L-1のFourier 変換をF_i,0,...,F_i,L-1とする。~f_m = (F_0,m,...,F_6,m)^Tおよびf_m = (F_x,m, F_y,m, F_z,m, F_0,m)^T の関係は

と書かれることから、~f_m= (F_0,m,...,F_6,m)^T の分散共分散行列をΨ_m とすると、f_mの分散共分散行列はBΨ_mB^T となる。従って、例えば空間的に無相関で等しいパワーの雑音を仮定する場合、Ψ_m は単位行列となるため、f_mの分散共分散行列W_m を

と置けば良い。

ある区域内で、エネルギー密度が一様でかつすべての方向に対するエネルギーの流れが等しい確率であるとみなせる分布をしている音場を拡散音場といい、残響環境の音場を良く近似的に表すことが知られている。拡散音場においては、２点間の空間相関係数が距離d にのみ依存し、

で与えられる。従って、拡散性雑音を仮定する場合、図2 のようなアレイ幾何の例では、~f_m =(F_0,m,...,F_6,m)^Tの分散共分散行列Ψ_m は

となる。これを用いて、f_mの分散共分散行列W_m をBΨ_mB^T と置けば良い。

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

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

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

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

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

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

時間周波数展開部２４は、空間差分算出部２２により得られた、基準点のマイクロホンにおける各時刻ｔ_lの観測信号f_0,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,m,n ²とを含む未知パラメータに基づいて、各音源ｋ、各時刻ｎ、及び各周波数ｍについて、上記式(40) によりη^(k) _m,nを計算する。

変数更新部２８は、期待値算出部２６によって算出されたη^(k) _m,nに基づいて、

ができるだけ大きくなるように、上記式（４７）〜式（５０）の何れかと、式（５５）〜式（５７）と、式（６２）と、式（６７）〜式（６８）と、式（７４）とに従って、各音源ｋの方向ベクトルn^(k),音源距離R^(k),パワーσ^(k) _0,m,n ²,混合比α^(k) _m,n,雑音共分散行列Γを更新する。

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

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

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

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

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

ステップＳ１２１では、上記ステップＳ１２０で得られた基準点のマイクロホンにおける各時刻ｔ_lの観測信号f_0,lから、各周波数ｍの観測時間周波数成分Ｆ_0,mを計算する。また、各時刻ｔ_lの各方向ｘ、ｙ、ｚの空間差分ｆ_x,l，ｆ_y,l，ｆ_z,lから、各周波数ｍの観測時間周波数成分Ｆ_x,m，Ｆ_y,m，Ｆ_z,mを計算する。

ステップＳ１２２では、各音源ｋの音源位置θ^(k)と各音源ｋの音源信号の成分エネルギーσ^(k) _0,m,n ²と各音源ｋの混合比α^(k) _m,nと雑音共分散行列Γとを含む未知パラメータに初期値を設定する。

そして、ステップＳ１２３では、上記ステップＳ１２１で取得した各周波数ｍの観測時間周波数成分Ｆ_x,m，Ｆ_y,m，Ｆ_z,m，Ｆ_0,mからなる観測周波数成分ｙと、上記ステップＳ１２２で初期設定された、又は後述するステップＳ１２４で前回更新された各音源ｋの音源位置θ^(k)、音源信号のパワーσ^(k) _0,m,n ²、混合比α^(k) _m,n、及び雑音共分散行列Γに基づいて、上記式(40) によりη^(k) _m,nを計算する。

ステップＳ１２４では、上記ステップＳ１２３で算出されたη^(k) _m,nに基づいて、

ができるだけ大きくなるように、上記式（４７）〜式（５０）の何れかと、式（５５）〜式（５７）と、式（６２）と、式（６７）〜式（６８）と、式（７４）とに従って、各音源ｋの方向ベクトルn^(k),音源距離R^(k),パワーσ^(k) _0,m,n ²,混合比α^(k) _m,n,雑音共分散行列Γを更新する。

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

ステップＳ１２６では、上記ステップＳ１２４で最終的に得られた各音源ｋの方向ベクトル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" はこの方法を意味する。）と比べ、高精度な定位が行えていることが分かった。

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

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

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

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

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

Claims

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

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

ただし、ｎが、時刻のインデックスであり、ｙ_ｍ,nが、時刻ｎの前記基準のマイクロホンの周波数ｍの観測時間周波数成分、及び時刻ｎの前記複数の方向の各々に対する周波数ｍの観測時間周波数成分を表し、ｚ_m,nは、時刻ｎの周波数ｍにおいて支配的な音源を表す変数であり、ｋは、音源のインデックスであり、Ｋ＋１は、雑音を表し、Γ_ｍは、前記加法雑音の周波数ｍの分散共分散行列であり、θ^(k)は、音源ｋの位置を表し、σ_x,m ^(k)2、σ_y,m ^(k)2、σ_z,m ^(k)2は、音源ｋからの音源信号の方向ｘ、ｙ、ｚの差分における周波数ｍの成分エネルギーを表し、σ_0,m,n ^(k)2は、音源ｋからの音源信号の前記基準のマイクロホンにおける時刻ｎの周波数ｍの成分エネルギーを表し、Ｒは、音源までの距離を表し、ｃは、音速を表し、Ｌは、周波数のインデッックスを規定するための定数であり、Ｔは、サンプリング周期を表し、ｎ_x、ｎ_y、ｎ_zは、音源へ向かう単位ベクトルの方向ｘ、ｙ、ｚの成分を表す。
前記音源位置推定部が推定することでは、ＥＭ（Expectation-Maximization）アルゴリズムにより、前記確率密度値が大きくなり、かつ、周波数領域表現における各時刻の各周波数成分において、高々１つの音源のみが支配的になるように、複数の音源ｋの各々までの距離Ｒ^(k)、前記複数の音源ｋの各々へ向かう単位ベクトルｎ^(k)、前記複数の音源ｋの各々の前記基準のマイクロホンにおける時刻ｎの周波数ｍの成分エネルギーσ_0,m,n ^(k)2、前記複数の音源ｋの各々についての時刻ｎの周波数ｍにおける混合比α_m,n ^(k)、及び前記加法雑音の周波数ｍの分散共分散行列Γ_ｍを繰り返し更新することにより、前記複数の音源の各々の位置を推定する請求項５記載の音源定位方法。
請求項１〜請求項３の何れか１項に記載の音源定位装置の各部としてコンピュータを機能させるためのプログラム。