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

Description

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

一般に、音響信号などの大規模なメディアデータを用いたアプリケーションにおいて、そのデータにいかに有効なラベルを付与するか、大量のデータをいかに計算効率良く処理するかが重要な課題となっている。
従来は、教師あり学習の下で、音響信号に含まれる音イベントを自動的にラベリングする技術が一般的であった。すなわち、学習データとして、ラベル付けされた音イベントの音源信号を事前に用意し、それを周波数分析して得られる音響的特徴と統計モデルを用いて、未知の音響信号の時系列データにラベル付けを行った（例えば、非特許文献１、非特許文献２、及び非特許文献３）。統計モデルでは、例えば、ガウス混合モデル（ＧＭＭ） を用いる場合、各音イベントを表現しうる音響的特徴の頻度分布が学習される。隠れマルコフモデル（ＨＭＭ）を用いる場合、音響的特徴の統計的な時間遷移が学習され、これを用いて未知の音響信号の時系列データにラベル付けを行う。これらのラベル付けの多段処理技術（音響信号が無音なのかそうでないか、音楽か話声かを段階的に識別する）も提案されている（例えば、非特許文献４）。また、音の重ね合わせを考慮して、音響信号をあらかじめ、非負値行列分解（ＮＭＦ）によって音源分離し、分離信号を学習データに用いる技術も提案されている（例えば、非特許文献５）。

J. Saunders," Real-time discrimination of broadcast speech/music, "in Proc. ICASSP 1996. T. Butko and C. Nadeu," Audio segmentation of broadcast news in the Albayzin-2010 evaluation: overview, results, and discussion, " EURASIP Journal onAudio, Speech, and Music Processing, 2011. A. Mesaros, T. Heittola, A. Eronen and T. Virtanen," Acoustic event detection in real life recordings, "in Proc. EUSIPCO 2010. T. Butko and C. Nadeu, " Audio segmentation of broadcast news:A hierarchical system with feature selection for the Albayzin-2010 evaluation, " in Proc.ICASSP 2011. T. Heittola, A. Mesaros, T. Virtanen and A. Eronen," Sound Event Detection in Multisource Environments Using Source Separation, "in Proc. CHiME 2011.

従来は、各音イベントの音響的特徴をあらかじめラベル付けされたデータから学習する必要があった。これには、人手で構築されたラベル付音響信号データが必要であり、その構築にはコストがかかる。また、そもそも観測される音響信号に、どのような音イベントが含まれているか事前知識がないことも多く、音イベントの総数や音イベントを表現するための音響的特徴はモデルが大量の音響信号データから自動的に決定してくれることが理想的である。

本発明は、上記の事情を考慮してなされたもので、音響信号の時系列データから、そこに含まれる音イベントの数、音イベントの状態数を自動的に決定して、全ての音イベントの音響的特徴とその発音区間を自律的に推定することができる音響信号分析装置、方法、及びプログラムを提供することを目的とする。

上記の目的を達成するために本発明に係る音響信号分析装置は、複数の音イベントが混在する音響信号の時系列データを入力として、観測時間周波数成分Ｙ_ω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙ＾を出力する時間周波数分解手段と、予め定められたＤ個の音イベントｄの各々における所定の基底の状態ｋの基底スペクトルを表わすパラメータＨ_ω,d ^(k)を要素にもつ三次元配列Ｈ^、前記Ｄ個の音イベントの各々に対する各時刻ｔの発音の有無を表すパラメータＵ_d,t（＝｛０、１｝）を要素にもつ二次元配列Ｕ^、及び前記Ｄ個の音イベントｄの各々に対して各時刻ｔに前記基底の状態が何れであるかを表わすパラメータＺ_ｄ,tを要素にもつ二次元配列Ｚ^の各々の初期値を設定すると共に、Ｄ^†個（Ｄ^†＝Ｄ＋１）の音イベントｄの各々について、前記パラメータＵ_d,tが０から１へ遷移する遷移確率ａ_d（ａ₁＞ａ₂＞・・・＞ａ_D-1＞ａ_D＞ａ_D†）、及び前記パラメータＵ_d,tが１から１へ遷移する遷移確率ｂ_dの各々の初期値を設定するパラメータ初期値設定手段と、０からＵ_d,t＝１となるｔをもつ音イベントｄの中で最も小さい値となる遷移確率ａ_dまでの一様分布に従って、補助変数ｓをサンプリングにより決定し、遷移確率ａ_D†及び補助変数ｓに応じて音イベントｄを追加する場合、Ｄ、Ｄ^†を更新すると共に、追加する音イベントｄに対するパラメータＵ_d,t及び前記所定個の基底の状態ｋを表わすパラメータＨ_ω,d ^(k)の初期値を設定する補助変数推論手段と、（ω、ｔ、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^に基づいてパラメータｐ_ω,t,dを算出し、算出した各パラメータｐ_ω,t,dと、前記二次元配列Ｙ^とをパラメータとする多項分布に従って、Ｄ^†個の音イベントｄの各々に対する時間周波数成分Ｃ_ω,t,ｄを要素にもつ三次元配列Ｃ＾をサンプリングにより決定する音源分離手段と、Ｄ^†個の音イベントｄの各々について、前記二次元配列Ｚ^、前記三次元配列Ｃ^，及び前記二次元配列Ｕ^に基づいて、各時刻ｔに前記基底の状態が各状態ｋである事後確率、及び各時刻ｔに前記基底の状態が新規の状態である事後確率を算出し、算出した事後確率に従って、前記二次元配列Ｚ^をサンプリングにより決定する状態系列推論手段と、（ｋ、ω、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、及び前記三次元配列Ｃ^に基づいて、パラメータφ_ω,d ^(k)、ψ_ω,d ^(k)を算出し、算出したパラメータφ_ω,d ^(k)、ψ_ω,d ^(k)をパラメータとする確率分布に従って、前記三次元配列Ｈ^をサンプリングにより決定する基底スペクトル推論手段と、（ｄ、ｔ）の全ての組み合わせについて、前記三次元配列Ｃ^、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記遷移行列ａ_d，ｂ_dに基づいて、パラメータＵ_d,tの事後分布を算出し、算出したパラメータＵ_d,Tの事後分布に従って、Ｕ_d,Tをサンプリングにより決定し、ｔ＝Ｔ−１，・・・，１について、算出したパラメータＵ_d,ｔの事後分布と、前記遷移行列ａ_d，ｂ_dに基づく事後分布ｐ（Ｕ_d,ｔ＋１｜Ｕ_d,ｔ）との積に従って、Ｕ_d,ｔをサンプリングにより決定することにより、前記二次元配列Ｕ^を求めるアクティベーション推論手段と、前記Ｄ^†個の音イベントｄの各々について、前記遷移確率ａ_d-1、ａ_d＋1及び前記二次元配列Ｕ^に基づく遷移確率ａ_dの確率分布に従って、遷移確率ａ_dをサンプリングにより決定し、前記遷移確率ａ_D†に基づく遷移確率ａ_dの確率分布に従って、遷移確率ａ_D†をサンプリングにより決定すると共に、前記Ｄ^†個の音イベントｄの各々について、前記二次元配列Ｕ^に基づく遷移確率ｂ_dの確率分布に従って、遷移確率ｂ_dをサンプリングにより決定する遷移確率推論手段と、予め定められた収束条件を満たすまで、前記補助変数推論手段、前記音源分離手段、前記状態系列推論手段、前記基底スペクトル推論手段、前記アクティベーション推論手段、及び前記遷移確率推論手段による各処理を繰り返し行い、前記三次元配列Ｃ＾、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^を出力する収束判定手段と、を含んで構成されている。

