JP6635903B2 - 音源位置推定装置、音源位置推定方法、及びプログラム - Google Patents

Description

本発明は、３次元空間に存在する音源の位置を推定する音源位置推定装置、音源位置推定方法、及びプログラムに関する。

従来の音源位置推定装置は、３次元空間に存在する音源の位置を推定するため、マイクロフォンアレイで取得したマルチチャネル信号に対して、時空間フーリエ変換を施して音場情報パラメータ群を算出する。そして、音場情報パラメータ群に対して部分空間に基づく手法（ＭＵＳＩＣ法：Multiple Signal Classification）を適用して音波の到来方向を推定する。最後に、異なる位置に設置された２つ以上のマイクロフォンアレイの推定結果に基づいて、音源の位置を推定する。上記に基づく音源位置推定装置は、例えば非特許文献１に開示されている。

D. Khaykin, et al., Coherent signals direction-of-arrival estimation using a spherical microphone array. in Proc. IEEE Workshop on Applications of Signal Processing to Audio and Acoustics(WASPAA), pp.221-224, 2009.

従来の方法は、残響音のない白色雑音の環境を仮定して音源位置の推定を行う。そのため、残響音（空間的に有色な雑音）が存在する実際の環境下では、音源位置の推定精度が低くなるという課題がある。

本発明は、この課題に鑑みてなされたものであり、残響音を含む雑音の構造を計算することで、音源位置の推定精度を向上させた音源位置推定装置、音源位置推定方法、及びプログラムを提供することを目的とする。

本実施形態の一態様に係る音源位置推定装置は、マイクロフォンアレイで取得したマルチチャネル信号から、音源の位置を推定する音源位置推定装置において、前記マルチチャネル信号を、周波数領域の信号に変換し、該信号に、ｎ次の球面調和関数の複素共役を乗じて音場の球面調和関数展開係数を計算する音場情報パラメータ群計算部と、前記球面調和関数展開係数から、前記音場の雑音の構造と複数の音波の到来方向を、ベイズの定理に基づいて計算する到来方向計算部とを備えることを要旨とする。

また、本実施形態の一態様に係る音源位置推定方法は、音源位置推定装置が行う音源位置推定方法であって、マイクロフォンアレイで取得したマルチチャネル信号を、周波数領域の信号に変換し、該信号に、ｎ次の球面調和関数の複素共役を乗じて音場の球面調和関数展開係数を計算する音場情報パラメータ群計算ステップと、前記球面調和関数展開係数から、前記音場の雑音の構造と複数の音波の到来方向を、ベイズの定理に基づいて計算する到来方向計算ステップとを行うことを要旨とする。

また、本実施形態の一態様に係るプログラムは、上記の音源位置推定装置としてコンピュータを機能させるためのプログラムである。

本発明によれば、音場の残響音を含む雑音の構造を計算することで、音源位置の推定精度を向上させた音源位置推定装置、音源位置推定方法、及びプログラムを提供することができる。

第１実施形態に係る音源位置推定装置の機能構成例を示す図である。 第１実施形態に係る音源位置推定装置の動作フローを示す図である。 球面座標系の例を示す図である。 第１実施形態に係る音源位置推定装置の音源位置計算部が、音源の位置を計算するのに用いる座標系を示す図である。 第２実施形態に係る音源位置推定装置の機能構成例を示す図である。 第２実施形態に係る音源位置推定装置が２組のマイクロフォンアレイで、音源の位置を計算する様子を概念的に示す図である。

以下、本発明の実施形態について図面を用いて説明する。複数の図面中同一のものに
は同じ参照符号を付し、説明は繰り返さない。

〔第１実施形態〕
図１に、第１実施形態に係る音源位置推定装置１の機能構成例を示す。図２に、その動作フローを示す。

