SPARSE SOUND FIELD DECOMPOSITION
WITH MULTICHANNEL EXTENSION OF COMPLEX NMF
1Naoki Murata, 1 Shoichi Koyama, 1,2 Hirokazu Kameoka, 1 Norihiro Takamune, and 1 Hiroshi Saruwatari1

Graduate School of Information Science and Technology, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
2
NTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation,
3-1 Morinosato Wakamiya, Atsugi-shi, Kanagawa 243-0198, Japan
ABSTRACTA sparse sound field decomposition method using prior information
on source signals in the time-frequency domain is proposed. Sparsesound field decomposition has been proved to be effective for various
acoustic signal processing applications. Current methods for sparsedecomposition are based only on the spatial sparsity of the source
distribution. However, it can be assumed that possible source signals
to be decomposed are approximately known in advance. To exploit
this prior information, we incorporated the complex nonnegative factorization model into sparse sound field decomposition. Since themagnitude spectrum of the possible source signals can be trainedin advance, accuracy of the sparse decomposition can be improved
even when the source signals are highly correlated and the sources
are in a highly noisy environment. In addition, the proposed decomposition algorithm is derived using the auxiliary function method.
Numerical experiments indicated that the sparse decomposition performance was significantly improved using the proposed method.
Index Terms— Sound field decomposition, source localization,
sparse representation, complex nonnegative matrix factorization,
auxiliary function method
1. INTRODUCTION
Sound field decomposition forms the foundation of various acoustic
signal processing applications, such as sound field analysis, reconstruction, and visualization. Its aim is to represent a sound field asa linear combination of fundamental solutions of the wave equation(or Helmholtz equation) from the results of the pressure measurements. Using this representation, the entire sound field can be estimated from the signals received by multiple microphones. It can beconsidered that this problem involves source localization or direction of arrival (DOA) estimation because the power distribution of
the pressure field indicates the locations or directions of the soundsources.
A conventional method for sound field decomposition is basedon spatial Fourier analysis of the sound field, which corresponds toplane wave decomposition [1]. In acoustic holography, the method
based on spatial Fourier analysis is known as near-field acoustic
holography (NAH). NAH is used to measure the pressure or velocity
distribution on a surface close to acoustic sources. In the context ofsound field recording and reproduction, which is targeted at highfidelity audio systems, a method of converting signals received by
microphones into driving signals of loudspeakers for reproduction isnecessary. The wave field reconstruction filtering method is an efficient and stable signal conversion method based on spatial Fourier

978-1-4799-9988-0/16/$31.00 ©2016 IEEE

345

analysis [2, 3]. Higher-order Ambisonics [4, 5] also enables signal
conversion by encoding and decoding processes in the spherical harmonic domain.
Source localization and DOA estimation have been developedin a different context from sound field decomposition [6, 7]. Delayand sum (DAS) is one of the most classical approaches [6]. Subspace methods, such as the multiple signal classification (MUSIC)
method [8], are also well-established methods based on the assumption of orthogonality of signal and noise subspaces. Since signal
reconstruction is not necessary in source localization, a wide rangeof methods have been proposed compared with for sound field decomposition.
In recent years, the sparse representation of a sound field has
attracted interest in various research fields involving acoustic signal processing. Sparse sound field decomposition has been applied
to acoustic holography [9] and sound field recording and reproduction [10, 11] to improve their spatial resolution limits. Several methods for source localization based on sparse representation have been
proposed [12–14], which also make higher spatial resolution possible. Generally, these methods are formulated in a common framework of an overcomplete linear equation of the observations and are
based on an assumption of a spatially sparse distribution of soundsources. Therefore, improving the performance of sparse sound field
representation has important implications for various applications.
In current methods for sparse sound field decomposition, no assumptions are imposed on the structure of the source signals in thetime-frequency domain. However, in some applications, it can beassumed that the source signals of interest are known in advance.
For instance, speech signals are the main target of sparse sound fielddecomposition in telecommunication systems. This assumption is
widely used in the blind source separation (BSS) problem. In nonnegative matrix factorization (NMF) [15], which is one of the wellknown methods for BSS, a time-frequency spectrogram of the source
signal is approximated as a product of two nonnegative low-rank matrices. This means that acoustic source signals are represented by alimited number of spectrogram components.
We propose a sparse sound field decomposition method thattakes the time-frequency-spectrum structures of the source signals
into consideration. The proposed method is derived by incorporatingcomplex NMF (CNMF) [16] into sparse sound field decomposition.
The source signal is modeled as a sum of products of a rank 1 magnitude spectrogram and its phase spectrogram. Therefore, a complex
spectrogram of the source signals can be incorporated into sparsedecomposition without transformation into nonnegative values. The
decomposition algorithm is obtained by an auxiliary function approach. By exploiting time-frequency-spectrum structures trained in

