JP6123995B2 - 音響信号分析装置及び音響信号分析プログラム - Google Patents

Description

本発明は、楽曲を表わす音響信号を入力して、前記入力した楽曲における拍点（拍のタイミング）及びテンポを検出する音響信号分析装置に関する。

従来から、例えば、下記特許文献１に示されているように、楽曲を表わす音響信号を入力して楽曲における拍点及びテンポを検出する音響信号分析装置は知られている。

特開２００９−２６５４９３号公報

上記特許文献１の音響信号分析装置は、まず、音響信号の強度（振幅）の変化に基づいて、拍点の候補として拍指標系列を計算する。次に、拍指標系列の計算結果に基づいて、楽曲のテンポを検出する。したがって、拍指標系列の検出精度が低い場合、テンポの検出精度も低くなる。

本発明は上記問題に対処するためになされたもので、その目的は、楽曲における拍点及びテンポの推移を高精度に検出できる音響信号分析装置を提供することにある。なお、下記本発明の各構成要件の記載においては、本発明の理解を容易にするために、実施形態の対応箇所の符号を括弧内に記載しているが、本発明の各構成要件は、実施形態の符号によって示された対応箇所の構成に限定解釈されるべきものではない。

上記目的を達成するために、本発明の特徴は、楽曲を表わす音響信号を入力する音響信号入力手段（Ｓ１２）と、前記楽曲の各区間における拍の存在に関する特徴を表わす第１特徴量（ＸＯ）及びテンポに関する特徴を表わす第２特徴量（ＸＢ）を計算する特徴量計算手段（Ｓ１６５，Ｓ１６７）と、前記楽曲の各区間における拍の存在に関する物理量（ｎ）及びテンポに関する物理量（ｂ）の組み合わせにより分類された状態（ｑ_ｂ，ｎ）の系列として記述された複数の確率モデルのうち、前記第１特徴量及び前記第２特徴量が前記楽曲の各区間において同時に観測される確率を表わす観測尤度（Ｌ）の系列が所定の基準を満たす確率モデルを選択することにより、前記楽曲における拍点及びテンポの推移を同時に推定する推定手段（Ｓ１７，Ｓ１８）と、を備えたことにある。

この場合、推定手段は、拍の存在に関する物理量に応じて規定された確率分布関数の確率変数として第１特徴量を代入して計算された確率を第１特徴量が観測される確率として出力する第１確率出力手段（Ｓ１７２）と、テンポに関する物理量の大きさに応じて規定された確率分布関数の確率変数として第２特徴量を代入して計算された確率を第２特徴量が観測される確率として出力する第２確率出力手段（Ｓ１７３）と、を備えるとよい。

この場合、第１確率出力手段は、前記拍の存在に関する物理量に応じて規定された正規分布の確率変数として第１特徴量を代入して計算された確率を前記第１特徴量が観測される確率として出力するとよい。

また、この場合、第２確率出力手段は、前記テンポに関する物理量に応じてそれぞれ設定された複数のテンプレートに対する第２特徴量の適合度を前記第２特徴量が観測される確率として出力するとよい。

また、この場合、前記楽曲の各区間は、前記入力した音響信号を所定の時間間隔をおいて区切ることにより形成された各フレームにそれぞれ相当し、特徴量計算手段は、フレームごとに振幅スペクトル（Ａ）を計算し、前記振幅スペクトルに周波数帯域（ｗ_ｋ）が互いに異なる複数の窓関数（ＢＰＦ）を掛けて前記周波数帯域ごとの振幅スペクトル（Ｍ）を生成し、フレーム間における前記周波数帯域ごとの振幅スペクトルの変化に基づいて計算した値を第１特徴量として出力する第１特徴量計算手段（Ｓ１６５）と、各フレームに対応する値を入力するごとに１つの値を出力するフィルタであって、出力した値を所定の時間だけ保持する保持手段（ｄ_ｂ）を備え、前記入力した値と前記所定の時間だけ保持された値とを所定の比率で加算して出力するフィルタ（ＦＢＢ）を備え、第１特徴量の系列をフィルタに入力して得られたデータ列の時系列を逆にしたデータ列を前記フィルタに再び入力して得られたデータ列を第２特徴量の系列として出力する第２特徴量計算手段（Ｓ１６７）と、を備えるとよい。

上記のように構成した音響信号分析装置によれば、拍の存在に関する特徴を表わす第１特徴量及びテンポに関する特徴を表わす第２特徴量を用いて計算された観測尤度の系列が所定の基準を満たす確率モデル（例えば、最も尤もらしい確率モデル、事後分布が最大となる確率モデルなど）が選択され、楽曲における拍点及びテンポの推移が同時に推定される。したがって、上記従来技術とは異なり、拍点及びテンポのうちの一方の推定精度が低いために他方の推定精度も低くなるという事態が生じない。よって、従来技術に比べて拍点及びテンポの推移の推定精度を向上させることができる。

また、本発明の他の特徴は、前記楽曲における拍点及びテンポの推移のいずれか一方又は両方の修正内容を表わす修正情報を入力する入力手段（１１，Ｓ２３）と、前記入力した修正情報に応じて前記観測尤度を修正する観測尤度修正手段（Ｓ２３）と、前記複数の確率モデルのうち、前記修正された観測尤度の系列が前記所定の基準を満たす確率モデルを、推定手段を用いて選択することにより、前記楽曲における拍点及びテンポの推移を同時に再推定する再推定手段（Ｓ２３，Ｓ１８）と、を備えたことにある。

これによれば、ユーザにより入力された修正情報に基づいて観測尤度が修正され、修正された観測尤度に基づいて楽曲における拍点及びテンポの推移が再推定される。つまり、修正されたフレームの前後にそれぞれ位置する１つ又は複数のフレームの状態が再計算（再選択）される。したがって、修正されたフレーム及びその前後に位置する１つ又は複数のフレームに亘り、拍の間隔（すなわち、テンポ）が滑らかに変化するような推定結果が得られる。

さらに、本発明の実施にあたっては、音響信号分析装置の発明に限定されることなく、同装置に適用されるコンピュータプログラムの発明としても実施し得るものである。

本発明の一実施形態に係る音響信号分析装置の全体構成を示すブロック図である。 確率モデルの概念図である。 音響信号分析プログラムを表わすフローチャートである。 特徴量計算プログラムを表わすフローチャートである。 分析対象の音響信号の波形を表わすグラフである。 １つのフレームを短時間フーリエ変換した音響スペクトル図である。 バンドパスフィルタの特性図である。 各周波数帯域の振幅の時間変化を示すグラフである。 オンセット特徴量の時間変化を示すグラフである。 コムフィルタのブロック図である。 ＢＰＭ特徴量の計算結果を示すグラフである。 対数観測尤度計算プログラムを表わすフローチャートである。 オンセット特徴量の観測尤度の計算結果を示す表である。 テンプレートの構成を示す表である。 ＢＰＭ特徴量の観測尤度の計算結果を示す表である。 拍・テンポ同時推定プログラムを表わすフローチャート 対数観測尤度の計算結果を示す表である。 先頭のフレームから各フレームまでオンセット特徴量及びＢＰＭ特徴量を観測したときに前記各フレームの各状態の尤度が最大となるような状態の系列を選択した場合の前記各状態の尤度の計算結果を示す表である。 遷移元の状態の計算結果を示す表である。 拍・テンポ情報リストの概略を示す概略図である。 テンポの推移の一例を示すグラフである。 テンポの推移の他の例を示すグラフである。 拍点を示すグラフである。

本発明の一実施形態に係る音響信号分析装置１０について説明する。音響信号分析装置１０は、以下説明するように、楽曲を表わす音響信号を入力して、その楽曲における拍点及びテンポの推移を検出する。音響信号分析装置１０は、図１に示すように、入力操作子１１、コンピュータ部１２、表示器１３、記憶装置１４、外部インターフェース回路１５及びサウンドシステム１６を備えており、これらがバスＢＳを介して接続されている。

