# Sparsity-based Correction of Exponential Artifacts††thanks: Preprint
submitted to Signal Processing (Manuscript)††thanks: Corresponding author: Yin
Ding<EMAIL_ADDRESS>
Yin Ding and Ivan W. Selesnick
New York University School of Engineering
6 MetroTech Center, Brooklyn, NY 11201, USA

###### Abstract

This paper describes an exponential transient excision algorithm (ETEA). In
biomedical time series analysis, e.g., in vivo neural recording and
electrocorticography (ECoG), some measurement artifacts take the form of
piecewise exponential transients. The proposed method is formulated as an
unconstrained convex optimization problem, regularized by smoothed
$\ell_{1}$-norm penalty function, which can be solved by majorization-
minimization (MM) method. With a slight modification of the regularizer, ETEA
can also suppress more irregular piecewise smooth artifacts, especially,
ocular artifacts (OA) in electroencephalography (EEG) data. Examples of
synthetic signal, EEG data, and ECoG data are presented to illustrate the
proposed algorithms.

##  1 Introduction

This work is motivated by the problem of suppressing various types of
artifacts in recordings of neural activity. In a recent study [25], typical
artifacts in in vivo neural recordings are classified into four types (Type 0
to 3, see Section 2.2 and Figure 3 in [25]). This classification covers many
artifacts in the scope of human brain activity recordings, e.g.,
electroencephalography (EEG) and electrocorticography (ECoG). In this paper,
we consider the suppression of Type 0 and Type 1 artifacts. For the purpose of
flexibility and generality, we redefine them in terms of morphological
characteristics:

  * •

Type 0: a smooth protuberance that can be modeled as $\hat{x}(t)=te^{-\alpha
t}$, when $t\geqslant t_{0}$.

  * •

Type 1: an abrupt jump followed by an exponential decay that can be modeled as
$\hat{x}(t)=e^{-\alpha t}$, when $t\geqslant t_{0}$.

Figure 1 shows examples of the two types of artifacts. We do not consider the
other two types in this work because our previous works have addressed
efficient algorithms to remove such artifacts. For instance, low-pass
filtering/total variation denoising ($\operatorname{\text{LPF/TVD}}$) [61]
suppresses Type 2 artifacts (Figure 1c), and lowpass filtering/compound sparse
denoising (LPF/CSD) [60, 61] can remove sparse and blocky spikes (Type 3 shown
in Figure 1d).

Figure 1: Examples of (a) Type 0 artifact, and (b) Type 1 artifact, and (c)
Type 2 artifact, and (d) Type 3 artifact.

The approach proposed in this paper is based on an optimization problem
intended to capture the primary morphological characteristics of the artifacts
using sparsity-inducing regularization. To formulate the problem, we model the
observed time series as

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

where $f$ is a lowpass signal, $x$ is a piecewise smooth transient signal
(i.e., Type 0 or Type 1 artifacts), and $w$ is stationary white Gaussian
noise. More specifically, $f$ is assumed to be restricted to a certain range
of low frequencies. In other words, $\mathbf{H}(f)\approx\mathbf{0}$, where
$\mathbf{H}$ is a high-pass filter. Note that in the signal model (1),
conventional LTI filtering is not suitable to estimate either $f$ or $x$ from
$y$, because component $x$, as a piecewise smooth signal comprised of
transients, is not band limited.

In order to estimate the components, we combine LTI filtering with sparsity-
based techniques. We formulate an optimization problem for both decomposition
and denoising. A computationally efficient algorithm is derived to solve the
optimization problem, based on the theory of majorization-minimization (MM)
[15, 24, 36].

In addition, this paper specifies how to generate a smoothed penalty function
and its majorizer from a non-smooth one, in order to overcome a numerical
issue that arises when the penalty function is not differentiable.

Table 1: Sparsity-promoting penalty functions ($a>0$).

Penalty | $\phi(u)$ | $\phi_{\epsilon}(u)$ | $\psi(u)=u/\phi_{\epsilon}^{\prime}(u)\quad$
---|---|---|---
abs | $|u|$ | $\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)$
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)$

###  1.1 Related works

Some recent works recover signals with transients by various algorithms. In
[19], a slowly varying signal is modeled as a local polynomial and an
optimization problem using Tikhonov regularization is formulated to capture
it. In [47], the slowly varying trend is modeled as a higher-order sparse-
derivative signal (e.g., the third-order derivative is sparse).

Instead of estimating the slowly varying component via regularization, the
$\operatorname{\text{LPF/TVD}}$ method [61] estimates a lowpass component by
LTI filtering and a piecewise constant component by optimization. In this
case, an optimization problem is formulated to estimate the piecewise constant
component. The approach proposed here uses a similar technique to recover the
lowpass component, but in contrast to $\operatorname{\text{LPF/TVD}}$, it is
more general — the regularization is more flexible with a tunable parameter,
so that $\operatorname{\text{LPF/TVD}}$ can be considered as a special case.

Another algorithm related to the approach taken in this paper is the transient
artifact reduction algorithm (TARA) [60] which is utilized to suppress
additive piecewise constant artifacts and spikes (similar to a hybrid of Type
2 and Type 3 artifact). The approaches proposed in this work target different
types of artifacts (Type 0 and Type 1) and applied in different applications.

The irregularity of Type 0 transients leads to a more complicated artifact
removal problem, where the artifact are irregular fluctuations. A typical
example in EEG is ocular artifacts (OA) caused by the blink and/or movement of
eyes. To suppress OA, there are approaches based on empirical mode
decomposition (EMD) [41, 42, 43, 69], and on independent component analysis
(ICA) methods [1, 12, 20, 26, 37, 48]. The concept of spatial-frequency in
acoustic analysis is also used to remove OA from multichannel signals [44,
64]. In this work, we present a new method to suppress ocular artifacts by
proposing a specific model and using sparse optimization.

This paper adopts a regularizer inspired by the generalized 1-D total
variation [27], wherein the derivative operator in conventional total
variation regularizer is generalized to a recursive filter. The regularizers
adopted in $\operatorname{\text{ETEA}}$ and second-order
$\operatorname{\text{ETEA}}$ coincide with first-order and second-order cases
of generalized 1-D total variation, respectively. Some differences to the
problem discussed in [27] are as follow. Firstly, the signal model (1) allows
a lowpass baseline as a component, hence, $\operatorname{\text{ETEA}}$ can be
seen as a combination of conventional LTI filtering and generalized 1-D total
variation. Secondly, we consider a formulation in terms of banded matrices,
for computational efficiency. Thirdly, we give optimality conditions of the
proposed problems, and use these conditions as a guide to set the
regularization parameters.

##  2 Preliminaries

###  2.1 Notation

We use bold uppercase letters for matrices, e.g., $\mathbf{A}$ and
$\mathbf{B}$, and bold lowercase letters for vectors, e.g., $\mathbf{x}$ and
$\mathbf{y}$. We use column vectors for one-dimensional series. For example, a
vector $\mathbf{x}\in{}^{N}$ is written as

$\displaystyle\mathbf{x}=\big{[}x(0),\ x(1),\ \cdots,\
x(N-1)\big{]}^{\mathsf{T}}$ (2)

where $[\ \cdot\ ]^{\mathsf{T}}$ denotes matrix transpose. The $\ell_{1}$-norm
and squared $\ell_{2}$-norm of $\mathbf{x}$ are defined as

$\displaystyle\displaystyle\lVert\mathbf{x}\rVert_{1}:=\sum_{n}|x(n)|,\qquad\displaystyle\lVert\mathbf{x}\rVert_{2}^{2}:=\sum_{n}{x^{2}(n)}.$
(3)

The inverse transpose of a matrix is denoted as $\mathbf{A}^{-\mathsf{T}}$.
The first-order difference operator $\mathbf{D}$ of size $(N-1)\times N$ is

$\mathbf{D}:=\begin{bmatrix}-1&1&&&\\\ &-1&1&&\\\ &&\ddots&\ddots&\\\
&&&-1&1\end{bmatrix}.$ (4)

We use $f(x;a)$ to denote a function of $x$ determined by parameter $a$, to
distinguish it from a function with two ordered arguments, e.g., $f(x,y)$.

###  2.2 LTI filter with matrix formulation

A discrete-time filter can be described as

$\displaystyle\sum_{k}a_{k}~{}y(n-k)=\sum_{k}b_{k}~{}x(n-k),$ (5)

with transfer function $H(z)=B(z)/A(z)$. For a finite length signal, (5) can
be rewritten as

$\displaystyle\mathbf{A}\mathbf{y}=\mathbf{B}\mathbf{x},$ (6)

where $\mathbf{A},\mathbf{B}\in{}^{N\times N}$ are both banded matrices. When
$\mathbf{H}$ is a zero-phase filter with $a_{-k}=a_{k}$ and $b_{-k}=b_{k}$,
then $\mathbf{A}$ and $\mathbf{B}$ are both symmetric Toeplitz matrices

