PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

Robust localization of early reflections in a room usingsemi real-valued EB-ESPRIT with three recurrence relations andLaplacian constraint
Byeongho JO1; Jung-Woo CHOI2
1, 2

Korea Advanced Institute of Science and Technology, Republic of Korea

ABSTRACT
To localize early reflections in a reverberant room, the direction-of-arrival (DOA) estimation technique has
been evolved in various ways. Among many techniques, the subspace-based localization techniques such as
eigenbeam estimation of signal parameters via rotational invariance techniques (EB-ESPRIT) have been
extensively researched for their ability to directly estimate DOA parameters from spherical microphonerecordings. The EB-ESPRIT estimates DOAs by computing the directional parameters from a single
recurrence relation of spherical harmonics. For estimating DOAs of early reflections, however, the
EB-ESPRIT is inappropriate due to the insufficient number of detectable echoes and its vulnerability to noise.
In this study, we attempt to overcome these issues by using semi real-valued processing of sphericalharmonics and extra constraints. The proposed technique can compute directional parameters more robustly
by deriving the real-valued forms of EB-ESPRIT equations and by imposing additional Laplacian constraint.
We show that the real-valued formula solved by the joint eigenvalue decomposition (JEVD) algorithm canenhance the accuracy of DOA estimation as compared to all existing EB-ESPRIT based techniques. In
simulations, we validate the robustness of the proposed technique in terms of the estimation error under
different signal-to-noise ratios (SNRs).
Keywords: Three Recurrence relations, Spherical microphone array, Early reflections

1. INTRODUCTION
There have been many studies to estimate the directions-of-arrival (DOAs) or localize the earlyreflections of a room using a spherical microphone array (SMA), which is applicable to the room
geometry inference (1), robust adaptive beamforming (2), de-reverberation (3), and sound field
reproduction with compensation of the reverberant components (4). Mabande et al. proposed theroom geometry inference method using the eigenbeam beamforming techniques with timedifferences of arrival (TDOAs) estimated by the cross-correlation analysis (1, 5). Morgenstern et al.
proposed to estimate DOAs of reflections using the de-correlation technique with a loudspeakerarray and a microphone array (6, 7). Tervo et al. introduced the maximum-likelihood (ML) method
based on the stochastic model to estimate the coherent early reflections (8).
These studies have limitations that the exhaustive grid-search of peaks should be accompaniedfrom the computed beamforming map. In contrast, the eigenbeam estimation of signal parameters via
rotational invariant techniques (EB-ESPRIT) has the advantage that it directly estimates DOAs in aparametric way (9). However, the limited number of detectable sound sources and the possibility ofmatrix ill-conditioning problems due to reflections impinging from the equator direction hinder theuse of EB-ESPRIT in highly reverberant conditions (10, 11). For these reasons, many alternative
types of EB-ESPRIT techniques have been proposed to increase the number of detectable sources byusing the different recurrence relations or to resolve the matrix ill-conditioning problem (12-15).
In this work, we introduce the novel EB-ESPRIT technique that combines three recurrence
relations and the semi real-valued processing. The semi real-valued processing is known to improve1
2

byeongho@kaist.ac.kr
jwoo@kaist.ac.kr

4949

the estimation performance albeit with the reduced computational complexity (16-19). In addition,
we estimate the directional parameter matrices by using iterative joint eigenvalue decomposition(JEVD) method (20-22) with the additional constraint based on the Helmholtz equation. By using theadditional constraint, the diagonalization of JVED is restricted by a physical constraint, which
reduces the effect of the noises. The estimation performance of the proposed method is validated by
the simulated spherical microphone array data installed in a box-shaped room under various noise
conditions.

2. SPHERICAL HARMONIC DOMAIN MODEL
2.1

Complex-valued data model

The 2D fields f (W) defined on a unit sphere ( S 2 ) can be converted into the spherical harmonic
domain (SHD) through the spherical Fourier transform (SFT) as (23)f nm := ò f (W)Ynm (W)* d W ,

(1)

S2

