(57)【要約】【課題】演奏された音楽に対する楽譜位置の推定を行う楽譜位置推定装置、楽譜位置推定方法、及び楽譜位置推定プログラムを提供する。【解決手段】音響特徴量生成部は入力された音響信号の特徴量を所定の区間ごとに生成し、楽譜位置推定部は前記音響信号の特徴量と楽譜情報の特徴量に基づき重み係数を算出し、前記重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定する。【選択図】図１
2011180590.tif 000002
【特許請求の範囲】
【請求項１】
入力された音響信号の特徴量を生成する音響特徴量生成部と、
前記音響信号の特徴量と楽譜情報の特徴量に基づき重み係数を算出し、前記重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定する楽譜位置推定部と、
を備えることを特徴とする楽譜位置推定装置。
【請求項２】
前記仮想テンポと第１の乱数に基づいて前記仮想楽譜位置を更新し、第２の乱数に基づいて前記仮想テンポを更新する状態遷移部を備えることを特徴とする請求項１に記載の楽譜位置推定装置。
【請求項３】
前記特徴量は、音階ごとの強度を複数個含むクロマベクトル又は周波数特性を含み、
前記楽譜位置推定部は、
前記重み係数を、前記音響信号によるクロマベクトルと前記楽譜情報によるクロマベクトルとの類似度を示すクロマベクトル重み係数、前記前記音響信号の周波数特性と前記楽譜情報の周波数特性との類似度を示す周波数特性重み係数、又は前記音響信号の周波数特性の自己相関に基づくテンポ重み係数を用いて算出すること、
を特徴とする請求項１又は２に記載の楽譜位置推定装置。
【請求項４】
前記仮想楽譜位置に基づいて前記推定した楽譜位置の信頼度を算出し、前記信頼度に基づいて前記推定した楽譜位置を出力するか否か判断する信頼度算出部を備えること、
を特徴とする請求項１乃至３のうちいずれか１項に記載の楽譜位置推定装置。
【請求項５】
前記推定した楽譜位置の信頼度を算出する信頼度算出部を備え、
前記楽譜位置推定部は、
前記算出された信頼度に基づき前記重み係数を算出する前記特徴量の区間を延長するか否か判断すること、
を特徴とする請求項１乃至３のうちいずれか１項に記載の楽譜位置推定装置。
【請求項６】
楽譜位置推定装置における楽譜位置推定方法において、
前記楽譜位置推定装置が、入力された音響信号の特徴量を所定の区間ごとに生成する第１の過程と、
前記楽譜位置推定装置が、前記音響信号の特徴量と楽譜情報の特徴量に基づき重み係数を算出し、前記重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定する第２の過程と
を有すること、
を特徴とする楽譜位置推定方法。
【請求項７】
楽譜位置推定装置のコンピュータに、
入力された音響信号の特徴量を生成する手順と、
前記音響信号の特徴量と楽譜情報の特徴量に基づき重み係数を算出し、前記重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定する手順と、
を実行させるための楽譜位置推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、楽譜位置推定装置、楽譜位置推定方法、及び楽譜位置推定プログラムに関する。
【背景技術】
【０００２】
近年、ロボットが人間社会においてますます利用されることが期待されている。このようなロボットは、人間とのインタラクション（相互作用；ｉｎｔｅｒａｃｔｉｏｎ）を行う手段を備えることが求められている。その手段を実現するためのメディア（媒体；ｍｅｄｉｕｍ）は、例えば言語による対話、音楽、アイコンタクト（視線を合わせること；ｅｙｅ  ｃｏｎｔａｃｔ）等である。とりわけ、音楽は、有望なメディアであり、人間社会において地域や人種を越えた重要かつ共通な要素である。例えば、ロボットが、音楽にあわせて人間と踊ったり、歌ったり、又は楽器を演奏することが期待されている。
【０００３】
例えば、特許文献１には楽譜を２次元画像データに変換し、この画像データから五線データ及び小節線データを除去し、音符および記号を認識する画像認識装置が記載されている。また、特許文献２には、楽譜データおよびグルーピングされた構造分析結果に基づき、楽曲中の各拍点の局所的な強さの値を算出し、各拍点に対応づけてその強さの値を記憶手段に記憶し、次の階層の拍節構造のすべての候補において、その候補に含まれる拍点の強さの値を記憶手段から読み出して総和を求め、それぞれの候補における拍点の強さの値の総和を比較し、総和が最も大きい候補を次の階層の拍節構造として選択する楽曲の拍節構造の自動分析方法が提案されている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特許第３１４７８４６号公報
【特許文献２】特開２００６−２０１２７８号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかし、特許文献２記載の方法では、収音した楽曲が楽譜上のどの位置（以下、楽譜位置（Ｓｃｏｒｅ  Ｐｏｓｉｔｉｏｎ）と呼ぶ）にあるかを認識できない。つまり、特許文献２記載の方法は、演奏されている曲のビート（拍；ｂｅａｔ）時間やテンポの抽出に失敗することがあるという問題点があった。
【０００６】
本発明は上記の点に鑑みてなされたものであり、演奏された音楽に対する楽譜位置を推定できる楽譜位置推定装置、楽譜位置推定方法、及び楽譜位置推定プログラムを提供することを課題としている。
【課題を解決するための手段】
【０００７】
（１）本発明は上記の課題を解決するためになされたものであり、本発明の一態様は、入力された音響信号の特徴量を生成する音響特徴量生成部と、前記音響信号の特徴量と楽譜情報の特徴量に基づき重み係数を算出し、前記重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定する楽譜位置推定部とを備えることを特徴とする楽譜位置推定装置である。
【０００８】
（２）本発明の他の態様は、前記仮想テンポと第１の乱数に基づいて前記仮想楽譜位置を前記所定の区間ごとに更新し、第２の乱数に基づいて前記仮想テンポを前記所定の区間ごとに更新する状態遷移部を備えることを特徴とする（１）の楽譜位置推定装置である。
【０００９】
（３）本発明の他の態様は、前記特徴量は、音階ごとの強度を複数個含むクロマベクトル又は周波数特性を含み、前記楽譜位置推定部は、前記重み係数を、所定の区間における前記音響信号によるクロマベクトルと前記楽譜情報によるクロマベクトルとの類似度を示すクロマベクトル重み係数、前記所定の区間における前記音響信号の周波数特性と前記楽譜情報の周波数特性との類似度を示す周波数特性重み係数、又は前記所定の区間における前記音響信号の周波数特性の自己相関に基づくテンポ重み係数を用いて算出することを特徴とする（１）又は（２）の楽譜位置推定装置である。
【００１０】
（４）本発明の他の態様は、前記仮想楽譜位置に基づいて前記推定した楽譜位置の信頼度を算出し、前記信頼度に基づいて前記推定した楽譜位置を出力するか否か判断する信頼度算出部を備えることを特徴とする（１）、（２）、又は（３）のいずれかの楽譜位置推定装置である。
【００１１】
（５）本発明の他の態様は、前記推定した楽譜位置の信頼度を算出する信頼度算出部を備え、前記楽譜位置推定部は、前記算出された信頼度に基づき前記重み係数を算出する前記特徴量の区間を延長するか否か判断することを特徴とする（１）、（２）、又は（３）のいずれかの楽譜位置推定装置である。
【発明の効果】
【００１２】
（１）に記載した態様によれば、入力された音響信号の特徴量を生成し、音響信号の特徴量と楽譜情報の特徴量に基づき重み係数を算出し、重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定するので、演奏された音楽に対する楽譜位置を精度よく推定することが可能になる。
（２）に記載した態様によれば、仮想テンポと第１の乱数に基づいて仮想楽譜位置を更新し、第２の乱数に基づいて仮想テンポを更新するため、演奏された音楽のテンポが変動しても楽譜位置を頑健に推定することが可能になる。
（３）に記載した態様によれば、クロマベクトル又は周波数特性を含む特徴量に基づき音響信号と楽譜情報との類似度を示す重み係数を算出するため、音響信号と楽譜情報が類似する仮想楽譜位置及び仮想テンポが重視される。そのため、当該発明は演奏された音楽に対する楽譜位置を精度よく推定することが可能になる。
（４）に記載した態様によれば、推定した楽譜位置の信頼度に基づいて推定した楽譜位置を出力するか否か判断するため、出力された楽譜位置の信頼度を確保することが可能になる。
（５）に記載した態様によれば、算出された信頼度に基づき重み係数を算出する特徴量の区間を延長するか否か判断するため、推定した楽譜位置の信頼度が確保される特徴量を含むように区間を延長して、出力された楽譜位置の信頼度を確保することが可能になる。
【図面の簡単な説明】
【００１３】
【図１】本発明の第１の実施形態に係る楽譜位置推定装置の構成を示す概略図である。
【図２】本実施形態に係る楽譜位置推定処理の概要を示す概念図である。
【図３】本実施形態に係る楽譜位置推定処理の概要を示すフローチャートである。
【図４】本実施形態に係る楽譜位置推定処理の一部を示すフローチャートである。
【図５】本実施形態に係る楽譜位置推定処理の他の一部を示すフローチャートである。
【図６】本実施形態に係る楽譜位置推定処理の残りの一部を示すフローチャートである。
【図７】本実施形態に係る楽譜位置の推定誤差の一例を示す図である。
【図８】本実施形態に係る信頼度の一例を示す図である。
【図９】本実施形態に係る推定誤差が生じる割合の一例を示す図である。
【図１０】本実施形態に係る推定誤差の一例を示す図である。
【図１１】本実施形態に係る推定誤差が生じる割合の他の例を示す図である。
【図１２】本実施形態に係る推定誤差の他の例を示す図である。
【図１３】本実施形態に係る推定誤差と粒子数の関係の一例を示す図である。
【図１４】本実施形態に係る精度と再現率の一例を示す図である。
【図１５】本発明の第２の実施形態に係る楽譜位置推定装置の構成を示す概略図である。
【図１６】本実施形態に係る観測時間を変更する処理の概要を示す概念図である。
【図１７】本実施形態に係る観測時間を変更する処理を示すフローチャートである。
【図１８】本実施形態に係る楽譜位置推定装置が行うＤＴＷ法を示すフローチャートである。
【図１９】本実施形態に係る楽譜特徴量と音響特徴量とのマッチング処理を示すフローチャートである。
【発明を実施するための形態】
【００１４】
（第１の実施形態）
以下、図面を参照しながら本発明の第１の実施形態について説明する。
図１は、本実施形態に係る楽譜位置推定装置１の構成を示す概略図である。楽譜位置（ｓｃｏｒｅ  ｐｏｓｉｔｉｏｎ）推定装置１は、音響信号入力部１０１、音響特徴量生成部１０２、楽譜情報記憶部１０３、楽譜情報入力部１０４、楽譜特徴量生成部１０５、楽譜位置推定部１０６、信頼度（ｃｏｎｆｉｄｅｎｃｅ）算出部１０７、楽譜位置出力部１０８、及びテンポ出力部１０９を含んで構成される。
【００１５】
音響信号入力部１０１は、演奏された音楽による音波をディジタル音響信号に変換し、変換したディジタル音響信号を音響特徴量生成部１０２に出力する。音響信号入力部１０１は、例えばマイクロホンとアナログ・ディジタル（Ａｎａｌｏｇ−ｔｏ−Ｄｉｇｉｔａｌ；Ａ／Ｄ）変換部（図示せず）とを含んで構成される。このマイクロホンは、人間が聴取することができる周波数帯域（例えば、２０Ｈｚ−２０ｋＨｚ）の音波を受信してアナログ音響信号を生成し、生成したアナログ音響信号をＡ／Ｄ変換部に出力する。Ａ／Ｄ変換部は、マイクロホンから入力されたアナログ音響信号をディジタル音響信号に変換し、変換したディジタル音響信号を音響特徴量生成部１０２に出力する。ここで、Ａ／Ｄ変換部は、入力されたアナログ音響信号を、例えば、サンプリング周波数４４．１ｋＨｚ、振幅を１６ビットの２進（バイナリ；ｂｉｎａｒｙ）データにパルス符号変調（Ｐｕｌｓｅ  Ｃｏｄｅ  Ｍｏｄｕｌａｔｉｏｎ；ＰＣＭ）することで、量子化されたディジタル音響信号に変換する。
【００１６】
音響特徴量生成部１０２は、音響信号入力部１０１から入力されたディジタル音響信号から、その特徴量（音響特徴量）を生成し、生成した音響特徴量を楽譜位置推定部１０６に出力する。
具体的には、音響特徴量生成部１０２は、音響特徴量として、例えば、スペクトログラム（ｓｐｅｃｔｒｏｇｒａｍ）及びクロマベクトル（ｃｈｒｏｍａ  ｖｅｃｔｏｒ）を生成する。音響特徴量生成部１０２は、周波数特性分析部１０２１及びクロマベクトル生成部１０２２を含んで構成される。
周波数特性分析部１０２１は、音響信号入力部１０１から入力されたディジタル音響信号を時間領域信号から周波数領域信号に変換する。ここで、周波数特性分析部１０２１は、例えば、２０４８個のサンプルからなるフレームごとにディジタル音響信号を高速フーリエ変換（Ｆａｓｔ  Ｆｏｕｒｉｅｒ  Ｔｒａｎｓｆｏｒｍ；ＦＦＴ）して周波数領域信号に変換する。このフレームは、ＦＦＴを行うディジタル音響信号の範囲である。周波数特性分析部１０２１は、予め設定された時刻ｔ毎（例えば、１０ｍｓ間隔）に新たな信号のサンプル（例えば、サンプリング周波数が４４．１ｋＨｚの場合、４４１サンプル）が含まれるようにフレームをシフトさせる。
周波数特性分析部１０２１は、フレーム時刻ｔ毎に変換した周波数領域信号の振幅（絶対値）をスペクトログラム（以下、音響スペクトログラムと呼ぶ）Ｘｆ，ｔとして算出し、算出した音響スペクトログラムＸｆ，ｔをクロマベクトル生成部１０２２及び楽譜位置推定部１０６に出力する。なお、ｆは周波数を表す。
【００１７】
クロマベクトル生成部１０２２は、周波数特性分析部１０２１から入力された音響スペクトログラムＸｆ，ｔに基づいて音響信号によるクロマベクトル（以下、音響クロマベクトルと呼ぶ）ｃｔａを算出する。クロマベクトルとは、音階（ｃｈｒｏｍａ）毎の成分のパワーを要素として含むベクトルである。クロマベクトルは、例えば、西洋音楽を構成する１２音階（Ｃ、Ｃ＃、Ｄ、Ｄ＃、Ｅ、Ｆ、Ｆ＃、Ｇ、Ｇ＃、Ａ、Ａ＃、Ｂ）における音階毎の成分の強度を要素とするベクトル（例えば、要素数が１２）である。
クロマベクトル生成部１０２２は、音響クロマベクトルｃｔａ＝［ｃｔａ（１），ｃｔａ（２），…，ｃｔａ（１２）］の音階ｎ毎の要素ｃｔａ（ｎ）を、例えば、次式を用いて算出する。
【００１８】
【数１】
2011180590.tif 000003
【００１９】
式（１）において、ＯｃｔＬは予め設定されたオクターブの下限（下限オクターブ）を表す。ＯｃｔＨは、予め設定されたオクターブの上限（上限オクターブ）を表す。ＢＰＦｎ，ｈ（ｆ）は、ｈ番目のオクターブにおける音階ｎの成分を抽出する帯域通過フィルタ（バンドパスフィルタ、Ｂａｎｄ  Ｐａｓｓ  Ｆｉｌｔｅｒ；ＢＰＦ）の入出力特性を表す。ｈ番目のオクターブにおける音階ｎの周波数ｆｎ，ｈにおいて、その値は最大となり、周波数ｆがゼロ又は無限大に近づくと、ゼロに近似する関数である。また、｜Ｘｆ，ｔ｜２は、パワースペクトルを示す。楽譜位置推定装置１は、このクロマベクトルを用いることで、オクターブが異なっても音階毎の強度を抽出でき、音響信号の調波構造（ｈａｒｍｏｎｉｃ  ｓｔｒｕｃｔｕｒｅ）を無視することができる。よって、楽譜位置推定装置１は、旋律を表すために重要な音階の計時変化を抽出でき、旋律を表すために余り重要ではないオクターブの差異を除外することができる。
クロマベクトル生成部１０２２は、算出した音響クロマベクトルｃｔａを楽譜位置推定部１０６に出力する。
【００２０】
楽譜情報記憶部１０３は、楽譜情報を楽曲毎に記憶する。楽譜情報は、楽譜フレーム（ｓｃｏｒｅ  ｆｒａｍｅ）ｋ（ｋは１又は１より大きい整数）毎に、楽曲を構成する音階を示す音階情報ｎｋｒを要素とするベクトルｎｋ＝［ｎｋ１，ｎｋ２，…，ｎｋｒｋ］Ｔで表される。ここで、ｒは１又は１より大きくｒｋ又はｒｋよりも小さい整数である。ｒｋは、楽譜フレームｋにおける音階の数を表す。Ｔは、ベクトルの転置を表す。楽譜フレームとは、１つの楽曲の楽譜情報において楽曲を時間分割する単位時間である。例えば、四分音符の長さが１２フレームである場合、楽譜情報の時間分解能（ｒｅｓｏｌｕｔｉｏｎ）は十六分音符の三分の一となる。例えば、楽譜フレーム２において、音階Ｃ、Ｅ及びＧを表す楽譜情報ｎ２は、要素数が３個のベクトル［Ｃ，Ｅ，Ｇ］となる。
【００２１】
楽譜情報入力部１０４は、楽譜情報記憶部１０３から、処理の対象とする楽曲の楽譜情報ｎｋを読み出す。楽譜情報入力部１０４は、読み出した楽譜情報ｎｋを楽譜特徴量生成部１０５に出力する。
【００２２】
楽譜特徴量生成部１０５は、楽譜情報入力部１０４から入力された楽譜情報ｎｋに基づいて、その特徴量（以下、楽譜特徴量と呼ぶ）を生成し、生成した楽譜特徴量を楽譜位置推定部１０６に出力する。楽譜特徴量生成部１０５は、楽譜特徴量として、例えば、スペクトログラム及びクロマベクトルを生成する。ここで、楽譜特徴量生成部１０５が生成したスペクトログラム及びクロマベクトルを、それぞれ楽譜スペクトログラム及び楽譜クロマベクトルと呼ぶ。
楽譜特徴量生成部１０５は、周波数特性合成部１０５１及びクロマベクトル生成部１０５２を含んで構成される。
周波数特性合成部１０５１は、楽譜情報入力部から入力された楽譜情報ｎｋに基づき、楽譜スペクトログラムＸ’ｆ，ｋを生成する。楽譜スペクトログラムＸ’ｆ，ｋとは、楽譜情報ｎｋに基づき演奏されるべき音楽を示す音響信号の周波数特性である。
周波数特性合成部１０５１は、例えば次式のように調和混合正規分布モデル（ｈａｒｍｏｎｉｃ  ｇａｕｓｓｉａｎ  ｍｉｘｔｕｒｅ  ｍｏｄｅｌ；ｈａｒｍｏｎｉｃ  ＧＭＭ）に基づいて楽譜スペクトログラムＸ’ｆ，ｋを生成する。
【００２３】
【数２】
2011180590.tif 000004
【００２４】
式（２）の右辺第１項において、ｇは、調波インデックス（ｈａｒｍｏｎｉｃ  ｉｎｄｅｘ）即ち、基本周波数（ｆｕｎｄａｍｅｎｔａｌ  ｆｒｅｑｕｅｎｃｙ）Ｆの整数倍（ｉｎｔｅｇｅｒ  ｍｕｌｔｉｐｌｅ）の周波数の成分を示すインデックスである。ｇ＝１は、基本波（ｆｕｎｄａｍｅｎｔａｌ  ｔｏｎｅ）を示す。Ｇは、調波の数（例えば１０）を示す。ｈ（ｇ）は、第ｇ高調波（ｇ−ｔｈ  ｈａｒｍｏｎｉｃ）の振幅（例えば、０．２のｇ乗）を表す。Ｎ（ｆ；ｇＦｎｋ，σ２）は、パワーの周波数特性が、中心周波数ｇＦｎｋ、分散σ２の正規分布（ｎｏｒｍａｌ  ｄｉｓｔｒｉｂｕｔｉｏｎ，ｇａｕｓｓｉａｎ  ｄｉｓｔｒｉｂｕｔｉｏｎ）関数を表す。即ち、式（２）の右辺第１項は純粋な調和混合正規分布モデルに基づき、各音階を構成する調波成分（ｈａｒｍｏｎｉｃ）を線形結合した成分を示す。なお、分散σ２は、例えば０．８である。
【００２５】
式（２）の右辺第２項において、Ｃ（ｆ）は、ピンクノイズ（ｐｉｎｋ  ｎｏｉｓｅ）成分のパワー（ｐｏｗｅｒ）の周波数特性を示す。ピンクノイズとは、オクターブ帯域毎のパワーが等しい周波数特性を有しランダムな波形を表す信号系列である。式（２）の右辺において、Ｃ（ｆ）を加算するのは、右辺第１項の成分が表すスペクトログラムを処理すると生じうるゼロ除算（ｚｅｒｏ  ｄｉｖｉｄｅ）を避けるためである。
Ｃ（ｆ）は、例えば次式（３）で示される周波数特性を示す。
【００２６】
【数３】
2011180590.tif 000005
【００２７】
Ａは、式（２）の右辺第２項の成分について右辺第１項の成分とのパワー比（例えば５％）を示す定数である。αは、パワーの減少率を示す定数である。αは、例えば、ｌｏｇ１０（Ｃ（ｆ＋Δｆ）／Ｃ（ｆ））＝−０．６となる値である。但し、Δｆは、帯域幅１ｋＨｚに対応する周波数領域信号を構成する周波数変数（ｆｒｅｑｕｅｎｃｙ  ｂｉｎ）の数である。
周波数特性合成部１０５１は、生成した楽譜スペクトログラムＸ’ｆ，ｋを楽譜位置推定部１０６に出力する。
【００２８】
クロマベクトル生成部１０５２は、楽譜情報入力部１０４から入力された楽譜情報ｎｋに基づき楽譜クロマベクトルｃｋｓを生成する。クロマベクトル生成部１０５２は、例えば、楽譜情報ｎｋに含まれる音階ｎに対応する要素ｃｋｓ（ｎ）を１とし、楽譜情報ｎｋに含まれない音階ｎに対応する要素ｃｋｓ（ｎ）をゼロとした楽譜クロマベクトルを生成する。例えば、ｎ２が［Ｃ，Ｅ，Ｇ］であるとき、クロマベクトル生成部１０５２は、楽譜クロマベクトルｃ２ｓ＝［１，０，０，０，１，０，０，１，０，０，０，０］を生成する。
クロマベクトル生成部１０５２は、生成した楽譜クロマベクトルｃｋｓを楽譜位置推定部１０６に出力する。
【００２９】
楽譜位置推定部１０６は、重み（ｗｅｉｇｈｔ）算出部１０６１、再標本化（ｒｅｓａｍｐｌｉｎｇ）部１０６２、状態遷移（ｓｔａｔｅ  ｔｒａｎｓｉｔｉｏｎ）部１０６３、楽譜位置算出部１０６４及びテンポ算出部１０６５を含んで構成される。
重み算出部１０６１は、音響特徴量生成部１０２から入力された音響特徴量と楽譜特徴量に基づいて後述する粒子毎に重み係数を算出する。重み算出部１０６１は、算出した粒子毎の重み係数を再標本化部１０６２に出力する。
【００３０】
具体的には、重み算出部１０６１は、クロマベクトル重み（ｃｈｒｏｍａ  ｖｅｃｔｏｒ  ｗｅｉｇｈｔ）算出部１０６１１、周波数特性重み（ｓｐｅｃｔｒｏｇｒａｍ  ｗｅｉｇｈｔ）算出部１０６１２、テンポ重み（ｔｅｍｐｏ  ｗｅｉｇｈｔ）算出部１０６１３、及び粒子重み算出部１０６１４を含んで構成される。
【００３１】
クロマベクトル重み算出部１０６１１は、クロマベクトル生成部１０２２から入力された音響クロマベクトルｃｔａとクロマベクトル生成部１０５２から入力された楽譜クロマベクトルｃｋｔｓに基づいてクロマベクトル重み係数ｗｃｈを算出する。クロマベクトル重み係数ｗｃｈは、音響クロマベクトルｃｔａと楽譜クロマベクトルｃｋｓとの類似度又は相関性を示す尺度である。クロマベクトル重み算出部１０６１１は、例えば、以下の式を用いてクロマベクトル重み係数ｗｃｈを算出する。
【００３２】
【数４】
2011180590.tif 000006
【００３３】
式（４）において、τは時刻を表す。Ｌは予め設定された時間を表す。ｔは最新のフレーム時刻を表す。ｋτは、観測時刻τにおける楽譜フレームｋτである。ｋτと観測時刻τは次式（５）の関係がある。
【００３４】
【数５】
2011180590.tif 000007
【００３５】
式（５）において、ｆｌｏｏｒ（．．．）は、ｆｌｏｏｒ関数（ｆｌｏｏｒ  ｆｕｎｃｔｉｏｎ；床関数）を表す。即ち、ｆｌｏｏｒ（．．．）は、実数．．．を、その整数部に切り捨てる関数である。ｋ’は、仮想楽譜位置である。ｂ’は、仮想拍間隔（ｂｅａｔ  ｉｎｔｅｒｖａｌ）である。ｋ’及びｂ’については、後述する。
即ち、式（４）は、クロマベクトル重み係数ｗｃｈが、音響クロマベクトルｃｔａと楽譜クロマベクトルｃｋｔｓとの内積を最新のフレームｔから予め設定された時間Ｌだけ遡って累積した値であることを表す。以下の説明では、この時間Ｌを観測時間（窓長；ｗｉｎｄｏｗ  ｌｅｎｇｔｈ）と呼ぶ。なお、Ｌは、例えば２．５秒である。
なお、式（４）によれば、クロマベクトル重み係数ｗｃｈは、音響クロマベクトルｃｔａと楽譜クロマベクトルｃｋｔｓが類似する（つまり、両ベクトルの方向が近くなる）ほど増加し、音響クロマベクトルｃｔａと楽譜クロマベクトルｃｋｔｓが一致するとき最大値１となる。
クロマベクトル重み算出部１０６１１は、算出したクロマベクトル重み係数ｗｃｈを粒子重み算出部１０６１４に出力する。
【００３６】
なお、クロマベクトル重み算出部１０６１１は、後述する粒子（ｐａｒｔｉｃｌｅ）ｉ毎に定められた仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’が、状態遷移部１０６３から入力される。クロマベクトル重み算出部１０６１１は、観測時間の始点ｔ−Ｌと終点ｔを、各々先行時間（ｌｏｏｋ  ａｈｅａｄ  ｔｉｍｅ）ΔＴだけ後の時刻に変更する。先行時間ΔＴとは、状態遷移部１０６３において１回の状態遷移毎に観測時刻を進める時間であり、例えば１秒である。これにより、クロマベクトル重み算出部１０６１１は、観測時間長Ｌを変更せずに観測時間をシフトさせる。また、クロマベクトル重み算出部１０６１１は、上述の処理において粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’に基づきクロマベクトル重み係数ｗｉｃｈを算出及び出力する。
【００３７】
周波数特性重み算出部１０６１２は、周波数特性分析部１０２１から入力された音響スペクトログラムＸｆ，ｔと周波数特性合成部１０５１から入力された楽譜スペクトログラムＸ’ｆ，ｋに基づいて周波数特性重み係数ｗｓｐを算出する。周波数特性重み係数ｗｓｐは、音響スペクトログラムＸｆ，ｔと楽譜スペクトログラムＸ’ｆ，ｋとの類似度又は相関性を示す尺度である。周波数特性重み算出部１０６１２は、例えば、音響スペクトログラムＸｆ，ｔと楽譜スペクトログラムＸ’ｆ，ｋのカルバック・ライブラー情報量（Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ  Ｄｉｖｅｒｇｅｎｃｅ）に基づいて周波数特性重み係数ｗｓｐと定める。ここで、周波数特性重み算出部１０６１２は、次式（６）を用いて周波数特性重み係数ｗｓｐを算出する。
【００３８】
【数６】
2011180590.tif 000008
【００３９】
式（６）において、Ｑは、以下の式で表される実数である。
【００４０】
【数７】
2011180590.tif 000009
【００４１】
但し、（７）式において、周波数特性重み算出部１０６１２は、次式（８）に示すように規格化した音響スペクトログラムＸｆ，ｔと楽譜スペクトログラムＸ’ｆ，ｋｔを用いる。
【００４２】
【数８】
2011180590.tif 000010
【００４３】
また、周波数特性重み算出部１０６１２においても、式（５）に示される時刻τと楽譜フレームｋτとの関係を用いる。
式（６）及び（７）によれば、周波数特性重み係数ｗｓｐは、音響スペクトログラムＸｆ，ｔと楽譜スペクトログラムＸ’ｆ，ｋｔが類似するほど増加し、両者が一致するとき、最大値１となる。このように、周波数特性重み算出部１０６１２は、周波数特性重み係数ｗｓｐを用いることで、クロマベクトルだけでは表現されない、入力音響信号や楽譜情報による周波数特性の広がりを考慮して確実に入力音響信号と楽譜情報との対応関係を把握することができる。
周波数特性重み算出部１０６１２は、算出した周波数特性重み係数ｗｓｐを粒子重み算出部１０６１４に出力する。
【００４４】
なお、周波数特性重み算出部１０６１２は、粒子ｉ毎に定められた仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’が、状態遷移部１０６３から入力される。このとき、周波数特性重み算出部１０６１２は、観測時間の始点ｔ−Ｌと終点ｔを、各々先行時間ΔＴだけ後の時刻に変更することで、観測時間を後にシフトさせる。また、周波数特性重み算出部１０６１２は、上述の処理において粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’に基づき周波数特性重み係数ｗｉｓｐを算出及び出力する。
【００４５】
テンポ重み算出部１０６１３は、周波数特性分析部１０２１から入力された音響スペクトログラムＸｆ，ｔからテンポ重み係数ｗｔを算出する。テンポ重み係数ｗｔは、入力音響信号が示す仮想拍間隔ｂ’の確からしさを示す尺度である。
テンポ重み算出部１０６１３は、例えば、音響スペクトログラムＸｆ，ｔと仮想拍間隔ｂ’だけずれた音響スペクトログラムＸｆ，ｔの正規化相互相関（ｎｏｒｍａｌｉｚｅｄ  ｃｒｏｓｓ  ｃｏｒｒｅｌａｔｉｏｎ）をテンポ重み係数ｗｔとして算出する。即ち、テンポ重み算出部１０６１３は、次式（９）を用いてテンポ重み係数ｗｔを算出する。
【００４６】
【数９】
2011180590.tif 000011
【００４７】
式（９）において、仮想拍間隔ｂ’が先行時間ΔＴと等しいか、より大きいとき、テンポ重み算出部１０６１３は、テンポ重み係数ｗｔをゼロと定める。
式（９）によれば、テンポ重み係数ｗｔは、現在の音響スペクトログラムと過去の音響スペクトログラムの相関性を表す。例えば、同じパワーの大きさで同じ音階が拍間隔ｂ’ごとに繰り返し演奏されている入力音響信号については、テンポ重み係数ｗｔは、１となる。
テンポ重み算出部１０６１３は、算出したテンポ重み係数ｗｔを粒子重み算出部１０６１４に出力する。
【００４８】
なお、テンポ重み算出部１０６１３は、後述する粒子毎に定められた仮想拍間隔ｂｉ’が、状態遷移部１０６３から入力される。このとき、テンポ重み算出部１０６１３は、観測時間の始点ｔ−Ｌと終点ｔを、各々先行時間ΔＴだけ後の時刻に変更することで、観測時間がシフトさせる。また、テンポ重み算出部１０６１３は、上述の処理において粒子ｉ毎の仮想拍間隔ｂｉ’に基づきテンポ重み係数ｗｉｔを算出及び出力する。
【００４９】
粒子重み算出部１０６１４は、クロマベクトル重み算出部１０６１１から入力されたクロマベクトル重み係数ｗｉｃｈ、周波数特性重み算出部１０６１２から入力された周波数特性重み係数ｗｉｓｐ、及びテンポ重み算出部１０６１３から入力されたテンポ重み係数ｗｉｔに基づき粒子ｉ毎の重み係数ｗｉを算出する。
粒子ｉとは、本実施形態において仮想的な変数である仮想楽譜位置ｋｉ’、仮想拍間隔ｂｉ’及び観測値である重み係数ｗｉ、の組である。これらの粒子は、粒子フィルタリング（ｐａｒｔｉｃｌｅ  ｆｉｌｔｅｒｉｎｇ）法により観測値から仮想的な変数を推定するために用いられる。例えば、本実施形態は、粒子ｉ毎の仮想楽譜位置ｋｉ’、仮想拍間隔ｂｉ’に基づき、楽譜位置ｋ’、拍間隔ｂ’の推定値を定める。
具体的には、粒子重み算出部１０６１４は、粒子ｉ毎の粒子重み係数ｗｉを次式（１０）に基づいて算出する。
【００５０】
【数１０】
2011180590.tif 000012
【００５１】
粒子重み算出部１０６１４は、算出した粒子ｉ毎の仮想楽譜位置ｋｉ’、仮想拍間隔ｂｉ’及び粒子重み係数ｗｉを再標本化部１０６２に出力する。
【００５２】
再標本化部１０６２は、粒子重み算出部１０６１４から粒子ｉ毎の仮想楽譜位置ｋｉ’、仮想拍間隔ｂｉ’及び粒子重み係数ｗｉが入力される。再標本化部１０６２は、粒子ｉ毎の粒子重み係数ｗｉに基づき、各粒子を再標本化する。再標本化部１０６２は、再標本化した粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’を状態遷移部１０６３に出力する。再標本化部１０６２は、再標本化に係る処理において粒子数を一定に保つために、棄却する粒子数と、分割又は複製によって新たに生成する粒子数を等しくする。
【００５３】
再標本化部１０６２は、例えばＳＩＲ（ｓａｍｐｌｉｎｇ  ｉｍｐｏｒｔａｎｃｅ  ｒｅｓａｍｐｌｉｎｇ）法を用いて粒子を再標本化する。
ここで、再標本化部１０６２は、次の処理を行って粒子重みｗｉが小さい粒子ｉを棄却する。このような粒子ｉに対応する、仮想楽譜位置ｋｉ’や仮想拍間隔ｂｉ’の信頼性が低いからである。
例えば、再標本化部１０６２は、粒子ｉ毎の粒子重み係数を、次式（１１）を用いて規格化し、規格化粒子重み係数ｐｉを算出する。
【００５４】
【数１１】
2011180590.tif 000013
【００５５】
式（１１）において、Ｎは粒子数（例えば、３００個）を示す。
再標本化部１０６２は、規格化粒子重み係数ｐｉが予め設定された規格化粒子重み係数の閾値ｐｔｈよりも小さいか否か判断する。再標本化部１０６２は、規格化粒子重み係数ｐｉが閾値ｐｔｈよりも小さいと判断した場合、対応する粒子ｉを棄却する。即ち、再標本化部１０６２は、粒子ｉの規格化粒子重み係数ｐｉ、仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’を削除する。再標本化部１０６２は、規格化粒子重み係数ｐｉが閾値ｐｔｈよりも大きい又は等しいと判断した場合、対応する粒子ｉを維持する。また、再標本化部１０６２は、棄却された粒子数をカウントする。
【００５６】
再標本化部１０６２は、次の処理を行って粒子重みｗｉが大きい粒子ｉを分割又は複製する。このような粒子ｉに対応する、仮想楽譜位置ｋｉ’や仮想拍間隔ｂｉ’の信頼性が高いからである。
例えば、再標本化部１０６２は、粒子ｉ毎の粒子重み係数を、次式（１２）を用いて再規格化し、再規格化粒子重み係数ｑｉを算出する。
【００５７】
【数１２】
2011180590.tif 000014
【００５８】
再標本化部１０６２は、再規格化粒子重み係数ｑｉが1よりも大きい粒子について整数部分の値ｎｑｉ（以下、再規格化整数と呼ぶ）と小数部分の値ｓｑｉ（以下、再規格化小数と呼ぶ）を算出する。
再標本化部１０６２は、再規格化整数ｎｑｉが２又は２より多い粒子ｉを分割して、新たにｎｑｉ−１個の粒子ｊ（ｊ≠ｉ）を各々生成する。即ち、再標本化部１０６２は、新たな粒子ｊに対応する仮想楽譜位置ｋｊ’及び仮想拍間隔ｂｊ’を分割直前の仮想楽譜位置ｋｉ‘及び仮想拍間隔ｂｉ’と定める。再標本化部１０６２は、棄却された粒子数から新たに生成した粒子数ｎｑｉ−１を順次差し引くことにより、所定の粒子数Ｎに不足する粒子数Ｎｓ（以下、不足粒子数と呼ぶ）を算出する。
不足粒子数がゼロに達した場合には、再標本化部１０６２は、後述する粒子を複製する処理を行わない。他方、不足粒子数が１以上である場合に、再標本化部１０６２は、粒子を複製する処理を行う。
そして、再標本化部１０６２は、新たに複製される粒子の総数が不足粒子数Ｎｓとなるように、再規格化小数ｓｑｉに比例する確率（以下、複製確率と呼ぶ）で、複製対象の粒子ｉを決定する。例えば、再標本化部１０６２は、次式（１３）で算出した複製確率ｓ’ｑｉが、乱数ｒｉ（但し、０＜ｒｉ＜１）よりも小さい粒子を複製対象の粒子と定める。
【００５９】
【数１３】
2011180590.tif 000015
【００６０】
再標本化部１０６２は、定めた粒子各々について、新たに１個の粒子（ｊ≠ｉ）を複製する。即ち、再標本化部１０６２は、新たな粒子ｊに対応する仮想楽譜位置ｋｊ’及び仮想拍間隔ｂｊ’を複製直前の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’と定める。
再標本化部１０６２は、上述の処理を完了して決定した粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’を状態遷移部１０６３に出力する。
【００６１】
状態遷移部１０６３は、再標本化部１０６２から粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’を入力する。状態遷移部１０６３は、仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’を、例えば式（１４）及び式（１５）を用いて更新する。
【００６２】
【数１４】
2011180590.tif 000016
【００６３】
【数１５】
2011180590.tif 000017
【００６４】
式（１４）及び（１５）において、ｕ、ｖは、各々平均値がゼロ、分散が各々σｕ、σｖとなるガウシアン乱数（ｇａｕｓｓｉａｎ  ｒａｎｄｏｍ  ｖａｒｉａｂｌｅ）である。例えば、σｕは１．０であり、σｖは０．２である。
つまり、式（１４）は、状態遷移部１０６３が、進行時間ΔＴ中の仮想拍間隔ｂｉ’（テンポ）に応じた楽譜位置の進行に一定範囲に分布するランダムな揺らぎを示す乱数ｕを加えて仮想楽譜位置ｋｉ‘を更新していることを表す。
他方、式（１５）は、状態遷移部１０６３が、一定範囲に分布するランダムな揺らぎを示す乱数ｖを加えて仮想拍間隔ｂｉ’を更新していることを表す。
このように、本実施形態は、粒子ｉごとに独立な状態遷移を行い、状態遷移において一定範囲内に分布する乱数を加えることで、演奏される音楽のテンポの揺らぎを考慮することができる。
状態遷移部１０６３は、更新した粒子ｉ毎の仮想楽譜位置ｋｉ’を重み算出部１０６１及び楽譜位置算出部１０６４に出力し、更新した粒子ｉ毎の仮想拍間隔ｂｉ’を重み算出部１０６１及びテンポ算出部１０６５に出力する。
【００６５】
楽譜位置算出部１０６４は、状態遷移部１０６３から入力された粒子ｉ毎の仮想楽譜位置ｋｉ’に基づいて推定楽譜位置μを算出する。楽譜位置算出部１０６４は、例えば、仮想楽譜位置ｋｉ’の単純平均をμとする。本実施形態では、これに限られず、例えば楽譜位置算出部１０６４は、仮想楽譜位置ｋｉ’の偏差が所定の閾値ｕＴＨよりも小さくなる粒子間の平均値を推定楽譜位置μと算出してもよい。これにより、楽譜位置算出部１０６４は、平均値を算出する際、稠密に分布している粒子だけを考慮し、極端に分散した粒子を除外できる。
楽譜位置算出部１０６４は、粒子ｉ毎の仮想楽譜位置ｋｉ’と算出した推定楽譜位置μを信頼度算出部１０７に出力し、算出した推定楽譜位置μを楽譜位置出力部１０８に出力する。
【００６６】
信頼度算出部１０７は、楽譜位置算出部１０６４から入力された粒子ｉ毎の仮想楽譜位置ｋｉ’と算出した推定楽譜位置μに基づき、推定楽譜位置μの信頼度係数を算出する。信頼度算出部１０７は、例えば、次式（１６）を用いて最新の時刻ｔにおける信頼度係数ｓ２（ｔ）を算出する。この信頼度係数ｓ２（ｔ）は、粒子ｉ間の仮想楽譜位置ｋｉ’の分散であり、値が大きくなるほど推定楽譜位置の信頼性が低下することを示す。
【００６７】
【数１６】
2011180590.tif 000018
【００６８】
ここで、信頼度算出部１０７は、最新の時刻における信頼度係数ｓ２（ｔ）と予め定めた遷移回数Ｉ（例えば５回）前の時刻ｔ−ＩΔＴにおける信頼度係数ｓ２（ｔ−ＩΔＴ）の差分（以下、信頼度差分と呼ぶ）が予め定めた閾値γｉｎｃＩを超えるか否かを判断する。また、信頼度算出部１０７は、信頼度差分が予め定めた閾値−γｄｅｃＩを下回るか否かを判断する。γｉｎｃ及びγｄｅｃ任意の実数であり、例えばともに４である。
信頼度算出部１０７は、信頼度差分が閾値γｉｎｃＩを超えると判断したとき、楽譜位置の出力を停止することを指示する楽譜位置出力停止信号を楽譜位置出力部１０８に出力する。信頼度算出部１０７は、信頼度差分が閾値−γｄｅｃＩを下回ると判断したとき、楽譜位置の出力を開始することを指示する楽譜位置出力開始信号を楽譜位置出力部１０８に出力する。
【００６９】
上述の例では、信頼度算出部１０７は、信頼度差分に基づき楽譜位置出力停止信号又は楽譜位置出力開始信号の出力を制御していたが、本実施形態ではこれに限らない。例えば、信頼度算出部１０７は、信頼度係数ｓ２（ｔ）自体が予め定めた閾値γｓを超えるか否かを判断し、信頼度係数ｓ２（ｔ）が閾値γｓを超えると判断したとき、楽譜位置出力停止信号を楽譜位置出力部１０８に出力し、信頼度係数ｓ２（ｔ）が閾値γｓを下回ると判断したとき、楽譜位置出力開始信号を楽譜位置出力部１０８に出力してもよい。
信頼度算出部１０７は、上記の信頼度係数の代わりに推定楽譜位置μの推定誤差ｅ（ｔ）の絶対値｜ｅ（ｔ）｜もしくは二乗誤差ｅ２（ｔ）を信頼度係数として算出し、算出した信頼度が予め定めた閾値γｅを超えるか否かを判断し、信頼度係数が閾値γｅを超えると判断したとき、楽譜位置出力停止信号を楽譜位置出力部１０８に出力し、信頼度係数が閾値γｅを下回ると判断したとき、楽譜位置出力開始信号を楽譜位置出力部１０８に出力してもよい。
ここで、推定誤差ｅ（ｔ）は、次式（１７）で表される。
【００７０】
【数１７】
2011180590.tif 000019
【００７１】
式（１７）において、ｓ（ｋ’（ｔ））は、フレームｋ’（ｔ）における真の（ｇｒｏｕｎｄ  ｔｒｕｔｈ）時刻である。ｋ’（ｔ）は、推定楽譜位置μである。
【００７２】
楽譜位置出力部１０８は、楽譜位置算出部１０６４から推定楽譜位置μが入力され、入力された推定楽譜位置μを出力する。
但し、楽譜位置出力部１０８は、信頼度算出部１０７から楽譜位置出力停止信号が入力されたとき、推定楽譜位置μの出力を停止する。他方、楽譜位置出力部１０８は、信頼度算出部１０７から楽譜位置出力開始信号が入力されたとき、停止していた推定楽譜位置μの出力を再開する。
これにより、本実施形態は、信頼性が低い推定楽譜位置μの出力を回避することができる。
【００７３】
テンポ算出部１０６５は、状態遷移部１０６３から入力された粒子ｉ毎の仮想拍間隔ｂｉ’に基づき推定拍間隔βを算出する。この推定拍間隔は、演奏されている音楽のテンポを示す尺度である。テンポ算出部１０６５は、例えば、仮想拍間隔ｂｉ’の単純平均をβとする。本実施形態では、これに限られず、例えばテンポ算出部１０６５は、仮想拍間隔ｂｉ’の偏差が所定の閾値νＴＨよりも小さくなる粒子間の平均値を推定拍間隔βと算出してもよい。これにより、テンポ算出部１０６５は、平均値を算出する際に稠密に分布している粒子だけを考慮し、極端に分散した粒子を除外できる。
テンポ算出部１０６５は、推定拍間隔βをテンポ出力部１０９に出力する。
テンポ出力部１０９は、テンポ算出部１０６５から推定拍間隔βが入力され、入力された推定拍間隔βを楽譜位置推定装置１の外部に出力する。
【００７４】
ロボットは、例えば出力された推定楽譜位置μ及び推定拍間隔βを歌声の合成や機構部の位置制御に用い、合成した歌声や機構部の動作を演奏されている音楽に同期させることができる。本実施形態の楽譜位置推定処理は、例えばロボットが演奏されている音楽に反応して歌声を合成したり、動作機構を制御する、等に応用することが期待されている。これらの応用において必然的に遅延が生じる問題を、ロボットが本実施形態により推定された楽譜位置を用いることで解決できる。
【００７５】
また、楽譜位置推定装置１は、当初、推定楽譜位置μと推定拍間隔βをともに出力するが、推定楽譜位置の信頼度増分が所定の閾値を越えると推定楽譜位置μの出力を停止する状態（リズムレベル（ｒｈｙｔｈｍ  ｌｅｖｅｌ）に遷移する。推定楽譜位置の信頼度増分が所定の閾値を下回ると推定楽譜位置μの出力を再開する状態（メロディレベル（ｍｅｌｏｄｙ  ｌｅｖｅｌ）に遷移する。これにより、本実施形態では信頼性が低い推定楽譜位置が出力されなくなるため、歌声の合成や機構部の位置制御等において、利用者の意図しない動作をすることを回避できる。
【００７６】
次に、本実施形態に係る楽譜位置推定処理の概要について図を用いて説明する。
図２は、本実施形態に係る楽譜位置推定処理の概要を示す概念図である。図２左列は観測（Ｉ．）、中列は再標本化（ＩＩ．）、右列は状態遷移（ＩＩＩ．）、の過程を示す。このような過程を含む一連の処理は、粒子フィルタリングに含まれる処理である。
図２左列は、上段から下段へ順に、音響信号、楽譜情報、粒子の分布を示す。各段ともに横軸は時刻である。上段は、縦軸は振幅であり音響信号の波形を表す。中段は、五線上に音符を用いて楽譜情報を表す。下段は、円を用いて仮想楽譜位置毎の粒子の分布を表す。円の半径は、各粒子の重み係数の大きさを表す。
上段と中段に示されている四角の枠は、音響信号と楽譜情報を観測対象とする観測時間を示す。
従って、図２左列は、観測（Ｉ．）の過程において、音響特徴量と楽譜特徴量に基づき、粒子毎に重み係数を算出することを示す。
【００７７】
図２中列は、中段に再標本化前の粒子の分布を、下段に再標本化後の粒子の分布を示す。ここで、再標本化前の粒子の分布は、同図左列下段に示されているものと同一である。再標本化前の粒子から下段に向かう破線の矢印と×印は、その粒子が棄却されることを示す。再標本化前の粒子から再標本化後の粒子に向かう細い実線の矢印は、その粒子が維持されることを示す。再標本化前の粒子から再標本化後の粒子に向かう太い実線の矢印は、その粒子が分割又は複製され、同一の仮想楽譜情報を有する粒子が生成されることを示す。
従って、図２中列は、再標本化において重み係数が小さい粒子が棄却され、重み係数が大きい粒子について分割又は複製されることにより棄却された数と同数の新たな粒子が生成されることを示す。
【００７８】
図２右列は、上段から下段へ順に、楽譜情報、遷移前の粒子の分布、遷移後の粒子の分布を示す。ここで、楽譜情報は、同図左列中段に示されているものと同一である。遷移前の粒子の分布は、同図中列下段に示されているものと同一である。中段において下に向かう矢印は、粒子の分布の中心が現在の楽譜位置であることを示す。下段において上に向かう矢印は、粒子の分布の中心が推定楽譜位置であることを示す。ここで、図２右列は、現在の楽譜位置と推定楽譜位置の間が、拍間隔であることを示す。中段の粒子から下段に向かう矢印は、その矢印が示す楽譜位置に粒子が遷移することを示す。
従って、図２右列は、状態遷移において粒子毎の仮想楽譜位置が更新され、更新後の仮想楽譜位置の分布に基づき楽譜位置が推定されることを示す。
【００７９】
図３は、本実施形態に係る楽譜位置推定処理の概要を示すフローチャートである。
（ステップＳ１）楽譜位置推定装置１の各構成部は、処理に用いる変数やデータを初期設定する。
（ステップＳ２）音響信号入力部１０１は、音響信号を入力し、音響信号の入力が継続するか否か判断する。音響信号入力部１０１が、音響信号の入力が継続すると判断したとき（ステップＳ２  Ｙ）、ステップＳ３に進む。音響信号入力部１０１が、音響信号の入力が継続していないと判断したとき（ステップＳ２  Ｎ）、処理を終了する。
（ステップＳ３）音響特徴量生成部１０２は、入力された音響信号に基づき特徴量を生成する。重み算出部１０６１は、生成した音響信号に基づく特徴量と楽譜情報に基づく特徴量から粒子毎の重み係数を算出する。（Ｉ．観測）
【００８０】
（ステップＳ４）再標本化部１０６２は、重み係数が所定の値より小さい粒子を棄却し、重み係数が所定の値より大きい粒子を分割又は複製する。（ＩＩ．再標本化）
（ステップＳ５）状態遷移部１０６３は、粒子毎の仮想楽譜位置及び仮想拍間隔を更新する。楽譜位置算出部１０６４は、更新後の仮想楽譜位置に基づき推定楽譜位置を算出する。テンポ算出部１０６５は、更新後の仮想拍間隔に基づき推定拍間隔を算出する。（ＩＩＩ．状態遷移）
【００８１】
次に、本実施形態に係る楽譜位置推定処理の詳細について図を用いて説明する。
図４は、本実施形態に係る楽譜位置推定処理の一部を示すフローチャートである。
図４は、初期設定（ステップＳ１）、音響信号入力の継続判断（ステップＳ２）及び観測（ステップＳ３）の各過程を含んで構成される。
初期設定（ステップＳ１）は、ステップＳ１０１−Ｓ１０４を含んで構成される。
（ステップＳ１０１）楽譜位置推定装置１の各構成部は、処理に用いる変数やデータを初期設定する。例えば、重み算出部１０６１は、粒子ｉ毎の仮想拍間隔ｂｉ‘、粒子ｉ毎の仮想楽譜位置ｋｉ‘、分散σｕ、σｖを初期設定する。
重み算出部１０６１は、粒子ｉ毎の仮想拍間隔ｂｉ‘の初期値を、例えば各々テンポｘに対応する拍間隔６０／ｘに定める。このテンポｘは、６０ｂｐｍ（ｂｅａｔ  ｐｅｒ  ｍｉｎｕｔｅ；１分あたりの拍数）から１８０ｂｐｍの間で均一に分布したものである。
重み算出部１０６１は、粒子ｉ毎の仮想楽譜位置ｋｉ‘の初期値を、例えば、ΔＴ／ｂｉ‘と定める。
重み算出部１０６１は、先行時間ΔＴを、例えば１秒と設定しておく。その後、ステップＳ１０２に進む。
（ステップＳ１０２）楽譜情報入力部１０４は、楽譜情報記憶部１０３から、処理の対象とする楽曲の楽譜情報ｎｋを読み出す。楽譜情報入力部１０４は、読み出した楽譜情報ｎｋを楽譜特徴量生成部１０５に出力する。その後、ステップＳ１０３に進む。
【００８２】
（ステップＳ１０３）楽譜特徴量生成部１０５が備える周波数特性合成部１０５１は、楽譜情報入力部１０４から入力された楽譜情報ｎｋに基づいて楽譜特徴量として、例えば楽譜スペクトログラムＸ’ｆ，ｋを算出する。周波数特性合成部１０５１は、楽譜スペクトログラムＸ’ｆ，ｋを算出する際、例えば式（２）に表されるｈａｒｍｏｎｉｃ  ＧＭＭを用いる。
周波数特性合成部１０５１は、その楽曲を構成する楽譜情報ｎｋに全体について楽譜スペクトログラムＸ’ｆ，ｋを算出しておいてもよい。また、予め算出した楽譜スペクトログラムＸ’ｆ，ｋを楽譜情報記憶部１０３に記憶しておき、周波数特性合成部１０５１は、楽譜情報記憶部１０３から楽譜スペクトログラムＸ’ｆ，ｋを読み出してもよい。
周波数特性合成部１０５１は、生成または読み出した楽譜スペクトログラムＸ’ｆ，ｋを楽譜位置推定部１０６に出力する。その後、ステップＳ１０４に進む。
【００８３】
（ステップＳ１０４）楽譜特徴量生成部１０５が備えるクロマベクトル生成部１０５２は、楽譜情報入力部１０４から入力された楽譜情報ｎｋに基づいて楽譜特徴量として、楽譜クロマベクトルｃｋｓを生成する。
クロマベクトル生成部１０５２は、その楽曲を構成する楽譜情報ｎｋ全てについて楽譜クロマベクトルｃｋｓを生成しておいてもよい。また、予め生成した楽譜クロマベクトルｃｋｓを楽譜情報記憶部１０３に記憶しておき、クロマベクトル生成部１０５２は、楽譜情報記憶部１０３から楽譜クロマベクトルｃｋｓを読み出してもよい。
クロマベクトル生成部１０５２は、生成又は読み出した楽譜クロマベクトルｃｋｓを楽譜位置推定部１０６に出力する。その後、ステップＳ２に進む。
【００８４】
（ステップＳ２）音響信号入力部１０１は、音響信号を入力し、音響信号の入力が継続するか否か判断する。音響信号入力部１０１が、音響信号の入力が継続すると判断したとき（ステップＳ２  Ｙ）、ステップＳ３０１に進む。音響信号入力部１０１が、音響信号の入力が継続していないと判断したとき（ステップＳ２  Ｎ）、処理を終了する。
【００８５】
観測（ステップＳ３）の過程は、ステップＳ３０１−Ｓ３０７を含んで構成される。
（ステップＳ３０１）音響信号入力部１０１は、アナログ音響信号を入力し、入力されたアナログ音響信号をディジタル信号に変換する。音響信号入力部１０１は、変換されたディジタル信号を音響特徴量生成部１０２に出力する。その後、ステップＳ３０２に進む。
（ステップＳ３０２）周波数特性分析部１０２１は、音響信号入力部１０１から入力されたディジタル音響信号を時間領域から周波数領域信号に変換する。周波数特性分析部１０２１は、変換された周波数領域信号の振幅を音響スペクトログラムＸｆ，ｔとして算出し、算出した音響スペクトログラムＸｆ，ｔをクロマベクトル生成部１０２２及び楽譜位置推定部１０６に出力する。その後、ステップＳ３０３に進む。
（ステップＳ３０３）クロマベクトル生成部１０２２は、周波数特性分析部１０２１から入力された音響スペクトログラムＸｆ，ｔから音響クロマベクトルを算出する。ここで、クロマベクトル生成部１０２２は、例えば、音響クロマベクトルｃｔａの音階ｎ毎の要素ｃｔａ（ｎ）を、式（１）を用いて算出する。
クロマベクトル生成部１０２２は、算出した音響クロマベクトルｃｔａを楽譜位置推定部１０６に出力する。その後、ステップＳ３０４に進む。
【００８６】
（ステップＳ３０４）クロマベクトル重み算出部１０６１１は、クロマベクトル生成部１０２２から入力された音響クロマベクトルｃｔａとクロマベクトル生成部１０５２から入力された楽譜クロマベクトルｃｋｔｓからクロマベクトル重み係数ｗｃｈを算出する。ここで、クロマベクトル重み算出部１０６１１は、例えば式（４）及び（５）を用い、状態遷移部１０６３から入力された粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’に基づき、粒子ｉ毎のクロマベクトル重み係数ｗｉｃｈを算出する。
クロマベクトル重み算出部１０６１１は、算出した粒子ｉ毎のクロマベクトル重み係数ｗｉｃｈを粒子重み算出部１０６１４に出力する。その後、ステップＳ３０５に進む。
【００８７】
（ステップＳ３０５）周波数特性重み算出部１０６１２は、周波数特性分析部１０２１から入力された音響スペクトログラムＸｆ，ｔと周波数特性合成部１０５１から入力された楽譜スペクトログラムＸ’ｆ，ｋから周波数特性重み係数ｗｓｐを算出する。ここで、周波数特性重み算出部１０６１２は、例えば式（５）−（８）を用い、状態遷移部１０６３から入力された粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’に基づき、粒子ｉ毎の周波数特性重み係数ｗｓｐを算出する。
周波数特性重み算出部１０６１２は、算出した粒子ｉ毎の周波数特性重み係数ｗｉｓｐを粒子重み算出部１０６１４に出力する。その後、ステップＳ３０６に進む。
【００８８】
（ステップＳ３０６）テンポ重み算出部１０６１３は、周波数特性分析部１０２１から入力された音響スペクトログラムＸｆ，ｔからテンポ重み係数ｗｔを算出する。ここで、テンポ重み算出部１０６１３は、例えば式（９）を用いて状態遷移部１０６３から入力された粒子ｉ毎の仮想拍間隔ｂｉ’に基づき、粒子ｉ毎のテンポ重み係数ｗｔを算出する。
テンポ重み算出部１０６１３は、算出した粒子ｉ毎のテンポ重み係数ｗｉｔを粒子重み算出部１０６１４に出力する。その後、ステップＳ３０７に進む。
【００８９】
（ステップＳ３０７）粒子重み算出部１０６１４は、クロマベクトル重み算出部１０６１１から入力されたクロマベクトル重み係数ｗｉｃｈ、周波数特性重み算出部１０６１２から入力された周波数特性重み係数ｗｉｓｐ、及びテンポ重み算出部１０６１３から入力されたテンポ重み係数ｗｉｔから粒子ｉ毎の重み係数ｗｉを算出する。ここで、粒子重み算出部１０６１４は、重み係数ｗｉを算出する際、例えば式（１０）を用いる。
粒子重み算出部１０６１４は、算出した粒子ｉ毎の仮想楽譜位置ｋｉ’、仮想拍間隔ｂｉ’及び粒子重み係数ｗｉを再標本化部１０６２に出力する。その後、ステップＳ４１１に進む。
【００９０】
図５は、本実施形態に係る楽譜位置推定処理の他の一部を示すフローチャートである。図５は、再標本化（ステップＳ４）の過程を含んで構成される。再標本化（ステップＳ４）は、棄却（ステップＳ４１）、分割（ステップＳ４２）及び複製（ステップＳ４３）の各過程を含んで構成される。
棄却（ステップＳ４１）の過程は、ステップＳ４１１−Ｓ４１４を含んで構成される。
（ステップＳ４１１）再標本化部１０６２は、粒子重み算出部１０６１４から粒子ｉ毎の仮想楽譜位置ｋｉ’、仮想拍間隔ｂｉ’及び粒子重み係数ｗｉを入力する。
再標本化部１０６２は、全ての粒子ｉについてステップＳ４１２−Ｓ４１４の処理を実行したか否か判断する。再標本化部１０６２は、全ての粒子ｉについてステップＳ４１２−Ｓ４１４の処理を実行したと判断したとき（ステップＳ４１１  Ｙ）、ステップＳ４２１に進む。再標本化部１０６２は、全ての粒子ｉについてステップＳ４１２−Ｓ４１４の処理を実行していないと判断したとき（ステップＳ４１１  Ｎ）、ステップＳ４１２に進む。
【００９１】
（ステップＳ４１２）再標本化部１０６２は、例えば、式（１１）を用いて、粒子ｉ毎の粒子重み係数ｗｉから規格化粒子重み係数ｐｉを算出する。再標本化部１０６２は、規格化粒子重み係数ｐｉが予め設定された閾値ｐｔｈよりも小さいか否か判断する。再標本化部１０６２は、規格化粒子重み係数ｐｉが閾値ｐｔｈよりも小さいと判断した場合（ステップＳ４１１  Ｙ）、ステップＳ４１３に進む。再標本化部１０６２は、規格化粒子重み係数ｐｉが閾値ｐｔｈよりも大きい又は等しいと判断した場合（ステップＳ４１１  Ｎ）、対応する粒子ｉを維持し、ステップＳ４１１に進む。
（ステップＳ４１３）再標本化部１０６２は、規格化粒子重み係数ｐｉが閾値ｐｔｈよりも小さいと判断された粒子ｉを棄却する。その後、ステップＳ４１４に進む。
（ステップＳ４１４）再標本化部１０６２は、棄却された粒子数を１個増加させ、棄却された粒子数をカウントする。その後、ステップＳ４１１に進む。
【００９２】
分割（ステップＳ４２）の過程は、ステップＳ４２１−Ｓ４２６を含んで構成される。
（ステップＳ４２１）再標本化部１０６２は、棄却されずに残った粒子ｉについてステップＳ４２２−Ｓ４２６の処理を実行したか否か判断する。再標本化部１０６２は、全ての残った粒子ｉについてステップＳ４２２−Ｓ４２６の処理を実行したと判断したとき（ステップＳ４２１  Ｙ）、ステップＳ４３１に進む。再標本化部１０６２は、全ての粒子ｉについてステップＳ４２２−Ｓ４２６の処理を実行していないと判断したとき（ステップＳ４２１  Ｎ）、ステップＳ４２２に進む。
（ステップＳ４２２）再標本化部１０６２は、例えば、式（１２）を用いて、粒子ｉ毎の規格化粒子重み係数ｐｉから再規格化粒子重み係数ｑｉを算出する。再標本化部１０６２は、再規格化粒子重み係数ｑｉが１よりも大きい粒子について再規格化整数ｎｑｉと再規格化小数ｓｑｉを算出する。その後、ステップＳ４２３に進む。
【００９３】
（ステップＳ４２３）再標本化部１０６２は、再規格化整数ｎｑｉが２又は２より多いか否か判断する。再標本化部１０６２は、再規格化整数ｎｑｉが２又は２より多いと判断したとき（ステップＳ４２３  Ｙ）、ステップＳ４２４に進む。再標本化部１０６２が、再規格化整数ｎｑｉが１であると判断したとき（ステップＳ４２３  Ｎ）、ステップＳ４２１に進む。
（ステップＳ４２４）再標本化部１０６２は、粒子ｉを分割して、新たにｎｑｉ−１個の粒子ｊ（ｊ≠ｉ）を生成する。即ち、再標本化部１０６２は、新たな粒子ｊに対応する仮想楽譜位置ｋｊ’及び仮想拍間隔ｂｊ’を分割前の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’と定める。その後、ステップＳ４２５に進む。
【００９４】
（ステップＳ４２５）再標本化部１０６２は、不足粒子数Ｎｓ（初期値は棄却された粒子数）から新たに生成した粒子数ｎｑｉ−１を差し引いて不足粒子数をカウントする。その後、ステップＳ４２６に進む。
（ステップＳ４２６）再標本化部１０６２は、不足粒子数Ｎｓがゼロに達したか否か判断する。再標本化部１０６２は、不足粒子数Ｎｓがゼロに達したと判断したとき（ステップＳ４２６  Ｙ）、再標本化に関する処理を終了し、粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’を状態遷移部１０６３に出力する。その後、ステップＳ５０１に進む。不足粒子数Ｎｓがゼロに達していないと判断したとき（ステップＳ４２６  Ｎ）、ステップＳ４２１に進む。
【００９５】
複製（ステップＳ４３）の過程は、ステップＳ４３１−Ｓ４３４を含んで構成される。
（ステップＳ４３１）再標本化部１０６２は、棄却されずに残った粒子ｉについてステップＳ４３２−Ｓ４３３の処理を実行したか否か判断する。再標本化部１０６２は、全ての残った粒子ｉについてステップＳ４３２−Ｓ４３３の処理を実行したと判断したとき（ステップＳ４３１  Ｙ）、ステップＳ４３４に進む。再標本化部１０６２は、全ての粒子ｉについてステップＳ４３２−Ｓ４３３の処理を実行していないと判断したとき（ステップＳ４３１  Ｎ）、ステップＳ４３２に進む。
（ステップＳ４３２）再標本化部１０６２は、例えば式（１３）を用いて再規格化小数ｓｑｉに比例する複製確率ｓ’ｑｉを算出する。その後、ステップＳ４３３に進む。
（ステップＳ４３３）再標本化部１０６２は、複製確率ｓ’ｑｉに基づいて新たに複製される粒子の総数が不足粒子数Ｎｓとなるように、複製対象の粒子ｉを決定する。その後、ステップＳ４３１に進む。
（ステップＳ４３４）再標本化部１０６２は、定めた粒子各々について、新たに１個の粒子（ｊ≠ｉ）を複製する。即ち、再標本化部１０６２は、新たな粒子ｊに対応する仮想楽譜位置ｋｊ’及び仮想拍間隔ｂｊ’を複製前の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’と定める。
再標本化部１０６２は、粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’を状態遷移部１０６３に出力する。その後、ステップＳ５０１に進む。
【００９６】
図６は、本実施形態に係る楽譜位置推定処理の残りの一部を示すフローチャートである。図６は、状態遷移（ステップＳ５）及び推定値出力（ステップＳ６）の過程を含んで構成される。状態遷移（ステップＳ５）の過程は、ステップＳ５０１−Ｓ５０４を含んで構成される。
（ステップＳ５０１）状態遷移部１０６３は、再標本化部１０６２から粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’を入力する。状態遷移部１０６３は、進行時間ΔＴ中に仮想拍間隔ｂｉ’（テンポ）で進行した楽譜位置に乱数ｕを加えて仮想楽譜位置ｋｉ’を更新する。状態遷移部１０６３は、仮想楽譜位置ｋｉ’を更新する際、例えば式（１４）を用いる。その後、ステップＳ５０２に進む。
（ステップＳ５０２）状態遷移部１０６３は、乱数ｖを加えて仮想拍間隔ｂｉ’を更新する。状態遷移部１０６３は、仮想拍間隔ｂｉ’を更新する際、例えば式（１５）を用いる。状態遷移部１０６３は、更新した粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’を重み算出部１０６１に出力する。状態遷移部１０６３は、更新した粒子ｉ毎の仮想楽譜位置ｋｉ‘を楽譜位置算出部１０６４に出力する。状態遷移部１０６３は、更新した粒子ｉ毎の仮想拍間隔ｂｉ’をテンポ算出部１０６５に出力する。その後、ステップＳ５０３に進む。
【００９７】
（ステップＳ５０３）楽譜位置算出部１０６４は、状態遷移部１０６３から入力された粒子ｉ毎の仮想楽譜位置ｋｉ‘に基づいて推定楽譜位置μを算出する。楽譜位置算出部１０６４は、例えば、仮想楽譜位置ｋｉ‘の単純平均をμとする。楽譜位置算出部１０６４は、粒子ｉ毎の仮想楽譜位置ｋｉ‘と算出した推定楽譜位置μを信頼度算出部１０７に出力する。楽譜位置算出部１０６４は、算出した推定楽譜位置μを楽譜位置出力部１０８に出力する。その後、ステップＳ５０３に進む。
（ステップＳ５０４）テンポ算出部１０６５は、状態遷移部１０６３から入力された粒子ｉ毎の仮想拍間隔ｂｉ’に基づいて推定拍間隔βを算出する。テンポ算出部１０６５は、例えば、仮想拍間隔ｂｉ’の単純平均をβとする。
テンポ算出部１０６５は、推定拍間隔βをテンポ出力部１０９に出力する。その後、ステップＳ６０１に進む。
【００９８】
推定値出力（ステップＳ６）の過程は、ステップＳ６０１−Ｓ６０６を含んで構成される。
（ステップＳ６０１）信頼度算出部１０７は、楽譜位置算出部１０６４から入力された粒子ｉ毎の仮想楽譜位置ｋｉ‘に基づき推定楽譜位置μを算出する。信頼度算出部１０７は、仮想楽譜位置ｋｉ‘と推定楽譜位置μに基づき、推定楽譜位置μの信頼度係数ｓ２（ｔ）を算出する。信頼度算出部１０７は、信頼度係数ｓ２（ｔ）を算出する際、例えば、式（１６）を用いる。その後、ステップＳ６０２に進む。
（ステップＳ６０２）信頼度算出部１０７は、信頼度差分が予め定めた閾値γｉｎｃＩを超えるか否かを判断する。信頼度算出部１０７は、信頼度差分が予め定めた閾値γｉｎｃＩを超えると判断したとき（ステップＳ６０２  Ｙ）、ステップＳ６０３に進む。信頼度算出部１０７は、信頼度差分が予め定めた閾値γｉｎｃＩを超えないと判断したとき（ステップＳ６０２  Ｎ）、ステップＳ６０４に進む。
（ステップＳ６０３）信頼度算出部１０７は、楽譜位置出力停止信号を楽譜位置出力部１０８に出力する。ここで、楽譜位置出力部１０８は、楽譜位置算出部１０６４から推定楽譜位置μが入力され、入力された推定楽譜位置μを出力する。楽譜位置出力部１０８は、信頼度算出部１０７から楽譜位置出力停止信号が入力されたとき、推定楽譜位置μの出力を停止する。その後、ステップＳ６０４に進む。
【００９９】
（ステップＳ６０４）信頼度算出部１０７は、信頼度算出部１０７は、信頼度差分が予め定めた閾値−γｄｅｃＩを下回るか否かを判断する。信頼度算出部１０７は、信頼度差分が予め定めた閾値−γｄｅｃＩを下回ると判断したとき（ステップＳ６０４  Ｙ）、ステップＳ６０５に進む。信頼度算出部１０７は、信頼度差分が予め定め閾値−γｄｅｃＩを下回らないと判断したとき（ステップＳ６０４  Ｎ）、ステップＳ６０６に進む。
（ステップＳ６０５）信頼度算出部１０７は、楽譜位置出力開始信号を楽譜位置出力部１０８に出力する。楽譜位置出力部１０８は、信頼度算出部１０７から楽譜位置出力開始信号が入力されたとき、停止していた推定楽譜位置μの出力を再開する。その後、ステップＳ６０６に進む。
（ステップＳ６０６）テンポ出力部１０９は、テンポ算出部１０６５から推定拍間隔βが入力され、入力された推定拍間隔βを出力する。その後、ステップＳ７に進む。
【０１００】
（ステップＳ７）重み算出部１０６１は、状態遷移部１０６３から粒子ｉ毎の仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’が入力される。クロマベクトル重み算出部１０６１１は、観測時間の始点ｔ−Ｌと終点ｔを、各々先行時間ΔＴだけ後の時刻にシフトさせて観測時間を変更する。
重み算出部１０６１は、入力された仮想楽譜位置ｋｉ’、仮想拍間隔ｂｉ’及び観測時間を、クロマベクトル重み算出部１０６１１及び周波数特性重み算出部１０６１２に出力する。重み算出部１０６１は、入力された仮想拍間隔ｂｉ’及び観測時間を、テンポ重み算出部１０６１３に出力する。重み算出部１０６１は、仮想楽譜位置ｋｉ’及び仮想拍間隔ｂｉ’を、粒子重み算出部１０６１４に出力する。その後、ステップＳ２に進む。
【０１０１】
次に、本実施形態に係る楽譜位置推定装置１によって推定された楽譜位置の例について図を用いて説明する。図７−１４は、特に断らない限り粒子数Ｎが３００の場合における楽譜位置推定装置１の動作結果を示す。
図７は、本実施形態に係る楽譜位置の推定誤差（ｅｓｔｉｍａｔｉｏｎ  ｅｒｒｏｒ、ｐｒｅｄｉｃｔｉｏｎ  ｅｒｒｏｒ）の一例を示す図である。図７において、縦軸は楽譜位置の推定誤差（秒）を示し、横軸は時刻（秒）を示す。
図７によれば、時刻０−１０秒又は２０−３０秒付近では、推定誤差は概ね±３秒未満  である。時刻３５−４０秒付近での推定誤差が１５秒前後と、より大きくなる。さらに、時刻４１−５０秒でも推定誤差は概ね±３秒未満となる。
【０１０２】
図８は、本実施形態に係る信頼度係数の一例を示す図である。図８において、縦軸は信頼度係数を示し、横軸は時刻（秒）を示す。図８によれば、時刻０−１０秒又は２０-３０秒付近では、信頼度係数が概ね２０未満である。図７の結果を考慮すると、粒子が正しい楽譜位置に収束していることを示す。他方、信頼度係数が時刻３５−４０秒付近で約６０から１００へと増加する。この結果は、信頼度係数の増加が推定誤差の累積によって生じることを示す。また、図８によれば、時刻４１−５０秒付近で信頼度係数が低下する。この低下は、粒子が正しい楽譜位置に収束することを意味する。
従って、本実施形態では、信頼度算出部１０７が、信頼度差分が所定の閾値以上であることを検知することで、顕著な推定誤差をもたらす推定楽譜位置の出力を停止することができる。また、信頼度算出部１０７が、信頼度差分が所定の閾値以下であることを検知することで、推定誤差が小さい推定楽譜位置の出力を開始することができる。これにより、推定誤差が著しい推定楽譜位置の出力を回避できることが裏付けられる。
【０１０３】
図９は、本実施形態に係る推定誤差が生じる割合の一例を示す図である。図９において、縦軸は、推定誤差の絶対値が１秒未満となる割合（以下、低推定誤差率（ｒａｔｅ  ｏｆ  ｌｏｗ  ｐｒｅｄｉｃｔｉｏｎ  ｅｒｒｏｒ）と呼ぶ）であり、横軸は、テンポ（ｂｐｍ）である。図９において、○及び×は、演奏された楽曲のテンポ毎の低推定誤差率である。○は、ドラム音を含む楽曲を示す。×は、ドラム音を含まない楽曲を示す。
図９は、テンポが１１０ｂｐｍ以下の楽曲のほとんどでは、低推定誤差率は０．９よりも高く、楽譜位置が精度よく推定されていることを示す。他方、図９は、テンポが１２０ｂｐｍ以上の楽曲では、低推定誤差率は０．１よりも低く、楽譜位置が良好に推定されていないことを示す。その理由は、１２０ｂｐｍ以上のテンポでは、拍間隔ｂが０．５秒以下となるため、先行時間ΔＴ（この例では１秒）内にテンポ重み係数ｗｔ（式（９）参照）のピークが２個以上存在するからである。つまり、一意にテンポ重み係数ｗｔが定まらないことが推定誤差を増加させる要因と考えられる。
【０１０４】
また、図９は、テンポが１１０ｂｐｍ以下において、ドラム音を含む楽曲よりもドラム音を含まない楽曲のほうが、楽譜位置が精度よく推定されていることを示す。これは、ドラム音の非調和性（ｉｎｈａｒｍｏｎｉｃ  ｎａｔｕｒｅ）がクロマベクトル重み係数ｗｉｃｈ、周波数特性重み算出部１０６１２から入力された周波数特性重み係数ｗｉｓｐを低下させる有力な原因となるからである。
なお、上述の低推定誤差率とは異なり、０．５秒未満又は２秒未満の推定誤差をとる割合を低推定誤差率と定めた場合でも、本実施形態では図９と同様な結果が得られる。
【０１０５】
図１０は、本実施形態に係る推定誤差の一例を示す図である。図１０において、縦軸は、推定誤差であり、横軸は、テンポ（ｂｐｍ）である。図１０において、○及び×は、演奏された楽曲のテンポ毎の平均推定誤差である。○は、ドラム音を含む楽曲を示す。×は、ドラム音を含まない楽曲を示す。また、○及び×の上下に伸びる棒線は、推定誤差の標準偏差を示す。
【０１０６】
図１０は、テンポが１１０ｂｐｍ以下の楽曲では、平均推定誤差が−１０秒から０秒の範囲であり、テンポが１２０ｂｐｍ以上の楽曲では、平均推定誤差が１０秒から２０秒に分布している。また、図１０は、テンポが１１０ｂｐｍ以下の楽曲のほとんどでは、標準偏差は７秒以内にとどまるのに対し、テンポが１２０ｂｐｍ以上の楽曲では、標準偏差は１０秒前後となる。即ち、図１０は、テンポが１１０ｂｐｍ以下の楽曲では、テンポが１２０ｂｐｍ以上の楽曲よりも、楽譜位置が精度よく推定されることを示す。
また、図１０は、テンポが１１０ｂｐｍ以下において、ドラム音を含まない楽曲よりもドラム音を含む楽曲のほうが、標準偏差が大きいことを示す。この点において、ドラム音を含む楽曲よりもドラム音を含まない楽曲のほうが、楽譜位置が精度よく推定されることを裏付ける。
【０１０７】
図１１は、本実施形態に係る推定誤差が生じる割合の他の例を示す図である。
この例は、楽譜位置出力部１０８が推定楽譜位置を出力する場合（メロディレベル）のみを考慮して算出した低推定誤差率を示す。
図１１において、縦軸は、低推定誤差率であり、横軸は、テンポ（ｂｐｍ）である。図１１において、○及び×は、演奏された楽曲のテンポ毎の低推定誤差率である。○は、ドラム音を含む楽曲を示す。×は、ドラム音を含まない楽曲を示す。
【０１０８】
図１１は、テンポが１１０ｂｐｍ以下の楽曲のほとんどでは、低推定誤差率は０．９よりも高く、楽譜位置が精度よく推定されていることを示す。他方、図１１は、テンポが１２０ｂｐｍ以上の楽曲では、低推定誤差率は０．２から０．８の間を分布し、楽譜位置が良好に推定されていないことを示す。しかし、図１１が示す低推定誤差率は、図９が示す低推定誤差率よりも高い。この結果は、信頼度算出部１０７が算出する信頼度に基づいて、楽譜位置出力部１０８からの推定楽譜位置の出力の有無を制御する構成が、楽譜位置の推定精度を向上させることを示す。
【０１０９】
図１２は、本実施形態に係る推定誤差の他の例を示す図である。
この例も、楽譜位置出力部１０８が推定楽譜位置を出力する場合（メロディレベル）のみを考慮して算出した低推定誤差率を示す。
図１２において、縦軸は、推定誤差であり、横軸は、テンポ（ｂｐｍ）である。図１２において、○及び×は、演奏された楽曲のテンポ毎の平均推定誤差である。○は、ドラム音を含む楽曲を示す。×は、ドラム音を含まない楽曲を示す。また、○及び×の上下に伸びる棒線は、推定誤差の標準偏差を示す。
【０１１０】
図１２は、テンポが１１０ｂｐｍ以下の楽曲では、平均推定誤差が−２秒から２秒の範囲であり、テンポが１２０ｂｐｍ以上のほとんどの楽曲では、平均推定誤差が１０秒から１５秒に分布している。また、図１２は、テンポが１１０ｂｐｍ以下の楽曲のほとんどでは、標準偏差は１秒以内にとどまるのに対し、テンポが１２０ｂｐｍ以上の楽曲では、標準偏差は７〜１０秒となる。図１０に示す結果と比較すると、図１２に示す平均推定誤差、標準偏差ともに低くなる。
この結果も、信頼度算出部１０７が算出する信頼度に基づいて、楽譜位置出力部１０８からの推定楽譜位置の出力の有無を制御する構成が、楽譜位置の推定精度を向上させることを示す。
【０１１１】
図１３は、本実施形態に係る推定誤差と粒子数の関係の一例を示す図である。
図１３において、縦軸は推定誤差を、横軸は粒子数を示す。◇、□、及び△は、テンポが９８、１２２、及び２００ｂｐｍの楽曲について、粒子数毎の平均推定誤差を示す。◇、□、及び△の上下に伸びる棒線は、粒子数毎の標準偏差を示す。
【０１１２】
図１３によれば、テンポが９８ｂｐｍの楽曲について、平均推定誤差は、粒子数１５０について−５秒であるのに対し、粒子数３００以上について−２〜−３秒に収束する。標準偏差は、粒子数１５０について上下４秒であるのに対し、粒子数３００以上について上下１〜２秒に収束する。
テンポが１２２ｂｐｍの楽曲について、平均推定誤差は、粒子数１５０について１４秒であるのに対し、粒子数３００以上について３〜５秒に収束する。標準偏差は、粒子数１５０について±８秒であるのに対し、粒子数３００以上について±３秒に収束する。
これに対し、テンポが２００ｂｐｍの楽曲について、粒子数に関わらず平均推定誤差は１５秒前後である。標準偏差は粒子数に関わらず±８秒である。
即ち、図１３の結果は、テンポが約１２０ｂｐｍ以下の楽曲について、粒子数が３００あれば楽譜位置を精度よく推定できることを示す。反面、テンポが高い楽曲については楽譜位置を精度よく推定できないことを示す。
【０１１３】
図１４は、本実施形態に係る精度（ｐｒｅｃｉｓｉｏｎ）と再現率（ｒｅｃａｌｌ）の一例を示す図である。
図１４において、縦軸は精度及び再現率を示し、横軸はテンポ（ｂｐｍ）を示す。
ここで、精度ξｐはｘｔ／ｙと、再現率ξｒはｘｔ／ｚである。但し、ｘｔは、推定誤差ｅ（ｔ）の絶対値が１秒を越える場合であって、リズムレベル（楽譜位置出力部１０８から推定楽譜位置を出力する場合）における出力の個数である。ｙは、リズムレベルにおける全ての出力の個数である。ｚは、推定誤差ｅ（ｔ）の絶対値が１秒を超える場合の出力の個数である。
【０１１４】
図１４において、×は演奏された楽曲のテンポ毎の精度を、○は演奏された楽曲の再現率を示す。図１４によれば、ほとんどの楽曲（特にテンポが１２０ｂｐｍ以上の全ての楽曲）について精度が概ね１である。これは、信頼度算出部１０７において信頼度に基づいて推定楽譜位置の出力を停止する場合は、概ね推定誤差が著しいことを示す。図１４によれば、再現率の平均は０．４３である。ここで、テンポが１１０ｂｐｍ以下の楽曲について、２曲を除いて再現率は０．０５未満となり、テンポが１２０ｂｐｍ以上の楽曲について、再現率は０．２〜０．８の間で分布する。テンポが１１０ｂｐｍ以下の楽曲について、そもそも推定誤差ｅ（ｔ）の絶対値が１秒を超える場合が少ないため、全体として推定精度の低下させずに済む。
【０１１５】
上述の結果からも、本実施形態は、従来困難だった拍構造が不明確な楽曲（例えば、ドラム音を含まない楽曲）であっても演奏されている楽譜位置を精度よく推定できる。また、上述の結果において、テンポが速い楽曲について推定誤差が増加する傾向が認められた。その場合でも、先行時間ΔＴをより短くすれば、テンポが速い楽曲でも楽譜位置を精度よく推定できることが示唆される。
このように、本実施形態は、入力された音響信号の特徴量と楽譜情報の特徴量に基づいて重み係数を算出し、その重み係数に対応する仮想楽譜位置及び仮想テンポを用いて楽譜位置を推定する。そのため、本実施形態によれば演奏された音楽に対する楽譜位置を精度よく推定することが可能になる。
【０１１６】
（第２の実施形態）
以下、図面を参照しながら本発明の第２の実施形態について説明する。
図１５は、本実施形態に係る楽譜位置推定装置２の構成を示す概略図である。
楽譜位置推定位置２は、楽譜位置推定装置１の楽譜位置推定部１０６及び信頼度算出部１０７の代わりに楽譜位置推定部２０６及び信頼度算出部２０７を備える。
楽譜位置推定部２０６は、楽譜位置推定部１０６の重み算出部１０６１の代わりに重み算出部２０６１を備える。重み算出部２０６１は、クロマベクトル重み算出部１０６１１、周波数特性重み算出部１０６１２、及びテンポ重み算出部１０６１３を備える点で重み算出部１０６１と同様である。
また、楽譜位置推定装置２のその他の機能及び構成は第１の実施形態に係る楽譜位置推定装置１と同様である。以下、本実施形態について、第１の実施形態との差異点を主に説明する。
【０１１７】
信頼度算出部２０７は、例えば、上述の方式の何れかによって算出した信頼度係数を楽譜位置推定部２０６に出力する点が、信頼度算出部１０７と異なる。信頼度算出部２０７のその他の機能及び構成は信頼度算出部１０７と同様である。
なお、信頼度算出部２０７は、楽譜位置推定部２０６に出力する信頼度係数として、楽譜位置出力停止信号及び楽譜位置出力開始信号の出力の要否を判断するために用いる信頼度係数とは異なる方式で算出してもよい。例えば、信頼度算出部２０７は、後述するＤＴＷ（ｄｙｎａｍｉｃ  ｔｉｍｅ  ｗａｒｐｉｎｇ；動的時間収縮）法を用いて距離を算出し、算出した距離を信頼度係数として楽譜位置推定部２０６に出力してもよい。後述するように、信頼度算出部２０７は、このＤＴＷ法を用いて音響特徴量と楽譜特徴量に基づき距離を算出する過程で、音響特徴量の区間と合致する楽譜特徴量の区間を定めることができる。
本実施形態では、この距離を算出する構成を重み算出部２０６１に備え、この距離を算出してもよい。
【０１１８】
重み算出部２０６１は、信頼度算出部２０７から入力された信頼度係数に基づき観測時間を変更するか否か判断する。ここで、観測時間を変更する処理について説明する。
図１６は、本実施形態に係る観測時間を変更する処理の概要を示す概念図である。
図１６上段において、縦軸は振幅、横軸は時刻である、曲線は入力音響信号の波形の一例を表す。図１６上段左側は、楽譜位置推定装置１の重み算出部１０６１と同様に観測時間をシフトする処理（Ａ．シフト）を示す。実線の四角形は、現時刻の観測時間を示し、破線の四角形は次時刻の観測時間を示す。即ち、図１６上段左側は、観測時間長Ｌを変更せずに先行時間ΔＴだけ後の時刻に観測時間が変更することを示す。
図１６上段右側は、観測時間を延長させる処理（Ｂ．延長）を示す。実線の四角形は、現時刻の観測時間を示し、破線の四角形は次時刻の観測時間を示す。即ち、図１６上段右側は、観測開始時刻を固定し観測終了時刻を先行時間ΔＴだけ後の時刻に変更することにより、観測時間長がＬからＬ＋ΔＴに延長することを示す。
【０１１９】
図１６中段及び下段は、観測時間を現時刻までに入力された音響信号全体に延長する処理（Ｃ．全区間に延長）を示す。図１６中段及び下段ともに、横軸は時刻を示す。図１６中段において、縦軸は振幅を示し、曲線は入力音響信号の波形の一例を表す。図１６下段は、楽譜情報を五線上に音階とその長さを示す音符を用いて示す。
図１６中段右側の実線の四角形は、現時刻の観測時間を示し、図１６下段右側の実線の四角形は同図中段の実線の四角形と対応する区間の楽譜情報を示す。図１６中段の破線の四角形は、次時刻の観測時間を示し、図１６下段の破線の四角形は同図中段の破線の四角形と対応する区間の楽譜情報を示す。即ち、図１６中段及び下段は、観測時間長Ｌ＋ΔＴの観測時間を現時刻までに入力された音響信号全体Ｌａｌｌに延長することを示す。
【０１２０】
次に本実施形態に係る重み算出部２０６１が観測時間を変更する処理について説明する。
図１７は、本実施形態に係る観測時間を変更する処理を示すフローチャートである。
重み算出部２０６１は、以下に説明する処理を、図６のステップＳ７の代わりに行う。
（ステップＳ７１）重み算出部２０６１は、信頼度算出部２０７から信頼度係数が入力される。その後、ステップＳ７２に進む。
（ステップＳ７２）重み算出部２０６１は、入力された信頼度係数が予め設定した閾値よりも大きい（即ち、推定楽譜時刻の信頼性が低い）か否か判断する。重み算出部２０６１は、信頼度が予め設定した閾値よりも小さいと判断したとき（ステップＳ７２  Ｎ）、ステップＳ７３に進む。重み算出部２０６１は、信頼度係数が予め設定した閾値よりも大きい又は等しいと判断したとき（ステップＳ７２  Ｙ）、ステップＳ７４に進む。
（ステップＳ７３）重み算出部２０６１は、観測時間を、観測時間長Ｌを変更せずに先行時間ΔＴだけシフトさせる（Ａ．シフト）。その後、図４のステップＳ２に進む。
（ステップＳ７４）重み算出部２０６１は、誤り回数を１加算してカウントする。誤り回数とは、楽譜位置の推定が誤った回数であり、初期値はゼロである。その後、ステップＳ７５に進む。
【０１２１】
（ステップＳ７５）重み算出部２０６１は、誤り回数が予め設定した閾値（例えば、５）よりも大きいか否か判断する。重み算出部２０６１は、誤り回数が予め設定した閾値よりも大きいと判断したとき（ステップＳ７５  Ｙ）、ステップＳ７７に進む。重み算出部２０６１は、誤り回数が予め設定した閾値よりも小さい又は等しいと判断したとき（ステップＳ７５  Ｎ）、ステップＳ７６に進む。
（ステップＳ７６）重み算出部２０６１は、観測時間の開始時刻を変更せず、終了時刻を先行時間ΔＴだけ延長することで観測時間長を延長する（Ｂ．延長）。誤り回数が予め設定した閾値（例えば、５）よりも大きいか否か判断する。その後、図４のステップＳ２に進む。
（ステップＳ７７）重み算出部２０６１は、観測時間を現時刻までに入力された音響信号全体に延長する（Ｃ．全区間に延長）。その後、図４のステップＳ２に進む。
【０１２２】
これにより、楽譜位置の推定誤差が著しい場合でも、本実施形態が音響信号（又は音響特徴量）と楽譜情報（又は楽譜特徴量）の観測時間を延長して、両特徴量が合致する部分が見出される可能性が高くなる。これにより、本実施形態に係る楽譜位置推定装置２は、推定楽譜位置の精度が低下した状態から精度がより向上した状態に回復することができる。ひいては、本実施形態は、推定楽譜位置の信頼性が低く出力されない状態から、信頼性の低下を回避して推定楽譜位置の出力を再開することができる。
【０１２３】
次に、上述のＤＴＷ法を用いて信頼度係数を算出する処理について説明する。
ＤＴＷ法では、重み算出部２０６１又は信頼度算出部２０７（以下、重み算出部２０６１等と呼ぶ）は、距離の算出対象とする音響特徴量の一部の区間と楽譜特徴量の一部の区間をマッチングさせ、全体の距離（累積距離）を最小化する音響特徴量の区間と楽譜特徴量の区間の組を決定する。即ち、重み算出部２０６１等は、この最小化された累積距離を信頼度係数と定める。以下の例では、このような音響特徴量の区間の始点を音響始点、終点を音響終点と呼び、このような楽譜特徴量区間の始点を楽譜始点、終点を楽譜終点と呼ぶ。
【０１２４】
図１８は、本実施形態に係る楽譜位置推定装置２が行うＤＴＷ法を示すフローチャートである。
（ステップＳ７１１）重み算出部２０６１等は、楽譜終点ｎ２を、例えば楽譜特徴量の時系列データ冒頭のフレームを含むサブ系列に初期設定する。
楽譜終点ｎ２は、ここでは楽譜特徴量の時系列データにおいて累積距離Ｄｎ２を算出するための区間の終点、つまり、その区間の楽譜特徴量を含むサブ系列（ｓｕｂｓｅｑｕｅｎｃｅ）である。サブ系列とは、楽譜特徴量の時系列データ及び音響特徴量の時系列データについて両者をマッチングすることにより両者が合致する区間を探索する単位となる小区間である。
本実施形態では、観測時間が複数個のサブ系列からなり、１個のサブ系列に整数個の楽譜フレームが含まれるように、サブ系列の長さ（例えば、２楽譜フレーム、１／２４秒に相当）を予め定めておく。その後、ステップＳ７１２に進む。
【０１２５】
（ステップＳ７１２）重み算出部２０６１等は、楽譜終点ｎ２が楽譜特徴量の時系列データの最初のサブ系列に達するまで、ステップＳ７１３及びステップＳ７１４を実行する。その後、ステップＳ７１５に進む。
（ステップＳ７１３）重み算出部２０６１等は、楽譜特徴量と音響特徴量をマッチング処理する。重み算出部２０６１等は、例えば後述するマッチング処理を行って、両特徴量間における累積距離Ｄｎ２を算出し、楽譜始点を示す対象楽譜サブ系列ｎ１及び音響始点を示す対象音響サブ系列ｍ１を決定する。決定された楽譜特徴量と音響特徴量とのマッチング処理の詳細については後述する。その後、ステップＳ７１４に進む。
（ステップＳ７１４）重み算出部２０６１等は、楽譜終点ｎ２を直前のサブ系列ｎ２−１にシフトさせる。その後、ステップＳ７１２に進む。
（ステップＳ７１５）重み算出部２０６１等は、累積距離Ｄｎ２が最小となる楽譜始点ｎ１から楽譜終点ｎ２までの楽譜特徴量を、音響特徴量と合致する区間と決定する。また、重み算出部１０６１等は、累積距離Ｄｎ２が最小となる音響始点ｍ１から音響終点ｍ２までの音響特徴量を、楽譜特徴量と合致する音響特徴量の区間と決定する。その後、図４のステップＳ２に進む。
これにより、重み算出部２０６１は、楽譜特徴量と合致する音響特徴量の時系列の区間を観測時間として重み係数ｗｉを算出することになる。
【０１２６】
次に、楽譜特徴量と音響特徴量とのマッチング処理について説明する。
図１９は、本実施形態に係る楽譜特徴量と音響特徴量とのマッチング処理を示すフローチャートである。
（ステップＳ７２１）重み算出部２０６１等は、音響終点ｍ２を、例えば音響特徴量の時系列データにおける最後のフレームを含むサブ系列に初期設定する。重み算出部２０６１等は、対象楽譜サブ系列ｎ１を楽譜終点ｎ２と初期設定し、対象音響サブ系列ｍ１を音響終点ｍ２と初期設定する。ここで、対象楽譜サブ系列ｎ１とは、楽譜特徴量の時系列データのうち処理対象となるサブ系列を示す。対象音響サブ系列ｍ１とは、音響特徴量の時系列データのうち処理対象となるサブ系列を示す。その後、ステップＳ７２２に進む。
このように最後のサブ系列からマッチング処理を行うのは、楽譜位置推定装置２が最新の音響信号に基づいて楽譜位置を推定できるようにするためである。
【０１２７】
（ステップＳ７２２）重み算出部２０６１等は、累積距離Ｄｎ２の初期値を対象楽譜サブ系列ｎ１における楽譜特徴量と対象音響サブ系列ｍ１における音響特徴量との距離ｄｎ１，ｍ１と初期設定する。距離ｄｎ１，ｍ１は、楽譜特徴量と音響特徴量の類似度又は相関性を示す尺度である。距離ｄｎ１，ｍ１は、両特徴量の類似度又は相関性が高いほど小さい値をとる関数であればよい。例えば、距離ｄｎ１，ｍ１は、クロマベクトル重み係数ｗｃｈに基づくエントロピー−ｗｃｈｌｏｇｗｃｈ、周波数特性重み係数ｗｉｓｐに基づくエントロピー−ｗｓｐｌｏｇｗｓｐ、又はこれらの合計値でもよい。ここで、重み算出部２０６１等は、例えば、クロマベクトル重み係数ｗｃｈを算出する際、対象楽譜サブ系列ｎ１又は対象音響サブ系列ｍ１に相当する時間を観測時間として式（４）及び（５）を用いてもよい。また、重み算出部１０６１等は、例えば、周波数特性重み係数ｗｉｓｐを算出する際、対象楽譜サブ系列ｎ１又は対象音響サブ系列ｍ１に相当する時間を観測時間として式（５）−（８）を用いてもよい。その後、ステップＳ７２３に進む。
【０１２８】
（ステップＳ７２３）重み算出部２０６１等は、対象楽譜サブ系列ｎ１が楽譜特徴量の時系列データにおける最初のサブ系列に達する、又は対象音響サブ系列ｍ１が音響特徴量の時系列データにおける最初のサブ系列に達するまで、ステップＳ７２４〜ステップＳ７２６の処理を繰り返す。
（ステップＳ７２４）重み算出部２０６１等は、次の３通り（１）−（３）のサブ系列の組それぞれについて距離を算出する。重み算出部２０６１等が、３通りの距離を算出するのは、演奏される音楽におけるテンポのゆらぎが、楽譜特徴量を相対的に遅延させる場合及び音響特徴量を相対的に遅延させる場合を考慮するためである。
（１）直前の対象楽譜サブ系列ｎ１−１における楽譜特徴量と現在の対象音響サブ系列ｍ１における音響特徴量との間の距離ｄｎ１−１，ｍ１、
（２）直前の対象楽譜サブ系列ｎ１−１における楽譜特徴量と直前の対象音響サブ系列ｍ１−１における音響特徴量との間の距離ｄｎ１−１，ｍ１−１、及び
（３）現在の対象楽譜サブ系列ｎ１における楽譜特徴量と直前の対象音響サブ系列ｍ１−１における音響特徴量との間の距離ｄｎ１，ｍ１−１。
このステップにおいて、距離を算出する方法はステップＳ７２２と同様な方法でよい。
重み算出部１０６１等は、ｄｎ１−１，ｍ１、ｄｎ１−１，ｍ１−１、及びｄｎ１，ｍ１−１のうち、最小の距離を選択する。その後、ステップＳ７２５に進む。
【０１２９】
（ステップＳ７２５）重み算出部２０６１等は、累積距離Ｄｎ２に選択された距離を加算する。その後、ステップＳ７２６に進む。
（ステップＳ７２６）重み算出部２０６１等は、対象楽譜サブ系列ｎ１及び対象音響サブ系列ｍ１を選択された距離に対応するサブ系列の組に変更する。つまり、距離ｄｎ１−１，ｍ１、が選択された場合、重み算出部２０６１等は、対象楽譜サブ系列ｎ１を直前の対象楽譜サブ系列ｎ１−１に変更し、対象音響サブ系列ｍ１を変更しない。
距離ｄｎ１−１，ｍ１−１、が選択された場合、重み算出部２０６１等は、対象楽譜サブ系列ｎ１を直前の対象楽譜サブ系列ｎ１−１に変更し、対象音響サブ系列ｍ１を直前の対象音響サブ系列ｍ１−１に変更する。
距離ｄｎ１−１，ｍ１−１、が選択された場合、重み算出部２０６１等は、対象楽譜サブ系列ｎ１を変更せず、対象音響サブ系列ｍ１を直前の対象音響サブ系列ｍ１−１に変更する。
その後、ステップＳ７２３に進む。
【０１３０】
これにより、重み算出部２０６１等は、音響特徴量と楽譜特徴量が合致する区間を特定して、重み算出部２０６１は、特定された区間に基づき重み係数ｗｉを算出する。よって、本実施形態に係る楽譜位置推定装置２は、さらに楽譜位置の推定誤りから容易に回復し、推定楽譜位置の出力を再開することができる。
【０１３１】
なお、上述の観測時間を変更する処理において、観測時間を全区間に延長（Ｃ．）した後にＤＴＷ法を用いて楽譜特徴量と音響特徴量が合致する区間を特定することを前提にした説明したが、本実施形態ではこれには限られない。本実施形態では、シフト（Ａ．）した後や、先行時間ΔＴだけ延長（Ｂ．）した後の観測時間における楽譜特徴量と音響特徴量の時系列データにＤＴＷ法を用いて合致する区間を特定してもよい。これにより、逐次に楽譜位置の推定誤差が精度よく算出できる状態を回復することができる。
本実施形態では、観測時間を延長（Ｂ．）する際の延長時間は先行時間ΔＴに限られず、任意の時間でよいし、始点を固定した時刻にせず先行させてもよい。
なお、本実施形態では、より大きいサブ系列（マッチング処理を行う区間数がより少ない）に基づきＤＴＷ法を実行して楽譜特徴量を基準に合致する音響特徴量の区間の候補を決定しておき、その決定された区間の候補に対し、より小さいサブ系列（マッチング処理を行う区間数がより多い）に基づきＤＴＷ法を実行してもよい。
また、上述のＤＴＷ法においては、楽譜特徴量を基準に合致する音響特徴量の区間を決定していたが、本実施形態では、音響特徴量を基準に合致する楽譜特徴量の区間を決定してもよい。
【０１３２】
本実施形態において、周波数特性分析部１０２１は、入力されたディジタル音響信号が周波数領域信号に変換できれば、ＦＦＴ以外の方式、例えばＤＦＴ（Ｄｉｓｃｒｅｔｅ  Ｆｏｕｒｉｅｒ  Ｔｒａｎｓｆｏｒｍ、離散フーリエ変換）を用いてもよい。
本実施形態において、再標本化部１０６２は、最も規格化粒子重み係数ｐｉが小さいものから予め定めた数Ｎｒ（例えば５０個）の粒子を棄却してもよい。
【０１３３】
なお、上述した実施形態における楽譜位置推定装置１又は２の一部、例えば、音響特徴量生成部１０２、楽譜特徴量生成部１０５、楽譜位置推定部１０６、２０６、及び信頼度算出部１０７、２０７をコンピュータで実現するようにしても良い。その場合、この制御機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現しても良い。なお、ここでいう「コンピュータシステム」とは、楽譜位置推定装置１又は２に内蔵されたコンピュータシステムであって、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含んでも良い。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。
また、上述した実施形態における楽譜位置推定装置１又は２の一部、または全部を、ＬＳＩ（Ｌａｒｇｅ  Ｓｃａｌｅ  Ｉｎｔｅｇｒａｔｉｏｎ）等の集積回路として実現しても良い。楽譜位置推定装置１又は２の各機能ブロックは個別にプロセッサ化してもよいし、一部、または全部を集積してプロセッサ化しても良い。また、集積回路化の手法はＬＳＩに限らず専用回路、または汎用プロセッサで実現しても良い。また、半導体技術の進歩によりＬＳＩに代替する集積回路化の技術が出現した場合、当該技術による集積回路を用いても良い。
【０１３４】
以上、図面を参照してこの発明の一実施形態について詳しく説明してきたが、具体的な構成は上述のものに限られることはなく、この発明の要旨を逸脱しない範囲内において様々な設計変更等をすることが可能である。
【符号の説明】
【０１３５】
１、２…楽譜位置推定装置、１０１…音響信号入力部、１０２…音響特徴量生成部、
１０２１…周波数特性分析部、１０２２…クロマベクトル生成部、
１０３…楽譜情報記憶部、１０４…楽譜情報入力部、
１０５…楽譜特徴量生成部、１０５１…周波数特性合成部、
１０５２…クロマベクトル生成部、
１０６、２０６…楽譜位置推定部、１０６１、２０６１…重み算出部、
１０６１１…クロマベクトル重み算出部、１０６１２…周波数特性重み算出部、
１０６１３…テンポ重み算出部、１０６１４…粒子重み算出部、
１０６２…再標本化部、１０６３…状態遷移部、１０６４…楽譜位置算出部、
１０６５…テンポ算出部、
１０７、２０７…信頼度算出部、１０８…楽譜位置出力部、１０９…テンポ出力部
【図１】
2011180590.tif 000020
【図２】
2011180590.tif 000021
【図３】
2011180590.tif 000022
【図４】
2011180590.tif 000023
【図５】
2011180590.tif 000024
【図６】
2011180590.tif 000025
【図７】
2011180590.tif 000026
【図８】
2011180590.tif 000027
【図９】
2011180590.tif 000028
【図１０】
2011180590.tif 000029
【図１１】
2011180590.tif 000030
【図１２】
2011180590.tif 000031
【図１３】
2011180590.tif 000032
【図１４】
2011180590.tif 000033
【図１５】
2011180590.tif 000034
【図１６】
2011180590.tif 000035
【図１７】
2011180590.tif 000036
【図１８】
2011180590.tif 000037
【図１９】
2011180590.tif 000038