Sparse models and recursive computations for determining arterial dynamics

Biomedical Signal Processing and Control 38 (2017) 9–21

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Sparse models and recursive computations for determining arterial dynamics Thendral Ganesh a , Jayaraj Joseph b , Bharath Bhikkaji a,∗ , Mohanasankar Sivaprakasam a,b a b

Indian Institute of Technology, Madras, India Healthcare Technology Innovation Centre, Chennai, India

a r t i c l e

i n f o

Article history: Received 7 June 2016 Received in revised form 8 December 2016 Accepted 22 February 2017 Keywords: Carotid artery Ultrasound Arterial stiffness Vascular compliance LTI system Sparsity

a b s t r a c t Arteries expand and contract in every cardiac cycle. Arteries of a healthy individual are elastic. Increased arterial stiffness is an established marker of the vascular health. An estimate of this vascular stiffness may be obtained by measuring the diameter of the Common Carotid Artery (CCA) in each cardiac cycle. This is typically done using image based systems. ARTSENS1 is a portable, image free, ultrasound device for evaluating the stiffness of the CCA. ARTSENS emits a sequence of ultrasound pulses and records the reflected echoes. These echoes are then used to identify the CCA and estimate its diameter, and thereby evaluate the arterial stiffness. This paper deals with development of algorithms for determining the echoes due to the CCA and the estimation of its diameter. Here, the propagation path of each ultrasound pulse is modeled as an FIR filter considering the Gaussian modulated sine (GMS) pulse as the input and its reflections from the walls of the artery and other anatomical structures as the output. The impulse response of the FIR filter is sparse as its output has only few significant echoes. The echoes are reconstructed using the estimated filter coefficients and observed that the reconstructed signal is noise free. This results in the reliable tracking of the artery walls and evaluating its lumen (inner) diameter. The filter coefficients (impulse response) are first determined using Matching Pursuit (MP) algorithms. Additionally, the MP algorithms are made recursive to enable online filtering of the data. The inner diameter of the CCA was calculated for twenty seven subjects using the reconstructed (filtered) data. The estimated diameters were compared with diameters obtained from a B-mode imaging system and was found to be in close match. Furthermore, it is found that for a subject, only the non-zero impulse responses and their sample numbers need to be stored to recover the filtered echoes. Leading to a significant data compression. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Coronary arteries (CoA) supply pure blood (oxygen-rich blood) from the aorta (largest artery that supplies oxygenated blood to the entire body) to the heart muscles. The inner walls of CoA are lined by a thin layer of endothelial cells. These cells preserve the CoA flexible and toned for a smooth flow of blood. Endothelial cells get damaged due to aging, high blood pressure and high cholesterol levels. Damage to the endothelial layer results in the accumulation of degenerative materials such as calcium, lipids, fat and other cellular waste on the inner walls of the artery, thereby narrowing the arterial diameter (lumen, see, Fig. 1) and stiffening the arterial walls.

This occludes the blood flow and diminishes the supply of oxygen and nutrients to the heart muscle causing a Coronary Heart Disease (CHD). Measurement of the arterial stiffness can be used to predict the CHD, [1]. Lumen diameter (along with systolic and diastolic blood pressure) reveals the stiffness ˇ2 of the coronary arteries, [2,3]. Coronary arteries are not available to non-invasive diagnostic tools (such as ultrasound), [4]. It has been noted that stiffness of the Common Carotid Artery (CCA) is a good indicator of the CHD, [5]. CCA runs along the neck region and supplies pure blood to the brain. And significantly, it is accessible for non-invasive diagnos-

p

2

∗ Corresponding author. 1 ARTSENS was developed in Healthcare Technology Innovation Centre. http://dx.doi.org/10.1016/j.bspc.2017.02.010 1746-8094/© 2017 Elsevier Ltd. All rights reserved.

ˇ=

log p s d

Ds −Dd Dd

, where ps , pd are the systolic and diastolic blood preassure

respectively,Ds and Dd are the systolic and diastolic diameteres of the artery.

10

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

Fig. 1. Layers of CCA.

tics. This paper deals with the estimation of the lumen diameter of the CCA. ARTSENS is an image free non-invasive hand held device that can be used for determining the stiffness of the CCA, [6,7]. The main purpose of ARTSENS is to screen large scale populations for early detection and prevention of cardio vascular diseases, [8]. When screening large scale populations the device used must be easily usable by a non-expert operator, portable and must provide results in very short time. Commercially available instruments for measuring CCA stiffness (such as Aloka e-tracking system) are expensive and require a licensed sonographer for operation, [9,10]. The complexity of the computations involved in generating the images and estimation of the arterial diameter from the ultrasound data (Radio Frequency ultrasound signals) makes their usage difficult when screening large populations. Furthermore, imaging systems are subject to the rules and regulations of the countries in use (for example, Preconception and Prenatal diagnostic techniques act (India)) due to their association with fetal sex determination and abortion. Their use is limited to licensed practitioners. While ARTSENS being image free and portable with automated algorithms for wall detection, tracking and stiffness measurements makes it amenable for health care workers to perform quick vascular health evaluation in field. This enables for early detection of at-risk subjects for timely intervention. The carotid Initima Media Thickness (IMT), aortic Pulse Wave Velocity (PWV) and CCA stiffness are established measurements for evaluating vascular health and risk, [11,12]. IMT measurements need a high fidelity imaging system [11,12] and a trained operator. Clinically relevant aortic PWV estimation requires measurement of carotid-femoral PWV [11,12], which is impractical in large scale screening. In many social contexts measurement of carotid-femoral PWV could be a very sensitive issue. However, CCA stiffness measurements could be performed either using ultrasound imaging systems or image free systems [13,14]. Most literature on measuring these indicators can be slotted into two classes, the ones using B-Mode images to measure the indicators, [9,15–18], and the others using Radio Frequency ultrasound signals (RF Signals) directly acquired from the ultrasound sensors (without generating an image), [13,14]. In many cases, the IMT and the lumen diameter of the CCA are estimated with manual intervention, [19,20]. B-Mode systems allow for marking the Intima, Media and Adventitia, see Fig. 1, manually on the images. The IMT and Lumen diameter are estimated by determining the distance between the marked points. In [15] an automated scheme without manual intervention is presented for determining the IMT. This method, based on edge detection schemes, though robust cannot be used for large scale screening due to the use of B-Mode imaging system and the computational complexities involved in generating the image and determining the IMT. In [16], an expert (physician) marks several points manually on the Intima, Media and Adventitia in an ultrasound image. Curves are fit for these marked points. The distance between two diametric points in parallel curves denote the dynamic Intima-Media thickness and lumen diameter. In [17], polynomials were fit to a sequence of images to determine the

Fig. 2. Data acquisition.

shape of the carotid artery, and the distension wave form was estimated from these curve fits. In [18], a sequence of B-mode images are considered and regions of interest are manually marked. Using multiresoultion analysis the contours of CCA are obtained from the marked regions. In [13], Radio Frequency (RF) ultrasound signals were used to assess the IMT, Lumen diameter and distension waveforms. However, RF ultrasound signals were acquired from an echo scanner (image based system) using a dedicated data acquisition system. The imaging system is used to position the ultrasound sensors to obtain RF signals that have significant reflections from the walls of the CCA. This is an important step as mis-positioning would lead to high noise signals with very-little contributions from the walls of CCA. The RF signals acquired were low-pass filtered to remove noise. An adaptive thresholding scheme was used to determine the intima, media and the lumen diameter. In [14], as in [13], RF ultrasound signals were acquired from an image based system using a dedicated data acquisition system. And positioning of the sensor were based on the imaging system. The acquired signal was filtered with an adaptive band-pass filter. Using an adaptive thresholding scheme the diameter of the lumen was estimated. The adaptive thresholding scheme presented in both [14,13] were both similar with a contour tracking around the peaks and an exponential decay between two peaks. The rate of exponential decay being changed adaptively depending on the size of the last peak encountered. Both filtering and wall detection in [14,13] involves making a certain educated guess on deciding the parameters involved. In summary, the three primary tasks, when using systems that acquire RF signals, are (1) Positioning of the sensor to record signals with significant contributions from the walls of CCA. (2) Filter the acquired signals to remove clutter. (3) Detect the walls of CCA and determine the lumen diameter from the filtered signal. Furthermore, all the three steps have to be done in real time. Unlike [13,14], ARTSENS is independent of any imaging system. Positioning of the ultrasound, Fig. 2, sensor to obtain signals with significant contributions from the CCA is explained in detail in Section 2. Here, a sparse signal processing approach is taken to filter the signal. A train of Gaussian Modulated Sinusoidal (GMS) pulses are fired at the CCA, Fig. 3. The reflections (echoes) due to each GMS pulse in the pulse train is acquired and recorded. These echoes are

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

