JP6915579B2 - 信号分析装置、信号分析方法および信号分析プログラム - Google Patents

Description

本発明は、信号分析装置、信号分析方法および信号分析プログラムに関する。

Ｎ´個（Ｎ´は０以上の整数）の音源信号が混在する状況において、それぞれ異なる位置でマイクロホンにより取得された複数の観測信号から、個々の音源信号を推定する音源分離技術がある。Ｎ´は真の音源数であり、Ｎは仮定された音源数であるとする。従来技術では、真の音源数Ｎ´が既知である状況を想定し、仮定された音源数をＮ＝Ｎ´と設定する。

T. Higuchi, N. Ito, S. Araki, T. Yoshioka, M. Delcroix, and T. Nakatani, "Online MVDR Beamformer Based on Complex Gaussian Mixture Model With Spatial Prior for Noise Robust ASR", IEEE/ACM Transactions on Audio, Speech, and Language Processing (ASLP), vol. 25, no. 4, pp. 780−793, Apr. 2017.

ここで、図６および図７を用いて、従来の音源分離装置の構成と処理について説明する。図６は、従来の音源分離装置の構成の一例を示す図である。図７は、従来の音源分離処理の処理手順の一例を示すフローチャートである。なお、ベクトル、行列又はスカラーであるＡに対し、“＾Ａ”と記載する場合は「“Ａ”の直上に“＾”が記された記号」と同じであるとする。また、ベクトル、行列又はスカラーであるＡに対し、“~Ａ”と記載する場合は「“Ａ”の直上に“~”が記された記号」と同じであるとする。

図６及び図７に示すように、従来の信号分析装置１Ｐは、観測信号ベクトル作成部１１Ｐ、初期化部（図示しない）、音源存在事後確率更新部１２Ｐ、記憶部１３Ｐ、音源存在事前確率更新部１４Ｐ、空間共分散行列更新部１５Ｐ、パワーパラメータ更新部１６Ｐ、収束判定部（図示しない）および音源信号成分推定部１７Ｐを有する。

観測信号ベクトル作成部１１Ｐは、まず、入力された観測信号ｙ_ｍ（τ）を取得し（ステップＳ４１）、短時間フーリエ変換などにより時間周波数領域の観測信号ｙ_ｍ（ｔ，ｆ）を計算する（ステップＳ４２）。ここで、ｔ＝１，・・・，Ｔはフレームのインデックスであり、ｆ＝１，・・・，Ｆは周波数ビンのインデックスであり、ｍ＝１，・・・，Ｍはマイクロホンのインデックスであり、τはサンプル点のインデックスである。Ｍ個のマイクロホンはそれぞれ異なる位置に配置されているとする。

次に、観測信号ベクトル作成部１１Ｐは、（１）式のように、取得されたＭ個すべての観測信号ｙ_ｍ（ｔ，ｆ）からなるＭ次元縦ベクトルである観測信号ベクトルｙ（ｔ，ｆ）を時間周波数点ごとに作成する（ステップＳ４３）。ここで、上付きのＴは転置を表す。

初期化部は、音源存在事前確率α_ｎ（ｆ）と、空間共分散行列Ｒ_ｎ（ｆ）と、パワーパラメータｖ_ｎ（ｔ，ｆ）と、の推定値の初期値を計算することでこれらのパラメータを初期化する（ステップＳ４４）。ただし、ｎ＝１，・・・，Ｎは音源のインデックスである。例えば、初期化部は、乱数に基づいてこれらの初期値を計算する。

音源存在事後確率更新部１２Ｐは、観測信号ベクトル作成部１１Ｐからの観測信号ベクトルｙ（ｔ，ｆ）と、音源存在事前確率更新部１４Ｐからの音源存在事前確率（ただし例外として、音源存在事後確率更新部１２Ｐにおける最初の処理の際には、初期化部からの音源存在事前確率の初期値）α_ｎ（ｆ）と、空間共分散行列更新部１５Ｐからの空間共分散行列（ただし例外として、音源存在事後確率更新部１２Ｐにおける最初の処理の際には、初期化部からの空間共分散行列の初期値）Ｒ_ｎ（ｆ）と、パワーパラメータ更新部からのパワーパラメータ（ただし例外として、音源存在事後確率更新部１２Ｐにおける最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖ_ｎ（ｔ，ｆ）と、を受け取って、音源存在事後確率λ_ｎ（ｔ，ｆ）を更新する（ステップＳ４５）。

記憶部１３Ｐは、各音源信号ｎおよび各周波数ビンｆに対する空間共分散行列の事前分布のパラメータを記憶する。

音源存在事前確率更新部１４Ｐは、音源存在事後確率更新部１２Ｐからの音源存在事後確率λ_ｎ（ｔ，ｆ）を受け取って、音源存在事前確率α_ｎ（ｆ）を更新する（ステップＳ４６）。

空間共分散行列更新部１５Ｐは、観測信号ベクトル作成部１１Ｐからの観測信号ベクトルｙ（ｔ，ｆ）と、音源存在事後確率更新部１２Ｐからの音源存在事後確率λ_ｎ（ｔ，ｆ）と、記憶部１３Ｐからの事前分布のパラメータと、パワーパラメータ更新部１６Ｐからのパワーパラメータ（ただし例外として、空間共分散行列更新部１５Ｐにおける最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖ_ｎ（ｔ，ｆ）と、を受け取って、空間共分散行列Ｒ_ｎ（ｆ）を更新する（ステップＳ４７）。

パワーパラメータ更新部１６Ｐは、観測信号ベクトル作成部１１Ｐからの観測信号ベクトルｙ（ｔ，ｆ）と、空間共分散行列更新部１５Ｐからの空間共分散行列Ｒ_ｎ（ｆ）と、を受け取って、パワーパラメータｖ_ｎ（ｔ，ｆ）を更新する（ステップＳ４８）。

収束判定部は、収束したかどうかの判定を行う（ステップＳ４９）。収束判定部によって収束していないと判定された場合（ステップＳ４９：Ｎｏ）、音源存在事後確率更新部１２Ｐでの処理（ステップＳ４５）に戻って、処理が継続される。一方、収束判定部によって収束したと判定された場合（ステップＳ４９：Ｙｅｓ）、音源信号成分推定部１７Ｐでの処理に進む。

音源信号成分推定部１７Ｐは、観測信号ベクトル作成部１１Ｐからの観測信号ベクトルｙ（ｔ，ｆ）と音源存在事後確率更新部１２Ｐからの音源存在事後確率λ_ｎ（ｔ，ｆ）とを受け取って、音源信号成分ｘ_ｎ（ｔ，ｆ）の推定値＾ｘ_ｎ（ｔ，ｆ）を計算し、出力する（ステップＳ５０）。

ここで、従来技術の特徴について説明する。観測信号ベクトル作成部１１Ｐにおいて作成された観測信号ベクトルｙ（ｔ，ｆ）は、Ｎ個の音源信号に由来する成分である音源信号成分ｘ_１（ｔ，ｆ），・・・，ｘ_Ｎ（ｔ，ｆ）の和として、（２）式で表される。

従来技術では、各音源信号は、時間周波数領域において、疎な点においてのみ有意なエネルギーを持つという性質（スパース性）を有すると仮定する。例えば、音声はこのスパース性を比較的よく満たすとされている。この仮定の下では、各時間周波数点では、観測信号ベクトルｙ（ｔ，ｆ）は、Ｎ個の音源信号成分ｘ_１（ｔ，ｆ），・・・，ｘ_Ｎ（ｔ，ｆ）のうちの一つだけからなると近似できる（（３）式）。