where W Î {(q , f ) | q Î [0, p ], f Î [0, 2p ]} denotes 2D angles in the spherical coordinates and (×)*
denotes the complex conjugate. Spherical harmonics is defined as Ynm (W) = N nm Pnm (cos q )eimf ,where N nm = {(2n + 1)(n - m)!/ 4p (n + m)!}1/ 2 and Pnm (×) is the associated Legendre function. For a unit
amplitude plane-wave with DOA W PW , the SHD signal is transformed from the spherical microphonearray (SMA) recording data as
f PW , nm = bn (kr )Ynm (W PW )* ,(2)
where k denotes the wavenumber, r is the radius of the spherical microphone array and bn (kr ) is amode-strength that depends on the radius and the boundary condition of the array surface. To get rid of
this dependency, the SHD signal is divided by the bn (kr ) of the respective orders ( n = 1,L , N ). Aftercompensation, the mode-compensated SHD (MC-SHD) signal of a unit amplitude plane-wave can be
written as aPW , nm = Ynm (W PW )* (24).
For D plane waves, the MC-SHD signal up to the truncated harmonics order N can be written as
the vector form:(3)
a = YH s ,where a = [a0,0 , a1, -1 , a1,0 , a1,1 , L , aN , N ]T , s = [ s1 ,L , sd , L , sD ]T , sd is the complex amplitude of d -th
Y H is the steering matrix in SHD whose d -th column is given byy (W d ) = [Y00 (W d ), Y1-1 (W d ), Y10 (W d ), Y11 (W d ), L , YNN (W d )]H . (×)T and (×) H denote the transpose andsource, andconjugate transpose, respectively. We assume that all MC-SHD signals are made by the short-time
Fourier transform (STFT), SFT, and mode-compensation in order. To simplify the notations, timeframe and frequency indices of MC-SHD signals are omitted. For EB-ESPRIT type algorithms, thecovariance matrix should be constructed for extracting the signal subspaces as follows,(4)
R c = E{aa H } = Y H E{ss H }Y = U s Σ s U s H ,where U s and Σ s are the signal subspace and the corresponding eigenvalue matrix, respectively.
E{}× denotes the expectation operator. The subscript c of R c represents the complex-valued matrix,
which is to distinguish the complex-valued covariance matrix from the real-valued covariance matrix
introduced in the following section.
2.2 Real-valued data model
For the unitary or semi real-valued processing (12, 18), the real-valued covariance matrix should beconstructed. First, the real spherical harmonics is defined as follows:ì 2 N n m Pn m (cos q ) sin( m f ), m < 0ïïYˆnm ( W ) = í N nm Pnm (cos q ), m=0.
ïm
ïî 2 N nm Pn (cos q ) cos(mf ) , m > 0

4950

(5)

Using the complex symmetry property of the spherical harmonics and Equation (5), the complex
spherical harmonic matrix can be transformed by a following (2n + 1) ´ (2n + 1) unitary matrix for ordern ³ 1:
é iI n
1 ê T
Qn =
ê 0n
2ê
êëG n

0n
2
0n

(-1) 2 iJ n ùú
0Tn ú ,J n G n úúû(6)where Q 0 = 1 , and G n is a n ´ n anti-diagonal matrix with ones in anti-diagonal positions and zeros
elsewhere, 0n represents a column vector containing of n zeros, and J n is a n ´ n anti-diagonal
matrix whose (i, n + 1 - i ) -th element is set to (-1)i . With Q n matrices ( n = 0,L , N ), one can constructa unitary matrix Q = blkdiag{Q 0 , L , Q N } . As shown in (12), a real-valued steering matrix can beobtained as QY H = Yˆ T = [yˆ 1 , L , yˆ D ] , where yˆ = [Yˆ00 (q , f ), Yˆ1-1 (q , f ), Yˆ10 (q , f ),L , YˆNN (q , f )]T .
Based on the proof in (12), the real-valued covariance matrix can be computed from thecomplex-valued covariance matrix as
ˆ T E{ss
ˆ =U
ˆ Σ
ˆ ˆH
ˆˆ H }Y(7)
R r = Re{QR c Q H } = Y
s s Uswhere sˆ = [Re{s}, Im{s}] , Uˆ s is the real-valued signal subspace, Σˆ s is the real-valued eigenvalue
matrix, and Re{}× extracts the real part of the argument. In the following section, the real-valued
signal subspace in Equation (7) is utilized for the recurrence relations of complex spherical harmonics.

