WO2006106946A1

WO2006106946A1 - 音高推定方法及び装置並びに音高推定用プログラム

Info

Publication number: WO2006106946A1
Application number: PCT/JP2006/306899
Authority: WO
Inventors: Masataka Goto
Original assignee: National Institute of Advanced Industrial Science and Technology AIST
Current assignee: National Institute of Advanced Industrial Science and Technology AIST
Priority date: 2005-04-01
Filing date: 2006-03-31
Publication date: 2006-10-12
Anticipated expiration: 2007-10-01
Also published as: GB2440079B; US20080312913A1; US7885808B2; GB0721502D0; JP2006285052A; JP4517045B2; GB2440079A

Description

明細書

音高推定方法及び装置並びに音高推定用プログラム

技術分野

[0001] 本発明は、混合音中の各構成音 (基本周波数)の音高と音量を推定する音高推定方法及び装置並びに音高推定用プログラムに関するものである。

背景技術

[0002] CD等による実世界の音響信号は、事前に音源数を仮定することが不可能な混合音である。このような混合音中では、周波数成分が頻繁に重複する上に、基本周波数成分が存在しないような音も存在する。しかし、従来の音高推定技術の多くは、少数の音源数を仮定し、周波数成分を局所的に追跡したり、基本周波数成分の存在に依存したり、していた。そのために、前述の実世界の混合音には適用できな力つた。

[0003] そこで、本発明者は、特許第 3413634号公報 (特許文献 1)に示される「音高推定方法及び装置」と題する発明を提案した。この発明においては、入力の混合音には、あらゆる基本周波数 (本願明細書中で抽象的に使用される「音高」に相当するもの）の音が様々な音量で同時に含まれていると考える。そしてこの発明では、統計的手法を利用するために、入力の周波数成分を確率密度関数 (観測した分布)で表現し、各音の高調波構造に対応する確率分布を音モデルとして導入する。そして、周波数成分の確率密度関数が、対象とするあらゆる基本周波数の音モデルの混合分布モデル (重み付き和のモデル)カゝら生成されたと考える。この混合分布中の各音モデルの重みは、各高調波構造が相対的にどれくらい優勢力を表すことから、基本周波数の確率密度関数と呼ぶ (混合分布中にお!、て音モデルが優勢になればなるほど、そのモデルの基本周波数の確率が高くなる)。この重みの値 (すなわち基本周波数の確率密度関数）は、 EM (Expectation - Maximization)ァノレゴリズム（Dempste r, A. P. , Laird, N.M and Rubin、 D. B. ： Maximum likeiinood from mc omplete data via the EM algorithm, J. Roy. Stat. Soc. B, Vol. 39, No . 1, pp. 1— 38 (1977) )を用いることで推定できる。こうして求めた基本周波数の確率密度関数は、混合音中の構成音が、どの音高でどれぐらいの音量で鳴っているかを表している。

[0004] また、本発明者は、この従来の「音高推定方法及び装置」の発明を発展または拡張させる技術を非特許文献 1及び 2の二つの論文に公表している。非特許文献 1は、 2 001年 5月に公表された「A PREDOMINANT- FO ESTIMATION METHOD FOR CD RECORDINGS: MAP ESTIMATION USING EM ALGORITHM FOR ADAPTIVE T ONE MODELSJと題する餘文で、 The 2001 IEEE International Conference on Acous tics, Speech, and Signal Processingの予稿集 Vの 3365- 3368頁に掲載された。また非特許文献 2は、 2004年 9月に公表された「A real-time music-scene-description syst em: predominant— FO estimation for detecting melody and bass lines in real-world au dio signalsjと題する論文で、 Speech Communication 43(2004)の第 311頁〜第 329 頁に掲載された。これら二つの非特許文献で提案された拡張は、音モデルの多重化と、音モデルのパラメータの推定と、モデルパラメータに関する事前分布の導入である。なおこれらの拡張は後に詳しく説明する。

特許文献 1：特許第 3413634号公報

非特許文献 1：「A PREDOMINANT- FO ESTIMATION METHOD FOR CD RECOR DINGS: MAP ESTIMATIONUSING EM ALGORITHM FOR ADAPTIVE TONE MO DELS」と題する餘文 (The2001 IEEE International Conference on Acoustics, Speech, and Signal Processingの予稿集 Vの 3365- 3368頁）

非特許文献 2：「A real-time music-scene-description system: predominant- FOestim ation for detecting melody and bass lines in real-world audio signalsjと題する文で (Speech Communication 43(2004)の第 311頁〜第 329頁）

発明の開示

発明が解決しょうとする課題

[0005] ところで、上記拡張された技術をコンピュータを用いて実現して、基本周波数の確率密度関数の重みと高調波成分の大きさとを推定するためには、演算回数が非常に多くなり、高速演算機能を有するコンピュータを用いなければ、短い時間で推定結果が得られな、と、う問題があった。

[0006] 本発明の目的は、従来よりも少ない演算回数で基本周波数の確率密度関数の重みと高調波成分の大きさとを推定できる音高推定方法及び装置並びに音高推定用プログラムを提供することにある。

課題を解決するための手段

[0007] 本発明の音高推定方法では、次のようにして基本周波数の確率密度関数の重みと高調波成分の大きさとを推定する。

[0008] まず、入力される混合音に含まれる周波数成分を観測し、観測した周波数成分を下記 (a)の対数周波数 X上の確率密度関数として表現する。

[数 1]

[0009] そして非特許文献 1及び 2に開示した技術 (音モデルの多重化、音モデルのパラメータの推定及びモデルパラメータに関する事前分布の導入)を、上記 (a)で表現される観測した周波数成分の確率密度関数から、下記 (b)で表現される基本周波数 Fの確率密度関数を求める過程にぉ、て採用する。

[数 2]

[0010] まず音モデルの多重化では、同一基本周波数に対して M種類の音モデルがあるものとして基本周波数力の m番目の音モデルの確率密度関数を p (x I F, m, ^(t) ( F, m) )と表現する。ただし、 μ ^ω (F, m)は、 m番目の音モデルの高調波成分の大きさの比率を表すモデルパラメータである。

[0011] また音モデルのパラメータの推定では、観測した周波数成分の確率密度関数が、下記 (c)で定義した混合分布モデル p (x I Θ ^(t))から生成されたものと考える。ただし ω ^(t) (F, m)は基本周波数力の m番目の音モデルの重みである。

[数 3] ρ(χ Ίη) ) dF

[0012] なお上記（c)の式において、 θ ^ωは音モデルの重み co^W(F, m)と音モデルの高調波成分の大きさの比率 /^(F, m)を包含したモデルパラメータ 0 ^(t) = ^(t), μ ^ω} であり、 co^(t) = {co^(t)(F, m) I Fl≤F≤Fh, m=l, ···, M}であり、 μ ^(t) = { ^(t) (F, m) I Fl≤F≤Fh, m=l, ···, M}であり、ここで Flは許容される基本周波数の下限であり、 Fhは許容される基本周波数の上限である。

そして上記 (b)の基本周波数 Fの確率密度関数を下記 (d)の解釈により重み ω ^(t) ( F, m)から求める。

画

[0013] さらにモデルパラメータに関する事前分布の導入では、モデルパラメータ Θ ^(t)〖こ関する事前分布に基づいてモデルパラメータ Θ ^(t)の最大事後確率推定量を EM(Expe ctation-Maximization)アルゴリズムを用いて推定する。そして、この推定から事前分布を考慮した上記 (b)の基本周波数 Fの確率密度関数と解釈できる重み ω ^(t) (F, m )と、すべての音モデルの確率密度関数 p (X I F, m, ^ ^(t) (F, m) )の ^(t) (F, m)が表す第 h次高調波成分の大きさ c^W(h I F, m) (h=l, ···, H)とを求めるために用いる下記 (e)及び (f)で表される二つのパラメータ推定値を求めるための式を定める。ただし Hは基本周波数の周波数成分も含めた高調波成分の数である。

[数 5]

丄十 p_t

[数 6]

(f)

[0014] 上記 (e)及び (f)中において、下記 (g)及び (h)は、下記 (i)及び (j)が 0になる無情報事前分布のときの最尤推定値である。

[数 7]

[数 8]

(h)

[数 9]

[数 10]

また上記 (e)及び (f )中にお、て、下記 (k)は重み ω ^(t) (F, m)の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、下記 (1)はモデルノラメータ ^(t) (F, m)の単峰性の事前分布を考えたときに最大値を取る最も起こりやす、パラメータであり、また上記 (i)は下記 (k)をどれぐら!/、重視した事前分布とするかを決めるパラメータであり、上記 (j)は下記 (1)をどれぐら、重視した事前分布とするかを決めるパラメータである。

[数 11]

m (I)

