Toward Musical-Noise-Free Blind Speech
Extraction: Concept and Its Applications
Ryoichi Miyazaki∗ , Hiroshi Saruwatari∗ , Satoshi Nakamura∗ , Kiyohiro Shikano∗ ,Kazunobu Kondo† , Jonathan Blanchette‡ , and Martin Bouchard‡∗ Graduate School of Information Science, Nara Institute of Science and Technology, Japan

E-mail: ryoichi-m@is.naist.jp† Corporate Research & Development Center, Yamaha Corporation, Japan‡ School of Information Technology and Engineering, University of Ottawa, Canada

Abstract—In this paper, we review a blind musical-noise-freespeech extraction method using a microphone array that canbe applied to nonstationary noise. In our previous study, it was
found that optimized iterative spectral subtraction (SS) resultsin speech enhancement with almost no musical noise generation,
but this method is valid only for stationary noise. The proposed
method consists of iterative blind dynamic noise estimation by,e.g., ICA or multichannel Wiener filtering, and musical-noisefree speech extraction by modified iterative SS, where multiple
iterative SS is applied to each channel while maintaining themultichannel property reused for the dynamic noise estimators.
Also, related to the proposed method, we discuss the justificationof applying ICA to such signals nonlinearly distorted by SS. From
objective and subjective evaluations simulating real-world handsfree speech communication system, we reveal that the proposed
method outperforms the conventional methods.

I. I NTRODUCTION
In the past few decades, many applications of speechcommunication systems, such as hearing aids and mobile
phones, have been investigated. It is, however, well known
that these systems always suffer from the deterioration of
speech quality under adverse noise conditions, and thus noise
reduction is a problem requiring urgent attention. Spectral
subtraction is a commonly used noise reduction method thathas high noise reduction performance with low computational
complexity [1], [2], [3], [4], [5]. However, in this method,
artificial distortion, so-called musical noise, arises owing tononlinear signal processing, leading to a serious deteriorationof sound quality.
To achieve high-quality noise reduction with low musicalnoise, an iterative SS method has been proposed [6], [7], [8].
This method is performed through signal processing in which
weak SS processes are iteratively applied to the input signal.
Also, some of the authors have reported the very interesting
phenomenon that this method with appropriate parameters
gives equilibrium behavior in the growth of higher-order
statistics with increasing number of iterations [9]. This meansthat almost no musical noise is generated even with high noise
reduction, which is one of the most desirable properties ofsingle-channel nonlinear noise reduction methods. Following
this finding, the authors have derived the optimal parameters satisfying the no musical noise generation condition by
analysis based on higher-order statistics. We have defined thismethod as musical-noise-free speech enhancement, where nomusical noise is generated even for a high SNR in iterative
SS [10].
In conventional iterative SS, however, it is assumed thatthe input noise signal is stationary, meaning that we canestimate the expectation of noise power spectral density froma time-frequency period of a signal that contains only noise.
In contrast, under real-world acoustical environments, e.g.,a nonstationary noise field, although it is necessary to dynamically estimate noise, this is very difficult. Therefore inthis paper, first, we propose a new iterative signal extraction
method using a microphone array that can be applied to
nonstationary noise [11], [12]. Our proposed method consistsof iterative blind dynamic noise estimation by independent
component analysis (ICA) [13], [14] and musical-noise-freespeech extraction by modified iterative SS, where multiple
iterative SS is applied to each channel while maintaining themulti-channel property reused for ICA.
Secondly, related to the proposed method, we discuss thejustification of applying ICA to such signals nonlinearly
distorted by SS. We theoretically clarify that the degradation
in ICA-based noise estimation obeys an amplitude variationin room transfer functions between the target user and microphones. Next, to reduce speech distortion, we introduce achannel selection strategy into ICA, where we automaticallychoose less varied inputs to maintain a high accuracy of thenoise estimation. Furthermore, we introduce a time-variant
noise PSD estimator [15] instead of ICA for improvement of
the noise estimation accuracy. From objective and subjective
evaluations, we reveal that the proposed method outperforms
the conventional method.
The rest of the paper is organized as follows. In Sect. II, we
describe related works on spectral subtraction and the musicalnoise metric. In Sect. III, new musical-noise-free blind speech
extraction method is proposed. In Sect. IV, an improvement
scheme for poor noise estimation is presented. In Sect. V,
objective and subjective evaluations are described. Following
a discussion on the results of the experiments, we present ourconclusions in Sect. VI.

β

Noise

|X(f,τ)|Short-time
FFT
Speech

1st iteration+ηSpectralsubtraction

Noise PSD
estimation E[|N 21 |]-C. Modeling of input signal
In this paper, we assume that the input signal x in the power
spectral domain is modeled using the gamma distribution as
P (x) =+
2nd iteration

Noise PSD
estimation E[|N 22 |]...
Final iteration

Spectralsubtraction-xα−1
exp(−x/θ),Γ(α)θα(2)where x ≥ 0, α > 0, and θ > 0. Here, α is the shapeparameter, θ is the scale parameter,and Γ(α) is the gamma∫∞function, defined as Γ(α) = 0 tα−1 exp(−t)dt.
D. Mathematical metric of musical noise generation via
higher-order statistics for non-iterative SS [16]+Spectralsubtraction Output

Noise PSD
estimation E[|NF2 |]-Fig. 1. Block diagram of iterative SS.

