JP6616472B2 - クラスタリング装置、クラスタリング方法及びクラスタリングプログラム - Google Patents

Description

本発明は、クラスタリング装置、クラスタリング方法及びクラスタリングプログラムに関する。

複数のマイクロホンにより取得した目的信号及び雑音を含む複数の観測信号から雑音を抑圧して目的信号を取得する技術として雑音除去技術がある。従来の雑音除去技術では、時間周波数領域（例えば短時間フーリエ変換領域）における目的信号のスパース性の仮定の下、「目的信号及び雑音からなる観測信号は、時間周波数領域の疎な点においてのみ目的信号を含み、他の点においては目的信号を含まず雑音のみを含む」というモデルを用いる。従来の雑音除去技術では、このモデルに基づき、観測信号の各時間周波数成分を、目的信号を含むものと含まないものとに分類する。

すなわち、従来の雑音除去技術では、上記モデルに基づき、「観測信号の各時間周波数成分から抽出した信号到来方向を特徴づける到来方向特徴量の分布は、目的信号を含む時間周波数成分においては目的信号の到来方向に対応する分布を持ち、目的信号を含まない時間周波数成分においては雑音の到来方向に対応する分布を持つ」と仮定する。そして、従来の雑音除去技術では、この仮定に基づき、到来方向特徴量に対してクラスタ数２のクラスタリングを適用することで、観測信号の時間周波数成分を、目的信号を含む時間周波数成分からなるクラスタ（以下「目的信号と対応するクラスタ」）と目的信号を含まない時間周波数成分からなるクラスタ（以下「目的信号と対応しないクラスタ」）とにクラスタリングする。

特開２０１３−０５４２５８号公報

しかしながら、上記技術では、目的信号と対応するクラスタと目的信号と対応しないクラスタへのクラスタリングはできても、得られた２つのクラスタのうちのどちらが目的信号と対応するクラスタであるかが不明である。そのため、各クラスタが、目的信号と対応するクラスタであるか、目的信号と対応しないクラスタであるかの判定を人為的に行う必要があるという問題がある。

本願が開示する実施形態の一例は、複数のマイクロホンにより取得した観測信号の時間周波数成分をクラスタリングすることで得られたクラスタの各々が、目的信号と対応するクラスタであるか、目的信号と対応しないクラスタであるかを判定することを目的とする。

上述した課題を解決し、目的を達成するために、本発明の本願のクラスタリング装置は、複数のマイクロホンで目的信号を含む音を観測した複数の観測信号の時間周波数成分を前記目的信号と対応するクラスタと前記目的信号と対応しないクラスタとにクラスタリングするクラスタリング装置であって、各観測信号の時間周波数成分からなる観測信号ベクトルに対応する特徴量ベクトルを抽出する特徴抽出部と、前記特徴量ベクトルに所定の混合モデルを当てはめることで前記クラスタ毎かつ時間周波数成分毎の事後確率を算出するクラスタリング部と、を含み、前記混合モデルは、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和で表され、前記各クラスタに関する特徴量ベクトルのモデル分布には、前記目的信号と対応しないクラスタに関する前記特徴量ベクトルのモデル分布は、前記目的信号と対応するクラスタに関する前記特徴量ベクトルのモデル分布よりも小さい到来方向局在度を持つという制約が課されていることを特徴とする。

本願が開示する実施形態の一例によれば、例えば、複数のマイクロホンにより取得した観測信号の時間周波数成分をクラスタリングすることで得られたクラスタの各々が、目的信号と対応するクラスタであるか、目的信号と対応しないクラスタであるかを判定することができる。

図１は、実施形態１（実施形態２）に係るクラスタ判定装置の一例を示す図である。 図２は、実施形態１（実施形態２）に係る第１例の到来方向局在度算出部を示す図である。 図３は、実施形態１（実施形態２）に係る第２例の到来方向局在度算出部を示す図である。 図４は、実施形態１（実施形態２）に係る第３例の到来方向局在度算出部を示す図である。 図５は、実施形態１に係る判定部の一例を示す図である。 図６は、実施形態１に係るクラスタ判定処理の一例を示すフローチャートである。 図７は、実施形態２に係るクラスタ判定処理の一例を示すフローチャートである。 図８は、実施形態３に係るクラスタリング装置の一例を示す図である。 図９は、実施形態３に係るクラスタリング処理の一例を示すフローチャートである。 図１０は、プログラムが実行されることにより、実施形態に係るクラスタ判定装置及びクラスタリング装置が実現されるコンピュータの一例を示す図である。

以下、本願の開示技術に係る実施形態の一例の説明に先立ち、開示技術の前提及び背景となる従来技術について説明する。従来技術の説明後、本願の開示技術に係る実施形態の一例について説明する。以下、重複する構成及び処理については、後出部分について同一符号を付与して説明を省略する。また、以下、下記（１−１）式〜（１−３）式に示すように、「“＾Ａ”を“Ａの直上に＾を付した記号”」、「“〜Ａ”を“Ａの直上に〜を付した記号”」、「“Ａ_C ^B”を“Ａに上付きＢを付し、下付きＣを付した記号”」と同等であるとする。また、以下、記号Ａがスカラーを示す場合には単に“Ａ”と表記し、記号Ａがベクトルを示す場合には“ベクトルＡ”と表記し、記号Ａが行列を示す場合には“行列Ａ”と表記し、記号Ａが集合を示す場合には“集合Ａ”と表記する。

［従来の雑音除去技術］
複数のマイクロホンにより取得した目的信号及び雑音を含む観測信号を用いて雑音を抑圧した強調信号を取得する従来の雑音除去技術について説明する。

従来の雑音除去技術では、「目的信号は時間周波数平面上の疎な点においてのみパワーを持つ」（スパース性）と仮定する。この仮定は音声信号に対して特に精度よく成立することが知られている。スパース性を仮定すれば、観測信号の各時間周波数成分は、下記（２）式に示す２つの状態のいずれかを取る。

状態Ｓ⁽⁰⁾では観測信号は雑音のみを含むことに着目すれば、雑音共分散行列＾Φ_f ^vは、下記（３）式により推定できる。

ここで、ベクトルｙ_tf:＝[ｙ_tf ⁽¹⁾・・・ｙ_tf ^(M)]^T∈Ｃ^M（Ｃ^MはＭ次元複素ベクトル空間）は、Ｍ個のマイクロホンにおける観測信号の短時間フーリエ変換（以下、観測信号ベクトルと呼ぶ）、γ_tf ⁽ⁿ⁾は状態Ｓ⁽ⁿ⁾（ｎ＝０，１）の検出により得られる各時間周波数点での状態Ｓ⁽ⁿ⁾の事後確率、ｔ∈｛１，・・・，Ｔ｝は時間フレームの番号、ｆ∈｛１，・・・，Ｆ｝は周波数ビンの番号、Ｔは時間フレームの数、Ｆは周波数ビンの数、上付きのＴは転置、上付きのＨはエルミート転置、変数の直上に付した＾は推定値であることを示す記号である。雑音共分散行列＾Φ_f ^vを用いれば、下記（４）式により、雑音を抑圧した強調信号である観測信号ベクトルｙ_tfへの目的信号の寄与ベクトルｓ_tf∈Ｃ^Mを推定できる。

