ROBUST SPEECH DEREVERBERATION BASED ON NON-NEGATIVITY
AND SPARSE NATURE OF SPEECH SPECTROGRAMS
Hirokazu Kameoka, Tomohiro Nakatani, Takuya Yoshioka,
NTT Communication Science Laboratories, NTT Corporation,
3-1 Morinosato Wakamiya, Atsugi, Kanagawa 243-0198, Japana model that implicitly assumes the additivity of power spectra, which holds only approximately, the performance may
be limited to some extent. However, we hypothesize that it
is advantageous in one sense. That is, the phase characteristic of a room impulse response is a characteristic that is especially sensitive to speaker movements. Therefore, the useof a convolution model of the subband envelopes, in which
only the time-invariance of the amplitude characteristic is assumed, may allow the system to be robust against the changes
in room characteristics.
On the basis of this hypothesis, we choose the latter approach in this paper. While most of the temporal envelope
ﬁltering techniques require strong assumptions about MTF(some require MTF to be measured prior to the dereverberation stage, while others employ a theoretically derived MTF
model with only a few parameters such as the reverberation
time RT60 ), the motivation for this work has been to develop
a blind temporal envelope ﬁltering method that only assumes
the non-negativity of the subband envelopes of speech and aroom impulse response, and the sparseness of speech. The
problem is thus formulated as a blind deconvolution problem
with non-negativity constraints, regularized by a sparsity cost.
For its optimization we derive an iterative algorithm that ensures a monotonic decrease in the objective function and thenon-negativity of the parameters.

ABSTRACTThis paper presents a blind dereverberation method designed
to recover the subband envelope of an original speech signal
from its reverberant version. The problem is formulated asa blind deconvolution problem with non-negative constraints,
regularized by the sparse nature of speech spectrograms. We
derive an iterative algorithm for its optimization, which canbe seen as a special case of the non-negative matrix factor
deconvolution. We conﬁrmed through experiments that thealgorithm is fast and robust to speaker movement.
Index Terms— Speech dereverberation, temporal envelope ﬁltering, non-negative blind deconvolution, sparseness
1. INTRODUCTION
Reverberation in a room severely degrades the characteristics
and quality of speech captured by distant microphones, thus
posing a severe problem for many speech applications. Theaim of speech dereverberation techniques is to recover a clean
speech signal from its reverberant version. When only a reverberant speech signal is accessible, dereverberation systems
must work ‘blind’.
A number of blind dereverberation systems have alreadybeen developed. While there are several viable solutions
in situations where multiple channels are available, single
channel systems still pose a formidable challenge. For single
channel systems two main approaches have been adopted.
One approach involves applying inverse ﬁltering techniques.
This approach requires assumptions about clean speech (e.g.
harmonicity, sparseness, non-stationarity, non-Gaussianity)
and/or the use of a speech model (e.g. dual excitation model,autoregressive model, codebook) to assess ‘speech naturalness’. Systems are then able to estimate the optimal inverse
ﬁlter such that the speech naturalness of the inverse-ﬁlteredsignal is maximized [1]–[9]. While inverse ﬁltering methods are usually performed in the time domain, the short-time
Fourier transform (STFT)-domain approach is more computationally efﬁcient [8, 9]. Within this framework, sinceroom impulse responses can vary rapidly according to suchfactors as the speaker’s position, the problem of how robustly
and adaptively a system can estimate an inverse ﬁlter from ashort-term observation has attracted particular attention.
With a different approach, an attempt is made to recover
the subband envelope (power envelope of a subband signal) of
the original speech by applying an inverse ﬁlter of the modulation transfer function (MTF) [10]–[17]. As it is based on

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

2. PRINCIPLE
2.1. Subband envelope model of reverberant speech
Let sk [t], hk [t] be the subband signals of speech and a room
impulse response at the kth subband where t is the time index.
The reverberant subband signal is then approximated by
xk [t] =sk [τ ]hk [t − τ ],(1)
τparticularly where each subband signal is obtained by STFT[8].
Then, the subband envelope can be written asxk [t]2 =s∗k [t − τ ]h∗k [τ ]sk [t − τ  ]hk [τ  ](2)=ττ

s∗k [t − τ ]sk [t − τ  ]hk [τ ]hk [τ  ]e−jφk [τ ] ejφk [τ ] ,τ,τ 

where hk [τ ] = |hk [τ ]|ejφk [τ ] . As hypohesized in 1, the phaseφk [τ ] can vary sensitively with respect to the reverberant conditions. We shall thus ﬁnd it convenient to treat the phase as a45

ICASSP 2009

