2-10-4
非同期録音ブラインド同期のための線形位相補償の効率的最尤解探索∗◎宮部滋樹 (筑波大), 小野順貴 (NII), 牧野昭二 (筑波大)1

はじめに非同期マイクロホンアレー [1, 2] は，複数の独立した録音機器を用いた非同期の多チャネル録音でアレー信号処理を行う枠組で，多チャネルの音声強調や音響モニタリングなどを，携帯電話やボイスレコーダといった汎用的な録音装置を組み合わせた簡易かつ柔軟な機器構成で行える利点がある．しかし，各チャネルの録音が同期していないために録音開始時刻やサンプリング周波数の不一致が起こり [1, 3]，アレー信号処理の性能が大幅に劣化するという問題がある．
我々は前回の発表 [4] において，ブラインド音源分離(BSS) [5] の前処理のために，サンプリング周波数のチャネル間ミスマッチを短時間フーリエ変換 (STFT)領域で最尤法により推定し補償する方法を提案した．
しかし最尤推定が解析的に解くことができないものとなるため，計算量の大きい離散値全探索を行う必要があった．本稿では，最適値付近では局所的に短報的になるという経験的に得られた性質を利用し，荒い離散値全探索による探索範囲の絞り込みと，黄金分割探索による詳細な探索を併用した，効率的な最尤解探索法について議論する．

2

時間周波数領域におけるサンプリング周波数ミスマッチの最尤推定

本節では [4] で提案したサンプリング周波数ミスマッチ補償のための最尤推定について述べる．ここでは 2 素子マイクロホンアレーについて議論するが，3素子以上への拡張は，基準チャネルを定めて基準以外のチャネルごとに補償することにより容易に行える．
サンプリング周波数ミスマッチの時間領域表現

2.1

同時刻における 2 つのマイクロホンの連続信号
x1 (t) , x2 (t)（t は連続時間）が別々の A/D 変換器でサンプリングされて離散信号 x1 (n) , x2 (n)（n はサンプル番号）が得られたとする．ここで x1 (n) のサンプリング周波数は fs ，x2 (n) のサンプリング周波数は未知の無次元数 ϵ（|ϵ| ≪ 1）により表される (1 + ϵ) fsであるとする．このとき各チャネルの離散信号と連続信号の関係は以下のように表される．
( )nx1 (n) = x1(1)fs()
nx2 (n) = x2+ T21(2)(1 + ϵ) fs
ここで連続時刻 t の原点は x1 (n) の録音開始時刻とし，T21 は x2 (n) の録音開始時刻とする．従って同じ
時刻 t を参照する第 i チャネル（i = 1, 2）の離散時∗

Window function

First sampleof the frame2


Frame centralsample ଵ

2
 1

Discrete time

Frame length Fig. 1 Relation between the frame central sample
m of a frame and the sample n1 inside the frame.
刻 ni の関係は

n2 = (1 + ϵ) n1 − fs T21(3)となり，長い時間が経過して n1 , n2 が大きくなるのに比例して参照する時刻の差が拡大する．

2.2

時間周波数領域におけるサンプリング周波数ミスマッチのモデル化

アレー信号処理の多くは STFT 領域で行われるため，サンプリング周波数ミスマッチ補償は STFT 領域での信号表現がよく近似する方法であれば十分であると考える．そのためにまず STFT のフレーム長を L，チャネル 1 のあるフレームの中心サンプルをm として，m − L/2 ≤ n1 ≤ m + L/2 − 1 という 1 フレーム内での時間の対応関係を考える（Fig. 1）．
式 (3) の関係から，サンプル m の近傍の n1 には

