(57)【要約】【課題】楽曲における拍点及びテンポの推移を高精度に検出できる音響信号分析装置を提供する。【解決手段】音響信号分析装置１０は、楽曲を表わす音響信号を入力する音響信号入力手段と、前記楽曲の各区間における前記拍の存在に関する特徴を表わす第１特徴量及びテンポに関する特徴を表わす第２特徴量を計算する特徴量計算手段と、前記楽曲の各区間における拍の存在に関する物理量及びテンポに関する物理量の組み合わせにより分類された状態ｑの系列として記述された複数の確率モデルのうち、前記第１特徴量及び前記第２特徴量が前記楽曲の各区間において同時に観測される確率を表わす観測尤度の系列が所定の基準を満たす確率モデルを選択することにより、前記楽曲における拍点及びテンポの推移を同時に推定する推定手段と、を備える。【選択図】図２
2014178394.tif 000002
【特許請求の範囲】
【請求項１】
楽曲を表わす音響信号を入力する音響信号入力手段と、
前記楽曲の各区間における拍の存在に関する特徴を表わす第１特徴量及びテンポに関する特徴を表わす第２特徴量を計算する特徴量計算手段と、
前記楽曲の各区間における拍の存在に関する物理量及びテンポに関する物理量の組み合わせにより分類された状態の系列として記述された複数の確率モデルのうち、前記第１特徴量及び前記第２特徴量が前記楽曲の各区間において同時に観測される確率を表わす観測尤度の系列が所定の基準を満たす確率モデルを選択することにより、前記楽曲における拍点及びテンポの推移を同時に推定する推定手段と、
を備えたことを特徴とする音響信号分析装置。
【請求項２】
請求項１に記載の音響信号分析装置において、
前記推定手段は、
前記拍の存在に関する物理量に応じて規定された確率分布関数の確率変数として前記第１特徴量を代入して計算された確率を前記第１特徴量が観測される確率として出力する第１確率出力手段と、
前記テンポに関する物理量に応じて規定された確率分布関数の確率変数として前記第２特徴量を代入して計算された確率を前記第２特徴量が観測される確率として出力する第２確率出力手段と、
を備えたことを特徴とする音響信号分析装置。
【請求項３】
請求項２に記載の音響信号分析装置において、
前記第１確率出力手段は、前記拍の存在に関する物理量に応じて規定された正規分布の確率変数として前記第１特徴量を代入して計算された確率を前記第１特徴量が観測される確率として出力することを特徴とする音響信号分析装置。
【請求項４】
請求項２に記載の音響信号分析装置において、
前記第２確率出力手段は、前記テンポに関する物理量に応じてそれぞれ設定された複数のテンプレートに対する前記第２特徴量の適合度を前記第２特徴量が観測される確率として出力することを特徴とする音響信号分析装置。
【請求項５】
請求項１乃至４のうちのいずれか１つに記載の音響信号分析装置において、
前記楽曲の各区間は、前記入力した音響信号を所定の時間間隔をおいて区切ることにより形成された各フレームにそれぞれ相当し、
前記特徴量計算手段は、
前記フレームごとに振幅スペクトルを計算し、前記振幅スペクトルに周波数帯域が互いに異なる複数の窓関数を掛けて前記周波数帯域ごとの振幅スペクトルを生成し、前記フレーム間における前記周波数帯域ごとの振幅スペクトルの変化に基づいて計算した値を前記第１特徴量として出力する第１特徴量計算手段と、
各フレームに対応する値を入力するごとに１つの値を出力するフィルタであって、出力した値を所定の時間だけ保持する保持手段を備え、前記入力した値と前記所定の時間だけ保持された値とを所定の比率で加算して出力するフィルタを備え、
前記第１特徴量の系列を前記フィルタに入力して得られたデータ列の時系列を逆にしたデータ列を前記フィルタに再び入力して得られたデータ列を前記第２特徴量の系列として出力する第２特徴量計算手段と、を備えたことを特徴とする音響信号分析装置。
【請求項６】
請求項１乃至５のうちのいずれか１つに記載の音響信号分析装置において、
前記楽曲における拍点及びテンポの推移のいずれか一方又は両方の修正内容を表わす修正情報を入力する修正情報入力手段と、
前記入力した修正情報に応じて前記観測尤度を修正する観測尤度修正手段と、
前記複数の確率モデルのうち、前記修正された観測尤度の系列が前記所定の基準を満たす確率モデルを、前記推定手段を用いて選択することにより、前記楽曲における拍点及びテンポの推移を同時に再推定する再推定手段と、を備えたことを特徴とする音響信号分析装置。
【請求項７】
コンピュータに、
楽曲を表わす音響信号を入力する音響信号入力ステップと、
前記楽曲の各区間における拍の存在に関する特徴を表わす第１特徴量及びテンポに関する特徴を表わす第２特徴量を計算する特徴量計算ステップと、
前記楽曲の各区間における拍の存在に関する物理量及びテンポに関する物理量の組み合わせにより分類された状態の系列として記述された複数の確率モデルのうち、前記第１特徴量及び前記第２特徴量が前記楽曲の各区間において同時に観測される確率を表わす観測尤度の系列が所定の基準を満たす確率モデルを選択することにより、前記楽曲における拍点及びテンポの推移を同時に推定する推定ステップと、
を実行させることを特徴とする音響信号分析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、楽曲を表わす音響信号を入力して、前記入力した楽曲における拍点（拍のタイミング）及びテンポを検出する音響信号分析装置に関する。
【背景技術】
【０００２】
従来から、例えば、下記特許文献１に示されているように、楽曲を表わす音響信号を入力して楽曲における拍点及びテンポを検出する音響信号分析装置は知られている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００９−２６５４９３号公報
【発明の概要】
【０００４】
上記特許文献１の音響信号分析装置は、まず、音響信号の強度（振幅）の変化に基づいて、拍点の候補として拍指標系列を計算する。次に、拍指標系列の計算結果に基づいて、楽曲のテンポを検出する。したがって、拍指標系列の検出精度が低い場合、テンポの検出精度も低くなる。
【０００５】
本発明は上記問題に対処するためになされたもので、その目的は、楽曲における拍点及びテンポの推移を高精度に検出できる音響信号分析装置を提供することにある。なお、下記本発明の各構成要件の記載においては、本発明の理解を容易にするために、実施形態の対応箇所の符号を括弧内に記載しているが、本発明の各構成要件は、実施形態の符号によって示された対応箇所の構成に限定解釈されるべきものではない。
【０００６】
上記目的を達成するために、本発明の特徴は、楽曲を表わす音響信号を入力する音響信号入力手段（Ｓ１２）と、前記楽曲の各区間における拍の存在に関する特徴を表わす第１特徴量（ＸＯ）及びテンポに関する特徴を表わす第２特徴量（ＸＢ）を計算する特徴量計算手段（Ｓ１６５，Ｓ１６７）と、前記楽曲の各区間における拍の存在に関する物理量（ｎ）及びテンポに関する物理量（ｂ）の組み合わせにより分類された状態（ｑｂ，ｎ）の系列として記述された複数の確率モデルのうち、前記第１特徴量及び前記第２特徴量が前記楽曲の各区間において同時に観測される確率を表わす観測尤度（Ｌ）の系列が所定の基準を満たす確率モデルを選択することにより、前記楽曲における拍点及びテンポの推移を同時に推定する推定手段（Ｓ１７，Ｓ１８）と、を備えたことにある。
【０００７】
この場合、推定手段は、拍の存在に関する物理量に応じて規定された確率分布関数の確率変数として第１特徴量を代入して計算された確率を第１特徴量が観測される確率として出力する第１確率出力手段（Ｓ１７２）と、テンポに関する物理量の大きさに応じて規定された確率分布関数の確率変数として第２特徴量を代入して計算された確率を第２特徴量が観測される確率として出力する第２確率出力手段（Ｓ１７３）と、を備えるとよい。
【０００８】
この場合、第１確率出力手段は、前記拍の存在に関する物理量に応じて規定された正規分布の確率変数として第１特徴量を代入して計算された確率を前記第１特徴量が観測される確率として出力するとよい。
【０００９】
また、この場合、第２確率出力手段は、前記テンポに関する物理量に応じてそれぞれ設定された複数のテンプレートに対する第２特徴量の適合度を前記第２特徴量が観測される確率として出力するとよい。
【００１０】
また、この場合、前記楽曲の各区間は、前記入力した音響信号を所定の時間間隔をおいて区切ることにより形成された各フレームにそれぞれ相当し、特徴量計算手段は、フレームごとに振幅スペクトル（Ａ）を計算し、前記振幅スペクトルに周波数帯域（ｗｋ）が互いに異なる複数の窓関数（ＢＰＦ）を掛けて前記周波数帯域ごとの振幅スペクトル（Ｍ）を生成し、フレーム間における前記周波数帯域ごとの振幅スペクトルの変化に基づいて計算した値を第１特徴量として出力する第１特徴量計算手段（Ｓ１６５）と、各フレームに対応する値を入力するごとに１つの値を出力するフィルタであって、出力した値を所定の時間だけ保持する保持手段（ｄｂ）を備え、前記入力した値と前記所定の時間だけ保持された値とを所定の比率で加算して出力するフィルタ（ＦＢＢ）を備え、第１特徴量の系列をフィルタに入力して得られたデータ列の時系列を逆にしたデータ列を前記フィルタに再び入力して得られたデータ列を第２特徴量の系列として出力する第２特徴量計算手段（Ｓ１６７）と、を備えるとよい。
【００１１】
上記のように構成した音響信号分析装置によれば、拍の存在に関する特徴を表わす第１特徴量及びテンポに関する特徴を表わす第２特徴量を用いて計算された観測尤度の系列が所定の基準を満たす確率モデル（例えば、最も尤もらしい確率モデル、事後分布が最大となる確率モデルなど）が選択され、楽曲における拍点及びテンポの推移が同時に推定される。したがって、上記従来技術とは異なり、拍点及びテンポのうちの一方の推定精度が低いために他方の推定精度も低くなるという事態が生じない。よって、従来技術に比べて拍点及びテンポの推移の推定精度を向上させることができる。
【００１２】
また、本発明の他の特徴は、前記楽曲における拍点及びテンポの推移のいずれか一方又は両方の修正内容を表わす修正情報を入力する入力手段（１１，Ｓ２３）と、前記入力した修正情報に応じて前記観測尤度を修正する観測尤度修正手段（Ｓ２３）と、前記複数の確率モデルのうち、前記修正された観測尤度の系列が前記所定の基準を満たす確率モデルを、推定手段を用いて選択することにより、前記楽曲における拍点及びテンポの推移を同時に再推定する再推定手段（Ｓ２３，Ｓ１８）と、を備えたことにある。
【００１３】
これによれば、ユーザにより入力された修正情報に基づいて観測尤度が修正され、修正された観測尤度に基づいて楽曲における拍点及びテンポの推移が再推定される。つまり、修正されたフレームの前後にそれぞれ位置する１つ又は複数のフレームの状態が再計算（再選択）される。したがって、修正されたフレーム及びその前後に位置する１つ又は複数のフレームに亘り、拍の間隔（すなわち、テンポ）が滑らかに変化するような推定結果が得られる。
【００１４】
さらに、本発明の実施にあたっては、音響信号分析装置の発明に限定されることなく、同装置に適用されるコンピュータプログラムの発明としても実施し得るものである。
【図面の簡単な説明】
【００１５】
【図１】本発明の一実施形態に係る音響信号分析装置の全体構成を示すブロック図である。
【図２】確率モデルの概念図である。
【図３】音響信号分析プログラムを表わすフローチャートである。
【図４】特徴量計算プログラムを表わすフローチャートである。
【図５】分析対象の音響信号の波形を表わすグラフである。
【図６】１つのフレームを短時間フーリエ変換した音響スペクトル図である。
【図７】バンドパスフィルタの特性図である。
【図８】各周波数帯域の振幅の時間変化を示すグラフである。
【図９】オンセット特徴量の時間変化を示すグラフである。
【図１０】コムフィルタのブロック図である。
【図１１】ＢＰＭ特徴量の計算結果を示すグラフである。
【図１２】対数観測尤度計算プログラムを表わすフローチャートである。
【図１３】オンセット特徴量の観測尤度の計算結果を示す表である。
【図１４】テンプレートの構成を示す表である。
【図１５】ＢＰＭ特徴量の観測尤度の計算結果を示す表である。
【図１６】拍・テンポ同時推定プログラムを表わすフローチャート
【図１７】対数観測尤度の計算結果を示す表である。
【図１８】先頭のフレームから各フレームまでオンセット特徴量及びＢＰＭ特徴量を観測したときに前記各フレームの各状態の尤度が最大となるような状態の系列を選択した場合の前記各状態の尤度の計算結果を示す表である。
【図１９】遷移元の状態の計算結果を示す表である。
【図２０】拍・テンポ情報リストの概略を示す概略図である。
【図２１】テンポの推移の一例を示すグラフである。
【図２２】テンポの推移の他の例を示すグラフである。
【図２３】拍点を示すグラフである。
【発明を実施するための形態】
【００１６】
本発明の一実施形態に係る音響信号分析装置１０について説明する。音響信号分析装置１０は、以下説明するように、楽曲を表わす音響信号を入力して、その楽曲における拍点及びテンポの推移を検出する。音響信号分析装置１０は、図１に示すように、入力操作子１１、コンピュータ部１２、表示器１３、記憶装置１４、外部インターフェース回路１５及びサウンドシステム１６を備えており、これらがバスＢＳを介して接続されている。
【００１７】
入力操作子１１は、オン・オフ操作に対応したスイッチ（例えば数値を入力するためのテンキー）、回転操作に対応したボリューム又はロータリーエンコーダ、スライド操作に対応したボリューム又はリニアエンコーダ、マウス、タッチパネルなどから構成される。これらの操作子は、演奏者の手によって操作されて、分析対象の楽曲の選択、音響信号の分析開始又は停止、楽曲の再生又は停止（後述するサウンドシステム１６からの出力又は停止）、音響信号の分析に関する各種パラメータの設定などに用いられる。入力操作子１１を操作すると、その操作内容を表す操作情報が、バスＢＳを介して、後述するコンピュータ部１２に供給される。
【００１８】
コンピュータ部１２は、バスＢＳにそれぞれ接続されたＣＰＵ１２ａ、ＲＯＭ１２ｂ及びＲＡＭ１２ｃからなる。ＣＰＵ１２ａは、詳しくは後述する音響信号分析プログラム及びそのサブルーチンをＲＯＭ１２ｂから読み出して実行する。ＲＯＭ１２ｂには、音響信号分析プログラム及びそのサブルーチンに加えて、初期設定パラメータ、表示器１３に表示される画像を表わす表示データを生成するための図形データ及び文字データなどの各種データが記憶されている。ＲＡＭ１２ｃには、音響信号分析プログラムの実行時に必要なデータが一時的に記憶される。
【００１９】
表示器１３は、液晶ディスプレイ（ＬＣＤ）によって構成される。コンピュータ部１２は、図形データ、文字データなどを用いて表示すべき内容を表わす表示データを生成して表示器１３に供給する。表示器１３は、コンピュータ部１２から供給された表示データに基づいて画像を表示する。例えば分析対象の楽曲の選択時には、楽曲のタイトルリストが表示される。また、例えば分析終了時には、拍点及びテンポの推移を表わす拍・テンポ情報リスト及びそのグラフ（図２０乃至図２３参照）が表示される。
【００２０】
また、記憶装置１４は、ＨＤＤ、ＦＤＤ、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤなどの大容量の不揮発性記録媒体と、同各記録媒体に対応するドライブユニットから構成されている。記憶装置１４には、複数の楽曲をそれぞれ表わす複数の楽曲データが記憶されている。楽曲データは、楽曲を所定のサンプリング周期（例えば４４．１ｋＨｚ）でサンプリングして得られた複数のサンプル値からなり、各サンプル値が記憶装置１４における連続するアドレスに順に記録されている。楽曲のタイトルを表わすタイトル情報、楽曲データの容量を表わすデータサイズ情報なども楽曲データに含まれている。楽曲データは予め記憶装置１４に記憶されていてもよいし、後述する外部インターフェース回路１５を介して外部から取り込んでもよい。記憶装置１４に記憶されている楽曲データは、ＣＰＵ１２ａによって読み込まれ、楽曲における拍点及びテンポの推移が分析される。
【００２１】
外部インターフェース回路１５は、音響信号分析装置１０を電子音楽装置、パーソナルコンピュータなどの外部機器に接続可能とする接続端子を備えている。音響信号分析装置１０は、外部インターフェース回路１５を介して、ＬＡＮ（Ｌｏｃａｌ  Ａｒｅａ  Ｎｅｔｗｏｒｋ）、インターネットなどの通信ネットワークにも接続可能である。
【００２２】
サウンドシステム１６は、楽曲データをアナログ音信号に変換するＤ／Ａ変換器、変換したアナログ音信号を増幅するアンプ、及び増幅されたアナログ音信号を音響信号に変換して出力する左右一対のスピーカを備えている。ユーザが入力操作子１１を用いて分析対象の楽曲の再生を指示すると、ＣＰＵ１２ａは、分析対象の楽曲データをサウンドシステム１６に供給する。これにより、ユーザは分析対象の楽曲を試聴できる。
【００２３】
つぎに、上記のように構成した音響信号分析装置１０の動作について説明する。まず、その概略について説明する。分析対象の楽曲は複数のフレームｔｉ｛ｉ＝０，１，・・・，ｌａｓｔ｝に分割される。そして、拍の存在に関する特徴を表すオンセット特徴量ＸＯ及びテンポに関する特徴を表すＢＰＭ特徴量ＸＢがフレームｔｉごとに計算される。各フレームｔｉにおける拍周期ｂの値（テンポの逆数に比例する値）及び次の拍までのフレーム数ｎの値の組み合わせに応じて分類された状態ｑｂ，ｎの系列として記述された確率モデル（隠れマルコフモデル）のうち、観測値としてのオンセット特徴量ＸＯ及びＢＰＭ特徴量ＸＢが同時に観測される確率を表わす観測尤度の系列が最も尤もらしい確率モデルが選択される（図２参照）。これにより、分析対象の楽曲における拍点及びテンポの推移が検出される。なお、拍周期ｂは、フレームの数によって表わされる。したがって、拍周期ｂの値は「１≦ｂ≦ｂｍａｘ」を満たす整数であり、拍周期ｂの値が「β」である状態では、フレーム数ｎの値は「０≦ｎ＜β」を満たす整数である。
【００２４】
つぎに、音響信号分析装置１０の動作について具体的に説明する。ユーザが音響信号分析装置１０の図示しない電源スイッチをオンにすると、ＣＰＵ１２ａは、図３に示す音響信号分析プログラムをＲＯＭ１２ｂから読み出して実行する。
【００２５】
ＣＰＵ１２ａは、ステップＳ１０にて音響信号分析処理を開始し、ステップＳ１１にて、記憶装置１４に記憶されている複数の楽曲データにそれぞれ含まれるタイトル情報を読み込んで、楽曲のタイトルをリスト形式で表示器１３に表示する。ユーザは、入力操作子１１を用いて、表示器１３に表示された楽曲の中から分析対象の楽曲データを選択する。なお、ステップＳ１１にて分析対象の楽曲データを選択する際、選択しようとする楽曲データが表す楽曲の一部又は全部を再生して楽曲データの内容を確認できるように構成してもよい。
【００２６】
つぎに、ＣＰＵ１２ａは、ステップＳ１２にて、音響信号分析のための初期設定を実行する。具体的には、前記選択された楽曲データのデータサイズ情報に応じた記憶領域をＲＡＭ１２ｃ内に確保し、前記確保した記憶領域に前記選択された楽曲データを読み込む。また、分析結果を表す拍・テンポ情報リスト、オンセット特徴量ＸＯ、ＢＰＭ特徴量ＸＢなどを一時的に記憶する領域をＲＡＭ１２ｃ内に確保する。
【００２７】
詳しくは後述するが、本プログラムによる分析結果は、記憶装置１４に保存される（ステップＳ２１）。前記選択された楽曲が本プログラムによって過去に分析されたことがあれば、記憶装置１４にその分析結果が保存されている。そこで、ＣＰＵ１２ａは、ステップＳ１３にて、前記選択された楽曲の分析に関する既存のデータ（以下、単に既存データと呼ぶ）を検索する。既存データが有れば、ＣＰＵ１２ａは、ステップＳ１４にて「Ｙｅｓ」と判定して、ステップＳ１５にて既存データをＲＡＭ１２ｃに読み込み、後述するステップＳ１９に処理を進める。一方、既存データが無ければ、ＣＰＵ１２ａは、ステップＳ１４にて「Ｎｏ」と判定して、その処理をステップＳ１６に進める。
【００２８】
ＣＰＵ１２ａは、ステップＳ１６にて、図４に示す特徴量計算プログラムをＲＯＭ１２ｂから読み出して実行する。特徴量計算プログラムは、音響信号分析プログラムのサブルーチンである。
【００２９】
ＣＰＵ１２ａは、ステップＳ１６１にて特徴量計算処理を開始する。つぎに、ＣＰＵ１２ａは、ステップＳ１６２にて、図５に示すように、前記選択された楽曲を所定の時間間隔をおいて区切り、複数のフレームｔｉ｛ｉ＝０，１，・・・，ｌａｓｔ｝に分割する。各フレームの長さは共通である。説明を簡単にするために、本実施形態では各フレームの長さを１２５ｍｓとする。上記のように、各楽曲のサンプリング周波数は４４．１ｋＨｚであるので、各フレームは、約５０００個のサンプル値から構成されている。そして、以下説明するように、フレームごとに、オンセット特徴量ＸＯ及びＢＰＭ（ｂｅａｔｓ  ｐｅｒ  ｍｉｎｕｔｅ（１分間あたりの拍数））特徴量ＸＢを計算する。
【００３０】
つぎに、ＣＰＵ１２ａは、ステップＳ１６３にて、フレームごとに短時間フーリエ変換を実行して、図６に示すように、各周波数ビンｆｊ｛ｊ＝１，２・・・｝の振幅Ａ（ｆｊ，ｔｉ）を計算する。そして、ＣＰＵ１２ａは、ステップＳ１６４にて、振幅Ａ（ｆ１，ｔｉ），Ａ（ｆ２，ｔｉ）・・・を周波数ビンｆｊごとに設けられたフィルタバンクＦＢＯｊによってフィルタ処理することにより、所定の周波数帯域ｗｋ｛ｋ＝１，２，・・・｝の振幅Ｍ（ｗｋ，ｔｉ）を計算する。周波数ビンｆｊ用のフィルタバンクＦＢＯｊは、図７に示すように、通過帯域の中心周波数が互いに異なる複数のバンドパスフィルタＢＰＦ（ｗｋ，ｆｊ）からなる。フィルタバンクＦＢＯｊを構成する各バンドパスフィルタＢＰＦ（ｗｋ，ｆｊ）の中心周波数は、対数周波数軸上で等間隔であり、かつ各バンドパスフィルタＢＰＦ（ｗｋ，ｆｊ）の通過帯域幅は、対数周波数軸上で共通である。各バンドパスフィルタＢＰＦ（ｗｋ，ｆｊ）は、通過帯域の中心周波数から通過帯域の下限周波数側及び上限周波数側へ向かうに従って徐々にゲインがそれぞれ小さくなるように構成されている。ＣＰＵ１２ａは、図４のステップＳ１６４に示すように、周波数ビンｆｊごとに振幅Ａ（ｆｊ，ｔｉ）とバンドパスフィルタＢＰＦ（ｗｋ，ｆｊ）のゲインとを積算する。そして、前記周波数ビンｆｊごとに計算した積算結果を全ての周波数ビンｆｊについて合算して振幅Ｍ（ｗｋ，ｔｉ）とする。上記のようにして計算された振幅Ｍの系列を図８に例示する。
【００３１】
つぎに、ＣＰＵ１２ａは、ステップＳ１６５にて、振幅Ｍの時間変化に基づいてフレームｔｉのオンセット特徴量ＸＯ（ｔｉ）を計算する。具体的には、図４のステップＳ１６５に示すように、周波数帯域ｗｋごとに、フレームｔｉ−１からフレームｔｉへの振幅Ｍの増加量Ｒ（ｗｋ，ｔｉ）を計算する。ただし、フレームｔｉ―１の振幅Ｍ（ｗｋ，ｔｉ−１）とフレームｔｉの振幅Ｍ（ｗｋ，ｔｉ）とが同じである場合、又はフレームｔｉの振幅Ｍ（ｗｋ，ｔｉ）がフレームｔｉ―１の振幅Ｍ（ｗｋ，ｔｉ−１）よりも小さい場合は、増加量Ｒ（ｗｋ，ｔｉ）は「０」とする。そして、周波数帯域ｗｋごとに計算した増加量Ｒ（ｗｋ，ｔｉ）を全ての周波数帯域ｗ１，ｗ２，・・・について合算してオンセット特徴量ＸＯ（ｔｉ）とする。上記のようにして計算されたオンセット特徴量ＸＯの系列を図９に例示する。一般に、楽曲においては、拍が存在する部分の音量が大きい。したがって、オンセット特徴量ＸＯ（ｔｉ）が大きいほど、フレームｔｉに拍が存在する可能性が高い。
【００３２】
つぎに、ＣＰＵ１２ａは、オンセット特徴量ＸＯ（ｔ０），ＸＯ（ｔ１）・・・を用いて、ＢＰＭ特徴量ＸＢをフレームｔｉごとに計算する。なお、フレームｔｉのＢＰＭ特徴量ＸＢ（ｔｉ）は、拍周期ｂごとに計算されたＢＰＭ特徴量ＸＢｂ＝１，２・・・（ｔｉ）の集合として表わされる（図１１参照）。まず、ＣＰＵ１２ａは、ステップＳ１６６にて、オンセット特徴量ＸＯ（ｔ０），ＸＯ（ｔ１）・・・をこの順にフィルタバンクＦＢＢに入力してフィルタ処理する。フィルタバンクＦＢＢは、拍周期ｂの値に応じてそれぞれ設けられた複数のコムフィルタＤｂからなる。コムフィルタＤｂ＝βは、フレームｔｉのオンセット特徴量ＸＯ（ｔｉ）を入力すると、前記入力したオンセット特徴量ＸＯ（ｔｉ）と「β」だけ先行するフレームｔｉ−βのオンセット特徴量ＸＯ（ｔｉ−β）に対する出力としてのデータＸＤｂ＝β（ｔｉ−β）とを所定の比率で加算してフレームｔｉのデータＸＤｂ＝β（ｔｉ）として出力する（図１０参照）。すなわち、コムフィルタＤｂ＝βは、データＸＤｂ＝βをフレーム数βに相当する時間だけ保持する保持手段としての遅延回路ｄｂ＝βを有する。上記のようにして、オンセット特徴量ＸＯの系列ＸＯ（ｔ）｛＝ＸＯ（ｔ０），ＸＯ（ｔ１）・・・｝をフィルタバンクＦＢＢに入力することにより、データＸＤｂの系列ＸＤｂ（ｔ）｛＝ＸＤｂ（ｔ０），ＸＤｂ（ｔ１）・・・｝が計算される。
【００３３】
つぎに、ＣＰＵ１２ａは、ステップＳ１６７にて、データＸＤｂの系列ＸＤｂ（ｔ）を時系列的に逆にしたデータ列をフィルタバンクＦＢＢに入力することにより、ＢＰＭ特徴量の系列ＸＢｂ（ｔ）｛＝ＸＢｂ（ｔ０），ＸＢｂ（ｔ１）・・・｝が得られる。これにより、オンセット特徴量ＸＯ（ｔ０），ＸＯ（ｔ１）・・・の位相とＢＰＭ特徴量ＸＢｂ（ｔ０），ＸＢｂ（ｔ１）・・・の位相のずれを「０」にすることができる。上記のようにして計算されたＢＰＭ特徴量ＸＢ（ｔｉ）を図１１に例示する。上記のように、ＢＰＭ特徴量ＸＢｂ（ｔｉ）は、オンセット特徴量ＸＯ（ｔｉ）と拍周期ｂの値に相当する時間（すなわち、フレーム数ｂ）だけ遅延させたＢＰＭ特徴量ＸＢｂ（ｔｉ―ｂ）とを所定の比率で加算して計算されるので、オンセット特徴量ＸＯ（ｔ０），ＸＯ（ｔ１）・・・が拍周期ｂの値に相当する時間間隔をおいてピークを有する場合、ＢＰＭ特徴量ＸＢｂ（ｔｉ）の値が大きくなる。楽曲のテンポは、１分間あたりの拍数で表されるから、拍周期ｂは１分間あたりの拍数の逆数に比例する。例えば、図１１に示す例では、拍周期ｂの値が「４」であるときのＢＰＭ特徴量ＸＢｂの値（ＢＰＭ特徴量ＸＢｂ＝４）が最も大きい。したがって、この例では拍が４フレームおきに存在する可能性が高い。本実施形態では、１フレームの時間の長さを１２５ｍｓとしたので、この場合の拍の間隔は０．５ｓである。すなわち、テンポは１２０ＢＰＭ（＝６０ｓ／０．５ｓ）である。
【００３４】
つぎに、ＣＰＵ１２ａは、ステップＳ１６８にて、特徴量計算処理を終了し、その処理を音響信号分析処理（メインルーチン）のステップＳ１７に進める。
【００３５】
ＣＰＵ１２ａは、ステップＳ１７にて、図１２に示す対数観測尤度計算プログラムをＲＯＭ１２ｂから読み出して実行する。対数観測尤度計算プログラムは、音響信号分析プログラムのサブルーチンである。
【００３６】
ＣＰＵ１２ａは、ステップＳ１７１にて対数観測尤度計算処理を開始する。そして、以下説明するように、オンセット特徴量ＸＯ（ｔｉ）の尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））、及びＢＰＭ特徴量ＸＢ（ｔｉ）の尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））を計算する。なお、上記の「Ｚｂ＝β，ｎ＝η（ｔｉ）」は、フレームｔｉにおいて、拍周期ｂの値が「β」であり、且つ次の拍までのフレーム数ｎの値が「η」である状態ｑｂ＝β，ｎ＝ηのみが生起していることを表わす。フレームｔｉにおいて状態ｑｂ＝β，ｎ＝ηと状態ｑｂ≠β，ｎ≠ηとが同時に生起することはない。したがって、尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ＝β，ｎ＝η（ｔｉ））は、フレームｔｉにおいて、拍周期ｂの値が「β」であり、且つ次の拍までのフレーム数ｎの値が「η」であるという条件のもとでオンセット特徴量ＸＯ（ｔｉ）が観測される確率を表わす。また、尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ＝β，ｎ＝η（ｔｉ））は、フレームｔｉにおいて、拍周期ｂの値が「β」であり、且つ次の拍までのフレーム数ｎの値が「η」であるという条件のもとでＢＰＭ特徴量ＸＢ（ｔｉ）が観測される確率を表わす。
【００３７】
まず、ＣＰＵ１２ａは、ステップＳ１７２にて、尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））を計算する。次の拍までのフレーム数ｎの値が「０」であるとき、オンセット特徴量ＸＯは、平均値が「３」であって、且つ分散が「１」である第１の正規分布に従って分布するものとする。すなわち、第１の正規分布の確率変数としてオンセット特徴量ＸＯ（ｔｉ）を代入した値を尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ，ｎ＝０（ｔｉ））として計算する。また、拍周期ｂの値が「β」であり、次の拍までのフレーム数ｎの値が「β／２」であるとき、オンセット特徴量ＸＯは、平均値が「１」であって、且つ分散が「１」である第２の正規分布に従って分布するものとする。すなわち、第２の正規分布の確率変数としてオンセット特徴量ＸＯ（ｔｉ）を代入した値を尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ＝β，ｎ＝β／２（ｔｉ））として計算する。また、次の拍までのフレーム数ｎの値が「０」及び「β／２」のうちのいずれの値とも異なるとき、オンセット特徴量ＸＯは、平均値が「０」であって、且つ分散が「１」である第３の正規分布に従って分布するものとする。すなわち、第３の正規分布の確率変数としてオンセット特徴量ＸＯ（ｔｉ）を代入した値を尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ，ｎ≠０，β／２（ｔｉ））として計算する。
【００３８】
オンセット特徴量ＸＯの系列が｛１０，２，０．５，５，１，０，３，４，２｝であるときの尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ＝６，ｎ（ｔｉ））の対数を計算した結果を、図１３に例示する。同図に示すように、オンセット特徴量ＸＯの値が大きいフレームｔｉほど、尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ，ｎ＝０（ｔｉ））が尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ，ｎ≠０（ｔｉ））に比べて大きい。このように、オンセット特徴量ＸＯの値が大きいフレームｔｉほど、フレーム数ｎの値が「０」であるときに拍が存在する可能性が高くなるように、確率モデル（第１乃至第３の正規分布、及びそれらのパラメータ（平均値及び分散））が設定されている。なお、第１乃至第３の正規分布のパラメータの値は、上記実施形態に限られない。これらのパラメータの値は、実験を繰り返して決定してもよいし、機械学習を用いて決定してもよい。なお、この例では、オンセット特徴量ＸＯの尤度Ｐを計算するための確率分布関数として正規分布を用いているが、確率分布関数として他の関数（例えば、ガンマ分布、ポアソン分布など）を用いても良い。
【００３９】
つぎに、ＣＰＵ１２ａは、ステップＳ１７３にて、尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））を計算する。尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ＝γ，ｎ（ｔｉ））は、図１４に示すテンプレートＴＰγ｛γ＝１，２・・・｝に対するＢＰＭ特徴量ＸＢ（ｔｉ）の適合度に相当する。具体的には、尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ＝γ，ｎ（ｔｉ））は、ＢＰＭ特徴量ＸＢ（ｔｉ）とテンプレートＴＰγ｛γ＝１，２・・・｝との内積に相当する（図１２のステップＳ１７３の演算式を参照）。なお、この演算式におけるκｂは、オンセット特徴量ＸＯに対するＢＰＭ特徴量ＸＢの重みを決定する係数である。つまり、κｂを大きく設定するほど、結果的に、後述する拍・テンポ同時推定処理においてＢＰＭ特徴量ＸＢが重視される。また、この演算式におけるＺ（κｂ）は、κｂに依存する正規化係数である。テンプレートＴＰγは、図１４に示すように、ＢＰＭ特徴量ＸＢ（ｔｉ）を構成するＢＰＭ特徴量ＸＢｂ（ｔｉ）にそれぞれ乗算される係数δγ，ｂからなる。テンプレートＴＰγは、その係数δγ，γが最大であり、係数δγ，２γ，係数δγ，３γ・・・，係数δγ，（「γ」の整数倍），・・・がそれぞれ極大となるように設定されている。すなわち、例えば、テンプレートＴＰγ＝２は、２フレームおきに拍が存在する楽曲に適合するように構成されている。なお、この例では、ＢＰＭ特徴量ＸＢの尤度Ｐを計算するためにテンプレートＴＰを用いているが、これに代えて確率分布関数（例えば、多項分布、ディリクレ分布、多次元正規分布、多次元ポアソン分布など）を用いても良い。
【００４０】
ＢＰＭ特徴量ＸＢ（ｔｉ）が図１１に示すような値であった場合に、図１４に示すテンプレートＴＰγ｛γ＝１，２・・・｝を用いて尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））を計算し、その対数を計算した結果を図１５に例示する。この例では、尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ＝４，ｎ（ｔｉ））が最も大きいので、ＢＰＭ特徴量ＸＢ（ｔｉ）は、テンプレートＴＰ４に最も適合している。
【００４１】
つぎに、ＣＰＵ１２ａは、ステップＳ１７４にて、尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））の対数と尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））の対数をそれぞれ加算し、その結果を対数観測尤度Ｌｂ，ｎ（ｔｉ）とする。なお、尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））と尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））とを積算した結果の対数を対数観測尤度Ｌｂ，ｎ（ｔｉ）としても同じ結果が得られる。つぎに、ＣＰＵ１２ａは、ステップＳ１７５にて、対数観測尤度計算処理を終了し、その処理を音響信号分析処理（メインルーチン）のステップＳ１８に進める。
【００４２】
つぎに、ＣＰＵ１２ａは、ステップＳ１８にて、図１６に示す拍・テンポ同時推定プログラムをＲＯＭ１２ｂから読み出して実行する。拍・テンポ同時推定プログラムは、音響信号分析プログラムのサブルーチンである。この拍・テンポ同時推定プログラムは、ビタビアルゴリズムを用いて最尤の状態の系列Ｑを計算するプログラムである。ここで、その概略について説明する。ＣＰＵ１２ａは、まず、フレームｔ０からフレームｔｉまでオンセット特徴量ＸＯ及びＢＰＭ特徴量ＸＢを観測したときにフレームｔｉの状態ｑｂ，ｎの尤度が最大となるような状態の系列を選択した場合の状態ｑｂ，ｎの尤度を尤度Ｃｂ，ｎ（ｔｉ）とするとともに、各状態ｑｂ，ｎに遷移する１つ前のフレームの状態（遷移元の状態）を状態Ｉｂ，ｎ（ｔｉ）として記憶する。つまり、遷移後の状態が状態ｑｂ＝βｅ，ｎ＝ηｅであって、遷移元の状態が状態ｑｂ＝βｓ，ｎ＝ηｓであるとき、状態Ｉｂ＝βｅ，ｎ＝ηｅ（ｔｉ）は、状態ｑｂ＝βｓ，ｎ＝ηｓである。ＣＰＵ１２ａは、上記のようにして尤度Ｃ及び状態Ｉをフレームｔｌａｓｔまで計算し、その結果を用いて最尤の状態の系列Ｑを選択する。
【００４３】
なお、後述する具体例では、その説明を簡単にするために、分析対象の楽曲の拍周期ｂの値が「３」、「４」及び「５」のうちのいずれかであるものとする。すなわち、対数観測尤度Ｌｂ，ｎ（ｔｉ）が図１７に例示するように計算された場合の拍・テンポ同時推定処理の手順を具体例として説明する。この例では、拍周期ｂの値が「３」、「４」及び「５」以外である状態の観測尤度が十分に小さいものとし、図１７乃至図１９では、拍周期ｂの値が「３」、「４」及び「５」以外である状態の観測尤度の図示を省略する。また、この例では、拍周期ｂの値が「βｓ」であり、且つフレーム数ｎの値が「ηｓ」である状態から、拍周期ｂの値が「βｅ」であり、且つフレーム数ｎの値が「ηｅ」である状態への対数遷移確率Ｔの値は、次のように設定されている。「ηｅ＝０」、「βｅ＝βｓ」、かつ「ηｅ＝βｅ−１」のとき、対数遷移確率Ｔの値は、「−０．２」である。また、「ηｓ＝０」、「βｅ＝βｓ＋１」、かつ「ηｅ＝βｅ−１」のとき、対数遷移確率Ｔの値は、「−０．６」である。また、「ηｓ＝０」、「βｅ＝βｓ−１」、かつ「ηｅ＝βｅ−１」のとき、対数遷移確率Ｔの値は、「−０．６」である。また、「ηｓ＞０」、「βｅ＝βｓ」、かつ「ηｅ＝ηｓ−１」のとき、対数遷移確率Ｔの値は、「０」である。上記以外の対数遷移確率Ｔの値は、「−∞」である。すなわち、フレーム数ｎの値が「０」である状態（ηｓ＝０）から次の状態へ遷移するとき、拍周期ｂの値は「１」だけ増減され得る。このとき、フレーム数ｎの値は、遷移後の拍周期ｂの値より「１」だけ小さい値に設定される。また、フレーム数ｎの値が「０」でない状態（ηｓ≠０）から次の状態へ遷移するとき、拍周期ｂの値は変更されず、フレーム数ｎの値が「１」だけ減少する。
【００４４】
以下、拍・テンポ同時推定処理について具体的に説明する。ＣＰＵ１２ａは、ステップＳ１８１にて拍・テンポ同時推定処理を開始する。つぎに、ユーザは、ステップＳ１８２にて、入力操作子１１を用いて、図１８に示すような、各状態ｑｂ，ｎに対応した、尤度Ｃの初期条件ＣＳｂ，ｎを入力する。なお、初期条件ＣＳｂ，ｎがＲＯＭ１２ｂに記憶されていて、ＣＰＵ１２ａがＲＯＭ１２ｂから初期条件ＣＳｂ，ｎを読み込むようにしてもよい。
【００４５】
つぎに、ＣＰＵ１２ａは、ステップＳ１８３にて、尤度Ｃｂ，ｎ（ｔｉ）及び状態Ｉｂ，ｎ（ｔｉ）を計算する。フレームｔ０において拍周期ｂの値が「βｅ」であって、フレーム数ｎの値が「ηｅ」である状態ｑｂ＝βｅ，ｎ＝ηｅの尤度Ｃｂ＝βｅ，ｎ＝ηｅ（ｔ０）は、初期条件ＣＳｂ＝βｅ，ｎ＝ηｅと対数観測尤度Ｌｂ＝βｅ，ｎ＝ηｅ（ｔ０）とを加算することにより計算される。
【００４６】
また、状態ｑｂ＝βｓ，ｎ＝ηｓから状態ｑｂ＝βｅ，ｎ＝ηｅに遷移したとき、尤度Ｃｂ＝βｅ，ｎ＝ηｅ（ｔｉ）｛ｉ＞０｝は次のように計算される。状態ｑｂ＝βｓ，ｎ＝ηｓのフレーム数ｎが「０」でないとき（すなわち、ηｓ≠０）、尤度Ｃｂ＝βｅ，ｎ＝ηｅ（ｔｉ）は、尤度Ｃｂ＝βｅ，ｎ＝ηｅ＋１（ｔｉ―１）と対数観測尤度Ｌｂ＝βｅ，ｎ＝ηｅ（ｔｉ）と対数遷移確率Ｔを加算して計算される。ただし、本実施形態では、遷移元の状態のフレーム数ｎが「０」でないときの対数遷移確率Ｔは「０」であるので、尤度Ｃｂ＝βｅ，ｎ＝ηｅ（ｔｉ）は、実質的には、尤度Ｃｂ＝βｅ，ｎ＝ηｅ＋１（ｔｉ―１）と対数観測尤度Ｌｂ＝βｅ，ｎ＝ηｅ（ｔｉ）とを加算することにより計算される（Ｃｂ＝βｅ，ｎ＝ηｅ（ｔｉ）＝Ｃｂ＝βｅ，ｎ＝ηｅ＋１（ｔｉ―１）＋Ｌｂ＝βｅ，ｎ＝ηｅ（ｔｉ））。また、この場合、状態Ｉｂ＝βｅ，ｎ＝ηｅ（ｔｉ）は、状態ｑβｅ，ηｅ＋１である。例えば、尤度Ｃが図１８に示すように計算された例では、尤度Ｃ４，１（ｔ２）の値は「２」であり、対数観測尤度Ｌ４，０（ｔ３）の値は「１」であるので、尤度Ｃ４，０（ｔ３）の値は「３」である。また、図１９に示すように、状態Ｉ４，０（ｔ３）は、状態ｑ４，１である。
【００４７】
また、状態ｑｂ＝βｓ，ｎ＝ηｓのフレーム数ｎが「０」のとき（ηｓ＝０）の尤度Ｃｂ＝βｅ，ｎ＝ηｅ（ｔｉ）は次のように計算される。この場合、状態の遷移に伴って拍周期ｂの値が増減され得る。そこで、まず、尤度Ｃβｅ−１，０（ｔｉ−１）、尤度Ｃβｅ，０（ｔｉ−１）、及び尤度Ｃβｅ＋１，０（ｔｉ−１）に対数遷移確率Ｔをそれぞれ加算し、そのうちの最大値に対数観測尤度Ｌｂ＝βｅ，ｎ＝ηｅ（ｔｉ）を加算した結果が尤度Ｃｂ＝βｅ，ｎ＝ηｅ（ｔｉ）である。また、状態Ｉｂ＝βｅ，ｎ＝ηｅ（ｔｉ）は、状態ｑβｅ−１，０、状態ｑβｅ，０、及び状態ｑβｅ＋１，０のうち、それらの尤度Ｃβｅ−１，０（ｔｉ−１）、尤度Ｃβｅ，０（ｔｉ−１）、及び尤度Ｃβｅ＋１，０（ｔｉ−１）に対数遷移確率Ｔをそれぞれ加算した値が最大となる状態ｑである。なお、厳密には、尤度Ｃｂ，ｎ（ｔｉ）は正規化される必要があるが、正規化されていなくても、拍点及びテンポの推移の推定に関しては、数理上同一の結果が得られる。
【００４８】
例えば、尤度Ｃ４，３（ｔ４）は、次のように計算される。遷移元の状態が状態ｑ３，０である場合、尤度Ｃ３，０（ｔ３）の値は「０．４」であり、対数遷移確率Ｔは「−０．６」であるので、尤度Ｃ３，０（ｔ３）と対数遷移確率Ｔとを加算した値は、「−０．２」である。また、遷移元の状態が状態ｑ４，０である場合、遷移元の尤度Ｃ４，０（ｔ３）の値は「３」であり、対数遷移確率Ｔは「−０．２」であるので、尤度Ｃ４，０（ｔ３）と対数遷移確率Ｔとを加算した値は、「２．８」である。また、遷移元の状態が状態ｑ５，０である場合、遷移元の尤度Ｃ５，０（ｔ３）の値は「１」であり、対数遷移確率Ｔは「−０．６」であるので、尤度Ｃ５，０（ｔ３）と対数遷移確率Ｔとを加算した値は、「０．４」である。したがって、尤度Ｃ４，０（ｔ３）に対数遷移確率Ｔを加算した値が最も大きい。また、対数観測尤度Ｌ４，３（ｔ４）の値は、「０」である。よって、尤度Ｃ４，３（ｔ４）の値は「２．８」（＝２．８＋０）であり、状態Ｉ４，３（ｔ４）は、状態ｑ４，０である。
【００４９】
上記のようにして、全てのフレームｔｉについて、全ての状態ｑｂ，ｎの尤度Ｃｂ，ｎ（ｔｉ）及び状態Ｉｂ，ｎ（ｔｉ）を計算し終えると、ＣＰＵ１２ａはステップＳ１８４にて、最尤の状態の系列Ｑ（＝｛ｑｍａｘ（ｔ０），ｑｍａｘ（ｔ１）・・・，ｑｍａｘ（ｔｌａｓｔ）｝）を次のようにして決定する。まず、ＣＰＵ１２ａは、フレームｔｌａｓｔにおける尤度Ｃｂ，ｎ（ｔｌａｓｔ）が最大である状態ｑｂ，ｎを、状態ｑｍａｘ（ｔｌａｓｔ）とする。ここで、状態ｑｍａｘ（ｔｌａｓｔ）の拍周期ｂの値を「βｍ」と表記し、フレーム数ｎの値を「ηｍ」と表記する。このとき、状態Ｉβｍ，ηｍ（ｔｌａｓｔ）がフレームｔｌａｓｔの１つ前のフレームｔｌａｓｔ−１の状態ｑｍａｘ（ｔｌａｓｔ−１）である。フレームｔｌａｓｔ−２、フレームｔｌａｓｔ−３、・・・の状態ｑｍａｘ（ｔｌａｓｔ−２）、状態ｑｍａｘ（ｔｌａｓｔ−３）・・・も状態ｑｍａｘ（ｔｌａｓｔ−１）と同様に決定される。すなわち、フレームｔｉ＋１の状態ｑｍａｘ（ｔｉ＋１）の拍周期ｂの値を「βｍ」と表記し、フレーム数ｎの値を「ηｍ」と表記したときの状態Ｉβｍ，ηｍ（ｔｉ＋１）がフレームｔｉ＋１の１つ前のフレームｔｉの状態ｑｍａｘ（ｔｉ）である。上記のようにして、ＣＰＵ１２ａは、フレームｔｌａｓｔ−１からフレームｔ０へ向かって順に状態ｑｍａｘを決定して、最尤の状態の系列Ｑを決定する。
【００５０】
例えば、図１８及び図１９に示す例では、フレームｔｌａｓｔ＝９においては、状態ｑ４，２の尤度Ｃ４，２（ｔｌａｓｔ＝９）が最大である。したがって、状態ｑｍａｘ（ｔｌａｓｔ＝９）は、状態ｑ４，２である。図１９によれば、状態Ｉ４，２（ｔ９）は状態ｑ４，３であるから、状態ｑｍａｘ（ｔ８）は状態ｑ４，３である。また、状態Ｉ４，３（ｔ８）は状態ｑ４，０であるから、状態ｑｍａｘ（ｔ７）は状態ｑ４，０である。状態ｑｍａｘ（ｔ６）乃至状態ｑｍａｘ（ｔ０）も状態ｑｍａｘ（ｔ８）及び状態ｑｍａｘ（ｔ７）と同様に決定する。このようにして図１８に矢印で示す最尤の状態の系列Ｑが決定される。この例では、拍の周期ｂの値はいずれのフレームｔｉにおいても「４」であると推定される。また、系列Ｑのうち、フレーム数ｎの値が「０」である状態ｑｍａｘ（ｔ１），ｑｍａｘ（ｔ５），ｑｍａｘ（ｔ８）に対応するフレームｔ１，ｔ５，ｔ８に拍が存在すると推定される。
【００５１】
つぎに、ＣＰＵ１２ａは、ステップＳ１８５にて、拍・テンポ同時推定処理を終了し、その処理を音響信号分析処理（メインルーチン）のステップＳ１９に進める。
【００５２】
ＣＰＵ１２ａは、ステップＳ１９にて、フレームｔｉごとに「ＢＰＭらしさ」、「観測に基づく確率」、「拍らしさ」、「拍が存在する確率」及び「拍が存在しない確率」を計算（図２０に示す演算式を参照）する。「ＢＰＭらしさ」は、フレームｔｉにおけるテンポの値が拍周期ｂに対応した値である確率を意味し、尤度Ｃｂ，ｎ（ｔｉ）を正規化するとともにフレーム数ｎについて周辺化することにより計算される。具体的には、拍周期ｂの値が「β」である場合の「ＢＰＭらしさ」は、フレームｔｉにおける全ての状態の尤度Ｃの合計に対する、拍周期ｂの値が「β」である状態の尤度Ｃの合計の割合である。また、「観測に基づく確率」は、観測値（すなわちオンセット特徴量ＸＯ）に基づいて計算された拍がフレームｔｉに存在する確率を意味する。具体的には、所定の基準値ＸＯｂａｓｅに対するオンセット特徴量ＸＯ（ｔｉ）の割合である。また、「拍らしさ」は、すべてのフレーム数ｎの値についてのオンセット特徴量ＸＯ（ｔｉ）の尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））を合算した値に対する尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ，０（ｔｉ））の割合である。また、「拍が存在する確率」及び「拍が存在しない確率」は、いずれも尤度Ｃｂ，ｎ（ｔｉ）を拍周期ｂについて周辺化することにより計算される。具体的には、「拍が存在する確率」は、フレームｔｉにおける全ての状態の尤度Ｃの合計に対する、フレーム数ｎの値が「０」である状態の尤度Ｃの合計の割合である。また、「拍が存在しない確率」は、フレームｔｉにおける全ての状態の尤度Ｃの合計に対する、フレーム数ｎの値が「０」でない状態の尤度Ｃの合計の割合である。
【００５３】
ＣＰＵ１２ａは、「ＢＰＭらしさ」、「観測に基づく確率」、「拍らしさ」、「拍が存在する確率」及び「拍が存在しない確率」を用いて、図２０に示す拍・テンポ情報リストを表示器１３に表示する。同リスト中の「テンポの推定値（ＢＰＭ）」の欄には、前記計算した「ＢＰＭらしさ」のうち最も確率の高い拍周期ｂに対応するテンポの値（ＢＰＭ）が表示される。また、前記決定した状態ｑｍａｘ（ｔｉ）のうちフレーム数ｎの値が「０」であるフレームの「拍の存在」の欄には「○」が表示され、その他のフレームの「拍の存在」の欄には「×」が表示される。また、ＣＰＵ１２ａは、テンポの推定値（ＢＰＭ）を用いて、図２１に示すようなテンポの推移を表わすグラフを表示器１３に表示する。図２１の例では、テンポの推移を棒グラフで表わしている。図１８及び図１９を用いて説明した例では、テンポの値が一定であるので図２１に示すような各フレームのテンポを表わすバーの高さが一定であるが、テンポが頻繁に変化する楽曲では、図２２に示すように、テンポの値に応じてバーの高さが異なる。これにより、ユーザは、テンポの推移を視覚的に認識することができる。また、ＣＰＵ１２ａは、前記計算した「拍が存在する確率」を用いて、図２３に示すような拍点を表わすグラフを表示器１３に表示する。
【００５４】
また、音響信号分析処理のステップＳ１３にて既存データを検索した結果、既存データが存在する場合には、ＣＰＵ１２ａは、ステップＳ１５にてＲＡＭ１２ｃに読み込んだ前回の分析結果に関する各種データを用いて、拍・テンポ情報リスト、テンポの推移を表わすグラフ、及び拍点を表わすグラフを表示器１３に表示する。
【００５５】
次に、ＣＰＵ１２ａは、ステップＳ２０にて、音響信号分析処理を終了するか否かを表すメッセージを表示器１３に表示して、ユーザからの指示を待つ。ユーザは入力操作子１１を用いて音響信号分析処理を終了するか、後述の拍・テンポ情報修正処理を実行するかのいずれかを指示する。例えば、マウスを用いて図示しないアイコンをクリックする。ユーザから音響信号分析処理を終了するよう指示された場合には、ＣＰＵ１２ａは「Ｙｅｓ」と判定してステップＳ２１にて尤度Ｃ、状態Ｉ、拍・テンポ情報リストなどの分析結果に関する各種データを楽曲のタイトルと関連付けて記憶装置１４に記憶して、ステップＳ２２にて音響信号分析処理を終了する。
【００５６】
一方、ステップＳ２０にて、音響信号分析処理を継続するように指示された場合には、ＣＰＵ１２ａは「Ｎｏ」と判定して、ステップＳ２３にて、テンポ情報修正処理を実行する。まず、ＣＰＵ１２ａは、ユーザが修正情報の入力を終了するまで待機する。ユーザは、入力操作子１１を用いて「ＢＰＭらしさ」、「拍が存在する確率」などの修正値を入力する。例えば、マウスを用いて修正するフレームを選択し、テンキーを用いて修正値を入力する。修正された項目の右側に配置された「Ｆ」の表示形態（例えば色）が変更され、その値が修正されたことが明示される。ユーザは、複数の項目について修正値を入力可能である。ユーザは修正値の入力を完了すると、入力操作子１１を用いて修正情報の入力を完了したことを指示する。例えば、マウスを用いて図示しない修正完了を表わすアイコンをクリックする。ＣＰＵ１２ａは、前記入力された修正値に応じて尤度Ｐ（ＸＯ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））及び尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ，ｎ（ｔｉ））のうちのいずれか一方又は両方を更新する。例えば、フレームｔｉにおける「拍が存在する確率」が高くなるように修正された場合であって、修正された値に関するフレーム数ｎの値が「ηｅ」であるときには、尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ，ｎ≠ηｅ（ｔｉ））を十分に小さい値に設定する。これにより、フレームｔｉでは、フレーム数ｎの値が「ηｅ」である確率が相対的に最も高くなる。また、例えば、フレームｔｉにおける「ＢＰＭらしさ」のうち、拍周期ｂの値が「βｅ」である確率が高くなるように修正された場合には、拍周期ｂの値が「βｅ」でない状態の尤度Ｐ（ＸＢ（ｔｉ）｜Ｚｂ≠βｅ，ｎ（ｔｉ））を十分に小さい値に設定する。これにより、フレームｔｉでは、拍周期ｂの値が「βｅ」である確率が相対的に最も高くなる。そして、ＣＰＵ１２ａは、拍・テンポ情報修正処理を終了して、その処理をステップＳ１８に進め、修正された対数観測尤度Ｌを用いて、拍・テンポ同時推定処理を再度実行する。
【００５７】
上記のように構成した音響信号分析装置１０によれば、拍点に関するオンセット特徴量ＸＯ及びテンポに関するＢＰＭ特徴量ＸＢを用いて計算された対数観測尤度Ｌの系列が最も尤もらしい確率モデルが選択され、楽曲における拍点及びテンポの推移が同時に推定される。したがって、上記従来技術とは異なり、拍点及びテンポのうちの一方の推定精度が低いために他方の推定精度も低くなるという事態が生じない。よって、従来技術に比べて楽曲における拍点及びテンポの推移の推定精度を向上させることができる。
【００５８】
また、本実施形態においては、フレーム数ｎの値が「０」である状態から、拍周期ｂの値が同じ状態又は拍周期ｂの値が「１」だけ異なる状態へのみ遷移可能に各状態間の遷移確率（対数遷移確率）が設定されている。これにより、テンポがフレーム間で急激に変化するような誤推定が防止される。したがって、楽曲として自然な拍点及びテンポの推移の推定結果を得ることができる。なお、テンポが急激に変化する楽曲に対しては、次の拍までのフレーム数ｎの値が「０」である状態から次の状態に遷移するとき、拍周期ｂの値が大きく異なる状態への遷移も可能なように各状態間の遷移確率（対数遷移確率）を設定すればよい。
【００５９】
また、拍・テンポ同時推定処理では、ビタビアルゴリズムを用いたので、他のアルゴリズム（例えば、「サンプリング法」、「前向き後向きアルゴリズム」など）を用いる場合に比べて計算量を削減できる。
【００６０】
また、ユーザにより入力された修正情報に基づいて対数観測尤度Ｌが修正され、修正された対数観測尤度Ｌに基づいて楽曲における拍点及びテンポの推移が再推定される。これにより、修正されたフレームの前後にそれぞれ位置する１つ又は複数のフレームの最尤の状態ｑｍａｘが再計算（再選択）される。したがって、修正されたフレーム及びその前後に位置する１つ又は複数のフレームに亘り、拍の間隔及びテンポが滑らかに変化するような推定結果が得られる。
【００６１】
上記のようにして推定された楽曲における拍点及びテンポの推移に関する情報は、例えば楽曲データの検索、伴奏を表わす伴奏データの検索などに利用される。また、分析対象とした楽曲に対する伴奏パートの自動生成、ハーモニーの自動付加などにも利用される。
【００６２】
さらに、本発明の実施にあたっては、上記実施形態に限定されるものではなく、本発明の目的を逸脱しない限りにおいて種々の変更が可能である。
【００６３】
例えば、上記実施形態では、観測値としてのオンセット特徴量ＸＯ及びＢＰＭ特徴量ＸＢが同時に観測される確率を表わす観測尤度の系列が最も尤もらしい確率モデルが選択される。しかし、確率モデルの選択基準は、上記実施形態に限られない。例えば、事後分布が最大となるような確率モデルを選択してもよい。
【００６４】
また、例えば、上記実施形態では、説明を簡単にするために、各フレームの長さを１２５ｍｓとしたが、より短く（例えば、５ｍｓ）してもよい。これによれば、拍点及びテンポの推定に関する分解能を向上させることができる。例えば、テンポを１ＢＰＭ刻みで推定できる。また、上記実施形態では、各フレームの長さを共通にしているが、各フレームの長さが異なっていてもよい。この場合であっても、オンセット特徴量ＸＯは、上記実施形態と同様にして計算できる。また、この場合、ＢＰＭ特徴量ＸＢの計算においては、コムフィルタの遅延量をフレームの長さに応じて変更すればよい。また、尤度Ｃの計算においては、各フレームの長さの最大公約数Ｆ（つまり、各フレームを構成するサンプル数の
最大公約数）を計算する。そして、フレームｔｉ（＝τ）の長さがＬ（τ）×Ｆと表わされたとき、状態ｑｂ，ｎ（ｎ≠０）から状態ｑｂ，ｎ−Ｌ（τ）へ遷移する確率を１００％とすればよい。
【００６５】
また、上記実施形態では、楽曲全体を分析対象としているが、楽曲の一部（例えば数小節）のみを分析対象としてもよい。この場合、入力した楽曲データのうち、分析対象とする部分を選択可能に構成するとよい。また、楽曲のうちの単一のパート（例えばリズムセクション）のみを分析対象としてもよい。
【００６６】
また、例えば、テンポの推定において、優先的に推定するテンポの範囲を指定可能に構成してもよい。具体的には、音響信号分析処理のステップＳ１２において、「Ｐｒｅｓｔｏ」、「Ｍｏｄｅｒａｔｏ」などのテンポを表わす用語を表示して、優先的に推定するテンポの範囲を選択可能に構成してもよい。例えば、「Ｐｒｅｓｔｏ」が選択された場合、ＢＰＭ＝１６０〜１９０の範囲以外の対数観測尤度Ｌを十分に小さく設定する。これにより、ＢＰＭ＝１６０〜１９０の範囲のテンポが優先的に推定される。これによれば、楽曲の大凡のテンポが既知である場合、テンポの推定精度を向上させることができる。
【００６７】
また、拍・テンポ情報修正処理（ステップＳ２３）では、ユーザは入力操作子１１を用いて修正内容を入力するように構成されている。これに代えて、又は加えて、外部インターフェース回路１５を介して接続された電子鍵盤楽器、電子打楽器などの操作子を用いて修正内容を入力可能に構成してもよい。例えば、ユーザが電子鍵盤楽器の鍵盤を数回打鍵すると、ＣＰＵ１２ａがその打鍵のタイミングからテンポを計算して、前記「ＢＰＭらしさ」の修正値として用いるように構成してもよい。
【００６８】
また、上記実施形態では、拍点及びテンポに関する修正値を何度でも入力可能に構成されている。しかし、例えば、「拍が存在する確率」の平均値が基準値（例えば８０％）に達した時点以降においては、拍点及びテンポに関する修正値を入力不可能としてもよい。
【００６９】
また、例えば、拍・テンポ情報修正処理（ステップＳ２３）において、ユーザによって指定されたフレームの拍・テンポ情報を入力された値に修正するとともに、そのフレームに近接するフレームの拍・テンポ情報を前記入力された値に応じて自動的に修正してもよい。例えば、連続する複数のフレームのテンポの推定値が同じ値であって、そのうちの１つのフレームのテンポの値が修正されたとき、前記複数のフレームのテンポの値を前記１つのフレームの修正値と同じ値に自動的に修正してもよい。
【００７０】
また、上記実施形態では、ステップＳ２３にて、ユーザが入力操作子１１を用いて修正値の入力を完了したことを指示すると、拍点及びテンポの同時推定が再び実行される。しかし、これに代えて、ユーザが少なくとも１つの修正値を入力した後、他の修正値が入力されないまま所定の時間（例えば１０秒）が経過したとき、自動的に拍点及びテンポの同時推定が再び実行されてもよい。
【００７１】
また、拍・テンポ情報リスト（図２０）の表示形態は、上記実施形態に限られない。例えば、上記実施形態では、「ＢＰＭらしさ」、「拍点らしさ」などを、確率（％）で表示しているが、これらを記号、文字列などを用いて表現してもよい。また、上記実施形態では、前記決定した状態ｑｍａｘ（ｔｉ）のうちフレーム数ｎの値が「０」であるフレームｔｉの「拍の存在」の欄には「○」が表示され、その他のフレームの「拍の存在」の欄には「×」が表示されるが、これに代えて、例えば、基準値（例えば８０％）以上であるとき「拍の存在」の欄に「○」が表示され、「拍点が存在する確率」が基準値未満であるとき、「拍の存在」の欄に「×」が表示されてもよい。また、この場合、複数の基準値を設けてもよい。例えば、第１基準値（＝８０％）と第２基準値（＝６０％）を設け、「拍点が存在する確率」が第１基準値以上であるとき、「拍の存在」の欄に「○」が表示され、第２基準値以上かつ第１基準値未満であるとき、「拍の存在」の欄に「△」が表示され、「拍点が存在する確率」が第２基準値未満であるとき、「拍の存在」の欄に「×」が表示されてもよい。また、テンポの推定値の欄には、「Ｐｒｅｓｔｏ」、「Ｍｏｄｅｒａｔｏ」などのテンポを表わす用語が表示されてもよい。
【符号の説明】
【００７２】
１０・・・音響信号分析装置、１１・・・入力操作子、ＸＯ・・・オンセット特徴量、ＸＢ・・・ＢＰＭ特徴量、ｂ・・・拍周期、ｎ・・・フレーム数、ＦＢＢ・・・フィルタバンク、ＴＰ・・・テンプレート
【図１】
2014178394.tif 000003
【図２】
2014178394.tif 000004
【図３】
2014178394.tif 000005
【図４】
2014178394.tif 000006
【図５】
2014178394.tif 000007
【図６】
2014178394.tif 000008
【図７】
2014178394.tif 000009
【図８】
2014178394.tif 000010
【図９】
2014178394.tif 000011
【図１０】
2014178394.tif 000012
【図１１】
2014178394.tif 000013
【図１２】
2014178394.tif 000014
【図１３】
2014178394.tif 000015
【図１４】
2014178394.tif 000016
【図１５】
2014178394.tif 000017
【図１６】
2014178394.tif 000018
【図１７】
2014178394.tif 000019
【図１８】
2014178394.tif 000020
【図１９】
2014178394.tif 000021
【図２０】
2014178394.tif 000022
【図２１】
2014178394.tif 000023
【図２２】
2014178394.tif 000024
【図２３】
2014178394.tif 000025