JP5924968B2 - 楽譜位置推定装置、及び楽譜位置推定方法 - Google Patents

Description

本発明は、楽譜位置推定装置、及び楽譜位置推定方法に関する。

従来から、楽曲を書き表した楽譜において演奏されている楽曲の位置（演奏位置）を推定する楽譜位置推定技術が開発されている。楽譜位置推定技術は、自動伴奏、自動合奏、楽譜への演奏位置の重畳表示を実現するための応用が期待されている。ここで、演奏されている楽曲の音響信号の周波数特性と楽譜情報に基づいて合成した楽曲の周波数特性を照合し、演奏されている楽曲について楽譜の位置を実時間で推定する技術が提案されている

例えば、特許文献１に記載の楽譜位置推定技術では、楽譜情報に基づく楽曲の周波数特性を調和混合正規分布モデル（ｈａｒｍｏｎｉｃ Ｇａｕｓｓｉａｎ Ｍｉｘｔｕｒｅ Ｍｏｄｅｌ；ｈａｒｍｏｎｉｃ ＧＭＭ、又は単にＧＭＭともいう。）に基づくテンプレートを用いて合成し、音響信号の周波数特性と照合していた。このテンプレートでは、高調波の振幅を、予め定めた実数のべき乗に定められていた。

特開２０１１−１８０５９０号公報

しかし、特許文献１記載の楽譜位置推定技術では、合成した楽曲の周波数特性と音響信号の周波数特性が合致せず、演奏位置を頑健に推定できないことがあった。
本発明は上記の点に鑑みてなされたものであり、演奏された楽曲の楽譜上の位置をより頑健に推定できる楽譜位置推定装置及び楽譜位置推定方法を提供することを課題としている。

（１）本発明は上記の課題を解決するためになされたものであり、本発明の一態様は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する周波数特性分析部と、楽譜情報が表す楽譜の位置毎の音階に基づく周波数特性であって少なくとも１つの調波構造を含む第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記調波構造に含まれる調波成分毎の前記調波構造への寄与を表す変数の確率分布と、前記調波構造毎の前記第２の周波数特性への寄与を表す変数の確率分布に基づいて前記第１の周波数特性についての尤度を最大にするように、前記楽譜の位置毎に算出する関連度算出部と、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する楽譜位置探索部と、を備えることを特徴とする楽譜位置推定装置である。

（２）本発明のその他の態様は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する周波数特性分析部と、楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記第２の周波数特性に含まれる周波数毎の振幅が表す度数が、対応する周波数における前記第１の周波数特性の振幅が表す度数だけ発生する確率分布に基づいて前記楽譜の位置毎に算出する関連度算出部と、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する楽譜位置探索部と、を備えることを特徴とする楽譜位置推定装置である。

（３）本発明のその他の態様は、前記楽譜位置探索部は、前記音階が変化する楽譜の位置に対応する時刻を含む観測区間における前記音響信号の自己相関に基づいて拍間隔を定め、定めた拍間隔に基づいて探索する楽譜の位置を更新することを特徴とする。

（４）本発明のその他の態様は、楽譜位置推定装置における方法であって、前記楽譜位置推定装置は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する過程と、前記楽譜位置推定装置は、楽譜情報が表す楽譜の位置毎の音階に基づく周波数特性であって少なくとも１つの調波構造を含む第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記調波構造に含まれる調波成分毎の前記調波構造への寄与を表す変数の確率分布と、前記調波構造毎の前記第２の周波数特性への寄与を表す変数の確率分布に基づいて前記第１の周波数特性についての尤度を最大にするように、前記楽譜の位置毎に算出する過程と、前記楽譜位置推定装置は、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する過程と、を有することを特徴とする楽譜位置推定方法。

（５）本発明のその他の態様は、楽譜位置推定装置における方法であって、前記楽譜位置推定装置は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する過程と、前記楽譜位置推定装置は、楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記第２の周波数特性に含まれる周波数毎の振幅が表す度数が、対応する周波数における前記第１の周波数特性の振幅が表す度数だけ発生する確率分布に基づいて前記楽譜の位置毎に算出する過程と、前記楽譜位置推定装置は、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する過程と、を有することを特徴とする楽譜位置推定方法である。

（１）、（２）、（４）、（５）に記載した態様によれば、音響信号に基づく第１の周波数特性と、楽譜情報における音階に基づく第２の周波数特性が完全に合致していなくとも、両者の関連性を検出することができる。そのため、演奏された音楽に対する楽譜の位置をより頑健に推定することが可能になる。
（１）、（４）に記載した態様によれば、楽譜情報に基づく第２の周波数成分に含まれる調波構造の周波数特性の変化、調波構造間の結合特性の変化に柔軟に対応できるので、音響信号に基づく第１の周波数特性との関連性をより的確に推定することができる。そのため、演奏された音楽に対する楽譜の位置をより正確に推定することが可能になる。
（２）、（５）に記載した態様によれば、音響信号に基づく第１の周波数特性の周波数毎の振幅と楽譜情報における音階に基づく第２の周波数特性の周波数毎の振幅が合致する度合いを鋭敏に検出することができる。そのため、演奏された音楽に対する楽譜の位置をより正確に推定することが可能になる。
（３）に記載した態様によれば、楽譜情報における音階が変化する楽譜の位置と、音響信号の振幅が変化する時刻が対応付けられる。そのため、周波数特性の周期性を表す拍間隔を正確に検知でき、ひいては演奏された音楽に対する楽譜の位置をより正確に推定することができる。

本発明の第１の実施形態に係る楽譜位置推定装置の構成を示す概略図である。 粒子フィルタリング法の処理を表す概念図である。 本実施形態に係る楽譜位置推定処理を表すフローチャートである。 楽譜情報と音響信号の関係の一例を表す概念図である。 本実施形態に係る粒子遷移処理を表すフローチャートである。 本実施形態に係る音源モデルの一例を表す概念図である。 本実施形態に係る調波構造の一例を表す図である。 本実施形態に係るヒストグラムの一例を表す図である。 本実施形態に係る周波数特性重み算出部の一構成例を表す概略図である。 本実施形態における周波数特性重み算出処理の一例を表すフローチャートである。 本発明の第２の実施形態に係る楽譜位置推定装置の構成を示す概略図である。 本実施形態に係る周波数特性重み算出部の構成を表す概略図である。 本実施形態における周波数特性重み算出処理を表すフローチャートである。 本発明の第３の実施形態に係る楽譜位置推定装置の構成を示す概略図である。 本実施形態に係る楽譜位置推定処理を表すフローチャートである。 楽譜位置の推定誤差の一例を表す図である。 楽譜位置の推定誤差のその他の例を表す図である。 楽譜位置の推定誤差のその他の例を表す図である。 音響信号と楽譜情報との関連性の一例を表す図である。

（第１の実施形態）
以下、図面を参照しながら本発明の第１の実施形態について説明する。
図１は、本実施形態に係る楽譜位置推定装置１の構成を示す概略図である。
楽譜位置推定装置１は、音響信号入力部１０１、音響特徴量生成部１０２、楽譜情報記憶部１１１、楽譜情報入力部１１２、関連度算出部１２１、楽譜位置探索部１３１、及び演奏情報出力部１４１を含んで構成される。

音響信号入力部１０１は、演奏された音楽を表す音波をディジタル音響信号に変換し、変換したディジタル音響信号を音響特徴量生成部１０２に出力する。音響信号入力部１０１は、例えばマイクロホンとアナログ・ディジタル（Ａｎａｌｏｇ−ｔｏ−Ｄｉｇｉｔａｌ；Ａ／Ｄ）変換部（図示せず）とを含んで構成される。このマイクロホンは、人間が聴取することができる周波数帯域（例えば、２０Ｈｚ−２０ｋＨｚ）の音波をアナログ音響信号に変換し、変換したアナログ音響信号をＡ／Ｄ変換部に出力する。Ａ／Ｄ変換部は、マイクロホンから入力されたアナログ音響信号をディジタル音響信号に変換し、変換したディジタル音響信号を音響特徴量生成部１０２に出力する。ここで、Ａ／Ｄ変換部は、入力されたアナログ音響信号を、例えば、サンプリング周波数４４．１ｋＨｚ、振幅を１６ビットの２進（バイナリ；ｂｉｎａｒｙ）データにパルス符号変調（Ｐｕｌｓｅ Ｃｏｄｅ Ｍｏｄｕｌａｔｉｏｎ；ＰＣＭ）して、サンプル時刻毎の量子化された振幅を表すディジタル音響信号に変換する。以下の説明では、入力されたディジタル音響信号を入力信号と呼ぶことがある。

音響特徴量生成部１０２は、音響信号入力部１０１から入力されたディジタル音響信号から、その物理的な特徴を表す音響特徴量を生成し、生成した音響特徴量を関連度算出部１２１及び楽譜位置探索部１３１に出力する。音響特徴量生成部１０２は、音響特徴量として、例えば、振幅周波数特性として周波数毎の振幅を表す音響スペクトログラム、自己相関、音響スペクトログラムの時間差分に基づく距離値を算出する。音響特徴量生成部１０２は、周波数特性分析部１０３及び相関算出部１０４を含んで構成される。

周波数特性分析部１０３は、音響信号入力部１０１から入力されたディジタル音響信号を時間領域信号から周波数領域信号に変換する。ここで、周波数特性分析部１０３は、例えば、予め定められた個数（例えば、２０４８個）サンプルからなるフレーム（以下、音響フレームとも呼ぶ）ごとにディジタル音響信号を高速フーリエ変換（Ｆａｓｔ Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ；ＦＦＴ）を行って周波数領域信号に変換する。周波数特性分析部１０３は、予め設定された音響フレーム時刻毎（音響フレーム間隔Δｔは、例えば１０ｍｓ）に音響フレームを移動させる。音響フレームを移動させることによって、所定数の新たな信号のサンプル（例えば、サンプリング周波数が４４．１ｋＨｚの場合、４４１サンプル）が音響フレームに含まれ、最も時刻が経過した同数の信号のサンプルが音響フレームから除外される。

周波数特性分析部１０３は、時刻ｔにおける周波数領域信号に基づいてスペクトログラム（以下、音響スペクトログラムと呼ぶ）Ａ_ｆ，ｔを算出する。音響スペクトログラムＡ_ｆ，ｔは、周波数領域信号の周波数ｆ毎の振幅（絶対値）Ａ_ｆ，ｔである。以下の説明では、音響スペクトログラムＡ_ｆ，ｔを、入力音響振幅と呼ぶことがある。周波数特性分析部１０３は、算出した音響スペクトログラムＡ_ｆ，ｔを相関算出部１０４、及び関連度算出部１２１の周波数特性重み算出部１２２に出力する。

相関算出部１０４は、周波数特性分析部１０３から入力された音響スペクトログラムＡ_ｆ，ｔに基づいて、拍間隔ｂについての自己相関を算出する。自己相関は、時刻ｔにおける音響スペクトログラムＡ_ｆ，ｔと、拍間隔ｂだけ過去の時刻ｔ−ｂにおける音響スペクトログラムＡ_{ｆ，ｔ−ｂ}の類似性を表す指標値である。拍間隔ｂとは、演奏される楽曲の時間単位である拍と、その拍に後続する拍との間の時間間隔を表す。相関算出部１０４は、自己相関として、例えば式（１）を用いてＲ（ｂ，｛Ａ_ｋ｝）を算出する。

式（１）において、｛Ａ_ｋ｝は、観測区間Ｔ_ｋにおける現フィルタリングステップｋまでの音響スペクトログラムＡ_ｆ，ｔの集合（セット）を表す。Ｔ_ｋは、観測区間｛τ｜ｋΔＴ−Ｗ＜τ≦ｋΔＴ｝を表す。τは、時刻を表す。ΔＴは、例えば後述する楽譜位置探索におけるフィルタリング間隔（例えば、１秒）である。Ｗは、観測区間長（例えば、２．５秒）である。フィルタリングステップｋとは、後述する粒子フィルタ法（パーティクルフィルタリング、ｐａｒｔｉｃｌｅ ｆｉｌｔｅｒｉｎｇ）を応用した楽譜位置推定における処理の繰り返しを表す整数である。
相関算出部１０４は、算出した自己相関を関連度算出部１２１の拍間隔重み算出部１２３及び楽譜位置探索部１３１の粒子遷移部１３２に出力する。相関算出部１０４は、音響スペクトログラムＡ_ｆ，ｔを粒子遷移部１３２に出力する。

相関算出部１０４は、自己相関Ｒ（ｂ，｛Ａ_ｋ｝）の代わりに、自己相関Ｒ（ｂ，｛Ξ_ｋ’｝）を算出してもよい。｛Ξ_ｋ’｝は、距離値Ξ_ｆ，ｔを聴感補正した距離値Ξ_ｆ，ｔ’の集合（セット）である。距離値Ξ_ｆ，ｔは、現時刻ｔにおける音響スペクトログラムＡ_ｆ，ｔと直前の音響フレームの時刻ｔ−Δｔにおける音響スペクトログラムＡ_{ｆ，ｔ―Δｔ}との差分の絶対値である。
相関算出部１０４は、距離値として、例えば式（２）で表されるΞ_ｆ，ｔを算出してもよい。

式（２）において、Δφ_ｆ，ｔは、式（３）で表される。

式（３）において、φ_ｆ，ｔ等は、対応する音響スペクトログラムＡ_ｆ，ｔの絶対位相（アンラップト位相［ｕｎｗｒａｐｐｅｄ ｐｈａｓｅ］）である。
相関算出部１０４は、線形領域で表された周波数ｆをメル尺度で表された周波数ｆ^ｍｅｌに変換する。線形領域で表された周波数ｆと、メル尺度で表された周波数ｆ^ｍｅｌとの関係は、式（４）で表される。

相関算出部１０４は、フレーム毎の周波数の次元数（例えば、１０２４）をＮ_ｍ（Ｎ_ｍは１よりも大きい整数、例えば、６４である。）に減少させる。そのため、相関算出部１０４は、距離値Ξ_ｆ，ｔを中心周波数ｆ_ｍ ^ｍｅｌ毎にフィルタリングして、聴感補正した距離値Ξ_ｆ，ｔ’を算出する。中心周波数ｆ_ｍ ^ｍｅｌは、式（５）に示されるように、０からナイキスト周波数ｆ_Ｎｙｑ ^ｍｅｌまでの間を等分する周波数である。

式（５）において、ｍは、０より大きくＮ_ｍと等しいか、Ｎ_ｍよりも小さい整数である。フィルタリングに用いる窓関数Ψ_ｍ（ｆ^ｍｅｌ）は、例えば、式（６）で表される。

式（６）は、窓関数Ψ_ｍ（ｆ^ｍｅｌ）の周波数特性は、各中心周波数ｆ_ｍ ^ｍｅｌにおいてピーク１をとり、隣接する中心周波数ｆ_ｍ−１ ^ｍｅｌ，ｆ_ｍ＋１ ^ｍｅｌにおいてゼロとなる三角形を表す。但し、ｍ＝Ｎ_ｍの場合、窓関数ψ_Ｎｍ（ｆ^ｍｅｌ）は、周波数ｆ_Ｎｍ−１ ^ｍｅｌ≦ｆ^ｍｅｌ＜ｆ_Ｎｍ ^ｍｅｌの区間のみ定められる。周波数ｆ_Ｎｍ ^ｍｅｌは、ナイキスト周波数ｆ_Ｎｙｑ ^ｍｅｌであるため、この周波数を越える周波数成分が存在しないからである。
相関算出部１０４は、例えば式（７）を用いて距離値Ξ_ｆ，ｔをフィルタリングする。

相関算出部１０４は、聴感補正した距離値Ξ_ｆ，ｔ’に基づいて、例えば式（８）を用いて自己相関Ｒ（ｂ，｛Ξ_ｋ’｝）を算出する。

距離値Ξ_ｆ，ｔ、又はΞ_ｆ，ｔ’は、音響スペクトログラムＡ_ｆ，ｔの拍間隔ｂに対する周期性を表す指標値である。例えば、定常的な音では、距離値Ξ_ｆ，ｔ、又はΞ_ｆ，ｔ’は、ゼロ又はゼロに近似する。そのため、自己相関Ｒ（ｂ，｛Ξ_ｋ’｝）を用いることで、立ち上がり部分（オンセット）のように、入力信号の振幅の変動が著しい部分に基づいて、振幅周波数特性の周期性を検知することが容易になる。
相関算出部１０４は、算出した距離値Ξ_ｆ，ｔ及び自己相関Ｒ（ｂ，｛Ξ_ｋ’｝）を粒子遷移部１３２に出力する。

楽譜情報記憶部１１１は、楽譜情報を楽曲毎に予め記憶されている。楽譜情報は、楽譜フレーム（ｓｃｏｒｅ ｆｒａｍｅ）ｐ（ｐは、整数）毎に、楽曲を構成する音階毎の基本周波数μ_ｐ ^ｌを要素とする基本周波数ベクトル{μ_ｐ}＝［μ_ｐ ^１，μ_ｐ ^２，…，μ_ｐ ^Ｌｐ］^Ｔで表される。ここで、Ｌ_ｐは、楽譜フレームｐにおける音階の数を表す。ｌは０より大きく、Ｌ_ｐと等しいかＬ_ｐよりも小さい整数である。Ｔは、ベクトル又は行列の転置を表す。つまり、基本周波数ベクトルμ_ｐは、同時に演奏される音階の周波数の組、即ちコード（和音、ｃｈｏｒｄ）を表す情報である。但し、本実施形態では、コードは、１個又は１個より多い数の音階を含んでいればよい。例えば、ある楽譜フレームｐ１において音階Ａ４、Ｃ４、及びＥ４を含むコードを表す基本周波数ベクトルμ_ｐ１は、要素数が３個のベクトル［４４０，５２３，６５９］となる。ここで、音階Ａ４、Ｃ４、Ｅ４の周波数は、それぞれ４４０、５２３、６５９Ｈｚである。