音源位置推定装置１は、マイクロフォンアレイ１１，１２、音場情報パラメータ群計算部２０、到来方向計算部３０、及び音源位置計算部４０を備える。なお、ここでは２つのマイクロフォンアレイ１１，１２を備える例を示しているが、マイクロフォンアレイは１つでも構わない。

マイクロフォンアレイ１１は、例えば同一平面上にＱ個のマイクロフォンが等間隔に配列されて構成される。マイクロフォンアレイ１１は、音場の音を取得してマルチチャネル信号ｓ１１_１（ｔ），…，Ｓ１１_Ｑ（ｔ）を出力する。マイクロフォンアレイ１２も同様である。

音場情報パラメータ群計算部２０は、マルチチャネル信号ｓ１１_１（ｔ），…，Ｓ１１_Ｑ（ｔ）を、周波数領域の信号ｓ１１_１（ω），…，Ｓ１１_Ｑ（ω）に変換し（ステップＳ１）、該信号に球面調和関数の複素共役を乗じて音場の球面調和関数展開係数を計算する（ステップＳ２）。球面調和関数とは、軌道角運動を表現する周知の固有関数のことである。

音場情報パラメータ群計算部２０は、例えばマイクロフォンアレイ１１，１２が空中に配置されている場合、次式で球面調和関数展開係数ｐ_ｎｍ（ｋ）を計算する。なお、ｋは波数であり、２πｆ＝ｃｋ（ｃは音速）である。

ここで、ｊ_ｎはｎ次の球ベッセル関数、ｐ（ｒ，θ，φ，ｋ）は周波数領域のマルチチャネル信号であり、ｒ，θ，φは図３に示す球面座標系の各パラメータである。また、Ｑはマイクロフォンアレイを構成するマイクロフォンの数である。また、Ｙ_ｎｍ（・）^＊はｎ次の球面調和関数、ｎは次数、ｍは度数であり、＊は複素共役を表す。また、α_ｑは重み係数である。重み係数について詳しくは後述する。

また、マイクロフォンアレイ１１，１２が、剛体（例えば半径ａの球体）の表面上に配置されている場合は、次式で球面調和関数展開係数ｐ_ｎｍ（ｋ）を計算する。

ここで、ｈ_ｎはｎ次の球ハンケル関数である。

到来方向計算部３０は、音場情報パラメータ群計算部２０で計算した球面調和関数展開係数ｐ_ｎｍ（ｋ）から、音場の雑音の構造と複数の音波の到来方向を、周知のベイズの定理に基づいて計算する。

ベイズの定理で球面調和関数展開係数ｐ_ｎｍ（ｋ）を求めるに当たって、本実施形態では、次式で球面調和関数展開係数ｐ_ｎｍ（ｋ）を定義する。

ここでＹ^Ｈ（Ψ）は、信号の到来方向（Ψ）に依存するサイズ（Ｎ＋１）^２×Ｌの行列である。ｓは、信号を表す長さＬの縦ベクトルである。ｎ_ｎｍは、音場における雑音ベクトルｎを短時間フーリエ変換及び球面調和関数変換した値、ｋは波数である。

そして、音場の雑音の構造は、球面調和関数領域において次式の共分散行列で定義する。

ここで、Ｅは統計的期待値、Ｈは複素共役転置である。

時間周波数領域の共分散行列Ｋの非対角成分が0の場合は、時間周波数領域において雑音が空間的に白色（残響が存在しない状態）であり、非対角成分に値を持つ場合は、雑音が空間的に有色（室内残響が有り）であることを表す。したがって、共分散行列Ｋ_ｎｍで音場の雑音の構造を表すことができる。

到来方向計算部３０は、音場の球面調和関数展開係数ｐ_ｎｍ（ｋ）に対してＭＣＭＣ法（マルコフ連鎖モンテカルロ法）によるベイズ推定を行う。