入力操作子１１は、オン・オフ操作に対応したスイッチ（例えば数値を入力するためのテンキー）、回転操作に対応したボリューム又はロータリーエンコーダ、スライド操作に対応したボリューム又はリニアエンコーダ、マウス、タッチパネルなどから構成される。これらの操作子は、演奏者の手によって操作されて、分析対象の楽曲の選択、音響信号の分析開始又は停止、楽曲の再生又は停止（後述するサウンドシステム１６からの出力又は停止）、音響信号の分析に関する各種パラメータの設定などに用いられる。入力操作子１１を操作すると、その操作内容を表す操作情報が、バスＢＳを介して、後述するコンピュータ部１２に供給される。

コンピュータ部１２は、バスＢＳにそれぞれ接続されたＣＰＵ１２ａ、ＲＯＭ１２ｂ及びＲＡＭ１２ｃからなる。ＣＰＵ１２ａは、詳しくは後述する音響信号分析プログラム及びそのサブルーチンをＲＯＭ１２ｂから読み出して実行する。ＲＯＭ１２ｂには、音響信号分析プログラム及びそのサブルーチンに加えて、初期設定パラメータ、表示器１３に表示される画像を表わす表示データを生成するための図形データ及び文字データなどの各種データが記憶されている。ＲＡＭ１２ｃには、音響信号分析プログラムの実行時に必要なデータが一時的に記憶される。

表示器１３は、液晶ディスプレイ（ＬＣＤ）によって構成される。コンピュータ部１２は、図形データ、文字データなどを用いて表示すべき内容を表わす表示データを生成して表示器１３に供給する。表示器１３は、コンピュータ部１２から供給された表示データに基づいて画像を表示する。例えば分析対象の楽曲の選択時には、楽曲のタイトルリストが表示される。また、例えば分析終了時には、拍点及びテンポの推移を表わす拍・テンポ情報リスト及びそのグラフ（図２０乃至図２３参照）が表示される。

また、記憶装置１４は、ＨＤＤ、ＦＤＤ、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤなどの大容量の不揮発性記録媒体と、同各記録媒体に対応するドライブユニットから構成されている。記憶装置１４には、複数の楽曲をそれぞれ表わす複数の楽曲データが記憶されている。楽曲データは、楽曲を所定のサンプリング周期（例えば４４．１ｋＨｚ）でサンプリングして得られた複数のサンプル値からなり、各サンプル値が記憶装置１４における連続するアドレスに順に記録されている。楽曲のタイトルを表わすタイトル情報、楽曲データの容量を表わすデータサイズ情報なども楽曲データに含まれている。楽曲データは予め記憶装置１４に記憶されていてもよいし、後述する外部インターフェース回路１５を介して外部から取り込んでもよい。記憶装置１４に記憶されている楽曲データは、ＣＰＵ１２ａによって読み込まれ、楽曲における拍点及びテンポの推移が分析される。

外部インターフェース回路１５は、音響信号分析装置１０を電子音楽装置、パーソナルコンピュータなどの外部機器に接続可能とする接続端子を備えている。音響信号分析装置１０は、外部インターフェース回路１５を介して、ＬＡＮ（Ｌｏｃａｌ Ａｒｅａ Ｎｅｔｗｏｒｋ）、インターネットなどの通信ネットワークにも接続可能である。

サウンドシステム１６は、楽曲データをアナログ音信号に変換するＤ／Ａ変換器、変換したアナログ音信号を増幅するアンプ、及び増幅されたアナログ音信号を音響信号に変換して出力する左右一対のスピーカを備えている。ユーザが入力操作子１１を用いて分析対象の楽曲の再生を指示すると、ＣＰＵ１２ａは、分析対象の楽曲データをサウンドシステム１６に供給する。これにより、ユーザは分析対象の楽曲を試聴できる。

つぎに、上記のように構成した音響信号分析装置１０の動作について説明する。まず、その概略について説明する。分析対象の楽曲は複数のフレームｔ_ｉ｛ｉ＝０，１，・・・，ｌａｓｔ｝に分割される。そして、拍の存在に関する特徴を表すオンセット特徴量ＸＯ及びテンポに関する特徴を表すＢＰＭ特徴量ＸＢがフレームｔ_ｉごとに計算される。各フレームｔ_ｉにおける拍周期ｂの値（テンポの逆数に比例する値）及び次の拍までのフレーム数ｎの値の組み合わせに応じて分類された状態ｑ_ｂ，ｎの系列として記述された確率モデル（隠れマルコフモデル）のうち、観測値としてのオンセット特徴量ＸＯ及びＢＰＭ特徴量ＸＢが同時に観測される確率を表わす観測尤度の系列が最も尤もらしい確率モデルが選択される（図２参照）。これにより、分析対象の楽曲における拍点及びテンポの推移が検出される。なお、拍周期ｂは、フレームの数によって表わされる。したがって、拍周期ｂの値は「１≦ｂ≦ｂ_ｍａｘ」を満たす整数であり、拍周期ｂの値が「β」である状態では、フレーム数ｎの値は「０≦ｎ＜β」を満たす整数である。

つぎに、音響信号分析装置１０の動作について具体的に説明する。ユーザが音響信号分析装置１０の図示しない電源スイッチをオンにすると、ＣＰＵ１２ａは、図３に示す音響信号分析プログラムをＲＯＭ１２ｂから読み出して実行する。

ＣＰＵ１２ａは、ステップＳ１０にて音響信号分析処理を開始し、ステップＳ１１にて、記憶装置１４に記憶されている複数の楽曲データにそれぞれ含まれるタイトル情報を読み込んで、楽曲のタイトルをリスト形式で表示器１３に表示する。ユーザは、入力操作子１１を用いて、表示器１３に表示された楽曲の中から分析対象の楽曲データを選択する。なお、ステップＳ１１にて分析対象の楽曲データを選択する際、選択しようとする楽曲データが表す楽曲の一部又は全部を再生して楽曲データの内容を確認できるように構成してもよい。

つぎに、ＣＰＵ１２ａは、ステップＳ１２にて、音響信号分析のための初期設定を実行する。具体的には、前記選択された楽曲データのデータサイズ情報に応じた記憶領域をＲＡＭ１２ｃ内に確保し、前記確保した記憶領域に前記選択された楽曲データを読み込む。また、分析結果を表す拍・テンポ情報リスト、オンセット特徴量ＸＯ、ＢＰＭ特徴量ＸＢなどを一時的に記憶する領域をＲＡＭ１２ｃ内に確保する。

詳しくは後述するが、本プログラムによる分析結果は、記憶装置１４に保存される（ステップＳ２１）。前記選択された楽曲が本プログラムによって過去に分析されたことがあれば、記憶装置１４にその分析結果が保存されている。そこで、ＣＰＵ１２ａは、ステップＳ１３にて、前記選択された楽曲の分析に関する既存のデータ（以下、単に既存データと呼ぶ）を検索する。既存データが有れば、ＣＰＵ１２ａは、ステップＳ１４にて「Ｙｅｓ」と判定して、ステップＳ１５にて既存データをＲＡＭ１２ｃに読み込み、後述するステップＳ１９に処理を進める。一方、既存データが無ければ、ＣＰＵ１２ａは、ステップＳ１４にて「Ｎｏ」と判定して、その処理をステップＳ１６に進める。

ＣＰＵ１２ａは、ステップＳ１６にて、図４に示す特徴量計算プログラムをＲＯＭ１２ｂから読み出して実行する。特徴量計算プログラムは、音響信号分析プログラムのサブルーチンである。

ＣＰＵ１２ａは、ステップＳ１６１にて特徴量計算処理を開始する。つぎに、ＣＰＵ１２ａは、ステップＳ１６２にて、図５に示すように、前記選択された楽曲を所定の時間間隔をおいて区切り、複数のフレームｔ_ｉ｛ｉ＝０，１，・・・，ｌａｓｔ｝に分割する。各フレームの長さは共通である。説明を簡単にするために、本実施形態では各フレームの長さを１２５ｍｓとする。上記のように、各楽曲のサンプリング周波数は４４．１ｋＨｚであるので、各フレームは、約５０００個のサンプル値から構成されている。そして、以下説明するように、フレームごとに、オンセット特徴量ＸＯ及びＢＰＭ（ｂｅａｔｓ ｐｅｒ ｍｉｎｕｔｅ（１分間あたりの拍数））特徴量ＸＢを計算する。

