# Blind Deconvolution of the Aortic Pressure Waveform Using the Malliavin Calculus

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

---

## Abstract

Multichannel Blind Deconvolution (MBD) is a powerful tool particularly for the identification and estimation of dynamical systems in which a sensor, for measuring the input, is difficult to place. This paper presents an MBD method, based on the Malliavin calculus MC (stochastic calculus of variations). The arterial network is modeled as a Finite Impulse Response (FIR) filter with unknown coefficients. The source signal central arterial pressure CAP is also unknown. Assuming that many coefficients of the FIR filter 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 proposed Malliavin calculus-based method. We have been able to deconvolve the measurements and obtain both the source signal and the arterial path or filter. The presented examples prove the superiority of the proposed method, as compared 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, 2]. A multichannel blind deconvolution problem can be considered as natural extension or generalization of instantaneous Blind Source Separation (BSS) problem [3, 4]. The problem of BSS has received wide attention in various fields such as signal analysis and processing of speech, image [5, 6], and biomedical signals, especially, signal extraction, enhancement, denoising, model reduction, and classification problems [7–9].The MBD is the technique that allows the estimation of both an unknown input and unknown channel dynamics from only channel outputs. Although one cannot place a sensor [10] 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 MBD technique distinguishes itself from other techniques that apply a predetermined transfer function [11, 12] 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.The physiologic state of the cardiovascular (CV) system can be most accurately assessed by using the aortic blood pressure or CAP [7, 13, 14] and flow. However, standard measurement of these signals, such as catheter, entails costly and risky surgical procedures. Therefore, most of the practically applicable methods aim to monitor the CV system based on peripheral circulatory signals, for example, arterial blood pressure at a distant site. Various methods have been developed to relate the peripheral signals to the CV state. These include blind deconvolution [15–17] methods for recovering the CAP signal from the upper-limb arterial blood pressure, and the estimation of CV parameters such as left ventricular elasticity, end diastolic volume, total peripheral resistance (TPR), and mean aortic flow from arterial BP measurement [18, 19]. A chronic challenge of these previous methods is that the dynamics of the CV system, which relates the aortic and peripheral signals, is unknown and time-varying as well. So this problem turns out to be an ill-posed system identification problem because we are asked to identify both the unknown system dynamics and input signal using the output signal measurement alone.Because of the practical difficulty in measuring the arterial pressure waveform near the heart [13, 15, 20], several mathematical transformation methods have been developed based on a generalized transfer function approach [21, 22]. Over the years many methods have been suggested for blind deconvolution for estimating central aortic pressure and flow [3, 9, 19, 23, 24]. One of the most popular and effective techniques is to assume an FIR model [1] and IIR model [18, 19, 25] for the modulating channels/paths and to estimate the coefficients [26] of this model. Through the inversion of the FIR filter, we get the original source signals. The principal assumption underlying these methods is that the arterial tree properties are constant over time and between individuals. A few methods have therefore been more recently developed towards “individualizing” the transfer function [1, 12, 21] relating peripheral artery pressure to central aortic pressure. These methods essentially involve (1) modeling the transfer function [11, 23] with physiologic parameters, (2) estimating a subset of the model parameters [7, 15] from the peripheral arterial pressure waveforms and/or other measurements from an individual while assuming values for the remaining parameters, (3) constructing a transfer function based on the estimated and assumed parameter values [26], and (4) applying the transfer function to the measured peripheral arterial pressure waveforms to predict the corresponding central aortic pressure waveform. While these methods attempt to determine a transfer function that is specific to an individual over a particular time period, only a partial individualization is actually obtained. Perhaps, as a result, these methods have found only limited success with results not much, if at all, better than the generalized transfer function approach.In this paper, we suggest characterizing the channels of the single-input, multioutput system model of the arterial tree by linear and time-variant FIR filters. If we make only one of the FIR filter parameters changing over time, then the problem is handled by the Ito calculus [27, 28], while, if we make more than one of the FIR filter parameters changing over time, then the problem could be handled by the Malliavin calculus [29] as we propose in this paper. This way, the ambiguity in the order of the impulse response is compensated by the time variations of the filter parameters [29–31].In this paper, we introduce a method that could estimate time-varying parameters/coefficients. It is based on the stochastic calculus of variations (Malliavin calculus) [28, 32–35]. We derive a closed-form expression for the unknown time varying FIR filter parameters by using the Clark-Ocone formula [30]. This will enable us to find a stochastic differential equation (SDE) for each unknown time-varying parameter of the FIR filter. Each SDE is function of PAP and some other unknown deterministic parameters. The unknown deterministic parameters are estimated through Monte Carlo simulation methods.The proposed method is then applied to noninvasive monitoring of the cardiovascular system of the swine. The arterial network is modeled as a multichannel system where the CAP is the input and pressure profiles measured at different branches of the artery, for example, radial and femoral arteries, are the outputs. The proposed method would allow us to estimate both the waveform of the input pressure and the arterial channel dynamics from outputs obtained with noninvasive sensors placed at different branches of the arterial network. Numerical examples verify the major theoretical results and the feasibility of the method. In Section2, we describe the blind deconvolution problem, conventional solution methods. In Section 3 we introduce the proposed method based on the Malliavin calculus. In Section 4, we present the results for the reconstruction of single-input CAP from two distant measures outputs PAP. Finally in Section 5, we provide summary and conclusions. The Appendix contains the technical derivations of the proposed method.


## 2. Problem Formulation

