(57)【特許請求の範囲】
【請求項１】
入力された音響信号の周波数特性を分析して第１の周波数特性を算出する周波数特性分析部と、
楽譜情報が表す楽譜の位置毎の音階に基づく周波数特性であって少なくとも１つの調波構造を含む第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記調波構造に含まれる調波成分毎の前記調波構造への寄与を表す変数の確率分布と、前記調波構造毎の前記第２の周波数特性への寄与を表す変数の確率分布に基づいて前記第１の周波数特性についての尤度を最大にするように、前記楽譜の位置毎に算出する関連度算出部と、
前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する楽譜位置探索部と、
を備えることを特徴とする楽譜位置推定装置。
【請求項２】
入力された音響信号の周波数特性を分析して第１の周波数特性を算出する周波数特性分析部と、
楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記第２の周波数特性に含まれる周波数毎の振幅が表す度数が、対応する周波数における前記第１の周波数特性の振幅が表す度数だけ発生する確率分布に基づいて前記楽譜の位置毎に算出する関連度算出部と、
前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する楽譜位置探索部と、
を備えることを特徴とする楽譜位置推定装置。
【請求項３】
前記楽譜位置探索部は、  
前記音階が変化する楽譜の位置に対応する時刻を含む観測区間における前記音響信号の自己相関に基づいて拍間隔を定め、定めた拍間隔に基づいて探索する楽譜の位置を更新することを特徴とする請求項１または請求項２に記載の楽譜位置推定装置。
【請求項４】
楽譜位置推定装置における方法であって、
前記楽譜位置推定装置は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する過程と、
前記楽譜位置推定装置は、楽譜情報が表す楽譜の位置毎の音階に基づく周波数特性であって少なくとも１つの調波構造を含む第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記調波構造に含まれる調波成分毎の前記調波構造への寄与を表す変数の確率分布と、前記調波構造毎の前記第２の周波数特性への寄与を表す変数の確率分布に基づいて前記第１の周波数特性についての尤度を最大にするように、前記楽譜の位置毎に算出する過程と、
前記楽譜位置推定装置は、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する過程と、
を有することを特徴とする楽譜位置推定方法。
【請求項５】
楽譜位置推定装置における方法であって、
前記楽譜位置推定装置は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する過程と、
前記楽譜位置推定装置は、楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記第２の周波数特性に含まれる周波数毎の振幅が表す度数が、対応する周波数における前記第１の周波数特性の振幅が表す度数だけ発生する確率分布に基づいて前記楽譜の位置毎に算出する過程と、
前記楽譜位置推定装置は、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する過程と、
を有することを特徴とする楽譜位置推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、楽譜位置推定装置、及び楽譜位置推定方法に関する。
【背景技術】
【０００２】
従来から、楽曲を書き表した楽譜において演奏されている楽曲の位置（演奏位置）を推定する楽譜位置推定技術が開発されている。楽譜位置推定技術は、自動伴奏、自動合奏、楽譜への演奏位置の重畳表示を実現するための応用が期待されている。ここで、演奏されている楽曲の音響信号の周波数特性と楽譜情報に基づいて合成した楽曲の周波数特性を照合し、演奏されている楽曲について楽譜の位置を実時間で推定する技術が提案されている
【０００３】
例えば、特許文献１に記載の楽譜位置推定技術では、楽譜情報に基づく楽曲の周波数特性を調和混合正規分布モデル（ｈａｒｍｏｎｉｃ  Ｇａｕｓｓｉａｎ  Ｍｉｘｔｕｒｅ  Ｍｏｄｅｌ；ｈａｒｍｏｎｉｃ  ＧＭＭ、又は単にＧＭＭともいう。）に基づくテンプレートを用いて合成し、音響信号の周波数特性と照合していた。このテンプレートでは、高調波の振幅を、予め定めた実数のべき乗に定められていた。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０１１−１８０５９０号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかし、特許文献１記載の楽譜位置推定技術では、合成した楽曲の周波数特性と音響信号の周波数特性が合致せず、演奏位置を頑健に推定できないことがあった。
本発明は上記の点に鑑みてなされたものであり、演奏された楽曲の楽譜上の位置をより頑健に推定できる楽譜位置推定装置及び楽譜位置推定方法を提供することを課題としている。
【課題を解決するための手段】
【０００６】
（１）本発明は上記の課題を解決するためになされたものであり、本発明の一態様は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する周波数特性分析部と、楽譜情報が表す楽譜の位置毎の音階に基づく周波数特性であって少なくとも１つの調波構造を含む第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記調波構造に含まれる調波成分毎の前記調波構造への寄与を表す変数の確率分布と、前記調波構造毎の前記第２の周波数特性への寄与を表す変数の確率分布に基づいて前記第１の周波数特性についての尤度を最大にするように、前記楽譜の位置毎に算出する関連度算出部と、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する楽譜位置探索部と、を備えることを特徴とする楽譜位置推定装置である。
【０００７】
（２）本発明のその他の態様は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する周波数特性分析部と、楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記第２の周波数特性に含まれる周波数毎の振幅が表す度数が、対応する周波数における前記第１の周波数特性の振幅が表す度数だけ発生する確率分布に基づいて前記楽譜の位置毎に算出する関連度算出部と、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する楽譜位置探索部と、を備えることを特徴とする楽譜位置推定装置である。
【０００８】
（３）本発明のその他の態様は、前記楽譜位置探索部は、前記音階が変化する楽譜の位置に対応する時刻を含む観測区間における前記音響信号の自己相関に基づいて拍間隔を定め、定めた拍間隔に基づいて探索する楽譜の位置を更新することを特徴とする。
【０００９】
（４）本発明のその他の態様は、楽譜位置推定装置における方法であって、前記楽譜位置推定装置は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する過程と、前記楽譜位置推定装置は、楽譜情報が表す楽譜の位置毎の音階に基づく周波数特性であって少なくとも１つの調波構造を含む第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記調波構造に含まれる調波成分毎の前記調波構造への寄与を表す変数の確率分布と、前記調波構造毎の前記第２の周波数特性への寄与を表す変数の確率分布に基づいて前記第１の周波数特性についての尤度を最大にするように、前記楽譜の位置毎に算出する過程と、前記楽譜位置推定装置は、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する過程と、を有することを特徴とする楽譜位置推定方法。
【００１０】
（５）本発明のその他の態様は、楽譜位置推定装置における方法であって、前記楽譜位置推定装置は、入力された音響信号の周波数特性を分析して第１の周波数特性を算出する過程と、前記楽譜位置推定装置は、楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、前記第２の周波数特性に含まれる周波数毎の振幅が表す度数が、対応する周波数における前記第１の周波数特性の振幅が表す度数だけ発生する確率分布に基づいて前記楽譜の位置毎に算出する過程と、前記楽譜位置推定装置は、前記重み値に基づいて前記音響信号に対応する楽譜の位置を探索する過程と、を有することを特徴とする楽譜位置推定方法である。
【発明の効果】
【００１１】
（１）、（２）、（４）、（５）に記載した態様によれば、音響信号に基づく第１の周波数特性と、楽譜情報における音階に基づく第２の周波数特性が完全に合致していなくとも、両者の関連性を検出することができる。そのため、演奏された音楽に対する楽譜の位置をより頑健に推定することが可能になる。
（１）、（４）に記載した態様によれば、楽譜情報に基づく第２の周波数成分に含まれる調波構造の周波数特性の変化、調波構造間の結合特性の変化に柔軟に対応できるので、音響信号に基づく第１の周波数特性との関連性をより的確に推定することができる。そのため、演奏された音楽に対する楽譜の位置をより正確に推定することが可能になる。
（２）、（５）に記載した態様によれば、音響信号に基づく第１の周波数特性の周波数毎の振幅と楽譜情報における音階に基づく第２の周波数特性の周波数毎の振幅が合致する度合いを鋭敏に検出することができる。そのため、演奏された音楽に対する楽譜の位置をより正確に推定することが可能になる。
（３）に記載した態様によれば、楽譜情報における音階が変化する楽譜の位置と、音響信号の振幅が変化する時刻が対応付けられる。そのため、周波数特性の周期性を表す拍間隔を正確に検知でき、ひいては演奏された音楽に対する楽譜の位置をより正確に推定することができる。
【図面の簡単な説明】
【００１２】
【図１】本発明の第１の実施形態に係る楽譜位置推定装置の構成を示す概略図である。
【図２】粒子フィルタリング法の処理を表す概念図である。
【図３】本実施形態に係る楽譜位置推定処理を表すフローチャートである。
【図４】楽譜情報と音響信号の関係の一例を表す概念図である。
【図５】本実施形態に係る粒子遷移処理を表すフローチャートである。
【図６】本実施形態に係る音源モデルの一例を表す概念図である。
【図７】本実施形態に係る調波構造の一例を表す図である。
【図８】本実施形態に係るヒストグラムの一例を表す図である。
【図９】本実施形態に係る周波数特性重み算出部の一構成例を表す概略図である。
【図１０】本実施形態における周波数特性重み算出処理の一例を表すフローチャートである。
【図１１】本発明の第２の実施形態に係る楽譜位置推定装置の構成を示す概略図である。
【図１２】本実施形態に係る周波数特性重み算出部の構成を表す概略図である。
【図１３】本実施形態における周波数特性重み算出処理を表すフローチャートである。
【図１４】本発明の第３の実施形態に係る楽譜位置推定装置の構成を示す概略図である。
【図１５】本実施形態に係る楽譜位置推定処理を表すフローチャートである。
【図１６】楽譜位置の推定誤差の一例を表す図である。
【図１７】楽譜位置の推定誤差のその他の例を表す図である。
【図１８】楽譜位置の推定誤差のその他の例を表す図である。
【図１９】音響信号と楽譜情報との関連性の一例を表す図である。
【発明を実施するための形態】
【００１３】
（第１の実施形態）
以下、図面を参照しながら本発明の第１の実施形態について説明する。
図１は、本実施形態に係る楽譜位置推定装置１の構成を示す概略図である。
楽譜位置推定装置１は、音響信号入力部１０１、音響特徴量生成部１０２、楽譜情報記憶部１１１、楽譜情報入力部１１２、関連度算出部１２１、楽譜位置探索部１３１、及び演奏情報出力部１４１を含んで構成される。
【００１４】
音響信号入力部１０１は、演奏された音楽を表す音波をディジタル音響信号に変換し、変換したディジタル音響信号を音響特徴量生成部１０２に出力する。音響信号入力部１０１は、例えばマイクロホンとアナログ・ディジタル（Ａｎａｌｏｇ−ｔｏ−Ｄｉｇｉｔａｌ；Ａ／Ｄ）変換部（図示せず）とを含んで構成される。このマイクロホンは、人間が聴取することができる周波数帯域（例えば、２０Ｈｚ−２０ｋＨｚ）の音波をアナログ音響信号に変換し、変換したアナログ音響信号をＡ／Ｄ変換部に出力する。Ａ／Ｄ変換部は、マイクロホンから入力されたアナログ音響信号をディジタル音響信号に変換し、変換したディジタル音響信号を音響特徴量生成部１０２に出力する。ここで、Ａ／Ｄ変換部は、入力されたアナログ音響信号を、例えば、サンプリング周波数４４．１ｋＨｚ、振幅を１６ビットの２進（バイナリ；ｂｉｎａｒｙ）データにパルス符号変調（Ｐｕｌｓｅ  Ｃｏｄｅ  Ｍｏｄｕｌａｔｉｏｎ；ＰＣＭ）して、サンプル時刻毎の量子化された振幅を表すディジタル音響信号に変換する。以下の説明では、入力されたディジタル音響信号を入力信号と呼ぶことがある。
【００１５】
音響特徴量生成部１０２は、音響信号入力部１０１から入力されたディジタル音響信号から、その物理的な特徴を表す音響特徴量を生成し、生成した音響特徴量を関連度算出部１２１及び楽譜位置探索部１３１に出力する。音響特徴量生成部１０２は、音響特徴量として、例えば、振幅周波数特性として周波数毎の振幅を表す音響スペクトログラム、自己相関、音響スペクトログラムの時間差分に基づく距離値を算出する。音響特徴量生成部１０２は、周波数特性分析部１０３及び相関算出部１０４を含んで構成される。
【００１６】
周波数特性分析部１０３は、音響信号入力部１０１から入力されたディジタル音響信号を時間領域信号から周波数領域信号に変換する。ここで、周波数特性分析部１０３は、例えば、予め定められた個数（例えば、２０４８個）サンプルからなるフレーム（以下、音響フレームとも呼ぶ）ごとにディジタル音響信号を高速フーリエ変換（Ｆａｓｔ  Ｆｏｕｒｉｅｒ  Ｔｒａｎｓｆｏｒｍ；ＦＦＴ）を行って周波数領域信号に変換する。周波数特性分析部１０３は、予め設定された音響フレーム時刻毎（音響フレーム間隔Δｔは、例えば１０ｍｓ）に音響フレームを移動させる。音響フレームを移動させることによって、所定数の新たな信号のサンプル（例えば、サンプリング周波数が４４．１ｋＨｚの場合、４４１サンプル）が音響フレームに含まれ、最も時刻が経過した同数の信号のサンプルが音響フレームから除外される。
【００１７】
周波数特性分析部１０３は、時刻ｔにおける周波数領域信号に基づいてスペクトログラム（以下、音響スペクトログラムと呼ぶ）Ａｆ，ｔを算出する。音響スペクトログラムＡｆ，ｔは、周波数領域信号の周波数ｆ毎の振幅（絶対値）Ａｆ，ｔである。以下の説明では、音響スペクトログラムＡｆ，ｔを、入力音響振幅と呼ぶことがある。周波数特性分析部１０３は、算出した音響スペクトログラムＡｆ，ｔを相関算出部１０４、及び関連度算出部１２１の周波数特性重み算出部１２２に出力する。
【００１８】
相関算出部１０４は、周波数特性分析部１０３から入力された音響スペクトログラムＡｆ，ｔに基づいて、拍間隔ｂについての自己相関を算出する。自己相関は、時刻ｔにおける音響スペクトログラムＡｆ，ｔと、拍間隔ｂだけ過去の時刻ｔ−ｂにおける音響スペクトログラムＡｆ，ｔ−ｂの類似性を表す指標値である。拍間隔ｂとは、演奏される楽曲の時間単位である拍と、その拍に後続する拍との間の時間間隔を表す。相関算出部１０４は、自己相関として、例えば式（１）を用いてＲ（ｂ，｛Ａｋ｝）を算出する。
【００１９】
【数１】
0005924968.tif 000002
【００２０】
式（１）において、｛Ａｋ｝は、観測区間Ｔｋにおける現フィルタリングステップｋまでの音響スペクトログラムＡｆ，ｔの集合（セット）を表す。Ｔｋは、観測区間｛τ｜ｋΔＴ−Ｗ＜τ≦ｋΔＴ｝を表す。τは、時刻を表す。ΔＴは、例えば後述する楽譜位置探索におけるフィルタリング間隔（例えば、１秒）である。Ｗは、観測区間長（例えば、２．５秒）である。フィルタリングステップｋとは、後述する粒子フィルタ法（パーティクルフィルタリング、ｐａｒｔｉｃｌｅ  ｆｉｌｔｅｒｉｎｇ）を応用した楽譜位置推定における処理の繰り返しを表す整数である。
相関算出部１０４は、算出した自己相関を関連度算出部１２１の拍間隔重み算出部１２３及び楽譜位置探索部１３１の粒子遷移部１３２に出力する。相関算出部１０４は、音響スペクトログラムＡｆ，ｔを粒子遷移部１３２に出力する。
【００２１】
相関算出部１０４は、自己相関Ｒ（ｂ，｛Ａｋ｝）の代わりに、自己相関Ｒ（ｂ，｛Ξｋ’｝）を算出してもよい。｛Ξｋ’｝は、距離値Ξｆ，ｔを聴感補正した距離値Ξｆ，ｔ’の集合（セット）である。距離値Ξｆ，ｔは、現時刻ｔにおける音響スペクトログラムＡｆ，ｔと直前の音響フレームの時刻ｔ−Δｔにおける音響スペクトログラムＡｆ，ｔ―Δｔとの差分の絶対値である。
相関算出部１０４は、距離値として、例えば式（２）で表されるΞｆ，ｔを算出してもよい。
【００２２】
【数２】
0005924968.tif 000003
【００２３】
式（２）において、Δφｆ，ｔは、式（３）で表される。
【００２４】
【数３】
0005924968.tif 000004
【００２５】
式（３）において、φｆ，ｔ等は、対応する音響スペクトログラムＡｆ，ｔの絶対位相（アンラップト位相［ｕｎｗｒａｐｐｅｄ  ｐｈａｓｅ］）である。
相関算出部１０４は、線形領域で表された周波数ｆをメル尺度で表された周波数ｆｍｅｌに変換する。線形領域で表された周波数ｆと、メル尺度で表された周波数ｆｍｅｌとの関係は、式（４）で表される。
【００２６】
【数４】
0005924968.tif 000005
【００２７】
相関算出部１０４は、フレーム毎の周波数の次元数（例えば、１０２４）をＮｍ（Ｎｍは１よりも大きい整数、例えば、６４である。）に減少させる。そのため、相関算出部１０４は、距離値Ξｆ，ｔを中心周波数ｆｍｍｅｌ毎にフィルタリングして、聴感補正した距離値Ξｆ，ｔ’を算出する。中心周波数ｆｍｍｅｌは、式（５）に示されるように、０からナイキスト周波数ｆＮｙｑｍｅｌまでの間を等分する周波数である。
【００２８】
【数５】
0005924968.tif 000006
【００２９】
式（５）において、ｍは、０より大きくＮｍと等しいか、Ｎｍよりも小さい整数である。フィルタリングに用いる窓関数Ψｍ（ｆｍｅｌ）は、例えば、式（６）で表される。
【００３０】
【数６】
0005924968.tif 000007
【００３１】
式（６）は、窓関数Ψｍ（ｆｍｅｌ）の周波数特性は、各中心周波数ｆｍｍｅｌにおいてピーク１をとり、隣接する中心周波数ｆｍ−１ｍｅｌ，ｆｍ＋１ｍｅｌにおいてゼロとなる三角形を表す。但し、ｍ＝Ｎｍの場合、窓関数ψＮｍ（ｆｍｅｌ）は、周波数ｆＮｍ−１ｍｅｌ≦ｆｍｅｌ＜ｆＮｍｍｅｌの区間のみ定められる。周波数ｆＮｍｍｅｌは、ナイキスト周波数ｆＮｙｑｍｅｌであるため、この周波数を越える周波数成分が存在しないからである。
相関算出部１０４は、例えば式（７）を用いて距離値Ξｆ，ｔをフィルタリングする。
【００３２】
【数７】
0005924968.tif 000008
【００３３】
相関算出部１０４は、聴感補正した距離値Ξｆ，ｔ’に基づいて、例えば式（８）を用いて自己相関Ｒ（ｂ，｛Ξｋ’｝）を算出する。
【００３４】
【数８】
0005924968.tif 000009
【００３５】
距離値Ξｆ，ｔ、又はΞｆ，ｔ’は、音響スペクトログラムＡｆ，ｔの拍間隔ｂに対する周期性を表す指標値である。例えば、定常的な音では、距離値Ξｆ，ｔ、又はΞｆ，ｔ’は、ゼロ又はゼロに近似する。そのため、自己相関Ｒ（ｂ，｛Ξｋ’｝）を用いることで、立ち上がり部分（オンセット）のように、入力信号の振幅の変動が著しい部分に基づいて、振幅周波数特性の周期性を検知することが容易になる。
相関算出部１０４は、算出した距離値Ξｆ，ｔ及び自己相関Ｒ（ｂ，｛Ξｋ’｝）を粒子遷移部１３２に出力する。
【００３６】
楽譜情報記憶部１１１は、楽譜情報を楽曲毎に予め記憶されている。楽譜情報は、楽譜フレーム（ｓｃｏｒｅ  ｆｒａｍｅ）ｐ（ｐは、整数）毎に、楽曲を構成する音階毎の基本周波数μｐｌを要素とする基本周波数ベクトル{μｐ}＝［μｐ１，μｐ２，…，μｐＬｐ］Ｔで表される。ここで、Ｌｐは、楽譜フレームｐにおける音階の数を表す。ｌは０より大きく、Ｌｐと等しいかＬｐよりも小さい整数である。Ｔは、ベクトル又は行列の転置を表す。つまり、基本周波数ベクトルμｐは、同時に演奏される音階の周波数の組、即ちコード（和音、ｃｈｏｒｄ）を表す情報である。但し、本実施形態では、コードは、１個又は１個より多い数の音階を含んでいればよい。例えば、ある楽譜フレームｐ１において音階Ａ４、Ｃ４、及びＥ４を含むコードを表す基本周波数ベクトルμｐ１は、要素数が３個のベクトル［４４０，５２３，６５９］となる。ここで、音階Ａ４、Ｃ４、Ｅ４の周波数は、それぞれ４４０、５２３、６５９Ｈｚである。
【００３７】
楽譜情報は、楽曲毎の演奏テンポｂｓ’を表すテンポ情報を含んでいでもよい。テンポ情報は、例えばテンポの逆数である拍間隔ｂｓ’で表される。
楽譜情報は、コードを示すコード情報ｎｐｌを要素とするコードベクトル{ｎｐ}＝［ｎｐ１，ｎｐ２，…，ｎｐＬｐ］Ｔを、基本周波数ベクトル{μｐ}と対応付けて含んでいてもよい。音階Ａ４、Ｃ４、及びＥ４を要素値として含むコードベクトルは、［Ａ４，Ｃ４，Ｅ４］となる。
楽譜フレームは、１つの楽曲の楽譜情報において楽曲を時間分割する単位時間である。例えば、１拍、即ち四分音符の長さが１２フレームである場合、楽譜情報の時間分解能（ｒｅｓｏｌｕｔｉｏｎ）は十六分音符の三分の一である。
【００３８】
楽譜情報入力部１１２は、楽譜情報記憶部１１１から、処理の対象とする楽曲の楽譜情報を読み出す。楽譜情報入力部１１２は、読み出した楽譜情報を関連度算出部１２１の周波数特性重み算出部１２２、及び楽譜位置探索部１３１の粒子遷移部１３２に出力する。
【００３９】
関連度算出部１２１は、音響特徴量生成部１０２から入力された音響特徴量と楽譜情報入力された楽譜情報が表す楽譜特徴量との関連度を表す重み値を楽譜位置毎に算出する。楽譜特徴量は、例えば、楽譜情報に基づいて演奏されたと仮定された場合における音響信号の振幅周波数特性である。楽譜位置とは、楽譜の冒頭からの拍数や楽譜フレーム数であり、演奏されている楽曲の楽譜上の拍数や楽譜フレーム数を指すこともある。関連度算出部１２１は、算出した楽譜位置毎の重み値を表す重み情報を楽譜位置探索部１３１及び演奏情報出力部１４１に出力する。
関連度算出部１２１は、例えば、周波数特性重み算出部１２２、拍間隔重み算出部１２３、及び粒子重み算出部１２４を含んで構成される。粒子重み算出部１２４は、上述の重み情報として、粒子重み情報を含んだ粒子情報を楽譜位置探索部１３１及び演奏情報出力部１４１に出力する。
【００４０】
周波数特性重み算出部１２２は、周波数特性分析部１０３から入力された音響スペクトログラムと楽譜情報入力部１１２から入力された楽譜情報に含まれる基本周波数ベクトルが表す周波数特性との楽譜位置毎の関連度（例えば、類似度）を表す周波数特性重み値を算出する。周波数特性重み算出部１２２は、算出した周波数特性重み値を表す周波数特性重み情報を粒子重み算出部１２４に出力する。周波数特性重み算出部１２２のより詳細な構成については後述する。
拍間隔重み算出部１２３は、相関算出部１０４から入力された自己相関に基づいて拍間隔毎に音響信号の振幅周波数特性の関連度を表す拍間隔重み値を算出する。拍間隔重み算出部１２３は、拍間隔重み値を表す拍間隔重み情報を粒子重み算出部１２４に出力する。拍間隔重み算出部１２３のより詳細な構成については後述する。
【００４１】
粒子重み算出部１２４は、後述する粒子遷移部１３２から入力された粒子ｉ毎の粒子情報、粒子分布値、状態遷移確率、周波数特性重み算出部１２２から入力された周波数特性重み情報と拍間隔重み算出部１２３から入力された拍間隔重み情報に基づいて、楽譜位置及び拍間隔毎の重み情報を粒子ｉ毎に生成する。
粒子とは、楽譜位置を表す楽譜位置情報、拍間隔を表す拍間隔情報及び入力音響信号に対する観測値である重み情報の組を表す情報の単位である。粒子重み算出部１２４が生成した重み情報を粒子重み情報と呼び、楽譜位置情報、拍間隔情報及び粒子重み情報の組を粒子情報と呼ぶ。粒子情報は、後述するように粒子フィルタリング法により、楽譜位置、拍間隔の推定値を定めるために用いられる。
粒子フィルタリング法を用いる場合、粒子重み算出部１２４は、例えば、式（９）を用いて粒子重み情報として粒子ｉ毎の重み値（以下、粒子重み値と呼ぶ）ｗｉ，ｋを算出する。
【００４２】
【数９】
0005924968.tif 000010
【００４３】
式（９）において、ｗｓｐｉ，ｋは、周波数特性重み値を表し、ｗｔｍｉ，ｋは、拍間隔重み値を表す。ｑ（ｐｉｋ，ｂｉｋ｜｛Ａｋ｝，ｂｓ’，ｏｐ）、ｑ（ｐｉｋ，ｂｉｋ｜ｐｉｋ−１，ｂｉｋ−１）は、それぞれ粒子ｉの粒子分布値、状態遷移確率を表す。粒子重み算出部１２４は、算出した粒子重み値ｗｉ，ｋを表す粒子重み情報を生成する。
粒子重み算出部１２４は、入力された粒子ｉ毎の粒子情報に含まれる粒子重み情報を、生成した粒子重み情報に置き換えることで粒子情報を更新する。
粒子重み算出部１２４は、更新した粒子情報を楽譜位置探索部１３１の再標本化部１３３、演奏情報出力部１４１の楽譜位置算出部１４２及び拍間隔算出部１４３に出力する。
【００４４】
楽譜位置探索部１３１は、粒子重み算出部１２４から入力された楽譜位置毎の重み情報に基づいて楽譜位置を探索する。楽譜位置探索部１３１は、例えば、重み情報が表す重み値が音響信号と楽譜情報との関連性が最も高い楽譜位置を探索する。楽譜位置探索部１３１は、探索した楽譜位置を表す楽譜位置情報を関連度算出部１２１に出力する。
【００４５】
楽譜位置探索部１３１は、例えば、粒子フィルタリング法を用いて楽譜位置を探索する。楽譜位置探索部１３１は粒子遷移部１３２及び再標本化部１３３を含んで構成される。ここで、粒子フィルタリング法について説明する。
図２は、粒子フィルタリング法の処理を表す概念図である。
粒子フィルタリング法は、Ｉ．遷移、ＩＩ．重み算出、ＩＩＩ．位置推定・再標本化、の各過程を有する。図２は、Ｉ．遷移、ＩＩ．重み算出、ＩＩＩ．位置推定・再標本化、の各過程における処理の概念を左から右へ順に表す。
【００４６】
図２のＩ．遷移において、横軸が時刻、即ち楽譜位置を表し、縦軸が拍間隔を表す。塗りつぶした円は、それぞれ粒子を表す。各矢印は、起点にある粒子が終点に遷移することを表す。ここで、粒子遷移部１３２が、各時刻、各拍間隔にわたって分布している複数の粒子から各粒子を１つずつ抽出し、抽出した粒子毎の楽譜位置を、その拍間隔に基づいて更新することを表す。
【００４７】
ＩＩ．重み算出では、上段に音響信号、中段に楽譜情報、下段に粒子をそれぞれ表す。各段ともに横軸は時刻、即ち楽譜位置を表す。中段において縦軸は音階を表す。上段と中段にそれぞれ表された四角形の幅は、観測区間を表す。上段の四角形で囲まれる音響信号が、中段の四角形で囲まれる楽譜情報に対して観測対象となる音響信号を表す。下段に表される塗りつぶしの円は、粒子を表す。円の半径は、観測対象となる音響信号と対応する楽譜情報に基づく粒子重みの大きさを表す。即ち、関連度算出部１２１は、粒子毎の楽譜位置及び拍間隔の妥当性を表す粒子重みを、観測区間内の音響信号の音響特徴量と楽譜情報に基づいて算出する。
【００４８】
ＩＩＩ．位置推定・再標本化では、上段に楽譜情報、中段に再標本化前の粒子、下段に再標本化後の粒子を表す。各段ともに横軸は時刻、即ち楽譜位置を表す。中段において、下向きの矢印は粒子の分布に基づいて推定された推定楽譜位置を表す。即ち、後述する楽譜位置算出部１４２が、粒子毎の楽譜位置を各自の粒子重みで重み付け加算して楽譜位置を推定することを表す。
中段の粒子を起点とし、下段の粒子を終点とする矢印は、中段の粒子を再標本化することを表す。濃く塗りつぶした粒子を終点とする太線の矢印は、起点の粒子の粒子重みが大きいために複数の粒子に分割・複製することを表す。×印を終点とする破線の矢印は、起点の粒子の粒子重みが小さいために、その粒子を棄却することを表す。即ち、再標本化部１３３が、粒子重みが大きい粒子を分割・複製し、粒子重みが小さい粒子ほど優先して棄却することを表す。これにより、楽譜位置及び拍間隔が妥当な粒子が残される。その後、Ｉ．遷移に戻り、Ｉ−ＩＩＩの過程を繰り返す。粒子遷移部１３２、再標本化部１３３の構成もしくは動作については後述する。
【００４９】
図１に戻り、演奏情報出力部１４１は、関連度算出部１２１から入力された関連度情報として粒子情報に基づいて楽譜位置を算出し、算出した楽譜位置を表す楽譜位置情報を生成する。演奏情報出力部１４１は、生成した楽譜位置情報を含む演奏情報を楽譜位置推定装置１の外部に出力する。
演奏情報出力部１４１は、例えば、楽譜位置算出部１４２、拍間隔算出部１４３、信頼度判定部１４４、及び楽譜位置出力部１４５を含んで構成される。
【００５０】
楽譜位置算出部１４２は、粒子重み算出部１２４から入力された粒子ｉ毎の粒子情報から楽譜位置情報と粒子重み情報を抽出する。楽譜位置算出部１４２は、粒子ｉ毎に抽出した楽譜位置情報が表す楽譜位置ｐｉｋに基づいて推定楽譜位置＜ｐｋ＞を算出する。楽譜位置算出部１４２は、例えば、楽譜位置ｐｉｋの単純平均を推定楽譜位置＜ｐｋ＞としてもよいし、粒子重み情報が表す粒子重み値ｗｉ，ｋによる加重平均を推定楽譜位置＜ｐｋ＞としてもよい。楽譜位置算出部１４２が加重平均を算出する際、必ずしも全ての粒子ｉを考慮する必要はなく、一部の粒子のみを考慮してもよい。一部の粒子とは、例えば、粒子重み値ｗｉ，ｋが最も大きい粒子から粒子重み値ｗｉ，ｋの大きい順に予め定めた個数Ｉｖ（例えば、粒子数Ｉの２％）までの粒子である。これにより、楽譜位置ｐｉｋの信頼性が低いことを表す粒子重み値ｗｉ，ｋがより小さい粒子が無視されるため、楽譜位置の推定精度が向上する。楽譜位置算出部１４２は、算出した推定楽譜位置＜ｐｋ＞を表す推定楽譜位置情報を信頼度判定部１４４に出力する。
楽譜位置算出部１４２は、一部の粒子を考慮して推定楽譜位置＜ｐｋ＞を算出する際、一部の粒子についての粒子重み値の総和の他に、さらに全部の粒子についての粒子重み値の総和を算出する。楽譜位置算出部１４２は、一部の粒子についての粒子重み値の総和、全部の粒子についての粒子重み値の総和を信頼度判定部１４４に出力する。
【００５１】
拍間隔算出部１４３は、粒子重み算出部１２４から入力された粒子ｉ毎の粒子情報から拍間隔情報と粒子重み情報を抽出する。拍間隔算出部１４３は、粒子ｉ毎に抽出した拍間隔情報が表す拍間隔ｂｉｋに基づいて推定拍間隔＜ｂｋ＞を算出する。拍間隔算出部１４３は、例えば、拍間隔ｂｉｋの単純平均を推定拍間隔＜ｂｋ＞としてもよいし、粒子重み情報が表す粒子重み値ｗｉ，ｋによる加重平均を推定拍間隔＜ｂｋ＞としてもよい。拍間隔算出部１４３が加重平均を算出する際、必ずしも全ての粒子ｉを考慮する必要はなく、楽譜位置算出部と同様に一部の粒子のみを考慮してもよい。これにより、推定拍間隔＜ｂｋ＞の信頼性が低いことを表す粒子重み値がより小さい粒子が無視され、楽譜位置の推定精度が向上する。拍間隔算出部１４３は、算出した推定拍間隔＜ｂｋ＞を表す推定拍間隔情報を拍間隔出力部１４６に出力する。
【００５２】
信頼度判定部１４４は、楽譜位置算出部１４２から入力された一部の粒子についての粒子重み値の総和、全部の粒子についての粒子重み値の総和に基づいて信頼度係数ｖｋを算出する。信頼度係数ｖｋは、例えば、式（１０）に示されるように、一部の粒子についての粒子重み値の総和を全部の粒子についての粒子重み値の総和で除算した値である。
【００５３】
【数１０】
0005924968.tif 000011
【００５４】
即ち、信頼度係数ｖｋは、粒子重み値ｗｉ，ｋが大きい粒子が一部の粒子に集中している度合いであり、推定した楽譜位置や拍間隔の信頼性を表す指標値である。
信頼度判定部１４４は、現在のフィルタリングステップｋにおける信頼度係数ｖｋから前フィルタリングステップｋ−１における信頼度係数ｖｋ−１の差分ｖｋ−ｖｋ−１を算出する。信頼度判定部１４４は、差分ｖｋ−ｖｋ−１が予め定めた閾値−ηｄｅｃ（例えば、−０．０７）よりも小さい場合、楽譜位置算出部１４２から入力された推定楽譜位置情報を楽譜位置出力部１４５へ出力することを停止する。この場合、信頼度係数ｖｋ−１が急激に低下するため、信頼性に欠ける楽譜位置情報の出力が回避される。
信頼度判定部１４４は、差分ｖｋ−ｖｋ−１が予め定めた閾値ηｉｎｃ（例えば、０．０８）よりも大きい場合、楽譜位置算出部１４２から入力された推定楽譜位置情報を楽譜位置出力部１４５へ出力することを再開する。
【００５５】
なお、信頼度判定部１４４は、上記の信頼度係数の代わりに推定楽譜位置情報が表す推定楽譜位置＜ｐｋ＞の推定誤差ｅ（ｔ）の絶対値｜ｅ（ｔ）｜もしくは二乗誤差ｅ２（ｔ）を信頼度係数として算出し、算出した信頼度が予め定めた閾値γｅを超えるか否かを判断してもよい。信頼度判定部１４４は、信頼度係数が閾値γｅを超えると判断したとき、楽譜位置算出部１４２から入力された推定楽譜位置情報の楽譜位置出力部１４５への出力を停止する。信頼度判定部１４４は、信頼度係数が閾値γｅを超えないと判断したとき、楽譜位置算出部１４２から入力された推定楽譜位置情報を楽譜位置出力部１４５へ出力する。  信頼度判定部１４４は、例えば式（１１）を用いて推定誤差ｅ（ｔ）を算出する。
【００５６】
【数１１】
0005924968.tif 000012
【００５７】
式（１１）において、ｓ（＜ｐｋ＞）は、推定楽譜位置＜ｐｋ＞の楽譜フレームに対応する時刻を表す。
【００５８】
楽譜位置出力部１４５は、信頼度判定部１４４から入力された推定楽譜位置情報を楽譜位置推定装置１の外部に出力する。
拍間隔出力部１４６は、拍間隔算出部１４３から入力された推定拍間隔情報を楽譜位置推定装置１の外部に出力する。拍間隔出力部１４６は、推定拍間隔情報自体の代わりに推定拍間隔情報が表す推定拍間隔＜ｂｋ＞の逆数である推定テンポを出力するようにしてもよい。
【００５９】
次に、本実施形態に係る楽譜位置推定処理について説明する。
図３は、本実施形態に係る楽譜位置推定処理を表すフローチャートである。
（ステップＳ１０１）楽譜位置推定装置１の各構成部は、処理に用いる変数やデータを初期設定する。例えば、粒子遷移部１３２は、粒子数Ｉ、粒子ｉ毎の拍間隔ｂｉｋ、楽譜位置ｐｉｋ、粒子分布値ｑ（ｐｉｋ，ｂｉｋ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）、状態遷移確率ｑ（ｐｉｋ，ｂｉｋ｜ｐｉｋ−１，ｂｉｋ−１）の初期値を定める。その後、ステップＳ１０２に進む。
【００６０】
（ステップＳ１０２）楽譜情報入力部１１２は、楽譜情報記憶部１１１から、処理の対象とする楽曲の楽譜情報を読み出す。楽譜情報入力部１１２は、読み出した楽譜情報を関連度算出部１２１の周波数特性重み算出部１２２、拍間隔重み算出部１２３、及び楽譜位置探索部１３１の粒子遷移部１３２に出力する。その後、ステップＳ１０３に進む。
【００６１】
（ステップＳ１０３）音響信号入力部１０１は、音波を表すアナログ音響信号をディジタル音響信号に変換し、変換したディジタル音響信号を周波数特性分析部１０３に出力する。その後、ステップＳ１０４に進む。
【００６２】
（ステップＳ１０４）周波数特性分析部１０３は、音響信号入力部１０１から入力されたディジタル音響信号を時間領域信号から周波数領域信号にフレーム時刻毎に変換する。周波数特性分析部１０３は、フレーム時刻毎の周波数領域信号に基づいて音響スペクトログラムＡｆ，ｔ（音響周波数特性）を算出する。周波数特性分析部１０３は、算出した音響スペクトログラムＡｆ，ｔを相関算出部１０４、及び関連度算出部１２１の周波数特性重み算出部１２２に出力する。その後、ステップＳ１０５に進む。
【００６３】
（ステップＳ１０５）相関算出部１０４は、周波数特性分析部１０３から入力された音響スペクトログラムＡｆ，ｔに基づいて、拍間隔ｂについての自己相関を、例えば式（１）を用いて算出する。相関算出部１０４は、算出した自己相関を関連度算出部１２１の拍間隔重み算出部１２３及び楽譜位置探索部１３１の粒子遷移部１３２に出力する。相関算出部１０４は、音響スペクトログラムＡｆ，ｔを粒子遷移部１３２に出力する。その後、ステップＳ１０６に進む。
【００６４】
（ステップＳ１０６）周波数特性重み算出部１２２は、周波数特性分析部１０３から入力された音響スペクトログラムと楽譜情報入力部１１２から入力された楽譜情報に含まれる基本周波数ベクトルが表す周波数特性との楽譜位置毎の関連度、即ち類似度を表す周波数特性重み情報を生成する。周波数特性重み算出部１２２は、生成した周波数特性重み情報を粒子重み算出部１２４に出力する。その後、ステップＳ１０７に進む。
【００６５】
（ステップＳ１０７）拍間隔重み算出部１２３は、相関算出部１０４から入力された自己相関に基づいて拍間隔毎に振幅特性の類似度を表す拍間隔重み情報を生成する。拍間隔重み算出部１２３は、生成した拍間隔重み情報を粒子重み算出部１２４に出力する。その後、ステップＳ１０８に進む。
【００６６】
（ステップＳ１０８）粒子重み算出部１２４は、周波数特性重み算出部１２２から入力された周波数特性重み情報、拍間隔重み算出部１２３から入力された拍間隔重み情報、粒子遷移部１３２から入力された粒子ｉ毎の粒子分布値、状態遷移確率、粒子情報に基づいて粒子重み情報を生成する。粒子重み算出部１２４は、例えば式（９）を用いて粒子重み値ｗｉ，ｋを算出し、算出した粒子重み値を表す粒子重み情報を生成する。粒子重み算出部１２４は、入力された粒子ｉ毎の粒子情報に含まれる粒子重み情報を、生成した粒子重み情報に置き換えることで粒子情報を更新する。
粒子重み算出部１２４は、更新した粒子情報を再標本化部１３３、楽譜位置算出部１４２及び拍間隔算出部１４３に出力する。その後、ステップＳ１０９に進む。
【００６７】
（ステップＳ１０９）楽譜位置算出部１４２は、粒子重み算出部１２４から入力された粒子ｉ毎の粒子情報から楽譜位置情報と粒子重み情報を抽出する。楽譜位置算出部１４２は、粒子ｉ毎に抽出した楽譜位置情報が表す楽譜位置ｐｉｋに基づいて推定楽譜位置＜ｐｋ＞を算出する。楽譜位置算出部１４２は、例えば、粒子重み情報が表す粒子重み値ｗｉ，ｋによる加重平均を推定楽譜位置＜ｐｋ＞とする。楽譜位置算出部１４２は、算出した推定楽譜位置＜ｐｋ＞を表す推定楽譜位置情報を信頼度判定部１４４に出力する。楽譜位置算出部１４２は、一部の粒子についての粒子重み値の総和と、全部の粒子についての粒子重み値の総和を算出する。楽譜位置算出部１４２は、算出した一部の粒子についての粒子重み値の総和、全部の粒子についての粒子重み値の総和を信頼度判定部１４４に出力する。その後、ステップＳ１１０に進む。
【００６８】
（ステップＳ１１０）拍間隔算出部１４３は、粒子重み算出部１２４から入力された粒子ｉ毎の粒子情報から拍間隔情報と粒子重み情報を抽出する。拍間隔算出部１４３は、例えば、粒子重み情報が表す粒子重み値ｗｉ，ｋによる加重平均を推定拍間隔＜ｂｋ＞として算出する。拍間隔算出部１４３は、算出した推定拍間隔＜ｂｋ＞を表す推定拍間隔情報を拍間隔出力部１４６に出力する。その後、ステップＳ１１１に進む。
【００６９】
（ステップＳ１１１）信頼度判定部１４４は、楽譜位置算出部１４２から入力された一部の粒子についての粒子重み値の総和、全部の粒子についての粒子重み値の総和に基づいて信頼度係数ｖｋを、例えば式（１０）を用いて算出する。
信頼度判定部１４４は、現在のフィルタリングステップｋにおける信頼度係数ｖｋから前フィルタリングステップｋ−１における信頼度係数ｖｋ−１の差分ｖｋ−ｖｋ−１を算出する。その後、ステップＳ１１２に進む。
【００７０】
（ステップＳ１１２）信頼度判定部１４４は、差分ｖｋ−ｖｋ−１が予め定めた閾値−ηｄｅｃ（閾値１）よりも小さいか否かを判断する。予め定めた閾値−ηｄｅｃよりも小さいと判断したとき（ステップＳ１１２  Ｙ）、ステップＳ１１３に進む。予め定めた閾値−ηｄｅｃよりも小さくないと判断したとき（ステップＳ１１２  Ｎ）、ステップＳ１１４に進む。
（ステップＳ１１３）信頼度判定部１４４は、楽譜位置算出部１４２から入力された推定楽譜位置情報＜ｐｋ＞の楽譜位置出力部１４５への出力を停止する。その後、ステップＳ１１４に進む。
【００７１】
（ステップＳ１１４）信頼度判定部１４４は、差分ｖｋ−ｖｋ−１が予め定めた閾値ηｉｎｃ（閾値２）よりも大きいか否かを判断する。予め定めた閾値ηｉｎｃよりも大きいと判断したとき（ステップＳ１１４  Ｙ）、ステップＳ１１５に進む。予め定めた閾値ηｉｎｃよりも大きくないと判断したとき（ステップＳ１１４  Ｎ）、ステップＳ１１６に進む。
（ステップＳ１１５）信頼度判定部１４４は、楽譜位置算出部１４２から入力された推定楽譜位置情報＜ｐｋ＞の楽譜位置出力部１４５への出力を再開する。その後、ステップＳ１１６に進む。
【００７２】
（ステップＳ１１６）再標本化部１３３は、入力された粒子情報に含まれる粒子重み情報が表す粒子重みに比例した確率で粒子数Ｉよりも多い予め定めた数の粒子を抽出する。再標本化部１３３は、増加した粒子数と同一の数だけ、粒子重み値ｗｉ，ｋがより小さい粒子を優先して抽出した粒子を棄却する。再標本化部１３３は、抽出された各粒子の楽譜位置情報及び拍間隔情報は、対応する抽出前の粒子の楽譜位置情報及び拍間隔情報と同一とする。再標本化部１３３は、棄却されずに残った粒子の粒子情報を粒子遷移部１３２に出力する。その後、ステップＳ１１７に進む。
【００７３】
（ステップＳ１１７）粒子遷移部１３２は、相関算出部１０４から距離値及び自己相関が、楽譜情報入力部１１２から入力された楽譜情報が、再標本化部１３３から粒子毎の粒子情報が入力される。粒子遷移部１３２は、入力された楽譜情報に基づいてオンセット情報を生成する。
粒子遷移部１３２は、楽譜情報に含まれるテンポ情報、自己相関、粒子情報に基づいて拍間隔分布値を算出する。粒子遷移部１３２は、算出した拍間隔分布値に比例する確率で粒子ｉの拍間隔ｂｉｋを定める。粒子遷移部１３２は、オンセット情報、距離値、及び抽出した拍間隔ｂｉｋに基づいて楽譜位置分布値を算出し、算出した楽譜位置分布値に比例する確率で粒子ｉの楽譜位置ｐｉｋを定める。
粒子遷移部１３２は、粒子ｉ毎の拍間隔分布値に、楽譜位置分布値を乗じて、粒子分布値ｑ（ｐｉｋ，ｂｉｋ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）を算出する。
【００７４】
粒子遷移部１３２は、粒子ｉ毎に前フィルタリングステップｋ−１における楽譜位置ｐｉｋ−１、拍間隔ｂｉｋ−１と、現フィルタリングステップｋで定めた楽譜位置ｐｉｋ、拍間隔ｂｉｋに基づいて、状態遷移確率ｑ（ｐｉｋ，ｂｉｋ｜ｐｉｋ−１，ｂｉｋ−１）を算出する。
粒子遷移部１３２は、粒子毎に粒子情報に含まれる拍間隔情報及び楽譜位置情報を、抽出した拍間隔を表す拍間隔情報及び楽譜位置を表す楽譜位置情報にそれぞれ置き換える。粒子遷移部１３２は、置き換えた粒子情報を関連度算出部１２１の周波数特性重み算出部１２２、拍間隔重み算出部１２３、粒子重み算出部１２４に出力する。
粒子遷移部１３２は、算出した粒子毎の粒子分布値と状態遷移確率を粒子重み算出部１２４に出力する。その後、処理を終了する。
【００７５】
ここで粒子遷移部１３２の構成について、より詳細に説明する。図１に戻り、粒子遷移部１３２は、相関算出部１０４から距離値及び自己相関が、楽譜情報入力部１１２から入力された楽譜情報が、再標本化部１３３から粒子毎の粒子情報が入力される。粒子遷移部１３２は、入力された楽譜情報からコードの変化（オンセット）を有するか否かを表すオンセット情報を生成する。コードの変化は、コードに含まれる音階の変化である。従って、コードの変化が楽譜位置では音響信号は、変化直後の音階の周波数成分が新たに発生し、振幅の立ち上がり（オンセット）として観測される。ここで、粒子遷移部１３２は、基本周波数ベクトル｛μｐ｝が、直前の楽譜フレームｐ−１における基本周波数ベクトル｛μｐ−１｝と異なる楽譜フレームｐがオンセットを有すると判断する。粒子遷移部１３２は、基本周波数ベクトル｛μｐ｝が、直前の楽譜フレームｐ−１における基本周波数ベクトル｛μｐ−１｝と同一ではない楽譜フレームｐはオンセットを有しないと判断する。粒子遷移部１３２は、楽譜フレームｐ毎にオンセットを有するか否かを表すオンセット情報を生成する。例えば、楽譜フレームｐがオンセットを有することを表すオンセット情報ｏｐは１であり、楽譜フレームｐがオンセットを有しないことを表すオンセット情報ｏｐは０である。
【００７６】
なお、本実施形態では、粒子遷移部１３２は、楽譜情報に対応するオンセット情報を楽譜情報入力部１１２から入力されるようにしてもよい。その場合、予め生成したオンセット情報を楽譜情報と対応付けて楽譜情報記憶部１１１に記憶させておき、楽譜情報入力部１１２が楽譜情報記憶部１１１から読み出したオンセット情報を粒子遷移部１３２に出力する。
【００７７】
粒子遷移部１３２は、生成したオンセット情報、入力された距離値、自己相関、及び楽譜情報に含まれるテンポ情報に基づいて粒子分布値を算出する。粒子遷移部１３２は、算出した粒子分布値が表す粒子分布で粒子を抽出することで粒子を遷移させる。粒子遷移部１３２は、遷移させた粒子毎の粒子情報を関連度算出部１２１の周波数特性重み算出部１２２及び拍間隔重み算出部１２４に出力する。
【００７８】
ここで、粒子遷移部１３２が、粒子分布値を算出する処理について説明する。
粒子遷移部１３２は、テンポ情報に基づいて拍間隔ｂに対する窓関数ψ（ｂ｜ｂｓ’）を算出する。窓関数ψ（ｂ｜ｂｓ’）は、例えば、テンポ情報が表すテンポ６０／ｂｓ’を中心とする予め定めた範囲内｜６０／ｂ−６０／ｂｓ’｜＜６０／ｂθでは１であり、その範囲外では０である。以下、６０／ｂθをテンポ窓長（ｔｅｍｐｏ  ｗｉｎｄｏｗ  ｌｅｎｇｔｈ）と呼ぶ。
【００７９】
粒子遷移部１３２は、自己相関Ｒ（ｂ，｛Ａｋ｝）に窓関数ψ（ｂ｜ｂｓ’）を乗じて拍間隔ｂに対する分布値ｑ（ｂ｜｛Ａｋ｝，ｂｓ’）を算出する。これにより、粒子遷移部１３２は、楽曲の拍間隔としてあり得ない拍間隔を排除し、拍構造の周期性が高い拍間隔をもつ粒子が抽出される確率を高くする。そして、粒子遷移部１３２は、入力された粒子情報から、算出した分布値ｑ（ｂ｜｛Ａｋ｝，ｂｓ’）に比例する確率で粒子ｉの拍間隔ｂｉｋを定める。
なお、粒子遷移部１３２は、自己相関Ｒ（ｂ，｛Ξｋ’｝）が相関算出部１０４から入力された場合には、Ｒ（ｂ，｛Ａｋ｝）の代わりにＲ（ｂ，｛Ξｋ’｝）を用いて分布値ｑ（ｂ｜｛Ａｋ｝，ｂｓ’）を算出してもよい。
【００８０】
粒子遷移部１３２は、生成したオンセット情報、距離値、及び抽出した拍間隔ｂｉｋに基づいて楽譜位置ｐに対する分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）を算出する。
ここで、粒子遷移部１３２は、例えば式（１２）を用いて分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）を算出する。
【００８１】
【数１２】
0005924968.tif 000013
【００８２】
式（１２）において、τ’（ｐ，ｂ）は、楽譜位置の探索範囲Ｐｉｋに含まれ、かつ時刻ｋΔＴにおける音響フレームと対応付けられた楽譜位置ｐにおけるオンセット情報ｏｐを表す。探索範囲Ｐｉｋは、前フィルタリングステップｋ−１における粒子ｉの楽譜位置ｐｉｋ−１から拍間隔ｂ＝ｂｉｋだけ進んだ楽譜位置ｐｉｋ−１’＝ｐｉｋ−１＋ΔＴ／ｂｉｋ−１を中心とし、予め定めた幅２δｐ（例えば、δｐ＝３）を有する。
即ち、式（１２）は、探索範囲Ｐｉｋにオンセットが含まれる場合、分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）は、オンセットを有する楽譜位置に対応する時刻における音響信号の振幅の累積値に比例する値であることを表す。探索範囲Ｐｉｋにオンセットが含まれない場合、１に比例する値であることを表す。楽譜位置ｐが、探索範囲Ｐｉｋの範囲外である場合には、分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）は、ゼロであることを表す。
【００８３】
ここで、楽譜位置と音響信号の振幅との関係について説明する。
図４は、楽譜情報と音響信号の関係の一例を表す概念図である。
図４の上段は、楽譜情報を表す。楽譜情報は、左から順に、１拍の音階Ａ４からなるコードｄ１、０．５拍のコードＢ４、Ｆ３からなるコードｄ２、０．５拍の音階Ｂ４、Ｇ３からなるコードｄ３を表す。この例では、音階Ｂ３を２つの０．５拍からなる音として扱う。図４の上下に延びる３本の破線は、各々オンセットを表す。図４では、オンセットを有する楽譜位置は、左端、左端から１拍、左端から１．５拍である。なお、各コードが継続する区間をセグメントといい、ｄ１、ｄ２、ｄ３は、各コードのセグメントを表す。
【００８４】
図４の下段は音響信号のスペクトログラムを表す画像である。水平方向は時刻を表し、上下方向は周波数を表す。この画像において濃淡は振幅を表す。明るい部分ほど振幅が大きく、暗い部分ほど振幅が小さいことを表す。図４の下段は、各周波数において振幅が一斉に増加し、その後、周波数が高い成分ほど早く振幅が減衰する時間変化が３回繰り返されることを表す。振幅が一斉に増加する部分が、音響信号の振幅が急激に増加する部分である。従って、図４において明るい部分を破線が通過していることは、音響信号の振幅が急激に増加する部分と楽譜情報に基づくオンセットと合致していることを表す。このとき、分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）は、大きな値をとる。言い換えれば、音響信号の振幅が急激に増加する部分が、楽譜情報に基づくオンセットと合致していない場合、分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）はより小さい値をとる。
これにより、粒子遷移部１３２は、分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）を用いることで、楽譜情報に含まれるオンセットと音響信号の振幅が増加する部分の楽譜位置ｐｉｋが推定される確率を高くすることができる。
【００８５】
図１に戻り、粒子遷移部１３２は、自己相関Ｒ（ｂ，｛Ξｋ’｝）が相関算出部１０４から入力された場合には、Ｒ（ｂ，｛Ａｋ｝）の代わりにＲ（ｂ，｛Ξｋ’｝）を用いてもよい。この場合において、粒子遷移部１３２は、例えば式（１３）を用いて分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）を算出してもよい。
【００８６】
【数１３】
0005924968.tif 000014
【００８７】
式（１３）において、ζｔは、周波数毎の距離値Ξｆ，ｔを累積した累積距離値ζｔを表す。例えば、累積距離値ζｔは、式（１４）を用いて算出される。
【００８８】
【数１４】
0005924968.tif 000015
【００８９】
式（１３）に戻り、ｐ’（τ）は、拍間隔ｂｉｋについての時刻τに対応する楽譜位置を表す。即ち、ｐ’（τ）＝ｐ−（ｔ−τ）／ｂｉｋという関係がある。この関係において、時刻ｔにおける楽譜位置がｐであると仮定されている。Ｐは、楽譜位置の探索範囲を表す。楽譜範囲Ｐは、前フィルタリングステップｋ−１における楽譜位置ｐｉｋ−１から拍間隔ｂｉｋだけ進んだ楽譜位置ｐｉｋ−１＋ΔＴ／ｂｉｋを中心とし、前後に予め定めた幅３δｐ（例えば、δｐ＝１）を有する。この楽譜位置は、粒子ｉの現フィルタリングステップｋにおける推定楽譜位置である。即ち、式（１３）は、探索範囲Ｐにオンセットが含まれる場合、分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）は、オンセットを有する楽譜位置に対応する時刻における累積距離値ζｔに比例する値であることを表す。探索範囲Ｐにオンセットが含まれない場合、１に比例する値であることを表す。楽譜位置ｐが、探索範囲Ｐの外である場合には、分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）は、ゼロであることを表す。
【００９０】
粒子遷移部１３２は、入力された粒子情報から算出した分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）に比例する確率で粒子ｉの楽譜位置ｐｉｋを定める。
なお、上述の探索範囲Ｐｉｋ−１（又はＰ）にオンセットが含まれない場合には、粒子遷移部１３２は、当該探索範囲からランダムに粒子ｉの楽譜位置ｐｉｋを抽出する。
粒子遷移部１３２は、粒子ｉ毎に算出した拍間隔の分布値（拍間隔分布値）ｑ（ｂｉｋ｜｛Ａｋ｝，ｂｓ’）に、楽譜位置の分布値（楽譜位置分布値）ｑ（ｐｉｋ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）を乗じて、粒子の分布値（粒子分布値）ｑ（ｐｉｋ，ｂｉｋ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）を算出する。
【００９１】
粒子遷移部１３２は、粒子ｉ毎にフィルタリングステップｋ−１からｋへの状態遷移確率ｑ（ｐｉｋ，ｂｉｋ｜ｐｉｋ−１，ｂｉｋ−１）を、例えば式（１５）を用いて算出する。
【００９２】
【数１５】
0005924968.tif 000016
【００９３】
式（１５）において、Ｎ（ａ｜ｂ，σ２）は、変数ａについての正規分布関数（ガウス関数、ガウシアンとも呼ばれる）を表す。ここで、原点がｂであり、分散がσ２である。σ２ｂは、拍間隔ｂの分散を表す実数値（例えば、０．２）である。ｐｉｋ−１＋ΔＴ／ｂｉｋ−１は、前フィルタリングステップｋ−１における楽譜位置ｐｉｋ−１から拍間隔ｂｉｋ−１だけ進んだ楽譜位置を表す。σ２ｐは、楽譜位置ｐの分散を表す実数値（例えば、１）である。なお、楽譜位置ｐｉｋ−１は、粒子ｉの楽譜位置情報が表す値であり、拍間隔ｂｉｋ−１は、粒子ｉの拍間隔情報が表す値である。これにより、実際の演奏によって生じる拍間隔や楽譜位置の揺らぎが考慮される。
粒子遷移部１３２は、粒子毎に粒子情報に含まれる拍間隔情報及び楽譜位置情報を、抽出した拍間隔を表す拍間隔情報及び楽譜位置を表す楽譜位置情報にそれぞれ置き換える。粒子遷移部１３２は、置き換えた粒子情報を関連度算出部１２１の周波数特性重み算出部１２２、拍間隔重み算出部１２３、粒子重み算出部１２４に出力する。
粒子遷移部１３２は、算出した粒子毎の粒子分布値ｑ（ｐｉｋ，ｂｉｋ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）と状態遷移確率ｑ（ｐｉｋ，ｂｉｋ｜ｐｉｋ−１，ｂｉｋ−１）を粒子重み算出部１２４に出力する。
【００９４】
初期（ｋ＝０）において、粒子遷移部１３２は、例えば、テンポ６０／ｂｓ’を中心とする予め定めた範囲内（ｂｓ’−６０／ｂθ＜ｂ＜ｂｓ’＋６０／ｂθ）の一様分布から粒子ｉ毎の拍間隔ｂｉｋを定める。また、粒子遷移部１３２は、粒子ｉ毎の楽譜位置ｐｉｋをゼロと定める。粒子遷移部１３２は、粒子分布値ｑ（ｐｉｋ，ｂｉｋ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）と状態遷移確率ｑ（ｐｉｋ，ｂｉｋ｜ｐｉｋ−１，ｂｉｋ−１）を１と定める。
【００９５】
次に、本実施形態に係る粒子遷移部１３２が行う粒子遷移処理について説明する。
図５は、本実施形態に係る粒子遷移処理を表すフローチャートである。
（ステップＳ２０１）粒子遷移部１３２は、相関算出部１０４から距離値及び自己相関が、楽譜情報入力部１１２から入力された楽譜情報が、再標本化部１３３から粒子毎の粒子情報が入力される。粒子遷移部１３２は、入力された楽譜情報に基づいてオンセット情報を生成する。その後、ステップＳ２０２に進む。
【００９６】
（ステップＳ２０２）粒子遷移部１３２は、入力された楽譜情報に含まれるテンポ情報に基づいて拍間隔ｂに対する窓関数ψ（ｂ｜ｂｓ’）を算出する。粒子遷移部１３２は、入力された自己相関Ｒ（ｂ，｛Ａｋ｝）に窓関数ψ（ｂ｜ｂｓ’）を乗じて拍間隔分布値ｑ（ｂ｜｛Ａｋ｝，ｂｓ’）を算出する。その後、ステップＳ２０３に進む。
（ステップＳ２０３）粒子遷移部１３２は、入力された粒子情報から、算出した拍間隔分布値ｑ（ｂ｜｛Ａｋ｝，ｂｓ’）に比例する確率で粒子ｉの拍間隔ｂｉｋを抽出する。その後、ステップＳ２０４に進む。
【００９７】
（ステップＳ２０４）粒子遷移部１３２は、生成したオンセット情報、距離値、及び抽出した拍間隔ｂｉｋに基づいて楽譜位置分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）を算出する。粒子遷移部１３２は、例えば式（１２）を用いて楽譜位置分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）を算出する。その後、ステップＳ２０５に進む。
【００９８】
（ステップＳ２０５）粒子遷移部１３２は、粒子遷移部１３２は、楽譜位置分布値ｑ（ｐ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）に比例する確率で粒子ｉの楽譜位置ｐｉｋを抽出する。但し、探索範囲Ｐｉｋ−１（又はＰ）にオンセットが含まれない場合には、粒子遷移部１３２は、探索範囲からランダムに粒子ｉの楽譜位置ｐｉｋを抽出する。その後、ステップＳ２０６に進む。
【００９９】
（ステップＳ２０６）粒子遷移部１３２は、粒子遷移部１３２は、粒子ｉ毎に拍間隔分布値ｑ（ｂｉｋ｜｛Ａｋ｝，ｂｓ’）に、楽譜位置分布値ｑ（ｐｉｋ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）を乗じて、粒子分布値ｑ（ｐｉｋ，ｂｉｋ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）を算出する。その後、ステップＳ２０７に進む。
【０１００】
（ステップＳ２０７）粒子遷移部１３２は、粒子ｉ毎に前フィルタリングステップｋ−１における楽譜位置ｐｉｋ−１、拍間隔ｂｉｋ−１、抽出された楽譜位置ｐｉｋ−１、拍間隔ｂｉｋ−１に基づいて、状態遷移確率ｑ（ｐｉｋ，ｂｉｋ｜ｐｉｋ−１，ｂｉｋ−１）を算出する。粒子遷移部１３２は、状態遷移確率ｑ（ｐｉｋ，ｂｉｋ｜ｐｉｋ−１，ｂｉｋ−１）を算出する際、例えば式（１５）を用いる。その後、ステップＳ２０８に進む。
【０１０１】
（ステップＳ２０８）粒子遷移部１３２は、粒子毎に粒子情報に含まれる拍間隔情報及び楽譜位置情報を、抽出した拍間隔を表す拍間隔情報及び楽譜位置を表す楽譜位置情報にそれぞれ置き換える。粒子遷移部１３２は、置き換えた粒子情報を関連度算出部１２１の周波数特性重み算出部１２２、拍間隔重み算出部１２３、粒子重み算出部１２４に出力する。
粒子遷移部１３２は、算出した粒子毎の粒子分布値ｑ（ｐｉｋ，ｂｉｋ｜｛Ａｋ｝，ｏｐ，ｂｉｋ）と状態遷移確率ｑ（ｐｉｋ，ｂｉｋ｜ｐｉｋ−１，ｂｉｋ−１）を粒子重み算出部１２４に出力する。その後、処理を終了する。
【０１０２】
図１に戻り、再標本化部１３３の動作について、より詳細に説明する。再標本化部１３３は、粒子重み算出部１２４から入力された粒子情報に基づいて各粒子を再標本化する。
再標本化部１３３は、入力された粒子情報に含まれる粒子重み情報が表す粒子重みに比例した確率で粒子を抽出する。再標本化部１３３は、粒子を抽出する際、例えばＳＩＲ（ｓａｍｐｌｉｎｇ  ｉｍｐｏｒｔａｎｃｅ  ｒｅｓａｍｐｌｉｎｇ）法を用いる。
再標本化部１３３は、粒子ｉ毎の粒子重み係数を、次式（１６）を用いて規格化し、規格化粒子重み係数ｗｉ，ｋ’を算出する。
【０１０３】
【数１６】
0005924968.tif 000017
【０１０４】
式（１６）において、Ｉは、予め定めた粒子数を示す整数（例えば、１５００）を表す。
再標本化部１３３は、算出した確率（規格化粒子重み係数ｗｉ，ｋ’）で次回の粒子ｉを抽出する。ここで、再標本化部１３３は、抽出された各粒子の粒子重みは、各々等しい値（例えば、１／Ｉ）とし、抽出された各粒子の楽譜位置情報及び拍間隔情報は、対応する抽出前の粒子の楽譜位置情報及び拍間隔情報と同一とする。再標本化部１３３は、抽出した粒子毎の楽譜位置情報、拍間隔情報及び粒子重みを含む粒子情報を生成する。
再標本化部１３３は、これにより、粒子数Ｉよりも多い予め定めた数の粒子を抽出する。その後、再標本化部１３３は、増加した粒子数と同一の数だけ、粒子を棄却する。ここで、再標本化部１３３は、抽出前の粒子重みが小さい粒子ほど優先して、その粒子に基づいて抽出した粒子情報を棄却する。これにより、再標本化部１３３は、粒子数（粒子情報の数）を一定数Ｉに保つ。再標本化部１３３は、棄却されずに残った粒子の粒子情報を粒子遷移部１３２に出力する。
【０１０５】
次に、拍間隔重み算出部１２３の動作について説明する。
拍間隔重み算出部１２３は、粒子遷移部１３２から入力された粒子情報に含まれる拍間隔情報を粒子ｉ毎に抽出する。拍間隔重み算出部１２３は、相関算出部１０４から入力された自己相関Ｒ（ｂ，｛Ａｋ｝）又はＲ（ｂ，｛Ξｋ’｝）から抽出した拍間隔情報が表す拍間隔ｂｉｋに対応する値Ｒ（ｂｉｋ，｛Ａｋ｝）又はＲ（ｂｉｋ，｛Ξｋ’｝）を拍間隔重み値ｗｔｍｉ，ｋと定める。拍間隔重み算出部１２３は、粒子ｉ毎の拍間隔重み値ｗｔｍｉ，ｋを表す拍間隔重み情報を粒子重み算出部１２４に出力する。
【０１０６】
なお、本実施形態では、拍間隔重み算出部１２３を省略してもよい。拍間隔重み算出部１２３を省略する場合、粒子遷移部１３２が相関算出部１０４から入力された自己相関Ｒ（ｂ，｛Ａｋ｝）又はＲ（ｂ，｛Ξｋ’｝）から拍間隔ｂｉｋに対応する値Ｒ（ｂｉｋ，｛Ａｋ｝）又はＲ（ｂｉｋ，｛Ξｋ’｝）を拍間隔重み値ｗｔｍｉ，ｋと定める。粒子遷移部１３２は、粒子ｉ毎に定めた拍間隔重み値ｗｔｍｉ，ｋを表す拍間隔重み情報を粒子重み算出部１２４に出力してもよい。
【０１０７】
次に、周波数特性重み算出部１２２について説明する。
本実施形態では、楽譜情報に基づく振幅周波数特性（例えば、楽譜スペクトログラム）が、既知の楽譜情報に含まれる楽譜位置毎の音階の基本周波数に基づく調波構造で表される成分を少なくとも１つ含むと仮定する。そこで、周波数特性重み算出部１２２は、観測値である音響信号に基づく振幅周波数特性（音響スペクトログラム）に対して尤度が最大となるように、調波構造毎の楽譜スペクトログラムへの寄与を表す変数の確率分布や、各調波構造に含まれる調波成分毎の当該調波構造への寄与を表す変数を変分パラメータとして算出する。周波数特性重み算出部１２２は、算出した変数が表す楽譜スペクトログラムが観測値として音響スペクトログラムが与えられているときの確率を算出し、算出した確率に基づいて周波数特性重み値ｗｓｐｉ，ｋを算出する。本実施形態では、上述の尤度が最大になるような周波数特性重み値ｗｓｐｉ，ｋを算出できればよいので、楽譜情報に基づく振幅周波数特性（楽譜スペクトログラム）を陽に合成しなくともよい。
【０１０８】
周波数特性重み算出部１２２は、例えば、ＬＨＡ（Ｌａｔｅｎｔ  Ｈａｒｍｏｎｉｃ  Ａｌｌｏｃａｔｉｏｎ；潜在的調波配分）法を用いて周波数特性分析部１０３から入力された音響スペクトログラムと楽譜情報入力部１１２から入力された楽譜情報に基づいて周波数特性重み情報を生成する。
ＬＨＡ法は、入力された音響信号のスペクトログラムを観測値とし、この観測値を調波構造と成分として含み、調波構造の楽譜スペクトログラムに対する寄与を表す変数の確率分布、調波成分の調波構造への寄与を表す係数の確率分布として表す音源モデルである。調波構造は、基本周波数の成分と基本周波数の整数倍の周波数の成分（高調波）を調波成分として含む周波数特性を表す。従って、基本周波数、各次数の振幅（もしくは、基準値に対する比率）が、調波構造を表す変数に含まれる。
【０１０９】
ここで、ＬＨＡ法に基づく音源モデルについて説明する。
図６は、本実施形態に係る音源モデルの一例を表す概念図である。
図６の右上に示されるｘｎは、周波数を表す。ｘｎを囲う二重円は、ｘｎが観測変数であることを表す。ｚｎからｘｎに向かう矢印、μｌからｘｎに向かう矢印及びΛｌからｘｎに向かう矢印は、周波数ｘｎの潜在変数（ｌａｔｅｎｔ  ｖａｒｉａｂｌｅ）ｚｎ、基本周波数μｌ及び精度（ｐｒｅｃｉｓｉｏｎ）Λｌに対する尤度ｑ（Ｘ｜Ｚ，｛μ｝，｛Λ｝）が仮定されていることを表す。Ｘ、Ｚ、｛μ｝、｛Λ｝は、それぞれ変数ｘｎ、ｚｎ、μｌ、Λｌの集合（セット）を表す。ｚｎ、Λｌを囲う円は、ｚｎ、Λｌがそれぞれ未知変数であることを示す。μｌが円で囲まれていないことは、基本周波数μｌが既知変数であることを表す。基本周波数μｌは、楽譜情報の基本周波数ベクトルに含まれる各音階の基本周波数であるためである。
πｄｌからｚｎに向かう矢印、θｌｍからｚｎに向かう矢印、潜在変数ｚｎの確率πｄｌ及び確率θｌｍに対する尤度ｑ（Ｚ｜｛π｝，｛θ｝）が仮定されていることを表す。ここで、｛π｝，｛θ｝は、それぞれ変数πｄｌ、θｌｍの集合（セット）を表す。πｄｌ、θｌｍを囲う円は、πｄｌ、θｌｍがそれぞれ未知変数であることを示す。末端に示されている未知変数πｄｌ、θｌｍ、Λｌに対して、それぞれ事前分布ｑ（｛π｝）、ｑ（｛θ｝）、ｑ（｛Λ｝）が仮定されていることを表す。従って、ＬＨＡ法では、確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝）は、ｑ（Ｘ｜Ｚ，｛μ｝，｛Λ｝）、ｑ（Ｘ｜Ｚ，｛μ｝，｛Λ｝）、ｑ（｛π｝）、ｑ（｛θ｝）、ｑ（｛Λ｝）の積で表される。
【０１１０】
図６において、変数を囲う四角形は、それぞれ尤度又は事前分布を与える単位、即ち、観測値を表現する確率分布を表現するための単位を表す。ｘｎとｚｎを囲う四角形は、ｘｎとｚｎが楽譜フレームｎ毎に与えられることを表す。この四角形の右下に表されているＮｄは、セグメントｄに含まれるフレーム数を表す。ｚｎ、ｘｎとともにπｄｌを囲う四角形は、πｄｌがセグメントｄ毎に与えられることを表す。この四角形の右下に表されているＤは、全セグメント数を表す。θｌｍとともにμｌ、Λｌを囲う四角形は、μｌ、Λｌが調波構造ｌ毎に与えられることを表す。この四角形の右下に表されているＬｄは、調波構造の総数を表す。θｌｍを囲う四角形は、θｌｍが調波構造ｌを構成する次数ｍの調波毎に与えられることを表す。但し、ｍ＝１は、基本周波数を表す。この四角形の右下に示されているＭは、各調波構造に含まれる調波の個数（最高次数）を表す。
【０１１１】
ＬＨＡ法では、調波構造ｌを、例えばＧＭＭを用いて表す。即ち、調波構造ｌは、次数ｍについて周波数ｍμｌを平均値とするガウス関数を確率θｌｍ（１≦ｍ≦Ｍ）で線形重み付け加算した確率分布で表される。ここで、調波構造ｌ毎の確率θｌｍの総和ΣｌＭθｌｍは１と規格化されている。即ち、確率θｌｍは、（Ｍは、各調波構造が含む調波の数を表す予め定められた整数である（例えば、Ｍ＝１０）。各次数ｍについて、ガウス関数の分散は精度Λｌである。即ち、確率θｌｍは、次数間の寄与の比率、言い換えれば次数ｍの調波成分の調波構造ｌへの寄与を表す変数である。精度Λｌは、基本周波数μｌの誤差の大きさを表す。
【０１１２】
図７は、本実施形態に係る調波構造の一例を表す図である。
図７において、縦軸は確率を表し、横軸は周波数を表す。図７において、細い破線は調波構造ｌが、周波数μｌ，２μｌ，…，Ｍμｌにおいて確率のピーク値θｌ１，θｌ２，…，θｌＭを有することを表す。一点破線は、ピーク値θｌ１，θｌ２，…，θｌＭを表す。周波数２μｌを中心とする確率値の曲線（ガウス関数）を挟む２つの矢印の間隔は、精度Λｌを表す。
【０１１３】
ＬＨＡ法では、観測値である音響スペクトログラムを周波数ｘｎの確率分布とみなす。本実施形態では、周波数ｘｎの確率分布を、調波構造ｌを確率πｄｌ（１≦ｌ≦ｌＬｄ）で線形重み付け加算した確率分布で近似する。即ち、確率πｄｌは、調波構造ｌの楽譜スペクトラムへの寄与を表す変数であり、この楽譜スペクトログラムは陽に算出されないがをも音響スペクトログラムを近似する。本実施形態では、音響スペクトログラム、調波構造ｌに含まれる確率値（即ち、振幅）を予め定めた量子化幅で量子化してもよい。これにより、確率分布は、振幅を度数とするヒストグラムとして表されるため演算量が低減する。
【０１１４】
図８は、本実施形態に係るヒストグラムの一例を表す図である。
図８の左上は、調波構造１を表すヒストグラムである。図８の左下は、調波構造ｌを現すヒストグラムである。図８の右側は、周波数ｘｎの確率分布を表すヒストグラムである。これらのヒストグラムにおいて、縦軸は振幅を表し、横軸は周波数を表す。調波構造１の振幅は、周波数μ１、２μ１、３μ１、４μ１、５μ１、６μ１においてピークを有する。各ピーク値の大きさは、楽器の音色などに応じて異なる。調波構造ｌの振幅は、周波数μｌ、２μｌにおいてピークを有する。図８に示す例では、周波数μｌは周波数μ１よりも高く、周波数μｌにおけるピーク値は、周波数μ１におけるピーク値よりも小さい。
調波構造１を表すヒストグラムの右側に表されているπｄ１は、調波構造１に対する確率を表す。調波構造ｌを表すヒストグラムの右側に表されているπｄｌは、調波構造ｌに対する確率を表す。図８の最左辺の上下に延びる括弧と「Ｌｄ個」の文字は、Ｌｄ個の調波構造ｌを確率πｄｌで重み付け加算して、図８の右側における周波数ｘｎの確率分布を近似することを表す。
【０１１５】
従って、本実施形態では観測値Ｘが与えられているときに、全ての未知変数の事後分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐｉｋ，｛μ｝）を粒子ｉ毎に算出する。全ての未知変数に対する事後分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐｉｋ，｛μ｝）を累積した累積確率値は、本音源モデルの適合性を表す。周波数特性重み算出部１２２は、例えば、式（１７）を用いて算出した累積確率値を周波数特性重み値の対数値ｌｏｇｗｓｐｉ，ｋとして算出する。
【０１１６】
【数１７】
0005924968.tif 000018
【０１１７】
しかし、式（１７）を用いた演算、とりわけＺに対する累算を解析的に行うことは困難である。そこで、本実施形態では、変分ベイズ法（ｖａｒｉａｔｉｏｎａｌ  Ｂａｙｅｓ［ＶＢ］  ｍｅｔｈｏｄ）を用いて、事後分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝｜Ｘ，ｐｉｋ，｛μ｝）を近似する。本実施形態では、解析的に求めることができる変分事後分布ｑ＊（Ｚ，｛π｝，｛θ｝，｛Λ｝）を仮定し、真の事後分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝｜Ｘ，ｐｉｋ，｛μ｝）との指標値（例えば、カルバック・ライブラー情報量［Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ  Ｄｉｖｅｒｇｅｎｃｅ；ＫＬ−ｄｉｖ］）を最小化する変分パラメータを算出する。ＫＬ−ｄｉｖ  ＤＫＬ（ｑ｜ｑ’）は、式（１８）に示すように、ある確率分布ｑ（ｊ）と他の確率分布ｑ’（ｊ）の差異を表す指標値である。
【０１１８】
【数１８】
0005924968.tif 000019
【０１１９】
変分ベイズ法では、全ての未知変数についての変分事後分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝）が、未知変数Ｚの変分事後分布ｑ（Ｚ）と未知変数｛π｝、｛θ｝、｛Λ｝の変分事後分布ｑ（｛π｝，｛θ｝，｛Λ｝）の積に因数分解できると仮定する。
後述するように、例えば、式（１９）で表される変分下限値（ｖａｒｉａｔｉｏｎａｌ  ｌｏｗｅｒ  ｂｏｕｎｄ）Ｌ（ｑ）が収束する変分パラメータを算出する。
【０１２０】
【数１９】
0005924968.tif 000020
【０１２１】
式（１９）において、ＥＺ，｛π｝，｛θ｝，｛Λ｝［…］は、変分事後分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝）に対する…の期待値を表す。
そして、変分下限値Ｌ（ｑ）を周波数特性重み値の対数値ｌｏｇｗｓｐｉ，ｋとして算出する。変分下限値Ｌ（ｑ）と変分事後分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝）と確率分布ｑ（Ｚ，｛π｝，｛θ｝，｛Λ｝｜Ｘ，ｐｉｋ，｛μ｝）の間のＫＬ−ｄｉｖの和は、常に式（１７）の右辺と一致するため、変分ベイズ法によりＫＬ−ｄｉｖを最小化することで、式（１７）の近似値として変分下限値Ｌ（ｑ）を代用する。
変分ベイズ法では、変分パラメータの算出において式（２０）、（２１）に示す関係が仮定されている。これらの関係は、ともに変分事後分布ｑ＊（｛Ｚ｝）、ｑ＊（｛π｝、｛θ｝、｛Λ｝）が最適であることを表す。
【０１２２】
【数２０】
0005924968.tif 000021
【０１２３】
式（２０）は、潜在変数｛Ｚ｝の変分事後分布ｑ＊（｛Ｚ｝）が、粒子ｉ毎の確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐｉｋ，｛μ｝）の、その他の未知変数｛π｝、｛θ｝、｛Λ｝の変分事後分布ｑ＊（｛π｝、｛θ｝、｛Λ｝）に対する期待値に比例することを表す。この関係を満足するように変分パラメータを算出する処理は、ＶＢ−Ｅ（Ｅｘｐｅｃｔａｔｉｏｎ）ステップと呼ばれる。
【０１２４】
【数２１】
0005924968.tif 000022
【０１２５】
式（２１）は、未知変数｛π｝、｛θ｝、｛Λ｝の変分事後分布ｑ＊（｛π｝、｛θ｝、｛Λ｝）が、粒子ｉ毎の確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐｉｋ，｛μ｝）の、潜在変数Ｚの変分事後分布ｑ＊（｛Ｚ｝）に対する期待値に比例することを表す。この関係を満足するように変分パラメータを算出する処理は、ＶＢ−Ｍ（Ｍａｘｉｍｉｚａｔｉｏｎ）ステップと呼ばれる。
【０１２６】
次に、ＬＨＡ法に変分ベイズ法を適用したときの周波数特性重み算出部１２２の構成について説明する。
図９は、本実施形態に係る周波数特性重み算出部１２２の構成を表す概略図である。
周波数特性重み算出部１２２は、データ入力部１２２１、第１変分パラメータ算出部１２２２、第２変分パラメータ算出部１２２３、変分下限値算出部１２２４、収束判定部１２２５、及び重み算出部１２２６を含んで構成される。
【０１２７】
データ入力部１２２１は、周波数特性分析部１０３から音響スペクトログラムＡｆ，ｔが、楽譜情報入力部１１２から楽譜情報が、粒子遷移部１３２から粒子毎の粒子情報がそれぞれ入力される。データ入力部１２２１は、音響スペクトログラムＡｆ，ｔを予め定めた量子化幅ΔＡで量子化して量子化音響スペクトログラム［Ａｆ，ｔ］を算出する。ここで、算出した量子化音響スペクトログラム［Ａｆ，ｔ］は、音響フレームｎにおける周波数ｘｎ毎の度数（確率）としたヒストグラムを表す。つまり、このヒストグラムは周波数ｘｎが発生する回数が量子化音響スペクトログラム［Ａｆ，ｔ］が表す振幅に比例することを表す。
【０１２８】
データ入力部１２２１は、粒子情報から粒子ｉ毎に楽譜位置情報ｐｉｋ、拍間隔情報ｂｉｋを抽出し、入力された楽譜情報から楽譜位置ｐｉｋから拍間隔ｂｉｋ前の楽譜位置ｐｉｋ−ｂｉｋまでの楽譜情報を抽出する。データ入力部１２２１は、抽出した楽譜情報からセグメントｄ毎の基本周波数ベクトルμｌを抽出する。
データ入力部１２２１は、算出した量子化音響スペクトログラム［Ａｆ，ｔ］から、粒子ｉ毎の楽譜位置ｐｉｋに対応する時刻ｋΔＴから楽譜位置ｐｉｋ−ｂｉｋに対応する時刻までの区間の量子化音響スペクトログラム［Ａｆ，ｔ］を抽出する。
データ入力部１２２１は、粒子ｉ毎に抽出した区間の量子化音響スペクトログラム［Ａｆ，ｔ］とセグメントｄ毎の基本周波数ベクトルμｐを、第１変分パラメータ算出部１２２２、第２変分パラメータ算出部１２２３、及び変分下限値算出部１２２４に出力する。
【０１２９】
第１変分パラメータ算出部１２２２は、データ入力部１２２１から、量子化音響スペクトログラム［Ａｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμｐが入力される。また、第２変分パラメータ算出部１２２３から変分パラメータαｄｌ、βｌｍ、ａｌ、ｂｌが入力される。
第１変分パラメータ算出部１２２２は、入力されたパラメータを用いて、上述のＶＢ−Ｅステップを実行する。即ち、第１変分パラメータ算出部１２２２は、例えば式（２２）を用いて、ｄ、ｎ、ｌ、ｍ毎にパラメータρｄｎｌｍを算出する。
【０１３０】
【数２２】
0005924968.tif 000023
【０１３１】
式（２２）において、ψ（…）は、…のディガンマ（ｄｉｇａｍｍａ）関数を表す。
αｄ・は、αｄｌのｌに対する総和Σｌαｄｌである。β・ｌは、βｄｌのｄに対する総和Σｄβｄｌである。
次に、第１変分パラメータ算出部１２２２は、式（２３）に基づいてｄ、ｎ、ｌ、ｍ毎にパラメータρｄｎｌｍをパラメータρｄｎ・・で除算して変分パラメータγｄｎｌｍを算出する。
【０１３２】
【数２３】
0005924968.tif 000024
【０１３３】
式（２３）において、パラメータρｄｎ・・は、パラメータρｄｎｌｍのｌ、ｍに対する総和ΣｌΣｍρｄｎｌｍである。
変分パラメータγｄｎｌｍは、負担率とも呼ばれ、観測した周波数ｘｎが調波構造ｌの第ｍ高調波に対応するガウス分布の事後確率を表す。
なお、式（２２）、（２３）は、潜在変数｛Ｚ｝の変分事後分布ｑ＊（｛Ｚ｝）が、変分パラメータγｄｎｌｍのｚｄｎｌｍ乗の積からなる多項分布（式（２４））に従うことを仮定して導出される。変分パラメータγｄｎｌｍのｚｄｎｌｍ乗の積は、変分パラメータγｄｎｌｍがｚｄｎｌｍ回発生する確率である。
【０１３４】
【数２４】
0005924968.tif 000025
【０１３５】
第１変分パラメータ算出部１２２２は算出した変分パラメータγｄｎｌｍを第２変分パラメータ算出部１２２３に出力する。
【０１３６】
第２変分パラメータ算出部１２２３は、データ入力部１２２１から、量子化音響スペクトログラム［Ａｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμｐが入力される。また、第１変分パラメータ算出部１２２２から変分パラメータγｄｎｌｍが入力される。
第２変分パラメータ算出部１２２３は、入力されたパラメータを用いて、上述のＶＢ−Ｍステップを実行する。ここで、第２変分パラメータ算出部１２２３は、例えば式（２５）を用いて、ｄ、ｌ毎にパラメータα０ｌにパラメータγｄ・ｌ・を加算して変分パラメータαｄｌを算出する。
【０１３７】
【数２５】
0005924968.tif 000026
【０１３８】
式（２５）において、パラメータγｄ・ｌ・は、変分パラメータγｄｎｌｍのｎ、ｍに対する総和ΣｌΣｍγｄｎｌｍである。パラメータα０ｌは、予め定められた実数を表す。
第２変分パラメータ算出部１２２３は、例えば式（２６）を用いて、ｌ、ｍ毎にパラメータβ０ｍにパラメータγ・・ｌｍを加算して変分パラメータβｌｍを算出する。
【０１３９】
【数２６】
0005924968.tif 000027
【０１４０】
式（２６）において、パラメータγ・・ｌｍは、変分パラメータγｄｎｌｍのｄ、ｎに対する総和ΣｄΣｎγｄｎｌｍである。パラメータβ０ｍは、予め定められた実数を表す。
第２変分パラメータ算出部１２２３は、例えば式（２７）を用いて、ｌ毎にパラメータａ０にパラメータγ・・ｌ・を加算して変分パラメータａｌを算出する。
【０１４１】
【数２７】
0005924968.tif 000028
【０１４２】
式（２７）において、パラメータγ・・ｌ・は、変分パラメータγｄｎｌｍのｄ、ｎ、ｍに対する総和ΣｄΣｎΣｍγｄｎｌｍである。ａ０は、予め定められた実数を表す。
第２変分パラメータ算出部１２２３は、例えば式（２８）を用いて、変分パラメータｂｌを算出する。
【０１４３】
【数２８】
0005924968.tif 000029
【０１４４】
式（２８）において、ｂ０は、予め定められた実数を表す。
なお、式（２５）は、未知変数｛πｄ｝の変分事後分布ｑ＊（｛π｝）が変分パラメータ｛αｄ｝に対するディリクレ分布の積（式（２９））に従うことを仮定して導出される。
【０１４５】
【数２９】
0005924968.tif 000030
【０１４６】
式（２９）において、Ｄｉｒ（｛πｄ｝｜｛αｄ｝）は、｛πｄ｝の｛αｄ｝に対するディリクレ分布を表す。｛πｄ｝は、変数πｄｌのｄに対する集合（セット）を表す。｛αｄ｝は、変数αｄｌのｄに対するセットを表す。式（２９）のディリクレ分布は、複数の事象が各々αｄｌ回発生する確率がπｄｌとなる確率を与える。
式（２６）は、未知変数｛θｌ｝の変分事後分布ｑ＊（｛θ｝）が変分パラメータ｛θｌ｝に対するディリクレ分布の積（式（３０））に従うことを仮定して導出される。
【０１４７】
【数３０】
0005924968.tif 000031
【０１４８】
式（３０）において、Ｄｉｒ（｛θｌ｝｜｛βｌ｝）は、｛θｌ｝の｛βｌ｝に対するディリクレ分布を表す。｛θｌ｝は、変数θｌｍのｌに対する集合（セット）を表す。｛βｌ｝は、変数βｌｍのｍに対する集合（セット）を表す。式（３０）のディリクレ分布は、複数の事象が各々βｄｌ回発生する確率がθｌｍとなる確率を与える。
式（２７）、（２８）は、未知変数｛Λ｝の変分事後分布ｑ＊（｛Λ｝）が変分パラメータａｌ、ｂｌに対するガンマ分布の積（式（３１））に従うことを仮定して導出される。
【０１４９】
【数３１】
0005924968.tif 000032
【０１５０】
式（３１）において、Ｇａｍ（Λｌ｜ａｌ，ｂｌ）は、変数Λｌのガンマ関数を表す。ａｌ、ｂｌは、それぞれガンマ関数の形状パラメータ、比率パラメータを表す。なお、変分事後分布ｑ＊（｛π｝、｛θ｝、｛Λ｝）は、式（２９）、（３０）、（３１）で各々表される変分事後分布の積ｑ＊（｛π｝）、ｑ＊（｛θ｝）ｑ＊（｛Λ｝）で因数分解できることが仮定されている。
第２変分パラメータ算出部１２２３は、算出した変分パラメータαｄｌ、βｌｍ、ａｌ、ｂｌと、第１変分パラメータ算出部１２２２から入力された変分パラメータγｄｎｌｍを変分下限値算出部１２２４に出力する。収束判定部１２２５から変分下限値Ｌ（ｑ）が収束したことを表す変分下限値信号が入力されたとき、第２変分パラメータ算出部１２２３は、第１変分パラメータ算出部１２２２への変分パラメータαｄｌ、βｌｍ、ａｌ、ｂｌの出力を停止する。
【０１５１】
変分下限値算出部１２２４は、第２変分パラメータ算出部１２２３から変分パラメータαｄｌ、βｌｍ、ａｌ、ｂｌ、γｄｎｌｍが入力される。変分下限値算出部１２２４は、データ入力部１２２１から入力された量子化音響スペクトログラム［Ａｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμｐに基づいて、例えば式（３２）に基づいて確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐｉｋ，｛μ｝）を算出する。
【０１５２】
【数３２】
0005924968.tif 000033
【０１５３】
式（３２）において、Ｘは入力された区間の量子化音響スペクトログラム［Ａｆ，ｔ］に基づく周波数ｘｎの集合（即ち、上述のヒストグラム）を表す。｛μ｝は入力されたセグメントｄ毎の基本周波数ベクトル｛μｐ｝（但し、楽譜フレームｐは、セグメントｄに属す）に含まれる基本周波数μｌの集合（セット）を表す。ｑ（Ｘ｜Ｚ，｛μ｝,Λ）は、例えば、式（３３）に示されるように、平均がｍ次の基本周波数μｌ、分散Λｌである、周波数ｘｎのガウス関数の積である。
【０１５４】
【数３３】
0005924968.tif 000034
【０１５５】
式（３２）において、ｑ（Ｚ｜｛π｝,｛θ｝）は、例えば、式（３４）に示されるようにπｄｌθｌｍのｚｄｎｌｍ乗の積からなる多項分布である。
【０１５６】
【数３４】
0005924968.tif 000035
【０１５７】
式（３２）において、ｑ（｛π｝）は、例えば、式（３５）に示されるように変数｛πｄ｝のパラメータセット｛α０｝に対するディリクレ分布の積である。
【０１５８】
【数３５】
0005924968.tif 000036
【０１５９】
パラメータセット｛α０｝は、上述のパラメータα０ｌの集合（セット）である。
式（３２）において、ｑ（｛θ｝）は、例えば、式（３６）に示されるように変数｛θｌ｝のパラメータセット｛β０｝に対するディリクレ分布の積である。
【０１６０】
【数３６】
0005924968.tif 000037
【０１６１】
パラメータセット｛β０｝は、上述のパラメータβ０ｍの集合（セット）である。
式（３２）において、ｑ（｛Λ｝）は、例えば、式（３７）に示されるように分散Λｌのパラメータａ０、ｂ０に対するガンマ分布の積である。
【０１６２】
【数３７】
0005924968.tif 000038
【０１６３】
パラメータａ０、ｂ０の値は、上述のパラメータａ０、ｂ０と同一の値である。
変分下限値算出部１２２４は、入力された変分パラメータαｄｌ、βｌｍ、ａｌ、ｂｌ、γｄｎｌｍに基づいて、式（２９）、（３０）、（３１）、（２４）を用いて、変分事後分布ｑ＊（｛π｝）、ｑ＊（｛θ｝）、ｑ＊（｛Λ｝）、ｑ＊（Ｚ）を算出する。変分下限値算出部１２２４は、算出したｑ＊（｛π｝）、ｑ＊（｛θ｝）、ｑ＊（｛Λ｝）、ｑ＊（Ｚ）の積をとり全未知変数についての変分事後分布ｑ＊（Ｚ，π，θ，Λ）を算出する。
【０１６４】
変分下限値算出部１２２４は、算出した確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐｉｋ，｛μ｝）と変分事後分布ｑ＊（Ｚ，π，θ，Λ）を、例えば式（１９）を用いて粒子ｉ毎の変分下限値Ｌ（ｑ）を算出する。変分下限値算出部１２２４は、式（１９）においてｑ（Ｚ，π，θ，Λ）の代わりに、算出したｑ＊（Ｚ，π，θ，Λ）を代入する。
なお、変分事後分布ｑ＊（｛π｝）、ｑ＊（｛θ｝）、ｑ＊（｛Λ｝）、ｑ＊（Ｚ）は、それぞれ式（２９）、（３０）、（３１）、（２４）で表され、確率分布ｑ（Ｘ｜Ｚ，｛μ｝,Λ）ｑ（Ｚ｜｛π｝,｛θ｝）、ｑ（｛π｝）、ｑ（｛θ｝）、ｑ（｛Λ｝）は、それぞれ式（３３）、（３４）、（３５）、（３６）、（３７）で表される関数を仮定しているため、変分下限値算出部１２２４は、式（１９）を用いて変分下限値Ｌ（ｑ）を解析的に算出することができる。かかる仮定により、事前分布と変分分布の関数形は相似、即ち共役になる。
変分下限値算出部１２２４は、算出した粒子ｉの変分下限値Ｌ（ｑ）を収束判定部１２２５に出力する。
【０１６５】
収束判定部１２２５は、変分下限値算出部１２２４から入力された粒子ｉ毎の変分下限値Ｌ（ｑ）が、収束したか否か判断する。収束判定部１２２５は、例えば、前回入力された変分下限値Ｌ（ｑ）との差分の絶対値が、予め定めた値よりも小さくなったとき、変分下限値Ｌ（ｑ）が収束したと判断する。
収束判定部１２２５は、変分下限値Ｌ（ｑ）が収束したと判断したとき、変分下限値Ｌ（ｑ）が収束したことを表す変分下限値信号を第２変分パラメータ算出部１２２３に出力する。また、収束判定部１２２５は、粒子ｉ毎の変分下限値Ｌ（ｑ）を重み算出部１２２６に出力する。
【０１６６】
重み算出部１２２６は、収束判定部１２２５から入力された粒子ｉの変分下限値Ｌ（ｑ）に基づいて粒子ｉの周波数特性重み値ｗｓｐｉ，ｋを算出する。重み算出部１２２６は、例えば、ｗｓｐｉ，ｋ＝ｅｘｐ（Ｌ（ｑ））と算出する。重み算出部１２２６は、算出した周波数特性重み値ｗｓｐｉ，ｋを表す周波数特性重み情報を粒子重み算出部１２４に出力する。
【０１６７】
次に、周波数特性重み算出部１２２が行う周波数特性重み算出処理について説明する。
図１０は、本実施形態における周波数特性重み算出処理の一例を表すフローチャートである。
（ステップＳ３０１）データ入力部１２２１は、周波数特性分析部１０３から入力された音響スペクトログラムＡｆ，ｔを予め定めた量子化幅ΔＡで量子化して量子化音響スペクトログラム［Ａｆ，ｔ］を算出する。
データ入力部１２２１は、粒子遷移部１３２から入力された粒子情報から粒子ｉ毎に楽譜位置情報、拍間隔情報を抽出する。データ入力部１２２１は、楽譜情報入力部１１２から入力された楽譜情報から楽譜位置ｐｉｋから拍間隔ｂｉｋ前の楽譜位置ｐｉｋ−ｂｉｋまでの楽譜情報を抽出し、抽出した楽譜情報からセグメントｄ毎の基本周波数ベクトルμｌを抽出する。
【０１６８】
データ入力部１２２１は、算出した量子化音響スペクトログラム［Ａｆ，ｔ］から、粒子ｉ毎の楽譜位置ｐｉｋに対応する時刻ｋΔＴから楽譜位置ｐｉｋ−ｂｉｋに対応する時刻までの区間の量子化音響スペクトログラム［Ａｆ，ｔ］を抽出する。
データ入力部１２２１は、粒子ｉ毎に抽出した区間の量子化音響スペクトログラム［Ａｆ，ｔ］とセグメントｄ毎の基本周波数ベクトルμｐを、第１変分パラメータ算出部１２２２、第２変分パラメータ算出部１２２３、及び変分下限値算出部１２２４に出力する。その後、ステップＳ３０２に進む。
【０１６９】
（ステップＳ３０２）第１変分パラメータ算出部１２２２は、データ入力部１２２１から、量子化音響スペクトログラム［Ａｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμｐが入力される。また、第１変分パラメータ算出部１２２２は、第２変分パラメータ算出部１２２３から変分パラメータαｄｌ、βｌｍ、ａｌ、ｂｌが入力される。
第１変分パラメータ算出部１２２２は、入力されたパラメータに基づいて、例えば式（２２）を用いて、パラメータρｄｎｌｍを算出し、式（２３）を用いて変分パラメータγｄｎｌｍを算出する。
第１変分パラメータ算出部１２２２は算出した変分パラメータγｄｎｌｍを第２変分パラメータ算出部１２２３に出力する。その後、ステップＳ３０３に進む。
【０１７０】
（ステップＳ３０３）第２変分パラメータ算出部１２２３は、データ入力部１２２１から、量子化音響スペクトログラム［Ａｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμｐが入力される。また、第１変分パラメータ算出部１２２２から変分パラメータγｄｎｌｍが入力される。第２変分パラメータ算出部１２２３は、入力された変分パラメータγｄｎｌｍに基づいて、例えば式（２５）、（２６）、（２７）、（２８）を用いて、変分パラメータαｄｌ、βｌｍ、ａｌ、ｂｌを算出する。第２変分パラメータ算出部１２２３は、変分パラメータαｄｌ、βｌｍ、ａｌ、ｂｌを第１変分パラメータ算出部１２２２と変分下限値算出部１２２４に出力する。第２変分パラメータ算出部１２２３は、第１変分パラメータ算出部１２２２から入力された変分パラメータγｄｎｌｍを変分下限値算出部１２２４に出力する。その後、ステップＳ３０４に進む。
【０１７１】
（ステップＳ３０４）変分下限値算出部１２２４は、データ入力部１２２１から入力された量子化音響スペクトログラム［Ａｆ，ｔ］、セグメントｄ毎の基本周波数ベクトルμｐに基づいて、例えば式（３２）に基づいて確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐｉｋ，｛μ｝）を算出する。変分下限値算出部１２２４は、第２変分パラメータ算出部１２２３から入力された変分パラメータαｄｌ、βｌｍ、ａｌ、ｂｌ、γｄｎｌｍに基づいて変分事後分布ｑ＊（Ｚ，π，θ，Λ）を算出する。
変分下限値算出部１２２４は、算出した確率分布ｑ（Ｘ，Ｚ，｛π｝，｛θ｝，｛Λ｝｜ｐｉｋ，｛μ｝）と変分事後分布ｑ＊（Ｚ，π，θ，Λ）を、例えば式（１９）を用いて粒子ｉ毎の変分下限値Ｌ（ｑ）を算出する
変分下限値算出部１２２４は、算出した粒子ｉの変分下限値Ｌ（ｑ）を収束判定部１２２５に出力する。その後、ステップＳ３０５に進む。
【０１７２】
（ステップＳ３０５）収束判定部１２２５は、変分下限値算出部１２２４から入力された粒子ｉ毎の変分下限値Ｌ（ｑ）が、収束したか否か判断する。ここで、収束判定部１２２５は、変分下限値算出部１２２４から前回入力された変分下限値Ｌ（ｑ）との差分の絶対値が、予め定めた値よりも小さくなったとき、変分下限値Ｌ（ｑ）が収束したと判断する。変分下限値Ｌ（ｑ）が収束したと判断したとき（ステップＳ３０５  Ｙ）、第２変分パラメータ算出部１２２３は、変分パラメータαｄｌ、βｌｍ、ａｌ、ｂｌの第１変分パラメータ算出部１２２２への出力を停止する。収束判定部１２２５は、粒子ｉ毎の変分下限値Ｌ（ｑ）を重み算出部１２２６に出力する。その後、ステップＳ３０６に進む。
変分下限値Ｌ（ｑ）が収束していないと判断したとき（ステップＳ３０５  Ｎ）、ステップＳ３０２に進む。
【０１７３】
（ステップＳ３０６）重み算出部１２２６は、算出した粒子ｉの変分下限値Ｌ（ｑ）に基づいて粒子ｉの周波数特性重み値ｗｓｐｉ，ｋを算出する。重み算出部１２２６は、算出した周波数特性重み値ｗｓｐｉ，ｋを表す周波数特性重み情報を粒子重み算出部１２４に出力する。その後、処理を終了する。
【０１７４】
以上、説明したように、本実施形態では、入力された音響信号の第１の周波数特性と楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性との関連度を表す重み値を、第２の周波数特性に含まれる成分毎の第１の周波数特性の確率分布に基づいて、楽譜の位置毎に算出する。また、本実施形態では、算出した重み値に基づいて入力された音響信号に対応する楽譜の位置を探索する。これにより、音響信号に基づく第１の周波数特性と、楽譜情報における音階に基づく第２の周波数特性が完全に合致していなくとも、両者の関連性を検出することができる。そのため、演奏された音楽に対する楽譜の位置をより頑健に推定することが可能になる。
また、本実施形態では、第２の周波数特性は少なくとも１つの調波構造を含み、調波構造に含まれる調波成分毎の調波構造への寄与を表す変数の確率分布と、調波構造毎の第２の周波数特性への寄与を表す変数の確率分布に基づいて第１の周波数特性についての尤度を最大にするように、重み値を算出する。そのため、楽譜情報に基づく第２の周波数成分に含まれる調波構造の周波数特性の変化、調波構造間の結合特性の変化に柔軟に対応できるので、音響信号に基づく第１の周波数特性との関連性をより的確に推定することができる。そのため、演奏された音楽に対する楽譜の位置をより正確に推定することが可能になる。
【０１７５】
（第２の実施形態）
次に、本発明の第２の実施形態について説明する。第１の実施形態と同一の構成又は処理については、同一の番号を付す。以下、第１の実施形態との差異点を主として説明する。
図１１は、本実施形態に係る楽譜位置推定装置２の構成を示す概略図である。
楽譜位置推定装置２は、楽譜位置推定装置１（図１）と同様に、音響信号入力部１０１、音響特徴量生成部１０２、楽譜情報記憶部１１１、楽譜情報入力部１１２、関連度算出部１２１、楽譜位置探索部１３１、及び演奏情報出力部１４１を含んで構成される。但し、関連度算出部１２１は、周波数特性重み算出部２２２、拍間隔重み算出部１２３、及び粒子重み算出部１２４を含んで構成される。即ち、楽譜位置推定装置２は、周波数特性重み算出部１２２の代わりに周波数特性重み算出部２２２を備える点で、楽譜位置推定装置１と異なる。
【０１７６】
本実施形態でも、楽譜情報に基づく振幅周波数特性（例えば、楽譜スペクトログラム）が、既知の楽譜情報に含まれる楽譜位置毎の音階の基本周波数に基づく調波構造で表される成分を少なくとも１つ含むと仮定する。そこで、周波数特性重み算出部２２２は、楽譜スペクトログラムに含まれる周波数毎の振幅が表す度数（振幅期待値）が、対応する周波数における音響スペクトログラムの振幅が表す度数だけ発生する確率分布に基づいて周波数特性重み値ｗｓｐｉ，ｋを算出する。周波数特性重み算出部２２２は、例えば、ＰＤＡ（Ｐｏｉｓｓｏｎ  Ｄｉｓｔｒｉｂｕｔｉｏｎ  Ａｍｐｌｉｔｕｄｅ；ポアソン分布振幅）を用いて周波数特性重み値ｗｉ，ｋｓｐを算出する。但し、周波数特性重み算出部２２２は、楽譜スペクトログラムの振幅の周波数方向に累積した総和が、予め定めた値（例えば、１）となるように各周波数の振幅を規格化する。
【０１７７】
ＰＤＡは、量子化音響スペクトログラムの振幅期待値（ｅｘｐｅｃｔｅｄ  ａｍｐｌｉｔｕｄｅ）に対するポアソン分布の関数値を周波数にわたり累積した累積値である。振幅期待値は、楽譜スペクトログラムに、量子化音響スペクトログラムを周波数にわたり累積した総和を乗じた値である。また、上述のように、量子化音響スペクトログラムは、周波数毎の度数を表すヒストグラムとみることもできる。上述のポアソン分布は、楽譜情報による周波数毎の平均度数だけ発生する事象が、入力音響信号による周波数毎の度数だけ発生する確率を表す。従って、ＰＤＡは、音響スペクトログラムと楽譜スペクトログラムとの差異の度合いを表す指標値を表す。例えば、ＰＤＡは、音響スペクトログラムのピーク周波数と合致する楽譜スペクトログラムのピーク周波数があれば大きい値をとる。ＰＤＡは、音響スペクトログラムのピーク周波数と合致する楽譜スペクトログラムがなければ小さい値をとる。そのため、ＰＤＡは、音響スペクトログラムと楽譜スペクトログラムの調波構造の差異に対して鋭敏に変化する。
【０１７８】
次に、周波数特性重み算出部２２２の構成について説明する。
図１２は、本実施形態に係る周波数特性重み算出部２２２の構成を表す概略図である。
周波数特性重み算出部１２２は、周波数特性量子化部２２２１、周波数特性合成部２２２２、振幅期待値算出部２２２３、振幅分布算出部２２２４及び重み算出部２２２５を含んで構成される。
【０１７９】
周波数特性量子化部２２２１は、周波数特性分析部１０３から入力された音響スペクトログラムＡｆ，ｔを予め定めた量子化幅ΔＡで量子化して量子化音響スペクトログラム［Ａｆ，ｔ］を算出する。周波数特性量子化部２２２１は、算出した量子化音響スペクトログラム［Ａｆ，ｔ］を振幅期待値算出部２２２３及び振幅分布算出部２２２４に出力する。
【０１８０】
周波数特性合成部２２２２は、楽譜情報入力部１１２から入力された楽譜情報から楽譜フレームｐ毎の基本周波数ベクトルμｐを抽出し、基本周波数ベクトルμｐに含まれる基本周波数μｐｌ毎に予め定めた音源モデル（例えば、ＧＭＭ）を用いて調波構造を生成する。周波数特性合成部２２２２は、生成した調波構造を基本周波数μｐｌにわたって累積して楽譜スペクトログラムＡ＾ｆ，ｐを生成する。生成した楽譜スペクトログラムＡ＾ｆ，ｐは、周波数にわたり累積した総和ΣｆＡ＾ｆ，ｐが１となるように規格化されている。
周波数特性合成部２２２２は、例えば式（３８）を用いて楽譜スペクトログラムＡ＾ｆ，ｐを生成する。
【０１８１】
【数３８】
0005924968.tif 000039
【０１８２】
式（３８）において、Ｃｈａｒｍは、例えば、式（３９）に示すように、楽譜スペクトログラムＡ＾ｆ，ｐの周波数毎の振幅値を周波数にわたり累積した総和が、０より大きく１より小さい実数（式（３９）では、例えば０．９）を与える定数である。
【０１８３】
【数３９】
0005924968.tif 000040
【０１８４】
式（３８）において、Ｍ、θｌｍ、Λｌは、例えば、それぞれ１０、０．２（ｍ−１）、０．８である。Ｃｆｌｏｏｒは、Ｃｈａｒｍとの和が１となる実数（例えば、０．１）である。Ｃｆｌｏｏｒを定めることにより、楽譜スペクトログラムＡ＾ｆ，ｐは、０より大きい実数値をとるため、後述する処理において零除算が回避される。
周波数特性合成部２２２２は、生成した楽譜スペクトログラムＡ＾ｆ，ｐを振幅期待値算出部２２２３に出力する。
【０１８５】
振幅期待値算出部２２２３は、周波数特性量子化部２２２１から量子化音響スペクトログラム［Ａｆ，ｔ］が、周波数特性合成部２２２２から楽譜スペクトログラムＡ＾ｆ，ｐが、粒子遷移部１３２から粒子毎の粒子情報がそれぞれ入力される。
振幅期待値算出部２２２３は、入力された粒子情報から粒子ｉ毎に楽譜位置情報ｐｉｋを抽出する。振幅期待値算出部２２２３は、抽出した楽譜位置ｐｉｋを基準（例えば、楽譜位置ｐｉｋに対応する時刻τがｋΔＴ）とした対応する時刻τ毎に、量子化音響スペクトログラム［Ａｆ，τ］の周波数ｆにわたった総和Ａ・τを算出する。
振幅期待値算出部２２２３は、入力された楽譜スペクトログラムＡ＾ｆ，ｐを、対応する楽譜位置ｐに対応する時刻τ毎に算出した総和Ａ・τを乗算して振幅期待値λｆ，τを算出する。振幅期待値算出部２２２３は、算出した振幅期待値λｆ，τを振幅分布算出部２２２４に出力する。
【０１８６】
振幅分布算出部２２２４は、周波数特性量子化部２２２１から量子化音響スペクトログラム［Ａｆ，ｔ］が、振幅期待値算出部２２２３から振幅期待値λｆ，τ、粒子遷移部１３２から粒子毎の粒子情報が、それぞれ入力される。
振幅分布算出部２２２４は、入力された粒子ｉ毎の粒子情報に基づき楽譜位置ｐｉｋを基準とした対応する時刻τ及び周波数ｆ毎に、振幅期待値λｆ，τの量子化音響スペクトログラム［Ａｆ，τ］に対するポアソン分布の関数値Ｐｏｉ（［Ａｆ，τ］｜λｆ，τ）を算出する。ポアソン分布は、式（４０）で表される確率密度関数である。
【０１８７】
【数４０】
0005924968.tif 000041
【０１８８】
式（４０）において、ｋｘは、０又は０よりも大きい整数、λは、０よりも大きい実数である。
算出した関数値は、ポアソン分布の性質により、周波数ｆにおける量子化音響スペクトログラム［Ａｆ，τ］の量子化した振幅が表す度数に対する楽譜スペクトログラムＡ＾ｆ，ｐ（τ）に含まれる振幅の分布関数の値を表す。即ち、算出した関数値は、周波数ｆにおける楽譜スペクトログラムＡ＾ｆ，ｐの振幅が表す度数が、量子化音響スペクトログラム［Ａｆ，τ］の量子化した振幅が表す度数だけ発生する確率値を表す。
振幅分布算出部２２２４は、算出した振幅分布値を重み算出部２２２５に出力する。
【０１８９】
重み算出部２２２５は、入力された振幅分布値を周波数ｆ及び時刻τにわたり累積してＰＤＡを算出する。累積する時刻τの範囲は、粒子毎の観測区間Ｔｋ｛τ｜ｋΔＴ−Ｗ＜τ≦ｋΔＴ｝である。重み算出部２２２５は、算出したＰＤＡの指数関数値を周波数特性重み値ｗｓｐｉ，ｋを算出する。ここで、重み算出部２２２５は、式（４１）を用いて周波数特性重み値ｗｓｐｉ，ｋを算出する。
【０１９０】
【数４１】
0005924968.tif 000042
【０１９１】
重み算出部２２２５は、算出した周波数特性重み値ｗｓｐｉ，ｋを表す周波数特性重み情報を粒子重み算出部１２４に出力する。
【０１９２】
次に、本実施形態に係る楽譜位置推定処理について説明する。本実施形態に係る楽譜位置推定処理は、図３に示す楽譜位置推定処理と同様である。但し、ステップＳ１０６において周波数特性重み算出部２２２は、次に説明する処理を実行する。
図１３は、本実施形態における周波数特性重み算出処理を表すフローチャートである。
（ステップＳ４０１）周波数特性量子化部２２２１は、周波数特性分析部１０３から入力された音響スペクトログラムＡｆ，ｔを予め定めた量子化幅ΔＡで量子化して量子化音響スペクトログラム［Ａｆ，ｔ］を算出する。周波数特性量子化部２２２１は、算出した量子化音響スペクトログラム［Ａｆ，ｔ］を振幅期待値算出部２２２３及び振幅分布算出部２２２４に出力する。その後、ステップＳ４０２に進む。
【０１９３】
（ステップＳ４０２）周波数特性合成部２２２２は、楽譜情報入力部１１２から入力された楽譜情報から楽譜フレームｐ毎の基本周波数ベクトル｛μｐ｝を抽出し、基本周波数ベクトル｛μｐ｝に含まれる基本周波数μｐｌ毎に予め定めた音源モデル（例えば、ＧＭＭ）を用いて調波構造を生成する。周波数特性合成部２２２２は、生成した調波構造を基本周波数μｐｌにわたって累積し、例えば式（３８）を用いて楽譜スペクトログラムＡ＾ｆ，ｐを生成する。
周波数特性合成部２２２２は、生成した楽譜スペクトログラムＡ＾ｆ，ｐを振幅期待値算出部２２２３に出力する。その後、ステップＳ４０３に進む。
【０１９４】
（ステップＳ４０３）振幅期待値算出部２２２３は、粒子遷移部１３２から入力された粒子情報から粒子ｉ毎に楽譜位置情報を抽出する。振幅期待値算出部２２２３は、抽出した楽譜位置情報が表す楽譜位置ｐｉｋを基準とした対応する時刻τ毎に、周波数特性量子化部２２２１から入力された量子化音響スペクトログラム［Ａｆ，τ］の周波数にわたった総和Ａ・τを算出する。振幅期待値算出部２２２３は、周波数特性合成部２２２２から入力された楽譜スペクトログラムＡ＾ｆ，ｐを、対応する楽譜位置ｐに対応する時刻τ毎に算出した総和Ａ・τを乗算して振幅期待値λｆ，τを算出する。振幅期待値算出部２２２３は、算出した振幅期待値λｆ，τを振幅分布算出部２２２４に出力する。その後、ステップＳ４０４に進む。
【０１９５】
（ステップＳ４０４）振幅分布算出部２２２４は、周波数特性量子化部２２２１から量子化音響スペクトログラム［Ａｆ，ｔ］が、振幅期待値算出部２２２３から振幅期待値λｆ，τ、粒子遷移部１３２から粒子毎の粒子情報がそれぞれ入力される。
振幅分布算出部２２２４は、粒子遷移部１３２から入力された粒子ｉ毎の粒子情報に含まれる楽譜位置情報を抽出する。振幅分布算出部２２２４は、抽出した楽譜位置情報が表す楽譜位置ｐｉｋを基準とした対応する時刻τ及び周波数ｆ毎に、振幅期待値λｆ，τの量子化音響スペクトログラム［Ａｆ，τ］に対するポアソン分布Ｐｏｉ（［Ａｆ，τ］｜λｆ，τ）の関数値を振幅分布値として算出する。
振幅分布算出部２２２４は、算出した振幅分布値を重み算出部２２２５に出力する。その後、ステップＳ４０５に進む。
【０１９６】
（ステップＳ４０５）重み算出部２２２５は、振幅分布算出部２２２４から入力された振幅分布値を周波数ｆ及び粒子毎の観測区間Ｔｋ内の時刻τにわたり累積してＰＤＡを算出する。重み算出部２２２５は、算出したＰＤＡの指数関数値を、例えば式（４１）を用いて周波数特性重み値ｗｓｐｉ，ｋを算出する。重み算出部２２２５は、算出した周波数特性重み値ｗｓｐｉ，ｋを表す周波数特性重み情報を粒子重み算出部１２４に出力する。その後、処理を終了する
【０１９７】
以上、説明したように、本実施形態では、第２の周波数特性は少なくとも１つの調波構造を含み、調波構造に含まれる調波成分毎の調波構造への寄与を表す変数の確率分布と、調波構造毎の第２の周波数特性への寄与を表す変数の確率分布に基づいて音響信号の第１の周波数特性についての尤度を最大にするように、楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を算出する。
これにより、音響信号に基づく第１の周波数特性の周波数毎の振幅と楽譜情報における音階に基づく第２の周波数特性の周波数毎の振幅が合致する度合いを鋭敏に検出することができる。そのため、演奏された音楽に対する楽譜の位置をより正確に推定することが可能になる。
【０１９８】
（第３の実施形態）
次に、本発明の第３の実施形態について説明する。第１の実施形態と同一の構成又は処理については、同一の番号を付す。以下、第１の実施形態との差異点を主として説明する。
図１４は、本実施形態に係る楽譜位置推定装置３の構成を示す概略図である。
楽譜位置推定装置３は、楽譜位置推定装置１（図１）と同様に、音響信号入力部１０１、音響特徴量生成部１０２、楽譜情報記憶部１１１、楽譜情報入力部１１２、関連度算出部１２１、楽譜位置探索部１３１、及び演奏情報出力部１４１を含んで構成される。
但し、音響特徴量生成部１０２は、周波数特性分析部１０３及び相関算出部１０４の他にクロマベクトル生成部３０５を備える点が楽譜位置推定装置１とは異なる。周波数特性分析部１０３は、算出した音響スペクトログラムＡｆ，ｔをクロマベクトル生成部３０５に出力する。
【０１９９】
クロマベクトル生成部３０５は、周波数特性分析部１０３から入力された音響スペクトログラムＡｆ，ｔに基づいて、音響信号によるクロマベクトル（以下、音響クロマベクトルと呼ぶ）ｃｔａを生成する。クロマベクトルとは、音階（ｃｈｒｏｍａ）毎の成分のパワーを要素として含むベクトルである。クロマベクトルは、例えば、西洋音楽を構成する１２音階（Ｃ、Ｃ＃、Ｄ、Ｄ＃、Ｅ、Ｆ、Ｆ＃、Ｇ、Ｇ＃、Ａ、Ａ＃、Ｂ）における音階毎の成分の強度を要素とするベクトル（例えば、要素数が１２）である。
クロマベクトル生成部３０５は、音響クロマベクトルｃｔａ＝［ｃｔ，１ａ，ｃｔ，２ａ，…，ｃｔ，１２ａ］の音階ｊ毎の要素ｃｔ，ｊａを、例えば、式（４２）を用いて算出する。
【０２００】
【数４２】
0005924968.tif 000043
【０２０１】
式（４２）において、ｏは、オクターブを表す整数値である。Ｏｃｔｌｏｗは予め設定されたオクターブの下限を表す整数値（下限オクターブ、例えば３）を表す。Ｏｃｔｈｉは、予め設定されたオクターブの上限を表す整数値（上限オクターブ、例えば６）を表す。Ｂｊ，ｏ（ｆ）は、オクターブｏにおける音階ｊの成分を抽出する帯域通過フィルタ（バンドパスフィルタ、ｂａｎｄｐａｓｓ  ｆｉｌｔｅｒ；ＢＰＦ）の入出力振幅周波数特性を表す。ｏ番目のオクターブにおける音階ｊの周波数ｆｊ，ｏにおいて、Ｂｊ，ｏ（ｆ）は、例えば、式（４３）に示されるように最大となり、周波数ｆがゼロ又は無限大に近づくと、ゼロに近似する関数である。
【０２０２】
【数４３】
0005924968.tif 000044
【０２０３】
式（４３）において、ｆｃｅｎｔ、ｆｊ，ｏｃｅｎｔは、それぞれ周波数ｆ、ｆｊ，ｏを対数領域で表したセント値である。セント値ｆｃｅｎｔと周波数ｆは、式（４４）に示す関係がある。
【０２０４】
【数４４】
0005924968.tif 000045
【０２０５】
式（４４）では、例えば、音階Ａ４（ｊ＝１、ｏ＝４）の基本周波数ｆが４４０Ｈｚであることが仮定されている。
クロマベクトル生成部３０５は、生成した音響クロマベクトルｃｔａを、後述するクロマベクトル重み算出部３２５に出力する。
【０２０６】
楽譜情報記憶部１１１に記憶されている楽譜情報には、楽曲を構成する各コードを構成する音階を要素として含むクロマベクトル（以下、楽譜クロマベクトルと呼ぶ）ｃｐｓが、基本周波数ベクトルμｐと対応付けられて含まれている。楽譜クロマベクトルの要素数は、１２である。つまり、楽譜クロマベクトルを構成する音階は、楽曲を表す音階のオクターブが区別されていない点で上述のコードベクトルを構成する音階と異なる。また、楽譜クロマベクトルｃｐｓを構成する要素ｃｐ，ｊｓは、０または１の値をとる。０とは、楽譜フレームｐにおいて音階ｊが存在しないことを表し、１とは、音階ｊが存在することを表す。例えば、音階Ａ４、Ｃ４、及びＥ４を要素値として含む楽譜クロマベクトルは、［１，０，０，１，０，０，０，１，０，０，０，０］となる。
楽譜情報入力部１１２は、楽譜情報記憶部１１１から読み出した楽譜情報に含まれる楽譜クロマベクトルｃｐｓをクロマベクトル重み算出部３２５に出力する。
【０２０７】
関連度算出部１２１は、周波数特性重み算出部１２２及び拍間隔重み算出部１２３の他にクロマベクトル重み算出部３２５を備え、粒子重み算出部１２４の代わりに粒子重み算出部３２４を備える。
なお、粒子遷移部１３２は、粒子ｉ毎の粒子情報をクロマベクトル重み算出部３２５に出力する。粒子遷移部１３２は、粒子ｉ毎の粒子分布値、状態遷移確率を粒子重み算出部３２４に出力する。
【０２０８】
クロマベクトル重み算出部３２５は、クロマベクトル生成部３０５から入力された音響クロマベクトルｃｔａを、時刻ｔ毎にノルムが１となるように各要素値を規格化する。クロマベクトル重み算出部３２５は、楽譜情報入力部１１２から入力された楽譜クロマベクトルｃｐｓを、楽譜フレームｐ毎にノルムが１となるように各要素値を規格化する。
クロマベクトル重み算出部３２５は、粒子遷移部１３２から入力された粒子情報から粒子ｉ毎に楽譜位置情報ｐｉｋを抽出する。クロマベクトル重み算出部３２５は、要素値を規格化した音響クロマベクトルｃτａと楽譜クロマベクトルｃｐ（τ）’ｓとの内積を、例えば式（４５）を用いて平均して、粒子ｉ毎のクロマベクトル重みｗｉ，ｋｃｈを算出する。内積は、ベクトル間の類似性を表す指標値である。
【０２０９】
【数４５】
0005924968.tif 000046
【０２１０】
式（４５）において、Ｔｋは、抽出した楽譜位置ｐｉｋを基準（時刻τがｋΔＴ）とした粒子毎の観測区間Ｔｋ｛τ｜ｋΔＴ−Ｗ＜τ≦ｋΔＴ｝を表す。
クロマベクトル重み算出部３２５は、算出した粒子ｉ毎のクロマベクトル重みｗｉ，ｋｃｈを表すクロマベクトル重み情報を粒子重み算出部３２４に出力する。
【０２１１】
粒子重み算出部３２４は、周波数特性重み算出部１２２から周波数特性重み情報が、拍間隔重み算出部１２３から拍間隔重み情報が、クロマベクトル重み算出部３２５からクロマベクトル重み情報が、粒子遷移部１３２から粒子ｉ毎の粒子分布値、状態遷移確率が入力される。
粒子重み算出部３２４は、式（４６）に示すように、周波数特性重み情報が表す周波数特性重み値、拍間隔重み情報が表す拍間隔重み値、クロマベクトル重み情報が表すクロマベクトル重み値、状態遷移確率を乗算し、さらに状態遷移確率で除算して、粒子ｉ毎の重み値ｗｉ，ｋを算出する。
【０２１２】
【数４６】
0005924968.tif 000047
【０２１３】
粒子重み算出部３２４は、算出した重み値ｗｉ，ｋを表す粒子重み情報を生成する。
粒子重み算出部３２４は、粒子遷移部１３２から入力された粒子ｉ毎の粒子情報に含まれる粒子重み情報を、生成した粒子重み情報に置き換えることで粒子情報を更新する。
粒子重み算出部３２４は、更新した粒子情報を再標本化部１３３、楽譜位置算出部１４２及び拍間隔算出部１４３に出力する。
なお、本実施形態に係る楽譜位置推定装置３は、周波数特性重み算出部１２２の代わりに周波数特性重み算出部２２２を備えてもよい。
【０２１４】
次に、本実施形態に係る楽譜位置推定処理について説明する。
図１５は、本実施形態に係る楽譜位置推定処理を表すフローチャートである。
図１５が示す楽譜位置推定処理は、ステップＳ１０１、Ｓ１０２−Ｓ１０５、Ｓ１０６、Ｓ１０７、Ｓ１０９−Ｓ１１７を有する点で、図３が示す楽譜位置推定処理と共通する。但し、ステップＳ１０１とステップＳ１０２の間でステップＳ４０１を実行し、ステップＳ１０５とステップＳ１０６の間でステップＳ４０２を実行し、ステップＳ１０８の代わりにステップＳ４０３、Ｓ４０４を実行する点が、図３が示す楽譜位置推定処理とは異なる。以下、処理が異なるステップについて説明する。
【０２１５】
（ステップＳ４０１）音源位置推定装置３では、楽曲を構成する各コードを構成する音階を要素として含む楽譜情報に基づいて楽譜クロマベクトルｃｐｓを予め生成しておく。音源位置推定装置３では、基本周波数ベクトルμｐと対応付けられ楽譜情報の一部として予め楽譜情報記憶部１１１に記憶させておく。これにより、楽譜情報入力部１１２が、楽譜情報記憶部１１１から読み出した楽譜情報から楽譜クロマベクトルｃｐｓをクロマベクトル重み算出部３２５に出力できるようにする。その後、ステップＳ１０２に進む。
【０２１６】
（ステップＳ４０２）クロマベクトル生成部３０５は、周波数特性分析部１０３から入力された音響スペクトログラムＡｆ，ｔに基づいて、例えば式（４２）を用いて算出される要素を有する音響クロマベクトルｃｔａを生成する。
クロマベクトル生成部３０５は、生成した音響クロマベクトルｃｔａをクロマベクトル重み算出部３２５に出力する。その後、ステップＳ１０６に進む。
【０２１７】
（ステップＳ４０３）クロマベクトル重み算出部３２５は、クロマベクトル生成部３０５から入力された音響クロマベクトルｃｔａを、時刻ｔ毎にノルムが１となるように各要素値を規格化する。クロマベクトル重み算出部３２５は、楽譜情報入力部１１２から入力された楽譜クロマベクトルｃｐｓを、楽譜フレームｐ毎にノルムが１となるように各要素値を規格化する。
クロマベクトル重み算出部３２５は、粒子遷移部１３２から入力された粒子情報から粒子ｉ毎に楽譜位置情報を抽出する。クロマベクトル重み算出部３２５は、要素値を規格化した音響クロマベクトルｃτａと楽譜クロマベクトルｃｐ（τ）ｓとの内積を、例えば式（４５）を用いて抽出した楽譜位置情報が表す楽譜位置を基準とした観測区間Ｔｋで平均する。これにより、クロマベクトル重み算出部３２５は、粒子ｉ毎のクロマベクトル重みｗｉ，ｋｃｈを算出する。クロマベクトル重み算出部３２５は、算出した粒子ｉ毎のクロマベクトル重みｗｉ，ｋｃｈを表すクロマベクトル重み情報を粒子重み算出部１２４に出力する。その後、ステップＳ４０４に進む。
【０２１８】
（ステップＳ４０４）粒子重み算出部３２４は、周波数特性重み算出部１２２から周波数特性重み情報が、拍間隔重み算出部１２３から拍間隔重み情報が、クロマベクトル重み算出部３２５からクロマベクトル重み情報が、粒子遷移部１３２から粒子ｉ毎の粒子分布値、状態遷移確率が入力される。
粒子重み算出部３２４は、周波数特性重み情報が表す周波数特性重み値、拍間隔重み情報が表す拍間隔重み値、クロマベクトル重み情報が表すクロマベクトル重み値、状態遷移確率を乗算し、さらに状態遷移確率で除算して、粒子ｉ毎の重み値ｗｉ，ｋを算出する。
粒子重み算出部３２４は、算出した重み値ｗｉ，ｋを表す粒子重み情報を生成する。
粒子重み算出部３２４は、粒子遷移部１３２から入力された粒子ｉ毎の粒子情報に含まれる粒子重み情報を、生成した粒子重み情報に置き換えることで粒子情報を更新する。
粒子重み算出部３２４は、更新した粒子情報を楽譜位置探索部１３１の再標本化部１３３、演奏情報出力部１４１の楽譜位置算出部１４２及び拍間隔算出部１４３に出力する。その後、ステップＳ１０９に進む。以下、図３に示す楽譜位置推定処理と同様にステップＳ１１０−１１７を実行する。
なお、本実施形態に係る楽譜位置推定処理では、ステップＳ１０６において周波数特性重み算出部１２２が図１０に示す周波数特性重み算出処理を実行する代わりに、周波数特性重み算出部２２２が図１３に示す周波数特性重み算出処理を実行してもよい。
【０２１９】
次に、本実施形態に係る楽譜位置推定装置２が推定した楽譜位置の検証例について説明する。検証においては特に断らない限り、周波数の上限ｆｍａｘを６０００Ｈｚ、テンポ窓長６０／ｂθを１５（ｂｐｍ）とした。
【０２２０】
図１６は、楽譜位置の推定誤差の一例を表す図である。
図１６において、縦軸は、推定誤差、横軸は、楽曲の番号を表す。図１５において、○印は、本実施形態による平均推定誤差、×印は、従来の楽譜位置推定システムであるＡｎｔｅｓｃｏｆｏによる平均推定誤差を表す（例えば、Ｃｏｎｔ．Ａ．：Ａ  Ｃｏｕｐｌｅｄ  Ｄｕｒａｔｉｏｎ−Ｆｏｃｕｓｅｄ  Ａｒｃｈｉｔｅｃｔｕｒｅ  ｆｏｒ  Ｒｅａｌｔｉｍｅ  Ｍｕｓｉｃ  ｔｏ  Ｓｃｏｒｅ  Ａｌｉｇｎｍｅｎｔ，ＩＥＥＥ  Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｎ  Ｐａｔｔｅｒｎ  Ａｎａｌｙｓｉｓ  ａｎｄ  Ｍａｃｈｉｎｅ  Ｉｎｔｅｌｌｉｇｅｎｃｅ，  Ｖｏｌ．３２，Ｎｏ．６，ｐｐ．９７４−９８７（２０１０））。Ａｎｔｅｓｃｏｆｏは、ＨＳＭＭ（Ｈｉｄｄｅｎ  Ｓｅｍｉ−Ｍａｒｋｏｖ  Ｍｏｄｅｌ；隠れセミマルコフモデル）に基づいて楽譜位置を推定する処理を行う。○印、×印を中心として上下に延びるバーの上限と下限は、それぞれ平均推定誤差から標準偏差だけ大きい推定誤差と、平均推定誤差から標準偏差だけ小さい推定誤差を表す。即ち、バーは、推定誤差の信頼区間を表す。楽曲１−２０は、各々テンポ及び演奏する楽器が異なる楽曲である。推定誤差は、式（１１）を用いて算出した値である。
【０２２１】
図１６によれば、本実施形態では、ほとんどの楽曲についてＡｎｔｅｓｃｏｆｏよりも推定誤差の平均値及び標準偏差が減少する。本実施形態の楽曲にわたる平均推定誤差の絶対値の平均値は、Ａｎｔｅｓｃｏｆｏの平均推定誤差の絶対値の平均値よりも６９％減少する。即ち、本実施形態の方が、Ａｎｔｅｓｃｏｆｏよりも推定精度が優れる。一部の楽曲、例えば楽曲６、１１では、本実施形態の平均推定誤差の絶対値の方が、Ａｎｔｅｓｃｏｆｏの平均推定誤差の絶対値よりも大きい。しかし、本実施形態の標準偏差の方が、Ａｎｔｅｓｃｏｆｏの標準偏差よりも小さい。即ち、本実施形態の方が、より頑健に楽譜位置を推定できることが示されている。
【０２２２】
図１７は、楽譜位置の推定誤差のその他の例を表す図である。
図１７において、横軸と縦軸の関係は、図１６と同様である。図１７において、○印、×印、△印は、テンポ窓長６０／ｂθが、それぞれ５、１５、３０（ｂｐｍ）である場合の平均推定誤差を表す。○印、×印、△印を中心として上下に延びるバーは、平均推定誤差を中心とし、その長さが標準偏差の２倍である信頼区間を表す。
図１７によれば、ほとんどの楽曲について、テンポ窓長が小さくなるほど平均推定誤差の絶対値が小さくなり、標準偏差も小さくなる傾向がある。即ち、テンポ窓長が小さいほど楽譜位置の推定精度が優れ、楽譜位置を安定して推定できることが示されている。
【０２２３】
図１８は、楽譜位置の推定誤差のその他の例を表す図である。
図１８において、縦軸は、推定誤差の絶対値を表し、横軸は、楽曲の番号を表す。棒グラフは、推定処理又はシステム毎の平均推定誤差を表す。最も右側の網掛けの棒グラフは、第１の実施形態においてＬＨＡ法を用いて周波数特性重み値ｗｉ，ｋｓｐを算出する処理を表す。右から２番目の縦縞で塗りつぶした棒グラフは、第２の実施形態においてＰＤＡを用いて周波数特性重み値ｗｉ，ｋｓｐを算出する処理を表す。左から２番目の横縞で塗りつぶした棒グラフは、第２の実施形態の構成においてＰＤＡの代わりにカルバック・ライブラー情報量（ＫＬ−ｄｉｖ）を用いて周波数特性重み値ｗｉ，ｋｓｐを算出する処理を表す。最も左側の斜めの縞模様で塗りつぶした棒グラフは、Ａｎｔｅｓｃｏｆｏが採用するＨＳＭＭを表す。また、各棒グラフの頂上から上に伸びるエラーバーの長さは、標準偏差を表す。
図１８によれば、楽曲１−１０のいずれにおいても、ＬＨＡ、ＰＤＡ、ＫＬ−ｄｉｖ、ＨＳＭＭの順で、平均推定誤差、標準偏差が、ともに小さくなる。即ち、この順序で、楽譜位置の推定精度が優れ、安定して楽譜位置を推定できることが示される。
【０２２４】
図１９は、音響信号と楽譜情報との関連性の一例を表す図である。
図１９において、縦軸は対数尤度を、横軸は楽譜位置偏差を表す。図１９が表す対数尤度は、ある区間の楽譜情報と、ある観測区間の音響信号との間の周波数特性重み値ｗｉ，ｋｓｐの対数値に相当する値である。即ち、対数尤度は、音響信号と楽譜情報との関連性を表す指標値である。但し、図１９に示す対数尤度は、計算精度を確保するために規格化されていない値である。楽譜位置偏差とは、楽譜情報の区間と対応する観測区間を基準とした、音響信号の観測区間のずれを、拍単位で表す値である。一点破線は、第１の実施形態においてＬＨＡ法を用いて周波数特性重み値ｗｉ，ｋｓｐを算出する処理を表す。破線は、第２の実施形態においてＰＤＡを用いて周波数特性重み値ｗｉ，ｋｓｐを算出する処理を表す。実線は、ＫＬ−ｄｉｖを用いて周波数特性重み値ｗｉ，ｋｓｐを算出する処理を表す。
表示形態が異なる３本の太線は、各々楽譜位置偏差毎の対数尤度を表し、３本の細線は、各々対数尤度の最大値を表す。いずれも、楽譜位置偏差が０の場合、対数尤度が最大になる。即ち、楽譜位置偏差が０とは、楽譜位置と音響信号の時刻の対応が取れている点であることを表す。
【０２２５】
図１９において、楽譜位置偏差による対数尤度の変化は、ＬＨＡ、ＰＤＡ、ＫＬ−ｄｉｖの順で著しい。とりわけ、第１の実施形態においてＬＨＡを用いた場合や、第２の実施形態においてＰＤＡを用いた場合の対数尤度の鋭敏な変化は、ＫＬ−ｄｉｖを用いた場合の緩慢な対数尤度の変化と対照的である。このことは、上述の実施形態において音源モデルとしてＬＨＡ又はＰＤＡを用い、最尤推定を行うことで楽譜位置が正確に推定できることを示す。
【０２２６】
以上、説明したように、上述した実施形態では、入力された音響信号に基づく第１の周波数特性と楽譜情報が表す楽譜の位置毎の音階に基づく第２の周波数特性と前記第１の周波数特性の関連度を表す重み値を、第２の周波数特性に含まれる成分毎の第１の周波数特性の確率分布に基づいて、楽譜の位置毎に算出する。また、上述した実施形態では、算出した重み値に基づいて音響信号に対応する楽譜の位置を探索する。
これにより、音響信号に基づく第１の周波数特性と、楽譜情報における音階に基づく第２の周波数特性が完全に合致していなくとも、両者の関連性を検出することができる。そのため、演奏された音楽に対する楽譜の位置をより頑健に推定することが可能になる。
【０２２７】
また、上述した実施形態では、音階が変化する楽譜の位置に対応する時刻を含む観測区間における音響信号の自己相関に基づいて拍間隔を定め、定めた拍間隔に基づいて探索する楽譜の位置を更新する。
これにより、楽譜情報における音階が変化する楽譜の位置と、音響信号の振幅が変化する時刻が対応付けられる。そのため、周波数特性の周期性を表す拍間隔を正確に検知でき、ひいては演奏された音楽に対する楽譜の位置をより正確に推定することができる。
【０２２８】
上述では、楽譜位置探索部１３１が、粒子フィルタリング法を用いて楽譜位置を探索する例について説明したが、上述した実施形態ではこれには限られない。上述した実施形態では、楽譜位置毎の楽譜情報と音響信号との関連度を表す重み情報に基づいて楽譜位置を定めることができる処理であれば、いかなる処理、例えば、ある探索区間内を候補となる楽譜位置を走査しながら重み情報を算出して、算出した重み情報に基づいて楽譜位置を探索してもよい。
【０２２９】
なお、上述した実施形態における楽譜位置推定装置１、２又は３の一部、例えば、音響特徴量生成部１０２、関連度算出部１２１、楽譜位置探索部１３１、及び演奏情報出力部１４１をコンピュータで実現するようにしても良い。その場合、この制御機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現しても良い。なお、ここでいう「コンピュータシステム」とは、楽譜位置推定装置１、２又は３に内蔵されたコンピュータシステムであって、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含んでも良い。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。
また、上述した実施形態における楽譜位置推定装置１、２又は３の一部、または全部を、ＬＳＩ（Ｌａｒｇｅ  Ｓｃａｌｅ  Ｉｎｔｅｇｒａｔｉｏｎ）等の集積回路として実現しても良い。楽譜位置推定装置１又は２の各機能ブロックは個別にプロセッサ化してもよいし、一部、または全部を集積してプロセッサ化しても良い。また、集積回路化の手法はＬＳＩに限らず専用回路、または汎用プロセッサで実現しても良い。また、半導体技術の進歩によりＬＳＩに代替する集積回路化の技術が出現した場合、当該技術による集積回路を用いても良い。
【０２３０】
以上、図面を参照してこの発明の一実施形態について詳しく説明してきたが、具体的な構成は上述のものに限られることはなく、この発明の要旨を逸脱しない範囲内において様々な設計変更等をすることが可能である。
【符号の説明】
【０２３１】
１、２…楽譜位置推定装置、１０１…音響信号入力部、１０２…音響特徴量生成部、
１０３…周波数特性分析部、１０４…相関算出部、３０５…クロマベクトル生成部
１１１…楽譜情報記憶部、１１２…楽譜情報入力部、
１２１…関連度算出部、１２２、２２２…周波数特性重み算出部、
１２２１…データ入力部、１２２２…第１変分パラメータ算出部、
１２２３…第２変分パラメータ算出部、１２２４…変分下限値算出部、
１２２５…収束判定部、１２２６…重み算出部、
２２２１…周波数特性量子化部、２２２２…周波数特性合成部、
２２２３…振幅期待値算出部、２２２４…振幅分布算出部、２２２５…重み算出部、
１２３…拍間隔重み算出部、１２４、３２４…粒子重み算出部、
３２５…クロマベクトル重み算出部
１３１…楽譜位置探索部、１３２…粒子遷移部、１３３…再標本化部、
１４１…演奏情報出力部、１４２…楽譜位置算出部、１４３…拍間隔算出部、
１４４…信頼度判定部、１４５…楽譜位置出力部、１４６…拍間隔出力部
【図１】
0005924968.tif 000048
【図２】
0005924968.tif 000049
【図３】
0005924968.tif 000050
【図５】
0005924968.tif 000051
【図６】
0005924968.tif 000052
【図７】
0005924968.tif 000053
【図８】
0005924968.tif 000054
【図９】
0005924968.tif 000055
【図１０】
0005924968.tif 000056
【図１１】
0005924968.tif 000057
【図１２】
0005924968.tif 000058
【図１３】
0005924968.tif 000059
【図１４】
0005924968.tif 000060
【図１５】
0005924968.tif 000061
【図１６】
0005924968.tif 000062
【図１７】
0005924968.tif 000063
【図１８】
0005924968.tif 000064
【図１９】
0005924968.tif 000065
【図４】
0005924968.tif 000066