$\displaystyle\mathbf{A}=\begin{bmatrix}a_{1}&a_{0}&&\\\ a_{0}&a_{1}&a_{0}&\\\
&\ddots&\ddots&\ddots\\\ &&a_{0}&a_{1}&a_{0}\\\ &&&a_{0}&a_{1}\end{bmatrix},$
(7a) $\displaystyle\mathbf{B}=\begin{bmatrix}b_{1}&b_{0}&\\\
b_{0}&b_{1}&b_{0}\\\ &\ddots&\ddots&\ddots\\\ &&b_{0}&b_{1}&b_{0}\\\
&&&b_{0}&b_{1}\end{bmatrix}.$ (7b)

In this case, disregarding a few samples on both edges, $\mathbf{y}$ can be
re-written as:

$\displaystyle\mathbf{y}=\mathbf{B}\mathbf{A}^{-1}\mathbf{x}=\mathbf{H}\mathbf{x}.$
(8)

A detailed description of such zero-phase filters is given in Section V of
Ref. [61], where two parameters, $d$ and $f_{c}$, determine the order of the
filter and the cut-off frequency, respectively. In this paper, $\mathbf{B}$ is
a square matrix, in contrast to [61]. The simplest case ($d=1$) is given by
(7).

###  2.3 Majorization-Minimization

Majorization-Minimization (MM) [15, 35, 57] is a procedure to replace a
difficult optimization problem by a sequence of simpler ones [15, 24, 36, 57].
Suppose a minimization problem has a objective function $J:{}^{N}\to\real$,
then the majorizer
$G(\mathbf{u},\mathbf{v}):\mathbb{R}^{N}\times\mathbb{R}^{N}\to\mathbb{R}$
satisfies

$\displaystyle G(\mathbf{u},\mathbf{v})$
$\displaystyle=J(\mathbf{u}),\quad\text{when}\ \mathbf{u}=\mathbf{v},$ (9a)
$\displaystyle G(\mathbf{u},\mathbf{v})$ $\displaystyle\geqslant
J(\mathbf{u}),\quad\text{when}\ \mathbf{u}\neq\mathbf{v}.$ (9b)

Then MM iteratively solves

$\displaystyle\mathbf{u}^{(k+1)}=\arg\min_{\mathbf{u}}G(\mathbf{u},\mathbf{u}^{(k)})$
(10)

until convergence, where $k$ is iteration index. The detailed derivation and
proof of convergence of MM are given in [35].

##  3 Majorization of smoothed penalty function

Figure 2: (a) Convex non-smooth penalty function and corresponding smoothed
penalty function. (b) Non-convex and smoothed penalty functions. (c) Majorizer
(gray) of the smoothed penalty function in (a). (d) Majorizer (gray) of the
smoothed non-convex penalty function in (b). For all the figures we set
$\epsilon=0.05,v=0.5,a=2$.

When using sparse-regularized optimization to recover a signal, the
$\ell_{1}$-norm is widely utilized. To further enhance sparsity, some methods
use iteratively re-weighting procedures [4, 34, 45, 63], or a non-convex
pseudo $\ell_{p}$-norm ($0<p<1$), or a mixed-norm in the regularizer [5, 68,
66, 58, 33]. Other non-convex penalty functions can also be used. In [46] a
logarithm penalty function is discussed (‘log’ function in Table 1), and in
[59], an arctangent function is used (‘atan’ function Table 1). However, each
of these functions is non-differentiable at zero, where
$\phi^{\prime}(0^{-})\neq\phi^{\prime}(0^{+})$. A method, to avoid the
discontinuity in the derivative, is to smooth the function.

Consider the smoothed function $\phi_{\epsilon}:\real\to\mathbb{R}_{+}$

$\displaystyle\phi_{\epsilon}(u):=\phi(\rho(u;\epsilon)),$ (11)

where function $\rho:\real\to\mathbb{R}_{+}$ is defined as

$\displaystyle\rho(u;\epsilon):=\sqrt{u^{2}+\epsilon},\quad\epsilon>0.$ (12)

Because $\sqrt{u^{2}+\epsilon}$ is always greater than zero, $\phi_{\epsilon}$
has a continuous first-order derivative

$\displaystyle\phi_{\epsilon}^{\prime}(u)=\frac{u}{\sqrt{u^{2}+\epsilon}}\phi^{\prime}(\sqrt{u^{2}+\epsilon}).$
(13)

Table 1 gives the smoothed penalty functions $\phi_{\epsilon}$ corresponding
to $\phi$ in the first column, and Figure 2(ab) illustrates the penalty
function and its smoothed version. The parameter $\epsilon$ controls the
similarity of $\phi_{\epsilon}$ to non-smooth function $\phi$. When
$\epsilon\to 0$, $\phi_{\epsilon}(u)\approx\phi(u)$. In practice, we set
$\epsilon$ to a very small number (e.g. $10^{-10}$).

###  3.1 Majorization of smoothed penalty function

We assume the non-smooth function $\phi$ satisfies the following conditions:

  1. 1.

$\phi(u)$ is continuous on .

  2. 2.

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

  3. 3.

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

  4. 4.

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

Under such assumptions, we find a quadratic function
$g:\real\times\real\to\real$,

$\displaystyle
g(u,v;\phi)=\frac{\phi^{\prime}(v)}{2v}u^{2}+\phi(v)-\frac{v}{2}\phi^{\prime}(v),$
(14)

majorizing $\phi$ when $v\neq 0$ [10, Lemma 1]. However, $\phi$ might not be
differentiable at zero, e.g., all functions $\phi$ in Table 1. To avoid this
issue, we use (11) to ensure the objective function is continuous and
differentiable. Then the MM method can be applied without the occurrence of
numerical issues (e.g., divide by zero). The following proposition indicates a
suitable majorizer of the smoothed penalty functions.

######  Proposition 1.

If function $\phi$ satisfies the listed four conditions in the previous
paragraph and $g$ in (14) is a majorizer of $\phi$, then

$\displaystyle
g_{\epsilon}(u,v;\phi_{\epsilon})=g(\rho(u;\epsilon),\rho(v;\epsilon);\phi),$
(15)

is majorizer of $\phi_{\epsilon}(u)=\phi(\rho(u;\epsilon))$ for $u,v\in\real$,
and if writing explicitly, then $g_{\epsilon}$ in (15) is

$\displaystyle
g_{\epsilon}(u,v;\phi_{\epsilon})=\frac{u^{2}}{2\psi(v)}+\phi_{\epsilon}(v)-\frac{v}{2}\phi_{\epsilon}^{\prime}(v),$
(16)

where $\psi(v)=v/\phi_{\epsilon}^{\prime}(v)\neq 0$.

The proof of Proposition 1 is given in A. Figure 2c and Figure 2d give
examples of using function (16) to majorize the smoothed penalty functions in
Figure 2a and Figure 2b, respectively.

##  4 Exponential transient excision algorithm

To formulate the problem of exponential transient excision, we define
$\mathbf{R}$ as

$\mathbf{R}:=\begin{bmatrix}-r&1&&\\\ &-r&1&&\\\ &&\ddots&\ddots&\\\
&&&-r&1\end{bmatrix}.$ (17)

The first-order derivative operator $\mathbf{D}$ is a special case of
$\mathbf{R}$ with $r=1$. In this paper we restrict $0<r<1$ to distinguish
them. As a filter, $\mathbf{R}$ can be seen to be a first-order filter
attenuating low frequencies, with a transfer function

$\displaystyle R(z):=1-rz^{-1}.$ (18)

Driving the system $R(z)$ with the step exponential

$x(n)=\begin{cases}0,&n<n_{0},\\\ r^{n-n_{0}},&n\geqslant n_{0},\end{cases}$
(19)

produces an impulse $\delta(n-n_{0})$ as an output. If input signal
$\mathbf{x}$ is a step exponential with rate $r$, i.e., a Type 1 artifact,
then

$\mathbf{v}=\mathbf{R}\mathbf{x}$ (20)

will be sparse.

Note that $r$ should be known beforehand. In practice, we can estimate it from
the time-constant of a Type 1 artifact. Using observation data $y$, if we can
measure the time $N_{0}$ over which a Type 1 transient decays to half its
initial height, then $r$ can be found by solving $r^{N_{0}}=0.5$. We ignore
the influence of the slowly varying baseline, as we assume the exponential
decays much faster. If the transients have an approximately equal decay rate,
then we can use an average measured form multiple transients.

###  4.1 Problem definition

In signal model (1), if an estimate $\hat{\mathbf{x}}$ is known, then we can
estimate $f$ by

$\displaystyle\hat{\mathbf{f}}=(\mathbf{y}-\hat{\mathbf{x}})-\mathbf{H}(\mathbf{y}-\hat{\mathbf{x}}),$
(21)