random variable and integrate it out of the model. Hence, byassuming that φk [t] is an independent random variable uniformly distributed on the interval D = [−π, π), the expectation of |xk [t]|2 leads to  
sk [t − τ ]2 hk [τ ]2 ,(3)E |xk [t]|2 =τwhich suggests that the subband envelope of reverberant
speech can be represented, in an expectation sense, as theconvolution of subband envelopes of clean speech and theroom impulse response.
2.2. Problem setting
For simplicity of notation, let Sk [t] ≡ |sk [t]|2 , Hk [t] ≡|hk [t]|2 . According to the discussion in 2.1, we model the kth
subband envelope of reverberant speech asSk [τ ]Hk [t − τ ],(4)Xk [t] ≡τand assume t Hk [t] = 1 in order to avoid an indeterminacy
in the scaling. Given an observed subband envelope, Yk [t],the goal is to ﬁnd an approximation such that Yk [t]  Xk [t]in which Sk [t] is sparse, based on the sparse nature of thespeech spectrogram. It should be noted that the ‘sparseness ofspeech’ is the only assumption we make about clean speech.
We now assume the following generative model
Yk [t] = Xk [t] + k [t].

(5)

The reconstruction error k [t] is assumed to include any errors resulting from the approximations in 2.1. Assuming k [t]is Gaussian white noise that follows N (0, σ 2 ), the likelihood of S ≡ {Sk [1], · · · , Sk [T ]}K
k=1 and H ≡ {Hk [1], · · · ,,givenY≡{Y[1],···, Yk [T ]}K
Hk [T ]}Kkk=1k=1 , is written as 1(Yk [t] − Xk [t])2√. (6)exp −P (Y |S, H) =2σ 2
2πσk,tWe use a generalized Gaussian prior for P (S)|Sk [t]|p1
exp−P (S) =bp2Γ 1 + p1 bk,t,(10)t2.3. Multiplicative update algorithm
Guided by the idea of NMF[18], we derive an efﬁcient iterative algorithm that ensures a monotonic decrease (convergence to a stationary point) in the objective function and, simultaneously, the non-negativity of the parameters.
First we derive the update formula for S. Let S  and H 
be the parameters at the previous iteration such that Sk [t] ≥ 0and Hk [t] ≥ 0. We then have an inequality
f (S, H ) ≤+ S  [τ ]H  [t − τ ] 
kk,t,τ

pSk [t]p−2kXk [t]Yk,t −Sk [τ ] X [t]Sk [τ ] k

2

Sk [t]2 + 2|Sk [t]|p − p|Sk [t]|p , (11)k,twhere the equality holds when Sk [t] = Sk [t]. Its proof isomitted owing space limitations. We write the right-hand sideof the inequality as f¯(S). It can be proved that the minimization of f¯(S) w.r.t. Sk [t] leads to a certain decrease off (S, H  ). Thus, differentiating f¯(S) partially w.r.t. Sk [t] andsetting it at 0, we obtain the update formula, which is often
referred to as the ‘multiplicative update rule’, for Sk [t]Sk [τ ] = Sk [τ ] t

Hk [t − τ ]Yk [t]Hk [t − τ ]Xk [t] + λp|S  k [τ ]|p−1, (12)where(7)(8)Regarding p, σ 2 and b as constant parameters, we arrive at thefollowing objective function:
2
Yk [t] − Xk [t] + 2λ
|Sk [t]|p . (9)f (S, H) ≡k,tminimize f (S, H) with respect to S and H
subject to
Hk [t] = 1, Hk [t] ≥ 0, Sk [t] ≥ 0.

t

and assume for convenience that P (Hk [1], · · · , Hk [T ]) is independent over k and follows a uniform distribution (or moreaccurately, a Dirichlet distribution
 whose parameters are allset at 1, as H is constrained to t Hk [t] = 1). When 0 <p < 2, P (S) becomes super-Gaussian and promotes sparsity
if the norm of S is bounded. The likely values of S and H
can thus be inferred from the posterior density
P (S, H|Y ) ∝ P (Y |S, H)P (S).

λ, determined by p, σ 2 and b, weighs the importance of thesparsity cost relative to the accurate reconstruction. We notice
that, as S and H are non-negative in nature, we must seek tominimize f (S, H) subject to Sk [t] ≥ 0 and Hk [t] ≥ 0. We
are therefore led to solve the following optimization problem:

k,t46Xk [t] =τSk [τ ]Hk [t − τ ].

(13)

We note that Eq. (12) comprises the product and sum of nonnegative entities, and so the non-negativity of the parameter
update is thus guaranteed. Next, we derive the update formulafor H. Similarly, we have an inequality
f (S  , H) ≤2
 S  [t − τ ]H  [τ ] 
