Modeling the Pulse Signal by Wave-Shape Function and Analyzing by Synchrosqueezing Transform

We apply the recently developed adaptive non-harmonic model based on the wave-shape function, as well as the time-frequency analysis tool called synchrosqueezing transform (SST) to model and analyze oscillatory physiological signals. To demonstrate how the model and algorithm work, we apply them to study the pulse wave signal. By extracting features called the spectral pulse signature, and based on functional regression, we characterize the hemodynamics from the radial pulse wave signals recorded by the sphygmomanometer. Analysis results suggest the potential of the proposed signal processing approach to extract health-related hemodynamics features.


Introduction
Hemodynamics is an essential ingredient of cardiovascular physiology which not only reflects the forces the heart needs to pump blood through the cardiovascular system, but also reflects the integrity of an individual's physiological system. While there are many aspects to hemodynamics [1][2][3][4][5], it is common to evaluate hemodynamics by assessing the pulse wave signals recorded at different locations. For example, the brachial cuff systolic and diastolic blood pressures derived from the pulse wave have important clinical applications, long established beginning 1870 [6]. In addition to the blood pressure, the pulse wave signal itself contains an abundance of clinical information. For example, in subjects with dysfunctional left ventricle, different carotid pulse patterns can be observed, such as hyperkinetic pulse, pulsus alternans, pulse bisferiens, pulsus parvus et tardus, etc [7,8]. However, it is difficult to obtain these carotid pulse wave signals non-invasively [9]. On the other hand, the non-invasive pulse diagnosis (pulsology) based on, for example, the brachial or radial pulse signal, provides several aspects of physiological information, such as the central pressure wave [10,11]. Other indices associated with the pulse wave signals recorded from different locations, including the augmentation index [5, Section 6.1] [9,12], are also available for clinical applications. Since there is a great amount of information within the pulse signal, a deeper and more extensive understanding of the pulse wave is undoubtedly important to better assess not only the cardiovascular but also the physiological integrity.
The common approaches to analyzing the pulse wave signals can be classified into two categories-time domain analysis and frequency domain analysis. In both categories, it is necessary to have either a representative pulse wave cycle, which reflects the interaction between the one time heart pump driving force and the impedance to blood flow, or a sequence of pulse wave cycles over a period of time, during which the heart rate is relatively constant. In the time domain analysis, different landmarks or features are identified from the pulse signal. For example, with the radial pulse signal is characterized by the height of the percussive wave, the height of the dicrotic notch, the length of the cardiac cycle, the length of acute ejection, etc. Different features with different physiological meanings can be obtained from the pulse signals recorded from other areas, and these features have been widely used in clinics for health evaluation [8]. Furthermore, new mathematical methods have been applied to investigate additional physiological information. For example, the central pulse wave could be reconstructed from the radial pulse wave by a generalized transfer function, and then be separated into inflection and reflection waves in order to evaluate markers of early vascular aging [13]. In the frequency domain analysis, the energy of different harmonic modes estimated by power spectrum analysis are also potential features for clinical applications. According to the resonance theory [14], different harmonic modes are associated with the integrity of different organs [15]. Several existing works show the potential of this spectral approach, such as the relationship between the internal organ disorders and the spectral content [16,17]. Although this approach is commonly applied, it has two significant limitations.
The first limitation is the over-simplification of the model underlying the analysis techniques. For example, heart rate variability (HRV) and pulsus alternans are often not used in the analysis. Though the time domain landmark features reflect the hemodynamics from different aspects, they might be confounded by HRV and pulsus alternans. Since HRV and pulsus alternans are inevitable physiological dynamics, without considering them in the model and analysis, both time domain and frequency domain analysis results may be inaccurate and the results might lead to incorrect interpretation. The second limitation is that it is not always straightforward to determine a representative pulse cycle or a period of time when the pulse wave signal is suitable for these analysis techniques, so that human intervention and subjective decision making are required. Indeed, to determine a pulse cycle, it is necessary to define what a pulse cycle means or, at least, to determine landmarks within the pattern to separate one pulse cycle from another. For example, with an electrocardiogram (ECG) signal, the repeating basic pattern is related to the electrophysiological activity of a normal heart beat, and a common chosen landmark for this is the R peak. As the R peak is usually dominant and significant, determining the oscillatory pattern for an ECG signal is not difficult. However, for pulse signal such as those acquired in this study, it is not always easy to define an oscillatory pattern or a landmark, in particular when the signal is recorded from an abnormal subject. Even more seriously, due to the possibility of a noisy record, the pulse cycle morphology can dramatically change from one pulse to another or from one subject to another (See Fig 1). Although this difficulty can be mitigated to some extent by noise reduction algorithms, due to the non-linear nature of the signal, their performance is not guaranteed.
To address these limitations, in this paper we apply a recently proposed and analyzed descriptive model, called the adaptive non-harmonic model, to describe the pulse wave signals. This model is characterized by what is referred to as the wave-shape function [18], which is defined in order to capture what a pulse wave cycle should look like, including the features in the time domain and frequency domain mentioned above. The main characteristic of the proposed model is that the HRV and pulsus alternans are decoupled from the wave-shape function. A companion algorithm, referred to as the synchrosqueezing transform (SST), is applied to provide an accurate estimation of the wave-shape function. SST is a recently proposed timefrequency analysis technique that has a mathematical foundation and robustness to different kinds of noise and different practical issues have been well established in the literature [19,20]. By the adaptive non-harmonic model and SST, the above-mentioned limitations can be alleviated. By a suitable manipulation of the estimated pulse wave cycle, a vector can be obtained as an associated feature, which further leads to our final index, referred to as the global pulse signature (GPS) [18,21]. This paper is organized as follows. In Section 2, we systematically model the pulse signal using the notion of wave-shape function and the adaptive non-harmonic model. In Section 3, we summarize the SST algorithm and discuss how the wave-shape function can be estimated based on the SST and functional regression technique and then generate the GPS index. In Section 4, we apply the proposed algorithm to study the radial pulse signal. The data collection procedure, the physiological background of pulsology and the analysis result are reported. In Section 5, we discuss our findings with their clinical interpretations and indicate directions for future research.