楽譜情報は、楽曲毎の演奏テンポｂ^ｓ’を表すテンポ情報を含んでいでもよい。テンポ情報は、例えばテンポの逆数である拍間隔ｂ^ｓ’で表される。
楽譜情報は、コードを示すコード情報ｎ_ｐ ^ｌを要素とするコードベクトル{ｎ_ｐ}＝［ｎ_ｐ ^１，ｎ_ｐ ^２，…，ｎ_ｐ ^Ｌｐ］^Ｔを、基本周波数ベクトル{μ_ｐ}と対応付けて含んでいてもよい。音階Ａ４、Ｃ４、及びＥ４を要素値として含むコードベクトルは、［Ａ４，Ｃ４，Ｅ４］となる。
楽譜フレームは、１つの楽曲の楽譜情報において楽曲を時間分割する単位時間である。例えば、１拍、即ち四分音符の長さが１２フレームである場合、楽譜情報の時間分解能（ｒｅｓｏｌｕｔｉｏｎ）は十六分音符の三分の一である。

楽譜情報入力部１１２は、楽譜情報記憶部１１１から、処理の対象とする楽曲の楽譜情報を読み出す。楽譜情報入力部１１２は、読み出した楽譜情報を関連度算出部１２１の周波数特性重み算出部１２２、及び楽譜位置探索部１３１の粒子遷移部１３２に出力する。

関連度算出部１２１は、音響特徴量生成部１０２から入力された音響特徴量と楽譜情報入力された楽譜情報が表す楽譜特徴量との関連度を表す重み値を楽譜位置毎に算出する。楽譜特徴量は、例えば、楽譜情報に基づいて演奏されたと仮定された場合における音響信号の振幅周波数特性である。楽譜位置とは、楽譜の冒頭からの拍数や楽譜フレーム数であり、演奏されている楽曲の楽譜上の拍数や楽譜フレーム数を指すこともある。関連度算出部１２１は、算出した楽譜位置毎の重み値を表す重み情報を楽譜位置探索部１３１及び演奏情報出力部１４１に出力する。
関連度算出部１２１は、例えば、周波数特性重み算出部１２２、拍間隔重み算出部１２３、及び粒子重み算出部１２４を含んで構成される。粒子重み算出部１２４は、上述の重み情報として、粒子重み情報を含んだ粒子情報を楽譜位置探索部１３１及び演奏情報出力部１４１に出力する。

周波数特性重み算出部１２２は、周波数特性分析部１０３から入力された音響スペクトログラムと楽譜情報入力部１１２から入力された楽譜情報に含まれる基本周波数ベクトルが表す周波数特性との楽譜位置毎の関連度（例えば、類似度）を表す周波数特性重み値を算出する。周波数特性重み算出部１２２は、算出した周波数特性重み値を表す周波数特性重み情報を粒子重み算出部１２４に出力する。周波数特性重み算出部１２２のより詳細な構成については後述する。
拍間隔重み算出部１２３は、相関算出部１０４から入力された自己相関に基づいて拍間隔毎に音響信号の振幅周波数特性の関連度を表す拍間隔重み値を算出する。拍間隔重み算出部１２３は、拍間隔重み値を表す拍間隔重み情報を粒子重み算出部１２４に出力する。拍間隔重み算出部１２３のより詳細な構成については後述する。

粒子重み算出部１２４は、後述する粒子遷移部１３２から入力された粒子ｉ毎の粒子情報、粒子分布値、状態遷移確率、周波数特性重み算出部１２２から入力された周波数特性重み情報と拍間隔重み算出部１２３から入力された拍間隔重み情報に基づいて、楽譜位置及び拍間隔毎の重み情報を粒子ｉ毎に生成する。
粒子とは、楽譜位置を表す楽譜位置情報、拍間隔を表す拍間隔情報及び入力音響信号に対する観測値である重み情報の組を表す情報の単位である。粒子重み算出部１２４が生成した重み情報を粒子重み情報と呼び、楽譜位置情報、拍間隔情報及び粒子重み情報の組を粒子情報と呼ぶ。粒子情報は、後述するように粒子フィルタリング法により、楽譜位置、拍間隔の推定値を定めるために用いられる。
粒子フィルタリング法を用いる場合、粒子重み算出部１２４は、例えば、式（９）を用いて粒子重み情報として粒子ｉ毎の重み値（以下、粒子重み値と呼ぶ）ｗ_ｉ，ｋを算出する。

式（９）において、ｗ^ｓｐ _ｉ，ｋは、周波数特性重み値を表し、ｗ^ｔｍ _ｉ，ｋは、拍間隔重み値を表す。ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜｛Ａ_ｋ｝，ｂ^ｓ’，ｏ_ｐ）、ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜ｐ^ｉ _ｋ−１，ｂ^ｉ _ｋ−１）は、それぞれ粒子ｉの粒子分布値、状態遷移確率を表す。粒子重み算出部１２４は、算出した粒子重み値ｗ_ｉ，ｋを表す粒子重み情報を生成する。
粒子重み算出部１２４は、入力された粒子ｉ毎の粒子情報に含まれる粒子重み情報を、生成した粒子重み情報に置き換えることで粒子情報を更新する。
粒子重み算出部１２４は、更新した粒子情報を楽譜位置探索部１３１の再標本化部１３３、演奏情報出力部１４１の楽譜位置算出部１４２及び拍間隔算出部１４３に出力する。

楽譜位置探索部１３１は、粒子重み算出部１２４から入力された楽譜位置毎の重み情報に基づいて楽譜位置を探索する。楽譜位置探索部１３１は、例えば、重み情報が表す重み値が音響信号と楽譜情報との関連性が最も高い楽譜位置を探索する。楽譜位置探索部１３１は、探索した楽譜位置を表す楽譜位置情報を関連度算出部１２１に出力する。

楽譜位置探索部１３１は、例えば、粒子フィルタリング法を用いて楽譜位置を探索する。楽譜位置探索部１３１は粒子遷移部１３２及び再標本化部１３３を含んで構成される。ここで、粒子フィルタリング法について説明する。
図２は、粒子フィルタリング法の処理を表す概念図である。
粒子フィルタリング法は、Ｉ．遷移、ＩＩ．重み算出、ＩＩＩ．位置推定・再標本化、の各過程を有する。図２は、Ｉ．遷移、ＩＩ．重み算出、ＩＩＩ．位置推定・再標本化、の各過程における処理の概念を左から右へ順に表す。

図２のＩ．遷移において、横軸が時刻、即ち楽譜位置を表し、縦軸が拍間隔を表す。塗りつぶした円は、それぞれ粒子を表す。各矢印は、起点にある粒子が終点に遷移することを表す。ここで、粒子遷移部１３２が、各時刻、各拍間隔にわたって分布している複数の粒子から各粒子を１つずつ抽出し、抽出した粒子毎の楽譜位置を、その拍間隔に基づいて更新することを表す。

ＩＩ．重み算出では、上段に音響信号、中段に楽譜情報、下段に粒子をそれぞれ表す。各段ともに横軸は時刻、即ち楽譜位置を表す。中段において縦軸は音階を表す。上段と中段にそれぞれ表された四角形の幅は、観測区間を表す。上段の四角形で囲まれる音響信号が、中段の四角形で囲まれる楽譜情報に対して観測対象となる音響信号を表す。下段に表される塗りつぶしの円は、粒子を表す。円の半径は、観測対象となる音響信号と対応する楽譜情報に基づく粒子重みの大きさを表す。即ち、関連度算出部１２１は、粒子毎の楽譜位置及び拍間隔の妥当性を表す粒子重みを、観測区間内の音響信号の音響特徴量と楽譜情報に基づいて算出する。

ＩＩＩ．位置推定・再標本化では、上段に楽譜情報、中段に再標本化前の粒子、下段に再標本化後の粒子を表す。各段ともに横軸は時刻、即ち楽譜位置を表す。中段において、下向きの矢印は粒子の分布に基づいて推定された推定楽譜位置を表す。即ち、後述する楽譜位置算出部１４２が、粒子毎の楽譜位置を各自の粒子重みで重み付け加算して楽譜位置を推定することを表す。
中段の粒子を起点とし、下段の粒子を終点とする矢印は、中段の粒子を再標本化することを表す。濃く塗りつぶした粒子を終点とする太線の矢印は、起点の粒子の粒子重みが大きいために複数の粒子に分割・複製することを表す。×印を終点とする破線の矢印は、起点の粒子の粒子重みが小さいために、その粒子を棄却することを表す。即ち、再標本化部１３３が、粒子重みが大きい粒子を分割・複製し、粒子重みが小さい粒子ほど優先して棄却することを表す。これにより、楽譜位置及び拍間隔が妥当な粒子が残される。その後、Ｉ．遷移に戻り、Ｉ−ＩＩＩの過程を繰り返す。粒子遷移部１３２、再標本化部１３３の構成もしくは動作については後述する。

図１に戻り、演奏情報出力部１４１は、関連度算出部１２１から入力された関連度情報として粒子情報に基づいて楽譜位置を算出し、算出した楽譜位置を表す楽譜位置情報を生成する。演奏情報出力部１４１は、生成した楽譜位置情報を含む演奏情報を楽譜位置推定装置１の外部に出力する。
演奏情報出力部１４１は、例えば、楽譜位置算出部１４２、拍間隔算出部１４３、信頼度判定部１４４、及び楽譜位置出力部１４５を含んで構成される。

楽譜位置算出部１４２は、粒子重み算出部１２４から入力された粒子ｉ毎の粒子情報から楽譜位置情報と粒子重み情報を抽出する。楽譜位置算出部１４２は、粒子ｉ毎に抽出した楽譜位置情報が表す楽譜位置ｐ_ｉ ^ｋに基づいて推定楽譜位置＜ｐ_ｋ＞を算出する。楽譜位置算出部１４２は、例えば、楽譜位置ｐ_ｉ ^ｋの単純平均を推定楽譜位置＜ｐ_ｋ＞としてもよいし、粒子重み情報が表す粒子重み値ｗ_ｉ，ｋによる加重平均を推定楽譜位置＜ｐ_ｋ＞としてもよい。楽譜位置算出部１４２が加重平均を算出する際、必ずしも全ての粒子ｉを考慮する必要はなく、一部の粒子のみを考慮してもよい。一部の粒子とは、例えば、粒子重み値ｗ_ｉ，ｋが最も大きい粒子から粒子重み値ｗ_ｉ，ｋの大きい順に予め定めた個数Ｉ_ｖ（例えば、粒子数Ｉの２％）までの粒子である。これにより、楽譜位置ｐ_ｉ ^ｋの信頼性が低いことを表す粒子重み値ｗ_ｉ，ｋがより小さい粒子が無視されるため、楽譜位置の推定精度が向上する。楽譜位置算出部１４２は、算出した推定楽譜位置＜ｐ_ｋ＞を表す推定楽譜位置情報を信頼度判定部１４４に出力する。
楽譜位置算出部１４２は、一部の粒子を考慮して推定楽譜位置＜ｐ_ｋ＞を算出する際、一部の粒子についての粒子重み値の総和の他に、さらに全部の粒子についての粒子重み値の総和を算出する。楽譜位置算出部１４２は、一部の粒子についての粒子重み値の総和、全部の粒子についての粒子重み値の総和を信頼度判定部１４４に出力する。

拍間隔算出部１４３は、粒子重み算出部１２４から入力された粒子ｉ毎の粒子情報から拍間隔情報と粒子重み情報を抽出する。拍間隔算出部１４３は、粒子ｉ毎に抽出した拍間隔情報が表す拍間隔ｂ_ｉ ^ｋに基づいて推定拍間隔＜ｂ_ｋ＞を算出する。拍間隔算出部１４３は、例えば、拍間隔ｂ_ｉ ^ｋの単純平均を推定拍間隔＜ｂ_ｋ＞としてもよいし、粒子重み情報が表す粒子重み値ｗ_ｉ，ｋによる加重平均を推定拍間隔＜ｂ_ｋ＞としてもよい。拍間隔算出部１４３が加重平均を算出する際、必ずしも全ての粒子ｉを考慮する必要はなく、楽譜位置算出部と同様に一部の粒子のみを考慮してもよい。これにより、推定拍間隔＜ｂ_ｋ＞の信頼性が低いことを表す粒子重み値がより小さい粒子が無視され、楽譜位置の推定精度が向上する。拍間隔算出部１４３は、算出した推定拍間隔＜ｂ_ｋ＞を表す推定拍間隔情報を拍間隔出力部１４６に出力する。

信頼度判定部１４４は、楽譜位置算出部１４２から入力された一部の粒子についての粒子重み値の総和、全部の粒子についての粒子重み値の総和に基づいて信頼度係数ｖ_ｋを算出する。信頼度係数ｖ_ｋは、例えば、式（１０）に示されるように、一部の粒子についての粒子重み値の総和を全部の粒子についての粒子重み値の総和で除算した値である。

即ち、信頼度係数ｖ_ｋは、粒子重み値ｗ_ｉ，ｋが大きい粒子が一部の粒子に集中している度合いであり、推定した楽譜位置や拍間隔の信頼性を表す指標値である。
信頼度判定部１４４は、現在のフィルタリングステップｋにおける信頼度係数ｖ_ｋから前フィルタリングステップｋ−１における信頼度係数ｖ_ｋ−１の差分ｖ_ｋ−ｖ_ｋ−１を算出する。信頼度判定部１４４は、差分ｖ_ｋ−ｖ_ｋ−１が予め定めた閾値−η_ｄｅｃ（例えば、−０．０７）よりも小さい場合、楽譜位置算出部１４２から入力された推定楽譜位置情報を楽譜位置出力部１４５へ出力することを停止する。この場合、信頼度係数ｖ_ｋ−１が急激に低下するため、信頼性に欠ける楽譜位置情報の出力が回避される。
信頼度判定部１４４は、差分ｖ_ｋ−ｖ_ｋ−１が予め定めた閾値η_ｉｎｃ（例えば、０．０８）よりも大きい場合、楽譜位置算出部１４２から入力された推定楽譜位置情報を楽譜位置出力部１４５へ出力することを再開する。

なお、信頼度判定部１４４は、上記の信頼度係数の代わりに推定楽譜位置情報が表す推定楽譜位置＜ｐ_ｋ＞の推定誤差ｅ（ｔ）の絶対値｜ｅ（ｔ）｜もしくは二乗誤差ｅ^２（ｔ）を信頼度係数として算出し、算出した信頼度が予め定めた閾値γ_ｅを超えるか否かを判断してもよい。信頼度判定部１４４は、信頼度係数が閾値γ_ｅを超えると判断したとき、楽譜位置算出部１４２から入力された推定楽譜位置情報の楽譜位置出力部１４５への出力を停止する。信頼度判定部１４４は、信頼度係数が閾値γ_ｅを超えないと判断したとき、楽譜位置算出部１４２から入力された推定楽譜位置情報を楽譜位置出力部１４５へ出力する。 信頼度判定部１４４は、例えば式（１１）を用いて推定誤差ｅ（ｔ）を算出する。

式（１１）において、ｓ（＜ｐ_ｋ＞）は、推定楽譜位置＜ｐ_ｋ＞の楽譜フレームに対応する時刻を表す。

楽譜位置出力部１４５は、信頼度判定部１４４から入力された推定楽譜位置情報を楽譜位置推定装置１の外部に出力する。
拍間隔出力部１４６は、拍間隔算出部１４３から入力された推定拍間隔情報を楽譜位置推定装置１の外部に出力する。拍間隔出力部１４６は、推定拍間隔情報自体の代わりに推定拍間隔情報が表す推定拍間隔＜ｂ_ｋ＞の逆数である推定テンポを出力するようにしてもよい。

次に、本実施形態に係る楽譜位置推定処理について説明する。
図３は、本実施形態に係る楽譜位置推定処理を表すフローチャートである。
（ステップＳ１０１）楽譜位置推定装置１の各構成部は、処理に用いる変数やデータを初期設定する。例えば、粒子遷移部１３２は、粒子数Ｉ、粒子ｉ毎の拍間隔ｂ^ｉ _ｋ、楽譜位置ｐ^ｉ _ｋ、粒子分布値ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）、状態遷移確率ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜ｐ^ｉ _ｋ−１，ｂ^ｉ _ｋ−１）の初期値を定める。その後、ステップＳ１０２に進む。

（ステップＳ１０２）楽譜情報入力部１１２は、楽譜情報記憶部１１１から、処理の対象とする楽曲の楽譜情報を読み出す。楽譜情報入力部１１２は、読み出した楽譜情報を関連度算出部１２１の周波数特性重み算出部１２２、拍間隔重み算出部１２３、及び楽譜位置探索部１３１の粒子遷移部１３２に出力する。その後、ステップＳ１０３に進む。

（ステップＳ１０３）音響信号入力部１０１は、音波を表すアナログ音響信号をディジタル音響信号に変換し、変換したディジタル音響信号を周波数特性分析部１０３に出力する。その後、ステップＳ１０４に進む。

（ステップＳ１０４）周波数特性分析部１０３は、音響信号入力部１０１から入力されたディジタル音響信号を時間領域信号から周波数領域信号にフレーム時刻毎に変換する。周波数特性分析部１０３は、フレーム時刻毎の周波数領域信号に基づいて音響スペクトログラムＡ_ｆ，ｔ（音響周波数特性）を算出する。周波数特性分析部１０３は、算出した音響スペクトログラムＡ_ｆ，ｔを相関算出部１０４、及び関連度算出部１２１の周波数特性重み算出部１２２に出力する。その後、ステップＳ１０５に進む。

