# Sparsity-based Algorithm for Detecting Faults in Rotating
Machines111Preprint submitted to Mechanical Systems and Signal Processing
(Elsevier).

Wangpeng He Email<EMAIL_ADDRESS>Tandon School of Engineering, New York
University, 6 Metrotech Center, Brooklyn, NY 11201, USA  State Key Laboratory
for Manufacturing and Systems Engineering, School of Mechanical Engineering,
Xi’an Jiaotong University, Xi’an 710049, P. R. China  Yin Ding Email:
<EMAIL_ADDRESS>Tandon School of Engineering, New York University, 6 Metrotech
Center, Brooklyn, NY 11201, USA  Yanyang Zi  State Key Laboratory for
Manufacturing and Systems Engineering, School of Mechanical Engineering, Xi’an
Jiaotong University, Xi’an 710049, P. R. China  Ivan W. Selesnick  Tandon
School of Engineering, New York University, 6 Metrotech Center, Brooklyn, NY
11201, USA

###### Abstract

This paper addresses the detection of periodic transients in vibration signals
for detecting faults in rotating machines. For this purpose, we present a
method to estimate periodic-group-sparse signals in noise. The method is based
on the formulation of a convex optimization problem. A fast iterative
algorithm is given for its solution. A simulated signal is formulated to
verify the performance of the proposed approach for periodic feature
extraction. The detection performance of comparative methods is compared with
that of the proposed approach via RMSE values and receiver operating
characteristic (ROC) curves. Finally, the proposed approach is applied to
compound faults diagnosis of motor bearings. The non-stationary vibration data
were acquired from a SpectraQuest’s machinery fault simulator. The processed
results show the proposed approach can effectively detect and extract the
useful features of bearing outer race and inner race defect.

##  1 Introduction

Rotating machinery is one of the most common types of mechanical equipment and
plays a significant role in industrial applications. Early detection of faults
developing in rotating machinery is of great importance to prevent economic
loss and personal casualties [1]. Rolling element bearings and gearboxes are
two kinds of widely used components in rotating machines and their failures
are among the most frequent reasons for machine breakdown.

Much attention has focused on vibration-based diagnosis of mechanical faults
in rotating machines [2]. The detection of periodically occurring transient
vibration signatures is of vital importance for vibration-based condition
monitoring and fault detection of rotating machinery [3]. However, these
useful transient features are usually buried in heavy background noise and
other irrelevant vibrations. To address this problem, many signal processing
methods have been introduced, such as singular value decomposition (SVD) [4],
empirical mode decomposition (EMD) [5], and methods based on different wavelet
transforms, e.g., dual-tree wavelet in [6], harmonic wavelet in [7], and
tunable Q-factor wavelet (TQWT) in [8]. These methods have achieved successful
applications in the field of machinery fault diagnosis.

Recently, an algorithm, called ‘overlapping group shrinkage’ (OGS) was
developed for estimating group-sparse signals in noise [9]. The OGS algorithm
was initially formulated as a convex optimization promoting group sparsity by
a convex regularization. In order to promote sparsity more strongly, an
improved version of OGS was developed, which utilizes a non-convex
regularization [10]. The superiority of denoising group-sparse signals using
the approach presented in [10] indicates its potential for effectively
extracting periodic transient pulses.

This paper aims to develop an approach for rotating machinery fault diagnosis
based on a periodic group-sparse signal representation. The signature of
localized faults of the gear teeth and bearing components generally exhibit
periodic transient pulses when a rotating machine is operated at a constant
speed [11]. Meanwhile, the large amplitudes of these useful features are not
only sparse but also tend to form groups. Several neighborhood-based denoising
methods have been developed for machinery fault diagnosis utilizing this
property [12, 13, 14]. Our approach is based on a signal model intended to
capture the useful impulsive features for machinery fault diagnosis. In
particular, this paper addresses the problem of estimating $x$ from a noisy
observation $y$. We model the measured discrete-time series, $y$, as

$y_{n}=x_{n}+w_{n},\qquad n=0,\ldots,N-1,$ (1)

where the signal $x$ is known to have a periodic group-sparse property and $w$
is white Gaussian noise. A group-sparse signal is one where large magnitude
signal values tend not to be isolated. Instead, these large magnitude values
tend to form groups.

Convex optimization is commonly used to estimate sparse vectors from noisy
signals, where we solve the optimization problem with the prototype

$x^{\mathsf{*}}=\arg\min_{x}\Big{\\{}F(x)=\displaystyle\frac{1}{2}\lVert
y-x\rVert_{2}^{2}+\lambda\Phi(x)\Big{\\}},$ (2)

where $\lambda>$ 0 is a regularization parameter and $\Phi$ : $\mathbb{R}^{N}$
$\to$ $\mathbb{R}$ is a sparsity-promoting penalty function (regularizer).
Conventionally, the regularizer $\Phi(x)$ is a convex function, e.g.
$\ell_{1}$-norm. In [15], an idea of using non-convex regularizer and keeping
the convexity of entire problem is used for signal denoising problem, where
the sparsity can be significantly promoted comparing to $\ell_{1}$-norm, and
the problem still has a unique solution due to the convexity.

In this paper, we adopt ideas from [15] and [9], and present a method for
estimating periodic-group-sparse signals in noise. We propose its use for
detecting faults in rotating machinery, where the fault characteristic
frequency (period of the group-sparse pulses) is used as prior knowledge.
Similar to the approach in [10], the non-convex regularization term in the
proposed method is properly chosen so as to ensure that the total objective
function $F$ is convex; however, in contrast to [10], where each group has to
be contiguous, we allow grouping with intervals, and furthermore periodically.

As a consequence, in this work, the regularization term $\Phi$ in (2) is
formulated specifically to utilize the periodicity of the impulsive fault
features. The aim of our approach is to capture the useful impulsive features
for the purpose of machinery fault diagnosis. Additionally, it has the
potential to separate compound fault features by utilizing different periods
of the periodic transient pulses corresponding to different fault frequencies
(e.g., various defect frequencies of rolling element bearings). The proposed
approach also reduces to a non-periodic group-sparse signal denoising method,
i.e., we can utilize the sparsity-based OGS approach [10], if we do not have
prior knowledge of the period of the group-sparse transients. Thus, the
proposed sparsity-based approach is a generalization of the non-convex
regularized OGS [10].

The paper is organized as follows. A brief review of OGS with convex and non-
convex regularization is given in Section 2. Section 3 presents a method for
denoising periodic group sparse signals. In Section 4 a simulation study is
performed to validate the effectiveness of the proposed method. Section 5
applies the proposed periodic group sparse denoising method to fault diagnosis
of motor bearings for further validation of its effectiveness. Finally,
conclusions are summarized in Section 6.

##  2 Review

In this section, we give short reviews of overlapping group shrinkage (OGS)
[9] and majorization-minimization (MM) [16].

###  2.1 Overlapping Group Shrinkage (OGS)