ICASSP 2016

advance, it is possible to improve the decomposition accuracy evenwhen the source signals are highly correlated and the sources are ina highly noisy environment. Numerical simulations are conducted
to evaluate the proposed method in terms of DOA estimation andsignal reconstruction.
2. SPARSE SOUND FIELD DECOMPOSITIONFirst, we formulate a signal model for sparse sound field decomposition and its optimization criterion for current methods in a generalized form [9–14]. Figure 1 shows this signal model in the special
case of a linear microphone array and plane wave decomposition.
The number of sources is assumed to be small in the region of interest. We consider M microphones receiving signals from I sources.
The locations (or directions) of the sources are denoted by position
vectors in an arbitrary number of dimensions vi (i ∈ {1, · · · , I}). A
model of observed signals in the time-frequency domain at a single
time frame can be represented asy(ω) = A(ω)s(ω) + z(ω),(1)

where ω denotes the frequency, y(ω) ∈ C and s(ω) ∈ C respectively denote the signals received by the microphones and thecomplex amplitudes of the source signals, and z(ω) ∈ CM is an unknown noise vector. Hereafter, we omit ω for notational simplicity.
The quantity
A = [d(v1 ) d(v2 ) · · · d(vI )]
(2)
M

I

is referred to as a steering matrix, where each column d(vi ) ∈ CM
is a vector whose elements consist of a fundamental solution of theHelmholtz equation defined by the locations of the sources and microphones, such as the Green’s function and plane wave function. In
the source localization problem, arbitrary transfer functions can beused for A. Our objective is to estimate both the parameters vi ands. Note that only vi is necessary in the source localization.
To solve (1) as a sparse representation problem, we generalize
the steering matrix A as an overcomplete dictionary matrix D comprising all possible source locations v̂n (n ∈ {1, · · · , N }), which
are referred to as grids, such thatD = [d(v̂1 ) d(v̂2 ) · · · d(v̂N )],(3)

where N denotes the number of grids and both N ≫ I and N ≫ M
are assumed. The source vector s can be accordingly modified tox ∈ CN , where the nth element xn has a nonzero value and equals
si when the ith source originates from v̂n , and is zero otherwise.
Using D and x, (1) can be reformulated as
y = Dx + z.

(4)

Fig. 1. Signal model for plane wave decomposition using linear array of microphones.

where 0 < p ≤ 1. When observations of multiple time frames areavailable and the source locations can be assumed to be static during
the observations, a matrix in which each column is an observation isconstructed, and a row-sparse penalty term of the matrix is usuallyused for decomposition, such as the ℓp/q -norm (q ≥ 2) [11, 18, 19].
This problem is generally referred to as the multiple measurement
vector (MMV) problem.
3. PROPOSED SIGNAL MODEL AND OBJECTIVE
FUNCTION
In (4) and (5), no assumptions are imposed on x for the temporal structure; therefore, the signal decomposition is based only on
the spatial sparsity of the sources. Although it is assumed that x
for multiple time frames has the same sparsity pattern in the MMV
problem, a time-frequency-spectrum structure that is unique in wideband acoustic signals is not exploited.
Kameoka et al. proposed CNMF [16], in which a complex spectrogram in the time-frequency domain is represented as a sum of
products of a rank 1 magnitude spectrogram and its phase spectrogram. We consider now to incorporate the signal model of CNMF
into sparse sound field decomposition since the complex spectrogram of the observed signals is necessary to model signals received
by multiple microphones. In CNMF, a complex spectrogram in thetime-frequency domain is modeled as∑
fk,t ≈
hk,l ul,t ejϕk,l,t ,(7)lWhen the source distribution is spatially sparse, x becomes a sparse
vector, i.e., x has a few nonzero elements. Therefore, sparse decomposition algorithms can be applied to estimate x from y [17].
The optimization criterion for the sparse sound field decomposition can be formulated as{}minimize ∥y − Dx∥22 + ηJ(x) ,(5)xwhere J(x) is a penalty term for inducing the sparsity of x and η
is a parameter that balances the approximation error and the penalty
J(x) [17]. In general, the ℓp -norm of x is used as J(x),J (x) = ∥x∥p =
p

