program

# Detection of Faults in Rotating Machinery Using
Periodic Time-Frequency Sparsity111pre-print submitted to Journal of Sound and
Vibration

Yin Ding email<EMAIL_ADDRESS>Tandon School of Engineering, New York
University, 6 Metrotech Center, Brooklyn, NY 11201, USA  Wangpeng He  School
of Aerospace Science and Technology, Xidian University, Xi’an, China  Binqiang
Chen  School of Aeronautics and Astronautics, Xiamen University, Xiamen, China
Yanyang Zi  School of Mechanical Engineering, Xi’an Jiaotong University,
Xi’an, China  Ivan W. Selesnick  Tandon School of Engineering, New York
University, 6 Metrotech Center, Brooklyn, NY 11201, USA

###### Abstract

This paper addresses the problem of extracting periodic oscillatory features
in vibration signals for detecting faults in rotating machinery. To extract
the feature, we propose an approach in the short-time Fourier transform (STFT)
domain where the periodic oscillatory feature manifests itself as a relatively
sparse grid. To estimate the sparse grid, we formulate an optimization problem
using customized binary weights in the regularizer, where the weights are
formulated to promote periodicity. In order to solve the proposed optimization
problem, we develop an algorithm called augmented Lagrangian majorization-
minimization algorithm, which combines the split augmented Lagrangian
shrinkage algorithm (SALSA) with majorization-minimization (MM), and is
guaranteed to converge for both convex and non-convex formulation. As
examples, the proposed approach is applied to simulated data, and used as a
tool for diagnosing faults in bearings and gearboxes for real data, and
compared to some state-of-the-art methods. The results show the proposed
approach can effectively detect and extract the periodical oscillatory
features.

##  1 Introduction

Condition monitoring and fault diagnosis of rotating machines are of great
importance to prevent machinery breakdown in numerous industrial applications.
Vibration signals generated by faults in rotating components have been widely
studied. Detecting and extracting transient vibration signatures is of vital
importance for vibration-based detection of faults in rotating machinery.
However, the observed vibration signals are usually corrupted by very heavy
background noise [59, 31]. Many diagnostic techniques have been developed to
extract the fault features based on different transforms, such as short time
Fourier transform (STFT) [80, 30, 43, 34, 58] and wavelet transform [79, 13,
78, 42, 48, 35, 38, 11, 41]. Some methods use morphological decomposition
methods, such as empirical mode decomposition (EMD) [46, 33], or demodulation
[68, 49, 75]. Some methods combine conventional time-frequency analysis with
other techniques suitable for transients detection, such as bispectrum-based
method [28, 19] and spectral kurtosis (SK) based methods [7, 6, 85, 12].

The use of sparse representations for fault detection was initially
illustrated in Ref. [80], where basis pursuit denoising (BPD) [17, 18] was
adopted to exploit the sparsity of vibration signals in different domains in
order to detect and extract fault features. In Ref. [58], the author
considered a morphological component analysis (MCA) problem [29, 70, 71],
which is a sparsity-based decomposition approach, to improve the extraction of
fault features. In Ref. [22] matching pursuit (MP) [53] was used on a
customized dictionary to extract oscillatory fault features. In Ref. [37], a
sparse representation using the tunable Q-factor wavelet transform [64, 63] is
utilized to extract oscillatory fault features. Most recently, sparse coding
techniques are introduced for fault feature extraction [73].

The periodicity (also called fault characteristic or fundamental frequency in
some cases) of potential fault features is very important information for
fault diagnosis. In many cases, this information can be simply obtained using
the geometry of the specific components under steady rotating speed, or
directly obtained from the user operation manual. It is already required as
necessary information in many fault diagnosis methods [79, 49, 72, 74]. Some
fault diagnosis methods use periodicity information as a verification step
after feature extraction (e.g. [50, 58]), and some methods use it to determine
an optimal frequency band when extracting the oscillation frequency of the
fault feature [49, 72, 74]. Only a few very recent works use the temporal
periodic structure to extract fault features [54, 36].

The overlapping group sparsity (OGS) approach has been discussed in [16, 15],
where a multi-dimensional group structure has been modeled in various domains
and convex optimization problems have been formulated, to take advantage of
group behavior (e.g. coefficients in STFT domain) to recover signals from
noisy observations. Using OGS to extract fault features was introduced and
formulated as an convex optimization problem in [36], in which the fault
features are assumed to be a periodic sequence of pulses. In addition, it was
shown that the objective function can be convex even if some specified non-
convex penalty functions are used to promote sparsity. Moreover, when the
periodicity information is known, a binary-weighted group structure was used
in the regularization, capturing the period of the potential fault features.

In this work, we consider the problem wherein the fault features are comprised
of a sequence of oscillatory transients, and sparsity should be promoted in
the time-frequency domain. The discrete-time vibration observation $y$ is
modeled by

$\displaystyle y=x+w,$ (1)

where $x$ is an approximately-periodic sequence of oscillatory transients
where the period is determined by the characteristic frequency of potential
faults, and $w$ is additive Gaussian noise. Moreover, we assume that the
oscillation frequencies of the transient component $x$ are unknown but
relatively sparse, and the periodicity is over the time domain.

Under the signal model, we seek a solution to the optimization problem

$c^{\star}=\arg\min_{c}\Big{\\{}F(c)=\frac{1}{2}\lVert
y-Ac\rVert_{2}^{2}+\lambda\Phi(c)\Big{\\}},$ (2)

for the purpose of extracting the fault feature $x\in{}^{N}$, where $c$
denotes the coefficients corresponding to an overcomplete transform $A$ of
$x$, where $x=Ac$, and $\lambda>$ 0 is a regularization parameter and $\Phi$
is a sparsity-promoting penalty function (regularizer) in the domain of $c$.
This paper aims to extract the oscillatory fault features as periodically
structured groups of coefficients in the time-frequency domain. We set the
operator $A$ to be the normalized inverse STFT, and the function $\Phi$ is
customized to capture the periodicity in the STFT domain.

To solve a relatively complicated sparsity-based optimization problem, where
an iterative scheme has to be adopted, methods are usually based on two
algorithms. One is based on majorization-minimization (MM) [32, 45, 39], such
as the algorithms given in [32, 56, 25, 66, 67, 65], and further the method of
iterated soft-thresholding algorithm (ISTA) [8, 24] can be considered as a
special case of MM asa well (majorizing data fidelity term based on Lipschitz
constant). Another trend is based on alternating direction method of
multipliers (ADMM) [10, 27], and its extension such as split augmented
Lagrangian shrinkage algorithm (SALSA) [2], and some recent applications using
such algorithm are presented in [55, 63].

Some algorithms can be explained as a combination of ADMM and MM as a special
case. For instance, a specific function of Bregman distance was described in
[57], which is similar to a majorizer, utilized to derive an iterative method
for total variation problem [62]. Such concept has been used on total
variation problem with multiple extensions such as combinations with wavelet-
based denoising, compressive sensing [77, 82, 52]. Moreover, the Bregman
distance function in [57] was used to in other iterative algorithms for convex
problems as well, such as the methods proposed in [83, 84]. Some recent work
in the field of optimization also consider this problem. In [47], a framework
is given that instead of using the Bregman distance function in [57], but
using proper ‘indefinite proximal terms’ to majorize augmented Lagrangian, the
resulting majorized ADMM-style algorithm will converge for convex problem.
Moreover, in [14], an ADMM-style algorithm with a embedded criteria checking
step is proposed based on majorizing the corresponding augmented Lagrangian as
well. Additionally, in [23], a scheme of majoring a cross term with the ADMM-
style iteration is proposed for convex problem.

In this work, we propose an algorithm that combines ADMM and MM with a more
general and simpler formulation, based on block successive upper-bound
minimization (BSUM) algorithm [61]. Moreover, the algorithm allows non-convex
regularization to promote the resulting sparsity, and it is guaranteed to
converge to a local minimum even though the objective function is not convex.

###  1.1 Notation

In this paper, the elements of a vector $x\in{}^{N}$ are denoted as $x_{n}$ or
$[x]_{n}$, and the norms are defined as

$\displaystyle\lVert x\rVert_{1}$ $\displaystyle:=\sum_{n}\left|x_{n}\right|,$
(3a) $\displaystyle\lVert x\rVert_{2}^{2}$
$\displaystyle:=\sum_{n}\left|x_{n}\right|^{2}.$ (3b)

The norms of a complex vector $z\in\mathbb{C}^{N}$ is defined as

$\displaystyle\lVert z\rVert_{1}$ $\displaystyle:=\sum_{n}\left|z_{n}\right|,$
(4a) $\displaystyle\lVert z\rVert_{2}^{2}$
$\displaystyle:=\sum_{n}\left|z_{n}\right|^{2},$ (4b)

where the definition of absolute value of a complex number $a\in\mathbb{C}$ is
given by

$\displaystyle\left|a\right|:=\sqrt{[\mathsf{Re}(a)]^{2}+[\mathsf{Im}(a)]^{2}},$
(5)

where $\mathsf{Re}(\cdot)$ and $\mathsf{Im}(\cdot)$ are operators taking the
real and imaginary part of a complex number respectively. As a consequence,
$\lVert z\rVert_{2}^{2}$ can be written explicitly as

$\displaystyle\lVert z\rVert_{2}^{2}$
$\displaystyle=\sum_{n}[\mathsf{Re}([z]_{n})]^{2}+[\mathsf{Im}([z]_{n})]^{2}$
(6a) $\displaystyle=\lVert z_{r}\rVert_{2}^{2}+\lVert z_{i}\rVert_{2}^{2},$
(6b)

where $z=z_{r}+\mathrm{i}z_{i}$, and $z_{r},z_{i}\in{}^{N}$ are vectors of the
real and imaginary part of $z$, respectively.

For multi-dimensional matrices, e.g. $B\in{}^{M\times N}$, we denote its
elements as $[B]_{m,n}$. We use bold $\mathbf{0}$ and $\mathbf{1}$ to express
all $0$’s and $1$’s matrices. More specifically, $\mathbf{0}(K,N)$ is a matrix
of size $K\times N$, and all its elements are zero, and $\mathbf{1}(K,N)$ is a
matrix, the elements of which are all one.

###  1.2 Review of Majorization-minimization

The majorization-minimization (MM) method [32, 45, 39] simplifies a
complicated optimization problem into a sequence of easier ones. To solve a
convex optimization problem

$\displaystyle u^{\star}=\arg\min_{u}F(u),$ (7)

the MM method is described by iteration

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

where $G:{}^{N}\times{}^{N}\to\real$ is a continuously differentiable
function, satisfying

$\displaystyle G(u,v)\geq F(u),\quad\text{ for all }u,v,$ (9a) $\displaystyle
G(u,u)=F(u).$ (9b)

The detailed proof of convergence for convex problems has been given in Ref.
[45, Chapter 10].

###  1.3 Review of BSUM

Block successive upper-bound minimization algorithm (BSUM) is described in
[61]. Its original goal is to overcome the uniqueness requirement of block
coordinate descent (BCD) method (also known as Gauss-Seidel method) [9, 51].
BSUM is an iterative method to solve optimization problems, which allows the
objective function to be non-convex and/or non-smooth. When the objective
function is neither convex nor smooth, the algorithm will at least converge to
a stationary point under some weak assumptions [61, Theorem 2]. This algorithm
has been applied to transceiver design in MIMO multi-cellular communication
system, which has a non-convex formulation [60]. In order to facilitate our
proposed algorithm, as a short review, we only rewrite a special case of BSUM
using our notation.

The problem is defined as

$\displaystyle x^{\star}=\arg\min_{x}Q(x),$ (10)

where $Q:{}^{N}\to\real$ is continuous and has at least one stationary point.
In this case, we slightly relax the definition of ‘optimal solution’, whereas
if the objective function is not convex, we allow the ‘optimum’ to be merely a
local minimizer. Note that, when function $Q$ is strictly convex, the unique
minimizer is the global minimizer, therefore the relaxation of optimal does
not affect convex problems, where $x^{\star}$ is still the true global
minimizer.

We assume the unknown $x\in{}^{N}$ can be divided into three blocks:
$x_{0}\in{}^{N_{0}}$, $x_{1}\in{}^{N_{1}}$ and $x_{2}\in{}^{N_{2}}$ with
$N=N_{0}+N_{1}+N_{2}$. Then the objective function $Q$ in (10) can be
rewritten using the expression of the three blocks, where the problem is

$\displaystyle\\{x_{0}^{\star},x_{1}^{\star},x_{2}^{\star}\\}=$
$\displaystyle\arg\min_{x\in{}^{N}}Q(x_{0},x_{1},x_{2}),\text{ where
}x=\begin{bmatrix}x_{0}\\\ x_{1}\\\ x_{2}\end{bmatrix},$ (11)

and if we can find a continuous function $G:{}^{N}\times{}^{N}\to\real$, which
is the block-wise upper-bound of $Q$, defined as

