# A Linear Transformation Approach for  Estimating Pulse Arrival Time

Authors: Dohyun Kim; Jong-Hoon Ahn; Jongshill Lee; Hoon Ki Park; In Young Kim
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

We propose a new mathematical framework for estimating pulse arrival time (PAT). Existing methods of estimating PAT rely on local characteristic points or global parametric models: local characteristic point methods detect points such as foot points, max points, or max slope points, while global parametric methods fit a parametric form to the anacrotic phase of pulse signals. Each approach has its strengths and weaknesses; we take advantage of the favorable properties of both approaches in our method. To be more precise, we transform continuous pulse signals into scalar timing codes through three consecutive transformations, the last of which is a linear transformation. By training the linear transformation method on a subset of data, the proposed method yields results that are robust to noise. We apply this method to real photoplethysmography (PPG) signals and analyze the agreement between our results and those obtained using a conventional approach.


---

## Body

## 1. Introduction

The importance of arterial stiffness as a cardiovascular disease index has been emphasized in recent years [1–6], because arterial stiffness can be acquired using inexpensive and noninvasive methods such as pulse wave velocity (PWV) [7, 8]. PWV is considered to be a good indicator for assessing arterial stiffness because it shows a strong correlation with cardiovascular events and mortality [1, 9–15]. Furthermore, PWV can be used for the continuous assessment of cardiovascular homeostasis and regulation [16].One approach to assess PWVin vivo relies on tracking pressure pulses that arise from the onset of left ventricular ejection. This is the common method for acquiring PWV in the arterial trees which uses ECG and two pressure pulses that are measured simultaneously. In general, pressure pulses are measured at the carotid and femoral arteries, respectively, and PWV is calculated as the distance between the two sites divided by the time for the pulse wave to travel that distance. The time that it takes the pulse pressure to travel from the carotid artery to the femoral artery is called pulse arrival time (PAT) [17].To measure pressure pulses at the carotid and femoral arteries, a catheter is generally used. However, it is difficult to measure pressure pulses without clinical assistance because this is an invasive method. For this reason, intensive efforts have been made to improve the performance of external skin transducers that can measure PWV in recent years. Several techniques have been developed to record pressure pulses. Among these, photoplethysmography (PPG) is particularly popular as a noninvasive, nonobstructive technique that is based on the temporal patterns of light absorption in living tissues because morphological characteristics of PPG are similar to pressure pulse, especially in the arteries [18].Pulse arrival time of PPG pulse is typically measured by detecting local characteristic points: the foot determined by the start point of the anacrotic phase (FOOT), the maximum-slope of the anacrotic phase (MS), and the maximum amplitude of the pulse (MAX) [15, 19, 20]. Unfortunately, however, these characteristic points often yield unstable and unreliable results when used to analyze PPG pulses with morphological variation due to underlying conditions [15]. Thus, the design of robust extraction techniques that are capable of estimating PAT from PPG pulses remains an unsolved problem. Solà and colleagues suggested that PAT could be estimated by parametric modeling of the anacrotic phase of pressure pulses in PPG. However, although their method produces robust and reliable results under noisy conditions, it is relatively computationally complex because of the need to fit a parametric function to every single pulse [21].Therefore, our aim was to develop a method to measure PAT with accuracy and reliability using simple operations. In the next section, we outline the mathematical framework that we developed to estimate PAT.


## 2. Methods

### 2.1. Representation of PPG Signals in Vector Space

LetΩ be the set of all continuous PPG signals measured from human arteries. We can define a sampling process ξM as follows:(2.1)ξM:Ω⟶RM,
where  ℝM is an M-dimensional Euclidean space. The mapping ξM reduces a continuous signal to an M-dimensional vector point. The vector point forms a lower-dimensional cluster in an M-dimensional space. Let us consider the cluster as being embedded by the manifold  Υ. Our goal in this paper is to find a mapping τ between this manifold and PAT:(2.2)τ:Υ⟶R1.
Then, the parametric estimation of PAT is given by the composition of two mappings: PAT=τ∘ξM.However, manifolds constructed through normal sampling at a constant frequency are highly curved. Training this type of manifold and mapping PAT using this manifold are challenging [22]. Let us consider simple translations of Gaussian peaks, as shown in Figure 1. The manifold constructed using simple translations is spirally curved. If slight variation is added to Gaussian peaks, it is not feasible to parameterize the manifold with well-defined functions.Figure 1
The estimation of PAT using a linear transformation approach.Now, suppose that we find the sampling processξM such that the manifold Υ can be flat and isometric along PAT. For instance, if three vectors f,g,h∈Υ are collinear and have isometric timing codestf, tg, th, that is,  (th-tg)f+(tf-th)g+(tg-tf)h=0, then we can always find the linear transformation ω such that tf=ωTf, tg=ωTg, and th=ωTh. This means that the special sampling process allows the mapping τ to be the simplest form by ω. However, we failed to find such a sampling process, regardless of the sampling frequencies applied to the continuous signals; convex combination of two different vectors f,g∈Υ cannot be used to represent the human artery PPG signal.In this context, we propose adding another transformation betweenξM and τ. By considering the new mapping, we intend to keep the mapping, τ, the linear transformation. If we denote the novel mapping as ζ, we can estimate the PAT of the PPG signals as(2.3)PAT=τ∘ζ∘ξM.We refer to this as a linear transformation approach for estimating PAT. In following subsections, we describe the new transformationζ in more detail.


