Modulation-Domain Kalman Filtering for Monaural Blind Speech Denoising and Dereverberation
Nikolaos Dionelis Mike Brookes





Modulation-Domain Kalman Filtering for Monaural Blind Speech Denoising
and Dereverberation 
    N. Dionelis, <https://www.commsp.ee.ic.ac.uk/ sap/people-nikolaos-dionelis/>,
M. Brookes, Member, IEEE
    December 30, 2023
==========================================================================================================


We describe a monaural speech enhancement algorithm based on modulation-domain
Kalman filtering to blindly track the time-frequency log-magnitude
spectra of speech and reverberation. We propose an adaptive algorithm
that performs blind joint denoising and dereverberation, while accounting
for the inter-frame speech dynamics, by estimating the posterior distribution
of the speech log-magnitude spectrum given the log-magnitude spectrum
of the noisy reverberant speech. The Kalman filter update step models
the non-linear relations between the speech, noise and reverberation
log-spectra. The Kalman filtering algorithm uses a signal model that
takes into account the reverberation parameters of the reverberation
time, T_60, and the direct-to-reverberant energy ratio (DRR)
and also estimates and tracks the T_60 and the DRR in every frequency
bin in order to improve the estimation of the speech log-magnitude
spectrum. The Kalman filtering algorithm is tested and graphs that
depict the estimated reverberation features over time are examined.
The proposed algorithm is evaluated in terms of speech quality, speech
intelligibility and dereverberation performance for a range of reverberation
parameters and SNRs, in different noise types, and is also compared
to competing denoising and dereverberation techniques. Experimental
results using noisy reverberant speech demonstrate the effectiveness
of the enhancement algorithm.













§ INTRODUCTION 


[lines=2]Nowadays, technology is ever evolving with
tremendous haste and the demand for speech enhancement systems is
evident. Speech enhancement in noisy reverberant environments, for
human listeners, is challenging. Speech is degraded by noise and reverberation
when captured using a near-field or far-field distant microphone <cit.>
<cit.>. A room impulse response (RIR) can include components
at long delays, hence resulting in reverberation and echoes <cit.>
<cit.>. Reverberation is a convolutive distortion that can be
quite long with a reverberation time, T_60, of more than 0.8
s. Due to convolution, reverberation induces long-term correlation
between consecutive observations. Reverberation and noise, which can
be stationary or non-stationary, have a detrimental impact on speech
quality and intelligibility. Reverberation, especially in the presence
of non-stationary noise, damages the intelligibility of speech.firstpage

The direct to reverberant energy ratio (DRR) and the reverberation
time, T_60, are the two main parameters of a reverberation model
<cit.> <cit.>. The DRR describes reverberation in the
space domain, depending on the positions of the sound source and the
receiver. The T_60 is the time interval required for a sound
level to decay 60 dB after ceasing its original stimulus. The reverberation
time, when measured in the diffuse sound field, is independent of
the source to microphone configuration and mainly depends on the room.
The impact of reverberation on auditory perception depends on the
T_60. If the T_60 is short, the environment reinforces the
sound which may enhance the sound perception <cit.> <cit.>.
On the contrary, if the T_60 is long, spoken syllables interfere
with future spoken syllables. Reverberation spreads energy over time
and this smearing across time has two effects: (a) the energy of individual
phonemes spreads out in time and, hence, plosives have a delayed decay
and fricatives are smoothed, and (b) preceding phonemes blur into
the current phonemes.

The aim of speech enhancement is to reduce and ideally eliminate the
effects of both noise and reverberation without distorting the speech
signal <cit.>. Enhancement algorithms typically aim to suppress
noise and late reverberation because early reverberation is not perceived
as separate sound sources and usually improves the quality and intelligibility
of speech. Noise is assumed to be uncorrelated with speech, early
reverberation is correlated with speech and late reverberation is
commonly assumed to be uncorrelated with speech <cit.>
<cit.>.

Speech enhancement can be performed in different domains. The ideal
domain should be chosen such that (a) good statistical models of speech
and noise exist in this domain, and (b) speech and noise are separable
in this domain. Speech and noise are additive in the time domain and
the Short Time Fourier Transform (STFT) domain <cit.> <cit.>.
The relation between speech and noise becomes progressively more complicated
in the amplitude, power and log-power spectral domains. Noise suppression
algorithms usually operate in a time-frequency STFT domain and these
techniques have been extended to address dereverberation. In <cit.>,
spectral enhancement methods based on a time-frequency gain, originally
developed for noise suppression, have been modified and employed for
dereverberation. Such algorithms suppress late reverberation assuming
that the early and late reverberation components are uncorrelated.
The spectral enhancement methods in <cit.> <cit.>
estimate the late reverberant spectral variance (LRSV) and use it
in the place of the noise spectral variance, reducing the problem
of late reverberation suppression to that of estimating the LRSV <cit.>.
In <cit.>, blind spectral weighting is employed to reduce the
overlap-masking effect of reverberation using an uncorrelated and
additive assumption for late reverberation.

Dereverberation algorithms that leave the phase unaltered and operate
in the amplitude, power or log-power spectral domains are relatively
insensitive to minor variations in the spatial placement of sources
<cit.>. Two criticisms of spectral enhancement algorithms based
on LRSV reverberationnoiseestimation are that they introduce musical
noise and suppress speech onsets when they over-estimate reverberation
<cit.>. The LRSV estimator in <cit.>, which is a continuation
of <cit.>, models the RIR in the STFT domain and not in
the time domain <cit.> <cit.>, using the same model
of the RIR that is attributed to J. Polack or J. Moorer <cit.>.
Reverberation is estimated in <cit.> considering the STFT energy
contribution of the direct path of speech and an external T_60
estimate.

Modelling the speech temporal dynamics is beneficial when the T_60
is long and the DRR is low <cit.> <cit.>. Joint denoising
and dereverberation using speech and noise tracking is performed in
<cit.>. The SPENDRED algorithm <cit.> <cit.>, which
is a model-based method with a convolution model for reverberation
based on the T_60 and the DRR, considers the speech temporal
dynamics. SPENDRED employs a parametric model of the RIR <cit.>
and performs frequency-dependent and time-varying T_60 and DRR
estimation. However, unless the source or the microphone are moving,
the T_60 and the DRR will be constant throughout the recording.
The SPENDRED algorithm assumes that DRR≥ 1 - (10^-6)^L/T_60
where L is the acoustic frame increment. For example, when T_60=0.4
s and L=5 ms, then DRR≥ -8 dB is assumed. In addition,
SPENDRED performs intra-frame correlation modelling, which can be
beneficial in adverse conditions, while typical algorithms decouple
different frequency dimensions <cit.>.

Statistical-based models, such as the SPENDRED algorithm, describe
reverberation by a convolution in the power spectral domain while
LRSV models describe reverberation as an additive distortion in the
power spectral domain <cit.> <cit.>. A model with an infinite
impulse response is used either with the two parameters of the T_60
and the DRR, as in <cit.> <cit.>, or with a finite number
of parameters. The infinite-order convolution model of reverberation
with the T_60 and the DRR is sparse and contrasts with the higher-order
autoregressive processes in the complex STFT domain, used in <cit.>
<cit.>.

The algorithms described in <cit.>, <cit.> and <cit.>
create non-linear observation models of noisy reverberant speech in
the log Mel-power spectral domain, using the reverberation-to-noise
ratio (RNR). As discussed in <cit.>, phase differences in Mel-frequency
bands have different properties from phase differences in STFT bins.
The phase factor between reverberant speech and noise is different
from that between speech and noise <cit.>. In <cit.>,
the phase factor between reverberant speech and noise in Mel-frequency
bands is examined.

In noisy reverberant conditions, finding the onset of speech phonemes
and determining which frames are unvoiced/silence is difficult, due
to the smearing across time, often leading to noise over-estimation.
The concatenation of different techniques for denoising and dereverberation
has lower performance than unified methods due to over-estimating
noise when estimating noise and reverberation separately <cit.>
<cit.>.