（ベイズ推定）
到来方向計算部３０は、音場の雑音の構造と複数の音波の到来方向を、ベイズ推定で求めるに当たり、最初に到来方向Ψ^（１）と共分散行列Ｋ_ｎｍ ^（１）の初期値を設定する。

次に、音場の球面調和関数展開係数ｓ_ｊ ^{（ｐ＋１）}を、条件付き確率分布ｐ（ｓ_ｊ｜Ｐ_ｎｍ，Ψ^（ｐ），Ｋ_ｎｍ ^（ｐ））からサンプリングする。サンプリングは、例えば一般的な棄却法に基づいて行う。

ここで、添え字の（ｐ）は、ＭＣＭＣ法の繰り返しのインデックスである。繰り返し回数は、例えば1000回程度である。また、添え字のｊは、フレーム番号である。

音場の球面調和関数展開係数ｓ_ｋの事前分布は、例えば複素ガウス分布Ν（０，Φ_０）で表されると仮定する。その場合、音場の球面調和関数展開係数ｓ_ｊの条件付き確率分布は、次式で表せる。

ただし、

である。なお、Ｐ_ｎｍはベクトルの集合、Ｋ_ｎｍ及びΦは行列、Ψはベクトルである。

音場の球面調和関数展開係数ｓ_ｋは共役性を考慮して選択する。ただし、Φ_０＝ｃＩであり、ｃ（一定値）は例えばｃ＝１０^{−４〜−６}×Ｙ^ＨＫ_ｎｍ ^−１Ｙの最大固有値とする。ここで、Ｉは単位行列である。

到来方向計算部３０は、次に共分散行列Ｋ^{（ｐ＋１）}を、条件付き確率分布ｐ（Ｋ｜Ｐ_ｎｍ，Ｓ^{（ｐ＋１）}，Ψ^（ｐ））からサンプリングする。

共分散行列Ｋ_ｎｍの事前分布ｐ（Ｋ_ｎｍ）は、次式で表す複素逆ウィシャート分布（共役性を考慮して選択）で表される。

ただし、ν_０は仮想的な標本数である。共分散行列Ｋ_ｎｍの条件付き確率分布は、次の複素逆ウィシャート分布で表される。

ここで、ｐ（Ｐ_ｎｍ｜Ψ，Ｓ，Ｋ_ｎｍ）は、音場の球面調和関数の尤度である。尤度については後述する。

到来方向計算部３０は、次に到来方向Ψ^＊を、次式に示す提案分布Ｊ（Ψ^＊｜Ψ^（ｐ））からサンプリングする。

ここでσ_Ψ ^２は探索範囲を制御する定数である。提案標本Ψ^＊は、次に示すようにして採択／棄却する。

式（１２）は、ｕを一様分布ｕ（０，１）からサンプリングし、Ｒ＞ｕの時にΨ^＊を採択することで実現できる。ただし、採択率Ｒは次式で定義する。なお、ｐ（Ψ^＊）/ｐ（Ψ^（ｐ））は１であると仮定する。

なお、音場の球面調和関数展開係数ｐ_ｎｍ（ｋ）の尤度は次式に示すように記述できる。

ただし、雑音の球面調和関数展開係数ｎ_ｎｍ（ｋ）が複素ガウス分布Ν（０，Ｋ_ｎｍ）に従うと仮定している。

音場の球面調和関数展開係数の尤度Ｚは、次式のように記述できる。

ただし、ｔｒ（・）は行列のトレースを表し、Ｃは次式で表される。

到来方向計算部３０は、上記の音波の到来方向（Ψ）と雑音の構造（共分散行列Ｋ_ｎｍ）を同時に求める計算を繰り返して、音波の到来方向を出力する（ステップＳ３、図２）。

音源位置計算部４０は、複数の音波の到来方向から、音源の位置を計算する（ステップＳ４）。この例では、マイクロフォンアレイ１１のマルチチャネル信号から計算した音波の到来方向Ψ_１（θ_１，φ_１）と、マイクロフォンアレイ１２のマルチチャネル信号から計算した音波の到来方向Ψ_２（θ_２，φ_２）の２つの到来方向から音源の位置を計算する。

