COMPLEX NMF: A NEW SPARSE REPRESENTATION FOR ACOUSTIC SIGNALS
Hirokazu Kameoka† , Nobutaka Ono‡ , Kunio Kashino† , Shigeki Sagayama‡† NTT Communication Science Laboratories, NTT Corporation,
3-1 Morinosato Wakamiya, Atsugi, Kanagawa 243-0198, Japan‡ Graduate School of Information Science and Technology, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan

ABSTRACTThis paper presents a new sparse representation for acoustic signals which is based on a mixing model defined in thecomplex-spectrum domain (where additivity holds), and allows us to extract recurrent patterns of magnitude spectra thatunderlie observed complex spectra and the phase estimatesof constituent signals. An efficient iterative algorithm is derived, which reduces to the multiplicative update algorithm
for non-negative matrix factorization developed by Lee undera particular condition.
Index Terms— Sparse signal representation, non-negative
matrix factorization, sparse coding, data-driven approach
1. INTRODUCTION
The use of sparse representations including sparse coding(SC) [1] and non-negative matrix factorization (NMF) [2]in acoustic signal processing has gradually increased in recent years. Given a set of observed vectors, the goal is tofind a set of basis vectors such that any observations canbe succinctly represented as a linear combination of a small
number of ‘active’ bases. While SC enforces the sparsenessof basis activations by incorporating a sparsity cost in theobjective function, NMF achieves sparseness as a result ofa side effect caused by non-negativity constraints. Unlike
certain other unsupervised basis learning techniques such as
principal component analysis (PCA) and independent component analysis (ICA), the particularity of SC and NMF is thatrecurrent patterns underlying the observations are consideredgood bases owing to the sparse nature of the decompositions.
[3]–[6] were among the first to apply this concept to an audio
processing problem, and since then, NMF, in particular, hasbeen applied extensively with considerable success to various
problems including automatic music transcription, monauralsource separation, speech denoising, bandwidth expansion,audio classification and speech dereverberation [7].
One reason for NMF being a powerful tool as regards extracting regularities or structural patterns from acoustic signals may be that the mixing model is defined in the magnitude spectrum domain. Because of the phase-invariant natureof magnitude spectra, NMF is able to project all signals thathave the same spectral shape onto a single basis. This allows
us to represent a variety of acoustic phenomena efficiently
using a very compact set of spectrum bases. However, themixing model implicitly assumes the additivity of magnitude
spectra, which holds only approximately. Although attempts

978-1-4244-2354-5/09/$25.00 ©2009 IEEE3437

are made to mitigate the non-additivity problem with respect
to NMF in certain recently published papers [8, 9], in these
papers the additivity of power spectra is only verified in a statistical sense on the assumption that the constituent signals
are Gaussian processes and perfectly independent. Anotherdrawback of NMF is that it cannot estimate the phase spectraof underlying constituent signals, which certainly limits its
range of applications.
Within the framework of SC or semi-NMF we can construct a mixing model defined in the time domain [10, 11] orin the complex-spectrum domain. Such a mixing model isvalid in terms of defining an additive quantity as the basis,
namely, a waveform signal. However, as pointed out in [6]because the phase coherence between frequency components
can be easily destroyed as a result of many factors, it is difficult to capture high-level structural elements from observations through the use of waveform/complex-spectrum bases.
In this paper, we present a new framework “ComplexNMF”, which offers the advantages of the NMF and SC
frameworks concurrently: Complex NMF, (1) is based ona mixing model defined in the complex-spectrum domain(where additivity holds), (2) can extract the recurrent patterns
of magnitude spectra that underlie the observed complex
spectra and the phase spectra of constituent signals, and (3)can be performed with an efficient iterative algorithm, which
reduces to the multiplicative update algorithm developed by
Lee [2] under a particular condition.
2. PRINCIPLE OF COMPLEX NMF
2.1. Sparse representation model
We start with the basic assumption that the short-term Fouriertransform (STFT) of an arbitrary acoustic signal, Fx,t ∈ C,
consists of K complex-valued elements
K

Fx,t =

