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

Abstract

【課題】音響信号の時系列データから、そこに含まれる音響イベントが、どのような音響要素で構成されているのかを精度よく推定することができるようにする。
【解決手段】音響特徴量成分Ｃ＾を決定し、各時刻ｔに基底の状態が何れであるか、または新規の状態であることを表わすＺ^を決定し、基底の各状態ｋの基底スペクトルを表わすＨ^を決定し、予め与えられたＬ個の音響イベントｌの各々に対する一部の時刻ｔnの各々の音響イベントの有無をＸ_ｌＮ^に基づいて、各時刻ｔの発音の有無を表すＵ^を決定し、遷移確率ａ_（d）を決定すると共に遷移確率ｂ_（d）を決定し、Ｌ個の音響イベントｌの各々に対するＤ個の音響要素ｄの各々の重みｗ_ｌ＾を求め、Ｌ個の音響イベントｌの各々に対する精度パラメータα_ｌを算出し、予め定められた収束条件を満たすまで、各処理を繰り返し行い、Ｃ＾、Ｚ^、Ｈ^、Ｕ^及びｗ_ｌ＾を出力する。
【選択図】図４

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 on Audio, 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.2 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,tが１から１へ遷移する遷移確率ｂ_dの各々の初期値を設定するステップと、前記音源分離手段によって、（ω、ｔ、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^に基づいてパラメータｐ_ω,t,dを算出し、算出した各パラメータｐ_ω,t,dと、前記二次元配列Ｙ^とをパラメータとする多項分布に従って、Ｄ個の音響要素ｄの各々に対する音響特徴量成分Ｃ_ω,t,ｄを要素にもつ三次元配列Ｃ＾をサンプリングにより決定するステップと、前記状態系列推論手段によって、前記Ｄ個の音響要素ｄの各々について、前記二次元配列Ｚ^、前記三次元配列Ｃ^，及び前記二次元配列Ｕ^に基づいて、各時刻ｔに前記基底の状態が各状態ｋである事後確率、及び各時刻ｔに前記基底の状態が新規の状態である事後確率を算出し、算出した事後確率に従って、前記二次元配列Ｚ^をサンプリングにより決定するステップと、前記基底スペクトル推論手段によって、（ｋ、ω、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^及び前記三次元配列Ｃ^に基づいて、パラメータφ_ω,d ^(k)、ψ_ω,d ^(k)を算出し、算出したパラメータφ_ω,d ^(k)、ψ_ω,d ^(k)をパラメータとする確率分布に従って、前記三次元配列Ｈ^をサンプリングにより決定するステップと、前記アクティベーション推論手段によって、（ｄ、ｔ）の全ての組み合わせについて、予め与えられた前記Ｌ個の音響イベントｌの各々に対する一部の時刻ｔnの各々の音響イベントの有無を表すパラメータＸ_ｌ,tｎ（＝｛０、１｝）を要素にもつ一次元配列Ｘ_ｌＮ^、前記三次元配列Ｃ^、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記遷移行列ａ_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をサンプリングにより決定するステップと、前記重みベクトル推論手段によって、前記Ｌ個の音響イベントｌの各々について、前記一次元配列Ｘ_ｌＮ^、前記二次元配列Ｕ^、前記一次元配列ｗ_ｌ＾、及び前記パラメータα_ｌに基づいて、事後確率が最大となるように前記一次元配列ｗ_ｌ＾を求めるステップと、前記精度パラメータ推論手段によって、前記Ｌ個の音響イベントｌの各々について、前記一次元配列ｗ_ｌ＾に基づく前記パラメータα_lの確率分布に従って、前記パラメータα_ｌをサンプリングにより決定するステップと、前記収束判定手段によって、予め定められた収束条件を満たすまで、前記音源分離手段、前記状態系列推論手段、前記基底スペクトル推論手段、前記アクティベーション推論手段、前記遷移確率推論手段、前記重みベクトル推論手段、及び前記精度パラメータ推論手段による各処理を繰り返し行い、前記三次元配列Ｃ＾、前記二次元配列Ｚ^、前記三次元配列Ｈ^、前記二次元配列Ｕ^、及び前記Ｌ個の音響イベントｌの各々の前記一次元配列ｗ_ｌ＾を出力するステップと、を含む。

本発明に係る音響信号分析装置は、 複数の音響要素が混在する音響信号の時系列データを入力として、音響特徴量Ｙ_ω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙ＾を出力する音響特徴量抽出手段と、予め定められたＤ個の音響要素ｄの各々における所定の基底の状態ｋの基底スペクトルを表わすパラメータＨ_ω,d ^(k)を要素にもつ三次元配列Ｈ^、前記Ｄ個の音響要素の各々に対する各時刻ｔの発音の有無を表すパラメータＵ_d,t（＝｛０、１｝）を要素にもつ二次元配列Ｕ^、及び前記Ｄ個の音響要素ｄの各々に対して各時刻ｔに前記基底の状態が何れであるかを表わすパラメータＺ_ｄ,tを要素にもつ二次元配列Ｚ^、予め定められたＬ個の音響イベントｌの各々に対する前記Ｄ個の音響要素ｄの各々の重みｗ_ｌ，ｄを要素にもつ一次元配列ｗ_ｌ＾、前記Ｌ個の音響イベントｌの各々に対する精度パラメータα_ｌの各々の初期値を設定すると共に、前記Ｄ個の音響要素ｄの各々について、前記パラメータＵ_d,tが０から１へ遷移する遷移確率ａ_d（ａ₁＞ａ₂＞・・・＞ａ_D-1＞ａ_D）、及び前記パラメータＵ_d,tが１から１へ遷移する遷移確率ｂ_dの各々の初期値を設定するパラメータ初期値設定手段と、（ω、ｔ、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^に基づいてパラメータｐ_ω,t,dを算出し、算出した各パラメータｐ_ω,t,dと、前記二次元配列Ｙ^とをパラメータとする多項分布に従って、Ｄ個の音響要素ｄの各々に対する音響特徴量成分Ｃ_ω,t,ｄを要素にもつ三次元配列Ｃ＾をサンプリングにより決定する音源分離手段と、前記Ｄ個の音響要素ｄの各々について、前記二次元配列Ｚ^、前記三次元配列Ｃ^，及び前記二次元配列Ｕ^に基づいて、各時刻ｔに前記基底の状態が各状態ｋである事後確率、及び各時刻ｔに前記基底の状態が新規の状態である事後確率を算出し、算出した事後確率に従って、前記二次元配列Ｚ^をサンプリングにより決定する状態系列推論手段と、（ｋ、ω、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^及び前記三次元配列Ｃ^に基づいて、パラメータφ_ω,d ^(k)、ψ_ω,d ^(k)を算出し、算出したパラメータφ_ω,d ^(k)、ψ_ω,d ^(k)をパラメータとする確率分布に従って、前記三次元配列Ｈ^をサンプリングにより決定する基底スペクトル推論手段と、（ｄ、ｔ）の全ての組み合わせについて、予め与えられた前記Ｌ個の音響イベントｌの各々に対する一部の時刻ｔnの各々の音響イベントの有無を表すパラメータＸ_ｌ,tｎ（＝｛０、１｝）を要素にもつ一次元配列Ｘ_ｌＮ^、前記三次元配列Ｃ^、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記遷移行列ａ_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をサンプリングにより決定する遷移確率推論手段と、前記Ｌ個の音響イベントｌの各々について、前記一次元配列Ｘ_ｌＮ^、前記二次元配列Ｕ^、前記一次元配列ｗ_ｌ＾、及び前記パラメータα_ｌに基づいて、事後確率が最大となるように前記一次元配列ｗ_ｌ＾を求める重みベクトル推論手段と、前記Ｌ個の音響イベントｌの各々について、前記一次元配列ｗ_ｌ＾に基づく前記パラメータα_lの確率分布に従って、前記パラメータα_ｌをサンプリングにより決定する精度パラメータ推論手段と、予め定められた収束条件を満たすまで、前記音源分離手段、前記状態系列推論手段、前記基底スペクトル推論手段、前記アクティベーション推論手段、前記遷移確率推論手段、前記重みベクトル推論手段、及び前記精度パラメータ推論手段による各処理を繰り返し行い、前記三次元配列Ｃ＾、前記二次元配列Ｚ^、前記三次元配列Ｈ^、前記二次元配列Ｕ^、及び前記Ｌ個の音響イベントｌの各々の前記一次元配列ｗ_ｌ＾を出力する収束判定手段と、を含んで構成されている。