本発明に係る音響信号分析方法は、時間周波数分解手段、パラメータ初期値設定手段、補助変数推論手段、音源分離手段、基底スペクトル推論手段、アクティベーション推論手段、遷移確率推論手段、及び収束判定手段を含む音響信号分析装置における音響信号分析方法であって、前記時間周波数分解手段によって、複数の音イベントが混在する音響信号の時系列データを入力として、観測時間周波数成分Ｙ_ω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙ＾を出力し、前記パラメータ初期値設定手段によって、予め定められたＤ個の音イベントｄの各々における所定の基底の状態ｋの基底スペクトルを表わすパラメータＨ_ω,d ^(k)を要素にもつ三次元配列Ｈ^、前記Ｄ個の音イベントの各々に対する各時刻ｔの発音の有無を表すパラメータＵ_d,t（＝｛０、１｝）を要素にもつ二次元配列Ｕ^、及び前記Ｄ個の音イベントｄの各々に対して各時刻ｔに前記基底の状態が何れであるかを表わすパラメータＺ_ｄ,tを要素にもつ二次元配列Ｚ^の各々の初期値を設定すると共に、Ｄ^†個（Ｄ^†＝Ｄ＋１）の音イベントｄの各々について、前記パラメータＵ_d,tが０から１へ遷移する遷移確率ａ_d（ａ₁＞ａ₂＞・・・＞ａ_D-1＞ａ_D＞ａ_D†）、及び前記パラメータＵ_d,tが１から１へ遷移する遷移確率ｂ_dの各々の初期値を設定し、前記補助変数推論手段によって、０からＵ_d,t＝１となるｔをもつ音イベントｄの中で最も小さい値となる遷移確率ａ_dまでの一様分布に従って、補助変数ｓをサンプリングにより決定し、遷移確率ａ_D†及び補助変数ｓに応じて音イベントｄを追加する場合、Ｄ、Ｄ^†を更新すると共に、追加する音イベントｄに対するパラメータＵ_d,t及び前記所定個の基底の状態ｋを表わすパラメータＨ_ω,d ^(k)の初期値を設定し、前記音源分離手段によって、（ω、ｔ、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^に基づいてパラメータｐ_ω,t,dを算出し、算出した各パラメータｐ_ω,t,dと、前記二次元配列Ｙ^とをパラメータとする多項分布に従って、Ｄ^†個の音イベントｄの各々に対する時間周波数成分Ｃ_ω,t,ｄを要素にもつ三次元配列Ｃ＾をサンプリングにより決定し、前記状態系列推論手段によって、Ｄ^†個の音イベントｄの各々について、前記二次元配列Ｚ^、前記三次元配列Ｃ^，及び前記二次元配列Ｕ^に基づいて、各時刻ｔに前記基底の状態が各状態ｋである事後確率、及び各時刻ｔに前記基底の状態が新規の状態である事後確率を算出し、算出した事後確率に従って、前記二次元配列Ｚ^をサンプリングにより決定するステップと、前記基底スペクトル推論手段によって、（ｋ、ω、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、及び前記三次元配列Ｃ^に基づいて、パラメータφ_ω,d ^(k)、ψ_ω,d ^(k)を算出し、算出したパラメータφ_ω,d ^(k)、ψ_ω,d ^(k)をパラメータとする確率分布に従って、前記三次元配列Ｈ^をサンプリングにより決定し、前記アクティベーション推論手段によって、（ｄ、ｔ）の全ての組み合わせについて、前記三次元配列Ｃ^、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記遷移行列ａ_d，ｂ_dに基づいて、パラメータＵ_d,tの事後分布を算出し、算出したパラメータＵ_d,Tの事後分布に従って、Ｕ_d,Tをサンプリングにより決定し、ｔ＝Ｔ−１，・・・，１について、算出したパラメータＵ_d,ｔの事後分布と、前記遷移行列ａ_d，ｂ_dに基づく事後分布ｐ（Ｕ_d,ｔ＋１｜Ｕ_d,ｔ）との積に従って、Ｕ_d,ｔをサンプリングにより決定することにより、前記二次元配列Ｕ^を求め、前記遷移確率推論手段によって、前記Ｄ^†個の音イベントｄの各々について、前記遷移確率ａ_d-1、ａ_d＋1及び前記二次元配列Ｕ^に基づく遷移確率ａ_dの確率分布に従って、遷移確率ａ_dをサンプリングにより決定し、前記遷移確率ａ_D†に基づく遷移確率ａ_dの確率分布に従って、遷移確率ａ_D†をサンプリングにより決定すると共に、前記Ｄ^†個の音イベントｄの各々について、前記二次元配列Ｕ^に基づく遷移確率ｂ_dの確率分布に従って、遷移確率ｂ_dをサンプリングにより決定し、前記収束判定手段によって、予め定められた収束条件を満たすまで、前記補助変数推論手段、前記音源分離手段、前記状態系列推論手段、前記基底スペクトル推論手段、前記アクティベーション推論手段、及び前記遷移確率推論手段による各処理を繰り返し行い、前記三次元配列Ｃ＾、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^を出力する。

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

以上説明したように、本発明の音響信号分析装置、方法、及びプログラムによれば、補助変数ｓをサンプリングにより決定して音イベントを追加し、音イベントｄの各々に対する時間周波数成分Ｃ_ω,t,ｄを要素にもつ三次元配列Ｃ＾をサンプリングにより決定し、音イベントｄの各々に対して各時刻ｔに基底の状態が何れであるか、または新規の状態であることを表わすパラメータＺ_ｄ,tを要素にもつ二次元配列Ｚ^をサンプリングにより決定し、音イベントｄの各々における基底の各状態ｋの基底スペクトルを表わすパラメータＨ_ω,d ^(k)を要素にもつ三次元配列Ｈ^をサンプリングにより決定し、音イベントｄの各々に対する各時刻ｔの発音の有無を表すパラメータＵ_d,tを要素にもつ二次元配列Ｕ^をサンプリングにより決定し、音イベントｄの各々について、パラメータＵ_d,tが０から１へ遷移する遷移確率ａ_dをサンプリングにより決定すると共にパラメータＵ_d,tが１から１へ遷移する遷移確率ｂ_dをサンプリングにより決定し、予め定められた収束条件を満たすまで、各処理を繰り返し行い、三次元配列Ｃ＾、二次元配列Ｚ^、三次元配列Ｈ^、及び二次元配列Ｕ^を出力することにより、音響信号の時系列データから、そこに含まれる音イベントの数、音イベントの状態数を自動的に決定して、全ての音イベントの音響的特徴とその発音区間を自律的に推定することができる、という効果が得られる。

音響信号への音イベントのマルチラベリングを示すイメージ図である。 音響信号のモデルパラメータの構造を示すイメージ図である。 本発明の実施の形態に係る提案モデルのグラフィカル表現図である。 本発明の実施の形態に係る音響信号分析装置の構成を示す概略図である。 本発明の実施の形態に係る音響信号分析装置における音響信号分析処理ルーチンの内容を示すフローチャートである。 本発明の実施の形態での予備実験に用いた音響信号のスペクトログラム（ピアノのＣ４、バイオリンのＥ４、オーボエのＧ４による演奏）を示す図である。 本発明の実施の形態での予備実験結果である、楽器音のスペクトルグラムへのマルチラベリングとセグメンテーションを示す図である。 本発明の実施の形態での予備実験結果である、学習された各楽器音のスペクトル集合を示す図である。 本発明の実施の形態におけるスライスサンプリングによるパラメータ更新（音イベント追加）のイメージ図である。 本発明の実施の形態におけるスライスサンプリングによるパラメータ更新（音イベント追加なし）のイメージ図である。

以下、図面を参照して本発明の実施の形態を詳細に説明する。以下では、まず本発明の実施の形態の概要及び原理について説明する。

＜発明の概要＞
本発明の実施の形態は、テレビやラジオ放送、ポッドキャスト、動画投稿サイトにおける動画などの音響信号の時系列データに対して、そこに含まれる様々な音イベント（音楽、会話、ナレーター、ベル音、ノイズなどの音のカテゴリを指す）のマルチラベリングを行うアルゴリズムに関する。具体的には、音響信号を周波数分析して得られるスペクトログラムから、そこに含まれる様々な音イベントの音響的特徴を表現するための基底となるスペクトルの集合（以降、基底スペクトル集合と呼ぶ）、各音イベントのスペクトルの状態遷移を表現する状態系列、そして各音イベントの発音区間（音が鳴っているか否か、ＯＮ／ＯＦＦの状態からなる）を表現するアクティベーション集合を抽出する音響信号分析装置、方法及びプログラムに関する。図１は、音響信号への音イベントのマルチラベリングの概略図を示す。また、図２は、観測スペクトログラムから抽出する基底スペクトル集合、状態系列、アクティベーション集合の概略図を示す。

本発明の実施の形態では、機械学習分野で注目されるノンパラメトリックベイズ法を利用して、大規模な音響信号の時系列データから、そこに含まれる全ての音イベントの音響的特徴とその発音区間を自律的に学習させるフレームワークを提供する。具体的には、音響信号に含まれる音イベントの数（図２に示すＤに相当する）が無限大の可能性を持つと仮定し、各音イベントの音源の持続性をマルコフモデルとして導入したMarkov Indian Buffet Process を用いて、各音イベントの音響的特徴およびそれらの発音区間を確率的に推定する。

＜原理＞
[１．音イベント検出のための非負値行列分解型スペクトログラムモデル]
非負値行列分解(ＮＭＦ)を音響信号に適用する場合、一般的には振幅スペクトログラム、またはパワースペクトログラムＹ＾＝（Ｙ_ω，ｔ）_Ω×T∈Ｒ^≧0,Ω×Ｔ（ただし、ω＝１,... ,Ωは周波数のインデックス、ｔ＝１,...,Ｔは時間のインデックス）を基底スペクトル集合Ｈ＾＝（Ｈ_ω，ｄ）_Ω×Ｄ∈Ｒ^≧0,Ω×Ｄと各基底スペクトルのアクティベーション集合Ｕ＾＝ （Ｕ_ｄ，ｔ）_Ｄ×T∈Ｒ^{≧0,Ｄ×Ｔ}の積で表現できるという仮定に基づいている（例えば、文献：P. Smaragdis and J. C. Brown. Nonnegativematrix factorization for polyphonic music transcription. In Proc. WASPAA 2003.や文献：T. Virtanen,“Monaural sound source separation by nonnegative matrix factorizationwith temporal continuity and sparseness criteria, ”IEEE Transactions on Audio, Speech,and Language Processing, vol. 15, pp. 1066-1074, Mar. 2007.）。すなわち、

のように観測スペクトログラムＹ＾をＤ個の頻出の基底スペクトルｈ_ｄ＾＝[Ｈ_１,ｄ,..., Ｈ_Ω,ｄ]^Ｔとそれぞれ基底スペクトルの音量変化を表すアクティベーションｕ_ｄ＾＝[Ｕ_ｄ,１,...,Ｕ_ｄ,Ｔ]で近似することに相当する。ここで、ｈ_ｄ＾とｕ_ｄ＾のペアを音イベントのコンポーネントと呼ぶ。このコンポーネントによって、一つの音イベントが表現されることが望ましいが、実際の音イベントのスペクトルは時間的に変化し、非定常であると言える。基底スペクトルが時間にともなって変化するように拡張したモデルが提案されており、そこでは各基底スペクトルは時刻ｔに、ある一つの状態Ｚ_ｄ,ｔ∈Ｎを取ると見なし、

