JP5134525B2 - 方向情報分布推定装置、音源数推定装置、音源方向測定装置、音源分離装置、それらの方法、それらのプログラム - Google Patents

Description

この発明は、音響信号処理に用いる方向情報分布推定装置、音源数推定装置、音源方向測定装置、音源分離装置、それらの方法、それらのプログラムに関する。

従来から、音源の方向情報の分布を推定する技術がある。この技術は、例えば、音響信号処理の分野において、同時に複数の人が発話した音声が混ざった信号が観測されている時に、各人の方向を推定したり、各人の音声を分離抽出する際に重要である。

図１に従来の方向情報分布推定装置１００の機能構成例を示し、図２Ａに音の方向情報について得られたヒストグラムＨの一例を示し、当該方向情報は方向情報分布推定装置１００に入力されるものである。方向情報分布推定装置１００の原理については非特許文献１に記載されている。図２Ａ、後述する図２Ｂ、図２Ｃの横軸は、音源の（音の到来方向）の角度を示す。方向情報分布推定装置１００の目的は、ヒストグラムＨに、正規分布モデルをフィッティングさせる（近似させる）ことである。特にこの技術では、ヒストグラムＨ中にある複数の分布の山それぞれに意味がある場合を考える。例えば、それぞれの分布の山が音声信号源や電波信号源などの推定方向情報を表している場合などである。図２Ａの例では、分布の山は４つ存在し、４つの分布の山をそれぞれａ〜ｄとする。分布の山ａは約−１１５度、分布の山ｂは約−２０度、分布の山ｃは約６０度、分布の山ｄは約１５０度に位置している。

ヒストグラムＨから音源方向を推定する場合には、分布の山の角度の平均値を求めるのであるが、この平均値を求めるためには、各分布の山をそれぞれ１個の確率分布モデルでフィッティングさせることが要求される。

従来の技術では、ヒストグラム全体を例えば混合正規分布モデル（Ｇａｕｓｓｉａｎ ｍｉｘｔｕｒｅ ｍｏｄｅｌ：ＧＭＭ）でモデル化していた。混合正規分布Ｇは以下の式（１）により表される。

図１を用いて、従来の方向情報分布推定装置１００の各構成部の処理を簡単に説明する。まず、方向情報（図２Ａに示すヒストグラムＨ）が事後確率計算部２に入力され、事後確率計算部２でＭ個の正規分布モデルごとに事後確率を求める。そして、Ｍ個の事後確率はパラメタ更新部４に入力される。平均更新手段４２、分散更新手段４４、混合重み更新手段４６がそれぞれ、Ｍ個の事後確率等を用いて、パラメタ保持部８に保持させつつ、μ_ｍ、σ_ｍ、α_ｍをＥＭアルゴリズムにより更新する。更新過程の詳細は省略する。

そして、例えば、更新回数が閾値Ｔを超えると、収束後パラメタθ_ｖ（θ_ｖ＝（μ_ｍ、σ_ｍ、α_ｍ） ｍ＝１、．．．、Ｍ）を出力する。出力された収束後パラメタθ_ｖを用いて、音源方向の測定や音源の分離を行う。
M.Mandel,D.Ellis,and T.Jebara,"An EM algorithm for localizing multiple sound sources in reverberant environments,"Proc.Neural Info.Proc.Sys.,2006.

図２Ａに示すヒストグラムＨに８つの正規分布モデルを用いて、方向情報分布推定装置１００によるフィッティングさせた結果を図２Ｂに示す。図２Ｂの−１１５度付近の箇所（図２Ｂ記載のＰの箇所）に注目されたい。Ｐの箇所では、−１１５度付近の分布も１つの正規分布モデルでフィッティングさせたいにも関わらず、２つの正規分布モデルがフィッティングしてしまう。そして、図２Ｂ記載の正規分布モデルをそれぞれ合計することで、図２Ｃ記載の混合正規分布モデルが求められ、求められた混合正規分布モデルの収束後パラメタθ_ｖが方向情報分布推定装置１００から出力される。この場合であると、Ｐに２つの正規分布モデルがフィッティングしてしまっていることから、正確な収束後パラメタθ_ｖを求めることができないという問題がある。その結果、正確な音源分離や音源方向推定を行うことができない。

本発明では、音源からの音情報として、複数のピークを持つ方向情報が与えられた場合に、各ピークにそれぞれ１つの確率分布モデルをフィッティングさせることのできる方向情報分布推定装置を提供することである。

この発明は音源からの音情報が複数のピークを持つ場合に、Ｍ（Ｍは１以上の整数）個の確率分布モデルを用いて、各ピークにそれぞれ１つの確率分布モデルをフィッティングさせる方向情報分布推定装置である。当該方向情報分布推定装置は、パラメタ保持部と、事後確率計算部と、更新部と、を備える。パラメタ保持部は、現在の確率分布モデルの各パラメタを保持している。事後確率計算部は、音情報と、現在の確率分布モデルの各パラメタを用いて、Ｍ個の確率分布モデルごとに事後確率を計算する。更新部は音情報と、Ｍ個の確率分布モデルごとの事後確率を用いて、現在の確率分布モデルの各パラメタを更新し、各パラメタ値が収束していると判断した場合には更新された各パラメタを出力し、各パラメタ値が収束していないと判断した場合には、更新された各パラメタをパラメタ保持部に現在の確率分布モデルの各パラメタとして保持させる。そして、各パラメタのうち、混合重みの事前分布にディリクレ分布を用いる。

この発明の方向情報分布推定装置では、確率分布モデルのパラメタである混合重みに事前分布としてディリクレ分布を与えることで、各ピークに対して少数の確率分布モデルのフィッティングが可能であり、結果として各ピークにそれぞれ１つの確率分布モデルをフィッティングさせることができる。

