INTRINSIC FREQUENCY ANALYSIS AND FAST ALGORITHMS

arXiv:1708.00465v1 [math.NA] 1 Aug 2017

PEYMAN TAVALLALI1,∗ , HANA KOOREHDAVOUDI2, JOANNA KRUPA3

Abstract. Intrinsic Frequency (IF) has recently been introduced as an ample signal processing
method for analyzing carotid and aortic pulse pressure tracings. The IF method has also been
introduced as an effective approach for the analysis of cardiovascular system dynamics. The physiological significance, convergence and accuracy of the IF algorithm has been established in prior
works. In this paper, we show that the IF method could be derived by appropriate mathematical approximations from the Navier-Stokes and elasticity equations. We further introduce a fast
algorithm for the IF method based on the mathematical analysis of this method. In particular,
we demonstrate that the IF algorithm can be made faster, by a factor or more than 100 times,
using a proper set of initial guesses based on the topology of the problem, fast analytical solution
at each point iteration, and substituting the brute force algorithm with a pattern search method.
Statistically, we observe that the algorithm presented in this article complies well with its bruteforce counterpart. Furthermore, we will show that on a real dataset, the fast IF method can draw
correlations between the extracted intrinsic frequency features and the infusion of certain drugs. In
general, this paper aims at a mathematical analysis of the IF method to show its possible origins
and also to present faster algorithms.

1 Division

of Engineering and Applied Sciences, California Institute of Technology, 1200 East Cali-

fornia Boulevard, MC 205-45, Pasadena, CA 91125, USA
2 Aerospace

and Mechanical Engineering, University of Southern California, Los Angeles, CA 90089-

1453, USA
3 Avicena

LLC, 2400 N Lincoln Ave, Altadena, CA 91001, USA

∗ Corresponding

Author, email: ptavalla@caltech.edu, tavallali@gmail.com

1. Introduction
Cardiovascular diseases (CVDs) and stroke are major causes of death in the United States. The
total cost related to CVDs and stroke was estimated to be more than $316 billion in 2011-2012
1

[20, 22]. Hence, clinical measurements of cardiovascular health indices are of great importance.
These methods and measurements are essential tools for monitoring cardiovascular health due to
their relative availability. For example, Left Ventricular Ejection Fraction (LVEF) is a measure of
left ventricular contractility [6] and Carotid-Femoral Pulse Wave Velocity (cfPWV) is a measure of
aortic stiffness [18].
However, current methods of measuring such indices are expensive, sometimes invasive, prone to
measurement errors, and not necessarily easy to use. For example, 2D LVEF echocardiography is not
accurate compared to more expensive and laborious gold standard cardiac MRI method [10, 12, 15,
14]. As another example, obtaining accurate cfPWV measurements often requires certain medical
devices and a well-trained staff within a clinical setting [29]. Consequently, continuous measurement
of these indices is not practical. These limitations emphasize the need for new cardiovascular
monitoring methods.
Intrinsic Frequency (IF) has been established as a new method of cardiovascular monitoring through
a novel signal processing methodology [26]. The IF method needs only an uncalibrated pulse pressure
[31] to extract pertinent information regarding the cardiovascular health of an individual [26]. The
IF method has also been shown to be capable of non-invasively measuring LVEF by means of an
iPhone camera [25]. We believe that methods like IF are of clinical and financial benefit in addressing
cardiovascular monitoring.
In this paper, at first, we provide an overview of the IF method. Next, we present an approximate
derivation of the IF model by combining Navier-Stokes equations and continuity with elasticity
equations. This helps to build a solid mathematical foundation for the IF method and the analysis
that follows. Later, we analyze the IF algorithm in the space of feasible solutions, and based on
that, we introduce a new version of the IF algorithm which is faster than the current brute-force IF
method [31] while maintaining the same accuracy. We then perform a case study on real pressure
waveforms drawn from canine data using our new algorithm. We will see that the fast IF algorithm
is capable of capturing the effects of different drug infusions on a canine subject.

2

2. Brief Overview of IF method
2.1. A History of Analyzing Cardiovascular Pulse Waveform. Blood pressure was first measured by Hales in 1735 [13]. In his measurements, he found that blood pressure is not constant in
the arterial system. He related these variations to the elasticity of the arteries [13]. Currently,
it is known that the shape of the arterial pulse wave is intimately related to the physiology and
pathology of the whole arterial system [1]. There has been much research on analyzing the dynamics
of blood pressure and flow in arterial systems [2, 9, 23, 24, 33]. Specifically, there are two main
approaches to analyzing cardiovascular pulse wave data. One approach is based on a systematic
mathematical framework for the cardiovascular system. The other is based on directly analyzing
the pulse pressure waveform using signal processing methods.
An example of the systematic framework can be seen with the set of Windkessel models [34]. The
formulation of a minimal lumped model of the arterial system was first presented by Westerhof et
al. [34]. Based on a Windkessel model, the arterial system dynamics have been modeled through a
combination of different elements such as resistance, compliance and impedance. In this simplified
model of the arterial system, the blood flow dynamics is modeled by the interaction between the
elements (assuming the blood flow acts as the current in the system). Because of the type of
modeling, the wave transmission of the blood flow is neglected. As a result, the Windkessel models
is not able to represent the entire dynamics of the blood flow in an arterial system accurately.
On the other hand, there are various methods for direct analysis of an arterial pulse waveform, in
both time and frequency domains [24]. For example, the impedance method, which is based on
Fourier transform, is a common method to analyze the pressure waveform in the frequency domain
[2]. As an example, Milnor has shown that the pressure and flow waveforms can be a superposition
of several harmonics using the Fourier method [21]. Another method to investigate the pressure
wave in the time domain is the wave intensity analysis which is based on wavelet transform [8].
These methods do not necessarily convey a physical understanding of the cardiovascular system.
The IF algorithm presented in [31] is analyzing a pulse waveform through a direct time-frequency
signal processing machinery setting, from a quantitative perspective. Although, in previous work
[26], we tried to qualitatively express a systems approach to the IF formulation, the quantitative

3

picture has not yet been clear. However, in this article, we show this connection from a quantitative
perspective.
2.2. IF Formulation. In the IF method, the aortic pressure waveform at time t ∈ [0, T ), for a
cardiac period T , can be represented as

(2.1)

S (ai , bi , p̄, ωi ; t) = (a1 cos ω1 t + b1 sin ω1 t + p̄) 1[0,T0 ) (t) +
(a2 cos ω2 t + b2 sin ω2 t + p̄) 1[T0 ,T ) (t) ,

with a continuity condition at T0 and periodicity at T . In this formulation, the indicator function
is defined as


 1, x ≤ t < y,
1[x,y) (t) =
 0,
else.

Here, a1 , b1 , a2 and b2 are the envelopes of the IF model fit. ω1 and ω2 are the Intrinsic Frequencies
(IFs) of the waveform. Further, p̄ is the mean pressure during the cardiac cycle. This type of
formulation embeds the coupling and decoupling of heart and aorta.
The goal of the IF model (2.1) is to extract a fit, called Intrinsic Mode Function (IMF), that carries
most of the energy (information) from an observed pressure waveform f (t) in one period. The latter
is done by solving the following optimization problem [31]:
minimize kf (t) − S (ai , bi , p̄, ωi ; t)k22

(2.2)

(2.3)

ai ,bi ,ωi ,p̄

subject to

a1 cos ω1 T0 + b1 sin ω1 T0 = a2 cos ω2 T0 + b2 sin ω2 T0 ,
a1

=

a2 cos ω2 T + b2 sin ω2 T.

Here, kk2 is the L2 -norm. The first linear condition in this optimization enforces the continuity of
the extracted IMF at the dicrotic notch. The second one imposes the periodicity. The mathematical
convergence and accuracy of the IF algorithm have been explained in a previous work [31]. In the
next sections, we explore the foundation of the IF algorithm and propose a faster IF algorithm.

4