where $\mathbf{H}=\mathbf{B}\mathbf{A}^{-1}$ is a highpass filter, where the
$\mathbf{A}$ and $\mathbf{B}$ matrices are both banded with a structure as
(7). Correspondingly, noise $w$ is
$\mathbf{y}-(\hat{\mathbf{x}}+\hat{\mathbf{f}})=\mathbf{H}(\mathbf{y}-\hat{\mathbf{x}})$,
which indicates a suitable data fidelity term
$\lVert\mathbf{H}(\mathbf{y}-\mathbf{x})\rVert_{2}^{2}$. Moreover, when the
component $x$ is modeled as (19), we can use $\mathbf{v}=\mathbf{R}\mathbf{x}$
as its sparse representation. Hence, we propose to formulate the estimation of
$x$ from noisy data $y$ as the unconstrained optimization problem:

$\displaystyle\mathbf{x}{{}^{*}}=\arg\min_{\mathbf{x}}\Big{\\{}P_{1}(\mathbf{x})=\lVert\mathbf{H}(\mathbf{y}-\mathbf{x})\rVert_{2}^{2}+\lambda\sum_{n}\phi_{\epsilon}([\mathbf{R}\mathbf{x}]_{n})\Big{\\}},$
(22)

where $\lambda>0$ is the regularization parameter related to noise level.

###  4.2 Algorithm derivation

Using the majorizer of smoothed penalty function (16), and recalling that
$\mathbf{H}=\mathbf{B}\mathbf{A}^{-1}$ illustrated in Section 2.2, we majorize
the objective function (22),

$\displaystyle G(\mathbf{x},\mathbf{z};P_{1})=$ $\displaystyle\
\lVert\mathbf{H}(\mathbf{y}-\mathbf{x})\rVert_{2}^{2}+\lambda\sum_{n}g_{\epsilon}([\mathbf{R}\mathbf{x}]_{n},[\mathbf{R}\mathbf{z}]_{n};\phi_{\epsilon})$
$\displaystyle=$ $\displaystyle\
\lVert\mathbf{B}\mathbf{A}^{-1}(\mathbf{y}-\mathbf{x})\rVert_{2}^{2}+\mathbf{x}^{\mathsf{T}}\mathbf{R}^{\mathsf{T}}\big{[}\mathbf{\Lambda}(\mathbf{R}\mathbf{z})\big{]}\mathbf{R}\mathbf{x}+C$
(23)

where $C$ does not depend on $\mathbf{x}$, and
$\big{[}\mathbf{\Lambda}(\mathbf{R}\mathbf{z})\big{]}\in{}^{(N-1)\times(N-1)}$
is a diagonal matrix

$\displaystyle\big{[}\mathbf{\Lambda}(\mathbf{R}\mathbf{z})\big{]}_{n,n}=\frac{\lambda}{2\psi([\mathbf{R}\mathbf{z}]_{n})}$
(24)

where $\psi(u):=u/\phi_{\epsilon}^{\prime}(u)$, which is listed in Table 1.
Then the minimizer of (4.2) is given by

$\displaystyle\mathbf{x}=\Big{[}\mathbf{A}^{-\mathsf{T}}\mathbf{B}^{\mathsf{T}}\mathbf{B}\mathbf{A}^{-1}+\mathbf{R}^{\mathsf{T}}\big{[}\mathbf{\Lambda}(\mathbf{R}\mathbf{z})\big{]}\mathbf{R}\Big{]}^{-1}\mathbf{A}^{-\mathsf{T}}\mathbf{B}^{\mathsf{T}}\mathbf{B}\mathbf{A}^{-1}\mathbf{y}.$
(25)

Note that, calculating (25) directly requires the solution to a large dense
system of equations. However, since $\mathbf{A}$ is invertable, we can factor
$\mathbf{A}$ out and write (25) as

$\displaystyle\mathbf{x}^{(k+1)}=\mathbf{A}\Big{[}\mathbf{B}^{\mathsf{T}}\mathbf{B}+\mathbf{A}^{\mathsf{T}}\mathbf{R}^{\mathsf{T}}\big{[}\mathbf{\Lambda}(\mathbf{R}\mathbf{x}^{(k)})\big{]}\mathbf{R}\mathbf{A}\Big{]}^{-1}\mathbf{B}^{\mathsf{T}}\mathbf{B}\mathbf{A}^{-1}\mathbf{y},$
(26)

where the MM procedure (10) is adopted with $\mathbf{z}=\mathbf{x}^{(k)}$.
Note that the matrix to be inverted in (26) is banded, so that the solution
can be computed efficiently by fast banded system solvers (e.g., [28, 29]).
Moreover, in (26), all matrix multiplications are between banded ones. Table 2
summarizes the algorithm to solve (22).

Table 2: Exponential Transient Excision Algorithm.
$\displaystyle\text{Input:}~{}\mathbf{y},~{}\lambda,~{}r,~{}\mathbf{A},~{}\mathbf{B}$
$\displaystyle\text{Initialization:}~{}\mathbf{x}\in\mathbb{R}^{N}$
$\displaystyle\mathbf{b}=~{}\mathbf{B}^{\mathsf{T}}\mathbf{B}\mathbf{A}^{-1}\mathbf{y}$
Repeat $\displaystyle\qquad\mathbf{v}=\ \mathbf{R}\mathbf{x}$
$\displaystyle\qquad[\mathbf{\Lambda}]_{n,n}=\frac{\lambda}{2\psi(v(n))}$
$\displaystyle\qquad\mathbf{Q}=\mathbf{B}^{\mathsf{T}}\mathbf{B}+\mathbf{A}^{\mathsf{T}}\mathbf{R}^{\mathsf{T}}\mathbf{\Lambda}\mathbf{R}\mathbf{A}$
$\displaystyle\qquad\mathbf{x}\ =\mathbf{A}\mathbf{Q}^{-1}\mathbf{b}$ Until
convergence
$\displaystyle\mathbf{f}=(\mathbf{y}-\mathbf{x})-\mathbf{H}(\mathbf{y}-\mathbf{x})$
$\displaystyle\text{Return: }\mathbf{x},~{}\mathbf{f}.$

###  4.3 Optimality condition and parameter selection

For problem (22), it is difficult to derive optimality conditions directly.
Here, we consider the optimality condition indirectly via an equivalent
problem, similar to the discussion of optimality condition for total variation
problems in [2, 17] and for $\operatorname{\text{LPF/TVD}}$ in [61].

To facilitate our derivation, we firstly denote the optimal solution of (22)
by $\mathbf{x}^{*}$. We then define $\mathbf{v}^{*}=\mathbf{R}\mathbf{x}^{*}$
where $\mathbf{R}$ is given by (17). We define
$\mathbf{G}\in{}^{N\times(N-1)}$,

$\mathbf{G}:=\left[\begin{array}[]{llccccc}0\\\ 1&0\\\ r&1&0&\\\
r^{2}&r&1&0&\\\ \vdots&&&\ddots&\ddots&\\\ r^{N-3}&\cdots&\cdots&r&1&0\\\
r^{N-2}&\cdots&\cdots&r^{2}&r&1\end{array}\right]$ (28)

which satisfies

$\displaystyle\mathbf{R}\mathbf{G}=\mathbf{I}.$ (29)

The operator $\mathbf{G}$ acts as an inverse filter of $\mathbf{R}$. Using the
above assumptions, we derive the following proposition regarding optimality
conditions for (22).

######  Proposition 2.

If $\mathbf{x}^{*}$ is an optimal solution to (22), then
$\mathbf{v}^{*}=\mathbf{R}\mathbf{x}^{*}$ is an optimal solution to problem

$\displaystyle\mathbf{v}^{*}=\arg\min_{\mathbf{v}}Q(\mathbf{v}),$ (30)

where

$\displaystyle
Q(\mathbf{v})=\lVert\mathbf{H}(\mathbf{y}_{0}-\mathbf{G}\mathbf{v})\rVert_{2}^{2}+\lambda\sum_{n}\phi_{\epsilon}([\mathbf{v}]_{n}).$
(31)

and
$\mathbf{y}_{0}=\mathbf{y}-(\mathbf{x}^{*}-\mathbf{G}\mathbf{R}\mathbf{x}^{*})$.

###### Proof.

We define $\mathbf{x}_{0}\in{}^{N}$ as

$\displaystyle\mathbf{x}_{0}:=\mathbf{G}\mathbf{v}^{*}=\mathbf{G}\mathbf{R}\mathbf{x}^{*}$
(32)

where $\mathbf{x}^{*}$ is the optimal solution to (22). Note that
$\mathbf{x}_{0}$ and $\mathbf{x}^{*}$ share an identical
$\mathbf{v}^{*}=\mathbf{R}\mathbf{x}_{0}=\mathbf{R}\mathbf{x}^{*}$. Moreover,
we define the difference between $\mathbf{x}^{*}$ and $\mathbf{x}_{0}$ as

$\displaystyle\mathbf{d}:=\mathbf{x}^{*}-\mathbf{x}_{0}.$ (33)

It follows that
$\mathbf{R}\mathbf{d}=\mathbf{R}\mathbf{x}^{*}-\mathbf{R}\mathbf{x}_{0}=\mathbf{0}$.
As a consequence, $\mathbf{x}_{0}$ is the optimal solution to problem

$\displaystyle\mathbf{x}_{0}=$
$\displaystyle\arg\min_{\mathbf{x}}P_{1}(\mathbf{u})$ (34a) $\displaystyle\
\text{s.t. }\mathbf{x}=\mathbf{u}-\mathbf{d}.$ (34b)