[0016] また上記 (g)及び (h)の式において、 m)及び/ ^W(F, m)は上記（e)及び (f)を反復計算する際の一つ前の古いパラメータ推定値であり、 ηは基本周波数であり、 υは何番目の音モデルであるかを示す。

[0017] 本発明が改良の対象とする音高推定方法では、上記 (e)及び (f )の二つのパラメ一タ推定値を求めるための式を用いた反復計算により、上記 (b)の基本周波数の確率密度関数と解釈できる重み co^(t)(F, m)とすべての音モデルの確率密度関数 p(x I F, m, μ ^(t)(F, m))のモデルパラメータ ^(t) (F, m)が表す第 h次高調波成分の大きさ c^(t)(h I F, m)とをコンピュータを禾 IJ用した演算により求めることにより基本周波数の音高を推定する。

[0018] 本発明では、上記 (e)及び (f)を上記 (g)及び (h)によってコンピュータを利用して演算するために、次のようにする。まず上記 (g)で表現される推定値を示す計算式中の分子を下記 (m)に示す Xの関数として展開する。但し下記 (m)に示す式において aT^(t)(F, m)は古い重みであり、 ^¾(11 I F, m)は古い第 h次高調波成分の大きさの比率であり、 Hは基本周波数の周波数成分を含めた高調波成分の数であり、 mは M種類の音モデルの中の何番目の音モデルかを示すものであり、 Wは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差である。

[数 13] f_{(t)f ir /{t)f l P} 、 1 ( ( ^ ( + 12001og₂ /i))²

(m)

[0019] 上記（m)中の 12001og hと exp [― (x— (F + 12001og h))²/2W²]を事前に計算し

2 2

てコンピュータのメモリに格納する。

[0020] 上記 (e)及び (f)の二つのパラメータ推定値を求めるための式を事前に定めた回数分だけ反復計算するために、上記 (g)及び (h)の式の計算では、観測した周波数成分の確率密度関数の離散化後の各周波数 Xに対して以下の第 1の演算処理を Nx回実行する。但し Nxは Xの定義域の範囲の離散化数である。

[0021] 第 1の演算処理では、 M種類の音モデルについてそれぞれ以下の第 2の演算処理を実行して、上記 (m)式の演算結果を求める。そして上記 (m)式の演算結果を基本周波数 Fと m番目の音モデルに関して積分して上記 (g)及び (h)の式の分母を求め、観測した周波数成分の確率密度関数を上記 (g)及び (h)の式に代入して上記 (g) 及び (h)の演算を行う。

[0022] また第 2の演算処理では、基本周波数を含む高調波成分の数 Hだけ、第 3の演算処理を実行して下記 (n)の式の演算結果を求め、下記 (n)の式の演算結果を H個加算して上記（m)の式の演算結果を求める。

[数 14]

[0023] さらに第 3の演算処理では、 X - (F+ 12001og h)が 0に近い基本周波数 Fに関して

2

第 4の演算処理を Na回実行して上記 (n)の式の演算結果を求める。但し Naは x— ( F + 12001og h)が充分 0に近い範囲の離散化後の Fの個数を表す小さい正の整数

2

である。

[0024] そして第 4の演算処理では、事前にメモリに格納した exp [—（x— (F+ 12001og h)

2

) 2W²]を用いて下記 (o)の式を求める。

[数 15]

[0025] 最後に、上記（o)の式に古、重み ω ' ^(t) (F, m)をかけて上記 (n)の式の演算結果を求める。

[0026] 本発明の方法によれば、事前にメモリに格納した exp [—（X— (F+ 12001og h) ) ²

2

/2W²]を用いることができるので、演算回数を減らすことができる。特に本発明では、第 4の演算処理の回数を Na回と少なくして、上記 (m)の式の演算結果を求めても、演算精度が低下することがな、ことを見、だしたことを根拠として、第 4の演算処理回数を制限する。その結果、従来よりも大幅に演算回数を少なくすることができて、演算時間の短縮ィ匕を可能にした。

[0027] なお対数周波数 Xと基本周波数 Fの離散化幅を dとしたときには、 (3W/d)より小さいもしくは近い正の整数 bを求めて、 Naを（2b + l)回と決定し、離散化して計算するときに X— (F + 12001og h)が一 b+ α, — b + l+ α, ···, 0+ α, ···, b— 1+ α, b

2

+ αの（2b + 1)通りの値を取るようにすればょ、。そしてメモリには、 X- (F+ 12001ο g h)カ b+ a, -b + l+ α, ···, 0+ α , ···, b— 1+ α , b+ αの（2b + l)通りの

2

値を取るときの exp [—（X— (F + 12001og h))²Z2W²]の値を事前に格納しておくの

2

が好まし!/、。ここで前述の wは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差である。またひは (F+12001og h)を離散化して表現したときの 0. 5以

2

下の小数である。なお、ここで（3WZd)の分子の 3は 3以外の任意の正の整数でもよぐ小さいほど演算回数が少なくなる。

[0028] より具体的には、対数周波数 Xと基本周波数 Fの離散化幅が 20cent (半音の音高差 lOOcentの五分の一)で Wが 17centのときは、第 4の演算処理を行う回数 Naは 5 回にすることが好ましい。この場合には、離散化して計算するときに X— (F+12001og h)が一 2+ α, -1+ α, 0+ α, 1+ α, 2+ αの 5通りの値を取る。そして、 αは（

2

F + 12001og h)を離散化して表現したときの 0. 5以下の小数である。このようにする

2

ことによって演算回数を大幅に少なくすることができる。なおこの場合においては、メモリには、 X— (F+12001og h)が— 2+α, — 0+ α, 1+ α, 2+αの値を

2

取るときの exp [—（X— (F+12001og h))²Z2W²]の値を事前に格納しておくのが好

2

ましい。また、 12001og hも事前に計算して格納しておくとよい。これらの格納によって

2

、演算回数を更に低減できる。

[0029] 本発明の音高推定装置では、前述の本発明の音高推定方法をコンピュータを用いて実施する。そのために本発明の音高推定装置は、上記 (g)で表現される推定値を示す計算式中の分子を上記 (m)に示す Xの関数として展開する手段と、上記）中の 12001og hと exp [—（x— (F + 12001og h))²/2W²]を事前に計算してコンビユー

2 2

タのメモリに格納する手段と、前述の第 1の演算処理を実行する第 1の演算処理手段と、前述の第 2の演算処理を実行する第 2の演算処理手段と、前述の第 3の演算処理を実行する第 3の演算処理手段と、前述の第 4の演算処理を実行する第 4の演算処理手段とを備えている。

[0030] また本発明の音高推定用プログラムは、本発明の音高推定方法をコンピュータを用いて実施するためにコンピュータにインストールされるものである。本発明の音高推定用プログラムは、上記 (_g)で表現される推定値を示す計算式中の分子を上記 (m) に示す Xの関数として展開する機能と、上記 (m)中の 12001og hと exp [—（x—（F+ 1

2

2001og h) ) ²Z2W²]を事前に計算してコンピュータのメモリに格納する機能と、前述

2

の第 1の演算処理を実行する機能と、前述の第 2の演算処理を実行する機能と、前述の第 3の演算処理を実行する機能と、前述の第 4の演算処理を実行する機能をコンピュータ内に実現するように構成されて、る。

発明の効果

[0031] 本発明によれば、音源数を仮定せず、周波数成分の局所的な追跡も行わず、また基本周波数成分の存在を前提とせずに、音高推定する際に、従来よりも大幅に演算回数を少なくすることができて、演算時間を短縮ィ匕することができる。

図面の簡単な説明

[0032] [図 1]図 1は音モデルのパラメータの推定を説明するために用いる図である。

[図 2]図 2は本発明のプログラムのアルゴリズムを示すフローチャートである。

[図 3]図 3は図 2のアルゴリズムの一部を詳細に示したフローチャートである。

発明を実施するための最良の形態

[0033] 以下図面を参照しながら本発明の音高推定方法及びプログラムの実施の形態の一例を詳細に説明する。まず本発明の方法の実施の形態の一例について説明する前提として、特許第 3413634号の発明を拡張した前述の非特許文献 1及び 2に提案された公知の三つの拡張方法につ!、て簡単に説明する。

[0034] [拡張方法 1]音モデルを多重化

特許第 3413634号公報に記載の発明では、同一基本周波数には一つの音モデルしか用意していな力つた。し力しながら実際には、ある基本周波数に、異なる高調波構造を持つ音が入れ替わり立ち替わり現れることがある。そこで、同一基本周波数に対して複数の音モデルを用意し、それらの混合分布でモデル化する。具体的な方法は、後に詳しく説明する。 [0035] [拡張方法 2]音モデルのパラメータを推定