と表現する（例えば、文献：A. Ozerov, C. F´evotte, and M. Charbit, “ Factorial scaledhidden Markov model for polyphonic audio representation and source separation, ”in Proc.WASPAA 2009. や文献：M. Nakano, J. Le Roux, H. Kameoka, N. Ono and S. Sagayama,“ Infinite-state spectrum model for music signal analysis, ”in Proc. ICASSP 2011.）。なお、記号に付された「＾」は、当該記号が行列または多次元配列またはベクトルであることを表わしている。本実施の形態はこのような基底スペクトルの状態遷移を考慮した非負値行列分解型のスペクトログラムモデルを土台とする。ただし、目的が音響信号に含まれる様々な音イベントの発音区間（ＯＮ／ＯＦＦ）の推定を含むため、アクティベーションは０（ＯＦＦ）もしくは１（ＯＮ）の値を取るものとする。すなわち、Ｕ_ｄ,ｔ∈｛０，１｝となる。したがって、これまでアクティベーションによって表現された音響信号の音量に関する情報はすべて、基底スペクトルに含まれて表現されることを想定する。
ＮＭＦは一般的に観測データとモデル間の距離尺度を用いて目的関数を設計し、それを最小化する制約付きの最適化問題として定式化される。この距離尺度の選び方は重要であり、従来さまざまな研究が行われてきた。よく用いられる尺度としては、Euclidean distance、一般化Kullback-Leibler divergence やItakura-Saito divergence などが挙げられるが、最近ではこれらを含むより広いクラスのβ-divergence が用いられることも多い（文献：D.FitzGerald,M.Cranitch and E.Coyle,“On the use of the Beta Divergence for Musical Source Separation,”in Proc. ISSC 2009.）。本実施の形態では距離尺度の選び方は中心的な話題ではないため、音源分離において性能が良いと報告されている振幅スペクトログラムに対する一般化Kullback-Leibler （ＫＬ） divergence（文献：A.T.Cemgil, “ Bayesian inference in non-negative matrixfactorisation models, ”in University of Cambridge, 2008）を用いた状況に限定して議論する。ただし、本実施の形態は距離尺度の選び方に依存したものではなく、軽微な修正によって他の尺度を用いることが可能である。観測スペクトログラムがＤ個の音イベントＣ＾＝（Ｃ_{ω，ｔ，ｄ}）_Ω×T×Ｄの重ね合わせとして表現されていると考えると、

のように書ける。上記（４）式で示されるポアソン分布の再生性より、

が成り立つ。これは前述したように、モデルと観測スペクトログラムの間の距離尺度をKLdivergence とした場合と等価である。二乗誤差規準やItakura-Saito divergence 規準の場合においては、上記（５）式 をガウス分布、複素ガウス分布とすることで同様に議論できる。

[２．スペクトログラムモデルのグラフィカル表現]
スペクトログラムモデルをベイズ的な枠組みの下でグラフィカル表現したものが図３である。スペクトログラムの生成過程を以下、順番に説明する。

[２．１ 基底スペクトルの生成（Ｈ＾の事前分布）]
基底スペクトルを生成するための事前分布として、上記（５）式で示されるポアソン分布の共役事前分布であるガンマ分布を利用する。実際、ガンマ分布の事前分布はスパースな解を導くという報告もある。そこで、

とする。ここで、ｋは基底スペクトルのインデックスであり、音イベントｄに対して合計Ｋ_ｄ個の基底スペクトルを準備する。このとき、基底スペクトル集合Ｈ＾の同時分布は

と書ける。上記（７）式で示される分布に従って、基底スペクトル集合の各要素を生成する。ここで、Θ_Ｈ＝｛φ_1,1,ψ_1,1,φ_2,1,ψ_2,1,..., φ_Ω,Ｄ,ψ_Ω,Ｄ}とする。

[２．２ アクティベーションの生成（Ｕ＾の事前分布）]
アクティベーションは各音イベントの発音区間を表現する。音源の持続性を表現するために、アクティベーションの系列のマルコフ性を仮定して、各音イベントごとに遷移行列

を用意する。つまり、０→０の遷移確率は１−ａ_ｄ，０→１の遷移確率はａ_ｄ，１→０の遷移確率は１−ｂ_ｄ，１→１ の遷移確率はｂ_ｄとし、これらの状態遷移によって、０と１からなる音イベントｄのアクティベーションを生成する。具体的には、Ｕ_ｄ，０＝０として、

のように上記（９）式で示されるベルヌーイ分布から生成される確率変数として記述できる。このとき、すべての音イベントのアクティベーション集合の同時分布ｐ（Ｕ＾｜ａ＾，ｂ＾）は

と書ける。ここで、ａ＾＝{ａ_１,...,ａ_Ｄ}、ｂ＾＝{ｂ_１,...,ｂ_Ｄ}とし、ｃ_ｄ ⁰⁰、ｃ_ｄ ⁰¹、ｃ_ｄ ¹⁰、ｃ_ｄ ¹¹はそれぞれ０→０,０→１,１→０,１→１に遷移する回数とする。さらに、ａ_ｄとｂ_ｄの事前分布をそれらの共役性から、

とおくと、Θ_Ｕ＝{α,γ,δ}に対するすべての音イベントのアクティベーションの同時分布は

となる。この分布に従って、アクティベーション集合を生成する。

[２．３ 基底スペクトルの状態系列の生成（Ｚ＾の事前分布）]
ｄ番目の音イベントの基底スペクトルの状態系列{Ｚ_ｄ,ｔ,...,Ｚ_ｄ,Ｔ}はそれぞれ離散的な値１,...,Ｋ_ｄ(状態のインデックス) を取る。

を各基底スペクトルの生起確率とすると、状態系列の同時分布は、

となる。ただし、δ（ｘ−ｋ）はクロネッカーのデルタであり、ｘ＝ｋのとき１、それ以外は０とする。ｎ_ｄ ^(ｋ)はＺ_ｄ,ｔ’＝ｋ（ｔ’＝１,...,Ｔ）を満たすｔ’の個数を表す。ここで、π_ｄ＾に対して、次のような事前分布を考える。

ただし、β_ｄは正のパラメータとする。生起確率π_ｄ＾を周辺化すると、β_ｄに対する状態系列の同時分布は

と書ける。したがって、状態系列集合の同時分布は

となる。上記（１６）式で示される分布に従って、状態系列集合を生成する。ここで、Θ_Ｚ＝{β_１,...,β_Ｄ} とする。

[２．４ モデルパラメータの同時分布]
図３のグラフィカル表現に基づいて、観測スペクトログラムとパラメータＣ＾,Ｚ＾,Ｈ＾,Ｕ＾の同時分布を書き起こすと、

となる。Θ_Ｚ,Θ_Ｈ,Θ_Ｕは超パラメータと呼び、Θ＝{Θ_Ｚ,Θ_Ｈ,Θ_Ｕ}とする。

[２．５ 観測スペクトログラムＹ＾および分離スペクトログラムＣ＾の対数尤度関数]
観測スペクトログラムＹ＾の対数尤度関数は、

となる。ここで、ｒ_ｄ,ｔ ^（ｋ）＝δ（Ｚ_ｄ,ｔ−ｋ）のようなインジケータを利用する。一方、観測スペクトログラムおよび音源分離された各音イベントのスペクトログラムからなる完全データＹ＾、Ｃ＾の対数尤度関数は

となる。

[３．音イベント数、基底スペクトルの状態数の決定方法]
２．で説明したスペクトログラムモデルを利用する上で、考えなければならない問題として、以下の２点を取り上げる。

１．観測スペクトログラムに含まれる音イベントの総数Ｄの決定方法
２．各音イベントを表現するための基底スペクトルの数（状態数）Ｋ_ｄの決定方法
従来は、音イベントの総数は事前に固定し、基底スペクトルの状態数はどの音イベントに対しても同数としたり、対象に応じて（例えば、話し声と音楽）、異なる状態数を与えることが多かった。しかし、一般的には観測される音響信号の事前知識がないことも多く、そこに含まれている音イベントの総数や必要な基底スペクトルの状態数はモデルが大量の音響信号データから自動的に決定されることが理想的である。

本実施の形態では、機械学習の分野で注目されているノンパラメトリックベイズ法を利用して、ＤとＫ_ｄの値を大規模データから自動的に決定する推論アルゴリズムを導出する。具体的には、以下の確率過程に基づいてアクティベーションＵ＾と状態系列Ｚ＾を生成する。

１．音イベント数を観測データから学習するために、Indian Buffet Process（T. L. Griffithsand Z. Ghahramani, “ Infinite latent feature models and the Indian buffet process,”In Proc. NIPS 2006.）にマルコフ性を導入したMarkov Indian Buffet Process（文献：J. V. Gael, Y. W. Teh, and Z. Ghahramani, “ The infinite factorial hidden Markov model, ”In Proc. NIPS 2008.）に基づいて、アクティベーションを生成する
２．各音イベントの基底スペクトルの状態数を観測データから学習するために、ChineseRestaurant Process（文献：Y. W. Teh, M. I. Jordan, “ Hierarchical Bayesian NonparametricModels with Applications,”in Bayesian Nonparametrics in Practice, Cambridge,UK: Cambridge University Press.）に基づいて、状態系列を生成する。ここで、Stick-breaking construction（文献：Y. W. Teh, D. G¨or¨ur, and Z. Ghahramani,“ Stick-breaking construction for the indian buffet process, ”in Proc. AISTAT 2007.）を利用して、Markov Indian Buffet Process を構成するため、遷移行列のパラメータａ＾とｂ＾の周辺化は行わない。具体的に、上記（１７）式 の観測データＹ＾とパラメータＣ＾,Ｚ＾,Ｈ＾,Ｕ＾,ａ＾,ｂ＾の同時分布は、