k=1

ak,x,t =|ak,x,t |ejφk,x,t ,(1)kwhere x and t are the frequency and frame indices, respectively. Here, we factorize the modulus, |ak,x,t |, into the product of nonnegative parameters, Hk,x and Uk,t ,|ak,x,t | = Hk,x Uk,t , Hk,x ≥ 0, Uk,t ≥ 0,and assume

Hk,x = 1 (k = 1, · · · , K),(2)(3)xICASSP 2009

in order to avoid an indeterminacy in the scaling. We therefore arrive at the following mixing modelFx,t =Hk,x Uk,t ejφk,x,t .
(4)k

In this model, each element comprises a static magnitude
spectral shape Hk,x , a time-varying activation coefficient
Uk,t and a time-varying phase spectrum φk,x,t . Given an observed spectrum Yx,t , the goal is to find a decomposition suchthat Yx,t  Fx,t in which the basis activations are sparse.
This basically means that any observed complex spectrumcan be well represented using only a few active magnitude
spectrum bases each of which is paired with an arbitrary
phase spectrum. Combining the goal of a small reconstruction error with that of sparseness, we consider that the model
will allow us to extract the recurrent patterns of magnitude
spectra underlying observed spectra as with NMF and, simultaneously, the phase spectra of constituent signals. As
seen from Eq. (4), the model cannot be expressed in matrix
notation as with NMF and SC, it can clearly be considered anew class of sparse representation model.
2.2. Problem setting
Given an observed complex spectrum, Yx,t ∈ C, we would
like to find the optimal estimates of Hk,x , Uk,t and φk,x,t . For
simplicity of notation let Y ≡ {Yx,t }X×T , F ≡ {Fx,t }X×T ,H ≡ {Hk,x }K×X , U ≡ {Uk,t }K×T , φ ≡ {φk,x,t }K×X×T .
Now, we assume the following generative modelYx,t = Fx,t + x,t ,(5)where the reconstruction error x,t is assumed to be complex
Gaussian white noise with mean 0 and variance σ 2 . The likelihood of θ = {H, U, φ} is thus given as follows: 1|Yx,t − Fx,t |2
P (Y |θ) =.
(6)
exp −πσ 2σ2x,tWe assume for convenience that the prior distributionsfor H, U and φ are independent, which yields P (θ) =P (H)P (U )P (φ), and that P (H) and P (φ) are uniform distributions. P (U ) corresponds to the sparsity cost, for which anatural choice is a generalized Gaussian prior|Uk,t |p1
 exp − p
P (U ) =,(7)b2Γ 1 + p1 bk,twhere p and b are the parameters that determine the shape ofthe distribution. When 0 < p < 2, P (U ) becomes superGaussian and promotes sparsity if the norm of U is bounded.
The likely values of H, U and φ can thus be inferred from theposterior density
P (θ|Y ) ∝ P (Y |θ)P (U ).

(8)

From Eqs. (3)–(8), we are thus led to solve the following
optimization problem:Yx,t − Fx,t 2 + 2λminimize f (θ) ≡|Uk,t |psubject to

x,tk,tHk,x = 1 (k = 1, · · · , K)

2.3. Iterative algorithm
We present an efficient algorithm that seeks to minimize f (θ).
The auxiliary function concept is utilized for our derivation,similar to [2]. We use G(θ) to denote an objective functionthat we want to minimize w.r.t. θ, and define an auxiliary
function of G(θ) as G+ (θ, θ̄) if it satisfies
G(z) = min G+ (z, z̄).