3. Approximate Derivation of the IF Model
As mentioned earlier, in a previous work [26], we tried to express a systems approach to the IF
formulation qualitatively. However, in this article, we show this connection from a quantitative
perspective. This section is devoted to this purpose.
In this paper, we assume that the Left Ventricle (LV), the aortic valve, aorta and the arterial system
can be represented by a simplified model as shown in Figure 1. Here, the LV and the aortic valve
are assumed to be the boundary condition at the entrance of the aortic tube and the arterial system
is the terminal boundary condition of the aortic tube. The boundary condition at the entrance of
the aortic tube changes from an LV boundary condition to a closed valve boundary condition, at the
dicrotic notch time T0 during a cardiac cycle [0, T ]. We further assume that blood is a Newtonian
incompressible fluid, the aorta is a straight and sufficiently long elastic tube with a constant circular
cross section and there is no external force causing flow rotation. These assumptions are not all
satisfied in a real cardiovascular system. However, they are useful in estimating the general behavior
of blood in aorta.
Combining the Navier-Stokes equations and continuity with the elasticity equation, we can drive a
model for the flow Q (x, t) and the pressure P (x, t) along the length x of an aorta as follow

(3.1)

(3.2)

−

∂P
∂Q
(x, t) = L
(x, t) + RQ (x, t) ,
∂x
∂t

−

∂Q
∂P
(x, t) = C
(x, t) .
∂x
∂t

The step by step derivation of these equations is presented in Appendix A. Parameters L, R, and
C represent inductance, resistance, and compliance of the blood in aorta. Here, 0 ≤ x ≤ h, where
h represents the aortic length. This model has also been discussed and simulated numerically in [3]
with a complex set of boundary conditions. Here, our main concentration will be on the aortic tube
oscillatory waveform solutions. Next, we will show that we can derive (2.1) from (3.1) and (3.2).

5

Since the input to the IF model (2.1) is a pressure waveform, we need to extract an equation for
the pressure P (x, t) from Equations (3.1) and (3.2) by eliminating the flow. Combining Equations
(3.1) and (3.2) results in

(3.3)

CL

∂P
∂2P
∂2P
(x,
t)
+
CR
(x,
t)
=
(x, t) .
∂t2
∂t
∂x2

Taking P (x, t) = K (t) p (x, t) + p̄, with p̄ as the the mean pressure, we can write equation (3.3) as

(3.4)





∂2p
CLK̈ (t) + CRK̇ (t) p (x, t) + 2CLK̇ (t) + CRK (t) ∂p
∂t (x, t) + CLK (t) ∂t2 (x, t)
2

∂ p
= K (t) ∂x
2 (x, t) .

Here, we have used the dot notation to represent the time derivative. We can simplify the term in
front of
R
− 2L
t

Ke

∂p
∂t

(x, t), in (3.4), by setting 2CLK̇ (t) + CRK (t) = 0. The latter has a solution K (t) =

for some constant K. This reduces Equation (3.4) into

(3.5)

CL

∂2p
CR2
∂2p
p
(x,
t)
=
(x,
t)
−
(x, t) .
∂t2
4L
∂x2

The solution of Equation (3.5) can be expressed in terms of eigenfunctions. In other words, using
the method of separation of the variables, one can express the solution of Equation (3.5) as
(3.6)

p (x, t) =

∞
X

Tn (t) Xn (x)

n=1

for

(3.7)

(3.8)

Tn (t) = αn sin (ωn t) + βn cos (ωn t) ,

Xn (x) = ζn sin

r

CR2
x
CL (ωn ) −
4L
2

6

!

+ ηn cos

r

!
CR2
x ,
CL (ωn ) −
4L
2

and some constants αn , βn , ζn and ηn . As a result, the solution of (3.3) can be expressed as
(
R P∞
P (x, t) = p̄ + Ke− 2L t n=1 (αn sin (ωn t) + βn cos (ωn t))
(3.9)

q
)

q
CL (ωn )2 −

ζn sin

CR2
4L x

CL (ωn )2 −

+ ηn cos

CR2
4L x

.

The variables ωn can be expressed based on the boundary conditions of the aortic tube. We need
to emphasize that for a period of the cardiac cycle [0, T ), the boundary conditions change before
and after the dicrotic notch T0 . Hence, for t ∈ [0, T ), Equation (3.9) can be written as

(3.10)



P (x, t) = p̄
(


R P∞
+1[0,T0 ) (t) K 1 e− 2L t n=1 α1n sin ωn1 t + βn1 cos ωn1 t

ζn1 sin

q

CL (ωn1 )2 −

CR2
4L

R

+1[T0 ,T ) (t) K 2 e− 2L t


ζn2 sin

q

CL (ωn2 )2

−



x + ηn1 cos
(

P∞



CL (ωn1 )2 −

CR2
4L

x

)



α2n sin ωn2 t + βn2 cos ωn2 t

n=1

CR2
4L x

q

+

ηn2 cos

q

CL (ωn2 )2

−

CR2
4L x

)
.

Here, the superscripts indicated with “1” belong to the form of the solution before the closure of
the aortic valve, and the superscripts indicated with “2” belong to the form of the solution after
the closure of the aortic valve. Constants K 1 , α1n , β 1n , ζ 1n , η 1n and ω 1n are found from the boundary
and initial conditions at systole. Similarly, constants K 2 , α2n , β 2n , ζ 2n , η 2n and ω 2n are found from the
boundary and initial conditions at diastole.
Equation (3.10) is explicitly showing the coupling and decoupling of heart and aorta before and
after the dicrotic notch. As the boundary conditions change during a cardiac cycle, the frequencies
of oscillation also change from ωn1 to ωn2 . Generally, Equation (3.9) can represent pressure waveform
for a Newtonian incompressible fluid in a straight and sufficiently long elastic tube with constant
circular cross section.

7

If the pressure is recorded at a specific point x0 on aorta, the terms containing the spacial variable
x would be fixed. In other words, Equation (3.10) would reduce to
P (x = x0 , t) = p̄
n
o

R P∞
+ K 1 e− 2L t n=1 κ1n α1n sin ωn1 t + βn1 cos ωn1 t
1[0,T0 ) (t)
n
o


P
R
2
2
2
2
2
+ K 2 e− 2L t ∞
1[T0 ,T ) (t) ,
n=1 κn αn sin ωn t + βn cos ωn t

(3.11)

for

ζn1 sin

r

κ2n = ζn2 sin

r

κ1n

=

and

CL (ωn1 )2

CR2
−
x0
4L

CR2
x0
CL (ωn2 )2 −
4L

!

!

ηn1 cos

r

+ ηn2 cos

r

+

CL (ωn1 )2

CR2
−
x0
4L

CR2
x0
CL (ωn2 )2 −
4L

!

!

.

Now, considering that the cardiac cycle length would be around 1.5 sec, at most, and taking into
R

account that R is smaller than L [3], one can use the approximation e− 2L t ≃ 1. Hence, Equation
(3.11) would become

(3.12)

P (x = x0 , t) ≈ p̄


1
1
1
1
1
+ K
1[0,T0 ) (t)
n=1 κn αn sin ωn t + βn cos ωn t
 2 P∞ 2 2


2
2
2
+ K
1[T0 ,T ) (t) .
n=1 κn αn sin ωn t + βn cos ωn t


P∞
1

Further, if most of the information, or energy, is carried out by the first terms in the series of the
solution, we can further write the approximated solution (3.12) as

(3.13)

Now, by relabeling

P (x = x0 , t) ≈ p̄


 1 1 1
1[0,T0 ) (t)
+ K κ1 α1 sin ω11 t + β11 cos ω11 t


 2 2 2
1[T0 ,T ) (t) .
+ K κ1 α1 sin ω12 t + β12 cos ω12 t

(3.14)

b1 = K 1 κ11 α11 ,

(3.15)

a1 = K 1 κ11 β11 ,

(3.16)

b2 = K 2 κ21 α21 ,

8

(3.17)

a2 = K 2 κ21 β12 ,

(3.18)

ω1 = ω11 ,

(3.19)

ω2 = ω12 ,