$\displaystyle G((x_{0},x_{1},x_{2}),(v_{0},x_{1},x_{2}))\geq
Q(x_{0},x_{1},x_{2}),\quad\text{ for all }x\text{ and }v_{0}\in{}^{N_{0}},$
(12a) $\displaystyle G((x_{0},x_{1},x_{2}),(x_{0},v_{1},x_{2}))\geq
Q(x_{0},x_{1},x_{2}),\quad\text{ for all }x\text{ and }v_{1}\in{}^{N_{1}},$
(12b) $\displaystyle G((x_{0},x_{1},x_{2}),(x_{0},x_{1},v_{2}))\geq
Q(x_{0},x_{1},x_{2}),\quad\text{ for all }x\text{ and }v_{2}\in{}^{N_{2}},$
(12c) $\displaystyle
G((x_{0},x_{1},x_{2}),(x_{0},x_{1},x_{2}))=Q(x_{0},x_{1},x_{2}),\quad\text{
for all }x,$ (12d)

then problem (11) can be solved by BSUM by the algorithm illustrated in [61,
Figure 2]. For the special case of three blocks, we can write the algorithm
explicitly as the steps in Table 1.

The condition (12) implies all the assumptions listed in [61, Assumption 2],
which ensures that the algorithm in Table 1 converges, and the mathematical
proof of convergence can be found in Section IV of [61]. It is worth noting
that the convergence behavior of BSUM does not rely on the convexity of the
problem.

Table 1: Block successive upper-bound minimization (BSUM) algorithm with three
blocks. $\displaystyle\text{initial:}~{}x\in{}^{N}$ repeat:
$\displaystyle\quad
x_{0}=\arg\min_{u_{0}}G((u_{0},x_{1},x_{2}),(x_{0},x_{1},x_{2}))$
$\displaystyle\quad
x_{1}=\arg\min_{u_{1}}G((x_{0},u_{1},x_{2}),(x_{0},x_{1},x_{2}))$
$\displaystyle\quad
x_{2}=\arg\min_{u_{2}}G((x_{0},x_{1},u_{2}),(x_{0},x_{1},x_{2}))$ end
$\displaystyle\text{return: }x=\\{x_{0},x_{1},x_{2}\\}$

##  2 SALMA

Split augmented Lagrangian shrinkage algorithm (SALSA) [2] is a recently
developed method for solving convex optimization problems. As an efficient
$\ell_{1}$-norm regularized convex problem solver, it has been widely used to
solve various signal processing problems. SALSA adopts variable splitting
technique and the concept of alternating direction method of multipliers
(ADMM) [10, 27], achieving a very simply structured and efficient iterative
algorithm. In this section, we propose an algorithm formulated very similar to
SALSA. However, in contrast to SALSA, the proposed method converges when used
to solve non-convex problems. Because we use the concept of majorization-
minimization (MM) in the proposed algorithm, it is termed split augmented
Lagrangian majorization-minimization algorithm (SALMA), and can be stated as
the following theorem.

######  Theorem 1.

For optimization problem

$\displaystyle x^{\star}=$
$\displaystyle\arg\min_{x\in{}^{N}}F_{0}(x)+F_{1}(x),$ (14)

where function $F_{0}:{}^{N}\to\real$ and $F_{1}:{}^{N}\to\real$ are both
finite and continuously differentiable, if a function
$G_{1}:{}^{N}\times{}^{N}\to\real$ majorizes $F_{1}$, i.e., satisfying

$\displaystyle G_{1}(u,v)\geq F_{1}(u),\quad\text{ for all }u,v\in{}^{N},$
(15) $\displaystyle G_{1}(u,u)=F_{1}(u),$ (16)

can be found, and it enables the following iterative algorithm

$\displaystyle u^{(i+1)}$
$\displaystyle=\arg\min_{u}G_{1}(u,u^{(i)})+\frac{\mu}{2}\lVert
u-x-d\rVert_{2}^{2},$ (17a) $\displaystyle x^{(i+1)}$
$\displaystyle=\arg\min_{x}F_{0}(x)+\frac{\mu}{2}\lVert u-x-d\rVert_{2}^{2},$
(17b) $\displaystyle d^{(i+1)}$ $\displaystyle=d^{(i)}-(u-x)$ (17c)

having unique solution at each step, then the algorithm (17) will converge,
and converges to a local minimizer of problem (14).

The proof of the above theorem in detail is given in Appendix A. Note that the
algorithm (17) has two steps identical to SALSA, and the difference is in the
sub-problem (17a), where the cost function is based on a continuously
differentiable majorizer of $F_{1}$. In contrast to SALSA, the proposed method
will converge when the objective function is not convex, because it is
actually a special case of BSUM.

In other words, the proposed algorithm utilizes a Gauss-Seidel procedure to
solve a majorization-minimization (MM) problem. This leads us to term the
algorithm SALMA (split augmented Lagrangian majorization-minimization
algorithm), which combines SALSA and majorization-minimization methods.

##  3 Proposed method

###  3.1 Smoothed penalty function

In this work, we use smoothed penalty function
$\phi_{\epsilon}:\real\to{}_{+}$ to induce sparsity, which is a smoothed
version of $\phi:\real\to{}_{+}$ in Table 2. It is defined as:

$\displaystyle\phi_{\epsilon}(u;a):=\phi(\sqrt{u^{2}+\epsilon};a),\quad\epsilon>0,$
(18)

where $\phi$ is the non-smooth penalty function satisfies the following
properties:

  1. 1.

$\phi(u;a)$ is continuous on .

  2. 2.

$\phi(u;a)$ is twice continuously differentiable on
$\mathbb{R}\backslash\\{0\\}$.

  3. 3.

$\phi(u;a)$ is even symmetric: $\phi(u)=\phi(-u)$.

  4. 4.

$\phi(u;a)$ is increasing and concave on +.

  5. 5.

$\phi(u;a)=\left|x\right|$, when $a=0$.

Several typical examples of $\phi$ are given in Table 2. Figure 1(a-b) gives
two specific examples of $\phi$ in dashed line, and corresponding smoothed
versions with $\epsilon=0.01$ in solid line.

In this work, we use the smoothed penalty function $\phi_{\epsilon}$ instead
of $\phi$, because after the smoothing, the function is continuously
differentiable on $u\in\real$, and we can easily find the corresponding
majorizer,

$\displaystyle
g(u,v):=\frac{u^{2}}{2\psi(v)}-\underbrace{\bigg{(}\frac{v^{2}}{2\psi(v)}-\phi_{\epsilon}(v;a)\bigg{)}}_{\text{only
depends on $v$}}.$ (19)

where $\psi(v):=v/\phi_{\epsilon}^{\prime}(v;a)$, and $g(u,v)$ is quadratic in
$u$. In Ref. [25], the same smoothed penalty functions are used and a detailed
proof is given to show that when $\phi$ satisfies the above 5 conditions, the
majorizer (19) satisfies

$\displaystyle g(u,v)\geq\phi_{\epsilon}(u;a),\quad\text{ for all }u,v,$ (20a)
$\displaystyle g(u,u)=\phi_{\epsilon}(u;a).$ (20b)

Figure 1(c-d) illustrates the majorizer $g$ of two penalty functions when
$v=0.5$. Note that $g$ never goes to extrema (i.e. $\pm\infty$) if $u$ and $v$
are both finite. Table 2 provides four examples of smoothed penalty functions.
Among all these functions, only the first one (‘abs’) is convex, and the rest
are all non-convex. In (18), $\phi$ is not differentiable at $0$, but
$\phi_{\epsilon}$ is differentiable, and when $\epsilon\to 0$, function
$\phi_{\epsilon}\approx\phi$, but preserving the differentiability at the
origin. In practice, we recommend to set $\epsilon$ very small, e.g.
$10^{-8}$.

Table 2: Sparsity-promoting penalty functions.

Penalty | $\phi(u;a)$ | $\phi_{\epsilon}(u;a)$ | $\psi(u)=u/\phi_{\epsilon}^{\prime}(u;a)\quad$
---|---|---|---
abs ($a=0$) | $\left|u\right|$ | $\sqrt{u^{2}+\epsilon}$ | $\sqrt{u^{2}+\epsilon}$
log | $\displaystyle\frac{1}{a}\log(1+a|u|)$ | $\displaystyle\frac{1}{a}\log(1+a\sqrt{u^{2}+\epsilon})$ | $\displaystyle\sqrt{u^{2}+\epsilon}\left(1+a\sqrt{u^{2}+\epsilon}\right)$
rat | $\displaystyle\frac{\left|u\right|}{1+a\left|u\right|/2}$ | $\displaystyle\frac{\sqrt{u^{2}+\epsilon}}{1+a\sqrt{u^{2}+\epsilon}/2}$ | $\displaystyle\sqrt{u^{2}+\epsilon}\left(1+a\sqrt{u^{2}+\epsilon}/2\right)^{2}$
atan | $\displaystyle\frac{2}{a\sqrt{3}}\left(\tan^{-1}\left(\frac{1+2a|u|}{\sqrt{3}}\right)-\frac{\pi}{6}\right)$ | $\displaystyle\frac{2}{a\sqrt{3}}\left(\tan^{-1}\left(\frac{1+2a\sqrt{u^{2}+\epsilon}}{\sqrt{3}}\right)-\frac{\pi}{6}\right)$ | $\displaystyle\sqrt{u^{2}+\epsilon}\left(1+a\sqrt{u^{2}+\epsilon}+a^{2}(u^{2}+\epsilon)\right)$

Figure 1: (a) $\ell_{1}$-norm penalty function and its smoothed version, and
(b) arctangent penalty function and its smoothed version, and (c) Majorizer
(gray) of the smoothed penalty function in (a). (d) Majorizer (gray) of the
smoothed penalty function in (b). For all the figures we set
$\epsilon=0.01,v=0.5,a=1$.

###  3.2 Problem definition

To estimate a sequence of oscillatory transients modeled by (1), through an
optimization problem prototyped by (2), we further formulate the optimization
problem explicitly as

$\displaystyle c^{\star}=\arg\min_{c}\Bigg{\\{}P(c)=\frac{1}{2}\lVert
y-Ac\rVert_{2}^{2}+\lambda\sum_{m_{1},m_{2}}\phi_{\epsilon}\big{(}\theta(c,B,m_{1},m_{2});a\big{)}\Bigg{\\}}$
(21)

where $y\in{}^{N}$, $c\in\mathbb{C}^{M_{1}\times M_{2}}$ are STFT
coefficients, and the matrix $A:\mathbb{C}^{M_{1}\times
M_{2}}\to\mathbb{C}^{N}$ here denotes the inverse short-time Fourier transform
operator. When a local optimal of problem (21) is obtained, component $x$
(periodic sequence of oscillatory transients) in signal model (1) is given by

$\displaystyle\hat{x}=Ac^{\star}.$ (22)

In problem (21), we defined the regularizer by function
$\theta:\mathbb{C}^{M_{1}\times M_{2}}\times{}^{K_{1}\times
K_{2}}\times\mathbb{Z}\times\mathbb{Z}\to\real$,

$\displaystyle\theta(c,B,m_{1},m_{2})$
$\displaystyle=\bigg{[}\sum_{k_{1}=0}^{K_{1}-1}\sum_{k_{2}=0}^{K_{2}-1}[B]_{k_{1},k_{2}}[c]_{m_{1}+k_{1},m_{2}+k_{2}}^{2}\bigg{]}^{1/2},$
(23)

which is similar to a 2-D convolution. Note that $B\in{}^{K_{1}\times K_{2}}$
is a two-dimensional binary-valued mask describing the group structure, which
is to be determined beforehand by the time-frequency information of $x$ and
the STFT.

We use $\odot$ to denote element-wise multiplication, so the function $\theta$
can be written as a weighted norm of a segment of STFT coefficient $c$,

$\displaystyle\theta(c,B,m_{1},m_{2}):=\lVert B\odot
S(c,m_{1},m_{2})\rVert_{2}.$ (24)

Here $S(c,m_{1},m_{2})$ denotes a $K_{1}\times K_{2}$ segment of matrix
$c\in\mathbb{C}^{M_{1}\times M_{2}}$, where the first element (up and left) is
$[c]_{m_{1},m_{2}}$. Moreover, since $c$ is a complex-valued array, where
$c=c_{r}+\mathrm{i}c_{i}$, it follows that,

$\displaystyle\lVert B\odot S(c,m_{1},m_{2})\rVert_{2}^{2}=\lVert B\odot
S(c_{r},m_{1},m_{2})\rVert_{2}^{2}+\lVert B\odot
S(c_{i},m_{1},m_{2})\rVert_{2}^{2}$ (25)

where $S(c_{r},m_{1},m_{2}),S(c_{i},m_{1},m_{2})\in{}^{K_{1}\times K_{2}}$ are
the corresponding arrays of the real and imaginary parts of
$S(c,m_{1},m_{2})$. Then, in implementation, the real and imaginary parts of
the problem can be calculated separately as real values.

For problem (21), in order to capture the known periodicity of $x$ as modeled
in (1), we further define binary weighting vectors $b_{1}\in\\{0,1\\}^{K_{1}}$
and $b_{2}\in\\{0,1\\}^{K_{2}}$ in the following structure,