There are several advantages to formulating sparse estimation as a convex
optimization problem. The most basic advantage is that the problem can then be
very reliably and efficiently solved using convex optimization methods [17].
Although a non-convex regularizer can promote sparsity more strongly, it
generally leads to a non-convex optimization problem with non-optimal local
minima [18]. To avoid the formulation of a non-convex optimization problem,
one may utilize a non-convex regularizer $\Phi$ designed so as to ensure the
total objective function is convex.

The problem of denoising a group sparse signal was addressed in [9] which
utilized convex optimization. An improved method was proposed in [10], which
utilized non-convex regularization designed to ensure convexity of the
objective function. The problem is solved efficiently by an iterative
algorithm based on majorization-minimization (MM) [16]. The objective function
for the OGS problem, with a group size of $K$ and convex or non-convex
regularization, is denoted as

$\displaystyle
x^{\mathsf{*}}=\arg\min_{x\in{}^{N}}\Bigg{\\{}P_{0}(x)=\displaystyle\frac{1}{2}\lVert
y-x\rVert_{2}^{2}+\lambda\sum_{n}\phi\bigg{(}\Big{[}\sum_{k\in\mathcal{K}}x_{n+k}^{2}\Big{]}^{1/2};a\bigg{)}\Bigg{\\}}$
(3)

where $\mathcal{K}:=\\{0,1,\ldots K-1\\}$, and $y\in{}^{N}$ is the noisy
observation. For $x\in{}^{N}$, we define $x_{n}=0$ for $n<0$ and $n\geqslant
N$. We assume the penalty function $\phi:\real\to{}_{+}$ satisfies the
properties:

  1. 1.

$\phi$ is continuous on .

  2. 2.

$\phi$ is twice differentiable on $\real\setminus\\{0\\}$.

  3. 3.

$\phi$ is increasing and concave on +.

  4. 4.

$\phi$ is symmetric, $\phi(-x;a)=\phi(x;a)$.

  5. 5.

$\phi^{\prime}(0^{+};a)=1$.

  6. 6.

$\phi^{\prime\prime}(x;a)\geqslant-a$ for all $x\neq 0$.

  7. 7.

$\phi(x;0)=|x|$.

which is used to induce the resulting sparsity in an optimization problem.

Figure 1: Examples of penalty functions. Table 1: Sparsity-promoting penalty
functions $(a\geqslant 0)$.

Penalty | $\phi(x;a)$ | $\psi(x)$
---|---|---
abs | $\left|x\right|$ | $\left|x\right|$
log | $\displaystyle\frac{1}{a}\log(1+a\left|x\right|)$ | $\displaystyle\left|x\right|(1+a\left|x\right|)$
rat | $\displaystyle\frac{\left|x\right|}{1+a\left|x\right|/2}$ | $\displaystyle\left|x\right|(1+a\left|x\right|)^{2}$
atan | $\displaystyle\frac{2}{a\sqrt{3}}\left(\tan^{-1}\left(\frac{1+2a\left|x\right|}{\sqrt{3}}\right)-\frac{\pi}{6}\right)$ | $\displaystyle\left|x\right|(1+a\left|x\right|+a^{2}\left|x\right|^{2})$

Note that the objective function $P_{0}$ in (3) may be convex even if the
regularizer (penalty) is not. When the penalty function $\phi$ satisfies the
conditions listed above, then the parameter ‘$a$’ can be chosen so that
$P_{0}$ is convex even if $\phi$ is not [15]. This approach is used to improve
OGS in [10] where it has been proved that the objective function $P_{0}$ in
(3) is strictly convex if

$\displaystyle 0\leqslant a<\frac{1}{\lambda K}.$ (4)

Table 1 gives some examples of penalty functions satisfying the above-listed
conditions, and these examples are illustrated in Fig. 1. In order to induce
the sparsity most strongly, the arctangent ‘atan’ function can be used among
all the three given non-convex functions. Fig. 1 also illustrates that under
the same ‘$a$’ parameter, ‘atan’ function adheres the most ‘concavity’. In
addition to these non-convex penalty functions, Fig. 1 also shows the
$\ell_{1}$-norm as a special case. Note that we allow $a=0$, which is the
extreme case of the three other penalties. In particular,
$\phi(x;0)=\left|x\right|$ is a convex function. The original OGS algorithm
given in [9] considers only the special case $a=0$, i.e., convex
regularization.

###  2.2 Majorization-minimization (MM)

A detailed derivation to solve problem (3), based on the majorization-
minimization (MM) method is given in [10]. The MM method simplifies a
complicated optimization problem into a sequence of easy ones, and is
described by the iteration

$\displaystyle u^{(i+1)}=\arg\min_{u}G(u,u^{(i)}),$ (5)

where $i$ denotes the number of iterations, and $G:{}^{N}\times{}^{N}\to\real$
is a majorizer of the objective function $J$, satisfying

$\displaystyle G(u,v)\geqslant J(u),\quad\text{ for all }u,v,$ (6a)
$\displaystyle G(u,u)=J(u).$ (6b)

A majorizer of $\phi$ is given by

$\displaystyle g(u,v)=\frac{(u^{2}-v^{2})}{2\psi(v)}+\phi(v;a),\quad\text{
when }v\neq 0,$ (7)

where $\psi$ is the function given in Table 1. For $v=0$, $g(u,0)$ is defined
by

$\displaystyle g(u,0):=\begin{cases}+\infty,&\text{ if }u\neq 0,\\\ 0,&\text{
if }u=0.\end{cases}$ (8)

As a consequence, the function $g:\real\times\real\to{}_{+}$ satisfies

$\displaystyle g(u,v)$ $\displaystyle\geqslant\phi(u;a),\quad\text{ when
}u\neq v,$ (9a) $\displaystyle g(u,u)$ $\displaystyle=\phi(u;a),$ (9b)

Note that $g(u,0)$ defined in (8) equals infinity except $u=0$. This forces
its minimizer to lock to $u=0$ in the MM iteration described in (5). This
issue in the OGS problem is discussed in [10], where it does not affect the
convergence when the algorithm is implemented with a non-zero initialization.

##  3 OGS with binary weights

To facilitate the following derivation, we define a binary sequence
$b=[b_{0},b_{1},\ldots b_{K-1}]$, with $b_{k}\in\\{0,1\\}$, and sets

$\displaystyle\mathcal{K}:=\\{0,1,\ldots,K-1\\},$ (10a)
$\displaystyle\mathcal{K}_{0}:=\\{k\in\mathcal{K}:b_{k}=0\\},$ (10b)
$\displaystyle\mathcal{K}_{1}:=\\{k\in\mathcal{K}:b_{k}=1\\}.$ (10c)

Since $b$ is a binary vector, we have
$\mathcal{K}=\mathcal{K}_{0}\cup\mathcal{K}_{1}$ and
$\mathcal{K}_{0}\cap\mathcal{K}_{1}=\emptyset$ is the empty set. We denote the
cardinality (size) of the sets $\mathcal{K}$, $\mathcal{K}_{1}$ and
$\mathcal{K}_{0}$ as $K$, $K_{1}$, and $K_{0}$, respectively, so that
$K=K_{0}+K_{1}$, and

