(57)【要約】【課題】  演奏音の基本周波数を正確に判定することができる演奏判定装置を提供する。【解決手段】  基本周波数推定処理１１５では、楽音の高調波構造を模した音モデル１１５Ｍを用い、収音部１０４により取得された演奏音の周波数成分が各種の基本周波数に対応した音モデルを重み付け加算した混合分布となるように各音モデルに対する重み値を最適化し、演奏音の基本周波数を推定する。類似度評価処理１１６は、この推定された基本周波数と楽曲データから得られた目標基本周波数との類似度を求め、一致判定処理１１７は、その類似度に基づいて演奏が正しく行われたか否かの判定を行う。【選択図】図１
2008209572.tif 000002
【特許請求の範囲】
【請求項１】
演奏音を外部から取得する演奏音取得手段と、
前記演奏音取得手段により取得される演奏音が本来有しているべき基本周波数である目標基本周波数を取得する目標基本周波数取得手段と、
楽器が発生する音の高調波構造を模した音モデルを用い、前記演奏音取得手段により取得された演奏音の周波数成分が各種の基本周波数に対応した音モデルを重み付け加算した混合分布となるように各音モデルに対する重み値を逐次更新して最適化し、前記演奏音取得手段により取得された演奏音の基本周波数を推定する基本周波数推定手段と、
前記目標周波数取得手段により取得された目標基本周波数と前記基本周波数推定手段における基本周波数の推定結果とに基づいて、前記演奏音取得手段により取得された演奏音の基本周波数に関する判定を行う判定手段と
を具備することを特徴とする演奏判定装置。
【請求項２】
前記基本周波数推定手段は、前記各種の基本周波数に対応した音モデルに対する重み値の最適化に際して、各音モデルに対する重み値に対し、前記目標基本周波数取得手段により取得された目標基本周波数における重み値を強調した事前分布を付与することを特徴とする請求項１に記載の演奏判定装置。
【請求項３】
前記基本周波数推定手段は、前記音モデルの高調波構造における各高調波成分の大きさの比率を変化させつつ、前記各音モデルに対する重み値を逐次更新して最適化する処理を行うことを特徴とする請求項１に記載の演奏判定装置。
【請求項４】
コンピュータを、
演奏音を外部から取得する演奏音取得手段と、
前記演奏音取得手段により取得される演奏音が本来有しているべき基本周波数である目標基本周波数を取得する目標基本周波数取得手段と、
楽器が発生する音の高調波構造を模した音モデルを用い、前記演奏音取得手段により取得された演奏音の周波数成分が各種の基本周波数に対応した音モデルを重み付け加算した混合分布となるように各音モデルに対する重み値を逐次更新して最適化し、前記演奏音取得手段により取得された演奏音の基本周波数を推定する基本周波数推定手段と、
前記目標周波数取得手段により取得された目標基本周波数と前記基本周波数推定手段における基本周波数の推定結果とに基づいて、前記演奏音取得手段により取得された演奏音の基本周波数に関する判定を行う判定手段として機能させることを特徴とするコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、楽譜等において指定された通りの音高の演奏音が発音されているか否かの判定を行う演奏判定装置および演奏判定プログラムに関する。
【背景技術】
【０００２】
従来から演奏の独習機能を備えた楽器が各種提供されている。鍵盤楽器を例に挙げると、この種の独習機能を備えた楽器は、表示装置への表示等により、押鍵すべき鍵をユーザに指示し、ユーザによって押鍵された鍵を検知して、正しい押鍵がなされているかをユーザに知らせ、ユーザに鍵盤演奏の独習を行わせる。このような独習機能を実現するためには、ユーザによって押鍵された鍵を検知する必要がある。従って、キースキャン機構を備えていない鍵盤楽器には、このような独習機能を持たせることができないという問題があった。そこで、鍵盤操作等を検知するのではなく、演奏音を収音して、その周波数解析を行い、楽譜通りの正しい音高の演奏音が発生しているか否かを判定する技術が提案されるに至った。例えば特許文献１に開示された技術では、各種の音高のピアノ音を収音して、それらの周波数解析を行い、各音高のピアノ音のパワースペクトルを求めて予め記憶しておく。そして、ピアノ演奏が行われる際には、演奏音を収音し、その周波数解析を行ってパワースペクトルを求め、この演奏音のパワースペクトルと予め記憶した各種の音高のピアノ音の各パワースペクトルとの類似度を求め、類似度の高さにより楽譜通りの演奏がなされているか否かを判定する。
【特許文献１】特開２００４−３４１０２６号公報
【特許文献２】特許第３４１３６３４号
【非特許文献１】後藤 真孝:"リアルタイム音楽情景記述システム: 全体構想と音高推定手法の拡張", 情報処理学会 音楽情報科学研究会 研究報告 2000-MUS-37-2,Vol.2000, No.94, pp.9-16, ２０００年１０月１６日
【発明の開示】
【発明が解決しようとする課題】
【０００３】
しかしながら、一般に楽音のパワースペクトルは多くの周波数位置に倍音成分を持っており、各成分の割合も多様であるため、対比される２つの楽音がある場合に、それらの基本周波数は異なっているがパワースペクトルの形状は類似している、ということが起こりうる。このため、上述した特許文献１の技術では、ある基本周波数の演奏音が収音された場合に、この演奏音とは基本周波数が異なるがパワースペクトルの形状が類似したピアノ音が選択され、演奏音の音高が誤って判定される可能性がある、という問題があった。また、特許文献１の技術では、収音した演奏音の基本周波数を求めないので、「本来、この音高であるべきものをこの音高で演奏した。」といった形式で演奏の誤りを指摘することができない、という問題があった。
【０００４】
この発明は、以上説明した事情に鑑みてなされたものであり、演奏音の基本周波数を正確に判定することができる演奏判定装置を提供することを目的としている。
【課題を解決するための手段】
【０００５】
この発明は、演奏音を外部から取得する演奏音取得手段と、前記演奏音取得手段により取得される演奏音が本来有しているべき基本周波数である目標基本周波数を取得する目標基本周波数取得手段と、楽器が発生する音の高調波構造を模した音モデルを用い、前記演奏音取得手段により取得された演奏音の周波数成分が各種の基本周波数に対応した音モデルを重み付け加算した混合分布となるように各音モデルに対する重み値を逐次更新して最適化し、前記演奏音取得手段により取得された演奏音の基本周波数を推定する基本周波数推定手段と、前記目標周波数取得手段により取得された目標基本周波数と前記基本周波数推定手段における基本周波数の推定結果とに基づいて、前記演奏音取得手段により取得された演奏音の基本周波数に関する判定を行う判定手段とを具備することを特徴とする演奏判定装置およびコンピュータを該演奏判定装置として機能させるコンピュータプログラムを提供する。
【０００６】
かかる発明によれば、楽器が発生する音の高調波構造を模した音モデルを用い、演奏音取得手段により取得された演奏音の周波数成分が各種の基本周波数に対応した音モデルを重み付け加算した混合分布となるように各音モデルに対する重み値を逐次更新して最適化し、演奏音取得手段により取得された演奏音の基本周波数を推定するようにしたので、演奏音の基本周波数を精度良く推定し、演奏音の基本周波数に関する判定を正確に行うことができる。
【発明を実施するための最良の形態】
【０００７】
以下、図面を参照し、この発明の実施の形態を説明する。
【０００８】
＜全体構成＞
図１は、この発明による演奏判定装置の一実施形態である教習伴奏システムの構成を示すブロック図である。この教習伴奏システムは、例えば鍵盤楽器等の楽器において稼動するシステムであり、ユーザに当該楽器の演奏の独習を行わせるシステムである。図１において、制御部１０１は、各種のプログラムを実行するＣＰＵと、このＣＰＵにワークエリアとして用いられるＲＡＭ等により構成されている。図１において、制御部１０１を表すボックス内には、制御部１０１のＣＰＵによって実行されるプログラムのうち、本実施形態による教習伴奏システムとしての機能を実現するためのプログラムの処理内容が図示されている。操作部１０２は、各種のコマンドや情報をユーザから受け取るための装置であり、楽器の本体に配備されたパネルスイッチ等の操作子群によって構成されている。ＨＤＤ（ハードディスク装置）１０３は、各種のプログラムやデータベースを記憶した記憶装置である。本実施形態による教習伴奏システムとしての機能を実現するためのプログラムもこのＨＤＤ１０３に記憶されている。制御部１０１のＣＰＵは、操作部１０２の操作により、教習伴奏システムとしての機能を起動すべき旨のコマンドが与えられたとき、この教習伴奏システムとしての機能を実現するためのプログラムをＲＡＭにロードして実行する。
【０００９】
収音部１０４は、外界の音を収音してアナログ音響信号を出力するマイクロフォンと、このアナログ音声信号をデジタル音響信号に変換するＡ／Ｄ変換器等により構成されている。本実施形態において、この収音部１０４は、演奏音を外部から取得する演奏音取得手段として用いられる。
【００１０】
楽曲記憶部１０５は、楽曲データを記憶するための記憶装置であり、例えばＲＡＭにより構成されている。ここで、楽曲データは、楽曲を構成するメロディパート、ベースパート等の各種のパートに対応した演奏データの集合体である。１つのパートに対応した演奏データは、演奏音の発生を指示するイベントデータと、演奏音の発生タイミングを指示するタイミングデータとにより構成された時系列データである。データ入力部１０６は、各種の楽曲の楽曲データを外部から取り込むための手段である。例えばＦＤやＩＣメモリ等の記憶媒体から楽曲データを読み出す装置やネットワークを介してサーバから楽曲データをダウンロードする通信装置等がこのデータ入力部１０６として用いられる。
【００１１】
サウンドシステム１０７は、デジタル音響信号をアナログ音響信号に変換するＤ／Ａ変換器と、アナログ音響信号を音として出力するスピーカ等により構成されている。表示部１０８は、例えば液晶パネルディスプレイであり、本実施形態では、演奏すべき楽曲の表示、押鍵すべき鍵をユーザに示すための鍵盤の画像の表示、ユーザにより行われた演奏が適切なものであったか否かの判定結果の表示等を行うための手段として用いられる。なお、判定結果は、表示に限らず、報知音や振動による提示などによりユーザに示してもよい。
【００１２】
次に、本実施形態における教習伴奏システムとしての機能を実現するためのプログラムの処理内容について説明する。まず、楽曲入力処理１１１は、操作部１０２を介して与えられるコマンドに従って、データ入力部１０６により楽曲データ１０５ａを取得し、楽曲記憶部１０５に格納する処理である。演奏位置制御処理１１２は、ユーザが演奏すべき位置を制御するとともに、その演奏位置に対応した演奏データを楽曲記憶部１０５内の楽曲データ１０５ａから抽出して出力し、かつ、この抽出した演奏データに基づき、ユーザが演奏すべき音の基本周波数である目標基本周波数を求めて出力する処理である。演奏位置制御処理１１２における演奏位置の制御の態様には２通りある。第１の態様は、ユーザがあるパートの楽器演奏をしており、この楽器演奏によりある演奏音が発生したとき、その演奏音が楽曲データ中の該当パートの演奏データにより指定された正しい音高を持った演奏音である場合に、その演奏音の次の演奏音の位置に演奏位置を進める態様である。第２の態様は、自動演奏、すなわち、パート毎にタイミングデータにより指定されたタイミングで各イベントデータを順次読み出し、この読み出しに連動させて演奏位置を進める態様である。演奏位置制御処理１１２にいずれの態様による演奏位置の制御を行わせるかは、操作部１０２を介して与えられるコマンドにより決定される。また、ユーザが楽曲データ１０５ａにおけるいずれのパートの演奏を行うかも、操作部１０２を介して与えられるコマンドにより決定される。
【００１３】
楽曲再生処理１１３は、演奏位置制御処理１１２により出力された演奏位置における演奏データの中からユーザの演奏パート以外のパートの演奏データを選択し、その演奏データが示す演奏音（すなわち、伴奏音）の波形のサンプルデータを生成してサウンドシステム１０７に供給する処理である。楽曲表示処理１１４は、ユーザが演奏を行うべき演奏位置および演奏音を示す情報を表示部１０８に表示させる処理である。この楽曲表示処理１１４の態様にも各種のものがあり得る。ある態様において、楽曲表示処理１１４では、楽曲データ１０５ａに従って、演奏すべき楽曲の楽譜を表示部１０８に表示させるとともに、演奏位置における演奏データに基づいて、ユーザが演奏すべき演奏位置を示すマークをその楽譜上に表示させる。他の態様における楽曲表示処理１１４では、例えば鍵盤の画像を表示部１０８に表示させるとともに、演奏位置における演奏データに基づいて、ユーザが押鍵すべき鍵を示すマークを表示させる。
【００１４】
基本周波数推定処理１１５は、楽器が発生する音の高調波構造を模した音モデル１１５Ｍを用い、収音部１０４によって収音された演奏音の周波数成分が各種の基本周波数に対応した音モデル１１５Ｍを重み付け加算した混合分布となるように各音モデル１１５Ｍに対する重み値を最適化し、最適化された各音モデルに対する重み値に基づいて、収音部１０４によって収音された演奏音の基本周波数を推定する処理である。本実施形態における基本周波数推定処理１１５では、演奏位置制御処理１１２から出力される目標基本周波数を事前知識として使用し、基本周波数の推定を行う。類似度評価処理１１６は、基本周波数推定処理１１５により推定された基本周波数と、演奏位置制御処理１１２により得られた目標基本周波数との類似度を算出する処理である。一致判定処理１１７は、類似度評価処理１１６により得られた類似度に基づき、基本周波数推定処理１１５により推定された基本周波数と演奏位置制御処理１１２により得られた目標基本周波数とが一致しているといえるか否かを判定する処理である。一致判定処理１１７の判定結果は、判定結果表示処理１１８に引き渡される他、上述した演奏位置制御処理１１２に引き渡される。演奏位置制御処理１１２では、上述した第１の態様が操作部１０２の操作により選択されている場合、一致判定処理１１７の判定結果が肯定的である場合に限り、演奏位置を次の演奏音の位置に進める制御が行われる。判定結果表示処理１１８は、一致判定処理１１７の判定結果、すなわち、ユーザが演奏データ通りの音高で演奏音を発生させたか否かを表示部１０８に表示させる処理である。
【００１５】
＜基本周波数推定処理１１５の内容＞
次に本実施形態における基本周波数推定処理１１５の処理内容の詳細について説明する。この基本周波数推定処理１１５は、特許文献２に開示された技術を基盤とし、これに非特許文献１に開示された改良を加えたものである。
【００１６】
特許文献２の技術では、メロディ音のものと考えられる帯域に属する周波数成分と、ベース音のものと考えられる帯域に属する周波数成分とを入力音響信号からＢＰＦにより別々に取り出し、それらの各帯域の周波数成分に基づき、メロディ音およびベース音の各々の基本周波数の推定を行う。
【００１７】
さらに詳述すると、特許文献２の技術では、音の高調波構造に対応した確率分布を持った音モデルを用意し、メロディ音の帯域の各周波数成分、ベース音の帯域の各周波数成分が、様々な基本周波数に対応した各音モデルを重み付け加算した混合分布であると考える。そして、各音モデルの重みの値をＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）アルゴリズムを用いて推定する。
【００１８】
このＥＭアルゴリズムは、隠れ変数を含む確率モデルに対して最尤推定を行うための反復アルゴリズムであり、局所最適解を求めることができる。ここで、最も大きな重みの値を持つ確率分布は、その時点で最も優勢な高調波構造であるとみなすことができるため、あとはその優勢な高調波構造における基本周波数を音高として求めればよい。この手法は基本周波数成分の存在に依存しないため、ミッシングファンダメンタル現象も適切に扱うことができ、基本周波数成分の存在に依存せずに、最も優勢な高調波構造を求めることができる。
【００１９】
非特許文献１は、特許文献２の技術に対して次の拡張を行ったものである。
＜拡張１：音モデルの多重化＞
特許文献２の技術では、同一基本周波数には１つの音モデルしか用意されていなかったが、実際には、ある基本周波数に、異なる高調波構造を持つ音が入れ替わり立ち替わり現れることがある。そこで、同一基本周波数に対して複数の音モデルを用意し、入力音響信号をそれらの混合分布としてモデル化した。
＜拡張２：音モデルのパラメータの推定＞
特許文献２の技術では、音モデルにおいて各高調波成分の大きさの比を固定していた（ある理想的な音モデルを仮定していた）。これは実世界の混合音中の高調波構造とは必ずしも一致しておらず、精度向上のためには洗練される余地が残されていた。そこで、音モデルの高調波成分の比率もモデルパラメータに加え、各時刻においてＥＭアルゴリズムにより推定するようにした。
＜拡張３：モデルパラメータに関する事前分布の導入＞
特許文献２の技術では、音モデルの重み（基本周波数の確率密度関数）に関する事前知識は仮定していなかった。しかし、この基本周波数の推定技術の用途によっては、たとえ事前に基本周波数がどの周波数の近傍にあるかを与えてでも、より誤検出の少ない基本周波数を求めたいというような要求も発生し得る。例えば、演奏分析やビブラート分析等の目的では、楽曲をヘッドホン聴取しながらの歌唱や楽器演奏によって、各時刻におけるおおよその基本周波数を事前知識として用意しておき、実際の楽曲中のより正確な基本周波数を得ることが求められている。そこで、特許文献２におけるモデルパラメータ（音モデルに対する重み値）の最尤推定の枠組みを拡張し、モデルパラメータに関する事前分布に基づいて最大事後確率推定（ＭＡＰ推定；Maximum A Posteriori Probability Estimation）を行うようにした。その際、＜拡張２＞においてモデルパラメータに加えた音モデルの高調波成分の大きさの比率に関する事前分布も導入した。
【００２０】
図２は、特許文献２の技術と非特許文献１の技術を組み合わせることにより構成された本実施形態における基本周波数推定処理１１５の処理内容を示す図である。基本周波数推定処理１１５ではメロディラインとベースラインを推定する。メロディは他よりも際立って聞こえる単音の系列、ベースはアンサンブル中で最も低い単音の系列であり、その時間的な変化の軌跡をそれぞれメロディラインＤｍ（ｔ）、ベースラインＤｂ（ｔ）と呼ぶ。時刻tにおける基本周波数Ｆ０をＦｉ（ｔ）（ｉ＝ｍ，ｂ）、振幅をＡｉ（ｔ）とすると、これらは以下のように表される。
【数１】
2008209572.tif 000003
【数２】
2008209572.tif 000004
【００２１】
このメロディラインＤｍ（ｔ）およびベースラインＤｂ（ｔ）を収音部１０４により収音された演奏音を示す入力音響信号から得るための手段として、基本周波数推定処理１１５は、瞬時周波数の算出１、周波数成分の候補の抽出２、周波数帯域の制限３、メロディラインの推定４ａおよびベースラインの推定４ｂの各処理を含む。また、メロディラインの推定４ａおよびベースラインの推定４ｂの各処理は、基本周波数の確率密度関数の推定４１およびマルチエージェントモデルによる基本周波数の継時的な追跡４２の各処理を各々含む。本実施形態では、ユーザの演奏パートがメロディパートである場合にはメロディラインの推定４ａが実行され、ユーザの演奏パートがベースパートである場合にはベースラインの推定４ｂが実行されるようになっている。
【００２２】
＜＜瞬時周波数の算出１＞＞
この処理では、入力音響信号を複数のＢＰＦからなるフィルタバンクに与え、フィルタバンクの各ＢＰＦの出力信号について、位相の時間微分である瞬時周波数（Flanagan, J.L. and Golden, R.M.: Phase Vocoder, The BellSystem
Technical J., Vol.45, pp.1493-1509 (1966)参照）を計算する。ここでは、上記Flanaganの手法を用い、短時間フーリエ変換(STFT)の出力をフィルタバンク出力と解釈して、効率良く瞬時周波数を計算する。入力音響信号ｘ(ｔ)に対する窓関数ｈ(ｔ)を用いたＳＴＦＴが式（３）および（４）により与えられるとき、瞬時周波数λ（ω，ｔ）は式（５）により求めることができる。
【００２３】
【数３】
2008209572.tif 000005
【数４】
2008209572.tif 000006
【数５】
2008209572.tif 000007
【００２４】
ここで、ｈ(ｔ)は時間周波数の局所化を与える窓関数である（例えば、最適な時間周波数の局所化を与えるガウス関数に２階のカーディナルＢ−スプライン関数を畳み込んで作成した時間窓など)。
【００２５】
この瞬時周波数を計算するのに、ウェーブレット変換を用いても良い。ここでは、計算量を減らすためにＳＴＦＴを用いるが、単一のＳＴＦＴのみを用いたのでは、ある周波数帯域における時間分解能や周波数分解能が悪くなってしまう。そこで、マルチレートフィルタバンク（Vetterli, M.: A Theory of Multirate Filter Banks, IEEE Trans. on
ASSP, Vol.ASSP-35, No.3, pp. 356-372 (1987)、参照）を構成し、リアルタイムに実行可能という制約のもとで、ある程度妥当な時間周波数分解能を得る。
【００２６】
＜＜周波数成分の候補の抽出２＞＞
この処理では、フィルタの中心周波数からその瞬時周波数への写像に基づいて、周波数成分の候補を抽出する（Charpentier, F.J.: Pitch detection using the short-termphase
spectrum, Proc. of ICASSP 86, pp.113-116 (1986)参照）。あるＳＴＦＴフィルタの中心周波数ωからその出力の瞬時周波数λ（ω，ｔ）への写像を考える。すると、もし周波数ψの周波数成分があるときには、ψがこの写像の不動点に位置し、その周辺の瞬時周波数の値はほぼ一定となる。つまり、全周波数成分の瞬時周波数Ψf(t)は、次式によって抽出することができる。
【００２７】
【数６】
2008209572.tif 000008
これらの周波数成分のパワーは、Ψf(t)の各周波数におけるＳＴＦＴパワースペクトルの値として得られるため、周波数成分のパワー分布関数Ψp(t)(ω)を次のように定義することができる。
【数７】
2008209572.tif 000009
【００２８】
＜＜周波数帯域の制限３＞＞
この処理では、抽出した周波数成分に重み付けすることで、周波数帯域を制限する。ここでは、メロディラインとベースライン用に、２種類のＢＰＦを用意する。メロディライン用ＢＰＦは、典型的なメロディラインの主要な基本波成分および高調波成分の多くを通過させることができ、かつ、基本周波数付近の重複が頻繁に起きる周波数帯域をある程度遮断する。一方、ベースライン用ＢＰＦは、典型的なベースラインの主要な基本周波数成分および高調波成分の多くを通過させることができ、かつ、他の演奏パートがベースラインよりも優勢になるような周波数帯域をある程度遮断する。
【００２９】
本実施形態では、以下、対数スケールの周波数をcentの単位(本来は音高差(音程)を表す尺度)で表し、Ｈｚで表された周波数ｆＨｚを、次のようにcentで表された周波数ｆｃｅｎｔに変換する。 
【数８】
2008209572.tif 000010
【数９】
2008209572.tif 000011
平均律の半音は１００ｃｅｎｔに、１オクターブは１２００ｃｅｎｔに相当する。
【００３０】
周波数ｘ ｃｅｎｔでのＢＰＦの周波数応答をＢＰＦｉ（ｘ）（ｉ＝ｍ，ｂ）とし、周波数成分のパワー分布関数をΨ’ｐ（ｔ）（ｘ）とすると、ＢＰＦを通過した周波数成分はＢＰＦｉ（ｘ）Ψ’ｐ（ｔ）（ｘ）と表すことができる。ただし、Ψ’ｐ（ｔ）（ｘ）は、周波数軸がｃｅｎｔで表されていることを除けばΨｐ（ｔ）（ω）と同じ関数である。ここで、次の段階の準備として、ＢＰＦを通過した周波数成分の確率密度関数ｐΨ（ｔ）（ｘ）を定義する。
【数１０】
2008209572.tif 000012
【００３１】
ここで、Ｐｏｗ（ｔ）は次式に示すようにＢＰＦを通過した周波数成分のパワーの合計である。
【数１１】
2008209572.tif 000013
【００３２】
＜＜基本周波数の確率密度関数の推定４１＞＞
この基本周波数の確率密度関数の推定４１では、ＢＰＦを通過した周波数成分の候補に対し、各高調波構造が相対的にどれくらい優勢かを表す基本周波数の確率密度関数を求める。この基本周波数の確率密度関数の推定４１の処理内容を、非特許文献１に開示されたものに改良を加えた内容となっている。
【００３３】
基本周波数の確率密度関数の推定４１では、上述した「拡張１」と「拡張２」を実現するために、同一基本周波数に対してＭｉ種類の音モデルがあるものとし（ｉはメロディ用（ｉ＝ｍ）かベース用（ｉ＝ｂ）かを示す）、基本周波数がＦであり、音モデルの種類がｍ番目の種類であり、モデルパラメータμ（ｔ）（Ｆ，ｍ）を持った音モデルｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））を次のように定義する。
【数１２】
2008209572.tif 000014
【数１３】
2008209572.tif 000015
【数１４】
2008209572.tif 000016
【数１５】
2008209572.tif 000017
【００３４】
この音モデルは、基本周波数がＦのときに、その高調波成分がどの周波数にどれくらい現れるかをモデル化したものである。Ｈｉは基本周波数成分も含めた高調波成分の数、Ｗｉ２はガウス分布Ｇ（ｘ；ｘ０，σ）の分散を表す。ｃ（ｔ）（ｈ｜Ｆ，ｍ）は、基本周波数がＦであるｍ番目の音モデルの第ｈ次調波成分の大きさを表し、次式を満たす。
【数１６】
2008209572.tif 000018
【００３５】
上記式（１６）に示すように、基本周波数Ｆの音モデルにおける重みｃ（ｔ）（ｈ｜Ｆ，ｍ）は、総和が１となるような予め定義された重みが用いられる。
【００３６】
基本周波数の確率密度関数の推定４１では、以上のような音モデルを使用し、周波数成分の確率密度関数ｐΨ（ｔ）（ｘ）が、次式で定義されるようなｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））の混合分布モデルｐ（ｘ｜θ（ｔ））から生成されたと考える。
【数１７】
2008209572.tif 000019
【数１８】
2008209572.tif 000020
【数１９】
2008209572.tif 000021
【数２０】
2008209572.tif 000022
ここで、ＦｈｉとＦｌｉは、許容される基本周波数の上限と下限であり、ｗ（ｔ）（Ｆ，ｍ）は次式を満たすような音モデルの重みである。
【数２１】
2008209572.tif 000023
【００３７】
実世界の混合音に対して事前に音源数を仮定することは不可能なため、式（１７）のように、あらゆる基本周波数の可能性を同時に考慮してモデル化することが重要となる。最終的に、モデルｐ（ｘ｜θ（ｔ））から、観測した確率密度関数ｐΨ（ｔ）（ｘ）が生成されたかのようにモデルパラメータθ（ｔ）を推定できれば、その重みｗ（ｔ）（Ｆ，ｍ）は各高調波構造が相対的にどれくらい優勢かを表すため、次式のように基本周波数の確率密度関数ｐＦ０（ｔ）（Ｆ）と解釈することができる。
【数２２】
2008209572.tif 000024
【００３８】
次に、上述した「拡張３」を実現するために、θ（ｔ）の事前分布ｐ０ｉ（θ（ｔ））を、式（２３）のように式（２４）と式（２５）の積で与える。
【数２３】
2008209572.tif 000025
【数２４】
2008209572.tif 000026
【数２５】
2008209572.tif 000027
【００３９】
ここで、ｐ０ｉ（ｗ（ｔ））とｐ０ｉ（μ（ｔ））は、最も起こりやすいパラメータをｗ０ｉ（ｔ）（Ｆ，ｍ）とμ０ｉ（ｔ）（Ｆ，ｍ）としたときに、そこで最大値を取るような単峰性の事前分布である。ただし、Ｚｗ、Ｚμは正規化係数、βｗｉ（ｔ）、βμｉ（ｔ）（Ｆ，ｍ）は、最大値をどれくらい重視した事前分布とするかを決めるパラメータで、０のときに無情報事前分布(一様分布）となる。また、Ｄｗ（ｗ０ｉ（ｔ）；ｗ（ｔ））、Ｄμ（μ０ｉ（ｔ）（Ｆ，ｍ）；μ（ｔ）（Ｆ，ｍ））は、次のようなＫ−Ｌ情報量（Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ’ｓ  ｉｎｆｏｒｍａｔｉｏｎ）である。
【数２６】
2008209572.tif 000028
【数２７】
2008209572.tif 000029
【００４０】
以上から、確率密度関数ｐΨ（ｔ）（ｘ）を観測したときに、そのモデルｐ（ｘ｜θ（ｔ））のパラメータθ（ｔ）を、事前分布ｐ０ｉ（θ（ｔ））に基づいて推定する問題を解けばよいことがわかる。この事前分布に基づくθ（ｔ）の最大事後確率推定量（ＭＡＰ推定量）は、次式を最大化することで得られる。
【数２８】
2008209572.tif 000030
【００４１】
この最大化問題は解析的に解くことが困難なため、前述のＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）アルゴリズムを用いてθ（ｔ）を推定する。ＥＭアルゴリズムは、Ｅステップ（ｅｘｐｅｃｔａｔｉｏｎ  ｓｔｅｐ）とＭステップ（ｍａｘｉｍｉｚａｔｉｏｎ  ｓｔｅｐ）を交互に繰返し適用することで、不完全な観測データ（この場合、ｐΨ（ｔ）（ｘ））から最尤推定をおこなうための反復アルゴリズムである。本実施形態では、ＥＭアルゴリズムを繰り返すことにより、ＢＰＦを通過した周波数成分の確率密度関数ｐΨ（ｔ）（ｘ）を、各種の基本周波数Ｆに対応した複数の音モデルｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））を重み付け加算した混合分布と考える場合において、最も尤もらしい重みのパラメータθ（ｔ）（＝｛ｗ（ｔ）（Ｆ，ｍ），μ（ｔ）（Ｆ，ｍ）｝）を求める。ここで、ＥＭアルゴリズムの各繰り返しでは、パラメータθ（ｔ）（＝｛ｗ（ｔ）（Ｆ，ｍ），μ（ｔ）（Ｆ，ｍ））に関して、古いパラメータ推定値θold（ｔ）（＝｛ｗｏｌｄ（ｔ）（Ｆ，ｍ），μｏｌｄ（ｔ）（Ｆ，ｍ）｝）を更新して新しい（より尤もらしい）パラメータ推定値θｎｅｗ（ｔ）（＝｛ｗｎｅｗ（ｔ）（Ｆ，ｍ），μｎｅｗ（ｔ）（Ｆ，ｍ）｝）を求めていく。θｏｌｄ（ｔ）の初期値には、１つ前の時刻ｔ−１における最終的な推定値を用いる。この古いパラメータ推定値θold（ｔ）から新しいパラメータ推定値θｎｅｗ（ｔ）を求める漸化式は、次のようになる。なお、この漸化式の導出過程は非特許文献１に詳細に説明されているので、そちらを参照されたい。
【数２９】
2008209572.tif 000031
【数３０】
2008209572.tif 000032
【００４２】
上記式（２９）および（３０）におけるｗＭＬ（ｔ）（Ｆ，ｍ）とｃＭＬ（ｔ）（ｈ｜Ｆ，ｍ）は、βｗｉ（ｔ）＝０、βμｉ（ｔ）（Ｆ，ｍ）＝０の無情報事前分布のとき、つまり、最尤推定の場合の推定値であり、次式により与えられる。
【数３１】
2008209572.tif 000033
【数３２】
2008209572.tif 000034
【００４３】
これらの反復計算により、事前分布を考慮した基本周波数の確率密度関数ｐＦ０（ｔ）（Ｆ）が、式（２３）によってｗ（ｔ）（Ｆ，ｍ）から求まる。さらに、すべての音モデルｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））の各高調波成分の大きさの比率ｃ（ｔ）（ｈ｜Ｆ，ｍ）も求まり、「拡張１」〜「拡張３」が実現される。
【００４４】
最も優勢な基本周波数Ｆｉ（ｔ）を決定するには、次式に示すように、基本周波数の確率密度関数ｐＦ０（ｔ）（Ｆ）（式（２２）より、式（２９）〜（３２）を反復計算した最終的な推定値として得られる）を最大にする周波数を求めればよい。
【数３３】
2008209572.tif 000035
こうして得られた周波数を音高とする。
【００４５】
＜＜マルチエージェントモデルによる基本周波数の継時的な追跡４２＞＞
基本周波数の確率密度関数において、同時に鳴っている音の基本周波数に対応する複数のピークが拮抗すると、それらのピークが確率密度関数の最大値として次々に選ばれてしまうことがあるため、このように単純に求めた結果は安定しないことがある。そこで、本実施形態では、大局的な観点から基本周波数を推定するために、基本周波数の確率密度関数の時間変化において複数のピークの軌跡を継時的に追跡し、その中で最も優勢で安定した基本周波数の軌跡を選択する。このような追跡処理を動的で柔軟に制御するために、マルチエージェントモデルを導入する。
【００４６】
マルチエージェントモデルは、１つの特徴検出器と複数のエージェントにより構成される（図３参照）。特徴検出器は、基本周波数の確率密度関数の中で目立つピークを拾い上げる。エージェントは基本的に、それらのピークに駆動されて軌跡を追跡していく。つまり、マルチエージェントモデルは、入力中で目立つ特徴を時間的に追跡する汎用の枠組みである。具体的には、各時刻において以下の処理がおこなわれる。
【００４７】
（１）基本周波数の確率密度関数が求まった後、特徴検出器は目立つピーク（最大ピークに応じて動的に変化する閾値を越えたピーク）を複数検出する。そして、目立つピークのそれぞれについて、周波数成分のパワーの合計Ｐｏｗ（ｔ）も考慮しながら、どれくらい将来有望なピークかを評価する。これは、現在時刻を数フレーム先の時刻とみなして、ピークの軌跡をその時刻まで先読みして追跡することで実現する。
【００４８】
（２）既に生成されたエージェントがあるときは、それらが相互作用しながら、目立つピークをそれに近い軌跡を持つエージェントへと排他的に割り当てる。複数のエージェントが割り当て候補に上がる場合には、最も信頼度の高いエージェントへと割り当てる。
【００４９】
（３）最も有望で目立つピークがまだ割り当てられていないときは、そのピークを追跡する新たなエージェントを生成する。
【００５０】
（４）各エージェントは累積ペナルティを持っており、それが一定の閾値を越えると消滅する。
【００５１】
（５）目立つピークが割り当てられなかったエージェントは、一定のペナルティを受け、基本周波数の確率密度関数の中から自分の追跡する次のピークを直接見つけようとする。もしそのピークも見つからないときは、さらにペナルティを受ける。さもなければ、ペナルティはリセットされる。
【００５２】
（６）各エージェントは、今割り当てられたピークがどれくらい有望で目立つかを表す度合いと、１つ前の時刻の信頼度との重み付き和によって、信頼度を自己評価する。
【００５３】
（７）時刻ｔにおける基本周波数Ｆｉ（ｔ）は、信頼度が高く、追跡しているピークの軌跡に沿ったパワーの合計が大きいエージェントに基づいて決定する。振幅Ａｉ（ｔ）は、基本周波数Ｆｉ（ｔ）の高調波成分等をΨｐ（ｔ）（ω）から抽出して決定する。
以上が本実施形態における基本周波数推定処理１１５の詳細である。
【００５４】
＜本実施形態の動作＞
次に本実施形態の動作を説明する。本実施形態において演奏位置制御処理１１２では、楽曲上においてユーザが演奏すべき位置を常時監視し、その演奏位置に対応した演奏データを楽曲記憶部１０５内の楽曲データ１０５ａから抽出して出力し、楽曲再生処理１１３および楽曲表示処理１１４に引き渡す。また、演奏位置制御処理１１２では、演奏位置に対応した演奏データに基づいてユーザの演奏パートの演奏音の目標基本周波数を求め、基本周波数推定処理１１５に引き渡す。
【００５５】
楽曲再生処理１１３では、演奏位置制御処理１１２から引き渡された演奏データに基づいて、ユーザの演奏パート以外のパートの演奏音（すなわち、伴奏音）の音響信号を生成し、サウンドシステム１０７にその再生を行わせる。また、楽曲表示処理１１４では、演奏位置制御処理１１２から引き渡された演奏データに基づいて、ユーザが演奏すべき演奏音を示す画像（例えば鍵盤において押鍵すべき鍵を示す画像）またはユーザが演奏すべき演奏位置を示す画像（楽譜上における演奏位置を示す画像）を表示部１０８に表示させる。
【００５６】
ユーザが楽器演奏を行って、その演奏音が収音部１０４により収音されると、その演奏音を示す入力音響信号が基本周波数推定処理１１５に引き渡される。基本周波数推定処理１１５では、楽器が発生する音の高調波構造を模した音モデル１１５Ｍを用い、入力音響信号の周波数成分が各種の基本周波数に対応した音モデル１１５Ｍを重み付け加算した混合分布となるように各音モデル１１５Ｍに対する重み値を最適化し、最適化された各音モデルに対する重み値に基づいて、入力音響信号が示す１または複数の演奏音の基本周波数を推定する。その際、本実施形態における基本周波数推定処理１１５では、演奏位置制御処理１１２から引き渡された目標基本周波数における重みが強調された事前分布ｐ０ｉ（θ（ｔ））を作成し、この事前分布ｐ０ｉ（θ（ｔ））用い、かつ、音モデルにおける各高調波成分の大きさの比を変化させつつＥＭアルゴリズムを実行し、演奏音の基本周波数を推定する。
【００５７】
そして、類似度評価処理１１６では、このようにして基本周波数推定処理１１５により推定された基本周波数と、演奏位置制御処理１１２により得られた目標基本周波数との類似度を算出する。類似度として何を用いるかについては各種の態様が考えられるが、例えば基本周波数推定処理１１５により推定された基本周波数の目標基本周波数に対する比（すなわち、対数表現された両周波数のずれであるセント値）を所定の値（例えば一音階分のセント値）によって割り算したものを類似度としてもよい。一致判定処理１１７では、この類似度評価処理１１６により得られた類似度に基づき、基本周波数推定処理１１５により推定された基本周波数と演奏位置制御処理１１２により得られた目標基本周波数とが一致しているといえるか否かを判定する。判定結果表示処理１１８では、一致判定処理１１７の判定結果、すなわち、ユーザが演奏データ通りの音高で演奏音を発生させたか否かを表示部１０８に表示させる。好ましい態様では、表示部１０８に楽譜を表示させ、判定結果表示処理１１８では、ユーザに演奏の誤りを適確に知らせるため、その楽譜上において、演奏位置に対応した演奏データにより指定された演奏音の音符（すなわち、目標基本周波数を示す音符）と、実際にユーザが発生させた演奏音の基本周波数を示す音符とを例えば異なる色で表示させる。
本実施形態では、演奏位置を進めつつ、以上の処理が繰り返される。
【００５８】
以上説明したように、本実施形態では、楽器が発生する音の高調波構造を模した音モデルを用い、収音部１０４により取得された演奏音の周波数成分が各種の基本周波数に対応した音モデルを重み付け加算した混合分布となるように各音モデルに対する重み値を最適化し、最適化された各音モデルに対する重み値に基づいて、演奏音の基本周波数を推定するようにしたので、演奏音の基本周波数を精度良く推定し、演奏音の基本周波数に関する判定を正確に行うことができる。また、本実施形態では、ユーザが発生させた演奏音の基本周波数を求めるので、「本来、この音高であるべきものをこの音高で演奏した。」といった形式で演奏の誤りをユーザに示すことができる。また、本実施形態では、演奏データから得られる目標基本周波数の位置において重みが強調された事前分布を作成し、この事前分布を用いて、基本周波数の推定を行うようにしているので、基本周波数の推定の精度を向上させることができる。また、本実施形態によれば、音モデルにおける高調波成分の大きさの比率を変えつつＥＭアルゴリズムを実行して演奏音の基本周波数の推定を行うようにしている。従って、ユーザが発生させる演奏音のスペクトル形状が演奏の強弱、タッチにより大きく変動する状況においても、そのスペクトル形状の変化に合わせて音モデルの高調波成分の大きさの比率を変化させ、精度よく演奏音の基本周波数の推定を行うことができる。
【００５９】
＜他の実施形態＞
以上、この発明の一実施形態について説明したが、この発明には他にも実施形態があり得る。例えば次の通りである。
【００６０】
（１）上記実施形態では、基本周波数推定処理１１５において、１または複数の基本周波数を推定結果として出力するようにしたが、演奏音の基本周波数の確率密度関数を推定結果として出力するようにしてもよい。この場合、類似度評価処理１１６では、目標基本周波数においてピークとなるガウス分布などの確率密度関数を生成し、この目標基本周波数の確率密度関数と基本周波数推定処理１１５により得られた基本周波数の確率密度関数の類似度を算出する。ここで、演奏位置において和音を演奏する場合には目標基本周波数は複数発生する。そのような場合には、複数の目標基本周波数において各々ピークとなる複数の確率密度関数を合成して目標基本周波数の確率密度関数を求める。演奏音の確率密度関数と目標基本周波数の確率密度関数とを用いた類似度の算出方法としては、例えば次のような各種の態様が考えられる。
【００６１】
（１−１）２つの確率密度関数の平均２乗誤差ＲＭＳ、すなわち、図４に示すように、演奏音の基本周波数の確率密度関数における確率密度と目標基本周波数の確率密度関数における確率密度との差分を２乗したものを全周波数帯域に亙って積分し、所定の定数Ｃにより除算したものの平方根の逆数を類似度とする。なお、この平方根の逆数の代わりに、所定の最大値からこの平方根を減算したものを類似度としてもよい。
【００６２】
（１−２）図５に示すように、周波数帯域を目標基本周波数の確率密度が高い有ピッチ領域と目標基本周波数の確率密度がほぼ０である無ピッチ領域とに分け、基本周波数推定処理１１５により得られた演奏音の基本周波数の確率密度関数において、有ピッチ領域に属する部分の確率密度の総和と、無ピッチ領域に属する部分の確率密度の総和とを求め、前者から後者を差し引いたものを類似度としてもよい。
【００６３】
（１−３）図６に示すように、演奏音の基本周波数の確率密度関数について、目標基本周波数を中心とした所定幅の周波数範囲における積分値を算出する。図示の例では、演奏位置において演奏すべき音が３個有り、その基本周波数がＦ１、Ｆ２、Ｆ３であるため、Ｆ１±ΔＦ、Ｆ２±ΔＦ、Ｆ３±ΔＦの各範囲について演奏音の確率密度関数の積分値（図示の斜線部分）を求めている。そして、各目標基本周波数を中心とした各範囲の各積分値を類似度として算出し、これが閾値を越えているか否かにより、各目標基本周波数の音を正しく演奏したか否かを判定する。この場合、演奏位置において演奏すべき音数が多いと、演奏音の確率密度関数において目標基本周波数の確率密度関数との類似度の低いピークが多数発生するため、実際に正しい演奏を行っていたとしても、正しい演奏が行われていないという誤った判定がなされるおそれがある。そこで、このような誤判定を防止するために、演奏位置において演奏すべき音の数がｋである場合に、目標基本周波数を中心とした範囲の積分値をｋ倍したものを類似度としてもよい。
【００６４】
（１−４）その他、演奏音の基本周波数の確率密度関数と目標基本周波数の確率密度関数から何らかの特徴量を抽出し、特徴量同士の積・べき乗・各種の数学関数などを用いた量を類似度とし、演奏音の基本周波数の確率密度関数と目標基本周波数の確率密度関数との差別化を容易にするようにしてもよい。
【００６５】
(１−５）上記のうちの例えば２種類の方法により２種類の類似度（第１および第２の類似度とする）を求め、第１および第２の類似度を線形結合した第３の類似度を演奏音が正しい音高であるか否かの判断のための類似度として採用してもよい。この場合、目標基本周波数通りの演奏音を発生させ、あるいは目標基本周波数からずれた基本周波数の演奏音を発生させる等、各種の条件で演奏音の発生および演奏音の基本周波数の推定を行わせる。そして、それら各種の条件において、第１の類似度と第２の類似度の重みを変えつつ、基本周波数の確率密度関数と目標基本周波数の確率密度関数とを用いた第３の類似度の算出を行う。そして、周知の判別分析の手法を用いることにより、演奏音の基本周波数と目標基本周波数とが一致しているのかずれているのかの差別化が容易な第３の類似度が得られるように、第１の類似度と第２の類似度の重みのバランスを定めればよい。なお、周知の判別分析の手法として知られる手法の他、ニューラルネットワークやＳＶＭ（サポートベクターマシン）として知られる手法も用いることができる。
【００６６】
（２）上記実施形態において、類似度評価処理１１６および一致判定処理１１７を実行する代わりに、基本周波数推定処理１１５において得られる基本周波数の確率密度関数の中から顕著なピークを選択し、そのピークが発生している基本周波数と目標基本周波数との一致度により、正しい音高で演奏が行われたか否かの判定を行うようにしてもよい。
【００６７】
（３）模範となる楽器演奏を録音することにより得られた音響信号のサンプルデータを楽曲データとして用い、この楽曲データに基本周波数推定処理を施すことにより、ユーザが発生すべき演奏音の目標基本周波数を求めるようにしてもよい。すなわち、図１において、収音部１０４により得られる演奏音の基本周波数を推定する基本周波数推定処理１１５の他に、演奏位置制御処理１１２により抽出される演奏位置の楽曲データ（模範演奏音のサンプルデータ）から模範演奏音の基本周波数を推定する基本周波数推定処理を設け、この基本周波数推定処理により推定された模範演奏音の基本周波数を目標基本周波数として用いるのである。この態様において、模範演奏の演奏音を収音部１０４により収音し、収音部１０４から得られる音響信号を模範演奏の楽曲データとして楽曲記憶部１０５に格納するようにしてもよい。
【図面の簡単な説明】
【００６８】
【図１】この発明による演奏判定装置の一実施形態である教習伴奏システムの構成を示すブロック図である。
【図２】同実施形態における基本周波数推定処理１１５の処理内容を示す図である。
【図３】同基本周波数推定処理１１５において行われるマルチエージェントモデルによる基本周波数の経時的な追跡を示す図である。
【図４】同実施形態における基本周波数の類似度の算出方法の変形例を示す図である。
【図５】同実施形態における基本周波数の類似度の算出方法の他の変形例を示す図である。
【図６】同実施形態における基本周波数の類似度の算出方法の他の変形例を示す図である。
【符号の説明】
【００６９】
１０１……制御部、１０２……操作部、１０３……ＨＤＤ、１０４……収音部、１０５……楽曲記憶部、１０６……データ入力部、１１１……楽曲入力処理、１１２……演奏位置制御処理、１１３……楽曲再生処理、１１４……楽曲表示処理、１１５……基本周波数推定処理、１１５Ｍ……音モデル、１１６……類似度評価処理、１１７……一致判定処理、１１８……判定結果表示処理、１……瞬時周波数の算出、２……周波数成分の候補の抽出、３……周波数帯域の制限、４ａ……メロディラインの推定、４ｂ……ベースラインの推定、４１……基本周波数の確率密度関数の推定、４２……マルチエージェントモデルによる基本周波数の継時的な追跡。
【図１】
2008209572.tif 000036
【図２】
2008209572.tif 000037
【図３】
2008209572.tif 000038
【図４】
2008209572.tif 000039
【図５】
2008209572.tif 000040
【図６】
2008209572.tif 000041