なお、上記（４）式において、ベクトル＾ｖ_tfは観測信号ベクトルｙ_tfに含まれる雑音の寄与ベクトルｖ_tfの推定値であるマルチチャネルウィーナフィルタ出力、γ_tf ⁽¹⁾は時間周波数マスク、行列＾Φ_f ^yは観測共分散行列である。本アプローチによる雑音除去の性能は事後確率γ_tf ⁽ⁿ⁾の精度に大きく依存する。

従来の雑音除去技術では、各時間周波数点にて観測信号ベクトルｙ_tfから抽出した信号到来方向を特徴づける到来方向特徴量ベクトルｚ_tfを、「状態Ｓ⁽⁰⁾と対応するクラスタ（目的信号と対応しないクラスタ）」と「状態Ｓ⁽¹⁾と対応するクラスタ（目的信号と対応するクラスタ）」の２つのクラスタにクラスタリングすることで、事後確率γ_tf ⁽ⁿ⁾を算出する。時間周波数成分（ｔ，ｆ）が状態Ｓ⁽ⁿ⁾（ｎ＝０，１）を取る場合の到来方向特徴量ベクトルｚ_tfの分布をｐ(ｚ_tf|Ｓ⁽ⁿ⁾，Θ)、状態の事前分布をＰ(Ｓ⁽ⁿ⁾|Θ)とすれば、到来方向特徴量ベクトルｚ_tfの分布は、下記（５）式のように混合分布でモデル化される。

ただし、上記（５）式における集合Θは、全てのモデルパラメータの集合である。ＥＭアルゴリズムなどにより、下記（６）式の尤度を最大化するパラメータ集合Θを求めることにより、上記（５）式の混合分布を、観測された到来方向特徴量ベクトルｚ_tfの分布に当てはめる。得られたパラメータ集合Θの推定値＾Θを用いれば、事後確率γ_tf ⁽ⁿ⁾を、下記（７）式により算出することができる。

到来方向特徴量ベクトルｚ_tfとしては、下記（８）式により定義される観測信号ベクトルｙ_tfの単位方向ベクトルを用いることができる。なお、下記（８）式において、||*||₂は、２ノルムを表す。

上記（８）式による到来方向特徴量ベクトルｚ_tfは、異なる周波数ビンでは異なる分布形状を持つことから、従来の雑音除去技術では、周波数ビン毎に独立にクラスタリングを行う。この場合、１つの状態が、異なる周波数ビンでは異なるクラスタ番号を持つというパーミュテーション問題が生じる。そこで、１つの状態に対応するクラスタを周波数ビン間で束ねる後処理を行う。

このように、従来の雑音除去技術では、状態Ｓ⁽ⁿ⁾（ｎ＝０，１）に対する到来方向特徴量ベクトルｚ_tfの分布形状に関する事前知識をモデル化していないため、得られた２つのクラスタの各々が状態Ｓ⁽ⁿ⁾（ｎ＝０，１）のどちらに対応するか不明であるというクラスタのラベルの曖昧性が生じる。このため、従来の雑音除去技術では、雑音除去を行うために、クラスタリング後に人手で各クラスタを各状態に対応付ける必要がある。

本発明は、このような背景に鑑みてなされたものであり、クラスタのラベルの曖昧性を解消し、状態Ｓ⁽ⁿ⁾（ｎ＝０，１）と紐づいたクラスタを得るための技術を提供する。これにより、クラスタリング後に人手で各クラスタを各状態に対応付ける必要がなくなり、雑音除去を自動で行うことができるようになる。

これを実現するために、本発明では、「目的信号は特定方向から到来するのに対し、雑音は多方向から一様に到来する」という性質に着目する。この性質により、状態Ｓ⁽⁰⁾では目的信号が存在せず雑音のみが存在するため、信号到来方向は一様である（信号到来方向の局在性が小さい）のに対し、状態Ｓ⁽¹⁾では目的信号が存在するため信号到来方向に偏りがある（信号到来方向の局在性が大きい）。本発明では、この点に着目し、クラスタ内の信号到来方向の局在性の大きさを表す到来方向局在度を用いて、「状態Ｓ⁽⁰⁾と対応するクラスタ（目的信号と対応しないクラスタ）」よりも「状態Ｓ⁽¹⁾と対応するクラスタ（目的信号と対応するクラスタ）」の方がより大きい到来方向局在度を持つと仮定することで、クラスタのラベルの曖昧性を解消する。

＜実施形態＞
以下、本願の開示技術に係る実施形態の一例について、図面を参照して説明する。なお、以下の実施形態により、本願が開示する技術が限定されるものではない。また、以下の各実施形態は、適宜組合せてもよい。

［実施形態１］
開示技術の実施形態１では、複数のマイクロホンで取得した観測信号と、クラスタのラベルの曖昧性が存在する２つのクラスタの事後確率とを用いて、各クラスタの到来方向局在度を周波数ビン毎に算出し、この到来方向局在度の周波数方向への集約値が最大となるクラスタを目的信号と対応するクラスタと判定することで、クラスタのラベルの曖昧性を解消する。

（実施形態１に係るクラスタ判定装置）
図１は、実施形態１に係るクラスタ判定装置の一例を示す図である。実施形態１に係るクラスタ判定装置１００は、Ｍ（≧２）個のマイクロホンで取得した観測信号の時間周波数変換ｙ_tf ⁽¹⁾，・・・，ｙ_tf ^(M)と、「各クラスタが状態Ｓ⁽ⁿ⁾（ｎ＝０，１）のどちらに対応するかが不明である」というクラスタのラベルの曖昧性が存在する２つのクラスタの事後確率β_tf ⁽ⁿ⁾とを受け取って、クラスタのラベルの曖昧性が解消された事後確率γ_tf ⁽ⁿ⁾を出力する。

なお、観測信号の時間周波数変換値ｙ_tf ⁽¹⁾，・・・，ｙ_tf ^(M)をベクトルとしてまとめたものを観測信号ベクトルｙ_tfとして、下記（９）式のように定義する。以下、クラスタのラベルの曖昧性が存在する事後確率をβ_tf ⁽ⁿ⁾、クラスタのラベルの曖昧性が存在しない事後確率をγ_tf ⁽ⁿ⁾と区別する。

クラスタ判定装置１００は、クラスタ内の信号到来方向の局在性の大きさを表す到来方向局在度κ_f ⁽ⁿ⁾を算出し、到来方向局在度κ_f ⁽ⁿ⁾を周波数方向に集約した集約値が最大となるクラスタを目的信号と対応するクラスタと判定する。

