2013.10.11EA@NTT CS Lab in Keihanna

時空間周波数領域クロストークキャンセラを用いた複数平行直線アレイによる音空間の収録と再生Sound ﬁeld recording and reproduction by multiple parallel linear arrays
using crosstalk cancellers in spatio-temporal frequency domain

岡本 拓磨，榎本 成悟，西村 竜一
情報通信研究機構

1

Presentation contents
Introduction
Purpose
複数平行直線マイクロホン / スピーカアレイを用いた大画面映像システムに最適な3次元音場収録・再現

Basic theory
Spectral division method

Proposed method
時空間周波数領域クロストークキャンセラの導出
計算機シミュレーション

Concluding remarks
2Introduction
大画面映像システム

200インチ裸眼立体システム
立体ディスプレイ・高精細ディスプレイ(4K，8K)
大画面：没入感
3D or 高精細：奥行き感，飛び出し感

大画面映像システムに最適な3次元音場収録・再現
映像と重らないスピーカ配置
没入感，奥行き感，飛び出し感を如何に制御するか!?

3

Purpose : Umekita acoustics
平行2直線スピーカアレイによる大画面システムに最適な音場再現

バーチャル音源の提示 (for rendered contents)
バーチャル音源の高さや奥行き，飛び出し音の制御
実環境音場の収録・再現 (for recorded contents)
最適なマイクロホン配置，信号処理

岡本ら，音講論(秋) 20134

How to arrange microphones?
裸眼立体ディスプレイの視聴位置の特徴

水平方向
複数人が同時に違う位置から見る立ち位置によって見え方がかわるため視聴者は左右に動く→ 最適視聴距離と平行に密に配置
仰角方向
頭の高さは基本的には画面中央付近→ 画面中央付近に密に配置

5

Proposed arrangement
複数平行直線アレイを用いた音空間の収録と再生
z

複数直線マイクロホンアレイ
水平方向は広く制御するために直線アレイを密に配置

Original sound ﬁeldx

仰角方向の音空間を収録するため複数の直線アレイを配置onhp
roic

rm
aen
Liz

複数直線スピーカアレイ
画面と重ならないように複数の直線スピーカアレイを配置

r loaLine

pe
uds

akrra
er asrayraeyys

Reproduced sound ﬁeld

xreen
csig

BCon

tr

intsopoly

6

Problem
一致zz
Linxe

ak
spe

udolr
ea

Original sound ﬁeld

ysar
rra
Reproduced sound ﬁeld

xn

ree
csg

Bi

Linn

horop

icmearrayr
aes
trnoCy(a) 音場収録部
収録信号

tsoinploy(b) 音場再現部
変換フィルタ↑
如何に計算するか!?

再生信号

7

Purpose従来法：最小二乗解(least squares method : LS method)：順問題↓
D(!) = G+ (!)P (!) ：逆問題

変換フィルタ
問題点
行列 G(!) は非常に悪条件→ 変換フィルタは極めて不安点
計算量がチャネル数の2乗オーダー

G(!)
..
.

D(!)

..
.

..
.

..
.

Microphones

P (!) = G(!)D(!)

loudspeakers

各スピーカ-マイクロホン間の伝達関数行列の逆特性を直接算出

P (!)

本研究の具体的な目的
複数平行直線アレイに最適な変換フィルタの設計
時空間周波数領域クロストークキャンセラの提案

8Basic theory (1)直線スピーカアレイを用いた水平面音場再現
波面合成法(Wave ﬁeld synthesis : WFS)
直線アレイによる音場(Single layer potential)
Z 1
P (x, !) =D(x0 , !)G3D (x x0 , !)dx01

3次元自由空間グリーン関数
G3D (xexp( jk|x x0 |)
x0 , !) =4⇡|x x0 |位置 xでの音圧 P (x, !)↓スピーカの駆動信号 D(x0 , !)と伝達関数 G3D (xx0 , !) との畳み込み

9Basic theory (2)直線マイクロホン / スピーカアレイを用いた水平面音場再現
Spectral division method (SDM)x 軸方向に空間フーリエ変換

J. Ahrens et al. 2010.

P̃ (kx , y, 0, !) = D̃(kx , !) · G̃(kx , y, 0, !)→ 時空間周波数領域ではかけ算