II. R ELATED WORKSIn this study, we apply the kurtosis ratio to a noise-only
time-frequency period of the subject signal for the assessmentof musical noise [16]. This measure is defined as
kurtosis ratio = kurtproc /kurtorg ,where kurtproc is the kurtosis of the processed signal andkurtorg is the kurtosis of the observed signal. Kurtosis isdefined as
kurt = µ4 /µ22 ,A. Conventional non-iterative SS [2]We apply short-time Fourier analysis to the observed signal,
which is a mixture of target speech and noise, to obtainthe time-frequency signal. We formulate conventional noniterative SS [2] in the time-frequency domain as follows: √ |x(f, τ )|2 − βE[|N |2 ] exp(jarg(x(f, τ )))y(f, τ ) =(if |x(f, τ )|2 > βE[|N |2 ]),ηx(f, τ ) (otherwise),(1)
where y(f, τ ) is the enhanced target speech signal, x(f, τ )is the observed signal, f denotes the frequency subband, τis the frame index, β is the oversubtraction parameter, andη is the flooring parameter. Here, E[|N |2 ] is the expectationof the random variable |N |2 corresponding to the noise power
spectra. In practice, we can approximate E[|N |2 ] by averaging
the observed noise power spectra |n(f, τ )|2 in the first Ksample frames, where we assume the absence of speech inthis period and noise stationarity. However, this often requires
high-accuracy voice activity detection.
B. Iterative SS [6], [7], [8]In an attempt to achieve high-quality noise reduction withlow musical noise, an improved method based on iterative SS
was proposed in previous studies [6], [7], [8]. This method is
performed through signal processing, in which the following
weak SS processes are recursively applied to the noise signal(see Fig. 1). (I) The average power spectrum of the input
noise is estimated. (II) The estimated noise prototype is thensubtracted from the input with the parameters specifically set
for weak subtraction, e.g., a large flooring parameter η anda small subtraction parameter β. (III) We then return to step(I) and substitute the resultant output (partially noise reduced
signal) for the input signal.

(3)

where µm is the mth-order moment, given by∫ ∞µm =xm P (x)dx,

(4)

(5)

0

and P (x) is the probability density function (p.d.f.) of a powerspectral-domain signal x. A kurtosis ratio of unity corresponds
to no musical noise. This measure increases as the amount ofgenerated musical noise increases.
The mth-order moment after SS, µm , is given by [9]µm = θnm M(αn , β, η, m),(6)

where θn is the noise scale parameter, αn is the noise shapeparameter, andM(αn , β, η, m) =S(αn , β, η) + η 2m F(αn , β, η),(7)m∑
Γ(m+1)Γ(αn +m−l, βαn )S(αn , β, m) =(−βαn )l,Γ(αn )Γ(l+1)Γ(m−l+1)l=0(8)γ(αn +m, βαn )F(αn , β, m) =.
(9)Γ(αn )
Γ(b, a) and γ(b, a) are the upper and∫ lower incomplete∞
gamma functions defined as Γ(b, a) = b ta−1 exp(−t)dt∫ b a−1and γ(b, a) = 0 t
exp(−t)dt, respectively. From (4), (6),and (7), the kurtosis after SS can be expressed as
kurt =M(αn , β, η, 4).
M2 (αn , β, η, 2)(10)

Using (3) and (10), we also express the kurtosis ratio as
kurtosis ratio =M(αn , β, η, 4)/M2 (αn , β, η, 2).
M(αn , 0, 0, 4)/M2 (αn , 0, 0, 2)(11)

Also, as a measure of noise reduction performance, the noise
reduction rate (NRR), the output signal-to-noise ratio (SNR)

Harmful

12

Kurtosis ratio
NaturalIII. P ROPOSED M ETHOD : E XTENSION TO M ICROPHONE
A RRAY S IGNAL P ROCESSING

non-iterative
iterative (η = 0.9)

108

A. Conventional blind spatial subtraction array

6

In the previous section, we assumed that the input noise signal is stationary, meaning that we can estimate the expectationof a noise signal from a time-frequency period of a signal that
contains only noise, i.e., speech absence. However, in actual
environments, e.g., a nonstationary noise field, it is necessary
to dynamically estimate the noise power spectral density.
To solve this problem, we previously proposed blind spatial
subtraction array (BSSA) [17], which involves accurate noise
estimation by ICA followed by a speech extraction procedure
based on SS (see Fig. 3). BSSA improves the noise reduction
performance, particularly in the presence of both of diffuse andnonstationary noises; thus, almost all the environmental noise
can be dealt with. However, BSSA always suffers from musicalnoise owing to SS. In addition, the output signal of BSSA
degenerates to a monaural (not multi-channel) signal, meaningthat ICA cannot be reapplied; thus, we cannot iteratively
estimate the noise power spectra. Therefore, it is impossible
to directly apply iterative SS to the conventional BSSA.

4
2
0
0

2

4

6

8

10

12

Noise reduction rate [dB]Fig. 2. Relation between NRR and kurtosis ratio obtained from theoreticalanalysis for Gaussian noise case.

minus the input SNR in dB, can be given in terms of a 1storder moment as [9]NRR = 10 log10αn.
M(αn , β, η, 1)(12)

E. Musical-noise-free speech enhancement [10]In [10], we have proposed musical-noise-free noise reduction, where no musical noise is generated even for a high
SNR in iterative SS. In the study, first, some of the authors
discovered an interesting phenomenon that the kurtosis ratio
sometimes does not change even after SS via mathematicalanalysis based on (11) [9]. This indicates that the kurtosis ratio can be maintained at unity even after iteratively
applying SS to improve the NRR, and thus no musicalnoise is generated. Following this finding, the authors havederived the optimal parameters satisfying the musical-noisefree condition by finding a fixed-point status in the kurtosis ratio, i.e., by solving M(αn , 0, 0, 4)/M2 (αn , 0, 0, 2) =M(αn , β, η, 4)/M2 (αn , β, η, 2) [10]. Given the noise shapeparameter αn , we can choose combinations of the oversubtraction parameter β and the flooring parameter η thatsimultaneously satisfy the musical-noise-free condition using
the following equation;

