2-1-5
Experimental Evaluation of Postfilter-Based Nonnegative Matrix Factorization
with Statistical Model Parameter Estimation ∗☆ Murota Yuki, Daichi Kitamura, Shunsuke Nakai, Hiroshi Saruwatari,
Satoshi Nakamura (Nara Institute of Science and Technology)
Kazunobu Kondo, Yu Takahashi (Yamaha Corporation)

1 Introduction

STFT

In recent years, music signal separation has been avery active area of signal-processing research. This technique is suitable for many potential applications, e.g.,controlling each source in a music signal in interactive
3D audio systems and realizing automatic music transcription for each instrument player.
Common methods used for audio signal separation,
which were mainly developed for speech enhancement,
are nonlinear filtering algorithms such as Wiener filtering (WF) [1] and the minimum mean-square error shorttime spectral amplitude (MMSE-STSA) estimator [2].
In particular, the MMSE-STSA estimator and its extended algorithms are optimal Bayesian estimators basedon the a priori speech statistical model, resulting in agreat improvement of sound quality. However, thesemethods cannot be eﬀectively applied to music signals
because it is diﬃcult to deal with nonstationary interference signals. Also, the performance has strong dependence on the selection of the a priori statistical model,
which cannot be determined automatically.
As another means of signal separation, nonnegative
matrix factorization (NMF) has been actively studied [3],in which an input spectrogram is decomposed into the
product of a spectral basis matrix and its activation matrix. In particular, supervised NMF (SNMF) [4],which
includes a priori training with some sample sounds of atarget instrument, can extract the target signal to someextent. Although SNMF can deal with nonstationary
source signals, it has an inherent drawback that the decomposition includes approximation only valid for linear
combinations of spectrograms, and is not strictly valid
for time- (or time-frequency-) domain complex-valued
mixtures. In addition, it is diﬃcult to implement the statistical model for time sequences.
Motivated by the complementarity between the properties of the MMSE-STSA estimator and SNMF, in thispaper, we propose a new approach based on the generalized Bayesian estimator with automatic prior estimation suitable for music signal separation. This method
consists of three parts, namely, the generalized MMSESTSA estimator [5] with a flexible target signal prior, the
SNMF-based interference spectrogram estimator, and a
new closed-form parameter estimation for the statistical model of the target signal time sequence based onhigher-order statistics.
Compared with the conventional methods, the proposed method has the following advantages: (I) The target signal extraction is carried out via the generalized
MMSE-STSA estimator so that the mixing of the timefrequency-domain complex-valued signals can be properly considered without any approximation. (II) Thanks
to the SNMF-based spectrogram estimator, we can dynamically estimate the nonstationary power spectra of
the interference signal. (III) The statistical model of thehidden target signal can be detected automatically only∗~S ∗ ( f ,τ )Y∗ ( f ,τ )Y∗ ( f ,τ )G ( f ,τ )Shape parameter
estimation
kurt target
Target kurtosis
estimation

ξ ( f ,τ )A posteriori
SNR estimation

PN~ ( f )Fig. 1~γ~ ( f ,τ )PN~ ( f )Interference signal
estimation by SNMF

ISTFT

G ( f ,τ )Spectral gain
calculation

A priori
SNR estimation

γ~ ( f ,τ )Block diagram of proposed method.

using the observable data, and can be used for optimal
Bayesian estimation with online target-signal prior adaptation.

2 Conventional method
2.1 Generalized MMSE-STSA estimator [5]We apply short-time Fourier analysis to the observed
signal, which is a mixture of target and interference
signals, to obtain the time-frequency-domain complexvalued signal
YR ( f, τ)+iYI ( f, τ) = (S R ( f, τ)+iS I ( f, τ))+(NR ( f, τ)+iNI ( f, τ)),(1)
where Y∗ ( f, τ) is the observed signal, S ∗ ( f, τ) is thetargetN∗ ( f, τ) is the interference signal, and{ signal,}∗= R, I denote the real and imaginary parts of the signal, respectively. Also, f is the frequency bin and τ isthe frame index. The MMSE-STSA estimator provides
the amplitude spectrum of the target signal based on the
MMSE criterion under the assumption that the amplitude
spectrum of the target signal obeys a Rayleigh distribution. The generalized MMSE-STSA estimator [5] provides the amplitude spectrum of the target signal under
the assumption of a chi distribution.
2.2 SNMF [4]The mixture model of NMF approximately assumes
the additivity of amplitude (or power) spectrograms as√√√
YR2 ( f, τ)+YI2 ( f, τ) ≈ S R2 ( f, τ)+S I2 ( f, τ)+ NR2 ( f, τ)+NI2 ( f, τ) .