$\displaystyle b_{1}$
$\displaystyle=\mathbf{1}{(1,K_{1})}=[\underbrace{1,1,\ldots,1}_{K_{1}}]$
(26a) $\displaystyle b_{2}$ $\displaystyle=[\mathbf{1}(1,N_{1}),\
\mathbf{0}(1,N_{0}),\ \mathbf{1}(1,N_{1}),\ \mathbf{0}(1,N_{0}),\ \ldots,\
\mathbf{1}(1,N_{1})]$ $\displaystyle=[\ \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}}\ ],$ (26b)

where $b_{2}$ captures the feature periodicity with a periodic binary
structure, so that the binary block $B$ spanning $M$ periods can be written
explicitly as

$\displaystyle B$ $\displaystyle=b_{1}^{\mathsf{T}}\times b_{2}=[\
\underbrace{\mathbf{1}(K_{1},N_{1}),\ \mathbf{0}(K_{1},N_{0}),\
\mathbf{1}(K_{1},N_{1}),\ \mathbf{0}(K_{1},N_{0}),\ \ldots,\
\mathbf{1}(K_{1},N_{1})}_{\text{spanning }M\text{ periods}}\ ],$ (27)

which is concatenating of $M$ all-$1$’s blocks and $M-1$ all-$0$’s blocks
intervally. For instance, when $M=4$, $K_{1}=4$, $N_{0}=4$ and $N_{1}=2$, the
block $B$ is a binary grid illustrated in Figure 2.

Figure 2: Example of binary block $B$ when $M=4$, $K_{1}=4$, $N_{0}=4$ and
$N_{1}=2$.

From (27) we can derive that the two-dimensional binary-weight block
$B\in{}^{K_{1}\times K_{2}}$ has a size determined by $K_{1}$ and

$\displaystyle K_{2}$ $\displaystyle=(N_{0}+N_{1})\times(M-1)+N_{1}.$ (28)

Using $B$ in the objective function of (21), the regularization term groups a
2-D structure of STFT coefficients with the respect to the periodicity
information of $x$ in (1).

###  3.3 Algorithm derivation

To solve the optimization problem (21), we split the variable $c$ by
introducing an extra equality constraint, and solve the equivalent problem

$\displaystyle\\{c^{\star},u^{\star}\\}=$
$\displaystyle\arg\min_{c,u}P_{0}(c)+P_{1}(u)$ (29a) $\displaystyle\text{
subject to }c-u=0$ (29b)

where $P_{0}:\mathbb{C}^{M_{1}\times M_{2}}\to\real$ and
$P_{1}:\mathbb{C}^{M_{1}\times M_{2}}\to\real$ are defined as,

$\displaystyle P_{0}(c)$ $\displaystyle:=\frac{1}{2}\lVert y-Ac\rVert_{2}^{2}$
(30a) $\displaystyle P_{1}(u)$
$\displaystyle:=\lambda\sum_{m_{1}}\sum_{m_{2}}\phi_{\epsilon}\big{(}\theta(u,B,m_{1},m_{2});a\big{)},$
(30b)

which leads to a same solution by $c=u$ to problem (21).

The majorizer of $P_{1}$ in (29) can be found by $g(u,v)$ in (19) as,

$\displaystyle G_{1}(u,v)$
$\displaystyle=\lambda\sum_{m_{1}}\sum_{m_{2}}g\big{(}\theta(u,B,m_{1},m_{2}),\theta(v,B,m_{1},m_{2})\big{)}$
(31a) $\displaystyle=\lambda\sum_{m_{1}}\sum_{m_{2}}g\big{(}\lVert B\odot
S(u,m_{1},m_{2})\rVert_{2},\lVert B\odot S(v,m_{1},m_{2})\rVert_{2}\big{)}$
(31b)
$\displaystyle=\frac{\lambda}{2}\sum_{m_{1}}\sum_{m_{2}}[r(v)]_{m_{1},m_{2}}[u]_{m_{1},m_{2}}^{2}+C$
(31c)

where $C$ is a constant, and $r(v)\in{}^{M_{1}\times M_{2}}$ is a matrix
having an explicit form

$\displaystyle[r(v)]_{m_{1},m_{2}}=\sum_{k_{1}=0}^{K_{1}-1}\sum_{k_{2}=0}^{K_{2}-1}\frac{[B]_{k_{1},k_{2}}^{2}}{\psi(\lVert
B\odot S(v,m_{1}-k_{1},m_{2}-k_{2})\rVert_{2})}.$ (32)

Then using SALMA (Theorem 1), problem (21) can be solved iteratively by

$\displaystyle u^{(i+1)}$
$\displaystyle=\arg\min_{u}G_{1}(u,u^{(i)})+\frac{\mu}{2}\lVert
u-c-d\rVert_{2}^{2}$ (33a) $\displaystyle c^{(i+1)}$
$\displaystyle=\arg\min_{c}P_{0}(c)+\frac{\mu}{2}\lVert u-c-d\rVert_{2}^{2}$
(33b) $\displaystyle d^{(i+1)}$ $\displaystyle=d^{(i)}-(u-c).$ (33c)

####  3.3.1 Sub-problem (33a)

The step (33a) can be written as

$\displaystyle u$ $\displaystyle=\arg\min_{u}G_{1}(u,v)+\frac{\mu}{2}\lVert
u-c-d\rVert_{2}^{2}$ (34)

where $G_{1}(u,v)$ is given by (31), therefore the problem has a more explicit
formulation

$\displaystyle u$ $\displaystyle=\arg\min_{u}\frac{\mu}{2}\lVert
u-(c+d)\rVert_{2}^{2}+\frac{\lambda}{2}\sum_{m_{1}}\sum_{m_{2}}[r(v)]_{m_{1},m_{2}}[u]_{m_{1},m_{2}}^{2},$
(35)

which is a least square problem with explicit solution

$\displaystyle u$
$\displaystyle=(c+d)./\Big{[}1+\frac{\lambda[r(v)]}{\mu}\Big{]},$ (36)

where $./$ denotes the point-wise division.

####  3.3.2 Sub-problem (33b)

The sub-problem (33b) can be written as

$\displaystyle c$ $\displaystyle=\arg\min_{c}\frac{1}{2}\lVert
y-Ac\rVert_{2}^{2}+\frac{\mu}{2}\lVert u-c-d\rVert_{2}^{2}$ (37)

which is also a least-square problem, and it has an explicit solution

$\displaystyle c=\Big{[}A^{\mathsf{H}}A+\mu
I\Big{]}^{-1}\Big{[}A^{\mathsf{H}}y+\mu(u-d)\Big{]}.$ (38)

The operator $A^{\mathsf{H}}$ is the complex conjugate transpose (Hermitian)
of $A$. Note that to compute $c$ by (38) we need to solve an inverse system
$\Big{[}A^{\mathsf{H}}A+\mu I\Big{]}^{-1}$. To overcome the computational
cost, we use the matrix inverse lemma (see the formulation in Appendix B),

$\displaystyle\Big{[}A^{\mathsf{H}}A+\mu I\Big{]}^{-1}$
$\displaystyle=\frac{1}{\mu}\bigg{\\{}I-A^{\mathsf{H}}\Big{[}\mu
I+AA^{\mathsf{H}}\Big{]}^{-1}A\bigg{\\}}$ (39a)
$\displaystyle=\frac{1}{\mu}\Big{[}I-\frac{1}{\mu+1}A^{\mathsf{H}}A\Big{]},$
(39b)

where since the operators $A^{\mathsf{H}}$ and $A$ here are the forward and
inverse STFT, where $x=AA^{\mathsf{H}}x$, then the property

$\displaystyle AA^{\mathsf{H}}=I$ (40)

holds. As a consequence, equation (38) can be written as,

$\displaystyle c$
$\displaystyle=\frac{1}{\mu}\Big{[}I-\frac{1}{\mu+1}A^{\mathsf{H}}A\Big{]}\Big{[}A^{\mathsf{H}}y+\mu(u-d)\Big{]}$
(41a)
$\displaystyle=\frac{1}{\mu}A^{\mathsf{H}}y+(u-d)-\frac{1}{\mu(\mu+1)}A^{\mathsf{H}}y-\frac{1}{1+\mu}A^{\mathsf{H}}A(u-d)$
(41b)
$\displaystyle=(u-d)+\frac{1}{\mu+1}A^{\mathsf{H}}\Big{[}y-A(u-d)\Big{]},$
(41c)

which is equivalent to (38), but does not require solving a linear system in
(38).

Using the results derived in (36) and (41), which solve the problems (33a) and
(33b) respectively, the algorithm (33) has explicit solutions of each step,

$\displaystyle u^{(i+1)}$ $\displaystyle=(c+d)./\Big{[}1+\frac{\lambda
r(u^{(i)})}{\mu}\Big{]}$ (42a) $\displaystyle c^{(i+1)}$
$\displaystyle=(u-d)+\frac{1}{\mu+1}A^{\mathsf{H}}\Big{[}y-A(u-d)\Big{]}$
(42b) $\displaystyle d^{(i+1)}$ $\displaystyle=d^{(i)}-(u-c).$ (42c)

The entire algorithm to solve problem (21) is summarized in Table 3. Note that
although the variables $u,c,d\in\mathbb{C}^{M_{1}\times M_{2}}$ in the
algorithm are all 2-D complex-valued matrices, in fact this does not affect
the feasibility of SALMA. In each step of the algorithm (Table 3), the real
and imaginary parts can be computed parallelly and separately, using real-
valued expressions. More specifically, it is equivalent to implement the
algorithm using 3-D matrices on ${}^{M_{1}\times M_{2}\times 2}$, and each
step leads to a same result. Hence, SALMA (Theorem 1) defined in real domain,
is still valid to solve the problem (21) (see Appendix C for more details).

Table 3: Solving problem (21) with SALMA.
$\displaystyle\text{initial:}~{}c,u,d\in\mathbb{C}^{M_{1}\times M_{2}},\mu>0,$
repeat:
$\displaystyle\quad[r]_{m_{1},m_{2}}=\sum_{k_{1}}\sum_{k_{2}}\frac{[B]_{k_{1},k_{2}}^{2}}{\psi(\lVert
B\odot S(u,m_{1}-k_{1},m_{2}-k_{2})\rVert_{2})}$ $\displaystyle\quad
u=(c+d)./\Big{[}1+\frac{\lambda}{\mu}r\Big{]}$ $\displaystyle\quad
c=(u-d)+\frac{1}{\mu+1}A^{\mathsf{H}}\Big{[}y-A(u-d)\Big{]}$
$\displaystyle\quad d=d-(u-c)$ end $\displaystyle\text{return: }c,\ x=Ac.$

##  4 Simulated data example

###  4.1 Example 1

Figure 3: (a) Simulated signal, and (b) noisy observation.

Figure 4: Result from convex penalty function, (a) estimated time domain
signal, (b) STFT coefficients, (c) detected oscillatory frequency components,
(d) Hilbert envelope spectrum of estimated signal. Figure 5: Converging
history of using different $\mu$’s.

Figure 6: Denoising results of (a) bandpass filtering, (b) wavelet-based
denoising.

Figure 7: Filtering result using MCKD of (a) test data, (b) test data with a
reduced noise level.

Figure 8: Enhanced result using non-convex penalty function, (a) estimated
time domain signal, (b) STFT coefficients, (c) detected frequency components,
(d) Hilbert envelope spectrum of estimated signal.

We apply the proposed method to the simulated data illustrated in Figure 3(a).
The simulated data is a 0.25 second signal with sampling rate $f_{s}=16$ kHz.
We simulate the vibration features caused by fault as a transient containing
two oscillation components at $1$ kHz and $2$ kHz. Each transient lasts for
about $4$ milliseconds with a random amplitude within a range, and the two
oscillatory components in one transient have independently random initial
phases. The detail of a simulated oscillatory fault is shown in Figure 3(a)
(inside the dashed line box). The clean data shown in Figure 3(a) is a
sequence of the simulated faults with a period $T=10$ milliseconds. Note that
since each transient has a random phase, in Figure 3(a), each transient is
different than the others. We add heavy noise to the test data, and the noisy
signal is shown in Figure 3(b).

In this example, we set the STFT to have a window size $R=32$ with $50\%$
window overlapping, and the length of Fourier transform is $L=256$. We can
establish the block $B$ in problem (21) by a condition of $b_{2}$ in (26)
using the periodicity information, where the length of one period in the STFT
domain is

$\displaystyle\frac{2Tf_{s}}{R}\approx N_{0}+N_{1}.$ (44)

In this example, we set $K_{1}=N_{1}=2$, and $N_{0}+N_{1}=2Tf_{s}/R=10$, and
$b_{2}$ as a binary weight vector is spanning $M=4$ periods. Hence, by (28),
the size of binary weight block $B$ is $2\times 32$. When using (27) to
express the block, $B$ consists of four blocks with $2\times 2$ ones and three
blocks with $2\times 8$ zeros.

