調波時間構造化クラスタリングによる CASA へのアプローチ
亀岡 弘和†

ルルー・ジョナトン†

小野

順貴†

嵯峨山茂樹†† 東京大学大学院情報理工学系研究科
E-mail: †{kameoka,leroux,onono,sagayama}@hil.t.u-tokyo.ac.jp
あらまし

人間の聴覚機能を計算機で実現しようという試みが活発に進められており，その枠組を総称して計算論的

聴覚情景分析 (Computational Auditory Scene Analysis: CASA) と呼ぶ。近年この研究分野における興味の対象は，Bregman が指摘した分凝要件 [1] に基づく混合音分離法の実現にある。CASA における多くの従来手法の間で共通するのは，各時刻で独立に調波成分を見つけ出すための処理 (周波数方向の群化) と，抽出された調波成分特徴量の時系
列を時間方向にスムージングする処理 (時間方向の群化) を多段処理的に行っている点である。しかしながら，より良い群化プロセスの実践のためには，これらは本来協調し合いながら行われるべきであり，個々の音源の時間周波数全域に渡ったスペクトル構造を一挙に推定できる方法論が不可欠であると我々は考える。本稿では，この観点から導かれる調波時間構造化クラスタリングと呼ぶ CASA のための新しいアプローチを提案する。
キーワード

計算論的聴覚情景分析, Bregman の分凝要件, 定 Q フィルタバンク, 調波時間構造化クラスタリング

Harmonic-Temporal Clustering:
A New Approach to CASA
Hirokazu KAMEOKA† , Jonathan LE ROUX† , Nobutaka ONO† , and Shigeki SAGAYAMA†† Graduate School of Information Science and Technology, The University of Tokyo
E-mail: †{kameoka,leroux,onono,sagayama}@hil.t.u-tokyo.ac.jp
Abstract There have been many studies on computational implementation of human auditory function. Sucha framework is referred to as computational auditory scene analysis (CASA). The recent interest in this area of
research is focused on the development of source separation algorithms based on the grouping cues suggested by
Bregman [1]. In most of the conventional methods for CASA, stream segregation is performed through two successive procedures: identification of the harmonic components and temporal smoothing of the extracted harmoniccomponent features. These procedures, however, should be done essentially in a cooperative way and this belief has
led us to formulate a unified estimation framework for the two dimensional structure of time-frequency components.
We propose in this paper a new approach for CASA called Harmonic Temporal Structured Clustering.
Key words CASA, Bregman’s grouping cues, constant Q filter bank, Harmonic Temporal Structured Clustering

1. 序

論

様々な音が混在する中で目的音の情報を計算機にうまく分離推定させることは，ロボット聴覚，音声認識や自動採譜などの様々な魅力的なアプリケーションを実現する上で工学的に重要な要素技術である一方で，数理的には大変難しい問題である。しかし，我々人間は，そのような環境の中でも容易に特定の音だけを選択的に聴取することができることから，聴覚は能動的に外界を把握するための優れた機能をもつと考えられている。このように音を通じて環境を把握することは聴覚情景分析 (Auditory Scene Analysis: ASA) と呼ばれ，Bregmanの著書 [1] により広く知られることとなった。この著書の中でBregman は，知識を利用しない聴覚の低次の分離能力に関して，

I. 音響信号はスペクトログラムに似た要素に「分解」される
II. 同じ音源に由来する要素は「群化」されて音脈を形成する
III. 群化のされやすさ (分凝要件) は，(a) 周期性 (調波構造)，(b) 調波成分の共通の立上り時刻，(c) 調波成分の共通の周波数および振幅変化，(d) 成分の連続性，(e) 時間周波数の近接性，(f) 共通の音源位置などに関係することを心理実験を通して示している [1]。近年このような聴覚機能を計算機で実現しようという試みが活発に進められており，その枠組を総称して計算論的聴覚情景分析 (Computational
Auditory Scene Analysis: CASA) と呼ぶ。
最近の CASA の研究における興味の対象は，以上の分凝要件に基づく混合音分離法の実現にある。すなわち，分解と群化のプロセスを，分凝要件に関係する物理量を用いてアルゴリズ—1—ムとして実現し，音脈の認識に有用な特徴量 (ピッチ周波数なする時変スペクトル構造を具体的にモデル化する。2. 5 節でモど) を抽出したり，目的音に相当する音脈の再構成を行うこと

デルパラメータを最適化する反復アルゴリズムを導出する。

が主な目的である。