と修正する。ここで、ａ＾＝{ａ_１,ａ_２,...,ａ_Ｄ}は、ａ_１＞ａ_２＞...＞ａ_Ｄのように順序付け、Stick-breaking construction より、Ｄ→∞とすることで、音イベント数を観測される音響信号の時系列データから自動的に学習させる。このときａ_ｄは、

によって生成される。ここで、Ｄ^†をアクティベーション集合Ｕ＾の行数に相当し（Ｕ＾はＤ^†×Ｔの行列をイメージされたい）、ａ＾＝{ａ_１,ａ_２,...,ａ_Ｄ†}とする。一方、ｄ＜Ｄ^†となるｄは、少なくとも１回はＵ_ｄ，：において１の値を持つ（アクティブとなる）とする。ｂ＾＝{ｂ_１,...,ｂ_Ｄ†}は、前述と同様に、ｂ_ｄ〜Ｂｅｔａ（γ,δ）に従って生成する。
基底スペクトルの状態数Ｋ_ｄをデータから学習させるためのChinese restaurant processについては、後述する４．３にて具体的に説明する。

[４．スライスサンプリングによるパラメータ推論]
スライスサンプリングと動的計画法を組み合わせて、提案モデルのパラメータを推論する。これは、Stick-breaking constructionの打ち切り数（音イベント数）を予め大きな値に固定するのではなく、スライスサンプリングによって適応的に打ち切り数を選択しながら、パラメータを推論する方法である。まず、スライス補助変数ｓを導入する。

ｓが与えられたとき、Ｕ＾の条件付き確率は、

となる。ここで、I(A) は，Aが真であるときI(A)＝１となり、それ以外で０となる関数とする。この式は、ａ_ｄ＜ｓとなるすべての音イベントｄのＵ_ｄ，：を０とする。ここで、Ｄ^＊をａ_Ｄ＊＞ｓとなる最大の音イベントのインデックスとする。
アクティブな音イベントｄ（すなわち，ｄ：∃ｔ，Ｕ_ｄ，ｔ＝１）に対して、ｄ＜Ｄ^†となるインデックスＤ^†を導入する（Ｄ^†自体は非アクティブな音イベントとする）。スライスサンプリングにおいて推論すべきパラメータは、スライス補助変数ｓとＤ^†番目までの音イベントとなる。すなわち，＜ｓ，Ｄ^＊，Ｄ^†，Ｃ_{：,：,１：Ｄ†}，Ｚ_{１：Ｄ†,：}，Ｈ^（：） _{：,１：Ｄ†}，Ｕ_{１：Ｄ†,：}，ａ_１：Ｄ†，ｂ_１：Ｄ†＞であり、尤度関数の値が収束するまで、これらを繰り返しサンプリングする。各パラメータのサンプリング方法を以下に示す。

[４．１ スライス補助変数ｓのサンプリング]
上記（２２）式から、スライス補助変数ｓをサンプリングする。

（ｓ＜ａ_Ｄ†の場合）
Ｄ^＊＝Ｄ^†と更新する。Ｄ^＊＜Ｄ^†である必要があるので、

にしたがって、ｄ＝Ｄ^＊＋１についてａ_ｄをサンプリングする（音イベントを追加する）。この分布はlogａ_ｄについて対数凹分布となるため、適応的棄却サンプリング（Adaptive rejection sampling, ARS）を利用できる。ｓ＞ａ_ｄとなるまで、ａ_ｄを繰り返しサンプリングし、ｓ＞ａ_ｄとなった初めてのｄをＤ^†とするように更新する。さらに、ここで追加された音イベントｄに関しては、Ｕ_ｄ，：＝０、Ｈ^（１） _ω，ｄ〜Ｇａｍｍａ（φ_ω，ｄ，ψ_ω，ｄ）としてパラメータを初期化する。

（ｓ≧ａ_Ｄ†の場合）
ｓ＜ａ_ｄとなる最大の音イベントのインデックスｄをＤ^＊＝ｄと更新する。Ｄ^†＝Ｄ^＊＋１と更新する。Ｄ^＊＜Ｄ^†なので音イベントを追加する必要はない。また、音イベントＤ^†よりもインデックスの大きいイベントｄは削除する。

［４．２ Ｃ_{ω，ｔ，１}，... ，Ｃ_{ω，ｔ，Ｄ†}のサンプリング］
Ｃ_{ω，ｔ，１}，... ，Ｃ_{ω，ｔ，Ｄ†}の条件付き確率は，

となる。ここで、ｒ_ｄ，ｔ ^（ｋ）＝δ（Ｚ_ｄ，ｔ−ｋ）とし、ｐ_{ω，ｔ，ｄ}は

とする。上記（２５）式で示される多項分布からＣ_{ω，ｔ，１}，... ，Ｃ_{ω，ｔ，Ｄ†}をサンプリングする。

［４．３ Ｚ_ｄ，ｔのサンプリング］
４．２の推論によって、音イベントごとの音源分離スペクトログラムＣ＾が得られる。次は、分離された音イベントｄの各時刻のスペクトルＣ_：,１,ｄ,...,Ｃ_：,Ｔ,ｄが、どの基底スペクトルによって表現されるか、状態系列Ｚ_ｄ,１,...,Ｚ_ｄ,Ｔを推論する。
ｐ（Ｚ_ｄ,１,...,Ｚ_ｄ,Ｔ｜Ｃ_：,１,ｄ,...,Ｃ_：,Ｔ,ｄ）を最大化する基底スペクトルの状態系列Ｚ_ｄ,１,...,Ｚ_ｄ,Ｔを求める際、単純に、Ｚ_ｄ,１,...,Ｚ_ｄ,Ｔの全ての可能な値で評価するのは計算量的に非現実である。そこで、ｐ（Ｚ_ｄ,１,...,Ｚ_ｄ,Ｔ｜Ｃ_：,１,ｄ,...,Ｃ_：,Ｔ,ｄ）からＺ_ｄ,１,...,Ｚ_ｄ,Ｔをサンプリングすることによって、ｐ（Ｚ_ｄ,１,...,Ｚ_ｄ,Ｔ｜Ｃ_：,１,ｄ,...,Ｃ_：,Ｔ,ｄ）の値の大きなＺ_ｄ,１,...,Ｚ_ｄ,Ｔを確率的に求める。このサンプリングを効率良く実現する方法として、Gibbsサンプリングを利用する。
Gibbsサンプリングでは、ｐ（Ｚ_ｄ,１,...,Ｚ_ｄ,Ｔ｜Ｃ_：,１,ｄ,...,Ｃ_：,Ｔ,ｄ）からＺ_ｄ,１,...,Ｚ_ｄ,Ｔを同時にサンプリングするのではなく、逐次的にサンプリングする。すなわち、Ｚ_ｄ,ｔをサンプリングする際、Ｚ_ｄ,ｔ以外の{Ｚ_ｄ,１,..., Ｚ_{ｄ,ｔ−１}, Ｚ_ｄ,ｔ+１,..., Ｚ_ｄ,Ｔ}（以降Ｚ_{ｄ, ＼ｔ}＾と記述する）の値は既知とした，ｐ（Ｚ_ｄ,ｔ｜Ｚ_ｄ,＼ｔ＾,Ｃ_：,ｔ,ｄ,Θ）より、Ｚ_ｄ,ｔをサンプリングする。

Ｚ_ｄ，ｔの条件付き確率は、Ｙ＾の尤度関数を考慮して，

と書ける。ここで、

とする。上記（２７）式に従って、既存の状態ｋの各々について条件付き確率を算出し、新規の状態ｋ^newについて条件付き確率を算出し、算出された各状態の条件付き確率により、Ｚ_ｄ，ｔをサンプリングして決定する。

上式から分かるように、状態数Ｋ_ｄを無限極限としたスペクトログラムモデルにおいて、各時刻に用いられる状態に着目すると、他の時刻に多く用いられている状態ほど使われやすくなる性質がある。また、新しい状態が用いられやすくなるか否かについてはパラメータβ_ｄが影響している。このような状態系列の生成方法はChinese restaurant process (CRP)と呼ばれ、ディリクレ混合過程の一構成を与える。CRPの重要な性質である交換可能性（Exchangeability）より、Ｚ_ｄ,ｔの任意のｔの入れ替えをおこなって出現順序を変えても結果は変わらない。

［４．４ Ｈ_ω，ｄ ^（ｋ）のサンプリング］
Ｈ_ω，ｄ ^（ｋ）の条件付き確率は、

となる。上記（２９）式で示されるガンマ分布からＨ_ω，ｄ ^（ｋ）をサンプリングする。

［４．５ Ｕ_ｄ，１,...,Ｕ_ｄ，Ｔのサンプリング］
Forward-filtering backward-sampling アルゴリズムを利用して、Ｕ_ｄ，１,...,Ｕ_ｄ，Ｔを推論する。
このとき、Ｕ＾のｄ行以外の値はすべて固定する。まず、ｔ＝１,...,Ｔに対して、