we can approximate the IF model (2.1). The continuity and periodicity conditions (2.3) can also be
approximated if we hold the assumption that most of the energy is carried out by the first terms in
the series of the solution (3.12).
In short, in this section, we have presented an approximate quantitative justification on the origins
of the IF method. In the next section, we move on with the analysis of the optimization problem
(2.2) subject to (2.3).

4. Analysis of The IF Algorithm
Practically, one must solve the discrete version of (2.2). We assume that the pressure waveform f (t)
is sampled uniformly. Also, we can simplify (2.2) by the fact that any sinusoid can be assumed to
start from time t = 0 with a compensation coming from a phase shift. In other words, any sinusoid
can be expressed as A cos ωt + B sin ωt, irrespective of whether the initial time is t = 0 or t = T0 .
Hence, the discrete format of (2.2) can be expressed as
kf − S (ai , bi , p̄, ωi ; t)k22

minimize
ai ,bi ,ωi ,p̄

(4.1)
subject to

a1 cos ω1 T0 + b1 sin ω1 T0 =

a2

,

= a2 cos ω2 (T − T0 ) + b2 sin ω2 (T − T0 ) ,

a1

for f = (f1 , . . . , fn+m )′ as the uniform sampling of the original cycle. Here, by taking
(4.2)

t = t′1 , t′2

′

= t11 , t21 , . . . , tn1 , t12 , t22 , . . . , tm
2

9

′

∈ R(n+m)×1

for t1 = (0, ∆t, 2∆t, . . . , T0 )′ ∈ Rn×1 and t2 = (∆t, 2∆t, . . . , T − T0 )′ ∈ Rm×1 , we have the discrete
form of S (ai , bi , p̄, ωi ; t) as


(4.3)

S (ai , bi , p̄, ωi ; t) = 

a1 cos ω1 t1 + b1 sin ω1 t1
a2 cos ω2 t2 + b2 sin ω2 t2



 + p̄1.

In this article, (.)′ denotes the transpose operator and the vector 1 = (1, 1, . . . , 1)′ ∈ R(n+m)×1 .
Also,
′
cos ω1 t1 = cos ω1 t11 , . . . , cos ω1 tn1 ,
′
sin ω1 t1 = sin ω1 t11 , . . . , sin ω1 tn1 ,
′
,
cos ω2 t2 = cos ω2 t12 , . . . , cos ω2 tm
2

′
.
sin ω2 t2 = sin ω2 t12 , . . . , sin ω2 tm
2

(4.4)

The constraints, in (4.1), can be written as

(4.5)



cos ω1 T0

−1

sin ω1 T0

1

− cos ω2 (T − T0 )

0






a
1 


  
 a2 
0
  0 

.

=
 b 
− sin ω2 (T − T0 )
0
 1 


b2

If we can solve for two, out of four, unknowns in (4.5), we would make (4.1) an unconstrained
optimization. However, it is important to check whether the matrix in (4.5) is of full rank or not.
In fact, the rows of this matrix are linearly independent except when
cos ω1 T0 cos ω2 (T − T0 ) = 1.

(4.6)
This will lead into two cases:

(1) Degenerate Case in which Equation (4.6) holds,
(2) General Case in which, it does not.
4.1. General Case (cos ω1 T0 cos ω2 (T − T0 ) 6= 1). One can solve the constraints in (4.1) for a1
and a2 to obtain
(4.7)

a1 =

b1 sin ω1 T0 cos ω2 (T − T0 ) + b2 sin ω2 (T − T0 )
,
1 − cos ω1 T0 cos ω2 (T − T0 )
10

(4.8)

a2 =

b1 sin ω1 T0 + b2 cos ω1 T0 sin ω2 (T − T0 )
.
1 − cos ω1 T0 cos ω2 (T − T0 )

Equations (4.7) and (4.8) would then simplify (4.3) into
(4.9)

S (ω1 , ω1 , b1 , b2 , p̄; t) = Q (ω1 , ω1 , b1 , b2 ; t) + p̄1,

where Q (ω1 , ω1 , b1 , b2 ; t) = b1 v1 (ω1 , ω2 ; t) + b2 v2 (ω1 , ω2 ; t) for
(4.10)



sin ω1 T0 cos ω2 (T −T0 )
1−cos ω1 T0 cos ω2 (T −T0 )

cos ω1 t1 + sin ω1 t1

sin ω1 T0
1−cos ω1 T0 cos ω2 (T −T0 )

cos ω2 t2



sin ω2 (T −T0 )
1−cos ω1 T0 cos ω2 (T −T0 ) cos ω1 t1
cos ω1 T0 sin ω2 (T −T0 )
1−cos ω1 T0 cos ω2 (T −T0 ) cos ω2 t2 +

v1 (ω1 , ω2 ; t) = 



,

and
(4.11)

v2 (ω1 , ω2 ; t) = 

sin ω2 t2



.

Using Equations (4.9)-(4.11), and dropping the dependencies in notation, simplifies (4.1) into
minimize kQ + p̄1 − f k22 .

(4.12)

ω1 ,ω2 ,b1 ,b2 ,p̄

This simplification has helped to eliminate the constraints in the optimization problem (4.1).
The minimization problem (4.12) is non-convex and non-linear in its parameters. So, in order to be
able to solve the problem, we can use the fact that the minimum of a function can first be found
over some variables and then over the remaining ones [7]. In other words, the optimization problem
in (4.12) can be written as
(4.13)

minimize
ω1 ,ω2



minimize kQ + p̄1 −
b1 ,b2 ,p̄

f k22



.

We call the inner optimization in (4.13) as P (ω1 , ω2 ). Solving for P (ω1 , ω2 ) is a classical least
squares problem. The solution existence and uniqueness of this optimization is mentioned in our
previous work [31]. To find the exact solution we simplify the objective function as

(4.14)

kQ + p̄1 − f k22 =

(Q + p̄1 − f )′ (Q + p̄1 − f )

= Q′ Q + 2p̄Q′ 1 − 2Q′ f − 2p̄f ′ 1 + p̄1′ 1 + f ′ f .

11

Substituting for Q = b1 v1 + b2 v2 , we convert (4.14) into
(4.15)

kQ + p̄1 − f k22 = b21 v1′ v1 + 2b1 b2 v1′ v2 + b22 v2′ v2 + 2p̄b1 v1′ 1 + 2p̄b2 v2′ 1
−2b1 v1′ f − 2b2 v2′ f − 2p̄f ′ 1 + p̄2 1′ 1 + f ′ f .

Since, in this part of the optimization, the values of ω1 and ω2 are fixed, we can find the optimal
values of b1 , b2 , and p̄ by setting the partial derivatives of (4.15) equal to zero. In other words, we
∂ (kQ+p̄1−f k22 )
∂ (kQ+p̄1−f k22 )
∂ (kQ+p̄1−f k22 )
set
=
0,
=
0,
and
= 0. Doing this, we find the optimal
∂b1
∂b2
∂ p̄
solution for b1 , b2 , and p̄, by

(4.16)



b∗1 (ω1 , ω2 )


 ∗
 b2 (ω1 , ω2 )

p̄∗ (ω1 , ω2 )





v1′ v1 v1′ v2 v1′ 1

 
  ′
 =  v1 v2 v2′ v2 v2′ 1
 
v1′ 1 v2′ 1 1′ 1

−1 





v1′ f


 ′
 v2 f

1′ f





.


Here, we have fulfilled the optimization part by solving a linear system. This could potentially
accelerate the IF algorithm. Finally, we only have to solve a minimization on
(4.17)

P (ω1 , ω2 ) = kQ (ω1 , ω2 , b∗1 (ω1 , ω2 ) , b∗2 (ω1 , ω2 ) ; t) + p̄∗ (ω1 , ω2 ) 1 − f k22 ,

which is
(4.18)

minimize P (ω1 , ω2 ) .
ω1 ,ω2

We note that a property of the function P (ω1 , ω2 ) is its differentiability, away from its singularities.
In fact, by definition, the function kQ + p̄1 − f k22 is directionally differentiable with respect to all
its variables. Hence, using the results in [5, 28], we can deduce that
(4.19)