Adaptive non-harmonic model for the pulse signal
In this section, we propose a mathematical model to describe the pulse wave signal. It is well known that hemodynamics is a highly complex subject [2-4, 22, 23]. Different signals recorded from different areas of the body by different instruments can shed light on different aspects of hemodynamics. The question of modeling hemodynamics has been discussed by several authors, who have proposed different models based on different physiological factors, such as An example of the radial pulse signal recorded from a subject with congestive heart failure. Note that at the 3 rd second there is an obvious artifact, indicated by the blue arrow, and at the 5.5 th second there are two shock noise events. The relatively small pulse wave at the 3 rd second is related to a premature atrial contraction. This is significant since patients with heart failure may have more premature beats than a normal control.
pressure gradient, resonance, vibration, etc, have been proposed; see [1,3,4,24] and the literatures cited therein for detailed discussions of these models. While these models shade a light on the understanding of the pulse wave signal, having a signal model to effectively quantify the pulse behavior remains a distant goal. In particular, inevitable physiological aspects, such as HRV and pulsus alternans, are commonly missed or disregarded in the existing models. Ignoring them in the model and analysis, however, might not only limit the analysis but also lead to incorrect findings. In this paper, we apply a different model and approach to resolve these issues.
The pulse wave signal is oscillatory in nature, and there are several well-known features that we consider here: how fast the signal oscillates, how much the oscillatory rate changes, how differently the signal behaves inside an oscillation, etc. Understanding these features plays a significant role in the data analysis, which can help to obtain more information about the physiological system. Based on these observations, we consider the following adaptive non-harmonic model [18] to model the pulse signal. This model differs from most existing models in that it is purely phenomenological; that is, the parameters in the model are solely driven by observations of the physiological signal, rather than by explicit and quantitative underlying mechanisms.
2.1.1 Wave shape function. Based on the periodic nature of an oscillatory function, we introduce the wave shape function, which is used to model how a signal oscillates over the period. This idea has been previously studied in [18,25,26]. To resolve the difficulty of defining a period, we shall call a function f τ-periodic, where τ > 0, if (1) the periodic condition is satisfied; that is, for all t 2 R and all k 2 Z, f(t + kτ) = f(t); and (2) for all 0 < τ 0 < τ the periodic condition is not satisfied. Take two positive values δ and θ and a positive integer D. A 1-periodic function s is called a wave shape function with parameters δ, D and θ if the following conditions are satisfied. First, the function is differentiable and its derivative is a Holder continuous function with the Holder coefficient α > 1/2. Denote b s to be the Fourier transform of the function s. The second condition that a wave shape function should satisfy is that s has the unit energy (unit L 2 -norm) and all the Fourier modes b sðkÞ, k 6 ¼ 1 are dominated by the product of δ and the first mode coefficient; that is, Furthermore, it is desirable for b s to be mostly concentrated in the low frequency region, which is quantified by θ by: X n>D jnb sðnÞj y: Conditions (1) and (2) reflect the practical finding of the spectral analysis of the pulse signals [27][28][29], and it can be observed that the amplitudes of tenth and higher order harmonics are negligible. Note that the commonly selected landmarks and lengths considered in the time domain analysis could be understood as morphological features describing the wave-shape function that is modeling how a pulse repeats itself. While there is a one-to-one relationship between these landmarks and Fourier coefficients b s, this relationship is nonlinear in nature. Thus, while we could view the Fourier coefficients of the wave shape function as the "features" of the pulse wave, their physiological interpretations are not directly related to the hemodynamic interpretations of these landmarks. See Fig 2 for an example of the wave-shape function describing the pulse wave signal. 2.1.2 Adaptive non-harmonic model. Using the notion of wave shape function, we now describe our phenomenological model to capture the recorded pulse signal. Fix parameters δ, θ, D for a wave shape function s. We consider the following Intrinsic Mode Type (IMT) functions to model the pulse wave signal. An IMF function, f, is a bounded and continuous function with a continuous derivative that satisfies the following format: where A is a positive differentiable function and ϕ is a monotonically differentiable function. Intuitively, A(t) describes how large the oscillation is at time t, and the positive function ϕ 0 (t), the first derivative of ϕ, describes how fast the oscillation is at time t. To see how it is interpreted in this way, consider a constant positive function A and a linear function with a positive slope as ϕ. Suppose ϕ(t) = ξ 0 t, where ξ 0 > 0. In this case, we know that f is a harmonic function with the frequency ξ 0 and the amplitude A. We could thus view Eq (3) as a generalization of the harmonic function. Though the heart rate is not constant, normally it does not change suddenly. Therefore we need following conditions to better quantify the pulse signal. Fix a small positive constant . Then we assume that A(t) does not change too fast in the sense that its derivative is bounded by ϕ 0 (t) and that ϕ 0 (t) does not change too fast in the sense that ϕ@(t) exists and is bounded by ϕ 0 (t); that is, we have jA 0 ðtÞj 0 ðtÞ and j 00 ðtÞj 0 ðtÞ for all time t: We would thus model the pulse wave signal by Eq (3) with the Condition (4), and call this model the adaptive non-harmonic model. We would call the monotonically increasing function ϕ(t) the phase function, ϕ 0 (t) the instantaneous frequency (IF), and A(t) the amplitude modulation (AM). An important issue regarding the identifiability issue of the phase function, the IF and the AM cannot be ignored if the discussion is to be rigorous. Indeed, there might be infinitely many different ways to represent a cosine function g 0 (t) = cos(2πt) in the format a(t) cos(2πϕ(t)) so that a > 0 and ϕ 0 > 0, even though it is well known that g 0 (t) is a harmonic function with amplitude 1 and frequency 1. Indeed, we could find infinitely many smooth functions α and β so that g 0 (t) = cos(t) = (1 + α(t)) cos(2π(t + β(t))), and in general there is no reason to favor α(t) = β(t) = 0. Clearly, before resolving this issue, amplitude 1 and frequency 1 cannot be taken as reliable features to quantify the signal g 0 . This identifiability issue has been well studied in [20,30], finding that if g 0 (t) = A(t)s(ϕ(t)) and g 0 (t) = (A(t) + α(t))s(ϕ(t) + β(t)) also holds, and both satisfy the Condition (4), then |α(t)| and |β 0 (t)| are both bounded by for all time t 2 R. Note that the IF and AM are always positive, but usually not constant. Clearly, when ϕ is a linear function with a positive slope and a is a positive constant, then the model is reduced to the harmonic model and the IF is equivalent to the notion frequency in the ordinary Fourier transform sense. The conditions |a 0 (t)| ϕ 0 (t) and |ϕ 00 (t)| ϕ 0 (t) force the signal to locally oscillate "regularly", that is, around time t 0 , the signal A(t)s(ϕ(t)) oscillates like A(t 0 )s(ϕ(t 0 ) − t 0 ϕ 0 (t 0 ) + ϕ 0 (t 0 )t), and hence the nominations of IF and AM. We mention that this model is a special case of a wider class of model composed of multiple oscillatory components in the sense that we only have one oscillatory component in the model. For details about a more general model, see [18][19][20]30].

