情報処理学会研究報告

Vol.2013-MUS-99 No.8
2013/5/11

IPSJ SIG Technical Report

スペクトル包絡と基本周波数の同時推定のための無限カーネル線形予測分析法
吉井 和佳1,a)

後藤 真孝1,b)

概要：本稿では，音声信号のスペクトル包絡と基本周波数とを同時に推定するための新しい線形予測分析

法について述べる．従来，ソース・フィルタ理論に基づく線形予測分析法では，所与の観測信号はガウス
性白色雑音を入力信号とする自己回帰系からの出力信号であると仮定して，全極型フィルタの係数を推定することが行われていた．しかし，声帯振動に起因する周期パルス（周波数領域では調波構造）が入力信号である場合には，推定されたスペクトル包絡（全極型フィルタの周波数応答）は，調波構造の倍音周波数で不必要に大きなピークをもつ問題があった．この問題を解決するため，入力信号と出力信号との関係をノンパラメトリックベイズガウス過程回帰モデルで表現する無限カーネル線形予測分析法を提案する．
本手法では，異なる基本周波数に対応する可算無限個のカーネルを考え，それらの凸結合で入力信号の周期性を表現する．ここで，無限個の非負の重みに対してガンマ過程事前分布を導入すると，スパースなマ
ルチカーネル学習を行うことができる．すなわち，全極型フィルタの係数を推定すると同時に，基本周波数に対応する優勢なカーネルを同定できる．本手法を話声・歌声信号に対して適用し，基本周波数をもつ
有声区間を同定しながら，基本周波数の影響を考慮したスペクトル包絡を推定できることを確かめた．

1. はじめには，音源信号がガウス性白色雑音であれば，全極型フィルタの係数を最尤推定の枠組みで精度良く求めることができ音響信号のスペクトル包絡推定は，話声・歌声信号分析る．しかし，観測信号が明確な音高をもっている，すなわの基礎をなす重要な技術である．音声分野におけるこれまち，音源信号が周期的であると，推定されたフィルタの周での研究により，人間の発声機構はソース・フィルタ理論

波数応答（スペクトル包絡）は調波構造のピークの位置ででよく説明できることが知られている（図 1）．具体的に不必要に鋭いピークをもつようなバイアスがかかる．

は，声帯振動に起因する音源信号が声道を通ることでそのこの問題を解決するため，これまで多くの研究が行われ音響特性が変化する過程を考える．時間領域では，音源信てきた．El-Jaroudi ら [2] は，全極型フィルタの周波数応答

号に声道フィルタのインパルス応答が畳みこまれた出力信を調波構造の離散的なピークに対してフィットさせる手法

号が観測信号であると解釈できる．一方，周波数領域では，(Discrete All-pole Filtering: DAP) を提案している．この観測スペクトルの微細構造と包絡構造がそれぞれ，音源信

手法は，Badeau ら [3] によって自己回帰・移動平均モデル

号とフィルタの周波数特性に対応していると仮定することにも適用可能なように拡張されている．一方，Oudot ら [4]が一般的である．本研究ではこのような仮定のもと，観測は，全極型フィルタがなめらかな周波数応答をもつような信号が与えられたときに，フィルタの周波数応答，すなわ

制約をかける手法を提案している．Villavicencio ら [5] は，

ちスペクトル包絡を推定する問題に取り組む．

ケプストラム平滑化を反復的に適用する手法を提案してい線形予測分析 (Linear Prediction: LP) [1] は，スペクトる．線形予測分析以外のスペクトル包絡推定法として，河ル包絡推定のための数学的に確立された方法のひとつであ原ら [6] は，音声スペクトルを周期・非周期成分とスペクる．線形予測分析では，観測信号は自己回帰過程に従う，

トル包絡に精度よく分離できる分析合成系 STRAIGHT を

すなわち，フィルタは全極型伝達関数で記述できると仮定

開発している．中野ら [7] は基本周波数の影響を回避する

する．声道は単純な音響管の連結であるとみなすと鼻子音

ため，隣接するフレームにわたってスペクトル包絡を平均以外には反共振は存在せず，人間の聴覚はスペクトルの化する手法を提案している．これらの手法はスペクトル包ピーク（フォルマント）に敏感であることから，この仮定

絡を精度良く求めることができるが，基本周波数の値が既は妥当であると考えられている．古典的な線形予測分析で知であることが前提であった．

1
a)
b)

