\frac{a'' - a \cdot {\left( \phi' \right)}^2}{{(\phi')}^2} \cdot \cos(\phi) +\\
&\hspace{0.45cm} \left[ \frac{-2 \cdot a'}{a \cdot {\left( \phi' \right)}^2} -
\frac{\phi''}{{(\phi')}^3} \right] \cdot
{\left( a' \cdot \cos(\phi) - a \cdot \phi' \cdot \sin(\phi) \right)} +\\
&\hspace{0.45cm} \left[ \frac{2 \cdot {(a')}^2 - a \cdot a''}
{a^2 \cdot {(\phi')}^2} + \frac{a' \cdot \phi''}{a \cdot {(\phi')}^3} + 1 \right] \cdot
{\left( a \cdot \cos(\phi) \right)}.
\end{align*}

Sine and cosine are separated and we show that their coefficients are zero, implying
that the entire expression is zero, as follows:

\begin{align*}
\mathcal{D}_{(a,\phi)} \I[a,\phi] &=
\left[ \frac{-2 \cdot a'}{\phi'} - \frac{a \cdot \phi''}{{(\phi')}^2}\right] \cdot
\sin(\phi) + \left[ \frac{a''}{{(\phi')}^2} - a \right] \cdot \cos(\phi) +\\
&\hspace{0.45cm} \left[ \frac{2 \cdot a'}{\phi'} + \frac{a \cdot \phi''}
{{(\phi')}^2} \right] \cdot \sin(\phi) +
\left[ \frac{-2 \cdot {\left(a'\right)}^2}{a \cdot {\left( \phi' \right)}^2} -
\frac{a' \cdot \phi''}{{(\phi')}^3} \right] \cdot \cos(\phi) +\\
&\hspace{0.45cm} \left[ \frac{2 \cdot {(a')}^2}{a \cdot {(\phi')}^2} - \frac{a''}
{{(\phi')}^2} + \frac{a' \cdot \phi''}{{(\phi')}^3} + a \right] \cdot \cos(\phi) \\
&= \left[ \frac{-2 \cdot a'}{\phi'} - \frac{a \cdot \phi''}{{(\phi')}^2} +
\frac{2 \cdot a'}{\phi'} + \frac{a \cdot \phi''}{{(\phi')}^2} \right] \cdot \sin(\phi) +\\
&\hspace{0.45cm} \left[ \frac{a''}{{(\phi')}^2} - a +
\frac{-2 \cdot {\left(a'\right)}^2}{a \cdot {\left( \phi' \right)}^2} -
\frac{a' \cdot \phi''}{{(\phi')}^3} +
\frac{2 \cdot {(a')}^2}{a \cdot {(\phi')}^2} - \frac{a''}
{{(\phi')}^2} + \frac{a' \cdot \phi''}{{(\phi')}^3} + a \right] \cdot \cos(\phi)\\
&= 0. \qedhere
\end{align*}

This shows that $\mathcal{D}_{(a,\phi)}$ is in fact an annihilating operator. When considering
the IMF differential operator from Definition <ref> again,
one can observe that the only way the
$a$ and $\phi$ interface with the operator is through
instantaneous envelope derivation operator $A[a]$
and inverse square continuous frequency operator $\Omega[\phi]$.
Given both are functions just like $a$ and $\phi$, we can
express the operator parametrized by $A$ and $\Omega$ instead of $a$
and $\phi$ and call it the modified IMF operator.

The modified IMF operator is defined as

\begin{align*}
\mathcal{\tilde{D}}_{(A,\Omega)} & := \Omega \cdot D^2 +\\
&\hspace{0.5cm} \left[ -2 \cdot \Omega \cdot A +
\frac{1}{2} \cdot \Omega' \right]
\cdot D^1 +\\
&\hspace{0.5cm} \left[ \Omega \cdot \left( {A}^2 - A' \right) -
\frac{1}{2} \cdot \Omega' \cdot A + 1 \right] \cdot D^0.
\end{align*}

How to determine $a$ from $A$ and $\phi$ from $\Omega$ shall not yet
be of concern here, but for instance in the case of $\Omega$,
the inverse square root of $\Omega$ yields $\phi'$ directly
(compare Definition <ref>).
It is now in our interest to examine the behaviour of the differential
operator under its parameters. As we vary $A$ and $\Omega$, we want to
know that if we found an annihilating pair we really obtained a unique
solution or not. We will approach this question just like the
cost functions in Section <ref> and consider the
$A$ and $\Omega$ to be spline functions that vary over their
B-spline basis coefficients $\boldsymbol{A} \in \R^n$ and
$\boldsymbol{\Omega} \in \R^n$. To give an example, we consider $A$
to be the spline function (using Definitions <ref> and

\begin{equation*}
\B_{k}(\boldsymbol{A}) = \sum_{i=0}^{n-1} A_i \cdot B_{i,k}
\end{equation*}

that varies over the vector entries $A_i$ of $\boldsymbol{A}$.
If we consider the norm of this function, we can for instance
calculate its partial derivative in $A_m$ for $m \in \{0,\dots,n-1\}$
using the chain rule as

\begin{equation*}
\frac{\partial {\| \B_{k}(\boldsymbol{A}) \|}_2^2}{\partial A_m} =
\frac{\partial}{\partial A_m} \left(
\int_{-\infty}^{\infty} {\left( \B_{k}(\boldsymbol{A}) \right)}^2 \mathrm{d}t
\right) =
\int_{-\infty}^{\infty}
\frac{\partial {\left( \B_{k}(\boldsymbol{A}) \right)}^2}
{\partial A_m} \mathrm{d}t =
\int_{-\infty}^{\infty}
2 \cdot \B_{k}(\boldsymbol{A}) \cdot B_{m,k}(t) \mathrm{d}t.
\end{equation*}

To go even further, we can of course also partially derive again,
this time in $A_p$ for $p \in \{0,\dots,n-1\}$. We obtain

\begin{equation*}
\frac{\partial^2 {\| \B_{k}(\boldsymbol{A}) \|}_2^2}
{\partial A_m \partial A_p} =
\frac{\partial}{\partial A_m} \left(
\int_{-\infty}^{\infty}
2 \cdot \B_{k}(\boldsymbol{A}) \cdot B_{m,k}(t) \mathrm{d}t
\right) =
\int_{-\infty}^{\infty}
2 \cdot B_{m,k}(t) \cdot B_{p,k}(t) \mathrm{d}t.
\end{equation*}

The result of this particular observation is that the covariation
of ${\| \B_{k}(\boldsymbol{A}) \|}_2^2$ in $A_m$ and $A_p$ is directly
related to the orthogonality of $B_{m,k}$ and $B_{p,k}$. If the
B-spline basis were truly orthogonal, the final integral will always
be zero. However, the B-spline basis is not orthogonal. Thus, if
$m$ and $p$ are close to each other or even equal, the
integral will be positive. Another interpretation is to consider
the Hessian matrix in partial derivates in $A_i$, which can be
used to prove that a function is convex in multiple variables.
After all, the function ${\| \B_{k}(\boldsymbol{A}) \|}_2^2$ is
a mapping $\R^n \to \R$. Making use of this general technique,
we prove the following

Let $f \in \mathcal{C}^2(\R,\R)$ and $k \ge 4$ (for derivability).
${\left\| \tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}
f \right\|}_2^2$ is strictly convex in $\boldsymbol{\Omega}$ but not
convex in $\boldsymbol{A}$.

By definition we obtain that

\begin{equation*}
{\left\| \tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}
f \right\|}_2^2 =
\int_{-\infty}^{\infty}
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
\right)}^2
\mathrm{d}t.
\end{equation*}

In particular we first consider the term within the norm itself and expand it
from its definition (see Definition <ref>)

\begin{align*}
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))} f &=
\B_{k}(\boldsymbol{\Omega}) \cdot f'' +\\
&\hspace{0.5cm} \bigg[ -2 \cdot \B_{k}(\boldsymbol{\Omega}) \cdot \B_{k}(\boldsymbol{A}) +
\frac{1}{2} \cdot {\left(\B_{k}(\boldsymbol{\Omega})\right)}' \bigg] \cdot f' +\\
&\hspace{0.5cm} \bigg[ \B_{k}(\boldsymbol{\Omega}) \cdot \left(
{\left(\B_{k}(\boldsymbol{A})\right)}^2 - {\left(\B_{k}(\boldsymbol{A})\right)}' \right) -
\frac{1}{2} \cdot {\left( \B_{k}(\boldsymbol{\Omega}) \right)}' \cdot \B_{k}(\boldsymbol{A}) +
1 \bigg] \cdot f.
\end{align*}

Making use of Definition <ref> to expand
the spline functions into B-spline expressions it follows

\begin{align*}
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))} f
&= \left(\sum_{i=0}^{n-1} \Omega_i \cdot B_{i,k}\right) \cdot f'' +\\
&\hspace{0.5cm} \Bigg[ -2 \cdot
\left(\sum_{i=0}^{n-1} \Omega_i \cdot B_{i,k}\right) \cdot
{\left(\sum_{i=0}^{n-1} A_i \cdot B_{i,k}\right)} +\\
\frac{1}{2} \cdot \left(\sum_{i=0}^{n-1} \Omega_i \cdot
B'_{i,k}\right) \Bigg] \cdot f' +\\
&\hspace{0.5cm} \Bigg[ \left(\sum_{i=0}^{n-1}
\Omega_i \cdot B_{i,k}\right) \cdot \Bigg(
{\left(\sum_{i=0}^{n-1} A_i \cdot B_{i,k}\right)}^2 -
{\left(\sum_{i=0}^{n-1} A_i \cdot B'_{i,k}\right)} \Bigg) -\\
\frac{1}{2} \cdot \left(\sum_{i=0}^{n-1} \Omega_i \cdot B'_{i,k}\right)
\cdot \left(\sum_{i=0}^{n-1} A_i \cdot B_{i,k}\right) +
1 \Bigg] \cdot f.
\end{align*}

We approach this proof checking if the requirements of
Theorem <ref> hold for distinct partial
derivatives for entries of $\boldsymbol{A}$ and $\boldsymbol{\Omega}$.
By applying partial derivatives we obtain that for $m,p \in \{ 0,\dots,n-1 \}$
it holds with the chain rule for partial derivatives in $\boldsymbol{A}$

\begin{align}
\frac{\partial^2 {\left\|\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),
\B_{k}(\boldsymbol{\Omega}))} f\right\|}_2^2}
{\partial A_m \partial A_p}(\boldsymbol{A},\boldsymbol{\Omega}) &=
\frac{\partial}{\partial A_p} \left(
\int_{-\infty}^{\infty}
\cdot
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
\cdot
\frac{
\partial
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
}{\partial A_m}
\mathrm{d}t
\right) \notag \\
&= \int_{-\infty}^{\infty}
\cdot
\frac{
\partial
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
}{\partial A_m}
\cdot
\frac{
\partial
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
}{\partial A_p} + \notag \\
\cdot
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
\cdot
\frac{
\partial^2
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
}{\partial A_m \partial A_p}
\mathrm{d}t \label{eq:proof-diffop-convex-A}
\end{align}

and analogously for partial derivatives in $\boldsymbol{\Omega}$

\begin{align}
\frac{\partial^2 {\left\|\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),
\B_{k}(\boldsymbol{\Omega}))} f\right\|}_2^2}
{\partial \Omega_m \partial \Omega_p}(\boldsymbol{A},\boldsymbol{\Omega}) &=
\int_{-\infty}^{\infty}
\cdot
\frac{
\partial
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
}{\partial \Omega_m}
\cdot
\frac{
\partial
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
}{\partial \Omega_p} + \notag \\
\cdot
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
\cdot
\frac{
\partial^2
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
}{\partial \Omega_m \partial \Omega_p}
\mathrm{d}t. \label{eq:proof-diffop-convex-Omega}
\end{align}

First for partial derivatives in $\boldsymbol{A}$, we consider the bare
derivatives of the operator (without the norm) that we found within the derivative
expressions of the operator within the norm. For the first order we find that

\begin{align*}
\frac{\partial \tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))} f }
{\partial A_m}(\boldsymbol{A},\boldsymbol{\Omega}) &= \left[ -2 \cdot
{\B_{k}(\boldsymbol{\Omega})} \cdot B_{m,k} \right] \cdot f' +\\
&\hspace{0.41cm} \left[ \B_{k}(\boldsymbol{\Omega}) \cdot \left(
2 \cdot \B_{k}(\boldsymbol{A}) \cdot B_{m,k} - B'_{m,k} \right) -
\frac{1}{2} \cdot {\left( \B_{k}(\boldsymbol{\Omega}) \right)}'
\cdot B_{m,k} \right] \cdot f\\
&= \left[ -2 \cdot
{\B_{k}(\boldsymbol{\Omega})} \cdot B_{m,k} \right] \cdot f' +\\
&\hspace{0.46cm} \Bigg[ \B_{k}(\boldsymbol{\Omega}) \cdot \left(
2 \cdot {\left(\sum_{i=0}^{n-1} A_i \cdot B_{i,k}\right)} \cdot
B_{m,k} - B'_{m,k} \right) -\\
&\hspace{0.7cm} \frac{1}{2} \cdot {\left( \B_{k}(\boldsymbol{\Omega}) \right)}'
\cdot B_{m,k} \Bigg] \cdot f,
\end{align*}

and for the second order we finally obtain

\begin{equation*}
\frac{\partial^2 \tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))} f }
{\partial A_m \partial A_p}(\boldsymbol{A},\boldsymbol{\Omega}) =
2 \cdot \B_{k}(\boldsymbol{\Omega}) \cdot B_{m,k} \cdot B_{p,k} \cdot f \not\equiv 0.
\end{equation*}

This implies that for $m = p$ the sign of the term in Equation (<ref>)
is not positive, as the bare operator is not zeroed out in the second summand.
Thus the necessary condition in Theorem <ref>
that all diagonal elements of the Hessian matrix must be positive is violated. It follows
that we can not show convexity for partial derivatives in $\boldsymbol{A}$, as the Hessian
matrix for derivatives in $\boldsymbol{A}$ can not be shown to be positive semidefinite.
For partial derivatives of the bare operator in $\boldsymbol{\Omega}$ we find
for the first order that

\begin{align*}
\frac{\partial \tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))} f }
{\partial \Omega_m}(\boldsymbol{A},\boldsymbol{\Omega}) &=
B_{m,k} \cdot f'' +\\
&\hspace{0.5cm} \bigg[ -2 \cdot \B_{k}(\boldsymbol{A}) \cdot B_{m,k} +
\frac{1}{2} \cdot B'_{m,k} \bigg] \cdot f' +\\
&\hspace{0.5cm} \bigg[ B_{m,k} \cdot \left(
{\left(\B_{k}(\boldsymbol{A})\right)}^2 - {\left(\B_{k}(\boldsymbol{A})\right)}' \right) -
\frac{1}{2} \cdot \B_{k}(\boldsymbol{A}) \cdot B'_{m,k} \bigg] \cdot f
\end{align*}

which yields for the second order that

\begin{equation*}
\frac{\partial^2 \tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))} f }
{\partial \Omega_m \partial \Omega_p}(\boldsymbol{A},\boldsymbol{\Omega}) \equiv 0.
\end{equation*}

We obtain from this result that the second summand in the integral in
Equation (<ref>) is zero and it holds that

\begin{equation*}
\frac{\partial^2 {\left\|\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),
\B_{k}(\boldsymbol{\Omega}))} f\right\|}_2^2}
{\partial \Omega_m \partial \Omega_p}(\boldsymbol{A},\boldsymbol{\Omega}) =
\int_{-\infty}^{\infty}
\cdot
\frac{
\partial
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
}{\partial \Omega_m}
\cdot
\frac{
\partial
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
}{\partial \Omega_p}
\mathrm{d}t.
\end{equation*}

We can immediately see that the Hessian matrix of
\B_{k}(\boldsymbol{\Omega}))} f\right\|}_2^2$ for derivatives in $\boldsymbol{\Omega}$
is symmetric. In particular, for $m = p$, it also holds

\begin{equation*}
\frac{\partial^2 {\left\|\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),
\B_{k}(\boldsymbol{\Omega}))} f\right\|}_2^2}
{\partial \Omega_m^2}(\boldsymbol{A},\boldsymbol{\Omega}) =
\int_{-\infty}^{\infty}
\cdot
\frac{
\partial
\tilde{\mathcal{D}}_{(\B_{k}(\boldsymbol{A}),\B_{k}(\boldsymbol{\Omega}))}f
}{\partial \Omega_m}
\right)}^2
\mathrm{d}t > 0,
\end{equation*}

which means that the diagonal entries of the Hessian matrix for
$\boldsymbol{\Omega}$ are strictly positive. It is also diagonally dominant
with the same argument as in the proof of Proposition <ref>,
namely due to the compact support and partial orthogonality of the
B-spline basis functions. It follows with Theorem <ref>
that the Hessian matrix for derivatives in $\boldsymbol{\Omega}$ is positive definite
and by Proposition <ref> that
\B_{k}(\boldsymbol{\Omega}))} f\right\|}_2^2$ is strictly convex in $\boldsymbol{\Omega}$.

Given this result we do not have the theoretical guarantee that our
operator gives us a unique and minimal solution, as it is not generally
convex. This is obviously undesirable as our primary motivation is to
find and examine methods that have a stronger theoretical foundation
than the classic heuristic EMD methods.

§.§ Simplification

The reassuring part of the result in Theorem <ref>
is that the function is strictly convex if we reduce the variation to
$\boldsymbol{\Omega}$ and keep $\boldsymbol{A}$ constant. Without loss of
generality, if we know that our input IMF has constant amplitude $1$,
namely that it only has the form $\cos(\phi)$, we can apply a simplified differential
operator to it of which we know that it is strictly convex. This assumption may
sound a bit too extravagant, but we will show in Chapter <ref>
that it is meaningful and use the results of the following subsection
to extract the instantaneous phase from IMFs with constant amplitude $1$.
If we know that our input IMFs will have the form $\cos(\phi)$, we might
wonder how our IMF differential operator changes under this assumption.

It holds

\begin{equation*}
\mathcal{D}_{(1,\phi)} = \Omega[\phi] \cdot D^2 +
\frac{1}{2} \cdot \Omega'[\phi] \cdot D^1 +
\end{equation*}

It follows directly from Definition <ref> and observing that

\begin{equation*}
A[1] = \frac{{(1)}'}{1} = 0.\qedhere
\end{equation*}

The great simplification of the operator is apparent. We now wonder how
the modified IMF operator behaves under the assumption that $a \equiv 1$.
Making the observation that the instantaneous envelope derivation operator
$A[a]$ (see Definition <ref>) vanishes
for $a \equiv 1$, as $A[1] = {(1)}' / 1 = 0$, we can
see that the parametrization for the modified IMF operator is
$(0,\Omega)$ and we can formulate the following

Let $f \in \mathcal{C}^2(\R,\R)$ and $k \ge 4$ (for derivability).
${\left\| \tilde{\mathcal{D}}_{(0,\B_{k}(\boldsymbol{\Omega}))}
f \right\|}_2^2$ is strictly convex in $\boldsymbol{\Omega}$.

This follows directly from Theorem <ref>.

Given this convexity property, we have a theoretical guarantee that we reach
a global minimum for a given input IMF and a unique $\boldsymbol{\Omega}$.
As already mentioned earlier, we obtain the desired instantaneous frequency
from $\boldsymbol{\Omega}$ by inversely applying the inverse square
continuous frequency operator in Definition <ref>.
This equates to inverting and taking the square root of $\Omega$, which
is a relatively simple operation.

§.§ Discretization

Given the results from this chapter and especially Theorem <ref>
we will as follows only consider the simple case with constant amplitude $a \equiv 1$
for discretization, as this will also be the only relevant case for the toolbox
presented in Chapter <ref> given we can't use the differential operator
to extract the amplitude anyway.
We have shown in Proposition <ref> that
${\left\| \tilde{\mathcal{D}}_{(0,\B_{k}(\boldsymbol{\Omega}))}
f \right\|}_2^2$ is strictly convex in $\boldsymbol{\Omega}$, however,
for the discretization we can make two observations to simplify it:
The first is that given $\boldsymbol{\Omega} \in \R^n$ we need a system of
at least $n$ samples and additional boundary conditions to solve the problem.
The second is that given we have a uniform grid we can,
instead of minimizing an integral function, minimize at least $n$ samples
f\right)(t_i)$ over $\boldsymbol{\Omega}$. We take this detour as
we can see that
$\tilde{\mathcal{D}}_{(0,\B_{k}(\boldsymbol{\Omega}))} f$ is linear
in $\boldsymbol{\Omega}$, yielding a least squares problem of
the form

\begin{equation*}
\left\| A \cdot \boldsymbol{\Omega} - b \right\|^2,
\end{equation*}

where $A$ is the matrix representing $\tilde{\mathcal{D}}_
{(0,\B_{k}(\boldsymbol{\Omega}))}$ and $b$ is the vector of
samples in $t_i$. This is better than the nonlinear problem
we would obtain by just using the squared integral.
Given the precomputation-concepts of the toolbox we make use of the
precomputed extended grid (see Section <ref>)
and evaluate $\tilde{\mathcal{D}}_{(0,\B_{k}(\boldsymbol{\Omega}))} f$ on the extended
grid, obtaining more than $n$ equations, one for each point on the extended
grid. This way we obtain implicit boundary conditions, saving us from proposing
possibly wrong or ill-chosen ones in the process. See Subsection <ref>
for more reflexions on boundary effects.
The instantaneous frequency is calculated from $\boldsymbol{\Omega}$
by evaluating $\B_k(\boldsymbol{\Omega})$, applying an inversion and
square root and running a B-splines-fit on the resulting data. Given
the nature of the transformation the inverse square
continuous frequency operator in
Definition <ref> specifies it is
most likely impossible to exploit any B-spline property to circumvent this
step and directly work with the B-spline coefficient vector $\boldsymbol{\Omega}$.

§ EXAMPLES

Following the previous theoretical perspective, this section gives
a few examples on the numerical behaviour of the simplified
IMF differential operator. For this purpose we restrict ourselves
to IMFs with constant instantaneous amplitude $a \equiv 1$ and known
analytical form and apply our operator to them. The question is
how well we manage to extract the instantaneous phase $\phi$,
which we can assess by comparing the results to the ground truth.
The examples were implemented using the ETHOS-toolbox and can be found
in Listing <ref>.
The parameters $(k,q,n)$ given in the figure captions refer to the
spline order $k$, in-fill-count $q$ (see Section <ref>)
and number of B-spline basis functions $n$. See Subsection <ref>
for a discussion on the boundary effects of these examples in the context
of information theory and other literature.
[Constant frequency]
As an introduction consider the simple IMF

\begin{equation}\label{eq:example_diffop-0-u}
u_0(t) := \cos(\phi_0(t)) := \cos(40 \cdot t),
\end{equation}

on the interval $[0,1]$ (see Figure <ref>).
The analytical instantaneous frequency $\phi'_0$ is
$40$, i.e. the IMF is of constant frequency.
Given the instantaneous amplitude is constantly $1$, we
can use the toolbox to fit the simple IMF differential operator to $u_0$ to calculate the
numerical instantaneous frequency $\tilde{\phi}'_0$.
The difference between $\phi'_0$ and $\tilde{\phi}'_0$
is too small to be visible in a normal plot and thus we examine the semi-log plot of the
relative error (see Figure <ref>). We can see that the
relative error is at most $0.1\%$ briefly at the beginning and stays below
$10^{-4}$ on the remaining interval.

[xlabel=$t$, ylabel=$u_0(t)$]
[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/regop.data/0-u.csv;

Plot of the IMF $u_0$ (see (<ref>)) from
Example <ref>.

[xlabel=$t$, ylabel near ticks,
ylabel=$\left| \frac{\phi'_0(t) -
\tilde{\phi}'_0(t)}{\phi'_0(t)} \right|$]
[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/regop.data/0-relerr.csv;

Semi-log plot of the relative error between the analytical
instantaneous frequency $\phi'_0=40$ and the solution
$\tilde{\phi}'_0$, obtained using the
IMF differential operator for $(k, q, n) = (4, 4, 180)$, from
Example <ref>.

[Harmonic peaks]
Consider the simple IMF

\begin{equation}\label{eq:example_diffop-1-u}
u_1(t) := \cos(\phi_1(t)) := \cos(3 \cdot \sin(3 \cdot \pi \cdot t) +
16 \cdot \pi \cdot t),
\end{equation}

on the interval $[0,1]$ (see Figure <ref>).
It is easy to analytically obtain the instantaneous
frequency $\phi'_1$ of the signal by calculation, namely

\begin{equation}\label{eq:example_diffop-1-freq-analytic}
\phi'_1(t) = 9 \cdot \pi \cdot \cos(3 \cdot \pi \cdot t) +
16 \cdot \pi,
\end{equation}

which you can find pictured in Figure <ref>.
Given the IMF has constant instantaneous amplitude $1$ we can use
the toolbox to fit the simple IMF differential operator to $u_1$ to calculate the numerical
instantaneous frequency $\tilde{\phi}'_1$. As can be seen the instantenous frequency
itself is a wave function too, explaining the irregular shape of the IMF.
The difference between $\phi'_1$ and $\tilde{\phi}'_1$
is too small to be visible in a normal plot and we thus examine the semi-log plot of the
relative error (see Figure <ref>). We can see that the
relative error is at most $1\%$ briefly at the beginning and between $1$ and $4$ orders
of magnitude lower on the remaining interval.

[xlabel=$t$, ylabel=$u_1(t)$]
[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/regop.data/1-u.csv;

Plot of the IMF $u_1$ (see (<ref>)) from
Example <ref>.

[xlabel=$t$, ylabel=$\phi'_1(t)$]
[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/regop.data/1-freq-analytic.csv;

Plot of the analytical instantaneous frequency
$\phi'_1$ (see (<ref>))
from Example <ref>.

[xlabel=$t$, ylabel near ticks,
ylabel=$\left| \frac{\phi'_1(t) -
\tilde{\phi}'_1(t)}{\phi'_1(t)} \right|$]
[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/regop.data/1-relerr.csv;

Semi-log plot of the relative error between the analytical
instantaneous frequency $\phi'_1$
(see (<ref>)) and the solution
$\tilde{\phi}'_1$, obtained using the
IMF differential operator for $(k, q, n) = (4, 4, 180)$, from
Example <ref>.

[Sigmoid up-chirp]
An aspect of interest is an IMF with rapidly increasing frequency in a short
timeframe. We want to know how well our differential operator handles such a case.
A signal whose frequency changes over time is called a chirp, and one
with increasing frequency over time an up-chirp.
Even with a rapid increase, as with all natural phenomena, we can reasonably expect
our frequency to still be smooth. This is best illustrated if we compare the
frequency with velocity. We can not have sudden changes in velocity of an object
either, as it would imply infinte acceleration in that moment.
To take the idea further, we can not have sudden changes in the acceleration either,
as it would imply infinite jerk (rate of change of acceleration) in that moment,
et cetera.
A good modelling function for this is a sigmoid
function, more precisely the logistic function, which we will make use of in this
Consider the simple IMF

\begin{equation}\label{eq:example_diffop-2-u}
u_2(t) := \cos(\phi_2(t)) := \cos\left( 40 \cdot t + \frac{100}{90} \cdot
\ln\!\left( 1 + \exp\!\left(90 \cdot (t - 0.5)\right)\right)\right)
\end{equation}

on the interval $[0,1]$ (see Figure <ref>). We calculate
the instantaneous frequency $\phi'_2$, which happens to be a transformation of the
logistic function, analytically as

\begin{equation}\label{eq:example_diffop-2-freq-analytic}
\phi'_2(t) = 40 + \frac{100}{1 + \exp\!\left( -90 \cdot (t - 0.5) \right)}.
\end{equation}

You can find it pictured in Figure <ref>. It
represents a sudden frequency increase from $40$ to $140$ in a very short
timeframe around the middle of the interval $[0,1]$.
As in Example <ref>, given the instantaneous amplitude is constantly $1$, we
can use the toolbox to fit the simple IMF differential operator to $u_2$ to calculate the
numerical instantaneous frequency $\tilde{\phi}'_2$.
The difference between $\phi'_2$ and $\tilde{\phi}'_2$
is too small to be visible in a normal plot and thus we examine the semi-log plot of the
relative error (see Figure <ref>). We can see that the
relative error is at most roughly $0.1\%$ and ranges between around
$2$ orders of magnitude below that.

[xlabel=$t$, ylabel=$u_2(t)$]
[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/regop.data/2-u.csv;

Plot of the IMF $u_2$ (see (<ref>)) from
Example <ref>.

[xlabel=$t$, ylabel=$\phi'_2(t)$]
[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/regop.data/2-freq-analytic.csv;

Plot of the analytical instantaneous frequency
$\phi'_2$ (see (<ref>))
from Example <ref>.

[xlabel=$t$, ylabel near ticks,
ylabel=$\left| \frac{\phi'_2(t) -
\tilde{\phi}'_2(t)}{\phi'_2(t)} \right|$]
[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/regop.data/2-relerr.csv;

Semi-log plot of the relative error between the analytical
instantaneous frequency $\phi'_2$
(see (<ref>)) and the solution
$\tilde{\phi}'_2$, obtained using the
IMF differential operator for $(k, q, n) = (4, 4, 180)$, from
Example <ref>.

§ DISCUSSION

In this chapter we have examined the IMF differential operator
(see Definition <ref>) as a possible means to extract
instantaneous amplitude and frequency from an IMF and to
regularize the EMD optimization problem introduced in
Definition <ref>. What we noticed in Theorem <ref> is that the IMF differential
operator is not convex in the parameters corresponding to amplitude
and frequency, which is why we modified it in
Definition <ref> to work only on IMFs with
constant amplitude $1$ and only extract the frequency, which we proved
in Corollary <ref>. Consequently,
we showed in the examples in Section <ref> that
the operator, as expected, successfully extracts the frequency from given
IMFs with constant amplitude $1$.
The limitatin to IMFs with constant amplitude $1$ appears to be a drastic
limitation, but we show in Chapter <ref> how to work around it and
extract the amplitude already during the EMD sifting process.
However, when considering the EMD optimization problem from Definition <ref> the general differential operator from Definition <ref> is unsuitable as a general regularization
term $R(a,\phi)$ from a theoretical perspective.
Given we have shown that it is not convex, one can even
consider it to be more of a heuristic tool.

CHAPTER: HYBRID OPERATOR-BASED METHODS

This chapter is the culmination of the results obtained in the previous
chapters. In Chapter <ref> we formulated and analyzed
the optimization problem (see Equation (<ref>))

\begin{equation*}
\begin{aligned}
\min_{u} \quad &
{\| s - \I[a,\phi] \|}_2^2 \\
\text{s.t.} \quad & (a,\phi) \in {\mathcal{S}}_{\mu_0,\mu_1,\mu_2}
\end{aligned}
\end{equation*}

In particular we have shown that this optimization problem can also be
expressed as a regularized optimization problem
(see Section <ref>)

\begin{align*}
\min_{(a,\phi) \in {\left( \mathcal{C}^0(\R,\R) \right)}^2} \quad &
c[s](a,\phi) + R(a,\phi).
\end{align*}

$R(a,\phi)$ is a regularization term that punishes solution candidates
of the optimization problem that are not an IMF soul
(see Definition <ref>).
In Chapter <ref> we introduced one possible way to define
this regularization term. We made use of a differential operator
$\mathcal{D}_{(\tilde{a},\tilde{\phi})}$ with parameters
$\tilde{a}$ and $\tilde{\phi}$ that annihilates IMF functions
$u(t) := a(t) \cdot \cos(\phi(t))$ when $\tilde{a} = a$ and
$\tilde{\phi} = \phi$. Unfortunately, this operator does not
yield unique results for a given input IMF
(see Theorem <ref>) and thus is
highly reliant on heuristics to work.
However, when reduced to input IMFs with constant amplitude $1$,
the differential operator is convex and the resulting frequency for
an IMF is unique (see Corollary <ref>). This result does not look very useful,
but can be leveraged when combined with the classic EMD method
proposed by <cit.>.
The result of this combination is a hybrid of classic and
modern methods and will be introduced later in this chapter. Before
considering this approach, we first introduce the classic EMD method.

§ CLASSIC EMD METHOD

The classic EMD method was first proposed in <cit.>
and will be described as follows. The EMD is a multistep
method, but we will without loss of generality
only consider a single extraction step.
In this step we separate a given multicomponent signal $s(t)$
into an IMF $u(t)$ and a residual $r(t)$.
This is without loss of generality, as subsequent extraction
steps are realized by considering the
residual of the previous step as the input signal for the current
step. Continuing this process, we sooner or later obtain a residual
that does not contain any more IMFs. The stopping criterion
might for instance be when the residual has no or at
most one local extremum, but this is not within the scope of this thesis.
After the separation of $s(t)$ into $u(t)$ and $r(t)$,
one can determine the instantaneous amplitude $a(t)$
and phase $\phi(t)$ of $u(t) := a(t) \cdot \cos(\phi(t))$ by
complexification of $u(t)$ using the Hilbert
transform, which will not be further elaborated here.
We note here though that this Hilbert transform
provides some numerical challenges. In particular,
it requires heuristics to work properly in the numerical
context, which is why alternatives to this approach are
desired and presented in this thesis.
The process of
separation is called sifting in the original
paper <cit.> and commonly referred to as the
empirical mode decomposition (EMD). The extraction of $a(t)$
and $\phi(t)$ from $u(t)$ is called the
Hilbert spectral-analysis (HSA).
The complete process of EMD and HSA is referred
to as the Hilbert-Huang-transform
(HHT). Of note here is though that because the
HHT describes a very specific approach using the
Hilbert transform, one finds that the term
EMD is often used to also include the
spectral analysis part that makes use of some
other method.
The sifting method of separating the signal $s(t)$ into
an IMF $u(t)$ and residual $r(t)$ can be separated into
three steps, illustrated in Figure <ref>
and given in Algorithm <ref>.

[domain=0:1, samples=200]
2 + 3 * (1 + x)^3 +
0.8 * (2 + cos(deg(4 * pi * x))) * sin(deg(30 * pi * x));
[domain=0:1, samples=200, densely dashed]
2 + 3 * (1 + x)^3 +
0.8 * (2 + cos(deg(4 * pi * x))) * 1;
[domain=0:1, samples=200, densely dashed]
2 + 3 * (1 + x)^3 +
0.8 * (2 + cos(deg(4 * pi * x))) * (-1);

envelope estimation

[domain=0:1, samples=200, densely dashed]
2 + 3 * (1 + x)^3;
[domain=0:1, samples=200]
2 + 3 * (1 + x)^3 +
0.8 * (2 + cos(deg(4 * pi * x))) * sin(deg(30 * pi * x));

envelope averaging

[domain=0:1, samples=200]
2 + 3 * (1 + x)^3;
[domain=0:1, samples=200]
0.8 * (2 + cos(deg(4 * pi * x))) * sin(deg(30 * pi * x));
[domain=0:1, samples=200, densely dashed]
0.8 * (2 + cos(deg(4 * pi * x))) * 1;

Visualization of the EMD sifting process of a $1$-component

multicomponent signal $s \in \mathcal{C}^2(\R,\R)$

intrinsic mode function $u \in \mathcal{C}^2(\R,\R)$

instantaneous amplitude $a \in \mathcal{C}^2(\R,\R)$

residual $r \in \mathcal{C}^2(\R,\R)$

$\ul{a} \longleftarrow \text{LowerEnvelope(s)}$
$\ol{a} \longleftarrow \text{UpperEnvelope}(s)$
$r \longleftarrow \frac{1}{2} \cdot (\ul{a} + \ol{a})$
$u \longleftarrow s - r$
$a \longleftarrow \ol{a} - r$

EMD sifting algorithm.

The first step is to estimate the lower and upper envelopes $\ul{a}(t)$
and $\ul{a}(t)$ of the input signal. What an envelope is exactly will
be defined later. The second is to take the average of $\ul{a}(t)$
and $\ul{a}(t)$, yielding the residual $r(t)$, and the third is to
separate the signal into residual and IMF $u(t)$ by subtracting $r(t)$
from $s(t)$. The instantaneous amplitude $a(t)$ of $u(t) :=
a(t) \cdot \cos(\phi(t))$ follows naturally by subtracting $r$ from
the upper envelope $\ol{a}$.
From this observation we can conclude two things: The first is that
the envelope estimation is central to the EMD method. The second
is that given we obtain the IMF $u(t)$ and its instantaneous amplitude
$a(t)$ naturally from the sifting process, we can make use of our
differential operator to extract the instantaneous phase $\phi(t)$.
This is because the IMF $\tilde{u}(t) := u(t) / a(t)$ has amplitude
$1$ and makes it possible to use the differential operator introduced
in Chapter <ref> in a theoretically meaningful way.
What is left to do is to analyze the envelope estimation method itself,
which we will do as follows.

§ ENVELOPE ESTIMATION

An envelope is not uniquely classified, but defined as a
function that encloses a function either from above (upper
envelope) or below (lower envelope).

Let $v,f \in \mathcal{C}^0(\R,\R)$. $v$ is a lower envelope
of $f$ if and only if

\begin{equation*}
v \preceq f.
\end{equation*}

$v$ is an upper envelope of $f$ if and only if

\begin{equation*}
v \succeq f.
\end{equation*}

As an example, an upper envelope for $\cos(t)$ is the constant function
$1$ and a lower envelope is the constant function $-1$, but we can also
choose $2$ and $-2$ or $\cos(t)$ for both as lower and upper amplitudes
(see Figure <ref>).

[domain=0:21, samples=200]cos(deg(x));
[densely dashed, domain=0:21, samples=200]0.2 + 0.8 * cos(deg(x)));
[densely dashed, domain=0:21, samples=200]0.5 + 0.5 * cos(deg(x)));
[densely dashed, domain=0:21, samples=20]1;
[densely dashed, domain=0:21, samples=20]2;

Examples (dashed) for upper envelopes of $\cos(t)$ (solid).

We can thus note that by far there is no unique choice for a lower
and upper envelope of a function and we will have to specify more
requirements the envelopes have to fulfill.
In the context of the empirical mode decomposition, determining the
lower and upper envelopes of an input signal is the central step to
obtain the residual and IMF, as explained in Section <ref>.
As follows, we will, without loss of generality, only consider the
upper envelope estimation. The procedure for the lower envelope follows
respectively, given the following

Let $f \in \mathcal{C}^0(\R,\R)$ and $v \in \mathcal{C}^0(\R,\R)$
be a lower envelope of $s$.
It holds that $-v$ is an upper envelope of $-f$.

It holds by Definition <ref> that $v \preceq f$ and

\begin{equation*}
v \preceq f
\Leftrightarrow -v \succeq -f. \qedhere
\end{equation*}

Thus, to determine the lower envelope we simply determine the negated
upper envelope of the negated input function.

§.§ Classic Envelope Estimation

Knowing the requirements for an upper envelope in the context of the
empirical mode decomposition listed previously, we now take a look at the
classic envelope estimation proposed in <cit.>.
When we reconsider the previous example $\cos(t)$ (which has upper envelope $1$)
we see that the function assumes the value
$1$ in its local maxima. Consequently, we can propose that an IMF assumes
the value of its upper envelope in its local maxima and we obtain the upper
envelope by interpolating them. The corresponding algorithm in pseudocode
can be found in Algorithm <ref>.
multicomponent signal $s \in \mathcal{C}^2(\R,\R)$
upper envelope $m \in \mathcal{C}^2(\R,\R)$
$P \longleftarrow \{ (t, s(t)) \in \R\times\R \mid s'(t) = 0 \land s''(t) < 0 \}$
$m \longleftarrow \text{Interpolate}(P)$

Classic upper envelope estimation algorithm.

The Interpolate-method in the Algorithm is left out by choice and
means the fitting of a B-spline-curve to each point in the set $p$.
The problem is that with varying amplitude the estimated envelope tends to
dip below the signal, thus violating the definition of an envelope not to
cut the signal at any moment. This is illustrated
in Figure <ref>.

[spy using outlines=rectangle,lens=scale=5,
size=2.5cm, connect spies]
[mark=none, smooth, very thin] table [x=x, y=y, col sep=comma] examples/envelope.data/x-s.csv;
[mark=none, smooth, ultra thin] table [x=x, y=y, col sep=comma] examples/envelope.data/x-a-classic.csv;
[mark=none, smooth, ultra thin, densely dashed] table [x=x, y=y, col sep=comma] examples/envelope.data/x-a-analytic.csv;
(spypoint_0) at (axis cs:0.321,7.5);
(glass_0) at (axis cs:0.16,18.3);
(spypoint_1) at (axis cs:0.927,21);
(glass_1) at (axis cs:0.84,4.2);

on (spypoint_0) in node[fill=white] at (glass_0);
on (spypoint_1) in node[fill=white] at (glass_1);

Multicomponent signal $s_0$ (thick) from Example <ref>
with analytical upper envelope (dashed) and estimated upper envelope using
the classic sifting algorithm (thin). Two sections where the latter cuts
the signal are enlarged.

As we can see, the classic method of interpolating the local maxima
reaches its limits very quickly and is in general not a very good
envelope estimation method, given it violates the definition.

§.§ Iterative Slope Envelope Estimation

There have been multiple approaches to the problem
with the intersection of envelope and signal that we described earlier.
<cit.> introduced an optimization scheme to obtain
the envelope, strictly enforcing the nature of the envelope definition,
but at the cost of the smoothness of the resulting amplitude estimation.
<cit.> approached the problem by analytically moving the interpolation
points from the local maxima to more fitting spots, with the disadvantage
that these approaches only work where it is at least possible to estimate
the current frequency. Additionally, it only allows to work with IMFs and
not a multicomponent signal, which we are relying on in the sifting process,
because it is not possible to easily find analytical results taking the
entire multicomponent signal into account.
If we take a step back and think how a human would draw an uppper envelope of
a signal $s(t)$ by hand, we see that the result $m_0(t)$ of the classic
sifting can be considered as a first step toward a better envelope
estimation which just needs some refining.
We do that by taking $m_0(t)$ and finding every point on the signal $s(t)$
where $m_0'(t) = s'(t)$ (matching slope) and $s''(t) < 0$ (negative curvature)
hold. We obtain the upper envelope $m_1(t)$ by interpolating these points.
Repeating this process yields a curve with a better fitting, as it
becomes by definition a tangential curve.
The algorithm describing this process can be found as pseudocode in
Algorithm <ref>. We begin
with a multicomponent signal $s(t)$ and a tolerance.
Our estimated upper envelope $m$ is first initialized to
the zero-function before entering the main loop, in which $m$
is copied to $\tilde{m}$ and $m$ set to the
next envelope estimate iterate. If the difference between
the previous and current envelope estimate iterate
is strictly smaller than our tolerance $\varepsilon$
in the supremum norm, we are done.
We make use of the supremum norm given it is easy
to calculate a close upper bound of it within the
well-conditioned B-spline basis (see
<ref>), which amounts
to just the supremum norm of the respective vector
of  B-spline basis coefficients.
multicomponent signal $s \in \mathcal{C}^2(\R,\R)$

tolerance $\varepsilon > 0$

upper envelope $m \in \mathcal{C}^2(\R,\R)$
$m \longleftarrow 0 \in \mathcal{C}^0(\R,\R)$
${\| m - \tilde{m} \|}_\infty < \varepsilon$
$\tilde{m} \longleftarrow m$
$p \longleftarrow \{ t \in \R \mid s'(t) = \tilde{m}'(t) \land s''(t) < 0 \}$
$m \longleftarrow \text{Interpolate}(p)$

Iterative slope upper envelope estimation algorithm.

Let us compare Algorithms <ref>
and <ref>.
We remind ourselves that to determine the upper envelope,
the classic method interpolates the local maxima.
The slope in the maxima is $0$ and the curvature is
negative. It is easy to see that the first iteration of the
iterative slope envelope algorithm is simply the classic
envelope estimation, because the slope of the $0$-function
is also zero. Thus, all slope matches in the signal are those
where the slope is zero.
Forcing the algorithm to finish after the first iteration by
setting $\varepsilon = \infty$ we obtain the classic
method. We can thus say that the proposed upper envelope
estimation algorithm is a generalization of the classic

Obtaining the lower envelope of a given input signal $s$ is analogous
to Proposition <ref> by
determining the negative upper envelope of the negated input
signal $-s$. Given these negations are linear time operations
there is no effect on the run-time of the algorithm regardless
of whether we estimate the upper or lower envelope.
An advantage of this algorithm over the method presented in <cit.>
is that we do not need to estimate the instantaneous frequency and do
not require the input signal to have any special form. Given our new
method is a generalization of the classic envelope estimation,
it fits more naturally into the existing methods. Moreover, we
solve the intersection problem as described in
Figure <ref> and obtain meaningful envelopes
that satisfy the definition.

§.§ Examples

The following examples were implemented using the ETHOS-toolbox developed
in the course of this thesis and can be found in Listing <ref>.
The parameters $(k,q,n,\varepsilon)$ given in the figure captions refer to the
spline order $k$, in-fill-count $q$ (see Section <ref>),
number of B-spline basis functions $n$ (see Definition <ref>) and envelope extraction
tolerance $\varepsilon$ (see Algorithm <ref>).
See Subsection <ref> for a discussion on the boundary
effects of these examples in the context of information theory and other literature.
Consider the composite signal

\begin{equation}\label{eq:example_envelope-0-s}
s_0(t) := 40 \cdot t + (20 + 10 \cdot \cos(5 \cdot \pi \cdot t))
\cdot \cos(25 \cdot \pi \cdot t)
\end{equation}

on the interval [0,1]. Beginning with the highest frequency
component, the first analytical envelope to be extracted by the sifting process

\begin{equation}\label{eq:example_envelope-0-m}
m_0(t) := 40 \cdot t + (20 + 10 \cdot \cos(5 \cdot \pi \cdot t)).
\end{equation}

In Figure <ref> you can see the result of
the proposed iterative slope sifting process compared with the
analytical envelope $m_0$.
Due to the little differences in most
parts we examine the semi-log plot of the relative error
(see Figure <ref>) for both the classic
and iterative slope sifting processes. We can see that the
relative error for the proposed iterative slope method is up to
an order of magnitude less in
some parts while staying equally good e.g. in the boundary regions,
which is more due to an information theoretical reason and not a
quality criterion of the sifting algorithm.
What is more important is that the iterative slope envelope is a true
envelope in that it does not cut the signal in any location like
the envelope obtained with the classic sifting algorithm. This is
due to the fact that the iterative slope envelope is a tangent
by construction.

[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/envelope.data/0-s.csv;
[mark=none, smooth, thin] table [x=x, y=y, col sep=comma] examples/envelope.data/0-a.csv;
[mark=none, smooth, dotted] table [x=x, y=y, col sep=comma] examples/envelope.data/0-a-analytic.csv;

Plot of the multicomponent signal $s_0$ (see (<ref>)), its analytical highest frequency
component envelope $m_0$ (dotted) (see (<ref>))
and calculated iterative slope envelope (thin) for
$(k, q, n, \varepsilon) = (4, 4, 180, 0.01)$ from
Example <ref>.

[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/envelope.data/0-relerr.csv;
[mark=none, smooth, dotted] table [x=x, y=y, col sep=comma] examples/envelope.data/0-relerr-classic.csv;

Semi-log plot of the relative errors between the analytical
highest frequency component envelope $m_0$
(see (<ref>)) and both the envelopes
obtained using the classic (dotted) and iterative slope
sifting algorithms for $(k,q,n,\varepsilon) = (4,4,180,0.01)$
from Example <ref>.

Consider the IMF

\begin{equation}\label{eq:example_envelope-1-s}
s_1(t) := \frac{1}{16} \cdot (t^2 + 2) \cdot
\cos(\pi \cdot \sin(8 \cdot t) + \pi)
\end{equation}

on the interval [-4,4]. The analytical envelope to be
extracted by the sifting process is

\begin{equation}\label{eq:example_envelope-1-m}
m_1(t) := t^2 + 2.
\end{equation}

In Figure <ref> you can see the result of
the proposed iterative slope sifting process compared with the
analytical envelope $m_1$.
Due to the little differences in most
parts we examine the semi-log plot of the relative error
(see Figure <ref>) for both the classic
and iterative slope sifting processes. We can see that the
relative error for the proposed iterative slope method is equal
to that of the classic sifting method and even up to an order of magnitude
lower in the increasing branch of $s_1$.
The reason the error is
not symmetric like that of the classic sifting method is because
even though the maxima are symmetrically distributed, $s_1$ itself
is not symmetric. The classic sifting method only considers the maxima
though and thus is oblivious to the shape of $s_1$ itself, unlike
the iterative slope method.

[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/envelope.data/1-s.csv;
[mark=none, smooth, thin] table [x=x, y=y, col sep=comma] examples/envelope.data/1-a.csv;
[mark=none, smooth, dotted] table [x=x, y=y, col sep=comma] examples/envelope.data/1-a-analytic.csv;

Plot of the multicomponent signal $s_1$ (see (<ref>)), its analytical highest frequency
component envelope $m_1$ (dotted) (see (<ref>))
and calculated iterative slope envelope (thin) for
$(k, q, n, \varepsilon) = (4, 4, 180, 0.1)$ from
Example <ref>.

[mark=none, smooth] table [x=x, y=y, col sep=comma] examples/envelope.data/1-relerr.csv;
[mark=none, smooth, dotted] table [x=x, y=y, col sep=comma] examples/envelope.data/1-relerr-classic.csv;

Semi-log plot of the relative errors between the analytical
highest frequency component envelope $m_1$
(see (<ref>)) and both the envelopes
obtained using the classic (dotted) and iterative slope
sifting algorithms for $(k,q,n,\varepsilon) = (4,4,180,0.1)$
from Example <ref>.

What we can clearly see is that the newly presented iterative slope sifting
algorithm provides a better envelope estimation than the classic sifting algorithm. Of note is especially the intuition behind it and the fact that
it is a generalization of the classic method. For this reason, we will
make use of it in our hybrid EMD algorithm presented in the next section.

§ HYBRID EMD ALGORITHM

This section presents a new EMD algorithm making use of
the new iterative slope sifting algorithm and the differential-
operator-based method presented earlier. Given the former is
considered a more classic approach compared to the operator-based
signal-separation and the latter operator-based method is a modern
concept it is fitting to call this algorithm a hybrid
algorithm. The entire procedure for a single step of the decomposition
is given in Algorithm <ref>, but we will construct the method
step by step in the following section for a given multicomponent signal $s$.
The complete decomposition is obtained by successive runs of the
algorithm with the residual subtracted from the signal as the input
for the next step.
multicomponent signal $s \in \mathcal{C}^2(\R,\R)$

tolerance $\varepsilon > 0$

intrinsic mode function $u \in \mathcal{C}^2(\R,\R)$

instantaneous amplitude $a \in \mathcal{C}^2(\R,\R)$

instantaneous frequency $\phi' \in \mathcal{C}^2(\R,\R)$

residual $r \in \mathcal{C}^2(\R,\R)$

$\ol{a} \leftarrow \text{UpperEnvelope}(s,\varepsilon)$
$\ul{a} \leftarrow -\text{UpperEnvelope}(-s,\varepsilon)$
$r \leftarrow \frac{1}{2} \cdot (\ol{a} + \ul{a})$
$u \leftarrow s - r$
$a \leftarrow \ol{a} - r$
$\Omega^\star \leftarrow \arg\min_{\Omega \in \mathcal{C}^2(\R,\R)}
\left(
{\left\|\tilde{\mathcal{D}}_{(0,\Omega)}(u / a)\right\|}_2^2
\right)$
$\phi' \leftarrow \frac{1}{\sqrt{\Omega^\star}}$

Hybrid operator-based empirical mode decomposition algorithm.

The first step is to determine the upper and lower envelopes $\ol{a}$
and $\ul{a}$ of $s$.
We make use of the fact that the lower envelope is just the negation
of the upper envelope of the negated signal $s$. This is why we
previously only considered the upper envelope estimation, as the lower
envelope estimation follows as a corollary.
The idea behind the following steps to obtain $r$, $u$ and $a$ were first
introduced in <cit.>. Once we've determined $\ol{a}$ and $\ul{a}$
we can calculate the residual $r$ as their mean.
The intrinsic mode function $u$ is obtained by subtracting $r$ from $s$
and the instantaneous amplitude $a$ is calculated by subtracting $r$ from
Next we use the simple case (see Subsection <ref>)
of the modified IMF differential operator
(see Definition <ref>)
to obtain our instantaneous frequency. We first solve the NSP optimization
problem for the inverse square continuous frequency operator $\Omega[\phi]$
(see Definition <ref>) and then calculate
the instantaneous frequency directly by applying the inverse square

§ ETHOS TOOLBOX

The central numerical piece of this thesis is the ETHOS toolbox.
It stands
for EMD Toolbox using Hybrid
Operator-Based Methods and B-splines and has been
developed in the course of this thesis to provide an implementation
for the new concepts presented in this work, making it possible
to do an empirical mode decomposition on a discrete input signal.
The implementation language is C99 (see <cit.>),
making use of the GNU Scientific Library (see <cit.>) for
the numerical backend (including B-splines).
All examples in this and previous chapters have been realized in
this toolbox (see Section <ref>) and the
entire source code is listed under Section <ref>.
The main header exposing the toolbox function is
(see Listing <ref>) with
the main datatype . Nearly all
functions take a  as input and it is
the main storage for system parameters and precomputed
As follows, we will take a look at the most important
functions with regard to the decomposition process. For
all exposed functions refer to
(see Listing <ref>).

§.§ Initialization and Precomputation

The main initialization and first step of any program
using the ETHOS toolbox is done by
(see Listing <ref>). It takes a
pointer to an ETHOS-struct , vector
of length  with spline
order , in-fill-count ,
density  and grid-type  as input and
fills the given ETHOS-struct with the necessary parameters
and precomputes data for later use. The parameters are
explained as follows.
The vector  contains all time-steps of the
discrete input signal, or comparable like a superset of
multiple possible interpolation areas, and it is our
interest to only take a certain subset of these steps for our spline
knot-vector. This is controlled by the parameters
, controlling the density
of the spline knot vector relative to the input vector
and residing in the interval $(0,1)$, and ,
controlling the way the selection is made
(uniformly or adaptively).
When the spline knot-vector is obtained, it is uniformly
in-filled with  points between each spline knot.
The motivation for this process is to be able to pre-evaluate
the splines and their derivatives during the initialization step
on this extended grid (i.e. the in-filled knot-vector).
The reason why the knot-vector is not just made denser is
because we want to limit the number of basis functions and
for the sake of plotting or general evaluation do not need so
many basis functions, because B-splines as is provide a great
amount of smoothness.
The evaluation happens only for the non-zero parts of each
basis function, resulting in linear memory complexity for
this precomputation step.
On this extended grid the B-splines and their first and
second derivatives are evaluated and the results stored
in  within the ETHOS struct, just like the
spline knot-vector in  and in-filled form
in .
The size of the B-spline basis  is of great
importance, as each function is stored internally as a
vector of this length, corresponding as coefficients
of the B-spline basis.

§.§ Data Filtering

After initialization, the next step is to filter the
discrete input data, with the goal of obtaining the
B-spline coefficients for this given function on the
initialized grid. The function for this purpose is
(see Listing <ref>)
and it takes an arbitrary discrete
input signal  with the time-steps ,
which do not have to agree completely with the
used in the initialization, but should agree on the
start- and endpoints.
During the fitting process, the B-splines of the
initialized basis are evaluated on all points in
and a weighted least-squares system solved with
a set of low-weight smoothness-terms of second order
besides the interpolation terms for each given datapoint.
It returns the B-spline coefficients in  best fitting the
given discrete input data as a vector, which is the
standard way of handling functions within the toolbox.
Even though not explicitly expressed, all of the
coefficient vectors have length  found
in the ETHOS struct.

§.§ Boundary Effects and Extension

The boundary effects we have seen in the previous
examples opens up some questions that will be
addressed here. As you can, for example, see
in Figure <ref>
the error goes up as it reaches the boundary.
A common countermeasure often (silently) employed
in the literature is to extend the signal beyond
the boundary, either by mirroring or other methods
(see <cit.> for further reading).
This way, the shock the algorithm is exposed to
is moved into the mirrored section or dampened, not
as heavily affecting the interior part one actually cares about.
When looking at this matter in an information
theoretical way this technique of extension is
rather dishonest about the performance of such
an algorithm and generates information where there
is none. It might be forgivable for applications
that care about a good represenation, but the real
challenge is to design robust algorithms and make them
comparable among each other without silent tricks like
this one. Moreover, there is not a canonical way to
extend beyond the boundary and it presents itself
more as its own field of research. This is the reason
why the author chose not to use boundary extension
methods for his examples and keep them
honest with regard to the boundary effects.
Despite these ethical concerns, the ETHOS toolbox includes
(see Listing <ref>)
which takes an arbitrary discrete input signal  with the
time-steps  and length  and calculates an extended signal
with the time-steps  and length
by mirroring the signal into the extended area.
This extension is parametrized by  between $0$ and $1$,
extending the signal by this fraction both on the left and right side.

§.§ Decomposition

The function to do the signal decomposition itself is
(see Listing <ref>).
It takes a pointer to an ETHOS-struct ,
B-spline coefficient vectors  of the output IMF,
of the amplitude of the output IMF,
of the frequency of the output IMF
and  of the input signal
of length  and tolerance .
The procedure aligns with Algorithm <ref>,
filling the input signal vector  with
the residual and ,
and  with the extracted IMF and its
amplitude and frequency respectively. The tolerance
is the tolerance for the iterative
slope algorithm (see Algorithm <ref>).
Subsequent invocations of
yield the complete decomposition.

§.§ Plotting

The ETHOS toolbox provides two ways of plotting data, either
as a CSV-output or output meant as input for the
command of the GNU plotting utilities. The input to those
plotting functions can either be a set of discrete points or
a spline function, represented with a coefficient vector.
The former is implemented as ,
the latter as  (see
Listing <ref>).

§.§ Envelope Estimation

The procedure to estimate the envelope is already used
in  described in Subsection <ref>
and not directly part of the decomposition path, however, it might
be of interest to test the envelope estimation itself separately.
This estimation is achieved with
(see Listing <ref>). It takes a pointer to an
ETHOS-struct , B-spline coefficient vectors
of the upper envelope and  of the input signal and
tolerance .
The procedure aligns with Algorithm <ref>
and stores the upper envelope estimation of  in
with the tolerance .
It can also be used to estimate the envelope using the classic
method by setting  to
(see Remark <ref>), defined
in <cit.> to represent infinity. This allows an easy comparison
of both methods, as done in Subsection <ref>.

§.§ IMF Characteristic

When we defined intrinsic mode functions in Definition <ref>
we parametrized the model with three parameters $(\mu_0,\mu_1,\mu_2)$ with
$\mu_0,\mu_1,\mu_2 > 0$, called the characteristic. We also established
the connection with the IMF accuracy presented in <cit.>
and the role and calculation of each parameter in Remark <ref>. The former especially underlines the relevance of these parameters and makes it interesting to further explore them instead of just
treating them as a theoretical tool.
For the purpose of determining the characteristic for a given IMF, the ETHOS toolbox offers the
(see Listing <ref>) function. It takes a pointer to an
ETHOS-struct , vectors  and
of length  and a vector  of length $3$,
filling  with the three characteristic values.
The vectors  and  are the B-spline
coefficients for the amplitude and frequency functions relative
to the current spline environment given with .
Even though it is not possible to control the characteristic of each
extracted IMF during the decomposition it is nevertheless possible
to ascertain the quality of the extraction afterwards using this tool.

§ EXAMPLES

These examples were implemented using the ETHOS-toolbox and can be found
in Listing <ref>.
The parameters $(k,q,n,\varepsilon)$ given in the figure captions refer to the
spline order $k$, in-fill-count $q$ (see Section <ref>),
number of B-spline basis functions $n$ and envelope extraction
tolerance $\varepsilon$ (see Algorithm <ref>).
See Subsection <ref> for a discussion on the boundary
effects of these examples in the context of information theory and other literature.

This example was inspired by <cit.>.
Consider the multicomponent signal

\begin{equation}\label{eq:example_emd-0-s}
s_0(t) := u_{0,0}(t) + u_{0,1}(t) + 20 \cdot (t + 1)
\end{equation}

with the first IMF component (characteristic
$\sim(6.60\cdot 10^1, 1.52\cdot 10^{-2}, 1.43 \cdot 10^0)$)

\begin{equation*}
u_{0,0}(t) := a_{0,0} \cdot \cos(\phi_{0,0}(t))
:= (t + 1) \cdot \cos((15 \cdot t + 21) \cdot \pi \cdot t)
\end{equation*}

and the second IMF component (characteristic
$\sim(1.57\cdot 10^1,1.91 \cdot 10^{-1}, 1.88 \cdot 10^{-12})$)

\begin{equation*}
u_{0,1}(t) := a_{0,1} \cdot \cos(\phi_{0,1}(t))
:= (3 \cdot t + 1) \cdot \cos(5 \cdot \pi \cdot t).
\end{equation*}

We calculate the instantaneous frequencies $\phi'_{0,0}$ and
$\phi'_{0,1}$ of both IMF components $u_{0,0}$ and $u_{0,1}$ as

\begin{equation*}
\phi'_{0,0} = 15 \cdot \pi \cdot t + (15 \cdot t + 21) \cdot \pi
\end{equation*}

\begin{equation*}
\phi'_{0,1} = 5 \cdot \pi.
\end{equation*}

Our objective is to run a full EMD on this input signal $s_0$, which
means that we, in each step, identify an IMF that we will further
analyze to obtain its instantaneous amplitude and frequency.
Due to the nature of our HOST-EMD
algorithm we first extract high-frequency components only to continue
to extract successively lower frequency components in subsequent steps,
corresponding to the target in this example to extract $u_{0,0}$ first
and then $u_{0,1}$.
In the first step we find the IMF $\tilde{u}_{0,0}$ with instantaneous
amplitude $\tilde{a}_{0,0}$ and frequency $\tilde{\phi}'_{0,0}$
(see Figure <ref>). Splitting $u_{0,0}$
analytically from the input signal $s_0$ we obtain our first residual

\begin{equation}\label{eq:example_emd-0-r-1}
r_{0,1}(t) := s_0 - u_{0,0}
\end{equation}

(see Figure <ref>), whose calculated form
$\tilde{r}_{0,1}$ we will use as the input signal for our second step.
Analogously, we find the IMF $\tilde{u}_{0,1}$ with instantaneous
amplitude $\tilde{a}_{0,1}$ and frequency $\tilde{\phi}'_{0,1}$
in the second
step (see Figure <ref>) and further
splitting $u_{0,1}$ from $r_{0,1}$ analytically yields the second

\begin{equation}\label{eq:example_emd-0-r-2}
r_{0,2}(t) := r_{0,1} - u_{0,1} = s_0 - u_{0,0} - u_{0,1}
\end{equation}

(see Figure <ref>) with its calculated form
$\tilde{r}_{0,2}$, which we identify
as the last residual given it obviously contains no further
IMF components.
Discussing errors is more difficult than in the other examples
presented in this thesis as we have two algorithms working in concert,
namely the proposed iterative slope sifting and differential operator
extraction algorithms. Nevertheless, they turn out to be working
independently and what we can note is that the only significant errors
are visible at the boundaries of the signal, which is to be expected
as it is blind for the analytical nature of the input signal.
In the interior of the signal, the relative errors range
between $10^{-2}$ and $10^{-6}$.

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

Plot of the multicomponent signal
$s_{0}$ (see (<ref>)) from
Example <ref>.

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

Plots of the first extracted
IMF-component $\tilde{u}_{0,0}$ of $s_0$
(see (<ref>))
and its respective instantaneous amplitude $\tilde{a}_{0,0}$ and
frequency $\tilde{\phi}'_{0,0}$ in the first row (with
analytical solutions $u_{0,0}, a_{0,0}$ and
$\phi'_{0,0}$ (dotted)) with
semi-log plots of the absolute (for $u_{0,0}(t)$) and
relative errors (for $a_{0,0}(t)$ and $\phi'_{0,0}(t)$) compared to the
analytical solutions in the second row for
$(k, q, n, \varepsilon) = (4, 4, 180, 0.01)$
from Example <ref>.

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

Plots of the calculated first residual $\tilde{r}_{0,1}$
(see (<ref>)) and its analytical solution
$r_{0,1}$ (dotted) on
the left and a semi-log plot of the relative error between both
on the right for $(k, q, n, \varepsilon) = (4, 4, 180, 0.01)$
from Example <ref>.

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

Plots of the second extracted
IMF-component $\tilde{u}_{0,1}$ of $s_0$
(see (<ref>))
and its respective instantaneous amplitude $\tilde{a}_{0,1}$ and
frequency $\tilde{\phi}'_{0,1}$ in the first row (with
analytical solutions $u_{0,1}, a_{0,1}$ and
$\phi'_{0,1}$ (dotted)) with
semi-log plots of the absolute (for $u_{0,1}(t)$) and
relative errors (for $a_{0,1}(t)$ and $\phi'_{0,1}(t)$) compared to the
analytical solutions in the second row for
$(k, q, n, \varepsilon) = (4, 4, 180, 0.01)$
from Example <ref>.

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

Plots of the calculated second residual $\tilde{r}_{0,2}$
(see (<ref>)) and its analytical solution
$r_{0,2}$ (dotted) on
the left and a semi-log plot of the relative error between both
on the right for $(k, q, n, \varepsilon) = (4, 4, 180, 0.01)$
from Example <ref>.

This example was inspired by <cit.>.
Consider the multicomponent signal

\begin{equation}\label{eq:example_emd-1-s}
s_1(t) := u_{1,0}(t) + u_{1,1}(t) + 25 \cdot t^3
\end{equation}

with the first IMF component (characteristic
$\sim(2.51\cdot 10^2, 1.19 \cdot 10^{-2}, 1.88 \cdot 10^{-12})$)

\begin{equation*}
u_{1,0}(t) := a_{1,0} \cdot \cos(\phi_{1,0}(t))
:= (t + 1) \cdot \cos((15 \cdot t + 21) \cdot \pi \cdot t)
\end{equation*}

and the second IMF component (characteristic
$\sim(4.71 \cdot 10^1, 4.24 \cdot 10^{-1}, 2.86 \cdot 10^{-12})$)

\begin{equation*}
u_{1,1}(t) := a_{1,1} \cdot \cos(\phi_{1,1}(t))
:= (3 \cdot t + 1) \cdot \cos(5 \cdot \pi \cdot t).
\end{equation*}

We calculate the instantaneous frequencies $\phi'_{1,0}$ and
$\phi'_{1,1}$ of both IMF components $u_{1,0}$ and $u_{1,1}$ as

\begin{equation*}
\phi'_{1,0} = 15 \cdot \pi \cdot t + (15 \cdot t + 21) \cdot \pi
\end{equation*}

\begin{equation*}
\phi'_{1,1} = 5 \cdot \pi.
\end{equation*}

Our objective is to run a full EMD on this input signal $s_0$, which
means that we, in each step, identify an IMF that we will further
analyze to obtain its instantaneous amplitude and frequency.
Due to the nature of our HOST-EMD
algorithm we first extract high-frequency components only to continue
to extract successively lower frequency components in subsequent steps,
corresponding to the target in this example to extract $u_{1,0}$ first
and then $u_{1,1}$.
In the first step we find the IMF $\tilde{u}_{1,0}$ with instantaneous
amplitude $\tilde{a}_{1,0}$ and frequency $\tilde{\phi}'_{1,0}$
(see Figure <ref>). Splitting $u_{1,0}$
analytically from the input signal $s_1$ we obtain our first residual

\begin{equation}\label{eq:example_emd-1-r-1}
r_{1,1}(t) := s_1 - u_{1,0}
\end{equation}

(see Figure <ref>), whose calculated form
$\tilde{r}_{1,1}$ we will use as the input signal for our second step.
Analogously, we find the IMF $\tilde{u}_{1,1}$ with instantaneous
amplitude $\tilde{a}_{1,1}$ and frequency $\tilde{\phi}'_{1,1}$
in the second
step (see Figure <ref>) and further
splitting $u_{1,1}$ from $r_{1,1}$ analytically yields the second

\begin{equation}\label{eq:example_emd-1-r-2}
r_{1,2}(t) := r_{1,1} - u_{1,1} = s_1 - u_{1,0} - u_{1,1}
\end{equation}

(see Figure <ref>) with its calculated form
$\tilde{r}_{1,2}$, which we identify
as the last residual given it obviously contains no further
IMF components.
Discussing errors is more difficult than in the other examples
presented in this thesis as we have two algorithms working in concert,
namely the proposed iterative slope sifting and differential operator
extraction algorithms. Nevertheless, they turn out to be working
independently and what we can note is that the only significant errors
are visible at the boundaries of the signal, which is to be expected
as it is blind for the analytical nature of the input signal.
In the interior of the signal, the relative errors range
between $10^{-2}$ and $10^{-6}$.

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

Plot of the multicomponent signal
$s_{1}$ (see (<ref>)) from
Example <ref>.

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

Plots of the first extracted
IMF-component $\tilde{u}_{1,0}$ of $s_1$
(see (<ref>))
and its respective instantaneous amplitude $\tilde{a}_{1,0}$ and
frequency $\tilde{\phi}'_{1,0}$ in the first row (with
analytical solutions $u_{1,0}, a_{1,0}$ and
$\phi'_{1,0}$ (dotted)) with
semi-log plots of the absolute (for $u_{1,0}(t)$) and
relative errors (for $a_{1,0}(t)$ and $\phi'_{1,0}(t)$) compared to the
analytical solutions in the second row for
$(k, q, n, \varepsilon) = (4, 4, 180, 0.01)$
from Example <ref>.

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

Plots of the calculated first residual $\tilde{r}_{1,1}$
(see (<ref>)) and its analytical solution
$r_{1,1}$ (dotted) on
the left and a semi-log plot of the relative error between both
on the right for $(k, q, n, \varepsilon) = (4, 4, 180, 0.01)$
from Example <ref>.

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

Plots of the second extracted
IMF-component $\tilde{u}_{1,1}$ of $s_1$
(see (<ref>))
and its respective instantaneous amplitude $\tilde{a}_{1,1}$ and
frequency $\tilde{\phi}'_{1,1}$ in the first row (with
analytical solutions $u_{1,1}, a_{1,1}$ and
$\phi'_{1,1}$ (dotted)) with
semi-log plots of the absolute (for $u_{1,1}(t)$) and
relative errors (for $a_{1,1}(t)$ and $\phi'_{1,1}(t)$) compared to the
analytical solutions in the second row for
$(k, q, n, \varepsilon) = (4, 4, 180, 0.01)$
from Example <ref>.

[mark=none, smooth] table
[x=x, y=y, col sep=comma]
[mark=none, smooth, dotted] table
[x=x, y=y, col sep=comma]

[mark=none, smooth] table
[x=x, y=y, col sep=comma]

Plots of the calculated second residual $\tilde{r}_{1,2}$
(see (<ref>)) and its analytical solution
$r_{1,2}$ (dotted) on
the left and a semi-log plot of the relative error between both
on the right for $(k, q, n, \varepsilon) = (4, 4, 180, 0.01)$
from Example <ref>.

CHAPTER: SUMMARY AND OUTLOOK

In the course of this thesis we started off with the construction
of an analytical model of the empirical mode decomposition and showed
that it satisfies strong duality, namely due to being Slater
regular, using B-spline properties and the theory of convex-like
optimization. This strong duality yielded a theoretical justification for
reformulating the constrained optimization problem into an unconstrained
optimization problem with a regularization term that enforces the
In the context of EMD, we examined one possible modern approach to such a
regularization term: The operator-based signal-separation (OSS)
null-space-pursuit (NSP) method that makes use of adaptive differential
operators to determine instantaneous amplitude and phase from a given IMF.
We observed that the operator can only yield unique results when the
IMF that is to be regularized has constant amplitude $1$, which is
a strong limitation.
We considered the classic EMD algorithm and noted the following:
Sifting, the process in which the signal is separated into an
IMF and a residual, is highly dependent on a good method for estimating
the upper envelope of a multicomponent signal. We identified the
weakness in the classic method that the envelope may intersect with
the signal itself, violating the definition of an envelope.
We presented a new approach called iterative slope envelope estimation
that solves this problem. Additionally, it is a generalization of the
classic envelope estimation method.
We also made the following observation: During sifting we also
obtain the amplitude of the IMF. This in turn meant that we can just
divide the IMF by its amplitude and obtain an IMF with constant
instantaneous amplitude $1$. This meant that we could use the
differential operator from the NSP method mentioned earlier and
be sure that it behaved properly.
Using both the new envelope estimation method and the differential
operator, we have obtained an approach that is a mix of classic
and modern methods. Consequently, we defined a hybrid EMD method
and examined it using multiple examples.
These examples were implemented in a toolbox called ETHOS
using the GNU Scientific Library, which was explained and documented
subsequently. One newly discovered approach in this process was to
evaluate the quality of each extracted IMF by calculating its
It is clear that only by building a strict theoretical foundation that
went further than previous works it was possible to obtain the results
laid out in this thesis. This foundation included convex constraints
instead of regularization terms, using B-splines for the theoretical and
practical modelling of functions and the theory of convex-like
The author suspects that given the underlying EMD optimization-problem is not
convex, previous attempts to show regularity were not followed through.
This is because it is widely assumed that convexity is a requirement for showing
Slater regularity. This is wrong,
as convex-likeness is sufficient, but generally not well-known.
The attractiveness of Slater regularity is due to the fact
that it applies to the entire optimization problem instead of just specific
points, and showing it for the EMD optimization problem yields the
regularity for any scenario.
A possible outlook for further works would be to expand the iterative slope
sifting algorithm to higher dimensions to explore their usefulness in
multidimensional EMD. In the context of the ETHOS toolbox, a possible
field of research would be to use more specific tools for
sparse optimization, possibly based on the modern GHOST
(General, Hybrid and Optimized Sparse Toolkit) sparse library (see <cit.>),
and refine the B-spline data fitting process, especially concerning
preventing overfitting and underfitting.
In the end, what is easy to see is that signal analysis as a whole and
EMD in particular are complex topics with many open questions, not only
in the theoretical sense but also in practical terms. Examples include
problems like mode mixing and noise distortion of an input signal.
What remains to be seen is how these problems can be solved most
effectively: In a preprocessing step before applying the
EMD method or as a part of a newly devised EMD method. Or maybe they are
simply unsolvable in terms of information theory. No matter the outcome,
given the far-reaching applications it has in many different fields,
every little problem solved in signal analysis may have far-reaching

CHAPTER: FUNCTION SPACE ORDER AND OPERATORS

This thesis makes use of operations on functions rather than scalars,
for instance as candidates in optimization problems. Because of that we
want to understand how function spaces work and how we can map them to
other vector spaces we can handle more easily. Even though the concepts
laid out as follows seem to be very intuitive and might not even need
further explanation beyond notation, the formal aspect of this topic
shall not be missed for completeness sake but also not unnecessarily
complicate the main matter, which is why this chapter is in the appendix.
The functions we are dealing with are relatively smooth functions
$\R \to \R$.
Formally speaking, we have a $p \geq 0$ such that our function can be
differentiated $p$ times and the $p$th derivative is continuous.
This set is commonly denoted as $\mathcal{C}^p(\R,\R)$, the set of
$p$-continuously differentiable functions. It is easy to see
that $\mathcal{C}^0(\R,\R) \supseteq \mathcal{C}^1(\R,\R) \supseteq
\mathcal{C}^2(\R,\R) \supseteq \dots \supseteq \mathcal{C}^\infty(\R,\R)$,
so if we define something for $\mathcal{C}^0(\R,\R)$ it automatically holds
for all functions from $\mathcal{C}^p(\R,\R)$ with $p \geq 0$. This is why,
as follows, we will only use the set of continuous functions
$\mathcal{C}^0(\R,\R)$ in our definitions to keep everything relatively
One aspect of interest is to be able to compare two functions in some
way. In other words, we want to find an order relation on
$\mathcal{C}^0(\R,\R)$. One way to do that is as follows: We say that
a function succeeds another function if and only if the former is
pointwise greater than or equal to the latter. We define
precedence respectively and note that, obviously, there are functions
which can not be compared this way. This concept
is illustrated in Figure <ref> and is the function space
order used in the course of this thesis.

[domain=0:1, samples=200]1 + 0.2 * x^3;
[dashed, domain=0:1, samples=200]
0.2 + 0.8 * cos(deg(x)));

$f$ succeeds $g$

[domain=0:1, samples=200]cos(6 * deg(x));
[dashed, domain=0:1, samples=200]
1.2 + cos(6 * deg(x));

$f$ precedes $g$

[domain=0:1, samples=200]x);
[dashed, domain=0:1, samples=200]

$f,g$ incomparable

Illustration of an intuitive concept for a
partial function space order
(see Proposition <ref>)
exemplified with two functions $f$ (solid)
and $g$ (dashed).

Before we are able to further formalize this idea, we first need
to introduce some preliminary definitions of what an order is and
how to define ordered vector spaces.

Let $S$ be a set and $\le$ a binary relation on $S$. $\le$
is a preorder on $S$ if and only if

\begin{equation*}
\forall x \in S \colon
x \le x,
\end{equation*}

\begin{equation*}
\forall x,y,z \in S \colon
a \le b \land b \le c \Rightarrow
a \le c.
\end{equation*}

A preorder is the weakest order we can find for a vector space.
If we can also show antisymmetry of a preorder, we obtain
a partial order, given in the following

Let $S$ be a set and $\le$ a preorder on $S$. $\le$
is a partial order on $S$ if and only if
$\le$ is antisymmetric, i.e.

\begin{equation*}
\forall a,b \in S \colon
a \le b \land b \le a \Rightarrow a = b.
\end{equation*}

One step higher would be a total order, which adds the
connex property meaning that all elements of the set are
comparable. As we've seen in Figure <ref> this
is not possible for vector spaces given the presence of
indeterminate relations.
If we find a binary relation on a vector space that is
either a preorder or partial order we would, as the next
step, hope that the orders are compatible with vector
space operations. In other words, this means that addition
and scalar multiplication preserve relations intuitively.
This is reflected in the following

Let $V$ be a real vector space and $\le$ a binary relation on $V$.
$V$ is a (pre)ordered vector space if and only if

$\le$ is a preorder or partial order respectively,

Addition compatibility:

\begin{equation*}
\forall x,y,z \in V \colon
x \le y \Rightarrow x + z \le y + z,
\end{equation*}

Scalar multiplication compatibility:

\begin{equation*}
\forall x,y \in V \colon
\forall \lambda \ge 0 \colon
x \le y \Rightarrow \lambda \cdot x \le
\lambda \cdot y.
\end{equation*}

With the definitions in place we can now formalize what has
been discussed at the beginning of the section and illustrated in
Figure <ref>.

Let $a,b \in \mathcal{C}^0(\R,\R)$.
$(\mathcal{C}^0(\R,\R), \preceq)$ is an ordered vector space

\begin{equation*}
a \preceq b
\forall t \in \R \colon a(t) \leq b(t).
\end{equation*}

Let $a,b,c \in \mathcal{C}^0(\R,\R)$. We first show that $\preceq$
is a partial order.

\begin{equation*}
\forall t \in \R \colon a(t) \leq a(t)
\Leftrightarrow
a \preceq a,
\end{equation*}

\begin{align*}
a \preceq b \land b \preceq c & \Leftrightarrow
\forall t \in \R \colon (a(t) \leq b(t) \land b(t) \leq c(t)) \\
& \Rightarrow
\forall t \in \R \colon a(t) \leq c(t) \\
& \Leftrightarrow
a \preceq c,
\end{align*}

\begin{align*}
a \preceq b \land b \preceq a & \Leftrightarrow
\forall t \in \R \colon (a(t) \leq b(t) \land b(t) \leq a(t)) \\
& \Leftrightarrow
\forall t \in \R \colon a(t) = b(t) \\
& \Leftrightarrow
a = b.
\end{align*}

Let $\lambda \ge 0$. We now show that $(\mathcal{C}^0(\R,\R), \preceq)$
satisfies the two axioms

of an ordered vector space such that the order is compatible
with the vector space operations.

Addition compatibility:

\begin{align*}
a \preceq b & \Leftrightarrow
\forall t \in \R \colon a(t) \leq b(t) \\
& \Leftrightarrow
\forall t \in \R \colon a(t) + c(t) \leq b(t) + c(t) \\
& \Leftrightarrow
\forall t \in \R \colon (a + c)(t) \leq (b + c)(t) \\
& \Leftrightarrow
a + c \preceq b + c,
\end{align*}

Scalar multiplication compatibility:

\begin{align*}
a \preceq b & \Leftrightarrow
\forall t \in \R \colon a(t) \leq b(t) \\
& \Leftrightarrow
\forall t \in \R \colon \lambda \cdot a(t)
\leq \lambda \cdot b(t) \\
& \Leftrightarrow
\forall t \in \R \colon (\lambda \cdot a)(t)
\leq (\lambda \cdot b)(t) \\
& \Leftrightarrow
\lambda \cdot a \preceq
\lambda \cdot b. \qedhere
\end{align*}

We have now shown that $\mathcal{C}^0(\R,\R)$ is an ordered vector
space with the order relation $\preceq$, but this is just one point of
interest. It is often useful to transform objects that are hard to work
with to a space where that is easy, perform operations on them and
then transform them back. One common example are polar coordinate transforms
which dramatically simplify many complicated integrals.
A requirement for such a transform in the general sense is that it is an
isomorphism and preserves the nature and/or relations of the objects.
In the context of optimization problems this means that a transform shall
preserve the relative order of objects. If a function precedes another its
transform shall do the same relative to the transform of the other function.
The motivation for this approach in this thesis is found within the strong
link between spline functions and their B-spline base coefficients, as introduced in
Chapter <ref>. The base coefficients of a spline function is
a vector in $\R^n$ and fully describes it. Due to the nature of the
B-splines the intuitive partial order $\leq$ on $\R^n$ with regard to the
base coefficients is equivalent to the partial order $\preceq$ on
the set of spline functions $\Sigma_k$. To put it differently, if
a spline function succeeds another, so do their base coefficients.
The advantage this brings is apparent: Vectors in $\R^n$ are much easier
to handle both theoretically and numerically, and any optimization problem
can be trivially transformed between both representations.
Now that we have understood the motivation behind this approach, we
can formalize such a transform between ordered vector spaces in the

Let $(V,\le_V), (W,\le_W)$ be ordered vector spaces and
$f\colon V \to W$. $f$ is an order isomorphism of
ordered vector-spaces if and only if

$f$ is an isomorphism,

$\forall x,y \in V \colon x \le_V y \Leftrightarrow
f(x) \le_W f(y)$.

As we can see, preserving the order structure opens up new
possibilities with regard to optimization problems. In our case,
even though we formulate it over a function space, we can
transform our candidates into $\R^n$ and examine the equivalent
optimization problem over $\R^n$, which is much more accessible
than the function space in many aspects.
As another remark we define the following operators on functions
that serve notational purposes but can be considered to be more
or less intuitive.

Let $a \in \mathcal{C}^0(\R,\R)$.
The modulus operator $|\cdot| : \mathcal{C}^0(\R,\R) \to
\mathcal{C}^0(\R,\R)$ is defined as

\begin{equation*}
|a|(t) := |a(t)|.
\end{equation*}

Let $p \geq 0$ and $a \in \mathcal{C}^0(\R,\R)$.
The power operator $\cdot^p : \mathcal{C}^0(\R,\R) \to
\mathcal{C}^0(\R,\R)$ is defined as

\begin{equation*}
a^p(t) := {\left(a(t)\right)}^p.
\end{equation*}

Let $p \in \N_0$ and $a \in \mathcal{C}^0(\R,\R)$. The differentiation operator
$D^p: \mathcal{C}^0(\R,\R) \to \mathcal{C}^0(\R,\R)$ is defined as

\begin{equation*}
\left(D^p a\right)(t) :=
\begin{cases}
a(t) & p = 0 \\
\frac{\mathrm{d}^p a}{\mathrm{d} t^p}(t) & p > 0
\end{cases}
\end{equation*}

What shall be apparent from the results of this chapter is that
this approach saves us from a lot of cumbersome notation and complexity
down the road. The alternative of always thinking of a function as a
set of samples may also solve the discretization problem, but requires
much more care in terms of parametrization and is much harder to access
theoretically. The latter is especially apparent when it comes to derivatives,
which are a crucial part of the theory of this thesis.

CHAPTER: CONVEXITY THEORY

The theory of convex sets and functions plays a central role in this
thesis and we will introduce it in this chapter. Even though most
of the given definitions may be known to the reader, the resulting
Theorem <ref> is used in all proofs
in the thesis for showing function convexity. Before we dive
in further, we introduce both concepts of convex sets and convex functions.

Let $V$ be a real or complex vector space and $C \subseteq V$.
$C$ is convex if and only if

\begin{equation*}
\forall a,b \in C \colon \forall \lambda \in [0,1] \colon
\lambda \cdot a + (1 - \lambda) \cdot b \in C.
\end{equation*}

Let $n \in \N$, $C \subseteq \R^n$ convex and $f \colon C \to \R$.
$f$ is a strictly convex function or convex function
if and only if

\begin{equation*}
\forall a,b \in C \colon \forall \lambda \in [0,1] \colon
f(\lambda \cdot a + (1 - \lambda) \cdot b) <
\lambda \cdot f(a) + (1 - \lambda) \cdot f(b)
\end{equation*}

\begin{equation*}
\forall a,b \in C \colon \forall \lambda \in [0,1] \colon
f(\lambda \cdot a + (1 - \lambda) \cdot b) \le
\lambda \cdot f(a) + (1 - \lambda) \cdot f(b)
\end{equation*}

For the purpose of determining convexity later the multidimensional
pendant of a second derivative is introduced as the Hessian matrix,
the matrix that contains all possible combinations of
mixed second partial derivatives.

Let $n \in \N$ and $f \in \mathcal{C}^2(\R^n,\R)$. The Hessian matrix
$H_f(x) \in \R^{n \times n}$ of $f$ in $x \in \R^n$ is defined as

\begin{equation*}
H_f(x) := {\left( \frac{\partial^2 f(x)}{\partial x_i \partial x_j} \right)}_
{i,j \in \{1,\dots,n\}}.
\end{equation*}

Another important concept is the positive definiteness that
can more or less be imagined to be a generalization of
positivity of scalars into the realm of matrices.

Let $n \in \N$ and $A \in \R^{n \times n}$. $A$ is positive definite or
positive semidefinite if and only if

\begin{equation*}
\forall v \in \R^n \colon v^T \cdot A \cdot v > 0
\end{equation*}

\begin{equation*}
\forall v \in \R^n \colon v^T \cdot A \cdot v \ge 0
\end{equation*}

One thing to note here though is that a matrix can in fact have
negative entries and be positive definite, but also have only
positive entries and still not be positive definite.
Using these definitions, we can give the following condition for a
sufficiently smooth function to be convex.

Let $n \in \N$ and $f \in \mathcal{C}^2(\R^n,\R)$. $f$ is
strictly convex or convex if
and only if for all $x \in \R^n$ the Hessian matrix $H_f(x)$ is
positive definite or positive semidefinite respectively.

See <cit.>.

The result is that object of interest is the Hessian matrix of
a given function we want to examine and its definiteness. The
direct approach to check definiteness is complicated, especially
for large systems like the ones we are dealing with in this thesis.
One useful result to simplify this process is the
Gershgorin-Hadamard-theorem that gives conditions
which are easy to check, yielding a positive definite matrix.
To formulate the theorem we must first consider so-called strictly
diagonally dominant matrices where the absolute value of a diagonal
entry is strictly larger than the sum of the absolute values of all
other entries in that row.

Let $n \in \N$ and $A \in \R^{n \times n}$. $A$ is
strictly diagonally dominant if	and only if

\begin{equation*}
\forall i \in \{1,\dots,n\} \colon |a_{ii}| >
\sum_{j \in \{1,\dots,n\} \setminus \{i\}} |a_{ij}|.
\end{equation*}

Let $n \in \N$ and $M \in \R^{n \times n}$ symmetric, diagonally
dominant and $\forall i \in \{1,\dots,n\} \colon a_{ii} > 0$.
Then $M$ is positive definite, non-singular and every eigenvalue of $M$
is positive.

See <cit.>.

Using this theorem it shall be possible to approach even large systems
in regard to positive definiteness. When we show positive definitness
of the Hessian matrix, convexity follows and we have obtained our desired

CHAPTER: NOTATION DIRECTORY

§ CHAPTER <REF>CH:INTRODUCTION: *CH:INTRODUCTION

p2cm p10cm
$\Phi(t)$     1-periodic function;
see Equation (<ref>)

$c_j$     Fourier coefficients;
see Equation (<ref>)

$(\mathcal{F}s)(f)$     Fourier transform;
see Equation (<ref>)

$\left\langle \cdot,\cdot \right\rangle$     standard inner product;
see Equation (<ref>)

$\psi_{j,k}$     wavelet basis function;
see Equation (<ref>)

$\psi$     mother wavelet;
see Equation (<ref>)

$c_{j,k}$     wavelet coefficients;
see Equation (<ref>)

${\| \cdot \|}_2$     $2$-norm based on the standard inner product

§ CHAPTER <REF>CH:B-SPLINES: *CH:B-SPLINES

p2cm p10cm
$\Pi_p$     set of polynomials of order $p$;
see Definition <ref>

$\Sigma_{k,T}$     spline function space;
see Definition <ref>

$\Sigma_{k}$     shorthand notation for $\Sigma_{k,T}$;
see Remark <ref>

$k$     spline function order;
see Definition <ref>

$\ell$     size of spline knot vector;
see Definition <ref>

$T$     spline knot vector;
see Definition <ref>

$\indicator(A)$     indicator function on set $A$;
see Definition <ref>

$B_{i,k,T}$     B-spline function;
see Definition <ref>

$B_{i,k}$     shorthand notation for $B_{i,k,\Delta_k(T)}$;
see Remark <ref>

$\Delta_k(T)$     extended knot vector;
see Definition <ref>

$n$     number of B-spline functions ($k + \ell - 2$);
see Definition <ref>

$\B_{k,T}$     coefficient spline mapping;
see Definition <ref>

$\B_{k}$     shorthand notation for $\B_{k,T}$;
see Remark <ref>

$c_{k,\infty}$     B-spline condition constant;
see Proposition <ref>

§ CHAPTER <REF>CH:EMD_ANALYSIS: *CH:EMD_ANALYSIS

p2cm p10cm
$a$     instantaneous amplitude;
see Defintion <ref>

$\phi$     instantaneous phase;
see Definition <ref>

$\phi'$     instantaneous frequency;
see Definition <ref>

${\mathcal{S}}_{\mu_0,\mu_1,\mu_2}$     set of intrinsic mode function souls (IMFS)
with characteristic $(\mu_0,\mu_1,\mu_2)$;
see Definition <ref>

$\boldsymbol{\mathcal{S}}_{\mu_0,\mu_1,\mu_2}$     set of intrinsic mode spline
function souls (IMSpFS);
see Definition <ref>

$\I[a,\phi]$     intrinsic mode function operator;
see Definition <ref>

$c_1[s](a,\phi)$     canonical EMD cost function;
see Definition <ref>

$\boldsymbol{c}_1[\boldsymbol{s}](\boldsymbol{a},\boldsymbol{\phi})$     canonical spline EMD cost function;
see Definition <ref>

$c_\ell[s](a,\phi)$     leakage factor EMD cost function;
see Definition <ref>

$\boldsymbol{c}_\ell[\boldsymbol{s}](\boldsymbol{a},\boldsymbol{\phi})$     leakage factor spline EMD cost function;
see Definition <ref>

$V^+$     positive cone of cone $V$;
see Definition <ref>

$V'$     dual cone of cone $V$;
see Definition <ref>

$\Lambda(x,\lambda)$     Lagrange function;
see Definition <ref>

$\ul{\Lambda}(\lambda)$     Lagrange dual function;
see Definition <ref>

§ CHAPTER <REF>CH:OSS: *CH:OSS

p2cm p10cm
$A[a]$     instantaneous envelope derivation operator;
see Definition <ref>

$\Omega[\phi]$     inverse square continuous frequency operator;
see Definition <ref>

$\mathcal{D}_{(a,\phi)}$     IMF differential operator;
see Definition <ref>

$\tilde{\mathcal{D}}_{(A,\Omega)}$     modified IMF differential operator;
see Definition <ref>

§ CHAPTER <REF>CH:FUNCSPACEOPT: *CH:FUNCSPACEOPT

p2cm p10cm
$\N$     set of natural numbers $\ge 1$

$\N_0$     set of natural numbers $\ge 0$

$\mathcal{C}^p(A,B)$     set of $p$-continuously differentiable
functions $A \to B$

$\preceq$     partial order on $\mathcal{C}^0(\R,\R)$;
see Proposition <ref>

$| \cdot |$     modulus operator;
see Definition <ref>

${(\cdot)}^p$     power operator;
see Definition <ref>

$D^p$     differentiation operator;
see Definition <ref>

§ CHAPTER <REF>CH:CONVEXITY_THEORY: *CH:CONVEXITY_THEORY

p2cm p10cm
$H_f(x)$     Hessian matrix of $f$ in $x$; see Definition <ref>

CHAPTER: CODE LISTINGS

§ ETHOS TOOLBOX

§.§ ethos.h

§.§ ethos.c

§.§ config.mk

§.§ Makefile

§ EXAMPLES

These programs expect libethos.a and ethos.h somewhere in
the environment. Set the variable  in config.mk
(see Listing <ref>) to point to the
directory containing both files. This in turn will properly
set the - and -flags in the preprocessor flags
and the linker flags  respectively.

§.§ emd.c

§.§ envelope.c

§.§ regop.c

§.§ util.h

§.§ util.c

§.§ config.mk

§.§ Makefile

§ LICENSE

This ISC license applies to all code listings in
Chapter <ref>.

Bourbaki, Nicolas  Topological Vector Spaces: Chapters 1–5,  1  Elements of
Springer-Verlag Berlin Heidelberg, Berlin,
Germany, 1, 2003,

Ciesielski, Krzysztof  Set Theory for the Working Mathematician,  39  London
Mathematical Society Student Texts.
Cambridge University Press, Cambridge, England,
UK, 1, 1997,

Daubechies, Ingrid  Ten Lectures on Wavelets,  61  CBMS-NSF Regional Conference Series
in Applied Mathematics.
Society for Industrial and Applied Mathematics
(SIAM), University City, Philadelphia, PA, USA, 1,

Boor, Carl-Wilhelm
Reinhold de A Practical Guide to Splines,
27  Applied Mathematical Sciences.
Springer-Verlag New York, New York City, NY, USA,

Dierckx, Paul Curve and Surface Fitting with Splines,  1  Numerical Mathematics and Scientific
Clarendon Press, Oxford, England, UK,

Daubechies, Ingrid, Jianfeng
Lu  HauTieng Wu  Synchrosqueezed wavelet transforms: An
empirical mode decomposition-like toolSynchrosqueezed wavelet transforms:
An empirical mode decomposition-like tool.
Applied and Computational Harmonic Analysis,
30(2):243–261, .032011long,

Fourier, Jean Baptiste
Joseph Théorie analytique de la
Ambroise Firmin Didot, père et fils, Paris,
France, 1, 1822.

Galassi, Mark, Jim
Davies, James Theiler,
Brian Gough, Gerard
Jungman, Patrick Alken,
Michael Booth, Fabrice
Rossi  Rhys
Ulerich GNU Scientific
Free Software Foundation, Boston, MA, USA, 2.5, .062018long.

Guo, Baokui, Silong
Peng, Xiyuan Hu  Pengcheng Xu Complex-valued
differential operator-based method for multi-component signal
separationComplex-valued differential operator-based method for
multi-component signal separation.
Signal Processing, 132:66–76,

Haar, Alfréd Zur Theorie der orthogonalen FunktionensystemeZur
Theorie der orthogonalen Funktionensysteme.
Mathematische Annalen, 69:331–371,

Horn, RogerAlan
Johnson Matrix Analysis.
Cambridge University Press, Cambridge, England,
UK, 2,

Huang, Boqiang  Angela Kunoth An optimization based empirical mode decomposition
schemeAn optimization based empirical mode decomposition scheme.
Journal of Computational and Applied Mathematics,
240:174–183, .032013long,
https://dx.doi.org/10.1016/j.cam.2012.07.012, MATA 2012.

Hu, Xiyuan, Silong
Peng  Wen-Liang Hwang EMD revisited: A new understanding of the envelope and
resolving the mode-mixing problem in AM-FM signalsEMD Revisited: A
New Understanding of the Envelope and Resolving the Mode-Mixing Problem in
AM-FM Signals.
IEEE Transactions on Signal Processing,
60(3):1075–1086, .032012long,

Hou, Thomas Yizhao   Zuoqiang Shi Adaptive data analysis via sparse time-frequency
representationAdaptive Data Analysis via Sparse Time-Frequency
Advances in Adaptive Data Analysis, 3(1&2):1–28,

Huang, Norden Eh, Zheng
Shen, Steven R. Long,
Manli C. Wu, Hsing H.
Shih, Quanan Zheng,
Nai-Chyuan Yen, Chi Chao
Tung  Henry
H. Liu The empirical mode decomposition and the hilbert spectrum for nonlinear and
non-stationary time series analysisThe empirical mode decomposition and the
Hilbert spectrum for nonlinear and non-stationary time series analysis.
Proceedings of the Royal Society of London. Series
A: Mathematical, Physical and Engineering Sciences, 454:903–995,

Huang, Chao, Lijun
Yang  Lihua
Yang $\epsilon$-Mono-Component: Its characterization and
construction$\epsilon$-Mono-Component: Its Characterization and
IEEE Transactions on Signal Processing,
63:234–243, .012015long,

ISOIEC JTC 1/SC 22  ISOIEC 9899:1999: Programming Languages — C.
International Organization for Standardization,
Geneva, Switzerland, 2,

Jahn, Johannes Introduction to the Theory of Nonlinear Optimization.
Springer-Verlag Berlin Heidelberg, Berlin,
Germany, 3, 2007,

Küpfmüller, Karl   Gerhard Kohn Theoretische Elektrotechnik und Elektronik: Eine Einführung.
Springer-Lehrbuch. Springer-Verlag Berlin
Heidelberg, Berlin, Germany, 15, 2000, .

Kreutzer, Moritz, Jonas
Thies, Melven Röhrig-Zöllner, Andreas Pieper,
Faisal Shahzad, Martin
Galgon, Achim Basermann, Holger Fehske, Georg Hager  Gerhard Wellein GHOST: Building blocks for high performance sparse
linear algebra on heterogeneous systemsGHOST: Building Blocks for High
Performance Sparse Linear Algebra on Heterogeneous Systems.
International Journal of Parallel Programming,
45(5):1046–1072, .102017long,

Körner, Thomas William  Fourier Analysis.
Cambridge University Press, Cambridge, England,
UK, 1, 1988,

Liu, Yanping, Yong
Wang  Zhen
Wang RBF prediction model based on EMD for forecasting GPS precipitable
water vapor and annual precipitationRBF Prediction Model Based on EMD
for Forecasting GPS Precipitable Water Vapor and Annual Precipitation.
Luo, Xun  (): 2nd International Conference On Systems
Engineering and Modeling (ICSEM-13),  35
Advances in Intelligent Systems
Research,  51–55, Paris, France,
.042013long. Atlantis Press.

Mallat, Stéphane Georges  A Wavelet Tour of Signal Processing.
Academic Press, Boston, MA, USA,
3, 2009, .

Niculescu, ConstantinP.
LarsErik Persson Convex Functions and Their
Applications: A Contemporary Approach,  24
CMS Books in Mathematics.
Springer-Verlag New York, New York City, NY, USA,
1, 2006, .

Peng, Silong  Wen-Liang Hwang Adaptive signal decomposition based on local narrow band
signalsAdaptive Signal Decomposition Based on Local Narrow Band Signals.
IEEE Transactions on Signal Processing,
56(7):2669–2676, .62008long,

Peng, Silong  Wen-Liang Hwang Null space pursuit: An operator-based approach to adaptive
signal separationNull Space Pursuit: An Operator-based Approach to Adaptive
Signal Separation.
IEEE Transactions on Signal Processing,
58(5):2475–2483, .52010long,

Schoenberg, Isaac Jacob  Contributions to the problem of
approximation of equidistant data by analytic functions. Part A. On the
problem of smoothing or graduation. A first class of analytic approximation
formulaeContributions to the problem of approximation of equidistant data
by analytic functions. Part A. On the problem of smoothing or
graduation. A first class of analytic approximation formulae.
Quarterly of Applied Mathematics, 4(1):45–99,

Schoenberg, Isaac Jacob  Contributions to the problem of
approximation of equidistant data by analytic functions. Part B. On the
problem of osculatory interpolation. A second class of analytic
approximation formulaeContributions to the problem of approximation of
equidistant data by analytic functions. Part B. On the problem of
osculatory interpolation. A second class of analytic approximation
Quarterly of Applied Mathematics, 4(2):112–141,

Schaefer, HelmutHeinrich
MichaelP. Wolff Topological Vector Spaces,
3  Graduate Texts in Mathematics.
Springer-Verlag New York, New York City, NY, USA,

Wu, Qin  Sherman Delbert Riemenschneider  Boundary extension and stop criteria for
empirical mode decompositionBoundary Extension and Stop Criteria for
Empirical Mode Decomposition.
Advances in Adaptive Data Analysis, 2(2):157–169,

Yong, Wang, Liu
Yanping  Yang
Jing Signal prediction based on empirical mode decomposition and artificial
neural networksSignal prediction based on empirical mode decomposition and
artificial neural networks.
Geodesy and Geodynamics, 3(1):52–56,

CHAPTER: EIGENSTÄNDIGKEITSERKLÄRUNG

Hiermit versichere ich an Eides statt, daß ich die vorliegende Arbeit
selbstständig und ohne die Benutzung anderer als der angegebenen
Hilfsmittel angefertigt habe. Alle Stellen, die wörtlich oder sinngemäß
aus veröffentlichten und nicht veröffentlichten Schriften entnommen wurden,
sind als solche kenntlich gemacht.
Die Arbeit ist in gleicher oder ähnlicher Form oder auszugsweise im
Rahmen einer anderen Prüfung noch nicht vorgelegt worden. Ich versichere,
daß die eingereichte elektronische Fassung der eingereichten Druckfassung
vollständig entspricht.

Laslo Hunhold