An extended Hilbert transform method for reconstructing the phase from an oscillatory signal

Rhythmic activity is ubiquitous in biological systems from the cellular to organism level. Reconstructing the instantaneous phase is the first step in analyzing the essential mechanism leading to a synchronization state from the observed signals. A popular method of phase reconstruction is based on the Hilbert transform, which can only reconstruct the interpretable phase from a limited class of signals, e.g., narrow band signals. To address this issue, we propose an extended Hilbert transform method that accurately reconstructs the phase from various oscillatory signals. The proposed method is developed by analyzing the reconstruction error of the Hilbert transform method with the aid of Bedrosian’s theorem. We validate the proposed method using synthetic data and show its systematically improved performance compared with the conventional Hilbert transform method with respect to accurately reconstructing the phase. Finally, we demonstrate that the proposed method is potentially useful for detecting the phase shift in an observed signal. The proposed method is expected to facilitate the study of synchronization phenomena from experimental data.


I. INTRODUCTION
Rhythmic activity is ubiquitous in biological systems, including cortical networks in the brain [1,2], human heart and respiratory system [3][4][5], circadian rhythm [6,7], gene expression [8], and animal gait [9][10][11].The phase description approach [12,13] describes the state of a multi-dimensional nonlinear oscillator using a variable called the phase and derives a reduced phase equation from a nonlinear dynamical system.This approach has promoted the understanding of how a population of nonlinear oscillatory elements can synchronize or form a cluster state.Theoretical studies based on the phase equation have been used to investigate the potential mechanisms underlying synchronization phenomena, including mutual coupling among the elements and the common inputs to the elements [14,15].
Fundamental questions in complex systems include how a system in the real-world achieves synchronization and what is the essential mechanism that leads to a synchronization state [16].While the theoretical studies provide potential explanations for the synchronization phenomena, they cannot directly answer these questions.It is essential to reconstruct the instantaneous phase from observed data (e.g., signals or time series) and to infer the phase equation from the reconstructed phase.Many studies have focused on the latter step, that is, they have developed the inference methods for the phase response curve [17][18][19][20][21] and the coupling function [22][23][24][25][26][27][28][29][30] from the phase (various reviews discuss this topic [31,32]).Conversely, a few studies [33] have focused on the former step, i.e., the reconstruction of the instantaneous phase from an observed signal.An accurate phase reconstruction is necessary to study the synchronization phenomena in data because these inference methods assume the perfect phase reconstruction.
There are two primary approaches to reconstructing the instantaneous phase from an oscillatory signal.One simple approach to reconstructing the phase is to use linear interpolation between the subsequent marker events.For example, the phase is defined as 0 or 2π at the time of the action potential (spike) for neuronal oscillators [17][18][19] or a heartbeat [3].This method can accurately reconstruct the phase when the noise level is not very high.However, this method is not applicable to signals without identifiable marker events, such as neuronal spikes.An alternative phase reconstruction approach is to apply the Hilbert transform to the observed signal [16,34,35].An advantage of the Hilbert transform method is that it is applicable even when there is no well-defined marker.Consequently, the Hilbert transform method has been applied to a variety of systems, e.g., the respiratory system in human [3,5], the gene expression in a cell [8], and the human brain activity [29,[36][37][38].The limitation of the Hilbert transform method is that it can reconstruct the physically interpretable phase from a limited class of signals, i.e., the narrow band signals [36,39].Therefore, it is necessary to carefully develop a pre-processing procedure via trial and error, which hinders the application of this method to oscillatory signals.Theoretical studies in signal processing have clarified the mathematical conditions of the signals on which the Hilbert transform method can reconstruct a meaningful phase [36,39,40].However, only a few attempts have been made to develop a method to reconstructing the phase from more general signals.
In this study, we propose an extension of the Hilbert transform method that can reconstruct the interpretable phase from a wider variety of signals.Here, we consider a new class of signals, called "weakly phase-modulated signals," which are an extension of the sinusoidal signals from which the conventional Hilbert transform method can reconstruct the phase.We first demonstrate that this conventional method cannot accurately extract the phase from these signals (Fig. 1).Then, we derive a new algorithm to reconstruct the phase from the phase-modulated signals and empirically show that the proposed method improves the reconstruction performance.This paper is organized as follows.We first review the conventional Hilbert transform method for reconstructing the instantaneous phase from data.In addition, we illustrate the limitation of the conventional method using an example.Second, we present the proposed method for reconstructing the instantaneous phase and examine the computational complexity of the algorithm.Third, we evaluate the performance of the phase reconstruction and compare its performance with that of the conventional method.Finally, we conclude this study and discuss future directions.

A. Estimating the instantaneous phase from an oscillatory signal
A standard method for reconstructing the instantaneous phase from an oscillatory signal is based on the Hilbert Transform (HT) [16,34,35].This method calculates the phase from the analytic signal, defined as where y(t) and H [y(t)] are the observed signal and its HT where P.V. refers to the Cauchy principal value.The HT method reconstructs the instantaneous phase by the argument of the analytic signal It is well-known that the HT method can reconstruct the interpretable phase from a particular class of signals.Let us consider the sinusoidal signal where ω is the effective frequency, and φ 0 is the initial phase.The HT method can perfectly reconstruct the interpretable phase from the signal: φ H (t) = ωt + φ 0 .Furthermore, it is possible to extend this result to signals with slow amplitude modulation where the amplitude A L (t) is the low-pass-filtered signal whose Fourier coefficients of the frequency higher than the effective frequency (f > ω) vanish.It can be shown [41] that the HT method can perfectly reconstruct the phase: φ H (t) = ωt + φ 0 .However, the HT method can only reconstruct the interpretable phase from a particular class of signals, i.e., the narrow band signals [36,39].
In this study, we extend the HT method for another type of signal, which we call "weakly phase-modulated signals" where u(t) is the phase-modulation from the sinusoidal signal.We set φ 0 = 0 without loss of generality by shifting the time.We applied the HT method to a phase-modulated signal (Fig. 1(a)).Figure 1(b) demonstrates that the HT method can accurately track the linear trend ωt and estimate the effective frequency ω even from a phase-modulated signal.Note that this method (Fig. 1(c), red) cannot accurately reconstruct the phase-modulation φ(t)− ωt.Then, we analyzed the power spectrum of the phase-modulation φ(t)− ωt to investigate the effect of the HT method.Figure 1(d) compares the power spectrum of the phase-modulation reconstructed using the HT method with that of the true phasemodulation.We plotted the frequency range of ω < f < 2ω because the phase-modulation is given by the sum of two sinusoidal functions (sin √ 2ωt and cos √ 3ωt) in this example.The result indicates that the HT method behaves like a low-pass filter, that is, it suppresses the spectral density of the peak frequencies (f = √ 2ω, √ 3ω).Motivated by this observation, we investigate how the HT method changes the power spectrum in the following subsection.Then we extend the HT method to reconstruct the instantaneous phase from an oscillatory signal; the proposed method preserves the power spectrum of the phase-modulation u(t).
Note that it is critical to reconstruct the phase-modulation u(t) accurately to study the synchronization mechanism [16], even though the modulation is small (Fig. 1(b), (c)).Many methods for inferring the phase coupling function rely on the assumption that the phase has been perfectly reconstructed; consequently, the bias in the phase reconstruction may induce a serious effects on the inference results.

B. Proposed method
As we observed in Fig. 1, the conventional HT method cannot reconstruct an interpretable phase from phasemodulated signals.In this subsection, we extend the HT method to include the phase-modulated signals (Eq. 6).
Let us assume that the signal is sampled at N time steps with a constant interval ∆t.We consider a phase modulated signal (Eq.6) sampled at time where A 0 is the amplitude and φ[k] := φ(k∆t) is the instantaneous phase at time t = k∆t.
We can analyze the effect of the phase-modulation on the phase reconstructed via the HT method with the aid of Bedrosian's theorem.The true phase-modulation u[k] := u(k∆t) and its reconstruction via the HT method u H [k] := φ H (k∆t) − ωk∆t can be represented as Fourier series: where ω n = 2nπ/N , and c n and c H n are given by the discrete Fourier transform of u[k] and u H [k], respectively.Assuming that the phase-modulation is small: 1, we can derive a formula that clarifies the relation between the Fourier coefficients (c n and c H n ) (see Methods for the derivation), where m := ωN ∆t/2π is a frequency index corresponding to the effective frequency ω, and z denotes the complex conjugate of a complex number z.In addition, it is assumed that the observation is sufficiently long to satisfy N > 4m, and that the number of data points N is even.If this number is odd, the term N/2 should be replaced with (N − 1)/2.This result (Eq.9) illustrates the effect of the HT method on the phase-modulation in the frequency domain.Eq. (9) shows that the phase reconstructed by the conventional HT method is inconsistent with the true phase for phase-modulated signals.This is because the Fourier coefficients reconstructed via the HT method c H n are not equal to those of the true phase-modulation c n .In addition, the result (Eq.9) implies that the HT method acts as a low-pass-like filter to the phase-modulation u(t).Let us consider the phase-modulated signal with a single frequency component j: Algorithm 1: Proposed method for reconstructing the instantaneous phase.
1: Calculate the initial guess of the phase φ H [k] by using the Hilbert transform method (Eq.3).2: Estimate the effective frequency ω = (φ where α is a non-zero complex value.If the phase-modulation frequency is lower than the effective frequency: j < m, the HT method perfectly reconstructs the true phase, i.e., c H n = c n for all the n.Conversely, when the phasemodulation frequency is higher than the effective frequency: j > m, the amplitude of the reconstructed phasemodulation is half of that of the true phase-modulation, i.e., c H n = c n /2 for all the n.Indeed, Fig. 1(d) shows that the Fourier coefficient of the reconstructed phase-modulation c H n is smaller than the true modulation c n near the dominant Fourier modes ( √ 2ω and √ 3ω).We can extend the HT method to accommodate phase-modulated signals.In the following, we describe the proposed method, which consists of five steps (Algorithm 1).First, we calculate the initial guess of the phase φ H [k] by using the conventional HT method (Eq.3).The Gibbs phenomenon dramatically impairs the phase reconstruction of the HT method when there is a large discrepancy between the values of the first and last point [38].To mitigate this phenomenon, we extract the peaks from the signal and restricted the analysis to be from the first peak to the last one before applying the HT method.Second, we estimate the effective frequency ω from the initial guess: , where T = (N − 1)∆t is the observation duration.Third, we calculate the discrete Fourier transform of the initial guess {c H n } (n = 1, 2, . . ., N ).Fourth, we correct the Fourier coefficient by inverting Eq. ( 9), where c P n is the corrected Fourier coefficient.The remaining coefficients c P n (N/2 < n ≤ N − 1) are calculated by using the formula c P N −n = cP n that reflects the fact that the phase-modulation u[k] is the real signal.Finally, in the fifth step, we reconstruct the phase-modulation by calculating the inverse Fourier transform of {c P n } and smoothing the phase signal.We identify outliers in the reconstructed phase using Median Absolute Deviation criteria [42] and replace the outliers with a linear interpolation of the nearest neighbors.
Finally, we compare the computational complexities of the conventional HT method and the proposed method for reconstructing the phase of a signal.Computational complexity, that is, the dependency of the computational time on the data length, is critical when analyzing a long signal.Let us denote the number of data points of the signal as N .The computational complexity of the HT method is O(N log N ), because we calculate the discrete Hilbert Transform (HT) by using the discrete Fourier transform (See Method).Next, we evaluate the computational complexity of the proposed method.First, the proposed method computes the HT: O(N log N ) (Step 1 in Algorithm 1).Next, the effective frequency is calculated: O(1) (Step 2).Then, the discrete Fourier transform is computed: O(N log N ) (Step 3) and the coefficients of the Fourier transform are corrected: O(N ) (Step 4).Finally, the method reconstructs the phase-modulation by calculating the inverse Fourier transform: O(N log N ) and smoothing it: O(N log N ) (Step 5).Therefore, the computational complexity of the proposed algorithm is O(N log N ), which is comparable to the conventional method.

C. Reconstruction performance of the proposed method
Here, we examine whether the proposed method can accurately reconstruct the instantaneous phase from an observed signal.First, we consider an oscillatory signal with a constant amplitude x(t) = cos (ωt + u(t)) , (11) where ω is the effective frequency and u(t) is the phase-modulation.The sampling time interval and the duration of the simulation are ∆t = 0.01 and T = 200, respectively, unless otherwise stated.We evaluated the performance of the phase reconstruction by analyzing the synthetic data based on two types of phase-modulated signals.The first signal is a quasi-periodic phase-modulation, where b is the amplitude of the phase-modulation.The second signal is the Ornstein-Uhalenbeck (OU) type phasemodulation where η(t) is the Gaussian white noise with zero mean and unit variance.
We applied the proposed method to a signal with quasi-periodic phase-modulation (Fig. 2(a)).Figure2(b) compares the phase reconstructed via the proposed method (blue) with that reconstructed via the conventional HT method (red).While the proposed method accurately reconstructs the phase-modulation, the conventional method cannot reconstruct it.In addition, we compared the power spectrum of the phase-modulation with that of the reconstructed phase-modulations (Fig. 2(c)).We found that the proposed method can reconstruct a phase-modulation whose power spectrum is consistent with the true power spectrum.Next, we applied the proposed method to a signal with the OU-type phase-modulation (Fig. 3(a)).Similar to the case of the quasi-periodic modulation, the proposed method can reconstruct the phase-modulation (Fig. 3(b)) and its power spectrum (Fig. 3(c)) accurately.While the conventional HT method can track the slow trend of the phase fluctuation, it cannot accurately reconstruct the phase-modulation.
Furthermore, we examined whether the proposed method can reconstruct the phase given a larger phase-modulation.We quantified the phase reconstruction performance based on the mean squared error of the phase-modulation.Figure 4(a) shows that the proposed method consistently performs better than the conventional HT method for a signal with a quasi-periodic phase-modulation for a range of the phase-modulation amplitude b.Nevertheless, the error of the proposed method increases with increasing phase-modulation amplitude.The performance deterioration may be due to the nonlinear effects, i.e., O 2 , which we neglected in the derivation of the method.Similarly, the proposed method consistently performs better than the HT method for a signal with an OU-type phase-modulation even when the phase-modulation is not small (Fig. 4 (b)).Even though we assumed a small phase-modulation to derive the proposed method, the results (Fig. 4) suggests that the proposed method provides better performance than the HT method for signals with moderate phase-modulations.Finally, we examined whether the proposed method can reconstruct the phase from a signal with the amplitude and phase modulations.We consider a signal whose amplitude is modulated by a sum of sinusoidal functions where ω is the effective frequency, Figure 5 demonstrates that while the proposed method can accurately reconstruct the instantaneous phase from the amplitude and phase modulated signal, the conventional HT method cannot.This result can be understood as follows.Bedrosian's theorem [41] states that the HT of the product of a high-pass and a low-pass signal with non-overlapping spectra is equal to that of the low-pass signal, and the HT of the high-pass signal.This theorem implies that the slow (or low-pass filtered) amplitude modulation (Eq.5) will not impair the phase reconstruction by the conventional HT method.Thus, it is natural to expect that the proposed method works even when we observe the weakly phase-modulated signals.
Furthermore, we examined whether the proposed method is robust against the fast amplitude modulation.We consider an amplitude and phase modulated signal (14) with Figure 6 shows the dependence of the reconstruction error on the amplitude r and frequency ν of the amplitude modulation.The error does not depend on the amplitude, which indicates that the proposed method works even for signals with moderate amplitude modulation (Fig. 6(a) ).While the error does not depend on the frequency ν in the range of ν < ω (Fig. 6(b) ), it increases when the frequency becomes larger than the effective frequency ω.
Nevertheless, the error of the proposed method is smaller than the HT method.Overall, these results suggest that the proposed method improves the phase reconstruction even for the signals with amplitude modulation.
D. Detecting a phase shift from an observed signal Biological oscillatory systems often exhibit "phase shifts," that is, a rapid change in the phase of a rhythm.For example, the phase of a circadian rhythm can change as a result of light exposure [43].It would be useful to develop a method for detecting the phase shifts in oscillatory signals.We examined whether the proposed method is potentially useful for detecting the phase shifts in data.As a minimal model, we consider a single oscillator exhibiting a phase shift: where the interval [T c , T c + ∆T ] represents the change period, i.e., the period in which the phase of the oscillator shifts with a frequency ω2 , and η(t) is the Gaussian white noise with zero mean and unit variance.The synthetic data were simulated with a total duration T = 10 and a change time T c = 5.0.
Figure 7(a) shows an observed signal before and after the change period, with the frequency increasing after t = 5.0.The proposed method can accurately track the change in the phase-modulation induced by the phase shift (Fig. 7(b): Blue).Conversely, it is difficult for the conventional HT method to infer the change period because of the smoothing effect in the reconstructed phase (Fig. 7(b): Red).Furthermore, we compared the phase reconstruction error of the proposed method with that of the conventional HT method.The proposed method achieved a smaller error than the HT method across a range of phase shift amplitude (ω 2 − ω1 )∆ T (Fig. 7(c)).Similar to the previous result concerning the reconstruction error (Fig. 4), the errors of these methods increase with increasing phase shift amplitude.In addition, we examined the dependency of the error on the duration of the phase shift ∆ T when the phase shift amplitude is fixed.The error of the HT method increases as the shift duration decreases (Fig. 7(d): Red).Conversely, the error of the proposed method is small even for a signal with a small duration ∆ T (Fig. 7(d): Blue).This result suggests that the proposed method is more suitable for detecting rapid phase shifts than the HT method.
is the n-th frequency component.For the derivation of Eq. ( 28), we used the formula c N −n = cn , which reflects the fact that the phase-modulation is a real signal.Here, the number of data points N is assumed to be even.If this number is odd, the term N/2 should be replaced with (N − 1)/2.Due to the linearity of the HT, the reconstructed phase-modulation can be written as We can further calculate Eq.( 29) for each term f (v n [k]) by dividing three cases based on a frequency index m = ωN ∆t/2π that corresponds to the effective frequency ω.
Case 1: Low frequency modulation: 0 ≤ n ≤ m − 1.Using Bedrosian's theorem (Eq.24) and Eq.( 22), we have Substituting Eq.( 30) into Eq.(27), we obtain Case 2: Middle frequency modulation: n = m.Using Bedrosian's theorem, we have where the observation is assumed to be sufficiently long: N > 4m.Substituting Eq.( 32) into Eq.(27), we obtain where c.c denotes the complex conjugate of the terms in the curly brackets.From Eq.( 36), we finally obtain Eq.( 9) by neglecting the higher order terms O( 2 ), which are expected to be small when the phase-modulation is weak.Note that we can also obtain a similar formula between c n and c H n (Eq.9) even when the amplitude A(t) is changes slowly, i.e., the amplitude is a low-pass signal whose spectrum does not overlap with the spectra of cos φ(t).

C. Phase reconstruction based on the discrete Hilbert transform
Here, we describe the conventional HT method for reconstructing the phase φ H [k] from a discrete oscillatory signal x[k].First, we pre-process the signal in order to mitigate the Gibbs phenomenon: we detect the first and the last peak points of the signal, and delete all of the data points prior to the first peak or after the last peak.Next, we calculate the Fourier transform of x[k] to obtain X n .Then, we calculate the discrete HT H d (x[k]) according to Eq. (20).Finally, we reconstruct the phase by calculating the argument of the analytic signal φ H [k] := arg (x[k] + iH d (x[k])).

V. DATA AVAILABILITY
The datasets used during this study are available from the corresponding author on reasonable request.

2 FIG. 2 . 2 FIG. 3 .FIG. 4 .
FIG.2.Reconstructing the instantaneous phase by the proposed method: a signal with the quasi-periodic phase-modulation.(a): Observed signal x(t) given by Eqs.(11) and(12).(b): phase-modulation u(t).(c): Power spectrum of the phasemodulation.Dashed lines represent the true phase-modulation u(t) in (b) and its power spectrum in (c).Red and blue lines represent the reconstructions by the conventional HT method and the proposed method, respectively.The dotted vertical lines in (c) represent the dominant frequencies of the phase-modulation: √ 2ω and √ 3ω.Parameters are ω = 2π and b = 0.2.

FIG. 5 .FIG. 6 .
FIG.5.Reconstructing the instantaneous phase by the proposed method: a signal with amplitude and phase-modulation.(a) Observed signal x(t) (Eq.14).(b) phase-modulation u(t).The dashed line represents the true phase-modulation u(t), and the red and blue line represents its reconstruction by the HT method and the proposed method, respectively.

FIG. 7 .
FIG. 7. Detecting the phase shift from an oscillatory signal.(a): Observed signal x(t) (Eq.15).(b): phase-modulation u(t).(c, d): Dependence of the phase reconstruction error on the phase shift amplitude in (c) and the phase shift duration in (d).We plotted the mean and standard deviation of the errors calculated from 100 trials in (c) and (d).Parameters were set as the effective frequency ω1 = 2π, the noise variance σ = 0.1, the shift amplitude (ω2 − ω1)∆T = 0.5 in (a), (b), and (d), and the the shift duration ∆T = 0.2 in (a), (b), and (c).