Despite the claim that it is inefficient to perform a two step procedure
that is comprised of denoising followed by dereverberation <cit.>
<cit.>, long-term linear prediction with pre-denoising can be
used to suppress noise and reverberation. With the weighted prediction
error (WPE) algorithm <cit.> <cit.>, reverberation is
represented as a one-dimensional convolution in each frequency bin.
In <cit.>, the WPE algorithm is discussed along with inter-frame
correlation. In <cit.>, the WPE algorithm is used in the complex
STFT domain performing batch processing and iteratively estimating,
first, the reverberation prediction coefficients and, then, the speech
spectral variance. The WPE linear filtering approach, which can be
employed in the power spectral domain <cit.> <cit.>, takes
into account past frames, from the 3-rd to the 40-th past frame
<cit.> <cit.>.

This paper presents an adaptive denoising and dereverberation Kalman
filtering framework that tracks the speech and reverberation spectral
log-magnitudes. In this paper, we extend the enhancer in <cit.>
to include dereverberation. Enhancement is performed using a Kalman
filter (KF) to model inter-frame correlations. We use an integrated
structure of two parallel signal models to track speech, reverberation
and the T_60 and DRR reverberation parameters. The T_60
and the DRR are updated in every frame to improve the estimation of
the speech log-magnitude spectrum. We create an observation model
and a series of non-linear KF update steps performing joint noise
and reverberation suppression by estimating the first two moments
of the posterior distribution of the speech log-spectrum given the
noisy reverberant log-spectrum. The log-spectral domain is chosen,
as in <cit.> <cit.>, because good speech models exist
in this domain. Modelling spectral log-amplitudes as Gaussian distributions
leads to good speech modelling in noisy reverberant environments since
super-Gaussian distributions that resemble the log-normal, such as
the Gamma <cit.> <cit.>, are used to model the speech
amplitude spectrum. Mean squared errors (MSEs) in the log-spectral
domain are a good measure to use for perceptual quality and speech
log-spectra are well modelled by Gaussian distributions, as in <cit.>
and <cit.>.

The structure of this paper is as follows. Section II describes the
signal model and Sec. III presents the enhancement algorithm and its
non-linear KF. The implementation and the validation of the algorithm
are in Sec. IV. The algorithm's evaluation is in Sec. V. Conclusions
are drawn in Sec. VI.



§ SIGNAL MODEL AND NOTATION


In the complex STFT domain, the noisy speech, Y_t(k), is given
by Y_t(k) = S_t(k) + R_t(k) + N_t(k) where S_t(k) is the direct
speech component, R_t(k) is the reverberant speech component and
N_t(k) is the noise, as for example in <cit.> <cit.>.
The time-frame index is t and the frequency bin index is k.
For clarity, we also define Z_t(k)=R_t(k)+N_t(k). We drop the time
and frequency indexes and we obtain Y = S + Z = S + R + N. We define
the log-magnitude spectrum of S as s=log(|S|) and we also define
r, n, y and z similarly.

In the signal model,british signal quantities with
capital letters, such as S_t, are complex numbers with magnitude
and phase values, |S_t| and ∠ S_t. In the complex STFT
domain, using a_t,b_t ∈ℝ, britishthe
reverberation signal model is given by

    R_t    =√(a_t) R_t-1exp(jθ_t)+√(b_t) S_t-1exp(jψ_t)
        =∑_τ=1^∞(∏_i=1^τ-1(√(a_t-i+1)exp(jθ_t-i+1)).
       ×.√(b_t-τ+1)S_t-τexp(jψ_t-τ+1)),


    a_t^T_60/L=10^-6,        b_t=1-a_t/DRR

where L is the acoustic frame increment. In (<ref>), the
factors exp(j θ_t) and exp(j ψ_t), where θ_t
and ψ_t are uniformly distributed phases, are used. In (2),
the DRR is defined in the power spectral domain <cit.> <cit.>
and the T_60 and the DRR are both time and frequency dependent,
as described in <cit.>.

The expression in (<ref>) is the convolution model for reverberation;
the most common reverberation model is this single-pole filter that
is described by the pole and zero positions that depend on the T_60
and the DRR <cit.> <cit.>. A convolution of infinite order
is used, with the two parameters of the T_60 and the DRR, to
describe reverberation <cit.> <cit.>. Models that describe
reverberation by a convolution are also discussed in <cit.>
<cit.>. britishThe signal model is defined
by (1) and by

    Y_t=S_t+Z_t,      Z_t=R_t+N_t,
       γ_t=0.5log(a_t),  δ_t=γ_t+r_t-1=log(√(a_t)|R_t-1|),
       β_t=0.5log(b_t),  ϵ_t=β_t+s_t-1=log(√(b_t)|S_t-1|)

britishwhere b_t > 0, 0 < a_t < 1 and γ_t < 0.

Figure 1 shows graphs of β against T_60 for a fixed DRR
and of β against DRR for a fixed T_60. If DRR = 0
dB, then b_t = 1 - a_t. If β_t = 0, then b_t=1 and DRR = 1-a_t.

Figure 2 illustrates the flowchart of the signal model. The britishreverberation
signal model in (<ref>) uses √(a_t) and √(b_t)
because the a_t and b_t reverberation parameters, in (<ref>),
and the DRR are defined in the power spectral domain, as in <cit.>
<cit.>. The a_t and b_t parameters are mapped to γ_t
and β_t using (<ref>) and (<ref>). The signals
z_t, δ_t and ϵ_t are the total distubance, the
old (decaying) reverberation and the new reverberation, respectively.
We note that z_t is defined in the first paragraph of this section
and that δ_t and ϵ_t are defined in (<ref>)
and (<ref>).

The signal model of how the reverberation parameters of γ_t
and β_t change over time is a random walk model. This is used
in the algorithm's KF prediction step for γ_t and β_t.

The signal model in Fig. 2 is directly linked to the alternating and
interacting KFs of the enhancement algorithm. The algorithm is
a collection of two KFs, the speech KF and the reverberation KF, that
estimate the speech and reverberation log-amplitude spectra and the
γ_t and β_t reverberation parameters. This KF algorithm
is described in detail in Sec. III.



§ THE SPEECH ENHANCEMENT ALGORITHM


The KF algorithm operates in the log-magnitude spectral domain, tracking
speech and reverberation. Figure <ref>
depicts the denoising and dereverberation algorithm that formulates
a model of reverberation as a first-order autoregressive process and
propagates the means and variances of the random variables. Almost
all the signals follow a Gaussian distribution and the distribution
of s_t conditioned on observations up to time τ is given
by p_s_t| τ(s) ≜ p(s_t | y_0, …, y_τ) = N(s_t| τ, Σ^(s)_t| τ).
In Fig. <ref>,
a Gaussian distribution is denoted by its mean, s_t|τ.

The core of the algorithm in Fig. <ref>
is the KF that is defined by the gray blocks in the flowchart diagram.
The non-linear KF estimates and tracks the posterior distributions
of the speech log-magnitude spectrum, s_t, the reverberation log-magnitude
spectrum, r_t, and the reverberation parameters, γ_t and
β_t.

The input to the algorithm in Fig. <ref>
is the noisy reverberant speech in the time domain. The algorithm's
first step is to perform a STFT and obtain the signal in the complex
STFT domain. The algorithm does not alter the noisy reverberant phase,
∠ Y, and uses the noisy reverberant amplitude spectrum, |Y|,
in three ways: in the speech KF prediction step, in the KF update
step and in the noise power modelling. The main part of the algorithm
is the KF and the speech KF state, s_t∈ℝ^p,
is the speech log-spectrum from the previous p frames,

    s_t=(s_t s_t-1 … s_t-p+1)^T.

The speech KF prediction step is based on autoregressive (AR) modelling
on the log-spectrum of pre-cleaned speech <cit.>. The reverberation
KF state is r_t and the KF states of the reverberation parameters
are γ_t and β_t. The KF observation is the noisy reverberant
speech log-spectrum, y_t, which is used in the KF update step to
compute the first two moments of the posterior of the speech log-spectrum.
The mean of the speech log-spectrum posterior is used together with
∠ Y to create the enhanced speech signal using the inverse
STFT (ISTFT).