（ステップＳ１０５）相関算出部１０４は、周波数特性分析部１０３から入力された音響スペクトログラムＡ_ｆ，ｔに基づいて、拍間隔ｂについての自己相関を、例えば式（１）を用いて算出する。相関算出部１０４は、算出した自己相関を関連度算出部１２１の拍間隔重み算出部１２３及び楽譜位置探索部１３１の粒子遷移部１３２に出力する。相関算出部１０４は、音響スペクトログラムＡ_ｆ，ｔを粒子遷移部１３２に出力する。その後、ステップＳ１０６に進む。

（ステップＳ１０６）周波数特性重み算出部１２２は、周波数特性分析部１０３から入力された音響スペクトログラムと楽譜情報入力部１１２から入力された楽譜情報に含まれる基本周波数ベクトルが表す周波数特性との楽譜位置毎の関連度、即ち類似度を表す周波数特性重み情報を生成する。周波数特性重み算出部１２２は、生成した周波数特性重み情報を粒子重み算出部１２４に出力する。その後、ステップＳ１０７に進む。

（ステップＳ１０７）拍間隔重み算出部１２３は、相関算出部１０４から入力された自己相関に基づいて拍間隔毎に振幅特性の類似度を表す拍間隔重み情報を生成する。拍間隔重み算出部１２３は、生成した拍間隔重み情報を粒子重み算出部１２４に出力する。その後、ステップＳ１０８に進む。

（ステップＳ１０８）粒子重み算出部１２４は、周波数特性重み算出部１２２から入力された周波数特性重み情報、拍間隔重み算出部１２３から入力された拍間隔重み情報、粒子遷移部１３２から入力された粒子ｉ毎の粒子分布値、状態遷移確率、粒子情報に基づいて粒子重み情報を生成する。粒子重み算出部１２４は、例えば式（９）を用いて粒子重み値ｗ_ｉ，ｋを算出し、算出した粒子重み値を表す粒子重み情報を生成する。粒子重み算出部１２４は、入力された粒子ｉ毎の粒子情報に含まれる粒子重み情報を、生成した粒子重み情報に置き換えることで粒子情報を更新する。
粒子重み算出部１２４は、更新した粒子情報を再標本化部１３３、楽譜位置算出部１４２及び拍間隔算出部１４３に出力する。その後、ステップＳ１０９に進む。

（ステップＳ１０９）楽譜位置算出部１４２は、粒子重み算出部１２４から入力された粒子ｉ毎の粒子情報から楽譜位置情報と粒子重み情報を抽出する。楽譜位置算出部１４２は、粒子ｉ毎に抽出した楽譜位置情報が表す楽譜位置ｐ_ｉ ^ｋに基づいて推定楽譜位置＜ｐ_ｋ＞を算出する。楽譜位置算出部１４２は、例えば、粒子重み情報が表す粒子重み値ｗ_ｉ，ｋによる加重平均を推定楽譜位置＜ｐ_ｋ＞とする。楽譜位置算出部１４２は、算出した推定楽譜位置＜ｐ_ｋ＞を表す推定楽譜位置情報を信頼度判定部１４４に出力する。楽譜位置算出部１４２は、一部の粒子についての粒子重み値の総和と、全部の粒子についての粒子重み値の総和を算出する。楽譜位置算出部１４２は、算出した一部の粒子についての粒子重み値の総和、全部の粒子についての粒子重み値の総和を信頼度判定部１４４に出力する。その後、ステップＳ１１０に進む。

（ステップＳ１１０）拍間隔算出部１４３は、粒子重み算出部１２４から入力された粒子ｉ毎の粒子情報から拍間隔情報と粒子重み情報を抽出する。拍間隔算出部１４３は、例えば、粒子重み情報が表す粒子重み値ｗ_ｉ，ｋによる加重平均を推定拍間隔＜ｂ_ｋ＞として算出する。拍間隔算出部１４３は、算出した推定拍間隔＜ｂ_ｋ＞を表す推定拍間隔情報を拍間隔出力部１４６に出力する。その後、ステップＳ１１１に進む。

（ステップＳ１１１）信頼度判定部１４４は、楽譜位置算出部１４２から入力された一部の粒子についての粒子重み値の総和、全部の粒子についての粒子重み値の総和に基づいて信頼度係数ｖ_ｋを、例えば式（１０）を用いて算出する。
信頼度判定部１４４は、現在のフィルタリングステップｋにおける信頼度係数ｖ_ｋから前フィルタリングステップｋ−１における信頼度係数ｖ_ｋ−１の差分ｖ_ｋ−ｖ_ｋ−１を算出する。その後、ステップＳ１１２に進む。

（ステップＳ１１２）信頼度判定部１４４は、差分ｖ_ｋ−ｖ_ｋ−１が予め定めた閾値−η_ｄｅｃ（閾値１）よりも小さいか否かを判断する。予め定めた閾値−η_ｄｅｃよりも小さいと判断したとき（ステップＳ１１２ Ｙ）、ステップＳ１１３に進む。予め定めた閾値−η_ｄｅｃよりも小さくないと判断したとき（ステップＳ１１２ Ｎ）、ステップＳ１１４に進む。
（ステップＳ１１３）信頼度判定部１４４は、楽譜位置算出部１４２から入力された推定楽譜位置情報＜ｐ_ｋ＞の楽譜位置出力部１４５への出力を停止する。その後、ステップＳ１１４に進む。

（ステップＳ１１４）信頼度判定部１４４は、差分ｖ_ｋ−ｖ_ｋ−１が予め定めた閾値η_ｉｎｃ（閾値２）よりも大きいか否かを判断する。予め定めた閾値η_ｉｎｃよりも大きいと判断したとき（ステップＳ１１４ Ｙ）、ステップＳ１１５に進む。予め定めた閾値η_ｉｎｃよりも大きくないと判断したとき（ステップＳ１１４ Ｎ）、ステップＳ１１６に進む。
（ステップＳ１１５）信頼度判定部１４４は、楽譜位置算出部１４２から入力された推定楽譜位置情報＜ｐ_ｋ＞の楽譜位置出力部１４５への出力を再開する。その後、ステップＳ１１６に進む。

（ステップＳ１１６）再標本化部１３３は、入力された粒子情報に含まれる粒子重み情報が表す粒子重みに比例した確率で粒子数Ｉよりも多い予め定めた数の粒子を抽出する。再標本化部１３３は、増加した粒子数と同一の数だけ、粒子重み値ｗ_ｉ，ｋがより小さい粒子を優先して抽出した粒子を棄却する。再標本化部１３３は、抽出された各粒子の楽譜位置情報及び拍間隔情報は、対応する抽出前の粒子の楽譜位置情報及び拍間隔情報と同一とする。再標本化部１３３は、棄却されずに残った粒子の粒子情報を粒子遷移部１３２に出力する。その後、ステップＳ１１７に進む。

（ステップＳ１１７）粒子遷移部１３２は、相関算出部１０４から距離値及び自己相関が、楽譜情報入力部１１２から入力された楽譜情報が、再標本化部１３３から粒子毎の粒子情報が入力される。粒子遷移部１３２は、入力された楽譜情報に基づいてオンセット情報を生成する。
粒子遷移部１３２は、楽譜情報に含まれるテンポ情報、自己相関、粒子情報に基づいて拍間隔分布値を算出する。粒子遷移部１３２は、算出した拍間隔分布値に比例する確率で粒子ｉの拍間隔ｂ^ｉ _ｋを定める。粒子遷移部１３２は、オンセット情報、距離値、及び抽出した拍間隔ｂ^ｉ _ｋに基づいて楽譜位置分布値を算出し、算出した楽譜位置分布値に比例する確率で粒子ｉの楽譜位置ｐ^ｉ _ｋを定める。
粒子遷移部１３２は、粒子ｉ毎の拍間隔分布値に、楽譜位置分布値を乗じて、粒子分布値ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）を算出する。

粒子遷移部１３２は、粒子ｉ毎に前フィルタリングステップｋ−１における楽譜位置ｐ^ｉ _ｋ−１、拍間隔ｂ^ｉ _ｋ−１と、現フィルタリングステップｋで定めた楽譜位置ｐ^ｉ _ｋ、拍間隔ｂ^ｉ _ｋに基づいて、状態遷移確率ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜ｐ^ｉ _ｋ−１，ｂ^ｉ _ｋ−１）を算出する。
粒子遷移部１３２は、粒子毎に粒子情報に含まれる拍間隔情報及び楽譜位置情報を、抽出した拍間隔を表す拍間隔情報及び楽譜位置を表す楽譜位置情報にそれぞれ置き換える。粒子遷移部１３２は、置き換えた粒子情報を関連度算出部１２１の周波数特性重み算出部１２２、拍間隔重み算出部１２３、粒子重み算出部１２４に出力する。
粒子遷移部１３２は、算出した粒子毎の粒子分布値と状態遷移確率を粒子重み算出部１２４に出力する。その後、処理を終了する。

ここで粒子遷移部１３２の構成について、より詳細に説明する。図１に戻り、粒子遷移部１３２は、相関算出部１０４から距離値及び自己相関が、楽譜情報入力部１１２から入力された楽譜情報が、再標本化部１３３から粒子毎の粒子情報が入力される。粒子遷移部１３２は、入力された楽譜情報からコードの変化（オンセット）を有するか否かを表すオンセット情報を生成する。コードの変化は、コードに含まれる音階の変化である。従って、コードの変化が楽譜位置では音響信号は、変化直後の音階の周波数成分が新たに発生し、振幅の立ち上がり（オンセット）として観測される。ここで、粒子遷移部１３２は、基本周波数ベクトル｛μ_ｐ｝が、直前の楽譜フレームｐ−１における基本周波数ベクトル｛μ_ｐ−１｝と異なる楽譜フレームｐがオンセットを有すると判断する。粒子遷移部１３２は、基本周波数ベクトル｛μ_ｐ｝が、直前の楽譜フレームｐ−１における基本周波数ベクトル｛μ_ｐ−１｝と同一ではない楽譜フレームｐはオンセットを有しないと判断する。粒子遷移部１３２は、楽譜フレームｐ毎にオンセットを有するか否かを表すオンセット情報を生成する。例えば、楽譜フレームｐがオンセットを有することを表すオンセット情報ｏ_ｐは１であり、楽譜フレームｐがオンセットを有しないことを表すオンセット情報ｏ_ｐは０である。

なお、本実施形態では、粒子遷移部１３２は、楽譜情報に対応するオンセット情報を楽譜情報入力部１１２から入力されるようにしてもよい。その場合、予め生成したオンセット情報を楽譜情報と対応付けて楽譜情報記憶部１１１に記憶させておき、楽譜情報入力部１１２が楽譜情報記憶部１１１から読み出したオンセット情報を粒子遷移部１３２に出力する。

粒子遷移部１３２は、生成したオンセット情報、入力された距離値、自己相関、及び楽譜情報に含まれるテンポ情報に基づいて粒子分布値を算出する。粒子遷移部１３２は、算出した粒子分布値が表す粒子分布で粒子を抽出することで粒子を遷移させる。粒子遷移部１３２は、遷移させた粒子毎の粒子情報を関連度算出部１２１の周波数特性重み算出部１２２及び拍間隔重み算出部１２４に出力する。

ここで、粒子遷移部１３２が、粒子分布値を算出する処理について説明する。
粒子遷移部１３２は、テンポ情報に基づいて拍間隔ｂに対する窓関数ψ（ｂ｜ｂ^ｓ’）を算出する。窓関数ψ（ｂ｜ｂ^ｓ’）は、例えば、テンポ情報が表すテンポ６０／ｂ^ｓ’を中心とする予め定めた範囲内｜６０／ｂ−６０／ｂ^ｓ’｜＜６０／ｂ^θでは１であり、その範囲外では０である。以下、６０／ｂ^θをテンポ窓長（ｔｅｍｐｏ ｗｉｎｄｏｗ ｌｅｎｇｔｈ）と呼ぶ。

粒子遷移部１３２は、自己相関Ｒ（ｂ，｛Ａ_ｋ｝）に窓関数ψ（ｂ｜ｂ^ｓ’）を乗じて拍間隔ｂに対する分布値ｑ（ｂ｜｛Ａ_ｋ｝，ｂ^ｓ’）を算出する。これにより、粒子遷移部１３２は、楽曲の拍間隔としてあり得ない拍間隔を排除し、拍構造の周期性が高い拍間隔をもつ粒子が抽出される確率を高くする。そして、粒子遷移部１３２は、入力された粒子情報から、算出した分布値ｑ（ｂ｜｛Ａ_ｋ｝，ｂ^ｓ’）に比例する確率で粒子ｉの拍間隔ｂ^ｉ _ｋを定める。
なお、粒子遷移部１３２は、自己相関Ｒ（ｂ，｛Ξ_ｋ’｝）が相関算出部１０４から入力された場合には、Ｒ（ｂ，｛Ａ_ｋ｝）の代わりにＲ（ｂ，｛Ξ_ｋ’｝）を用いて分布値ｑ（ｂ｜｛Ａ_ｋ｝，ｂ^ｓ’）を算出してもよい。

粒子遷移部１３２は、生成したオンセット情報、距離値、及び抽出した拍間隔ｂ^ｉ _ｋに基づいて楽譜位置ｐに対する分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）を算出する。
ここで、粒子遷移部１３２は、例えば式（１２）を用いて分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）を算出する。

式（１２）において、τ’（ｐ，ｂ）は、楽譜位置の探索範囲Ｐ^ｉ _ｋに含まれ、かつ時刻ｋΔＴにおける音響フレームと対応付けられた楽譜位置ｐにおけるオンセット情報ｏ_ｐを表す。探索範囲Ｐ^ｉ _ｋは、前フィルタリングステップｋ−１における粒子ｉの楽譜位置ｐ^ｉ _ｋ−１から拍間隔ｂ＝ｂ^ｉ _ｋだけ進んだ楽譜位置ｐ^ｉ _ｋ−１’＝ｐ^ｉ _ｋ−１＋ΔＴ／ｂ^ｉ _ｋ−１を中心とし、予め定めた幅２δ_ｐ（例えば、δ_ｐ＝３）を有する。
即ち、式（１２）は、探索範囲Ｐ^ｉ _ｋにオンセットが含まれる場合、分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）は、オンセットを有する楽譜位置に対応する時刻における音響信号の振幅の累積値に比例する値であることを表す。探索範囲Ｐ^ｉ _ｋにオンセットが含まれない場合、１に比例する値であることを表す。楽譜位置ｐが、探索範囲Ｐ^ｉ _ｋの範囲外である場合には、分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）は、ゼロであることを表す。

ここで、楽譜位置と音響信号の振幅との関係について説明する。
図４は、楽譜情報と音響信号の関係の一例を表す概念図である。
図４の上段は、楽譜情報を表す。楽譜情報は、左から順に、１拍の音階Ａ４からなるコードｄ１、０．５拍のコードＢ４、Ｆ３からなるコードｄ２、０．５拍の音階Ｂ４、Ｇ３からなるコードｄ３を表す。この例では、音階Ｂ３を２つの０．５拍からなる音として扱う。図４の上下に延びる３本の破線は、各々オンセットを表す。図４では、オンセットを有する楽譜位置は、左端、左端から１拍、左端から１．５拍である。なお、各コードが継続する区間をセグメントといい、ｄ１、ｄ２、ｄ３は、各コードのセグメントを表す。

図４の下段は音響信号のスペクトログラムを表す画像である。水平方向は時刻を表し、上下方向は周波数を表す。この画像において濃淡は振幅を表す。明るい部分ほど振幅が大きく、暗い部分ほど振幅が小さいことを表す。図４の下段は、各周波数において振幅が一斉に増加し、その後、周波数が高い成分ほど早く振幅が減衰する時間変化が３回繰り返されることを表す。振幅が一斉に増加する部分が、音響信号の振幅が急激に増加する部分である。従って、図４において明るい部分を破線が通過していることは、音響信号の振幅が急激に増加する部分と楽譜情報に基づくオンセットと合致していることを表す。このとき、分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）は、大きな値をとる。言い換えれば、音響信号の振幅が急激に増加する部分が、楽譜情報に基づくオンセットと合致していない場合、分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）はより小さい値をとる。
これにより、粒子遷移部１３２は、分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）を用いることで、楽譜情報に含まれるオンセットと音響信号の振幅が増加する部分の楽譜位置ｐ^ｉ _ｋが推定される確率を高くすることができる。

図１に戻り、粒子遷移部１３２は、自己相関Ｒ（ｂ，｛Ξ_ｋ’｝）が相関算出部１０４から入力された場合には、Ｒ（ｂ，｛Ａ_ｋ｝）の代わりにＲ（ｂ，｛Ξ_ｋ’｝）を用いてもよい。この場合において、粒子遷移部１３２は、例えば式（１３）を用いて分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）を算出してもよい。

式（１３）において、ζ_ｔは、周波数毎の距離値Ξ_ｆ，ｔを累積した累積距離値ζ_ｔを表す。例えば、累積距離値ζ_ｔは、式（１４）を用いて算出される。

式（１３）に戻り、ｐ’（τ）は、拍間隔ｂ^ｉ _ｋについての時刻τに対応する楽譜位置を表す。即ち、ｐ’（τ）＝ｐ−（ｔ−τ）／ｂ^ｉ _ｋという関係がある。この関係において、時刻ｔにおける楽譜位置がｐであると仮定されている。Ｐは、楽譜位置の探索範囲を表す。楽譜範囲Ｐは、前フィルタリングステップｋ−１における楽譜位置ｐ^ｉ _ｋ−１から拍間隔ｂ^ｉ _ｋだけ進んだ楽譜位置ｐ^ｉ _ｋ−１＋ΔＴ／ｂ^ｉ _ｋを中心とし、前後に予め定めた幅３δ_ｐ（例えば、δ_ｐ＝１）を有する。この楽譜位置は、粒子ｉの現フィルタリングステップｋにおける推定楽譜位置である。即ち、式（１３）は、探索範囲Ｐにオンセットが含まれる場合、分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）は、オンセットを有する楽譜位置に対応する時刻における累積距離値ζ_ｔに比例する値であることを表す。探索範囲Ｐにオンセットが含まれない場合、１に比例する値であることを表す。楽譜位置ｐが、探索範囲Ｐの外である場合には、分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）は、ゼロであることを表す。