本発明に係るプログラムは、上記の音響信号分析方法を構成する各ステップをコンピュータに実行させるためのプログラムである。

以上説明したように、本発明の音響信号分析方法、装置、及びプログラムによれば、音響要素ｄの各々に対する音響特徴量成分Ｃ_ω,t,ｄを要素にもつＣ＾をサンプリングにより決定し、音響要素ｄの各々に対して各時刻ｔに基底の状態が何れであるか、または新規の状態であることを表わすパラメータＺ_ｄ,tを要素にもつＺ^をサンプリングにより決定し、音響要素ｄの各々における基底の各状態ｋの基底スペクトルを表わすパラメータＨ_ω,d ^(k)を要素にもつＨ^をサンプリングにより決定し、予め与えられたＬ個の音響イベントｌの各々に対する一部の時刻ｔ_nの各々の音響イベントの有無を表すパラメータＸ_ｌ,tｎ（＝｛０、１｝）を要素にもつＸ_ｌＮ^を用いて、音響要素ｄの各々に対する各時刻ｔの発音の有無を表すパラメータＵ_d,tを要素にもつＵ^をサンプリングにより決定し、音響要素ｄの各々について、パラメータＵ_d,tが０から１へ遷移する遷移確率ａ_（d）をサンプリングにより決定すると共にパラメータＵ_d,tが１から１へ遷移する遷移確率ｂ_（d）をサンプリングにより決定し、予め与えられたＸ_ｌＮ^を用いて、音響イベントｌの各々に対するＤ個の音響要素ｄの各々の重みｗ_ｌ，ｄを要素にもつｗ_ｌ＾を求め、音響イベントｌの各々に対する精度パラメータα_ｌを算出することを、予め定められた収束条件を満たすまで繰り返し行い、Ｃ＾、Ｚ^、Ｈ^、Ｕ^及びｗ_ｌ＾を出力することにより、音響信号の時系列データから、そこに含まれる音響イベントが、どのような音響要素で構成されているのかを精度よく推定することができる、という効果が得られる。

音響信号への音響イベントのマルチラベリングを示すイメージ図である。 音響イベントの音響要素を構成する基底スペクトル集合、状態系列集合、及びアクティベーション集合を示すイメージ図である。 本発明の実施の形態に係る提案モデルのグラフィカル表現図である。 本発明の実施の形態に係る音響信号分析装置の構成を示す概略図である。 本発明の実施の形態に係る音響信号分析装置における音響信号分析処理ルーチンの内容を示すフローチャートである。 本発明の実施の形態の実験結果を示す図である。 本発明の実施の形態における学習データの長さに対する実験結果を示す図である。 本発明の実施の形態と従来技術との性能を比較した実験結果を示す図である。 アクティベーション行列に推定される未知の音響イベントを示す図である。

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

＜発明の概要＞
本発明の実施の形態は、テレビやラジオ放送、ポッドキャスト、動画投稿サイトにおける音響信号の時系列データに対して、そこに含まれる、人に認識されうる音の事象“音響イベント”（話声や歌声、笑い声やあいづちをはじめ、動物の鳴き声、楽器音、環境音、効果音などの音のカテゴリを指す）を書き起こすアルゴリズムに関する。具体的には、音響信号を周波数分析して得られるスペクトログラムをフィルタバンク処理し、その出力値行列に含まれる、様々な音響要素（音声における音素、音楽における楽器音など）の音響的特徴を表現するための基底スペクトル集合、各音響要素の基底スペクトルの状態遷移を表現する状態系列集合、そして各音響要素の発音区間（鳴っているか否か、ＯＮ／ＯＦＦの状態からなる）を表現するアクティベーション集合を抽出し、抽出されたアクティベーションのベイズロジスティック回帰に基づいて、フレームごとに音響イベントをラベル付けする信号解析装置、信号解析方法及び信号解析プログラムに関する。図１は、音響信号から音響イベントを書き起こすための概略図を示す。また、図２は、フィルタバンク出力値を分解して得られる音響要素の基底スペクトル集合、状態系列集合、アクティベーション集合の概略図を示す。

非特許文献６（大石康智、持橋大地、松井知子、中野允裕、亀岡弘和、泉谷知範、柏野邦夫、「ノンパラメトリックベイズアプローチに基づく音響イベント検出」、電子情報通信学会技術研究報告 Vol.112,No.111、2012年6月22日、p.37-42）では、機械学習分野で注目されるノンパラメトリックベイズ法を利用して、大規模な音響信号の時系列データから、そこに含まれる全ての音響要素の音響的特徴とその発音区間を自律的に学習させるフレームワークを提供した（音響イベントはこの音響要素を組み合わせることによって構成されると想定する）。具体的には、音響信号に含まれる音響要素の数（上記図２に示すＤに相当する）が無限個の可能性を持つと仮定し、各音響要素の継続時間を表現するためにマルコフ性を導入したMarkov Indian Buffet Process を用いて、各音響要素の音響的特徴およびそれらのアクティベーションを確率的に推論する。本実施の形態では、推定された音響要素のアクティベーションを利用して、ベイズロジスティック回帰によって、各フレームに音響イベントをラベル付けすることを提供する。

＜原理＞
[１. 音響イベント転写モデル]
音響イベントの重ねあわせをモデル化するために、可変基底型ＮＭＦ（例えば、参考文献１（A. Ozerov, C. F´evotte, and M. Charbit,“ Factorial scaled hidden 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.）を参照）を利用する。これは各音響要素の基底スペクトルが時間にともなって変化するように拡張したＮＭＦである。各音響要素の音色を特徴付けるため、振幅スペクトログラムをメルフィルタバンク処理して得られた出力値行列である音響特徴量行列Ｙ＾＝（Ｙ_ω，ｔ）_Ω×T∈Ｒ^≧0,Ω×Ｔに、このＮＭＦを適用する。なお、記号に付された「＾」は、当該記号が行列または多次元配列またはベクトルである。ここで、ω＝１,...,Ωはメルフィルタバンクの中心周波数のインデックス、ｔ＝１,...,Ｔは時間のインデックス、ｄ＝１,...,Ｄは音響要素のインデックスを表す。このとき、各音響要素の基底スペクトルは時刻t に、ある一つの状態Ｚ_ｄ,ｔ∈Ｎを取ると見なし、

と表現する。これは図２に示すような要素からなる生成モデルである。まず、Ｈ_ｄ＾＝（Ｈ^（ｋ） _ω，ｄ）_Ω×Ｋｄは音響要素ｄのＫ_ｄ個の基底スペクトルを意味する。Ｕ＾＝（Ｕ_ｄ，ｔ）_Ｄ×Tは、音響要素のｏｎ／ｏｆｆを表現する２値系列からなるアクティベーション集合を意味する。このとき、音響信号の音量に関する情報はすべて、基底スペクトルに含まれて表現されると想定する。そして、Ｃ_ｄ＾＝（Ｃ_ω，t，d）_Ω×Ｔは、音響要素ｄのメルフィルタバンク出力値成分である音響特徴量成分に相当し、Ｃ_ω，t，dはＨ^（ｋ） _ω，ｄ，Ｕ_ｄ，ｔ，Ｚ_ｄ,ｔをパラメータとするポアソン分布から生成される。最終的に音響信号はこれらの音響要素の和で表現される。先に述べたように、音響イベントはこれらの音響要素の組み合わせで表現されることを想定している。
さらに、Ｕ＾とＺ＾の事前分布としてｍＩＢＰとＣＲＰを導入し、無限個の音響要素、無限個の基底スペクトルを仮定するモデルへと拡張する。また、半教師あり学習の下で、ベイズロジスティック回帰を利用して、音響要素の組み合わせに音響イベントをラベル付けるアルゴリズムを導入する。以下にその導入方法を説明する。

[１．１ ｍＩＢＰによって生成されるアクティベーション行列]
Markov Indian buffet process(ｍＩＢＰ)と呼ばれるノンパラメトリックベイジアン因子モデル（参考文献３（J. V. Gael et al., “ The infinite factorial hidden Markov model, ”in Proc. NIPS 2008）を参照）は、時間的な依存性を考慮するようにIndian buffet process (IBP)（参考文献４（T. L. Griffiths et al.,“ Infinite latent feature models and the Indian buffet process, ”in Proc. NIPS 2006.）を参照）を拡張したものである。これは各フレームにおける音響要素がｏｎかｏｆｆかどうかをモデル化するためにアクティベーション行列Ｕ＾に対する事前分布を定義する。次の２つの特性を満たす。
１.行数（音響要素の数に相当する）は任意に大きいサイズを想定する
２.各列（フレームに相当する）のｏｎ／ｏｆｆは、音響要素の継続時間（発音時間）を表現するよう、マルコフ過程に従って生起されるこの構成方法が時系列データに対する因子表現を学習することを可能にする。このｍＩＢＰをアクティベーション行列Ｕ＾の事前分布に導入する。具体的には、音響要素の発音時間を表現するために、各音響要素ごとに遷移行列

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

