Measuring the effect of aging on vibrations of the carotid artery wall using empirical mode decomposition method.

Quantification of arterial elasticity and its dependency to age is considered in this paper. We use radiofrequency (RF) signals from carotid artery ultasonography to evaluate this dependency. Blood pressure, blood flow, and tethering to surrounding tissue are the main causes of the motion of the carotid wall. Tracking carotid artery wall motion from a series of ultrasound B-mode images is challenging due to the presence of noise and variable contrast. Moreover, the process of converting RF signals into the B-mode images causes some information to be lost. Hence, our goal is to extract the carotid wall motions and vibrations from RF signals. After extraction and removing the wall motion by using the phased tracking method combined with continuous wavelet transform, the vibrations of carotid inner wall in different subjects in different ages are compared with each other. Empirical mode decomposition method is used for extracting the first intrinsic mode function for different subjects' vibration and then their zero-crossing rates are compared. The results show the vibrations of the carotid inner wall are clearly decreased by age.

INTRODUCTION

Among patients with carotid atherosclerosis as the major cause of stroke and the third most common cause of death in most industrialized countries, only a few have warning symptoms. Athermanous plaques with relatively low stenosis degree may produce symptoms while high stenosis plaques may remain asymptomatic.[1] Age, systolic blood pressure and hyperlipidemia were found to be the important parameters associated with carotid artery elasticity values independently.[2]

The elasticity of the carotid wall is changed in the presence of atherosclerosis as a consequence of carotid wall composition and decreased by age. Blood pressure, blood flow, and tethering to surrounding tissue induce carotid wall motion. Compressive stresses in the radial direction and tensile stresses in the longitudinal and circumferential directions are caused by blood pressure.[3]

The carotid wall motion can be estimated using a number of techniques to quantify the carotid wall elasticity and contractility. A recognized indicator of the mechanical properties of carotid wall is radial movement, which is used in arterial disease examinations.[4]

Carotid wall motion can be quantitatively estimated by using image based methods that use temporal sequences of ultrasound B-mode images. Tracking carotid wall motion in real-time based on B-mode images is challenging since the motion is non-rigid and the rapidly acquired B-mode images contains noise. Moreover, the contrast and intensity of the B-mode images may change over time due to flowing blood.

Using different motion estimation methods leads to different results and interpretation. Although different suitable methods for carotid artery wall motion estimation are proposed, there is no quantitative evaluation of these methods. In the previous studies, the utilized methods were apparently selected randomly among a large number of image analysis methods.[3] Moreover, converting radiofrequency (RF) signals (echoes) into the B-mode images causes some information to be lost. Due to these reasons and more accurate view of the works[56] on RF signals, our goal is to extract the carotid wall motions and vibrations from RF signals. In this paper, a framework for quantitative evaluation of wall motion analysis and vibration extraction from RF ultrasound signals is introduced.

Related Works

Following, we briefly explain the existing related works on RF signals of arteries and heart wall.

Hasegawa et al. evaluated the elasticity of a local region of the arterial wall in.[7] Their proposed phased-tracking method detects the small deformation of the arterial wall by a precise estimation of its velocity.

Based on the change of the time delay for two-way propagation of ultrasound between an ultrasonic probe and the arterial wall, the phase shift of echoes is calculated to estimate the velocity. Although the phase shift can be estimated at each frequency within the frequency bandwidth of an echo, they showed that the velocity should be estimated at the centre frequency of an echo to prevent from being biased.[8]

This error during a frame interval is not so large, but it becomes significant as the result of temporal integration of the estimated velocity to track the position of the arterial wall.

In this paper, a new method is proposed to extract the vibrations of the inner carotid wall. The wall tracking method extracts the vibrations by using the phase change of Hilbert transforms of the RF echoes received from two corresponding layers in consecutive frames after eliminating the displacement of the wall.[9]

The carotid wall motion due to heart-beat is calculated and removed to extract wall vibrations. In order to prove the hypothesis that wall vibrations decrease while age increases, the vibrations of different subjects in different ages should be compared with each other. Empirical mode decomposition (EMD) method is used for extraction of vibrations.

EMD is a relatively new method to analyze the data by using an iterative process to obtain intrinsic mode functions (IMFs). Based on the analysis of the extracted IMFs, the carotid inner wall vibrations can be compared and evaluated in their related frequency bands. According to the EMD analysis, the first intrinsic mode function (IMF1), with the highest frequency and lowest amplitude, contains the carotid wall vibrations. The first IMFs obtained from different subjects' vibrations are then compared by their zero-crossing rates.