2. 2 擬似周期信号の定 Q フィルタバンク出力の導出
音源 k の信号モデルとして n 次調波成分の瞬時位相が nθk (u)+ϕk,n ，瞬時振幅が w
ek,n (u) の擬似周期信号の解析信号表現

このような背景のもと，Cooke [2]，Brown ら [3]，Ellis [4]，Fishbach [5]，中谷ら [6] は分凝要件に基づく混合音分離の手法を提案している。これらの手法の多くは規則記述的や人工知能的なアプローチをとるため，トップダウン式にさまざまな制約条件を導入しやすい点では優れるが，細分化された閾値判定処理で構成される煩雑なアルゴリズムにならざるを得ないのが難点である。一方，西ら [7]，鵜木ら [8]，安部ら [9], [10]，Wu
ら [11] は，分凝要件を数学的な制約として定義し，それに基づ
く制約つき最適化問題として定式化を行っている。これらに代表される CASA の多くの従来手法の間で共通するのは，時間周波数分解による分解プロセスを経て，群化プロセスが，瞬時特徴量の尤度を計算するステップと，その尤度に基づいて瞬時
特徴量の時間軌跡を Kalman フィルタなどのダイナミカルシス
テム，隠れマルコフモデル，マルチエージェントシステムなどを用いて逐次的に追跡するステップの多段処理により構成される点にある。上記の群化プロセスのうちの前半のステップは，各時刻で調波構造をなしている周波数成分を見つけ出すための処理であり，隣接する時刻と独立な周波数方向の群化に相当する。そして，後半のステップは，抽出された成分あるいはそれに相当する特徴量の時系列を時間方向にスムージングする処理であり，時間方向の群化に相当する。
しかしながら，周波数方向と時間方向の群化を順次行うこと
が群化プロセスの最適な実践方法であるとは必ずしも言えない。なぜなら，周波数方向の群化の精度は隣接する時刻間での成分の連続性を考慮することではるかに向上されうるし，時間方向の群化の精度もまた周波数方向の群化が高精度であればあるほど高くなるはずであることから，まさに両者は，以上の意
味で鶏と卵の関係にあるからである。我々は以上の問題意識のもと，より良い群化プロセスの実践のためには，個々の音源の時間周波数全域に渡ったスペクトル構造を一挙に推定できる方法論が不可欠であると考え，この観点から問題の定式化を目指した。本稿では，調波時間構造化クラスタリングと呼ぶ CASAのための新しいアプローチを提案する。

2. 調波時間構造化クラスタリングの定式化
2. 1 方法論の概要と本章の構成
入力された音響信号の時間周波数成分を音脈に対応する時間
周波数成分に分解する問題を定式化するにあたり，まず知覚上の単位である音脈の概念を計測可能な物理量として具現化する必要があるが，ここで言う知覚と物理量の橋渡し役を担うのがまさに Bregman の分凝要件である。そして，この分凝要件か
ら逸脱しない範囲の自由度をもった時変スペクトルを直接的にモデル化し，これを混合したもので対象の時変スペクトルを説明しようとすることが，本稿で提案する調波時間構造化クラス
タリングのアイディアの要点である。すなわち，分解された要素 (時間周波数成分) を音脈に対応する時変スペクトル成分に文字通りクラスタ化することが狙いである。
定 Q フィルタバンクが聴覚末梢系を良く模擬したモデルであることが知られることから，次節で，調波構造と共通周波数変化をなす擬似周期信号モデルの定 Q フィルタバンク出力を導出し，続く 2. 3 節ではさらに共通の立上りと共通振幅変化の制
約，2. 4 節では成分の連続性の仮定を新たに設け，音脈に対応

fk (u) =N
Xw
ek,n (u)ej (nθk (u)+ϕk,n )

(1)

n=1

を考え，この信号の定 Q フィルタバンク出力を導出する。u は時刻を表す。まず，wavelet 基底関数を

µ

u−t
1
ψ
ψa,t (u) , √
a
2π|a|

¶
(2)

と定義する。ただし，a はスケールパラメータ，t はシフトパラ
メータ，ψ(u) はアドミッシブル条件を満たす周波数 1 の任意のアナライジング wavelet である。fk (u) の連続 wavelet 変換は

¡

D

¢

E

Wk log a1 , t , fk (u), ψa,t (u)
=

Z ∞X
N

u∈R

