JP2012027196A - 信号分析装置、方法、及びプログラム - Google Patents

Abstract

【課題】時系列信号から、頻出のパターンとパターンの出現の時間的な関連とを精度良く抽出することができるようにする。
【解決手段】信号分析部３によって、パラメータＨ_ω,φ,dとパラメータＰ_φ,t,dとパラメータＵ_d,tとの積を、（ｄ，φ）の全ての組み合わせについて足し合わせたモデルについて、時間周波数解析部１によって出力された二次元配列Ｙとモデルとの距離をβダイバージェンスで表わした目的関数の値が小さくなるように、パラメータＨ_ω,φ,d、パラメータＰ_φ,t,d、及びパラメータＵ_d,tの各々の値を更新する。収束判定部３４は、繰り返し回数Ｓに到達するまで、補助変数更新部３１による計算、スケジューリング係数更新部３２による計算、及びパラメータ更新部３３による更新を繰り返し行う。
【選択図】図１

Description

本発明は、信号分析装置、方法、及びプログラムに係り、特に、時系列信号から、信号パラメータを分析する信号分析装置、方法、及びプログラムに関する。

従来、信号の分析に関しては、信号の時間周波数表現に基づく各種の分析方法が広く知られている。信号内に繰り返して現れるパターン等のような高度な構造を抽出する分析手法として非負値行列分解（ＮＭＦ：Ｎｏｎ‐ｎｅｇａｔｉｖｅ Ｍａｔｒｉｘ Ｆａｃｔｏｒｉｚａｔｉｏｎ）が知られている（例えば、非特許文献１）。

ＮＭＦは、自動採譜や、モノラル混合信号からの音源の分離に適用されている（例えば、非特許文献２、非特許文献３）。

また、ポリフォニックオーディオデータを確率モデルで表現することにより、信号分離を行う技術も知られている（非特許文献４）。これは、Factorial Scaled Hidden markov Modelを利用したものであり、距離尺度を板倉齊藤距離とした場合のＮＭＦモデルを、時間的な関連が表現できるように拡張したものに相当する。

D. D. Lee, H. S. Seung, "Learning the part of objects by non-negative matrix factorization," Nature Vol.401, pp. 788-791, 1999. P. Smaragdis, J. C. Brown, "Non-Negative Matrix Factorization for Music Transcription,"ln Proc. 2003 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA2003), pp. 177-180,2003. T. Virtanen, "Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria," IEEE Transactions on Audio, Speech, and Language Processing, vol. 15, pp. 1066-1074, 2007. A. Ozerov, C. Fevotte and M. Charbit, "Factorial Scaled Hidden Markov Model for Polyphonic Audio Representation and Source Separation", IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, pp.121-124, 2009.

従来のＮＭＦによる信号の分析は、信号内に頻出するパターンを表出するものであるが、それらのパターンの出現に関する時間的な関連を表現することができない、という問題がある。例えば、ピアノのような打鍵楽器は、音の立ち上がりにおける、アタックといわれる広帯域なエネルギーを持つ状態から、サステイン、ディケイ、リリースなどのいくつかの状態を経て消音する。また、歌声や楽器のビブラートでは、一音の中でも基本周波数が変化する。これらの楽器音の詳細なスペクトルの変化は、単にＮＭＦの分解性能を上げるだけでは表現することができない。

また、上記の非特許文献４に記載の技術によれば、ＮＭＦを、ある特定の条件の場合について時間的な関連を表現したモデルに拡張することができ、スペクトルの時間的な変化を分析することができる。板倉齊藤距離は、β-divergenceと呼ばれる距離尺度の特殊な場合（β＝０）である。

しかしながら、ＮＭＦにおける板倉齊藤距離は、局所解に捕まりやすい距離尺度として知られている。そのため、音源の種類によっては精度の良い音源分離を行うことができない、という問題がある。また、板倉齊藤距離（β＝０）以外の距離尺度に基づいて、スペクトルの時間的な変化を分析できるようにＮＭＦを拡張する方法は知られていなかった。

本発明は、上記の課題を解決するためになされたもので、時系列信号から、頻出のパターンとパターンの出現の時間的な関連とを精度良く抽出することができる信号分析装置、方法、及びプログラムを提供することを目的とする。

上記の目的を達成するために本発明に係る信号分析装置は、時系列信号を入力として、観測時間周波数成分Ｙ_ω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙを出力する時間周波数分解手段と、複数の基底ｄの各々における全ての基底の状態φを表わすパラメータＨ_ω,φ,dを要素にもつ三次元配列Ｈ、前記複数の基底ｄの各々に対して各時刻ｔに前記基底の状態が前記三次元配列Ｈの要素のうちの何れであるかを表わすパラメータＰ_φ,t,dを要素にもつ三次元配列Ｐ、及び前記複数の基底ｄの各々に対する各時刻ｔのゲインを表すパラメータＵ_d,tを要素にもつ二次元配列Ｕの各々の初期値を設定するパラメータ初期値設定手段と、全てのパラメータＨ_ω,φ,dが非負値であり、全てのパラメータＵ_d,tが非負値であり、かつ、前記三次元配列Ｐが、前記複数の基底ｄの各々について予め許可された前記基底の状態の状態遷移に基づくものである、という条件の下で、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を、（ｄ，φ）の全ての組み合わせについて足し合わせたモデルについて、前記時間周波数分解手段によって出力された二次元配列Ｙと前記モデルとの距離をβダイバージェンスで表わした目的関数の値が小さくなるように、前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの各々の値を更新する信号分析手段と、を含んで構成されている。

本発明に係る信号分析装置において、前記信号分析手段は、前記時間周波数分解手段によって出力された二次元配列Ｙの観測時間周波数成分Ｙ_ω,tと、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を（ｄ，φ）の全ての組み合わせについて足し合わせたものとの差分に、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を全ての基底の状態φについて合計したものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算した商μ_{ω，ｔ，ｄ}を、乗算したものに、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて合計したものを加算して、補助変数＾Ｙ_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算する補助変数計算手段と、全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算するスケジューリング係数計算手段と、前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記パラメータＰ_φ,t,dと前記補助変数＾Ｙ_{ω，ｔ，ｄ}との積から、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を減算したものを２乗して、前記商μ_{ω，ｔ，ｄ}で除算した商を、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、前記補助変数＾Ｙ_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値の２乗と前記パラメータＵ_d,tとを乗算した積を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積の２乗を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたもので除算して、前記パラメータＨ_ω,φ,dの更新値を（ω，φ，ｄ）の全ての組み合わせについて算出し、前記パラメータＨ_ω,φ,dの更新値と前記補助変数＾Ｙ_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値の２乗とを乗算した積を、前記商μ_{ω，ｔ，ｄ}で除算した商を、（ω，ｔ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積の２乗を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたもので除算して、前記パラメータＵ_d,tの更新値を（ｄ，ｔ）の全ての組み合わせについて算出するパラメータ更新手段と、前記繰り返し回数Ｓに到達するまで、前記補助変数計算手段による計算、前記スケジューリング係数計算手段による計算、及び前記パラメータ更新手段による更新を繰り返し行う収束判定手段と、を備えるようにすることができる。

本発明に係る信号分析装置において、前記信号分析手段は、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて足し合わせたものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算して、補助変数η_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算する補助変数計算手段と、全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算するスケジューリング係数計算手段と、前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記補助変数η_{ω，ｔ，ｄ}と、前記パラメータＰ_φ,t,dと、前記観測時間周波数成分Ｙ_ω,tと、前記補助変数η_{ω，ｔ，ｄ}及び前記観測時間周波数成分Ｙ_ω,tの積を前記パラメータＨ_ω,φ,d及び前記パラメータＵ_d,tの積で除算した商について対数をとったものと、を乗算した積を、前記観測時間周波数成分Ｙ_ω,tを前記基底の数Ｄと前記状態数Φとの積で除算した商で減算し、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を加算したものを、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、前記補助変数η_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積を各時刻ｔについて足し合わせたもので除算して、前記パラメータＨ_ω,φ,dの更新値を（ω，φ，ｄ）の全ての組み合わせについて算出し、前記補助変数η_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を（ω，φ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積を（ω，φ）の全ての組み合わせについて足し合わせたもので除算して、前記パラメータＵ_d,tの更新値を（ｄ，ｔ）の全ての組み合わせについて算出するパラメータ更新手段と、前記繰り返し回数Ｓに到達するまで、前記補助変数計算手段による計算、前記スケジューリング係数計算手段による計算、及び前記パラメータ更新手段による更新を繰り返し行う収束判定手段と、を備えるようにすることができる。