クラスタ判定装置１００は、到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）、判定部２０を有する。到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）は、観測信号ベクトルｙ_tfと事後確率β_tf ⁽ⁿ⁾を受け取って、到来方向局在度κ_f ⁽ⁿ⁾を算出し、判定部２０に受け渡す。判定部２０は、事後確率β_tf ⁽ⁿ⁾と到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）からの到来方向局在度κ_f ⁽ⁿ⁾を受け取って、クラスタのラベルの曖昧性が解消された事後確率γ_tf ⁽ⁿ⁾を算出して出力する。

到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）において算出する到来方向局在度としては、特徴量局在度、固有値スパース度、異方度等を用いることができる。以下では、特徴量局在度を用いる到来方向局在度算出部１０Ａを第１例とし、固有値スパース度を用いる到来方向局在度算出部１０Ｂを第２例とし、異方度を用いる到来方向局在度算出部１０Ｃを第３例として説明する。

（第１例の到来方向局在度算出部）
図２は、実施形態１に係る第１例の到来方向局在度算出部を示す図である。第１例の到来方向局在度算出部１０Ａが到来方向局在度として用いる特徴量局在度は、各時間周波数点で抽出された信号到来方向を特徴づける到来方向特徴量の局在性の大きさである。

到来方向局在度算出部１０Ａは、到来方向特徴量算出部１０Ａ−１、特徴量局在度算出部１０Ａ−２を有する。到来方向特徴量算出部１０Ａ−１は、観測信号ベクトルｙ_tfを受け取って、到来方向特徴量ベクトルｚ_tfを算出し、特徴量局在度算出部１０Ａ−２に受け渡す。到来方向特徴量としては、マイクロホン間の位相差や振幅比、マイクロホン間到来時間差、ＤＯＡ（Direction Of Arrival）、観測信号ベクトルｙ_tfの単位方向ベクトル等がある。例えば、観測信号ベクトルｙ_tfの単位方向ベクトルは、下記（１０）式により算出できる。

特徴量局在度算出部１０Ａ−２は、事後確率β_tf ⁽ⁿ⁾と、到来方向特徴量算出部１０Ａ−１からの到来方向特徴量ベクトルｚ_tfを受け取って、特徴量局在度κ_f ⁽ⁿ⁾を算出し、特徴量局在度κ_f ⁽ⁿ⁾を到来方向局在度として出力する。

特徴量局在度κ_f ⁽ⁿ⁾としては、到来方向特徴量ベクトルｚ_tfのクラスタ内分散の逆数、到来方向特徴量ベクトルｚ_tfのクラスタ内エントロピーの符号を反転させたものなどを用いることができる。また、到来方向特徴量ベクトルｚ_tfとして観測信号ベクトルｙ_tfの単位方向ベクトルを用いた場合には、観測された到来方向特徴量ベクトルｚ_tfの分布に複素ワトソン分布を当てはめたときの集中パラメータを特徴量局在度κ_f ⁽ⁿ⁾として用いることもできる。

集中パラメータκ_f ⁽ⁿ⁾は、以下の手順で算出できる。まず、単位方向ベクトルｚ_tfのクラスタ内共分散行列Ｒ_f ⁽ⁿ⁾を、下記（１１）式により算出する。

次に、クラスタ内共分散行列Ｒ_f ⁽ⁿ⁾の固有値分解により最大固有値μ_f ⁽ⁿ⁾を取得する。最後に、集中パラメータκ_f ⁽ⁿ⁾を、下記（１２）式により算出する。

（第２例の到来方向局在度算出部）
図３は、実施形態１に係る第２例の到来方向局在度算出部を示す図である。第２例の到来方向局在度算出部１０Ｂが到来方向局在度として用いる固有値スパース度は、観測信号ベクトルのクラスタ内共分散行列の固有値分布のスパース性の大きさである。一般に、ある信号の共分散行列の各固有値は対応する固有ベクトル方向におけるその信号の分散を表す。したがって、観測信号ベクトルのクラスタ内共分散行列の固有値分布がスパースであるということは、そのクラスタ内において観測信号ベクトルの方向が特定の方向に局在していることを意味する。

観測信号ベクトルの方向は信号到来方向と対応するから、観測信号ベクトルの方向が局在しているということは、信号到来方向が局在していることを意味する。したがって、固有値スパース度を到来方向局在度として用いることができる。

到来方向局在度算出部１０Ｂは、共分散行列算出部１０Ｂ−１、固有値分布算出部１０Ｂ−２、固有値スパース度算出部１０Ｂ−３を有する。

共分散行列算出部１０Ｂ−１は、事後確率β_tf ⁽ⁿ⁾と、観測信号ベクトルｙ_tfを受け取って、クラスタ内共分散行列Φ_f ⁽ⁿ⁾を算出し、固有値分布算出部１０Ｂ−２に受け渡す。クラスタ内共分散行列Φ_f ⁽ⁿ⁾は、下記（１３）式により算出できる。

固有値分布算出部１０Ｂ−２は、共分散行列算出部１０Ｂ−１からクラスタ内共分散行列Φ_f ⁽ⁿ⁾を受け取って、クラスタ内共分散行列Φ_f ⁽ⁿ⁾の固有値λ_f1 ⁽ⁿ⁾，・・・，λ_fM ⁽ⁿ⁾を算出し、下記（１４）式のように定義される固有値分布ベクトルλ_f ⁽ⁿ⁾を固有値スパース度算出部１０Ｂ−３に受け渡す。ここで、λ_f1 ⁽ⁿ⁾，λ_f2 ⁽ⁿ⁾，・・・，λ_fM-1 ⁽ⁿ⁾，λ_fM ⁽ⁿ⁾は、λ_f1 ⁽ⁿ⁾≧λ_f2 ⁽ⁿ⁾≧・・・≧λ_fM-1 ⁽ⁿ⁾≧λ_fM ⁽ⁿ⁾の大小関係が成り立つように降順に並べるものとする。

固有値スパース度算出部１０Ｂ−３は、固有値分布算出部１０Ｂ−２から固有値分布ベクトルλ_f ⁽ⁿ⁾を受け取って、固有値スパース度κ_f ⁽ⁿ⁾を算出して到来方向局在度として出力する。固有値スパース度κ_f ⁽ⁿ⁾としては、最大固有値λ_f1 ⁽ⁿ⁾と最小固有値λ_fM ⁽ⁿ⁾の比（すなわちクラスタ内共分散行列Φ_f ⁽ⁿ⁾の条件数）、零固有値の個数（すなわちクラスタ内共分散行列Φ_f ⁽ⁿ⁾の零空間の次元）、固有値分布ベクトルλ_f ⁽ⁿ⁾の１ノルムと２ノルムの比（すなわちクラスタ内共分散行列Φ_f ⁽ⁿ⁾のトレースノルムとフロベニウスノルムの比）、固有値分布ベクトルλ_f ⁽ⁿ⁾のエントロピーの符号を反転させたもの等を用いることができる。