のようなベルヌーイ分布から生成される確率変数として記述できる。このとき、すべての音響要素のアクティベーションからなるアクティベーション行列Ｕ＾は

に従う。ここで、ａ＾＝{ａ_１，．．．，ａ_Ｄ}、ｂ＾＝{ｂ_１，．．．，ｂ_Ｄ}とし、ｃ_ｄ ^００、ｃ_ｄ ^０１、ｃ_ｄ ^１０、ｃ_ｄ ^１１はそれぞれ０→０、０→１、１→０、１→１に遷移する回数とする。さらに、ａ_ｄとｂ_ｄの事前分布をそれらの共役性から、

とする。式（４）における有限モデルをＤ→∞として無限モデルに拡張するために、stick breaking construction（参考文献５（Y. W. Teh et al.,“ Stick-breaking construction for the Indian buffet process, ” in Proc. NIPS 2007）参照）を適用する。ここで、ａ＾をａ₍₁₎＞ａ₍₂₎＞．．．＞ａ_(D)のように順序付けると、Ｄ→∞におけるａ_（ｄ）の生成過程は以下の法則に従う。

一方、ｂ_ｄはベータ分布Ｂｅｔａ（θ_ｂ ^（０），θ_ｂ ^（１））から独立にサンプルされる。したがって、ｄ番目に大きいａ_ｄに対応する変数をｂ_（ｄ）とすると、ｂ_（ｄ）〜Ｂｅｔａ（θ_ｂ ^（０），θ_ｂ ^（１））となる。ｍＩＢＰでは、頻繁にアクティベートされる音響要素ほど、小さいインデックスｄに割り当てられる。

[１．２ ＣＲＰによって生成される時変なスペクトル特性]
音響イベントを表現するために必要な、基底スペクトルの状態数は音響イベントごとに異なる。例えば、ピアノの単音であれば、“attack”、“decay”、“sustain”、“release”と呼ばれる状態があるため、４つの基底スペクトルを用いて表現されるだろう。同様のことが音声を構成する音素にも言える。ドアの開閉音や食器音は、楽器の単音に比べて突発音であるため、基底スペクトルの状態数は少ないかもしれない。交通騒音（ノイズ）は常に1 つの基底スペクトルで表現されるかもしれない。このように、基底スペクトルの状態数は固定するのではなく、音響信号から自動的に決定されることが望ましい。本発明者らは可変基底型ＮＭＦにディリクレ過程を導入した無限状態スペクトルモデルを提案し、音楽音響信号を楽器の単音ごとに分解できることを示した（参考文献２（M. Nakano et al.,“ Infinite-state spectrum model for music signal analysis, ” in Proc. ICASSP 2011）参照）。同様の枠組みを音響イベント検出のために導入する。
音響要素ｄの基底スペクトルの状態系列{Ｚ_ｄ,１，．．．，Ｚ_ｄ,Ｔ}はそれぞれ離散的な値１，．．．，Ｋ_ｄ（すなわち状態のインデックス）をとる。これらの状態系列の同時分布は

を用いて、

となる。

であり、δ（・）はクロネッカーのデルタ関数である。ここで、π＾_ｄの事前分布として、多項分布の共役事前分布である対称なディリクレ分布

を考える。ただし、θ_β ^（ｄ）は正のパラメータとする。このように事前分布を定めると、θ_β ^（ｄ）に対する状態系列の同時分布は

と書ける。このとき、｛Ｚ_ｄ，１,...,Ｚ_{ｄ，ｔ−１}，Ｚ_{ｄ，ｔ＋１},...,Ｚ_ｄ，Ｔ｝ (以降ではＺ＾_ｄ，＼ｔと表現する) が与えられたときのＺ_ｄ，ｔの条件付き確率は次のように与えられる。

ただし、ｎ_ｄ，＼ｔ ^(ｋ)はＺ_ｄ，ｔ’＝ｋ（Ｚ_ｄ，ｔ’∈Ｚ＾_ｄ，＼ｔ）を満たすｔ′の個数を表す。さらに、音響要素を表現するスペクトルの状態数をデータから学習する枠組みを導入する。ＣＲＰ（参考文献６（D. J. Aldous,“ Representations for partially exchangeable arrays of random variables, ”Journal of Multivariate Analysis, vol. 11, pp. 581−598, 1981.）や参考文献７（Y.W.Teh et al., “ Hierarchical Bayesian nonparametric models with applications, ”Cambridge University Press, 2010.）を参照）では、有限として考えた状態数をＫ_ｄ→∞とすることで、Ｚ_ｄ，ｔの条件付き確率を次のように表す。

ここで、Ｋ_ｔ，＼＋はｎ_ｄ，＼ｔ ^(ｋ)＞０となるクラスの個数である。上式から分かるように、無限状態スペクトルモデルにおいて、各時刻に用いられる状態に着目すると、他の時刻に多く用いられている状態ほど使われやすくなる性質がある．また、新しい状態が用いられやすくなるか否かはパラメータθ_β ^（ｄ）に影響する。

[１．３ ベイズロジスティック回帰に基づくマルチラベリング]
アクティベーション行列から、各フレームで発音する音響イベントを推定する（マルチラベリングする）ために、ベイズロジスティック回帰（参考文献８（D. J. C. MacKay,“ The evidence framework applied to classification networks, ”Neural Computation, vol. 4,pp. 448−472, 1992.）を参照）を利用する。これはトピックモデルを教師あり学習できるよう拡張したsLDA（参考文献９（D. M. Blei and J. D. McAuliffe,“ Supervised topic models, ”in Proc. NIPS 2007.）を参照）と類似した適用方法である。まず、フレームt における音響要素のアクティベーションをまとめてＵ_ｔ＾＝［Ｕ_１，ｔ，Ｕ_２，ｔ，．．．，Ｕ_Ｄ，ｔ］^Ｔと表現する。また、音響イベントl に対するラベルをＸ_ｌＮ＾＝｛Ｘ_ｌ，ｔ１，Ｘ_ｌ，ｔ２，．．．，Ｘ_ｌ，ｔＮ｝（Ｘ_ｌ，ｔｎ∈｛０，１｝） と表現すると、その尤度関数は、

と書ける。ここで、ｗ_ｌ＾＝［ｗ_ｌ，１，ｗ_ｌ，２，．．．，ｗ_ｌ，Ｄ］^Ｔは重みベクトルであり、σ（・） はロジスティックシグモイド関数

である。さらに、次の形の単純な等方ガウス事前分布

とα_ｌに対する共役超事前分布をガンマ分布

を導入する。これは超パラメータθ_α ^（０）とθ_α ^（１）で定まる。Ｉ＾_ＤはＤ×Ｄの単位行列を表す。このとき、アクティベーションベクトルＵ_ｔ＾が与えられた下で、音響イベントラベルＸ_ｌ，ｔの予測分布は、

となる。上式では、ロジスティックシグモイド関数とガウス分布のたたみ込み積分を近似するために、プロビット関数に基づく近似予測分布を利用する。ｗ_ｌＭＡＰ＾はｗ＾のＭＡＰ（最大事後確率）解を表し、実際はラベル付けされた少量の学習データを利用して推定される。

音響イベント書き起こしモデルのグラフィカル表現を図3 に示す。これは音響特徴行列Ｙ＾とラベル行列Ｘ＾を共に生成するモデルある。ここで、計算の都合上、基底スペクトル集合Ｈ＾の事前分布として、ガンマ分布Ｈ_ω，ｄ ^（ｋ）〜Ｇａｍｍａ（θ_φ ^（ω，d），θ_ψ ^（ω，d））を仮定する。また、θ_ａ，θ_ｂ ^（０），θ_ｂ ^（１），θ_α ^（０），θ_α ^（１），θ_β ^（ｄ），θ_φ ^（ω，d），θ_ψ ^（ω，d）はすべて固定の超パラメータである。
また、音響特徴量行列Ｙ＾と、ある特定のフレームだけにラベルが付与された音響イベントラベル行列Ｘ＾＝[Ｘ_１Ｎ＾，Ｘ_２Ｎ＾，．．．，Ｘ_ＬＮ＾]^Ｔとの対数尤度関数は、