Theorem 1. G(θ) is non-increasing under the updates, θ̄ ←argminθ̄ G+ (θ, θ̄) and θ ← argminθ G+ (θ, θ̄).
Proof: Assume that we set θ at an arbitrary value θ . Letθ̄+1 = argminθ̄ G+ (θ , θ̄) and θ+1 = argminθ G+ (θ, θ̄+1 ).
It is obvious that G(θ ) = G+ (θ , θ̄+1 ). We deducefrom θ+1 = argminθ G+ (θ, θ̄+1 ) that G+ (θ , θ̄+1 ) ≥G+ (θ+1 , θ̄+1 ) . By definition, from G+ (θ+1 , θ̄+1 ) ≥G(θ+1 ) we verify that G(θ ) ≥ G(θ+1 ).
By iteratively updating θ and θ̄, G(θ) will converge to astationary point. We can apply this concept to the minimization of f (θ) using the following theorem:
Theorem 2 (Auxiliary function). The function+f (θ, θ̄) ≡+λ3438

 Ȳk,x,t − Hk,x Uk,t ejφk,x,t 2βk,x,tk,x,tp|Ūk,t |p−2 Uk,t 2 + 2|Ūk,t |p − p|Ūk,t |p , (10)k,twith θ̄ = {Ȳ , Ū }, Ȳ ≡ {Ȳk,x,t }K×X×T , Ū ≡ {Ūk,t }K×T ,is an auxiliary function for f (θ), if k Ȳk,x,t = Yx,t and0 < p ≤ 2. βk,x,t can be any positive number satisfying+k βk,x,t = 1. f (θ, θ̄) is minimized w.r.t. θ̄ whenȲk,x,t = Hk,x Uk,t ejφk,x,t + βk,x,t Yx,t − Fx,t ,(11)Ūk,t = Uk,t .

(12)

Proof: This follows from the following lemmas.
Lemma 1. With k Ȳk,x,t = Yx,t and for any βk,x,t > 0such that k βk,x,t = 1, Ȳk,x,t − Hk,x Uk,t ejφk,x,t 2
Yx,t − Fx,t 2 ≤. (13)βk,x,tk

Proof: By adding the Lagrange multiplier term to the right∗
andhand side, and then taking the derivative w.r.t. Ȳk,x,t
setting it at zero, we obtain Eq. (11). Substituting this to
the right-hand side, we determine the minimum value of theright-hand side, which is equal to the left-hand side of theinequality.
Lemma 2. If 0 < p ≤ 2, for any Uk,t ∈ R and Ūk,t ∈ R,|Uk,t |p ≤x(9)z̄

p|Ūk,t |p−2
p|Ūk,t |pUk,t 2 + |Ūk,t |p −.
22(14)

Proof: Regarding |Uk,t |p as the function of Uk,t , the righthand side of the inequality amounts to a convex quadratic
function that is tangent to |Uk,t |p at argument Uk,t = ±Ūk,t .

Substituting this into Eq. (11),
Hk,x Uk,t+ βk,x,t 1 −Ȳk,x,t = Yx,t|Yx,t |As the update rule for θ̄ is shown in Eqs. (11), (12), weonly need to derive the update rule for θ. Although it is necessary to take account of the normalization condition for Hk,xand solve the Lagrange dual problem, here we simply describe a convenient approach, which consists of solving theunconstrained minimization of f + (θ, θ̄) w.r.t. Hk,x and projecting its solution onto the constraint space. Differentiatingf + (θ, θ̄) partially w.r.t. Hk,x and Uk,t , and setting them atzero, we obtain update rules for Hk,x and Uk,t :Substituting this result and Eq. (18) into Eq. (15),2 Hk,x Uk,t+ Uk,t |Yx,t | −Hn,x Un,tβk,x,ttnHk,x ←.
 Uk,t 2

 Uk,tHk,x =βk,x,t Uk,t 2t Hk,xUk,t =∗Re Ȳk,x,tejφk,x,ttx(15)βk,x,tβk,x,t.
+ λp|Ūk,t |k,x,twhere Ak,x,t = Ȳk,x,t Hk,x Uk,t /βk,x,t andRe Ȳk,x,tIm Ȳk,x,tcos Ck,x,t =, sin Ck,x,t =.
|Ȳk,x,t ||Ȳk,x,t |f + (θ, θ̄) is obviously minimized when cos(φk,x,t −Ck,x,t ) =cos φk,x,t cos Ck,x,t + sin φk,x,t sin Ck,x,t = 1, that is,cos φk,x,t = cos Ck,x,t and sin φk,x,t = sin Ck,x,t . This
leads to the update formula for ejφk,x,teȲk,x,t=.
|Ȳk,x,t |(17)

