1-Q-25
Convolutive HMM に基づく歌声の基本周波数制御モデルの提案とそのパラメータ最尤推定∗◎大石 康智 (名大), 亀岡 弘和, 柏野 邦夫 (NTT CS 研), 武田 一哉 (名大)1

はじめに歌声の基本周波数 (F0) 軌跡は，階段状の旋律概形に様々な動的変動成分が複雑に重ね合わさった状態で観測され，この変動成分は，歌声の個人性知覚に影響を与えるという知見が報告されている [1]．つまり，旋律に相当する譜面情報に加え，譜面には載っていない歌唱者の個性や感情のような非譜面情報が，F0軌跡に含まれていると考えられる．本研究では，文献 [2] の F0 制御モデルを参考にして，F0 軌跡から旋律概形を表す階段状の入力信号と動的変動因子を表すインパルス応答を同時推定 (分離) することを考える．ただ，譜面情報無しの状態で，観測される F0 軌跡だけから，それぞれの信号を推定する逆問題はそもそも不良設定問題であるため，従来のモデルだけではこの問題を解くことができない．そこで，我々は階段状の拘束をもつ特殊な入力信号を隠れマルコフ
モデル (HMM) でモデル化し，システムの伝達関数を全極モデルで表現することによる制約条件を設けた上で，Viterbi 学習と線形予測分析 (LPC) 的な解法の反復により，モデルパラメータを効率的に推定するアルゴリズムを提案する．その定式化と実装を行い，
F0 軌跡から旋律概形と動的変動成分をともに推定できること，さらに推定されたモデルパラメータによって，F0 軌跡を生成可能であることを確認する．
歌声の F0 軌跡から旋律概形と動的変動成分を分離する技術は様々な応用に貢献できる．例えば，歌声素片の連結や HMM を利用した歌声合成手法 [3, 4] によって，表情豊かな歌声合成音が実現されているが，知覚に影響を与える動的変動成分を精緻にモデル化する提案手法によって，さらに人間らしい合成音に近づけると考えられる．その他，歌唱スタイルの転写や歌唱力評価などへの応用が考えられる．

2数波
周u
時間

アイデア1

アイデア2

数波
周

HMMのパラメータ

U (z)

=

ここで，状態系列

y時間

アイデア3a0 +∑ az
Mj =1−jj

Y ( z ) = Z {y n }, U ( z ) = Z {u n }

aj

F0軌跡

インパルス応答j

最尤推定に基づくパラメータ最適化

ステップ1 Viterbi学習
HMMのパラメータ，状態系列の推定

反復推定

ステップ2 LPC的解法(a0 ,K, a ) の推定
M

Fig. 1 提案する歌声の F0 制御モデルの概略図
ここで，階段状の入力信号 un と観測される F0 軌跡
yn の z 変換を U (z)，Y (z)，伝達関数の次数を M と表す．音声分析の用途においては，零平均の Gauss 性雑音の入力信号を仮定し，フィルタと入力のスケールの任意性を除く目的で通常 a0 = 1 と置くが，これに対し，ここで我々が考えたい F0 制御モデルには，動的変動成分のスケールを入力音高信号に依らず自由に調節できる機構が備わっていることが望ましい．
そこで，a0 をパラメータ化することで，ビブラート等の動的変動成分のスケールを入力信号とは独立に調節できるようにする．すなわち，音声分析の用途では通常定数として扱われる a0 を，自由パラメータとして扱うことがここでは本質的に重要な意味をもつ．
式 (1) を逆 z 変換すると，

a0 yn +MX

aj yn−j = un

(2)

j=1

提案する歌声の F0 制御モデルは，Fig. 1 に示すような 3 つのアイデアに基づくモデルである．
アイデア 1 旋 律 概 形 を 表 す 階 段 状 の 入 力 信 号 は ，HMM によってモデル化される．
アイデア 2 F0 制御モデルは，全極モデルで表される伝達関数 (式 (1)) によってモデル化される．
アイデア 3 各モデルパラメータは，最尤推定に基づ
く以下の 2 つのステップの反復推定 によって最適化される．
ステップ 1 状態系列と HMM のパラメータはViterbi 学習によって推定される．
ステップ 2 インパルス応答 a0 , . . . , aM は，LPC的な解法によって推定される．
以下の節で各アイデアの詳細を述べる．

