"入力された音響信号の特徴量を生成する音響特徴量生成部と、前記音響信号の特徴量と楽譜情報の特徴量に基づき重み係数を算出し、前記重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定する楽譜位置推定部と、を備えることを特徴とする楽譜位置推定装置。前記仮想テンポと第1の乱数に基づいて前記仮想楽譜位置を更新し、第2の乱数に基づいて前記仮想テンポを更新する状態遷移部を備えることを特徴とする請求項1に記載の楽譜位置推定装置。前記特徴量は、音階ごとの強度を複数個含むクロマベクトル又は周波数特性を含み、前記楽譜位置推定部は、前記重み係数を、前記音響信号によるクロマベクトルと前記楽譜情報によるクロマベクトルとの類似度を示すクロマベクトル重み係数、前記前記音響信号の周波数特性と前記楽譜情報の周波数特性との類似度を示す周波数特性重み係数、又は前記音響信号の周波数特性の自己相関に基づくテンポ重み係数を用いて算出すること、を特徴とする請求項1又は2に記載の楽譜位置推定装置。前記仮想楽譜位置に基づいて前記推定した楽譜位置の信頼度を算出し、前記信頼度に基づいて前記推定した楽譜位置を出力するか否か判断する信頼度算出部を備えること、を特徴とする請求項1乃至3のうちいずれか1項に記載の楽譜位置推定装置。前記推定した楽譜位置の信頼度を算出する信頼度算出部を備え、前記楽譜位置推定部は、前記算出された信頼度に基づき前記重み係数を算出する前記特徴量の区間を延長するか否か判断すること、を特徴とする請求項1乃至3のうちいずれか1項に記載の楽譜位置推定装置。楽譜位置推定装置における楽譜位置推定方法において、前記楽譜位置推定装置が、入力された音響信号の特徴量を所定の区間ごとに生成する第1の過程と、前記楽譜位置推定装置が、前記音響信号の特徴量と楽譜情報の特徴量に基づき重み係数を算出し、前記重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定する第2の過程とを有すること、を特徴とする楽譜位置推定方法。楽譜位置推定装置のコンピュータに、入力された音響信号の特徴量を生成する手順と、前記音響信号の特徴量と楽譜情報の特徴量に基づき重み係数を算出し、前記重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定する手順と、を実行させるための楽譜位置推定プログラム。しかし、特許文献2記載の方法では、収音した楽曲が楽譜上のどの位置以下、楽譜位置ScorePositionと呼ぶにあるかを認識できない。つまり、特許文献2記載の方法は、演奏されている曲のビート拍;beat時間やテンポの抽出に失敗することがあるという問題点があった。本発明は上記の点に鑑みてなされたものであり、演奏された音楽に対する楽譜位置を推定できる楽譜位置推定装置、楽譜位置推定方法、及び楽譜位置推定プログラムを提供することを課題としている。1本発明は上記の課題を解決するためになされたものであり、本発明の一態様は、入力された音響信号の特徴量を生成する音響特徴量生成部と、前記音響信号の特徴量と楽譜情報の特徴量に基づき重み係数を算出し、前記重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定する楽譜位置推定部とを備えることを特徴とする楽譜位置推定装置である。2本発明の他の態様は、前記仮想テンポと第1の乱数に基づいて前記仮想楽譜位置を前記所定の区間ごとに更新し、第2の乱数に基づいて前記仮想テンポを前記所定の区間ごとに更新する状態遷移部を備えることを特徴とする1の楽譜位置推定装置である。3本発明の他の態様は、前記特徴量は、音階ごとの強度を複数個含むクロマベクトル又は周波数特性を含み、前記楽譜位置推定部は、前記重み係数を、所定の区間における前記音響信号によるクロマベクトルと前記楽譜情報によるクロマベクトルとの類似度を示すクロマベクトル重み係数、前記所定の区間における前記音響信号の周波数特性と前記楽譜情報の周波数特性との類似度を示す周波数特性重み係数、又は前記所定の区間における前記音響信号の周波数特性の自己相関に基づくテンポ重み係数を用いて算出することを特徴とする1又は2の楽譜位置推定装置である。4本発明の他の態様は、前記仮想楽譜位置に基づいて前記推定した楽譜位置の信頼度を算出し、前記信頼度に基づいて前記推定した楽譜位置を出力するか否か判断する信頼度算出部を備えることを特徴とする1、2、又は3のいずれかの楽譜位置推定装置である。5本発明の他の態様は、前記推定した楽譜位置の信頼度を算出する信頼度算出部を備え、前記楽譜位置推定部は、前記算出された信頼度に基づき前記重み係数を算出する前記特徴量の区間を延長するか否か判断することを特徴とする1、2、又は3のいずれかの楽譜位置推定装置である。1に記載した態様によれば、入力された音響信号の特徴量を生成し、音響信号の特徴量と楽譜情報の特徴量に基づき重み係数を算出し、重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定するので、演奏された音楽に対する楽譜位置を精度よく推定することが可能になる。2に記載した態様によれば、仮想テンポと第1の乱数に基づいて仮想楽譜位置を更新し、第2の乱数に基づいて仮想テンポを更新するため、演奏された音楽のテンポが変動しても楽譜位置を頑健に推定することが可能になる。3に記載した態様によれば、クロマベクトル又は周波数特性を含む特徴量に基づき音響信号と楽譜情報との類似度を示す重み係数を算出するため、音響信号と楽譜情報が類似する仮想楽譜位置及び仮想テンポが重視される。そのため、当該発明は演奏された音楽に対する楽譜位置を精度よく推定することが可能になる。4に記載した態様によれば、推定した楽譜位置の信頼度に基づいて推定した楽譜位置を出力するか否か判断するため、出力された楽譜位置の信頼度を確保することが可能になる。5に記載した態様によれば、算出された信頼度に基づき重み係数を算出する特徴量の区間を延長するか否か判断するため、推定した楽譜位置の信頼度が確保される特徴量を含むように区間を延長して、出力された楽譜位置の信頼度を確保することが可能になる。第1の実施形態以下、図面を参照しながら本発明の第1の実施形態について説明する。図1は、本実施形態に係る楽譜位置推定装置1の構成を示す概略図である。楽譜位置scoreposition推定装置1は、音響信号入力部101、音響特徴量生成部102、楽譜情報記憶部103、楽譜情報入力部104、楽譜特徴量生成部105、楽譜位置推定部106、信頼度confidence算出部107、楽譜位置出力部108、及びテンポ出力部109を含んで構成される。音響信号入力部101は、演奏された音楽による音波をディジタル音響信号に変換し、変換したディジタル音響信号を音響特徴量生成部102に出力する。音響信号入力部101は、例えばマイクロホンとアナログディジタルAnalog-to-Digital;A/D変換部図示せずとを含んで構成される。このマイクロホンは、人間が聴取することができる周波数帯域例えば、20Hz-20kHzの音波を受信してアナログ音響信号を生成し、生成したアナログ音響信号をA/D変換部に出力する。A/D変換部は、マイクロホンから入力されたアナログ音響信号をディジタル音響信号に変換し、変換したディジタル音響信号を音響特徴量生成部102に出力する。ここで、A/D変換部は、入力されたアナログ音響信号を、例えば、サンプリング周波数44.1kHz、振幅を16ビットの2進バイナリ;binaryデータにパルス符号変調PulseCodeModulation;PCMすることで、量子化されたディジタル音響信号に変換する。音響特徴量生成部102は、音響信号入力部101から入力されたディジタル音響信号から、その特徴量音響特徴量を生成し、生成した音響特徴量を楽譜位置推定部106に出力する。具体的には、音響特徴量生成部102は、音響特徴量として、例えば、スペクトログラムspectrogram及びクロマベクトルchromavectorを生成する。音響特徴量生成部102は、周波数特性分析部1021及びクロマベクトル生成部1022を含んで構成される。周波数特性分析部1021は、音響信号入力部101から入力されたディジタル音響信号を時間領域信号から周波数領域信号に変換する。ここで、周波数特性分析部1021は、例えば、2048個のサンプルからなるフレームごとにディジタル音響信号を高速フーリエ変換FastFourierTransform;FFTして周波数領域信号に変換する。このフレームは、FFTを行うディジタル音響信号の範囲である。周波数特性分析部1021は、予め設定された時刻t毎例えば、10ms間隔に新たな信号のサンプル例えば、サンプリング周波数が44.1kHzの場合、441サンプルが含まれるようにフレームをシフトさせる。周波数特性分析部1021は、フレーム時刻t毎に変換した周波数領域信号の振幅絶対値をスペクトログラム以下、音響スペクトログラムと呼ぶXf,tとして算出し、算出した音響スペクトログラムXf,tをクロマベクトル生成部1022及び楽譜位置推定部106に出力する。なお、fは周波数を表す。クロマベクトル生成部1022は、周波数特性分析部1021から入力された音響スペクトログラムXf,tに基づいて音響信号によるクロマベクトル以下、音響クロマベクトルと呼ぶctaを算出する。クロマベクトルとは、音階chroma毎の成分のパワーを要素として含むベクトルである。クロマベクトルは、例えば、西洋音楽を構成する12音階C、C、D、D、E、F、F、G、G、A、A、Bにおける音階毎の成分の強度を要素とするベクトル例えば、要素数が12である。クロマベクトル生成部1022は、音響クロマベクトルcta=[cta1,cta2,...,cta12]の音階n毎の要素ctanを、例えば、次式を用いて算出する。式1において、OctLは予め設定されたオクターブの下限下限オクターブを表す。OctHは、予め設定されたオクターブの上限上限オクターブを表す。BPFn,hfは、h番目のオクターブにおける音階nの成分を抽出する帯域通過フィルタバンドパスフィルタ、BandPassFilter;BPFの入出力特性を表す。h番目のオクターブにおける音階nの周波数fn,hにおいて、その値は最大となり、周波数fがゼロ又は無限大に近づくと、ゼロに近似する関数である。また、|Xf,t|2は、パワースペクトルを示す。楽譜位置推定装置1は、このクロマベクトルを用いることで、オクターブが異なっても音階毎の強度を抽出でき、音響信号の調波構造harmonicstructureを無視することができる。よって、楽譜位置推定装置1は、旋律を表すために重要な音階の計時変化を抽出でき、旋律を表すために余り重要ではないオクターブの差異を除外することができる。クロマベクトル生成部1022は、算出した音響クロマベクトルctaを楽譜位置推定部106に出力する。楽譜情報記憶部103は、楽譜情報を楽曲毎に記憶する。楽譜情報は、楽譜フレームscoreframekkは1又は1より大きい整数毎に、楽曲を構成する音階を示す音階情報nkrを要素とするベクトルnk=[nk1,nk2,...,nkrk]Tで表される。ここで、rは1又は1より大きくrk又はrkよりも小さい整数である。rkは、楽譜フレームkにおける音階の数を表す。Tは、ベクトルの転置を表す。楽譜フレームとは、1つの楽曲の楽譜情報において楽曲を時間分割する単位時間である。例えば、四分音符の長さが12フレームである場合、楽譜情報の時間分解能resolutionは十六分音符の三分の一となる。例えば、楽譜フレーム2において、音階C、E及びGを表す楽譜情報n2は、要素数が3個のベクトル[C,E,G]となる。楽譜情報入力部104は、楽譜情報記憶部103から、処理の対象とする楽曲の楽譜情報nkを読み出す。楽譜情報入力部104は、読み出した楽譜情報nkを楽譜特徴量生成部105に出力する。楽譜特徴量生成部105は、楽譜情報入力部104から入力された楽譜情報nkに基づいて、その特徴量以下、楽譜特徴量と呼ぶを生成し、生成した楽譜特徴量を楽譜位置推定部106に出力する。楽譜特徴量生成部105は、楽譜特徴量として、例えば、スペクトログラム及びクロマベクトルを生成する。ここで、楽譜特徴量生成部105が生成したスペクトログラム及びクロマベクトルを、それぞれ楽譜スペクトログラム及び楽譜クロマベクトルと呼ぶ。楽譜特徴量生成部105は、周波数特性合成部1051及びクロマベクトル生成部1052を含んで構成される。周波数特性合成部1051は、楽譜情報入力部から入力された楽譜情報nkに基づき、楽譜スペクトログラムX’f,kを生成する。楽譜スペクトログラムX’f,kとは、楽譜情報nkに基づき演奏されるべき音楽を示す音響信号の周波数特性である。周波数特性合成部1051は、例えば次式のように調和混合正規分布モデルharmonicgaussianmixturemodel;harmonicGMMに基づいて楽譜スペクトログラムX’f,kを生成する。式2の右辺第1項において、gは、調波インデックスharmonicindex即ち、基本周波数fundamentalfrequencyFの整数倍integermultipleの周波数の成分を示すインデックスである。g=1は、基本波fundamentaltoneを示す。Gは、調波の数例えば10を示す。hgは、第g高調波g-thharmonicの振幅例えば、0.2のg乗を表す。Nf;gFnk,σ2は、パワーの周波数特性が、中心周波数gFnk、分散σ2の正規分布normaldistribution,gaussiandistribution関数を表す。即ち、式2の右辺第1項は純粋な調和混合正規分布モデルに基づき、各音階を構成する調波成分harmonicを線形結合した成分を示す。なお、分散σ2は、例えば0.8である。式2の右辺第2項において、Cfは、ピンクノイズpinknoise成分のパワーpowerの周波数特性を示す。ピンクノイズとは、オクターブ帯域毎のパワーが等しい周波数特性を有しランダムな波形を表す信号系列である。式2の右辺において、Cfを加算するのは、右辺第1項の成分が表すスペクトログラムを処理すると生じうるゼロ除算zerodivideを避けるためである。Cfは、例えば次式3で示される周波数特性を示す。Aは、式2の右辺第2項の成分について右辺第1項の成分とのパワー比例えば5を示す定数である。αは、パワーの減少率を示す定数である。αは、例えば、log10CfΔf/Cf=-0.6となる値である。但し、Δfは、帯域幅1kHzに対応する周波数領域信号を構成する周波数変数frequencybinの数である。周波数特性合成部1051は、生成した楽譜スペクトログラムX’f,kを楽譜位置推定部106に出力する。クロマベクトル生成部1052は、楽譜情報入力部104から入力された楽譜情報nkに基づき楽譜クロマベクトルcksを生成する。クロマベクトル生成部1052は、例えば、楽譜情報nkに含まれる音階nに対応する要素cksnを1とし、楽譜情報nkに含まれない音階nに対応する要素cksnをゼロとした楽譜クロマベクトルを生成する。例えば、n2が[C,E,G]であるとき、クロマベクトル生成部1052は、楽譜クロマベクトルc2s=[1,0,0,0,1,0,0,1,0,0,0,0]を生成する。クロマベクトル生成部1052は、生成した楽譜クロマベクトルcksを楽譜位置推定部106に出力する。楽譜位置推定部106は、重みweight算出部1061、再標本化resampling部1062、状態遷移statetransition部1063、楽譜位置算出部1064及びテンポ算出部1065を含んで構成される。重み算出部1061は、音響特徴量生成部102から入力された音響特徴量と楽譜特徴量に基づいて後述する粒子毎に重み係数を算出する。重み算出部1061は、算出した粒子毎の重み係数を再標本化部1062に出力する。具体的には、重み算出部1061は、クロマベクトル重みchromavectorweight算出部10611、周波数特性重みspectrogramweight算出部10612、テンポ重みtempoweight算出部10613、及び粒子重み算出部10614を含んで構成される。クロマベクトル重み算出部10611は、クロマベクトル生成部1022から入力された音響クロマベクトルctaとクロマベクトル生成部1052から入力された楽譜クロマベクトルcktsに基づいてクロマベクトル重み係数wchを算出する。クロマベクトル重み係数wchは、音響クロマベクトルctaと楽譜クロマベクトルcksとの類似度又は相関性を示す尺度である。クロマベクトル重み算出部10611は、例えば、以下の式を用いてクロマベクトル重み係数wchを算出する。式4において、τは時刻を表す。Lは予め設定された時間を表す。tは最新のフレーム時刻を表す。kτは、観測時刻τにおける楽譜フレームkτである。kτと観測時刻τは次式5の関係がある。式5において、floor...は、floor関数floorfunction;床関数を表す。即ち、floor...は、実数...を、その整数部に切り捨てる関数である。k’は、仮想楽譜位置である。b’は、仮想拍間隔beatintervalである。k’及びb’については、後述する。即ち、式4は、クロマベクトル重み係数wchが、音響クロマベクトルctaと楽譜クロマベクトルcktsとの内積を最新のフレームtから予め設定された時間Lだけ遡って累積した値であることを表す。以下の説明では、この時間Lを観測時間窓長;windowlengthと呼ぶ。なお、Lは、例えば2.5秒である。なお、式4によれば、クロマベクトル重み係数wchは、音響クロマベクトルctaと楽譜クロマベクトルcktsが類似するつまり、両ベクトルの方向が近くなるほど増加し、音響クロマベクトルctaと楽譜クロマベクトルcktsが一致するとき最大値1となる。クロマベクトル重み算出部10611は、算出したクロマベクトル重み係数wchを粒子重み算出部10614に出力する。なお、クロマベクトル重み算出部10611は、後述する粒子particlei毎に定められた仮想楽譜位置ki’及び仮想拍間隔bi’が、状態遷移部1063から入力される。クロマベクトル重み算出部10611は、観測時間の始点t-Lと終点tを、各々先行時間lookaheadtimeΔTだけ後の時刻に変更する。先行時間ΔTとは、状態遷移部1063において1回の状態遷移毎に観測時刻を進める時間であり、例えば1秒である。これにより、クロマベクトル重み算出部10611は、観測時間長Lを変更せずに観測時間をシフトさせる。また、クロマベクトル重み算出部10611は、上述の処理において粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’に基づきクロマベクトル重み係数wichを算出及び出力する。周波数特性重み算出部10612は、周波数特性分析部1021から入力された音響スペクトログラムXf,tと周波数特性合成部1051から入力された楽譜スペクトログラムX’f,kに基づいて周波数特性重み係数wspを算出する。周波数特性重み係数wspは、音響スペクトログラムXf,tと楽譜スペクトログラムX’f,kとの類似度又は相関性を示す尺度である。周波数特性重み算出部10612は、例えば、音響スペクトログラムXf,tと楽譜スペクトログラムX’f,kのカルバックライブラー情報量Kullback-LeiblerDivergenceに基づいて周波数特性重み係数wspと定める。ここで、周波数特性重み算出部10612は、次式6を用いて周波数特性重み係数wspを算出する。式6において、Qは、以下の式で表される実数である。但し、7式において、周波数特性重み算出部10612は、次式8に示すように規格化した音響スペクトログラムXf,tと楽譜スペクトログラムX’f,ktを用いる。また、周波数特性重み算出部10612においても、式5に示される時刻τと楽譜フレームkτとの関係を用いる。式6及び7によれば、周波数特性重み係数wspは、音響スペクトログラムXf,tと楽譜スペクトログラムX’f,ktが類似するほど増加し、両者が一致するとき、最大値1となる。このように、周波数特性重み算出部10612は、周波数特性重み係数wspを用いることで、クロマベクトルだけでは表現されない、入力音響信号や楽譜情報による周波数特性の広がりを考慮して確実に入力音響信号と楽譜情報との対応関係を把握することができる。周波数特性重み算出部10612は、算出した周波数特性重み係数wspを粒子重み算出部10614に出力する。なお、周波数特性重み算出部10612は、粒子i毎に定められた仮想楽譜位置ki’及び仮想拍間隔bi’が、状態遷移部1063から入力される。このとき、周波数特性重み算出部10612は、観測時間の始点t-Lと終点tを、各々先行時間ΔTだけ後の時刻に変更することで、観測時間を後にシフトさせる。また、周波数特性重み算出部10612は、上述の処理において粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’に基づき周波数特性重み係数wispを算出及び出力する。テンポ重み算出部10613は、周波数特性分析部1021から入力された音響スペクトログラムXf,tからテンポ重み係数wtを算出する。テンポ重み係数wtは、入力音響信号が示す仮想拍間隔b’の確からしさを示す尺度である。テンポ重み算出部10613は、例えば、音響スペクトログラムXf,tと仮想拍間隔b’だけずれた音響スペクトログラムXf,tの正規化相互相関normalizedcrosscorrelationをテンポ重み係数wtとして算出する。即ち、テンポ重み算出部10613は、次式9を用いてテンポ重み係数wtを算出する。式9において、仮想拍間隔b’が先行時間ΔTと等しいか、より大きいとき、テンポ重み算出部10613は、テンポ重み係数wtをゼロと定める。式9によれば、テンポ重み係数wtは、現在の音響スペクトログラムと過去の音響スペクトログラムの相関性を表す。例えば、同じパワーの大きさで同じ音階が拍間隔b’ごとに繰り返し演奏されている入力音響信号については、テンポ重み係数wtは、1となる。テンポ重み算出部10613は、算出したテンポ重み係数wtを粒子重み算出部10614に出力する。なお、テンポ重み算出部10613は、後述する粒子毎に定められた仮想拍間隔bi’が、状態遷移部1063から入力される。このとき、テンポ重み算出部10613は、観測時間の始点t-Lと終点tを、各々先行時間ΔTだけ後の時刻に変更することで、観測時間がシフトさせる。また、テンポ重み算出部10613は、上述の処理において粒子i毎の仮想拍間隔bi’に基づきテンポ重み係数witを算出及び出力する。粒子重み算出部10614は、クロマベクトル重み算出部10611から入力されたクロマベクトル重み係数wich、周波数特性重み算出部10612から入力された周波数特性重み係数wisp、及びテンポ重み算出部10613から入力されたテンポ重み係数witに基づき粒子i毎の重み係数wiを算出する。粒子iとは、本実施形態において仮想的な変数である仮想楽譜位置ki’、仮想拍間隔bi’及び観測値である重み係数wi、の組である。これらの粒子は、粒子フィルタリングparticlefiltering法により観測値から仮想的な変数を推定するために用いられる。例えば、本実施形態は、粒子i毎の仮想楽譜位置ki’、仮想拍間隔bi’に基づき、楽譜位置k’、拍間隔b’の推定値を定める。具体的には、粒子重み算出部10614は、粒子i毎の粒子重み係数wiを次式10に基づいて算出する。粒子重み算出部10614は、算出した粒子i毎の仮想楽譜位置ki’、仮想拍間隔bi’及び粒子重み係数wiを再標本化部1062に出力する。再標本化部1062は、粒子重み算出部10614から粒子i毎の仮想楽譜位置ki’、仮想拍間隔bi’及び粒子重み係数wiが入力される。再標本化部1062は、粒子i毎の粒子重み係数wiに基づき、各粒子を再標本化する。再標本化部1062は、再標本化した粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’を状態遷移部1063に出力する。再標本化部1062は、再標本化に係る処理において粒子数を一定に保つために、棄却する粒子数と、分割又は複製によって新たに生成する粒子数を等しくする。再標本化部1062は、例えばSIRsamplingimportanceresampling法を用いて粒子を再標本化する。ここで、再標本化部1062は、次の処理を行って粒子重みwiが小さい粒子iを棄却する。このような粒子iに対応する、仮想楽譜位置ki’や仮想拍間隔bi’の信頼性が低いからである。例えば、再標本化部1062は、粒子i毎の粒子重み係数を、次式11を用いて規格化し、規格化粒子重み係数piを算出する。式11において、Nは粒子数例えば、300個を示す。再標本化部1062は、規格化粒子重み係数piが予め設定された規格化粒子重み係数の閾値pthよりも小さいか否か判断する。再標本化部1062は、規格化粒子重み係数piが閾値pthよりも小さいと判断した場合、対応する粒子iを棄却する。即ち、再標本化部1062は、粒子iの規格化粒子重み係数pi、仮想楽譜位置ki’及び仮想拍間隔bi’を削除する。再標本化部1062は、規格化粒子重み係数piが閾値pthよりも大きい又は等しいと判断した場合、対応する粒子iを維持する。また、再標本化部1062は、棄却された粒子数をカウントする。再標本化部1062は、次の処理を行って粒子重みwiが大きい粒子iを分割又は複製する。このような粒子iに対応する、仮想楽譜位置ki’や仮想拍間隔bi’の信頼性が高いからである。例えば、再標本化部1062は、粒子i毎の粒子重み係数を、次式12を用いて再規格化し、再規格化粒子重み係数qiを算出する。再標本化部1062は、再規格化粒子重み係数qiが1よりも大きい粒子について整数部分の値nqi以下、再規格化整数と呼ぶと小数部分の値sqi以下、再規格化小数と呼ぶを算出する。再標本化部1062は、再規格化整数nqiが2又は2より多い粒子iを分割して、新たにnqi-1個の粒子jj≠iを各々生成する。即ち、再標本化部1062は、新たな粒子jに対応する仮想楽譜位置kj’及び仮想拍間隔bj’を分割直前の仮想楽譜位置ki‘及び仮想拍間隔bi’と定める。再標本化部1062は、棄却された粒子数から新たに生成した粒子数nqi-1を順次差し引くことにより、所定の粒子数Nに不足する粒子数Ns以下、不足粒子数と呼ぶを算出する。不足粒子数がゼロに達した場合には、再標本化部1062は、後述する粒子を複製する処理を行わない。他方、不足粒子数が1以上である場合に、再標本化部1062は、粒子を複製する処理を行う。そして、再標本化部1062は、新たに複製される粒子の総数が不足粒子数Nsとなるように、再規格化小数sqiに比例する確率以下、複製確率と呼ぶで、複製対象の粒子iを決定する。例えば、再標本化部1062は、次式13で算出した複製確率s’qiが、乱数ri但し、0<ri<1よりも小さい粒子を複製対象の粒子と定める。再標本化部1062は、定めた粒子各々について、新たに1個の粒子j≠iを複製する。即ち、再標本化部1062は、新たな粒子jに対応する仮想楽譜位置kj’及び仮想拍間隔bj’を複製直前の仮想楽譜位置ki’及び仮想拍間隔bi’と定める。再標本化部1062は、上述の処理を完了して決定した粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’を状態遷移部1063に出力する。状態遷移部1063は、再標本化部1062から粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’を入力する。状態遷移部1063は、仮想楽譜位置ki’及び仮想拍間隔bi’を、例えば式14及び式15を用いて更新する。式14及び15において、u、vは、各々平均値がゼロ、分散が各々σu、σvとなるガウシアン乱数gaussianrandomvariableである。例えば、σuは1.0であり、σvは0.2である。つまり、式14は、状態遷移部1063が、進行時間ΔT中の仮想拍間隔bi’テンポに応じた楽譜位置の進行に一定範囲に分布するランダムな揺らぎを示す乱数uを加えて仮想楽譜位置ki‘を更新していることを表す。他方、式15は、状態遷移部1063が、一定範囲に分布するランダムな揺らぎを示す乱数vを加えて仮想拍間隔bi’を更新していることを表す。このように、本実施形態は、粒子iごとに独立な状態遷移を行い、状態遷移において一定範囲内に分布する乱数を加えることで、演奏される音楽のテンポの揺らぎを考慮することができる。状態遷移部1063は、更新した粒子i毎の仮想楽譜位置ki’を重み算出部1061及び楽譜位置算出部1064に出力し、更新した粒子i毎の仮想拍間隔bi’を重み算出部1061及びテンポ算出部1065に出力する。楽譜位置算出部1064は、状態遷移部1063から入力された粒子i毎の仮想楽譜位置ki’に基づいて推定楽譜位置μを算出する。楽譜位置算出部1064は、例えば、仮想楽譜位置ki’の単純平均をμとする。本実施形態では、これに限られず、例えば楽譜位置算出部1064は、仮想楽譜位置ki’の偏差が所定の閾値uTHよりも小さくなる粒子間の平均値を推定楽譜位置μと算出してもよい。これにより、楽譜位置算出部1064は、平均値を算出する際、稠密に分布している粒子だけを考慮し、極端に分散した粒子を除外できる。楽譜位置算出部1064は、粒子i毎の仮想楽譜位置ki’と算出した推定楽譜位置μを信頼度算出部107に出力し、算出した推定楽譜位置μを楽譜位置出力部108に出力する。信頼度算出部107は、楽譜位置算出部1064から入力された粒子i毎の仮想楽譜位置ki’と算出した推定楽譜位置μに基づき、推定楽譜位置μの信頼度係数を算出する。信頼度算出部107は、例えば、次式16を用いて最新の時刻tにおける信頼度係数s2tを算出する。この信頼度係数s2tは、粒子i間の仮想楽譜位置ki’の分散であり、値が大きくなるほど推定楽譜位置の信頼性が低下することを示す。ここで、信頼度算出部107は、最新の時刻における信頼度係数s2tと予め定めた遷移回数I例えば5回前の時刻t-IΔTにおける信頼度係数s2t-IΔTの差分以下、信頼度差分と呼ぶが予め定めた閾値γincIを超えるか否かを判断する。また、信頼度算出部107は、信頼度差分が予め定めた閾値-γdecIを下回るか否かを判断する。γinc及びγdec任意の実数であり、例えばともに4である。信頼度算出部107は、信頼度差分が閾値γincIを超えると判断したとき、楽譜位置の出力を停止することを指示する楽譜位置出力停止信号を楽譜位置出力部108に出力する。信頼度算出部107は、信頼度差分が閾値-γdecIを下回ると判断したとき、楽譜位置の出力を開始することを指示する楽譜位置出力開始信号を楽譜位置出力部108に出力する。上述の例では、信頼度算出部107は、信頼度差分に基づき楽譜位置出力停止信号又は楽譜位置出力開始信号の出力を制御していたが、本実施形態ではこれに限らない。例えば、信頼度算出部107は、信頼度係数s2t自体が予め定めた閾値γsを超えるか否かを判断し、信頼度係数s2tが閾値γsを超えると判断したとき、楽譜位置出力停止信号を楽譜位置出力部108に出力し、信頼度係数s2tが閾値γsを下回ると判断したとき、楽譜位置出力開始信号を楽譜位置出力部108に出力してもよい。信頼度算出部107は、上記の信頼度係数の代わりに推定楽譜位置μの推定誤差etの絶対値|et|もしくは二乗誤差e2tを信頼度係数として算出し、算出した信頼度が予め定めた閾値γeを超えるか否かを判断し、信頼度係数が閾値γeを超えると判断したとき、楽譜位置出力停止信号を楽譜位置出力部108に出力し、信頼度係数が閾値γeを下回ると判断したとき、楽譜位置出力開始信号を楽譜位置出力部108に出力してもよい。ここで、推定誤差etは、次式17で表される。式17において、sk’tは、フレームk’tにおける真のgroundtruth時刻である。k’tは、推定楽譜位置μである。楽譜位置出力部108は、楽譜位置算出部1064から推定楽譜位置μが入力され、入力された推定楽譜位置μを出力する。但し、楽譜位置出力部108は、信頼度算出部107から楽譜位置出力停止信号が入力されたとき、推定楽譜位置μの出力を停止する。他方、楽譜位置出力部108は、信頼度算出部107から楽譜位置出力開始信号が入力されたとき、停止していた推定楽譜位置μの出力を再開する。これにより、本実施形態は、信頼性が低い推定楽譜位置μの出力を回避することができる。テンポ算出部1065は、状態遷移部1063から入力された粒子i毎の仮想拍間隔bi’に基づき推定拍間隔βを算出する。この推定拍間隔は、演奏されている音楽のテンポを示す尺度である。テンポ算出部1065は、例えば、仮想拍間隔bi’の単純平均をβとする。本実施形態では、これに限られず、例えばテンポ算出部1065は、仮想拍間隔bi’の偏差が所定の閾値νTHよりも小さくなる粒子間の平均値を推定拍間隔βと算出してもよい。これにより、テンポ算出部1065は、平均値を算出する際に稠密に分布している粒子だけを考慮し、極端に分散した粒子を除外できる。テンポ算出部1065は、推定拍間隔βをテンポ出力部109に出力する。テンポ出力部109は、テンポ算出部1065から推定拍間隔βが入力され、入力された推定拍間隔βを楽譜位置推定装置1の外部に出力する。ロボットは、例えば出力された推定楽譜位置μ及び推定拍間隔βを歌声の合成や機構部の位置制御に用い、合成した歌声や機構部の動作を演奏されている音楽に同期させることができる。本実施形態の楽譜位置推定処理は、例えばロボットが演奏されている音楽に反応して歌声を合成したり、動作機構を制御する、等に応用することが期待されている。これらの応用において必然的に遅延が生じる問題を、ロボットが本実施形態により推定された楽譜位置を用いることで解決できる。また、楽譜位置推定装置1は、当初、推定楽譜位置μと推定拍間隔βをともに出力するが、推定楽譜位置の信頼度増分が所定の閾値を越えると推定楽譜位置μの出力を停止する状態リズムレベルrhythmlevelに遷移する。推定楽譜位置の信頼度増分が所定の閾値を下回ると推定楽譜位置μの出力を再開する状態メロディレベルmelodylevelに遷移する。これにより、本実施形態では信頼性が低い推定楽譜位置が出力されなくなるため、歌声の合成や機構部の位置制御等において、利用者の意図しない動作をすることを回避できる。次に、本実施形態に係る楽譜位置推定処理の概要について図を用いて説明する。図2は、本実施形態に係る楽譜位置推定処理の概要を示す概念図である。図2左列は観測I.、中列は再標本化II.、右列は状態遷移III.、の過程を示す。このような過程を含む一連の処理は、粒子フィルタリングに含まれる処理である。図2左列は、上段から下段へ順に、音響信号、楽譜情報、粒子の分布を示す。各段ともに横軸は時刻である。上段は、縦軸は振幅であり音響信号の波形を表す。中段は、五線上に音符を用いて楽譜情報を表す。下段は、円を用いて仮想楽譜位置毎の粒子の分布を表す。円の半径は、各粒子の重み係数の大きさを表す。上段と中段に示されている四角の枠は、音響信号と楽譜情報を観測対象とする観測時間を示す。従って、図2左列は、観測I.の過程において、音響特徴量と楽譜特徴量に基づき、粒子毎に重み係数を算出することを示す。図2中列は、中段に再標本化前の粒子の分布を、下段に再標本化後の粒子の分布を示す。ここで、再標本化前の粒子の分布は、同図左列下段に示されているものと同一である。再標本化前の粒子から下段に向かう破線の矢印と×印は、その粒子が棄却されることを示す。再標本化前の粒子から再標本化後の粒子に向かう細い実線の矢印は、その粒子が維持されることを示す。再標本化前の粒子から再標本化後の粒子に向かう太い実線の矢印は、その粒子が分割又は複製され、同一の仮想楽譜情報を有する粒子が生成されることを示す。従って、図2中列は、再標本化において重み係数が小さい粒子が棄却され、重み係数が大きい粒子について分割又は複製されることにより棄却された数と同数の新たな粒子が生成されることを示す。図2右列は、上段から下段へ順に、楽譜情報、遷移前の粒子の分布、遷移後の粒子の分布を示す。ここで、楽譜情報は、同図左列中段に示されているものと同一である。遷移前の粒子の分布は、同図中列下段に示されているものと同一である。中段において下に向かう矢印は、粒子の分布の中心が現在の楽譜位置であることを示す。下段において上に向かう矢印は、粒子の分布の中心が推定楽譜位置であることを示す。ここで、図2右列は、現在の楽譜位置と推定楽譜位置の間が、拍間隔であることを示す。中段の粒子から下段に向かう矢印は、その矢印が示す楽譜位置に粒子が遷移することを示す。従って、図2右列は、状態遷移において粒子毎の仮想楽譜位置が更新され、更新後の仮想楽譜位置の分布に基づき楽譜位置が推定されることを示す。図3は、本実施形態に係る楽譜位置推定処理の概要を示すフローチャートである。ステップS1楽譜位置推定装置1の各構成部は、処理に用いる変数やデータを初期設定する。ステップS2音響信号入力部101は、音響信号を入力し、音響信号の入力が継続するか否か判断する。音響信号入力部101が、音響信号の入力が継続すると判断したときステップS2Y、ステップS3に進む。音響信号入力部101が、音響信号の入力が継続していないと判断したときステップS2N、処理を終了する。ステップS3音響特徴量生成部102は、入力された音響信号に基づき特徴量を生成する。重み算出部1061は、生成した音響信号に基づく特徴量と楽譜情報に基づく特徴量から粒子毎の重み係数を算出する。I.観測ステップS4再標本化部1062は、重み係数が所定の値より小さい粒子を棄却し、重み係数が所定の値より大きい粒子を分割又は複製する。II.再標本化ステップS5状態遷移部1063は、粒子毎の仮想楽譜位置及び仮想拍間隔を更新する。楽譜位置算出部1064は、更新後の仮想楽譜位置に基づき推定楽譜位置を算出する。テンポ算出部1065は、更新後の仮想拍間隔に基づき推定拍間隔を算出する。III.状態遷移次に、本実施形態に係る楽譜位置推定処理の詳細について図を用いて説明する。図4は、本実施形態に係る楽譜位置推定処理の一部を示すフローチャートである。図4は、初期設定ステップS1、音響信号入力の継続判断ステップS2及び観測ステップS3の各過程を含んで構成される。初期設定ステップS1は、ステップS101-S104を含んで構成される。ステップS101楽譜位置推定装置1の各構成部は、処理に用いる変数やデータを初期設定する。例えば、重み算出部1061は、粒子i毎の仮想拍間隔bi‘、粒子i毎の仮想楽譜位置ki‘、分散σu、σvを初期設定する。重み算出部1061は、粒子i毎の仮想拍間隔bi‘の初期値を、例えば各々テンポxに対応する拍間隔60/xに定める。このテンポxは、60bpmbeatperminute;1分あたりの拍数から180bpmの間で均一に分布したものである。重み算出部1061は、粒子i毎の仮想楽譜位置ki‘の初期値を、例えば、ΔT/bi‘と定める。重み算出部1061は、先行時間ΔTを、例えば1秒と設定しておく。その後、ステップS102に進む。ステップS102楽譜情報入力部104は、楽譜情報記憶部103から、処理の対象とする楽曲の楽譜情報nkを読み出す。楽譜情報入力部104は、読み出した楽譜情報nkを楽譜特徴量生成部105に出力する。その後、ステップS103に進む。ステップS103楽譜特徴量生成部105が備える周波数特性合成部1051は、楽譜情報入力部104から入力された楽譜情報nkに基づいて楽譜特徴量として、例えば楽譜スペクトログラムX’f,kを算出する。周波数特性合成部1051は、楽譜スペクトログラムX’f,kを算出する際、例えば式2に表されるharmonicGMMを用いる。周波数特性合成部1051は、その楽曲を構成する楽譜情報nkに全体について楽譜スペクトログラムX’f,kを算出しておいてもよい。また、予め算出した楽譜スペクトログラムX’f,kを楽譜情報記憶部103に記憶しておき、周波数特性合成部1051は、楽譜情報記憶部103から楽譜スペクトログラムX’f,kを読み出してもよい。周波数特性合成部1051は、生成または読み出した楽譜スペクトログラムX’f,kを楽譜位置推定部106に出力する。その後、ステップS104に進む。ステップS104楽譜特徴量生成部105が備えるクロマベクトル生成部1052は、楽譜情報入力部104から入力された楽譜情報nkに基づいて楽譜特徴量として、楽譜クロマベクトルcksを生成する。クロマベクトル生成部1052は、その楽曲を構成する楽譜情報nk全てについて楽譜クロマベクトルcksを生成しておいてもよい。また、予め生成した楽譜クロマベクトルcksを楽譜情報記憶部103に記憶しておき、クロマベクトル生成部1052は、楽譜情報記憶部103から楽譜クロマベクトルcksを読み出してもよい。クロマベクトル生成部1052は、生成又は読み出した楽譜クロマベクトルcksを楽譜位置推定部106に出力する。その後、ステップS2に進む。ステップS2音響信号入力部101は、音響信号を入力し、音響信号の入力が継続するか否か判断する。音響信号入力部101が、音響信号の入力が継続すると判断したときステップS2Y、ステップS301に進む。音響信号入力部101が、音響信号の入力が継続していないと判断したときステップS2N、処理を終了する。観測ステップS3の過程は、ステップS301-S307を含んで構成される。ステップS301音響信号入力部101は、アナログ音響信号を入力し、入力されたアナログ音響信号をディジタル信号に変換する。音響信号入力部101は、変換されたディジタル信号を音響特徴量生成部102に出力する。その後、ステップS302に進む。ステップS302周波数特性分析部1021は、音響信号入力部101から入力されたディジタル音響信号を時間領域から周波数領域信号に変換する。周波数特性分析部1021は、変換された周波数領域信号の振幅を音響スペクトログラムXf,tとして算出し、算出した音響スペクトログラムXf,tをクロマベクトル生成部1022及び楽譜位置推定部106に出力する。その後、ステップS303に進む。ステップS303クロマベクトル生成部1022は、周波数特性分析部1021から入力された音響スペクトログラムXf,tから音響クロマベクトルを算出する。ここで、クロマベクトル生成部1022は、例えば、音響クロマベクトルctaの音階n毎の要素ctanを、式1を用いて算出する。クロマベクトル生成部1022は、算出した音響クロマベクトルctaを楽譜位置推定部106に出力する。その後、ステップS304に進む。ステップS304クロマベクトル重み算出部10611は、クロマベクトル生成部1022から入力された音響クロマベクトルctaとクロマベクトル生成部1052から入力された楽譜クロマベクトルcktsからクロマベクトル重み係数wchを算出する。ここで、クロマベクトル重み算出部10611は、例えば式4及び5を用い、状態遷移部1063から入力された粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’に基づき、粒子i毎のクロマベクトル重み係数wichを算出する。クロマベクトル重み算出部10611は、算出した粒子i毎のクロマベクトル重み係数wichを粒子重み算出部10614に出力する。その後、ステップS305に進む。ステップS305周波数特性重み算出部10612は、周波数特性分析部1021から入力された音響スペクトログラムXf,tと周波数特性合成部1051から入力された楽譜スペクトログラムX’f,kから周波数特性重み係数wspを算出する。ここで、周波数特性重み算出部10612は、例えば式5-8を用い、状態遷移部1063から入力された粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’に基づき、粒子i毎の周波数特性重み係数wspを算出する。周波数特性重み算出部10612は、算出した粒子i毎の周波数特性重み係数wispを粒子重み算出部10614に出力する。その後、ステップS306に進む。ステップS306テンポ重み算出部10613は、周波数特性分析部1021から入力された音響スペクトログラムXf,tからテンポ重み係数wtを算出する。ここで、テンポ重み算出部10613は、例えば式9を用いて状態遷移部1063から入力された粒子i毎の仮想拍間隔bi’に基づき、粒子i毎のテンポ重み係数wtを算出する。テンポ重み算出部10613は、算出した粒子i毎のテンポ重み係数witを粒子重み算出部10614に出力する。その後、ステップS307に進む。ステップS307粒子重み算出部10614は、クロマベクトル重み算出部10611から入力されたクロマベクトル重み係数wich、周波数特性重み算出部10612から入力された周波数特性重み係数wisp、及びテンポ重み算出部10613から入力されたテンポ重み係数witから粒子i毎の重み係数wiを算出する。ここで、粒子重み算出部10614は、重み係数wiを算出する際、例えば式10を用いる。粒子重み算出部10614は、算出した粒子i毎の仮想楽譜位置ki’、仮想拍間隔bi’及び粒子重み係数wiを再標本化部1062に出力する。その後、ステップS411に進む。図5は、本実施形態に係る楽譜位置推定処理の他の一部を示すフローチャートである。図5は、再標本化ステップS4の過程を含んで構成される。再標本化ステップS4は、棄却ステップS41、分割ステップS42及び複製ステップS43の各過程を含んで構成される。棄却ステップS41の過程は、ステップS411-S414を含んで構成される。ステップS411再標本化部1062は、粒子重み算出部10614から粒子i毎の仮想楽譜位置ki’、仮想拍間隔bi’及び粒子重み係数wiを入力する。再標本化部1062は、全ての粒子iについてステップS412-S414の処理を実行したか否か判断する。再標本化部1062は、全ての粒子iについてステップS412-S414の処理を実行したと判断したときステップS411Y、ステップS421に進む。再標本化部1062は、全ての粒子iについてステップS412-S414の処理を実行していないと判断したときステップS411N、ステップS412に進む。ステップS412再標本化部1062は、例えば、式11を用いて、粒子i毎の粒子重み係数wiから規格化粒子重み係数piを算出する。再標本化部1062は、規格化粒子重み係数piが予め設定された閾値pthよりも小さいか否か判断する。再標本化部1062は、規格化粒子重み係数piが閾値pthよりも小さいと判断した場合ステップS411Y、ステップS413に進む。再標本化部1062は、規格化粒子重み係数piが閾値pthよりも大きい又は等しいと判断した場合ステップS411N、対応する粒子iを維持し、ステップS411に進む。ステップS413再標本化部1062は、規格化粒子重み係数piが閾値pthよりも小さいと判断された粒子iを棄却する。その後、ステップS414に進む。ステップS414再標本化部1062は、棄却された粒子数を1個増加させ、棄却された粒子数をカウントする。その後、ステップS411に進む。分割ステップS42の過程は、ステップS421-S426を含んで構成される。ステップS421再標本化部1062は、棄却されずに残った粒子iについてステップS422-S426の処理を実行したか否か判断する。再標本化部1062は、全ての残った粒子iについてステップS422-S426の処理を実行したと判断したときステップS421Y、ステップS431に進む。再標本化部1062は、全ての粒子iについてステップS422-S426の処理を実行していないと判断したときステップS421N、ステップS422に進む。ステップS422再標本化部1062は、例えば、式12を用いて、粒子i毎の規格化粒子重み係数piから再規格化粒子重み係数qiを算出する。再標本化部1062は、再規格化粒子重み係数qiが1よりも大きい粒子について再規格化整数nqiと再規格化小数sqiを算出する。その後、ステップS423に進む。ステップS423再標本化部1062は、再規格化整数nqiが2又は2より多いか否か判断する。再標本化部1062は、再規格化整数nqiが2又は2より多いと判断したときステップS423Y、ステップS424に進む。再標本化部1062が、再規格化整数nqiが1であると判断したときステップS423N、ステップS421に進む。ステップS424再標本化部1062は、粒子iを分割して、新たにnqi-1個の粒子jj≠iを生成する。即ち、再標本化部1062は、新たな粒子jに対応する仮想楽譜位置kj’及び仮想拍間隔bj’を分割前の仮想楽譜位置ki’及び仮想拍間隔bi’と定める。その後、ステップS425に進む。ステップS425再標本化部1062は、不足粒子数Ns初期値は棄却された粒子数から新たに生成した粒子数nqi-1を差し引いて不足粒子数をカウントする。その後、ステップS426に進む。ステップS426再標本化部1062は、不足粒子数Nsがゼロに達したか否か判断する。再標本化部1062は、不足粒子数Nsがゼロに達したと判断したときステップS426Y、再標本化に関する処理を終了し、粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’を状態遷移部1063に出力する。その後、ステップS501に進む。不足粒子数Nsがゼロに達していないと判断したときステップS426N、ステップS421に進む。複製ステップS43の過程は、ステップS431-S434を含んで構成される。ステップS431再標本化部1062は、棄却されずに残った粒子iについてステップS432-S433の処理を実行したか否か判断する。再標本化部1062は、全ての残った粒子iについてステップS432-S433の処理を実行したと判断したときステップS431Y、ステップS434に進む。再標本化部1062は、全ての粒子iについてステップS432-S433の処理を実行していないと判断したときステップS431N、ステップS432に進む。ステップS432再標本化部1062は、例えば式13を用いて再規格化小数sqiに比例する複製確率s’qiを算出する。その後、ステップS433に進む。ステップS433再標本化部1062は、複製確率s’qiに基づいて新たに複製される粒子の総数が不足粒子数Nsとなるように、複製対象の粒子iを決定する。その後、ステップS431に進む。ステップS434再標本化部1062は、定めた粒子各々について、新たに1個の粒子j≠iを複製する。即ち、再標本化部1062は、新たな粒子jに対応する仮想楽譜位置kj’及び仮想拍間隔bj’を複製前の仮想楽譜位置ki’及び仮想拍間隔bi’と定める。再標本化部1062は、粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’を状態遷移部1063に出力する。その後、ステップS501に進む。図6は、本実施形態に係る楽譜位置推定処理の残りの一部を示すフローチャートである。図6は、状態遷移ステップS5及び推定値出力ステップS6の過程を含んで構成される。状態遷移ステップS5の過程は、ステップS501-S504を含んで構成される。ステップS501状態遷移部1063は、再標本化部1062から粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’を入力する。状態遷移部1063は、進行時間ΔT中に仮想拍間隔bi’テンポで進行した楽譜位置に乱数uを加えて仮想楽譜位置ki’を更新する。状態遷移部1063は、仮想楽譜位置ki’を更新する際、例えば式14を用いる。その後、ステップS502に進む。ステップS502状態遷移部1063は、乱数vを加えて仮想拍間隔bi’を更新する。状態遷移部1063は、仮想拍間隔bi’を更新する際、例えば式15を用いる。状態遷移部1063は、更新した粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’を重み算出部1061に出力する。状態遷移部1063は、更新した粒子i毎の仮想楽譜位置ki‘を楽譜位置算出部1064に出力する。状態遷移部1063は、更新した粒子i毎の仮想拍間隔bi’をテンポ算出部1065に出力する。その後、ステップS503に進む。ステップS503楽譜位置算出部1064は、状態遷移部1063から入力された粒子i毎の仮想楽譜位置ki‘に基づいて推定楽譜位置μを算出する。楽譜位置算出部1064は、例えば、仮想楽譜位置ki‘の単純平均をμとする。楽譜位置算出部1064は、粒子i毎の仮想楽譜位置ki‘と算出した推定楽譜位置μを信頼度算出部107に出力する。楽譜位置算出部1064は、算出した推定楽譜位置μを楽譜位置出力部108に出力する。その後、ステップS503に進む。ステップS504テンポ算出部1065は、状態遷移部1063から入力された粒子i毎の仮想拍間隔bi’に基づいて推定拍間隔βを算出する。テンポ算出部1065は、例えば、仮想拍間隔bi’の単純平均をβとする。テンポ算出部1065は、推定拍間隔βをテンポ出力部109に出力する。その後、ステップS601に進む。推定値出力ステップS6の過程は、ステップS601-S606を含んで構成される。ステップS601信頼度算出部107は、楽譜位置算出部1064から入力された粒子i毎の仮想楽譜位置ki‘に基づき推定楽譜位置μを算出する。信頼度算出部107は、仮想楽譜位置ki‘と推定楽譜位置μに基づき、推定楽譜位置μの信頼度係数s2tを算出する。信頼度算出部107は、信頼度係数s2tを算出する際、例えば、式16を用いる。その後、ステップS602に進む。ステップS602信頼度算出部107は、信頼度差分が予め定めた閾値γincIを超えるか否かを判断する。信頼度算出部107は、信頼度差分が予め定めた閾値γincIを超えると判断したときステップS602Y、ステップS603に進む。信頼度算出部107は、信頼度差分が予め定めた閾値γincIを超えないと判断したときステップS602N、ステップS604に進む。ステップS603信頼度算出部107は、楽譜位置出力停止信号を楽譜位置出力部108に出力する。ここで、楽譜位置出力部108は、楽譜位置算出部1064から推定楽譜位置μが入力され、入力された推定楽譜位置μを出力する。楽譜位置出力部108は、信頼度算出部107から楽譜位置出力停止信号が入力されたとき、推定楽譜位置μの出力を停止する。その後、ステップS604に進む。ステップS604信頼度算出部107は、信頼度算出部107は、信頼度差分が予め定めた閾値-γdecIを下回るか否かを判断する。信頼度算出部107は、信頼度差分が予め定めた閾値-γdecIを下回ると判断したときステップS604Y、ステップS605に進む。信頼度算出部107は、信頼度差分が予め定め閾値-γdecIを下回らないと判断したときステップS604N、ステップS606に進む。ステップS605信頼度算出部107は、楽譜位置出力開始信号を楽譜位置出力部108に出力する。楽譜位置出力部108は、信頼度算出部107から楽譜位置出力開始信号が入力されたとき、停止していた推定楽譜位置μの出力を再開する。その後、ステップS606に進む。ステップS606テンポ出力部109は、テンポ算出部1065から推定拍間隔βが入力され、入力された推定拍間隔βを出力する。その後、ステップS7に進む。ステップS7重み算出部1061は、状態遷移部1063から粒子i毎の仮想楽譜位置ki’及び仮想拍間隔bi’が入力される。クロマベクトル重み算出部10611は、観測時間の始点t-Lと終点tを、各々先行時間ΔTだけ後の時刻にシフトさせて観測時間を変更する。重み算出部1061は、入力された仮想楽譜位置ki’、仮想拍間隔bi’及び観測時間を、クロマベクトル重み算出部10611及び周波数特性重み算出部10612に出力する。重み算出部1061は、入力された仮想拍間隔bi’及び観測時間を、テンポ重み算出部10613に出力する。重み算出部1061は、仮想楽譜位置ki’及び仮想拍間隔bi’を、粒子重み算出部10614に出力する。その後、ステップS2に進む。次に、本実施形態に係る楽譜位置推定装置1によって推定された楽譜位置の例について図を用いて説明する。図7-14は、特に断らない限り粒子数Nが300の場合における楽譜位置推定装置1の動作結果を示す。図7は、本実施形態に係る楽譜位置の推定誤差estimationerror、predictionerrorの一例を示す図である。図7において、縦軸は楽譜位置の推定誤差秒を示し、横軸は時刻秒を示す。図7によれば、時刻0-10秒又は20-30秒付近では、推定誤差は概ね±3秒未満である。時刻35-40秒付近での推定誤差が15秒前後と、より大きくなる。さらに、時刻41-50秒でも推定誤差は概ね±3秒未満となる。図8は、本実施形態に係る信頼度係数の一例を示す図である。図8において、縦軸は信頼度係数を示し、横軸は時刻秒を示す。図8によれば、時刻0-10秒又は2030秒付近では、信頼度係数が概ね20未満である。図7の結果を考慮すると、粒子が正しい楽譜位置に収束していることを示す。他方、信頼度係数が時刻35-40秒付近で約60から100へと増加する。この結果は、信頼度係数の増加が推定誤差の累積によって生じることを示す。また、図8によれば、時刻41-50秒付近で信頼度係数が低下する。この低下は、粒子が正しい楽譜位置に収束することを意味する。従って、本実施形態では、信頼度算出部107が、信頼度差分が所定の閾値以上であることを検知することで、顕著な推定誤差をもたらす推定楽譜位置の出力を停止することができる。また、信頼度算出部107が、信頼度差分が所定の閾値以下であることを検知することで、推定誤差が小さい推定楽譜位置の出力を開始することができる。これにより、推定誤差が著しい推定楽譜位置の出力を回避できることが裏付けられる。図9は、本実施形態に係る推定誤差が生じる割合の一例を示す図である。図9において、縦軸は、推定誤差の絶対値が1秒未満となる割合以下、低推定誤差率rateoflowpredictionerrorと呼ぶであり、横軸は、テンポbpmである。図9において、○及び×は、演奏された楽曲のテンポ毎の低推定誤差率である。○は、ドラム音を含む楽曲を示す。×は、ドラム音を含まない楽曲を示す。図9は、テンポが110bpm以下の楽曲のほとんどでは、低推定誤差率は0.9よりも高く、楽譜位置が精度よく推定されていることを示す。他方、図9は、テンポが120bpm以上の楽曲では、低推定誤差率は0.1よりも低く、楽譜位置が良好に推定されていないことを示す。その理由は、120bpm以上のテンポでは、拍間隔bが0.5秒以下となるため、先行時間ΔTこの例では1秒内にテンポ重み係数wt式9参照のピークが2個以上存在するからである。つまり、一意にテンポ重み係数wtが定まらないことが推定誤差を増加させる要因と考えられる。また、図9は、テンポが110bpm以下において、ドラム音を含む楽曲よりもドラム音を含まない楽曲のほうが、楽譜位置が精度よく推定されていることを示す。これは、ドラム音の非調和性inharmonicnatureがクロマベクトル重み係数wich、周波数特性重み算出部10612から入力された周波数特性重み係数wispを低下させる有力な原因となるからである。なお、上述の低推定誤差率とは異なり、0.5秒未満又は2秒未満の推定誤差をとる割合を低推定誤差率と定めた場合でも、本実施形態では図9と同様な結果が得られる。図10は、本実施形態に係る推定誤差の一例を示す図である。図10において、縦軸は、推定誤差であり、横軸は、テンポbpmである。図10において、○及び×は、演奏された楽曲のテンポ毎の平均推定誤差である。○は、ドラム音を含む楽曲を示す。×は、ドラム音を含まない楽曲を示す。また、○及び×の上下に伸びる棒線は、推定誤差の標準偏差を示す。図10は、テンポが110bpm以下の楽曲では、平均推定誤差が-10秒から0秒の範囲であり、テンポが120bpm以上の楽曲では、平均推定誤差が10秒から20秒に分布している。また、図10は、テンポが110bpm以下の楽曲のほとんどでは、標準偏差は7秒以内にとどまるのに対し、テンポが120bpm以上の楽曲では、標準偏差は10秒前後となる。即ち、図10は、テンポが110bpm以下の楽曲では、テンポが120bpm以上の楽曲よりも、楽譜位置が精度よく推定されることを示す。また、図10は、テンポが110bpm以下において、ドラム音を含まない楽曲よりもドラム音を含む楽曲のほうが、標準偏差が大きいことを示す。この点において、ドラム音を含む楽曲よりもドラム音を含まない楽曲のほうが、楽譜位置が精度よく推定されることを裏付ける。図11は、本実施形態に係る推定誤差が生じる割合の他の例を示す図である。この例は、楽譜位置出力部108が推定楽譜位置を出力する場合メロディレベルのみを考慮して算出した低推定誤差率を示す。図11において、縦軸は、低推定誤差率であり、横軸は、テンポbpmである。図11において、○及び×は、演奏された楽曲のテンポ毎の低推定誤差率である。○は、ドラム音を含む楽曲を示す。×は、ドラム音を含まない楽曲を示す。図11は、テンポが110bpm以下の楽曲のほとんどでは、低推定誤差率は0.9よりも高く、楽譜位置が精度よく推定されていることを示す。他方、図11は、テンポが120bpm以上の楽曲では、低推定誤差率は0.2から0.8の間を分布し、楽譜位置が良好に推定されていないことを示す。しかし、図11が示す低推定誤差率は、図9が示す低推定誤差率よりも高い。この結果は、信頼度算出部107が算出する信頼度に基づいて、楽譜位置出力部108からの推定楽譜位置の出力の有無を制御する構成が、楽譜位置の推定精度を向上させることを示す。図12は、本実施形態に係る推定誤差の他の例を示す図である。この例も、楽譜位置出力部108が推定楽譜位置を出力する場合メロディレベルのみを考慮して算出した低推定誤差率を示す。図12において、縦軸は、推定誤差であり、横軸は、テンポbpmである。図12において、○及び×は、演奏された楽曲のテンポ毎の平均推定誤差である。○は、ドラム音を含む楽曲を示す。×は、ドラム音を含まない楽曲を示す。また、○及び×の上下に伸びる棒線は、推定誤差の標準偏差を示す。図12は、テンポが110bpm以下の楽曲では、平均推定誤差が-2秒から2秒の範囲であり、テンポが120bpm以上のほとんどの楽曲では、平均推定誤差が10秒から15秒に分布している。また、図12は、テンポが110bpm以下の楽曲のほとんどでは、標準偏差は1秒以内にとどまるのに対し、テンポが120bpm以上の楽曲では、標準偏差は7~10秒となる。図10に示す結果と比較すると、図12に示す平均推定誤差、標準偏差ともに低くなる。この結果も、信頼度算出部107が算出する信頼度に基づいて、楽譜位置出力部108からの推定楽譜位置の出力の有無を制御する構成が、楽譜位置の推定精度を向上させることを示す。図13は、本実施形態に係る推定誤差と粒子数の関係の一例を示す図である。図13において、縦軸は推定誤差を、横軸は粒子数を示す。◇、□、及び△は、テンポが98、122、及び200bpmの楽曲について、粒子数毎の平均推定誤差を示す。◇、□、及び△の上下に伸びる棒線は、粒子数毎の標準偏差を示す。図13によれば、テンポが98bpmの楽曲について、平均推定誤差は、粒子数150について-5秒であるのに対し、粒子数300以上について-2~-3秒に収束する。標準偏差は、粒子数150について上下4秒であるのに対し、粒子数300以上について上下1~2秒に収束する。テンポが122bpmの楽曲について、平均推定誤差は、粒子数150について14秒であるのに対し、粒子数300以上について3~5秒に収束する。標準偏差は、粒子数150について±8秒であるのに対し、粒子数300以上について±3秒に収束する。これに対し、テンポが200bpmの楽曲について、粒子数に関わらず平均推定誤差は15秒前後である。標準偏差は粒子数に関わらず±8秒である。即ち、図13の結果は、テンポが約120bpm以下の楽曲について、粒子数が300あれば楽譜位置を精度よく推定できることを示す。反面、テンポが高い楽曲については楽譜位置を精度よく推定できないことを示す。図14は、本実施形態に係る精度precisionと再現率recallの一例を示す図である。図14において、縦軸は精度及び再現率を示し、横軸はテンポbpmを示す。ここで、精度ξpはxt/yと、再現率ξrはxt/zである。但し、xtは、推定誤差etの絶対値が1秒を越える場合であって、リズムレベル楽譜位置出力部108から推定楽譜位置を出力する場合における出力の個数である。yは、リズムレベルにおける全ての出力の個数である。zは、推定誤差etの絶対値が1秒を超える場合の出力の個数である。図14において、×は演奏された楽曲のテンポ毎の精度を、○は演奏された楽曲の再現率を示す。図14によれば、ほとんどの楽曲特にテンポが120bpm以上の全ての楽曲について精度が概ね1である。これは、信頼度算出部107において信頼度に基づいて推定楽譜位置の出力を停止する場合は、概ね推定誤差が著しいことを示す。図14によれば、再現率の平均は0.43である。ここで、テンポが110bpm以下の楽曲について、2曲を除いて再現率は0.05未満となり、テンポが120bpm以上の楽曲について、再現率は0.2~0.8の間で分布する。テンポが110bpm以下の楽曲について、そもそも推定誤差etの絶対値が1秒を超える場合が少ないため、全体として推定精度の低下させずに済む。上述の結果からも、本実施形態は、従来困難だった拍構造が不明確な楽曲例えば、ドラム音を含まない楽曲であっても演奏されている楽譜位置を精度よく推定できる。また、上述の結果において、テンポが速い楽曲について推定誤差が増加する傾向が認められた。その場合でも、先行時間ΔTをより短くすれば、テンポが速い楽曲でも楽譜位置を精度よく推定できることが示唆される。このように、本実施形態は、入力された音響信号の特徴量と楽譜情報の特徴量に基づいて重み係数を算出し、その重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定する。そのため、本実施形態によれば演奏された音楽に対する楽譜位置を精度よく推定することが可能になる。第2の実施形態以下、図面を参照しながら本発明の第2の実施形態について説明する。図15は、本実施形態に係る楽譜位置推定装置2の構成を示す概略図である。楽譜位置推定位置2は、楽譜位置推定装置1の楽譜位置推定部106及び信頼度算出部107の代わりに楽譜位置推定部206及び信頼度算出部207を備える。楽譜位置推定部206は、楽譜位置推定部106の重み算出部1061の代わりに重み算出部2061を備える。重み算出部2061は、クロマベクトル重み算出部10611、周波数特性重み算出部10612、及びテンポ重み算出部10613を備える点で重み算出部1061と同様である。また、楽譜位置推定装置2のその他の機能及び構成は第1の実施形態に係る楽譜位置推定装置1と同様である。以下、本実施形態について、第1の実施形態との差異点を主に説明する。信頼度算出部207は、例えば、上述の方式の何れかによって算出した信頼度係数を楽譜位置推定部206に出力する点が、信頼度算出部107と異なる。信頼度算出部207のその他の機能及び構成は信頼度算出部107と同様である。なお、信頼度算出部207は、楽譜位置推定部206に出力する信頼度係数として、楽譜位置出力停止信号及び楽譜位置出力開始信号の出力の要否を判断するために用いる信頼度係数とは異なる方式で算出してもよい。例えば、信頼度算出部207は、後述するDTWdynamictimewarping;動的時間収縮法を用いて距離を算出し、算出した距離を信頼度係数として楽譜位置推定部206に出力してもよい。後述するように、信頼度算出部207は、このDTW法を用いて音響特徴量と楽譜特徴量に基づき距離を算出する過程で、音響特徴量の区間と合致する楽譜特徴量の区間を定めることができる。本実施形態では、この距離を算出する構成を重み算出部2061に備え、この距離を算出してもよい。重み算出部2061は、信頼度算出部207から入力された信頼度係数に基づき観測時間を変更するか否か判断する。ここで、観測時間を変更する処理について説明する。図16は、本実施形態に係る観測時間を変更する処理の概要を示す概念図である。図16上段において、縦軸は振幅、横軸は時刻である、曲線は入力音響信号の波形の一例を表す。図16上段左側は、楽譜位置推定装置1の重み算出部1061と同様に観測時間をシフトする処理A.シフトを示す。実線の四角形は、現時刻の観測時間を示し、破線の四角形は次時刻の観測時間を示す。即ち、図16上段左側は、観測時間長Lを変更せずに先行時間ΔTだけ後の時刻に観測時間が変更することを示す。図16上段右側は、観測時間を延長させる処理B.延長を示す。実線の四角形は、現時刻の観測時間を示し、破線の四角形は次時刻の観測時間を示す。即ち、図16上段右側は、観測開始時刻を固定し観測終了時刻を先行時間ΔTだけ後の時刻に変更することにより、観測時間長がLからLΔTに延長することを示す。図16中段及び下段は、観測時間を現時刻までに入力された音響信号全体に延長する処理C.全区間に延長を示す。図16中段及び下段ともに、横軸は時刻を示す。図16中段において、縦軸は振幅を示し、曲線は入力音響信号の波形の一例を表す。図16下段は、楽譜情報を五線上に音階とその長さを示す音符を用いて示す。図16中段右側の実線の四角形は、現時刻の観測時間を示し、図16下段右側の実線の四角形は同図中段の実線の四角形と対応する区間の楽譜情報を示す。図16中段の破線の四角形は、次時刻の観測時間を示し、図16下段の破線の四角形は同図中段の破線の四角形と対応する区間の楽譜情報を示す。即ち、図16中段及び下段は、観測時間長LΔTの観測時間を現時刻までに入力された音響信号全体Lallに延長することを示す。次に本実施形態に係る重み算出部2061が観測時間を変更する処理について説明する。図17は、本実施形態に係る観測時間を変更する処理を示すフローチャートである。重み算出部2061は、以下に説明する処理を、図6のステップS7の代わりに行う。ステップS71重み算出部2061は、信頼度算出部207から信頼度係数が入力される。その後、ステップS72に進む。ステップS72重み算出部2061は、入力された信頼度係数が予め設定した閾値よりも大きい即ち、推定楽譜時刻の信頼性が低いか否か判断する。重み算出部2061は、信頼度が予め設定した閾値よりも小さいと判断したときステップS72N、ステップS73に進む。重み算出部2061は、信頼度係数が予め設定した閾値よりも大きい又は等しいと判断したときステップS72Y、ステップS74に進む。ステップS73重み算出部2061は、観測時間を、観測時間長Lを変更せずに先行時間ΔTだけシフトさせるA.シフト。その後、図4のステップS2に進む。ステップS74重み算出部2061は、誤り回数を1加算してカウントする。誤り回数とは、楽譜位置の推定が誤った回数であり、初期値はゼロである。その後、ステップS75に進む。ステップS75重み算出部2061は、誤り回数が予め設定した閾値例えば、5よりも大きいか否か判断する。重み算出部2061は、誤り回数が予め設定した閾値よりも大きいと判断したときステップS75Y、ステップS77に進む。重み算出部2061は、誤り回数が予め設定した閾値よりも小さい又は等しいと判断したときステップS75N、ステップS76に進む。ステップS76重み算出部2061は、観測時間の開始時刻を変更せず、終了時刻を先行時間ΔTだけ延長することで観測時間長を延長するB.延長。誤り回数が予め設定した閾値例えば、5よりも大きいか否か判断する。その後、図4のステップS2に進む。ステップS77重み算出部2061は、観測時間を現時刻までに入力された音響信号全体に延長するC.全区間に延長。その後、図4のステップS2に進む。これにより、楽譜位置の推定誤差が著しい場合でも、本実施形態が音響信号又は音響特徴量と楽譜情報又は楽譜特徴量の観測時間を延長して、両特徴量が合致する部分が見出される可能性が高くなる。これにより、本実施形態に係る楽譜位置推定装置2は、推定楽譜位置の精度が低下した状態から精度がより向上した状態に回復することができる。ひいては、本実施形態は、推定楽譜位置の信頼性が低く出力されない状態から、信頼性の低下を回避して推定楽譜位置の出力を再開することができる。次に、上述のDTW法を用いて信頼度係数を算出する処理について説明する。DTW法では、重み算出部2061又は信頼度算出部207以下、重み算出部2061等と呼ぶは、距離の算出対象とする音響特徴量の一部の区間と楽譜特徴量の一部の区間をマッチングさせ、全体の距離累積距離を最小化する音響特徴量の区間と楽譜特徴量の区間の組を決定する。即ち、重み算出部2061等は、この最小化された累積距離を信頼度係数と定める。以下の例では、このような音響特徴量の区間の始点を音響始点、終点を音響終点と呼び、このような楽譜特徴量区間の始点を楽譜始点、終点を楽譜終点と呼ぶ。図18は、本実施形態に係る楽譜位置推定装置2が行うDTW法を示すフローチャートである。ステップS711重み算出部2061等は、楽譜終点n2を、例えば楽譜特徴量の時系列データ冒頭のフレームを含むサブ系列に初期設定する。楽譜終点n2は、ここでは楽譜特徴量の時系列データにおいて累積距離Dn2を算出するための区間の終点、つまり、その区間の楽譜特徴量を含むサブ系列subsequenceである。サブ系列とは、楽譜特徴量の時系列データ及び音響特徴量の時系列データについて両者をマッチングすることにより両者が合致する区間を探索する単位となる小区間である。本実施形態では、観測時間が複数個のサブ系列からなり、1個のサブ系列に整数個の楽譜フレームが含まれるように、サブ系列の長さ例えば、2楽譜フレーム、1/24秒に相当を予め定めておく。その後、ステップS712に進む。ステップS712重み算出部2061等は、楽譜終点n2が楽譜特徴量の時系列データの最初のサブ系列に達するまで、ステップS713及びステップS714を実行する。その後、ステップS715に進む。ステップS713重み算出部2061等は、楽譜特徴量と音響特徴量をマッチング処理する。重み算出部2061等は、例えば後述するマッチング処理を行って、両特徴量間における累積距離Dn2を算出し、楽譜始点を示す対象楽譜サブ系列n1及び音響始点を示す対象音響サブ系列m1を決定する。決定された楽譜特徴量と音響特徴量とのマッチング処理の詳細については後述する。その後、ステップS714に進む。ステップS714重み算出部2061等は、楽譜終点n2を直前のサブ系列n2-1にシフトさせる。その後、ステップS712に進む。ステップS715重み算出部2061等は、累積距離Dn2が最小となる楽譜始点n1から楽譜終点n2までの楽譜特徴量を、音響特徴量と合致する区間と決定する。また、重み算出部1061等は、累積距離Dn2が最小となる音響始点m1から音響終点m2までの音響特徴量を、楽譜特徴量と合致する音響特徴量の区間と決定する。その後、図4のステップS2に進む。これにより、重み算出部2061は、楽譜特徴量と合致する音響特徴量の時系列の区間を観測時間として重み係数wiを算出することになる。次に、楽譜特徴量と音響特徴量とのマッチング処理について説明する。図19は、本実施形態に係る楽譜特徴量と音響特徴量とのマッチング処理を示すフローチャートである。ステップS721重み算出部2061等は、音響終点m2を、例えば音響特徴量の時系列データにおける最後のフレームを含むサブ系列に初期設定する。重み算出部2061等は、対象楽譜サブ系列n1を楽譜終点n2と初期設定し、対象音響サブ系列m1を音響終点m2と初期設定する。ここで、対象楽譜サブ系列n1とは、楽譜特徴量の時系列データのうち処理対象となるサブ系列を示す。対象音響サブ系列m1とは、音響特徴量の時系列データのうち処理対象となるサブ系列を示す。その後、ステップS722に進む。このように最後のサブ系列からマッチング処理を行うのは、楽譜位置推定装置2が最新の音響信号に基づいて楽譜位置を推定できるようにするためである。ステップS722重み算出部2061等は、累積距離Dn2の初期値を対象楽譜サブ系列n1における楽譜特徴量と対象音響サブ系列m1における音響特徴量との距離dn1,m1と初期設定する。距離dn1,m1は、楽譜特徴量と音響特徴量の類似度又は相関性を示す尺度である。距離dn1,m1は、両特徴量の類似度又は相関性が高いほど小さい値をとる関数であればよい。例えば、距離dn1,m1は、クロマベクトル重み係数wchに基づくエントロピー-wchlogwch、周波数特性重み係数wispに基づくエントロピー-wsplogwsp、又はこれらの合計値でもよい。ここで、重み算出部2061等は、例えば、クロマベクトル重み係数wchを算出する際、対象楽譜サブ系列n1又は対象音響サブ系列m1に相当する時間を観測時間として式4及び5を用いてもよい。また、重み算出部1061等は、例えば、周波数特性重み係数wispを算出する際、対象楽譜サブ系列n1又は対象音響サブ系列m1に相当する時間を観測時間として式5-8を用いてもよい。その後、ステップS723に進む。ステップS723重み算出部2061等は、対象楽譜サブ系列n1が楽譜特徴量の時系列データにおける最初のサブ系列に達する、又は対象音響サブ系列m1が音響特徴量の時系列データにおける最初のサブ系列に達するまで、ステップS724~ステップS726の処理を繰り返す。ステップS724重み算出部2061等は、次の3通り1-3のサブ系列の組それぞれについて距離を算出する。重み算出部2061等が、3通りの距離を算出するのは、演奏される音楽におけるテンポのゆらぎが、楽譜特徴量を相対的に遅延させる場合及び音響特徴量を相対的に遅延させる場合を考慮するためである。1直前の対象楽譜サブ系列n1-1における楽譜特徴量と現在の対象音響サブ系列m1における音響特徴量との間の距離dn1-1,m1、2直前の対象楽譜サブ系列n1-1における楽譜特徴量と直前の対象音響サブ系列m1-1における音響特徴量との間の距離dn1-1,m1-1、及び3現在の対象楽譜サブ系列n1における楽譜特徴量と直前の対象音響サブ系列m1-1における音響特徴量との間の距離dn1,m1-1。このステップにおいて、距離を算出する方法はステップS722と同様な方法でよい。重み算出部1061等は、dn1-1,m1、dn1-1,m1-1、及びdn1,m1-1のうち、最小の距離を選択する。その後、ステップS725に進む。ステップS725重み算出部2061等は、累積距離Dn2に選択された距離を加算する。その後、ステップS726に進む。ステップS726重み算出部2061等は、対象楽譜サブ系列n1及び対象音響サブ系列m1を選択された距離に対応するサブ系列の組に変更する。つまり、距離dn1-1,m1、が選択された場合、重み算出部2061等は、対象楽譜サブ系列n1を直前の対象楽譜サブ系列n1-1に変更し、対象音響サブ系列m1を変更しない。距離dn1-1,m1-1、が選択された場合、重み算出部2061等は、対象楽譜サブ系列n1を直前の対象楽譜サブ系列n1-1に変更し、対象音響サブ系列m1を直前の対象音響サブ系列m1-1に変更する。距離dn1-1,m1-1、が選択された場合、重み算出部2061等は、対象楽譜サブ系列n1を変更せず、対象音響サブ系列m1を直前の対象音響サブ系列m1-1に変更する。その後、ステップS723に進む。これにより、重み算出部2061等は、音響特徴量と楽譜特徴量が合致する区間を特定して、重み算出部2061は、特定された区間に基づき重み係数wiを算出する。よって、本実施形態に係る楽譜位置推定装置2は、さらに楽譜位置の推定誤りから容易に回復し、推定楽譜位置の出力を再開することができる。なお、上述の観測時間を変更する処理において、観測時間を全区間に延長C.した後にDTW法を用いて楽譜特徴量と音響特徴量が合致する区間を特定することを前提にした説明したが、本実施形態ではこれには限られない。本実施形態では、シフトA.した後や、先行時間ΔTだけ延長B.した後の観測時間における楽譜特徴量と音響特徴量の時系列データにDTW法を用いて合致する区間を特定してもよい。これにより、逐次に楽譜位置の推定誤差が精度よく算出できる状態を回復することができる。本実施形態では、観測時間を延長B.する際の延長時間は先行時間ΔTに限られず、任意の時間でよいし、始点を固定した時刻にせず先行させてもよい。なお、本実施形態では、より大きいサブ系列マッチング処理を行う区間数がより少ないに基づきDTW法を実行して楽譜特徴量を基準に合致する音響特徴量の区間の候補を決定しておき、その決定された区間の候補に対し、より小さいサブ系列マッチング処理を行う区間数がより多いに基づきDTW法を実行してもよい。また、上述のDTW法においては、楽譜特徴量を基準に合致する音響特徴量の区間を決定していたが、本実施形態では、音響特徴量を基準に合致する楽譜特徴量の区間を決定してもよい。本実施形態において、周波数特性分析部1021は、入力されたディジタル音響信号が周波数領域信号に変換できれば、FFT以外の方式、例えばDFTDiscreteFourierTransform、離散フーリエ変換を用いてもよい。本実施形態において、再標本化部1062は、最も規格化粒子重み係数piが小さいものから予め定めた数Nr例えば50個の粒子を棄却してもよい。なお、上述した実施形態における楽譜位置推定装置1又は2の一部、例えば、音響特徴量生成部102、楽譜特徴量生成部105、楽譜位置推定部106、206、及び信頼度算出部107、207をコンピュータで実現するようにしても良い。その場合、この制御機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現しても良い。なお、ここでいうコンピュータシステムとは、楽譜位置推定装置1又は2に内蔵されたコンピュータシステムであって、OSや周辺機器等のハードウェアを含むものとする。また、コンピュータ読み取り可能な記録媒体とは、フレキシブルディスク、光磁気ディスク、ROM、CD-ROM等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらにコンピュータ読み取り可能な記録媒体とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含んでも良い。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。また、上述した実施形態における楽譜位置推定装置1又は2の一部、または全部を、LSILargeScaleIntegration等の集積回路として実現しても良い。楽譜位置推定装置1又は2の各機能ブロックは個別にプロセッサ化してもよいし、一部、または全部を集積してプロセッサ化しても良い。また、集積回路化の手法はLSIに限らず専用回路、または汎用プロセッサで実現しても良い。また、半導体技術の進歩によりLSIに代替する集積回路化の技術が出現した場合、当該技術による集積回路を用いても良い。以上、図面を参照してこの発明の一実施形態について詳しく説明してきたが、具体的な構成は上述のものに限られることはなく、この発明の要旨を逸脱しない範囲内において様々な設計変更等をすることが可能である。"