ここで、ｎ（ｔ，ｆ）は、時間周波数点（ｔ，ｆ）において存在する音源信号のインデックスであり、１以上Ｎ以下の整数の値を取る。

（３）式のモデルの下では、各時間周波数点（ｔ，ｆ）において存在する音源信号のインデックスｎ（ｔ，ｆ）の推定値＾ｎ（ｔ，ｆ）が得られれば、音源分離を実現できる。すなわち、一旦＾ｎ（ｔ，ｆ）が得られれば、次の（４）式のように、ｎ番目の音源信号が存在する時間周波数点以外の音のエネルギーを遮断するかまたは減衰させることにより、ｎ番目の音源信号成分ｘ_ｎ（ｔ，ｆ）の推定値＾ｘ_ｎ（ｔ，ｆ）を得ることができる、すなわち、音源分離が実現できる。

従来技術では、観測信号ベクトルｙ（ｔ，ｆ）の確率分布を次の（５）式の混合複素ガウス分布でモデル化し、観測信号ベクトルｙ（ｔ，ｆ）にこのモデルを当てはめることにより、ｎ（ｔ，ｆ）の推定を実現する。

ここで、ｐ_Ｇは複素ガウス分布を表す（Ｇはガウス（Gauss）の頭文字である）。Ｒ_ｎ（ｆ）は、各音源の空間的特性（音響伝達特性）を表すパラメータである空間共分散行列であり、ｖ_ｎ（ｔ，ｆ）は、各音源のパワースペクトルをモデル化するパラメータであるパワーパラメータである。α_ｎ（ｆ）は、（６）式を満たす混合重みであり、本明細書では音源存在事前確率とも呼ぶ。

また、Θは、すべての未知パラメータをまとめて表したものであり、具体的には、音源存在事前確率α_ｎ（ｆ）、空間共分散行列Ｒ_ｎ（ｆ）、パワーパラメータｖ_ｎ（ｔ，ｆ）からなる。ひとたびパラメータΘが推定できれば、観測信号ベクトルｙ（ｔ，ｆ）が与えられた下での音源インデックスｎ（ｔ，ｆ）の事後確率を、次の（７）式により求めることができる。

これを用いて、次の（８）式のように音源インデックスｎ（ｔ，ｆ）を推定することができる。

この音源インデックスの推定値を用いれば、（４）式に従って、音源分離を実現できる。

このアプローチに基づいて高精度な音源分離を実現するためにはパラメータΘの正確な推定が鍵となる。一般に、与えられる観測信号の長さが長ければ長いほどパラメータΘの正確な推定が容易になり、与えられる観測信号の長さが短ければ短いほどパラメータΘの正確な推定が困難になる。そこで、与えられる観測信号の長さが短くなった場合におけるパラメータΘの推定精度の劣化を防ぐために、パラメータΘに関する事前知識を表す事前分布を適切に定めることが重要である。事前分布を適切に定めれば、与えられる観測信号の長さが短くなった場合でも、パラメータΘに関する事前知識に基づいて、パラメータΘをある程度正確に推定できるため、パラメータΘの推定精度の大幅な低下を防ぐことができる。事前分布はまた、オンライン処理における音源信号が鳴り始めた直後におけるパラメータの推定精度の劣化防止や、パーミュテーション問題の回避のためにも重要である。

ここで、パーミュテーション問題について説明する。観測信号ベクトルｙ（ｔ，ｆ）は、周波数ビンごとに異なる分布に従う。このため、（５）式のような混合モデルを用いた音源インデックスｎ（ｔ，ｆ）の推定（クラスタリング）に基づく音源分離アプローチでは、一般に、各周波数ビン内に限定した音源の分類（クラスタリング）はできても、異なる周波数間で音源の対応をとることはできない。これが、パーミュテーション問題と呼ばれている。

従来技術では、各音源の音源位置が既知であるという仮定の下、各音源信号の空間的特性をモデル化するパラメータである空間共分散行列Ｒ_ｎ（ｆ）の事前分布ｐ（Ｒ_ｎ（ｆ））を設計していた。具体的には、従来技術では、空間共分散行列Ｒ_ｎ（ｆ）の事前分布ｐ（Ｒ_ｎ（ｆ））を、次の（９）式の逆ウィシャート分布によりモデル化する。

ここで、ＩＷは、逆ウィシャート分布を表す（「ＩＷ」は、「Inverse Wishart（逆ウィシャート）」の頭文字である）。~Ψ_ｎ（ｆ）は事前分布ｐ（Ｒ_ｎ（ｆ））の山（モード）の位置をモデル化するスケール行列であり、~ν_ｎ（ｆ）は事前分布ｐ（Ｒ_ｎ（ｆ））の山の広がりをモデル化する自由度である。以下、自由度~ν_ｎ（ｆ）は音源および周波数ビンに依らず一定であると仮定し、単に~νと書く。事前分布ｐ（Ｒ_ｎ（ｆ））のパラメータであるスケール行列~Ψ_ｎ（ｆ）および自由度~νは、パラメータＲ_ｎ（ｆ）をモデル化するパラメータであり、その意味でハイパーパラメータと呼ばれる。

（９）式より、すべての周波数ビンにおける空間共分散行列Ｒ_ｎ（１），・・・，Ｒ_ｎ（Ｆ）の事前分布ｐ（Ｒ_ｎ（１），・・・，Ｒ_ｎ（Ｆ））は、次の（１０）式のようになる。

ここで周波数間の独立性を仮定した。

従来技術では、各音源の音源位置が既知であるという仮定の下、事前分布ｐ（Ｒ_ｎ（ｆ））のハイパーパラメータであるスケール行列~Ψ_ｎ（ｆ）および自由度~νを既知であると仮定していた。これらのハイパーパラメータは、学習用データに基づいて、事前に学習することができる。すなわち、各音源の音源位置が既知の場合には、音源ごとに既知である音源位置から音源信号が到来する場合の観測信号を実測し、これを学習用データとして用いることにより、事前分布ｐ（Ｒ_ｎ（ｆ））のハイパーパラメータであるスケール行列~Ψ_ｎ（ｆ）および自由度~νを事前に学習することができる。

従来技術では、この事前分布に基づき、以下の（１１）式〜（１４）式に示す更新則を交互に繰り返し適用することにより、パラメータΘを推定する。

（１１）式の処理は音源存在事後確率更新部１２Ｐにおいて、（１２）式の処理は音源存在事前確率更新部１４Ｐにおいて、（１３）式の処理は空間共分散行列更新部１５Ｐにおいて、（１４）式の処理はパワーパラメータ更新部１６Ｐにおいて、それぞれ行われる。音源信号成分推定部１７Ｐは、上記の処理により得られた音源存在事後確率更新部１２Ｐからの音源存在事後確率λ_ｎ（ｔ，ｆ）に基づいて、（８）式により音源インデックスの推定値＾ｎ（ｔ，ｆ）を計算し、さらに（４）式により音源信号成分の推定値＾ｘ_ｎ（ｔ，ｆ）を計算する。

しかしながら、従来技術では、各音源信号に対する音源位置が既知であると仮定しており、各音源信号に対する音源位置が未知である場合には適用できなかった。

本発明は、上記に鑑みてなされたものであって、各音源信号に対する音源位置が未知である場合にも、各音源信号の空間的特性をモデル化するパラメータである空間パラメータ（例えば、空間共分散行列）の事前分布に基づいて音源分離などの信号分析を行うことができる信号分析装置、信号分析方法および信号分析プログラムを提供することを目的とする。

