1-1-11
複数次数の高次モーメントを用いた非線形 MUSIC ∗☆杉本侑哉, 宮部滋樹, 山田武志, 牧野昭二 (筑波大),Fred Juang (ジョージア工科大)

1

はじめにマイクロホンアレーを用いた信号処理において音源の方向推定（DOA）は，話者追跡や音声強調などの様々な応用技術の前処理として必要不可欠である．高分解能に方向推定を行う代表的な手法として，観測信号の部分空間分析に基づく MUSIC（MUltiple SIgnalClassification）法 [1] が知られている．しかし MUSICは，観測信号の伝達関数ベクトルと直交する雑音部分空間を推定する必要があり，共分散行列の次元数に拘束される．N 個の音源方向を推定するためにはM (> N ) 個のマイクロホンが必要であり，N が Mに近付くほど推定性能は低下する．したがって，多数の音源を推定する場合，アレーの大規模が避けられないといった問題がある．
このような背景から，音源数がマイクロホン数よりも多い劣決定条件においても高分解能に方向推定を行うために，いくつかの MUSIC 拡張手法が提案されている．例えば，90 年代初めの 4 次キュムラント分析に基づく手法 [2] や，近年提案された上手法を任意の偶数次キュムラント分析へと一般化した 2q-MUSIC[3]などが挙げられる．これらの手法は，高次クロスキュムラントを用いてバーチャルに観測信号の次元数を増大させることで，MUSIC よりも高い解像度と劣決定条件での方向推定を達成している．また，異なったアプローチでの MUSIC 拡張手法として，我々は写像
MUSIC[4] を提案している．写像 MUSIC は劣決定条
件で音声強調を行う非線形ビームフォーマ [5, 6, 7] のアイディアを MUSIC に応用した手法であり，任意の偶数次モーメント分析に相当する．推定性能は従来のキュムラント分析に基づく手法とほぼ同等であるが，観測の次元数を増大させるために高次元写像を用いることで大きく計算量が削減されている．
しかし従来の MUSIC 拡張手法には，[3] に示されているように，推定の分解能を高めるためには統計量の次数を上げる必要があるが，同時に統計的バイアスにより安定性が損われるといったトレードオフがあった．本稿ではこれを解決し，さらに分解能を高めた方向推定を目的とした，複数次数の高次モーメント分析に基づく方向推定手法を提案する．提案手法
は，単一の偶数次モーメント分析のための写像を，直和により簡便に拡張した写像を用いて複数次数モーメントの分析を行う．単一の高い次数のモーメント分∗

析において劣化要因となる統計的バイアスの影響を，低次モーメントとの同時分析によって緩和することで，従来の MUSIC 拡張手法に比べて極めて高い分解能を達成する．

2

問題設定

本稿を通して，信号は複素振幅として表現する．観測信号を次のようにモデル化する．

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

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

(3)

i=1

T

ただし，t = 1, · · · , L はある処理区間における時間フレームインデックス，ω は周波数，M は素子数，Nは音源数，si (ω, t) は i 番目の音源信号の複素振幅，xj (ω, t) は j 番目の素子における観測信号の複素振幅，nj (ω, t) は j 番目の素子における観測雑音の複素振幅，aj,i (ω) は i 番目の音源から j 番目のマイクロホンへの伝達関数を表す．(1) に示す信号モデルにおいて，各マイクでは複数の音源から到来する信号と，雑音が線形に混ざり合って観測される．本稿における問題設定を，このような混合信号の分析により各音源 s1 (ω, t), · · · , sN (ω, t) の到来方向を推定することとし，次式に示す任意の方向 θ のステアリングベクトル b(ω; θ) (||b(ω; θ)||2 = 1) から，伝達関数ベクトル a1 (ω), · · · , aN (ω) との誤差がより小さいステアリ
ングベクトル b(ω; θ1 ), · · · , b(ω; θN ) を探索することにより行う．