本発明に係る信号分析装置において、前記信号分析手段は、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を（ｄ，φ）の全ての組み合わせについて足し合わせて、補助変数Ｚ_ω，ｔを（ω，ｔ）の全ての組み合わせについて計算すると共に、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて足し合わせたものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算して、補助変数η_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算する補助変数計算手段と、全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算するスケジューリング係数計算手段と、前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記パラメータＰ_φ,t,d、前記補助変数η_{ω，ｔ，ｄ}の２乗、及び前記観測時間周波数成分Ｙ_ω,tの積を、前記パラメータＨ_ω,φ,d及び前記パラメータＵ_d,tの積で除算した商と、前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの積を前記補助変数Ｚ_ω，ｔで除算した商と、前記補助変数Ｚ_ω，ｔを前記観測時間周波数成分Ｙ_ω,tで除算した商について対数をとったものから２を減算して前記基底の数Ｄで除算した商と、を加算したものを、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、前記パラメータＰ_φ,t,dの更新値と前記補助変数η_{ω，ｔ，ｄ}の２乗と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を、前記パラメータＵ_d,tで除算した商を、各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積を、前記補助変数Ｚ_ω，ｔで除算した商を、各時刻ｔについて足し合わせたもので除算した商の平方根を、前記パラメータＨ_ω,φ,dの更新値として（ω，φ，ｄ）の全ての組み合わせについて算出し、前記パラメータＰ_φ,t,dの更新値と前記補助変数η_{ω，ｔ，ｄ}の２乗と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を、前記パラメータＨ_ω,φ,dの更新値で除算した商を、（ω，φ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積を前記補助変数Ｚ_ω，ｔで除算した商を、（ω，φ）の全ての組み合わせについて足し合わせたもので除算した商の平方根を、前記パラメータＵ_d,tの更新値として（ｄ，ｔ）の全ての組み合わせについて算出するパラメータ更新手段と、前記繰り返し回数Ｓに到達するまで、前記補助変数計算手段による計算、前記スケジューリング係数計算手段による計算、及び前記パラメータ更新手段による更新を繰り返し行う収束判定手段と、を備えるようにすることができる。

本発明に係る信号分析方法は、時系列信号を入力として、観測時間周波数成分Ｙ_ω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙを出力し、複数の基底ｄの各々における全ての基底の状態φを表わすパラメータＨ_ω,φ,dを要素にもつ三次元配列Ｈ、前記複数の基底ｄの各々に対して各時刻ｔに前記基底の状態が前記三次元配列Ｈの要素のうちの何れであるかを表わすパラメータＰ_φ,t,dを要素にもつ三次元配列Ｐ、及び前記複数の基底ｄの各々に対する各時刻ｔのゲインを表すパラメータＵ_d,tを要素にもつ二次元配列Ｕの各々の初期値を設定し、全てのパラメータＨ_ω,φ,dが非負値であり、全てのパラメータＵ_d,tが非負値であり、かつ、前記三次元配列Ｐが、前記複数の基底ｄの各々について予め許可された前記基底の状態の状態遷移に基づくものである、という条件の下で、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を、（ｄ，φ）の全ての組み合わせについて足し合わせたモデルについて、前記時間周波数分解手段によって出力された二次元配列Ｙと前記モデルとの距離をβダイバージェンスで表わした目的関数の値が小さくなるように、前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの各々の値を更新することを特徴とする。

本発明に係る信号分析方法において、前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの各々の値を更新することは、前記出力された二次元配列Ｙの観測時間周波数成分Ｙ_ω,tと、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を（ｄ，φ）の全ての組み合わせについて足し合わせたものとの差分に、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を全ての基底の状態φについて合計したものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算した商μ_{ω，ｔ，ｄ}を、乗算したものに、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて合計したものを加算して、補助変数＾Ｙ_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算し、全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記パラメータＰ_φ,t,dと前記補助変数＾Ｙ_{ω，ｔ，ｄ}との積から、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を減算したものを２乗して、前記商μ_{ω，ｔ，ｄ}で除算した商を、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、前記補助変数＾Ｙ_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値の２乗と前記パラメータＵ_d,tとを乗算した積を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積の２乗を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたもので除算して、前記パラメータＨ_ω,φ,dの更新値を（ω，φ，ｄ）の全ての組み合わせについて算出し、前記パラメータＨ_ω,φ,dの更新値と前記補助変数＾Ｙ_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値の２乗とを乗算した積を、前記商μ_{ω，ｔ，ｄ}で除算した商を、（ω，ｔ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積の２乗を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたもので除算して、前記パラメータＵ_d,tの更新値を（ｄ，ｔ）の全ての組み合わせについて算出し、前記繰り返し回数Ｓに到達するまで、前記補助変数＾Ｙ_{ω，ｔ，ｄ}の計算、前記スケジューリング係数ｋ_φ ^（ｓ）の計算、前記パラメータＨ_ω,φ,dの更新値の算出、前記パラメータＰ_φ,t,dの更新値の算出、及び前記パラメータＵ_d,tの更新値の算出を繰り返し行うことを特徴とすることができる。

本発明に係る信号分析方法は、前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの各々の値を更新することは、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて足し合わせたものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算して、補助変数η_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算し、全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記補助変数η_{ω，ｔ，ｄ}と、前記パラメータＰ_φ,t,dと、前記観測時間周波数成分Ｙ_ω,tと、前記補助変数η_{ω，ｔ，ｄ}及び前記観測時間周波数成分Ｙ_ω,tの積を前記パラメータＨ_ω,φ,d及び前記パラメータＵ_d,tの積で除算した商について対数をとったものと、を乗算した積を、前記観測時間周波数成分Ｙ_ω,tを前記基底の数Ｄと前記状態数Φとの積で除算した商で減算し、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を加算したものを、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、前記補助変数η_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積を各時刻ｔについて足し合わせたもので除算して、前記パラメータＨ_ω,φ,dの更新値を（ω，φ，ｄ）の全ての組み合わせについて算出し、前記補助変数η_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を（ω，φ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積を（ω，φ）の全ての組み合わせについて足し合わせたもので除算して、前記パラメータＵ_d,tの更新値を（ｄ，ｔ）の全ての組み合わせについて算出し、前記繰り返し回数Ｓに到達するまで、前記補助変数η_{ω，ｔ，ｄ}の計算、前記スケジューリング係数ｋ_φ ^（ｓ）の計算、前記パラメータＨ_ω,φ,dの更新値の算出、前記パラメータＰ_φ,t,dの更新値の算出、及び前記パラメータＵ_d,tの更新値の算出を繰り返し行うことを特徴とすることができる。

本発明に係る信号分析方法において、前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの各々の値を更新することは、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を（ｄ，φ）の全ての組み合わせについて足し合わせて、補助変数Ｚ_ω，ｔを（ω，ｔ）の全ての組み合わせについて計算すると共に、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて足し合わせたものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算して、補助変数η_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算し、全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記パラメータＰ_φ,t,d、前記補助変数η_{ω，ｔ，ｄ}の２乗、及び前記観測時間周波数成分Ｙ_ω,tの積を、前記パラメータＨ_ω,φ,d及び前記パラメータＵ_d,tの積で除算した商と、前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの積を前記補助変数Ｚ_ω，ｔで除算した商と、前記補助変数Ｚ_ω，ｔを前記観測時間周波数成分Ｙ_ω,tで除算した商について対数をとったものから２を減算して前記基底の数Ｄで除算した商と、を加算したものを、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、前記パラメータＰ_φ,t,dの更新値と前記補助変数η_{ω，ｔ，ｄ}の２乗と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を、前記パラメータＵ_d,tで除算した商を、各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積を、前記補助変数Ｚ_ω，ｔで除算した商を、各時刻ｔについて足し合わせたもので除算した商の平方根を、前記パラメータＨ_ω,φ,dの更新値として（ω，φ，ｄ）の全ての組み合わせについて算出し、前記パラメータＰ_φ,t,dの更新値と前記補助変数η_{ω，ｔ，ｄ}の２乗と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を、前記パラメータＨ_ω,φ,dの更新値で除算した商を、（ω，φ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積を前記補助変数Ｚ_ω，ｔで除算した商を、（ω，φ）の全ての組み合わせについて足し合わせたもので除算した商の平方根を、前記パラメータＵ_d,tの更新値として（ｄ，ｔ）の全ての組み合わせについて算出し、前記繰り返し回数Ｓに到達するまで、前記補助変数Ｚ_ω，ｔ及び前記補助変数η_{ω，ｔ，ｄ}の計算、前記スケジューリング係数ｋ_φ ^（ｓ）の計算、前記パラメータＨ_ω,φ,dの更新値の算出、前記パラメータＰ_φ,t,dの更新値の算出、及び前記パラメータＵ_d,tの更新値の算出を繰り返し行うことを特徴とすることができる。

本発明に係るプログラムは、上記の信号分析装置の各手段としてコンピュータを機能させるためのプログラムである。

以上説明したように、本発明の信号分析装置、方法、及びプログラムによれば、時系列信号の観測時間周波数成分Ｙ_ω,tを要素とする二次元配列Ｙとの距離をβダイバージェンスで表わした目的関数の値が小さくなるように、各基底ｄにおける各基底の状態φを表わすパラメータＨ_ω,φ,d、各時刻ｔにおいて基底の状態がどの状態であるかを表わすパラメータＰ_φ,t,d、及び各基底ｄに対する各時刻ｔのゲインを表すパラメータＵ_d,tを更新することにより、時系列信号から、頻出のパターンとパターンの出現の時間的な関連とを精度良く抽出することができる、という効果が得られる。

