本発明の一実施形態に係る音響信号分析装置10について説明する。音響信号分析装置10は、以下説明するように、楽曲を表わす音響信号を入力して、その楽曲における拍点及びテンポの推移を検出する。音響信号分析装置10は、図1に示すように、入力操作子11、コンピュータ部12、表示器13、記憶装置14、外部インターフェース回路15及びサウンドシステム16を備えており、これらがバスBSを介して接続されている。
入力操作子11は、オン・オフ操作に対応したスイッチ(例えば数値を入力するためのテンキー)、回転操作に対応したボリューム又はロータリーエンコーダ、スライド操作に対応したボリューム又はリニアエンコーダ、マウス、タッチパネルなどから構成される。これらの操作子は、演奏者の手によって操作されて、分析対象の楽曲の選択、音響信号の分析開始又は停止、楽曲の再生又は停止(後述するサウンドシステム16からの出力又は停止)、音響信号の分析に関する各種パラメータの設定などに用いられる。入力操作子11を操作すると、その操作内容を表す操作情報が、バスBSを介して、後述するコンピュータ部12に供給される。
コンピュータ部12は、バスBSにそれぞれ接続されたCPU12a、ROM12b及びRAM12cからなる。CPU12aは、詳しくは後述する音響信号分析プログラム及びそのサブルーチンをROM12bから読み出して実行する。ROM12bには、音響信号分析プログラム及びそのサブルーチンに加えて、初期設定パラメータ、表示器13に表示される画像を表わす表示データを生成するための図形データ及び文字データなどの各種データが記憶されている。RAM12cには、音響信号分析プログラムの実行時に必要なデータが一時的に記憶される。
表示器13は、液晶ディスプレイ(LCD)によって構成される。コンピュータ部12は、図形データ、文字データなどを用いて表示すべき内容を表わす表示データを生成して表示器13に供給する。表示器13は、コンピュータ部12から供給された表示データに基づいて画像を表示する。例えば分析対象の楽曲の選択時には、楽曲のタイトルリストが表示される。また、例えば分析終了時には、拍点及びテンポの推移を表わす拍・テンポ情報リスト及びそのグラフ(図20乃至図23参照)が表示される。
また、記憶装置14は、HDD、FDD、CD−ROM、MO、DVDなどの大容量の不揮発性記録媒体と、同各記録媒体に対応するドライブユニットから構成されている。記憶装置14には、複数の楽曲をそれぞれ表わす複数の楽曲データが記憶されている。楽曲データは、楽曲を所定のサンプリング周期(例えば44.1kHz)でサンプリングして得られた複数のサンプル値からなり、各サンプル値が記憶装置14における連続するアドレスに順に記録されている。楽曲のタイトルを表わすタイトル情報、楽曲データの容量を表わすデータサイズ情報なども楽曲データに含まれている。楽曲データは予め記憶装置14に記憶されていてもよいし、後述する外部インターフェース回路15を介して外部から取り込んでもよい。記憶装置14に記憶されている楽曲データは、CPU12aによって読み込まれ、楽曲における拍点及びテンポの推移が分析される。
外部インターフェース回路15は、音響信号分析装置10を電子音楽装置、パーソナルコンピュータなどの外部機器に接続可能とする接続端子を備えている。音響信号分析装置10は、外部インターフェース回路15を介して、LAN(Local Area Network)、インターネットなどの通信ネットワークにも接続可能である。
サウンドシステム16は、楽曲データをアナログ音信号に変換するD/A変換器、変換したアナログ音信号を増幅するアンプ、及び増幅されたアナログ音信号を音響信号に変換して出力する左右一対のスピーカを備えている。ユーザが入力操作子11を用いて分析対象の楽曲の再生を指示すると、CPU12aは、分析対象の楽曲データをサウンドシステム16に供給する。これにより、ユーザは分析対象の楽曲を試聴できる。
つぎに、上記のように構成した音響信号分析装置10の動作について説明する。まず、その概略について説明する。分析対象の楽曲は複数のフレームti{i=0,1,・・・,last}に分割される。そして、拍の存在に関する特徴を表すオンセット特徴量XO及びテンポに関する特徴を表すBPM特徴量XBがフレームtiごとに計算される。各フレームtiにおける拍周期bの値(テンポの逆数に比例する値)及び次の拍までのフレーム数nの値の組み合わせに応じて分類された状態qb,nの系列として記述された確率モデル(隠れマルコフモデル)のうち、観測値としてのオンセット特徴量XO及びBPM特徴量XBが同時に観測される確率を表わす観測尤度の系列が最も尤もらしい確率モデルが選択される(図2参照)。これにより、分析対象の楽曲における拍点及びテンポの推移が検出される。なお、拍周期bは、フレームの数によって表わされる。したがって、拍周期bの値は「1≦b≦bmax」を満たす整数であり、拍周期bの値が「β」である状態では、フレーム数nの値は「0≦n<β」を満たす整数である。
つぎに、音響信号分析装置10の動作について具体的に説明する。ユーザが音響信号分析装置10の図示しない電源スイッチをオンにすると、CPU12aは、図3に示す音響信号分析プログラムをROM12bから読み出して実行する。
CPU12aは、ステップS10にて音響信号分析処理を開始し、ステップS11にて、記憶装置14に記憶されている複数の楽曲データにそれぞれ含まれるタイトル情報を読み込んで、楽曲のタイトルをリスト形式で表示器13に表示する。ユーザは、入力操作子11を用いて、表示器13に表示された楽曲の中から分析対象の楽曲データを選択する。なお、ステップS11にて分析対象の楽曲データを選択する際、選択しようとする楽曲データが表す楽曲の一部又は全部を再生して楽曲データの内容を確認できるように構成してもよい。
つぎに、CPU12aは、ステップS12にて、音響信号分析のための初期設定を実行する。具体的には、前記選択された楽曲データのデータサイズ情報に応じた記憶領域をRAM12c内に確保し、前記確保した記憶領域に前記選択された楽曲データを読み込む。また、分析結果を表す拍・テンポ情報リスト、オンセット特徴量XO、BPM特徴量XBなどを一時的に記憶する領域をRAM12c内に確保する。
詳しくは後述するが、本プログラムによる分析結果は、記憶装置14に保存される(ステップS21)。前記選択された楽曲が本プログラムによって過去に分析されたことがあれば、記憶装置14にその分析結果が保存されている。そこで、CPU12aは、ステップS13にて、前記選択された楽曲の分析に関する既存のデータ(以下、単に既存データと呼ぶ)を検索する。既存データが有れば、CPU12aは、ステップS14にて「Yes」と判定して、ステップS15にて既存データをRAM12cに読み込み、後述するステップS19に処理を進める。一方、既存データが無ければ、CPU12aは、ステップS14にて「No」と判定して、その処理をステップS16に進める。
CPU12aは、ステップS16にて、図4に示す特徴量計算プログラムをROM12bから読み出して実行する。特徴量計算プログラムは、音響信号分析プログラムのサブルーチンである。
CPU12aは、ステップS161にて特徴量計算処理を開始する。つぎに、CPU12aは、ステップS162にて、図5に示すように、前記選択された楽曲を所定の時間間隔をおいて区切り、複数のフレームti{i=0,1,・・・,last}に分割する。各フレームの長さは共通である。説明を簡単にするために、本実施形態では各フレームの長さを125msとする。上記のように、各楽曲のサンプリング周波数は44.1kHzであるので、各フレームは、約5000個のサンプル値から構成されている。そして、以下説明するように、フレームごとに、オンセット特徴量XO及びBPM(beats per minute(1分間あたりの拍数))特徴量XBを計算する。
つぎに、CPU12aは、ステップS163にて、フレームごとに短時間フーリエ変換を実行して、図6に示すように、各周波数ビンfj{j=1,2・・・}の振幅A(fj,ti)を計算する。そして、CPU12aは、ステップS164にて、振幅A(f1,ti),A(f2,ti)・・・を周波数ビンfjごとに設けられたフィルタバンクFBOjによってフィルタ処理することにより、所定の周波数帯域wk{k=1,2,・・・}の振幅M(wk,ti)を計算する。周波数ビンfj用のフィルタバンクFBOjは、図7に示すように、通過帯域の中心周波数が互いに異なる複数のバンドパスフィルタBPF(wk,fj)からなる。フィルタバンクFBOjを構成する各バンドパスフィルタBPF(wk,fj)の中心周波数は、対数周波数軸上で等間隔であり、かつ各バンドパスフィルタBPF(wk,fj)の通過帯域幅は、対数周波数軸上で共通である。各バンドパスフィルタBPF(wk,fj)は、通過帯域の中心周波数から通過帯域の下限周波数側及び上限周波数側へ向かうに従って徐々にゲインがそれぞれ小さくなるように構成されている。CPU12aは、図4のステップS164に示すように、周波数ビンfjごとに振幅A(fj,ti)とバンドパスフィルタBPF(wk,fj)のゲインとを積算する。そして、前記周波数ビンfjごとに計算した積算結果を全ての周波数ビンfjについて合算して振幅M(wk,ti)とする。上記のようにして計算された振幅Mの系列を図8に例示する。
つぎに、CPU12aは、ステップS165にて、振幅Mの時間変化に基づいてフレームtiのオンセット特徴量XO(ti)を計算する。具体的には、図4のステップS165に示すように、周波数帯域wkごとに、フレームti−1からフレームtiへの振幅Mの増加量R(wk,ti)を計算する。ただし、フレームti―1の振幅M(wk,ti−1)とフレームtiの振幅M(wk,ti)とが同じである場合、又はフレームtiの振幅M(wk,ti)がフレームti―1の振幅M(wk,ti−1)よりも小さい場合は、増加量R(wk,ti)は「0」とする。そして、周波数帯域wkごとに計算した増加量R(wk,ti)を全ての周波数帯域w1,w2,・・・について合算してオンセット特徴量XO(ti)とする。上記のようにして計算されたオンセット特徴量XOの系列を図9に例示する。一般に、楽曲においては、拍が存在する部分の音量が大きい。したがって、オンセット特徴量XO(ti)が大きいほど、フレームtiに拍が存在する可能性が高い。
つぎに、CPU12aは、オンセット特徴量XO(t0),XO(t1)・・・を用いて、BPM特徴量XBをフレームtiごとに計算する。なお、フレームtiのBPM特徴量XB(ti)は、拍周期bごとに計算されたBPM特徴量XBb=1,2・・・(ti)の集合として表わされる(図11参照)。まず、CPU12aは、ステップS166にて、オンセット特徴量XO(t0),XO(t1)・・・をこの順にフィルタバンクFBBに入力してフィルタ処理する。フィルタバンクFBBは、拍周期bの値に応じてそれぞれ設けられた複数のコムフィルタDbからなる。コムフィルタDb=βは、フレームtiのオンセット特徴量XO(ti)を入力すると、前記入力したオンセット特徴量XO(ti)と「β」だけ先行するフレームti−βのオンセット特徴量XO(ti−β)に対する出力としてのデータXDb=β(ti−β)とを所定の比率で加算してフレームtiのデータXDb=β(ti)として出力する(図10参照)。すなわち、コムフィルタDb=βは、データXDb=βをフレーム数βに相当する時間だけ保持する保持手段としての遅延回路db=βを有する。上記のようにして、オンセット特徴量XOの系列XO(t){=XO(t0),XO(t1)・・・}をフィルタバンクFBBに入力することにより、データXDbの系列XDb(t){=XDb(t0),XDb(t1)・・・}が計算される。
つぎに、CPU12aは、ステップS167にて、データXDbの系列XDb(t)を時系列的に逆にしたデータ列をフィルタバンクFBBに入力することにより、BPM特徴量の系列XBb(t){=XBb(t0),XBb(t1)・・・}が得られる。これにより、オンセット特徴量XO(t0),XO(t1)・・・の位相とBPM特徴量XBb(t0),XBb(t1)・・・の位相のずれを「0」にすることができる。上記のようにして計算されたBPM特徴量XB(ti)を図11に例示する。上記のように、BPM特徴量XBb(ti)は、オンセット特徴量XO(ti)と拍周期bの値に相当する時間(すなわち、フレーム数b)だけ遅延させたBPM特徴量XBb(ti―b)とを所定の比率で加算して計算されるので、オンセット特徴量XO(t0),XO(t1)・・・が拍周期bの値に相当する時間間隔をおいてピークを有する場合、BPM特徴量XBb(ti)の値が大きくなる。楽曲のテンポは、1分間あたりの拍数で表されるから、拍周期bは1分間あたりの拍数の逆数に比例する。例えば、図11に示す例では、拍周期bの値が「4」であるときのBPM特徴量XBbの値(BPM特徴量XBb=4)が最も大きい。したがって、この例では拍が4フレームおきに存在する可能性が高い。本実施形態では、1フレームの時間の長さを125msとしたので、この場合の拍の間隔は0.5sである。すなわち、テンポは120BPM(=60s/0.5s)である。
つぎに、CPU12aは、ステップS168にて、特徴量計算処理を終了し、その処理を音響信号分析処理(メインルーチン)のステップS17に進める。
CPU12aは、ステップS17にて、図12に示す対数観測尤度計算プログラムをROM12bから読み出して実行する。対数観測尤度計算プログラムは、音響信号分析プログラムのサブルーチンである。
CPU12aは、ステップS171にて対数観測尤度計算処理を開始する。そして、以下説明するように、オンセット特徴量XO(ti)の尤度P(XO(ti)|Zb,n(ti))、及びBPM特徴量XB(ti)の尤度P(XB(ti)|Zb,n(ti))を計算する。なお、上記の「Zb=β,n=η(ti)」は、フレームtiにおいて、拍周期bの値が「β」であり、且つ次の拍までのフレーム数nの値が「η」である状態qb=β,n=ηのみが生起していることを表わす。フレームtiにおいて状態qb=β,n=ηと状態qb≠β,n≠ηとが同時に生起することはない。したがって、尤度P(XO(ti)|Zb=β,n=η(ti))は、フレームtiにおいて、拍周期bの値が「β」であり、且つ次の拍までのフレーム数nの値が「η」であるという条件のもとでオンセット特徴量XO(ti)が観測される確率を表わす。また、尤度P(XB(ti)|Zb=β,n=η(ti))は、フレームtiにおいて、拍周期bの値が「β」であり、且つ次の拍までのフレーム数nの値が「η」であるという条件のもとでBPM特徴量XB(ti)が観測される確率を表わす。
まず、CPU12aは、ステップS172にて、尤度P(XO(ti)|Zb,n(ti))を計算する。次の拍までのフレーム数nの値が「0」であるとき、オンセット特徴量XOは、平均値が「3」であって、且つ分散が「1」である第1の正規分布に従って分布するものとする。すなわち、第1の正規分布の確率変数としてオンセット特徴量XO(ti)を代入した値を尤度P(XO(ti)|Zb,n=0(ti))として計算する。また、拍周期bの値が「β」であり、次の拍までのフレーム数nの値が「β/2」であるとき、オンセット特徴量XOは、平均値が「1」であって、且つ分散が「1」である第2の正規分布に従って分布するものとする。すなわち、第2の正規分布の確率変数としてオンセット特徴量XO(ti)を代入した値を尤度P(XO(ti)|Zb=β,n=β/2(ti))として計算する。また、次の拍までのフレーム数nの値が「0」及び「β/2」のうちのいずれの値とも異なるとき、オンセット特徴量XOは、平均値が「0」であって、且つ分散が「1」である第3の正規分布に従って分布するものとする。すなわち、第3の正規分布の確率変数としてオンセット特徴量XO(ti)を代入した値を尤度P(XO(ti)|Zb,n≠0,β/2(ti))として計算する。
オンセット特徴量XOの系列が{10,2,0.5,5,1,0,3,4,2}であるときの尤度P(XO(ti)|Zb=6,n(ti))の対数を計算した結果を、図13に例示する。同図に示すように、オンセット特徴量XOの値が大きいフレームtiほど、尤度P(XO(ti)|Zb,n=0(ti))が尤度P(XO(ti)|Zb,n≠0(ti))に比べて大きい。このように、オンセット特徴量XOの値が大きいフレームtiほど、フレーム数nの値が「0」であるときに拍が存在する可能性が高くなるように、確率モデル(第1乃至第3の正規分布、及びそれらのパラメータ(平均値及び分散))が設定されている。なお、第1乃至第3の正規分布のパラメータの値は、上記実施形態に限られない。これらのパラメータの値は、実験を繰り返して決定してもよいし、機械学習を用いて決定してもよい。なお、この例では、オンセット特徴量XOの尤度Pを計算するための確率分布関数として正規分布を用いているが、確率分布関数として他の関数(例えば、ガンマ分布、ポアソン分布など)を用いても良い。
つぎに、CPU12aは、ステップS173にて、尤度P(XB(ti)|Zb,n(ti))を計算する。尤度P(XB(ti)|Zb=γ,n(ti))は、図14に示すテンプレートTPγ{γ=1,2・・・}に対するBPM特徴量XB(ti)の適合度に相当する。具体的には、尤度P(XB(ti)|Zb=γ,n(ti))は、BPM特徴量XB(ti)とテンプレートTPγ{γ=1,2・・・}との内積に相当する(図12のステップS173の演算式を参照)。なお、この演算式におけるκbは、オンセット特徴量XOに対するBPM特徴量XBの重みを決定する係数である。つまり、κbを大きく設定するほど、結果的に、後述する拍・テンポ同時推定処理においてBPM特徴量XBが重視される。また、この演算式におけるZ(κb)は、κbに依存する正規化係数である。テンプレートTPγは、図14に示すように、BPM特徴量XB(ti)を構成するBPM特徴量XBb(ti)にそれぞれ乗算される係数δγ,bからなる。テンプレートTPγは、その係数δγ,γが最大であり、係数δγ,2γ,係数δγ,3γ・・・,係数δγ,(「γ」の整数倍),・・・がそれぞれ極大となるように設定されている。すなわち、例えば、テンプレートTPγ=2は、2フレームおきに拍が存在する楽曲に適合するように構成されている。なお、この例では、BPM特徴量XBの尤度Pを計算するためにテンプレートTPを用いているが、これに代えて確率分布関数(例えば、多項分布、ディリクレ分布、多次元正規分布、多次元ポアソン分布など)を用いても良い。
BPM特徴量XB(ti)が図11に示すような値であった場合に、図14に示すテンプレートTPγ{γ=1,2・・・}を用いて尤度P(XB(ti)|Zb,n(ti))を計算し、その対数を計算した結果を図15に例示する。この例では、尤度P(XB(ti)|Zb=4,n(ti))が最も大きいので、BPM特徴量XB(ti)は、テンプレートTP4に最も適合している。
つぎに、CPU12aは、ステップS174にて、尤度P(XO(ti)|Zb,n(ti))の対数と尤度P(XB(ti)|Zb,n(ti))の対数をそれぞれ加算し、その結果を対数観測尤度Lb,n(ti)とする。なお、尤度P(XO(ti)|Zb,n(ti))と尤度P(XB(ti)|Zb,n(ti))とを積算した結果の対数を対数観測尤度Lb,n(ti)としても同じ結果が得られる。つぎに、CPU12aは、ステップS175にて、対数観測尤度計算処理を終了し、その処理を音響信号分析処理(メインルーチン)のステップS18に進める。
つぎに、CPU12aは、ステップS18にて、図16に示す拍・テンポ同時推定プログラムをROM12bから読み出して実行する。拍・テンポ同時推定プログラムは、音響信号分析プログラムのサブルーチンである。この拍・テンポ同時推定プログラムは、ビタビアルゴリズムを用いて最尤の状態の系列Qを計算するプログラムである。ここで、その概略について説明する。CPU12aは、まず、フレームt0からフレームtiまでオンセット特徴量XO及びBPM特徴量XBを観測したときにフレームtiの状態qb,nの尤度が最大となるような状態の系列を選択した場合の状態qb,nの尤度を尤度Cb,n(ti)とするとともに、各状態qb,nに遷移する1つ前のフレームの状態(遷移元の状態)を状態Ib,n(ti)として記憶する。つまり、遷移後の状態が状態qb=βe,n=ηeであって、遷移元の状態が状態qb=βs,n=ηsであるとき、状態Ib=βe,n=ηe(ti)は、状態qb=βs,n=ηsである。CPU12aは、上記のようにして尤度C及び状態Iをフレームtlastまで計算し、その結果を用いて最尤の状態の系列Qを選択する。
なお、後述する具体例では、その説明を簡単にするために、分析対象の楽曲の拍周期bの値が「3」、「4」及び「5」のうちのいずれかであるものとする。すなわち、対数観測尤度Lb,n(ti)が図17に例示するように計算された場合の拍・テンポ同時推定処理の手順を具体例として説明する。この例では、拍周期bの値が「3」、「4」及び「5」以外である状態の観測尤度が十分に小さいものとし、図17乃至図19では、拍周期bの値が「3」、「4」及び「5」以外である状態の観測尤度の図示を省略する。また、この例では、拍周期bの値が「βs」であり、且つフレーム数nの値が「ηs」である状態から、拍周期bの値が「βe」であり、且つフレーム数nの値が「ηe」である状態への対数遷移確率Tの値は、次のように設定されている。「ηe=0」、「βe=βs」、かつ「ηe=βe−1」のとき、対数遷移確率Tの値は、「−0.2」である。また、「ηs=0」、「βe=βs+1」、かつ「ηe=βe−1」のとき、対数遷移確率Tの値は、「−0.6」である。また、「ηs=0」、「βe=βs−1」、かつ「ηe=βe−1」のとき、対数遷移確率Tの値は、「−0.6」である。また、「ηs>0」、「βe=βs」、かつ「ηe=ηs−1」のとき、対数遷移確率Tの値は、「0」である。上記以外の対数遷移確率Tの値は、「−∞」である。すなわち、フレーム数nの値が「0」である状態(ηs=0)から次の状態へ遷移するとき、拍周期bの値は「1」だけ増減され得る。このとき、フレーム数nの値は、遷移後の拍周期bの値より「1」だけ小さい値に設定される。また、フレーム数nの値が「0」でない状態(ηs≠0)から次の状態へ遷移するとき、拍周期bの値は変更されず、フレーム数nの値が「1」だけ減少する。
以下、拍・テンポ同時推定処理について具体的に説明する。CPU12aは、ステップS181にて拍・テンポ同時推定処理を開始する。つぎに、ユーザは、ステップS182にて、入力操作子11を用いて、図18に示すような、各状態qb,nに対応した、尤度Cの初期条件CSb,nを入力する。なお、初期条件CSb,nがROM12bに記憶されていて、CPU12aがROM12bから初期条件CSb,nを読み込むようにしてもよい。
つぎに、CPU12aは、ステップS183にて、尤度Cb,n(ti)及び状態Ib,n(ti)を計算する。フレームt0において拍周期bの値が「βe」であって、フレーム数nの値が「ηe」である状態qb=βe,n=ηeの尤度Cb=βe,n=ηe(t0)は、初期条件CSb=βe,n=ηeと対数観測尤度Lb=βe,n=ηe(t0)とを加算することにより計算される。
また、状態qb=βs,n=ηsから状態qb=βe,n=ηeに遷移したとき、尤度Cb=βe,n=ηe(ti){i>0}は次のように計算される。状態qb=βs,n=ηsのフレーム数nが「0」でないとき(すなわち、ηs≠0)、尤度Cb=βe,n=ηe(ti)は、尤度Cb=βe,n=ηe+1(ti―1)と対数観測尤度Lb=βe,n=ηe(ti)と対数遷移確率Tを加算して計算される。ただし、本実施形態では、遷移元の状態のフレーム数nが「0」でないときの対数遷移確率Tは「0」であるので、尤度Cb=βe,n=ηe(ti)は、実質的には、尤度Cb=βe,n=ηe+1(ti―1)と対数観測尤度Lb=βe,n=ηe(ti)とを加算することにより計算される(Cb=βe,n=ηe(ti)=Cb=βe,n=ηe+1(ti―1)+Lb=βe,n=ηe(ti))。また、この場合、状態Ib=βe,n=ηe(ti)は、状態qβe,ηe+1である。例えば、尤度Cが図18に示すように計算された例では、尤度C4,1(t2)の値は「2」であり、対数観測尤度L4,0(t3)の値は「1」であるので、尤度C4,0(t3)の値は「3」である。また、図19に示すように、状態I4,0(t3)は、状態q4,1である。
また、状態qb=βs,n=ηsのフレーム数nが「0」のとき(ηs=0)の尤度Cb=βe,n=ηe(ti)は次のように計算される。この場合、状態の遷移に伴って拍周期bの値が増減され得る。そこで、まず、尤度Cβe−1,0(ti−1)、尤度Cβe,0(ti−1)、及び尤度Cβe+1,0(ti−1)に対数遷移確率Tをそれぞれ加算し、そのうちの最大値に対数観測尤度Lb=βe,n=ηe(ti)を加算した結果が尤度Cb=βe,n=ηe(ti)である。また、状態Ib=βe,n=ηe(ti)は、状態qβe−1,0、状態qβe,0、及び状態qβe+1,0のうち、それらの尤度Cβe−1,0(ti−1)、尤度Cβe,0(ti−1)、及び尤度Cβe+1,0(ti−1)に対数遷移確率Tをそれぞれ加算した値が最大となる状態qである。なお、厳密には、尤度Cb,n(ti)は正規化される必要があるが、正規化されていなくても、拍点及びテンポの推移の推定に関しては、数理上同一の結果が得られる。
例えば、尤度C4,3(t4)は、次のように計算される。遷移元の状態が状態q3,0である場合、尤度C3,0(t3)の値は「0.4」であり、対数遷移確率Tは「−0.6」であるので、尤度C3,0(t3)と対数遷移確率Tとを加算した値は、「−0.2」である。また、遷移元の状態が状態q4,0である場合、遷移元の尤度C4,0(t3)の値は「3」であり、対数遷移確率Tは「−0.2」であるので、尤度C4,0(t3)と対数遷移確率Tとを加算した値は、「2.8」である。また、遷移元の状態が状態q5,0である場合、遷移元の尤度C5,0(t3)の値は「1」であり、対数遷移確率Tは「−0.6」であるので、尤度C5,0(t3)と対数遷移確率Tとを加算した値は、「0.4」である。したがって、尤度C4,0(t3)に対数遷移確率Tを加算した値が最も大きい。また、対数観測尤度L4,3(t4)の値は、「0」である。よって、尤度C4,3(t4)の値は「2.8」(=2.8+0)であり、状態I4,3(t4)は、状態q4,0である。
上記のようにして、全てのフレームtiについて、全ての状態qb,nの尤度Cb,n(ti)及び状態Ib,n(ti)を計算し終えると、CPU12aはステップS184にて、最尤の状態の系列Q(={qmax(t0),qmax(t1)・・・,qmax(tlast)})を次のようにして決定する。まず、CPU12aは、フレームtlastにおける尤度Cb,n(tlast)が最大である状態qb,nを、状態qmax(tlast)とする。ここで、状態qmax(tlast)の拍周期bの値を「βm」と表記し、フレーム数nの値を「ηm」と表記する。このとき、状態Iβm,ηm(tlast)がフレームtlastの1つ前のフレームtlast−1の状態qmax(tlast−1)である。フレームtlast−2、フレームtlast−3、・・・の状態qmax(tlast−2)、状態qmax(tlast−3)・・・も状態qmax(tlast−1)と同様に決定される。すなわち、フレームti+1の状態qmax(ti+1)の拍周期bの値を「βm」と表記し、フレーム数nの値を「ηm」と表記したときの状態Iβm,ηm(ti+1)がフレームti+1の1つ前のフレームtiの状態qmax(ti)である。上記のようにして、CPU12aは、フレームtlast−1からフレームt0へ向かって順に状態qmaxを決定して、最尤の状態の系列Qを決定する。
例えば、図18及び図19に示す例では、フレームtlast=9においては、状態q4,2の尤度C4,2(tlast=9)が最大である。したがって、状態qmax(tlast=9)は、状態q4,2である。図19によれば、状態I4,2(t9)は状態q4,3であるから、状態qmax(t8)は状態q4,3である。また、状態I4,3(t8)は状態q4,0であるから、状態qmax(t7)は状態q4,0である。状態qmax(t6)乃至状態qmax(t0)も状態qmax(t8)及び状態qmax(t7)と同様に決定する。このようにして図18に矢印で示す最尤の状態の系列Qが決定される。この例では、拍の周期bの値はいずれのフレームtiにおいても「4」であると推定される。また、系列Qのうち、フレーム数nの値が「0」である状態qmax(t1),qmax(t5),qmax(t8)に対応するフレームt1,t5,t8に拍が存在すると推定される。
つぎに、CPU12aは、ステップS185にて、拍・テンポ同時推定処理を終了し、その処理を音響信号分析処理(メインルーチン)のステップS19に進める。
CPU12aは、ステップS19にて、フレームtiごとに「BPMらしさ」、「観測に基づく確率」、「拍らしさ」、「拍が存在する確率」及び「拍が存在しない確率」を計算(図20に示す演算式を参照)する。「BPMらしさ」は、フレームtiにおけるテンポの値が拍周期bに対応した値である確率を意味し、尤度Cb,n(ti)を正規化するとともにフレーム数nについて周辺化することにより計算される。具体的には、拍周期bの値が「β」である場合の「BPMらしさ」は、フレームtiにおける全ての状態の尤度Cの合計に対する、拍周期bの値が「β」である状態の尤度Cの合計の割合である。また、「観測に基づく確率」は、観測値(すなわちオンセット特徴量XO)に基づいて計算された拍がフレームtiに存在する確率を意味する。具体的には、所定の基準値XObaseに対するオンセット特徴量XO(ti)の割合である。また、「拍らしさ」は、すべてのフレーム数nの値についてのオンセット特徴量XO(ti)の尤度P(XO(ti)|Zb,n(ti))を合算した値に対する尤度P(XO(ti)|Zb,0(ti))の割合である。また、「拍が存在する確率」及び「拍が存在しない確率」は、いずれも尤度Cb,n(ti)を拍周期bについて周辺化することにより計算される。具体的には、「拍が存在する確率」は、フレームtiにおける全ての状態の尤度Cの合計に対する、フレーム数nの値が「0」である状態の尤度Cの合計の割合である。また、「拍が存在しない確率」は、フレームtiにおける全ての状態の尤度Cの合計に対する、フレーム数nの値が「0」でない状態の尤度Cの合計の割合である。
CPU12aは、「BPMらしさ」、「観測に基づく確率」、「拍らしさ」、「拍が存在する確率」及び「拍が存在しない確率」を用いて、図20に示す拍・テンポ情報リストを表示器13に表示する。同リスト中の「テンポの推定値(BPM)」の欄には、前記計算した「BPMらしさ」のうち最も確率の高い拍周期bに対応するテンポの値(BPM)が表示される。また、前記決定した状態qmax(ti)のうちフレーム数nの値が「0」であるフレームの「拍の存在」の欄には「○」が表示され、その他のフレームの「拍の存在」の欄には「×」が表示される。また、CPU12aは、テンポの推定値(BPM)を用いて、図21に示すようなテンポの推移を表わすグラフを表示器13に表示する。図21の例では、テンポの推移を棒グラフで表わしている。図18及び図19を用いて説明した例では、テンポの値が一定であるので図21に示すような各フレームのテンポを表わすバーの高さが一定であるが、テンポが頻繁に変化する楽曲では、図22に示すように、テンポの値に応じてバーの高さが異なる。これにより、ユーザは、テンポの推移を視覚的に認識することができる。また、CPU12aは、前記計算した「拍が存在する確率」を用いて、図23に示すような拍点を表わすグラフを表示器13に表示する。
また、音響信号分析処理のステップS13にて既存データを検索した結果、既存データが存在する場合には、CPU12aは、ステップS15にてRAM12cに読み込んだ前回の分析結果に関する各種データを用いて、拍・テンポ情報リスト、テンポの推移を表わすグラフ、及び拍点を表わすグラフを表示器13に表示する。
次に、CPU12aは、ステップS20にて、音響信号分析処理を終了するか否かを表すメッセージを表示器13に表示して、ユーザからの指示を待つ。ユーザは入力操作子11を用いて音響信号分析処理を終了するか、後述の拍・テンポ情報修正処理を実行するかのいずれかを指示する。例えば、マウスを用いて図示しないアイコンをクリックする。ユーザから音響信号分析処理を終了するよう指示された場合には、CPU12aは「Yes」と判定してステップS21にて尤度C、状態I、拍・テンポ情報リストなどの分析結果に関する各種データを楽曲のタイトルと関連付けて記憶装置14に記憶して、ステップS22にて音響信号分析処理を終了する。
一方、ステップS20にて、音響信号分析処理を継続するように指示された場合には、CPU12aは「No」と判定して、ステップS23にて、テンポ情報修正処理を実行する。まず、CPU12aは、ユーザが修正情報の入力を終了するまで待機する。ユーザは、入力操作子11を用いて「BPMらしさ」、「拍が存在する確率」などの修正値を入力する。例えば、マウスを用いて修正するフレームを選択し、テンキーを用いて修正値を入力する。修正された項目の右側に配置された「F」の表示形態(例えば色)が変更され、その値が修正されたことが明示される。ユーザは、複数の項目について修正値を入力可能である。ユーザは修正値の入力を完了すると、入力操作子11を用いて修正情報の入力を完了したことを指示する。例えば、マウスを用いて図示しない修正完了を表わすアイコンをクリックする。CPU12aは、前記入力された修正値に応じて尤度P(XO(ti)|Zb,n(ti))及び尤度P(XB(ti)|Zb,n(ti))のうちのいずれか一方又は両方を更新する。例えば、フレームtiにおける「拍が存在する確率」が高くなるように修正された場合であって、修正された値に関するフレーム数nの値が「ηe」であるときには、尤度P(XB(ti)|Zb,n≠ηe(ti))を十分に小さい値に設定する。これにより、フレームtiでは、フレーム数nの値が「ηe」である確率が相対的に最も高くなる。また、例えば、フレームtiにおける「BPMらしさ」のうち、拍周期bの値が「βe」である確率が高くなるように修正された場合には、拍周期bの値が「βe」でない状態の尤度P(XB(ti)|Zb≠βe,n(ti))を十分に小さい値に設定する。これにより、フレームtiでは、拍周期bの値が「βe」である確率が相対的に最も高くなる。そして、CPU12aは、拍・テンポ情報修正処理を終了して、その処理をステップS18に進め、修正された対数観測尤度Lを用いて、拍・テンポ同時推定処理を再度実行する。
上記のように構成した音響信号分析装置10によれば、拍点に関するオンセット特徴量XO及びテンポに関するBPM特徴量XBを用いて計算された対数観測尤度Lの系列が最も尤もらしい確率モデルが選択され、楽曲における拍点及びテンポの推移が同時に推定される。したがって、上記従来技術とは異なり、拍点及びテンポのうちの一方の推定精度が低いために他方の推定精度も低くなるという事態が生じない。よって、従来技術に比べて楽曲における拍点及びテンポの推移の推定精度を向上させることができる。
また、本実施形態においては、フレーム数nの値が「0」である状態から、拍周期bの値が同じ状態又は拍周期bの値が「1」だけ異なる状態へのみ遷移可能に各状態間の遷移確率(対数遷移確率)が設定されている。これにより、テンポがフレーム間で急激に変化するような誤推定が防止される。したがって、楽曲として自然な拍点及びテンポの推移の推定結果を得ることができる。なお、テンポが急激に変化する楽曲に対しては、次の拍までのフレーム数nの値が「0」である状態から次の状態に遷移するとき、拍周期bの値が大きく異なる状態への遷移も可能なように各状態間の遷移確率(対数遷移確率)を設定すればよい。
また、拍・テンポ同時推定処理では、ビタビアルゴリズムを用いたので、他のアルゴリズム(例えば、「サンプリング法」、「前向き後向きアルゴリズム」など)を用いる場合に比べて計算量を削減できる。
また、ユーザにより入力された修正情報に基づいて対数観測尤度Lが修正され、修正された対数観測尤度Lに基づいて楽曲における拍点及びテンポの推移が再推定される。これにより、修正されたフレームの前後にそれぞれ位置する1つ又は複数のフレームの最尤の状態qmaxが再計算(再選択)される。したがって、修正されたフレーム及びその前後に位置する1つ又は複数のフレームに亘り、拍の間隔及びテンポが滑らかに変化するような推定結果が得られる。
上記のようにして推定された楽曲における拍点及びテンポの推移に関する情報は、例えば楽曲データの検索、伴奏を表わす伴奏データの検索などに利用される。また、分析対象とした楽曲に対する伴奏パートの自動生成、ハーモニーの自動付加などにも利用される。
さらに、本発明の実施にあたっては、上記実施形態に限定されるものではなく、本発明の目的を逸脱しない限りにおいて種々の変更が可能である。
例えば、上記実施形態では、観測値としてのオンセット特徴量XO及びBPM特徴量XBが同時に観測される確率を表わす観測尤度の系列が最も尤もらしい確率モデルが選択される。しかし、確率モデルの選択基準は、上記実施形態に限られない。例えば、事後分布が最大となるような確率モデルを選択してもよい。
また、例えば、上記実施形態では、説明を簡単にするために、各フレームの長さを125msとしたが、より短く(例えば、5ms)してもよい。これによれば、拍点及びテンポの推定に関する分解能を向上させることができる。例えば、テンポを1BPM刻みで推定できる。また、上記実施形態では、各フレームの長さを共通にしているが、各フレームの長さが異なっていてもよい。この場合であっても、オンセット特徴量XOは、上記実施形態と同様にして計算できる。また、この場合、BPM特徴量XBの計算においては、コムフィルタの遅延量をフレームの長さに応じて変更すればよい。また、尤度Cの計算においては、各フレームの長さの最大公約数F(つまり、各フレームを構成するサンプル数の
最大公約数)を計算する。そして、フレームti(=τ)の長さがL(τ)×Fと表わされたとき、状態qb,n(n≠0)から状態qb,n−L(τ)へ遷移する確率を100%とすればよい。
また、上記実施形態では、楽曲全体を分析対象としているが、楽曲の一部(例えば数小節)のみを分析対象としてもよい。この場合、入力した楽曲データのうち、分析対象とする部分を選択可能に構成するとよい。また、楽曲のうちの単一のパート(例えばリズムセクション)のみを分析対象としてもよい。
また、例えば、テンポの推定において、優先的に推定するテンポの範囲を指定可能に構成してもよい。具体的には、音響信号分析処理のステップS12において、「Presto」、「Moderato」などのテンポを表わす用語を表示して、優先的に推定するテンポの範囲を選択可能に構成してもよい。例えば、「Presto」が選択された場合、BPM=160〜190の範囲以外の対数観測尤度Lを十分に小さく設定する。これにより、BPM=160〜190の範囲のテンポが優先的に推定される。これによれば、楽曲の大凡のテンポが既知である場合、テンポの推定精度を向上させることができる。
また、拍・テンポ情報修正処理(ステップS23)では、ユーザは入力操作子11を用いて修正内容を入力するように構成されている。これに代えて、又は加えて、外部インターフェース回路15を介して接続された電子鍵盤楽器、電子打楽器などの操作子を用いて修正内容を入力可能に構成してもよい。例えば、ユーザが電子鍵盤楽器の鍵盤を数回打鍵すると、CPU12aがその打鍵のタイミングからテンポを計算して、前記「BPMらしさ」の修正値として用いるように構成してもよい。
また、上記実施形態では、拍点及びテンポに関する修正値を何度でも入力可能に構成されている。しかし、例えば、「拍が存在する確率」の平均値が基準値(例えば80%)に達した時点以降においては、拍点及びテンポに関する修正値を入力不可能としてもよい。
また、例えば、拍・テンポ情報修正処理(ステップS23)において、ユーザによって指定されたフレームの拍・テンポ情報を入力された値に修正するとともに、そのフレームに近接するフレームの拍・テンポ情報を前記入力された値に応じて自動的に修正してもよい。例えば、連続する複数のフレームのテンポの推定値が同じ値であって、そのうちの1つのフレームのテンポの値が修正されたとき、前記複数のフレームのテンポの値を前記1つのフレームの修正値と同じ値に自動的に修正してもよい。
また、上記実施形態では、ステップS23にて、ユーザが入力操作子11を用いて修正値の入力を完了したことを指示すると、拍点及びテンポの同時推定が再び実行される。しかし、これに代えて、ユーザが少なくとも1つの修正値を入力した後、他の修正値が入力されないまま所定の時間(例えば10秒)が経過したとき、自動的に拍点及びテンポの同時推定が再び実行されてもよい。
また、拍・テンポ情報リスト(図20)の表示形態は、上記実施形態に限られない。例えば、上記実施形態では、「BPMらしさ」、「拍点らしさ」などを、確率(%)で表示しているが、これらを記号、文字列などを用いて表現してもよい。また、上記実施形態では、前記決定した状態qmax(ti)のうちフレーム数nの値が「0」であるフレームtiの「拍の存在」の欄には「○」が表示され、その他のフレームの「拍の存在」の欄には「×」が表示されるが、これに代えて、例えば、基準値(例えば80%)以上であるとき「拍の存在」の欄に「○」が表示され、「拍点が存在する確率」が基準値未満であるとき、「拍の存在」の欄に「×」が表示されてもよい。また、この場合、複数の基準値を設けてもよい。例えば、第1基準値(=80%)と第2基準値(=60%)を設け、「拍点が存在する確率」が第1基準値以上であるとき、「拍の存在」の欄に「○」が表示され、第2基準値以上かつ第1基準値未満であるとき、「拍の存在」の欄に「△」が表示され、「拍点が存在する確率」が第2基準値未満であるとき、「拍の存在」の欄に「×」が表示されてもよい。また、テンポの推定値の欄には、「Presto」、「Moderato」などのテンポを表わす用語が表示されてもよい。