### 2.2. Conjugate Transformations

In the previous subsection, we framed a set of three transformations to change continuous pulse functions into scalar timing codes. The first transformation,ξM, is needed to reduce continuous pulse functions into M-dimensional vector points, while the third transformation τ is the linear transformation. The function of the second transformation,  ζ, will be fully discussed in this subsection, in which we consider a well-known transformation called the convex conjugate or Legendre transformation.


#### 2.2.1. Convex Conjugate

The convex conjugate is a transformation that maps a convex function onto another convex function [23, 24]. A convex function always has its conjugate function: the conjugate function is also a convex function. First, we outline why we need convex functions. A Gaussian function, which was exemplified as a pulse in the previous subsection, is the starting point for developing our idea. Gaussian functions have a single peak and are nonnegative over the entire region. Such a Gaussian function can be derived from a convex function by differentiating the convex function twice. We can therefore consider convex functions instead of Gaussian functions. A general type of pulse function that has mixed-signed values, unlike Gaussian functions, will be discussed later.Let us consider two arbitrary convex functions. When their first derivatives become inverses of each other, two functions are referred to as “convex conjugate”. If two functionsf(t) and f̃(t̃) have such a relation, then(2.4)dfdt∘df̃dt̃=I=df̃dt̃∘dfdt,
where I is an identity function, that is, I∘t=t and I∘t̃=t̃. To find an explicit expression for f̃, we assume(2.5)t̃=dfdt∘t.
From (2.4) and (2.5), we obtain(2.6)I∘t=df̃dt̃∘dfdt∘t⟹t=df̃dt̃∘t̃.
Conversely, when we assume t=(df̃/dt̃)∘t̃, we also obtain t̃=(df/dt)∘t. Thus, (2.4) gives the reciprocal expressions (2.5) and (2.6), which are referred to as variable change.When we assume a finite domainΩ on which a convex function is defined, the independent variable t on the domain can be changed into its conjugate variable t̃ through convex conjugate. Then, the function form f can be changed into the form f̃ by replacing t with t̃. To be precise, the explicit form of the convex conjugate from above relations is(2.7)f̃(t̃)=sup⁡⁡{tt̃-f(t)∣t∈Ω},
where the conjugate variable t̃ is expressed as the gradient of f at t. By differentiating tt̃-f(t) with regard to t and equating this result to zero, we can confirm that t̃ is expressed as the gradient of f. Conversely, the new convex function can be converted back to the original function in the same manner:(2.8)f(t)=sup⁡⁡{t̃t-f(t)∣t̃∈Ω̃}.
In this case, the original variable t is expressed as the gradient of f̃ at t̃. Variables t and t̃ are basically conservative fields with regard to each other. Convex conjugation was originally derived from duality relationship between points and lines. The functional relationship specified by f(t) can be represented equally as well as a set of points t, or as a set of tangent lines specified by their gradients and intercept values.


#### 2.2.2. Nonnegative Conjugate

Now, we introduce a new conjugate transformation termednonnegative conjugate. This transformation is closely related to the former convex conjugate. If f is twice continuously differentiable and the domain is Ω, then we can characterize a convex function as follows:(2.9)fis convex ifff′′(t)≥0forΩ.
This is a link between convex conjugate and nonnegative conjugate based on the following definition.Definition 2.1.
Suppose that two convex functionsf(t) and f̃(t̃) are in convex conjugate for t∈Ω and t̃∈Ω̃ and their second derivativesf′′(t) and f̃′′(t̃) are denoted as I(t) and Ĩ(t̃), respectively. Then I(t) and Ĩ(t̃) are said to be nonnegative conjugate of each other on domains Ω and Ω̃.The two-dimensional conjugate transform that is analogous to this nonnegative conjugate has been applied to image morphing [25].Let us calculate the second derivatives directly. As mentioned in (2.5) and (2.6), the first derivatives represent variable change between t and t̃. The second derivative of f(t) is given as(2.10)f′′(t)=d2f(t)dt2=dt̃dt.
Similarly, the second derivative of f̃′′(t̃) is given as(2.11)f̃′′(t̃)=d2f̃(t̃)dt̃2=dtdt̃.
From (2.10) and (2.11), we obtain the following reciprocal relation between two second derivatives:(2.12)f′′(t)f̃′′(t̃)=1.
If we denote f′′(t) and f̃′′(t̃) as I(t) and Ĩ(t̃), the expression can be rewritten as I(t)Ĩ(t̃)=1. Then, the nonnegative conjugate of  I(t) can be expressed as(2.13)Ĩ(t̃)=1I(t)=1I(df̃(t̃)/dt̃)=1I((d/dt̃)sup⁡⁡{tt̃-f(t)∣t∈Ω})=1I((d/dt̃)sup⁡⁡{tt̃-∫∫I(t)dt∣t∈Ω}).
Like variable change in the convex conjugate transformation, the variable t of the nonnegative function I(t) on the domain Ω can be formally changed by using the nonnegative conjugate. This yields another nonnegative function Ĩ(t̃) on the domain Ω̃ when the variable t is replaced with its conjugate variable t̃. Equation (2.13) has a very complex form, but the variable change between t and t̃ has the following concise forms:(2.14)t̃=∫I(t)dt,t=∫Ĩ(t̃)dt̃.Alternatively, the nonnegative conjugate can also be derived from the equidistribution principle. First, the conjugate variablet̃ is introduced such that a nonnegative distribution I(t) becomes constant with 1 in the conjugate coordinate t̃: I(t)dt=dt̃. The conjugate function Ĩ(t̃) also becomes constant with 1 in the original coordinate t by the same form: Ĩ(t̃)dt̃=dt. As a result, we can obtain the reciprocal relation between I(t) and Ĩ(t̃) and biconjugacy from the equidistribution principle, that is, I(t)Ĩ(t̃)=1 and Ĩ̃=I. This approach is equivalent to solving the Jacobian equation:
(2.15)I(t)=dt̃dt,Ĩ(t̃)=dtdt̃.
Note that equations in (2.15) are the same as (2.10) and (2.11), respectively.