（第３例の到来方向局在度算出部）
図４は、実施形態１に係る第３例の到来方向局在度算出部を示す図である。第３例の到来方向局在度算出部１０Ｃが到来方向局在度として用いる異方度は、参考文献１及び参考文献２に開示されている行列部分空間の概念に基づき、以下のように定義される。参考文献１は、「N. Ito, E. Vincent, N. Ono, and S. Sagayama, “General Algorithms for Estimating Spectrogram and Transfer Functions of Target Signal for Blind Suppression of Diffuse Noise,” Proc. IEEE International Workshop on Machine Learning for Signal Processing (MLSP), Sep. 2013.」である。参考文献２は、「N. Ito, E. Vincent, T. Nakatani, N. Ono, S. Araki, and S. Sagayama, “Blind Suppression of Nonstationary Diffuse Acoustic Noise Based on Spatial Covariance Matrix Decomposition,” Journal of Signal Processing Systems, vol. 79, no. 2, pp. 145-157, May 2015.」である。

すなわち、Ｍ次エルミート行列全体の集合を実ベクトル空間（行列のなすベクトル空間なので行列空間と呼ぶ）とみなし、あらゆる方向から到来する等方性信号の共分散行列を前述の行列空間の部分空間（以下、等方性行列部分空間と呼ぶ）としてモデル化することで、観測信号のクラスタ内共分散行列Φ_f ⁽ⁿ⁾は、等方性行列部分空間に属する成分行列Φ_f ^||(n)と、異方性行列部分空間（等方性行列部分空間の直交補空間）に属する成分行列Φ_f ^⊥(n)とに直交分解できる。成分行列Φ_f ^⊥(n)及び成分行列Φ_f ^||(n)のノルム比（以下、異方度と呼ぶ）は、下記（１５）式で与えられる。異方度は、各クラスタ内における信号到来方向の異方性の大きさ（等方性の小ささ）を表すので、これを到来方向局在度として用いることができる。

到来方向局在度算出部１０Ｃは、共分散行列算出部１０Ｃ−１、直交分解部１０Ｃ−２、異方度算出部１０Ｃ−３を有する。共分散行列算出部１０Ｃ−１は、事後確率β_tf ⁽ⁿ⁾と、観測信号ベクトルｙ_tfを受け取って、観測信号のクラスタ内共分散行列Φ_f ⁽ⁿ⁾を算出し、直交分解部１０Ｃ−２に受け渡す。クラスタ内共分散行列Φ_f ⁽ⁿ⁾は、下記（１６）式により算出できる。

直交分解部１０Ｃ−２は、共分散行列算出部１０Ｃ−１からクラスタ内共分散行列Φ_f ⁽ⁿ⁾を受け取って、等方性行列部分空間に属する成分行列Φ_f ^||(n)と、異方性行列部分空間に属する成分行列Φ_f ^⊥(n)とを算出し、異方度算出部１０Ｃ−３に受け渡す。等方性行列部分空間に属する成分行列Φ_f ^||(n)と、異方性行列部分空間に属する成分行列Φ_f ^⊥(n)は、それぞれ下記（１７−１）式及び（１７−２）式により算出できる。

ここで、｛Ｑ_f ^(d)｝_d=1 ^Dは等方性行列部分空間の直交基底を表し、｛Ｑ_f ^(d)｝_d=D+1 ^M2（M2はM²を表す）は異方性行列部分空間の直交基底を表し、<Φ_f ⁽ⁿ⁾，Ｑ_f ^(d)>は下記（１８）式により定義されるクラスタ内共分散行列Φ_f ⁽ⁿ⁾と行列Ｑ_f ^(d)の内積を表し、Ｄは等方性行列部分空間の次元を表し、||*||_Fはフロベニウスノルムを表し、trはトレースを表す。

直交基底｛Ｑ_f ^(d)｝_d=1 ^Dとしては、例えば、下記（１９）式により定義されるｂ_f ^(m,n)を（ｍ，ｎ）成分とする行列Ｂ_fのみからなる、下記（２０）式で表現される基底を用いることができる。ここで、sinc(*)はsinc(x):＝sin(x)/xで定義されるsinc関数を表し、ω_fはｆ番目の周波数ビンに対応する角周波数を表し、ｄ_mnはマイクロホンｍとマイクロホンｎとの間の距離を表し、ｃは音速を表す。直交基底｛Ｑ_f ^(d)｝_d=1 ^Dの他の例は、上記の参考文献１及び参考文献２に開示されている。

異方度算出部１０Ｃ−３は、直交分解部１０Ｃ−２から等方性行列部分空間に属する成分行列Φ_f ^||(n)と、異方性行列部分空間に属する成分行列Φ_f ^⊥(n)とを受け取って、到来方向局在度である異方度κ_f ⁽ⁿ⁾を下記（２１）式により算出して出力する。

（判定部）
図５は、実施形態１に係る判定部の一例を示す図である。判定部２０は、周波数方向集約部２０−１、クラスタ置換部２０−２を有する。

周波数方向集約部２０−１は、到来方向局在度κ_f ⁽ⁿ⁾を受け取って、その周波数方向への集約値〜κ⁽ⁿ⁾を算出してクラスタ置換部２０−２に受け渡す。集約値〜κ⁽ⁿ⁾としては、集合｛κ_f ⁽ⁿ⁾｝_f=1 ^Fに属する到来方向局在度κ_f ⁽ⁿ⁾の大きさを表す代表値（相加平均などの平均、中央値、最頻値、最大値、最小値、四分位点など）や、周波数ビン毎にκ_f ⁽ⁿ⁾が最大となる番号ｎに１票を与える投票における得票数等を用いることができる。全周波数ビンにわたる集合｛κ_f ⁽ⁿ⁾｝_f=1 ^Fに対する集約値の代わりに集合｛κ_f ⁽ⁿ⁾|ｆ∈Ｓ｝に対する集約値を用いてもよい。ここで、集合Ｓは、集合｛１，・・・，Ｆ｝の任意の部分集合である。例えば、集約値〜κ⁽ⁿ⁾が得票数である場合には、集約値〜κ⁽ⁿ⁾は、下記（２２）式により計算される。ここで、|*|は集合の元の個数である。

また、例えば、集約値〜κ⁽ⁿ⁾が全周波数ビンにわたる相加平均である場合には、集約値〜κ⁽ⁿ⁾は、下記（２３）式により計算される。

クラスタ置換部２０−２は、事後確率β_tf ⁽ⁿ⁾と、周波数方向集約部２０−１からの集約値〜κ⁽ⁿ⁾を受け取って、クラスタのラベルの曖昧性が解消された事後確率γ_tf ⁽ⁿ⁾を作成し出力する。事後確率γ_tf ⁽ⁿ⁾の作成は、下記（２４）式による。