### 2.1. Multichannel Dynamic Systems

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. Figure 1 illustrates, in a block diagram form, the relation between the central aortic pressure waveform u(t) and the peripheral arterial pressure waveforms yi(t). The arterial channels hi(t) relating the common input to each output represent the vascular dynamic properties of different arterial tree paths as characterized by finite impulse responses (FIRs). The main idea is therefore to determine the absolute central aortic pressure waveform within an arbitrary scale factor by mathematically analyzing two or more PAP waveforms or related signals so as to extract their common features. An ancillary idea is to then determine the parameters of the determined central aortic pressure waveform.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)].We estimate the FIR filter coefficients by the conventional methods in Section2.2. We introduce our proposed method that is based on Malliavin calculus in Section 3. We will be working in the probability space (Ω,ℱ,𝒫). To simplify the exposure, we shall assume that we have only two measurements outputs of a modulated version of the source signal that are given as(2.1)y1(t)=h1(t)*u(t)+ε1(t),(2.2)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 arterial paths, “*” is the convolution operation, y1(t) (femoral AP) and y2(t) (radial AP) are the observed measurements, and ε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.3)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) and h2(t) and on the right-hand side we get: (2.4)h2(t)*y1(t)=h1(t)*(h2(t)*u(t))+h2(t)*ε1(t)=h1(t)*y2(t)-h1(t)*ε2(t)+h2(t)*ε2(t).
Thus, (2.5)h2(t)*y1(t)=h1(t)*y2(t)-h1(t)*ε2(t)+h2(t)*ε1(t).
Note that this equation does not include the unknown 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: (2.6)h1(z)=1+β1z-1+β2z-2,
that is, (2.7)y1(k)=u(k)+β1u(k-1)+β2u(k-2)+ε1(k),
and in matrix format for N data points:(2.8)[y1(2)y(3)⋯y(N-1)]=[β2β110⋯00β2β11⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯0β2β11][u(0)u(1)⋯u(N-1)]+[ε1(2)ε1(3)⋯ε1(N-1)],
that is,
(2.9)Y̲1=H1U̲+ε̲1,
where(2.10)Y̲1=[y1(2)y1(3)⋯y1(N-1)]T,(2.11)U̲=[u(0)u(1)⋯u(N-1)]T,(2.12)ε̲1=[ε1(2)ε1(3)⋯ε1(N-1)]T,(2.13)H1=[β2β110⋯00β2β11⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯0β2β11],
and “T” stands for transpose.Similarly,(2.14)h2(z)=1+α1z-1+α2z-2,
that is, (2.15)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.16)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.17)h2(z)y1(z)=h1(z)y2(z)-h1(z)ε2(z)+h2(z)ε1(z).
In the time domain, we get the equation:
(2.18)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.2. A Conventional Method for the Estimation of the System

The familiar scalar regression format as shown in (2.18) 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.Rearrange (2.18), we gets
(2.19)a[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 are I and J, we get
(2.20)[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(2.21)[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 [36](2.22)AICc=n(lnσ2+1)+2n(q+1)(n-q-2),
where n is the number of observations, q=I+J is a number of unknown, and σ2 is an estimate of the error variance. We choose the order q such that AICc is minimized.Once the coefficients of the FIR filter are estimated, we use inverse filtering to find an estimate for the source signalU as follows: (2.23)Û̲=Ĥ1T(Ĥ1Ĥ1T)-1Y̲1,
where (2.24)Y̲1=[y1(2)y1(3)⋯y1(N-1)]T,(2.25)Û̲=[û(0)û(1)⋯û(N-1)]T,(2.26)Ĥ1=[β̂2β̂110⋯00β̂2β̂11⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯0β̂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 [19].


## 2.1. Multichannel Dynamic Systems

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. Figure 1 illustrates, in a block diagram form, the relation between the central aortic pressure waveform u(t) and the peripheral arterial pressure waveforms yi(t). The arterial channels hi(t) relating the common input to each output represent the vascular dynamic properties of different arterial tree paths as characterized by finite impulse responses (FIRs). The main idea is therefore to determine the absolute central aortic pressure waveform within an arbitrary scale factor by mathematically analyzing two or more PAP waveforms or related signals so as to extract their common features. An ancillary idea is to then determine the parameters of the determined central aortic pressure waveform.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)].We estimate the FIR filter coefficients by the conventional methods in Section2.2. We introduce our proposed method that is based on Malliavin calculus in Section 3. We will be working in the probability space (Ω,ℱ,𝒫). To simplify the exposure, we shall assume that we have only two measurements outputs of a modulated version of the source signal that are given as(2.1)y1(t)=h1(t)*u(t)+ε1(t),(2.2)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 arterial paths, “*” is the convolution operation, y1(t) (femoral AP) and y2(t) (radial AP) are the observed measurements, and ε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.3)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) and h2(t) and on the right-hand side we get: (2.4)h2(t)*y1(t)=h1(t)*(h2(t)*u(t))+h2(t)*ε1(t)=h1(t)*y2(t)-h1(t)*ε2(t)+h2(t)*ε2(t).
Thus, (2.5)h2(t)*y1(t)=h1(t)*y2(t)-h1(t)*ε2(t)+h2(t)*ε1(t).
Note that this equation does not include the unknown 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: (2.6)h1(z)=1+β1z-1+β2z-2,
that is, (2.7)y1(k)=u(k)+β1u(k-1)+β2u(k-2)+ε1(k),
and in matrix format for N data points:(2.8)[y1(2)y(3)⋯y(N-1)]=[β2β110⋯00β2β11⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯0β2β11][u(0)u(1)⋯u(N-1)]+[ε1(2)ε1(3)⋯ε1(N-1)],
that is,
(2.9)Y̲1=H1U̲+ε̲1,
where(2.10)Y̲1=[y1(2)y1(3)⋯y1(N-1)]T,(2.11)U̲=[u(0)u(1)⋯u(N-1)]T,(2.12)ε̲1=[ε1(2)ε1(3)⋯ε1(N-1)]T,(2.13)H1=[β2β110⋯00β2β11⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯0β2β11],
and “T” stands for transpose.Similarly,(2.14)h2(z)=1+α1z-1+α2z-2,
that is, (2.15)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.16)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.17)h2(z)y1(z)=h1(z)y2(z)-h1(z)ε2(z)+h2(z)ε1(z).
In the time domain, we get the equation:
(2.18)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.2. A Conventional Method for the Estimation of the System