2.1 動的変動因子を表すインパルス応答のモデル化
F0 制御モデルの動的特性を表すインパルス応答は，以下の伝達関数によってモデル化される．

∗F0制御
F0制御モデル制御モデル

数波 減衰，定常振動
時間周全極モデルによる
HMMを利用した階段状の入力信号の 伝達関数のモデル化
モデル化 Y ( z )1

歌声の F0 制御モデルの提案

1Y (z)=PMU (z)a0 + j=1 aj z −j＊音高目標））
旋律概形（（音高目標
旋律概形(1)となり，入出力関係を表す差分方程式が導かれる．いま，F0 軌跡と入力信号の関係はこの差分方程式に従うと想定するが，実際にはこれにさらに微細変動成分 [2] が加わる．この微細変動成分を平均 0，分散 σ 2の正規分布に従う Gauss 性白色雑音 ε とし，ûn =PM
j=0 aj yn−j と un との間に，ûn − un = ε,(ε ∼ N (0, σ 2 ))(3)を仮定する．

2.2 階段状の入力信号のモデル化
階段状の拘束をもつ入力信号 un は，Fig. 2 に示す状態集合 S = {S1 , . . . , SI } からなる HMM によってモデル化される．1 回の状態遷移により実数が 1 個生成されるモデルであり，観測される F0 軌跡から状態遷移 (どの音高間でどのタイミングで遷移したか) を一意に決定できない．つまり，
un = mqn , (mSi ∈ R, qn ∈ {S1 , . . . , SI })(4)と表し，ここで mSi は状態 Si における出力確率分布(式 (3) で仮定する正規分布に相当する．) の平均を表す．一様な Markov 連鎖を想定し，状態 Sj から Si への遷移確率は P(Si |Sj ) と表す．qn は状態集合 S の要
素 (HMM の状態番号) の中のいずれかの値をとる．

Maximum Likelihood Parameter Estimation of F0 Control Model for Singing Voices based on Convolutive
HMM. by Yasunori OHISHI, Hirokazu KAMEOKA, Kunio KASHINO, Kazuya TAKEDA

日本音響学会講演論文集- 379 -2008年9月

