Intrinsic Frequency Analysis and Fast Algorithms

08/01/2017 ∙ by Peyman Tavallali, et al. ∙ California Institute of Technology 0

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 brute-force 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.



There are no comments yet.


page 21

page 22

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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

[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. 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 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 , for a cardiac period , can be represented as


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

Here, , , and are the envelopes of the IF model fit. and are the Intrinsic Frequencies (IFs) of the waveform. Further, 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 in one period. The latter is done by solving the following optimization problem [31]:


Here, is the -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.

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 during a cardiac cycle . 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 and the pressure along the length of an aorta as follow


The step by step derivation of these equations is presented in Appendix A. Parameters , , and represent inductance, resistance, and compliance of the blood in aorta. Here, , where 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).

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


Taking , with as the the mean pressure, we can write equation (3.3) as


Here, we have used the dot notation to represent the time derivative. We can simplify the term in front of , in (3.4), by setting . The latter has a solution for some constant . This reduces Equation (3.4) into


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




and some constants , , and . As a result, the solution of (3.3) can be expressed as


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


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 , , , , and are found from the boundary and initial conditions at systole. Similarly, constants , , , , and 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 to . 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.

If the pressure is recorded at a specific point on aorta, the terms containing the spacial variable would be fixed. In other words, Equation (3.10) would reduce to




Now, considering that the cardiac cycle length would be around , at most, and taking into account that is smaller than [3], one can use the approximation . Hence, Equation (3.11) would become


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


Now, by relabeling


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 is sampled uniformly. Also, we can simplify (2.2) by the fact that any sinusoid can be assumed to start from time with a compensation coming from a phase shift. In other words, any sinusoid can be expressed as , irrespective of whether the initial time is or . Hence, the discrete format of (2.2) can be expressed as


for as the uniform sampling of the original cycle. Here, by taking


for and , we have the discrete form of as


In this article,

denotes the transpose operator and the vector

. Also,


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


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


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 ()

One can solve the constraints in (4.1) for and to obtain


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


where for




Using Equations (4.9)-(4.11), and dropping the dependencies in notation, simplifies (4.1) into


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


We call the inner optimization in (4.13) as . Solving for 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


Substituting for , we convert (4.14) into


Since, in this part of the optimization, the values of and are fixed, we can find the optimal values of , , and by setting the partial derivatives of (4.15) equal to zero. In other words, we set , , and . Doing this, we find the optimal solution for , , and , by


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


which is


We note that a property of the function is its differentiability, away from its singularities. In fact, by definition, the function is directionally differentiable with respect to all its variables. Hence, using the results in [5, 28], we can deduce that


is directionally differentiable with respect to and . This property can be exploited if one tries to solve (4.18) using a gradient based optimization method [4].

4.2. Degenerate Case ()

The solution of (4.6) can be expressed as nodes of a lattice in plane. To be more specific, we have






If , from (4.5) we have . On the other hand, if , from (4.5) we have . In both of these cases, we can express (4.3) as


where , for . If ,


Similarly, if , we have


In both of the cases, we have




Here, and are zero vectors in and , respectively. It is clear, from (4.26) and (4.27), that . Using (4.23), and a similar approach we employed in (4.15) and (4.16), we find the optimal solution for , , , and , by


for . Hence, similar to (4.17), for or , we only have to solve a minimization on


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, and 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 . 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 domain


is made discrete, namely . The second step is a minimization to find , see (4.18). The final step is a brute-force search on to find the minimum of .

All three steps can be optimized to make the IF algorithm faster. Regarding the domain of optimization , defined in (5.1), we know from our previous work in [25] that the average IF solution, for a physiological pulse waveform recording, is confined to a smaller domain expressed as


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 .

  1. Make discrete for a uniform mesh , ,

  2. For all solve


    and store for minimizers .

  3. Find the intrinsic frequencies (IFs)

Algorithm 1 Intrinsic Frequency

Before moving on, we show the topology of the function and also its minima locations in and 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 . The complex nature of can be seen in these figures. We purposefully plotted in the dimensionless coordinates and to show the behavior of this function with respect to the lattice node locations defined in (4.20)-(4.22). To have a better view and understanding of the topology, a contour of is shown in those figures. The general topology of , 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 functions construct two different areas in the dimensionless coordinates and . We call these areas as the upper lobe and lower lobe. The upper lobe is an area, in , confined above the line . The lower lobe is an area, in , confined below the line . 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].


Let be given.

Let be the initial guess.

Let be the tolerance used to test for convergence.

Let be the initial value of the step length control parameter.

Algorithm. For each iteration

Step 1. Let be the set of coordinate directions , where is the ith unit coordinate vector in .

Step 2. If there exists such that , then do the following:

  • Set .

  • Set .

Step 3. Otherwise, for all , so do the following:

  • Set .

  • Set .

  • If , then terminate.

Algorithm 2 Pattern Search [17]

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.

Figure 2, bottom right, shows the results of Algorithm 3. In this figure, when using Algorithm 3, we have used two initial guesses and