となる。

[２ パラメータの推論]
ｍＩＢＰのstick breaking constructionは、スライスサンプリングと動的計画法を組み合わせることで、モデルパラメータを推論できる（参考文献３（J. V. Gael et al.,“The infinite factorial hidden Markov model, ”in Proc. NIPS 2008）を参照）。スライスサンプラーが適応的に、無限次元モデルの打ち切りを行い、その後、動的計画法によって各パラメータが推論される。しかしながら、本実施の形態はロジスティック回帰に基づくマルチラベリング処理を含むため（計算コストを軽くするため）、打ち切りされたstick breaking construction の下で（Ｄを大きな値に固定して）、ギブスサンプリングによってパラメータを推論する。以下に各パラメータの推論方法を示す。

[２．１ Ｃ_{ω，ｔ，１}，．．．，Ｃ_{ω，ｔ，Ｄ}の推論]
Ｃ_{ω，ｔ，１}，．．．，Ｃ_{ω，ｔ，Ｄ}の条件付き確率は、

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

とする。この多項分布からＣ_{ω，ｔ，１}，．．．，Ｃ_{ω，ｔ，Ｄ}をサンプリングする。

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

と書ける。ここで、

とする。これらの確率を用いて、Ｚ_ｄ，ｔの状態割り当てを行う。

[２．３ Ｈ_ω，ｄ ^（ｋ）の推論]
Ｈ_ω，ｄ ^（ｋ）の条件付き確率は、

となる。このガンマ分布からＨ_ω，ｄ ^（ｋ）をサンプリングする。

[２．４ Ｕ_ｄ，１，，．．．，Ｕ_ｄ，Ｔの推論]
Forward-filtering backward-sampling アルゴリズム（参考文献１０（S. L. Scott, “Bayesian methods for hidden Markov models: Recursive computing in the 21st century, ”JASA, vol. 97, pp. 337−351, 2002.）を参照）を利用して、Ｕ_ｄ，１，，．．．，Ｕ_ｄ，Ｔを推論する。このとき、Ｕ＾のｄ行以外の値はすべて固定する。まず、ｔ＝１，．．．，Ｔに対して、

を再帰的に計算する。ここで、Ｘ_ｌ，ｔの尤度関数は、ラベルが付与されているフレームｔだけ計算される。次に、

からＵ_ｄ，Ｔをサンプリングする。そして、ｔ＝Ｔ−１，．．．，１に対して、Ｕ_{ｄ，ｔ＋１}が与えられた下で、

に従って、Ｕ_ｄ，ｔを後方から順番にサンプリングすれば、Ｕ_ｄ，１，，．．．，Ｕ_ｄ，Ｔが求まる。

[２．５ ａ_（ｄ）の推論]
ｄ＝１，．．．，Ｄ に対して、ａ_（ｄ） の条件付き確率は、

となり、このベータ分布からａ_（ｄ）をサンプリングする。

[２．６ ｂ_（ｄ）の推論]
ｄ＝１，．．．，Ｄに対して、ｂ_（ｄ）の条件付き確率は、

となり、このベータ分布からｂ_（ｄ）をサンプリングする。

[２．７ ｗ_ｌ＾の推論]
音響信号の一部のＮフレームには音響イベントのｏｎ／ｏｆｆに関するラベルが付いているとする。ラベルが付けられたフレームのインデックスをｔ_１，ｔ_２，．．．，ｔ_Ｎとし、これらのフレームの音響要素のアクティベーション｛Ｕ_ｔ１＾，Ｕ_ｔ２＾，．．．，Ｕ_ｔＮ＾｝ と音響イベントｌに対するラベルＸ_ｌＮ＾＝［Ｘ_ｌ，ｔ１，Ｘ_ｌ，ｔ２，．．．，Ｘ_ｌ，ｔＮ］^Ｔを用いて、ｗ_ｌを推論する。ロジスティックシグモイド関数の積である尤度関数と事前確率分布が共役な関係ではないため、ここでは、ニュートン-ラフソン法を利用してＭＡＰ（最大事後確率）推定値ｗ_ｌＭＡＰ＾を求める。更新式は、

と導出される。ここで、

とする。ｗ_ｌ ^{（ｎｅｗ）}＾とｗ_ｌ ^{（ｏｌｄ）}＾の差が10⁻⁸ 以下になったとき、ｗ_ｌＭＡＰ＾＝ｗ_ｌ ^{（ｎｅｗ）}＾とする。

[２．８ α_ｌの推論]
ｌ＝１，．．．，Ｌに対して、α_ｌの条件付き確率は、

と近似され、このガンマ分布からα_ｌをサンプリングする。

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

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

本実施の形態に係る音響信号分析装置は、入力部１と、記憶部２と、演算部３と、出力部４とを備えている。本実施の形態は、分析対象である音響信号の時系列の一部分に対してＬ個の音響イベントの値をラベル付けしたデータを学習データとして入力する。そして、前述したパラメータ推定アルゴリズムを用いて信号解析を行い、分析対象である音響信号の時系列の全ての部分に対して、推定したパラメータを適用する構成である。

入力部１には、分析対象である音響信号の時系列と、当該音響信号の時系列について、ある特定のフレームだけにラベルが付与された音響イベントラベル行列Ｘ＾＝[Ｘ_１Ｎ＾，Ｘ_２Ｎ＾，．．．，Ｘ_ＬＮ＾]^Ｔとが入力される。ここで、Ｘ＾は、音響信号の時系列の一部分（ある特定のフレーム）について予め付与された教師ラベル行列である。

記憶部２は、入力部１に入力された音響信号の時系列と、Ｘ＾とを記憶する。また、記憶部２には、予め定められた超パラメータθ_ａ＝１，θ_ｂ ^（０）＝１０，θ_ｂ ^（１）＝１，θ_β ^（ｄ）＝１，θ_φ ^（ω，d）＝１，θ_ψ ^（ω，d）＝１ （ｄ ＝１，．．．，Ｄ，ω＝１，．．．，Ω），θ_α ^（０）＝１，θ_α ^（１）＝１０００が記憶されている。これらの超パラメータの各値は予備実験によって決定した。また、記憶部２は、後述する各処理での結果等を記憶する。

演算部３は、音響特徴抽出部３０と、パラメータ初期値生成部３２と、モデルパラメータ更新部３４と、音響イベントラベル推論部３６とから構成されている。

音響特徴抽出部３０は、記憶部２に記憶された音響信号の時系列を読み出し、当該音響信号に基づいて、音響特徴量行列Ｙ＾（以下、Ｙ＾と称する。Ｙ＾は二次元配列である。）を作成する。本実施の形態では、音響イベントの音色に着目し、フレーム長１００ｍｓ、フレームシフト長１００ｍｓとして、フレームごとに短時間フーリエ変換を行い、その振幅スペクトルをメルフィルタバンク処理して得られる出力値であるＹ_ω,tを各要素にもつ二次元配列Ｙ＾（ω＝１,...,Ω、ｔ＝１,...,Ｔ）を生成し、記憶部２へ格納する。このとき、例えば、メルフィルタの数を２４個とすれば、Ω＝２４であり、Ｔはフレームの総数に相当する。また、上記（１）式に示すように、本実施の形態では、ポアソン分布に従う確率変数として、音響特徴量行列をモデル化するため、Ｙ＾の全ての要素の値を整数値に丸め込む。

パラメータ初期値生成部３２は、基底スペクトル集合Ｈ＾（以下、Ｈ＾と称する。）、アクティベーション集合Ｕ＾（以下、Ｕ＾と称する。）、状態系列Ｚ＾（以下、Ｚ＾と称する。）、重みベクトルｗ_ｌ＾（以下、ｗ_ｌ＾と称する。）、精度パラメータα_ｌ、及び遷移確率ａ＾並びにｂ＾の各パラメータの初期値を生成し、記憶部２へ格納する。
なお、Ｈ＾は三次元配列、Ｕ＾とＺ＾とは二次元配列、ｗ_ｌ＾とａ＾とｂ＾とは一次元配列である。

パラメータ初期値生成部３２は、音響特徴抽出部３０により作成されたＹ＾に基づいて、上記（１）式に従って、通常のＮＭＦを適用し、推定されたＨ＾とＵ＾を初期値とし、記憶部２に格納する。ここで、推定されたＨ_{１：Ω，ｄ}はｄ番目の音響要素の初期基底スペクトルとする（初期状態数はＫ_ｄ＝１とする）。

