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

Authors: Ahmed S. Abutaleb; M. El-Sayed Waheed; Nermeen M. Elhamy
Journal: Mathematical Problems in Engineering
 (2010)

---

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


---

## Body

## 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 Figure3) 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 Figure1) 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
TheM (>1) measured and sampled peripheral AP waveforms [yi(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 Section2.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(Ω,ℱ,P). 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)y1(t)=h1(t)*u(t)+ε1(t),y2(t)=h2(t)*u(t)+ε2(t),
where u(t) is the unknown source signal (central AP), h1(t) and h2(t)are unknown filters (hemodynamic response at time t) or modulating paths, “*” is the convolution operation, y1(t) (femoral AP) and y2(t) (radial AP) are the observed measurements, ε1(t) and ε2(t) are the measurements noise. The objective is to deconvolve y1(t) and y2(t) to estimate u(t).  If we convolve y1(t) with h2(t), we will get(2.2)h2(t)*y1(t)=h2(t)*(h1(t)*u(t))+h2(t)*ε1(t),
since the convolution is a commutative operation, then exchanging h1(t)andh2(t) and on the right hand side, we get(2.3)h2(t)*y1(t)=h1(t)*(h2(t)*u(t))+h2(t)*ε1(t)=h1(t)*y2(t)-h1(t)*ε2(t)+h2(t)*ε1(t).
Thus,(2.4)h2(t)*y1(t)=h1(t)*y2(t)-h1(t)*ε2(t)+h2(t)*ε1(t).
Note that this equation does not include the input u(t). It represents the constraints among the channel dynamics or filters and observed output. Substituting a measured time series of output data for y1(t) and y2(t), the above equation can be solved for the unknown parameters involved in h1(t) and h2(t).  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 Z transforms as follows:(2.5)h1(z)=1+β1z-1+β2z-1,
that is,(2.6)y1(k)=u(k)+β1u(k-1)+β2u(k-2)+ε1(k),
and in matrix format for N data points(2.7)[y1(2)y1(3)⋮y1(N-1)]=[β2β110⋯00β2β11⋯0⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮00β2β11][u(0)u(1)⋮u(N-1)]+[ε1(2)ε1(3)⋮ε1(N-1)],
that is,
(2.8)Y1̲=H1U̲+ε1̲,
where(2.9)Y1=[y1(2)y1(3)⋯y1(N-1)]T,U̲=[u(0)u(1)⋯u(N-1)]T,ε1̲=[ε1(2)ε1(3)⋯ε1(N-1)]T,(2.10)H1=[β2β110⋯00β2β11⋯0⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮00β2β11],
and “T” stands for transpose.Similarly,(2.11)h2(z)=1+α1z-1+α2z-2,
that is,(2.12)y2(k)=u(k)+α1u(k-1)+α2u(k-2)+ε2(k),
where h1(z) and h2(z) are the Z transforms of the discrete versions of h1(t) and h2(t), respectively. Since,(2.13)h2(t)*y1(t)=h1(t)*y2(t)-h1(t)*ε2(t)+h2(t)*ε1(t),
taking the Z transform of the discrete version of both sides, we get(2.14)h2(z)y1(z)=h1(z)y2(z)-h1(z)ε2(z)+h2(z)ε1(z).
In the time domain, we get the following equation:
(2.15)y1(k)+α1y1(k-1)+α2y1(k-2)=y2(k)+β1y2(k-1)+β2y2(k-2)+[-ε2(k)-β1ε2(k-1)-β2ε2(k-2)+ε1(k)+α1ε1(k-1)-α2ε1(k-2)].


### 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 α1,α2, β1, and β2. Rearranging (2.15), we get
(2.16)[y1(k)-y2(k)]=-α1y1(k-1)-α2y1(k-2)+β1y2(k-1)+β2y2(k-2)+[-ε2(k)-β1ε2(k-1)-β2ε2(k-2)+ε1(k)+α1ε1(k-1)-α2ε1(k-2)].
For the general case where the order of the FIR filters is I and J, we get(2.17)[y1(k)-y2(k)]=-∑i=1Iαiy1(k-i)+∑j=1Jβjy2(k-j)+[-∑j=1Jβjε2(k-j)+∑i=1Iαiε1(k-i)]+[ε1(k)-ε2(k)].
This could be approximated as follows:(2.18)[y1(k)-y2(k)]≈-∑i=1Iαiy1(k-i)+∑j=1Jβjy2(k-j)+ε(k).
Since the filter orders are unknown, one could use the corrected Akaike information criterion AICc to determine both “I” and “J”. Assume that the error termε(k) is zero mean Gaussian with variance σ2, the AICc is defined as follows [17]:(2.19)AICc=n(lnσ2+1)+2n(p+1)(n-p-2),
where n is the number of observations, p=I+J is a number of unknown, and σ2 is an estimate of  the error variance. We choose the order p such that AICc  is minimized.Once the coefficients of the FIR filter are estimated, we use inverse filtering to find and estimate for the source signalU as follows: (2.20)Û̲=Ĥ1T(Ĥ1Ĥ1T)-1Y1̲,
where
(2.21)Y1̲=[y1(2)y1(3)⋯y1(N-1)]T,Û̲=[û(0)û(1)⋯û(N-1)]T,(2.22)Ĥ1=[β̂2β̂110⋯00β̂2β̂11⋯0⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮00β̂2β̂11].
The symbol “∧” on top of the variable refers to estimation. For example, β̂2  and β̂1  are the estimates for  β2 and β1, 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,(α1(t),α2(t),β1(t),β2(t)) or in discrete form (α1(k),α2(k),β1(k),β2(k)). 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, Py(k) with (d-1) components, could be modeled as follows:(2.23)Py(k)=∑i=2dηi(k)Si(k),
where ηi(k), i>1, is the ith unknown time-varying coefficient, and the random error has been included in each ηi(k) as follows:
(2.24)py(k)=[y1(k)-y2(k)],η2(k)=-α1(k),η3(k)=-α2(k),η4(k)=β1(k),η5(k)=β2(k),S2(k)=y1(k-1),S3(k)=y1(k-2),S4(k)=y2(k-1),S5(k)=y2(k-2).
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 processesy1(k-1), y1(k-2), y2(k-1), and y2(k-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, y1(k-1),y1(k-2),y2(k-1),andy2(k-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 β̂2 has SDE of the OU form(2.25)dβ2(t)=c2(β2(0)-β2(t))dt+e2dW2(t),
where c2, e2 are unknown constants to be estimated, W2(t) is a Wiener process, and β2(0) is the estimated value through the constant coefficient model of the conventional method.Using the model of Section2.1, we get(2.26)[y1(k)-y2(k)]≈-α1y1(k-1)-α2y1(k-2)+β1y2(k-1)+β2y2(k-2).
Rearranging to separate the measurements of β2(k), we get
(2.27)β2(k)≈[y1(k)-y2(k)]-[α̂1y1(k-1)-α̂2y1(k-2)+β̂1y2(k-1)]y2(k-2),
where α̂1,α̂2, and β̂1 are the estimates of α1,α2, and β1, respectively. These estimators could be obtained through the least square method or any other method.β2(k) 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 β2(k)  of  (2.27), we get spikes and erroneous estimates of the input AP. As such an estimate for β2(k) is needed. We do this by modeling the time-varying coefficient β2(k) as an OU process with unknown coefficients. Estimating the coefficients of β2(k) will yield an estimate for β2(k).Specifically we have discrete measurements of the stochastic processβ2(t). We need to estimate the unknown deterministic parameters of this process; mainly c2 and e2 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,T], squaring both sides of (2.25) we get
(2.28)[dβ2(t)]2=e2dt,
where we use the properties of the Ito calculus, mainly,(2.29)dtdt=0,dtdW2(t)=0,dW2(t)dW2(t)=dt.
Thus,(2.30)ê2=1T∫0T[dβ2(t)]2dt.


#### 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:(2.31)ĉ2=-∫0Tβ2(t)dβ2(t)∫0Tβ22(t)dt-β2(0)∫0Tβ2(t)dt.
Thus, an estimate for β2(t) is obtained by substituting (2.30) and (2.31) in (2.25). Once the parameters are estimated, we use inverse filtering to find U (the aortic pressure waveform) as follows: (2.32)Û̲=Ĥ1T(Ĥ1Ĥ1T)-1Y1̲,
where
(2.33)Y1̲=[y1(2)y1(3)⋯y1(N-1)]T,Û̲=[û(0)û(1)⋯û(N-1)]T,(2.34)Ĥ1=[β̂2(t1)β̂110⋯00β̂2(t2)β̂11⋯0⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮00β̂2(tN-2)β̂11].
We note from matrix Ĥ1 that the parameter β̂2(t) 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 Figure1 are linear time variant.(2)
The system is represented by FIR. The filtersh1(t) and h2(t) are second order and contain many parameters some of which are time varying.(3)
The parameterβ2(t) is varying with time. Other parameters α1, α2, and β1 are constant.(4)
An OU model is assumed for the time varying coefficientβ2(t).(5)
β̂2(t) 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 β̂2(t) 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 Û̲ is calibrated to absolute pressure by scaling it to have the same mean value as the measured peripheral artery pressures as follows. 										(2.35)Ûs(t)=Û(t)∑t=0N-1((Y1(t)+Y2(t))/2)∑t=0N-1Û(t),t∈[0,N-1],
where Ûs(t) is the final absolute (scaled) estimated aortic pressure, that is,
(2.36)Ûs(t)=Û(t)×Scale factor,Scale factor=Average Femoral Pressure/2+Average  Radial Pressure/2Average estimated Aortic Pressure.


### 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α̂1, α̂2, β̂1, and β̂2 by regression (2.26).(4)
β2(0) is the estimated value through the constant coefficient model of the conventional method (2.18), where β2(0) equal β̂2 constant.(5)
Measureβ2 that is being varied with time by (2.27).(6)
Estimate the diffusion parametere2 and the drift parameter c2 in OU process by (2.30) and (2.31), respectively, and then calculate β̂2(t) (estimated β2 varying with time) by (2.25).(7)
Calculate the matrixĤ1 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.


## 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 α1,α2, β1, and β2. Rearranging (2.15), we get
(2.16)[y1(k)-y2(k)]=-α1y1(k-1)-α2y1(k-2)+β1y2(k-1)+β2y2(k-2)+[-ε2(k)-β1ε2(k-1)-β2ε2(k-2)+ε1(k)+α1ε1(k-1)-α2ε1(k-2)].
For the general case where the order of the FIR filters is I and J, we get(2.17)[y1(k)-y2(k)]=-∑i=1Iαiy1(k-i)+∑j=1Jβjy2(k-j)+[-∑j=1Jβjε2(k-j)+∑i=1Iαiε1(k-i)]+[ε1(k)-ε2(k)].
This could be approximated as follows:(2.18)[y1(k)-y2(k)]≈-∑i=1Iαiy1(k-i)+∑j=1Jβjy2(k-j)+ε(k).
Since the filter orders are unknown, one could use the corrected Akaike information criterion AICc to determine both “I” and “J”. Assume that the error termε(k) is zero mean Gaussian with variance σ2, the AICc is defined as follows [17]:(2.19)AICc=n(lnσ2+1)+2n(p+1)(n-p-2),
where n is the number of observations, p=I+J is a number of unknown, and σ2 is an estimate of  the error variance. We choose the order p such that AICc  is minimized.Once the coefficients of the FIR filter are estimated, we use inverse filtering to find and estimate for the source signalU as follows: (2.20)Û̲=Ĥ1T(Ĥ1Ĥ1T)-1Y1̲,
where
(2.21)Y1̲=[y1(2)y1(3)⋯y1(N-1)]T,Û̲=[û(0)û(1)⋯û(N-1)]T,(2.22)Ĥ1=[β̂2β̂110⋯00β̂2β̂11⋯0⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮00β̂2β̂11].
The symbol “∧” on top of the variable refers to estimation. For example, β̂2  and β̂1  are the estimates for  β2 and β1, 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,(α1(t),α2(t),β1(t),β2(t)) or in discrete form (α1(k),α2(k),β1(k),β2(k)). 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, Py(k) with (d-1) components, could be modeled as follows:(2.23)Py(k)=∑i=2dηi(k)Si(k),
where ηi(k), i>1, is the ith unknown time-varying coefficient, and the random error has been included in each ηi(k) as follows:
(2.24)py(k)=[y1(k)-y2(k)],η2(k)=-α1(k),η3(k)=-α2(k),η4(k)=β1(k),η5(k)=β2(k),S2(k)=y1(k-1),S3(k)=y1(k-2),S4(k)=y2(k-1),S5(k)=y2(k-2).
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 processesy1(k-1), y1(k-2), y2(k-1), and y2(k-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, y1(k-1),y1(k-2),y2(k-1),andy2(k-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 β̂2 has SDE of the OU form(2.25)dβ2(t)=c2(β2(0)-β2(t))dt+e2dW2(t),
where c2, e2 are unknown constants to be estimated, W2(t) is a Wiener process, and β2(0) is the estimated value through the constant coefficient model of the conventional method.Using the model of Section2.1, we get(2.26)[y1(k)-y2(k)]≈-α1y1(k-1)-α2y1(k-2)+β1y2(k-1)+β2y2(k-2).
Rearranging to separate the measurements of β2(k), we get
(2.27)β2(k)≈[y1(k)-y2(k)]-[α̂1y1(k-1)-α̂2y1(k-2)+β̂1y2(k-1)]y2(k-2),
where α̂1,α̂2, and β̂1 are the estimates of α1,α2, and β1, respectively. These estimators could be obtained through the least square method or any other method.β2(k) 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 β2(k)  of  (2.27), we get spikes and erroneous estimates of the input AP. As such an estimate for β2(k) is needed. We do this by modeling the time-varying coefficient β2(k) as an OU process with unknown coefficients. Estimating the coefficients of β2(k) will yield an estimate for β2(k).Specifically we have discrete measurements of the stochastic processβ2(t). We need to estimate the unknown deterministic parameters of this process; mainly c2 and e2 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,T], squaring both sides of (2.25) we get
(2.28)[dβ2(t)]2=e2dt,
where we use the properties of the Ito calculus, mainly,(2.29)dtdt=0,dtdW2(t)=0,dW2(t)dW2(t)=dt.
Thus,(2.30)ê2=1T∫0T[dβ2(t)]2dt.


### 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:(2.31)ĉ2=-∫0Tβ2(t)dβ2(t)∫0Tβ22(t)dt-β2(0)∫0Tβ2(t)dt.
Thus, an estimate for β2(t) is obtained by substituting (2.30) and (2.31) in (2.25). Once the parameters are estimated, we use inverse filtering to find U (the aortic pressure waveform) as follows: (2.32)Û̲=Ĥ1T(Ĥ1Ĥ1T)-1Y1̲,
where
(2.33)Y1̲=[y1(2)y1(3)⋯y1(N-1)]T,Û̲=[û(0)û(1)⋯û(N-1)]T,(2.34)Ĥ1=[β̂2(t1)β̂110⋯00β̂2(t2)β̂11⋯0⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮00β̂2(tN-2)β̂11].
We note from matrix Ĥ1 that the parameter β̂2(t) 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 Figure1 are linear time variant.(2)
The system is represented by FIR. The filtersh1(t) and h2(t) are second order and contain many parameters some of which are time varying.(3)
The parameterβ2(t) is varying with time. Other parameters α1, α2, and β1 are constant.(4)
An OU model is assumed for the time varying coefficientβ2(t).(5)
β̂2(t) 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 β̂2(t) 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 Û̲ is calibrated to absolute pressure by scaling it to have the same mean value as the measured peripheral artery pressures as follows. 										(2.35)Ûs(t)=Û(t)∑t=0N-1((Y1(t)+Y2(t))/2)∑t=0N-1Û(t),t∈[0,N-1],
where Ûs(t) is the final absolute (scaled) estimated aortic pressure, that is,
(2.36)Ûs(t)=Û(t)×Scale factor,Scale factor=Average Femoral Pressure/2+Average  Radial Pressure/2Average estimated Aortic Pressure.


## 2.2.1. Estimation of the Diffusion Parameter of (2.25)

For an observation period[0,T], squaring both sides of (2.25) we get
(2.28)[dβ2(t)]2=e2dt,
where we use the properties of the Ito calculus, mainly,(2.29)dtdt=0,dtdW2(t)=0,dW2(t)dW2(t)=dt.
Thus,(2.30)ê2=1T∫0T[dβ2(t)]2dt.


## 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:(2.31)ĉ2=-∫0Tβ2(t)dβ2(t)∫0Tβ22(t)dt-β2(0)∫0Tβ2(t)dt.
Thus, an estimate for β2(t) is obtained by substituting (2.30) and (2.31) in (2.25). Once the parameters are estimated, we use inverse filtering to find U (the aortic pressure waveform) as follows: (2.32)Û̲=Ĥ1T(Ĥ1Ĥ1T)-1Y1̲,
where
(2.33)Y1̲=[y1(2)y1(3)⋯y1(N-1)]T,Û̲=[û(0)û(1)⋯û(N-1)]T,(2.34)Ĥ1=[β̂2(t1)β̂110⋯00β̂2(t2)β̂11⋯0⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮00β̂2(tN-2)β̂11].
We note from matrix Ĥ1 that the parameter β̂2(t) 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 Figure1 are linear time variant.(2)
The system is represented by FIR. The filtersh1(t) and h2(t) are second order and contain many parameters some of which are time varying.(3)
The parameterβ2(t) is varying with time. Other parameters α1, α2, and β1 are constant.(4)
An OU model is assumed for the time varying coefficientβ2(t).(5)
β̂2(t) 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 β̂2(t) 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 Û̲ is calibrated to absolute pressure by scaling it to have the same mean value as the measured peripheral artery pressures as follows. 										(2.35)Ûs(t)=Û(t)∑t=0N-1((Y1(t)+Y2(t))/2)∑t=0N-1Û(t),t∈[0,N-1],
where Ûs(t) is the final absolute (scaled) estimated aortic pressure, that is,
(2.36)Ûs(t)=Û(t)×Scale factor,Scale factor=Average Femoral Pressure/2+Average  Radial Pressure/2Average estimated Aortic Pressure.


## 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α̂1, α̂2, β̂1, and β̂2 by regression (2.26).(4)
β2(0) is the estimated value through the constant coefficient model of the conventional method (2.18), where β2(0) equal β̂2 constant.(5)
Measureβ2 that is being varied with time by (2.27).(6)
Estimate the diffusion parametere2 and the drift parameter c2 in OU process by (2.30) and (2.31), respectively, and then calculate β̂2(t) (estimated β2 varying with time) by (2.25).(7)
Calculate the matrixĤ1 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 Figure3). 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:
(3.1)SNRE=10log∑ku2(k)∑i[u(i)-û(i)]2,
where û(i) is the estimated value of the pressure at instant “i” (see Table 1).(2)
The mean absolute percent error (MAPE) was taken as another performance measure for this evaluation (see Table2). It is defined as follows:
(3.2)MAPE=(1N∑i=1N|û(i)-u(i)|u(i))×100.Table 1
The first performance measure (SNRE).
MethodsSNRE0.0SNRE0.012SNRE0.024SNRE0.037SNRE0.049SNRE0.0615FIR-2 order9.006 db7.061 db6.730 db6.380 db6.073 db5.995 dbFIR-4 order12.349 db11.248 db11.076 db10.974 db10.582 db10.108 dbStochastic Calculus23.296 db22.049 db21.047 db20.727 db19.400 db18.199 dbThe multichannel blind deconvolution (MBD) technique was applied by using the conventional method (FIR-2 model with two orders and FIR-4 model with 4 orders) and the proposed method (OU model or stochastic calculus). We compared between these methods by using performance measure SNRE (3.1) at variant noise. SNRE0.0 is the signal to noise ratio of the estimate of Aortic AP with free noise. SNRE0.012 is the signal to noise ratio of the estimate with noise variance = 0.012, SNRE0.024 with noise variance = 0.024, and so forth.Table 2
The second performance measure (MAPE).
MethodsMAPE0.0MAPE0.012MAPE0.024MAPE0.037MAPE0.049MAPE0.0615FIR-2 order21.969%28.853%31.050%32.274%33.084%34.423%FIR-4 order12.941%14.061%17.120%17.966%18.090%20.801%Stochastic Calculus2.615%3.122%3.218%3.670%4.134%4.754%The performance measure for the estimated femoral AP dataMAPE (3.2) was used to evaluate the cardiovascular system at variant noise for three methods (FIR-2 order and FIR-4 order and Stochastic Calculus). MAPE0.0is the mean absolute percent error with free noise of the estimated central AP waveform. MAPE0.012 is the mean absolute percent error of the estimate with noise variance = 0.012, MAPE0.024 with noise variance = 0.024, and so forth.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,y1(k) (Femoral AP) and y2(k) (Radial AP) that are assumed to be represented by the following equations:(3.3)y1(k)=u(k)+β1u(k-1)+β2u(k-2)+β3u(k-3)+β4u(k-4),(3.4)y2(k)=u(k)+α1u(k-1)+α2u(k-2).In the conventional method, the order of the filters was assumed to be just two, that is,β3=β4=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 β2(t) has the following form:
(3.5)dβ2(t)=c2(β2(0)-β2(t))dt+e2dW2(t),
where c2, e2 are unknown constants that were estimated as explained in Sections 2.2.2 and 2.2.3. Remember that in the conventional method
(3.6)[y1(k)-y2(k)]=α1y1(k-1)-α2y1(k-2)+β1y2(k-1)+β2y2(k-2).
We estimate the unknown coefficients α1, α2, β1, and β2 through regression analysis.Figure4 shows the estimated pressure at the root of the aorta AP from the received signal y1 (Femoral AP) using the proposed method based on the stochastic calculus (OU) (SNRE0.0=23.296 db) (see Table 1). SNRE0.0 is the signal to noise ratio of the estimate with free noise. The mean absolute percent error with free noise (MAPE0.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 (SNRE0.0=9.006 db) and the proposed method OU (SNRE0.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 β̂2(t) compared to the noisy measurements of β2(k) is shown in Figure 5.Figure 5
A typical estimate forβ̂2(t) (2.25) compared to the noisy measurements of β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 I=4, J=2. Tables 1 and 2 show the estimated values of central AP using both methods, where SNRE0.012 is the signal to noise ratio of the estimate with noise variance = 0.012, SNRE0.024 with noise variance = 0.024, and so forth, and so as MAPE.Figure8 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 MAPE0.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.


## 3.1. Example  1, Noise-Free Measurements

We have two measurements,y1(k) (Femoral AP) and y2(k) (Radial AP) that are assumed to be represented by the following equations:(3.3)y1(k)=u(k)+β1u(k-1)+β2u(k-2)+β3u(k-3)+β4u(k-4),(3.4)y2(k)=u(k)+α1u(k-1)+α2u(k-2).In the conventional method, the order of the filters was assumed to be just two, that is,β3=β4=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 β2(t) has the following form:
(3.5)dβ2(t)=c2(β2(0)-β2(t))dt+e2dW2(t),
where c2, e2 are unknown constants that were estimated as explained in Sections 2.2.2 and 2.2.3. Remember that in the conventional method
(3.6)[y1(k)-y2(k)]=α1y1(k-1)-α2y1(k-2)+β1y2(k-1)+β2y2(k-2).
We estimate the unknown coefficients α1, α2, β1, and β2 through regression analysis.Figure4 shows the estimated pressure at the root of the aorta AP from the received signal y1 (Femoral AP) using the proposed method based on the stochastic calculus (OU) (SNRE0.0=23.296 db) (see Table 1). SNRE0.0 is the signal to noise ratio of the estimate with free noise. The mean absolute percent error with free noise (MAPE0.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 (SNRE0.0=9.006 db) and the proposed method OU (SNRE0.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 β̂2(t) compared to the noisy measurements of β2(k) is shown in Figure 5.Figure 5
A typical estimate forβ̂2(t) (2.25) compared to the noisy measurements of β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 I=4, J=2. Tables 1 and 2 show the estimated values of central AP using both methods, where SNRE0.012 is the signal to noise ratio of the estimate with noise variance = 0.012, SNRE0.024 with noise variance = 0.024, and so forth, and so as MAPE.Figure8 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 MAPE0.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].



---

*Source: 602373-2010-08-16.xml*