3-1-6
高次相関を用いた非線形 MUSIC による高分解能方位推定∗☆杉本侑哉, 宮部滋樹, 山田武志, 牧野昭二 (筑波大)

1

はじめに

音源信号の存在する空間を分析するマイクロホン
アレー信号処理において，音源の方位推定は，事前情報に音源の方位情報を要する音源分離や雑音抑圧など，幅広い分野で求められ，様々な研究が行われている．空間共分散行列の高位主成分ベクトルとステアリ
ングベクトルの相似性を利用した MUSIC(MUltiple
SIgnal Classification)[1] は，高分解能に方位推定を行う手法として知られている．しかし，MUSIC は共分散行列の次元数に拘束される．N 個の音源の方位を推定するためには，M (> N ) 個の素子を用いる必要があり，また N が M に近付くと推定精度が低下する．従って，多数の音源の方位を推定する場合，アレーの大規模化が避けられないという問題がある．このような背景から，統計量として高次キュムラントを用い，次元数を増やした分析へと MUSIC を拡張する
いくつかの手法が提案されてきた．例えば，統計量に4 次キュムラントを用いる 4-MUSIC[2] や，4-MUSICをさらに高次元に拡張する 2q-MUSIC[3] が知られている．これらの手法は，音源数が素子数よりも多い劣決定条件の下での方位推定が可能であり，同時に方位推定の分解能を向上させる．
本稿では，4-MUSIC などとは異なったアプローチ
で劣決定条件に MUSIC を拡張する手法を提案する．
パターン認識の分野では，カーネル関数を用いた高次元ヒルベルト空間への非線形写像により，識別器の線
形識別性を高めるカーネル法が知られている．また，カーネル法を用いて観測信号を高次元化することで，ビームフォーミングの性能を高める手法 [4] が提案されている．カーネル法により処理の次元数を増やす
アプローチは，共分散行列の次元数に拘束を受ける
MUSIC にとって極めて有効であると考えられる．そこで我々は，カーネル法を用いて高次元ヒルベルト空間で MUSIC を行うカーネル MUSIC を提案し，定式化を行う．どのような条件のカーネル関数が MUSICの分析に適しているかを定め，その中で最も単純な二次カーネルに注目する．統計的・幾何学的な視点からの議論を通じて，二次カーネルを用いた二次カー
ネル MUSIC が 4 次モーメントの分析であることと，(M 2 − 1) までの音源を高分解能に方位推定可能であることを示す．さらに，観測信号の時間長が長い場合に計算が効率的な，二次カーネルの与える非線形写像を解析的に表現して用いたエルミート MUSIC を提案し，従来手法に比べてより少ない計算量で方位推定が可能であることを示す．実験により，音源の方位推定性能が従来手法とほぼ同等で，残響に対してやや頑健であることを確認する．

2x(ω, t) = [x1 (ω, t), · · · , xM (ω, t)]T=N∑ai (ω)si (ω, t) + n(ω, t),

(1)

(2)

ai (ω) = [a1,i (ω), · · · , aM,i (ω)] .

(3)

T

ここで，t = 1, · · · , L はある処理区間における時間フレームインデックス，ω は周波数，M は素子数，N
は音源数，si (ω, t) は i 番目の音源信号の複素振幅，
xj (ω, t) は j 番目の素子における観測信号の複素振幅，nj (ω, t) は j 番目の素子における観測雑音の複素振幅，
aj,i (ω) は i 番目の音源から j 番目のマイクロホンへの伝達関数を表す．ai (ω) は，各音源の到来方位を表現するベクトルであることから，本稿ではこれを i 番
目の音源の方向ベクトルと呼ぶ．
式 (1) に示される信号モデルにおいて，各素子では複数の音源から到来する信号と，雑音が線形に混ざり合って観測される．方位推定技術の目標は，このような混合信号から各音源の方向ベクトル a1 (ω), · · · , aN (ω)
を個別に同定することである．従って，本稿での問題
設定を観測信号ベクトル x(ω, t) を用いて，真の音源方位に対応する方向ベクトルと類似のステアリング
ベクトルを推定することと定める．

3

