(57)【要約】      （修正有）【課題】拍点、テンポ、小節線の位置及びコード進行の推定精度を向上させた音響信号分析装置及び音響信号分析プログラムを提供する。【解決手段】音響信号分析装置は、分析対象としての楽曲の演奏音を表わす音響信号を取り込む。次に楽曲の各区間において発音される和音に関する特徴を表わす第１特徴量、拍の存在に関する特徴を表わす第２特徴量及びテンポに関する特徴を表わす第３特徴量を計算し、前記第１乃至第３特徴量の観測尤度を計算する。次に、前記各区間が属する小節内における前記各区間の位置に関連して各状態間の遷移確率が設定されている第１モデルの事後分布と、楽曲の各区間における拍の存在に関する物理量、テンポに関する物理量及び各区間が属する小節内における各区間の位置に関する物理量の組み合わせの系列として第２モデルの事後分布とを、観測尤度及び所定の事前分布を用いて同時（一体的）に計算する。【選択図】図２Ｂ
2015200803.tif 000002
【特許請求の範囲】
【請求項１】
分析対象としての楽曲の演奏音を表わす音響信号を取り込む音響信号取得手段と、
前記取り込んだ音響信号に基づいて、前記楽曲の各区間において発音される和音に関する特徴を表わす第１特徴量、拍の存在に関する特徴を表わす第２特徴量及びテンポに関する特徴を表わす第３特徴量を計算する特徴量計算手段と、
前記第１乃至第３特徴量の観測尤度を計算する観測尤度計算手段と、
前記楽曲の各区間において発音される和音に関する物理量の系列として記述された第１モデルであって前記各区間が属する小節内における前記各区間の位置に関連して各状態間の遷移確率が設定されている第１モデルの事後分布と、前記楽曲の各区間における拍の存在に関する物理量、テンポに関する物理量及び前記各区間が属する小節内における前記各区間の位置に関する物理量の組み合わせの系列として記述された第２モデルの事後分布とを、前記計算した観測尤度及び所定の事前分布を用いて同時に計算する事後分布計算手段と、
を備えた音響信号分析装置。
【請求項２】
請求項１に記載の音響信号分析装置において、
前記事後分布計算手段は、変分ベイズ推定法を用いて、前記第１モデル及び前記第２モデルの事後分布をそれぞれ近似的に計算する、音響信号分析装置。
【請求項３】
請求項２に記載の音響信号分析装置において、
前記第１モデル及び前記第２モデルの各状態の遷移確率の事前分布は、ディリクレ分布に従う、音響信号分析装置。
【請求項４】
音響信号分析装置が備えるコンピュータに、
分析対象としての楽曲の演奏音を表わす音響信号を取り込む音響信号取得ステップと、
前記取り込んだ音響信号に基づいて、前記楽曲の各区間において発音される和音に関する特徴を表わす第１特徴量、拍の存在に関する特徴を表わす第２特徴量及びテンポに関する特徴を表わす第３特徴量を計算する特徴量計算ステップと、
前記第１乃至第３特徴量の観測尤度を計算する観測尤度計算ステップと、
前記楽曲の各区間において発音される和音に関する物理量の系列として記述された第１モデルであって前記各区間が属する小節内における前記各区間の位置に関連して各状態間の遷移確率が設定されている第１モデルの事後分布と、前記楽曲の各区間における拍の存在に関する物理量、テンポに関する物理量及び前記各区間が属する小節内における前記各区間の位置に関する物理量の組み合わせの系列として記述された第２モデルの事後分布とを、前記計算した観測尤度及び所定の事前分布を用いて計算する事後分布計算ステップと、
を実行させる、コンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、楽曲を表わす音響信号を分析して、楽曲における拍点（拍のタイミング）、テンポ、小節線の位置、及び楽曲の各区間で発音されるコード（和音）を推定する音響信号分析装置及び音響信号分析プログラムに関する。
【背景技術】
【０００２】
従来から、例えば、下記非特許文献１に記載されているように、楽曲における拍点、及び楽曲の各区間で発音されるコード（和音）を推定する音響信号分析装置は知られている。この音響信号分析装置は、まず、音響信号を分析して、楽曲の拍点を推定している。そして、前記推定した拍点においてコード変化が生起し、かつ小節の先頭でコード変化が生起するという仮定の下で、楽曲の各区間のコード及び小節線の位置を推定している。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】Ｍ．Ｇｏｔｏ  ｅｔ  ａｌ．  、“ＳＯＮＧＬＥ：  Ａ  ＷＥＢ  ＳＥＲＶＩＣＥ  ＦＯＲ  ＡＣＴＩＶＥ  ＭＵＳＩＣ  ＬＩＳＴＥＮＩＮＧ  ＩＭＰＲＯＶＥＤ  ＢＹ  ＵＳＥＲ  ＣＯＮＴＲＩＢＵＴＩＯＮＳ”、ＩＳＭＩＲ、２０１１、ｐ．３１１−３１６
【発明の概要】
【０００４】
上記非特許文献１に記載の音響信号分析装置においては、拍点の推定精度が低い場合には、コードの推定精度及び小節線の位置の推定精度も低くなる。
【０００５】
本発明は上記問題に対処するためになされたもので、その目的は、拍点、テンポ、小節線の位置及びコード進行の推定精度を向上させた音響信号分析装置及び音響信号分析プログラムを提供することにある。なお、下記本発明の各構成要件の記載においては、本発明の理解を容易にするために、実施形態の対応箇所の符号を括弧内に記載しているが、本発明の各構成要件は、実施形態の符号によって示された対応箇所の構成に限定解釈されるべきものではない。
【０００６】
上記目的を達成するために、本発明の特徴は、分析対象としての楽曲の演奏音を表わす音響信号を取り込む音響信号取得手段（Ｓ１１）と、前記取り込んだ音響信号に基づいて、前記楽曲の各区間において発音される和音に関する特徴を表わす第１特徴量（ｘ（ｃ））、拍の存在に関する特徴を表わす第２特徴量（ｘ（ｏ））及びテンポに関する特徴を表わす第３特徴量（ｘ（ｂ））を計算する特徴量計算手段（Ｓ１４）と、前記第１乃至第３特徴量の観測尤度（Ｌ（ｃ），Ｌ（ｏ），Ｌ（ｂ），Ｌ（ｏ，ｂ））を計算する観測尤度計算手段（Ｓ１５）と、前記楽曲の各区間において発音される和音に関する物理量の系列として記述された第１モデル（Ｚ（ｃ））であって前記各区間が属する小節内における前記各区間の位置に関連して各状態間の遷移確率（τｉ，ｊ（ｃ））が設定されている第１モデルの事後分布と、前記楽曲の各区間における拍の存在に関する物理量（ｎ）、テンポに関する物理量（ｂ）及び前記各区間が属する小節内における前記各区間の位置に関する物理量（ｓ）の組み合わせの系列として記述された第２モデル（Ｚ（ｏ））の事後分布とを、前記計算した観測尤度及び所定の事前分布を用いて同時（一体的）に計算する事後分布計算手段（Ｓ１６〜Ｓ２６）と、を備えた音響信号分析装置（１０）としたことにある。
【０００７】
この場合、前記事後分布計算手段は、変分ベイズ推定法を用いて、前記第１モデル及び前記第２モデルの事後分布をそれぞれ近似的に計算するとよい。
【０００８】
また、この場合、前記第１モデル及び前記第２モデルの各状態の遷移確率の事前分布は、ディリクレ分布に従うとよい。
【０００９】
一般に、コードの変化は、拍点において生起する可能性が高い。また、遷移元のコードと遷移先のコードは、遷移先（又は遷移元）の区間が属する小節内における前記区間の位置（言い換えれば、前記区間からみた直前の小節線から数えた拍数。）に依存する。以下、「遷移先（又は遷移元）の区間が属する小節内における前記区間の位置」を「拍子位置」と呼ぶ。そこで、本発明に係る音響信号分析装置においては、各区間のコードの系列として記述された第１モデルの各状態間の遷移確率が、拍子位置に関連して設定されている。そして、第１モデルの事後分布と、前記楽曲の各区間における拍の存在に関する物理量、テンポに関する物理量及び前記各区間が属する小節内における前記各区間の位置（つまり拍子位置）に関する物理量の組み合わせの系列として記述された第２モデルの事後分布とが、同時（一体的）に計算される。第１モデルと第２モデルとは、拍子位置に関して互いに依存しており、この両モデルの事後分布を同時（一体的）に計算することにより、楽曲における拍点、テンポ、拍子位置、及びコードに関する情報が精度良く得られる。言い換えれば、一方のモデルの推定結果が他方のモデルの推定結果にフィードバックされて推定結果が更新される。これにより、両モデルの推定精度が高められる。したがって、本発明に係る音響信号分析装置によれば、拍点、テンポ、小節線の位置及びコード進行の推定精度を従来よりも向上させることができる。
【００１０】
また、本発明は、音響信号分析装置が備えるコンピュータに適用されるコンピュータプログラムとしても実施可能である。
【図面の簡単な説明】
【００１１】
【図１】本発明の一実施形態に係る音響信号分析装置の構成を表わすブロック図である。
【図２Ａ】音響信号分析処理の前半部分を表わすフローチャートである。
【図２Ｂ】音響信号分析処理の後半部分を表わすフローチャートである。
【図３】分析対象の楽曲を表わす音響信号の波形を表わすグラフである。
【図４】コード特徴量の概念図である。
【図５】コムフィルタのブロック図である。
【図６】ＢＰＭ特徴量の計算結果を示すグラフである。
【図７】テンプレートの構成を示す表である。
【図８】推定結果の表示例である。
【発明を実施するための形態】
【００１２】
本発明の一実施形態に係る音響信号分析装置１０について説明する。音響信号分析装置１０は、以下説明するように、楽曲を表わす音響信号を取り込んで、その楽曲における拍点、テンポ、小節線の位置及びコード進行を検出する。音響信号分析装置１０は、図１に示すように、入力操作子１１、コンピュータ部１２、表示器１３、記憶装置１４、外部インターフェース回路１５及びサウンドシステム１６を備えており、これらがバスＢＳを介して接続されている。
【００１３】
入力操作子１１は、オン・オフ操作に対応したスイッチ（例えば数値を入力するためのテンキー）、回転操作に対応したボリューム又はロータリーエンコーダ、スライド操作に対応したボリューム又はリニアエンコーダ、マウス、タッチパネルなどから構成される。これらの操作子は、演奏者の手によって操作されて、分析対象の楽曲の選択、音響信号の分析開始又は停止、楽曲の再生又は停止（後述するサウンドシステム１６からの出力又は停止）、音響信号の分析に関する各種パラメータの設定などに用いられる。入力操作子１１を操作すると、その操作内容を表す操作情報が、バスＢＳを介して、後述するコンピュータ部１２に供給される。
【００１４】
コンピュータ部１２は、バスＢＳにそれぞれ接続されたＣＰＵ１２ａ、ＲＯＭ１２ｂ及びＲＡＭ１２ｃからなる。ＣＰＵ１２ａは、詳しくは後述する音響信号分析プログラム及びそのサブルーチンをＲＯＭ１２ｂから読み出して実行する。ＲＯＭ１２ｂには、音響信号分析プログラム及びそのサブルーチンに加えて、初期設定パラメータ、表示器１３に表示される画像を表わす表示データを生成するための図形データ及び文字データなどの各種データが記憶されている。ＲＡＭ１２ｃには、音響信号分析プログラムの実行時に、各種データが一時的に記憶される。
【００１５】
表示器１３は、液晶ディスプレイ（ＬＣＤ）によって構成される。コンピュータ部１２は、図形データ、文字データなどを用いて表示すべき内容を表わす表示データを生成して表示器１３に供給する。表示器１３は、コンピュータ部１２から供給された表示データに基づいて画像を表示する。例えば分析対象の楽曲の選択時には、楽曲のタイトルリストが表示される。また、例えば分析終了時には、拍点及び小節線を表わすグラフ、テンポの推移を表わすグラフ、コード進行を表わすコード名の系列などが表示される。
【００１６】
また、記憶装置１４は、ＨＤＤ、ＦＤＤ、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤなどの大容量の不揮発性記録媒体と、同各記録媒体に対応するドライブユニットから構成されている。記憶装置１４には、複数の楽曲をそれぞれ表わす複数の楽曲データが記憶されている。楽曲データは、楽曲を所定のサンプリング周期（例えば１／４４１００秒）でサンプリングして得られた複数のサンプル値からなり、各サンプル値が記憶装置１４における連続するアドレスに順に記録されている。楽曲のタイトルを表わすタイトル情報、楽曲データの容量を表わすデータサイズ情報なども楽曲データに含まれている。楽曲データは予め記憶装置１４に記憶されていてもよいし、後述する外部インターフェース回路１５を介して外部機器から取り込んでもよい。
【００１７】
外部インターフェース回路１５は、音響信号分析装置１０を電子音楽装置、パーソナルコンピュータなどの外部機器に接続可能とする接続端子を備えている。音響信号分析装置１０は、外部インターフェース回路１５を介して、ＬＡＮ（Ｌｏｃａｌ  Ａｒｅａ  Ｎｅｔｗｏｒｋ）、インターネットなどの通信ネットワークにも接続可能である。
【００１８】
サウンドシステム１６は、楽曲データをアナログ音信号に変換するＤ／Ａ変換器、変換したアナログ音信号を増幅するアンプ、及び増幅されたアナログ音信号を音響信号に変換して出力する左右一対のスピーカを備えている。ユーザが入力操作子１１を用いて分析対象の楽曲の再生を指示すると、ＣＰＵ１２ａは、分析対象の楽曲データをサウンドシステム１６に供給する。これにより、ユーザは分析対象の楽曲を試聴できる。
【００１９】
つぎに、音響信号分析装置１０の動作の概略を説明する。本実施形態においては、分析対象の楽曲は複数のフレームｔ（＝０，１，・・・）に分割される。各フレームｔにおける、直前の小節線からの拍数として表わされた拍子位置ｓ（例えば、４拍子の楽曲においては、「ｓ」は、「１」、「２」、「３」又は「４」）、次の拍までのフレーム数ｎ及び拍周期ｂの値（テンポの逆数に比例する値）の値の組み合わせの系列によって表わされたモデルＺ（ｏ）（図８参照）と、各フレームｔにおいて発音されているコード（和音）ｉ（＝「Ｃ」、「Ｃｍ」、「Ｄ」、・・・）の系列によって表されたモデルＺ（ｃ）（図８参照）と、に基づいて、拍点、テンポ、小節線の位置（具体的には各フレームの拍子位置）及びコード進行が、ベイズ推定法を用いて同時（一体的）に推定される。なお、詳しくは後述するように、コードの遷移は、拍子位置に関連している。したがって、上記の両モデルは互いに関連付けられている。ただし、ベイズ推定を解析的に実行することは困難であるので、本実施形態では、変分ベイズ法を用いる。具体的には、下記の式（１）に示すような変分事後分布の積として表わされる分布を真の事後分布に近似させる。つまり、前記変分事後分布の積として表わされる分布と真の事後分布との差が収束するまで、各変分事後分布を反復的に更新する。なお、拍周期ｂは、フレームの数によって表わされる。したがって、拍周期ｂの値は「１≦ｂ≦ｂｍａｘ」を満たす整数であり、拍周期ｂの値が「η」である状態では、フレーム数ｎの値は「０≦ｎ＜η」を満たす整数である。
【数１】
2015200803.tif 000003
【００２０】
また、上記の式（１）における「θ（ｃ）」は、コード進行のモデル（つまり、モデルＺ（ｃ））に関連するパラメータである。具体的には、「θ（ｃ）」は、詳しくは後述する、各コードｉに対応したコード特徴量ｘ（ｃ）を規定するパラメータ、及び各コードｉに対応したテンプレートである。また、「θ（ｏ）」は、拍点のモデル（つまり、モデルＺ（ｏ））に関連するパラメータである。具体的には、「θ（ｏ）」は、詳しくは後述する、オンセット特徴量ｘ（ｏ）の分布（正規分布）を規定するパラメータである。また、「θ（ｂ）」は、テンポのモデル（つまり、モデルＺ（ｏ））に関連するパラメータである。具体的には、「θ（ｂ）」は、詳しくは後述する、拍周期ｂに対応したテンプレートである。
【００２１】
つぎに、音響信号分析装置１０の動作について具体的に説明する。ユーザが音響信号分析装置１０の図示しない電源スイッチをオンにすると、ＣＰＵ１２ａは、図２Ａ及び図２Ｂに示す音響信号分析プログラムをＲＯＭ１２ｂから読み出して実行する。なお、図２Ｂにおいては、「判断」のステップを六角形で示す。
【００２２】
ＣＰＵ１２ａは、ステップＳ１０にて音響信号分析処理を開始し、ステップＳ１１にて、記憶装置１４に記憶されている複数の楽曲データにそれぞれ含まれるタイトル情報を読み込んで、楽曲のタイトルをリスト形式で表示器１３に表示する。ユーザは、入力操作子１１を用いて、表示器１３に表示された楽曲の中から分析対象の楽曲データを選択する。なお、ステップＳ１１にて分析対象の楽曲データを選択する際、選択しようとする楽曲データが表す楽曲の一部又は全部を再生して楽曲データの内容を確認できるように構成してもよい。
【００２３】
つぎに、ＣＰＵ１２ａは、ステップＳ１２にて、音響信号分析のための初期設定処理を実行する。具体的には、前記選択された楽曲データのデータサイズ情報に応じた記憶領域をＲＡＭ１２ｃ内に確保し、前記確保した記憶領域に前記選択された楽曲データを読み込む。また、後述するコード特徴量ｘ（ｃ）、オンセット特徴量ｘ（ｏ）、ＢＰＭ特徴量ｘ（ｂ）などを記憶する記憶領域をＲＡＭ１２ｃ内に確保する。また、後述する変分事後分布を計算するために一時的に利用される変数を記憶する記憶領域をＲＡＭ１２ｃ内に確保する。また、同ステップＳ１２において、ユーザは、前記選択した楽曲の調ｋｅｙを、入力操作子１１を用いて入力する。つまり、本実施形態においては、前記選択した楽曲の調ｋｅｙが既知である。
【００２４】
ＣＰＵ１２ａは、ステップＳ１３にて、図３に示すように、前記選択された楽曲を所定の時間間隔をおいて区切り、複数のフレームｔ（＝０，１，・・・）に分割する。各フレームの長さは共通である。
【００２５】
次に、ＣＰＵ１２ａは、ステップＳ１４にて、各フレームｔの各特徴量を計算する。具体的には、ＣＰＵ１２ａは、コード（和音）に関する特徴を表わすコード特徴量ｘ（ｃ）、拍の存在に関する特徴を表すオンセット特徴量ｘ（ｏ）、及びテンポに関する特徴を表すＢＰＭ（ｂｅａｔｓ  ｐｅｒ  ｍｉｎｕｔｅ（１分間あたりの拍数））特徴量ｘ（ｂ）をフレームｔごとに計算する。
【００２６】
フレームｔのコード特徴量ｘ（ｃ）（ｔ）は次のようにして計算される。まず、ＣＰＵ１２ａは、各フレームｔの各周波数ビンのパワーを、その周波数に最も近い音高の周波数（例えば平均律における各音高の基本周波数）にマッピングする。次に、上記のようにして各音高にマッピングされたパワーのうち、低音域（例えば「Ｂ１」以下）に属するパワーをピッチクラス（Ｃ，Ｃ＃，Ｄ，・・・，Ｂ＃）ごとに加算（又は積算）する。このようにして計算された各ピッチクラスのパワーからなる１２次元の特徴量をベース特徴量ＨＰＣＰ（Ｂ）と呼ぶ（図４参照）。また、各音高にマッピングされたパワーのうち、高音域（例えば「Ｃ２」以上）に属するパワーをピッチクラス（Ｃ，Ｃ＃，Ｄ，・・・，Ｂ＃）ごとに加算（又は積算）する。このようにして計算された各ピッチクラスのパワーからなる１２次元の特徴量をトレブル特徴量ＨＰＣＰ（Ｔ）と呼ぶ。各フレームｔに関するベース特徴量ＨＰＣＰ（Ｂ）及びトレブル特徴量ＨＰＣＰ（Ｔ）からなる２４次元の特徴量がコード特徴量ｘ（ｃ）（ｔ）である。
【００２７】
フレームｔのオンセット特徴量ｘ（ｏ）（ｔ）は、次のようにして計算される。ＣＰＵ１２ａは、まず、フレームｔに対して短時間フーリエ変換を実行し、各周波数ビンの信号強度を計算する。次に、ＣＰＵ１２ａは、メルフィルタバンクを用いて、各周波数帯域ｆｂｙ（例えば、ｙ＝１，２，・・・，２０）の信号強度Ｍ（ｆｂｙ，ｔ）を計算する。次に、ＣＰＵ１２ａは、フレーム間における各周波数帯域の信号強度の増加量Ｒ（ｆｂｙ，ｔ）を計算する。下記の式（２）に示すように、フレーム間における前記各周波数帯域の信号強度の増加量の総和がオンセット特徴量ｘ（ｏ）（ｔ）である。
【数２】
2015200803.tif 000004
【００２８】
フレームｔのＢＰＭ特徴量ｘ（ｂ）（ｔ）は、次のようにして計算される。ＣＰＵ１２ａは、まず、オンセット特徴量ｘ（ｏ）（０），ｘ（ｏ）（１）・・・をこの順にフィルタバンクＦＢＢ（図５参照）に入力する。フィルタバンクＦＢＢは、拍周期ｂの値に応じてそれぞれ設けられた複数のコムフィルタＣＦｂからなる。コムフィルタＣＦｂは、１つのデータが入力される度に１つのデータを出力する。コムフィルタＣＦｂは、過去の出力データを拍周期ｂの値に応じた個数だけ記憶するＦＩＦＯ（＝Ｆｉｒｓｔ  Ｉｎ  Ｆｉｒｓｔ  Ｏｕｔ）メモリを有しており、入力されたデータと前記記憶手段に記憶されているデータのうちの最古のデータを所定の比率（例えば、１：１（すなわち、λ＝０．５））で加算して出力する。オンセット特徴量ｘ（ｏ）の系列ｘ（ｏ）（ｔ）｛＝ｘ（ｏ）（０），ｘ（ｏ）（１）・・・｝をフィルタバンクＦＢＢに入力することにより得られたデータｘｂ（Ｄ）の系列ｘｂ（Ｄ）（ｔ）｛＝ｘｂ（Ｄ）（０），ｘｂ（Ｄ）（１）・・・の時系列を逆にして、フィルタバンクＦＢＢに再度入力することにより、拍周期ｂに関するＢＰＭ特徴量の系列ｘｂ（ｂ）（ｔ）｛＝ｘｂ（ｂ）（０），ｘｂ（ｂ）（１）・・・｝が得られる。フレームｔのＢＰＭ特徴量ｘ（ｂ）（ｔ）は、拍周期ｂごとに計算されたＢＰＭ特徴量ｘｂ＝１，２・・・（ｂ）（ｔ）の集合として表わされる（図６参照）。
【００２９】
次に、ＣＰＵ１２ａは、ステップＳ１５にて、各フレームｔの各特徴量の観測尤度を計算する。具体的には、ＣＰＵ１２ａは、コード特徴量ｘ（ｃ）（ｔ）の観測尤度Ｌ（ｃ）（ｔ）、オンセット特徴量ｘ（ｏ）（ｔ）の観測尤度Ｌ（ｏ）（ｔ）及びＢＰＭ特徴量ｘ（ｂ）（ｔ）の観測尤度Ｌ（ｂ）（ｔ）を計算する。
【００３０】
コード特徴量ｘ（ｃ）（ｔ）の観測尤度Ｌ（ｃ）（ｔ）は下記の式（３）のように表わされる。
【数３】
2015200803.tif 000005
【００３１】
対数観測尤度Ｌｉ（ｃ）（ｔ）は、トレブル特徴量ＨＰＣＰ（Ｔ）及びベース特徴量ＨＰＣＰ（Ｂ）のテンプレートＴＭＰｉ（Ｔ）及びテンプレートＴＭＰｉ（Ｂ）に対する適合度に相当する。テンプレートＴＭＰｉ（Ｔ）及びテンプレートＴＭＰｉ（Ｂ）は、コードｉに対応する係数列である。例えば、「Ｃｍａｊ」に対応する係数列ＴＭＰｉ＝Ｃｍａｊ（Ｔ）は、｛１，０，０，０，１，０，０，１，０，０，０，０｝と表わされる。また、「Ｃｍｉｎ」に対応する係数列ＴＭＰｉ＝Ｃｍｉｎ（Ｔ）は、｛１，０，０，１，０，０，０，１，０，０，０，０｝と表わされる。なお、テンプレートＴＭＰｉ（Ｂ）は、テンプレートＴＭＰｉ（Ｔ）よりもルート音に対応する成分が強調されている。トレブル特徴量ＨＰＣＰ（Ｔ）と係数列ＴＭＰｉ（Ｔ）とのコサイン距離、及びベース特徴量ＨＰＣＰ（Ｂ）と係数列ＴＭＰｉ（Ｂ）とのコサイン距離の重み付け和が、対数観測尤度Ｌｉ（ｃ）（ｔ）である。例えば、トレブル特徴量ＨＰＣＰ（Ｔ）と係数列ＴＭＰｉ（Ｔ）とのコサイン距離、及びベース特徴量ＨＰＣＰ（Ｂ）と係数列ＴＭＰｉ（Ｂ）とのコサイン距離の重み係数をそれぞれ「１．０」及び「１．０」に設定する。ベース音を重視する場合には、前記重み係数をそれぞれ「１．０」及び「２．０」に設定すればよい。なお、「ｚｉ（ｃ）（ｔ）」は、次に説明するような２値変数である。すなわち、この変数は、フレームｔにおいて、コードが「ｉ」である場合に「１」であり、それ以外の場合に「０」である。
【００３２】
オンセット特徴量ｘ（ｏ）（ｔ）及びＢＰＭ特徴量ｘ（ｂ）（ｔ）の観測尤度Ｌ（ｏ）（ｔ）及び観測尤度Ｌ（ｂ）（ｔ）は下記の式（４）及び式（５）のように表わされる。なお、式（４）及び式（５）における「：」は、その次元に関する総和を意味する（下記の式（６）参照）。
【数４】
2015200803.tif 000006
【数５】
2015200803.tif 000007
【数６】
2015200803.tif 000008
【００３３】
本実施形態においては、オンセット特徴量ｘ（ｏ）（ｔ）は、次の拍点までのフレーム数ｎの値に応じて設定された正規分布に従うものとする。つまり、オンセット特徴量ｘ（ｏ）（ｔ）の観測尤度Ｌ（ｏ）（ｔ）は、次の拍点までのフレーム数ｎの値に応じて設定された正規分布の確率変数としてオンセット特徴量ｘ（ｏ）（ｔ）を代入することにより計算される。例えば、フレーム数ｎの値が「０」であるときは、平均値が「３」であって、且つ分散が「１」である正規分布が用いられる。また、フレーム数ｎの値が「０」でないとき、平均値が「０」であって、且つ分散が「１」である正規分布が用いられる。上記の両正規分布を規定する平均値及び分散が上記式（１）における「θ（ｏ）」である。フレーム数ｎの値が「０」であるときの正規分布を規定するパラメータが式（４）における「θ０（ｏ）」である。フレーム数ｎの値が「０」でないときの正規分布を規定するパラメータが式（４）における「θ１（ｏ）」である。なお、「ｚｓ，ｎ，ｂ（ｏ）（ｔ）」は、次に説明するような２値変数である。すなわち、この変数は、フレームｔにおいて、拍子位置が「ｓ」であり、次の拍点までのフレーム数が「ｎ」であって、かつ拍周期が「ｂ」である場合に「１」であり、それ以外の場合に「０」である。
【００３４】
また、ＢＰＭ特徴量ｘ（ｂ）（ｔ）の観測尤度Ｌ（ｂ）（ｔ）は、拍周期ｂごとに設けられたテンプレートＴＭＰ（ｂ）に対するＢＰＭ特徴量ｘ（ｂ）（ｔ）の適合度に相当する。テンプレートＴＭＰ（ｂ）は、ＢＰＭ特徴量ｘ（ｂ）（ｔ）を構成するＢＰＭ特徴量ｘｂ（ｂ）（ｔ）にそれぞれ乗算される係数ζｂ，γ｛＝１，２・・・｝の系列からなる（図７参照）。テンプレートＴＭＰｂを構成する係数ζｂ，γのうち、インデックスγが拍周期ｂに等しい係数及び拍周期ｂの整数倍に等しい係数が極大となるように、テンプレートＴＭＰｂが設定されている。式（５）に示すように、テンプレートＴＭＰｂとＢＰＭ特徴量ｘ（ｂ）（ｔ）の内積が、観測尤度Ｌ（ｂ）（ｔ）である。なお、この演算式における「μｂ」は、オンセット特徴量ｘ（ｏ）（ｔ）に対するＢＰＭ特徴量ｘ（ｂ）（ｔ）の重みを決定する係数である。つまり、「μｂ」を大きく設定するほど、結果的に、ＢＰＭ特徴量ｘ（ｂ）（ｔ）が重視される。また、この演算式におけるＺ（μｂ）は、「μｂ」に依存する正規化係数である。なお、上記のテンプレートＴＭＰｂが上記式（１）における「θ（ｂ）」に相当する。
【００３５】
さらに、ＣＰＵ１２ａは、同ステップＳ１５において、式（７）に示す、観測尤度Ｌ（ｏ，ｂ）（ｔ）を計算する。観測尤度Ｌ（ｏ，ｂ）（ｔ）は、観測尤度Ｌ（ｏ）（ｔ）と観測尤度Ｌ（ｂ）（ｔ）との積である
【数７】
2015200803.tif 000009
【００３６】
次に、ＣＰＵ１２ａはステップＳ１６〜ステップＳ２７にて拍点、テンポ、小節線の位置及びコード進行の推定処理を実行する。ここで、コードの遷移（コード進行）のモデルＺ（ｃ）は、下記の式（８）に示すように定義される。
【数８】
2015200803.tif 000010
【００３７】
式（８）における、「τｉ，ｊ（ｃ，０）」は、いずれか一方のフレームが拍点ではないフレーム間においてコードｉからコードｊへ遷移する確率を表わす。一方、「τｉ，ｊ（ｃ，ｓ）」は、拍子位置ｓへの遷移が生起するフレーム間（つまり、遷移先のフレームに拍点が存在する場合）において、コードｉからコードｊへ遷移する確率を表わす。コードの遷移は拍点において生起する確率が高い。言い換えれば、拍点でないフレーム間においては、同一のコードに遷移する確率が高い。したがって、「τｉ，ｊ＝ｉ（ｃ，０）」は「τｉ，ｊ≠ｉ（ｃ，０）」よりも大きい。
【００３８】
一般に、コード間の遷移確率は、拍点の拍子位置ｓに依存する。例えば、４分の４拍子の楽曲において、４拍目（つまりｓ＝４）の和音が「Ｇ７」であるとき、次の小節の１拍目（つまりｓ＝１）の和音は「Ｃ」である可能性が高い（ドミナントモーション）。また、コード間の遷移確率は、楽曲の調ｋｅｙに依存する。例えば、コード「Ｃ」からコード「Ｆ」への遷移はハ長調の楽曲において生起する可能性が高い。そこで、予め種々の楽曲におけるコード間の遷移回数をカウントし、前記カウントしたコード間の遷移回数が、データベースとしてＲＯＭ１２ｂに記憶されている。データベースは、拍子ごとに構成されている。つまり、音響信号分析装置１０は、例えば、４分の３拍子の楽曲を分析するときに用いるデータベース、４分の４拍子の楽曲を分析するときに用いるデータベース、８分の６拍子の楽曲を分析するときに用いるデータベースなどを備える。そして、各データベースには、コード間の遷移回数が、調ｋｅｙ及び拍子位置ｓに関連づけて記憶されている。
【００３９】
また、拍点、テンポ及び拍子位置の遷移のモデルＺ（ｏ）は、下記の式（９）に示すように表わされる。
【数９】
2015200803.tif 000011
【００４０】
式（９）における、「τｓ０，ｎ０，ｂ０，ｓ１，ｎ１，ｂ１（ｏ）」は、拍子位置ｓ０、フレーム数ｎ０、且つ拍周期ｂ０である状態から、拍子位置ｓ１、フレーム数ｎ１、且つ拍周期ｂ１である状態へ遷移する確率を表わす。本実施形態では、フレーム数ｎの値が「０」でない状態（ｎ０≠０）から次の状態へ遷移するとき、拍子位置ｓ及び拍周期ｂの値は変更されず、フレーム数ｎの値が「１」だけ減少する。一方、フレーム数ｎの値が「０」である状態（ｎ０＝０）から次の状態へ遷移するとき、拍子位置ｓ及び拍周期ｂの値は変更され得る。このとき、フレーム数ｎの値は、遷移後の拍周期ｂの値より「１」だけ小さい値に設定される。以下、拍子位置ｓの遷移確率を「τ（ｓ）」と表記し、拍周期ｂの遷移確率を「τ（ｂ）」と表記する。
【００４１】
また、各遷移確率の事前分布は下記の式（１０）乃至（１２）に示すように定義される。
【数１０】
2015200803.tif 000012
【数１１】
2015200803.tif 000013
【数１２】
2015200803.tif 000014
【００４２】
各遷移確率は多項分布に従うので、本実施形態では、多項分布の共役事前分布であるディリクレ分布を採用した。式（１０）における「νｉ，ｊ（ｃ）」は、上記のコード進行に関するデータベースに相当する。つまり、例えば、「νｉ，ｊ（ｃ）」は、「４分の４拍子の楽曲であって、且つハ長調の楽曲において、４拍目から１拍目への遷移するとき、「Ｇ７」から「Ｃｍａｊ」への遷移がＭ回あった」という情報に相当する。また、式（１１）における「νｓ０，ｓ１（ｓ）」は、拍子位置ｓ０から拍子位置ｓ１へ遷移する度合いを設定するパラメータである。例えば、「ν１，２（ｓ）」、「ν２，３（ｓ）」、「ν３，４（ｓ）」及び「ν４，１（ｓ）」が「１」に近い値（例えば「０．９９」）に設定され、その他の遷移に関するパラメータは充分に小さい値に設定される。また、式（１２）における「νｂ０，ｂ１（ｂ）」は、拍周期ｂ０から拍周期ｂ１へ遷移する度合いを設定するパラメータである。例えば、「νｂ０，ｂ１（ｂ）」を、所定の係数ρ，ωを用いて表わされた関数ρ×ｅｘｐ（−ω（ｂ０−ｂ１）２）と設定することにより、テンポが急激に変化することを抑制できる。係数ρ及び係数ωは、例えば、「１０．０」及び「５．０」にそれぞれ設定される。
【００４３】
モデルＺ（ｃ）の変分事後分布ｑ（Ｚ（ｃ））は、下記の式（１３）のように表わされる。なお、以降、「ｘ」の分布が「ｐ（ｘ）」に従う場合の「ｆ（ｘ）」の期待値を「〈ｆ（ｘ）〉ｐ（ｘ）」と表記する。ただし、標記が煩雑になるのを防止するため、「ｐ（ｘ）」を省略する場合もある。
【数１３】
2015200803.tif 000015
【００４４】
式（１３）における「〈ｌｏｇτｉ，ｊ（ｃ，０）〉」は、式（１４）に示すように表される。ただし、「ψ（ｘ）」は、ｄｉｇａｍｍａ関数である。
【数１４】
2015200803.tif 000016
【００４５】
また、ここで、下記の式（１５）及び式（１６）に示す量を定義しておく。式（１５）は、モデルＺ（ｃ）の各フレームの対数観測尤度に相当し、式（１６）は、モデルＺ（ｃ）の各状態間の対数遷移確率に相当する。
【数１５】
2015200803.tif 000017
【数１６】
2015200803.tif 000018
【００４６】
一方、モデルＺ（ｏ）の変分事後分布ｑ（Ｚ（ｏ））は、下記の式（１７）のように表わされる。
【数１７】
2015200803.tif 000019
【００４７】
なお、式（１７）における「ξｉ，ｊ（ｃ）（ｔ）」は、下記の式（１８）に示すように、フレームｔへの遷移において、コードが「ｉ」から「ｊ」へ遷移した回数の期待値を意味する。
【数１８】
2015200803.tif 000020
【００４８】
ここで、下記の式（１９）及び式（２０）に示す量を定義しておく。式（１９）は、モデルＺ（ｏ）の各フレームの対数観測尤度に相当し、式（２０）は、モデルＺ（ｏ）の各状態間の対数遷移確率に相当する。
【数１９】
2015200803.tif 000021
【数２０】
2015200803.tif 000022
【００４９】
なお、式（２０）における「δ（ｎ）」は、Ｋｒｏｎｅｃｋｅｒのデルタ関数である。
【００５０】
上記のように、変分事後分布ｑ（Ｚ（ｃ））は、モデルＺ（ｏ）の期待値に依存し、変分事後分布ｑ（Ｚ（ｏ））は、モデルＺ（ｃ）の期待値に依存している。したがって、変分事後分布ｑ（Ｚ（ｃ））を更新した後、その期待値に依存する変分事後分布ｑ（Ｚ（ｏ））を更新する必要があり、変分事後分布ｑ（Ｚ（ｏ））を更新した後、その期待値に依存する変分事後分布ｑ（Ｚ（ｃ））を更新する必要がある。このように、変分事後分布ｑ（Ｚ（ｃ））及び変分事後分布ｑ（Ｚ（ｏ））を交互に更新することにより、両変分事後分布を収束させる。
【００５１】
具体的には、まず、ＣＰＵ１２ａは、ステップＳ１６にて、推論対象をモデルＺ（ｃ）に設定する。つぎに、ＣＰＵ１２ａは、ステップＳ１７にて、現在の推論対象のモデルがモデルＺ（ｃ）か否かを判定する。現在の推論対象のモデルがモデルＺ（ｃ）である場合には、ＣＰＵ１２ａは「Ｙｅｓ」と判定して、ステップＳ１８にて、式（１５）及び式（１６）に基づいて、対数観測尤度Ｏｉ（ｃ）（ｔ）及び対数遷移確率Ｔｉ，ｊ（ｃ）（ｔ）を計算する。つぎに、ＣＰＵ１２ａは、ステップＳ１９にて、前向き後ろ向きアルゴリズムを用いて、前向き変数αｉ（ｃ）（ｔ）及び後ろ向き変数βｉ（ｃ）（ｔ）を計算する。具体的には、ＣＰＵ１２ａは、下記の式（２１）及び式（２２）に基づいて、前向き変数αｉ（ｃ）（ｔ）及び後ろ向き変数βｉ（ｃ）（ｔ）を計算する。なお、式（２１）における「Ｎｉ（ｃ）（ｔ）」は、式（２３）に示すような正規化係数である。
【数２１】
2015200803.tif 000023
【数２２】
2015200803.tif 000024
【数２３】
2015200803.tif 000025
【００５２】
次に、ＣＰＵ１２ａは、ステップＳ２０にて、下記の式（２４）に基づいて、各フレームｔのモデルＺ（ｃ）の期待値を計算する。
【数２４】
2015200803.tif 000026
【００５３】
また、ＣＰＵ１２ａは、同ステップＳ２０にて、下記の式（２５）に基づいて、フレームｔにおいてコードｉからコードｊに遷移した回数の期待値を計算する。
【数２５】
2015200803.tif 000027
【００５４】
次に、ＣＰＵ１２ａは、ステップＳ２１にて、変分事後分布ｑ（Ｚ（ｏ））及び変分事後分布ｑ（Ｚ（ｃ））が収束したか否かを判定する。変分事後分布ｑ（Ｚ（ｏ））及び変分事後分布ｑ（Ｚ（ｃ））のうち少なくとも一方が収束していない場合には、ＣＰＵ１２ａは「Ｎｏ」と判定して、ステップＳ２２にて、推論対象を切り替え、ステップＳ１７に処理を進める。すなわち、ステップＳ２２において、現在の推論対象がモデルＺ（ｃ）である場合には、ＣＰＵ１２ａは、推論対象をモデルＺ（ｏ）に設定する。一方、ステップＳ２２において、現在の推論対象がモデルＺ（ｏ）である場合には、ＣＰＵ１２ａは、推論対象をモデルＺ（ｃ）に設定する。
【００５５】
ステップＳ１７において、推論対象がモデルＺ（ｏ）である場合には、ＣＰＵ１２ａは「Ｎｏ」と判定して、ステップＳ２３にて、式（１９）及び式（２０）に基づいて、対数観測尤度Ｏｓ，ｎ，ｂ（ｏ）（ｔ）」及び対数遷移確率Ｔｓ０，ｎ０，ｂ０ｓ１，ｎ１，ｂ１（ｏ）（ｔ）を計算する。つぎに、ＣＰＵ１２ａは、ステップＳ２４にて、前向き後ろ向きアルゴリズムを用いて、前向き変数αｓ，ｎ，ｂ（ｏ）（ｔ）及び後ろ向き変数βｓ，ｎ，ｂ（ｏ）（ｔ）を計算する。具体的には、ＣＰＵ１２ａは、下記の式（２６）及び式（２７）に基づいて、前向き変数αｓ，ｎ，ｂ（ｏ）（ｔ）及び後ろ向き変数βｓ，ｎ，ｂ（ｏ）（ｔ）を計算する。なお、式（２６）における「Ｎｓ，ｎ，ｂ（ｏ）（ｔ）」は、式（２８）に示すような正規化係数である。
【数２６】
2015200803.tif 000028
【数２７】
2015200803.tif 000029
【数２８】
2015200803.tif 000030
【００５６】
次に、ＣＰＵ１２ａは、ステップＳ２５にて、下記の式（２９）に基づいて、各フレームｔのモデルＺ（ｏ）の期待値を計算し、ステップＳ２１に処理を進める。
【数２９】
2015200803.tif 000031
【００５７】
ステップＳ２１において、変分事後分布ｑ（Ｚ（ｏ））及び変分事後分布ｑ（Ｚ（ｃ））のいずれもが収束した場合には、ＣＰＵ１２ａは「Ｙｅｓ」と判定して、ステップＳ２６にて、各フレームｔにおけるモデルＺ（ｏ）及びモデルＺ（ｃ）の期待値がそれぞれ最大になる系列を選択することにより、分析対象の楽曲における拍点、テンポ、小節線の位置及びコード進行を推定する。ＣＰＵ１２ａは、前記推定結果を例えば図８に示すような形式で表示する。そして、ＣＰＵ１２ａは、ステップＳ２７にて、音響信号分析処理を終了する。
【００５８】
音響信号分析装置１０においては、モデルＺ（ｃ）の各状態間の遷移確率τｉ、ｊ（ｃ）が、拍子位置ｓに関連して設定されている。そして、モデルＺ（ｃ）の事後分布と、モデルＺ（ｏ）の事後分布とが、同時（一体的）に計算される。モデルＺ（ｃ）とモデルＺ（ｏ）とは、拍子位置ｓに関して互いに依存しており、この両モデルの事後分布を同時（一体的）に計算することにより、楽曲における拍点、テンポ、拍子位置、及びコードに関する情報が精度良く得られる。言い換えれば、一方のモデルの推定結果が他方のモデルの推定結果にフィードバックされて推定結果が更新される（Ｓ１６〜Ｓ２６）。これにより、両モデルの推定精度が高められる。したがって、音響信号分析装置１０によれば、拍点、テンポ、小節線の位置及びコード進行の推定精度を従来よりも向上させることができる。
【００５９】
さらに、本発明の実施にあたっては、上記実施形態に限定されるものではなく、本発明の目的を逸脱しない限りにおいて種々の変更が可能である。
【００６０】
例えば、上記実施形態では、楽曲全体を分析対象としているが、楽曲の一部（例えば数小節）のみを分析対象としてもよい。この場合、入力した楽曲データのうち、分析対象とする部分を選択可能に構成するとよい。また、楽曲のうちの単一のパート（例えばリズムセクション）のみを分析対象としてもよい。
【００６１】
また、上記実施形態では、楽曲の調が既知であると仮定しているが、楽曲の調も、拍点、テンポ、小節線の位置及びコード進行と同時に推定しても良い。例えば、コードの遷移確率τ（ｃ）を２４個のブロックから構成されるブロック対角行列とし，各ブロックにそれぞれの調におけるコード進行のデータを対応させるとよい。具体的には、コードの種類を４つとすると、ブロック対角行列は、９６×９６（９６＝４（コードの種類）×１２×２（１２個の調のルート×（長調または短調）））の要素を有する。ブロック対角行列の要素（１，１）から要素（４，４）は、調が「Ｃｍａｊ」である場合のコードの遷移を表わす。要素（５，５）から要素（８，８）は、調が「Ｃｍｉｎ」である場合のコードの遷移を表わす。要素（９，９）から要素（１２，１２）は、調が「Ｃ＃ｍａｊ」である場合のコードの遷移を表わす。要素（１３，１３）から要素（１６，１６）は、調が「Ｃ＃ｍｉｎ」である場合のコードの遷移を表わす。他の要素についても、上記の要素と同様に設定することで、すべての調のコード進行データを網羅できる。この場合、ブロックをまたがるような遷移は許されないので、単一の調におけるコード進行を推定することに相当する。
【００６２】
また、例えば、ユーザによって指定されたフレームの拍点、テンポ、拍子位置、及びコードのいずれか又は複数の要素を修正可能に構成しても良い。この場合、ユーザが入力操作子１１を用いて修正値を入力すると、ＣＰＵ１２ａは、前記指定されたフレームの各種要素の値を前記入力された値に修正するとともに、そのフレームに近接するフレームの各種要素を前記入力された値に応じて自動的に修正するように構成するとよい。例えば、連続する複数のフレームのテンポの推定値が同じ値であって、そのうちの１つのフレームのテンポの値が修正されたとき、前記複数のフレームのテンポの値を前記１つのフレームの修正値と同じ値に自動的に修正してもよい。
【００６３】
また、例えば、コード進行の推定において、優先的に推定するコードを指定可能に構成してもよい。例えば、メジャーコードを優先的に推定するように設定可能に構成してもよい。この場合、メジャーコード以外の対数観測尤度を十分に小さく設定する。これにより、メジャーコードが優先的に推定される。これによれば、楽曲の大凡のコードが既知である場合、コードの推定精度を向上させることができる。
【００６４】
また、例えば、テンポの推定において、優先的に推定するテンポの範囲を指定可能に構成してもよい。具体的には、「Ｐｒｅｓｔｏ」、「Ｍｏｄｅｒａｔｏ」などのテンポを表わす用語を表示して、優先的に推定するテンポの範囲を選択可能に構成してもよい。例えば、「Ｐｒｅｓｔｏ」が選択された場合、ＢＰＭ＝１６０〜１９０の範囲以外の対数観測尤度を十分に小さく設定する。これにより、ＢＰＭ＝１６０〜１９０の範囲のテンポが優先的に推定される。これによれば、楽曲の大凡のテンポが既知である場合、テンポの推定精度を向上させることができる。
【符号の説明】
【００６５】
１０・・・音響信号分析装置、ＨＰＣＰ（Ｂ）・・・ベース特徴量、ＨＰＣＰ（Ｔ）・・・トレブル特徴量、Ｌｉ（ｃ），Ｌ（ｏ），Ｌ（ｂ），Ｌ（ｏ，ｂ）・・・観測尤度、Ｏｉ（ｃ），Ｏｓ，ｎ，ｂ（ｏ）・・・対数観測尤度、Ｔｉ，ｊ（ｃ），Ｔｓ０，ｎ０，ｂ０ｓ１，ｎ１，ｂ１（ｏ）・・・対数遷移確率、Ｚ（ｃ），Ｚ（ｏ）・・・モデル、ｂ・・・拍周期、ｉ・・・コード、ｋｅｙ・・・調、ｎ・・・フレーム、ｑ（Ｚ（ｃ）），ｑ（Ｚ（ｏ））・・・変分事後分布、ｓ・・・拍子位置、ｔ・・・フレーム、ｘ（ｏ）・・・オンセット特徴量、ｘ（ｃ）・・・コード特徴量
【図１】
2015200803.tif 000032
【図２Ａ】
2015200803.tif 000033
【図２Ｂ】
2015200803.tif 000034
【図３】
2015200803.tif 000035
【図４】
2015200803.tif 000036
【図５】
2015200803.tif 000037
【図６】
2015200803.tif 000038
【図７】
2015200803.tif 000039
【図８】
2015200803.tif 000040