以下に、発明を実施するための最良の形態を示す。なお、同じ機能を持つ構成部や同じ処理を行う過程には同じ番号を付し、重複説明を省略する。

図３に実施例１の方向情報分布推定装置２００の機能構成例を示し、図４に処理フローを示す。図５Ａに入力されるヒストグラムＨを示し、図５Ｂ、Ｃにそれぞれ、方向情報分布推定装置２００で得られる分布の一例を示す。また、図５ＡのヒストグラムＨは図２Ａと同様であるとする。

この実施例１では、用いるＭ個の確率分布モデルとしてＭ個の正規分布を用いる例を示し、入力される音情報を方向情報Ｄとし、方向情報Ｄの一例としてヒストグラムＨである場合を示す。音情報が方向情報である場合には、横軸は角度、縦軸は度数を表す。そして、方向情報分布推定装置２００はヒストグラムＨが複数のピークを持つ場合に、各ピークにそれぞれ１つの確率分布モデルをフィッティングさせる。正規分布モデルでなくとも、確率分布モデルであれば、どのモデルを用いても良い。ここで、各ピークは、１つの音源方向に対応するものである。
通常、方向情報Ｄには、２πのｋ（ｋは整数）倍の不定性が含まれるため、ここではそれを許すＷｒａｐｐｅｄ ＧＭＭを用いる。Ｗｒａｐｐｅｄ ＧＭＭであるＧは以下の式（２）で表すことができる。

θを混合正規分布の平均μ、分散σ、混合重みαをまとめたものを示し、つまり、θ＝（μ_ｍ、σ_ｍ、α_ｍ）＝（μ_１、σ_１、α_１、．．．、μ_ｍ、σ_ｍ、α_ｍ、．．．、μ_Ｍ、σ_Ｍ、α_Ｍ）となる。また、ｔを更新回数（時刻）とし、θに更新回数の概念を付与したもの、つまり、ｔ回更新したθをθ^ｔとすると、θ^ｔ＝（μ_１ ^ｔ、σ_１ ^ｔ、α_１ ^ｔ、．．．、μ_ｍ ^ｔ、σ_ｍ ^ｔ、α_ｍ ^ｔ、．．．、μ_Ｍ ^ｔ、σ_Ｍ ^ｔ、α_Ｍ ^ｔ）となる。また、記憶部１６には予め用いる正規分布モデルのモデル数Ｍと混合正規分布モデルの各パラメタの初期値θ^０が記憶されている。事前分布情報保持部１１０には、後述のハイパーパラメタφ、重みパラメタｃが保持されている。

方向情報分布推定装置２００には、Ｎ個の方向情報Ｄ＝｛ｄ_１、．．．、ｄ_ｎ、．．．、ｄ_Ｎ｝と重み係数Ａ＝｛ａ_１、．．．、ａ_ｎ、．．．、ａ_Ｎ｝が入力される。重み係数Ａは、方向情報の各要素ｄ_ｎ（ｎ＝１、．．．、Ｎ）に対する重み係数である。この重み係数は例えば、方向情報Ｄが得られる頻度や方向情報Ｄが得られた時の信頼度（取得信号のパワーや信号の瞬時的ＳＮ比など）により与えることができる。または全てのｎについてａ_ｎ＝１としても良い。

まず、ｔ＝０と設定し（つまり更新回数が０）、ｔ＝０のときの混合正規分布のパラメタθ^０の値を設定し、用いる正規分布のモデル数Ｍ、ｋの範囲であるＫ、更新回数閾値Ｔまたは差閾値Δ（後述する）を設定する。更新回数閾値Ｔまたは差閾値Δは、後述する収束判定処理の際に用いられる（ステップＳ２）。

事後確率計算部１２は、音情報（この実施例１では方向情報Ｄ）と、現在の混合正規分布のパラメタθ^ｔ（＝（μ_ｍ ^ｔ、σ_ｍ ^ｔ、α_ｍ ^ｔ ｍ＝１、．．．、Ｍ））から、Ｍ個の正規分布ごとに事後確率ｐ（ｍ、ｋ│ｄ_ｎ、θ^ｔ）を計算する（ステップＳ６）。またパラメタ保持部１８には、現在の混合正規分布のθ^ｔが保持されている。事後確率計算部１２は具体的には例えば、以下の式（５）により計算する。

式（５）の右辺の分子「ｐ（ｍ、ｋ、ｄ_ｎ│θ^ｔ）は上記式（３）（４）で表される「ｐ（ｍ、ｋ、ｄ_ｎ│θ）」内のθに更新回数ｔの概念を付与したものである。

次に、更新部１４は、方向情報Ｄと事後確率ｐ（ｍ、ｋ│ｄ_ｎ、θ^ｔ）を用いて、現在の混合正規分布の各パラメタθ^ｔを更新する（ステップＳ８）。以下、更新処理について詳細に説明する。更新部１４は更新処理の際に、ハイパーパラメタφ、重みパラメタｃを事前分布情報保持部１１０から取り出す。この発明では、パラメタθの更新処理は、正規分布のパラメタθの混合重みα_ｍに適切な事前分布を与え、例えばＥＭアルゴリズムにて行う。この実施例１では、混合重みα_ｍの事前分布として、ディリクレ分布を考える。ディリクレ分布の詳細は、参考文献１である「Ｃ．Ｍ．ビショップ著（元田、栗田他訳） 「パターン認識と機械学習（上）」、シュプリンガー・ジャパン２００７年 ｐ．７４−ｐ．７７」等に記載されている。ディリクレ分布は例えば以下の式（６）で表される。