粒子遷移部１３２は、入力された粒子情報から算出した分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）に比例する確率で粒子ｉの楽譜位置ｐ^ｉ _ｋを定める。
なお、上述の探索範囲Ｐ^ｉ _ｋ−１（又はＰ）にオンセットが含まれない場合には、粒子遷移部１３２は、当該探索範囲からランダムに粒子ｉの楽譜位置ｐ^ｉ _ｋを抽出する。
粒子遷移部１３２は、粒子ｉ毎に算出した拍間隔の分布値（拍間隔分布値）ｑ（ｂ^ｉ _ｋ｜｛Ａ_ｋ｝，ｂ^ｓ’）に、楽譜位置の分布値（楽譜位置分布値）ｑ（ｐ^ｉ _ｋ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）を乗じて、粒子の分布値（粒子分布値）ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）を算出する。

粒子遷移部１３２は、粒子ｉ毎にフィルタリングステップｋ−１からｋへの状態遷移確率ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜ｐ^ｉ _ｋ−１，ｂ^ｉ _ｋ−１）を、例えば式（１５）を用いて算出する。

式（１５）において、Ｎ（ａ｜ｂ，σ^２）は、変数ａについての正規分布関数（ガウス関数、ガウシアンとも呼ばれる）を表す。ここで、原点がｂであり、分散がσ^２である。σ^２ _ｂは、拍間隔ｂの分散を表す実数値（例えば、０．２）である。ｐ^ｉ _ｋ−１＋ΔＴ／ｂ^ｉ _ｋ−１は、前フィルタリングステップｋ−１における楽譜位置ｐ^ｉ _ｋ−１から拍間隔ｂ^ｉ _ｋ−１だけ進んだ楽譜位置を表す。σ^２ _ｐは、楽譜位置ｐの分散を表す実数値（例えば、１）である。なお、楽譜位置ｐ^ｉ _ｋ−１は、粒子ｉの楽譜位置情報が表す値であり、拍間隔ｂ^ｉ _ｋ−１は、粒子ｉの拍間隔情報が表す値である。これにより、実際の演奏によって生じる拍間隔や楽譜位置の揺らぎが考慮される。
粒子遷移部１３２は、粒子毎に粒子情報に含まれる拍間隔情報及び楽譜位置情報を、抽出した拍間隔を表す拍間隔情報及び楽譜位置を表す楽譜位置情報にそれぞれ置き換える。粒子遷移部１３２は、置き換えた粒子情報を関連度算出部１２１の周波数特性重み算出部１２２、拍間隔重み算出部１２３、粒子重み算出部１２４に出力する。
粒子遷移部１３２は、算出した粒子毎の粒子分布値ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）と状態遷移確率ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜ｐ^ｉ _ｋ−１，ｂ^ｉ _ｋ−１）を粒子重み算出部１２４に出力する。

初期（ｋ＝０）において、粒子遷移部１３２は、例えば、テンポ６０／ｂ^ｓ’を中心とする予め定めた範囲内（ｂ^ｓ’−６０／ｂ^θ＜ｂ＜ｂ^ｓ’＋６０／ｂ^θ）の一様分布から粒子ｉ毎の拍間隔ｂ^ｉ _ｋを定める。また、粒子遷移部１３２は、粒子ｉ毎の楽譜位置ｐ^ｉ _ｋをゼロと定める。粒子遷移部１３２は、粒子分布値ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）と状態遷移確率ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜ｐ^ｉ _ｋ−１，ｂ^ｉ _ｋ−１）を１と定める。

次に、本実施形態に係る粒子遷移部１３２が行う粒子遷移処理について説明する。
図５は、本実施形態に係る粒子遷移処理を表すフローチャートである。
（ステップＳ２０１）粒子遷移部１３２は、相関算出部１０４から距離値及び自己相関が、楽譜情報入力部１１２から入力された楽譜情報が、再標本化部１３３から粒子毎の粒子情報が入力される。粒子遷移部１３２は、入力された楽譜情報に基づいてオンセット情報を生成する。その後、ステップＳ２０２に進む。

（ステップＳ２０２）粒子遷移部１３２は、入力された楽譜情報に含まれるテンポ情報に基づいて拍間隔ｂに対する窓関数ψ（ｂ｜ｂ^ｓ’）を算出する。粒子遷移部１３２は、入力された自己相関Ｒ（ｂ，｛Ａ_ｋ｝）に窓関数ψ（ｂ｜ｂ^ｓ’）を乗じて拍間隔分布値ｑ（ｂ｜｛Ａ_ｋ｝，ｂ^ｓ’）を算出する。その後、ステップＳ２０３に進む。
（ステップＳ２０３）粒子遷移部１３２は、入力された粒子情報から、算出した拍間隔分布値ｑ（ｂ｜｛Ａ_ｋ｝，ｂ^ｓ’）に比例する確率で粒子ｉの拍間隔ｂ^ｉ _ｋを抽出する。その後、ステップＳ２０４に進む。

（ステップＳ２０４）粒子遷移部１３２は、生成したオンセット情報、距離値、及び抽出した拍間隔ｂ^ｉ _ｋに基づいて楽譜位置分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）を算出する。粒子遷移部１３２は、例えば式（１２）を用いて楽譜位置分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）を算出する。その後、ステップＳ２０５に進む。

（ステップＳ２０５）粒子遷移部１３２は、粒子遷移部１３２は、楽譜位置分布値ｑ（ｐ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）に比例する確率で粒子ｉの楽譜位置ｐ^ｉ _ｋを抽出する。但し、探索範囲Ｐ^ｉ _ｋ−１（又はＰ）にオンセットが含まれない場合には、粒子遷移部１３２は、探索範囲からランダムに粒子ｉの楽譜位置ｐ^ｉ _ｋを抽出する。その後、ステップＳ２０６に進む。

（ステップＳ２０６）粒子遷移部１３２は、粒子遷移部１３２は、粒子ｉ毎に拍間隔分布値ｑ（ｂ^ｉ _ｋ｜｛Ａ_ｋ｝，ｂ^ｓ’）に、楽譜位置分布値ｑ（ｐ^ｉ _ｋ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）を乗じて、粒子分布値ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）を算出する。その後、ステップＳ２０７に進む。

（ステップＳ２０７）粒子遷移部１３２は、粒子ｉ毎に前フィルタリングステップｋ−１における楽譜位置ｐ^ｉ _ｋ−１、拍間隔ｂ^ｉ _ｋ−１、抽出された楽譜位置ｐ^ｉ _ｋ−１、拍間隔ｂ^ｉ _ｋ−１に基づいて、状態遷移確率ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜ｐ^ｉ _ｋ−１，ｂ^ｉ _ｋ−１）を算出する。粒子遷移部１３２は、状態遷移確率ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜ｐ^ｉ _ｋ−１，ｂ^ｉ _ｋ−１）を算出する際、例えば式（１５）を用いる。その後、ステップＳ２０８に進む。

（ステップＳ２０８）粒子遷移部１３２は、粒子毎に粒子情報に含まれる拍間隔情報及び楽譜位置情報を、抽出した拍間隔を表す拍間隔情報及び楽譜位置を表す楽譜位置情報にそれぞれ置き換える。粒子遷移部１３２は、置き換えた粒子情報を関連度算出部１２１の周波数特性重み算出部１２２、拍間隔重み算出部１２３、粒子重み算出部１２４に出力する。
粒子遷移部１３２は、算出した粒子毎の粒子分布値ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜｛Ａ_ｋ｝，ｏ_ｐ，ｂ^ｉ _ｋ）と状態遷移確率ｑ（ｐ^ｉ _ｋ，ｂ^ｉ _ｋ｜ｐ^ｉ _ｋ−１，ｂ^ｉ _ｋ−１）を粒子重み算出部１２４に出力する。その後、処理を終了する。

図１に戻り、再標本化部１３３の動作について、より詳細に説明する。再標本化部１３３は、粒子重み算出部１２４から入力された粒子情報に基づいて各粒子を再標本化する。
再標本化部１３３は、入力された粒子情報に含まれる粒子重み情報が表す粒子重みに比例した確率で粒子を抽出する。再標本化部１３３は、粒子を抽出する際、例えばＳＩＲ（ｓａｍｐｌｉｎｇ ｉｍｐｏｒｔａｎｃｅ ｒｅｓａｍｐｌｉｎｇ）法を用いる。
再標本化部１３３は、粒子ｉ毎の粒子重み係数を、次式（１６）を用いて規格化し、規格化粒子重み係数ｗ_ｉ，ｋ’を算出する。

式（１６）において、Ｉは、予め定めた粒子数を示す整数（例えば、１５００）を表す。
再標本化部１３３は、算出した確率（規格化粒子重み係数ｗ_ｉ，ｋ’）で次回の粒子ｉを抽出する。ここで、再標本化部１３３は、抽出された各粒子の粒子重みは、各々等しい値（例えば、１／Ｉ）とし、抽出された各粒子の楽譜位置情報及び拍間隔情報は、対応する抽出前の粒子の楽譜位置情報及び拍間隔情報と同一とする。再標本化部１３３は、抽出した粒子毎の楽譜位置情報、拍間隔情報及び粒子重みを含む粒子情報を生成する。
再標本化部１３３は、これにより、粒子数Ｉよりも多い予め定めた数の粒子を抽出する。その後、再標本化部１３３は、増加した粒子数と同一の数だけ、粒子を棄却する。ここで、再標本化部１３３は、抽出前の粒子重みが小さい粒子ほど優先して、その粒子に基づいて抽出した粒子情報を棄却する。これにより、再標本化部１３３は、粒子数（粒子情報の数）を一定数Ｉに保つ。再標本化部１３３は、棄却されずに残った粒子の粒子情報を粒子遷移部１３２に出力する。

次に、拍間隔重み算出部１２３の動作について説明する。
拍間隔重み算出部１２３は、粒子遷移部１３２から入力された粒子情報に含まれる拍間隔情報を粒子ｉ毎に抽出する。拍間隔重み算出部１２３は、相関算出部１０４から入力された自己相関Ｒ（ｂ，｛Ａ_ｋ｝）又はＲ（ｂ，｛Ξ_ｋ’｝）から抽出した拍間隔情報が表す拍間隔ｂ_ｉ ^ｋに対応する値Ｒ（ｂ^ｉ _ｋ，｛Ａ_ｋ｝）又はＲ（ｂ^ｉ _ｋ，｛Ξ_ｋ’｝）を拍間隔重み値ｗ^ｔｍ _ｉ，ｋと定める。拍間隔重み算出部１２３は、粒子ｉ毎の拍間隔重み値ｗ^ｔｍ _ｉ，ｋを表す拍間隔重み情報を粒子重み算出部１２４に出力する。

なお、本実施形態では、拍間隔重み算出部１２３を省略してもよい。拍間隔重み算出部１２３を省略する場合、粒子遷移部１３２が相関算出部１０４から入力された自己相関Ｒ（ｂ，｛Ａ_ｋ｝）又はＲ（ｂ，｛Ξ_ｋ’｝）から拍間隔ｂ^ｉ _ｋに対応する値Ｒ（ｂ^ｉ _ｋ，｛Ａ_ｋ｝）又はＲ（ｂ^ｉ _ｋ，｛Ξ_ｋ’｝）を拍間隔重み値ｗ^ｔｍ _ｉ，ｋと定める。粒子遷移部１３２は、粒子ｉ毎に定めた拍間隔重み値ｗ^ｔｍ _ｉ，ｋを表す拍間隔重み情報を粒子重み算出部１２４に出力してもよい。

次に、周波数特性重み算出部１２２について説明する。
本実施形態では、楽譜情報に基づく振幅周波数特性（例えば、楽譜スペクトログラム）が、既知の楽譜情報に含まれる楽譜位置毎の音階の基本周波数に基づく調波構造で表される成分を少なくとも１つ含むと仮定する。そこで、周波数特性重み算出部１２２は、観測値である音響信号に基づく振幅周波数特性（音響スペクトログラム）に対して尤度が最大となるように、調波構造毎の楽譜スペクトログラムへの寄与を表す変数の確率分布や、各調波構造に含まれる調波成分毎の当該調波構造への寄与を表す変数を変分パラメータとして算出する。周波数特性重み算出部１２２は、算出した変数が表す楽譜スペクトログラムが観測値として音響スペクトログラムが与えられているときの確率を算出し、算出した確率に基づいて周波数特性重み値ｗ^ｓｐ _ｉ，ｋを算出する。本実施形態では、上述の尤度が最大になるような周波数特性重み値ｗ^ｓｐ _ｉ，ｋを算出できればよいので、楽譜情報に基づく振幅周波数特性（楽譜スペクトログラム）を陽に合成しなくともよい。

周波数特性重み算出部１２２は、例えば、ＬＨＡ（Ｌａｔｅｎｔ Ｈａｒｍｏｎｉｃ Ａｌｌｏｃａｔｉｏｎ；潜在的調波配分）法を用いて周波数特性分析部１０３から入力された音響スペクトログラムと楽譜情報入力部１１２から入力された楽譜情報に基づいて周波数特性重み情報を生成する。
ＬＨＡ法は、入力された音響信号のスペクトログラムを観測値とし、この観測値を調波構造と成分として含み、調波構造の楽譜スペクトログラムに対する寄与を表す変数の確率分布、調波成分の調波構造への寄与を表す係数の確率分布として表す音源モデルである。調波構造は、基本周波数の成分と基本周波数の整数倍の周波数の成分（高調波）を調波成分として含む周波数特性を表す。従って、基本周波数、各次数の振幅（もしくは、基準値に対する比率）が、調波構造を表す変数に含まれる。

ここで、ＬＨＡ法に基づく音源モデルについて説明する。
図６は、本実施形態に係る音源モデルの一例を表す概念図である。
図６の右上に示されるｘ_ｎは、周波数を表す。ｘ_ｎを囲う二重円は、ｘ_ｎが観測変数であることを表す。ｚ_ｎからｘ_ｎに向かう矢印、μ_ｌからｘ_ｎに向かう矢印及びΛ_ｌからｘ_ｎに向かう矢印は、周波数ｘ_ｎの潜在変数（ｌａｔｅｎｔ ｖａｒｉａｂｌｅ）ｚ_ｎ、基本周波数μ_ｌ及び精度（ｐｒｅｃｉｓｉｏｎ）Λ_ｌに対する尤度ｑ（Ｘ｜Ｚ，｛μ｝，｛Λ｝）が仮定されていることを表す。Ｘ、Ｚ、｛μ｝、｛Λ｝は、それぞれ変数ｘ_ｎ、ｚ_ｎ、μ_ｌ、Λ_ｌの集合（セット）を表す。ｚ_ｎ、Λ_ｌを囲う円は、ｚ_ｎ、Λ_ｌがそれぞれ未知変数であることを示す。μ_ｌが円で囲まれていないことは、基本周波数μ_ｌが既知変数であることを表す。基本周波数μ_ｌは、楽譜情報の基本周波数ベクトルに含まれる各音階の基本周波数であるためである。
π_ｄｌからｚ_ｎに向かう矢印、θ_ｌｍからｚ_ｎに向かう矢印、潜在変数ｚ_ｎの確率π_ｄｌ及び確率θ_ｌｍに対する尤度ｑ（Ｚ｜｛π｝，｛θ｝）が仮定されていることを表す。ここで、｛π｝，｛θ｝は、それぞれ変数π_ｄｌ、θ_ｌｍの集合（セット）を表す。π_ｄｌ、θ_ｌｍを囲う円は、π_ｄｌ、θ_ｌｍがそれぞれ未知変数であることを示す。末端に示されている未知変数π_ｄｌ、θ_ｌｍ、Λ_ｌに対して、それぞれ事前分布ｑ（｛π｝）、ｑ（｛θ｝）、ｑ（｛Λ｝）が仮定されていることを表す。従って、ＬＨＡ法では、確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝）は、ｑ（Ｘ｜Ｚ，｛μ｝，｛Λ｝）、ｑ（Ｘ｜Ｚ，｛μ｝，｛Λ｝）、ｑ（｛π｝）、ｑ（｛θ｝）、ｑ（｛Λ｝）の積で表される。

図６において、変数を囲う四角形は、それぞれ尤度又は事前分布を与える単位、即ち、観測値を表現する確率分布を表現するための単位を表す。ｘ_ｎとｚ_ｎを囲う四角形は、ｘ_ｎとｚ_ｎが楽譜フレームｎ毎に与えられることを表す。この四角形の右下に表されているＮ_ｄは、セグメントｄに含まれるフレーム数を表す。ｚ_ｎ、ｘ_ｎとともにπ_ｄｌを囲う四角形は、π_ｄｌがセグメントｄ毎に与えられることを表す。この四角形の右下に表されているＤは、全セグメント数を表す。θ_ｌｍとともにμ_ｌ、Λ_ｌを囲う四角形は、μ_ｌ、Λ_ｌが調波構造ｌ毎に与えられることを表す。この四角形の右下に表されているＬ_ｄは、調波構造の総数を表す。θ_ｌｍを囲う四角形は、θ_ｌｍが調波構造ｌを構成する次数ｍの調波毎に与えられることを表す。但し、ｍ＝１は、基本周波数を表す。この四角形の右下に示されているＭは、各調波構造に含まれる調波の個数（最高次数）を表す。

ＬＨＡ法では、調波構造ｌを、例えばＧＭＭを用いて表す。即ち、調波構造ｌは、次数ｍについて周波数ｍμ_ｌを平均値とするガウス関数を確率θ_ｌｍ（１≦ｍ≦Ｍ）で線形重み付け加算した確率分布で表される。ここで、調波構造ｌ毎の確率θ_ｌｍの総和Σ_ｌ ^Ｍθ_ｌｍは１と規格化されている。即ち、確率θ_ｌｍは、（Ｍは、各調波構造が含む調波の数を表す予め定められた整数である（例えば、Ｍ＝１０）。各次数ｍについて、ガウス関数の分散は精度Λ_ｌである。即ち、確率θ_ｌｍは、次数間の寄与の比率、言い換えれば次数ｍの調波成分の調波構造ｌへの寄与を表す変数である。精度Λ_ｌは、基本周波数μ_ｌの誤差の大きさを表す。