(2)
SNMF [4]consists of two processes, namely, training of
the target sound and separation of the observed signal. In
the training process, the supervised bases are trained as adictionary of the target sound. Then, the observed spectrogram is decomposed into the target and other spectrograms using the supervised bases in the separation process.

3 Proposed method
3.1 Generalized MMSE-STSA estimator
Figure 1 depicts a block diagram of the proposed
method. In the next part, we describe these processes
in detail.

“統計モデルパラメータ推定を用いたポストフィルタに基づく非負値行列因子分解の実験的評価,” by 室田勇騎,
北村大地, 中井駿介, 猿渡洋, 中村哲 (奈良先端大), 近藤多伸, 高橋祐 (ヤマハ株式会社).

日本音響学会講演論文集- 677 -

2014年3月

In the generalized MMSE-STSA estimator, the amplitude spectrum of the target signal is estimated on thebasis of the MMSE criterion under a certain target prior.
The processed signal S̃ ( f, τ) via the generalized MMSESTSA estimator is given byS˜∗ ( f, τ) = G( f, τ)Y∗ ( f, τ),(3)√
ν( f, τ) Γ(ρ+0.5) Φ(0.5−ρ, 1, −ν( f, τ))G( f, τ) =··,γ( f, τ)Γ(ρ)Φ(1−ρ, 1, −ν( f, τ))(4)where Γ(·) is the gamma function, Φ(a, b; k) = F 1 (a, b; k)is the confluent hypergeometric function, and()−1
ν( f, τ) = γ̃( f, τ)ξ̃( f, τ) 1 + ξ̃( f, τ) .

(5)

Here, ξ̃( f, τ) and γ̃( f, τ) are the estimated a priori and aposteriori SNRs, respectively, which are defined asξ̃( f, τ) = αγ̃( f, τ − 1)G2 ( f, τ) + (1 − α)max[γ( f, τ) − 1, 0],(6)γ̃( f, τ) =(YR2 + YI2 )/PÑ ( f ),(7)

where PÑ ( f ) is the estimated interference power spectral
density and α is the forgetting factor.
In the generalized MMSE-STSA estimator, the a priori statistical model of the target signal amplitude spectrum is set to
(
)ρ
(
)
ρ
ρ
2
2ρ−1
2
x
exp
x
,
(8)
p(x) =
Γ(ρ) E[x2 ]E [x2 ]
where p(x) is the p.d.f. of signal x in the amplitude domain and ρ is the shape parameter. Here, ρ = 1 gives
a Rayleigh distribution that corresponds to a Gaussian
distribution in the time domain, and a smaller value of ρ
corresponds to a super-Gaussian distribution signal.
In the generalized MMSE-STSA estimator, to calculate γ̃( f, τ), dynamic estimation is required if the interference signal is nonstationary, and estimation of the
shape parameter ρ, which depends on the type of target signal, is also required. To solve these problems, we
propose the use of SNMF as the interference signal estimator and estimate the shape parameter ρ using higherorder statistics of the target signal.
3.2 Interference signal estimation by SNMF
The following equation represents the decomposition
model of SNMF using the trained supervision components F( f, k):√A( f, τ) = YR2 ( f, τ) + YI2 ( f, τ)∑∑≈ k F( f, k)V(k, τ) + n H( f, n)U(n, τ), (9)where F( f, k) is a nonnegative element of the supervised basis matrix trained in advance, which involves
spectral patterns of the target signal as column vectors,
V( f, k) is a nonnegative element of an activation matrix
that corresponds to F( f, k), H( f, n) represents a nonnegative element of the other basis matrix, which involves∑ residual spectral patterns that cannot be expressedby k F( f, k)V(k, τ), and U(n, τ) is a nonnegative element of the activation matrix that corresponds to H( f, n).
Moreover, k is the basis index of F( f, k), and n is the basis index of H( f, n). The supervised basis matrix can betrained using sample sounds of the target signal in the