を再帰的に計算する。次に、ｐ（Ｕ_ｄ，Ｔ｜Ｙ_{：，１：Ｔ}，Ｃ_{：，１：Ｔ，ｄ}，Ｚ_ｄ，：，Ｈ^（：） _：，ｄ）からＵ_ｄ，Ｔをサンプリングする。そして、ｔ＝Ｔ−１,...,１に対して、Ｕ_{ｄ，ｔ＋１}が与えられた下で、

に従って、Ｕ_ｄ，ｔを後方から順番にサンプリングする。

［４．６ ａ_ｄのサンプリング］
ｄ＝１,...,Ｄ^†−１対して、ａ_ｄの条件付き確率は、

となり、上記（３２）式で示されるベータ分布からａ_ｄをサンプリングする。ｄ＝Ｄ^†の音イベントは非アクティブとなるが、スライスサンプリングの打ち切り数のために、ａ_Ｄ†を計算する必要がある。ａ_Ｄ†の条件付き確率は，上記（２４）式をｄ＝Ｄ^†としたときであり、適応的棄却サンプリングを用いて、ａ_Ｄ†をサンプリングする。

［４．７ ｂ_ｄのサンプリング］
ｄ＝１,...,Ｄ^†に対して、ｂ_ｄの条件付き確率は、

となり、上記（３３）式で示されるベータ分布からｂ_ｄをサンプリングする。

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

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

本実施の形態に係る音響信号分析装置は、入力部１と、記憶部２と、演算部３と、出力部４とを備えている。また、演算部３は、短時間フーリエ変換部１０と、パラメータ初期値生成部１２と、スライス補助変数推論部１６と、音源分離推論部１８と、状態系列推論部２０と、基底スペクトル推論部２２と、アクティベーション推論部２４と、遷移確率推論部２６と、収束判定部２８とを備えている。また、スライス補助変数推論部１６、音源分離推論部１８、状態系列推論部２０、基底スペクトル推論部２２、アクティベーション推論部２４、遷移確率推論部２６、及び収束判定部２８は、モデルパラメータ更新部１４を構成する。この実施形態は、前述したパラメータ推定アルゴリズムを用いて信号解析を行う構成である。

入力部１には、分析する対象である音響信号の時系列が入力される。記憶部２は、入力部１に入力された音響信号の時系列を記憶する。また、記憶部２は、後述する各処理での結果等を記憶する。

短時間フーリエ変換部１０は、記憶部２に記憶された音響信号の時系列を読み出して、音響信号の振幅スペクトログラムを計算する。このとき、音響信号の時系列に対して時間フレームを設定し、時間フレーム長６４ｍｓ、時間フレームシフト長３２ｍｓとして、時間フレームごとに短時間フーリエ変換を行い、観測時間周波数成分Ｙ_ω,tを各（ω，ｔ）の要素にもつ二次元配列Ｙ＾（＝観測スペクトログラムＹ＾）を生成し、記憶部２に格納する。このとき、周波数のインデックスの最大値Ω＝５１２であり、時間のインデックスの最大値Ｔはフレームの総数に相当する。また、上記（５）式 に示すように、本実施の形態ではポアソン分布に従う確率変数として、スペクトログラムをモデル化するため、Ｙ＾の全ての要素を整数値に丸め込む。

また、記憶部２には、予め定められた超パラメータα＝４、β_１,...,β_Ｄ†＝１、γ＝５００、δ＝１、φ_1,1,ψ_1,1,φ_2,1,ψ_2,1,...,φ_Ω,Ｄ†,ψ_Ω,Ｄ†＝１からなるパラメータが記憶されている。

パラメータ初期値生成部１２は、基底スペクトル集合Ｈ＾、アクティベーション集合Ｕ＾、状態系列Ｚ＾、音イベント数Ｄ、及び遷移確率ａ＾並びにｂ＾の各パラメータの初期値を生成し、記憶部２へ格納する。
基底スペクトル集合Ｈ＾及びアクティベーション集合Ｕ＾については、記憶部２に記憶された二次元配列Ｙ＾に対し、上記（１）式に従って、通常のＮＭＦを適用し、推定されたＨ＾とＵ＾を初期値とし、記憶部２に格納する。具体的には、基底の状態ｋを１のみとし、基底スペクトル集合Ｈ＾の要素である全てのパラメータＨ_ω,d ⁽¹⁾が非負値であり、アクティベーション集合Ｕ＾の要素である全てのパラメータＵ_d,tが非負値である、という条件の下で、パラメータＨ_ω,d ⁽¹⁾とパラメータＵ_d,tとの積を、全てのｄについて足し合わせたモデルについて、二次元配列Ｙ＾とモデルとの距離を表わした目的関数の値が小さくなるように、パラメータＨ_ω,d ⁽¹⁾及びパラメータＵ_d,tの各々を推定して、基底スペクトル集合Ｈ＾及びアクティベーション集合Ｕ＾の初期値を設定する。
この部分のＮＭＦは、周知技術により実現でき、例えば、文献：A.T.Cemgil,“ Bayesian inference in non-negative matrix factorisation models, ” in University of Cambridge, 2008. や文献：M. Hoffman, D. Blei, and P. Cook,“ Bayesian nonparametric matrix factorization for recorded music, ”in Proc. ICML, 2010.で提案されるＮＭＦ手法を利用する。ただし、Ｕ＾に関しては、その平均値よりも大きい要素は１に、平均値よりも小さい要素は０に二値化して初期値とする。
状態系列Ｚ＾については、全ての要素の初期値としては１が設定され、記憶部２に格納される。
音イベント数Ｄについては、適切な初期値が設定され、記憶部２に格納される。例えば、初期値として音イベント数Ｄ＝１０とすることができる。
また、遷移確率ａ＾については、上記（２１）式に従って、ａ＾＝{ａ_１,... ,ａ_Ｄ,ａ_Ｄ†}についての初期値が生成され、そして遷移確率ｂ＾については、上記（１１）式に従って、初期値が生成され、各々の初期値が記憶部２に格納される。

出力部４は、モデルパラメータ更新部１４で求めた各パラメータを出力する。

次に、モデルパラメータ更新部１４の具体処理について説明する。

スライス補助変数推論部１６は、上記（２２）式に従って、スライスサンプリングのためのスライス補助変数ｓを生成し、これに従って、図９及び図１０に示すように、ａ_ｄ、ｂ_ｄの追加・削減をおこなって、Ｄ^†とＤ^＊を更新し、Ｕ＾のサイズを再構成する。ここでＤ^†はアクティベーションＵ＾の行数に相当する。
具体的には、記憶部２に記憶された遷移確率ａ＾、及びアクティベーション集合Ｕ＾に基づいて、上記（２２）式に従って、スライス補助変数ｓを生成し、記憶部２へ格納する。そして、生成されたスライス補助変数ｓと、遷移確率ａ_Ｄ†を比較し、ｓ＜ａ_Ｄ†の場合には、図９に示すように、上記（２４）式に従ってａ_ｄをサンプリングして、新たな音イベントを追加すると共に、Ｄ^†とＤ^＊を更新する。ｓ≧ａ_Ｄ†の場合には、図１０に示すように、ｓ＜ａ_ｄとなる最大の音イベントのインデックスｄをＤ^＊＝ｄと更新する。Ｄ^†＝Ｄ^＊＋１と更新する。Ｄ^＊＜Ｄ^†なので音イベントを追加する必要はない。また、音イベントＤ^†よりもインデックスの大きいイベントｄは削除する。
また新たに音イベントが追加された場合には、新たに追加された音イベントに対応する基底スペクトルＨ^（１） _{１,ｄ，・・・，}Ｈ^（１） _Ω,ｄ、及びアクティベーションＵ_ｄ，：も生成し、記憶部２に格納する。具体的には、追加された音イベントｄに関して、Ｕ_ｄ，：＝０、Ｈ^（１） _ω，ｄ〜Ｇａｍｍａ（φ_ω，ｄ，ψ_ω，ｄ）としてパラメータを初期化し、記憶部２に格納する。

音源分離推論部１８は、上記（２５）式及び（２６）式から、数値的サンプリングによってＣ＾のすべての要素を推論する。具体的には、記憶部２に記憶された二次元配列Ｙ＾、状態系列Ｚ＾、基底スペクトル集合Ｈ＾、及びアクティベーション集合Ｕ＾に基づいて、上記（２５）式及び（２６）式に従って、Ｃ＾のすべての要素をサンプリングして推論し、記憶部２に格納する。

状態系列推論部２０は、上記（２７）式及び（２８）式から、ギブスサンプリングによってＺ＾のすべての要素の推論を行う。具体的には、記憶部２に記憶された二次元配列Ｙ＾、超パラメータβ_１,...,β_Ｄ†、φ_1,1,ψ_1,1,φ_2,1,ψ_2,1,...,φ_Ω,Ｄ†,ψ_Ω,Ｄ†、音源分離スペクトログラムＣ＾、状態系列Ｚ＾、及びアクティベーション集合Ｕ＾に基づいて、上記（２７）式及び（２８）式に従って、ｐ（Ｚ_ｄ,１,...,Ｚ_ｄ,Ｔ｜Ｃ_・,１,ｄ,...,Ｃ_・,Ｔ,ｄ）の値の大きなＺ_ｄ,１,...,Ｚ_ｄ,Ｔを確率的に求め、記憶部２に格納する。なお、Ｚ_ｄ,ｔを求める際の状態系列Ｚ＾については、Ｚ_ｄ,ｔ以外の｛Ｚ_ｄ,１,...,Ｚ_{ｄ,ｔ−１},Ｚ_{ｄ,ｔ＋１},...,Ｚ_ｄ,ｔ｝（＝Ｚ_ｄ,＼ｔ）の値を既知として、Ｚ_ｄ,ｔをサンプリングして推論し、記憶部２に格納する。