Apart from the speech log-spectrum, the non-linear KF also tracks
the reverberation log-spectrum, r_t, and the γ_t and β_t
reverberation parameters. The KF, as defined by the gray blocks in
Fig. <ref>,
has a speech KF prediction step, a reverberation KF prediction step
and a series of KF update steps. The reverberation KF is comprised
of the blocks “Reverberation KF prediction”, “KF Update” and
“γ_t, β_t KF Update”. These three blocks perform
joint denoising and dereverberation and estimate s_t|t and Σ_t|t^(s)
to enhance noisy reverberant speech.







The structure of the rest of this algorithm description section is
as follows. Sections III.A and III.B present the speech and reverberation
KF prediction steps, respectively. Section III.C describes the KF
update step and Sec. III.D the priors for the γ_t and β_t
parameters that are needed so that the KF (a) distinguishes between
speech and reverberation, and (b) does not diverge to non-realistic
T_60 and DRR estimates. Section III.E describes the unshaded
peripheral blocks in Fig. <ref>.



 §.§ The Speech KF Prediction Step


The speech KF prediction step is linear and is related to the “Speech
KF prediction”, “Decorrelate” and “Recorrelate” blocks in
Fig. <ref>.
The speech KF prediction step is described in <cit.> and in
<cit.> <cit.> and is based on conditional distributions
to model short-term dependencies. Decorrelation and recorrelation
of the speech KF state in (<ref>) are performed after and
before the speech KF prediction step, respectively. The decorrelation
and recorrelation operations in Fig. <ref>,
which are performed so that the non-linear KF update step can be applied,
perform vector-matrix and matrix-matrix multiplications for the speech
KF state mean and its covariance matrix, respectively, using B_t ∈ℝ^p × p
<cit.> <cit.>. The outputs of the “Decorrelate” block
are: (a) the first element of the speech KF state, and (b) the rest
elements of the speech KF state.

The KF prediction step propagates the first and second moments of
the speech KF state <cit.> <cit.>. Inter-frame linear
relationships are used for the speech KF prediction step that uses
AR modelling in the log-magnitude spectral domain. In the speech KF
prediction step, s_t is predicted as a linear combination
of s_t-1 using the speech AR coefficients that are obtained
from the “Speech AR(p)” block in Fig. <ref>,
which uses pre-cleaned speech as an input. After the speech KF prediction
step, s_t is correlated; we decorrelate the speech KF
state with a linear transformation (using B_t) to simplify
the KF update step and impose the observation constraint <cit.>.
The KF update step changes only the first element of the speech KF
state and after the KF update step, recorrelation is applied with
a linear transformation (using B_t^-1) to continue the
KF recursion.



 §.§ The Reverberation KF Prediction Step


The presented algorithm uses a KF prediction step for γ_t
and β_t that assumes that the variance of γ_t and β_t
increases over time, preserving their mean. The KF algorithm implements
a random-walk prediction step, performing the operations of

    γ_t|t-1^'=γ_t-1|t-1,     Σ_t|t-1^(γ^')=Σ_t-1|t-1^(γ)+Q_γ,
       β_t|t-1^'=β_t-1|t-1,     Σ_t|t-1^(β^')=Σ_t-1|t-1^(β)+Q_β

where Q_γ is a fixed error variance for γ and Q_β
is a fixed error variance for β. The values used for the prediction
error variances, Q_γ and Q_β, depend on the rate
at which the T_60 and the DRR are likely to change in a real
situation.

After the reverberation KF prediction step, the algorithm computes
and imposes priors on γ_t and β_t using Gaussian-Gaussian
multiplication. The internally computed priors for γ_t and
β_t in the “γ_t, β_t priors” block in Fig.
<ref>
are explained in Sec. III.D. After imposing the priors, the outputs
are γ_t|t-1, Σ_t|t-1^(γ) and β_t|t-1,
Σ_t|t-1^(β). We note that a prime diacritic, ',
is used in (7) and (8) to denote quantities before the priors.


The “Reverberation KF prediction” block in englishFig.
<ref>
estimates the first two moments of the prior distribution of the reverberation
spectral log-amplitude, i.e. r_t|t-1 and its variance.english
The algorithm performs a reverberation KF prediction step based on
the previous posterior of both speech and reverberation using the
signal model in (<ref>), where a is less than unity and this
makes the reverberation KF prediction step stable.


From (<ref>) and Fig. 2, the STFT-domain reverberation is the
sum of two components arising, respectively, from the reverberation
and speech components of the previous frame. The old reverberation,
δ_t, and the new reverberation, ϵ_t, are defined
in (<ref>) and (<ref>), respectively. The KF algorithm
calculates the prior distributions of these two components in the
log-amplitude spectral domain using δ_t|t-1 = γ_t|t-1 + r_t-1|t-1
and ϵ_t|t-1 = β_t|t-1 + s_t-1|t-1. These equations
are based on (<ref>) and (<ref>) with a common condition
added to all terms. Assuming that r_t-1 and s_t-1 are uncorrelated
with γ_t and β_t, respectively, britishthe
means and variances of the two Gaussian distributions are added. The
variances therefore add,

    Σ_t|t-1^(δ)=Σ_t|t-1^(γ)+Σ_t-1|t-1^(r),
       Σ_t|t-1^(ϵ)=Σ_t|t-1^(β)+Σ_t-1|t-1^(s).


As shown in Fig. <ref>,
the final operation of tbritishhe “Reverberation
KF prediction” block is to compute the prior distribution, r_t|t-1.
The addition in the cbritishomplex STFT domain
of two random variables in the log-spectral domain is modelled. The
reverberation log-amplitude spectrum is estimated by modelling the
addition in the STFT domain of two random variables in the log-amplitude
spectral domain. Given two disturbance sources, we combine them into
a single disturbance source.

From this point onwards in Sec. III.B, the time-frame subscript t|t-1
is omitted for clarity. For example, p(δ) is used instead
of p_δ_t|t-1(δ), which is defined in Sec. III.

From (<ref>), the reverberation component is the STFT-domain
sum of two elements arising, respectively, from the reverberation
and speech components in the previous frame. The log-amplitude spectral
domain distributions of these two elements, δ_t|t-1 and
ϵ_t|t-1, were calculated in the preceding paragraphs.
A two-dimensional Gaussian distribution is used for p(δ,ϵ) = p(δ) p(ϵ)
assuming independence between δ and ϵ. We assume
that the phase difference, η, between the two disturbance sources,
δ and ϵ, is uniformly distributed, i.e. η∼ U(-π, π),
and independent of their magnitudes. We write e^2r=e^2δ+e^2ϵ+2cos(η)e^δ+ϵ
and r=0.5log(e^2δ+e^2ϵ+2cos(η)e^δ+ϵ),
which takes account of η <cit.> <cit.>. Next, we
calculate

    𝔼{ r}    =∫_η=-π^πp(η) ∫_δ,ϵr p(δ,ϵ) dδ dϵ dη
        =∫_η=0^π1π ∫_δ,ϵ0.5log(e^2δ+e^2ϵ+2cos(η)e^δ+ϵ)
       × p(δ) p(ϵ) dδ dϵ dη,
    𝔼{ r^2}    =∫_η=0^π1π ∫_δ,ϵ(0.5log(e^2δ+e^2ϵ..
       ..+2cos(η)e^δ+ϵ))^2× p(δ) p(ϵ) dδ dϵ dη

where K_(δ,ϵ) sigma points are used to evaluate the
inner integral over (δ,ϵ) <cit.> and K_η
the outer integral over η.

Equations (<ref>) and (<ref>) estimate the first two
moments of the prior distribution of reverberation, r_t|t-1.
In this algorithm description section, we provide expressions for
either the second moment or the variance. We convert implicitly between
them using, for example, Σ_t^(r) = 𝔼{ r_t^2 }  - ( 𝔼{ r_t } )^2
for r_t.



 §.§ The Non-Linear KF Update Step