1b(ω; θ) = √ [exp(−jωτ1 ), · · · , exp(−jωτM )] (4)M
ただし，τi (i = 1, · · · , M ) は方向 θ から到来する信号が，アレーの基準点と各素子の間に作る時間遅延量を表す．

3

写像 MUSIC

3.1

写像 MUSIC による方向推定のアルゴリズム

写像 MUSIC[4] は観測信号 x(ω, t) を非線形関数
0φ : CM → CM (M 0 ≥ M ) により写像し，M 0 次元

Non-linear MUSIC based on the analysis of multiple higher order moments. by Yuya, SUGIMOTO,Shigeki, MIYABE, Takeshi, YAMADA, Shoji, MAKINO (University of Tsukuba), Fred, Jduang (Georgia
Institute of Technology)

日本音響学会講演論文集- 601 -2013年9月

ユークリッド空間においてベクトル φ(x(ω, t)) の共分散行列を用いて MUSIC と同様の分析を行う方向推定手法である．本章では，写像 MUSIC による推定のアルゴリズムについて述べる．まず，写像 MUSICが正しく方向推定を行うためには，写像の前後で各素子間の相関情報が保持される必要がある．この条
件は式 (5)–(7) により表される．

1. ノルムの大小関係を維持する．
||x||2 ≥ ||y||2 → ||φ(x)||2 ≥ ||φ(y)||2(5)2. 原点を移動させない．
x = 0 → φ(x) = 0(7)(8)

ただし，E[·] は期待値演算，[·]H は複素共役転置を表す．共分散行列 R(ω) の固有値分解より，次式を得る．

R(ω) = V(ω)E(ω)VH (ω)(9)V(ω) = [v1 (ω), · · · , vM 0 (ω)]
VH (ω)V(ω) = IM 0(10)E(ω) = diag[e1 (ω), · · · , eM 0 (ω)]
e1 (ω) ≥ · · · ≥ eM 0 (ω)
0

M = dim[φ(x(ω, t))](11)

span[φ(ai (ω))|i = 1, · · · , N ](13)(ただし，span[·] は引数ベクトルの張る空間を表す)
写像空間における雑音部分空間 N (ω) は次式で定義される．
N (ω) , span[vk (ω)|k = N 0 + 1, · · · , M 0 ](14)

このとき N 0 は，写像関数 φ と音源信号数 N に依存して，N から増加した数となるが，伝達関数ベクトルの写像 φ(ai (ω)) と雑音部分空間 N (ω) の間の直交

日本音響学会講演論文集

k = N 0 + 1, · · · , M 0 )(15)

ステアリングベクトルが音源の伝達関数ベクトルを模擬する仮定 (b(ω; θi ) ∼= ai (ω) (i = 1, · · · , N )) と (15)より，雑音部分空間へのステアリングベクトルの写像
φ(b(ω; θ)) の直交射影を評価することで音源方向探索が可能である．以上より，ステアリングベクトルの方向 θ が音源方向 θi と一致した時にピークを形成する，MUSIC スコア f (ω; θ) が次式により定義できる．

1f (ω; θ) , ∑M 0H (b(ω; θ))v (ω)|2k(16)なお，(16) に定義した f (ω; θ) は狭帯域での MUSICスコアであり，広帯域の推定結果を得るためにはこれらの結果を周波数方向に平均する必要がある．周波数平均の方法としてはいくつか考えられるが [8]，本稿では以下の幾何平均を採用する．
[] K1∏
f (θ) ,f (ω; θ)(17)ω

ただし，K は平均演算に用いる周波数ビンの数を表す．推定音源方向は f (θ) のピーク探索により与えられる.

3.2