METHODS

In order to extract the carotid wall vibrations, carotid inner wall is more appropriate to analyze than the outer one. Because of the echogenic interface between the blood-intima-media layers, the carotid inner wall better reflects the ultrasound signal.[10]

Extraction of Vibrations on a Moving Object

The method for extraction of the low-amplitude vibrations from the motion of the carotid wall due to heart-beat is explained in this section.

Ultrasonic transducer transmits the RF pulses with a centre frequency of f0 and a time interval, T. The reflected echoes from the scatterers and/or different tissue interfaces are received by the same transducer [Figures 1 and 2]. The propagation delay, τ1(n), between the transducer and the reflector (scatterer or interface) at the nth transmission (frame) depends on the sound wave velocity, c0, and the distance, d1(n), between the transducer and the reflector:

Figure 1

Schematic of a new method for measurement of the vessel wall motion using an average phase of point scatterers in layers[11]

Figure 2

Estimation of the motion of the arterial wall by the phase tracking method in combination with Hilbert transforms. MT: Total number of assigned combinations; ML: Number of sampled points between two points of an assigned combination; T: Frame interval; Δx: Spacing of sampled points; d1 (n) and d1i(n): Depths of the lumen-intima and media-adventitia interfaces in the nth frame; xm (n) and x' m (n): mth combination of two points in the nth frame; νm (n) and ν' m (n): Velocities at xm (n) and x' m (n), respectively; Δh (n, x): Change in thickness at depth x in the nth frame [11]

Phase, θ1 (n) of a received echo is expressed by θ1 (n) =2πτ1 (n) f0 at the centre frequency, f0.[11] Therefore, distances d1 (n) and d1 (n + 1) between ultrasonic transducer and the reflector in the n-th and (n + 1)th frames are respectively defined by.

The average velocity of the reflector, ν1 (n) during a frame interval, T, is defined as:

The Hilbert transform is applied to the windowed and sampled RF echoes, in order to obtain phase, θ1(n), of an RF echo and the sampled complex signal, z (n; x), along an ultrasonic beam (x-axis).[9]

As shown in Figure 2, by referring to an M-mode image, the initial depth is manually assigned as the position of the reflector at the first frame (x1(1)”d1(1)). The phase,

x1(1) falls in the duration of the echo from d1(1).

Following describes the way we propose to obtain the phase shift, Δθ1(n), between two corresponding windows in consecutive frames.

The mother wavelet function of the continuous wavelet transform is first built based on four periods of a sine wave with frequency equal to the centre frequency of the transmitted ultrasound pulse multiplied by a Hanning window. The continuous wavelet transform is then applied on the two consecutive frames of RF data received by the ultrasonic transducer.

If we set the transmitted pulse as the mother wavelet, the time lag, that maximizes the cross correlation between absolute values of continuous wavelet transform of signals in the nth and (n + 1)th frames, can be computed. One of the signals is displaced by this amount of lag, and then the phase change, Δθ1(n), for two corresponding windows within a period equal to ΔT can be calculated by:

where θHIL (n + 1) and θHIL (n) are the phases of Hilbert transforms of two corresponding layers in successive frames which are unwrapped to reduce the large differences due to skipping phases from P to –π The average velocity, ν1 (n), in Eq. 4 is given by:

The next object position, , is estimated by multiplying the resultant velocity by the period ΔT:

Given , and the velocity , the object position as well as the motion velocity on a large amplitude motion are determined at the same time.

EMD

EMD was recently proposed to help analyzing the nonlinear and non-stationary time series. Previous signal processing methods used for medical ultrasound have been based on the assumption of a linear time-invariant system instead.[12]

Since real signals often consist of multiple frequency components and the Hilbert transform only gives a single value at each time instant, the instantaneous frequency is ambiguous. To prevent such and ambiguity and make the definition of instantaneous frequency more meaningful, the real part of the Fourier transform of the signal should only consist of positive frequencies. However, this restriction may make these solutions useless for non-stationary data. Huang[13] developed a modified version of this restriction that can be applied locally. They defined a class of functions, called IMFs, for which the instantaneous frequency can be defined anywhere. They also proposed EMD as a technique to break a signal down into a set of IMFs. Each IMF corresponds to a single oscillatory mode of the original data that achieves the requirement for the signal to be narrowband.

The intrinsic oscillatory modes can be identified based on their characteristic time scales in two below ways:

From the time between alternations of local maxima and minima

From the time between the zero crossings.

Since measuring the time between the zero crossings requires the mean of the data to be zero, EMD is based on the first way.

