(57)【要約】【課題】複数のマイクロホンにより取得した観測信号の時間周波数成分をクラスタリングすることで得られたクラスタの各々が、目的信号と対応するクラスタであるか、目的信号と対応しないクラスタであるかを判定する。【解決手段】クラスタリング装置２００は、特徴抽出部２１０およびクラスタリング部２２０を有する。特徴抽出部２１０は、各観測信号の時間周波数成分からなる観測信号ベクトルに対応する特徴量ベクトルを抽出する。クラスタリング部２２０は、特徴量ベクトルに所定の混合モデルを当てはめることでクラスタ毎かつ時間周波数成分毎の事後確率を算出する。【選択図】図８
2018189994.tif 000002
【特許請求の範囲】
【請求項１】
複数のマイクロホンで目的信号を含む音を観測した複数の観測信号の時間周波数成分を前記目的信号と対応するクラスタと前記目的信号と対応しないクラスタとにクラスタリングするクラスタリング装置であって、
各観測信号の時間周波数成分からなる観測信号ベクトルに対応する特徴量ベクトルを抽出する特徴抽出部と、
前記特徴量ベクトルに所定の混合モデルを当てはめることで前記クラスタ毎かつ時間周波数成分毎の事後確率を算出するクラスタリング部と、
を含み、
前記混合モデルは、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和で表され、
前記各クラスタに関する特徴量ベクトルのモデル分布には、前記目的信号と対応しないクラスタに関する前記特徴量ベクトルのモデル分布は、前記目的信号と対応するクラスタに関する前記特徴量ベクトルのモデル分布よりも小さい到来方向局在度を持つという制約が課されている
ことを特徴とするクラスタリング装置。
【請求項２】
複数のマイクロホンで目的信号を含む音を観測した複数の観測信号の時間周波数成分を前記目的信号と対応するクラスタと前記目的信号と対応しないクラスタとにクラスタリングするクラスタリング装置が実行するクラスタリング方法であって、
各観測信号の時間周波数成分からなる観測信号ベクトルに対応する特徴量ベクトルを抽出する特徴抽出ステップと、
前記特徴量ベクトルに所定の混合モデルを当てはめることで前記クラスタ毎かつ時間周波数成分毎の事後確率を算出するクラスタリングステップと、
を含み、
前記混合モデルは、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和で表され、
前記各クラスタに関する特徴量ベクトルのモデル分布には、前記目的信号と対応しないクラスタに関する前記特徴量ベクトルのモデル分布は、前記目的信号と対応するクラスタに関する前記特徴量ベクトルのモデル分布よりも小さい到来方向局在度を持つという制約が課されている
ことを特徴とするクラスタリング方法。
【請求項３】
請求項１に記載のクラスタリング装置としてコンピュータを機能させるためのクラスタリングプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、クラスタリング装置、クラスタリング方法及びクラスタリングプログラムに関する。
【背景技術】
【０００２】
複数のマイクロホンにより取得した目的信号及び雑音を含む複数の観測信号から雑音を抑圧して目的信号を取得する技術として雑音除去技術がある。従来の雑音除去技術では、時間周波数領域（例えば短時間フーリエ変換領域）における目的信号のスパース性の仮定の下、「目的信号及び雑音からなる観測信号は、時間周波数領域の疎な点においてのみ目的信号を含み、他の点においては目的信号を含まず雑音のみを含む」というモデルを用いる。従来の雑音除去技術では、このモデルに基づき、観測信号の各時間周波数成分を、目的信号を含むものと含まないものとに分類する。
【０００３】
すなわち、従来の雑音除去技術では、上記モデルに基づき、「観測信号の各時間周波数成分から抽出した信号到来方向を特徴づける到来方向特徴量の分布は、目的信号を含む時間周波数成分においては目的信号の到来方向に対応する分布を持ち、目的信号を含まない時間周波数成分においては雑音の到来方向に対応する分布を持つ」と仮定する。そして、従来の雑音除去技術では、この仮定に基づき、到来方向特徴量に対してクラスタ数２のクラスタリングを適用することで、観測信号の時間周波数成分を、目的信号を含む時間周波数成分からなるクラスタ（以下「目的信号と対応するクラスタ」）と目的信号を含まない時間周波数成分からなるクラスタ（以下「目的信号と対応しないクラスタ」）とにクラスタリングする。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０１３−０５４２５８号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、上記技術では、目的信号と対応するクラスタと目的信号と対応しないクラスタへのクラスタリングはできても、得られた２つのクラスタのうちのどちらが目的信号と対応するクラスタであるかが不明である。そのため、各クラスタが、目的信号と対応するクラスタであるか、目的信号と対応しないクラスタであるかの判定を人為的に行う必要があるという問題がある。
【０００６】
本願が開示する実施形態の一例は、複数のマイクロホンにより取得した観測信号の時間周波数成分をクラスタリングすることで得られたクラスタの各々が、目的信号と対応するクラスタであるか、目的信号と対応しないクラスタであるかを判定することを目的とする。
【課題を解決するための手段】
【０００７】
上述した課題を解決し、目的を達成するために、本発明の本願のクラスタリング装置は、複数のマイクロホンで目的信号を含む音を観測した複数の観測信号の時間周波数成分を前記目的信号と対応するクラスタと前記目的信号と対応しないクラスタとにクラスタリングするクラスタリング装置であって、各観測信号の時間周波数成分からなる観測信号ベクトルに対応する特徴量ベクトルを抽出する特徴抽出部と、前記特徴量ベクトルに所定の混合モデルを当てはめることで前記クラスタ毎かつ時間周波数成分毎の事後確率を算出するクラスタリング部と、を含み、前記混合モデルは、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和で表され、前記各クラスタに関する特徴量ベクトルのモデル分布には、前記目的信号と対応しないクラスタに関する前記特徴量ベクトルのモデル分布は、前記目的信号と対応するクラスタに関する前記特徴量ベクトルのモデル分布よりも小さい到来方向局在度を持つという制約が課されていることを特徴とする。
【発明の効果】
【０００８】
本願が開示する実施形態の一例によれば、例えば、複数のマイクロホンにより取得した観測信号の時間周波数成分をクラスタリングすることで得られたクラスタの各々が、目的信号と対応するクラスタであるか、目的信号と対応しないクラスタであるかを判定することができる。
【図面の簡単な説明】
【０００９】
【図１】図１は、実施形態１（実施形態２）に係るクラスタ判定装置の一例を示す図である。
【図２】図２は、実施形態１（実施形態２）に係る第１例の到来方向局在度算出部を示す図である。
【図３】図３は、実施形態１（実施形態２）に係る第２例の到来方向局在度算出部を示す図である。
【図４】図４は、実施形態１（実施形態２）に係る第３例の到来方向局在度算出部を示す図である。
【図５】図５は、実施形態１に係る判定部の一例を示す図である。
【図６】図６は、実施形態１に係るクラスタ判定処理の一例を示すフローチャートである。
【図７】図７は、実施形態２に係るクラスタ判定処理の一例を示すフローチャートである。
【図８】図８は、実施形態３に係るクラスタリング装置の一例を示す図である。
【図９】図９は、実施形態３に係るクラスタリング処理の一例を示すフローチャートである。
【図１０】図１０は、プログラムが実行されることにより、実施形態に係るクラスタ判定装置及びクラスタリング装置が実現されるコンピュータの一例を示す図である。
【発明を実施するための形態】
【００１０】
以下、本願の開示技術に係る実施形態の一例の説明に先立ち、開示技術の前提及び背景となる従来技術について説明する。従来技術の説明後、本願の開示技術に係る実施形態の一例について説明する。以下、重複する構成及び処理については、後出部分について同一符号を付与して説明を省略する。また、以下、下記（１−１）式〜（１−３）式に示すように、「“＾Ａ”を“Ａの直上に＾を付した記号”」、「“〜Ａ”を“Ａの直上に〜を付した記号”」、「“ＡCB”を“Ａに上付きＢを付し、下付きＣを付した記号”」と同等であるとする。また、以下、記号Ａがスカラーを示す場合には単に“Ａ”と表記し、記号Ａがベクトルを示す場合には“ベクトルＡ”と表記し、記号Ａが行列を示す場合には“行列Ａ”と表記し、記号Ａが集合を示す場合には“集合Ａ”と表記する。
【００１１】
【数１】
2018189994.tif 000003
【００１２】
［従来の雑音除去技術］
複数のマイクロホンにより取得した目的信号及び雑音を含む観測信号を用いて雑音を抑圧した強調信号を取得する従来の雑音除去技術について説明する。
【００１３】
従来の雑音除去技術では、「目的信号は時間周波数平面上の疎な点においてのみパワーを持つ」（スパース性）と仮定する。この仮定は音声信号に対して特に精度よく成立することが知られている。スパース性を仮定すれば、観測信号の各時間周波数成分は、下記（２）式に示す２つの状態のいずれかを取る。
【００１４】
【数２】
2018189994.tif 000004
【００１５】
状態Ｓ(0)では観測信号は雑音のみを含むことに着目すれば、雑音共分散行列＾Φfvは、下記（３）式により推定できる。
【００１６】
【数３】
2018189994.tif 000005
【００１７】
ここで、ベクトルｙtf:＝[ｙtf(1)・・・ｙtf(M)]T∈ＣM（ＣMはＭ次元複素ベクトル空間）は、Ｍ個のマイクロホンにおける観測信号の短時間フーリエ変換（以下、観測信号ベクトルと呼ぶ）、γtf(n)は状態Ｓ(n)（ｎ＝０，１）の検出により得られる各時間周波数点での状態Ｓ(n)の事後確率、ｔ∈｛１，・・・，Ｔ｝は時間フレームの番号、ｆ∈｛１，・・・，Ｆ｝は周波数ビンの番号、Ｔは時間フレームの数、Ｆは周波数ビンの数、上付きのＴは転置、上付きのＨはエルミート転置、変数の直上に付した＾は推定値であることを示す記号である。雑音共分散行列＾Φfvを用いれば、下記（４）式により、雑音を抑圧した強調信号である観測信号ベクトルｙtfへの目的信号の寄与ベクトルｓtf∈ＣMを推定できる。
【００１８】
【数４】
2018189994.tif 000006
【００１９】
なお、上記（４）式において、ベクトル＾ｖtfは観測信号ベクトルｙtfに含まれる雑音の寄与ベクトルｖtfの推定値であるマルチチャネルウィーナフィルタ出力、γtf(1)は時間周波数マスク、行列＾Φfyは観測共分散行列である。本アプローチによる雑音除去の性能は事後確率γtf(n)の精度に大きく依存する。
【００２０】
従来の雑音除去技術では、各時間周波数点にて観測信号ベクトルｙtfから抽出した信号到来方向を特徴づける到来方向特徴量ベクトルｚtfを、「状態Ｓ(0)と対応するクラスタ（目的信号と対応しないクラスタ）」と「状態Ｓ(1)と対応するクラスタ（目的信号と対応するクラスタ）」の２つのクラスタにクラスタリングすることで、事後確率γtf(n)を算出する。時間周波数成分（ｔ，ｆ）が状態Ｓ(n)（ｎ＝０，１）を取る場合の到来方向特徴量ベクトルｚtfの分布をｐ(ｚtf|Ｓ(n)，Θ)、状態の事前分布をＰ(Ｓ(n)|Θ)とすれば、到来方向特徴量ベクトルｚtfの分布は、下記（５）式のように混合分布でモデル化される。
【００２１】
【数５】
2018189994.tif 000007
【００２２】
ただし、上記（５）式における集合Θは、全てのモデルパラメータの集合である。ＥＭアルゴリズムなどにより、下記（６）式の尤度を最大化するパラメータ集合Θを求めることにより、上記（５）式の混合分布を、観測された到来方向特徴量ベクトルｚtfの分布に当てはめる。得られたパラメータ集合Θの推定値＾Θを用いれば、事後確率γtf(n)を、下記（７）式により算出することができる。
【００２３】
【数６】
2018189994.tif 000008
【００２４】
【数７】
2018189994.tif 000009
【００２５】
到来方向特徴量ベクトルｚtfとしては、下記（８）式により定義される観測信号ベクトルｙtfの単位方向ベクトルを用いることができる。なお、下記（８）式において、||*||2は、２ノルムを表す。
【００２６】
【数８】
2018189994.tif 000010
【００２７】
上記（８）式による到来方向特徴量ベクトルｚtfは、異なる周波数ビンでは異なる分布形状を持つことから、従来の雑音除去技術では、周波数ビン毎に独立にクラスタリングを行う。この場合、１つの状態が、異なる周波数ビンでは異なるクラスタ番号を持つというパーミュテーション問題が生じる。そこで、１つの状態に対応するクラスタを周波数ビン間で束ねる後処理を行う。
【００２８】
このように、従来の雑音除去技術では、状態Ｓ(n)（ｎ＝０，１）に対する到来方向特徴量ベクトルｚtfの分布形状に関する事前知識をモデル化していないため、得られた２つのクラスタの各々が状態Ｓ(n)（ｎ＝０，１）のどちらに対応するか不明であるというクラスタのラベルの曖昧性が生じる。このため、従来の雑音除去技術では、雑音除去を行うために、クラスタリング後に人手で各クラスタを各状態に対応付ける必要がある。
【００２９】
本発明は、このような背景に鑑みてなされたものであり、クラスタのラベルの曖昧性を解消し、状態Ｓ(n)（ｎ＝０，１）と紐づいたクラスタを得るための技術を提供する。これにより、クラスタリング後に人手で各クラスタを各状態に対応付ける必要がなくなり、雑音除去を自動で行うことができるようになる。
【００３０】
これを実現するために、本発明では、「目的信号は特定方向から到来するのに対し、雑音は多方向から一様に到来する」という性質に着目する。この性質により、状態Ｓ(0)では目的信号が存在せず雑音のみが存在するため、信号到来方向は一様である（信号到来方向の局在性が小さい）のに対し、状態Ｓ(1)では目的信号が存在するため信号到来方向に偏りがある（信号到来方向の局在性が大きい）。本発明では、この点に着目し、クラスタ内の信号到来方向の局在性の大きさを表す到来方向局在度を用いて、「状態Ｓ(0)と対応するクラスタ（目的信号と対応しないクラスタ）」よりも「状態Ｓ(1)と対応するクラスタ（目的信号と対応するクラスタ）」の方がより大きい到来方向局在度を持つと仮定することで、クラスタのラベルの曖昧性を解消する。
【００３１】
＜実施形態＞
以下、本願の開示技術に係る実施形態の一例について、図面を参照して説明する。なお、以下の実施形態により、本願が開示する技術が限定されるものではない。また、以下の各実施形態は、適宜組合せてもよい。
【００３２】
［実施形態１］
開示技術の実施形態１では、複数のマイクロホンで取得した観測信号と、クラスタのラベルの曖昧性が存在する２つのクラスタの事後確率とを用いて、各クラスタの到来方向局在度を周波数ビン毎に算出し、この到来方向局在度の周波数方向への集約値が最大となるクラスタを目的信号と対応するクラスタと判定することで、クラスタのラベルの曖昧性を解消する。
【００３３】
（実施形態１に係るクラスタ判定装置）
図１は、実施形態１に係るクラスタ判定装置の一例を示す図である。実施形態１に係るクラスタ判定装置１００は、Ｍ（≧２）個のマイクロホンで取得した観測信号の時間周波数変換ｙtf(1)，・・・，ｙtf(M)と、「各クラスタが状態Ｓ(n)（ｎ＝０，１）のどちらに対応するかが不明である」というクラスタのラベルの曖昧性が存在する２つのクラスタの事後確率βtf(n)とを受け取って、クラスタのラベルの曖昧性が解消された事後確率γtf(n)を出力する。
【００３４】
なお、観測信号の時間周波数変換値ｙtf(1)，・・・，ｙtf(M)をベクトルとしてまとめたものを観測信号ベクトルｙtfとして、下記（９）式のように定義する。以下、クラスタのラベルの曖昧性が存在する事後確率をβtf(n)、クラスタのラベルの曖昧性が存在しない事後確率をγtf(n)と区別する。
【００３５】
【数９】
2018189994.tif 000011
【００３６】
クラスタ判定装置１００は、クラスタ内の信号到来方向の局在性の大きさを表す到来方向局在度κf(n)を算出し、到来方向局在度κf(n)を周波数方向に集約した集約値が最大となるクラスタを目的信号と対応するクラスタと判定する。
【００３７】
クラスタ判定装置１００は、到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）、判定部２０を有する。到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）は、観測信号ベクトルｙtfと事後確率βtf(n)を受け取って、到来方向局在度κf(n)を算出し、判定部２０に受け渡す。判定部２０は、事後確率βtf(n)と到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）からの到来方向局在度κf(n)を受け取って、クラスタのラベルの曖昧性が解消された事後確率γtf(n)を算出して出力する。
【００３８】
到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）において算出する到来方向局在度としては、特徴量局在度、固有値スパース度、異方度等を用いることができる。以下では、特徴量局在度を用いる到来方向局在度算出部１０Ａを第１例とし、固有値スパース度を用いる到来方向局在度算出部１０Ｂを第２例とし、異方度を用いる到来方向局在度算出部１０Ｃを第３例として説明する。
【００３９】
（第１例の到来方向局在度算出部）
図２は、実施形態１に係る第１例の到来方向局在度算出部を示す図である。第１例の到来方向局在度算出部１０Ａが到来方向局在度として用いる特徴量局在度は、各時間周波数点で抽出された信号到来方向を特徴づける到来方向特徴量の局在性の大きさである。
【００４０】
到来方向局在度算出部１０Ａは、到来方向特徴量算出部１０Ａ−１、特徴量局在度算出部１０Ａ−２を有する。到来方向特徴量算出部１０Ａ−１は、観測信号ベクトルｙtfを受け取って、到来方向特徴量ベクトルｚtfを算出し、特徴量局在度算出部１０Ａ−２に受け渡す。到来方向特徴量としては、マイクロホン間の位相差や振幅比、マイクロホン間到来時間差、ＤＯＡ（Direction  Of  Arrival）、観測信号ベクトルｙtfの単位方向ベクトル等がある。例えば、観測信号ベクトルｙtfの単位方向ベクトルは、下記（１０）式により算出できる。
【００４１】
【数１０】
2018189994.tif 000012
【００４２】
特徴量局在度算出部１０Ａ−２は、事後確率βtf(n)と、到来方向特徴量算出部１０Ａ−１からの到来方向特徴量ベクトルｚtfを受け取って、特徴量局在度κf(n)を算出し、特徴量局在度κf(n)を到来方向局在度として出力する。
【００４３】
特徴量局在度κf(n)としては、到来方向特徴量ベクトルｚtfのクラスタ内分散の逆数、到来方向特徴量ベクトルｚtfのクラスタ内エントロピーの符号を反転させたものなどを用いることができる。また、到来方向特徴量ベクトルｚtfとして観測信号ベクトルｙtfの単位方向ベクトルを用いた場合には、観測された到来方向特徴量ベクトルｚtfの分布に複素ワトソン分布を当てはめたときの集中パラメータを特徴量局在度κf(n)として用いることもできる。
【００４４】
集中パラメータκf(n)は、以下の手順で算出できる。まず、単位方向ベクトルｚtfのクラスタ内共分散行列Ｒf(n)を、下記（１１）式により算出する。
【００４５】
【数１１】
2018189994.tif 000013
【００４６】
次に、クラスタ内共分散行列Ｒf(n)の固有値分解により最大固有値μf(n)を取得する。最後に、集中パラメータκf(n)を、下記（１２）式により算出する。
【００４７】
【数１２】
2018189994.tif 000014
【００４８】
（第２例の到来方向局在度算出部）
図３は、実施形態１に係る第２例の到来方向局在度算出部を示す図である。第２例の到来方向局在度算出部１０Ｂが到来方向局在度として用いる固有値スパース度は、観測信号ベクトルのクラスタ内共分散行列の固有値分布のスパース性の大きさである。一般に、ある信号の共分散行列の各固有値は対応する固有ベクトル方向におけるその信号の分散を表す。したがって、観測信号ベクトルのクラスタ内共分散行列の固有値分布がスパースであるということは、そのクラスタ内において観測信号ベクトルの方向が特定の方向に局在していることを意味する。
【００４９】
観測信号ベクトルの方向は信号到来方向と対応するから、観測信号ベクトルの方向が局在しているということは、信号到来方向が局在していることを意味する。したがって、固有値スパース度を到来方向局在度として用いることができる。
【００５０】
到来方向局在度算出部１０Ｂは、共分散行列算出部１０Ｂ−１、固有値分布算出部１０Ｂ−２、固有値スパース度算出部１０Ｂ−３を有する。
【００５１】
共分散行列算出部１０Ｂ−１は、事後確率βtf(n)と、観測信号ベクトルｙtfを受け取って、クラスタ内共分散行列Φf(n)を算出し、固有値分布算出部１０Ｂ−２に受け渡す。クラスタ内共分散行列Φf(n)は、下記（１３）式により算出できる。
【００５２】
【数１３】
2018189994.tif 000015
【００５３】
固有値分布算出部１０Ｂ−２は、共分散行列算出部１０Ｂ−１からクラスタ内共分散行列Φf(n)を受け取って、クラスタ内共分散行列Φf(n)の固有値λf1(n)，・・・，λfM(n)を算出し、下記（１４）式のように定義される固有値分布ベクトルλf(n)を固有値スパース度算出部１０Ｂ−３に受け渡す。ここで、λf1(n)，λf2(n)，・・・，λfM-1(n)，λfM(n)は、λf1(n)≧λf2(n)≧・・・≧λfM-1(n)≧λfM(n)の大小関係が成り立つように降順に並べるものとする。
【００５４】
【数１４】
2018189994.tif 000016
【００５５】
固有値スパース度算出部１０Ｂ−３は、固有値分布算出部１０Ｂ−２から固有値分布ベクトルλf(n)を受け取って、固有値スパース度κf(n)を算出して到来方向局在度として出力する。固有値スパース度κf(n)としては、最大固有値λf1(n)と最小固有値λfM(n)の比（すなわちクラスタ内共分散行列Φf(n)の条件数）、零固有値の個数（すなわちクラスタ内共分散行列Φf(n)の零空間の次元）、固有値分布ベクトルλf(n)の１ノルムと２ノルムの比（すなわちクラスタ内共分散行列Φf(n)のトレースノルムとフロベニウスノルムの比）、固有値分布ベクトルλf(n)のエントロピーの符号を反転させたもの等を用いることができる。
【００５６】
（第３例の到来方向局在度算出部）
図４は、実施形態１に係る第３例の到来方向局在度算出部を示す図である。第３例の到来方向局在度算出部１０Ｃが到来方向局在度として用いる異方度は、参考文献１及び参考文献２に開示されている行列部分空間の概念に基づき、以下のように定義される。参考文献１は、「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(n)は、等方性行列部分空間に属する成分行列Φf||(n)と、異方性行列部分空間（等方性行列部分空間の直交補空間）に属する成分行列Φf⊥(n)とに直交分解できる。成分行列Φf⊥(n)及び成分行列Φf||(n)のノルム比（以下、異方度と呼ぶ）は、下記（１５）式で与えられる。異方度は、各クラスタ内における信号到来方向の異方性の大きさ（等方性の小ささ）を表すので、これを到来方向局在度として用いることができる。
【００５８】
【数１５】
2018189994.tif 000017
【００５９】
到来方向局在度算出部１０Ｃは、共分散行列算出部１０Ｃ−１、直交分解部１０Ｃ−２、異方度算出部１０Ｃ−３を有する。共分散行列算出部１０Ｃ−１は、事後確率βtf(n)と、観測信号ベクトルｙtfを受け取って、観測信号のクラスタ内共分散行列Φf(n)を算出し、直交分解部１０Ｃ−２に受け渡す。クラスタ内共分散行列Φf(n)は、下記（１６）式により算出できる。
【００６０】
【数１６】
2018189994.tif 000018
【００６１】
直交分解部１０Ｃ−２は、共分散行列算出部１０Ｃ−１からクラスタ内共分散行列Φf(n)を受け取って、等方性行列部分空間に属する成分行列Φf||(n)と、異方性行列部分空間に属する成分行列Φf⊥(n)とを算出し、異方度算出部１０Ｃ−３に受け渡す。等方性行列部分空間に属する成分行列Φf||(n)と、異方性行列部分空間に属する成分行列Φf⊥(n)は、それぞれ下記（１７−１）式及び（１７−２）式により算出できる。
【００６２】
【数１７】
2018189994.tif 000019
【００６３】
ここで、｛Ｑf(d)｝d=1Dは等方性行列部分空間の直交基底を表し、｛Ｑf(d)｝d=D+1M2（M2はM2を表す）は異方性行列部分空間の直交基底を表し、<Φf(n)，Ｑf(d)>は下記（１８）式により定義されるクラスタ内共分散行列Φf(n)と行列Ｑf(d)の内積を表し、Ｄは等方性行列部分空間の次元を表し、||*||Fはフロベニウスノルムを表し、trはトレースを表す。
【００６４】
【数１８】
2018189994.tif 000020
【００６５】
直交基底｛Ｑf(d)｝d=1Dとしては、例えば、下記（１９）式により定義されるｂf(m,n)を（ｍ，ｎ）成分とする行列Ｂfのみからなる、下記（２０）式で表現される基底を用いることができる。ここで、sinc(*)はsinc(x):＝sin(x)/xで定義されるsinc関数を表し、ωfはｆ番目の周波数ビンに対応する角周波数を表し、ｄmnはマイクロホンｍとマイクロホンｎとの間の距離を表し、ｃは音速を表す。直交基底｛Ｑf(d)｝d=1Dの他の例は、上記の参考文献１及び参考文献２に開示されている。
【００６６】
【数１９】
2018189994.tif 000021
【００６７】
【数２０】
2018189994.tif 000022
【００６８】
異方度算出部１０Ｃ−３は、直交分解部１０Ｃ−２から等方性行列部分空間に属する成分行列Φf||(n)と、異方性行列部分空間に属する成分行列Φf⊥(n)とを受け取って、到来方向局在度である異方度κf(n)を下記（２１）式により算出して出力する。
【００６９】
【数２１】
2018189994.tif 000023
【００７０】
（判定部）
図５は、実施形態１に係る判定部の一例を示す図である。判定部２０は、周波数方向集約部２０−１、クラスタ置換部２０−２を有する。
【００７１】
周波数方向集約部２０−１は、到来方向局在度κf(n)を受け取って、その周波数方向への集約値〜κ(n)を算出してクラスタ置換部２０−２に受け渡す。集約値〜κ(n)としては、集合｛κf(n)｝f=1Fに属する到来方向局在度κf(n)の大きさを表す代表値（相加平均などの平均、中央値、最頻値、最大値、最小値、四分位点など）や、周波数ビン毎にκf(n)が最大となる番号ｎに１票を与える投票における得票数等を用いることができる。全周波数ビンにわたる集合｛κf(n)｝f=1Fに対する集約値の代わりに集合｛κf(n)|ｆ∈Ｓ｝に対する集約値を用いてもよい。ここで、集合Ｓは、集合｛１，・・・，Ｆ｝の任意の部分集合である。例えば、集約値〜κ(n)が得票数である場合には、集約値〜κ(n)は、下記（２２）式により計算される。ここで、|*|は集合の元の個数である。
【００７２】
【数２２】
2018189994.tif 000024
【００７３】
また、例えば、集約値〜κ(n)が全周波数ビンにわたる相加平均である場合には、集約値〜κ(n)は、下記（２３）式により計算される。
【００７４】
【数２３】
2018189994.tif 000025
【００７５】
クラスタ置換部２０−２は、事後確率βtf(n)と、周波数方向集約部２０−１からの集約値〜κ(n)を受け取って、クラスタのラベルの曖昧性が解消された事後確率γtf(n)を作成し出力する。事後確率γtf(n)の作成は、下記（２４）式による。
【００７６】
【数２４】
2018189994.tif 000026
【００７７】
（実施形態１に係るクラスタ判定処理）
図６は、実施形態１に係るクラスタ判定処理の一例を示すフローチャートである。先ず、到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）が、観測信号ベクトルｙtfと事後確率βtf(n)を受け取って、クラスタ毎かつ周波数ビン毎の到来方向局在度κf(n)を算出する（ステップＳ１１）。次に、判定部２０の周波数方向集約部２０−１が、到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）から到来方向局在度κf(n)を受け取り、クラスタ毎の集約値〜κ(n)を算出する（ステップＳ１２）。次に、判定部２０のクラスタ置換部２０−２が、事後確率βtf(n)と、周波数方向集約部２０−１からの集約値〜κ(n)を受け取って、クラスタのラベルの曖昧性が解消されたクラスタ毎かつ時間周波数成分毎の事後確率γtf(n)を作成し出力する（クラスタ置換、ステップＳ１３）。
【００７８】
（実施形態１の変形例１）
到来方向局在度κf(n)が既知の場合には、実施形態１に係るクラスタ判定装置は、前記到来方向局在度κf(n)を入力として用い、到来方向局在度算出部１０Ａ（１０Ｂ，１０Ｃ）を含まない構成としても良い。すなわち、実施形態１の変形例１に係るクラスタ判定装置は、判定部２０のみからなる。前記判定部２０は、観測信号の時間周波数成分を複数のクラスタへクラスタリングすることで得られたクラスタのラベルの曖昧性が存在する２つのクラスタの事後確率βtf(n)と到来方向局在度κf(n)を受け取って、クラスタのラベルの曖昧性が解消された事後確率γtf(n)を作成し出力する。判定部２０における処理は、既に述べた通りであるため詳細な説明を省略する。
【００７９】
なお、到来方向局在度κf(n)が既知の場合とは、例えば、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和である混合モデルを前記特徴量ベクトルに当てはめることによって前記クラスタリングを行なう場合に、前記混合モデルを前記特徴量ベクトルに当てはめたときの前記混合モデルのパラメータを用いて到来方向局在度κf(n)が得られる場合などである。
【００８０】
（実施形態１の変形例２）
実施形態１に係るクラスタ判定装置は、クラスタのラベルの曖昧性が解消された事後確率γtf(n)を出力する代わりに、γtf(n)=  βtf(Π(n))となる集合｛０，１｝上の置換Πを出力する構成としても良い。前記置換Πは、クラスタのラベルの曖昧性が存在する２つのクラスタのそれぞれが、目的信号と対応するクラスタであるか、目的信号と対応しないクラスタであるかを指定する対応関係である。すなわち、実施形態１の変形例２に係るクラスタ判定装置は、到来方向局在度算出部と判定部からなり、前記判定部は、周波数方向集約部とクラスタ置換部からなり、前記クラスタ置換部は、周波数方向集約部からの集約値〜κ(n)を受け取って前記置換Πを決定し出力する。前記置換Πは、「〜κ(0)＜〜κ(1)」ならば「Π(０)＝０、Π(１)＝１」、「〜κ(0)≧〜κ(1)」ならば「Π(０)＝１、Π(１)＝０」となるように決定すればよい。なお、前記到来方向局在度算出部と前記周波数方向集約部における処理は、既に述べた通りであるため詳細な説明を省略する。
【００８１】
［実施形態２］
実施形態２では、複数のマイクロホンで取得した観測信号とクラスタのラベルの曖昧性が存在する２つのクラスタの事後確率とを用いて各クラスタの到来方向局在度を周波数ビン毎に算出し、周波数ビン毎に到来方向局在度が最大となるクラスタを目的信号と対応するクラスタと判定することで、クラスタのラベルの曖昧性を解消する。
【００８２】
実施形態２では、クラスタのラベルの曖昧性を解消するという効果に加えて、パーミュテーションの曖昧性を解消するという効果もある。事後確率を周波数ビン毎のクラスタリングにより算出する場合、１つの状態が、異なる周波数ビンでは異なるクラスタ番号を持つというパーミュテーションの曖昧性がある。実施形態２では、周波数ビン毎にクラスタのラベルの曖昧性を解消するため、上記パーミュテーションの曖昧性も解消することができる（パーミュテーション解決）。
【００８３】
実施形態２では、実施形態１の場合と同様に、観測信号ベクトルｙtfと、クラスタのラベルの曖昧性が存在する２つのクラスタの事後確率βtf(n)とを受け取って、クラスタのラベルの曖昧性が解消された事後確率γtf(n)を出力する。
【００８４】
実施形態２のポイントは、周波数ビン毎に到来方向局在度κf(n)が最大となるクラスタを目的信号と対応するクラスタと判定することである。実施形態２は、周波数ビン毎に判定を行う点で、実施形態１とは異なる。これにより、クラスタのラベルの曖昧性を解消する効果に加え、パーミュテーションの曖昧性を解消する効果も得ることができる。
【００８５】
（実施形態２に係るクラスタ判定装置）
図１を参照して、実施形態２に係るクラスタ判定装置の一例を説明する。実施形態２に係るクラスタ判定装置１００Ｄは、到来方向局在度算出部１０Ｄ、判定部２０Ｄを有する。
【００８６】
到来方向局在度算出部１０Ｄは、観測信号ベクトルｙtfと、事後確率βtf(n)を受け取って到来方向局在度κf(n)を算出し、判定部２０Ｄに受け渡す。到来方向局在度算出部１０Ｄにおける処理は、実施形態１の場合と同様であるから、詳細な説明は省略する。
【００８７】
判定部２０Ｄは、事後確率βtf(n)と、到来方向局在度算出部１０Ｄからの到来方向局在度κf(n)を受け取って、クラスタのラベルの曖昧性が解消された事後確率γtf(n)を作成し出力する。事後確率γtf(n)は、下記（２５）式のように、周波数ビン毎に到来方向局在度κf(n)を比較して作成する。
【００８８】
【数２５】
2018189994.tif 000027
【００８９】
（実施形態２に係るクラスタ判定処理）
図７は、実施形態２に係るクラスタ判定処理の一例を示すフローチャートである。先ず、到来方向局在度算出部１０Ｄが、観測信号ベクトルｙtfと、事後確率βtf(n)を受け取ってクラスタ毎かつ周波数ビン毎の到来方向局在度κf(n)を算出する（ステップＳ２１）。次に判定部２０Ｄが、事後確率βtf(n)と、到来方向局在度算出部１０Ｄからの到来方向局在度κf(n)を受け取って、クラスタのラベルの曖昧性が解消されたクラスタ毎かつ時間周波数成分毎の事後確率γtf(n)を作成し出力する（判定処理、ステップＳ２２）。
【００９０】
（実施形態２の変形例１）
到来方向局在度κf(n)が既知の場合には、実施形態２に係るクラスタ判定装置は、前記到来方向局在度κf(n)を入力として用い、到来方向局在度算出部１０Ｄを含まない構成としても良い。すなわち、実施形態２の変形例１に係るクラスタ判定装置は、判定部２０Ｄのみからなる。前記判定部２０Ｄは、観測信号の時間周波数成分を複数のクラスタへクラスタリングすることで得られたクラスタのラベルの曖昧性が存在する２つのクラスタの事後確率βtf(n)と到来方向局在度κf(n)を受け取って、クラスタのラベルの曖昧性が解消された事後確率γtf(n)を作成し出力する。判定部２０Ｄにおける処理は、既に述べた通りであるため詳細な説明を省略する。
【００９１】
（実施形態２の変形例２）
実施形態２に係るクラスタ判定装置は、クラスタのラベルの曖昧性が解消された事後確率γtf(n)を出力する代わりに、γtf(n)=  βtf(Πｆ(n))となる周波数ビン毎の集合｛０，１｝上の置換Πｆを出力する構成としても良い。ただし、前記ΠｆはΠｆを表す。前記置換Πｆは、各周波数ビンにおいて、クラスタのラベルの曖昧性が存在する２つのクラスタのそれぞれが、目的信号と対応するクラスタであるか、目的信号と対応しないクラスタであるかを指定する対応関係である。すなわち、実施形態２の変形例２に係るクラスタ判定装置は、到来方向局在度算出部と判定部からなり、前記判定部は到来方向局在度算出部からの到来方向局在度κf(n)を受け取って前記置換Πｆを作成し出力する。前記置換Πｆの作成は、「κf(０)＜κf(１)」ならば「Πｆ(０)＝０、Πｆ(１)＝１」とし、「κf(０)≧κf(１)」ならば「Πｆ(０)＝１、Πｆ(１)＝０」とすればよい。なお、前記到来方向局在度算出部における処理は、既に述べた通りであるため詳細な説明を省略する。
【００９２】
［実施形態３］
実施形態３では、複数のマイクロホンで取得した観測信号を用いて、「目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布よりも、目的信号と対応するクラスタに関する特徴量ベクトルのモデル分布の方が、より大きい到来方向局在度を持つ」という制約を課された、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和である混合モデルを特徴量ベクトルに当てはめることにより、クラスタのラベルの曖昧性を生じさせずにクラスタリングを行なう。
【００９３】
実施形態３では、各周波数ビンにおいてクラスタのラベルの曖昧性を生じさせずにクラスタリングを行なうことができるという効果があり、したがってまた、パーミュテーションの曖昧性を生じさせずにクラスタリングを行なうことができるという効果もある。また、「目的信号は特定方向から到来するのに対し、雑音は多方向から一様に到来する」という目的信号と雑音の性質をモデル化しているため、この性質を事前知識として用いることで、クラスタリングの精度自体を改善する効果もある。
【００９４】
実施形態３では、観測信号ベクトルｙtf:＝[ｙtf(1)・・・ｙtf(M)]Tを受け取ってクラスタのラベルの曖昧性がない事後確率γtf(n)を出力する。実施形態３のポイントは、目的信号と雑音の前記性質をモデル化することにより、クラスタのラベルの曖昧性を生じさせずにクラスタリングを行なう点にある。より具体的には、目的信号と雑音の前記性質に基づき、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布よりも、目的信号と対応するクラスタに関する特徴量ベクトルのモデル分布の方が、より大きい到来方向局在度を持つという制約を課された、各クラスタに関する特徴量ベクトルのモデル分布の重み付き和である混合モデルを特徴量ベクトルに当てはめることによりクラスタリングを行なう。
【００９５】
実施形態３では、例えば、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布よりも、目的信号と対応するクラスタに関する特徴量ベクトルのモデル分布の方が、より大きい到来方向局在度を持つように、前記混合モデルのパラメータに制約を課す。例えば、目的信号と対応しないクラスタの到来方向局在度が、その取りうる値の中の最小値を取るように、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布のパラメータに制約を課す。
【００９６】
すなわち、例えば、前記特徴量ベクトルとして観測信号ベクトルの単位方向ベクトルを用い、前記混合モデルとして下記参考文献３において開示されている複素ワトソン分布の重み付き和である混合ワトソン分布を用いる場合、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布である複素ワトソン分布のパラメータである集中パラメータを０に固定することにより、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布の到来方向局在度（特徴量局在度）である集中パラメータが最小値０を取るようにする。また、例えば、前記特徴量ベクトルとして観測信号ベクトルを用い、前記混合モデルとして下記参考文献４において開示されている複素ガウス分布の重み付き和である混合ガウス分布を用いる場合、（イ）目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布である複素ガウス分布のパラメータである共分散行列が単位行列の定数倍であるという制約を課すことにより、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布の到来方向局在度（固有値スパース度）であるクラスタ内共分散行列の最大固有値と最小固有値の比が最小値１を取るようにする、もしくは、（ロ）目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布である複素ガウス分布のパラメータである共分散行列が前記等方性行列部分空間に属するという制約を課すことにより、目的信号と対応しないクラスタに関する特徴量ベクトルのモデル分布の到来方向局在度（異方度）が最小値０を取るようにする。参考文献３は、「伊藤信貴，荒木章子，木下慶介，中谷智広，“音源位置情報に基づく劣決定ブラインド音源分離のためのパーミュテーションフリークラスタリング法”，電子情報通信学会論文誌Ａ，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を算出し、クラスタリング部２２０に受け渡す。
【０１０２】
【数２６】
2018189994.tif 000028
【０１０３】
クラスタリング部２２０は、特徴抽出部２１０から到来方向特徴量ベクトルｚtfを受け取り、クラスタのラベルの曖昧性がない事後確率（以下、単に「事後確率」と呼ぶ）γtf(n)を出力する。
【０１０４】
クラスタリング部２２０は、事後確率更新部２２１、パラメータ更新部２２２、パラメータ記憶部２２３を有する。事後確率更新部２２１は、特徴抽出部２１０からの到来方向特徴量ベクトルｚtfと、パラメータ記憶部２２３からのパラメータΘを受け取り、事後確率γtf(n)を更新してパラメータ更新部２２２に受け渡すとともに、クラスタリング部２２０での処理が終了した際に事後確率γtf(n)を出力する。パラメータΘは、下記（２７）式により定義される。また、事後確率γtf(n)は、下記（２８）式により更新される。
【０１０５】
【数２７】
2018189994.tif 000029
【０１０６】
【数２８】
2018189994.tif 000030
【０１０７】
ここで、ｐ(ｚtf|Ｓ(n)，Θ)（ｎ＝０，１）は、状態Ｓ(n)に対する到来方向特徴量ベクトルｚtfのモデル分布であり、下記（２９−１）式及び（２９−２）式により計算される。ただし、下記（２９−１）式における状態Ｓ(0)に対する分布ｐ(ｚtf|Ｓ(0)，Θ)は、複素ワトソン分布の集中パラメータを０と置くことで得られる分布（単位超球面上の一様分布）である。また、下記（２９−２）式におけるＫ(＊，＊，＊)は、クンマー関数を表す。
【０１０８】
【数２９】
2018189994.tif 000031
【０１０９】
パラメータ更新部２２２は、事後確率更新部２２１からの事後確率γtf(n)と、特徴抽出部２１０からの到来方向特徴量ベクトルｚtfを受け取ってパラメータΘを更新し、パラメータ記憶部２２３に受け渡す。パラメータ更新部２２２の詳細は後述する。パラメータ記憶部２２３は、パラメータ更新部２２２からパラメータΘを受け取って記憶し、事後確率更新部２２１での次の処理の際にパラメータΘを受け渡す。
【０１１０】
パラメータ更新部２２２は、混合重み更新部２２２−１、特徴量共分散行列算出部２２２−２、平均方向ベクトル更新部２２２−３、集中パラメータ更新部２２２−４を有する。
【０１１１】
混合重み更新部２２２−１は、事後確率更新部２２１から事後確率γtf(n)を受け取って混合重みαf(n)を更新し、パラメータ記憶部２２３に受け渡す。混合重みαf(n)は、下記（３０）式により更新できる。
【０１１２】
【数３０】
2018189994.tif 000032
【０１１３】
特徴量共分散行列更新部２２２−２は、特徴抽出部２１０からの到来方向特徴量ベクトルｚtfと、事後確率更新部２２１からの事後確率γtf(n)とを受け取って、目的信号に対応するクラスタの特徴量共分散行列Ｒf(1)を更新し、平均方向ベクトル更新部２２２−３と、集中パラメータ更新部２２２−４に受け渡す。特徴量共分散行列Ｒf(1)は、下記（３１）式により更新できる。
【０１１４】
【数３１】
2018189994.tif 000033
【０１１５】
平均方向ベクトル更新部２２２−３は、特徴量共分散行列更新部２２２−２から特徴量共分散行列Ｒf(1)を受け取って目的信号に対応するクラスタの平均方向ベクトルａf(1)を更新し、パラメータ記憶部２２３に受け渡す。平均方向ベクトルａf(1)は、特徴量共分散行列Ｒf(1)の固有値分解により得られた最大固有値に対応するノルム１の固有ベクトルにより更新すればよい。一方、上記のように、目的信号に対応しないクラスタは単位超球面上の一様分布によりモデル化しているため、目的信号に対応しないクラスタの平均方向ベクトルは更新する必要がない。
【０１１６】
集中パラメータ更新部２２２−４は、特徴量共分散行列更新部２２２−２から特徴量共分散行列Ｒf(1)を受け取って目的信号に対応するクラスタの集中パラメータκf(1)を更新し、パラメータ記憶部２２３に受け渡す。集中パラメータκf(1)は、特徴量共分散行列Ｒf(1)の固有値分解により得られた最大固有値μf(1)を用いて、下記（３２）式により更新すればよい。一方、上記のように、目的信号に対応しないクラスタは単位超球面上の一様分布によりモデル化しているため目的信号に対応しないクラスタの集中パラメータは更新する必要がない。
【０１１７】
【数３２】
2018189994.tif 000034
【０１１８】
（実施形態３に係るクラスタリング処理）
図９は、実施形態３に係るクラスタリング処理の一例を示すフローチャートである。先ず、特徴抽出部２１０が、観測信号ベクトルｙtfを受け取って時間周波数成分毎の到来方向特徴量ベクトルｚtfを抽出する（ステップＳ３１）。次に、事後確率更新部２２１が、特徴抽出部２１０からの到来方向特徴量ベクトルｚtfと、パラメータ記憶部２２３からのパラメータΘとを受け取って、クラスタ毎かつ時間周波数成分毎の事後確率γtf(n)を更新する（ステップＳ３２）。
【０１１９】
次に、混合重み更新部２２２−１が、事後確率更新部２２１から事後確率γtf(n)を受け取って、クラスタ毎かつ周波数ビン毎の混合重みを更新する（ステップＳ３３）。次に、特徴量共分散行列更新部２２２−２が、事後確率更新部２２１からの事後確率γtf(n)と、特徴抽出部２１０からの到来方向特徴量ベクトルｚtfとを受け取って、目的信号に対応するクラスタの周波数ビン毎の特徴量共分散行列Ｒf(1)を算出する（ステップＳ３４）。
【０１２０】
次に、平均方向ベクトル更新部２２２−３が、特徴量共分散行列更新部２２２−２から特徴量共分散行列Ｒf(1)を受け取って、目的信号に対応するクラスタの周波数ビン毎の平均方向ベクトルａf(1)を更新する（ステップＳ３５）。次に、集中パラメータ更新部２２２−４が、特徴量共分散行列更新部２２２−２から特徴量共分散行列Ｒf(1)を受け取って、目的信号に対応するクラスタの周波数ビン毎の集中パラメータκf(1)を更新する（ステップＳ３６）。次に、図示しない終了条件判定部が、終了条件が満たされているかを判定し（ステップＳ３７）、満たされている場合（ステップＳ３７：Ｙｅｓ）は事後確率更新部２２１が事後確率γtf(n)を出力し、満たされていない場合（ステップＳ３７：Ｎｏ）は、ステップＳ３２に戻って処理を続ける。なお、終了条件は、例えば、ステップＳ３２〜Ｓ３７のループ処理を規定回数繰り返した、事後確率γtf(n)が所定条件を満たした等である。
【０１２１】
以上の各実施形態は、任意のクラスタリング技術と組み合わせて利用することができる。
【０１２２】
また、クラスタ判定装置１００，１００Ｄ、クラスタリング装置２００において行われる各処理は、全部又は任意の一部が、ＣＰＵ等の処理装置及び処理装置により解析実行されるプログラムにて実現されてもよい。また、クラスタ判定装置１００，１００Ｄ、クラスタリング装置２００において行われる各処理は、ワイヤードロジックによるハードウェアとして実現されてもよい。
【０１２３】
また、実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的に行うこともできる。もしくは、実施形態において説明した各処理のうち、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。この他、上述及び図示の処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて適宜変更することができる。
【０１２４】
（プログラムについて）
図１０は、プログラムが実行されることにより、実施形態に係るクラスタ判定装置及びクラスタリング装置が実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ（Central  Processing  Unit）１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。コンピュータ１０００において、これらの各部はバス１０８０によって接続される。
【０１２５】
メモリ１０１０は、ＲＯＭ（Read  Only  Memory）１０１１及びＲＡＭ（Random  Access  Memory）１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ（Basic  Input  Output  System）等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０３１に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１０４１に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１０４１に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１０５１、キーボード１０５２に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１０６１に接続される。
【０１２６】
ハードディスクドライブ１０３１は、例えば、ＯＳ１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわち、クラスタ判定装置１００，１００Ｄ、クラスタリング装置２００の各処理を規定するプログラムは、コンピュータ１０００によって実行される指令が記述されたプログラムモジュール１０９３として、例えばハードディスクドライブ１０３１に記憶される。例えば、クラスタ判定装置１００，１００Ｄ、クラスタリング装置２００における機能構成と同様の情報処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０３１に記憶される。
【０１２７】
また、実施形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０３１に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０３１に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。
【０１２８】
なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０３１に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１０４１等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワーク（ＬＡＮ（Local  Area  Network）、ＷＡＮ（Wide  Area  Network）等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。
【０１２９】
上記の各実施形態は、本願が開示する技術に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【０１３０】
１０Ａ、１０Ｂ、１０Ｃ、１０Ｄ  到来方向局在度算出部
１０Ａ−１  到来方向特徴量算出部
１０Ａ−２  特徴量局在度算出部
１０Ｂ−１、１０Ｃ−１  共分散行列算出部
１０Ｂ−２  固有値分布算出部
１０Ｂ−３  固有値スパース度算出部
１０Ｃ−２  直交分解部
１０Ｃ−３  異方度算出部
２０、２０Ｄ  判定部
２０−１  周波数方向集約部
２０−２  クラスタ置換部
１００、１００Ｄ  クラスタ判定装置
２００  クラスタリング装置
２１０  特徴抽出部
２２０  クラスタリング部
２２１  事後確率更新部
２２２  パラメータ更新部
２２２−１  混合重み更新部
２２２−２  特徴量共分散行列更新部
２２２−３  平均方向ベクトル更新部
２２２−４  集中パラメータ更新部
２２３  パラメータ記憶部
１０００  コンピュータ
１０１０  メモリ
１０２０  ＣＰＵ
【図１】
2018189994.tif 000035
【図２】
2018189994.tif 000036
【図３】
2018189994.tif 000037
【図４】
2018189994.tif 000038
【図５】
2018189994.tif 000039
【図６】
2018189994.tif 000040
【図７】
2018189994.tif 000041
【図８】
2018189994.tif 000042
【図９】
2018189994.tif 000043
【図１０】
2018189994.tif 000044