従来手法

3.1 MUSIC
ここでは，MUSIC について述べる．MUSIC は，空間共分散行列の上位主成分ベクトルとステアリング
ベクトルの相似性を利用し，マイクロホンアレーで観測された混合信号から各音源の方位推定を行う手法である．マイクロホンアレーの素子数が M の場合，
(M − 1) までの音源方位を高い分解能で推定できる．
前節でモデル化された観測系において，n(ω, t) の各要素が互いに無相関であると仮定すると，観測信号の共分散行列は以下のように表される．
R(ω) = E[x(ω, t)xH (ω, t)]
=

N
∑

2
Pi (ω)ai (ω)aHi (ω) + σ (ω)IM .

(4)

i=1

ここで，Pi (ω) = E[|si (ω, t)|2 ] は，i 番目の音源の信号パワー，σ 2 (ω) は雑音パワー，Ij は j 次元の単位行列を表す．共分散行列 R(ω) の固有値分解により，次式を得る.

R(ω) = V(ω)E(ω)VH (ω),(5)

V(ω) = [v1 (ω), · · · , vM (ω)], V (ω)V(ω) = IM ,(6)
H

E(ω) = diag[e1 (ω), · · · , eM (ω)] + σ 2 (ω)IM ,e1 (ω) ≥ · · · ≥ eN (ω),
eN +1 (ω), · · · , eM (ω) = 0.

問題設定

本稿を通して信号は周波数領域で表されることとする．観測信号は次のようにモデル化できる．

n(ω, t) = [n1 (ω, t), · · · , nM (ω, t)]T ,(7)

こ こ で ，e1 (ω), · · · , eM (ω) は R(ω) の 固 有 値 ，v1 (ω), · · · , vM (ω) は 各 固 有 値 に 対 応 す る 固 有 ベ
ク ト ル と し ，diag[· · · ] は 引 数 の 要 素 か ら な る 対角行列を表す．音源信号に対応する固有ベクトルi=1∗

High-resolution DOA estimation with non-linear MUSIC based on higher-order corvariance. by Yuya
SUGIMOTO, Shigeki MIYABE, Takeshi YAMADA, and Shoji MAKINO. (University of Tsukuba)

日本音響学会講演論文集- 763 -

2012年3月

v1 (ω), · · · , vN (ω) によって張られる部分空間は，各音源の方向ベクトル a1 (ω), · · · , aN (ω) が張る部分空間と一致する．この部分空間を信号部分空間と呼び，S(ω) と表す．
N
S(ω) , span[ai (ω)]Ni=1 = span[vi (ω)]i=1 .

(8)さらに，span[vi (ω)]Mi=1 において信号部分空間の直交補空間を雑音部分空間と呼ぶと，雑音部分空間の基
底ベクトル vN +1 (ω), · · · , vM (ω) と各音源の方向ベクトルの間には以下の関係が成り立つ．

aH
i (ω)vj (ω) = 0,(i, j) ∈ {i = 1, · · · , N j = N + 1, · · · , M }.

(9)

式 (9) の直交性を利用して，次式に示す MUSIC スコ
ア f (ω; θ) を考える．

f (ω; θ) , ∑M1

H
i=N +1 |b (ω; θ)vi (ω)|

.

(10)

ここで，b(ω; θ) (||b(ω; θ)|| = 1) は方向 θ のステアリ
ングベクトルを表す．この時，i = 1, · · · , N いずれかについて bH (ω; θ)ai (ω) = ||ai (ω)|| に近づく θ に対して，f (ω; θ) は大きな値を取る．MUSIC は，この性質を利用して，f (ω; θ) が大きな値を取る方向を音源方向として推定する手法である．なお，式 (10) によって求まるのはある狭帯域での方位推定結果であるが，所望の周波数帯域における単一の推定結果は，各周波数ビンでの推定結果を周波数平均することで求められる．
MUSIC の注意すべき性質として，式 (4) の第二項が単位行列の定数倍でない時に固有ベクトルの値が変わることから，各素子の観測雑音に相関がある場
合推定性能が低下する．また，雑音部分空間を表現するために少なくとも一次元の部分空間が必要となるため，推定する音源数 N は素子数 M 未満に拘束される．高分解能な方位推定を行うためには雑音部分空間の高い表現能力が求められ，M は N より十分大
きくなくてはならない．