本発明の第１の実施の形態に係る音響信号パラメータ分析装置の構成を示す概略図である。 提案モデルのイメージ図である。 本発明の第１の実施の形態に係る音響信号パラメータ分析装置における音響信号分析処理ルーチンの内容を示すフローチャートである。 （Ａ）分析結果として得られた基底の状態を示すグラフ、（Ｂ）入力された音響信号から計算された振幅スペクトログラムを示すグラフ、及び（Ｃ）分析結果として得られたアクティベーションを示すグラフである。 従来法と提案法について得られたＳＮＲとエラー率を示す図である。 （Ａ）３音源を混合した観測スペクトログラムを示すグラフ、（Ｂ）分離した１音源のスペクトログラムを示すグラフ、及び（Ｃ）基底の状態遷移を示すグラフである。 反復回数と目的関数の値との関係を示すグラフである。

以下、図面を参照して本発明の実施の形態を詳細に説明する。本発明で提案する手法では、入力信号を時間周波数にしたスペクトログラムを、限られた数のソースの和によって近似することにより、信号の分析を行う。ここで、一つのソースでは、非負値な基底スペクトルパラメータが、限られた数の状態を持ち、各時刻にその中の一つの状態が非負値なアクティビティ係数パラメータによって出現している、としてモデル化する。このとき、状態の出現には事前に拘束を加えることができる。入力信号のスペクトログラムとモデルとの間のβ-divergenceを目的関数とし、モデルパラメータの推定を、目的関数の値を最小化する問題として定式化し、反復的に一つの解を探索するアルゴリズムを導出する。

まず、本発明の原理について説明する。まず、信号のモデルについて説明する。

本発明で提案するモデルは、ＮＭＦによるスペクトログラムの分解表現を拡張したモデルである。ＮＭＦによるスペクトログラムの分解表現は、観測スペクトログラムを非負値行列

と見立て、

となるような基底

とアクティベーション

の２つの非負値行列を決めることによって得られる。ここで、ω（ω＝１，２，…，Ω），ｔ（ｔ＝１，２，・・・，Ｔ）はそれぞれ周波数と時刻に対応するインデックスである。また、Ｙ_ω，ｔは、行列Ｙの（ω，ｔ）成分であり、Ｒ^{≧０，Ω×Ｔ}は、各成分が０以上の実数からなるΩ×Ｔ行列の集合を表す。

このモデルは、観測スペクトログラムがＤ個のスペクトル（基底）とアクティベーションとの積によって表現される、というモデルになっている。スペクトルとアクティベーションのペアが、楽器音１音高となることを狙ったモデルである。

実際の楽器音（信号）は、時間に伴ってスペクトルが変化することを考慮すると、スペクトルが各時刻ｔにおいて一つの状態であるとみなし、時刻ｔにおけるｄ番目の基底の状態をＨ_ｄ，ｔ ^{（φｄ，ｔ）}で表すことによって、上記（１）式を、以下の（２）式に拡張する。

ここで、基底の状態も限られたΦ個しか存在せず、時間に伴ってそれらの基底の状態を遷移するものと仮定し、Ｈ＝（Ｈ_{ω，φ，ｄ}）_Ω×Φ×Ｄ∈Ｒ^{≧０，Ω×Φ×Ｄ}を用いてモデルを表現し直す。各ｄ＝１，・・・，Ｄに対するΦ個の基底の状態ｈ_φ ^（ｄ）＝（Ｈ_{１，φ，ｄ}，・・・，Ｈ_{Ω，φ，ｄ}）^Ｔ（φ＝１，・・・，Φ）に関して、許される状態遷移を隣接行列Ａ（ｄ）＝（Ａ_φ，φ ^（ｄ））_Φ×Φを用いて表す。状態ｉから状態ｊへの遷移が許される場合、Ａ_ｉ，ｊ ^（ｄ）＝１とし、遷移が許されない場合、Ａ_ｉ，ｊ ^（ｄ）＝０とする。例えば，以下の（３）式で表される隣接行列Ａは、状態のインデックス１→２、２→３、３→１の遷移のみが許されることを表している。

ここで、Ｐ＝（Ｐ_{φ，ｔ，ｄ}）_Φ×Ｔ×Ｄ∈Ｒ^{≧０，Φ×Ｔ×Ｄ}を考え、ｄ＝１，・・・，Ｄに対して各時刻ｔに基底の状態がｈ_φ ^（ｄ）のときにＰ_{φ，ｔ，ｄ}＝１とし、そうでないときにＰ_{φ，ｔ，ｄ}＝０とすることで、提案モデルは、図２に示すように、（Ｈ_{ω，φ，ｄ}）_Ω×Φ×Ｄと、（Ｐ_{φ，ｔ，ｄ}）_Φ×Ｔ×Ｄと、（Ｕ_t，ｄ）_T×Dと用いて表わされ、以下の（４）式で表現することができる。

パラメータＨ_ω,φ,dは、各基底ｄにおける全ての基底の状態φを表わす。パラメータＰ_φ,t,dは、各基底ｄに対して各時刻ｔに基底の状態が、パラメータＨ_ω,φ,dを要素にもつ三次元配列Ｈの要素のうちの何れであるかを表わす。パラメータＵ_d,tは、各基底ｄに対する各時刻ｔのゲインを表わす。

今、各ｄ＝１，・・・，Ｄに対するＡ^（ｄ）を固定した時にとりうるＰの集合をＧとする。すると、本発明により解くべき問題は、観測されたＹから、三次元行列Ｈ，二次元行列Ｕ，及び三次元行列Ｐ∈Ｇを推定することに帰着される。

つまり、行列Ｙに対して、上記（４）式を満たす各Ｈ_{ω，φ，ｄ}，Ｕ_ｄ，ｔ，Ｐ_{φ，ｔ，ｄ}を求める。具体的には、観測結果Ｙとモデル（上記（４）式の右辺）との距離（近さ）を目的関数とし、所定の条件の下で目的関数の値を最小化することにより、各Ｈ_{ω，φ，ｄ}，Ｕ_ｄ，ｔ，Ｐ_{φ，ｔ，ｄ}を求める。観測結果Ｙとの距離尺度としては、β-divergenceを用いる。具体的には、観測されたＹから、

という条件のもとで、以下の（５）式で表される目的関数Ｊ（θ）の値を最小化するθ＝｛Ｈ，Ｕ，Ｐ｝を求める。

ここで、関数Ｄ_βは以下の（６）式で定義される。

次に、音響信号の信号パラメータを分析して出力する音響信号パラメータ分析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。

図１に示すように、第１の実施の形態に係る音響信号パラメータ分析装置は、ＣＰＵと、ＲＡＭと、後述する音響信号分析処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。

音響信号パラメータ分析装置は、時間周波数解析部１と、パラメータ初期値設定部２と、信号分析部３と、記憶部４と、出力部５とを備えている。また、信号分析部３は、補助変数更新部３１と、スケジューリング係数更新部３２と、パラメータ更新部３３と、収束判定部３４とを備えている。

時間周波数解析部１は、時系列信号としての観測された音響信号ｙを入力として、時間周波数成分（観測時間周波数成分）Ｙ_ω,t（ω＝１,・・・,Ω,ｔ＝１,・・・,Ｔは、それぞれ周波数、時刻に対応するインデックスを示す。）を各（ω，ｔ）の要素にもつ二次元配列Ｙを計算し、信号分析部３に出力する。また、計算した時間周波数成分Ｙ_ω,tを、記憶部４に記憶しておく。より詳細には、時間周波数解析部１は、時系列信号ｙを入力として、短時間フーリエ変換（Short-Time Fourier Transform；ＳＴＦＴ）を用いて時間周波数解析を行うことにより時間周波数成分Ｙ_ω,tを計算し、時間周波数成分Ｙ_ω，ｔを格納した行列（振幅スペクトログラム）Ｙ＝（Ｙ_ω，ｔ）_Ω×Ｔを出力する。なお、時間周波数成分Ｙ_ω,tは、ウェーブレット変換を用いて計算してもよい。

パラメータ初期値設定部２は、信号分析部３で用いる各パラメータＨ_{ω，φ，ｄ}，Ｕ_ｄ，ｔ，Ｐ_{φ，ｔ，ｄ}の初期値Ｈ_{ω，φ，ｄ} ^（０），Ｕ_ｄ，ｔ ^（０），Ｐ_{φ，ｔ，ｄ} ^（０）を設定する。初期値として、乱数を用いて適当な値を設定すればよい。

出力部５は、信号分析部３で求めた各パラメータＨ_{ω，φ，ｄ}，Ｕ_ｄ，ｔ，Ｐ_{φ，ｔ，ｄ}を出力する。

次に、信号分析部３の具体処理について説明する。本実施の形態では、上記（５）式の目的関数を構成する距離尺度として、ユークリッド距離を利用した場合（β＝２）について説明する。