図７は、本実施形態に係る調波構造の一例を表す図である。
図７において、縦軸は確率を表し、横軸は周波数を表す。図７において、細い破線は調波構造ｌが、周波数μ_ｌ，２μ_ｌ，…，Ｍμ_ｌにおいて確率のピーク値θ_ｌ１，θ_ｌ２，…，θ_ｌＭを有することを表す。一点破線は、ピーク値θ_ｌ１，θ_ｌ２，…，θ_ｌＭを表す。周波数２μ_ｌを中心とする確率値の曲線（ガウス関数）を挟む２つの矢印の間隔は、精度Λ_ｌを表す。

ＬＨＡ法では、観測値である音響スペクトログラムを周波数ｘ_ｎの確率分布とみなす。本実施形態では、周波数ｘ_ｎの確率分布を、調波構造ｌを確率π_ｄｌ（１≦ｌ≦ｌＬ_ｄ）で線形重み付け加算した確率分布で近似する。即ち、確率π_ｄｌは、調波構造ｌの楽譜スペクトラムへの寄与を表す変数であり、この楽譜スペクトログラムは陽に算出されないがをも音響スペクトログラムを近似する。本実施形態では、音響スペクトログラム、調波構造ｌに含まれる確率値（即ち、振幅）を予め定めた量子化幅で量子化してもよい。これにより、確率分布は、振幅を度数とするヒストグラムとして表されるため演算量が低減する。

図８は、本実施形態に係るヒストグラムの一例を表す図である。
図８の左上は、調波構造１を表すヒストグラムである。図８の左下は、調波構造ｌを現すヒストグラムである。図８の右側は、周波数ｘ_ｎの確率分布を表すヒストグラムである。これらのヒストグラムにおいて、縦軸は振幅を表し、横軸は周波数を表す。調波構造１の振幅は、周波数μ_１、２μ_１、３μ_１、４μ_１、５μ_１、６μ_１においてピークを有する。各ピーク値の大きさは、楽器の音色などに応じて異なる。調波構造ｌの振幅は、周波数μ_ｌ、２μ_ｌにおいてピークを有する。図８に示す例では、周波数μ_ｌは周波数μ_１よりも高く、周波数μ_ｌにおけるピーク値は、周波数μ_１におけるピーク値よりも小さい。
調波構造１を表すヒストグラムの右側に表されているπ_ｄ１は、調波構造１に対する確率を表す。調波構造ｌを表すヒストグラムの右側に表されているπ_ｄｌは、調波構造ｌに対する確率を表す。図８の最左辺の上下に延びる括弧と「Ｌ_ｄ個」の文字は、Ｌ_ｄ個の調波構造ｌを確率π_ｄｌで重み付け加算して、図８の右側における周波数ｘ_ｎの確率分布を近似することを表す。

従って、本実施形態では観測値Ｘが与えられているときに、全ての未知変数の事後分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐ^ｉ _ｋ，｛μ｝）を粒子ｉ毎に算出する。全ての未知変数に対する事後分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐ^ｉ _ｋ，｛μ｝）を累積した累積確率値は、本音源モデルの適合性を表す。周波数特性重み算出部１２２は、例えば、式（１７）を用いて算出した累積確率値を周波数特性重み値の対数値ｌｏｇｗ^ｓｐ _ｉ，ｋとして算出する。

しかし、式（１７）を用いた演算、とりわけＺに対する累算を解析的に行うことは困難である。そこで、本実施形態では、変分ベイズ法（ｖａｒｉａｔｉｏｎａｌ Ｂａｙｅｓ［ＶＢ］ ｍｅｔｈｏｄ）を用いて、事後分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝｜Ｘ，ｐ^ｉ _ｋ，｛μ｝）を近似する。本実施形態では、解析的に求めることができる変分事後分布ｑ^＊（Ｚ，｛π｝，｛θ｝，｛Λ｝）を仮定し、真の事後分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝｜Ｘ，ｐ^ｉ _ｋ，｛μ｝）との指標値（例えば、カルバック・ライブラー情報量［Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ Ｄｉｖｅｒｇｅｎｃｅ；ＫＬ−ｄｉｖ］）を最小化する変分パラメータを算出する。ＫＬ−ｄｉｖ Ｄ_ＫＬ（ｑ｜ｑ’）は、式（１８）に示すように、ある確率分布ｑ（ｊ）と他の確率分布ｑ’（ｊ）の差異を表す指標値である。

変分ベイズ法では、全ての未知変数についての変分事後分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝）が、未知変数Ｚの変分事後分布ｑ（Ｚ）と未知変数｛π｝、｛θ｝、｛Λ｝の変分事後分布ｑ（｛π｝，｛θ｝，｛Λ｝）の積に因数分解できると仮定する。
後述するように、例えば、式（１９）で表される変分下限値（ｖａｒｉａｔｉｏｎａｌ ｌｏｗｅｒ ｂｏｕｎｄ）Ｌ（ｑ）が収束する変分パラメータを算出する。

式（１９）において、Ｅ_{Ｚ，｛π｝，｛θ｝，｛Λ｝}［…］は、変分事後分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝）に対する…の期待値を表す。
そして、変分下限値Ｌ（ｑ）を周波数特性重み値の対数値ｌｏｇｗ^ｓｐ _ｉ，ｋとして算出する。変分下限値Ｌ（ｑ）と変分事後分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝）と確率分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝｜Ｘ，ｐ^ｉ _ｋ，｛μ｝）の間のＫＬ−ｄｉｖの和は、常に式（１７）の右辺と一致するため、変分ベイズ法によりＫＬ−ｄｉｖを最小化することで、式（１７）の近似値として変分下限値Ｌ（ｑ）を代用する。
変分ベイズ法では、変分パラメータの算出において式（２０）、（２１）に示す関係が仮定されている。これらの関係は、ともに変分事後分布ｑ^＊（｛Ｚ｝）、ｑ^＊（｛π｝、｛θ｝、｛Λ｝）が最適であることを表す。

式（２０）は、潜在変数｛Ｚ｝の変分事後分布ｑ^＊（｛Ｚ｝）が、粒子ｉ毎の確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐ^ｉ _ｋ，｛μ｝）の、その他の未知変数｛π｝、｛θ｝、｛Λ｝の変分事後分布ｑ^＊（｛π｝、｛θ｝、｛Λ｝）に対する期待値に比例することを表す。この関係を満足するように変分パラメータを算出する処理は、ＶＢ−Ｅ（Ｅｘｐｅｃｔａｔｉｏｎ）ステップと呼ばれる。

式（２１）は、未知変数｛π｝、｛θ｝、｛Λ｝の変分事後分布ｑ^＊（｛π｝、｛θ｝、｛Λ｝）が、粒子ｉ毎の確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐ^ｉ _ｋ，｛μ｝）の、潜在変数Ｚの変分事後分布ｑ^＊（｛Ｚ｝）に対する期待値に比例することを表す。この関係を満足するように変分パラメータを算出する処理は、ＶＢ−Ｍ（Ｍａｘｉｍｉｚａｔｉｏｎ）ステップと呼ばれる。

次に、ＬＨＡ法に変分ベイズ法を適用したときの周波数特性重み算出部１２２の構成について説明する。
図９は、本実施形態に係る周波数特性重み算出部１２２の構成を表す概略図である。
周波数特性重み算出部１２２は、データ入力部１２２１、第１変分パラメータ算出部１２２２、第２変分パラメータ算出部１２２３、変分下限値算出部１２２４、収束判定部１２２５、及び重み算出部１２２６を含んで構成される。

データ入力部１２２１は、周波数特性分析部１０３から音響スペクトログラムＡ_ｆ，ｔが、楽譜情報入力部１１２から楽譜情報が、粒子遷移部１３２から粒子毎の粒子情報がそれぞれ入力される。データ入力部１２２１は、音響スペクトログラムＡ_ｆ，ｔを予め定めた量子化幅ΔＡで量子化して量子化音響スペクトログラム［Ａ_ｆ，ｔ］を算出する。ここで、算出した量子化音響スペクトログラム［Ａ_ｆ，ｔ］は、音響フレームｎにおける周波数ｘ_ｎ毎の度数（確率）としたヒストグラムを表す。つまり、このヒストグラムは周波数ｘ_ｎが発生する回数が量子化音響スペクトログラム［Ａ_ｆ，ｔ］が表す振幅に比例することを表す。

データ入力部１２２１は、粒子情報から粒子ｉ毎に楽譜位置情報ｐ^ｉ _ｋ、拍間隔情報ｂ^ｉ _ｋを抽出し、入力された楽譜情報から楽譜位置ｐ^ｉ _ｋから拍間隔ｂ^ｉ _ｋ前の楽譜位置ｐ^ｉ _ｋ−ｂ^ｉ _ｋまでの楽譜情報を抽出する。データ入力部１２２１は、抽出した楽譜情報からセグメントｄ毎の基本周波数ベクトルμ_ｌを抽出する。
データ入力部１２２１は、算出した量子化音響スペクトログラム［Ａ_ｆ，ｔ］から、粒子ｉ毎の楽譜位置ｐ^ｉ _ｋに対応する時刻ｋΔＴから楽譜位置ｐ^ｉ _ｋ−ｂ^ｉ _ｋに対応する時刻までの区間の量子化音響スペクトログラム［Ａ_ｆ，ｔ］を抽出する。
データ入力部１２２１は、粒子ｉ毎に抽出した区間の量子化音響スペクトログラム［Ａ_ｆ，ｔ］とセグメントｄ毎の基本周波数ベクトルμ_ｐを、第１変分パラメータ算出部１２２２、第２変分パラメータ算出部１２２３、及び変分下限値算出部１２２４に出力する。

第１変分パラメータ算出部１２２２は、データ入力部１２２１から、量子化音響スペクトログラム［Ａ_ｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμ_ｐが入力される。また、第２変分パラメータ算出部１２２３から変分パラメータα_ｄｌ、β_ｌｍ、ａ_ｌ、ｂ_ｌが入力される。
第１変分パラメータ算出部１２２２は、入力されたパラメータを用いて、上述のＶＢ−Ｅステップを実行する。即ち、第１変分パラメータ算出部１２２２は、例えば式（２２）を用いて、ｄ、ｎ、ｌ、ｍ毎にパラメータρ_ｄｎｌｍを算出する。

式（２２）において、ψ（…）は、…のディガンマ（ｄｉｇａｍｍａ）関数を表す。
α_ｄ・は、α_ｄｌのｌに対する総和Σ_ｌα_ｄｌである。β_・ｌは、β_ｄｌのｄに対する総和Σ_ｄβ_ｄｌである。
次に、第１変分パラメータ算出部１２２２は、式（２３）に基づいてｄ、ｎ、ｌ、ｍ毎にパラメータρ_ｄｎｌｍをパラメータρ_ｄｎ・・で除算して変分パラメータγ_ｄｎｌｍを算出する。

式（２３）において、パラメータρ_ｄｎ・・は、パラメータρ_ｄｎｌｍのｌ、ｍに対する総和Σ_ｌΣ_ｍρ_ｄｎｌｍである。
変分パラメータγ_ｄｎｌｍは、負担率とも呼ばれ、観測した周波数ｘ_ｎが調波構造ｌの第ｍ高調波に対応するガウス分布の事後確率を表す。
なお、式（２２）、（２３）は、潜在変数｛Ｚ｝の変分事後分布ｑ^＊（｛Ｚ｝）が、変分パラメータγ_ｄｎｌｍのｚ_ｄｎｌｍ乗の積からなる多項分布（式（２４））に従うことを仮定して導出される。変分パラメータγ_ｄｎｌｍのｚ_ｄｎｌｍ乗の積は、変分パラメータγ_ｄｎｌｍがｚ_ｄｎｌｍ回発生する確率である。

第１変分パラメータ算出部１２２２は算出した変分パラメータγ_ｄｎｌｍを第２変分パラメータ算出部１２２３に出力する。

第２変分パラメータ算出部１２２３は、データ入力部１２２１から、量子化音響スペクトログラム［Ａ_ｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμ_ｐが入力される。また、第１変分パラメータ算出部１２２２から変分パラメータγ_ｄｎｌｍが入力される。
第２変分パラメータ算出部１２２３は、入力されたパラメータを用いて、上述のＶＢ−Ｍステップを実行する。ここで、第２変分パラメータ算出部１２２３は、例えば式（２５）を用いて、ｄ、ｌ毎にパラメータα_０ｌにパラメータγ_ｄ・ｌ・を加算して変分パラメータα_ｄｌを算出する。

式（２５）において、パラメータγ_ｄ・ｌ・は、変分パラメータγ_ｄｎｌｍのｎ、ｍに対する総和Σ_ｌΣ_ｍγ_ｄｎｌｍである。パラメータα_０ｌは、予め定められた実数を表す。
第２変分パラメータ算出部１２２３は、例えば式（２６）を用いて、ｌ、ｍ毎にパラメータβ_０ｍにパラメータγ_・・ｌｍを加算して変分パラメータβ_ｌｍを算出する。

式（２６）において、パラメータγ_・・ｌｍは、変分パラメータγ_ｄｎｌｍのｄ、ｎに対する総和Σ_ｄΣ_ｎγ_ｄｎｌｍである。パラメータβ_０ｍは、予め定められた実数を表す。
第２変分パラメータ算出部１２２３は、例えば式（２７）を用いて、ｌ毎にパラメータａ_０にパラメータγ_・・ｌ・を加算して変分パラメータａ_ｌを算出する。

式（２７）において、パラメータγ_・・ｌ・は、変分パラメータγ_ｄｎｌｍのｄ、ｎ、ｍに対する総和Σ_ｄΣ_ｎΣ_ｍγ_ｄｎｌｍである。ａ_０は、予め定められた実数を表す。
第２変分パラメータ算出部１２２３は、例えば式（２８）を用いて、変分パラメータｂ_ｌを算出する。

式（２８）において、ｂ_０は、予め定められた実数を表す。
なお、式（２５）は、未知変数｛π_ｄ｝の変分事後分布ｑ^＊（｛π｝）が変分パラメータ｛α_ｄ｝に対するディリクレ分布の積（式（２９））に従うことを仮定して導出される。

式（２９）において、Ｄｉｒ（｛π_ｄ｝｜｛α_ｄ｝）は、｛π_ｄ｝の｛α_ｄ｝に対するディリクレ分布を表す。｛π_ｄ｝は、変数π_ｄｌのｄに対する集合（セット）を表す。｛α_ｄ｝は、変数α_ｄｌのｄに対するセットを表す。式（２９）のディリクレ分布は、複数の事象が各々α_ｄｌ回発生する確率がπ_ｄｌとなる確率を与える。
式（２６）は、未知変数｛θ_ｌ｝の変分事後分布ｑ^＊（｛θ｝）が変分パラメータ｛θ_ｌ｝に対するディリクレ分布の積（式（３０））に従うことを仮定して導出される。

式（３０）において、Ｄｉｒ（｛θ_ｌ｝｜｛β_ｌ｝）は、｛θ_ｌ｝の｛β_ｌ｝に対するディリクレ分布を表す。｛θ_ｌ｝は、変数θ_ｌｍのｌに対する集合（セット）を表す。｛β_ｌ｝は、変数β_ｌｍのｍに対する集合（セット）を表す。式（３０）のディリクレ分布は、複数の事象が各々β_ｄｌ回発生する確率がθ_ｌｍとなる確率を与える。
式（２７）、（２８）は、未知変数｛Λ｝の変分事後分布ｑ^＊（｛Λ｝）が変分パラメータａ_ｌ、ｂ_ｌに対するガンマ分布の積（式（３１））に従うことを仮定して導出される。

式（３１）において、Ｇａｍ（Λ_ｌ｜ａ_ｌ，ｂ_ｌ）は、変数Λ_ｌのガンマ関数を表す。ａ_ｌ、ｂ_ｌは、それぞれガンマ関数の形状パラメータ、比率パラメータを表す。なお、変分事後分布ｑ^＊（｛π｝、｛θ｝、｛Λ｝）は、式（２９）、（３０）、（３１）で各々表される変分事後分布の積ｑ^＊（｛π｝）、ｑ^＊（｛θ｝）ｑ^＊（｛Λ｝）で因数分解できることが仮定されている。
第２変分パラメータ算出部１２２３は、算出した変分パラメータα_ｄｌ、β_ｌｍ、ａ_ｌ、ｂ_ｌと、第１変分パラメータ算出部１２２２から入力された変分パラメータγ_ｄｎｌｍを変分下限値算出部１２２４に出力する。収束判定部１２２５から変分下限値Ｌ（ｑ）が収束したことを表す変分下限値信号が入力されたとき、第２変分パラメータ算出部１２２３は、第１変分パラメータ算出部１２２２への変分パラメータα_ｄｌ、β_ｌｍ、ａ_ｌ、ｂ_ｌの出力を停止する。

変分下限値算出部１２２４は、第２変分パラメータ算出部１２２３から変分パラメータα_ｄｌ、β_ｌｍ、ａ_ｌ、ｂ_ｌ、γ_ｄｎｌｍが入力される。変分下限値算出部１２２４は、データ入力部１２２１から入力された量子化音響スペクトログラム［Ａ_ｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμ_ｐに基づいて、例えば式（３２）に基づいて確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐ^ｉ _ｋ，｛μ｝）を算出する。