Note $x_{0}(0)=0$ must be satisfied in (34) because of the definition of
$\mathbf{G}$. Thus, we can write an equivalent problem to (34) using (32),

$\displaystyle\\{\mathbf{x}_{0},\mathbf{v}^{*}\\}=\arg\min_{\mathbf{x},\mathbf{v}}P_{1}(\mathbf{u})$
(35a) $\displaystyle\text{s.t. }\mathbf{x}=\mathbf{u}-\mathbf{d}$ (35b)
$\displaystyle\mathbf{x}=\mathbf{G}\mathbf{v}\quad$ (35c)

where $\\{\mathbf{x}_{0},\mathbf{v}^{*}\\}$ must be the optimal solution. In
problem (35), $\mathbf{u}$ is uniquely determined by linear functions of
$\mathbf{v}$ and $\mathbf{d}$, so problem (35) can be simplified by
substituting the variables,

$\displaystyle\\{\mathbf{x}_{0},\mathbf{v}^{*}\\}=\arg\min_{\mathbf{x},\mathbf{v}}Q(\mathbf{v})$
(36a) $\displaystyle\text{s.t. }\mathbf{x}=\mathbf{G}\mathbf{v}$ (36b)

where $Q(\mathbf{v})=P_{1}(\mathbf{u})$, with
$\mathbf{u}=\mathbf{G}\mathbf{v}+\mathbf{d}$, and can be written explicitly,

$\displaystyle Q(\mathbf{v})$
$\displaystyle=\lVert\mathbf{H}(\mathbf{y}-\mathbf{G}\mathbf{v}-\mathbf{d})\rVert_{2}^{2}+\lambda\sum_{n}\phi_{\epsilon}([\mathbf{R}\mathbf{G}\mathbf{v}+\mathbf{R}\mathbf{d}]_{n}).$
(37)

Because $\mathbf{R}\mathbf{G}=\mathbf{I}$ and
$\mathbf{R}\mathbf{d}=\mathbf{0}$, we simplify (37) as

$\displaystyle Q(\mathbf{v})$
$\displaystyle=\lVert\mathbf{H}(\mathbf{y}_{0}-\mathbf{G}\mathbf{v})\rVert_{2}^{2}+\lambda\sum_{n}\phi_{\epsilon}([\mathbf{v}]_{n})$
(38)

where $\mathbf{y}_{0}=\mathbf{y}-\mathbf{d}$.

In this case, the equality constraint in (36) is redundant. We can solve the
unconstrained problem (30) first, and then compute $\mathbf{x}_{0}$ by (32).
This implies that $\mathbf{v}^{*}$ is the optimal solution to (30) and
satisfies $\nabla Q(\mathbf{v}^{*})=\mathbf{0}$. ∎

Using Proposition 2, instead of problem (22), we alternatively consider the
optimality condition of problem (30) where as long as $\mathbf{x}^{*}$ is an
optimal solution of (22), $\mathbf{v}^{*}$ is an optimal solution of (30).
Moreover, we can rewrite the optimality condition of (30) as

$\displaystyle
p(n)=\lambda\phi_{\epsilon}^{\prime}([\mathbf{v}^{*}]_{n})=\lambda\phi_{\epsilon}^{\prime}([\mathbf{R}\mathbf{x}^{*}]_{n}),$
(39)

where
$\mathbf{p}=2\mathbf{G}^{\mathsf{T}}\mathbf{H}^{\mathsf{T}}\mathbf{H}(\mathbf{y}_{0}-\mathbf{x}_{0})$,
which can be rewritten as

$\displaystyle\mathbf{p}=2\mathbf{G}^{\mathsf{T}}\mathbf{H}^{\mathsf{T}}\mathbf{H}(\mathbf{y}-\mathbf{x}^{*}).$
(40)

Note that, writing $\mathbf{p}$ as (40) does not require the solution to
problem (30) or to compute $\mathbf{d}$ and $\mathbf{y}_{0}$. Therefore,
although we derived an indirect way to verify the optimality condition for
(22), the final procedure is direct.

Setting parameter $\lambda$. The equation (39) can be used as a guide to set
the regularization parameter $\lambda$. Suppose the observation data is
composed of Gaussian noise only, then a proper value of $\lambda$ should make
the solution of (22) almost identically zero. Thus, its sparse representation
$\mathbf{v}^{*}$, given by (20), should be all zero as well. In this case, we
can calculate a vector $\mathbf{q}$, which depends only on noise,

$\displaystyle\mathbf{q}=2\mathbf{G}^{\mathsf{T}}\mathbf{H}^{\mathsf{T}}\mathbf{H}\mathbf{w}.$
(41)

When $\mathbf{p}$ is calculated by (40), we find a constraint that
$-\lambda<p(n)<\lambda$, by the property of $\phi_{\epsilon}$. Furthermore,
since $\mathbf{x}^{*}\approx\mathbf{0}$ is expected when the observation is
pure noise (i.e., $\mathbf{y}=\mathbf{w}$), the values of
$\mathbf{q}\approx\mathbf{p}$ can be considered bounded

$\displaystyle q(n)\in[-\lambda,+\lambda],\quad\text{when
}|v(n)|\leqslant\beta_{\epsilon},\text{ for }n\in\mathbb{Z}_{N}$ (42)

with high probability. The value $\beta_{\epsilon}$ here is related to
$\epsilon$, and since $\epsilon$ is extremely small, it can be simply assumed
to be 0. As a consequence, the value of $\lambda$ needs to satisfy

$\displaystyle\lambda\geqslant\max\big{\\{}|2[\mathbf{G}^{\mathsf{T}}\mathbf{H}^{\mathsf{T}}\mathbf{H}\mathbf{w}]_{n}|,n\in\mathbb{Z}_{N}\big{\\}}.$
(43)

Further, if the statistical property of the noise can be exploited, e.g.,
$w(n)\sim\mathcal{N}(0,\sigma_{w}^{2})$, $\lambda$ can be set statistically,
such that

$\displaystyle\lambda=2.5\sigma_{w}\lVert 2\mathbf{h}_{1}\rVert_{2}$ (44)

where $\sigma_{w}$ is the standard deviation of the noise and $\mathbf{h}_{1}$
is the impulse response corresponding to system
$\mathbf{G}^{\mathsf{T}}\mathbf{H}^{\mathsf{T}}\mathbf{H}$. Note that,
$\mathbf{G}^{\mathsf{T}}\mathbf{H}^{\mathsf{T}}\mathbf{H}$ is a filter with a
transfer function $P_{1}(z)=H^{2}(z)/{R(z)}$. Hence, although we give the
matrix $\mathbf{G}$ to derive the approach to set $\lambda$, it is not
required to compute $\mathbf{h}_{1}$ because the filter $P_{1}(z)$ can be
implemented directly using $H(z)$ in (8), and $R(z)$ in (18).

###  4.4 Example: synthetic signal

Figure 3: Example 1: Simulated data comprising a lowpass component and
exponential transients.

Figure 4: Example 1: The denoising results of (a) lowpass filtering, (b)
$\operatorname{\text{LPF/TVD}}$, (c) $\operatorname{\text{ETEA}}$ (proposed
method).

Figure 5: Example 1: Comparison of denoising and decomposition results using
(a) Wavelet-based method and, (b) $\operatorname{\text{ETEA}}$ with non-convex
regularization.

To illustrate $\operatorname{\text{ETEA}}$, the algorithm is tested on a
simulated signal and compared to other methods. Figure 3 shows the test signal
and its components: a lowpass signal $f$, a transient component $x$, composed
of three exponential decays, and white Gaussian noise with $\sigma_{w}=0.20$.
The noisy observation is shown in Figure 3 in gray.

Denoising. The result of conventional lowpass filtering is shown in Figure 4a.
The result is smooth, but the step edges of each pulse are diminished.
Moreover, if compared to the true low-pass baseline, the result is also
distorted. The result of using $\operatorname{\text{LPF/TVD}}$ [61] is
illustrated in Figure 4b. Although the abrupt jumps have been preserved,
serious staircase effects appear. $\operatorname{\text{LPF/TVD}}$ uses the
first-order derivative operator $\mathbf{D}$ in regularization. The derivative
of an exponential is another exponential, therefore the derivative is unlikely
to have a sufficiently sparse representation. A compensative method is to
widen the bandwidth of the lowpass filter, however, this leads to a noisy
result. The result presented in Figure 4b is selected by tuning all necessary
parameters of $\operatorname{\text{LPF/TVD}}$ to optimize RMSE (root-mean-
square error) value.

The result of the proposed ETEA method is shown in Figure 4c. The step edges
and decay behavior are well preserved (RMSE = 0.057). We use a smoothed
$\ell_{1}$-norm penalty function. In this example, we set $\epsilon=10^{-10}$,
$r=0.94$, filter order parameter $d=1$, cut-off frequency
$f_{c}=0.013\mathrm{~{}cycle/sample}$, and regularization parameter is
calculated by (44) based on the noise principle.