上述した課題を解決し、目的を達成するために、本発明の信号分析装置は、Ｎ個（Ｎは２以上の整数）の信号源からの信号の空間的特性をモデル化するパラメータを空間パラメータとする場合、空間パラメータの各信号源に対する事前分布を、空間パラメータのＫ個（Ｋは２以上の整数）の各信号源位置候補に対する事前分布の線型結合である混合分布によりモデル化するときの混合重みであり、信号源ごとの各信号源位置候補から信号が到来する確率である、信号源位置事前確率を推定する推定部を有することを特徴とする。

本発明によれば、各音源信号に対する音源位置が未知である場合にも、空間パラメータの事前分布に基づいて音源分離などの信号分析を行うことができる。

図１は、第１の実施形態に係る信号分析装置の構成の一例を示す図である。 図２は、第１の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。 図３は、第１の実施形態の変形例４に係る信号分析装置の構成の一例を示す図である。 図４は、第１の実施形態の変形例４に係る信号分析処理の処理手順の一例を示すフローチャートである。 図５は、プログラムが実行されることにより、信号分析装置が実現されるコンピュータの一例を示す図である。 図６は、従来の音源分離装置の構成の一例を示す図である。 図７は、従来の音源分離処理の処理手順の一例を示すフローチャートである。

以下に、本願に係る信号分析装置、信号分析方法および信号分析プログラムの実施形態を図面に基づいて詳細に説明する。また、本発明は、以下に説明する実施形態により限定されるものではない。なお、以下では、ベクトル、行列又はスカラーであるＡに対し、“＾Ａ”と記載する場合は「“Ａ”の直上に“＾”が記された記号」と同じであるとする。また、ベクトル、行列又はスカラーであるＡに対し、“~Ａ”と記載する場合は「“Ａ”の直上に“~”が記された記号」と同じであるとする。

［第１の実施形態］
まず、第１の実施形態に係る信号分析装置について説明する。なお、第１の実施形態においては、Ｎ´個（Ｎ´は０以上の整数）の音源信号が混在する状況において、それぞれ異なる位置でマイクロホンにより取得されたＭ個（Ｍは２以上の整数）の観測信号ｙ_ｍ（τ）（ｍ＝１，・・・，Ｍはマイクロホンのインデックス、τはサンプル点のインデックス）が信号分析装置に入力されるものとする。Ｎ´は真の音源数であり、Ｎは仮定された音源数であるとする。第１の実施形態では、真の音源数Ｎ´が既知である状況を想定し、仮定された音源数をＮ＝Ｎ´と設定する。なお、本第１の実施形態における「音源信号」は、目的信号（例えば、音声）であってもよいし、特定の音源位置から到来する雑音である方向性雑音（例えば、テレビから流れる音楽）であってもよい。また、様々な音源位置から到来する雑音である拡散性雑音を、まとめて１つの「音源信号」とみなしてもよい。拡散性雑音の例としては、雑踏やカフェ等における大勢の人々の話し声、駅や空港における足音、空調による雑音などが挙げられる。

図１および図２を用いて、第１の実施形態の構成と処理について説明する。図１は、第１の実施形態に係る信号分析装置の構成の一例を示す図である。図２は、第１の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。第１の実施形態に係る信号分析装置１は、例えば、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、ＣＰＵ（Central Processing Unit）等を含むコンピュータ等に所定のプログラムが読み込まれて、ＣＰＵが所定のプログラムを実行することで実現される。

図１に示すように、信号分析装置１は、観測信号ベクトル作成部１１、初期化部（図示しない）、推定部１０、記憶部１３、パワーパラメータ更新部１８、パーミュテーション解決部（図示しない）、収束判定部（図示しない）、音源信号成分推定部１９を有する。

まず、信号分析装置１の各部の概要について説明する。観測信号ベクトル作成部１１は、まず、入力された観測信号ｙ_ｍ（τ）を取得し（ステップＳ１）、短時間フーリエ変換などにより時間周波数領域の観測信号ｙ_ｍ（ｔ，ｆ）を計算する（ステップＳ２）。ここで、ｔ＝１，・・・，Ｔはフレームのインデックスであり、ｆ＝１，・・・，Ｆは周波数ビンのインデックスである。

次に、観測信号ベクトル作成部１１は、取得されたＭ個すべての観測信号ｙ_ｍ（ｔ，ｆ）からなるＭ次元縦ベクトルである観測信号ベクトルｙ（ｔ，ｆ）、すなわち（１５）式で表される観測信号ベクトルｙ（ｔ，ｆ）、を時間周波数点ごとに作成する（ステップＳ３）。ここで、上付きのＴは転置を表す。

本第１の実施形態では、各音源信号は、Ｋ個の音源位置の候補のいずれかから到来すると仮定し、それらの音源位置候補をインデックス（以下、「音源位置インデックス」）１，・・・，Ｋで表す。例えば、音源が円卓の周りに着席して会話している複数の話者であり、Ｍ個のマイクロホンが円卓の中央の数ｃｍ四方程度の小領域内に置かれており、音源位置として円卓の中央から見たときの音源の方位角のみに注目するとき、０°〜３６０°をＫ等分したＫ個の方位角Δφ，２Δφ，・・・，ＫΔφ（Δφ＝３６０°／Ｋ）を音源位置候補とすることができる。この例に限らず、一般に任意の所定のＫ点を、音源位置候補として指定することができる。また、音源位置候補は、拡散性雑音を表す音源位置候補でもよい。拡散性雑音は、１つの音源位置から到来するのではなく、多数の音源位置から到来する。このような拡散性雑音も「多数の音源位置から到来する」という１つの音源位置候補とみなすことにより、拡散性雑音が存在する状況でも正確な推定が可能になる。

初期化部は、音源存在事前確率α_ｎ（ｆ）と、音源位置事前確率β_ｋｎと、空間共分散行列Ｒ_ｎ（ｆ）と、パワーパラメータｖ_ｎ（ｔ，ｆ）と、の推定値の初期値を計算する（ステップＳ４）。ただし、ｎ＝１，・・・，Ｎは音源のインデックス、ｋ＝１，・・・・，Ｋは音源位置インデックスである。例えば、初期化部は、乱数に基づいてこれらの初期値を計算する。

推定部１０は、音源位置事前確率を推定する。本第１の実施形態では、Ｎ個の音源の位置からの信号の空間的特性をモデル化するパラメータである空間パラメータとして、空間共分散行列を用いる。音源位置事前確率は、空間共分散行列（空間パラメータ）の各音源に対する事前分布を、空間共分散行列（空間パラメータ）のＫ個（Ｋは２以上の整数）の各音源位置候補に対する事前分布の線型結合である混合分布によりモデル化するときの混合重みであり、音源ごとの各音源位置候補から信号が到来する確率である。推定部１０は、音源存在事後確率更新部１２、音源位置事後確率更新部１４、音源存在事前確率更新部１５、音源位置事前確率更新部１６および空間共分散行列更新部１７を有する。

音源存在事後確率更新部１２は、観測信号ベクトル作成部１１からの観測信号ベクトルｙ（ｔ，ｆ）と、音源存在事前確率更新部１５からの音源存在事前確率（ただし例外として、音源存在事後確率更新部１２における最初の処理の際には、初期化部からの音源存在事前確率の初期値）α_ｎ（ｆ）と、空間共分散行列更新部１７からの空間共分散行列（ただし例外として、音源存在事後確率更新部１２における最初の処理の際には、初期化部からの空間共分散行列の初期値）Ｒ_ｎ（ｆ）と、パワーパラメータ更新部１８からのパワーパラメータ（ただし例外として、音源存在事後確率更新部１２における最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖ_ｎ（ｔ，ｆ）と、を受け取って、音源存在事後確率λ_ｎ（ｔ，ｆ）を更新する（ステップＳ５）。