The EMD algorithm extracts the oscillatory mode that exhibits the highest local frequency of the data and leaves the remainder as a "residual." The complete decomposition of the data is obtained by applying the algorithm successively on the sequence of residuals. The final residual is a constant, a monotone trend, or a curve with a single extremum.[14] Given a signal, x (n), the IMFs are obtained using the algorithm called sifting. The sifting process can be summarized as follows:

The upper and the lower envelope are determined based on separately interpolating the local maxima and local minima points with a cubic spline, respectively. The mean envelope, as the half sum of the upper and the lower envelopes, is subtracted from the initial signal, and the same interpolation scheme is applied on the remainder. The sifting process terminates when the mean envelope is reasonably zero everywhere. The resultant signal is the first IMF and the same procedure after removing the previous IMFs iteratively yields the higher order IMFs.[15]

After the IMFs, c(n), 1 ≤ j ≤ M, have been calculated, decomposition of the original signal, x (n), into M empirical modes and a residue, r(n), is then achieved:

Zero-crossing Rate

Zero-crossing rate is one of the most straightforward methods used for computing the fundamental frequency. The zero-crossing rate is an indicator of the frequency at which the energy is concentrated in the signal spectrum. As the name suggests, the algorithm works by stepping through the waveform and measuring the rate of zero-crossing.[9]

At the last step of the proposed method, the zero-crossing rate of the first IMF1 is computed. The IMF1, as EMD method suggests, has the maximum frequency and minimum amplitude components. All steps are shown in the block diagram in Figure 3.

Figure 3

Block diagram of the vibration extraction and related evaluation

MATERIALS

Data

The probe and the scanner should be set properly during the study. Moreover, subjects and probe setup should ensure that the recorded arterial motion is only due to the hemodynamic forces. Deep breathing and increased pressure on the probe as sources of unwanted tissue displacement may affect tissue compliance. Hence, the sequences should be recorded with the minimal pressure on the probe. Furthermore, as the motion analysis algorithms estimate the relative motion between the artery and probe, it is important to minimize the probe motion to record the absolute arterial motion.