図４に、音波の到来方向を計算するのに用いる座標系を示す。本実施形態では、一方のマイクロフォンアレイ１１の中心Ｏ_１と、他方のマイクロフォンアレイ１２の中心Ｏ_２とを直線上に配置し、中心Ｏ_１を含む平面ｘ_１ｙ_１に対して中心Ｏ_２を含む平面ｘ_２ｙ_２が180°回転した座標系を用いる。

音源位置計算部４０は、２つのマイクロフォンアレイ１１，１２のそれぞれの中心Ｏ_１，Ｏ_２を原点として直線上に配置し、一方のマイクロフォンアレイ１１の第１の三次元座標（ｘ_１，ｙ_１，ｚ_１）と、該三次元座標に対してｘｙ平面を180°回転させた他方のマイクロフォンアレイ１２の第２の三次元座標（ｘ_２，ｙ_２，ｚ_２）の２つの座標系において、マイクロフォンアレイ１１と１２の中心間の距離をｃとし、第１の三次元座標で表した第１到来方向をΨ_１（θ_１，φ_１）とし、前記第２の三次元座標で表した第２到来方向をΨ_２（θ_２，φ_２）とした場合に、音源ｒ_ｓの位置を（ｘ_ｓ，ｙ_ｓ，ｚ_ｓ）として次式に基づいて計算する。

ここで、式（２）と式（４）に示した重み係数α_ｑの具体例を示す。

等角度サンプリングの場合の重み係数α_ｑは、次式で表せる。

ただし

また、ガウシアンサンプリングの場合は、