記憶部１３は、各音源位置候補ｋ、各周波数ビンｆに対する空間共分散行列の事前分布のパラメータを記憶する。

音源位置事後確率更新部１４は、記憶部１３からの事前分布のパラメータと、音源位置事前確率更新部１６からの音源位置事前確率（ただし例外として、音源位置事後確率更新部１４における最初の処理の際には、初期化部からの音源位置事前確率の初期値）β_ｋｎと、空間共分散行列更新部１７からの空間共分散行列（ただし例外として、音源位置事後確率更新部１４における最初の処理の際には、初期化部からの空間共分散行列の初期値）Ｒ_ｎ（ｆ）と、を受け取って、音源位置事後確率μ_ｋｎを更新する（ステップＳ６）。

音源存在事前確率更新部１５は、音源存在事後確率更新部１２からの音源存在事後確率λ_ｎ（ｔ，ｆ）を受け取って、音源存在事前確率α_ｎ（ｆ）を更新する（ステップＳ７）。

音源位置事前確率更新部１６は、音源位置事後確率更新部１４からの音源位置事後確率μ_ｋｎを受け取って、音源位置事前確率β_ｋｎを更新する（ステップＳ８）。

空間共分散行列更新部１７は、観測信号ベクトル作成部１１からの観測信号ベクトルｙ（ｔ，ｆ）と、音源存在事後確率更新部１２からの音源存在事後確率λ_ｎ（ｔ，ｆ）と、記憶部１３からの事前分布のパラメータと、音源位置事後確率更新部１４からの音源位置事後確率μ_ｋｎと、パワーパラメータ更新部１８からのパワーパラメータ（ただし例外として、空間共分散行列更新部１７における最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖ_ｎ（ｔ，ｆ）と、を受け取って、空間共分散行列Ｒ_ｎ（ｆ）を更新する（ステップＳ９）。

パワーパラメータ更新部１８は、観測信号ベクトル作成部１１からの観測信号ベクトルｙ（ｔ，ｆ）と、空間共分散行列更新部１７からの空間共分散行列Ｒ_ｎ（ｆ）と、を受け取って、パワーパラメータｖ_ｎ（ｔ，ｆ）を更新する（ステップＳ１０）。

パーミュテーション解決部は、音源存在事前確率更新部１５からの音源存在事前確率α_ｎ（ｆ）と、空間共分散行列更新部１７からの空間共分散行列Ｒ_ｎ（ｆ）と、パワーパラメータ更新部１８からのパワーパラメータｖ_ｎ（ｔ，ｆ）と、を受け取り、音源存在事前確率α_ｎ（ｆ）と、空間共分散行列Ｒ_ｎ（ｆ）と、パワーパラメータｖ_ｎ（ｔ，ｆ）と、を更新することでパーミュテーション問題を解決する（ステップＳ１１）。具体的には、パーミュテーション解決部は、尤度または対数尤度または補助関数などの評価値が最大となるように、音源インデックスｎを周波数ビンごとに付け替えることにより、これらのパラメータを更新する。すなわち、周波数ビンｆにおける音源インデックスｎの付け替えを全単射σ_f：｛１，・・・，Ｎ｝→｛１，・・・，Ｎ｝で表すとき、各周波数ビンｆにおいてこれらのパラメータの音源インデックスｎをσ_f（ｎ）に付け替えたときの尤度または対数尤度または補助関数などの評価値が最大になるように全単射σ_fを求め、求めた全単射σ_fを用いて各周波数ビンｆにおいてこれらのパラメータの音源インデックスｎをσ_f（ｎ）に付け替えることにより、これらのパラメータを更新する。なお、パーミュテーション解決部は、音源存在事前確率α_ｎ（ｆ）と、空間共分散行列Ｒ_ｎ（ｆ）と、パワーパラメータｖ_ｎ（ｔ，ｆ）と、のすべてを更新する代わりに、その一部のみ（例えば、空間共分散行列Ｒ_ｎ（ｆ）のみ）を更新してもよい。なお、パーミュテーション解決部での処理は必須ではない。

続いて、収束判定部は、収束したかどうかの判定を行う（ステップＳ１２）。収束判定部が収束していないと判定した場合（ステップＳ１２：Ｎｏ）、音源存在事後確率更新部１２での処理（ステップＳ５）に戻って、以降の処理が継続される。一方、収束判定部が収束したと判定した場合（ステップＳ１２：Ｙｅｓ）、音源信号成分推定部１９における処理（ステップＳ１３）に進む。

音源信号成分推定部１９は、観測信号ベクトル作成部１１からの観測信号ベクトルｙ（ｔ，ｆ）と音源存在事後確率更新部１２からの音源存在事後確率λ_ｎ（ｔ，ｆ）とを受け取って、音源信号成分ｘ_ｎ（ｔ，ｆ）の推定値＾ｘ_ｎ（ｔ，ｆ）を計算し、出力する（ステップＳ１３）。

次に、第１の実施形態の特徴について、従来技術と対比しながら説明する。前述の通り、従来技術では、すべての周波数ビンにおける空間共分散行列Ｒ_ｎ（１），・・・，Ｒ_ｎ（Ｆ）の事前分布ｐ（Ｒ_ｎ（１），・・・，Ｒ_ｎ（Ｆ））を、次の（１６）式（（１０）式を再掲）によりモデル化していた。

しかしながら、従来技術では、各音源の音源位置が既知であると仮定しており、各音源の音源位置が未知の場合には適用できないという問題があった。

これに対し、本第１の実施形態では、すべての周波数ビンにおける空間共分散行列Ｒ_ｎ（１），・・・，Ｒ_ｎ（Ｆ）の事前分布ｐ（Ｒ_ｎ（１），・・・，Ｒ_ｎ（Ｆ））を、次の（１７）式の混合複素逆ウィシャート分布でモデル化する。

これは、音源位置候補ｋに対する事前分布を、音源ｎが音源位置候補ｋにある確率β_ｋｎを重みとして平均した形になっている。本第１の実施形態では各音源の音源位置が未知であると仮定しているから、β_ｋｎは未知の確率である。ただし、β_ｋｎは確率であるから、次の（１８）式を満たすものとする。

このように、未知の確率β_ｋｎによる重み付き和に基づくことで、各音源の音源位置が未知の場合でも、空間共分散行列の事前分布を設計することができる。β_ｋｎは未知であるが、これも未知パラメータとみなし、他の未知パラメータと同時に推定することができる。

本第１の実施形態では、各音源位置候補ｋ、各周波数ビンｆに対する複素逆ウィシャート分布のパラメータΨ_ｋ（ｆ），ν_ｋ（ｆ）は、事前に準備され、記憶部１３に記憶されているものとする。これらのパラメータは、マイクロホン配置の情報に基づいて事前に準備してもよいし、実測データから事前に学習してもよい。

例えば、マイクロホン配置の情報に基づいて事前に準備する場合には、各マイクロホンｍの直交座標をｒ_ｍとして、各音源位置候補ｋに対応する平面波のステアリングベクトルを（１９）式により計算し、Ψ_ｋ（ｆ），ν_ｋ（ｆ）を次の（２０）式および（２１）式により計算すればよい。

