Multichannel Blind Deconvolution Using the Stochastic Calculus for the Estimation of the Central Arterial Pressure

Abstract

A new tool for estimation of both the central arterial pressure and the unknown channel dynamics has been developed. Given two peripheral waveform measurements, this new signal processing algorithm generates two models that represent the distinct branch dynamic behavior associated with the measured signals. The framework for this methodology is based on a Multichannel Blind Deconvolution (MBD) technique that has been reformulated to use Stochastic Calculus (SC). The technique is based on MBD of dynamic system are mathematically analyzed, in order to reconstruct the common unobserved input within an arbitrary scale factor. The convolution process is modeled as a Finite Impulse Response (FIR) filter with unknown coefficients. The source signal is also unknown. Assuming that one of the FIR filter coefficients are time varying, we have been able to get accurate estimation results for the source signal, even though the filter order is unknown. The time varying filter coefficients have been estimated through the SC algorithm, and we have been able to deconvolve the measurements and obtain both the source signal and the convolution path. The positive results demonstrate that the SC approach is superior to conventional methods.

1. Introduction

In this paper, we present a new approach to monitor central arterial pressure using the Multichannel Blind Deconvolution (MBD) [1]. A multichannel blind deconvolution problem can be considered as natural extension or generalization of instantaneous Blind Source Separation (BSS) problem [2]. The MBD is the technique that allows the estimation of both an unknown input and unknown channel dynamics from only channel outputs. It is a powerful tool particularly for the estimation of system dynamics in which a sensor, for measuring the input, is difficult to place. Therefore, the outputs can be explored to find both the channel dynamics and the input. Although one cannot place a sensor to directly measure the input, yet, it may be recovered from the outputs that are simultaneously measured at the multiple branches of the system. The salient feature of MBD is that the channel dynamics and the common input can be identified in real-time based on output data alone. This technique distinguishes itself from other techniques that apply a predetermined transfer function to interpret sensor data. The other techniques cannot account for individual differences nor can they account for dynamic changes in the subject’s physiologic state.

In many situations, one is able to obtain measurements of a modulated signal. The original source is unknown but usually is of great value. To complicate the matter further, the modulating path is usually partially or completely unknown. This situation happens, for example, when measuring the blood pressure at a distance from the heart while one is interested in the blood pressure at the root of the aorta or central arterial pressure (AP). However, in most cases, the modulating path or the filter is slowly varying, that is, the measured signal does not suffer from abrupt changes. Blind deconvolution is the way to estimate the source signal from the received modulated signal given that the modulating path or the filter is unknown IEEE [3].

Central aortic pressure and flow tracings convey important information about the cardiovascular status and the cardiac performance. Traditionally, aortic flow and pressure recording have required invasive catheterization, which may lead to bleeding, infection, and sometimes thromboembolic complications. The challenge in central hemodynamic monitoring is that the central flow and pressure, the input to the arterial system, are not easily measurable.

Over the years, many methods have been suggested for blind deconvolution for estimating central aortic pressure and flow [2–5]. One of the most popular and effective techniques is to assume an FIR model [1], IIR model [4, 6, 7], or generalized FIR models [8, 9] for the modulating channels/paths and to estimate the coefficients of this model. Through the inversion of the FIR filter, we get the original source signals. A major drawback of this model and the others is that if the filter order is incorrect, the estimated source signals will be far from the true values. Hahn et al. [8] suggested the use of a Laguerre series data-compression technique to obtain the central flow from peripheral arterial pressure. The cardiovascular (CV) dynamics is considered to be unknown and time-varying process. However, within a short time window, they assumed that the system can be approximated as a time-invariant linear system. Sugimachi et al.[10] needed the radial arterial pressure and the timing of the central pressure upstroke (for the determination of the wave transmission delay parameter) to reconstruct central aortic pressure. They individualize the transfer function by changing the wave transmission delay parameter. Swamy et al. [1] introduced a technique to estimate the central AP signals as well as estimating other key (CV) physiologic measures from couple of peripheral AP measurements based on a standard Multichannel Blind System Identification (MBSI) method. Xu et al. [11] developed a technique to monitor cardiac output (CO) and left arterial pressure (LAP) by mathematical analysis of the pulmonary artery pressure (PAP) waveform. Porta et al. [12] reviewed cardiovascular variability (CV) regulation systems by pointing out the role of the main rhythms and the various control and functional systems involved.

It is desired that the number of coefficients to reproduce the true impulse response be small, as shown in this paper, because a large number of coefficients is difficult to identify unless the system input signal is rich enough [1, 13] as shown in Figures 6 and 7. In [10], they depend on physical analysis by the collection of measured reference parameters, varying with time, with using the measured transfer function for the estimation of central aortic pressure. In [10], they depend on prior knowledge of the system input.

In this paper, we suggest making some of the FIR filter parameters changing over time. This way, the ambiguity in the filter order is compensated by the time variations of some of the parameters [14]. The time-varying parameters are estimated through the stochastic calculus [15, 16]. The proposed method needs only two distant measurements of the modulated AP. In this paper, we used radial arterial pressure and femoral arterial pressure, and we did not require prior knowledge of the system input to construct central aortic pressure AP. In Section 2, we describe the blind deconvolution problem, conventional solution methods, and the proposed method based on the Ito calculus. In Section 3, we present the results for the reconstruction of single input (the aortic pressure waveform) from two distant measures outputs (peripheral artery pressure waveforms). The same approach could be used in other signals as well. Finally in Section 4, we provide summary and conclusions. The Appendix contains the technical derivation of the proposed method.

2. System Identification and Methodology

The cardiovascular system is topologically analogous to a multichannel dynamic system. Pressure wave emanating from a common source, the heart, is broadcast and transmitted through the many vascular pathways. Therefore, noninvasive circulatory measurements taken at different locations (as shown in Figure 3) can be treated as multichannel data and processed with an MBD algorithm.

