(57)【要約】【課題】各音源信号に対する音源位置が未知である場合にも、空間パラメータの事前分布に基づいて音源分離などの信号分析を行うことができる。【解決手段】信号分析装置１は、Ｎ個（Ｎは２以上の整数）の信号源からの信号の空間的特性をモデル化するパラメータを空間パラメータとする場合、空間パラメータの各信号源に対する事前分布を、空間パラメータのＫ個（Ｋは２以上の整数）の各信号源位置候補に対する事前分布の線型結合である混合分布によりモデル化するときの混合重みであり、信号源ごとの各信号源位置候補から信号が到来する確率である、信号源位置事前確率を推定する推定部１０を有する。【選択図】図１
2019184773.tif 000002
【特許請求の範囲】
【請求項１】
Ｎ個（Ｎは２以上の整数）の信号源からの信号の空間的特性をモデル化するパラメータを空間パラメータとする場合、前記空間パラメータの各信号源に対する事前分布を、前記空間パラメータのＫ個（Ｋは２以上の整数）の各信号源位置候補に対する事前分布の線形結合である混合分布によりモデル化するときの混合重みであり、前記信号源ごとの前記各信号源位置候補から信号が到来する確率である、信号源位置事前確率を推定する推定部を有することを特徴とする信号分析装置。
【請求項２】
前記空間パラメータは、空間共分散行列であり、
前記混合分布は混合複素逆ウィシャート分布であることを特徴する請求項１に記載の信号分析装置。
【請求項３】
前記推定部は、未知のパラメータの事後確率を最大化するための目的関数についての補助関数であり、前記目的関数に含まれる前記線形結合における和演算が対数演算の中に含まれない補助関数を用いた補助関数法により前記信号源位置事前確率を推定することを特徴とする請求項１または２に記載の信号分析装置。
【請求項４】
前記推定部は、実際の信号源の数Ｎ´に対し十分に大きい数で仮定する信号源の数をＮとしたとき、各ｎ（ｎは１以上Ｎ以下の整数）に対し、前記信号源位置事前確率が最大となる信号源の位置候補を信号源位置の推定値とし、得られたＮ個の信号源の位置を、階層クラスタリングによりクラスタリングし、得られたクラスタの個数を、実際の音源数Ｎ´の推定値とすることを特徴とする請求項１〜３のいずれか一つに記載の信号分析装置。
【請求項５】
信号分析装置が実行する信号分析方法であって、
Ｎ個（Ｎは２以上の整数）の信号源からの信号の空間的特性をモデル化するパラメータを空間パラメータとする場合、前記空間パラメータの各信号源に対する事前分布を、前記空間パラメータのＫ個（Ｋは２以上の整数）の各信号源位置候補に対する事前分布の線型結合である混合分布によりモデル化するときの混合重みであり、前記信号源ごとの前記各信号源位置候補から信号が到来する確率である、信号源位置事前確率を推定する工程
を含んだことを特徴とする信号分析方法。
【請求項６】
コンピュータを、請求項１〜４のいずれか一つに記載の信号分析装置として機能さ
せるための信号分析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号分析装置、信号分析方法および信号分析プログラムに関する。
【背景技術】
【０００２】
Ｎ´個（Ｎ´は０以上の整数）の音源信号が混在する状況において、それぞれ異なる位置でマイクロホンにより取得された複数の観測信号から、個々の音源信号を推定する音源分離技術がある。Ｎ´は真の音源数であり、Ｎは仮定された音源数であるとする。従来技術では、真の音源数Ｎ´が既知である状況を想定し、仮定された音源数をＮ＝Ｎ´と設定する。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】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.
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ここで、図６および図７を用いて、従来の音源分離装置の構成と処理について説明する。図６は、従来の音源分離装置の構成の一例を示す図である。図７は、従来の音源分離処理の処理手順の一例を示すフローチャートである。なお、ベクトル、行列又はスカラーであるＡに対し、“＾Ａ”と記載する場合は「“Ａ”の直上に“＾”が記された記号」と同じであるとする。また、ベクトル、行列又はスカラーであるＡに対し、“~Ａ”と記載する場合は「“Ａ”の直上に“~”が記された記号」と同じであるとする。
【０００５】
図６及び図７に示すように、従来の信号分析装置１Ｐは、観測信号ベクトル作成部１１Ｐ、初期化部（図示しない）、音源存在事後確率更新部１２Ｐ、記憶部１３Ｐ、音源存在事前確率更新部１４Ｐ、空間共分散行列更新部１５Ｐ、パワーパラメータ更新部１６Ｐ、収束判定部（図示しない）および音源信号成分推定部１７Ｐを有する。
【０００６】
観測信号ベクトル作成部１１Ｐは、まず、入力された観測信号ｙｍ（τ）を取得し（ステップＳ４１）、短時間フーリエ変換などにより時間周波数領域の観測信号ｙｍ（ｔ，ｆ）を計算する（ステップＳ４２）。ここで、ｔ＝１，・・・，Ｔはフレームのインデックスであり、ｆ＝１，・・・，Ｆは周波数ビンのインデックスであり、ｍ＝１，・・・，Ｍはマイクロホンのインデックスであり、τはサンプル点のインデックスである。Ｍ個のマイクロホンはそれぞれ異なる位置に配置されているとする。  
【０００７】
次に、観測信号ベクトル作成部１１Ｐは、（１）式のように、取得されたＭ個すべての観測信号ｙｍ（ｔ，ｆ）からなるＭ次元縦ベクトルである観測信号ベクトルｙ（ｔ，ｆ）を時間周波数点ごとに作成する（ステップＳ４３）。ここで、上付きのＴは転置を表す。
【０００８】
【数１】
2019184773.tif 000003
【０００９】
初期化部は、音源存在事前確率αｎ（ｆ）と、空間共分散行列Ｒｎ（ｆ）と、パワーパラメータｖｎ（ｔ，ｆ）と、の推定値の初期値を計算することでこれらのパラメータを初期化する（ステップＳ４４）。ただし、ｎ＝１，・・・，Ｎは音源のインデックスである。例えば、初期化部は、乱数に基づいてこれらの初期値を計算する。
【００１０】
音源存在事後確率更新部１２Ｐは、観測信号ベクトル作成部１１Ｐからの観測信号ベクトルｙ（ｔ，ｆ）と、音源存在事前確率更新部１４Ｐからの音源存在事前確率（ただし例外として、音源存在事後確率更新部１２Ｐにおける最初の処理の際には、初期化部からの音源存在事前確率の初期値）αｎ（ｆ）と、空間共分散行列更新部１５Ｐからの空間共分散行列（ただし例外として、音源存在事後確率更新部１２Ｐにおける最初の処理の際には、初期化部からの空間共分散行列の初期値）Ｒｎ（ｆ）と、パワーパラメータ更新部からのパワーパラメータ（ただし例外として、音源存在事後確率更新部１２Ｐにおける最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖｎ（ｔ，ｆ）と、を受け取って、音源存在事後確率λｎ（ｔ，ｆ）を更新する（ステップＳ４５）。
【００１１】
記憶部１３Ｐは、各音源信号ｎおよび各周波数ビンｆに対する空間共分散行列の事前分布のパラメータを記憶する。
【００１２】
音源存在事前確率更新部１４Ｐは、音源存在事後確率更新部１２Ｐからの音源存在事後確率λｎ（ｔ，ｆ）を受け取って、音源存在事前確率αｎ（ｆ）を更新する（ステップＳ４６）。
【００１３】
空間共分散行列更新部１５Ｐは、観測信号ベクトル作成部１１Ｐからの観測信号ベクトルｙ（ｔ，ｆ）と、音源存在事後確率更新部１２Ｐからの音源存在事後確率λｎ（ｔ，ｆ）と、記憶部１３Ｐからの事前分布のパラメータと、パワーパラメータ更新部１６Ｐからのパワーパラメータ（ただし例外として、空間共分散行列更新部１５Ｐにおける最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖｎ（ｔ，ｆ）と、を受け取って、空間共分散行列Ｒｎ（ｆ）を更新する（ステップＳ４７）。
【００１４】
パワーパラメータ更新部１６Ｐは、観測信号ベクトル作成部１１Ｐからの観測信号ベクトルｙ（ｔ，ｆ）と、空間共分散行列更新部１５Ｐからの空間共分散行列Ｒｎ（ｆ）と、を受け取って、パワーパラメータｖｎ（ｔ，ｆ）を更新する（ステップＳ４８）。
【００１５】
収束判定部は、収束したかどうかの判定を行う（ステップＳ４９）。収束判定部によって収束していないと判定された場合（ステップＳ４９：Ｎｏ）、音源存在事後確率更新部１２Ｐでの処理（ステップＳ４５）に戻って、処理が継続される。一方、収束判定部によって収束したと判定された場合（ステップＳ４９：Ｙｅｓ）、音源信号成分推定部１７Ｐでの処理に進む。
【００１６】
音源信号成分推定部１７Ｐは、観測信号ベクトル作成部１１Ｐからの観測信号ベクトルｙ（ｔ，ｆ）と音源存在事後確率更新部１２Ｐからの音源存在事後確率λｎ（ｔ，ｆ）とを受け取って、音源信号成分ｘｎ（ｔ，ｆ）の推定値＾ｘｎ（ｔ，ｆ）を計算し、出力する（ステップＳ５０）。
【００１７】
ここで、従来技術の特徴について説明する。観測信号ベクトル作成部１１Ｐにおいて作成された観測信号ベクトルｙ（ｔ，ｆ）は、Ｎ個の音源信号に由来する成分である音源信号成分ｘ１（ｔ，ｆ），・・・，ｘＮ（ｔ，ｆ）の和として、（２）式で表される。
【００１８】
【数２】
2019184773.tif 000004
【００１９】
従来技術では、各音源信号は、時間周波数領域において、疎な点においてのみ有意なエネルギーを持つという性質（スパース性）を有すると仮定する。例えば、音声はこのスパース性を比較的よく満たすとされている。この仮定の下では、各時間周波数点では、観測信号ベクトルｙ（ｔ，ｆ）は、Ｎ個の音源信号成分ｘ１（ｔ，ｆ），・・・，ｘＮ（ｔ，ｆ）のうちの一つだけからなると近似できる（（３）式）。
【００２０】
【数３】
2019184773.tif 000005
【００２１】
ここで、ｎ（ｔ，ｆ）は、時間周波数点（ｔ，ｆ）において存在する音源信号のインデックスであり、１以上Ｎ以下の整数の値を取る。
【００２２】
（３）式のモデルの下では、各時間周波数点（ｔ，ｆ）において存在する音源信号のインデックスｎ（ｔ，ｆ）の推定値＾ｎ（ｔ，ｆ）が得られれば、音源分離を実現できる。すなわち、一旦＾ｎ（ｔ，ｆ）が得られれば、次の（４）式のように、ｎ番目の音源信号が存在する時間周波数点以外の音のエネルギーを遮断するかまたは減衰させることにより、ｎ番目の音源信号成分ｘｎ（ｔ，ｆ）の推定値＾ｘｎ（ｔ，ｆ）を得ることができる、すなわち、音源分離が実現できる。
【００２３】
【数４】
2019184773.tif 000006
【００２４】
従来技術では、観測信号ベクトルｙ（ｔ，ｆ）の確率分布を次の（５）式の混合複素ガウス分布でモデル化し、観測信号ベクトルｙ（ｔ，ｆ）にこのモデルを当てはめることにより、ｎ（ｔ，ｆ）の推定を実現する。
【００２５】
【数５】
2019184773.tif 000007
【００２６】
ここで、ｐＧは複素ガウス分布を表す（Ｇはガウス（Gauss）の頭文字である）。Ｒｎ（ｆ）は、各音源の空間的特性（音響伝達特性）を表すパラメータである空間共分散行列であり、ｖｎ（ｔ，ｆ）は、各音源のパワースペクトルをモデル化するパラメータであるパワーパラメータである。αｎ（ｆ）は、（６）式を満たす混合重みであり、本明細書では音源存在事前確率とも呼ぶ。
【００２７】
【数６】
2019184773.tif 000008
【００２８】
また、Θは、すべての未知パラメータをまとめて表したものであり、具体的には、音源存在事前確率αｎ（ｆ）、空間共分散行列Ｒｎ（ｆ）、パワーパラメータｖｎ（ｔ，ｆ）からなる。ひとたびパラメータΘが推定できれば、観測信号ベクトルｙ（ｔ，ｆ）が与えられた下での音源インデックスｎ（ｔ，ｆ）の事後確率を、次の（７）式により求めることができる。
【００２９】
【数７】
2019184773.tif 000009
【００３０】
これを用いて、次の（８）式のように音源インデックスｎ（ｔ，ｆ）を推定することができる。
【００３１】
【数８】
2019184773.tif 000010
【００３２】
この音源インデックスの推定値を用いれば、（４）式に従って、音源分離を実現できる。
【００３３】
このアプローチに基づいて高精度な音源分離を実現するためにはパラメータΘの正確な推定が鍵となる。一般に、与えられる観測信号の長さが長ければ長いほどパラメータΘの正確な推定が容易になり、与えられる観測信号の長さが短ければ短いほどパラメータΘの正確な推定が困難になる。そこで、与えられる観測信号の長さが短くなった場合におけるパラメータΘの推定精度の劣化を防ぐために、パラメータΘに関する事前知識を表す事前分布を適切に定めることが重要である。事前分布を適切に定めれば、与えられる観測信号の長さが短くなった場合でも、パラメータΘに関する事前知識に基づいて、パラメータΘをある程度正確に推定できるため、パラメータΘの推定精度の大幅な低下を防ぐことができる。事前分布はまた、オンライン処理における音源信号が鳴り始めた直後におけるパラメータの推定精度の劣化防止や、パーミュテーション問題の回避のためにも重要である。
【００３４】
ここで、パーミュテーション問題について説明する。観測信号ベクトルｙ（ｔ，ｆ）は、周波数ビンごとに異なる分布に従う。このため、（５）式のような混合モデルを用いた音源インデックスｎ（ｔ，ｆ）の推定（クラスタリング）に基づく音源分離アプローチでは、一般に、各周波数ビン内に限定した音源の分類（クラスタリング）はできても、異なる周波数間で音源の対応をとることはできない。これが、パーミュテーション問題と呼ばれている。
【００３５】
従来技術では、各音源の音源位置が既知であるという仮定の下、各音源信号の空間的特性をモデル化するパラメータである空間共分散行列Ｒｎ（ｆ）の事前分布ｐ（Ｒｎ（ｆ））を設計していた。具体的には、従来技術では、空間共分散行列Ｒｎ（ｆ）の事前分布ｐ（Ｒｎ（ｆ））を、次の（９）式の逆ウィシャート分布によりモデル化する。
【００３６】
【数９】
2019184773.tif 000011
【００３７】
ここで、ＩＷは、逆ウィシャート分布を表す（「ＩＷ」は、「Inverse  Wishart（逆ウィシャート）」の頭文字である）。~Ψｎ（ｆ）は事前分布ｐ（Ｒｎ（ｆ））の山（モード）の位置をモデル化するスケール行列であり、~νｎ（ｆ）は事前分布ｐ（Ｒｎ（ｆ））の山の広がりをモデル化する自由度である。以下、自由度~νｎ（ｆ）は音源および周波数ビンに依らず一定であると仮定し、単に~νと書く。事前分布ｐ（Ｒｎ（ｆ））のパラメータであるスケール行列~Ψｎ（ｆ）および自由度~νは、パラメータＲｎ（ｆ）をモデル化するパラメータであり、その意味でハイパーパラメータと呼ばれる。
【００３８】
（９）式より、すべての周波数ビンにおける空間共分散行列Ｒｎ（１），・・・，Ｒｎ（Ｆ）の事前分布ｐ（Ｒｎ（１），・・・，Ｒｎ（Ｆ））は、次の（１０）式のようになる。
【００３９】
【数１０】
2019184773.tif 000012
【００４０】
ここで周波数間の独立性を仮定した。
【００４１】
従来技術では、各音源の音源位置が既知であるという仮定の下、事前分布ｐ（Ｒｎ（ｆ））のハイパーパラメータであるスケール行列~Ψｎ（ｆ）および自由度~νを既知であると仮定していた。これらのハイパーパラメータは、学習用データに基づいて、事前に学習することができる。すなわち、各音源の音源位置が既知の場合には、音源ごとに既知である音源位置から音源信号が到来する場合の観測信号を実測し、これを学習用データとして用いることにより、事前分布ｐ（Ｒｎ（ｆ））のハイパーパラメータであるスケール行列~Ψｎ（ｆ）および自由度~νを事前に学習することができる。
【００４２】
従来技術では、この事前分布に基づき、以下の（１１）式〜（１４）式に示す更新則を交互に繰り返し適用することにより、パラメータΘを推定する。
【００４３】
【数１１】
2019184773.tif 000013
【数１２】
2019184773.tif 000014
【数１３】
2019184773.tif 000015
【数１４】
2019184773.tif 000016
【００４４】
（１１）式の処理は音源存在事後確率更新部１２Ｐにおいて、（１２）式の処理は音源存在事前確率更新部１４Ｐにおいて、（１３）式の処理は空間共分散行列更新部１５Ｐにおいて、（１４）式の処理はパワーパラメータ更新部１６Ｐにおいて、それぞれ行われる。音源信号成分推定部１７Ｐは、上記の処理により得られた音源存在事後確率更新部１２Ｐからの音源存在事後確率λｎ（ｔ，ｆ）に基づいて、（８）式により音源インデックスの推定値＾ｎ（ｔ，ｆ）を計算し、さらに（４）式により音源信号成分の推定値＾ｘｎ（ｔ，ｆ）を計算する。
【００４５】
しかしながら、従来技術では、各音源信号に対する音源位置が既知であると仮定しており、各音源信号に対する音源位置が未知である場合には適用できなかった。
【００４６】
本発明は、上記に鑑みてなされたものであって、各音源信号に対する音源位置が未知である場合にも、各音源信号の空間的特性をモデル化するパラメータである空間パラメータ（例えば、空間共分散行列）の事前分布に基づいて音源分離などの信号分析を行うことができる信号分析装置、信号分析方法および信号分析プログラムを提供することを目的とする。
【課題を解決するための手段】
【００４７】
上述した課題を解決し、目的を達成するために、本発明の信号分析装置は、Ｎ個（Ｎは２以上の整数）の信号源からの信号の空間的特性をモデル化するパラメータを空間パラメータとする場合、空間パラメータの各信号源に対する事前分布を、空間パラメータのＫ個（Ｋは２以上の整数）の各信号源位置候補に対する事前分布の線型結合である混合分布によりモデル化するときの混合重みであり、信号源ごとの各信号源位置候補から信号が到来する確率である、信号源位置事前確率を推定する推定部を有することを特徴とする。
【発明の効果】
【００４８】
本発明によれば、各音源信号に対する音源位置が未知である場合にも、空間パラメータの事前分布に基づいて音源分離などの信号分析を行うことができる。
【図面の簡単な説明】
【００４９】
【図１】図１は、第１の実施形態に係る信号分析装置の構成の一例を示す図である。
【図２】図２は、第１の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【図３】図３は、第１の実施形態の変形例４に係る信号分析装置の構成の一例を示す図である。
【図４】図４は、第１の実施形態の変形例４に係る信号分析処理の処理手順の一例を示すフローチャートである。
【図５】図５は、プログラムが実行されることにより、信号分析装置が実現されるコンピュータの一例を示す図である。
【図６】図６は、従来の音源分離装置の構成の一例を示す図である。
【図７】図７は、従来の音源分離処理の処理手順の一例を示すフローチャートである。
【発明を実施するための形態】
【００５０】
以下に、本願に係る信号分析装置、信号分析方法および信号分析プログラムの実施形態を図面に基づいて詳細に説明する。また、本発明は、以下に説明する実施形態により限定されるものではない。なお、以下では、ベクトル、行列又はスカラーであるＡに対し、“＾Ａ”と記載する場合は「“Ａ”の直上に“＾”が記された記号」と同じであるとする。また、ベクトル、行列又はスカラーであるＡに対し、“~Ａ”と記載する場合は「“Ａ”の直上に“~”が記された記号」と同じであるとする。
【００５１】
［第１の実施形態］
まず、第１の実施形態に係る信号分析装置について説明する。なお、第１の実施形態においては、Ｎ´個（Ｎ´は０以上の整数）の音源信号が混在する状況において、それぞれ異なる位置でマイクロホンにより取得されたＭ個（Ｍは２以上の整数）の観測信号ｙｍ（τ）（ｍ＝１，・・・，Ｍはマイクロホンのインデックス、τはサンプル点のインデックス）が信号分析装置に入力されるものとする。Ｎ´は真の音源数であり、Ｎは仮定された音源数であるとする。第１の実施形態では、真の音源数Ｎ´が既知である状況を想定し、仮定された音源数をＮ＝Ｎ´と設定する。なお、本第１の実施形態における「音源信号」は、目的信号（例えば、音声）であってもよいし、特定の音源位置から到来する雑音である方向性雑音（例えば、テレビから流れる音楽）であってもよい。また、様々な音源位置から到来する雑音である拡散性雑音を、まとめて１つの「音源信号」とみなしてもよい。拡散性雑音の例としては、雑踏やカフェ等における大勢の人々の話し声、駅や空港における足音、空調による雑音などが挙げられる。
【００５２】
図１および図２を用いて、第１の実施形態の構成と処理について説明する。図１は、第１の実施形態に係る信号分析装置の構成の一例を示す図である。図２は、第１の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。第１の実施形態に係る信号分析装置１は、例えば、ＲＯＭ（Read  Only  Memory）、ＲＡＭ（Random  Access  Memory）、ＣＰＵ（Central  Processing  Unit）等を含むコンピュータ等に所定のプログラムが読み込まれて、ＣＰＵが所定のプログラムを実行することで実現される。
【００５３】
図１に示すように、信号分析装置１は、観測信号ベクトル作成部１１、初期化部（図示しない）、推定部１０、記憶部１３、パワーパラメータ更新部１８、パーミュテーション解決部（図示しない）、収束判定部（図示しない）、音源信号成分推定部１９を有する。
【００５４】
まず、信号分析装置１の各部の概要について説明する。観測信号ベクトル作成部１１は、まず、入力された観測信号ｙｍ（τ）を取得し（ステップＳ１）、短時間フーリエ変換などにより時間周波数領域の観測信号ｙｍ（ｔ，ｆ）を計算する（ステップＳ２）。ここで、ｔ＝１，・・・，Ｔはフレームのインデックスであり、ｆ＝１，・・・，Ｆは周波数ビンのインデックスである。
【００５５】
次に、観測信号ベクトル作成部１１は、取得されたＭ個すべての観測信号ｙｍ（ｔ，ｆ）からなるＭ次元縦ベクトルである観測信号ベクトルｙ（ｔ，ｆ）、すなわち（１５）式で表される観測信号ベクトルｙ（ｔ，ｆ）、を時間周波数点ごとに作成する（ステップＳ３）。ここで、上付きのＴは転置を表す。
【００５６】
【数１５】
2019184773.tif 000017
【００５７】
本第１の実施形態では、各音源信号は、Ｋ個の音源位置の候補のいずれかから到来すると仮定し、それらの音源位置候補をインデックス（以下、「音源位置インデックス」）１，・・・，Ｋで表す。例えば、音源が円卓の周りに着席して会話している複数の話者であり、Ｍ個のマイクロホンが円卓の中央の数ｃｍ四方程度の小領域内に置かれており、音源位置として円卓の中央から見たときの音源の方位角のみに注目するとき、０°〜３６０°をＫ等分したＫ個の方位角Δφ，２Δφ，・・・，ＫΔφ（Δφ＝３６０°／Ｋ）を音源位置候補とすることができる。この例に限らず、一般に任意の所定のＫ点を、音源位置候補として指定することができる。また、音源位置候補は、拡散性雑音を表す音源位置候補でもよい。拡散性雑音は、１つの音源位置から到来するのではなく、多数の音源位置から到来する。このような拡散性雑音も「多数の音源位置から到来する」という１つの音源位置候補とみなすことにより、拡散性雑音が存在する状況でも正確な推定が可能になる。
【００５８】
初期化部は、音源存在事前確率αｎ（ｆ）と、音源位置事前確率βｋｎと、空間共分散行列Ｒｎ（ｆ）と、パワーパラメータｖｎ（ｔ，ｆ）と、の推定値の初期値を計算する（ステップＳ４）。ただし、ｎ＝１，・・・，Ｎは音源のインデックス、ｋ＝１，・・・・，Ｋは音源位置インデックスである。例えば、初期化部は、乱数に基づいてこれらの初期値を計算する。
【００５９】
推定部１０は、音源位置事前確率を推定する。本第１の実施形態では、Ｎ個の音源の位置からの信号の空間的特性をモデル化するパラメータである空間パラメータとして、空間共分散行列を用いる。音源位置事前確率は、空間共分散行列（空間パラメータ）の各音源に対する事前分布を、空間共分散行列（空間パラメータ）のＫ個（Ｋは２以上の整数）の各音源位置候補に対する事前分布の線型結合である混合分布によりモデル化するときの混合重みであり、音源ごとの各音源位置候補から信号が到来する確率である。推定部１０は、音源存在事後確率更新部１２、音源位置事後確率更新部１４、音源存在事前確率更新部１５、音源位置事前確率更新部１６および空間共分散行列更新部１７を有する。
【００６０】
音源存在事後確率更新部１２は、観測信号ベクトル作成部１１からの観測信号ベクトルｙ（ｔ，ｆ）と、音源存在事前確率更新部１５からの音源存在事前確率（ただし例外として、音源存在事後確率更新部１２における最初の処理の際には、初期化部からの音源存在事前確率の初期値）αｎ（ｆ）と、空間共分散行列更新部１７からの空間共分散行列（ただし例外として、音源存在事後確率更新部１２における最初の処理の際には、初期化部からの空間共分散行列の初期値）Ｒｎ（ｆ）と、パワーパラメータ更新部１８からのパワーパラメータ（ただし例外として、音源存在事後確率更新部１２における最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖｎ（ｔ，ｆ）と、を受け取って、音源存在事後確率λｎ（ｔ，ｆ）を更新する（ステップＳ５）。
【００６１】
記憶部１３は、各音源位置候補ｋ、各周波数ビンｆに対する空間共分散行列の事前分布のパラメータを記憶する。
【００６２】
音源位置事後確率更新部１４は、記憶部１３からの事前分布のパラメータと、音源位置事前確率更新部１６からの音源位置事前確率（ただし例外として、音源位置事後確率更新部１４における最初の処理の際には、初期化部からの音源位置事前確率の初期値）βｋｎと、空間共分散行列更新部１７からの空間共分散行列（ただし例外として、音源位置事後確率更新部１４における最初の処理の際には、初期化部からの空間共分散行列の初期値）Ｒｎ（ｆ）と、を受け取って、音源位置事後確率μｋｎを更新する（ステップＳ６）。
【００６３】
音源存在事前確率更新部１５は、音源存在事後確率更新部１２からの音源存在事後確率λｎ（ｔ，ｆ）を受け取って、音源存在事前確率αｎ（ｆ）を更新する（ステップＳ７）。
【００６４】
音源位置事前確率更新部１６は、音源位置事後確率更新部１４からの音源位置事後確率μｋｎを受け取って、音源位置事前確率βｋｎを更新する（ステップＳ８）。
【００６５】
空間共分散行列更新部１７は、観測信号ベクトル作成部１１からの観測信号ベクトルｙ（ｔ，ｆ）と、音源存在事後確率更新部１２からの音源存在事後確率λｎ（ｔ，ｆ）と、記憶部１３からの事前分布のパラメータと、音源位置事後確率更新部１４からの音源位置事後確率μｋｎと、パワーパラメータ更新部１８からのパワーパラメータ（ただし例外として、空間共分散行列更新部１７における最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖｎ（ｔ，ｆ）と、を受け取って、空間共分散行列Ｒｎ（ｆ）を更新する（ステップＳ９）。
【００６６】
パワーパラメータ更新部１８は、観測信号ベクトル作成部１１からの観測信号ベクトルｙ（ｔ，ｆ）と、空間共分散行列更新部１７からの空間共分散行列Ｒｎ（ｆ）と、を受け取って、パワーパラメータｖｎ（ｔ，ｆ）を更新する（ステップＳ１０）。
【００６７】
パーミュテーション解決部は、音源存在事前確率更新部１５からの音源存在事前確率αｎ（ｆ）と、空間共分散行列更新部１７からの空間共分散行列Ｒｎ（ｆ）と、パワーパラメータ更新部１８からのパワーパラメータｖｎ（ｔ，ｆ）と、を受け取り、音源存在事前確率αｎ（ｆ）と、空間共分散行列Ｒｎ（ｆ）と、パワーパラメータｖｎ（ｔ，ｆ）と、を更新することでパーミュテーション問題を解決する（ステップＳ１１）。具体的には、パーミュテーション解決部は、尤度または対数尤度または補助関数などの評価値が最大となるように、音源インデックスｎを周波数ビンごとに付け替えることにより、これらのパラメータを更新する。すなわち、周波数ビンｆにおける音源インデックスｎの付け替えを全単射σf：｛１，・・・，Ｎ｝→｛１，・・・，Ｎ｝で表すとき、各周波数ビンｆにおいてこれらのパラメータの音源インデックスｎをσf（ｎ）に付け替えたときの尤度または対数尤度または補助関数などの評価値が最大になるように全単射σfを求め、求めた全単射σfを用いて各周波数ビンｆにおいてこれらのパラメータの音源インデックスｎをσf（ｎ）に付け替えることにより、これらのパラメータを更新する。なお、パーミュテーション解決部は、音源存在事前確率αｎ（ｆ）と、空間共分散行列Ｒｎ（ｆ）と、パワーパラメータｖｎ（ｔ，ｆ）と、のすべてを更新する代わりに、その一部のみ（例えば、空間共分散行列Ｒｎ（ｆ）のみ）を更新してもよい。なお、パーミュテーション解決部での処理は必須ではない。
【００６８】
続いて、収束判定部は、収束したかどうかの判定を行う（ステップＳ１２）。収束判定部が収束していないと判定した場合（ステップＳ１２：Ｎｏ）、音源存在事後確率更新部１２での処理（ステップＳ５）に戻って、以降の処理が継続される。一方、収束判定部が収束したと判定した場合（ステップＳ１２：Ｙｅｓ）、音源信号成分推定部１９における処理（ステップＳ１３）に進む。
【００６９】
音源信号成分推定部１９は、観測信号ベクトル作成部１１からの観測信号ベクトルｙ（ｔ，ｆ）と音源存在事後確率更新部１２からの音源存在事後確率λｎ（ｔ，ｆ）とを受け取って、音源信号成分ｘｎ（ｔ，ｆ）の推定値＾ｘｎ（ｔ，ｆ）を計算し、出力する（ステップＳ１３）。
【００７０】
次に、第１の実施形態の特徴について、従来技術と対比しながら説明する。前述の通り、従来技術では、すべての周波数ビンにおける空間共分散行列Ｒｎ（１），・・・，Ｒｎ（Ｆ）の事前分布ｐ（Ｒｎ（１），・・・，Ｒｎ（Ｆ））を、次の（１６）式（（１０）式を再掲）によりモデル化していた。
【００７１】
【数１６】
2019184773.tif 000018
【００７２】
しかしながら、従来技術では、各音源の音源位置が既知であると仮定しており、各音源の音源位置が未知の場合には適用できないという問題があった。
【００７３】
これに対し、本第１の実施形態では、すべての周波数ビンにおける空間共分散行列Ｒｎ（１），・・・，Ｒｎ（Ｆ）の事前分布ｐ（Ｒｎ（１），・・・，Ｒｎ（Ｆ））を、次の（１７）式の混合複素逆ウィシャート分布でモデル化する。
【００７４】
【数１７】
2019184773.tif 000019
【００７５】
これは、音源位置候補ｋに対する事前分布を、音源ｎが音源位置候補ｋにある確率βｋｎを重みとして平均した形になっている。本第１の実施形態では各音源の音源位置が未知であると仮定しているから、βｋｎは未知の確率である。ただし、βｋｎは確率であるから、次の（１８）式を満たすものとする。
【００７６】
【数１８】
2019184773.tif 000020
【００７７】
このように、未知の確率βｋｎによる重み付き和に基づくことで、各音源の音源位置が未知の場合でも、空間共分散行列の事前分布を設計することができる。βｋｎは未知であるが、これも未知パラメータとみなし、他の未知パラメータと同時に推定することができる。
【００７８】
本第１の実施形態では、各音源位置候補ｋ、各周波数ビンｆに対する複素逆ウィシャート分布のパラメータΨｋ（ｆ），νｋ（ｆ）は、事前に準備され、記憶部１３に記憶されているものとする。これらのパラメータは、マイクロホン配置の情報に基づいて事前に準備してもよいし、実測データから事前に学習してもよい。
【００７９】
例えば、マイクロホン配置の情報に基づいて事前に準備する場合には、各マイクロホンｍの直交座標をｒｍとして、各音源位置候補ｋに対応する平面波のステアリングベクトルを（１９）式により計算し、Ψｋ（ｆ），νｋ（ｆ）を次の（２０）式および（２１）式により計算すればよい。
【００８０】
【数１９】
2019184773.tif 000021
【００８１】
【数２０】
2019184773.tif 000022
【数２１】
2019184773.tif 000023
【００８２】
ここで、ｄｋはｋ番目の音源位置候補に対応する音源信号の到来方向を表す単位ベクトル、ｃは音速、ωｆは周波数ビンｆに対応する角周波数、（２１−１）式に示すｊは虚数単位、上付きのＨはエルミート転置である。
【００８３】
【数２２】
2019184773.tif 000024
【００８４】
ここで、本第１の実施形態における事前分布（（１７）式）の導出について説明する。各音源の音源位置は未知であると仮定し、各音源ｎの音源位置に対応する音源位置インデックスｋｎは、（２２）式に示す未知の確率分布に従うと仮定する。βｋｎは、音源ごとの音源位置インデックスの確率分布である音源位置事前確率である。
【００８５】
【数２３】
2019184773.tif 000025
【００８６】
さらに、本第１の実施形態では、音源ｎに対する音源位置インデックスがｋｎ＝ｋであるという条件の下で、音源ｎの空間共分散行列Ｒｎ（１），・・・，Ｒｎ（Ｆ）が、互いに独立に確率分布（（２３）式）に従うものとする。
【００８７】
【数２４】
2019184773.tif 000026
【００８８】
ここで、Ψｋ（ｆ）は、各音源位置候補に対する空間共分散行列の事前分布の山（モード）の位置を表すパラメータ（スケール行列）であり、νｋ（ｆ）は、各音源位置候補に対する空間共分散行列の事前分布の山の広がり（自由度）を表すパラメータである。また、ＩＷＣ（Σ；Ψ，ν）は、（２４）式に示すものであり、スケール行列がΨ、自由度がνである複素逆ウィシャート分布である。
【００８９】
【数２５】
2019184773.tif 000027
【００９０】
（２２）式および（２３）式のモデル化の下では、音源ｎの空間共分散行列Ｒｎ（１），・・・，Ｒｎ（Ｆ）の確率分布は、次の（２５）式〜（２８）式で与えられる。
【００９１】
【数２６】
2019184773.tif 000028
【００９２】
本実施形態では、事前分布（（１７）式）に基づき、パラメータを推定する。以下、本実施形態におけるパラメータ推定アルゴリズムについて説明する。なお、以下では簡単のため、複素逆ウィシャート分布「ＩＷＣ」を、添え字Ｃを省略して単に「ＩＷ」と表す。空間共分散行列Ｒｎ（ｆ）以外の未知パラメータの事前分布は一様分布であると仮定すると、パラメータΘの事前分布は次の（２９）式および（３０）式で与えられる。
【００９３】
【数２７】
2019184773.tif 000029
【００９４】
なお、本第１の実施形態におけるパラメータΘは、音源存在事前確率αｎ（ｆ）、パワーパラメータｖｎ（ｔ，ｆ）、空間共分散行列Ｒｎ（ｆ）および音源位置事前確率βｋｎからなる。
【００９５】
一方、パラメータΘが与えられた下で、各時間周波数点における観測信号ベクトルｙ（ｔ，ｆ）が互いに独立であると仮定すると、尤度が次の（３１）式および（３２）式で与えられる。
【００９６】
【数２８】
2019184773.tif 000030
【００９７】
ここで、Ｙは、すべての時間周波数点における観測信号ベクトルｙ（ｔ，ｆ）をまとめて表したものである。
【００９８】
本第１の実施形態では、パラメータΘの事後確率ｐ（Θ｜Ｙ）を最大化することにより、パラメータΘを推定する。ベイズの定理より、この事後確率は（３３）式のように表せ、両辺の対数を取ると、（３４）式となる。
【００９９】
【数２９】
2019184773.tif 000031
【０１００】
【数３０】
2019184773.tif 000032
【０１０１】
ｌｎｐ（Ｙ）はパラメータΘに依らないから、事後確率ｐ（Θ｜Ｙ）のΘに関する最大化は、次の（３５）式のΘに関する最大化と等価であり、したがって次の（３６）式に示す目的関数Ｊ（Θ）のΘに関する最大化と等価である。
【０１０２】
【数３１】
2019184773.tif 000033
【０１０３】
ここで、＝の上に“ｃ”が記された記号は、パラメータΘに依存しない定数の差を除いて両辺が等しいことを表す記号である。また、「Ａ＝：Ｂ」は、ＢをＡによって定義することを表す。
【０１０４】
上式の目的関数Ｊ（Θ）の最大化は、補助関数法に基づいて行うことができる。補助関数法では、パラメータΘと補助変数と呼ばれる変数Φとの関数である補助関数Ｑ（Θ，Φ）に基づいて、以下の２つのステップを交互に反復する。
１．補助関数Ｑ（Θ，Φ）を補助変数Φに関して最大化することにより、補助変数Φを更新するステップ
２．補助関数Ｑ（Θ，Φ）が減少しないようにパラメータΘを更新するステップ
【０１０５】
ただし、補助関数Ｑ（Θ，Φ）は、次の（３７）式に示す条件を満たすものとする。
【０１０６】
【数３２】
2019184773.tif 000034
【０１０７】
この補助関数法により、目的関数Ｊ（Θ）を単調増加させることができる。すなわち、ｉ回目の反復の結果得られたパラメータΘの推定値をΘ（ｉ）として、（３８）式が成り立つ。
【０１０８】
【数３３】
2019184773.tif 000035
【０１０９】
実際、ｉ回目の反復の結果得られた補助変数Φの値をΦ（ｉ）とすると、（３７）式より、（３９）式および（４０）式が成り立つ。
【０１１０】
【数３４】
2019184773.tif 000036
【０１１１】
しかるに、以下の（４１）式が成り立つから、（３８）式が得られる。
【０１１２】
【数３５】
2019184773.tif 000037
【０１１３】
補助関数法においては、（３７）式を満たすような補助関数Ｑ（Θ，Φ）を設計する必要がある。そのために、本第１の実施形態では、イェンセンの不等式を用いる。ｆを凸関数とし、ｗ１，・・・，ｗＬを（４２）式を満たす非負の数とし、ｘ１，・・・，ｘＬを実数とするとき、（４３）式が成り立つ（等号成立条件はｘ１＝・・・＝ｘＬ）ことが知られている。
【０１１４】
【数３６】
2019184773.tif 000038
【０１１５】
【数３７】
2019184773.tif 000039
【０１１６】
これは、イェンセンの不等式と呼ばれる。特に、ｆ（ｘ）＝−ｌｎｘとすると、（４４）式を得る。
【０１１７】
【数３８】
2019184773.tif 000040
【０１１８】
λ１（ｔ，ｆ），・・・，λＮ（ｔ，ｆ）を（４５）式を満たす非負の数とすると（４４）式より（４６）式および（４７）式が得られる。
【０１１９】
【数３９】
2019184773.tif 000041
【０１２０】
【数４０】
2019184773.tif 000042
【０１２１】
また、μ１ｎ，・・・，μＫｎを（４８）式を満たす非負の数とすると、（４４）式より（４９）式および（５０）式が得られる。
【０１２２】
【数４１】
2019184773.tif 000043
【０１２３】
【数４２】
2019184773.tif 000044
【０１２４】
（４７）式および（５０）式より、（５１）式が得られる。
【０１２５】
【数４３】
2019184773.tif 000045
【０１２６】
よって、（５１）式の右辺を、（５２）式とおくと、（３６）式および（５１）式より、（５３）式が成り立つ。
【０１２７】
【数４４】
2019184773.tif 000046
【０１２８】
【数４５】
2019184773.tif 000047
【０１２９】
ただし、補助変数Φは、λｎ（ｔ，ｆ）とμｋｎとからなるものとする。
【０１３０】
（５１）式の等号成立条件は、（５４）式および（５５）式である。
【０１３１】
【数４６】
2019184773.tif 000048
【０１３２】
これは、次の（５６）式および（５７）式と等価である。
【０１３３】
【数４７】
2019184773.tif 000049
【０１３４】
したがって、（５８）式が成り立つ。
【０１３５】
【数４８】
2019184773.tif 000050
【０１３６】
（５３）式および（５８）式より、（５２）式のＱ（Θ，Φ）が（３７）式を満たすことが分かる。これで、目的関数Ｊ（Θ）に対する補助関数が設計できた。
【０１３７】
本第１の実施形態では、(５２)式の補助関数Ｑ（Θ，Φ）に基づいて、補助変数ΦおよびパラメータΘを次のようにして更新する。まず、補助変数Φの更新は、（５６）式および（５７）式により行えばよい。また、パラメータΘの更新は、次の（５９）式〜（６２）式を用いて行えばよい。
【０１３８】
【数４９】
2019184773.tif 000051
【数５０】
2019184773.tif 000052
【数５１】
2019184773.tif 000053
【数５２】
2019184773.tif 000054
【０１３９】
このように、本第１の実施形態では、（３６）式の目的関数を直接最大化する代わりに、補助関数Ｑ（Θ，Φ）に基づいて、補助関数Ｑ（Θ，Φ）を補助変数Φに関して最大化することによりΦを更新するステップと、補助関数Ｑ（Θ，Φ）が減少しないようにパラメータΘを更新するステップと、を交互に反復することにより、（３６）式の目的関数を間接的に最大化する。（３６）式の目的関数においては、対数ｌｎの中にｋに関する和Σｋ＝１Ｋが含まれており、（３６）式の目的関数の各パラメータに関する微分が煩雑な形になるため、（３６）式の目的関数を勾配法などにより直接最大化しようとすると、更新則が煩雑な形になる。これに対し、補助関数Ｑ（Θ，Φ）では、ｋに関する和Σｋ＝１Ｋが対数ｌｎの外に出た形になっており、補助関数Ｑ（Θ，Φ）の各パラメータに関する微分が単純な形になる。また、勾配法では、反復ごとのパラメータの更新量を定めるステップサイズを調整する必要があるが、補助関数法では、ステップサイズが不要であるため、ステップサイズを調整する必要がない。
【０１４０】
（５６）式により更新されたλｎ（ｔ，ｆ）は、観測信号ベクトルｙ（ｔ，ｆ）が観測された「後」の音源存在確率に他ならない。実際、ベイズの定理より、（５６）式は（６３）式とも書ける。
【０１４１】
【数５３】
2019184773.tif 000055
【０１４２】
そこで、λｎ（ｔ，ｆ）を音源存在事後確率と呼ぶ。これに対し、αｎ（ｆ）（（６４）式））は、観測信号ベクトルｙ（ｔ，ｆ）が観測される「前」の音源存在確率であるから、音源存在事前確率と呼ぶ。
【０１４３】
【数５４】
2019184773.tif 000056
【０１４４】
また、（５７）式により更新されたμｋｎは、空間共分散行列Ｒｎ（１），・・・，Ｒｎ（Ｆ）が与えられた「後」の音源位置確率に他ならない。実際、（５７）は、（６５）式とも書ける。
【０１４５】
【数５５】
2019184773.tif 000057
【０１４６】
そこで、μｋｎを音源位置事後確率と呼ぶ。これに対し、βｋｎ（（６６）式）は、空間共分散行列Ｒｎ（１），・・・，Ｒｎ（Ｆ）が与えられる「前」の音源位置確率であるため、音源位置事前確率と呼ぶ。
【０１４７】
【数５６】
2019184773.tif 000058
【０１４８】
（５６）式の処理は音源存在事後確率更新部１２において、（５７）式の処理は音源位置事後確率更新部１４において、（５９）式の処理は音源存在事前確率更新部１５において、（６０）式の処理は音源位置事前確率更新部１６において、（６１）式の処理は空間共分散行列更新部１７において、（６２）式の処理はパワーパラメータ更新部１８において、それぞれ行われる。
【０１４９】
ここで、上述のパラメータΘの更新則（５９）式〜（６２）式の導出について説明する。まず、（５２）式の補助関数は次の（６７）式および（６８）式のように計算できる。ここで、ＣはパラメータΘに依らない定数である。
【０１５０】
【数５７】
2019184773.tif 000059
【０１５１】
音源存在事前確率αｎ（ｆ）の更新則（５９）式を導出するために、拘束条件（６）式に注意して、ξをラグランジュの未定乗数として、（６９）式をαｎ（ｆ）で微分したものを０と置くと、（７０）式となる。
【０１５２】
【数５８】
2019184773.tif 000060
【０１５３】
【数５９】
2019184773.tif 000061
【０１５４】
（７０）式をαｎ（ｆ）について解くと、（７１）式となる。
【０１５５】
【数６０】
2019184773.tif 000062
【０１５６】
（７１）式に含まれるラグランジュの未定乗数ξの値を決定するために、（７１）式を拘束条件（６）式に代入すると、（７２）式〜（７４）式となる。
【０１５７】
【数６１】
2019184773.tif 000063
【０１５８】
よって、ξ＝Ｔであるから、音源存在事前確率αｎ（ｆ）の更新則（５９）式が得られる。音源位置事前確率βｋｎの更新則（６０）式も同様にして導出できるから、説明を省略する。
【０１５９】
空間共分散行列Ｒｎ（ｆ）の更新則（６１）式を導出するために、（６８）式をＲｎ（ｆ）で微分したものを０と置くと、（７５）式となる。
【０１６０】
【数６２】
2019184773.tif 000064
【０１６１】
上式の両辺に対し、左右からそれぞれＲｎ（ｆ）を掛けると、（７６）式となる。これをＲｎ（ｆ）について解けば、空間共分散行列Ｒｎ（ｆ）の更新則（６１）式が得られる。
【０１６２】
【数６３】
2019184773.tif 000065
【０１６３】
パワーパラメータｖｎ（ｔ，ｆ）の更新則（６２）式を導出するために、（６８）式をｖｎ（ｔ，ｆ）で微分したものを０と置くと、（７７）式となる。
【０１６４】
【数６４】
2019184773.tif 000066
【０１６５】
これをｖｎ（ｔ，ｆ）について解けば、パワーパラメータｖｎ（ｔ，ｆ）の更新則（６２）式が得られる。以上で、上述のパラメータΘの更新則（５９）式〜（６２）式が導出できた。
【０１６６】
本第１の実施形態においては、複素ガウス分布のパラメータである空間共分散行列Ｒｎ（ｆ）の事前分布が、複素逆ウィシャート分布に基づく事前分布であるというモデル化に基づいている。このように、複素ガウス分布と複素逆ウィシャート分布とを組み合わせて用いることにより、補助関数Ｑ（Θ，Φ）が、その空間共分散行列Ｒｎ（ｆ）に関する微分を０と置いた式がＲｎ（ｆ）について解ける（上述）ような形になる。これは、複素逆ウィシャート分布が複素ガウス分布の共役事前分布であることに起因する。共役事前分布については、参考文献２「C.M.  Bishop,“Pattern  Recognition  and  Machine  Learning”,  Springer,  2006.」を参照されたい。
【０１６７】
［第１の実施形態の効果］
このように、本第１の実施形態では、空間共分散行列の各信号源に対する事前分布を、空間共分散行列の複数の各信号源位置候補に対する事前分布の線型結合である混合分布によりモデル化するときの混合重みであり、信号源ごとの各信号源位置候補から信号が到来する確率である、信号源位置事前確率を推定する。具体的には、本第１の実施形態では、空間共分散行列の各信号源に対する事前分布を（１７）式のようにモデル化している。そして、本第１の実施形態では、未知の確率である音源位置事前確率βｋｎによる重み付き和に基づくことによって、各音源の音源位置が未知の場合でも、空間共分散行列の事前分布を設計することができる。したがって、本第１の実施形態では、各音源信号に対する音源位置が未知である場合にも、空間共分散行列の事前分布に基づいて音源分離を行うことができる。
【０１６８】
また、本第１の実施形態では、（５２）式に示すように、対数ｌｎの中にｋに関する和がない補助関数を用いるため、補助関数の各パラメータに関する微分が単純になり、パラメータの更新演算が煩雑ではなくなる。
【０１６９】
また、本第１の実施形態では、空間共分散行列の事前分布が、複素逆ウィシャート分布に基づく事前分布であるというモデル化に基づいている。このように、本第１の実施形態では、複素ガウス分布と複素逆ウィシャート分布とを組み合わせて用いることにより、補助関数Ｑ（Θ，Φ）が、その空間共分散行列Ｒｎ（ｆ）に関する微分を０と置いた式がＲｎ（ｆ）について解ける。
【０１７０】
［第１の実施形態の変形例１］
本第１の実施形態では、観測データとして観測信号ベクトルｙ（ｔ，ｆ）を用いたが、他の特徴ベクトルまたは特徴量を観測データとして用いてもよい。例えば、観測信号ベクトルｙ（ｔ，ｆ）に基づいて、（７８）式および（７９）式で定義される特徴ベクトルｚ（ｔ，ｆ）を用いてもよい。
【０１７１】
【数６５】
2019184773.tif 000067
【数６６】
2019184773.tif 000068
【０１７２】
また、観測データとして、マイクロホン間の位相差、振幅比や、音源信号の到来時間差、到来方向などの特徴量を用いてもよい。
【０１７３】
また、本第１の実施形態では、特徴ベクトルである観測信号ベクトルに当てはめる混合モデルとして、混合複素ガウス分布を用いたが、利用される特徴ベクトルに応じて、様々な混合モデル（例えば、混合ガウス分布、混合ラプラス分布、混合複素ワトソン分布、混合複素ビンガム分布、混合複素角度中心ガウス分布、フォンミーゼス分布など）を用いることができる。また、混合モデルに限らず、複素ガウス分布などのモデルを特徴ベクトルである観測信号ベクトルに当てはめてもよい。
【０１７４】
また、本第１の実施形態では、空間共分散行列の事前分布を混合複素逆ウィシャート分布によりモデル化したが、混合複素ウィシャート分布などの他のモデルによりモデル化してもよい。
【０１７５】
また、本第１の実施形態では、モデルを観測データに当てはめるために、パラメータΘの事後確率を最大化する方法を採用したが、他の方法によりモデルを観測データに当てはめてもよい。
【０１７６】
また、本第１の実施形態では、補助関数法により最適化を行ったが、勾配法などの他の方法により最適化を行ってもよい。その場合、音源存在事後確率更新部１２および音源位置事後確率更新部１４は必須ではない。
【０１７７】
［第１の実施形態の変形例２］
真の音源数Ｎ´が未知の場合に、真の音源数Ｎ´の推定や音源分離を行う第１の実施形態の変形例２について説明する。本変形例では、仮定された音源数ＮはＮ≧Ｎ´となるように十分大きく設定されているものとする。例えば、想定される音源数が高々６個であると分かっているような場合には、仮定された音源数はＮ＝６と設定すればよい。なお、実際の音源数は４個である場合には、Ｎ´＝４となる。
【０１７８】
推定部１０は、各ｎ（ｎは１以上Ｎ以下の整数）に対し、音源位置事前確率更新部１６からの音源位置事前確率βｋｎが最大となるｋに対応する音源位置候補を音源位置の推定値とする。そして、信号分析装置１は、このようにして得られたＮ個の音源位置を、階層クラスタリングなどによりクラスタリングし、得られたクラスタの個数を、実際の音源数Ｎ´の推定値＾Ｎ´とする。
【０１７９】
クラスタリングにより得られた＾Ｎ´個の各クラスタは、＾Ｎ´個の実際の音源に対応するものとみなされる。従ってこのクラスタリングにより、Ｎ個の仮定する各音源ｎが、＾Ｎ´個の実際の音源のうちのどれに対応するか、が分かる。音源分離を行う場合には、この対応関係を利用して、推定部１０が以降の処理も行う。
【０１８０】
推定部１０は更に、得られた＾Ｎ´個の各クラスタｎ´（ｎ´は１以上＾Ｎ´以下の整数であるクラスタのインデックス）に対し、Ｎ個の仮定する音源の音源存在事後確率λｎ（ｔ，ｆ）のうち該クラスタに対応するものを加算することにより、ｎ´番目の実際の音源の音源存在事後確率λ´ｎ´（ｔ，ｆ）を計算する。推定部１０は更に、式（８）と同様に、各時間周波数点（ｔ，ｆ）に対し、実際の音源の音源存在事後確率λ´ｎ´（ｔ，ｆ）が最大となる番号ｎ´に対応する実際の音源からの信号が（ｔ，ｆ）にて鳴っていると判定する。推定部１０は更に、（４）式と同様に、実際の音源の音源信号成分の推定値＾ｘ´ｎ´（ｔ，ｆ）を、（ｔ，ｆ）においてｎ´番目の実際の音源が鳴っていると判定された場合にはｙ（ｔ，ｆ）とし、そうでないと判定された場合には０とすることにより、音源分離を行う。
【０１８１】
［第１の実施形態の変形例３］
本第１の実施形態は、音信号に限らず、他の信号（脳波、脳磁図、無線信号など）に対して適用してもよい。本第１の実施形態における観測信号は、複数のマイクロホン（マイクロホンアレイ）により取得された観測信号に限らず、脳波計、脳磁計、アンテナアレイなどの他のセンサアレイ（複数のセンサ）により取得された、空間上の位置から時系列として発生する信号からなる観測信号であってもよい。
【０１８２】
［第１の実施形態の変形例４］
観測信号ベクトルｙ（ｔ，ｆ）の確率分布を次の（８０）式の複素ガウス分布によりモデル化する例を第１の実施形態の変形例４として説明する。この場合のパラメータΘの更新則は、第１の実施形態の（５６）、（５７）、（５９）、（６０）、（６１）、（６２）式に代えて、（８１）式〜（８６）式のようになる。
【０１８３】
【数６７】
2019184773.tif 000069
【０１８４】
図３および図４を用いて、第１の実施形態の変形例４の構成と処理について説明する。図３は、第１の実施形態の変形例４に係る信号分析装置の構成の一例を示す図である。図４は、第１の実施形態の変形例４に係る信号分析処理の処理手順の一例を示すフローチャートである。
【０１８５】
図３に示すように、本第１の実施形態の変形例４に係る信号分析装置２０１は、観測信号ベクトル作成部１１、初期化部（図示しない）、記憶部１３、推定部２１０、パワーパラメータ更新部２１８、収束判定部（図示しない）を有する。推定部２１０は、音源位置事後確率更新部２１２、音源信号事後確率更新部２１３、音源位置事前確率更新部２１４、空間共分散行列更新部２１７を有する。
【０１８６】
観測信号ベクトル作成部１１は、第１の実施形態と同様に、観測信号ベクトルｙ（ｔ，ｆ）を（１）式により作成する（ステップＳ２１〜ステップＳ２３）。
【０１８７】
初期化部は、音源位置事前確率βｋｎと、空間共分散行列Ｒｎ（ｆ）と、パワーパラメータｖｎ（ｔ，ｆ）と、の推定値の初期値を計算する（ステップＳ２４）。ただし、ｎ＝１，・・・，Ｎは音源のインデックス、ｋ＝１，・・・，Ｋは音源位置候補のインデックスである。例えば初期化部は、乱数に基づいてこれらの初期値を計算する。また、初期化部は、ｎを初期化する（ステップＳ２５）。
【０１８８】
なお、記憶部１３は、各音源位置候補ｋ、各周波数ビンｆに対する空間共分散行列の事前分布のパラメータであるΨｋ（ｆ）およびνｋ（ｆ）を記憶する。
【０１８９】
続いて、信号分析装置２０１は、ｎに１を加算して（ステップＳ２６）、ステップＳ２７〜ステップＳ３１の処理を行う。
【０１９０】
音源位置事後確率更新部２１２は、記憶部１３からの事前分布のパラメータであるΨｋ（ｆ）およびνｋ（ｆ）と、音源位置事前確率更新部２１４からの音源位置事前確率（ただし例外として、音源位置事後確率更新部２１２における最初の処理の際には、初期化部からの音源位置事前確率の初期値）βｋｎと、空間共分散行列更新部２１７からの空間共分散行列（ただし例外として、音源位置事後確率更新部２１２における最初の処理の際には、初期化部からの空間共分散行列の初期値）Ｒｎ（ｆ）と、を受け取って、音源位置事後確率μｋｎを（８１）式により更新する（ステップＳ２７）。
【０１９１】
【数６８】
2019184773.tif 000070
【０１９２】
音源信号事後確率更新部２１３は、観測信号ベクトル作成部１１からの観測信号ベクトルｙ（ｔ，ｆ）と、パワーパラメータ更新部２１８からのパワーパラメータ（ただし例外として、音源信号事後確率更新部２１３における最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖｎ（ｔ，ｆ）と、空間共分散行列更新部２１７からの空間共分散行列（ただし例外として、音源信号事後確率更新部２１３における最初の処理の際には、初期化部からの空間共分散行列の初期値）Ｒｎ（ｆ）と、を受け取って、音源信号成分ｘｎ（ｔ，ｆ）の事後確率の平均ξｎ（ｔ，ｆ）および共分散行列Σｎ（ｔ，ｆ）を、（８２）式および（８３）式により更新する（ステップＳ２８）。
【０１９３】
【数６９】
2019184773.tif 000071
【数７０】
2019184773.tif 000072
【０１９４】
音源位置事前確率更新部２１４は、音源位置事後確率更新部２１２からの音源位置事後確率μｋｎを受け取って、音源位置事前確率βｋｎを（８４）式により更新する（ステップＳ２９）。
【０１９５】
【数７１】
2019184773.tif 000073
【０１９６】
空間共分散行列更新部２１７は、記憶部１３からの事前分布のパラメータであるΨｋ（ｆ）およびνｋ（ｆ）と、音源位置事後確率更新部２１２からの音源位置事後確率μｋｎと、音源信号事後確率更新部２１３からの事後確率の平均ξｎ（ｔ，ｆ）および共分散行列Σｎ（ｔ，ｆ）と、パワーパラメータ更新部２１８からのパワーパラメータ（ただし例外として、空間共分散行列更新部２１７における最初の処理の際には、初期化部からのパワーパラメータの初期値）ｖｎ（ｔ，ｆ）と、を受け取って、空間共分散行列Ｒｎ（ｆ）を（８５）式により更新する（ステップＳ３０）。
【０１９７】
【数７２】
2019184773.tif 000074
【０１９８】
パワーパラメータ更新部２１８は、空間共分散行列更新部２１７からの空間共分散行列Ｒｎ（ｆ）と、音源信号事後確率更新部２１３からの事後確率の平均ξｎ（ｔ，ｆ）および共分散行列Σｎ（ｔ，ｆ）と、を受け取って、パワーパラメータｖｎ（ｔ，ｆ）を（８６）式により更新する（ステップＳ３１）。
【０１９９】
【数７３】
2019184773.tif 000075
【０２００】
そして、信号分析装置２０１は、ｎ＝Ｎか否かを判定する（ステップＳ３２）。信号分析装置２０１は、ｎ＝Ｎでないと判定した場合（ステップＳ３２：Ｎｏ）、ステップＳ２６に戻る。これに対し、信号分析装置２０１は、ｎ＝Ｎであると判定した場合（ステップＳ３２：Ｙｅｓ）、収束判定部による判定処理に進む。
【０２０１】
収束判定部は、収束したかどうかの判定を行う（ステップＳ３３）。信号分析装置２０１は、収束していないと収束判定部が判定した場合（ステップＳ３３：Ｎｏ）、ステップＳ２５に戻って、処理を継続する。一方、収束したと収束判定部が判定した場合（ステップＳ３３：Ｙｅｓ）、音源信号事後確率更新部２１３は、事後確率の平均ξｎ（ｔ，ｆ）を、音源信号成分ｘｎ（ｔ，ｆ）の推定値＾ｘｎ（ｔ，ｆ）として出力し（ステップＳ３４）、信号分析装置２０１での処理が終了する。
【０２０２】
［第１の実施形態の変形例５］
第１の実施形態では、空間共分散行列により音源信号の空間的特性をモデル化したが、他のパラメータにより音源信号の空間的特性をモデル化してもよい。音源信号の空間的特性をモデル化するパラメータを、ここでは空間パラメータと呼ぶ。
【０２０３】
例えば、空間パラメータとしてステアリングベクトルを用い、これにより音源信号の空間的特性をモデル化してもよい。この場合、観測信号ベクトルｙ（ｔ，ｆ）の確率分布は、例えば次の（８７）式の複素ガウス分布によりモデル化できる。
【０２０４】
【数７４】
2019184773.tif 000076
【０２０５】
ここで、ｈｎ（ｆ）は、音源信号ｎの空間的特性をモデル化する空間パラメータであるステアリングベクトルであり、σ１２は正則化のための正数である。この場合、ｈｎ（ｆ）の事前分布は次の（８８）式で与えられる。但し、（８８）式における「ｐ」は、複素ガウス分布「ｐＧ」を表す。
【０２０６】
【数７５】
2019184773.tif 000077
【０２０７】
ここで、ｇｋ（ｆ）とσ２２はハイパーパラメータである。ｇｋ（ｆ）はｋ番目の音源位置候補に対するステアリングベクトルであり、σ２２は正則化のための正数である。以上のモデル化に基づいて、第１の実施形態と同様にパラメータΘを推定すればよい。
【０２０８】
［システム構成等］
また、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示のように構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部又は一部を、各種の負荷や使用状況等に応じて、任意の単位で機能的又は物理的に分散・統合して構成することができる。さらに、各装置にて行われる各処理機能は、その全部又は任意の一部が、ＣＰＵおよび当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。
【０２０９】
また、本実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的に行うこともでき、あるいは、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。この他、上記文書中や図面中で示した処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて任意に変更することができる。すなわち、上記学習方法および音声認識方法において説明した処理は、記載の順にしたがって時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。
【０２１０】
［プログラム］
図５は、プログラムが実行されることにより、信号分析装置１，２０１が実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。これらの各部は、バス１０８０によって接続される。
【０２１１】
メモリ１０１０は、ＲＯＭ１０１１およびＲＡＭ１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ（Basic  Input  Output  System）等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０９０に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１１００に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１１００に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１１１０、キーボード１１２０に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１１３０に接続される。
【０２１２】
ハードディスクドライブ１０９０は、例えば、ＯＳ（Operating  System）１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわち、信号分析装置１，２０１の各処理を規定するプログラムは、コンピュータ１０００により実行可能なコードが記述されたプログラムモジュール１０９３として実装される。プログラムモジュール１０９３は、例えばハードディスクドライブ１０９０に記憶される。例えば、信号分析装置１，２０１における機能構成と同様の処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０９０に記憶される。なお、ハードディスクドライブ１０９０は、ＳＳＤ（Solid  State  Drive）により代替されてもよい。
【０２１３】
また、上述した実施形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０９０に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０９０に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。
【０２１４】
なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０９０に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１１００等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３およびプログラムデータ１０９４は、ネットワーク（ＬＡＮ（Local  Area  Network）、ＷＡＮ（Wide  Area  Network）等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３およびプログラムデータ１０９４は、他のコンピュータから、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。
【０２１５】
以上、本発明者によってなされた発明を適用した実施形態について説明したが、本実施形態による本発明の開示の一部をなす記述および図面により本発明は限定されることはない。すなわち、本実施形態に基づいて当業者等によりなされる他の実施形態、実施例および運用技術等はすべて本発明の範疇に含まれる。
【符号の説明】
【０２１６】
１，２０１，１Ｐ  信号分析装置
１０  推定部
１１，１１Ｐ  観測信号ベクトル作成部
１２，１２Ｐ  音源存在事後確率更新部
１３，１３Ｐ  記憶部
１４，２１２  音源位置事後確率更新部
１４Ｐ  音源存在事前確率更新部
１５  音源存在事前確率更新部
１６，２１４  音源位置事前確率更新部
１７，２１７，１５Ｐ  空間共分散行列更新部
１８，２１８，１６Ｐ  パワーパラメータ更新部
１９，１７Ｐ  音源信号成分推定部
２１３  音源信号事後確率更新部
【図１】
2019184773.tif 000078
【図２】
2019184773.tif 000079
【図３】
2019184773.tif 000080
【図４】
2019184773.tif 000081
【図５】
2019184773.tif 000082
【図６】
2019184773.tif 000083
【図７】
2019184773.tif 000084