日本音響学会講演論文集∑
training process. Hence, ideally, k F(∑f, k)V(k, τ) represents the target signal components and n H( f, n)U(n, τ)
represents the other components diﬀerent from the target
signals after the decomposition.
∑A( f, τ) −n H( f, n)U(n, τ) (or∑ We can usek F( f, k)V(k, τ)) as the estimated amplitude spectrogram of the interference signal for∑ the generalized
MMSE-STSA estimator. Thus, ( n H( f, n)U(n, τ))2
is regarded as a good estimate of PÑ ( f ) in (7) in thetime-frequency grids even if the interference sounds are
nonstationary, which is common in actual music signals.
3.3 Target signal prior estimation
3.3.1 Shape parameter and kurtosis
Generally, we cannot obtain any a priori statisticalmodel (8) from the training data (e.g., several octave
notes of the target instrument) in SNMF because the statistical time structure is quite diﬀerent from that of thetargetsignal S ∗ ( f, τ). Also, the target signal component∑k F( f, k)V(k, τ) in SNMF cannot be used because its
accuracy is not sometimes enough. Therefore, we inversely calculate the parameter of the target amplitude
spectrogram in a data-driven manner, utilizing two observable statistics of the input signal and interference
spectrogram estimated by SNMF.
Regarding the chi distribution p(x) in (8), the shapeparameter ρ can be written as follows:ρ = (µ4 /µ22 − 1)−1 ,(10)

where µ4 /µ22 is called the kurtosis and µm is the mthorder moment of the amplitude spectrum. µm is defined
as
∫ ∞
µm =
xm p(x)dx.
(11)0

From this relation, the shape parameter of the subjective
target signal can be estimated by obtaining its kurtosis
value. In general, however, it is diﬃcult to directly estimate the kurtosis of such a hidden target signal becauseof the contamination by additive interference signals. In
the following subsections, a new algorithm of target kurtosis estimation is proposed for estimating of the shapeparameter of the target p.d.f.
To cope with the mathematical problem that the mixing of the target and interference signals is additive but
generally their higher-order moments are not additive,
we introduce the cumulant, which holds the additivity for
additive variables. Meanwhile, in transformation froma waveform to its amplitude spectrum, the exponentiation operation is conducted but the cumulant does nothave a straightforward relationship. In this case, we use
the moment instead of the cumulant. Thus, we propose
moment-cumulant transformation.
3.3.2 Moment-cumulant transformation
In this section, we derive some formulas regarding
moment-cumulant transformation. They explicitly represent the relations between the moment and cumulant
in each order, which are useful for estimating the kurtosis of the target amplitude spectrum.
First, the characteristic function ϕ x (it) of the randomvariable x is defined as∫ ∞ϕ x (it) =eitx P(x)dx.
(12)−∞Then, we can define the mth-order moment µm (x) and- 678 -

2014年3月

the mth-order cumulant κm (x) of x as follows:µm (x) =∂(m) ϕ x (it),∂it(m) t=0(13)κm (x) =∂(m) log ϕ x (it).
∂it(m)
t=0(14)

In the same manner, we can estimate the statistics of thesquared variable of S I . Given µm (S R2 ) and µm (S I2 ), we
can calculate the cumulant of the power spectrum S R2 +S I2
asκm (S R2 + S I2 )= κm (S R2 ) + κm (S I2 )∑∏= π(m) (−1)|π(m)|−1 (|π(m)| − 1)! B∈π(m) µ|B| (S R2 )∑∏+ π(m) (−1)|π(m)|−1 (|π(m)| − 1)! B∈π(m) µ|B| (S I2 ),(24)

Next, polynomial forms of interrelations between the
moment and cumulant are derived below. From (13), themth-order moment µm (x) can be rewritten as∂(m) exp(log ϕ x (it))∂it(m)t=0∑ ∏= π(m) B∈π(m) κ|B| (x),and the mth-order moment of the power spectrum isgiven byµm (x) =where we use a combinational form of Faà di Bruno’s
formula,∑∏∂(m) f (g(x))= π(m) f (|π(m)|) (g(x)) B∈π(m) [g(x))](|B|) ,
∂x(m)

=

|π(m)|−1(16)∏(|π(m)| − 1)! B∈π(m) µ|B| (x).
(17)

3.3.3 Estimation of target signal kurtosis from observations
In this section, we estimate the amplitude-domain kurtosis of the target signal. First, we express the target
kurtosis by using complex-valued variables as intrinsic
parameters. After that, we show that the kurtosis can be
represented in the amplitude-spectrogram domain. Hereafter, we ignore the indexes f, τ, k, and n of each signal.
Only the statistics of (YR +iYI ) and (NR2 + NI2 )1/2 obtainedby SNMF are observable, and (S R2 + S I2 )1/2 is a hidden
value to be estimated. First, we assume to obtain thefollowing mth-order moments from the data:µm (YR ) = E[YRm ],µm (YI ) = E[YIm ],µm (NR ) = E[NRm ],µm (NI ) = E[NIm ].