$\displaystyle\sum_{k\in\mathcal{K}}b_{k}=\sum_{k\in\mathcal{K}_{1}}b_{k}$
$\displaystyle=K_{1}.$ (11)

The following proposition is straightforward.

######  Proposition 1.

Let $\phi:\real\to\real$ satisfy the conditions listed in Section 2. When
$\gamma>0$ and $\lambda>0$, the function $p:\real\to\real$,

$\displaystyle p(v)=\frac{\gamma}{2}v^{2}+\lambda\phi(v\;a)$ (12)

is strictly convex if

$\displaystyle\phi^{\prime\prime}(v;a)>-\frac{\gamma}{\lambda}.$ (13)

Detailed proof of this proposition can be found in Appendix A.

###  3.1 Problem definition

We define the objective function $P_{1}:{}^{N}\to\real$ as

$\displaystyle
x^{\mathsf{*}}=\arg\min_{x}\Bigg{\\{}P_{1}(x)=\displaystyle\frac{1}{2}\lVert
y-x\rVert_{2}^{2}+\lambda\sum_{n}\phi\big{(}\theta(x,b,n);a\big{)}\Bigg{\\}}$
(14)

where the binary-weighted grouping function
$\theta:{}^{N}\times{}^{K}\times\mathbb{Z}\to\real$ is defined as

$\displaystyle\theta(x,b,n):=\Big{[}\sum_{k=0}^{K-1}b_{k}x_{n+k}^{2}\Big{]}^{1/2},$
(15)

which is the Euclidean norm of a binary weighted block. For $x\in{}^{N}$, we
define $B_{n}(x)\in{}^{K}$ as

$\displaystyle B_{n}(x):=[x_{n},x_{n+1},\ldots x_{n+K-1}],$ (16)

i.e., a $K$-point subvector of $x$, starting at index $n$. Hence
$\theta(x,b,n)$ can be written as

$\displaystyle\theta(x,b,n)=\lVert b\odot B_{n}(x)\rVert_{2},$ (17)

where $\odot$ denotes element-wise multiplication.

Note that in (14) if $b_{k}=1$ for all $k\in\mathcal{K}$, then
$\theta^{2}(x,b,n)=\sum_{k}x_{n+k}^{2}$, and problem (14) reduces to (3).
Therefore (3) is a special case of (14). Moreover, if $K=K_{1}=1$, this
problem further reduces to scalar (i.e., non-group) sparse denoising. In the
following discussion, we consider the group sparse case: $N\gg K\geqslant
K_{1}>0$. The case $K_{1}=0$ is trivial.

In the following section, we exploit the convexity condition of (14), to
constrain the non-convex penalty function $\phi$, so as to ensure that the
objective function $P_{1}$ is convex. Therefore, the problem formulation we
ultimately propose is a convex optimization problem. Hence, it will have no
non-optimal local minima in which an iterative optimization algorithm can be
trapped.

###  3.2 Convexity conditions

######  Proposition 2.

Let $b\in\\{0,1\\}^{K}$ be a binary vector with $\sum_{k}b_{k}=K_{1}$. The
function $P:{}^{K}\to\real$ defined as

$\displaystyle
P(u)=\frac{1}{2K_{1}}\sum_{k=0}^{K-1}b_{k}u_{k}^{2}+\lambda\phi\bigg{(}\Big{[}\sum_{k=0}^{K-1}b_{k}u_{k}^{2}\Big{]}^{1/2};a\bigg{)},$
(18)

is strictly convex if

$\displaystyle\phi^{\prime\prime}(x;a)>-\frac{1}{K_{1}\lambda}\qquad\text{ for
all }x\neq 0.$ (19)

This Proposition is proved in detail in Appendix B. Utilizing the above
results, we find a range of $a$ ensuring that, even though $\phi$ is non-
convex, the objective function $P_{1}$ in (14) is strictly convex.

######  Theorem 1.

Assume that problem $P_{1}$ (14) has $K_{1}=\sum_{k}b_{k}$ non-zero weights in
one binary weighted group, then the objective function is strictly convex if
the parameter ‘$a$’ of the penalty function $\phi(\cdot;a)$ satisfies

$\displaystyle 0\leqslant a<\frac{1}{K_{1}\lambda}.$ (20)

###### Proof.

First, since $\sum_{k}b_{k}=K_{1}$ and $b_{k}\geqslant 0$, $x_{n}^{2}\geqslant
0$, it follows that

$\displaystyle K_{1}\sum_{n}x_{n}^{2}$
$\displaystyle=\sum_{k}b_{k}\sum_{n}x_{n}^{2}=\sum_{n}\sum_{k}b_{k}x_{n}^{2}.$
(21)

Using (21), we write
$\displaystyle\frac{1}{2}\sum_{n}x_{n}^{2}=\frac{1}{2K_{1}}\sum_{n}\sum_{k}b_{k}x_{n}^{2}$.
Therefore, the data fidelity term in problem (14) can be written as

$\displaystyle F(x)$ $\displaystyle=\displaystyle\frac{1}{2}\lVert
y-x\rVert_{2}^{2}=\displaystyle\frac{1}{2}\sum_{n}x_{n}^{2}+L(x)$
$\displaystyle=\frac{1}{2K_{1}}\sum_{n=0}^{N-K}\Big{(}\sum_{k}b_{k}x_{n+k}^{2}\Big{)}+L(x),$
(22)

where $L(x)$ is linear in $x$. Adding $L(x)$ to a strictly convex function
yields a strictly convex function.

Using the above results, the objective function in problem (14) can be
reorganized as

$\displaystyle P_{1}(x)$
$\displaystyle=\sum_{n=0}^{N-K}\Bigg{[}\frac{1}{2K_{1}}\sum_{k}b_{k}x_{n+k}^{2}+\lambda\phi\bigg{(}\Big{[}\sum_{k}b_{k}x_{n+k}^{2}\Big{]}^{1/2};a\bigg{)}\Bigg{]}+L(x),$
$\displaystyle=\sum_{n=0}^{N-K}P(B_{n}(x))+L(x),$ (23)

where $B_{n}(x)$ is defined in (16). As a consequence, if $P$ is strictly
convex, then $P_{1}$ is strictly convex. The condition for convexity of $P$ is
given in Proposition 2. Hence, as long as the inequality condition (19) is
satisfied, $P$ is strictly convex. Moreover, $\phi$ satisfies condition
$\phi^{\prime\prime}(x;a)\geqslant-a$ (condition 6 in Section 2). This implies
that when (20) is satisfied, $P$ is strictly convex, and the entire objective
function $P_{1}$ is convex. ∎

