Evaluating Physiological Dynamics via Synchrosqueezing: Prediction of Ventilator Weaning

Oscillatory phenomena abound in many types of signals. Identifying the individual oscillatory components that constitute an observed biological signal leads to profound understanding about the biological system. The instantaneous frequency (IF), the amplitude modulation (AM), and their temporal variability are widely used to describe these oscillatory phenomena. In addition, the shape of the oscillatory pattern, repeated in time for an oscillatory component, is also an important characteristic that can be parametrized appropriately. These parameters can be viewed as phenomenological surrogates for the hidden dynamics of the biological system. To estimate jointly the IF, AM, and shape, this paper applies a novel and robust time-frequency analysis tool, referred to as the synchrosqueezing transform (SST). The usefulness of the model and SST are shown directly in predicting the clinical outcome of ventilator weaning. Compared with traditional respiration parameters, the breath-to-breath variability has been reported to be a better predictor of the outcome of the weaning procedure. So far, however, all these indices normally require at least \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\hbox{20}$\end{document}min of data acquisition to ensure predictive power. Moreover, the robustness of these indices to the inevitable noise is rarely discussed. We find that based on the proposed model, SST and only \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\hbox{3}$\end{document} min of respiration data, the ROC area under curve of the prediction accuracy is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\hbox{0.76}$\end{document}. The high predictive power that is achieved in the weaning problem, despite a shorter evaluation period, and the stability to noise suggest that other similar kinds of signal may likewise benefit from the proposed model and SST.