#### 2.2.3. Nonnegative Conjugate of a Nonnegative Vector

In the previous section, we described a method of transformation based on the convex conjugate. However, although the nonnegative conjugate transforms a continuous function into another continuous function, the transformationζ should map an M-dimensional vector onto another M-dimensional vector. Thus, we require a discrete version of the nonnegative conjugate.Let us denote anM-dimensional column vector with nonnegative components as I, that is, I=(I1⋯IM)T and  Ii≥0. Then its nonnegative conjugate is denoted as Ĩ. To transform I into Ĩ, we have to link I with a continuous function I(t) by(2.16)I(t)≡I[t],
where the function  [t]  is a ceiling function that gives the smallest integer not less than t. Then I(t) is a continuous and nonnegative function defined from t=0 tot=M. If we denote a cumulative distribution ∫0tI(t)dt as s(t), we obtain t̃=s(t) from (2.14), and t is given by its inverse:(2.17)t=s-1(t̃).
Then, from (2.13), we obtain(2.18)Ĩ(t̃)=1I(t)=1I[t]=1I[s-1(t̃)].
Finally, we can change it into the M-dimensional vector Ĩ by(2.19)Ĩi=∫i-1iĨ(t̃)dt.̃
Applying the same procedure to Ĩ, we can transform it back to the original vector I. However, this vector is not exactly same as the original vector. As the dimensionality of M increases, the error, I-Ĩ̃, converges to a zero vector.


### 2.3. Application to PPG Signals

All experiments and analyses were performed using ECG and PPG signals extracted from the publically available MIMIC database that contains data from intensive care unit patients admitted to Boston Beth Israel Hospital. ECG and PPG signals were measured to 500 and 125 samples per second, respectively. First,R peaks were detected from ECG based on the assumption that the R peak represents the onset time of left ventricular ejection. Therefore, the position of R peaks was used to segment single PPG pulses from the full PPG signals. Raw PPG signals were low-pass filtered at 15 Hz, then single PPG pulses were separated by synchronized R peaks. The extracted single PPG pulses were resampled to 500 Hz to improve accuracy, and then FOOT and MAX points of single PPG pulses were detected by the traditional method that detects characteristic points [26].Each single PPG pulse was divided into two parts by the FOOT point: the front part from theR peak to the FOOT point, and the rear part from the FOOT point to the next R peak. The time difference at each part was calculated in different ways than that used to estimate PAT. The time difference ta at the front part was derived by a simple translation to change the number of samples into time (seconds), and the time difference tb at the rear part was calculated by the nonnegative conjugate transformation and linear projection, as shown in Figure 2.Figure 2
Single PPG pulse processing to estimate PAT.Various single PPG pulses with different amplitudes, shapes, or pulse widths were represented as single points in anM-dimensional vector space after the nonnegative conjugate transformation. The points that corresponded to nonnegative conjugate vectors, Ĩ, were located on a same line in an M-dimensional vector space. This characteristic is referred to as collinearity.


#### 2.3.1. Training a Projection MatrixW

To derive the matrixW that projects collinear points in an M-dimensional vector space into a one-dimensional time space, 10,000 different PPG pulses were extracted as the training set, and the linear projection matrix was trained according to the MAX point of the PPG pulse, because this is the most obvious characteristic point. Only rear parts of PPG pulses were used to train the linear projection matrix, which we derived by pseudoinverse operation between nonnegative conjugate transformed pulses and the known time information of MAX points as follows:(2.20)tb=WĨor(tb1⋯tbN)=W(Ĩ1⋯ĨN).
The matrix W was calculated from the training samples by using the pseudo-inverse relationship(2.21)W=tbĨT(ĨĨT)-1,
where W=(w1⋯wM)∈ℝ1×M, tb is the known time set from the R peaks to the MAX points, Ĩ is nonnegative conjugate transformed PPG pulse, and W is the derived linear projecting matrix.