Note that Theorem 1 generalizes the convexity condition of OGS in [10,
Corollary 2]. When every element in binary vector $b$ equals $1$ then
$K_{1}=K$ and Theorem 1 reduces to Corollary 2 in [10] as a special case.

We have proved that under a more flexible group structure (binary weights),
non-convex penalty functions can be utilized to promote structured sparsity,
and the convexity of the objective function is preserved when the
regularization parameter $a$ is suitably set. Moreover, the result also shows
that when maximizing the non-convexity of the penalty function, only the the
nonzero weights matter for the selection of the parameter $a$.

###  3.3 Algorithm derivation

To minimize $P_{1}$ using the MM procedure, we define a majorizer
$G:{}^{N}\times{}^{N}\to\real$, defined by

$\displaystyle G(x,v)$ $\displaystyle=\displaystyle\frac{1}{2}\lVert
y-x\rVert_{2}^{2}+\lambda\sum_{n}g\big{(}\theta(x,b,n),\theta(v,b,n)\big{)}.$
(24)

We can write $G$ as

$\displaystyle G(x,v)$ $\displaystyle=\displaystyle\frac{1}{2}\lVert
y-x\rVert_{2}^{2}+\frac{\lambda}{2}\sum_{n}\frac{1}{\psi(\theta(v,b,n))}\theta^{2}(x,b,n)+C$
(25a) $\displaystyle=\displaystyle\frac{1}{2}\lVert
y-x\rVert_{2}^{2}+\frac{\lambda}{2}\sum_{n}\sum_{k}\frac{b_{k}}{\psi(\theta(v,b,n))}x^{2}_{n+k}+C$
(25b) $\displaystyle=\displaystyle\frac{1}{2}\lVert
y-x\rVert_{2}^{2}+\frac{\lambda}{2}\sum_{n}[r(v)]_{n}x_{n}^{2}+C$ (25c)

where $r(v)\in{}^{N}$ is defined as

$\displaystyle[r(v)]_{n}:=\displaystyle\sum_{j=0}^{K-1}\frac{b_{j}}{\psi\big{(}\theta(v,b,n-j)\big{)}}.$
(26)

Then $G(x,v)$ can be written as

$\displaystyle G(x,v)$
$\displaystyle=\displaystyle\frac{1}{2}\sum_{n}x_{n}^{2}+\frac{\lambda}{2}\sum_{n}r_{n}x_{n}^{2}-\sum_{n}y_{n}x_{n}+C(y)$
(27)
$\displaystyle=\sum_{n}\Big{(}\displaystyle\frac{1}{2}+\frac{\lambda[r(v)]_{n}}{2}\Big{)}x_{n}^{2}-y_{n}x_{n}+C(y)$
(28)

which has an explicit minimizer $x_{n}={y_{n}}/({1+\lambda[r(v)]_{n}})$.
Hence, the MM iteration (5) is given by

$\displaystyle x_{n}^{(i+1)}=\frac{y_{n}}{1+\lambda[r(x_{n}^{(i)})]_{n}}.$
(29)

Table 2 gives the explicit steps to solve (14), assuming the penalty function
is chosen from Table 1 and satisfies (20). This guarantees the problem (14) is
strictly convex and consequently MM procedure (5) will converge to the unique
global minimizer.

Note that zero-locking might occur when using quadratic function to majorize
non-smooth function [19]. For the OGS problem, initializing the algorithm by
$x^{(0)}=y$ avoids this issue; a detailed proof is given in [10, Lemma B].
This lemma is not affected by introducing binary group weights. In other
words, if the function $F$ in [10, Lemma B] is substituted with $P_{1}$ (14),
the derivation is still valid with an almost identical proof. Moreover, since
we allow $\psi(x)$ to be 0, Equation (26) may lead to a ‘divide-by-zero’. The
work of [10] contains a sufficient discussion that this problem is avoided by
the initialization $x^{(0)}=y$, based on a same lemma. Hence, there is no
zero-locking or ‘divide-by-zero’ issue when solving the problem (14) by the
algorithm in Table 2.

Table 2: OGS with binary weights.
$\displaystyle\text{Input:}~{}y\in{}^{N},~{}\lambda,~{}b\in\\{0,1\\}^{K}$
$\displaystyle\text{Initialization:}~{}x=y,~{}S=\\{n:y_{n}\neq 0\\}$ Repeat
for $n\in S$: $\displaystyle\qquad
r_{n}=\sum_{j=0}^{K-1}\frac{b_{j}}{\psi\big{(}\theta(x,b,n-j)\big{)}}$
$\displaystyle\qquad x_{n}=\frac{y_{n}}{1+\lambda r_{n}}$ $\displaystyle\qquad
S=\\{n:\left|x_{n}\right|>\epsilon\\}$ Until convergence
$\displaystyle\text{Return: }x$

###  3.4 Periodicity-induced OGS (POGS)

In the previous sections, we have given a method for group sparse denoising
with binary weights within the group. In signal model (1), since the
periodicity of impulsive faults in $x$ is assumed to be approximately
consistent over a reasonable duration, the time interval between two
consecutive faults can be considered identical within the support of a group.
Moreover, when the period $T$ of a potential fault is known or predictable
from the knowledge of the machinery, we can select the group with a length $K$
and its zero and non-zero entries by

$\displaystyle N_{0}+N_{1}$ $\displaystyle\approx f_{s}T,$ (31a)
$\displaystyle N_{0}+N_{1}$ $\displaystyle=K/M,$ (31b)

where $N_{1}$ and $N_{0}$ are the estimated length (in samples) of impulsive
transients and the time interval between them, and integer $M\geqslant 2$ is
the number of periods contained in one group, and $f_{s}$ is sampling rate.
Thus, in one group, the numbers of zero and non-zero entries are
$K_{0}=MN_{0}$ and $K_{1}=MN_{1}$, respectively.

Moreover, when the transient sequence is periodical, the binary weight $b_{k}$
are chosen according to a periodic group structure. Specifically, under the
constraint (31), $b\in\\{0,1\\}^{K}$ is given by

$\displaystyle b=[\ \underbrace{1,1,\ldots,1}_{N_{1}},\
\underbrace{0,0,\ldots,0}_{N_{0}},\ \underbrace{1,1,\ldots,1}_{N_{1}},\
\underbrace{0,0,\ldots,0}_{N_{0}},\ \ldots,\ \
\underbrace{1,1,\ldots,1}_{N_{1}},\ \underbrace{0,0,\ldots,0}_{N_{0}}\ ],$
(32)

where in each period, there are $N_{1}$ non-zero entries grouped according to
the impulsive signal, and the entire group comprises $M$ periods. In this
case, the last $N_{0}$ zeros in (32) has no effect in problem (14) by the
definition (15). In practice, we omit trailing zeroes and the actual length of
$b$ involved in computation is $K-N_{0}$. Note that although we use the
parameters such as $K,K_{0},K_{1},N_{0},N_{1}$ and $M$ to illustrate the
binary pattern structure, in fact given $f_{s}$ and $T$, we only need to
select $N_{1}$ and $M$, then use (31), the pattern of $b$ (32) can be
determined.