ここで、αは混合重み行列であり、α＝｛α_１、．．．、α_ｍ、．．．、α_Ｍ｝で表され、Σ_ｍ ^Ｍα_ｍ＝１、０≦α_ｍ≦１という条件を満たす。これは混合正規分布のパラメタである混合重みの条件と同じであることに注意されたい。またβ（φ）は正規化項（ベータ分布）であり、ここで、ハイパーパラメタφを１より小さい正の値（例えば、０．９）に設定すると、α_ｍのごく少数のみが十分に大きな値を持ち、残りは０に近い値をとるようになる。この性質を式（１）で表される混合正規分布Ｇの混合重みα_ｍに対して適用することで、混合正規分布Ｇのうちの少数の正規分布のみに十分大きな混合重みがかかり、その他の正規分布の混合重みは０に近くなる。結果として、なるべく少数の正規分布によるフィッティングが可能である。

次に、この事前分布を含みながら、パラメタ更新を行うためのＥＭアルゴリズムを導出する。まず、最尤推定のためのコスト関数Ｌ（θ）は次のように与えられる。

また、重みパラメタｃは、式（９）の第１項と第２項の重みをコントロールするパラメタである。

となる。ここで、式（１１）のE［H］は式Hの期待値を示し、式（１２）中のｐ（ｍ、ｋ│ｄ_ｎ、θ^ｔ）は式（５）で表される事後確率分布である。ここで、従来法のＥＭアルゴリズムでは式（１２）中のｌｏｇ（ｐ（α））がないことに注意されたい。

また上述の通り、この場合には、音情報は音の到来方向を示すＮ個の方向情報ｄ_ｎ（ｎ＝１、．．．、Ｎ）であり、当該方向情報ｄ_ｎには２ｋπ倍（ｋは整数）の不定性が含まれているとし、混合正規分布は、ラップＧＭＭであり、ｃは重みパラメタであり、φはハイパーパラメタであり、Ｋはｋの範囲を示す。

図３中の更新部１４中の平均更新手段１４２が式（１３）より現在の平均μ_ｍ ^ｔを更新することで更新後の平均μ_ｍ ^ｔ＋１を出力する。分散更新手段１４４が式（１４）より分散σ_ｍを更新することで更新後の分散σ_ｍ ^ｔ＋１を出力する。混合重み更新手段１４６が式（１５）により混合重みα_ｍを更新することで更新後の混合重みα_ｍ ^ｔ＋１を出力する。パラメタ算出手段が、更新後の平均μ_ｍ ^ｔ＋１、分散σ_ｍ ^ｔ＋１、混合重みα_ｍ ^ｔ＋１についての更新後のパラメタθ^ｔ＋１を算出する。

各パラメタの更新処理が数回行われ（ステップＳ４）、更新部１４内の収束判定手段１５０は、更新されたθ^ｔ＋１に対して、予め定められた規則により、各パラメタ値が収束しているか否かの収束判定を行う（ステップＳ１０）。各パラメタ値が収束していると判断した場合には、更新されたパラメタθ^ｔ＋１を出力する。また、各パラメタ値が収束していないと判断した場合には、更新されたパラメタθ^ｔ＋１を現在の確率分布モデルの平均、分散、混合重みとしてパラメタ保持部１８に保持させる。そして、収束判定手段１５０が、各パラメタ値が収束していると判断するまで、ステップＳ４〜ステップＳ１０の処理を繰り返す。

ここで収束判定に用いる予め定められた規則の例を説明する。更新回数閾値Ｔを用いる例を説明すると、更新部１４内のカウント手段（図示せず）は更新回数ｔをカウントし、更新回数ｔが更新回数閾値Ｔを超えた場合には、十分更新しており、収束していると判断して、更新後のパラメタθを出力する。また、差閾値Δを用いる例を説明すると、以下の式（１６）の式を満たす場合には、収束していると判断して、更新後のパラメタθを出力する。
│Ｑ（θ│θ^ｔ＋１）−Ｑ（θ│θ^ｔ）│＜Δ （１６）
パラメタ算出手段１４８により算出されるθは図５Ｃの混合正規分布の各パラメタである。

また、この実施例１において、Ｋ＝０とすれば、ラップＧＭＭではなく、通常のＧＭＭによるフィッティングを行うことができる。この場合は、音情報Ｄとして方向情報である必要はない。例えば、音源からの音をＪ個のマイクロホン２０_ｊ（ｊ＝１、．．．、Ｊ）で収音した場合であると、マイクロホン２０_ｊとマイクロホン２０_ｊ’ （ｊ’＝１、．．．、Ｊであり、ｊ≠ｊ’）とのマイクロホン間位相差ｑ’_ｊｊ’を音情報Ｄとしても良い。またこの実施例１では混合重みα_ｍのみに事前分布を導入したが、各ガウス分布の平均μ_ｍと分散σ_ｍに対しても事前分布を導入することで、より精度の高いＧＭＭフィッティング（方向情報分布推定処理）を実現できる。また、各ガウス分布の各パラメタである平均μ_ｍ、分散σ_ｍ、混合重みα_ｍ、に事前分布を導入した場合にのＧＭＭフィッティングには、ＥＭアルゴリズムのほかにもベイズ推定などの様々なアルゴリズムが知られている。これらの拡張は当業者であれば、上記参考文献１などを参照すれば、容易に実現できるため、ここでは省略する。

この実施例１で説明したように、式（５）のハイパーパラメタφを１より小さな正の値（例えば、０．９）に設定すると、ディリクレ分布の性質からα_ｍのごく少数のみが十分に大きな値を持ち、残りは０に近い値をとるようになる。上記式（１）に示すＧＭＭの少数の正規分布のみに十分大きな混合重みα_ｍがかかり、そのほかの正規分布の重みは０に近くなる。この性質を用いることにより、なるべく少数の正規分布によるモデルフィッティングが可能になる。