となる．ここで各フレームの F0 軌跡は独立に観2測 さ れ る も の と 想 定 す る ．ま た ，σ (t) は フ レ ーム t に お い て 仮 定 す る 式 (3) の 正 規 分 布 の 分 散値に相当する．パラメータ集合 Θ の事後確率は，P(Θ|y(1) , . . . , y(T ) ) ∝ P(y(1) , . . . , y(T ) |Θ)P(Θ) なので，U (Θ) , P(y(1) , . . . , y(T ) |Θ)P(Θ) とおくと，cent ][P (S i +1 S i +1 )mSi+1qn+2 = Si+1

Si +1qn = Si

mS i

Si

qn+1 = Si

qn−2 = Si −1 qn−1 = Si−1

mSi−1
S i −1n−2n −1nn +1n + 2 [time]

Fig. 2 階段状の入力信号のモデル化
log U (Θ) =F0軌跡 yS3

入力信号 u

S2S1

フレーム tインパルス応答 a (jt )y (t ) , u (t )Fig. 3

j2.4 時変な F0 制御モデルへの拡張
Fig. 3 より，旋律の変化とともに，ビブラートやオーバーシュートのような動的変動成分が次々と変化することがわかる．これは，図 1 の F0 制御モ
デルのインパルス応答が，時々刻々と変化するためであると考えられる．そこで，長さ N の F0 軌跡
y = {y1 , . . . , yN } と入力信号 u = {mq1 , . . . , mqN }をフレームに分割し，フレームごとにインパルス応答を推定する．長さ L のフレーム t における F0(t)(t)
軌跡と入力信号をそれぞれ y(t) = {y1 , . . . , yL }，u(t) = {mq(t) , . . . , mq(t) }，インパルス応答を θ (t) =L(t){a0 , . . . , aM } と表す．ここで入力信号を生成する
HMM のパラメータ mS1 , . . . , mSI はすべてのフレーム (フレームの総数は T ) にわたって共通するものと考える．この状態系列と HMM のパラメータを合わせて，ω = {q1 , . . . , qN , mS1 , . . . , mSI } と表すと，推定しなければならないパラメータの集合は，
Θ = {θ (1) , . . . , θ (T ) , ω} となる．

F0 制御モデルのパラメータ最尤推定

3

F0 軌跡の集合 {y(1) , . . . , y(T ) } が与えられたとき，パラメータ集合 Θ を最尤推定する方法を提案する．
{y(1) , . . . , y(T ) } の対数尤度は，式 (3) より，log P(y(1) , . . . , y(T ) |Θ) =TXT nXJ ≡ −(L − M )log P(y(t) |θ (t) , ω)−2σ (t)2

l=M +1

j=0l

LMN³X´2 XX(t) (t)aj yl−j − mq(t) +Pqn ,qn−11

t=1 2σ(t) 2

l=M +1lj=0n=1がパラメータの集合 Θ に関して最大化したい目的関数である．しかし，式 (8) を最大化する Θ は解析的に求めることができない．そこで，Θ の各要素に関して，他の要素を固定したもとで式 (8) を最大化する
ステップを式 (8) の値が収束するまで繰り返す (Fig.
4)．

3.1 インパルス応答 θ (t) の推定
パラメータ集合 ω を固定して，式 (8) を最大にす(t)る θ (t) を求める．J を aj ′ , (j ′ = 0, . . . , M ) に関して偏微分して 0 とおくと，j ′ = 0 のとき(t)a0LM³X´X2(t) (t)(t)aj yl−j − mq(t) yl + σ (t) (L−M ) = 0
l=M +1lj=0(9)
j ′ = 1, . . . , M のとき
LM³X´
X(t) (t)(t)
aj yl−j − mq(t) yl−j ′ = 0(10)lj=0となる．J を σ (t) に関して偏微分して 0 とおくと，2(5)日本音響学会講演論文集

TX

σ (t) =LM
³X
´2 o
X(t) (t)aj yl−j − mq(t)(t)log σ (t) a0(8)

t=11

TXt=1

√
(t)
− (L − M ) log( 2πσ (t) a0 )−t=1log P(q1 , . . . , qN ) = log P(q1 )P(q2 |q1 ) . . . P(qN |qN −1 )(7)となる．遷移確率 P(Si |Sj ) は事前に決定する定数である．簡単のため，PSi ,Sj ≡ log P(Si |Sj ) とおく．以上より，式 (5) と式 (7) を式 (6) に代入し，定数項を除いた以下の式

l=M +1t=1=log P(θ (t) )と 記 述 で き る ．こ こ で 事 前 確 率 P(θ (t) ) とP(mS1 , . . . , mSI ) は一様分布を想定し，P(q1 , . . . , qN )は，一様な Markov 連鎖を想定するため，

j時変な F0 制御モデルへの拡張

1T
X+ log P(mS1 , . . . , mSI ) + log P(q1 , . . . , qN )(6)

2.3 各モデルパラメータの解釈
インパルス応答 a0 , . . . , aM は，歌声の動的変動因子を表すパラメータである．音高が安定するときの振動や，音高目標に到るまでの連続的な音高遷移が表現される．状態系列 q1 , . . . , qN は，音長を決定するパラメータである．これは，必ずしも譜面に記される音符の音長に対応するわけではなく，歌唱者の意図やスタイルに基づいて生成される運動指令の長さを表現するものと想定する．最後に，HMM のパラ
メータ mS1 , . . . , mSI は，歌唱者が意図する旋律の音高 (音高目標値) に対応するパラメータである．これは，必ずしも譜面に記される音符の音高 (客観的に定まっている音高値) には対応しない．

(t)log P(y(t) |θ (t) , ω) +t=1

ビブラート
オーバーシュート
プレパレーションa (jt +1 )フレーム t+1 y (t +1) , u (t +1)

TX

1L−MLM³X´2
X(t) (t)aj yl−j − mq(t)(11)
l=M +1j=0lとなる．式 (5) の正規化項に動的変動成分のスケー(t)ル a0 を含むため，j ′ が 0 とそれ以外で，偏微分の式が異なる．したがって，以上の連立方程式を解析- 380 -

2008年9月

START

初期値 q ,K , q の設定
初期値 m ,K, m の設定
θ( ) ,K, θ( ) の推定
q , K , q の推定
m , K , m の推定Noinit
1init
NinitS1init
SI1T

1N

S1SI

a0 の推定(t )K

a1(t ) , , aM(t ) の推定σ (t ) の推定
2

J の収束性の判定

J の収束性の判定

Fig. 4K

θ(t ) の推定 (t = 1, , T )

NoYes

Yes
RETURN
反復法によるモデルパラメータの推定手順的に解くことはできない．そこで，まず，式 (9) の(t)(t)(t)a1 , . . . , aM , σ (t) を固定して，a0 に関する 2 次方(t)程式を解析的に解く．次に，式 (10) の a0 を固定し(t)(t)て，a1 , . . . , aM に関する連立方程式を解析的に解
2
く．最後に式 (11) の σ (t) を計算する．この計算過程を目的関数 J が収束するまで順番に繰り返して，(t)(t)a0 , . . . , aM , σ (t) を更新させる．

3.2 状態系列 q1 , . . . , qN の推定
インパルス応答の集合 {θ (1) , . . . , θ (T ) } と HMMのパラメータ mS1 , . . . , mSI を固定して，状態系列
q1 , . . . , qN に関して，式 (8) の最大化を考える．この問題は，Viterbi アルゴリズム (動的計画法) により効率的に解ける．最初から時刻 k に状態 Si に至るまでの最適な状態系列は漸化式を利用して，hδk (Si ) = max δk−1 (Sj )Sj−X

1(t,l)∈Ck 2σ(t)(t)22 (ûl − mSi ) + PSi ,Sji(Ck = {(t, l)|yk ∈ y(t) , 1 ≤ t ≤ T, 1 ≤ l ≤ L})(12)と記述できる．ここで，Ck は，1 から T までのフレー(t)ムの中で，yk に相当する yl のフレーム番号 t とインデクス l の組み合わせからなる集合である．式 (12)を k = N まで計算すれば，最適な経路が求められる．

3.3 HMM のパラメータの推定
インパルス応答の集合 {θ (1) , . . . , θ (T ) } と状態系列
q1 , . . . , qN を固定して，式 (8) が最大になるように
HMM のパラメータ mS1 , . . . , mSI を推定する．J をmSi に関して偏微分して，0 とおくと，X (t)1 X 1ûl(13)mSi =|Di ||Cn |n∈Di(t,l)∈Cnを得る．ここで，集合 Di = {n|qn = Si } とし，|Di | をその要素数とする．最終的に，3.2 節と本節で推定される状態系列 q1 , . . . , qN と HMM のパラ
メータ mS1 , . . . , mSI から，階段状の入力信号 un =mqn , (n = 1, 2, . . . , N ) が決定される．
以上のように階段状の入力信号を HMM によってモデル化し，そのパラメータ学習方法が入力信号と F0軌跡に基づく逆畳み込みと Viterbi 学習から成ることから，提案モデルを畳み込み HMM(ConvolutiveHMM) に基づく歌声の F0 制御モデルと呼ぶ．

日本音響学会講演論文集

3.4 初期値の設定局所解に収束することを防ぐために，2 つの初期値設定を行う (Fig. 4)．まず，観測系列 y に 3.2 節
で述べた Viterbi アルゴリズムを適用して初期状態init
系列 q1init , . . . , qNを決定する．次に，3.3 節と同様に，初期状態系列から初期の HMM のパラメータ
initminit
S1 , . . . , mSI を決定する．
3.5 F0 軌跡の生成
フレーム t において推定されたインパルス応答(t)(t)θ (t) と状態系列 q1 , . . . , qL ，HMM のパラメータ
mS1 , . . . , mSI から，F0 軌跡を生成する．式 (2) より，フレーム t において生成される F0 軌跡 ŷ(t) は，
1 ¡ (t) X (t) (t) ¢ul −aj ŷl−j , (l ≥ M + 1) (14)(t)
a0j=1M(t)ŷl=(t)となる．ただし，ŷl4(t)= yl , (l ≤ M ) とする．

評価実験

畳み込み HMM に基づく F0 制御モデルのパラメー
タ学習アルゴリズムによって，各フレームにおいて推定される階段状の入力信号 u(t) とインパルス応答
θ (t) ，それらから生成される F0 軌跡 ŷ(t) の性能を評価する．

4.1 歌声データベース
クラシックの声楽家，プロのポップス歌手，音楽的な訓練を受けていない素人 (それぞれ男女 1 名ずつの計 6 名) の歌声からなるデータベースを構築した．歌
唱者は，ヘッドフォンで歌唱曲の旋律 (ガイドトーン)を聴きながら伴奏なしで歌唱した．歌唱曲は，「きらきら星」，「喜びの歌 (Beethoven の交響曲第 9 番第 4楽章の歌の部分を岩佐東一郎氏によって作詞されたもの)」である．日本語歌詞による歌唱，ハミングによる歌唱の 2 パターンを収録した．
4.2 実験条件
F0 は，de Cheveigné らの提案した YIN[5] を利用して 10ms ごとに推定される．なお，Hz で表される周波数 yHz を，cent で表される対数スケールの周波数 ycent に変換する．200ms 以内の F0 が推定されない区間は，無声音の区間とみなし，F0 を線形補間する．また，ガイドトーンからも F0 を推定し，提案手法により推定される階段状の入力信号の評価のために利用する．
推定された F0 軌跡は，500ms ずつシフトさせて得られる長さ N の素片に分割される．本実験ではN = 2s とする．この素片が 2.4 節で述べた y に相当し，素片ごとに提案手法を適用する．素片に F0 が推定されない区間 (休符区間) が含まれる場合は，今回は分析を行わない．各素片におけるフレームは 100msずつシフトさせる．また，各素片の HMM の状態数 Iは相平面を利用して決定する．相平面に描かれる F0軌跡のアトラクタの分布を混合ガウス分布 (GMM) によって学習し，その極大値の数を HMM の状態数 Iとする [6]．各状態の初期状態確率は 1/I ，自己遷移確率は 0.9，状態間の遷移確率は 0.1/(I − 1) とした．
各フレームにおいて仮定する正規分布の分散の初期値は 2500 とした．
4.3 フレーム長と伝達関数の次数の評価ある素片の推定結果を Fig. 5 に示す．フレーム長
L は 250ms ，伝達関数の次数 M は 3 とした．音高遷
移の部分を除いて，入力信号 u は対応する区間のガイドトーンの F0 軌跡 g と近いものが推定される．し- 381 -

2008年9月

[cent ]6000%][率解正階段状の入力信号 u
ガイドトーン g生成されるF0軌跡 ŷ ( )][%率解正]%[率解正][%率解正t5800

5600

観測されるF0軌跡 y
02 [sec ]1数波
周

伝達関数の振幅応答zH5Fig. 5 推定結果の例
かし，F0 軌跡 y(t) と生成される F0 軌跡 ŷ(t) との間に誤差が生じる箇所がいくつか存在する．素片によっては ŷ(t) が発散してしまうフレームも存在した．
以上を踏まえて，L と M を変化させたときの各フレームの u(t) と ŷ(t) の推定性能を評価した．フレー
ムごとに以下の平均二乗誤差 (RMS)v
uLX
u 1(t)(t)RMSu(t) = t(ul − gl )2
L−M
l=M +1(15)v
uLX
u 1(t)(t)(ŷl − yl )2RMSŷ(t) = t
L−M

65.6

ŷ (t ) の評価

u (t )の評価

83.5

83.7

62.2

M =327.2

M =683.2

29.3

M =383.7

83.2

83.7

2.4
M =34.4
M =653.9

14.9

ms

84.7

M = 12L: 750

M =910.6
M =683.9

M = 12
L: 50044.0 84.0

83.8

10.1
M =369.3

M =9

M =683.2

フレーム長 L: 250ms
68.8 84.1

ms

84.0

19.984.2

M = 12
L: 1

M =9s

84.0

84.7
7.9
M = 12

5.7
M =9Fig. 6 フレーム長と伝達関数の次数に関する評価
0.2

クラシックの声楽家
（女性）

度頻
対相

000.2

2

000.2

2

4

4

度頻
対相

002

4

クラシックの声楽家
（男性）

000.2

2

6 [Hz]00

2

素人 （女性）0.2

6 [Hz]

00

6 [Hz]

プロのポップス歌手（女性）

度頻対相0.2

4

6 [Hz]

プロのポップス歌手（男性）
4

6 [Hz]

素人 （男性）24

6 [Hz]

l=M +1を計算し，分析を行ったすべてのフレームのうち正しく推定された割合を正解率として以下に定義する．

RMSu(t) < ξ となるフレーム数
× 100
すべての素片の総フレーム数
RMSŷ(t) < 10 となるフレーム数
正解率ŷ(t) =× 100
すべての素片の総フレーム数(16)

正解率u(t) =ここで，ξ は各素片の y と g の RMS とする．ガイド
トーンを聴きながら歌唱しても，歌唱者によってはガイドトーンからずれて歌唱するため，ξ 以内で u(t) が推定されれば正解とした．ŷ(t) に関しては，RMSŷ(t)が 10cent 以内で生成されれば正解とした．
推定性能を Fig. 6 に示す．L を短くするにつれて，ŷ(t) の正解率が向上した．つまり，時変な F0 制御モデルとして，素片をフレームに分割し，フレームごとに伝達関数を推定することの有効性を確認できた．
また，M の増加によって伝達関数の自由度が増すた
め，正解率が向上した．しかし，総フレーム数の 3 割は，伝達関数が適切に推定されないため，F0 軌跡を生成させることが困難であった．また，L を短く，Mを大きくしすぎると，素片内で扱うモデルパラメータ数が増えてしまう問題もある．したがって，情報量基準などを利用して L や M を変化させながら推定を行うなどの提案手法の改善が必要とされる．
伝達関数の振幅応答に基づく歌唱者の個性に関する考察
Fig. 5 の下段に式 (1) から計算される伝達関数の振幅応答を示す．音高が安定するフレームでは共振周波数をもち，そのいくつかはビブラート (4〜7Hz 程度の周期的変動 [1]) に対応すると考えられる．そこで，
式 (16) で正解とされたフレームの共振周波数の分布を Fig. 7 に示す．歌唱者ごとに分布の形状が異なり，特に女性の声楽家は 2.5Hz あたりに分布が集中した．
この歌手は音高が安定するときに顕著に振動させる歌い方をすることを確認した．一方，素人の歌唱者ほど平坦な分布となり，推定された伝達関数の振幅応答から個性や歌唱力の違いを確認することができた．

4.4

日本音響学会講演論文集

Fig. 7

5

歌唱者ごとの共振周波数の分布

まとめと今後の展開

歌声の F0 軌跡から，階段状の入力信号と動的変動因子を表すインパルス応答を同時推定するためのF0 制御モデルとパラメータ学習アルゴリズムを提案した．実験結果より，分析を行った総フレーム数の69.3%は正確にそれぞれの信号を推定できることを確認した．また，推定される伝達関数の振幅応答に，歌唱者の個性や歌い方の違いが表れることを確認した．
性能を向上させるためには，フレーム長や伝達関数の次数を可変にしてパラメータを推定することが必要とされる．また，F0 軌跡の符号化，歌唱者の歌い方を反映した歌声合成，歌唱力評価など本研究の応用先について検討する必要がある．楽器音や生体信号などの時系列信号への適用や，提案モデルを多変量化し，MFCC ベクトルの時系列などを動的モデリングして音声分析に利用することも今後の展開である．
謝辞 本研究は日本学術振興会特別研究員 (DC2) 科研費の補助を受けた．また，これまで本研究に対し，有益なご助言を頂いた後藤真孝氏 (産総研)，伊藤克亘氏 (法政大) に感謝致します．また，
F0 制御モデルに関して有益なご意見を頂いた南泰浩氏 (NTTCS研) に感謝致します．

参考文献[1] 齋藤 毅ほか：歌声の個人性知覚に寄与する音響特徴の検討，
音講論集, 2-Q-26，pp. 601–602 (2007).
[2] 齋藤 毅ほか：SingBySpeaking:歌声知覚に重要な音響特徴を制御して話声を歌声に変換するシステム，情処研報音楽情報科学，Vol. 2008, No. 74, pp. 25–32 (2008).
[3] Bonada, J. et al.: Synthesis of the Singing Voice by Performance Sampling and Spectral Models, IEEE Signal
Processing Magazine, Vol. 24, No. 2, pp. 67–79 (2007).
[4] 酒向慎司ほか：声質と歌唱スタイルを自動学習可能な歌声合成システム，情処研報音楽情報科学，Vol. 2008, No. 74, pp.
51–58 (2008).
[5] de Cheveigné, A. et al.: YIN, a fundamental frequency
estimator for speech and music, Journal of the AcousticalSociety of America, Vol. 111, No. 4, pp. 1917–1930 (2002).
[6] 大石康智ほか：相平面を利用した歌声の F0 軌跡の新しい表現方法，電子情報通信学会総合大会，pp. S–51–S–52 (2008).

- 382 -2008年9月