本実施の形態では上記（５）式の目的関数の値を最小化するにあたり、補助関数法を用いる。目的関数Ｊ（θ）に対して、以下の（７）式が成り立つとき、Ｊ^（＋）（θ，＾θ）を補助関数、＾θを補助変数と呼ぶ。

本実施の形態における補助関数Ｊ^（＋）（θ，＾θ）は、∀ω,t,d,μ_ω,t,d＞０，Σ_ｄμ_{ω，ｔ，ｄ}＝１，Σ_ｄ＾Ｙ_{ω，ｔ，ｄ}＝Ｙ_ω，ｔという条件の下で、＾θ＝（＾Ｙ_{ω，ｔ，ｄ}）_Ω×Ｔ×Ｄとすると、以下の（８）式に示すように設計することができる。

ここで、Ｊ（θ）≦Ｊ^（＋）（θ，＾θ）における等号が成立する条件は、以下の（９）式で表される。

μ_{ω，ｔ，ｄ}の与え方は自由であるが、ここでは、より効率的に計算できるようにするため、以下の（１０）式に示すように、μ_{ω，ｔ，ｄ}を与える。

信号分析部３は、収束判定部３４で収束したと判定されるまで、補助変数更新部３１、スケジューリング係数更新部３２、パラメータ更新部３３の処理を繰り返すことにより、各パラメータＨ_{ω，φ，ｄ}，Ｕ_ｄ，ｔ，Ｐ_{φ，ｔ，ｄ}の値を更新していく。以下では、繰り返し回数が第ｓ回目（ｓ＝１，２，・・・，Ｓ）の実行における各部の処理について説明をする。

補助変数更新部３１は、まず、記憶部４に記憶されている、前回（第ｓ−１回目）の実行で得られたパラメータＨ_{ω，φ，ｄ} ^{（ｓ−１）}，Ｕ_ｄ，ｔ ^{（ｓ−１）}，Ｐ_{φ，ｔ，ｄ} ^{（ｓ−１）}を用いて、以下の（９´）式に従って、（ω，ｔ，ｄ）の組み合わせの各々について補助変数＾Ｙ_{ω，ｔ，ｄ} ^{（ｓ−１）}の値を算出する。算出した値＾Ｙ_{ω，ｔ，ｄ} ^{（ｓ−１）}、μ_{ω，ｔ，ｄ} ^{（ｓ−１）}を記憶部４に記憶しておく。

スケジューリング係数更新部３２は、第ｓ回目（s=1,2,…,S）の実行におけるスケジューリング係数ｋ_φ ^（ｓ）を以下の（１１）式に従って計算する。

上記（１１）式に示すようなスケジューリング係数を与えることにより、各実行ステップにおいて、目的関数の値を単調収束させることが可能となる。

また、パラメータ更新部３３は、計算されたスケジューリング係数ｋ_φ ^（ｓ）を用いて、以下の（１３）式に従って、第s回目の実行におけるパラメータＰ_{φ，ｔ，ｄ}の更新値Ｐ_{φ，ｔ，ｄ} ^（ｓ）を、（φ，ｔ，ｄ）の組み合わせの各々について求める。

ここで、変数＾Ｐ_{φ，ｔ，ｄ}は、記憶部４に記憶されている、パラメータＨ_{ω，φ，ｄ} ^{（ｓ−１）}，Ｕ_ｄ，ｔ ^{（ｓ−１）}，Ｐ_{φ，ｔ，ｄ} ^{（ｓ−１）}，補助変数＾Ｙ_{ω，ｔ，ｄ} ^{（ｓ−１）}、μ_{ω，ｔ，ｄ} ^{（ｓ−１）}を用いて、以下の（１４）式に従って、（φ，ｔ，ｄ）の組み合わせの各々について求められる。

ただし、Ｇは、対象の音響信号に応じて各ｄ＝１，・・・，Ｄに対する隣接行列Ａ^（ｄ）を定めた時にとりうるＰの集合であり、予め定められている。

そして、パラメータ更新部３３は、求めた更新値Ｐ_{φ，ｔ，ｄ} ^（ｓ）、及び記憶部４に記憶されている＾Ｙ_{ω，ｔ，ｄ} ^{（ｓ−１）}，Ｕ_ｄ，ｔ ^{（ｓ−１）}，μ_{ω，ｔ，ｄ} ^{（ｓ−１）}を用いて、以下の（１５）式に従って、第s回目の実行におけるパラメータＨ_{ω，φ，ｄ}の更新値Ｈ_{ω，φ，ｄ} ^（ｓ）を、（ω，φ，ｄ）の組み合わせの各々について求める。

さらに、パラメータ更新部３３は、求めた更新値Ｐ_{φ，ｔ，ｄ} ^（ｓ）と更新値Ｈ_{ω，φ，ｄ} ^（ｓ）、及び記憶部４に記憶されている＾Ｙ_{ω，ｔ，ｄ} ^{（ｓ−１）}，Ｕ_ｄ，ｔ ^{（ｓ−１）}，μ_{ω，ｔ，ｄ} ^{（ｓ−１）}を用いて、以下の（１６）式に従って、第s回目の実行におけるパラメータＵ_ｄ，ｔの更新値Ｕ_ｄ，ｔ ^（ｓ）を、（ｄ，ｔ）の組み合わせの各々について求める。

収束判定部３４は、各パラメータが収束したか否かを判定し、収束していない場合には、繰り返し回数ｓを１つ増加させて、補助変数更新部３１、スケジューリング係数更新部３２、及びパラメータ更新部３３の各処理を繰り返す。収束判定部３４は、収束したと判定した場合には、最後に更新されたパラメータ値Ｐ_{φ，ｔ，ｄ} ^（ｓ）、Ｈ_{ω，φ，ｄ} ^（ｓ）、Ｕ_ｄ，ｔ ^（ｓ）を、目的関数を最小化するパラメータＰ_{φ，ｔ，ｄ}、Ｈ_{ω，φ，ｄ}、Ｕ_ｄ，ｔとして出力部５へ出力する。

収束判定部３４において収束したか否かを判定する方法としては、繰り返し回数ｓが予め定めた回数Ｓに達したときに収束したと判定する方法を用いればよい。なお、s-1回目のパラメータを用いたときの目的関数の値とs回目のパラメータを用いたときの目的関数の値との差が、予め定めた閾値よりも小さくなったときに収束したと判定するようにしてもよい。

次に、本実施の形態に係る音響信号パラメータ分析装置の作用について説明する。まず、分析対象の時系列信号として音響信号が音響信号パラメータ分析装置に入力され、図示しないメモリに格納される。そして、音響信号パラメータ分析装置において、図３に示す音響信号分析処理ルーチンが実行される。

まず、ステップ１００において、メモリから、あるフレーム内の音響信号を読み込み、音響信号に対して、短時間フーリエ変換を用いた時間周波数分析を行った結果から、観測時間周波数成分Ｙ_ω,tを各（ω，ｔ）の要素にもつ二次元配列Ｙを生成して、記憶部４に記憶する。

そして、ステップ１０２において、乱数を用いて、パラメータθ＝｛Ｈ、Ｕ、Ｐ｝の初期値を設定して、記憶部４に記憶する。

次にステップ１０４では、上記ステップ１０２で設定されたパラメータθ、又は後述するステップ１０８で更新されたパラメータθと、上記ステップ１００で生成された二次元行列Ｙとに基づいて、上記（９´）式、（１０´）式に従って、補助係数＾Ｙ_{ω，ｔ，ｄ}を各（ω，ｔ，ｄ）の組み合わせについて算出して、記憶部４に記憶する。また、上記ステップ１０４で算出される変数μ_{ω，ｔ，ｄ}も、記憶部４に記憶する。

そして、ステップ１０６では、予め設定された基底の状態数Φ、予め設定された繰り返し回数Ｓ、及び現在の繰り返し回数ｓに基づいて、上記（１１）式に従って、スケジューリング係数ｋ_φを、全ての基底の状態φについて計算して、記憶部４に記憶する。

ステップ１０８では、上記ステップ１０２で設定されたパラメータθ、又は前回のステップ１０８で更新されたパラメータθと、上記ステップ１０４で計算された補助係数＾Ｙ_{ω，ｔ，ｄ}及び変数μ_{ω，ｔ，ｄ}と、予め設定された行列の集合Ｇと用いて、上記（１４）式に従って、変数＾Ｐ_{φ，ｔ，ｄ}を（φ，ｔ，ｄ）の組み合わせの各々について求める。そして、求められた変数＾Ｐ_{φ，ｔ，ｄ}と上記ステップ１０６で計算されたスケジューリング係数ｋ_φとに基づいて、パラメータＰ_{φ，ｔ，ｄ}の更新値を、上記（１３）式に従って、（φ，ｔ，ｄ）の組み合わせの各々について計算して、記憶部４に記憶する。

