MAXIMUM LIKELIHOOD APPROACH TO SPEECH ENHANCEMENT FORNOISY REVERBERANT SIGNALS
Takuya Yoshioka, Tomohiro Nakatani, Takafumi Hikichi, and Masato Miyoshi
NTT Communication Science Laboratories, NTT Corporation
2-4, Hikari-dai, Seika-cho, Soraku-gun, Kyoto 619-0237, Japan
Noise
Dt ,ω

ABSTRACTThis paper proposes a speech enhancement method for signals contaminated by room reverberation and additive background noise.
The following conditions are assumed: (1) The spectral componentsof speech and noise are statistically independent Gaussian randomvariables. (2) The convolutive distortion channel is modeled as anauto-regressive system in each frequency bin. (3) The power spectral density of speech is modeled as an all-pole spectrum, while thatof noise is assumed to be stationary and given in advance. Under
these conditions, the proposed method estimates the parameters of
the channel and those of the all-pole speech model based on the maximum likelihood estimation method. Experimental results showedthat the proposed method successfully suppressed the reverberation
and additive noise from three-second noisy reverberant signals when
the reverberation time was 0.5 seconds and the reverberant signal tonoise ratio was 10 dB.
Index Terms— Speech enhancement, Reverberation, Additivenoise, Maximum likelihood estimation
1. INTRODUCTION
Speech signals captured by microphones in rooms are often corrupted by reverberation and background noise. The distortion system
of interest is modeled as Fig. 1. Recovering the clean speech signals
from the distorted signals observed at microphones is indispensable
for many audio applications working in such realistic environments.
Thus far, it has been a difﬁcult problem to recover clean speech
signals from their noisy and reverberated version. Conventional
speech dereverberation techniques based on a subspace method arequite sensitive to additive noise [1]. In recent years, we have developed dereverberation methods that are based on statistical models ofspeech [2, 3]. Those methods calculate maximum likelihood (ML)estimates of convolutive distortion channels, where the likelihood
functions are deﬁned based on the speech models. They exhibit robustness against noise with a relatively low noise level. However,recovering the clean speech signals in highly noisy reverberant environments remains a challenging problem [4].
In this paper, we propose a speech enhancement method for signals corrupted by both convolutive distortion channels and additive
stationary noise. The proposed method is based on the ML estimation methodology. An observed signal is ﬁrst analyzed with a shorttime Fourier transform (STFT), and the resultant complex spectrogram is then processed. We assume that the spectral components ofspeech and noise signals are statistically independent Gaussian random variables. The convolutive distortion channel is modeled as anauto-regressive (AR) system in each frequency bin. The power spectral density (PSD) of speech is modeled by an all-pole spectrum. The

1-4244-1484-9/08/$25.00 ©2008 IEEE

4585

Clean
speechSt ,ω

Convolutive
distortion
channel Reverberant
speechX t ,ω

+Noisy reverberant
speech observed ata microphone
Yt ,ω

Fig. 1. Distortion system

PSD of noise is assumed to be known in advance. Under these conditions, we can deﬁne the likelihood function of the parameters of thechannel and those of the all-pole speech model. By maximizing thislikelihood function, the proposed method estimates the parameters.
The proposed method was tested on signals corrupted by a channel
with a reverberation time of 0.5 seconds and by additive noise witha reverberant signal to noise ratio (RSNR) of 10 dB. The proposed
method successfully suppressed the reverberation and noise from theobserved three-second signals.
2. TASK AND APPROACH
2.1. Task formulation
Let Xt,ω and Dt,ω , respectively, denote the complex spectral components of a reverberant speech signal and a noise signal at the t-thtime frame and a frequency bin with frequency ω. The spectral component Yt,ω of an observed signal is given by
Yt,ω = Xt,ω + Dt,ω .

(1)

Let St,ω denote the complex spectral component of a clean speech
signal. Provided that a convolutive distortion channel is modeled byan AR system in each frequency bin, the reverberant speech spectral
component Xt,ω is then represented as
K

Xt,ω =∗gk,ω
Xt−k,ω + St,ω ,

(2)

k=1

where gk,ω is the k-th regression coefﬁcient of the distortion channel
at frequency ω, K is the regression order of the distortion channel,
and superscript ∗ denotes a complex conjugate operation. Note thefollowing three points as regards the convolutive distortion model(2).
(p1) It is assumed that the convolutive distortion channel has nocross-band components [5].

ICASSP 2008