Our technique applies a novel MBD method to two peripheral AP waveforms (outputs) in order to reconstruct the central AP waveform (input) within an arbitrary scale factor. The channels relating the common input to each output represent the vascular dynamic properties of different arterial tree paths (see Figure 1) and are assumed to be characterized by finite impulse responses (FIRs). The filters contain many parameters. We estimate the coefficients by the conventional method (in Section 2.1), then we assume that one of the coefficients is varying with time. This way we will be able to compensate for the small number of FIR filter parameters and for the time variation of the channel.

Figure 1 

The $M$ ($>$1) measured and sampled peripheral AP waveforms $[y_i(t),1<i<M]$ are modeled as outputs of $M$ unknown channels driven by the common input $[u(t)]$.

There are many methods for solving problems of estimation time-varying coefficients. But we introduce a method that is based on stochastic calculus (in Section 2.2). Assuming a slow time-varying regression coefficient, we assume that it is evolving according to the Ornstein-Uhlenbeck (OU) process. The unknown parameters of the OU process are estimated by the maximum likelihood method (defined in the Appendix). Finally, through the inversion of the FIR filter, we get the original source signal (Central/Aortic AP) within an arbitrary scale factor.

We will be working in the probability space (Ω,ℱ,𝑃). To simplify the exposure, we will assume that we have only two measurements outputs of a modulated version of the source signal that are given as follows: 𝑦₁(𝑡)=ℎ₁(𝑡)∗𝑢(𝑡)+𝜀₁𝑦(𝑡),₂(𝑡)=ℎ₂(𝑡)∗𝑢(𝑡)+𝜀₂(𝑡), (2.1) where 𝑢(𝑡) is the unknown source signal (central AP), ℎ₁(𝑡) and ℎ₂(𝑡)are unknown filters (hemodynamic response at time 𝑡) or modulating paths, “∗” is the convolution operation, 𝑦₁(𝑡) (femoral AP) and 𝑦₂(𝑡) (radial AP) are the observed measurements, 𝜀₁(𝑡) and 𝜀₂(𝑡) are the measurements noise. The objective is to deconvolve 𝑦₁(𝑡) and 𝑦₂(𝑡) to estimate 𝑢(𝑡). If we convolve 𝑦₁(𝑡) with ℎ₂(𝑡), we will getℎ₂(𝑡)∗𝑦₁(𝑡)=ℎ₂(𝑡)∗(ℎ₁(𝑡)∗𝑢(𝑡))+ℎ₂(𝑡)∗𝜀₁(𝑡), (2.2) since the convolution is a commutative operation, then exchanging ℎ₁(𝑡)and ℎ₂(𝑡) and on the right hand side, we getℎ₂(𝑡)∗𝑦₁(𝑡)=ℎ₁(𝑡)∗(ℎ₂(𝑡)∗𝑢(𝑡))+ℎ₂(𝑡)∗𝜀₁(𝑡)=ℎ₁(𝑡)∗𝑦₂(𝑡)−ℎ₁(𝑡)∗𝜀₂(𝑡)+ℎ₂(𝑡)∗𝜀₁(𝑡). (2.3) Thus,ℎ₂(𝑡)∗𝑦₁(𝑡)=ℎ₁(𝑡)∗𝑦₂(𝑡)−ℎ₁(𝑡)∗𝜀₂(𝑡)+ℎ₂(𝑡)∗𝜀₁(𝑡). (2.4) Note that this equation does not include the input 𝑢(𝑡). It represents the constraints among the channel dynamics or filters and observed output. Substituting a measured time series of output data for 𝑦₁(𝑡) and 𝑦₂(𝑡), the above equation can be solved for the unknown parameters involved in ℎ₁(𝑡) and ℎ₂(𝑡). Once the filters are obtained, we will use their inverses to find an estimate for the source signal. To simplify the exposure further, assume that the modulating filters, that represent the signal paths or channel dynamics, are second-order linear time invariant and have the 𝑍 transforms as follows:ℎ₁(𝑧)=1+𝛽₁𝑧⁻¹+𝛽₂𝑧⁻¹, (2.5) that is,𝑦₁(𝑘)=𝑢(𝑘)+𝛽₁𝑢(𝑘−1)+𝛽₂𝑢(𝑘−2)+𝜀₁(𝑘), (2.6) and in matrix format for 𝑁 data points𝑦₁𝑦(2)₁⋮𝑦(3)₁=𝛽(𝑁−1)²𝛽₁10⋯00𝛽₂𝛽₁1⋯0⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮00𝛽₂𝛽₁1⋮+𝑢(0)𝜀𝑢(1)𝑢(𝑁−1)¹𝜀(2)₁⋮𝜀(3)₁(𝑁−1), (2.7) that is, 𝑌₁=𝐻₁𝑈+𝜀₁, (2.8) where𝑌₁=𝑦₁(2)𝑦₁(3)⋯𝑦₁(𝑁−1)^𝑇,𝑈=𝑢(0)𝑢(1)⋯𝑢(𝑁−1)^𝑇,𝜀₁=𝜀₁(2)𝜀₁(3)⋯𝜀₁(𝑁−1)^𝑇,𝐻₁=𝛽₂𝛽₁10⋯00𝛽₂𝛽₁1⋯0⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮00𝛽₂𝛽₁1, (2.9)(2.10) and “𝑇” stands for transpose.

