{'fraction_non_alphanumeric': 0.0436506452216196, 'fraction_numerical': 0.009912100243126987, 'mean_word_length': 4.662007623888183, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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.

II. RESULTS
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 t = k∆t 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: := max k |u[k]| 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ω, andz 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: 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: ∆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 =c P 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.  (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 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 A(t) = 1 + r cos νt, u(t) = 0.2 sin √ 2ωt + cos √ 3ωt . 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.

III. DISCUSSION
We proposed an extension of the Hilbert Transform (HT) method for reconstructing the phase from an observed signal. We addressed a limitation of the conventional HT method, that is, the conventional method has been proven to work only for the narrow band signals. We demonstrated that the conventional HT method cannot accurately reconstruct the interpretable phase from phase-modulated signals. Conversely, our method can extract the phase from these types of signals (Figs. 2-4). Furthermore, we have demonstrated the performance of the proposed method by using the simulated data with the amplitude-and phase-modulated signals (Figs. 5 and 6). Consequently, the extended HT method is a promising tool for investigating synchronization phenomena through analyzing oscillatory signals in biological systems.
There are limitations to the proposed method, offering the opportunity for future improvements. First, we assumed that the phase-modulation u(t) is so small that higher-order terms are negligible. In addition, we assumed that the amplitude of the signal is approximately constant. Despite of these assumptions, the numerical results (Figs. 4, 5 and 6) indicate that the proposed method is, to some extent, robust against violations of these assumptions. As shown in Fig. 6, an amplitude modulation of higher frequency can impair the performance of the proposed method. Future work is required to develop a phase reconstruction method that is robust to rapid amplitude modulations. Next, our method assumed a specific form of the observation signals, s(t) = A(t) cos(ωt + u(t)). This assumption can fail in real-world systems. For example, the signals generated by the van der Pol oscillator [21] and signals from ECG measurement [24] are distorted from the sinusoidal oscillations. To address this issue, a transformation method from the protophase to phase [24,44] might improve the performance of the proposed method.
Another future research direction is to apply the proposed method to a real-world dataset. We have shown that the proposed method can accurately reconstruct the high-frequency component compared with the conventional HT method. Therefore, the proposed method should be able to improve the estimation performances of the phase equations from the observed signals in biological systems. Moreover, the proposed method can identify the phase shifts in oscillatory signals more accurately than the conventional method (Fig. 7). It would be an interesting future study to identify the phase shifts in the circadian rhythms using the proposed method.

A. Discrete Hilbert transform
We used the Hilbert Transform (HT) for discrete signals, i.e., we used the discrete Hilbert Transform (HT) to analyze the signals. In this subsection, we describe the definition and properties of the discrete HT. Let X(ω) be the Fourier transform of a continuous signal x(t): The HT of x(t) (Eq. 2) can be written by using its Fourier transform [35] H where sgn denotes the sign function defined as The discrete HT is defined to satisfy the property similar to Eq. (17). We consider the inverse discrete-time Fourier transform of a signal (sequence) x[k]: where X n is the discrete Fourier transform of x[k], ω n = 2nπ/N is the n-th frequency. The discrete HT of x[k] is defined as where S n is for n = 0, 1 for 0 < n ≤ N/2, −1 for N/2 < n ≤ N − 1.

(21)
Again, we assume 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. From this definition (Eqs.20 and 21), we can derive the formula implying that the discrete HT can also reconstruct the phase of the sinusoidal wave, H d (cos(ωk)) = sin(ωk), H d (sin(ωk)) = − cos(ωk).
where 0 < ω < π is the frequency parameter. Next, we introduce Bedrosian's theorem, which states that the HT of the product of a high-pass signal and a lowpass signal with non-overlapping spectra is given by the product of the low-pass signal and the HT of the high-pass signal. Formally, it is written as follows: Theorem (Bedrosian's theorem for the discrete Hilbert transform [35,45]). Let x[k] and y[k] (k = 0, 1, . . . , N − 1) be sequences with their Fourier transform X n and Y n (n = 0, 1, . . . , N − 1), respectively. If there exists an integer 0 < m < N/2 such that B. Analysis of phase-modulation reconstructed via the conventional Hilbert Transform method Here, we analyze the phase-modulation reconstructed via the conventional HT method. The aim of this subsection is to derive a formula (Eq. 9) that characterizes the relationship between the phase-modulation and its reconstruction. Let us assume that we observe a weakly phase-modulated signal where φ[k] =ωk∆t + u[k] is the instantaneous phase at time k∆t,ω is the effective frequency, and u[k] is a small phase-modulation: = max k |u[k]| 1. Note that the amplitude can be set as 1 (A 0 = 1) without loss of generality. The HT method reconstructs the phase via the argument of the analytic signal (Eq. 3). Substituting Eq. (25) To analyze the discrete HT in Eq.
where ω n = 2πn/N is a frequency of the n-th Fourier component, v n [k] = c n e ikωn for n = 0 or n = N/2, c n e ikωn +c n e −ikωn otherwise, 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.

VI. CODE AVAILABILITY
The code of the extended Hilbert transform is available from the corresponding author on reasonable request.