#### 2.3.2. PAT Estimation

Single PPG pulses extracted according toR peaks of ECG were divided into front and rear sections using the FOOT point, and time values were calculated for each section separately. First, the time of the front part, ta, was calculated using the number of samples between the R peak and FOOT point divided by the sampling rate. Second, the time of the rear part, tb, was acquired by linear operation between the linear projecting matrix (W) and the nonnegative conjugate transformed pulse (Ĩ). Finally, the PAT was obtained by simple summation of ta and tb:(2.22)ta=Number of samplesSampling rate(seconds),tb=WI,̃tb=(w1⋯wM)(Ĩ1⋮ĨM)(seconds),PAT=ta+tb.


## 2.1. Representation of PPG Signals in Vector Space

LetΩ be the set of all continuous PPG signals measured from human arteries. We can define a sampling process ξM as follows:(2.1)ξM:Ω⟶RM,
where  ℝM is an M-dimensional Euclidean space. The mapping ξM reduces a continuous signal to an M-dimensional vector point. The vector point forms a lower-dimensional cluster in an M-dimensional space. Let us consider the cluster as being embedded by the manifold  Υ. Our goal in this paper is to find a mapping τ between this manifold and PAT:(2.2)τ:Υ⟶R1.
Then, the parametric estimation of PAT is given by the composition of two mappings: PAT=τ∘ξM.However, manifolds constructed through normal sampling at a constant frequency are highly curved. Training this type of manifold and mapping PAT using this manifold are challenging [22]. Let us consider simple translations of Gaussian peaks, as shown in Figure 1. The manifold constructed using simple translations is spirally curved. If slight variation is added to Gaussian peaks, it is not feasible to parameterize the manifold with well-defined functions.Figure 1
The estimation of PAT using a linear transformation approach.Now, suppose that we find the sampling processξM such that the manifold Υ can be flat and isometric along PAT. For instance, if three vectors f,g,h∈Υ are collinear and have isometric timing codestf, tg, th, that is,  (th-tg)f+(tf-th)g+(tg-tf)h=0, then we can always find the linear transformation ω such that tf=ωTf, tg=ωTg, and th=ωTh. This means that the special sampling process allows the mapping τ to be the simplest form by ω. However, we failed to find such a sampling process, regardless of the sampling frequencies applied to the continuous signals; convex combination of two different vectors f,g∈Υ cannot be used to represent the human artery PPG signal.In this context, we propose adding another transformation betweenξM and τ. By considering the new mapping, we intend to keep the mapping, τ, the linear transformation. If we denote the novel mapping as ζ, we can estimate the PAT of the PPG signals as(2.3)PAT=τ∘ζ∘ξM.We refer to this as a linear transformation approach for estimating PAT. In following subsections, we describe the new transformationζ in more detail.


## 2.2. Conjugate Transformations

In the previous subsection, we framed a set of three transformations to change continuous pulse functions into scalar timing codes. The first transformation,ξM, is needed to reduce continuous pulse functions into M-dimensional vector points, while the third transformation τ is the linear transformation. The function of the second transformation,  ζ, will be fully discussed in this subsection, in which we consider a well-known transformation called the convex conjugate or Legendre transformation.


### 2.2.1. Convex Conjugate

The convex conjugate is a transformation that maps a convex function onto another convex function [23, 24]. A convex function always has its conjugate function: the conjugate function is also a convex function. First, we outline why we need convex functions. A Gaussian function, which was exemplified as a pulse in the previous subsection, is the starting point for developing our idea. Gaussian functions have a single peak and are nonnegative over the entire region. Such a Gaussian function can be derived from a convex function by differentiating the convex function twice. We can therefore consider convex functions instead of Gaussian functions. A general type of pulse function that has mixed-signed values, unlike Gaussian functions, will be discussed later.Let us consider two arbitrary convex functions. When their first derivatives become inverses of each other, two functions are referred to as “convex conjugate”. If two functionsf(t) and f̃(t̃) have such a relation, then(2.4)dfdt∘df̃dt̃=I=df̃dt̃∘dfdt,
where I is an identity function, that is, I∘t=t and I∘t̃=t̃. To find an explicit expression for f̃, we assume(2.5)t̃=dfdt∘t.
From (2.4) and (2.5), we obtain(2.6)I∘t=df̃dt̃∘dfdt∘t⟹t=df̃dt̃∘t̃.
Conversely, when we assume t=(df̃/dt̃)∘t̃, we also obtain t̃=(df/dt)∘t. Thus, (2.4) gives the reciprocal expressions (2.5) and (2.6), which are referred to as variable change.When we assume a finite domainΩ on which a convex function is defined, the independent variable t on the domain can be changed into its conjugate variable t̃ through convex conjugate. Then, the function form f can be changed into the form f̃ by replacing t with t̃. To be precise, the explicit form of the convex conjugate from above relations is(2.7)f̃(t̃)=sup⁡⁡{tt̃-f(t)∣t∈Ω},
where the conjugate variable t̃ is expressed as the gradient of f at t. By differentiating tt̃-f(t) with regard to t and equating this result to zero, we can confirm that t̃ is expressed as the gradient of f. Conversely, the new convex function can be converted back to the original function in the same manner:(2.8)f(t)=sup⁡⁡{t̃t-f(t)∣t̃∈Ω̃}.
In this case, the original variable t is expressed as the gradient of f̃ at t̃. Variables t and t̃ are basically conservative fields with regard to each other. Convex conjugation was originally derived from duality relationship between points and lines. The functional relationship specified by f(t) can be represented equally as well as a set of points t, or as a set of tangent lines specified by their gradients and intercept values.