Figure 4 shows the results obtained using the proposed method with convex
(smoothed $\ell_{1}$-norm) penalty function. Figure 4(b) is the STFT
coefficients achieved by solving problem (21), and Figure 4(a) is the
estimated $\hat{x}$ calculated from the STFT coefficients in Figure 4(b). The
proposed method recovers the oscillatory transients at accurate time and
frequency. Although the amplitude of the estimated $\hat{x}$ is slightly
diminished, all the ‘fault features’ (transients) are recovered at accurate
locations, and in Figure 4(b), a grid of sparse blocks is distributed at the
correct frequencies (about $1$ and $2$ kHz). In this example, we set
$\lambda=18$ (the method to select the parameters will be discussed in the
following section). Note that, similar to SALSA, the proposed method requires
a parameter $\mu>0$ which effects the convergence rate. We use several values
of $\mu$ for this example, the comparison of convergence speed is shown in
Figure 5. Through numerical experiments, we find that setting $\mu=1$ gives a
reasonable convergence rate, not only in this example, but also for the
applications illustrated in the following sections. To identify the frequency
components, we sum the obtained STFT spectrum over time, to obtain a frequency
distribution in Figure 4(c), where two peaks indicate the true oscillatory
frequencies. Meanwhile, the Hilbert envelope spectrum of the estimated signal
is illustrate in Figure 4(d). The useful fault frequency (inter-transient
period) and its harmonic components are clearly revealed in the Hilbert
envelope spectrum.

For comparison, we utilize bandpass filtering. Bandpass filtering is a
conventional and effective method in fault diagnosis field to extract
oscillatory features or emphasize specific frequency components, and is used
in conjunction with some other techniques [40, 76, 68, 21]. For instance, in
[76], the author uses bandpass filtering after removing detected harmonics,
and in [68, 69], the author demodulates vibration signals at characteristic
frequency after using specially designed bandpass filters as pre-processing,
and in [21], bandpass filtering is used as a pre-processing for a singular
value decomposition (SVD) based analysis.

In this example, we use two third-order Butterworth bandpass filters with
center frequencies at $1$ kHz and $2$ kHz respectively with ‘forward-backward’
zero-phase filtering, and the passband is about $400$ Hz for both filters.
Because the the passbands of the two filters do not overlap, we simply add the
output signals from these two filters, and the result is shown in Figure 6(a).
This result shows that linear time-invariant (LTI) filtering technique hardly
recovers a zero baseline. Moreover, although some transients with large
amplitude can be visually distinguished in Figure 6(a), it is hard to
distinguish the transients with small amplitudes, e.g., at about $t=0.08$ to
$0.12$ second. Furthermore, the result in Figure 6(a) is obtained by knowing
the oscillatory frequencies exactly, but in practice, especially in the area
of fault diagnosis, it is rarely possible. In this example, even though we use
the accurate frequencies as prior-knowledge for the bandpass filters, the
proposed method still attains a significantly better root-mean-square error
(RMSE) value.

Wavelet-based denoising has also been used to extracted fault features [79,
81, 1, 48]. In this example, we adopt a conventional wavelet-coefficient
thresholding method to denoise the test signal. More specifically, a 5-scale
undecimated wavelet transform [20, 44] with 3 vanishing moments is used, and
the result is shown in Figure 6(b). For denoising, we apply hard-thresholding
and chose the threshold value by $3\sigma$-rule for each subband. In Figure
6(b), although some large amplitude transients can be recovered at correct
locations, they lose the oscillatory structure, where most recovered
transients have irregular waveforms [see the cropped example inside the dashed
line box in Figure 6(b)]. In the time range from about $0.10$ to $0.13$
seconds, the transients are almost completely diminished.

As an application of fault diagnosis, we apply maximum correlated Kurtosis
deconvolution (MCKD) [54] to the test data222 Matlab implementation is
available online:
http://www.mathworks.com/matlabcentral/fileexchange/31326-maximum-correlated-
kurtosis-deconvolution--mckd-. This method deconvolves the signal from an
estimated FIR filter that maximizes the correlated Kurtosis value locally. It
requires the period of potential fault features. To proceed with the
algorithm, we set the period length to be the true period (160 samples), and
the FIR filter length to be 100, and the M-shift value to be 5 as suggested in
the Matlab implementation available online. The result of MCKD is shown in
Figure 7(a), where we found that the algorithm is not able to give a
deconvolution result that properly indicates most of the fault features at the
noise level of our test signal. When comparing to the noisy data in Figure
3(b), we can see the algorithm mistakenly promotes the impulsive features in
noise. When we reduce the noise level, the algorithm is capable to properly
indicate periodic fault features as promoted impulses after deconvolution. In
Figure 7(b), we show an example of applying MCKD to the data in Figure 3(b),
but with a reduced noise level from $\sigma=150$ to $50$.

###  4.2 Enhancement by non-convex penalty

Non-convex penalty functions can be utilized to further enhance the
effectiveness of our proposed formulation. When the objective function (21)
has a non-convex regularizer, the problem is generally non-convex, and when
minimizing a non-convex function, the algorithm may get trapped in a local
minima. To reduce this problem, we use the solution from the convex
formulation to initialize the algorithm. Through numerical experiments, we
suggest a range of values for parameter $a$, which controls the ‘non-
convexity’ of the penalty functions in Table 2,

$\displaystyle 0\leq a\leq\frac{1}{\lambda K},\quad
K=\sum_{k_{1},k_{2}}[B]_{k_{1},k_{2}},$ (45)

where $K$ is also the number of $1$’s when we restrict $B$ to be a binary
block as (27). This range does not guarantee the problem is convex, but in
practice, using the non-convex penalty regularizers in this range
significantly enhances the result and gives a stable convergence behavior.

To further reduce the problem of local minima, we can gradually increase the
parameter $a$ to its maximum value. For instance, when the upper bound of $a$
is $0.25$ in (45), we can run the algorithm (in Table 3) starting with $a=0$
obtaining a result with convex penalty function, and then run the algorithm
with non-convex penalty function five more times initialized by previous
iteration, gradually increasing $a$ by $0.05$ each time. This simultaneously
helps to restrict the non-convex solution close to the convex solution and
promote stronger sparsity.

In this example, we use the arctangent function (‘atan’ in Table 2), and the
result is shown in Figure 8(a). There is a significant improvement reflected
by the RMSE value compared to the convex penalty function in Figure 4. We also
show the absolute value of STFT coefficients in Figure 8(b), and the frequency
component indication in Figure 8(c), where significant enhancement of sparsity
can be seen. Meanwhile, the Hilbert envelope spectrum of estimated signals is
illustrate in Figure 8(d). The useful fault frequency and its harmonic
components are clearly detected in the Hilbert envelope spectrum.

###  4.3 Parameter selection

Setting STFT parameters. Setting STFT parameters. To use the proposed method,
we suggest choosing a relatively small window length because, assuming the
length of the periodic transient is unknown, a short window length can
preserve the resolution of the time-frequency spectrum. Moreover, in order to
promote the sparsity effectively, the transform A should be over-complete. We
suggest choosing the length of the FFT for each window (the value L in the
examples) 4 or 8 times the window length (the value R in the examples). If $L$
is too large, it will affect the computation efficiency. In all the tests and
examples, we consistently used the STFT with a squared-sine window with 50 %
overlapping.

Setting binary weighting grid $B$. Presumably, other than the information of
periodicity, any other facts about the oscillatory feature is not given in the
proposed method. Consequently, as a lower limit of grouping scheme, we can
chose the grid on the time-frequency domain adhering $K_{1}=N_{1}=2$, so that
a neighboring coefficients on both time and frequency domain are considered.
Moreover, the value of $M$, which is the number spanning in $B$, has no
specific limits, and in all tests and examples, we consistently set the value
of $M$ to 4 or 8.

Setting regularization parameter $\lambda$. Based on signal model (1), we
assume the noise (excluding the periodic oscillatory transients) is white
noise. Consequently, when the observation $y$ is noise only as $w$, the
solution of problemn(21) should be almost all zero. In this case, if the
variance of noise $\sigma_{w}^{2}$ is known, a signal with only noise
$\hat{w}_{n}\sim\mathcal{N}(0,\sigma_{w})$ can be built, and the
regularization parameter $\lambda$ in (21) can be trained by letting
$x^{\star}=Ac^{\star}\approx\mathbf{0}$ through (21) when $y=\hat{w}$. It is
well-known that in most signal denoising methods using sparsity-based
optimization the regularization parameter is proportional to the noise level,
wherein the heavier the noise is, the larger regularization parameter should
be utilized. Through numerical experiments, we found that for problem (21),
the optimal lambda is an approximately linear function of $\sigma_{w}$.
Moreover, it is certain that when there is no noise in the observation signal,
$\lambda$ should be 0, therefore other than building a noise signal to select
$\lambda$ for (21), alternatively it is only necessary to determine a factor
$\eta$ and set $\lambda\approx\eta\sigma_{w}$. For guidance, we provide a
table as a reference to select $\lambda$ at the variance of noise $\sigma=1$,
under several typical settings of STFT parameters.

Table 4: Factor $\eta$ for selecting $\lambda=\eta\sigma$, when fixing $K_{1}=N_{1}=2,\sigma=1$. Parameters of STFT |  $\eta$ ($M=4$) |  $\eta$ ($M=8$)
---|---|---
$R=16,L=64$ | 0.120 | 0.060
$R=16,L=128$ | 0.085 | 0.060
$R=32,L=128$ | 0.120 | 0.065
$R=32,L=256$ | 0.090 | 0.060

Note that in practice when the deviation of noise $\sigma_{w}$ is not unity,
it is just proportional to the case of $\sigma=1$, so it is enough to provide
a suggestion of $\eta$ for $\lambda=\eta\sigma$ when $\sigma=1$ in Table 4.
Using the table, we can set $\lambda$ by

$\displaystyle\lambda=\eta\cdot\alpha=\eta\cdot\sigma_{w}.$ (46)

For instance, in the example shown in Figure 4 and Figure 8, since we know the
noise level is $\sigma_{w}=150$, then when the STFT has the parameters
$R=16,L=64$ and the binary grid has $M=4$, $\lambda$ can be selected by
$\lambda=0.120\times 150=18$.

Moreover, when the noise level is not known, we suggest to use a method given
in [26] to estimate the noise level, which is

$\displaystyle\hat{\sigma}_{w}=\mathsf{MAD}(y)/0.6745$ (47)

which is a conventional estimator of noise level dependent on the noisy
observation only, used for wavelet-based denoising, where $\mathsf{MAD}$ is
the median absolute deviation operator defined as

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

Note that all the examples in this paper use the above method to select
parameters, including the examples of synthetic data before this section and
the examples of experimental data (Example 2 and Example 3) following this
section, and all the specific parameters used in the examples are listed in
the figures.

##  5 Experiment and engineering validation

In this section, we illustrate the proposed method by applying it to data
recorded from two different machines. One is from a public dataset concerning
bearing faults. The other data was collected by State Key Laboratory for
Manufacturing and Systems Engineering of Xi’an Jiaotong University for fault
diagnosis of gearbox [13].

###  5.1 Example 2: experimental data validation

Table 5: Defect frequencies of drive end bearing (multiple by running speed) Inner race | Outer race | Cage train | Rolling element
---|---|---|---
5.4152 | 3.5848 | 0.39828 | 4.7135
Figure 9: Example 2: Test stand for vibration monitoring.

Figure 10: Example 2: (a) test data, (b) estimated signal, (c) STFT
coefficients, (d) smoothed profile.

Figure 11: Example 2: (a) test data of normal setup, (b) estimated signal

Figure 12: Example 2: (a) Kurtogram, (b) envelop.

In this example, we utilized the proposed method to detect potential inner
race fault of the drive end bearing from an experimental setup. The
experimental vibration data was acquired from the Case Western Reserve
University Bearing Data Center333The data is available online at:
http://csegroups.case.edu/bearingdatacenter/pages/download-data-file.

As shown in Figure 9, the test stand consists of a 2 hp motor (left), a torque
transducer/encoder (center), and a dynamometer (right). The bearing installed
on the drive end is a 6205-2RS-JEM SKF deep-groove ball bearing. The potential
defect frequencies of drive end bearing are listed in Table 5, dependent on
the location of the fault. More details about the experiment setup can be
found at the website
http://csegroups.case.edu/bearingdatacenter/pages/apparatus-procedures.

A measured vibration signal is shown in Figure 10(a), with a sampling rate
$f_{s}=12$ kHz. The signal was collected using accelerometers, which were
attached to the housing with magnetic bases. The current rotational speed of
the motor is approximately 1730 r/min. Thus, according to Table 5, the fault
characteristic frequency of the drive end bearing inner race is about
$f_{i}=1730/60\times 5.4152\approx 156$ Hz. To run the proposed algorithm, we
set the STFT window length to $R=16$ samples, and the number of the frequency
bins to $128$ for each window. For the binary block in the STFT domain, we set
that each block spans $M=4$ periods, and $K_{1}=2$, $N_{1}=2$, then the block
$B$ can be determined by (26), (27) and (44). Moreover, the regularization
parameter $\lambda$ is selected using Table 4 and formula (46) with an
estimated $\hat{\sigma}_{w}=0.225$ calculated by (47).