Decomposition. Wavelet-based methods for artifact correction have been
described in [1, 25, 40, 56]. In Figure 5a, we illustrate denoising and
decomposition results obtained using the stationary (un-decimated) wavelet
transform [11] with Haar wavelet using hard-threshold determined by the
$\sqrt{2\log{N}}\sigma_{w}$ thresholding scheme. The denoised result is
further enhanced by the artifact-free denoising method [13, 14] which uses
total variation minimization to overcome pseudo-Gibbs oscillations. Although
the denoising output is relatively smooth and captures the discontinuities,
the decomposed components are both distorted. Compared to the true components
in Figure 3, the lowpass subband deviates from true lowpass component
(especially at about $n=100$). The transient component, reconstructed from
other subbands after denoising, has to compensate for the error. Hence, the
estimated transient component does not have a zero baseline (see bottom of
Figure 5a).

Non-convex penalty functions can be used to improve the result of ETEA. Here
we use the smoothed non-convex logarithm penalty function in Table 1 ($a=2$).
The filtering result is shown in Figure 5b, where discontinuities are more
accurately preserved, and the RMSE is reduced to 0.043 compared with the
result in Figure 4c. We also illustrate the decomposed $f$ and $x$ components
in Figure 5b. $\operatorname{\text{ETEA}}$ recovers both of the components
accurately. The algorithm run-time is not affected by the choice of penalty
function. In this example, 50 iterations of the algorithm takes about 25 ms on
a MacBook Pro 2012 with a 2.7 GHz CPU, implemented in Matlab 8.4.

###  4.5 Example: artifact removal of ECoG data

Conventional EEG studies and most clinical EEG applications are restricted
below 75 Hz [51]. Advanced measuring methods such as ECoG records multichannel
cortical potentials from the micro-electrodes located inside human skull with
a higher sampling rate. In this example, we use a 5-second ECoG signal epoch,
with a sampling rate of 2713 Hz, recorded from a mesial temporal lobe epilepsy
(MTLE) patient.

Figure 6a shows the raw data in gray and the corrected data by ETEA in black.
There are two Type 1 artifacts identified in this 5-second epoch. One is at
about $t=1.20$ second, and the other is at about $t=3.60$ second. In signal
model (1), suppressing component x which represents artifacts, the corrected
data $(f+w)$ should preserve the components. We show the corrected data in
Figure 6a, where the two Type 1 spikes are correctly removed. The decomposed
artifact is illustrated in Figure 6a as well, which adheres a zero-baseline.

Since wavelet-based methods have been successfully applied to suppress
artifacts [1, 25, 40]. a comparison with a wavelet-based method is shown in
Figure 6b. We use the stationary (un-decimated) wavelet transform [11] with
Haar wavelet filter and the non-negative garrote threshold function [18],as
recommended in [25]. The thresholding has been applied to all the subbands
except the lowpass band, and the artifact component is obtained by subtracting
the corrected data from the raw data. As shown, this approach estimates
transient pulses but the estimated pulse adheres to the shape of the Haar
wavelet filter: a positive-negative pulse (see the bottom square box in Figure
6b), but not the true Type 1 artifact in Figure 1b.

To more clearly illustrate the estimated results by the two methods, we show
the details of the corrected data and the artifact from $t=3.50$ to $3.65$
second in dashed line boxes in Figure 1a and Figure 1b, respectively. Since
the wavelet-based method cannot correctly estimate Type 1 artifact in this
case, a ‘bump’ can be observed in Figure 6b to compensate the error. Moreover,
this inappropriate estimation will cause the corrected data to change before
Type 1 artifact actually occurs (see Figure 1b). In contrast, ETEA estimates
the artifact as the abrupt drift with a decay, then it exicises the artifact
without influencing the data before the transient occurs (see Figure 1a).

In this example, for the signal with a length of $13565$ samples, the proposed
algorithm converges within 40 iterations, and takes about 330 ms on a MacBook
Pro 2012 with 2.7 GHz CPU, implemented in Matlab 8.4. The cost function
history is shown in Figure 7.

Figure 6: Example 2: Comparison of artifacts removal for ECoG data. (a) ETEA
(proposed method), (b) wavelet-based method.
Figure 7: Example 2: Cost function history.

##  5 Higher-order ETEA

Figure 8: Example 3: The simulated test data, generated by driving an AR
filter with a sparse sequence and adding noise.

Figure 9: Example 3: Denoising and decomposition results by (a) EMD based
denoising (EMD-CIIT [32]) and (b) $\operatorname{\text{ETEA}}$.

In the previous section, we have presented $\operatorname{\text{ETEA}}$ based
on signal model (1), suitable for Type 1 artifacts, where the component $x$
has discontinuous step exponential transients. Here we illustrate another
version of $\operatorname{\text{ETEA}}$ based on signal model (1), where the
component $x$ models Type 0 artifacts.

First, we consider the operator

$\mathbf{R}_{2}:=\begin{bmatrix}r^{2}&-2r&1&&\\\ &r^{2}&-2r&1&&\\\
&&\ddots&\ddots&\ddots&\\\ &&&r^{2}&-2r&1\end{bmatrix},$ (45)

where $\mathbf{R}_{2}\in{}^{(N-2)\times N}$ has three non-zero coefficients in
each row. As a special case, when $r=1$, $\mathbf{R}_{2}$ is the second-order
difference operator, and $\lVert\mathbf{R}_{2}\mathbf{x}\rVert_{1}$ is the
same as the regularizer used in [30] for $\ell_{1}$ detrending.

Taking $\mathbf{R}_{2}$ as a filter, the transfer function is

$\displaystyle R_{2}(z)=(1-rz^{-1})^{2},$ (46)

which has a double zero at $z=r$. The impulse response of system
$R_{2}^{-1}(z)$ is $(n+1)r^{n}$, when $n\geqslant 0$. It has the same shape as
Figure 1a, which is a suitable model for the Type 0 artifacts in [25].
Therefore, $\mathbf{R}_{2}\mathbf{x}$ is sparse when $x$ is composed of such
piecewise smooth transients. Similar to (22), we formulate the problem

$\displaystyle\mathbf{x}^{*}=\arg\min_{\mathbf{x}}\Big{\\{}P_{2}(\mathbf{x})=\lVert\mathbf{H}(\mathbf{y}-\mathbf{x})\rVert_{2}^{2}+\lambda\sum_{n}\phi_{\epsilon}([\mathbf{R}_{2}\mathbf{x}]_{n})\Big{\\}}$
(47)

to estimate the transient component $x$, so that the corrected data is
$y-\hat{x}$, and the low-pass trend $f$ can be estimated by (21).
Additionally, (47) can be easily solved by ETEA (in Table 2) substituting
$\mathbf{R}$ by $\mathbf{R}_{2}$.

An example of complicated and irregular artifacts is the eye blink/movement
artifacts in EEG data, which may have various morphologies and durations. In
this case, we assume that the artifacts can be estimated by a continuous
piecewise smooth waveform, generated by applying a certain sparse impulse
sequence to $R_{2}^{-1}(z)$ (46). In other words, we broaden our signal model
so that the artifact is not only equivalent to an isolated Type 0 transient as
in [25], but also a superposition of multiple such transients, with some
freedom of scaling and shifting. The problem of estimating $x$ can be solved
by (47) in this case as well.

###  5.1 Example: Simulated data

Figure 8 shows the simulated data and its lowpass and transient components.
The transient component has several pulses with different heights and widths.
They are obtained by feeding a sequence of impulses into system
$R_{2}^{-1}(z)$. The filter is given in (5) with $a_{k}=[1,-2r,r^{2}]$,
$b_{0}=1$, and $r=0.950$.

Empirical mode decomposition (EMD) [22, 23], a powerful method for analyzing
signals, has been successfully utilized in different fields, including
neuroscience, biometrics, speech recognition, electrocardiogram (ECG)
analysis, and fault detection [8, 16, 21, 39, 54, 62, 65, 67].

As a comparison to the proposed approach, we use the EMD based denoising
algorithm in [32], which uses wavelet coefficient thresholding techniques on
decomposed IMFs. More specifically, we use clear first iterative interval
thresholding (CIIT) with smoothly clipped absolute deviation (SCAD) penalty
thresholding [31, 32], with 20 iterations, and the result is shown in Figure
9b. In this example, in order to perform decomposition, among the entire eight
IMFs, the thresholding is performed on IMFs 1-5, their summation is considered
as the transients illustrated in Figure 9a, and IMFs 6-8 are considered as the
estimation of the lowpass component. The EMD based method achieves a decent
denoising performance, but does not accurately estimate the components. For
instance, the simulated data has a smooth dip at about $n=700$ (circled in
Figure 9a). EMD decomposes it into higher IMFs since they are varying slowly,
which degrades the estimation. We may group the lowpass and transient
components differently to avoid this problem, for instance, grouping IMF 1-6
together in order to include more oscillations into transient component, but
this causes other distortion, where the decomposed transient component
contains a lowpass signal and does not adhere to a baseline of zero.