∗
(u)du
w
ek,n (u)ej (nθk (u)+ϕk,n ) ψa,t(3)−∞ n=1∗で与えられるが，ここで，一般に ψa,t(u) は時刻 t に局在すること，瞬時位相 θk (u) と調波成分の瞬時振幅 w
ek,n (u) は緩やかに変化することを踏まえ，時刻 t 周辺で θk (u) と w
ek,n (u) を¡¢θk (u) ' θk (t) + θk0 (t) u − tw
ek,n (u) ' w
ek,n (t),(4)と 0 次および 1 次近似する。瞬時位相の導関数は瞬時周波数なので，θk0 (t) は瞬時基本周波数を表し，これを µk (t) と置くことにすると，式 (4) より式 (3) は，¡¢

Wk log a1 , t 'NXw
ek,n (t)ej (nθk (t)+ϕk,n )Z ∞∗ejnµk (t)(u−t) ψa,t
(u)du
−∞

n=1

と近似できる。また，一般化 Parseval の等式より

Z ∞∗ejnµk (t)(u−t) ψa,t
(u)du
−∞

µ

¿
ejnµk (t)(u−t) , √

=

¿
=

√

¡

u−t
1
ψ
a
2π|a|

¶À

¢

À

u∈R

¡
¢
1
= Ψ∗ anµk (t)2πδ ω − nµk (t) , √ Ψ(aω)
2πω∈Rであるので，結局¡¢

Wk log a1 , t '

N
X

¡

¢

w
ek,n (t)Ψ∗ anµk (t) ej (nθk (t)+ϕk,n )n=1である。あとはこれを K 個分足したもの¡¢

W log a1 , t =KX¡

Wk log a1 , t

¢

(5)

k=1

が混合音信号の時刻 t における定 Q フィルタバンク出力となる。
ところで，a の次元は周期であるので，対数周波数 x = log a1—2—

に変数変換し，Ωk (t) , log µk (t) と置くと

XX
k

³

´

w
ek,n (t)Ψ∗ ne−x+Ωk (t) ej (nθk (t)+ϕk,n )Ψ(ω)

W (x, t) =

1

(6)

0.5

n

と書ける。ここで，周波数特性 Ψ(ω) が，次のような ω = 1 で最大値をとる単峰的な実関数 (図 1 参照)

0

0
0

図1

(7)

(ω 5 0)となるようなアナライジング wavelet を選ぶと，式 (6) は，Ã ¡w
ek,n (t) exp−

x − Ωk (t) − log n

ω

6

8

10

周波数特性 Ψ(ω)(式 (7)) の σ = 21 のときの概形。

n=1

n=2
n=3

ej (nθk (t)+ϕk,n )Ω k+log2

Ωk

y
rg
e
n
e

¯
¯
¯W (x, t)¯2 =
¯
¯2
Ã ¡
¢2 !
¯
X ¯¯
x − Ωk (t) − log nj (nθk (t)+ϕk,n ) ¯eek,n (t) exp −
¯w
¯
4σ 2
¯
¯
k,n
Ã ¡
¢2 !
XX
x − Ωk (t) − log nw
ek,n (t)w
ek0 ,n0 (t) exp −+2
Ã ¡x − Ωk0 (t) − log n0Ω k +log3 ...

log-freq.

y =0

y =1

y =2 y = 3
y =4

φk

φk

τk

4σ

| 0
k=k
| 0 n=n...

図 2 式 (8) の概形。

となる。以上より，式 (5) のパワースペクトルは

−

4

¢2 !

4σ 2

k,n

exp

2

n=4

W (x, t) =

X

1

power

Ã ¡

¢2 !
log
ω

exp −(ω > 0)4σ 2Ψ(ω) = Ψ∗ (ω) =φkφ k φk

y =5

φk

φk

φk

time

図 3 式 (11) の概形。

¢2 !
j(nθk (t)+n0 θk0 (t)+ϕk,n +ϕk0 ,n0 )

e

4σ 2

を満たすものとすると，パラメータ evk,n は音源 k の n 次調波成分のエネルギーと書けるが，ここで，上式第二項において k =| k0 または n =| n0の exp 項同士の重なりがほとんどない（注 1） と仮定すれば，第二R∞−∞

wk,n (t)dt に対応する。また，必ずしも必要ではないが，便宜的に evk,n , wk vk,n と分解し，vk,n は，P

nvk,n = 1 を満たすものとする。正規化共通瞬時振幅 uk (t)項は第一項に比べると無視できるほど小さいと見なせるため，√
wk,n (t) , 2πσ|w
ek,n (t)|2 と置くと，時刻 t における混合音信はノンパラメトリックに定めることもできるが，ここではパラ号モデルの定 Q フィルタバンク出力パワーは各 Gauss 分布関