The KF algorithm decomposes the noisy reverberant observation, y_t,
into its component parts using distributions in the log-magnitude
spectral domain. britishThe decompositions are
based on Fig. 2british and the signal model in
(<ref>)british and (<ref>british).
The KF algorithm performs a series of low-dimensional operations
instead of a high-dimensional one in the KF update step.british
The adaptive KF algorithm propagates backwards through Fig. 2 and
decomposesbritish: (a) y_tbritish
into speech, s_t,british and into reverberation
and noise, z_tbritish, (b) the reverberation
and noise, z_tbritish, into reverberation,
r_tbritish, and into noise, n_tbritish,
and (c) the reverberation, r_tbritish, into
“old reverberation”, δ_t, and britishinto
british“new reverberation”, ϵ_tbritish.
The reverberation and noise log-spectrum, z_t, is a variable of
the “KF Update” block inbritish Fig. <ref>.

The KF algorithm in Fig. <ref>
uses the noisy reverberant observation, y_t, to first update s_t
and z_t and then update r_t. The posterior z_t|t is computed
in the “KF Update” block in Fig. <ref>.
In the proposed KF algorithm, the observed y_t affects z_t|t
directly and, in turn, z_t|t affects r_t|t. We hence divide
the observation update into two steps: (a) we use the log-spectrum
observation, y_t, to estimate the posterior distributions s_t|t
and z_t|t because Y_t = S_t + Z_t in (<ref>), and (b)
we use z_t|t as an “observation”
to obtain the posterior distributions r_t|t and n_t|t because
Z_t = R_t + N_t in (<ref>). In (b), we calculate an updated
version of r_t by using the posterior z_t|t as a KF observation
constraint. Hence, according to (a) and (b), the log-spectrum observation,
y_t, provides new information about s_t and r_t.

The sequence of operations involved in britishthe
“Reverberation KF prediction”, “KF update” and “γ_t,
β_t KF update” blocks are listed in Table <ref>.british
The “Reverberation KF prediction” block britishin
Fig. <ref>,
which was britishpresented in Sec. III.B, performs
the first five operations in Table <ref>. britishThe
“KF update” and “γ_t, β_t KF update” blocks
in Fig. <ref>british
perform the next seven operations in Table <ref>british,
i.e. steps 6-12. The bottom “z^-1” block in Fig. <ref>
britishperforms step 13. These 13 steps constitute
the dereverberation KF update step. The non-linear dereverberation
KF update step computes the first two moments of the posterior distributions
for most signal quantities and, moreover, includes the prediction
step as well for some quantities. Both means and variances are computed
for the tracked Gaussian signals; for clarity, the variances, such
as Σ_t|t^(r) for r_t, are not included in Table
<ref>.









[caption = The operations performed in the “Reverberation KF prediction”, “KF Update” and “γ_t, β_t KF Update” blocks in Fig. <ref>. Steps 7-12 perform specific signal decompositions propagating backwards through Fig. 2.,   
label = tab:ppssaaa,   
pos = hp,   
doinside=, width=.492]p0.2cm p7.8cm 
    Inputs: (a) s_t|t-1 from the speech KF prediction step, (b) n_t|- from external noise estimation, (c) y_t from observation, and (d) s_t-1|t from s_t|t in step 7 and B_t^-1 from Sec. III.A.
1:     γ_t|t-1 ← γ_t-1|t-1 from step 13 2:     β_t|t-1 ← β_t-1|t-1 from step 13 3:     δ_t|t-1 ← γ_t|t-1, r_t-1|t-1 from steps 1, 13 4:     ϵ_t|t-1 ← β_t|t-1, s_t-1|t-1 from steps 2, 13 5:     r_t|t-1 ← δ_t|t-1, ϵ_t|t-1 from steps 3, 4 6:     z_t|t-1 ← r_t|t-1, n_t|- from step 5 and input (b) 7:     s_t|t, z_t|t ← s_t|t-1, z_t|t-1, y_t from 6 and inputs (a), (c) 8:     r_t|t, [n_t|t]  ← r_t|t-1, n_t|-, z_t|t from 5, 7 and input (b) 9:     ϵ_t|t^' ← β_t|t-1, s_t-1|t from step 2 and input (d) 10:     δ_t|t, ϵ_t|t ← δ_t|t-1, ϵ_t|t^', r_t|t from steps 3, 8, 9 11:     γ_t|t, [r_t-1|t]  ← γ_t|t-1, r_t-1|t-1, δ_t|t from steps 1, 10, 13 12:     β_t|t, [s_t-1|t]  ← β_t|t-1, s_t-1|t, ϵ_t|t from 2, 10 and input (d)13:     r_t-1|t-1, γ_t-1|t-1, β_t-1|t-1 ← r_t|t, γ_t|t, β_t|t
 




For clarity, from this point onwards in Sec. III.C, the time subscript
is included only if it differs from t|t-1. Table <ref>
shows the time subscripts. We also denote n_t|- by n_t|t-1.

The KF algorithm computes the total distubance, z_t, from (<ref>)
in step 6 using similar equations to step 5. A two-dimensional Gaussian
distribution is used for p(r,n) and independence is assumed between
r and n. Hence, p(r,n) = p(r) p(n). The phase-sensitive KF
algorithm assumes that the phase difference, ζ, between the
two disturbance sources, r and n, is uniformly distributed and
independent of their magnitudes. From (<ref>), we write e^2z=e^2r+e^2n+2cos(ζ)e^r+n
and thus z=0.5log(e^2r+e^2n+2cos(ζ)e^r+n).
Next, using p(ζ)=1/2π, the first two moments of z
are given by

    𝔼{ z}    =∫_ζ=0^π1π∫_r,n0.5log(e^2r+e^2n+2cos(ζ)e^r+n)
       × p(r) p(n) dr dn dζ,
    𝔼{ z^2}    =∫_ζ=0^π1π∫_r,n(0.5log(e^2r+e^2n..
       ..+2cos(ζ)e^r+n))^2× p(r) p(n) dr dn dζ

where K_(r,n) sigma points are used to evaluate the inner integral
over (r,n) <cit.> and K_ζ the outer integral over
ζ.

The KF algorithm performs noise suppression with steps 6 and 7. Step
7 decomposes y_t into s_t and z_t, as shown in the signal
model in Fig. 2, estimating both s_t|t and z_t|t.

Step 7 performs the first signal decomposition, y_t into s_t
and z_t, when propagating backwards through the signal model in
Fig. 2. Step 7 applies the observation constraint, y_t, according
to (<ref>). As in <cit.> <cit.>, the variables are
first transformed according to (s,z,λ) ⇒ (u,y,λ)
where u = z-s and λ = ∠ S - ∠ Z. This variable
transformation is performed to allow the imposition of the scalar
KF observation, y_t.

The noisy reverberant log-amplitude spectrum, y, is given by y = s + 0.5 log(1 + exp(2(z-s)) + 2 cos(λ) exp(z-s)).
The KF update step assumes that ∠ Z is uniformly distributed,
∠ Z ∼ U(-π,π). Therefore, λ∼ U(-π, π).
The first two moments of the posteriors of s_t and z_t are computed
using

    𝔼{ s_t^m_1z_t^m_2 | y_0,y_1,…,y_t}
        =∫_λ=-π^π∫_u=-∞^∞s^m_1z^m_2 p(u,λ|y_t) du dλ
       ∝1|Δ|∫_λ=-π^πp(λ)∫_u=-∞^∞s^m_1z^m_2 p(s,z) du dλ
       ∝∫_λ=0^π∫_u=-∞^∞s^m_1z^m_2 p(s) p(z) du dλ

where the Jacobian determinant is Δ = -1 and the moment indexes,
m_1 and m_2, are integers, 0 ≤ m_1, m_2 ≤ 2. We denote
the variables for two moment indexes by m_1 and m_2.

The first two moments of the posterior distributions of s_t and
z_t are estimated. In (<ref>), the priors of the speech
and of the noise and late reverberation are assumed to be independent,
i.e. p(s,z) = p(s) p(z). In addition, in (<ref>), K_λ
weighted sigma points are used to evaluate the outer integral over
λ <cit.>.

