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

Description

この発明は、音源位置を推定する技術に関する。

音響ビームフォーミングを用いてTV会議室等の空間における音源位置を特定する音源定位技術が知られている。音響ビームフォーミングは、マイクロホンアレイを空間的に異なる位置に配置し、それらの出力信号に時空間フィルタを施すことで指向性を制御し、特定方向からの音響信号のみを得る技術である。非特許文献１には、１個のマイクロホンアレイを用いて音響ビームフォーミングを行うことで特定の領域に存在する音源位置を特定する技術が知られている（例えば、非特許文献１参照。）。

WOLFGANG H. KUMMER, "Basic Array Theory", PROCEEDINGS OF THE IEEE, VOL.80, NO.1, JANUARY 1922

しかしながら、非特許文献１の技術では、音響ビームフォーミングの指向方向の共線上に複数の音源がある場合には、各音源を区別することができない。

この発明の目的は、音響ビームフォーミングの指向方向の共線上に複数の音源がある場合であっても、これらの複数の音源の位置を推定することができる音源位置推定装置、方法及びプログラムを提供することである。

この発明の１つの態様による音源位置推定装置によれば、Nを２以上の整数とし、Lを実際の音源の数とし、JをL≦J<Nの関係を満たす整数とし、Kを２以上の整数とし、S_iをi番目の仮想的な音源の位置とし、n番目の検査ポイントP_nを所定の位置とし、φ_k,n,iをマイクロホンアレイkから見た仮想的な音源の位置S_iを基準とする検査ポイントP_nの方向として、J個の検査ポイントP_n (n=1,2,…,J)の各々について、各マイクロホンアレイkで、ある仮想的な音源の位置S_iを基準として各検査ポイントP_nの方向をビームフォーミング測定することにより各周波数ωについての観測信号のパワースペクトルY_k(ω,φ_k,n,i)を取得し、これらの取得されたY_k(ω,φ_k,n,i)に基づく検査ポイントP_nについての観測信号のクロススペクトルを要素とするJ次元の観測結果ベクトル~yを生成する観測結果ベクトル生成部と、各マイクロホンアレイkについての伝達関数H_k(ω,φ_k,n,i)が記憶されている記憶部と、記憶部から読み込んだ伝達関数H_k(ω,φ_k,n,i)を用いて以下の式により定義される観測行列Bを生成する観測行列生成部と、

観測結果ベクトル~y及び観測行列Bを用いて~y=Bxという関係を満たす未知ベクトルxの中でノルムが最小となる未知ベクトルxを計算し、そのノルムが最小となる未知ベクトルxの要素の中で絶対値が大きいものからL個の要素に対応するL個の仮想的な音源の位置S_iを推定された音源の位置として出力する音源位置推定部と、を含む。

音響ビームフォーミングの指向方向の共線上に複数の音源がある場合であっても、これらの複数の音源の位置を推定することができる。

音源位置推定装置の例を説明するためのブロック図。 音源位置推定方法の例を説明するための流れ図。 発明の技術的背景を説明するための流れ図。 発明の技術的背景を説明するための流れ図。 発明の技術的背景を説明するための流れ図。

［技術的背景］
まず、この発明の技術的背景について説明する。

Kをマクロホンアレイの数、Lを実際の音源の数とする。また、Nをビームフォーミング測定を行うための当初の検査ポイントの数（従来法において必要な検査ポイントの数）とする。Kは２以上の正の整数であり、Lは１以上の正の整数であり、NはL<Nを満たす正整数である。すなわち、Nは２以上の正の整数である。

なお、JをL≦J<<Nという関係を満たす整数として、ビームフォーミング測定を行うための実際の検査ポイントの数をJ個にすることができるが、これについては後述する。

k番目（k=1,2,…,K）のマイクロホンアレイの伝達関数Ｈ_k（ω，θ）を以下のように表すものとする。Mはk番目のマイクロホンアレイを構成するマイクの数である。

ここで、a_k(ω,θ)は以下のように表されるものとする。

ここで、θは指向方向、ωは周波数、cは音速、jは虚数単位、dはk番目のマイクロホンアレイにおける各マイク間の距離、Λは行列又はベクトルの複素共役転置、Tは行列又はベクトルの転置を表す。なお、指向方向θは、図３に示すように、例えば各マイクロホンアレイを構成するマイクの配列方向に対する垂直方向とのなす角度で定義される。

