伊藤 信貴 
  
エマニュエル ・ ヴァンソン 
  
小野 順貴 
  
レミ ・ グリボンヴァル 
  
† 嵯峨山茂樹 
  
東京大学大学院 情報理工学系研究科 
  
Nobutaka Ito 
email: †ito@hil.t.u-tokyo.ac.jp
  
Emmanuel Vincent 
email: †emmanuel.vincent@inria.fr
  
Nobutaka Ono 
  
Rémi Gribonval 
email: remi.gribonval@inria.fr
  
Shigeki Sagayama 
email: sagayama@hil.t.u-tokyo.ac.jp
  
  
INFORMATION AND COMMUNICATION ENGINEERS 信学技報

Keywords: Source localization, covariance matrix, microphone arrays, MUSIC, diffuse noise 1

In this paper, we propose a method for estimating the azimuths of multiple sound sources accurately even in the presence of diffuse noise. MUSIC (MUltiple SIgnal Classification) for the estimation of the azimuths of multiple sources is robust against spatially white noise but the estimation performance degrades in the presence of diffuse noise. Based on a low-rank assumption on the covariance matrix of directional signals and the assumption that the covariance matrix of diffuse noise belongs to a subspace in a matrix space, the proposed method estimates the covariance matrix of the directional signals and applies MUSIC to the estimated matrix. The subspace model on the covariance matrix of diffuse noise includes as special cases noise models such as spatially uncorrelated noise, noise with a given coherence matrix, and isotropic noise observed with a crystal array. We showed through experiments with real-world noise recordings that the proposed method estimated the azimuths of multiple sources more accurately than the conventional MUSIC.

拡散性雑音環境における複数音源定位の問題 設定