The result obtained using second-order $\operatorname{\text{ETEA}}$ is
illustrated in Figure 9b. $\operatorname{\text{ETEA}}$ estimates both the low-
pass and transient components well, and recovers the signal by $x+f$ precisely
with RMSE = 0.87 (with the smoothed $\ell_{1}$ penalty function). The
regularization parameter $\lambda$ for problem (47) was similar to (44),

$\displaystyle\lambda\approx 2.5\sigma_{w}\lVert 2\mathbf{h}_{2}\rVert_{2},$
(48)

where $\mathbf{h}_{2}$ is the impulse response of system
$H^{2}(z)/{R_{2}(z)}$, and $R_{2}(z)$ is defined in (46). The decomposition is
accurate. There are no compensating waveforms between the estimated $x$ and
$f$ at $n=700$, comparing to the estimation in Figure 9a.

Through numerical experiments, we found that second-order ETEA is not very
sensitive to parameter $r$. Figure 9c shows the RMSE of denoising the data in
Figure 8a, using $r\in[0.85,0.99]$. In this test, all filter parameters
($f_{c}$ and $d$) are the same and $\lambda$ is set by (48). In most cases,
the results are better than EMD-CIIT. In addition, $r$ should not be too small
or very close to 1. If $r$ has to be very small to yield a good estimation of
component $x$, it must fluctuate extremely rapidly, then it must be closer to
a sequence of sparse spikes (i.e., Type 3 artifacts in [25]), which differs
from the signal model. For such a signal, other algorithms are more suitable,
e.g., LPF/CSD [60, 61]. Similarly, if $r$ has to be very close to 1 to fit the
transients, then the transients must be very close to a piecewise linear
signal, which is also not how we model the signal initially. As a consequence,
we suggest to set $r$ in the range $0.90<r<0.98$.

###  5.2 Example: ocular artifacts suppression

Figure 10: Example 4: (a) EEG data with ocular artifacts, and corrected data
using (b) $\operatorname{\text{ETEA}}$ and (c) MEMD.

In this example, we use $\operatorname{\text{ETEA}}$ with $\mathbf{R}_{2}$ to
correct EEG with eye blink/movement artifacts. Figure 10a shows a 10 second
signal from channel Fp1, with sampling rate $f_{s}=256$ Hz, downloaded from
[38]. As a channel located on forehead, Fp1 is very sensitive to the motion of
eyes and eyebrows, and in this example, eye movement artifacts with large
amplitudes are present through the entire signal. Applying second-order
$\operatorname{\text{ETEA}}$ with $r=0.95$, the results for corrected data and
extracted OA are illustrated in Figure 10b.

As a comparison, we use multivariate empirical mode decomposition (MEMD) [52]
to correct the data. MEMD is a recently developed algorithm extending
conventional EMD. It has been used in different aspects of EEG signal analysis
and applications [9, 50, 49], including removing the ocular artifacts (OA)
from multichannel EEG data [53, 42]. In this example, we use 4 EEG channels
(Fp1, Fp2, C3, C4) measured simultaneously as the input, and decompose the
higher-index IMFs (low-frequency subbands) considered to be artifacts [53,
Section V]. More specifically, among all 14 IMFs decomposed in this example,
we use IMF 1-4 as the corrected data, and the rest as artifacts. The corrected
data and estimated OA are shown in Figure 10c.

From the results in Figure 10b and Figure 10c, $\operatorname{\text{ETEA}}$
estimates artifacts more clearly than MEMD. In the MEMD result, some small-
amplitude higher frequency oscillations leak into the artifact (e.g., about
$t=5.5$ and $9.0$ second). Moreover, some artifacts are introduced after
applying MEMD method to correct the data, (e.g., about $t=3.0$ and $9.5$
second in Figure 10c). Some oscillations are generated as transients where the
abrupt artifacts occur. In contrast, in Figure 10b, there are no oscillations
introduced either in the estimated artifacts or the corrected data.

Computational Efficiency. Figure 11 shows the average computation time as a
function of the signal length. In this experiment, we calculate the time of
computation of ETEA (22) and second-order ETEA (47) with different filter
settings (controlled by parameter $d$), using input signal lengths from 5000
to $10^{5}$. For each length, we average the computation time over 10 trials.
For each trial, run 40 iterations of each algorithm. The experiment is
implemented in Matlab 8.4 on a MacBook Pro 2012 with 2.7 GHz CPU.

As shown, the proposed algorithms have a run time of order $O(N)$. Most of the
computation time is consumed by the step of solving the linear system
$\mathbf{Q}^{-1}\mathbf{b}$ in Table 2, where $\mathbf{Q}$ is a banded matrix
and we use a fast banded solver.

Additionally, fast iterative shrinkage-thresholding algorithm (FISTA) [3, 7],
which is an acceleration scheme for iterative shrinkage-thresholding algorithm
(ISTA) [6], (a special formulation of MM) may be used to further accelerate
the algorithm.

Figure 11: Comparison of computation time with different signal length.

##  6 Conclusion

This paper proposes a new algorithm for denoising and artifact removal for
signals comprising artifacts arising in measured data, e.g., neural time-
series recordings, using sparse optimization. The first algorithm,
$\operatorname{\text{ETEA}}$, assumes the signal is composed of a lowpass
signal and an exponential transients (Type 1). It is formulated as an
optimization problem regularized by differentiable and smooth penalty
function. The second algorithm is an extension of
$\operatorname{\text{ETEA}}$, using a higher-order recursive filter, which is
applicable for correction of continuous protuberance transients (Type 0), and
more irregular artifacts. As applications, we have shown that
$\operatorname{\text{ETEA}}$ with different regularizers ($\mathbf{R}$ and
$\mathbf{R}_{2}$) are suitable for the suppression of Type 1 artifact in ECoG
data and ocular artifacts (OA) (as sequential Type 0 artifacts) in
conventional EEG data, with detailed comparisons to some state-of-the-art
methods. Both of the above algorithms are computationally efficient because
they are formulated in terms of banded matrices. A promising future work is to
extend the above data correcting methods to multichannel data.

## Acknowledgments

The authors would like to thank Jonathan Viventi of the Department of
Biomedical Engineering of Duke University, and Justin Blanco of the Electrical
and Computer Engineering Department of United States Naval Academy, for
providing the data and giving useful comments.

##  Appendix A Proof of Proposition 1

###### Proof.

Substitute the variables $u$ and $v$ in (14) by

$\displaystyle u=\sqrt{x^{2}+\epsilon},\text{ and }v=\sqrt{z^{2}+\epsilon}.$
(49)

Therefore, for $x,z\in\mathbb{R}$, the inequality holds:

$\displaystyle\frac{\phi^{\prime}(\sqrt{z^{2}+\epsilon})}{2\sqrt{z^{2}+\epsilon}}(x^{2}+\epsilon)+\phi(\sqrt{z^{2}+\epsilon})-\frac{\sqrt{z^{2}+\epsilon}}{2}\phi^{\prime}(\sqrt{z^{2}+\epsilon})\geqslant\phi(\sqrt{x^{2}+\epsilon}).$
(50)

The right of the inequality is the majorizer of the smoothed penalty function
$\phi(\sqrt{x^{2}+\epsilon})$.
1) If $z\neq 0$, we can multiply $z$ to the nominator and denominator of the
first term on the left side of (50),

$\displaystyle\frac{\phi^{\prime}(\sqrt{z^{2}+\epsilon})}{2\sqrt{z^{2}+\epsilon}}(x^{2}+\epsilon)=\frac{x^{2}}{2z}\bigg{(}\phi^{\prime}(\sqrt{z^{2}+\epsilon})\frac{z}{\sqrt{z^{2}+\epsilon}}\bigg{)}+\frac{\epsilon}{2z}\bigg{(}\phi^{\prime}(\sqrt{z^{2}+\epsilon})\frac{z}{\sqrt{z^{2}+\epsilon}}\bigg{)}=\frac{x^{2}}{2z}\phi_{\epsilon}^{\prime}(z)+\frac{\epsilon}{2z}\phi_{\epsilon}^{\prime}(z).$
(51)

Multiplying the nominator and denominator of the third term on the left side
of (50) by $\sqrt{z^{2}+\epsilon}$,

$\displaystyle\frac{\sqrt{z^{2}+\epsilon}}{2}\phi^{\prime}(\sqrt{z^{2}+\epsilon})=$
$\displaystyle\frac{z^{2}}{2\sqrt{z^{2}+\epsilon}}\phi^{\prime}(\sqrt{z^{2}+\epsilon})+\frac{\epsilon}{2\sqrt{z^{2}+\epsilon}}\phi^{\prime}(\sqrt{z^{2}+\epsilon})$
$\displaystyle=$
$\displaystyle\frac{z}{2}\bigg{(}\phi^{\prime}(\sqrt{z^{2}+\epsilon})\frac{z}{\sqrt{z^{2}+\epsilon}}\bigg{)}+\frac{\epsilon}{2z}\bigg{(}\phi^{\prime}(\sqrt{z^{2}+\epsilon})\frac{z}{\sqrt{z^{2}+\epsilon}}\bigg{)}$
$\displaystyle=$
$\displaystyle\frac{z}{2}\phi_{\epsilon}^{\prime}(z)+\frac{\epsilon}{2z}\phi_{\epsilon}^{\prime}(z).$
(52)