そして、上記ステップ１０２で設定されたパラメータＵ，Ｈ、又は前回のステップ１０８で更新されたパラメータＵ，Ｈと、上記ステップ１０４で計算された補助係数＾Ｙ_{ω，ｔ，ｄ}及びμ_{ω，ｔ，ｄ}と、上記で計算されたパラメータＰ_{φ，ｔ，ｄ}の更新値とに基づいて、上記（１５）式に従って、パラメータＨ_{ω，φ，ｄ}の更新値を、（ω，φ，ｄ）の組み合わせの各々について計算して、記憶部４に記憶する。

次に、上記ステップ１０２で設定されたパラメータＵ、又は前回のステップ１０８で更新されたパラメータＵと、上記ステップ１０４で計算された補助係数＾Ｙ_{ω，ｔ，ｄ}及びμ_{ω，ｔ，ｄ}と、上記で計算されたパラメータＰ_{φ，ｔ，ｄ}の更新値及びパラメータＨ_{ω，φ，ｄ}の更新値とに基づいて、上記（１６）式に従って、パラメータＵ_ｄ，ｔの更新値を、（ｄ，ｔ）の組み合わせの各々について計算して、記憶部４に記憶する。

次のステップ１１０では、所定の収束条件として、繰り返し回数ｓが、Ｓに到達したか否かを判定し、繰り返し回数ｓがＳに到達していない場合には、所定の収束条件が成立していないと判断して、上記ステップ１０４へ戻り、上記ステップ１０８で更新したパラメータθを用いて、上記ステップ１０４〜ステップ１０８の処理を繰り返す。一方、繰り返し回数ｓがＳに到達した場合には、所定の収束条件が成立したと判断し、ステップ１１２で、上記ステップ１０８で最終的に更新されたパラメータθを出力して、音響信号分析処理ルーチンを終了する。

次に、第２の実施の形態について説明する。なお、第２の実施の形態に係る音響信号パラメータ分析装置の構成は、第１の実施の形態と同様であるため、同一符号を付して説明を省略する。

第２の実施の形態では、目的関数の距離尺度に、ＫＬ ｄｉｖｅｒｇｅｎｃｅを用いている点が、第１の実施の形態と異なっている。

第２の実施の形態に係る音響信号パラメータ分析装置の信号分析部３の具体処理について説明する。本実施の形態では、上記（５）式の目的関数を構成する距離尺度として、Ｋullback-Leibler (KL) divergenceを利用した場合（β＝１）について説明する。

本実施の形態では、以下の（１７）式で表される補助関数Ｊ^（＋）（θ，＾θ）を用いる。

ここで、Ｊ（θ）≦Ｊ^（＋）（θ，＾θ）における等号が成立する条件は、以下の（１８）式で表される。

信号分析部３は、収束判定部３４で収束したと判定されるまで、補助変数更新部３１、スケジューリング係数更新部３２、パラメータ更新部３３の処理を繰り返すことにより、各パラメータＨ_{ω，φ，ｄ}，Ｕ_ｄ，ｔ，Ｐ_{φ，ｔ，ｄ}の値を更新していく。以下では、繰り返し回数が第ｓ回目（ｓ＝１，２，・・・，Ｓ）の実行における各部の処理について説明をする。なお、スケジューリング係数更新部３２及び収束判定部３４の処理は、第１の実施の形態と同様であるため、説明を省略する。

補助変数更新部３１は、まず、記憶部４に記憶されている、前回（第ｓ−１回目）の実行で得られたパラメータＨ_{ω，φ，ｄ} ^{（ｓ−１）}，Ｕ_ｄ，ｔ ^{（ｓ−１）}，Ｐ_{φ，ｔ，ｄ} ^{（ｓ−１）}を用いて、以下の（１８´）式に従って、（ω，ｔ，ｄ）の組み合わせの各々について補助変数η_{ω，ｔ，ｄ} ^{（ｓ−１）}の値を算出する。算出した値η_{ω，ｔ，ｄ} ^{（ｓ−１）}を記憶部４に記憶しておく。

また、パラメータ更新部３３は、計算されたスケジューリング係数ｋ_φ ^（ｓ）を用いて、以下の（１９）式に従って、第s回目の実行におけるパラメータＰ_{φ，ｔ，ｄ}の更新値Ｐ_{φ，ｔ，ｄ} ^（ｓ）を、（φ，ｔ，ｄ）の組み合わせの各々について求める。

ここで、＾Ｐ_{φ，ｔ，ｄ}は、記憶部４に記憶されている、観測時間周波数成分Ｙ_ω，ｔ、パラメータＨ_{ω，φ，ｄ} ^{（ｓ−１）}，Ｕ_ｄ，ｔ ^{（ｓ−１）}，Ｐ_{φ，ｔ，ｄ} ^{（ｓ−１）}，η_{ω，ｔ，ｄ} ^{（ｓ−１）}と、予め定められた基底の状態数Φ及び基底の数Ｄと、予め定められた行列の集合Ｇとを用いて、以下の式に従って、（φ，ｔ，ｄ）の組み合わせの各々について求められる。

そして、パラメータ更新部３３は、求めた更新値Ｐ_{φ，ｔ，ｄ} ^（ｓ）、及び記憶部４に記憶されている補助変数η_{ω，ｔ，ｄ} ^{（ｓ−１）}，パラメータＵ_ｄ，ｔ ^{（ｓ−１）}、観測時間周波数成分Ｙ_ω，ｔを用いて、以下の（２０）式に従って、第s回目の実行におけるパラメータＨ_{ω，φ，ｄ}の更新値Ｈ_{ω，φ，ｄ} ^（ｓ）を、（ω，φ，ｄ）の組み合わせの各々について求める。

さらに、パラメータ更新部３３は、求めた更新値Ｐ_{φ，ｔ，ｄ} ^（ｓ）と更新値Ｈ_{ω，φ，ｄ} ^（ｓ）、及び記憶部４に記憶されている補助変数η_{ω，ｔ，ｄ} ^{（ｓ−１）}，観測時間周波数成分Ｙ_ω，ｔを用いて、以下の（２１）式に従って、第s回目の実行におけるパラメータＵ_ｄ，ｔの更新値Ｕ_ｄ，ｔ ^（ｓ）を、（ｄ，ｔ）の組み合わせの各々について求める。

なお、第２の実施の形態に係る音響信号パラメータ分析装置の他の構成及び作用については、第１の実施の形態と同様であるため、説明を省略する。

次に、第３の実施の形態について説明する。なお、第３の実施の形態に係る音響信号パラメータ分析装置の構成は、第１の実施の形態と同様であるため、同一符号を付して説明を省略する。

第３の実施の形態では、目的関数の距離尺度に、板倉齊藤距離を用いている点が、第１の実施の形態と異なっている。

第３の実施の形態に係る音響信号パラメータ分析装置の信号分析部３の具体処理について説明する。本実施の形態では、上記（５）式の目的関数を構成する距離尺度として、板倉齊藤距離を利用した場合（β＝０）について説明する。

本実施の形態では、以下の（２２）式で表される補助関数Ｊ^（＋）（θ，＾θ）を用いる。

ここで、Ｊ（θ）≦Ｊ^（＋）（θ，＾θ）における等号が成立する条件は、以下の（２３）式、（２４）式で表される。

信号分析部３は、収束判定部３４で収束したと判定されるまで、補助変数更新部３１、スケジューリング係数更新部３２、パラメータ更新部３３の処理を繰り返すことにより、各パラメータＨ_{ω，φ，ｄ}，Ｕ_ｄ，ｔ，Ｐ_{φ，ｔ，ｄ}の値を更新していく。以下では、繰り返し回数が第ｓ回目（ｓ＝１，２，・・・，Ｓ）の実行における各部の処理について説明をする。なお、スケジューリング係数更新部３２及び収束判定部３４の処理は、第１の実施の形態と同様であるため、説明を省略する。

補助変数更新部３１は、まず、記憶部４に記憶されている、前回（第ｓ−１回目）の実行で得られたパラメータＨ_{ω，φ，ｄ} ^{（ｓ−１）}，Ｕ_ｄ，ｔ ^{（ｓ−１）}，Ｐ_{φ，ｔ，ｄ} ^{（ｓ−１）}を用いて、以下の（２３´）式に従って、（ω，ｔ）の組み合わせの各々について補助変数Ｚ_ω，ｔ ^{（ｓ−１）}の値を算出すると共に、以下の（２４´）式に従って、（ω，ｔ，ｄ）の組み合わせの各々について補助変数η_{ω，ｔ，ｄ} ^{（ｓ−１）}の値を算出する。算出した値Ｚ_ω，ｔ ^{（ｓ−１）}、η_{ω，ｔ，ｄ} ^{（ｓ−１）}を記憶部４に記憶しておく。

また、パラメータ更新部３３は、計算されたスケジューリング係数ｋ_φ ^（ｓ）を用いて、以下の（２５）式に従って、第s回目の実行におけるパラメータＰ_{φ，ｔ，ｄ}の更新値Ｐ_{φ，ｔ，ｄ} ^（ｓ）を、（φ，ｔ，ｄ）の組み合わせの各々について求める。