（実施形態１に係るクラスタ判定処理）
図６は、実施形態１に係るクラスタ判定処理の一例を示すフローチャートである。先ず、到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）が、観測信号ベクトルｙ_tfと事後確率β_tf ⁽ⁿ⁾を受け取って、クラスタ毎かつ周波数ビン毎の到来方向局在度κ_f ⁽ⁿ⁾を算出する（ステップＳ１１）。次に、判定部２０の周波数方向集約部２０−１が、到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）から到来方向局在度κ_f ⁽ⁿ⁾を受け取り、クラスタ毎の集約値〜κ⁽ⁿ⁾を算出する（ステップＳ１２）。次に、判定部２０のクラスタ置換部２０−２が、事後確率β_tf ⁽ⁿ⁾と、周波数方向集約部２０−１からの集約値〜κ⁽ⁿ⁾を受け取って、クラスタのラベルの曖昧性が解消されたクラスタ毎かつ時間周波数成分毎の事後確率γ_tf ⁽ⁿ⁾を作成し出力する（クラスタ置換、ステップＳ１３）。

（実施形態１の変形例１）
到来方向局在度κ_f ⁽ⁿ⁾が既知の場合には、実施形態１に係るクラスタ判定装置は、前記到来方向局在度κ_f ⁽ⁿ⁾を入力として用い、到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）を含まない構成としても良い。すなわち、実施形態１の変形例１に係るクラスタ判定装置は、判定部２０のみからなる。前記判定部２０は、観測信号の時間周波数成分を複数のクラスタへクラスタリングすることで得られたクラスタのラベルの曖昧性が存在する２つのクラスタの事後確率β_tf ⁽ⁿ⁾と到来方向局在度κ_f ⁽ⁿ⁾を受け取って、クラスタのラベルの曖昧性が解消された事後確率γ_tf ⁽ⁿ⁾を作成し出力する。判定部２０における処理は、既に述べた通りであるため詳細な説明を省略する。

なお、到来方向局在度κ_f ⁽ⁿ⁾が既知の場合とは、例えば、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和である混合モデルを前記特徴量ベクトルに当てはめることによって前記クラスタリングを行なう場合に、前記混合モデルを前記特徴量ベクトルに当てはめたときの前記混合モデルのパラメータを用いて到来方向局在度κ_f ⁽ⁿ⁾が得られる場合などである。

（実施形態１の変形例２）
実施形態１に係るクラスタ判定装置は、クラスタのラベルの曖昧性が解消された事後確率γ_tf ⁽ⁿ⁾を出力する代わりに、γ_tf ⁽ⁿ⁾= β_tf ^(Π(n))となる集合｛０，１｝上の置換Πを出力する構成としても良い。前記置換Πは、クラスタのラベルの曖昧性が存在する２つのクラスタのそれぞれが、目的信号と対応するクラスタであるか、目的信号と対応しないクラスタであるかを指定する対応関係である。すなわち、実施形態１の変形例２に係るクラスタ判定装置は、到来方向局在度算出部と判定部からなり、前記判定部は、周波数方向集約部とクラスタ置換部からなり、前記クラスタ置換部は、周波数方向集約部からの集約値〜κ⁽ⁿ⁾を受け取って前記置換Πを決定し出力する。前記置換Πは、「〜κ⁽⁰⁾＜〜κ⁽¹⁾」ならば「Π(０)＝０、Π(１)＝１」、「〜κ⁽⁰⁾≧〜κ⁽¹⁾」ならば「Π(０)＝１、Π(１)＝０」となるように決定すればよい。なお、前記到来方向局在度算出部と前記周波数方向集約部における処理は、既に述べた通りであるため詳細な説明を省略する。

［実施形態２］
実施形態２では、複数のマイクロホンで取得した観測信号とクラスタのラベルの曖昧性が存在する２つのクラスタの事後確率とを用いて各クラスタの到来方向局在度を周波数ビン毎に算出し、周波数ビン毎に到来方向局在度が最大となるクラスタを目的信号と対応するクラスタと判定することで、クラスタのラベルの曖昧性を解消する。

実施形態２では、クラスタのラベルの曖昧性を解消するという効果に加えて、パーミュテーションの曖昧性を解消するという効果もある。事後確率を周波数ビン毎のクラスタリングにより算出する場合、１つの状態が、異なる周波数ビンでは異なるクラスタ番号を持つというパーミュテーションの曖昧性がある。実施形態２では、周波数ビン毎にクラスタのラベルの曖昧性を解消するため、上記パーミュテーションの曖昧性も解消することができる（パーミュテーション解決）。

実施形態２では、実施形態１の場合と同様に、観測信号ベクトルｙ_tfと、クラスタのラベルの曖昧性が存在する２つのクラスタの事後確率β_tf ⁽ⁿ⁾とを受け取って、クラスタのラベルの曖昧性が解消された事後確率γ_tf ⁽ⁿ⁾を出力する。

実施形態２のポイントは、周波数ビン毎に到来方向局在度κ_f ⁽ⁿ⁾が最大となるクラスタを目的信号と対応するクラスタと判定することである。実施形態２は、周波数ビン毎に判定を行う点で、実施形態１とは異なる。これにより、クラスタのラベルの曖昧性を解消する効果に加え、パーミュテーションの曖昧性を解消する効果も得ることができる。

（実施形態２に係るクラスタ判定装置）
図１を参照して、実施形態２に係るクラスタ判定装置の一例を説明する。実施形態２に係るクラスタ判定装置１００Ｄは、到来方向局在度算出部１０Ｄ、判定部２０Ｄを有する。

到来方向局在度算出部１０Ｄは、観測信号ベクトルｙ_tfと、事後確率β_tf ⁽ⁿ⁾を受け取って到来方向局在度κ_f ⁽ⁿ⁾を算出し、判定部２０Ｄに受け渡す。到来方向局在度算出部１０Ｄにおける処理は、実施形態１の場合と同様であるから、詳細な説明は省略する。

判定部２０Ｄは、事後確率β_tf ⁽ⁿ⁾と、到来方向局在度算出部１０Ｄからの到来方向局在度κ_f ⁽ⁿ⁾を受け取って、クラスタのラベルの曖昧性が解消された事後確率γ_tf ⁽ⁿ⁾を作成し出力する。事後確率γ_tf ⁽ⁿ⁾は、下記（２５）式のように、周波数ビン毎に到来方向局在度κ_f ⁽ⁿ⁾を比較して作成する。

（実施形態２に係るクラスタ判定処理）
図７は、実施形態２に係るクラスタ判定処理の一例を示すフローチャートである。先ず、到来方向局在度算出部１０Ｄが、観測信号ベクトルｙ_tfと、事後確率β_tf ⁽ⁿ⁾を受け取ってクラスタ毎かつ周波数ビン毎の到来方向局在度κ_f ⁽ⁿ⁾を算出する（ステップＳ２１）。次に判定部２０Ｄが、事後確率β_tf ⁽ⁿ⁾と、到来方向局在度算出部１０Ｄからの到来方向局在度κ_f ⁽ⁿ⁾を受け取って、クラスタのラベルの曖昧性が解消されたクラスタ毎かつ時間周波数成分毎の事後確率γ_tf ⁽ⁿ⁾を作成し出力する（判定処理、ステップＳ２２）。