I. INTRODUCTION
B IOLOGICAL signals can contain a wealth of information.
In particular, to evaluate a person's physiological condition we can extract information from a variety of biological signals such as ECG, respiratory signals, blood pressure, and circadian rhythm [1]- [8]. In some cases, this information is easy to read and interpret, in others, it is less accessible, and more sophisticated approaches are needed to extract the information. Many of the measured signals are oscillatory, and one particular and common technique is to focus on the oscillatory features. The fundamental quantity describing the oscillation is its period, which is defined to be the time needed for an observer to observe a "complete and intact oscillation"; this can be expressed by the frequency, which qualitatively is the inverse of the period, that is, it gives the number of oscillations per unit time period. For example, the periods of the ECG signal, the respiratory signal, and the circadian rhythm are about 1 s, 5 s, and 24 h, corresponding to 1 Hz, 0.2 Hz, and 11.6 mHz, respectively [2]- [5]. In recent years, growing evidence suggests that information extracted from biological signals with oscillatory features has diagnostic and prognostic value in various diseases [1], [2], [4]- [6], [9], [10].
Mathematically, frequency analysis is typically studied via the Fourier transform when the signal can be assumed to be stationary. However, it has been long observed that the stationarity assumption is too restrictive for physiological signals, and more information of the physiological system can be extracted if one allows time dependence in the frequency or period. For example, the variability of the time intervals between sequential heart beats, or heart rate variability (HRV) [2], [4], observed in ECG signals, and the variability of the time intervals between sequential breath intakes, or respiratory rate variability (RRV) [3], [10]- [16] are well known to be related to physiological dynamics. Accurate extraction of this type of time-varying information improves diagnostic accuracy and treatment quality [3], [4], [9], [10] and much effort has been put on this direction. In general, time-varying frequency is not measured directly, but is inferred from the behavior, in time, of the oscillation-to-oscillation intervals. A well-known example is the analysis of R-peak to R-peak intervals (RRI) to reveal HRV information [4], [9], [17]. Many techniques have been introduced, including spectral methods and nonlinear dynamical analysis, such as Poincare map, entropy analysis, fractal analysis, to analyze these oscillation-tooscillation intervals [10]- [17].
These established analysis techniques have at least the following three limitations, however, in their use for the study of Fig. 1. Typical recorded respiratory signal from an intensive care unit patient with support from a ventilator. The arrows with a mark R indicate times in the signal where it would be difficult to identify a basic oscillation if only peaks were taken into consideration; the arrow with a mark B indicates a brief machine recalibration; the arrow with a mark G indicates an invalid respiratory trigger. This patient succeeded in ventilator weaning. physiological signals. A first limitation is the (relatively) large number of oscillations that must be observed. For example, for analysis of respiratory signals, at least 300 and 100-1000 oscillations are needed for methods that use a Poincare map [10], [14], [15] or approximate entropy [12], [18], respectively. It is feasible to collect many points for the ECG signal, since this requires only a continuous recording for about 5 min, a reasonable length for a bedside observation. For signals that function on larger time scales than the ECG signal, the story is different. For example, for respiratory signals, we normally need at least 20 min or longer to collect the necessary amount of data, which is usually difficult in certain clinical settings such as patients in intensive care unit [10], [12], [14]- [16] or newborns [18]. The situation would be even worse if we wanted to analyze larger time-scale physiological signals such as circadian rhythm [5], [7].
The second limitation is that it is not always straightforward to determine the oscillation-to-oscillation time series from the given oscillatory signal. Recall that this determination depends on being able to isolate individual oscillations, which requires the identification of the complete repeating basic pattern and, possibly, landmarks within the pattern. Given a suitable definition of this repeating pattern, the oscillation-to-oscillation time series is determined by finding the landmarks for each oscillation. For example, for the ECG signal, the pattern is related to the electrophysiological activity of a normal heart beat and the landmark is defined to be the R peak, and the RRI-time series is based on the R-peak detection [4], [9], [17]. For other physiological signals, it is not always easy to define a basic "oscillation" or a "landmark," even for healthy subjects. For example, although we can provide a definition for respiratory signals, in practice determining the landmarks is not easy, specifically when there is invalid or doubly triggered respiration (see Fig. 1). This difficulty can be mitigated to some extent by noise removal or noise reduction algorithms [17], but even then no reliable determination of the "true" landmarks can be guaranteed. Sometimes it is hard to even provide a universally accepted definition of a landmark, e.g., for the electroesophagographic signal.
The third limitation is the overreduction of the information inside the physiological signals by retaining only the oscillationto-oscillation time series. For example, respiratory signal information is hidden in the time varying amplitude of the ECG signal [19], which is lost in the RRI time series. Indeed, the basic pattern of the oscillation of the ECG signal itself also varies according to time due to the cardiac axis rotation induced by the respiration and other effects.
To address these limitations, we propose in this paper a descriptive model featuring a more fine-grained description of the oscillatory physiological signal than given by the time intervals between sequential oscillations. The model is characterized by the wave-shape function, which is defined to replace the definition of an oscillation and the landmark, the instantaneous frequency, the variability of which is defined to be a proxy for the physiological dynamics, and the amplitude modulation, which is aimed to capture more physiological information. The companion algorithm, referred to as the synchrosqueezing transform (SST), is introduced to provide an accurate estimation of the instantaneous frequency and the amplitude modulation [20], [21].
We recently reported [14], [15] that small variabilities of respiratory parameters including rate and flow are associated with a high incidence of weaning failure in intensive care unit patients, and these variabilities may serve as reliable predictors for weaning patients from mechanical ventilation. As an testbed of the proposed model and algorithm, we analyze the respiratory signals collected from patients in one of our recent studies [15]. We show that the variation of the recorded respiratory signal, monitored for as brief a period as 3 min, together with the tidal volume information, can be used to define a weaning index (WIN) that predicts weaning process outcome with a success rate quantitatively as high as with the area under curve (AUC) of 0.76 when analyzed by the receiver operation characteristic (ROC).