P (ω1 , ω2 ) = minimize kQ + p̄1 − f k22
b1 ,b2 ,p̄

is directionally differentiable with respect to ω1 and ω2 . This property can be exploited if one tries
to solve (4.18) using a gradient based optimization method [4].
4.2. Degenerate Case (cos ω1 T0 cos ω2 (T − T0 ) = 1). The solution of (4.6) can be expressed as
nodes of a lattice N in ω1 ω2 plane. To be more specific, we have
(4.20)

N = Γ1 ∪ Γ2 ,
12

where
(4.21)

Γ1 = {(ω1 , ω2 ) |ω1 T0 = (2k1 + 1) π, ω2 (T − T0 ) = (2k2 + 1) π, k1 ∈ Z, k2 ∈ Z } ,

and
(4.22)

Γ2 = {(ω1 , ω2 ) |ω1 T0 = 2k1 π, ω2 (T − T0 ) = 2k2 π, k1 ∈ Z, k2 ∈ Z } .

If (ω1 , ω2 ) ∈ Γ1 , from (4.5) we have a1 = −a2 . On the other hand, if (ω1 , ω2 ) ∈ Γ2 , from (4.5) we
have a1 = a2 . In both of these cases, we can express (4.3) as
(4.23)

S (ω1 , ω1 , a1 , b1 , b2 , p̄; t) = Q (ω1 , ω1 , a1 , b1 , b2 ; t) + p̄1,

where Q (ω1 , ω1 , a1 , b1 , b2 ; t) = a1 w0Γi (ω1 , ω2 ; t) + b1 w1 (ω1 , ω2 ; t) + b2 w2 (ω1 , ω2 ; t), for i = 1, 2. If
(ω1 , ω2 ) ∈ Γ1 ,

(4.24)



w0Γ1 = 



cos ω1 t1

.

− cos ω2 t2

Similarly, if (ω1 , ω2 ) ∈ Γ2 , we have

(4.25)



cos ω1 t1

w0Γ2 = 

cos ω2 t2



sin ω1 t1



02



.

In both of the cases, we have

(4.26)

w1 = 

01



,

and
(4.27)

w2 = 

sin ω2 t2



.

Here, 01 and 02 are zero vectors in Rm×1 and Rn×1 , respectively. It is clear, from (4.26) and (4.27),
that w1′ w2 = w2′ w1 = 0. Using (4.23), and a similar approach we employed in (4.15) and (4.16),

13

we find the optimal solution for a1 , b1 , b2 , and p̄, by
(4.28)


a∗1,i (ω1 , ω2 )



 b∗ (ω1 , ω2 )
 1,i

 b∗ (ω , ω )
 2,i 1 2

p̄∗,i (ω1 , ω2 )




′

′

′
′
w0Γi w0Γi
w0Γi 1
w0Γi w2
w0Γi w1

′
w0Γi w1
0
w1′ 1
w1′ w1

′
w0Γi w2
w2′ 1
0
w2′ w2

′
w0Γi 1
w1′ 1
w2′ 1
1′ 1

 

 
 
 
 
=
 
 
 

−1  
′
Γi
w
f
0
 
 
  w′ f
 
1
 
  w′ f
 
2
 
1′ f







,




for i = 1, 2. Hence, similar to (4.17), for (ω1 , ω2 ) ∈ Γ1 or (ω1 , ω2 ) ∈ Γ2 , we only have to solve a
minimization on
(4.29)


P (ω1 , ω2 ) = Q ω1 , ω1 , a∗1,i (ω1 , ω2 ) , b∗1,i (ω1 , ω2 ) , b∗2,i (ω1 , ω2 ) ; t + p̄∗,i (ω1 , ω2 ) 1 − f

2
.
2

Note that, from a machine learning perspective, the nodes specified in (4.20) do not have important
information physiologically as they could be inferred from the systolic and diastolic parts of a waveform alone. In other words, even if these points present a global minima, they are not informative
as we already know the systolic and diastolic inverses,

1
T0

and

1
T −T0

respectively, as possible inputs

to any machine learning algorithm. Hence, these points could possibly be ignored in a search for an
optimum point of (4.1).

5. Fast IF Algorithms
In this section, we present a fast IF algorithm which is based on the results presented in the previous
section and the topology of the solution space for P (ω1 , ω2 ). In order to keep the fluency of this
section, we mention the original IF algorithm (see Algorithm 1) as presented in [31].
Algorithm 1 has three major steps. In the first step, the (ω1 , ω2 ) domain
(5.1)

Df r = {(ω1 , ω2 ) |0 < ω1 ≤ C, 0 < ω2 ≤ C }

is made discrete, namely D f r . The second step is a minimization to find P (ω1 , ω2 ), see (4.18). The
final step is a brute-force search on D f r to find the minimum of P (ω1 , ω2 ).
All three steps can be optimized to make the IF algorithm faster. Regarding the domain of optimization Df r , defined in (5.1), we know from our previous work in [25] that the average IF solution,
14

for a physiological pulse waveform recording, is confined to a smaller domain D expressed as
(5.2)



ω1 T0
ω2 (T − T0 )
D = (ω1 , ω2 ) 0.5 6
6 1.5, 0.5 6
63 .
π
π

This will make the first step search area more well-defined and optimized. In the previous section,
we have been able to find some analytic solutions (see (4.16)) for the inner optimization part of
problem (4.13). This will help us to substitute an analytic solution instead of an iterative [11] or
QR decomposition [32] solution for (5.3). Finally, the brute-force part can be substituted with
an appropriate direct search algorithm [17], e.g. pattern search algorithm [16]. It can even be
substituted with an appropriate gradient based algorithm [4, 7], e.g. gradient descent, as we know
the differentiability of P (ω1 , ω2 ).
Algorithm 1 Intrinsic Frequency
(1) Make Df r discrete for a uniform r × r mesh D f r , r ∈ N,



m
l
l
m
ω1 , ω2 ω1 = C, ω2 = C; l, m ∈ {0, 1, . . . , r} .
Df r =
r
r

(2) For all l, m ∈ {0, 1, . . . , r} solve
2
Pn
m
l
minimize
j=1 f (tj ) − S ai , bi , p̄, ω1 , ω2 ; tj
ai ,bi ,p̄

(5.3)

a1 cos ω1 T0 + b1 sin ω1 T0 = a2 cos ω2 T0 + b2 sin ω2 T0 ,
a1
= a2 cos ω2 T + b2 sin ω2 T.

2
P
and store P ω1l , ω2m = nj=1 f (tj ) − S a∗i , b∗i , p̄∗ , ω1l , ω2m ; tj
for minimizers a∗i , b∗i , p̄∗ .
(3) Find the intrinsic frequencies (IFs)

 
(ω1∗ , ω2∗ ) = argmin P ω1l , ω2m .
subject to

l,m

Before moving on, we show the topology of the P (ω1 , ω2 ) function and also its minima locations
in ω1 and ω2 space. These will provide useful insights on where to set the initialization point(s)
of a possible fast IF algorithm. The data description is provided in the next section. In Figures 2
and 3, we have presented two different dog aortic pressure cycles with the IMF extracted by the
means of the brute-force IF Algorithm 1. Figures 2 and 3, top right, show the heat-map plots of


P ω1πT0 , ω2 (Tπ−T0 ) . The complex nature of P (ω1 , ω2 ) can be seen in these figures. We purposefully
plotted P in the dimensionless coordinates

ω1 T0
π

and

ω2 (T −T0 )
π

to show the behavior of this function

with respect to the lattice node locations N defined in (4.20)-(4.22). To have a better view and


understanding of the P (ω1 , ω2 ) topology, a contour of P ω1πT0 , ω2 (Tπ−T0 ) is shown in those figures.
15

The general topology of P



ω1 T0 ω2 (T −T0 )
π ,
π


, for all aortic or carotid pulse waveforms, is similar to

the ones presented in Figures 2 and 3. However, the location of the minimizer is not similar.