Physiological interpretation.
The main reason we consider a time-varying frequency and amplitude model like Eq (3) is to capture the physiological facts of HRV and pulsus alternans.
Since the heart does not beat at a constant rate [31], the pulse signal should not oscillate with a constant frequency. The non-constant heart rate is modeled by IF and hence the HRV could be evaluated from analyzing IF. In the past decade, due to the health information embedding in HRV and the trend towards personalized health care requirements, estimating the HRV from the pulse signal extracted from different resources has attracted significant research [32][33][34]. However, the existence of HRV is commonly ignored in the pulse analysis literature. For example, in the traditional spectral analysis approach to the radial pulse, researchers need to analyze the pulse signal over the interval where the heart rate closely resembles a constant [15,29].
The possible pulsus alternans is captured by the AM. This model is particularly important in patients susceptible to heart failure since pulsus alternans is a manifestation of severe impairment of the left ventricular systolic function [7]. An example of a wave-shape function and an IMT function representing the pulse signal are shown in Fig 2. 2.1.4 Recorded pulse signal. In practice, noise is inevitable when recording signals. Thus, we model the recorded pulse signal by where A(t)s(ϕ(t)) is the adaptive non-hamornic model for the pulse signal, σ is a slowly varying smooth function quantifying the possible non-stationarity in the data collection process, and F is a stationary random process with unit standard deviation describing the noise. See, for example, [20] for a discussion of the noise. Here, we do not assume that the noise is Gaussian or white, while the commonly encountered Gaussian white noise is a specific example. Since in general the noise might be complicated with time-dependence and be far from Gaussian, we take this noise model into account. While the pulse signal could be recorded by different ways, like tonometer, photoplethymography, video [33], etc, the wave-shape function could be different for different types of pulse signals.