実施例１では、Ｍ個の確率分布モデルとして、Ｍ個の正規分布モデルを用いたが、実施例２では、Ｍ個のフォン・ミーゼス（ｖｏｎ Ｍｉｓｅｓ）分布モデルを用いる。フォン・ミーゼス分布は角度の分布を表す関数であり、フォン・ミーゼス分布モデルの詳細は、参考文献２「Ｋ．Ｖ．Ｍａｒｄｉａ、”Ｓｔａｔｉｓｔｉｃｓ ｏｆ Ｄｉｒｅｃｔｉｏｎａｌ Ｄａｔａ”、Ａｃａｄｅｍｉｃ Ｐｒｅｓｓ、１９７２、３．４．９節」などに記載されている。フォン・ミーゼス分布を用いる効果は、正規分布モデルを用いた場合と比較すると、ｋおよびＫの値を考慮する必要がないため、演算処理が削減されることである。

この実施例２の方向情報分布推定装置３００の機能構成例、処理フローは図３、図４とほぼ同様であるが、図３中の分散更新手段１４４が拡散パラメタ更新手段１６０に代替されている点が異なる。以下、詳細に説明する。また、フォン・ミーゼス分布モデルのパラメタθをθ＝｛μ_ｍ、к_ｍ、α_ｍ｝とし、к_ｍは拡散パラメタである。

まず、事後確率計算部１２は音情報Ｄ（例えば方向情報Ｄ）とパラメタ保持部１８に保持されている現在のパラメタθ^ｔ＝｛α_ｍ ^ｔ、μ_ｍ ^ｔ、к_ｍ ^ｔ｝からＭ個それぞれのフォン・ミーゼス分布モデルに関する事後確率ｐ（ｍ│ｄ_ｎ、θ^ｔ）を求める。

この式（１７）は式（４）と対応しているものであり、式（１７）中の右辺の分子ｐ（ｍ、ｄ_ｎ│θ^ｔ）は、フォン・ミーゼス分布ｇ（ｄ_ｎ；μ_ｍ、к_ｍ）である。

ここで、−π＜ｄ_ｎ≦π、−π＜μ_ｍ≦πとし、また、к_ｍ＞０である。また、Ｉ_０（ｘ）は０次の第１種の変形されたベッセル関数である。

次に、更新部１４は、音情報Ｄと事後確率ｐ（ｍ│ｄ_ｎ、θ^ｔ）を用いて、フォン・ミーゼス分布のパラメタθ、つまり、平均μ_ｍ ^ｔ、拡散パラメタк_ｍ ^ｔ、混合重みα_ｍ ^ｔを更新する。以下、詳細に説明する。

平均更新手段１４２は平均μ_ｍ ^ｔを例えば以下の式（２０）により更新する。

ここで、ａｒｃｔａｎ（ｘ）は−π／２＜μ_ｍ＜π／２の値を返すのが一般的であるから、−π＜μ_ｍ＜πのデータを扱うには、以下の演算も行う。
式（２０）の値が負の場合、μ_ｍ ^ｔとμ_ｍ ^ｔ＋πの両方について、式（２１）に示すＱ関数の２次導関数を計算し、式（２１）の値が負になるほうをμ_ｍ ^ｔ＋１とする。

式（２０）の値が正の場合、μ_ｍ ^ｔ、μ_ｍ ^ｔ−πについて式（２１）を計算し、これが負になる方をμ_ｍ ^ｔとして保存する。
拡散パラメタ更新手段１６０は例えば以下の式（２２）により更新する。

ここで、Ｉ（к_ｍ ^ｔ＋１）を拡散パラメタ関数とする。к_ｍ ^ｔ＋１は解析的に得られないが、次のように得ることができる。拡散パラメタ関数Ｉ（к_ｍ ^ｔ＋１）は、単調増加関数である。そこで、ある範囲のк（例えば、０≦к≦１００）について、「к_ｍ ^ｔ＋１」と「Ｉ（к_ｍ ^ｔ＋１）」とを対応させたルックアップテーブルを用意しておく。当該ルックアップテーブルは、拡散パラメタ更新手段１６０中の記憶部（図示せず）に記憶させておけばよい。そして、Ｉ（к_ｍ ^ｔ＋１）が求まると、ルックアップテーブルを参照して、Ｉ（к_ｍ ^ｔ＋１）に対応するк_ｍ ^ｔ＋１を出力する。
混合重み更新手段１４６は、例えば、以下の式（２３）により混合重みα_ｍを更新する。

このようにして、更新部１４は分布パラメタθ^ｔ（＝｛α_ｍ、μ_ｍ、к_ｍ｝）を更新する。
この実施例２の方向情報分布推定装置３００のように、フォン・ミーゼス分布を用いることで、ｋに関する推定操作が不要であるため、実施例１の方向情報分布推定装置２００と比較して、計算コストやパラメタθの収束時間を削減できる。

［実験結果１］
図５を用いて、実施例１で説明した方向情報分布推定装置２００によるフィッティングの実験結果について説明する。実験条件として８（＝Ｍ）つの正規分布からなる混合正規分布をフィッティングさせ、ハイパーパラメタφを０．９とする。上述のように図５Ａに入力される方向情報ｄ_ｎについての図２Ａと同様のヒストグラムＨを示し、図５Ｂに方向情報分布推定装置２００のフィッティング処理による正規分布の結果を示し、図５Ｃに図５Ｂの正規分布を合計した混合正規分布（ＧＭＭ）を示す。図５Ｂ記載のＰの箇所（−１１５度付近）に注目すると、１つの正規分布でフィッティングできていることが理解されよう。従って、図５Ｃに示す求められる混合正規分布は、正確なものである。従って、実施例３〜５で説明する音源数推定処理、音源方向測定処理、音源分離処理も正確に行うことができる。