11

Fig. 5. Positioning of the ultrasound sensor.

Fig. 3. (a) Train of M GMS pulses. (b) Reflection due to M = 10th GMS pulse. (c) Reflection due to M = 150th GMS pulse. (d) Reflection due to M = 250th GMS pulse.

2. Hardware specification and data acquisition Noise (w) Input (u)

LTI system (Propagation path)

h

+

Output (y)

+

Fig. 4. Model of the propagation path.

due to reflections form CCA, anatomic structures in the vicinity of CCA and other noise sources. The propagation path of each GMS pulse is interpreted as an LTI system, Fig. 4. In other words, the recorded echo is interpreted as an output of an LTI system (FIR -filter), having sparse coefficients (most filter coefficients being zero), with the GMS pulse as the input. Details on estimation of the filter coefficients, both in real time and offline, are presented in Section 3. Once the coefficients are determined the filtered signal can be obtained by convolving the GMS with the filter coefficients. This not only denoises the signal but also leads to a significant compression in the data that needs to be stored. Here, filtering of the signals does not involve any heuristic parameter choices. The filtered signals, which has only the significant echoes, is thresholded and stacked up in a matrix form with each row corresponding to a reflection of a single ultrasound pulse. An algorithms is then presented to detect walls of the CCA and estimate the diameter and the distension waveform in real time. This paper is constructed as follows. In Section 2 a brief description of ARTSENS and positioning of the ultrasound sensor is presented. In Section 3, modelling the propagation path as an LTI system is described. Estimating its filter coefficients using Basis Pursuit Denoising (BPD) is explained in Section 3.1. The computationally simpler Matching Pursuit (MP) is used to estimate filter coefficients in Section 3.2. Also MP is made recursive in Section 3.2. In Section 4 an algorithm is devised for locating the near Wall (also known as Proximal Wall (PW)) and the far Wall (also known as Distal Wall (D)) of the CCA, and estimate its lumen diameter. The walls of the CCA are detected and the inner diameter is estimated using the filtered signals. In Section 5, Hit Rate of the algorithm is discussed in detail. The paper concludes in Section 6.

In this section a concise summary of ARTSENS is provided along with details on the positioning of ultrasound sensor to obtain meaningful signals. The ultrasound transducer is the primary component of the ARTSENS. It emits GMS pulses with center frequency 5 MHz and receives ultrasonic signals. Here, the transducer records the reflected echoes. A captured echo is first filtered by a high pass filter with a cut off frequency of 2.56 MHz and then amplified with a gain of 40 dB. Further, the echoes are digitized with a sampling frequency of 100 MHz. A detailed description of the design can be found in [6,21]. The ultrasound transducer used in ARTSENS has a maximum sensitivity along its central axis. The beam spread of this transducer is just ˛ = 1.3◦ . That is echoes that are not within a cone, with its axis along the transducers central axis, and of solid angle ˛ = 1.3◦ would not be captured. Assuming the CCA to be a cylindrical tube, if the transducers axis aligns with the diameter of the CCA then a maximum magnitude reflection is captured. If the alignment does not happen, i.e., if the axis aligns with a chord instead of a diameter, then the echoes, being the normal to the point of incidence, would not be captured due to the low sensitivity. Thus all echoes recorded by ARTSENS are due to echoes when the probe is aligned with the CCAs diameter, see Fig. 5. 3. Modeling In order to track the walls of CCA, 800 Gaussian Modulated Sinusoidal pulses (GMS), Fig. 3 were sent with a time interval of 10 ms between pulses. The reflections due to each sinusoidal pulse were recorded. Acquired signals for M = 10, 150 and 250 are shown in Fig. 3b–d respectively. As seen in these figures, there are only three significant echoes per reflection of which the two adjacent echoes represent the echoes from the Near Wall and Far Wall of CCA. Fig. 1 shows the inner layers of CCA. In Fig. 3b–d, either the first and second or the second and third echoes represent reflections from PW and DW respectively. For a single pulse, the propagation path can be modeled as a Linear Time Invariant (LTI) system. Here the GMS pulse is taken as the input and the resultant echoes as the output, Fig. 4. Mathematically stated, y(i) (n)  u(n)  h(i) (n) + w(i) (n), n = 1, 2, . . ., N,

(1)

12

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

where u(n) represents the ith GMS pulse, h(i) (n) and y(i) (n) represent the filter coefficients and reflection that corresponds to the ith pulse respectively, w(n) represents the noise added during the measurement and  represents the linear convolution operator. The superscript (i) is not attached to the input, u(n) as all GMS pulses are same. Note (1) can be rewritten as, (i)

y(i)  Uh

+ w(i) , i = 1, 2, . . ., M

u(1)

⎢ u(2) ⎢ U ⎢. ⎣ .. u(N)

0

···

0

0

u(1)

0

···

0

.. .

.. .

.. .

.. .

u(N − 1)

u(N − 2)

···

u(1)

y(i) = [ y(1) y(2) h

= [ h(1) h(2)

Amplitude

⎤ ⎥ ⎥ ⎥ ⎦ N×N

···

−0.1

···

w(2) · · ·

Sample 2000No.(N)

1900

Fig. 6. Channel coefficient h(1) obtained from (3).

algorithm detects the walls of the CCA. The HR is calculated for all 27 subjects. All the algorithms mentioned in paper are coded in MATLAB and downloaded using SIMULATION in a standard laptop.3 The following subsections explains the BPD and MP algorithms in detail.



y(N) ] , 3.1. Basis pursuit denoising 

h(N) ]

A popular method to solve for h(i) from (2) is to minimize the cost function,

and w(i) = [ w(1)

0

−0.2 1 1000

is the N × N convolution matrix,

(i)

0.1

(2)

where

⎡

0.2



w(N) ]

are N × 1 vectors. In general h(i) , (2), can be estimated using least squares. Note, as there are only few significant echoes in the reflected signal, h(i) will be sparse. However, least squares does not result in a sparse h(i) , [22]. In order to obtain a sparse vector, BPD is normally used, [23,24]. In an earlier paper, [25], the authors had used BPD for estimating h(i) . Solving BPD involves the use of search algorithms that are computationally intensive. In [25], BPD was solved using CVX a matlab toolbox, [26]. A practitioner or an user would like to observe the filtered signals almost immediately after placing the transducer on the neck. When screening large populations, an estimate of the inner diameter is desired while attending to a patient. This is not possible in the case of BPD. Orthogonal Matching Pursuit (OMP) is a computationally efficient tool to solve the sparse signal recovery problem. In [27] it was shown that OMP was as good as BPD in recovering a sparse signal. In this paper, the OMP and a simplification of the OMP known as Matching Pursuit (MP) are used to estimate the filter coefficients. And the results are compared. The MP algorithm used is similar to the one presented in [28]. OMP involves a matrix inversion at every iteration, while, MP has just a single matrix inversion at the end. In addition, the MP algorithm is made recursive using Recursive Least Squares (RLS). This avoids any matrix inversion. And makes their implementation simpler. The estimates of h(i) generated using OMP, MP and Recursive Matching Pursuits (RMP) are convolved with the input to obtain the filtered data. Finally, an algorithm is devised for detecting the walls of the CCA from the noised signal. The above mentioned algorithms, in particular the recursive algorithms, are computationally simple and easy to implement. These enable a user to observe filtered data and the diameter waveform when attending to a patient. In this paper, the ultrasonic data is obtained from 27 subjects using ARTSENS with M being equal to 800, (2). They are used to validate the above mentioned algorithms. For the sake of brevity, plots comparing the reconstructed signal (filtered signal) is with the data y(i) is presented only for i = 1 and M = 1. The inner diameter is calculated for all the 27 subjects and compared with inner diameter obtained from a B-mode imaging system. Also, a parameter called as Hit Rate (HR) which evaluates the ability of the algorithm to detect the walls of CCA is presented. A hit is when, the wall detection