ここで、＾Ｐ_{φ，ｔ，ｄ}は、記憶部４に記憶されている、観測時間周波数成分Ｙ_ω，ｔ、パラメータＨ_{ω，φ，ｄ} ^{（ｓ−１）}，Ｕ_ｄ，ｔ ^{（ｓ−１）}，Ｐ_{φ，ｔ，ｄ} ^{（ｓ−１）}，Ｚ_ω，ｔ ^{（ｓ−１）}、η_{ω，ｔ，ｄ} ^{（ｓ−１）}と、予め定められた基底の数Ｄと、予め設定された行列の集合Ｇとを用いて、以下の式に従って、（φ，ｔ，ｄ）の組み合わせの各々について求められる。

そして、パラメータ更新部３３は、求めた更新値Ｐ_{φ，ｔ，ｄ} ^（ｓ）、及び記憶部４に記憶されている補助変数Ｚ_ω，ｔ ^{（ｓ−１）}、η_{ω，ｔ，ｄ} ^{（ｓ−１）}，パラメータＵ_ｄ，ｔ ^{（ｓ−１）}、観測時間周波数成分Ｙ_ω，ｔを用いて、以下の（２６）式に従って、第s回目の実行におけるパラメータＨ_{ω，φ，ｄ}の更新値Ｈ_{ω，φ，ｄ} ^（ｓ）を、（ω，φ，ｄ）の組み合わせの各々について求める。

さらに、パラメータ更新部３３は、求めた更新値Ｐ_{φ，ｔ，ｄ} ^（ｓ）と更新値Ｈ_{ω，φ，ｄ} ^（ｓ）、及び記憶部４に記憶されている補助変数Ｚ_ω，ｔ ^{（ｓ−１）}、η_{ω，ｔ，ｄ} ^{（ｓ−１）}，観測時間周波数成分Ｙ_ω，ｔを用いて、以下の（２７）式に従って、第s回目の実行におけるパラメータＵ_ｄ，ｔの更新値Ｕ_ｄ，ｔ ^（ｓ）を、（ｄ，ｔ）の組み合わせの各々について求める。

なお、第３の実施の形態に係る音響信号パラメータ分析装置の他の構成及び作用については、第１の実施の形態と同様であるため、説明を省略する。

次に、楽器単音による音響信号を分析対象として、上述した第２の実施の形態に係る手法を適用し、シミュレーション実験を行った結果について説明する。

ここで、収束判定部３４における収束条件として、「反復回数Ｓ＝２００回に到達したか否か」を用い、提案手法の基礎的な動作確認として、楽器単音からスペクトルの状態変化を分析した。実験に用いた信号は、ピアノＣ３の４つの音長の入った音響信号をＭＩＤＩによって合成したものとした。振幅スペクトログラムＹは、短時間フーリエ変換（サンプリング数１６ｋＨｚ，フレーム長６４ｍｓ、フレームシフト３２ｍｓ、Ｈａｎｎｉｎｇ窓）により計算した（図４（Ｂ）参照）。ピアノにおける基底状態の状態遷移はleft-to-rightに起こることを想定し、隣接行列Ａ^（ｄ）は、以下の（３２）式のように予め定めた。

また、基底数Ｄ＝１、状態数Φ＝４を予め定め、パラメータの初期値を乱数により与えた。信号分析を行った結果、例えば、図４（Ａ）、（Ｃ）に示すような、４つの基底の状態Ｈ_{ω，１，１}と〜Ｈ_{ω，４，１}とアクティベーションＵ_１，ｔが得られた。得られた結果から、既定の状態遷移に着目すると、音の立ち上がりから２状態は音長によらず同程度の長さを持ち、減衰の長さが異なっている様子が確認できた。また基底として得られたスペクトルパターンに着目すると、各倍音ごとに異なる減衰をしていることも分かった。

次に、提案法による音源分離の性能評価を行うために、従来のＮＭＦと提案法との間のＳＮＲ（Ｓｉｇｎａｌ−ｔｏ−Ｎｏｉｓｅ Ｒａｔｉｏ）を比較した結果について説明する。提案法として、上述した第１の実施の形態の手法（ＥＵＣ−ｐｒｏｐｏｓｅｄ）、第２の実施の形態の手法（ＫＬ−ｐｒｏｐｏｓｅｄ）、及び第３の実施の形態の手法（ＩＳ−ｐｒｏｐｏｓｅｄ）をそれぞれ採用した。振幅スペクトログラム上で、各音源の元信号Ｉ_ω，ｔと推定した各音源＾Ｉ_ω，ｔとの間のＳＮＲを、以下の式により計算した。

観測信号として、ヴァイオリンによるＤ♭、Ｆ、Ａ♭の混合音を用いた。ビブラートによる基底の状態遷移がleft-to-rightに起こると想定した場合と、ergodicに起こると想定して隣接行列を以下の（３３）式のように設定した場合とについて、それぞれ実験を行った。

基底数Ｄ＝３、各基底の状態数Φ＝５を予め設定した。比較対象とする従来のＮＭＦでは、基底数Ｄ＝３，Ｄ＝１５として分解を行った。また、分離した各信号は元信号のどの音源に属するか分らないため、元信号の全てとのＳＮＲを求め、最も大きな値を持つ音源への割り当てを行った。もし割り当ての存在しない音源があった場合は、その音源はＳＮＲの算出には用いず、代わりに、割り当ての存在しなかった音源数の割合をエラー率として計算した。

従来法及び提案法の各々において、パラメータの初期値を乱数によって与え、１０回の試行を行ったところ、図５に示すような、ＳＮＲとエラー率が得られた。また、上述した第２の実施の形態の手法を用いて音源分離した結果、図６（Ａ）に示すような、３音源を混合した観測スペクトログラムから、図６（Ｂ）、（Ｃ）に示すように、分離した１音源のスペクトログラムと、その基底の状態遷移とが得られた。上記図６（Ｃ）に示すように、基底の状態変化から、ヴァイオリンのビブラートが５つの状態遷移として学習できていることが確認できた。

また、ＳＮＲに関して、第３の実施の形態の手法を除いて、提案法では、基底数が１５のときの従来法と同程度の性能を示している。従来法では、各音源がいくつかの信号に分解され、それらの統合に元信号を利用しているが、提案法においては、音源ごとに分解された信号が階層的にまとまった形で得られる。このことから、提案法は、元信号の事前知識を必要としないで音源分離を行うことのできる手法としても有効であることが確認できた。

また、上記図５に示す結果から、提案法において、距離尺度として板倉齊藤距離（ＩＳ−ｐｒｏｐｏｓｅｄ）を用いた場合よりも、ユークリッド距離（ＥＵＣ−ｐｒｏｐｏｓｅｄ）とＫＬダイバージェンス（ＫＬ−ｐｒｏｐｏｓｅｄ）を距離尺度にして用いた場合の結果の方が、性能が良いことが分かった。このように、距離尺度を変えることで、より精度良い音源分離を行うことができる。

次に、上述した第１の実施の形態の手法を用いた場合における目的関数の値の変化を、スケジューリング係数を用いない手法と比較した結果について説明する。状態遷移を表すＰの更新に、スケジューリング係数を導入することによる効果を検証するための実験を行った。観測信号は、上述した音源分離の性能評価のための実験と同じものを用いた。第１実施の形態の方法を用い、収束判定部３４の収束条件として「反復回数ｓ＝１００回に到達したか否か」を用いた。スケジューリング係数を導入しない場合には、ｋ_１ ^（ｓ）＝１とし、残りのｋ_ｉ ^（ｓ）（ｉ＝２，・・・，Φ）を０にして計算した。つまり、上記（１３）式の代わりに、以下の式に従って、パラメータＰ_{φ，ｔ，ｄ}の更新値Ｐ_{φ，ｔ，ｄ} ^（ｓ）を計算した。

図７に示すように、スケジューリング係数を用いた場合（図７の実線を参照）と用いなかった場合（図７の破線を参照）をそれぞれ３回ずつ試行したときの目的関数の変化の様子が得られた。スケジューリング係数を導入することにより、局所解が回避できていることが確認できた。

以上説明したように、本発明の実施の形態に係る音響信号パラメータ分析装置によれば、時系列信号の観測時間周波数成分Ｙ_ω,tを要素とする二次元配列Ｙとの距離をβダイバージェンスで表わした目的関数の値が小さくなるように、各基底ｄにおける各基底の状態φを表わすパラメータＨ_ω,φ,d、状態遷移を表わすパラメータＰ_φ,t,d、及びアクティベーションＵ_d,tを更新することにより、時系列信号から、頻出のパターンとパターンの出現の時間的な関連とを同時に精度良く抽出することができる。

また、楽器単音のスペクトルが、限られた数のスペクトルパターンの状態遷移によって構成されるという階層構造を利用し、階層的にスパース表現をする枠組みを提案することにより、信号から頻出のパターンとそれらのパターンの出現の時間的な関連とを同時に抽出可能とした。