We recorded the ultrasonic data of the right common carotid artery of 12 healthy subjects (in two age groups (age diversity = 43.27 ± 17.39 years) with no history of cardiovascular disease, hypertension, or diabetes. All the subjects gave informed written consent for the examination. Their information is given in Table 1.

Table 1

Zero-crossing rates in subjects of all ages

Before ultrasonography, the subjects rested for at least 10 min in the supine position until their heart rate and blood pressure reached a steady state. In the present study, we used a MyLab™60 (ESAOTE, Italy) ultrasound system equipped with a LA532 linear array probe 4-13 MHz linear transducer with a frame rate of 60 frames per second that was connected to a PC for post processing. Each record contained approximately 6 cardiac cycles, and the right common carotid artery was scanned in the longitudinal view.

The RF signals are recorded and converted to data matrices using MATLAB (version 7.14.0.739 the MathWorks Inc., 2012). Matrices have dimensions of 2526 × 129 × 355 in average. The number of sample points in the depth direction is 2526 and the number of RF lines is 129, and there are 355 frames in 6 s of data acquisition for every case.

The next step is to extract the RF signals of the inner carotid wall that are used for phase tracking step. The inner and outer walls of the carotid artery are shown in Figure 4. A sample RF signal of the inner wall on the B-mode image is also shown in Figure 4 to better demonstrate the scanning protocol.

Figure 4

The inner and outer carotid wall (the B -mode image of a frame is shown to clarify). A sample RF signal of the inner wall on the B-scan is shown to better demonstrate the scanning protocol

RESULTS

Extracting the Vibrations of Carotid Wall

As shown in Figure 5, vibrations are carried by the inner wall motion that is larger in amplitude. The wall vibrations also have phases different from the ones of the wall motion. As it is mentioned in section "Methods", the vibrations are extracted after the wall motion is eliminated. In this paper, we consider the phase of wall vibrations to evaluate the effect of age on them. The inner wall motion and the corresponding extracted vibrations are shown in Figure 5a and b respectively.

Figure 5

(a) Vibrations and motion of a healthy data (b) Extracted vibrations

EMD Based Vibration Extraction

The extracted wall motion of a subject along with its IMFs and residual are shown in Figure 6. The first IMF appears to not follow the form of the original signal [Figure 6a], but the peaks corresponding to the wall motion are clear [Figure 6b]. The result from the second IMF shows a more likely trend to the original signal [Figure 6c]. In spite of losing some parts of the original form of the signal at this stage, six large wall motions resulting from heartbeat in the third IMF are still obvious [Figure 6d]. The other IMFs of the signal [Figure 6e–g] are not of interest to be considered because they do not provide us with more vibration related information.

Figure 6

Intrinsic mode functions (IMFs) obtained from empirical mode decomposition applying on the carotid wall. (a) Derived wall motion and vibrations on it, (b) IMF1, (c) IMF2, (d) IMF3, (e) IMF4, (f) IMF5, (g) Residual

Zero-Crossing Rate of the First IMF

The zero-crossing rates of the first IMFs for twelve subjects in different age and gender are presented in Table 1. The change of zero-crossing rate by age is illustrated in Figure 7.

Figure 7

Zero-crossing rates of the intrinsic mode function obtained via applying empirical mode decomposition on the interior wall motions of the carotid in subjects of all ages. The general trend graph coincides with the basic premise of vibration less at older ages. The horizontal axis is the subject's age. The linear fit shows the general trend

As it can be inferred from Figure 7 and Table 1, the subjects with label GO (77-year-old), S (59-year-old) and G (62-year-old) are old and the zero-crossing rates of their IMF1 are lower than the nine younger subjects.

DISCUSSION

In this paper, the phase tracking method in combination with Hilbert and wavelet transforms is first used to extract and eliminate the motion of the carotid inner wall. At the next step, the EMD method extracts the inner wall vibrations carried on the motion of the inner carotid wall caused by heartbeat. We then aim to measure the aging effect on the vibrations by comparing the zero-crossing rates of the obtained IMFs. It has been shown that, using this method, vibrations of carotid wall can be extracted with a respectable spatial resolution and comparing the zero-crossing rates of the IMFs proves the aging effect on carotid wall vibrations.

The major novelties of the proposed method are:

Using ultrasound RF signals that contain more carotid wall information in comparison to B-mode images

Extracting the carotid wall vibrations from carotid wall motion to consider them separately

Tracking the carotid wall using phase based method coupled with Hilbert transform and wavelet transform

Using EMD method to extract the wall vibrations.

In future investigations using the other properties of these vibrations, the aging effect will be determined more preciously. Furthermore, we are going to apply this method on cases with cardiovascular diseases for evaluating its ability to detecting the abnormalities.

BIOGRAPHIES

Fereshteh Yousefi Rizi received the B.S. degree in Biomedical Engineering from Shahed University, Tehran, Iran, in 2004 and the M.S. degree in Biomedical Engineering from Medical Physics and Biomedical Engineering Department, Faculty of Medicine, Tehran University of Medical Sciences, Tehran, Iran in 2008. Now she is a Ph.D. candidate in Biomedical Engineering, School of Electrical and Computer Engineering, University of Tehran, Tehran, Iran. Her research interests are Medical Signal and Image processing and Vascular Ultrasonography.

E-mail: f.yousefirizi@ut.ac.ir

Seyed Kamaledin Setarehdan received the B.Sc. degree in Electronic Engineering from the University of Tehran, Tehran, Iran, in 1989 and the M.Sc. degree in Biomedical Engineering from Sharif University of Technology, Tehran, Iran in 1992. He received the Ph.D. degree from the Electrical and Electronic Engineering Department of the University of Strathclyde in Glasgow, UK in 1998 in the field of medical image and signal processing. From March 1998 to December 2000 he was a postdoctoral research fellow in the Signal Processing Division of the Electrical and Electronic Engineering Department, University of Strathclyde in Glasgow, UK. In January 2000, he joined the Faculty of Electrical and Computer Engineering of the University of Tehran, Tehran, Iran. Dr. Setarehdan's main research interests are medical signal and image processing in general, medical ultrasound, EEG based brain computer interfacing and medical applications of the Near-Infrared Spectroscopy.

E-mail: ksetareh@ut.ac.ir

Hamid Behnam received the B.S. degree in Electrical Engineering from Iran University of Science and Technology, Tehran, Iran, in 1988, the M.S. degree in Medical Engineering from Sharif University of Technology, Tehran, Iran, in 1992, the Ph.D. degree in Applied Electronics from Tokyo Institute of Technology, Tokyo, Japan in 1998. Since 1998 till 2004, he was a researcher at Iran Research Organization for Science and Technology and from 2004 he has been a faculty member at Iran University of Science and Technology, in Tehran, Iran. Currently he is an Assistant Professor of Biomedical Engineering at the IUST. His research interests are Ultrasound in Medicine, Medical Image processing and Medical signal processing.

E-mail: behnam@iust.ac.ir