Consequently, we propose to recover a periodic impulsive signature from a
noisy observation by solving problem (14) with $b$ defined by (31) and (32).
We refer to this method as Periodicity-induced OGS (POGS). This is an
extension of OGS accounting for the periodicity of the sparse signal.

As a special case, if the period $T$ is not known, we may use conventional OGS
(3) to detect faults. If a period $T$ can be determined by inspecting the
output produced by OGS, then an enhanced result with better accuracy may be
achieved by POGS using the determined periodicity.

##  4 Simulation validation

Figure 2:  Example 1: Simulated signal: (a) clean data and (b) noisy data.

Figure 3:  Example 1: Denoising results. (a) Proposed method with
$N_{1}=4,M=4$. (b) Proposed method with $N_{1}=2,M=4$. (c) Conventional OGS.
(d) Wavelet-based denoising.

To test the proposed method, we apply it to the simulated data illustrated in
Fig. 2(a). The signal is a 1-second signal with sampling rate $f_{s}=6400$ Hz,
and is composed of a periodic sequence of transients occurring with $80$ Hz.

We simulate the vibration signal containing features caused by machinery
defect as a sequence of impulsive transients, and each transient consists for
10 samples (1.6 ms when $f_{s}=6400$ Hz). In this example, each transient is
composed of a random number of sinusoidal components, each with random
frequency and random initial phase. More specifically, each transient can be
written as

$\displaystyle v(n)=\sum_{i=1}^{U}A_{i}\sin(\omega_{i}n+\beta_{i}),\quad
n=0,1,2\ldots 9,$ (33)

where $1\leqslant U\leqslant 10$ is a random integer, and for each $i$,
$A_{i}$ is a random amplitude, and $\omega_{i}$ is a random frequency, and
$\beta_{i}$ is a random initial phase. The sequence of transients is shown in
Fig. 2(a), and we show the detail of one transient at about $t=0.40$ second in
the box. The generated test signal is multiplied by a constant, so that it has
unit standard deviation.

In order to evaluate the false alarm rate, the first part of the test signal
contains no transient. Fig. 2(a) shows the clean test signal, where there are
50 faults starting at approximately $t=0.36$ second with a period $T=1/80$
second. White Gaussian noise with standard deviation $\sigma=2.5$ is added to
the simulated fault sequence, as illustrated in Fig. 2(b).

Figure 4: Example 1: Binary pattern $b$ when $M=4$, $N_{1}=4$.

When the periodicity of the faults is known, we can use POGS, where (31) and
(32) are used to define the binary weight vector $b$. In this example, we set
$N_{1}=4$ and $M=4$. Since $T=1/80$ seconds, we can calculate
$N_{0}=6400/80-N_{1}=80-4=76$ according to conditions (31), then the explicit
pattern of $b$ can be determined as Fig. 4. The denoising result is shown in
Fig. 3(a), and the transients can be easily identified with an almost pure
zero baseline. In this example, we use OGS and POGS with the atan penalty
function and non-convexity parameter $a$ set to its maximum value of
$1/(K_{1}\lambda)$, so as to maximally induce sparsity subject to the
constraint that the objective function $P_{1}$ is convex.

Fig. 3(b) shows another example using $N_{1}=2$ to determine the pattern $b$.
The result is slightly worse in RMSE than Fig. 3(a), because $N_{1}=2$ in this
example does not match the simulated data. However, $N_{1}=2$ is the lower
limit of a realistic value, because in practice it is very rare that the
transients are all single-sample spikes, when the data is properly measured.
As a consequence, the result in Fig. 3(b) can be understand as the worst case
of choosing an inappropriate $b$.

Fig. 3(c) shows the result of denoising using conventional non-periodic OGS
(3) with group size $K=8$ and the arctangent (atan) non-convex penalty
function. The regularization parameter $\lambda$ is chosen by the look-up
table in [9], which sets $\lambda$ proportional to the noise $\sigma$. In this
example, we set $\lambda=0.52\sigma$. The result in Fig. 3(c) misses some
faults and yields several false detections, e.g., the ones at about $0.82$
second, and some false transients appear before $t=0.36$ second.

As a comparison, we adopt conventional wavelet-based denoising method to the
test signal. More specifically, a 5-scale undecimated wavelet transform [20]
with 3 vanishing moments is used, and the result is shown in Fig. 3(d). For
denoising, we apply hard-thresholding and chose the threshold value by
$3\sigma$-rule for each subband. In Fig. 3(d), although some large amplitude
transients can be recovered at correct locations, they exhibit the same shape
as the chosen wavelet [see the box in Fig. 3(d)].

Figure 5:  Example 1: Results of fast spectral kurtosis. (a) Kurtogram. (b)
Envelope of extracted transients.

We compare the performance with fast spectral kurtosis [21] This method
produces the kurtogram in Fig. 5(a), where the kurtosis maximum is at the
third level with an estimated ‘optimal carrier frequency’ at 1000 Hz [see the
bright yellow area in Fig. 5(a)]. The corresponding amplitude of the extracted
transients is shown in Fig. 5(b), with an automatically generated threshold
shown as a gray dashed line. The peaks after $t=0.36$ seconds have a greater
density, which indicates that it is more likely that faults occur after
$t=0.36$ seconds. However, the useful repetitive transients are surrounded by
strong irrelevant noise.

###  4.1 Parameter selection

Figure 6:  Example 1: Optimal $\lambda$ at different noise level. (Note that
the vertical axes of the sub-figures are different.)

Setting pattern $b$. In many cases of fault diagnosis, it is feasible to
estimate the period of the transients based by component geometry and rolling
speed, but it might be difficult to estimate the duration of the transient.
However, since the binary structure $b$ is acting as a sliding window
capturing the global periodicity structure, it is not necessary to match the
length of a transient exactly. As a consequence, we suggest to constrain the
value of $N_{1}$ relatively small, e.g. $2\leqslant N_{1}\leqslant 4$.
However, if the sampling rate is quite high and a transient may contain more
samples, then the value of $N_{1}$ can be specified greater than 4. The value
of $M$ determines the number of non-zero entries (the value of $K_{1}$) in
$b$, and $K_{1}$ effects the non-convexity of the regularizer in the proposed
problem (14). In our experiments, we keep use $M=4$.

Setting regularization parameter $\lambda$. In order to explore the
correlation of the regularization parameter $\lambda$ to the binary pattern
$b$ (32), we show the optimal $\lambda$ as a function of the standard
deviation $\sigma$ of noise using different binary patterns in Fig. 6. We
define the optimal $\lambda$ as the value minimizing the RMSE for each fixed
binary pattern $b$.