Step 8 performs the second signal decomposition, z_t into r_t
and n_t, when propagating backwards through the proposed signal
model in Fig. 2. Step 8 decomposes z_t|t into r_t and n_t,
according to (<ref>), estimating both r_t|t and n_t|t.
britishStep 8 performs an integral over britishz_t
where the integrand is similar to step 7 and (<ref>) and to
the KF update step in <cit.>british <cit.>british.
Instead of a scalar observation, as in step 7, the observation in
step 8 is a distribution; step 8 performs an outer integral over the
observation distribution and the integrand is similar to step 7. Step
8 uses 𝔼{ r^m} = ∫_z_t𝔼{ r^m | z_t} dz_t
where m ∈ℕ_0 and 0 ≤ m ≤ 2british.
Step 7 computes a two-dimensional integral over the variables of (z-s)
and of the phase difference between S and Z, using y_t that
reduces the probability space from three to two dimensions. In step
8, the KF algorithm calculates a three-dimensional integral over:
(a) (r-n), (b) the phase difference between R and N, i.e.
ι = ∠ R - ∠ N, and (c) the posterior of z_t.
Assuming ∠ N ∼ U(-π,π), step 8 computes

    𝔼{ r_t^m_1n_t^m_2 | y_0,y_1,…,y_t}∝∫_z_t=-∞^∞p_z_t|t(z)
       ×∫_ι=0^π ∫_u'=-∞^∞r^m_1n^m_2 p(r) p(n) du' dι dz_t

where u' = (r-n) and ι∼ U(-π,π). In (<ref>),
K_ι weighted sigma points are used to evaluate the integral
over ι and K_z_t|t sigma points are used to evaluate
the integral over z_t <cit.>.


Ienglishn Table <ref>, the signals
in square brackets are calculated but are not used in the KF recursion.
In step 8, the posterior of n_t is estimated using (<ref>)
but is not used in the KF recursion.