In order to validate the reliability of the proposed algorithm, we also test
it on a normal (no fault) setup under the same rotation speed. The test signal
is shown in Figure 11(a), and using the same parameters, there is no
oscillatory feature extracted, and the output of the proposed method is almost
zero [see Figure 11(b)].

To achieve a more visualizable oscillatory feature detection signal, we sum
the absolute values of STFT coefficients $c^{\star}$ through frequency domain,
and then the obtained result is lowpass filtered,

$\displaystyle
s(m_{2})=\textsf{LPF}\Bigg{\\{}\sum_{m_{1}}\left|[c^{\star}]_{m_{1},m_{2}}\right|\Bigg{\\}},$
(49)

where $m_{1}$ denotes the frequency index and LPF denotes a lowpass filter.
Then we can obtain a smooth time series $s\in{}^{M_{2}}$, indicating the time
and strength of extracted oscillatory features caused by faults [see Figure
10(c)]. For comparison, we utilize fast spectral kurtosis (SK) [7, 6] 444The
matlab implementation is available online:
http://www.mathworks.com/matlabcentral/fileexchange/48912-fast-kurtogram to
analyze the same vibration data. This method extracts the ‘envelope’ of the
optimal subband indicating potential fault features. The Kurtogram and the
envelope of the optimal subband is shown in Figure 12, where we use the
optimal subband SK detected automatically at $3$ kHz and level $2.2$. In
Figure 12(b) the large spikes indicate the fault features. The profile in
Figure 10(d) more clearly indicates the fault features than the one in Figure
12(b). We can apply lowpass filtering to the ‘noisy’ envelop extracted by SK
as well, but in this case, it tends to eliminate some spikes with small
amplitudes, such as the one at about $t=0.03$ second.

Figure 13: Example 3: Schematic of the oxygen separation and compression unit
and the sensor locations (left), and Illustration of rub-impact between the
end face of the anti-thrust plate on the pinion and that of the bull gear
(right).

Figure 14: Example 3: (a) measured test data. (b) frequency domain.

###  5.2 Example 3: rub-impact fault detection in the gearbox of an oil
refinery

In this engineering application, we applied the proposed method to detect a
rub-impact fault developed in the gearbox of an oil refinery to further verify
its effectiveness. One large oxygen separation and compression unit of this
oil refinery was operating with severely abnormal sounds after an overhaul and
thus suspected to have faults induced in its components.

The schematic of the oxygen separation and compression unit is illustrated in
Figure 13. The characteristic frequencies of interest are listed in Table 6.
The sampling frequency of the acquired vibration signals is 15 kHz, and each
epoch of signal is of length 1024. One acceleration signal measured from
Sensor 5 and its Fourier spectrum are displayed in Figure 14(a) and Figure
14(b), respectively. Peak frequencies of 1479 Hz, 2947 Hz, and 4219 Hz can be
observed in Figure 14(b). In addition, the peak frequencies are all surrounded
by sidebands spaced at about 213 Hz. Note that 213 Hz exactly corresponds to
the rotating frequency of shaft II, as shown in Table 6. Thus, the rotating
frequency of shaft II ($f_{2}=$ 213 Hz) can be utilized as prior knowledge of
fault features.

In the practical shutdown inspection, a rub-impact fault between the end face
of the anti-thrust plate on the pinion and that of the bull gear was found.
The fault was due to the inappropriate installation of the pinion [13]. The
end face of the anti-thrust plate and that of the bull gear were not strictly
parallel due to non-perpendicularity between the pinion and Shaft II. Figure
13 illustrates the rubbing fault occurred between the end face of the anti-
thrust plate on the pinion and that of the bull gear.

In order to capture the potential fault feature in Shaft II, the formula (44)
is implemented with $T=1/f_{2}$ to build the binary weight block (27). More
specifically, in this example, we set $M=8$, $K_{1}=2$, $N_{1}=2$, and
$N_{0}+N_{1}=9$, to establish the the binary weight block $B$ with (26) and
(27), and the regularization parameter $\lambda$ is determined by (46) with an
estimated $\hat{\sigma}_{w}=167$ calculated by (47). Figure 15 shows the
result using proposed method with non-convex (‘atan’) penalty function, where
a periodic sequence of oscillatory transients is obtained with a baseline
almost zero. In the STFT domain, a sequence of periodic clusters can be
observed at about $4$ kHz, which coincides to the peak signature frequency
around $4$ kHz in Figure 14(b). The proposed method not only extracts the
periodic oscillatory fault features, but also gives the oscillation frequency
precisely.

Figure 15: Example 3: Result using non-convex penalty function, (a) estimated
time domain signal, (b) STFT coefficients.

As a comparison, we tested the cyclic spectral analysis method proposed in [3,
4, 5]555The matlab implementation is available online at
http://www.mathworks.com/matlabcentral/fileexchange/48909-cyclic-spectral-
analysis, which is one of the ‘state-of-the-art’ analysis method for vibration
signal based on cyclostationary processes. The results are shown in Figure 16.
More specifically, Figure 16(a) shows the cyclic spectral density, where
several suspicious frequency components can be observed on the density map. In
Figure 16(a), at the location where the period (or cyclic frequency) is close
to the truth (213 Hz), two main components at about 1.5 kHz and 4 kHz
respectively can be observed, which coincides to the result given in Figure
15(c), however there exists several other suspicious components on the map as
well, such as the one close to 3.5 kHz at cyclic frequency close to 1050 Hz
and 1250 Hz.

Note that the cyclic spectral analysis method does not require a period (or
cyclic frequency), wherein a suitable range of possible cyclic frequency is
needed to be defined (here we use 0 to 1500 Hz), therefore as a more fair
compassion, we can assume that all the components observed at 213 Hz are the
trustable results, and ignore other suspicious high density points. In this
case, due to the lack of denoising procedure, the density map are quite
blurred comparing the a clear sparse indication given in Figure 15(b) and (c).
In addition, the proposed method in this work, is able to achieve a sparse
time-frequency spectrum indicating not only the frequency but also the time
information, where the cyclostationary process based method does not have such
feasibility. Moreover, we found in this example, the cyclic spectral coherence
[see Figure 16(b)] fails to indicate any frequency information.

Figure 16: Example 3: (a) Cyclic Spectral density and, (b) Cyclic spectral coherence of the test signal given in Figure 14. Table 6: The relevant characteristic frequencies of oxygen separation and compression unit Rotating speed of motor (r/m) | 2985
---|---
Rotating frequency of shaft I (Hz) | 49.75
Rotating frequency of shaft II (Hz) | 213
Meshing frequency of gearbox (Hz) | 6815.75
Rotating frequency of fan blades (Hz) | 3620.86/4472.83

##  6 Conclusion

This paper proposes an approach using structured sparsity in the STFT domain
to extract oscillatory features caused by faults in rotating machinery. We
model the feature caused by a potential fault as a quasi-periodic sequence of
oscillatory transients, where each transient may contain multiple consistent
frequency components. The proposed approach uses smoothed (convex or non-
convex) penalty functions to promote the sparsity of specially grouped STFT
domain coefficients in an optimization problem. This approach detects sparse
and periodically distributed coefficients, and the periodically oscillatory
fault features can be extracted. Using the estimated STFT coefficients, the
signature frequency of the fault feature can be directly observed due to the
nature of sparsity in the time-frequency domain. To solve the proposed
optimization problem efficiently, an iterative algorithm which combines
majorization-minimization and split augmented Lagrangian shrinkage algorithm
has been illustrated. When the problem is not convex, this algorithm still
converges to a local optimum. The effectiveness of the proposed approach is
verified by simulation data, experimental data, and engineering data. The
results demonstrate the notable effectiveness of the proposed approach in
extracting transients for detecting faults in rotating machines.

##  Appendix A Proof of Theorem 1

###### Proof.

By introducing a variable $u\in{}^{N}$ and adding an equality constraint, we
can split the variable $x$, and then problem (14) is equivalent to

$\displaystyle\\{x^{\star},u^{\star}\\}=$
$\displaystyle\arg\min_{x,u}\Big{\\{}F(x,u)=F_{0}(x)+F_{1}(u)\Big{\\}}$ (A.1)
$\displaystyle\text{subject to }x-u=0.$

This problem has an augmented Lagrangian

$\displaystyle
L_{A}(x,u,\beta)=F(x,u)+\beta^{\mathsf{T}}(u-x)+\frac{\mu}{2}\lVert
u-x\rVert_{2}^{2},$ (A.2)

where $\beta\in{}^{N}$, and $\mu>0$. Moreover, we assume that $L_{A}$ in (A.2)
is finite when $x,u,\beta$ are all finite, and the saddle point
$(x^{\star},u^{\star},\beta^{\star})$ exists and should be finite. This
implies that the saddle point $(x^{\star},u^{\star},\beta^{\star})$ is at
least a local minimizer of problem (A.1). If function
$G_{1}:{}^{N}\times{}^{N}\to\real$ is continuously differentiable and
majorizes [is a upper-bound satisfying (15)] of $F_{1}$, then a upper-bound of
$L_{A}$ can be found as

$\displaystyle
G_{A}((x,u,\beta),(x,v,\gamma))=F_{0}(x)+G_{1}(u,v)+G_{2}((x,u,\beta),\gamma)+\frac{\mu}{2}\lVert
u-x\rVert_{2}^{2},$ (A.3)

where $G_{2}((x,u,\beta),\gamma)$ is a upper-bound of
$\beta^{\mathsf{T}}(u-x)$, defined as

$\displaystyle
G_{2}((x,u,\beta),\gamma):=\beta^{\mathsf{T}}(u-x)+\frac{1}{2\mu}\lVert\beta-\gamma\rVert_{2}^{2}\geq\beta^{\mathsf{T}}(u-x),\quad\text{for
all }\beta,\gamma\in{}^{N}.$ (A.4)

Note that, different to the algorithms proposed in [47, 14, 23], since the
augmented Lagrangian is also a function of the introduced multiplier [denoted
as $\beta$ in (A.2)], BSUM can be utilized directly, so that $G_{A}$ in (A.3)
satisfies the block-wise upper-bound conditions listed in (12), when
vectorizing $\\{x,u,\beta\\}$ as a single variable. Hence, BSUM for three
blocks (Table 1) can be directly applied to find the stationary point of
$L_{A}$. Moreover, since $L_{A}$ goes to extrema only when $\\{x,u,\beta\\}$
goes to extrema, BSUM will converge to a non-extremal stationary point of
$L_{A}$, which is the saddle point (see convergence analysis in Ref. [61]).

As a consequence, we can write the algorithm in Table 1 explicitly
‘minimizing’ $L_{A}$ in (A.2), and it will converge to a saddle point of
$L_{A}$,

$\displaystyle u^{(i+1)}$
$\displaystyle=\arg\min_{u}G_{1}(u,u^{(i)})+\beta^{\mathsf{T}}(u-x)+\frac{\mu}{2}\lVert
u-x\rVert_{2}^{2},$ (A.5a) $\displaystyle x^{(i+1)}$
$\displaystyle=\arg\min_{x}F_{0}(x)+\beta^{\mathsf{T}}(u-x)+\frac{\mu}{2}\lVert
u-x\rVert_{2}^{2},$ (A.5b) $\displaystyle\beta^{(i+1)}$
$\displaystyle=\arg\min_{\beta}\beta^{\mathsf{T}}(u-x)+\frac{1}{2\mu}\lVert\beta-\beta^{(i)}\rVert_{2}^{2}$
(A.5c)

where the step (A.5c) can be solved explicitly as

$\displaystyle\beta^{(i+1)}=\beta^{(i)}+\mu(x-u).$ (A.6)

Similar to SALSA, we can reorganize the steps by introducing a variable
$d\in{}^{N}$, and then after simple variable substitutions (same to change
ALM/MM version I to version II in Section II of Ref. [2]), solving (A.5) for
$\\{x^{\star},u^{\star}\\}$ turns out to be iteratively computing the steps
illustrated in (17). ∎

##  Appendix B Matrix inverse lemma

The matrix inverse lemma has several forms. A common form is

$\displaystyle\left(A+BD^{-1}C\right)^{-1}$
$\displaystyle=A^{-1}-A^{-1}B\left(D+CA^{-1}B\right)^{-1}CA^{-1}.$ (B.1)

##  Appendix C Equivalent formulation of real values

To rewrite the algorithm into real value formulation, we consider the complex
variables as an extended real value array, such that

$\displaystyle\bar{c}=\begin{bmatrix}c_{r}\\\
c_{i}\end{bmatrix}=\begin{bmatrix}\mathsf{Re}(c)\\\
\mathsf{Im}(c)\end{bmatrix}\in{}^{M_{1}\times M_{2}\times 2},$ (C.1)

and correspondingly the original variable $c$ can be recovered by
$c=c_{r}+c_{i}\mathrm{i}$. Additionally, to facilitate further derivations, we
also define all the other complex variables into a formulation of real values,

$\displaystyle\bar{u}=\begin{bmatrix}u_{r}\\\
u_{i}\end{bmatrix},\quad\bar{d}=\begin{bmatrix}d_{r}\\\
d_{i}\end{bmatrix},\quad\bar{v}=\begin{bmatrix}v_{r}\\\ v_{i}\end{bmatrix}.$
(C.2)