Our investigations show that the locations of the minimizers of all P functions construct two different
areas in the dimensionless coordinates

ω1 T0
π

and

ω2 (T −T0 )
.
π

We call these areas as the upper lobe and

lower lobe. The upper lobe is an area, in D, confined above the line
is an area, in D, confined below the line

ω2 (T −T0 )
π

ω2 (T −T0 )
π

= 1. The lower lobe

= 1. This is also the case for human subject data

[25]. This type of topology suggests two critical initial guess areas for any non-brute-force algorithm
solving (4.1): one set of points in the upper lobe, the other in the lower. In the remaining part of
this section, we introduce a fast IF algorithm based on the pattern search method [17].
5.1. Pattern Search IF. The pattern search algorithm (or sometimes called the compass search
algorithm) is explained in detail in [17]. For completeness, we have summarized the pattern search
algorithm in Algorithm 2. The convergence analysis of this method is expressed in [17].
Algorithm 2 Pattern Search [17]
Initialization.
Let f : Rn → R be given.
Let x0 ∈ Rn be the initial guess.
Let △tol > 0 be the tolerance used to test for convergence.
Let △0 > △tol be the initial value of the step length control parameter.
Algorithm. For each iteration k = 1, 2, . . .
Step 1. Let D⊕ be the set of coordinate directions {±ei |i = 1, . . . , n }, where ei is the i th unit
coordinate vector in Rn .
Step 2. If there exists dk ∈ D⊕ such that f (xk + △k dk ) < f (xk ), then do the following:
• Set xk+1 = xk + △k dk .
• Set △k+1 = △k .
Step 3. Otherwise, f (xk + △k dk ) > f (xk ) for all dk ∈ D⊕ , so do the following:
• Set xk+1 = xk .
• Set △k+1 = 12 △k .
• If △k+1 < △tol , then terminate.

The fast IF algorithm, without considering the nodes (4.20), is expressed in Algorithm 3. As
mentioned before, what makes Algorithm 3 fast is embedded in three different objects:
(1) The initial guess set up in the initialization part of the algorithm.
(2) The fast analytic solution at each point iteration defined by (4.17) and (4.16).
(3) The pattern search part which is a substitute for the brute force algorithm.
16

Figure 2, bottom right, shows the results of Algorithm 3. In this figure, when using Algorithm




3, we have used two initial guesses ω1πT0 = 1, ω2 (Tπ−T0 ) = 2 and ω1πT0 = 1, ω2 (Tπ−T0 ) = 0.9 . As

depicted on the figure, the initial guess located in the upper lobe has converged towards the true

minimizer in D. On a PC having 8 threads, Intel® Core™ i7-4700MQ CPU @ 2.40GHz × 8, running
a Matlab implementation of the brute-force Algorithm 1 in parallel takes roughly 85 seconds. On
the other hand, achieving the same minimizer, using a sequential version of the fast Algorithm 3,
takes approximately 0.5 seconds.
The same test was done for another aortic cycle presented in Figure 3. We used the same initial
guesses as before. This time, on the same PC, using the same implementations, the brute-force
Algorithm 1 took roughly 80 seconds and the fast Algorithm 3 took approximately 0.5 seconds.
These two examples show a speed up of almost 160 times. In the next section we present more
about the statistical accuracy of Algorithm 3 and its physiological capabilities.
Algorithm 3 Fast IF
Initialization.
Let f ∈ Rnn+m be a given discrete aortic/carotid signal with
o specified T0 and T − T0 .
ω2 (T −T0 )
ω1 T0
63 .
Let D = (ω1 , ω2 ) 0.5 6 π 6 1.5, 0.5 6
π
SM
Let G = l=1 {(ω1 , ω2 )l ∈ D} be the set of M random initial guesses excluding the nodes (4.20).
Let △ωtol > 0 be the convergence tolerance.
Let △ω0 > △ωtolbe the initial step length.
Let ω̃ k = ω1k , ω2k and
 
  
 

 
2
P ω̃ k = Q ω̃ k , b∗1 ω̃ k , b∗2 ω̃ k ; t + p̄∗ ω̃ k 1 − f
2

for the k th iteration, defined by (4.17), which is solved using (4.16).
Let D⊕ = {±ej |j = 1, 2 }, where ej is the j th unit coordinate vector in R2 .
Algorithm. For each initial guess ω̃ i ∈ G, i = 1, . . . , M , andfor each iteration
ki = 0, 1, . . .

k
k
i
i
Step 1. If there exists dki ∈ D⊕ such that P ω̃ + △ωki dki < P ω̃ , then:
• ω̃ ki +1 = ω̃ ki + △ωki dki .
• △ωki +1 = △ωki .


Step 2. Otherwise, if P ω̃ ki + △ωki dki > P ω̃ ki for all dki ∈ D⊕ , then:

• ω̃ ki +1 = ω̃ ki .
• △ωki +1 = 21 △ωki .
• If △ωki +1 < △ωtol , then terminate and ω̃i∗ = ω̃ ki +1 .
Step 3. The solution isω̃ ∗ = arg min P (ω̃i∗ ).
i∈{1,...,M }

17

6. Real Data Example
The real dog data used in this manuscript is well described in [30]. Since, at the time of the the data
retrieval, the data was downloaded with different sampling rates, we re-sampled all six dog data
at 500 Hz. We used a modified version of the automatic cycle selection introduced in [36] to pick
cycles. Dicrotic notch locations were then found from the picked cycles [19]. We totally extracted
59384 acceptable aortic cycles form those six dogs.
6.1. Statistical Accuracy. To check the statistical accuracy of the fast IF algorithm versus the
brute-force IF algorithm, we compared the results of these two algorithms on the extracted 59384
dog aortic cycles. The brute-force IF algorithm (Algorithm 1) was run over the sample set with a


mesh size (min ω1l − ω1m = min ω2l − ω2m ) of 0.02π. Algorithm 3 was run on the same sample
l6=m

l6=m

set of 59384 aortic cycles with △ωtol = 0.001, and △ω0 = 0.1, comprising a mesh size of

0.1
.
26

The

brute-force algorithm has a larger mesh size due to heavy computational cost of this algorithm.
The maximum average difference between the IFs found by these two algorithms was found to be
less than 0.0475. This difference is smaller than both mesh sizes used for the brute-force and fast
IF algorithms. This shows that, on average, the fast IF algorithm (Algorithm 3) reaches the same
minima as the brute-force algorithm (Algorithm 1).
6.2. Physiological Observations. To evaluate the new fast IF algorithm (Algorithm 3), we applied the algorithm on the measured aortic pressure signal from one dog experiencing various
pharmacological interventions, see Figure 4. During the experiment, the dog was under the following pharmacological influences: infusion of dobutamine (5-20 µg/kg/min), phenylephrine (2-8
µg/kg/min) and nitroglycerin (4 µg/kg/min) during different time intervals.
The third panel, in Figure 4, shows the dosage and duration of each drug in the experiment. In the
first phase of the experiment dobutamine has been injected at a low dosage followed by a fluctuation
in the dosage of injection. The effect of dobutamine on the cardiovascular system is to increase the
strength and force of the heartbeat. Consequently, it forces more blood to circulate throughout
the body. In previous works [25, 26], we hypothesized that ω1 would be a representative of heart
functionality. We also hypothesized that ω1 and ω2 would try to keep a balance during changes.
These hypotheses can be seen during the injection of dobutamine in this figure.
18

Next, phenylephrine has been injected at a low dosage and the dosage is then increased over time.
Phenylephrine is a decongestant which affects the cardiovascular system by shrinking blood vessels.
ω2 shows an almost monotone decrease during the infusion of phenylephrine. This is again in
qualitative accord with what we presented in [26].
Lastly, nitroglycerin has been injected at a constant dosage. Nitroglycerin helps to dilate the blood
vessels. This dilation can be captured with ω2 , as can be seen from the figure. Generally, based
on this figure, IFs are able to capture changes in the dynamics of the system under the effects of
different drugs.