### 2.2.2. Nonnegative Conjugate

Now, we introduce a new conjugate transformation termednonnegative conjugate. This transformation is closely related to the former convex conjugate. If f is twice continuously differentiable and the domain is Ω, then we can characterize a convex function as follows:(2.9)fis convex ifff′′(t)≥0forΩ.
This is a link between convex conjugate and nonnegative conjugate based on the following definition.Definition 2.1.
Suppose that two convex functionsf(t) and f̃(t̃) are in convex conjugate for t∈Ω and t̃∈Ω̃ and their second derivativesf′′(t) and f̃′′(t̃) are denoted as I(t) and Ĩ(t̃), respectively. Then I(t) and Ĩ(t̃) are said to be nonnegative conjugate of each other on domains Ω and Ω̃.The two-dimensional conjugate transform that is analogous to this nonnegative conjugate has been applied to image morphing [25].Let us calculate the second derivatives directly. As mentioned in (2.5) and (2.6), the first derivatives represent variable change between t and t̃. The second derivative of f(t) is given as(2.10)f′′(t)=d2f(t)dt2=dt̃dt.
Similarly, the second derivative of f̃′′(t̃) is given as(2.11)f̃′′(t̃)=d2f̃(t̃)dt̃2=dtdt̃.
From (2.10) and (2.11), we obtain the following reciprocal relation between two second derivatives:(2.12)f′′(t)f̃′′(t̃)=1.
If we denote f′′(t) and f̃′′(t̃) as I(t) and Ĩ(t̃), the expression can be rewritten as I(t)Ĩ(t̃)=1. Then, the nonnegative conjugate of  I(t) can be expressed as(2.13)Ĩ(t̃)=1I(t)=1I(df̃(t̃)/dt̃)=1I((d/dt̃)sup⁡⁡{tt̃-f(t)∣t∈Ω})=1I((d/dt̃)sup⁡⁡{tt̃-∫∫I(t)dt∣t∈Ω}).
Like variable change in the convex conjugate transformation, the variable t of the nonnegative function I(t) on the domain Ω can be formally changed by using the nonnegative conjugate. This yields another nonnegative function Ĩ(t̃) on the domain Ω̃ when the variable t is replaced with its conjugate variable t̃. Equation (2.13) has a very complex form, but the variable change between t and t̃ has the following concise forms:(2.14)t̃=∫I(t)dt,t=∫Ĩ(t̃)dt̃.Alternatively, the nonnegative conjugate can also be derived from the equidistribution principle. First, the conjugate variablet̃ is introduced such that a nonnegative distribution I(t) becomes constant with 1 in the conjugate coordinate t̃: I(t)dt=dt̃. The conjugate function Ĩ(t̃) also becomes constant with 1 in the original coordinate t by the same form: Ĩ(t̃)dt̃=dt. As a result, we can obtain the reciprocal relation between I(t) and Ĩ(t̃) and biconjugacy from the equidistribution principle, that is, I(t)Ĩ(t̃)=1 and Ĩ̃=I. This approach is equivalent to solving the Jacobian equation:
(2.15)I(t)=dt̃dt,Ĩ(t̃)=dtdt̃.
Note that equations in (2.15) are the same as (2.10) and (2.11), respectively.


### 2.2.3. Nonnegative Conjugate of a Nonnegative Vector

In the previous section, we described a method of transformation based on the convex conjugate. However, although the nonnegative conjugate transforms a continuous function into another continuous function, the transformationζ should map an M-dimensional vector onto another M-dimensional vector. Thus, we require a discrete version of the nonnegative conjugate.Let us denote anM-dimensional column vector with nonnegative components as I, that is, I=(I1⋯IM)T and  Ii≥0. Then its nonnegative conjugate is denoted as Ĩ. To transform I into Ĩ, we have to link I with a continuous function I(t) by(2.16)I(t)≡I[t],
where the function  [t]  is a ceiling function that gives the smallest integer not less than t. Then I(t) is a continuous and nonnegative function defined from t=0 tot=M. If we denote a cumulative distribution ∫0tI(t)dt as s(t), we obtain t̃=s(t) from (2.14), and t is given by its inverse:(2.17)t=s-1(t̃).
Then, from (2.13), we obtain(2.18)Ĩ(t̃)=1I(t)=1I[t]=1I[s-1(t̃)].
Finally, we can change it into the M-dimensional vector Ĩ by(2.19)Ĩi=∫i-1iĨ(t̃)dt.̃
Applying the same procedure to Ĩ, we can transform it back to the original vector I. However, this vector is not exactly same as the original vector. As the dimensionality of M increases, the error, I-Ĩ̃, converges to a zero vector.