n2 = (1 + ϵ) (n1 − m) + (1 + ϵ) m − fs T21⇋ (n2 − m) = (1 + ϵ) (n1 − m) + ϵm − fs T21(4)という関係が成り立ち，(n2 − m) と (n1 − m) の対応関係は，m とともに ϵm だけ拡大していくことがわ
かる．|ϵ (n1 − m)| は第 2.3 節の処理によって 0 に近い値となるため，ここでは ϵ (n1 − m) を無視すると，(n2 − m) ≈ (n1 − m) + ϵm − fs T21= (n1 − m) + τ21 (m; ϵ)τ21 (m; ϵ) = ϵ (m − M )fs T21M=ϵ(5)(6)(7)のように，フレーム内で時間差が n1 に依存せず一定と仮定したモデルを得る．ここで τ21 (m; ϵ) はフレームの時刻オフセットであり，M は両チャネルが同じ連続時刻を指して n1 = M のとき n2 = M となる離散時刻で，以下ではオフセット原点と呼ぶこととする．
もし全てのフレームの中心サンプル m に対して

|τ21 (m; ϵ)| = |ϵ (m − M )| ≪ L(8)

Eﬃcient Maximum Likelihood Solution Search of Linear Phase Compensation for Blind Synchronizationof Asynchronized Recording. by MIYABE, Shigeki (University of Tsukuba), ONO, Nobutaka (NII),
MAKINO, Shoji (University of Tsukuba)

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

2013年3月

が満たされるとき，オフセット τ21 (m; ϵ) の補償はフレーム内の時間シフトとみなすことができ，これはSTFT 領域では線形位相により近似的に補償することができる．観測信号 x1 (n) , x2 (n) を一律でフレーム分析に施して得られる第 m サンプルを中心とするフレーム波形の短時間フーリエ変換 Xi (k, m) , i =1, 2, −L/2 < k ≤ L/2 を以下のように表される．

Xi (k, m) =L−1∑(w (l) xi

l=0Ll+m−2)(2πȷklexp −L)(9)
ここで w (l) は長さ L の再合成可能な窓関数，k =−L/2 + 1, . . . , L/2 は離散周波数インデックス，ȷ =√−1 を表す．第 2 チャネルの信号 X (k, m) にオフセッ
ト τ21 (m; ϵ) を補償する線形位相を与えることにより，補償した信号 X̂2 (k, m) は以下のように得られる．
()
2πȷkτ21 (m; ϵ)X̂2 (k, m) = X2 (k, m) exp
L()
2πȷkϵ (m − M )= X2 (k, m) exp
(10)
L

2.3

観測信号全体の相互相関を用いた録音開始時刻
ずれの大まかな補償式 (2) に現れる録音開始時刻の差 T21 は，観測信号x1 (n) と x2 (n) のみを用いて推定することは難しい．
しかしながら，BSS のような音源方位情報を明に利用しないクラスのアレー信号処理では，小さなチャネル間の固定時間差は許容され，正確な録音開始時刻差 T21 の推定は不要となる．そのためここでは，第2.2 節で述べた線形位相補償モデルが成り立つための条件である，|ϵ (n1 − m)| が全てのフレーム中心サン
プル m おいて無視できるようにするための，大まかな録音時刻差 T21 の補償を行う．
サンプリング周波数のミスマッチが十分小さく |ϵ| ≪1 ととなる場合は，サンプリング周波数のミスマッチの影響を受けていても x1 (n) と x2 (n) の相関は高いと仮定できる．そこで以下のように相関を最大にする x2 (n) の遅延量 δ12 を求める．

∑min(N1 ,N2 +δ)−1

δ12 = arg max−N2 <δ<N1x1 (n) x2 (n − δ)

x2 (n) ← x2 (n − δ12 )(12)とすることにより，x1 (n1 ) と x2 (n2 ) のサンプリング周波数ミスマッチ原点を信号のオーバーラップするサンプル区間の中央付近に移動する．また式 (10) 中のオフセット原点 M には，以下のようにオーバーラッ
プの中央付近のサンプル番号を推定として与える．

min (N1 − δ12 , N2 ) − max (0, δ21 ) − 1⌋
2(13)

日本音響学会講演論文集

2.4