式（３２）において、Ｘは入力された区間の量子化音響スペクトログラム［Ａ_ｆ，ｔ］に基づく周波数ｘ_ｎの集合（即ち、上述のヒストグラム）を表す。｛μ｝は入力されたセグメントｄ毎の基本周波数ベクトル｛μ_ｐ｝（但し、楽譜フレームｐは、セグメントｄに属す）に含まれる基本周波数μ_ｌの集合（セット）を表す。ｑ（Ｘ｜Ｚ，｛μ｝,Λ）は、例えば、式（３３）に示されるように、平均がｍ次の基本周波数μ_ｌ、分散Λ_ｌである、周波数ｘ_ｎのガウス関数の積である。

式（３２）において、ｑ（Ｚ｜｛π｝,｛θ｝）は、例えば、式（３４）に示されるようにπ_ｄｌθ_ｌｍのｚ_ｄｎｌｍ乗の積からなる多項分布である。

式（３２）において、ｑ（｛π｝）は、例えば、式（３５）に示されるように変数｛π_ｄ｝のパラメータセット｛α_０｝に対するディリクレ分布の積である。

パラメータセット｛α_０｝は、上述のパラメータα_０ｌの集合（セット）である。
式（３２）において、ｑ（｛θ｝）は、例えば、式（３６）に示されるように変数｛θ_ｌ｝のパラメータセット｛β_０｝に対するディリクレ分布の積である。

パラメータセット｛β_０｝は、上述のパラメータβ_０ｍの集合（セット）である。
式（３２）において、ｑ（｛Λ｝）は、例えば、式（３７）に示されるように分散Λ_ｌのパラメータａ_０、ｂ_０に対するガンマ分布の積である。

パラメータａ_０、ｂ_０の値は、上述のパラメータａ_０、ｂ_０と同一の値である。
変分下限値算出部１２２４は、入力された変分パラメータα_ｄｌ、β_ｌｍ、ａ_ｌ、ｂ_ｌ、γ_ｄｎｌｍに基づいて、式（２９）、（３０）、（３１）、（２４）を用いて、変分事後分布ｑ^＊（｛π｝）、ｑ^＊（｛θ｝）、ｑ^＊（｛Λ｝）、ｑ^＊（Ｚ）を算出する。変分下限値算出部１２２４は、算出したｑ^＊（｛π｝）、ｑ^＊（｛θ｝）、ｑ^＊（｛Λ｝）、ｑ^＊（Ｚ）の積をとり全未知変数についての変分事後分布ｑ^＊（Ｚ，π，θ，Λ）を算出する。

変分下限値算出部１２２４は、算出した確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐ^ｉ _ｋ，｛μ｝）と変分事後分布ｑ^＊（Ｚ，π，θ，Λ）を、例えば式（１９）を用いて粒子ｉ毎の変分下限値Ｌ（ｑ）を算出する。変分下限値算出部１２２４は、式（１９）においてｑ（Ｚ，π，θ，Λ）の代わりに、算出したｑ^＊（Ｚ，π，θ，Λ）を代入する。
なお、変分事後分布ｑ^＊（｛π｝）、ｑ^＊（｛θ｝）、ｑ^＊（｛Λ｝）、ｑ^＊（Ｚ）は、それぞれ式（２９）、（３０）、（３１）、（２４）で表され、確率分布ｑ（Ｘ｜Ｚ，｛μ｝,Λ）ｑ（Ｚ｜｛π｝,｛θ｝）、ｑ（｛π｝）、ｑ（｛θ｝）、ｑ（｛Λ｝）は、それぞれ式（３３）、（３４）、（３５）、（３６）、（３７）で表される関数を仮定しているため、変分下限値算出部１２２４は、式（１９）を用いて変分下限値Ｌ（ｑ）を解析的に算出することができる。かかる仮定により、事前分布と変分分布の関数形は相似、即ち共役になる。
変分下限値算出部１２２４は、算出した粒子ｉの変分下限値Ｌ（ｑ）を収束判定部１２２５に出力する。

収束判定部１２２５は、変分下限値算出部１２２４から入力された粒子ｉ毎の変分下限値Ｌ（ｑ）が、収束したか否か判断する。収束判定部１２２５は、例えば、前回入力された変分下限値Ｌ（ｑ）との差分の絶対値が、予め定めた値よりも小さくなったとき、変分下限値Ｌ（ｑ）が収束したと判断する。
収束判定部１２２５は、変分下限値Ｌ（ｑ）が収束したと判断したとき、変分下限値Ｌ（ｑ）が収束したことを表す変分下限値信号を第２変分パラメータ算出部１２２３に出力する。また、収束判定部１２２５は、粒子ｉ毎の変分下限値Ｌ（ｑ）を重み算出部１２２６に出力する。

重み算出部１２２６は、収束判定部１２２５から入力された粒子ｉの変分下限値Ｌ（ｑ）に基づいて粒子ｉの周波数特性重み値ｗ^ｓｐ _ｉ，ｋを算出する。重み算出部１２２６は、例えば、ｗ^ｓｐ _ｉ，ｋ＝ｅｘｐ（Ｌ（ｑ））と算出する。重み算出部１２２６は、算出した周波数特性重み値ｗ^ｓｐ _ｉ，ｋを表す周波数特性重み情報を粒子重み算出部１２４に出力する。

次に、周波数特性重み算出部１２２が行う周波数特性重み算出処理について説明する。
図１０は、本実施形態における周波数特性重み算出処理の一例を表すフローチャートである。
（ステップＳ３０１）データ入力部１２２１は、周波数特性分析部１０３から入力された音響スペクトログラムＡ_ｆ，ｔを予め定めた量子化幅ΔＡで量子化して量子化音響スペクトログラム［Ａ_ｆ，ｔ］を算出する。
データ入力部１２２１は、粒子遷移部１３２から入力された粒子情報から粒子ｉ毎に楽譜位置情報、拍間隔情報を抽出する。データ入力部１２２１は、楽譜情報入力部１１２から入力された楽譜情報から楽譜位置ｐ^ｉ _ｋから拍間隔ｂ^ｉ _ｋ前の楽譜位置ｐ^ｉ _ｋ−ｂ^ｉ _ｋまでの楽譜情報を抽出し、抽出した楽譜情報からセグメントｄ毎の基本周波数ベクトルμ_ｌを抽出する。

データ入力部１２２１は、算出した量子化音響スペクトログラム［Ａ_ｆ，ｔ］から、粒子ｉ毎の楽譜位置ｐ^ｉ _ｋに対応する時刻ｋΔＴから楽譜位置ｐ^ｉ _ｋ−ｂ^ｉ _ｋに対応する時刻までの区間の量子化音響スペクトログラム［Ａ_ｆ，ｔ］を抽出する。
データ入力部１２２１は、粒子ｉ毎に抽出した区間の量子化音響スペクトログラム［Ａ_ｆ，ｔ］とセグメントｄ毎の基本周波数ベクトルμ_ｐを、第１変分パラメータ算出部１２２２、第２変分パラメータ算出部１２２３、及び変分下限値算出部１２２４に出力する。その後、ステップＳ３０２に進む。

（ステップＳ３０２）第１変分パラメータ算出部１２２２は、データ入力部１２２１から、量子化音響スペクトログラム［Ａ_ｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμ_ｐが入力される。また、第１変分パラメータ算出部１２２２は、第２変分パラメータ算出部１２２３から変分パラメータα_ｄｌ、β_ｌｍ、ａ_ｌ、ｂ_ｌが入力される。
第１変分パラメータ算出部１２２２は、入力されたパラメータに基づいて、例えば式（２２）を用いて、パラメータρ_ｄｎｌｍを算出し、式（２３）を用いて変分パラメータγ_ｄｎｌｍを算出する。
第１変分パラメータ算出部１２２２は算出した変分パラメータγ_ｄｎｌｍを第２変分パラメータ算出部１２２３に出力する。その後、ステップＳ３０３に進む。

（ステップＳ３０３）第２変分パラメータ算出部１２２３は、データ入力部１２２１から、量子化音響スペクトログラム［Ａ_ｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμ_ｐが入力される。また、第１変分パラメータ算出部１２２２から変分パラメータγ_ｄｎｌｍが入力される。第２変分パラメータ算出部１２２３は、入力された変分パラメータγ_ｄｎｌｍに基づいて、例えば式（２５）、（２６）、（２７）、（２８）を用いて、変分パラメータα_ｄｌ、β_ｌｍ、ａ_ｌ、ｂ_ｌを算出する。第２変分パラメータ算出部１２２３は、変分パラメータα_ｄｌ、β_ｌｍ、ａ_ｌ、ｂ_ｌを第１変分パラメータ算出部１２２２と変分下限値算出部１２２４に出力する。第２変分パラメータ算出部１２２３は、第１変分パラメータ算出部１２２２から入力された変分パラメータγ_ｄｎｌｍを変分下限値算出部１２２４に出力する。その後、ステップＳ３０４に進む。

（ステップＳ３０４）変分下限値算出部１２２４は、データ入力部１２２１から入力された量子化音響スペクトログラム［Ａ_ｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμ_ｐに基づいて、例えば式（３２）に基づいて確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐ^ｉ _ｋ，｛μ｝）を算出する。変分下限値算出部１２２４は、第２変分パラメータ算出部１２２３から入力された変分パラメータα_ｄｌ、β_ｌｍ、ａ_ｌ、ｂ_ｌ、γ_ｄｎｌｍに基づいて変分事後分布ｑ^＊（Ｚ，π，θ，Λ）を算出する。
変分下限値算出部１２２４は、算出した確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐ^ｉ _ｋ，｛μ｝）と変分事後分布ｑ^＊（Ｚ，π，θ，Λ）を、例えば式（１９）を用いて粒子ｉ毎の変分下限値Ｌ（ｑ）を算出する
変分下限値算出部１２２４は、算出した粒子ｉの変分下限値Ｌ（ｑ）を収束判定部１２２５に出力する。その後、ステップＳ３０５に進む。

（ステップＳ３０５）収束判定部１２２５は、変分下限値算出部１２２４から入力された粒子ｉ毎の変分下限値Ｌ（ｑ）が、収束したか否か判断する。ここで、収束判定部１２２５は、変分下限値算出部１２２４から前回入力された変分下限値Ｌ（ｑ）との差分の絶対値が、予め定めた値よりも小さくなったとき、変分下限値Ｌ（ｑ）が収束したと判断する。変分下限値Ｌ（ｑ）が収束したと判断したとき（ステップＳ３０５ Ｙ）、第２変分パラメータ算出部１２２３は、変分パラメータα_ｄｌ、β_ｌｍ、ａ_ｌ、ｂ_ｌの第１変分パラメータ算出部１２２２への出力を停止する。収束判定部１２２５は、粒子ｉ毎の変分下限値Ｌ（ｑ）を重み算出部１２２６に出力する。その後、ステップＳ３０６に進む。
変分下限値Ｌ（ｑ）が収束していないと判断したとき（ステップＳ３０５ Ｎ）、ステップＳ３０２に進む。

（ステップＳ３０６）重み算出部１２２６は、算出した粒子ｉの変分下限値Ｌ（ｑ）に基づいて粒子ｉの周波数特性重み値ｗ^ｓｐ _ｉ，ｋを算出する。重み算出部１２２６は、算出した周波数特性重み値ｗ^ｓｐ _ｉ，ｋを表す周波数特性重み情報を粒子重み算出部１２４に出力する。その後、処理を終了する。

以上、説明したように、本実施形態では、入力された音響信号の第１の周波数特性と楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性との関連度を表す重み値を、第２の周波数特性に含まれる成分毎の第１の周波数特性の確率分布に基づいて、楽譜の位置毎に算出する。また、本実施形態では、算出した重み値に基づいて入力された音響信号に対応する楽譜の位置を探索する。これにより、音響信号に基づく第１の周波数特性と、楽譜情報における音階に基づく第２の周波数特性が完全に合致していなくとも、両者の関連性を検出することができる。そのため、演奏された音楽に対する楽譜の位置をより頑健に推定することが可能になる。
また、本実施形態では、第２の周波数特性は少なくとも１つの調波構造を含み、調波構造に含まれる調波成分毎の調波構造への寄与を表す変数の確率分布と、調波構造毎の第２の周波数特性への寄与を表す変数の確率分布に基づいて第１の周波数特性についての尤度を最大にするように、重み値を算出する。そのため、楽譜情報に基づく第２の周波数成分に含まれる調波構造の周波数特性の変化、調波構造間の結合特性の変化に柔軟に対応できるので、音響信号に基づく第１の周波数特性との関連性をより的確に推定することができる。そのため、演奏された音楽に対する楽譜の位置をより正確に推定することが可能になる。

（第２の実施形態）
次に、本発明の第２の実施形態について説明する。第１の実施形態と同一の構成又は処理については、同一の番号を付す。以下、第１の実施形態との差異点を主として説明する。
図１１は、本実施形態に係る楽譜位置推定装置２の構成を示す概略図である。
楽譜位置推定装置２は、楽譜位置推定装置１（図１）と同様に、音響信号入力部１０１、音響特徴量生成部１０２、楽譜情報記憶部１１１、楽譜情報入力部１１２、関連度算出部１２１、楽譜位置探索部１３１、及び演奏情報出力部１４１を含んで構成される。但し、関連度算出部１２１は、周波数特性重み算出部２２２、拍間隔重み算出部１２３、及び粒子重み算出部１２４を含んで構成される。即ち、楽譜位置推定装置２は、周波数特性重み算出部１２２の代わりに周波数特性重み算出部２２２を備える点で、楽譜位置推定装置１と異なる。

本実施形態でも、楽譜情報に基づく振幅周波数特性（例えば、楽譜スペクトログラム）が、既知の楽譜情報に含まれる楽譜位置毎の音階の基本周波数に基づく調波構造で表される成分を少なくとも１つ含むと仮定する。そこで、周波数特性重み算出部２２２は、楽譜スペクトログラムに含まれる周波数毎の振幅が表す度数（振幅期待値）が、対応する周波数における音響スペクトログラムの振幅が表す度数だけ発生する確率分布に基づいて周波数特性重み値ｗ^ｓｐ _ｉ，ｋを算出する。周波数特性重み算出部２２２は、例えば、ＰＤＡ（Ｐｏｉｓｓｏｎ Ｄｉｓｔｒｉｂｕｔｉｏｎ Ａｍｐｌｉｔｕｄｅ；ポアソン分布振幅）を用いて周波数特性重み値ｗ_ｉ，ｋ ^ｓｐを算出する。但し、周波数特性重み算出部２２２は、楽譜スペクトログラムの振幅の周波数方向に累積した総和が、予め定めた値（例えば、１）となるように各周波数の振幅を規格化する。

ＰＤＡは、量子化音響スペクトログラムの振幅期待値（ｅｘｐｅｃｔｅｄ ａｍｐｌｉｔｕｄｅ）に対するポアソン分布の関数値を周波数にわたり累積した累積値である。振幅期待値は、楽譜スペクトログラムに、量子化音響スペクトログラムを周波数にわたり累積した総和を乗じた値である。また、上述のように、量子化音響スペクトログラムは、周波数毎の度数を表すヒストグラムとみることもできる。上述のポアソン分布は、楽譜情報による周波数毎の平均度数だけ発生する事象が、入力音響信号による周波数毎の度数だけ発生する確率を表す。従って、ＰＤＡは、音響スペクトログラムと楽譜スペクトログラムとの差異の度合いを表す指標値を表す。例えば、ＰＤＡは、音響スペクトログラムのピーク周波数と合致する楽譜スペクトログラムのピーク周波数があれば大きい値をとる。ＰＤＡは、音響スペクトログラムのピーク周波数と合致する楽譜スペクトログラムがなければ小さい値をとる。そのため、ＰＤＡは、音響スペクトログラムと楽譜スペクトログラムの調波構造の差異に対して鋭敏に変化する。

次に、周波数特性重み算出部２２２の構成について説明する。
図１２は、本実施形態に係る周波数特性重み算出部２２２の構成を表す概略図である。
周波数特性重み算出部１２２は、周波数特性量子化部２２２１、周波数特性合成部２２２２、振幅期待値算出部２２２３、振幅分布算出部２２２４及び重み算出部２２２５を含んで構成される。

周波数特性量子化部２２２１は、周波数特性分析部１０３から入力された音響スペクトログラムＡ_ｆ，ｔを予め定めた量子化幅ΔＡで量子化して量子化音響スペクトログラム［Ａ_ｆ，ｔ］を算出する。周波数特性量子化部２２２１は、算出した量子化音響スペクトログラム［Ａ_ｆ，ｔ］を振幅期待値算出部２２２３及び振幅分布算出部２２２４に出力する。

周波数特性合成部２２２２は、楽譜情報入力部１１２から入力された楽譜情報から楽譜フレームｐ毎の基本周波数ベクトルμ_ｐを抽出し、基本周波数ベクトルμ_ｐに含まれる基本周波数μ_ｐ ^ｌ毎に予め定めた音源モデル（例えば、ＧＭＭ）を用いて調波構造を生成する。周波数特性合成部２２２２は、生成した調波構造を基本周波数μ_ｐ ^ｌにわたって累積して楽譜スペクトログラムＡ＾_ｆ，ｐを生成する。生成した楽譜スペクトログラムＡ＾_ｆ，ｐは、周波数にわたり累積した総和Σ_ｆＡ＾_ｆ，ｐが１となるように規格化されている。
周波数特性合成部２２２２は、例えば式（３８）を用いて楽譜スペクトログラムＡ＾_ｆ，ｐを生成する。

式（３８）において、Ｃ_ｈａｒｍは、例えば、式（３９）に示すように、楽譜スペクトログラムＡ＾_ｆ，ｐの周波数毎の振幅値を周波数にわたり累積した総和が、０より大きく１より小さい実数（式（３９）では、例えば０．９）を与える定数である。