(18)
(19)
(20)
(21)

Generally, the cumulant has the property of additivity for
additive independent variables, i.e., κm (a + b) = κm (a) +κm (b). Using this relation and (17), we can estimate thecumulant of the real part of the target signal S R =YR − NRasκm (S R ) = κm (YR ) − κm (NR )∏∑= π(m) (−1)|π(m)|−1 (|π(m)| − 1)! B∈π(m) µ|B| (YR )∏∑− π(m) (−1)|π(m)|−1 (|π(m)| − 1)! B∈π(m) µ|B| (NR ).
(22)

日本音響学会講演論文集

∑

π(2m)

∏

B∈π(2m) κ|B| (S R ).

2
2
B∈π(m) κ|B| (S R + S I ).

(25)

(26)

Using (18)–(26), we can estimate the resultant kurtosisof the target amplitude spectrum as
1kurttarget =µ4 ((S R2 + S I2 ) 2 )
1µ22 ((S R2 + S I2 ) 2 )N (µm (YR ), µm (YI ), µm (NR ), µm (NI )), (27)=
D (µm (YR ), µm (YI ), µm (NR ), µm (NI ))where
N (µm (YR ), µm (YI ), µm (NR ), µm (NI ))= µ4 (YR ) + µ4 (YI ) − µ4 (NR ) − µ4 (NI )+ 6µ22 (NR ) + 6µ22 (NI ) + 2µ2 (YR )µ2 (YI ) + 2µ2 (NR )µ2 (NI )− 6µ2 (YR )µ2 (NR ) − 6µ2 (YI )µ2 (NI )− 2µ2 (YR )µ2 (NI ) − 2µ2 (YI )µ2 (NR ),(28)
D (µm (YR ), µm (YI ), µm (NR ), µm (NI ))= µ22 (YR ) + µ22 (YI ) + µ22 (NR ) + µ22 (NI ) + 2µ2 (YR )µ2 (YI )− 2µ2 (YR )µ2 (NR ) − 2µ2 (YR )µ2 (NI ) − 2µ2 (YI )µ2 (NR )− 2µ2 (YI )µ2 (NI ) + 2µ2 (NR )µ2 (NI ).
(29)
Next, the sums of the 4th-order moments µ4 (YR ) +µ4 (YI ) and µ4 (NR ) + µ4 (NI ) are represented by the
amplitude-domain kurtosis of the observed signal spectrum and the interference signal spectrum as(µ4 (YR ) + µ4 (YI ) = µ22 (YR ) + µ22 (YI ) + 2µ2 (YR )µ2 (YI )) µ4 (A)µ22 (A)−2µ2 (YR )µ2 (YI ),(30)−2µ2 (NR )µ2 (NI ).

(31)) µ (∑ HU )µ4 (NR ) + µ4 (NI ) = µ22 (NR ) + µ22 (NI ) + 2µ2 (NR )µ2 (NI ) 24 (∑n )µ2 n HU(

If we assume that the real and imaginary parts are i.i.d.,µ2 (YR ) equals µ2 (YI ) and µ2 (NR ) equals µ2 (NI ). Under
this assumption and (15) and (17), we obtain the following relation for the SNMF output:

The statistics of the squared variable of S R is given byµm (S R2 ) = µ2m (S R ) =∏µm ((S R2 + S I2 ) 2 ) = µ m2 (S R2 + S I2 ).

∏(|π(m)|)(ϕ x (it)) B∈π(m) [ϕ x (it))](|B|) t=0π(m) log

π(m) (−1)π(m)

Furthermore, the mth-order moment of the amplitude
spectrum (S R2 + S I2 )1/2 is

∑
∑

∑

1

where π(m) runs through the list of all partitions of a setof size m, B ∈ π(m) means that B is one of the blocks
into which the set is partitioned, and |B| is the size of theset B.
In the same manner, from (14), the mth-order cumulant κm (x) is given byκm (x) =µm (S R2 + S I2 ) =

(15)

(23)

- 679 -1µ2 (YR ) = µ2 (YI ) = µ2 (A),21 ∑µ2 (NR ) = µ2 (NI ) = µ2 ( n HU).
2

(32)
(33)

2014年3月

12

(a)

11

11

10

10

SDR [dB]

SDR [dB]

12

WF + SNMF

98
7
6

Fig. 2

MMSE-STSA + SNMF

16

(b)




(c)

14

9
8

12
10

7

8

6

6

Average SDRs for extraction of (a) clarinet signal, (b) oboe signal, and (c) cello signal for each method.