Similarly,ℎ₂(𝑧)=1+𝛼₁𝑧⁻¹+𝛼₂𝑧⁻², (2.11) that is,𝑦₂(𝑘)=𝑢(𝑘)+𝛼₁𝑢(𝑘−1)+𝛼₂𝑢(𝑘−2)+𝜀₂(𝑘), (2.12) where ℎ₁(𝑧) and ℎ₂(𝑧) are the 𝑍 transforms of the discrete versions of ℎ₁(𝑡) and ℎ₂(𝑡), respectively. Since,ℎ₂(𝑡)∗𝑦₁(𝑡)=ℎ₁(𝑡)∗𝑦₂(𝑡)−ℎ₁(𝑡)∗𝜀₂(𝑡)+ℎ₂(𝑡)∗𝜀₁(𝑡), (2.13) taking the 𝑍 transform of the discrete version of both sides, we getℎ₂(𝑧)𝑦₁(𝑧)=ℎ₁(𝑧)𝑦₂(𝑧)−ℎ₁(𝑧)𝜀₂(𝑧)+ℎ₂(𝑧)𝜀₁(𝑧). (2.14) In the time domain, we get the following equation: 𝑦₁(𝑘)+𝛼₁𝑦₁(𝑘−1)+𝛼₂𝑦₁(𝑘−2)=𝑦₂(𝑘)+𝛽₁𝑦₂(𝑘−1)+𝛽₂𝑦₂(𝑘−2)+[−𝜀₂(𝑘)−𝛽₁𝜀₂(𝑘−1)−𝛽₂𝜀₂(𝑘−2)+𝜀₁(𝑘)+𝛼₁𝜀₁(𝑘−1)−𝛼₂𝜀₁(𝑘−2)]. (2.15)

2.1. A Conventional Method for the Estimation of the System

The familiar scalar regression format as shown in (2.15) is the shape of a regression equation with a correlated error (colored noise). Unless we take this into consideration, the ordinary least square (OLS) method will yield biased estimates for the unknown coefficients 𝛼₁,𝛼₂, 𝛽₁, and 𝛽₂. Rearranging (2.15), we get [𝑦₁(𝑘)−𝑦₂(𝑘)]=−𝛼₁𝑦₁(𝑘−1)−𝛼₂𝑦₁(𝑘−2)+𝛽₁𝑦₂(𝑘−1)+𝛽₂𝑦₂(𝑘−2)+[−𝜀₂(𝑘)−𝛽₁𝜀₂(𝑘−1)−𝛽₂𝜀₂(𝑘−2)+𝜀₁(𝑘)+𝛼₁𝜀₁(𝑘−1)−𝛼₂𝜀₁(𝑘−2)]. (2.16) For the general case where the order of the FIR filters is 𝐼 and 𝐽, we get[𝑦₁(𝑘)−𝑦₂(𝑘)]=−^𝐼_𝑖=1𝛼_𝑖𝑦₁(𝑘−𝑖)+^𝐽_𝑗=1𝛽_𝑗𝑦₂(𝑘−𝑗)+[−^𝐽_𝑗=1𝛽_𝑗𝜀₂(𝑘−𝑗)+^𝐼_𝑖=1𝛼_𝑖𝜀₁(𝑘−𝑖)]+[𝜀₁(𝑘)−𝜀₂(𝑘)]. (2.17) This could be approximated as follows:[𝑦₁(𝑘)−𝑦₂(𝑘)]≈−^𝐼_𝑖=1𝛼_𝑖𝑦₁(𝑘−𝑖)+^𝐽_𝑗=1𝛽_𝑗𝑦₂(𝑘−𝑗)+𝜀(𝑘). (2.18) Since the filter orders are unknown, one could use the corrected Akaike information criterion AIC_c to determine both “𝐼” and “𝐽”. Assume that the error term𝜀(𝑘) is zero mean Gaussian with variance 𝜎², the AIC_c is defined as follows [17]:AIC_c=𝑛(ln𝜎²+1)+2𝑛(𝑝+1)(𝑛−𝑝−2), (2.19) where 𝑛 is the number of observations, 𝑝=𝐼+𝐽 is a number of unknown, and 𝜎² is an estimate of the error variance. We choose the order 𝑝 such that AIC_c is minimized.

Once the coefficients of the FIR filter are estimated, we use inverse filtering to find and estimate for the source signal 𝑈 as follows: 𝑈=𝐻^𝑇₁(𝐻₁𝐻^𝑇₁)⁻¹𝑌₁, (2.20) where 𝑌₁=𝑦₁(2)𝑦₁(3)⋯𝑦₁(𝑁−1)^𝑇,𝑈=𝑢(0)𝑢(1)⋯𝑢(𝑁−1)^𝑇,𝐻₁=𝛽₂𝛽₁0𝛽10⋯0₂𝛽₁𝛽1⋯0⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮00₂𝛽₁1. (2.21)(2.22) The symbol “∧” on top of the variable refers to estimation. For example, 𝛽₂ and 𝛽₁ are the estimates for 𝛽₂ and 𝛽₁, respectively. Due to its simplicity, the above-mentioned method is the one that is commonly used [4].

2.2. Ito Calculus for the Estimation of the Unknown Time-Varying Coefficient

The linear filter assumption is just an approximation to reality. Sometimes the media is nonlinear, time varying, random, or all. The measured signals are usually noisy. The filter order is usually unknown. All the factors suggest that the FIR filter model is an approximation. To compensate for these assumptions, we suggest to make one or some of the unknown coefficients varying with time, that is, (𝛼₁(𝑡),𝛼₂(𝑡),𝛽₁(𝑡),𝛽₂(𝑡)) or in discrete form (𝛼₁(𝑘),𝛼₂(𝑘),𝛽₁(𝑘),𝛽₂(𝑘)). In this paper, however, we restrict ourselves to only one time-varying parameter. Now the problem becomes that of the estimation of the unknown time-varying coefficients. The details of the estimation procedure are given in this section.

