(57)【要約】【課題】  記憶する音モデルの個数が比較的少なくて済み、かつ、高い推定精度で基本周波数を推定できるようにする。【解決手段】  音モデル補間処理５では、記憶装置に記憶された複数種類の音モデルに補間処理を施し、各音モデルの中間の基本周波数に対応した複数種類の音モデルを生成する。基本周波数の確率密度関数の推定４１では、記憶された音モデルおよび補間処理により得られた音モデルを用いて、様々な基本周波数の音モデルの混合分布を構成し、この混合分布が入力音響信号の周波数成分の分布となるように、各音モデルに対する重み値を最適化し、最適化された各音モデルの重み値を音源の音の基本周波数の確率密度関数とする。【選択図】図１
2008058753.tif 000002
【特許請求の範囲】
【請求項１】
楽器から発音される複数種類の音の高調波構造を各々定義した複数種類の音モデルを記憶する記憶手段と、
前記記憶手段に記憶された複数種類の音モデルを各々の基本周波数に従って序列化し、序列化された複数種類の音モデルに対して、基本周波数に基づく補間処理を施し、序列化された各音モデルの中間の基本周波数に対応した複数種類の音モデルを生成する補間手段と、
前記記憶手段に記憶された複数種類の音モデルおよび前記補間手段により生成された複数種類の音モデルを用いて、各種の高調波構造および基本周波数を有する複数の音モデルを重み付け加算した混合分布を構成し、この混合分布が入力音響信号の周波数成分の分布となるように、各音モデルに対する重み値を最適化し、最適化された各音モデルの重み値を前記入力音響信号の元である音源の音の基本周波数の確率密度関数として推定する確率密度関数推定手段と、
前記基本周波数の確率密度関数に基づいて前記入力音響信号における１または複数の音源の音の基本周波数を推定して出力する基本周波数推定手段と
を具備することを特徴とする音分析装置。
【請求項２】
前記確率密度推定手段は、各音モデルの種類毎に、当該音モデルが本来有していた基本周波数に合わせて、基本周波数の範囲を設定し、この設定した基本周波数の範囲外の周波数における当該音モデルに対する重み値を制限して、前記各音モデルに対する重み値の最適化を行うことを特徴とする請求項１に記載の音分析装置。
【請求項３】
前記基本周波数推定手段は、前記確率密度関数におけるピークを複数検出し、各ピークの信頼度に基づき、信頼度が高くパワーの大きい基本周波数を出力することを特徴とする請求項１に記載の音分析装置。
【請求項４】
コンピュータを、
楽器から発音される複数種類の音の高調波構造を各々定義した複数種類の音モデルを記憶する記憶手段と、
前記記憶手段に記憶された複数種類の音モデルを各々の基本周波数に従って序列化し、序列化された複数種類の音モデルに対して、基本周波数に基づく補間処理を施し、序列化された各音モデルの中間の基本周波数に対応した複数種類の音モデルを生成する補間手段と、
前記記憶手段に記憶された複数種類の音モデルおよび前記補間手段により生成された複数種類の音モデルを用いて、各種の高調波構造および基本周波数を有する複数の音モデルを重み付け加算した混合分布を構成し、この混合分布が入力音響信号の周波数成分の分布となるように、各音モデルに対する重み値を最適化し、最適化された各音モデルの重み値を前記入力音響信号の元である音源の音の基本周波数の確率密度関数として推定する確率密度関数推定手段と、
前記基本周波数の確率密度関数に基づいて前記入力音響信号における１または複数の音源の音の基本周波数を推定して出力する基本周波数推定手段と
して機能させることを特徴とするコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、市販のＣＤ（ｃｏｍｐａｃｔ  ｄｉｓｃ）などに収録されている、歌声や複数種類の楽器音を同時に含む音楽音響信号を対象に、メロディ音やベース音の音高（本明細書では基本周波数の意味で用いる）を推定する音分析装置およびプログラムに関する。
【背景技術】
【０００２】
多数の音源の音が混ざり合ったモノラルの音響信号中から、ある特定の音源の音高を推定することは、非常に困難である。混合音に対して音高推定することが難しい本質的な理由の１つに、時間周波数領域において、ある音の周波数成分が同時に鳴っている他の音の周波数成分と重複することが挙げられる。例えば、歌声、鍵盤楽器（ピアノ等）、ギター、ベースギター、ドラムス等で演奏される典型的なポピュラー音楽では、メロディを担う歌声の高調波構造の一部（特に基本周波数成分）は、鍵盤楽器、ギターの高調波成分やベースギターの高次の高調波成分、スネアドラム等の音に含まれるノイズ成分などと頻繁に重複する。そのため、各周波数成分を局所的に追跡するような手法は、複雑な混合音に対しては安定して機能しない。基本周波数成分が存在することを前提に高調波構造を推定する手法もあるが、そのような手法は、ミッシングファンダメンタル（ｍｉｓｓｉｎｇ  ｆｕｎｄａｍｅｎｔａｌ）現象を扱えないという大きな欠点を持つ。さらに、同時に鳴っている他の音の周波数成分が基本周波数成分と重複すると、有効に機能しない。
【０００３】
以上のような理由により、従来、単一音のみか、非周期的な雑音を伴った単一音を収録した音響信号を対象とした音高の推定技術はあったが、市販のＣＤに記録された音響信号のように複数の音が混ざり合ったものについて音高を推定する技術はなかった。
【０００４】
しかしながら、近年、統計的手法を利用することにより、混合音に含まれる各音の音高を適切に推定する技術が提案されるに至った。特許文献１の技術である。
【０００５】
この特許文献１の技術では、メロディ音のものと考えられる帯域に属する周波数成分と、ベース音のものと考えられる帯域に属する周波数成分とを入力音響信号からＢＰＦにより別々に取り出し、それらの各帯域の周波数成分に基づき、メロディ音およびベース音の各々の基本周波数の推定を行う。
【０００６】
さらに詳述すると、特許文献１の技術では、音の高調波構造に対応した確率分布を持った音モデルを用意し、メロディ音の帯域の各周波数成分、ベース音の帯域の各周波数成分が、様々な基本周波数に対応した各音モデルを重み付け加算した混合分布であると考える。そして、各音モデルの重みの値をＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）アルゴリズムを用いて推定する。
【０００７】
このＥＭアルゴリズムは、隠れ変数を含む確率モデルに対して最尤推定を行うための反復アルゴリズムであり、局所最適解を求めることができる。ここで、最も大きな重みの値を持つ確率分布は、その時点で最も優勢な高調波構造であるとみなすことができるため、あとはその優勢な高調波構造における基本周波数を音高として求めればよい。この手法は基本周波数成分の存在に依存しないため、ミッシングファンダメンタル現象も適切に扱うことができ、基本周波数成分の存在に依存せずに、最も優勢な高調波構造を求めることができる。
【０００８】
非特許文献１は、特許文献１の技術に対して次の拡張を行った技術を開示している。
【０００９】
＜拡張１：音モデルの多重化＞
特許文献１の技術では、同一基本周波数には１つの音モデルしか用意されていなかったが、実際には、ある基本周波数に、異なる高調波構造を持つ音が入れ替わり立ち替わり現れることがある。そこで、同一基本周波数に対して複数の音モデルを用意し、入力音響信号をそれらの混合分布としてモデル化した。
【００１０】
＜拡張２：音モデルのパラメータの推定＞
特許文献１の技術では、音モデルにおいて各高調波成分の大きさの比を固定していた（ある理想的な音モデルを仮定していた）。これは実世界の混合音中の高調波構造とは必ずしも一致しておらず、精度向上のためには洗練される余地が残されていた。そこで、音モデルの高調波成分の比率もモデルパラメータに加え、各時刻においてＥＭアルゴリズムにより推定するようにした。
【００１１】
＜拡張３：モデルパラメータに関する事前分布の導入＞
特許文献１の技術では、音モデルの重み（基本周波数の確率密度関数）に関する事前知識は仮定していなかった。しかし、この基本周波数の推定技術の用途によっては、たとえ事前に基本周波数がどの周波数の近傍にあるかを与えてでも、より誤検出の少ない基本周波数を求めたいというような要求も発生し得る。例えば、演奏分析やビブラート分析等の目的では、楽曲をヘッドホン聴取しながらの歌唱や楽器演奏によって、各時刻におけるおおよその基本周波数を事前知識として用意しておき、実際の楽曲中のより正確な基本周波数を得ることが求められている。そこで、特許文献１におけるモデルパラメータ（音モデルに対する重み値）の最尤推定の枠組みを拡張し、モデルパラメータに関する事前分布に基づいて最大事後確率推定（ＭＡＰ推定；Maximum A Posteriori Probability Estimation）を行うようにした。その際、＜拡張２＞においてモデルパラメータに加えた音モデルの高調波成分の大きさの比率に関する事前分布も導入した。
【特許文献１】特許第３４１３６３４号
【非特許文献１】後藤 真孝:"リアルタイム音楽情景記述システム: 全体構想と音高推定手法の拡張", 情報処理学会 音楽情報科学研究会 研究報告 2000-MUS-37-2,Vol.2000, No.94, pp.9-16, ２０００年１０月１６日
【発明の開示】
【発明が解決しようとする課題】
【００１２】
非特許文献１に開示された技術によれば、上記「拡張１」を導入したことにより、例えば音源が高調波構造の異なった複数の音を発生し得るような場合に、それらの各高調波構造に対応した複数の音モデルを用意しておくことで、各音の基本周波数の推定精度が向上することが期待される。しかしながら、基本周波数の推定精度を高めるために、多数の音モデルを用いるとなると、そのような多数の音モデルを作成するのに多大な工数を要し、また、多数の音モデルを記憶させるための記憶容量を音分析装置内に確保しなければならないという問題がある。
【００１３】
この発明は、以上説明した事情に鑑みてなされたものであり、記憶する音モデルの個数が比較的少なくて済み、かつ、高い推定精度で基本周波数の推定を行うことができる音分析装置およびプログラムを提供することを目的としている。
【課題を解決するための手段】
【００１４】
この発明は、楽器から発音される複数種類の音の高調波構造を各々定義した複数種類の音モデルを記憶する記憶手段と、前記記憶手段に記憶された複数種類の音モデルを各々の基本周波数に従って序列化し、序列化された複数種類の音モデルに対して、基本周波数に基づく補間処理を施し、序列化された各音モデルの中間の基本周波数に対応した複数種類の音モデルを生成する補間手段と、前記記憶手段に記憶された複数種類の音モデルおよび前記補間手段により生成された複数種類の音モデルを用いて、各種の高調波構造および基本周波数を有する複数の音モデルを重み付け加算した混合分布を構成し、この混合分布が入力音響信号の周波数成分の分布となるように、各音モデルに対する重み値を最適化し、最適化された各音モデルの重み値を前記入力音響信号の元である音源の音の基本周波数の確率密度関数として推定する確率密度関数推定手段と、前記基本周波数の確率密度関数に基づいて前記入力音響信号における１または複数の音源の音の基本周波数を推定して出力する基本周波数推定手段とを具備することを特徴とする音分析装置並びにコンピュータを該音分析装置として機能させるコンピュータプログラムを提供する。
【００１５】
かかる発明によれば、基本周波数に基づく音モデルの補間が行われ、音モデルが補充された状態で、音モデルを用いた基本周波数の推定が行われるので、記憶手段に記憶させる音モデルの個数が比較的少ない場合であっても高い精度で基本周波数の推定を行うことができる。
【発明を実施するための最良の形態】
【００１６】
以下、図面を参照し、この発明の実施の形態を説明する。
【００１７】
＜Ａ．第１実施形態＞
＜全体構成＞
図１は、この発明の第１実施形態による音分析プログラムの処理内容を示す図である。この音分析プログラムは、自然界から音響信号を取得する収音機能、ＣＤ等の記録媒体から音楽の音響信号を再生する再生機能またはネットワークを介して音楽の音響信号を取得する通信機能等の音響信号取得機能を備えたパーソナルコンピュータ等のコンピュータにインストールされて実行される。本実施形態による音分析プログラムを実行するコンピュータは、本実施形態による音分析装置として機能する。
【００１８】
本実施形態による音分析プログラムは、音響信号取得機能を介して取得されたモノラルの音楽音響信号に対し、その中のある音源の音高を推定する。その最も重要な例として、ここではメロディラインとベースラインを推定する。メロディは他よりも際立って聞こえる単音の系列、ベースはアンサンブル中で最も低い単音の系列であり、その時間的な変化の軌跡をそれぞれメロディラインＤｍ（ｔ）、ベースラインＤｂ（ｔ）と呼ぶ。時刻tにおける基本周波数Ｆ０をＦｉ（ｔ）（ｉ＝ｍ，ｂ）、振幅をＡｉ（ｔ）とすると、これらは以下のように表される。
【数１】
2008058753.tif 000003
【数２】
2008058753.tif 000004
【００１９】
このメロディラインＤｍ（ｔ）およびベースラインＤｂ（ｔ）を入力音響信号から得るための手段として、音分析プログラムは、瞬時周波数の算出１、周波数成分の候補の抽出２、周波数帯域の制限３、メロディラインの推定４ａおよびベースラインの推定４ｂ並びに音モデル補間処理５の各処理を含む。また、メロディラインの推定４ａおよびベースラインの推定４ｂの各処理は、基本周波数の確率密度関数の推定４１およびマルチエージェントモデルによる基本周波数の継時的な追跡４２の各処理を各々含む。本実施形態において、瞬時周波数の算出１、周波数成分の候補の抽出２、周波数帯域の制限３、メロディラインの推定４ａおよびベースラインの推定４ｂの処理内容は、前掲特許文献１および非特許文献１に開示されたものと基本的に同様である。本実施形態の特徴は、音モデル補間処理５を追加した点にある。以下、本実施形態による音分析プログラムを構成する各処理の内容を説明する。
【００２０】
＜瞬時周波数の算出１＞
この処理では、入力音響信号を複数のＢＰＦからなるフィルタバンクに与え、フィルタバンクの各ＢＰＦの出力信号について、位相の時間微分である瞬時周波数（Flanagan, J.L. and Golden, R.M.: Phase Vocoder, The BellSystem
Technical J., Vol.45, pp.1493-1509 (1966)参照）を計算する。ここでは、上記Flanaganの手法を用い、短時間フーリエ変換(STFT)の出力をフィルタバンク出力と解釈して、効率良く瞬時周波数を計算する。入力音響信号ｘ(ｔ)に対する窓関数ｈ(ｔ)を用いたＳＴＦＴが式（３）および（４）により与えられるとき、瞬時周波数λ（ω，ｔ）は式（５）により求めることができる。
【００２１】
【数３】
2008058753.tif 000005
【数４】
2008058753.tif 000006
【数５】
2008058753.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)は、次式によって抽出することができる。
【００２５】
【数６】
2008058753.tif 000008
これらの周波数成分のパワーは、Ψf(t)の各周波数におけるＳＴＦＴパワースペクトルの値として得られるため、周波数成分のパワー分布関数Ψp(t)(ω)を次のように定義することができる。
【数７】
2008058753.tif 000009
【００２６】
＜周波数帯域の制限３＞
この処理では、抽出した周波数成分に重み付けすることで、周波数帯域を制限する。ここでは、メロディラインとベースライン用に、２種類のＢＰＦを用意する。メロディライン用ＢＰＦは、典型的なメロディラインの主要な基本波成分および高調波成分の多くを通過させることができ、かつ、基本周波数付近の重複が頻繁に起きる周波数帯域をある程度遮断する。一方、ベースライン用ＢＰＦは、典型的なベースラインの主要な基本周波数成分および高調波成分の多くを通過させることができ、かつ、他の演奏パートがベースラインよりも優勢になるような周波数帯域をある程度遮断する。
【００２７】
本実施形態では、以下、対数スケールの周波数をcentの単位(本来は音高差(音程)を表す尺度)で表し、Ｈｚで表された周波数ｆＨｚを、次のようにcentで表された周波数ｆｃｅｎｔに変換する。 
【数８】
2008058753.tif 000010
【数９】
2008058753.tif 000011
平均律の半音は１００ｃｅｎｔに、１オクターブは１２００ｃｅｎｔに相当する。
【００２８】
周波数ｘ ｃｅｎｔでのＢＰＦの周波数応答をＢＰＦｉ（ｘ）（ｉ＝ｍ，ｂ）とし、周波数成分のパワー分布関数をΨ’ｐ（ｔ）（ｘ）とすると、ＢＰＦを通過した周波数成分はＢＰＦｉ（ｘ）Ψ’ｐ（ｔ）（ｘ）と表すことができる。ただし、Ψ’ｐ（ｔ）（ｘ）は、周波数軸がｃｅｎｔで表されていることを除けばΨｐ（ｔ）（ω）と同じ関数である。ここで、次の段階の準備として、ＢＰＦを通過した周波数成分の確率密度関数ｐΨ（ｔ）（ｘ）を定義する。
【数１０】
2008058753.tif 000012
【００２９】
ここで、Ｐｏｗ（ｔ）は次式に示すようにＢＰＦを通過した周波数成分のパワーの合計である。
【数１１】
2008058753.tif 000013
【００３０】
＜基本周波数の確率密度関数の推定４１＞
この基本周波数の確率密度関数の推定４１では、ＢＰＦを通過した周波数成分の候補に対し、各高調波構造が相対的にどれくらい優勢かを表す基本周波数の確率密度関数を求める。この基本周波数の確率密度関数の推定４１の処理内容は、非特許文献１に開示された内容となっている。
【００３１】
基本周波数の確率密度関数の推定４１では、上述した「拡張１」と「拡張２」を実現するために、同一基本周波数に対してＭｉ種類の音モデルがあるものとし（ｉはメロディ用（ｉ＝ｍ）かベース用（ｉ＝ｂ）かを示す）、基本周波数がＦであり、音モデルの種類がｍ番目の種類であり、モデルパラメータμ（ｔ）（Ｆ，ｍ）を持った音モデルｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））を次のように定義する。
【数１２】
2008058753.tif 000014
【数１３】
2008058753.tif 000015
【数１４】
2008058753.tif 000016
【数１５】
2008058753.tif 000017
【００３２】
この音モデルは、基本周波数がＦのときに、その高調波成分がどの周波数にどれくらい現れるかをモデル化したものである。Ｈｉは基本周波数成分も含めた高調波成分の数、Ｗｉ２はガウス分布Ｇ（ｘ；ｘ０，σ）の分散を表す。ｃ（ｔ）（ｈ｜Ｆ，ｍ）は、基本周波数がＦであるｍ番目の音モデルの第ｈ次調波成分の大きさを表し、次式を満たす。
【数１６】
2008058753.tif 000018
【００３３】
ｍ番目の音モデルがある基本周波数Ｆの音モデルとして使用され、基本周波数の確率密度関数の推定に用いられる場合、その基本周波数Ｆの音モデルにおける重みｃ（ｔ）（ｈ｜Ｆ，ｍ）として、上記式（１６）に示すように、総和が１となるように予め定義された重みｃ（ｔ）（ｈ｜Ｆ，ｍ）が用いられる。
【００３４】
基本周波数の確率密度関数の推定４１では、以上のような音モデルを使用し、周波数成分の確率密度関数ｐΨ（ｔ）（ｘ）が、次式で定義されるようなｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））の混合分布モデルｐ（ｘ｜θ（ｔ））から生成されたと考える。
【数１７】
2008058753.tif 000019
【数１８】
2008058753.tif 000020
【数１９】
2008058753.tif 000021
【数２０】
2008058753.tif 000022
ここで、ＦｈｉとＦｌｉは、許容される基本周波数の上限と下限であり、ｗ（ｔ）（Ｆ，ｍ）は次式を満たすような音モデルの重みである。
【数２１】
2008058753.tif 000023
【００３５】
実世界の混合音に対して事前に音源数を仮定することは不可能なため、式（１７）のように、あらゆる基本周波数の可能性を同時に考慮してモデル化することが重要となる。最終的に、モデルｐ（ｘ｜θ（ｔ））から、観測した確率密度関数ｐΨ（ｔ）（ｘ）が生成されたかのようにモデルパラメータθ（ｔ）を推定できれば、その重みｗ（ｔ）（Ｆ，ｍ）は各高調波構造が相対的にどれくらい優勢かを表すため、次式のように基本周波数の確率密度関数ｐＦ０（ｔ）（Ｆ）と解釈することができる。
【数２２】
2008058753.tif 000024
【００３６】
次に、上述した「拡張３」を実現するために、θ（ｔ）の事前分布ｐ０ｉ（θ（ｔ））を、式（２３）のように式（２４）と式（２５）の積で与える。
【数２３】
2008058753.tif 000025
【数２４】
2008058753.tif 000026
【数２５】
2008058753.tif 000027
【００３７】
ここで、ｐ０ｉ（ｗ（ｔ））とｐ０ｉ（μ（ｔ））は、最も起こりやすいパラメータをｗ０ｉ（ｔ）（Ｆ，ｍ）とμ０ｉ（ｔ）（Ｆ，ｍ）としたときに、そこで最大値を取るような単峰性の事前分布である。ただし、Ｚｗ、Ｚμは正規化係数、βｗｉ（ｔ）、βμｉ（ｔ）（Ｆ，ｍ）は、最大値をどれくらい重視した事前分布とするかを決めるパラメータで、０のときに無情報事前分布(一様分布）となる。また、Ｄｗ（ｗ０ｉ（ｔ）；ｗ（ｔ））、Ｄμ（μ０ｉ（ｔ）（Ｆ，ｍ）；μ（ｔ）（Ｆ，ｍ））は、次のようなＫ−Ｌ情報量（Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ’ｓ  ｉｎｆｏｒｍａｔｉｏｎ）である。
【数２６】
2008058753.tif 000028
【数２７】
2008058753.tif 000029
【００３８】
以上から、確率密度関数ｐΨ（ｔ）（ｘ）を観測したときに、そのモデルｐ（ｘ｜θ（ｔ））のパラメータθ（ｔ）を、事前分布ｐ０ｉ（θ（ｔ））に基づいて推定する問題を解けばよいことがわかる。この事前分布に基づくθ（ｔ）の最大事後確率推定量（ＭＡＰ推定量）は、次式を最大化することで得られる。
【数２８】
2008058753.tif 000030
【００３９】
この最大化問題は解析的に解くことが困難なため、前述のＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）アルゴリズムを用いてθ（ｔ）を推定する。ＥＭアルゴリズムは、Ｅステップ（ｅｘｐｅｃｔａｔｉｏｎ  ｓｔｅｐ）とＭステップ（ｍａｘｉｍｉｚａｔｉｏｎ  ｓｔｅｐ）を交互に繰返し適用することで、不完全な観測データ（この場合、ｐΨ（ｔ）（ｘ））から最尤推定をおこなうための反復アルゴリズムである。本実施形態では、ＥＭアルゴリズムを繰り返すことにより、ＢＰＦを通過した周波数成分の確率密度関数ｐΨ（ｔ）（ｘ）を、各種の基本周波数Ｆに対応した複数の音モデルｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））を重み付け加算した混合分布と考える場合において、最も尤もらしい重みのパラメータθ（ｔ）（＝｛ｗ（ｔ）（Ｆ，ｍ），μ（ｔ）（Ｆ，ｍ）｝）を求める。ここで、ＥＭアルゴリズムの各繰り返しでは、パラメータθ（ｔ）（＝｛ｗ（ｔ）（Ｆ，ｍ），μ（ｔ）（Ｆ，ｍ））に関して、古いパラメータ推定値θold（ｔ）（＝｛ｗｏｌｄ（ｔ）（Ｆ，ｍ），μｏｌｄ（ｔ）（Ｆ，ｍ）｝）を更新して新しい（より尤もらしい）パラメータ推定値θｎｅｗ（ｔ）（＝｛ｗｎｅｗ（ｔ）（Ｆ，ｍ），μｎｅｗ（ｔ）（Ｆ，ｍ）｝）を求めていく。θｏｌｄ（ｔ）の初期値には、１つ前の時刻ｔ−１における最終的な推定値を用いる。この古いパラメータ推定値θold（ｔ）から新しいパラメータ推定値θｎｅｗ（ｔ）を求める漸化式は、次のようになる。なお、この漸化式の導出過程は非特許文献１に詳細に説明されているので、そちらを参照されたい。
【数２９】
2008058753.tif 000031
【数３０】
2008058753.tif 000032
【００４０】
上記式（２９）および（３０）におけるｗＭＬ（ｔ）（Ｆ，ｍ）とｃＭＬ（ｔ）（ｈ｜Ｆ，ｍ）は、βｗｉ（ｔ）＝０、βμｉ（ｔ）（Ｆ，ｍ）＝０の無情報事前分布のとき、つまり、最尤推定の場合の推定値であり、次式により与えられる。
【数３１】
2008058753.tif 000033
【数３２】
2008058753.tif 000034
【００４１】
これらの反復計算により、事前分布を考慮した基本周波数の確率密度関数ｐＦ０（ｔ）（Ｆ）が、式（２２）によってｗ（ｔ）（Ｆ，ｍ）から求まる。さらに、すべての音モデルｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））の各高調波成分の大きさの比率ｃ（ｔ）（ｈ｜Ｆ，ｍ）も求まり、「拡張１」〜「拡張３」が実現される。
【００４２】
最も優勢な基本周波数Ｆｉ（ｔ）を決定するには、次式に示すように、基本周波数の確率密度関数ｐＦ０（ｔ）（Ｆ）（式（２２）より、式（２９）〜（３２）を反復計算した最終的な推定値として得られる）を最大にする周波数を求めればよい。
【数３３】
2008058753.tif 000035
こうして得られた周波数を音高とする。これが、本実施形態において、基本周波数の確率密度関数に基づいて入力音響信号における１または複数の音源の音の基本周波数を推定して出力する基本周波数推定手段としての処理である。
【００４３】
＜マルチエージェントモデルによる基本周波数の継時的な追跡４２＞
基本周波数の確率密度関数において、同時に鳴っている音の基本周波数に対応する複数のピークが拮抗すると、それらのピークが確率密度関数の最大値として次々に選ばれてしまうことがあるため、このように単純に求めた結果は安定しないことがある。そこで、本実施形態における基本周波数推定手段としての処理では、大局的な観点から基本周波数を推定するために、基本周波数の確率密度関数の時間変化において複数のピークの軌跡を継時的に追跡し、その中で最も優勢で安定した基本周波数の軌跡を選択する。このような追跡処理を動的で柔軟に制御するために、マルチエージェントモデルを導入する。
【００４４】
マルチエージェントモデルは、１つの特徴検出器と複数のエージェントにより構成される（図２参照）。特徴検出器は、基本周波数の確率密度関数の中で目立つピークを拾い上げる。エージェントは基本的に、それらのピークに駆動されて軌跡を追跡していく。つまり、マルチエージェントモデルは、入力中で目立つ特徴を時間的に追跡する汎用の枠組みである。具体的には、各時刻において以下の処理がおこなわれる。
【００４５】
（１）基本周波数の確率密度関数が求まった後、特徴検出器は目立つピーク（最大ピークに応じて動的に変化する閾値を越えたピーク）を複数検出する。そして、目立つピークのそれぞれについて、周波数成分のパワーの合計Ｐｏｗ（ｔ）も考慮しながら、どれくらい将来有望なピークかを評価する。これは、現在時刻を数フレーム先の時刻とみなして、ピークの軌跡をその時刻まで先読みして追跡することで実現する。
【００４６】
（２）既に生成されたエージェントがあるときは、それらが相互作用しながら、目立つピークをそれに近い軌跡を持つエージェントへと排他的に割り当てる。複数のエージェントが割り当て候補に上がる場合には、最も信頼度の高いエージェントへと割り当てる。
【００４７】
（３）最も有望で目立つピークがまだ割り当てられていないときは、そのピークを追跡する新たなエージェントを生成する。
【００４８】
（４）各エージェントは累積ペナルティを持っており、それが一定の閾値を越えると消滅する。
【００４９】
（５）目立つピークが割り当てられなかったエージェントは、一定のペナルティを受け、基本周波数の確率密度関数の中から自分の追跡する次のピークを直接見つけようとする。もしそのピークも見つからないときは、さらにペナルティを受ける。さもなければ、ペナルティはリセットされる。
【００５０】
（６）各エージェントは、今割り当てられたピークがどれくらい有望で目立つかを表す度合いと、１つ前の時刻の信頼度との重み付き和によって、信頼度を自己評価する。
【００５１】
（７）時刻ｔにおける基本周波数Ｆｉ（ｔ）は、信頼度が高く、追跡しているピークの軌跡に沿ったパワーの合計が大きいエージェントに基づいて決定する。振幅Ａｉ（ｔ）は、基本周波数Ｆｉ（ｔ）の高調波成分等をΨｐ（ｔ）（ω）から抽出して決定する。
【００５２】
＜本実施形態の改良点（音モデル補間処理５）＞
一般に楽器から発音される音のスペクトル形状は音高（基本周波数）に依存して変化する。従って、基本周波数の推定精度を高めるためには、様々な基本周波数を持った音を楽器から収音し、これらの各音から作成した多くの音モデルを用いて、基本周波数の確率密度関数の推定４１を実行した方が好ましい。しかし、そのような多数の音モデルを基本周波数の確率密度関数の推定４１に用いるとなると、それらの多数の音モデルを作成するのには多大な工数を要し、また、多数の音モデルを記憶させるための記憶容量を音分析装置内に確保しなければならないという問題がある。そこで、本実施形態では、次のような改良が行われている。すなわち、音分析装置の記憶装置には、様々な基本周波数に対応した比較的少数の代表的な音モデルのみを各々の基本周波数と対応付けて記憶させ、音分析プログラムの実行時に、この記憶装置に記憶された比較的少数の代表的な音モデルから多数の音モデルを生成し、基本周波数の確率密度関数の推定４１に引き渡すのである。
【００５３】
そして、本実施形態では、記憶装置に記憶された比較的少数の代表的な音モデルから多数の音モデルを生成し、基本周波数の確率密度関数の推定４１に引き渡すための手段として、図１に示す音モデル補間処理５が音分析プログラムに追加されている。この音モデル補間処理は、記憶装置に記憶された複数種類の音モデルを各々の基本周波数に従って序列化し、序列化された複数種類の音モデルに対して、基本周波数に基づく補間処理を施し、序列化された各音モデルの中間の基本周波数に対応した複数種類の音モデルを生成する処理である。本実施形態による音分析プログラムは、その実行開始の初期に、この音モデル補間処理５を実行し、記憶装置に記憶された代表的な音モデルとこの音モデル補間処理５により得られた音モデルとを基本周波数の確率密度関数の推定４１に引き渡すように構成されている。
【００５４】
図３は代表的な音モデルの選出と音モデル補間処理５の具体例を示すものである。この例では、ギターの全フレットにおいて５フレット毎に代表フレットを選び、それらの各代表フレットを指で押さえたときのギター音の音モデルを作成し、代表的な音モデルとして音分析装置の記憶装置に記憶させる。そして、各代表フレットに挟まれた中間フレットに対応したギター音の音モデルは、音モデル補間処理５に生成させる。音モデル補間処理５では、中間フレットに対応した音モデルのｈ次倍音成分（ｈ＝１〜Ｈｉ）を、その中間フレットの低音側の代表フレットに対応した音モデルのｈ次倍音成分（ｈ＝１〜Ｈｉ）とその中間フレットの高音側の代表フレットに対応した音モデルのｈ次倍音成分（ｈ＝１〜Ｈｉ）とから生成する。この音モデル補間処理５に関しては各種の態様が考えられる。ある好ましい態様では、低音側の代表フレットに対応した音モデルの基本周波数をＦａ、ｈ次倍音成分をｃ（ｈ｜Ｆａ、ｍａ）、高音側の代表フレットに対応した音モデルの基本周波数をＦｂ、ｈ次倍音成分をｃ（ｈ｜Ｆｂ、ｍｂ）、中間フレットに対応した音モデルの基本周波数をＦｃ、ｈ次倍音成分をｃ（ｈ｜Ｆｃ、ｍｃ）とした場合に、次式に示す１次補間により中間フレットに対応した音モデルのｈ次倍音成分ｃ（ｈ｜Ｆｃ、ｍｃ）を求める。
【数３４】
2008058753.tif 000036
【００５５】
基本周波数の確率密度関数の推定４１では、このようにして得られる中間フレットに対応した音モデルと元々記憶装置に記憶されていた代表フレットに対応した代表的な音モデルの両方が用いられる。
【００５６】
代表的な音モデルを作成する基本周波数は、楽器の構造により定まる倍音特性の特徴に着目して選出するのが好ましい。具体的には、倍音構造が急激に変化する基本周波数領域においては、密に代表的な音モデルを作成して記憶装置に記憶させるのが効果的である。
例えばギターでは、ある弦の最高フレットまでは倍音構造が連続的に変化するが、それより半音高い音を出すには、より高い音を発音しうる別の弦を奏さねばならず、ここで倍音構造が不連続となる。さらにいえば、弦の開放弦および低いフレット付近では倍音構造の変化は緩やかであり、高いフレット付近では１フレットの違いでも大きく倍音構造が変化するから、これを反映し、低域では粗い間隔で、高域になるほど密に、代表的な音モデルを作成する基本周波数を選ぶと有効である。また、ピアノでは最低音域は１本の弦、低音域は２本の弦、高音域は３本の弦、最高音域は、３本の弦であるが他と異なりミュート機構がない、といったように周波数帯域により弦構造が異なり、これに呼応して音色すなわち倍音構造も特定の周波数で不連続に変化する。そのような不連続点では代表的な音モデルを密に配置すれば、少ない数の音モデルでも、基本周波数の推定精度を高めることが可能となる。
【００５７】
以上説明した本実施形態によれば、記憶する音モデルのデータ量を削減しつつ、音域ごとに異なる音源特性をより詳細に音モデルとして表現し、さらには少数のパラメータ調整により実際の入力音に最適に音モデルの形状を調整することが可能となる。
【００５８】
＜Ｂ．第２実施形態＞
上記第１実施形態では、基本周波数によって音のスペクトル形状が異なることを考慮し、比較的少ない代表的な音モデルからより多くの種類の基本周波数に対応した音モデルを補間処理により生成した。本実施形態では、基本周波数の確率密度関数の推定４１において、各音モデル（代表的な音モデルおよび音モデル補間処理５により得られた音モデル）の種類毎に、当該音モデルが本来有していた基本周波数に合わせて、基本周波数の範囲を設定し、この設定した基本周波数の範囲外の周波数における当該音モデルに対する重み値を制限して、各音モデルに対する重み値の最適化を行う。さらに詳述すると、次の通りである。
【００５９】
まず、ＥＭアルゴリズムによる基本周波数の確率密度関数の推定に関して、音モデルの種類毎に基本周波数の適用範囲を定める。基本周波数の適用範囲の下限Ｆｌｍおよび上限Ｆｈｍの定め方には各種の方法が考えられる。例えばｍ番目の種類の音モデルｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））があるフレット位置におけるギター音の音モデルである場合に、そのフレット位置での基本周波数と隣接する低音側のフレット位置での基本周波数との中間の周波数をＦｌｍとし、そのフレット位置での基本周波数と隣接する高音側のフレット位置での基本周波数との中間の周波数をＦｈｍとしてもよい。あるいはＦｌｍとＦｈｍとの間をもっと広くとり、基本周波数の隣接した各音モデル間で基本周波数の適用範囲をオーバラップさせてもよい。
【００６０】
そして、本実施形態では、ＥＭアルゴリズムの過程において、ｍ番目の種類の音モデルｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））は、その適用範囲から外れるような基本周波数（すなわち、Ｆ＜ＦｌｍまたはＦ＞Ｆｈｍであるような基本周波数）の確率密度の推定に使われないようにする。具体的には、ｍ番目の種類の音モデルｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））は、何ら策を講じないとすると、ＥＭアルゴリズムにおいて多くの種類の基本周波数Ｆに対応した音モデルとして使用され得るが、式（２９）〜（３２）の漸化式を繰り返す際、各基本周波数Ｆに対応した各音モデルｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））のうちＦ＜ＦｌｍまたはＦ＞Ｆｈｍであるような基本周波数Ｆに対応した音モデルについては、それらに対する重み値ｗ（Ｆ，ｍ）の初期値を０にするのである。
【００６１】
このようにすることで、ｍ番目の種類の音モデルｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））は、Ｆ＜ＦｌｍまたはＦ＞Ｆｈｍであるような基本周波数Ｆの確率密度の推定には一切使われないようになる。このような処理を全ての種類の音モデルについて行う。
【００６２】
この態様によれば、音源から発生し得る各音の基本周波数の範囲に合わせて、その音に対応した音モデルの適用可能な基本周波数の範囲（Ｆｌｍ，Ｆｈｍ）を個別的に定義しておくことで、各音の基本周波数の範囲を考慮した適切な基本周波数の推定を行うことができる。
【００６３】
＜他の実施形態＞
以上、この発明の一実施形態について説明したが、この発明には他にも実施形態があり得る。例えば次の通りである。
【００６４】
（１）音モデル補間処理５には、１次補間に限らず、０次補間、スプライン補間など一般に知られる補間法を広く用いることができる。さらに、代表的な音モデルにも誤差が含まれ得るので、自己回帰を用いることにより、代表的な音モデルの各間の音モデルを求めてもよい。
【００６５】
（２）形状の異なった多くの種類の音モデルを用いることは基本周波数の推定精度の向上に寄与するが、同じような形状の音モデルを多数用いたとしても基本周波数の推定精度の向上は期待できない。そこで、音モデル補間処理５では、代表的な音モデルの各間の全てではなく、音モデルの形状がある程度以上変化する一部区間のみについて補間を行うようにしてもよい。
【００６６】
（３）音モデル補間処理５において、周波数領域により異なる補間演算方法で音モデルの補間を行ってもよい。例えば基本周波数の変化に対する音モデルの形状の変化が緩やかな周波数領域では１次補間により音モデルを求め、基本周波数の変化に対する音モデルの形状の変化が比較的急激な周波数領域ではより高次の補間により音モデルを求める、といった態様が考えられる。あるいは補間演算方法自体を変える代わりに、補間演算用のパラメータを周波数領域に応じて変えてもよい。
【００６７】
（４）上記各実施形態では、基本周波数の確率密度関数の推定４１により得られる基本周波数をマルチエージェントに追跡させることにより、最終的な基本周波数を決定したが、基本周波数の確率密度関数の推定４１において誤推定の確率が低く、信頼性の高い推定結果が得られる場合には、マルチエージェントによる追跡を省略してもよい。
【００６８】
（５）上記各実施形態では、音分析装置に「拡張１」（音モデルの多重化）に加えて、「拡張２」（音モデルのパラメータの推定)を導入したが、「拡張２」は導入せず、基本周波数の確率密度関数の推定では、例えば漸化式（２９）および（３０）のうち漸化式（２９）のみを逐次演算し、音モデルに対する重みｗ（Ｆ，ｍ）の更新のみを行うようにしてもよい。
【００６９】
（６）上記各実施形態では、音分析装置に「拡張３」（モデルパラメータに関する事前分布の導入）を導入したが、音分析装置はこれを導入しない構成としてもよい。
【図面の簡単な説明】
【００７０】
【図１】この発明の第１実施形態である音分析プログラムの処理内容を示す図である。
【図２】１つの特徴検出器と複数のエージェントにより構成されるマルチエージェントモデルによる基本周波数の経時的な追跡を示す図である。
【図３】同実施形態における代表的な音モデルの選出例と音モデル補間処理の内容を示す図である。
【符号の説明】
【００７１】
１……瞬時周波数の算出、２……周波数成分の候補の算出、３……周波数帯域の制限、４ａ……メロディラインの推定、４ｂ……ベースラインの推定、４１……基本周波数の確率密度関数の推定、４２……マルチエージェントモデルによる基本周波数の継時的な追跡、５……音モデル補間処理。
【図１】
2008058753.tif 000037
【図２】
2008058753.tif 000038
【図３】
2008058753.tif 000039