A. Model
In this section, we provide a phenomenological model describing general oscillatory physiological signals. A physiological system is closely linked with a variety of other physiological systems that interact in complex ways; it is well known that, for example, chemical set points and metabolic demand play a role in respiration patterns [3], [15]. Our treatment of these signals will be purely phenomenological; that is, the parameters and indices we will derive from observations of the physiological signal will be based solely on these signals themselves, and not on explicit, quantitative models of the underlying mechanisms. We will show by example in the next section that these parameters and indices in the model contain information that can provide insight into the functioning of the underlying physiology.
The major characteristic pattern of an oscillatory physiological signal is that it is a (fairly) periodic phenomenon; we therefore model it (without noise) as where s(·) is a continuously differentiable periodic function we call the wave shape function, s(t + 2π) = s(t) for all t; it is an oscillating function that satisfies some mild technical conditions [22] (Note that to make the discussion clear, we assume that the signal has just one component, unlike [22], where superpositions of several components were considered.) We call the derivative φ (·) of the phase function φ(t) the instantaneous frequency of f ; we require it to be positive, but it need not be constant; we allow it to vary in time, as long as the variations are slight from one period to the next, i.e., |φ (t)| ≤ φ (t) for all t, where is some small, preassigned number. Likewise, the amplitude A(t) should be positive, but is allowed to vary slightly as well, i.e., |A (t)| ≤ φ (t). In summary, we have the following conditions for all t ∈ R: (2) For the identifiability problem raised in this model and hence the terminologies instantaneous frequency and amplitude modulation, we refer the reader to [21] for the discussion.
In reality, the recorded physiological signal g(t) is contaminated by noise or measurement error, and we model the deviated physiological signal as [21], [23] where Φ(t) is a generalized stationary random process with finite variance and σ(t) is a slowly varying smooth function which capture the heteroscedasticity of the error. Although the possible noise appearing in the medical signal is versatile, our model covers a large portion of it, for example, the time-dependent noise, the Poisson noise and even "slightly" nonstationary noise. Note that the commonly used Gaussian white noise model is when σ(t) = 1 and Φ(t) is the derivative of the Brownian motion.

B. Methodology
Given the model (1) for the oscillatory physiological signal f (t), we want to capture the time-varying quantities of the signal, including the instantaneous frequency φ (t) and the amplitude modulation A(t), when the signal is contaminated by noise as the model (3). It is well known that the continuous wavelet transform (CWT) and short time Fourier transform (STFT) provide profound information about these time-varying quantities, in particular the instantaneous frequency, but accurate extraction via these methods remains an issue, even after many years of research. Reallocation is a technique widely employed in order to get accurate estimates of the instantaneous frequency [24]- [26]. In general, these methods reallocate the wavelet coefficients or STFT coefficients according to some "regrouping" rule, making it possible to read the instantaneous frequency from the resulting time-frequency plane representation.
Synchrosqueezing transform (SST) is a recently introduced novel reallocation technique introduced in [27] in order to analyze speech signals; it was theoretically proved to enjoy several nice properties, useful in our analysis of [20]- [23]. Specifically, the instantaneous frequency φ (t) and the amplitude modulation A(t) can be accurately estimated and the estimation does not depend on whether or not the wave-shape function is a cosine [22]; moreover, the SST is robust to several different types of noise, like the white or colored Gaussian noise [23] or almost stationary generalized random process [21]. Furthermore, the analysis result is adaptive to the data in the sense that the error is dependent only on the first three moments of the chosen mother wavelet and its derivative instead of the profile of the mother wavelet. We summarize the reallocation technique and SST in Appendix A and its numerical implementation in Appendix B. 1 For the sake of convenience, in what follows we shall use the acronyms SSTIF (for SynchroSqueezing Transform-derived Instantaneous Frequency) to refer to the SST-estimated instantaneous frequency, and SSTAM (SynchroSqueezing Transformderived Amplitude Modulation) to refer to the SST-estimated amplitude modulation.

