(57)【要約】【課題】音響信号の時系列データから、そこに含まれる全ての音イベントの音響的特徴とその発音区間を自律的に推定することができるようにする。【解決手段】スライス補助変数推論部１６により補助変数ｓをサンプリングにより決定する。音源分離推論部１８により時間周波数成分Ｃω,t,ｄを要素にもつ三次元配列Ｃ＾を決定する。状態系列推論部２０により基底の状態が何れであるかを表わすパラメータである二次元配列Ｚ^を決定する。基底スペクトル推論部２２により、基底の状態ｋの基底スペクトルを表わすパラメータＨω,d(k)を要素にもつ三次元配列Ｈ^を決定する。アクティベーション推論部２４により、発音の有無を表すパラメータＵd,ｔを要素にもつ二次元配列Ｕ^を求める。遷移確率推論部２６により、パラメータＵd,tが０から１へ遷移する遷移確率ａd及びａD†を決定すると共に、前記パラメータＵd,tが１から１へ遷移する遷移確率ｂdを決定する。【選択図】図４
2014044296.tif 000002
【特許請求の範囲】
【請求項１】
複数の音イベントが混在する音響信号の時系列データを入力として、観測時間周波数成分Ｙω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙ＾を出力する時間周波数分解手段と、
予め定められたＤ個の音イベントｄの各々における所定の基底の状態ｋの基底スペクトルを表わすパラメータＨω,d(k)を要素にもつ三次元配列Ｈ^、前記Ｄ個の音イベントの各々に対する各時刻ｔの発音の有無を表すパラメータＵd,t（＝｛０、１｝）を要素にもつ二次元配列Ｕ^、及び前記Ｄ個の音イベントｄの各々に対して各時刻ｔに前記基底の状態が何れであるかを表わすパラメータＺｄ,tを要素にもつ二次元配列Ｚ^の各々の初期値を設定すると共に、Ｄ†個（Ｄ†＝Ｄ＋１）の音イベントｄの各々について、前記パラメータＵd,tが０から１へ遷移する遷移確率ａd（ａ1＞ａ2＞・・・＞ａ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をサンプリングにより決定する遷移確率推論手段と、
予め定められた収束条件を満たすまで、前記補助変数推論手段、前記音源分離手段、前記状態系列推論手段、前記基底スペクトル推論手段、前記アクティベーション推論手段、及び前記遷移確率推論手段による各処理を繰り返し行い、前記三次元配列Ｃ＾、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^を出力する収束判定手段と、
を含む音響信号分析装置。
【請求項２】
前記補助変数推論手段は、０からＵd,t＝１となるｔをもつ音イベントｄの中で最も小さい値となる遷移確率ａdまでの一様分布に従って、補助変数ｓをサンプリングにより決定し、
前記遷移確率ａD†が、前記補助変数ｓより大きい場合、追加する音イベントｄの遷移確率ａdを、前記遷移確率ａD†に基づく遷移確率ａdの確率分布に従ってサンプリングにより決定すると共に、前記追加する音イベントｄに対するパラメータＵd,t及び前記所定個の基底の状態ｋを表わすパラメータＨω,d(k)の初期値を設定し、ＤをＤ†に更新すると共に、サンプリングにより決定され、かつ、前記補助変数ｓより小さくなる遷移確率ａｄにおけるｄに、Ｄ†を更新し、
前記遷移確率ａD†が、前記補助変数ｓ以下である場合、ｓ＜ａｄとなる最大の音イベントのインデックスｄをＤ＊＝ｄと更新し、Ｄ†＝Ｄ＊＋１と更新し、音イベントＤ†よりもインデックスの大きいイベントｄは削除する請求項１記載の音響信号分析装置。
【請求項３】
初期値設定手段は、
前記予め定められたＤ個の音イベントｄの各々における基底の状態数を１とし、全てのパラメータＨω,d(1)が非負値であり、全てのパラメータＵd,tが非負値である、という条件の下で、前記パラメータＨω,d(1)と前記パラメータＵd,tとの積を、全てのｄについて足し合わせたモデルについて、前記時間周波数分解手段によって出力された観測時間周波数成分Ｙω,tと前記モデルとの距離を表わした目的関数の値が小さくなるように、前記パラメータＨω,d(1)及び前記パラメータＵd,tの各々を推定し、前記推定された前記パラメータＵd,tの各々を二値化することにより、前記三次元配列Ｈ^及び前記二次元配列Ｕ^の初期値を設定すると共に、前記パラメータＺｄ,tの各々を１とした二次元配列Ｚ^の初期値を設定する請求項１又は２記載の音響信号分析装置。
【請求項４】
時間周波数分解手段、パラメータ初期値設定手段、補助変数推論手段、音源分離手段、状態系列推論手段、基底スペクトル推論手段、アクティベーション推論手段、遷移確率推論手段、及び収束判定手段を含む音響信号分析装置における音響信号分析方法であって、
前記音響信号分析装置は、
前記時間周波数分解手段によって、複数の音イベントが混在する音響信号の時系列データを入力として、観測時間周波数成分Ｙω,t（ωは周波数、ｔは時刻のインデックスである。）を要素にもつ二次元配列Ｙ＾を出力するステップと、
前記パラメータ初期値設定手段によって、予め定められたＤ個の音イベントｄの各々における所定の基底の状態ｋの基底スペクトルを表わすパラメータＨω,d(k)を要素にもつ三次元配列Ｈ^、前記Ｄ個の音イベントの各々に対する各時刻ｔの発音の有無を表すパラメータＵd,t（＝｛０、１｝）を要素にもつ二次元配列Ｕ^、及び前記Ｄ個の音イベントｄの各々に対して各時刻ｔに前記基底の状態が何れであるかを表わすパラメータＺｄ,tを要素にもつ二次元配列Ｚ^の各々の初期値を設定すると共に、Ｄ†個（Ｄ†＝Ｄ＋１）の音イベントｄの各々について、前記パラメータＵd,tが０から１へ遷移する遷移確率ａd（ａ1＞ａ2＞・・・＞ａ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をサンプリングにより決定するステップと、
前記収束判定手段によって、予め定められた収束条件を満たすまで、前記補助変数推論手段、前記音源分離手段、前記状態系列推論手段、前記基底スペクトル推論手段、前記アクティベーション推論手段、及び前記遷移確率推論手段による各処理を繰り返し行い、前記三次元配列Ｃ＾、前記二次元配列Ｚ^、前記三次元配列Ｈ^、及び前記二次元配列Ｕ^を出力するステップと、
を含む音響信号分析方法。
【請求項５】
コンピュータを、請求項１〜請求項３の何れか１項記載の音響信号分析装置の各手段として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音響信号分析装置、方法、及びプログラムに係り、特に、音響信号の時系列から、信号パラメータを分析する音響信号分析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
一般に、音響信号などの大規模なメディアデータを用いたアプリケーションにおいて、そのデータにいかに有効なラベルを付与するか、大量のデータをいかに計算効率良く処理するかが重要な課題となっている。
従来は、教師あり学習の下で、音響信号に含まれる音イベントを自動的にラベリングする技術が一般的であった。すなわち、学習データとして、ラベル付けされた音イベントの音源信号を事前に用意し、それを周波数分析して得られる音響的特徴と統計モデルを用いて、未知の音響信号の時系列データにラベル付けを行った（例えば、非特許文献１、非特許文献２、及び非特許文献３）。統計モデルでは、例えば、ガウス混合モデル（ＧＭＭ） を用いる場合、各音イベントを表現しうる音響的特徴の頻度分布が学習される。隠れマルコフモデル（ＨＭＭ）を用いる場合、音響的特徴の統計的な時間遷移が学習され、これを用いて未知の音響信号の時系列データにラベル付けを行う。これらのラベル付けの多段処理技術（音響信号が無音なのかそうでないか、音楽か話声かを段階的に識別する）も提案されている（例えば、非特許文献４）。また、音の重ね合わせを考慮して、音響信号をあらかじめ、非負値行列分解（ＮＭＦ）によって音源分離し、分離信号を学習データに用いる技術も提案されている（例えば、非特許文献５）。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】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（ａ1＞ａ2＞・・・＞ａ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（ａ1＞ａ2＞・・・＞ａ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.）。すなわち、
【００１５】
【数１】
2014044296.tif 000003
のように観測スペクトログラムＹ＾をＤ個の頻出の基底スペクトルｈｄ＾＝[Ｈ１,ｄ,..., ＨΩ,ｄ]Ｔとそれぞれ基底スペクトルの音量変化を表すアクティベーションｕｄ＾＝[Ｕｄ,１,...,Ｕｄ,Ｔ]で近似することに相当する。ここで、ｈｄ＾とｕｄ＾のペアを音イベントのコンポーネントと呼ぶ。このコンポーネントによって、一つの音イベントが表現されることが望ましいが、実際の音イベントのスペクトルは時間的に変化し、非定常であると言える。基底スペクトルが時間にともなって変化するように拡張したモデルが提案されており、そこでは各基底スペクトルは時刻ｔに、ある一つの状態Ｚｄ,ｔ∈Ｎを取ると見なし、
【００１６】
【数２】
2014044296.tif 000004
と表現する（例えば、文献：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×Ｄの重ね合わせとして表現されていると考えると、
【００１７】
【数３】
2014044296.tif 000005
【００１８】
【数４】
2014044296.tif 000006
のように書ける。上記（４）式で示されるポアソン分布の再生性より、
【００１９】
【数５】
2014044296.tif 000007
が成り立つ。これは前述したように、モデルと観測スペクトログラムの間の距離尺度をKLdivergence とした場合と等価である。二乗誤差規準やItakura-Saito divergence 規準の場合においては、上記（５）式 をガウス分布、複素ガウス分布とすることで同様に議論できる。
【００２０】
[２．スペクトログラムモデルのグラフィカル表現]
スペクトログラムモデルをベイズ的な枠組みの下でグラフィカル表現したものが図３である。スペクトログラムの生成過程を以下、順番に説明する。
【００２１】
[２．１  基底スペクトルの生成（Ｈ＾の事前分布）]
基底スペクトルを生成するための事前分布として、上記（５）式で示されるポアソン分布の共役事前分布であるガンマ分布を利用する。実際、ガンマ分布の事前分布はスパースな解を導くという報告もある。そこで、
【００２２】
【数６】
2014044296.tif 000008
とする。ここで、ｋは基底スペクトルのインデックスであり、音イベントｄに対して合計Ｋｄ個の基底スペクトルを準備する。このとき、基底スペクトル集合Ｈ＾の同時分布は
【００２３】
【数７】
2014044296.tif 000009
と書ける。上記（７）式で示される分布に従って、基底スペクトル集合の各要素を生成する。ここで、ΘＨ＝｛φ1,1,ψ1,1,φ2,1,ψ2,1,..., φΩ,Ｄ,ψΩ,Ｄ}とする。
【００２４】
[２．２  アクティベーションの生成（Ｕ＾の事前分布）]
アクティベーションは各音イベントの発音区間を表現する。音源の持続性を表現するために、アクティベーションの系列のマルコフ性を仮定して、各音イベントごとに遷移行列
【００２５】
【数８】
2014044296.tif 000010
を用意する。つまり、０→０の遷移確率は１−ａｄ，０→１の遷移確率はａｄ，１→０の遷移確率は１−ｂｄ，１→１ の遷移確率はｂｄとし、これらの状態遷移によって、０と１からなる音イベントｄのアクティベーションを生成する。具体的には、Ｕｄ，０＝０として、
【００２６】
【数９】
2014044296.tif 000011
のように上記（９）式で示されるベルヌーイ分布から生成される確率変数として記述できる。このとき、すべての音イベントのアクティベーション集合の同時分布ｐ（Ｕ＾｜ａ＾，ｂ＾）は
【００２７】
【数１０】
2014044296.tif 000012
と書ける。ここで、ａ＾＝{ａ１,...,ａＤ}、ｂ＾＝{ｂ１,...,ｂＤ}とし、ｃｄ00、ｃｄ01、ｃｄ10、ｃｄ11はそれぞれ０→０,０→１,１→０,１→１に遷移する回数とする。さらに、ａｄとｂｄの事前分布をそれらの共役性から、
【００２８】
【数１１】
2014044296.tif 000013
とおくと、ΘＵ＝{α,γ,δ}に対するすべての音イベントのアクティベーションの同時分布は
【００２９】
【数１２】
2014044296.tif 000014
となる。この分布に従って、アクティベーション集合を生成する。
【００３０】
[２．３  基底スペクトルの状態系列の生成（Ｚ＾の事前分布）]
ｄ番目の音イベントの基底スペクトルの状態系列{Ｚｄ,ｔ,...,Ｚｄ,Ｔ}はそれぞれ離散的な値１,...,Ｋｄ(状態のインデックス) を取る。
【００３１】
【数１３】
2014044296.tif 000015
を各基底スペクトルの生起確率とすると、状態系列の同時分布は、
【００３２】
【数１４】
2014044296.tif 000016
となる。ただし、δ（ｘ−ｋ）はクロネッカーのデルタであり、ｘ＝ｋのとき１、それ以外は０とする。ｎｄ(ｋ)はＺｄ,ｔ’＝ｋ（ｔ’＝１,...,Ｔ）を満たすｔ’の個数を表す。ここで、πｄ＾に対して、次のような事前分布を考える。
【００３３】
【数１５】
2014044296.tif 000017
ただし、βｄは正のパラメータとする。生起確率πｄ＾を周辺化すると、βｄに対する状態系列の同時分布は
【００３４】
【数１６】
2014044296.tif 000018
と書ける。したがって、状態系列集合の同時分布は
【００３５】
【数１７】
2014044296.tif 000019
となる。上記（１６）式で示される分布に従って、状態系列集合を生成する。ここで、ΘＺ＝{β１,...,βＤ} とする。
【００３６】
[２．４  モデルパラメータの同時分布]
図３のグラフィカル表現に基づいて、観測スペクトログラムとパラメータＣ＾,Ｚ＾,Ｈ＾,Ｕ＾の同時分布を書き起こすと、
【００３７】
【数１８】
2014044296.tif 000020
となる。ΘＺ,ΘＨ,ΘＵは超パラメータと呼び、Θ＝{ΘＺ,ΘＨ,ΘＵ}とする。
【００３８】
[２．５  観測スペクトログラムＹ＾および分離スペクトログラムＣ＾の対数尤度関数]
観測スペクトログラムＹ＾の対数尤度関数は、
【００３９】
【数１９】
2014044296.tif 000021
となる。ここで、ｒｄ,ｔ（ｋ）＝δ（Ｚｄ,ｔ−ｋ）のようなインジケータを利用する。一方、観測スペクトログラムおよび音源分離された各音イベントのスペクトログラムからなる完全データＹ＾、Ｃ＾の対数尤度関数は
【００４０】
【数２０】
2014044296.tif 000022
となる。
【００４１】
[３．音イベント数、基底スペクトルの状態数の決定方法]
２．で説明したスペクトログラムモデルを利用する上で、考えなければならない問題として、以下の２点を取り上げる。
【００４２】
１．観測スペクトログラムに含まれる音イベントの総数Ｄの決定方法
２．各音イベントを表現するための基底スペクトルの数（状態数）Ｋｄの決定方法
従来は、音イベントの総数は事前に固定し、基底スペクトルの状態数はどの音イベントに対しても同数としたり、対象に応じて（例えば、話し声と音楽）、異なる状態数を与えることが多かった。しかし、一般的には観測される音響信号の事前知識がないことも多く、そこに含まれている音イベントの総数や必要な基底スペクトルの状態数はモデルが大量の音響信号データから自動的に決定されることが理想的である。
【００４３】
本実施の形態では、機械学習の分野で注目されているノンパラメトリックベイズ法を利用して、ＤとＫｄの値を大規模データから自動的に決定する推論アルゴリズムを導出する。具体的には、以下の確率過程に基づいてアクティベーションＵ＾と状態系列Ｚ＾を生成する。
【００４４】
１．音イベント数を観測データから学習するために、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 を構成するため、遷移行列のパラメータａ＾とｂ＾の周辺化は行わない。具体的に、上記（１７）式 の観測データＹ＾とパラメータＣ＾,Ｚ＾,Ｈ＾,Ｕ＾,ａ＾,ｂ＾の同時分布は、
【００４５】
【数２１】
2014044296.tif 000023
と修正する。ここで、ａ＾＝{ａ１,ａ２,...,ａＤ}は、ａ１＞ａ２＞...＞ａＤのように順序付け、Stick-breaking construction より、Ｄ→∞とすることで、音イベント数を観測される音響信号の時系列データから自動的に学習させる。このときａｄは、
【００４６】
【数２２】
2014044296.tif 000024
によって生成される。ここで、Ｄ†をアクティベーション集合Ｕ＾の行数に相当し（Ｕ＾はＤ†×Ｔの行列をイメージされたい）、ａ＾＝{ａ１,ａ２,...,ａＤ†}とする。一方、ｄ＜Ｄ†となるｄは、少なくとも１回はＵｄ，：において１の値を持つ（アクティブとなる）とする。ｂ＾＝{ｂ１,...,ｂＤ†}は、前述と同様に、ｂｄ〜Ｂｅｔａ（γ,δ）に従って生成する。
基底スペクトルの状態数Ｋｄをデータから学習させるためのChinese restaurant processについては、後述する４．３にて具体的に説明する。
【００４７】
[４．スライスサンプリングによるパラメータ推論]
スライスサンプリングと動的計画法を組み合わせて、提案モデルのパラメータを推論する。これは、Stick-breaking constructionの打ち切り数（音イベント数）を予め大きな値に固定するのではなく、スライスサンプリングによって適応的に打ち切り数を選択しながら、パラメータを推論する方法である。まず、スライス補助変数ｓを導入する。
【００４８】
【数２３】
2014044296.tif 000025
ｓが与えられたとき、Ｕ＾の条件付き確率は、
【００４９】
【数２４】
2014044296.tif 000026
となる。ここで、I(A) は，Aが真であるときI(A)＝１となり、それ以外で０となる関数とする。この式は、ａｄ＜ｓとなるすべての音イベントｄのＵｄ，：を０とする。ここで、Ｄ＊をａＤ＊＞ｓとなる最大の音イベントのインデックスとする。
アクティブな音イベントｄ（すなわち，ｄ：∃ｔ，Ｕｄ，ｔ＝１）に対して、ｄ＜Ｄ†となるインデックスＤ†を導入する（Ｄ†自体は非アクティブな音イベントとする）。スライスサンプリングにおいて推論すべきパラメータは、スライス補助変数ｓとＤ†番目までの音イベントとなる。すなわち，＜ｓ，Ｄ＊，Ｄ†，Ｃ：,：,１：Ｄ†，Ｚ１：Ｄ†,：，Ｈ（：）：,１：Ｄ†，Ｕ１：Ｄ†,：，ａ１：Ｄ†，ｂ１：Ｄ†＞であり、尤度関数の値が収束するまで、これらを繰り返しサンプリングする。各パラメータのサンプリング方法を以下に示す。
【００５０】
[４．１  スライス補助変数ｓのサンプリング]
上記（２２）式から、スライス補助変数ｓをサンプリングする。
【００５１】
（ｓ＜ａＤ†の場合）
Ｄ＊＝Ｄ†と更新する。Ｄ＊＜Ｄ†である必要があるので、
【００５２】
【数２５】
2014044296.tif 000027
にしたがって、ｄ＝Ｄ＊＋１についてａｄをサンプリングする（音イベントを追加する）。この分布はlogａｄについて対数凹分布となるため、適応的棄却サンプリング（Adaptive rejection sampling, ARS）を利用できる。ｓ＞ａｄとなるまで、ａｄを繰り返しサンプリングし、ｓ＞ａｄとなった初めてのｄをＤ†とするように更新する。さらに、ここで追加された音イベントｄに関しては、Ｕｄ，：＝０、Ｈ（１）ω，ｄ〜Ｇａｍｍａ（φω，ｄ，ψω，ｄ）としてパラメータを初期化する。
【００５３】
（ｓ≧ａＤ†の場合）
ｓ＜ａｄとなる最大の音イベントのインデックスｄをＤ＊＝ｄと更新する。Ｄ†＝Ｄ＊＋１と更新する。Ｄ＊＜Ｄ†なので音イベントを追加する必要はない。また、音イベントＤ†よりもインデックスの大きいイベントｄは削除する。
【００５４】
［４．２  Ｃω，ｔ，１，... ，Ｃω，ｔ，Ｄ†のサンプリング］
Ｃω，ｔ，１，... ，Ｃω，ｔ，Ｄ†の条件付き確率は，
【００５５】
【数２６】
2014044296.tif 000028
となる。ここで、ｒｄ，ｔ（ｋ）＝δ（Ｚｄ，ｔ−ｋ）とし、ｐω，ｔ，ｄは
【００５６】
【数２７】
2014044296.tif 000029
とする。上記（２５）式で示される多項分布からＣω，ｔ，１，... ，Ｃω，ｔ，Ｄ†をサンプリングする。
【００５７】
［４．３  Ｚｄ，ｔのサンプリング］
４．２の推論によって、音イベントごとの音源分離スペクトログラムＣ＾が得られる。次は、分離された音イベントｄの各時刻のスペクトルＣ：,１,ｄ,...,Ｃ：,Ｔ,ｄが、どの基底スペクトルによって表現されるか、状態系列Ｚｄ,１,...,Ｚｄ,Ｔを推論する。
ｐ（Ｚｄ,１,...,Ｚｄ,Ｔ｜Ｃ：,１,ｄ,...,Ｃ：,Ｔ,ｄ）を最大化する基底スペクトルの状態系列Ｚｄ,１,...,Ｚｄ,Ｔを求める際、単純に、Ｚｄ,１,...,Ｚｄ,Ｔの全ての可能な値で評価するのは計算量的に非現実である。そこで、ｐ（Ｚｄ,１,...,Ｚｄ,Ｔ｜Ｃ：,１,ｄ,...,Ｃ：,Ｔ,ｄ）からＺｄ,１,...,Ｚｄ,Ｔをサンプリングすることによって、ｐ（Ｚｄ,１,...,Ｚｄ,Ｔ｜Ｃ：,１,ｄ,...,Ｃ：,Ｔ,ｄ）の値の大きなＺｄ,１,...,Ｚｄ,Ｔを確率的に求める。このサンプリングを効率良く実現する方法として、Gibbsサンプリングを利用する。
Gibbsサンプリングでは、ｐ（Ｚｄ,１,...,Ｚｄ,Ｔ｜Ｃ：,１,ｄ,...,Ｃ：,Ｔ,ｄ）からＺｄ,１,...,Ｚｄ,Ｔを同時にサンプリングするのではなく、逐次的にサンプリングする。すなわち、Ｚｄ,ｔをサンプリングする際、Ｚｄ,ｔ以外の{Ｚｄ,１,..., Ｚｄ,ｔ−１, Ｚｄ,ｔ+１,..., Ｚｄ,Ｔ}（以降Ｚｄ, ＼ｔ＾と記述する）の値は既知とした，ｐ（Ｚｄ,ｔ｜Ｚｄ,＼ｔ＾,Ｃ：,ｔ,ｄ,Θ）より、Ｚｄ,ｔをサンプリングする。
【００５８】
Ｚｄ，ｔの条件付き確率は、Ｙ＾の尤度関数を考慮して，
【００５９】
【数２８】
2014044296.tif 000030
と書ける。ここで、
【００６０】
【数２９】
2014044296.tif 000031
とする。上記（２７）式に従って、既存の状態ｋの各々について条件付き確率を算出し、新規の状態ｋnewについて条件付き確率を算出し、算出された各状態の条件付き確率により、Ｚｄ，ｔをサンプリングして決定する。
【００６１】
上式から分かるように、状態数Ｋｄを無限極限としたスペクトログラムモデルにおいて、各時刻に用いられる状態に着目すると、他の時刻に多く用いられている状態ほど使われやすくなる性質がある。また、新しい状態が用いられやすくなるか否かについてはパラメータβｄが影響している。このような状態系列の生成方法はChinese restaurant process (CRP)と呼ばれ、ディリクレ混合過程の一構成を与える。CRPの重要な性質である交換可能性（Exchangeability）より、Ｚｄ,ｔの任意のｔの入れ替えをおこなって出現順序を変えても結果は変わらない。
【００６２】
［４．４  Ｈω，ｄ（ｋ）のサンプリング］
Ｈω，ｄ（ｋ）の条件付き確率は、
【００６３】
【数３０】
2014044296.tif 000032
となる。上記（２９）式で示されるガンマ分布からＨω，ｄ（ｋ）をサンプリングする。
【００６４】
［４．５  Ｕｄ，１,...,Ｕｄ，Ｔのサンプリング］
Forward-filtering backward-sampling アルゴリズムを利用して、Ｕｄ，１,...,Ｕｄ，Ｔを推論する。
このとき、Ｕ＾のｄ行以外の値はすべて固定する。まず、ｔ＝１,...,Ｔに対して、
【００６５】
【数３１】
2014044296.tif 000033
を再帰的に計算する。次に、ｐ（Ｕｄ，Ｔ｜Ｙ：，１：Ｔ，Ｃ：，１：Ｔ，ｄ，Ｚｄ，：，Ｈ（：）：，ｄ）からＵｄ，Ｔをサンプリングする。そして、ｔ＝Ｔ−１,...,１に対して、Ｕｄ，ｔ＋１が与えられた下で、
【００６６】
【数３２】
2014044296.tif 000034
に従って、Ｕｄ，ｔを後方から順番にサンプリングする。
【００６７】
［４．６  ａｄのサンプリング］
ｄ＝１,...,Ｄ†−１対して、ａｄの条件付き確率は、
【００６８】
【数３３】
2014044296.tif 000035
となり、上記（３２）式で示されるベータ分布からａｄをサンプリングする。ｄ＝Ｄ†の音イベントは非アクティブとなるが、スライスサンプリングの打ち切り数のために、ａＤ†を計算する必要がある。ａＤ†の条件付き確率は，上記（２４）式をｄ＝Ｄ†としたときであり、適応的棄却サンプリングを用いて、ａＤ†をサンプリングする。
【００６９】
［４．７  ｂｄのサンプリング］
ｄ＝１,...,Ｄ†に対して、ｂｄの条件付き確率は、
【００７０】
【数３４】
2014044296.tif 000036
となり、上記（３３）式で示されるベータ分布からｂｄをサンプリングする。
【００７１】
＜システム構成＞
次に、音響信号の信号パラメータを分析して出力する音響信号分析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【００７２】
図４に示すように、本実施の形態に係る音響信号分析装置は、ＣＰＵと、ＲＡＭと、後述する音響信号分析処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【００７３】
本実施の形態に係る音響信号分析装置は、入力部１と、記憶部２と、演算部３と、出力部４とを備えている。また、演算部３は、短時間フーリエ変換部１０と、パラメータ初期値生成部１２と、スライス補助変数推論部１６と、音源分離推論部１８と、状態系列推論部２０と、基底スペクトル推論部２２と、アクティベーション推論部２４と、遷移確率推論部２６と、収束判定部２８とを備えている。また、スライス補助変数推論部１６、音源分離推論部１８、状態系列推論部２０、基底スペクトル推論部２２、アクティベーション推論部２４、遷移確率推論部２６、及び収束判定部２８は、モデルパラメータ更新部１４を構成する。この実施形態は、前述したパラメータ推定アルゴリズムを用いて信号解析を行う構成である。
【００７４】
入力部１には、分析する対象である音響信号の時系列が入力される。記憶部２は、入力部１に入力された音響信号の時系列を記憶する。また、記憶部２は、後述する各処理での結果等を記憶する。
【００７５】
短時間フーリエ変換部１０は、記憶部２に記憶された音響信号の時系列を読み出して、音響信号の振幅スペクトログラムを計算する。このとき、音響信号の時系列に対して時間フレームを設定し、時間フレーム長６４ｍｓ、時間フレームシフト長３２ｍｓとして、時間フレームごとに短時間フーリエ変換を行い、観測時間周波数成分Ｙω,tを各（ω，ｔ）の要素にもつ二次元配列Ｙ＾（＝観測スペクトログラムＹ＾）を生成し、記憶部２に格納する。このとき、周波数のインデックスの最大値Ω＝５１２であり、時間のインデックスの最大値Ｔはフレームの総数に相当する。また、上記（５）式 に示すように、本実施の形態ではポアソン分布に従う確率変数として、スペクトログラムをモデル化するため、Ｙ＾の全ての要素を整数値に丸め込む。
【００７６】
また、記憶部２には、予め定められた超パラメータα＝４、β１,...,βＤ†＝１、γ＝５００、δ＝１、φ1,1,ψ1,1,φ2,1,ψ2,1,...,φΩ,Ｄ†,ψΩ,Ｄ†＝１からなるパラメータが記憶されている。
【００７７】
パラメータ初期値生成部１２は、基底スペクトル集合Ｈ＾、アクティベーション集合Ｕ＾、状態系列Ｚ＾、音イベント数Ｄ、及び遷移確率ａ＾並びにｂ＾の各パラメータの初期値を生成し、記憶部２へ格納する。
基底スペクトル集合Ｈ＾及びアクティベーション集合Ｕ＾については、記憶部２に記憶された二次元配列Ｙ＾に対し、上記（１）式に従って、通常のＮＭＦを適用し、推定されたＨ＾とＵ＾を初期値とし、記憶部２に格納する。具体的には、基底の状態ｋを１のみとし、基底スペクトル集合Ｈ＾の要素である全てのパラメータＨω,d(1)が非負値であり、アクティベーション集合Ｕ＾の要素である全てのパラメータＵd,tが非負値である、という条件の下で、パラメータＨω,d(1)とパラメータＵd,tとの積を、全てのｄについて足し合わせたモデルについて、二次元配列Ｙ＾とモデルとの距離を表わした目的関数の値が小さくなるように、パラメータＨω,d(1)及びパラメータＵ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をサンプリングにより決定し、予め定められた収束条件を満たすまで、各処理を繰り返し行い、三次元配列Ｃ＾、二次元配列Ｚ^、三次元配列Ｈ^、及び二次元配列Ｕ^を出力することにより、音響信号の時系列データから、そこに含まれる音イベントの数、音イベントの状態数を自動的に決定して、全ての音イベントの音響的特徴とその発音区間を自律的に推定することができる。
【０１０３】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１０４】
また、上述の音響信号分析装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０１０５】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【０１０６】
１     入力部
３     演算部
４     出力部
１０   短時間フーリエ変換部
１２   パラメータ初期値生成部
１４   モデルパラメータ更新部
１６   スライス補助変数推論部
１８   音源分離推論部
２０   状態系列推論部
２２   基底スペクトル推論部
２４   アクティベーション推論部
２６   遷移確率推論部
２８   収束判定部
【図１】
2014044296.tif 000037
【図２】
2014044296.tif 000038
【図３】
2014044296.tif 000039
【図４】
2014044296.tif 000040
【図５】
2014044296.tif 000041
【図６】
2014044296.tif 000042
【図７】
2014044296.tif 000043
【図８】
2014044296.tif 000044
【図９】
2014044296.tif 000045
【図１０】
2014044296.tif 000046