特許第 3413634号公報に記載の従来の音モデルでは、各高調波成分の大きさの比率を固定して、た (ある理想的な音モデルを仮定して、た)。し力しながらこれは実世界の混合音中の高調波構造とは必ずしも一致しておらず、精度向上のためには、さらに改善をする余地が残されていた。そこで拡張方法 2では、音モデルの高調波成分の大きさの比率もモデルパラメータに力卩ぇ音モデルのパラメータを各時刻において EMアルゴリズムで推定することとした。具体的な方法は後に説明する。

[0036] [拡張方法 3]モデルパラメータに関する事前分布を導入

特許第 3413634号公報に記載の従来の方法では、音モデルの重み (基本周波数の確率密度関数）に関する事前知識は仮定していな力つた。しかし本発明を様々な用途に適用していく上で、たとえ事前に基本周波数がどの周波数の近傍にあるかを与えてでも、より誤検出の少な、基本周波数を求めた、と、うような応用が考えられる。例えば、演奏分析やビブラート分析等の目的では、楽曲をヘッドホン力も聴取しながらの歌唱や楽器演奏によって、各時刻におけるおおよその基本周波数を事前知識として用意しておき、実際の楽曲中のより正確な基本周波数を得ることが求められている。そこで、従来のモデルパラメータの最尤推定の枠組みを拡張し、モデルパラメータに関する事前分布に基づヽて最大事後確率推定 (MAP推定： Maximum A Posteriori Probability Estimation)を行う。その際、 [拡張方法 2]でモデルパラメータに加えた、音モデルの高調波成分の大きさの比率に関する事前分布も導入する。具体的な方法は後に説明する。

[0037] 拡張方法 1〜3について、式を用いて更に具体的に説明する。まず入力される混合音 (入力音響信号)に含まれる観測した周波数成分の確率密度関数を下記（1)で表現する。

[数 16]

[0038] そして上記（1)式の周波数成分の確率密度関数から、下記（2)式で表現される基本周波数 Fの確率密度関数を求める過程における、具体的な拡張方法を述べる。 [数 17]

[0039] 上記（1)式の観測した周波数成分の確率密度関数は、混合音 (入力音響信号)から例えばマルチレートフィルタバンク（Vetterli,M.: A Theory of Multirate Filter Banks , IEEE Trans, on ASSP, Vol.ASSP- 35, No.3, pp.356- 372 (1987)参照）によって求めることができる。なおこのマルチレートフィルタバンクについては、特許第 3413634号公報の図 2及び前述の非特許文献 2に示された Fig . 3にバイナリ一ツリー状のフィルタバンクの構成の一例とその詳細が説明されている。ここで、（1)式及び（2)式中の t はフレームシフト（10msec)を時間単位とする時刻であり、 xと Fは centの単位で表された対数スケールの対数周波数と基本周波数である。なお、 Hzで表された周波数 f

H

は、下記（3)式により centで表された周波数 f に変換されるものとする。

z cent

[数 18] f_cen. — 一一- ₍₃₎

■/ cent = 1200 log &o 2 _{440 χ 2} 5

[0040] そこで、前述の [拡張方法 1]と [拡張方法 2]を実現するために、同一基本周波数に対して Μ種類の音モデルがあるものとし、基本周波数が Fの m番目の音モデルの確率密度関数 p (X I F、 m、 ^ω (Fゝ m) )にモデルパラメータ μ ^(t) (Fゝ m)を導入する。

[0041] なお以下に説明する（4)式から（51)式までについては、前述の非特許文献 1にお V、て（2)式〜（36)式としてすでに公表されて、るので参照された！、。

[0042] 基本周波数が Fの m番目の音モデルの確率密度関数 p (x | F、 m、 ^(t) (F、 m) ) は、次式のように表されるものとする。

[数 19]

H

p(a;| F, m， i ( (F， ) ) - ^ ( ， | F， τη, μ ) (F， m) )

-- - --… (4)

[数 20]

—-- -(5)

[数 21]

μ (り（F，m)

1，.·.，Η}

"(6)

[数 22]

上記 (4)〜（7)の式は、基本周波数力のときに、その高調波成分がどの周波数にどれくらい現れるかをモデル化したものである（図 1)。上記式においては、 Ηは基本周波数 Fの周波数成分も含めた高調波成分の数、 Wはガウス分布 G(x;x、 σ)の標

0

準偏差を表している。また c(t) (h I F、 m)は、第 h次高調波成分の大きさを表し、次式を満たすものとする。

[数 23]

そして、上記（1)式で表される観測した周波数成分の確率密度関数が、次式で定義されるような、 p(x I F、 m、 ^ω (Fゝ m))の混合分布モデル p(x | θ ^ω)から生成されたと考える。

[数 24]

"(9) [数 25]

b O "(10) l r

fv (F, m) I Fl < < Fh₇m = 1,...,M}

-(11)

[数 27] 1/ μ (

(12)

[0045] 上記（11)及び（12)式において、 Fhと Flは、許容される基本周波数の上限と下限でであありり、、 w^(t) (F、 m)は、次式を満たすような、音モデルの重みである。

[数 28]

[0046] 実世界の混合音に対して事前に音源数を仮定することは不可能なため、上記 (6) 式のように、あらゆる基本周波数の可能性を同時に考慮してモデルィ匕することが重要となる。最終的に、モデル p(x I Θ ^(t))から、観測した確率密度関数 [上記（1)式]が生成されたかのようにモデルパラメータ Θ ^(t)を推定できれば、その重み w^(t) (F、 m)は各高調波構造が相対的にどれくらい優勢力を表すことになる。そのため、次式のように基本周波数 Fの確率密度関数を解釈することができる。

[数 29]

M

p 、()(F) = w (F, m) (Fl < < Fh) -(14)

[0047] 次に、前述の [拡張方法 3]の事前分布の導入を行う。 [拡張方法 3]を実現するために、 0 (t)の事前分布 p (0 ^(t))は、下記の式（19)のように下記式（20)と下記式（

Oi

21)の積で与えられる。下記の式（19)〜（21)に示される p (co^(t))と p (/x ^(t))は、最

Oi Oi

も起りやすいパラメータを

[数 30]

w^⁽ (F,m) -……〖¹⁵) と

[数 31]

Mo ^m) …-- - としたとき（ただし、式（16)は

[数 32]

—(〗7) である)に、そこで最大値を取るような単峰性の事前分布である。ただし、 Z 、 Z は正 β 規化係数であり、

[数 33]

β^ , , --- (₁₈, の二つは、最大値をどれくらい重視した事前分布とするかを決めるパラメータであり、 0のときに無情報事前分布 (一様分布)となる。

[数 34]

Poi(0^(t)) =Ροί(^^(ί)) ροΐ( ^(ί)) ——— (¹⁹)

[数 35]

P i(w^(t)) = exp(— ） D_w(w^;w^(t)))

w

(20) [数 36] ρ_0ί(μ^(ί)) = - (- )

μ ^j广 F

F^l

- -…- (21) 上記（20)式中の下記の式

[数 37]

D_w(u^; …― と上記（21)式中の下記の式

[数 38]

は、次のような K— L情報量（Kullback— Leibler's information)である。

[数 39]

D dF

—— f24)

[数 40]

—— 25) 以上の説明から、上記（1)式の確率密度関数を観測したときに、そのモデル p (X I Θ ^(t))のパラメータ θ ^(t)を、事前分布 p ( Θ ^(t))に基づいて推定する問題を解けばよい

Oi

ことがわかる。この事前分布 p ( Θ ^(t))に基づくパラメータ θ ^(t)の最大事後確率推定

Oi

量 (MAP推定量）は、次式を最大化することで得られる。

[数 41] )十 log p„ dx

-…… (26)

[0049] し力しながらこの最大化問題は解析的に解くことが困難なため、 EMアルゴリズム ( Dempster, A. P. , Laird, N.M and Rubin、 D. B. ： Maximum likelihood f rom incomplete data via the EM algorithm, J. Roy. Stat. Soc. B, Vol . 39, No. 1, pp. 1 38 (1977) )を用! /、て 0 ^ωを推定する。 ΕΜァノレゴ!;ズムは、不完全な観測データ力最尤推定を行うために用いられることが多いが、最大事後確率推定の場合にも適用できる。最尤推定では、平均対数尤度の条件付き期待値を求める Εステップ（expectation step)とその最大化をおこなう Mステップ（maximiz ation step)を交互に繰り返す。し力最大事後確率推定の場合には、条件付き期待値に事前分布の対数を加えたものの最大化を繰り返す。ここでは各繰り返しにお V、て、古、パラメータ推定値 0 ' ^(t) = {w' ^(t\ μ ' ^(t) }を更新して新、パラメータ推定値 (下記式 (27) )

を求めていく。

[数 42] —— (27)

[0050] 対数周波数 Xにおヽて観測した各周波数成分が、どの基本周波数のどの音モデルのどの倍音力生成されたのかを表す隠れ変数 F、 m、 hを導入して、 EMァルゴリズムを以下のように定式ィ匕することができる。

[0051] (Eステップ）

最尤推定の場合には、平均対数尤度の条件付き期待値 <3 ( θ ^ω I θ ^))を求めるが、最大事後確率推定の場合には、それに log p ( θ ^ω

Oi ^ω)を加えた Q ( θ

MAP I 9 ^(t

))を求める。

[数 43]

…… 128)