束つき混合正規分布モデル (図 3 参照)

メトリックモデルを仮定することにし，uk (t) を次のような拘数が x = Ωk (t) + log n でピークをとるような混合正規分布と
同形の関数

¯
¯
¯W (x, t)¯2 'X wk,n (t)k,nX uk,y
√

y=0

Ã

√
exp
2πσ

¡
−

x − Ωk (t) − log n
2σ 2

µ

Y −1

uk (t) =

2πφk

exp −

(t − τk − yφk )2
2φk 2

P

¢2 !

とする。従って，式 (10) を満たすには，

(8)

y

¶
(11)

uk,y = 1 でありさえすれば良い。この選択は，非負制約や一定の滑らかさの制約を導入しやすくするだけでなく，のちに述べる効果的な最適化アルゴリズムの導出に大変都合が良い。τk は先頭の Gaussとなる (図 2 参照)。

2. 3 共通振幅変化の制約
音源 k のモデルにおいて，各調波成分の共通の立上りと共通の振幅変化を仮定することは，n 次調波成分の瞬時振幅を
wk,n (t) = e
vk,n uk (t)(9)

するパラメータと見なせる。式 (9)，(11) を式 (8) に代入すると，具体的にK X
N Y−1(x−Ωk (t)−log n)2(t−τk −yφk )2
XX−
wk vk,n uk,y −2σ 2
2φ2
k=1 n=1 y=0と置くことに相当する。ここで，uk (t) が

2πσφk

e

k

(12)

Z ∞
uk (t)dt = 1分布関数の中心を表すため，これは音源 k の立上り時刻に相当

(10)

−∞

（注 1）：この仮定は，音声スパース性と呼ばれる混合音声信号の時間周波数成分
はまばらにしか存在しないという性質に基づくものであり，多チャンネルブライを得る (図 4 参照)。

2. 4 成分の連続性の制約
未知の連続関数であるピッチ軌跡関数 Ωk (t) は，t の多項式
Ωk (t) , Ωk,0 + Ωk,1 t + Ωk,2 t2 + Ωk,3 t3 + · · ·(13)

ンド音源分離においては多用される [12]。ここでは，パワースペクトルの加法性の正当化のためにスパース性の仮定を利用している。

と，より安定な曲線を形成する区分的多項式 (3 次 Spline 関数)—3—Power density

I(Θ) 5 I + (Θ, m) =wk,n (t)ZZ µ X
time

¯
¯2
¯
¯2
mk,n,y (x, t)¯Y (x, t)¯¯¯mk,n,y (x, t) Y (x, t) log
Wk,n,y (x, t)

k,n,y

X³
−

log-frequency

図 4 分凝要件に基づく時変パワースペクトルモデルの概形。

¯¯2

mk,n,y (x, t)¯Y (x, t)¯ − Wk,n,y (x, t)

´¶
dxdt

(17)

k,n,y

のような凸不等式が成り立つ。ただし，Ωk (b)

Wk,n,y (x, t) ,Ωk,i-1

Ωk,i

b i-1

bi

Ωk,i+1b³−