B. Iterative blind spatial subtraction array [11]In this section, we propose a new multi-iterative blind signal
extraction method integrating iterative blind noise estimation
by ICA and iterative noise reduction by SS. As mentioned
previously, the conventional BSSA cannot iteratively andaccurately estimate noise by ICA because the conventional
BSSA performs a delay and sum (DS) operation before SS.
To solve this problem, we propose a new BSSA structurethat performs multiple independent SS in each channel before
DS; we call this structure channel-wise SS [18], [19], [20].
Using this structure, we can equalize the number of channelsof the observed signal to that of the signals after channelwise SS. Therefore, we can iteratively apply noise estimation
by ICA and speech extraction by SS (see Fig. 4). Also,
the advantage of the proposed structure is that ICA has thepossibility of adaptively estimating the distorted wavefront ofa speech signal to some extent even after SS, because ICA
is a blind signal identification method that does not requireη 4 ={F(αn , β, 4)(αn +1)αn −F 2 (αn , β, 2)(αn +3)(αn +2)}−1 knowledge of the target signal direction. Details of this issue[
will be discussed in Sect. III-C. Hereafter, we refer to this
proposed BSSA as iterative BSSA.
S(αn , β, 2)F(αn , β, 2)(αn +3)(αn +2)
We conduct iterative BSSA in the following manner, where[2
thesuperscript [i] represents the value in the ith iteration of± {S(αn , β, 2)F(αn , β, 2)(αn +3)(αn +2)}
SS
(initiallyi = 0).
{}− F(αn , β, 4)(αn +1)αn −F 2 (αn , β, 2)(αn +3)(αn +2)
The observed signal vector of the K-channel array
] (I)
] 12
{
}
in the time-frequency domain, x[0] (f, τ ), is given by
S(αn , β, 4)(αn +1)αn −S 2 (αn , β, 2)(αn +3)(αn +2).
x[0] (f, τ ) = h(f )s(f, τ ) + n(f, τ ),
(14)
(13)
where h(f ) = [h1 (f ), h2 (f ) . . . , hK (f )]T is a colFigure 2 shows an example of the kurtosis ratio in optimized
umn vector of the transfer functions from the target
iterative SS, where Gaussian noise is assumed. We can confirm
signal position to each microphone, s(f, τ ) is thethe flat trace of the kurtosis, indicating no musical noise
target speech signal, and n(f, τ ) is a column vector
generation.
of the additive noise.

(II)Next, we perform signal separation using ICA as [13][i]o (f, τ ) =W ICA (f )x[i] (f, τ ),(15)[i][p+1][i][i]
H
W ICA (f ) =µ[I − ⟨φ(o (f, τ ))(o (f, τ )) ⟩τ ][i][p][i][p]· W ICA (f ) + W ICA (f ),
(16)

F
F
T

User’sspeech

[i]

θU

Phase
compensation

Noise

θU
FDICA

+

Σ

Spectralsubtraction-

θU

0 Projection
back

θU

Σ

[i][p]

where W ICA (f ) is a demixing matrix, µ is the stepsize parameter, [p] is used to express the value of thepth step in the ICA iterations, I is the identity matrix,⟨·⟩τ denotes a time-averaging operator, and φ(·) is
an appropriate nonlinear vector function. Then, weconstruct a noise-only vector,
[i]

[i]

[i]

onoise (f, τ ) =[o1 (f, τ ), . . . , oU −1 , 0,
[i]

[i]

oU +1 (f, τ ), . . . , oK (f, τ )]T ,Fig. 3. Block diagram of conventional BSSA [17].

the microphone position, c is the sound velocity, andθU is the estimated direction of arrival of the target
speech. Moreover, [A]lj represents the entry of A in
the lth row and jth column.

(17)
C. Accuracy of wavefront estimated by ICA after SS

where U is the signal number for speech, and weapply the projection back operation to remove theambiguity of the amplitude and construct the estimated noise signal, z [i] (f, τ ), as

(III)

(IV)

In this subsection, we discuss the accuracy of the estimatednoise signal in each iteration of iterative BSSA. In actual
environments, not only point-source noise but also non-pointsource (e.g., diffuse) noise often exists. It is known that ICA[i][i]z [i] (f, τ ) = W ICA (f )−1 onoise (f, τ ).
(18) is proficient in noise estimation rather than speech estimation
under such a noise condition [17]. This is because the target
Next, we perform SS independently in each in- speech can be regarded as a point-source signal (thus, theput channel and derive the multiple target-speech- wavefront is static in each subband) and ICA acts as anenhanced signals. This procedure can be given byeffective blocking filter of the speech wavefront even in a time[i+1]invariant manner, resulting in good noise estimation. However,xk (f, τ ) =in iterative BSSA, we should address the inherent question of
 √
[i]
[i]
[i]

 |xk (f, τ )|2 − β|zk (f, τ )|2 exp(j arg(xk (f, τ )))whether the distorted speech wavefront after nonlinear noise[i][i]reduction such as SS can be blocked by ICA or not; thus, the(if |xk (f, τ )|2 > β|zk (f, τ )|2 ),
speech component after channel-wise SS can become a point[i]ηxk (f, τ ) (otherwise),source again or not.
(19)
Hereafter, we quantify the degree of point-source-likeness
[i+1]
for
SS-applied speech signals. For convenience of discussion,
where xk (f, τ ) is the target-speech-enhanced siga
simpletwo-channel array model is assumed. First, we define
nal obtained by SS at a specific channel k. Then
[i+1]
the
speechcomponent in each channel after channel-wise SS
we return to step (II) with x(f, τ ). When weasobtain sufficient noise reduction performance, go tostep (IV).
ŝ1 (f, τ ) = h1 (f )s(f, τ ) + ∆s1 (f, τ ),(24)
Finally, we obtain the resultant target-speechŝ2 (f, τ ) = h2 (f )s(f, τ ) + ∆s2 (f, τ ),(25)
enhanced signal by applying DS to x[∗] (f, τ ), where∗ is the number of iterations after which sufficient where s(f, τ ) is the original point-source speech signal,
noise reduction performance is obtained. This proceŝk (f, τ ) is the speech component after channel-wise SS at thedure can be expressed bykth channel, and ∆s (f, τ ) is the speech component distortedk[∗]y(f, τ ) = wT
DS (f )x (f, τ ),
(DS)