具体的には、パラメータ初期値生成部３２は、Ｈ＾について、各音響要素ｄの基底の状態ｋを１のみとする。そして、Ｈ＾の要素である全てのパラメータＨ_ω,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.）を参照）ため、詳細な説明を省略する。

また、Ｕ＾については、その中央値よりも大きい要素は１に、中央値よりも小さい要素は０に二値化して初期値とする。

Ｚ＾については、全ての要素の初期値として１が設定され、記憶部２に格納される。

音響要素数Ｄについては、適切な初期値が設定され、記憶部２に格納される。本実施の形態では、打ち切りstick breaking construcitonをＵ＾に適用するため、音響要素数Ｄの値はできる限り大きい値が望ましい。後述する評価実験では、Ｄ＝１００とした。

遷移確率ａ＾については、上記（６）式に従って、ａ＾＝{ａ_１,...,ａ_Ｄ}についての初期値が生成され、そして遷移確率ｂ＾については、上記（５）式に従って、初期値が生成され、各々の初期値が記憶部２に格納される。

ｗ_ｌ＾、精度パラメータα_ｌについては、音響イベントｌ＝１，．．．，Ｌに対して、初期値として適切な値が設定され、各々の初期値が記憶部２に格納される。

モデルパラメータ更新部３４は、前述したパラメータ推定アルゴリズムを用いて、各パラメータを更新して各パラメータを推定する。モデルパラメータ更新部３４は、音源分離推論部３４０と、状態系列推論部３４１と、基底スペクトル推論部３４２と、アクティベーション推論部３４３と、遷移確率推論部３４４と、重みベクトル推論部３４５と、精度パラメータ推論部３４６と、収束判定部３４７とから構成されている。

音源分離推論部３４０は、上記（１８）式及び（１９）式に従って、数値的サンプリングによって、各音響要素ｄについての音響特徴量成分を表す行列Ｃ＾（以下、Ｃ＾と称する。Ｃ＾は三次元配列である。）のすべての要素を推論する。具体的には、記憶部２に記憶されたＹ＾、Ｚ＾、Ｈ＾、及びＵ＾に基づいて、上記（１８）式及び（１９）式に従って、Ｃ＾のすべての要素をサンプリングして推論し、記憶部２に格納する。

状態系列推論部３４１は、上記（２０）式及び（ｃ）式から、ギブスサンプリングによってＺ＾のすべての要素の推論を行う。具体的には、記憶部２に記憶されたＹ＾、超パラメータθ_φ ^（ω，d），θ_ψ ^（ω，d）（ｄ ＝１，．．．，Ｄ，ω＝１，．．．，Ω）、Ｃ＾、Ｚ＾、及びＵ＾に基づいて、上記（２０）式及び（ｃ）式に従って、ｐ（Ｚ_ｄ,１,...,Ｚ_ｄ,Ｔ｜Ｃ_・,１,ｄ,...,Ｃ_・,Ｔ,ｄ）の値の大きなＺ_ｄ,１,...,Ｚ_ｄ,Ｔを確率的に求め、記憶部２に格納する。なお、Ｚ_ｄ,ｔを求める際のＺ＾については、Ｚ_ｄ,ｔ以外の｛Ｚ_ｄ,１,...,Ｚ_{ｄ,ｔ−１},Ｚ_{ｄ,ｔ＋１},...,Ｚ_ｄ,ｔ｝（＝Ｚ_ｄ,＼ｔ）の値を既知として、Ｚ_ｄ,ｔをサンプリングして推論し、記憶部２に格納する。

基底スペクトル推論部３４２は、上記（２１）式から、数値的サンプリングによってＨ＾のすべての要素を推論する。具体的には、記憶部２に記憶された超パラメータθ_φ ^（ω，d），θ_ψ ^（ω，d）（ｄ ＝１，．．．，Ｄ，ω＝１，．．．，Ω）、Ｃ＾、Ｚ＾、及びＵ＾に基づいて、上記（２１）式に従って、数値的サンプリングによってＨ＾のすべての要素をサンプリングして推論し、記憶部２に格納する。

アクティベーション推論部３４３は、上記（２２）式及び（２３）式からなるForward filtering backward sampling アルゴリズムに基づいて、Ｕ＾のすべての要素を推論する。具体的には、記憶部２に記憶されたＹ＾、Ｃ＾、Ｚ＾、Ｈ＾、Ｘ＾及びＵ＾に基づいて、まずForward filteringによって、上記（２２）式に従って、各ｔに対するＵ_ｄ,ｔの事後分布を再帰的に計算した後、Ｂackward samplingによって、上記（２３）式に従って、各ｔに対するＵ_ｄ,ｔを数値的にサンプリングして求め、記憶部２に格納する。なお、Forward filteringでは、Ｕ_ｄ,ｔの事後分布を求める際に、Ｕ_ｄ,ｔ-1の事後分布を用い、Ｂackward samplingでは、Ｕ_ｄ,ｔを求める際に、Ｕ_ｄ,ｔ＋1を用いる。

遷移確率推論部３４４は、記憶部２に記憶された超パラメータθ_ｂ ^（０）、θ_ｂ ^（１）及びＵ＾に基づいて、上記（２４）式及び（２５）式に従って、遷移確率ａ_（ｄ）及びｂ_（ｄ）を推論する。具体的には、記憶部２に記憶された遷移確率ａ＾、及びＵ＾に基づいて、上記（２４）式に従って、ｄ＝１,...,Ｄに対するａ_（ｄ）をサンプリングして推論し、記憶部２に格納する。また、ｂ_（ｄ）については、記憶部２に記憶された超パラメータθ_ｂ ^（０）、θ_ｂ ^（１）、及びＵ＾に基づいて、上記（２５）式に従って、ｄ＝１,...,Ｄに対するｂ_（ｄ）をサンプリングにより推論し、記憶部２に格納する。

重みベクトル推論部３４５は、音響イベントｌ＝１，．．．，Ｌに対して、ニュートン-ラフソン法を利用して、最大事後確率推定値ｗ_ｌＭＡＰ＾を求める。具体的には、音響イベントｌ＝１，．．．，Ｌに対して、記憶部２に記憶されたＵ＾、Ｘ_ｌＮ＾、精度パラメータα_ｌ、及びｗ_ｌ ^{（ｏｌｄ）}＾に基づいて、上記（２６）式に従って、重みベクトルｗ_ｌ ^{（ｎｅｗ）}＾を更新することを繰り返し、更新したｗ_ｌ ^{（ｎｅｗ）}＾とｗ_ｌ ^{（ｏｌｄ）}＾の差が１０^ー８以下になったとき、ｗ_ｌＭＡＰ＾＝ｗ_ｌ ^{（ｎｅｗ）}＾として最大事後確率推定値ｗ_ｌＭＡＰ＾を更新し、記憶部２に格納する。

精度パラメータ推論部３４６は、音響イベントｌ＝１，．．．，Ｌに対して、記憶部２に記憶された超パラメータθ_α ^（０）、θ_α、及びｗ_ｌＭＡＰ＾に基づいて、上記（２８）に従って、数値的サンプリングによって、α_ｌを推論し、記憶部２に格納する。

収束判定部３４７は、記憶部２に記憶されたＹ＾及びＸ＾と、推論されたＺ＾、Ｈ＾、Ｕ＾、ｗ＾とを用いて、上記（１７）式の対数尤度関数を計算する。更新前のモデルパラメータを用いて計算した対数尤度関数の値と更新後の値との誤差が、所定の閾値以下であれば、モデルパラメータの更新が収束したと判定する。一方、誤差が所定の閾値を超えるようであれば、モデルパラメータの更新は収束していないと判定し、音源分離推論部３４０の処理に戻る。本実施の形態ではこの誤差の閾値を実験的にε＝１．０×１０^−５とした。

音響イベントラベル推論部３６は、Ｌ個の音響イベントｌの各々について、記憶部２に記憶されたＵ＾、Ｘ_ｌＮ^、精度パラメータα_ｌ及びｗ_ｌＭＡＰ＾に基づいて、上記（１６）式に従って、各時刻ｔのパラメータＸ_ｌ,tの予測確率を算出し、時刻ｔ毎に、時刻ｔの予測確率が、０．５より大きい場合には、パラメータＸ_ｌ,tの値を１とし、時刻ｔの予測確率が、０．５以下である場合には、パラメータＸ_ｌ,tの値を０とする。

出力部４は、記憶部２に記憶されているモデルパラメータＣ＾、Ｚ＾、Ｈ＾、Ｕ＾、ａ＾、ｂ＾、ｗ_ｌＭＡＰ＾をすべて出力すると共に、音響イベントラベル推論部３６により求められた各音響イベントｌに対する時刻ｔ毎のパラメータＸ_ｌ,tを出力する。

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

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