また、dは、全てのマイクロホンアレイにおいて同じであってよいし、マイクロホンアレイごとに異なってもよい。このことを明らかにするために、図３では、dのことをd_kと表記している。同様に、Mは、全てのマイクロホンアレイにおいて同じであってよいし、マイクロホンアレイごとに異なってもよい。このことを明らかにするために、図３では、MのことをM_kと表記している。

今、位置Sに音源が存在していると仮定すると、各マイクロホンアレイで位置Sの方向をビームフォーミング測定することにより得られる観測信号のクロスパワースペクトル^Y(ω,Ｓ)は、以下のように表すことができる。なお、より正確には、^Y(ω,Ｓ)はクロスパワースペクトルの周波数帯域ωの成分である。

ここで、^Y(ω,Ｓ),^H(ω,Ｓ),^X(ω,Ｓ)は以下の式により表わされる。

θ_kは、k番目のマイクロホンアレイから見た位置Sの音源の方向（角度）であり、Y_k(ω,θ_k)は、k番目のマイクロホンアレイでθ_k方向をビームフォーミング測定して得られる観測信号のパワースペクトルである。また、X(ω,S)は位置Sの音源の信号のパワースペクトルである。

X(ω,S)は位置Sに音源が存在していれば大きな値を取り、位置Sに音源が存在していなければX(ω,S)=0となり得る。言い換えれば、位置Sの近くに音源が存在しているほど大きな値を取る。ゆえに、^Y(ω,S)は、位置Sの近くに音源が存在しているほどその値が大きくなる。

２次元平面上での音源の位置を推定するためには、例えば図４のように、対象とする２次元平面上の領域を格子で分割し、格子の交点（以下、格子点とする。）のそれぞれを検査ポイントP_n(n=1,2,…,N)として順番にビームフォーミング測定を行う。すなわち、Nは音源位置推定の対象となる空間全体を網羅する数が必要となる。図４の例では、左上の格子点からスタートし、（１）順次右側の格子点についてのビームフォーミング測定が行われ、（２）一番右側の格子点についてのビームフォーミング測定が終わったら一列下の一番左側の格子点についてのビームフォーミング測定が行われる。この（１）（２）を繰り返すことにより、最終的には右下の格子点についてビームフォーミング測定が行われる。その後、各検査ポイントについての観測信号のクロスパワースペクトルを求め、その値が最大となる位置に音源が存在すると推定すればよい。

仮に、位置S_i(i=1,2,…,N)にそれぞれ音源が存在すると仮定すると、観測信号のクロスパワースペクトル^Y(ω,S_i)は以下のように表せる。

ここで、^H(ω,P_n,S_i)は以下のように表すことができる。

ここで、θ_k(P_n)はk番目のマイクロホンアレイからみた検査ポイントP_nの方向（角度）であり、θ_k(S_i)はk番目のマイクロホンアレイからみた位置S_iの方向（角度）である。もちろん本当の音源の位置は分からないので、上述の式では、全ての検査ポイントP_nの位置のそれぞれに音源が存在すると仮定している。つまり、S_iとP_nは同数と仮定している。位置S_iに本当に音源が存在していれば、^Y(ω,S_i)は大きな値を取る。よって、観測信号により計算される^Y(ω,S_i)が大きな値を取る位置に音源が存在すると推定することができる。

ここで、計算量を少なくするために、圧縮センシングの原理を応用することにより、図５のように例えばランダムに選ばれた少数の検査ポイントについてビームフォーミング測定を行った結果のみに基づいて音源位置を推定してもよい。

圧縮センシングは、観測行列と観測結果ベクトルから未知ベクトルを推定する問題として定式化される。そこで、式(1)を行列形式で表すと以下のように書ける。
y=Hx …(3)

ここで、y,H,xは以下のように表される。

観測結果ベクトルyは、各検査ポイントについてビームフォーミング測定をして得られた観測信号により得られたクロススペクトルを要素とするベクトルであるとする。未知ベクトルxは、各位置S_iに音源が存在すると仮定したときの音源のパワースペクトルを要素とするベクトルとする。観測行列Hは、検査ポイントP_nの方向と音源が存在すると仮定される位置S_iの方向とに依存する伝達関数を要素とする行列であるとする。

