VOL. E104-D NO. 3
MARCH 2021
The usage of this PDF file must comply with the IEICE Provisions
on Copyright.
The author(s) can distribute this PDF file for research andeducational (nonprofit) purposes only.
Distribution by anyone other than the author(s) is prohibited.

IEICE TRANS. INF. & SYST., VOL.E104–D, NO.3 MARCH 2021

441

PAPER

Noise Robust Acoustic Anomaly Detection System with Nonnegative
Matrix Factorization Based on Generalized Gaussian Distribution∗
Akihito AIBA†a) , Minoru YOSHIDA† , Nonmembers, Daichi KITAMURA††,††† , Member,
Shinnosuke TAKAMICHI††† , Nonmember, and Hiroshi SARUWATARI††† , Member

SUMMARY We studied an acoustic anomaly detection system forequipments, where the outlier detection method based on recorded sounds
is used. In a real environment, the SNR of the target sound against background noise is low, and there is the problem that it is necessary to catchslight changes in sound buried in noise. In this paper, we propose a system in which a sound source extraction process is provided at the preliminary stage of the outlier detection process. In the proposed system, nonnegative matrix factorization based on generalized Gaussian distribution(GGD-NMF) is used as a sound source extraction process. We evaluatedthe improvement of the anomaly detection performance in a low-SNR environment. In this experiment, SNR capable of detecting an anomaly was
greatly improved by providing GGD-NMF for preprocessing.
key words: nonnegative matrix factorization, generalized Gaussian distribution, anomaly detection, outlier detection

1.

Introduction

Several systems have been proposed with which abnormal
operation is detected from sensor signals by constantly sensing the state of facilities such as water supply and drainage
and air conditioning in place of human beings [1], [2]. Here,we consider an anomaly detection system based on the operating sound recorded by a microphone.
One of the major problems in anomaly detection applications is the diﬃculty in obtaining samples at anomalous times. In particular, it is diﬃcult to intentionally createan anomalous state of a facility that is expensive and important for building management, for the purpose of sample
collection. Also, waiting for a naturally occurring anomaly
is impractical as some anomalies occur only once every few
months or years. For this reason, outlier detection is often adopted in anomaly detection applications [3], [4]. In
outlier detection, only easily available normal samples arecollected in advance. Then their features are modeled, andsamples that deviate from the model are judged to be anomalies when obtained. In this research, we also use this outlier
detection approach.
Manuscript received May 18, 2020.
Manuscript revised October 14, 2020.
Manuscript publicized December 18, 2020.
†
The authors are with Ricoh Company Ltd., Ebina-shi, 243–
0460 Japan.
††
The author is with the National Institute of Technology,
Kagawa College, Takamatsu-shi, 761–8058 Japan.
†††
The authors are with the University of Tokyo, Tokyo, 113–8656 Japan.
∗
This is a paper on system development.
a) E-mail: akihito.aiba@jp.ricoh.com
DOI: 10.1587/transinf.2020EDK0002

When sound data is used as input, the background noise
in a real environment also becomes a problem. We often
see a situation where multiple equipment units are installed
in one room and generate very loud operation noise. In
such a situation, even if the operation sound of one deviceshows an abnormal feature, its detection is very diﬃcult because the signal power of the anomalous sound is lower
than that of the operation sound of other devices. This is
a common problem not only in outlier detection but also insound discrimination. There are existing methods for improving the SNR by extracting sound sources before discrimination [5], [6].
In this paper, a method based on nonnegative matrix
factorization (NMF) [7] is used for sound source extraction
processing. Regarding the sound source extraction, deepneural network (DNN)-based processing [8]–[10] has been
proposed for general use in recent years. However, it cannot be introduced in our anomaly detection application because DNN requires an enormous number of clean objective sounds but we cannot obtain them in advance. On theother hand, the proposed NMF-based method has an advantage to requiring no clean data. In particular, we extend the NMF based on generalized Gaussian distribution(GGD-NMF) [11] proposed by Kitamura to a semisupervised NMF (SS-NMF) [12], [13] and confirm its eﬃcacy inthe anomaly detection system. GGD-NMF incorporates a
generative model based on GGD into NMF to control the
low-rank property during decomposition. The evaluation is
performed using data recorded in a real environment.
This paper is organized as follows. In Sect. 2, we introduce the basis theory of NMF and its application to audio
signals. In Sect. 3, the details of the proposed anomaly detection system are explained. The experimental evaluationof the proposed system is described in Sect. 4. The paper isconcluded in Sect. 5.
2. Conventional NMF
In this section, the NMF algorithm is described first. Next,
Itakura–Saito NMF (ISNMF) [14], which has a Gaussiandistribution-based generative model, is introduced.
2.1

NMF

NMF is a mathematical algorithm for extracting a limited
number of nonnegative features from a nonnegative matrix.