[数 44] p^{( >} (x)E_F,_n,_h[logp(x_J F, m,

dx

-oo

(29)

[0052] 上記式において、条件付き期待値 E [a I b]は、条件 bにより決定される確率分

F、 m、 h

布を持つ隠れ変数 F、 m、 hに関する、 aの期待値を意味する。

[0053] (Mステップ）

Q Θ

MAP ( Θ ^(t) I θ ' ^ω)を ^(t)の関数として最大化して、更新後の新 U、推定値

[数 45]

Θ —— (30) を下記の式で得る。

[数 46]

W) - argmax

^(ί)) ---OD

、て、上記（29)式は下記のようになる。

[数 47]

(32) 上記式中の完全データの対数尤度は

[数 48]

—— -(33) で与えられる。また、 log ρ ( 0 ^(t))は、 [数 49]

-

(34) となる。

次に、 Mステップに関しては、上記（31)式が、上記（8)式と上記（13)式を条件とする条件付き変分問題となっている。この問題は、 Lagrangeの乗数 w、 λ μを導入し、次の Euler - Lagrangeの微分方程式によって解くことができる。

[数 50]

(t)_{f r} 、

= 0

M(Fh-Fl) '

•(35)

[数 51]

― ¾)(ί ) 0

"(36) これより、

[数 52]

m)

(37)

[数 53] dx + ¾)(F，m)4? ( |F,m)

(38) が得られる。これらの式において、 Lagrangeの乗数は（8)式、（13)式から

[数 54]

\ ― 1 I り (39)

八 ― 丄 I wi

[数 55]

140) と定まり、 p (F、 m、 h I χ、 θ ' ^(t))、 p (F、 m I x、 Θ ' ^(t))はベイズの定理から、

[数 56]

141 )

[数 57]

42) となる。以上から、新しいパラメータ推定値

[数 58]

と

[数 59]

を求める式は次のようになる。

[数 60] w(t (F, rn) = (45)

1十 3^{T)

[数 61]

一-—一-… (46) 式中の

[数 62]

と

[数 63]

は、

[数 64] 0 ---(49)

の無情報事前分布のとき、つまり最尤推定の場合の推定値である。

[数 65]

"(50)

[数 66] ,

--——- (51 )

[0055] これらの反復計算により、事前分布を考慮した上記（2)式の基本周波数の確率密度関数が、上記（14)式によって重み w^(t) (F、 m)から求まる。さらに、すべての音モデル p (x I Fゝ m、 ^ω (F、 m) )の各高調波成分の大きさの比率 c^(t) (h | F、 m)も求まる。これにより [拡張方法 1]〜 [拡張方法 3]が実現される。

[0056] 上記のように拡張した音高推定手法をコンピュータを用いて実行するためには、上記 (45)式と上記 (46)式の反復計算が必要となる。しかし、これらの式の反復計算では、上記（50)式と上記（51)式の計算量が多いため、計算能力の限られた (計算速度が遅い)コンピュータでは、この式をそのまま計算しょうとすると計算時間が非常に長くなつてしまうという問題が生じる。

[0057] 計算時間が非常に長くなる理由を説明する。最初に、上記 (50)式を用いて計算結果を普通に求めるときには、どのような計算が必要かを説明する。まず、上記（50)式の計算では、対象となる範囲の Fと mに関して、上記（50)式の右辺の積分内の分子 [数 67]

Η

二 ιν'^ (F, m)∑ p(x, h\F, m,„ m)) G(x; F十 1200 log₂ h, W) — _exp(- 1 土

(52) を xの関数として計算する (上記 (4)式〜（7)式により展開する)。ここで、説明のための計算例として、 Xの定義域の範囲を 360個（N個）に離散化し、 F1から Fhまでの範

X

囲の Fを 300 (N個）に離散化して計算するものと仮定する。また、音モデルの数 M

F

は 3、高調波成分の数 Hは 16とする。このとき、上記（52)式を計算するには、下記の (53)式の計算を 16回繰り返すことになる。

[数 68] t_{t) ( — （F ₊ 1200 1og₂ /i)

C exiH———— ———リ

--— -(53)

[0058] 上記（50)式の右辺の積分内の分子を求めるためには、ある Xに関する上記（52) 式の計算を一回行う。そして上記（50)式の右辺の積分内の分母を求めるためには、 Fと mに関して、 300 X 3回(N X M回）、（52)式の計算を繰り返す必要がある。

F

[0059] さらに、その Xを定義域範囲内で 360通り変化させて積分するため、

[数 69]

を求めるには、分母は 16 X (300 X 3) X 360回、分子は 16 X 360回、上記（53)式の計算を繰り返す必要がある。分母は Fと mを変化させても共通なので再度計算する必要はな、が、分子はすべての F (300通り）と m (3通り）につ!/、て求める必要がある。そのため、最終的に、分母と分子のそれぞれで 16 X (300 X 3) X 360回（H X N

F

X M X N回、計 5184000回）、上記（53)式の計算を繰り返すことになる。ここで、

X

分子を分母よりも前に計算しておけば、分母はその分子を合計すれば求めることができる。よって、分子と分母を共に求めたとしても、上記（53)式の計算は 5184000 回繰り返すことになる。

[0060] そこで本発明は、この計算時間を以下のようにして大幅に削減して高速ィ匕する。以下本発明の方法により、上記の通常計算方法を高速化した高速計算方法を図 2及び図 3に示した本発明のプログラムのアルゴリズムを示すフローチャートを参照しながら説明する。まず、上記（50)式の計算では、対象となる範囲の Fと mに関して、上記（5 0)式の右辺の積分内の分子を、 Xの関数として上記（52)式により計算する。

[0061] 図 2に示すように、上記（52)式中の 12001og hと exp [—（x— (F+ 12001og h) )

2 2

2W²]を事前に計算してコンピュータのメモリに格納する。そして上記 (45)式及び (4 6)式の二つのパラメータ推定値を求めるための式を事前に定めた回数分だけ (もしくは収束するまで)反復計算するために、上記（50)式及び（51)式の計算では、図 3に示すように、上記 (47)式及び (48)式を 0で初期化した後、観測した周波数成分の確率密度関数の各対数周波数 Xに対して以下の第 1の演算処理を Nx回実行する。伹し Nxは Xの定義域の範囲の離散化数である。

[0062] 第 1の演算処理では、 M種類の音モデルについてそれぞれ以下の第 2の演算処理を実行して、上記（52)式の演算結果を求める。そして上記（52)式の演算結果を基本周波数 Fと m番目の音モデルに関して積分して上記（50)式及び（51)式の分母を求め、観測した周波数成分の確率密度関数を上記（50)式及び (51)式に代入して上記（50)式及び（51)式の演算を行う。

[0063] 第 2の演算処理では、基本周波数の周波数成分を含む高調波成分の数 Hだけ、第 3の演算処理を実行して後に説明する下記の（55)式の演算結果を求める。そして（ 55)式の演算結果を H個加算して上記（52)式の演算結果を求める。

[数 70] w'(t、（F, m) p(x, h\ F, m， μ' （F, m))

- w'(t) (F, m (り (h\ F, ?n) G(x F十 1200 log₂ h, W)

._{m /} 、 _{/M /} , 、 1 , r — (F + 1200 1og₉ /i) ) ²、 = w m)c (h\ F, m) exp( —— ―—二—一.）

- --… (55)

[0064] 上記（55)式は、（51)式の右辺の積分内の分子を、対象となる範囲の Fと m、 hに関して Xの関数として計算するものである。 (55)式は（52)式から

[数 71]

、H

∑^丄 — … を取り除いたものである。

[0065] 前述の第 3の演算処理では、 X— (F+ 12001og h)が 0に近い基本周波数 Fに関し

2

て第 4の演算処理を Na回実行して上記（55)式の演算結果を求める。但し本発明では、 Naは X— (F+ 12001og h)が充分 0に近い範囲の Fの個数を表す小さい正の整数としている。後に説明するように好ましくは、対数周波数 Xと基本周波数 Fの離散化幅 dが 20cent (半音の音高差 lOOcentの五分の一)で、前述のガウス分布の標準偏差 Wが 17centのときは、この Naは 5である。

[0066] 第 4の演算処理では、事前にメモリに格納した exp [—（x— (F+ 12001og h) ) ²/2

2

W²]を用いて上記（53)式の演算を行う。そして上記（53)式に古い重み α ^ω (Ρ, m )をかけて上記（55)式の演算結果を求める。このようにして本発明では、音高を推定する。