まず、ステップＳ１００において、記憶部２から、各フレーム内の音響信号の時系列を読み込む。そして、音響信号の時系列に対して、短時間フーリエ変換を用いた時間周波数分析を行い、振幅スペクトルを結果として得る。そして、その振幅スペクトルに対し、メルフィルタバンク処理して得られる音響特徴量Ｙ_ω,tを各（ω，ｔ）の要素にもつ二次元配列Ｙ＾を生成し、記憶部２へ格納する。

そして、ステップＳ１０２において、記憶部２に記憶されている超パラメータθ_ａ、θ_ｂ ^（０）、θ_ｂ ^（１）、θ_β ^（ｄ）、θ_φ ^（ω，d）、θ_ψ ^（ω，d）（ｄ ＝１，．．．，Ｄ，ω＝１，．．．，Ω）、θ_α ^（０）、θ_α ^（１）の値を読み込む。

次にステップＳ１０６において、パラメータ初期値生成部３２によって、Ｈ＾、Ｕ＾、Ｚ＾、音響要素数Ｄ、遷移確率ａ＾、ｂ＾、ｗ_ｌ＾、及び精度パラメータα_ｌの各パラメータの初期値を生成する。
具体的には、上記（１）式に従って、上記ステップＳ１００において生成されたＹ＾に基づいて、Ｈ＾とＵ＾を推定し、初期値として記憶部２へ記憶する。Ｚ＾については、全ての要素の初期値として１を設定し、音響要素数Ｄについては、初期値として例えばＤ＝１００を設定する。また遷移確率ａ＾については、上記（６）式に従って、ａ＾＝{ａ_１,...,ａ_Ｄ}についての初期値が生成され、遷移確率ｂ＾については、上記（５）式に従って初期値が生成される。ｗ_ｌ＾、精度パラメータα_ｌについては、音響イベントｌ＝１，．．．，Ｌに対して、初期値として適切な値が設定される。そして、生成された各々の初期値を記憶部２へ記憶する。

ステップＳ１０８では、音源分離推論部３４０によって、Ｃ＾の全ての要素を推論する。具体的には、上記ステップＳ１００において生成されたＹ＾、並びに上記ステップＳ１０６で設定されたＺ＾、Ｈ＾及びＵ＾、又は、後述するステップＳ１１０〜ステップＳ１１４で前回推論されたＺ＾、Ｈ＾及びＵ＾に基づいて、上記（１８）式及び（１９）式に従って、Ｃ＾のすべての要素をサンプリングして推論し、記憶部２へ記憶する。

ステップＳ１１０では、状態系列推論部３４１によって、ギブスサンプリングによって、Ｚ＾のすべての要素の推論を行う。具体的には、上記ステップＳ１００で生成されたＹ＾、上記ステップＳ１０２で読み込まれた超パラメータθ_φ ^（ω，d），θ_ψ ^（ω，d）（ｄ ＝１，．．．，Ｄ，ω＝１，．．．，Ω）、上記ステップＳ１０６で設定されたＺ＾及びＵ＾、又は本ステップＳ１１０及び後述するステップＳ１１４で前回推論されたＺ＾及びＵ＾、並びに上記ステップＳ１０８で推論されたＣ＾に基づいて、上記（２０）式及び（ｃ）式に従って、Ｚ_ｄ,１,...,Ｚ_ｄ,Ｔをサンプリングして推論する。
なお、Ｚ_ｄ,ｔを求める際の状態系列Ｚ＾については、Ｚ_ｄ,ｔ以外の｛Ｚ_ｄ,１,...,Ｚ_{ｄ,ｔ−１},Ｚ_{ｄ,ｔ＋１},...,Ｚ_ｄ,ｔ｝（＝Ｚ_ｄ,＼ｔ）の値を既知として、Ｚ_ｄ,ｔをサンプリングして推論する。ここで、Ｚ_ｄ,１,...,Ｚ_{ｄ,ｔ−１}までは、本ステップＳ１１０で推論された直前の値を用い、Ｚ_{ｄ,ｔ＋１},...,Ｚ_ｄ,ｔについては、上記ステップＳ１０６で生成されたＺ＾又は本ステップＳ１１０で前回推論されたＺ＾を用いて、Ｚ_ｄ,ｔをサンプリングして推論し、記憶部２へ記憶する。

ステップＳ１１２では、基底スペクトル推論部３４２によって、Ｈ＾の全ての要素の推論を行う。具体的には、上記ステップＳ１０２で読み込まれた超パラメータθ_φ ^（ω，d），θ_ψ ^（ω，d）（ｄ ＝１，．．．，Ｄ，ω＝１，．．．，Ω）、上記ステップＳ１０６で設定されたＵ＾、又は後述するステップＳ１１４で前回推論されたＵ＾、上記ステップＳ１０８で推論されたＣ＾、及び上記ステップＳ１１０で推論されたＺ＾に基づいて、上記（２１）式に従って、Ｈ＾のすべての要素をサンプリングして推論し、記憶部２へ記憶する。

ステップＳ１１４では、アクティベーション推論部３４３によって、Ｕ＾の全ての要素の推論を行う。具体的には、上記ステップＳ１００において生成されたＹ＾、記憶部２に記憶されたＸ_ｌＮ＾、上記ステップＳ１０８で推論されたＣ＾、上記ステップＳ１１０で推論されたＺ＾、上記ステップＳ１１２で推論されたＨ＾、及び上記ステップＳ１０６で設定されたＵ＾又は本ステップＳ１１４で前回推論されたＵ＾に基づいて、上記（２２）式及び（２３）式に従って、Ｕ＾の全ての要素の推論を行う。詳細には、まずForward filteringによって、上記（２２）式に従って、各ｔに対するＵ_ｄ,ｔの事後分布を再帰的に計算した後、Ｂackward samplingによって、上記（２３）式に従ってＵ_ｄ,ｔを数値的にサンプリングして求めることにより、Ｕ＾の全ての要素を推論し、記憶部２へ記憶する。

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

ステップＳ１１８において、重みベクトル推論部３４５によって、ｗ_ｌＭＡＰ＾の推論を行う。具体的には、音響イベントｌ＝１，．．．，Ｌに対して、上記ステップＳ１１４で推論されたＵ＾、記憶部２に記憶されたＸ_ｌＮ＾、上記ステップＳ１０６で設定されたか、又は後述するステップＳ１２０で前回推論された精度パラメータα_ｌ、及び上記ステップＳ１０６で設定されたか、又は本ステップＳ１１８で前回推論されたｗ_ｌ ^{（ｏｌｄ）}＾に基づいて、上記（２６）式に従って、重みベクトルｗ_ｌ ^{（ｎｅｗ）}＾を更新することを繰り返し、更新したｗ_ｌ ^{（ｎｅｗ）}＾とｗ_ｌ ^{（ｏｌｄ）}＾の差が１０^ー８以下になったとき、ｗ_ｌＭＡＰ＾＝ｗ_ｌ ^{（ｎｅｗ）}＾として最大事後確率推定値ｗ_ｌＭＡＰ＾を更新し、記憶部２に格納する。

ステップＳ１２０において、精度パラメータ推論部３４６によって、α_ｌを推論する。
具体的には、音響イベントｌ＝１，．．．，Ｌに対して、上記ステップＳ１０２で読み込まれた超パラメータθ_α ^（０）、θ_α ^（１）、及び上記ステップＳ１１８で推論されたｗ_ｌＭＡＰ＾に基づいて、上記（２８）に従って、数値的サンプリングによって、α_ｌを推論する。そして、推論されたα_ｌを記憶部２へ記憶する。

ステップＳ１２２では、収束判定部３４７によって、上記ステップＳ１００で得られたＹ＾及び記憶部２に記憶されたＸ＾と、上記ステップＳ１１０〜Ｓ１１８で推論されたＺ＾、Ｈ＾、Ｕ＾、及びｗ_ｌ＾とを用いて、上記（１７）式の対数尤度関数を計算する。そして、更新前のモデルパラメータを用いて計算した対数尤度関数の値と更新後の値との誤差が、所定の閾値以下であれば、収束したと判定する。収束していないと判定された場合には、ステップＳ１０８へ移行し、上記ステップＳ１０８〜ステップＳ１２０で推論したパラメータを用いて、上記ステップＳ１０８〜ステップＳ１２０の処理を繰り返す。収束したと判定された場合には、ステップＳ１２４へ移行する。