基底スペクトル推論部２２は、上記（２９）式から、数値的サンプリングによってＨ＾のすべての要素を推論する。具体的には、記憶部２に記憶された超パラメータφ_1,1,ψ_1,1,φ_2,1,ψ_2,1,...,φ_Ω,Ｄ†,ψ_Ω,Ｄ†、音源分離スペクトログラムＣ＾、状態系列Ｚ＾、及びアクティベーション集合Ｕ＾に基づいて、上記（２９）式に従って、数値的サンプリングによってＨ＾のすべての要素をサンプリングして推論し、記憶部２に格納する。

アクティベーション推論部２４は、上記（３０）式及び（３１）式からなるForward filtering backward sampling アルゴリズムに基づいて、Ｕ＾のすべての要素を推論する。具体的には、記憶部２に記憶された二次元配列Ｙ＾、音源分離スペクトログラムＣ＾、状態系列Ｚ＾、基底スペクトル集合Ｈ＾、及びアクティベーション集合Ｕ＾に基づいて、まずForward filteringによって、上記（３０）式を再帰的に計算した後、Ｂackward samplingによって、上記（３１）式に従ってＵ_ｄ,ｔを数値的にサンプリングして求め、記憶部２に格納する。なお、Ｕ_ｄ,ｔを求める際のアクティベーション集合Ｕ＾については、Forward filteringでは、Ｕ_ｄ,ｔ-1を用い、Ｂackward samplingでは、Ｕ_ｄ,ｔ＋1を用いて、Ｕ＾の全ての要素をサンプリングして推論し、記憶部２に格納する。

遷移確率推論部２６は、記憶部２に記憶された超パラメータα、γ、δ、及びアクティベーション集合Ｕ＾に基づいて、上記（２４）式、上記（３２）式、及び上記（３３）式に従って、遷移確率ａ_ｄ及びｂ_ｄを推論する。具体的には、記憶部２に記憶された遷移確率ａ＾、及びアクティベーション集合Ｕ＾に基づいて、上記（３２）式に従って、ｄ＝１,...,Ｄ^†−１に対するａ_ｄをサンプリングして推論し、記憶部２に格納する。一方、ｄ＝Ｄ^†の場合には、超パラメータα及びａ_Ｄ†−１に基づいて、上記（２４）式に従って、ａ_Ｄ†をサンプリングして推論し、記憶部２に格納する。また、ｂ_ｄについては、記憶部２に記憶された超パラメータγ、δ、及びアクティベーション集合Ｕ＾に基づいて、上記（３３）式に従って、ｄ＝１,...,Ｄ^†に対するｂ_ｄをサンプリングにより推論し、記憶部２に格納する。

収束判定部２８は、記憶部２に記憶された、二次元配列Ｙ＾と推論されたＺ＾、Ｈ＾、Ｕ＾を用いて、上記（１８）式の対数尤度関数を計算する。更新前のモデルパラメータを用いて計算した対数尤度関数の値と更新後の値との誤差が、所定の閾値ε以下であれば、収束したと判定する。

出力部４は、収束判定部２８で収束したと判定された場合には、記憶部２に記憶されているモデルパラメータＣ＾、Ｚ＾、Ｈ＾、Ｕ＾をすべて出力する。なお、本実施形態ではこの誤差を実験的にε＝１．０×１０^−５とした。

なお、収束したか否かを判定する方法としては、対数尤度関数を用いる方法以外に、モデルパラメータ各々の値を更新前と更新後とで比較しても良いし、予め定めた繰り返し回数に到達したか否かで判定を行っても良い。本実施形態ではモデルパラメータ各々の値を更新前と更新後とで比較する場合、この誤差がε＝１．０×１０^−５であれば良好な結果であることを実験的に確認した。また、予め定めた繰り返し回数に到達したか否かで判定をする場合、１０００回の繰り返し回数が必要であることも実験的に確認している。

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

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

そして、ステップＳ１０２において、記憶部２に記憶されている超パラメータα、β_１,...,β_Ｄ†、γ、δ、φ_1,1,ψ_1,1,φ_2,1,ψ_2,1,...,φ_Ω,Ｄ†,ψ_Ω,Ｄ†（＝Θ）の値を読み込む。

次にステップＳ１０４では、パラメータ初期値生成部１２が、基底スペクトル集合Ｈ＾、アクティベーション集合Ｕ＾、状態系列Ｚ＾、音イベント数Ｄ、及び遷移確率ａ＾並びにｂ＾の各パラメータの初期値を生成する。
具体的には、上記（１）式に従って、上記ステップＳ１００において生成された二次元配列Ｙ＾に基づいて、基底スペクトル集合Ｈ＾とアクティベーション集合Ｕ＾を推定し、初期値として記憶部２へ記憶する。状態系列Ｚ＾については、全ての要素の初期値として１を設定し、音イベント数Ｄについては、初期値として例えばＤ＝１０を設定する。また遷移確率ａ＾については、上記（２１）式に従って、ａ＾＝{ａ_１,...,ａ_Ｄ†}についての初期値が生成され、遷移確率ｂ＾については、上記（１１）式に従って初期値が生成される。そして、生成された各々の初期値を記憶部２へ記憶する。

ステップＳ１０６では、スライス補助変数ｓの推論をする。具体的には、上記ステップＳ１０４で生成された遷移確率ａ＾又は後述するステップＳ１１６で前回推論された遷移確率ａ＾、及び上記ステップＳ１０４で生成されたアクティベーション集合Ｕ＾又は後述するステップＳ１１４で前回推論されたアクティベーション集合Ｕ＾に基づいて、上記（２２）式に従って、スライス補助変数ｓを生成する。そして、生成されたスライス補助変数ｓと、遷移確率ａ_Ｄ†を比較し、ｓ＜ａ_Ｄ†の場合には、上記図９に示すように、上記（２４）式に従ってａ_ｄをサンプリングして、新たな音イベントを追加すると共に、Ｄ^＊、Ｄ^†を更新する。ｓ≧ａ_Ｄ†の場合には、上記図１０に示すように、ｓ＜ａ_ｄとなる最大の音イベントのインデックスｄをＤ^＊＝ｄと更新する。Ｄ^†＝Ｄ^＊＋１と更新する。Ｄ^＊＜Ｄ^†なので音イベントを追加する必要はない。また、音イベントＤ^†よりもインデックスの大きいイベントｄは削除する。
また新たに音イベントが追加された場合には、新たに追加された音イベントに対応する基底スペクトルＨ^（１） _{１,ｄ，・・・，}Ｈ^（１） _Ω,ｄ、及びアクティベーションＵ_ｄ，：も生成する。具体的には、追加された音イベントｄに関して、Ｕ_ｄ，：＝０、Ｈ^（１） _ω，ｄ〜Ｇａｍｍａ（φ_ω，ｄ，ψ_ω，ｄ）としてパラメータを初期化する。

ステップＳ１０８では、音源分離スペクトログラムＣ＾の全ての要素を推論する。具体的には、上記ステップＳ１００において生成された二次元配列Ｙ＾、及び上記ステップＳ１０４で生成された状態系列Ｚ＾、基底スペクトル集合Ｈ＾並びにアクティベーション集合Ｕ＾、又は、後述するステップＳ１１０〜ステップＳ１１４で前回推論された状態系列Ｚ＾、基底スペクトル集合Ｈ＾並びにアクティベーション集合Ｕ＾に基づいて、上記（２５）式及び（２６）式に従って、Ｃ＾のすべての要素をサンプリングして推論し、記憶部２へ記憶する。

ステップＳ１１０では、ギブスサンプリングによって、状態系列Ｚ＾のすべての要素の推論を行う。具体的には、上記ステップＳ１００で生成された二次元配列Ｙ＾、上記ステップＳ１０２で読み込まれた超パラメータβ_１,...,β_Ｄ†、φ_1,1,ψ_1,1,φ_2,1,ψ_2,1,...,φ_Ω,Ｄ†,ψ_Ω,Ｄ†、上記ステップＳ１０４で生成された状態系列Ｚ＾並びにアクティベーション集合Ｕ＾、又は本ステップＳ１１０及び後述するステップＳ１１４で前回推論された状態系列Ｚ＾並びにアクティベーション集合Ｕ＾、及び上記ステップＳ１０８で推論された音源分離スペクトログラムＣ＾に基づいて、上記（２７）式及び（２８）式に従って、Ｚ_ｄ,１,...,Ｚ_ｄ,Ｔをサンプリングして推論する。
なお、Ｚ_ｄ,ｔを求める際の状態系列Ｚ＾については、Ｚ_ｄ,ｔ以外の｛Ｚ_ｄ,１,...,Ｚ_{ｄ,ｔ−１},Ｚ_{ｄ,ｔ＋１},...,Ｚ_ｄ,ｔ｝（＝Ｚ_ｄ,＼ｔ）の値を既知として、Ｚ_ｄ,ｔをサンプリングして推論する。ここで、Ｚ_ｄ,１,...,Ｚ_{ｄ,ｔ−１}までは、本ステップＳ１１０で推論された直前の値を用い、Ｚ_{ｄ,ｔ＋１},...,Ｚ_ｄ,ｔについては、上記ステップＳ１０４で生成された状態系列Ｚ＾又は本ステップＳ１１０で前回推論された状態系列Ｚ＾を用いて、Ｚ_ｄ,ｔをサンプリングして推論し、記憶部２へ記憶する。