式（３８）において、Ｍ、θ_ｌｍ、Λ_ｌは、例えば、それぞれ１０、０．２^{（ｍ−１）}、０．８である。Ｃ_{ｆｌｏｏｒ}は、Ｃ_ｈａｒｍとの和が１となる実数（例えば、０．１）である。Ｃ_{ｆｌｏｏｒ}を定めることにより、楽譜スペクトログラムＡ＾_ｆ，ｐは、０より大きい実数値をとるため、後述する処理において零除算が回避される。
周波数特性合成部２２２２は、生成した楽譜スペクトログラムＡ＾_ｆ，ｐを振幅期待値算出部２２２３に出力する。

振幅期待値算出部２２２３は、周波数特性量子化部２２２１から量子化音響スペクトログラム［Ａ_ｆ，ｔ］が、周波数特性合成部２２２２から楽譜スペクトログラムＡ＾_ｆ，ｐが、粒子遷移部１３２から粒子毎の粒子情報がそれぞれ入力される。
振幅期待値算出部２２２３は、入力された粒子情報から粒子ｉ毎に楽譜位置情報ｐ^ｉ _ｋを抽出する。振幅期待値算出部２２２３は、抽出した楽譜位置ｐ^ｉ _ｋを基準（例えば、楽譜位置ｐ^ｉ _ｋに対応する時刻τがｋΔＴ）とした対応する時刻τ毎に、量子化音響スペクトログラム［Ａ_ｆ，τ］の周波数ｆにわたった総和Ａ_・τを算出する。
振幅期待値算出部２２２３は、入力された楽譜スペクトログラムＡ＾_ｆ，ｐを、対応する楽譜位置ｐに対応する時刻τ毎に算出した総和Ａ_・τを乗算して振幅期待値λ_ｆ，τを算出する。振幅期待値算出部２２２３は、算出した振幅期待値λ_ｆ，τを振幅分布算出部２２２４に出力する。

振幅分布算出部２２２４は、周波数特性量子化部２２２１から量子化音響スペクトログラム［Ａ_ｆ，ｔ］が、振幅期待値算出部２２２３から振幅期待値λ_ｆ，τ、粒子遷移部１３２から粒子毎の粒子情報が、それぞれ入力される。
振幅分布算出部２２２４は、入力された粒子ｉ毎の粒子情報に基づき楽譜位置ｐ^ｉ _ｋを基準とした対応する時刻τ及び周波数ｆ毎に、振幅期待値λ_ｆ，τの量子化音響スペクトログラム［Ａ_ｆ，τ］に対するポアソン分布の関数値Ｐｏｉ（［Ａ_ｆ，τ］｜λ_ｆ，τ）を算出する。ポアソン分布は、式（４０）で表される確率密度関数である。

式（４０）において、ｋ_ｘは、０又は０よりも大きい整数、λは、０よりも大きい実数である。
算出した関数値は、ポアソン分布の性質により、周波数ｆにおける量子化音響スペクトログラム［Ａ_ｆ，τ］の量子化した振幅が表す度数に対する楽譜スペクトログラムＡ＾_{ｆ，ｐ（τ）}に含まれる振幅の分布関数の値を表す。即ち、算出した関数値は、周波数ｆにおける楽譜スペクトログラムＡ＾_ｆ，ｐの振幅が表す度数が、量子化音響スペクトログラム［Ａ_ｆ，τ］の量子化した振幅が表す度数だけ発生する確率値を表す。
振幅分布算出部２２２４は、算出した振幅分布値を重み算出部２２２５に出力する。

重み算出部２２２５は、入力された振幅分布値を周波数ｆ及び時刻τにわたり累積してＰＤＡを算出する。累積する時刻τの範囲は、粒子毎の観測区間Ｔ_ｋ｛τ｜ｋΔＴ−Ｗ＜τ≦ｋΔＴ｝である。重み算出部２２２５は、算出したＰＤＡの指数関数値を周波数特性重み値ｗ^ｓｐ _ｉ，ｋを算出する。ここで、重み算出部２２２５は、式（４１）を用いて周波数特性重み値ｗ^ｓｐ _ｉ，ｋを算出する。

重み算出部２２２５は、算出した周波数特性重み値ｗ^ｓｐ _ｉ，ｋを表す周波数特性重み情報を粒子重み算出部１２４に出力する。

次に、本実施形態に係る楽譜位置推定処理について説明する。本実施形態に係る楽譜位置推定処理は、図３に示す楽譜位置推定処理と同様である。但し、ステップＳ１０６において周波数特性重み算出部２２２は、次に説明する処理を実行する。
図１３は、本実施形態における周波数特性重み算出処理を表すフローチャートである。
（ステップＳ４０１）周波数特性量子化部２２２１は、周波数特性分析部１０３から入力された音響スペクトログラムＡ_ｆ，ｔを予め定めた量子化幅ΔＡで量子化して量子化音響スペクトログラム［Ａ_ｆ，ｔ］を算出する。周波数特性量子化部２２２１は、算出した量子化音響スペクトログラム［Ａ_ｆ，ｔ］を振幅期待値算出部２２２３及び振幅分布算出部２２２４に出力する。その後、ステップＳ４０２に進む。

（ステップＳ４０２）周波数特性合成部２２２２は、楽譜情報入力部１１２から入力された楽譜情報から楽譜フレームｐ毎の基本周波数ベクトル｛μ_ｐ｝を抽出し、基本周波数ベクトル｛μ_ｐ｝に含まれる基本周波数μ_ｐ ^ｌ毎に予め定めた音源モデル（例えば、ＧＭＭ）を用いて調波構造を生成する。周波数特性合成部２２２２は、生成した調波構造を基本周波数μ_ｐ ^ｌにわたって累積し、例えば式（３８）を用いて楽譜スペクトログラムＡ＾_ｆ，ｐを生成する。
周波数特性合成部２２２２は、生成した楽譜スペクトログラムＡ＾_ｆ，ｐを振幅期待値算出部２２２３に出力する。その後、ステップＳ４０３に進む。

（ステップＳ４０３）振幅期待値算出部２２２３は、粒子遷移部１３２から入力された粒子情報から粒子ｉ毎に楽譜位置情報を抽出する。振幅期待値算出部２２２３は、抽出した楽譜位置情報が表す楽譜位置ｐ^ｉ _ｋを基準とした対応する時刻τ毎に、周波数特性量子化部２２２１から入力された量子化音響スペクトログラム［Ａ_ｆ，τ］の周波数にわたった総和Ａ_・τを算出する。振幅期待値算出部２２２３は、周波数特性合成部２２２２から入力された楽譜スペクトログラムＡ＾_ｆ，ｐを、対応する楽譜位置ｐに対応する時刻τ毎に算出した総和Ａ_・τを乗算して振幅期待値λ_ｆ，τを算出する。振幅期待値算出部２２２３は、算出した振幅期待値λ_ｆ，τを振幅分布算出部２２２４に出力する。その後、ステップＳ４０４に進む。

（ステップＳ４０４）振幅分布算出部２２２４は、周波数特性量子化部２２２１から量子化音響スペクトログラム［Ａ_ｆ，ｔ］が、振幅期待値算出部２２２３から振幅期待値λ_ｆ，τ、粒子遷移部１３２から粒子毎の粒子情報がそれぞれ入力される。
振幅分布算出部２２２４は、粒子遷移部１３２から入力された粒子ｉ毎の粒子情報に含まれる楽譜位置情報を抽出する。振幅分布算出部２２２４は、抽出した楽譜位置情報が表す楽譜位置ｐ^ｉ _ｋを基準とした対応する時刻τ及び周波数ｆ毎に、振幅期待値λ_ｆ，τの量子化音響スペクトログラム［Ａ_ｆ，τ］に対するポアソン分布Ｐｏｉ（［Ａ_ｆ，τ］｜λ_ｆ，τ）の関数値を振幅分布値として算出する。
振幅分布算出部２２２４は、算出した振幅分布値を重み算出部２２２５に出力する。その後、ステップＳ４０５に進む。

（ステップＳ４０５）重み算出部２２２５は、振幅分布算出部２２２４から入力された振幅分布値を周波数ｆ及び粒子毎の観測区間Ｔ_ｋ内の時刻τにわたり累積してＰＤＡを算出する。重み算出部２２２５は、算出したＰＤＡの指数関数値を、例えば式（４１）を用いて周波数特性重み値ｗ^ｓｐ _ｉ，ｋを算出する。重み算出部２２２５は、算出した周波数特性重み値ｗ^ｓｐ _ｉ，ｋを表す周波数特性重み情報を粒子重み算出部１２４に出力する。その後、処理を終了する

以上、説明したように、本実施形態では、第２の周波数特性は少なくとも１つの調波構造を含み、調波構造に含まれる調波成分毎の調波構造への寄与を表す変数の確率分布と、調波構造毎の第２の周波数特性への寄与を表す変数の確率分布に基づいて音響信号の第１の周波数特性についての尤度を最大にするように、楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を算出する。
これにより、音響信号に基づく第１の周波数特性の周波数毎の振幅と楽譜情報における音階に基づく第２の周波数特性の周波数毎の振幅が合致する度合いを鋭敏に検出することができる。そのため、演奏された音楽に対する楽譜の位置をより正確に推定することが可能になる。

（第３の実施形態）
次に、本発明の第３の実施形態について説明する。第１の実施形態と同一の構成又は処理については、同一の番号を付す。以下、第１の実施形態との差異点を主として説明する。
図１４は、本実施形態に係る楽譜位置推定装置３の構成を示す概略図である。
楽譜位置推定装置３は、楽譜位置推定装置１（図１）と同様に、音響信号入力部１０１、音響特徴量生成部１０２、楽譜情報記憶部１１１、楽譜情報入力部１１２、関連度算出部１２１、楽譜位置探索部１３１、及び演奏情報出力部１４１を含んで構成される。
但し、音響特徴量生成部１０２は、周波数特性分析部１０３及び相関算出部１０４の他にクロマベクトル生成部３０５を備える点が楽譜位置推定装置１とは異なる。周波数特性分析部１０３は、算出した音響スペクトログラムＡ_ｆ，ｔをクロマベクトル生成部３０５に出力する。

クロマベクトル生成部３０５は、周波数特性分析部１０３から入力された音響スペクトログラムＡ_ｆ，ｔに基づいて、音響信号によるクロマベクトル（以下、音響クロマベクトルと呼ぶ）ｃ_ｔ ^ａを生成する。クロマベクトルとは、音階（ｃｈｒｏｍａ）毎の成分のパワーを要素として含むベクトルである。クロマベクトルは、例えば、西洋音楽を構成する１２音階（Ｃ、Ｃ＃、Ｄ、Ｄ＃、Ｅ、Ｆ、Ｆ＃、Ｇ、Ｇ＃、Ａ、Ａ＃、Ｂ）における音階毎の成分の強度を要素とするベクトル（例えば、要素数が１２）である。
クロマベクトル生成部３０５は、音響クロマベクトルｃ_ｔ ^ａ＝［ｃ_ｔ，１ ^ａ，ｃ_ｔ，２ ^ａ，…，ｃ_ｔ，１２ ^ａ］の音階ｊ毎の要素ｃ_ｔ，ｊ ^ａを、例えば、式（４２）を用いて算出する。

式（４２）において、ｏは、オクターブを表す整数値である。Ｏｃｔ_ｌｏｗは予め設定されたオクターブの下限を表す整数値（下限オクターブ、例えば３）を表す。Ｏｃｔ_ｈｉは、予め設定されたオクターブの上限を表す整数値（上限オクターブ、例えば６）を表す。Ｂ_ｊ，ｏ（ｆ）は、オクターブｏにおける音階ｊの成分を抽出する帯域通過フィルタ（バンドパスフィルタ、ｂａｎｄｐａｓｓ ｆｉｌｔｅｒ；ＢＰＦ）の入出力振幅周波数特性を表す。ｏ番目のオクターブにおける音階ｊの周波数ｆ_ｊ，ｏにおいて、Ｂ_ｊ，ｏ（ｆ）は、例えば、式（４３）に示されるように最大となり、周波数ｆがゼロ又は無限大に近づくと、ゼロに近似する関数である。

式（４３）において、ｆ^ｃｅｎｔ、ｆ_ｊ，ｏ ^ｃｅｎｔは、それぞれ周波数ｆ、ｆ_ｊ，ｏを対数領域で表したセント値である。セント値ｆ^ｃｅｎｔと周波数ｆは、式（４４）に示す関係がある。

式（４４）では、例えば、音階Ａ４（ｊ＝１、ｏ＝４）の基本周波数ｆが４４０Ｈｚであることが仮定されている。
クロマベクトル生成部３０５は、生成した音響クロマベクトルｃ_ｔ ^ａを、後述するクロマベクトル重み算出部３２５に出力する。

楽譜情報記憶部１１１に記憶されている楽譜情報には、楽曲を構成する各コードを構成する音階を要素として含むクロマベクトル（以下、楽譜クロマベクトルと呼ぶ）ｃ_ｐ ^ｓが、基本周波数ベクトルμ_ｐと対応付けられて含まれている。楽譜クロマベクトルの要素数は、１２である。つまり、楽譜クロマベクトルを構成する音階は、楽曲を表す音階のオクターブが区別されていない点で上述のコードベクトルを構成する音階と異なる。また、楽譜クロマベクトルｃ_ｐ ^ｓを構成する要素ｃ_ｐ，ｊ ^ｓは、０または１の値をとる。０とは、楽譜フレームｐにおいて音階ｊが存在しないことを表し、１とは、音階ｊが存在することを表す。例えば、音階Ａ４、Ｃ４、及びＥ４を要素値として含む楽譜クロマベクトルは、［１，０，０，１，０，０，０，１，０，０，０，０］となる。
楽譜情報入力部１１２は、楽譜情報記憶部１１１から読み出した楽譜情報に含まれる楽譜クロマベクトルｃ_ｐ ^ｓをクロマベクトル重み算出部３２５に出力する。

関連度算出部１２１は、周波数特性重み算出部１２２及び拍間隔重み算出部１２３の他にクロマベクトル重み算出部３２５を備え、粒子重み算出部１２４の代わりに粒子重み算出部３２４を備える。
なお、粒子遷移部１３２は、粒子ｉ毎の粒子情報をクロマベクトル重み算出部３２５に出力する。粒子遷移部１３２は、粒子ｉ毎の粒子分布値、状態遷移確率を粒子重み算出部３２４に出力する。

クロマベクトル重み算出部３２５は、クロマベクトル生成部３０５から入力された音響クロマベクトルｃ_ｔ ^ａを、時刻ｔ毎にノルムが１となるように各要素値を規格化する。クロマベクトル重み算出部３２５は、楽譜情報入力部１１２から入力された楽譜クロマベクトルｃ_ｐ ^ｓを、楽譜フレームｐ毎にノルムが１となるように各要素値を規格化する。
クロマベクトル重み算出部３２５は、粒子遷移部１３２から入力された粒子情報から粒子ｉ毎に楽譜位置情報ｐ^ｉ _ｋを抽出する。クロマベクトル重み算出部３２５は、要素値を規格化した音響クロマベクトルｃ_τ ^ａと楽譜クロマベクトルｃ_{ｐ（τ）’} ^ｓとの内積を、例えば式（４５）を用いて平均して、粒子ｉ毎のクロマベクトル重みｗ_ｉ，ｋ ^ｃｈを算出する。内積は、ベクトル間の類似性を表す指標値である。

式（４５）において、Ｔ_ｋは、抽出した楽譜位置ｐ^ｉ _ｋを基準（時刻τがｋΔＴ）とした粒子毎の観測区間Ｔ_ｋ｛τ｜ｋΔＴ−Ｗ＜τ≦ｋΔＴ｝を表す。
クロマベクトル重み算出部３２５は、算出した粒子ｉ毎のクロマベクトル重みｗ_ｉ，ｋ ^ｃｈを表すクロマベクトル重み情報を粒子重み算出部３２４に出力する。

粒子重み算出部３２４は、周波数特性重み算出部１２２から周波数特性重み情報が、拍間隔重み算出部１２３から拍間隔重み情報が、クロマベクトル重み算出部３２５からクロマベクトル重み情報が、粒子遷移部１３２から粒子ｉ毎の粒子分布値、状態遷移確率が入力される。
粒子重み算出部３２４は、式（４６）に示すように、周波数特性重み情報が表す周波数特性重み値、拍間隔重み情報が表す拍間隔重み値、クロマベクトル重み情報が表すクロマベクトル重み値、状態遷移確率を乗算し、さらに状態遷移確率で除算して、粒子ｉ毎の重み値ｗ_ｉ，ｋを算出する。

粒子重み算出部３２４は、算出した重み値ｗ_ｉ，ｋを表す粒子重み情報を生成する。
粒子重み算出部３２４は、粒子遷移部１３２から入力された粒子ｉ毎の粒子情報に含まれる粒子重み情報を、生成した粒子重み情報に置き換えることで粒子情報を更新する。
粒子重み算出部３２４は、更新した粒子情報を再標本化部１３３、楽譜位置算出部１４２及び拍間隔算出部１４３に出力する。
なお、本実施形態に係る楽譜位置推定装置３は、周波数特性重み算出部１２２の代わりに周波数特性重み算出部２２２を備えてもよい。

次に、本実施形態に係る楽譜位置推定処理について説明する。
図１５は、本実施形態に係る楽譜位置推定処理を表すフローチャートである。
図１５が示す楽譜位置推定処理は、ステップＳ１０１、Ｓ１０２−Ｓ１０５、Ｓ１０６、Ｓ１０７、Ｓ１０９−Ｓ１１７を有する点で、図３が示す楽譜位置推定処理と共通する。但し、ステップＳ１０１とステップＳ１０２の間でステップＳ４０１を実行し、ステップＳ１０５とステップＳ１０６の間でステップＳ４０２を実行し、ステップＳ１０８の代わりにステップＳ４０３、Ｓ４０４を実行する点が、図３が示す楽譜位置推定処理とは異なる。以下、処理が異なるステップについて説明する。