Substituting Eq. (17) into Eqs. (15), (16), we notice that the
non-negativity of Hk,x and Uk,t is preserved if we start with
non-negative initial conditions for Hk,x and Uk,t .
2.4. Condition for equivalence to NMF
The iterative algorithm presented in 2.3 reduces to the Lee’s
multiplicative update algorithm [2] under a particular condition. Let us consider a particular situation where we fix thevalue of φ, from the beginning of the iteration, atejφk,x,t =Yx,t.
|Yx,t |

(18)

3439

.

βk,x,tt

As βk,x,t is a parameter that can be chosen arbitrarily subject
to βk,x,t > 0 and k βk,x,t = 1, the convergence of theiterative algorithm is still guaranteed even if we change its
value in parallel with H and U in each iteration. If we decideto replace the value of βk,x,t at each iteration by(19)

n(16)

p−2

We note that these updates are guaranteed to provide the minimum value for each coordinate because it is obvious from aninspection of their second derivatives. Next, we derive the update rule for φk,x,t . By using c to denote the terms that do notdepend on φk,x,t , f + (θ, θ̄) can be simply written as follows
|Ak,x,t | cos(φk,x,t − Ck,x,t ),f + (θ, θ̄) = c − 2

jφk,x,t|Yx,t |Hk,x Uk,tβk,x,t = ,Hn,x Un,t∗
Re Ȳk,x,tejφk,x,tβk,x,t Hk,x 2x,n Hn,x Un,tthe update formula for Hk,x leads eventually to
Uk,t |Yx,t |Hk,x ← Hk,x ttUk,tHn,x Un,t,(20)

nwhich amounts to Lee’s multiplicative update rule for theFrobenius norm criterion. We obtain the same conclusion for
Uk,t , if λ = 0. Based on this analysis, we can conjecture that(1) the present algorithm works as efficiently and stably as
the Lee’s NMF algorithm,(2) updating βk,x,t according to Eq. (19) is somewhat more
effective than fixing it at a constant value, and(3) we can stabilize the algorithm by running NMF at thebeginning of the iteration, which can be performed simply by fixing the value of φ at Eq. (18).
The iterative algorithm is summarized as follows:
1. Initialize H, U and φ.
2. Update θ̄ = {Ȳ , Ū } according to Eqs. (11), (12).
3. Update θ = {H, U, φ} according to Eqs. (15), (16),(17) and Hk,x ← Hk,x / n Hn,x .
4. Update β according to Eq. (19) and return to 2.
3. EXPERIMENTS
In this section we report some results for speech data excerpted from the ATR B-set speech database. The aim of theexperiments was to determine whether Complex NMF has aneffect as with NMF to extract the underlying patterns of magnitude spectra from audio data. All speech data were monaural and sampled at 16kHz. STFT was computed using a Hanning window that was 32ms long with a 16ms overlap. p andλ were set at p = 1.2，λ = x,t |Yx,t |2 /K 1−p/2 ×10−5 . The

250

Frequency bin #Frequency bin #250

200200

150

150

100

100
505050100

150200250
Frame #

300

350

400

50450

Fig. 1. Magnitude spectrogram of speech A100150200250
Frame #300250
Frame #300350400450Fig. 3. Magnitude spectrogram of mixed speech
Frequency bin #2502001501005050100150200350400450Fig. 4. Synthesized signal corresponding to speech A
4. CONCLUDING REMARKS

