(57)【要約】【課題】時系列信号から、頻出のパターンとパターンの出現の時間的な関連とを精度良く抽出することができるようにする。【解決手段】信号分析部３によって、パラメータＨω,φ,dとパラメータＰφ,t,dとパラメータＵd,tとの積を、（ｄ，φ）の全ての組み合わせについて足し合わせたモデルについて、時間周波数解析部１によって出力された二次元配列Ｙとモデルとの距離をβダイバージェンスで表わした目的関数の値が小さくなるように、パラメータＨω,φ,d、パラメータＰφ,t,d、及びパラメータＵd,tの各々の値を更新する。収束判定部３４は、繰り返し回数Ｓに到達するまで、補助変数更新部３１による計算、スケジューリング係数更新部３２による計算、及びパラメータ更新部３３による更新を繰り返し行う。【選択図】図１
2012027196.tif 000002
【特許請求の範囲】
【請求項１】
時系列信号を入力として、観測時間周波数成分Ｙω,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の更新値の算出を繰り返し行う
ことを特徴とする請求項５記載の信号分析方法。
【請求項９】
請求項１〜請求項４の何れか１項に記載の信号分析装置の各手段としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号分析装置、方法、及びプログラムに係り、特に、時系列信号から、信号パラメータを分析する信号分析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
従来、信号の分析に関しては、信号の時間周波数表現に基づく各種の分析方法が広く知られている。信号内に繰り返して現れるパターン等のような高度な構造を抽出する分析手法として非負値行列分解（ＮＭＦ：Ｎｏｎ‐ｎｅｇａｔｉｖｅ  Ｍａｔｒｉｘ  Ｆａｃｔｏｒｉｚａｔｉｏｎ）が知られている（例えば、非特許文献１）。
【０００３】
ＮＭＦは、自動採譜や、モノラル混合信号からの音源の分離に適用されている（例えば、非特許文献２、非特許文献３）。
【０００４】
また、ポリフォニックオーディオデータを確率モデルで表現することにより、信号分離を行う技術も知られている（非特許文献４）。これは、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を目的関数とし、モデルパラメータの推定を、目的関数の値を最小化する問題として定式化し、反復的に一つの解を探索するアルゴリズムを導出する。
【００２２】
まず、本発明の原理について説明する。まず、信号のモデルについて説明する。
【００２３】
本発明で提案するモデルは、ＮＭＦによるスペクトログラムの分解表現を拡張したモデルである。ＮＭＦによるスペクトログラムの分解表現は、観測スペクトログラムを非負値行列
【００２４】
【数１】
2012027196.tif 000003
【００２５】
と見立て、
【００２６】
【数２】
2012027196.tif 000004
【００２７】
となるような基底
【００２８】
【数３】
2012027196.tif 000005
【００２９】
とアクティベーション
【００３０】
【数４】
2012027196.tif 000006
【００３１】
の２つの非負値行列を決めることによって得られる。ここで、ω（ω＝１，２，…，Ω），ｔ（ｔ＝１，２，・・・，Ｔ）はそれぞれ周波数と時刻に対応するインデックスである。また、Ｙω，ｔは、行列Ｙの（ω，ｔ）成分であり、Ｒ≧０，Ω×Ｔは、各成分が０以上の実数からなるΩ×Ｔ行列の集合を表す。
【００３２】
このモデルは、観測スペクトログラムがＤ個のスペクトル（基底）とアクティベーションとの積によって表現される、というモデルになっている。スペクトルとアクティベーションのペアが、楽器音１音高となることを狙ったモデルである。
【００３３】
実際の楽器音（信号）は、時間に伴ってスペクトルが変化することを考慮すると、スペクトルが各時刻ｔにおいて一つの状態であるとみなし、時刻ｔにおけるｄ番目の基底の状態をＨｄ，ｔ（φｄ，ｔ）で表すことによって、上記（１）式を、以下の（２）式に拡張する。
【００３４】
【数５】
2012027196.tif 000007
【００３５】
ここで、基底の状態も限られたΦ個しか存在せず、時間に伴ってそれらの基底の状態を遷移するものと仮定し、Ｈ＝（Ｈω，φ，ｄ）Ω×Φ×Ｄ∈Ｒ≧０，Ω×Φ×Ｄを用いてモデルを表現し直す。各ｄ＝１，・・・，Ｄに対するΦ個の基底の状態ｈφ（ｄ）＝（Ｈ１，φ，ｄ，・・・，ＨΩ，φ，ｄ）Ｔ（φ＝１，・・・，Φ）に関して、許される状態遷移を隣接行列Ａ（ｄ）＝（Ａφ，φ（ｄ））Φ×Φを用いて表す。状態ｉから状態ｊへの遷移が許される場合、Ａｉ，ｊ（ｄ）＝１とし、遷移が許されない場合、Ａｉ，ｊ（ｄ）＝０とする。例えば，以下の（３）式で表される隣接行列Ａは、状態のインデックス１→２、２→３、３→１の遷移のみが許されることを表している。
【００３６】
【数６】
2012027196.tif 000008
【００３７】
ここで、Ｐ＝（Ｐφ，ｔ，ｄ）Φ×Ｔ×Ｄ∈Ｒ≧０，Φ×Ｔ×Ｄを考え、ｄ＝１，・・・，Ｄに対して各時刻ｔに基底の状態がｈφ（ｄ）のときにＰφ，ｔ，ｄ＝１とし、そうでないときにＰφ，ｔ，ｄ＝０とすることで、提案モデルは、図２に示すように、（Ｈω，φ，ｄ）Ω×Φ×Ｄと、（Ｐφ，ｔ，ｄ）Φ×Ｔ×Ｄと、（Ｕt，ｄ）T×Dと用いて表わされ、以下の（４）式で表現することができる。
【００３８】
【数７】
2012027196.tif 000009
【００３９】
パラメータＨω,φ,dは、各基底ｄにおける全ての基底の状態φを表わす。パラメータＰφ,t,dは、各基底ｄに対して各時刻ｔに基底の状態が、パラメータＨω,φ,dを要素にもつ三次元配列Ｈの要素のうちの何れであるかを表わす。パラメータＵd,tは、各基底ｄに対する各時刻ｔのゲインを表わす。
【００４０】
今、各ｄ＝１，・・・，Ｄに対するＡ（ｄ）を固定した時にとりうるＰの集合をＧとする。すると、本発明により解くべき問題は、観測されたＹから、三次元行列Ｈ，二次元行列Ｕ，及び三次元行列Ｐ∈Ｇを推定することに帰着される。
【００４１】
つまり、行列Ｙに対して、上記（４）式を満たす各Ｈω，φ，ｄ，Ｕｄ，ｔ，Ｐφ，ｔ，ｄを求める。具体的には、観測結果Ｙとモデル（上記（４）式の右辺）との距離（近さ）を目的関数とし、所定の条件の下で目的関数の値を最小化することにより、各Ｈω，φ，ｄ，Ｕｄ，ｔ，Ｐφ，ｔ，ｄを求める。観測結果Ｙとの距離尺度としては、β-divergenceを用いる。具体的には、観測されたＹから、
【００４２】
【数８】
2012027196.tif 000010
【００４３】
という条件のもとで、以下の（５）式で表される目的関数Ｊ（θ）の値を最小化するθ＝｛Ｈ，Ｕ，Ｐ｝を求める。
【００４４】
【数９】
2012027196.tif 000011
【００４５】
ここで、関数Ｄβは以下の（６）式で定義される。
【００４６】
【数１０】
2012027196.tif 000012
【００４７】
次に、音響信号の信号パラメータを分析して出力する音響信号パラメータ分析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【００４８】
図１に示すように、第１の実施の形態に係る音響信号パラメータ分析装置は、ＣＰＵと、ＲＡＭと、後述する音響信号分析処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【００４９】
音響信号パラメータ分析装置は、時間周波数解析部１と、パラメータ初期値設定部２と、信号分析部３と、記憶部４と、出力部５とを備えている。また、信号分析部３は、補助変数更新部３１と、スケジューリング係数更新部３２と、パラメータ更新部３３と、収束判定部３４とを備えている。
【００５０】
時間周波数解析部１は、時系列信号としての観測された音響信号ｙを入力として、時間周波数成分（観測時間周波数成分）Ｙω,t（ω＝１,・・・,Ω,ｔ＝１,・・・,Ｔは、それぞれ周波数、時刻に対応するインデックスを示す。）を各（ω，ｔ）の要素にもつ二次元配列Ｙを計算し、信号分析部３に出力する。また、計算した時間周波数成分Ｙω,tを、記憶部４に記憶しておく。より詳細には、時間周波数解析部１は、時系列信号ｙを入力として、短時間フーリエ変換（Short-Time Fourier Transform；ＳＴＦＴ）を用いて時間周波数解析を行うことにより時間周波数成分Ｙω,tを計算し、時間周波数成分Ｙω，ｔを格納した行列（振幅スペクトログラム）Ｙ＝（Ｙω，ｔ）Ω×Ｔを出力する。なお、時間周波数成分Ｙω,tは、ウェーブレット変換を用いて計算してもよい。
【００５１】
パラメータ初期値設定部２は、信号分析部３で用いる各パラメータＨω，φ，ｄ，Ｕｄ，ｔ，Ｐφ，ｔ，ｄの初期値Ｈω，φ，ｄ（０），Ｕｄ，ｔ（０），Ｐφ，ｔ，ｄ（０）を設定する。初期値として、乱数を用いて適当な値を設定すればよい。
【００５２】
出力部５は、信号分析部３で求めた各パラメータＨω，φ，ｄ，Ｕｄ，ｔ，Ｐφ，ｔ，ｄを出力する。
【００５３】
次に、信号分析部３の具体処理について説明する。本実施の形態では、上記（５）式の目的関数を構成する距離尺度として、ユークリッド距離を利用した場合（β＝２）について説明する。
【００５４】
本実施の形態では上記（５）式の目的関数の値を最小化するにあたり、補助関数法を用いる。目的関数Ｊ（θ）に対して、以下の（７）式が成り立つとき、Ｊ（＋）（θ，＾θ）を補助関数、＾θを補助変数と呼ぶ。
【００５５】
【数１１】
2012027196.tif 000013
【００５６】
本実施の形態における補助関数Ｊ（＋）（θ，＾θ）は、∀ω,t,d,μω,t,d＞０，Σｄμω，ｔ，ｄ＝１，Σｄ＾Ｙω，ｔ，ｄ＝Ｙω，ｔという条件の下で、＾θ＝（＾Ｙω，ｔ，ｄ）Ω×Ｔ×Ｄとすると、以下の（８）式に示すように設計することができる。
【００５７】
【数１２】
2012027196.tif 000014
【００５８】
ここで、Ｊ（θ）≦Ｊ（＋）（θ，＾θ）における等号が成立する条件は、以下の（９）式で表される。
【００５９】
【数１３】
2012027196.tif 000015
【００６０】
μω，ｔ，ｄの与え方は自由であるが、ここでは、より効率的に計算できるようにするため、以下の（１０）式に示すように、μω，ｔ，ｄを与える。
【００６１】
【数１４】
2012027196.tif 000016
【００６２】
信号分析部３は、収束判定部３４で収束したと判定されるまで、補助変数更新部３１、スケジューリング係数更新部３２、パラメータ更新部３３の処理を繰り返すことにより、各パラメータＨω，φ，ｄ，Ｕｄ，ｔ，Ｐφ，ｔ，ｄの値を更新していく。以下では、繰り返し回数が第ｓ回目（ｓ＝１，２，・・・，Ｓ）の実行における各部の処理について説明をする。
【００６３】
補助変数更新部３１は、まず、記憶部４に記憶されている、前回（第ｓ−１回目）の実行で得られたパラメータＨω，φ，ｄ（ｓ−１），Ｕｄ，ｔ（ｓ−１），Ｐφ，ｔ，ｄ（ｓ−１）を用いて、以下の（９´）式に従って、（ω，ｔ，ｄ）の組み合わせの各々について補助変数＾Ｙω，ｔ，ｄ（ｓ−１）の値を算出する。算出した値＾Ｙω，ｔ，ｄ（ｓ−１）、μω，ｔ，ｄ（ｓ−１）を記憶部４に記憶しておく。
【００６４】
【数１５】
2012027196.tif 000017
【００６５】
スケジューリング係数更新部３２は、第ｓ回目（s=1,2,…,S）の実行におけるスケジューリング係数ｋφ（ｓ）を以下の（１１）式に従って計算する。
【００６６】
【数１６】
2012027196.tif 000018
【００６７】
上記（１１）式に示すようなスケジューリング係数を与えることにより、各実行ステップにおいて、目的関数の値を単調収束させることが可能となる。
【００６８】
また、パラメータ更新部３３は、計算されたスケジューリング係数ｋφ（ｓ）を用いて、以下の（１３）式に従って、第s回目の実行におけるパラメータＰφ，ｔ，ｄの更新値Ｐφ，ｔ，ｄ（ｓ）を、（φ，ｔ，ｄ）の組み合わせの各々について求める。
【００６９】
【数１７】
2012027196.tif 000019
【００７０】
ここで、変数＾Ｐφ，ｔ，ｄは、記憶部４に記憶されている、パラメータＨω，φ，ｄ（ｓ−１），Ｕｄ，ｔ（ｓ−１），Ｐφ，ｔ，ｄ（ｓ−１），補助変数＾Ｙω，ｔ，ｄ（ｓ−１）、μω，ｔ，ｄ（ｓ−１）を用いて、以下の（１４）式に従って、（φ，ｔ，ｄ）の組み合わせの各々について求められる。
【００７１】
【数１８】
2012027196.tif 000020
【００７２】
ただし、Ｇは、対象の音響信号に応じて各ｄ＝１，・・・，Ｄに対する隣接行列Ａ（ｄ）を定めた時にとりうるＰの集合であり、予め定められている。
【００７３】
そして、パラメータ更新部３３は、求めた更新値Ｐφ，ｔ，ｄ（ｓ）、及び記憶部４に記憶されている＾Ｙω，ｔ，ｄ（ｓ−１），Ｕｄ，ｔ（ｓ−１），μω，ｔ，ｄ（ｓ−１）を用いて、以下の（１５）式に従って、第s回目の実行におけるパラメータＨω，φ，ｄの更新値Ｈω，φ，ｄ（ｓ）を、（ω，φ，ｄ）の組み合わせの各々について求める。
【００７４】
【数１９】
2012027196.tif 000021
【００７５】
さらに、パラメータ更新部３３は、求めた更新値Ｐφ，ｔ，ｄ（ｓ）と更新値Ｈω，φ，ｄ（ｓ）、及び記憶部４に記憶されている＾Ｙω，ｔ，ｄ（ｓ−１），Ｕｄ，ｔ（ｓ−１），μω，ｔ，ｄ（ｓ−１）を用いて、以下の（１６）式に従って、第s回目の実行におけるパラメータＵｄ，ｔの更新値Ｕｄ，ｔ（ｓ）を、（ｄ，ｔ）の組み合わせの各々について求める。
【００７６】
【数２０】
2012027196.tif 000022
【００７７】
収束判定部３４は、各パラメータが収束したか否かを判定し、収束していない場合には、繰り返し回数ｓを１つ増加させて、補助変数更新部３１、スケジューリング係数更新部３２、及びパラメータ更新部３３の各処理を繰り返す。収束判定部３４は、収束したと判定した場合には、最後に更新されたパラメータ値Ｐφ，ｔ，ｄ（ｓ）、Ｈω，φ，ｄ（ｓ）、Ｕｄ，ｔ（ｓ）を、目的関数を最小化するパラメータＰφ，ｔ，ｄ、Ｈω，φ，ｄ、Ｕｄ，ｔとして出力部５へ出力する。
【００７８】
収束判定部３４において収束したか否かを判定する方法としては、繰り返し回数ｓが予め定めた回数Ｓに達したときに収束したと判定する方法を用いればよい。なお、s-1回目のパラメータを用いたときの目的関数の値とs回目のパラメータを用いたときの目的関数の値との差が、予め定めた閾値よりも小さくなったときに収束したと判定するようにしてもよい。
【００７９】
次に、本実施の形態に係る音響信号パラメータ分析装置の作用について説明する。まず、分析対象の時系列信号として音響信号が音響信号パラメータ分析装置に入力され、図示しないメモリに格納される。そして、音響信号パラメータ分析装置において、図３に示す音響信号分析処理ルーチンが実行される。
【００８０】
まず、ステップ１００において、メモリから、あるフレーム内の音響信号を読み込み、音響信号に対して、短時間フーリエ変換を用いた時間周波数分析を行った結果から、観測時間周波数成分Ｙω,tを各（ω，ｔ）の要素にもつ二次元配列Ｙを生成して、記憶部４に記憶する。
【００８１】
そして、ステップ１０２において、乱数を用いて、パラメータθ＝｛Ｈ、Ｕ、Ｐ｝の初期値を設定して、記憶部４に記憶する。
【００８２】
次にステップ１０４では、上記ステップ１０２で設定されたパラメータθ、又は後述するステップ１０８で更新されたパラメータθと、上記ステップ１００で生成された二次元行列Ｙとに基づいて、上記（９´）式、（１０´）式に従って、補助係数＾Ｙω，ｔ，ｄを各（ω，ｔ，ｄ）の組み合わせについて算出して、記憶部４に記憶する。また、上記ステップ１０４で算出される変数μω，ｔ，ｄも、記憶部４に記憶する。
【００８３】
そして、ステップ１０６では、予め設定された基底の状態数Φ、予め設定された繰り返し回数Ｓ、及び現在の繰り返し回数ｓに基づいて、上記（１１）式に従って、スケジューリング係数ｋφを、全ての基底の状態φについて計算して、記憶部４に記憶する。
【００８４】
ステップ１０８では、上記ステップ１０２で設定されたパラメータθ、又は前回のステップ１０８で更新されたパラメータθと、上記ステップ１０４で計算された補助係数＾Ｙω，ｔ，ｄ及び変数μω，ｔ，ｄと、予め設定された行列の集合Ｇと用いて、上記（１４）式に従って、変数＾Ｐφ，ｔ，ｄを（φ，ｔ，ｄ）の組み合わせの各々について求める。そして、求められた変数＾Ｐφ，ｔ，ｄと上記ステップ１０６で計算されたスケジューリング係数ｋφとに基づいて、パラメータＰφ，ｔ，ｄの更新値を、上記（１３）式に従って、（φ，ｔ，ｄ）の組み合わせの各々について計算して、記憶部４に記憶する。
【００８５】
そして、上記ステップ１０２で設定されたパラメータＵ，Ｈ、又は前回のステップ１０８で更新されたパラメータＵ，Ｈと、上記ステップ１０４で計算された補助係数＾Ｙω，ｔ，ｄ及びμω，ｔ，ｄと、上記で計算されたパラメータＰφ，ｔ，ｄの更新値とに基づいて、上記（１５）式に従って、パラメータＨω，φ，ｄの更新値を、（ω，φ，ｄ）の組み合わせの各々について計算して、記憶部４に記憶する。
【００８６】
次に、上記ステップ１０２で設定されたパラメータＵ、又は前回のステップ１０８で更新されたパラメータＵと、上記ステップ１０４で計算された補助係数＾Ｙω，ｔ，ｄ及びμω，ｔ，ｄと、上記で計算されたパラメータＰφ，ｔ，ｄの更新値及びパラメータＨω，φ，ｄの更新値とに基づいて、上記（１６）式に従って、パラメータＵｄ，ｔの更新値を、（ｄ，ｔ）の組み合わせの各々について計算して、記憶部４に記憶する。
【００８７】
次のステップ１１０では、所定の収束条件として、繰り返し回数ｓが、Ｓに到達したか否かを判定し、繰り返し回数ｓがＳに到達していない場合には、所定の収束条件が成立していないと判断して、上記ステップ１０４へ戻り、上記ステップ１０８で更新したパラメータθを用いて、上記ステップ１０４〜ステップ１０８の処理を繰り返す。一方、繰り返し回数ｓがＳに到達した場合には、所定の収束条件が成立したと判断し、ステップ１１２で、上記ステップ１０８で最終的に更新されたパラメータθを出力して、音響信号分析処理ルーチンを終了する。
【００８８】
次に、第２の実施の形態について説明する。なお、第２の実施の形態に係る音響信号パラメータ分析装置の構成は、第１の実施の形態と同様であるため、同一符号を付して説明を省略する。
【００８９】
第２の実施の形態では、目的関数の距離尺度に、ＫＬ  ｄｉｖｅｒｇｅｎｃｅを用いている点が、第１の実施の形態と異なっている。
【００９０】
第２の実施の形態に係る音響信号パラメータ分析装置の信号分析部３の具体処理について説明する。本実施の形態では、上記（５）式の目的関数を構成する距離尺度として、Ｋullback-Leibler (KL) divergenceを利用した場合（β＝１）について説明する。
【００９１】
本実施の形態では、以下の（１７）式で表される補助関数Ｊ（＋）（θ，＾θ）を用いる。
【００９２】
【数２１】
2012027196.tif 000023
【００９３】
ここで、Ｊ（θ）≦Ｊ（＋）（θ，＾θ）における等号が成立する条件は、以下の（１８）式で表される。
【００９４】
【数２２】
2012027196.tif 000024
【００９５】
信号分析部３は、収束判定部３４で収束したと判定されるまで、補助変数更新部３１、スケジューリング係数更新部３２、パラメータ更新部３３の処理を繰り返すことにより、各パラメータＨω，φ，ｄ，Ｕｄ，ｔ，Ｐφ，ｔ，ｄの値を更新していく。以下では、繰り返し回数が第ｓ回目（ｓ＝１，２，・・・，Ｓ）の実行における各部の処理について説明をする。なお、スケジューリング係数更新部３２及び収束判定部３４の処理は、第１の実施の形態と同様であるため、説明を省略する。
【００９６】
補助変数更新部３１は、まず、記憶部４に記憶されている、前回（第ｓ−１回目）の実行で得られたパラメータＨω，φ，ｄ（ｓ−１），Ｕｄ，ｔ（ｓ−１），Ｐφ，ｔ，ｄ（ｓ−１）を用いて、以下の（１８´）式に従って、（ω，ｔ，ｄ）の組み合わせの各々について補助変数ηω，ｔ，ｄ（ｓ−１）の値を算出する。算出した値ηω，ｔ，ｄ（ｓ−１）を記憶部４に記憶しておく。
【００９７】
【数２３】
2012027196.tif 000025
【００９８】
また、パラメータ更新部３３は、計算されたスケジューリング係数ｋφ（ｓ）を用いて、以下の（１９）式に従って、第s回目の実行におけるパラメータＰφ，ｔ，ｄの更新値Ｐφ，ｔ，ｄ（ｓ）を、（φ，ｔ，ｄ）の組み合わせの各々について求める。
【００９９】
【数２４】
2012027196.tif 000026
【０１００】
ここで、＾Ｐφ，ｔ，ｄは、記憶部４に記憶されている、観測時間周波数成分Ｙω，ｔ、パラメータＨω，φ，ｄ（ｓ−１），Ｕｄ，ｔ（ｓ−１），Ｐφ，ｔ，ｄ（ｓ−１），ηω，ｔ，ｄ（ｓ−１）と、予め定められた基底の状態数Φ及び基底の数Ｄと、予め定められた行列の集合Ｇとを用いて、以下の式に従って、（φ，ｔ，ｄ）の組み合わせの各々について求められる。
【０１０１】
【数２５】
2012027196.tif 000027
【０１０２】
そして、パラメータ更新部３３は、求めた更新値Ｐφ，ｔ，ｄ（ｓ）、及び記憶部４に記憶されている補助変数ηω，ｔ，ｄ（ｓ−１），パラメータＵｄ，ｔ（ｓ−１）、観測時間周波数成分Ｙω，ｔを用いて、以下の（２０）式に従って、第s回目の実行におけるパラメータＨω，φ，ｄの更新値Ｈω，φ，ｄ（ｓ）を、（ω，φ，ｄ）の組み合わせの各々について求める。
【０１０３】
【数２６】
2012027196.tif 000028
【０１０４】
さらに、パラメータ更新部３３は、求めた更新値Ｐφ，ｔ，ｄ（ｓ）と更新値Ｈω，φ，ｄ（ｓ）、及び記憶部４に記憶されている補助変数ηω，ｔ，ｄ（ｓ−１），観測時間周波数成分Ｙω，ｔを用いて、以下の（２１）式に従って、第s回目の実行におけるパラメータＵｄ，ｔの更新値Ｕｄ，ｔ（ｓ）を、（ｄ，ｔ）の組み合わせの各々について求める。
【０１０５】
【数２７】
2012027196.tif 000029
【０１０６】
なお、第２の実施の形態に係る音響信号パラメータ分析装置の他の構成及び作用については、第１の実施の形態と同様であるため、説明を省略する。
【０１０７】
次に、第３の実施の形態について説明する。なお、第３の実施の形態に係る音響信号パラメータ分析装置の構成は、第１の実施の形態と同様であるため、同一符号を付して説明を省略する。
【０１０８】
第３の実施の形態では、目的関数の距離尺度に、板倉齊藤距離を用いている点が、第１の実施の形態と異なっている。
【０１０９】
第３の実施の形態に係る音響信号パラメータ分析装置の信号分析部３の具体処理について説明する。本実施の形態では、上記（５）式の目的関数を構成する距離尺度として、板倉齊藤距離を利用した場合（β＝０）について説明する。
【０１１０】
本実施の形態では、以下の（２２）式で表される補助関数Ｊ（＋）（θ，＾θ）を用いる。
【０１１１】
【数２８】
2012027196.tif 000030
【０１１２】
ここで、Ｊ（θ）≦Ｊ（＋）（θ，＾θ）における等号が成立する条件は、以下の（２３）式、（２４）式で表される。
【０１１３】
【数２９】
2012027196.tif 000031
【０１１４】
信号分析部３は、収束判定部３４で収束したと判定されるまで、補助変数更新部３１、スケジューリング係数更新部３２、パラメータ更新部３３の処理を繰り返すことにより、各パラメータＨω，φ，ｄ，Ｕｄ，ｔ，Ｐφ，ｔ，ｄの値を更新していく。以下では、繰り返し回数が第ｓ回目（ｓ＝１，２，・・・，Ｓ）の実行における各部の処理について説明をする。なお、スケジューリング係数更新部３２及び収束判定部３４の処理は、第１の実施の形態と同様であるため、説明を省略する。
【０１１５】
補助変数更新部３１は、まず、記憶部４に記憶されている、前回（第ｓ−１回目）の実行で得られたパラメータＨω，φ，ｄ（ｓ−１），Ｕｄ，ｔ（ｓ−１），Ｐφ，ｔ，ｄ（ｓ−１）を用いて、以下の（２３´）式に従って、（ω，ｔ）の組み合わせの各々について補助変数Ｚω，ｔ（ｓ−１）の値を算出すると共に、以下の（２４´）式に従って、（ω，ｔ，ｄ）の組み合わせの各々について補助変数ηω，ｔ，ｄ（ｓ−１）の値を算出する。算出した値Ｚω，ｔ（ｓ−１）、ηω，ｔ，ｄ（ｓ−１）を記憶部４に記憶しておく。
【０１１６】
【数３０】
2012027196.tif 000032
【０１１７】
また、パラメータ更新部３３は、計算されたスケジューリング係数ｋφ（ｓ）を用いて、以下の（２５）式に従って、第s回目の実行におけるパラメータＰφ，ｔ，ｄの更新値Ｐφ，ｔ，ｄ（ｓ）を、（φ，ｔ，ｄ）の組み合わせの各々について求める。
【０１１８】
【数３１】
2012027196.tif 000033
【０１１９】
ここで、＾Ｐφ，ｔ，ｄは、記憶部４に記憶されている、観測時間周波数成分Ｙω，ｔ、パラメータＨω，φ，ｄ（ｓ−１），Ｕｄ，ｔ（ｓ−１），Ｐφ，ｔ，ｄ（ｓ−１），Ｚω，ｔ（ｓ−１）、ηω，ｔ，ｄ（ｓ−１）と、予め定められた基底の数Ｄと、予め設定された行列の集合Ｇとを用いて、以下の式に従って、（φ，ｔ，ｄ）の組み合わせの各々について求められる。
【０１２０】
【数３２】
2012027196.tif 000034
【０１２１】
そして、パラメータ更新部３３は、求めた更新値Ｐφ，ｔ，ｄ（ｓ）、及び記憶部４に記憶されている補助変数Ｚω，ｔ（ｓ−１）、ηω，ｔ，ｄ（ｓ−１），パラメータＵｄ，ｔ（ｓ−１）、観測時間周波数成分Ｙω，ｔを用いて、以下の（２６）式に従って、第s回目の実行におけるパラメータＨω，φ，ｄの更新値Ｈω，φ，ｄ（ｓ）を、（ω，φ，ｄ）の組み合わせの各々について求める。
【０１２２】
【数３３】
2012027196.tif 000035
【０１２３】
さらに、パラメータ更新部３３は、求めた更新値Ｐφ，ｔ，ｄ（ｓ）と更新値Ｈω，φ，ｄ（ｓ）、及び記憶部４に記憶されている補助変数Ｚω，ｔ（ｓ−１）、ηω，ｔ，ｄ（ｓ−１），観測時間周波数成分Ｙω，ｔを用いて、以下の（２７）式に従って、第s回目の実行におけるパラメータＵｄ，ｔの更新値Ｕｄ，ｔ（ｓ）を、（ｄ，ｔ）の組み合わせの各々について求める。
【０１２４】
【数３４】
2012027196.tif 000036
【０１２５】
なお、第３の実施の形態に係る音響信号パラメータ分析装置の他の構成及び作用については、第１の実施の形態と同様であるため、説明を省略する。
【０１２６】
次に、楽器単音による音響信号を分析対象として、上述した第２の実施の形態に係る手法を適用し、シミュレーション実験を行った結果について説明する。
【０１２７】
ここで、収束判定部３４における収束条件として、「反復回数Ｓ＝２００回に到達したか否か」を用い、提案手法の基礎的な動作確認として、楽器単音からスペクトルの状態変化を分析した。実験に用いた信号は、ピアノＣ３の４つの音長の入った音響信号をＭＩＤＩによって合成したものとした。振幅スペクトログラムＹは、短時間フーリエ変換（サンプリング数１６ｋＨｚ，フレーム長６４ｍｓ、フレームシフト３２ｍｓ、Ｈａｎｎｉｎｇ窓）により計算した（図４（Ｂ）参照）。ピアノにおける基底状態の状態遷移はleft-to-rightに起こることを想定し、隣接行列Ａ（ｄ）は、以下の（３２）式のように予め定めた。
【０１２８】
【数３５】
2012027196.tif 000037
【０１２９】
また、基底数Ｄ＝１、状態数Φ＝４を予め定め、パラメータの初期値を乱数により与えた。信号分析を行った結果、例えば、図４（Ａ）、（Ｃ）に示すような、４つの基底の状態Ｈω，１，１と〜Ｈω，４，１とアクティベーションＵ１，ｔが得られた。得られた結果から、既定の状態遷移に着目すると、音の立ち上がりから２状態は音長によらず同程度の長さを持ち、減衰の長さが異なっている様子が確認できた。また基底として得られたスペクトルパターンに着目すると、各倍音ごとに異なる減衰をしていることも分かった。
【０１３０】
次に、提案法による音源分離の性能評価を行うために、従来のＮＭＦと提案法との間のＳＮＲ（Ｓｉｇｎａｌ−ｔｏ−Ｎｏｉｓｅ  Ｒａｔｉｏ）を比較した結果について説明する。提案法として、上述した第１の実施の形態の手法（ＥＵＣ−ｐｒｏｐｏｓｅｄ）、第２の実施の形態の手法（ＫＬ−ｐｒｏｐｏｓｅｄ）、及び第３の実施の形態の手法（ＩＳ−ｐｒｏｐｏｓｅｄ）をそれぞれ採用した。振幅スペクトログラム上で、各音源の元信号Ｉω，ｔと推定した各音源＾Ｉω，ｔとの間のＳＮＲを、以下の式により計算した。
【０１３１】
【数３６】
2012027196.tif 000038
【０１３２】
観測信号として、ヴァイオリンによるＤ♭、Ｆ、Ａ♭の混合音を用いた。ビブラートによる基底の状態遷移がleft-to-rightに起こると想定した場合と、ergodicに起こると想定して隣接行列を以下の（３３）式のように設定した場合とについて、それぞれ実験を行った。
【０１３３】
【数３７】
2012027196.tif 000039
【０１３４】
基底数Ｄ＝３、各基底の状態数Φ＝５を予め設定した。比較対象とする従来のＮＭＦでは、基底数Ｄ＝３，Ｄ＝１５として分解を行った。また、分離した各信号は元信号のどの音源に属するか分らないため、元信号の全てとのＳＮＲを求め、最も大きな値を持つ音源への割り当てを行った。もし割り当ての存在しない音源があった場合は、その音源はＳＮＲの算出には用いず、代わりに、割り当ての存在しなかった音源数の割合をエラー率として計算した。
【０１３５】
従来法及び提案法の各々において、パラメータの初期値を乱数によって与え、１０回の試行を行ったところ、図５に示すような、ＳＮＲとエラー率が得られた。また、上述した第２の実施の形態の手法を用いて音源分離した結果、図６（Ａ）に示すような、３音源を混合した観測スペクトログラムから、図６（Ｂ）、（Ｃ）に示すように、分離した１音源のスペクトログラムと、その基底の状態遷移とが得られた。上記図６（Ｃ）に示すように、基底の状態変化から、ヴァイオリンのビブラートが５つの状態遷移として学習できていることが確認できた。
【０１３６】
また、ＳＮＲに関して、第３の実施の形態の手法を除いて、提案法では、基底数が１５のときの従来法と同程度の性能を示している。従来法では、各音源がいくつかの信号に分解され、それらの統合に元信号を利用しているが、提案法においては、音源ごとに分解された信号が階層的にまとまった形で得られる。このことから、提案法は、元信号の事前知識を必要としないで音源分離を行うことのできる手法としても有効であることが確認できた。
【０１３７】
また、上記図５に示す結果から、提案法において、距離尺度として板倉齊藤距離（ＩＳ−ｐｒｏｐｏｓｅｄ）を用いた場合よりも、ユークリッド距離（ＥＵＣ−ｐｒｏｐｏｓｅｄ）とＫＬダイバージェンス（ＫＬ−ｐｒｏｐｏｓｅｄ）を距離尺度にして用いた場合の結果の方が、性能が良いことが分かった。このように、距離尺度を変えることで、より精度良い音源分離を行うことができる。
【０１３８】
次に、上述した第１の実施の形態の手法を用いた場合における目的関数の値の変化を、スケジューリング係数を用いない手法と比較した結果について説明する。状態遷移を表すＰの更新に、スケジューリング係数を導入することによる効果を検証するための実験を行った。観測信号は、上述した音源分離の性能評価のための実験と同じものを用いた。第１実施の形態の方法を用い、収束判定部３４の収束条件として「反復回数ｓ＝１００回に到達したか否か」を用いた。スケジューリング係数を導入しない場合には、ｋ１（ｓ）＝１とし、残りのｋｉ（ｓ）（ｉ＝２，・・・，Φ）を０にして計算した。つまり、上記（１３）式の代わりに、以下の式に従って、パラメータＰφ，ｔ，ｄの更新値Ｐφ，ｔ，ｄ（ｓ）を計算した。
【０１３９】
【数３８】
2012027196.tif 000040
【０１４０】
図７に示すように、スケジューリング係数を用いた場合（図７の実線を参照）と用いなかった場合（図７の破線を参照）をそれぞれ３回ずつ試行したときの目的関数の変化の様子が得られた。スケジューリング係数を導入することにより、局所解が回避できていることが確認できた。
【０１４１】
以上説明したように、本発明の実施の形態に係る音響信号パラメータ分析装置によれば、時系列信号の観測時間周波数成分Ｙω,tを要素とする二次元配列Ｙとの距離をβダイバージェンスで表わした目的関数の値が小さくなるように、各基底ｄにおける各基底の状態φを表わすパラメータＨω,φ,d、状態遷移を表わすパラメータＰφ,t,d、及びアクティベーションＵd,tを更新することにより、時系列信号から、頻出のパターンとパターンの出現の時間的な関連とを同時に精度良く抽出することができる。
【０１４２】
また、楽器単音のスペクトルが、限られた数のスペクトルパターンの状態遷移によって構成されるという階層構造を利用し、階層的にスパース表現をする枠組みを提案することにより、信号から頻出のパターンとそれらのパターンの出現の時間的な関連とを同時に抽出可能とした。
【０１４３】
また、任意の実数βについて、β-divergenceを距離尺度として用いてスペクトルの時間的な変化を分析できるようにＮＭＦを拡張することができる。これにより、入力される信号に適した距離尺度を用いることができるので、様々な信号に対して、音の特性に関する事前知識や特殊な前処理を必要とせずに、スペクトルの時間的な変化を精度良く分析することができる。
【０１４４】
また、上述した背景技術の説明において、上記の非特許文献４の技術は、距離尺度を板倉齊藤距離とした場合のＮＭＦモデルを拡張したものに相当すると述べたが、本発明の実施の形態と比較すると、拡張したモデルのパラメータを求める具体的な方法が全く異なる。また、本発明によれば、距離尺度として板倉齊藤距離以外の距離尺度を用いることも可能であると共に、非特許文献４の手法よりも早く収束させることができる、という利点がある。また、非特許文献４に記載の技術でも、ヒューリスティックな方法（非特許文献４におけるＥＭ−ＭＵ）を導入することで収束を早くしようとする試みはされているものの、そのアルゴリズムが必ず収束することが証明されていない。一方、本発明の手法によれば、必ず収束させることができ、実装も容易である。
【０１４５】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１４６】
例えば、上述の音響信号パラメータ分析装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０１４７】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【産業上の利用可能性】
【０１４８】
本発明による信号の分析は、多重音解析のための分析手法として、潜在的に、音楽音響信号からの自動採譜や、音源分離、楽器音認識、音楽加工などに応用されることが考えられる。例えば、ピアノ曲に本発明の分析を適用した場合には、各音高のアタックやサステインやリリースといった状態とそれらの出現とを同時に推定することが期待される。これは各音のオンセットについて、より正確に手掛かりを与え、自動採譜の性能を高めると考えられる。また、弦楽器や歌声によるビブラートを、ピアノ曲と同じようにいくつかの区分定常なスペクトルの状態遷移として表現出来ることが期待でき、多重音からの音源分離の性能を向上させることが可能であり、より精密な音楽加工に応用することが可能であると考えられる。
【符号の説明】
【０１４９】
１     時間周波数解析部
２     パラメータ初期値設定部
３     信号分析部
４     記憶部
５     出力部
３１   補助変数更新部
３２   スケジューリング係数更新部
３３   パラメータ更新部
３４   収束判定部
【図１】
2012027196.tif 000041
【図２】
2012027196.tif 000042
【図３】
2012027196.tif 000043
【図４】
2012027196.tif 000044
【図５】
2012027196.tif 000045
【図７】
2012027196.tif 000046
【図６】
2012027196.tif 000047