ここで、実際の音源の数がLがNよりも小さな数（L<N）とすると、未知ベクトルxの要素中で実際には音源が存在しない位置S_iに関する^X(ω,S_i)の値は０となり、未知ベクトルxが疎なベクトルとなるので、圧縮センシングの理論を応用することで未知ベクトルxを推定することができる。

そこで、観測数である検査ポイントの数をNではなく、検査ポイントの数をJとする。ただし、L≦J<<Nとする。もちろん、L≦J<Nであってもよい。yからJ個の検査ポイントについての観測信号をランダムに抽出するための変換行列Aを考える。Aは、Hと無相関なJ×Nの行列である。例えば、J=5, N=10としたとき、

といったような行列を用いればよい。そして、
B=AH
とおくと、式(3)の両辺に左からAを乗算することで、
~y=Bx …(5)
が得られる。~yはJ個の検査ポイントについての観測信号から計算されるクロススペクトルを要素とするベクトル、すなわち、実際のビームフォーミング測定により得られる観測信号に基づくクロススペクトルの列である。よって、~yを観測結果ベクトル、Bを観測行列、xを未知ベクトルとして、圧縮センシングの理論を応用すれば、式(5)の疎な未知ベクトルxを推定することができる。すなわち、下記式を最小化する問題を解くことにより、未知ベクトルxを推定すればよい。

このような~xを求める方法としては、周知の基底追跡法(basis pursuit method)やmatching pursuit法及びその派生形、さらに不動点法，EMアルゴリズム（確率伝搬法）などを用いればよい（例えば、参考文献１参照。）。これによってL₁ノルムを最小化あるいは近似的に最小化する方法がある。

〔参考文献１〕Richard Baraniuk, “Compressing Sensing”, Lecture Notes in IEEE Signal Processing Magazine, Volume 24, July, 2007

~xが求まると、~xの要素の中で大きな値を取る方からL個の要素に対応するL個の位置Ｓ_ｉが音源位置として特定される。

［実施形態］
音源位置推定装置の構成例を図１に示す。また、この音源位置推定装置による音源位置推定方法の流れを図２に示す。

音源位置推定装置は、観測結果ベクトル生成部１、記憶部２、観測行列生成部３及び音源位置推定部４を例えば備えている。

Nを２以上の整数とし、i=1,2,…,Nとして、S_iをi番目の仮想的な音源の位置とする。当初の仮想的な音源の数はN個であるとする。

Lを実際の音源の数として、JをL≦J<Nの関係を満たす整数とし、n=1,2,…,Jとして、n番目の検査ポイントP_nを所定の位置とする。検査ポイントP_nの数はJ個である。

Kを２以上の整数とし、k=1,2,…,Kとする。マイクロホンアレイkの数はK個である。マイクロホンアレイkとは、k番目のマイクロホンアレイのことである。

φ_k,n,iを、マイクロホンアレイkから見た仮想的な音源の位置S_iを基準とする検査ポイントP_nの方向とする。φ_k,n,iは、例えば式(2)のように定義される。

観測結果ベクトル生成部１は、各マイクロホンアレイkで、ある仮想的な音源の位置S_iを基準として各検査ポイントP_nの方向をビームフォーミング測定することにより各周波数ωについての観測信号のパワースペクトルY_k(ω,φ_k,n,i)を取得し、これらの取得されたY_k(ω,φ_k,n,i)に基づいて以下の式により定義される観測結果ベクトル~yを生成する（ステップＳ１）。ここで、各検査ポイントにおいて測定対象とする仮想的な音源の位置は１つ(S_i)とする。つまり、当初のN個の仮想的な音源位置のうち、J個の仮想的な音源位置を対象としてビームフォーミング測定を行う。生成された観測結果ベクトル~yは、音源位置推定部４に提供される。

Y⁽ⁿ⁾(ω,S_i)は、観測結果ベクトル~yのn番目の要素を意味する。つまり、~yは、J個の検査ポイントP_n(n=1,2,…,J)の各々について、ある仮想的な音源の位置S_i(i∈{1,2,…,N})を基準として当該検査ポイントP_nの方向を各マイクロホンアレイでビームフォーミング測定することにより得た観測信号Y_k(ω,φ_k,n,i)のクロススペクトルY^(ｎ)(ω，S_i)を要素とするJ次元のベクトルである。