In this test, using the data in Fig. 2(a), we search for $\lambda$ minimizing
an average RMSE (20 trials) at each noise level generated by different random
seeds. We present the results of select $M$ from 1 to 4 in Fig. 6(a) to (d),
respectively. In each figure, we show the optimal $\lambda$ as a function of
$\sigma$ under different $N_{1}$. Note that, in Fig. 6(a), when $M=1$ and
$N_{1}=1$, the proposed method adheres neither grouping nor periodic
structure, in which case, the problem (14) coincides to the BPD problem with
non-convex regularizer.

Note that since we simulate the test signal with a unit standard deviation,
the horizontal axis in Fig. 6 is also the ratio of deviation from noise to
data. In other word, we can use Fig. 6 as a look-up table, when the input
noise-to-signal ratio (SNR) is known.

Moreover, all the 16 curves in Fig. 6 show that the optimal $\lambda$ varies
approximately linearly with noise level. In practice, we suggest to chose
$\lambda$ in (14) proportional to the noise level as $\lambda=r\sigma$.
Through further experiments, we provide Table 3 as a guide for choosing the
multiplier $r$.

Table 3: Selection of $r$ for setting $\lambda$. $M$ $N_{1}$ | 1 | 2 | 3 | 4
---|---|---|---|---
1 | 3.700 | 1.700 | 1.150 | 0.925
2 | 1.700 | 0.850 | 0.625 | 0.475
3 | 1.150 | 0.625 | 0.450 | 0.375
4 | 0.925 | 0.475 | 0.375 | 0.325

###  4.2 Receiver operating characteristic (ROC) evaluation

Figure 7:  Example 1: ROC curves of faults detection. Figure 8:  Example of
classification rule.

The above comparative evaluation of OGS, POGS, wavelet, and fast spectral
kurtosis uses RMSE as an indicator of denoising performance. However, the RMSE
by itself is not a sufficient indicator of fault detection accuracy. In the
following, we use a receiver operating characteristic (ROC) based approach to
evaluate the relative accuracy of the methods.

The receiver operating characteristic (ROC) curve is a well known detection
performance evaluation methodology [22]. ROC curves are well suited to the
problem of vibration-based diagnosis applications [23]. An ROC curve is
generated by plotting the probability of a false alarm against the probability
of detection as the threshold level is varied. Since the POGS approach is
focused on machinery fault detection, the ROC curve is utilized to validate
the superior detection performance of POGS compared to other methods. We
define the classification rule as: if one sample in a transient [generated by
(33)] is detected as positive, then this entire transient a fault feature is
detected and all the remaining samples are all assigned to be positive.

This rule is slightly different to the one used in [23], since using the
sample-wise decision rule in [23] may cause problem of overweighting sample
recovery, but neglecting detecting fault features as transients of a cluster
of samples. We use a simple example in Fig. 8 to illustrate the problem and
our classification rule. Suppose that the true signal has 8 samples consists 2
periods (4 samples for each period), and each period has 2 samples positive
(labeled as circles). If Method 1 detects 2 samples belonging to two different
transients respectively, and method 2 detects 2 samples belonging to only one
transients, then the rule of [23] will give that the two methods have
identical accuracy, because if merely counting the samples, they have a same
number of samples detected. This is undesirable, since Method 2 misses an
entire transient. To overcome this issue, we re-label the detect result, where
if any sample in a transient of each period is detected, we re-label all the
samples in the entire transient to be detected. In the example of Fig. 8,
after the re-labeling, Method 1 has a better accuracy because it detects both
of the transients.

Fig. 7 shows the ROC curves from using OGS, POGS and wavelet-based method
respectively. POGS achieves an almost perfect detection result. Also, OGS is
better than wavelet-based method. In this example, because the results from
POGS with different parameter settings [in Fig. 3(a) and Fig. 3(b)] obtain
almost identical ROC curves, we chose to show the one in Fig. 3(a) where
$N_{1}=4$ and $M=4$.

We also use the extracted envelope from fast spectral kurtosis method to
perform a similar ROC analysis. Note that although we plot it together with
other ROC curves in Fig. 7, since its envelope has a different length than the
other results, the ROC curve is generated by a different number of samples.

##  5 Experimental and engineering data validation

###  5.1 Example 2: Rolling bearing with defect on outer race

Table 4: The parameters of the tested bearing Inner Race | Outer Race | Roller | Number of roller | Contact angle
---|---|---|---|---
160 mm | 290 mm | 34 mm | 17 | $0^{\degree}$
Figure 9: Example 2: Fault on the outer race of the testing bearing.

Figure 10: Example 2: (a) Test data. (b) Result of proposed method with
$N_{1}=2,M=4$. (c) Result of proposed method with $N_{1}=4,M=4$. (d) Result of
AR-MED.

In this example, the proposed approach is applied to a vibration signal
collected from a locomotive rolling element bearing with defect. The testing
locomotive rolling bearing with a slight scrape on its outer race is shown in
Fig. 9. The vibration signal is measured from acceleration sensors, using SONY
Ex data acquisition system when the electric locomotive was running. The
bearing type is 552732QT and its specification is shown in Table 4. The
sampling rate is 12.8 kHz and the current rotating speed is approximately 481
rpm. Thus, based on the geometric parameters and rotational speed, the
characteristic frequency of the outer race defect is calculated to be
$f_{o}=57.8$ Hz.

To adopt the proposed POGS (14), if the duration of a transient is uncertain,
we can directly chose the pattern $b$ (32) with $M=4$ and $N_{1}=2$, then by
Table 3 we can chose $\lambda$ with the respect to the noise level. Note that
since the ‘noise’ in the vibration data for fault detection is the background
out of the transient sequence, in practice it can be easily estimated using
healthy data.

In this example, we illustrate how to determine $\lambda$ without healthy
data. We estimate the deviation of background by the formula

$\displaystyle\hat{\sigma}=\mathsf{MAD}(y)/0.6745$ (34)

which is a conventional estimator of noise level used for wavelet-based
denoising [24], where $\mathsf{MAD}$ in (34) stands for median absolute
deviation, defined as

$\displaystyle\mathsf{MAD}(y):=\mathsf{median}(\left|y_{n}-\mathsf{median}(y)\right|).$
(35)

In this example, the estimated deviation of background is
$\hat{\sigma}=0.1606$ obtained by (34). The parameter $\lambda$ can be
obtained using $\lambda=r\hat{\sigma}$ where $r=0.475$ (from Table 3). The
result is shown in Fig. 10(b), where a periodical phenomenon can be easily
identified. Moreover, we also test POGS with $M=4$ and $N_{1}=4$, using the
same method to set $\lambda$. The result is shown in Fig. 10(c).

As a comparison, we also test the data with autoregression model assisted
minimum entropy deconvolution (AR-MED) [25]222 Implementation available online
http://www.mathworks.com/matlabcentral/fileexchange/41614-ar-filter-+-minimum-
entropy-deconvolution-for-bearing-fault-diagnosis. The order of AR filter is
chosen by maximizing Kurtosis value as the implementation suggested, and the
estimated filter in the MED step has a length of 50. The result is shown in
Fig. 10(d), where most of the impulses are promoted after deconvolution.
However, since the baseline is still relatively noisy, the regularity of
periodicity is not very clear.