本稿を通じて次の記法を用いる。上付き記号 * および H はそ れぞれ複素共役とエルミート転置を表す。信号は時間周波数領 域における零平均の確率変数で表され、τ および ω はそれぞれ フレーム番号と角周波数を表す。ベクトルの信号 ¸(τ, ω) の共 分散行列を Φ¸¸(τ, ω) , E[¸(τ, ω)¸H(τ, ω)] (1) により表す。但し、E[•] は期待値演算である。 M 素子のマイクロフォンアレイが拡散性雑音の存在下で、 到来方向未知の L(L < M )個の方向性信号を受信するとす る。ここで、これらの方向性信号のいくつかは方向性の干渉 音であっても良い。s(τ, ω) ∈ C L を基準点(アレイ配置の重 心など)で観測された L 個の方向性信号からなるベクトル、 x(τ, ω) ∈ C M を各マイクロフォンにおける観測信号からなる ベクトル、v(τ, ω) ∈ C M を各マイクロフォンにおける拡散性 雑音からなるベクトルとする。方位角 θ から水平に入射する平 面波のステアリングベクトルを d(ω; θ) , h e -jωδ 1 (θ) e -jωδ 2 (θ) . . . e -jωδ M (θ) i T (2) と記す。ここで、δm(θ) は方位角 θ の平面波が基準点から m 番 目のマイクロフォンに到達するまでに掛かる時間であり、マイ クロフォン座標を用いて幾何学的に計算できる θ の既知の関数 である。以上より、θ l を l 番目の音源の到来方位角とし、行列 D(ω) = h d(ω; θ1) d(ω; θ2) • • • d(ω; θL) i (3) を定義すると、観測信号のモデルは次のように書ける。 x(τ, ω) , D(ω)s(τ, ω) + v(τ, ω) (4) s(τ, ω) と v(τ, ω) が無相関であるとすると、共分散行列の間に 次の関係が成立する。 Φxx(τ, ω) = D(ω)Φss(τ, ω)D H (ω) + Φvv(τ, ω) (5) 本 稿 で は Φxx(τ, ω) を 観 測 さ れ た 共 分 散 行 列 と 呼 び 、 D(ω)Φss(τ, ω)D H (ω) を方向性信号の共分散行列、Φvv(τ, ω) を拡散性雑音の共分散行列と呼ぶ。 拡散性雑音環境下での複数音源定位の目的は、x(τ, ω) が与 えられたときに {θ l |l = 1, . . . , L} を推定することである。 3. MUSIC 3. 1 MUSIC による複数音源定位の原理 MUSIC による音源定位は、アレイ多様体と信号部分空間の 交点の検出に基づく。アレイ多様体とは θ を変化させたときに ステアリングベクトル d(ω; θ) が描く C M 内の曲線 M(ω) , {d(ω; θ)|θ ∈ [0, 2π)} (6) である。信号部分空間は真の音源方向に対応するステアリング ベクトルが張る C M の部分空間 S(ω) , span{d(ω; θ l )} L l=1 (7) である。θ = θ1, θ2, . . . , θL において曲線 M(ω) は S(ω) と交わ るから、このとき S(ω) の直交補空間である雑音部分空間 N (ω) への d(ω; θ) の正射影は零となる。従って、正射影の 2 ノルム の 2 乗の逆数 fN(ω; θ) , 1 V H (ω)d(ω; θ) 2 2 (8) は θ = θ1, θ2, . . . , θL においてピーク値をとる。ここで、V (ω) は N (ω) の正規直交基底を列として持つ M × (M -L) 行列で ある。本稿では fN(ω; θ) を狭帯域 MUSIC スペクトルと呼ぶ。 単一の方位角推定結果を得るためには、狭帯域 MUSIC スペ クトルを周波数方向に平均して単一の広帯域 MUSIC スペクト ルを計算する必要がある。平均の取り方は複数考えられるが、 例えば算術、幾何、調和平均を用いた場合はそれぞれ以下のよ うになる [3], [4]。 fW,A(θ) , 1 |Ω| X ω∈Ω fN(ω; θ) (算術平均) (9) fW,G(θ) , » Y ω∈Ω fN(ω; θ) -1/|Ω| (幾何平均) (10) fW,H(θ) , |Ω| P ω∈Ω 1/fN(ω; θ) (調和平均) (11) ここで Ω は平均に用いられる ω の値全体の集合である。方位 角の推定値は広帯域 MUSIC スペクトルのピークを検出するこ とにより得られる。 3. 2 雑音部分空間の同定 狭帯域 MUSIC スペクトル (8) の計算においては、雑音部分 空間 N (ω) の正確な同定が重要である。N (ω) は方向性信号の 共分散行列の零空間に一致するから、方向性信号のみが存在す る場合にはそれは容易に求まる。一般には雑音が存在するため に方向性信号の共分散行列自体は観測できないが、空間的に白 色な雑音は共分散行列の固有空間を変えないため、この場合に は観測された共分散行列の最小固有値に属する固有空間とし て N (ω) を求めることができる。しかしながら、拡散性雑音は マイクロフォン間で相関を持つため、共分散行列の固有ベクト ルの方向を変え、拡散性雑音の存在下では観測信号から直接 N (ω) を正確に同定することが困難である。 4. 共分散行列の雑音除去のための提案法 本節では、MUSIC を拡散性雑音に対して頑健にすることを 目指して、観測された共分散行列から方向性信号の共分散行列 を復元する方法を示す。これができれば、雑音部分空間 N (ω) は復元された行列の零空間として同定することができる。 拡散性雑音の共分散行列は一般に必ず観測できるとは限らな いが、空間的な無相関性など、なんらかのモデルが仮定できる 場合も多い。本研究では、こうしたモデルの多くが行列空間内 の部分空間として表現できることに着目し、まず観測された共 分散行列をこの雑音共分散行列部分空間の直交補空間に正射影 することにより、拡散性雑音の寄与を含まない成分を得る。こ の成分は、理想的には方向性信号の共分散行列のみの寄与で決 まるが、一方、雑音共分散行列部分空間に含まれる信号共分散 成分は得られない。そこで次に、方向性信号の共分散行列が低 ランク(音源数に一致)であることを仮定し、行列の要素(も しくはそれらの線形結合)の一部分から全体を復元する行列補 完と呼ばれる数理的技術を用いて、方向性信号の共分散行列を 推定する。 4. 1 拡散性雑音の共分散行列の部分空間モデル 本稿では、拡散性雑音の共分散行列が、M × M エルミート 行列の全体がなす R 上のベクトル空間 H ⊂ C M ×M の或る部 分空間 V(ω) に属すると仮定する。H には内積 A, B , M X m=1 M X n=1 amnb * mn (12) = tr(AB) (13) とノルム A , p A, A (14) 
が付随する。以下に示すように、この部分空間モデルは、空間 的に無相関な雑音 [5] やコヒーレンス行列が既知の雑音 [6], [7] といった雑音モデルを特別な場合として含む。 4. 1. 1 空間的に無相関な雑音 マイクロフォン同士が、波長と比較して十分に離れていれば、 拡散性雑音はマイクロフォン間で無相関であると近似できる。 このとき拡散性雑音の共分散行列は対角だから、このモデルは 次の部分空間に対応する。 V = {A ∈ H|A: 対角行列 } (15) このとき、V およびその直交補空間 V ⊥ への射影演算子 P[A] および P ⊥ [A] は以下で与えられる。 P[A] = diag(A) (16) P ⊥ [A] = off-diag(A) (17) ここで、diag[•] は非対角成分を 0 で置き換える演算、off-diag[•] は対角成分を 0 で置き換える演算である。 4. 1. 2 コヒーレンス行列が既知である雑音 球面対称な雑音場など、理想的な拡散性雑音場ではそのコ ヒーレンス行列(正規化された共分散行列)Γ(ω) が分かって いることがある。このような場合に対応する部分空間は次で与 えられる。 V(ω) = {kΓ(ω)|k ∈ R} (18) 
射影演算子は