Synchrosqueezing transform
It is a common practice to apply the Fourier transform to study oscillatory signals, in particular the pulse wave signal in which we are interested [15][16][17]. As useful as the spectral analysis is, however, it is well known that when the signal is not composed of harmonic functions, the power spectrum determined by the Fourier transform might be misleading. For spectral pulse wave analysis, since HRV and pulsus alternans are inevitable, in practice analysts must carefully choose the signal. But there is still no guarantee that the HRV and pulsus alternans influence could be eliminated, and this could become a confounder in the analysis.
While the oscillatory signals with time-varying AM and IF are ubiquitous, much effort has been extended over the past few decades to address this problem, and several approaches have been proposed, ranging from empirical mode decomposition and its variations [35][36][37], to optimized window and the approximation theory approach [38,39], to name but a few. Moreover, research interest has recently been extended to multivariate time series analysis [40][41][42][43][44] to further take the spatial information into account. Among the different approaches, one active field regarding this issue is time-frequency (TF) analysis. Based on different principles, several TF analysis techniques have been proposed. Typical examples include linear methods such as the short time Fourier transform (STFT) and the continuous wavelet transform (CWT), or the quadratic methods such as the Wigner-Ville distribution or the Cohen class. We refer the reader to a standard textbook [45] for more information. While these methods have attracted a great deal of attention in different fields other than signal processing, they are limited by either the Heisenberg uncertainty principal or the mixing issue. We demonstrate this idea by discussing STFT. The main idea forming the basis for STFT is dividing the signal into overlapping small pieces, and then studying the spectral behaviors over these small pieces. Mathematically, for a function f, its STFT associated with a window function h(t) can be defined as where t 2 R is the time, Z 2 R þ is the frequency, h is the window function determined by the user-for which a common choice is the Gaussian function with kernel bandwidth σ > 0, i.e. h(t) = (2πσ) −1/2 e −t 2 /σ2 . However, the Heisenberg uncertainty principal limits how well the spectrum could be estimated over each piece, thereby limiting the STFT. Similar discussions hold for CWT and other linear TF analysis techniques, and we refer the reader to [45] for details. In the current work, in order to capture the hemodynamics, which have oscillations on the order of 1 second, we run STFT with σ = 0.5 so that the window is not too long to lose the dynamic information, and is not too short to cause a numerical artifact.
To handle these fundamental issues, a nonlinear TF analysis technique, synchrosqueezing transform (SST) [19,20], which is a special reassignment technique (RM) [46], is proposed to obtain a sharper time-frequency representation. Since these techniques were introduced, they have been widely applied in different fields; see [47] for a summary of its applications in different fields.
Here we briefly summarize SST, and refer the reader to [46] for RM. We mention that the SST could be applied to different linear TF analysis, like STFT, CWT, wave packet transform [25] or S-transform [48] and theoretically the results do not depend significantly on the chosen linear TF analysis. But to simplify the discussion, we choose to do the analysis with STFT. The Matlab code for the SST algorithm based on both CWT and STFT can be downloaded from https://sites.google.com/site/hautiengwu/home/download. In brief, the SST sharpens the TF representation determined by STFT by reallocating the STFT coefficients along the frequency axis to the "correct" frequency slot which represents the IF of an oscillatory component. Mathematically, The SST of f is defined as where g α is an approximate δ-function in the sense that g is a fast decaying smooth function with R g(x)dx = 1, so that g a ðtÞ :¼ In plain language, by reading o ðhÞ f ðt; ZÞ, we collect all STFT coefficients indicating the existence of an oscillatory component with frequency ξ to the slot ξ. Note that compared with RM, in SST, the coefficients are reallocated along the frequency axis, so the causality is preserved; second, in SST we reallocate the STFT coefficient instead of the spectrogram coefficient. These two facts allow us to reconstruct the oscillatory components of interest, in particular when the signal is noisy or is composed of several oscillatory components. It has been clearly found that for the phenomenological model of the pulse wave signal of interest, the IF ϕ 0 (t) and the AM A(t) can be accurately estimated from the recorded pulse signal [19,20]. Precisely, we could prove that at time t, the coefficients in S ðhÞ f ðt; xÞ are dominant when ξ % ϕ 0 (t). This property allows us an accurate estimate of the IF ϕ 0 by, for example, the curve extraction technique.
Denote the estimated ϕ 0 by e 0 . We can then estimate the amplitude modulation A(t) and the phase function ϕ(t) by building e RðtÞ :¼ hð0Þ The estimator of A(t) is thus defined as e AðtÞ :¼ j e RðtÞj, and hence an estimator for ϕ(t), denoted as e ðtÞ, can be obtained by unwrapping the phase of the complex-valued signal e RðtÞ= e AðtÞ. We refer readers interested in SST to [19,20] for the detailed numerical algorithms and the theory beyond them.
Here, we show the results of pulse analysis by applying the SST in Fig 3. To demonstrate its benefit on a noisy signal, we artificially add noise on the second half of the signal, and show the result in Fig 4. In this example, the noise is an ARMA(1,1), where ARMA stands for autoregressive and moving averaging, time series determined by the autoregression polynomial a(z) = 0.5z + 1 and the moving averaging polynomial b(z) = −0.3z + 1, with the innovation process taken as i.i.d. student t 3 random variables so that the signal to noise ratio, defined as 20 log stdðsignalÞ stdðnoiseÞ , is 0 dB. It is clear that though the signal to noise is low, the noise is non-stationary since it exists only over a finite interval, and the noise has the fat tail behavior described by the student t 3 random variables. Thus, the SST algorithm could reliably extract the instantaneous frequency, so that we could obtain reliable HRV information.