産業技術総合研究所
Umezono 1-1-1, Tsukuba, Ibaraki 305–8568, Japank.yoshii(at)aist.go.jp
m.goto(at)aist.go.jpⓒ 2013 Information Processing Society of Japan

近年，基本周波数とスペクトル包絡をモデル化するうえで，確率モデルに基づくアプローチが有望視されている．
佐宗ら [8] は自己回帰隠れマルコフモデル (AR-HMM) と1

情報処理学会研究報告

Vol.2013-MUS-99 No.8
2013/5/11

IPSJ SIG Technical Report
周期信号

時間

周波数

ソース

白色雑音

2.1 確率モデルの定式化

自己回帰系
インパルス応答の畳みこみ

まず，観測信号が自己回帰過程に従うと仮定する基本的
観測信号

時間

フィルタ

時間

周波数応答の乗算
時間

周波数

図 1

周波数
周波数

発声機構に対するソース・フィルタモデル呼ばれる，状態遷移が循環するように拘束をかけた HMMを用いて音源信号の周期性を表現する手法を提案している．
戸田ら [9] は，あらかじめ基本周波数の値を与える必要があるものの，音源信号に起因する調波構造の時間的なダイナミクスをトラジェクトリ HMM を用いて表現する手法を提案している．亀岡ら [10, 11] は，基本周波数とスペク
トル包絡を同時に推定するための先駆的な研究を行っている．例えば，異なる基本周波数に対応する多数の周期カーネルを内包するガウス過程に基づくマルチカーネル線形回
帰分析法 (Multiple Kernel Linear Prediction: MKLP) を提案している [10]．優勢なカーネルを決定（基本周波数を推定）するには，スペクトル包絡を推定すると同時にカー
ネルの重みに対する事後分布最大化推定が行われる．我々は，この研究をさらに発展させることを試みる．
本稿では，近年着目されているノンパラメトリックベイズ
理論とカーネル法の強みを同時に取り入れた無限カーネル
線形予測分析法 (Inﬁnite Kernel Linear Prediction: IKLP)を提案する．我々は，MKLP におけるカーネルの個数が無限に発散した極限を考え，それら無限個のカーネルの重みに対してガンマ過程事前分布を仮定する．変分ベイズ法を用いて事後分布推定を行うことで，ほとんど全てのカーネルの重みがゼロにほぼ等しくなり，基本周波数に対応するカーネルの重みのみが有意な値をもつようなスパースな学
習をスペクトル包絡推定と同時に行うことができる．本研
究は，線形予測分析法の背後にある観測信号の生成過程に対する仮定を明らかにすることで信号処理分野に貢献するだけではなく，マルチカーネル学習 [12] に対する効率的な収束保証付きの最適化アルゴリズムを提案することで機械学習分野にも貢献することができると考える．

2. 線形予測分析
本章では，音声・歌声などの観測信号からスペクトル包絡（全極型フィルタの係数）を推定するための確率モデル
について概観し，その最新の手法として亀岡らの手法 [10]な定式化について説明する．あるフレームに含まれる M
個の連続するサンプルを x = (x1 , x2 , · · · , xM )T とする．
いま，x が P 次の自己回帰過程に従うと仮定すると，

xm =

P

p=1 ap xm−p + m(1)と書くことができる．ここで，a = (a1 , · · · , aP )T は全極型フィルタの P 個の係数であり，線形予測係数と呼ばれる．

 = (1 , · · · , M )T は誤差項である．ソース・フィルタ理論の観点からは， が声帯振動に起因する音源信号（ソース）に対応し，a が声道（フィルタ）の反響特性を表していると解釈できる．この自己回帰モデルは，m を入力として xmを出力する線形系とみなすことができ，その全極型伝達関
数は A(z) = 1/(1 − a1 z −1 · · · − aP z −P ) で与えられる．すなわち，x および  の Z 変換を X(z) および E(z) とすると，X(z) = E(z)A(z) が成り立つ．全極型フィルタの周波数応答（スペクトル包絡）は，m を周波数ビンのインデッ
クすると考えると |A(e2πim/M )|2 で与えられる（図 1）．
我々の目的は，観測信号 x が与えられたときに，確率的な枠組みのもとでフィルタ係数 a を求めることである．この問題はフィルタ係数 a および音源信号  が両方とも未知のもとでは不良設定問題であるため，音源信号  の性質に関するなんらかの仮定が必要になる．一般的には， はガ
ウス性白色雑音であると仮定する．

 ∼ N (0, νI)(2)