（ステップＳ４０１）音源位置推定装置３では、楽曲を構成する各コードを構成する音階を要素として含む楽譜情報に基づいて楽譜クロマベクトルｃ_ｐ ^ｓを予め生成しておく。音源位置推定装置３では、基本周波数ベクトルμ_ｐと対応付けられ楽譜情報の一部として予め楽譜情報記憶部１１１に記憶させておく。これにより、楽譜情報入力部１１２が、楽譜情報記憶部１１１から読み出した楽譜情報から楽譜クロマベクトルｃ_ｐ ^ｓをクロマベクトル重み算出部３２５に出力できるようにする。その後、ステップＳ１０２に進む。

（ステップＳ４０２）クロマベクトル生成部３０５は、周波数特性分析部１０３から入力された音響スペクトログラムＡ_ｆ，ｔに基づいて、例えば式（４２）を用いて算出される要素を有する音響クロマベクトルｃ_ｔ ^ａを生成する。
クロマベクトル生成部３０５は、生成した音響クロマベクトルｃ_ｔ ^ａをクロマベクトル重み算出部３２５に出力する。その後、ステップＳ１０６に進む。

（ステップＳ４０３）クロマベクトル重み算出部３２５は、クロマベクトル生成部３０５から入力された音響クロマベクトルｃ_ｔ ^ａを、時刻ｔ毎にノルムが１となるように各要素値を規格化する。クロマベクトル重み算出部３２５は、楽譜情報入力部１１２から入力された楽譜クロマベクトルｃ_ｐ ^ｓを、楽譜フレームｐ毎にノルムが１となるように各要素値を規格化する。
クロマベクトル重み算出部３２５は、粒子遷移部１３２から入力された粒子情報から粒子ｉ毎に楽譜位置情報を抽出する。クロマベクトル重み算出部３２５は、要素値を規格化した音響クロマベクトルｃ_τ ^ａと楽譜クロマベクトルｃ_ｐ（τ） ^ｓとの内積を、例えば式（４５）を用いて抽出した楽譜位置情報が表す楽譜位置を基準とした観測区間Ｔ_ｋで平均する。これにより、クロマベクトル重み算出部３２５は、粒子ｉ毎のクロマベクトル重みｗ_ｉ，ｋ ^ｃｈを算出する。クロマベクトル重み算出部３２５は、算出した粒子ｉ毎のクロマベクトル重みｗ_ｉ，ｋ ^ｃｈを表すクロマベクトル重み情報を粒子重み算出部１２４に出力する。その後、ステップＳ４０４に進む。

（ステップＳ４０４）粒子重み算出部３２４は、周波数特性重み算出部１２２から周波数特性重み情報が、拍間隔重み算出部１２３から拍間隔重み情報が、クロマベクトル重み算出部３２５からクロマベクトル重み情報が、粒子遷移部１３２から粒子ｉ毎の粒子分布値、状態遷移確率が入力される。
粒子重み算出部３２４は、周波数特性重み情報が表す周波数特性重み値、拍間隔重み情報が表す拍間隔重み値、クロマベクトル重み情報が表すクロマベクトル重み値、状態遷移確率を乗算し、さらに状態遷移確率で除算して、粒子ｉ毎の重み値ｗ_ｉ，ｋを算出する。
粒子重み算出部３２４は、算出した重み値ｗ_ｉ，ｋを表す粒子重み情報を生成する。
粒子重み算出部３２４は、粒子遷移部１３２から入力された粒子ｉ毎の粒子情報に含まれる粒子重み情報を、生成した粒子重み情報に置き換えることで粒子情報を更新する。
粒子重み算出部３２４は、更新した粒子情報を楽譜位置探索部１３１の再標本化部１３３、演奏情報出力部１４１の楽譜位置算出部１４２及び拍間隔算出部１４３に出力する。その後、ステップＳ１０９に進む。以下、図３に示す楽譜位置推定処理と同様にステップＳ１１０−１１７を実行する。
なお、本実施形態に係る楽譜位置推定処理では、ステップＳ１０６において周波数特性重み算出部１２２が図１０に示す周波数特性重み算出処理を実行する代わりに、周波数特性重み算出部２２２が図１３に示す周波数特性重み算出処理を実行してもよい。

次に、本実施形態に係る楽譜位置推定装置２が推定した楽譜位置の検証例について説明する。検証においては特に断らない限り、周波数の上限ｆ_ｍａｘを６０００Ｈｚ、テンポ窓長６０／ｂ^θを１５（ｂｐｍ）とした。

図１６は、楽譜位置の推定誤差の一例を表す図である。
図１６において、縦軸は、推定誤差、横軸は、楽曲の番号を表す。図１５において、○印は、本実施形態による平均推定誤差、×印は、従来の楽譜位置推定システムであるＡｎｔｅｓｃｏｆｏによる平均推定誤差を表す（例えば、Ｃｏｎｔ．Ａ．：Ａ Ｃｏｕｐｌｅｄ Ｄｕｒａｔｉｏｎ−Ｆｏｃｕｓｅｄ Ａｒｃｈｉｔｅｃｔｕｒｅ ｆｏｒ Ｒｅａｌｔｉｍｅ Ｍｕｓｉｃ ｔｏ Ｓｃｏｒｅ Ａｌｉｇｎｍｅｎｔ，ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｐａｔｔｅｒｎ Ａｎａｌｙｓｉｓ ａｎｄ Ｍａｃｈｉｎｅ Ｉｎｔｅｌｌｉｇｅｎｃｅ， Ｖｏｌ．３２，Ｎｏ．６，ｐｐ．９７４−９８７（２０１０））。Ａｎｔｅｓｃｏｆｏは、ＨＳＭＭ（Ｈｉｄｄｅｎ Ｓｅｍｉ−Ｍａｒｋｏｖ Ｍｏｄｅｌ；隠れセミマルコフモデル）に基づいて楽譜位置を推定する処理を行う。○印、×印を中心として上下に延びるバーの上限と下限は、それぞれ平均推定誤差から標準偏差だけ大きい推定誤差と、平均推定誤差から標準偏差だけ小さい推定誤差を表す。即ち、バーは、推定誤差の信頼区間を表す。楽曲１−２０は、各々テンポ及び演奏する楽器が異なる楽曲である。推定誤差は、式（１１）を用いて算出した値である。

図１６によれば、本実施形態では、ほとんどの楽曲についてＡｎｔｅｓｃｏｆｏよりも推定誤差の平均値及び標準偏差が減少する。本実施形態の楽曲にわたる平均推定誤差の絶対値の平均値は、Ａｎｔｅｓｃｏｆｏの平均推定誤差の絶対値の平均値よりも６９％減少する。即ち、本実施形態の方が、Ａｎｔｅｓｃｏｆｏよりも推定精度が優れる。一部の楽曲、例えば楽曲６、１１では、本実施形態の平均推定誤差の絶対値の方が、Ａｎｔｅｓｃｏｆｏの平均推定誤差の絶対値よりも大きい。しかし、本実施形態の標準偏差の方が、Ａｎｔｅｓｃｏｆｏの標準偏差よりも小さい。即ち、本実施形態の方が、より頑健に楽譜位置を推定できることが示されている。

図１７は、楽譜位置の推定誤差のその他の例を表す図である。
図１７において、横軸と縦軸の関係は、図１６と同様である。図１７において、○印、×印、△印は、テンポ窓長６０／ｂ^θが、それぞれ５、１５、３０（ｂｐｍ）である場合の平均推定誤差を表す。○印、×印、△印を中心として上下に延びるバーは、平均推定誤差を中心とし、その長さが標準偏差の２倍である信頼区間を表す。
図１７によれば、ほとんどの楽曲について、テンポ窓長が小さくなるほど平均推定誤差の絶対値が小さくなり、標準偏差も小さくなる傾向がある。即ち、テンポ窓長が小さいほど楽譜位置の推定精度が優れ、楽譜位置を安定して推定できることが示されている。

図１８は、楽譜位置の推定誤差のその他の例を表す図である。
図１８において、縦軸は、推定誤差の絶対値を表し、横軸は、楽曲の番号を表す。棒グラフは、推定処理又はシステム毎の平均推定誤差を表す。最も右側の網掛けの棒グラフは、第１の実施形態においてＬＨＡ法を用いて周波数特性重み値ｗ_ｉ，ｋ ^ｓｐを算出する処理を表す。右から２番目の縦縞で塗りつぶした棒グラフは、第２の実施形態においてＰＤＡを用いて周波数特性重み値ｗ_ｉ，ｋ ^ｓｐを算出する処理を表す。左から２番目の横縞で塗りつぶした棒グラフは、第２の実施形態の構成においてＰＤＡの代わりにカルバック・ライブラー情報量（ＫＬ−ｄｉｖ）を用いて周波数特性重み値ｗ_ｉ，ｋ ^ｓｐを算出する処理を表す。最も左側の斜めの縞模様で塗りつぶした棒グラフは、Ａｎｔｅｓｃｏｆｏが採用するＨＳＭＭを表す。また、各棒グラフの頂上から上に伸びるエラーバーの長さは、標準偏差を表す。
図１８によれば、楽曲１−１０のいずれにおいても、ＬＨＡ、ＰＤＡ、ＫＬ−ｄｉｖ、ＨＳＭＭの順で、平均推定誤差、標準偏差が、ともに小さくなる。即ち、この順序で、楽譜位置の推定精度が優れ、安定して楽譜位置を推定できることが示される。

図１９は、音響信号と楽譜情報との関連性の一例を表す図である。
図１９において、縦軸は対数尤度を、横軸は楽譜位置偏差を表す。図１９が表す対数尤度は、ある区間の楽譜情報と、ある観測区間の音響信号との間の周波数特性重み値ｗ_ｉ，ｋ ^ｓｐの対数値に相当する値である。即ち、対数尤度は、音響信号と楽譜情報との関連性を表す指標値である。但し、図１９に示す対数尤度は、計算精度を確保するために規格化されていない値である。楽譜位置偏差とは、楽譜情報の区間と対応する観測区間を基準とした、音響信号の観測区間のずれを、拍単位で表す値である。一点破線は、第１の実施形態においてＬＨＡ法を用いて周波数特性重み値ｗ_ｉ，ｋ ^ｓｐを算出する処理を表す。破線は、第２の実施形態においてＰＤＡを用いて周波数特性重み値ｗ_ｉ，ｋ ^ｓｐを算出する処理を表す。実線は、ＫＬ−ｄｉｖを用いて周波数特性重み値ｗ_ｉ，ｋ ^ｓｐを算出する処理を表す。
表示形態が異なる３本の太線は、各々楽譜位置偏差毎の対数尤度を表し、３本の細線は、各々対数尤度の最大値を表す。いずれも、楽譜位置偏差が０の場合、対数尤度が最大になる。即ち、楽譜位置偏差が０とは、楽譜位置と音響信号の時刻の対応が取れている点であることを表す。

図１９において、楽譜位置偏差による対数尤度の変化は、ＬＨＡ、ＰＤＡ、ＫＬ−ｄｉｖの順で著しい。とりわけ、第１の実施形態においてＬＨＡを用いた場合や、第２の実施形態においてＰＤＡを用いた場合の対数尤度の鋭敏な変化は、ＫＬ−ｄｉｖを用いた場合の緩慢な対数尤度の変化と対照的である。このことは、上述の実施形態において音源モデルとしてＬＨＡ又はＰＤＡを用い、最尤推定を行うことで楽譜位置が正確に推定できることを示す。

以上、説明したように、上述した実施形態では、入力された音響信号に基づく第１の周波数特性と楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、第２の周波数特性に含まれる成分毎の第１の周波数特性の確率分布に基づいて、楽譜の位置毎に算出する。また、上述した実施形態では、算出した重み値に基づいて音響信号に対応する楽譜の位置を探索する。
これにより、音響信号に基づく第１の周波数特性と、楽譜情報における音階に基づく第２の周波数特性が完全に合致していなくとも、両者の関連性を検出することができる。そのため、演奏された音楽に対する楽譜の位置をより頑健に推定することが可能になる。

また、上述した実施形態では、音階が変化する楽譜の位置に対応する時刻を含む観測区間における音響信号の自己相関に基づいて拍間隔を定め、定めた拍間隔に基づいて探索する楽譜の位置を更新する。
これにより、楽譜情報における音階が変化する楽譜の位置と、音響信号の振幅が変化する時刻が対応付けられる。そのため、周波数特性の周期性を表す拍間隔を正確に検知でき、ひいては演奏された音楽に対する楽譜の位置をより正確に推定することができる。

上述では、楽譜位置探索部１３１が、粒子フィルタリング法を用いて楽譜位置を探索する例について説明したが、上述した実施形態ではこれには限られない。上述した実施形態では、楽譜位置毎の楽譜情報と音響信号との関連度を表す重み情報に基づいて楽譜位置を定めることができる処理であれば、いかなる処理、例えば、ある探索区間内を候補となる楽譜位置を走査しながら重み情報を算出して、算出した重み情報に基づいて楽譜位置を探索してもよい。

なお、上述した実施形態における楽譜位置推定装置１、２又は３の一部、例えば、音響特徴量生成部１０２、関連度算出部１２１、楽譜位置探索部１３１、及び演奏情報出力部１４１をコンピュータで実現するようにしても良い。その場合、この制御機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現しても良い。なお、ここでいう「コンピュータシステム」とは、楽譜位置推定装置１、２又は３に内蔵されたコンピュータシステムであって、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含んでも良い。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。
また、上述した実施形態における楽譜位置推定装置１、２又は３の一部、または全部を、ＬＳＩ（Ｌａｒｇｅ Ｓｃａｌｅ Ｉｎｔｅｇｒａｔｉｏｎ）等の集積回路として実現しても良い。楽譜位置推定装置１又は２の各機能ブロックは個別にプロセッサ化してもよいし、一部、または全部を集積してプロセッサ化しても良い。また、集積回路化の手法はＬＳＩに限らず専用回路、または汎用プロセッサで実現しても良い。また、半導体技術の進歩によりＬＳＩに代替する集積回路化の技術が出現した場合、当該技術による集積回路を用いても良い。

以上、図面を参照してこの発明の一実施形態について詳しく説明してきたが、具体的な構成は上述のものに限られることはなく、この発明の要旨を逸脱しない範囲内において様々な設計変更等をすることが可能である。

１、２…楽譜位置推定装置、１０１…音響信号入力部、１０２…音響特徴量生成部、
１０３…周波数特性分析部、１０４…相関算出部、３０５…クロマベクトル生成部
１１１…楽譜情報記憶部、１１２…楽譜情報入力部、
１２１…関連度算出部、１２２、２２２…周波数特性重み算出部、
１２２１…データ入力部、１２２２…第１変分パラメータ算出部、
１２２３…第２変分パラメータ算出部、１２２４…変分下限値算出部、
１２２５…収束判定部、１２２６…重み算出部、
２２２１…周波数特性量子化部、２２２２…周波数特性合成部、
２２２３…振幅期待値算出部、２２２４…振幅分布算出部、２２２５…重み算出部、
１２３…拍間隔重み算出部、１２４、３２４…粒子重み算出部、
３２５…クロマベクトル重み算出部
１３１…楽譜位置探索部、１３２…粒子遷移部、１３３…再標本化部、
１４１…演奏情報出力部、１４２…楽譜位置算出部、１４３…拍間隔算出部、
１４４…信頼度判定部、１４５…楽譜位置出力部、１４６…拍間隔出力部

Claims (5)

1. 入力された音響信号の周波数特性を分析して第１の周波数特性を算出する周波数特性分析部と、
   楽譜情報が表す楽譜の位置毎の音階に基づく周波数特性であって少なくとも１つの調波構造を含む第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記調波構造に含まれる調波成分毎の前記調波構造への寄与を表す変数の確率分布と、前記調波構造毎の前記第２の周波数特性への寄与を表す変数の確率分布に基づいて前記第１の周波数特性についての尤度を最大にするように、前記楽譜の位置毎に算出する関連度算出部と、
   前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する楽譜位置探索部と、
   を備えることを特徴とする楽譜位置推定装置。
2. 入力された音響信号の周波数特性を分析して第１の周波数特性を算出する周波数特性分析部と、
   楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記第２の周波数特性に含まれる周波数毎の振幅が表す度数が、対応する周波数における前記第１の周波数特性の振幅が表す度数だけ発生する確率分布に基づいて前記楽譜の位置毎に算出する関連度算出部と、
   前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する楽譜位置探索部と、
   を備えることを特徴とする楽譜位置推定装置。
3. 前記楽譜位置探索部は、
   前記音階が変化する楽譜の位置に対応する時刻を含む観測区間における前記音響信号の自己相関に基づいて拍間隔を定め、定めた拍間隔に基づいて探索する楽譜の位置を更新することを特徴とする請求項１または請求項２に記載の楽譜位置推定装置。
4. 楽譜位置推定装置における方法であって、
   前記楽譜位置推定装置は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する過程と、
   前記楽譜位置推定装置は、楽譜情報が表す楽譜の位置毎の音階に基づく周波数特性であって少なくとも１つの調波構造を含む第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記調波構造に含まれる調波成分毎の前記調波構造への寄与を表す変数の確率分布と、前記調波構造毎の前記第２の周波数特性への寄与を表す変数の確率分布に基づいて前記第１の周波数特性についての尤度を最大にするように、前記楽譜の位置毎に算出する過程と、
   前記楽譜位置推定装置は、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する過程と、
   を有することを特徴とする楽譜位置推定方法。
5. 楽譜位置推定装置における方法であって、
   前記楽譜位置推定装置は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する過程と、
   前記楽譜位置推定装置は、楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記第２の周波数特性に含まれる周波数毎の振幅が表す度数が、対応する周波数における前記第１の周波数特性の振幅が表す度数だけ発生する確率分布に基づいて前記楽譜の位置毎に算出する過程と、
   前記楽譜位置推定装置は、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する過程と、
   を有することを特徴とする楽譜位置推定方法。