(57)【要約】【課題】反射や残響の影響が大きい場合でも、信号源数を正確に推定する。【解決手段】周波数領域変換部１０は、観測信号を時間周波数領域に変換し、周波数ビンごとに複数の分離信号に分離する。アクティビティ計算部１２は、分離信号からアクティビティを求める。クラスタリング部１４は、アクティビティを複数のアクティビティ・クラスタにクラスタリングする。信号源数決定部１６は、アクティビティ・クラスタそれぞれに含まれるアクティビティの数に基づいて信号源数の推定値を求める。【選択図】図３
2015055843.tif 000002
【特許請求の範囲】
【請求項１】
観測信号を周波数ビンごとに複数の分離信号に分離したときの、前記分離信号の各々から取得した分離信号のアクティビティを複数のアクティビティ・クラスタにクラスタリングするクラスタリング部と、
前記アクティビティ・クラスタのうち信号源に対応すると判定されるアクティビティ・クラスタの個数を信号源数の推定値として求める信号源数決定部と、
を含む信号源数推定装置。
【請求項２】
請求項１に記載の信号源数推定装置であって、
前記信号源数決定部は、前記アクティビティ・クラスタごとのアクティビティの大きさに対応する特徴量に基づいて前記アクティビティ・クラスタをクラスタリングして得られる二つのクラスタのうち前記アクティビティの大きさに対応する特徴量が大きい方のクラスタに含まれるアクティビティ・クラスタの総数を前記信号源数の推定値とするものである
信号源数推定装置。
【請求項３】
請求項１に記載の信号源数推定装置であって、
前記信号源数決定部は、前記アクティビティ・クラスタごとのアクティビティの大きさに対応する特徴量が所定の閾値以上であるアクティビティ・クラスタの総数を前記信号源数の推定値とするものである
信号源数推定装置。
【請求項４】
請求項１から３のいずれかに記載の信号源数推定装置であって、
前記アクティビティは、時間区間ごとの前記分離信号の寄与度である
信号源数推定装置。
【請求項５】
請求項１から３のいずれかに記載の信号源数推定装置であって、
前記アクティビティは、前記分離信号のパワーを用いて求めた時間周波数マスクである
信号源数推定装置。
【請求項６】
請求項１から５のいずれかに記載の信号源数推定装置であって、
前記クラスタリング部は、前記アクティビティを、アクティビティの分布を表す確率モデルに当てはめ、前記確率モデルを評価する評価関数を用いて前記確率モデルのモデルパラメータを計算することにより、複数のアクティビティ・クラスタにクラスタリングするものである
信号源数推定装置。
【請求項７】
クラスタリング部が、観測信号を周波数ビンごとに複数の分離信号に分離したときの、前記分離信号の各々から取得した分離信号のアクティビティを複数のアクティビティ・クラスタにクラスタリングするクラスタリングステップと、
信号源数決定部が、前記アクティビティ・クラスタのうち信号源に対応すると判定されるアクティビティ・クラスタの個数を信号源数の推定値として求める信号源数決定ステップと、
を含む信号源数推定方法。
【請求項８】
請求項１から６に記載の信号源数推定装置としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、観測信号を用いて信号源の数を推定する技術に関する。
【背景技術】
【０００２】
残響環境において、N個の信号源からの信号が混合されてM個のセンサで観測されたとする（Nは1以上の整数）。このとき、j（1≦j≦M）番目のセンサで観測される信号の時間周波数変換xj(τ,f)は、式（１）でモデル化される。
【０００３】
【数１】
2015055843.tif 000003
【０００４】
ここで、τは時間フレーム番号を表し、fは周波数ビン番号を表し、sk(τ,f)はk番目の信号源の時間周波数変換を表す。また、hjk(f)はk番目の信号源からj番目のセンサへの伝達関数を表し、信号源とセンサの距離による信号の遅延や減衰及び壁などにおける反射や残響の影響をモデル化している。
【０００５】
この発明で扱う信号源数推定の問題は、M個のセンサによる観測信号xj(τ,f)（j=1,…,M）を用いて、信号源の個数Nを推定する問題である。特に、クラスタリングに基づく信号源数推定の枠組みでは、信号源ごとに固有の値を取る特徴量fζ（ζ=1,…,L）を用い、その特徴量のヒストグラムが各信号源に対応するピークを持つことに基づき、クラスタリングとモデル選択により各ピークに一つの確率モデルを当てはめ、確率モデルの数を数えることにより信号源数を推定する。
【０００６】
従来の信号源数推定技術として、音の到来方向のクラスタリングによる音源数推定技術が知られている（例えば、特許文献１）。図１に、特許文献１に記載されたクラスタリングによる音源数推定装置の一例を示す。これは、上述のクラスタリングに基づく信号源数推定の一般的な枠組みにおいて、特徴量として到来方向を用いて、モデル選択機構として混合重みの事前分布にディリクレ分布を与えた場合に相当する。
【０００７】
図１に示す音源数推定装置４００は、J（J≧2）個の収音手段２０ｊ（例えば、マイクロフォン）（j=1,…,J）、周波数領域変換部３０、パワー推定部３２、到来方向推定部３４、方向情報分布推定装置２００、音源数推定部３６からなる。J個の収音手段２０ｊにより観測された観測信号は、周波数領域変換部３０により周波数領域へ変換される。到来方向推定部３４は、周波数領域に変換された観測信号を用いて到来方向を計算する。パワー推定部３２は、クラスタリングの際に尤度関数に与える重みであるパワーを計算する。この重みは、パワーが小さく信頼性の小さいデータがクラスタリングへ与える影響を小さくするために与えられる。方向情報分布推定装置２００は、パワーと到来方向とを用いて、到来方向を各音源に対応するクラスタにクラスタリングする。音源数推定部３６は、クラスタの個数を数えることにより音源数の推定値を計算する。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開２０１０−１４５８３６号公報
【発明の概要】
【発明が解決しようとする課題】
【０００９】
観測信号を用いて計算した音の到来方向の推定値は、壁などにおける反射や残響の影響を受けやすいという問題がある。無響室のように壁などにおける反射や残響の影響が小さく、音の伝搬が平面波伝搬と近似できる場合には、音の到来方向を正確に推定することができるが、壁における反射や残響が大きい場合には正確に推定することが困難である。そのため、特許文献１の音源数推定装置では、このような場合には音源数を正確に推定することが困難であった。
【００１０】
この発明の目的は、壁などにおける反射や残響の影響が大きい場合でも、信号源数を正確に推定することである。
【課題を解決するための手段】
【００１１】
上記の課題を解決するために、この発明の信号源数推定装置は、クラスタリング部と信号源数決定部を含む。クラスタリング部は、観測信号を周波数ビンごとに複数の分離信号に分離したときの、前記分離信号の各々から取得した分離信号のアクティビティを複数のアクティビティ・クラスタにクラスタリングする。信号源数決定部は、アクティビティ・クラスタのうち信号源に対応すると判定されるアクティビティ・クラスタの個数を信号源数の推定値として求める。
【発明の効果】
【００１２】
この発明の信号源数推定技術は、クラスタリングに基づく信号源数推定の枠組みにおいて、特徴量としてアクティビティを用いて信号源数を推定する。アクティビティは、壁などにおける反射や残響の有無にかかわらず、信号ごとに固有のパターンを示す特徴量である。したがって、この発明の信号源数推定技術によれば、壁などにおける反射や残響の影響が大きい場合でも、信号源数を正確に推定することができる。
【図面の簡単な説明】
【００１３】
【図１】図１は、従来の音源数推定装置の機能構成を例示する図である。
【図２】図２は、発明の原理を説明するための図である。
【図３】図３は、信号源数推定装置の機能構成を例示する図である。
【図４】図４は、信号源数推定方法の処理フローを例示する図である。
【図５】図５は、クラスタリング部の機能構成を例示する図である。
【図６】図６は、信号源数決定部の機能構成を例示する図である。
【図７】図７は、実験環境を説明するための図である。
【発明を実施するための形態】
【００１４】
以下、この発明の実施の形態について詳細に説明する。なお、図面中において同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【００１５】
[発明の原理]
実施形態の説明に先立ち、この発明の原理について説明する。
【００１６】
この発明では、従来のクラスタリングにおける信号源数推定の枠組みにおいて、特徴量としてアクティビティを用い、アクティビティをクラスタリングすることで信号源数を推定する。信号のアクティビティとは、分離信号の振幅の時間変化に対応する時系列情報であり、例えば、時間区間ごとの分離信号の振幅値の時系列、時間区間ごとの分離信号の寄与度の時系列などである。信号のアクティビティは、壁などにおける反射や残響の有無にかかわらず、信号ごとに固有のパターンを示すという性質を持つ。この性質に基づいて信号のアクティビティをクラスタリングすることで、壁などにおける反射や残響の影響が大きい場合でも、信号源数を正確に推定することが可能となる。
【００１７】
信号のアクティビティを得るために、従来から知られている音源分離方法を応用する。ただし、従来の音源分離方法は、音源数が既知であることが前提とされるが、この発明では信号源数が未知の状況である。そのため、想定される信号源数よりも多い数を仮定して、従来の音源分離方法により分離信号を得る。
【００１８】
こうして得られた分離信号には、信号源に対応しない雑音などの信号が含まれる。そのため、信号源に対応しない分離信号を信号源に対応するものとして扱うことで信号源数を過大推定してしまうことを防ぐ必要がある。これを実現するためには、信号源に対応しない分離信号を除去する、もしくは、信号源に対応しない分離信号が信号源数推定に与える影響を小さくするなどの対策が考えられる。
【００１９】
そこで、この発明では、信号源に対応しない分離信号のアクティビティが、信号源に対応する分離信号のアクティビティと比較して小さな値をとることが多いという性質を利用して、信号源に対応しない分離信号の影響を排除する。すなわち、分離信号のアクティビティをクラスタリングする際に、小さいアクティビティをクラスタリングに利用しない、あるいは、小さいアクティビティからなるクラスタを排除する。
【００２０】
同一の信号源に由来する分離信号は、そのアクティビティが互いに類似するため、同じクラスタに分類されることが期待される。また、信号源に対応しない分離信号のアクティビティは小さい値を持つので、クラスタリングに用いられない、または、クラスタリング結果から排除される。したがって、最終的にアクティビティをクラスタリングして得られるクラスタの数が信号源の数と高い精度で一致することが期待される。なお、後述の実施形態ではアクティビティをクラスタリングして得られるクラスタを、アクティビティ・クラスタと呼ぶ。
【００２１】
図２に上述した発明の原理を概念的に示す。図中左側の「Speech activity sequences」は観測信号を周波数ビンごとに分離した分離信号から求めたアクティビティである。図中右側の「Speech activity clusters」はアクティビティをクラスタリングしたアクティビティ・クラスタである。詳細は後述するが、アクティビティの時系列をクラスタリングし、アクティビティの代表値ρiが大きいアクティビティ・クラスタ（図中では「Active」と表す）は信号源に対応するものと判定する。一方、アクティビティの代表値ρiが小さいアクティビティ・クラスタ（図中では「Garbage」と表す）は信号源に対応しないものと判定する。そして、信号源に対応するものと判定された上記アクティビティ・クラスタの数を信号源数として決定する。
【００２２】
[発明のポイント]
この発明の第一のポイントは、従来のクラスタリングによる信号源数推定の枠組みにおいて、特徴量としてアクティビティを用い、アクティビティをクラスタリングする点である。
【００２３】
第二のポイントは、信号源数が事前に分かってない状況で、観測信号からアクティビティを求めるために、実際よりも多い信号源数を仮定して信号源分離を行うことにより、アクティビティを求めることである。
【００２４】
第三のポイントは、信号に対応しないアクティビティを除去するか、これによる影響を小さくする手段を導入したことである。これは、実際よりも多い信号源数を仮定したことにより信号に対応しないアクティビティが生じ、信号に対応しないアクティビティを信号に対応するアクティビティと誤認識することで信号源数を過大に推定してしまうことを防ぐためである。
【００２５】
［第一実施形態］
図３を参照して、第一実施形態に係る信号源数推定装置１００の機能構成の一例を説明する。信号源数推定装置１００は、周波数領域変換部１０、アクティビティ計算部１２、クラスタリング部１４、信号源数決定部１６を含む。信号源数推定装置１００は、例えば、中央演算処理装置（Central Processing Unit、CPU）、主記憶装置（Random Access Memory、RAM）などを有する公知又は専用のコンピュータに特別なプログラムが読み込まれて構成された特別な装置である。信号源数推定装置１００は、例えば、中央演算処理装置の制御のもとで各処理を実行する。信号源数推定装置１００に入力されたデータや各処理で得られたデータは、例えば、主記憶装置に格納され、主記憶装置に格納されたデータは必要に応じて読み出されて他の処理に利用される。
【００２６】
図４を参照して、信号源数推定装置１００の実行する信号源数推定方法の処理フローの一例を、実際に行われる手続きの順に従って説明する。
【００２７】
ステップＳ１において、周波数領域変換部１０へM個のマイクロフォン２０1,…,２０Mにより観測された時間領域の混合信号~x(t)が入力される。時間領域の混合信号~x(t)は式（２）で定義される。
【００２８】
【数２】
2015055843.tif 000004
【００２９】
ここで、~xm(t)はm番目のマイクロフォン２０mで観測された時間領域の混合信号を表し、tは時間インデックスを表し、上付き文字のTはベクトルの転置を表す。
【００３０】
周波数領域変換部１０は、入力された時間領域の混合信号~x(t)に対して、短時間フーリエ変換などの時間周波数変換を行い、時間周波数領域の混合信号x(τ,f)を生成し、出力する。出力した時間周波数領域の混合信号x(τ,f)はアクティビティ計算部１２へ入力される。時間周波数領域の混合信号x(τ,f)は式（３）で定義される。
【００３１】
【数３】
2015055843.tif 000005
【００３２】
ここで、xm(τ,f)はm番目のマイクロフォン２０mで観測された時間領域の混合信号~xm(t)の時間周波数変換を表し、τは時間フレーム番号を表し、fは周波数ビン番号を表す。
【００３３】
ステップＳ２において、アクティビティ計算部１２は、入力された時間周波数領域の混合信号x(τ,f)に対し、アクティビティy(ζ)を計算し出力する。出力したアクティビティy(ζ)はクラスタリング部１４へ入力される。
【００３４】
第一実施形態では、分離信号のアクティビティとして、時間区間ごとの分離信号の寄与度（事後確率）の時系列、すなわち、時間周波数領域の混合信号x(τ,f)にk（1≦k≦K）番目のクラスタが寄与する事後確率ηkτfを用いる。事後確率ηkτfを求める方法は、例えば、「特開２００９―０５３３４９号公報（参考文献１）」に開示されている方法が利用できる。
【００３５】
クラスタの総数Kは想定される信号源数の上限以上に設定する。各周波数ビンfにおいて各クラスタkに対して事後確率の時系列[ηk1f…ηkTf]Tが定義できる。ここで、下付き文字に現れるTはフレームの総数、上付き文字のTはベクトルの転置を表す。これらにζ=k+(f-1)Kで定義される通し番号を振って、アクティビティy(ζ)を式（４）により定義する。
【００３６】
【数４】
2015055843.tif 000006
【００３７】
アクティビティの総数はL:=K×Fとなる。ここで、Fは周波数ビンの総数である。
【００３８】
なお、事後確率は上述の例のようにフレームごとに求めてもよいし、数フレームごとに求めたものを用いてもよい。
【００３９】
ステップＳ３からＳ９において、クラスタリング部１４は、入力されたアクティビティy(ζ)をアクティビティの分布を表す確率モデルに当てはめ、確率モデルを評価する所定の評価関数を用いて確率モデルのモデルパラメータを計算することにより、クラスタリングを行う。
【００４０】
信号の振幅の立ち上がり及び立下がりは、信号源ごとに固有のパターンを示すため、同じ信号源に対応する分離信号のアクティビティは互いに類似する。この発明では、この性質を利用して、アクティビティをクラスタリングする。同じ信号源に由来する分離信号のアクティビティを一つのクラスタにまとめることにより、アクティビティのクラスタ（以下、アクティビティ・クラスタとも言う）の個数に基づいて信号源数を推定することができる。
【００４１】
アクティビティの分布を表す確率モデルは、例えば、以下のようにモデル化されたものである。i番目のクラスタに関するアクティビティの分布を、例えば、式（５）に示すようにワトソン（Watson）分布でモデル化する。
【００４２】
【数５】
2015055843.tif 000007
【００４３】
ここで、Γはガンマ関数である。また、z(ζ)は、前処理として、アクティビティy(ζ)をノルムが1となるように式（６）に示すように正規化したものである。
【００４４】
【数６】
2015055843.tif 000008
【００４５】
ここで、||・||はユークリッドノルムを表す。wiはクラスタiに関する正規化後のアクティビティz(ζ)の分布の中心を表し、平均方向（mean orientation）と呼ばれ、||wi||=1なる制約条件を満たす。кiはクラスタiに関する正規化後のアクティビティz(ζ)の分布の集中度合いを表し、密度パラメータ（concentration parameter）と呼ばれる。M(a,b,x)はクンマー（Kummer）関数である。クンマー関数についての詳細は、例えば、「S. Sra and D. Karp, “The multivariate Watson distribution: Maximum-likelihood estimation and other aspects”, Journal of Multivariate Analysis, vol. 114, pp. 256-269, Feb. 2013.（参考文献２）」を参照されたい。
【００４６】
以上より、正規化後のアクティビティz(ζ)の尤度関数は、式（７）で表される混合モデルで与えられる。
【００４７】
【数７】
2015055843.tif 000009
【００４８】
ここで、αiは混合重みであり、αiはΣi=1Kαi=1なる制約条件を満たす。また、Θは式（８）に示すパラメータ集合である。
【００４９】
【数８】
2015055843.tif 000010
【００５０】
ここで、例えば{αi}iのような記法は、{αi|∀i}なる集合を表すものとする。なお、iはアクティビティ・クラスタのインデックスを表しており、時間周波数成分のクラスタのインデックスであるkと区別されたい。
【００５１】
第一実施形態においては、簡単のため、アクティビティ・クラスタの数を、時間周波数成分のクラスタ数と同数のKと仮定しているが、これらは必ずしも同数である必要はない。アクティビティ・クラスタの数は、存在しうる信号源の上限以上となっていれば、自由に設定することができる。
【００５２】
すでに述べたように、アクティビティは信号源ごとに固有のパターンを示すため、アクティビティの分布は、各信号源に対応する複数のピークをもつ。アクティビティ・クラスタの個数に基づいて信号源数を推定するためには、アクティビティの分布の各ピークに一つの確率分布モデルをフィッティングさせることが重要である。そこで、混合重みに、事前分布として式（９）で与えられるディリクレ（Dirichlet）事前分布を与える。
【００５３】
【数９】
2015055843.tif 000011
【００５４】
ここで、Γはガンマ関数であり、φはハイパーパラメータである。φを0＜φ＜1の範囲に定めることにより、少数の混合重みのみが大きな値を取りやすいようにすることができる。混合重み以外のパラメータについては一様な事前分布を仮定する。したがって、パラメータ集合Θの事前分布は、式（１０）で表される。
【００５５】
【数１０】
2015055843.tif 000012
【００５６】
ここで、記号∝は、比例を表す。
【００５７】
クラスタリング部１４では、正規化後のアクティビティz(ζ)を、以上のようにモデル化された確率モデルに当てはめ、確率モデルを評価する所定の評価関数を用いて、事後確率及び最適なパラメータ集合Θを求め、その結果に基づいて、アクティビティをクラスタリングする。
【００５８】
以下、クラスタリング部１４の処理をより具体的に説明する。
【００５９】
クラスタリング部１４は、例えば、図５に示すとおり、正規化部１４０、事後確率計算部１４２、パラメータ更新部１４４、パラメータ保持部１４６、クラスタ決定部１４８から構成される。パラメータ更新部１４４は、さらに、混合重み更新手段１４４０、相関行列更新手段１４４２、平均方向更新手段１４４４、密度パラメータ更新手段１４４６、重み更新手段１４４８から構成される。
【００６０】
図４に示すステップＳ３において、クラスタリング部１４は、混合重みαi、平均方向wi、密度パラメータкiからなるパラメータ集合Θの初期値をパラメータ保持部１４６に記憶する。この初期値は、例えば、式（１１）に示すように設定すればよい。
【００６１】
【数１１】
2015055843.tif 000013
【００６２】
また、wiは{z(ζ)}ζから無作為に選べばよい。
【００６３】
ステップＳ４において、正規化部１４０は、クラスタリングに先立つ前処理として、入力されたアクティビティy(ζ)を、式（１２）に示すようにユークリッドノルムが1となるように正規化し、その計算結果を正規化後のアクティビティz(ζ)として出力する。
【００６４】
【数１２】
2015055843.tif 000014
【００６５】
ステップＳ５において、事後確率計算部１４２は、正規化部１４０が出力した正規化後のアクティビティz(ζ)と、パラメータ保持部１４６に記憶されたパラメータ集合Θ’:={{α’i}i,{w’i}i,{к’i}i}とを用いて、正規化後のアクティビティz(ζ)に関する事後確率γi(ζ)を式（１３）により計算し、出力する。ここで、事後確率γi(ζ):=p(i|z(ζ),Θ’)は、正規化後のアクティビティz(ζ)がi番目のクラスタに属する事後確率である。
【００６６】
【数１３】
2015055843.tif 000015
【００６７】
ステップＳ６において、パラメータ更新部１４４の混合重み更新手段１４４０は、式（１４）に示すように、事後確率計算部１４２の出力する事後確率γi(ζ)とパラメータ保持部１４６に記憶された重みλ’とを用いて、混合重みαiを更新する。
【００６８】
【数１４】
2015055843.tif 000016
【００６９】
ステップＳ６において、パラメータ更新部１４４の相関行列更新手段１４４２は、事後確率計算部１４２の出力する事後確率γi(ζ)と正規化部１４０の出力するアクティビティz(ζ)とを用いて、クラスタiに関する相関行列Riを計算する。相関行列Riは式（１５）により定義される。
【００７０】
【数１５】
2015055843.tif 000017
【００７１】
ステップＳ６において、パラメータ更新部１４４の平均方向更新手段１４４４は、相関行列Riの最大固有値に対応する固有ベクトルであって、ユークリッドノルムが1であるものとして、平均方向wiを更新する。
【００７２】
ステップＳ６において、パラメータ更新部１４４の密度パラメータ更新手段１４４６は、式（１６）に示すように、平均方向wiと相関行列Riとを用いて、密度パラメータкiを更新する。
【００７３】
【数１６】
2015055843.tif 000018
【００７４】
ここで、ri:=wiTRiwiである。
【００７５】
ステップＳ６において、パラメータ更新部１４４の重み更新手段１４４８は、式（１７）に示すように、混合重みαi、平均方向wi、密度パラメータкiを用いて、重みλを更新する。
【００７６】
【数１７】
2015055843.tif 000019
【００７７】
以下、パラメータ更新部１４４における各パラメータの更新式の導出根拠を説明する。パラメータの更新はEMアルゴリズムを導出してそれに基づき行う。なお、ワトソン分布のインデックスiは、EMアルゴリズムにおける隠れ変数として扱う。まず、MAP（maximum a posteriori）推定のための評価関数G(Θ)は、式（１８）（１９）により与えられる。
【００７８】
【数１８】
2015055843.tif 000020
【００７９】
ここで、すでに定義したようにp(z(ζ)│i,wi,кi)は式（２０）で与えられ、const.はパラメータ集合Θによらない定数を表す。
【００８０】
【数１９】
2015055843.tif 000021
【００８１】
EMアルゴリズムで用いる評価関数であるQ関数は、式（２１）で与えられる。
【００８２】
【数２０】
2015055843.tif 000022
【００８３】
式（２１）の第二項は、ディリクレ事前分布に由来する項であり、少数の混合重みαiのみが大きな値を取るようにする効果がある。さらに、マイクや信号源の配置などの条件が変わっても、より適するφの値を適応的に決定するために、「S. Medasani and R. Krishnapuram, “Categorization of image databases for efficient retrieval using robust mixture decomposition”, Computer Vision and Image Understanding, vol. 83, pp. 216-235, 2001.（参考文献３）」にならい、式（２１）を変更し、式（２２）に示すようにハイパーパラメータφの代わりに、適応的なハイパーパラメータλを導入する。
【００８４】
【数２１】
2015055843.tif 000023
【００８５】
ハイパーパラメータλは、EMアルゴリズムの反復ごとに、式（２３）により更新される。
【００８６】
【数２２】
2015055843.tif 000024
【００８７】
各パラメータの更新式は、上記の修正されたQ関数~Q(Θ)のパラメータによる偏微分を0と置くことにより導かれる。ただし、制約条件が存在する場合には、ラグランジュ（Lagrange）の未定乗数法により導かれる。
【００８８】
まず、混合重みの更新式は、制約条件Σi=1Kαi=1に留意して、式（２４）に対し、式（２５）の連立方程式を解くことにより得られる。ここで、μはラグランジュの未定乗数である。
【００８９】
【数２３】
2015055843.tif 000025
【００９０】
平均方向の更新式は、制約条件||wi||=1に留意して、式（２６）に対し、式（２７）の連立方程式を解くことにより得られる。
【００９１】
【数２４】
2015055843.tif 000026
【００９２】
密度パラメータの更新式は、まず、Q関数~Q(Θ)を密度パラメータкiで偏微分し、式(28)を得る。
【００９３】
【数２５】
2015055843.tif 000027
【００９４】
ここで、M'(1/2,T/2,кi)は、M(1/2,T/2,кi)のкiによる偏微分であり、Riはすでに定義したように、式（２９）により定義される。また、ri:=wiTRiwiである。
【００９５】
【数２６】
2015055843.tif 000028
【００９６】
「S. Sra and D. Karp, “The multivariate Watson distribution: Maximum-likelihood estimation and other aspects”, Journal of Multivariate Analysis, vol. 114, pp. 256-269, Feb. 2013.（参考文献４）」に記載の近似法により、式（３０）に示すように、κiの近似解が得られる。
【００９７】
【数２７】
2015055843.tif 000029
【００９８】
ステップＳ７において、クラスタリング部１４は、事前に設定したパラメータ更新回数の最大値max_iter（例えば200）に達しているか、または式（３１）に示す収束条件を満たしているかを判定する。
【００９９】
【数２８】
2015055843.tif 000030
【０１００】
ここで、thr_convは収束判定のための閾値であり、例えば10-10とする。
【０１０１】
ステップＳ８において、パラメータ保持部１４６は、パラメータ更新部１４４での更新処理により得られたパラメータ集合Θ及び重みλを記憶する。パラメータ保持部１４６は、事後確率計算部１４２及びパラメータ更新部１４４での次回の処理の際に、保持しているパラメータ集合Θ及び重みλをパラメータ集合Θ'及び重みλ'として提供する。
【０１０２】
クラスタリング部１４は、ステップＳ７においてパラメータ更新回数の最大値に達するか式（３１）に示した収束条件を満たすまで、事後確率計算部１４２及びパラメータ更新部１４４における処理並びにパラメータ保持部１４６への更新データの読み書きを繰り返し実行する。反復終了後の事後確率γi(ζ)は、クラスタ決定部１４８に入力される。
【０１０３】
ステップＳ９において、クラスタ決定部１４８は、正規化する前のアクティビティy(ζ)と反復終了後の事後確率γi(ζ)とを用い、アクティビティ・クラスタCi（1≦i≦K）を式（３２）により決定する。
【０１０４】
【数２９】
2015055843.tif 000031
【０１０５】
上述のクラスタリング部１４の例では、アクティビティの分布の各ピークに一つの確率分布モデルをフィッティングするためのモデル選択の機構を備えたクラスタリング法として、混合重みに適応的なハイパーパラメータを持つディリクレ事前分布を与えた混合モデルを用いていた。ただし、モデル選択の機構を備えたクラスタリング法は、これに限らず様々な方法を用いてもよく、例えば、固定されたハイパーパラメータを持つディリクレ事前分布を用いてもよいし、ノンパラメトリックベイズモデル（例えば、ディリクレ過程混合モデル）、情報量基準（例えば、赤池情報量基準やベイズ情報量基準）、交叉検証（cross validation）、仮説検定などを用いてもよい。
【０１０６】
また、上述のクラスタリング部１４の例では、すべての周波数ビンにおけるアクティビティを用いてクラスタリングを行う例を説明したが、一部の周波数ビンにおけるアクティビティのみを用いてクラスタリングを行ってもよい。用いる周波数ビンは、好ましくは4kHz〜5.56kHzの周波数帯域に含まれる200個の周波数ビンであるが、これに限らず任意の周波数ビンを選択することができる。
【０１０７】
信号源数決定部１６は、クラスタリング部１４で得たアクティビティ・クラスタを用いて信号源数Nを決定し、出力する。具体的には、小さいアクティビティからなるアクティビティ・クラスタ（すなわち、信号源に対応しないと想定されるアクティビティ・クラスタ）を除く残りのアクティビティ・クラスタの数を信号源数Nとして推定する。
【０１０８】
以下、信号源数決定部１６の処理をより具体的に説明する。
【０１０９】
信号源数決定部１６は、例えば、図６に示すとおり、アクティビティ総和計算部１６０、アクティビティ総和クラスタリング部１６２、要素数計数部１６４から構成される。
【０１１０】
図４に示すステップＳ１０において、信号源数決定部１６のアクティビティ総和計算部１６０は、入力されたアクティビティ・クラスタCi（1≦i≦K）を用いて、式（３３）によりアクティビティ総和ρi（1≦i≦K）を計算し、出力する。
【０１１１】
【数３０】
2015055843.tif 000032
【０１１２】
ここで、sum(・)はベクトルの要素の総和を表す。
【０１１３】
ステップＳ１１において、信号源数決定部１６のアクティビティ総和クラスタリング部１６２は、入力されたアクティビティ総和ρiに基づいて、アクティビティ・クラスタCiを二つのクラスタに分類する。クラスタリングは、例えば、「C.M. Bishop, “Pattern Recognition and Machine Learning,” Springer, 2006.（参考文献５）」に記載のk-meansクラスタリングにおいてクラスタ数を2としたものにより行い、二つのクラスタのうちセントロイド値が小さくない方を信号源に対応するクラスタDactiveとして得る。信号源に対応するアクティビティ・クラスタは信号源に対応しないアクティビティ・クラスタよりも大きなアクティビティ総和を持つと考えられる。したがって、信号源に対応するアクティビティ総和と信号源に対応しないアクティビティ総和とは、それぞれがクラスタをなし、セントロイド値が大きいクラスタが信号源に対応すると考えられるためである。
【０１１４】
ステップＳ１２において、信号源数決定部１６の要素数計数部１６４は、入力されたDactiveの要素数を数え、信号源数Nとして出力する。
【０１１５】
上述の信号源数決定部１６の例では、各アクティビティ・クラスタのアクティビティ総和間の距離に基づいて二つのクラスタにクラスタリングする処理を説明した。ただし、アクティビティ総和はアクティビティの大きさに対応した特徴量の一例であり、アクティビティ総和に限らず、他の特徴量を用いても良い。例えば、各アクティビティ・クラスタに含まれるアクティビティのいずれか一つの大きさ（ノルム）や、各クラスタ内のアクティビティの平均値の大きさ（ノルム）などの、アクティビティ・クラスタ内に含まれるアクティビティの代表値を用いて計算した特徴量であってもよい。
【０１１６】
また、アクティビティ・クラスタ間の距離に基づいて二つのクラスタに分類してから信号源に対応するクラスタを得る代わりに、各アクティビティ・クラスタの代表値（または代表値に当該アクティビティ・クラスタの大きさに対応した重みを乗じた値）の大きさが所定の閾値以上であるか否かを判定し、所定の閾値以上であると判定されたアクティビティ・クラスタの数を信号源数Nとして出力する構成としてもよい。
【０１１７】
［変形例］
第一実施形態では、クラスタリング部１４でアクティビティのクラスタリングを行った後に、信号源数決定部１６で信号源に対応しないアクティビティ・クラスタ（つまり、小さいアクティビティからなるアクティビティ・クラスタ）を排除する処理を行ったが、この順序を逆にしてもよい。すなわち、まず、周波数ビンごとの分離信号のアクティビティを、その大きさ（ノルム）に基づいて二つのクラスタにクラスタリングし、そのうちノルムが小さいアクティビティからなる方のクラスタを、信号源に対応しないクラスタとして除外し、ノルムが大きいアクティビティからなる方のクラスタを、信号源に対応するクラスタとして保持し、前記信号源に対応するクラスタに属するアクティビティのみを、クラスタリング部１４の処理によりアクティビティ・クラスタにクラスタリングしてもよい。この場合、信号源に対応しないアクティビティはすでに除外されているため、信号源数決定部１６においては、単に、混合重みαiが所定の閾値以上であるアクティビティ・クラスタの個数を信号源数Nとして出力すればよい。
【０１１８】
また、第一実施形態では、信号源に対応しないアクティビティは信号源に対応するアクティビティよりも小さいという性質を利用して信号源に対応しないアクティビティを除外している。この性質を利用する他の方法によって信号源に対応しないアクティビティを除外してもよいことは言うまでもない。例えば、まず、周波数ビンごとの分離信号のアクティビティのノルム||y(ζ)||を計算し、前記ノルムが所定の閾値未満であるアクティビティy(ζ)を信号源に対応しないアクティビティとして除外し、前記ノルムが所定の閾値以上であるアクティビティy(ζ)を信号源に対応するアクティビティとして保持し、前記信号源に対応するアクティビティについてのみクラスタリング部１４の処理によりアクティビティ・クラスタにクラスタリングしてもよい。この場合も、信号源に対応しないアクティビティはすでに除外されているため、信号源数決定部１６においては、単に、混合重みαiが所定の閾値以上であるアクティビティ・クラスタの個数を信号源数Nとして出力すればよい。
【０１１９】
あるいは、信号源に対応しないアクティビティを陽に除外する代わりに、クラスタリング部１４におけるクラスタリングにおいて、信号に対応しないアクティビティを無視してもよい。すなわち、クラスタリング部１４における評価関数~Q(Θ)に対し、アクティビティの大きさに基づいて式（３４）に示すように重みづけを行う。
【０１２０】
【数３１】
2015055843.tif 000033
【０１２１】
ここで、ν(ζ)はアクティビティの大きさに基づく重みであり、式（３５）により定義される。
【０１２２】
【数３２】
2015055843.tif 000034
【０１２３】
ここで、指数aは所定の正の数である。式（３４）に基づくパラメータ更新式も、第一実施形態と同様の方法で導出することができる。このように、アクティビティの大きさに基づく重みを導入すれば、クラスタリング部１４におけるクラスタリングは、大きいν(ζ)を持つ（すなわち、信号源に対応する）アクティビティの影響を大きく受け、小さいν(ζ)を持つ（すなわち、信号源に対応しない）アクティビティの影響はあまり受けない。したがって、信号源に対応しないアクティビティを陽に除外しなくても、それらをほとんど無視してクラスタリングを行うことができる。この場合、信号源に対応しないアクティビティはクラスタリングにおいてほとんど無視されるため、信号源数決定部１６においては、単に、混合重みαiが所定の閾値以上であるアクティビティ・クラスタの個数を信号源数Nとして出力すればよい。
【０１２４】
［第二実施形態］
第一実施形態に係る信号源数推定装置１００では、アクティビティ計算部１２において、アクティビティとして事後確率の時系列を用いた。しかし、この発明で利用するアクティビティは事後確率の時系列に限定されず、各周波数ビンにおける時間周波数マスクを用いてもよい。なお、各周波数ビンにおける時間周波数マスクは、分離信号のパワーを用いて計算できる。以下では、時間周波数マスクの一例として、ウィーナー（Wiener）マスクを用いる場合について説明するが、どのような時間周波数マスクを用いてもよい。例えば、「Y. Ephraim and D. Malah, “Speech enhancement using a minimum-mean square error short-time spectral amplitude estimator”, IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-32, pp. 1109-1121, Dec. 1984.（参考文献６）」に開示されているMMSE-STSAマスクを用いてもよいし、「Y. Ephraim and D. Malah, “Speech enhancement using a minimum mean-square error log-spectral amplitude estimator”, IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-33, no. 2, pp. 443-445, Apr. 1985.（参考文献７）」に開示されているMMSE-logSTSAマスクを用いてもよい。
【０１２５】
以下、アクティビティとしてウィーナーマスクを用いる場合のアクティビティ計算部の処理について説明する。
【０１２６】
第二実施形態に係るアクティビティ計算部１３は、入力された時間周波数領域の混合信号x(τ,f)に対し、任意の信号源分離技術を適用してK個の分離信号ξkτf∈Cを求める。ここで、kは分離信号の番号、τはフレーム番号、fは周波数ビン番号、Cは複素数全体を表す。この信号源分離技術は、クラスタリングに基づく信号源分離技術でもよいが、他の任意の信号源分離技術を用いてもよい。例えば、独立成分分析に基づく音源分離技術などを適用することができる。ここで、分離信号の個数Kは、想定される信号源数の上限以上に設定するものとする。
【０１２７】
続いて、アクティビティ計算部１３は、分離信号ξkτfを用いて、ウィーナーマスクηkτfを式（３６）により計算する。
【０１２８】
【数３３】
2015055843.tif 000035
【０１２９】
これにより、各周波数ビンfにおいて、各分離信号の番号kに対して、アクティビティ[ηk1f…ηkTf]Tが定義できる。ここで、下付き文字に現れるTはフレームの総数、上付き文字のTはベクトルの転置を表す。これらにζ=k+(f-1)Kで定義される通し番号を振って、アクティビティy(ζ)を式（３７）により定義する。
【０１３０】
【数３４】
2015055843.tif 000036
【０１３１】
［実験結果］
この発明の信号源数推定技術による信号源数の推定精度を評価するために実験を行った。観測信号は、計測した室内インパルス応答を音源信号に畳み込むことにより生成した。
【０１３２】
図７を参照して実験環境における音源とマイクロフォンの配置を説明する。直方体の空間の中央付近に正三角形を描くように三つのマイクロフォンを配置し、それらのマイクロフォンを取り囲む円を描くように四つの音源を配置した。直方体の空間の大きさは、幅4.45メートル、奥行3.55メートル、高さ2.5メートルである。三つのマイクロフォンの描く正三角形は一辺4センチメートルである。四つの音源の描く円の半径は1.2メートルである。四つの音源は図７の紙面下方向を0°として反時計回りに70°、150°、245°、315°の位置に配置されている。三つのマイクロフォンと四つの音源は水平に設置され、その床面からの高さは1.2メートルとした。
【０１３３】
表１にその他の実験条件を示す。
【０１３４】
【表１】
2015055843.tif 000037
【０１３５】
予備実験により、アクティビティ・クラスタリングに用いる周波数ビンの帯域を4kHz〜5.56kHz（200の周波数ビンを含む）に限定した場合に特に信号源数推定精度が良いことが分かったため、この帯域に限定した。
【０１３６】
時間周波数成分のクラスタリングは、「特開２００９−０５３３４９号公報（参考文献８）」に記載の方法により行った。
【０１３７】
実験では、信号源数、マイクロフォン数、残響時間、クラスタ数（K）の組み合わせの各々に対し、音源信号の組み合わせを変えて8回の試行を行ったうち、正しく信号源数が推定された回数の割合を正解率として算出した。
【０１３８】
表２は、信号源数を2個とし、マイクロフォン数を3個とした場合の実験結果である。使用した音源は図７に示した音源１、２であり、使用したマイクロフォンは図７に示したマイクロフォン１、２、３である。
【０１３９】
【表２】
2015055843.tif 000038
【０１４０】
表３は、信号源数を2個とし、マイクロフォン数を2個とした場合の実験結果である。使用した音源は図７に示した音源１、２であり、使用したマイクロフォンは図７に示したマイクロフォン１、２である。
【０１４１】
【表３】
2015055843.tif 000039
【０１４２】
表４は、信号源数を3個とし、マイクロフォン数を3個とした場合の実験結果である。使用した音源は図７に示した音源１、２、３であり、使用したマイクロフォンは図７に示したマイクロフォン１、２、３である。
【０１４３】
【表４】
2015055843.tif 000040
【０１４４】
表５は、信号源数を3個とし、マイクロフォン数を2個とした場合の実験結果である。使用した音源は図７に示した音源１、２、３であり、使用したマイクロフォンは図７に示したマイクロフォン１、２である。
【０１４５】
【表５】
2015055843.tif 000041
【０１４６】
表６は、信号源数を4個とし、マイクロフォン数を3個とした場合の実験結果である。使用した音源は図７に示した音源１、２、３、４であり、使用したマイクロフォンは図７に示したマイクロフォン１、２、３である。
【０１４７】
【表６】
2015055843.tif 000042
【０１４８】
表２から表４に示すように、信号源数≦マイクロフォン数の場合には、クラスタ数（K）や残響時間によらず100％の正解率が得られた。さらに、表５及び表６に示すように、より困難な条件である信号源数＞マイクロフォン数の場合にも、高い正解率が得られた。表５の3音源、2マイクロフォンの場合、K=12のとき、残響時間によらず100％の正解率が得られた。表６の4音源、3マイクロフォンという最も困難な場合でも、K=8,10,12のとき、残響時間によらず75％以上という高い正解率を達成した。
【０１４９】
以上の実験結果より、この発明の信号源数推定技術により十分な効果が得られることが確認された。
【０１５０】
［プログラム、記録媒体］
この発明は上述の実施形態に限定されるものではなく、この発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。上記実施例において説明した各種の処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。
【０１５１】
また、上記実施形態で説明した各装置における各種の処理機能をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記各装置における各種の処理機能がコンピュータ上で実現される。
【０１５２】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。
【０１５３】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【０１５４】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０１５５】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【符号の説明】
【０１５６】
１００  信号源数推定装置
１０  周波数領域変換部
１２、１３  アクティビティ計算部
１４  クラスタリング部
１４０  正規化部
１４２  事後確率計算部
１４４  パラメータ更新部
１４６  パラメータ保持部
１４８  クラスタ決定部
１６  信号源数決定部
１６０  アクティビティ総和計算部
１６２  アクティビティ総和クラスタリング部
１６４  要素数計数部
４００  音源数推定装置
２０  収音手段
３０  周波数領域変換部
３２  パワー推定部
３４  到来方向推定部
３６  音源数推定部
２００  方向情報分布推定装置
【図１】
2015055843.tif 000043
【図２】
2015055843.tif 000044
【図３】
2015055843.tif 000045
【図４】
2015055843.tif 000046
【図５】
2015055843.tif 000047
【図６】
2015055843.tif 000048
【図７】
2015055843.tif 000049