III. TESTBED: VENTILATOR WEANING PROBLEM
Making a weaning decision for a patient on a ventilator is clinically an important issue. The RRV has proved to be helpful in predicting the outcome of weaning intubated patients from the ventilator [14]- [16]. Extended intubation has many negative side effects, such as an increased risk for infection [28], [29]; ideally physicians seek to extubate as soon as medically possible. Yet, weaning too early carries risk as well: reintubation leads to stress to patients or a higher mortality rate [28], [29]). It is thus important to decide accurately when patients can be weaned from the ventilator. To increase the weaning success, the present common practice is to conduct spontaneous breathing trials before a weaning attempt; the final weaning decision is based on the patient's performance during spontaneous breathing trials, characterized by parameters derived directly from the respiration signal, such as the rapid shallow breathing index (RSBI) [30], and subjective evaluation by the clinician. Unfortunately, weaning failure still occurs in a significant percentage of patients who are judged ready-to-wean [28], [29]. Recently, several RRV-based predictors have been proposed to increase the rate of success weaning in this context; these newer predictors are reported to have a higher accuracy than RSBI in predicting weaning success or failure [10], [14]- [16]. Since the oscillation-to-oscillation time series is the focus of RRV analysis, observation of the respiration signal usually takes at least 20 min to guarantee highly significant prediction accuracy [13]- [15]. In addition, these more accurate RRV-based predictors rely on accurate timing of breath intake "peaks", and are thus likely subject to stability issues caused by the inevitable noise.
To mitigate these limitations, we analyze the respiratory signal by the approach proposed in Section II. First, we model the respiratory signal R(t) (without noise) as where s(·), A(t), and φ(t) satisfy (2). Fig. 2 illustrates the role of the different constituents of f(t) modeling the respiratory signal.
In reality, the recorded respiration signal Resp(t) is noisy, and we model it as where σ(t)Φ(t) satisfies the same conditions as that of (3). We apply the methodology described in Section II to directly  analyze Resp(t) signals measured continuously during clinically feasible time intervals, far shorter than 20 min. The SST result of the recorded respiratory signal demonstrated in Fig. 1 is shown in Fig. 3.

A. Study Material
To validate the combination of the model and the SST algorithm in deriving the dynamics of a physiological signal with large oscillatory scale, we consider the following database collected in a recent study [15] for the purpose of studying the ventilator weaning problem.
All protocols in that study [15] were approved by the Institutional Review Board of Taipei Veterans General Hospital, Taipei, Taiwan, and written informed consent was obtained from patients. The study subjects were 68 ready-for-weaning intubated patients collected in the intensive care unit of Taipei Veterans General Hospital, Taipei, Taiwan. In particular, all subjects are with RSBI ≤110 breaths/min/L since we excluded patients not ready for weaning with RSBI >110 breaths/min/L before spontaneous breathing trial (SBT) due to the restriction from the Institutional Review Board. For each subject, we continuously recorded a 30 min flow signal at the sampling rate 100 Hz during SBT under the T-piece ventilator mode. The characteristics of these patients and the protocol to perform SBT are described in detail in [15]. These 68 subjects are divided into weaning success (n = 45) and failure (n = 23) groups, based upon their extubation outcomes. Extubation was defined to be successful if patients did not need the ventilator again for at least 48 h after extubation. Reinstitution of either noninvasive or invasive mechanical ventilation within 48 h of extubation was considered an extubation failure. In [15], the data show no difference between the success and failure groups in the mean values of six clinically used weaning predictors measured before subject inclusion, and also in the mean values of three breathing pattern parameters measured after subject inclusion; in other words, the mean values of these nine clinically used weaning predictors could not allow us to discriminate success cases from failures.

B. WIN Index
We now define the WIN index capturing the breathing pattern variability. Given the respiratory signal Resp(t), we apply SST to get the SSTIF and SSTAM. In the respiratory signal, the SSTAM can be understood as the instantaneous tidal volume. Then, we define the WIN index  shown in Fig. 5. As shown, the cut-off value of WIN index could allow us to separate a majority of patients with success weaning procedure from those with failure weaning procedure. Accordingly, this cut-off value of WIN index might be helpful for the physicians to make decision which patients merit weaning from mechanical ventilation. Notice that our method requires only three consecutive minutes of respiratory signal observation. It is important to note that the occurrence, during the observation windows, of machine calibration or some breathing irregularities, such as coughs, invalid trigger, etc., which normally require special attention in the existing analyses, is not an impediment for our method. These properties of the SST for a patient in the success group are shown in Fig. 3 for demonstration.