spectrogram Y and the parameter set Θ. Therefore, if the true values of Θ were known, clean speech spectrogram S could then be
estimated by using this distribution based on (2). At this point, ourinterest is now directed at estimating the unknown parameter set Θ.
To estimate Θ, we employ the maximum likelihood (ML) estimation method. With this method, once the spectrogram Y of thenoisy reverberant speech is observed, Θ̂ that maximizes the log likelihood function log p(Y|Θ) is calculated. In the next section, we
speciﬁcally derive the ML estimator of Θ.

(p2) The convolutive distortion channel is assumed to be invertiblein each frequency bin. By working in the frequency domain,
the noninvertibility problem of a non-minimum phase channel in the full band may be mitigated [6]. In fact, the proposed
method was tested by using a non-minimum phase channel1 .
(p3) (2) ignores the part of the convolutive distortion that is shorter
than the frame size.
The task to be solved is to estimate the clean speech spectrogram
S = {St,ω }t,ω from the observed spectrogram Y = {Yt,ω }t,ω .
{St,ω }t,ω means the set of St,ω over all the frame indices and discrete frequencies. Hereafter, we use the same notation. We also represent the spectrogram of the reverberant speech as X = {Xt,ω }t,ω .
To make the estimation task tractable, we introduce statisticalmodels of the speech and noise spectra as constraints on the speechand noise. Under the abovementioned Gaussian assumptions regarding the speech and noise spectral components, the probability density functions (PDFs) of these components are given byp(St,ω | s λt (ω)) =p(Dt,ω | d λt (ω)) =1|St,ω |2exp −π s λt (ω)
s λt (ω)

3.1. Expectation maximization algorithm
To maximize the log likelihood function log p(Y|Θ), we employ the
expectation maximization (EM) algorithm because the model (7) involves the latent variable S. The EM algorithm iteratively maximizes the expected log likelihood function of complete data {Y, X }(often called an auxiliary function) as follows.

(3)





1
|Dt,ω |2exp −,π d λt (ω)
d λt (ω)

(4)





where s λt (ω) and d λt (ω) are the speech and noise PSDs at thet-th frame, respectively. St,ω and Dt ,ω are considered to be statistically independent for any t, ω, t , and ω  . Furthermore, if(t, ω) = (t , ω  ), St,ω and St ,ω (and Dt,ω and Dt ,ω ) are assumed to be statistically independent, although they are dependent
in an exact sense owing to the window function and frame overlap. The speech PSD is, for simplicity, assumed to be of the all-pole
form, though more complicated models may be incorporated. Thus,
we haves λt (ω) =3. DERIVATION OF MAXIMUM LIKELIHOODESTIMATOR

2
s σt
jω
|At (e )|2

E-step Calculate p(X |Y, Θ(i) ), i.e. the conditional posterior of thereverberant speech spectrogram, where Θ(i) is the estimateof Θ in the i-th iteration.
M-step Update the estimate of Θ by Θ(i+1) = argmaxΘ Q(Θ|Θ(i) ),(5)

At (z) =1 − at,1 z −1 − · · · − at,P z −P ,(6)

where P is the number of poles and {at,k , s σt2 }k are the all-polemodel parameters at the t-th frame. The noise PSD is assumed tobe known in advance, and might be obtained from intervals without
speech.
By putting (1) to (6) together, we have the PDF p(Y, X |Θ) of
the observed spectrogram Y and the reverberant speech spectrogram
X given parameter set Θ = {gk,ω , at,l , s σt2 }k,ω,l,t as|At (ejω )|2|Yt,ω − Xt,ω |2p(Y, X |Θ) ∝2 exp −λ (ω) s σt
d λt (ω)
t,ω d t
t,ω






where
Q(Θ|Θ(i) ) = log p(Y, X |Θ)p(X |Y,Θ(i) ) .

Below we embody the E-step and M-step, and then summarize anddiscuss the algorithm of the proposed method.
3.2. Conditional posterior of reverberant speech spectrogram
In the E-step, we calculate the conditional posterior p(X |Y, Θ) of
the reverberant speech spectrogram. Here, we only show the derived p(X |Y, Θ) owing to limited space. Let us deﬁne ω and ω
∗
∗
∗
∗ H
as ω = [XT,ω
, · · · , X1,ω]H and ω = [YT,ω
, · · · , Y1,ω] , respectively, where T is the number of time frames and superscript H
denotes a conjugate transpose operation. Then, we havep(X |Y, Θ) =