c 2021 The Institute of Electronics, Information and Communication Engineers
Copyright ⃝IEICE TRANS. INF. & SYST., VOL.E104–D, NO.3 MARCH 2021

442
I×J
be a nonnegative observation matrix to be
Let X ∈ R≥0analyzed. In the field of acoustics, the complex spectrogram C ∈ CI×J is obtained by applying short-time Fouriertransform (STFT) to the time-domain signal. Then, the nonnegative spectrogram X is calculated from C, e.g., the amplitude spectrogram (X = |C|.1 ) or the power spectrogram(X = |C|.2 ), where | · |.p returns a matrix with the elementwise absolute value and the pth power of the input matrix.
I and J are the numbers of rows and columns in X, respectively. In the case of a spectrogram, I and J correspond to
the numbers of frequency bins and time frames, respectively.
NMF is an approximate decomposition as follows:

X ≈ X̂ = FG∑=f k gTk ,

(1)
(2)

k
I×K
is the basis matrix, andwhere F = ( f 1 . . . f K ) ∈ R≥0
K×J
T
G = (g1 . . . gK ) ∈ R≥0 is the coeﬃcient matrix or activation matrix. k = 1, 2, . . . , K is the index of bases. K is thenumber of bases and is generally set to a value suﬃciently
smaller than I and J. Therefore, Eq. (1) is interpreted as a
low-rank approximation that represents the observation matrix X with a limited number of bases. The basis vectorsI
f k ∈ R≥0
represent the frequently appearing spectral patterns (parts) in the spectrogram X, and the activation vectors
J
represent their time-varying gains [11].
gk ∈ R≥0
The matrix variables F and G can be estimated by solving the following optimization problem:

min D(X|FG) s.t. fik , gk j ≥ 0 ∀i, j, k,
F,G

(3)

where i = 1, 2, . . . , I is the index of frequency bins andj = 1, 2, . . . , J is the index of time frames. fik and gk j are
the elements of F and G, respectively. D(·|·) is a similarity
function between two input matrices. For example, squared
Euclidean distance [15], generalized Kullbak–Leibler (KL)divergence [15], and Itakura–Saito (IS) divergence [14] areoften used.
2.2 ISNMF
In this section, we consider NMF based on IS divergence forthe similarity function D in Eq. (3) (ISNMF). IS divergenceis a similarity function defined as
|c|2
|c|2
−
log
− 1.
(4)
σ2
σ2
∑
By setting |c|2 = xi j and σ2 = k fik gk j , we can rewrite theoptimization problem of ISNMF as∑  xi j
∑



+ log
min
fik gk j  ∑F,Gk fik gk ji, jkDIS (c|σ) =s.t. fik , gk j ≥ 0 ∀i, j, k,(5)where xi j is the element of the complex spectrogram X andthe constant terms are omitted.

In ISNMF, a generative model of the complex spectrogram C is assumed as described below. Let us assume thateach time-frequency element of the complex spectrogram C,
ci j , can∑be decomposed into K spectral components ci j,k as
ci j = k ci j,k , and ci j,k is assumed to be generated from azero-mean circularly symmetric complex Gaussian distribution with the variance σ2i j,k > 0 [14]. Then, the generativemodel of the complex spectrogram C becomes∏∏ ∑ 2 σi j,k  .
p(ci j ) =NC 0,
C∼
(6)
i, j

i, j

k

On the basis of this generative model, we consider
the maximum likelihood estimation problem of the varianceσ2i j,k . The negative log-likelihood function of C is∏ ∑ 2 σi j,k − log L = − log
NC 0,
i, j

k



∑  |ci j |2
∑


2

=
σi j,k + log π .
 ∑ 2 + log
σ
k

i, j

i j,k

(7)

k

The maximum likelihood estimate of the variance σ2i j,k can
be obtained by minimizing Eq. (7). Here, when comparing
Eq. (7) with the ISNMF minimization problem (5), if xi j =|ci j |2 and σ2i j,k = fik gk j , both equations become equivalent up
to the constant term.
From the above, it is seen that applying ISNMF to thepower spectrogram X = |C|.2 is equivalent to the maximum
likelihood estimation of the variance σ2i j,k based on the generative model (6).
3. Proposed System
In this section, GGD-NMF is described first. Then, we describe the proposed system that extends GGD-NMF to SSNMF and incorporates it into the anomaly detection system.
3.1

GGD-NMF

In ISNMF, a complex Gaussian distribution is assumed asa generative model of the observation complex spectrogram
C. GGD-NMF is a generalization of ISNMF, which assumes
the complex GGD as a generative model of C.
By denoting the circularly symmetric complex GGD
with the mean value 0, shape parameter ρ > 0, and scaleparameters σi j that fluctuate along time and frequency as
GC (ci j ; 0, ρ, σi j ), we can describe the generative model in
GGD-NMF as∏C ∼GC (ci j ; 0, ρ, σi j )
i, j