7. Conclusion
In this paper, we provided a mathematical foundation for the IF model [31]. We showed how to
derive an estimation of the IF model (2.1) by considering basic physics principles. More precisely,
we showed that the IF model can be estimated from Navier-Stokes and elasticity equations.
We further analysed the IF model (4.1). This helped to introduce a fast algorithm for the IF method
(Algorithm 3). What made this algorithm fast was embedded in the proper set up of the initial
guesses based on the topology of the problem, fast analytic solution at each point iteration, and
substituting the brute force algorithm with a pattern search method. These changes would convert
an iterative and brute-force method (Algorithm 1) into an algebraic and iterative method (Algorithm
3). The presented fast algorithm, in this article, has a speed up of more than 100 times compared
to the brute-force algorithm provided in [31]. From a statistical perspective, we have also shown
that the algorithm presented in this article complies well with the brute-force implementations of
this method.
We also showed, on a real dataset, that the fast IF Algorithm 3 can depict correlations between its
outputs and infusion of certain drugs. This part of our paper can be subject to further physiological
and clinical investigations in a future work.

19

8. Authors’ Contributions
P.T. conceived of the mathematical and numerical methods of the study, carried out the modeling,
programmed the initial code of the method, and drafted the manuscript; H.K. helped with the
mathematical and numerical derivations, helped draft and revise the manuscript, and conducted
the real data case example; J.K. conducted the brute-force simulations on the real data case example,
and helped draft the manuscript. All authors gave final approval for publication.

9. Acknowledgement
We would like to thank Mr. Sean Brady for constructive discussions and editorial comments.

10. Research Ethics
All experiments and procedures were reviewed and approved by the MSU All-University Committee
on Animal Use and Care [30].

11. Permission to Carry Out Fieldwork
This study did not have fieldwork.

12. Funding
This work was not funded.

13. Figures

Figure 1. Simplified cardiovascular system model schematic
20

160

Pressure

140

120

100

0

0.2

0.4

0.6

Time

Figure 2. Up-Left: A dog aortic pressure cycle (in blue), its dicrotic
notch (black


ω1 T0 ω2 (T −T0 )
dot), and the IMF (in red). Up-Right: heat-map plot of P
,
for
π
π
the cycle in left with
of the solution marked with red dot. Down-Left:
 the location 
ω1 T0 ω2 (T −T0 )
Contour plot of P
. The location of the minimizer of P is shown by
π ,
π


a red dot. Down-Right: Contour plot of P ω1πT0 , ω2 (Tπ−T0 ) and the location of the
minimizer of P tracked by the pattern search. The true optimum point is marked
with a red asterisk. The upper pattern search set (in black) has converged towards
the correct optimum. The lower pattern search set (in magenta) has converged to a
local minima near the node.

21

145

140

135

Pressure

130

125

120

115

110

105
0

0.1

0.2

0.3

0.4

0.5

0.6

Time

Figure 3. Up-Left: A dog aortic pressure cycle (in blue), its dicrotic
notch (black


ω1 T0 ω2 (T −T0 )
dot), and the IMF (in red). Up-Right: heat-map plot of P
,
for
π
π
the cycle in left with
of the solution marked with red dot. Down-Left:
 the location 
ω1 T0 ω2 (T −T0 )
Contour plot of P
. The location of the minimizer of P is shown by
π ,
π


a red dot. Down-Right: Contour plot of P ω1πT0 , ω2 (Tπ−T0 ) and the location of the
minimizer of P tracked by the pattern search. The true optimum point is marked
with a red asterisk. The lower pattern search set (in black) has converged towards
the correct optimum. The upper pattern search set (in magenta) has converged to
a local minima near the node.

22

Figure 4. Drug effects on ω1 and ω2 . First Panel: The measured aortic pressure waveform recorded in time. Second Panel: Dosage of dobutamine (in red),
phenylephrine (in purple), and nitroglycerin (in green) during the aortic pressure
measurement. Third Panel: Changes of ω1 in units of bit per minute (bpm) over
the measurement time. Each drug effect is projected with its corresponding color.
No drug areas are in blue. Fourth Panel: Changes of ω2 in units of bpm over the
measurement time. Each drug effect is projected with its corresponding color. No
drug areas are in blue.

References
[1] Jordi Alastruey, Kim H Parker, and Spencer J Sherwin. Arterial pulse wave haemodynamics. In 11th International
Conference on Pressure Surges, pages 401–442. Virtual PiE Led t/a BHR Group: Lisbon, Portugal, 2012.
[2] Alberto P Avolio, Mark Butlin, and Andrew Walsh. Arterial blood pressure measurement and pulse wave analysistheir role in enhancing cardiovascular assessment. Physiological measurement, 31(1):R1, 2009.
[3] DAVID S Berger, JK Li, and ABRAHAM Noordergraaf. Differential effects of wave reflections and peripheral
resistance on aortic blood pressure: a model-based study. American Journal of Physiology-Heart and Circulatory
Physiology, 266(4):H1626–H1642, 1994.
[4] D.P. Bertsekas. Nonlinear programming. Athena Scientific, 1999.
[5] OF Borisenko and LI Minchenko. Directional derivatives of the maximum function. Cybernetics and Systems
Analysis, 28(2):309–312, 1992.

23

[6] Barry A Borlaug, Carolyn SP Lam, Véronique L Roger, Richard J Rodeheffer, and Margaret M Redfield.
Contractility and ventricular systolic stiffening in hypertensive heart disease: insights into the pathogenesis of
heart failure with preserved ejection fraction. Journal of the American College of Cardiology, 54(5):410–418,
2009.
[7] S. Boyd and L. Vandenberghe. Convex optimization. Cambridge university press, 2004.
[8] Mirko De Melis, Umberto Morbiducci, Ernst R Rietzschel, Marc De Buyzere, Ahmad Qasem, Luc Van Bortel,
Tom Claessens, Franco M Montevecchi, Albert Avolio, and Patrick Segers. Blood pressure waveform analysis by
means of wavelet transform. Medical & biological engineering & computing, 47(2):165–173, 2009.
[9] S.J. Denardo, R. Nandyala, G.L. Freeman, G.L. Pierce, and W.W. Nichols. Pulse wave analysis of the aortic
pressure waveform in severe left ventricular systolic dysfunctionclinical perspective. Circulation: Heart Failure,
3(1):149–156, 2010.
[10] Marc Dewey, Mira Müller, Stephan Eddicks, Dirk Schnapauff, Florian Teige, Wolfgang Rutsch, Adrian C Borges,
and Bernd Hamm. Evaluation of global and regional left ventricular function with 16-slice computed tomography,
biplane cineventriculography, and two-dimensional transthoracic echocardiography: comparison with magnetic
resonance imaging. Journal of the American College of Cardiology, 48(10):2034–2044, 2006.
[11] Nick Gould and Philippe L Toint. Preprocessing for quadratic programming. Mathematical Programming,
100(1):95–132, 2004.
[12] Johannes Greupner, Elke Zimmermann, Andrea Grohmann, Hans-Peter Dübel, Till Althoff, Adrian C Borges,
Wolfgang Rutsch, Peter Schlattmann, Bernd Hamm, and Marc Dewey. Head-to-head comparison of left ventricular function assessment with 64-row computed tomography, biplane left cineventriculography, and both 2-and
3-dimensional transthoracic echocardiography: comparison with magnetic resonance imaging as the reference
standard. Journal of the American College of Cardiology, 59(21):1897–1907, 2012.
[13] W Dallas Hall. Stephen hales: theologian, botanist, physiologist, discoverer of hemodynamics. Clinical cardiology,
10(8):487–489, 1987.
[14] Rainer Hoffmann, Giuseppe Barletta, Stephan von Bardeleben, Jean Louis Vanoverschelde, Jaroslaw Kasprzak,
Christian Greis, and Harald Becher. Analysis of left ventricular volumes and function: a multicenter comparison of cardiac magnetic resonance imaging, cine ventriculography, and unenhanced and contrast-enhanced
two-dimensional and three-dimensional echocardiography. Journal of the American Society of Echocardiography,
27(3):292–301, 2014.
[15] Rainer Hoffmann, Stephan von Bardeleben, Folkert ten Cate, Adrian C Borges, Jaroslaw Kasprzak, Christian
Firschke, Stephane Lafitte, Nidal Al-Saadi, Stefanie Kuntz-Hehner, Marc Engelhardt, et al. Assessment of systolic
left ventricular function: a multi-centre comparison of cineventriculography, cardiac magnetic resonance imaging,
unenhanced and contrast-enhanced echocardiography. European Heart Journal, 26(6):607–616, 2005.
[16] Robert Hooke and Terry A Jeeves. “direct search”solution of numerical and statistical problems. Journal of the
ACM (JACM), 8(2):212–229, 1961.