## 2.2.1. Convex Conjugate

The convex conjugate is a transformation that maps a convex function onto another convex function [23, 24]. A convex function always has its conjugate function: the conjugate function is also a convex function. First, we outline why we need convex functions. A Gaussian function, which was exemplified as a pulse in the previous subsection, is the starting point for developing our idea. Gaussian functions have a single peak and are nonnegative over the entire region. Such a Gaussian function can be derived from a convex function by differentiating the convex function twice. We can therefore consider convex functions instead of Gaussian functions. A general type of pulse function that has mixed-signed values, unlike Gaussian functions, will be discussed later.Let us consider two arbitrary convex functions. When their first derivatives become inverses of each other, two functions are referred to as “convex conjugate”. If two functionsf(t) and f̃(t̃) have such a relation, then(2.4)dfdt∘df̃dt̃=I=df̃dt̃∘dfdt,
where I is an identity function, that is, I∘t=t and I∘t̃=t̃. To find an explicit expression for f̃, we assume(2.5)t̃=dfdt∘t.
From (2.4) and (2.5), we obtain(2.6)I∘t=df̃dt̃∘dfdt∘t⟹t=df̃dt̃∘t̃.
Conversely, when we assume t=(df̃/dt̃)∘t̃, we also obtain t̃=(df/dt)∘t. Thus, (2.4) gives the reciprocal expressions (2.5) and (2.6), which are referred to as variable change.When we assume a finite domainΩ on which a convex function is defined, the independent variable t on the domain can be changed into its conjugate variable t̃ through convex conjugate. Then, the function form f can be changed into the form f̃ by replacing t with t̃. To be precise, the explicit form of the convex conjugate from above relations is(2.7)f̃(t̃)=sup⁡⁡{tt̃-f(t)∣t∈Ω},
where the conjugate variable t̃ is expressed as the gradient of f at t. By differentiating tt̃-f(t) with regard to t and equating this result to zero, we can confirm that t̃ is expressed as the gradient of f. Conversely, the new convex function can be converted back to the original function in the same manner:(2.8)f(t)=sup⁡⁡{t̃t-f(t)∣t̃∈Ω̃}.
In this case, the original variable t is expressed as the gradient of f̃ at t̃. Variables t and t̃ are basically conservative fields with regard to each other. Convex conjugation was originally derived from duality relationship between points and lines. The functional relationship specified by f(t) can be represented equally as well as a set of points t, or as a set of tangent lines specified by their gradients and intercept values.


## 2.2.2. Nonnegative Conjugate

Now, we introduce a new conjugate transformation termednonnegative conjugate. This transformation is closely related to the former convex conjugate. If f is twice continuously differentiable and the domain is Ω, then we can characterize a convex function as follows:(2.9)fis convex ifff′′(t)≥0forΩ.
This is a link between convex conjugate and nonnegative conjugate based on the following definition.Definition 2.1.
Suppose that two convex functionsf(t) and f̃(t̃) are in convex conjugate for t∈Ω and t̃∈Ω̃ and their second derivativesf′′(t) and f̃′′(t̃) are denoted as I(t) and Ĩ(t̃), respectively. Then I(t) and Ĩ(t̃) are said to be nonnegative conjugate of each other on domains Ω and Ω̃.The two-dimensional conjugate transform that is analogous to this nonnegative conjugate has been applied to image morphing [25].Let us calculate the second derivatives directly. As mentioned in (2.5) and (2.6), the first derivatives represent variable change between t and t̃. The second derivative of f(t) is given as(2.10)f′′(t)=d2f(t)dt2=dt̃dt.
Similarly, the second derivative of f̃′′(t̃) is given as(2.11)f̃′′(t̃)=d2f̃(t̃)dt̃2=dtdt̃.
From (2.10) and (2.11), we obtain the following reciprocal relation between two second derivatives:(2.12)f′′(t)f̃′′(t̃)=1.
If we denote f′′(t) and f̃′′(t̃) as I(t) and Ĩ(t̃), the expression can be rewritten as I(t)Ĩ(t̃)=1. Then, the nonnegative conjugate of  I(t) can be expressed as(2.13)Ĩ(t̃)=1I(t)=1I(df̃(t̃)/dt̃)=1I((d/dt̃)sup⁡⁡{tt̃-f(t)∣t∈Ω})=1I((d/dt̃)sup⁡⁡{tt̃-∫∫I(t)dt∣t∈Ω}).
Like variable change in the convex conjugate transformation, the variable t of the nonnegative function I(t) on the domain Ω can be formally changed by using the nonnegative conjugate. This yields another nonnegative function Ĩ(t̃) on the domain Ω̃ when the variable t is replaced with its conjugate variable t̃. Equation (2.13) has a very complex form, but the variable change between t and t̃ has the following concise forms:(2.14)t̃=∫I(t)dt,t=∫Ĩ(t̃)dt̃.Alternatively, the nonnegative conjugate can also be derived from the equidistribution principle. First, the conjugate variablet̃ is introduced such that a nonnegative distribution I(t) becomes constant with 1 in the conjugate coordinate t̃: I(t)dt=dt̃. The conjugate function Ĩ(t̃) also becomes constant with 1 in the original coordinate t by the same form: Ĩ(t̃)dt̃=dt. As a result, we can obtain the reciprocal relation between I(t) and Ĩ(t̃) and biconjugacy from the equidistribution principle, that is, I(t)Ĩ(t̃)=1 and Ĩ̃=I. This approach is equivalent to solving the Jacobian equation:
(2.15)I(t)=dt̃dt,Ĩ(t̃)=dtdt̃.
Note that equations in (2.15) are the same as (2.10) and (2.11), respectively.