ここで、ｄ_ｋはｋ番目の音源位置候補に対応する音源信号の到来方向を表す単位ベクトル、ｃは音速、ω_ｆは周波数ビンｆに対応する角周波数、（２１−１）式に示すｊは虚数単位、上付きのＨはエルミート転置である。

ここで、本第１の実施形態における事前分布（（１７）式）の導出について説明する。各音源の音源位置は未知であると仮定し、各音源ｎの音源位置に対応する音源位置インデックスｋ_ｎは、（２２）式に示す未知の確率分布に従うと仮定する。β_ｋｎは、音源ごとの音源位置インデックスの確率分布である音源位置事前確率である。

さらに、本第１の実施形態では、音源ｎに対する音源位置インデックスがｋ_ｎ＝ｋであるという条件の下で、音源ｎの空間共分散行列Ｒ_ｎ（１），・・・，Ｒ_ｎ（Ｆ）が、互いに独立に確率分布（（２３）式）に従うものとする。

ここで、Ψ_ｋ（ｆ）は、各音源位置候補に対する空間共分散行列の事前分布の山（モード）の位置を表すパラメータ（スケール行列）であり、ν_ｋ（ｆ）は、各音源位置候補に対する空間共分散行列の事前分布の山の広がり（自由度）を表すパラメータである。また、ＩＷ_Ｃ（Σ；Ψ，ν）は、（２４）式に示すものであり、スケール行列がΨ、自由度がνである複素逆ウィシャート分布である。

（２２）式および（２３）式のモデル化の下では、音源ｎの空間共分散行列Ｒ_ｎ（１），・・・，Ｒ_ｎ（Ｆ）の確率分布は、次の（２５）式〜（２８）式で与えられる。

本実施形態では、事前分布（（１７）式）に基づき、パラメータを推定する。以下、本実施形態におけるパラメータ推定アルゴリズムについて説明する。なお、以下では簡単のため、複素逆ウィシャート分布「ＩＷ_Ｃ」を、添え字Ｃを省略して単に「ＩＷ」と表す。空間共分散行列Ｒ_ｎ（ｆ）以外の未知パラメータの事前分布は一様分布であると仮定すると、パラメータΘの事前分布は次の（２９）式および（３０）式で与えられる。

なお、本第１の実施形態におけるパラメータΘは、音源存在事前確率α_ｎ（ｆ）、パワーパラメータｖ_ｎ（ｔ，ｆ）、空間共分散行列Ｒ_ｎ（ｆ）および音源位置事前確率β_ｋｎからなる。

一方、パラメータΘが与えられた下で、各時間周波数点における観測信号ベクトルｙ（ｔ，ｆ）が互いに独立であると仮定すると、尤度が次の（３１）式および（３２）式で与えられる。

ここで、Ｙは、すべての時間周波数点における観測信号ベクトルｙ（ｔ，ｆ）をまとめて表したものである。

本第１の実施形態では、パラメータΘの事後確率ｐ（Θ｜Ｙ）を最大化することにより、パラメータΘを推定する。ベイズの定理より、この事後確率は（３３）式のように表せ、両辺の対数を取ると、（３４）式となる。

ｌｎｐ（Ｙ）はパラメータΘに依らないから、事後確率ｐ（Θ｜Ｙ）のΘに関する最大化は、次の（３５）式のΘに関する最大化と等価であり、したがって次の（３６）式に示す目的関数Ｊ（Θ）のΘに関する最大化と等価である。

ここで、＝の上に“ｃ”が記された記号は、パラメータΘに依存しない定数の差を除いて両辺が等しいことを表す記号である。また、「Ａ＝：Ｂ」は、ＢをＡによって定義することを表す。

上式の目的関数Ｊ（Θ）の最大化は、補助関数法に基づいて行うことができる。補助関数法では、パラメータΘと補助変数と呼ばれる変数Φとの関数である補助関数Ｑ（Θ，Φ）に基づいて、以下の２つのステップを交互に反復する。
１．補助関数Ｑ（Θ，Φ）を補助変数Φに関して最大化することにより、補助変数Φを更新するステップ
２．補助関数Ｑ（Θ，Φ）が減少しないようにパラメータΘを更新するステップ

ただし、補助関数Ｑ（Θ，Φ）は、次の（３７）式に示す条件を満たすものとする。

この補助関数法により、目的関数Ｊ（Θ）を単調増加させることができる。すなわち、ｉ回目の反復の結果得られたパラメータΘの推定値をΘ^（ｉ）として、（３８）式が成り立つ。

実際、ｉ回目の反復の結果得られた補助変数Φの値をΦ^（ｉ）とすると、（３７）式より、（３９）式および（４０）式が成り立つ。

しかるに、以下の（４１）式が成り立つから、（３８）式が得られる。

補助関数法においては、（３７）式を満たすような補助関数Ｑ（Θ，Φ）を設計する必要がある。そのために、本第１の実施形態では、イェンセンの不等式を用いる。ｆを凸関数とし、ｗ_１，・・・，ｗ_Ｌを（４２）式を満たす非負の数とし、ｘ_１，・・・，ｘ_Ｌを実数とするとき、（４３）式が成り立つ（等号成立条件はｘ_１＝・・・＝ｘ_Ｌ）ことが知られている。

これは、イェンセンの不等式と呼ばれる。特に、ｆ（ｘ）＝−ｌｎｘとすると、（４４）式を得る。

λ_１（ｔ，ｆ），・・・，λ_Ｎ（ｔ，ｆ）を（４５）式を満たす非負の数とすると（４４）式より（４６）式および（４７）式が得られる。

また、μ_１ｎ，・・・，μ_Ｋｎを（４８）式を満たす非負の数とすると、（４４）式より（４９）式および（５０）式が得られる。

（４７）式および（５０）式より、（５１）式が得られる。

よって、（５１）式の右辺を、（５２）式とおくと、（３６）式および（５１）式より、（５３）式が成り立つ。

ただし、補助変数Φは、λ_ｎ（ｔ，ｆ）とμ_ｋｎとからなるものとする。

（５１）式の等号成立条件は、（５４）式および（５５）式である。

これは、次の（５６）式および（５７）式と等価である。

したがって、（５８）式が成り立つ。

（５３）式および（５８）式より、（５２）式のＱ（Θ，Φ）が（３７）式を満たすことが分かる。これで、目的関数Ｊ（Θ）に対する補助関数が設計できた。

本第１の実施形態では、(５２)式の補助関数Ｑ（Θ，Φ）に基づいて、補助変数ΦおよびパラメータΘを次のようにして更新する。まず、補助変数Φの更新は、（５６）式および（５７）式により行えばよい。また、パラメータΘの更新は、次の（５９）式〜（６２）式を用いて行えばよい。

このように、本第１の実施形態では、（３６）式の目的関数を直接最大化する代わりに、補助関数Ｑ（Θ，Φ）に基づいて、補助関数Ｑ（Θ，Φ）を補助変数Φに関して最大化することによりΦを更新するステップと、補助関数Ｑ（Θ，Φ）が減少しないようにパラメータΘを更新するステップと、を交互に反復することにより、（３６）式の目的関数を間接的に最大化する。（３６）式の目的関数においては、対数ｌｎの中にｋに関する和Σ_ｋ＝１ ^Ｋが含まれており、（３６）式の目的関数の各パラメータに関する微分が煩雑な形になるため、（３６）式の目的関数を勾配法などにより直接最大化しようとすると、更新則が煩雑な形になる。これに対し、補助関数Ｑ（Θ，Φ）では、ｋに関する和Σ_ｋ＝１ ^Ｋが対数ｌｎの外に出た形になっており、補助関数Ｑ（Θ，Φ）の各パラメータに関する微分が単純な形になる。また、勾配法では、反復ごとのパラメータの更新量を定めるステップサイズを調整する必要があるが、補助関数法では、ステップサイズが不要であるため、ステップサイズを調整する必要がない。