The familiar scalar regression format as shown in (2.18) 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.Rearrange (2.18), we gets
(2.19)a[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 are I and J, we get
(2.20)[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(2.21)[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 [36](2.22)AICc=n(lnσ2+1)+2n(q+1)(n-q-2),
where n is the number of observations, q=I+J is a number of unknown, and σ2 is an estimate of the error variance. We choose the order q such that AICc is minimized.Once the coefficients of the FIR filter are estimated, we use inverse filtering to find an estimate for the source signalU as follows: (2.23)Û̲=Ĥ1T(Ĥ1Ĥ1T)-1Y̲1,
where (2.24)Y̲1=[y1(2)y1(3)⋯y1(N-1)]T,(2.25)Û̲=[û(0)û(1)⋯û(N-1)]T,(2.26)Ĥ1=[β̂2β̂110⋯00β̂2β̂11⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯0β̂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 [19].


## 3. The Malliavin Calculus and the Generalized Clark-Ocone Formula for the Estimation of the Unknown Time-Varying Coefficients: The Proposed Approach

The linear filter assumption is just an approximation to reality. Sometimes the media (arterial tree) is nonlinear, time varying, random, or all. Moreover and above the measured signals are usually noisy. The filter order is usually unknown. All the factors suggest that the FIR filter model is just an approximation. To compensate for these approximations, we suggest making some of the unknown filter coefficients varying with time, for example, (α1(t),α2(t),β1(t),β2(t)). Notice that in this case, the filter coefficients do not have any physical meaning and we are only concerned with the estimation of their values. Now the problem becomes that of the estimation of the unknown time-varying coefficients. The details of the estimation procedure are, briefly, given in this section. More details could be found in the appendix.We now recast the problem in the format that could be handled by the Malliavin calculus. Using (2.21), the observed signal/regressor, py(k)=y1(k)-y2(k), is the sum of stochastic processes, y1(k-1),y1(k-2),…,y2(k-1),y2(k-2),… weighted by the unknown time-varying parameters, α1(k),α2(k),…,β1(k),β2(k),…. The objective is to estimate, from the observations, the unknown time-varying coefficients α1(k),α2(k),…,β1(k),β2(k),… .


### 3.1. The Estimates for the Time-Varying Filter Parameters

To solve the estimation of the time-varying parameters problem, we imbed the sum of the processes (the observations) into another signal,p(t), which we call the augmented observed signal. The augmented observations are the original observations plus a deterministic component. The addition of the known deterministic component is needed to facilitate the analysis.The augmented observed signal,p(t), with (d) components, could be modeled, in the continuous time, as follows:(3.1)p(t)=η1(t)S1(t)+∑i=2dηi(t)Si(t)t∈[0,T]=η1(t)S1(t)+py(t),
where py(t)=[y1(t)-y2(t)]; see (2.21). ηi(t),i>1, is the ith unknown time-varying coefficient and are defined as follows for second-order filters:η2(t)=-α1(t),η3(t)=-α2(t),η4(t)=β1(t),η5(t)=β2(t),S2(t)=y1(t-Δ),S3(t)=y1(t-2Δ),S4(t)=y2(t-Δ),S5(t)=y2(t-2Δ),Δ is the sampling interval, η1(t) = known constant, with (3.2)η1(t)=knownconstant,(3.3)dS1(t)=λ1S1(t)dt,
where λ1 and S1(t) are known.In this analysis, the presence of a deterministic component will make the exposition easier. In this case, It isS1(t). This component acts as a reference signal or a numeraire. Since the observed signal does not usually come with a known deterministic component, we should add a known deterministic component and proceed with the analysis. We will then try to find an expression for the SDE that describes the evolution of p(t). This expression will be a function of the unknown stochastic coefficients ηi(t). Assuming a shape for the final value p(T), as a function of the Wiener processes W̲(t), and using the generalized Clark-Ocone formula [30, 37, 38], we will be able to find a closed-form expression for the coefficients ηi(t) as a function of the Wiener processes W̲(t). The details of the derivations are given in the appendix and are also given in [29].Assume that the unknown coefficients,ηi(t), are varying slowly with time, that is(3.4)dηi(t)→0,
or more precisely, dηi(t)≪dSi(t).There are several commonly models that represent different physical situations. In our case, the signals,Si(t), each has an SDE of the form of an Ornstein-Uhlenbeck (OU) process with trend(3.5)dSi(t)=ci(t)[Oi(t)-Si(t)]dt+ei(t)dWi(t),i>1.
This form represents an Ornstein-Uhlenbeck (OU) process and has the solution [12, 39](3.6)Si(t)=e-Ci(t)[Si(0)+∫0tci(u)Oi(u)eCi(u)du]+e-Ci(t)∫0tei(u)eCi(u)dWi(u),i>1,
where Ci(t)=∫0tci(u)du.The coefficientsci(t),ei(t), and the trend Oi(t) are estimated from the observed data.  For constant quantities ci,Oi,ei, and after some manipulations we end up with the estimate of the unknown time-varying parameters (see the appendix):(3.7)ηi(t)=eλ1teiEQ{[DtiF(T)-F(T)∫tTDtiθi(s,Wi(s))dW̃i(s)]ℱt},
where
(3.8)θi(u,Wi(u))=ci-λ1ei[Si(0)eciu+∫0ueieci(u-r)dWi(r)],(3.9)F(T)=e-λ1Tp(T),EQ{X/ℱt} is the conditional expectation of the random variable X under the probability measure Q [37] and X=[DtiF(T)-F(T)∫tTDtiθi(s,Wi(s))dW̃i(s)]. This measure Q is related to the old probability measure 𝒫 through a function of the stochastic process θi(u,Wi(u)).ℱt is the filtration generated by the Wiener vector W̲(t). t∈[0,T]. ℱt contains exactly the information learned by observing the Wiener vector processes up to time t. Dtiθi(s,Wi(s))is the Malliavin derivative of θi(s) with respect to the ith Wiener process Wi(s)·F(T) is, in general, a function of the Wiener vector W̲(T). EQ{X/ℱt} is usually difficult to calculate. ηi(t) is in general a function of W̲(t), and it could be approximated by the SDE:(3.10)dηi(t)=μi(W̲(t),t)dt+∑jκij(W̲(t),t)dWj(t),
and eiηi(t) could be represented by the SDE(3.11)d[eiηi(t)]=eiμi(W̲(t),t)dt+ei∑jκij(W̲(t),t)dWj(t).
Notice that the drift and the diffusion are related to each other through θi(t,Wi(t)) [29, 32] see the appendix.  Let
(3.12)eiμi(W̲(t),t)=∑jϕij(t,W̲(t))θj(t,W̲(t)),
then
(3.13)eiκij(W̲(t),t)=ϕij(t,W̲(t)).
The last equality comes from the solution of the CO formula θi(t,Wi(t)) (see the appendix). An expression for ϕij(t,W̲(t)) is needed to solve for the unknowns ηi(t). An exact expression, however, is difficult to obtain. Instead, we assume a polynomial shape with unknown parameters (δij,γkij,λklij), namely,(3.14)ϕij(t,W̲(t))=δij+∑kγkijWk(t)+∑k∑lλklijWk(t)Wl(t)+⋯.
The unknown parameters are estimated by minimizing the squared error, e2(T), between the observed, p(t), and the estimated observations p̂(t), where(3.15)p̂(t)=∑iη̂i(t)Si(t),i≥1,(3.16)d[eiη̂i(t)]=eiμ̂i(W̲(t),t)dt+ei∑jκ̂ij(W̲(t),t)dWj(t),(3.17)eiμ̂i(W̲(t),t)=∑jϕ̂ij(t,W̲(t))θj(t,W̲(t)),(3.18)eiκ̂ij(W̲(t),t)=ϕ̂ij(t,W̲(t)),(3.19)ϕ̂ij(t,W̲(t))=δij+∑kγkijWk(t)+∑k∑lλklijWk(t)Wl(t),(3.20)e2(T)=∫0T[p(t)-p̂(t)]2dt.Once the parametersηi(t) are estimated, we use inverse filtering to find U (the aortic pressure waveform). In the case that the filter order is two, that is, only two time-varying parameters, η̂2(t)=α̂1(t),η̂3(t)=α̂2(t),η̂4(t)=β̂1(t), and η̂5(t)=β̂2(t). In this case, the aortic pressure Û is estimated as
(3.21)Û̲=Ĥ1T(Ĥ1Ĥ1T)-1Y̲1,
where
(3.22)Y̲1=[y1(2)y1(3)⋯y1(N-1)]T,Û̲=[û(0)û(1)⋯û(N-1)]T,Ĥ1=[β̂2(t1)β̂1(t2)10⋯00β̂2(t2)β̂1(t3)1⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯0β̂2(tN-2)β̂1(tN-1)1].
We note from matrix that the parameters β̂1(t) and β̂2(t) have values that are changing across the sample.In order to further reduce the mathematical complexity, one could restrict the analysis to the estimation of only two coefficientsα1(t) and β1(t). The rest of filters coefficients are considered constants and the OLS estimates of them are used in the equations.


### 3.2. Summary of the Proposed Algorithm and its Assumptions

Our MBD technique is based on a set of assumptions as follows.(1)
The common input is obtained at the output of two sensors. The channels relating the common input to each distinct output, as in Figure1, are linear time variant.(2)
The system is represented by FIR filters. The filtersh1(t) and h2(t) are second-order and contain many parameters some of which are time-varying.(3)
The FIR filter parametersα1(t),α2(t),β1(t),andβ2(t) are varying slowly with time.(4)
An OU model is assumed to describe the signals/regressorsSi(t) which represent the processes y1(t-Δ),y1(t-2Δ),y2(t-Δ), and y2(t-2Δ)); see (3.1) and (3.6).(5)
The time-varying FIR filter parametersηi(t) are estimated by using Malliavin calculus and the generalization Clark-Ocone formula; see (3.7).(6)
Equataion (3.21) is used with the time-varying parameters/coefficients η̂i(t) to find an estimate for the input AP.An absolute central aortic pressure waveform is determined by scaling the estimated input based on the measured waveforms. In certain embodiments, the reconstructed waveform is calibrated to absolute pressure based on the measured peripheral artery pressure waveforms. For example, the reconstructed waveform is scaled to have a mean value similar to the mean value of the measured waveforms. This scaling step is well justified, since the paths from the central aorta to peripheral arteries offer very little resistance to blood flow due to Poiseuille's law [1]. Certain embodiments scale the reconstructed waveform to have a mean value specifically equal to that of one of the measured peripheral artery pressure waveforms or the waveform with the largest mean value plus a constant (whose value may be between, e.g., 0 and 3 mmHg). Certain alternative embodiments scale the reconstructed waveform to have a mean value equal to the mean (or medium) of the mean values of some or all of the measured peripheral artery pressure waveforms plus a constant (whose value may be between, e.g., 0 and 3 mmHg).The reconstructed waveformÛ̲ is calibrated to absolute pressure by scaling it to have the same mean value as the measured peripheral artery pressures as follows:(3.23)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, ÛS(t)=Û(t) × Scale factor, (3.24)Scalefactor=averagefemoralpressure/2+averageradialpressure/2averageestimatedaorticpressure.


### 3.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 (3.5).(3)
Estimate the parametersηi(t) by the following steps for each unknown time-varying coefficient ηi(t).
(i)
Initiate an estimate for the unknown parametersδij,γkij,λklij.(ii)
Use Monte Carlo simulation method to generate Wiener processesW̲(t).(iii)
Use (3.19) to an estimate for ϕ̂ij(t,W̲(t)).(iv)
Use (3.16) to find an estimate for η̂i(t).(v)
Calculate the summed squared error of (3.20).(vi)
Change the estimate for the unknown parametersδij,γkij,andλklij using the least square method, for example, such that the error is reduced.(vii)
Repeat steps (ii)–(vi) and stop when the error is minimum or the number of iterations is exceeded.(4)
Calculate the matrixĤ1 in (3.25).(5)
Estimate the source signal by (3.21).(6)
Calculate the estimated aortic pressure within scale factor by (3.23).Figure 2
Segments of measured central arterial pressure (AP), femoral AP, and radial AP waveforms from one swine dataset.Figure 3
The arterial pressure measurements at displayed points in the figure.


## 3.1. The Estimates for the Time-Varying Filter Parameters

To solve the estimation of the time-varying parameters problem, we imbed the sum of the processes (the observations) into another signal,p(t), which we call the augmented observed signal. The augmented observations are the original observations plus a deterministic component. The addition of the known deterministic component is needed to facilitate the analysis.The augmented observed signal,p(t), with (d) components, could be modeled, in the continuous time, as follows:(3.1)p(t)=η1(t)S1(t)+∑i=2dηi(t)Si(t)t∈[0,T]=η1(t)S1(t)+py(t),
where py(t)=[y1(t)-y2(t)]; see (2.21). ηi(t),i>1, is the ith unknown time-varying coefficient and are defined as follows for second-order filters:η2(t)=-α1(t),η3(t)=-α2(t),η4(t)=β1(t),η5(t)=β2(t),S2(t)=y1(t-Δ),S3(t)=y1(t-2Δ),S4(t)=y2(t-Δ),S5(t)=y2(t-2Δ),Δ is the sampling interval, η1(t) = known constant, with (3.2)η1(t)=knownconstant,(3.3)dS1(t)=λ1S1(t)dt,
where λ1 and S1(t) are known.In this analysis, the presence of a deterministic component will make the exposition easier. In this case, It isS1(t). This component acts as a reference signal or a numeraire. Since the observed signal does not usually come with a known deterministic component, we should add a known deterministic component and proceed with the analysis. We will then try to find an expression for the SDE that describes the evolution of p(t). This expression will be a function of the unknown stochastic coefficients ηi(t). Assuming a shape for the final value p(T), as a function of the Wiener processes W̲(t), and using the generalized Clark-Ocone formula [30, 37, 38], we will be able to find a closed-form expression for the coefficients ηi(t) as a function of the Wiener processes W̲(t). The details of the derivations are given in the appendix and are also given in [29].Assume that the unknown coefficients,ηi(t), are varying slowly with time, that is(3.4)dηi(t)→0,
or more precisely, dηi(t)≪dSi(t).There are several commonly models that represent different physical situations. In our case, the signals,Si(t), each has an SDE of the form of an Ornstein-Uhlenbeck (OU) process with trend(3.5)dSi(t)=ci(t)[Oi(t)-Si(t)]dt+ei(t)dWi(t),i>1.
This form represents an Ornstein-Uhlenbeck (OU) process and has the solution [12, 39](3.6)Si(t)=e-Ci(t)[Si(0)+∫0tci(u)Oi(u)eCi(u)du]+e-Ci(t)∫0tei(u)eCi(u)dWi(u),i>1,
where Ci(t)=∫0tci(u)du.The coefficientsci(t),ei(t), and the trend Oi(t) are estimated from the observed data.  For constant quantities ci,Oi,ei, and after some manipulations we end up with the estimate of the unknown time-varying parameters (see the appendix):(3.7)ηi(t)=eλ1teiEQ{[DtiF(T)-F(T)∫tTDtiθi(s,Wi(s))dW̃i(s)]ℱt},
where
(3.8)θi(u,Wi(u))=ci-λ1ei[Si(0)eciu+∫0ueieci(u-r)dWi(r)],(3.9)F(T)=e-λ1Tp(T),EQ{X/ℱt} is the conditional expectation of the random variable X under the probability measure Q [37] and X=[DtiF(T)-F(T)∫tTDtiθi(s,Wi(s))dW̃i(s)]. This measure Q is related to the old probability measure 𝒫 through a function of the stochastic process θi(u,Wi(u)).ℱt is the filtration generated by the Wiener vector W̲(t). t∈[0,T]. ℱt contains exactly the information learned by observing the Wiener vector processes up to time t. Dtiθi(s,Wi(s))is the Malliavin derivative of θi(s) with respect to the ith Wiener process Wi(s)·F(T) is, in general, a function of the Wiener vector W̲(T). EQ{X/ℱt} is usually difficult to calculate. ηi(t) is in general a function of W̲(t), and it could be approximated by the SDE:(3.10)dηi(t)=μi(W̲(t),t)dt+∑jκij(W̲(t),t)dWj(t),
and eiηi(t) could be represented by the SDE(3.11)d[eiηi(t)]=eiμi(W̲(t),t)dt+ei∑jκij(W̲(t),t)dWj(t).
Notice that the drift and the diffusion are related to each other through θi(t,Wi(t)) [29, 32] see the appendix.  Let
(3.12)eiμi(W̲(t),t)=∑jϕij(t,W̲(t))θj(t,W̲(t)),
then
(3.13)eiκij(W̲(t),t)=ϕij(t,W̲(t)).
The last equality comes from the solution of the CO formula θi(t,Wi(t)) (see the appendix). An expression for ϕij(t,W̲(t)) is needed to solve for the unknowns ηi(t). An exact expression, however, is difficult to obtain. Instead, we assume a polynomial shape with unknown parameters (δij,γkij,λklij), namely,(3.14)ϕij(t,W̲(t))=δij+∑kγkijWk(t)+∑k∑lλklijWk(t)Wl(t)+⋯.
The unknown parameters are estimated by minimizing the squared error, e2(T), between the observed, p(t), and the estimated observations p̂(t), where(3.15)p̂(t)=∑iη̂i(t)Si(t),i≥1,(3.16)d[eiη̂i(t)]=eiμ̂i(W̲(t),t)dt+ei∑jκ̂ij(W̲(t),t)dWj(t),(3.17)eiμ̂i(W̲(t),t)=∑jϕ̂ij(t,W̲(t))θj(t,W̲(t)),(3.18)eiκ̂ij(W̲(t),t)=ϕ̂ij(t,W̲(t)),(3.19)ϕ̂ij(t,W̲(t))=δij+∑kγkijWk(t)+∑k∑lλklijWk(t)Wl(t),(3.20)e2(T)=∫0T[p(t)-p̂(t)]2dt.Once the parametersηi(t) are estimated, we use inverse filtering to find U (the aortic pressure waveform). In the case that the filter order is two, that is, only two time-varying parameters, η̂2(t)=α̂1(t),η̂3(t)=α̂2(t),η̂4(t)=β̂1(t), and η̂5(t)=β̂2(t). In this case, the aortic pressure Û is estimated as
(3.21)Û̲=Ĥ1T(Ĥ1Ĥ1T)-1Y̲1,
where
(3.22)Y̲1=[y1(2)y1(3)⋯y1(N-1)]T,Û̲=[û(0)û(1)⋯û(N-1)]T,Ĥ1=[β̂2(t1)β̂1(t2)10⋯00β̂2(t2)β̂1(t3)1⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯0β̂2(tN-2)β̂1(tN-1)1].
We note from matrix that the parameters β̂1(t) and β̂2(t) have values that are changing across the sample.In order to further reduce the mathematical complexity, one could restrict the analysis to the estimation of only two coefficientsα1(t) and β1(t). The rest of filters coefficients are considered constants and the OLS estimates of them are used in the equations.


## 3.2. Summary of the Proposed Algorithm and its Assumptions

Our MBD technique is based on a set of assumptions as follows.(1)
The common input is obtained at the output of two sensors. The channels relating the common input to each distinct output, as in Figure1, are linear time variant.(2)
The system is represented by FIR filters. The filtersh1(t) and h2(t) are second-order and contain many parameters some of which are time-varying.(3)
The FIR filter parametersα1(t),α2(t),β1(t),andβ2(t) are varying slowly with time.(4)
An OU model is assumed to describe the signals/regressorsSi(t) which represent the processes y1(t-Δ),y1(t-2Δ),y2(t-Δ), and y2(t-2Δ)); see (3.1) and (3.6).(5)
The time-varying FIR filter parametersηi(t) are estimated by using Malliavin calculus and the generalization Clark-Ocone formula; see (3.7).(6)
Equataion (3.21) is used with the time-varying parameters/coefficients η̂i(t) to find an estimate for the input AP.An absolute central aortic pressure waveform is determined by scaling the estimated input based on the measured waveforms. In certain embodiments, the reconstructed waveform is calibrated to absolute pressure based on the measured peripheral artery pressure waveforms. For example, the reconstructed waveform is scaled to have a mean value similar to the mean value of the measured waveforms. This scaling step is well justified, since the paths from the central aorta to peripheral arteries offer very little resistance to blood flow due to Poiseuille's law [1]. Certain embodiments scale the reconstructed waveform to have a mean value specifically equal to that of one of the measured peripheral artery pressure waveforms or the waveform with the largest mean value plus a constant (whose value may be between, e.g., 0 and 3 mmHg). Certain alternative embodiments scale the reconstructed waveform to have a mean value equal to the mean (or medium) of the mean values of some or all of the measured peripheral artery pressure waveforms plus a constant (whose value may be between, e.g., 0 and 3 mmHg).The reconstructed waveformÛ̲ is calibrated to absolute pressure by scaling it to have the same mean value as the measured peripheral artery pressures as follows:(3.23)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, ÛS(t)=Û(t) × Scale factor, (3.24)Scalefactor=averagefemoralpressure/2+averageradialpressure/2averageestimatedaorticpressure.


## 3.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 (3.5).(3)
Estimate the parametersηi(t) by the following steps for each unknown time-varying coefficient ηi(t).
(i)
Initiate an estimate for the unknown parametersδij,γkij,λklij.(ii)
Use Monte Carlo simulation method to generate Wiener processesW̲(t).(iii)
Use (3.19) to an estimate for ϕ̂ij(t,W̲(t)).(iv)
Use (3.16) to find an estimate for η̂i(t).(v)
Calculate the summed squared error of (3.20).(vi)
Change the estimate for the unknown parametersδij,γkij,andλklij using the least square method, for example, such that the error is reduced.(vii)
Repeat steps (ii)–(vi) and stop when the error is minimum or the number of iterations is exceeded.(4)
Calculate the matrixĤ1 in (3.25).(5)
Estimate the source signal by (3.21).(6)
Calculate the estimated aortic pressure within scale factor by (3.23).Figure 2
Segments of measured central arterial pressure (AP), femoral AP, and radial AP waveforms from one swine dataset.Figure 3
The arterial pressure measurements at displayed points in the figure.


## 4. Computed Examples

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, radial artery pressure 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 Excel by using Risk Solver Platform V9.0 [40, 41]. Fourth, we evaluated the proposed method by the performance measure SNRE. The signal-to-noise ratio of the estimates (SNRE) was taken as the measure of performance for this evaluation. It is defined as(4.1)SNRE=10log∑ku2(k)∑k[u(i)-û(i)]2db,
where û(i) is the estimated value of the pressure at instant “i” (see Table 1).Table 1
The performance measure as the signal-to-noise ratio of the aortic pressure estimates (SNRE) using the conventional FIR-based method (FIR-2 model with two orders and FIR-4 model with 4 orders) and the proposed Malliavin-calculus-based method.
MethodsFIR-2FIR-4Malliavin calculusSNRE15.125 db18.754 db25.111 dbIn 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 was assumed noisy and the order of each of the filters was unknown. The filter order was estimated using a corrected Akaike information criterion (AIC) [36]. In both examples, by using the proposed method, all coefficients of the two filters were assumed to be unknown slowly time-varying, η2(t),…,ηd(t). It was assumed that the signals/regressors Si(t) follow an Ornstein-Uhlenbeck (OU) process [31] with some unknown constant parameters see; (3.5).In both cases, the proposed approach outperformed the conventional method. The pressure at the root of the aorta (central AP) was successfully estimated.


### 4.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 equations (4.2)y1(k)=u(k)+β1u(k-1)+β2u(k-2),(4.3)y2(k)=u(k)+α1u(k-1)+α2u(k-2)+α3u(k-3)+α4u(k-4).
In the conventional method, the order of the filters was assumed to be just two, that is, α3=α3=0. In the proposed method, we made the same order assumption of two but all coefficients are slowly time-varying. Specifically, we assumed that the SDE of the signals, Si(t), of each has the form of the Ornstein-Uhlenbeck (OU) process (4.4)dSi(t)=ci[Oi(t)-Si(t)]dt+eidWi(t),i>1,
where ci,ei are unknown constants, Oi(t) is a polynomial in time that represents the trend or the base data, and Si(t) bounces around Oi(t). Remember that in the conventional method (4.5)[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. In the proposed method, however, the slowly time-varying coefficients η̂i(t), i=2,3,4,5; are estimated by the (3.16) in Figure 5 and assumed η̂1(t)=1.Figure4 shows the estimated pressure at the root of the aorta AP using the proposed method based on the Malliavin calculus (Section 3.2) and measured aorta or central AP, (SNRE = 25.111 db) (see Table 1). In Figure 7, we compared the estimated central arterial pressure using the conventional method (Section 2.2) with second-order FIR-2 (SNRE = 15.125 db) and the proposed Malliavin calculus-based method (SNRE = 25.111 db). So applying the conventional method with fourth-order FIR-4, we obtained a better result than by using FIR-2. But still the Malliavin calculus (MC) proved to be the best (see Table 1 and Figure 7). We note the difference between three methods in Figures 6 and 7. Therefore, typical estimates for η̂2(t),η̂3(t),η̂4(t),η̂5(t) using the Malliavin calculus are shown in Figure 5, where η̂2(t)=α̂1(t),η̂3(t)=α̂2(t),η̂4(t)=β̂1(t), and η̂5(t)=β̂2(t).Figure 4
The estimated aorta AP using Malliavin calculus and measured aorta or central AP waveforms.Typical estimate ofηi, i>1, using Malliavin calculus method.
(a)(b)(c)(d)Figure 6
The performace measure SNRE of the estimated noise aorta AP at variant noise by using Malliavin and FIR_4 and FIR-2 methods.Figure 7
Estimated aorta AP waveform using the conventional methods (FIR-2 model with two orders and FIR-4 model with four orders) and proposed method (Malliavin calculus).


### 4.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 (4.2) and (4.3). 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. Table 2 shows the estimated values of central AP using both methods.Table 2
The performance measure (SNRE) after added white Gaussian noise.
FIR-2FIR-4Malliavin calculusNoise varianceSNRSNRESNRESNRE138.199 db13.8976 db18.599 db24.927 db1028.217 db13.479 db18.356 db23.443 db7019.833 db12.748 db16.548 db18.808 db10018.312 db12.441 db15.767 db17.664 db30013.705 db10.851 db13.044 db13.768 db50011.633 db9.7677 db11.466 db11.803 db70010.309 db8.9428 db10.374 db10.523 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 (Malliavin calculus). We compared between these methods by using the performance measure SNRE (3.1) with the signal-to-noise ratio SNR at variant noise.Figure9 shows the measured pressure at the root of the aorta (with added high noise) and the estimated pressure using the proposed method that is based on the Malliavin calculus. By using the proposed method based on Malliavin calculus (stochastic calculus of variations), the SNRE (as shown in Table 2) at noise variance 70 = 18.808 db, where the signal-to-noise ratio (SNR = 19.833 db). While, in Figure 8, the estimated central arterial pressure (with added the same noise variance) using the conventional method with second-order FIR-2 (SNRE = 12.748 db) and with four order FIR-4 (SNRE = 16.548 db). The performance measure (as shown in Figure 6) of the FIR filter-based methods and the proposed method are compared in Table 2 for different noise levels. The results show the superior performance of the proposed method.Figure 8
Estimated noisy aorta AP waveforms using the conventional methods FIR-4 model  with 4 orders and FIR-2 Model with 2 Orders. (SNR = 19.833 db, SNRE = 16.548 db)Figure 9
The estimated aorta AP using Malliavin calculus at high noise variance and measured aorta or central AP waveforms. (SNR = 19.833 db, SNRE = 18.808 db).


## 4.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 equations (4.2)y1(k)=u(k)+β1u(k-1)+β2u(k-2),(4.3)y2(k)=u(k)+α1u(k-1)+α2u(k-2)+α3u(k-3)+α4u(k-4).
In the conventional method, the order of the filters was assumed to be just two, that is, α3=α3=0. In the proposed method, we made the same order assumption of two but all coefficients are slowly time-varying. Specifically, we assumed that the SDE of the signals, Si(t), of each has the form of the Ornstein-Uhlenbeck (OU) process (4.4)dSi(t)=ci[Oi(t)-Si(t)]dt+eidWi(t),i>1,
where ci,ei are unknown constants, Oi(t) is a polynomial in time that represents the trend or the base data, and Si(t) bounces around Oi(t). Remember that in the conventional method (4.5)[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. In the proposed method, however, the slowly time-varying coefficients η̂i(t), i=2,3,4,5; are estimated by the (3.16) in Figure 5 and assumed η̂1(t)=1.Figure4 shows the estimated pressure at the root of the aorta AP using the proposed method based on the Malliavin calculus (Section 3.2) and measured aorta or central AP, (SNRE = 25.111 db) (see Table 1). In Figure 7, we compared the estimated central arterial pressure using the conventional method (Section 2.2) with second-order FIR-2 (SNRE = 15.125 db) and the proposed Malliavin calculus-based method (SNRE = 25.111 db). So applying the conventional method with fourth-order FIR-4, we obtained a better result than by using FIR-2. But still the Malliavin calculus (MC) proved to be the best (see Table 1 and Figure 7). We note the difference between three methods in Figures 6 and 7. Therefore, typical estimates for η̂2(t),η̂3(t),η̂4(t),η̂5(t) using the Malliavin calculus are shown in Figure 5, where η̂2(t)=α̂1(t),η̂3(t)=α̂2(t),η̂4(t)=β̂1(t), and η̂5(t)=β̂2(t).Figure 4
The estimated aorta AP using Malliavin calculus and measured aorta or central AP waveforms.Typical estimate ofηi, i>1, using Malliavin calculus method.
(a)(b)(c)(d)Figure 6
The performace measure SNRE of the estimated noise aorta AP at variant noise by using Malliavin and FIR_4 and FIR-2 methods.Figure 7
Estimated aorta AP waveform using the conventional methods (FIR-2 model with two orders and FIR-4 model with four orders) and proposed method (Malliavin calculus).


## 4.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 (4.2) and (4.3). 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. Table 2 shows the estimated values of central AP using both methods.Table 2
The performance measure (SNRE) after added white Gaussian noise.
FIR-2FIR-4Malliavin calculusNoise varianceSNRSNRESNRESNRE138.199 db13.8976 db18.599 db24.927 db1028.217 db13.479 db18.356 db23.443 db7019.833 db12.748 db16.548 db18.808 db10018.312 db12.441 db15.767 db17.664 db30013.705 db10.851 db13.044 db13.768 db50011.633 db9.7677 db11.466 db11.803 db70010.309 db8.9428 db10.374 db10.523 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 (Malliavin calculus). We compared between these methods by using the performance measure SNRE (3.1) with the signal-to-noise ratio SNR at variant noise.Figure9 shows the measured pressure at the root of the aorta (with added high noise) and the estimated pressure using the proposed method that is based on the Malliavin calculus. By using the proposed method based on Malliavin calculus (stochastic calculus of variations), the SNRE (as shown in Table 2) at noise variance 70 = 18.808 db, where the signal-to-noise ratio (SNR = 19.833 db). While, in Figure 8, the estimated central arterial pressure (with added the same noise variance) using the conventional method with second-order FIR-2 (SNRE = 12.748 db) and with four order FIR-4 (SNRE = 16.548 db). The performance measure (as shown in Figure 6) of the FIR filter-based methods and the proposed method are compared in Table 2 for different noise levels. The results show the superior performance of the proposed method.Figure 8
Estimated noisy aorta AP waveforms using the conventional methods FIR-4 model  with 4 orders and FIR-2 Model with 2 Orders. (SNR = 19.833 db, SNRE = 16.548 db)Figure 9
The estimated aorta AP using Malliavin calculus at high noise variance and measured aorta or central AP waveforms. (SNR = 19.833 db, SNRE = 18.808 db).


## 5. Discussion and Conclusions

We have developed a new technique to estimate the aortic pressure waveform from multiple measured peripheral artery pressure waveforms. The technique is based on multichannel blind deconvolution in which two or more measured outputs (peripheral artery pressure waveforms) of a single input, multioutput system (arterial tree) are mathematically analyzed so as to reconstruct the common unobserved input (aortic pressure waveform). Each channel is modeled as an FIR filter. We assumed that all of the FIR filter parameters/coefficients are varying slowly with time. Their values were estimated using methods based on the Malliavin calculus. By this assumption, time-varying parameters, we were able to compensate for the wrong FIR filter order and the possible time variations/nonlinearities of the channels. The blind deconvolution problem was reformulated as a regression problem with unknown time-varying regression coefficients. We introduced a new method for the estimation of these slowly time-varying regression coefficients. The method relies heavily on the Malliavin calculus (stochastic calculus of variations) and the generalized Clark-Ocone formula to find a closed-form expression for the estimates of the unknown time-varying coefficients. Some approximations were needed to find mathematically tractable estimates. While the approach is quite general and can be applied to any signal model, we present only one signal model, the Ornstein Uhlenbeck process. Other models could have been used as well.We tested the proposed technique in swine experiments, and our results showed superior performance for our proposed approach compared to conventional methods. 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. We suggest studying the case where the time-varying parameters/coefficients are rapidly changing over time. Another extension is the estimation of the flow at the root of the Aorta from peripheral measurements of the flow/pressure.



---

*Source: 102581-2010-09-23.xml*