（実施形態２の変形例１）
到来方向局在度κ_f ⁽ⁿ⁾が既知の場合には、実施形態２に係るクラスタ判定装置は、前記到来方向局在度κ_f ⁽ⁿ⁾を入力として用い、到来方向局在度算出部１０Ｄを含まない構成としても良い。すなわち、実施形態２の変形例１に係るクラスタ判定装置は、判定部２０Ｄのみからなる。前記判定部２０Ｄは、観測信号の時間周波数成分を複数のクラスタへクラスタリングすることで得られたクラスタのラベルの曖昧性が存在する２つのクラスタの事後確率β_tf ⁽ⁿ⁾と到来方向局在度κ_f ⁽ⁿ⁾を受け取って、クラスタのラベルの曖昧性が解消された事後確率γ_tf ⁽ⁿ⁾を作成し出力する。判定部２０Ｄにおける処理は、既に述べた通りであるため詳細な説明を省略する。

（実施形態２の変形例２）
実施形態２に係るクラスタ判定装置は、クラスタのラベルの曖昧性が解消された事後確率γ_tf ⁽ⁿ⁾を出力する代わりに、γ_tf ⁽ⁿ⁾= β_tf ^(Πｆ(n))となる周波数ビン毎の集合｛０，１｝上の置換Π_ｆを出力する構成としても良い。ただし、前記ΠｆはΠ_ｆを表す。前記置換Π_ｆは、各周波数ビンにおいて、クラスタのラベルの曖昧性が存在する２つのクラスタのそれぞれが、目的信号と対応するクラスタであるか、目的信号と対応しないクラスタであるかを指定する対応関係である。すなわち、実施形態２の変形例２に係るクラスタ判定装置は、到来方向局在度算出部と判定部からなり、前記判定部は到来方向局在度算出部からの到来方向局在度κ_f ⁽ⁿ⁾を受け取って前記置換Π_ｆを作成し出力する。前記置換Π_ｆの作成は、「κ_f ^(０)＜κ_f ^(１)」ならば「Π_ｆ(０)＝０、Π_ｆ(１)＝１」とし、「κ_f ^(０)≧κ_f ^(１)」ならば「Π_ｆ(０)＝１、Π_ｆ(１)＝０」とすればよい。なお、前記到来方向局在度算出部における処理は、既に述べた通りであるため詳細な説明を省略する。

［実施形態３］
実施形態３では、複数のマイクロホンで取得した観測信号を用いて、「目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布よりも、目的信号と対応するクラスタに関する特徴量ベクトルのモデル分布の方が、より大きい到来方向局在度を持つ」という制約を課された、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和である混合モデルを特徴量ベクトルに当てはめることにより、クラスタのラベルの曖昧性を生じさせずにクラスタリングを行なう。

実施形態３では、各周波数ビンにおいてクラスタのラベルの曖昧性を生じさせずにクラスタリングを行なうことができるという効果があり、したがってまた、パーミュテーションの曖昧性を生じさせずにクラスタリングを行なうことができるという効果もある。また、「目的信号は特定方向から到来するのに対し、雑音は多方向から一様に到来する」という目的信号と雑音の性質をモデル化しているため、この性質を事前知識として用いることで、クラスタリングの精度自体を改善する効果もある。

実施形態３では、観測信号ベクトルｙ_tf:＝[ｙ_tf ⁽¹⁾・・・ｙ_tf ^(M)]^Tを受け取ってクラスタのラベルの曖昧性がない事後確率γ_tf ⁽ⁿ⁾を出力する。実施形態３のポイントは、目的信号と雑音の前記性質をモデル化することにより、クラスタのラベルの曖昧性を生じさせずにクラスタリングを行なう点にある。より具体的には、目的信号と雑音の前記性質に基づき、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布よりも、目的信号と対応するクラスタに関する特徴量ベクトルのモデル分布の方が、より大きい到来方向局在度を持つという制約を課された、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和である混合モデルを特徴量ベクトルに当てはめることによりクラスタリングを行なう。

実施形態３では、例えば、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布よりも、目的信号と対応するクラスタに関する特徴量ベクトルのモデル分布の方が、より大きい到来方向局在度を持つように、前記混合モデルのパラメータに制約を課す。例えば、目的信号と対応しないクラスタの到来方向局在度が、その取りうる値の中の最小値を取るように、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布のパラメータに制約を課す。

すなわち、例えば、前記特徴量ベクトルとして観測信号ベクトルの単位方向ベクトルを用い、前記混合モデルとして下記参考文献３において開示されている複素ワトソン分布の重み付き和である混合ワトソン分布を用いる場合、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布である複素ワトソン分布のパラメータである集中パラメータを０に固定することにより、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布の到来方向局在度（特徴量局在度）である集中パラメータが最小値０を取るようにする。また、例えば、前記特徴量ベクトルとして観測信号ベクトルを用い、前記混合モデルとして下記参考文献４において開示されている複素ガウス分布の重み付き和である混合ガウス分布を用いる場合、（イ）目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布である複素ガウス分布のパラメータである共分散行列が単位行列の定数倍であるという制約を課すことにより、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布の到来方向局在度（固有値スパース度）であるクラスタ内共分散行列の最大固有値と最小固有値の比が最小値１を取るようにする、もしくは、（ロ）目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布である複素ガウス分布のパラメータである共分散行列が前記等方性行列部分空間に属するという制約を課すことにより、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布の到来方向局在度（異方度）が最小値０を取るようにする。参考文献３は、「伊藤信貴，荒木章子，木下慶介，中谷智広，“音源位置情報に基づく劣決定ブラインド音源分離のためのパーミュテーションフリークラスタリング法”，電子情報通信学会論文誌Ａ，vol，J97-A，no．4，pp.234-246，Apr，2014.」である。参考文献４は、「N. Ito, S. Araki, T. Yoshioka, and T. Nakatani, “Relaxed Disjointness Based Clustering for Joint Blind Source Separation and Dereverberation,” Proc. International Workshop on Acoustic Signal Enhancement (IWAENC), Sep. 2014.」である。

また、実施形態３では、例えば、混合モデルを用いた反復アルゴリズムに基づくクラスタリングにおいて、各反復において、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布よりも、目的信号と対応するクラスタに関する特徴量ベクトルのモデル分布の方が、より小さい到来方向局在度を持つ場合に、「（イ）目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布のパラメータを調整して、当該モデル分布の到来方向局在度を減少させる」、もしくは、「（ロ）目的信号と対応するクラスタに関する特徴量ベクトルのモデル分布のパラメータを調整して当該モデル分布の到来方向局在度を増加させる」ことにより、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布よりも、目的信号と対応するクラスタに関する特徴量ベクトルのモデル分布の方が、より大きい到来方向局在度を持つようにする。