3.2 4-MUSIC
4-MUSIC は，高次統計量として 4 次キュムラント
を用いた音源の方位推定手法である．マイクロホンア
レーの素子数が M の場合，音源数 (M 2 − 1) までの方向推定が可能であり，方位推定の分解能も MUSIC
より高い．4-MUSIC では，各要素が 4 次キュムラン
トよりなる M 2 × M 2 行列 C(ω) を固有値分解する．
4 次キュムラントは式 (11) のように計算される．

り方位推定の分解能が向上し，また推定する音源数が素子数よりも多い劣決定条件でも方位推定が可能である．写像される空間の次元はカーネルの選択に依存し，無限次元への写像も可能であるが，本稿ではパワー推定型のビームフォーマ [4] で用いられるもっとも単純な二次カーネルによる M 2 次元空間への写像に焦点を当てて，幾何学的・統計的な視点から議論する．さらに，二次カーネルを用いたカーネル MUSICと等価で，M 2 < L (L はある処理区間における時間
フレームインデックスの総数) の場合に計算が効率的なエルミート MUSIC アルゴリズムを提案する．

4.1 一般化カーネル MUSIC の定式化
正定値カーネル関数 : CM × CM → C は，カーネル関数が決定する写像 φ : CM → H で与えられるヒ
ルベルト空間 H における内積となることが知られている．
k(x, y) = hφ(x), φ(y)iH .
(12)
ここで，x, y ∈ CM ，h·, ·iH はヒルベルト空間 H
における内積を表す．非線形なカーネル関数により与えられるヒルベルト空間 H は一般に CM よりも高次元であるため，カーネル関数は識別問題において特徴ベクトルの線形識別性を高めるために頻繁に用いられる．MUSIC における共分散行列 R(ω) =E[x(ω, t)xH (ω, t)] の固有値問題をヒルベルト空間
H における相関線形作用素 Γω : H → H, Γω =E[h·, φ(x(ω, t))iH hφ(x(ω, t)), ·iH ] の固有値問題に，
ステアリングベクトル b(ω; θ) を写像 φ(b(ω; θ)) に置き換えることにより， 高次元写像を用いて分解能を高めたカーネル MUSIC が実現される．
一般に写像関数 φ は解析的に記述できないため，相関線形作用素 Γω もまた解析的に記述できない．しかし，Γω の固有値問題

Γω (υi (ω)) = λi (ω)υi (ω),
hυi (ω), υj (ω)iH = δij , λ1 (ω) ≥ λ2 (ω) ≥ · · · ≥ 0.
(13)
(δij は，クロネッカ―のデルタ) における固有値 λi (ω)と固有ベクトル υi (ω) ∈ H を，観測信号 x(ω, t), t =1, · · · , L の L × L グラム行列 K(ω) の固有値分解を用いて以下のように記述する方法がカーネル PCA として提案されている [5]．
λi (ω) = di (ω)/L,1
υi (ω) = √uji (ω)φ(x(ω, j)),di (ω) j=1= U(ω)D(ω)UH (ω),

E[Xi Xj Xk Xl ] − E[Xi Xj ]E[Xk Xl ]

H

(11)

方位推定は固有ベクトルの分析によって MUSIC と同様の手順で行われる．

4

提案手法

パターン認識の分野では，カーネル関数を用いた
高次元のヒルベルト空間への写像により，識別器の線
形識別性を高めるカーネル法が知られている．本節ではカーネル法を MUSIC に応用したカーネル MUSICを提案する．雑音部分空間の次元数を増やすことによ日本音響学会講演論文集(15)

K(ω) = [k(x(ω, i), x(ω, j))]ij

Cum[Xi , Xj , Xk , Xl ] =− E[Xi Xk ]E[Xj Xl ] − E[Xi Xl ]E[Xj Xk ].

(14)

L
∑

(16)

U(ω) = [uji ]ji , U (ω)U(ω) = IL ,(17)