ここで，ν はノイズの分散であり，I は単位行列である．すなわち， の M 個の成分は N (0, ν) に従って独立同分布することを示している．ここで，行列 Ψ ∈ RM ×M および行
列 X ∈ RM ×P を
⎡ 1

⎢ −a . . .
⎢ 1⎢ .
⎢ . ..
.
.
Ψ =⎢
⎢ −aP
⎢
⎢
..
⎣
.

⎤

..

0
.
..

0 −a

P

.

..

.

⎡

⎥
⎢
⎥
⎢
⎥
⎢
⎥
⎢
⎥, X = ⎢
⎥
⎢
⎥
⎢
⎥
⎢
⎦
⎣

· · · −a1 1

0
x1
.
.
.
.
.
.
.
.
.

···
..
.
..

.

0
.
.
.
0
x1
.
.
.

⎤
⎥
⎥
⎥
⎥
⎥(3)
⎥
⎥
⎥
⎦

xM −1 · · · xM −Pとすると，式 (1) は  = Ψx と簡潔に書けて，

x = Ψ−1 (4)を得る．式 (2) および (4) を用いると，観測信号 x の尤度はx ∼ N (0, νΨ−1 Ψ−T )(5)を紹介する．まず，音源信号がガウス性白色雑音に従うとで与えられる．式 (5) が古典的な線形予測分析の確率モデ仮定する古典的な自己回帰モデルについて説明する．次ルである．この確率モデルに対して最尤推定を行う場合にに，音源信号が周期性をもつ場合でも精度のよい推定が可は，最適なフィルタ係数 a は正規方程式 X T Xa = X T x

能なカーネル化されたモデルについて紹介する．

の解として求めることができる．

ⓒ 2013 Information Processing Society of Japan2

情報処理学会研究報告

Vol.2013-MUS-99 No.8
2013/5/11

IPSJ SIG Technical Report

音源信号  が式 (2) で与えられる等方的なガウス分布に2.3 マルチカーネル学習

従う場合には，このスペクトル包絡推定法は有効に機能す

次に，音源信号  の性質を反映するようなグラム行列 Kる．しかし，観測信号が音声信号や歌声信号である場合にの設計法について述べる．本稿では，観測信号 x（あるいは，声帯の周期振動に起因する明確な周期性をもっているは音源信号 ）が基本周波数をもつ周期信号である場合をため，なんらかの対策が必要になる．

考えているので，K として周期カーネルを利用するのが自2然である．例えば，k(t, t ) = exp(−2 sin2 (π t−t
T )/l ) は良


2.2 ガウス過程に基づくカーネル化
亀岡ら [10] は，式 (5) で与えられる確率モデルをガウス
過程回帰の観点でカーネル拡張する方法を提案している．
これまで，音源信号  はガウス性白色雑音であると仮定してきた．一方，ここでは，音源信号 (t) は時間 t 上に定義された連続関数であると考え，時刻 t から関数 (t) を線形回帰する問題について考える．いま，{m }M
m=1 は，時刻{tm }M
m=1 における関数 (t) の出力値であるとみなす．我々の目的は，連続関数 (t) を J 個の基底関数 {φj (t)}Jj=1 の和で表現することである．
J
(t) = j=1 wj φj (t) + η(t) = φ(t)T w + η(t)(6)

φ(t) = (φ1 (t), · · · , φJ (t))T とした．時刻 {tm }M
m=1 における関数 η(t) の出力値を並べたものを η = (η1 , · · · , ηM )T とし，計画行列を Φ = (φ(t1 ), · · · , φ(tM ))T ∈ RM ×J とすると，式 (6) の時刻 {tm }M
m=1 における回帰モデルは(7)で与えられる．いま，重み ω および誤差信号 η が等方的なガウス分布に従うと仮定すると

w ∼ N (0, νw I),η ∼ N (0, νe I)(8)と書ける．ここで，νw および νe はガウス分布の分散を表がガウス分布であることから，関数 (t) はガウス過程 [13]に従うことが分かる．その振る舞いはグラム行列（カーネル）K = ΦΦT によって規定される．カーネル K の各要
素は基底関数の内積として定義されている．

(10)