Feature extraction by estimation of the wave shape functions
As discussed in Section 2.1, the main feature of interest for the pulse analysis is the wave-shape function in the adaptive non-harmonic model, in particular in the pulse spectral analysis. In this section, we describe an algorithm to extract the wave shape function. Note that since the waveshape function can be expanded by its Fourier coefficients, the pulse signal f(t) = A(t)s(ϕ(t)) can be represented as where a ' 2 R and b ' 2 R are the Fourier coefficients of the shape function s. In this study, we estimate the wave shape functions based on Eq (10) using the standard functional regression [21,49]. Consider the wave shape functions s with parameters δ, D, θ in Eq (5). To simplify our discussion, we assume that θ = 0 and that the noise is stationary, that is σ = 1. We would choose D = 6 in this study, based on the discussion in Section 2. Thus, the pulse signal Eq (10) becomes After discretization by the sampling period Δt > 0 over time interval [Δt, NΔt], the recorded pulse signal is saved digitally as a N-dim vector, Y 2 R N , so that its l-th entry is Y l = f(lΔt) + F l , where l = 1, . . ., N and F is a random vector satisfying var(F l ) = 1 for all l, which might not be Gaussian and the covariant matrix might not be the identity. Denote the discretized estimators of A(t) and ϕ(t) by e A 2 R N and e 2 R N . Note that it has been well established the discretized estimation of A and ϕ are accurate with error of order [20]. To simplify the discussion, we assume that the estimates e A and e are precise without error; that is, e AðlÞ ¼ AðlDtÞ and e ðlÞ ¼ ðlDtÞ for all l = 1, . . ., N. In the general case, the analysis result will deviate by an error of order of . We thus construct the following "functional vectors" where ℓ = 0, . . ., D and c ℓ , d ℓ are N-dim vector whose k-th entries are  (a) The radial pulse signal R(t) recorded from a normal subject that is contaminated by the autoregressive and moving averaged noise generated from the student t 3 i.i.d. random process from the 2.5-th second to the 5.5-th second, which is denoted as Y(t). The radial pulse signal is the same as that shown in Fig 3. Note that there are several spikes in Y(t), which are generated by the fat tail natural of the student t 3 random variable. (b) The time-frequency (TF) representation of Y(t) is determined by the synchrosqueezing transform (SST). The dominant curve in the TF representation is associated with the IF induced by the heart rate variability (HRV). It is clear that the TF representation determined by SST is robust to the noise even if the noise is time-dependent with a fat tail distribution, while the dominant curve representing the instantaneous frequency is slightly deformed due to the noise. where k = 1, . . ., N. As a result, the recorded pulse signal satisfies To estimate the parameters α ℓ and β ℓ from the functional vectors c, observe that the (2D + 1) × (2D + 1) matrix cc T is diagonal dominant. Since EðDtFc T Þ ¼ 0 and var(ΔtFc T ) = O(LΔt). Thus we can estimate γ by where e g ¼ ½e a 0 ; e a 1 ; . . . ; e a D ; e b 1 ; . . . ; e b D T 2 R 2Dþ1 , and hence estimate the oscillatory signal by e f ðtÞ :¼ e g T c: ð16Þ When e f contains an accurate estimation of the wave-shape function s with L 2 error of order , we could take the Fourier coefficients e g into account as the feature of the pulse signal. We call the 2D + 1 dimensional vector e g the spectral pulse signature (SPS) for the recorded pulse signal.
Note that the ℓ-th component of the power spectrum of s could be estimated by ðe a 2 ' þ e b 2 ' Þ=4. The main benefit of this approach is that the influence of HRV, which is modeled by IF, is eliminated automatically and the wave-shape function is better estimated.

