Phase demodulation with iterative Hilbert transform embeddings

TL;DR

This paper introduces an iterative Hilbert transform embedding method for more accurate phase demodulation of complex, fast-modulated signals, surpassing traditional single-application approaches.

Contribution

The authors develop an iterative approach that significantly improves phase demodulation accuracy for complex waveforms, with theoretical convergence analysis for simple cases.

Findings

Iterative method enhances phase accuracy beyond traditional approaches.

Applicable to complex and fast modulated waveforms.

Convergence proven for simple cosine waveforms.

Abstract

We propose an efficient method for demodulation of phase modulated signals via iterated Hilbert transform embeddings. We show that while a usual approach based on one application of the Hilbert transform provides only an approximation to a proper phase, with iterations the accuracy is essentially improved, up to precision limited mainly by the discretization effects. We demonstrate that the method is applicable to arbitrarily complex waveforms, and to modulations fast compared to the basic frequency. Furthermore, we develop a perturbative theory applicable to simple cosine waveforms, showing convergence of the technique.

Equations

$$x(t)=S(\phi (t)).$$

$$S_1(\phi )=\cos (\phi )$$

$$S_2(\phi )=\cos (\phi )-0.7\sin (2\phi )+\cos (3\phi )$$

$$\dot{\vec{y}}=\vec{F}(\vec{y})$$

$$\dot{\vec{y}}=\vec{F}(\vec{y})+\epsilon \vec{P}(\vec{y},t),$$

$$\dot{\phi }=\omega +\epsilon Q(\phi ,t)$$

$$x(t)=x({\vec{y}}_0(\phi (t))).$$

$$\psi =\Psi (\phi ),\Psi (s+2\pi )=\Psi (s)+2\pi ,{\Psi }^{\prime }>0,$$

$$x(t)=\tilde{S}(\psi (t)),\tilde{S}=S({\Psi }^{-1}(\psi (t))$$

$$x(t+L(t))=x(t).$$

$$y(t)=\hat{H}[x(t)]=\frac{1}{\pi }P.V.{\int }_{-\infty }^{\infty }\frac{x(t^{\prime })}{t-t^{\prime }}d{t}^{\prime }.$$

$${\theta }^{(a)}(t)=\text{arg}[x(t)+iy(t)].$$

$${\theta }^{(b)}(T)={\int }_0^{T}dt\sqrt{{\left(\frac{dx}{dt}\right)}^2+{\left(\frac{dy}{dt}\right)}^2}.$$

$$\tilde{\epsilon }(L)=\frac{{\int }_{{\theta }_{min}}^{{\theta }_{max}}(x({\theta }^{\prime }+L)-x({\theta }^{\prime }){)}^{2}d{\theta }^{\prime }}{{\int }_{{\theta }_{min}}^{{\theta }_{max}}{x}^2({\theta }^{\prime })d{\theta }^{\prime }}.$$

$$\hat{L}=\text{argmin }\tilde{\epsilon }(L),\epsilon =\underset{L}{\min }\tilde{\epsilon }(L)=\tilde{\epsilon }(\hat{L}).$$

$$\theta (t)=\hat{P}[x(t)]$$

$${\theta }_{n+1}({\theta }_n)=\hat{P}[x({\theta }_n)]$$

$$\begin{array}{c}y({\theta }_n)=\hat{H}[x({\theta }_n)]=\frac{1}{\pi }\text{PV}{\int }_{-\infty }^{\infty }\frac{x({\theta }_n^{\prime })}{{\theta }_n-{\theta }_n^{\prime }}d{\theta }_n^{\prime },\\ {\theta }_{n+1}=\text{arg }[x({\theta }_n)+iy({\theta }_n)].\end{array}$$

$$\phi (t)=t+ϵq(t),ϵ\ll 1,$$

$$x(t)=$$

$$\hat{H}[x](t)\approx$$

$${\theta }^{(a)}(t)$$

$${\theta }^{(b)}(t)$$

$$+ϵ{\int }_0^t\frac{S^{\prime }(\tau )(q^{\prime }(\tau )S^{\prime }(\tau )+S^{\prime \prime }(\tau )q(\tau ))+\hat{H}[S^{\prime }(\tau )](\hat{H}[q^{\prime }(\tau )S^{\prime }(\tau )]+\hat{H}[q(\tau )S^{\prime \prime }(\tau )])}{\sqrt{(S(\tau ){)}^2+\hat{H}[S(\tau ){]}^2}}d\tau$$

$${\Delta }^{(a)}(t)=ϵ\left(q(t)+\cos (t)\hat{H}[q(t)\sin (t)]-q(t){\sin }^2(t)\right),$$

$${\Delta }^{(b)\prime }(t)=\frac{ϵ}{2}{q}^{\prime }(t)+\frac{ϵ}{2}[\cos (2t)q^{\prime }(t)-\sin (2t)q(t)+2\cos (t)\hat{H}[q^{\prime }(t)\sin (t)]+2\cos (t)\hat{H}[q(t)\cos (t)]].$$

$$\hat{H}[Q_L(t)Q_H(t)]=Q_L(t)\hat{H}[Q_H(t)].$$

$$\hat{H}[q(t)\sin t]=\hat{H}[(\ell (t)+h(t))\sin (t)]=-\ell (t)\cos (t)+\hat{H}[h(t)]\sin (t),$$

$${\Delta }^{(a)}(t)=\frac{ϵ}{2}h(t)+\frac{ϵ}{2}\left(h(t)\cos (2t)+\sin (2t)\hat{H}[h(t)]\right),$$

$${\Delta }^{(b)\prime }(t)=\frac{ϵ}{2}{h}^{\prime }(t)+\frac{ϵ}{2}\left(\hat{H}[h(t)]+\cos (2t)h^{\prime }(t)-\sin (2t)h(t)+\sin (2t)\hat{H}[h^{\prime }(t)]+\cos (2t)\hat{H}[h(t)]\right).$$

$${\mathcal{F}}_{{\Delta }^{(a)}}(\omega )=$$

$${\mathcal{F}}_{{\Delta }^{(b)}}(\omega )=$$

$$\frac{\omega +1}{\omega }(1+\text{sgn}(\omega +2)){\mathcal{F}}_h(\omega +2))].$$

