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.
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 () modeling the transfer function [11, 23] with physiologic parameters, () 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, () constructing a transfer function based on the estimated and assumed parameter values [26], and () 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 Section 2, 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 Figure 3) 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 and the peripheral arterial pressure waveforms . The arterial channels 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.
We estimate the FIR filter coefficients by the conventional methods in Section 2.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 where is the unknown source signal (central AP), and are unknown filters (hemodynamic response at time ) or arterial paths, “” is the convolution operation, (femoral AP) and (radial AP) are the observed measurements, and and are the measurements noise. The objective is to deconvolve and to estimate . If we convolve with , we will get since the convolution is a commutative operation, then exchanging and and on the right-hand side we get: Thus, Note that this equation does not include the unknown input . It represents the constraints among the channel dynamics or filters and observed output. Substituting a measured time series of output data for and , the above equation can be solved for the unknown parameters involved in and . Once the filters are obtained, we will use their inverses to find an estimate for the source signal. To simplify the exposure further, assume that the modulating filters, that represent the signal paths or channel dynamics, are second-order linear time invariant and have the transforms: that is, and in matrix format for data points: that is, where and “” stands for transpose.
Similarly, that is, where and are the transforms of the discrete versions of and , respectively. Since, taking the transform of the discrete version of both sides, we get In the time domain, we get the equation:
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 .
Rearrange (2.18), we gets For the general case where the order of the FIR filters are and , we get This could be approximated as Since the filter orders are unknown, one could use the corrected Akaike information criterion to determine both “” and “”. Assume that the error term, , is zero mean Gaussian with variance , the is defined as [36] where is the number of observations, is a number of unknown, and is an estimate of the error variance. We choose the order such that is minimized.
Once the coefficients of the FIR filter are estimated, we use inverse filtering to find an estimate for the source signal as follows: where The symbol “” on top of the variable refers to estimation. For example, and are the estimates for and respectively. Due to its simplicity, the above mentioned method is the one that is commonly used [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, (). 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, , is the sum of stochastic processes, weighted by the unknown time-varying parameters, . The objective is to estimate, from the observations, the unknown time-varying coefficients .
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, , 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, , with components, could be modeled, in the continuous time, as follows: where ; see (2.21). , is the unknown time-varying coefficient and are defined as follows for second-order filters:
is the sampling interval, = known constant, with where and are known.
In this analysis, the presence of a deterministic component will make the exposition easier. In this case, It is . 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 . This expression will be a function of the unknown stochastic coefficients . Assuming a shape for the final value , as a function of the Wiener processes , and using the generalized Clark-Ocone formula [30, 37, 38], we will be able to find a closed-form expression for the coefficients as a function of the Wiener processes . The details of the derivations are given in the appendix and are also given in [29].
Assume that the unknown coefficients, , are varying slowly with time, that is or more precisely, .
There are several commonly models that represent different physical situations. In our case, the signals, , each has an SDE of the form of an Ornstein-Uhlenbeck (OU) process with trend This form represents an Ornstein-Uhlenbeck (OU) process and has the solution [12, 39] where .
The coefficients , and the trend are estimated from the observed data. For constant quantities , and after some manipulations we end up with the estimate of the unknown time-varying parameters (see the appendix): where is the conditional expectation of the random variable X under the probability measure [37] and . This measure is related to the old probability measure through a function of the stochastic process is the filtration generated by the Wiener vector . . contains exactly the information learned by observing the Wiener vector processes up to time . is the Malliavin derivative of with respect to the Wiener process is, in general, a function of the Wiener vector . is usually difficult to calculate. is in general a function of , and it could be approximated by the SDE: and could be represented by the SDE Notice that the drift and the diffusion are related to each other through [29, 32] see the appendix. Let then The last equality comes from the solution of the CO formula (see the appendix). An expression for is needed to solve for the unknowns . An exact expression, however, is difficult to obtain. Instead, we assume a polynomial shape with unknown parameters (), namely, The unknown parameters are estimated by minimizing the squared error, , between the observed, , and the estimated observations , where
Once the parameters are estimated, we use inverse filtering to find (the aortic pressure waveform). In the case that the filter order is two, that is, only two time-varying parameters, , and . In this case, the aortic pressure is estimated as where We note from matrix that the parameters and 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 and . 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 Figure 1, are linear time variant. (2)The system is represented by FIR filters. The filters and are second-order and contain many parameters some of which are time-varying.(3)The FIR filter parameters are varying slowly with time.(4)An OU model is assumed to describe the signals/regressors which represent the processes , and ); see (3.1) and (3.6).(5)The time-varying FIR filter parameters 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 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 waveform is calibrated to absolute pressure by scaling it to have the same mean value as the measured peripheral artery pressures as follows: where is the final absolute (scaled) estimated aortic pressure, that is, × Scale factor,
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 by the following steps for each unknown time-varying coefficient .(i)Initiate an estimate for the unknown parameters .(ii)Use Monte Carlo simulation method to generate Wiener processes .(iii)Use (3.19) to an estimate for .(iv)Use (3.16) to find an estimate for .(v)Calculate the summed squared error of (3.20).(vi)Change the estimate for the unknown parameters 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 matrix in (3.25).(5)Estimate the source signal by (3.21).(6)Calculate the estimated aortic pressure within scale factor by (3.23).
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 Figure 3). Third, we demonstrate the ability of the proposed approach to extract aorta AP waveform from multichannel. All the calculations in the algorithm were performed under 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 where is the estimated value of the pressure at instant “” (see Table 1).
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 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, . It was assumed that the signals/regressors 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, (femoral AP) and (radial AP) that are assumed to be represented by the equations In the conventional method, the order of the filters was assumed to be just two, that is, . 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, , of each has the form of the Ornstein-Uhlenbeck (OU) process where are unknown constants, is a polynomial in time that represents the trend or the base data, and bounces around . Remember that in the conventional method We estimate the unknown coefficients and , through regression analysis. In the proposed method, however, the slowly time-varying coefficients , are estimated by the (3.16) in Figure 5 and assumed .
Figure 4 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 using the Malliavin calculus are shown in Figure 5, where and .
(a)
(b)
(c)
(d)
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 , . Table 2 shows the estimated values of central AP using both methods.
Figure 9 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.
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.
Appendix
A. Derivation of the Estimates [37, 38, 42]
Let the augmented signal be defined as It is assumed that we have slowly varying unknown coefficients , that is; or more precisely The signals, , each has an SDE of the form: which has a solution: and assume that with known and known .
Under the slowly varying-coefficient assumptions we have: Substituting for we get an expression for the SDE of as Rearrange (A.1), we get then, which is reduced to: We need to find an expression for that is function only of and ; more specifically as function of and the parameters of . The first step is to find as function of and or . This will be obtained through the generalized Clark-Ocone formula. To do this, we need to do some mathematical manipulations such as the change of the probability measure. This will enable us to find an expression for as function of and .
Let , , , , , a -dimensional Brownian motion, , is the accompanying filtration, is a -dimensional adapted process and .
Define that is, and define the new probability measure where is the old probability measure and is the new probability measure.
Notice that, for a random variable , where the expected value with respect to the new probability measure Q, and is the expected value with respect to the old probability measure .
Girsanov's theorem states that is a Wiener process with respect to the probability measure . In addition is a martingale with respect to [37, 43]. We choose such that the drift coefficient is function only of . Other options or shapes could also be used. Thus, let Notice that Substitute (A.13), (A.16a), and (A.16b) into (A.11), we get: this has the form where
A.1. The Estimates for the Stochastic Time-Varying Coefficients [30, 34]
The generalized version of the Clark-Ocone formula will now be used to find an estimate for the unknown time-varying coefficients as function of and .
Let applying Ito lemma, we get Or, equivalently, Define the random variable We shall now use the generalized Clark-Ocone formula [30] that represents the random variable as: where is the Malliavin derivative of [34] and it is a vector, is the Malliavin derivative of and it is a diagonal matrix since . Equating (A.22) and (A.23) we get Equating the coefficients of on both sides, we get hence,Remember that , , and . After some manipulations, we get () equations in the unknowns . Using the probability measure , each equation is of the form: Or, using the probability measure where .
In our case, , and . Substituting these values and setting , we get
Equation (A.28) when solved yields an estimate for the unknown time-varying coefficient . This requires a knowledge of the shapes of , , , , and . Since a mathematical model for is available, it is easy to get the desired values.
A.2. An Approximate Expression for
Most of the time, it is difficult to find a closed-form expression for this equation. Instead we shall try to put it as a summation of three terms: () constant, () integration w.r.t. time, and () Integration w.r.t. a Wiener process. This way we will be able, in our case, to find a closed-form solution to this formula. We shall present an approximate solution to the above formula.
Since is in general function of , it could be represented by the SDE In our analysis we are more concerned with . Thus we need an SDE for . This has the form Let Using the generalized Clark-Ocone formula for the vector case and the fact that “,” one obtains an expression for the random variable as: where, from the generalized Clark-Ocone formula, In (A.32), the term “” could be approximated as a constant value. Thus, (A.32) or the random variable could be approximated as An approximate SDE for could thus be obtained as Since Then Substituting the different expressions for and for and equating them, we get that is,
Notice that the drift and the diffusion coefficients of are related to each other. This came from the approximate solution of the CO formula.
Without an expression for , we would not be able to solve for the unknowns . Instead of an exact expression, which is difficult to obtain, we shall assume a polynomial shape with unknown parameters , , and , namely, The unknown parameters are estimated by minimizing the squared error, , between the observed, , and estimated observations , where
Abbreviations
AIC: | Akaike information criterion |
AP: | Arterial pressure |
BP: | Blood pressure |
CO: | Clark-Ocone |
CV: | Cardiovascular |
FIR: | Finite impulse response |
OLS: | Ordinary least square |
PAP: | Peripheral arterial pressure |
SDE: | Stochastic differential equation |
SNR: | Signal-to-noise ratio |
SNRE: | Signal-to-noise ratio of the estimate |
TPR: | Total peripheral resistance. |