(20)

(DS)

wDS (f ) = [w1 (f ), . . . , wK (f )],
(21)
1
(DS)
wk (f ) =exp(−2j(f /N )fs dk sin θU /c),
K
(22)
)
([
]

by channel-wise SS. Also, we assume that s(f, τ ), ∆s1 (f, τ ),and ∆s2 (f, τ ) are uncorrelated with each other. Obviously,ŝ1 (f, τ ) and ŝ2 (f, τ ) can be regarded as being generated by apoint source if ∆s1 (f, τ ) and ∆s2 (f, τ ) are zero, i.e., a valid
static blocking filter can be obtained by ICA as, (23)[W ICA (f )]11 ŝ1 (f, τ ) + [W ICA (f )]12 ŝ2 (f, τ )= ([W ICA (f )]11 h1 (f ) + [W ICA (f )]12 h2 (f ))s(f, τ )
= 0,
(26)

where y(f, τ ) is the final output signal of iterative
BSSA, wDS is the filter coefficient vector of DS, N
is the DFT size, fs is the sampling frequency, dk iswhere we assume U = 1 and, e.g., [W ICA (f )]11 = h2 (f )and [W ICA (f )]12 = −h1 (f ). However, if ∆s1 (f, τ ) and∆s2 (f, τ ) become nonzero as a result of SS, ICA does not

W ICA (f )−1] kU[∗]
W ICA (f )−1[∗]

arg
θU = sin−1