Hk [τ ] kk
Y[t]−[t]XkXk [t]Hk [τ ] kk,t,τ+ 2λ|Sk [t]|p . (14)k,t

In the same way, differentiating the right-hand side of thisinequality partially w.r.t. Hk [t] and setting it at 0, we obtainthe update formula for Hk [t]tHk [τ ] = Hk [τ ] tSk [t − τ ]Yk [t]Sk [t − τ ]Xk [t].

(15)Note, however,that as we did not take into account the con
straint t Hk [t] = 1 in the above, a convenient update procedure for Hk [t] shall consist of computing Eq. (15) andthen projecting
it onto the constraint space, that is, Hk [t] ←Hk [t]/ t Hk [t ].
2.4. Modulation frequency domain processing
As can be seen from an inspection of Eqs. (12), (15) there
is a need to compute the convolution of S  and H  , and thecross-correlations of H  and Y , H  and X  , S  and Y , and S and X  . Fortunately, they can be computed very efﬁciently
using the Fast Fourier Transform (FFT) as described below.
Let F be the Fourier transform operator such that Sk [u] =F {Sk [t]}u , Hk [u] = F {Hk [t]}u and Yk [u] = F {Yk [t]}u ,where u corresponds to the modulation frequency index andHk [u] is therefore nothing else than the modulation transfer
function (MTF). Then, the update formulae for H and S canbe rewritten as
Hk [t] = Hk [t]= Hk [t]= Hk [t]F −1 {F {Sk [t]}∗u F {Yk [t]}u }t
F −1 {F {Sk [t]}∗u F {Xk [t]}u }t

F −1 {F {Sk [t]}∗u F {Yk [t]}u }t
F −1 {F {Sk [t]}∗u F {Sk [t]}u F {Hk [t]}u }t
F −1 {Sk ∗ [u]Yk [u]}t,F −1 {|Sk [u]|2 Hk [u]}t(16)F −1 {Hk ∗ [u]Yk [u]}t
Sk [t] = Sk [t] −1 , (17)F {Sk [u]|Hk [u]|2 }t + λp|Sk [t]|p−1
from which we see that the computation of the multiplicative factors can be performed in the modulation frequency
domain, which is fast and also easy to implement.
2.5. Interpretation as ‘Diagonal’ Sparse NMFD
The algorithm presented here turns out to be a special caseof the non-negative matrix factor deconvolution (NMFD)[19]
when λ = 0. While NMF uses a model of the form Y W U , NMFD employs an extended modelY J

j→Wj U ,(18)j=1where Y is the matrix to be decomposed, and W t and H are
the bases and weight matrices. The j → operator shifts thecolumns of its argument by j −1 positions to the right, i.e.,
1→1 2 3 4
1 2 3 4
A=, A =
(19)
5 6 7 85 6 7 8

2→3→0 1 2 3
0 0 1 2A =, A =.
(20)
0 5 6 70 0 5 647

We will write Eq. (4) using the notation deﬁned above. LetH t be the diagonal matrix whose entries of the main diagonal
are H1 [t], H2 [t], · · · , HK [t] such that⎞⎛H1 [t]
0..
⎠,Ht ≡ ⎝(21).
0HK [t]and S be the ‘spectrogram’ matrix of clean speech such that⎞⎛S1 [1] · · · S1 [T ].. ⎠ .
S ≡ ⎝ ...
(22).
SK [1] · · · SK [T ]Eq. (4) can then be written in a matrix notation
X=T


t→

Ht S ,(23)t=1which has the same form as the NMFD model. This method
can thus be explained as a particular case of NMFD where thebases matrices H t are constrained to diagonal matrices.
By analogy with NMFD, the temporal envelope of thespectrogram blurring effect caused by reverberation corresponds to the temporal evolution of the basis component andinterestingly, the clean speech spectrogram corresponds tothe activation matrix. The reverberant spectrogram is thus
considered to comprise blurring envelopes ‘activated’ by thecomponents of a clean speech spectrogram.
3. EXPERIMENTS
In this section we report some results for speech data recorded
under several different reverberant conditions. All the speech
data were monaural and sampled at 16kHz. The STFT wascomputed using a Hanning window that was 64ms long with
2−pa 32ms overlap.
 p and λ were set at p = 1.2 and λ = E