サンプリング周波数ミスマッチ推定を評価する
尤度の定式化観測されるすべての音源は定常かつ位置の移動が無いと仮定すると，正確な ϵ の推定を用いてサンプリング周波数のミスマッチを補償した観測信号[]TX̂ (k, m; ϵ) = X1 (k, m) , X̂2 (k, m; ϵ)(14)は離散周波数 k 毎に定常となる．正確な ϵ の推定による補償で定常性を回復した観測信号 X̂ (k, m; ϵ) が共分散行列 V (k) の零平均多変量複素正規分布に従うと仮定すると，その対数尤度は∑(H−1J (V, ϵ) =− X̂ (k, m; ϵ) V (k) X̂ (k, m; ϵ)k,m)− log 2π 2 − log det V (k)(15)Hと表せる．ここで {·} は複素共役転置を表し，また共分散行列 V (k) は尤度関数を決定するパラメタであり，以下のように X̂ (k, m; ϵ) を用いた標本推定により推定する．
1 ∑HV (k) ←X̂ (k, m; ϵ) X̂ (k, m; ϵ)
(16)
|∀m|
∀m

これを式 (15) に代入すると第 1 項が定数となるため、対数尤度関数は定数項を除いて以下のように単純化される。
∑∑
Hlog detX̂ (k, m; ϵ) X̂ (k, m; ϵ)J (ϵ) = −km(17)サンプリング周波数ミスマッチ ϵ の補償が不正確であれば，時刻のドリフトによる定常性の低下により対数尤度 J (ϵ) が小さくなるため，J (ϵ) の最大化によりϵ を推定することができる．しかしこの尤度最大化は解析的に解くことができないため，次節では効率的な ϵ の最尤解の探索方法について議論する．

n=max(0,δ)(11)
ここで Ni , i = 1, 2 は xi (n) のサンプル数を表す．そして x2 (n) を δ12 だけ遅延させてM ←⌊ここで ⌊·⌋ は床関数を表す．この全サンプルのシフトとオフセット原点の推定により，式 (8), (10) のϵ (m − M ) は全ての m において可能な限り 0 に近い値となる．

3黄金分割法による最尤解の効率的探索

前節まででは [4] で提案したサンプリング周波数のミスマッチ ϵ を評価するための対数尤度関数 J (ϵ) を議論したが，この対数尤度関数を最大化する ϵ は解析的に求めることができない．推定するパラメタはϵ のみであり，その最適化には離散値全探索を行うことも考えられるが，一つの ϵ の評価のために全帯域での共分散行列とその逆行列の計算が必要であるため，高い解像度の離散値全探索を行うためには計算量が膨大になる．本節ではこの最尤推定問題の効率的な解探索法について議論する．
この最尤推定問題で求めるべきパラメタは ϵ のみであるため，一次元の最適化問題の代表的な解法で- 734 -

2013年3月

x 10

4Step 3: J (p) ≤ J (q) なら-3-4a = p, p = q, q = a + φ (b − a)(21)

Log likelihood function-5とし，そうでなければ以下のように更新する．

-6b = q, q = p, p = b − φ (b − a)-7

Step 4: (b − a) が十分に小さくなければ step 2 に戻り，十分に小さければ以下のようにとして ϵ の最尤推定値を求め，終了する．

-8-9-10-11-5(22)x 10-4-3-2-101Mismatch

234-4Fig. 2 An example of the log likelihood function
J (ϵ). We can see that the function is locally convex
in the neighbor of its maximum.
ある黄金分割探索 [6] を用いることが考えられる．黄金分割探索は関数の最大値または最小値を繰り返し探索により探索範囲を狭めながら求める手法で，関数が局所的に単峰的な範囲では最適解に一意に収束する．図 2 に示す例のように，対数尤度関数 J (ϵ) は最大値周辺では単峰性関数となることが経験的にわかっているため，適切に探索範囲を絞り込むことで黄金分割探索の利用が可能となる．
黄金分割探索の前にまず刻みの荒い離散値全探索により探索範囲を絞り込む．サンプリング周波数ミスマッチ ϵ を範囲 [−E, E] で等間隔に D 分割した