24

[17] T.G. Kolda, R.M. Lewis, and V. Torczon. Optimization by direct search: New perspectives on some classical
and modern methods. SIAM review, 45(3):385–482, 2003.
[18] Stephane Laurent, John Cockcroft, Luc Van Bortel, Pierre Boutouyrie, Cristina Giannattasio, Daniel Hayoz,
Bruno Pannier, Charalambos Vlachopoulos, Ian Wilkinson, and Harry Struijker-Boudier. Expert consensus document on arterial stiffness: methodological issues and clinical applications. European heart journal, 27(21):2588–
2605, 2006.
[19] Bing Nan Li, Ming Chui Dong, and Mang I Vai. On an automatic delineator for arterial blood pressure waveforms.
Biomedical Signal Processing and Control, 5(1):76–81, 2010.
[20] D. Lloyd-Jones, R.J. Adams, T.M. Brown, M. Carnethon, S. Dai, G. De Simone, T.B. Ferguson, E. Ford,
K. Furie, C. Gillespie, et al. Heart disease and stroke statistics 2010 update. Circulation, 121(7):e46–e215, 2010.
[21] W.R. Milnor. Hemodynamics. Williams & Wilkins, 1989.
[22] Dariush Mozaffarian, Emelia J Benjamin, Alan S Go, Donna K Arnett, Michael J Blaha, Mary Cushman,
Sandeep R Das, Sarah de Ferranti, Jean-pierre Després, Heather J Fullerton, et al. Executive summary: Heart
disease and stroke statistics-2016 update: A report from the american heart association. Circulation, 133(4):447,
2016.
[23] W.W. Nichols and M.F. O’Rourke. McDonald’s blood flow in arteries: theoretical, experimental, and clinical
principles. CRC Press, 2011.
[24] W.W. Nichols, M.F. O’Rourke, C. Hartley, et al. McDonald’s blood flow in arteries: theoretical, experimental
and clinical principles, volume 340809418. Hodder Arnold London, 2005.
[25] Niema M Pahlevan, Derek G Rinderknecht, Peyman Tavallali, Marianne Razavi, Thao T Tran, Michael W Fong,
Robert A Kloner, Marie Csete, and Morteza Gharib. Noninvasive iphone measurement of left ventricular ejection
fraction using intrinsic frequency methodology. Critical care medicine, 2017.
[26] Niema M Pahlevan, Peyman Tavallali, Derek G Rinderknecht, Danny Petrasek, Ray V Matthews, Thomas Y
Hou, and Morteza Gharib. Intrinsic frequency for a systems approach to haemodynamic waveform analysis with
clinical applications. Journal of The Royal Society Interface, 11(98):20140617, 2014.
[27] Kim H Parker. A brief history of arterial wave mechanics. Medical and Biological Engineering and Computing,
47(2):111–118, 2009.
[28] BN Pshenichnyi. Necessary conditions for an extremum. New York, 1971.
[29] Marek W Rajzer, Wiktoria Wojciechowska, Marek Klocek, Ilona Palka, Malgorzata Brzozowska-Kiszka, and
Kalina Kawecka-Jaszcz. Comparison of aortic pulse wave velocity measured by three techniques: Complior,
sphygmocor and arteriograph. Journal of hypertension, 26(10):2001–2007, 2008.
[30] Gokul Swamy, Jacob Kuiper, Madhu SR Gudur, N Bari Olivier, and Ramakrishna Mukkamala. Continuous left
ventricular ejection fraction monitoring by aortic pressure waveform analysis. Annals of biomedical engineering,
37(6):1055, 2009.
[31] Peyman Tavallali, Thomas Y Hou, Derek G Rinderknecht, and Niema M Pahlevan. On the convergence and
accuracy of the cardiovascular intrinsic frequency method. Royal Society Open Science, 2(12):150475, 2015.

25

[32] Lloyd N Trefethen and David Bau III. Numerical linear algebra, volume 50. Siam, 1997.
[33] Yuh-Ying Lin Wang, Wah-Keung Sze, Jian-Guo Bau, Sheng-Hung Wang, Ming-Yie Jan, Tse-Lin Hsu, and WeiKung Wang. The ventricular-arterial coupling system can be analyzed by the eigenwave modes of the whole
arterial system. Applied Physics Letters, 92(15):153901, 2008.
[34] Nicolaas Westerhof, Frederik Bosman, Cornelis J De Vries, and Abraham Noordergraaf. Analog studies of the
human systemic arterial tree. Journal of biomechanics, 2(2):121–143, 1969.
[35] Mair Zamir and M Zamir. The physics of pulsatile flow. Springer, 2000.
[36] W Zong, T Heldt, GB Moody, and RG Mark. An open-source algorithm to detect onset of arterial blood pressure
pulses. In Computers in Cardiology, 2003, pages 259–262. IEEE, 2003.

Appendix A

In this appendix, we show how one can derive (3.2) and (3.1), from the Navier-Stokes and elasticity
equations. Having x, r and θ as the cylindrical coordinate system, with x in the direction of the
aortic length, the momentum and the continuity equations are

(13.1)

ρ

∂u
∂t

µ

(13.2)

(13.3)

(13.4)



∂u
+ u ∂u
∂x + v ∂r +
∂2u
∂x2

+

∂2u
∂r 2

+

ω ∂u
r ∂θ

1 ∂u
r ∂r

+



+

∂P
∂x

1 ∂2u
r 2 ∂θ 2



=
,



∂v
∂v
w ∂v
w2
+ ∂P
+
u
+
v
+
−
ρ ∂v
∂t
∂x
∂r
r ∂θ
r
∂r =
 2

∂2v
1 ∂v
1 ∂2v
2 ∂w
∂ v
v
µ ∂x
,
2 + ∂r 2 + r ∂r − r 2 + r 2 ∂θ 2 − r 2 ∂θ
 1 ∂P
∂w
∂w
w ∂w
vw
ρ ∂w
+ r ∂θ =
∂t + u ∂x + v ∂r + r ∂θ + r
 2

2
1 ∂2w
2 ∂v
w
+
+
µ ∂∂xw2 + ∂∂rw2 + 1r ∂w
−
,
∂r
r2
r 2 ∂θ 2
r 2 ∂θ
∂u ∂v v 1 ∂w
+
+ +
= 0.
∂x ∂r
r r ∂θ

Here, u, v and w are velocity vector components in x, r and θ directions, respectively. We assume
that aorta is a straight and sufficiently long tube with constant circular cross section with the tube
wall following classical elasticity theory dynamics and blood is considered to be an incompressible
Newtonian fluid with the velocity field being axisymmetric. In the absents of any external forces
causing flow rotation, the assumption that the flow field is symmetrical about the longitudinal axis
of the tube is justified. This means w =

∂w
∂θ

=

∂v
∂θ

26

=

∂u
∂θ

=

∂P
∂θ

= 0. Hence, Equations (13.1)-(13.4)

will be simplified as
(13.5)

ρ

(13.6)



ρ



∂u
∂u
∂u
+u
+v
∂t
∂x
∂r

∂v
∂v
∂v
+u
+v
∂t
∂x
∂r





∂P
=µ
+
∂x

∂P
=µ
+
∂r





∂ 2 u ∂ 2 u 1 ∂u
+ 2 +
∂x2
∂r
r ∂r