すなわち、例えば、前記特徴量ベクトルとして観測信号ベクトルの単位方向ベクトルを用い、前記混合モデルとして混合ワトソン分布を用いる場合、各反復において、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布である複素ワトソン分布よりも、目的信号と対応するクラスタに関する特徴量ベクトルのモデル分布である複素ワトソン分布の方が、より小さい到来方向局在度（特徴量局在度）である集中パラメータを持つ場合に、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布である複素ワトソン分布のパラメータである集中パラメータを減少させることにより、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布よりも、目的信号と対応するクラスタに関する特徴量ベクトルのモデル分布の方が、より大きい到来方向局在度である集中パラメータを持つようにする。

以下では、観測信号ベクトルの単位方向ベクトルを特徴量ベクトルとして用い、目的信号と対応しないクラスタに関する前記特徴量ベクトルのモデル分布である複素ワトソン分布の集中パラメータを０に固定した混合ワトソン分布を混合モデルとして用いる場合について、より詳しく説明する。

（実施形態３に係るクラスタリング装置）
図８は、実施形態３に係るクラスタリング装置の一例を示す図である。実施形態３に係るクラスタリング装置２００は、特徴抽出部２１０、クラスタリング部２２０を有する。

特徴抽出部２１０は、観測信号ベクトルｙ_tfを受け取って各時間周波数点における信号到来方向を特徴づける到来方向特徴量ベクトルである下記（２６）式の観測信号ベクトルｙ_tfの単位方向ベクトルｚ_tfを算出し、クラスタリング部２２０に受け渡す。

クラスタリング部２２０は、特徴抽出部２１０から到来方向特徴量ベクトルｚ_tfを受け取り、クラスタのラベルの曖昧性がない事後確率（以下、単に「事後確率」と呼ぶ）γ_tf ⁽ⁿ⁾を出力する。

クラスタリング部２２０は、事後確率更新部２２１、パラメータ更新部２２２、パラメータ記憶部２２３を有する。事後確率更新部２２１は、特徴抽出部２１０からの到来方向特徴量ベクトルｚ_tfと、パラメータ記憶部２２３からのパラメータΘを受け取り、事後確率γ_tf ⁽ⁿ⁾を更新してパラメータ更新部２２２に受け渡すとともに、クラスタリング部２２０での処理が終了した際に事後確率γ_tf ⁽ⁿ⁾を出力する。パラメータΘは、下記（２７）式により定義される。また、事後確率γ_tf ⁽ⁿ⁾は、下記（２８）式により更新される。

ここで、ｐ(ｚ_tf|Ｓ⁽ⁿ⁾，Θ)（ｎ＝０，１）は、状態Ｓ⁽ⁿ⁾に対する到来方向特徴量ベクトルｚ_tfのモデル分布であり、下記（２９−１）式及び（２９−２）式により計算される。ただし、下記（２９−１）式における状態Ｓ⁽⁰⁾に対する分布ｐ(ｚ_tf|Ｓ⁽⁰⁾，Θ)は、複素ワトソン分布の集中パラメータを０と置くことで得られる分布（単位超球面上の一様分布）である。また、下記（２９−２）式におけるＫ(＊，＊，＊)は、クンマー関数を表す。

パラメータ更新部２２２は、事後確率更新部２２１からの事後確率γ_tf ⁽ⁿ⁾と、特徴抽出部２１０からの到来方向特徴量ベクトルｚ_tfを受け取ってパラメータΘを更新し、パラメータ記憶部２２３に受け渡す。パラメータ更新部２２２の詳細は後述する。パラメータ記憶部２２３は、パラメータ更新部２２２からパラメータΘを受け取って記憶し、事後確率更新部２２１での次の処理の際にパラメータΘを受け渡す。

パラメータ更新部２２２は、混合重み更新部２２２−１、特徴量共分散行列算出部２２２−２、平均方向ベクトル更新部２２２−３、集中パラメータ更新部２２２−４を有する。

混合重み更新部２２２−１は、事後確率更新部２２１から事後確率γ_tf ⁽ⁿ⁾を受け取って混合重みα_f ⁽ⁿ⁾を更新し、パラメータ記憶部２２３に受け渡す。混合重みα_f ⁽ⁿ⁾は、下記（３０）式により更新できる。

特徴量共分散行列更新部２２２−２は、特徴抽出部２１０からの到来方向特徴量ベクトルｚ_tfと、事後確率更新部２２１からの事後確率γ_tf ⁽ⁿ⁾とを受け取って、目的信号に対応するクラスタの特徴量共分散行列Ｒ_f ⁽¹⁾を更新し、平均方向ベクトル更新部２２２−３と、集中パラメータ更新部２２２−４に受け渡す。特徴量共分散行列Ｒ_f ⁽¹⁾は、下記（３１）式により更新できる。

平均方向ベクトル更新部２２２−３は、特徴量共分散行列更新部２２２−２から特徴量共分散行列Ｒ_f ⁽¹⁾を受け取って目的信号に対応するクラスタの平均方向ベクトルａ_f ⁽¹⁾を更新し、パラメータ記憶部２２３に受け渡す。平均方向ベクトルａ_f ⁽¹⁾は、特徴量共分散行列Ｒ_f ⁽¹⁾の固有値分解により得られた最大固有値に対応するノルム１の固有ベクトルにより更新すればよい。一方、上記のように、目的信号に対応しないクラスタは単位超球面上の一様分布によりモデル化しているため、目的信号に対応しないクラスタの平均方向ベクトルは更新する必要がない。

集中パラメータ更新部２２２−４は、特徴量共分散行列更新部２２２−２から特徴量共分散行列Ｒ_f ⁽¹⁾を受け取って目的信号に対応するクラスタの集中パラメータκ_f ⁽¹⁾を更新し、パラメータ記憶部２２３に受け渡す。集中パラメータκ_f ⁽¹⁾は、特徴量共分散行列Ｒ_f ⁽¹⁾の固有値分解により得られた最大固有値μ_f ⁽¹⁾を用いて、下記（３２）式により更新すればよい。一方、上記のように、目的信号に対応しないクラスタは単位超球面上の一様分布によりモデル化しているため目的信号に対応しないクラスタの集中パラメータは更新する必要がない。

（実施形態３に係るクラスタリング処理）
図９は、実施形態３に係るクラスタリング処理の一例を示すフローチャートである。先ず、特徴抽出部２１０が、観測信号ベクトルｙ_tfを受け取って時間周波数成分毎の到来方向特徴量ベクトルｚ_tfを抽出する（ステップＳ３１）。次に、事後確率更新部２２１が、特徴抽出部２１０からの到来方向特徴量ベクトルｚ_tfと、パラメータ記憶部２２３からのパラメータΘとを受け取って、クラスタ毎かつ時間周波数成分毎の事後確率γ_tf ⁽ⁿ⁾を更新する（ステップＳ３２）。