1Ωk,i (bi+1 − t) + Ωk,i+1 (t − bi )bi+1 − bi´¤£1(t − bi )(bi+1 − t) (bi+2 − t)Ω00k,i + (t − bi−1 )Ω00k,i+1 ,6t ∈ [bi , bi+1 [(14)(図 5 参照) の両者を検討する。ただし，上記 3 次 Spline ピッチ軌跡関数において，時間軸上のノード b1 , · · · , bI は定数であり，各ノード上でのピッチ周波数値 Ωk,1 , · · · , Ωk,I がパラメータである。また，各ノード上でのピッチ軌跡関数の 2 次微分係00
数 Ω00
k,i は，Ωk,i ，Ωk,i をそれぞれ要素にもつベクトルを Ω，Ω00 とすると，Ω00 = M Ω で求められる。M は，b1 , · · · , bI が定数 (既知) であれば陽に決定される定数行列である。
Y −1
以上より，{Ωk,i }Ii=1 ，{vk,n }N
n=1 ，{uk,y }y=0 ，wk ，τk ，φkがすべて求まれば，Ωk (t) および wk,n (t) が決まるため，式 (1)より，音源 k の信号を再構成することができる。次節より，これらのパラメータを推定するための反復アルゴリズムを導く。
2. 5 パラメータ推定アルゴリズム
Y −1
K
パラメータ {{Ωk,i }Ii=1 ,{vk,n }N
n=1 ,{uk,y }y=0 ,wk ,τk ,φk }k=1をまとめて Θ と表記することとし，対象の信号を最も良く説明する Θ を求めることがここでの目標である。そのためにはまず，対象とする信号の定 Q フィルタバンク出力パワーを |Y (x, t)|2とすると，|Y (x, t)|2 と |W (x, t)|2 との近さを表す尺度を導入する必要がある。|Y (x, t)|2 と |W (x, t)|2 はいずれも非負関数であるが，非負関数間の近さを測る尺度として I ダイバージェ
ンスと呼ばれるものがあり，

ZZ Ã
I(Θ) ,(x−Ωk (t)−log n)2(t−τk −yφk )2
−
2σ 2
2φ2
k

であるが，これは具体的には時刻 τk + yφk 近傍にのみ分布
b i+1

図 5 3 次 Spline ピッチ軌跡関数 (式 (14))Ωk (t) ,wk vk,n uk,y −
e
2πσφk

¯
¯
¯Y (x, t)¯2¯
¯2
¯Y (x, t)¯ log ¯¯¯W (x, t)¯2

!
³¯
¯2 ¯
¯2 ´
− ¯Y (x, t)¯ − ¯W (x, t)¯
dxdt

(15)

する，n 次調波成分のパワースペクトルの一部を表す。このWk,n,y (x, t) を，音脈を構成する「要素」に対応するものと考え，以後サブクラスタモデルと呼ぶ。さて，式 (17) を注意深く見てみると，右辺は実は mk,n,y (x, t)|Y (x, t)|2 と Wk,n,y (x, t)との間の I ダイバージェンスの総和に他ならないことに気づく。CASA やブラインド音源分離の分野ではバイナリマスクを用いて時間周波数の各成分を排他的に 1 つの音源にのみ帰属させることによる音源分離法 [12] が知られているが，これに対し，
mk,n,y (x, t) は時間周波数成分パワー |Y (x, t)|2 を要素に対応する成分に分解するためのいわば連続値マスクと見なせる。この不等式は，サブクラスタモデル Wk,n,y (x, t) と，連続値マス
クにより分解された成分との I ダイバージェンスの総和をとったものが，全体同士の I ダイバージェンスより決して小さくな
ることはないことを示唆している。
この事実は，次のような反復推定アルゴリズムの導出へとつながる。まず，I + (Θ, m) を任意の Θ のもとで mk,n,y (x, t)
に関して最小化すると，不等式より I + (Θ, m) は I(Θ) より小さくなることなく等号成立状態に至る。導出は省略するが，I + (Θ, m) を最小化する mk,n,y (x, t) は，変分法に基づき，Wk,n,y (x, t)
mk,n,y (x, t) = X
Wk,n,y (x, t)

(18)

k,n,y

と解析的に得られる。この結果を式 (17) 右辺に代入すると，たしかに I + (Θ, m) = I(Θ) となることが分かる。続いて，この状態から I + (Θ, m) を最小化する，あるいは減少させる Θ を求めることができるならば，不等式 (17) より，減少させられた I + (Θ, m) よりさらに I(Θ) は小さくなっているはずである(図 6 参照)。従って，次のような 2 ステップの反復計算
Step 0 初期値 Θ0 を選び，` = 1 とおく。
Step 1 m(`) = argmin I + (Θ(`−1) , m)m
Step 2 I + (Θ, m(`) ) 5 I + (Θ(`−1) , m(`) ) なる Θ をΘ(`) ，` ← ` + 1 として Step 1 へ戻る。
を行うことで I(Θ) を単調減少させていくことができる。I(Θ)と定義される。以後これを最小化する Θ を求めるためのアルゴリズムについて議論する。
まず，負の対数関数が凸関数であることに注目すると，X

mk,n,y (x, t) = 1, ∀k, n, y : mk,n,y (x, t) ∈ (0, 1) (16)k,n,yを満たすような任意の重み関数 mk,n (x, t) を導入しては下に有界なので，上記の反復計算は収束性が保証される。
この反復計算は，I(Θ) の最小化が解析的に行えない代わりに I + (Θ, m) を最小化あるいは減少させられる Θ の更新式が解析的に得られる場合に特に効果的である。ところで，この反復計算は，実は，EM(Expectation-Maximization) アルゴリズ
ム [13] と形式上は等価である。しかしながら，文献 [13] の導出に直接従えば，|Y (x, t)|2 と |W (x, t)|2 がいずれも確率密度関数でなければならず，パワースペクトルに対して Bayes 則 (条—4—

件つき確率の公式や周辺化操作) が適用できるかどうかは自明ではない。その意味で，上記のように Bayes 則を適用せずとも

EM アルゴリズムと形式的に等価な反復アルゴリズムが導ける
ことは，本問題設定上においては重要な意味をもつ。
Y −1
Step 2 における {Ωk,i }Ii=1 ,{vk,n }N
n=1 ,{uk,y }y=0 ,wk ,τk ,φk の更新式は，導出は省略するが，
(`)

wk =

XZZ

(`)

τk

=

(`)

wk

(`)

(`)

wk

(`)
φk =

I(Θ

(`)

wk

(`)

(`−1)

Φk,n,y (x, t)(t − yφk

)dxdt

(`)

Φk,n,y (x, t)dxdt(`)

Φk,n,y (x, t)

n

³¡

1

(`)

(`)

¢1

(`) 2

αk 2 + 4βk wk

(`)
2wk

(`)

´

− αk

図 7 混合音声信号のスペクトログラムと推定されたピッチ軌跡関数のように解析的に得られる。ただし，¯¯

推定最大パラメータを同等な反復アルゴリズムで求めることが

¯
¯

2
(`)
(i)
Φk,n,y (x, t) = mk,n,y (x, t) Y (x, t)

できる。これにより，他のさまざまな先験的制約をスムーズにである。なお，上記の Ωk,0 は，ピッチ周波数軌跡が多項式 (式(13)) の場合の 0 次項の係数である。1 次以上の項の係数の更新式も同様に解析的に得られるが，ここでは紙面の都合上省略する。一方で，ピッチ周波数軌跡が 3 次 Spline 関数 (式 (14))の場合の各係数の更新式は，

XZZ ³
n,y

¡

ek,n,i t; Ωk,i
x−Ω
(`)

(`)

¢´ ∂Ωk (t)X Z Z µ ∂Ωk (t) ¶2
n,y

∂Ωk,i

∂Ωk,i

k,i

k,i

は t と定数行列 M にのみ依る項である。ところで，σ は定 Q
フィルタバンクの解析条件で決まる定数であったが，これも音源 k ごとの変数と見なして推定することもできる。その場合，
σk の更新式は同様にして

XZZ

(`)

wk

聴覚においてひとまとりと知覚される単位が物理的な意味での一連の音響信号と対応するとは限らないが，実際のアプリかを評価実験により確認することは，方法論自体の工学的な実と結びついているかを示す 1 つの目安にもなると考えられる。

e13. 性能の定量的評価

用可能性だけでなく，Bregman の分凝要件がいかに物理的事実(`)
Φk,n,y (x, t)dxdt

(`)
(`)
(`)
(`−1)
(`−1)
, Ωk,i+1 ,
(`−1)
(`)
(`)
(`)
. . . , Ωk,I ) と す る 。ま た ，Ωk,n,i (t; Ωk,i ) = Ωk (t; Ωk,i ) −∂Ωk (t) (`)
k (t)Ωk,i + log n であり，この項は Ωk,i に依らない。 ∂Ω
∂Ω
∂Ω

Ã

導入できる点もここで強調しておく (文献 [15], [16] 参照)。

ケーションを想定して提案する方法論がどの程度の性能を示す(`)

Φk,n,y (x, t)dxdt

で与えられる。ただし，Ωk,i = (Ωk,1 , . . . , Ωk,i−1 , Ωk,i

(`)
σk =

)