つぎに、ＣＰＵ１２ａは、ステップＳ１６３にて、フレームごとに短時間フーリエ変換を実行して、図６に示すように、各周波数ビンｆ_ｊ｛ｊ＝１，２・・・｝の振幅Ａ（ｆ_ｊ，ｔ_ｉ）を計算する。そして、ＣＰＵ１２ａは、ステップＳ１６４にて、振幅Ａ（ｆ_１，ｔ_ｉ），Ａ（ｆ_２，ｔ_ｉ）・・・を周波数ビンｆ_ｊごとに設けられたフィルタバンクＦＢＯ_ｊによってフィルタ処理することにより、所定の周波数帯域ｗ_ｋ｛ｋ＝１，２，・・・｝の振幅Ｍ（ｗ_ｋ，ｔ_ｉ）を計算する。周波数ビンｆ_ｊ用のフィルタバンクＦＢＯ_ｊは、図７に示すように、通過帯域の中心周波数が互いに異なる複数のバンドパスフィルタＢＰＦ（ｗ_ｋ，ｆ_ｊ）からなる。フィルタバンクＦＢＯ_ｊを構成する各バンドパスフィルタＢＰＦ（ｗ_ｋ，ｆ_ｊ）の中心周波数は、対数周波数軸上で等間隔であり、かつ各バンドパスフィルタＢＰＦ（ｗ_ｋ，ｆ_ｊ）の通過帯域幅は、対数周波数軸上で共通である。各バンドパスフィルタＢＰＦ（ｗ_ｋ，ｆ_ｊ）は、通過帯域の中心周波数から通過帯域の下限周波数側及び上限周波数側へ向かうに従って徐々にゲインがそれぞれ小さくなるように構成されている。ＣＰＵ１２ａは、図４のステップＳ１６４に示すように、周波数ビンｆ_ｊごとに振幅Ａ（ｆ_ｊ，ｔ_ｉ）とバンドパスフィルタＢＰＦ（ｗ_ｋ，ｆ_ｊ）のゲインとを積算する。そして、前記周波数ビンｆ_ｊごとに計算した積算結果を全ての周波数ビンｆ_ｊについて合算して振幅Ｍ（ｗ_ｋ，ｔ_ｉ）とする。上記のようにして計算された振幅Ｍの系列を図８に例示する。

つぎに、ＣＰＵ１２ａは、ステップＳ１６５にて、振幅Ｍの時間変化に基づいてフレームｔ_ｉのオンセット特徴量ＸＯ（ｔ_ｉ）を計算する。具体的には、図４のステップＳ１６５に示すように、周波数帯域ｗ_ｋごとに、フレームｔ_ｉ−１からフレームｔ_ｉへの振幅Ｍの増加量Ｒ（ｗ_ｋ，ｔ_ｉ）を計算する。ただし、フレームｔ_ｉ―１の振幅Ｍ（ｗ_ｋ，ｔ_ｉ−１）とフレームｔ_ｉの振幅Ｍ（ｗ_ｋ，ｔ_ｉ）とが同じである場合、又はフレームｔ_ｉの振幅Ｍ（ｗ_ｋ，ｔ_ｉ）がフレームｔ_ｉ―１の振幅Ｍ（ｗ_ｋ，ｔ_ｉ−１）よりも小さい場合は、増加量Ｒ（ｗ_ｋ，ｔ_ｉ）は「０」とする。そして、周波数帯域ｗ_ｋごとに計算した増加量Ｒ（ｗ_ｋ，ｔ_ｉ）を全ての周波数帯域ｗ_１，ｗ_２，・・・について合算してオンセット特徴量ＸＯ（ｔ_ｉ）とする。上記のようにして計算されたオンセット特徴量ＸＯの系列を図９に例示する。一般に、楽曲においては、拍が存在する部分の音量が大きい。したがって、オンセット特徴量ＸＯ（ｔ_ｉ）が大きいほど、フレームｔ_ｉに拍が存在する可能性が高い。

つぎに、ＣＰＵ１２ａは、オンセット特徴量ＸＯ（ｔ_０），ＸＯ（ｔ_１）・・・を用いて、ＢＰＭ特徴量ＸＢをフレームｔ_ｉごとに計算する。なお、フレームｔ_ｉのＢＰＭ特徴量ＸＢ（ｔ_ｉ）は、拍周期ｂごとに計算されたＢＰＭ特徴量ＸＢ_{ｂ＝１，２・・・}（ｔ_ｉ）の集合として表わされる（図１１参照）。まず、ＣＰＵ１２ａは、ステップＳ１６６にて、オンセット特徴量ＸＯ（ｔ_０），ＸＯ（ｔ_１）・・・をこの順にフィルタバンクＦＢＢに入力してフィルタ処理する。フィルタバンクＦＢＢは、拍周期ｂの値に応じてそれぞれ設けられた複数のコムフィルタＤ_ｂからなる。コムフィルタＤ_ｂ＝βは、フレームｔ_ｉのオンセット特徴量ＸＯ（ｔ_ｉ）を入力すると、前記入力したオンセット特徴量ＸＯ（ｔ_ｉ）と「β」だけ先行するフレームｔ_ｉ−βのオンセット特徴量ＸＯ（ｔ_ｉ−β）に対する出力としてのデータＸＤ_ｂ＝β（ｔ_ｉ−β）とを所定の比率で加算してフレームｔ_ｉのデータＸＤ_ｂ＝β（ｔ_ｉ）として出力する（図１０参照）。すなわち、コムフィルタＤ_ｂ＝βは、データＸＤ_ｂ＝βをフレーム数βに相当する時間だけ保持する保持手段としての遅延回路ｄ_ｂ＝βを有する。上記のようにして、オンセット特徴量ＸＯの系列ＸＯ（ｔ）｛＝ＸＯ（ｔ_０），ＸＯ（ｔ_１）・・・｝をフィルタバンクＦＢＢに入力することにより、データＸＤ_ｂの系列ＸＤ_ｂ（ｔ）｛＝ＸＤ_ｂ（ｔ_０），ＸＤ_ｂ（ｔ_１）・・・｝が計算される。

つぎに、ＣＰＵ１２ａは、ステップＳ１６７にて、データＸＤ_ｂの系列ＸＤ_ｂ（ｔ）を時系列的に逆にしたデータ列をフィルタバンクＦＢＢに入力することにより、ＢＰＭ特徴量の系列ＸＢ_ｂ（ｔ）｛＝ＸＢ_ｂ（ｔ_０），ＸＢ_ｂ（ｔ_１）・・・｝が得られる。これにより、オンセット特徴量ＸＯ（ｔ_０），ＸＯ（ｔ_１）・・・の位相とＢＰＭ特徴量ＸＢ_ｂ（ｔ_０），ＸＢ_ｂ（ｔ_１）・・・の位相のずれを「０」にすることができる。上記のようにして計算されたＢＰＭ特徴量ＸＢ（ｔ_ｉ）を図１１に例示する。上記のように、ＢＰＭ特徴量ＸＢ_ｂ（ｔ_ｉ）は、オンセット特徴量ＸＯ（ｔ_ｉ）と拍周期ｂの値に相当する時間（すなわち、フレーム数ｂ）だけ遅延させたＢＰＭ特徴量ＸＢ_ｂ（ｔ_ｉ―ｂ）とを所定の比率で加算して計算されるので、オンセット特徴量ＸＯ（ｔ_０），ＸＯ（ｔ_１）・・・が拍周期ｂの値に相当する時間間隔をおいてピークを有する場合、ＢＰＭ特徴量ＸＢ_ｂ（ｔ_ｉ）の値が大きくなる。楽曲のテンポは、１分間あたりの拍数で表されるから、拍周期ｂは１分間あたりの拍数の逆数に比例する。例えば、図１１に示す例では、拍周期ｂの値が「４」であるときのＢＰＭ特徴量ＸＢ_ｂの値（ＢＰＭ特徴量ＸＢ_ｂ＝４）が最も大きい。したがって、この例では拍が４フレームおきに存在する可能性が高い。本実施形態では、１フレームの時間の長さを１２５ｍｓとしたので、この場合の拍の間隔は０．５ｓである。すなわち、テンポは１２０ＢＰＭ（＝６０ｓ／０．５ｓ）である。

