MULTI-RESOLUTION SIGNAL DECOMPOSITION WITH
TIME-DOMAIN SPECTROGRAM FACTORIZATION
Hirokazu Kameoka
The University of Tokyo / Nippon Telegraph and Telephone Corporation

ABSTRACTThis paper proposes a novel framework that makes it possible to realize non-negative matrix factorization (NMF)-likesignal decompositions in the time-domain. This new formulation also allows for an extension to multi-resolution signal
decomposition, which was not possible with the conventional
NMF framework.
Index Terms— Audio source separation, non-negative
matrix factorization (NMF), majorization-minimization (MM),auxiliary function, multi-resolution representation
1. INTRODUCTION
Many sound recordings are mixtures of multiple soundsources. Audio source separation, i.e., the process by which
individual sound sources are separated from a mixture signal,
has long been a formidable challenge in the field of audio
signal processing.
In recent years, non-negative matrix factorization (NMF)
has attracted a lot of attention after being proposed as a powerful approach for music transcription [1] and audio source
separation [2]. With this approach, the magnitude (or power)
spectrogram of a mixture signal, interpreted as a non-negative
matrix Y , is factorized into the product of two non-negative
matrices H and U . This can be interpreted as approximating
the observed spectra at each time frame as a linear sum of basis spectra scaled by time-varying amplitudes, and amounts todecomposing the observed spectrogram into the sum of rank1 spectrograms. An important feature of NMF is that its nonnegativity constraint usually induces sparse representations,
i.e., U with a relatively large number of zero entries. Thismeans that each observed spectrum is parsimoniously represented using only a few active basis spectra. In such situations, the sequence of observed spectra can be approximatedreasonably well when each basis spectrum expresses the spectrum of an underlying audio event that occurs frequently in theentire observed range. Thus, with music signals, each basis
spectrum usually becomes the spectrum of a frequently used
pitch in the music piece.
Although the concept of the NMF-based audio source separation approach has been shown to be successful, one limitation is that it does not take account of phase information: theadditivity of magnitude (or power) spectra is assumed, whichholds only approximately. To overcome this limitation, we
have previously proposed a framework called the “complex
NMF” [3], where the observed complex spectrum at each timeframe is modeled as the sum of components, each of which is
described by the multiplication of a static basis spectrum, atime-varying amplitude and a time-varying phase spectrum.

2. FORMULATIONLet us denote an observed signal at time tn by y[n], andthe signal of the entire period by y = (y[1], . . . , y[N ])T ∈RN . While the NMF approach considers decomposing an observed magnitude spectrogram into the sum of rank-1 spectro-

This work was supported by JSPS KAKENHI Grant Number 26730100.

978-1-4673-6997-8/15/$31.00 ©2015 IEEE

Unlike NMF, this model allows the components to cancel
each other out, and so without any constraints, it does not naturally produce sparse representations. However, it has beenshown that an additional sparsity constraint yields sparse representations similar to NMF. With a similar motivation, twogroups (Parry et al. and Févotte et al.) have independently
proposed a generative model of complex-valued coefficientsof the short-time Fourier transform (STFT) of a mixture signal, where the power spectrogram of each underlying component is modeled as a rank-1 matrix (similarly to NMF) and thephase spectrogram is treated as uniformly distributed latentvariables [4, 5]. They showed that the maximum likelihood
estimation of the power spectral density parameters amounts
to fitting the NMF model to an observed power spectrogramusing the Itakura-Saito (IS) divergence as a goodness-of-fit
criterion. This approach is called IS-NMF.
To the best of our knowledge, complex NMF and ISNMF are among the first phase-aware NMF variants to be
proposed. Although both approaches treat each element of
the phase spectrogram as an independent parameter (or latentvariable), the phases of time-frequency components are, infact, constrained and dependent on each other. This is because the spectrograms obtained with typical time-frequency
transforms (such as the STFT and the wavelet transform)
are redundant representations (see Fig. 1 regarding STFT
spectrograms). For example, the STFT spectrogram is computed by concatenating the Fourier transforms of overlapping
short-time frames of the signal. Hence, all the elements of
the STFT spectrogram must satisfy a certain condition toensure that the waveforms within the overlapping segment ofconsecutive frames are consistent [6]. The shortcomings ofthe complex NMF and IS-NMF frameworks are that they failto take account of this kind of redundancy. One possible way
to further develop improved variants of these models would
be to incoporate into the models the explicit condition thatspectrograms must satisfy. However, this would make thesemodels overcomplicated and difficult to optimize.
Instead of using mixing models defined in the timefrequency domain, this paper proposes introducing a timedomain model that makes it possible to realize NMF-likesignal decompositions. Moreover, we will show that the
time-domain formulation also allows for an extension tomulti-resolution signal decompositions, that was not possible
with the conventional frameworks.