###  5.2 Example 3: Motor bearing with multiple faults

Figure 11: Example 3: Spectra Quest’s machinery fault simulator. Figure 12:
Example 3: Test data. (a) The measured vibration signal of the motor housing.
(b) Fourier spectrum. (c) Hilbert envelope spectrum.  Figure 13: Example 3:
Outer race faults detection ($N_{1}=2,M=4$): (a) Extracted transients, (b)
Fourier spectrum and (c) Hilbert envelope spectrum. Figure 14: Example 3:
Inner race faults detection ($N_{1}=2,M=4$): (a) Extracted transients, (b)
Fourier spectrum and (c) Hilbert envelope spectrum.

Figure 15:  Example 3: Results of fast spectral kurtosis. (a) Kurtogram. (b) Envelope of extracted transients. Table 5: Fault frequencies for MFS motor bearing  Component | FTF | BPFO | BPFI | BSF
---|---|---|---|---
Motor Bearings | 0.384 | 3.066 | 4.932 | 2.03

To further demonstrate the effectiveness of the proposed method for machinery
fault detection, applications of motor bearing fault diagnosis are studied in
this section. The experiment is performed on a Spectra Quest’s machinery fault
simulator (MFS) illustrated in Fig. 11. The user does not need to make any
modifications to the motor provided. The simulation setup consists of a motor
with intentionally faulted bearings: one bearing with an inner race fault, and
one bearing with an outer race fault. Therefore, the fault diagnosis of the
motor bearings is equivalent to a compound faults detection case. The motor
bearing fault frequencies for MFS components are given in Table 5.
Accelerometers were mounted on the motor housing. The vibration signals were
measured at a sampling rate $f_{s}=6400$ Hz. The rotating speed of the motor
is 1433 rpm (23.89 Hz). Hence, the characteristic frequencies of the inner
race and outer race of the motor bearings are calculated to be $f_{i}\approx
117.8$ Hz and $f_{o}\approx 73.2$ Hz, respectively.

An observed vibration signal with a duration of 1 second is illustrated in
Fig. 12(a). However, the useful periodic pulses are buried in strong
background noise and irrelevant interference. The frequency spectrum and the
Hilbert envelope spectrum of the signal are shown in Fig. 12(b) and (c),
respectively. It can be observed from Fig. 12(b) that the energy of the
spectrum is distributed along the whole frequency range. Peaks at 24 Hz and
its harmonics can be observed in the low frequency band, which correspond to
the rotating frequency and its harmonics. The useful characteristic
frequencies used to monitor the health status of motor bearings can not be
observed in Fig. 12(c).

The proposed POGS approach is utilized to process the vibration signal. The
results are shown in Fig. 13 and Fig. 14. We run POGS twice with two group
structures, determined by the inner race period $T_{i}$ and the outer race
period $T_{o}$, for the purpose of separating the useful impulsive fault
features. In particular, keeping $N_{1}=2$ and $M=4$, we define two different
group structures using (31) and (32) based on the inner race period
$T_{i}=1/f_{i}$ and the outer race period $T_{o}=1/f_{o}$, respectively. We
use the atan penalty function with $a=1/(K_{1}\lambda)$ to ensure convexity of
the objective function. The value of $\lambda$ obtained using healthy data, is
set so as to diminish the healthy data to almost all zeros.

The two periodic-related values correspond to the outer race defect and the
inner race defect frequencies respectively. Strong periodic impulses with
intervals of approximately 0.0133 second (75 Hz) are clearly revealed Fig.
13(a), which is exactly in accordance with the outer race characteristic
frequency of 73.2 Hz. Similarly, periodic transient features with the period
0.0085 second (118 Hz) can be observed in Fig. 14(a), which is approximately
the inner race characteristic frequency of 117.8 Hz. To further reveal the
characteristic frequencies, the frequency spectrum and the Hilbert envelope
spectrum of the processed signals are shown in Fig. 13(b) and (c) and Fig.
14(b) and (c). We also present the smoothed profiles of the Hilbert envelop
spectrum to indicate the characteristic frequencies more clearly. The Hilbert
spectrum of the processed result illustrated in Fig. 13(c) is obtained using
the proposed approach with prior knowledge of the outer race characteristic
frequency 73.2 Hz. Apparently, the characteristic frequencies of outer race 73
Hz and its harmonic components are clearly revealed, as shown in Fig. 13(c).
Similarly, Fig. 14(c) is obtained utilizing the proposed approach with prior
knowledge of the outer race characteristic frequency 117.8 Hz. The
characteristic frequencies of inner race 117 Hz and its harmonic components
can be clearly observed in Fig. 14(c). Thus, the proposed periodic non-convex
regularized OGS approach successfully detects the compound faults of the motor
bearings. More specifically, the fault features of outer race defect and inner
race defect are clearly separated utilizing the proposed approach.

To further demonstrate the effectiveness of the proposed approach, we also
processed the vibration signal using spectral kurtosis and the results are
presented in Fig. 15. Fig. 15(a) is the fast kurtogram, where the optimal
carrier is detected at $2050$ Hz. Under this frequency, periodic transients
can be observed in the envelope of the filtered signal, in accordance with the
rotating speed of motor (23.89 Hz), as shown in Fig.15(b). No further fault-
related information can be observed in Fig. 15(b).

##  6 Conclusion

This paper proposes a periodic group sparsity approach for the purpose of
detecting faults in rotating machinery. The approach uses non-convex penalty
functions to promote periodic group sparsity. We show how to constrain the
non-convex penalty functions to ensure that the objective function is convex.
The OGS method was introduced in [9] and extended to non-convex regularized
OGS in [10]. A novelty of the proposed approach is that the proposed penalty
function models the periodicity of the sparse groups, making it suitable
specifically for feature extraction in machinery fault diagnosis. The period
of the sparse pulses is chosen based on prior knowledge of machine geometry
under inspect. Moreover, the proposed approach is able to separate compound
fault features by utilizing different periods of the periodic pulses
corresponding to different fault frequencies (e.g., outer race and inner race
characteristic frequencies of rolling element bearings) as demonstrated in
Section 5. The effectiveness of the proposed approach is verified by
simulation and experimental data. The processed results demonstrate that the
proposed approach outperforms other methods.

## Appendices

##  Appendix A Proof of Proposition 1

Proof.  This proposition can be proven by taking the second-order derivative,
where $\phi(v;a)$ is twice differentiable when $v\neq 0$. The second-order
derivative of (12) is

$\displaystyle
p^{\prime\prime}(v)=\gamma+\lambda\phi^{\prime\prime}(v;a),\quad v\neq 0.$
(36)