[
(
)ρ ]
2 |ci j |
=
exp −
,
1− 2
2
ρ σi ji, j 2 ρ πσi j Γ (2/ρ)∑
σipj =
fik gk j ,
∏

k

2

ρ1− ρ

(8)
(9)

AIBA et al.: NOISE ROBUST ACOUSTIC ANOMALY DETECTION SYSTEM WITH NMF BASED ON GENERALIZED GAUSSIAN DISTRIBUTION

443

where Γ(·) is the gamma function and p is a domain parameter corresponding to the domain of the low-rank approximation. In the case of p = 1, the amplitude spectrogram |C|.1
is modeled by NMF variables (F and G), and in the case of
p = 2, the power spectrogram |C|.2 is modeled. The complex
GGD GC (ci j ; 0, ρ, σi j ) becomes equivalent to the complex
Gaussian distribution and the complex Laplace distribution
when ρ = 2 and ρ = 1, respectively. Also, the complex
GGD becomes sub- and super-Gaussian when ρ > 2 andρ < 2, respectively.
In the generative model (8), the complex GGD is independently defined in each time-frequency slot, and its scaleparameter σi j can fluctuate along time and frequency axes.
In this case, the macro model that includes all the timefrequency slots, i.e., the generative model of the spectrogram C, always approaches to a super-Gaussian distribution
when σi j dynamically varies along i and j [16]. For thisreason, the generative model (8) with sub-Gaussian GGD(ρ > 2) has a versatility as depicted in Fig. 1 in [16], namely,
the macro model of Eq. (8) with ρ > 2 becomes either subGaussian or super-Gaussian. This versatile model is desirable and often provides improved performance in NMFbased acoustic modeling [16], [17].
Let us clarify the relationship between the generativemodel (8) and the similarity function in NMF. The divergence based on the complex GGD is derived from the difference in log-likelihood (deviance). The divergence can becalculated as follows [11]:
( )ρ
2
2 |c|
DGGD (c||σ) = −2 log |c| − + 2 log σ +
ρ
ρ σ
[( )ρ
( )ρ
]
|c|
2 |c|
(10)
− log
−1 .
=
ρ σ
σ
This divergence (10) can be interpreted as a generalizationof IS divergence (4) for the shape parameter ρ. From theabove, the function to be minimized in GGD-NMF is as follows (except for the constant term):∑∑ 
|ci j |ρ
DGGD (ci j |σi j ) = (∑
) ρp
i, j
i, j
k fik gk j

∑

ρ
(11)
+ log
fik gk j  .
pkThus, GGD-NMF (Eq. (11)) and ISNMF (Eq. (5)) are equivalent algorithms when ρ = p and xi j = |c| p .
The variables fik and gk j that minimize Eq. (11) can becalculated by the auxiliary function technique. The update
rules are derived as follows [11], [18]:∑ pz j ∑ i j 2 gk j  ρ+p( k′ fik′ gk′ j ) ,
fik ← fik  ∑
(12)

1

j ∑k′ fik′ gk′ j gk j∑ pz
 i ∑ i j 2 fik  ρ+p
 ( k′ fik′ gk′ j )
gk j ← gk j  ∑ , i∑ 1fik 
k′ fik′ gk′ j(13)

zi j =(

ρ)
1− ρ pci j p σi j p

.

(14)

In the update rules, zi j is regarded as a virtual observationconsisting of the weighted geometric mean of the observation |ci j | and its low-rank model σi j , and GGD-NMF is interpreted as ISNMF on this virtual observation. When ρ = p,
the virtual observation becomes zi j = |ci j | p , which coincides
with the original ISNMF. When ρ < p, since the virtual
observation zi j is the weighted geometric mean of |ci j | andσi j , GGD-NMF tends to prevent σi j from an overfitting tothe observation |ci j |. This is because a heavy-tailed (superGaussian) distribution is assumed and it allows outliers in|ci j |, resulting in low-rank-enhanced NMF modeling. On theother hand, when ρ > p, the geometric mean corresponds tothe point externally dividing |ci j | and σi j , and the error withthe current approximation is further emphasized, which mitigates the excessive low-rank modeling [11]. Thus, in GGDNMF, it is possible to control the low-rank property of theapproximated model by adjusting the parameters ρ and p.
3.2

GGD-Based SS-NMF

In SS-NMF [12], [13], we have a sample sound (training
dataset) of one source, and its basis matrix is trained in advance by decomposing the sample sound using simple NMF.
Then, the observed mixture spectrogram X is divided into
the components with the pretrained (supervised) basis matrix and the other basis matrix asX ≈ ∑FG + HU ∑=f k gk +hl ul ,k(15)lwhere F is the pretrained basis matrix that has spectral patI×L
isterns of the target source, and H = (h1 . . . hL ) ∈ R≥0
the other basis matrix estimated in the separation stage andrepresents the spectral patterns of the other sources in X. G
L×J
are the coeﬃcient matrices forand U = (u1 . . . uL )T ∈ R≥0
F and H, respectively, and both are also estimated in theseparation stage. L is the number of other bases in H, andl = 1, 2, . . . , L is their index.
The matrices G, H, and U can be estimated by the following optimization problem, which is the separation stageof SS-NMF.
min