つぎに、ＣＰＵ１２ａは、ステップＳ１６８にて、特徴量計算処理を終了し、その処理を音響信号分析処理（メインルーチン）のステップＳ１７に進める。

ＣＰＵ１２ａは、ステップＳ１７にて、図１２に示す対数観測尤度計算プログラムをＲＯＭ１２ｂから読み出して実行する。対数観測尤度計算プログラムは、音響信号分析プログラムのサブルーチンである。

ＣＰＵ１２ａは、ステップＳ１７１にて対数観測尤度計算処理を開始する。そして、以下説明するように、オンセット特徴量ＸＯ（ｔ_ｉ）の尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））、及びＢＰＭ特徴量ＸＢ（ｔ_ｉ）の尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））を計算する。なお、上記の「Ｚ_{ｂ＝β，ｎ＝η}（ｔ_ｉ）」は、フレームｔ_ｉにおいて、拍周期ｂの値が「β」であり、且つ次の拍までのフレーム数ｎの値が「η」である状態ｑ_{ｂ＝β，ｎ＝η}のみが生起していることを表わす。フレームｔ_ｉにおいて状態ｑ_{ｂ＝β，ｎ＝η}と状態ｑ_{ｂ≠β，ｎ≠η}とが同時に生起することはない。したがって、尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_{ｂ＝β，ｎ＝η}（ｔ_ｉ））は、フレームｔ_ｉにおいて、拍周期ｂの値が「β」であり、且つ次の拍までのフレーム数ｎの値が「η」であるという条件のもとでオンセット特徴量ＸＯ（ｔ_ｉ）が観測される確率を表わす。また、尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_{ｂ＝β，ｎ＝η}（ｔ_ｉ））は、フレームｔ_ｉにおいて、拍周期ｂの値が「β」であり、且つ次の拍までのフレーム数ｎの値が「η」であるという条件のもとでＢＰＭ特徴量ＸＢ（ｔ_ｉ）が観測される確率を表わす。

まず、ＣＰＵ１２ａは、ステップＳ１７２にて、尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））を計算する。次の拍までのフレーム数ｎの値が「０」であるとき、オンセット特徴量ＸＯは、平均値が「３」であって、且つ分散が「１」である第１の正規分布に従って分布するものとする。すなわち、第１の正規分布の確率変数としてオンセット特徴量ＸＯ（ｔ_ｉ）を代入した値を尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_{ｂ，ｎ＝０}（ｔ_ｉ））として計算する。また、拍周期ｂの値が「β」であり、次の拍までのフレーム数ｎの値が「β／２」であるとき、オンセット特徴量ＸＯは、平均値が「１」であって、且つ分散が「１」である第２の正規分布に従って分布するものとする。すなわち、第２の正規分布の確率変数としてオンセット特徴量ＸＯ（ｔ_ｉ）を代入した値を尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_{ｂ＝β，ｎ＝β／２}（ｔ_ｉ））として計算する。また、次の拍までのフレーム数ｎの値が「０」及び「β／２」のうちのいずれの値とも異なるとき、オンセット特徴量ＸＯは、平均値が「０」であって、且つ分散が「１」である第３の正規分布に従って分布するものとする。すなわち、第３の正規分布の確率変数としてオンセット特徴量ＸＯ（ｔ_ｉ）を代入した値を尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_{ｂ，ｎ≠０，β／２}（ｔ_ｉ））として計算する。

オンセット特徴量ＸＯの系列が｛１０，２，０．５，５，１，０，３，４，２｝であるときの尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_{ｂ＝６，ｎ}（ｔ_ｉ））の対数を計算した結果を、図１３に例示する。同図に示すように、オンセット特徴量ＸＯの値が大きいフレームｔ_ｉほど、尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_{ｂ，ｎ＝０}（ｔ_ｉ））が尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_{ｂ，ｎ≠０}（ｔ_ｉ））に比べて大きい。このように、オンセット特徴量ＸＯの値が大きいフレームｔ_ｉほど、フレーム数ｎの値が「０」であるときに拍が存在する可能性が高くなるように、確率モデル（第１乃至第３の正規分布、及びそれらのパラメータ（平均値及び分散））が設定されている。なお、第１乃至第３の正規分布のパラメータの値は、上記実施形態に限られない。これらのパラメータの値は、実験を繰り返して決定してもよいし、機械学習を用いて決定してもよい。なお、この例では、オンセット特徴量ＸＯの尤度Ｐを計算するための確率分布関数として正規分布を用いているが、確率分布関数として他の関数（例えば、ガンマ分布、ポアソン分布など）を用いても良い。

つぎに、ＣＰＵ１２ａは、ステップＳ１７３にて、尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））を計算する。尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_{ｂ＝γ，ｎ}（ｔ_ｉ））は、図１４に示すテンプレートＴＰ_γ｛γ＝１，２・・・｝に対するＢＰＭ特徴量ＸＢ（ｔ_ｉ）の適合度に相当する。具体的には、尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_{ｂ＝γ，ｎ}（ｔ_ｉ））は、ＢＰＭ特徴量ＸＢ（ｔ_ｉ）とテンプレートＴＰ_γ｛γ＝１，２・・・｝との内積に相当する（図１２のステップＳ１７３の演算式を参照）。なお、この演算式におけるκ_ｂは、オンセット特徴量ＸＯに対するＢＰＭ特徴量ＸＢの重みを決定する係数である。つまり、κ_ｂを大きく設定するほど、結果的に、後述する拍・テンポ同時推定処理においてＢＰＭ特徴量ＸＢが重視される。また、この演算式におけるＺ（κ_ｂ）は、κ_ｂに依存する正規化係数である。テンプレートＴＰ_γは、図１４に示すように、ＢＰＭ特徴量ＸＢ（ｔ_ｉ）を構成するＢＰＭ特徴量ＸＢ_ｂ（ｔ_ｉ）にそれぞれ乗算される係数δ_γ，ｂからなる。テンプレートＴＰ_γは、その係数δ_γ，γが最大であり、係数δ_γ，２γ，係数δ_γ，３γ・・・，係数δ_{γ，（「γ」の整数倍）}，・・・がそれぞれ極大となるように設定されている。すなわち、例えば、テンプレートＴＰ_γ＝２は、２フレームおきに拍が存在する楽曲に適合するように構成されている。なお、この例では、ＢＰＭ特徴量ＸＢの尤度Ｐを計算するためにテンプレートＴＰを用いているが、これに代えて確率分布関数（例えば、多項分布、ディリクレ分布、多次元正規分布、多次元ポアソン分布など）を用いても良い。

ＢＰＭ特徴量ＸＢ（ｔ_ｉ）が図１１に示すような値であった場合に、図１４に示すテンプレートＴＰ_γ｛γ＝１，２・・・｝を用いて尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））を計算し、その対数を計算した結果を図１５に例示する。この例では、尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_{ｂ＝４，ｎ}（ｔ_ｉ））が最も大きいので、ＢＰＭ特徴量ＸＢ（ｔ_ｉ）は、テンプレートＴＰ_４に最も適合している。

つぎに、ＣＰＵ１２ａは、ステップＳ１７４にて、尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））の対数と尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））の対数をそれぞれ加算し、その結果を対数観測尤度Ｌ_ｂ，ｎ（ｔ_ｉ）とする。なお、尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））と尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））とを積算した結果の対数を対数観測尤度Ｌ_ｂ，ｎ（ｔ_ｉ）としても同じ結果が得られる。つぎに、ＣＰＵ１２ａは、ステップＳ１７５にて、対数観測尤度計算処理を終了し、その処理を音響信号分析処理（メインルーチン）のステップＳ１８に進める。