## 2.2.3. Nonnegative Conjugate of a Nonnegative Vector

In the previous section, we described a method of transformation based on the convex conjugate. However, although the nonnegative conjugate transforms a continuous function into another continuous function, the transformationζ should map an M-dimensional vector onto another M-dimensional vector. Thus, we require a discrete version of the nonnegative conjugate.Let us denote anM-dimensional column vector with nonnegative components as I, that is, I=(I1⋯IM)T and  Ii≥0. Then its nonnegative conjugate is denoted as Ĩ. To transform I into Ĩ, we have to link I with a continuous function I(t) by(2.16)I(t)≡I[t],
where the function  [t]  is a ceiling function that gives the smallest integer not less than t. Then I(t) is a continuous and nonnegative function defined from t=0 tot=M. If we denote a cumulative distribution ∫0tI(t)dt as s(t), we obtain t̃=s(t) from (2.14), and t is given by its inverse:(2.17)t=s-1(t̃).
Then, from (2.13), we obtain(2.18)Ĩ(t̃)=1I(t)=1I[t]=1I[s-1(t̃)].
Finally, we can change it into the M-dimensional vector Ĩ by(2.19)Ĩi=∫i-1iĨ(t̃)dt.̃
Applying the same procedure to Ĩ, we can transform it back to the original vector I. However, this vector is not exactly same as the original vector. As the dimensionality of M increases, the error, I-Ĩ̃, converges to a zero vector.


## 2.3. Application to PPG Signals

All experiments and analyses were performed using ECG and PPG signals extracted from the publically available MIMIC database that contains data from intensive care unit patients admitted to Boston Beth Israel Hospital. ECG and PPG signals were measured to 500 and 125 samples per second, respectively. First,R peaks were detected from ECG based on the assumption that the R peak represents the onset time of left ventricular ejection. Therefore, the position of R peaks was used to segment single PPG pulses from the full PPG signals. Raw PPG signals were low-pass filtered at 15 Hz, then single PPG pulses were separated by synchronized R peaks. The extracted single PPG pulses were resampled to 500 Hz to improve accuracy, and then FOOT and MAX points of single PPG pulses were detected by the traditional method that detects characteristic points [26].Each single PPG pulse was divided into two parts by the FOOT point: the front part from theR peak to the FOOT point, and the rear part from the FOOT point to the next R peak. The time difference at each part was calculated in different ways than that used to estimate PAT. The time difference ta at the front part was derived by a simple translation to change the number of samples into time (seconds), and the time difference tb at the rear part was calculated by the nonnegative conjugate transformation and linear projection, as shown in Figure 2.Figure 2
Single PPG pulse processing to estimate PAT.Various single PPG pulses with different amplitudes, shapes, or pulse widths were represented as single points in anM-dimensional vector space after the nonnegative conjugate transformation. The points that corresponded to nonnegative conjugate vectors, Ĩ, were located on a same line in an M-dimensional vector space. This characteristic is referred to as collinearity.


### 2.3.1. Training a Projection MatrixW

To derive the matrixW that projects collinear points in an M-dimensional vector space into a one-dimensional time space, 10,000 different PPG pulses were extracted as the training set, and the linear projection matrix was trained according to the MAX point of the PPG pulse, because this is the most obvious characteristic point. Only rear parts of PPG pulses were used to train the linear projection matrix, which we derived by pseudoinverse operation between nonnegative conjugate transformed pulses and the known time information of MAX points as follows:(2.20)tb=WĨor(tb1⋯tbN)=W(Ĩ1⋯ĨN).
The matrix W was calculated from the training samples by using the pseudo-inverse relationship(2.21)W=tbĨT(ĨĨT)-1,
where W=(w1⋯wM)∈ℝ1×M, tb is the known time set from the R peaks to the MAX points, Ĩ is nonnegative conjugate transformed PPG pulse, and W is the derived linear projecting matrix.


### 2.3.2. PAT Estimation

Single PPG pulses extracted according toR peaks of ECG were divided into front and rear sections using the FOOT point, and time values were calculated for each section separately. First, the time of the front part, ta, was calculated using the number of samples between the R peak and FOOT point divided by the sampling rate. Second, the time of the rear part, tb, was acquired by linear operation between the linear projecting matrix (W) and the nonnegative conjugate transformed pulse (Ĩ). Finally, the PAT was obtained by simple summation of ta and tb:(2.22)ta=Number of samplesSampling rate(seconds),tb=WI,̃tb=(w1⋯wM)(Ĩ1⋮ĨM)(seconds),PAT=ta+tb.


## 2.3.1. Training a Projection MatrixW