(i) 2

arg minˇy(i) − Uh 2 + ˛h(i) 2 ,

(3)

h(i)

where ||·||2 represents the 2-norm of a vector, ˛ and ˇ are the parameters that control the trade-off between the residual error and the 2-norm of h(i) . It has been shown that the minimizer of (3) is h(i) = [(U ∗ U + (˛/ˇ) I)

−1

U ∗ y],

(4) h(i) ,

elesee [22]. As the minimization in (3) is over the 2-norm of ments with larger magnitudes dominate over the smaller ones. As a result minimization happens over the larger magnitude components with little effect on the smaller values. Therefore, the solution h(i) , (4) would not be a sparse vector. In Fig. 6, h(1) obtained from least squares (4) is shown. Note that (1) h is not sparse. In order to overcome this problem, cost function (3) is replaced by (i)

arg minˇy(i) − Uh 22 + ˛h(i) 1 .

(5)

h(i)

Here ||·||1 denotes the 1-norm of a vector, denotes the 1-norm. Cost function (5) is referred   to as the BPD. Unlike (3), here just the absolute values h(i) (m) are pruned. Therefore, elements with less significant magnitudes are also minimized (made to zero), [22]. This results in a sparse vector with a few significant non-zero values. 3.2. Matching pursuit algorithms Unlike least squares there exists no analytic expression (4), for the minimize of BPD, [22]. However, the cost function (5) is convex in h(i) . Hence, it has a unique global minimizer. The global minimizer is numerically determined using various search algorithms [29–33]. In [25], BPD problem was solved using the CVX toolbox, [26]. These search algorithms are computationally intensive and would require a large number of iterations to reach the minima. Matching Pursuit algorithms are fast and efficient algorithms to recover a sparse signal [27].

3 Dell,Processor:Intel(R)Core(TM)[email protected]. ory:8GBRAM.64bitoperatingsystem.

Mem-

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

3.2.1. Orthogonal matching pursuit algorithm Denote the columns of U as u(1) , u(2) , . . ., u(N) . Since h(i) is a sparse vector, it would at most have k non-zero coefficients with k  N. That is, the data vector y(i) is a linear combination of only k columns of U. Orthogonal Matching Pursuit (OMP) algorithm finds these k columns of U that significantly contributes in forming the data vector y(i) . For the sake of completion, the steps involved in OMP is provided below. 1. At the first iteration t = 1, a residual vector r0 is set to r0 = y(i) and an index set 0 is equated to the null set. 2. Compute the magnitudes of the inner products

    < r0 , u(m) > = r0 u(m)  m = 1, 2, . . ., N

(6)



1  arg max < r0 , u(m) > .

(7)

m=1,...,N

3. Update the index set to 1 = 0 4. A least squares problem (i)



{1 } and set U1

(i)

h1 = arg miny(i) − U1 h1 2

= [u(1 ) ] .

(8)

(i) h 1

is solved to obtain the estimate of h(i) due to U1 . 5. The new estimate of the data vector y(i) and the new residual (i) (i) (i) are calculated by setting yˆ 1 = U1 h1 and r1 = y(i) − yˆ 1 respectively. 6. At the second iteration t = 2, the new residual r1 is used to determine

  2 = arg max < r1 , u(m) >

Method

Time(s) – MATLAB (standard laptop)

BPD MP MPSR

11.1923 0.2519 0.351

3.2.2. Matching pursuit algorithm 1. Initialize the residual r0 = y(i) , the index set ␹ iteration counter t = 1. 2. Find the index t such that



7. The index set is updated to 2 = 1 matrix is set to U2 = [U1 , u(2 ) ]. 8. As before, a least squares estimate





(ut ) rt−1  (ut ) ut

u(t ) .

5. Increment t and go to step 2. 6. Repeat this process k times such that rk+1 is insignificant.Let Uk denote N × k matrix with the k chosen columns. The k × 1 coefficient (i) vector hk is taken to be minimize of J(h) = y(i) − Uk h2

(14)

over all k × 1 vectors h. The coefficient vector ting



(i)

h (m) =

0

if m = / t , t = 1, 2, . . ., k

(i) hk (m)

if m = t , t = 1, 2, . . ., k

h(i)

is formed by set-

.

(15)

Note (4) will give the minimize of (14). The MP algorithm has two elements:

(10) The number of multiplications involved in MP is,

of h(i) due to the updated column matrix U2 is obtained. (i) (i) 9. A new estimate yˆ 2 = U2 h2 of the data vector and the residual (i)

r2 = y(i) − yˆ 2 are computed. 10. This process is repeated until t = k, where k is such that residual rk+1 is numerically insignificant. h (i)



0

if = / t , t = 1, 2, . . ., k

(i) hk (m)

if m = t , t = 1, 2, . . ., k

.

(11)

In the above mentioned OMP algorithm a least squares problem involving a matrix inversion is solved for each iteration. This is computationally expensive. It can be shown that the number of multiplications involved in OMP is, 2

(13)

1. Determining the columns of U that contribute to the output. 2. Determination of the non-zero elements of h(i) .

(i)

kN + (N + c)



4. Set the residual rt = rt−1 − rˆt−1 , where rˆt−1 =

(i) h 2

h (m) =

=  and the

the matrix of chosen columns: 3. Augment  the index set and t and Ut = [Ut−1 , u(1 ) ] . t = t−1

{2 } and the column

h2 = arg miny(i) − U2 h2 2

(i)

0

t = arg max < rt−l , u(m) > .

(9)

m=1,...,N

(i)

Table 1 Time taken by the algorithms to Denoise one Frame.

m=1,···,N

and determine the maximize, i. e ., set



13



k(k + 1) 2

2

 +N

k(k + 1) 2

 +

k(k + 1)(2k + 1) 6

(12)

where k is the number of iterations, N is the length of data and c is constant

due  to the inversion of matrix (inversion of an k × k matrix is O k3 ). In all the plots presented here the data length was taken to be N = 1900 and the number of iterations were found to be k = 18. Setting c = 1, the number of multiplications involved in OMP is 120894150. The constant c depends on the entry of the matrix. Refer to the Appendix for a derivation of (12). In the next subsection the OMP algorithm is modified to avoid the matrix inversion at every iteration.

2

3

(kN + 3kN + k2 N + k2 + ck ).

(16)

A derivation of the above is presented in the Appendix. For N = 1900 and k = 18. Setting c = 1, the number of multiplications involved in MP is 67504356. The time taken to obtain the filtered signal is tabulated in Table 1. The time was estimated using the tic and toc commands in Matlab. In order to keep the exposition simple, first a semi-recursive algorithm, wherein, only the least squares (14) is made recursive. In the semi-recursive algorithm it is assumed that the columns of U that contribute to the output are known before having to solve the least squares. Finally, in the fully recursive scheme the choice of the columns are also made recursive. 3.2.3. MP semi-recursive (MPSR) Note solving (14) requires inversion of the matrix Uk Uk . In the following, it is assumed that the matrix with k columns that combine to give coefficient vector h(i) are known. Without loss of generality they are taken to be u(1) , u(2) , . . ., u(k) . Here, the com(i) putation of hk from (14) is made recursive using Recursive Least Squares (RLS), [34]. At iteration t, t ≤ N, let



Uk,t−1  u(1) t−1 and (i)



yt−1  y(i) (1) (i)

(2)

ut−1

y(i) (2)

···

···

(k)

ut−1



(t−1)×k

y(i) (t − 1)