We now recast the problem in the format that could be handled by the Ito calculus. Using (2.16), the observed signal, 𝑃_𝑦(𝑘) with (𝑑−1) components, could be modeled as follows:𝑃_𝑦(𝑘)=^𝑑_𝑖=2𝜂_𝑖(𝑘)𝑆_𝑖(𝑘), (2.23) where 𝜂_𝑖(𝑘), 𝑖>1, is the 𝑖th unknown time-varying coefficient, and the random error has been included in each 𝜂_𝑖(𝑘) as follows: 𝑝_𝑦(𝑘)=[𝑦₁(𝑘)−𝑦₂(𝑘)],𝜂₂(𝑘)=−𝛼₁(𝑘),𝜂₃(𝑘)=−𝛼₂𝜂(𝑘),₄(𝑘)=𝛽₁(𝑘),𝜂₅(𝑘)=𝛽₂(𝑘),𝑆₂(𝑘)=𝑦₁(𝑘−1),𝑆₃(𝑘)=𝑦₁𝑆(𝑘−2),₄(𝑘)=𝑦₂(𝑘−1),𝑆₅(𝑘)=𝑦₂(𝑘−2). (2.24) This equation will be used later in another report, if we allow all the coefficients to be time varying [18].

In the proposed approach, it is assumed that the stochastic processes 𝑦₁(𝑘−1), 𝑦₁(𝑘−2), 𝑦₂(𝑘−1), and 𝑦₂(𝑘−2) are independent. If they are correlated, the correlation coefficients will be known. It is also assumed that a stochastic differential equation (SDE) for each process is known. Usually, but not necessary, an Ornstein-Uhlenbeck (OU) process is assumed to describe the evolution of the processes, 𝑦₁(𝑘−1),𝑦₁(𝑘−2),𝑦₂(𝑘−1),and𝑦₂(𝑘−2) [19].

We propose to model the unknown time-varying coefficients as OU processes. The OU models are used when the trend in the time-varying parameter is known or could be guessed. The OU model represents a signal that is bouncing around its trend [20]. In our case, we assume that all the coefficients are constants and only 𝛽₂ has SDE of the OU form𝑑𝛽₂(𝑡)=𝑐₂(𝛽₂(0)−𝛽₂(𝑡))𝑑𝑡+𝑒₂𝑑𝑊₂(𝑡), (2.25) where 𝑐₂, 𝑒₂ are unknown constants to be estimated, 𝑊₂(𝑡) is a Wiener process, and 𝛽₂(0) is the estimated value through the constant coefficient model of the conventional method.

Using the model of Section 2.1, we get[𝑦₁(𝑘)−𝑦₂(𝑘)]≈−𝛼₁𝑦₁(𝑘−1)−𝛼₂𝑦₁(𝑘−2)+𝛽₁𝑦₂(𝑘−1)+𝛽₂𝑦₂(𝑘−2). (2.26) Rearranging to separate the measurements of 𝛽₂(𝑘), we get 𝛽₂(𝑘)≈[𝑦₁(𝑘)−𝑦₂(𝑘)]−[𝛼₁𝑦₁(𝑘−1)−𝛼₂𝑦₁(𝑘−2)+𝛽₁𝑦₂(𝑘−1)]𝑦₂(𝑘−2), (2.27) where 𝛼₁,𝛼₂, and 𝛽₁ are the estimates of 𝛼₁,𝛼₂, and 𝛽₁, respectively. These estimators could be obtained through the least square method or any other method.

𝛽₂(𝑘) of (2.27) is a noisy measurement, not real, because the model order is not known and we have approximated the order by 2 for simplicity. If we use 𝛽₂(𝑘) of (2.27), we get spikes and erroneous estimates of the input AP. As such an estimate for 𝛽₂(𝑘) is needed. We do this by modeling the time-varying coefficient 𝛽₂(𝑘) as an OU process with unknown coefficients. Estimating the coefficients of 𝛽₂(𝑘) will yield an estimate for 𝛽₂(𝑘).

Specifically we have discrete measurements of the stochastic process 𝛽₂(𝑡). We need to estimate the unknown deterministic parameters of this process; mainly 𝑐₂ and 𝑒₂ of (2.25). We use the maximum likelihood method to achieve this objective (see the appendix). Other methods could be used [18] as well.

2.2.1. Estimation of the Diffusion Parameter of (2.25)

For an observation period [0,𝑇], squaring both sides of (2.25) we get [𝑑𝛽₂(𝑡)]²=𝑒₂𝑑𝑡, (2.28) where we use the properties of the Ito calculus, mainly,𝑑𝑡𝑑𝑡=0,𝑑𝑡𝑑𝑊₂(𝑡)=0,𝑑𝑊₂(𝑡)𝑑𝑊₂(𝑡)=𝑑𝑡. (2.29) Thus,𝑒₂̂=1𝑇^𝑇₀[𝑑𝛽₂(𝑡)]²𝑑𝑡. (2.30)

2.2.2. Estimation of the Drift Parameter of (2.25)

As explained in Appendix and following [15, 16], the maximum likelihood estimate of the drift parameter is given as follows:𝑐₂̂∫=−^𝑇0𝛽₂(𝑡)𝑑𝛽₂(𝑡)∫^𝑇₀𝛽²₂(𝑡)𝑑𝑡−𝛽₂(0)∫^𝑇₀𝛽₂(𝑡)𝑑𝑡. (2.31) Thus, an estimate for 𝛽₂(𝑡) is obtained by substituting (2.30) and (2.31) in (2.25). Once the parameters are estimated, we use inverse filtering to find 𝑈 (the aortic pressure waveform) as follows: 𝑈=𝐻^𝑇₁(𝐻₁𝐻^𝑇₁)⁻¹𝑌₁, (2.32) where 𝑌₁=𝑦₁(2)𝑦₁(3)⋯𝑦₁(𝑁−1)^𝑇,𝑈=𝑢(0)𝑢(1)⋯𝑢(𝑁−1)^𝑇,𝐻₁=𝛽₂(𝑡₁)𝛽₁0𝛽10⋯0₂(𝑡₂)𝛽₁𝛽1⋯0⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮00₂(𝑡_𝑁−2)𝛽₁1. (2.33)(2.34) We note from matrix 𝐻₁ that the parameter 𝛽₂(𝑡) has values that are changing across the sample.