ステップＳ１２４において、音響イベントラベル推論部３６によって、Ｌ個の音響イベントｌの各々について、記憶部２に記憶されたＵ＾、Ｘ_ｌＮ^、精度パラメータα_ｌ及びｗ_ｌＭＡＰ＾に基づいて、上記（１６）式に従って、各時刻ｔのパラメータＸ_ｌ,tの予測確率を算出する。そして、Ｌ個の音響イベントｌの各々について、時刻ｔ毎に、時刻ｔの予測確率が、０．５より大きい場合には、パラメータＸ_ｌ,tの値を１とし、時刻ｔの予測確率が、０．５以下である場合には、パラメータＸ_ｌ,tの値を０とする。

ステップＳ１２６では、出力部４によって、推論された各パラメータ（Ｃ＾、Ｚ＾、Ｈ＾、Ｕ＾、ａ＾、ｂ＾、ｗ_ｌＭＡＰ＾）を結果として出力すると共に、音響イベントラベル推論部３６により求められた各音響イベントｌの時刻ｔ毎のパラメータＸ_ｌ,tを出力し、音響信号分析処理ルーチンが終了する。

＜実験結果＞
英語学習用ポッドキャストの音響信号（計３５分のうち、はじめの５分）を用いて、本技術の性能を評価する。音響信号は１６ｋＨｚサンプリングで１６ビット量子化されたものである。

図６の上段の図は、１００ｍｓごとに手動でラベル付けされた音響イベントラベルを示す。黒色が音響イベントのｏｎを、白色がｏｆｆを表す。この音響信号は“音楽”、“効果音”、“電話のベル音”、５名の“音声”からなる８種類の音響イベントを含む。音響信号は、フレームシフト長１００ｍｓ、フレーム長１００ｍｓ、ハニング窓を用いてフレームに分割され、短時間フーリエ変換によって振幅スペクトログラムに変換される。そして、各フレームの振幅スペクトルをフィルタバンク処理して得られる出力値（２４個の値）を、音響特徴量行列とする（すなわち、Ｙ＾は２４×３０００の行列である）。

図６の上段の図に示されるように、学習データ（あらかじめラベル付けされたデータ）として３ 種類（（１）５０秒、（２）１００秒、（３）１５０秒）を用意した。そして、後半の１５０秒に対する推定結果に基づいて、本技術の性能を評価する。評価尺度として、各イベントのラベル予測確率｛ｐ（Ｘ_{｛ｌ，１５０１｝}｜Ｕ＾_ｔ，Ｘ＾_ｌＮ）， ｐ（Ｘ_{｛ｌ，１５０２｝}｜Ｕ＾_ｔ，Ｘ＾_ｌＮ），…ｐ（Ｘ_{｛ｌ，３０００｝}｜Ｕ＾_ｔ,Ｘ＾_ｌＮ）｝と正解ラベルとを比較することで計算されるＲＯＣ曲線の下側面積ＡＵＣの値を利用する。

図６の中段及び下段の図は、アクティベーション行列Ｕ＾と学習データ（３）を利用して推定されたラベル付け結果を示す。
Ｘ_ｌ，ｔ＝Ｉ（ｐ（Ｘ_{｛ｌ,ｔ｝}|Ｕ＾_ｔ，Ｘ＾_ｌＮ)＞０．５）
の基準に基づいて、各フレームに付与する音響イベントｌのラベルの値が０であるか１であるかを判定した。ここで、ｐ（Ｘ_{｛ｌ,ｔ｝}｜Ｕ＾_ｔ，Ｘ＾_ｌＮ)は予測確率である。また、Ｉ（Ａ）はインジケータ関数であり、条件Ａが真であれば１を、偽であれば０を返す関数である。アクティベーション行列のロジスティック回帰によって推定されるラベル付結果が、正解ラベルと全体的に近い結果が得られたことがわかる。また、女性の話者を識別することが難しいこともわかる。

図７は、学習データの長さに対するラベル付性能を評価した図である。ここで、“電話のベル音” と“男性Ａの音声” は評価データに含まれていないため、これらについては評価しない。学習データの長さが増えるにつれて、各イベントのラベル付性能を表すＡＵＣの値が向上することがわかる。

図８は、ベースライン手法と本実施の形態の性能を比較する図である。どちらも学習データ（３）を使用した結果である。ベースライン手法では、音響イベントのラベルが付与されるフレームの音響特徴量（２４個のメルフィルタバンク出力値）の分布をガウス混合モデル（ＧＭＭ）を利用して学習し、評価データに対して、各音響イベントラベルの事後確率を計算する。そして、５点移動平均フィルタを利用して、事後確率を平滑化し、平均ＡＵＣの値を計算した。本技術はベースライン手法と同等以上の性能が得られていることがわかった。ただし注目すべきは、（１）や（２）のように学習データが少ない場合、ベースライン手法はＧＭＭを効果的に学習できない。しかしながら、本技術は、たとえ少量の学習データであったとしても、ラベル付けされていないデータをうまく活用することによって、全体のイベントラベルを推定できることが分かった。

また、識別モデルにとって、学習データに含まれていない未知の音響イベントを検出することは難しい。しかし、生成的なアプローチである本技術が未知の音響イベントを検出できる可能性があることを示す。ここでは、“効果音”に着目する。実際、この効果音は様々な種類の音に対してラベル付けされているため、性能は低い（上記図７参照）。しかし、推定されたアクティベーション行列を見てみると興味深い（図９参照）。“効果音” のラベルが付与される区間において、特定の音響要素がアクティブとなっていることがわかる。これらの音響要素が新しい音響イベントとして検出されていることを意味する。

以上説明したように、本発明の実施の形態に係る音響信号分析装置によれば、音響要素ｄの各々に対する音響特徴量成分Ｃ_ω,t,ｄを要素にもつＣ＾をサンプリングにより決定し、音響要素ｄの各々に対して各時刻ｔに基底の状態が何れであるか、または新規の状態であることを表わすパラメータＺ_ｄ,tを要素にもつＺ^をサンプリングにより決定し、音響要素ｄの各々における基底の各状態ｋの基底スペクトルを表わすパラメータＨ_ω,d ^(k)を要素にもつＨ^をサンプリングにより決定し、予め与えられたＬ個の音響イベントｌの各々に対する一部の時刻ｔ_nの各々の音響イベントの有無を表すパラメータＸ_ｌ,tｎ（＝｛０、１｝）を要素にもつＸ_ｌＮ^を用いて、音響要素ｄの各々に対する各時刻ｔの発音の有無を表すパラメータＵ_d,tを要素にもつＵ^をサンプリングにより決定し、音響要素ｄの各々について、パラメータＵ_d,tが０から１へ遷移する遷移確率ａ_（d）をサンプリングにより決定すると共にパラメータＵ_d,tが１から１へ遷移する遷移確率ｂ_（d）をサンプリングにより決定し、予め与えられたＸ_ｌＮ^を用いて、音響イベントｌの各々に対するＤ個の音響要素ｄの各々の重みｗ_ｌ，ｄを要素にもつｗ_ｌ＾を求め、音響イベントｌの各々に対する精度パラメータα_ｌを算出することを、予め定められた収束条件を満たすまで繰り返し行い、Ｃ＾、Ｚ^、Ｈ^、Ｕ^及びｗ_ｌ＾を出力することにより、音響信号の時系列データから、そこに含まれる音響イベントが、どのような音響要素で構成されているのかを精度よく推定することができる。

また、本発明の実施の形態に係る音響信号分析装置によれば、音響信号の時系列データ及び一部のフレームに対する音響イベントのラベル行列から、全ての音響イベントの音響的特徴とその発音区間を自律的に推定することができる。

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

例えば、本実施の形態における音響特徴抽出部３０では、音響信号の時系列に対し、フレームごとに短時間フーリエ変換を行い、その振幅スペクトルをメルフィルタバンク処理してＹ＾を生成する場合を例に説明したが、これに限定されるものではなく、振幅スペクトログラムそのものを音響特徴量行列Ｙ＾としてもよい。

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

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

１ 入力部
２ 記憶部
３ 演算部
４ 出力部
３０ 音響特徴抽出部
３２ パラメータ初期値生成部
３４ モデルパラメータ更新部
３６ 音響イベントラベル推論部
３４０ 音源分離推論部
３４１ 状態系列推論部
３４２ 基底スペクトル推論部
３４３ アクティベーション推論部
３４４ 遷移確率推論部
３４５ 重みベクトル推論部
３４６ 精度パラメータ推論部
３４７ 収束判定部