（５６）式により更新されたλ_ｎ（ｔ，ｆ）は、観測信号ベクトルｙ（ｔ，ｆ）が観測された「後」の音源存在確率に他ならない。実際、ベイズの定理より、（５６）式は（６３）式とも書ける。

そこで、λ_ｎ（ｔ，ｆ）を音源存在事後確率と呼ぶ。これに対し、α_ｎ（ｆ）（（６４）式））は、観測信号ベクトルｙ（ｔ，ｆ）が観測される「前」の音源存在確率であるから、音源存在事前確率と呼ぶ。

また、（５７）式により更新されたμ_ｋｎは、空間共分散行列Ｒ_ｎ（１），・・・，Ｒ_ｎ（Ｆ）が与えられた「後」の音源位置確率に他ならない。実際、（５７）は、（６５）式とも書ける。

そこで、μ_ｋｎを音源位置事後確率と呼ぶ。これに対し、β_ｋｎ（（６６）式）は、空間共分散行列Ｒ_ｎ（１），・・・，Ｒ_ｎ（Ｆ）が与えられる「前」の音源位置確率であるため、音源位置事前確率と呼ぶ。

（５６）式の処理は音源存在事後確率更新部１２において、（５７）式の処理は音源位置事後確率更新部１４において、（５９）式の処理は音源存在事前確率更新部１５において、（６０）式の処理は音源位置事前確率更新部１６において、（６１）式の処理は空間共分散行列更新部１７において、（６２）式の処理はパワーパラメータ更新部１８において、それぞれ行われる。

ここで、上述のパラメータΘの更新則（５９）式〜（６２）式の導出について説明する。まず、（５２）式の補助関数は次の（６７）式および（６８）式のように計算できる。ここで、ＣはパラメータΘに依らない定数である。

音源存在事前確率α_ｎ（ｆ）の更新則（５９）式を導出するために、拘束条件（６）式に注意して、ξをラグランジュの未定乗数として、（６９）式をα_ｎ（ｆ）で微分したものを０と置くと、（７０）式となる。

（７０）式をα_ｎ（ｆ）について解くと、（７１）式となる。

（７１）式に含まれるラグランジュの未定乗数ξの値を決定するために、（７１）式を拘束条件（６）式に代入すると、（７２）式〜（７４）式となる。

よって、ξ＝Ｔであるから、音源存在事前確率α_ｎ（ｆ）の更新則（５９）式が得られる。音源位置事前確率β_ｋｎの更新則（６０）式も同様にして導出できるから、説明を省略する。

空間共分散行列Ｒ_ｎ（ｆ）の更新則（６１）式を導出するために、（６８）式をＲ_ｎ（ｆ）で微分したものを０と置くと、（７５）式となる。

上式の両辺に対し、左右からそれぞれＲ_ｎ（ｆ）を掛けると、（７６）式となる。これをＲ_ｎ（ｆ）について解けば、空間共分散行列Ｒ_ｎ（ｆ）の更新則（６１）式が得られる。

パワーパラメータｖ_ｎ（ｔ，ｆ）の更新則（６２）式を導出するために、（６８）式をｖ_ｎ（ｔ，ｆ）で微分したものを０と置くと、（７７）式となる。

これをｖ_ｎ（ｔ，ｆ）について解けば、パワーパラメータｖ_ｎ（ｔ，ｆ）の更新則（６２）式が得られる。以上で、上述のパラメータΘの更新則（５９）式〜（６２）式が導出できた。

本第１の実施形態においては、複素ガウス分布のパラメータである空間共分散行列Ｒ_ｎ（ｆ）の事前分布が、複素逆ウィシャート分布に基づく事前分布であるというモデル化に基づいている。このように、複素ガウス分布と複素逆ウィシャート分布とを組み合わせて用いることにより、補助関数Ｑ（Θ，Φ）が、その空間共分散行列Ｒ_ｎ（ｆ）に関する微分を０と置いた式がＲ_ｎ（ｆ）について解ける（上述）ような形になる。これは、複素逆ウィシャート分布が複素ガウス分布の共役事前分布であることに起因する。共役事前分布については、参考文献２「C.M. Bishop,“Pattern Recognition and Machine Learning”, Springer, 2006.」を参照されたい。

［第１の実施形態の効果］
このように、本第１の実施形態では、空間共分散行列の各信号源に対する事前分布を、空間共分散行列の複数の各信号源位置候補に対する事前分布の線型結合である混合分布によりモデル化するときの混合重みであり、信号源ごとの各信号源位置候補から信号が到来する確率である、信号源位置事前確率を推定する。具体的には、本第１の実施形態では、空間共分散行列の各信号源に対する事前分布を（１７）式のようにモデル化している。そして、本第１の実施形態では、未知の確率である音源位置事前確率β_ｋｎによる重み付き和に基づくことによって、各音源の音源位置が未知の場合でも、空間共分散行列の事前分布を設計することができる。したがって、本第１の実施形態では、各音源信号に対する音源位置が未知である場合にも、空間共分散行列の事前分布に基づいて音源分離を行うことができる。

また、本第１の実施形態では、（５２）式に示すように、対数ｌｎの中にｋに関する和がない補助関数を用いるため、補助関数の各パラメータに関する微分が単純になり、パラメータの更新演算が煩雑ではなくなる。

また、本第１の実施形態では、空間共分散行列の事前分布が、複素逆ウィシャート分布に基づく事前分布であるというモデル化に基づいている。このように、本第１の実施形態では、複素ガウス分布と複素逆ウィシャート分布とを組み合わせて用いることにより、補助関数Ｑ（Θ，Φ）が、その空間共分散行列Ｒ_ｎ（ｆ）に関する微分を０と置いた式がＲ_ｎ（ｆ）について解ける。

［第１の実施形態の変形例１］
本第１の実施形態では、観測データとして観測信号ベクトルｙ（ｔ，ｆ）を用いたが、他の特徴ベクトルまたは特徴量を観測データとして用いてもよい。例えば、観測信号ベクトルｙ（ｔ，ｆ）に基づいて、（７８）式および（７９）式で定義される特徴ベクトルｚ（ｔ，ｆ）を用いてもよい。

また、観測データとして、マイクロホン間の位相差、振幅比や、音源信号の到来時間差、到来方向などの特徴量を用いてもよい。

また、本第１の実施形態では、特徴ベクトルである観測信号ベクトルに当てはめる混合モデルとして、混合複素ガウス分布を用いたが、利用される特徴ベクトルに応じて、様々な混合モデル（例えば、混合ガウス分布、混合ラプラス分布、混合複素ワトソン分布、混合複素ビンガム分布、混合複素角度中心ガウス分布、フォンミーゼス分布など）を用いることができる。また、混合モデルに限らず、複素ガウス分布などのモデルを特徴ベクトルである観測信号ベクトルに当てはめてもよい。

また、本第１の実施形態では、空間共分散行列の事前分布を混合複素逆ウィシャート分布によりモデル化したが、混合複素ウィシャート分布などの他のモデルによりモデル化してもよい。

また、本第１の実施形態では、モデルを観測データに当てはめるために、パラメータΘの事後確率を最大化する方法を採用したが、他の方法によりモデルを観測データに当てはめてもよい。