2.2.3. Summary of the Proposed Algorithm and Its Assumptions

Our MBD technique is based on a set of assumptions.(1)The common input is obtained at the output of two sensors. The channels relating the common input to reach distinct output in Figure 1 are linear time variant.(2)The system is represented by FIR. The filters ℎ₁(𝑡) and ℎ₂(𝑡) are second order and contain many parameters some of which are time varying.(3)The parameter 𝛽₂(𝑡) is varying with time. Other parameters 𝛼₁, 𝛼₂, and 𝛽₁ are constant.(4)An OU model is assumed for the time varying coefficient 𝛽₂(𝑡).(5)𝛽₂(𝑡) has SDE as shown in (2.25). The parameters, in this equation, are estimated by maximum likelihood method.(6)Equation (2.32) is used with the time-varying parameter 𝛽₂(𝑡) to find an estimate for the input AP.

The aortic pressure waveform input is reconstructed to within an arbitrary scale factor by deconvolving the estimated FIRs from the measured peripheral artery pressure. The technique then calibrates the reconstructed waveform to absolute pressure by using Poiseuille’s law [1]. The reconstructed waveform 𝑈 is calibrated to absolute pressure by scaling it to have the same mean value as the measured peripheral artery pressures as follows. 𝑈_𝑠∑(𝑡)=𝑈(𝑡)^{𝑁−1𝑡=0}((𝑌₁(𝑡)+𝑌₂(𝑡))/2)∑^𝑁−1_𝑡=0𝑈(𝑡),𝑡∈[0,𝑁−1], (2.35) where 𝑈_𝑠(𝑡) is the final absolute (scaled) estimated aortic pressure, that is, 𝑈_𝑠(𝑡)=𝑈(𝑡)×Scalefactor,Scalefactor=AverageFemoralPressure/2+AverageRadialPressure/2.AverageestimatedAorticPressure (2.36)

2.3. The Proposed Algorithm

We now give a brief summary in an algorithmic form to describe the methodology as follows.(1)Insert data in [1] as shown in Figure 2.(2)Describe the form of the stochastic process OU using the formula in (2.25).(3)Estimate the parameters 𝛼₁, 𝛼₂, 𝛽₁, and 𝛽₂ by regression (2.26).(4)𝛽₂(0) is the estimated value through the constant coefficient model of the conventional method (2.18), where 𝛽₂(0) equal 𝛽₂ constant.(5)Measure 𝛽₂ that is being varied with time by (2.27).(6)Estimate the diffusion parameter 𝑒₂ and the drift parameter 𝑐₂ in OU process by (2.30) and (2.31), respectively, and then calculate 𝛽₂(𝑡) (estimated 𝛽₂ varying with time) by (2.25).(7)Calculate the matrix 𝐻₁ in (2.32).(8)Estimate the source signal by (2.32).(9)Calculate the estimated aortic pressure within scale factor by (2.35).

Figure 2 

Segments of measured central arterial pressure (AP), femoral AP, and radial AP waveforms from one swine dataset.

Figure 3 

The aorta and arteries. Solid gold dots indicate pulse points in arteries. These are areas in which the pulse (expansion and contraction of a superficial artery) can be felt.

3. Results

To test the proposed approach, first, we took real data from [1] by inserting the graph that contains the measured data on software CURVESCAN and extracting points (see Figure 2). Second, we simulated a set of 300 data points on computer measured data.

Multichannel blind deconvolution was experimentally evaluated with respect to measured data in which femoral artery pressure (AP), radial AP waveforms, and aortic pressure waveform were simultaneously measured (see Figure 3). Third, we demonstrate the ability of the proposed approach to extract Aorta AP waveform from multichannel. All the calculations in the algorithm were performed under MATLAB 7.2 [20, 21]. Fourth, we evaluated the proposed method by two performance measures.(1)The signal to noise ratio of the estimates (SNRE) was taken as the measure of performance for this evaluation. It is defined as follows: ∑SNRE=10log^𝑘𝑢²(𝑘)∑_𝑖[𝑢(𝑖)−𝑢(𝑖)]², (3.1) where 𝑢(𝑖) is the estimated value of the pressure at instant “𝑖” (see Table 1).(2)The mean absolute percent error (MAPE) was taken as another performance measure for this evaluation (see Table 2). It is defined as follows: 1MAPE=(𝑁^𝑁𝑖=1|𝑢(𝑖)−𝑢(𝑖)|𝑢(𝑖))×100. (3.2)

Table 1 

The first performance measure (SNRE).

Table 2 

The second performance measure (MAPE).

In the first example, it was assumed that the data were noise-free but the order of each of the filters was unknown. In the second example, the data were assumed noisy and the order of each of the filters was unknown. The filter order was estimated using a corrected Akaike information criterion (AIC) [17]. In both examples, by using the proposed method, only one coefficient of the two filters was assumed to be unknown and time varying. The rest of the coefficients were time invariant. It was assumed that this coefficient follows an Ornstein-Uhlenbeck (OU) process with some unknown constant parameters (see (2.25)). The rest of the filter coefficients were assumed unknown deterministic constants.