一方、上述のように、図１Ｂに示す従来の方向情報分布推定装置１００の実験結果については、Ｐの箇所では、２つの正規分布がフィッティングしてしまい、図１Ｃに示すＧＭＭは、不正確なものとなってしまう。

この実施例３では、実施例１、２で説明した方向情報分布推定装置２００、３００を用いた音源数測定装置４００について説明する。図６に音源数測定装置４００の機能構成例を示す。この実施例３の音源数測定装置４００は、Ｊ（Ｊは２以上の整数）個の収音手段２０_ｊ（例えば、マイクロホン ｊ＝１、．．．、Ｊ）に接続されている場合を説明する。そして、ある収録時間内（例えば５秒間など）に複数の音源から音が発せられた場合に、当該音をＪ個の収音手段２０_ｊで収録したとする（以下、状況Ｘという。）。この実施例３の音源数測定装置４００は、収録音のみを用いて、音を発した音源の数を推定する。

収音手段２０_ｊから入力された音信号をｘ_ｉ（ｓ）とし、ｓを離散時刻とする。周波数領域変換部３０は音信号ｘ_ｉ（ｓ）を周波数領域音信号Ｘ_ｊ（ｆ、τ）に変換する。ｆは周波数、τは時間フレーム番号である。また、この実施例３ではｎ＝τＦ＋ｆを考える。ただしＦは周波数領域の数である。
パワー推定部３２は周波数領域音信号から音のパワーを求める。求め方の一例として、パワー推定部３２は、各時間周波数（ｆ、τ）における周波数領域音信号Ｘ_ｊ（ｆ、τ）の信号パワー│Ｘ_ｊ（ｆ、τ）│^２を演算し出力する。出力された信号パワー│Ｘ_ｊ（ｆ、τ）│^２が、上述した重み係数ａ_ｎとして、以後用いられる。

また、到来方向推定部３４は周波数領域音信号から音の到来方向情報を求める。求め方の例を詳細に説明する。到来方向推定部３４は、収音手段間位相差演算手段３４２、到来方向情報生成手段３４４とで構成されている。まず、収音手段間位相差演算手段３４２が各フレームτ、各周波数ｆにおいて、各収音手段の全ての組み合わせ（マイクロホンペア）について収音手段間位相差ｑ’_ｊｊ’（ｆ、τ）を以下の式（２４）により求める。ただし、ｊ＝１、．．．、Ｊであり、ｊ’＝１、．．．、Ｊであり、ｊ≠ｊ’とする。
ｑ’_ｊｊ’（ｆ、τ）＝｛ａｒｇ［Ｘ_ｊ（ｆ、τ）Ｘ^＊ _ｊ’（ｆ、τ）］｝／２πｆ
（２４）

ただし、「^＊」は複素共役であることを示す。そして、全てのｑ’_ｊｊ’（ｆ、τ）を並べたベクトルをＱ’（ｆ、τ）とする。音の到来方向情報Ｑ（ｆ、τ）は音速Ｃと、各収音手段の座標系Ｄを用いて、以下の式（２５）により求められる。
Ｑ（ｆ、τ）＝ＣＤ^＋Ｑ’（ｆ、τ） （２５）

ここでＣは音速であり、「^＋」は、Ｍｏｏｒｅ−Ｐｅｎｒｏｓｅの擬似逆行列を表し、Ｄ＝［Ｄ_１−Ｄ_Ｌ、．．．、Ｄ_ｊ−Ｄ_Ｌ、．．．、Ｄ_Ｊ−Ｄ_Ｌ］^Ｔであり、Ｄ_ｊは収音手段２０_ｊの座標（ｘ、ｙ、ｚ）と並べたベクトルであり、ＬはＪ個の収音手段のうち代表として選ばれた代表収音手段のインデックスである。到来方向情報Ｑ（ｆ、τ）のｘｙｚ座標（ｘ_Ｑ、ｙ_Ｑ、ｚ_Ｑ）は、到来方向水平角（以下、単に「水平角」という。）をΨ（ｆ、τ）とし、到来方向仰角（以下、単に「仰角」という。）をΩ（ｆ、τ）とすると、以下の式（２６）で表すことができる。
Ｑ（ｆ、τ）＝（ｘ_Ｑ、ｙ_Ｑ、ｚ_Ｑ）
＝（ｃｏｓΨ（ｆ、τ）ｃｏｓΩ（ｆ、τ）、
ｓｉｎΨ（ｆ、τ）ｃｏｓΩ（ｆ、τ）、
ｓｉｎΩ（ｆ、τ）） （２６）

この実施例では、水平角Ψ（ｆ、τ）のみを用いる。求められた到来方向情報Ｑ（ｆ、τ）を方向情報ｄ_ｎとして用いる。また、方向情報ｄ_ｎについてヒストグラムを作成すると、図２Ａに示すヒストグラムＨが得られる。次に、方向情報ｄ_ｎ、方向情報ａ_ｎは方向情報分布推定装置２００（または３００）に入力され、実施例１（または実施例２）で説明した処理により、パラメタθが出力される。以降、出力されたパラメタθを決定後パラメタθとする。