H,G,U

D(X|FG + HU)
s.t. gk j , hil , ul j ≥ 0 ∀i, j, k, l

(16)

Since F contains the spectral patterns of the target source
and is fixed in the optimization of Eq. (16), FG and HU
correspond to the estimated nonnegative spectrograms of thetarget source and the other source components, respectively.
Note that the similarity function D(·|·) must be the same asthat used in the training of F.
In the proposed system, we introduce GGD-NMF into
the framework of SS-NMF. The function to be minimizedchanges from Eq. (11) as follows:IEICE TRANS. INF. & SYST., VOL.E104–D, NO.3 MARCH 2021

444

The details of each step are described below.
3.4

Basis Training Step

In the proposed system, first, we obtain the basis matrix
of background noise. A time-domain sound signal x0 is
recorded when the target machine is not operating and only
background noise exists. The spectrogram of x0 , X0 ∈0RΩ×T
is obtained by STFT, where Ω and T 0 are the num≥0
bers of frequency bins and time frames, respectively. It can0be considered as a noise sample spectrogram N0 ∈ RΩ×T
≥0
because X0 does not include the operation sound of the target machine. Then, N0 is decomposed by GGD-NMF as
X0 = N0 ≈ F0W 0 ,Fig. 1∑

DGGD (ci j |σi j ) =i, jSystem overview.

∑ [( |ci j | )ρi, j

σipj =∑kσi j
fik gk j +basis matrix that represents only background noise sepectra
is obtained in this training process.

]+ ρ log σi j ,
∑

(23)

Ω×K
is a basis matrix for background noise andwhere F0 ∈ R≥0
K×T 0
W 0 ∈ R≥0 is a coeﬃcient matrix for F0 . As a result, the

(17)
3.5

hil ul j .

(18)

l

The variables that minimize Eq. (17) can be calculated by
the auxiliary function technique similarly to the derivationof Eqs. (12) and (13) [18]. The update rules for solving
GGD-based SS-NMF are derived as
p
∑
 ρ+p
zi j
 i ∑∑
2 fik 
 ( k′ fik′ gk′ j + l′ hil′ ul′ j )
gk j ← gk j  ∑
(19)
 ,
1∑
 i ∑
f
ik
′
′
′
′
f
g
+
h
u
′
′
( k ik k j l il l j )
p
∑
 ρ+p
zi j
 j ∑∑
2 ul j 


′
′
′
′
f
g
+
h
u
′
′
( k ik k j l il l j ) ,hil ← hil  ∑
(20)

1∑
 j ∑
u
l
j
( k′ fik′ gk′ j + l′ hil′ ul′ j )
p
∑
 ρ+p
zi j
 i ∑∑
2 hil  ( k′ fik′ gk′ j + l′ hil′ ul′ j ) ul j ← ul j  ∑
(21)
 ,
1∑
 i ∑
( k′ fik′ gk′ j + l′ hil′ ul′ j ) hil()ρ
1− ρ pzi j = ci j p σi j p .
(22)
3.3 System Overview
Figure 1 shows the overview of the proposed system using
GGD-NMF for anomaly detection. GGD-NMF is expected
to be eﬀective in improving the anomaly detection performance even in a low-SNR environment. The proposed system consists of the following three major steps.
1. Basis training step: training the basis matrix F0 for
background noise.
2. Detector training step: training to learn the normal operation sound of a target machine, which is extracted
by GGD-based SS-NMF.
3. Detection step: detecting anomalies from the extracted
machine-operation sound.

Detector Training Step

In the detector training step, the operation sound of the target machine is extracted from the noisy mixture sound, andthe extracted sound is used to train the anomaly detector. To
extract the operation sound of the target machine, the basis
matrix F0 obtained in the basis training step is used.
We prepare a noisy sound dataset {x1(1) , x1(2) , . . . , x1(N) }recorded when the target machine was operating normallyin a noisy environment. Their nonnegative spectrograms
(n)
(N)
(2)
{X(1)
1 , X 1 , . . . , X 1 } are calculated by STFT, where X 1 ∈Ω×T (n)R≥0 1 , T 1(n) is the number of time frames, N is the numberof samples in the dataset, and n = 1, 2, . . . , N is their index.
We assume that X(n)1 can be approximated as the sum of twononnegative spectrograms as
(n)
(n)
X(n)
1 ≃ S1 + N 1 ,Ω×T (n)(24)Ω×T (n)11
and N(n)are the nonnegawhere S(n)1 ∈ R≥0
1 ∈ R≥0tive spectrograms of the clean operation sound of the target
machine and the background noise, respectively.
In the detector training step, first, the observed noisy
mixture X(n)1 is decomposed into two estimates, i.e., the(n)background noise Ñ1 and the operation sound of the target(n)machine S̃1 , by GGD-based SS-NMF with the pretrained
basis matrix F0 as follows:

(n)

(n)

X(n)
1 ≃ S̃1 + Ñ 1 ,

(25)

(n)
(n)
S̃1 = H(n)1 U1 ,(n)Ñ1 = F0 G(n)
1 ,

(26)
(27)
L×T (n)(n)Ω×L1
are the basis and cowhere H(n)1 ∈ R≥0 and U 1 ∈ R≥0
eﬃcient matrices that correspond to the operation sound of
K×T (n)1
the target machine as Eq. (26) and G(n)is the co1 ∈ R≥0
eﬃcient matrix of F0 , which corresponds to the backgroundAIBA et al.: NOISE ROBUST ACOUSTIC ANOMALY DETECTION SYSTEM WITH NMF BASED ON GENERALIZED GAUSSIAN DISTRIBUTION

445

noise as Eq. (27).
Next, the anomaly detector is trained using the ex(n)tracted normal operation sound S̃1 as the training data. In
the proposed system, we use an outlier detection methodwith Autoencoder (AE) [4]. AE is a neural network trained
to reconstruct input data through dimensional compression.
For this reason, the reconstruction error of the trained AE
becomes smaller for known inputs, whereas the reconstruction error becomes larger for unknown inputs. Thus, it ispossible to detect outliers on the basis of the reconstruction
error.
A parameter of AE, Θ, is optimized by minimizing aloss function J:[{} 2]− fAE S′(n,m)J(Θ) = E S′(n,m)
,
(28)
Θ
1
1
2 n,m

p

AE
(n,m) . pS′(n,m)
= S̃1
1

,

(29)

AE
∈ RΩ×Tis a spectrogram piece cut out from≥0(n)S̃1 every time frame number T AE , m = 1, 2, . . . , Mn is the(n,m)

Trainingbasis matrix F
Training
detector
Testing
detection

Table 1 Dataset
Items

SNR

Time

Background noise

-

0.4 h

Normal bearing sound+ background noise
Normal bearing sound+ background noise
Abnormal bearing sound+ background noise−24 dB

2.0 h−24 dB

5.2 h−12 dB−18 dB−24 dB

5.2 h
5.2 h
5.2 h

proposed system is evaluated. To confirm the validity ofthe target source extraction process, a system without GGDNMF (omitting basis training step and GGD-based SS-NMF
in the other steps from the proposed system) and the proposed system with various parameter settings were compared.

where S̃1

4.1

index of cut spectrogram pieces, and Mn is the number ofcut spectrogram pieces for the nth sample sound. Function
fAE {· | Θ} represents the processing of AE with the given
parameter Θ and pAE is a domain parameter used in AE.

Two types of data, normal and abnormal, were used as samples of the sound to be detected from a dataset published byCase Western Reserve University Bearing Data Center [19].
This dataset consists of vibration signal data from motor
bearings acquired using an acceleration sensor and also includes vibration data of normal bearings and abnormal bearings with scratches. The data used as normal were ”Normal 2” of ”Normal Baseline Data,” and the data used asabnormal were ”OR007@6 2” of ”48k Drive End Bearing
Fault Data.”
As the background noise, the sound recorded in theequipment room at the Ebina Oﬃce of Ricoh Co., Ltd., wasused, where several machines such as pumps were operatingin the same room. The total recording time was about 5.2 h(18770 s).
These data were mixed, and the experimental data wereprepared as shown in Table 1. To train the basis matrix F0 ,we used 1440-s-long background noise, which was a diﬀerent time period from the background noise sound for training the detector (AE). The background noise for the detection test includes the time periods used for training the basis
matrix and the detector. The volume of background noise isconstant in each dataset, namely, when the SNR is changed,
only the volume of the bearing is changed. Since the bearing vibration data have a length of only 10 s for both normal
and anomaly cases, the background noise data was divided
every 10 s, and each divided time period was mixed with thebearing vibration data.
AE for anomaly detection consists of four fully connected layers of the encoder and decoder. The input andoutput spectrograms are treated as one-dimensional vectors.
The number of nodes in each layer of the encoder is halved
in order from the input layer, and the decoder has the opposite configuration. Batch normalization [20] is performed on
the output of each layer. The ReLU function [21] was usedas the activation function. Other experimental conditions areshown in Table 2.

3.6 Detection Step
In the detection step, a sound signal x2 is recorded when
the target machine is operating in a noisy environment. Its
nonnegative spectrogram X2 is obtained by STFT, where2
X2 ∈ RΩ×T≥0 , and T 2 is the number of time frames. The
operation sound of the target machine is extracted from X2
by GGD-based SS-NMF with F0 , as follows:
X2 ≃ F0 G2 + H2 U2 ,S̃2 = H2 U2 ,