図 6 パラメータ反復推定アルゴリズムの概念図。

n,y

(`)

(l)

(`)

Φk,n,y (x, t)(x − log n)dxdt

X Z Z (`)
(`)
(`)


Φk,n,y (x, t)y(t − τk )dxdt αk ,n,y
X Z Z (`)
(`)
(`)


Φk,n,y (x, t)(t − τk )2 dxdt
 βk ,Ωk,i =

step 2

y

ZZ
1 X

(`)

uk,y =

I + (Θ (l-1), m ( l ))

n,y

ZZ
1 X

(`)

vk,n =XZZ

1)

I + (Θ ( l ,) m( l ))

n,y

wk

I(Θ

(`)

ZZ
1 X

(`)

step 1(l-1)

Φk,n,y (x, t)dxdt

n,yΩk,0 =

I + (Θ (l-1), m (l-1))

!12
(`)
(`)
Φk,n,y (x, t)(x − Ωk (t) − log n)2 dxdt

n,y

と求まる。
また，本稿では紙面の都合上割愛せざるを得なかったが，式(15) の最小化は，回帰分析の問題と捉えると，|W (x, t)|2 をパラメータとする |Y (x, t)|2 の連続 Poisson 分布の結合分布を尤度とした最尤推定と等価であることが分かっており，それにより，適当なパラメータの事前分布を定めることにより事後確率

3. 1 混合音声のピッチ周波数推定精度
まず，提案法の動作確認のため，単一話者音声のピッチ周波数の推定精度を評価した。女性と男性話者による計 100 個の英文の読み上げ音声の Laryngograph 信号 [14] を実験データとした。ここでは，HTC のピッチ周波数軌跡は 3 次スプライン関数を仮定することとした。分析条件やピッチ周波数の正解データの作成方法などに関する詳細は文献 [16] に示されている。推定
精度は，サンプリング周期ごとに推定ピッチ周波数と正解ピッ
チ周波数を照らし合せ，正解からの誤差が 20%以内だったもの
を正解として算出した。表 1 に，現在最先端のピッチ周波数推定手法である YIN [17] やそれ以外の手法との誤り率の比較を示す。YIN 以外の略称が指す手法名については文献 [16] に詳し
く示されている。単ピッチ推定に特化したこれらの従来手法とHTC の性能が遜色がないことは特筆すべき結果と言える。
続いて，同じ分析条件で，2 話者 (男性と女性，男性同士，女性同士) の音声を平均パワーが等しくなるように混合した 150
個の混合音声信号を対象とし，混合音声信号のピッチ推定精度を評価した。CASA の最先端の手法の 1 つである Wu らの手—5—Energy density

na
cf
ad
dit
iv
TE e
MP
O
YI
N
HT
C

s
ac
f

sh

cc

Method

pd
a
fxa
c
fxc
ep
ac

表 1 単ピッチ推定精度の種々の従来手法との比較

Gross
19.016.8 15.8 9.2 6.8 12.8 1.9 1.7 3.6 3.2 1.4 3.5
error (%)

表 2 混合音声信号のピッチ推定精度の WWB 法との比較

81.5
69.0
90.8
83.5

time (s)Energy density

86.8
87.9
95.6
90.2

nt)