[0067] 上記をより具体的な例を用いて説明する。

まず上記（52)式中の

[数 72] e_xp (― ί£ ί£± ? ¾ ) -… -… ）の演算を、 Xと (F+ 12001og h)の差が大きくなると急速に 0に近づくことを利用し、そ

2

の差が一定範囲内 (本実施例では、 Xと Fの離散化幅が 20cent (半音の音高差 100c entの五分の一）で Wが 17centのときは、例えば ± 2以内の 5回（ = Na)とする)だけ計算する。

[0068] ここで、ある対数周波数 Xに関して、（50)式の右辺の積分内の分母を計算することを考える。上記の計算範囲の制限により、 (F+ 12001og h)の近傍の Xでのみ上記（5

2

7)式を計算し、それ以外では 0とみなして計算しないことにする。このようにすると、ある対数周波数 Xを起点に考えると、（50)式の右辺の積分内の分母を求めるのに、 (5 3)式の計算を 16 X 300 X 3回繰り返す必要はなぐ 16 X 5 X 3回（H X Na X M回）繰り返すだけでよいことになる。つまり、（50)式の右辺の積分内の分母の基本周波数 7?に関する積分は、その基本周波数 r?がほぼ Xに等しいとき、第 2倍音 7? + 12001 og 2がほぼ Xに等しいとき、第 3倍音 7? + 12001og 3がほぼ Xに等しいとき、…、第 1

2 2

6倍音 7? + 12001og 16がほぼ Xに等しいときの僅か 16 X 5箇所に関する（53)式の積

2

分になる。

[0069] そして、その Xを定義域範囲内で 360通り変化させて積分するため、分母は（53)式の計算を 16 X 5 X 3 X 360回（H X Na X M X Nx回）繰り返すことになる。これは、 [数 73] wl⁽ (F₇ m) —一- をすベての基本周波数 F (300通り）と音モデルの数 m (3通り）につ!/、て求める場合にも共通に用いることができるので、上記は一度だけ計算すればよい。一方、（50) 式の右辺の積分内の分子に関しては、ある Xにおける計算に関連している基本周波数 Fは、その値の範囲の 300箇所より大幅に少なぐ 16 X 5箇所となる。これは、分母のときと同様に考えると、基本周波数 Fがほぼ Xに等しいとき、 5箇所の Fに対して分子を計算するだけでよいからである。同様に、第 2〜16倍音 F + 12001og hがほぼ x

2 に等しいときも分子を計算する必要があるので、合計 16 X 5回の（53)式の計算が必要になる。つまり、ある Xにおける分子の計算結果は、 80箇所の Fにだけ影響を与え、残り 220箇所には影響を与えない。 m (3通り）についても求めるので、最終的に、分母と分子のそれぞれで 16 X 5 X 3 X 360回（H X Na X M X Nx回、計 86400回）、 (53)式の計算を繰り返すことになる。ここで、分子を分母よりも前に計算しておけば、分母はその分子を合計すれば求めることができる。よって、分子と分母を共に求めたとしても、（53)式の計算は、 86400回繰り返せばよいことがわかる。この演算回数は、前述の高速ィ匕を実施しない場合の演算回数の 1Z60の演算回数であり、この程度の演算回数であれば普通に利用できるパーソナルコンピュータでも短時間で演算が可能である。

さらに、（53)式の計算自体を高速ィ匕することも考える。（57)式の計算に着目すると、 X- (F + 12001og h)の差が一定範囲内（ここでは Xと Fの離散化幅が 20cent (半音

2

の音高差 lOOcentの五分の一）で Wが 17centのときを前提に ± 2以内の 5回とする) になる場合にだけ計算する場合、

[数 74] exp(- —— (⁵⁹)

の yは、離散化して計算するときに常に y=— 2+ α、 - 1 + 0+ 1 + 2 + αの 5通りしか取らないことがわかる。ここで、 αは、（F+ 12001og h)を離散化して表現したときの、 0. 5以下の小数である。したがって、上記の 5通りについて（59)式を事前に計算して格納しておけば、実際の推定時には、それを読み出して掛け算するだけで等価な計算をおこなうことができ、非常に高速ィ匕できる。また、 12001og hも事

2 前に計算して格納しておくとよい。なお、この高速ィ匕は、 Xと Fの離散化幅を dとしたときに、（3WZd)より小さいもしくは近い正の整数 b (上記では 2)を求め、 Naを（2b+l )回と設定するように一般ィ匕できる。その場合、 x-(F+12001og h)が— b+ α, — b

2

+1+ α, ···, 0+ α, ···, b-l+ α, b+ αの（2b+ 1)通りの値を取る。なお、ここで（3WZd)の分子の 3は 3以外の任意の正の整数でもよぐ小さいほど演算回数が少なくなる。

[0071] 次に、（51)式の計算では、（50)式と右辺の積分内の分母は共通である。（51)式の右辺の積分内の分子は、対象となる範囲の Fと m、 hに関して、前述の（55)式を X の関数として計算すればよい。前述のとおり、これは、（52)式から（56)式を取り除いたものであり、上記で述べた高速ィヒ手法によって同様に（51)式の計算は、高速化できる。

[0072] 上記の計算の流れを整理すると次のようになる。

1. 12001og hと exp[— (X— (F + 12001og h))²/2W²]を事前に計算してメモリに

2 2

格納する。

2.以下の処理を収束するまで、もしくは、事前に定めた回数だけ繰り返す。

3.入力音響信号の周波数成分の確率密度関数（（1)式)の各周波数 Xに対して、 Nx回、以下の処理をおこなう（例えば、定義域の範囲を 360個に離散化していたら 3 60回おこなう）。

4.事前に計算した結果を利用して、 X— (F + 12001og h)がほぼ 0の Fに関して、 m

2

ごとに（51)式の右辺の積分内の分子

[数 75]

- --一 (60) を M回求める。また、これから（50)式の右辺の積分内の分子（（52)式)も求める。 5.上記の結果を利用して、（50)式と（51)式の右辺の積分内の分母を求める。

6. (50)式と（51)式の右辺の積分内の分数の値が決まるので、それを、現在の Xの計算に関連して、る基本周波数 F (300箇所中 16 X 5 (H X Na)箇所）に限定して、（ 47)式、（48)式に加算していく。

[0073] こうして、 Xを変化させながら順に加算していくことで、（50)式と（51)式の右辺の積分が実現できる。

[0074] 本発明の方法を実施する図 2及び図 3に示したアルゴリズムを実施するプログラムをコンピュータで実行することにより、上記の各演算を行う手段がコンピュータ内に実現されて本発明の音高推定装置が構成される。したがって本発明の音高推定装置は、本発明のプログラムをコンピュータで実行した結果物である。

[0075] 上記のようにして基本周波数の確率密度関数と解釈できる重み ω ^(t) (F, m)とすべての音モデルの確率密度関数 p (x I F, m, ^ ^(t) (F, m) )の第 h次高調波成分の大きさ c^(t) (h I F, m)とをコンピュータを利用した演算により求めることにより、少なくとも従来の 60倍以上の速さで演算を完了することができるので、高速のコンピュータを用いなくてもリアルタイムでの音高推定が可能になる。

[0076] なお基本周波数の確率密度関数と解釈できる重みが求められた後の処理は、特許第 3413634号公報に記載のようにマルチエージェントモデルを導入して、確率密度関数の中で所定の基準を満たすピークの軌跡が異なるエージェントを追跡し、信頼度が高くパワーの大き!/、エージェントが持つ基本周波数の軌跡を採択すればょ、。なおこの点は特許第 3413634号公報及び前述の非特許文献 1及び 2に詳しく説明されているので省略する。

Claims

請求の範囲

入力される混合音に含まれる周波数成分を観測し、観測した周波数成分を下記 (a )の対数周波数 X上の確率密度関数として表現するものとし、

[数 76] ） …- ）前記観測した周波数成分の確率密度関数から、下記 (b)で表現される基本周波数 Fの確率密度関数を求める過程において、音モデルの多重化、音モデルのパラメ一タの推定及びモデルパラメータに関する事前分布の導入を採用し、

[数 77]

前記音モデルの多重化では、同一基本周波数に対して M種類の音モデルがあるものとして前記基本周波数カ^の m番目の音モデルの確率密度関数を p(x I F, m, μ ^(t) (F, m) )と表現し、ただし、 μ ^(t) (F, m)は前記 m番目の音モデルの高調波成分の大きさの比率を表すモデルパラメータであり、

前記音モデルのパラメータの推定では、前記観測した周波数成分の確率密度関数が、下記 (c)で定義した混合分布モデル p(x I Θ ^(t))から生成されたものと考え、ただし0> ^ω (F, m)は基本周波数が Fの m番目の音モデルの重みであり、

[数 78]