86

ICASSP 2015

Complex
spectrogram

En re signal...

Frame 1

Frame 1...

Frame 2

En re signal ofl-th component
Basis func onwith frequency kand me m

Frame 2

...

..
.
...

Fig. 2. Illustration of time-frequency basis functions.

..
.

3. OPTIMIZATION ALGORITHM
3.1. General principle of auxiliary function approach
Although it is difficult to solve the above optimization problem analytically, we can develop a computationally efficientalgorithm for finding a locally optimal solution based on theauxiliary function concept. Here we introduce the generalprinciple of the auxiliary function approach (also called the
majorization-minimization approach).
Let us use D(θ) to denote an objective function that wewant to minimize with respect to θ. G(θ, α) is defined as an
auxiliary function for D(θ) if it satisfies

Fig. 1. Redundancy of STFT spectrograms.
grams, we consider decomposing the observed time-domain
signal y into the sum of L signal components:
y=

L
∑

sl ,

(1)

l=1

such that the magnitude spectrogram of each component is
as close to a rank-1 structure as possible. Here, let us use
ψ k,m ∈ CN to denote arbitrary basis functions for timefrequency analysis, where k and m are the frequency and time
indices, respectively. See Fig. 2 for a graphical illustration.
With STFT, ψ k,m is a windowed complex sinusoid. By using ψ k,m , the magnitude spectrogram of sl can be written as|ψ Hk,m sl |. Thus, the problem of interest can be cast as the
optimization problem of minimizing
∑∑
2
(|ψ H
I(θ) =
k,m sl | − Hk,l Ul,m ) + R(U ),lsubject toD(θ) = min G(θ, α).
αWe call α an auxiliary variable. By using G(θ, α), D(θ) canbe iteratively decreased according to the following theorem:
Theorem 1. D(θ) is non-increasing under the updates, θ ←argminθ G(θ, α) and α ← argminα G(θ, α).
3.2. Designing auxiliary function
When applying the auxiliary function approach to a certain
minimization problem, the first step is to design an auxiliary
function that bounds the objective function from above. The
main difficulty with the present optimization problem lies inp
the discontinuity of the gradients of |ψ H
k,m sl | and |Ul,m | .
We can design an auxiliary function for I(θ) by invoking thefollowing two inequalities.

k,m∑sl = y,

(5)

(2)

l