一方，任意の正定値行列はグラム行列として有効であることが知られており，基底関数を明示せずに Km,m = k(tm , tm )として直接 K を計算することもできる（カーネルトリック）．ここで，k(t, t ) はカーネル関数と呼ばれる．このとき，式 (4) および (9) を用いると，観測信号 x の尤度はx ∼ N (0, Ψ−1 (νw K + νe I)Ψ−T )つ正弦波の足し合わせで表現する方法を提案している．

φj (t) =H

h=1 sint−c

2πh T j(cj は位相)(12)

このとき，グラム行列 K は式 (10) から求められる．しかし，実際には音源信号  の周期 T は未知であり，観測信号
この問題に対する強力な解法のひとつに，マルチカーネ
ル学習 [12] が知られている．具体的には，グラム行列 K
を I 個の異なるグラム行列の重みつけ和として表現する．

K=

I

i=1 θi K i(13)

ここで，θ = {θi }Ii=1 はカーネルの重みであり，K i は周期

Ti をもつ周期カーネルである．カーネルの重みは，観測信号 x において各周期カーネルがどの程度優勢であるか，すなわち，基本周波数がどのあたりに存在するかを示していx ∼ N 0, Ψ−1 νw

I

i=1 θi K i + νe I

Ψ−T

(14)

フィルタ係数 a およびカーネルの重み θ は，EM アルゴ
リズムを用いて推定することができる．文献 [10] において(9)を得る．任意の時刻 {tm }M
m=1 における関数 (t) の周辺分布

Km,m = φ(tm )T φ(tm )亀岡ら [10] は，各基底関数を H 個の等しいパワーを持る．最終的に，観測信号 x の尤度は以下で与えられる．

す．式 (7) および (8) を用いると，最終的に

 ∼ N (0, νw ΦΦT + νe I)

き，全ての基底関数が周期 T （基本周波数は 1/T ）の周期関数であることが暗黙的に仮定されている．

x が与えられたときに，周期 T を推定する必要がある．

ここで，η(t) は誤差関数，ω ∈ RJ は基底関数の重みであり，

 = Φw + η

く知られた周期カーネルで，その周期は T である．このと(11)は，異なる周期 {Ti }Ii=1 をもつ数百個のカーネルを考え，カーネルの重み θ に対する事前分布として一般化ガウス分
布をおくことで，θ を事後確率最大化 (MAP) 推定により求めている．しかし，MAP 推定の枠組みでは，θ を完全にスパースに導くことは原理的にできない．

3. 無限カーネル線形予測分析
本章では，音声・歌声信号のスペクトル包絡と基本周波数を確率的な枠組みで同時に推定するため，ノンパラメト
リックベイズモデルに基づく無限カーネル線形予測分析
法 (IKLP) を提案する．まず，式 (14) で与えられるマルチカーネル線形予測分析法 (MKLP) の確率モデルにおいて，カーネルの個数を無限に発散させたときの極限（I → ∞）について考える．次に，確率モデルのベイズ的な取り扱いで与えられる．これはガウス過程回帰モデル [10] であり，を可能にするため，未知のパラメータに対して適切な事前式 (5) をその特別な場合として含んでいる．実際，ΦΦ = I

分布を設計する．最後に，未知のパラメータの事後分布をであれば，すなわち，J 個の基底関数が互いに独立であれ計算するため，変分ベイズ法に基づく効率的で収束の保証

ば，ν = νw + νe とすると式 (5) が得られる．

された最適化アルゴリズムの導出を行う．

Tⓒ 2013 Information Processing Society of Japan3

情報処理学会研究報告

Vol.2013-MUS-99 No.8
2013/5/11

IPSJ SIG Technical Report

3.1 ノンパラメトリックベイズモデル

変分ベイズ法では，変分事後分布の真の事後分布に対す式 (14) において I → ∞ とすると，観測信号 x の尤度はx ∼ N 0, Ψ−1 νw

I→∞
i=1

θi K i + νe I Ψ−T(15)で与えられる．まず，無限次元の非負ベクトル θ に対する事前分布として，ガンマ過程を利用する．具体的には，

θi ∼ Gamma(α/I, α)(16)

として，打ち切りレベル I を無限大に発散させたときに，

θ は，集中度 α をもつガンマ過程から得られる無限次元の非負系列となる．このとき，任意の正整数  > 0 に対して，

θi >  となる要素の個数 I + はほとんど確実に有限である