P(ω)[A] = tr(AΓ(ω)) tr(Γ 2 (ω)) Γ(ω) (19) 
P ⊥ (ω)[A] = A - tr(AΓ(ω)) tr(Γ 2 (ω)) Γ(ω) (20) により与えられる。 4. 1. 3 結晶型アレイで観測された等方的雑音 我々は論文 [8] において、結晶型アレイ(図 1 参照)と呼ば れる一定のクラスの対称なアレイ配置を用いるならば、等方的 雑音の任意の共分散行列 Φvv(τ, ω) は定ユニタリ行列 P によ り対角化されることを示した。ここで、等方的雑音とは次の 2 つの条件を満たす雑音のことと定義する。 • パワースペクトログラムは全てのマイクロフォンにおい て同一である: [Φvv]11(τ, ω) = [Φvv]22(τ, ω) = • • • = [Φvv]MM (τ, ω) (21) • チャネル間のクロススペクトログラムは、等距離のマイ クロフォン対に対しては等しくなる: rmn = rpq ⇒ [Φvv]mn(τ, ω) = [Φvv]pq(τ, ω) (22) ここで rmn は m 番目および n 番目のマイクロフォン間の距離 である。 図 1 結晶型アレイの例 上記のモデルは H の部分空間として下記のように書きなおせ る。Φvv(τ, ω) が P により対角化されるというのは、Φvv(τ, ω) が部分空間 V = {A ∈ H|P H AP : 対角 } (23) に属することと同値である。V はより具体的に V = span{pmp H m } M m=1 (24) とも書ける。ここで、pm は P の第 m 列である。 射影演算子 P および P ⊥ は下記のように導くことができる。 {pmp H m } M m=1 が V の正規直交基底をなすことに注意すると、 P[A] = M X m=1 A, pmp H m pmp H m (25) = M X m=1 tr(Apmp H m )pmp H m (26) = M X m=1 (p H m Apm)pmp H m ( 27 
)
= P diag(P H AP )P H (28) 