To derive the matrixW that projects collinear points in an M-dimensional vector space into a one-dimensional time space, 10,000 different PPG pulses were extracted as the training set, and the linear projection matrix was trained according to the MAX point of the PPG pulse, because this is the most obvious characteristic point. Only rear parts of PPG pulses were used to train the linear projection matrix, which we derived by pseudoinverse operation between nonnegative conjugate transformed pulses and the known time information of MAX points as follows:(2.20)tb=WĨor(tb1⋯tbN)=W(Ĩ1⋯ĨN).
The matrix W was calculated from the training samples by using the pseudo-inverse relationship(2.21)W=tbĨT(ĨĨT)-1,
where W=(w1⋯wM)∈ℝ1×M, tb is the known time set from the R peaks to the MAX points, Ĩ is nonnegative conjugate transformed PPG pulse, and W is the derived linear projecting matrix.


## 2.3.2. PAT Estimation

Single PPG pulses extracted according toR peaks of ECG were divided into front and rear sections using the FOOT point, and time values were calculated for each section separately. First, the time of the front part, ta, was calculated using the number of samples between the R peak and FOOT point divided by the sampling rate. Second, the time of the rear part, tb, was acquired by linear operation between the linear projecting matrix (W) and the nonnegative conjugate transformed pulse (Ĩ). Finally, the PAT was obtained by simple summation of ta and tb:(2.22)ta=Number of samplesSampling rate(seconds),tb=WI,̃tb=(w1⋯wM)(Ĩ1⋮ĨM)(seconds),PAT=ta+tb.


## 3. Results and Discussion

To assess the agreement between the traditional method and our novel PAT estimation method, we evaluated a subset of data from the MIMIC database. This database is part of the Physionet platform and contains data from over 72 intensive care unit patients at the Boston Beth Israel Hospital, but we selected only those records for which nonsaturated and nonmissing ECG and PPG signals were simultaneously measured and available [27, 28].We adopted the agreement analysis proposed by Bland and   Altman: given two different estimating methods, their agreement can be assessed by computing the standard deviation of two sets of estimates. The proposed strategy computes the differences between measurements provided by two methods and then computes their dispersion. Two methods have good agreement if dispersion is minimal [29].To acquire PAT estimates using our method, a linear projecting matrix was derived using a training process, and the derived linear projection matrix was applied to two test sets consisting of 2947 and 2890 PPG pulses, respectively. The mean difference and standard deviations of two sets of PAT estimates were calculated. Results of the agreement analysis are shown in Figure3. The 95% limits of agreement were calculated as mean difference  ±1.96*  standard deviation at each set.Agreement analysis between the two different PAT estimations.
(a)(b)For the first test set, 98.9% of pulses were located between −34.3 ms and 28.3 ms as the limits of agreement, while in the second test set, 96.2% of the pulses were located between −33.4 ms and 25.3 ms.The typical PAT estimation method, which detects characteristic points, is not accurate when applied to PPG pulse types with different morphologies. We therefore proposed a novel PAT estimation method that provides robust results by considering the morphological characteristics of PPG pulses according to the properties of blood vessels. We initially attempted to find a relation between a projecting factor and original PPG pulses,f, but were unsuccessful because of the broad dispersion of pulses in an M-dimensional vector space. We therefore decided to transform original PPG pulses into another form. We used the nonnegative conjugate transformation, because the nonnegative conjugate transformed signal,  Ĩ, had the property of collinearity in an M-dimensional vector space and could be used to estimate the PAT by projection onto a one-dimensional time space. We derived a linear projection matrix through the training process for linear projection and applied this matrix to two different morphological PPG pulse sets. PAT values estimated from annotated MAX points and the proposed linear projecting method were in good agreement; over 95% of the data were included within the 95% limits of agreement.Although our method provides results that appear to be highly accurate, it can show different results according to the linear projection matrix that is derived by different numbers of pulses and dimensionality. Therefore, an optimal combination that is applicable to a variety of morphological PPG pulse types should be determined. For instance, the size of the projecting matrix needs to be adjusted and the time delay caused by the nonnegative conjugate transformation needs to be addressed. Our novel approach still has some limitations in terms of its clinical application for real-time continuous monitoring of PAT as well as stiffness and blood pressure assessment using PPG; the linear projecting matrix needs to be optimized and the time delay caused by the nonnegative conjugate transformation needs to be addressed. Once these issues are addressed, however, our method has great potential in clinical practice to precisely assess cardiovascular risk associated with blood vessels.


## 4. Conclusions

Various PAT estimation methods exist, most of which are based on unsupervised extraction of characteristic points in PPG signals. Despite the good performance of these PAT estimation techniques when applied to clean PPG signals, they are less reliable when used to analyze morphologically variable PPG signals. Thus, we designed a novel, simple linear model based on the nonnegative conjugate transformation. This easy, stable PAT estimation method relies on training of the linear model using various samples. Because our method extracts information from various pressure pulse, it can be applied to different morphological signals without special conditions. In conclusion, we developed a novel method that can be used to estimate PAT robustly for a variety of PPG signals with different morphological characteristics.



---

*Source: 643653-2012-04-11.xml*