(c) ただし上記 (c)にお、て、 θ ^ωは前記音モデルの重み ω ^(t) (F, m)と前記音モデノレの高調波成分の大きさの比率 μ ^(t) (F, m)を包含したモデルパラメータ Θ ^W = { ω ^(t) , (^t)}であり、 ω^ω = {ω^ω(Ρ, m) | Fl≤F≤Fh, m=l, ···, M}であり、 μ ^ω = {μ ⁽' )(F, m) I Fl≤F≤Fh, m=l, ···, M}であり、ここで Flは許容される基本周波数の下限であり、 Fhは許容される基本周波数の上限であり、前記基本周波数 Fの確率密度関数を下記 (d)の解釈により前記重み ω (t) (F, m) から求め、

[数 79]

M

p ₀(F) = w^(F m) (Fl < F < Fh) (d) m=l

前記モデルパラメータに関する事前分布の導入では、前記モデルパラメータ Θ ^(t)に関する事前分布に基づいて前記モデルパラメータ Θ ^(t)の最大事後確率推定量を E M (Expectation-Maximization)アルゴリズムを用いて推定し、この推定から事前分布を考慮した上記 (b)の基本周波数 Fの確率密度関数と解釈できる前記重み ω ^(t) (F, m)とすべての音モデルの確率密度関数 p (X I F, m, ^ ^(t) (F, m) )の ^(t) (F, m)が表す第 h次高調波成分の大きさ c^W (h I F, m) (h= l, · ··, H)とを求めるために用いる下記 (e)及び (f)で表される二つのパラメータ推定値を求めるための式を定め、ただし Hは前記基本周波数の周波数成分も含めた高調波成分の数であり、

[数 80]

[数 81] w^(F, m) 1( I F, m) + ¾) (F, m)c (h\ F, m

«S ( f ( ! ™ m、)丄十 ^ ( ( m)

-- --— (f) また上記 (e)及び (f)中にお、て、下記 (g)及び (h)は、下記 (i)及び (j)が 0になる無情報事前分布のときの最尤推定値であり、

[数 82] ' MlA , ^m)一 / ΡΦ "

-- (g) [数 83]

( άη

…… -- Ν

[数 84] β^) (i)

ド Wl

[数 85] β ( 71、 —一 - < ) 上記 (e)及び (f)中にお、て、下記 (k)は前記重み ω ^(t) (F, m)の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、下記 (1)は前記モデルパラメータ ^(t) (F, m)の単峰性の事前分布を考えたときに最大値を取る最も起こりやす、パラメータであり、また上記 (i)は下記 (k)をどれぐら!/、重視した事前分布とするかを決めるパラメータであり、上記 (j)は下記 (1)をどれぐら、重視した事前分布とするかを決めるパラメータであり、

[数 86] 、 — -- W

[数 87] c^⁽ (h\ F m) (') また上記 (g)及び (h)にお、て、 ω ' ^(t) (F, m)及び μ ' ^(t) (F, m)は上記 (e)及び (f )を反復計算する際の一つ前の古いパラメータ推定値であり、 ηは基本周波数であり、 Vは何番目の音モデルであるかを示し、

上記 (e)及び (f)の二つのパラメータ推定値を求めるための式を用いた反復計算により、上記 (b)の基本周波数の確率密度関数と解釈できる前記重み co ^(t) (F, m)とすベての音モデルの確率密度関数 p (x I F, m, μ ^(t) (F, m) )の前記モデルパラメータ μ ^ω (F, m)が表す第 h次高調波成分の大きさ c^(t) (h | F, m)とをコンピュータを利用した演算により求めることにより基本周波数の音高を推定する音高推定方法であって上記 (e)及び (f)を上記 (g)及び (h)によって前記コンピュータを利用して演算するために、上記 (g)で表現される推定値を示す計算式中の分子を下記 (m)に示す Xの関数として展開し、但し下記 (m)に示す式において aT ^(t) (F, m)は古い重みであり、 c' ^(t) (h | F, m)は古い第 h次高調波成分の大きさの比率であり、 Hは前記基本周波数の周波数成分も含めた高調波成分の数であり、 mは M種類の音モデルの中の何番目の音モデルかを示すものであり、 Wは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、

[数 88]

„ ) -

上記（m)中の 12001og hと exp [― (x— (F + 12001og h) ) ²/2W²]を事前に計算し

2 2

て前記コンピュータのメモリに格納し、

上記 (e)及び (f)の二つのパラメータ推定値を求めるための式を事前に定めた回数分だけ反復計算するために、上記 (g)及び (h)の式の計算では、前記観測した周波数成分の確率密度関数の離散化後の各周波数 Xに対して以下の第 1の演算処理を Nx回実行し、但し Nxは Xの定義域の範囲の離散化数であり、

前記第 1の演算処理では、 M種類の音モデルについてそれぞれ以下の第 2の演算処理を実行して、上記 (m)式の演算結果を求め、上記 (m)式の演算結果を基本周波数 Fと m番目の音モデルに関して積分して上記 (g)及び (h)の式の分母を求め前記観測した周波数成分の確率密度関数を上記 (g)及び (h)の式に代入して上記 (g) 及び (h)の演算を行い、

前記第 2の演算処理では、基本周波数を含む高調波成分の数 Hだけ、第 3の演算処理を実行して下記 (n)の式の演算結果を求め、下記 (n)の式の演算結果を H個加算して上記（m)の式の演算結果を求め、 [数 89]

前記第 3の演算処理では、 X- (F+12001og h)が 0に近い基本周波数 Fに関して

2

第 4の演算処理を Na回実行して上記 (n)の式の演算結果を求め、但し Naは x—（F + 12001og h)が充分 0に近、範囲の離散化後の Fの個数を表す小さ!/、正の整数で

2

あり、

前記第 4の演算処理では、事前に前記メモリに格納した exp [— (x- (F+12001og h) ) ²/2W²]を用いて下記（o)の式を求め、

[数 90]

上記（o)の式に古、重み ω ' ^(t) (F, m)をかけて上記 (n)の式の演算結果を求めることを特徴とする音高推定方法。

[2] 前記対数周波数 Xと前記基本周波数 Fの離散化幅を dとしたときに、（3WZd)より小さいもしくは近い正の整数 bを求めて、前記 Naを（2b +1)回と決定し、離散化して計算するときに X— (F+12001og h)が— b+ α, — b+l+ α, ···, 0+ α, ···, b— 1

2

+ a , b+ aの（2b + 1)通りの値を取り、ここで Wは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、 aは (F+12001og h)を離散化して表現

2

したときの 0. 5以下の小数である請求項 1に記載の音高推定方法。

[3] 前記対数周波数 Xと前記基本周波数 Fの離散化幅を dとしたときに、（3WZd)より小さいもしくは近い正の整数 bを求めて、前記 Naを（2b +1)回と決定し、離散化して計算するときに X— (F+12001og h)が— b+ α, — b+l+ α, ···, 0+ α, ···, b— 1

2

+ a , b+ aの（2b + 1)通りの値を取り、ここで Wは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、 aは (F+12001og h)を離散化して表現したときの 0. 5以下の小数であり、

前記メモリには、 X— (F+12001og h)が— b+ α, — b + 1+ α, ···, 0+ α, ···, b

2

-1+ a, b+ αの（2b + l)通りの値を取るときの exp[— (x— (F + 12001og h))

2

2W²]の値が事前に格納されて!、る請求項 1に記載の音高推定方法。

[4] 前記対数周波数 Xと前記基本周波数 Fの離散化幅が 20centで前記 Wが 17cent のときに、前記 Naを 5と定め、離散化して計算するときに X— (F + 12001og h)がー 2

2

+ α，一 1+ α， 0+ α， 1+ α， 2+ αの値を取り、ここで αは（F+12001og h)を離

2 散化して表現したときの 0. 5以下の小数である請求項 1に記載の音高推定方法。

[5] 前記対数周波数 Xと前記基本周波数 Fの離散化幅が 20centで前記 Wが 17cent のときに、前記 Naを 5と定め、離散化して計算するときに X— (F + 12001og h)がー 2

2

2 散化して表現したときの 0. 5以下の小数であり、

前記メモリに ίま、 X— (F+12001og h)力

2 ー 2+α, — 0+ α, 1+ α, 2+ a の値を取るときの exp [—（X— (F+12001og h))²Z2W²]の値が事前に格納されて

2

V、る請求項 1に記載の音高推定方法。

[6] 入力される混合音に含まれる周波数成分を観測し、観測した周波数成分を下記 (a )の対数周波数 X上の確率密度関数として表現するものとし、

[数 91] ρ⁽^(^χ) ^(a) 前記観測した周波数成分の確率密度関数から、下記 (b)で表現される基本周波数

Fの確率密度関数を求める過程において、音モデルの多重化、音モデルのパラメ一タの推定及びモデルパラメータに関する事前分布の導入を採用し、

[数 92]