IV. DISCUSSION
In this paper, we introduced a phenomenological model and the synchrosqueezing transform (SST) to alleviate the limitations of traditional methods for the analysis of oscillatory physiological signals. The usefulness of this combination is shown by applying it to the study of the ventilator weaning problem. Our results show that the SSTIF and SSTAM and their derived quantity WIN from the respiratory signal provide a suitable criterion for clinical use to predict weaning outcome, with an ROC-AUC of 0.76.
At first sight, this result does not improve upon what has been reported for RRV-based predictions [14]- [16]. One should take into account, however, the following two points. First, the WIN requires only 3 min of consecutively recorded data rather than the at least 20 min typically required for the RRV-based indices; since a 3 min-observation is entirely feasible in clinical practice, whereas a 20 min-observation is much less so this is a crucial difference. Second, the definition of an oscillation and its landmark are not critical in the analysis. Note that even if the definition of an oscillation and its landmark are precisely given, the landmark detection is typically not robust under perturbations by noise, e.g., invalid breathing in the respiratory signal. Although the necessarily noisy signal can be "denoised," existing noise removal or noise reduction algorithms cannot guarantee a reliable determination of the "true" landmarks. Even if "oscillation-based segmentation" could be carried out perfectly on a long enough recorded physiological signal so that enough oscillations are collected, it is not straightforward to obtain such an uninterrupted signal, in practice, especially when the period of each oscillation is large than, say, 2 s. Different kinds of interruptions during the signal record, for example, machine calibration, coughs, suction, irritation, and so on, require the signal analyst to artificially cut-and-stitch together pieces of signal so that the oscillation-to-oscillation intervals can be defined, as needed for the traditional method, even though there is no theoretical support to show that the dynamical variabilities evaluated from such manipulated time series are still as trustworthy. Similar potential disruptions typically occur in other oscillatory physiological signals. We therefore expect that the combination of the proposed phenomenological model and the analysis method, SST, have great potential in dealing with broader physiological signals than the traditional approaches, especially when the scale of the signal is large and the determination of the landmarks is not easy.
In the introduction, we identified three drawbacks of presently used methods. We have already shown that the method proposed here does not suffer from the first two drawbacks: we can use much shorter time observation time observation windows, and we need not identify precisely the oscillating pattern to be detected. The third drawback we formulated was the relative poverty of the traditional landmark-to-landmark time interval series. The importance of keeping a richer description is illustrated by our method's use of variation in both frequency and amplitude in the ventilator weaning example.
Physiologically, the respiration is not only controlled by the neural respiratory center, but also controlled by the arterial chemoreceptors, lung vagal sensory receptors, lung mechanics, etc. In a normal subject, these control factors are integrated in a complex way which leads to the RRV [12], [14], [15]. Although we do not have a definite evidence, based on the reasoning about the relationship between the decreased HRV (or RRV) and the severity of disease [2], [4], [14]- [16], we hypothesize that the distinguishing power of WIN is a consequence of the possible disintegrity of the respiratory control factors, and this disintegrity leads to the decreased variability in the breathing pattern. This hypothesis indicates that the subject with decreased WIN is not completely ready for weaning, and explains why we observe the different patterns in the two groups. Also note that RRV is a different notion compared with the descriptive breathing parameters, like total volume, peak respiratory flow, total breath duration, and so on, which quantify the average behavior of one breath but not the complex integrity of the control factors.
This discussion would not be complete without listing the shortcomings of our approach. First, the phenomenological nature of our analysis limits the possibility to extract detailed underlying mechanisms leading to the variability and hence the prediction outcome. For example, we are not able (nor did we attempt) to distinguish if the observed variability is purely neural in origin, or if mechanical factors also play a role; finer modeling is needed to decide this, if possible. Second, from the viewpoint of estimating physiological dynamics, monitoring variability from the signal of one physiological subsystem, for example, the respiratory signal, is not sufficient: the physiological dynamics is the outcome of the complicated interactions between different physiological subsystems. Incorporating different simultaneously recorded biomedical signals, such as electrocardiograms, respiratory signal, blood pressure, and so on, should lead to a more informative description of the systematic dynamics. A study of how to extract the information on the interaction between different subsystems via the combination of the existing model and the SST algorithm is now ongoing. Specifically, although we now have a suitable method to extract information from different oscillatory physiological signals of different scales, finding their interactions and inferring more information remains an open problem. Third, the limitation inherent to the SST itself cannot be overlooked. Indeed, it is an intrinsic limitation of the SST that its estimation of the instantaneous frequency is less reliable when this instantaneous frequency has large local variations. To address this, we need a better theoretically rigorous approach to estimate the instantaneous frequency for such signals; work on this is ongoing. Finally, the study in the testbed has two limitations: it is retrospective, and the data have been collected for a small clinical population only, although the group of patients is homogeneous, in the sense that all patients are confirmed to be ready for weaning based on the RSBI. Yet it certainly warrants a follow-up prospective and large-scale clinical study to investigate its clinical applicability.
Despite these shortcomings, our example supports that the criterion we propose (in the form of the WIN) has the potential to assist physicians in assessing weaning readiness. The result encourages its application to other different kinds of oscillatory signals, in particular to those with (relatively) long periods.
In conclusion, the proposed model and algorithm together efficiently relieve the difficulty shared by the traditional methods presently used to analyze physiological dynamics or more general oscillatory signals-the data length needed for analyzing the dynamics is significantly shortened, the effect of the inevitable noise is reduced, and the pattern of the oscillatory phenomenon does not play a significant role.
Besides the traditional analytic approach, for example, the Hilbert transform, there are many existing TF analysis method available to estimate φ (t) of the signal expressed in (9) [35]. Due to the Heisenberg uncertainty principal, the ambiguity in the TF representation is inevitable, and the reallocation technique was proposed for the sake of sharpening the TF representation [24]- [26], [34]- [36]. SST is a special reallocation technique. In the following, we start from giving the precise conditions and definitions about the functional class modeling the oscillatory physiological signal, and then discuss the reallocation technique and its special case SST. We mention that we can also consider STFT, but here we focus on the CWT to simplify the discussion.
Fix a Schwartz function ψ so that supp ψ ⊂ [1 − Δ, 1 + Δ], where 0 < Δ < 1 which is commonly called the mother wavelet. Recall that the CWT [37] of a given f (t) ∈ S is defined by where a > 0 and b ∈ R. Here, we follow the convention in the wavelet literature that a means scale and b means time.
To ease the notation, the moments of ψ are denoted as I