つぎに、ＣＰＵ１２ａは、ステップＳ１８にて、図１６に示す拍・テンポ同時推定プログラムをＲＯＭ１２ｂから読み出して実行する。拍・テンポ同時推定プログラムは、音響信号分析プログラムのサブルーチンである。この拍・テンポ同時推定プログラムは、ビタビアルゴリズムを用いて最尤の状態の系列Ｑを計算するプログラムである。ここで、その概略について説明する。ＣＰＵ１２ａは、まず、フレームｔ_０からフレームｔ_ｉまでオンセット特徴量ＸＯ及びＢＰＭ特徴量ＸＢを観測したときにフレームｔ_ｉの状態ｑ_ｂ，ｎの尤度が最大となるような状態の系列を選択した場合の状態ｑ_ｂ，ｎの尤度を尤度Ｃ_ｂ，ｎ（ｔ_ｉ）とするとともに、各状態ｑ_ｂ，ｎに遷移する１つ前のフレームの状態（遷移元の状態）を状態Ｉ_ｂ，ｎ（ｔ_ｉ）として記憶する。つまり、遷移後の状態が状態ｑ_{ｂ＝βｅ，ｎ＝ηｅ}であって、遷移元の状態が状態ｑ_{ｂ＝βｓ，ｎ＝ηｓ}であるとき、状態Ｉ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）は、状態ｑ_{ｂ＝βｓ，ｎ＝ηｓ}である。ＣＰＵ１２ａは、上記のようにして尤度Ｃ及び状態Ｉをフレームｔ_ｌａｓｔまで計算し、その結果を用いて最尤の状態の系列Ｑを選択する。

なお、後述する具体例では、その説明を簡単にするために、分析対象の楽曲の拍周期ｂの値が「３」、「４」及び「５」のうちのいずれかであるものとする。すなわち、対数観測尤度Ｌ_ｂ，ｎ（ｔ_ｉ）が図１７に例示するように計算された場合の拍・テンポ同時推定処理の手順を具体例として説明する。この例では、拍周期ｂの値が「３」、「４」及び「５」以外である状態の観測尤度が十分に小さいものとし、図１７乃至図１９では、拍周期ｂの値が「３」、「４」及び「５」以外である状態の観測尤度の図示を省略する。また、この例では、拍周期ｂの値が「βｓ」であり、且つフレーム数ｎの値が「ηｓ」である状態から、拍周期ｂの値が「βｅ」であり、且つフレーム数ｎの値が「ηｅ」である状態への対数遷移確率Ｔの値は、次のように設定されている。「ηｅ＝０」、「βｅ＝βｓ」、かつ「ηｅ＝βｅ−１」のとき、対数遷移確率Ｔの値は、「−０．２」である。また、「ηｓ＝０」、「βｅ＝βｓ＋１」、かつ「ηｅ＝βｅ−１」のとき、対数遷移確率Ｔの値は、「−０．６」である。また、「ηｓ＝０」、「βｅ＝βｓ−１」、かつ「ηｅ＝βｅ−１」のとき、対数遷移確率Ｔの値は、「−０．６」である。また、「ηｓ＞０」、「βｅ＝βｓ」、かつ「ηｅ＝ηｓ−１」のとき、対数遷移確率Ｔの値は、「０」である。上記以外の対数遷移確率Ｔの値は、「−∞」である。すなわち、フレーム数ｎの値が「０」である状態（ηｓ＝０）から次の状態へ遷移するとき、拍周期ｂの値は「１」だけ増減され得る。このとき、フレーム数ｎの値は、遷移後の拍周期ｂの値より「１」だけ小さい値に設定される。また、フレーム数ｎの値が「０」でない状態（ηｓ≠０）から次の状態へ遷移するとき、拍周期ｂの値は変更されず、フレーム数ｎの値が「１」だけ減少する。

以下、拍・テンポ同時推定処理について具体的に説明する。ＣＰＵ１２ａは、ステップＳ１８１にて拍・テンポ同時推定処理を開始する。つぎに、ユーザは、ステップＳ１８２にて、入力操作子１１を用いて、図１８に示すような、各状態ｑ_ｂ，ｎに対応した、尤度Ｃの初期条件ＣＳ_ｂ，ｎを入力する。なお、初期条件ＣＳ_ｂ，ｎがＲＯＭ１２ｂに記憶されていて、ＣＰＵ１２ａがＲＯＭ１２ｂから初期条件ＣＳ_ｂ，ｎを読み込むようにしてもよい。

つぎに、ＣＰＵ１２ａは、ステップＳ１８３にて、尤度Ｃ_ｂ，ｎ（ｔ_ｉ）及び状態Ｉ_ｂ，ｎ（ｔ_ｉ）を計算する。フレームｔ_０において拍周期ｂの値が「βｅ」であって、フレーム数ｎの値が「ηｅ」である状態ｑ_{ｂ＝βｅ，ｎ＝ηｅ}の尤度Ｃ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_０）は、初期条件ＣＳ_{ｂ＝βｅ，ｎ＝ηｅ}と対数観測尤度Ｌ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_０）とを加算することにより計算される。

また、状態ｑ_{ｂ＝βｓ，ｎ＝ηｓ}から状態ｑ_{ｂ＝βｅ，ｎ＝ηｅ}に遷移したとき、尤度Ｃ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）｛ｉ＞０｝は次のように計算される。状態ｑ_{ｂ＝βｓ，ｎ＝ηｓ}のフレーム数ｎが「０」でないとき（すなわち、ηｓ≠０）、尤度Ｃ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）は、尤度Ｃ_{ｂ＝βｅ，ｎ＝ηｅ＋１}（ｔ_ｉ―１）と対数観測尤度Ｌ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）と対数遷移確率Ｔを加算して計算される。ただし、本実施形態では、遷移元の状態のフレーム数ｎが「０」でないときの対数遷移確率Ｔは「０」であるので、尤度Ｃ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）は、実質的には、尤度Ｃ_{ｂ＝βｅ，ｎ＝ηｅ＋１}（ｔ_ｉ―１）と対数観測尤度Ｌ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）とを加算することにより計算される（Ｃ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）＝Ｃ_{ｂ＝βｅ，ｎ＝ηｅ＋１}（ｔ_ｉ―１）＋Ｌ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ））。また、この場合、状態Ｉ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）は、状態ｑ_{βｅ，ηｅ＋１}である。例えば、尤度Ｃが図１８に示すように計算された例では、尤度Ｃ_４，１（ｔ_２）の値は「２」であり、対数観測尤度Ｌ_４，０（ｔ_３）の値は「１」であるので、尤度Ｃ_４，０（ｔ_３）の値は「３」である。また、図１９に示すように、状態Ｉ_４，０（ｔ_３）は、状態ｑ_４，１である。

また、状態ｑ_{ｂ＝βｓ，ｎ＝ηｓ}のフレーム数ｎが「０」のとき（ηｓ＝０）の尤度Ｃ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）は次のように計算される。この場合、状態の遷移に伴って拍周期ｂの値が増減され得る。そこで、まず、尤度Ｃ_{βｅ−１，０}（ｔ_ｉ−１）、尤度Ｃ_βｅ，０（ｔ_ｉ−１）、及び尤度Ｃ_{βｅ＋１，０}（ｔ_ｉ−１）に対数遷移確率Ｔをそれぞれ加算し、そのうちの最大値に対数観測尤度Ｌ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）を加算した結果が尤度Ｃ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）である。また、状態Ｉ_{ｂ＝βｅ，ｎ＝ηｅ}（ｔ_ｉ）は、状態ｑ_{βｅ−１，０}、状態ｑ_βｅ，０、及び状態ｑ_{βｅ＋１，０}のうち、それらの尤度Ｃ_{βｅ−１，０}（ｔ_ｉ−１）、尤度Ｃ_βｅ，０（ｔ_ｉ−１）、及び尤度Ｃ_{βｅ＋１，０}（ｔ_ｉ−１）に対数遷移確率Ｔをそれぞれ加算した値が最大となる状態ｑである。なお、厳密には、尤度Ｃ_ｂ，ｎ（ｔ_ｉ）は正規化される必要があるが、正規化されていなくても、拍点及びテンポの推移の推定に関しては、数理上同一の結果が得られる。