81.8
83.4
95.8
87.0

(ce

93.3
96.1
98.9
96.1

cy

男性と女性
男性同士
女性同士
Total

10%
WWB

en

HTC

qu

20%
WWB

fre
g-

HTC

Lo

Gross error thresholdmethods

法 (以後，WWB 法) [11] との性能比較を表 2 に示す。ここでは，正解と見なす許容誤差を ±10% と ±20% とした場合の両者の精度を評価し，HTC がいずれにおいても WWB 法を上回る成績を得た。提案法の動作を直感的に確認するため，混合音
声信号のスペクトログラムと推定されたピッチ軌跡関数を重ね

Lo
c
en
nt)

ce

y(

3. 2 音楽音響信号の音高推定精度
次に，音楽音響信号を対象とし，音名の推定精度を評価する
ため，RWC 研究用音楽データベース [18] よりクラシックとジャズの 8 楽曲を実験データに選び，PreFEst [19] と呼ぶ最先端の音楽の音高推定手法と性能比較を行った。ここでの実装では，
HTC のピッチ軌跡関数は 0 次多項式とした。なお，その他の分析条件や正解データの作成方法などに関する詳細は文献 [15]に示されている。実験の結果，PreFEst は平均音名推定精度が62.4%であったのに対し，HTC は 70.4%という成績を得，音楽の音高推定においても HTC が効果的であることを確認した。
時変スペクトルモデルがどの程度正確に推定されているかを視覚的に確認するため，対象とする信号のスペクトログラムと並べて最適推定された時変スペクトルモデルを図 8 に示す。

qu

fre

g-

て表示したものを図 7 に示す。

time (s)図 8 対象のスペクトログラム (上) と推定された時変スペクトルモデル (下)

[8]
[9]

[10]

[11]

[12]

4. ま と め
本稿では，Bregman の分凝要件に立脚し，個々の音源の時間[13]

周波数全域に渡ったスペクトル構造を一挙に推定できる新しい方法論を提案し，これを調波時間構造化クラスタリング (HTC)と名づけた。混合音声信号および音楽音響信号のピッチ周波数[14]推定精度の評価実験を通し，それぞれの分野における最先端の従来法の精度を上回る結果を得た。
文
[1]
[2]

[3]
[4]
[5]
[6]

[7]

献

A. S. Bregman, Auditory Scene Analysis, MIT Press, Cambridge, 1990.
M. Cooke, G. J. Brown, M. Crawford and P. Green, “Computational Auditory Scene Analysis – Listening to Several Thingsat Once,” Endevour New Series, Vol. 17, No. 4, pp. 186–190,
1993.
G. J. Brown and M. Cooke, “Computational Auditory Scene
Analysis,” Comp. Speech & Lang., No. 8, pp. 297–336, 1994.
D. P. W. Ellis, “A Computer Implementation of Psychoacoustic
Grouping Rules,” In Proc. IEEE ICPR’94, pp. 108–112, 1994.
A. Fishbach, “Primary Segmentation of Auditory Scene,” In
Proc. IEEE ICPR’94, pp. 113–117, 1994.
T. Nakatani, M. Goto and H. G. Okuno, “Localization by Harmonic Structure and Its Application to Harmonic Sound Segregation,” In Proc. IEEE ICASSP’96, pp. 653–656, 1996.
西, 安部, 安藤, “聴覚情景分析のための多重ピッチ追跡と調波分離アルゴ
リズム,” 計測自動制御学会, Vol. 34, No. 6, pp. 483–490, 1998.

鵜木, 赤木, “聴覚の情景分析に基づいた雑音下の調波複合音の一抽出法,”
電子情報通信学会論文誌, Vol. J82-A, No. 10, pp. 1497–1507, 1999.

[15]

[16]

[17]

[18]

[19]

安部, 安藤, “共有 FM-AM の時間周波数統合に基づく聴覚情景分析(I)—Lagrange 微分特徴量とその周波数統合—,” 電子情報通信学会論文誌, Vol. J83-D-II, No. 2, pp. 458–467, 2000.
安部, 安藤, “共有 FM-AM の時間周波数統合に基づく聴覚情景分析(II)—最適な時間軸統合とストリーム音の再合成—,” 電子情報通信学会
論文誌, Vol. J83-D-II, No. 2, pp. 468–477, 2000.
M. Wu, D. L. Wang and G. J. Brown, “A Multipitch Tracking
Algorithm for Noisy Speech,” IEEE Trans., Speech and Audio
Process., Vol. 11, pp. 229–241, 2003.
Ö. Yilmaz and S. Rickard, “Blind Separation of Speech Mixtures via Time-Frequency Masking,” IEEE Trans., Signal Process., Vol. 52, No. 7, pp. 1830–1847, 2004.
A. P. Dempster, N. M. Laird and D. B. Rubin, “Maximum
Likelihood from Incomplete Data via the EM Algorithm,” J.
of Royal Statistical Society Series B, Vol. 39, pp. 1–38, 1977.
P. C. Bagshaw, S. M. Hiller and M. A. Jack, “Enhanced Pitch
Tracking and the Processing of F0 Contours for Computer andIntonation Teaching,” In Proc. Eurospeech’93, pp. 1003–1006,
1993.
H. Kameoka, T. Nishimoto and S. Sagayama, “A Multipitch
Analyzer Based on Harmonic Temporal Structured Clustering,”
IEEE Trans., Audio, Speech, Language Process., in Press, 2006.
J. Le Roux, H. Kameoka, N. Ono, A. de Cheveigné and
S. Sagayama, “Single and Multiple Pitch Contour Estimation
through Parametric Spectrogram Modeling of Speech in Noisy
Environments,” IEEE Trans., Audio, Speech, Language Process., submitted in 2006.
A. de Cheveigné and H. Kawahara, “YIN, a Fundamental Frequency Estimator for Speech and Music,” J. Acoust. Soc. Am.,
111(4), pp. 1917–1930, 2002.
後藤, 橋口, 西村, 岡, “RWC 研究用音楽データベース: 研究目的で利用可能な著作権処理済み楽曲・楽器音データベース,” 情報処理学会論文誌,
Vol.45, No.3, pp.728–738, 2004.
M. Goto, “A Real-Time Music-Scene-Description System:
Predominant-F0 Estimation for Detecting Melody and Bass
Lines in Real-World Audio Signals,” ISCA Journal, Vol. 43,
No. 4, pp. 311–329, 2004.

—6—