JP2017151220A

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

Info

Publication number: JP2017151220A
Application number: JP2016032364A
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: JP6529451B2

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) より、

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

により得られる。

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

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

とする。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 ステップ）
θをθ´に代入し、式(39) により

を計算する。

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

＜Mステップ更新式＞
A(θ^(k)) はC₁(θ^(k)),...,C_L(θ^(k)) を対角に並べたブロック対角行列なので、

と書かれる。ただし、

である。また、Φ^(k) _m はΦ^(k)の4×4のブロック対角成分である。

M ステップではQ(θ,θ´) ができるだけ大きくなるようにn^(k), R^(k),σ^(k) _0,m ²,Γ を更新する。

Q(θ,θ´) を最大にするn^(k), R^(k),σ^(k) _0,m ²,Γの同時最適解を解析的に求めることは難しいが、座標勾配法によりそれぞれの変数に関してQ(θ,θ´) が最大となるように反復更新することでp(y｜θ) を局所最大化することができる（EM アルゴリズムでは、M ステップで補助関数が単調に増大することが保証されていれば収束性は保証される）。以下に、M ステップの更新方法を1 例示す。

＜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 と置くことにより、

を得る。ただし、[・]_i,j は行列のi 行j 列目の成分を表す。

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

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

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

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

を0 と置くことにより、

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

となる。ただし、

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

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

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

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

を得る。

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

のように、正規化分散共分散行列モデルW_m と周波数mの成分エネルギーν² _mの積で表し、ν² _mを推定すべき変数とする。後述するがW_m は空間無相関モデルや拡散音場モデルなどから導かれる定数行列である。Q(θ,θ´) のν² に関する偏微分を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 ²に基づいて、上記式(39) により

を計算する。

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

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

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

収束判定条件を満たしたときに、最終的に得られた各音源ｋの方向ベクトル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 ²とに初期値を設定する。

そして、ステップＳ１２３では、上記ステップＳ１２１で取得した各周波数ｍの観測時間周波数成分Ｆ_x,m，Ｆ_y,m，Ｆ_z,m，Ｆ_0,mからなる観測周波数成分ｙと、上記ステップＳ１２２で初期設定された、又は後述するステップＳ１２４で前回更新された各音源ｋの音源位置θ^(k)及び音源信号の成分エネルギーσ^(k) _0,m ²とに基づいて、上記式(39) により

を計算する。

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

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

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

ステップＳ１２６では、上記ステップＳ１２４で最終的に得られた各音源ｋの方向ベクトル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

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

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

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