3. SEMI REAL-VALUED EB-ESPRITThe proposed semi real-valued (SR)-EB-ESPRIT utilizes three recurrence relations of sphericalharmonics whose directional parameters correspond to x, y , and z directions in the Cartesian
coordinates (25, 26):
sin q d cos fd Ynm (W d )* = ( wnm+-11Ynm+1-1 (W d )* - wn- mYnm-1-1 (W d )* - wn-+( 1m +1)Ynm+1+1 (W d )* + wnmYnm-1+1 (W d )* ) / 2
sin q d sin fd Ynm (W d )* = ( wnm+-11Ynm+1-1 (W d )* - wn- mYnm-1-1 (W d )* + wn-+( 1m +1)Ynm+1+1 (W d )* - wnmYnm-1+1 (W d )* ) / 2i ,
m
d n

m m
n n -1*mm
n +1 n +1

*

(8)

*

cos q Y (W d ) = v Y (W d ) + v Y (W d )mnwhere w = (n - m - 1)(n - m) / (2n - 1)(2n + 1) and vnm = (n - m)(n + m) / (2n - 1)(2n + 1) .
The above three recurrence relations can be re-written as the matrix form with the shifted andorder-reduced ( N - 1 ) harmonics matrix Y ( m ,n ) as
Y (0,0) H Φ (1) = ( Wnm+-11Y ( +1, -1) H - Wn- m Y ( -1, -1) H - Wn-+(1m +1) Y ( +1, +1) H + Wnm Y ( -1, +1) H ) / 2
Y (0,0) H Φ (2) = ( Wnm+-11Y ( +1, -1) H - Wn- m Y ( -1, -1) H + Wn-+(1m +1) Y ( +1, +1) H - Wnm Y ( -1, +1) H ) / 2i ,Y(1)(0,0) H
(2)

Φ

(3)

m
n

=V Y

( -1,0) H

m
n +1

+V Y

(9)

( +1,0) H(3)where Φ , Φ , and Φ
are D ´ D real-valued diagonal matrices containing directionalparameters sin q d cos fd , sin q d sin fd , and cos q d , respectively. Wnm++mn and Vnm++mn are N 2 ´ N 2
real-valued diagonal matrices containing wnm++mn and vnm++mn , respectively. The d -th column of Y ( m ,n ) H
is defined as ynm (W d ) = [Y00++mn (W d ), Y1-+1m+n (W d ), Y10++mn (W d ), Y11++mn (W d ), L , Y((NN--1)1)++mn (W d )]H .
As the directional parameter matrices ( Φ(1) , Φ(2) , and Φ(3) ) contain DOAs ( W d ), one has to solve
the Equation (9) for the matrices ( Φ(1) , Φ(2) , and Φ(3) ). However, Y ( m ,n ) matrices cannot becomputed directly from the MC-SHD signals. We utilize a relationship between the real-valued signal
subspace ( Uˆ s ) and the complex-valued steering matrix ( Y H ) in SHD (12), given byˆ = Y H Tˆ ,Us = QH U
s(10)where T̂ is a D ´ D real-valued nonsingular matrix. Using Equation (10), we can transform areal-valued signal subspace ( Uˆ s ) obtained from a real-valued covariance matrix ( R r ) into acomplex-valued signal subspace ( U s ) that can be applied to recurrence relations of complex sphericalharmonics (Equation (8)), and obtain a real-valued common transformation matrix T̂ as well. From
Equations (7) and (10), and by multiplying T̂ matrix, Equation (9) can be re-written as4951

U (0,0)Ψ (1) = ( Wnm+-11U (s+1, -1) - Wn- m U (s-1, -1) - Wn-+(1m +1) U (s+1, +1) + Wnm U (s-1, +1) ) / 2sU (0,0)Ψ (2) = ( Wnm+-11U (s+1, -1) - Wn- m U (s-1, -1) + Wn-+(1m +1) U (s+1, +1) - Wnm U (s-1, +1) ) / 2i ,
s
U

(0,0)
s

Ψ

(3)

m
n

=V U

( -1,0)
s

m
n +1

+V U

(11)

( +1,0)s