(30)
(31)

2
is the nonnegative spectrogram of a clean
where S̃2 ∈ RΩ×T
≥0
Ω×L
and U2 ∈operation sound of the target machine, H2 ∈ R≥0
L×T 2
R≥0 are respectively the basis and coeﬃcient matrices that
correspond to the operation sound of the target machine as
K×T 2
Eq. (31), and G2 ∈ R≥0
is the coeﬃcient matrix of F0 ,which corresponds to the background noise.
The anomaly score α is calculated as follows:M{} 2
1 ∑ ′(m)S2 − fAE S′(m)
α =
Θ ,
2
2
M m=1
(m) .

S′(m)
= S̃2
2

pAE
p

,

(32)
(33)

AE
where S̃2 ∈ RΩ×T
is the spectrogram piece cut out from≥0S̃2 every time frame number T AE . The threshold value ϕ isdetermined in advance, and if α > ϕ, it is judged that an
anomaly has occurred on the target machine.

(m)

4.

Evaluation

In this section, the anomaly detection performance of theExperimental Conditions

IEICE TRANS. INF. & SYST., VOL.E104–D, NO.3 MARCH 2021

446

4.2 Experimental Results and Considerations
Figure 2 shows the performance evaluation results of eachsystem. Here, the threshold ϕ for anomaly detection is experimentally determined to be the highest F-measure. In thecomparison of the cases without and with NMF preprocessTable 2

Experimental conditions

Sampling rate
Window length / shift length
FFT length
NMF
Number of bases
in F0 (K)
Number of basesin H1 and H2 (L)FeatureAE
Structure
Nodes
Epoch

32000 Hz10 ms / 5 ms
512 points (16 ms)
256
1
Power spectrogram
256 dimensions × 32 frames(pAE = 2, Ω = 256, T AE = 32)
Autoencoder
8192-4096-2048-1024-1024-2048-4096-8192
200

ing, the case with preprocessing shows a higher F-measure(with suitable parameter settings). Therefore, it is considered that noise reduction by NMF is an eﬀective preprocessing approach for anomaly detection.
For each case of without NMF, with ISNMF (p = 0.5),
and with GGD-NMF (ρ = 4.0 and p = 0.5), the precision–
recall curves under each SNR condition are shown in Fig. 3.
GGD-NMF has better characteristics than ISNMF, and it isconsidered that a superior anomaly detection model is generated. Since high performance is obtained in the case of
ρ > p, it is considered that the mitigation of the low-rank
property in NMF modeling is eﬀective for anomaly detection. This is because the mitigation of the low-rank property
in the NMF decomposition does not tolerate spectral outliers in the spectrogram, namely, GGD-NMF captures thedetailed spectral diﬀerences between the normal and abnormal bearing sounds. The mitigation of the low-rank property is derived from the generative model of sub-Gaussian
distribution, which corresponds to ρ > 2. To the best of
our knowledge, no sub-Gaussian-model-based NMF exists
other than GGD-NMF. Thus, the proposed system can beconsidered as the first practical application of sub-GaussianFig. 2F-measure values with various GGD-NMF parameters.

Fig. 3

Precision–Recall curves under each SNR condition.

AIBA et al.: NOISE ROBUST ACOUSTIC ANOMALY DETECTION SYSTEM WITH NMF BASED ON GENERALIZED GAUSSIAN DISTRIBUTION

447

Fig. 6

Mean power spectrum of sounds extracted by ISNMF (p = 0.5).

Fig. 4 Examples of power spectrograms of (a) normal bearing vibration,(b) mixed sound of (a) and background noise (SNR: −24 dB), (c) soundextracted by ISNMF (p = 0.5) from (b), and (d) sound extracted by GGDNMF (ρ = 4.0, p = 0.5) from (b).

Fig. 7 Mean power spectrum of sounds extracted by GGD-NMF (ρ =4.0, p = 0.5).

Fig. 5 Examples of power spectrograms of (a) abnormal bearing vibration, (b) mixed sound of (a) and background noise (SNR: −24 dB), (c)sound extracted by ISNMF (p = 0.5) from (b), and (d) sound extracted by
GGD-NMF (ρ = 4.0, p = 0.5) from (b).