また、任意の実数βについて、β-divergenceを距離尺度として用いてスペクトルの時間的な変化を分析できるようにＮＭＦを拡張することができる。これにより、入力される信号に適した距離尺度を用いることができるので、様々な信号に対して、音の特性に関する事前知識や特殊な前処理を必要とせずに、スペクトルの時間的な変化を精度良く分析することができる。

また、上述した背景技術の説明において、上記の非特許文献４の技術は、距離尺度を板倉齊藤距離とした場合のＮＭＦモデルを拡張したものに相当すると述べたが、本発明の実施の形態と比較すると、拡張したモデルのパラメータを求める具体的な方法が全く異なる。また、本発明によれば、距離尺度として板倉齊藤距離以外の距離尺度を用いることも可能であると共に、非特許文献４の手法よりも早く収束させることができる、という利点がある。また、非特許文献４に記載の技術でも、ヒューリスティックな方法（非特許文献４におけるＥＭ−ＭＵ）を導入することで収束を早くしようとする試みはされているものの、そのアルゴリズムが必ず収束することが証明されていない。一方、本発明の手法によれば、必ず収束させることができ、実装も容易である。

なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

例えば、上述の音響信号パラメータ分析装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。

また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。

本発明による信号の分析は、多重音解析のための分析手法として、潜在的に、音楽音響信号からの自動採譜や、音源分離、楽器音認識、音楽加工などに応用されることが考えられる。例えば、ピアノ曲に本発明の分析を適用した場合には、各音高のアタックやサステインやリリースといった状態とそれらの出現とを同時に推定することが期待される。これは各音のオンセットについて、より正確に手掛かりを与え、自動採譜の性能を高めると考えられる。また、弦楽器や歌声によるビブラートを、ピアノ曲と同じようにいくつかの区分定常なスペクトルの状態遷移として表現出来ることが期待でき、多重音からの音源分離の性能を向上させることが可能であり、より精密な音楽加工に応用することが可能であると考えられる。

１ 時間周波数解析部
２ パラメータ初期値設定部
３ 信号分析部
４ 記憶部
５ 出力部
３１ 補助変数更新部
３２ スケジューリング係数更新部
３３ パラメータ更新部
３４ 収束判定部

Claims (9)

1. 時系列信号を入力として、観測時間周波数成分Ｙ_ω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙを出力する時間周波数分解手段と、
   複数の基底ｄの各々における全ての基底の状態φを表わすパラメータＨ_ω,φ,dを要素にもつ三次元配列Ｈ、前記複数の基底ｄの各々に対して各時刻ｔに前記基底の状態が前記三次元配列Ｈの要素のうちの何れであるかを表わすパラメータＰ_φ,t,dを要素にもつ三次元配列Ｐ、及び前記複数の基底ｄの各々に対する各時刻ｔのゲインを表すパラメータＵ_d,tを要素にもつ二次元配列Ｕの各々の初期値を設定するパラメータ初期値設定手段と、
   全てのパラメータＨ_ω,φ,dが非負値であり、全てのパラメータＵ_d,tが非負値であり、かつ、前記三次元配列Ｐが、前記複数の基底ｄの各々について予め許可された前記基底の状態の状態遷移に基づくものである、という条件の下で、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を、（ｄ，φ）の全ての組み合わせについて足し合わせたモデルについて、前記時間周波数分解手段によって出力された二次元配列Ｙと前記モデルとの距離をβダイバージェンスで表わした目的関数の値が小さくなるように、前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの各々の値を更新する信号分析手段と、
   を含む信号分析装置。
2. 前記信号分析手段は、
   前記時間周波数分解手段によって出力された二次元配列Ｙの観測時間周波数成分Ｙ_ω,tと、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を（ｄ，φ）の全ての組み合わせについて足し合わせたものとの差分に、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を全ての基底の状態φについて合計したものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算した商μ_{ω，ｔ，ｄ}を、乗算したものに、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて合計したものを加算して、補助変数＾Ｙ_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算する補助変数計算手段と、
   全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、
   全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算するスケジューリング係数計算手段と、
   前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記パラメータＰ_φ,t,dと前記補助変数＾Ｙ_{ω，ｔ，ｄ}との積から、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を減算したものを２乗して、前記商μ_{ω，ｔ，ｄ}で除算した商を、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、
   前記補助変数＾Ｙ_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値の２乗と前記パラメータＵ_d,tとを乗算した積を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積の２乗を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたもので除算して、前記パラメータＨ_ω,φ,dの更新値を（ω，φ，ｄ）の全ての組み合わせについて算出し、
   前記パラメータＨ_ω,φ,dの更新値と前記補助変数＾Ｙ_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値の２乗とを乗算した積を、前記商μ_{ω，ｔ，ｄ}で除算した商を、（ω，ｔ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積の２乗を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたもので除算して、前記パラメータＵ_d,tの更新値を（ｄ，ｔ）の全ての組み合わせについて算出するパラメータ更新手段と、
   前記繰り返し回数Ｓに到達するまで、前記補助変数計算手段による計算、前記スケジューリング係数計算手段による計算、及び前記パラメータ更新手段による更新を繰り返し行う収束判定手段と、を備えた請求項１記載の信号分析装置。
3. 前記信号分析手段は、
   前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて足し合わせたものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算して、補助変数η_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算する補助変数計算手段と、
   全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、
   全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算するスケジューリング係数計算手段と、
   前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記補助変数η_{ω，ｔ，ｄ}と、前記パラメータＰ_φ,t,dと、前記観測時間周波数成分Ｙ_ω,tと、前記補助変数η_{ω，ｔ，ｄ}及び前記観測時間周波数成分Ｙ_ω,tの積を前記パラメータＨ_ω,φ,d及び前記パラメータＵ_d,tの積で除算した商について対数をとったものと、を乗算した積を、前記観測時間周波数成分Ｙ_ω,tを前記基底の数Ｄと前記状態数Φとの積で除算した商で減算し、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を加算したものを、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、
   前記補助変数η_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積を各時刻ｔについて足し合わせたもので除算して、前記パラメータＨ_ω,φ,dの更新値を（ω，φ，ｄ）の全ての組み合わせについて算出し、
   前記補助変数η_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を（ω，φ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積を（ω，φ）の全ての組み合わせについて足し合わせたもので除算して、前記パラメータＵ_d,tの更新値を（ｄ，ｔ）の全ての組み合わせについて算出するパラメータ更新手段と、
   前記繰り返し回数Ｓに到達するまで、前記補助変数計算手段による計算、前記スケジューリング係数計算手段による計算、及び前記パラメータ更新手段による更新を繰り返し行う収束判定手段と、を備えた請求項１記載の信号分析装置。
4. 前記信号分析手段は、
   前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を（ｄ，φ）の全ての組み合わせについて足し合わせて、補助変数Ｚ_ω，ｔを（ω，ｔ）の全ての組み合わせについて計算すると共に、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて足し合わせたものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算して、補助変数η_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算する補助変数計算手段と、
   全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、
   全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算するスケジューリング係数計算手段と、
   前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記パラメータＰ_φ,t,d、前記補助変数η_{ω，ｔ，ｄ}の２乗、及び前記観測時間周波数成分Ｙ_ω,tの積を、前記パラメータＨ_ω,φ,d及び前記パラメータＵ_d,tの積で除算した商と、前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの積を前記補助変数Ｚ_ω，ｔで除算した商と、前記補助変数Ｚ_ω，ｔを前記観測時間周波数成分Ｙ_ω,tで除算した商について対数をとったものから２を減算して前記基底の数Ｄで除算した商と、を加算したものを、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、
   前記パラメータＰ_φ,t,dの更新値と前記補助変数η_{ω，ｔ，ｄ}の２乗と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を、前記パラメータＵ_d,tで除算した商を、各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積を、前記補助変数Ｚ_ω，ｔで除算した商を、各時刻ｔについて足し合わせたもので除算した商の平方根を、前記パラメータＨ_ω,φ,dの更新値として（ω，φ，ｄ）の全ての組み合わせについて算出し、
   前記パラメータＰ_φ,t,dの更新値と前記補助変数η_{ω，ｔ，ｄ}の２乗と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を、前記パラメータＨ_ω,φ,dの更新値で除算した商を、（ω，φ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積を前記補助変数Ｚ_ω，ｔで除算した商を、（ω，φ）の全ての組み合わせについて足し合わせたもので除算した商の平方根を、前記パラメータＵ_d,tの更新値として（ｄ，ｔ）の全ての組み合わせについて算出するパラメータ更新手段と、
   前記繰り返し回数Ｓに到達するまで、前記補助変数計算手段による計算、前記スケジューリング係数計算手段による計算、及び前記パラメータ更新手段による更新を繰り返し行う収束判定手段と、を備えた請求項１記載の信号分析装置。