例えば、尤度Ｃ_４，３（ｔ_４）は、次のように計算される。遷移元の状態が状態ｑ_３，０である場合、尤度Ｃ_３，０（ｔ_３）の値は「０．４」であり、対数遷移確率Ｔは「−０．６」であるので、尤度Ｃ_３，０（ｔ_３）と対数遷移確率Ｔとを加算した値は、「−０．２」である。また、遷移元の状態が状態ｑ_４，０である場合、遷移元の尤度Ｃ_４，０（ｔ_３）の値は「３」であり、対数遷移確率Ｔは「−０．２」であるので、尤度Ｃ_４，０（ｔ_３）と対数遷移確率Ｔとを加算した値は、「２．８」である。また、遷移元の状態が状態ｑ_５，０である場合、遷移元の尤度Ｃ_５，０（ｔ_３）の値は「１」であり、対数遷移確率Ｔは「−０．６」であるので、尤度Ｃ_５，０（ｔ_３）と対数遷移確率Ｔとを加算した値は、「０．４」である。したがって、尤度Ｃ_４，０（ｔ_３）に対数遷移確率Ｔを加算した値が最も大きい。また、対数観測尤度Ｌ_４，３（ｔ_４）の値は、「０」である。よって、尤度Ｃ_４，３（ｔ_４）の値は「２．８」（＝２．８＋０）であり、状態Ｉ_４，３（ｔ_４）は、状態ｑ_４，０である。

上記のようにして、全てのフレームｔ_ｉについて、全ての状態ｑ_ｂ，ｎの尤度Ｃ_ｂ，ｎ（ｔ_ｉ）及び状態Ｉ_ｂ，ｎ（ｔ_ｉ）を計算し終えると、ＣＰＵ１２ａはステップＳ１８４にて、最尤の状態の系列Ｑ（＝｛ｑ_ｍａｘ（ｔ_０），ｑ_ｍａｘ（ｔ_１）・・・，ｑ_ｍａｘ（ｔ_ｌａｓｔ）｝）を次のようにして決定する。まず、ＣＰＵ１２ａは、フレームｔ_ｌａｓｔにおける尤度Ｃ_ｂ，ｎ（ｔ_ｌａｓｔ）が最大である状態ｑ_ｂ，ｎを、状態ｑ_ｍａｘ（ｔ_ｌａｓｔ）とする。ここで、状態ｑ_ｍａｘ（ｔ_ｌａｓｔ）の拍周期ｂの値を「βｍ」と表記し、フレーム数ｎの値を「ηｍ」と表記する。このとき、状態Ｉ_{βｍ，ηｍ}（ｔ_ｌａｓｔ）がフレームｔ_ｌａｓｔの１つ前のフレームｔ_{ｌａｓｔ−１}の状態ｑ_ｍａｘ（ｔ_{ｌａｓｔ−１}）である。フレームｔ_{ｌａｓｔ−２}、フレームｔ_{ｌａｓｔ−３}、・・・の状態ｑ_ｍａｘ（ｔ_{ｌａｓｔ−２}）、状態ｑ_ｍａｘ（ｔ_{ｌａｓｔ−３}）・・・も状態ｑ_ｍａｘ（ｔ_{ｌａｓｔ−１}）と同様に決定される。すなわち、フレームｔ_ｉ＋１の状態ｑ_ｍａｘ（ｔ_ｉ＋１）の拍周期ｂの値を「βｍ」と表記し、フレーム数ｎの値を「ηｍ」と表記したときの状態Ｉ_{βｍ，ηｍ}（ｔ_ｉ＋１）がフレームｔ_ｉ＋１の１つ前のフレームｔ_ｉの状態ｑ_ｍａｘ（ｔ_ｉ）である。上記のようにして、ＣＰＵ１２ａは、フレームｔ_{ｌａｓｔ−１}からフレームｔ_０へ向かって順に状態ｑ_ｍａｘを決定して、最尤の状態の系列Ｑを決定する。

例えば、図１８及び図１９に示す例では、フレームｔ_{ｌａｓｔ＝９}においては、状態ｑ_４，２の尤度Ｃ_４，２（ｔ_{ｌａｓｔ＝９}）が最大である。したがって、状態ｑ_ｍａｘ（ｔ_{ｌａｓｔ＝９}）は、状態ｑ_４，２である。図１９によれば、状態Ｉ_４，２（ｔ_９）は状態ｑ_４，３であるから、状態ｑ_ｍａｘ（ｔ_８）は状態ｑ_４，３である。また、状態Ｉ_４，３（ｔ_８）は状態ｑ_４，０であるから、状態ｑ_ｍａｘ（ｔ_７）は状態ｑ_４，０である。状態ｑ_ｍａｘ（ｔ_６）乃至状態ｑ_ｍａｘ（ｔ_０）も状態ｑ_ｍａｘ（ｔ_８）及び状態ｑ_ｍａｘ（ｔ_７）と同様に決定する。このようにして図１８に矢印で示す最尤の状態の系列Ｑが決定される。この例では、拍の周期ｂの値はいずれのフレームｔ_ｉにおいても「４」であると推定される。また、系列Ｑのうち、フレーム数ｎの値が「０」である状態ｑ_ｍａｘ（ｔ_１），ｑ_ｍａｘ（ｔ_５），ｑ_ｍａｘ（ｔ_８）に対応するフレームｔ_１，ｔ_５，ｔ_８に拍が存在すると推定される。

つぎに、ＣＰＵ１２ａは、ステップＳ１８５にて、拍・テンポ同時推定処理を終了し、その処理を音響信号分析処理（メインルーチン）のステップＳ１９に進める。

ＣＰＵ１２ａは、ステップＳ１９にて、フレームｔ_ｉごとに「ＢＰＭらしさ」、「観測に基づく確率」、「拍らしさ」、「拍が存在する確率」及び「拍が存在しない確率」を計算（図２０に示す演算式を参照）する。「ＢＰＭらしさ」は、フレームｔ_ｉにおけるテンポの値が拍周期ｂに対応した値である確率を意味し、尤度Ｃ_ｂ，ｎ（ｔ_ｉ）を正規化するとともにフレーム数ｎについて周辺化することにより計算される。具体的には、拍周期ｂの値が「β」である場合の「ＢＰＭらしさ」は、フレームｔ_ｉにおける全ての状態の尤度Ｃの合計に対する、拍周期ｂの値が「β」である状態の尤度Ｃの合計の割合である。また、「観測に基づく確率」は、観測値（すなわちオンセット特徴量ＸＯ）に基づいて計算された拍がフレームｔ_ｉに存在する確率を意味する。具体的には、所定の基準値ＸＯ_ｂａｓｅに対するオンセット特徴量ＸＯ（ｔ_ｉ）の割合である。また、「拍らしさ」は、すべてのフレーム数ｎの値についてのオンセット特徴量ＸＯ（ｔ_ｉ）の尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））を合算した値に対する尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_ｂ，０（ｔ_ｉ））の割合である。また、「拍が存在する確率」及び「拍が存在しない確率」は、いずれも尤度Ｃ_ｂ，ｎ（ｔ_ｉ）を拍周期ｂについて周辺化することにより計算される。具体的には、「拍が存在する確率」は、フレームｔ_ｉにおける全ての状態の尤度Ｃの合計に対する、フレーム数ｎの値が「０」である状態の尤度Ｃの合計の割合である。また、「拍が存在しない確率」は、フレームｔ_ｉにおける全ての状態の尤度Ｃの合計に対する、フレーム数ｎの値が「０」でない状態の尤度Ｃの合計の割合である。

ＣＰＵ１２ａは、「ＢＰＭらしさ」、「観測に基づく確率」、「拍らしさ」、「拍が存在する確率」及び「拍が存在しない確率」を用いて、図２０に示す拍・テンポ情報リストを表示器１３に表示する。同リスト中の「テンポの推定値（ＢＰＭ）」の欄には、前記計算した「ＢＰＭらしさ」のうち最も確率の高い拍周期ｂに対応するテンポの値（ＢＰＭ）が表示される。また、前記決定した状態ｑ_ｍａｘ（ｔ_ｉ）のうちフレーム数ｎの値が「０」であるフレームの「拍の存在」の欄には「○」が表示され、その他のフレームの「拍の存在」の欄には「×」が表示される。また、ＣＰＵ１２ａは、テンポの推定値（ＢＰＭ）を用いて、図２１に示すようなテンポの推移を表わすグラフを表示器１３に表示する。図２１の例では、テンポの推移を棒グラフで表わしている。図１８及び図１９を用いて説明した例では、テンポの値が一定であるので図２１に示すような各フレームのテンポを表わすバーの高さが一定であるが、テンポが頻繁に変化する楽曲では、図２２に示すように、テンポの値に応じてバーの高さが異なる。これにより、ユーザは、テンポの推移を視覚的に認識することができる。また、ＣＰＵ１２ａは、前記計算した「拍が存在する確率」を用いて、図２３に示すような拍点を表わすグラフを表示器１３に表示する。

