(57)【要約】【課題】  入力音響信号の波形がアタック区間で不安定なものとなっても、全体としての基本周波数の推定精度を高める。【解決手段】  アタック検出１ａでは、入力音響信号をフレームに分割してアタック検出を行う。基本周波数の確率密度関数の推定４１では、フレーム毎に、ＥＭアルゴリズムにより基本周波数の確率密度関数を推定する。マルチエージェントモデルによる基本周波数の継続的な追跡４２では、この確率密度関数から基本周波数を推定する。その際、処理対象となるフレームがアタック区間のものか否かにより、基本周波数を得るための処理４１または４２の演算の態様を切り換える制御が行われる。【選択図】図１
2008209579.tif 000002
【特許請求の範囲】
【請求項１】
入力音響信号を所定時間長のフレームに分割し、フレーム毎に入力音響信号がアタック区間の信号であるか否かを判定するアタック検出手段と、
フレーム毎に、各々音源の音の高調波構造に対応した構造を有する確率密度関数である音モデルを使用して、各種の基本周波数に対応した複数の音モデルを重み付け加算した混合分布を構成し、この混合分布が入力音響信号の周波数成分の分布となるように、各音モデルに対する重み値を逐次更新して最適化し、最適化された各音モデルの重み値を前記音源の音の基本周波数の確率密度関数として推定する確率密度関数推定手段と、
前記基本周波数の確率密度関数に基づいて前記入力音響信号に含まれる１または複数の音源の音の基本周波数を推定して出力する基本周波数推定手段と、
前記確率密度関数推定手段の処理対象となるフレームがアタック区間のものであるか否かにより前記確率密度関数推定手段における基本周波数の確率密度関数の推定のための演算または前記基本周波数推定手段における基本周波数の推定のための演算の態様を切り換える演算制御手段と
を具備することを特徴とする音分析装置。
【請求項２】
前記演算制御手段は、前記確率密度関数推定手段の処理対象であるフレームがアタック区間のものである場合には、当該フレームにおける重み値の逐次更新が所定の初期値から開始されるように前記確率密度関数推定手段における基本周波数の確率密度関数の推定のための演算の制御を行い、前記確率密度関数推定手段の処理対象であるフレームがアタック区間のものでない場合には、当該フレームにおける重み値の逐次更新が前フレームにおける重み値の最終値を初期値として開始されるように前記確率密度関数推定手段における基本周波数の確率密度関数の推定のための演算の制御を行うことを特徴とする請求項１に記載の音分析装置。
【請求項３】
前記演算制御手段は、前記確率密度関数推定手段の処理対象であるフレームがアタック区間のものである場合には、当該フレームにおける重み値の逐次更新が所定の初期値から開始されるように前記確率密度関数推定手段における基本周波数の確率密度関数の推定のための演算の制御を行い、前記確率密度関数推定手段の処理対象であるフレームがアタック区間のものでない場合には、当該フレームにおける重み値の逐次更新が前フレームにおける重み値の最終値と所定の初期値とをミキシングした重み値を初期値として開始されるように前記確率密度関数推定手段における基本周波数の確率密度関数の推定のための演算の制御を行い、かつ、直前のアタック区間の終了時からの経過時間が長くなるに従って、前フレームでの重み値の最終値が強調されるように、前フレームにおける重み値の最終値と所定の初期値とのミキシング比を制御することを特徴とする請求項１に記載の音分析装置。
【請求項４】
前記演算制御手段は、前記確率密度関数推定手段の処理対象であるフレームがアタック区間のものであるか否かにより、前記基本周波数の確率密度関数の推定に用いられる音モデルの切り換えを行うことを特徴とする請求項１に記載の音分析装置。
【請求項５】
前記演算制御手段は、前記確率密度関数推定手段の処理対象であるフレームがアタック区間のものでない場合には前記基本周波数推定手段に当該フレームについての基本周波数の推定結果の出力を行わせ、前記確率密度関数推定手段の処理対象であるフレームがアタック区間のものである場合には前記基本周波数推定手段に当該フレームについての基本周波数の推定結果の出力を行わせないことを特徴とする請求項１に記載の音分析装置。
【請求項６】
コンピュータを、
入力音響信号を所定時間長のフレームに分割し、フレーム毎に入力音響信号がアタック区間の信号であるか否かを判定するアタック検出手段と、
フレーム毎に、各々音源の音の高調波構造に対応した構造を有する確率密度関数である音モデルを使用して、各種の基本周波数に対応した複数の音モデルを重み付け加算した混合分布を構成し、この混合分布が入力音響信号の周波数成分の分布となるように、各音モデルに対する重み値を逐次更新して最適化し、最適化された各音モデルの重み値を前記音源の音の基本周波数の確率密度関数として推定する確率密度関数推定手段と、
前記基本周波数の確率密度関数に基づいて前記入力音響信号に含まれる１または複数の音源の音の基本周波数を推定して出力する基本周波数推定手段と、
前記確率密度関数推定手段の処理対象となるフレームがアタック区間のものであるか否かにより前記確率密度関数推定手段における基本周波数の確率密度関数の推定のための演算または前記基本周波数推定手段における基本周波数の推定のための演算の態様を切り換える演算制御手段と
して機能させることを特徴とするコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、市販のＣＤ（ｃｏｍｐａｃｔ  ｄｉｓｃ）などに収録されている、歌声や複数種類の楽器音を同時に含む音楽音響信号を対象に、メロディ音やベース音の音高（本明細書では基本周波数の意味で用いる）を推定する音分析装置およびプログラムに関する。
【背景技術】
【０００２】
多数の音源の音が混ざり合ったモノラルの音響信号中から、ある特定の音源の音高を推定することは、非常に困難である。混合音に対して音高推定することが難しい本質的な理由の１つに、時間周波数領域において、ある音の周波数成分が同時に鳴っている他の音の周波数成分と重複することが挙げられる。例えば、歌声、鍵盤楽器（ピアノ等）、ギター、ベースギター、ドラムス等で演奏される典型的なポピュラー音楽では、メロディを担う歌声の高調波構造の一部（特に基本周波数成分）は、鍵盤楽器、ギターの高調波成分やベースギターの高次の高調波成分、スネアドラム等の音に含まれるノイズ成分などと頻繁に重複する。そのため、各周波数成分を局所的に追跡するような手法は、複雑な混合音に対しては安定して機能しない。基本周波数成分が存在することを前提に高調波構造を推定する手法もあるが、そのような手法は、ミッシングファンダメンタル（ｍｉｓｓｉｎｇ  ｆｕｎｄａｍｅｎｔａｌ）現象を扱えないという大きな欠点を持つ。さらに、同時に鳴っている他の音の周波数成分が基本周波数成分と重複すると、有効に機能しない。
【０００３】
以上のような理由により、従来、単一音のみか、非周期的な雑音を伴った単一音を収録した音響信号を対象とした音高の推定技術はあったが、市販のＣＤに記録された音響信号のように複数の音が混ざり合ったものについて音高を推定する技術はなかった。
【０００４】
しかしながら、近年、統計的手法を利用することにより、混合音に含まれる各音の音高を適切に推定する技術が提案されるに至った。特許文献１の技術である。
【０００５】
この特許文献１の技術では、メロディ音のものと考えられる帯域に属する周波数成分と、ベース音のものと考えられる帯域に属する周波数成分とを入力音響信号からＢＰＦにより別々に取り出し、それらの各帯域の周波数成分に基づき、メロディ音およびベース音の各々の基本周波数の推定を行う。
【０００６】
さらに詳述すると、特許文献１の技術では、音の高調波構造に対応した確率分布を持った音モデルを用意し、メロディ音の帯域の各周波数成分、ベース音の帯域の各周波数成分が、様々な基本周波数に対応した各音モデルを重み付け加算した混合分布であると考える。そして、各音モデルの重みの値をＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）アルゴリズムを用いて推定する。
【０００７】
このＥＭアルゴリズムは、隠れ変数を含む確率モデルに対して最尤推定を行うための反復アルゴリズムであり、局所最適解を求めることができる。ここで、最も大きな重みの値を持つ確率分布は、その時点で最も優勢な高調波構造であるとみなすことができるため、あとはその優勢な高調波構造における基本周波数を音高として求めればよい。この手法は基本周波数成分の存在に依存しないため、ミッシングファンダメンタル現象も適切に扱うことができ、基本周波数成分の存在に依存せずに、最も優勢な高調波構造を求めることができる。
【特許文献１】特許第３４１３６３４号
【特許文献２】特許第３６６０５９９号
【発明の開示】
【発明が解決しようとする課題】
【０００８】
ところで、上述した従来の音分析装置では、入力音響信号を一定時間長のフレームに分割し、フレーム単位でＥＭアルゴリズムを実行し、音源の音の基本周波数の推定を行っていた。また、各フレームにおいて、ＥＭアルゴリズムの繰り返しにより各種の基本周波数の音モデルに対する重み値を更新して最適化する際には、前フレームにおいて推定された重み値の最終値を引き継ぎ、これを初期状態として当該フレームにおけるＥＭアルゴリズムを実行していた。しかしながら、一般に楽音は、アタック区間において波形が不安定なものとなり易い。このため、従来の音分析装置は、波形が不安定なアタック区間の入力音響信号の基本周波数の推定処理を行うときに、推定処理が不安定な状態に陥り、基本周波数の誤推定が連続して発生し易いという問題があった。
【０００９】
この発明は、以上説明した事情に鑑みてなされたものであり、入力音響信号の波形がアタック区間において不安定なものとなる場合においても、全体としての基本周波数の推定精度を高めることができる音分析装置および音分析プログラムを提供することを目的としている。
【課題を解決するための手段】
【００１０】
この発明は、入力音響信号を所定時間長のフレームに分割し、フレーム毎に入力音響信号がアタック区間の信号であるか否かを判定するアタック検出手段と、フレーム毎に、各々音源の音の高調波構造に対応した構造を有する確率密度関数である音モデルを使用して、各種の基本周波数に対応した複数の音モデルを重み付け加算した混合分布を構成し、この混合分布が入力音響信号の周波数成分の分布となるように、各音モデルに対する重み値を逐次更新して最適化し、最適化された各音モデルの重み値を前記音源の音の基本周波数の確率密度関数として推定する確率密度関数推定手段と、前記基本周波数の確率密度関数に基づいて前記入力音響信号に含まれる１または複数の音源の音の基本周波数を推定して出力する基本周波数推定手段と、前記確率密度関数推定手段の処理対象となるフレームがアタック区間のものであるか否かにより前記確率密度関数推定手段における基本周波数の確率密度関数の推定のための演算または前記基本周波数推定手段における基本周波数の推定のための演算の態様を切り換える演算制御手段とを具備することを特徴とする音分析装置並びにコンピュータを前記音分析装置として機能させるコンピュータプログラムを提供する。
【００１１】
かかる発明によれば、確率密度関数推定手段の処理対象となるフレームがアタック区間のものであるか否かにより前記確率密度関数推定手段における基本周波数の確率密度関数の推定のための演算または前記基本周波数推定手段における基本周波数の推定のための演算の態様を切り換える演算制御手段とを設けたので、全体としての基本周波数の推定精度を高めるのに適した適切な演算の態様を選択し、確率密度関数推定手段または基本周波数推定手段に実行させることができ、例えば曲全体としての基本周波数の精度を高めることができる。
【発明を実施するための最良の形態】
【００１２】
以下、図面を参照し、この発明の実施の形態を説明する。
【００１３】
＜全体構成＞
図１は、この発明の一実施形態による音分析プログラムの処理内容を示す図である。この音分析プログラムは、自然界から音響信号を取得する収音機能、ＣＤ等の記録媒体から音楽の音響信号を再生する再生機能またはネットワークを介して音楽の音響信号を取得する通信機能等の音響信号取得機能を備えたパーソナルコンピュータ等のコンピュータにインストールされて実行される。本実施形態による音分析プログラムを実行するコンピュータは、本実施形態による音分析装置として機能する。
【００１４】
本実施形態による音分析プログラムは、音響信号取得機能を介して取得されたモノラルの音楽音響信号に対し、その中のある音源の音高を推定する。その最も重要な例として、ここではメロディラインとベースラインを推定する。メロディは他よりも際立って聞こえる単音の系列、ベースはアンサンブル中で最も低い単音の系列であり、その時間的な変化の軌跡をそれぞれメロディラインＤｍ（ｔ）、ベースラインＤｂ（ｔ）と呼ぶ。時刻tにおける基本周波数Ｆ０をＦｉ（ｔ）（ｉ＝ｍ，ｂ）、振幅をＡｉ（ｔ）とすると、これらは以下のように表される。
【数１】
2008209579.tif 000003
【数２】
2008209579.tif 000004
【００１５】
このメロディラインＤｍ（ｔ）およびベースラインＤｂ（ｔ）を入力音響信号から得るための手段として、音分析プログラムは、瞬時周波数の算出１、アタック検出１ａ、周波数成分の候補の抽出２、周波数帯域の制限３、メロディラインの推定４ａおよびベースラインの推定４ｂの各処理を含む。また、メロディラインの推定４ａおよびベースラインの推定４ｂの各処理は、基本周波数の確率密度関数の推定４１およびマルチエージェントモデルによる基本周波数の継時的な追跡４２の各処理を各々含む。本実施形態において、瞬時周波数の算出１、周波数成分の候補の抽出２、周波数帯域の制限３の処理内容は、前掲特許文献１に開示されたものと基本的に同様である。本実施形態の特徴は、アタック検出１ａを設けた点並びにこのアタック検出１ａの処理結果に基づいて制御されるメロディラインの推定４ａおよびベースラインの推定４ｂの処理内容にある。以下、本実施形態による音分析プログラムを構成する各処理の内容を説明する。
【００１６】
＜瞬時周波数の算出１＞
この瞬時周波数の算出１と、アタック検出１ａと、周波数成分の候補の抽出２と、周波数帯域の制限３と、メロディラインの推定４ａおよびベースラインの推定４ｂにおける基本周波数の確率密度関数の推定４１の各処理は、時間軸上において音響信号を分割した一定時間長のフレームを単位として実行される。以下において時刻ｔは、具体的にはフレームの番号である。瞬時周波数の算出１では、入力音響信号を複数のＢＰＦからなるフィルタバンクに与え、フィルタバンクの各ＢＰＦの出力信号について、位相の時間微分である瞬時周波数（Flanagan, J.L. and Golden, R.M.: Phase Vocoder, The BellSystem
Technical J., Vol.45, pp.1493-1509 (1966)参照）を計算する。ここでは、上記Flanaganの手法を用い、短時間フーリエ変換(STFT)の出力をフィルタバンク出力と解釈して、効率良く瞬時周波数を計算する。入力音響信号ｘ(ｔ)に対する窓関数ｈ(ｔ)を用いたＳＴＦＴが式（３）および（４）により与えられるとき、瞬時周波数λ（ω，ｔ）は式（５）により求めることができる。
【００１７】
【数３】
2008209579.tif 000005
【数４】
2008209579.tif 000006
【数５】
2008209579.tif 000007
【００１８】
ここで、ｈ(ｔ)は時間周波数の局所化を与える窓関数である（例えば、最適な時間周波数の局所化を与えるガウス関数に２階のカーディナルＢ−スプライン関数を畳み込んで作成した時間窓など)。
【００１９】
この瞬時周波数を計算するのに、ウェーブレット変換を用いても良い。ここでは、計算量を減らすためにＳＴＦＴを用いるが、単一のＳＴＦＴのみを用いたのでは、ある周波数帯域における時間分解能や周波数分解能が悪くなってしまう。そこで、マルチレートフィルタバンク（Vetterli, M.: A Theory of Multirate Filter Banks, IEEE Trans. on ASSP,
Vol.ASSP-35, No.3, pp. 356-372 (1987)、参照）を構成し、リアルタイムに実行可能という制約のもとで、ある程度妥当な時間周波数分解能を得る。
【００２０】
＜アタック検出１ａ＞
この処理では、時間軸上において入力音響信号を分割した各フレームが入力音響信号のアタック区間内のフレームであるか否かの判定を行い、フレーム毎にそのフレームがアタック区間のものか否かを示す情報をメロディラインの推定４ａおよびベースラインの推定４ｂに引き渡す。フレームが入力音響信号のアタック区間のものであるか否かの判定方法には周知の各種の方法があるが、例えば特許文献２に開示されているように、各フレームをより時間長の短い複数の解析区間に分割し、これらの複数の解析区間における音響信号のエネルギーの変動を解析することにより、当該フレームがアタック区間のものか否かを判定することが可能である。
【００２１】
＜周波数成分の候補の抽出２＞
この処理では、フィルタの中心周波数からその瞬時周波数への写像に基づいて、周波数成分の候補を抽出する（Charpentier, F.J.: Pitch detection using the short-termphase
spectrum, Proc. of ICASSP 86, pp.113-116 (1986)参照）。あるＳＴＦＴフィルタの中心周波数ωからその出力の瞬時周波数λ（ω，ｔ）への写像を考える。すると、もし周波数ψの周波数成分があるときには、ψがこの写像の不動点に位置し、その周辺の瞬時周波数の値はほぼ一定となる。つまり、全周波数成分の瞬時周波数Ψf(t)は、次式によって抽出することができる。
【数６】
2008209579.tif 000008
これらの周波数成分のパワーは、Ψf(t)の各周波数におけるＳＴＦＴパワースペクトルの値として得られるため、周波数成分のパワー分布関数Ψp(t)(ω)を次のように定義することができる。
【数７】
2008209579.tif 000009
【００２２】
＜周波数帯域の制限３＞
この処理では、抽出した周波数成分に重み付けすることで、周波数帯域を制限する。ここでは、メロディラインとベースライン用に、２種類のＢＰＦを用意する。メロディライン用ＢＰＦは、典型的なメロディラインの主要な基本波成分および高調波成分の多くを通過させることができ、かつ、基本周波数付近の重複が頻繁に起きる周波数帯域をある程度遮断する。一方、ベースライン用ＢＰＦは、典型的なベースラインの主要な基本周波数成分および高調波成分の多くを通過させることができ、かつ、他の演奏パートがベースラインよりも優勢になるような周波数帯域をある程度遮断する。
【００２３】
本実施形態では、以下、対数スケールの周波数をcentの単位(本来は音高差(音程)を表す尺度)で表し、Ｈｚで表された周波数ｆＨｚを、次のようにcentで表された周波数ｆｃｅｎｔに変換する。 
【数８】
2008209579.tif 000010
【数９】
2008209579.tif 000011
平均律の半音は１００ｃｅｎｔに、１オクターブは１２００ｃｅｎｔに相当する。
【００２４】
周波数ｘ ｃｅｎｔでのＢＰＦの周波数応答をＢＰＦｉ（ｘ）（ｉ＝ｍ，ｂ）とし、周波数成分のパワー分布関数をΨ’ｐ（ｔ）（ｘ）とすると、ＢＰＦを通過した周波数成分はＢＰＦｉ（ｘ）Ψ’ｐ（ｔ）（ｘ）と表すことができる。ただし、Ψ’ｐ（ｔ）（ｘ）は、周波数軸がｃｅｎｔで表されていることを除けばΨｐ（ｔ）（ω）と同じ関数である。ここで、次の段階の準備として、ＢＰＦを通過した周波数成分の確率密度関数ｐΨ（ｔ）（ｘ）を定義する。
【数１０】
2008209579.tif 000012
【００２５】
ここで、Ｐｏｗ（ｔ）は次式に示すようにＢＰＦを通過した周波数成分のパワーの合計である。
【数１１】
2008209579.tif 000013
【００２６】
＜基本周波数の確率密度関数の推定４１＞
この処理では、ＢＰＦを通過した周波数成分の候補に対し、各高調波構造が相対的にどれくらい優勢かを表す基本周波数の確率密度関数を求める。そのために、本実施形態では、周波数成分の確率密度関数ｐΨ（ｔ）（ｘ）が、高調波構造を持つ音をモデル化した確率分布（音モデル）の混合分布モデル（重み付き和のモデル）から生成されたと考える。基本周波数がＦの音モデルの確率密度関数をｐ（ｘ｜Ｆ）とすると、その混合分布モデルｐ（ｘ；θ（ｔ））は、次式により定義することができる。
【数１２】
2008209579.tif 000014
【数１３】
2008209579.tif 000015
【００２７】
ここで、ＦｈｉとＦｌｉは、許容される基本周波数の上限と下限であり、ＢＰＦの通過帯域により決定される。また、ｗ（ｔ）（Ｆ）は、次式を満たすような、音モデルｐ（ｘ｜Ｆ）の重みである。
【数１４】
2008209579.tif 000016
【００２８】
ＣＤ等による実世界の音響信号に対して事前に音源数を仮定することは不可能なため、このように、あらゆる基本周波数の可能性を同時に考慮してモデル化することが重要となる。もし、観測した周波数成分ｐΨ（ｔ）（ｘ）がモデルｐ（ｘ；θ（ｔ））から生成されたかのようにモデルパラメータθ（ｔ）を推定できれば、ｐΨ（ｔ）（ｘ）は個々の音モデルへと分解されたとみなすことができ、次式に示すように、各基本周波数Ｆの音モデルに対する重みｗ（ｔ）（Ｆ）を、基本周波数Ｆの確率密度関数ｐＦＯ（ｔ）（Ｆ）と解釈することができる。
【数１５】
2008209579.tif 000017
【００２９】
つまり、混合分布中において、ある音モデルｐ（ｘ｜Ｆ）が優勢になればなるほど（すなわち、ｗ（ｔ）（Ｆ）が大きくなるほど）、ｐＦＯ（ｔ）（Ｆ）において、そのモデルの基本周波数Ｆの確率が高くなる。
【００３０】
以上から、確率密度関数ｐΨ（ｔ）（ｘ）を観測したときに、そのモデルｐ（ｘ；θ（ｔ））のパラメータθ（ｔ）を推定する問題を解けばよいことがわかる。θ（ｔ）の最尤推定量は、次式で定義される平均対数尤度を最大化することで得られる。
【数１６】
2008209579.tif 000018
【００３１】
この最大化問題は解析的に解くことが困難なため、前述のＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）アルゴリズムを用いてθ（ｔ）を推定する。ＥＭアルゴリズムは、Ｅステップ（ｅｘｐｅｃｔａｔｉｏｎ  ｓｔｅｐ）とＭステップ（ｍａｘｉｍｉｚａｔｉｏｎ  ｓｔｅｐ）を交互に繰返し適用することで、不完全な観測データ（この場合、ｐΨ（ｔ）（ｘ））から最尤推定をおこなうための反復アルゴリズムである。本実施形態では、ＥＭアルゴリズムを繰り返すことにより、ＢＰＦを通過した周波数成分の確率密度関数ｐΨ（ｔ）（ｘ）を、各種の基本周波数Ｆに対応した複数の音モデルｐ（ｘ｜Ｆ）を重み付け加算した混合分布と考える場合において、最も尤もらしい重みのパラメータθ（ｔ）（＝｛ｗ（ｔ）（Ｆ）｜Ｆｌｉ≦Ｆ≦Ｆｈｉ｝を求める。ここで、ＥＭアルゴリズムの各繰り返しでは、パラメータθ（ｔ）（＝｛ｗ（ｔ）（Ｆ）｜Ｆｌｉ≦Ｆ≦Ｆｈｉ｝）に関して、古いパラメータ推定値θold（ｔ）（＝｛ｗold（ｔ）（Ｆ）｜Ｆｌｉ≦Ｆ≦Ｆｈｉ｝）を更新して新しい（より尤もらしい）パラメータ推定値θｎｅｗ（ｔ）（＝｛ｗｎｅｗ（ｔ）（Ｆ）｜Ｆｌｉ≦Ｆ≦Ｆｈｉ｝）を求めていく。この古いパラメータ推定値θold（ｔ）から新しいパラメータ推定値θｎｅｗ（ｔ）を求める漸化式は、次のようになる。なお、この漸化式の導出過程は特許文献１に詳細に説明されているので、そちらを参照されたい。
【数１７】
2008209579.tif 000019
【００３２】
図２は、本実施形態において音モデルｐ（ｘ｜Ｆ）に対する重みのパラメータθ（ｔ）（＝｛ｗ（ｔ）（Ｆ）｜Ｆｌｉ≦Ｆ≦Ｆｈｉ｝がＥＭアルゴリズムにより更新される過程を例示したものである。なお、図２では、図示を簡素化するために、周波数成分の数が４個である音モデルを使用した場合の例が示されている。
【００３３】
本実施形態におけるＥＭアリゴリズムでは、各基本周波数Ｆに対応した音モデルｐ（ｘ｜Ｆ）と現状における各音モデルに対する重み値ｗold（ｔ）（Ｆ）とに基づき、次式に従い、周波数ｘ毎に各音モデルに対応したスペクトル分配比を求める。
【数１８】
2008209579.tif 000020
【００３４】
上記式（１８）に示すように、ある周波数ｘにおける各音モデルｐ（ｘ｜Ｆ）に対応したスペクトル分配比（ｘ｜Ｆ）は、重み値ｗold（Ｆ）（ｔ）の乗算された各音モデルｐ（ｘ｜Ｆ）の周波数ｘにおける各振幅値ｗold（Ｆ）（ｔ）ｐ（ｘ｜Ｆ）の総和（式（１８）における分母の積分値に相当）を求め、その総和により各振幅値ｗold（Ｆ）（ｔ）ｐ（ｘ｜Ｆ）を除算することにより得られる。式（１８）から明らかなように、各周波数ｘにおいて、各音モデルｐ（ｘ｜Ｆ）に対応した各スペクトル分配比（ｘ｜Ｆ）は、総和が１となるように正規化されたものとなる。
【００３５】
そして、本実施形態では、各周波数ｘにおいて、その周波数ｘにおける確率密度関数ｐΨ（ｔ）（ｘ）の関数値をその周波数ｘにおける各音モデルｐ（ｘ｜Ｆ）のスペクトル分配比により分配し、音モデルｐ（ｘ｜Ｆ）毎に、このようにして分配された確率密度関数ｐΨ（ｔ）（ｘ）の関数値を総計して、各音モデルｐ（ｘ｜Ｆ）の取り分とする。そして、全音モデルの取り分を合計し、その合計値によって各音モデルの取り分を除算し、総和が１となるように正規化された各音モデルｐ（ｘ｜Ｆ）の取り分を新たな重みパラメータｗｎｅｗ（ｔ）（Ｆ）とするのである。以上の処理が繰り返されることにより、基本周波数Ｆの異なった各音モデルｐ（ｘ｜Ｆ）のうち混合音の周波数成分の確率密度関数ｐΨ（ｔ）（ｘ）によって支持されている確率の高いものに対する重みパラメータｗ（ｔ）（Ｆ）が次第に強調されてゆく。この結果、重みパラメータｗ（ｔ）（Ｆ）は、ＢＰＦを通過した混合音における基本周波数の確率密度関数を示すに至る。
【００３６】
＜マルチエージェントモデルによる基本周波数の継時的な追跡４２（基本周波数推定手段としての処理）＞
本実施形態による音分析プログラムは、以上のようにして得られる基本周波数の確率密度関数に基づいて、入力音響信号に含まれる１または複数の音源の音の基本周波数を推定して出力する基本周波数推定手段としての処理を含む。この処理では、最も優勢な基本周波数Ｆｉ（ｔ）を決定するために、次式に示すように、基本周波数の確率密度関数ｐＦ０（ｔ）（Ｆ）（式（１５）より、式（１７）を反復計算した最終的な推定値として得られる）を最大にする周波数を基本周波数の推定値として求める。
【数１９】
2008209579.tif 000021
【００３７】
ところで、基本周波数の確率密度関数において、同時に鳴っている音の基本周波数に対応する複数のピークが拮抗すると、それらのピークが確率密度関数の最大値として次々に選ばれてしまうことがあるため、このように単純に求めた結果は安定しないことがある。そこで、本実施形態における基本周波数推定手段としての処理では、大局的な観点から基本周波数を推定するために、基本周波数の確率密度関数の時間変化において複数のピークの軌跡を継時的に追跡し、その中で最も優勢で安定した基本周波数の軌跡を選択する。このような追跡処理を動的で柔軟に制御するために、マルチエージェントモデルを導入する。
【００３８】
マルチエージェントモデルは、１つの特徴検出器と複数のエージェントにより構成される（図３参照）。特徴検出器は、基本周波数の確率密度関数の中で目立つピークを拾い上げる。エージェントは基本的に、それらのピークに駆動されて軌跡を追跡していく。つまり、マルチエージェントモデルは、入力中で目立つ特徴を時間的に追跡する汎用の枠組みである。具体的には、各時刻において以下の処理がおこなわれる。
【００３９】
（１）基本周波数の確率密度関数が求まった後、特徴検出器は目立つピーク（最大ピークに応じて動的に変化する閾値を越えたピーク）を複数検出する。そして、目立つピークのそれぞれについて、周波数成分のパワーの合計Ｐｏｗ（ｔ）も考慮しながら、どれくらい将来有望なピークかを評価する。これは、現在時刻を数フレーム先の時刻とみなして、ピークの軌跡をその時刻まで先読みして追跡することで実現する。
【００４０】
（２）既に生成されたエージェントがあるときは、それらが相互作用しながら、目立つピークをそれに近い軌跡を持つエージェントへと排他的に割り当てる。複数のエージェントが割り当て候補に上がる場合には、最も信頼度の高いエージェントへと割り当てる。
【００４１】
（３）最も有望で目立つピークがまだ割り当てられていないときは、そのピークを追跡する新たなエージェントを生成する。
【００４２】
（４）各エージェントは累積ペナルティを持っており、それが一定の閾値を越えると消滅する。
【００４３】
（５）目立つピークが割り当てられなかったエージェントは、一定のペナルティを受け、基本周波数の確率密度関数の中から自分の追跡する次のピークを直接見つけようとする。もしそのピークも見つからないときは、さらにペナルティを受ける。さもなければ、ペナルティはリセットされる。
【００４４】
（６）各エージェントは、今割り当てられたピークがどれくらい有望で目立つかを表す度合いと、１つ前の時刻の信頼度との重み付き和によって、信頼度を自己評価する。
【００４５】
（７）時刻ｔにおける基本周波数Ｆｉ（ｔ）は、信頼度が高く、追跡しているピークの軌跡に沿ったパワーの合計が大きいエージェントに基づいて決定する。振幅Ａｉ（ｔ）は、基本周波数Ｆｉ（ｔ）の高調波成分等をΨｐ（ｔ）（ω）から抽出して決定する。
【００４６】
＜＜特許文献１の技術に対する本実施形態の改良点＞＞
図４は本実施形態における基本周波数の確率密度関数の推定４１の処理内容を示している。図４に示すように、基本周波数の確率密度関数の推定４１においては、ＥＭアルゴリズムのＥステップおよびＭステップ４１１と、収束判定４１２とを繰り返す。
【００４７】
まず、ＥステップおよびＭステップ４１１では、前掲式（１７）の漸化式に従い、基本周波数の確率密度関数、すなわち、各種の基本周波数Ｆに対応した音モデルの重み値θ＝θｎｅｗ（ｔ）（＝｛ｗnew（ｔ）（Ｆ）｜Ｆｌｉ≦Ｆ≦Ｆｈｉ｝）を求める。
【００４８】
次に収束判定４１２では、今回のＥステップおよびＭステップ４１１において得られた各種の基本周波数Ｆに対応した音モデルの重み値θ＝θｎｅｗ（ｔ）とその前の重み値θ＝θold（ｔ）とを比較し、重み値θの変化分が許容範囲内に収まったか否かを判定する。そして、重み値θの変化分が許容範囲内に収まったと判定した場合には、基本周波数の確率密度関数の推定４１の処理を終了し、基本周波数の確率密度関数の最終値をマルチエージェントモデルによる基本周波数の継時的な追跡４２に引き渡す。
【００４９】
本実施形態による音分析プログラムには、アタック検出１ａから出力される情報に基づいて、基本周波数の確率密度関数の推定４１または基本周波数推定手段であるマルチエージェントモデルによる基本周波数の継時的な追跡４２の処理の演算態様を制御する演算制御手段が設けられている。これが特許文献１の技術に対する本実施形態の改良点である。基本周波数の確率密度関数の推定４１または基本周波数推定手段としてのマルチエージェントモデルによる基本周波数の継時的な追跡４２の処理の演算態様の制御には、以下の４態様がある。ユーザは、図示しない操作部の操作により、音分析プログラムの演算制御手段にいずれの態様で演算態様の制御を行わせるかを指定することができる。
【００５０】
＜＜＜第１の態様＞＞＞
この第１の態様は、処理対象であるフレームがアタック区間のものである場合には、当該フレームにおける重み値ｗ（ｔ）（Ｆ）の逐次更新が所定の初期値ｗflat（Ｆ）から開始されるように、基本周波数の確率密度関数の推定４１の演算制御を行い、処理対象であるフレームがアタック区間のものでない場合には、当該フレームにおける重み値ｗ（ｔ）（Ｆ）の逐次更新が前フレームにおける重み値ｗ（ｔ-1）（Ｆ）の最終値を初期値として開始されるように基本周波数の確率密度関数の推定４１のための演算制御を行う態様である。
【００５１】
従来の技術の下では、フレーム毎に、上述した漸化式（１７）を繰り返して基本周波数の確率密度関数の推定を行う場合に、ｗ(ｔ)(Ｆ)の初期値として、１つ前の時刻ｔ−１（１つ前のフレーム）における重み値の最終値ｗ(t-1)(Ｆ)を用いた。しかし、このように前フレームにおける基本周波数の確率密度関数の最終状態を初期値として用いると、波形が不安定なアタック区間のフレームの基本周波数の推定を行う際に、推定の処理が不安定になり、誤推定に陥り易い。そこで、この第１の態様では、各フレームに関して重み値ｗ（ｔ）（Ｆ）の逐次更新を開始する際、そのフレームがアタック区間以外のものである場合には前フレームにおける重み値の最終値ｗ(t-1)(Ｆ)を初期値とし、アタック区間のものである場合には、例えば全周波数帯域においてフラットな重み値を持った所定の初期値ｗflat（Ｆ）を初期値とするのである。より具体的には、次の通りである。
【００５２】
まず、本実施形態では、各フレームでの重み値ｗ(t)(Ｆ)の逐次更新を開始するとき、図４に示すように、前フレームでの重み値の最終値ｗ(t-1)(Ｆ)に係数ｒを乗算したものと、所定の初期値ｗflat（Ｆ）に係数１−ｒを乗算したものとを加算し、その加算結果を当該フレームにおける重み値ｗ(t)(Ｆ)の初期値とする。
【００５３】
そして、図５に示すように、アタック区間に属しないフレームの処理時には、ｒの値を１とすることにより、前フレームにおける重み値の最終値ｗ(t-1)(Ｆ)を重み値ｗ(t)(Ｆ)の初期値とし、アタック区間に属するフレームの処理時には、ｒの値を０とすることにより、所定の初期値ｗflat（Ｆ）を重み値ｗ(t)(Ｆ)の初期値とするのである。
【００５４】
以上のように、この態様によれば、アタック区間のフレームの処理時には、当該フレームでの重み値ｗ(t)(Ｆ)の逐次更新の際に、前フレームにおける最終的な重み値ｗ(t-1)(Ｆ)が初期値として採用されない。従って、アタック区間に入力音響信号の波形が不安定となり、基本周波数の推定処理が不安定になる場合でも、連続して基本周波数の誤推定が発生するのを回避することができ、全体としての基本周波数の推定精度を高めることができる。
【００５５】
＜＜＜第２の態様＞＞＞
強いタッチで楽器演奏が行われた等の場合には、音響信号のアタック区間が終了しても暫くの間は波形の不安定な状態が続くことがある。そのような場合、アタック区間が終了した後のフレームであっても、前フレームの最終的な重み値ｗ(t-1)(Ｆ)を初期値として使用して、重み値ｗ(t)(Ｆ)の逐次更新を行うと、誤った基本周波数において重み値がピークとなり、基本周波数の誤推定が発生するおそれがある。
【００５６】
そこで、この第２の態様では、次のようにして、重み値ｗ(ｔ)(Ｆ)の初期値を制御する。まず、第２の態様でも、上記第１の態様と同様、各フレームでの重み値ｗ(t)(Ｆ)の逐次更新を開始するとき、図４に示すように、前フレームでの重み値の最終値ｗ(t-1)(Ｆ)に係数ｒを乗算したものと、所定の初期値ｗflat（Ｆ）に係数１−ｒを乗算したものとを加算し、その加算結果を当該フレームにおける重み値ｗ(t)(Ｆ)の初期値とする。
【００５７】
そして、図６に示すように、アタック区間に属するフレームの処理時には、ｒの値を０とすることにより、所定の初期値ｗflat（Ｆ）を重み値ｗ(t)(Ｆ)の初期値とする。
【００５８】
そして、アタック区間が終了した後は、フレームが切り換わる都度、ｒの値を０から１に向けて徐々に高めてゆく。すなわち、第２の態様では、アタック区間に属しない各フレームの処理時には、当該フレームにおける重み値ｗ(t)(Ｆ)の逐次更新が前フレームにおける重み値の最終値ｗ(t-1)(Ｆ)と所定の初期値ｗflat（Ｆ）とをミキシングした重み値を初期値として開始されるように基本周波数の確率密度関数の推定４１の演算制御を行い、かつ、直前のアタック区間の終了時からの経過時間が長くなるに従って、前フレームでの重み値の最終値ｗ(t-1)(Ｆ)が強調されるように、前フレームにおける重み値の最終値ｗ(t-1)(Ｆ)と所定の初期値ｗflat（Ｆ）とのミキシング比を制御するのである。
【００５９】
この態様によれば、アタック区間が終了しても暫くの間は音響信号の波形が不安定になる状況でも、基本周波数の誤推定が連続して発生するのを回避し、基本周波数の推定の精度を高めることができる。
【００６０】
＜＜＜第３の態様＞＞＞
上記第１の態様および第２の態様では、アタック検出１ａから引き渡される情報に従い、音モデルに対する重み値ｗ(t)(Ｆ)の初期値の制御を行った。これに対し、第３の態様では、図７に示すように、アタック区間以外の区間では、通常の音モデルが用いられ、アタック区間では、アタック区間用の音モデルが用いられるように、ＥステップおよびＭステップ４１１に用いられる音モデルの切り換えを行う。
【００６１】
ここで、アタック区間用の音モデルとしては、実際の楽音の高調波構造よりも顕著なピークが少なく、周波数軸上において各高調波成分の振幅値が緩やかなカーブを描いて変化する高調波構造を持った音モデルを用いる。アタック区間においてこのような音モデルを用いることにより、入力音響信号の波形の変化に対して安定性の良い基本周波数の推定を行うことが可能となる。
【００６２】
＜＜＜第４の態様＞＞＞
上記第１〜第３の態様では、基本周波数の確率密度関数の推定４１が、アタック検出１ａの処理結果に基づく制御の対象となった。これに対し、この第４の態様では、基本周波数推定手段であるマルチエージェントモデルによる基本周波数の継時的な追跡４２が、アタック検出１ａの処理結果に基づく制御の対象となる。すなわち、この第４の態様において音分析プログラムの演算制御手段は、アタック区間では、基本周波数の確率密度関数の推定４１により得られる基本周波数の確率密度関数が得られたとしても、その確率密度関数に基づく基本周波数の推定および出力を行わないように、マルチエージェントモデルによる基本周波数の継時的な追跡４２を制御する。すなわち、基本周波数の誤推定が発生するアタック区間では、基本周波数の推定および出力を行わず、出力される基本周波数についてのみ推定の精度を高める趣旨である。
【００６３】
＜他の実施形態＞
以上、この発明の一実施形態について説明したが、この発明には他にも実施形態があり得る。例えば次の通りである。
（１）第１の態様または第２の態様の一方と、第３の態様とを併用し得るように音分析プログラムの演算制御手段を構成してもよい。
（２）第３の態様において、アタック区間のみならず、アタック区間が終了してから所定時間が経過するまでの期間についても、アタック区間用の音モデルを使用して、基本周波数の確率密度関数の推定４１を実行するように制御してもよい。
（３）第４の態様において、アタック区間のみならず、アタック区間が終了してから所定時間が経過するまでの期間についても、基本周波数の推定および出力を停止させてもよい。
【図面の簡単な説明】
【００６４】
【図１】この発明の一実施形態である音分析プログラムの処理内容を示す図である。
【図２】同実施形態において音モデルに対する重みのパラメータがＥＭアルゴリズムにより更新される過程を例示した図である。
【図３】１つの特徴検出器と複数のエージェントにより構成されるマルチエージェントモデルによる基本周波数の経時的な追跡を示す図である。
【図４】同実施形態における基本周波数の確率密度関数の推定４１の処理内容を示す図である。
【図５】同実施形態による音分析プログラムの演算制御手段により実行される演算制御の第１の態様を示すタイムチャートである。
【図６】同演算制御手段により実行される演算制御の第２の態様を示すタイムチャートである。
【図７】同演算制御手段により実行される演算制御の第３の態様を示すタイムチャートである。
【符号の説明】
【００６５】
１……瞬時周波数の算出、１ａ……アタック検出、２……周波数成分の候補の抽出、３……周波数帯域の制限、４ａ……メロディラインの推定、４ｂ……ベースラインの推定、４１……基本周波数の確率密度関数の推定、４２……マルチエージェントモデルによる基本周波数の継時的な追跡。
【図１】
2008209579.tif 000022
【図２】
2008209579.tif 000023
【図３】
2008209579.tif 000024
【図４】
2008209579.tif 000025
【図５】
2008209579.tif 000026
【図６】
2008209579.tif 000027
【図７】
2008209579.tif 000028