In both cases, the proposed approach outperformed the conventional method and, therefore, the method was used in the paper of [9]. The (MAPE) value for the estimated central AP data in this paper was 2.615% (see Table 2), and the (MAPE) value in the paper of [9] was 3.2%. The pressure at the root of the aorta (central AP) was successfully estimated.

3.1. Example  1, Noise-Free Measurements

We have two measurements, 𝑦₁(𝑘) (Femoral AP) and 𝑦₂(𝑘) (Radial AP) that are assumed to be represented by the following equations:𝑦₁(𝑘)=𝑢(𝑘)+𝛽₁𝑢(𝑘−1)+𝛽₂𝑢(𝑘−2)+𝛽₃𝑢(𝑘−3)+𝛽₄𝑦𝑢(𝑘−4),₂(𝑘)=𝑢(𝑘)+𝛼₁𝑢(𝑘−1)+𝛼₂𝑢(𝑘−2). (3.3)(3.4)

In the conventional method, the order of the filters was assumed to be just two, that is, 𝛽₃=𝛽₄=0. In the proposed method, we made the same order assumption of two, but we made one of the coefficients time varying. Specifically, we assumed that the SDE of 𝛽₂(𝑡) has the following form: 𝑑𝛽₂(𝑡)=𝑐₂(𝛽₂(0)−𝛽₂(𝑡))𝑑𝑡+𝑒₂𝑑𝑊₂(𝑡), (3.5) where 𝑐₂, 𝑒₂ are unknown constants that were estimated as explained in Sections 2.2.2 and 2.2.3. Remember that in the conventional method [𝑦₁(𝑘)−𝑦₂(𝑘)]=𝛼₁𝑦₁(𝑘−1)−𝛼₂𝑦₁(𝑘−2)+𝛽₁𝑦₂(𝑘−1)+𝛽₂𝑦₂(𝑘−2). (3.6) We estimate the unknown coefficients 𝛼₁, 𝛼₂, 𝛽₁, and 𝛽₂ through regression analysis.

Figure 4 shows the estimated pressure at the root of the aorta AP from the received signal 𝑦₁ (Femoral AP) using the proposed method based on the stochastic calculus (OU) (SNRE_0.0=23.296 db) (see Table 1). SNRE_0.0 is the signal to noise ratio of the estimate with free noise. The mean absolute percent error with free noise (MAPE_0.0) of the estimated central AP waveform was 2.615%.

Figure 4 

The estimated Aorta AP using OU model and measured Aorta or Central AP waveforms.

We compared the estimated central arterial pressure using the conventional method with second order FIR-2 (SNRE_0.0=9.006 db) and the proposed method OU (SNRE_0.0=23.296 db). We note the difference between them in Figure 6, where we estimated the source signal (Aorta AP) from the received signal (Radial AP), see (3.4). While in Figure 7, after applying the conventional method with fourth-order FIR-4, see (3.3), we obtained a better result than by using FIR-2, because the order of filter is increased and the source signal (Aorta AP) is estimated from the received signal (Femoral AP). Scientifically, Estimation central AP from femoral AP is better than radial AP [1, 11]. But still the stochastic calculus (OU) proved to be the best (see Tables 1 and 2). A typical estimate for 𝛽₂(𝑡) compared to the noisy measurements of 𝛽₂(𝑘) is shown in Figure 5.

Figure 5 

A typical estimate for ${\widehat{\beta }}_2(t)$ (2.25) compared to the noisy measurements of ${\beta }_2(k)$ (2.27).

Figure 6 

Estimated AP waveforms from radial AP using the conventional method (FIR-2 model with two orders), compared to estimated AP waveforms from femoral AP using stochastic calculus-based method (OU).

Figure 7 

Estimated Aorta AP waveforms from femoral AP using the conventional method (FIR-4 model with 4 orders), and stochastic calculus-based method (OU).

3.2. Example  2, Noisy Measurements

In this example, white Gaussian noise was added to the original values of the measured pressures. It is assumed that the channels are represented by a second-order FIR filter and fourth-order FIR filter ((3.3) and (3.4)). However, in this case, for conventional method, the order of the filters was estimated by using the corrected AIC. The estimated orders turned out to be 𝐼=4, 𝐽=2. Tables 1 and 2 show the estimated values of central AP using both methods, where SNRE_0.012 is the signal to noise ratio of the estimate with noise variance = 0.012, SNRE_0.024 with noise variance = 0.024, and so forth, and so as MAPE.

Figure 8 shows the measured pressure at the root of the aorta (with added noise) and the estimated pressure using the proposed method that is based on the stochastic calculus. By using the proposed method based on the stochastic calculus, the MAPE_0.012=3.122%. The performances of the FIR filter-based methods and the proposed method are compared in Tables 1 and 2 for different noise levels. The results show the superior performance of the proposed method.

Figure 8 

The estimated Aorta AP using OU model at high noise variance and measured Aorta or Central AP waveforms.

4. Summary and Conclusions

In this report we presented a novel technique to deconvolve the aortic pressure waveform from multiple peripheral artery pressure waveform measurements, using multichannel blind deconvolution. We applied the technique to femoral and radial AP waveforms measured in the swine over 2-minute intervals of a peripheral AP waveform. We assumed that one of the FIR filter coefficients is time varying. Its values were estimated using methods based on the Ito calculus. By this assumption, we were able to compensate for the wrong FIR filter order and the possible time variations of the channels. The results showed superior performance for our proposed approach compared to conventional methods.

In this study, only one unknown time-varying coefficient was assumed to follow the Ornstein-Uhlenbeck stochastic process [20]. Other models could have been used as well. The Ito calculus techniques were used to estimate the coefficients of this Ornstein-Uhlenbeck model. We tested the proposed technique in swine experiments, and our results showed that the MAPE value for the estimated femoral AP data was 2.615%. Our way to reconstructed AP is simple and straightforward. Our method needs only the calculation of pressure wave components in the time domain and does not need calculations in the frequency domain and no need to large computer time. Because of this simplicity, it is quite possible to implement this method in monitoring central pressure AP on line.