ϵd = −E +2dE
, d = 0, . . . , D − 1
D−1(18)について，すべての J (ϵd ) を比較して最大値を与える分割番号 d∗ を求める．

d∗ = arg max J (ϵd )(19)

d=0,...,D−1

この離散値全探索の範囲 [−E, E] は，機器間のサンプリング周波数ミスマッチとして妥当な範囲に設定すればよい．一般的な機器のサンプリング周波数ミスマッチは 10−5 オーダーであるといわれているため，E は 10−4 またはその数倍が妥当といえる．また，探索範囲の分割数 D は，音声を観測する場合なら 10 から 100 程度で適切に単峰な範囲を絞り込むことができる．
次に，[ϵd∗ −1 , ϵd∗ +1 ] を探索範囲とした ϵ についての J (ϵ) の黄金分割探索を以下のように行う．

Step 1: 探索範囲 [a, b] の初期値を a = ϵd∗ −1 , b =ϵd∗ +1 と定める．
Step 2: p = b − φ (b − a), q = a + φ (b − a) としてJ (p)，J (q) を求める．ここで φ は以下で与えられる．
√5−1φ=(20)2日本音響学会講演論文集

ϵ=5a+b2(23)

黄金分割探索は探索範囲を繰り返し毎に φ ≈ 0.62倍に狭めることができるため，離散値全探索で P 回のlog 1
関数評価が必要な P 点等分割の精度の探索は，log Pφ ≈2.08 log P 回の関数評価で達成することができる．

4

実験

提案手法の性能と非同期アレー信号処理への有効性を検証するため，複数話者の音声の混合のマイクロホンアレーによる観測信号に人工的にサンプリング周波数のミスマッチを与え，提案手法にるサンプリング周波数ミスマッチ補償の計算効率，補償精度とブラインド音源分離の性能への寄与を評価する．

4.1

実験条件

使用した観測信号は，２人の話者による発話に実測したインパルス応答を畳み込んで混合したマイクロホンアレー観測信号で，片方のチャネルのサンプリング周波数を 100 タップのポリフェーズフィルタを用いて人為的に変更してサンプリング周波数のミスマッチを模擬した．変更前のサンプリング周波数は16,000 Hz で，サンプリング周波数の変更は 16, 000 ±0.5，16, 000 ± 1，16,000±1.5 Hz の 6 種類とした．これらはそれぞれ ±31.25，±62.5，±93.75 ppm に相当
し，別々の A/D 変換器を用いた場合のサンプリング周波数のミスマッチとして現実的な大きさのものである．音源分離評価のための分離手法には補助関数
法独立ベクトル分析 [7] を用いた．その他の実験条件
を Table 1 に示す．

4.2

計算量の評価Table 2 に，10 秒の観測信号を提案手法の E =5 × 10−4 , D = 10 とした分割の荒い離散値全探索の範囲絞り込みの後に，絞り込まれた範囲を P 点の解候補に等分割して比較する離散値全探索と，この P点解候補探索と同等の精度の黄金分割探索の計算量を比較する．表より，提案手法は高い解像度でも計算量をあまり増やすことなく探索を行えていることがわかる．
以下の実験では，E = 5 × 10−4 , D = 10 の離散値全探索による探索範囲絞り込みと，Step 4 における収束条件を (b − a) < 10−9 とした，P = 105 等分割の解候補探索に相当する黄金分割探索用いて評価を行う．

- 735 -

2013年3月

2.50E-06

Experimental conditions

Speech signal

Concatenated word
uttrances by 4 speakers
3, 5, 10, 20, 30 seconds
T60 of 130 ms
4,096 samples
2,048 samples
1.5 m[-50◦ , 30◦ ], [-60◦ , -10◦ ]2 cm
Opteron 2.8 GHz

Table 2102

103

104

105

Discrete searchGolden section

10.17
2.52

85.52
2.96

339.70
3.31

3337
3.74

4.3

サンプリング周波数ミスマッチの推定精度