また、音響信号分析処理のステップＳ１３にて既存データを検索した結果、既存データが存在する場合には、ＣＰＵ１２ａは、ステップＳ１５にてＲＡＭ１２ｃに読み込んだ前回の分析結果に関する各種データを用いて、拍・テンポ情報リスト、テンポの推移を表わすグラフ、及び拍点を表わすグラフを表示器１３に表示する。

次に、ＣＰＵ１２ａは、ステップＳ２０にて、音響信号分析処理を終了するか否かを表すメッセージを表示器１３に表示して、ユーザからの指示を待つ。ユーザは入力操作子１１を用いて音響信号分析処理を終了するか、後述の拍・テンポ情報修正処理を実行するかのいずれかを指示する。例えば、マウスを用いて図示しないアイコンをクリックする。ユーザから音響信号分析処理を終了するよう指示された場合には、ＣＰＵ１２ａは「Ｙｅｓ」と判定してステップＳ２１にて尤度Ｃ、状態Ｉ、拍・テンポ情報リストなどの分析結果に関する各種データを楽曲のタイトルと関連付けて記憶装置１４に記憶して、ステップＳ２２にて音響信号分析処理を終了する。

一方、ステップＳ２０にて、音響信号分析処理を継続するように指示された場合には、ＣＰＵ１２ａは「Ｎｏ」と判定して、ステップＳ２３にて、テンポ情報修正処理を実行する。まず、ＣＰＵ１２ａは、ユーザが修正情報の入力を終了するまで待機する。ユーザは、入力操作子１１を用いて「ＢＰＭらしさ」、「拍が存在する確率」などの修正値を入力する。例えば、マウスを用いて修正するフレームを選択し、テンキーを用いて修正値を入力する。修正された項目の右側に配置された「Ｆ」の表示形態（例えば色）が変更され、その値が修正されたことが明示される。ユーザは、複数の項目について修正値を入力可能である。ユーザは修正値の入力を完了すると、入力操作子１１を用いて修正情報の入力を完了したことを指示する。例えば、マウスを用いて図示しない修正完了を表わすアイコンをクリックする。ＣＰＵ１２ａは、前記入力された修正値に応じて尤度Ｐ（ＸＯ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））及び尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_ｂ，ｎ（ｔ_ｉ））のうちのいずれか一方又は両方を更新する。例えば、フレームｔ_ｉにおける「拍が存在する確率」が高くなるように修正された場合であって、修正された値に関するフレーム数ｎの値が「ηｅ」であるときには、尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_{ｂ，ｎ≠ηｅ}（ｔ_ｉ））を十分に小さい値に設定する。これにより、フレームｔ_ｉでは、フレーム数ｎの値が「ηｅ」である確率が相対的に最も高くなる。また、例えば、フレームｔ_ｉにおける「ＢＰＭらしさ」のうち、拍周期ｂの値が「βｅ」である確率が高くなるように修正された場合には、拍周期ｂの値が「βｅ」でない状態の尤度Ｐ（ＸＢ（ｔ_ｉ）｜Ｚ_{ｂ≠βｅ，ｎ}（ｔ_ｉ））を十分に小さい値に設定する。これにより、フレームｔ_ｉでは、拍周期ｂの値が「βｅ」である確率が相対的に最も高くなる。そして、ＣＰＵ１２ａは、拍・テンポ情報修正処理を終了して、その処理をステップＳ１８に進め、修正された対数観測尤度Ｌを用いて、拍・テンポ同時推定処理を再度実行する。

上記のように構成した音響信号分析装置１０によれば、拍点に関するオンセット特徴量ＸＯ及びテンポに関するＢＰＭ特徴量ＸＢを用いて計算された対数観測尤度Ｌの系列が最も尤もらしい確率モデルが選択され、楽曲における拍点及びテンポの推移が同時に推定される。したがって、上記従来技術とは異なり、拍点及びテンポのうちの一方の推定精度が低いために他方の推定精度も低くなるという事態が生じない。よって、従来技術に比べて楽曲における拍点及びテンポの推移の推定精度を向上させることができる。

また、本実施形態においては、フレーム数ｎの値が「０」である状態から、拍周期ｂの値が同じ状態又は拍周期ｂの値が「１」だけ異なる状態へのみ遷移可能に各状態間の遷移確率（対数遷移確率）が設定されている。これにより、テンポがフレーム間で急激に変化するような誤推定が防止される。したがって、楽曲として自然な拍点及びテンポの推移の推定結果を得ることができる。なお、テンポが急激に変化する楽曲に対しては、次の拍までのフレーム数ｎの値が「０」である状態から次の状態に遷移するとき、拍周期ｂの値が大きく異なる状態への遷移も可能なように各状態間の遷移確率（対数遷移確率）を設定すればよい。

また、拍・テンポ同時推定処理では、ビタビアルゴリズムを用いたので、他のアルゴリズム（例えば、「サンプリング法」、「前向き後向きアルゴリズム」など）を用いる場合に比べて計算量を削減できる。

また、ユーザにより入力された修正情報に基づいて対数観測尤度Ｌが修正され、修正された対数観測尤度Ｌに基づいて楽曲における拍点及びテンポの推移が再推定される。これにより、修正されたフレームの前後にそれぞれ位置する１つ又は複数のフレームの最尤の状態ｑ_ｍａｘが再計算（再選択）される。したがって、修正されたフレーム及びその前後に位置する１つ又は複数のフレームに亘り、拍の間隔及びテンポが滑らかに変化するような推定結果が得られる。

上記のようにして推定された楽曲における拍点及びテンポの推移に関する情報は、例えば楽曲データの検索、伴奏を表わす伴奏データの検索などに利用される。また、分析対象とした楽曲に対する伴奏パートの自動生成、ハーモニーの自動付加などにも利用される。

さらに、本発明の実施にあたっては、上記実施形態に限定されるものではなく、本発明の目的を逸脱しない限りにおいて種々の変更が可能である。

例えば、上記実施形態では、観測値としてのオンセット特徴量ＸＯ及びＢＰＭ特徴量ＸＢが同時に観測される確率を表わす観測尤度の系列が最も尤もらしい確率モデルが選択される。しかし、確率モデルの選択基準は、上記実施形態に限られない。例えば、事後分布が最大となるような確率モデルを選択してもよい。

また、例えば、上記実施形態では、説明を簡単にするために、各フレームの長さを１２５ｍｓとしたが、より短く（例えば、５ｍｓ）してもよい。これによれば、拍点及びテンポの推定に関する分解能を向上させることができる。例えば、テンポを１ＢＰＭ刻みで推定できる。また、上記実施形態では、各フレームの長さを共通にしているが、各フレームの長さが異なっていてもよい。この場合であっても、オンセット特徴量ＸＯは、上記実施形態と同様にして計算できる。また、この場合、ＢＰＭ特徴量ＸＢの計算においては、コムフィルタの遅延量をフレームの長さに応じて変更すればよい。また、尤度Ｃの計算においては、各フレームの長さの最大公約数Ｆ（つまり、各フレームを構成するサンプル数の
最大公約数）を計算する。そして、フレームｔ_ｉ（＝τ）の長さがＬ（τ）×Ｆと表わされたとき、状態ｑ_{ｂ，ｎ（ｎ≠０）}から状態ｑ_{ｂ，ｎ−Ｌ（τ）}へ遷移する確率を１００％とすればよい。

また、上記実施形態では、楽曲全体を分析対象としているが、楽曲の一部（例えば数小節）のみを分析対象としてもよい。この場合、入力した楽曲データのうち、分析対象とする部分を選択可能に構成するとよい。また、楽曲のうちの単一のパート（例えばリズムセクション）のみを分析対象としてもよい。

また、例えば、テンポの推定において、優先的に推定するテンポの範囲を指定可能に構成してもよい。具体的には、音響信号分析処理のステップＳ１２において、「Ｐｒｅｓｔｏ」、「Ｍｏｄｅｒａｔｏ」などのテンポを表わす用語を表示して、優先的に推定するテンポの範囲を選択可能に構成してもよい。例えば、「Ｐｒｅｓｔｏ」が選択された場合、ＢＰＭ＝１６０〜１９０の範囲以外の対数観測尤度Ｌを十分に小さく設定する。これにより、ＢＰＭ＝１６０〜１９０の範囲のテンポが優先的に推定される。これによれば、楽曲の大凡のテンポが既知である場合、テンポの推定精度を向上させることができる。