P ⊥ [A] = A -P[A] (29) = P off-diag(P H AP )P H (30) 4. 2 方向性信号の共分散行列の復元アルゴリズム Φxx をノイズフリーな空間 V ⊥ に射影するために、V の 正規直交基底を Q1, • • • , Qp と書き、V ⊥ の正規直交基底を Qp+1, • • • , Q M 2 と書く。但し、p , dimV である。Φxx の {Qi} M 2 i=1 による展開を Φxx = M 2 X i=1 φiQi (31) と書くと、Φxx の V ⊥ 成分は P ⊥ [Φxx] = M 2 X i=p+1 φiQi (32) である。但し、 φi , Φxx, Qi (33) である。 射影 P ⊥ [Φxx] は理想的にはノイズフリーであり、方向性信 号の共分散行列の V ⊥ 成分に一致するが、一方、Φxx の V 成 分は拡散性雑音の共分散行列の寄与も含む。そこで、V 成分を 未知数 {zi} p i=1 を用いて表した行列 Y , p X i=1 ziQi + M 2 X i=p+1 φiQi (34) を考え、これを Y = Θ + E (35) のように、真の方向性信号の共分散行列 Θ の誤差 E を含む観 測であるとみなし、低ランク制約 Θ ∈ Ψ の下、Θ を最尤推定 する。ここで、Ψ はランク L 以下の M × M 半正定値エルミー ト行列全体の集合である。 いま、 Θ = M 2 X i=1 θiQi (36) E = M 2 X i=1 ẽiQi (37) により係数 { θi}、{ẽi} を定義すると、式 (35) は zi = θi + ẽi (i = 1, . . . , p) (38) φi = θi + ẽi (i = p + 1, . . . , M 2 ) (39) と書き直せる。最尤法を定式化するために、誤差 ẽi は i.i.d. でお のおの複素ガウス分布 N (0, σ 2 ) に従うと仮定する。z1, • • • , zp は未知であるから隠れ変数とみなし、{ φp+1, • • • , φM 2 } を観測 (不完全)データ、{z1, • • • , zp, φp+1, • • • , φM 2 } を完全データ とみなし、制約条件 Θ = M 2 X i=1 θiQi ∈ Ψ (40) の下、観測データの尤度を最大化する { θi} M 2 i=1 を EM アルゴリ ズムにより探索する。 完全データの対数尤度は ln p({zi} p i=1 , { φi} M 2 i=p+1 ; { θi} M 2 i=1 , σ 2 ) , -M 2 ln σ 2 - p X i=1 |zi -θi| 2 σ 2 - M 2 X i=p+1 | φi -θi| 2 σ 2 (41) である。これを隠れ変数に関して平均化して、Q 関数は、 Q({ θi} M 2 i=1 , σ 2 ; { θ(t) i } M 2 i=1 , σ2(t) ) = -M 2 ln σ 2 -p σ2(t) σ 2 - p X i=1 | θ(t) i -θi| 2 σ 2 (42) - M 2 X i=p+1 | φi -θi| 2 σ 2 となる。これを最大化する { θi} M 2 i=1 (つまり、{ θ(t+1) i } M 2 i=1 )は、 次式を最小化する。 p X i=1 | θ(t) i -θi| 2 + M 2 X i=p+1 | φi -θi| 2 (43) = Ŷ (t+1) -Θ 2 F (44) ここで、 Ŷ (t+1) , p X i=1 θ(t) i Qi + M 2 X i=p+1 φiQi (45) = P[ Θ(t) ] + P ⊥ [Φxx] (46) 
Θ(t) , M 2 X i=1 θ(t) i Qi (47) である。よって、結局 Θ(t+1) = arg max Θ∈Ψ Ŷ (t+1) -Θ 2 F (48) とすれば良い。 紙面の都合上証明は省略するが、式 (48) の解は以下のよう に Ŷ (t+1) のある種の打ち切り固有値分解により与えられる。 Ŷ (t+1) の固有値分解を Ŷ (t+1) = U (t+1) Λ (t+1) U H(t+1) (49) 
とおく。ここで U (t+1) はユニタリ行列である。式 (48) の解は Θ(t+1) = U (t+1) Λ (t+1) T U H(t+1) (50) である。ここで、Λ (t+1) T は Λ (t+1) の対角要素( Ŷ (t+1) の固有 値)を大きい方から L 番目まで残して、残りは 0 で置き換える 打ち切りを行った後、更に負の固有値があればそれを 0 で置き 換えることにより得られる行列である。 以下にアルゴリズムをまとめる。 • Ŷ (0) = Θ(0) = Φxx • t = 0, 1, . . . に対し、以下をループ停止条件が満たされ るまで反復する。 -Ŷ (t+1) , P[ Θ(t) ] + P ⊥ [Φxx] -Y (t+1) の固有値分解 Ŷ (t+1) = U (t+1) Λ (t+1) U H(t+1) を 計算する。 -Θ(t+1) = U (t+1) Λ (t+1) T U H(t+1) なお、観測された共分散行列のみから信号共分散行列を一意 に復元できるためには、信号共分散行列のランクが一定以下で ある必要があると考えられる。そのようなランクの上限につい ては今後の検討課題としたい。