ただし、Ｐ_Ｎ＋１(cosθq）=0,0≦q≦Nを満たすN+1個のθを選択する。φは等角度サンプリングと同様に選択する。

また、一様サンプリングの場合は、

である。Ｑは、マイクロフォンアレイを構成するマイクロフォンの数である。

以上説明したように、音場の球面調和関数展開係数を算出し、音波の到来方向と雑音の共分散行列をＭＣＭＣ法によるベイズ推定によって推定することで、ある点で録音した音場ψ_ｓ（・）を、球面調和関数展開係数ｐ_ｎｍを用いて次式によって表現することが出来る。

本実施形態によれば、残響音を含む雑音の構造と音波の到来方向を同時に計算して音源位置の推定を行う。したがって、音源位置の推定精度を向上させた音源位置推定装置を提供することができる。

〔第２実施形態〕
図５に、第２実施形態に係る音源位置推定装置２の機能構成例を示す。音源位置推定装置２は、マイクロフォンアレイ１１，１２，１３、音場情報パラメータ群計算部２０２、到来方向計算部３０２、及び音源位置計算部４０２を備える。

音源位置推定装置２は、マイクロフォンアレイを３つ以上備える点で、音源位置推定装置１（図１）と異なる。マイクロフォンアレイは、４つ以上の偶数個を備えても良い。

音場情報パラメータ群計算部２０２は、３つのマイクロフォンアレイ１１，１２，１３で取得したマルチチャネル信号から、音場の球面調和関数展開係数ｐ_ｎｍ（ｋ）をそれぞれ計算する。音場の球面調和関数展開係数ｐ_ｎｍ（ｋ）の計算方法は、第１実施形態と同じである。

到来方向計算部３０２は、音場情報パラメータ群計算部２０２で計算したマイクロフォンアレイ１１，１２，１３に対応する音場の球面調和関数展開係数ｐ_ｎｍ（ｋ）から、音波の到来方向をそれぞれ計算する。マイクロフォンアレイ１１から求めた音場の球面調和関数展開係数ｐ_ｎｍ（ｋ）からは、第１到来方向Ψ_１（θ_１，φ_１）を計算する。マイクロフォンアレイ１２から求めた音場の球面調和関数展開係数ｐ_ｎｍ（ｋ）からは、第２到来方向Ψ_２（θ_２，φ_２）を計算する。マイクロフォンアレイ１３から求めた音場の球面調和関数展開係数ｐ_ｎｍ（ｋ）からは、第３到来方向Ψ_３（θ_３，φ_３）を計算する。到来方向Ψ_１，Ψ_２，Ψ_３の計算方法は、第１実施形態と同じである。

音源位置計算部４０２は、隣接するマイクロフォンアレイから取得した音波の到来方向から求めた複数の音源の位置を平均して出力する。２つの到来方向から音源の位置を計算する方法は、第１実施形態と同じである。

図６に、音源位置推定装置２が２以上の音波の到来方向から、１つの音源の位置を計算する方法を概念的に示す。Ｏ_１はマイクロフォンアレイ１１の中心であり、Ｏ_２はマイクロフォンアレイ１２の中心であり、Ｏ_３はマイクロフォンアレイ１３の中心である。

音源位置計算部４０２は、第１到来方向Ψ_１（θ_１，φ_１）と第２到来方向Ψ_２（θ_２，φ_２）とから音源の位置ｒ_１２を計算する。また、第２到来方向Ψ_２（θ_２，φ_２）と第３到来方向Ψ_３（θ_３，φ_３）とから音源の位置ｒ_２３を計算する。

そして、音源位置計算部４０２は、音源の位置ｒ_１２とｒ_２３とを平均した位置を、音源の位置として出力する。このように動作することで、音源位置推定装置２は、音源の位置に含まれる誤差を少なくすることができる。

以上説明したように本実施形態の音源位置推定装置１，２によれば、音波の到来方向と音場の雑音の構造を同時に計算することで、残響音が有る環境下においても音源の位置の推定精度を向上させることができる。

なお、マイクロフォンアレイの数を２つ又は３つの例で説明を行ったが、本発明はこの例に限定されない。例えば、１つのマイクロフォンアレイが１０個のマイクロフォンで構成される場合は、その内の５個を一方のマイクロフォンアレイ、また、他の５個を他方のマイクロフォンアレイと見立てれば良い。

また、３つ以上のマイクロフォンアレイを備えても良い。その場合は、隣接するマイクロフォンアレイから取得した音波の到来方向から求めた複数の音源の位置を平均することで、音源の位置の推定精度を更に向上させることが可能である。

このように本発明は、上記の実施形態に限定されるものではなく、その要旨の範囲内で変形が可能である。

なお、上記装置における処理部をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、各装置における処理部がコンピュータ上で実現される。

また、各処理部は、コンピュータ上で所定のプログラムを実行させることにより構成することにしてもよいし、これらの処理内容の少なくとも一部をハードウェア的に実現することとしても良い。

１、２：音源位置推定装置
１１、１２、１３：マイクロフォンアレイ
２０、２０２：音場情報パラメータ群計算部
２１：第１時間フーリエ変換部
２２：第１球面調和関数変換部
２３：第２時間フーリエ変換部
２４：第２球面調和関数変換部
３０、３０２：到来方向計算部
３１：第１ベイズ推定部
３２：第２ベイズ推定部
４０、４０２：音源位置計算部

Claims (8)

1. マイクロフォンアレイで取得したマルチチャネル信号から、音源の位置を推定する音源位置推定装置において、
   前記マルチチャネル信号を、周波数領域の信号に変換し、該信号に、ｎ次の球面調和関数の複素共役を乗じて音場の球面調和関数展開係数を計算する音場情報パラメータ群計算部と、
   前記球面調和関数展開係数から、前記音場の雑音の構造と複数の音波の到来方向を、ベイズの定理に基づいて計算する到来方向計算部と
   を備えることを特徴とする音源位置推定装置。
2. 請求項１に記載した音源位置推定装置において、
   前記球面調和関数展開係数をｐ_ｎｍ、信号の到来方向に依存するサイズ（Ｎ＋１）^２×Ｌの行列をＹ^Ｈ（Ψ）、Ｈは複素共役転置、前記信号を表す長さＬの縦ベクトルをｓ、前記音場の雑音ベクトルを短時間フーリエ変換及び球面調和関数変換した値をｎ_ｎｍ、及び（ｋ）を波数とした場合に、前記球面調和関数展開係数ｐ_ｎｍ（ｋ）は、次式で表されると仮定する
   

   

   ことを特徴とする音源位置推定装置。
3. 請求項１又は２に記載した音源位置推定装置において、
   前記音場の雑音の構造は、球面調和関数領域において次式の共分散行列で定義し、
   

   

   ここで、ｎ_ｎｍは前記音場における雑音ベクトルを短時間フーリエ変換及び球面調和関数変換した値、Ｈは複素共役転置、Ｅは統計的期待値であり、
   前記到来方向計算部は、前記共分散行列を、音場の球面調和関数展開係数の集合、信号の球面調和関数展開係数の集合、及び音波の到来方向の条件付き確率分布からサンプリングして更新する
   ことを特徴とする音源位置推定装置。
4. 請求項１乃至３の何れかに記載した音源位置推定装置において、
   複数の前記音波の到来方向から、音源の位置を計算する音源位置計算部を備え、
   前記音源位置計算部は、２つの前記マイクロフォンアレイのそれぞれの中心を原点として直線上に配置し、一方の前記マイクロフォンアレイの第１の三次元座標と、前記第１の三次元座標に対してｘｙ平面を180°回転させた他方のマイクロフォンアレイの第２の三次元座標の２つの座標系において、前記マイクロフォンアレイの中心間の距離をｃとし、前記第１の三次元座標で表した第１到来方向をΨ_１（θ_１，φ_１）とし、前記第２の三次元座標で表した第２到来方向をΨ_２（θ_２，φ_２）とした場合に、前記音源の位置を（ｘ_ｓ，ｙ_ｓ，ｚ_ｓ）として次式に基づいて計算する
   

   

   ことを特徴とする音源位置推定装置。
5. 請求項４に記載した音源位置推定装置において、
   ３つ以上のマイクロフォンアレイを備え、
   前記音源位置計算部は、隣接する前記マイクロフォンアレイから取得した音波の到来方向から２個以上の前記音源の位置を求め、２個以上の前記音源の位置を平均して出力する
   ことを特徴とする音源位置推定装置。
6. 音源位置推定装置が行う音源位置推定方法であって、
   マイクロフォンアレイで取得したマルチチャネル信号を、周波数領域の信号に変換し、該信号に、ｎ次の球面調和関数の複素共役を乗じて音場の球面調和関数展開係数を計算する音場情報パラメータ群計算ステップと、
   前記球面調和関数展開係数から、前記音場の雑音の構造と複数の音波の到来方向を、ベイズの定理に基づいて計算する到来方向計算ステップと
   を行うことを特徴とする音源位置推定方法。
7. 請求項６に記載した音源位置推定方法において、
   前記球面調和関数展開係数をｐ_ｎｍ、信号の到来方向に依存するサイズ（Ｎ＋１）^２×Ｌの行列をＹ^Ｈ（Ψ）、Ｈは複素共役転置、前記信号を表す長さＬの縦ベクトルをｓ、前記音場の雑音ベクトルを短時間フーリエ変換及び球面調和関数変換した値をｎ_ｎｍ、及び（ｋ）を波数とした場合に、前記球面調和関数展開係数ｐ_ｎｍ（ｋ）は、次式で表されると仮定する
   

   

   ことを特徴とする音源位置推定方法。
8. 請求項１乃至５の何れかに記載した音源位置推定装置としてコンピュータを機能させるためのプログラム。

Images