また、本第１の実施形態では、補助関数法により最適化を行ったが、勾配法などの他の方法により最適化を行ってもよい。その場合、音源存在事後確率更新部１２および音源位置事後確率更新部１４は必須ではない。

［第１の実施形態の変形例２］
真の音源数Ｎ´が未知の場合に、真の音源数Ｎ´の推定や音源分離を行う第１の実施形態の変形例２について説明する。本変形例では、仮定された音源数ＮはＮ≧Ｎ´となるように十分大きく設定されているものとする。例えば、想定される音源数が高々６個であると分かっているような場合には、仮定された音源数はＮ＝６と設定すればよい。なお、実際の音源数は４個である場合には、Ｎ´＝４となる。

推定部１０は、各ｎ（ｎは１以上Ｎ以下の整数）に対し、音源位置事前確率更新部１６からの音源位置事前確率β_ｋｎが最大となるｋに対応する音源位置候補を音源位置の推定値とする。そして、信号分析装置１は、このようにして得られたＮ個の音源位置を、階層クラスタリングなどによりクラスタリングし、得られたクラスタの個数を、実際の音源数Ｎ´の推定値＾Ｎ´とする。

クラスタリングにより得られた＾Ｎ´個の各クラスタは、＾Ｎ´個の実際の音源に対応するものとみなされる。従ってこのクラスタリングにより、Ｎ個の仮定する各音源ｎが、＾Ｎ´個の実際の音源のうちのどれに対応するか、が分かる。音源分離を行う場合には、この対応関係を利用して、推定部１０が以降の処理も行う。

推定部１０は更に、得られた＾Ｎ´個の各クラスタｎ´（ｎ´は１以上＾Ｎ´以下の整数であるクラスタのインデックス）に対し、Ｎ個の仮定する音源の音源存在事後確率λ_ｎ（ｔ，ｆ）のうち該クラスタに対応するものを加算することにより、ｎ´番目の実際の音源の音源存在事後確率λ´_ｎ´（ｔ，ｆ）を計算する。推定部１０は更に、式（８）と同様に、各時間周波数点（ｔ，ｆ）に対し、実際の音源の音源存在事後確率λ´_ｎ´（ｔ，ｆ）が最大となる番号ｎ´に対応する実際の音源からの信号が（ｔ，ｆ）にて鳴っていると判定する。推定部１０は更に、（４）式と同様に、実際の音源の音源信号成分の推定値＾ｘ´_ｎ´（ｔ，ｆ）を、（ｔ，ｆ）においてｎ´番目の実際の音源が鳴っていると判定された場合にはｙ（ｔ，ｆ）とし、そうでないと判定された場合には０とすることにより、音源分離を行う。

［第１の実施形態の変形例３］
本第１の実施形態は、音信号に限らず、他の信号（脳波、脳磁図、無線信号など）に対して適用してもよい。本第１の実施形態における観測信号は、複数のマイクロホン（マイクロホンアレイ）により取得された観測信号に限らず、脳波計、脳磁計、アンテナアレイなどの他のセンサアレイ（複数のセンサ）により取得された、空間上の位置から時系列として発生する信号からなる観測信号であってもよい。

［第１の実施形態の変形例４］
観測信号ベクトルｙ（ｔ，ｆ）の確率分布を次の（８０）式の複素ガウス分布によりモデル化する例を第１の実施形態の変形例４として説明する。この場合のパラメータΘの更新則は、第１の実施形態の（５６）、（５７）、（５９）、（６０）、（６１）、（６２）式に代えて、（８１）式〜（８６）式のようになる。

図３および図４を用いて、第１の実施形態の変形例４の構成と処理について説明する。図３は、第１の実施形態の変形例４に係る信号分析装置の構成の一例を示す図である。図４は、第１の実施形態の変形例４に係る信号分析処理の処理手順の一例を示すフローチャートである。

図３に示すように、本第１の実施形態の変形例４に係る信号分析装置２０１は、観測信号ベクトル作成部１１、初期化部（図示しない）、記憶部１３、推定部２１０、パワーパラメータ更新部２１８、収束判定部（図示しない）を有する。推定部２１０は、音源位置事後確率更新部２１２、音源信号事後確率更新部２１３、音源位置事前確率更新部２１４、空間共分散行列更新部２１７を有する。

観測信号ベクトル作成部１１は、第１の実施形態と同様に、観測信号ベクトルｙ（ｔ，ｆ）を（１）式により作成する（ステップＳ２１〜ステップＳ２３）。

初期化部は、音源位置事前確率β_ｋｎと、空間共分散行列Ｒ_ｎ（ｆ）と、パワーパラメータｖ_ｎ（ｔ，ｆ）と、の推定値の初期値を計算する（ステップＳ２４）。ただし、ｎ＝１，・・・，Ｎは音源のインデックス、ｋ＝１，・・・，Ｋは音源位置候補のインデックスである。例えば初期化部は、乱数に基づいてこれらの初期値を計算する。また、初期化部は、ｎを初期化する（ステップＳ２５）。

なお、記憶部１３は、各音源位置候補ｋ、各周波数ビンｆに対する空間共分散行列の事前分布のパラメータであるΨ_ｋ（ｆ）およびν_ｋ（ｆ）を記憶する。

続いて、信号分析装置２０１は、ｎに１を加算して（ステップＳ２６）、ステップＳ２７〜ステップＳ３１の処理を行う。

音源位置事後確率更新部２１２は、記憶部１３からの事前分布のパラメータであるΨ_ｋ（ｆ）およびν_ｋ（ｆ）と、音源位置事前確率更新部２１４からの音源位置事前確率（ただし例外として、音源位置事後確率更新部２１２における最初の処理の際には、初期化部からの音源位置事前確率の初期値）β_ｋｎと、空間共分散行列更新部２１７からの空間共分散行列（ただし例外として、音源位置事後確率更新部２１２における最初の処理の際には、初期化部からの空間共分散行列の初期値）Ｒ_ｎ（ｆ）と、を受け取って、音源位置事後確率μ_ｋｎを（８１）式により更新する（ステップＳ２７）。

音源信号事後確率更新部２１３は、観測信号ベクトル作成部１１からの観測信号ベクトルｙ（ｔ，ｆ）と、パワーパラメータ更新部２１８からのパワーパラメータ（ただし例外として、音源信号事後確率更新部２１３における最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖ_ｎ（ｔ，ｆ）と、空間共分散行列更新部２１７からの空間共分散行列（ただし例外として、音源信号事後確率更新部２１３における最初の処理の際には、初期化部からの空間共分散行列の初期値）Ｒ_ｎ（ｆ）と、を受け取って、音源信号成分ｘ_ｎ（ｔ，ｆ）の事後確率の平均ξ_ｎ（ｔ，ｆ）および共分散行列Σ_ｎ（ｔ，ｆ）を、（８２）式および（８３）式により更新する（ステップＳ２８）。

音源位置事前確率更新部２１４は、音源位置事後確率更新部２１２からの音源位置事後確率μ_ｋｎを受け取って、音源位置事前確率β_ｋｎを（８４）式により更新する（ステップＳ２９）。

空間共分散行列更新部２１７は、記憶部１３からの事前分布のパラメータであるΨ_ｋ（ｆ）およびν_ｋ（ｆ）と、音源位置事後確率更新部２１２からの音源位置事後確率μ_ｋｎと、音源信号事後確率更新部２１３からの事後確率の平均ξ_ｎ（ｔ，ｆ）および共分散行列Σ_ｎ（ｔ，ｆ）と、パワーパラメータ更新部２１８からのパワーパラメータ（ただし例外として、空間共分散行列更新部２１７における最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖ_ｎ（ｔ，ｆ）と、を受け取って、空間共分散行列Ｒ_ｎ（ｆ）を（８５）式により更新する（ステップＳ３０）。