with respect to θ = {H, U , S} where H = {Hk,l }, U ={Ul,m } and S = {sl }. Note that Hk,l ≥ 0 and Ul,m ≥ 0 are
analogous to the basis spectrum and the time-varying amplitude in the NMF model. The first term of the above objectivefunction becomes 0 when the magnitude spectrograms of all
the members of S have exactly rank-1 structures. The second
term R(U ) is a regularization term for U . It is important tonote that as with complex NMF, this model allows the components to cancel each other out, and so some constraint isneeded to induce the sparsity of U . For this purpose, we define R(U ) using the `p norm∑
R(U ) = 2λ|Ul,m |p ,(3)

Lemma 1. For any complex number z and any complex number c satisfying |c| = 1, we have−|z| ≤ −Re(c∗ z).

(6)

Equality holds when c = z/|z|.
Lemma 2. When 0 < p < 2, for any real (or complex) number x, we have2|x|p ≤ p|v|p−2 |x|2 + 2 − p|v|p .

l,mEquality holds when v = x.

where λ > 0 weighs the importance of the sparsity cost relative to the fitting cost. When 0 < p < 2, R(U ) promotes
sparsity if the norm of U is bounded. To bound U , we assume
∑
2
Hk,l
= 1.
(4)

By aplying (6) and (7) to I(θ), we obtain∑∑
2
I(θ) ≤|ψ H
k,m sl − Hk,l Ul,m cl,k,m |
l

k

87

k,m

(7)

+λ

∑(

)
2
p|Vl,m |p−2 Ul,m
+ (2 − p)|Vl,m |p .

(8)

..
.

l,m

...

The right-hand side of this inequality can be used as an auxiliary function for I(θ). Here, C = {cl,k,m } and V = {Vl,m }are the auxiliary variables. This auxiliary function is minimized with respect to C and V (equality of the above inequality holds) whenH
cl,k,m = ψ Hk,m sl /|ψ k,m sl |,Vl,m = Ul,m .

...

(9)
(10)

...

It should be noted that (9) indicates the phase spectrogram ofthe signal estimate sl , updated at the previous iteration.
As the update rule for the auxiliary variables is given by(9) and (10), all we need is to derive the update rule for θ.

..
.

∑

3.3. Update rulesHm ψ k,m ψ k,m can be set at a
circulant matrix in both cases of STFT and wavelet transform.

We can use the method of Lagrange multipliers to derive theupdate equations for S and H. By setting the partial derivative of the Lagrangian
)
(∑
∑
2
T
|ψ H
s
−
H
U
c
|
s
−
y
,
+
γ
l
k,l
l,m
l,k,m
l
k,m
s

By setting the partial derivatives of the auxiliary function
with respect to Hk,l and Ul,m at zeros, the update rules for H
and U are given as

Fig. 3. Graphical illustration of Ψ.

l,k,m

∑

l

with respect to sl at zero, we obtain sl = Ψ
where
∑∑
Ψ=
Re(ψ k,m ψ Hk,m ),kdl =(dl − γ s ),Re(Hk,l Ul,m cl,k,m ψ k,m ).

Ul,m |ψ H
k,m sl |Hk,l = √∑ m
∑
k(

∑

H
2
m Ul,m |ψ k,m sl |),Hk Hk,l |ψ k,m sl |.
2p−2
k Hk,l + λp|Vl,m |Ul,m = ∑

(11)

m

∑∑
k

−1

(15)

(16)

The iterative algorithm is summarized as follows:

(12)

m

1. Initialize H, U and S.
2. Update C and V using (9) and (10).
3. Update H, U and S using (14)–(16) and return to 2.

(12) amounts to performing time-frequency synthesis to obtain the l-th signal component, using the current estimates of
the time-frequency coefficients, Hk,l Ul,m cl,k,m . This process
corresponds to the inverse STFT when ψ k,m are defined as
the STFT basis functions. By substituting this result into (1)and solving λs , we can write the update equation for sl as:
{
∑ )}
1(
dl .
(13)
Ψy −
sl = Ψ−1 dl +L

For the initialization (step 1), we can use conventional NMF
algorithms followed by Wiener filtering to obtain the estimates of H, U and S.
3.4. Relation to complex NMF algorithm
It is interesting to compare the update rules of the present algorithm with those of the complex NMF algorithm [3]. The
aim with Complex NMF is to approximate an observed complex spectrogram Yk,m with a model of the form∑
(17)
Fk,m =
Hk,m Ul,m ejφl,k,m ,l

It may appear that Ψ ∈ RN ×N must be inverted to compute (13). However, this can be avoided under the following
setting. ∑First, the time-frequency basis function can be set
so that m ψ k,m ψ Hk,m becomes a circulant matrix (see Fig.
3). Thus, Ψ can be diagonalized using the discrete Fouriertransform matrix: Ψ = F F H ΨF F H = F DF H , where
F denotes the discrete Fourier tranform matrix and D denotes a diagonal matrix whose diagonal elements are given as
the sum of the power spectra of ψ 1,m , . . . , ψ K,m . It can beshown (without proof, owing to space limitations), for example, that with STFT using the square-root Hanning window,
D becomes an identity matrix, in which case Ψ also becomes
an identity matrix. Thus, (13) can be simplified to∑ )1(
sl = dl +
y−
dl .
(14)
L

l

where φl,k,m denotes the phase spectrogram of the l-th signal component. Under a certain condition, the complex NMF
algorithm can be described as follows:
1. Initialize H, U and φ.
2. Update X = {Xl,k,m } and V = {Vl,m } using(18)
Xl,k,m =Hk,l Ul,m ejφl,k,m
(
)
∑
1
jφ
+ L Yk,m − l0 Hk,l0 Ul0 ,m e l0 ,k,m ,Vl,m =Ul,m .
(19)

l

88

3. Update H, U and φ using∑m Ul,m |Xl,k,m |Hk,l = √∑ ∑,2k(m Ul,m |Xl,k,m |)∑Hk,l |Xl,k,m |Ul,m = ∑ k2
,
p−2
H
k
k,l + λp|Vl,m |ejφl,k,m = Xl,k,m /|Xl,k,m |.

Table 1. SNR improvements (dB) obtained by NMF and TSF.
(20)

Track #
1
2
3
4
5

(21)
(22)

and return to 2.
Interestingly, there is a similarity between the present andcomplex NMF algorithms, even though the models and theobjective functions are different. Specifically, Xl,k,m andejφl,k,m of the complex NMF algorithm are analogous to sl
and cl,k,m of the present algorithm, respectively.
Xl,k,m can be viewed as an estimate of the complex
spectrogram of the l-th signal component. In step 2, Xl,k,m
is updated by adding the portion of the error between theobserved spectrogram and the model to the current estimate of Hk,l Ul,m ejφl,k,m . φ is then updated at its argumentarg(Xl,k,m ), and H and U are updated using its magnitude|Xl,k,m |. Similarly, according to (14), sl is updated by adding
the portion of the error between an observed signal and thesum of dl to the current estimate of dl , where dl is the signal
converted from the set {Hk,l Ul,m cl,k,m }k,m . According to(9), (15) and (16), cl,k,m is then updated at the argument ofψHk,m sl (the complex spectrogram of sl ), and H and U are
updated using the magnitude of ψ H
k,m sl .
One drawback with complex NMF is that Xl,k,m is notguaranteed to satisfy the explicit condition that complex spectrograms must satisfy. This implies that the complex NMF
algorithm is designed to search an unnessarily large solution
space for the optimal parameters. By contrast, the present
algorithm always ensures that the estimate of the complex
spectrogram of the l-th latent component, ψ H
k,m sl , is associated with a time-domain signal sl , keeping the search withina proper solution space.

NMF
8.53
4.96
12.27
10.41
10.19

TSF
11.03
8.20
13.05
11.23
11.53

∑
subject to r,l sr,l = y. Here, each r indicates a different time-frequency resolution. The above formulation can be
explained as follows. All the signal components are divided
into R groups, each of which consists of Lr members, andsr,l denotes the l-th signal component within the group r.
Since ψ r,k,m is indexed by r, the spectrograms of the signal components in a different group r0 will have a different
time-frequency resolution.
It should be noted that the parameter update rules are derived in the same way as in the previous section.
5. EXPERIMENTS
We quantitatively compared the source separation performance of NMF (followed by Wiener filtering) and the proposed method (hereafter, time-domain spectrogram factorization (TSF)) by conducting supervised source separation experiments. We used professionally produced music recordings from the SiSEC 2013 database, available athttps://sisec.wiki.irisa.fr/, as the experimental data. Each recording is a mixture of multiple tracks, eachof which is produced by a single instrument or singer. The
separated tracks are also available. We divided each recording
into two segments, namely a test data segment and a training
data segment. With all these methods, the basis spectra werepretrained using the individual tracks of the training data, andthen source separation was performed on the test data. All
the audio samples were monaural and sampled at 22.05kHz.
STFT was computed using a square-root Hanning windowthat was 32ms long with a 16ms overlap. With both methods, 6 basis spectra were assigned to each track. Thus, for
5-track recordings, a total of 30 basis spectra were used forthe separation. Tab. 1 shows the signal-to-noise Ratio (SNR)
improvements after the separations with the two methods.
From these results, we confirmed that TSF performed better
than NMF.

4. MULTI-RESOLUTION DECOMPOSITION
Another important advantage of the present model is that it
enables multi-resolution signal decompositions.
The optimal time-frequency resolution may depend on thekinds of sound sources. For example, temporal resolution
is more important than frequency resolution for percussive
sounds, with the opposite being true for low-pitched sounds.
Thus, it is desirable to be able to arbitrarily set different timefrequency resolutions for individual sources. Of course, thiswas not possible with the NMF framework and its variants
including complex NMF and IS-NMF, since in the conventional frameworks, the decompositions are carried out in thetime-frequency domain with a particular resolution. By contrast, the proposed framework allows us to model the timefrequency structure of each component with a different resolution, thanks to the time-domain formulation.
The multi-resolution extension is straightforward. By additionally introducing a “resolution index” r, we obtain an
objective function for the multi-resolution signal decomposition problem:
∑∑
2
(|ψ H
I(θ) =
r,k,m sr,l | − Hr,k,l Ur,l,m ) + R(U ),6. CONCLUSIONSThis paper proposed time-domain spectrogram factorization(TSF), which makes it possible to realize NMF-like signal decompositions in the time domain. This new formulation also
allows for an extension to multi-resolution signal separation,
which was not possible with the conventional NMF framework. We confirmed through supervised source separation
experiments that the proposed method outperformed NMF.

r,l k,m89

7. REFERENCES[1] P. Smaragdis and J. C. Brown, “Non-negative matrix factorization for music transcription,” in Proceedings of theIEEE Workshop on Applications of Signal Processing to
Audio and Acoustics (WASPAA), 2003, pp. 177–180.
[2] T. Virtanen, “Sound source separation using sparse coding with temporal continuity objective,” in Proceedings of the International Computer Music Conference(ICMC), 2003, pp. 231–234.
[3] H. Kameoka, N. Ono, K. Kashino, and S. Sagayama,“Complex NMF: A new sparse representation for acoustic signals,” in Proceedings of 2009 IEEE International
Conference on Acoustics, Speech and Signal Processing(ICASSP2009), 2009, pp. 3437–3440.
[4] R. M. Parry and I. Essa, “Phase-aware non-negative spectrogram factorization,” in Proceedings of the International Conference on Independent Component Analysis
and Signal Separation (ICA), 2007, pp. 536–543.
[5] C. Févotte, N. Bertin, and J. -L. Durrieu, “Nonnegative
matrix factorization with the Itakura-Saito divergence.
With application to music analysis,” Neural Computation, vol. 21, no. 3, 2009.
[6] J. Le Roux, H. Kameoka, N. Ono, and S. Sagayama,“Fast signal reconstruction from magnitude stft spectrogram based on spectrogram consistency,” in Proceedingsof The 13th International Conference on Digital Audio
Effects (DAFx-10), 2010, pp. 397–403.

90