Using the results in (51) and (A), the inequality (50) can be rewritten as

$\displaystyle\frac{x^{2}}{2z}\phi_{\epsilon}^{\prime}(z)+\frac{\epsilon}{2z}\phi_{\epsilon}^{\prime}(z)+\phi_{\epsilon}(z)-\frac{z}{2}\phi_{\epsilon}^{\prime}(z)-\frac{\epsilon}{2z}\phi_{\epsilon}^{\prime}(z)\geqslant\phi_{\epsilon}(x),$
(53)

which can be reorganize into

$\displaystyle g_{\epsilon}(x,z;\phi_{\epsilon})$
$\displaystyle=\frac{\phi_{\epsilon}^{\prime}(z)}{2z}x^{2}+\phi_{\epsilon}(z)-\frac{z}{2}\phi_{\epsilon}^{\prime}(z)\geqslant\phi_{\epsilon}(x).$
(54)

2) If $z=0$ and $x\neq 0$, by Lagrange’s Mean Value Theorem [55, Theorem
5.10], since function $\phi(x)$ is continuous and
$\phi^{\prime\prime}(x)\leqslant 0$ on $x>0$, there exists a value $\xi$ in
the range $\sqrt{\epsilon}<\xi<\sqrt{x^{2}+\epsilon}$ that
$\phi^{\prime}(\sqrt{\epsilon})\geqslant\phi^{\prime}(\xi)\geqslant\phi^{\prime}(\sqrt{x^{2}+\epsilon})$
satisfying

$\displaystyle\phi^{\prime}(\xi)\left(\sqrt{x^{2}+\epsilon}-\sqrt{\epsilon}\right)=\phi(\sqrt{x^{2}+\epsilon})-\phi(\sqrt{\epsilon}).$
(55)

Moreover, consider the square that is always positive

$\displaystyle(\sqrt{x^{2}+\epsilon}-\sqrt{\epsilon})^{2}=x^{2}+\epsilon-2\sqrt{\epsilon(x^{2}+\epsilon)}+\epsilon>0,$
(56)

which implies the inequality

$\displaystyle x^{2}>2\sqrt{\epsilon(x^{2}+\epsilon)}-2\epsilon.$ (57)

Furthermore, we can multiply both sides of (57) by a positive term
$\displaystyle\frac{\phi^{\prime}(\sqrt{\epsilon})}{2\sqrt{\epsilon}}$, and
then adopt the result from (55), so that:

$\displaystyle\frac{\phi^{\prime}(\sqrt{\epsilon})}{2\sqrt{\epsilon}}x^{2}$
$\displaystyle>\frac{\phi^{\prime}(\sqrt{\epsilon})}{2\sqrt{\epsilon}}\left(2\sqrt{\epsilon(x^{2}+\epsilon)}-2\epsilon\right)$
(58a)
$\displaystyle=\phi^{\prime}(\sqrt{\epsilon})\left(\sqrt{x^{2}+\epsilon}-\sqrt{\epsilon}\right)$
(58b)
$\displaystyle>\phi^{\prime}(\xi)\left(\sqrt{x^{2}+\epsilon}-\sqrt{\epsilon}\right)$
(58c) $\displaystyle=\phi(\sqrt{x^{2}+\epsilon})-\phi(\sqrt{\epsilon}),$ (58d)

which leads to

$\displaystyle\frac{\phi^{\prime}(\sqrt{\epsilon})}{2\sqrt{\epsilon}}x^{2}+\phi(\sqrt{\epsilon})>\phi(\sqrt{x^{2}+\epsilon}).$
(59)

Because $\phi_{\epsilon}$ in (13) is differentiable on , we can find its
second-order derivative at zero,

$\displaystyle\lim_{z\to
0}\phi_{\epsilon}^{\prime\prime}(z)=\frac{\phi^{\prime}(\sqrt{\epsilon})}{\sqrt{\epsilon}},$
(60)

and by L’H$\hat{\textrm{o}}$pital’s rule [55, Theorem 5.13], we have

$\displaystyle\lim_{z\to 0}\frac{\phi_{\epsilon}^{\prime}(z)}{2z}=\lim_{z\to
0}\frac{\phi_{\epsilon}^{\prime\prime}(z)}{2}=\frac{\phi^{\prime}(\sqrt{\epsilon})}{2\sqrt{\epsilon}},$
(61)

which implies that (59) is in the same form of the majorizer (54) at $z=0$.
3) If $x=z=0$, then the condition (9) follows immediately. ∎