ステップＳ１１２では、基底スペクトル集合Ｈ＾の全ての要素の推論を行う。具体的には、上記ステップＳ１０２で読み込まれた超パラメータφ_1,1,ψ_1,1,φ_2,1,ψ_2,1,...,φ_Ω,Ｄ†,ψ_Ω,Ｄ†、上記ステップＳ１０４で生成されたアクティベーション集合Ｕ＾、又は後述するステップＳ１１２で前回推論されたアクティベーション集合Ｕ＾、上記ステップＳ１０８で推論された音源分離スペクトログラムＣ＾、及びステップＳ１１０で推論された状態系列Ｚ＾に基づいて、上記（２９）式に従って、Ｈ＾のすべての要素をサンプリングして推論し、記憶部２へ記憶する。

ステップＳ１１４では、アクティベーション集合Ｕ＾の全ての要素の推論を行う。具体的には、上記ステップＳ１００において生成された二次元配列Ｙ＾、上記ステップＳ１０６で推論された音源分離スペクトログラムＣ＾、上記ステップＳ１１０で推論された状態系列Ｚ＾、上記ステップＳ１１２で推論された基底スペクトル集合Ｈ＾、及び上記ステップＳ１０４で生成されたアクティベーション集合Ｕ＾又は本ステップＳ１１４で前回推論されたアクティベーション集合Ｕ＾に基づいて、上記（３０）式及び（３１）式に従って、Ｕ＾の全ての要素の推論を行う。詳細には、まずForward filteringによって、上記（３０）式を再帰的に計算した後、Ｂackward samplingによって、上記（３１）式に従ってＵ_ｄ,ｔを数値的にサンプリングして求めることにより、Ｕ＾の全ての要素を推論し、記憶部２へ記憶する。

ステップＳ１１６では、遷移確率ａ_ｄ、ｂ_ｄの推論を行う。具体的には、本ステップＳ１０４で今回または前回生成された遷移確率ａ、上記ステップＳ１０２で読み込まれた超パラメータα、γ、δ、及び上記ステップＳ１１４で推論されたＵ_ｄ,１,..., Ｕ_ｄ,Ｔに基づいて、上記（２４）式、（３２）式、及び（３３）式に従って、遷移確率ａ_ｄ、ｂ_ｄの推論を行う。まず、上記（３２）式に基づいて、ｄ＝１,...,Ｄ^†−１に対するａ_ｄを推論する。一方、ｄ＝Ｄ^†の場合には、上記（２４）式に従って、ａ_Ｄ†をサンプリングして推論する。また、ｂ_ｄについては、上記（３３）式に基づいて、ｄ＝１,...,Ｄ^†に対するｂ_ｄを推論する。そして、推論された遷移確率ａ_ｄ、ｂ_ｄを記憶部２へ記憶する。

ステップＳ１１８では、二次元配列Ｙ＾と推論されたＺ＾、Ｈ＾、Ｕ＾を用いて、対数尤度関数ｌｏｇｐ（Ｙ＾｜Ｚ＾，Ｈ＾，Ｕ＾，Θ）を計算する。そして、更新前のモデルパラメータを用いて計算した対数尤度関数の値と更新後の値との誤差が、所定の閾値以下であれば、収束したと判定する。収束していないと判定された場合には、ステップＳ１０６へ移行し、上記ステップＳ１０６〜ステップＳ１１６で推論したパラメータを用いて、上記ステップＳ１０６〜ステップＳ１１６の処理を繰り返す。収束したと判定された場合には、ステップＳ１２０へ移行する。

ステップＳ１２０では、推論されたパラメータが結果として出力部４に出力され、音響信号分析処理ルーチンが終了する。

＜実験結果＞
予備実験として、オーボエ（Ｇ４、ソ音）、バイオリン（Ｅ４、ミ音）、ピアノ（Ｃ４、ド音）の音が混合された音響信号に対して、上記のアルゴリズムの適用を試みた。図６は、サンプリング周波数１６ｋＨｚ、量子化ビット数１６の音響信号をフレーム長６４ｍｓ、フレームシフト長３２ｍｓで周波数分析したときの振幅スペクトログラムである。パラメータの初期化方法については前述したとおりである。図７にマルチラベリング、およびセグメンテーションの結果を示す。この図から、オーボエ、ピアノ、バイオリンにうまく分離して特徴を捉えていることがわかり、本アルゴリズムの有効性が期待できる。図８は学習された各楽器音のスペクトル集合を示す。各楽器の音を表現しうるスペクトル数（バイオリンは６個、オーボエは５個、ピアノは７個）が、音響信号データから推定される。それぞれのスペクトル形状が各楽器に対応していることを定性的に確認した。

以上説明したように、本発明の実施の形態に係る音響信号分析装置によれば、補助変数ｓをサンプリングにより決定して音イベントを追加し、音イベントｄの各々に対する時間周波数成分Ｃ_ω,t,ｄを要素にもつ三次元配列Ｃ＾をサンプリングにより決定し、音イベントｄの各々に対して各時刻ｔに基底の状態が何れであるか、または新規の状態であることを表わすパラメータＺ_ｄ,tを要素にもつ二次元配列Ｚ^をサンプリングにより決定し、音イベントｄの各々における基底の各状態ｋの基底スペクトルを表わすパラメータＨ_ω,d ^(k)を要素にもつ三次元配列Ｈ^をサンプリングにより決定し、音イベントｄの各々に対する各時刻ｔの発音の有無を表すパラメータＵ_d,tを要素にもつ二次元配列Ｕ^をサンプリングにより決定し、音イベントｄの各々について、パラメータＵ_d,tが０から１へ遷移する遷移確率ａ_dをサンプリングにより決定すると共にパラメータＵ_d,tが１から１へ遷移する遷移確率ｂ_dをサンプリングにより決定し、予め定められた収束条件を満たすまで、各処理を繰り返し行い、三次元配列Ｃ＾、二次元配列Ｚ^、三次元配列Ｈ^、及び二次元配列Ｕ^を出力することにより、音響信号の時系列データから、そこに含まれる音イベントの数、音イベントの状態数を自動的に決定して、全ての音イベントの音響的特徴とその発音区間を自律的に推定することができる。

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

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

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

１ 入力部
３ 演算部
４ 出力部
１０ 短時間フーリエ変換部
１２ パラメータ初期値生成部
１４ モデルパラメータ更新部
１６ スライス補助変数推論部
１８ 音源分離推論部
２０ 状態系列推論部
２２ 基底スペクトル推論部
２４ アクティベーション推論部
２６ 遷移確率推論部
２８ 収束判定部

Claims (5)

1. 複数の音イベントが混在する音響信号の時系列データを入力として、観測時間周波数成分Ｙ_ω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙ＾を出力する時間周波数分解手段と、
   予め定められたＤ個の音イベントｄの各々における所定の基底の状態ｋの基底スペクトルを表わすパラメータＨ_ω,d ^(k)を要素にもつ三次元配列Ｈ^、前記Ｄ個の音イベントの各々に対する各時刻ｔの発音の有無を表すパラメータＵ_d,t（＝｛０、１｝）を要素にもつ二次元配列Ｕ^、及び前記Ｄ個の音イベントｄの各々に対して各時刻ｔに前記基底の状態が何れであるかを表わすパラメータＺ_ｄ,tを要素にもつ二次元配列Ｚ^の各々の初期値を設定すると共に、Ｄ^†個（Ｄ^†＝Ｄ＋１）の音イベントｄの各々について、前記パラメータＵ_d,tが０から１へ遷移する遷移確率ａ_d（ａ₁＞ａ₂＞・・・＞ａ_D-1＞ａ_D＞ａ_D†）、及び前記パラメータＵ_d,tが１から１へ遷移する遷移確率ｂ_dの各々の初期値を設定するパラメータ初期値設定手段と、
   ０からＵ_d,t＝１となるｔをもつ音イベントｄの中で最も小さい値となる遷移確率ａ_dまでの一様分布に従って、補助変数ｓをサンプリングにより決定し、遷移確率ａ_D†及び補助変数ｓに応じて音イベントｄを追加する場合、Ｄ、Ｄ^†を更新すると共に、追加する音イベントｄに対するパラメータＵ_d,t及び前記所定個の基底の状態ｋを表わすパラメータＨ_ω,d ^(k)の初期値を設定する補助変数推論手段と、
   （ω、ｔ、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^に基づいてパラメータｐ_ω,t,dを算出し、算出した各パラメータｐ_ω,t,dと、前記二次元配列Ｙ^とをパラメータとする多項分布に従って、Ｄ^†個の音イベントｄの各々に対する時間周波数成分Ｃ_ω,t,ｄを要素にもつ三次元配列Ｃ＾をサンプリングにより決定する音源分離手段と、
   Ｄ^†個の音イベントｄの各々について、前記二次元配列Ｚ^、前記三次元配列Ｃ^，及び前記二次元配列Ｕ^に基づいて、各時刻ｔに前記基底の状態が各状態ｋである事後確率、及び各時刻ｔに前記基底の状態が新規の状態である事後確率を算出し、算出した事後確率に従って、前記二次元配列Ｚ^をサンプリングにより決定する状態系列推論手段と、
   （ｋ、ω、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、及び前記三次元配列Ｃ^に基づいて、パラメータφ_ω,d ^(k)、ψ_ω,d ^(k)を算出し、算出したパラメータφ_ω,d ^(k)、ψ_ω,d ^(k)をパラメータとする確率分布に従って、前記三次元配列Ｈ^をサンプリングにより決定する基底スペクトル推論手段と、
   （ｄ、ｔ）の全ての組み合わせについて、前記三次元配列Ｃ^、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記遷移行列ａ_d，ｂ_dに基づいて、パラメータＵ_d,tの事後分布を算出し、算出したパラメータＵ_d,Tの事後分布に従って、Ｕ_d,Tをサンプリングにより決定し、ｔ＝Ｔ−１，・・・，１について、算出したパラメータＵ_d,ｔの事後分布と、前記遷移行列ａ_d，ｂ_dに基づく事後分布ｐ（Ｕ_d,ｔ＋１｜Ｕ_d,ｔ）との積に従って、Ｕ_d,ｔをサンプリングにより決定することにより、前記二次元配列Ｕ^を求めるアクティベーション推論手段と、
   前記Ｄ^†個の音イベントｄの各々について、前記遷移確率ａ_d-1、ａ_d＋1及び前記二次元配列Ｕ^に基づく遷移確率ａ_dの確率分布に従って、遷移確率ａ_dをサンプリングにより決定し、前記遷移確率ａ_D†に基づく遷移確率ａ_dの確率分布に従って、遷移確率ａ_D†をサンプリングにより決定すると共に、前記Ｄ^†個の音イベントｄの各々について、前記二次元配列Ｕ^に基づく遷移確率ｂ_dの確率分布に従って、遷移確率ｂ_dをサンプリングにより決定する遷移確率推論手段と、
   予め定められた収束条件を満たすまで、前記補助変数推論手段、前記音源分離手段、前記状態系列推論手段、前記基底スペクトル推論手段、前記アクティベーション推論手段、及び前記遷移確率推論手段による各処理を繰り返し行い、前記三次元配列Ｃ＾、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^を出力する収束判定手段と、
   を含む音響信号分析装置。