音源数測定部３６は、決定後パラメタθの混合重みα_ｍ（ｍ＝１、．．．、Ｍ）のうち、混合重みが予め定められた第１閾値ε１（例えば１０^−６）よりも大きな値である方向情報分布モデルの個数Ｍ’を測定する。測定された個数Ｍ’を音源数として出力する。何故なら、方向情報分布推定装置２００（３００）の演算が十分収束している場合には、決定後パラメタθの中の混合重みα_ｍのうち十分大きな値を持つ個数はヒストグラム中の分布の山の数と等しくなるからである。以下の説明では、音源と認められたものについての方向情報分布モデルを音源該当方向情報分布モデルという。

また、方向情報分布推定装置２００（３００）の演算が十分に収束していない場合は、音源数測定部３６は、次のような推定処理を行うことが好ましい。まず、音源数測定装置４００内で、方向情報分布推定装置２００（実施例１で説明）を用いた場合には、音源数測定部３６は、混合重みα_ｍが第１閾値ε１よりも大きく、かつ分散σ_ｍが予め定められた第２閾値ε２（例えば１５度）よりも小さい方向情報分布モデルを音源該当方向情報分布モデルとして検出し、これら検出された音源該当方向情報分布モデルの個数Ｍ’を測定すればよい。また、音源数測定装置４００内で、方向情報分布推定装置３００（実施例２で説明）を用いた場合には、音源数測定部３６は、混合重みα_ｍが第１閾値ε１よりも大きく、かつ拡散パラメタк_ｍが第３閾値（例えば１０）よりも大きい方向情報分布モデルを音源該当方向情報分布モデルとして検出し、これら検出された音源該当方向情報分布モデルの個数Ｍ’を測定すればよい。

従来の方向情報分布推定装置１００は、ヒストグラムの各ピークに対して、正規分布をフィッティングさせると、図２Ｂに示すように、１つのピークに対して、２つの正規分布をフィッティングさせる場合がある。従って、方向情報分布推定装置１００を適用した音源数測定装置であると、誤った音源数測定をしてしまう。しかし、実施例１、２で説明した方向情報分布推定装置２００（または３００）により、図５Ｂに示すように、１つのピークに対して１つの確率分布モデル（例えば、正規分布モデルやフォン・ミーゼス分布）をフィッティングさせることができるので、正確な音源数を測定できる。

この実施例４では、音源方向測定装置５００について説明する。音源方向測定装置５００は、状況Ｘの場合に、収録音のみを用いて音源の方向を推定する。図７に音源方向測定装置５００の機能構成例を示す。図７の例では、音源方向測定装置５００は音源数測定装置４００（実施例３で説明）と音源方向測定部３８とで構成されている。

音源数測定装置４００の処理が終了すると、音源方向測定部３８は、音源該当方向情報分布モデルのインデックスｍ’｛ｍ’＝１、．．．、Ｍ’｝に対応する平均パラメタμ_ｍ’を方向情報分布推定装置２００から取り出し、推定すべき音源方向として当該平均パラメタμ_ｍ’を出力する。

この実施例４のように、音源方向測定装置５００内に具備する方向情報分布推定装置２００（３００）により、正確な方向情報分布処理がされることから、音源方向測定装置５００は正確な音源方向測定を行うことができる。

この実施例５では、音源分離装置６００について説明する。音源分離装置６００は、状況Ｘの場合に、収録音のみを用いて音源からの音信号を分離抽出する。図８に音源分離装置６００の機能構成例を示す。図８では、音源分離装置６００は音源数測定装置４００（実施例３で説明）と分離部４０、時間領域変換部４１とで構成されている。

音源数測定装置４００の処理が終了すると、分離部４０は、音源数測定装置４００で定められた音源該当方向情報分布モデルのインデックスｍ’について以下の処理を行う。

音源数測定装置４００が、方向情報分布推定装置２００を具備している場合には、分離部４０は、以下の式（２７）により、Ｍ’個の正規分布(式(５)参照)に関する事後確率ｐ（ｍ’、ｋ│ｄ_ｎ、θ^ｔ）を周辺化することで、周辺化事後確率ｐ（ｍ’│ｄ_ｎ、θ^ｔ）を求める。
ｐ（ｍ’│ｄ_ｎ、θ）＝Σ_ｋ＝−Ｋ ^Ｋｐ（ｍ’、ｋ│ｄ_ｎ、θ^ｔ） （２７）
また、音源数測定装置４００が、方向情報分布推定装置３００を具備している場合には、周辺化処理を行わず、上記式（１７）の演算結果を用いる。

また、周波数領域変換部３０からの周波数領域音信号Ｘ_ｊ（ｆ、τ）は分離部４０に入力される。分離部４０は周辺化事後確率と周波数領域音信号とを掛け合わせる。つまり、以下の式（２８）を演算することで、ｍ’番目の信号の推定に対応する周波数領域目的信号（分離信号）を出力する。
Ｙ_ｎｍ’＝Ｘ_ｎｐ（ｍ’│ｄ_ｎ、θ） （２８）
ここでＸ_ｎはＸ_１（ｆ、τ）を上述のｎ＝τＦ＋ｆで変形したものである。出力される周波数領域目的信号は時間周波数表現（ｆ、τ）を用いると、ｎ＝τＦ＋ｆより以下の式（２９）で表される事に留意されたい。
Ｙ_ｍ’（ｆ、τ）＝Ｘ_１（ｆ、τ）ｐ（ｍ’│Ψ（ｆ、τ）、θ） （２９）

そして、時間領域変換部４１は周波数領域目的信号Ｙ_ｍ’（ｆ、τ）を時間領域に変換することで、目的信号ｙ_ｍ’（ｔ）を求め、出力する。
この実施例５のように、音源分離装置６００内に具備する方向情報分布推定装置２００（３００）により、正確な方向情報分布処理がされることから、音源分離装置６００は正確な信号分離を行うことができる。