ことが証明されており，無限次元の空間においてスパースな学習が可能である根拠となっている．現実的には，打ちるカルバック・ライブラー (KL) ダイバージェンスを単調減少させるように，各因子を反復的に最適化する．これは，対数周辺尤度 log p(x) の変分下限 L を単調増加させることと等価である．ここで，L は以下の通り与えられる．

log p(x) ≥ E[log p(x|θ, a, νw , νe )](20)+ E[log p(θ)] + E[log p(a)] + E[log p(νw )] + E[log p(νe )]− E[log q(θ)] − E[log q(a)] − E[log q(νw )] − E[log q(νe )] ≡ Lしかし，右辺の第一項（対数尤度関数の期待値）は依然として解析的に計算ができないため，L ≥ L となるようなさらなる変分下限 L を構成し，L を逐次最大化することを考える．このとき，各因子の更新則は切りレベル I を集中度 α より十分大きくしておけば，θ 中

q(θ) ∝ p(θ) exp(Eq(a,νw ,νe ) [log q(x|θ, a, νw , νe )])のいくつかの要素のみがゼロより大きな有意な値を持つこ

q(νw ) ∝ p(νw ) exp(Eq(θ ,a,νe ) [log q(x|θ, a, νw , νe )]) (21)とが期待できる．
完全なベイズ的な取り扱いのため，ガウス分布の分散 νw

および νe に対する事前分布に，ガンマ分布を利用する．

νw ∼ Gamma(aw , bw ),νe ∼ Gamma(ae , be )(17)

ここで，a∗ および b∗ はガンマ分布の形状パラメータおよびレートパラメータである．さらに，フィルタ係数 a に対する事前分布として，ガウス分布を利用する．

a ∼ N (0, λI)q(νe ) ∝ p(νe ) exp(Eq(θ ,a,νw ) [log q(x|θ, a, νw , νe )])で与えられる．ここで，q(x|θ, a, νw , νe ) は p(x|θ, a, νw , νe )の変分下限であり，式 (24) で与えられる．ただし，共役性の問題から a の完全なベイズ的な取り扱いは困難であるため，q(a) = δa∗ (a) であると仮定する．ここで，δa∗ は，ある a∗ において関数値は無限大となり，それ以外はゼロとなるようなディラックのデルタ関数である．

3.2.1 行列に関する不等式(18)

解析的に計算可能な変分下限 L を導出するためには，行列に関する 2 つの不等式を用いる必要がある．まず，行列ここで，λ は超パラメータである．
基本周波数を推定するには，重みの期待値 E[θi ] が最大となるカーネル K i を同定すればよい．このとき，基本周や行列変数関数に関する重要な概念を整理しておく．
定義 1 (行列の半正定値性) ある実対称行列 A が半正定

波数は 1/Ti となる．スペクトル包絡は，フィルタ係数の期値性を満たすとは，任意の実ベクトル z に対して z T Az ≥ 0

して，MKLP [10] と同様に，分散の比率 E[νw ]/E[νw + νe ]る，あるいは A = Z T Z となるような実行列 Z が存在す

待値 E[a] を用いて計算できる．このモデルの副次的効果とが成立する，あるいは A の全ての固有値がゼロ以上であから有声・無声の判定が可能である．

ることを言う．これらの条件はすべて等価である．

3.2 変分ベイズ法

f (·) が凸であるとは，任意の実数 0 ≤ λ ≤ 1 に対して

定義 2 (関数の凸性・凹性) 行列変数スカラー値関数

我々の目的は，観測信号 x が与えられたときに，未知パ
ラメータの同時的な事後分布 p(θ, a, νw , νe |x) をベイズのλf (A) + (1 − λ)f (B) ≥ f (λA + (1 − λ)B) が成り立つことを言う．関数 f (·) が凹であるとは，λf (A) + (1 − λ)f (B) ≤定理 p(θ, a, νw , νe |x) = p(x, θ, a, νw , νe )/p(x) に従って計f (λA + (1 − λ)B) が成り立つことを言う．

的に計算することは困難であるが，変分ベイズ法を用いれであるとすると，以下の 2 つの補題が得られる．

算することである．正規化項である周辺尤度 p(x) を解析

ば真の事後分布の近似を効率よく求めることができる．具
体的には，真の事後分布 p(θ, a, νw , νe |x) を，以下のように因子分解可能な変分事後分布q(θ, a, νw , νe ) = q(a)q(νw )q(νe )補題 1 関数 f (V ) = log |V | は凹関数である.