2. 前記補助変数推論手段は、０からＵ_d,t＝１となるｔをもつ音イベントｄの中で最も小さい値となる遷移確率ａ_dまでの一様分布に従って、補助変数ｓをサンプリングにより決定し、
   前記遷移確率ａ_D†が、前記補助変数ｓより大きい場合、追加する音イベントｄの遷移確率ａ_dを、前記遷移確率ａ_D†に基づく遷移確率ａ_dの確率分布に従ってサンプリングにより決定すると共に、前記追加する音イベントｄに対するパラメータＵ_d,t及び前記所定個の基底の状態ｋを表わすパラメータＨ_ω,d ^(k)の初期値を設定し、ＤをＤ^†に更新すると共に、サンプリングにより決定され、かつ、前記補助変数ｓより小さくなる遷移確率ａ_ｄにおけるｄに、Ｄ^†を更新し、
   前記遷移確率ａ_D†が、前記補助変数ｓ以下である場合、ｓ＜ａ_ｄとなる最大の音イベントのインデックスｄをＤ^＊＝ｄと更新し、Ｄ^†＝Ｄ^＊＋１と更新し、音イベントＤ^†よりもインデックスの大きいイベントｄは削除する請求項１記載の音響信号分析装置。
3. 初期値設定手段は、
   前記予め定められたＤ個の音イベントｄの各々における基底の状態数を１とし、全てのパラメータＨ_ω,d ⁽¹⁾が非負値であり、全てのパラメータＵ_d,tが非負値である、という条件の下で、前記パラメータＨ_ω,d ⁽¹⁾と前記パラメータＵ_d,tとの積を、全てのｄについて足し合わせたモデルについて、前記時間周波数分解手段によって出力された観測時間周波数成分Ｙ_ω,tと前記モデルとの距離を表わした目的関数の値が小さくなるように、前記パラメータＨ_ω,d ⁽¹⁾及び前記パラメータＵ_d,tの各々を推定し、前記推定された前記パラメータＵ_d,tの各々を二値化することにより、前記三次元配列Ｈ^及び前記二次元配列Ｕ^の初期値を設定すると共に、前記パラメータＺ_ｄ,tの各々を１とした二次元配列Ｚ^の初期値を設定する請求項１又は２記載の音響信号分析装置。
4. 時間周波数分解手段、パラメータ初期値設定手段、補助変数推論手段、音源分離手段、状態系列推論手段、基底スペクトル推論手段、アクティベーション推論手段、遷移確率推論手段、及び収束判定手段を含む音響信号分析装置における音響信号分析方法であって、
   前記音響信号分析装置は、
   前記時間周波数分解手段によって、複数の音イベントが混在する音響信号の時系列データを入力として、観測時間周波数成分Ｙ_ω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙ＾を出力するステップと、
   前記パラメータ初期値設定手段によって、予め定められたＤ個の音イベントｄの各々における所定の基底の状態ｋの基底スペクトルを表わすパラメータＨ_ω,d ^(k)を要素にもつ三次元配列Ｈ^、前記Ｄ個の音イベントの各々に対する各時刻ｔの発音の有無を表すパラメータＵ_d,t（＝｛０、１｝）を要素にもつ二次元配列Ｕ^、及び前記Ｄ個の音イベントｄの各々に対して各時刻ｔに前記基底の状態が何れであるかを表わすパラメータＺ_ｄ,tを要素にもつ二次元配列Ｚ^の各々の初期値を設定すると共に、Ｄ^†個（Ｄ^†＝Ｄ＋１）の音イベントｄの各々について、前記パラメータＵ_d,tが０から１へ遷移する遷移確率ａ_d（ａ₁＞ａ₂＞・・・＞ａ_D-1＞ａ_D＞ａ_D†）、及び前記パラメータＵ_d,tが１から１へ遷移する遷移確率ｂ_dの各々の初期値を設定するステップと、
   前記補助変数推論手段によって、０からＵ_d,t＝１となるｔをもつ音イベントｄの中で最も小さい値となる遷移確率ａ_dまでの一様分布に従って、補助変数ｓをサンプリングにより決定し、遷移確率ａ_D†及び補助変数ｓに応じて音イベントｄを追加する場合、Ｄ、Ｄ^†を更新すると共に、追加する音イベントｄに対するパラメータＵ_d,t及び前記所定個の基底の状態ｋを表わすパラメータＨ_ω,d ^(k)の初期値を設定するステップと、
   前記音源分離手段によって、（ω、ｔ、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^に基づいてパラメータｐ_ω,t,dを算出し、算出した各パラメータｐ_ω,t,dと、前記二次元配列Ｙ^とをパラメータとする多項分布に従って、Ｄ^†個の音イベントｄの各々に対する時間周波数成分Ｃ_ω,t,ｄを要素にもつ三次元配列Ｃ＾をサンプリングにより決定するステップと、
   前記状態系列推論手段によって、Ｄ^†個の音イベントｄの各々について、前記二次元配列Ｚ^、前記三次元配列Ｃ^，及び前記二次元配列Ｕ^に基づいて、各時刻ｔに前記基底の状態が各状態ｋである事後確率、及び各時刻ｔに前記基底の状態が新規の状態である事後確率を算出し、算出した事後確率に従って、前記二次元配列Ｚ^をサンプリングにより決定するステップと、
   前記基底スペクトル推論手段によって、（ｋ、ω、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、及び前記三次元配列Ｃ^に基づいて、パラメータφ_ω,d ^(k)、ψ_ω,d ^(k)を算出し、算出したパラメータφ_ω,d ^(k)、ψ_ω,d ^(k)をパラメータとする確率分布に従って、前記三次元配列Ｈ^をサンプリングにより決定するステップと、
   前記アクティベーション推論手段によって、（ｄ、ｔ）の全ての組み合わせについて、前記三次元配列Ｃ^、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記遷移行列ａ_d，ｂ_dに基づいて、パラメータＵ_d,tの事後分布を算出し、算出したパラメータＵ_d,Tの事後分布に従って、Ｕ_d,Tをサンプリングにより決定し、ｔ＝Ｔ−１，・・・，１について、算出したパラメータＵ_d,ｔの事後分布と、前記遷移行列ａ_d，ｂ_dに基づく事後分布ｐ（Ｕ_d,ｔ＋１｜Ｕ_d,ｔ）との積に従って、Ｕ_d,ｔをサンプリングにより決定することにより、前記二次元配列Ｕ^を求めるステップと、
   前記遷移確率推論手段によって、前記Ｄ^†個の音イベントｄの各々について、前記遷移確率ａ_d-1、ａ_d＋1及び前記二次元配列Ｕ^に基づく遷移確率ａ_dの確率分布に従って、遷移確率ａ_dをサンプリングにより決定し、前記遷移確率ａ_D†に基づく遷移確率ａ_dの確率分布に従って、遷移確率ａ_D†をサンプリングにより決定すると共に、前記Ｄ^†個の音イベントｄの各々について、前記二次元配列Ｕ^に基づく遷移確率ｂ_dの確率分布に従って、遷移確率ｂ_dをサンプリングにより決定するステップと、
   前記収束判定手段によって、予め定められた収束条件を満たすまで、前記補助変数推論手段、前記音源分離手段、前記状態系列推論手段、前記基底スペクトル推論手段、前記アクティベーション推論手段、及び前記遷移確率推論手段による各処理を繰り返し行い、前記三次元配列Ｃ＾、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^を出力するステップと、
   を含む音響信号分析方法。
5. コンピュータを、請求項１〜請求項３の何れか１項記載の音響信号分析装置の各手段として機能させるためのプログラム。