D(ω) = diag [d1 (ω), · · · , dL (ω)] ,d1 (ω) ≥ · · · ≥ dL (ω) ≥ 0.

(18)

ここで，[x]ij は第 i 行 j 列要素が x である行列を表す．固有値 λi (ω) の比零の値の数を r(ω) とすると，
r(ω)
span[υi (ω)]i=N +1 が雑音部分空間に相当するため，ヒルベルト空間 H における MUSIC スコア gH (ω) は以下のように求められる．

gH (ω) , ∑r(ω)- 764 -

1i=N +1 |hφ(b(ω; θ)), υi (ω)iH |2012年3月

=

1
∑r(ω)

∑L

i=N +1√uji (ω)j=1

di (ω).
k (b(ω; θ), x(ω, j))(19)

4.2 カーネル関数の適切な選択
カーネル MUSIC はヒルベルト空間 H における相関の分析であり，写像 φ の選択によってその性質が大幅に変化する．そのためカーネル関数を適切に選ぶ必要がある．まず信号の写像は，1) ノルムの大小関係が維持され，kxk > kyk のときに kφ(x)kH > kφ(y)kH（k·kH はヒルベルト空間におけるノルム）となり，2)写像が原点を移動させず，φ(0) が H における原点に一致し，かつ 3) 信号間の直交性が維持されて xH y = 0のときに hφ(x), φ(y)iH = 0 となるという 3 つの条件を満たす必要がある．また，ヒルベルト空間への写像のためにはカーネル関数が正定値カーネルである必要がある．これらの条件を満たすカーネル関数の例
として，次式が考えられる．
k(x, y) =∞∑pi |x y| x y +
H

2i H

i=0

∞
∑

qi |x y| ,
H

2i

i=1

pi , qi ≥ 0 for ∀i.

(20)

本稿では，式 (20) を満たすカーネル関数の中で最も単純で，またパワー推定型の非線形ビームフォーマに用いられる以下の二次カーネル k2 (·, ·) に焦点を当てる．
2
k2 (x, y) , xH y .
(21)
二次カーネル k2 (·, ·) によって再生されるヒルベルト2
空間 H2 は RM であり，与えられる写像 φ2 : CM →2
RM は以下のようにユークリッドベクトル関数として解析的に表すことができる．
TTT
φ2 (x) , [φT
abs (x), φre (x), φim (x)] ,(22)

φabs (x) , [∀|xi | |1 ≤ i ≤ M ] ,(23)√
φre (x) , 2[∀Re[xi x∗j ]|2 ≤ i ≤ M, 1 ≤ j ≤ i − 1]T ,(24)√∗
φim (x) , 2[∀Im[xi xj ]|2 ≤ i ≤ M, 1 ≤ j ≤ i − 1]T .
(25)
2

T

本稿では写像 φ2 をエルミート展開写像と呼ぶこととする．
以下では，二次カーネルを用いたカーネル MUSICの幾何学的・統計的性質について議論する．エルミー2ト展開写像 φ2 : CM → RM の定義より，特徴ベク
トルから高次ユークリッド空間への写像によって次元数が M から M 2 へと増加していることが分かる．このとき，エルミート展開写像 φ2 の与える相関線形作
用素 Γω を用いた MUSIC は，次元数の増加による雑音部分空間の表現能力の向上によって (M 2 − 1) までの音源の方位推定が可能である．また，式 (22)–(25)から分かるように，エルミート展開写像は観測信号のx(ω, t) の二次関数による写像であり，エルミート展
開写像の与える相関線形作用素を用いた二次カーネル MUSIC は 4 次のモーメントの分析となっていることが分かる．同様に，式 (20) においてより高次の係
数 pi , qi が非零となるカーネル関数を用いたカーネル
MUSIC は，より高い偶数次モーメントを用いた分析となることが分かる．