［実験結果２］
次に、方向情報分布推定装置２００を用いた、音源数測定装置４００（実施例３で説明）と音源分離装置６００（実施例５で説明）と（以下、「発明法」という。）、従来の方向情報分布推定装置１００を用いた音源数測定装置、音源分離装置（以下、「従来法」という。）とを比較した実験結果について説明する。まず図９を用いて、実験条件について説明する。長手方向４．４５ｍ（＝Ｌｂ）、短手方向３．５５ｍ（＝Ｌａ）、高さ２．５ｍの室内に、３つのマイクロホンＺ_１、Ｚ_２、Ｚ_３が、正三角形の各頂点に配置される。隣接するマイクロホン同士の間隔は４ｃｍであり、３つのマイクロホンの収音面がそれぞれ外側に向けられる。３つのマイクロホンＺ_１、Ｚ_２、Ｚ_３がなす正三角形の重心は、図９の左下の頂点Ｘから長手方向に２．５６ｍ（＝Ｌｄ）であり、短手方向に１．８ｍ（＝Ｌｃ）の箇所に位置する。また、３つのマイクロホンＺ_１、Ｚ_２、Ｚ_３を囲むように２〜４つのスピーカ（図９の例では４つのスピーカＳ_１、Ｓ_２、Ｓ_３、Ｓ_４とする。）が円周Ｒの方向に配置され、音を発しているとする。当該円周Ｒの半径は、５０ｃｍまたは１１０ｃｍであり、音の反響時間は１２８ｍｓであるとする。マイクロホンＺ_１、Ｚ_２、Ｚ_３、スピーカＳ_１、Ｓ_２、Ｓ_３、Ｓ_４の高さは全て１．２ｍとする。

実験項目については、（１）音源（スピーカ）が２個、３個、４個の場合について音源の数を測定できるか（音源数測定処理）、（２）音源からの音信号を分離できるか（音源分離処理）、である。これらの項目について、スピーカから発せられる音の音質を変えたり、スピーカがなす円周Ｒの半径を変えるなどして、２０通りの組み合わせについて実験を行った。

このような条件下で、図１０に実験結果を示す。図１０では、音源数測定処理については、２０通りのうち、どの程度の確率で音源数Ｗを判定しているか評価し、音源分離処理については信号対妨害音比（Ｓｉｇｎａｌ ｔｏ Ｉｎｔｅｒｆｅｒｅｎｃｅ Ｒａｔｉｏ：ＳＩＲ）を評価した。図１０からも理解されるように、音源数処理、音源分離処理については従来法では、誤った結果を出力しているが、発明法ではほぼ正確な結果を出していることが理解されよう。

＜ハードウェア構成＞
本発明は上述の実施の形態に限定されるものではない。また、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
また、上述の構成をコンピュータによって実現する場合、方向情報分布推定装置２００（３００）、音源数推定装置４００、音源方向測定装置５００、音源分離装置６００が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、処理機能がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。

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

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

また、本実施例で説明した方向情報分布推定装置２００（３００）、音源数推定装置４００、音源方向測定装置５００、音源分離装置６００は、ＣＰＵ（Central Processing Unit）、入力部、出力部、補助記憶装置、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）及びバスを有している（何れも図示せず）。

ＣＰＵは、読み込まれた各種プログラムに従って様々な演算処理を実行する。補助記憶装置は、例えば、ハードディスク、ＭＯ（Magneto-Optical disc）、半導体メモリ等であり、ＲＡＭは、ＳＲＡＭ(Static Random Access Memory)、ＤＲＡＭ (Dynamic Random Access Memory)等である。また、バスは、ＣＰＵ、入力部、出力部、補助記憶装置、ＲＡＭ及びＲＯＭを通信可能に接続している。

＜ハードウェアとソフトウェアとの協働＞
本発明の方向情報分布推定装置、音源数推定装置、音源方向測定装置、音源分離装置は、コンピュータの記録部に、本発明の各構成部としてを動作させるプログラムを読み込ませ、処理部、入力部、出力部などを動作させることで実現できる。また、コンピュータに読み込ませる方法としては、プログラムをコンピュータ読み取り可能な記録媒体に記録しておき、記録媒体からコンピュータに読み込ませる方法、サーバ等に記録されたプログラムを、電気通信回線等を通じてコンピュータに読み込ませる方法などがある。

従来の方向情報分布推定装置の機能構成例を示したブロック図。 Ａは従来の方向情報分布推定装置に入力されるヒストグラムであり、Ｂは従来の方向情報分布推定装置により正規分布をフィッティング処理の結果であり、Ｃはこれら正規分布についての混合正規分布を示す。 本実施例の方向情報分布推定装置を示した図。 本実施例の方向情報分布推定装置の処理フローを示した図。 Ａは本実施例の方向情報分布推定装置に入力されるヒストグラムであり、Ｂは本実施例の方向情報分布推定装置により正規分布をフィッティング処理の結果であり、Ｃはこれら正規分布についての混合正規分布を示す。 本実施例の音源数測定装置の機能構成例を示したブロック図。 本実施例の音源方向測定装置の機能構成例を示したブロック図。 本実施例の音源分離装置の機能構成例を示したブロック図。 実験２の実験条件を示した図。 実験２の結果を示す図。

Claims (13)