また、拍・テンポ情報修正処理（ステップＳ２３）では、ユーザは入力操作子１１を用いて修正内容を入力するように構成されている。これに代えて、又は加えて、外部インターフェース回路１５を介して接続された電子鍵盤楽器、電子打楽器などの操作子を用いて修正内容を入力可能に構成してもよい。例えば、ユーザが電子鍵盤楽器の鍵盤を数回打鍵すると、ＣＰＵ１２ａがその打鍵のタイミングからテンポを計算して、前記「ＢＰＭらしさ」の修正値として用いるように構成してもよい。

また、上記実施形態では、拍点及びテンポに関する修正値を何度でも入力可能に構成されている。しかし、例えば、「拍が存在する確率」の平均値が基準値（例えば８０％）に達した時点以降においては、拍点及びテンポに関する修正値を入力不可能としてもよい。

また、例えば、拍・テンポ情報修正処理（ステップＳ２３）において、ユーザによって指定されたフレームの拍・テンポ情報を入力された値に修正するとともに、そのフレームに近接するフレームの拍・テンポ情報を前記入力された値に応じて自動的に修正してもよい。例えば、連続する複数のフレームのテンポの推定値が同じ値であって、そのうちの１つのフレームのテンポの値が修正されたとき、前記複数のフレームのテンポの値を前記１つのフレームの修正値と同じ値に自動的に修正してもよい。

また、上記実施形態では、ステップＳ２３にて、ユーザが入力操作子１１を用いて修正値の入力を完了したことを指示すると、拍点及びテンポの同時推定が再び実行される。しかし、これに代えて、ユーザが少なくとも１つの修正値を入力した後、他の修正値が入力されないまま所定の時間（例えば１０秒）が経過したとき、自動的に拍点及びテンポの同時推定が再び実行されてもよい。

また、拍・テンポ情報リスト（図２０）の表示形態は、上記実施形態に限られない。例えば、上記実施形態では、「ＢＰＭらしさ」、「拍点らしさ」などを、確率（％）で表示しているが、これらを記号、文字列などを用いて表現してもよい。また、上記実施形態では、前記決定した状態ｑ_ｍａｘ（ｔ_ｉ）のうちフレーム数ｎの値が「０」であるフレームｔ_ｉの「拍の存在」の欄には「○」が表示され、その他のフレームの「拍の存在」の欄には「×」が表示されるが、これに代えて、例えば、基準値（例えば８０％）以上であるとき「拍の存在」の欄に「○」が表示され、「拍点が存在する確率」が基準値未満であるとき、「拍の存在」の欄に「×」が表示されてもよい。また、この場合、複数の基準値を設けてもよい。例えば、第１基準値（＝８０％）と第２基準値（＝６０％）を設け、「拍点が存在する確率」が第１基準値以上であるとき、「拍の存在」の欄に「○」が表示され、第２基準値以上かつ第１基準値未満であるとき、「拍の存在」の欄に「△」が表示され、「拍点が存在する確率」が第２基準値未満であるとき、「拍の存在」の欄に「×」が表示されてもよい。また、テンポの推定値の欄には、「Ｐｒｅｓｔｏ」、「Ｍｏｄｅｒａｔｏ」などのテンポを表わす用語が表示されてもよい。

１０・・・音響信号分析装置、１１・・・入力操作子、ＸＯ・・・オンセット特徴量、ＸＢ・・・ＢＰＭ特徴量、ｂ・・・拍周期、ｎ・・・フレーム数、ＦＢＢ・・・フィルタバンク、ＴＰ・・・テンプレート

Claims (7)

1. 楽曲を表わす音響信号を入力する音響信号入力手段と、
   前記楽曲の各区間における拍の存在に関する特徴を表わす第１特徴量及びテンポに関する特徴を表わす第２特徴量を計算する特徴量計算手段と、
   前記楽曲の各区間における拍の存在に関する物理量及びテンポに関する物理量の組み合わせにより分類された状態の系列として記述された複数の確率モデルのうち、前記第１特徴量及び前記第２特徴量が前記楽曲の各区間において同時に観測される確率を表わす観測尤度の系列が所定の基準を満たす確率モデルを選択することにより、前記楽曲における拍点及びテンポの推移を同時に推定する推定手段と、
   を備えたことを特徴とする音響信号分析装置。
2. 請求項１に記載の音響信号分析装置において、
   前記推定手段は、
   前記拍の存在に関する物理量に応じて規定された確率分布関数の確率変数として前記第１特徴量を代入して計算された確率を前記第１特徴量が観測される確率として出力する第１確率出力手段と、
   前記テンポに関する物理量に応じて規定された確率分布関数の確率変数として前記第２特徴量を代入して計算された確率を前記第２特徴量が観測される確率として出力する第２確率出力手段と、
   を備えたことを特徴とする音響信号分析装置。
3. 請求項２に記載の音響信号分析装置において、
   前記第１確率出力手段は、前記拍の存在に関する物理量に応じて規定された正規分布の確率変数として前記第１特徴量を代入して計算された確率を前記第１特徴量が観測される確率として出力することを特徴とする音響信号分析装置。
4. 請求項２に記載の音響信号分析装置において、
   前記第２確率出力手段は、前記テンポに関する物理量に応じてそれぞれ設定された複数のテンプレートに対する前記第２特徴量の適合度を前記第２特徴量が観測される確率として出力することを特徴とする音響信号分析装置。
5. 請求項１乃至４のうちのいずれか１つに記載の音響信号分析装置において、
   前記楽曲の各区間は、前記入力した音響信号を所定の時間間隔をおいて区切ることにより形成された各フレームにそれぞれ相当し、
   前記特徴量計算手段は、
   前記フレームごとに振幅スペクトルを計算し、前記振幅スペクトルに周波数帯域が互いに異なる複数の窓関数を掛けて前記周波数帯域ごとの振幅スペクトルを生成し、前記フレーム間における前記周波数帯域ごとの振幅スペクトルの変化に基づいて計算した値を前記第１特徴量として出力する第１特徴量計算手段と、
   各フレームに対応する値を入力するごとに１つの値を出力するフィルタであって、出力した値を所定の時間だけ保持する保持手段を備え、前記入力した値と前記所定の時間だけ保持された値とを所定の比率で加算して出力するフィルタを備え、
   前記第１特徴量の系列を前記フィルタに入力して得られたデータ列の時系列を逆にしたデータ列を前記フィルタに再び入力して得られたデータ列を前記第２特徴量の系列として出力する第２特徴量計算手段と、を備えたことを特徴とする音響信号分析装置。
6. 請求項１乃至５のうちのいずれか１つに記載の音響信号分析装置において、
   前記楽曲における拍点及びテンポの推移のいずれか一方又は両方の修正内容を表わす修正情報を入力する修正情報入力手段と、
   前記入力した修正情報に応じて前記観測尤度を修正する観測尤度修正手段と、
   前記複数の確率モデルのうち、前記修正された観測尤度の系列が前記所定の基準を満たす確率モデルを、前記推定手段を用いて選択することにより、前記楽曲における拍点及びテンポの推移を同時に再推定する再推定手段と、を備えたことを特徴とする音響信号分析装置。
7. コンピュータに、
   楽曲を表わす音響信号を入力する音響信号入力ステップと、
   前記楽曲の各区間における拍の存在に関する特徴を表わす第１特徴量及びテンポに関する特徴を表わす第２特徴量を計算する特徴量計算ステップと、
   前記楽曲の各区間における拍の存在に関する物理量及びテンポに関する物理量の組み合わせにより分類された状態の系列として記述された複数の確率モデルのうち、前記第１特徴量及び前記第２特徴量が前記楽曲の各区間において同時に観測される確率を表わす観測尤度の系列が所定の基準を満たす確率モデルを選択することにより、前記楽曲における拍点及びテンポの推移を同時に推定する推定ステップと、
   を実行させることを特徴とする音響信号分析プログラム。

Images