前記音モデルの多重化では、同一基本周波数に対して M種類の音モデルがあるものとして前記基本周波数カ^の m番目の音モデルの確率密度関数を p(x I F, m, μ ^ω (F, m) )と表現し、ただし、 μ ^ω (F, m)は前記 m番目の音モデルの高調波成分の大きさの比率を表すモデルパラメータであり、

前記音モデルのパラメータの推定では、前記観測した周波数成分の確率密度関数が、下記 (c)で定義した混合分布モデル p (x I Θ ^(t))から生成されたものと考え、ただし0> ^ω (F, m)は基本周波数が Fの m番目の音モデルの重みであり、

[数 93] dF

ただし上記 (c)にお、て、 θ ^ωは前記音モデルの重み ω ^(t) (F, m)と前記音モデノレの高調波成分の大きさの比率 μ ^(t) (F, m)を包含したモデルパラメータ Θ ^W = { ω ^(t) , (^t) }であり、 ω ^ω = { ω ^ω (Ρ, m) | Fl≤F≤Fh, m= l, · ··, M}であり、 _μ ^ω = { μ ⁽' ) (F, m) I Fl≤F≤Fh, m= l, · ··, M}であり、ここで Flは許容される基本周波数の下限であり、 Fhは許容される基本周波数の上限であり、

前記基本周波数 Fの確率密度関数を下記 (d)の解釈により前記重み ω (t) (F, m) から求め、

[数 94] pP⁽ ₀ (F) 二 (Fl < F < Fh) (d)

前記モデルパラメータに関する事前分布の導入では、前記モデルパラメータ Θ ^(t)に関する事前分布に基づいて前記モデルパラメータ Θ ^(t)の最大事後確率推定量を E M (Expectation-Maximization)アルゴリズムを用いて推定し、この推定から事前分布を考慮した上記 (b)の基本周波数 Fの確率密度関数と解釈できる前記重み ω ^(t) (F, m)とすべての音モデルの確率密度関数 p (X I F, m, ^ ^(t) (F, m) )の ^(t) (F, m)が表す第 h次高調波成分の大きさの c^(t) (h I F, m) (h= l, · ··, H)とを求めるために用いる下記 (e)及び (f)で表される二つのパラメータ推定値を求めるための式を定め、ただし Hは前記基本周波数の周波数成分も含めた高調波成分の数であり、 [数 95] (e)

[数 96]

(f) また上記 (e)及び (f)中にお！/、て、下記 (g)及び (h)は、下記 (i)及び (j)が 0になる無情報事前分布のときの最尤推定値であり、

[数 97]

" Fm―、 -广 "(り (x) ^)(F'^m) ?¾^^^(ί)^ 2 _dx

^ML(F，）—厶 ^φ ( ^ Σ^ ) レ) レ (り ("，レ))

—— (g)

[数 98]

w^_h(F,m) ^J-°° /FI∑di ノ ^(f)(", ），レ， ^(ί)(", )）

(h)

[数 99]

) ― —- (i) [数 100]

上記 (e)及び (f)中にお!/、て、下記 (k)は前記重み ω ^(t) (F, m)の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、下記 (1)は前記モデルパラメータ w _(F, _m)の単峰性の事前分布を考えたときに最大値を取る最も起こりやすヽパラメータであり、また上記 (i)は下記 (_k)をどれぐらヽ重視した事前分布とするかを決めるパラメータであり、上記 (j)は下記 (1)をどれぐら、重視した事前分布とするかを決めるパラメータであり、

[数 101]

[数 102]

また上記 (g)及び (h)にお、て、 ω ' ^(t) (F, m)及び μ ' ^(t) (F, m)は上記 (e)及び (f )を反復計算する際の一つ前の古いパラメータ推定値であり、 ηは基本周波数であり、 Vは何番目の音モデルであるかを示し、

上記 (e)及び (f)の二つのパラメータ推定値を求めるための式を用いた反復計算により、上記 (b)の基本周波数の確率密度関数と解釈できる前記重み co ^(t) (F, m)とすベての音モデルの確率密度関数 p (x I F, m, μ ^(t) (F, m) )の前記モデルパラメータ μ ^(t) (F, m)が表す第 h次高調波成分の大きさ c^(t) (h I F, m)とをコンピュータ内に以下の機能を実現する手段を構成して演算により求めることにより基本周波数の音高を推定する音高推定装置であって、

上記 (e)及び (f)を上記 (g)及び (h)によって前記コンピュータを利用して演算するために、上記 (g)で表現される推定値を示す計算式中の分子を下記 (m)に示す Xの関数として展開する手段を備え、但し下記 (m)に示す式にお、て ω ' ^(t) (F, m)は古い重みであり、 c' ^(t) (h I F, m)は古い第 h次高調波成分の大きさの比率であり、 Hは前記基本周波数の周波数成分も含めた高調波成分の数であり、 mは M種類の音モデルの中の何番目の音モデルかを示すものであり、 Wは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、

[数 103] (x (F + 1200 1og₂ /z,) )²

- ^ ^——

^rW²

上記（m)中の 12001og hと exp [― (x— (F + 12001og h) ) ²Z2W²]を事前に計算し

2 2

て前記コンピュータのメモリに格納する手段を備え、

上記 (e)及び (f)の二つのパラメータ推定値を求めるための式を事前に定めた回数分だけ反復計算するために、上記 (g)及び (h)の式の計算では、前記観測した周波数成分の確率密度関数の離散化後の各周波数 Xに対して以下の第 1の演算処理を Nx回実行する第 1の演算処理手段を備え、但し Nxは Xの定義域の範囲の離散化数であり、

前記第 1の演算処理手段では、 M種類の音モデルについてそれぞれ以下の第 2の演算処理を実行して、上記 (m)式の演算結果を求め、上記 (m)式の演算結果を基本周波数 Fと m番目の音モデルに関して積分して上記 (g)及び (h)の式の分母を求め前記観測した周波数成分の確率密度関数を上記 (g)及び (h)の式に代入して上記 (g)及び (h)の演算を行う第 2の演算処理手段を備え、

前記第 2の演算処理手段では、基本周波数を含む高調波成分の数 Hだけ、第 3の演算処理を実行して下記 (n)の式の演算結果を求め、下記 (n)の式の演算結果を H 個加算して上記 (m)の式の演算結果を求める第 3の演算処理手段を備え、

[数 104]

(x 一 (F + 1200 1ο¾ / ))²

前記第 3の演算処理手段では、 X - (F+ 12001og h)が 0に近い基本周波数 Fに関

2

して第 4の演算処理を Na回実行して上記 (n)の式の演算結果を求める第 4の演算処理手段を備え、但し Naは X— (F+ 12001og h)が充分 0に近い範囲の離散化後の F

2

の個数を表す小さ、正の整数であり、

前記第 4の演算処理手段では、事前に前記メモリに格納した exp [— (x- (F+ 120 Olog h))²Z²W²]を用いて下記（o)の式を求め、

2

[数 105]

）

(o) 上記（o)の式に古、重み ω ' ^(t) (F, m)をかけて上記 (n)の式の演算結果を求めることを特徴とする音高推定装置。

[7] 前記対数周波数 Xと前記基本周波数 Fの離散化幅を dとしたときに、（3WZd)より小さいもしくは近い正の整数 bを求めて、前記 Naを（2b +1)回と決定し、離散化して計算するときに X— (F+12001og h)が— b+ α, — b+l+ α, ···, 0+ α, ···, b— 1

2

したときの 0. 5以下の小数である請求項 6に記載の音高推定装置。

[8] 前記対数周波数 Xと前記基本周波数 Fの離散化幅を dとしたときに、（3WZd)より小さいもしくは近い正の整数 bを求めて、前記 Naを（2b +1)回と決定し、離散化して計算するときに X— (F+12001og h)が— b+ α, — b+l+ α, ···, 0+ α, ···, b— 1

2

したときの 0. 5以下の小数であり、

前記メモリには、 X— (F+12001og h)が— b+ a, — b + l+ a, ···, 0+ a, ···, b

2

-1+ a, b+ aの（2b + l)通りの値を取るときの exp[— (x— (F + 12001og h))

2

2W²]の値が事前に格納されて、る請求項 6に記載の音高推定装置。

[9] 前記対数周波数 Xと前記基本周波数 Fの離散化幅が 20centで前記 Wが 17cent のときに、前記 Naを 5と定め、離散化して計算するときに X— (F + 12001og h)がー 2

2

2 散化して表現したときの 0. 5以下の小数である請求項 6に記載の音高推定装置。

[10] 前記対数周波数 Xと前記基本周波数 Fの離散化幅が 20centで前記 Wが 17cent のときに、前記 Naを 5と定め、離散化して計算するときに X— (F + 12001og h)がー 2

2 + α, -1+ α, 0+ α, 1+ α, 2+ αの値を取り、ここで αは（F+ 12001og h)を離