Step 9 determines a preliminary estimate of the posterior distribution
of new reverberation component, ϵ_t in (<ref>). This
preliminary estimate, denoted ϵ_t|t^', combines
an updated estimate of the previous frame's speech, s_t-1 with
the prior estimate of the reverberation parameter, β_t. In
step 9,british two random variables in the log-amplitude
spectral domain are added; the means and variances of the two Gaussian
variables are: ϵ_t|t^' = β_t|t-1 + s_t-1|tbritish
and Σ_t|t^(ϵ^') = Σ_t|t-1^(β) + Σ_t-1|t^(s).
In this addition, β_t and s_t-1 are assumed to be independent.

Step 10 decomposes the reverberation, r_t, into a new reverberation
component, ϵ_t and an old decaying reverberation component,
δ_t, using (1), (4) and (5). Step 10 uses the prior distributions
for these two components from steps 3 and 9. britishStep
10 performs the same operation as step 8. In analogy to step 7, the
variable transformations in steps 8 and 10 are (r,n,ι) ⇒ (u',z,ι)
and (δ,ϵ,ξ) ⇒ (u”,r,ξ), respectively.
In step 10, the KF algorithm performs an integral over britishr
where the integrand is similar to the KF update step in <cit.>british.
Step 10 estimates the posterior distributions of δ and ϵ;
it estimates both δ_t|t and ϵ_t|t using δ_t|t-1
and ϵ_t|t'. Step 10 computes

    𝔼{δ_t^m_1ϵ_t^m_2 | y_0,y_1,…,y_t}∝∫_r_t=-∞^∞p_r_t|t(r)
       ×∫_ξ=0^π ∫_u”=-∞^∞δ^m_1 ϵ^m_2 p(δ) p_ϵ_t|t^'(ϵ) du” dξ dr_t

where u” = (δ - ϵ) and ξ∼ U(-π, π) is
the phase difference of the STFT-domain signals of δ and ϵ.
Using (<ref>), we decompose r_t|t into old reverberation
and new reverberation, estimating δ_t|t and ϵ_t|t.
In (<ref>), K_ξ sigma points are used to evaluate the
integral over ξ and K_r_t|t the integral over r_t.

Steps 10-12 perform the final signal decomposition, r_t into γ_t
and β_t, when propagating backwards through the signal model
in Fig. 2.british In steps 11 and 12, the KF algorithm
computes the first two moments of the posterior distributions of γ_t
and β_t. In step 11, the algorithm performs an integral over
britishδ_t using weighted sigma points
where the integrand models the britishaddition
of two random variables in the log-amplitude spectral domain, γ_t
and r_t-1. Likewise, in step 12, the algorithm performs an integral
over britishϵ_t using sigma points where
the integrand models the britishaddition of two
random variables in the log-amplitude spectral domain, β_t
and s_t-1. Steps 11 and 12 perform a linear KF update and impose
a straight line observation constraint because the signal model is
that γ_t and r_t-1 are additive and that β_t and
s_t-1 are additive, according to (<ref>) and (<ref>).

In step 11, the KF algorithm decomposes the old reverberation component,
δ_t, according to (<ref>) into the sum of a reverberation
parameter, γ_t, and the previous frame's reverberation, r_t-1.
We define x_1 ∼N(x_1;  μ_1, S_1)
where x_1 = (γ_t,   r_t-1,   w_t)^T ∈ℝ^3,
w_t is zero-mean Gaussian with variance Σ_t^(δ),
μ_1 = (γ_t|t-1,   r_t-1|t-1,   0)^T ∈ℝ^3
and S_1 = diag((Σ_t|t-1^(γ),  Σ_t-1|t-1^(r),  Σ_t|t^(δ))^T) ∈ℝ^3 × 3
where diag(j) is a diagonal matrix with the elements
of j on the main diagonal of the matrix. If x_1 ∼N(x_1; μ_1, S_1),
then

    (x_1|g^Tx_1=δ_t|t)    ∼N(x_1; (I-H_1g^T)μ_1+H_1δ_t|t,.
          .(I-H_1g^T)S_1)

where H_1 = S_1 g (g^T S_1 g)^-1
and g = (1,   1,   -1)^T ∈ℝ^3 <cit.>.
We compute (I - H_1g^T)μ_1+H_1 δ_t|t∈ℝ^3
and (I-H_1g^T)S_1 ∈ℝ^3 × 3
and set γ_t|t equal to the first element of this computed
mean and Σ_t|t^(γ) equal to the first element of
this variance.

Likewise, for step 12, we define x_2 ∼N(x_2;  μ_2, S_2)
where x_2 = (β_t,   s_t-1,   q_t)^T ∈ℝ^3,
q_t is zero-mean Gaussian with variance Σ_t^(ϵ),
μ_2 = (β_t|t-1,   s_t-1|t,   0)^T ∈ℝ^3
and, moreover, S_2 = diag((Σ_t|t-1^(β),  Σ_t-1|t^(s),  Σ_t|t^(ϵ))^T) ∈ℝ^3 × 3.
Next, we use

    (x_2|g^Tx_2=ϵ_t|t)    ∼N(x_2; (I-H_2g^T)μ_2+H_2ϵ_t|t,.
          .(I-H_2g^T)S_2)

where H_2 = S_2 g (g^T S_2 g)^-1
<cit.>. The KF algorithm computes (I - H_2g^T)μ_2+H_2 ϵ_t|t∈ℝ^3
and (I-H_2g^T)S_2 ∈ℝ^3 × 3
and sets β_t|t equal to the first element of this computed
mean and Σ_t|t^(β) equal to the first element of this
computed variance.

In steps 11 and 12, the presented KF algorithm computes γ_t|t,
Σ_t|t^(γ) and β_t|t, Σ_t|t^(β),
respectively. Steps 11 and 12 use p(γ_t, r_t-1 | δ_t|t + w_t)
and p(β_t, s_t-1 | ϵ_t|t + q_t), respectively.

Step 13 applies a one-frame delay, i.e. britishz^-1,
to continue the KF recursion britishas shown in
Fig. <ref>.
In summary, the 13 steps have the four main operations of: (a) step
5, (b) step 7, (c) step 8, and (d) step 11. The operations performed
in the other steps are either simple or identical to these operations.


According to the 13 stepsenglish in Table <ref>,
the γ_t and β_t reverberation parameters are affected
by the observed y_t through a series of operations that estimate
the first two moments of the posterior distributions. The estimate
of γ_t is affected by the previous estimate of β_t
because of the sequence of operations: β_t-1|t-1⇒β_t|t-1⇒ϵ_t|t^'⇒δ_t|t⇒γ_t|t.
The estimate of β_t depends on the previous estimate of γ_t
because of the sequence of operations: γ_t-1|t-1⇒γ_t|t-1⇒δ_t|t-1⇒ϵ_t|t⇒β_t|t.

The proposed 13 steps do not include the speech KF prediction step,
which is shown in englishFig. <ref>,
that is used to calculate s_t-1|t that is needed in steps 9 and
12. englishAfter step 7, according to englishFig.
<ref>
and Sec. III.B, recorrelation of the speech KF state is performed
with B_t^-1. Using s_t|t, from the recorrelation
operation, s_t-1|t is obtained. In steps 9 and 12, s_t-1|t
is used as a better estimate of s_t-1 than s_t-1|t-1.

In step 9, we note that two sub-indexes of t|t-1 and t-1|t give
rise to a sub-index of t|t: ϵ_t|t^' = f(β_t|t-1, s_t-1|t).
In step 9, we introduceenglish the notation ϵ_t|t^'
to avoid using the same symbol for two different posterior distributions,
englishϵ_t and englishϵ_t^'.




 §.§ The Priors for the Reverberation Parameters 


This section describes the priors for the γ_t and β_t
reverberation parameters, which are based on <cit.> <cit.>
and <cit.>. The priors for γ_t and β_t are imposed
using Gaussian-Gaussian multiplication; γ_t is modelled with
a Gaussian distribution and its internal prior is also a Gaussian.
Likewise, β_t is modelled with a Gaussian and its internal
prior is also a Gaussian.

The priors for γ_t and β_t are estimated from spectral
log-amplitude observations in the free decay region (FDR), which is
comprised of M_t consecutive frames with decreasing energy. We
define the look-ahead factor C and the frame index l=t-M_t+C+1,t-M_t+C+2,...,t+C.
The least squares (LS) fit to the FDR is found using r_l = θ_1 x_l + θ_2
where x_l is the time index in seconds and depends on L. The
parameters of the straight line are the slope, θ_1, and the
y-intercept, θ_2. For clarity, the frame subscript t is
omitted from θ_1 and θ_2. We define θ=(θ_1 θ_2)^T∈ℝ^2.
Its Gaussian distribution is

    N(θ;μ_θ,Σ_θ) ∝ |Σ_θ|^-0.5
       ×exp(-0.5(θ-μ_θ)^TΣ_θ^-1(θ-μ_θ))

where μ_θ∈ℝ^2 and Σ_θ∈ℝ^2 × 2.

The log-likelihood of θ is given by l(θ) = c_1 - 0.5 (θ - μ_θ)^T Σ_θ^-1 (θ - μ_θ) = c_2 - 0.5 ∑_τ=t-M_t+C+1^t+C ((θ_1 x_τ + θ_2 - r_τ)^T (Σ_τ^(r))^-1 (θ_1 x_τ + θ_2 - r_τ))
where c_1 and c_2 are constants.

We define the vectors x∈ℝ^M_t and r∈ℝ^M_t
from x_l and r_l, respectively. We define Σ_r∈ℝ^M_t × M_t
as the covariance matrix of r. The regression coefficients
as a Gaussian distribution are N(θ;μ_θ,Σ_θ)
where μ_θ=A^-1b and Σ_θ = A^-1,

    A=([ x^TΣ_r^-1x; 1^TΣ_r^-1x ] [ x^TΣ_r^-11; 1^TΣ_r^-11 ]),
       b=([ x^TΣ_r^-1r; 1^TΣ_r^-1r ])

where 1∈ℝ^M_t is a vector of ones. Figure
<ref> shows the log-likelihood of the Gaussian distribution
of θ with mean μ_θ and variance Σ_θ
when the correlation between θ_1 and θ_2 is strong.
We note that the maximum likelihood solution for the mean of θ_1,
when θ_2 is not considered, can be found in <cit.>.



The noisy reverberation, z_t, is the observation and, in this case,
our aim is to find r_t and Σ_t^(r). We denote the noise
floor in the power spectral domain by |N|^2. In a FDR, at high
RNRs and when |N|^2 is low, we have r_t ≈ z_t because:
z = r + 0.5log(1 + RNR^-1 + 2 cos(ζ)  RNR^-0.5),
which was also used in (<ref>) and (<ref>), where
RNR = exp(2(r-n)).

To solve the problem of finding r_t and Σ_t^(r) given
the observed z_t, we first employ a minimum MSE (MMSE) algorithm,
such as the traditional Log-MMSE estimator <cit.>, to
remove the noise and estimate the signal's spectral amplitude, as
in the end of Sec. 3.1 in <cit.>, and then perform a KF update
step and use z_t as the KF observation. For the linear KF <cit.>
<cit.>, when the observation noise covariance matrix tends to
zero at high RNRs, then the mean of the KF state goes to the KF observation
and the variance of the KF state goes to zero. Therefore, a solution
to finding r_t and Σ_t^(r) is to set r_t = z_t and
Σ_t^(r) equal to a small value, such as 1 dB^2.

The T_60 can be estimated from θ_1 using T_60 = - 3 log(10)/θ_1,
where θ_1 is computed from FDR data in the log-magnitude spectral
domain <cit.>. The KF algorithm, which models γ_t with
a Gaussian distribution, uses a Gaussian prior for γ with
mean L μ_θ_1 and variance L^2 Σ_θ_1 because
γ = L θ_1. Using γ = L θ_1, θ_1
is directly mapped to γ_t, without going via the T_60.

The algorithm estimates a prior for β directly, at the same
time as the prior for γ, rather than first estimating the
T_60 and the DRR. Figure <ref> depicts y_t, s_t
and r_t over time, at 1 kHz, when the T_60 is 0.68 s,
the DRR is -2.52 dB, the noise type is white and the SNR is 20
dB. In Fig. <ref>, s_t and r_t are the true speech
and reverberation powers; the latter is computed by convolving the
RIR without its first 30 ms with the speech signal <cit.>.
If we choose an appropriate FDR, then we can see that (<ref>)
is related to the FDR and the first frame of the FDR to the DRR. The
β reverberation parameter is computed from the difference between
the log-amplitude of the first frame and the y-intercept of the LS
fit, i.e. θ_2. For the LS fit to the FDR, we choose to not
use the first frame of the FDR. For the subtraction of two independent
random variables, their means subtract and their variances add; θ_2
has a mean and a variance and the log-amplitude of the first frame
of the FDR has a mean and a fixed small variance. The prior for β
is given by

    β=r_t-M_t+C+1-θ_2.


We denote the mean of the internally computed Gaussian γ_t
prior by γ_t|t-1^'' and its variance by Σ_t|t-1^(γ^'').
We use

    γ_t|t-1=γ_t|t-1^'Σ_t|t-1^(γ^'')+γ_t|t-1^''Σ_t|t-1^(γ^')Σ_t|t-1^(γ^')+Σ_t|t-1^(γ^''),
       Σ_t|t-1^(γ)=Σ_t|t-1^(γ^')Σ_t|t-1^(γ^'')Σ_t|t-1^(γ^')+Σ_t|t-1^(γ^'').

Likewise, we denote the mean of the Gaussian β_t prior by β_t|t-1^''
and its variance by Σ_t|t-1^(β^''). For
β_t, we use

    β_t|t-1=β_t|t-1^'Σ_t|t-1^(β^'')+β_t|t-1^''Σ_t|t-1^(β^')Σ_t|t-1^(β^')+Σ_t|t-1^(β^''),
       Σ_t|t-1^(β)=Σ_t|t-1^(β^')Σ_t|t-1^(β^'')Σ_t|t-1^(β^')+Σ_t|t-1^(β^'').


Assuming that the RIR is known, it is straightforward to compute the
T_60: compute the energy decay curve, plot it on a log scale
and estimate the T_60 from its slope. On the contrary, estimating
the T_60 blindly is not trivial. The KF algorithm estimates the
T_60, applying internally estimated priors using the decay rate
of the LS fit to the FDR <cit.>, so that the KF does not diverge
and does not treat reverberation as speech.



 §.§ The Peripheral Blocks of the Algorithm


This section describes the unshaded peripheral blocks of the KF algorithm
in Fig. <ref>.
The algorithm uses pre-cleaning before performing speech AR modelling,
as in <cit.> and <cit.>. Pre-cleaning has also been used
in <cit.> <cit.>, in <cit.> <cit.> and in
<cit.> <cit.>. The “Speech pre-cleaning” block in
Fig. <ref>
affects the “Speech KF prediction” block and not the observation
of the non-linear “KF Update” block, i.e. y_t, used in step
7.

For the “Noise power estimation” block, an external noise power
estimator, such as <cit.>, is used. External noise estimation
and log-normal noise power modelling, as in <cit.> <cit.>,
are used for n_t|-, which is then used in step 6 of Table I.

In summary, the KF algorithm tracks speech and reverberation in the
spectral log-magnitude domain along with γ_t and β_t,
as described in Fig. <ref>
and in the signal model in Fig. 2.



§ IMPLEMENTATION, TESTING AND VALIDATION


We use acoustic frames of length 32 ms and an acoustic frame increment
of L=8 ms in (2). We also use modulation frames of 64 ms and
a modulation frame increment of 8 ms <cit.> <cit.>.
In Fig. <ref>,
for speech amplitude spectrum pre-cleaning, we use the Log-MMSE estimator
<cit.> followed by the WPE dereverberation algorithm
<cit.> <cit.>. The dimension of the speech KF state, s_t,
in (<ref>) is p=2. In the “Noise power estimation”
block of Fig. <ref>,
we use external noise power estimation from <cit.> using the
implementation in <cit.>.

The outer integrals in Secs. III.B and III.C are performed using sigma
points, as in <cit.>. The numbers of sigma points used in (<ref>)-(<ref>)
are K_(δ,ϵ) = K_(r,n) = K_z_t|t = K_r_t|t = 3
and K_η = K_ζ = K_λ = K_ι = K_ξ = 6.
For the latter cases, for ξ, the sigma points are at ξ = ((1:K_ξ)-0.5) π/K_ξ
and the weights are all equal to 1/K_ξ <cit.>
<cit.>. With this choice of sigma points for ξ, the integral
will be exact for an integrand that is a sum of terms of the form
cos(n ξ) for 0 ≤ n ≤ 2K_ξ+1.

In Sec. III.D, for the γ_t and β_t priors, the look-ahead
factor is C=3 frames. For the FDR that is comprised of frames with
decreasing energy, M_t is computed in every frame.



Figures 6(a)-(c) illustrate s_t, z_t and r_t at 1 kHz
over time, as Fig. <ref>
in <cit.>. Figure 6(a) shows the predicted and true speech power,
i.e. s_t|t and the true s_t, and Fig. 6(b) the predicted and
true noise and reverberation power, i.e. z_t|t and the true z_t.
The correlation coefficient between the true z_t and z_t|t
is 0.89. Figure 6(c) shows the predicted and true reverberation
power, i.e. r_t|t and the true r_t. The correlation coefficient
between the true r_t and r_t|t is 0.77 (<0.89). The
noise type is white, the SNR is 20 dB, T_60 = 0.61 s and DRR = -1.74
dB.

The ordering of the graphs in Figs. 6(a)-(c) matches the ordering
of the signal decompositions in the KF algorithm, in Table <ref>.
The noisy reverberant observation is first decomposed into s_t
and z_t with step 7 in Table <ref>; Figs. 6(a) and
(b) illustrate s_t|t and z_t|t, respectively, over time.
Then, z_t|t is decomposed into r_t and n_t with step 8;
Fig. 6(c) depicts r_t|t over time.



The T_60 and DRR estimates should converge to their true constant
values when the talker and the microphone are not moving and the frequency
variations in the reflection coefficients are not modelled <cit.>
<cit.>. Internally estimated priors for γ_t and β_t
make the reverberation parameters converge over time. Figure 6(d)
shows the predicted and true T_60 and DRR reverberation parameters
over time. The dashed lines are the signals' standard deviations,
computed from γ_t and β_t.



§ EVALUATION OF THE ALGORITHM


The proposed KF algorithm is evaluated in terms of the perceptual
evaluation of speech quality (PESQ) <cit.>, the cepstral
distance (CD) spectral divergence metric <cit.>, the
reverberation decay tail (RDT) dereverberation metric <cit.>
and the STOI speech intelligibility metric <cit.> <cit.>.
The ideal values of CD and RDT, which have been also used in <cit.>,
are zero. For evaluation, the TIMIT database <cit.>, sampled
at 16 kHz, and the RSG-10 noise database <cit.> are used.
From the TIMIT database, 52 clean speech utterances are chosen.
We use artificially-created reverberation with the image method <cit.>
using the implementation in <cit.>. The wall reflection coefficient
is adjusted to vary the T_60 and hence also the DRR. The KF algorithm
is evaluated with noisy reverberant speech signals at various SNRs,
from 5 to 20 dB, using random noise segments and the noise types
of white, babble and factory.





The KF algorithm is compared to the SPENDRED algorithm <cit.>
<cit.>, which jointly performs blind denoising and dereverberation,
and to algorithms that sequentially perform denoising and dereverberation,
specifically to the Log-MMSE estimator <cit.> followed
by the WPE algorithm <cit.> <cit.> and to the Log-MMSE
estimator followed by an inverse-filter (IF) dereverberation method
that is based on T_60 and DRR estimates <cit.>. For the
IF, we blindly estimate the T_60 and the DRR for the entire speech
utterance using the algorithm in <cit.> <cit.>, whose
implementation was generously provided by the author. As found in
the Acoustic Characterisation of Environments (ACE) challenge <cit.>,
the T_60 estimator in <cit.> had the best performance of
the examined T_60 estimators.

For evaluation, the examined acoustic conditions are shown in Table
<ref>. Two different rooms with dimensions 5 × 4 × 4
m and 10 × 7 × 3 m are used. The distance between the
microphone and the talker is 1.5 m. Table <ref>
is sorted with respect to the T_60 and the DRR from A to I and
from J to V. We plot the improvement in PESQ, ΔPESQ, against
the index of the acoustic environment to evaluate the KF algorithm.
The top axis shows the raw PESQ of the noisy reverberant speech and
is monotonic from A to I and from J to V. Similar graphs are also
plotted for the other metrics. The ordering of the legends matches
that of the algorithms at low T_60 values.

Table <ref> also shows the a and b reverberation
values of the RIRs and we note that one of the baselines that we use,
i.e. the SPENDRED algorithm <cit.> <cit.>, assumes that
b ≤ 1. This is valid in Table <ref> because
b << 1 and its range is small.





 §.§ Overall Performance Against the T_60 and the DRR


This section investigates the overall performance of the KF algorithm
against the T_60 and the DRR. Figure 7 examines: (a) the ΔPESQ,
where higher scores signify better speech quality, (b) the ΔCD,
where lower values signify better speech quality, (c) the ΔRDT,
where lower values signify better dereverberation, and (d) the ΔSTOI,
where higher scores signify better intelligibility. Figure 7 first
presents the results that are related to speech quality, in (a) and
(b). Then, it shows the dereverberation results, in (c), and the speech
intelligibility results, in (d). The graphs are against the index
of the acoustic environment and, hence, against the T_60 and
the DRR. The average over the noise types of white, babble and factory
at 20 dB SNR is shown. We note that noisy reverberant speech with
T_60≈ 0.6 s has a raw PESQ score of 2 in Fig. 7(a).

In Fig. 7(a), the top axis is monotonic from A to I and from J to
V and there is no transition from index I to index J.





From Fig. 7(a), in terms of PESQ, the proposed algorithm consistently
yields improved speech quality in challenging environments compared
to the examined baselines. Compared to the unprocessed noisy reverberant
speech, the algorithm shows improved performance for all the examined
T_60 range from 0.18 to 1.05 s and DRR range from 8.43
to -3.33 dB. Compared to the unprocessed speech, for a T_60
of 0.3 and 1 s, the algorithm has a ΔPESQ of 0.35
and 0.25, respectively, for the SNR of 20 dB averaged over the
examined noises.

From Fig. 7(b), in terms of CD, the KF algorithm yields improved speech
quality in acoustic environments with a T_60 from 0.18 to
1.05 s. The KF algorithm shows a deteriorating CD improvement with
increasing T_60. Compared to the unprocessed signal, for a T_60
of 0.3 and 1 s, the algorithm has a CD improvement of approximately
-1.1 and -0.6, respectively, for the SNR of 20 dB averaged
over the tested noise types.

From Fig. 7(c), we observe that the KF algorithm yields improved speech
dereverberation in challenging acoustic conditions and that the ΔRDT
improves with increasing T_60 values. From the raw RDT of noisy
reverberant speech in Fig. 7(c), we observe that the indexes A and
J have very low reverberation. The indexes A and J have raw RDT scores
0.4 and 0.7, respectively. The indexes G-I and T-V have high
reverberation with a high raw RDT score. The proposed algorithm in
these high reverberation cases achieves a large RDT improvement decreasing
the RDT metric from approximately 2.6 to 0.8.

From Fig. 7(d), the KF algorithm shows marginally improved STOI performance
for the examined T_60 range compared to the unprocessed speech
and to the examined baselines. For a T_60 of 0.3 and 1
s, the algorithm has a ΔSTOI of 0.03 and 0.07, respectively.
With increasing T_60, the ΔSTOI scores slightly increase.
The baselines have negative ΔSTOI.



We use the noise types of white, babble and factory at 10 dB SNR
and we obtain Fig. 8. Figure 8 shows similar graphs to Fig. 7 but
for a SNR of 10 dB. For all the examined T_60 values, the
KF algorithm shows a consistent improvement in the evaluation metrics
of PESQ, CD, RDT and STOI.

We now use the noise types of white, babble and factory at 5 dB
SNR and we obtain Fig. 9. Figure 9 shows similar graphs to Figs. 7
and 8. For all the examined SNRs, the KF algorithm shows a consistent
improvement in the examined metrics, depending on the T_60. From
Fig. 9(a), we observe that the ΔPESQ scores of the KF algorithm
at 5 dB SNR are similar to the ΔPESQ scores at 10 dB
SNR in Fig. 8(a).

The KF algorithm shows improved PESQ performance for the examined
SNRs from 5 dB to 20 dB. The algorithm has better speech quality
performance compared to the baselines for all the examined T_60
values. Compared to the unprocessed noisy reverberant speech, for
the SNR of 5 dB, the algorithm has a ΔPESQ of about 0.45
for a T_60 of 0.6 s and a raw PESQ of 1.6. For the SNR
of 20 dB, the algorithm has a ΔPESQ of approximately 0.3
for a T_60 of 0.6 s and a raw PESQ of 2. Comparing Figs.
7-9, we observe that the raw PESQ decreases with decreasing SNR while
the ΔPESQ increases.

The KF algorithm yields improved dereverberation performance in terms
of RDT in adverse conditions in Fig. 9(c). The algorithm's performance
improves with decreasing SNR. We observe that the indexes A-C and
J-K have a high raw RDT despite that they have low reverberation at
5 dB SNR, when the effect of noise is higher than that of reverberation.

The results of the log-likelihood ratio (LLR) speech quality metric
<cit.> <cit.>, which is used as the main evaluation
metric in the dereverberation algorithm in <cit.>, resemble
the results of CD in Fig. 10. Both CD and LLR have been used in the
REVERB challenge <cit.> <cit.>. It was found in <cit.>
that LLR correlates well with speech quality although slightly less
well than PESQ. Figure 10 examines ΔLLR when white noise is
used at 20 and 10 dB SNR. From Fig. 10, we observe that decreasing
the SNR from 20 to 10 dB increases the raw LLR of noisy reverberant
speech and deteriorates the ΔLLR.



 §.§ Overall Performance Against the SNR


This section investigates the overall performance of the KF algorithm
against the SNR. Figure 11 presents the algorithm's PESQ performance
compared to the baselines. Figures 11(a) and 11(b) depict the PESQ
improvement, ΔPESQ, against the SNR for white noise when:
(a) T_60 = 0.40 s and DRR = 0.17 dB, which is case L
in Table <ref>, and (b) T_60 = 0.73 s and DRR = -2.1
dB, which is case R in Table <ref>. The indexes
L and R were chosen because they are not extreme cases and they have
different T_60 and DRR values. The ordering of the legends in
Fig. 11 matches that of the algorithms at high SNRs.



Figures 11(c) and 11(d) show the ΔPESQ against the SNR for
babble and factory noises, for cases L and R in Table <ref>
respectively, for the KF algorithm compared to the baselines. Babble
noise is used for the solid lines, the top axis and the legends while
factory noise is used for the dashed lines. The KF algorithm has a
ΔPESQ that is dependent on the SNR and noise type and is,
most of the times, increasing with decreasing SNRs. Compared to the
unprocessed speech in Fig. 11(c), the algorithm has a ΔPESQ
of about 0.35 for factory noise for SNRs from 5 to 20 dB,
while decreasing for lower SNRs.

Figure 12 examines the RDT improvement, ΔRDT, against the
SNR for white noise, for cases L and R in Table <ref>,
in (a) and (b). Figure 12 also depicts the ΔRDT against the
SNR for babble and factory noises, for cases L and R in Table <ref>,
in (c) and (d). Babble noise is used for the solid lines, the top
axis and the legends. Factory noise is used for the dashed lines.
Figures 12(a)-(d) show that the ΔRDT improves and the raw
RDT of noisy reverberant speech increases with decreasing SNR. The
ΔRDT scores are better for stationary white noise than for
non-stationary factory and babble noises.

Figure 13 shows the ΔPESQ and the ΔRDT against the
SNR, from 5 to 20 dB, for the KF algorithm for white, factory
and babble noises. Figure 13 examines the acoustic conditions J to
V that correspond to the room 10 × 7 × 3 m in Table <ref>.
Figure 13 differs from Figs. 11 and 12 in presenting the indexes J
to V and not only L and R. From Fig. 13(a), we observe that the ΔPESQ
decreases as the T_60 increases. The ΔPESQ improves with
decreasing SNR and the ΔPESQ is higher for white noise than
for factory and babble noises. In Fig. 13(b), the ΔRDT improves
as the T_60 increases and, moreover, the ΔRDT improves
with decreasing SNR.

In summary, in this evaluation section, we have tested the KF algorithm
in different SNRs, noise types and acoustic conditions with different
T_60 and DRR values. The algorithm is effective in enhancing
distorted speech, decomposing noisy reverberant speech into speech,
reverberation and noise. Regarding noise robustness, Figs. 7-13 show
that the proposed KF algorithm achieves a significant performance
gain over different noise types and SNRs, compared to the unprocessed
noisy reverberant speech and to the examined baselines.



§ CONCLUSION


In this paper, we present a monaural speech enhancement algorithm
based on Kalman filtering in the log-magnitude spectral domain to
blindly suppress noise and reverberation while accounting for inter-frame
speech dynamics. The first two moments of the posterior distribution
of the speech log-magnitude spectrum are estimated in noisy reverberant
environments using a model that adaptively updates the T_60 and
DRR reverberation parameters. The non-linear KF algorithm updates
and tracks the two reverberation parameters of γ_t and β_t
to further improve the estimation of the speech log-magnitude spectrum.
In this paper, we show by means of theoretical and experimental analyses
that Kalman filtering can be performed in challenging conditions by
performing the proposed signal decompositions into speech, reverberation
and noise, propagating backwards through the proposed signal model.
Experimental results using instrumental measures show improved performance
against both the T_60 and the SNR compared to the unprocessed
noisy reverberant speech and to alternative competing techniques that
perform blind denoising and dereverberation either in concatenation
or jointly.




IEEEtran
tocsection



