Intrinsic Frequency Analysis and Fast Algorithms

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.

capable of non-invasively measuring LVEF by means of an iPhone camera 14 . 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 13 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.

Brief Overview of IF Method
A History of Analyzing Cardiovascular Pulse Waveform. Blood pressure was first measured by Hales in 1735 15 . 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 15 . 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 16 . There has been much research on analyzing the dynamics of blood pressure and flow in arterial systems [17][18][19][20][21] . 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 22 . The formulation of a minimal lumped model of the arterial system was first presented by Westerhof et al. 22 . 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 20 . For example, the impedance method, which is based on Fourier transform, is a common method to analyze the pressure waveform in the frequency domain 17 . As an example, Milnor has shown that the pressure and flow waveforms can be a superposition of several harmonics using the Fourier method 23 . Another method to investigate the pressure wave in the time domain is the wave intensity analysis which is based on wavelet transform 24 . These methods do not necessarily convey a physical understanding of the cardiovascular system.
The IF algorithm presented in 13 is analyzing a pulse waveform through a direct time-frequency signal processing machinery setting, from a quantitative perspective. Although, in previous work 12 , 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.

IF
with a continuity condition at T 0 and periodicity at T. In this formulation, the indicator function is defined as Here, a 1 , b 1 , a 2 and b 2 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 13 : Here, 2 is the L 2 -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 13 . In the next sections, we explore the foundation of the IF algorithm and propose a faster IF algorithm.

Approximate Derivation of the IF Model
As mentioned earlier, in a previous work 12 , 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.
Here we use a simplified model to address our approach. However, for a more general modeling, analysis and estimation of the blood flow and pressure estimation one could see 25,26 .
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 Fig. 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 T 0 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 .
The step by step derivation of these equations is presented in the supplementary material. 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 27 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 P(x, t) from Equations (3.1) and (3.2) by eliminating the flow. Combining Equations (3.1) and (3.2) results in A similar expression could also be found for the flow field Q(x, t). ) , with p 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 ∂ ∂ x t ( , ) p t , in (3.4), by setting n n n n n and some constants α n , β n , ζ n and η n . As a result, the solution of (3.3) can be expressed as 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 T 0 . Hence, for ∈ t T [0, ), 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.
, η n 1 and ω n 1 are found from the boundary and initial conditions at systole. Similarly, constants K 2 , α n 2 , β n 2 , ζ n 2 , η n 2 and ω n 2 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 ω n 1 to ω n 2 . 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 x 0 on aorta, the terms containing the spacial variable x would be fixed. In other words, Equation (3.10) would reduce to Now, considering that the cardiac cycle length would be around 1.5 sec, at most, and taking into account that R is smaller than L 27   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) 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).

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 = T 0 . Hence, the discrete format of (2.2) can be expressed as In this article, (.)′ denotes the transpose operator and the vector 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
General Case (cos ω 1 T 0 cos ω 2 (T − T 0 ) ≠ 1). One can solve the constraints in (4.1) for a 1 and a 2 to obtain where Q(ω 1 , 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 28 . In other words, the optimization problem in (4.12) can be written as 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 13 . To find the exact solution we simplify the objective function as  Since, in this part of the optimization, the values of ω 1 and ω 2 are fixed, we can find the optimal values of b 1 , b 2 , and p by setting the partial derivatives of (4.15) equal to zero. In other words, we set . Doing this, we find the optimal solution for b 1 , b 2 , and p, 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 ( , , ( , ), ( , ); ) ( , ) , (4 17) We note that a property of the function P(ω 1 , ω 2 ) is its differentiability, away from its singularities. In fact, by definition, the function + − p Q 1 f 2 2 is directionally differentiable with respect to all its variables. Hence, using the results in 29,30 , we can deduce that 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 31 .
Degenerate Case (cos ω 1 T 0 cos ω 2 (T − T 0 ) = 1). The solution of (4.6) can be expressed as nodes of a lattice  in ω 1 ω 2 plane. To be more specific, we have If (ω 1 , ω 2 ) ∈ Γ 1 , from (4.5) we have a 1 = −a 2 . On the other hand, if (ω 1 , ω 2 ) ∈ Γ 2 , from (4.5) we have a 1 = a 2 . In both of these cases, we can express (4.3) as where ω ω a b b Q t ( , , , , ; ) Similarly, if (ω 1 , ω 2 ) ∈ Γ 2 , we have In both of the cases, we have Here, 0 1 and 0 2 are zero vectors in  × m 1 and  × n 1 , respectively. It is clear, from (4.26) and (4.27), that ′ ′ = = w w w w 0 . Using (4.23), and a similar approach we employed in (4.15) and (4.16), we find the optimal solution for a 1 , b 1 , b 2 , and p, by Hence, similar to (4.17), for (ω 1 , ω 2 ) ∈ Γ 1 or (ω 1 , ω 2 ) ∈ Γ 2 , 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,

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 13 .
Algorithm 1 has three major steps. In the first step, the (ω 1 ,ω 2 ) domain is made discrete, namely  fr . The second step is a minimization to find P(ω 1 , ω 2 ), see (4.18). The final step is a brute-force search on fr  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  fr , defined in (5.1), we know from our previous work in 14 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 32 or QR decomposition 33 solution for (5.3). Finally, the brute-force part can be substituted with an appropriate direct search algorithm 34 , e.g. pattern search algorithm 35 . It can even be substituted with an appropriate gradient based algorithm 28,31 , e.g. gradient descent, as we know the differentiability of P(ω 1 , ω 2 ). 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 Figs 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 . This is also the case for human subject data 14 . 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 34 .
Pattern Search IF. The pattern search algorithm (or sometimes called the compass search algorithm) is explained in detail in 34 . For completeness, we have summarized the pattern search algorithm in Algorithm 2. The convergence analysis of this method is expressed in 34 .

Real Data Example
Data Description. The real dog data used in this manuscript is well described in 36 . Briefly, six normal adult beagle dogs had undergone the data collection experiment. One dog was involved in a sterile surgical procedure for implanting chronic recording transducers. All the dogs had undergone general anesthesia with propofol and maintained it with inhaled isoflorane. For the aortic pressure waveform measurements, a micro-manometer-tipped catheter was inserted into a femoral artery and guided into the descending thoracic aorta. The transducer outputs were transfered to a personal computer through an A/D conversion system. The measurements were collected over a period of 50-170 minutes. Some pharmacological influences and total intra-vascular volume changes were imposed on the dogs: dobutamine, esmolol, verapamil, phenylephrine, nitroprusside, saline and progressive hemorrhage. For a detailed description, see 36 .
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 37 to pick cycles. Dicrotic notch locations were then found from the picked cycles 38 . We totally extracted 59384 acceptable aortic cycles form those six dogs.

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

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 Fig. 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 second panel, in Fig. 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 third and fourth panels, in Fig. 4, show the trends of ω 1 and ω 2 for both the fast IF algorithm (Algorithm 3 in blue, red, purple and green dots) and the brute-force algorithm (Algorithm 1 in dashed black line). Both algorithms follow the same trends showing the same precision. In other words, if the force-brute IF is used as a control, the fast IF is exactly duplicating the trend.
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 12, 14 , 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.
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 12 .
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.

Conclusion
In this paper, we provided a mathematical foundation for the IF model 13 . 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 13 . 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. Data availability. All data generated or analysed during this study are included in the Supplementary Information files. The datasets generated during and analysed during the current study are also available from the corresponding author on reasonable request.

Research ethics.
All experiments and procedures were reviewed and approved by the MSU All-University Committee on Animal Use and Care 36 .