Fig. 2. 10 patterns of magnitude spectra and their activation patterns
algorithm was run for 30 iterations. H and U were initializedwith random values and φk,x,t was initially set at arg(Yx,t ).
For the first experiment we tested Complex NMF, with
K = 10, on a single voice signal from a female speaker(“speech A”). Fig. 1 shows the spectrogram of speech A. The
patterns of magnitude spectra discovered by using Complex
NMF can be seen in Fig. 2 along with the corresponding activation patterns. As shown in Fig. 2, harmonic structures arevividly captured in H even without any prior knowledge.
For the second experiment we tested Complex NMF with
K = 30 on a mixed voice signal (see Fig. 3). The mixed
speech data was created by adding speech data from anotherfemale speaker (“speech B”) to speech A. In this experiment,
the aim was to ascertain whether each of the estimated patterns of magnitude spectra corresponds to a single-voice spectrum. We selected the pattern of magnitude spectrum that wasclosest to the true spectrum of each speaker per frame andthen concatenated the framewise signals, each of which weconstructed using the selected pattern and the corresponding
activation coefficient and phase spectrum to synthesize thewhole signal stream. If the synthesized signal contains twovoices, it may suggest that the model is overfitting the observation. In contrast, if the synthesized signal does not soundlike speech, it may suggest that the model is underfitting theobservation. Seen in this light, we calculated SNR to see how
well each speech signal could be restored based on the above
procedure. As a result, we obtained 4.3dB and 4.0dB for thetwo speakers, which constituted an improvement from 0.87dB
and -0.87dB, respectively. This result suggests that each estimated pattern of magnitude spectrum corresponds fairly appropriately to a single voice spectrum. The restored speech
signal corresponding to speech A can be seen in Fig. 4.

3440

We developed a new framework for the sparse representationof acoustic signals. We named it “Complex NMF” (whichmay appear contradictory since the term “complex” conflicts with the term “non-negative” and the decomposition is
no longer a matrix factorization), because, in common withNMF, it has the ability to generate non-negative matrices H
and U , and is different in that the input matrix Y is assumed
to be a complex matrix and it also generates a third-rankcomplex-valued tensor ejφk,x,t . As Complex NMF is capableof estimating the phase spectra of constituent signals, future
work will include its extension to the formulation for microphone array processing. Motivated by [12], we also want toinvestigate a general class of the prior structure suited to themodel introduced in this paper.
5. REFERENCES[1] B.A. Olshausen and D.J. Field, “Emergence of simple-cell receptive

field properties by learning a sparse code for natural images,” Nature,381, 607–609, 1996.
[2] D.D. Lee and H.S. Seung, “Algorithms for nonnegative matrix factorization,” Proc. NIPS’00, 556–562, 2000.
[3] M.D. Plumbley, S.A. Abdallah, J.P. Bello, M.E. Davies, G. Monti andM.B. Sandler, “Automatic music transcription and audio source separation,” Cybernetics and Systems, 33(6), 603–627, 2002.
[4] P. Smaragdis and J.C. Brown, “Non-negative matrix factorization formusic transcription,” Proc. WASPAA’03, 177–180, 2003.
[5] T. Virtanen, “Sound source separation using sparse coding with temporal continuity objective,” Proc. ICMC’03, 231–234, 2003.
[6] S.A. Abdallah, M.D. Plumbley, “Unsupervised analysis of polyphonicmusic using sparse coding,” IEEE Trans. NN, 17(1), 179–196, 2006.
[7] H. Kameoka T. Nakatani and T. Yoshioka, “Robust speech dereverberation based on non-negativity and sparse nature of speech spectrograms,” Proc. ICASSP’09, to appear.
[8] R.M. Parry and I. Essa, “Phase-Aware Non-negative Spectrogram Factorization,” Proc. ICA’07, 536–543, 2007.
[9] C. Févotte, N. Bertin and J.-L. Durrieu, “Nonnegative matrix factorization with the Itakura-Saito divergence with application to music
analysis,” Technical Report TELECOM ParisTech 2008D006, 2008.
[10] T. Blumensath and M.E. Davies, “Sparse and shift-invariant representations of music,” IEEE Trans. ASLP, 14(1), 50–57, 2006.
[11] J. Le Roux, A. de Cheveigne and L.C. Parra, “Adaptive templatematching with shift-invariant semi-NMF,” Proc. NIPS’08, 2008.
[12] A.T. Cemgil, “Bayesian inference for nonnegative matrix factorisation
models,” Technical Report CUED/F-INFENG/TR.609, University ofCambridge, 2008.