In the future, we suggest expanding this method by applying it to real data taken from human cardiovascular simulator. In the presented study, only one parameter was varying with time, while in the future we may use more than just one parameter varying with time. We might study the general case, where no assumptions are imposed on the speed of variations of the time-varying parameters and their numbers. This generalization may improve the accuracy of the estimates. The method of the Malliavin Calculus is proposed to solve this problem [18].

Appendix

Maximum Likelihood Estimate of the Drift Parameters [16, 19, 22]

The maximum likelihood estimate of the drift part of an SDE is discussed in this appendix. Let the SDE of a scalar process 𝑋(𝑡) be given as follows: 𝑑𝑋(𝑡)=𝑏(𝑡)𝑑𝑡+𝜎(𝑋)𝑑𝑊(𝑡). (A.1) The drift, 𝑏(𝑥) could be completely unknown or partially known except for some unknown parameters 𝜃 to be estimated, that is, we have 𝑏(𝑋,𝜃). There are several methods to estimate these unknown parameters. The maximum likelihood method is the most powerful and there are other methods that we will not present. The probability density function 𝑓(𝑋(𝑡)) of the state 𝑋(𝑡) (or the process at time 𝑡) is given as follows: 1𝑓(𝑋(𝑡))=𝐺(𝑏)𝜎²(𝑋(𝑡))exp{2₀^𝑋(𝑡)𝑏(𝑢)𝜎²(𝑢)𝑑𝑢}, (A.2) where 𝐺(𝑏) is the normalization factor such that ∫𝑓(𝑋(𝑡))𝑑𝑥(𝑡)=1 and is given as 𝐺(𝑏)=^+∞_−∞1𝜎²(𝑋(𝑡))exp{2₀^𝑋(𝑡)𝑏(𝑢)𝜎²(𝑢)𝑑𝑢}𝑑𝑋(𝑡). (A.3)