Therefore, when $v\neq 0$, it is sufficient that, if
$\gamma+\lambda\phi^{\prime\prime}(v;a)>0$, then $p^{\prime\prime}(v)>0$. When
$v=0$, Lemma A in [10] can be utilized directly, to show that, since
$p^{\prime}(0^{-})<p^{\prime}(0^{+})$, under the condition (13), function $p$
in (12) is strictly convex on .

##  Appendix B Proof of Proposition 2

Proof.  We rewrite $\sum_{k}b_{k}u_{k}^{2}$ with the respect to (10) as

$\displaystyle\sum_{k\in\mathcal{K}}b_{k}u_{k}^{2}=\sum_{k\in\mathcal{K}_{1}}b_{k}u_{k}^{2}+\sum_{k\in\mathcal{K}_{0}}b_{k}u_{k}^{2}.$
(37)

Then (18) is given by

$\displaystyle P(u)$
$\displaystyle=\frac{1}{2K_{1}}\Big{[}\sum_{k\in\mathcal{K}_{1}}b_{k}u_{k}^{2}+\underbrace{\sum_{k\in\mathcal{K}_{0}}b_{k}u_{k}^{2}}_{0}\Big{]}+\lambda\phi\bigg{(}\Big{[}\sum_{k\in\mathcal{K}_{1}}b_{k}u_{k}^{2}+\underbrace{\sum_{k\in\mathcal{K}_{0}}b_{k}u_{k}^{2}}_{0}\Big{]}^{1/2};a\bigg{)}$
(38a)
$\displaystyle=\frac{1}{2K_{1}}\sum_{k\in\mathcal{K}_{1}}u_{k}^{2}+\lambda\phi\bigg{(}\Big{[}\sum_{k\in\mathcal{K}_{1}}u_{k}^{2}\Big{]}^{1/2};a\bigg{)}.$
(38b)

Using Proposition 1 with $v=[\sum_{k\in\mathcal{K}_{1}}u_{k}^{2}]^{1/2}$ and
$\gamma=1/K_{1}$, it follows that, if (19) is satisfied, then $P$ is strictly
convex.

## References

  * [1] A. Heng, S. Zhang, A. C. Tan, J. Mathew, Rotating machinery prognostics: State of the art, challenges and opportunities, Mech. Syst. Signal Process. 23 (3) (2009) 724–739.
  * [2] R. Yan, R. X. Gao, X. Chen, Wavelets for fault diagnosis of rotary machines: A review with applications, Signal Process. 96 (2014) 1–15.
  * [3] R. B. Randall, J. Antoni, Rolling element bearing diagnostics – a tutorial, Mech. Syst. Signal Process. 25 (2) (2011) 485–520.
  * [4] F. Cong, J. Chen, G. Dong, F. Zhao, Short-time matrix series based singular value decomposition for rolling bearing fault diagnosis, Mech. Syst. Signal Process. 34 (1) (2013) 218–230.
  * [5] Y. Lei, J. Lin, Z. He, M. J. Zuo, A review on empirical mode decomposition in fault diagnosis of rotating machinery, Mech. Syst. Signal Process. 35 (1) (2013) 108–126.
  * [6] C. Zhang, B. Li, B. Chen, H. Cao, Y. Zi, Z. He, Periodic impulsive fault feature extraction of rotating machinery using dual-tree rational dilation complex wavelet transform, ASME J. Manuf. Sci. Eng. 136 (5) (2014) 051011.
  * [7] R. Yan, R. X. Gao, Harmonic wavelet-based data filtering for enhanced machine defect identification, J. Sound Vib. 329 (15) (2010) 3203–3217.
  * [8] W. He, Y. Zi, B. Chen, F. Wu, Z. He, Automatic fault feature extraction of mechanical anomaly on induction motor bearing using ensemble super-wavelet transform, Mech. Syst. Signal Process. 54-55 (2015) 457–480.
  * [9] P.-Y. Chen, I. W. Selesnick, Translation-invariant shrinkage/thresholding of group sparse signals, Signal Processing 94 (2014) 476–489.
  * [10] P.-Y. Chen, I. W. Selesnick, Group-sparse signal denoising: Non-convex regularization, convex optimization, IEEE Trans. Signal Process. 62 (13) (2014) 3464–3478.
  * [11] Q. He, X. Wang, Time–frequency manifold correlation matching for periodic fault identification in rotating machines, J. Sound Vib. 332 (10) (2013) 2611–2626.
  * [12] B. Chen, Z. Zhang, C. Sun, B. Li, Y. Zi, Z. He, Fault feature extraction of gearbox by using overcomplete rational dilation discrete wavelet transform on signals measured from vibration sensors, Mech. Syst. Signal Process. 33 (2012) 275–298.
  * [13] Z. Li, Z. He, Y. Zi, Y. Wang, Customized wavelet denoising using intra-and inter-scale dependency for bearing fault detection, J. Sound Vib. 313 (1) (2008) 342–359.
  * [14] W. He, Y. Zi, B. Chen, S. Wang, Z. He, Tunable Q-factor wavelet transform denoising with neighboring coefficients and its application to rotating machinery fault diagnosis, Sci. China Technol. Sci. 56 (8) (2013) 1956–1965.
  * [15] I. Selesnick, I. Bayram, Sparse signal estimation by maximally sparse convex optimization, IEEE Trans. Signal Process. 62 (5) (2014) 1078–1092.
  * [16] K. Lange, Optimization, Springer New York, 2004.
  * [17] S. Boyd, L. Vandenberghe, Convex optimization, Cambridge university press, 2009\.
  * [18] I. W. Selesnick, A. Parekh, I. Bayram, Convex 1-d total variation denoising with non-convex regularization, IEEE Signal Processing Letters 22 (2) (2015) 141–144.
  * [19] M. Figueiredo, J. Bioucas-Dias, R. Nowak, Majorization-minimization algorithms for wavelet-based image restoration, IEEE Trans. Image Process. 16 (12) (2007) 2980–2991.
  * [20] R. R. Coifman, D. L. Donoho, Translation-invariant de-noising, in: Wavelet and statistics, Springer-Verlag, 1995, pp. 125–150.
  * [21] J. Antoni, R. Randall, The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines, Mech. Syst. Signal Process. 20 (2) (2006) 308–331.
  * [22] S. Lee, P. White, Higher-order time–frequency analysis and its application to fault detection in rotating machinery, Mech. Syst. Signal Process. 11 (4) (1997) 637–650.
  * [23] J. Nichols, S. Trickey, M. Seaver, S. Motley, Using ROC curves to assess the efficacy of several detectors of damage-induced nonlinearities in a bolted composite structure, Mech. Syst. Signal Process. 22 (7) (2008) 1610–1622.
  * [24] D. Donoho, I. Johnstone, I. M. Johnstone, Ideal spatial adaptation by wavelet shrinkage, Biometrika 81 (1993) 425–455.
  * [25] H. Endo, R. Randall, Enhancement of autoregressive model based gear tooth fault detection technique by the use of minimum entropy deconvolution filter, Mechanical Systems and Signal Processing 21 (2) (2007) 906–919.