A continuous function is called an nonharmonic intrinsic mode function if it satisfies
where and is a small parameter chosen in the model; II) s : [0, 1] → R is C 1,α , where α > 1/2, and 1-periodic function with unit L 2 norm, | s(k)| ≤ δ| s(1)| for all k = 1, where δ ≥ 0 is a small parameter, and n> D |n s(n)| ≤ θ for some small parameter θ ≥ 0 and D ∈ N. Note that we assume that the signal has just one component, unlike the case in previous studies wherein the signals were considered to include several components [22]. We shall call φ the phase function of the signal f (t) and the derivative φ (t) of the phase function the instantaneous frequency (IF) of f (t). Next, we model the measured physiological signal as where s(t), A(t), and φ(t) satisfy (I) and (II) and the following condition is satisfied: III) T : R → R is in C 1 (R) so that its Fourier transform exists in the distribution sense, and |T (ψ a,b )|, |T (ψ a,b )| ≤ C T for all b ∈ R and a ∈ (0, 1+Δ c 1 ], for some C T ≥ 0 [21]. IV) Φ(t) is a stationary generalized random process (GRP) or "almost" stationary GRP [21] independent of A(t)s(φ(t)), which is introduced to model the measurement error. To be more specific, the power spectrum dη of the given GRP Φ satisfies (1 + |ξ|) −2l dη < ∞ for some l > 0. Also assume σ ∈ C ∞ so that σ L ∞ 1 and σ := max =1,...,max{1,l} { σ ( ) L ∞ } 1, and varΦ(ψ) = 1. Now we discuss the reallocation technique and the SST. Take the CWT of a given observation Y in (9) Note that (10) is well defined since σ(t)Φ(t) is GRP and by the assumption of T . The reallocation technique "sharpens" W Y by "reallocating" the value at (t, ξ) to a different point (t , ξ ) according to some reassignment rules [25], [26], [34], [35], where t might be different from t. In contrast to the reallocation technique, in the SST, W Y is reallocated from (t, ξ) to a different point (t, ξ ) according to a different reassignment rules which only reallocates the frequency axis and preserves the time information. Note that it is important in biomedical application, in particular when prediction is the purpose, since in general we do not have the future information.
Here we detail the SST, which is separated into three steps. First, calculate W Y (a, b). Second, calculate the reallocation rule ω Y defined on R + × R: Third, the SST of Y (t) is defined by reallocating the coefficients of W Y (a, b) according to the reallocation rule ω f (a, b): where (b, ξ) ∈ R × R + , Γ > 0 is the threshold chosen by the user, h is a kernel function which is smooth enough and decays fast enough. Intuitively, at each time point t, we collect all CWT coefficients with scales a at which ω Y (a, b) is close to ξ and put them in the (b, ξ) slot. As is shown in [20]- [22], S Γ Y (b, ξ) will only have dominant values around φ (b). This property allows us an accurate estimate of φ by, for example, the curve extraction technique. The estimated φ is denoted by φ .
With the estimated φ , we can estimate A(t) and φ(t) in (9) by the following reconstruction formula. Define where R ψ := ψ (ζ ) ζ dζ and Re means taking the real part. Based on the theorem in [20]- [22], the estimator of A(t) is defined as A(t) := | R(t)|. Then an estimator for φ(t), denoted as φ(t), can be obtained by unwrapping the phase of the complexvalued signal R(t)/ A(t). As is shown in [21,Th. 3.1], these estimators are accurate and are robust to the existence of the trend T (t) and noise σ(t)Φ(t).
Notice that we can interpret (13) as an adaptive band-pass filter. Indeed, at each time t, the estimation of φ (t) by SST can be interpreted as the main frequency region corresponding to the physiological signal at time t. Then, the reconstruction of the physiological signal at time t is based on the chosen frequency band with bandwidth 1/3 .
We summarize the theoretical results of SST relevant to this paper, and refer the reader to [20], [21], [23] for the precise statement of the theorem. P 1 : SST is robust to the several different kinds of noise, which might be slightly nonstationary. Thus, we are able to accurately estimate the IF and AM [21], [23]. P 2 : Since SST is local in nature, we are able to detect components that do not exist all the time and hence the dynamical behavior of the signal [20], [21]. P 3 : The time-frequency representation S Γ Y (b, ξ) is visually informative. See Fig. 3 for example. P 4 : SST is "adaptive" to the data in the sense that the error in the estimation depends only on the first three moments of the mother wavelet instead of the profile of the mother wavelet. See [23, Fig. 6] for a visual demonstration. P 5 : The existence of the trend modeled in (9) do not interfere with the SST. Thus, we are able to separate the trend and periodic components. Note that a smooth function T ∈ C ∞ ∩ S so that its Fourier transform T is compactly supported in a small interval around 0 is a special case.