[

k′ U

2πfs c−1 (dk − dk′ )1st iteration
F
F
T

User’sspeech2nd iteration

Channelwise
spectralsubtraction

Final iteration

Channelwise
spectralsubtraction

Channelwise
spectralsubtraction

θU

0 Projection
back

θU

Noise
FDICA

FDICA

0 Projection
back

FDICA

0 Projection
back

Σ
Phase
compensation

Fig. 4. Block diagram of proposed iterative BSSA.

have a valid speech blocking filter with a static (time-invariant)form.
Second, the cosine distance between speech power spectra|ŝ1 (f, τ )|2 and |ŝ2 (f, τ )|2 is introduced in each frequency
subband to indicate the degree of point-source-likeness, as∑|ŝ1 (f, τ )|2 |ŝ2 (f, τ )|2√∑
COS(f ) = √∑ τ
.
(27)
4
4
τ |ŝ1 (f, τ )|
τ |ŝ2 (f, τ )|From (27), the cosine distance reaches its maximum value ofunity if and only if ∆s1 (f, τ ) = ∆s2 (f, τ ) = 0, regardless ofthe values of h1 (f ) and h2 (f ), meaning that the SS-applied
speech signals ŝ1 (f, τ ) and ŝ2 (f, τ ) can be assumed to be
produced by the point source. The value of COS(f ) decreases
with increasing magnitudes of ∆s1 (f, τ ) and ∆s2 (f, τ ) as
well as the difference between h1 (f ) and h2 (f ); this indicatesthe non-point-source state.
Third, we evaluate the degree of point-source-likeness ineach iteration of iterative BSSA by using COS(f ). We statistically estimate the distorted speech component of the enhanced
signal in each iteration. Here, we assume that the original
speech power spectrum |s(f, τ )|2 obeys a gamma distribution
with a shape parameter of 0.1 (this is a typical value forspeech) as|s(f, τ )|2 ∼x−0.9
exp(−x/θs ),Γ(0.1)θs0.1[i+1]|ŝk (f, τ )|2 =

[i]
2
i
1−i

 |ŝk (f, τ )| − βθn M (αn , β, η, 1)αn[i]
(if |ŝk (f, τ )|2 > βθn Mi (αn , β, η, 1)αn1−i ), 2 [i]η |ŝk (f, τ )|2 (otherwise).

(29)

and thus, the amount of residual noise after the ith iteration
is given byHere, we define the input SNR as the average of both channel
SNRs,
(
)
0.1|h2 (f )|2 θs
1 0.1|h1 (f )|2 θs+
ISNR(f ) =
2
αn θn
αn θn
0.1θs
(|h1 (f )|2 + |h2 (f )|2 ).
(32)=
2αn θnIf we normalize the speech scale parameter θs to unity, from(32), the noise scale parameter θn is given byθn =iµ1 = θn αn {M(αn , β, η, 1)/αn }= θn Mi (αn , β, η, 1)αn1−i .

0.1(|h1 (f )|2 + |h2 (f )|2 ),2αn ISNR(f )(33)

and using (33), we can reformulate (31) as|ŝk (f, τ )|2 =[i]0.1(|h1 (f )|2+|h2 (f )|2 )2Mi (αn , β, η, 1)αn−i |ŝk (f, τ )| −β
2ISNR(f )
[i]

(30)

Next, we assume that speech and noise are disjoint, i.e.,there are no overlaps in the time-frequency domain, and thatspeech distortion is caused by subtracting the average noise

0.1(|h1 (f )|2+|h2 (f )|2 )(if |ŝk (f, τ )| > β
2ISNR(f ) η 2 |ŝ[i] (f, τ )|2 (otherwise).
2

Mi (αn , β, η, 1)αn−i ),k(34)Furthermore, we define the transfer function ratio (TFR) as
TFR(f ) = |h1 (f )/h2 (f )|2 ,(35)

and if we normalize |h1 (f )|2 to unity in each frequency
subband, |h1 (f )|2 +|h2 (f )|2 becomes 1+1/TFR(f ). Finally,
we express (34) in terms of the input SNR ISNR(f ) and thetransfer function ratio TFR(f ) as
[i+1]

[i]

(31)

[i+1]

(28)

where θs is the speech scale parameter. Regarding the amountof noise to be subtracted, the 1st-order moment of the noise
power spectra is equal to θn αn when the number of iterations,i, equals zero. Also, the value of αn does not change in each
iteration when we use the specific parameters β and η thatsatisfy the musical-noise-free condition because the kurtosis
ratio does not change in each iteration. If we perform SS only
once, the rate of noise decrease is given byM(αn , β, η, 1)/αn ,from the pure speech component. Thus, the speech component[i+1]|ŝk (f, τ )|2 at the kth channel after the ith iteration isrepresented by subtracting the amount of residual noise (30)as|ŝk (f, τ )|2 = [i]0.1(1+1/TFR(f ))2i−i |ŝk (f, τ )| −β 2ISNR(f ) M (αn , β, η, 1)αn[i]
0.1(1+1/TFR(f))(if |ŝk (f, τ )|2 > β 2ISNR(f ) Mi (αn , β, η, 1)αn−i ), 2 [i]η |ŝk (f, τ )|2 (otherwise).
(36)

2

Value of TFR(f) (|h1(f)/h2(f)| )1.0(a) Input SNR is set to 10 dB.

0.9

0.8

0.7

0.6

0.5

(b) Input SNR is set to 5 dB.

10.4(c) Input SNR is set to 0 dB.

1

1

0.99

0.98

0.97

3

i=1

i=2

4

6

5

i=3

i=4

i=5

i=6

Cosine distance

Cosine distance

Cosine distance

0.9
0.95

0.9

i=7

7 8 9 10 11 12 13 14 15
Noise reduction rate [dB]

0.85

3

i=1

i=2

4

6

5

i=3

i=4

i=5

i=6

i=7

7 8 9 10 11 12 13 14 15
Noise reduction rate [dB]

0.8
0.7
0.6
0.5
0.4

3

i=1

i=2

4

6

5

i=3

i=4

i=5

i=6

i=7

7
8 9 10 11 12 13 14 15
Noise reduction rate [dB]

As can be seen, the speech component is subjected to greater
subtraction and distortion as ISNR(f ) and/or TFR(f ) decrease.
Figure 5 shows the relation between the TFR and thecorresponding value of COS(f ) calculated by (27) and (36).
In Fig. 5, we plot the average of COS(f ) over the whole
frequency subbands. The noise shape parameter αn is set to 0.2
with the assumption of super-Gaussian noise (this corresponds
to the real noises used in Sect. V), the input SNR is set
to 10 dB, 5 dB, or 0 dB, and the noise scale parameter θn
is uniquely determined by (33) and the previous parametersettings. The TFR is set from 0.4 to 1.0 (|h1 (f )| is fixed
to 1.0). Note that the TFR is highly correlated to the room
reverberation and the interelement spacing of the microphonearray; we determined the range of the TFR by simulatinga typical moderately reverberant room and the array with
2.15 cm interelement spacing used in Sect. V (see the exampleof the TFR in Fig. 6). For the internal parameters used in
iterative BSSA in this simulation, β and η are 8.5 and 0.9,
respectively, which satisfy the musical-noise-free condition.
In addition, the smallest value on the horizontal axis is 3 dB
in Fig. 5 because DS is still performed even when i = 0.
From Figs. 5(a) and (b), which correspond to relatively highinput SNRs, we can confirm that the degree of point-sourcelikeness, i.e., COS(f ), is almost maintained when the TFR
is close to 1 even if the speech components are distorted by
iterative BSSA. Also, it is worth mentioning that the degreeof point-source-likeness is still above 0.9 even when the TFR
is decreased to 0.4 and i is increased to 6. This means that
almost 90% of the speech components can be regarded as apoint source and thus can be blocked by ICA. In contrast, from
Fig. 5(c), which shows the case of a low input SNR, when the
TFR is dropped to 0.4 and i is more than 3, the degree ofpoint-source-likeness is lower than 0.6. Thus, less than 60%of the speech components can be regarded as a point source.
However, this is a worst-case scenario; actually when the TFR
is dropped to 0.4 and i is more than 3, the degree of pointsource-likeness is lower than 0.6. Thus, more than 40% of the
speech components cannot be regarded as a point source, andthis leads to poor noise estimation.

TFR(f) (|h1(f)/h2(f)|)

Fig. 5. Relation between number of iterations of iterative BSSA and cosine distance. Input SNR is (a) 10 dB, (b) 5 dB, and (c) 0 dB.

1
0.8
0.6
0.4
0.2
0

0

1000 2000 3000 4000 5000 6000 7000 8000
Frequency [Hz]Fig. 6. Typical examples of TFR(f ) (|h1 (f )/h2 (f )|2 ) in each frequency
subband, where solid and broken lines are different combinations of microphones.

IV. I MPROVEMENT SCHEME FOR POOR NOISE ESTIMATION
A. Channel selection in ICA
In this subsection, we propose a channel selection strategy
in ICA for achieving high accuracy of noise estimation.
As mentioned previously, speech distortion is subjected to
ISNR(f ) and TFR(f ), and the accuracy of noise estimation
is degraded along with its speech distortion. Figure 6 shows atypical example of the TFR. From Fig. 6, we can confirm that
the TFRs in different combinations of microphones are not
the same in each frequency subband; in a specific frequency,
one microphone pair has higher TFR(f ) than another pair, andvice versa in another frequency. Thus, we are able to select theappropriate combination of the microphones to obtain higher
TFR.
Therefore, we introduce the channel selection method into
ICA in each frequency subband, where we automaticallychoose less varied inputs to maintain high accuracy of noise
estimation. Hereafter, we describe the detail of the channel
selection method. First, we calculate the average power of theobserved signal xk (f, τ ) at the kth channel asEτ [|xk (f, τ )|2 ] = Eτ [|s(f, τ )|2 ]|hk (f )|2 +Eτ [|nk (f, τ )|2 ].
(37)
Here, Eτ [|s(f, τ )|2 ] is a constant, and if we assume the diffuse
noise field, Eτ [|nk (f, τ )|2 ] is also a constant. Thus, we can

1.2

1.2

(a) 4-microphone case

0.8
0.6
0.4
0.2
0.0
-0.2

(b) 3-microphone case

1.0
Log kurtosis ratio

Log kurtosis ratio

1.0

0.8
0.6
0.4
0.2
0.0

-5 dB

0 dB

1.2

5 dB

-0.2

10 dB

(c) 2-microphone case

-5 dB

0 dB

5 dB

10 dB

MMSE-STSA estimator (commonly used speech enhancement method)

Log kurtosis ratio

1.0
Iterative spectral subtraction (single channel musical-noise-free speech enhancement method)
0.8
BSSA (conventional blind speech extraction method)
0.6
Iterative BSSA using ICA
0.4
Iterative BSSA with channel selection using ICA
0.2
Iterative BSSA using time-variant nonlinear noise estimator
0.0
Iterative BSSA with channel selection using time-variant nonlinear noise estimator
-0.2

-5 dB

0 dB

5 dB

10 dB

Fig. 7. Kurtosis ratio obtained from experiment for traffic noise under 10-dB NRR condition.

5
4
3
2
1
0

-5 dB

0 dB

8

Cepstral distortion [dB]6(a) 4-microphone case
Cepstral distortion [dB]

Cepstral distortion [dB]

6

5 dB

4
3
2
1
0

10 dB

(b) 3-microphone case

5

-5 dB

0 dB

5 dB

10 dB

(c) 2-microphone case

7

MMSE-STSA estimator (commonly used speech enhancement method)

6

Iterative spectral subtraction (single channel musical-noise-free speech enhancement method)

5

BSSA (conventional blind speech extraction method)

4

Iterative BSSA using ICA

3

Iterative BSSA with channel selection using ICA

2Iterative BSSA using time-variant nonlinear noise estimator

10

Iterative BSSA with channel selection using time-variant nonlinear noise estimator
-5 dB

0 dB

5 dB

10 dB

-5 dB

0 dB

5 dB

10 dB

Fig. 8. Cepstral distortion obtained from experiment for traffic noise under 10-dB NRR condition.

estimate the relative order of |hk (f )|2 by comparing (37) for
every k.
Next, we sort Eτ [|xk (f, τ )|2 ] in descending order and selectthe channels k corresponding to high amplitude of |hk (f )|2
satisfying the following condition:
max Eτ [|xk (f, τ )|2 ] · ξ ≤ Eτ [|xk (f, τ )|2 ],k(38)

where ξ(< 1) is the threshold for the selection.
Finally, we perform noise estimation based on ICA using
the selected channels in each frequency subband, and we apply
the projection back operation to remove the ambiguity of theamplitude and construct the estimated noise signal.

B. Time-variant noise PSD estimator
In the previous section, we reveal that the speech components cannot be regarded as a point source, and this leads topoor noise estimation in iterative BSSA. To solve this problem,
we introduce a time-variant noise PSD estimator [15] insteadof ICA for improvement of the noise estimation accuracy.
This method has been developed for future high-end binauralhearing aids and performs a prediction on the left noisy signal
from the right noisy signal via an Wiener filter, followed byan auto-PSD of the difference between the left noisy signal
and the prediction. By applying the noise PSD estimated from

1.2

1.2

(a) 4-microphone case

0.8
0.6
0.4
0.2
0.0
-0.2

(b) 3-microphone case

1.0
Log kurtosis ratio

Log kurtosis ratio

1.0

0.8
0.6
0.4
0.2
0.0

-5 dB

0 dB

1.2

5 dB

-0.2

10 dB

(c) 2-microphone case

-5 dB

0 dB

5 dB

10 dB

MMSE-STSA estimator (commonly used speech enhancement method)

Log kurtosis ratio

1.0
Iterative spectral subtraction (single channel musical-noise-free speech enhancement method)
0.8
BSSA (conventional blind speech extraction method)
0.6
Iterative BSSA using ICA
0.4
Iterative BSSA with channel selection using ICA
0.2
Iterative BSSA using time-variant nonlinear noise estimator
0.0
Iterative BSSA with channel selection using time-variant nonlinear noise estimator
-0.2

-5 dB

0 dB

5 dB

10 dB

Fig. 9. Kurtosis ratio obtained from experiment for railway station noise under 10-dB NRR condition.

5
4
3
2
1
0

-5 dB

0 dB

6
Cepstral distortion [dB]6(a) 4-microphone case
Cepstral distortion [dB]

Cepstral distortion [dB]65 dB(c) 2-microphone case

4
3
2
1
0

10 dB

(b) 3-microphone case

5

-5 dB

0 dB

5 dB

10 dB

MMSE-STSA estimator (commonly used speech enhancement method)5
Iterative spectral subtraction (single channel musical-noise-free speech enhancement method)
4
BSSA (conventional blind speech extraction method)3
Iterative BSSA using ICA
2Iterative BSSA with channel selection using ICA
1Iterative BSSA using time-variant nonlinear noise estimator
0

-5 dB

0 dB

5 dB

10 dB

Iterative BSSA with channel selection using time-variant nonlinear noise estimator

Fig. 10. Cepstral distortion obtained from experiment for railway station noise under 10-dB NRR condition.

this estimator to (19), we can perform the speech extraction.
The procedure of this noise PSD estimator is described inAppendix.
V. E XPERIMENT IN R EAL W ORLD

were sampled at 16 kHz with 16-bit accuracy. The FFT size
was 1024, and the frame shift length was 256. We use 10
speakers (5 males and 5 females) as sources of the original
target speech signal. The input SNR was -5, 0, 5, and 10 dB.

A. Experimental conditions

B. Objective evaluation

We conducted objective and subjective evaluation experiments to confirm the validity of the proposed methods under
the diffuse and nonstationary noise condition. The size of theexperimental room was 4.2×3.5×3.0 m3 and the reverberation
time was approximately 200 ms. We used a two-, three- or
four-element microphone array with an interelement spacingof 2.15 cm, and the direction of the target speech was set to
be normal to the array. All the signals used in this experiment

We conducted an objective experiment evaluation under thesame NRR condition. Figures 7, 8, 9, and 10 show the kurtosis
ratio and cepstral distortion obtaind from the experimentswith real traffic noise and railway station noise, where weevaluate 10-dB NRR (i.e., output SNRs = 5, 10, 15, and20 dB) signals processed by three conventional methods,
namely, the minimum mean-square error (MMSE) short-time
spectral amplitude (STSA) estimator [21], simple combination

100

100

80

80

Preference score [%]Preference score [%]Conventional BSSA
95% confidence interval
Iterative BSSA using ICA
Iterative BSSA with channel selection using ICA
Iterative BSSA using time-variant nonlinear noise estimator
Iterative BSSA with channel selection using time-variant nonlinear noise estimator

60
40
20
0

60
40
20
0

(a) Traffic noise case(b) Railway station noise caseFig. 11. Subjective evaluation results for (a) traffic noise and (b) railway station noise.

of BSSA and musical-noise-free iterative spectral subtraction,and our proposed methods, iterative BSSA (using ICA ortime-variant noise estimator with/without channel selection),
respectively. Here, we did not apply the channel selection
method to the two-microphone case because ICA or timevariant noise estimation need at least two-channel signals.
Also, we applied minimum statistics noise PSD estimator [5]to the MMSE STSA estimator and musical-noise-free iterative
spectral subtraction, and we use the decision-directed approach
for a priori SNR estimation in the MMSE STSA estimator.
From Figs. 7 and 9, we can confirm that iterative BSSA
methods outperform the MMSE STSA estimator and theconventional BSSA in the kurtosis ratio. In particular, thekurtosis ratios of the proposed methods are mostly close
to 1.0. This means that the proposed iterative methods did
not generate any musical noise. However, iterative BSSA
methods leads to greater speech distortion compared with theconventional BSSA (see Figs. 8 and 10). Therefore, a tradeoff exists between the amount of musical noise generationand speech distortion in the conventional BSSA and iterative
BSSA methods.
C. Subjective evaluation
Since we found the above-mentioned trade-off, we next
conducted a subjective evaluation for setting the performancecompetition. In the evaluation, we presented a pair of 10-dB
NRR signals processed by the conventional BSSA and four
our proposed iterative BSSAs (using ICA or time-variant noise
estimator with/without channel selection) in random order to
10 examinees, who selected which signal they preferred fromthe viewpoint of total sound quality, e.g., less musical noise,less speech distortion, etc.
The result of experiment is shown in Fig. 11 for (a) traffic
noise and (b) railway station noise. It is found that the output

signals of some iterative BSSAs are preferred to that of theconventional BSSA, indicating the higher sound quality of the
proposed method in terms of human perception. This result isplausible because humans are often more sensitive to musicalnoise than to speech distortion as indicated in past studies,e.g., [22].
VI. C ONCLUSION
In this paper, we first proposed iterative BSSA using a new
BSSA structure, which generates almost no musical noise evenwith increasing noise reduction performance. Our theoreticalanalysis indicates that the accuracy of noise estimation isdegraded along with its speech distortion. Next, we conducteda channel selection strategy in ICA for achieving high accuracyof noise estimation. From the evaluation experiments, it was
shown that there is a trade-off between the amount of musicalnoise generation and speech distortion in both the conventional
BSSA and iterative BSSA. However, in a subjective preference
test, iterative BSSA obtained a higher preference score thanthe conventional BSSA. Thus, iterative BSSA is advantageous
to the conventional BSSA in terms of sound quality.
A PPENDIX
This appendix provides a brief review of the time-variant
nonlinear noise estimator. For more detailed information,
Ref. [15] can be available.
Let xL (f, τ ) and xR (f, τ ) be noisy signals received at
the left and right microphones in the time-frequency domain,
defined as
xL (f, τ ) = hL (f )s(f, τ ) + nL (f, τ ),xR (f, τ ) = hR (f )s(f, τ ) + nR (f, τ ),

(39)
(40)

where hL (f ) and hR (f ) are the left and right transfer functions, respectively. Next, the left and right auto-power spectraldensities, ΓLL (f ) and ΓRR (f ), can be expressed as follows:

R EFERENCES[1] P. C. Loizou, Speech Enhancement Theory and Practice CRC Press,
Taylor & Francis Group FL, 2007.
[2] S. F. Boll, “Suppression of acoustic noise in speech using spectral subΓRR (f, τ ) = |HR (f )|2 ΓSS (f, τ ) + ΓNN (f, τ ),(42)
traction,” IEEE Transactions on Acoustics, Speech, Signal Processing,vol.27, no.2, pp.113–120, 1979.
where ΓSS (f, τ ) is the power spectral density of the target[3] M. Berouti, R. Schwartz, and J. Makhoul, “Enhancement of speechcorrupted by acoustic noise,” Proc. ICASSP79, pp.208–211, 1979.
speech signal, and ΓNN (f, τ ) is the power spectral density of[4] R. McAulay and M. Malpass, “Speech enhancement using a softthe noise signal. In this paper, we assume that the left and rightdecision noise suppression filter,” IEEE Transactions on Acoustics,
noise power spectral densities are approximately the same, i.e.,
Speech and Signal Processing, vol.28, no.2, pp.137–145, 1980.
[5] R. Martin, “Spectral subtraction based on minimum statistics,” Proc.
ΓNL NL (f, τ ) ≃ ΓNR NR (f, τ ) ≃ ΓNN (f, τ ).
EUSIPCO94, pp.1182–1185, 1994.
Next, we consider the Wiener solution between the left and[6] K. Yamashita, S. Ogata, T. Shimamura, “Spectral subtraction iterated
right transfer functions, which is defined aswith weighting factors,” Proc. IEEE Speech Coding Workshop, pp.138–
140, 2002.
[7] M. R. Khan, T. Hansen, “Iterative noise power subtraction technique forΓLR (f, τ )HW (f, τ ) =,(43)
improved speech quality,” Proc. ICECE2008, pp.391–394, 2008.
ΓRR (f, τ )[8] S. Li, J.-Q. Wang, M. Niu, X.-J. Jing, T. Liu, “Iterative spectral subtraction method for millimeter-wave conducted speech enhancement,”
where ΓLR (f ) is the cross-power spectral density between
Journal of Biomedical Science and Engineering, vol.2010, no.3, pp.187–
192, 2010.
the left and the right noisy signals. The cross-power spectral[9] T. Inoue, H. Saruwatari, Y. Takahashi, K. Shikano, K. Kondo, “Theodensity expression then becomesretical analysis of iterative weak spectral subtraction via higher-order
statistics,” Proc. MLSP2010, pp.220–225, 2010.
∗
ΓLR (f, τ ) = ΓSS (f, τ )HL (f )HR (f ).
(44) [10] R. Miyazaki, H Saruwatari, T Inoue, Y. Takahashi, K. Shikano,
K. Kondo, “Musical-noise-free speech enhancement based on optimized
iterative spectral subtraction,” IEEE Transactions on Audio, Speech and
Therefore, substituting (44) into (43) yields
Language Processing, vol.20, no.7, pp.2080–2094, 2012.
[11] R. Miyazaki, H. Saruwatari, K. Shikano, K. Kondo, “Musical-noise-freeΓSS (f, τ )HL (f )HR∗ (f )blind speech extraction using ICA-based noise estimation and iterative
HW (f, τ ) =.
(45)
ΓRR (f, τ )spectral subtraction,” Proc. ISSPA2012, pp.322–327, 2012.
[12] R. Miyazaki, H. Saruwatari, K. Shikano, K. Kondo, “Musical-noise-freeblind speech extraction using ICA-based noise estimation and iterative
Furthermore, using (41) and (42), the squared magnitude
spectral subtraction with channel selection,” Proc. IWAENC2012, 2012.
response of the Wiener solution in (45) can be also expressed[13] P. Comon, “Independent component analysis, a new concept?,” Signal
asProcessing, vol.36, pp.287–314, 1994.
[14] H. Saruwatari, T. Kawamura, T. Nishikawa, A. Lee, K. Shikano, “Blind(ΓLL (f, τ )−ΓNN (f, τ ))(ΓRR (f, τ )−ΓNN (f, τ ))source separation based on a fast-convergence algorithm combining ICA|HW (f, τ )|2 =.
and beamforming,” IEEE Transactions on Speech and Audio Processing,Γ2RR (f, τ )vol.14, no.2, pp.666–678, 2006.
(46)[15] A. Homayoun, M. bouchard, “Improved noise power spectrum density
estimation for binaural hearing aids operating in a diffuse noise field
Equation (46) is rearranged into a quadratic equation as in theenvironment,” IEEE Transactions on Audio, Speech, and Language
Processing, vol.17, no.4, pp.521–533, 2009.
following:[16] Y. Uemura, Y. Takahashi, H. Saruwatari, K. Shikano, K. Kondo, “Au2tomatic optimization scheme of spectral subtraction based on musicalΓNN (f, τ ) − ΓNN (f, τ ) (ΓLL (f, τ ) + ΓRR (f, τ ))noise assessment via higher-order statistics,” Proc. IWAENC2008, 2008.
+ ΓEE (f, τ )ΓRR (f, τ ) = 0,(47) [17] Y. Takahashi, T. Takatani, K. Osako, H. Saruwatari, K. Shikano, “Blindspatial subtraction array for speech enhancement in noisy environment,”
IEEE Transactions on Audio, Speech, and Language Processing, vol.17,
whereno.4, pp.650–664, 2009.
[18] Y. Takahashi, H. Saruwatari, K. Shikano, K. Kondo, “Musical-noise2ΓEE (f, τ ) = ΓLL (f, τ ) − ΓRR (f, τ )|HW (f )| .
(48)
analysis in methods of integrating microphone array and spectral subtraction based on higher-order statistics,” EURASIP Jouranl on Advancesin Signal Processing, vol.2010, Article ID 431347, 25 pages, 2010.
Consequently, the noise power spectral density ΓNN (f ) can be[19] H. Saruwatari, Y. Ishikawa, Y. Takahashi, T. Inoue, K. Shikano,estimated by solving the quadratic equation in (47) as follows:
K. Kondo, “Musical noise controllable algorithm of channelwise spectral
subtraction and adaptive beamforming based on higher-order statistics,”
1
IEEE Transactions on Audio, Speech and Language Processing, vol.19,
ΓNN (f, τ ) = (ΓLL (f, τ ) + ΓRR (f, τ )) − ΓLRavg (f, τ ),no.6, pp.1457–1466, 2011.
2(49) [20] R. Miyazaki, H. Saruwatari, K. Shikano, “Theoretical analysis ofamounts of musical noise and speech distortion in structure-generalized
1
parametric spatial subtraction array,” IEICE Transactions Fundamentals,
ΓLRavg (f, τ ) = {(ΓLL (f, τ ) + ΓRR (f, τ ))2
vol.95-A, no.2, pp.586–590, 2012.
2− 4ΓEE (f, τ )ΓRR (f, τ )}0.5 .
(50) [21] Y. Ephraim, and D. Malah, “Speech enhancement using a minimum
mean-square error short-time spectral amplitude estimator,” IEEE Transactions on Acoustics, Speech and Signal Processing, vol.32, no.6,
ACKNOWLEDGMENT
pp.1109–1121, 1984.
[22] Y. Uemura, Y. Takahashi, H. Saruwatari, K. Shikano, K. Kondo,“Musical noise generation analysis for noise reduction methods based
This work was supported by JST Core Research of Evalon spectral subtraction and MMSE STSA estimation,” Proc. ICASSP,
uational Science and Technology (CREST) and Grant-in-Aid
pp.4433–4436, 2009.

ΓLL (f, τ ) = |HL (f )|2 ΓSS (f, τ ) + ΓNN (f, τ ),for JSPS Fellows.

(41)