2 散化して表現したときの 0. 5以下の小数であり、

前記メモリには、 X— (F+12001og h)が— 2+ α, — 0+ α, 1+ α, 2+ a

2

の値を取るときの exp [—（X— (F+12001og h))²Z2W²]の値が事前に格納されて

2

いる請求項 6に記載の音高推定装置。

[数 106]

ΡΦ^} (^) (^A) 前記観測した周波数成分の確率密度関数から、下記 (b)で表現される基本周波数

[数 107] ) -…… ）前記音モデルの多重化では、同一基本周波数に対して M種類の音モデルがあるものとして前記基本周波数カ^の m番目の音モデルの確率密度関数を p(x I F, m, μ ^(t) (F, m) )と表現し、ただし、 μ ^(t) (F, m)は前記 m番目の音モデルの高調波成分の大きさの比率を表すモデルパラメータであり、

[数 108]

(c) ただし上記 (c)にお、て、 θ ^(t)は前記音モデルの重み ω ^(t) (F, m)と前記音モデノレの高調波成分の大きさの比率 μ ^(t) (F, m)を包含したモデルパラメータ Θ ^W = { ω ^(t) ^(t)}であり、 ω^ω = {ω^ω(Ρ, m) | Fl≤F≤Fh, m=l, ···, M}であり、 μ^ω = {μ⁽' )(F, m) I Fl≤F≤Fh, m=l, ···, M}であり、ここで Flは許容される基本周波数の下限であり、 Fhは許容される基本周波数の上限であり、

[数 109]

M

p¾( ) 二 ' ( (F,m) (Fl < < Fh) (d)

m=l

前記モデルパラメータに関する事前分布の導入では、前記モデルパラメータ Θ ^(t)に関する事前分布に基づいて前記モデルパラメータ Θ ^(t)の最大事後確率推定量を E M (Expectation-Maximization)アルゴリズムを用いて推定し、この推定から事前分布を考慮した上記 (b)の基本周波数 Fの確率密度関数と解釈できる前記重み ω ^(t) (F, m)とすべての音モデルの確率密度関数 p (X I F, m, ^ ^(t) (F, m) )の ^(t) (F, m)が表す第 h次高調波成分の大きさの c^(t)(h I F, m) (h=l, ···, H)とを求めるために用いる下記 (e)及び (f)で表される二つのパラメータ推定値を求めるための式を定め、ただし Hは前記基本周波数の周波数成分も含めた高調波成分の数であり、

[数 110] w^Hr.m) =————————^^————

1十^ )

[数 111]

(f) また上記 (e)及び (f)中にお、て、下記 (g)及び (h)は、下記 (i)及び (j)が 0になる無情報事前分布のときの最尤推定値であり、 [数 112]

(g)

[数 113]

(h)

[数 114]

[数 115]

上記 (e)及び (f)中にお、て、下記 (k)は前記重み ω ^ω (F, m)の単峰性の事前分布を考えたときに最大値を取る最も起こりやすいパラメータであり、下記 (1)は前記モデルパラメータ ^(t) (F, m)の単峰性の事前分布を考えたときに最大値を取る最も起こりやす、パラメータであり、また上記 (i)は下記 (k)をどれぐら!/、重視した事前分布とするかを決めるパラメータであり、上記 (j)は下記 (1)をどれぐら、重視した事前分布とするかを決めるパラメータであり、

[数 116]

[数 117]

また上記 (g)及び (h)にお、て、 ω ' ^(t) (F, m)及び μ " ^(t) (F, m)は上記 (e)及び (f )を反復計算する際の一つ前の古いパラメータ推定値であり、 ηは基本周波数であり、 Vは何番目の音モデルであるかを示し、

上記 (e)及び (f)の二つのパラメータ推定値を求めるための式を用いた反復計算により、上記 (b)の基本周波数の確率密度関数と解釈できる前記重み co ^(t) (F, m)とすベての音モデルの確率密度関数 p (x I F, m, μ ^(t) (F, m) )の前記モデルパラメータ μ ^(t) (F, m)が表す第 h次高調波成分の大きさ c^(t) (h I F, m)とをコンピュータを利用した演算により求めるため前記コンピュータにインスト一ノレされて前記コンピュータ内に必要な機能を実現するための音高推定用プログラムであって、

上記 (e)及び (f)を上記 (g)及び (h)によって前記コンピュータを利用して演算するために、上記 (g)で表現される推定値を示す計算式中の分子を下記 (m)に示す Xの関数として展開する機能、但し下記 (m)に示す式にお、て ω ' ^(t) (F, m)は古、重みであり、 c' ^(t) (h I F, m)は古い第 h次高調波成分の大きさの比率であり、 Hは前記基本周波数の周波数成分も含めた高調波成分の数であり、 mは M種類の音モデルの中の何番目の音モデルかを示すものであり、 Wは各高調波成分をガウス分布で表現するときのガウス分布の標準偏差であり、

[数 118]

(a: - ( + 12001og₂ ?,))²

(り (F, m) 5 ) (/i|F，m) exp( m 上記（m)中の 12001og hと exp [― (x— (F + 12001og h) ) ²Z2W²]を事前に計算し

2 2

て前記コンピュータのメモリに格納する機能、

上記 (e)及び (f)の二つのパラメータ推定値を求めるための式を事前に定めた回数分だけ反復計算するために、上記 (g)及び (h)の式の計算では、前記観測した周波数成分の確率密度関数の離散化後の各周波数 Xに対して以下の第 1の演算処理を Nx回実行する機能、但し Nxは Xの定義域の範囲の離散化数であり、

前記第 1の演算処理では、 M種類の音モデルについてそれぞれ以下の第 2の演算処理を実行して、上記 (m)式の演算結果を求め、上記 (m)式の演算結果を基本周波数 Fと m番目の音モデルに関して積分して上記 (g)及び (h)の式の分母を求め前記観測した周波数成分の確率密度関数を上記 (g)及び (h)の式に代入して上記 (g) 及び (h)の演算を行う機能、

前記第 2の演算処理では、基本周波数を含む高調波成分の数 Hだけ、第 3の演算処理を実行して下記 (n)の式の演算結果を求め、下記 (n)の式の演算結果を H個加算して上記 (m)の式の演算結果を求める機能、

[数 119]

' 、、

„ )ご( ( |

2

第 4の演算処理を Na回実行して上記 (n)の式の演算結果を求める機能、但し Naは x (F + 12001og h)が充分 0に近い範囲の離散化後の Fの個数を表す小さい正の整

2

数であり、

前記第 4の演算処理では、事前に前記メモリに格納した exp [— (x- (F+12001og

2 h) ) ²/2W²]を用いて下記（o)の式を求め、

[数 120] 、 1 ₍ ( + 12001og₂ /i))²

， ^m) ^exp(一^—— —)

上記（o)の式に古、重み ω ' ^(t) (F, m)をかけて上記 (n)の式の演算結果を求める機能を前記コンピュータ内に実現するための音高推定用プログラム。

[12] 前記対数周波数 Xと前記基本周波数 Fの離散化幅を dとしたときに、 (3W/d)より小さいもしくは近い正の整数 bを求めて、前記 Naを（2b +1)回と決定し、離散化して計算するときに X— (F+12001og h)が— b+ α, b+l+ α, ···, 0+ α, ···, b— 1

2

したときの 0.5以下の小数である請求項 11に記載の音高推定用プログラム。

[13] 前記対数周波数 Xと前記基本周波数 Fの離散化幅を dとしたときに、 (3W/d)より小さいもしくは近い正の整数 bを求めて、前記 Naを（2b +1)回と決定し、離散化して計算するときに X— (F+12001og h)が— b+ α, — b+l+ α, ···, 0+ α, ···, b— 1

2

したときの 0. 5以下の小数であり、

前記メモリには、 X— (F+12001og h)が一 b+ひ， b + l+ひ， …， 0+ひ， …， b

2

2W²]の値が事前に格納されて、る請求項 11に記載の音高推定用プログラム。

[14] 前記対数周波数 Xと前記基本周波数 Fの離散化幅が 20centで前記 Wが 17cent のときに、前記 Naを 5と定め、離散化して計算するときに X (F + 12001og h)がー 2

2

2 散化して表現したときの 0. 5以下の小数である請求項 11に記載の音高推定用プログラム。

[15] 前記対数周波数 Xと前記基本周波数 Fの離散化幅が 20centで前記 Wが 17cent のときに、前記 Naを 5と定め、離散化して計算するときに X (F + 12001og h)がー 2

2

2 散化して表現したときの 0. 5以下の小数であり、

前記メモリに ίま、 X— (F+12001og h)力ー 2+α, — 0+ α, 1+ α, 2+ a

2

の値を取るときの exp [ （X— (F+12001og h))²Z2W²]の値が事前に格納されて

2

V、る請求項 11に記載の音高推定用プログラム。