where yt−1 is a (t − 1) × 1 vector.



(17)

(18)

14

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

1

Amplitude

Amplitude

1

0

−1 BPD using CVX Toolbox MP Semi−Recursive

−2 1

0

−1 OMP MP Semi−Recursive

−2 1000 1

No. (N) 2500 1500 Sample 2000

1900

1900

Sample No.(N)

(a)

(a)

1

Amplitude

Amplitude

1

0

0

−1

−1

−21

OMP MP Semi−Recursive

−21 0

BPD using CVX Toolbox MP Semi−Recursive

500

1000 Sample No. (N)1500

1900

Sample No.(N)

1900

(b)

(b)

(1)

Fig. 8. (a) Channel coefficient h(1) . (b) Noised signal yd . (1)

Fig. 7. (a) Channel coefficient h(1) . (b) Denoised signal yd .

1

denotes the mth column In ((17)) t − 1 elements, i.e.,



(m)

ut−1  u(m) (1)

u(m) (2)

···

u(m) (t − 1)

u(m)



with only the first

.

Amplitude

(m) ut−1

(19)

0

−1

Alternatively, Uk,t−1 denotes Uk truncated by retaining only the first (i) t − 1 rows. Similarly, yt−1

denotes truncated y(i)

elements. Let



(i)

ht−1  h(i) (1) t−1

(i)

ht−1 (2)

···

(i)

ht−1 (k)

with only first t − 1



Acquired Echo Denoised Signal

−2 0 1

Sample No.(N) 1500 1000

500

1900 (1)

k×1

(20)

Fig. 9. The echo from CCA (red) using ARTSENS and reconstructed signal yd using MP Semi-Recursive scheme(black). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

denotes the minimize of the cost function (i)

Jt−1 (h) = yt−1 − Ut−1 h2

(21)

over all k × 1 vectors h. Augmenting the latest values, set



Uk,t  and (i) yt





Uk,t−1



(22)

ut

(i)

yt−1

 (23)

y(i) (t)

where ut denotes the tth row of Uk and y(i) (t) is the tth element of y(i) . The recursions, (i)

(Pt−1 u t )

(i)

ht = ht−1 +

(1 + ut Pt−1 u t )

(y(i) (t) − u t h(t−1) )

 Pt = Pt−1 − (Pt−1 u t (1 + ut Pt−1 ut )

(i)

−1

(24)

ut Pt−1 )

(i)

generate the estimate ht , the minimize of Jt (h), ((21)) at iteration (i) h−1

t. Initialize = 0 and P−1 = I. The number of multiplications involved in MPSR is 2

2

(kN + 4kN + k N).

from both of these methods are sparse. Fig. 7b shows the recon(1) structed signal yd obtained by using BPD solved by CVX toolbox and MP Semi-Recursive scheme. It is apparent that the noised sig(1) nal yd obtained from both these methods are of a reasonably good match. Fig. 8a compares the h(1) obtained from OMP and MP SemiRecursive scheme after N iterations, where N = 1900 is the total data (1) length. Similarly, Fig. 8b compares the noised signal yd obtained from OMP and MP Semi-Recursive scheme. Fig. 9 shows the signal y(1) acquired using ARTSENS and noised (1) signal yd . It is apparent from the figure that MP semi recursive scheme retains the significant echoes and suppresses the noise in the signal y(1) . 3.2.4. MP fully recursive In this sub section, the process of recursively finding the matrix Uk is presented. (i) 1. At the first iteration t = 1, a residual vector r0 is set to r0 = yp where p is an integer much smaller than N and reasonably greater than one. Set (i)

(25)



yp  y(i) (1)

y(i) (2)

···

y(i) (p)



(26)

p×1 (i)

For N = 1900 and k = 18. The number of multiplications involved in MPSR is 67532400. As there are no matrix inversions, the constant c is not involved in (25). Hence the number of computations doesn’t depend on data. The time taken is tabulated in Table 1. Thus, all these methods are easily implementable for real time applications and accurate estimation of inner diameter as only multiplications are involved. Fig. 7a shows the sparse h(1) obtained by using BPD solved by CVX toolbox and MP Semi-Recursive scheme. Note h(1) obtained

and equate the index set 0 to null set. Note yp contains the first p elements of y(i) . 2. Compute the magnitudes of the inner products

       (m)  (m) < rt , up > = rt up  m = 1, 2, · · ·, N where (m)

up



 u(m) (1)

u(m) (2)

···

u(m) (p)

 p×1

,

(27)

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

1 Amplitude

Amplitude

1

0

−1 BPD using CVX Toolbox MP Fully Recursive

−2 01

(1)

(m)

m=1,...,N

1900

1

 

> .



4. Set the residual rt = rt−1 − rˆt−1 , where rˆt−1 =

   (up t ) up t

up(t ) .

1

(i)

J(h) = yp − Uk1 ,p hk1 ,p 2 .

(28)

Note (28) can be solved using the recursions (24). For the next iteration, next block of elements in y(i) and next block of rows in U are added in (26) and (27) respectively. i.e.,

⎤

⎢ (i) ⎥ ⎢ y (p + 1) ⎥ ⎢ ⎥ (i) yp+c  ⎢ ⎥ ⎢ .. ⎥ ⎣. ⎦

0

−1



(up t ) rt−1

5. Increment t and go to step 2. 6. Repeat this process k1 times such that rk1 +1 is insignificant. Note by construction k1 < k, where k is the sparsity of h(i) . 7. Let Uk1 ,p denote the p × k1 matrix with k1 chosen columns considering only the first p elements in y(i) and first p rows of U. (i) The k1 × 1 coefficient vector hk ,p is taken to be minimize of

(29)

(1)

MPSR and RMP involves only matrix multiplications. Hence, they are easily implementable for real time applications. (i) The denoised signal yd is reconstructed using the channel coef-

ficient h(i) that was estimated by, (i)

(i)

yd = Uh .

where c is the block size. (i) Then set t = 1 and r0 = yp+c for m = 1, 2, . . ., N. Now repeat the steps 1 to 7 to get Uk2 ,p+c and new hk2 ,p+c . Note the recursions in (24) will run p times at the first iteration, p + c times at the second iteration, p + 2c times in the third iteration and so on. To illustrate the scheme in 3.2.4, p was set to p = 600 and c was set to c = 650. Fig. 10 shows the reconstructed signal when the first block of y(1) (i.e., first p elements of y(1) ) and U (first p rows of U) were used and k1 was found to be k1 = 6. Similarly, Fig. 11 shows the reconstructed signal when the second block of y(1) (i.e., p + 1th to p + cth elements of y(1) ) and U (p + 1th row to p + cth row) were augmented and k2 was found to be k2 = 12. Fig. 12 shows the reconstructed signal when the third block of y(1) (i.e., p + c + 1th to p + 2cth elements of y(1) ) and U (i.e., p + c + 1th row to p + 2cth row of U) were augmented and k3 was found to be k3 = 18. Note, p + 2c = N.

(31)

4. Wall detection algorithm In order to find the inner diameter of CCA, walls of the artery must be located first. This section presents an algorithm to locate the artery walls and then to find its inner diameter. Reconstruct the noised signal from (31) for all i = 1, 2, . . ., 800. A binary signal, b

up+c

1900

Fig. 12. Reconstructed signal yd using MP Fully Recursive method when third block

(i)

(30)

No. (N) 1500 500 Sample1000

of y(1) and U were augmented and k3 = 18.



⎤

BPD using CVX Tolbox MP Fully Recursive

−2 01

y(i) (p + c) Uk,p

1000No. (N)1500 Sample

(1)