4.3 エルミート MUSIC
4.1 節で議論したカーネル MUSIC に二次カーネル
k2 (·, ·) を用いる場合は L 次元の固有値問題を解く必要があるが，L > M 2 となる条件では，エルミート
展開写像の与える M 2 × M 2 共分散行列 R2 (ω):][(27)
R2 (ω) = E φ2 (x(ω, t))φT
2 (x(ω, t))の固有値問題を解く方が計算が効率的である．本節では，R2 (ω) を用いるエルミート MUSIC アルゴリズムを提案し，4-MUSIC より計算が効率的であることを示す．共分散行列 R2 (ω) の計算のためには，式 (22)–(25) の写像 φ2 (x(ω, t)), t = 1, · · · , L を求めるための2M 2 L 回の実数乗算と，φ2 (x(ω, t)) を R2 (ω) に変換するための M 2 (M 2 + 1)L/2 回の実数乗算のみが必要
であり，4-MUSIC におけるキュムラント行列 C(ω)
の計算よりも乗算の回数が少ない．また，4-MUSICにおける M 2 × M 2 次元キュムラント行列 C(ω) が複素行列であるのに対し，R2 (ω) は同じ M 2 × M 2 行
列であるが実数行列となるため，固有値問題を解くための計算が大きく単純化される．
エルミート展開写像を直接用いたエルミート MUSIC の MUSIC スコア g2 (ω; θ) は以下のように求められる．

g2 (ω; θ) , ∑M 2

T

ここで x = [x1 , · · · , xM ] , xi ∈ C である. この等価性は以下の内積の一致により証明される．

k2 (x, y) = xH y
=

M
∑

2

2

2

|xi | |yi | + 2

i=1

+2

M ∑
i−1
∑

[
]
[
]
Re xi x∗j Re yi yj∗

i=2 j=1
M ∑
i−1
∑

[
]
[
]
Im xi x∗j Im yi yj∗

1T
i=N +1 φ2 (b(ω; θ))wi (ω)
T

5

,

(28)

R2 (ω) = W(ω)Λ(ω)W (ω),W(ω) = [w1 (ω), · · · , wM 2 (ω)] ,(29)

WT (ω)W(ω) = IM 2 ,Λ(ω) = diag [λ1 (ω), · · · , λM 2 (ω)] ,λ1 (ω) ≥ · · · ≥ λM 2 (ω) ≥ 0.

(30)
(31)

実験

i=2 j=1= φT
2 (x)φ2 (y).

(26)

ここで式 (22)–(25) に注目する．φabs (x) は，エルミー
ト行列 xxH の対角成分よりなるべクトルとなっている．同様に，φre (x)，φim (x) は xxH の対角成分を除いた下三角要素の実部と虚部よりなるベクトルとなっている．入力ベクトル x より得られるエルミート行
列 xxH の成分に対応した写像となっていることから，

日本音響学会講演論文集

本章では，提案するエルミート MUSIC の有効性を検証するために，計算機によって行ったシミュレー
ション実験について述べる．比較のための従来手法として，3 章に示した MUSIC，4-MUSIC に加えて，
遅延和アレー，GCC-Eckart[6] を用いた．

5.1 実験条件
提案するエルミート MUSIC の方位推定性能の評
価を，音声混合のシミュレーションを用いて行った．

- 765 -

2012年3月

Experimental location

Estimation result

5

1.1

30°

1

4

90°
3

300°

2

240°
150°
1

Estimation score

0.9
0.8

Delay and sum
GCC-Eckart
MUSIC
4-MUSIC
Proposed method

0.7
0.6
0.5
0.4
0.3

1

2

3

4

5

90

180

270

360

Direction [degree]

Fig. 1

Experimental results of estimation performance.

Table 1
アレー形状

鏡像法 [7] によりインパルス応答を与え，異なる 3 つの残響時間で音声混合した．観測信号には，拡散性雑音を模擬した 3 つの異なるピンクノイズを重畳した[8]．Table 1 に各実験条件について示す．性能評価には，{0◦ , 30◦ , 60◦ , · · · , 330◦ } の 30◦ 間隔 12 方向における全組み合わせで求めた，真の音源方位と各手法の推定方位の間の最小二乗誤差平方根平均（RMSE）を用いた．また，音源数分のピーク探索結果が得られず，正しく方位推定が行えなかった場合は，全ての推定方位のエラーを平均したペナルティを与えた．
なお，M ≤ N の場合における MUSIC の方位推定
は，各帯域で最小固有値に対応する一次元部分空間
を雑音部分空間とみなし，広帯域のピークを検出することにより行った．

Delay and sum

4-MUSIC

Proposed method

101
10

1010
T60 = 0 s T60 = 0.12 s T60 = 0.3 s T60 = 0 s T60 = 0.12 s T60 = 0.3 s T60 = 0 s T60 = 0.12 s T60 = 0.3 s

SNR = 0dB

Fig. 2

SNR = 10dB

SNR = 20dB

Experimental results of RMSE.

が効率的なエルミート MUSIC を提案し，二次カー
ネルの幾何学的・統計的な性質について述べた．従来手法との性能比較から，ほぼ同等の推定性能であり，高残響下においては提案手法の推定性能が優れていることを示した．計算が効率的であることからも，提案手法は十分有効な手法であることを確認した．
今後，より高次に拡張したカーネル関数での性能評価に取り組む予定である．
謝辞 この研究は，科研費たものである．

23760325 の支援を受け

参考文献

おわりに

本稿では，劣決定条件で高分解能に複数音源の方位を推定する手法として，カーネル法と MUSIC を組
み合わせたカーネル MUSIC を提案した．また，二次カーネル MUSIC と等価で，M 2 < L の場合に計算

日本音響学会講演論文集

MUSIC

1002
10

0.1-1
10

5.2 結果と考察
Fig. 1 にある方位の組み合わせのもとで，残響時間が 0.12 s の場合の方位推定の結果を示す．左側が，音源とマイクロホンアレーの配置図，右側が推定結果のグラフとなっている．Fig. 2 に RMSE の全試行における平均を示す．結果から，提案手法と 4-MUSICがほぼ同等の優れた方位推定性能を示し，また高残響化において提案手法が優れていることが分かる．計算が効率的であることからも，提案手法は十分有効な手法であることが確認できる．

6

GCC-Eckart

1000
103

RMSE [degree]

Experimental conditions
円形アレー
アレーから距離 1.5 m の
音源
点音源
素子数
4
音源数
5
アレー半径
0.2 m
部屋サイズ
5 m×5 m
残響時間 (T60 )0, 0.12, 0.3 [s]SNR0, 10, 20 [dB]
サンプリング周波数
16 kHz
512 for T60 of0s
フレーム長1024 for T60 of 0.12 s
2048 for T60 of 0.3 s
シフト長
1/4 フレーム
窓関数
ハミング窓- 766 -

[1] R. O. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Trans. Antennas and Propagation, vol.AP–34, no.3, 1986.
[2] B. Porat and B. Friedlander, “Direction ﬁnding algorithms based on higher order statistics,” IEEE Trans.
Signal Processing, vol.39, no.9, pp.2016–2024, 1991.
[3] P. Chevalier, A. Ferreol and L. Albera, “High resolution direction ﬁnding from higher order statistics: The
2q-MUSIC algorithm,” IEEE Trans. Signal Processing,vol.54, no.8, pp.2986–2997, 2006.
[4] S. Miyabe, et al., “Kernel-based nonlinear independent
component analysis for underdetermined blind source
separation,” Proc. ICASSP, pp.1641–1644, 2008.
[5] B. Scholkopf, A. Smola and K. R. Muller, “Support vector methods in learning and feature extraction,” Australian J. of Intelligent Information Processing Systems,
vol.10, no.5, pp.1299–1319, 1998.
[6] C. H. Knapp and G. C. Carter, “The generalized correlation method for estimation of time delay,” IEEE Trans.
Acoustics, Speech and Signal Processing, vol.24, pp.320–
327, 1976.
[7] J. A. Allen and D. A. Berkley, “Image method foreﬃciently simulating small-room acoustics,” J. of the
Acoustical Society of America, vol.65, no.4, pp.943–950,
1979.
[8] I. A. McCowan and H. Bourlard, “Microphone arraypost-ﬁlter based on noise ﬁeld coherence,” IEEE Trans.
Speech and Audio Processing, Vol.11, No.6, pp.709–716,
2003.

2012年3月