N {




ω;

ω , Σω }(9)ωH
H H −1(Bω BωH )ω =(Bω Bω + Gω Aω Aω Gω )H −1
Σω =(Bω BωH + Gω Aω AH,
ω Gω )



(8)

ω

(10)
(11)



+

|At (ejω )|2 |Xt,ω −


2
s σt

K
∗
2
k=1 gk,ω Xt−k,ω |

.

(7)





We can also obtain the PDF p(Y|Θ) of the observed spectrogram Y
given Θ by integrating out X from p(Y, X |Θ).
2.2. Maximum likelihood estimation approach
As shown later, we can derive the posterior distribution p(X |Y, Θ)of the reverberant speech spectrogram X conditioned on observed
1 For simplicity, we limit the derivation to a single-microphone case.
However, the generalization of the proposed method to a multi-microphonecase, where the channel is always invertible [7], is strainghtforward.

4586

where N { ; , Σ} denotes the complex Gaussian PDF ofwith meanand covariance matrix
random variableΣ.
In (11), Gω is a T -dimensional Toeplitz matrix hav∗∗
, · · · , −gK,ω
, 0, · · · , 0]H as its ﬁrst column anding [1, −g1,ω
[1, 0, · · · , 0] as its ﬁrst row. Aω and Bω are T -dimensional
diagonal matrices whose t-th diagonals are 1/ s λT −t+1 (ω) and1/ d λT −t+1 (ω), respectively.
(11) indicates that we can calculate ω independently for eachfrequency and that the mean, or MMSE estimate, of ω is the output
of a Wiener ﬁlter applied to ω . In particular, if Gω is an identity
matrix (i.e. a reverberation-free condition), (11) reduces to Xt,ω ={ s λt (ω)/( s λt (ω) + d λt (ω))}Yt,ω , which corresponds to the wellknown Wiener-ﬁlter-based noise suppression [8].


3.3. Maximization of auxiliary function(i)− log( s σt2 )Q(Θ|Θ ) =
t=1

−

3000

0.4
0.2
0
−0.2

ω

2



2
s σt

K
∗
2
k=1 gk,ω Xt−k,ω |p(X |Y,Θ(i) ). (12)

1500

0.5

1

1.5

2

2.5

0

3

0

0.5

1

1.5

2

2.5

2

2.5

2

2.5

2

2.5

Time (sec)Time (sec)∗∗
To derive a condition that ω = [g1,ω
, · · · , gK,ω]H maximizing(12) must satisfy, let us take the gradient of (12) with respect to ωand set it at zero. See [9] for the gradient of a function of complex
variables. Organizing the resultant equation leads to

2000

500

−0.6
−0.8
0

2500

1000

−0.4

|At (e )| |Xt,ω −
jω

3500

Frequency (Hz)

T

4000

0.6

Amplitude

The M-step maximizes the auxiliary function Q(Θ|Θ(i) ) deﬁned by(8). The auxiliary function is represented as

1
0.8

1

4000

0.8

3500



0.6

3000



ω = x ω.
Amplitude

0.2
0
−0.2

(13)



Frequency (Hz)

x Rω

0.4

2000
1500
1000

−0.4

In (13), x Rω and x ω are deﬁned as

2500

500

−0.6



−0.8
0

T

|At (ejω )|2
x Rω =
2
s σt
t=1

1

1.5

2

2.5

0

3

0

0.5

1

t−1,ω


H
t−1,ω p(X |Y,Θ(i) )


(14)

1

4000

0.8

3500

0.6

∗
∗
, · · · , Xt−K+1,ω
]H .
where t,ω = [Xt,ωOn the other hand, by taking the gradient of (12) with respect toHand setting it at zero, we obtain the condit = [at,1 , · · · , at,P ]tion for t as

Frequency (Hz)(15)

Amplitude





3000

0.4

∗
t−1,ω Xt,ω p(X |Y,Θ(i) ) ,1.5Time (sec)Time (sec)T|At (ejω )|2

x ω =
2
s σt
t=1

0.5

0.2
0
−0.2

2500
2000
1500
1000



−0.4

500

−0.6



−0.8
0

0.5

1

t = s t.

(16)In (16), s Rt and s t are deﬁned as(P × P matrix)

(17)
(18)



|St,ω |2 p(X |Y,Θ(i) ) ejωk ,

(19)

4000

0.6

3500

0.4

3000

0

0.5

1

2500
2000
1500

500

−0.6

0.5

1

1.5

2

2.5

0

3

0

0.5

1

|At (ejω )|2 |St,ω |2 p(X |Y,Θ(i) ) .

(20)

1.5

Time (sec)Time (sec)where St,ω is given by (2) and N is the frame size. Also, the condition for s σt2 can be easily derived as

1.5

1000

−0.4

−0.8
0

ω

0

Time (sec)

0.8

−0.2

H
s t =[ s rt (1), · · · , s rt (P )]2s σt =

Amplitudes Rt =[ s rt (l − m)]1≤l,m≤P

0

3

0.2



1
s rt (k) =N

2.5

Frequency (Hz)

2Time (sec)



s Rt

1.5

Fig. 2. Waveforms (left panels) and spectrograms (right panels) of,
from top to bottom, clean speech, reverberant speech, noisy reverberant speech, and recovered speech.

ωNote that St,ω , and hence, t and s σt2 are dependent on ω .
We ﬁnd that the solutions for (13) and for the simultaneous equation of (16) and (20) are mutually dependent. Hence, there is no
closed-form solution for the simultaneous equation of (13), (16) and(20). We therefore propose updating the estimates of { ω }ω and{ t , s σt2 }t by using the following rules:

3.4. Summary and discussion
To summarize Sections 3.1 to 3.3, the algorithm for estimating theparameter set Θ is drawn as follows.







(i+1)
(i) −1(i)= s Rt
s t
t

(21)



2 (i+1)
=
s σt

(i)

|At (ejω )(i+1) |2 |St,ω |2 p(X |Y,Θ(i) )

(22)



(23)



(i)

(i)

(i)

(i)
, and x ω(i)In (21) to (23), s Rt , s t , At (ejω )(i) , St,ω , x Rωare calculated based on (17), (18), (6), (2), (14), and (15), respectively, where the values of Θ are set at Θ(i) . This variant of the EM
algorithm, where the M-step is replaced by several conditional maximization steps, is called the expectation conditional maximization(ECM) algorithm [10]. The ECM algorithm is proven to converge toa local optimum.
(s2) Calculate the meansp(X |Y, Θ(i) ) of (11).
(i+1)(s3) Calculate { tto (23).
ω(i+1)(i+1) −1
(i+1)
= x Rω.
x ωω(s1) Set initial parameter estimate Θ(0) and iteration index i at 0.

, s σt2

and

(i+1)

covariances(i+1)}t and { ωthatspecify}ω by using (21)(s4) Increment i and return to (s2) unless convergence is reached.



4587

It is noteworthy that the proposed method is an extension of thenoise suppression method based on the all-pole speech model [11] aswell as that of the dereverberation method proposed in [2]. In fact, if
we ignore the convolutive distortion channel by setting its order K
at 0, the proposed method becomes identical to the noise suppression method described in [11]. On the other hand, if we ignore the5. CONCLUSION

Table 1. SASNR improvement in dB.
Condition ( d: enabled, ×: disabled)ddNoise suppression
×dd
Dereverberation
×
Male
4.25
1.80
7.77
4.67
1.17
7.67
Femaleadditive noise by setting the noise variance d λt (ω) at 0 for any ω,
the proposed method is reduced to the frequency domain version of
the dereverberation method described in [2].
4. EXPERIMENTAL RESULTS
We conducted experiments to evaluate the performance of the proposed method. Japanese utterances of 10 males and 10 females weretaken from the JNAS database. All the utterances were three seconds
long. The signals were sampled at 8 kHz and quantized at 16-bit resolution. The speech signals of these utterances were convolved witha non-minimum phase room impulse response measured in a roomwith a reverberation time of 0.5 seconds to synthesize reverberant
signals. The reverberant signals were then contaminated by additive white Gaussian noises with a reverberant signal to noise ratio(RSNR) of 10 dB. The system paramters were set as follows: the
STFT frame size was 256 samples, the STFT frame shift was 128
samples, the window function was a Hanning window, the order ofthe convolutive distortion channel K was 30 for all frequency bins,
the assumed number of poles of the speech signals P was 12, andthe number of iteratons for the EM algorithm was 5. The proposed
method was implemented as a Matlab program. Each iteration was
performed with a real time factor of 0.67 on a Linux PC equippedwith a 3.6 GHz Pentium 4 processor.
The speech enhancement performance was evaluated by using asegmental amplitude signal to noise ratio (SASNR) deﬁned as
T

SASNR =

1
10 log10
T t=1




2
ω |St,ω |2
ω ||Ŝt,ω | − |St,ω ||

(dB).

(24)

We ﬁrst provide an example of the experimental results. Fig. 2
shows the waveforms and spectrograms of clean speech, reverberant speech, noisy reverberant speech, and recovered speech. It canbe seen that the proposed method successfully suppressed both thereverberant components and the noise, although the speech spectral components were underestimated in the higher frequency band,
where the signal to noise ratios were very low.
Table 1 summarizes the improvement in SASNR for each gender. The proposed method improved the SASNRs by an averageof 7.72 dB. When the dereverberation function is disabled by setting
the assumed convolutive distortion channel order K at 0, the average
improvement in the SASNRs was 4.46 dB. On the other hand, when
the noise suppression function is disabled by setting d λt (ω) at 0, theaverage improvement of the SASNRs fell down to 1.49 dB. Therefore, it can be concluded that the proposed method could enhance
the quality of the noisy reverberant speech signals by integrating the
dereverberation and noise suppression.
The residual distortion is attributed partly to the part of the convolutive distortion shorter than the frame size, which the proposed
method ignores. The inaccuracy of the convolutive distortion model(2) may be another reason for the residual distortion.

4588

In this paper, we described an ML estimation method for the regression coefﬁcients of a convolutive distortion channel in each frequency bin and the parameters of an all-pole speech model. An EM
algorithm was employed, in which the estimates of all the parameters
and those of the reverberant speech spectrograms are updated alternately. The derived algorithm is an extension of the noise suppression method described in [8] and the speech dereverberation method
reported in [2], both of which are based on the all-pole speech model.
Experimental results showed that the proposed method could suppress the reverberation and additive noise from three-second noisy
reverberant signals.
Future work includes adaptive estimation of the parameters of
the channel and the speech model in order to cope with a timevarying convolutive distortion channel. Adaptive estimation of thenoise PSD should also be addressed.
6. REFERENCES[1] T. Hikichi, M. Delcroix, and M. Miyoshi, “Blind dereverberation based on estimates of signal transmission channels without precise information on channel order,” in Proc. Int’l Conf.
Acoust., Speech, Signal Process., 2005, pp. 1069–1072.
[2] T. Yoshioka, T. Hikichi, and M. Miyoshi, “Dereverberationby using time-variant nature of speech production system,”
EURASIP J. Adv. Signal Process., vol. 2007, 2007, Article
ID 65698, 15 pages, doi:10.1155/2007/65698.
[3] T. Nakatani, B. H. Juang, T. Yoshioka, K. Kinoshita, and
M. Miyoshi, “Importance of energy and spectral features inGaussian source model for speech dereverberation,” in Proc.
IEEE Worksh. Appl. Signal Process. Audio, Acoust., 2007, pp.
299–302.
[4] K. Kinoshita, T. Nakatani, M. Delcroix, and M. Miyoshi,“Multi-step linear prediction based speech dereverberation innoisy reverberant environment,” in Proc. Interspeech, 2007,
pp. 854–857.
[5] Y. Avargel and I. Cohen, “System identiﬁcation in the shorttime Fourier transform domain with crossband ﬁltering,” IEEE
Trans. Audio, Speech, Lang. Process., vol. 15, no. 4, pp. 1305–
1319, 2007.
[6] H. Wang and F. Itakura, “Dereverberation of speech signals
based on sub-band envelope estimation,” IEICE Trans. Fund.,vol. E74-A, no. 11, pp. 3576–3583, 1991.
[7] M. Miyoshi and Y. Kaneda, “Inverse ﬁltering of room acoustics,” IEEE Trans. Acoust., Speech, Signal Process., vol. 36,
no. 2, pp. 145–152, 1988.
[8] J. S. Lim and A. V. Oppenheim, “Enhancement and bandwidthcompression of noisy speech,” Proc. IEEE, vol. 67, no. 12, pp.
1586–1604, 1979.
[9] S. Haykin, Adaptive Filter Theory, Prentice Hall, fourth edition, 2001.
[10] X. L. Meng and D. B. Rubin, “Maximum likelihood estimation
via the ECM algorithm: a general framework,” Biometrika,
vol. 80, no. 2, pp. 267–278, 1993.
[11] J. S. Lim and A. V. Oppenheim, “All-pole modeling of degraded speech,” IEEE Trans. Acoust. Speech, Signal Process.,
vol. 26, no. 3, pp. 197–210, 1978.