補題 2 関数 g(V ) = z T V −1 z は凸関数である.
紙面の都合上これらの証明は省略するが，定義 1 および定義 2 に従えば，簡単に確認することができる．

iq(θi )(19)で近似することを考える．これは，事後分布において各パ
ラメータの独立性を仮定していることを意味しており，真の事後分布との間にはいくらかの乖離が存在する．
ⓒ 2013 Information Processing Society of Japan

次に，V を任意の半正定値行列，z を任意の実ベクトル

いま，各補題に関して不等式を導くことを考える．まず，
凹関数 f (V ) に関して，任意の半正定値行列 Ω を展開点とした 1 次のテイラー展開を考えると以下を得る．

log |V | ≤ log |Ω| + tr(Ω−1 V ) − M

(22)

4

情報処理学会研究報告

Vol.2013-MUS-99 No.8
2013/5/11

IPSJ SIG Technical Report

ここで，M は行列 V のサイズであり，等号成立条件は

Ω = V である．次に，凸関数 g(V ) に関して，澤田ら [14]4. 実験
本章では，無限カーネル線形予測分析法 (IKLP) の基本

によって提案された行列不等式を適用すると以下を得る．
−1
I
I
zT
z ≤ i=1 z T ΥTi V −1
(23)i Υi zi=1 V i

的な振る舞いを確認するために行った実験について述べる．

ここで，{V i }Ii=1 は任意の半正定値行列の集合であり，

4.1 実験条件{Υi }Ii=1 は 足 し 合 わ せ る と 単 位 行 列 に な る よ う な 任 意の 行 列 の 集 合（ 補 助 変 数 ）で あ る ．等 号 成 立 条 件 は ，I
Υi = V i ( i =1 V i )−1 である．この不等式を証明するには，ラグランジュの未定乗数法を用いて，右辺の最小値が左辺に等しくなることを確かめればよい．

実験には，16kHz でサンプリングされた 2 種類の音響信号を用いた．ひとつは，RWC 研究用音楽データベース：
音楽ジャンル RWC-MDB-G-2001 [16] の楽曲 No.91 の男性の無伴奏歌唱のうち冒頭部分 6.62s である．基本周波数
軌跡の正確なアノテーションを行い，STRAIGHT [6] を用

3.2.2 変分下限と反復最適化

いて分析・再合成を行った信号を観測信号として用いた．

まず，解析的に計算可能な変分下限 L の導出を行う．半正定値行列 K を K = νw i θi K i + νe I とし，式 (22) おいる女性話者の 4s の音声信号 FSUSA101 である．これらえられる L の第一項) の変分下限は以下で与えられる．

線形予測分析法への悪影響が大きいと考えられる．