なお、ここでの検査ポイントP_nのインデックスnは、技術的背景の欄の記載で言うところのN個の検査ポイントのうちの先頭からＪ個を使うという意味ではなく、N個の検査ポイントの中からランダムに選んだＪ個の検査ポイント又はN個の検査ポイントの中の所定のＪ個の検査ポイントについて１から順にインデックスを振り直したものに相当する点に注意されたい。

また、仮想的な音源の位置S₁, …, S_Jは、S₁, …, S_J∈{S₁, …, S_N}であるが、{S₁,…,S_N}のうちの先頭からJ個の音源の位置を使うという意味ではなく、N個の仮想的な音源の位置の中からランダムに選んだＪ個の仮想的な音源の位置又はN個の仮想的な音源の位置の中の所定のＪ個の仮想的な音源の位置について１から順にインデックスを振り直したものに相当する。

検査ポイントP_nは、技術的背景の欄で述べたように、仮想的な音源の位置S_iと同じ位置にあってもよいし、異なる位置にあってもよい。また、技術的背景の欄で述べたように、格子点上にあってもよいし、格子点上になくてもよい。

記憶部２には、各マイクロホンアレイkについての伝達関数H_k(ω,φ_k,n,i)が予め記憶されているとする。

観測行列生成部３は、記憶部２から読み込んだ伝達関数H_k(ω,φ_k,n,i)を用いて以下の式により定義される観測行列Bを生成する（ステップＳ２）。生成された観測行列Bは、音源位置推定部４に提供される。

観測行列Bは、検査ポイントP_n( n=1,2,…, J)についての伝達関数を要素に持つ行列である。なお、観測行列Bは、N個の検査ポイントについての伝達関数を要素に持つ行列H（式(4)参照)の左から変換行列Aを乗じる(B=AH)ことにより生成しても良い。このとき変換行列AはJ×N行列であって、例えば、当初のN個の検査ポイントのうち実際にビームフォーミング測定に用いたJ個の検査ポイントに対応する行以外の行の要素は全て0であり、J個の検査ポイントに対応するJ個の行がJ個の一次独立な単位ベクトルとなるように構成すればよい。単位ベクトル中の1となる要素の位置は、観測結果ベクトル~yにおける検査ポイントの順序に依存して決まる。

音源位置推定部４は、観測結果ベクトル~y及び上記観測行列Bを用いて~y=Bxという関係を満たす未知ベクトルxの中でノルムが最小となる未知ベクトルxを計算し、そのノルムが最小となる未知ベクトルxの要素の中で絶対値が大きい要素に対応する仮想的な音源の位置S_iを推定された音源の位置として出力する（ステップＳ３）。

言い換えれば、音源位置推定部４は、観測結果ベクトル~y及び上記観測行列Bを用いて、例えば式(6)の最小化問題を解くことにより、未知ベクトルxを求める。

~y=Bxという関係を満たす未知ベクトルxの中でノルムが最小となる未知ベクトルxの計算は、技術的背景の欄で記載したように、例えば参考文献１に記載された手法や、周知の基底追跡法(basis pursuit method)やmatching pursuit法及びその派生形、さらに不動点法，EMアルゴリズム（確率伝搬法）などを用いればよい。

Cを1≦C≦Lの関係を満たす正の整数として、未知ベクトルxの要素の中で絶対値が大きい要素とは、未知ベクトルxのJ個の要素を大きい順に並べたときの上位C個の要素のことである。

このように、複数のマイクロホンアレイを用いることにより、何れかのマイクロホンアレイの音響ビームフォーミングの指向方向の共線上に複数の音源がある場合であっても、これらの複数の音源の位置を推定することができる。

また、圧縮センシングの原理を応用し、測定する検査ポイントP_nの数を、当初の仮想的な音源の位置S_iの数であるNよりも少なくした場合には、少ない演算量で音源位置を推定することができる。

［変形例等］
式（6）は疎な未知ベクトルxを求めるため一手段である。式（6）は下記のようなL₀擬ノルム最小化問題の緩和問題であるので、直接的にL₀擬ノルム最小化を行ってもよい。これは組み合わせ最適化問題を解くことによって実行可能である。

また、近似的に最小化を行う方法が提案されている。さらに、L_p（擬）ノルム（0<p<1）の最小化問題に書き換えることも可能であり、これを最小化又は近似的に最小化することによっても疎なxが求められることが知られている。