実験による性能評価

Θ(t) -Θ(t-1) F Θ(t) F < 10 -4 (51) とした。なお、目的音声は ATR 日本語データベース [11] から 取った。 方向性信号の共分散行列推定の精度を確認するために、図 2 にフロベニウス・ノルムで測った真の方向性信号の共分散行列 Θ に対する相対誤差 Θ -Θ F Θ F ( 

  位実験の結果を示す。直径 5 cm の正方形アレイを用いて新宿駅 で収録した雑音[START_REF] Ito | Designing the Wiener post-filter for diffuse noise suppression using imaginary parts of inter-channel cross-spectra[END_REF] に、鏡像法[START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF] m×3 m とし、アレイは部屋の中央に水平に配置 するものとし、音源はアレイと同一水平面内、1 m の距離にあ るものとし、音源数は 2 とした。1 番目のマイクロフォンにおけ る観測信号の SN 比は、いずれの音源についても -15 dB とし た。観測信号の信号長は 10 秒、サンプリング周波数は 16 kHz であった。観測信号の分析には短時間フーリエ変換を用い、フ レーム長を 512 点、フレームシフトを 64 点とし、ハミング窓 を用いた。Φxx の計算には全フレームを用い、x(τ, ω)x H (τ, ω) を時間平均することにより求めた。このように計算された Φxx に MUSIC を適用した場合(以下、従来法)と、雑音除去され た共分散行列に MUSIC を適用した場合(以下、提案法)で音 源定位の性能を比較する。拡散性雑音の共分散行列のモデル部 分空間としては、コヒーレンス行列が既知の雑音 (18) を用い、 コヒーレンス行列は球面等方場に対するものを用いた [7]。反復 の停止条件は

図 3

 3 相対誤差である。音源の方位角は 60 • 、330 • であり、残響時間 は 130 ms であった。図より、共分散行列の雑音除去により、真 の方向性信号の共分散行列との誤差が著しく減少していること が分かる。 において、従来法と提案法による広帯域 MUSIC スペ 図 3 広帯域 MUSIC スペクトルの例。破線:従来の MUSIC、実線: 提案法、直線:真の音源の方位角。 図 4 平均二乗誤差の平方根。破線:従来の MUSIC、実線:提案法。 クトルを比較する。諸条件は図 2 と同一である。SN 比の著し く低い 3 番目の帯域(63 Hz)までは捨て、4 番目以降の帯域 における狭帯域 MUSIC スペクトルを幾何平均した。図中の直 線は真の音源の方位角を示す。従来法(破線)は 60 • の音源に 対応するピークが鈍く、ピーク位置が大きくずれているのに加 え、偽のピークも存在する。それに対し、提案法(実線)によ る MUSIC スペクトルにおいては、比較的正確な位置に鋭い ピークがあり、偽のピークも抑制されていることが分かる。 最後に、従来法と提案法による方位推定精度の統計的な比較 を行った。図 4 は平均二乗誤差の平方根(root mean square error)を残響時間に対してプロットしたものである。集合 {0 • , 60 • , . . . , 300 • } から抽出した任意の 2 つの方位角に対して 実験を行った結果を平均して計算した。