## References

  * [1] M. T. Akhtar, W. Mitsuhashi, and C. J. James.  Employing spatially constrained ICA and wavelet denoising, for automatic removal of artifacts from multichannel EEG data.  Signal Processing, 92(2):401–416, 2012.
  * [2] F. Bach, R. Jenatton, J. Mairal, and G. Obozinski.  Optimization with sparsity-inducing penalties.  Foundations and Trends in Machine Learning, 4(1):1–106, 2012.
  * [3] A. Beck and M. Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM J. Imag. Sci., 2(1):183–202, 2009.
  * [4] E. J. Candès, M. B. Wakin, and S. Boyd.  Enhancing sparsity by reweighted l1 minimization.  J. Fourier Anal. Appl., 14(5):877–905, December 2008.
  * [5] M. Cetin and W.C. Karl.  Feature-enhanced synthetic aperture radar image formation based on nonquadratic regularization.  IEEE Trans. Image Process., 10(4):623–631, April 2001.
  * [6] A. Chambolle, R. A. De Vore, N.-Y. Lee, and B. J. Lucier.  Nonlinear wavelet image processing: variational problems, compression, and noise removal through wavelet shrinkage.  IEEE Trans. Image Process., 7(3):319–335, March 1998.
  * [7] A. Chambolle and V. R. Dossal.  On the convergence of the iterates of FISTA.  hal-01060130, September 2014.  preprint.
  * [8] C.-P. Chang, J.-C. Lee, Y. Su, P. S. Huang, and T.-M. Tu.  Using empirical mode decomposition for iris recognition.  Computer Standards & Interfaces, 31(4):729–739, 2009.
  * [9] H.-C. Chang, P.-L. Lee, M.-T. Lo, Y.-T. Wu, K.-W. Wang, and G.-Y. Lan.  Inter-trial analysis of post-movement beta activities in EEG signals using multivariate empirical mode decomposition.  IEEE Trans. Neural Systems and Rehabilitation Engineering, 21(4):607–615, July 2013.
  * [10] 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.
  * [11] R. R. Coifman and D. L. Donoho.  Translation-invariant de-noising.  In A. Antoniadis, editor, Wavelets and Statistics. Springer-Verlag Lecture Notes, 1995.
  * [12] J. Dammers, M. Schiek, F. Boers, C. Silex, M. Zvyagintsev, U. Pietrzyk, and K. Mathiak.  Integration of amplitude and phase statistics for complete artifact removal in independent components of neuromagnetic recordings.  IEEE Trans. Biomed. Eng., 55(10):2353–2362, October 2008.
  * [13] S. Durand and J. Froment.  Artifact free signal denoising with wavelets.  In Proc. ICASSP, 2001.
  * [14] S. Durand and J. Froment.  Reconstruction of wavelet coefficients using total variation minimization.  SIAM J. Sci. Comput., 24(5):1754–1767, 2003.
  * [15] 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.
  * [16] J. Fleureau, A. Kachenoura, L. Albera, J.-C. Nunes, and L. Senhadji.  Multivariate empirical mode decomposition and application to multichannel filtering.  Signal Processing, 91(12):2783–2792, 2011.
  * [17] J.-J. Fuchs.  On sparse representations in arbitrary redundant bases.  IEEE Trans. Inform. Theory, 50(6):1341–1344, 2004.
  * [18] H. Gao.  Wavelet shrinkage denoising using the non-negative garrote.  Journal of Computational and Graphical Statistics, 7(4):pp. 469–488, 1998.
  * [19] A. Gholami and S. M. Hosseini.  A balanced combination of Tikhonov and total variation regularizations for reconstruction of piecewise-smooth signals.  Signal Processing, 93(7):1945–1960, 2013.
  * [20] C. Guerrero-Mosquera and A. Navia-Vázquez.  Automatic removal of ocular artefacts using adaptive filtering and independent component analysis for electroencephalogram data.  IET Signal Processing, 6(2):99–106, 2012.
  * [21] H. Huang and J. Pan.  Speech pitch determination based on Hilbert-Huang transform.  Signal Processing, 86(4):792–803, April 2006.
  * [22] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N. C. Yen, C. C. Tung, and H. H. Liu.  The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis.  Proc. Roy. Soc. Lon. A, 454:903–995, March 1998.
  * [23] N. E. Huang, M.-L. C. Wu, S. R. Long, S. S. P. Shen, W. Qu, P. Gloersen, and K. L. Fan.  A confidence limit for the empirical mode decomposition and Hilbert spectral analysis.  Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 459(2037):2317–2345, September 2003.
  * [24] D. R. Hunter and K. Lange.  A tutorial on MM algorithms.  Amer. Statist., 58:30–37, 2004.
  * [25] M. K. Islam, A. Rastegarnia, A. T. Nguyen, and Z. Yang.  Artifact characterization and removal for in vivo neural recording.  Journal of Neuroscience Methods, 226(0):110–123, 2014.
  * [26] C. A. Joyce, I. F. Gorodnitsky, and M. Kutas.  Automatic removal of eye movement and blink artifacts from EEG data using blind component separation.  Psychophysiology, 41(2):313–325, 2004.
  * [27] F. I. Karahanoglu, I. Bayram, and D. Van De Ville.  A signal processing approach to generalized 1-d total variation.  IEEE Trans. Signal Process., 59(11):5265–5274, November 2011.
  * [28] E. Kilic.  Explicit formula for the inverse of a tridiagonal matrix by backward continued fractions.  Applied Mathematics and Computation, 197(1):345–357, 2008.
  * [29] E. Kilic and P. Stanica.  The inverse of banded matrices.  Journal of Computational and Applied Mathematics, 237(1):126–135, 2013.
  * [30] S. Kim, K. Koh, S. Boyd, and D. Gorinevsky.  $\ell_{1}$ trend filtering.  SIAM Review, 51(2):339–360, 2009.
  * [31] Y. Kopsinis and S. McLaughlin.  Empirical mode decomposition based soft-thresholding.  In 16th European Signal Processing Conference, 2008.
  * [32] Y. Kopsinis and S. McLaughlin.  Development of EMD-based denoising methods inspired by wavelet thresholding.  IEEE Trans. Signal Process., 57(4):1351–1362, April 2009.
  * [33] M. Kowalski.  Sparse regression using mixed norms.  Applied and Computational Harmonic Analysis, 27(3):303 – 324, 2009\.
  * [34] I. Kozlov and A. Petukhov.  Sparse solutions of underdetermined linear systems.  In W. Freeden et al., editor, Handbook of Geomathematics. Springer, 2010.
  * [35] K. Lange.  Optimization.  Springer New York, 2004.
  * [36] K. Lange, D. Hunter, and I. Yang.  Optimization transfer using surrogate objective functions.  J. of Comp. Graph. Statist., 9:1–20, 2000.
  * [37] N. Mammone, F. La Foresta, and F. C. Morabito.  Automatic artifact rejection from multichannel scalp EEG by wavelet ICA.  IEEE J. Sensors, 12(3):533–542, March 2012.
  * [38] D. Mandic.  Empirical mode decomposition, multivariate EMD, matlab and data sources.
  * [39] D. P. Mandic, N. U. Rehman, Z. Wu, and N. E. Huang.  Empirical mode decomposition-based time-frequency analysis of multivariate signals: The power of adaptive data analysis.  Signal Processing Magazine, IEEE, 30(6):74–86, November 2013.
  * [40] B. Molavi and G. A. Dumont.  Wavelet-based motion artifact removal for functional near-infrared spectroscopy.  Physiological Measurement, 33(2):259, 2012.
  * [41] M. K. I. Molla, M. R. Islam, T. Tanaka, and T. M. Rutkowski.  Artifact suppression from EEG signals using data adaptive time domain filtering.  Neurocomputing, 97:297–308, 2012.
  * [42] M. K. I. Molla, T. Tanaka, and T. M. Rutkowski.  Multivariate EMD based approach to EOG artifacts separation from EEG.  In Proc. ICASSP 2012, pages 653–656, 2012.
  * [43] M. K. I. Molla, T. Tanaka, T. M. Rutkowski, and A. Cichocki.  Separation of EOG artifacts from EEG signals using bivariate EMD.  In Proc. ICASSP 2010, pages 562–565, 2010.
  * [44] K. Nazarpour, Y. Wongsawat, S. Sanei, J. A. Chambers, and S. Oraintara.  Removal of the eye-blink artifacts from EEGs via STF-TS modeling and robust minimum variance beamforming.  IEEE Trans. Biomed. Eng., 55(9):2221–2231, 2008.
  * [45] D. Needell.  Noisy signal recovery via iterative reweighted l1-minimization.  In Proc. Forty-Third Asilomar Conference on Signals, Systems and Computers, pages 113–117, 2009.
  * [46] M. Nikolova.  Analysis of the recovery of edges in images and signals by minimizing nonconvex regularized least-squares.  Multiscale Modeling and Simulation, 4(3):960–991, 2005.
  * [47] 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.
  * [48] B. Noureddin, P. D. Lawrence, and G. E. Birch.  Online removal of eye movement and blink EEG artifacts using a high-speed eye tracker.  IEEE Trans. Biomed. Eng., 59(8):2103–2110, 2012.
  * [49] A. Omidvarnia, G. Azemi, P. B. Colditz, and B. Boashash.  A time-frequency based approach for generalized phase synchrony assessment in nonstationary multivariate signals.  Digital Signal Processing, 23(3):780–790, 2013.
  * [50] C. Park, M. Plank, J. Snider, S. Kim, H. C. Huang, S. Gepshtein, T. P. Coleman, and H. Poizner.  EEG gamma band oscillations differentiate the planning of spatially directed movements of the arm versus eye: Multivariate empirical mode decomposition analysis.  IEEE Trans. Neural Systems and Rehabilitation Engineering, 22(5):1083–1096, September 2014.
  * [51] R. M. Rangayyan.  Biomedical Signal Analysis - A Case-study Approach.  IEEE and Wiley, New York, NY, 2002.
  * [52] N. U. Rehman and D. P. Mandic.  Multivariate empirical mode decomposition.  Proceedings of the Royal Society A, 466(2117):1291–1302, 2010.
  * [53] N. U. Rehman and D. P. Mandic.  Filter bank property of multivariate empirical mode decomposition.  IEEE Trans. Signal Process., 59(5):2421–2426, May 2011.
  * [54] G. Rilling and P. Flandrin.  One or two frequencies? The empirical mode decomposition answers.  IEEE Trans. Signal Process., 56(1):85–95, 2008.
  * [55] W. Rudin.  Principles of mathematical analysis. McGraw-Hill, 1976.
  * [56] H. Sato, N. Tanaka, M. Uchida, Y. Hirabayashi, M. Kanai, T. Ashida, I. Konishi, and A. Maki.  Wavelet analysis for detecting body-movement artifacts in optical topography signals.  NeuroImage, 33(2):580–587, 2006.
  * [57] E. D. Schifano, R. L. Strawderman, and M. T. Wells.  Majorization-minimization algorithms for nonsmoothly penalized objective functions.  Electron. J. Statist., 4:1258–1299, 2010.
  * [58] I. W. Selesnick, S. Arnold, and V. Dantham.  Polynomial smoothing of time series with additive step discontinuities.  IEEE Trans. Signal Process., 60(12):6305–6318, December 2012.
  * [59] I. W. Selesnick and I. Bayram.  Sparse signal estimation by maximally sparse convex optimization.  IEEE Trans. Signal Process., 62(5):1078–1092, March 2014.
  * [60] 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.
  * [61] I. W. Selesnick, H. L. Graber, S. Douglas, S. Pfeil, and R. L. Barbour.  Simultaneous low-pass filtering and total variation denoising.  IEEE Trans. Signal Process., 62(5):1109–1124, March 2014.
  * [62] B. Tang, S. Dong, and T. Song.  Method for eliminating mode mixing of empirical mode decomposition based on the revised blind source separation.  Signal Processing, 92(1):248–258, 2012.
  * [63] D. Wipf and S. Nagarajan.  Iterative reweighted $\ell_{1}$ and $\ell_{2}$ methods for finding sparse solutions.  IEEE. J. Sel. Top. Signal Processing, 4(2):317–329, April 2010\.
  * [64] Y. Wongsawat.  Efficient implementation of RMVB for eyeblink artifacts removal of EEG via STF-TS modeling.  In Proc. ROBIO 2008, pages 1567–1572, 2008.
  * [65] X. Xie.  Illumination preprocessing for face images based on empirical mode decomposition.  Signal Processing, 103(0):250–257, 2014.
  * [66] M. Yaghoobi, T. Blumensath, and M. E. Davies.  Dictionary learning for sparse approximations with the majorization method.  IEEE Trans. Signal Process., 57(6):2178–2191, June 2009.
  * [67] J. Yan and L. Lu.  Improved Hilbert–Huang transform based weak signal detection methodology and its application on incipient fault diagnosis and ECG signal analysis.  Signal Processing, 98:74–87, 2014.
  * [68] S. Yu, A. S. Khwaja, and J. Ma.  Compressed sensing of complex-valued data.  Signal Processing, 92(2):357–362, 2012.
  * [69] H. Zeng, A. Song, R. Yan, and H. Qin.  EOG artifact correction from EEG recording using stationary subspace analysis and empirical mode decomposition.  Sensors, 13(11):14839–14859, 2013.