where Ψ ( I ) = Tˆ -1Φ( I ) Tˆ ("I Î {1, 2,3}) . The least-squares solution of Equation (11) yield Ψ ( I ) matrices.
The least-squares solution can be unique when the pseudo-inverse matrix U (0,0)
is the overdeterminedssystem, which corresponds N 2 ³ D . From this inequality condition, the maximum number ofdetectable sources is determined to N 2 .
To estimate simultaneously Φ( I ) matrices from the estimated Ψ ( I ) matrices, one should computejoint-diagonalizing matrix. In the following section, we show how to compute the joint-diagonalizing
matrix A by using a block coordinate procedure and additional constraints. From the estimated A
matrix that corresponds to the inverse of the real-valued common transformation matrix ( Τˆ -1 ), theˆ ( I ) matrices can be computed from Φˆ ( I ) = AΨ ( I ) A -1 . Using the elements of the Φˆ ( I ) matrices, DOA
Φof d -th source ( W d ) can be estimated with a four-quadrant inverse tangent ( atan{}× ) asˆ (1) ] }}2 + {Re{[Φˆ (2) ] }}2 , Re{[Φˆ (3) ] }}qˆd = atan2{ {Re{[Φ
d ,d
d ,d
d ,d
ˆ (1) ] }}
ˆ (2) ] }, Re{[Φfˆd = atan2{Re{[Φ
d ,d
d ,d

,

(12)

4. JOINT DIAGONALIZATION WITH CONSTRAINT
4.1 A block coordinate procedure
As mentioned in (15), Ψ ( I ) matrices are jointly diagonalized by using the joint eigenvalue
decomposition algorithms, which can be achieved by minimizing
3

J ( A) = å ZDiag {AΨ ( I ) A -1}
I =1

2
F

,

(13)

where ZDiag{}× sets to zero the diagonal of the matrix element, and × F denotes the Frobenius normof the matrix. In this work, we utilize a block coordinate procedure of the JAPAM algorithm (Joint
eigenvalue decomposition algorithms using a Parameterized Matrix) (22). The block coordinate
procedure divides the problem of estimating A by minimizing J ( A) into D( D - 1) / 2 successive
simpler optimization problems by using the following matrices,D -1

DA = Õ i =1 Õ j = i +1 A ( i , j ) ,(14)

where matrices A (i , j ) are equal to the identity matrix that depends on four unknown parameters:
A1(i , j ) , A2(i , j ) , A3(i , j ) , A4(i , j ) :ì"(i ¢, j ¢) Ï {(i, i ), (i, j ), ( j , i ), ( j , j )}, [ A (i , j ) ]i ¢, j ¢ = d i ¢, j ¢ï (i , j )(i , j )(i , j )(i , j ),í[ A ]i ,i = A1 , [ A ]i , j = A2
ï[ A (i , j ) ] = A(i , j ) ,[ A (i , j ) ] = A(i , j )
j ,i
3
j, j
4
î

(15)

where d is the Kronecker delta function. For (i, j ) step, the following matrices are updated asìï A ¬ A (i , j ) A.
í (I )(i , j )(I )( i , j ) -1
ïîΨ ¬ A Ψ ( A ) , "I Î {1, 2,3}

(16)

In (i, j ) -update, we define the following matricesΘ ( I ) = A (i , j ) Ψ ( I ) ( A (i , j ) ) -1 , "I Î {1, 2,3} .

(17)Then, the cost function (Equation (13)) can be re-written as
3

C ( A (i , j ) ) = å I =1 (| [Θ ( I ) ]i , j |2 + | [Θ ( I ) ] j ,i |2 ) .

4952

(18)

As we focus on Θi(,I j) and Θ (jI,i) matrices, Equation (18) can be transformed to the simpler costfunction form with 2 ´ 2 matrices as% (i , j ) ) = 3 (| [ΘC (Aå I =1 % ( I ) ]1,2 |2 + | [Θ% ( I ) ]2,1 |2 ) ,

(19)

where
% ( I ),(i , j ) = A% (i , j ) Ψ% ( I ),(i , j ) ( A% (i , j ) ) -1 =
Θ
-1

.
(20)
[Ψ ( I ) ]i , j ù é A1(i , j ) A2(i , j ) ù, "I Î {1, 2,3}ú ê (i , j )(I )(i , j ) ú[Ψ ] j , j ûú ë A3
A4 û
% ( I ) ] and [Θ% ( I ) ] are computed as
In this work, without assumptions in (22), two elements [Θé [Θ ( I ) ]i ,iê (I )
ëê[Θ ] j ,i[Θ ( I ) ]i , j ù é A1(i , j )ú=ê[Θ ( I ) ] j , j ûú ë A3(i , j )(I )A2(i , j ) ù é [Ψ ]i ,iú ê (I )A4(i , j ) û ëê[Ψ ] j ,i

1,2

2,1

% ( I ) ] = A(i , j ) A(i , j ) ( Y ( I ) - Y ( I ) ) + ( A(i , j ) ) 2 Y ( I ) - ( A(i , j ) ) 2 Y ( I )
[Θ
1,2
1
2
j. j
i ,i
1
i, j
2
j ,i
% ( I ) ] = A(i , j ) A(i , j ) ( Y ( I ) - Y ( I ) ) + ( A(i , j ) ) 2 Y ( I ) - ( A(i , j ) ) 2 Y ( I )
[Θ
1,2
3
4
i, j
j, j
4
j ,i
3
i, j

.

(21)

4.2 JEVD with Laplacian constraint
The sound field satisfies the Helmholtz equation (27). Therefore, the following equation always
holds for the wavenumber components ( k x , k y , k z ) in the Cartesian coordinates:
k x2 + k y2 + k z2 = k 2 .
(1)

(2)

(22)

(3)

In an ideal case without noises, Φ , Φ , and Φ matrices obtained from the acquired soundfield are the real-valued diagonal matrices with ( cos fd sin q d , sin fd sin q d , and cos q d ) as their(d , d ) -th element, respectively. Therefore, the following expression should be satisfied:Φ (1) Φ (1) + Φ (2) Φ (2) + Φ (3) Φ (3)ˆ =Iˆ -1 (Ψ (1) Ψ (1) + Ψ (2) Ψ (2) + Ψ (3) Ψ (3) )T
=T

.

(23)

D´ D

However, in practical situations, it is impossible to satisfy Equation (23) because of the noises. In
this work, we propose to use a robust JEVD algorithm by adding the iterative JEVD algorithm (22) into
the above Laplacian constraint (L-constraint). The cost function of the proposed JEVD algorithm ischanged from the cost function of a block coordinate procedure (Equation (19)) as% (i , j ) ) = (1 - l )Cconst ( Aå I =1 ( [Θ% ( I ) ]1,2 + [Θ% ( I ) ]2,1 ) + l ( [Ξ% ]1,2 + [Ξ% ]2,1 + [Ξ% ]1,1 - 1 + [Ξ% ]2,2 - 1 ) ,
2

3

2

2

2

2

2

(24)

% (i , j ) (Ψ% (1)(i , j ) Ψ% (1)(i , j ) + Ψ% (2)(i , j ) Ψ% (2)(i , j ) + Ψ% (3)(i , j ) Ψ% (3)(i , j ) )( A
% (i , j ) ) -1 Î ¡ 2´2 , and elements of Ξ
% (i , j )where Ξ% (i , j ) = A
are computed as follows,% ] = A(i , j ) A(i , j ) (Y (i , j ) - Y (i , j ) ) + ( A(i , j ) ) 2 Y (i , j ) - ( A(i , j ) ) 2 Y (i , j )
[Ξ
1,2
1
2
(4)
(1)
1
(2)
2
(3)
% ] = A(i , j ) A(i , j ) (Y (i , j ) - Y (i , j ) ) + ( A(i , j ) ) 2 Y (i , j ) - ( A(i , j ) ) 2 Y (i , j )
[Ξ
2,1
3
4
(1)
(4)
4
(3)
3
(2)
% ] = A(i , j ) A(i , j ) Y (i , j ) - A(i , j ) A(i , j ) Y (i , j ) + A(i , j ) A(i , j ) Y (i , j ) - A(i , j ) A(i , j ) Y (i , j )
[Ξ
1,1
1
4
(1)
1
3
(2)
2
4
(3)
2
3
(4)

,

(25)

% ] = A(i , j ) A(i , j ) Y (i , j ) + A(i , j ) A(i , j ) Y (i , j ) - A(i , j ) A(i , j ) Y (i , j ) - A(i , j ) A(i , j ) Y (i , j )
[Ξ
1,2
1
4
(4)
1
3
(2)
2
4
(3)
2
3
(1)

where
3

Y ((1)i , j ) = å I =1 ([Ψ ( I ) ]i ,i [Ψ ( I ) ] i ,i + [Ψ ( I ) ]i , j [Ψ ( I ) ] j ,i )3

i, j )Y ((2)= å I =1 ([Ψ ( I ) ]i ,i [Ψ ( I ) ]i , j + [Ψ ( I ) ]i , j [Ψ ( I ) ] j , j )3

Y ((3)i , j ) = å I =1 ([Ψ ( I ) ] j ,i [Ψ ( I ) ]i ,i + [Ψ ( I ) ] j , j [Ψ ( I ) ] j ,i )

.

(26)

3

i, j )Y ((4)= å I =1 ([Ψ ( I ) ] j ,i [Ψ ( I ) ]i , j + [Ψ ( I ) ] j , j [Ψ ( I ) ] j , j )In Equation (24), l is the parameter that determines the weight between the original cost function(Equation (19)) and the L-constraint, ranging from 0 to 1. In this work, for each update of (i, j ) , we4953

solve the unconstraint optimization problem with the cost function in Equation (24) by using the
Quasi-Newton method (28, 29).

5. EVALUATION
5.1 Setup
To investigate the estimation performance for the early reflections, we conducted simulations withsynthesized room impulse responses. The impulse responses of the room were generated with aspherical microphone array impulse response generator (30) and convolved with the white Gaussian
noise to generate MC-SHD data. The rigid body SMA with radius 0.042 m that composed of 32
uniformly sampled microphones was used. The SMA signals were sampled at 16 kHz. The STFT wasdone by Rectangular window of 480 samples (30 ms) in length, a hop size of 120 samples (7.5 ms), andDFT size of 512.
The room setup is described in the left side of Figure 1. The height of the room is 3.03 m, and theheight of the center of the SMA and loudspeaker unit are 1.41 m. The sound field was decomposed up
to the fourth-order of spherical harmonics ( N = 4). To de-correlate the highly correlated earlyreflections, the real-valued covariance matrix ( R r ) was averaged over 500 time frames and frequency
range within 2.275–3.9 kHz ( kr = 1.75 ~ 3 ) (31, 32). The signal-to-noise ratio (SNR) is defined to thetotal power of the sound sources to the total power of the self-microphone noises. The self-microphone
noises were generated with the uncorrelated white Gaussian noises. The ground truth DOAs of direct
sound and early reflections are shown in the right side of Figure 1. The estimation errors between thetrue DOAs and estimated DOAs were computed as DW d = cos -1 (cos q d cos qˆd + cos(fd - fˆd ) sin q d sin qˆd ) .
K

D

2The root-mean-square error (RMSE) is defined as RMSE(W) = {(1/ KD)å k =1 å d =1 DW d }1/ 2 with
K = 300 independent trials. The reverberation time of the simulated room was set to 300 ms.

Figure 1 – The room setup and ground truth DOAs of the direct sound and early reflections in a room.
The baselines of this work were DOA-vector EB-ESPRIT (15) and the SR-EB-ESPRIT using
JAPAM-5 algorithm without L-constraint (22). The diagonalizing matrix A in the first iteration was
initialized as the inverse of eigenvector matrix of Ψ (1) for JAPAM-5 algorithm and the proposed
optimization algorithm. The maximum iteration number was set to 150, and the stopping criterion for
JAPAM-5 was set as same in (22). In the proposed method, l was set to 0.05. The optimization withthe proposed cost function (Equation (24)) was solved by using the unconstraint multivariable andnonlinear optimization solver (fminunc) of the MATLAB® . The optimization algorithm was selectedas ‘Quasi-Newton’ algorithm (29), and the initial points of A1(i , j ) , A2(i , j ) , A3(i , j ) , and A4(i , j ) were selected% (i , j ) ) ) for each (i, j ) step. The total
from {0.3, 0.65, 1} that minimize the cost function output ( Cconst ( A
iteration number of the proposed algorithm was set to 1.

5.2 Results
In Figure 2 (a), RMSEs evaluated in various SNR conditions (SNR = 10 dB ~ 40 dB) are shown.
RMSEs for all algorithms change slightly as SNR decreases. This is because the artifacts of the late4954

reverberant components have dominant influences on the performance of algorithms. The RMSEs of
SR-EB-ESPRIT with JAPAM-5 and with L-constraint are only one-fourth compared to the case of
DOA-vector EB-ESPRIT in all SNR regions. The proposed SR-EB-ESPRIT with L-constraint
estimates slightly better than the SR-EB-ESPRIT with JAPAM-5, which is caused by adding thephysical constraints based on the Helmholtz equation. To validate the performance for each reflection,
the RMSE for each wall was examined in SNR = 10 dB condition (Figure 2 (b)). Without direct sound,
RMSEs of DOA-vector EB-ESPRIT are quite high. Unlike the RMSEs of DOA-vector EB-ESPRIT,
both RMSEs of SR-EB-ESPRIT with JAPAM-5 and SR-EB-ESPRIT with L-constraint are reasonably
small (under 4 degrees), but their performance gap is ignorable. To summarize, both SR-EB-ESPRIT
techniques proposed in this work show outstanding performance compared to the state-of-the art
technique. The L-constraint brings a slight performance gain in overall RMSEs, but the differences
depend on the wall directions.

(a)

(b)

Figure 2 – (a) RMSEs for various SNRs (10 ~ 40 dB)(b) RMSEs of the direct sound and the early reflections with SNR = 10 dB.

6. CONCLUSIONSWe proposed the DOA estimation method for direct and early reflection sounds in a room by using
the semi real-valued processing and joint eigenvalue decomposition with constraint. The proposed
method can reduce the computational complexity by using the real-valued covariance matrix andreal-valued common transformation matrix. In addition, by using the L-constraint based on theHelmholtz equation, the robustness of the proposed method is slightly improved. Through thesimulation results in various SNRs, we showed that the estimation performance of the proposed
method outperforms other baseline algorithms in a room.

ACKNOWLEDGEMENTS
This work was supported by the Korea Research Foundation (KRF) under Contract No.
2019R1A2C1007393 and by the BK21 Plus program through the National Research Foundation (NRF)
funded by the Ministry of Education of Korea.

REFERENCES
1. Mabande E, Kowalczyk K, Sun H, Kellermann W. Room geometry inference based on sphericalmicrophone array eigenbeam processing. J Acoust Soc Am. 2013; 134(4) 2773–2789.
2. Sun H, Yan S, Svensson U. P. Space domain optimal beamforming for spherical microphone arrays.
Proc IEEE Int Conf Acoust Speech Signal Process; 2010. p. 117–120.
3. Peled Y, Rafaely B. Method for dereverberation and noise reduction using spherical microphone arrays.
Proc IEEE Int Conf Acoust Speech Signal Process; 2010. p. 113–116.
4. Betlehem T, Poletti M. A. Two dimensional sound field reproduction using higher order sources to
exploit room reflections. J Acoust Soc Am. 2014; 135(4) 1820–1833.

4955

5. Sun H, Mabande E, Kowalczyk K, Kellermann W. Joint DOA and TDOA estimation for 3D
localization of reflective surfaces using eigenbeam MVDR and spherical microphone arrays. Proc
IEEE Int Conf Acoust Speech Signal Process. 2011; p. 113–116.
6. Morgenstern H, Zotter F, Rafaely B. Joint spherical beam forming for directional analysis of reflections
in rooms. J. Acoust. Soc. Am. 2012; 131(4) 3207.
7. Morgenstern H, Rafaely B. Modal smoothing for analysis of room reflections measured with sphericalmicrophone and loudspeaker arrays. J Acoust Soc Am. 2018; 143(2) 1008–1018.
8. Tervo S, Politis A. Direction of arrival estimation of reflections from room impulse responses using a
spherical microphone array. IEEE/ACM Trans Audio Speech Lang Process. 2015; 23(10) 1539–1551.
9. Teutsch H, Kellermann W. Eigen-beam processing for direction-of-arrival estimation using sphericalapertures. Proc. Joint Workshop Hands-Free Commun and Microphone Arrays; 2005.
10. Sun H, Teutsch H, Mabande E, Kellermann W. Robust localization of multiple sources in reverberant
environments using EB-ESPRIT with spherical microphone arrays. Proc IEEE Int Conf Acoust Speech
Signal Process; 2011. p. 117–120.
11. Sun, H, Mabande E, Kowalczyk K, Kellermann W. Localization of distinct reflections in rooms using
spherical microphone array eigenbeam processing. J Acoust Soc Am. 2012; 131(4) 2828–2840.
12. Huang Q, Zhang L, Yong F. Two Step Spherical Harmonics ESPRIT-Type Algorithms and Performance
Analysis. IEEE/ACM Trans Audio Speech Lang Process. 2018; 26(9) 1684–1697.
13. Jo B, Choi J. W. Direction of arrival estimation using nonsingular spherical ESPRIT. J Acoust Soc Am.
2018; 143(3) EL181–EL187.
14. Jo B, Choi J. W. Parametric direction-of-arrival estimation with three recurrence relations of sphericalharmonics. J Acoust Soc Am. 2019; 145(1) 480–488.
15. Herzog A, Habets E. A. Eigenbeam-ESPRIT for DOA-Vector Estimation. IEEE Signal Process. Lett.
2019; 26(4) 572–576.
16. Huarng K. C, Yeh C. C. A unitary transformation method for angle-of-arrival estimation. IEEE Trans
Signal Process. 1991; 39(4) 975–977.
17. Haardt M, Nossek J. A. Unitary ESPRIT: How to obtain increased estimation accuracy with a reduced
computational burden. IEEE Trans signal process. 1995; 43(5) 1232–1242.
18. Goossens R, Rogier H. Unitary spherical ESPRIT: 2D-angle estimation with spherical arrays for scalar
fields. IET signal process. 2009; 3(3) 221–231.
19. Krishnaveni V, Kesavamurthy T, Aparna B. Beamforming for direction-of-arrival (DOA) estimation-asurvey. Int J Comput Appl. 2013; 61(11).
20. Fu T, Gao X. Simultaneous diagonalization with similarity transformation for non-defective matrices.
Proc IEEE Int Conf Acoust Speech Signal Process; 2006. p. IV.
21. Balda E. R, Cheema S. A, Weiss A, Yeredor A, Haardt M. Perturbation analysis of joint eigenvalue
decomposition algorithms. Proc IEEE Int Conf Acoust Speech Signal Process; 2017. p. 3101–3105.
22. André R, Luciani X, Moreau E. A new class of block coordinate algorithms for the joint eigenvalue
decomposition of complex matrices. Signal Process. 2018; 145 78–90.
23. Rafaely B. Fundamentals of spherical array processing; 2015.
24. Jarrett D, Habets E, Naylor P. Theory and applications of spherical microphone array processing; 2017.
25. Choi C, Ivanic J, Gordon M. S, Ruedenberg K. Rapid and stable determination of rotation matrices
between spherical harmonics by direct recursion. J Chem Phys. 1999; 111(19) 8825–8831.
26. Gumerov N. A, Duraiswami R. Fast multipole methods for the Helmholtz equation in three
dimensions; 2005.
27. Kinsler L. E, Frey A. R, Coppens A. B, Sanders J. V. Fundamentals of Acoustics. 4th Edition; 1999.
28. Chong E. K, Zak S. H. An introduction to optimization; 2013.
29. Shanno D. F. Conditioning of quasi-Newton methods for function minimization. Math of comput.
1970; 24(111) 647–656.
30. Jarrett D, Habets E. A. P, Thomas M. R. P, Naylor P. A. Rigid sphere room impulse response
simulation: Algorithm and applications. J Acoust Soc Am. 2012; 132(3) 1462–1472.
31. Khaykin D, Rafaely B. Coherent signals direction-of-arrival estimation using a spherical microphonearray: Frequency smoothing approach. Proc IEEE Workshop on Appl of Signal Process Audio andAcoustics; 2009. p. 221–224.
32. Huleihel N, Rafaely B. Spherical array processing for acoustic analysis using room impulse responses
and time-domain smoothing. J Acoust Soc Am. 2013; 133(6) 3995–4007.

4956