Comparing the present method to previous ones
Here we summarize the difference between our approach and the commonly applied spectral analysis [27][28][29]. Recall that in spectral analysis, the power spectrum of a selected recorded pulse signal is evaluated and the energy of different harmonic modes are considered as features of the subject, indicating aspects of physiological health. The selection criteria for the pulse signal interval is that the heart rate is almost a constant [15,29], since the power spectrum approach is sensitive to the non-constant frequency and non-constant amplitude signal. The first difference is that based on the adaptive non-harmonic model and SST, the instantaneous frequency and amplitude modulation can be accurately obtained, so we do not need to select an interval from a recorded pulse signal. Second, since the HRV and pulsus alternans are physiologically inevitable, the proposed approach is more physiologically feasible.
Third, note that the information obtained from the spectral analysis is different from the SPS. Indeed, under the assumption for Eq (11), we have the following direct expansion: where a 0 ¼ŝð0Þ and for ℓ = 1, . . ., D,ŝð'Þ ¼ a ' e iy ' , a ℓ ! 0, θ ℓ 2 [0, 2π), α ℓ = 2a ℓ cos(θ ℓ ) and β ℓ = −2a ℓ sinθ ℓ . Here, θ ℓ reflects the phase of the ℓ-th harmonic component hidden inside the wave shape signal. Clearly, since power spectrum analysis commonly takes jŝð'Þj 2 ¼ a 2 ' , the phase information of the ℓ-th harmonic component is missed. However, in SPS the phase information is preserved and the hemodynamics could be more faithfully captured.

Testbed: the radial pulse diagnosis
To test how well the proposed adaptive non-harmonic model and the algorithm work in practice, in this section we study the radial pulse wave signal on congestive heart failure subjects. The signals are recorded from three landmarks. First, the radial styloid process; second the middle position between the radial styloid process and the palm and the proximal point with the same distance from the first to the second location. We mention that these three locations are commonly recognized in the pulse diagnosis; the first one is called guan, the second one is called chun and the third one is called chi. We thus use guan, chun and chi to refer to these three locations in this study.