よび (23) を用いると，E[log p(x|θ, a, νw , νe )] (式 (20) で与M
1
1
log(2π) − E[log |K|] − E[xT ΨT K −1 Ψx]
2
2
2
1
1
1
−1
−1
E[ν
θ
]
tr(Ω
K
)
−
) + const.
E[ν≥ − log |Ω| −w iie ] tr(Ω
2
2 i
2




1
1
−1
−
E 1 xT ΨT ΥT
E 1 xT ΨT ΥTi K i Υi Ψx −0 Υ0 Ψx2 i νw θi
2 νe
(24)

E[log p(x|·)] = −ここで，Ω は任意の半正定値行列であり，Υ = {Υi }I→∞i=0

は足すと単位行列になるような補助変数である．右辺を最大化する，すなわち，等号が成立するときの条件は
Ω = E[νw ] i E[θi ]K i + E[νe ]I
(25)
−1

−1

Υ0 = E ν1e
S −1
(26)

−1
で与えられる．ここで，S = i E νw1θi
K i + E[ ν1e ]−1 I
Υi = E νw1θi

K i S −1 ,とした．いま，式 (24) はパラメータ自身とその逆数に関する期待値計算を含んでいることに注意する．すなわち，十分統計量は x および 1/x である．一方，ガンマ分布の十分
統計量は log(x) および x であるので，変分事後分布は一般
化逆正規分布 (GIG) 分布で与えられる [15]．

もうひとつは，ATR 音声データベース [17] に収録されての信号は基本周波数が 300Hz 以上の帯域に存在し，従来のIKLP は各短時間フレームごとに独立して行う．フレー
ム長は 2048 点 (M = 2048)・シフト長は 160 点とし，窓関数にはガウス窓を用いた．超パラメータは α = 1.0,

aw = bw = ae = be = 1.0, P = 30, λ = 0.1 とした．また，

100Hz から 400Hz まで 6 セント間隔で 400 個の基本周波数に対応するカーネル {K i }I=400を準備した．
i=1

4.2 実験結果

図 2 および図 3 に示した実験結果から，IKLP はスペク

トル包絡と基本周波数を同時に推定することができるだけではなく，MKLP と同様に有声・無声区間の検出を行える
ことが確認できた．図 2 において，4.2s から 4.6s 付近に存在する有声区間の検出が不安定になっている原因は，調波構造の高次倍音成分が比較的弱く，音源信号の性質を周期
カーネルよりも単位行列カーネル（白色雑音）を用いて表現する方が適切であると判断されたからである．6.0s 付近に存在するビブラートにおいても，基本周波数を精度よく(27)推定することができた．図 2 および図 3 の 2 列目および 3q(νw ) = GIG(νw |γw , ρw , τw ), q(νe ) = GIG(νe |γe , ρe , τe )列目を比較すると分かる通り，推定されたスペクトル包絡q(θi ) = GIG(θi |γi , ρi , τi )(ρ/τ )γ/2

ただし，GIG(x|γ, ρ, τ ) = 2Kγ (√ρτ ) xγ−1 e−(ρx+τ /x)/2 である．このとき，変分パラメータの更新則は以下となる．

γi =α/I, ρi = 2α + E[νw ] tr(Ω−1 K i ) 
τi =E ν1w xT ΨT ΥTi K −1i Υi Ψxγw =aw , ρw = 2bw + i E[θi ] tr(Ω−1 K i )  τw = i E θ1i xT ΨT ΥTi K −1i Υi Ψxは，調波構造の倍音ピークの影響をほとんど受けていない．

IKLP は理論的に妥当ではあるものの，基本周波数推定においてしばしば半ピッチ誤りを起こすことがあった．この理由は，式 (15) で与えられる尤度関数は，モデルの共分
散 Ψ−1 KΨ−T が観測した共分散 xxT を過小評価する場合(28)

γe =ae , ρe = 2be + tr(Ω−1 ), τe = xT ΨT ΥT0 Υ0 Ψxには大きなペナルティをかけるが，過大評価しても小さなぺナルティしか与えないからである．このことは，式 (15)が板倉・齋藤 (IS) ダイバージェンス [1] と密接に関係していることを示唆する．実際，K が単位行列であれば，IKLP

フィルタ係数 a の MAP 推定値は，L の偏微分をゼロとは IS ダイバージェンス基準の自己回帰モデルの最適化問おくことで求まる．具体的には，正則化付きの正規方程式

題に帰着する．IS ダイバージェンスは KL ダイバージェン

X + λI)a = X Σ x の解で与えられる．ただ

T
−1
T
1
1
i E νw θi Υi K i Υi + E νe Υ0 Υ0 とした．

スのような凸性を持たないため，初期値依存性が高く，局

T