APPENDIX B NUMERICAL IMPLEMENTATION OF SST
In this section, we provide the numerical implementation detail of SST. The MATLAB code is available in http://sites. google.com/site/hautiengwu/ and we refer the readers to [23] for more implementation details.
Take a discretization of (9), denoted as Y ∈ R N , with the sampling period τ > 0 from time τ to time Nτ , that is, Y (n) = Y (nτ ). Here, we have to be careful about the meaning of Y (nτ ) when Φ is a GRP. Theoretically, when Φ is in the general sense, for example, the differentiation of the standard Brownian motion, it does not make sense to directly evaluate Φ at a particular time [38]. Also, as is discussed in [21] and the references therein, the relationship between the continuous random process and discrete time series is not always one to one. For example, not every autoregressive and moving average time series (ARMA) can be embedded into a continuous ARMA random process. Thus, in the discrete case, sometimes the random error term needs to be modeled separately, that is, where n = 1, . . . , N, f, T and σ satisfy conditions (I), (II) and (III), and Ψ is a zero-mean stationary time series, which can be taken as an ARAM time series, generalized autoregressive conditional heteroskedasticity time series, etc. In the following, to simplify the discussion, we assume that the discretization of (9) can be carried out directly, and refer the reader to [21] for more theoretical results about discretization.
In practice, to prevent boundary effects, we pad Y on both sides by reflecting the signal near the boundary so that its length is N = 2 L +1 , where L is the minimal integer such that N > N. Although it works well in practice, we emphasize that doing so is not the optimal solution. To ease the notation, in the following, we use the same notation Y to denote the padded signal and N to indicate the length of the padded signal.