In the sequel, the focus will be on parametric estimation for the vector 𝜃. The likelihood function could be obtained for the process 𝑋(𝑡) when the drift parameters are unknown and the diffusion part is completely known. If the diffusion part of the SDE is partially unknown, we cannot obtain the likelihood function. Let us consider the following three models of the same process 𝑋(𝑡) with three probability measures; namely, 𝑑𝑋(𝑡)=𝑏₁(𝑋)𝑑𝑡+𝜎(𝑋)𝑑𝑊(𝑡),𝑋⁽¹⁾(0), (A.4) with initial conditions 𝑋⁽¹⁾(0) and with probability measure 𝑃₁𝑑𝑋(𝑡)=𝑏₂(𝑋)𝑑𝑡+𝜎(𝑋)𝑑𝑊(𝑡),𝑋⁽²⁾(0), (A.5) with initial conditions 𝑋⁽²⁾(0) and with probability measure 𝑃₂𝑑𝑋(𝑡)=𝜎(𝑋)𝑑𝑊(𝑡),𝑋(0), (A.6) with initial conditions 𝑋(0) and with probability measure 𝑃. Under general conditions, the probabilities measures 𝑃₁, 𝑃₂, and 𝑃 are equivalent and the corresponding Radon-Nikodym derivatives are given as follows: 𝑑𝑃₁=𝑓𝑑𝑃¹(𝑋(0))𝑓(𝑋(0))exp{^𝑇₀𝑏₁(𝑋(𝑡))𝜎²1(𝑋(𝑡))𝑑𝑋(𝑡)−2^𝑇₀𝑏₁²(𝑋(𝑡))𝜎²(𝑋(𝑡))𝑑𝑡},𝑑𝑃₂𝑑𝑃₁=𝑓₂(𝑋(0))𝑓₁(𝑋(0))exp{^𝑇₀𝑏₂(𝑋(𝑡))−𝑏₁(𝑋(𝑡))𝜎²1(𝑋(𝑡))𝑑𝑋(𝑡)−2^𝑇₀𝑏₂²(𝑋(𝑡))−𝑏₁²(𝑋(𝑡))𝜎²(𝑋(𝑡))𝑑𝑡}. (A.7) We will use the equivalence of these measures to find the likelihood function. Let 𝑑𝑋(𝑡)=𝑏(𝑋,𝜃)𝑑𝑡+𝜎(𝑋)𝑑𝑊(𝑡),𝑋(0), (A.8) the measures {𝑃_𝜃,𝜃∈Θ} are equivalent. The likelihood ratio 𝐿(𝜃,𝜃₁,𝑋) becomes 𝐿(𝜃,𝜃₁,𝑋)=𝑑𝑃_𝜃𝑑𝑃_𝜃1=𝐺(𝜃₁)𝐺(𝜃)exp{2₀^𝑋(0)𝑏(𝜃,𝑢)−𝑏(𝜃₁,𝑢)𝜎²(𝑢)𝑑𝑢}×exp{^𝑇₀𝑏(𝜃,𝑋(𝑡))−𝑏(𝜃₁,𝑋(𝑡))𝜎²1(𝑋(𝑡))𝑑𝑋(𝑡)−2^𝑇₀𝑏²(𝜃,𝑋(𝑡))−𝑏²(𝜃₁,𝑋(𝑡))𝜎²(𝑋(𝑡))𝑑𝑡}, (A.9) where 𝜃₁ is a known assumed quantity that has no effect on the maximum likelihood estimates. This is because later on we will maximize with respect to 𝜃. This will result in zero value for the derivatives depending on 𝜃₁ not 𝜃. The likelihood ratio will have another expression that is easier to simulate if we replace the stochastic integral with respect to 𝑋(𝑡) with a Riemann integral, namely, ^𝑇₀𝑏(𝜃,𝑋(𝑡))−𝑏(𝜃₁,𝑋(𝑡))𝜎²(𝑋(𝑡))𝑑𝑋(𝑡)=₀^𝑋(𝑡)𝑏(𝜃,𝑢)−𝑏(𝜃₁,𝑢)𝜎²−1(𝑢)𝑑𝑢2^𝑇₀𝜕[𝜕𝑋(𝑡)𝑏(𝜃,𝑋(𝑡))−𝑏(𝜃₁,𝑋(𝑡))𝜎²(𝑋(𝑡))]𝜎²(𝑋(𝑡))𝑑𝑡. (A.10) The maximum likelihood estimate 𝜃 is defined as the solution of the equation 𝜃=argsup_𝜃𝐿(𝜃,𝜃₁,𝑋), (A.11) or 𝜃=argsup_𝜃log𝐿(𝜃,𝜃₁,𝑋). (A.12) In the case that the likelihood function has a continuous derivative with respect to 𝜃∈Θ, the maximum likelihood estimate 𝜃 is obtained by solving the log likelihood function as follows: 𝜕log𝐿(𝜃,𝜃₁,𝑋)𝜕𝜃|_𝜃=𝜃=0, (A.13) that is, −𝜕𝐺/𝜕𝜃𝜕𝐺(𝜃)+2₀^𝑋(0)𝜕𝑏(𝑢,𝜃)/𝜕𝜃𝜎²(𝑢)𝑑𝑢+^𝑇₀𝜕𝑏(𝑋(𝑡),𝜃)/𝜕𝜃𝜎²(𝑋(𝑡))[𝑑𝑋(𝑡)−𝑏(𝑋(𝑡),𝜃)𝑑𝑡]=0. (A.14) This formula has an Ito integral with respect to 𝑑𝑋(𝑡). Instead, there is another easier formula 0=−𝜕𝐺/𝜕𝜃𝜕𝐺(𝜃)+2₀^𝑋(0)𝜕𝑏(𝑢,𝜃)/𝜕𝜃𝜎²(𝑢)𝑑𝑢+₀^𝑋(𝑇)𝜕𝑏(𝑢,𝜃)/𝜕𝜃𝜎²−(𝑢)𝑑𝑢_𝑇0[𝑏(𝑋(𝑡),𝜃)𝜕𝑏(𝑋(𝑡),𝜃)/𝜕𝜃𝜎²+𝜕(𝑋(𝑡))(𝜕𝑋(𝑡)𝜕𝑏(𝑋(𝑡),𝜃)/𝜕𝜃2𝜎²(𝑋(𝑡)))𝜎²(𝑋(𝑡))]𝑑𝑡. (A.15) In the case that 𝑋(0) is independent of 𝜃, we have a reduced expression as follows: ^𝑇₀𝜕𝑏(𝑋(𝑡),𝜃)𝜕𝜃𝜎²(𝑋(𝑡))[𝑑𝑋(𝑡)−𝑏(𝑋(𝑡),𝜃)𝑑𝑡]=0 (A.16) or ₀^𝑋(𝑇)𝜕𝑏(𝑢,𝜃)𝜕𝜃𝜎²−(𝑢)𝑑𝑢_𝑇0[𝑏(𝑋(𝑡),𝜃)𝜕𝑏(𝑋(𝑡),𝜃)/𝜕𝜃𝜎²+𝜕(𝑋(𝑡))(𝜕𝑋(𝑡)𝜕𝑏(𝑋(𝑡),𝜃)/𝜕𝜃2𝜎²(𝑋(𝑡)))𝜎²(𝑋(𝑡))]𝑑𝑡=0. (A.17) It is this expression that we used to find the estimates of the drift parameters of (Section 2.2.2). The previous analysis is general case. In this paper, assume that we have Ornstein_Uhlenbeck (OU) process given by the SDE equation 𝑑𝛽₂(𝑡)=𝑐₂(𝛽₂(0))−𝛽₂(𝑡)𝑑𝑡+𝑒₂𝑑𝑊₂(𝑡). (A.18) Comparing to (A.8), We have 𝑏(𝑋(𝑡),𝜃)=𝑏(𝛽₂(𝑡),𝑐₂)=𝑐₂(𝛽₂(0)−𝛽₂(𝑡)),𝜕𝑏(𝑋(𝑡),𝜃)/𝜕𝜃𝜃=𝜕𝑏(𝛽₂(𝑡),𝑐₂)/𝜕𝑐₂𝑐₂=−𝛽₂(𝑡). (A.19) Substituting in (A.16), we get ^𝑇₀𝜕𝑏(𝛽₂(𝑡),𝑐₂)/𝜕𝑐₂𝑐₂[𝑑𝛽₂(𝑡)−𝑐₂(𝛽(0))−𝛽₂(𝑡)𝑑𝑡]=0,_𝑇0−𝛽₂(𝑡)𝑒²𝑑𝛽₂(𝑡)+𝑐₂^𝑇₀𝛽₂(𝑡)[𝛽₂(0)−𝛽₂(𝑡)]𝑑𝑡𝑒²=0,_𝑇0−𝛽₂(𝑡)𝑒²𝑑𝛽₂(𝑡)=−𝑐₂^𝑇₀𝛽₂(𝑡)[𝛽₂(0)−𝛽₂(𝑡)]𝑑𝑡𝑒². (A.20) This is reduced to 𝑐₂̂∫=−^𝑇0𝛽₂(𝑡)𝑑𝛽₂(𝑡)∫^𝑇₀𝛽²₂(𝑡)𝑑𝑡−𝛽₂(0)∫^𝑇₀𝛽₂(𝑡)𝑑𝑡. (A.21)