(57)【要約】      （修正有）【課題】従来よりも少ない演算回数で基本周波数の確率密度関数の重みと高調波成分の大きさとを推定する。【解決手段】2006285052.tif 000123中の1200ｌｏｇ２ｈとｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を事前に計算してコンピュータのメモリに格納する。ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が０に近い基本周波数Ｆに関してだけ上記式の演算を実行する。これにより従来よりも大幅に演算回数を少なくすることができて、演算時間の短縮化を可能にする。【選択図】  図３
2006285052.tif 000002
【特許請求の範囲】
【請求項１】
入力される混合音に含まれる周波数成分を観測し、観測した周波数成分を下記（ａ）の対数周波数ｘ上の確率密度関数として表現するものとし、
【数１】
2006285052.tif 000078
前記観測した周波数成分の確率密度関数から、下記（ｂ）で表現される基本周波数Ｆの確率密度関数を求める過程において、音モデルの多重化、音モデルのパラメータの推定及びモデルパラメータに関する事前分布の導入を採用し、
【数２】
2006285052.tif 000079
前記音モデルの多重化では、同一基本周波数に対してＭ種類の音モデルがあるものとして前記基本周波数がＦのｍ番目の音モデルの確率密度関数をｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））と表現し、ただし、μ（ｔ）（Ｆ，ｍ）は前記ｍ番目の音モデルの高調波成分の大きさの比率を表すモデルパラメータであり、
前記音モデルのパラメータの推定では、前記観測した周波数成分の確率密度関数が、下記（ｃ）で定義した混合分布モデルｐ（ｘ｜θ（ｔ））から生成されたものと考え、ただしω（ｔ）（Ｆ，ｍ）は基本周波数がＦのｍ番目の音モデルの重みであり、
【数３】
2006285052.tif 000080
ただし上記（ｃ）において、θ（ｔ）は前記音モデルの重みω（ｔ）（Ｆ，ｍ）と前記音モデルの高調波成分の大きさの比率μ（ｔ）（Ｆ，ｍ）を包含したモデルパラメータθ（ｔ）＝｛ω（ｔ），μ（ｔ）｝であり、ω（ｔ）＝｛ω（ｔ）（Ｆ，ｍ）｜Ｆｌ≦Ｆ≦Ｆｈ，ｍ＝１，…，Ｍ｝であり、μ（ｔ）＝｛μ（ｔ）（Ｆ，ｍ）｜Ｆｌ≦Ｆ≦Ｆｈ，ｍ＝１，…，Ｍ｝であり、ここでＦｌは許容される基本周波数の下限であり、Ｆｈは許容される基本周波数の上限であり、
前記基本周波数Ｆの確率密度関数を下記（ｄ）の解釈により前記重みω（ｔ）（Ｆ，ｍ）から求め、
【数４】
2006285052.tif 000081
前記モデルパラメータに関する事前分布の導入では、前記モデルパラメータθ（ｔ）に関する事前分布に基づいて前記モデルパラメータθ（ｔ）の最大事後確率推定量をＥＭ（Expectation-Maximization）アルゴリズムを用いて推定し、この推定から事前分布を考慮した上記（ｂ）の基本周波数Ｆの確率密度関数と解釈できる前記重みω（ｔ）（Ｆ，ｍ）とすべての音モデルの確率密度関数ｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））のμ（ｔ）（Ｆ，ｍ）が表す第ｈ次高調波成分の大きさｃ（ｔ）（ｈ｜Ｆ，ｍ）（ｈ＝１，…，Ｈ）とを求めるために用いる下記（ｅ）及び（ｆ）で表される二つのパラメータ推定値を求めるための式を定め、ただしＨは前記基本周波数の周波数成分も含めた高調波成分の数であり、
【数５】
2006285052.tif 000082
【数６】
2006285052.tif 000083
また上記（ｅ）及び（ｆ）中において、下記（ｇ）及び（ｈ）は、下記（ｉ）及び（ｊ）が０になる無情報事前分布のときの最尤推定値であり、
【数７】
2006285052.tif 000084
【数８】
2006285052.tif 000085
【数９】
2006285052.tif 000086
【数１０】
2006285052.tif 000087
上記（ｅ）及び（ｆ）中において、下記（ｋ）は前記重みω（ｔ）（Ｆ，ｍ）の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、下記（ｌ）は前記モデルパラメータμ（ｔ）（Ｆ，ｍ）の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、また上記（ｉ）は下記（ｋ）をどれぐらい重視した事前分布とするかを決めるパラメータであり、上記（ｊ）は下記（ｌ）をどれぐらい重視した事前分布とするかを決めるパラメータであり、
【数１１】
2006285052.tif 000088
【数１２】
2006285052.tif 000089
また上記（ｇ）及び（ｈ）において、ω´（ｔ）（Ｆ，ｍ）及びμ´（ｔ）（Ｆ，ｍ）は上記（ｅ）及び（ｆ）を反復計算する際の一つ前の古いパラメータ推定値であり、ηは基本周波数であり、υは何番目の音モデルであるかを示し、
上記（ｅ）及び（ｆ）の二つのパラメータ推定値を求めるための式を用いた反復計算により、上記（ｂ）の基本周波数の確率密度関数と解釈できる前記重みω（ｔ）（Ｆ，ｍ）とすべての音モデルの確率密度関数ｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））の前記モデルパラメータμ（ｔ）（Ｆ，ｍ）が表す第ｈ次高調波成分の大きさｃ（ｔ）（ｈ｜Ｆ，ｍ）とをコンピュータを利用した演算により求めることにより基本周波数の音高を推定する音高推定方法であって、
上記（ｅ）及び（ｆ）を上記（ｇ）及び（ｈ）によって前記コンピュータを利用して演算するために、上記（ｇ）で表現される推定値を示す計算式中の分子を下記（ｍ）に示すｘの関数として展開し、但し下記（ｍ）に示す式においてω´（ｔ）（Ｆ，ｍ）は古い重みであり、ｃ´（ｔ）（ｈ｜Ｆ，ｍ）は古い第ｈ次高調波成分の大きさの比率であり、Ｈは前記基本周波数の周波数成分も含めた高調波成分の数であり、ｍはＭ種類の音モデルの中の何番目の音モデルかを示すものであり、Ｗは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、
【数１３】
2006285052.tif 000090
上記（ｍ）中の1200ｌｏｇ２ｈとｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を事前に計算して前記コンピュータのメモリに格納し、
上記（ｅ）及び（ｆ）の二つのパラメータ推定値を求めるための式を事前に定めた回数分だけ反復計算するために、上記（ｇ）及び（ｈ）の式の計算では、前記観測した周波数成分の確率密度関数の離散化後の各周波数ｘに対して以下の第１の演算処理をＮｘ回実行し、但しＮｘはｘの定義域の範囲の離散化数であり、
前記第１の演算処理では、Ｍ種類の音モデルについてそれぞれ以下の第２の演算処理を実行して、上記（ｍ）式の演算結果を求め、上記（ｍ）式の演算結果を基本周波数Ｆとｍ番目の音モデルに関して積分して上記（ｇ）及び（ｈ）の式の分母を求め前記観測した周波数成分の確率密度関数を上記（ｇ）及び（ｈ）の式に代入して上記（ｇ）及び（ｈ）の演算を行い、
前記第２の演算処理では、基本周波数を含む高調波成分の数Ｈだけ、第３の演算処理を実行して下記（ｎ）の式の演算結果を求め、下記（ｎ）の式の演算結果をＨ個加算して上記（ｍ）の式の演算結果を求め、
【数１４】
2006285052.tif 000091
前記第３の演算処理では、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が０に近い基本周波数Ｆに関して第４の演算処理をＮａ回実行して上記（ｎ）の式の演算結果を求め、但しＮａはｘ−（Ｆ＋1200ｌｏｇ２ｈ）が充分０に近い範囲の離散化後のＦの個数を表す小さい正の整数であり、
前記第４の演算処理では、事前に前記メモリに格納したｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を用いて下記（ｏ）の式を求め、
【数１５】
2006285052.tif 000092
上記（ｏ）の式に古い重みω´（ｔ）（Ｆ，ｍ）をかけて上記（ｎ）の式の演算結果を求めることを特徴とする音高推定方法。
【請求項２】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅をｄとしたときに、（３Ｗ／ｄ）より小さいもしくは近い正の整数ｂを求めて、前記Ｎａを（２ｂ＋１）回と決定し、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取り、ここでＷは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、αは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数である請求項１に記載の音高推定方法。
【請求項３】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅をｄとしたときに、（３Ｗ／ｄ）より小さいもしくは近い正の整数ｂを求めて、前記Ｎａを（２ｂ＋１）回と決定し、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取り、ここでＷは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、αは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数であり、
前記メモリには、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取るときのｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］の値が事前に格納されている請求項１に記載の音高推定方法。
【請求項４】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅が２０ｃｅｎｔで前記Ｗが１７ｃｅｎｔのときに、前記Ｎａを５と定め、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−２＋α，−１＋α，０＋α，１＋α，２＋αの値を取り、ここでαは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数である請求項１に記載の音高推定方法。
【請求項５】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅が２０ｃｅｎｔで前記Ｗが１７ｃｅｎｔのときに、前記Ｎａを５と定め、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−２＋α，−１＋α，０＋α，１＋α，２＋αの値を取り、ここでαは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数であり、
前記メモリには、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−２＋α，−１＋α，０＋α，１＋α，２＋αの値を取るときのｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］の値が事前に格納されている請求項１に記載の音高推定方法。
【請求項６】
入力される混合音に含まれる周波数成分を観測し、観測した周波数成分を下記（ａ）の対数周波数ｘ上の確率密度関数として表現するものとし、
【数１６】
2006285052.tif 000093
前記観測した周波数成分の確率密度関数から、下記（ｂ）で表現される基本周波数Ｆの確率密度関数を求める過程において、音モデルの多重化、音モデルのパラメータの推定及びモデルパラメータに関する事前分布の導入を採用し、
【数１７】
2006285052.tif 000094
前記音モデルの多重化では、同一基本周波数に対してＭ種類の音モデルがあるものとして前記基本周波数がＦのｍ番目の音モデルの確率密度関数をｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））と表現し、ただし、μ（ｔ）（Ｆ，ｍ）は前記ｍ番目の音モデルの高調波成分の大きさの比率を表すモデルパラメータであり、
前記音モデルのパラメータの推定では、前記観測した周波数成分の確率密度関数が、下記（ｃ）で定義した混合分布モデルｐ（ｘ｜θ（ｔ））から生成されたものと考え、ただしω（ｔ）（Ｆ，ｍ）は基本周波数がＦのｍ番目の音モデルの重みであり、
【数１８】
2006285052.tif 000095
ただし上記（ｃ）において、θ（ｔ）は前記音モデルの重みω（ｔ）（Ｆ，ｍ）と前記音モデルの高調波成分の大きさの比率μ（ｔ）（Ｆ，ｍ）を包含したモデルパラメータθ（ｔ）＝｛ω（ｔ），μ（ｔ）｝であり、ω（ｔ）＝｛ω（ｔ）（Ｆ，ｍ）｜Ｆｌ≦Ｆ≦Ｆｈ，ｍ＝１，…，Ｍ｝であり、μ（ｔ）＝｛μ（ｔ）（Ｆ，ｍ）｜Ｆｌ≦Ｆ≦Ｆｈ，ｍ＝１，…，Ｍ｝であり、ここでＦｌは許容される基本周波数の下限であり、Ｆｈは許容される基本周波数の上限であり、
前記基本周波数Ｆの確率密度関数を下記（ｄ）の解釈により前記重みω（ｔ）（Ｆ，ｍ）から求め、
【数１９】
2006285052.tif 000096
前記モデルパラメータに関する事前分布の導入では、前記モデルパラメータθ（ｔ）に関する事前分布に基づいて前記モデルパラメータθ（ｔ）の最大事後確率推定量をＥＭ（Expectation-Maximization）アルゴリズムを用いて推定し、この推定から事前分布を考慮した上記（ｂ）の基本周波数Ｆの確率密度関数と解釈できる前記重みω（ｔ）（Ｆ，ｍ）とすべての音モデルの確率密度関数ｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））のμ（ｔ）（Ｆ，ｍ）が表す第ｈ次高調波成分の大きさのｃ（ｔ）（ｈ｜Ｆ，ｍ）（ｈ＝１，…，Ｈ）とを求めるために用いる下記（ｅ）及び（ｆ）で表される二つのパラメータ推定値を求めるための式を定め、ただしＨは前記基本周波数の周波数成分も含めた高調波成分の数であり、
【数２０】
2006285052.tif 000097
【数２１】
2006285052.tif 000098
また上記（ｅ）及び（ｆ）中において、下記（ｇ）及び（ｈ）は、下記（ｉ）及び（ｊ）が０になる無情報事前分布のときの最尤推定値であり、
【数２２】
2006285052.tif 000099
【数２３】
2006285052.tif 000100
【数２４】
2006285052.tif 000101
【数２５】
2006285052.tif 000102
上記（ｅ）及び（ｆ）中において、下記（ｋ）は前記重みω（ｔ）（Ｆ，ｍ）の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、下記（ｌ）は前記モデルパラメータμ（ｔ）（Ｆ，ｍ）の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、また上記（ｉ）は下記（ｋ）をどれぐらい重視した事前分布とするかを決めるパラメータであり、上記（ｊ）は下記（ｌ）をどれぐらい重視した事前分布とするかを決めるパラメータであり、
【数２６】
2006285052.tif 000103
【数２７】
2006285052.tif 000104
また上記（ｇ）及び（ｈ）において、ω´（ｔ）（Ｆ，ｍ）及びμ´（ｔ）（Ｆ，ｍ）は上記（ｅ）及び（ｆ）を反復計算する際の一つ前の古いパラメータ推定値であり、ηは基本周波数であり、υは何番目の音モデルであるかを示し、
上記（ｅ）及び（ｆ）の二つのパラメータ推定値を求めるための式を用いた反復計算により、上記（ｂ）の基本周波数の確率密度関数と解釈できる前記重みω（ｔ）（Ｆ，ｍ）とすべての音モデルの確率密度関数ｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））の前記モデルパラメータμ（ｔ）（Ｆ，ｍ）が表す第ｈ次高調波成分の大きさｃ（ｔ）（ｈ｜Ｆ，ｍ）とをコンピュータ内に以下の機能を実現する手段を構成して演算により求めることにより基本周波数の音高を推定する音高推定装置であって、
上記（ｅ）及び（ｆ）を上記（ｇ）及び（ｈ）によって前記コンピュータを利用して演算するために、上記（ｇ）で表現される推定値を示す計算式中の分子を下記（ｍ）に示すｘの関数として展開する手段を備え、但し下記（ｍ）に示す式においてω´（ｔ）（Ｆ，ｍ）は古い重みであり、ｃ´（ｔ）（ｈ｜Ｆ，ｍ）は古い第ｈ次高調波成分の大きさの比率であり、Ｈは前記基本周波数の周波数成分も含めた高調波成分の数であり、ｍはＭ種類の音モデルの中の何番目の音モデルかを示すものであり、Ｗは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、
【数２８】
2006285052.tif 000105
上記（ｍ）中の1200ｌｏｇ２ｈとｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を事前に計算して前記コンピュータのメモリに格納する手段を備え、
上記（ｅ）及び（ｆ）の二つのパラメータ推定値を求めるための式を事前に定めた回数分だけ反復計算するために、上記（ｇ）及び（ｈ）の式の計算では、前記観測した周波数成分の確率密度関数の離散化後の各周波数ｘに対して以下の第１の演算処理をＮｘ回実行する第１の演算処理手段を備え、但しＮｘはｘの定義域の範囲の離散化数であり、
前記第１の演算処理手段では、Ｍ種類の音モデルについてそれぞれ以下の第２の演算処理を実行して、上記（ｍ）式の演算結果を求め、上記（ｍ）式の演算結果を基本周波数Ｆとｍ番目の音モデルに関して積分して上記（ｇ）及び（ｈ）の式の分母を求め前記観測した周波数成分の確率密度関数を上記（ｇ）及び（ｈ）の式に代入して上記（ｇ）及び（ｈ）の演算を行う第２の演算処理手段を備え、
前記第２の演算処理手段では、基本周波数を含む高調波成分の数Ｈだけ、第３の演算処理を実行して下記（ｎ）の式の演算結果を求め、下記（ｎ）の式の演算結果をＨ個加算して上記（ｍ）の式の演算結果を求める第３の演算処理手段を備え、
【数２９】
2006285052.tif 000106
前記第３の演算処理手段では、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が０に近い基本周波数Ｆに関して第４の演算処理をＮａ回実行して上記（ｎ）の式の演算結果を求める第４の演算処理手段を備え、但しＮａはｘ−（Ｆ＋1200ｌｏｇ２ｈ）が充分０に近い範囲の離散化後のＦの個数を表す小さい正の整数であり、
前記第４の演算処理手段では、事前に前記メモリに格納したｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を用いて下記（ｏ）の式を求め、
【数３０】
2006285052.tif 000107
上記（ｏ）の式に古い重みω´（ｔ）（Ｆ，ｍ）をかけて上記（ｎ）の式の演算結果を求めることを特徴とする音高推定装置。
【請求項７】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅をｄとしたときに、（３Ｗ／ｄ）より小さいもしくは近い正の整数ｂを求めて、前記Ｎａを（２ｂ＋１）回と決定し、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取り、ここでＷは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、αは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数である請求項６に記載の音高推定装置。
【請求項８】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅をｄとしたときに、（３Ｗ／ｄ）より小さいもしくは近い正の整数ｂを求めて、前記Ｎａを（２ｂ＋１）回と決定し、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取り、ここでＷは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、αは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数であり、
前記メモリには、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取るときのｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］の値が事前に格納されている請求項６に記載の音高推定装置。
【請求項９】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅が２０ｃｅｎｔで前記Ｗが１７ｃｅｎｔのときに、前記Ｎａを５と定め、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−２＋α，−１＋α，０＋α，１＋α，２＋αの値を取り、ここでαは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数である請求項６に記載の音高推定装置。
【請求項１０】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅が２０ｃｅｎｔで前記Ｗが１７ｃｅｎｔのときに、前記Ｎａを５と定め、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−２＋α，−１＋α，０＋α，１＋α，２＋αの値を取り、ここでαは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数であり、
前記メモリには、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−２＋α，−１＋α，０＋α，１＋α，２＋αの値を取るときのｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］の値が事前に格納されている請求項６に記載の音高推定装置。
【請求項１１】
入力される混合音に含まれる周波数成分を観測し、観測した周波数成分を下記（ａ）の対数周波数ｘ上の確率密度関数として表現するものとし、
【数３１】
2006285052.tif 000108
前記観測した周波数成分の確率密度関数から、下記（ｂ）で表現される基本周波数Ｆの確率密度関数を求める過程において、音モデルの多重化、音モデルのパラメータの推定及びモデルパラメータに関する事前分布の導入を採用し、
【数３２】
2006285052.tif 000109
前記音モデルの多重化では、同一基本周波数に対してＭ種類の音モデルがあるものとして前記基本周波数がＦのｍ番目の音モデルの確率密度関数をｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））と表現し、ただし、μ（ｔ）（Ｆ，ｍ）は前記ｍ番目の音モデルの高調波成分の大きさの比率を表すモデルパラメータであり、
前記音モデルのパラメータの推定では、前記観測した周波数成分の確率密度関数が、下記（ｃ）で定義した混合分布モデルｐ（ｘ｜θ（ｔ））から生成されたものと考え、ただしω（ｔ）（Ｆ，ｍ）は基本周波数がＦのｍ番目の音モデルの重みであり、
【数３３】
2006285052.tif 000110
ただし上記（ｃ）において、θ（ｔ）は前記音モデルの重みω（ｔ）（Ｆ，ｍ）と前記音モデルの高調波成分の大きさの比率μ（ｔ）（Ｆ，ｍ）を包含したモデルパラメータθ（ｔ）＝｛ω（ｔ），μ（ｔ）｝であり、ω（ｔ）＝｛ω（ｔ）（Ｆ，ｍ）｜Ｆｌ≦Ｆ≦Ｆｈ，ｍ＝１，…，Ｍ｝であり、μ（ｔ）＝｛μ（ｔ）（Ｆ，ｍ）｜Ｆｌ≦Ｆ≦Ｆｈ，ｍ＝１，…，Ｍ｝であり、ここでＦｌは許容される基本周波数の下限であり、Ｆｈは許容される基本周波数の上限であり、
前記基本周波数Ｆの確率密度関数を下記（ｄ）の解釈により前記重みω（ｔ）（Ｆ，ｍ）から求め、
【数３４】
2006285052.tif 000111
前記モデルパラメータに関する事前分布の導入では、前記モデルパラメータθ（ｔ）に関する事前分布に基づいて前記モデルパラメータθ（ｔ）の最大事後確率推定量をＥＭ（Expectation-Maximization）アルゴリズムを用いて推定し、この推定から事前分布を考慮した上記（ｂ）の基本周波数Ｆの確率密度関数と解釈できる前記重みω（ｔ）（Ｆ，ｍ）とすべての音モデルの確率密度関数ｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））のμ（ｔ）（Ｆ，ｍ）が表す第ｈ次高調波成分の大きさのｃ（ｔ）（ｈ｜Ｆ，ｍ）（ｈ＝１，…，Ｈ）とを求めるために用いる下記（ｅ）及び（ｆ）で表される二つのパラメータ推定値を求めるための式を定め、ただしＨは前記基本周波数の周波数成分も含めた高調波成分の数であり、
【数３５】
2006285052.tif 000112
【数３６】
2006285052.tif 000113
また上記（ｅ）及び（ｆ）中において、下記（ｇ）及び（ｈ）は、下記（ｉ）及び（ｊ）が０になる無情報事前分布のときの最尤推定値であり、
【数３７】
2006285052.tif 000114
【数３８】
2006285052.tif 000115
【数３９】
2006285052.tif 000116
【数４０】
2006285052.tif 000117
上記（ｅ）及び（ｆ）中において、下記（ｋ）は前記重みω（ｔ）（Ｆ，ｍ）の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、下記（ｌ）は前記モデルパラメータμ（ｔ）（Ｆ，ｍ）の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、また上記（ｉ）は下記（ｋ）をどれぐらい重視した事前分布とするかを決めるパラメータであり、上記（ｊ）は下記（ｌ）をどれぐらい重視した事前分布とするかを決めるパラメータであり、
【数４１】
2006285052.tif 000118
【数４２】
2006285052.tif 000119
また上記（ｇ）及び（ｈ）において、ω´（ｔ）（Ｆ，ｍ）及びμ´（ｔ）（Ｆ，ｍ）は上記（ｅ）及び（ｆ）を反復計算する際の一つ前の古いパラメータ推定値であり、ηは基本周波数であり、υは何番目の音モデルであるかを示し、
上記（ｅ）及び（ｆ）の二つのパラメータ推定値を求めるための式を用いた反復計算により、上記（ｂ）の基本周波数の確率密度関数と解釈できる前記重みω（ｔ）（Ｆ，ｍ）とすべての音モデルの確率密度関数ｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））の前記モデルパラメータμ（ｔ）（Ｆ，ｍ）が表す第ｈ次高調波成分の大きさｃ（ｔ）（ｈ｜Ｆ，ｍ）とをコンピュータを利用した演算により求めるため前記コンピュータにインストールされて前記コンピュータ内に必要な機能を実現するための音高推定用プログラムであって、
上記（ｅ）及び（ｆ）を上記（ｇ）及び（ｈ）によって前記コンピュータを利用して演算するために、上記（ｇ）で表現される推定値を示す計算式中の分子を下記（ｍ）に示すｘの関数として展開する機能、但し下記（ｍ）に示す式においてω´（ｔ）（Ｆ，ｍ）は古い重みであり、ｃ´（ｔ）（ｈ｜Ｆ，ｍ）は古い第ｈ次高調波成分の大きさの比率であり、Ｈは前記基本周波数の周波数成分も含めた高調波成分の数であり、ｍはＭ種類の音モデルの中の何番目の音モデルかを示すものであり、Ｗは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、
【数４３】
2006285052.tif 000120
上記（ｍ）中の1200ｌｏｇ２ｈとｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を事前に計算して前記コンピュータのメモリに格納する機能、
上記（ｅ）及び（ｆ）の二つのパラメータ推定値を求めるための式を事前に定めた回数分だけ反復計算するために、上記（ｇ）及び（ｈ）の式の計算では、前記観測した周波数成分の確率密度関数の離散化後の各周波数ｘに対して以下の第１の演算処理をＮｘ回実行する機能、但しＮｘはｘの定義域の範囲の離散化数であり、
前記第１の演算処理では、Ｍ種類の音モデルについてそれぞれ以下の第２の演算処理を実行して、上記（ｍ）式の演算結果を求め、上記（ｍ）式の演算結果を基本周波数Ｆとｍ番目の音モデルに関して積分して上記（ｇ）及び（ｈ）の式の分母を求め前記観測した周波数成分の確率密度関数を上記（ｇ）及び（ｈ）の式に代入して上記（ｇ）及び（ｈ）の演算を行う機能、
前記第２の演算処理では、基本周波数を含む高調波成分の数Ｈだけ、第３の演算処理を実行して下記（ｎ）の式の演算結果を求め、下記（ｎ）の式の演算結果をＨ個加算して上記（ｍ）の式の演算結果を求める機能、
【数４４】
2006285052.tif 000121
前記第３の演算処理では、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が０に近い基本周波数Ｆに関して第４の演算処理をＮａ回実行して上記（ｎ）の式の演算結果を求める機能、但しＮａはｘ−（Ｆ＋1200ｌｏｇ２ｈ）が充分０に近い範囲の離散化後のＦの個数を表す小さい正の整数であり、
前記第４の演算処理では、事前に前記メモリに格納したｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を用いて下記（ｏ）の式を求め、
【数４５】
2006285052.tif 000122
上記（ｏ）の式に古い重みω´（ｔ）（Ｆ，ｍ）をかけて上記（ｎ）の式の演算結果を求める機能を前記コンピュータ内に実現するための音高推定用プログラム。
【請求項１２】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅をｄとしたときに、（３Ｗ／ｄ）より小さいもしくは近い正の整数ｂを求めて、前記Ｎａを（２ｂ＋１）回と決定し、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取り、ここでＷは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、αは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数である請求項１１に記載の音高推定用プログラム。
【請求項１３】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅をｄとしたときに、（３Ｗ／ｄ）より小さいもしくは近い正の整数ｂを求めて、前記Ｎａを（２ｂ＋１）回と決定し、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取り、ここでＷは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、αは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数であり、
前記メモリには、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取るときのｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］の値が事前に格納されている請求項１１に記載の音高推定用プログラム。
【請求項１４】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅が２０ｃｅｎｔで前記Ｗが１７ｃｅｎｔのときに、前記Ｎａを５と定め、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−２＋α，−１＋α，０＋α，１＋α，２＋αの値を取り、ここでαは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数である請求項１１に記載の音高推定用プログラム。
【請求項１５】
前記対数周波数ｘと前記基本周波数Ｆの離散化幅が２０ｃｅｎｔで前記Ｗが１７ｃｅｎｔのときに、前記Ｎａを５と定め、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−２＋α，−１＋α，０＋α，１＋α，２＋αの値を取り、ここでαは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数であり、
前記メモリには、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−２＋α，−１＋α，０＋α，１＋α，２＋αの値を取るときのｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］の値が事前に格納されている請求項１１に記載の音高推定用プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、混合音中の各構成音（基本周波数）の音高と音量を推定する音高推定方法及び装置並びに音高推定用プログラムに関するものである。
【背景技術】
【０００２】
ＣＤ等による実世界の音響信号は、事前に音源数を仮定することが不可能な混合音である。このような混合音中では、周波数成分が頻繁に重複する上に、基本周波数成分が存在しないような音も存在する。しかし、従来の音高推定技術の多くは、少数の音源数を仮定し、周波数成分を局所的に追跡したり、基本周波数成分の存在に依存していた。そのために、前述の実世界の混合音には適用できなかった。
【０００３】
そこで発明者は、特許第３４１３６３４号公報（特許文献１）に示される「音高推定方法及び装置」と題する発明を提案した。この従来の発明では、入力の混合音には、あらゆる基本周波数（本願明細書中で抽象的に使用される「音高」に相当するもの）の音が様々な音量で同時に含まれていると考える。そしてこの発明では、統計的手法を利用するために、入力の周波数成分を確率密度関数（観測した分布）で表現し、各音の高調波構造に対応する確率分布を音モデルとして導入する。そして、周波数成分の確率密度関数が、対象とするあらゆる基本周波数の音モデルの混合分布モデル（重み付き和のモデル）から生成されたと考える。この混合分布中の各音モデルの重みは、各高調波構造が相対的にどれくらい優勢かを表すことから、基本周波数の確率密度関数と呼ぶ（混合分布中において音モデルが優勢になればなるほど、そのモデルの基本周波数の確率が高くなる）。この重みの値（すなわち基本周波数の確率密度関数）は、ＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）アルゴリズム（Ｄｅｍｐｓｔｅｒ，Ａ．Ｐ．，Ｌａｉｒｄ,Ｎ.Ｍ  ａｎｄ  Ｒｕｂｉｎ、Ｄ．Ｂ．：Ｍａｘｉｍｕｍ  ｌｉｋｅｌｉｈｏｏｄ  ｆｒｏｍ  ｉｎｃｏｍｐｌｅｔｅ  ｄａｔａ  ｖｉａ  ｔｈｅ  ＥＭ  ａｌｇｏｒｉｔｈｍ，Ｊ．Ｒｏｙ．Ｓｔａｔ．Ｓｏｃ．Ｂ，Ｖｏｌ．３９，Ｎｏ．１，ｐｐ．１−３８（１９７７））を用いることで推定できる。こうして求めた基本周波数の確率密度関数は、混合音中の構成音が、どの音高でどれぐらいの音量で鳴っているかを表している。
【０００４】
また発明者は、この従来の発明を発展または拡張させる技術を非特許文献１及び２の二つの論文に公表している。まず非特許文献１は、２００１年５月に公表された「A PREDOMINANT-FO ESTIMATION METHOD FOR CD RECORDINGS: MAP ESTIMATION USING EM ALGORITHM FOR ADAPTIVE TONE MODELS」と題する論文で、The 2001 IEEE International Conference on Acoustics, Speech, and Signal Processingの予稿集Ｖの3365-3368頁に掲載された。また非特許文献２は、２００４年９月に公表された「A real-time music-scene-description system: predominant-FO estimation for detecting melody and bass lines in real-world audio signals」と題する論文で、Speech Communication 43(2004)の第３１１頁〜第３２９頁に掲載された。これら二つの非特許文献で提案された拡張は、音モデルの多重化と、音モデルのパラメータの推定と、モデルパラメータに関する事前分布の導入である。なおこれらの拡張は後に詳しく説明する。
【特許文献１】特許第３４１３６３４号公報
【非特許文献１】「APREDOMINANT-FO ESTIMATION METHOD FOR CD RECORDINGS: MAP ESTIMATION USING EM ALGORITHMFOR ADAPTIVE TONE MODELS」と題する論文（The 2001 IEEE International Conference onAcoustics, Speech, and Signal Processingの予稿集Ｖの3365-3368頁）
【非特許文献２】「Areal-time music-scene-description system: predominant-FO estimation fordetecting melody and bass lines in real-world audio signals」と題する論文で（Speech Communication43(2004)の第３１１頁〜第３２９頁）
【発明の開示】
【発明が解決しようとする課題】
【０００５】
しかしながら上記拡張された技術をコンピュータを用いて実現して、基本周波数の確率密度関数の重みと高調波成分の大きさとを推定するためには、演算回数が非常に多くなり、高速演算機能を有するコンピュータを用いなければ、短い時間で推定結果が得られない問題があった。
【０００６】
本発明の目的は、従来よりも少ない演算回数で基本周波数の確率密度関数の重みと高調波成分の大きさとを推定できる音高推定方法及び装置並びに音高推定用プログラムを提供することにある。
【課題を解決するための手段】
【０００７】
本発明の音高推定方法では、次のようにして基本周波数の確率密度関数の重みと高調波成分の大きさとを推定する。
【０００８】
まず入力される混合音に含まれる周波数成分を観測し、観測した周波数成分を下記（ａ）の対数周波数ｘ上の確率密度関数として表現する。
【数４６】
2006285052.tif 000003
【０００９】
そして非特許文献１及び２に開示した技術（音モデルの多重化、音モデルのパラメータの推定及びモデルパラメータに関する事前分布の導入）を、上記（ａ）で表現される観測した周波数成分の確率密度関数から、下記（ｂ）で表現される基本周波数Ｆの確率密度関数を求める過程においてを採用する。
【数４７】
2006285052.tif 000004
【００１０】
まず音モデルの多重化では、同一基本周波数に対してＭ種類の音モデルがあるものとして基本周波数がＦのｍ番目の音モデルの確率密度関数をｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））と表現する。ただし、μ（ｔ）（Ｆ，ｍ）は、ｍ番目の音モデルの高調波成分の大きさの比率を表すモデルパラメータである。
【００１１】
また音モデルのパラメータの推定では、観測した周波数成分の確率密度関数が、下記（ｃ）で定義した混合分布モデルｐ（ｘ｜θ（ｔ））から生成されたものと考える。ただしω（ｔ）（Ｆ，ｍ）は基本周波数がＦのｍ番目の音モデルの重みである。
【数４８】
2006285052.tif 000005
【００１２】
なお上記（ｃ）の式において、θ（ｔ）は音モデルの重みω（ｔ）（Ｆ，ｍ）と音モデルの高調波成分の大きさの比率μ（ｔ）（Ｆ，ｍ）を包含したモデルパラメータθ（ｔ）＝｛ω（ｔ），μ（ｔ）｝であり、ω（ｔ）＝｛ω（ｔ）（Ｆ，ｍ）｜Ｆｌ≦Ｆ≦Ｆｈ，ｍ＝１，…，Ｍ｝であり、μ（ｔ）＝｛μ（ｔ）（Ｆ，ｍ）｜Ｆｌ≦Ｆ≦Ｆｈ，ｍ＝１，…，Ｍ｝であり、ここでＦｌは許容される基本周波数の下限であり、Ｆｈは許容される基本周波数の上限である。
そして上記（ｂ）の基本周波数Ｆの確率密度関数を下記（ｄ）の解釈により重みω（ｔ）（Ｆ，ｍ）から求める。
【数４９】
2006285052.tif 000006
【００１３】
さらにモデルパラメータに関する事前分布の導入では、モデルパラメータθ（ｔ）に関する事前分布に基づいてモデルパラメータθ（ｔ）の最大事後確率推定量をＥＭ（Expectation-Maximization）アルゴリズムを用いて推定する。そして、この推定から事前分布を考慮した上記（ｂ）の基本周波数Ｆの確率密度関数と解釈できる重みω（ｔ）（Ｆ，ｍ）と、すべての音モデルの確率密度関数ｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））のμ（ｔ）（Ｆ，ｍ）が表す第ｈ次高調波成分の大きさｃ（ｔ）（ｈ｜Ｆ，ｍ）（ｈ＝１，…，Ｈ）とを求めるために用いる下記（ｅ）及び（ｆ）で表される二つのパラメータ推定値を求めるための式を定める。ただしＨは基本周波数の周波数成分も含めた高調波成分の数である。
【数５０】
2006285052.tif 000007
【数５１】
2006285052.tif 000008
【００１４】
上記（ｅ）及び（ｆ）中において、下記（ｇ）及び（ｈ）は、下記（ｉ）及び（ｊ）が０になる無情報事前分布のときの最尤推定値である。
【数５２】
2006285052.tif 000009
【数５３】
2006285052.tif 000010
【数５４】
2006285052.tif 000011
【数５５】
2006285052.tif 000012
【００１５】
また上記（ｅ）及び（ｆ）中において、下記（ｋ）は重みω（ｔ）（Ｆ，ｍ）の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、下記（ｌ）はモデルパラメータμ（ｔ）（Ｆ，ｍ）の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、また上記（ｉ）は下記（ｋ）をどれぐらい重視した事前分布とするかを決めるパラメータであり、上記（ｊ）は下記（ｌ）をどれぐらい重視した事前分布とするかを決めるパラメータである。
【数５６】
2006285052.tif 000013
【数５７】
2006285052.tif 000014
【００１６】
また上記（ｇ）及び（ｈ）の式において、ω´（ｔ）（Ｆ，ｍ）及びμ´（ｔ）（Ｆ，ｍ）は上記（ｅ）及び（ｆ）を反復計算する際の一つ前の古いパラメータ推定値であり、ηは基本周波数であり、υは何番目の音モデルであるかを示す。
【００１７】
本発明が改良の対象とする音高推定方法では、上記（ｅ）及び（ｆ）の二つのパラメータ推定値を求めるための式を用いた反復計算により、上記（ｂ）の基本周波数の確率密度関数と解釈できる重みω（ｔ）（Ｆ，ｍ）とすべての音モデルの確率密度関数ｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））のモデルパラメータμ（ｔ）（Ｆ，ｍ）が表す第ｈ次高調波成分の大きさｃ（ｔ）（ｈ｜Ｆ，ｍ）とをコンピュータを利用した演算により求めることにより基本周波数の音高を推定する。
【００１８】
本発明では、上記（ｅ）及び（ｆ）を上記（ｇ）及び（ｈ）によってコンピュータを利用して演算するために、次のようにする。まず上記（ｇ）で表現される推定値を示す計算式中の分子を下記（ｍ）に示すｘの関数として展開する。但し下記（ｍ）に示す式においてω´（ｔ）（Ｆ，ｍ）は古い重みであり、ｃ´（ｔ）（ｈ｜Ｆ，ｍ）は古い第ｈ次高調波成分の大きさの比率であり、Ｈは基本周波数の周波数成分を含めた高調波成分の数であり、ｍはＭ種類の音モデルの中の何番目の音モデルかを示すものであり、Ｗは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差である。
【数５８】
2006285052.tif 000015
【００１９】
上記（ｍ）中の1200ｌｏｇ２ｈとｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を事前に計算してコンピュータのメモリに格納する。
【００２０】
上記（ｅ）及び（ｆ）の二つのパラメータ推定値を求めるための式を事前に定めた回数分だけ反復計算するために、上記（ｇ）及び（ｈ）の式の計算では、観測した周波数成分の確率密度関数の離散化後の各周波数ｘに対して以下の第１の演算処理をＮｘ回実行する。但しＮｘはｘの定義域の範囲の離散化数である。
【００２１】
第１の演算処理では、Ｍ種類の音モデルについてそれぞれ以下の第２の演算処理を実行して、上記（ｍ）式の演算結果を求める。そして上記（ｍ）式の演算結果を基本周波数Ｆとｍ番目の音モデルに関して積分して上記（ｇ）及び（ｈ）の式の分母を求め、観測した周波数成分の確率密度関数を上記（ｇ）及び（ｈ）の式に代入して上記（ｇ）及び（ｈ）の演算を行う。
【００２２】
また第２の演算処理では、基本周波数を含む高調波成分の数Ｈだけ、第３の演算処理を実行して下記（ｎ）の式の演算結果を求め、下記（ｎ）の式の演算結果をＨ個加算して上記（ｍ）の式の演算結果を求める。
【数５９】
2006285052.tif 000016
【００２３】
さらに第３の演算処理では、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が０に近い基本周波数Ｆに関して第４の演算処理をＮａ回実行して上記（ｎ）の式の演算結果を求める。但しＮａはｘ−（Ｆ＋1200ｌｏｇ２ｈ）が充分０に近い範囲の離散化後のＦの個数を表す小さい正の整数である。
【００２４】
そして第４の演算処理では、事前にメモリに格納したｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を用いて下記（ｏ）の式を求める。
【数６０】
2006285052.tif 000017
【００２５】
最後に、上記（ｏ）の式に古い重みω´（ｔ）（Ｆ，ｍ）をかけて上記（ｎ）の式の演算結果を求める。
【００２６】
本発明の方法によれば、事前にメモリに格納したｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を用いることができるので、演算回数を減らすことができる。特に本発明では、第４の演算処理の回数をＮａ回と少なくして、上記（ｍ）の式の演算結果を求めても、演算精度が低下することがないことを見いだしたことを根拠として、第４の演算処理回数を制限する。その結果、従来よりも大幅に演算回数を少なくすることができて、演算時間の短縮化を可能にした。
【００２７】
なお対数周波数ｘと基本周波数Ｆの離散化幅をｄとしたときには、（３Ｗ／ｄ）より小さいもしくは近い正の整数ｂを求めて、Ｎａを（２ｂ＋１）回と決定し、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取るようにすればよい。そしてメモリには、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取るときのｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］の値を事前に格納しておくのが好ましい。ここで前述のＷは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差である。またαは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数である。なお、ここで（３Ｗ／ｄ）の分子の３は３以外の任意の正の整数でもよく、小さいほど演算回数が少なくなる。
【００２８】
より具体的には、対数周波数ｘと基本周波数Ｆの離散化幅が２０ｃｅｎｔ(半音の音高差１００ｃｅｎｔの五分の一)でＷが１７ｃｅｎｔのときは、第４の演算処理を行う回数Ｎａは５回にすることが好ましい。この場合には、離散化して計算するときにｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−２＋α，−１＋α，０＋α，１＋α，２＋αの５通りの値を取る。そして、αは（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの０．５以下の小数である。このようにすることによって演算回数を大幅に少なくすることができる。なおこの場合においては、メモリには、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−２＋α，−１＋α，０＋α，１＋α，２＋αの値を取るときのｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］の値を事前に格納しておくのが好ましい。また、1200ｌｏｇ２ｈも事前に計算して格納しておくとよい。これらの格納によって、演算回数を更に低減できる。
【００２９】
本発明の音高推定装置では、前述の本発明の音高推定方法をコンピュータを用いて実施する。そのために本発明の音高推定装置は、上記（ｇ）で表現される推定値を示す計算式中の分子を上記（ｍ）に示すｘの関数として展開する手段と、上記（ｍ）中の1200ｌｏｇ２ｈとｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を事前に計算してコンピュータのメモリに格納する手段と、前述の第１の演算処理を実行する第１の演算処理手段と、前述の第２の演算処理を実行する第２の演算処理手段と、前述の第３の演算処理を実行する第３の演算処理手段と、前述の第４の演算処理を実行する第４の演算処理手段とを備えている。
【００３０】
また本発明の音高推定用プログラムは、本発明の音高推定方法をコンピュータを用いて実施するためにコンピュータにインストールされるものである。本発明の音高推定用プログラムは、上記（ｇ）で表現される推定値を示す計算式中の分子を上記（ｍ）に示すｘの関数として展開する機能と、上記（ｍ）中の1200ｌｏｇ２ｈとｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を事前に計算してコンピュータのメモリに格納する機能と、前述の第１の演算処理を実行する機能と、前述の第２の演算処理を実行する機能と、前述の第３の演算処理を実行する機能と、前述の第４の演算処理を実行する機能をコンピュータ内に実現するように構成されている。
【発明の効果】
【００３１】
本発明によれば、音源数を仮定せず、周波数成分の局所的な追跡も行わず、また基本周波数成分の存在を前提とせずに、音高推定する際に、従来よりも大幅に演算回数を少なくすることができて、演算時間を短縮化することができる。
【発明を実施するための最良の形態】
【００３２】
以下図面を参照しながら本発明の音高推定方法及びプログラムの実施の形態の一例を詳細に説明する。まず本発明の方法の実施の形態の一例について説明する前提として、特許第３４１３６３４号の発明を拡張した前述の非特許文献１及び２に提案された公知の三つの拡張方法について簡単に説明する。
【００３３】
［拡張方法１］音モデルを多重化
特許第３４１３６３４号公報に記載の発明では、同一基本周波数には一つの音モデルしか用意していなかった。しかしながら実際には、ある基本周波数に、異なる高調波構造を持つ音が入れ替わり立ち替わり現れることがある。そこで、同一基本周波数に対して複数の音モデルを用意し、それらの混合分布でモデル化する。具体的な方法は、後に詳しく説明する。
【００３４】
［拡張方法２］音モデルのパラメータを推定
特許第３４１３６３４号公報に記載の従来の音モデルでは、各高調波成分の大きさの比率を固定していた（ある理想的な音モデルを仮定していた）。しかしながらこれは実世界の混合音中の高調波構造とは必ずしも一致しておらず、精度向上のためには、さらに改善をする余地が残されていた。そこで拡張方法２では、音モデルの高調波成分の大きさの比率もモデルパラメータに加え音モデルのパラメータを各時刻においてＥＭアルゴリズムで推定することとした。具体的な方法は後に説明する。
【００３５】
［拡張方法３］モデルパラメータに関する事前分布を導入
特許第３４１３６３４号公報に記載の従来の方法では、音モデルの重み（基本周波数の確率密度関数）に関する事前知識は仮定していなかった。しかし本発明を様々な用途に適用していく上で、たとえ事前に基本周波数がどの周波数の近傍にあるかを与えてでも、より誤検出の少ない基本周波数を求めたいというような応用が考えられる。例えば、演奏分析やビブラート分析等の目的では、楽曲をへッドホンから聴取しながらの歌唱や楽器演奏によって、各時刻におけるおおよその基本周波数を事前知識として用意しておき、実際の楽曲中のより正確な基本周波数を得ることが求められている。そこで、従来のモデルパラメータの最尤推定の枠組みを拡張し、モデルパラメータに関する事前分布に基づいて最大事後確率推定（ＭＡＰ推定：Ｍａｘｉｍｕｍ  Ａ  Ｐｏｓｔｅｒｉｏｒｉ  Ｐｒｏｂａｂｉｌｉｔｙ  Ｅｓｔｉｍａｔｉｏｎ）を行う。その際、［拡張方法２］でモデルパラメータに加えた、音モデルの高調波成分の大きさの比率に関する事前分布も導入する。具体的な方法は後に説明する。
【００３６】
拡張方法１〜３について、式を用いて更に具体的に説明する。まず入力される混合音（入力音響信号）に含まれる観測した周波数成分の確率密度関数を下記（１）で表現する。
【数６１】
2006285052.tif 000018
【００３７】
そして上記（１）式の周波数成分の確率密度関数から、下記（２）式で表現される基本周波数Ｆの確率密度関数を求める過程における、具体的な拡張方法を述べる。
【数６２】
2006285052.tif 000019
【００３８】
上記（１）式の観測した周波数成分の確率密度関数は、混合音（入力音響信号）から例えばマルチレートフィルタバンク（Vetterli,M.: A Theory of Multirate Filter Banks, IEEE Trans. on ASSP, Vol.ASSP-35, No.3, pp.356-372 (1987) 参照）によって求めることができる。なおこのマルチレートフィルタバンクについては、特許第３４１３６３４号公報の図２及び前述の非特許文献２に示されたＦｉｇ．３にバイナリーツリー状のフィルタバンクの構成の一例とその詳細が説明されている。ここで、（１）式及び（２）式中のｔはフレームシフト（１０ｍｓｅｃ）を時間単位とする時刻であり、ｘとＦはｃｅｎｔの単位で表された対数スケールの対数周波数と基本周波数である。なお、Ｈｚで表された周波数ｆＨｚは、下記（３）式によりｃｅｎｔで表された周波数ｆｃｅｎｔに変換されるものとする。
【数６３】
2006285052.tif 000020
【００３９】
そこで前述の［拡張方法１］と［拡張方法２］を実現するために、同一基本周波数に対してＭ種類の音モデルがあるものとし、基本周波数がＦのｍ番目の音モデルの確率密度関数ｐ（ｘ｜Ｆ、ｍ、μ（ｔ）（Ｆ、ｍ））にモデルパラメータμ（ｔ）（Ｆ、ｍ）を導入する。
【００４０】
なお以下に説明する（４）式から（５１）式までについては、前述の非特許文献１において（２）式〜（３６）式としてすでに公表されているので参照されたい。
【００４１】
基本周波数がＦのｍ番目の音モデルの確率密度関数ｐ（ｘ｜Ｆ、ｍ、μ（ｔ）（Ｆ、ｍ））は、次式のように表されるものとする。
【数６４】
2006285052.tif 000021
【数６５】
2006285052.tif 000022
【数６６】
2006285052.tif 000023
【数６７】
2006285052.tif 000024
【００４２】
上記（４）〜（７）の式は、基本周波数がＦのときに、その高調波成分がどの周波数にどれくらい現れるかをモデル化したものである（図１）。上記式においては、Ｈは基本周波数Ｆの周波数成分も含めた高調波成分の数、Ｗはガウス分布Ｇ（ｘ；ｘ０、σ）の標準偏差を表している。またｃ（ｔ）（ｈ｜Ｆ、ｍ）は、第ｈ次高調波成分の大きさを表し、次式を満たすものとする。
【数６８】
2006285052.tif 000025
【００４３】
そして、上記（１）式で表される観測した周波数成分の確率密度関数が、次式で定義されるような、ｐ（ｘ｜Ｆ、ｍ、μ（ｔ）（Ｆ、ｍ））の混合分布モデルｐ（ｘ｜θ（ｔ））から生成されたと考える。
【数６９】
2006285052.tif 000026
【数７０】
2006285052.tif 000027
【数７１】
2006285052.tif 000028
【数７２】
2006285052.tif 000029
【００４４】
上記（１１）及び（１２）式において、ＦｈとＦｌは、許容される基本周波数の上限と下限であり、ｗ（ｔ）（Ｆ、ｍ）は、次式を満たすような、音モデルの重みである。
【数７３】
2006285052.tif 000030
【００４５】
実世界の混合音に対して事前に音源数を仮定することは不可能なため、上記（６）式のように、あらゆる基本周波数の可能性を同時に考慮してモデル化することが重要となる。最終的に、モデルｐ（ｘ｜θ（ｔ））から、観測した確率密度関数［上記（１）式］が生成されたかのようにモデルパラメータθ（ｔ）を推定できれば、その重みｗ（ｔ）（Ｆ、ｍ）は各高調波構造が相対的にどれくらい優勢かを表すことになる。そのため、次式のように基本周波数Ｆの確率密度関数を解釈することができる。
【数７４】
2006285052.tif 000031
【００４６】
次に前述の［拡張方法３］の事前分布の導入を行う。［拡張方法３］を実現するために、θ（ｔ）の事前分布ｐ０ｉ（θ（ｔ））は、下記の式（１９）のように下記式（２０）と下記式（２１）の積で与えられる。下記の式（１９）〜（２１）に示されるｐ０ｉ（ω（ｔ））とｐ０ｉ（μ（ｔ））は、最も起りやすいパラメータを
【数７５】
2006285052.tif 000032
と
【数７６】
2006285052.tif 000033
としたとき（ただし、式（１６）は
【数７７】
2006285052.tif 000034
である)に、そこで最大値を取るような単峰性の事前分布である。ただし、Ｚｗ、Ｚμは正規化係数であり、
【数７８】
2006285052.tif 000035
の二つは、最大値をどれくらい重視した事前分布とするかを決めるパラメータであり、０のときに無情報事前分布（一様分布）となる。
【数７９】
2006285052.tif 000036
【数８０】
2006285052.tif 000037
【数８１】
2006285052.tif 000038
上記（２０）式中の下記の式
【数８２】
2006285052.tif 000039
と上記（２１）式中の下記の式
【数８３】
2006285052.tif 000040
は、次のようなＫ−Ｌ情報量（Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ’ｓ  ｉｎｆｏｒｍａｔｉｏｎ）である。
【数８４】
2006285052.tif 000041
【数８５】
2006285052.tif 000042
【００４７】
以上の説明から、上記（１）式の確率密度関数を観測したときに、そのモデルｐ（ｘ｜θ（ｔ））のパラメータθ（ｔ）を、事前分布ｐ０ｉ（θ（ｔ））に基づいて推定する問題を解けばよいことがわかる。この事前分布ｐ０ｉ（θ（ｔ））に基づくパラメータθ（ｔ）の最大事後確率推定量（ＭＡＰ推定量）は、次式を最大化することで得られる。
【数８６】
2006285052.tif 000043
【００４８】
しかしながらこの最大化問題は解析的に解くことが困難なため、ＥＭアルゴリズム（Ｄｅｍｐｓｔｅｒ，Ａ．Ｐ．，Ｌａｉｒｄ,Ｎ.Ｍ  ａｎｄ  Ｒｕｂｉｎ、Ｄ．Ｂ．：Ｍａｘｉｍｕｍ  ｌｉｋｅｌｉｈｏｏｄ  ｆｒｏｍ  ｉｎｃｏｍｐｌｅｔｅ  ｄａｔａ  ｖｉａ  ｔｈｅ  ＥＭ  ａｌｇｏｒｉｔｈｍ，Ｊ．Ｒｏｙ．Ｓｔａｔ．Ｓｏｃ．Ｂ，Ｖｏｌ．３９，Ｎｏ．１，ｐｐ．１−３８（１９７７））を用いてθ（ｔ）を推定する。ＥＭアルゴリズムは、不完全な観測データから最尤推定を行うために用いられることが多いが、最大事後確率推定の場合にも適用できる。最尤推定では、平均対数尤度の条件付き期待値を求めるＥステップ（ｅｘｐｅｃｔａｔｉｏｎ  ｓｔｅｐ）とその最大化をおこなうＭステップ（ｍａｘｉｍｉｚａｔｉｏｎ  ｓｔｅｐ）を交互に繰り返す。しかし最大事後確率推定の場合には、条件付き期待値に事前分布の対数を加えたものの最大化を繰り返す。ここでは各繰り返しにおいて、古いパラメータ推定値θ’（ｔ）＝｛ｗ’（ｔ）、μ’（ｔ）｝を更新して新しいパラメータ推定値（下記式（２７））
を求めていく。
【数８７】
2006285052.tif 000044
【００４９】
対数周波数ｘにおいて観測した各周波数成分が、どの基本周波数のどの音モデルのどの倍音から生成されたのかを表す隠れ変数Ｆ、ｍ、ｈを導入して、ＥＭアルゴリズムを以下のように定式化することができる。
【００５０】
（Ｅステップ）
最尤推定の場合には、平均対数尤度の条件付き期待値Ｑ（θ（ｔ）｜θ´（ｔ））を求めるが、最大事後確率推定の場合には、それにｌｏｇ  ｐ０ｉ（θ（ｔ））を加えたＱＭＡＰ（θ（ｔ）｜θ（ｔ））を求める。
【数８８】
2006285052.tif 000045
【数８９】
2006285052.tif 000046
【００５１】
上記式において、条件付き期待値ＥＦ、ｍ、ｈ［ａ｜ｂ］は、条件ｂにより決定される確率分布を持つ隠れ変数Ｆ、ｍ、ｈに関する、ａの期待値を意味する。
【００５２】
（Ｍステップ）
ＱＭＡＰ（θ（ｔ）｜θ´（ｔ））をθ（ｔ）の関数として最大化して、更新後の新しい推定値
【数９０】
2006285052.tif 000047
を下記の式で得る。
【数９１】
2006285052.tif 000048
Ｅステップにおいて、上記（２９）式は下記のようになる。
【数９２】
2006285052.tif 000049
上記式中の完全データの対数尤度は
【数９３】
2006285052.tif 000050
で与えられる。また、ｌｏｇ  ｐ０ｉ（θ（ｔ））は、
【数９４】
2006285052.tif 000051
となる。
【００５３】
次に、Ｍステップに関しては、上記（３１）式が、上記（８）式と上記（１３）式を条件とする条件付き変分問題となっている。この問題は、Ｌａｇｒａｎｇｅの乗数λｗ、λμを導入し、次のＥｕｌｅｒ−Ｌａｇｒａｎｇｅの微分方程式によって解くことができる。
【数９５】
2006285052.tif 000052
【数９６】
2006285052.tif 000053
これより、
【数９７】
2006285052.tif 000054
【数９８】
2006285052.tif 000055
が得られる。これらの式において、Ｌａｇｒａｎｇｅの乗数は（８）式、（１３）式から
【数９９】
2006285052.tif 000056
【数１００】
2006285052.tif 000057
と定まり、ｐ（Ｆ、ｍ、ｈ｜ｘ、θ’（ｔ））、ｐ（Ｆ、ｍ｜ｘ、θ’（ｔ））はベイズの定理から、
【数１０１】
2006285052.tif 000058
【数１０２】
2006285052.tif 000059
となる。以上から、新しいパラメータ推定値
【数１０３】
2006285052.tif 000060
と
【数１０４】
2006285052.tif 000061
を求める式は次のようになる。
【数１０５】
2006285052.tif 000062
【数１０６】
2006285052.tif 000063
式中の
【数１０７】
2006285052.tif 000064
と
【数１０８】
2006285052.tif 000065
は、
【数１０９】
2006285052.tif 000066
の無情報事前分布のとき、つまり最尤推定の場合の推定値である。
【数１１０】
2006285052.tif 000067
【数１１１】
2006285052.tif 000068
【００５４】
これらの反復計算により、事前分布を考慮した上記（２）式の基本周波数の確率密度関数が、上記（１４）式によって重みｗ（ｔ）（Ｆ、ｍ）から求まる。さらに、すべての音モデルｐ（ｘ｜Ｆ、ｍ、μ（ｔ）（Ｆ、ｍ））の各高調波成分の大きさの比率ｃ（ｔ）（ｈ｜Ｆ、ｍ）も求まる。これにより［拡張方法１］〜［拡張方法３］が実現される。
【００５５】
上記のように拡張した音高推定手法をコンピュータを用いて実行するためには、上記（４５）式と上記（４６）式の反復計算が必要となる。しかし、これらの式の反復計算では、上記（５０）式と上記（５１）式の計算量が多いため、計算能力の限られた（計算速度が遅い）コンピュータでは、この式をそのまま計算しようとすると計算時間が非常に長くなってしまうという問題が生じる。
【００５６】
計算時間が非常に長くなる理由を説明する。最初に、上記（５０）式を用いて計算結果を普通に求めるときには、どのような計算が必要かを説明する。まず、上記（５０）式の計算では、対象となる範囲のＦとｍに関して、上記（５０）式の右辺の積分内の分子
【数１１２】
2006285052.tif 000069
をｘの関数として計算する（上記（４）式〜（７）式により展開する）。ここで、説明のための計算例として、ｘの定義域の範囲を３６０個（ＮＸ個）に離散化し、ＦｌからＦｈまでの範囲のＦを３００（ＮＦ個）に離散化して計算するものと仮定する。また、音モデルの数Ｍは３、高調波成分の数Ｈは１６とする。このとき、上記（５２）式を計算するには、下記の（５３）式の計算を１６回繰り返すことになる。
【数１１３】
2006285052.tif 000070
【００５７】
上記（５０）式の右辺の積分内の分子を求めるためには、あるｘに関する上記（５２）式の計算を一回行う。そして上記（５０）式の右辺の積分内の分母を求めるためには、Ｆとｍに関して、３００×３回（ＮＦ×Ｍ回）、（５２）式の計算を繰り返す必要がある。
【００５８】
さらに、そのｘを定義域範囲内で３６０通り変化させて積分するため、
【数１１４】
2006285052.tif 000071
を求めるには、分母は１６×（３００×３）×３６０回、分子は１６×３６０回、上記（５３）式の計算を繰り返す必要がある。分母はＦとｍを変化させても共通なので再度計算する必要はないが、分子はすべてのＦ（３００通り）とｍ（３通り）について求める必要がある。そのため、最終的に、分母と分子のそれぞれで１６×（３００×３）×３６０回（Ｈ×ＮＦ×Ｍ×ＮＸ回、計５１８４０００回）、上記（５３）式の計算を繰り返すことになる。ここで、分子を分母よりも前に計算しておけば、分母はその分子を合計すれば求めることができる。よって、分子と分母を共に求めたとしても、上記（５３）式の計算は５１８４０００回繰り返すことになる。
【００５９】
そこで本発明は、この計算時間を以下のようにして大幅に削減して高速化する。以下本発明の方法により、上記の通常計算方法を高速化した高速計算方法を図２及び図３に示した本発明のプログラムのアルゴリズムを示すフローチャートを参照しながら説明する。まず、上記（５０）式の計算では、対象となる範囲のＦとｍに関して、上記（５０）式の右辺の積分内の分子を、ｘの関数として上記（５２）式により計算する。
【００６０】
図２に示すように、上記（５２）式中の1200ｌｏｇ２ｈとｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を事前に計算してコンピュータのメモリに格納する。そして上記（４５）式及び（４６）式の二つのパラメータ推定値を求めるための式を事前に定めた回数分だけ（もしくは収束するまで）反復計算するために、上記（５０）式及び（５１）式の計算では、図３に示すように、上記（４７）式及び（４８）式を０で初期化した後、観測した周波数成分の確率密度関数の各対数周波数ｘに対して以下の第１の演算処理をＮｘ回実行する。但しＮｘはｘの定義域の範囲の離散化数である。
【００６１】
第１の演算処理では、Ｍ種類の音モデルについてそれぞれ以下の第２の演算処理を実行して、上記（５２）式の演算結果を求める。そして上記（５２）式の演算結果を基本周波数Ｆとｍ番目の音モデルに関して積分して上記（５０）式及び（５１）式の分母を求め、観測した周波数成分の確率密度関数を上記（５０）式及び（５１）式に代入して上記（５０）式及び（５１）式の演算を行う。
【００６２】
第２の演算処理では、基本周波数の周波数成分を含む高調波成分の数Ｈだけ、第３の演算処理を実行して後に説明する下記の（５５）式の演算結果を求める。そして（５５）式の演算結果をＨ個加算して上記（５２）式の演算結果を求める。
【数１１５】
2006285052.tif 000072
【００６３】
上記（５５）式は、（５１）式の右辺の積分内の分子を、対象となる範囲のＦとｍ、ｈに関してｘの関数として計算するものである。（５５）式は（５２）式から
【数１１６】
2006285052.tif 000073
を取り除いたものである。
【００６４】
前述の第３の演算処理では、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が０に近い基本周波数Ｆに関して第４の演算処理をＮａ回実行して上記（５５）式の演算結果を求める。但し本発明では、Ｎａはｘ−（Ｆ＋1200ｌｏｇ２ｈ）が充分０に近い範囲のＦの個数を表す小さい正の整数としている。後に説明するように好ましくは、対数周波数ｘと基本周波数Ｆの離散化幅ｄが２０ｃｅｎｔ(半音の音高差１００ｃｅｎｔの五分の一)で、前述のガウス分布の標準偏差Ｗが１７ｃｅｎｔのときは、このＮａは５である。
【００６５】
第４の演算処理では、事前にメモリに格納したｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を用いて上記（５３）式の演算を行う。そして上記（５３）式に古い重みω´（ｔ）（Ｆ，ｍ）をかけて上記（５５）式の演算結果を求める。このようにして本発明では、音高を推定する。
【００６６】
上記をより具体的な例を用いて説明する。
まず上記（５２）式中の
【数１１７】
2006285052.tif 000074
の演算を、ｘと（Ｆ＋1200ｌｏｇ２ｈ）の差が大きくなると急速に０に近づくことを利用し、その差が一定範囲内（本実施例では、ｘとＦの離散化幅が２０ｃｅｎｔ(半音の音高差１００ｃｅｎｔの五分の一）でＷが１７ｃｅｎｔのときは、例えば±２以内の５回（＝Ｎａ）とする)だけ計算する。
【００６７】
ここで、ある対数周波数ｘに関して、（５０）式の右辺の積分内の分母を計算することを考える。上記の計算範囲の制限により、（Ｆ＋1200ｌｏｇ２ｈ）の近傍のｘでのみ上記（５７）式を計算し、それ以外では０とみなして計算しないことにする。このようにすると、ある対数周波数ｘを起点に考えると、（５０）式の右辺の積分内の分母を求めるのに、（５３）式の計算を１６×３００×３回繰り返す必要はなく、１６×５×３回（Ｈ×Ｎａ×Ｍ回）繰り返すだけでよいことになる。つまり、（５０）式の右辺の積分内の分母の基本周波数ηに関する積分は、その基本周波数ηがほぼｘに等しいとき、第２倍音η＋1200ｌｏｇ２２がほぼｘに等しいとき、第３倍音η＋1200ｌｏｇ２３がほぼｘに等しいとき、・・・、第１６倍音η＋1200ｌｏｇ２１６がほぼｘに等しいときの僅か１６×５箇所に関する（５３）式の積分になる。
【００６８】
そして、そのｘを定義域範囲内で３６０通り変化させて積分するため、分母は（５３）式の計算を１６×５×３×３６０回（Ｈ×Ｎａ×Ｍ×Ｎｘ回）繰り返すことになる。これは、
【数１１８】
2006285052.tif 000075
をすべての基本周波数Ｆ（３００通り）と音モデルの数ｍ（３通り）について求める場合にも共通に用いることができるので、上記は一度だけ計算すればよい。一方、（５０）式の右辺の積分内の分子に関しては、あるｘにおける計算に関連している基本周波数Ｆは、その値の範囲の３００箇所より大幅に少なく、１６×５箇所となる。これは、分母のときと同様に考えると、基本周波数Ｆがほぼｘに等しいとき、５箇所のＦに対して分子を計算するだけでよいからである。同様に、第２〜１６倍音Ｆ＋1200ｌｏｇ２ｈがほぼｘに等しいときも分子を計算する必要があるので、合計１６×５回の（５３）式の計算が必要になる。つまり、あるｘにおける分子の計算結果は、８０箇所のＦにだけ影響を与え、残り２２０箇所には影響を与えない。ｍ（３通り）についても求めるので、最終的に、分母と分子のそれぞれで１６×５×３×３６０回（Ｈ×Ｎａ×Ｍ×Ｎｘ回、計８６４００回）、（５３）式の計算を繰り返すことになる。ここで、分子を分母よりも前に計算しておけば、分母はその分子を合計すれば求めることができる。よって、分子と分母を共に求めたとしても、（５３）式の計算は、８６４００回繰り返せばよいことがわかる。この演算回数は、前述の高速化を実施しない場合の演算回数の１／６０の演算回数であり、この程度の演算回数であれば普通に利用できるパーソナルコンピュータでも短時間で演算が可能である。
【００６９】
さらに、（５３）式の計算自体を高速化することも考える。（５７）式の計算に着目すると、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）の差が一定範囲内（ここではｘとＦの離散化幅が２０ｃｅｎｔ(半音の音高差１００ｃｅｎｔの五分の一）でＷが１７ｃｅｎｔのときを前提に±２以内の５回とする)になる場合にだけ計算する場合、
【数１１９】
2006285052.tif 000076
のｙは、離散化して計算するときに常にｙ＝−２＋α、−１＋α、０＋α、１＋α、２＋αの５通りしか取らないことがわかる。ここで、αは、（Ｆ＋1200ｌｏｇ２ｈ）を離散化して表現したときの、０．５以下の小数である。したがって、上記の５通りについて（５９）式を事前に計算して格納しておけば、実際の推定時には、それを読み出して掛け算するだけで等価な計算をおこなうことができ、非常に高速化できる。また、1200ｌｏｇ２ｈも事前に計算して格納しておくとよい。なお、この高速化は、ｘとＦの離散化幅をｄとしたときに、（３Ｗ／ｄ）より小さいもしくは近い正の整数ｂ（上記では２）を求め、Ｎａを（２ｂ＋１）回と設定するように一般化できる。その場合、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）が−ｂ＋α，−ｂ＋１＋α，…，０＋α，…，ｂ−１＋α，ｂ＋αの（２ｂ＋１）通りの値を取る。なお、ここで（３Ｗ／ｄ）の分子の３は３以外の任意の正の整数でもよく、小さいほど演算回数が少なくなる。
【００７０】
次に、（５１）式の計算では、（５０）式と右辺の積分内の分母は共通である。（５１）式の右辺の積分内の分子は、対象となる範囲のＦとｍ、ｈに関して、前述の（５５）式をｘの関数として計算すればよい。前述のとおり、これは、（５２）式から（５６）式を取り除いたものであり、上記で述べた高速化手法によって同様に（５１）式の計算は、高速化できる。
【００７１】
上記の計算の流れを整理すると次のようになる。
１．1200ｌｏｇ２ｈとｅｘｐ［−（ｘ−（Ｆ＋1200ｌｏｇ２ｈ））２／２Ｗ２］を事前に計算してメモリに格納する。
２．以下の処理を収束するまで、もしくは、事前に定めた回数だけ繰り返す。
３．入力音響信号の周波数成分の確率密度関数（（１）式）の各周波数ｘに対して、Ｎｘ回、以下の処理をおこなう（例えば、定義域の範囲を３６０個に離散化していたら３６０回おこなう）。
４．事前に計算した結果を利用して、ｘ−（Ｆ＋1200ｌｏｇ２ｈ）がほぼ０のＦに関して、ｍごとに（５１）式の右辺の積分内の分子
【数１２０】
2006285052.tif 000077
をＭ回求める。また、これから（５０）式の右辺の積分内の分子（（５２）式）も求める。
５．上記の結果を利用して、（５０）式と（５１）式の右辺の積分内の分母を求める。
６．（５０）式と（５１）式の右辺の積分内の分数の値が決まるので、それを、現在のｘの計算に関連している基本周波数Ｆ（３００箇所中１６×５（Ｈ×Ｎａ）箇所）に限定して、（４７）式、（４８）式に加算していく。
【００７２】
こうして、ｘを変化させながら順に加算していくことで、（５０）式と（５１）式の右辺の積分が実現できる。
【００７３】
本発明の方法を実施する図２及び図３に示したアルゴリズムを実施するプログラムをコンピュータで実行することにより、上記の各演算を行う手段がコンピュータ内に実現されて本発明の音高推定装置が構成される。したがって本発明の音高推定装置は、本発明のプログラムをコンピュータで実行した結果物である。
【００７４】
上記のようにして基本周波数の確率密度関数と解釈できる重みω（ｔ）（Ｆ，ｍ）とすべての音モデルの確率密度関数ｐ（ｘ｜Ｆ，ｍ，μ（ｔ）（Ｆ，ｍ））の第ｈ次高調波成分の大きさｃ（ｔ）（ｈ｜Ｆ，ｍ）とをコンピュータを利用した演算により求めることにより、少なくとも従来の６０倍以上の速さで演算を完了することができるので、高速のコンピュータを用いなくてもリアルタイムでの音高推定が可能になる。
【００７５】
なお基本周波数の確率密度関数と解釈できる重みが求められた後の処理は、特許第３４１３６３４号公報に記載のようにマルチエージェントモデルを導入して、確率密度関数の中で所定の基準を満たすピークの軌跡が異なるエージェントを追跡し、信頼度が高くパワーの大きいエージェントが持つ基本周波数の軌跡を採択すればよい。なおこの点は特許第３４１３６３４号公報及び前述の非特許文献１及び２に詳しく説明されているので省略する。
【図面の簡単な説明】
【００７６】
【図１】音モデルのパラメータの推定を説明するために用いる図である。
【図２】本発明のプログラムのアルゴリズムを示すフローチャートである。
【図３】図２のアルゴリズムの一部を詳細に示したフローチャートである。
【図１】
2006285052.tif 000124
【図２】
2006285052.tif 000125
【図３】
2006285052.tif 000126