Material
All protocols in this study were approved by the Institutional Review Board of Chang Gung Memorial Hospital, Linkou, Taoyuan, Taiwan (93-6288), Taiwan and written informed consent was obtained from all patients. Nineteen normal subjects without history of heart disease are included in the control group, and 17 subjects with congestive heart failure (CHF) are included in the study group. The diagnosis of CHF subjects is based on the criteria indicated by the Framingham heart study. The participants were invited for pulse examination in a quiet and temperature-controlled room in Chang Gung Memorial Hospital, Linkou branch, Taiwan. Pulse wave signals were recorded from chun, guan and chi positions of both hands by a tonometer (Wang's sphygmometer, PDS-2000). The sampling rate of the signal is at 100Hz. For each subject, we collect 10 seconds signal for each position on both hands, and repeat 2 or 3 times. The pulse wave was recorded in sitting position with the wrist comfortably resting on a small pillow at the level of the heart.
We recruited 17 patients with CHF for the study group and 19 normal individuals for the controls. The age of the study and control group are 64.3 ± 23.7 and 63.2 ± 15.8 respectively. The male/female ratio were 15/7 in the study group and 10/10 in the control group. There were no significant differences in age or sex. As a result, we obtain 53 (respectively 39) pulse wave signals recorded from chun from the normal (respectively CHF) group, 55 (respectively 34) pulse wave signals recorded from guan from the normal (respectively CHF) group and 50 (respectively 41) pulse wave signals recorded from chi from the normal (respectively CHF) group.

Statistical Analysis and Global pulse signature (GPS)
To test if the SPS indices of the normal group differ from those of the CHF group, we apply the currently developed one-way ANOVA for functional data, called the globalized pointwise F (GPF) test [50]. For readers having interest in this technique, we refer them to [50]. In this study, we consider p values less than 0.01 as significant.
As the SPS index is a high dimensional vector, it is not easy to visualize. To provide an easyto-use index for the pulse diagnosis, we consider the following approach to integrate the information in the SPS index. Suppose from a fixed position of the i-th subject we obtain a SPS index b g i . The associated outcome of this SPS index is denoted as y i , and y i = 1 (respectively y i = 0) means the subject is in the CHF group (respectively control group). Thus we have the dataset fb g i ; y i g N i¼1 .
When the SPS index is located in the 2D + 1 Euclidean space and the sampling size is limited, we apply the partial least squares (PLS) regression to find a linear regression model by projecting SPS and the response variables to a new space. Here we briefly recall the PLS regression. PLS regression finds components from X ¼ fb g i g N i¼1 that are relevant for the outcome Y ¼ fy i g N i¼1 by seeking a set of components that performs a simultaneous decomposition X and Y with the constraint that these components explain as much as possible of the covariance between X and Y. Then the decomposition of X is applied to predict the group. For details about PLS regression, see [51]. Suppose the PLS regression coefficient is b 2 R ð2Dþ1ÞÂ1 . Then the prediction result under the linear regression model for b g i , denoted as b y i :¼ ½1 b g i b 2 R, is referred to as the global pulse signature (GPS) index. The GPS index is integrated information derived from the SPS index, which reflects the subject's condition of interest.
For the purpose of prediction based on GPS, we can further apply the receiver operating characteristic (ROC) to determine the threshold to classify the subjects into two groups. We report the sensitivity, specificity, accuracy and the area under curve (AUC) to evaluate the classification result. The confidence interval (CI) of AUC is evaluated by 1000 bootstrap replicas. To assess how the results based on PLS and ROC will generalize to an independent data set, we run leave-one-out cross validation (LOOCV) 200 times for each position of both hands, and report the accuracy.

Results
First, we show the synchrosqueezing transform of the pulse wave signal from a subject with CHF and the associated estimated wave shape function in Fig 5. Note that due to the inevitable deviation, for example the one at the 4-th second, and the HRV, the power spectrum estimated from the pulse wave signal is spreading. Note that estimating the wave shape function could be viewed as the power spectrum analysis of the pulse wave signal after correcting the IF and AM. In other words, as we could estimate IF and AM accurately from the pulse wave signal, we could resample the signal according to the estimated IF and then normalize the signal by the estimated AM. Fig 5. (a) The pulse signal recorded from a subject with congestive heart failure (CHF), with the y-axis as the arbitrary unit. (b) The estimated wave shape function for the subject with CHF. (c) The time-frequency (TF) representation determined by the synchrosqueezing transform. It is clear that the pulse around the 2.5 th second takes a longer time to finish, which leads to the slower instantaneous frequency (the dominant curve in the TF representation is around 0.9 Hz at time2.5 th second). It is also clear that the instantaneous frequency is not constant due to the inevitable heart rate variability (HRV). Note that while there is a significant deviation at the 2.5-th second, the estimation result catches most of the shape information. Then we apply PLS to obtain the GPS index to distinguish the two groups. The histogram of GPS indices determined from different positions on both hands are shown in Fig 7. It is clear that the GPS of subjects in the CHF group is smaller than that of the control group. The ROC analysis results, including the sensitivity and specificity and AUC, from different positions are shown in Fig 8. The sensitivity, specificity, accuracy, AUC and the accuracy of LOOCV of the GPS determined from different positions on both hands are summarized in Table 1.

Discussion
In summary, in this paper we study the physiological signal by the adaptive non-harmonic model, the SST and the functional regression technique. The usefulness of the proposed scheme is supported by an encouraging analysis result of the radial pulse signal. Although we analyze the radial pulse signal as the test case for this study, it should be noted that the proposed model and analysis technique could be applied to other pulse signals obtained by different instruments. For example, contact photoplethysmogram (PPG) measurement [52] or PhysioCam non-contact PPG measurement [33], which represent the changes of blood volume in the vessel obtained through an optical transmission measurement or real-time camera images, and hence reflect the change in vascular blood volume associated with the cardiac beat. While these signals represent different aspects of the hemodynamics than the radial pulse signal we study in this paper, we could expect to obtain a broader angle of view about human health if information obtained from these signals could be combined. We will report the research progress in the future work. Since the potential of the SST and other nonlinear TF analysis techniques have been shown in this study and other clinical problems, for example,    [53][54][55], by taking the wave-shape model into account, we could expect a broader application and better analysis results. We also point out the relationship of the pulse signal analysis results with the traditional Chinese medicine (TCM) theory. In short, the results we obtained in this study could lead to a potential means to help establish the foundation of TCM theory. See Fig 9 for an illustration of the summary of the pulse diagnosis theory in the TCM theory. Note that while TCM has been widely applied in the eastern culture, up to now, a systematical/scientific theory understanding the mechanism beyond pulse diagnosis is not yet well established [56], but its usefulness has been demonstrated [58]. There are several studies based on hemodynamics aiming to understand the mechanism; for example, Wang et. al. [24,59] proposed that the pulse wave consists of numerous harmonic waves, and each harmonic wave is associated with an internal organ and carries information of different meridians over the body. That study also noted that the harmonic waves of the upper, middle and lower section of the body may correspond to the three sections of pulse in timeline and that the waveform on both hands are not totally the same. On the other hand, our approach is purely from the phenomenological viewpoint and adaptive harmonic analysis. In our results, the classification result of the CHF group and the normal group is better on the right hand side, especially the right chi. Since the right chi position demonstrates the most significant difference between groups, it may reveal the importance of the right chi (kidney yang) signal in evaluating the clinical condition of CHF. According to TCM theory, kidney yang is the fundamental support to maintain the human life. In Western medicine, the hemodynamics of renal artery, such as the blood pressure, plays an important role in regulating the overall whole cardiovascular function. On the other hand, a possible cause for this observation is that the pulse pressure of the right arm is usually higher than the An illustration of anatomical locations associated with the terminologies used in pulse diagnosis. In the TCM pulse diagnosis, it is stated that the pulse on the right hand manifests the condition of Qi [56,57], while the pulse on the left hand reflects the blood. For different pulse positions, left chun, guan and chi reflect the heart, the liver, and the kidney respectively; the right chun, guan and chi reflect the lung, the spleen, the kidney respectively. The left chi manifests mainly the kidney yin, and the right chi manifests mainly the kidney yang (the Life Gate in TCM).
doi:10.1371/journal.pone.0157135.g009 pulse pressure of the left arm. It is presumed that the pulse wave signal of the right arm may have a higher signal-to-noise ratio than the pulse wave signal of the left arm. The left chun and right guan position are the two second important roles in our analysis. According to TCM theory, the left chun position (heart) manifests the general condition of the cardiovascular system, hence our results are consistent with clinical experience. Thus, our findings partially support the pulse diagnosis theory in TCM that the waveforms on different positions of radial artery contain different information. A similar finding was reported by Young et. al., who found that the although the augmentation indices determined from the radial pulse waves recorded from chun, guan and chi are not significantly different, the estimated aortic augmentation indices determined from the radial pulse waves recorded from chun, guan and chi are not identical [60]. In conclusion, since the study of pulse diagnosis from different aspects is an active field, we would expect our proposed model and method could help to further study the experiences and working practice of pulse diagnosis, e.g. the nature and dynamic of disease, and its relationship to modern hemodynamics.
Limitations of this study should also be mentioned. First, the tonometer (PDS-2000) we applied in this study records only the two-dimensional data (pressure-time) of the pulse. Since more advanced instruments now available can obtain three dimensional data [61], further study should be undertaken. Second, although the phenomenological model we propose is capable of capturing IF and AM as well as the wave shape function, it is clearly not the optimal solution. It is clear that there is a more complex interaction between IF, AM and the wave shape function than what we consider in the adaptive non-harmonic model. On the one hand, a more general model, like the time-varying wave shape function, could be considered based on the physiological background. On the other hand, we conjecture that this relationship might be better captured by combining the existing hemodynamic models with the proposed phenomenological model. This finer model might better capture the physiological information hidden inside the pulse wave signal and lead to a better algorithm to better study the recorded pulse wave signal. A systematic study and its application of this issue will be reported in future work. Moreover, from the clinical viewpoint, the sample size in this study is limited, and the interesting clinical problems, like outcome prediction or early CHF diagnosis are not discussed. Also, to simplify the discussion and avoid possible confounders, in this study we limit our analysis to subjects with CHF. Thus, we could not make the final conclusion about the pulse diagnosis. To use the research results in clinics, a larger scale clinical study with CHF and other diseases is needed to conclude the current findings, and we will report our continuing research in the future.
Supporting Information S1 Code. The Matlab code. The Matlab code used to analyze the pulse wave signal. (ZIP)