3. Augment matrix of chosen columns:  the index set and the t) t = t−1 t and Ut,p = [Ut−1,p , u( p ].

⎢u ⎥ ⎢ p+1 ⎥ ⎥ Uk,p+c  ⎢ ⎢. ⎥ ⎣ .. ⎦

500

Fig. 11. Reconstructed signal yd using MP Fully Recursive method when second

Amplitude

 

t = arg max < rt−1 , up

⎡

BPD using CVX Toolbox

block of y(1) and U were augmented k2 = 12.

is the mth column of the matrix U with first p elements.Determine the maximize,

(i)

−1 MP Fully Recursive

Fig. 10. Reconstructed signal yd using MP Fully Recursive method for first block of

yp

0

−21 0

No. (N) 1500 1900 500 Sample1000

y(1) and U and k1 = 6.

⎡

15

=

(i)

(i)

1,

yd (n) ≥ max (yd )/100, i = 1, 2, . . ., 800

0,

otherwise (i)

is constructed by thresholding each yd .

A matrix D is constructed by stacking the binary signal b(i) for i = 1, 2, . . ., 800 as rows, 

D  [b(1) (n) b(2) (n) · · · b(800) (n)] ,

(32)

where  denotes the transpose. Matrix D, (32), is plotted as a 2D contour plot, see Fig. 13a. In this plot the rows of D(i, n) for i = 1, 2, . . ., 200 is plotted along the y axis with the x-axis being the sample number (or the column index n) of the matrix D respectively. Pulsation of the artery walls in opposite directions is evident from the above plot. The blank white space indicates the zeros in matrix D and the red colored blocks indicate the ones in the matrix D. From Fig. 13a, it can be observed that the adjacent blocks R1 and R2 pulsate in opposite direction where as the blocks R2 and R3 pulsate in the same direction, with R3 having very limited or no pulsation. This is because the walls of CCA pulsate in opposite direction as they expand (diastole) and contract (systole) in the process of a cardiac cycle. Thus it is obvious that the adjacent blocks R1 and R2 represent the reflections from the walls of CCA.

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

No. of Frames

16

and Eq. (35) reduces to,

200 150 100 50

R1

R2

k1

0

500

e1(B2 ) = arg max{x} + k1

R3

k2

1000

x

e1(B2 ) = 425 + 600 = 1025

k3

1500

Thus s1(B2 ) = 946 and e1(B2 ) = 1025 for row m = 1 in the block B2 . This procedure is repeated for rows m = 2, 3, . . ., 800 and block B3 . Using (34) and (35), a 1 × 6 matrix Xm is generated with its mth row being,

200 150 100 50



(B1 ) Xm = sm

01

No. (N) 1500 500 Sample 1000

1900

Using this fact, an algorithm is built to automate the detection of the walls of CCA. As it is apparent in the plot that each blocks of ones are partitioned by huge columns of zeros (>100), matrix D can be markedly partitioned into blocks as (33)

where B1 , B2 and B3 are 800 × k1 , 800 × k2 and 800 × k3 matrices respectively and k1 + k2 + k3 = N. Here k1 = 600, k2 = 600 and k3 = 700. These numbers are chosen based on the fact that there are large number of zeros between the reflections. Let (Bp ) sm = arg min{x} +

p−1 

x Bp (m,x)=1

(Bp ) em

= arg max{x} + x Bp (m,x)=1

kj , p = 1, 2, 3



X = X1

···

X2

X800



(B2 ) (B ) Xmd = sm − em 1

(35)

s1(B1 ) = arg min{x} + k0 x



(36)

,

(B )

(B )

sm 3 − em 2



,

(37)

(B1 ) − em



(B )

(B )

(38) (B1 )

{s1 1 , s1 1 + 1, · · ·, sF (B1 )

{sF

(B1 )

, sF

(B )

(B1 )

s1 1 > sF

(B )

(B )

and (B1 ) S = {s2(B1 ) , s3(B1 ) , · · · sF−1 }.

(39)

where Q is the set of F consecutive rows of matrix X which gives the column numbers corresponding to where the 1 starts in the matrix D, for the rows m = 1 to F, R is the range vector that either contains the numbers starting from s1(B1 ) to sF(B1 ) or the numbers starting from sF(B1 ) to s1(B1 ) depending on whether s1(B1 ) < sF(B1 ) or s1(B1 ) > sF(B1 ) respectively and S is a subset of Q, that contains all the elements from Q except s1(B1 ) and sF(B1 ) . In this example, Q = {135, 132, 129, 129, 129, 128, 128, 127, 127, 126}, R = {135, 134, 133, 132, 131, 130, 129, 128, 127, 126} and S = {132, 129, 129, 129, 128, 128, 127, 127}. If

then

e1(B1 ) = arg max{x} + k0

direction =

Thus s1(B1 ) = 135 and e1(B1 ) = 196 for row m = 1 in the block B1 . Now for row m = 1 in the block B2 , Eq. (34) reduces to,

(B1 )

− 1, · · ·, s1 1 − 1, s1 1 },

and the (35) reduces to,

e1(B1 ) = 196 + 0 = 196

− 1, sF

(B )

s1 1 < sF

S˜ = {si(B1 ) ∈ S |si(B1 ) ∈ / R} = 

x

(B1 )

},

s1(B1 ) = 135 + 0 = 135

s1(B2 ) = 346 + 600 = 946



Q = {s1(B1 ) , s2(B1 ) , · · · sF(B1 ) },

(Bp ) and e(Bp ) represents the column number of where k0 = 0 . Here sm m the mth row where the one starts to appear and ends in the block Bp respectively. Eq. (34) finds the column number, in the matrix D, where the 1 starts for the mth row and in the Bp th block. Similarly, Eq. (35) finds the column number, in the matrix D, where the 1 ends for the mth row and in the Bp th block. For example, from Fig. 13 it can be see that, for the row m = 1 in the block B1 , the 1 starts at 135th column and ends at 196th column of B1 . Thus Eq. (34) reduces to,

x

(B )

em 3

in (37) will give the no. of for m = 1, 2 · · ·, 800. The term columns between the 1 that ends in the first block B1 and the 1 that starts in the block B2 in the matrix D. Matrix Xmd is constructed to know the distance between the reflections in Fig. 7b (Black dotted waveform). For example, for row m = 1, X1d = [750 350]. At least 2 rows of the matrix X is needed to identify the walls of the artery. Pick F (<20) rows of matrix X. The direction that the ones move along the rows of the block B1 can be identified by defining the sets

(34)

j=0

s1(B2 ) = arg min{x} + k1

(B )

sm 3

(B2 ) sm

R= kj , p = 1, 2, 3

(B )

em 2

is constructed. For example, for row m = 1, X1 = [135 196 946 1025 1375 1443]. A distance matrix Xmd is constructed using matrix Xm with its mth row being,

j=0 p−1 

(B )

sm 2

Stacking up matrix Xm for m = 1, 2 · · ·, 800, a matrix X,

Fig. 13. (a) 2 − D plot of the matrix D(i, n), (32), with i = 1, 2, . . ., 200 of the 800 rows along y axis and the corresponding sample indices n in the x axis. (b) 2 − D of the matrix D once the artery walls are located. It can be observed that R1 and R2 pulsate in the opposite direction R3 showing very limited pulsation along the same direction as R2 .

D800×N = [B1 |B2 |B3 ],

(B )

em 1



(40)

L,

s1 1 − sF

(B )

(B1 )

>0

R,

(B ) s1 1

(B ) − sF 1

<0

.

(41)

where  represents the null set. Eq. (40) finds the elements in the set S and not in set R. This is done to check whether the elements in S lie with in the range s1(B1 ) : sF(B1 ) i.e., 126 : 135. This corresponds to artery wall either moving from contraction to expansion state or from expansion to contraction state. Refer Fig. 14a and 14b. In the example considered, S˜ =  and the direction is found to be L as s1(B1 ) − sF(B1 ) = 135 − 126 > 0. If / R} = / , S˜ = {si(B1 ) ∈ S |si(B1 ) ∈

(42)

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

17

Fig. 14. (a) and (b) Artery wall movement either from expansion to contraction or from contraction to expansion. (c) and (d) Artery wall movement either expansion to contraction to expansion or contraction to expansion to contraction.

define









˜ s(B1 ) > s(B1 ) smax = max si(B1 ) ∈ S| 1 i and

(43)

˜ s(B1 ) < s(B1 ) . smin = min si(B1 ) ∈ S| 1 i

(44)

The direction is given by,



direction =

(B )

(B )

(B )

(B )

R → L,

smax > s1 1 , smin > s1 1

L → R,

smax < s1 1 , smin < s1 1

.

(45)

This corresponds to artery wall either moving from contraction to expansion and again to contraction state or moving from expansion to contraction and again to expansion state. Refer Fig. 14c and d. The same procedure is carried out for blocks B2 and B3 . After finding the direction of each block,locate the two consecutive blocks moving in opposite direction. In this illustration (see Fig. 13) the blocks moving in opposite direction is found to be B1 and B2 . Therefore these blocks correspond to the walls of artery. However, Block B3 corresponds to some other anatomical structure, it can be removed. This process of locating the artery walls can also be made recursive in time for each m. The sample indices of the ones in Xm can be recursively updated for each m as in (34) and (35). Once the indices are found for F < 20, the direction of each block can be found using (38) to (45). The whole process of determining the lumen diameter and distension waveform is summarized as a flow chart in Fig. 15 The lumen diameter of CCA is evaluated as,



Ld = (s(B2 ) − e(B1 ) )c/2fs 1 1

···

(B )

(B )

(sM 2 − eM 1 )c/2fs

T

where c = 1540 m/s speed of sound in tissues and fs = 100 MHz is the sampling frequency. Figs. 16 and 17 shows the plot of the distension waveform estimated using the reconstructed signal from BPD and MP SemiRecursive scheme, respectively, of one subject. In both the methods, the inner diameter shows good repeatability (between 5.9 and 6 mm over the 10 cycles plotted) with a precision of ±0.5%. Also refer to Fig. 18 for a Bland-Altman plot, a commonly used plot in clinical literature, [35]. In the Bland-Altman plot, the differences between the inner diameter from BPD and MPSR are plotted against the mean of the two measurements for all the 27 subjects. The blue ¯ taken line in Fig. 18 represents the mean of the differences (m) over all the 27 subjects. And the red dotted lines represents the ¯ − 1.96sd ) and (m ¯ + 1.96sd ), where sd is the agreement intervals (m standard deviation of the differences, with in which the 95% of the differences are expected to fall for a good match. It can be noted that all the points are within the two read lines suggesting a good match. Fig. 19 shows the comparison of the Inner diameter obtained from a B-mode Imaging system with the inner diameters estimated using BPD and MPSR for all the 27 subjects. It is evident form Fig. 19 that the values obtained form B-mode and proposed method are in a reasonably close match for most values. As above, a BlandAltman plot, Fig. 20, is presented with differences between the

Fig. 15. Summary of the process to determine the CCA walls and estimation of its lumen diameter and distension waveform.

inner diameter obtained from B-Mode Imaging system and MPSR plotted against the mean of the two measurements for all 27 subjects. Note all the 27 points are within the two red lines suggesting a good match between diameter measurements of B-Mode images and MPSR values. 5. Hit Rate As mentioned above, ARTSENS was used to record the ultrasound data from 27 human subjects. Train of 800 ultrasound pulses

HIT RATE(%)

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

Inner diameter (mm)

18

6.3 6.2 6.1 6 5.9

100 80 60 Volunteer 1 Volunteer 2 Volunteer 3

40 0

0

2

4 6 Time (s)

8

10 20 No. of frames

30

10 Fig. 21. Hit Rate (%) vs frames using BPD.

Inner diameter(mm)

Fig. 16. Inner diameter – BPD.

6.3 6.2 6.1 6 5.9 0

2

4 6 Time(s)

8

10

Fig. 22. Hit Rate (%) vs frames using MP semi-recursive scheme.

Fig. 17. Inner diameter – MP Semi-Recursive scheme.

Table 2 Time taken by the algorithms to detect the walls of CCA.

0.240492 1.96s

Method

Time(s) – MAT LAB (standard laptop)

BPD MP MPSR

55.9615 1.2595 1.755

-0.0107407 Mean

0

FOR F = 5 10

-0.26197 -1.96s

-0.5 4

4.5

5

5.5

6

6.5

7

7.5

Frequency

Difference of BPD and MPSR inner diameters

0.5

8

Average of two BPD and MPSR inner diameters

Inner Diameter (mm)

Fig. 18. Bland-Altman plot. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

10

6 4 2 0 0

5

10

15

20

25

60

70 80 Hit Rate (%)

90

100

Fig. 23. Distribution of the Hit Rate (%) for F = 5.

30

Subjects Fig. 19. Comparison of inner diameter. Difference of B-mode and MPSR inner diameter

0 50

B−mode imaging system MPSR BPD

8

5

were sent and its corresponding echoes were captured for each subject. Each reflected echo is considered as a frame. Typically, F(<20) number of frames were enough to identify the CCA walls. Hit Rate (HR) determines the efficiency of the algorithm to detect the artery walls. A hit is counted when the PW and DW of CCA is detected by the algorithm. For every trial, the algorithm uses F consecutive frames, with the first frame being randomly picked from 800 − F frames, to detect the walls of CCA. In order to calculate HR nearly 100 trials were carried out.

2

HR =

1.5 1 0.41832 1.96s

0.5 0

-0.4884 Mean

-0.5 -1

-1.395 -1.96s

-1.5 -2 4.4

4.6

4.8

5

5.2

5.4

5.6

5.8

6

6.2

Average of B-mode and MPSR inner diameters

Fig. 20. Bland-Altman plot. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

No. of Hits × 100. No. of Trials

(46)

HR Vs different number of frames (F), obtained from reconstructed signal using BPD is plotted in the Fig. 21. Similarly, the plot of HR vs different number of frames (F), obtained from reconstructed signal using MP Semi-Recursive scheme is shown in Fig. 22. The algorithm were able to locate the artery walls for more than 80% of the times using just F ≥ 5 frames. In Table 2, the time time taken by the algorithms BPD, MP and MPSR, to estimate filter coefficients for 5 frames and detecting the walls using these 5 frames are tabulated. The Figs. 23–25 show the distribution of HR for F = 5, F = 15 and F = 25 respectively for all 27 subjects. It is evident from the plots that HR improves when the value of F is increased. For F = 5, the

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

FOR F = 15

Frequency

10

5

0 50

60

70 80 Hit Rate (%)

90

100

Fig. 24. Distribution of the Hit Rate (%) for F = 15.

FOR F = 25

Frequency

10 8 6 4 2 0 50

60

70 80 Hit Rate (%)

90

100

Fig. 25. Distribution of the Hit Rate (%) for F = 25.

number of subjects having HR more than 80% is 12. Similarly for F = 15 it is 14 subjects and for F = 25 it is 18 subjects. 6. Conclusion This paper considers the detection of the Common Carotid Artery (CCA) and the estimation of its inner diameter from ultrasound measurements. Blood flow causes the CCA to pulsate. Tracking the inner diameter during pulsation is a good indicator of its stiffness. ARTSENS is a non-invasive image-less technology to track the dynamics of the walls of CCA and estimate the arterial stiffness for screening large populations. ARTSENS has an ultrasound probe, which when kept on the neck of a subject, and triggered, emits a train of GMS pulses. The resulting reflections were recorded by ARTSENS. The reflected signal, corresponding to every single pulse, would contain echoes from the walls of CCA and other anatomical structures in the vicinity. The main difference between the reflected signals would be in the temporal position of the echoes due to the walls of the CCA. The temporal position of the echoes due to other anatomical structures would not change as they remain static. In this paper, each reflected signal is modeled as the output of an FIR filter with the input being the Gaussian modulated sinusoidal pulse. The impulse response of the filter was estimated from the input and the output data. As there are only a few anatomical structures that reflect the Gaussian modulated pulses, the reflected signals and the impulse response tend to be sparse. Reconstructing the signal by involving the estimated impulse response with the pulse results in a noise free data. Basis Pursuit Denoising (BPD) was used to estimate the filter coefficients (impulse response) in [25]. The Hit Rate using BPD filtered data was found to be above 80% which just five frames or echoes. However, solving BPD involves a computationally intensive search procedure. A physician or a user would not get to see a filtered signal “as soon gas’s the ultrasound probe is placed on the neck. This is necessary to adjust the probe position to capture a good response. Furthermore, the diameter waveform cannot be observed when attending to a patient. This is necessary when screening large patients.

19

Orthogonal Matching Pursuit (OMP) is a computationally simpler alternative to estimate the sparse filter coefficients. OMP involves solving a least squares problem (Matrix inversion) at every iteration. Here the OMP algorithm is simplified to Matching Pursuit (MP) that involves solving just a single least squares (a single matrix inversion) at the end. In order to avoid even this, the MP was made recursive (MP semi recursive and MP fully recursive). The recursive versions were then used to the estimate filter coefficients and reconstruct a clean signal online. These recursive algorithms enable observing a filtered signal “immediately” after placing the probe on the neck. They also allow for the user to observe the diameter waveform in quick time when attending to a patient. An algorithm was developed to locate the walls of CCA and then estimates its lumen diameter from the reconstructed signals. Reconstructed data obtained from both BPD and Matching Pursuit Semi-Recursive (MPSR) was used to evaluate the performance of the algorithm. Both resulted in a Hit Rate greater than 80% with just five reflected signals. A stable Tracking of the inner diameter of CCA was achieved for MP, MPSR and RMP. Besides robust detection and tracking, the above approach of modeling the reflected signals as outputs of a filter with only a few non-zero sparse coefficients leads to a significant data compression. As it is adequate to store only the non-zero impulse response and its corresponding sample number (indices). In [36] a study was carried to determine the match between ultrasound measurements (B-Mode) of IMT, lumen diameter of the CCA with the equivalent pathological measurements. The measurements were carried out on 66 terminally ill neurological patients with the ultrasound measurements taken a few days before their death and pathological measurements, by autopsy, after their death. It was concluded therein that the B-mode measurements provided reliable estimates of the lumen diameter of the CCA. From the Bland-Altman plot, Fig. 20, it is clear that the lumen-diameter estimated using algorithms presented here a close match to the corresponding ones obtained using B-Mode images. Thus it can be concluded that the estimates of the distension waveform and lumen-diameter presented here are clinically reliable. Among all the references mentioned in this paper, only [14,13] use RF use signals, the rest use B-Mode images. In [14] results for two subjects were presented. One of 40 years having an average diameter of 7.7 mm and another having a average distension of 500 ␮m, while the other a 70 year old with a distension of 400 ␮m (no diameter values mentioned). In [13] average diameter and the distension of a young healthy subject was presented. The values were 7 mm and 689 ␮m. As the data set presented is limited a direct comparison with the results in this paper is not possible. However, it can be noted that all the values presented were in the expected range. The inter-operator reproducibility and intra-operator variability studies have been performed for ARTSENS [37–39]. The intraoperator and inter-operator variability was found to be less than 15%, which is acceptable and comparable to results reported in other literature [40,41]. ARTSENS was found to provide measures of stiffness with sufficient accuracy and repeatability, as compared to other imaging based systems. ARTSENS also demonstrated capability to detect the age-associated increase in vascular stiffness [38] thereby indicating potential to detect physiological changes in stiffness. By virtue of its automated measurement algorithms which eliminate the need of a trained sonographer, and its inherent portability it is a promising alternative to current systems for vascular stiffness evaluation in large scale screening. An extensive study to quantify the ability of ARTSENS to detect changes in arterial stiffness in subjects with established risk markers (such as low HDL, hypertension, diabetes, central obesity etc) is currently ongoing.

20

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

Appendix A

requires 2N computations per iteration. Therefore, for k iterations 2kN computations. And lastly, the Equation in Section III B. 2

1. Derivation of the expressions (12), (16) and (25)

J(h) = y(i) − Uk h2

Note that BPD, OMP and MPSR attempt to solve the cost function, (i)

argminh(i) ˇy(i) − Uh 22 + ˛h(i) 1 .

(47)

(a) MOP The computational complexity of OMP algorithm is calculated here.

    r0 , u(m)  = r0 u(m)  m = 1, 2, . . ., N

(48)

requires N2 multiplication for every iteration. As the algorithm is found to run for k iterations, total number of multiplication required due to (48) is kN2 . And note in step 4 of Section III B. 1 (i)

h1 = argmin

(i)

(i) 1

h

y(i) − U1 h1 2

(49)

a least squares problem is solved for which the solution is

h1 = U1 U1 (i)

−1

U1 y(i) .

(50)



In (50), to obtain A = U1 U1 matrix, N multiplications are required as U1 is N × 1 matrix at the 1st iteration. Similarly, to obtain B = U1 y(i) matrix, N multiplications are required. Finally to obtain A−1 B matrix, N multiplications are required. Note the OMP algoTherefore, in general for jth iteration, to rithm run for  k iterations. Uj Uj matrix, j3 N multiplications are required. Sim-

obtain Aj =

Uj y(i) ,

ilarly, to obtain Bj = obtain

A−1 Bj , j2 j



jN multiplications are required and to

multiplications are required. And for the matrix

−1

 3

(56)

is a least squares problem whose solution is obtained using

h = Uk Uk

−1

Uk y(i) .

(57)

In (57) to obtain Ak = Uk Uk



matrix, k2 N multiplications are

required. Similarly, to obtain Bk = Uk y(i) matrix, kN multiplications

are required. And finally to obtain A−1 × Bk matrix, k2 multik plications are required. To obtain the inverse of matrix Ak , ck3 multiplications are required computations as it involves matrix inversion. Thus the computational complexity for MP algorithm is given by 2

3

(kN + 3kN + k2 N + k2 + ck ).

(58)

It can be seen from (53) and (58) that the computational complexity is significantly reduced. The complexity of MP

computational  2 algorithm is in the order of O kN . Further the computations of MP algorithm is simplified using Recursive Least Squares (LS). (c) MPS The computational complexity for MPS is calculated here. Eqs. (54) and (55) are used in MPS which requires kN2 and 2kN computations for k iterations respectively. In the Recursive equations in Section III B.3, (i)

(i)

ht = ht−1 + Pt =

(Pt−1 u t )

(1 + ut Pt−1 u t )

(y(i) (t) − u t h(t−1) ) (i)

(59)

 −1 Pt−1 − (Pt−1 u t (1 + ut Pt−1 ut ) ut Pt−1 ),



multiplications are required. Thus for

2 to perform D = Pt−1 u t , k multiplications are required. And to

k iterations, total number of multiplications required is calculated as

perform (ut D) and u t h(t−1) , 2k multiplications are required. Note

inversion Uj Uj

(N + c)

k 

j3 + N

j=1

, O cj

j=k 

j+

j=1

k 

j2

(51)

j=1

Solving (51), (N + c)

 k(k + 1) 2 2

+N

 k(k + 1)  k(k + 1)(2k + 1) +

2

6

2



k(k + 1) 2

2



+N

k(k + 1) 2

 +

k(k + 1)(2k + 1) 6

.

(52)

.

(53)

In order to avoid matrix inversion at every iteration, OMP is modified in such a way that least squares had to be solved once at the end. This is known as Matching Pursuit (MP) algorithm, [28]. It is worth noting that w.r.t to the current case where N = 1900 and k = 18, k3 and k4 is much larger than N. Note Eq. (53) is in the O(kN2 + k4 N). (b) MP The computational complexity of MP algorithm is calculated here. In step 2 of Section III B. 2





t = argmaxm=1,...,N rt−l , u(m)  .

(54)

requires N2 multiplication for every iteration. As the algorithm is found to run for k iterations. Total number of multiplication required in Eq. (12) is kN2 . It can be noted in that step 4 of Section III B. 2, 

rˆt−1 =

(ut ) rt−1 

(ut ) ut

u(t )

(i)



these recursive equations run N times. In total, the computational complexity of MPS is, 2

(kN + 4kN + k2 N).

(60)

References

Summing up the number of computations, we get, kN + (N + c)



(55)

[1] J. Blacher, R. Asmar, S. Djane, G.M. London, M.E. Safar, Aortic pulse wave velocity as a marker of cardiovascular risk in hypertensive patients, Hypertension 33 (1999) 1111–1117. [2] J.J. Oliver, D.J. Webb, Noninvasive assessment of arterial stiffness and risk of atherosclerotic events, Arterioscler. Thromb. Vasc. Biol. 23 (2003) 554–566. [3] M. Cecelja, P. Chowienczyk, Role of arterial stiffness in cardiovascular disease, JRSM Cardiovasc. Dis. 1 (2012) 11. [4] Z.S. Hosseini, E. Zahedi, H.M. Attar, H. Fakhrzadeh, M.H. Parsafar, Discrimination between different degrees of coronary artery disease using time-domain features of the finger photoplethysmogram in response to reactive hyperemia, Biomed. Sig. Process. Control 18 (2015) 282–292. [5] D. Baldassarre, A. Hamsten, F. Veglia, U. de Faire, S.E. Humphries, A.J. Smit, P. Giral, S. Kurl, R. Rauramaa, E. Mannarino, et al., Measurements of carotid intima-media thickness and of interadventitia common carotid diameter improve prediction of cardiovascular events: results of the improve (carotid intima media thickness [imt] and imt-progression as predictors of vascular events in a high risk European population) study, J. Am. Coll. Cardiol. 60 (2012) 1489–1499. [6] J. Joseph, V. Jayashankar, A virtual instrument for automated measurement of arterial compliance, J. Med. Dev. 4 (2010) 045004. [7] A.K. Sahani, J. Joseph, M. Sivaprakasam, Automated system for imageless evaluation of arterial compliance, Engineering in Medicine and Biology Society (EMBC), 2012 Annual International Conference of the IEEE, IEEE, 227–231. [8] https://htic.iitm.ac.in/?q=artsens, 2016. [9] J. Stoitsis, S. Golemati, S. Kendros, K. Nikita, Automated detection of the carotid artery wall in b-mode ultrasound images using active contours initialized by the Hough transform, in: Engineering in Medicine and Biology Society. EMBS 2008. 30th Annual International Conference of the IEEE, IEEE, 2008, pp. 3146–3149. [10] D.E. Ilea, C. Duffy, L. Kavanagh, A. Stanton, P.F. Whelan, Fully automated segmentation and tracking of the intima media thickness in ultrasound video

T. Ganesh et al. / Biomedical Signal Processing and Control 38 (2017) 9–21

[11]

[12] [13]

[14]

[15]

[16] [17]

[18]

[19] [20] [21]

[22] [23] [24] [25]

[26]

sequences of the common carotid artery, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 60 (2013). Stephane Laurent, et al., Expert consensus document on arterial stiffness: methodological issues and clinical applications, Eur. Heart J. 27 (2006) 2588–2605. Stephane Laurent, et al., Abridged version of the expert consensus document on arterial stiffness, Artery Res. 1 (2007) 1–12. P.J. Brands, A.P. Hoeks, J. Willigers, C. Willekes, R.S. Reneman, An integrated system for the non-invasive assessment of vessel wall and hemodynamic properties of large arteries by means of ultrasound, Eur. J. Ultrasound 9 (1999) 257–266. G. Bambi, T. Morganti, S. Ricci, E. Boni, F. Guidi, C. Palombo, P. Tortoli, A novel ultrasound instrument for investigation of arterial mechanics, Ultrasonics 42 (2004) 731–737. M. Bastida-Jumilla, R.-M. Menchón-Lara, J. Morales-Sánchez, R. Verdú-Monedero, J. Larrey-Ruiz, J.-L. Sancho-Gómez, Frequency-domain active contours solution to evaluate intima-media thickness of the common carotid artery, Biomed. Sig. Process. Control 16 (2015) 68–79. D. Jegeleviˇcius, A. Lukoˇseviˇcius, Ultrasonic measurements of human carotid artery wall intima-media thickness, Ultrasound 43 (2014) 43–47. Y. Kawamura, Y. Yokota, F. Nogata, Estimation of carotid diameter with heartbeat on longitudinal b-mode ultrasonic images, 2008 30th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, IEEE, 3142–3145. M. Gutierrez, P. Pilon, S. Lage, L. Kopel, R. Carvalho, S. Furuie, Automatic measurement of carotid diameter and wall thickness in ultrasound images, in: Computers in Cardiology, IEEE, 2002, pp. 359–362. P.J. Touboul, M. Hennerici, et al., Mannheim carotid intima-medima thickness consensus, Cerebrovasc. Dis. 23 (2007) 75–80. R.S.J. Salonen, Ultrasound b-mode imaging in observational studies of atherosclerotic progression, Circulation 87 (Suppl. II) (1993) 56–65. J. Joseph, R. Radhakrishnan, S. Kusmakar, A.S. Thrivikraman, M. Sivaprakasam, Technical validation of artsens-an image free device for evaluation of vascular stiffness, IEEE J. Transl. Eng. Health Med. 3 (2015) 1–13. I. Selesnick, Introduction to Sparsity in Signal Processing, Connexions, 2012. F. Santosa, W.W. Symes, Linear inversion of band-limited reflection seismograms, SIAM J. Sci. Stat. Comput. 7 (1986) 1307–1330. E.J. Candes, J.K. Romberg, T. Tao, Stable signal recovery from incomplete and inaccurate measurements, Commun. Pure Appl. Math. 59 (2006) 1207–1223. T. Ganesh, J. Joseph, B. Bhikkaji, M. Sivaprakasam, Sparse models for determining arterial dynamics, Acoustics, Speech and Signal Processing (ICASSP), 2015 IEEE International Conference on, IEEE, 852–856. http://cvxr.com/cvx/, 2009.

21

[27] J.A. Tropp, A.C. Gilbert, Signal recovery from random measurements via orthogonal matching pursuit, IEEE Trans. Inf. Theory 53 (2007) 4655–4666. [28] S.G. Mallat, Z. Zhang, Matching pursuits with time-frequency dictionaries, IEEE Trans. Sig. Process. 41 (1993) 3397–3415. [29] S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004. [30] M. Aharon, M. Elad, A.M. Bruckstein, On the uniqueness of overcomplete dictionaries, and a practical way to retrieve them, Linear Algeb. Appl. 416 (2006) 48–67. [31] B. Efron, T. Hastie, I. Johnstone, R. Tibshirani, et al., Least angle regression, Ann. Stat. 32 (2004) 407–499. [32] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, Distributed optimization and statistical learning via the alternating direction method of multipliers, Found. Trends® Mach. Learn. 3 (2011) 1–122. [33] M.A. Figueiredo, R.D. Nowak, S.J. Wright, Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems, IEEE J. Sel. Top. Sig. Process. 1 (2007) 586–597. [34] T. Kailath, A.H. Sayed, B. Hassibi, Linear Estimation, vol. 1, Prentice Hall Upper Saddle River, NJ, 2000. [35] P. Myles, J. Cui, I. using the bland-altman method to measure agreement with repeated measures, Br. J. Anaesth. 99 (2007) 309–311. [36] Gernot Schulte-Altedorneburg, et al., Accuracy of in vivo carotid b-mode ultrasound compared with pathological analysis intima-media thickening, lumen diameter, and cross-sectional area, Stroke (2001) 1520–1526. [37] Jayaraj Joseph, et al., Artsens-an image-free system for noninvasive evaluating arterial compliance, 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) (2013) 4054–4057. [38] Jayaraj Joseph, et al., Artsens pen: a portable, image-free device for automated evaluation of vascular stiffness, IEEE International Symposium on Medical Measurements and Applications (MeMeA) (2016) 1–6. [39] Jayaraj Joseph, et al., Image-free evaluation of carotid artery stiffness using artsens: a repeatability study, 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (2014) 1957–1960. [40] S.L. Magda, et al., Comparative reproducibility of the noninvasive ultrasound methods for the assessment of vascular function, Heart Vessels 28 (2013) 143–150. [41] E.C. Godia, et al., Carotid artery distensibility: a reliability study, J. Ultrasound Med. 26 (2007) 1157–1165.