Claims (5)

1. 音響特徴量抽出手段、パラメータ初期値設定手段、音源分離手段、状態系列推論手段、基底スペクトル推論手段、アクティベーション推論手段、遷移確率推論手段、重みベクトル推論手段、精度パラメータ推論手段、及び収束判定手段を含む音響信号分析装置における音響信号分析方法であって、
   前記音響特徴量抽出手段によって、複数の音響要素が混在する音響信号の時系列データを入力として、音響特徴量Ｙ_ω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙ＾を出力するステップと、
   前記パラメータ初期値設定手段によって、予め定められたＤ個の音響要素ｄの各々における所定の基底の状態ｋの基底スペクトルを表わすパラメータＨ_ω,d ^(k)を要素にもつ三次元配列Ｈ^、前記Ｄ個の音響要素の各々に対する各時刻ｔの発音の有無を表すパラメータＵ_d,t（＝｛０、１｝）を要素にもつ二次元配列Ｕ^、及び前記Ｄ個の音響要素ｄの各々に対して各時刻ｔに前記基底の状態が何れであるかを表わすパラメータＺ_ｄ,tを要素にもつ二次元配列Ｚ^、予め定められたＬ個の音響イベントｌの各々に対する前記Ｄ個の音響要素ｄの各々の重みｗ_ｌ，ｄを要素にもつ一次元配列ｗ_ｌ＾、前記Ｌ個の音響イベントｌの各々に対する精度パラメータα_ｌの各々の初期値を設定すると共に、前記Ｄ個の音響要素ｄの各々について、前記パラメータＵ_d,tが０から１へ遷移する遷移確率ａ_d（ａ₁＞ａ₂＞・・・＞ａ_D-1＞ａ_D）、及び前記パラメータＵ_d,tが１から１へ遷移する遷移確率ｂ_dの各々の初期値を設定するステップと、
   前記音源分離手段によって、（ω、ｔ、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^に基づいてパラメータｐ_ω,t,dを算出し、算出した各パラメータｐ_ω,t,dと、前記二次元配列Ｙ^とをパラメータとする多項分布に従って、Ｄ個の音響要素ｄの各々に対する音響特徴量成分Ｃ_ω,t,ｄを要素にもつ三次元配列Ｃ＾をサンプリングにより決定するステップと、
   前記状態系列推論手段によって、前記Ｄ個の音響要素ｄの各々について、前記二次元配列Ｚ^、前記三次元配列Ｃ^，及び前記二次元配列Ｕ^に基づいて、各時刻ｔに前記基底の状態が各状態ｋである事後確率、及び各時刻ｔに前記基底の状態が新規の状態である事後確率を算出し、算出した事後確率に従って、前記二次元配列Ｚ^をサンプリングにより決定するステップと、
   前記基底スペクトル推論手段によって、（ｋ、ω、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^及び前記三次元配列Ｃ^に基づいて、パラメータφ_ω,d ^(k)、ψ_ω,d ^(k)を算出し、算出したパラメータφ_ω,d ^(k)、ψ_ω,d ^(k)をパラメータとする確率分布に従って、前記三次元配列Ｈ^をサンプリングにより決定するステップと、
   前記アクティベーション推論手段によって、（ｄ、ｔ）の全ての組み合わせについて、予め与えられた前記Ｌ個の音響イベントｌの各々に対する一部の時刻ｔnの各々の音響イベントの有無を表すパラメータＸ_ｌ,tｎ（＝｛０、１｝）を要素にもつ一次元配列Ｘ_ｌＮ^、前記三次元配列Ｃ^、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記遷移行列ａ_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をサンプリングにより決定するステップと、
   前記重みベクトル推論手段によって、前記Ｌ個の音響イベントｌの各々について、前記一次元配列Ｘ_ｌＮ^、前記二次元配列Ｕ^、前記一次元配列ｗ_ｌ＾、及び前記パラメータα_ｌに基づいて、事後確率が最大となるように前記一次元配列ｗ_ｌ＾を求めるステップと、
   前記精度パラメータ推論手段によって、前記Ｌ個の音響イベントｌの各々について、前記一次元配列ｗ_ｌ＾に基づく前記パラメータα_lの確率分布に従って、前記パラメータα_ｌをサンプリングにより決定するステップと、
   前記収束判定手段によって、予め定められた収束条件を満たすまで、前記音源分離手段、前記状態系列推論手段、前記基底スペクトル推論手段、前記アクティベーション推論手段、前記遷移確率推論手段、前記重みベクトル推論手段、及び前記精度パラメータ推論手段による各処理を繰り返し行い、前記三次元配列Ｃ＾、前記二次元配列Ｚ^、前記三次元配列Ｈ^、前記二次元配列Ｕ^、及び前記Ｌ個の音響イベントｌの各々の前記一次元配列ｗ_ｌ＾を出力するステップと、
   を含む音響信号分析方法。
2. 音響イベントラベル推論手段によって、前記Ｌ個の音響イベントｌの各々について、前記二次元配列Ｕ^、前記一次元配列Ｘ_ｌＮ^、前記精度パラメータα_ｌ、及び前記一次元配列ｗ_ｌ＾に基づいて、各時刻ｔのパラメータＸ_ｌ,tの予測確率を算出し、前記予測確率に基づいて、各時刻ｔのパラメータＸ_ｌ,tの値を求めるステップを更に含む請求項１記載の音響信号分析方法。
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,tが１から１へ遷移する遷移確率ｂ_dの各々の初期値を設定するパラメータ初期値設定手段と、
   （ω、ｔ、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^に基づいてパラメータｐ_ω,t,dを算出し、算出した各パラメータｐ_ω,t,dと、前記二次元配列Ｙ^とをパラメータとする多項分布に従って、Ｄ個の音響要素ｄの各々に対する音響特徴量成分Ｃ_ω,t,ｄを要素にもつ三次元配列Ｃ＾をサンプリングにより決定する音源分離手段と、
   前記Ｄ個の音響要素ｄの各々について、前記二次元配列Ｚ^、前記三次元配列Ｃ^，及び前記二次元配列Ｕ^に基づいて、各時刻ｔに前記基底の状態が各状態ｋである事後確率、及び各時刻ｔに前記基底の状態が新規の状態である事後確率を算出し、算出した事後確率に従って、前記二次元配列Ｚ^をサンプリングにより決定する状態系列推論手段と、
   （ｋ、ω、ｄ）の全ての組み合わせについて、前記二次元配列Ｚ^及び前記三次元配列Ｃ^に基づいて、パラメータφ_ω,d ^(k)、ψ_ω,d ^(k)を算出し、算出したパラメータφ_ω,d ^(k)、ψ_ω,d ^(k)をパラメータとする確率分布に従って、前記三次元配列Ｈ^をサンプリングにより決定する基底スペクトル推論手段と、
   （ｄ、ｔ）の全ての組み合わせについて、予め与えられた前記Ｌ個の音響イベントｌの各々に対する一部の時刻ｔnの各々の音響イベントの有無を表すパラメータＸ_ｌ,tｎ（＝｛０、１｝）を要素にもつ一次元配列Ｘ_ｌＮ^、前記三次元配列Ｃ^、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記遷移行列ａ_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をサンプリングにより決定する遷移確率推論手段と、
   前記Ｌ個の音響イベントｌの各々について、前記一次元配列Ｘ_ｌＮ^、前記二次元配列Ｕ^、前記一次元配列ｗ_ｌ＾、及び前記パラメータα_ｌに基づいて、事後確率が最大となるように前記一次元配列ｗ_ｌ＾を求める重みベクトル推論手段と、
   前記Ｌ個の音響イベントｌの各々について、前記一次元配列ｗ_ｌ＾に基づく前記パラメータα_lの確率分布に従って、前記パラメータα_ｌをサンプリングにより決定する精度パラメータ推論手段と、
   予め定められた収束条件を満たすまで、前記音源分離手段、前記状態系列推論手段、前記基底スペクトル推論手段、前記アクティベーション推論手段、前記遷移確率推論手段、前記重みベクトル推論手段、及び前記精度パラメータ推論手段による各処理を繰り返し行い、前記三次元配列Ｃ＾、前記二次元配列Ｚ^、前記三次元配列Ｈ^、前記二次元配列Ｕ^、及び前記Ｌ個の音響イベントｌの各々の前記一次元配列ｗ_ｌ＾を出力する収束判定手段と、
   を含む音響信号分析装置。
5. 請求項１〜請求項３の何れか１項記載の音響信号分析方法を構成する各ステップをコンピュータに実行させるためのプログラム。