3次元グリーン関数の空間フーリエ変換
8
⇣p⌘
(2)
j22y>
Hkk>
x04>
>
>
>
for 0  |kx |  |k|<G̃(kx , y, 0, !) =⌘⇣p>
>
1
22y>
K
kk>
0x2⇡
>
>
:for 0 < |k|  |kx |時空間周波数領域の各スピーカの駆動信号P̃ (kx , yref , 0, !)D̃(kx , !) =G̃(kx , yref , 0, !)変換フィルタ10

Basic theory (3)直線マイクロホン / スピーカアレイを用いた水平面音場再現の実装
実装のためには離散化P̃ (kx,i , y, !) = D̃(kx,i , !) · G̃(kx,i , y, !)(i = 1, 2, · · · , K)

アレイの打ち切り誤差 → 窓関数の導入

D̃(kx , !)

D(!)

d(t)

..
.

Loudspeaker array

Tapering window

IDFT

Spatial IDFT

SDM
D̃(kx , !) =P̃ (kx , yref , !)G̃(kx , yref , !)Spatial DFT.
..

.
..

DFT.
..

P̃ (kx , !)

.
..

P (!)

.
..

.
..
p(t)

.
..

.
..

Microphone array

実装例

11

Strategy for proposed method
複数平行直線アレイの時空間領域信号( M = 2, N = 2 )各アレイ間の違い=距離の違いz
D̃1 (kx,i )Linear loudspeaker array 1( K ch)

r1,1r1,2⊗
x r2,1D̃2 (kx,i )Linear mic. array 1( K ch)P̃1 (kx,i )yr2,2Linear mic. array 2( K ch)P̃2 (kx,i )Linear loudspeaker array 2( K ch)P̃ (kx,i , y, !) = D̃(kx,i , !) · G̃(kx,i , y, !) の重ね合わせP̃1 (kx,i ) = G̃(kx,i , r1,1 )D̃1 (kx,i ) + G̃(kx,i , r2,1 )D̃2 (kx,i )P̃2 (kx,i ) = G̃(kx,i , r1,2 )D̃1 (kx,i ) + G̃(kx,i , r2,2 )D̃2 (kx,i )12

Proposed method
時空間周波数領域クロストークキャンセラ
複数平行直線アレイの時空間領域信号の行列表記P̃ (kx,i ) = [P̃1 (kx,i ) P̃2 (kx,i )]TG̃(kx,i , r1,1 ) G̃(kx,i , r2,1 )G̃(kx,i ) =G̃(kx,i , r1,2 ) G̃(kx,i , r2,2 )P̃ (kx,i ) = G̃(kx,i )D̃(kx,i )D̃(kx,i ) = [D̃1 (kx,i ) D̃2 (kx,i )]T

クロストークキャンセラD̃(kx,i ) = G̃(kx,i ) 1 P̃ (kx,i )F̃1,1 (kx,i )1G̃(kx,i ) = F̃ (kx,i ) =F̃1,2 (kx,i )F̃2,1 (kx,i )F̃2,2 (kx,i )

時空間周波数領域の各スピーカの駆動信号
D̃1 (kx,i ) = F̃1,1 (kx,i )P̃1 (kx,i ) + F̃2,1 (kx,i )P̃2 (kx,i )D̃2 (kx,i ) = F̃1,2 (kx,i )P̃1 (kx,i ) + F̃2,2 (kx,i )P̃2 (kx,i )13Generalized proposed method
時空間周波数領域クロストークキャンセラ(一般化)
複数平行直線アレイの時空間領域信号の行列表記P̃ (kx,i ) = [P̃1 (kx,i ) P̃2 (kx,i ) · · · P̃N (kx,i )]T
3
2
G̃(kx,i , r1,1 ) · · · G̃(kx,i , rM,1 )
7
6.
.
.
.
.
.
G̃(kx,i ) = 4
5.
.
.
G̃(kx,i , r1,N )···P̃ (kx,i ) = G̃(kx,i )D̃(kx,i )G̃(kx,i , rM,N )D̃(kx,i ) = [D̃1 (kx,i ) D̃2 (kx,i ) · · · D̃M (kx,i )]T

クロストークキャンセラD̃(kx,i ) = G̃(kx,i )+ P̃ (kx,i )2F̃1,1 (kx,i )6..
+G̃(kx,i ) = F̃ (kx,i ) = 4.

F̃1,M (kx,i )

···..
.
···

3F̃N,1 (kx,i )7..
5.

F̃N,M (kx,i )14

Block diagram of proposed method
Multiple linear microphone arrays