( N
∑

)1/p
|xn |

p

,

(6)

where k, t, and l respectively denote the frequency, time frame, andbasis index, fk,t is the complex spectrogram of a source signal, hk,lis a static magnitude spectrum, ul,t is a time-varying activation coefficient, and ϕn,k,t is a time-varying phase spectrum. This modelallows any observed complex spectrum to be accurately represented
only using a few active magnitude spectrum bases, each of which ispaired with an arbitrary phase spectrum.
The signal model (7) can be introduced into the source component x in (4). Again, n denotes the index of the grids. When xn,k,t
represents the complex spectrogram in the time-frequency domain ateach grid, xn,k,t can be described as∑
xn,k,t ≈hk,l un,l,t ejϕn,k,t .
(8)
l

n=1

346

We assume that the magnitude spectrum hk,l is static and trained inadvance using possible monaural source signals. Although ϕn,k,tshould have the basis index l when the CNMF model is straightforwardly applied, we omit this index and assume that each grid has its
own phase spectrogram. Denoting the microphone index as m, theobserved signal of the mth microphone, ym,k,t , can be represented
as∑∑ym,k,t ≈dm,n,khk,l un,l,t ejϕn,k,t ,(9)nf + (U, Φ, Y, U)∑ |y m,n,k,l,t − dm,n,k hk,l un,l,t ejϕn,k,t |2=βm,n,k,l,t
m,n,k,l,t)∑{ (∑ 2 )p/2−1 ( ∑ 2+λpun,l,tun,l,t − u2n,l,tlwhere dm,n,k is the element of the dictionary matrix between thenth grid and the mth microphone in the kth frequency bin.
Finally, we formulate the objective function under the assumption that dm,n,k is known and the basis magnitude spectrum hk,l istrained in advance. We define the modeling error of (9) as
R(U, Φ) =∑∑∑hk,l un,l,t ejϕn,k,t |2 ,|ym,k,t −dm,n,knm,k,tupper bound function of f (U, Φ, Y, U), f + (U, Φ, Y, U), can be
defined as

(10)

l

n

l,t

∑(∑
n

u2n,l,tThis function satisfies f (U, Φ) ≤ f + (U, Φ, Y, U), and the equality is satisfied when
y m,n,k,l,t = dm,n,k hk,l un,l,t ejϕn,k,t)(∑
dm,n,k hk,l un,l,t ejϕn,k,t+βm,n,k,l,t ym,k,t −.

un,l,t = un,l,t .

l,t

(12)

U,Φ

where λ is a parameter used to balance the modeling error R(U, Φ)and the penalty J p,2 (U).
4. OPTIMIZATION ALGORITHM BASED ON AUXILIARYFUNCTION METHOD
We derive an optimization algorithm using the auxiliary function
method [16], which can give stable and fast update rules. First, wedevelop an auxiliary function of the objective function using the following inequalities:
dm,n,kn≤∑hk,l un,l,t ejϕn,k,t |2l∑ |y m,n,k,l,t − dm,n,k hk,l un,l,t ejϕn,k,t |2βm,n,k,l,t

n,l

,

(13)

and
∑ 2
2(
un,l,t )p/2
l,t

≤ p(

∑
l,t

(16)

(11)

{
}
minimize R(U, Φ) + 2λJ p,2 (U) ,

∑

n,l

and

)p/2

Therefore, the optimization criterion is formulated as|ym,k,t −

(15)

l,t

where U and Φ represent third-rank tensors whose elements are
un,l,t and ϕn,k,t , respectively. To obtain a spatially sparse solutionof U, the following penalty term is introduced:
J p,2 (U) =l,t(∑)p/2 }+2u2n,l,t.

(17)

Therefore, the conditions required for f + (U, Φ, Y, U) to be an
auxiliary function of f (U, Φ) are satisfied [16]. Here, Y and U are
tensors whose elements consist of y m,n,k,l,t and un,l,t , respectively.
Using the auxiliary function (15), we can derive iterative update
rules for U and Φ. Although un,l,t should be updated under thenonnegative constraint, we here update un,l,t without this constraint.
First, by solving ∂f + /∂un,l,t = 0, we obtain the following update
rule:
∑

∗
jϕn,k,t]/βm,n,k,l,t
m,k ℜ[y m,n,k,l,t dm,n,k hk,l e∑.
2
2 h2 /β
p/2−1
|d
|
+
λp(
m,n,k
m,n,k,l,tk,lm,k
l,t un,l,t )

un,l,t ← ∑(18)Here, ℜ[·] denotes the operator that returns the real part of a complex
number and ∗ denotes the operator that returns the complex conjugate value. To obtain the update rule for Φ, we use the following∑∑′
properties. When θ′ satisfies ejθ = i ai b∗i /| i ai b∗i |, θ′ mini∑
mizes i |ai − bi ejθ |2 , where ai , bi ∈ C and θ ∈ R. Then, theupdate rule is obtained as∑∗m,l y m,n,k,l,t dm,n,k hk,l un,l,t /βm,n,k,l,t. (19)
ejϕn,k,t ← ∑| m,l y m,n,k,l,t d∗m,n,k hk,l un,l,t /βm,n,k,l,t |We here set the arbitrary variable βm,n,k,l,t to∑ 2∑ 2
u2n,l,t )p/2−1 (
un,l,t − u2n,l,t ) + 2(
un,l,t )p/2 .
l,tβm,n,k,l,t = ∑|dm,n,k |hk,l un,l,t
n,l |dm,n,k |hk,l un,l,t

.

(20)

l,t

(14)
Here, y m,n,k,l,t and un,l,t are auxiliary variables and βm,n,k,l,t
is∑ an arbitrary constant. The inequality (13) is satisfied when(14) is satisfied whenn,l y m,n,k,l,t = ym,k,t , and the inequality∑
0 < p ≤ 2, 0 ≤ βm,n,k,l,t ≤ 1, and n,l βm,n,k,l,t = 1. Then, an347

We summarize the proposed algorithm as follows. In the training stage, the magnitude spectrum bases hk,l are trained in advanceusing possible monaural source signals and dm,n,k is obtained by
setting the locations of microphones and grids. In the decomposition stage, U and Φ are first initialized. Then, Y, U, U, and Φ areiteratively updated using (16), (17), (18), and (19), respectively.

5. EXPERIMENTSl,t

For a quantitative evaluation, we define the F -measure as
Fmsr = 2|supp{Pest (θn )} ∩ supp{Ptrue (θn )}|,|supp{Pest (θn )}| + |supp{Ptrue (θn )}|(22)

where Pest (θ) and Ptrue (θ) are the estimated and true DOA spectra,
respectively. The operator supp(·) extracts a set of indexes suchthat the value of each element of the solution P (θn ) is larger than athreshold value µ,supp{P (θn )} = {n ∈ {1, ..., N } | P (θn ) > µ}.

(23)

Since there exists no optimal threshold value for all the methods, wehere changed µ from –15 dB to –1 dB in this experiment. The trueindexes of θn were set to those for the two grids nearest the true
angle. Fmsr was averaged over 100 trials, i.e., 100 pairs of DOAs
randomly selected.
We also evaluated the signal reconstruction performance, using
the signal-to-distortion ratio (SDR) defined as∑
2
m,k,t |yorg,m,k,t |SDR = 10 log10 ∑,(24)
2
m,k,t |yorg,k,m,t − yest,m,k,t |

348

-5
-10
Proposed
M-FOCUSS
DAS
Capon
MUSIC

-15
-20

-80

-60

-40

-20

0

20

40

60

80

Angle [degrees]Fig. 2. Examples of DOA spectra for each method when the true
angles were set to –43.5◦ and 54.2◦ .
0.6

Proposed
M-FOCUSS
DAS
Capon
MUSIC

0.4

Fmsr

Numerical experiments on DOA estimation and signal reconstruction in terms of plane wave decomposition were conducted to evaluate the proposed method. We compared the proposed method (Proposed), M-FOCUSS [18], which is a sparse decomposition method
for the MMV problem based only on the simultaneous spatial sparsity of the sources, and conventional DOA estimation methods, i.e.,
DAS, the Capon beamformer (CB) [7], and MUSIC.
Omnidirectional microphones were linearly aligned along the xaxis with their center at the origin. The number of microphones wasfour and they were arranged at intervals of 0.24 m; therefore, thearray length was 0.72 m. Two plane wave sources, whose angles ofarrival were randomly determined within the range of (−90◦ , 90◦ ),impinged on the microphones. Identical source signals were generated by artificial MIDI sounds, i.e., the two source signals were
the same. The musical instrument used to generate the source signals as above was the oboe. The sampling frequency was 16 kHz.
The short-time Fourier transform (STFT) was performed to obtaintime-frequency spectrograms of the observed signals. A square-root
Hanning window of 32 ms length with a 16 ms overlap was usedin the STFT. The training data were two octaves of notes that covered all the notes of the source signal. They were also generated byartificial MIDI sounds. The number of bases was 25. A white Gaussian noise was added to the signals received by the microphones sothat the signal-to-noise ratio became 20 dB. The dictionary matrix,
i.e., the steering vector for the conventional DOA estimation methods, was generated so that the intervals between the grids were 1.0◦
within the range of (−90◦ , 90◦ ) as shown in Fig. 1. The angles ofeach grid is denoted as θn .
In M-FOCUSS, the parameters for the row-sparse penalty (p andq in [18]) were set to p = 1 and q = 2. The maximum number of
iterations was 400. The parameter corresponding to η in (5) was
set to 1.0. In Proposed, λ in (12) was set to 10−5 in the first five
iterations, and 2.0 × 102 in the remaining 95 iterations.
To evaluate the source localization performance, we used the
DOA spectrum [7], which indicates the power distribution of planewaves from each angle θn . For example, in Proposed, the DOA
spectrum P (θn ) was obtained as∑ 2
P (θn ) = 10 log10 un,,l,t  .
(21)

DOA spectrum [dB]

0

0.2

0
-15

-13

-11

-9

-5
-3
-1
µ -7
Fig. 3. Relation between Fmsr and threshold value µ.

where yorg,m,k,t is the original signal in the observations andyest,m,k,t is the signal in the observations reconstructed using thedecomposed signals.
Figure 2 shows examples of the DOA spectra of Proposed, MFOCUSS, DAS, CB, and MUSIC when the true angles were –43.5◦
and 54.2◦ . The DOA spectra of the conventional methods (DAS, CB,
and MUSIC) had a false peak around 0◦ . This is because these methods were derived under the assumption that the source signals areuncorrelated. Since the source signals were identical in this experimental setup, it is very difficult to accurately estimate the true source
angles. In the DOA spectra of M-FOCUSS and Proposed, the spatial
resolution was significantly improved by assuming the spatial sparsity of the sources. However, the DOA spectrum of M-FOCUSS stillcontained multiple false peaks. In contrast, that of Proposed had two
peaks at accurate angles. The average values of Fmsr were shown in
Figure 3. When the true angles were –43.5◦ and 54.2◦ , the SDRs of
Proposed and M-FOCUSS were 5.27 dB and 2.13 dB, respectively.
This indicates the efficacy of the proposed method, particularly forthe decomposition of multiple coherent sources.
6. CONCLUSION
We proposed a sparse sound field decomposition method using priorinformation on source signals in the time-frequency domain. The
proposed method was derived by incorporating the CNMF model ofsource signals into sparse sound field decomposition. Whereas current methods are based only on the spatial sparsity of the source distribution, the proposed method exploits magnitude spectrum bases
pretrained using possible monaural source signals. The decomposition algorithm was derived using the auxiliary function method.
Numerical simulation results indicated that the DOA estimation andsignal reconstruction performance of the proposed method was significantly improved from that of current methods, particularly whenthe multiple source signals were highly correlated.

7. REFERENCES[1] E. G. Williams, Fourier Acoustics: Sound Radiation and
Nearfield Acoustical Holography, Academic Press, New York,1999.
[2] S. Koyama, K. Furuya, Y. Hiwasaki, and Y. Haneda, “Analytical approach to wave field reconstruction filtering in spatiotemporal frequency domain,” IEEE Trans. Audio, Speech,
Lang. Process., vol. 21, no. 4, pp. 685–696, 2013.
[3] S. Koyama, K. Furuya, Y. Hiwasaki, Y. Haneda, and Y. Suzuki,“Wave field reconstruction filtering in cylindrical harmonic domain for with-height recording and reproduction,” IEEE/ACM
Trans. Audio, Speech, Lang. Process., vol. 22, no. 10, pp.
1546–1557, 2014.
[4] J. Daniel, “Spatial sound encoding including near field effect:
Introducing distance coding filters and a viable, new ambisonics format,” in Proc. 23rd AES Int. Conf., Copenhagen, May
2003.
[5] M. Poletti, “Three-dimensional surround sound systems basedon spherical harmonics,” J. Audio Eng. Soc., vol. 53, no. 11,
pp. 1004–1025, 2005.
[6] M. Brandstein and W. Darren, Microphone arrays: signal processing techniques and applications, Springer Science & Business, 2001.
[7] V. Trees and L. Harry, Optimum array processing: part IV ofdetection, estimation, and modulation, John Wiley & Sons,
2002.
[8] R. O. Schmidt, “Multiple emitter location and signal parameter
estimation,” IEEE Trans. Antennas Propag., vol. 34, no. 3, pp.
276–280, 1986.
[9] G. Chardon, L. Daudet, A. Peillot, F. Ollivier, N. Bertin,
and N. Gribonval, “Near-field acoustic holography using
sparse regularization and compressive sampling principles,” J.
Acoust. Soc. Am., vol. 132, no. 3, pp. 1521–1534, 2012.
[10] S. Koyama, S. Shimauchi, and H. Ohmuro, “Sparse soundfield representation in recording and reproduction for reducing
spatial aliasing artifacts,” in Proc. IEEE Int. Conf. Acoust.,
Speech, Signal Process. (ICASSP), Florence, May 2014, pp.
4476–4480.
[11] S. Koyama, N. Murata, and H. Saruwatari, “Structured
sparse signal models and decomposition algorithm for superresolution in sound field recording and reproduction,” in Proc.
IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP),
Brisbane, Apr. 2015, pp. 619–623.
[12] D. Malioutov, M. Çetin, and A. S. Willsky, “A sparse signal
reconstruction perspective for source localization with sensor
arrays,” IEEE Trans. Signal Process., vol. 53, no. 8, pp. 3010–
3022, 2005.
[13] A. Asaei, H. Bourlard, M. Taghizadeh, and V. Cevher, “Modelbased sparse component analysis for reverberant speech localization,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), Florence, May 2014, pp. 1453–1457.
[14] T. Noohi, N. Epain, and C. T. Jin, “Super-resolution acoustic
imaging using sparse recovery with spatial priming,” in Proc.
IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP),
Apr. 2015, pp. 2414–2418.

349[15] P. Smaragdis and J. C. Brown, “Non-negative matrix factorization for polyphonic music transcription,” in Proc. IEEE Int.
Workshop Appl. Signal Process. Audio Acoust. (WASPAA), Oct.
2003, pp. 177–180.
[16] H. Kameoka, N. Ono, K. Kashino, and S. Sagayama, “Complex NMF: A new sparse representation for acoustic signals,”in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process.
(ICASSP), Taipei, Apr. 2009, pp. 3437–3440.
[17] M. Elad, Sparse and Redundant Representations: From Theory
to Applications in Signal and Image Processing, Springer, New
York, 2010.
[18] S. F. Cotter, B. D. Rao, K. Engan, and K. Kreutz-Delgado,“Sparse solutions to linear inverse problems with multiple measurement vectors,” IEEE Trans. Signal Process., vol. 53, no. 7,
pp. 2477–2488, 2005.
[19] A. Rakotomamonjy, “Surveying and computing simultaneous
sparse approximation (or group-lasso) algorithms,” Signal Process., vol. 91, pp. 1505–1526, 2011.