$${\mathcal{F}}_{\Delta }(\omega )=\{\begin{array}{ll}A(\omega ){\mathcal{F}}_h(\omega )+B(\omega ){\mathcal{F}}_h(\omega +2)& \text{if }\omega >1,\\ B(\omega ){\mathcal{F}}_h(\omega +2)+A(\omega ){\mathcal{F}}_h(\omega -2)& \text{if }-1<\omega <1,\\ B(\omega ){\mathcal{F}}_h(\omega )+A(\omega ){\mathcal{F}}_h(\omega -2)& \text{if }\omega <-1,\end{array}$$

Full text

1 Introduction

Demodulation of nearly oscillatory signals is a standard problem in signal processing. In classical communication schemes, one uses amplitude, phase, and frequency modulations to store information in an analog signal. At the receiving side, demodulation is needed to extract this information. In a more general context, when one tries to characterize the dynamics of an oscillatory system based on some observable, modulation carries important information about the system operation. For example, the electrocardiogram is an important signal characterizing the human cardio-vascular system. One of the widely used tools is demodulation of this signal via extracting information about the instantaneous periods (RR intervals), with further analysis of demodulated heart rate variability [1]. In many situations one considers a demodulation task for harmonic carrying waveforms only (see, e.g., [2]), but in fact one quite often needs to demodulate signals with unknown, potentially rather complex waveform (like the electrocardiogram signal).

In this paper we focus on the problem of demodulation of phase modulated signals with arbitrary waveforms and arbitrary phase modulations. This task is essential for characterization of oscillatory processes resulting from forced and/or coupled self-sustained oscillators, e.g. for synchronization analysis [3, 4, 5] and for reconstruction of the phase dynamics [6, 7, 8]. In particular, in Refs. [9, 10] a continuous phase was extracted from the electrocardiogram signal, and further used to build a phase model for the cardio-respiratory interaction.

The Hilbert transform (HT) is a widely used tool in oscillatory data analysis and signal processing. It lies at the heart of the analytic signal approach [11, 12, 13, 14, 15], which provides an effective processing of nearly harmonic signals. The Hilbert transform is also extensively used for the phase demodulation in synchronization studies [16, 3], see also recent review [17]. One uses it to perform a two-dimensional embedding of the observed signal; this embedding then is usually treated as an analytic signal and the phase is extracted as the argument of the complex state (or, more generally, as an angle of rotation with respect to some point on the embedding plane). This approach gives a good approximation, but not the exact phase demodulation. The quality of the analytic signal based phase demodulation is better for slow modulation, although also in this situation, as will be illustrated below, it is not exact.

One can attribute non-exactness of the HT embedding to the fact, that the HT generally mixes phase and amplitude modulations: if one starts with a purely phase modulated signal, then the HT embedded signal will have amplitude modulation as well. In signal processing, workarounds exist [18, 14] that introduce filtering before application of the HT and/or decomposition of the signal into several functions to avoid the mixing.

In this paper we propose a phase demodulation method based on iterative HT embeddings. The first step is the same as in the traditional approach, it yields an approximation, but the signal is not fully demodulated. We use the phase obtained in the first step as a new time variable to perform the next embedding, with better demodulation. The repetition of this procedure provides better and better approximations. We will demonstrate that the iterations result in a very accurate demodulation with errors related to numerical errors in calculating the Hilbert transform integral. Noteworthy, our method is not restricted to simple waveforms like cosine function; we will present its implementation for arbitrary (but smooth) waveforms. Thus, no filtering is involved in the approach we describe.

For completeness, we mention another class of approaches to phase demodulation, where not the phase itself is determnined, but the instanteneous frequency is evaluated via a time-frequency decomposition (see, e.g., [19]). Recently, a variant of this approach – wavelet-based synchrosqueezing – have been suggested [20, 15, 21]. Comparison with these methods is a subject of future research.

The paper is organized as follows. First, in section 2, we elaborate on our definition of the waveform and the phase for a signal, and discuss connections to dynamical systems analysis. Then we present our method in sections 3,4. Analytic results on the convergence of the method, based on a perturbation analysis, are described in section 5. Numerical results, confirming accuracy of the technique for a range of strongly modulated signals and for complex waveforms, are presented in section 6. In section 7, a transformation from a phase variable, obtained in the demodulation procedure, to another phase is discussed. We end with conclusions in section 8.

2 Phase modulated signals

2.1 Definition of a phase modulated signal

The subject of our study are phase-modulated signals. First we define a waveform $S(\varphi)=S(\varphi+2\pi)$ as a $2\pi$-periodic function of its argument. The argument, interpreted as a phase, is a monotonous function of time $\dot{\varphi}>0$. If $\dot{\varphi}\neq\text{const}$, one speaks of a phase modulated signal

[TABLE]

Our task is, given the signal $x(t)$, to reconstruct the phase $\varphi(t)$ and the waveform $S(\varphi)$.

Our method below is not restricted to simple waveforms: $S(\varphi)$ can be a rather complex function with several maxima and minima. Below in the numerical examples and in theoretical considerations we will use

[TABLE]

as a simple waveform, and

[TABLE]

as a complex waveform.

2.2 Phase modulation in driven oscillators

The motivation for this setup comes from the inverse problems for nonlinear dynamical systems. Phase modulated signals naturally appear if a dynamical system with stable self-sustained oscillations (i.e. with a limit cycle) is driven by a weak external force (or interacts with another dynamical system) [22, 23, 16, 24], as we briefly outline below.

Suppose we have an autonomous dynamical system

[TABLE]

with a stable limit cycle $\vec{y}_{0}(t)$. This means that the system generates periodic oscillations and a generic scalar observable $x(\vec{y})$ is a periodic function of time $x(\vec{y}_{0}(t))$. On the limit cycle and in its basin of attraction one can introduce the phase variable $\varphi=\Phi(\vec{y})$, such that the limit cycle is parametrized by this phase $\vec{y}_{0}(\varphi+2\pi)=\vec{y}_{0}(\varphi)$, and the phase grows linearly in time with the natural frequency of the limit cycle $\dot{\varphi}=\omega$.

If there is a weak external force acting on the system

[TABLE]

then, in the first approximation in $\varepsilon$, a phase-modulated dynamics on the limit cycle represents the full oscillation of $\vec{y}(t)$ close to $\vec{y}_{0}(\varphi(t))$, where the phase obeys the equation

[TABLE]

and $Q=\text{grad}\Phi|_{\vec{y}_{0}(\varphi)}\cdot\vec{P}(\vec{y}_{0}(\varphi),t)$ is determined by the phase response curve of the limit cycle and the force $\vec{P}$ [25]. Accordingly, the scalar observable $x$ is the phase-modulated signal as described above:

[TABLE]

Here the waveform is determined by the form of the limit cycle and of the scalar observable $S(\varphi)=x(\vec{y}_{0}(\varphi))$. The inverse problem in the context of the theory above is formulated as follows: From scalar observations $x(t)$ of a driven dynamical system, one seeks to reconstruct the dynamics of the phase $\varphi(t)$. The first step here is the reconstruction of the phase $\varphi(t)$, what is the topic of the present study.

2.3 Phase vs protophase

In this section we discuss whether the decomposition of a signal (1) into the phase and the waveform is unique. The answer is obvious: No. Indeed, one can perform a transformation

[TABLE]

to a new phase variable $\psi$, in terms of which the signal is represented as

[TABLE]

with the new $2\pi$-periodic waveform $\tilde{S}(\psi)=\tilde{S}(\psi+2\pi)$. Generally, there is no particular reason to prefer the representation in terms of the waveform $S$ and the phase $\varphi$ to the representation in terms of the waveform $\tilde{S}$ and of the phase $\psi$.

In applications, “preferable” phase variables can exist. For example, for a dynamical system where the equation for the phase in form (4) is used, the “true” phase has the property of a uniform (on average) rotation, so that the distribution density of this phase does not depend on the phase itself. See more discussion in Refs. [6], where general possible phase variables have been named “protophases”, and a transformation of form (5) from a protophase to the phase having uniform distribution density has been suggested as a step in the phase reconstruction. However, in other applications other conditions on the phase may be imposed. On the other hand, there can be a “preferable” waveform. For example, a waveform with one maximum and one minimum over the period, can be transformed to a cosine waveform. The ambiguity of the phase definition results also in the ambiguity of the instantaneous frequency, if the latter is defined as the time derivative of the phase.

Summarizing, we state that the decomposition of a signal (1) is not unique, and in different methods, described below, different protophases and correspondingly different waveforms will be obtained. We will, however, consider demodulation as successful, if at least one particular protophase and the corresponding waveform can be found. The success of demosulation means that a protophase $\psi(t)$ is found, such that values of $|x(\psi)-x(\psi+2\pi)|$ are very small. Quantitatively, we define the squared integrated error of demodulation in Eq. (9) below.

2.4 Link function approach

There is another formulation of the demodulation problem, based on the link function approach [M. Holschneider, private communication]. Given a modulated signal (1), one looks for a function of time $L(t)$, called link function, such that

[TABLE]

If a protophase $\psi(t)$ of the signal is known, then one can define the link function as $\psi(t+L(t))=\psi(t)+2\pi$, what can be resolved via the inverse function $t(\psi)$ as $L(t)=t(\psi(t)+2\pi)-t$. One can easily see that the link function is independent on the choice of the protophase, and in this sense is invariant. For the dynamical systems described in section 2.2, the link function gives the return time for the Poincaré map from the value of the phase to the next value shifted by $2\pi$. The task of finding a protophase if the link function is known, is, to the best of our knowledge, not solved.

3 Proxi-phase extraction

Our main method is based on an iteration procedure. At each iteration stage, an approximation to the phase is extracted from the phase-modulated signal $x(t)$; we call the result of this extraction “proxi-phase” below. In fact, this extraction step is performed only once in most existing approaches to extract the phase dynamics from data (for example, see [26, 27, 6, 10]).

The usual approach is to perform a two-dimensional embedding $x(t)\to(x(t),y(t))$, so that a phase-modulated signal is represented by “rotations” or “loops” on the plane $(x,y)$. In the field of dynamical system reconstruction, different choices for $y$ have been discussed, two popular are a delayed signal $y=x(t-\tau)$ and the derivative $y=\frac{dx}{dt}$ [28, 29, 30]. However, the mostly stable and reliable results are achieved with the Hilbert transform [3]:

[TABLE]

The resulting trajectory on the plane $(x(t),y(t))$ we call HT-embedding.

Below we will consider two basic ways to extract a proxi-phase from the HT embedding. In the simplest case, where the waveform is close to a $\cos(\varphi)$ function, the HT embedding is close to a circle, surrounding the origin (see Fig. 1(a) below), and one can define the proxi-phase $\theta^{(a)}(t)$ as the angle from the origin (or another point close to the origin) to the position of the signal:

[TABLE]

This definition corresponds to the proxi-phase as the argument of the analytic signal $z(t)=x(t)+iy(t)$. Therefore we will call $\theta^{(a)}$ the analytic proxi-phase. The proxi-phase $\theta^{(a)}$ should be additionally unwrapped, so that it is a motonously growing function of time.

This simple definition of the proxi-phase will not work for a complex waveform with several maxima and minima on the period. In this case the embedding on the plane $(x,y)$ has many loops (see Fig. 1(e) below), and the evaluation of the proxi-phase according to (7) is not possible. Even for the simple waveform $S(\varphi)=\cos(\varphi)$, if the phase modulation is strong and the modulation frequency is high, a situation may occur where the analytic proxi-phase according to (7) is non-monotonous. This makes application of the analytic proxi-phase impossible for such cases.

A more universal proxi-phase definition was suggested in Ref. [6] and adopted for the analysis of an electrocardiogram signal in Ref. [9]. This proxi-phase is based on the curve length, to avoid the mentioned problems. The length-based definition of a proxi-phase $\theta^{(b)}(t)$ is

[TABLE]

For any sufficiently smooth signal, this definition provides a monotonously growing length proxi-phase, even if the embedded trajectory has loops.

A small drawback of this definition is that the proxi-phase $\theta^{(b)}(t)$ is not normalized to intervals of $2\pi$ at each rotation, contrary to the analytic proxi-phase (7). This does not affect further iteration steps to be described below, because the HT is invariant under rescaling of the function argument. At the very end of the procedure, one can perform the normalization by multiplying $\theta^{(b)}$ with a factor $2\pi\hat{L}^{-1}$, if needed. Here $\hat{L}$ is the average period of the signal in terms of the length (8) after the iteration procedure. To define it, we first define the normalized periodicity error of a function $x(\theta)$, with respect to a period guess $L$, as

[TABLE]

In case of the analytic proxi-phase, we can set $L=2\pi$, then expression (9) provides the quality of the phase reconstruction. $\varepsilon=\tilde{\varepsilon}(2\pi)$ only weakly dependent on the protophase definition. If it tends to zero for one protophase, then it tends to zero for other protophases as well.

For the length proxi-phase, the value of $L$ is not known a priori. Here, we obtain an “optimal” period by minimizing the error (9) with respect to $L$:

[TABLE]

This error will be used in figures 4 and 5 below to quantify quality of demodulation.

The proxi-phase extraction can be viewed as a (nonlinear) operator $\hat{P}$, providing for a given signal $x(t)$ a proxi-phase

[TABLE]

with one of the methods above.

4 Iterated Hilbert transform embeddings

Application of the proxi-phase extraction by one of the methods described in section 3 above shows that this proxi-phase is not the phase, because in the function $x(\theta)$ some phase modulation is still present such that it is not exactly periodic.

We suggest to proceed in an iterative manner: We consider the obtained proxi-phase (which we supply with an index that counts the steps in the iteration process) $\theta_{1}$ as a new time variable, and apply the proxi-phase extraction to the signal $x(\theta_{1})$. Namely, we calculate $\theta_{2}(\theta_{1})=\hat{P}[x(\theta_{1})]$, and so on. In other words, we perform iterations

[TABLE]

The original time can be treated as the zero-index proxi-phase $\theta_{0}=t$. To illustrate, we write explicitly the steps for the analytic signal proxi-phase $\theta^{(a)}$:

[TABLE]

In numerical implementation of the iterations, one has to perform the HT for a function defined on a non-uniform grid; see A for the implementation we used.

Our numerical simulations show that this iterative procedure effectively reduces the errors (9), and, with high accuracy, $\theta_{n}$ converges to a correct protophase for large $n$. Just visually one can check this by comparison of embedded plots $(x(\theta_{n}),y(\theta_{n}))$ for different iteration steps $n$, see Fig. 1. If the signal $x(\theta_{n})$ is a periodic function of $\theta_{n}$, then $y(\theta_{n})$ is also a periodic function, and on the embedding plane one observes just a closed curve. Contrary to this, if $\theta_{n}$ is just an approximation and $x(\theta_{n})$ is not periodic but contains rest modulation, then one observes a non-closed trajectory forming a sequence of non-equal loops.

In Fig. 1 one can see that just one step in the phase extraction does not provide a good phase – the embeddings in panels (b),(f) look like quite wide bands. In course of iterations, these bands become narrower, and already at the $10$-th step one observes very narrow lines (panels (d),(h)), what means that the demodulation procedure is successful.

Below we first present a theory, showing that the proposed iteration method indeed converges for a cosine waveform and a weak modulation. Then, in section 6, we will first check predictions of the theory; and then explore convergence of the method for strong modulations, by virtue of the error analysis according to expression (9).

5 Theory

Our strategy here is to show that at each iteration (10) the residue between the true phase and the calculated phase decreases. We represent the modulation of the true phase as

[TABLE]

and restrict our analysis to terms of first order in $\epsilon$. We have fixed the basic frequency of the signal to one, without loss of generality. Thus, the signal and its HT are represented as

[TABLE]

where prime denotes derivative. We calculate both proxi-phases in the first order in $\epsilon$ based on (7) and (8) respectively, and obtain general expressions, valid for arbitrary waveforms $S$ and modulations $q(t)$:

[TABLE]

We are able to give compact tractable expressions for the simplest waveform $S(\varphi)=\cos\varphi$ only. This case is treated below.

The difference $\Delta^{(a,b)}(t)=\varphi(t)-\theta^{(a,b)}(t)$ between the true phase and the reconstructed one should be interpreted as the rest modulation after the first iteration. For the analytic proxi-phase (7) and the length based proxi-phase (8), in the case $S(\varphi)=\cos(\varphi)$, the expressions for the rest modulation, obtained through simple calculations from (13), read

[TABLE]

Here, to get rid of the integral in the expression for $\Delta^{(b)}$, we take derivative with respect to time $t$.

We now use the Bedrosian’s theorem [31], which states that the HT of a product of a low-frequency $Q_{L}(t)$ and a high-frequency $Q_{H}(t)$ signals and can be represented as

[TABLE]

The condition for the Bedrosian’s theorem is that the spectra of $Q_{L}$ and $Q_{H}$ do not overlap: the whole spectrum of $Q_{L}$ lies left to the spectrum of $Q_{H}$.

Thus, to evaluate the terms $\hat{H}[q(t)\sin(t)]$, $\hat{H}[q(t)\cos(t)]$, and $\hat{H}[q^{\prime}(t)\sin(t)]$, we represent $q(t)$ as a sum of two components $q(t)=\ell(t)+h(t)$, where the Fourier spectrum of the slow component $\ell(t)$ contains frequencies $|\omega|<1$, while the Fourier spectrum of the fast component $h(t)$ contains frequencies that obey $|\omega|>1$. Then, we obtain

[TABLE]

and similar for $\hat{H}[q(t)\cos(t)]$ and $\hat{H}[q^{\prime}(t)\sin(t)]$. Substituting this in (14),(15), we see that the component $\ell(t)$ does not contribute to the rest modulation. Only the contribution of the high-frequency component remains:

[TABLE]

At this stage, it is more convenient to operate with the Fourier spectra $\mathcal{F}_{h}(\omega),\mathcal{F}_{\Delta}(\omega)$ of functions $\epsilon h(t),\Delta(t)$. The Hibert transform in the Fourier space has the form $\mathcal{F}(\hat{H}[h(t)])(\omega)=-i\text{\ sgn}(\omega)\mathcal{F}_{h}(\omega)$. Application of the Fourier transform to expressions (16) and (17) yields:

[TABLE]

We now use that $\mathcal{F}_{h}(\omega)=0$ for $|\omega|<1$, what allows us to simplify expressions (18) to

[TABLE]

where

[TABLE]

In (20), the upper sign is used for $\omega>1$ and the lower sign is used for $\omega<-1$.

Relation (19) can be considered as a transformation of the rest phase modulation at a step of our iteration procedure. Denoting the spectrum of the high-frequency modulation $\epsilon h(t)$ at the $n$-th step of the iteration procedure as $\mathcal{F}_{n}$, we can write a general recursion formula for the evolution of this spectrum under iterations:

[TABLE]

This relation is the main result of our weak modulation analysis. We now discuss its meaning for the iteration procedure of demodulation (where we focus, due to the evident symmetry, on positive frequencies only).

The first observation is that in both cases of proxi-phase definition, slow modulation (i.e. with frequencies lower than the basic oscillation frequency) is resolved exactly already in the first iteration. This is a well-known empirical fact that the traditional HT phase extraction (where only the first iteration is used) works well for slow modulation, but is not so good for fast modulation. 2. 2.

The high-frequency components do not disappear immediately, but are, on one hand, damped by factor $A(\omega)$, and on the other hand, produce components shifted to lower frequencies and damped by factor $B(\omega)$. Schematically, this is illustrated in Fig. 2. This scheme shows that for all spectra decaying at infinity, i.e. with $\lim_{\omega\to\infty}|\mathcal{F}_{0}(\omega)|=0$, the modulation eventually disappears, i.e. $\lim_{n\to\infty}|\mathcal{F}_{n}(\omega)|=0$. For some important classes of spectra, one can show that the eigenvalues of the transformation operator (21) are smaller than one (see  B), so that the iteration procedure converges exponentially. In particular, for spectra $\mathcal{F}_{h}(\omega)$ which decay exponentially $\mathcal{F}_{h}(\omega)\sim\exp[-\alpha\omega]$ (what corresponds to a smooth modulation $h(t)$), we have $|\mathcal{F}_{n}|\sim\left(\nicefrac{{1+\exp[-2\alpha]}}{{2}}\right)^{n}\underset{n\to\infty}{\to}0$. 3. 3.

In the analysis above we have not considered components of modulation having frequencies $1,3,5,...$, i.e. harmonics of the basic oscillation frequency. Such components cannot be demodulated, because one cannot distinguish them from the original phase – in fact, a modulation with these frequencies is equivalent to changing the shape of the waveform $S(\varphi)$.

6 Numerical tests

6.1 Testing theoretical relations

Our first numerical example is intended to illustrate the theory of Section 5. Therefore, we use here the simple waveform $S(\varphi)=\cos(\varphi)$, and the modulation is a function of time

[TABLE]

possesing two frequencies $f$ and $g$. We performed iterative HT embeddings as described above, calculated at each step the analytic proxi-phase, and found the spectral components of the differences $\theta^{(a)}_{n}(t)-\varphi(t)$. These spectral components are presented in Fig. 3 for four numerical setups.

• Case (a).

Panel (a) of Fig. 3 shows the case $a=0$, $b=0.2$, $g=3.73$ of a single-harmonic modulation with high frequency. Here the modulation is relatively weak, and we expect relations (19) to hold. One can indeed see that, in agreement with the theory, $\theta^{(a)}_{1}-\varphi$ contains harmonics at frequencies $g$ and $g-2$ with the same amplitude. At the next iteration, the component at $g$ decreases by factor 2, while the component at $g-2$ remains the same; additionally a component with frequency $|g-4|$ appears. Also at iterations 3 and 4, relation (21) is well fulfilled.

• Case (b).

Panel (b) illustrates nonlinear effects at a strong low-frequency single-harmonic modulation. Here $a=1$, $f=0.41$, and $b=0$. One can see that although the frequency of modulation is less than one, this component at $f$ does not disappear completely in $\theta^{(a)}_{1}-\varphi$ as predicted by the linear theory, but is finite (already rather small). In the first iteration higher harmonics $nf$ and the combinational frequencies $nf-2$ appear, and they are relatively large. During next iterations the amplitudes of these components decrease with a factor $\approx 1/2$.

• Case (c).

Panel (c) shows the case of relatively strong two-frequency modulation, where both frequencies are low: $a=b=0.3$, $f=0.41$, $g=0.73$. Due to nonlinearity, many combinational frequencies $kf+lg$, $kl+fg-2$ with integers $k$ and $l$ are exited. It appears that all the new frequencies, like in case (b), are created in the first iteration, further iterations follow roughly the law (21).

• Case (d).

Panel (d) shows, like case (c), a strong two-frequency modulation, but with high basic frequencies $a=b=0.1$, $f=2.41$, $g=3.73$. One can see that at some combinational frequencies, the level even initially increases due to the cascade process of amplitude shifting (see, e.g., components with $f+g-4$), and one needs more iterations to reduce error at these components.

In conclusion, the presented numerical tests show that the theory based on the linear approximation works well for weak modulation, but for strong modulation, there are essential nonlinear effects at the first few iterations. Because the effective modulations become weaker in the course of iterations, the theory better describes higher iterations.

6.2 Accuracy of the iterative procedure

In the following, we present numerical results for the signals $x(t)=S_{1,2}(\varphi(t))$, where the waveforms are given by expressions (2),(3), and the phase modulation is a quasiperiodic function of time

[TABLE]

We have chosen $\omega_{1}=\sqrt{2},\;\omega_{2}=\sqrt{3},\;a=b=0.3$. Parameter $R$ allows us to vary the basic frequencies, keeping the range of the phase modulation (which is quite large, $\text{min}(\dot{\varphi})=1-a(\omega_{1}+\omega_{2})\approx 0.056$, so the first-order theoretical analysis of section 5 is not applicable) constant. In all runs we used discretisation step $\Delta t=0.002$, the length of the time series was about 95 basic periods. Convergence was characterized by evaluating error (9), where, to discard boundar effects, only central 80% of the time series where used.

Our goal is to quantify performance of the iterative procedure described, in dependence on the characteristic time scale defined by parameter $R$, on the properties of the waveform (simple vs complex), and on the choice of proxi-phase in the iterative process (analytic vs length-based).

6.2.1 Simple waveform: comparison of analytic and length-based proxi-phases

For the simple waveform (2), one can use both proxi-phases - the analytic one and the length-based one. The results are shown in Fig. 4. One can see that both proxi-phases ensure demodulation with a very small rest error, if the modulation is slow ($R\leq 1$). For the lowest-frequency case, the length method performs even significantly better at high iterations.

For larger values of $R\geq 2$, the length method appears to be the only applicable one. Here, the analytic proxi-phase fails because at one of the first iterations the analytic proxi-phase becomes a non-monotonous function of time, and iterations break. In contradistinction, the length-based phase always provides a monotonous estimate of the phase. One can see that for the fastest modulation explored ($R=8$), the error decreases not monotonously with $n$ (it increases at $n\gtrsim 10$), but even in this case, the error of the iteration procedure eventually becomes rather small.

6.2.2 Complex waveform

For the complex waveform (3), only the length-based phase can be used. In Fig. 5 we show the errors vs the iteration number $n$, for the same values of the parameter $R$ as in Fig. 4. In all cases the final error is very small, close to $10^{-7}-10^{-8}$. Similarly to the case of the simple waveform, here at the fastest example ($R=8$), the demodulation is not monotonous.

The presented results confirm that the described method provides an effective demodulation of the complex signal. We attribute the rest error to numerical inaccuracies in calculation of the HT, of the length of the trajectory on the plane $(x,y)$, and in the evaluation of the error itself. This is supported by the observation that the asymptotic errors grow significantly if a larger discretization step is used.

It is important to mention, that contrary to the analytic proxi-phase (7) which is $2\pi$-periodic, the length-based proxi-phase (8) does not have a “natural” period: strictly speaking, a proper normalization to period $2\pi$ is only possible at the end stage of the iteration procedure, where the periodicity of $x(\theta^{(b)})$ is well-defined.

6.3 Arbitrary phase modulation

Above we considered phase modulations as quasiperiodic deviations from the linear phase growth. The method is, however, not restricted to such cases. As an example, we present here results for the complex waveform $S_{2}$ (Eq. (3)) with phase modulation defined as $\dot{\varphi}=\Omega(t)$, with a rather arbitrarily chosen instantaneous frequency

[TABLE]

see Fig. 6(b). For the sampling with $\Delta t=0.005$, the demodulation error $\approx 3\cdot 10^{-8}$ was achieved after 20 iterations (see Fig. 6(c)). In Fig. 6(a) the quality of demodulation is illustrated visually, like in Fig. 1.

7 Transformation of a protophase to true phase

As discussed in section 2.3 above, there is an intrinsic ambiguity in the phase demodulation problem: one can perform a transformation of a waveform together with a transformation of an obtained protophase, such that the signal remains the same. Thus, one needs additional conditions to specify a particular phase and a particular waveform uniquely. In the context of the dynamical system applications (section 2.2), a natural condition is that the phase should grow on average uniformly in time. This condition was adopted in [6], where a procedure for the transformation $\theta\to C(\theta)$ from a protophase to the uniformly on average growing phase was described. Roughly speaking, for a wrapped protophase one determines the probability distribution density on the interval $[0,2\pi)$, and performes a transformation such that the new phase has the uniform distribution density. To illustrate this approach, we take an example from section 6.2.2 with $x(t)=S_{2}(\varphi(t))$, where $S_{2}$ is given by Eq. (3) and $\varphi(t)$ is given by Eq. (22) with $R=8$. In Fig. 7 we show the protophase obtained via the length proxi-phase method, and its transformation according to method [6]. The transformed phase is close to the original one, up to an unavoidable, but unimportant, shift. We stress here, that the accuracy of the transformation $\theta\to C(\theta)$ is rather weak, because it relies on the empirical probability density, which one obtains based on finite-size data.

8 Conclusion

In summary, we have proposed a method of phase demodulation based on the iterative HT embeddings. In all tested cases it leads to an effective accurate reconstruction of the waveform and of the phase of the signal, up to unavoidable ambiguity due to protophase variety. Our method can be applied to simple and complex waveforms. In all cases we have found that the proxi-phase extraction based on the length of the embedded trajectory is superior to a traditional definition based on the analytic signal.

We have performed theoretical analysis valid for a weak modulation of the simplest cosine waveform. This analysis shows that while slow modulation is detected exactly already within the first iteration, for fast components the iterations are indeed needed. Furthermore, this analysis yields exponential (in number of iterations) convergence rate for modulation signals which spectra decay exponentially at high frequencies. We have shown also that nonlinear effects in the case of not so weak modulation lead to appearance of combinational frequencies in the modulation spectrum, what nevertheless does not prevent eventual convergence of the procedure.

In this study we considered only regular signals and used a quite high sampling rate. This is mainly to achieve the highest possible performance, because larger discretization steps lead to errors in calculations of the Hilbert transforms. The analysis of different factors affecting the accuracy will be subject of further research. Also applications of the method to practical situations of signals from forced dynamical systems are now under consideration and will be reported elsewhere.

Acknowledgement

A. P. was supported by Russian Science Foundation (Grant No. 17-12-01534). E. G. was supported by Friedrich-Ebert Stiftung. We thank Michael Rosenblum, Michael Feldman, and Mathias Holschneider for their advices and fruitful discussions.