複数の長さのデータに対するサンプリング周波数の補正精度を平均二乗誤差 (RMSE) で評価した結果を図 3 に示す．最も短い 3 秒の観測信号でも RMSE
は元のサンプリング周波数ミスマッチ ϵ の 10 分の 1以下に収束し，データが増えるにつれて急速に小さくなる．従って観測信号の定常性を仮定した尤度がサンプリング周波数のミスマッチの評価尺度として有効であるということがわかる．

4.4

音源分離への寄与

サンプリング周波数のミスマッチを提案手法により補償することで，音源分離の性能が回復することを確かめるための分離精度評価を行った．分離フィル
タの学習には与えられた観測信号全体を用いた．評価尺度には，非目的成分の抑圧の強さを表す信号対干渉比 (signal-to-interference ratio; SDR) を用いた
[8]．また，SDR を算出するための参照信号としては，サンプリング周波数の変更を施していないマイクロホンにおける各音源の音像を用いた．
実験結果を図 4 に示す．まず Unprocessed が非常に低い値を示していることから，この条件ではサンプリング周波数ミスマッチの補償をしなければ音源分離ができない厳しい条件であるということがわかる．
手動でパラメタ ϵ の正しい値を与えた位相補償はサンプリング周波数ミスマッチがない場合よりも SDR が2 dB 程度低いだけであり，BSS のためのサンプリング周波数ミスマッチ補償に STFT 領域における位相補償が有効であるということがわかる．また，サンプ
リング周波数ミスマッチをブラインドに推定して補償する提案手法は，手動でパラメタ ϵ を与えた場合とほとんど性能差がなく，提案手法の最尤推定は性能限界に近い高い精度でサンプリング周波数ミスマッチ ϵ

日本音響学会講演論文集

1.00E-06
5.00E-07

051015202530

Signal length [s]Fig. 3 Root mean squared error (RMSE) of estimated sampling frequency mismatch ϵ.
20
No mismatch

Comparison of computation time in secondP

1.50E-06

0.00E+00

SDR [dB]

Signal length
Rverberation time
Frame length L
Frame shift M
Source distanceSource directions
Microphone spacingCPU

2.00E-06

RMSE ofTable 115

Unprocessed10

Linear phase compensationof ground truth of5Linear phase compensation
with maximum likelihood
estimate of001020
Time [s]30Fig. 4 Signal to distortion ratios (SDRs) of BSSresults.
を推定できているということを示している．以上より，音源分離のためのサンプリング周波数ミスマッチ補償として提案手法が有効であるということが確認された．

5

おわりに

本稿では，文献 [4] において提案した，非同期マイクロホンアレーのための観測チャネル間サンプリン
グ周波数ミスマッチのブラインド補償を，効率的な計算で達成する手法について述べた．STFT 領域のサンプリング周波数ミスマッチの影響をフレームの位置に線形な位相シフトとモデル化し，観測される音源が定常で移動がないと仮定した最尤推定は，解が解析的に得られないために計算量の大きな離散値全探索が必要であった．そこで尤度関数が大域最適解周辺で単峰性を示すことに着目し，荒い離散値全探索による範囲の絞り込みと黄金分割探索による高速な推定を用いる効率的な探索法を提案した．黄金分割探索の計算効率を評価し，またサンプリング周波数ミスマッチ補償のブラインド音源分離への寄与を確認した．

参考文献[1] Liu, Proc. IWAENC, 2008.
[2] Ono et al., Proc. WASPAA, 161-164, 2009.
[3] Markovich-Golan et al., Proc. IWAENC, 2012.
[4] 宮部ほか, 音講論 (秋)，689-670，2005．
[5] Makino et al., “Blind Speech Separation,”
Springer, 2007.
[6] Gill et al., “Practical Optimization,” Academic
Press, 1981.
[7] Ono, Proc. WASPAA, 189-192, 2011.
[8] Vincent et al., Proc. ICA, 552-559, 2007.

- 736 -2013年3月