小節の先頭拍推定方法および先頭拍推定システム
【課題】楽曲の途中で拍数が変化した場合であっても各小節の先頭拍を推定できるようにした先頭拍の推定方法を提供する。
【解決手段】楽曲となる音響波形の入力を受け付ける楽曲受付手段１と、その受け付けた楽曲から拍の位置を推定する拍時刻推定手段２と、推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出する変動値算出手段３と、その変動値について４拍内の大小をスコアリングするスコアリング手段４と、すべての拍について一定周期の拍ごとに前記スコアリングされた値の時間変化を算出する時間変化算出手段５と、当該時間変化が最大となった位置を先頭拍の変動した位置として推定する変動位置推定手段６とを備える。
【解決手段】楽曲となる音響波形の入力を受け付ける楽曲受付手段１と、その受け付けた楽曲から拍の位置を推定する拍時刻推定手段２と、推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出する変動値算出手段３と、その変動値について４拍内の大小をスコアリングするスコアリング手段４と、すべての拍について一定周期の拍ごとに前記スコアリングされた値の時間変化を算出する時間変化算出手段５と、当該時間変化が最大となった位置を先頭拍の変動した位置として推定する変動位置推定手段６とを備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、楽曲から拍を推定し、その推定された拍からそれぞれの小節の先頭拍となる位置を推定できるようにした先頭拍推定方法に関するものである。
【背景技術】
【０００２】
従来より、入力された楽曲の中から拍の位置を推定するとともに、その推定された拍の中から小節の先頭拍の位置を推定できるようにした技術が提案されている。このような小節の先頭拍の位置を推定することができれば、楽曲同士をシンクロさせながら遷移させたり、あるいは、楽曲と映像とを同期させて映像を表示させたりするようなアプリケーションに適用することができる（特許文献１）。
【０００３】
このような小節の先頭拍の位置を推定する方法に関して、従来の推定方法について説明する。
【０００４】
まず、楽曲の小節の先頭拍を推定する場合、入力された楽曲の拍の位置を推定する。この拍の位置を推定する方法としては、例えば、ＰＣＭデータから得られた楽曲波形信号に対してＦＦＴを用いてパワースペクトル変化率を計算し、音が出力されていると推定された位置を拍時刻とする方法や（特許文献２）、ＦＦＴを用いてパワースペクトルを計算し、自己相関関数および相互相関関数を用いることにより拍時刻を推定するようにした方法（特許文献３）などがある。
【０００５】
次に、例えば、４／４拍子の場合、この推定された拍を４拍ごとに、小節内における拍の順序を表す１〜４のＩＤを暫定的に割り当て、各拍時刻前後の隣接する分析フレームのスペクトルを求めて、そのスペクトルの変動量を算出する。そして、すべての小節における第m（＝１〜４）拍目におけるスペクトルの変動量を加算する。すなわち、第１小節から最終小節における第１拍目のスペクトルの変動量を加算し、第１小節から最終小節における第２拍目のスペクトルの変動量を加算し、第１小節から最終小節における第３拍目のスペクトルの変動量を加算し、第１小節から最終小節における第４拍目のスペクトルの変動量を加算する。そして、その加算されたそれぞれのスペクトルの変動量のうち、最も変動量の大きな拍の位置を先頭拍として推定し、このように推定された先頭拍から４拍ごとの拍の位置を、各小節における先頭拍として推定するようにしている（非特許文献１）。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】Matthew E. P. Davies and Mark D. Plumbley, "A spectral difference approach to downbeat extraction in musical audio", EUSIPCO 2006, Florence, Italy, September 4-8, (2006).
【特許文献】
【０００７】
【特許文献１】特開２００２−２６０３４０号公報
【特許文献２】特開２００７−１９９３０６号公報
【特許文献３】特開２００７−２４８８９５号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、上述のように先頭拍を推定する方法では、次に示すような問題点を生ずる。
【０００９】
すなわち、上述のような先頭拍を推定する方法では、楽曲の途中で小節内の拍数が変化した場合には対応することができない。具体的には、４／４拍子の楽曲の途中で１小節内の拍数が５拍もしくは６拍に変化した場合、従来の方法では、その変更された位置からの先頭拍の位置がずれてしまう。また、一部の小節の拍数が変更される場合に限らず、楽曲の途中から一定区間３／４拍子に変更されるような場合にも、同様にその位置から先頭拍の位置がずれてしまう。
【００１０】
そこで、本発明は、楽曲の途中で拍数が変化した場合であっても各小節の先頭拍を推定できるようにした先頭拍の推定方法を提供することを目的とする。
【課題を解決するための手段】
【００１１】
すなわち、本発明は上記課題を解決するために、楽曲の入力を受け付けるステップと、当該受け付けた楽曲から拍の位置を推定するステップと、当該推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出するステップと、すべての拍について一定周期の拍ごとに前記変動値の時間変化を算出するステップと、当該時間変化の変化量が最大となった位置、もしくは、時間変化の変化量が所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するステップとを備えるようにしたものである。
【００１２】
このようにすれば、楽曲の途中で小節内の拍数が変化した場合、その時点で変動値の時間変化が大きく変動するため、先頭拍の変化した位置を推定することができるようになる。これにより、その変化した位置から新たに先頭拍を推定すれば、楽曲の途中で小節内の拍数が変化した場合であっても、すべての小節の先頭拍を推定することができるようになる。
【００１３】
また、このような発明において、当該推定された拍に対して、周期的な識別符号を付与するステップと、識別符号の付与されたすべての拍における前後の分析フレームを周波数解析してスペクトルを算出し、当該算出されたスペクトルの変動値を算出するステップと、前記周期的に付与された識別符号の範囲内で前記変動値を順位付けするステップと、当該順位付けされた値について、前記識別符号ごとに時間軸方向に加算していくステップと、当該時間軸方向に加算された値に対して時間変化を算出するステップと、当該時間変化の変化量が最大となった位置、もしくは、時間変化の変化量が所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するステップとを備えるようにする。
【００１４】
このようにすれば、変動値をスコアに置き換えて時間変化をみるようにするので、先頭拍以外の拍における突発的に発生する大きなスペクトルの変動の影響を低減することができるようになる。
【００１５】
さらには、前記時間軸方向に加算された値に対して、回帰直線に関する値を算出するステップと、当該算出された回帰直線に関する値と前記時間軸方向に加算された値との差分を算出するステップを備え、当該算出された差分の値が最大となった位置、もしくは、所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定する。
【００１６】
このようにすれば、先頭拍の位置が変動した場合、回帰直線との差分が大きくなるため、これによって先頭拍の変動した位置を容易に推定することができるようになる。
【発明の効果】
【００１７】
本発明によれば、楽曲の入力を受け付けるステップと、当該受け付けた楽曲から拍の位置を推定するステップと、当該推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出するステップと、すべての拍について一定周期の拍ごとに前記変動値の時間変化を算出するステップと、当該時間変化の変化量が最大となった位置、もしくは、時間変化の変化量が所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するステップとを備えるようにしたので、楽曲の途中で先頭拍の変化した位置を推定することができるようになる。これにより、その変化した位置から新たに先頭拍を推定すれば、楽曲の途中で小節内の拍数が変化した場合であっても、すべての小節の先頭拍を推定することができるようになる。
【図面の簡単な説明】
【００１８】
【図１】本発明の一実施の形態を示す先頭拍推定システムの機能ブロック図
【図２】推定された拍時刻上でスペクトルが大きく変動する周期（PLSV）を示す図
【図３】同形態における拍内のスペクトラルフラックスの順位付けスコアを付与する図
【図４】同形態における順位付けスコアの部分和を時系列で示した図
【図５】同形態における回帰直線との差を示す図
【図６】同形態における処理の流れを示すフローチャート
【発明を実施するための最良の形態】
【００１９】
以下、本発明の一実施の形態について図面を参照しながら説明する。この実施の形態における先頭拍推定システム１００は、一台のコンピューターもしくは複数のコンピューターによって構成されるものであって、図１に示すように、楽曲の入力を受け付ける楽曲受付手段１と、この受け付けた楽曲から拍の位置を推定する拍時刻推定手段２と、その推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出する変動値算出手段３とを備えてなるもので、特徴的には、すべての拍について、例えば、４／４拍子の場合は４拍目ごとの変動値の時間変化を算出する時間変化算出手段５と、その時間変化が最大となった位置を先頭拍（Downbeat）が変動した位置として推定する変動位置推定手段６とを備えるようにしたものである。そして、これによって、先頭拍の変動した位置を推定し、そこから新たに先頭拍を推定できるようにしている。以下、本実施の形態の各機能実現手段について詳細に説明する。なお、これらの各機能実現手段は、コンピューターを構成するCPU、メモリやハードディスクなどの記憶デバイス、マイクやスピーカー、ディスプレイなどのハードウェアと、そのハードディスクに格納されたプログラムとを協働させて実現される。
【００２０】
まず、楽曲受付手段１は、ユーザーから楽曲の入力を受け付ける。ユーザーから楽曲を受け付ける場合は、インターネットなどを介して楽曲をダウンロードして受け付けるようにしてもよく、あるいは、ＣＤなどの記憶媒体から楽曲を受け付けるようにしてもよい。このように楽曲を受け付ける場合は、音響波形として楽曲を受け付ける。
【００２１】
拍時刻推定手段２は、このように受け付けた楽曲の中から拍となる位置の時刻を推定する。この拍となる位置を推定するには、既存の種々の手法を用いることができるが、ここでは一例として、次のような処理を行うようにする。
【００２２】
まず、受け付けた楽曲から拍時刻の推定する場合、まず受け付けた音響波形を平滑化し、その平滑化された波形信号からアクセントとなる時刻を抽出する。この平滑化された振幅包絡には、アクセントとは考えられない微小なピークが多数存在するため、アクセントの時刻を正しく検出できない場合がある。このため、所定の閾値以下の微小なパワーを除いておおまかなアクセントの時刻を抽出する。この閾値を用いる場合、例えば、非常に激しい演奏と非常に穏やか演奏が混在している場合、非常に穏やかな演奏が行われている区間のすべてが閾値以下となってしまい、その区間のアクセントを検出することができなくなる。このため、微小なパワーの除去する場合は、各区間ごとに最適な値に変化する閾値を用いるようにする。次に、平滑化された波形情報をフーリエ変換してパワースペクトルが最大となる周期を求め、先にアクセント時刻と推定された時刻からこの周期時刻ごとに、その周期の時刻の近傍に他のアクセントが存在しているかどうかを検出する。そして、その周期の時刻の近傍に他のアクセントが存在している場合は、その時刻を拍時刻として推定する。また、これ以外の方法として、特開2008−107569号に示されるように、音響波形に対して全波整流を行った後、フィルタによって所定の周波数帯域の信号を抽出し、ピークの存在する時刻を拍時刻として推定する方法なども用いることができる。そして、このように拍時刻を推定した後、これらすべての拍についてBｊ（ｊ＝0,1,2,…）とラベリングする。
【００２３】
次に、このように楽曲の拍時刻を推定した後、この拍から小節における先頭拍を推定する。なお、ここでは楽曲として４／４拍子の楽曲を想定して説明するが、３／４拍子など他の拍子についても同様の処理を適用することができる。
【００２４】
先頭拍を推定する場合、まず、変動値算出手段３によって、各拍時刻における前後の微小時間における周波数の変動を算出する。このような処理を行うのは、各小節の先頭拍の時刻上では和音やメロディが大きく変化し、スペクトルが大きく変動する傾向が強いからである。ここでは、まず、先の拍時刻推定手段２で推定された拍Bｊ（ｊ＝0,1,2,…）について、周期的な４拍ごとに小節内のＩＤを割り当て、例えば、第k小節（k＝0,1,2,…）のm番目（m＝1,2,3,4）の拍をBm（k)とする。そして、Bm（k）と割り当てられた拍の前後のフレームからＦＦＴを用いてパワースペクトルを算出し、その前後におけるフレームのパワースペクトルの差分であるスペクトラムフラックスを算出する。なお、このスペクトラムフラックスが、本発明の変動値に相当する。
【００２５】
次に、時間変化算出手段５では、この変動値算出手段３によって算出された各拍の前後におけるスペクトルの変動値について各mの値毎に時間変化を算出する。すなわち、各小節における１拍目のスペクトラムフラックスの時間変化、各小節における２拍目のスペクトラムフラックスの時間変化、各小節における３拍目のスペクトラムフラックスの時間変化、各小節における４拍目のスペクトラムフラックスの時間変化を算出する。ここで、暫定的に決定されたk番目の小節におけるB1（k）、B2（k）、B3（k）、B4（k）のスペクトルの変動量を比較すると、多少の順位の入れ替わりはあるものの、各小節における先頭拍に対応するmのＩＤを持つ拍が常に最も大きいスペクトルの変動量を示すと考えられる。ここで、人間が主観的に感じる先頭拍（Downbeat）の周期をPSD（Periodicity of subjective downbeats）、音響波形上のスペクトルが大きく変動する周期をPLSV（Periodicity of the Large Spectral Variation）と呼ぶことにし、図２にPSDが変化しない場合、および、PSDが変化する場合の各Bm（k）におけるスペクトルの変動量の部分和（すなわち、各小節におけるk拍目の和）を示す。図２の上図は、拍時刻Bm（k）のスペクトルの変動量を算出し、mのＩＤ別の部分和をプロットしたグラフである。一方、図２の下図は、上図と同一の楽曲で、PSDが一時的に変化する場合の各Bm（k）におけるスペクトルの変動量からmのＩＤ別の部分和をプロットしグラフである。ここで、スペクトルの変動量の指標には、数１におけるスペクトラルフラックスを用いることとした。
【００２６】
【数１】
【００２７】
数１におけるNは、ＦＦＴに用いたサンプル点数で、8192点（すなわち、約185.8msec）を用いている。hは周波数解析を行う際のフレームのＩＤ、Lh[n]はh番目のフレームにおける周波数nのパワーを表している。図２の上図に示すように、特定の拍（図２ではm＝2が先頭拍となっている）におけるスペクトルの変動量が突出して増加しているのが確認できる。そして、図２のグラフの作成に用いた楽曲では、図２の上図で最大値を示しているm＝2と対応する拍B2（k）が先頭拍と一致していることを確認することができた。このように、PSDが変化しない楽曲に対しては、mのＩＤ別にスペクトルの変動量の合計値を求めることで先頭拍をある程度推定可能であることが分かる。ここで、図２の上図において、PSDが一時的に変化した場合について考察すると、図２の上図におけるm＝2の拍時刻においてスペクトルの変動量が減少し、先頭拍と対応する新たなmのＩＤを持つ拍時刻において、スペクトルの変動量が増加すると考えられる。図２の下図は、図２の上図で用いた楽曲と同一の楽曲に対して、本来存在しない拍の音響波形を意図的に挿入し、図２の上図と同様の方法でグラフを作成した結果である。意図的に拍の音響波形を挿入した時刻は、楽曲のちょうど半分の時刻（図２下図のChanging time）である。意図的に挿入された拍の音響波形が原因となり、PSDが変化し、本来２拍目であった先頭拍が、３拍目となっていることが確認できる。これは、意図的に挿入された拍の音響波形が原因で、PLSVが一時的に変化し、その時刻の前半部分までは、m＝2の拍時刻と実際の先頭拍が対応していたが、後半部分では、先頭拍と対応するmが3へと変化したためであると言える。このように、楽曲中においてPSDが一時的に変化した場合、該当する時刻を検出し、その時刻を基準とした区間ごとに先頭拍を推定する必要がある。
【００２８】
ところで、先頭拍以外の拍において、突発的に大きなスペクトルの変動が発生した場合、該当する拍におけるmのスペクトルの変動量の増加傾向に影響を与える可能性がある。そこで、本実施の形態では、この対策として、拍時刻Bm（k）におけるスペクトルの変動量Fm（k）へBm（k−1）およびBm（k+1）におけるスペクトルの変動量Fm（k−1）およびFm(k+1)を加算し、平均を求めることで、先頭拍以外の拍における突発的に発生する大きなスペクトルの変動の影響を低減することとした。また、この処理によって得られる拍時刻Bm（k）におけるスペクトルの変動量の平均値をAm（k）と表現する。Am（k）の算出式を数２に示す。
【００２９】
【数２】
【００３０】
ここで、Kは、暫定的に決定された小節における最終小節のＩＤを表している。
【００３１】
また、Am（k）は、フペクトラルフラックスに基づいて算出されているため、例えば、スネアドラムなどの全周波数帯域にパワーをもつ音がいくつかの小節にわたって発音されていた場合に、連続して非常に大きな値を示す場合がある。したがって、スペクトルの変動量のおおまかな増加傾向を、楽曲の局所的な物理的特徴に依存するAm(k)から直接求めることはあまり望ましくない。これは、PSDが変化せず、先頭拍と対応する拍において、スペクトルの変動量が常に最も大きい値を示す場合でも、スペクトルの変動量の増加傾向が一定になるとは限らないためである。この問題は、PSDが一時的に変化した時刻の特徴を得る際に悪影響を与え、最終的にはPSDが一時的に変化した時刻の誤検出の原因になると考えられる。この問題を解決するため、スペクトルの変動量の増加傾向を求める際に、Am（k）をそのまま用いるのではなく、暫定的に決定された各小節内のAm（k）の値を、スコアリング手段４によって順位別のスコア（4,3,2,1）に置き換え、その順序に基づいたスコアを付与することで、スペクトルの変動量の差を無視することとした。この順位スコアを先頭拍の可能性を表す指標と解釈し、「Downbeatness」と呼ぶとともに、拍時刻Bm（k）におけるDownbeatnessの値をDm（k）と表現する。
【００３２】
図３にDownbeatnessを算出する際の流れを示す。図３に示すように、暫定的に決定されたk番目の小節におけるA1（k）、A2（k）、A3（k）、A4（k）を比較し、最も大きい値を持つm番目の拍時刻Bm（k)から順に、4,3,2,1と順位スコアの値を付与していく。こうすることで、スペクトルの変動の急激な増加、あるいは減少の影響を低減することができる。
【００３３】
次に、変動位置推定手段６は、このように時系列に並べられたスペクトルの変動値に基づいて、先頭拍の変動した位置を推定する。
【００３４】
この変動位置推定手段６では、図２で用いた楽曲からDownbeatnessを算出し、mのＩＤ別の部分和を算出する。ここで得られた、各拍時刻B m（k）におけるmのＩＤ別の部分和の時間系列を、先頭拍らしさ関数と呼び、Sm（k）と表現する。先頭拍らしさ関数Sm（k）の算出式を数３に示す。
【００３５】
【数３】
【００３６】
図４に、図２で用いた楽曲から先頭拍らしさ関数を算出した結果を示す。図４に示すように、Downbeatnessは、各mにおけるスペクトルの変動量の増加傾向を良く表していることが確認できる。特に、図４の下図における破線の円内において、S3（k）が急激に増加していることが確認できる。そして、Sm（k）の局所的な傾きが急激に変化した地点は、Sm（k）の回帰直線とSm（k）の距離が遠い地点と見なすことができる。Sm（k）の回帰直線をCm（k）とすると、暫定的に決定されたk番目の小節におけるm番目のSm（k）と回帰直線Cm（k）との距離Rm（k）は数４で求めることができる。
【００３７】
【数４】
【００３８】
したがって、Sm（k）の局所的な傾きが急激に変化した地点、すなわちPSDが一時的に変化した可能性が高い時刻は、Rm（k）のピークから得られると期待できる。このとき、PSDが一時的に変化した時刻においては、特定のmのみではなく、全mにおけるSm（k）において、スペクトルの変動量の増加傾向が変化する可能性が高い。したがって、PSDが一時的に変化した時刻を得るためにRm（k）のピークを求める場合は、m＝1,2,3,4についてそれぞれRm（k）を加算し、ピークとなる時刻を求める処理が妥当であると言える。
【００３９】
図５に、図４のPSDが変化する場合の楽曲からRm（k）を抽出した結果を示す。図５に示すように、m＝1,2,3,4のRm（k）を加算した時間関数（図５の細い破線）における中間部分のピークと、楽曲のPSDが変化した時刻が概ね一致していることが確認できる。
【００４０】
ここでは、Rm（k）を加算した時間関数を先頭拍差異検出関数と呼ぶこととし、Rsum（k）と表現する。
【００４１】
図５から確認できるように、Rsum（k）は、PSDが変化した時刻と良い対応を持つ。そこで、楽曲中の拍時刻に基づいて、このRsum（k）を算出し、その最大値であるピークを検出することで、PSDが変化した時刻を検出し、この時刻を先頭拍の変化した時刻と推定する。そして、このように先頭拍の変化した時刻を推定できた場合、その拍から新たに先頭拍の推定処理を同様に行っていくようにする。
【００４２】
次に、このように構成された先頭拍推定システム１００を用いて楽曲の先頭拍を推定する方法について、図６のフローチャートを用いて説明する。
【００４３】
まず、先頭拍を推定する楽曲の音響波形の入力を受け付け（ステップＳ１）、これを記憶部に格納しておく。
【００４４】
次に、その受け付けた音響波形から拍となる位置の時刻を推定する（ステップＳ２）。この受け付けた音響波形から拍時刻を推定する場合は、入力された音響波形を平滑化し、所定の閾値以下の微小なパワーを除いたおおまかなアクセント時刻を移動平均によって求める。一方、これとは別に、平滑化された波形情報をフーリエ変換してパワースペクトルが最大となる周期を求め、先におおまかに推定されたアクセント時刻からこの周期時刻ごとに、その周期時刻の近傍に他のアクセントが存在しているかどうかを検出する。そして、その周期時刻の近傍に他のアクセントが存在している場合は、その時刻を拍時刻として推定する（ステップＳ２）。
【００４５】
次に、このように推定されたすべての拍についてBｊ(ｊ＝0,1,2,…)とラベリングするとともに、周期的な拍（ここでは４拍）ごとにＩＤ（Bm(k)）を割り当てるラベリング処理を行う（ステップＳ３）。
【００４６】
そして、このようにBm（k）とラベリングされたすべての拍について、その拍時刻の前後のフレームからＦＦＴを用いてパワースペクトルを算出し、そのパワースペクトルの差分であるスペクトルフラックス（すなわち、スペクトルの変動値）を算出する（ステップＳ４）。
【００４７】
次に、このスペクトルフラックスを各ＩＤごとに時系列に並べるようにするが、ここでは、突発的なスペクトルの変動を無視できるように、スコアリング手段４によって周期的に割り当てられたＩＤの範囲内（すなわち、k＝1,2,3,4の範囲内）でスペクトルフラックスの大きいものから順に4,3,2,1とスコアを与える（ステップＳ５）。そして、このように与えられたスコアを各ＩＤごとに時系列に加算して先頭拍らしさ関数Sm（k）を算出し、図４に対応する数値列を取得する（ステップＳ６）。
【００４８】
そして、この先頭拍らしさ関数Sm（k）から、スペクトルの変動量が急激に変化した位置を算出すべく、Sm（k）に対して単一の回帰直線Cm（k）を算出し（ステップＳ７）、その回帰直線Rm（k）との距離を数４により求めて、その値を加算した先頭拍差異検出関数を取得する。そして、その最大値であるピークを抽出して、その時刻を先頭拍の変化した時刻と推定する（ステップＳ８）。
【００４９】
このように上記実施の形態によれば、楽曲の入力を受け付ける楽曲受付手段１と、当該受け付けた楽曲から拍の位置を推定する拍時刻推定手段２と、当該推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出する変動値算出手段３と、すべての拍について一定周期の拍ごとに前記変動値の時間変化を算出する時間変化算出手段５と、当該時間変化が最大となった位置を先頭拍の変動した位置として推定する変動位置推定手段６とを備えるようにしたので、先頭拍の変化した位置を推定することができるようになる。これにより、その変化した位置から新たに先頭拍を推定すれば、楽曲の途中で小節内の拍数が変化した場合であっても、すべての小節の先頭拍を推定することができるようになる。
【００５０】
また、周期的に付与されたＩＤ（k＝1,2,3,4）の範囲内で変動値を順位付けし、その順位付けされた値に基づいて時間変化を算出するようにしたので、先頭拍以外の拍における突発的に発生する大きなスペクトルの変動の影響を低減することができるようになる。
【００５１】
さらには、時間軸方向に加算された値に対して、回帰直線に関する値を算出し、その回帰直線に関する値と時間軸方向に加算された値との差分を算出して、その差分が最大となった位置を先頭拍の変動した位置と推定するようにしたので、これによって先頭拍の変動した位置を容易に推定することができるようになる。
【００５２】
なお、本発明は上記実施の形態に限定されることなく種々の態様で実施することができる。
【００５３】
例えば、上記実施の形態では、スペクトルフラックスの時間変化を算出する場合、スコアリング処理を行った後に時間変化を算出するようにしたが、スペクトルフラックスの状態で時間変化を算出するようにしてもよい。
【００５４】
また、上記実施の形態では、時系列にスペクトルの変動値を加算して回帰直線との差を求めるようにしたが、時系列にスペクトルの変動値を加算するのではなく、それぞれの変動値そのものを時系列に並べ、スコアリングされた各ＩＤの時系列の値との大小関係をみることによって変動位置を推定するようにしてもよい。
【００５５】
さらには、上記実施の形態では、周期的なＩＤとして４／４拍子に対応した４拍ごとにＩＤを付与するようにしたが、３／４拍子の場合は、３拍ごとにＩＤを付与するようにしてもよい。
【００５６】
また、上記実施の形態では、回帰直線との距離Rm（k）が最大となる位置を変動位置と推定するようにしたが、複数回変動位置が生ずる可能性もあるため、所定の閾値を設定しておき、その閾値よりもRm（k）の値が大きい場合は、先頭拍が変動した位置と推定するようにしてもよい。
【実施例１】
【００５７】
上記手法によって、推定された先頭拍が正しいかどうかを評価した。上記手法の評価基準の一つとして、PSDが変化しない場合の楽曲において、提案手法と既存手法とを比較した際の推定精度への影響が挙げられる。これは、PSDが変化しない楽曲において、提案手法を用いて先頭拍を推定した場合の誤推定が起こる割合を評価することができる。この評価を行うためには、楽曲に対して人手で拍時刻および先頭拍の時刻のラベリングデータを作成し、ラベリングされた拍時刻に基づいて、既存手法および提案手法によって先頭拍を推定し、正解率を比較すれば良いと言える。以降この評価基準の条件を、S（Simple） 条件と呼ぶ。
【００５８】
もう一つの評価基準として、PSDが変化する場合の楽曲および推定された拍時刻上においてPLSVが変化する場合の楽曲において、提案手法と既存手法とを比較した際の推定精度の違いが挙げられる。すなわち、既存手法では推定が困難な楽曲に対する提案手法の有効性を評価することができる。
【００５９】
PSDが変化する状況としては、楽曲の小節の長さが一時的に変化する状況が考えられる。また、推定された拍時刻上においてPLSVが変化する状況としては、拍時刻の誤推定が発生した状況が挙げられる。これらの状況下において、既存手法および提案手法によって先頭拍を推定し、正解率を比較すれば良いと言える。本実施例では、これらの状況を意図的に再現することで、既存手法では推定が困難な楽曲に対する提案手法の有効性を評価することとした。具体的には、楽曲に対して人手で拍時刻および先頭拍の時刻のラベリングデータを作成し、ある単一の小節に対して、本来存在しない拍の音響波形を挿入することで、小節の長さが一時的に変化する状況を再現することとした。さらに、ラベリングされた拍時刻に対して、意図的に誤った拍時刻をラベリングすることで、拍時刻の誤推定が発生した状況を再現することとした。以降この評価基準の条件を、C（Challenge）条件と呼び、特に小節の長さが一時的に変化する条件をC m （measure）条件、拍時刻の誤推定が発生する条件をC b （beat）条件と呼ぶこととする。
【００６０】
上記各条件で評価を行うために、まず聴取者にRWC データベースから抜粋した楽曲を聴取させ、一定の間隔で発音される正弦波の時間間隔を拍時刻および先頭拍と同期するようにそれぞれ調整させた。そして、楽曲の拍時刻および先頭拍の時刻において発音される正弦波の時間系列データのwaveファイルをそれぞれ作成させ、それらのデータを正解データとした。
【００６１】
この正解データと、楽曲から推定された先頭拍の時刻を比較する必要がある。そこで、推定された先頭拍の時刻で正弦波が発音されるよう合成音を作成し、正解データとのパワーの比率を正解率Qとして求めることで評価を行うこととした。この正解率Qの算出方法を数５に示す。
【００６２】
【数５】
【００６３】
ここで、giは正解データの音響信号におけるパワー関数、eiは推定された先頭拍の時刻で正弦波が発音される音響信号におけるパワー関数を表している。
【００６４】
表1に、評価結果を示す。表1におけるEは既存手法、Pは提案手法を表している。表1のS条件において、提案手法と比較して、既存手法の方が僅かに高い推定精度であることについて言及すると、評価に用いた楽曲がPSDの変化しない楽曲であったためであると考えられる。次に、表1のC条件に示す通り、楽曲の中盤でPSDが変化したり（Cm条件）、拍の誤検出が発生すると（Cｂ条件）、既存手法の推定率が大きく低下しているのが確認できる。それに対して、提案手法を用いた場合、PSDが変化しない場合と同等の推定精度で先頭拍を推定できていることが確認できる（表1内の下線部）。したがって、PSDが変化する楽曲、あるいは拍時刻の誤推定の発生により、推定された拍時刻上においてPLSVが変化するような楽曲に対して、提案手法が有効であることを確認できた。このため、提案手法を用いることで、小節の長さが変化する楽曲や、本来拍が存在しない時刻に拍時刻が挿入されている場合でも、それらの影響を受けずに先頭拍を推定できることが確認できた。
【００６５】
【表１】
【符号の説明】
【００６６】
１・・・楽曲受付手段
２・・・拍時刻推定手段
３・・・変動値算出手段
４・・・スコアリング手段
５・・・時間変化算出手段
６・・・変動位置推定手段
１００・・・先頭拍推定システム
【技術分野】
【０００１】
本発明は、楽曲から拍を推定し、その推定された拍からそれぞれの小節の先頭拍となる位置を推定できるようにした先頭拍推定方法に関するものである。
【背景技術】
【０００２】
従来より、入力された楽曲の中から拍の位置を推定するとともに、その推定された拍の中から小節の先頭拍の位置を推定できるようにした技術が提案されている。このような小節の先頭拍の位置を推定することができれば、楽曲同士をシンクロさせながら遷移させたり、あるいは、楽曲と映像とを同期させて映像を表示させたりするようなアプリケーションに適用することができる（特許文献１）。
【０００３】
このような小節の先頭拍の位置を推定する方法に関して、従来の推定方法について説明する。
【０００４】
まず、楽曲の小節の先頭拍を推定する場合、入力された楽曲の拍の位置を推定する。この拍の位置を推定する方法としては、例えば、ＰＣＭデータから得られた楽曲波形信号に対してＦＦＴを用いてパワースペクトル変化率を計算し、音が出力されていると推定された位置を拍時刻とする方法や（特許文献２）、ＦＦＴを用いてパワースペクトルを計算し、自己相関関数および相互相関関数を用いることにより拍時刻を推定するようにした方法（特許文献３）などがある。
【０００５】
次に、例えば、４／４拍子の場合、この推定された拍を４拍ごとに、小節内における拍の順序を表す１〜４のＩＤを暫定的に割り当て、各拍時刻前後の隣接する分析フレームのスペクトルを求めて、そのスペクトルの変動量を算出する。そして、すべての小節における第m（＝１〜４）拍目におけるスペクトルの変動量を加算する。すなわち、第１小節から最終小節における第１拍目のスペクトルの変動量を加算し、第１小節から最終小節における第２拍目のスペクトルの変動量を加算し、第１小節から最終小節における第３拍目のスペクトルの変動量を加算し、第１小節から最終小節における第４拍目のスペクトルの変動量を加算する。そして、その加算されたそれぞれのスペクトルの変動量のうち、最も変動量の大きな拍の位置を先頭拍として推定し、このように推定された先頭拍から４拍ごとの拍の位置を、各小節における先頭拍として推定するようにしている（非特許文献１）。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】Matthew E. P. Davies and Mark D. Plumbley, "A spectral difference approach to downbeat extraction in musical audio", EUSIPCO 2006, Florence, Italy, September 4-8, (2006).
【特許文献】
【０００７】
【特許文献１】特開２００２−２６０３４０号公報
【特許文献２】特開２００７−１９９３０６号公報
【特許文献３】特開２００７−２４８８９５号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、上述のように先頭拍を推定する方法では、次に示すような問題点を生ずる。
【０００９】
すなわち、上述のような先頭拍を推定する方法では、楽曲の途中で小節内の拍数が変化した場合には対応することができない。具体的には、４／４拍子の楽曲の途中で１小節内の拍数が５拍もしくは６拍に変化した場合、従来の方法では、その変更された位置からの先頭拍の位置がずれてしまう。また、一部の小節の拍数が変更される場合に限らず、楽曲の途中から一定区間３／４拍子に変更されるような場合にも、同様にその位置から先頭拍の位置がずれてしまう。
【００１０】
そこで、本発明は、楽曲の途中で拍数が変化した場合であっても各小節の先頭拍を推定できるようにした先頭拍の推定方法を提供することを目的とする。
【課題を解決するための手段】
【００１１】
すなわち、本発明は上記課題を解決するために、楽曲の入力を受け付けるステップと、当該受け付けた楽曲から拍の位置を推定するステップと、当該推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出するステップと、すべての拍について一定周期の拍ごとに前記変動値の時間変化を算出するステップと、当該時間変化の変化量が最大となった位置、もしくは、時間変化の変化量が所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するステップとを備えるようにしたものである。
【００１２】
このようにすれば、楽曲の途中で小節内の拍数が変化した場合、その時点で変動値の時間変化が大きく変動するため、先頭拍の変化した位置を推定することができるようになる。これにより、その変化した位置から新たに先頭拍を推定すれば、楽曲の途中で小節内の拍数が変化した場合であっても、すべての小節の先頭拍を推定することができるようになる。
【００１３】
また、このような発明において、当該推定された拍に対して、周期的な識別符号を付与するステップと、識別符号の付与されたすべての拍における前後の分析フレームを周波数解析してスペクトルを算出し、当該算出されたスペクトルの変動値を算出するステップと、前記周期的に付与された識別符号の範囲内で前記変動値を順位付けするステップと、当該順位付けされた値について、前記識別符号ごとに時間軸方向に加算していくステップと、当該時間軸方向に加算された値に対して時間変化を算出するステップと、当該時間変化の変化量が最大となった位置、もしくは、時間変化の変化量が所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するステップとを備えるようにする。
【００１４】
このようにすれば、変動値をスコアに置き換えて時間変化をみるようにするので、先頭拍以外の拍における突発的に発生する大きなスペクトルの変動の影響を低減することができるようになる。
【００１５】
さらには、前記時間軸方向に加算された値に対して、回帰直線に関する値を算出するステップと、当該算出された回帰直線に関する値と前記時間軸方向に加算された値との差分を算出するステップを備え、当該算出された差分の値が最大となった位置、もしくは、所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定する。
【００１６】
このようにすれば、先頭拍の位置が変動した場合、回帰直線との差分が大きくなるため、これによって先頭拍の変動した位置を容易に推定することができるようになる。
【発明の効果】
【００１７】
本発明によれば、楽曲の入力を受け付けるステップと、当該受け付けた楽曲から拍の位置を推定するステップと、当該推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出するステップと、すべての拍について一定周期の拍ごとに前記変動値の時間変化を算出するステップと、当該時間変化の変化量が最大となった位置、もしくは、時間変化の変化量が所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するステップとを備えるようにしたので、楽曲の途中で先頭拍の変化した位置を推定することができるようになる。これにより、その変化した位置から新たに先頭拍を推定すれば、楽曲の途中で小節内の拍数が変化した場合であっても、すべての小節の先頭拍を推定することができるようになる。
【図面の簡単な説明】
【００１８】
【図１】本発明の一実施の形態を示す先頭拍推定システムの機能ブロック図
【図２】推定された拍時刻上でスペクトルが大きく変動する周期（PLSV）を示す図
【図３】同形態における拍内のスペクトラルフラックスの順位付けスコアを付与する図
【図４】同形態における順位付けスコアの部分和を時系列で示した図
【図５】同形態における回帰直線との差を示す図
【図６】同形態における処理の流れを示すフローチャート
【発明を実施するための最良の形態】
【００１９】
以下、本発明の一実施の形態について図面を参照しながら説明する。この実施の形態における先頭拍推定システム１００は、一台のコンピューターもしくは複数のコンピューターによって構成されるものであって、図１に示すように、楽曲の入力を受け付ける楽曲受付手段１と、この受け付けた楽曲から拍の位置を推定する拍時刻推定手段２と、その推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出する変動値算出手段３とを備えてなるもので、特徴的には、すべての拍について、例えば、４／４拍子の場合は４拍目ごとの変動値の時間変化を算出する時間変化算出手段５と、その時間変化が最大となった位置を先頭拍（Downbeat）が変動した位置として推定する変動位置推定手段６とを備えるようにしたものである。そして、これによって、先頭拍の変動した位置を推定し、そこから新たに先頭拍を推定できるようにしている。以下、本実施の形態の各機能実現手段について詳細に説明する。なお、これらの各機能実現手段は、コンピューターを構成するCPU、メモリやハードディスクなどの記憶デバイス、マイクやスピーカー、ディスプレイなどのハードウェアと、そのハードディスクに格納されたプログラムとを協働させて実現される。
【００２０】
まず、楽曲受付手段１は、ユーザーから楽曲の入力を受け付ける。ユーザーから楽曲を受け付ける場合は、インターネットなどを介して楽曲をダウンロードして受け付けるようにしてもよく、あるいは、ＣＤなどの記憶媒体から楽曲を受け付けるようにしてもよい。このように楽曲を受け付ける場合は、音響波形として楽曲を受け付ける。
【００２１】
拍時刻推定手段２は、このように受け付けた楽曲の中から拍となる位置の時刻を推定する。この拍となる位置を推定するには、既存の種々の手法を用いることができるが、ここでは一例として、次のような処理を行うようにする。
【００２２】
まず、受け付けた楽曲から拍時刻の推定する場合、まず受け付けた音響波形を平滑化し、その平滑化された波形信号からアクセントとなる時刻を抽出する。この平滑化された振幅包絡には、アクセントとは考えられない微小なピークが多数存在するため、アクセントの時刻を正しく検出できない場合がある。このため、所定の閾値以下の微小なパワーを除いておおまかなアクセントの時刻を抽出する。この閾値を用いる場合、例えば、非常に激しい演奏と非常に穏やか演奏が混在している場合、非常に穏やかな演奏が行われている区間のすべてが閾値以下となってしまい、その区間のアクセントを検出することができなくなる。このため、微小なパワーの除去する場合は、各区間ごとに最適な値に変化する閾値を用いるようにする。次に、平滑化された波形情報をフーリエ変換してパワースペクトルが最大となる周期を求め、先にアクセント時刻と推定された時刻からこの周期時刻ごとに、その周期の時刻の近傍に他のアクセントが存在しているかどうかを検出する。そして、その周期の時刻の近傍に他のアクセントが存在している場合は、その時刻を拍時刻として推定する。また、これ以外の方法として、特開2008−107569号に示されるように、音響波形に対して全波整流を行った後、フィルタによって所定の周波数帯域の信号を抽出し、ピークの存在する時刻を拍時刻として推定する方法なども用いることができる。そして、このように拍時刻を推定した後、これらすべての拍についてBｊ（ｊ＝0,1,2,…）とラベリングする。
【００２３】
次に、このように楽曲の拍時刻を推定した後、この拍から小節における先頭拍を推定する。なお、ここでは楽曲として４／４拍子の楽曲を想定して説明するが、３／４拍子など他の拍子についても同様の処理を適用することができる。
【００２４】
先頭拍を推定する場合、まず、変動値算出手段３によって、各拍時刻における前後の微小時間における周波数の変動を算出する。このような処理を行うのは、各小節の先頭拍の時刻上では和音やメロディが大きく変化し、スペクトルが大きく変動する傾向が強いからである。ここでは、まず、先の拍時刻推定手段２で推定された拍Bｊ（ｊ＝0,1,2,…）について、周期的な４拍ごとに小節内のＩＤを割り当て、例えば、第k小節（k＝0,1,2,…）のm番目（m＝1,2,3,4）の拍をBm（k)とする。そして、Bm（k）と割り当てられた拍の前後のフレームからＦＦＴを用いてパワースペクトルを算出し、その前後におけるフレームのパワースペクトルの差分であるスペクトラムフラックスを算出する。なお、このスペクトラムフラックスが、本発明の変動値に相当する。
【００２５】
次に、時間変化算出手段５では、この変動値算出手段３によって算出された各拍の前後におけるスペクトルの変動値について各mの値毎に時間変化を算出する。すなわち、各小節における１拍目のスペクトラムフラックスの時間変化、各小節における２拍目のスペクトラムフラックスの時間変化、各小節における３拍目のスペクトラムフラックスの時間変化、各小節における４拍目のスペクトラムフラックスの時間変化を算出する。ここで、暫定的に決定されたk番目の小節におけるB1（k）、B2（k）、B3（k）、B4（k）のスペクトルの変動量を比較すると、多少の順位の入れ替わりはあるものの、各小節における先頭拍に対応するmのＩＤを持つ拍が常に最も大きいスペクトルの変動量を示すと考えられる。ここで、人間が主観的に感じる先頭拍（Downbeat）の周期をPSD（Periodicity of subjective downbeats）、音響波形上のスペクトルが大きく変動する周期をPLSV（Periodicity of the Large Spectral Variation）と呼ぶことにし、図２にPSDが変化しない場合、および、PSDが変化する場合の各Bm（k）におけるスペクトルの変動量の部分和（すなわち、各小節におけるk拍目の和）を示す。図２の上図は、拍時刻Bm（k）のスペクトルの変動量を算出し、mのＩＤ別の部分和をプロットしたグラフである。一方、図２の下図は、上図と同一の楽曲で、PSDが一時的に変化する場合の各Bm（k）におけるスペクトルの変動量からmのＩＤ別の部分和をプロットしグラフである。ここで、スペクトルの変動量の指標には、数１におけるスペクトラルフラックスを用いることとした。
【００２６】
【数１】
【００２７】
数１におけるNは、ＦＦＴに用いたサンプル点数で、8192点（すなわち、約185.8msec）を用いている。hは周波数解析を行う際のフレームのＩＤ、Lh[n]はh番目のフレームにおける周波数nのパワーを表している。図２の上図に示すように、特定の拍（図２ではm＝2が先頭拍となっている）におけるスペクトルの変動量が突出して増加しているのが確認できる。そして、図２のグラフの作成に用いた楽曲では、図２の上図で最大値を示しているm＝2と対応する拍B2（k）が先頭拍と一致していることを確認することができた。このように、PSDが変化しない楽曲に対しては、mのＩＤ別にスペクトルの変動量の合計値を求めることで先頭拍をある程度推定可能であることが分かる。ここで、図２の上図において、PSDが一時的に変化した場合について考察すると、図２の上図におけるm＝2の拍時刻においてスペクトルの変動量が減少し、先頭拍と対応する新たなmのＩＤを持つ拍時刻において、スペクトルの変動量が増加すると考えられる。図２の下図は、図２の上図で用いた楽曲と同一の楽曲に対して、本来存在しない拍の音響波形を意図的に挿入し、図２の上図と同様の方法でグラフを作成した結果である。意図的に拍の音響波形を挿入した時刻は、楽曲のちょうど半分の時刻（図２下図のChanging time）である。意図的に挿入された拍の音響波形が原因となり、PSDが変化し、本来２拍目であった先頭拍が、３拍目となっていることが確認できる。これは、意図的に挿入された拍の音響波形が原因で、PLSVが一時的に変化し、その時刻の前半部分までは、m＝2の拍時刻と実際の先頭拍が対応していたが、後半部分では、先頭拍と対応するmが3へと変化したためであると言える。このように、楽曲中においてPSDが一時的に変化した場合、該当する時刻を検出し、その時刻を基準とした区間ごとに先頭拍を推定する必要がある。
【００２８】
ところで、先頭拍以外の拍において、突発的に大きなスペクトルの変動が発生した場合、該当する拍におけるmのスペクトルの変動量の増加傾向に影響を与える可能性がある。そこで、本実施の形態では、この対策として、拍時刻Bm（k）におけるスペクトルの変動量Fm（k）へBm（k−1）およびBm（k+1）におけるスペクトルの変動量Fm（k−1）およびFm(k+1)を加算し、平均を求めることで、先頭拍以外の拍における突発的に発生する大きなスペクトルの変動の影響を低減することとした。また、この処理によって得られる拍時刻Bm（k）におけるスペクトルの変動量の平均値をAm（k）と表現する。Am（k）の算出式を数２に示す。
【００２９】
【数２】
【００３０】
ここで、Kは、暫定的に決定された小節における最終小節のＩＤを表している。
【００３１】
また、Am（k）は、フペクトラルフラックスに基づいて算出されているため、例えば、スネアドラムなどの全周波数帯域にパワーをもつ音がいくつかの小節にわたって発音されていた場合に、連続して非常に大きな値を示す場合がある。したがって、スペクトルの変動量のおおまかな増加傾向を、楽曲の局所的な物理的特徴に依存するAm(k)から直接求めることはあまり望ましくない。これは、PSDが変化せず、先頭拍と対応する拍において、スペクトルの変動量が常に最も大きい値を示す場合でも、スペクトルの変動量の増加傾向が一定になるとは限らないためである。この問題は、PSDが一時的に変化した時刻の特徴を得る際に悪影響を与え、最終的にはPSDが一時的に変化した時刻の誤検出の原因になると考えられる。この問題を解決するため、スペクトルの変動量の増加傾向を求める際に、Am（k）をそのまま用いるのではなく、暫定的に決定された各小節内のAm（k）の値を、スコアリング手段４によって順位別のスコア（4,3,2,1）に置き換え、その順序に基づいたスコアを付与することで、スペクトルの変動量の差を無視することとした。この順位スコアを先頭拍の可能性を表す指標と解釈し、「Downbeatness」と呼ぶとともに、拍時刻Bm（k）におけるDownbeatnessの値をDm（k）と表現する。
【００３２】
図３にDownbeatnessを算出する際の流れを示す。図３に示すように、暫定的に決定されたk番目の小節におけるA1（k）、A2（k）、A3（k）、A4（k）を比較し、最も大きい値を持つm番目の拍時刻Bm（k)から順に、4,3,2,1と順位スコアの値を付与していく。こうすることで、スペクトルの変動の急激な増加、あるいは減少の影響を低減することができる。
【００３３】
次に、変動位置推定手段６は、このように時系列に並べられたスペクトルの変動値に基づいて、先頭拍の変動した位置を推定する。
【００３４】
この変動位置推定手段６では、図２で用いた楽曲からDownbeatnessを算出し、mのＩＤ別の部分和を算出する。ここで得られた、各拍時刻B m（k）におけるmのＩＤ別の部分和の時間系列を、先頭拍らしさ関数と呼び、Sm（k）と表現する。先頭拍らしさ関数Sm（k）の算出式を数３に示す。
【００３５】
【数３】
【００３６】
図４に、図２で用いた楽曲から先頭拍らしさ関数を算出した結果を示す。図４に示すように、Downbeatnessは、各mにおけるスペクトルの変動量の増加傾向を良く表していることが確認できる。特に、図４の下図における破線の円内において、S3（k）が急激に増加していることが確認できる。そして、Sm（k）の局所的な傾きが急激に変化した地点は、Sm（k）の回帰直線とSm（k）の距離が遠い地点と見なすことができる。Sm（k）の回帰直線をCm（k）とすると、暫定的に決定されたk番目の小節におけるm番目のSm（k）と回帰直線Cm（k）との距離Rm（k）は数４で求めることができる。
【００３７】
【数４】
【００３８】
したがって、Sm（k）の局所的な傾きが急激に変化した地点、すなわちPSDが一時的に変化した可能性が高い時刻は、Rm（k）のピークから得られると期待できる。このとき、PSDが一時的に変化した時刻においては、特定のmのみではなく、全mにおけるSm（k）において、スペクトルの変動量の増加傾向が変化する可能性が高い。したがって、PSDが一時的に変化した時刻を得るためにRm（k）のピークを求める場合は、m＝1,2,3,4についてそれぞれRm（k）を加算し、ピークとなる時刻を求める処理が妥当であると言える。
【００３９】
図５に、図４のPSDが変化する場合の楽曲からRm（k）を抽出した結果を示す。図５に示すように、m＝1,2,3,4のRm（k）を加算した時間関数（図５の細い破線）における中間部分のピークと、楽曲のPSDが変化した時刻が概ね一致していることが確認できる。
【００４０】
ここでは、Rm（k）を加算した時間関数を先頭拍差異検出関数と呼ぶこととし、Rsum（k）と表現する。
【００４１】
図５から確認できるように、Rsum（k）は、PSDが変化した時刻と良い対応を持つ。そこで、楽曲中の拍時刻に基づいて、このRsum（k）を算出し、その最大値であるピークを検出することで、PSDが変化した時刻を検出し、この時刻を先頭拍の変化した時刻と推定する。そして、このように先頭拍の変化した時刻を推定できた場合、その拍から新たに先頭拍の推定処理を同様に行っていくようにする。
【００４２】
次に、このように構成された先頭拍推定システム１００を用いて楽曲の先頭拍を推定する方法について、図６のフローチャートを用いて説明する。
【００４３】
まず、先頭拍を推定する楽曲の音響波形の入力を受け付け（ステップＳ１）、これを記憶部に格納しておく。
【００４４】
次に、その受け付けた音響波形から拍となる位置の時刻を推定する（ステップＳ２）。この受け付けた音響波形から拍時刻を推定する場合は、入力された音響波形を平滑化し、所定の閾値以下の微小なパワーを除いたおおまかなアクセント時刻を移動平均によって求める。一方、これとは別に、平滑化された波形情報をフーリエ変換してパワースペクトルが最大となる周期を求め、先におおまかに推定されたアクセント時刻からこの周期時刻ごとに、その周期時刻の近傍に他のアクセントが存在しているかどうかを検出する。そして、その周期時刻の近傍に他のアクセントが存在している場合は、その時刻を拍時刻として推定する（ステップＳ２）。
【００４５】
次に、このように推定されたすべての拍についてBｊ(ｊ＝0,1,2,…)とラベリングするとともに、周期的な拍（ここでは４拍）ごとにＩＤ（Bm(k)）を割り当てるラベリング処理を行う（ステップＳ３）。
【００４６】
そして、このようにBm（k）とラベリングされたすべての拍について、その拍時刻の前後のフレームからＦＦＴを用いてパワースペクトルを算出し、そのパワースペクトルの差分であるスペクトルフラックス（すなわち、スペクトルの変動値）を算出する（ステップＳ４）。
【００４７】
次に、このスペクトルフラックスを各ＩＤごとに時系列に並べるようにするが、ここでは、突発的なスペクトルの変動を無視できるように、スコアリング手段４によって周期的に割り当てられたＩＤの範囲内（すなわち、k＝1,2,3,4の範囲内）でスペクトルフラックスの大きいものから順に4,3,2,1とスコアを与える（ステップＳ５）。そして、このように与えられたスコアを各ＩＤごとに時系列に加算して先頭拍らしさ関数Sm（k）を算出し、図４に対応する数値列を取得する（ステップＳ６）。
【００４８】
そして、この先頭拍らしさ関数Sm（k）から、スペクトルの変動量が急激に変化した位置を算出すべく、Sm（k）に対して単一の回帰直線Cm（k）を算出し（ステップＳ７）、その回帰直線Rm（k）との距離を数４により求めて、その値を加算した先頭拍差異検出関数を取得する。そして、その最大値であるピークを抽出して、その時刻を先頭拍の変化した時刻と推定する（ステップＳ８）。
【００４９】
このように上記実施の形態によれば、楽曲の入力を受け付ける楽曲受付手段１と、当該受け付けた楽曲から拍の位置を推定する拍時刻推定手段２と、当該推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出する変動値算出手段３と、すべての拍について一定周期の拍ごとに前記変動値の時間変化を算出する時間変化算出手段５と、当該時間変化が最大となった位置を先頭拍の変動した位置として推定する変動位置推定手段６とを備えるようにしたので、先頭拍の変化した位置を推定することができるようになる。これにより、その変化した位置から新たに先頭拍を推定すれば、楽曲の途中で小節内の拍数が変化した場合であっても、すべての小節の先頭拍を推定することができるようになる。
【００５０】
また、周期的に付与されたＩＤ（k＝1,2,3,4）の範囲内で変動値を順位付けし、その順位付けされた値に基づいて時間変化を算出するようにしたので、先頭拍以外の拍における突発的に発生する大きなスペクトルの変動の影響を低減することができるようになる。
【００５１】
さらには、時間軸方向に加算された値に対して、回帰直線に関する値を算出し、その回帰直線に関する値と時間軸方向に加算された値との差分を算出して、その差分が最大となった位置を先頭拍の変動した位置と推定するようにしたので、これによって先頭拍の変動した位置を容易に推定することができるようになる。
【００５２】
なお、本発明は上記実施の形態に限定されることなく種々の態様で実施することができる。
【００５３】
例えば、上記実施の形態では、スペクトルフラックスの時間変化を算出する場合、スコアリング処理を行った後に時間変化を算出するようにしたが、スペクトルフラックスの状態で時間変化を算出するようにしてもよい。
【００５４】
また、上記実施の形態では、時系列にスペクトルの変動値を加算して回帰直線との差を求めるようにしたが、時系列にスペクトルの変動値を加算するのではなく、それぞれの変動値そのものを時系列に並べ、スコアリングされた各ＩＤの時系列の値との大小関係をみることによって変動位置を推定するようにしてもよい。
【００５５】
さらには、上記実施の形態では、周期的なＩＤとして４／４拍子に対応した４拍ごとにＩＤを付与するようにしたが、３／４拍子の場合は、３拍ごとにＩＤを付与するようにしてもよい。
【００５６】
また、上記実施の形態では、回帰直線との距離Rm（k）が最大となる位置を変動位置と推定するようにしたが、複数回変動位置が生ずる可能性もあるため、所定の閾値を設定しておき、その閾値よりもRm（k）の値が大きい場合は、先頭拍が変動した位置と推定するようにしてもよい。
【実施例１】
【００５７】
上記手法によって、推定された先頭拍が正しいかどうかを評価した。上記手法の評価基準の一つとして、PSDが変化しない場合の楽曲において、提案手法と既存手法とを比較した際の推定精度への影響が挙げられる。これは、PSDが変化しない楽曲において、提案手法を用いて先頭拍を推定した場合の誤推定が起こる割合を評価することができる。この評価を行うためには、楽曲に対して人手で拍時刻および先頭拍の時刻のラベリングデータを作成し、ラベリングされた拍時刻に基づいて、既存手法および提案手法によって先頭拍を推定し、正解率を比較すれば良いと言える。以降この評価基準の条件を、S（Simple） 条件と呼ぶ。
【００５８】
もう一つの評価基準として、PSDが変化する場合の楽曲および推定された拍時刻上においてPLSVが変化する場合の楽曲において、提案手法と既存手法とを比較した際の推定精度の違いが挙げられる。すなわち、既存手法では推定が困難な楽曲に対する提案手法の有効性を評価することができる。
【００５９】
PSDが変化する状況としては、楽曲の小節の長さが一時的に変化する状況が考えられる。また、推定された拍時刻上においてPLSVが変化する状況としては、拍時刻の誤推定が発生した状況が挙げられる。これらの状況下において、既存手法および提案手法によって先頭拍を推定し、正解率を比較すれば良いと言える。本実施例では、これらの状況を意図的に再現することで、既存手法では推定が困難な楽曲に対する提案手法の有効性を評価することとした。具体的には、楽曲に対して人手で拍時刻および先頭拍の時刻のラベリングデータを作成し、ある単一の小節に対して、本来存在しない拍の音響波形を挿入することで、小節の長さが一時的に変化する状況を再現することとした。さらに、ラベリングされた拍時刻に対して、意図的に誤った拍時刻をラベリングすることで、拍時刻の誤推定が発生した状況を再現することとした。以降この評価基準の条件を、C（Challenge）条件と呼び、特に小節の長さが一時的に変化する条件をC m （measure）条件、拍時刻の誤推定が発生する条件をC b （beat）条件と呼ぶこととする。
【００６０】
上記各条件で評価を行うために、まず聴取者にRWC データベースから抜粋した楽曲を聴取させ、一定の間隔で発音される正弦波の時間間隔を拍時刻および先頭拍と同期するようにそれぞれ調整させた。そして、楽曲の拍時刻および先頭拍の時刻において発音される正弦波の時間系列データのwaveファイルをそれぞれ作成させ、それらのデータを正解データとした。
【００６１】
この正解データと、楽曲から推定された先頭拍の時刻を比較する必要がある。そこで、推定された先頭拍の時刻で正弦波が発音されるよう合成音を作成し、正解データとのパワーの比率を正解率Qとして求めることで評価を行うこととした。この正解率Qの算出方法を数５に示す。
【００６２】
【数５】
【００６３】
ここで、giは正解データの音響信号におけるパワー関数、eiは推定された先頭拍の時刻で正弦波が発音される音響信号におけるパワー関数を表している。
【００６４】
表1に、評価結果を示す。表1におけるEは既存手法、Pは提案手法を表している。表1のS条件において、提案手法と比較して、既存手法の方が僅かに高い推定精度であることについて言及すると、評価に用いた楽曲がPSDの変化しない楽曲であったためであると考えられる。次に、表1のC条件に示す通り、楽曲の中盤でPSDが変化したり（Cm条件）、拍の誤検出が発生すると（Cｂ条件）、既存手法の推定率が大きく低下しているのが確認できる。それに対して、提案手法を用いた場合、PSDが変化しない場合と同等の推定精度で先頭拍を推定できていることが確認できる（表1内の下線部）。したがって、PSDが変化する楽曲、あるいは拍時刻の誤推定の発生により、推定された拍時刻上においてPLSVが変化するような楽曲に対して、提案手法が有効であることを確認できた。このため、提案手法を用いることで、小節の長さが変化する楽曲や、本来拍が存在しない時刻に拍時刻が挿入されている場合でも、それらの影響を受けずに先頭拍を推定できることが確認できた。
【００６５】
【表１】
【符号の説明】
【００６６】
１・・・楽曲受付手段
２・・・拍時刻推定手段
３・・・変動値算出手段
４・・・スコアリング手段
５・・・時間変化算出手段
６・・・変動位置推定手段
１００・・・先頭拍推定システム
【特許請求の範囲】
【請求項１】
楽曲の入力を受け付けるステップと、
当該受け付けた楽曲から拍の位置を推定するステップと、
当該推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出するステップと、
すべての拍について一定周期の拍ごとに前記変動値の時間変化を算出するステップと、
当該時間変化の変化量が最大となった位置、もしくは、時間変化の変化量が所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するステップと、
を備えたことを特徴とする小節における先頭拍推定方法。
【請求項２】
前記時間変化を算出するステップが、
前記推定された拍に対して、周期的な識別符号を付与し、当該識別符号の付与されたすべての拍における前後の分析フレームを周波数解析してスペクトルを算出し、当該算出されたスペクトルの変動値を算出するとともに、各識別符号が付与された拍について、前記変動値を楽曲の時間軸方向に加算するようにしたものである請求項１に記載の先頭拍推定方法。
【請求項３】
前記変動した位置と推定するステップが、
時間軸方向に加算された値に対して回帰直線に関する値を算出し、
当該算出された回帰直線に関する値と前記時間軸方向に加算された値との差分を算出し、
当該算出された差分の値が最大となった位置、もしくは、所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するようにしたものである請求項１に記載の先頭拍推定方法。
【請求項４】
楽曲の入力を受け付ける楽曲受付手段と、
当該受け付けた楽曲から拍の位置を推定する拍時刻推定手段と、
当該推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出する変動値算出手段と、
すべての拍について一定周期の拍ごとに前記変動値の時間変化を算出する時間変化算出手段と、
当該時間変化が最大となった位置、もしくは、時間変化が所定の閾値よりも大きくなった位置を、先頭拍の変動した位置として推定する変動位置推定手段と、
を備えたことを特徴とする先頭拍推定システム。
【請求項５】
前記時間変化算出手段が、
当該推定された拍に対して、周期的な識別符号を付与し、当該識別符号の付与されたすべての拍における前後の分析フレームを周波数解析してスペクトルを算出し、当該算出されたスペクトルの変動値を算出するとともに、各識別符号が付与された拍について、前記変動値を楽曲の時間軸方向に加算するように構成されたものである請求項４に記載の先頭拍推定システム。
【請求項６】
前記変動位置推定手段が、
時間軸方向に加算された値に対して回帰直線に関する値を算出し、
当該算出された回帰直線に関する値と前記時間軸方向に加算された値との差分を算出し、
当該算出された差分の値が最大となった位置、もしくは、所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するように構成されたものである請求項４に記載の先頭拍推定システム。
【請求項１】
楽曲の入力を受け付けるステップと、
当該受け付けた楽曲から拍の位置を推定するステップと、
当該推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出するステップと、
すべての拍について一定周期の拍ごとに前記変動値の時間変化を算出するステップと、
当該時間変化の変化量が最大となった位置、もしくは、時間変化の変化量が所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するステップと、
を備えたことを特徴とする小節における先頭拍推定方法。
【請求項２】
前記時間変化を算出するステップが、
前記推定された拍に対して、周期的な識別符号を付与し、当該識別符号の付与されたすべての拍における前後の分析フレームを周波数解析してスペクトルを算出し、当該算出されたスペクトルの変動値を算出するとともに、各識別符号が付与された拍について、前記変動値を楽曲の時間軸方向に加算するようにしたものである請求項１に記載の先頭拍推定方法。
【請求項３】
前記変動した位置と推定するステップが、
時間軸方向に加算された値に対して回帰直線に関する値を算出し、
当該算出された回帰直線に関する値と前記時間軸方向に加算された値との差分を算出し、
当該算出された差分の値が最大となった位置、もしくは、所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するようにしたものである請求項１に記載の先頭拍推定方法。
【請求項４】
楽曲の入力を受け付ける楽曲受付手段と、
当該受け付けた楽曲から拍の位置を推定する拍時刻推定手段と、
当該推定された拍における前後の分析フレームからスペクトルの差分である変動値を算出する変動値算出手段と、
すべての拍について一定周期の拍ごとに前記変動値の時間変化を算出する時間変化算出手段と、
当該時間変化が最大となった位置、もしくは、時間変化が所定の閾値よりも大きくなった位置を、先頭拍の変動した位置として推定する変動位置推定手段と、
を備えたことを特徴とする先頭拍推定システム。
【請求項５】
前記時間変化算出手段が、
当該推定された拍に対して、周期的な識別符号を付与し、当該識別符号の付与されたすべての拍における前後の分析フレームを周波数解析してスペクトルを算出し、当該算出されたスペクトルの変動値を算出するとともに、各識別符号が付与された拍について、前記変動値を楽曲の時間軸方向に加算するように構成されたものである請求項４に記載の先頭拍推定システム。
【請求項６】
前記変動位置推定手段が、
時間軸方向に加算された値に対して回帰直線に関する値を算出し、
当該算出された回帰直線に関する値と前記時間軸方向に加算された値との差分を算出し、
当該算出された差分の値が最大となった位置、もしくは、所定の閾値よりも大きくなった位置を先頭拍の変動した位置と推定するように構成されたものである請求項４に記載の先頭拍推定システム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図２】
【図３】
【図４】
【図５】
【図６】
【公開番号】特開２０１３−２０００７（Ｐ２０１３−２０００７Ａ）
【公開日】平成２５年１月３１日（２０１３．１．３１）
【国際特許分類】
【出願番号】特願２０１１−１５１９８４（Ｐ２０１１−１５１９８４）
【出願日】平成２３年７月８日（２０１１．７．８）
【出願人】（５９７０６５３２９）学校法人 龍谷大学 (120)
【Ｆターム（参考）】
【公開日】平成２５年１月３１日（２０１３．１．３１）
【国際特許分類】
【出願日】平成２３年７月８日（２０１１．７．８）
【出願人】（５９７０６５３２９）学校法人 龍谷大学 (120)
【Ｆターム（参考）】
[ Back to top ]