次に、混合重み更新部２２２−１が、事後確率更新部２２１から事後確率γ_tf ⁽ⁿ⁾を受け取って、クラスタ毎かつ周波数ビン毎の混合重みを更新する（ステップＳ３３）。次に、特徴量共分散行列更新部２２２−２が、事後確率更新部２２１からの事後確率γ_tf ⁽ⁿ⁾と、特徴抽出部２１０からの到来方向特徴量ベクトルｚ_tfとを受け取って、目的信号に対応するクラスタの周波数ビン毎の特徴量共分散行列Ｒ_f ⁽¹⁾を算出する（ステップＳ３４）。

次に、平均方向ベクトル更新部２２２−３が、特徴量共分散行列更新部２２２−２から特徴量共分散行列Ｒ_f ⁽¹⁾を受け取って、目的信号に対応するクラスタの周波数ビン毎の平均方向ベクトルａ_f ⁽¹⁾を更新する（ステップＳ３５）。次に、集中パラメータ更新部２２２−４が、特徴量共分散行列更新部２２２−２から特徴量共分散行列Ｒ_f ⁽¹⁾を受け取って、目的信号に対応するクラスタの周波数ビン毎の集中パラメータκ_f ⁽¹⁾を更新する（ステップＳ３６）。次に、図示しない終了条件判定部が、終了条件が満たされているかを判定し（ステップＳ３７）、満たされている場合（ステップＳ３７：Ｙｅｓ）は事後確率更新部２２１が事後確率γ_tf ⁽ⁿ⁾を出力し、満たされていない場合（ステップＳ３７：Ｎｏ）は、ステップＳ３２に戻って処理を続ける。なお、終了条件は、例えば、ステップＳ３２〜Ｓ３７のループ処理を規定回数繰り返した、事後確率γ_tf ⁽ⁿ⁾が所定条件を満たした等である。

以上の各実施形態は、任意のクラスタリング技術と組み合わせて利用することができる。

また、クラスタ判定装置１００，１００Ｄ、クラスタリング装置２００において行われる各処理は、全部又は任意の一部が、ＣＰＵ等の処理装置及び処理装置により解析実行されるプログラムにて実現されてもよい。また、クラスタ判定装置１００，１００Ｄ、クラスタリング装置２００において行われる各処理は、ワイヤードロジックによるハードウェアとして実現されてもよい。

また、実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的に行うこともできる。もしくは、実施形態において説明した各処理のうち、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。この他、上述及び図示の処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて適宜変更することができる。

（プログラムについて）
図１０は、プログラムが実行されることにより、実施形態に係るクラスタ判定装置及びクラスタリング装置が実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ（Central Processing Unit）１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。コンピュータ１０００において、これらの各部はバス１０８０によって接続される。

メモリ１０１０は、ＲＯＭ（Read Only Memory）１０１１及びＲＡＭ（Random Access Memory）１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ（Basic Input Output System）等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０３１に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１０４１に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１０４１に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１０５１、キーボード１０５２に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１０６１に接続される。

ハードディスクドライブ１０３１は、例えば、ＯＳ１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわち、クラスタ判定装置１００，１００Ｄ、クラスタリング装置２００の各処理を規定するプログラムは、コンピュータ１０００によって実行される指令が記述されたプログラムモジュール１０９３として、例えばハードディスクドライブ１０３１に記憶される。例えば、クラスタ判定装置１００，１００Ｄ、クラスタリング装置２００における機能構成と同様の情報処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０３１に記憶される。

また、実施形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０３１に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０３１に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。

なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０３１に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１０４１等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワーク（ＬＡＮ（Local Area Network）、ＷＡＮ（Wide Area Network）等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。

上記の各実施形態は、本願が開示する技術に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。

１０Ａ、１０Ｂ、１０Ｃ、１０Ｄ 到来方向局在度算出部
１０Ａ−１ 到来方向特徴量算出部
１０Ａ−２ 特徴量局在度算出部
１０Ｂ−１、１０Ｃ−１ 共分散行列算出部
１０Ｂ−２ 固有値分布算出部
１０Ｂ−３ 固有値スパース度算出部
１０Ｃ−２ 直交分解部
１０Ｃ−３ 異方度算出部
２０、２０Ｄ 判定部
２０−１ 周波数方向集約部
２０−２ クラスタ置換部
１００、１００Ｄ クラスタ判定装置
２００ クラスタリング装置
２１０ 特徴抽出部
２２０ クラスタリング部
２２１ 事後確率更新部
２２２ パラメータ更新部
２２２−１ 混合重み更新部
２２２−２ 特徴量共分散行列更新部
２２２−３ 平均方向ベクトル更新部
２２２−４ 集中パラメータ更新部
２２３ パラメータ記憶部
１０００ コンピュータ
１０１０ メモリ
１０２０ ＣＰＵ

Claims (3)

1. 複数のマイクロホンで目的信号を含む音を観測した複数の観測信号の時間周波数成分を前記目的信号と対応するクラスタと前記目的信号と対応しないクラスタとにクラスタリングするクラスタリング装置であって、
   各観測信号の時間周波数成分からなる観測信号ベクトルに対応する特徴量ベクトルを抽出する特徴抽出部と、
   前記特徴量ベクトルに所定の混合モデルを当てはめることで前記クラスタ毎かつ時間周波数成分毎の事後確率を算出するクラスタリング部と、
   を含み、
   前記混合モデルは、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和で表され、
   前記各クラスタに関する特徴量ベクトルのモデル分布は、前記目的信号と対応しないクラスタに関する前記特徴量ベクトルのモデル分布が、前記目的信号と対応するクラスタに関する前記特徴量ベクトルのモデル分布よりも小さい到来方向局在度を持つという制約が課されて作成されたものである
   ことを特徴とするクラスタリング装置。
2. 複数のマイクロホンで目的信号を含む音を観測した複数の観測信号の時間周波数成分を前記目的信号と対応するクラスタと前記目的信号と対応しないクラスタとにクラスタリングするクラスタリング装置が実行するクラスタリング方法であって、
   各観測信号の時間周波数成分からなる観測信号ベクトルに対応する特徴量ベクトルを抽出する特徴抽出ステップと、
   前記特徴量ベクトルに所定の混合モデルを当てはめることで前記クラスタ毎かつ時間周波数成分毎の事後確率を算出するクラスタリングステップと、
   を含み、
   前記混合モデルは、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和で表され、
   前記各クラスタに関する特徴量ベクトルのモデル分布は、前記目的信号と対応しないクラスタに関する前記特徴量ベクトルのモデル分布が、前記目的信号と対応するクラスタに関する前記特徴量ベクトルのモデル分布よりも小さい到来方向局在度を持つという制約が課されて作成されたものである
   ことを特徴とするクラスタリング方法。
3. 請求項１に記載のクラスタリング装置としてコンピュータを機能させるためのクラスタリングプログラム。

Images