パワーパラメータ更新部２１８は、空間共分散行列更新部２１７からの空間共分散行列Ｒ_ｎ（ｆ）と、音源信号事後確率更新部２１３からの事後確率の平均ξ_ｎ（ｔ，ｆ）および共分散行列Σ_ｎ（ｔ，ｆ）と、を受け取って、パワーパラメータｖ_ｎ（ｔ，ｆ）を（８６）式により更新する（ステップＳ３１）。

そして、信号分析装置２０１は、ｎ＝Ｎか否かを判定する（ステップＳ３２）。信号分析装置２０１は、ｎ＝Ｎでないと判定した場合（ステップＳ３２：Ｎｏ）、ステップＳ２６に戻る。これに対し、信号分析装置２０１は、ｎ＝Ｎであると判定した場合（ステップＳ３２：Ｙｅｓ）、収束判定部による判定処理に進む。

収束判定部は、収束したかどうかの判定を行う（ステップＳ３３）。信号分析装置２０１は、収束していないと収束判定部が判定した場合（ステップＳ３３：Ｎｏ）、ステップＳ２５に戻って、処理を継続する。一方、収束したと収束判定部が判定した場合（ステップＳ３３：Ｙｅｓ）、音源信号事後確率更新部２１３は、事後確率の平均ξ_ｎ（ｔ，ｆ）を、音源信号成分ｘ_ｎ（ｔ，ｆ）の推定値＾ｘ_ｎ（ｔ，ｆ）として出力し（ステップＳ３４）、信号分析装置２０１での処理が終了する。

［第１の実施形態の変形例５］
第１の実施形態では、空間共分散行列により音源信号の空間的特性をモデル化したが、他のパラメータにより音源信号の空間的特性をモデル化してもよい。音源信号の空間的特性をモデル化するパラメータを、ここでは空間パラメータと呼ぶ。

例えば、空間パラメータとしてステアリングベクトルを用い、これにより音源信号の空間的特性をモデル化してもよい。この場合、観測信号ベクトルｙ（ｔ，ｆ）の確率分布は、例えば次の（８７）式の複素ガウス分布によりモデル化できる。

ここで、ｈ_ｎ（ｆ）は、音源信号ｎの空間的特性をモデル化する空間パラメータであるステアリングベクトルであり、σ_１ ^２は正則化のための正数である。この場合、ｈ_ｎ（ｆ）の事前分布は次の（８８）式で与えられる。但し、（８８）式における「ｐ」は、複素ガウス分布「ｐ_Ｇ」を表す。

ここで、ｇ_ｋ（ｆ）とσ_２ ^２はハイパーパラメータである。ｇ_ｋ（ｆ）はｋ番目の音源位置候補に対するステアリングベクトルであり、σ_２ ^２は正則化のための正数である。以上のモデル化に基づいて、第１の実施形態と同様にパラメータΘを推定すればよい。

［システム構成等］
また、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示のように構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部又は一部を、各種の負荷や使用状況等に応じて、任意の単位で機能的又は物理的に分散・統合して構成することができる。さらに、各装置にて行われる各処理機能は、その全部又は任意の一部が、ＣＰＵおよび当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。

また、本実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的に行うこともでき、あるいは、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。この他、上記文書中や図面中で示した処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて任意に変更することができる。すなわち、上記学習方法および音声認識方法において説明した処理は、記載の順にしたがって時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。

［プログラム］
図５は、プログラムが実行されることにより、信号分析装置１，２０１が実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。これらの各部は、バス１０８０によって接続される。

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

ハードディスクドライブ１０９０は、例えば、ＯＳ（Operating System）１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわち、信号分析装置１，２０１の各処理を規定するプログラムは、コンピュータ１０００により実行可能なコードが記述されたプログラムモジュール１０９３として実装される。プログラムモジュール１０９３は、例えばハードディスクドライブ１０９０に記憶される。例えば、信号分析装置１，２０１における機能構成と同様の処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０９０に記憶される。なお、ハードディスクドライブ１０９０は、ＳＳＤ（Solid State Drive）により代替されてもよい。

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

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

以上、本発明者によってなされた発明を適用した実施形態について説明したが、本実施形態による本発明の開示の一部をなす記述および図面により本発明は限定されることはない。すなわち、本実施形態に基づいて当業者等によりなされる他の実施形態、実施例および運用技術等はすべて本発明の範疇に含まれる。

１，２０１，１Ｐ 信号分析装置
１０ 推定部
１１，１１Ｐ 観測信号ベクトル作成部
１２，１２Ｐ 音源存在事後確率更新部
１３，１３Ｐ 記憶部
１４，２１２ 音源位置事後確率更新部
１４Ｐ 音源存在事前確率更新部
１５ 音源存在事前確率更新部
１６，２１４ 音源位置事前確率更新部
１７，２１７，１５Ｐ 空間共分散行列更新部
１８，２１８，１６Ｐ パワーパラメータ更新部
１９，１７Ｐ 音源信号成分推定部
２１３ 音源信号事後確率更新部

Claims (6)

1. Ｎ個（Ｎは２以上の整数）の信号源からの信号の空間的特性をモデル化するパラメータを空間パラメータとする場合、前記空間パラメータの各信号源に対する事前分布を、前記空間パラメータのＫ個（Ｋは２以上の整数）の各信号源位置候補に対する事前分布の線形結合である混合分布によりモデル化するときの混合重みであり、前記信号源ごとの前記各信号源位置候補から信号が到来する確率である、信号源位置事前確率を推定する推定部を有することを特徴とする信号分析装置。
2. 前記空間パラメータは、空間共分散行列であり、
   前記混合分布は混合複素逆ウィシャート分布であることを特徴する請求項１に記載の信号分析装置。
3. 前記推定部は、未知のパラメータの事後確率を最大化するための目的関数についての補助関数であり、前記目的関数に含まれる前記線形結合における和演算が対数演算の中に含まれない補助関数を用いた補助関数法により前記信号源位置事前確率を推定することを特徴とする請求項１または２に記載の信号分析装置。
4. 前記推定部は、実際の信号源の数Ｎ´に対し十分に大きい数で仮定する信号源の数をＮとしたとき、各ｎ（ｎは１以上Ｎ以下の整数）に対し、前記信号源位置事前確率が最大となる信号源の位置候補を信号源位置の推定値とし、得られたＮ個の信号源の位置を、階層クラスタリングによりクラスタリングし、得られたクラスタの個数を、実際の音源数Ｎ´の推定値とすることを特徴とする請求項１〜３のいずれか一つに記載の信号分析装置。
5. 信号分析装置が実行する信号分析方法であって、
   Ｎ個（Ｎは２以上の整数）の信号源からの信号の空間的特性をモデル化するパラメータを空間パラメータとする場合、前記空間パラメータの各信号源に対する事前分布を、前記空間パラメータのＫ個（Ｋは２以上の整数）の各信号源位置候補に対する事前分布の線型結合である混合分布によりモデル化するときの混合重みであり、前記信号源ごとの前記各信号源位置候補から信号が到来する確率である、信号源位置事前確率を推定する工程
   を含んだことを特徴とする信号分析方法。
6. コンピュータを、請求項１〜４のいずれか一つに記載の信号分析装置として機能さ
   せるための信号分析プログラム。

Images