where E = k,t Yk [t] × 10−8 . The algorithm was run for
20 iterations. Sk [t] was initially set equal to Yk [t]. Hk [t]were initially set at a decaying exponential envelope. Theﬁnal reconstruction of the speech waveform was performed
using the original phase function of the observation.
For the ﬁrst experiment we tested our method on a synthesized reverberant speech signal. The test data were created byconvolving a clean speech signal from a female speaker, excerpted from the ATR speech database, with a room impulse
response, measured in a room with an RT60 of 0.5s. Fig.
1 (a), (b) show the spectrograms of the clean speech signal
and the test data, respectively. The spectrograms of the dereverberated signal obtained with the conventional and present
methods can be seen in Fig. 1 (c), (d). Here and subsequently,
when we refer to the conventional method, we are referring tothe single channel version of [8]. We used signal-to-noise
ratio (SNR) as a quantitative measure of the dereverberation
performance. As a result, the present method improved theSNR from 2.30 to 2.94dB, while the conventional method wasonly able to improve it to 2.53dB.
For the second exeperiment we used a speech signal
recorded in a reverberant room. The clean speech signal was

250

Frequency bin #Frequency bin #250
200150
100
50

200
150
100
50

10203040

Frame #

50607020(a) Clean speech (female speaker)
Frequency bin #Frequency bin #2001501005080

Frame #

100120200
150
10050

10203040

Frame #

50607080

20(b) Synthesized reverberant speech (RT60 =0.5s)

406080

Frame #

100120(b) Speech data recorded by a distant microphone in a reverberant room

250

250

Frequency bin #Frequency bin #

60

250

250

200
150
100
50200
150
100
50

10203040

Frame #

5060708020406080

Frame #

100120(c) Dereverberated signal obtained with the 1ch version of [8](c) Dereverberated signal obtained with the 1ch version of [8]

250

250

Frequency bin #Frequency bin #40(a) Clean speech data (female speaker)

200150100502001501005010203040

Frame #

5060708020406080

Frame #100120(d) Dereverberated signal obtained with the present method(d) Dereverberated signal obtained with the present method

Fig. 1. Test on synthesized reverberant speech data

Fig. 2. Test on synthesized reverberant speech data

played by a loudspeaker and recorded by a distant microphone. The loudspeaker was carried by a participant, whowas asked to move around the room. Fig. 2 (a), (b) show
the spectrograms of the clean speech signal and the test data,respectively. The spectrograms of the dereverberated signal
obtained with the conventional and present methods can beseen in Fig. 1 (c), (d). The present method improved the SNRfrom 1.29 to 1.68dB, while the conventional method wasonly able to improve it to 1.52dB. The workstation used toperform the experiments had a Core 2 Duo processor with a
2.66GHz clock speed and a 1.99GB memory. The algorithm
was implemented in Matlab on a Windows platform. Thealgorithm usually converged within fewer than 20 iteration
cycles at a real time factor of around 1/3.

non-negativity constraints, and an iterative algorithm for its
optimization, which can be seen as a special case of NMFD.
The algorithm presented in this paper is fast, easy to implement, and robust to speaker movement.

4. CONCLUDING REMARKSIn this paper we developed a new dereverberation method designed to recover the subband envelope of an original speech
signal from its reverberant version. The discussions in thispaper mainly focused on the validity of the subband envelope model of reverberant speech, problem setting based on
the non-negativity and the sparse nature of speech spectrograms leading to a regularized optimization problem with48

5. REFERENCES[1] B. Yegnanarayana, Proc. ICASSP’98, 1, 405–408, 1998.
[2] M.S. Brandstein, Proc. ICASSP’98, 6, 3613–3616, 1998.
[3] B. W. Gillespie et al., ICASSP’01, 6, 3701–3704, 2001.
[4] T. Nakatani et al., IEEE Trans. ASLP, 15(1), 80–95, 2007.
[5] K. Kinoshita et al., Proc. ICASSP’06, 1, 817–820, 2006.
[6] T. Yoshioka et al., Proc. IWAENC’06, 2006.
[7] T. Nakatani et al., Proc. ICASSP’07, 193–197, 2007.
[8] T. Nakatani et al., Proc. ICASSP’08, 85–88, 2008.
[9] T. Yoshioka et al., Proc. ICASSP’08, 4585–4588, 2008.
[10] T. Langhans et al., ICASSP’82, 156–159, 1982.
[11] J. Mourjopoulos et al., ICASSP’83, 1144–1147, 1983.
[12] H. G. Hirsch, EURASIP, 1177–1180, 1988.
[13] H. Hermansky et al., ICASSP’95, 405–408, 1995.
[14] C. Avendano et al., Proc. ICSLP’96, 2, 889–892, 1996.
[15] S. Hirobayashi et al., IEICE Trans. A, J81A(10), 1323–1330,
1998.

[16] M. Unoki et al., ICASSP’03, 1, 840–843, 2003.
[17] E. A. P. Habets, Proc. ProRISC’04, 250–254, 2004.
[18] D. D. Lee et al., Proc. NIPS’00, 556–562, 2000.
[19] P. Smaragdis et al., Proc. ICA’04, 494–499, 2004.