5. 時系列信号を入力として、観測時間周波数成分Ｙ_ω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙを出力し、
   複数の基底ｄの各々における全ての基底の状態φを表わすパラメータＨ_ω,φ,dを要素にもつ三次元配列Ｈ、前記複数の基底ｄの各々に対して各時刻ｔに前記基底の状態が前記三次元配列Ｈの要素のうちの何れであるかを表わすパラメータＰ_φ,t,dを要素にもつ三次元配列Ｐ、及び前記複数の基底ｄの各々に対する各時刻ｔのゲインを表すパラメータＵ_d,tを要素にもつ二次元配列Ｕの各々の初期値を設定し、
   全てのパラメータＨ_ω,φ,dが非負値であり、全てのパラメータＵ_d,tが非負値であり、かつ、前記三次元配列Ｐが、前記複数の基底ｄの各々について予め許可された前記基底の状態の状態遷移に基づくものである、という条件の下で、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を、（ｄ，φ）の全ての組み合わせについて足し合わせたモデルについて、前記時間周波数分解手段によって出力された二次元配列Ｙと前記モデルとの距離をβダイバージェンスで表わした目的関数の値が小さくなるように、前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの各々の値を更新する
   ことを特徴とする信号分析方法。
6. 前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの各々の値を更新することは、
   前記出力された二次元配列Ｙの観測時間周波数成分Ｙ_ω,tと、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を（ｄ，φ）の全ての組み合わせについて足し合わせたものとの差分に、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を全ての基底の状態φについて合計したものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算した商μ_{ω，ｔ，ｄ}を、乗算したものに、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて合計したものを加算して、補助変数＾Ｙ_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算し、
   全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、
   全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、
   前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記パラメータＰ_φ,t,dと前記補助変数＾Ｙ_{ω，ｔ，ｄ}との積から、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を減算したものを２乗して、前記商μ_{ω，ｔ，ｄ}で除算した商を、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、
   前記補助変数＾Ｙ_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値の２乗と前記パラメータＵ_d,tとを乗算した積を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積の２乗を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたもので除算して、前記パラメータＨ_ω,φ,dの更新値を（ω，φ，ｄ）の全ての組み合わせについて算出し、
   前記パラメータＨ_ω,φ,dの更新値と前記補助変数＾Ｙ_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値の２乗とを乗算した積を、前記商μ_{ω，ｔ，ｄ}で除算した商を、（ω，ｔ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積の２乗を、前記商μ_{ω，ｔ，ｄ}で除算した商を、各時刻ｔについて足し合わせたもので除算して、前記パラメータＵ_d,tの更新値を（ｄ，ｔ）の全ての組み合わせについて算出し、
   前記繰り返し回数Ｓに到達するまで、前記補助変数＾Ｙ_{ω，ｔ，ｄ}の計算、前記スケジューリング係数ｋ_φ ^（ｓ）の計算、前記パラメータＨ_ω,φ,dの更新値の算出、前記パラメータＰ_φ,t,dの更新値の算出、及び前記パラメータＵ_d,tの更新値の算出を繰り返し行う
   ことを特徴とする請求項５記載の信号分析方法。
7. 前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの各々の値を更新することは、
   前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて足し合わせたものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算して、補助変数η_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算し、
   全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、
   全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、
   前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記補助変数η_{ω，ｔ，ｄ}と、前記パラメータＰ_φ,t,dと、前記観測時間周波数成分Ｙ_ω,tと、前記補助変数η_{ω，ｔ，ｄ}及び前記観測時間周波数成分Ｙ_ω,tの積を前記パラメータＨ_ω,φ,d及び前記パラメータＵ_d,tの積で除算した商について対数をとったものと、を乗算した積を、前記観測時間周波数成分Ｙ_ω,tを前記基底の数Ｄと前記状態数Φとの積で除算した商で減算し、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を加算したものを、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、
   前記補助変数η_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積を各時刻ｔについて足し合わせたもので除算して、前記パラメータＨ_ω,φ,dの更新値を（ω，φ，ｄ）の全ての組み合わせについて算出し、
   前記補助変数η_{ω，ｔ，ｄ}と前記パラメータＰ_φ,t,dの更新値と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を（ω，φ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積を（ω，φ）の全ての組み合わせについて足し合わせたもので除算して、前記パラメータＵ_d,tの更新値を（ｄ，ｔ）の全ての組み合わせについて算出し、
   前記繰り返し回数Ｓに到達するまで、前記補助変数η_{ω，ｔ，ｄ}の計算、前記スケジューリング係数ｋ_φ ^（ｓ）の計算、前記パラメータＨ_ω,φ,dの更新値の算出、前記パラメータＰ_φ,t,dの更新値の算出、及び前記パラメータＵ_d,tの更新値の算出を繰り返し行う
   ことを特徴とする請求項５記載の信号分析方法。
8. 前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの各々の値を更新することは、
   前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を（ｄ，φ）の全ての組み合わせについて足し合わせて、補助変数Ｚ_ω，ｔを（ω，ｔ）の全ての組み合わせについて計算すると共に、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとを乗算した積を全ての基底の状態φについて足し合わせたものを、前記パラメータＨ_ω,φ,dと前記パラメータＰ_φ,t,dと前記パラメータＵ_d,tとの積を（ｄ，φ）の全ての組み合わせについて足し合わせたもので除算して、補助変数η_{ω，ｔ，ｄ}を（ω，ｔ，ｄ）の全ての組み合わせについて計算し、
   全ての基底の状態φのうちの１つ目の状態φについて、予め定められた基底の状態数Φの２倍から２を減算した値を、予め定められた繰り返し回数Ｓと前記状態数Φとの積で除算した商に、現時点での繰り返し回数ｓを乗算し、２を前記状態数Φで除算した商を加算し１を減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、
   全ての基底の状態φのうちの残りの状態φの各々について、２を前記状態数Φで除算した商を、２を前記繰り返し回数Ｓと前記状態数Φとの積で除算した商に前記繰り返し回数ｓを乗算したもので減算したもの、及び前記基底の状態数Φの逆数の何れか大きい方を、スケジューリング係数ｋ_φ ^（ｓ）として計算し、
   前記予め許可された前記基底の状態の状態遷移に基づく前記パラメータＰ_φ,t,dのうち、前記パラメータＰ_φ,t,d、前記補助変数η_{ω，ｔ，ｄ}の２乗、及び前記観測時間周波数成分Ｙ_ω,tの積を、前記パラメータＨ_ω,φ,d及び前記パラメータＵ_d,tの積で除算した商と、前記パラメータＨ_ω,φ,d、前記パラメータＰ_φ,t,d、及び前記パラメータＵ_d,tの積を前記補助変数Ｚ_ω，ｔで除算した商と、前記補助変数Ｚ_ω，ｔを前記観測時間周波数成分Ｙ_ω,tで除算した商について対数をとったものから２を減算して前記基底の数Ｄで除算した商と、を加算したものを、（ω，ｔ，ｄ，φ）の全ての組み合わせについて足し合わせたものが最小となる、前記パラメータＰ_φ,t,dを（φ，ｔ，ｄ）の全ての組み合わせについて求め、変数＾Ｐ_φ,t,dとし、（φ，ｔ，ｄ）の全ての組み合わせについて、前記スケジューリング係数ｋ_ｎ ^（ｓ）と前記変数＾Ｐ_φ-(n-1),t,dとの積を、ｎ＝１からｎ＝Φまで足し合わせたものを、前記パラメータＰ_φ,t,dの更新値として算出し、
   前記パラメータＰ_φ,t,dの更新値と前記補助変数η_{ω，ｔ，ｄ}の２乗と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を、前記パラメータＵ_d,tで除算した商を、各時刻ｔについて足し合わせたものを、前記パラメータＰ_φ,t,dの更新値と前記パラメータＵ_d,tと乗算した積を、前記補助変数Ｚ_ω，ｔで除算した商を、各時刻ｔについて足し合わせたもので除算した商の平方根を、前記パラメータＨ_ω,φ,dの更新値として（ω，φ，ｄ）の全ての組み合わせについて算出し、
   前記パラメータＰ_φ,t,dの更新値と前記補助変数η_{ω，ｔ，ｄ}の２乗と前記観測時間周波数成分Ｙ_ω,tとを乗算した積を、前記パラメータＨ_ω,φ,dの更新値で除算した商を、（ω，φ）の全ての組み合わせについて足し合わせたものを、前記パラメータＨ_ω,φ,dの更新値と前記パラメータＰ_φ,t,dの更新値とを乗算した積を前記補助変数Ｚ_ω，ｔで除算した商を、（ω，φ）の全ての組み合わせについて足し合わせたもので除算した商の平方根を、前記パラメータＵ_d,tの更新値として（ｄ，ｔ）の全ての組み合わせについて算出し、
   前記繰り返し回数Ｓに到達するまで、前記補助変数Ｚ_ω，ｔ及び前記補助変数η_{ω，ｔ，ｄ}の計算、前記スケジューリング係数ｋ_φ ^（ｓ）の計算、前記パラメータＨ_ω,φ,dの更新値の算出、前記パラメータＰ_φ,t,dの更新値の算出、及び前記パラメータＵ_d,tの更新値の算出を繰り返し行う
   ことを特徴とする請求項５記載の信号分析方法。
9. 請求項１〜請求項４の何れか１項に記載の信号分析装置の各手段としてコンピュータを機能させるためのプログラム。