1. 音源からの音情報の分布が複数のピークを持つ場合に、Ｍ（Ｍは１以上の整数）個の確率分布モデルを用いて、当該確率分布の各パラメタを更新することで、各ピークにそれぞれ１つの確率分布モデルをフィッティングさせる方向情報分布推定装置であって、
   現在の確率分布モデルの各パラメタを保持しているパラメタ保持部と、
   前記音情報と、前記現在の確率分布モデルの各パラメタを用いて、Ｍ個の確率分布モデルごとに事後確率を計算する事後確率計算部と、
   前記音情報と、前記Ｍ個の確率分布モデルごとの事後確率を用いて、前記現在の確率分布モデルの各パラメタを前記更新し、各パラメタ値が収束していると判断した場合には更新された各パラメタを出力し、各パラメタ値が収束していないと判断した場合には、更新された各パラメタを前記パラメタ保持部に前記現在の確率分布モデルの各パラメタとして保持させる更新部と、を備え、
   前記更新部は、前記各パラメタのうち、混合重みの事前分布としてハイパーパラメタを１より小さい正の値に設定したディリクレ分布を用いることを特徴とする方向情報分布推定装置。
2. 請求項１記載の方向情報分布推定装置であって、
   前記確率分布モデルは、正規分布モデルであり、
   前記正規分布モデルの各パラメタは、混合重み、平均、分散、であることを特徴とする方向情報分布推定装置。
3. 請求項１記載の方向情報分布推定装置であって、
   前記確率分布モデルは、フォン・ミーゼス分布モデルであり、
   前記フォン・ミーゼス分布モデルの各パラメタは、混合重み、平均、拡散パラメタ、であることを特徴とする方向情報分布推定装置。
4. 複数の収音手段で入力された音信号を周波数領域に変換することで、周波数領域音信号を求める周波数領域変換部と、
   前記周波数領域音信号から音の到来方向情報を求める到来方向推定部と、
   前記周波数領域音信号のパワーを求めるパワー推定部と、
   前記音の到来方向情報を音情報とし、前記パワーを重み係数として、方向情報分布モデルを求める請求項１〜３何れかに記載の方向情報分布推定装置と、
   混合重みが予め定められた第１閾値よりも大きな値である音源該当方向情報分布モデルの個数Ｍ’を測定することで、音源数を求める音源数測定部と、を備える音源数推定装置。
5. 請求項４記載の音源数推定装置と、
   各音源該当方向情報分布モデルの各パラメタのうち、平均を音源方向として出力する音源方向測定部と、を備える音源方向測定装置。
6. 請求項４記載の音源数推定装置と、
   前記Ｍ’個の音源該当方向情報分布モデルごとの周辺化事後確率を求め、当該周辺化事後確率と前記周波数領域音信号とを掛け合わせることで、周波数領域目的信号を求める分離部と、
   前記周波数領域目的信号を時間領域に変換することで、目的信号を求める時間領域変換部と、を備える音源分離装置。
7. 音源からの音情報の分布が複数のピークを持つ場合に、Ｍ（Ｍは１以上の整数）個の確率分布モデルを用いて、各ピークにそれぞれ１つの確率分布モデルをフィッティングさせる方向情報分布推定方法であって、
   現在の確率分布モデルの各パラメタを保持しているパラメタ保持過程と、
   前記音情報と、前記現在の確率分布モデルの各パラメタを用いて、Ｍ個の確率分布モデルごとに事後確率を計算する事後確率計算過程と、
   前記音情報と、前記Ｍ個の確率分布モデルごとの事後確率を用いて、前記現在の確率分布モデルの各パラメタを更新し、各パラメタ値が収束していると判断した場合には更新された各パラメタを出力し、各パラメタ値が収束していないと判断した場合には、更新された各パラメタを前記パラメタ保持過程に前記現在の確率分布モデルの各パラメタとして保持させる更新過程と、を有し、
   前記更新過程は、前記各パラメタのうち、混合重みの事前分布としてハイパーパラメタを１より小さい正の値に設定したディリクレ分布を用いることを特徴とする方向情報分布推定方法。
8. 請求項７記載の方向情報分布推定方法であって、
   前記確率分布モデルは、正規分布モデルであり、
   前記正規分布モデルの各パラメタは、混合重み、平均、分散、であることを特徴とする方向情報分布推定方法。
9. 請求項７記載の方向情報分布推定方法であって、
   前記確率分布モデルは、フォン・ミーゼス分布モデルであり、
   前記フォン・ミーゼス分布モデルの各パラメタは、混合重み、平均、拡散パラメタ、であることを特徴とする方向情報分布推定方法。
10. 複数の収音手段で入力された音信号を周波数領域に変換することで、周波数領域音信号を求める周波数領域変換過程と、
    前記周波数領域音信号から音の到来方向情報を求める到来方向推定過程と、
    前記周波数領域音信号のパワーを求めるパワー推定過程と、
    前記音の到来方向情報を音情報とし、前記パワーを重み係数として、方向情報分布モデルを求める請求項７〜９何れかに記載の方向情報分布推定方法の各過程と、
    混合重みが予め定められた第１閾値よりも大きな値である音源該当方向情報分布モデルの個数Ｍ’を測定することで、音源数を求める音源数測定過程と、を有する音源数推定方法。
11. 請求項１０記載の音源数推定方法の各過程と、
    各音源該当方向情報分布モデルの各パラメタのうち、平均を音源方向として出力する音源方向測定過程と、を有する音源方向測定方法。
12. 請求項１０記載の音源数推定方法の各過程と、
    前記Ｍ’個の音源該当方向情報分布モデルごとの周辺化事後確率を求め、当該周辺化事後確率と前記周波数領域音信号とを掛け合わせることで、周波数領域目的信号を求める分離過程と、
    前記周波数領域目的信号を時間領域に変換することで、目的信号を求める時間領域変換過程と、を有する音源分離方法。
13. 請求項７〜９何れかに記載の方向情報分布推定方法、または請求項１０記載の音源数推定方法、または請求項１１記載の音源方向測定方法、または請求項１２記載の音源分離方法、の各過程をコンピュータに実行させるためのプログラム。

Images