,

∂ 2 v ∂ 2 v 1 ∂v
v
− 2
+ 2 +
2
∂x
∂r
r ∂r r



,

∂u ∂v v
+
+ = 0.
∂x ∂r
r

(13.7)

Since the radius of the tube a is smaller than the length of the tube, a ≪ h, and also the average
velocity of the blood in aorta is smaller than the speed of wave propagation [35], we have
(13.8)

u

∂u
∂u
≪
,
∂x
∂t

(13.9)

v

∂u
∂u
≪
,
∂r
∂t

(13.10)

u

∂v
∂v
≪
,
∂x
∂t

(13.11)

v

∂v
∂v
≪
,
∂r
∂t

(13.12)

∂2u
∂2u
≪
,
∂x2
∂r 2

(13.13)

∂2v
∂2v
≪
.
∂x2
∂r 2

Using these, Equations (13.5)-(13.7) will reduce to
(13.14)

(13.15)

(13.16)

ρ

∂u ∂P
+
=µ
∂t
∂x

∂v ∂P
+
=µ
ρ
∂t
∂r





∂ 2 u 1 ∂u
+
∂r 2
r ∂r



∂ 2 v 1 ∂v
v
+
− 2
2
∂r
r ∂r r

∂u ∂v v
+
+ =0
∂x ∂r
r
27



It is important to mention that the velocity vector is a function of time and location. In other
words, we have u (x, r, t) and v (x, r, t). Based on the characteristic length of the problem, pressure
P can be assumed to be a function of x and t and not r, i.e. P (x, t) (See Chapter 5 of [35]). Hence
we can set

∂P
∂r

≈ 0 in (13.15).

Momentum Equations: Considering (13.14), one can integrate both sides with respect to the
differential element of the area 2πrdr.
(13.17)

a(x,t)
ˆ

∂u
2πrρ dr +
∂t

a(x,t)
ˆ

0

a(x,t)
ˆ

∂P
dr =
2πr
∂x

2πrµ

0

0



∂ 2 u 1 ∂u
+
∂r 2
r ∂r



dr.

The upper boundary of this integral is the radius a (x, t) of the tube. We know that the flow Q (x, t)
is defined as
(13.18)

Q (x, t) =

a(x,t)
ˆ

2πru (x, r, t) dr.

0

Hence, using Leibniz rule we can find the derivative of the flow with respect to time as

(13.19)

∂Q
(x, t) =
∂t

a(x,t)
ˆ

2πr

∂a
∂u
(x, r, t) dr + 2πa (x, t)
(x, t) u (x, a (x, t) , t) .
∂t
∂t

0

Considering no slip boundary condition u (x, a (x, t) , t) = 0 on the tube wall, Equation (13.19)
reduces to
∂Q
(x, t) =
∂t

(13.20)

a(x,t)
ˆ

2πr

∂u
(x, r, t) dr.
∂t

0

Therefore Equation (13.17) will be simplified to

(13.21)

∂Q
∂P
ρ
+ πa2
= 2πµ
∂t
∂x

a(x,t)
ˆ

r

0



∂ 2 u 1 ∂u
+
∂r 2
r ∂r



dr.

Since we have
(13.22)

a(x,t)
ˆ

r

0



∂ 2 u 1 ∂u
+
∂r 2
r ∂r



dr =

Equation (13.21) would become
28

a(x,t)
ˆ
0


∂ r ∂u
∂r
dr,
∂r

(13.23)

ρ

∂P
∂u a(x,t)
∂Q
+ πa2
= 2πµr
|
.
∂t
∂x
∂r 0

Using separation of variables u (x, r, t) = U (r) ū (x, t), for some function U (r) and ū (x, t), Equation
(13.23) will be simplified to

(13.24)

ρ

∂Q
∂P
+ πa2
= 2πµa (x, t) ū (x, t)
∂t
∂x




dU
|a(x,t) .
dr

From Equation (13.18) and u (x, r, t) = U (r) ū (x, t), we have



(13.25)


ū (x, t) = Q (x, t) 2π


rU (r) dr 

0

Therefore equation (13.24) will be simplified to

(13.26)

ρ

−1

a(x,t)
ˆ

∂P
∂Q
(x, t) + πa2 (x, t)
(x, t) = µa (x, t)
∂t
∂x







dU

|

dr a(x,t)

.

−1

a(x,t)
ˆ
0


rU (r) dr 

Q (x, t) .

From (13.26), we can relabel some terms and introduce the inductance L and resistance R as

L (x, t) =

(13.27)

(13.28)

R (x, t) = −µ

These will convert (13.26) into
(13.29)

−





ρ
πa2 (x, t)



dU

|a(x,t) πa (x, t)
dr

,

a(x,t)
ˆ
0

−1


rU (r) dr 

∂Q
∂P
(x, t) = L (x, t)
(x, t) + R (x, t) Q (x, t) .
∂x
∂t
29

.

Continuity Equation: Again, applying Leibniz rule to equation (13.18) we can find the derivative
of the flow with respect to x

(13.30)

∂Q
(x, t) =
∂x

a(x,t)
ˆ

2πr

∂a
∂u
(x, r, t) dr + 2πa (x, t)
(x, t) u (x, a (x, t) , t) .
∂x
∂x

0

Considering the no slip boundary condition u (x, a(x, t), t) = 0, equation (13.30) reduces to
∂Q
(x, t) =
∂x

(13.31)

a(x,t)
ˆ

2πr

∂u
(x, r, t) dr.
∂x

0

Now, we can rewrite the equation of continuity (13.16) as

(13.32)

a(x,t)
ˆ

∂u
2πr dr +
∂x

a(x,t)
ˆ


∂v
+ v dr = 0.
2π r
∂r

0

0



This equation, using (13.31), will result in
∂Q
+
∂x

(13.33)

a(x,t)
ˆ

2π

∂ (rv)
dr = 0.
∂r

0

Simplifying the latter would show that
∂Q
+ 2πa (x, t) v (a (x, t)) = 0.
∂x

(13.34)
We note that v (a (x, t)) =

∂a
∂t

(x, t).

Hence, having A = πa2 (x, t), we can conclude

∂A
∂t

=

2πa (x, t) v (a (x, t)). Consequently, Equation (13.34) can be written as
(13.35)

∂Q ∂A
+
= 0.
∂x
∂t

Using the chain rule we have
(13.36)

∂A
∂A ∂P
=
.
∂t
∂P ∂t

Considering the wave speed c0 of an incompressible fluid in an elastic tube we have
Matching this with Equations (13.35) and (13.33) would result in
(13.37)

∂Q
πa2 (x, t) ∂P
(x, t) +
(x, t) = 0.
∂x
∂t
ρc20
30

∂A
∂P

=

A
ρc20

[27].

In this equation, we can relabel

πa2 (x,t)
ρc20

−

(13.38)

as the compliance C (x, t). Hence, (13.37) would become

∂Q
∂P
(x, t) = C (x, t)
(x, t) .
∂x
∂t

Wave Equations: As depicted so far, using the approximations in this appendix, and considering
the mentioned assumptions, we can characterize the wave dynamics of the blood flow in aorta using
the hyperbolic equations (13.38) and (13.29), namely
−

(13.39)

(13.40)

−

∂Q
∂P
(x, t) = C (x, t)
(x, t) ,
∂x
∂t

∂P
∂Q
(x, t) = L (x, t)
(x, t) + R (x, t) Q (x, t) .
∂x
∂t

The coefficients C (x, t), L (x, t) and R (x, t), in Equations (13.38) and (13.29), are all positive and
functions of a (x, t). However, since a (x, t) is not changing drastically with respect to x and t, we
can approximate all these coefficients with their corresponding constant mean values C, L and R.
Using this approximation, we can rewrite the wave equations as
−

(13.41)

(13.42)

−

∂P
∂Q
(x, t) = C
(x, t) ,
∂x
∂t

∂P
∂Q
(x, t) = L
(x, t) + RQ (x, t) .
∂x
∂t

31