12Mzp(t)Original sound ﬁeld

Discrete Fourier transformxP (!)

Spatial DFTsicrorm
inea

LyarraenphoP̃ (kx , !)Crosstalk cancellers
in wavenumber domainyD̃(kx,i ) = F̃ (kx,i )P̃ (kx,i )
zuloearLinereakdsp

D̃(kx , !)ysarraSpatial IDFT

Reproduced sound ﬁeldD(!)

x

Inverse DFT

reencsig

Bd(t)
lp
ntro

Co

Tapering window

soint

y

Multiple linear loudspeaker arrays
12

N

15

Computer simulation
シミュレーション条件
解析周波数：500 Hz
音速：343.26 m/s
音源：Monopole
音源位置：(0, ­0.5, 0.5)
スピーカ間隔：0.0.7 m
スピーカ数：96 2 = 192 ch
マイクロホン数：96 x 2 =192 ch
時間平均二乗誤差
P2|p̂(t,x)
p(t,x)|iiE(x) = 10log10 i P
2|p(t,x)|ii16

Simulation results (z=0)
Z = 0.0 [m]Z = 0.0 [m]Z = 0.0 [m]66111200−2−0.5

0.5

X axis [m]0.54200−2−0.5−440.5

X axis [m]4−103

Y axis [m]3.5

4−2−4−1−62

2.5

3Y axis [m]

3.5

4−62

2.5

Z = 0.0 [m]3Y axis [m]

3.5

401

1−5−50.5Error

0.5−100−15−0.5−100−15−0.5−20−1−20−12Original−6

Z = 0.0 [m]0X axis [m]2.50−0.5−122−4X axis [m]X axis [m]62.5

3Y axis [m]3.5

LS method

4−25[dB]

22.5

3

Y axis [m]3.5

4

Proposed method−25[dB]

17

Simulation results (x=0)
X = 0.0 [m]X = 0.0 [m]X = 0.0 [m]66

111200−2−0.50.5

Z axis [m]0.54
200−2−0.5−440.5

Z axis [m]4−12.5

3Y axis [m]3.5

400−2−0.5−4−1−6−122.5

3Y axis [m]

3.5

4−62

X = 0.0 [m]2.5

3Y axis [m]

3.5

40

1

1−5−50.5Error

0.5−100−15−0.5−100−15−0.5−20−1−20−12

Original−6X = 0.0 [m]0Z axis [m]22−4Z axis [m]Z axis [m]62.5

3Y axis [m]3.5

LS method

4−25[dB]

22.5

3Y axis [m]3.5

4

Proposed method−25[dB]

18

Simulation results (y=3)
Y = 3.0 [m]Y = 3.0 [m]Y = 3.0 [m]66111200−2−0.50.5

Z axis [m]0.54200−2−0.5−440.5

Z axis [m]4−1−0.50

X axis [m]0.5100−2−0.5−4−1−6−1−1−0.50

X axis [m]

0.5

1−6−1

Y = 3.0 [m]−0.50

X axis [m]

0.5

101

1−5−50.5

Error

0.5
−100−15−0.5−100−15−0.5−20−1−20−1−1

Original−6Y = 3.0 [m]0Z axis [m]−12−4Z axis [m]Z axis [m]6−0.50X axis [m]0.5LS method

1−25[dB]−1−0.50

X axis [m]0.51

Proposed method−25[dB]

19

Conventional LS vs Proposed350

変換フィルタ行列の大きさ
従来法： M K ⇥ N K
提案法： M ⇥ N

300

条件数
従来法： 20log10 {cond(G(!))}
提案法： 20log10 {cond(F̃ (kx,i ))}→ 提案法が格段に小さい

Condition number [dB]

変換フィルタの安定性

Nyquist frequency

250
200
150
1005000LS method
Regularized LSMaximum (proposed)
Average (proposed)

12Frequency [kHz]34

変換フィルタ導出の計算量
直線アレイの数 M, N ⌧ 各アレイの素子数 K
従来法：逆行列の演算 O(K 2 )提案法：空間フーリエ変換(FFT) O(Klog2 K)→ K が大きい場合は提案法の方が演算量が小さい
20

Concluding remarks
複数平行直線アレイを用いた音空間の収録と再生
時空間領域クロストークキャンセラの提案
定式化
計算機シミュレーション
従来法よりも安定かつ高速な変換フィルタを実現

今後の予定
実環境での測定
音像定位実験

21