(X Σ

−1

T

−1

し，Σ−1 =

ⓒ 2013 Information Processing Society of Japan

所解に陥りやすいと考えられている [18]．この問題を解決

5

情報処理学会研究報告

Vol.2013-MUS-99 No.8
2013/5/11

IPSJ SIG Technical Report

8kHz

8kHz

8kHz

400Hz

1.0

観測スペクトログラム

8kHz

再構成スペクトログラム

8kHz

推定されたスペクトル包絡8kHz推定された基本周波数 (最も優勢なカーネル)

正解のF0軌跡400Hz

有声・無声の度合い
1.0
1s

2s

3s

4s

5s

6s

[5]

[6]

本稿では，マルチカーネル学習の枠組みで，スペクトル
包絡と基本周波数を同時に推定できるノンパラメトリック
ベイズモデルについて述べた．実験の結果，調波構造を考[7]慮したなめらかなスペクトル包絡が推定できることが確認[8]できた．従来法との詳細な比較実験は課題である．
今後の展開には，いくつかの興味深い方向性が考えられ[9]る．まず，推定精度を向上させるためには，音源信号  に対する精緻なモデル [19, 20] に基づいてカーネル関数 Φ を設計し，カーネルの周期パラメータ T 自体を経験ベイズ法の枠組みで最尤推定する方向性が考えられる．また，提案

[10]
[11]

手法は IS ダイバージェンス基準の線形予測分析 [1] の自然な拡張になっていることから，同様のマルチカーネル学[12]

習の枠組みで，IS ダイバージェンス基準の非負値行列因子[13]

分解 (NMF) [15, 21] の本質的な拡張が可能になるはずである．我々はすでに，従来の NMF における基底ベクトルの[14]

要素間の相関構造を考慮した新しい音響信号分解法について研究を進めており，優れた音源分離結果を得ている．詳[15]

細については稿を改めて報告したい．

[16]

謝辞: 本研究の一部は，JSPS 科研費 23700184 および

JST OngaCREST プロジェクトの支援を受けた．
参考文献
[1]
[2]
[3]

[4]

F. Itakura and S. Saito. Analysis synthesis telephonybased on the maximum likelihood method. ICA, 1968.
A. El-Jaroudi and J. Makhoul. Discrete all-pole modeling. IEEE Trans. on SP, 39(2):411–423, 1991.
R. Badeau and B. David. Weighted maximum likelihood
autoregressive and moving average spectrum modeling.
ICASSP, pp.3761–3764, 2008.
M. Oudot et al. Estimation of the spectral envelopeof voiced sounds using a penalized likelihood approach.
IEEE Trans. on SAP, 9(5):469–481, 2001.

ⓒ 2013 Information Processing Society of Japan推定されたスペクトル包絡推定された基本周波数 (最も優勢なカーネル)

有声・無声の度合い
2s

3s

4s

図 3 女性話声に対する推定結果

り，基本周波数の初期化を工夫する必要がある．

5. おわりに再構成スペクトログラム

1s

図 2 男性歌唱に対する推定結果

するには，基本周波数の時間的なダイナミクスを考慮した観測スペクトログラム

[17]

[18]

[19]

[20]
[21]

F. Villavicencio et al. Improving LPC spectral envelope
extraction of voiced speech by true-envelope estimation.
ICASSP, pp.869–872, 2006.
H. Kawahara et al. Restructuring speech representations
using a pitch-adaptive time-frequency smoothing and an
instantaneous-frequency-based F0 extraction: Possible
role of a repetitive structure in sounds. Speech Communication, 27(3–4):187–207, 1999.
T. Nakano and M. Goto. A spectral envelope estimation
method based on F0. SAPA-SCALE, pp.11–16, 2012.
A. Sasou and K. Tanaka. Robust LP analysis using glottal source HMM with application to high-pitched andnoise corrupted speech. Eurospeech, 2001.
T. Toda and K. Tokuda. Statistical approach to vocal
tract transfer function estimation based on factor analyzed trajectory HMM. ICASSP, pp.3925–3928, 2008.
亀岡弘和 et al. マルチカーネル線形予測モデルによる音声分析. 日本音響学会春季研究発表会, 2010.
H. Kameoka et al. Speech spectrum modeling for joint estimation of spectral envelope and fundamental frequency.
IEEE Trans. on ASLP, 18(6):1507–1516, 2010.
G. Lanckriet et al. Learning the kernel matrix withsemideﬁnite programming. JMLR, 5:27–72, 2004.
C. E. Rasmussen and C. K. I. Williams, editors. Gaussian Processes for Machine Learning. MIT Press, 2006.
H. Sawada et al. Eﬃcient algorithms for multichannel
extensions of Itakura-Saito nonnegative matrix factorization. ICASSP, pp.261–264, 2012.
M. Hoﬀman et al. Bayesian nonparametric matrix factorization for recorded music. ICML, pp.439–446, 2010.
M. Goto et al. RWC music database: Popular, classical,
and jazz music database. ISMIR, pp.287–288, 2002.
A. Kuramatsu et al. ATR Japanese speech database asa tool of speech recognition and synthesis. Speech Communication, 9(4):357–363, 1990.
N. Bertin et al. A tempering approach for Itakura-Saito
non-negative matrix factorization. with application tomusic transcription. ICASSP, pp.1545–1548, 2009.
D. H. Klatt and L. C. Klatt. Analysis, synthesis andperception of voice quality variations among female andmale talkers. JASA, 87(2):820–857, 1990.
G. Fant et al. A four-parameter model of glottal ﬂow.
STL-QPSR, 26(4):1–13, 1985.
C. Févotte et al. Nonnegative matrix factorization with
the Itakura-Saito divergence: With application to music
analysis. Neural Computation, 21(3):793–830, 2009.

6