A. Step 1: Numerically Implement the CWT
For the sake of implementing (10), we take the scales a j = 2 j /n v τ, j = 1, . . . , Ln v , where n v is the "voice number" chosen by the user. In practice, n v = 32 performs well. We denote the implemented CWT as a N × n a matrix W Y . We may directly follow the code available from wavelab 2 which imple-ments CWT in the Fourier domain, or directly implement it by convolution in the time domain.

B. Step 2: Numerically Implement ω Y (a, b)
To numerically implement the reallocation rule ω Y (a, b) (11), we have to discretize ∂ b W Y (a, b). It is evaluated directly by finite difference at the time axis b and we denote the result as a N × n a matrix ∂ b W Y . The ω Y (a, b) is then implemented as a N × n a matrix ω Y by the following entrywise calculation: Note that although W Y (i, j) = 0 is numerically unstable, but this unstability will be taken care in the next step.

C.
Step 3: Numerically Implement S Γ Y (b, ξ). To implement the Synchrosqueezing transform S Γ Y (b, ξ) (12), we discretize the frequency domain [ 1 N τ , 1 2τ ] by equally spaced intervals of length Δ ξ = 1 N τ . Here, 1 N τ and 1 2τ are the minimal and maximal frequencies detectable by the Fourier transform theorem. Note that the dc (direct current) term is not considered in SST. Denote n ξ = , which is the number of the discretization of the frequency axis. Fix Γ > 0, the discretized S Γ Y (b, ξ), denoted by a N × n ξ matrix S Y , is evaluated by where i = 1, . . . , N and j = 1, . . . , n ξ . Notice that the number Γ is a hard thresholding parameter, which is chosen to reduce the influence of noise and numerical error encountered when W Y (i, j) is small. We choose Γ = 1e − 3 in this study. If the error is Gaussian white noise, the choice of Γ is suggested in [23]. In general, more study is needed to adaptively choose Γ from a given time series.
To visually see the time frequency representation of the SST, we may directly plot the N × n ξ matrix R Y defined by R Y (i, j) := |S Y (i, j)| 2 for all i = 1, . . . , N and j = 1, . . . , n ξ , and observe its behavior.

D. Step 4: Estimate IF, AM, and Trend From S Y
According to the theoretical statement, the SST time frequency representation will be dominant in the IF [21], [23]. Thus, we estimate IF by fitting a discretized curve c * ∈ Z N n ξ , where Z n ξ = {1, . . . , n ξ } indexes the discretized frequency axis, to the dominant area of S Y . Precisely, we maximize the following functional over c ∈ Z N n ξ : where λ > 0 determines the regularity of the estimated curve. The first term is aimed to capture the maximal value of S Y at each time and the second term is aimed to impose regularity of the extracted curve. In this study, we simply choose λ = 1.
Denote the maximizer of the functional in (14) as c * ∈ R N . Then SSTIF, denoted as φ ∈ R N , is given by where n = 1, . . . , N, f 1 (n) ∈ C, and is the real part. The phase function φ is then estimated by unwrapping f 1 (n)/|f 1 (n)|. Note that estimating the phase by integrating the estimated IF is not recommended since the error might be accumulated. Lastly, the trend T (t) at time t = nτ can be estimated by where T ∈ R N and ξ l is chosen by the user, if recovering the trend is necessary.