NMF, which gives superior performance in an actual task.
In addition, the macro model of Eq. (8) with ρ > 2 has aversatility as described in Sect. 3.1. This property also contributed to the improvement of the proposed anomaly detection performance.
The eﬀect of GGD-NMF is analyzed from the actual
processed signals. Figure 4 shows the examples of normal spectrogram and Fig. 5 shows the examples of abnormal
spectrogram. It can be seen that GGD-NMF extracts features of the bearing vibration from the background noise as
well as ISNMF does; however, there exists slight diﬀerence
around 3 kHz between GGD-NMF and ISNMF in Fig. 5(d).
To analyze the diﬀerence in detail, we next show Figs. 6 and7 that depict the mean power spectra of the signals extracted
by ISNMF (p = 0.5) and GGD-NMF (ρ = 4.0, p = 0.5),
respectively. The larger the diﬀerence between the normal and anomaly spectra becomes, the more eﬀective thesource extraction process is for anomaly detection. The
spectra extracted by GGD-NMF show a greater diﬀerence
between the normal and abnormal cases in the band of ap-

proximately 3.0–3.3 kHz than those of ISNMF. This difference is likely to aﬀect the F-measure values. At the3187.5 Hz frequency bin, the diﬀerences between the normal and abnormal mean power spectra extracted by GGDNMF are greater than those of ISNMF by 5.5 dB (under
SNR = −12 dB), 2.4 dB (under SNR = −18 dB), and 0.6 dB(under SNR = −24 dB).
5. Conclusion
In this study, we aimed to improve the performance of theanomaly detection system in a low-SNR environment. For
this purpose, we proposed an anomaly detection system
based on the outlier detection method with target sound extraction based on GGD-NMF.
From the results of evaluation experiments, it was confirmed that the proposed system can improve the detection
performance in low-SNR environments compared with asystem without preprocessing. In addition, the mitigation of
the low-rank property in NMF resulted in a decomposition
more suitable for anomaly detection.
Acknowledgments
This work was partly supported by JSPS KAKENHI Grant
Number JP19H01116.
References[1] V. Chandola, A. Banerjee, and V. Kumar, “Anomaly detection: A
survey,” ACM Comput. Surv., vol.41, no.3, July 2009.

IEICE TRANS. INF. & SYST., VOL.E104–D, NO.3 MARCH 2021

448

[2] D. Ramotsoela, A. Abu-Mahfouz, and G. Hancke, “A survey ofanomaly detection in industrial wireless sensor networks with critical water system infrastructure as a case study,” Sensors, vol.18,
no.8, p.2491, 2018.
[3] K. Singh and S. Upadhyaya, “Outlier detection: applications andtechniques,” Int. J. Computer Science Issues, vol.9, no.1, pp.307–
323, 2012.
[4] J. An and S. Cho, “Variational autoencoder based anomaly detection
using reconstruction probability,” SNU Data Mining Center, 2015-2
Special Lecture on IE, 2015.
[5] H. Qiu, J. Lee, J. Lin, and G. Yu, “Wavelet filter-based weak signature detection method and its application on rolling element bearing
prognostics,” J. Sound and Vibration, vol.289, no.4, pp.1066–1090,
Feb. 2006.
[6] N. Kitaoka, I. Akahori, and S. Nakagawa, “Speech recognition
under noisy environments using spectral subtraction with smoothing of time direction and real-time cepstral mean normalization,”
Proc. International Workshop on Hands-Free Speech Communication, pp.159–162, 2001.
[7] D.D. Lee and H.S. Seung, “Learning the parts of objects by nonnegative matrix factorization,” Nature, vol.401, no.6755, p.788, Oct.
1999.
[8] X. lu, Y. Tsao, S. Matsuda, and C. Hori, “Speech enhancement basedon deep denoising auto-encoder,” Proc. Interspeech, pp.436–440,
2013.
[9] D. Wang and J. Chen, “Supervised speech separation based on deeplearning: An overview,” IEEE/ACM Trans. Audio, Speech, Language Process., vol.26, no.10, pp.1702–1726, Oct. 2018.
[10] H. Shao, H. Jiang, F. Wang, and H. Zhao, “An enhancement
deep feature fusion method for rotating machinery fault diagnosis,”
Knowledge-Based Systems, vol.119, pp.200–220, March 2017.
[11] D. Kitamura, “Nonnegative matrix factorization based on complex
generative model,” Acoustical Science and Technology, vol.40, no.3,
pp.155–161, 2019.
[12] P. Smaragdis, B. Raj, and M. Shashanka, “Supervised and semisupervised separation of sounds from single-channel mixtures,”Proc. the 7th Int. Conf. Independent Component Analysis and Signal
Separation, pp.414–421, 2007.
[13] D. Kitamura, H. Saruwatari, K. Yagi, K. Shikano, Y. Takahashi,
and K. Kondo, “Music signal separation based on supervised nonnegative matrix factorization with orthogonality and maximumdivergence penalties,” IEICE Trans. Fundamentals, vol.97-A, no.5,
pp.1113–1118, May 2014.
[14] C. Févotte, N. Bertin, and J.L. Durrieu, “Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis,” Neural Computation, vol.21, no.3, pp.793–830, March
2009.
[15] D.D. Lee and H.S. Seung, “Algorithms for non-negative matrix factorization,” Advances in Neural Information Processing Systems,
vol.13, pp.556–562, 2000.
[16] S. Mogami, N. Takamune, D. Kitamura, H. Saruwatari, Y.
Takahashi, K. Kondo, and N. Ono, “Independent low-rank matrix
analysis based on time-variant sub-gaussian source model for determined blind source separation,” IEEE/ACM Trans. Audio, Speech,
Language Process., vol.28, pp.503–518, Jan. 2019.
[17] K. Kamo, Y. Kubo, N. Takamune, D. Kitamura, H. Saruwatari,
Y. Takahashi, and K. Kondo, “Joint-diagonalizability-constrained
multichannel nonnegative matrix factorization based on multivariate complex sub-gaussian distribution,” Proc. EUSIPCO2020, 2020.
(arXiv preprint arXiv:2007.00416).
[18] D. Kitamura, S. Mogami, Y. Mitsui, N. Takamune, H. Saruwatari, N.
Ono, Y. Takahashi, and K. Kondo, “Generalized independent lowrank matrix analysis using heavy-tailed distributions for blind source
separation,” EURASIP J. Adv. Signal Process., vol.2018, no.1, Article No. 28, May 2018.
[19] Case Western Reserve University, “Bearing Data Center,” http://
csegroups.case.edu/bearingdatacenter/home, accessed May. 12. 2017.