method is better than those of the other methods. This
result indicates the eﬃcacy of introducing the flexible apriori statistical model of the target signal. The simple
MMSE-STSA estimator also assumes the fixed a priori
model of the Gaussian distribution but the assumption is
not appropriate for a music target signal. In contrast, theproposed method almost always uses more spiky p.d.f.
(ρ < 1), which enhances the true sparseness of the target
music signal.

1

 

Proposed method

SDR [dB]

SNMF¶ DE? E E E DE E E EE E EE E DD¶
´EEEEE E D
Fig. 3

Scores of each part.

Finally, using (30)–(33), (27) is rewritten as follows:∑∑∑µ (A)−µ ( n HU )+4µ22 ( n HU)−4µ2 (A)µ2 ( n HU)
kurttarget = 4 µ42 (A)+µ.
2 (∑ HU)−2µ (A)µ (∑ HU)
2
2
n
n
2
2

(34)
The shape parameter of the target p.d.f. can be estimated using (10) and (34). Note that all the estimates
can be obtained from the result of SNMF without using
any waveforms; this implies a good applicability to the
combination with SNMF.

4 Experiments
4.1 Experimental conditions
In this experiment, we compared four methods, i.e.,
simple SNMF [4], WF with SNMF-based interference
estimation (WF+SNMF) [6], the MMSE-STSA estimator with SNMF-based interference estimation (MMSESTSA+SNMF), and the proposed method. WF and the
MMSE-STSA estimator utilized the interference spectrogram estimated by SNMF. We used three instrumental signals, namely, an oboe, clarinet, and cello, as thetarget sounds (each melody part is depicted in Fig. 3).
These signals were artificially generated by a MIDI synthesizer, and the observed signals were produced by mixing two sources selected from three signals with the same
power. In estimation of the interference signal using
SNMF, we used the same MIDI sounds of the target instruments as supervision for the training process. The
training sounds contained two octave notes that cover all
the notes of the target signal in the observed signal. The
spectrograms were computed using a 92-ms-long rectangular window with a 46-ms overlap shift. Moreover,the number of trained bases was 100 and the number ofother bases was 50. In the proposed method, the forgetting factor α and amplitude compression parameter βwere set to 0.1 and 1.0, respectively.
4.2 Experimental results and discussion
We used the signal-to-distortion ratio (SDR) definedin [7] as the evaluation score. SDR indicates the overall quality of the separated target sound, showing highseparation and less artificial distortion.
Figure 2 shows the average SDRs for each method andeach target instrument. From the SDR results, we canconfirm that the separation performance of the proposed

日本音響学会講演論文集

5 Conclusions
In this paper, we propose a new approach for addressing music signal separation based on the generalized Bayesian estimator with automatic prior adaptation. From the experimental evaluation, it is found thatthe proposed method outperforms competitive methods,
namely, simple NMF, WF, and the MMSE-STSA estimator with a fixed Gaussian prior.

References[1] P. C. Loizou, “Speech Enhancement Theory andPractice” CRC Press, Taylor & Francis Group,FL, 2007.
[2] Y. Ephraim, D. Malah, “Speech enhancement using a minimum mean-square error short-time spectral amplitude estimator,” IEEE Trans. Acoustics, Speech and Signal Processing, vol.32, no.6,
pp.1109–1121, 1984.
[3] D. D. Lee, H. S. Seung, “Algorithms for nonnegative matrix factorization,” Proc. Advances in
Neural Information Processing Systems, vol.13,
pp.556–562, 2001. 2001.
[4] P. Smaragdis, B. Raj, M. Shashanka, “Supervisedand semi-supervised separation of sounds fromsingle-channel mixtures,” Proc. 7th International
Conference on Independent Component Analysis
and Signal Separation, pp.414–421, 2007.
[5] C. Breithaupt, M. Krawczyk, R. Martin, “Parameterized MMSE spectral magnitude estimation for
the enhancement of noisy speech,” Proc. IEEE International Conference on Acoustics, Speech andSignal Processing, pp.4037–4040, 2008.
[6] E. M. Grais, H. Erdogan, “Single channel speechmusic separation using nonnegative matrix factorization and spectral masks,” Proc. IEEE 17th International Conference on Digital Signal Processing,
pp.1–6, 2011[7] E. Vincent, R. Gribonval, C. Fevotte, “Performance measurement in blind audio source separation,” IEEE Trans. Audio, Speech and Language
Processing, vol.14, no.4, pp.1462–1469, 2006.

- 680 -

2014年3月