Moreover, we also consider the short-time Fourier transform operator
$A^{\mathsf{H}}$ and its inverse operator $A$ in (21) into a formulation that
works on real domain. We define $\bar{A}:{}^{M_{1}\times M_{2}\times
2}\to{}^{N}$ is a operator working equivalently to $A$ when realigning $c$ as
$\bar{c}$ in (C.1), where

$\displaystyle\bar{A}\bar{c}=x,$ (C.3)

and its inverse is to do the transform and then align the real and imaginary
parts of the output,

$\displaystyle\bar{A}^{\mathsf{T}}x=\begin{bmatrix}\mathsf{Re}(A^{\mathsf{H}}x)\\\
\mathsf{Im}(A^{\mathsf{H}}x)\end{bmatrix},$ (C.4)

Hence, the property (40) is still persevered as
$\bar{A}\bar{A}^{\mathsf{T}}=I$. As a consequence, solving the problem
sequence (33) is equivalent to compute the following steps,

$\displaystyle u^{(i+1)}=$
$\displaystyle\left\\{\begin{array}[]{lr}\displaystyle\bar{u}^{(i+1)}=(\bar{c}+\bar{d})./\Big{[}1+\frac{\lambda\bar{r}(\bar{u}^{(i)})}{\mu}\Big{]}\\\\[8.00003pt]
u^{(i+1)}=u_{r}^{(i+1)}+\mathrm{i}u_{i}^{(i+1)}\end{array}\right.$ (C.5c)
$\displaystyle c^{(i+1)}=$
$\displaystyle\left\\{\begin{array}[]{lr}\displaystyle\bar{c}^{(i+1)}=(\bar{u}-\bar{d})+\frac{1}{\mu+1}\bar{A}^{\mathsf{T}}\Big{[}y-\bar{A}(\bar{u}-\bar{d})\Big{]}\\\\[8.00003pt]
c^{(i+1)}=c_{r}^{(i+1)}+\mathrm{i}c_{i}^{(i+1)}\end{array}\right.$ (C.5f)
$\displaystyle d^{(i+1)}=$
$\displaystyle\left\\{\begin{array}[]{lr}\displaystyle\bar{d}^{(i+1)}=\bar{d}^{(i)}-(\bar{u}-\bar{c})\\\\[3.99994pt]
d^{(i+1)}=d_{r}^{(i+1)}+\mathrm{i}d_{i}^{(i+1)}\end{array}\right.$ (C.5i)

where $\bar{r}(\bar{v})\in{}^{M_{1}\times M_{2}\times 2}$ in (C.5c) is defined
as

$\displaystyle[\bar{r}(\bar{v})]_{m_{1},m_{2},j}=\sum_{k_{1}=0}^{K_{1}-1}\sum_{k_{2}=0}^{K_{2}-1}\frac{[B]_{k_{1},k_{2}}^{2}}{\psi(p(\bar{v},B,m_{1}-k_{1},m_{2}-k_{2}))},\quad\text{for
}j=1,2.$ (C.6)

where

$\displaystyle p(\bar{v},B,j_{1},j_{2})=\Big{[}\lVert B\odot
S(v_{r},j_{1},j_{2})\rVert_{2}^{2}+\lVert B\odot
S(v_{i},j_{1},j_{2})\rVert_{2}^{2}\Big{]}^{1/2}.$ (C.7)

Note that when comparing (32) to (C.6),

$\displaystyle[\bar{r}(\bar{v})]_{m_{1},m_{2},j}=[r(v)]_{m_{1},m_{2}},\quad\text{for
}j=1,2.$ (C.8)

Therefore, each step of (42) can be understood as computing corresponding real
and imaginary parts as real values, and then recovering the corresponding
complex results afterwards. When we denote the steps on complex domain as
(42), it does not affect the independency of real and imaginary parts. Hence,
the convergence of using SALMA for this problem is ensured as each step solved
on real domain with proper vectorization equivalently.

## References

  * [1] S. Abbasion, A. Rafsanjani, A. Farshidianfar, and N. Irani.  Rolling element bearings multi-fault classification based on the wavelet denoising and support vector machine.  Mech. Syst. Signal Process., 21(7):2933 – 2945, 2007.
  * [2] M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo.  Fast image recovery using variable splitting and constrained optimization.  IEEE Trans. Image Process., 19(9):2345–2356, September 2010.
  * [3] J. Antoni.  Cyclic spectral analysis in practice.  Mech. Syst. Signal Process., 21(2):597–630, 2007.
  * [4] J. Antoni.  Cyclic spectral analysis of rolling-element bearing signals: Facts and fictions.  J. Sound Vib., 304(3–5):497–529, 2007.
  * [5] J. Antoni.  Cyclostationarity by examples.  Mech. Syst. Signal Process., 23(4):987–1036, 2009.
  * [6] Jerome Antoni.  Fast computation of the kurtogram for the detection of transient faults.  Mech. Syst. Signal Process., 21(1):108–124, 2007.
  * [7] Jérôme Antoni and RB Randall.  The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines.  Mech. Syst. Signal Process., 20(2):308–331, 2006.
  * [8] A. Beck and M. Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM J. Imag. Sci., 2(1):183–202, 2009.
  * [9] A. Beck and L. Tetruashvili.  On the convergence of block coordinate descent type methods.  SIAM Journal on Optimization, 23(4):2037–2060, 2013.
  * [10] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends in Machine Learning, 3(1):1–122, 2011.
  * [11] I. Soltani Bozchalooi and M. Liang.  A smoothness index-guided approach to wavelet parameter selection in signal de-noising and fault detection.  J. Sound Vib., 308(1–2):246–267, 2007.
  * [12] B. Chen, Z. Zhang, Y. Zi, Z. He, and C. Sun.  Detecting of transient vibration signatures using an improved fast spatial–spectral ensemble kurtosis kurtogram and its applications to mechanical signature analysis of short duration data from rotating machinery.  Mech. Syst. Signal Process., 40(1):1–37, 2013.
  * [13] Binqiang Chen, Zhousuo Zhang, Chuang Sun, Bing Li, Yanyang Zi, and Zhengjia 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:275–298, 2012.
  * [14] L. Chen, D. Sun, and K.-C. Toh.  An Efficient Inexact Symmetric Gauss-Seidel Based Majorized ADMM for High-Dimensional Convex Composite Conic Programming.  preprint available at: http://adsabs.harvard.edu/abs/2015arXiv150600741C, June 2015.
  * [15] P.-Y. Chen and I. W. Selesnick.  Group-sparse signal denoising: Non-convex regularization, convex optimization.  IEEE Trans. Signal Process., 62(13):3464–3478, July 2014.
  * [16] P.-Y. Chen and I. W. Selesnick.  Translation-invariant shrinkage/thresholding of group sparse signals.  Signal Processing, 94:476–489, 2014.
  * [17] S. Chen and D. L. Donoho.  Basis pursuit.  In 1994 Conference Record of the Twenty-Eighth Asilomar, Conference on Signals, Systems and Computers, 1994., volume 1, pages 41–44, October 1994.
  * [18] S. Chen, D. L. Donoho, and M. A. Saunders.  Atomic decomposition by basis pursuit.  SIAM J. Sci. Comput., 20(1):33–61, 1998.
  * [19] T.W.S. Chow and F. Gou.  Three phase induction machines asymmetrical faults identification using bispectrum.  IEEE Trans. Energy Conversion, 10(4):688–693, December 1995.
  * [20] R. R. Coifman and D. L. Donoho.  Translation-invariant de-noising.  In Wavelet and statistics, pages 125–150. Springer-Verlag, 1995\.
  * [21] F. Cong, W. Zhong, S. Tong, N. Tang, and J. Chen.  Research of singular value decomposition based on slip matrix for rolling bearing fault diagnosis.  J. Sound Vib., 344(0):447–463, 2015.
  * [22] L. Cui, J. Wang, and S. Lee.  Matching pursuit of an adaptive impulse dictionary for bearing fault diagnosis.  J. Sound Vib., 333(10):2840–2862, 2014.
  * [23] Y. Cui, X. Li, D. Sun, and K.-C. Toh.  On the convergence properties of a majorized ADMM for linearly constrained convex optimization problems with coupled objective functions.  preprint available at: http://adsabs.harvard.edu/abs/2015arXiv150200098C, January 2015.
  * [24] I. Daubechies, M. Defriese, and C. De Mol.  An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.  Commun. Pure Appl. Math, LVII:1413–1457, 2004.
  * [25] Y. Ding and I. W. Selesnick.  Sparsity-based correction of exponential artifacts.  Signal Process., 120:236–248, March 2016.
  * [26] D. Donoho, I. Johnstone, and I. M. Johnstone.  Ideal spatial adaptation by wavelet shrinkage.  Biometrika, 81:425–455, 1993.
  * [27] J. Eckstein and D. P. Bertsekas.  On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators.  Math. Program., 55(3):293–318, June 1992.
  * [28] M. El Hachemi Benbouzid.  A review of induction motors signature analysis as a medium for faults detection.  IEEE Trans. Industrial Electronics, 47(5):984–993, October 2000\.
  * [29] M. Elad, J. Starck, P. Querre, and D. Donoho.  Simultaneous cartoon and texture image inpainting using morphological component analysis (MCA).  J. of Appl. and Comp. Harm. Analysis, 19(3):340–358, November 2005\.
  * [30] Zhipeng Feng, Ming Liang, and Fulei Chu.  Recent advances in time–frequency analysis methods for machinery fault diagnosis: A review with application examples.  Mech. Syst. Signal Process., 38(1):165–205, 2013.
  * [31] Zhipeng Feng and Ming J Zuo.  Vibration signal models for fault diagnosis of planetary gearboxes.  J. Sound Vib., 331(22):4919–4939, 2012.
  * [32] M. Figueiredo, J. Bioucas-Dias, and R. Nowak.  Majorization-minimization algorithms for wavelet-based image restoration.  IEEE Trans. Image Process., 16(12):2980–2991, December 2007.
  * [33] Wei Guo and W Tse Peter.  A novel signal compression method based on optimal ensemble empirical mode decomposition for bearing vibration signals.  J. Sound Vib., 332(2):423–441, 2013.
  * [34] Q. He and X. Wang.  Time-frequency manifold correlation matching for periodic fault identification in rotating machines.  J. Sound Vib., 332(10):2611–2626, 2013.
  * [35] Qingbo He.  Vibration signal classification by wavelet packet energy flow manifold learning.  J. Sound Vib., 332(7):1881–1894, 2013.
  * [36] W. He, Y. Ding, Y. Zi, and I. W. Selesnick.  Sparsity-based algorithm for detecting faults in rotating machines.  Mech. Syst. Signal Process., 72–73:46–64, May 2016.
  * [37] WangPeng He, YanYang Zi, BinQiang Chen, Shuai Wang, and ZhengJia He.  Tunable Q-factor wavelet transform denoising with neighboring coefficients and its application to rotating machinery fault diagnosis.  Sci. China Technol. Sci., 56(8):1956–1965, 2013.
  * [38] Wangpeng He, Yanyang Zi, Binqiang Chen, Feng Wu, and Zhengjia He.  Automatic fault feature extraction of mechanical anomaly on induction motor bearing using ensemble super-wavelet transform.  Mech. Syst. Signal Process., 54-55:457–480, 2015.
  * [39] D. R. Hunter and K. Lange.  A tutorial on MM algorithms.  Amer. Statist., 58:30–37, 2004.
  * [40] W Jiang, S. K. Spurgeon, J. A. Twiddle, F. S. Schlindwein, Y. Feng, and S. Thanagasundram.  A wavelet cluster-based band-pass filtering and envelope demodulation approach with application to fault diagnosis in a dry vacuum pump.  Proc. the Institution of Mechanical Engineers, Part C: J. Mechanical Engineering Science., 221(11):1279–1286, 2007.
  * [41] P. K. Kankar, S. C. Sharma, and S. P. Harsha.  Rolling element bearing fault diagnosis using wavelet transform.  Neurocomputing, 74(10):1638–1645, 2011.
  * [42] C. Kar and A. R. Mohanty.  Monitoring gear vibrations through motor current signature analysis and wavelet transform.  Mech. Syst. Signal Process., 20(1):158–187, 2006.
  * [43] C. Kar and A. R. Mohanty.  Vibration and current transient monitoring for gearbox fault detection using multiresolution fourier transform.  J. Sound Vib., 311(1–2):109 – 132, 2008.
  * [44] M. Lang, H. Guo, J. E. Odegard, C. S. Burrus, and R. O. Wells, Jr.  Noise reduction using an undecimated discrete wavelet transform.  IEEE Signal Processing Letters, 3(1):10–12, January 1996.
  * [45] K. Lange, D. Hunter, and I. Yang.  Optimization transfer using surrogate objective functions.  J. of Comp. Graph. Statist., 9:1–20, 2000.
  * [46] Yaguo Lei, Jing Lin, Zhengjia He, and Ming J Zuo.  A review on empirical mode decomposition in fault diagnosis of rotating machinery.  Mech. Syst. Signal Process., 35(1):108–126, 2013.
  * [47] M. Li, D. Sun, and K.-C. Toh.  A Majorized ADMM with Indefinite Proximal Terms for Linearly Constrained Convex Composite Optimization.  December 2014.
  * [48] Z. Li, Z. He, Y. Zi, and Y. Wang.  Customized wavelet denoising using intra- and inter-scale dependency for bearing fault detection.  J. Sound Vib., 313(1–2):342–359, 2008.
  * [49] M. Liang and I. Soltani Bozchalooi.  An energy operator approach to joint application of amplitude and frequency-demodulations for bearing fault detection.  Mech. Syst. Signal Process., 24(5):1473 – 1494, 2010.
  * [50] M. Liang and H. Faghidi.  Intelligent bearing fault detection by enhanced energy operator.  Expert Syst. Appl., 41(16):7223–7234, 2014.
  * [51] Z.Q. Luo and P. Tseng.  On the convergence of the coordinate descent method for convex differentiable minimization.  Journal of Optimization Theory and Applications, 72(1):7–35, 1992\.
  * [52] J. Ma.  Compressed sensing by inverse scale space and curvelet thresholding.  Applied Mathematics and Computation, 206(2):980–988, December 2008\.
  * [53] S. Mallat and Z. Zhang.  Matching pursuits with time-frequency dictionaries.  IEEE Trans. Signal Process., 41:3397–3415, December 1993.
  * [54] G. L. McDonald, Q. Zhao, and M. J. Zuo.  Maximum correlated kurtosis deconvolution and application on gear tooth chip fault detection.  Mech. Syst. Signal Process., 33(0):237–255, 2012.
  * [55] X. Ning and I. W. Selesnick.  ECG enhancement and QRS detection based on sparse derivatives.  Biomedical Signal Processing and Control, 8(6):713–723, 2013.
  * [56] J. P. Oliveira, J. M. Bioucas-Dias, and M. A. T. Figueiredo.  Adaptive total variation image deblurring: A majorization–minimization approach.  Signal Processing, 89(9):1683–1693, 2009.
  * [57] S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin.  An iterative regularization method for total variation-based image restoration.  Multiscale Modeling & Simulation, 4(2):460–489, 2005.
  * [58] Y. Qin, Y. Mao, and B. Tang.  Vibration signal component separation by iteratively using basis pursuit and its application in mechanical fault detection.  J. Sound Vib., 332(20):5217–5235, 2013.
  * [59] Robert B Randall and Jerome Antoni.  Rolling element bearing diagnostics – a tutorial.  Mech. Syst. Signal Process., 25(2):485–520, 2011.
  * [60] M. Razaviyayn, H. Baligh, A. Callard, and Z.-Q. Luo.  Joint user grouping and transceiver design in a MIMO interfering broadcast channel.  IEEE Trans. Signal Process., 62(1):85–94, January 2014.
  * [61] M. Razaviyayn, M. Hong, and Z. Luo.  A unified convergence analysis of block successive minimization methods for nonsmooth optimization.  SIAM Journal on Optimization, 23(2):1126–1153, 2013.
  * [62] L. Rudin, S. Osher, and E. Fatemi.  Nonlinear total variation based noise removal algorithms.  Physica D, 60:259–268, 1992.
  * [63] I. W. Selesnick.  Sparse signal representations using the tunable Q-factor wavelet transform.  In Proc. SPIE, volume 8138 (Wavelets and Sparsity XIV), 2011.
  * [64] I. W. Selesnick.  Wavelet transform with tunable Q-factor.  IEEE Trans. Signal Process., 59(8):3560–3575, August 2011.
  * [65] I. W. Selesnick.  Generalized total variation: Tying the knots.  IEEE Signal Processing Letters, 22(11):2009–2013, November 2015\.
  * [66] I. W. Selesnick, H. L. Graber, Y. Ding, T Zhang, and R. L. Barbour.  Transient artifact reduction algorithm (TARA) based on sparse optimization.  IEEE Trans. Signal Process., 62(24):6596–6611, December 2014.
  * [67] I. W. Selesnick, A. Parekh, and I. Bayram.  Convex 1-d total variation denoising with non-convex regularization.  IEEE Signal Processing Letters, 22(2):141–144, February 2015.
  * [68] Y.-T. Sheen.  A complex filter for vibration signal demodulation in bearing defect diagnosis.  J. Sound Vib., 276(1–2):105–119, 2004.
  * [69] Y.-T. Sheen.  An analysis method for the vibration signal with amplitude modulation in a bearing system.  J. Sound Vib., 303(3–5):538–552, 2007.
  * [70] J.-L. Starck, M. Elad, and D. Donoho.  Redundant multiscale transforms and their application for morphological component analysis.  Advances in Imaging and Electron Physics, 132:287–348, 2004.
  * [71] J.-L. Starck, M. Elad, and D. Donoho.  Image decomposition via the combination of sparse representation and a variational approach.  IEEE Trans. Image Process., 14(10):1570–1582, October 2005.
  * [72] W. Su, F. Wang, H. Zhu, Z. Zhang, and Z. Guo.  Rolling element bearing faults diagnosis based on optimal morlet wavelet filter and autocorrelation enhancement.  Mech. Syst. Signal Process., 24(5):1458–1472, 2010.
  * [73] Haifeng Tang, Jin Chen, and Guangming Dong.  Sparse representation based latent components analysis for machinery weak fault detection.  Mechanical Systems and Signal Processing, 46(2):373–388, 2014.
  * [74] D. Wang and Q. Miao.  Smoothness index-guided bayesian inference for determining joint posterior probability distributions of anti-symmetric real laplace wavelet parameters for identification of different bearing faults.  J. Sound Vib., 345(0):250 – 266, 2015.
  * [75] J. Wang and Q. He.  Exchanged ridge demodulation of time-scale manifold for enhanced fault diagnosis of rotating machinery.  J. Sound Vib., 333(11):2450–2464, 2014.
  * [76] W. Wang.  Early detection of gear tooth cracking using the resonance demodulation technique.  Mech. Syst. Signal Process., 15(5):887–903, 2001.
  * [77] J. Xu and S. Osher.  Iterative regularization and nonlinear inverse scale space applied to wavelet-based denoising.  IEEE Trans. Image Process., 16(2):534–544, February 2007.
  * [78] Ruqiang Yan and Robert X Gao.  Harmonic wavelet-based data filtering for enhanced machine defect identification.  J. Sound Vib., 329(15):3203–3217, 2010.
  * [79] Ruqiang Yan, Robert X Gao, and Xuefeng Chen.  Wavelets for fault diagnosis of rotary machines: A review with applications.  Signal Process., 96:1–15, 2014.
  * [80] H. Yang, J. Mathew, and L. Ma.  Fault diagnosis of rolling element bearings using basis pursuit.  Mech. Syst. Signal Process., 19(2):341–356, 2005.
  * [81] H.-T. Yang and C.-C. Liao.  A de-noising scheme for enhancing wavelet-based power quality monitoring system.  IEEE Trans. Power Delivery, 16(3):353–360, July 2001.
  * [82] W. Yin, S. Osher, D. Goldfarb, and J. Darbon.  Bregman iterative algorithms for $\ell_{1}$-minimization with applications to compressed sensing.  SIAM Journal on Imaging Sciences, 1(1):143–168, 2008.
  * [83] X. Zhang, M. Burger, X. Bresson, and S. Osher.  Bregmanized nonlocal regularization for deconvolution and sparse reconstruction.  SIAM Journal on Imaging Sciences, 3(3):253–276, 2010.
  * [84] X. Zhang, M. Burger, and S. Osher.  A unified primal-dual algorithm framework based on bregman iteration.  Journal of Scientific Computing, 46(1):20–46, 2010.
  * [85] Yongxiang Zhang and RB Randall.  Rolling element bearing fault diagnosis based on the combination of genetic algorithms and fast kurtogram.  Mech. Syst. Signal Process., 23(5):1509–1517, 2009.


as shown in eq. (49). Again using standard thermodynamic relations, one finds
the energy and entropy are given by

$\displaystyle E^{\rm hot}$ $\displaystyle=$ $\displaystyle-T_{\textrm{\tiny
DCFT}}^{2}\frac{\partial I^{\textrm{\tiny hot}}_{\textrm{\tiny E}}}{\partial
T_{\textrm{\tiny DCFT}}}=\frac{\pi\,c}{6}\,T_{\textrm{\tiny
DCFT}}^{2}\,\left(L_{1}+L_{2}\right)\,,$ (310) $\displaystyle S^{\rm hot}$
$\displaystyle=$ $\displaystyle-\frac{\partial(T_{\textrm{\tiny
DCFT}}\,I^{\textrm{\tiny hot}}_{\textrm{\tiny E}})}{\partial T_{\textrm{\tiny
DCFT}}}=\frac{\pi\,c}{3}\,T_{\textrm{\tiny
DCFT}}\,\left(L_{1}+L_{2}\right)+4\log g\,.$

Hence we see that both regions contribute a thermal energy and entropy
(proportional to $L_{1}$ and $L_{2}$) that is characteristic of the deconfined
phase of the CFT. The two defects do not contribute to the energy but each
contribute a constant defect entropy $S_{\rm def}=2\log g$ – recall eq. (7).

## References

  * (1) S. Ryu and T. Takayanagi, _Holographic derivation of entanglement entropy from AdS/CFT_ , _Phys. Rev. Lett._ 96 (2006) 181602 [hep-th/0603001].
  * (2) S. Ryu and T. Takayanagi, _Aspects of Holographic Entanglement Entropy_ , _JHEP_ 08 (2006) 045 [hep-th/0605073].
  * (3) V. E. Hubeny, M. Rangamani and T. Takayanagi, _A Covariant holographic entanglement entropy proposal_ , _JHEP_ 07 (2007) 062 [0705.0016].
  * (4) A. Almheiri, R. Mahajan, J. Maldacena and Y. Zhao, _The Page curve of Hawking radiation from semiclassical geometry_ , _JHEP_ 03 (2020) 149 [1908.10996].
  * (5) N. Engelhardt and A. C. Wall, _Quantum Extremal Surfaces: Holographic Entanglement Entropy beyond the Classical Regime_ , _JHEP_ 01 (2015) 073 [1408.3203].
  * (6) T. Faulkner, A. Lewkowycz and J. Maldacena, _Quantum corrections to holographic entanglement entropy_ , _JHEP_ 11 (2013) 074 [1307.2892].
  * (7) A. Almheiri, N. Engelhardt, D. Marolf and H. Maxfield, _The entropy of bulk quantum fields and the entanglement wedge of an evaporating black hole_ , _JHEP_ 12 (2019) 063 [1905.08762].
  * (8) G. Penington, _Entanglement Wedge Reconstruction and the Information Paradox_ , 1905.08255.
  * (9) G. Penington, S. H. Shenker, D. Stanford and Z. Yang, _Replica wormholes and the black hole interior_ , 1911.11977.
  * (10) A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian and A. Tajdini, _Replica Wormholes and the Entropy of Hawking Radiation_ , _JHEP_ 05 (2020) 013 [1911.12333].
  * (11) V. Balasubramanian, A. Kar and T. Ugajin, _Entanglement between two gravitating universes_ , _Class. Quant. Grav._ 39 (2022) 174001 [2104.13383].
  * (12) A. Miyata and T. Ugajin, _Entanglement between two evaporating black holes_ , _JHEP_ 09 (2022) 009 [2111.11688].
  * (13) L. Anderson, O. Parrikar and R. M. Soni, _Islands with gravitating baths: towards ER = EPR_ , _JHEP_ 21 (2020) 226 [2103.14746].
  * (14) H. Geng, A. Karch, C. Perez-Pardavila, S. Raju, L. Randall, M. Riojas et al., _Information Transfer with a Gravitating Bath_ , 2012.04671.
  * (15) H. Geng, S. Lüst, R. K. Mishra and D. Wakeham, _Holographic BCFTs and Communicating Black Holes_ , 2104.07039.
  * (16) V. Balasubramanian, B. Craps, M. Khramtsov and E. Shaghoulian, _Submerging islands through thermalization_ , _JHEP_ 10 (2021) 048 [2107.14746].
  * (17) Y. Liu, Z.-Y. Xian, C. Peng and Y. Ling, _Black holes entangled by radiation_ , _JHEP_ 09 (2022) 179 [2205.14596].
  * (18) V. Balasubramanian, Y. Nomura and T. Ugajin, _de Sitter space is sometimes not empty_ , 2308.09748.
  * (19) V. Balasubramanian, A. Kar and T. Ugajin, _Entanglement between two disjoint universes_ , _JHEP_ 02 (2021) 136 [2008.05274].
  * (20) V. Balasubramanian, A. Kar and T. Ugajin, _Islands in de Sitter space_ , _JHEP_ 02 (2021) 072 [2008.05275].
  * (21) J. Maldacena and L. Susskind, _Cool horizons for entangled black holes_ , _Fortsch. Phys._ 61 (2013) 781 [1306.0533].
  * (22) H. Z. Chen, R. C. Myers, D. Neuenfeld, I. A. Reyes and J. Sandor, _Quantum Extremal Islands Made Easy, Part I: Entanglement on the Brane_ , _JHEP_ 10 (2020) 166 [2006.04851].
  * (23) H. Z. Chen, R. C. Myers, D. Neuenfeld, I. A. Reyes and J. Sandor, _Quantum Extremal Islands Made Easy, Part II: Black Holes on the Brane_ , _JHEP_ 12 (2020) 025 [2010.00018].
  * (24) J. Hernandez, R. C. Myers and S.-M. Ruan, _Quantum extremal islands made easy. Part III. Complexity on the brane_ , _JHEP_ 02 (2021) 173 [2010.16398].
  * (25) G. Grimaldi, J. Hernandez and R. C. Myers, _Quantum extremal islands made easy. Part IV. Massive black holes on the brane_ , _JHEP_ 03 (2022) 136 [2202.00679].
  * (26) H. Z. Chen, Z. Fisher, J. Hernandez, R. C. Myers and S.-M. Ruan, _Information Flow in Black Hole Evaporation_ , _JHEP_ 03 (2020) 152 [1911.03402].
  * (27) M. Rozali, J. Sully, M. Van Raamsdonk, C. Waddell and D. Wakeham, _Information radiation in BCFT models of black holes_ , 1910.12836.
  * (28) V. Balasubramanian, A. Kar, O. Parrikar, G. Sárosi and T. Ugajin, _Geometric secret sharing in a model of Hawking radiation_ , 2003.05448.
  * (29) H. Z. Chen, Z. Fisher, J. Hernandez, R. C. Myers and S.-M. Ruan, _Evaporating Black Holes Coupled to a Thermal Bath_ , _JHEP_ 01 (2021) 065 [2007.11658].
  * (30) E. Caceres, A. Kundu, A. K. Patra and S. Shashi, _Warped Information and Entanglement Islands in AdS/WCFT_ , 2012.05425.
  * (31) C. Krishnan, _Critical Islands_ , 2007.06551.
  * (32) J. Sully, M. Van Raamsdonk and D. Wakeham, _BCFT entanglement entropy at large central charge and the black hole interior_ , 2004.13088.
  * (33) H. Omiya and Z. Wei, _Causal structures and nonlocality in double holography_ , _JHEP_ 07 (2022) 128 [2107.01219].
  * (34) D. Neuenfeld, _Homology conditions for RT surfaces in double holography_ , _Class. Quant. Grav._ 39 (2022) 075009 [2105.01130].
  * (35) D. Neuenfeld, _The Dictionary for Double Holography and Graviton Masses in d Dimensions_ , 2104.02801.
  * (36) J. Chu, F. Deng and Y. Zhou, _Page curve from defect extremal surface and island in higher dimensions_ , _JHEP_ 10 (2021) 149 [2105.09106].
  * (37) T. Li, M.-K. Yuan and Y. Zhou, _Defect extremal surface for reflected entropy_ , _JHEP_ 01 (2022) 018 [2108.08544].
  * (38) Z. Wang, Z. Xu, S. Zhou and Y. Zhou, _Partial reduction and cosmology at defect brane_ , _JHEP_ 05 (2022) 049 [2112.13782].
  * (39) H. Geng, Y. Nomura and H.-Y. Sun, _An Information Paradox and Its Resolution in de Sitter Holography_ , 2103.07477.
  * (40) Y. Ling, P. Liu, Y. Liu, C. Niu, Z.-Y. Xian and C.-Y. Zhang, _Reflected entropy in double holography_ , _JHEP_ 02 (2022) 037 [2109.09243].
  * (41) H. Geng, A. Karch, C. Perez-Pardavila, S. Raju, L. Randall, M. Riojas et al., _Jackiw-Teitelboim Gravity from the Karch-Randall Braneworld_ , _Phys. Rev. Lett._ 129 (2022) 231601 [2206.04695].
  * (42) K. Suzuki and T. Takayanagi, _BCFT and Islands in two dimensions_ , _JHEP_ 06 (2022) 095 [2202.08462].
  * (43) A. Karch, H. Sun and C. F. Uhlemann, _Double holography in string theory_ , _JHEP_ 10 (2022) 012 [2206.11292].
  * (44) F. Deng, Y.-S. An and Y. Zhou, _JT gravity from partial reduction and defect extremal surface_ , _JHEP_ 02 (2023) 219 [2206.09609].
  * (45) T. Anous, M. Meineri, P. Pelliconi and J. Sonner, _Sailing past the End of the World and discovering the Island_ , _SciPost Phys._ 13 (2022) 075 [2202.11718].
  * (46) D. Neuenfeld and M. Srivastava, _On the causality paradox and the Karch-Randall braneworld as an EFT_ , _JHEP_ 10 (2023) 164 [2307.10392].
  * (47) J.-C. Chang, S. He, Y.-X. Liu and L. Zhao, _Island formula in Planck brane_ , _JHEP_ 11 (2023) 006 [2308.03645].
  * (48) Y. Liu, Q. Chen, Y. Ling, C. Peng, Y. Tian and Z.-Y. Xian, _Entanglement of defect subregions in double holography_ , 2312.08025.
  * (49) J. K. Basak, D. Basu, V. Malvimat, H. Parihar and G. Sengupta, _Holographic Reflected Entropy and Islands in Interface CFTs_ , 2312.12512.
  * (50) H.-S. Jeong, K.-Y. Kim and Y.-W. Sun, _Entanglement entropy analysis of dyonic black holes using doubly holographic theory_ , _Phys. Rev. D_ 108 (2023) 126016 [2305.18122].
  * (51) T. Kawamoto, S.-M. Ruan and T. Takayanagia, _Gluing AdS/CFT_ , _JHEP_ 07 (2023) 080 [2303.01247].
  * (52) C. Bachas and V. Papadopoulos, _Phases of Holographic Interfaces_ , _JHEP_ 04 (2021) 262 [2101.12529].
  * (53) M. Freedman and M. Headrick, _Bit threads and holographic entanglement_ , _Commun. Math. Phys._ 352 (2017) 407 [1604.00354].
  * (54) B. Yoshida and A. Kitaev, _Efficient decoding for the Hayden-Preskill protocol_ , 1710.03363.
  * (55) I. Affleck and A. W. W. Ludwig, _Universal noninteger ’ground state degeneracy’ in critical quantum systems_ , _Phys. Rev. Lett._ 67 (1991) 161.
  * (56) T. Takayanagi, _Holographic Dual of BCFT_ , _Phys. Rev. Lett._ 107 (2011) 101602 [1105.5165].
  * (57) M. Fujita, T. Takayanagi and E. Tonni, _Aspects of AdS/BCFT_ , _JHEP_ 11 (2011) 043 [1108.5152].
  * (58) W. Israel, _Singular hypersurfaces and thin shells in general relativity_ , _Nuovo Cim._ B44S10 (1966) 1.
  * (59) R. Emparan, C. V. Johnson and R. C. Myers, _Surface terms as counterterms in the AdS / CFT correspondence_ , _Phys. Rev. D_ 60 (1999) 104001 [hep-th/9903238].
  * (60) G. Hayward, _Gravitational action for space-times with nonsmooth boundaries_ , _Phys. Rev. D_ 47 (1993) 3275.
  * (61) D. Brill and G. Hayward, _Is the gravitational action additive?_ , _Phys. Rev. D_ 50 (1994) 4914 [gr-qc/9403018].
  * (62) C. Bachas and I. Brunner, _Fusion of conformal interfaces_ , _JHEP_ 02 (2008) 085 [0712.0076].
  * (63) T. Hartman, C. A. Keller and B. Stoica, _Universal Spectrum of 2d Conformal Field Theory in the Large c Limit_ , _JHEP_ 09 (2014) 118 [1405.5137].
  * (64) S. Cooper, M. Rozali, B. Swingle, M. Van Raamsdonk, C. Waddell and D. Wakeham, _Black Hole Microstate Cosmology_ , _JHEP_ 07 (2019) 065 [1810.10601].
  * (65) P. Hayden and J. Preskill, _Black holes as mirrors: Quantum information in random subsystems_ , _JHEP_ 09 (2007) 120 [0708.4025].
  * (66) Y. Nakayama, A. Miyata and T. Ugajin, _The Petz (lite) recovery map for the scrambling channel_ , _PTEP_ 2023 (2023) 123B04 [2310.18991].
  * (67) B. Yoshida, _Unpublished note_ , .
  * (68) K. Jensen and A. Karch, _Holographic Dual of an Einstein-Podolsky-Rosen Pair has a Wormhole_ , _Phys. Rev. Lett._ 111 (2013) 211602 [1307.1132].
  * (69) J. Sonner, _Holographic Schwinger Effect and the Geometry of Entanglement_ , _Phys. Rev. Lett._ 111 (2013) 211603 [1307.6850].
  * (70) V. E. Hubeny and G. W. Semenoff, _String worldsheet for accelerating quark_ , _JHEP_ 10 (2015) 071 [1410.1171].
  * (71) T. Ugajin, _Two dimensional quantum quenches and holography_ , 1311.2562.
  * (72) J. Erdmenger, D. Fernandez, M. Flory, E. Megias, A.-K. Straub and P. Witkowski, _Time evolution of entanglement for holographic steady state formation_ , _JHEP_ 10 (2017) 034 [1705.04696].
  * (73) S. S. Gubser, I. R. Klebanov and A. M. Polyakov, _A Semiclassical limit of the gauge / string correspondence_ , _Nucl. Phys. B_ 636 (2002) 99 [hep-th/0204051].
  * (74) J. Kim and M. Porrati, _Long string dynamics in pure gravity on AdS 3_, _J. Exp. Theor. Phys._ 120 (2015) 477 [1410.3424].
  * (75) J. Kim and M. Porrati, _More on long string dynamics in gravity on AdS 3 : Spinning strings and rotating BTZ black holes_, _Phys. Rev. D_ 91 (2015) 124061 [1503.06875].
  * (76) C. Bachas, Z. Chen and V. Papadopoulos, _Steady states of holographic interfaces_ , _JHEP_ 11 (2021) 095 [2107.00965].
  * (77) H. Maxfield and Z. Wang, _Gravitating spinning strings in AdS 3_, _JHEP_ 07 (2022) 075 [2203.02492].
  * (78) V. Papadopoulos, _Membranes, holography, and quantum information_ , 2310.18521.
  * (79) I. Akal, Y. Kusuki, N. Shiba, T. Takayanagi and Z. Wei, _Holographic moving mirrors_ , _Class. Quant. Grav._ 38 (2021) 224001 [2106.11179].
  * (80) I. Akal, Y. Kusuki, T. Takayanagi and Z. Wei, _Codimension two holography for wedges_ , _Phys. Rev. D_ 102 (2020) 126007 [2007.06800].
  * (81) R. Emparan, G. T. Horowitz and R. C. Myers, _Exact description of black holes on branes_ , _JHEP_ 01 (2000) 007 [hep-th/9911043].
  * (82) R. Emparan, G. T. Horowitz and R. C. Myers, _Exact description of black holes on branes. 2. Comparison with BTZ black holes and black strings_ , _JHEP_ 01 (2000) 021 [hep-th/9912135].
  * (83) R. Emparan, R. Gregory and C. Santos, _Black holes on thick branes_ , _Phys. Rev. D_ 63 (2001) 104022 [hep-th/0012100].
  * (84) P. Kanti, _Black holes in theories with large extra dimensions: A Review_ , _Int. J. Mod. Phys. A_ 19 (2004) 4899 [hep-ph/0402168].
  * (85) A. S. Majumdar and N. Mukherjee, _Braneworld black holes in cosmology and astrophysics_ , _Int. J. Mod. Phys. D_ 14 (2005) 1095 [astro-ph/0503473].
  * (86) A. L. Fitzpatrick, L. Randall and T. Wiseman, _On the existence and dynamics of braneworld black holes_ , _JHEP_ 11 (2006) 033 [hep-th/0608208].
  * (87) S. Creek, R. Gregory, P. Kanti and B. Mistry, _Braneworld stars and black holes_ , _Class. Quant. Grav._ 23 (2006) 6633 [hep-th/0606006].
  * (88) R. Gregory, _Braneworld black holes_ , _Lect. Notes Phys._ 769 (2009) 259 [0804.2595].
  * (89) P. Figueras and T. Wiseman, _Gravity and large black holes in Randall-Sundrum II braneworlds_ , _Phys. Rev. Lett._ 107 (2011) 081101 [1105.2558].
  * (90) I. Akal, T. Kawamoto, S.-M. Ruan, T. Takayanagi and Z. Wei, _Page curve under final state projection_ , _Phys. Rev. D_ 105 (2022) 126026 [2112.08433].
  * (91) T. Kawamoto, S.-M. Ruan, Y.-k. Suzuki and T. Takayanagi, _A half de Sitter holography_ , _JHEP_ 10 (2023) 137 [2306.07575].
  * (92) V. Balasubramanian and P. Kraus, _A Stress tensor for Anti-de Sitter gravity_ , _Commun. Math. Phys._ 208 (1999) 413 [hep-th/9902121].