偶数次モーメント分析のための写像 φd[4] において，我々は任意の 2d 次クロスモーメントd分析のための写像 φd : CM → CM を提案している．
]T[ dd∏∏l⊗l⊗
xc1l (ω, t), · · · ,xc d (ω, t)
φd (x(ω, t)) ,M ll=1(12)た だ し ，e1 (ω), · · · , eM 0 (ω)
は 固 有 値 ，v1 (ω), · · · , vM 0 (ω) は 各 固 有 値 に 対 応 す る 固 有ベクトル，IM 0 は M 0 次元の単位行列，diag[·] は引数を対角成分とする対角行列，dim[·] は次元数を表す．
このとき，伝達関数 a に係る方向性の音源信号の写像と対応する固有ベクトルの本数を N 0 で表すと，' span[vj (ω)|j = 1, · · · , N 0 ](i = 1, · · · , N

k=N 0 +1 |φ

ただし，[·]H は複素共役転置を表す．条件 (5)–(7) を満たす φ により観測信号 x(ω, t) を写像し，ベクトルφ(x(ω, t)) から共分散行列 R(ω) を得る．

R(ω) = E[φ(x(ω, t))φH (x(ω, t))]|φH (ai (ω))vk (ω)|2  1(6)3. 信号間の直交性を維持する．
xH y = 0 → φH (x)φ(y) = 0性は条件 (5)–(7) によって保証されるため，次の条件を満たす．

{xl⊗ ,l=1(18)xx∗(if l is odd)
(if l is even)(19)

ただし，[·]∗ は複素共役を表す．ここで，ckl は φd の計算に用いる x の要素番号を指定するインデックスであり，{1, · · · , M } から d 要素選択する全 M d 通りの重複順列を，任意の順番で行ベクトルとして並べた M d × d インデックス行列 C の (k, l) 成分として与えられる．

11 ··· 1
1 ··· 2  1
C , [ckl ]kl = (20)..
.

M M ··· M
ただし，[·]ij は (i, j) 要素が引数により与えられる行列を表す．写像 φd を用いた写像 MUSIC は，キュムラント分析に基づく従来の MUSIC 拡張手法と比べて，少ない計算量でほぼ同程度の推定性能を得ることができる．

- 602 -

2013年9月

提案手法

φd1 ,··· ,dm , φd1 ⊕ · · · ⊕ φdm

Rotation direction

Window

1070

より高い分解能での方向推定のために，複数次数のクロスモーメントを同時分析するように拡張した
写像 MUSIC を提案する．提案手法は，3 章で述べた
写像 MUSIC に，以下の写像 φd1 ,··· ,dm を適用することで実現される．

170

Door
2095

150

20340

4200

Room size
Reference point of angle 340 x 400 x 270
is displaced by 5∘(21){φd1 , · · · , φdm } = {φd | d = d1 , · · · , dm }400

di+1 > di(i = 1, · · · , m − 1)

(22)
Fig. 1

ただし，⊕ は直和を表す．φd1 ,··· ,dm は，それぞれ次数の異なる 2di 次モーメント分析のための写像 φdi (i =1, · · · , m) の直和として与えられ，φd1 ,··· ,dm の共分散行列 R⊕
d1 ,··· ,dm (ω) は次式のように複数次数のモーメントを含むことから，φd1 ,··· ,dm を用いた写像 MUSICはこれらの複数次数モーメント同時分析に相当する．

R⊕
d1 ,··· ,dm (ω)[]0= E φd1 ,··· ,dm (x0 (ω, t))φHd1 ,··· ,dm (x (ω, t))
Rd1 ,d1 (ω) Rd1 ,d2 (ω) · · · Rd1 ,dm (ω)
 Rd2 ,d1 (ω) Rd2 ,d2 (ω) · · · Rd2 ,dm (ω) 


=
..
..
..
..
.

.
.
.
Rdm ,d1 (ω) Rdm ,d2 (ω) · · · Rdm ,dm (ω)(23)[]0Rdi ,dj (ω) , E φdi (x0 (ω, t))φH
dj (x (ω, t)) = [rpq ]pq(i, j = 1, · · · , m

p = 1, · · · , M

q = 1, · · · , M )∗ ()  djdi∏∏00 x cl⊗  rpq = E x cl⊗
plql
l=1

一方，異なる次数のモーメントを同時分析する場合は観測信号のスケーリングが問題となる．(23)–(25)に示すように，R⊕d1 ,··· ,dm (ω) は様々な次数のモーメントを含んでいるが，モーメント次数の高い小行列の要素ほど低次小行列のそれに対して，観測信号ベクトルのノルム ||x(ω, t)|| の大小に伴って極端に値が増減する．観測信号をそのまま用いると，大きい信号成分が含まれる帯域とそれ以外の帯域では，複数次数モーメントの分析の基準が異なったものとなるため，モーメント同時分析を正しく行うためには異なる次数の要素全てを適切に重み付ける必要がある．本稿
では各周波数ビンにおいて，φd1 ,··· ,dm の要素中で最大次数となる dm を基準として，L2dm -ノルムによりこの重み付けを行う．

x0 (ω, t) =didj

(24)
(25)

l=1

ただし，x0 (ω, t) = [x01 (ω, t), · · · , x0M (ω, t)] は正規化された観測信号（正規化の規則は後述），小行列
Rdi ,dj (ω) は (di + dj ) 次クロスモーメント行列を表す．このとき行列 R⊕
d1 ,··· ,dm (ω) 全体と，単一で最大次数の小行列 Rdm ,dm (ω) との次元数を比較すると，そ∑m
れぞれ i=1 M di と M dm となる．次元数の差がどの程度の大きさかは φd1 ,··· ,dm の設計に依存するが，これは単一次数モーメントを分析する写像 MUSIC と比べて，提案手法がより大きな観測の次元拡張を行っていることを示している．さらに R⊕
d1 ,··· ,dm (ω) の分
析では，低次のモーメントを高次のモーメントと同時に分析することで，特に高次モーメントで問題となる統計的バイアスの影響が抑えられる．このように，φd1 ,··· ,dm を用いた提案手法による複数次数の高次モーメント同時分析は，従来の単一次数の高次統計量に基づいた分析に対して，より大きな次元拡張と統計的バイアスの抑圧を達成している．

Recording environment.

wx(ω, t)
√ []
∑M
2dm
12dm
E M
i=1 |xi (ω, t)|(26)ただし，w は L2dm -ノルム基準の正規化を調整するためのパラメータである. 提案手法による分析では，wが大きい場合は高次モーメントを重視したものとなり，w が小さい場合はその逆となる．

T

日本音響学会講演論文集5

実験

本章では，提案手法の有効性を検証するために行った実験について述べる．推定性能の比較のための従来手法として，写像 MUSIC (d = 2, 3) と 2qMUSIC (q = 2, 3) を用いた．

5.1

実験条件

図 1 に示す環境で収録した伝達関数による日本語音声の畳み込み混合を用いて，提案手法の性能評価を行った．観測信号には，拡散性雑音 [9] を模擬した
ピンクノイズを異なる 3 つの SNR で重畳した．優決定条件と劣決定条件それぞれの下での推定性能を検証するために，異なる 2 つの音源数で評価を行った．
他の実験条件については表 1 に示す．
方向推定の性能評価には，{0◦ , 30◦ · · · , 330◦ } の 12方向にランダムに配置した 3 または 5 音源を 100 通- 603 -2013年9月

20

Table 1 Experimental conditions# of sound sources
Sampling frequency
Room size
Reverberation
Speech duration
Frame length
Frame shift length
Window function

3, 5
16 kHzり用意し，真の音源方向と推定方向の間の 2 乗平均平方根誤差 (RMSE) を用いた．音源数分のピーク探索
結果が得られず正しく推定方向を得られなかった場合
は，全ての推定方向のエラーを平均したペナルティを与えた．また，(13)–(16) に示す写像空間の信号部分空間の次元数 N 0 は，3 音源では N 0 = 15，5 音源では N 0 = 30 を用い，(26) に示す正規化パラメータ wは w = 2 を用いた．なお，M ≤ N における MUSICの方向推定は，最小固有値に対応する一次元部分空間を雑音部分空間とみなして行った．

5.2

実験結果

図 2–3 に音源数ごとに得られた各手法の RMSE を示す．全ての条件の下で提案手法が最も優れており，
写像 MUSIC は 2q-MUSIC に対してやや良好な結果となっている．写像 MUSIC の次数 d = 2 と d = 3の結果を比較すると，前者の方が優れた性能を示し，2q-MUSIC の q = 2 と q = 3 の比較についても同様の傾向が現れている．これらの結果は，[3] で検証されているように，次数の高い統計量ほどより大きなバイアスを受けることによるもので，今回の実験に用いた時間長ではスナップショットの数が不十分であったと考えられる．他方，提案手法と写像 MUSIC の比較は，複数次数のモーメント同時分析の有効性を示唆している．今回の条件では，6 次モーメントのみを用いた分析 (d = 3) はやや劣った結果となっているが，提案手法の結果からこれを 4 次モーメント (d = 2) の補佐的な情報として利用することで，推定性能が向上することを確認できる．また，提案手法は次数の高いモーメントの方が優れた推定結果を示す十分なスナップショット数の信号を分析する場合でも，安定化した高次モーメントとともに低次モーメントの情報を利用することとなるため，やはり単一次数のモーメント分析よりも優れた推定結果が期待できる．

6おわりに

本稿では，劣決定条件に適用可能な新たな高分解能音源方向推定手法として，写像 MUSIC の拡張手法

日本音響学会講演論文集

RMSE [degree]

Sound sources

Circular array withradius of 0.1 m
Speakers 1.5 m
apart from array# of sensors4SNR20 dB
3.4 m × 4 m × 2.7 m
T60 of 0.3 s1s512 samples
256 samples
Hamming window

Mapped MUSIC⣔ิ2(d=2)

Mapped MUSIC
26.5(d=3)
2q-MUSIC
27.4(q=2)
2q-MUSIC
29.9(q=3)

105

Proposed method⣔ิ60Fig. 2

Experimental results for 3 sources within 1 sec.

40MUSIC⣔ิ1
Mapped MUSIC30
RMSE [degree]

Sensor array configuration

MUSICMUSIC15

22.6(d=2)

Mapped MUSIC

26.5(d=3)

202q-MUSIC

27.4(q=2)

2q-MUSIC

29.9(q=3)10

Proposed method⣔ิ60Fig. 3

Experimental results for 5 sources within 1 sec.

を提案した．従来の単一偶数次モーメント分析のための写像について示し，提案手法の方向推定アルゴリズムがこれを直和して得られる新たな写像に基づくことを述べた．提案手法による分析は，観測信号の複数次数の高次モーメント同時分析に相当しており，高次モーメントの分解能の高さと低次モーメントの推定における安定性の両方の利点を備えたものであることを述べた．単一の高次統計量分析に基づいた従来手法との推定性能比較実験によって，複数次数の高次統計量の同時分析が有効であり，提案手法により極めて高い推定性能が得られることを確認した．
謝辞 この研究は，科研費 23760325 の支援を受けたものである．

参考文献[1] R. O. Schmidt, IEEE Trans. Antennas Propag.,
vol.AP–34, pp.276–280, 1986.
[2] B. Porat et al., IEEE Trans. Signal Process.,
vol.39, pp.2016–2024, 1991.
[3] P. Chevalier et al., IEEE Trans. Signal Process.,
vol.54, pp.2986–2997, 2006.
[4] Y. Sugimoto et al., Proc. IWAENC, 2012.
[5] H. Saruwatari et al., IEICE Trans. Fundam.,
vol.E82-A, pp.1501–1510, 1999.
[6] H. Kamiyanagida et al., IEICE Trans. Fundam.,
vol.E84-A, pp.999-1010, 2001.
[7] S. Miyabe et al., Proc. ICASSP, pp.1641–1644,
2009.
[8] M. Wax et al., IEEE Trans. Acoust. Speech Signal Process., vol.32, pp.817–827, 1984.
[9] I. A. McCowan et al., IEEE Trans. Speech Audio
Process., vol.11, pp.709–716, 2003.

- 604 -

2013年9月