したがって、音源位置推定部４で用いられる未知ベクトルxのノルムは、L₁ノルムであってもよいし、L₀擬ノルムであってもよいし、L_p（擬）ノルム（0<p<1）であってもよい。

上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。

上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体の例は、非一時的な（non-transitory）記録媒体である。このような記録媒体の例は、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等である。

このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。処理の実行時、このコンピュータは、自己の記録装置に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。

上記実施形態では、コンピュータ上で所定のプログラムを実行させて本装置の処理機能が実現されたが、これらの処理機能の少なくとも一部がハードウェアで実現されてもよい。

１ 観測結果ベクトル生成部
２ 記憶部
３ 観測行列生成部
４ 音源位置推定部

Claims (2)

1. Nを２以上の整数とし、Lを実際の音源の数とし、JをL≦J<Nの関係を満たす整数とし、Kを２以上の整数とし、S_iをi番目の仮想的な音源の位置とし、n番目の検査ポイントP_nを所定の位置とし、φ_k,n,iをマイクロホンアレイkから見た仮想的な音源の位置S_iを基準とする検査ポイントP_nの方向として、
   J個の検査ポイントP_n(n=1,2,…,J)の各々について、各マイクロホンアレイkで、ある仮想的な音源の位置S_iを基準として各検査ポイントP_nの方向をビームフォーミング測定することにより各周波数ωについての観測信号のパワースペクトルY_k(ω,φ_k,n,i)を取得し、これらの取得されたY_k(ω,φ_k,n,i)に基づく検査ポイントP_nについての観測信号のクロススペクトルを要素とするJ次元の観測結果ベクトル~yを生成する観測結果ベクトル生成部と、
   各上記マイクロホンアレイkについての伝達関数H_k(ω,φ_k,n,i)が記憶されている記憶部と、
   上記記憶部から読み込んだ伝達関数H_k(ω,φ_k,n,i)を用いて以下の式により定義される観測行列Bを生成する観測行列生成部と、
   

   

   
   上記観測結果ベクトル~y及び上記観測行列Bを用いて~y=Bxという関係を満たす未知ベクトルxの中でノルムが最小となる未知ベクトルxを計算し、そのノルムが最小となる未知ベクトルxの要素の中で絶対値が大きいものからL個の要素に対応するL個の仮想的な音源の位置S_iを推定された音源の位置として出力する音源位置推定部と、
   を含む音源位置推定装置。
2. Nを２以上の整数とし、Lを実際の音源の数とし、JをL≦J<Nの関係を満たす整数とし、Kを２以上の整数とし、S_iをi番目の仮想的な音源の位置とし、n番目の検査ポイントP_nを所定の位置とし、φ_k,n,iをマイクロホンアレイkから見た仮想的な音源の位置S_iを基準とする検査ポイントP_nの方向とし、
   記憶部には、各上記マイクロホンアレイkについての伝達関数H_k(ω,φ_k,n,i)が記憶されているとして、
   観測ベクトル生成部が、J個の検査ポイントP_n(n=1,2,…,J)の各々について、各マイクロホンアレイkで、ある各仮想的な音源の位置S_iを基準として各検査ポイントP_nの方向をビームフォーミング測定することにより各周波数ωについての観測信号のパワースペクトルY_k(ω,φ_k,n,i)を取得し、これらの取得されたY_k(ω,φ_k,n,i)に基づく検査ポイントP_nについての観測信号のクロススペクトルを要素とするJ次元の観測結果ベクトル~yを生成する観測結果ベクトル生成ステップと、
   観測行列生成部が、上記記憶部から読み込んだ伝達関数H_k(ω,φ_k,n,i)を用いて以下の式により定義される観測行列Bを生成する観測行列生成ステップと、
   

   

   
   音源位置推定部が、上記観測結果ベクトル~y及び上記観測行列Bを用いて~y=Bxという関係を満たす未知ベクトルxの中でノルムが最小となる未知ベクトルxを計算し、そのノルムが最小となる未知ベクトルxの要素の中で絶対値が大きいものからL個の要素に対応するL個の仮想的な音源の位置S_iを推定された音源の位置として出力する音源位置推定ステップと、
   を含む音源位置推定装置。

Images