[20] S. Ioﬀe and C. Szegedy, “Batch normalization: Accelerating deep
network training by reducing internal covariate shift,” arXiv preprint
arXiv:1502.03167, 2015.
[21] X. Glorot, A. Bordes, and Y. Bengio, “Deep sparse rectifier neuralnetworks,” Proc. the Fourteenth Int. Conf. Artificial Intelligence andStatistics, pp.315–323, 2011.

Akihito Aiba
received the B.E. and M.E. degrees from Tohoku University in 2008 and 2010,
respectively. He joined Ricoh Company, Ltd.
in 2010. His research interests include acoustic signal processing and machine learning. He
is a member of the Acoustical Society of Japan.

Minoru Yoshidareceived the M.E. degreefrom Tokyo Denki University in acoustical digital signal processing in 1993. He joined Ricoh
Company, Ltd. in 2015. His research interests
are digital signal processing for outlier detection.

Daichi Kitamura
received the Ph.D. degreefrom SOKENDAI, Hayama, Japan. He joined
the University of Tokyo, Tokyo, Japan, in 2017as a Research Associate, and he moved to the
National Institute of Technology, Kagawa Collage, Takamatsu, Japan, as an Assistant Professor, in 2018. His research interests include audio
source separation, statistical signal processing,and machine learning. He was the recipient of
the Awaya Prize Young Researcher Award fromThe Acoustical Society of Japan (ASJ) in 2015,
Ikushi Prize from Japan Society for the Promotion of Science in 2017,
Best Paper Award from IEEE Signal Processing Society Japan in 2017,
and Itakura Prize Innovative Young Researcher Award from ASJ in 2018.

Shinnosuke Takamichi
received the B.E.
degree from Nagaoka University of Technology, Nagaoka, Japan, in 2011, and the M.E. and
Ph.D. degrees from the Graduate School of Information Science, Nara Institute of Science andTechnology, Ikoma, Japan, in 2013 and 2016,
respectively. He is currently an Assistant Professor at The University of Tokyo. He has received more than ten paper/achievement awards
including the 3rd IEEE Signal Processing Society Japan Young Author Best Paper Award. He
is a member of ASJ, IPSJ, ISCA, and IEEE SPS.

AIBA et al.: NOISE ROBUST ACOUSTIC ANOMALY DETECTION SYSTEM WITH NMF BASED ON GENERALIZED GAUSSIAN DISTRIBUTION

449

Hiroshi Saruwatari
Hiroshi Saruwatarireceived the B.E., M.E., and Ph.D. degrees
from Nagoya University, Nagoya, Japan, in
1991, 1993, and 2000, respectively. He joined
SECOM Intelligent Systems Laboratory, Tokyo,
Japan, in 1993, and the Nara Institute of Science
andTechnology, Ikoma, Japan, in 2000. Since
2014, he has been a Professor with the University of Tokyo, Tokyo. His research interests include statistical audio signal processing, blindsource separation, and speech enhancement. He
was the recipient of the paper awards from IEICE in 2001 and 2006,
fromTAF in 2004, 2009, 2012, and 2018, from IEEE-IROS2005 in 2006,
and from APSIPA in 2013 and 2018. He was also the recipient of the
DOCOMO Mobile Science Award in 2011, Ichimura Award in 2013, The
Commendation for Science and Technology by the Minister of Education in2015, Achievement Award from IEICE in 2017, and Hoko-Award in 2018.
He has been professionally involved in various volunteer works for IEEE,EURASIP, IEICE, and ASJ.