A novel power spectrum calculation method using phase-compensation and weighted averaging for the estimation of ultrasound attenuation

Ultrasonics 50 (2010) 592–599

Contents lists available at ScienceDirect

Ultrasonics journal homepage: www.elsevier.com/locate/ultras

A novel power spectrum calculation method using phase-compensation and weighted averaging for the estimation of ultrasound attenuation Seo Weon Heo a, Hyungsuk Kim b,* a b

School of Electronic and Electrical Engineering, Hongik University, Seoul 121-791, Republic of Korea Department of Electrical Engineering, Kwangwoon University, Seoul 139-701, Republic of Korea

a r t i c l e

i n f o

Article history: Received 14 September 2009 Received in revised form 15 December 2009 Accepted 16 December 2009 Available online 30 December 2009 Keywords: Ultrasound Attenuation Phase-compensation Block power spectrum Short-time Fourier analysis

a b s t r a c t An estimation of ultrasound attenuation in soft tissues is critical in the quantitative ultrasound analysis since it is not only related to the estimations of other ultrasound parameters, such as speed of sound, integrated scatterers, or scatterer size, but also provides pathological information of the scanned tissue. However, estimation performances of ultrasound attenuation are intimately tied to the accurate extraction of spectral information from the backscattered radiofrequency (RF) signals. In this paper, we propose two novel techniques for calculating a block power spectrum from the backscattered ultrasound signals. These are based on the phase-compensation of each RF segment using the normalized cross-correlation to minimize estimation errors due to phase variations, and the weighted averaging technique to maximize the signal-to-noise ratio (SNR). The simulation results with uniform numerical phantoms demonstrate that the proposed method estimates local attenuation coefficients within 1.57% of the actual values while the conventional methods estimate those within 2.96%. The proposed method is especially effective when we deal with the signal reflected from the deeper depth where the SNR level is lower or when the gated window contains a small number of signal samples. Experimental results, performed at 5 MHz, were obtained with a one-dimensional 128 elements array, using the tissue-mimicking phantoms also show that the proposed method provides better estimation results (within 3.04% of the actual value) with smaller estimation variances compared to the conventional methods (within 5.93%) for all cases considered. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction Medical ultrasound systems have been widely utilized for visualizing the internal organs of human body as well as for evaluating the types and pathological states of an object scanned. Various medical ultrasound parameters including the speed of sound [1–3], integrated backscatter [4–7], scatterer size [8–10] have been investigated in terms of the pathological information of patients. Ultrasound attenuation among medical ultrasound parameters provides not only important information for diseases in liver [11–13], breast [14], bone [15,16], intercostal tissue [17], and the myocardium [18], but also crucial backgrounds for estimating other parameters by compensating the amount of reflected signal decays along depths. Since attenuation in soft tissues generally demonstrate linear frequency-dependence [19,20], many attenuation estimation methods utilize the spectral information of reflected radiofrequency (RF) signals using the finite length of a gated window at each depth. Two major approaches for estimating the attenuation * Corresponding author. Tel.: +82 2 940 8373; fax: +82 2 940 5141. E-mail address: [email protected] (H. Kim). 0041-624X/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.ultras.2009.12.004

parameter in a spectral domain are spectral difference [21] and spectral shift [22]. While the spectral difference approach chiefly calculates the amplitude decay of the backscattered RF signal [19,23,24], the spectral shift approach estimates the downshift of power spectra along the propagation depth [25–28]. The spectral shift approach is generally more robust than the spectral difference approach at organ boundaries, but more sensitive to local spectral noise and has difficulty in compensating for diffraction effects. Most spectral domain algorithms, however, commonly utilize the spectral information of short-gated RF signals backscattered at different depths. Since the attenuation in soft tissues is frequency-dependent, the power spectrum at a deeper depth demonstrates downshift to a lower frequency range compared with that at a shallower depth. The short-time Fourier analysis method [29], widely used in the spectral analysis for estimating ultrasound parameters, calculates the block power spectrum at each depth by averaging several short-gated RF segments. Therefore, an improved estimation algorithm of the power spectrum from the observed RF signals is required. In this paper, we propose two novel techniques of calculating a block power spectrum in the short-time Fourier analysis. The methods are based on observing the signal characteristics of the

593

S.W. Heo, H. Kim / Ultrasonics 50 (2010) 592–599

2. Algorithms of extracting the spectral information from reflected RF signals 2.1. Phase and amplitude variations of reflected RF signals In the uniform medium, RF signals reflected from the same depth should be phase-aligned under the assumptions of Rayleigh scattering, the constant speed of sound, and the Born approximation. However, in analyzing the captured signal, we noticed that there existed small but evident phase and magnitude discrepancy among the adjacent data segments reflected from the same depth. One of the convincing explanations about the cause of those phase/ amplitude variations is that the tissue is often modeled as an aggregate of small sub-wavelength point scatters. Hence, the scattered signal from such an object is then the sum of many small reflections [30]. Then, by the central limit theorem, the reflected signal can be modeled as that a complex Gaussian random variable is multiplied to the original signal. The resulting signal pattern is represented by the speckle noise in the B-mode imaging. Wagner et al., similarly, showed the effect of coherent phase combining of multiple scattered signals, which described speckle noises in ultrasound images [31]. While there exist several works to reduce the speckle noise in the B-mode imaging using the linear/nonlinear signal processing techniques (such as [30,31]), as far as authors know, little work has been reported to consider those effects in the power spectrum estimation from the RF signal and that is one of the motivations of this paper. The effect of phase variation of the backscattered RF signal becomes more critical as the size of the gated window decreases since the number of acoustic waves in the gated window declines. When a taper function such as hamming (or hanning) window is applied to minimize spectral leakage at both ends of a reflected RF data segment, the effective number of acoustic waves becomes much smaller. To minimize these effects, the phase difference between the adjacent signal segments within a block should be compensated before the gating operation. The combination of randomly backscattered signals also induces the amplitude fluctuation in addition to the phase variation. Since the SNRs of RF data segments within the same block are different due to the amplitude variation, the weighted averaging could minimize the effect of spectral noise. This technique known as the maximal ratio combining (MRC) [32] has been widely adopted in the wireless communication system to improve the SNR of multi-path reflected signals in the Rayleigh fading channel where the signal is corrupted by multiplicative noises (very similar to the mechanism of speckle noise generation). Because backscattered RF signals from a deeper depth have a lower SNR than those from a shallower depth, the weighted combining scheme would

improve estimation performances, especially for the signals reflected back from the deeper depth. 2.2. Representation of reflected RF signals Under the assumptions of the identical acoustic path property for closely penetrated sound beams as well as fine lateral resolution compared with backscattering characteristics in the modern medical ultrasound transducer design, adjacent RF signals reflected from the same depth should be phase-aligned in uniform medium. In other words, the received RF signals closely located in parallel at the same depth could be assumed to bear no phase difference. In practice, however, backscattered RF data segments from the same depth exhibit random phase/amplitude variations as shown in Fig. 1, which shows three adjacent RF data segments of 1 mm long around the beam focus. These phase/amplitude variations do not come from the sources (or scatterers) because we assume the Rayleigh scattering in soft tissues. It is well-known that when the wavelength is large compared to a scattering object, an individual reflection from roughness features on the surface of the object fails to cause any noticeable interference effect [33]. Therefore, we model phase/ amplitude variations as the constructive and/or destructive combining of multiple unresolved reflected signals due to randomlydistributed small scatterers. For the linear array transducers, which are widely used in current ultrasound pulse-echo systems, the received signal (exactly, one A-line signal) consists of several responses from active piezoelectric elements for the beam focus and the apodization in the transmitting and receiving processes, respectively. These functions enable array transducers not only to enhance the imaging quality around the focal region but also effectively to reduce sidelobes. If we assume that L piezoelectric elements are excited when transmitting and receiving, the received signal (again, one A-line signal), r(t) can be written as [30]

rðtÞ ¼

L X i¼1

Ari

L X

Atj pðt  sri  stj Þ;

ð1Þ

j¼1

where p(t) is the transmitted pulse, and A and s are the weighting factors and time delays for beam focus (when transmitting) and apodization (when receiving), respectively. The subscripts t and r represent the transmitting and receiving processes, respectively.

1

Normalized Receiving Amplitude

received RF signal, which shows both the phase and signal level variations between RF segments. The proposed method includes phase-compensation using a cross-correlation between adjacent RF signals and the weighted averaging technique proportional to their signal-to-noise ratios (SNRs), and the result of which would extract spectral information with smaller mean estimation error and variance. The paper is organized as follows. The next section presents a theoretical derivation and details of the proposed techniques. What follows in the next section are an ultrasound simulation program and the overall procedures for estimating the attenuation coefficients. The next section accounts for simulation and experimental results using tissue-mimicking phantoms that illustrate the estimation performance with the proposed method. The final section discusses several aspects of the proposed techniques and summarizes the contributions of this paper.

0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8

(n-1)th segment n-th segment (n+1)-th segment

-1 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Relative Depth (mm) Fig. 1. Three adjacent RF segments reflected from the same depth (around the focus region). The amplitude and phase are changed due to randomly-distributed scatterers on their paths even in the uniform medium.

S.W. Heo, H. Kim / Ultrasonics 50 (2010) 592–599

Tissues are generally modeled as an aggregate of small subwavelength point scatterers. Therefore, the backscattered signal from a tissue is then the sum of many small reflections. The Eq. (1), however, does not reflect a combined effect of those point scatterers, so we modify the equation to consider them as follows. If we assume that there exist N scatterers at any given depth, the received signal from those scatteres, denoted by s(t), can be written as

sðtÞ ¼

N X

Bk rðt  sk Þ;

ð2Þ

k¼1

where Bk and sk are the reflection coefficient and time delay, respectively, and r(t) is given in Eq. (1). Since A and s in Eq. (1), which are design parameters of a transducer to achieve the maximum power transmission and reception, do not affect phase/amplitude variations in the received signal, we simply assume that s(t) is the transmitted signal to ease the explanation about the cases of phase/ amplitude variations. If we assume that the transmitted pulse is A cosðxc tÞ (the transmitted signal, in fact, is modulated around the center frequency with the Gaussian-like spectrum, but here we assume that only a pure sinusoidal signal is transmitted for simplicity), the backscattered signal, composed of several signals scattered from the small sub-wavelength point scatterers around the target region, can be described by

sðtÞ ¼

N X

Ak cosðxc ðt  sk ÞÞ

k¼1

¼

N X

Ak ðcos xc t  cos xc sk þ sin xc t  sin xc sk Þ

k¼1

¼

N X

! Ak cos xc sk cos xc t þ

n¼1

¼ R cosðxc ðt  /ÞÞ:

N X

! Ak sin xc sk sin xc t

n¼1

ð3Þ

If we assume that N is large, then the central limit theorem leads that R and / are modeled as random variables having the Rayleigh and uniform distributions, respectively. In fact, this phenomenon is observed in several other areas of research including the speckle noise pattern in ultrasound imaging, signal fading in wireless communication, and synthetic aperture radar (SAR) image processing [34]. As is clear from Eq. (3), the reflected signal suffers from the signal delay as well as the random amplitude variation. The amplitude variation due to the random variable R also affects the SNR of each of RF data segments. These phase/amplitude variations consequently affect the accuracy of block power spectra calculated. 2.3. Phase-compensation and weighted averaging of the block power spectrum In order to improve the accuracy of an averaged power spectrum, we propose two techniques accounting for the phase and SNR variations in computing the block power spectrum. Firstly, the phase mismatch (variation) of backscattered RF data segments within a block is compensated by comparing the normalized crosscorrelation between adjacent RF data segments. Fig. 1 shows phase differences between adjacent RF segments. In order to clearly show the phase difference rather than amplitude distortion due to a small SNR, three adjacent RF segments are selected from the middle of the lateral dimension at the beam focus region (4 cm) where the maximum power occurs. The length of each RF segment is 1 mm in the axial direction. As shown in Fig. 1, three adjacent segments look similar to one another because of the uniform medium, but present small variations in amplitude and phase since they

experience different scatterers and path characteristics. The normalized cross-correlation between adjacent RF segments is calculated across 10% of the window size at both ends, and the phasecompensated RF segment is extracted where the maximum cross-correlation occurs before computing the Fourier transform. The phase-compensation improves the quality of the estimated block power spectrum, especially when the size of the gated window becomes smaller because it contains relatively the small number of (effective) acoustic waves. Secondly, instead of the uniform averaging, the weighted averaging technique with respect to the total received RF signal power is applied to combine the individual Fourier spectrum of each RF data segment. Partly due to the random variation of amplitudes by the combining of several unresolved backscattered signals as well as due to the inhomogeneity of backscattered signal paths corresponding to each RF signal, the received signal levels (consequently, the SNRs) are different within the same block. Fig. 2 shows the average receiving power of RF data segments along depths, which are diffraction-compensated by the reference signal. Each block contains 19 RF segments (represented a 4 mm lateral block size), with a gated window length of 2 mm (in the axial direction). The error-bar represents the standard deviation in the lateral direction at each depth. As shown in Fig. 2, the average receiving signal powers gradually decrease as the acoustic beam propagates. However, the fluctuations of the average power (standard deviations) increase along depths since the signal powers become lower while the statistical noise levels are almost the same (because of uniform medium). Since we can reasonably assume that the power spectrum of RF segments with a low signal power (or low SNR) is less reliable, the estimation performance would be improved if they are properly weighted and averaged. The calculation of the optimal weighting factor in terms of the SNR is known as the maximal ratio combining (MRC), where the weighting factor is proportional to the SNR (or for simplicity, the received signal power) of an individual signal [32]. Therefore, the weighting averaging technique would maximize the overall SNR within each block. Now let us summarize our proposed methods with simple equations in the following. Assume that the transmitted signal is a single sinusoidal wave with a center frequency of xc and there are L reflected RF data segments in parallel (the lateral direction),

7

3

x 10

2.5

Average Receiving Power

594

2

1.5

1

0.5

0

0

1

2

3

4

5

6

7

8

Depth (cm) Fig. 2. Average receiving powers along depths. Each block contains 19 RF segments, which are gated with a window length of 4 mm. The error-bar represents the standard deviation in lateral direction at each depth.

595

S.W. Heo, H. Kim / Ultrasonics 50 (2010) 592–599

which are gated by a small window, within a block. The lth received signal (or lth A-line signal within a block) is given by

sðt; lÞ ¼ Rl cosðxc ðt  /l ÞÞ:

ð4Þ

The original (conventional) method calculates the Fourier spectrum of each data segment without consideration of phase and amplitude differences (Rl and /l) between data segments. Then non-coherent uniform-weighted combining of the Fourier spectra is applied only to improve the statistical accuracy. If we denote the Fourier spectrum estimator as P[x], the block power spectrum at depth z using the original (uniform averaging) method, Runi (f, z), can be described by

Runi ðf ; zÞ ¼

L 1X P½sðt; lÞ: L l¼1

ð5Þ

However, our approach compensates the phase (delay) and applies weighted averaging of each data segments, and the block power spectrum at the depth z using the proposed method, Rwgt (f, z), can be simply described by

Rwgt ðf ; zÞ ¼

L X

wl P½sðt  ^tl ; lÞ;

ð6Þ

l¼1

where ^tl is the estimated delay and wl is the weighting factor of each RF segment computing based on the observed signal level of lth segment. As shown in Eq. (6), the proposed method applies weighting factors proportional to the SNR after compensating the phase variation in the time domain while the original method simply averages without phase compensation, as shown in Eq. (5). The estimated delay, ^t l , is computed using a normalized cross-correlation between adjacent RF segments in the time domain and given by

^t l ¼ max cov½sðt þ Dt; lÞsðt; l þ 1Þ ; Dt

rl rlþ1

ð7Þ

where s(t, l) represents the lth RF segment in a block, cov[x] is a covariance operator, and rl and rl+1 are the standard deviations of s(t, l) and s(t, l + 1), respectively. As mentioned above, the maximum of Dt is set to 10% of the gated window length. 3. Procedures and methods for attenuation estimation 3.1. Overall procedures of spectral domain analysis Since most methods that estimate attenuation coefficients assume linear frequency-dependence, the spectral information of a reflected RF signal is widely utilized to estimate attenuation in soft tissues [23,25,35]. Both spectral difference and spectral shift approaches in estimating the ultrasound attenuation property, commonly, calculate the block power spectrum at each depth before applying their own estimation algorithms, and the short-time Fourier analysis is widely used to calculate block power spectrum using gated window on the reflected RF signal. In the short-time Fourier analysis, every RF signal in the axial direction is divided into relatively short segments using a gated window such as a hamming function to reduce spectral leakage at both boundaries. The Fourier spectra of several gated RF segments at the same depth (in the lateral direction) are averaged to improve the statistical accuracy, and finally various estimation algorithms in the spectral domain are applied using calculated block power spectra. As for the quality of a calculated block power spectrum, there is a tradeoff between accuracy and spatial resolution depending on the size of a gated window [35]. Briefly, the size of a gated window should be large enough to extract accurate (which means small bias and variance in block power spectrum estimation) spectral

information from the backscattered RF signal while it should also be small enough to satisfy the stationary assumption and provide sufficient spatial resolution of the final attenuation image estimated. The number of acoustic waves within a gated window determines the performance of a block power spectrum estimator. When the taper function such as a hamming window is applied, the effective number of acoustic waves for calculating the Fourier spectrum also decreases due to the slow-varying property at both ends. After selecting an appropriate size of gated window, each RF segment gated by a taper window is transformed by the fast Fourier transform (FFT) to represent spectral information. As mentioned above, several Fourier spectra at the same depth are simply averaged to obtain a block power spectrum to improve the statistical accuracy. Due to the inhomogeneous natures of soft tissues, however, adjacent RF segments in the lateral direction experience different acoustic path and scatterer characteristics. When more than 10 RF data segments (under the assumptions of 0.2 mm element-spacing in parallel and a minimum lateral block size of 2 mm) are averaged to compute the average block power spectrum, individual RF segments might have different phase delays and SNRs due to the tissue-inhomogeneity and random scattering. The phase variation and SNR differences would affect the accuracy of the averaged block power spectrum, especially when the size of the gated window decreases and the RF signals are received from a deeper depth. 3.2. RF data simulation and processing RF echo signals for both the reference and sample phantom with various attenuation coefficients and transducer properties were generated using an ultrasound simulation program [36]. This program was a frequency domain model based on the classical diffraction theory for continuous wave propagation [37]. In this program, the beam pattern of an ultrasound pulse transmitted by a series of sources (such as array elements) was computed using the superposition principle. For this study, a linear array transducer consisting of elements with a size of 0.2 mm by 10 mm and with a center-to-center distance of 0.2 mm was utilized. Each beam line was formed using 128 consecutive elements. The elevational focus was set to be the same as the lateral transmit focus. A dynamic receiving focus and a dynamic aperture were utilized such that the F-number was fixed at 2. The ultrasound field varied with the axial distance from the transducer. A Gaussian-shaped pulse with a center frequency of 5 MHz and a bandwidth of 80% was used as the incident pulse. The parameters of a simulated transducer are summarized in Table 1. Three uniform phantoms were simulated by randomly distributing 50 lm polystyrene beads and 25 lm glass beads in a medium that had a sound speed of 1540 m/s. The scatterer number density was set at 9.686 per cubic millimeter, which simulates the Rayleigh scattering. The modeled phantom dimensions were 40 mm (width) by 80 mm (height) by 10 mm (thickness), and the distance between adjacent beam lines was fixed at 0.2 mm. The

Table 1 Simulation parameters of a transducer. Transducer type

Linear array

Element size Number of elements Element spacing F-number Center frequency Bandwidth Beam focus

0.2 mm  10 mm 128 0.2 mm 2 5 MHz 80% 40 mm

596

S.W. Heo, H. Kim / Ultrasonics 50 (2010) 592–599

center of the phantom was placed at a 40 mm depth from the transducer, which coincided with the transmit and receive foci. The backscattered signals were computed by first applying a Gaussian filter centered at 5 MHz with the specified frequency bandwidth of the transducer, and then computing the inverse Fourier transform. This yielded the RF echo data for each beam line that was used to form the simulated ultrasound images. The reference phantom contained 50 lm polystyrene beads whose attenuation slope was 0.5 dB/cm/MHz, while two sample phantoms contained 25 lm glass beads which had the attenuation slope of 0.3 dB/cm/MHz and 0.7 dB/cm/MHz, respectively. After generating simulated RF data with a specified value of the attenuation coefficient, the entire RF data frame was divided into smaller overlapping 2-D blocks of sufficient sizes to obtain a consistent power spectrum. As mentioned above, the 2-D block size should be small enough to satisfy the stationary assumption, namely no variation in the attenuation, as well as to provide sufficient spatial resolution for the attenuation estimate. However, the block size has to be large enough to generate an accurate and robust power spectrum of backscattered RF signals. The criteria based on the evaluation of the Full-Width-Half-Maximum (FWHM) of the power spectrum were used to determine the block size that would provide a consistent power spectrum [35]. The block power spectrum was calculated using a short-time Fourier analysis by averaging the Fourier spectra from windowed segments within the block [29]. In this procedure, the phase-compensation and the weighted average techniques were applied. The gated window size was set to half the axial length of the block and applied a 50% overlapping in the axial direction. Each windowed RF segment was gated by a hanning window to minimize the spectral leakage artifact. A 50% overlapping of the 2-D blocks was also applied along the axial and lateral directions, respectively. Frequency smoothing where adjacent frequency estimates are averaged using a moving average window, was utilized to further reduce spectral noise artifacts in the power spectra [38]. 3.3. Experimental RF data acquisition Two tissue-mimicking (TM) phantoms were also used to evaluate the performance of the proposed method. While the reference phantom was uniform and consisted of 45–53 lm glass beads in a gelatin for which the attenuation slope was 0.5 dB/cm/MHz [39], the sample phantom consisted of 20 lm glass beads, and its attenuation slope was 0.8 dB/cm/MHz. The TM phantoms were scanned using a Siemens Antares ultrasound system (Siemens Medical Systems, Issaquah, WA, USA) using a linear array transducer operated at a 5 MHz center frequency with an 80% bandwidth. The transmit focus was set to the center of the scanning depth, 4 cm, and the region-of-interests (ROIs) were selected before and after the focus, 3 cm and 5 cm, respectively. RF data from these ROIs were acquired at five independent locations to improve the statistical accuracy. The average power spectra of ROIs were calculated using the proposed technique as well as the original method. 3.4. Attenuation estimation algorithms Although various algorithms of estimating attenuation coefficients have been proposed in the spectral domain, there are two major approaches. One is the spectral difference algorithm and the other is the spectral shift algorithm. To compare the performances of the proposed phase-compensated and weighted averaging method with the original uniform averaging method, we instead compare the performances of three different attenuation coefficient estimators, which all utilize the calculated block power

spectrum. They are centroid downshift [25], reference phantom [23] and hybrid method [40], respectively. In the centroid downshift method, the spectral shift fc(z) is estimated by calculating the centroid of the spectrum as a function of depth, and then the attenuation slope is estimated using

bs ðdB=cm=MHzÞ ¼ 

8:686 dfc ðzÞ  ; 4r2 dz

ð8Þ

where z is the depth of the region of interest from the transducer and r2 is the variance of the transmit pulse. In the reference phantom method, after compensating the diffraction and system-dependent effects using Eq. (9), the attenuation coefficient is calculated by taking the logarithm of this ratio of power spectra, RS(f, z), and a linear regression along the depth.

RSðf ; zÞ ¼

Ss ðf ; zÞ ¼ expf4ðbs  br Þfzg; Sr ðf ; zÞ

ð9Þ

where S(f, z) is the intensity of the backscattered RF signal received at the ultrasound transducer while the subscripts, r and s, represent the reference phantom and sample respectively. In the hybrid method, the Gaussian filtering at the transmit center frequency is applied to the ratio of power spectra between the reference and the sample, which are the same as Eq. (9). And, the attenuation slope is estimated using Eq. (8), similar to the centroid downshift method. Since the attenuation coefficient is linearly proportional to the derivative of center frequencies with respect to the distance (depth) as shown in Eq. (8), the absolute values of center frequencies are not required in the hybrid method.

4. Results and discussion 4.1. Simulation results with numerical phantoms In this section, we compare the performances of the three attenuation estimation methods described in the previous section using the original (uniform averaging) and the proposed (phase-compensated and weighted averaging) block power spectra. Since it is hard to directly compare the quality of the two block power spectra, the estimation bias and variance of attenuation estimation algorithms are presented in order to compare the estimation performances of two block power spectra. All the estimation parameters including block size, gated window, and linear fitting window are identical for each estimation method. Since the numerical phantoms are uniform, the estimated attenuation coefficients for each method are averaged in the entire region, and the estimation variances are computed for all the blocks to compare the estimation precision. Table 2 shows the estimation results for the three methods using a 4 mm gated window. In Table 2, the ‘RPM’ in the column of methods represents the reference phantom method and the standard deviations of the estimated attenuation coefficients are noted with parenthesis for each case. Although all the attenuation coefficients estimated are around the real value of the sample phantoms, 0.3 dB/cm/MHz and 0.7 dB/cm/MHz, respectively, the proposed method estimates the local attenuation coefficients within 1.57% of the actual values, while the conventional methods estimate those within 2.96%. Since simulated numerical phantoms contain (statistically) well-defined noises and have no physical errors such as analog-to-digital conversion and quantization processes, the estimated average values of the estimator are quite accurate and the estimation variances are also stable for all cases when they are calculated in the entire (relatively large) region. However, the proposed method achieves relatively better estimation performance for all cases, and it could mean the proposed method is a better power spectrum estimator.

597

S.W. Heo, H. Kim / Ultrasonics 50 (2010) 592–599 Table 2 Estimation results of simulated numerical phantoms (unit: dB/cm/MHz). Reference: 0.5/sample: 0.3

Reference: 0.5/sample: 0.7

Original

Proposed

Original

Proposed

RPM Centroid Hybrid

0.3076 (0.0688) 0.2863 (0.1149) 0.3091 (0.0743)

0.3040 (0.0547) 0.3006 (0.0864) 0.3084 (0.0536)

0.6952 (0.0599) 0.6577 (0.1247) 0.7063 (0.0683)

0.6997 (0.0425) 0.6695 (0.0939) 0.7051 (0.0449)

4.2. Experimental results using TM phantoms

Estimated Attenuation Coefficient (dB/cm/MHz)

The estimation performances of the proposed techniques are also verified with uniform TM phantoms mentioned in the

(a)

original proposed

0.712

0.71

0.708

0.706

0.704

0.702 0.7

2mm

(b)

4mm

5mm

6mm

0.8 0.78 0.76 0.74 0.72 0.7 0.68 0.66 0.64 original proposed

0.62 2mm

0.4

3mm

Window Size

Estimated Attenuation Coefficient (dB/cm/MHz)

Note that the centroid downshift method exhibits relatively larger estimation errors among estimation methods because its performance is closely related to the compensation of diffraction effects due to the beam focusing. Although the estimation variance of the centroid downshift method is still larger than other estimation methods, the proposed method makes the estimation bias and variance smaller in this case, too. Fig. 3 shows the estimation bias and variance of the reference phantom method along the axial direction. The attenuation coefficient of the sample was 0.3 dB/cm/MHz and a 4 mm gated window was utilized. The lateral averaging was applied to calculate the attenuation coefficients at each depth in order to present the estimation performance as the acoustic beam propagates. Since the SNR of backscattered RF signals generally decreases due to the accumulated attenuation and longer travelling path, the estimated attenuation coefficients in deeper depths deviate from the real value. However, the proposed technique gives smaller estimation bias and variance in all depths. The estimation performances with various sizes of gated windows, from 2 mm to 6 mm, using the hybrid estimation method are shown in Fig. 4. As shown in Fig. 4a, the estimation accuracy is improved as the gated window becomes longer in both block power spectra while the block power spectrum calculated by the proposed method demonstrates better estimation performance for all gated window cases. For the standard deviations of estimated attenuation coefficients shown in Fig. 4b, the proposed block power spectrum also achieves 32.26% smaller standard deviation (on average) compared to the original block power spectrum for all window cases.

Estimated Attenuation Coefficient (dB/cm/MHz)

Method

3mm

4mm

5mm

6mm

Window Size original proposed

0.38

Fig. 4. Estimation performance with various sizes of gated windows using the hybrid method. (a) Estimated attenuation coefficients. (b) Estimated attenuation coefficients with standard deviations.

0.36 0.34 0.32 0.3 0.28 0.26 0.24 0.22

0

1

2

3

4

5

6

7

8

9

Depth (cm) Fig. 3. Estimated attenuation coefficients along depths using the reference phantom method. The attenuation coefficient of the sample is 0.3 dB/cm/MHz, and a 4 mm gated window is utilized. The estimation bias and variance of the proposed method are improved for the entire depth.

previous section. The attenuation coefficients of reference and sample phantoms were 0.5 dB/cm/MHz and 0.8 dB/cm/MHz, respectively, and the transmit center frequency was 5 MHz. Two ROIs were selected before and after the focus (3 cm and 5 cm), while the transmit focus was set to 4 cm. A 5 mm gated window was utilized and 50% overlaps in the axial and the lateral directions were allowed for all block power spectra to improve the statistical accuracy. Table 3 shows the estimation results using the original (uniform averaging) and the proposed (phase-compensated and weighted averaging) block power spectra for three estimation methods. As shown in Table 2, the ‘RPM’ represents the reference phantom method and the standard deviations of estimation values are noted in parenthesis. Similar to simulation results, the proposed block power spectra achieve better results that estimate the attenuation coefficients within 3.04% of the actual value with smaller estima-

598

S.W. Heo, H. Kim / Ultrasonics 50 (2010) 592–599

Table 3 Estimation results of tissue-mimicking phantoms (unit: dB/cm/MHz). Method

Original

Proposed

3 cm

5 cm

3 cm

5 cm

RPM Centroid Hybrid

0.8351 (0.1148) 0.9014 (0.1542) 0.8507 (0.1224)

0.7856 (0.1339) 0.8679 (0.1682) 0.8153 (0.1351)

0.8127 (0.0894) 0.8863 (0.1337) 0.8110 (0.1195)

0.8011 (0.0963) 0.8273 (0.1394) 0.8078 (0.1308)

tion variances, while the conventional methods estimate those within 5.93%. The overall estimation results at a deeper depth (5 cm) have relatively larger standard deviations than those at a shallower depth (3 cm) because the SNR quickly decreases after the beam focus. Note that the estimated attenuation coefficients at 3 cm are overestimated than those at 5 cm since the diffraction effects due to beam focusing is not correctly compensated. The imperfect compensation of the diffraction effect could provide a bias to the overall estimation procedures. Similar to the simulation results, the distinct improvement of estimation performances for the reference phantom method is achieved among other estimation methods since an accurate calculation of the block power spectrum is the most important factor when the frequency compounding is not applied.

[3]

[4]

[5]

[6]

[7]

[8]

5. Conclusions The Quantitative Ultrasound (QUS) analysis has been intensively studied to provide quantitative information from the reflected RF signals. Estimation of the ultrasound attenuation property is one of the crucial parameters since it not only provides pathological information on the tissues scanned but also enables to estimate other ultrasound parameters accurately by compensating the accumulated signal decay. Most of the QUS analyses, however, utilize spectral information obtained by short-gated RF signals such as the short-time Fourier analysis. Therefore, a correct estimation of spectral property in all depths from the backscattered RF signals provides a better fundamental background for further analysis. In this paper, we propose two techniques for using the phasecompensation and the weighted averaging within a block to calculate a more accurate block power spectrum in the short-time Fourier analysis method. These techniques reduce estimation errors of the block power spectrum caused by random phase variations of the RF data segments reflected from the same depth, and improve the estimation performance of the block power spectrum by properly combining the spectral information of individual RF data segments. Ultrasound simulation results with uniform numerical phantoms demonstrate that the estimation bias using the proposed techniques (1.57% of the actual value) is smaller than that using the conventional methods (2.96%) for all three different estimation algorithms considered. Experimental results with TM phantoms also show that the proposed techniques provide smaller estimation bias (3.04% of the actual value) than the conventional methods (5.93%). In addition, as the SNR of the reflected signal becomes lower (generally at deeper depths) and the gated window gets smaller, the proposed techniques achieve better estimation performances in terms of bias and variance for all attenuation estimation methods. References [1] D.E. Robinson, J.O. Ophir, L.S. Wilson, C.F. Chen, Pulse-echo ultrasound speed measurements: progress and prospects, Ultrasound Med. Biol. 17 (1991) 633– 646. [2] U. Techavipoo, T. Varghese, Q. Chen, T.A. Stiles, J.A. Zagzebski, G.R. Frank, Temperature dependence of ultrasonic propagation speed and attenuation in

[9]

[10]

[11] [12]

[13]

[14]

[15] [16]

[17]

[18]

[19] [20]

[21] [22] [23]

[24]

[25]

[26]

[27]

excised canine liver tissue measured using transmitted and reflected pulses, J. Acoust. Soc. Am. 115 (6) (2004) 2859–2865. Y. Levy, Y. Agnon, H. Azhari, Measurement of speed of sound dispersion in soft tissues using a double frequency continuous wave method, Ultrasound Med. Biol. 32 (7) (2006) 1065–1071. Z. Vered, G.A. Mohr, B. Barzilai, C.J. Gessler, S.A. Wickline, K.A. Wear, T.A. Shoup, A.N. Weiss, B.E. Sobel, J.G. Miller, Ultrasound integrated backscatter tissue characterization of remote myocardial infarction in human subjects, J. Am. Coll. Cardiol. 13 (1989) 84–91. K.A. Wear, T.A. Stiles, G.R. Frank, E.L. Madsen, F. Cheng, E.J. Feleppa, C.S. Hall, B.S. Kim, P. Lee, W.D. O’Brien Jr., M.L. Oelze, B.I. Raju, K.K. Shung, T.A. Wilson, J.R. Yuan, Interlaboratory comparison of ultrasonic backscatter coefficient measurements from 2 to 9 MHz, J. Ultrasound Med. 24 (9) (2005) 1235–1250. S.L. Bridal, C. Fournier, A. Coron, I. Leguerney, P. Laugier, Ultrasonic backscatter and attenuation (11–27 MHz) variation with collagen fiber distribution in ex vivo human dermis, Ultrason. Imaging 28 (1) (2006) 23–40. L.R. Taggart, R.E. Baddour, A. Giles, G.J. Czarnota, M.C. Kolios, Ultrasonic characterization of whole cells and isolated nuclei, Ultrasound Med. Biol. 33 (3) (2007) 389–401. M.F. Insana, R.F. Wagner, Describing small-scale structure in random media using pulse-echo ultrasound, J. Acoust. Soc. Am. 87 (1990) 179–192. T.J. Hall, M.F. Insana, L.A. Harrison, G.G. Cox, Ultrasonic measurement of glomerular diameters in normal adult humans, Ultrasound Med. Biol. 22 (8) (1996) 987–997. W. Liu, J.A. Zagzebski, T. Varghese, A.L. Gerig, T.J. Hall, Spectral and scatterersize correlation during angular compounding: simulations and experimental studies, Ultrason. Imaging 28 (4) (2006) 230–244. P.A. Narayana, J. Ophir, On the frequency dependence of attenuation in normal and fatty liver, IEEE Trans. Sonics Ultrason. 30 (1983) 379–383. Y. Fujii, N. Taniguchi, K. Itoh, K. Shigeta, Y. Wang, J.W. Tsao, K. Kumasaki, T. Itoh, A new method for attenuation coefficient measurement in the liver: comparison with the spectral shift central frequency method, J. Ultrasound Med. 21 (7) (2002) 783–788. M. Meziri, W.C. Pereira, A. Abdelwahab, C. Degott, P. Laugier, In vitro chronic hepatic disease characterization with a multiparametric ultrasonic approach, Ultrasonics 43 (5) (2005) 305–313. S.W. Huang, P.C. Li, Ultrasonic computed tomography reconstruction of the attenuation coefficient using a linear array, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 52 (11) (2005) 2011–2022. K.A. Wear, Ultrasonic attenuation in human calcaneus from 0.2 to 1.7 MHz, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 48 (2001) 602–608. K.A. Wear, Characterization of trabecular bone using the backscattered spectral centroid shift, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 50 (4) (2003) 402–407. R.T. Towa, R.J. Miller, L.A. Frizzell, J.F. Zachary, W.D. O’Brien Jr., Attenuation coefficient and propagation speed estimates of rat and pig intercostal tissue as a function of temperature, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 49 (10) (2002) 1411–1420. S.L. Baldwin, K.R. Marutyan, M. Yang, K.D. Wallace, M.R. Holland, J.G. Miller, Measurements of the anisotropy of ultrasonic attenuation in freshly excised myocardium, J. Acoust. Soc. Am. 119 (2006) 3130–3139. R. Kuc, Bounds on estimating the acoustic attenuation of small tissue regions from reflected ultrasound, Proc. IEEE 73 (1985) 1159–1168. S.W. Flax, N.J. Pelc, G.H. Glover, F.D. Gutmann, M. McLachlan, Spectral characterization and attenuation measurements in ultrasound, Ultrason. Imaging 5 (1983) 95–116. B. Zhao, O.A. Basir, G.S. Mittal, Estimation of ultrasound attenuation and dispersion using short time Fourier transform, Ultrasonics 43 (2004) 375–381. R. Kuc, H. Li, Reduced-order autoregressive modeling for center-frequency estimation, Ultrason. Imaging 7 (1985) 244–251. L.X. Yao, J.A. Zagzebski, E.L. Madsen, Backscatter coefficient measurements using a reference phantom to extract depth-dependent instrumentation factors, Ultrason. Imaging 12 (1990) 58–70. L.S. Wilson, D.E. Robinson, B.D. Doust, Frequency domain processing for ultrasonic attenuation measurement in liver, Ultrason. Imaging 6 (1984) 278– 292. M. Fink, F. Hottier, J.F. Cardoso, Ultrasonic signal processing for in vivo attenuation measurement: short time Fourier analysis, Ultrason. Imaging 5 (1983) 117–135. C. Kasai, K. Namekawa, A. Koyano, R. Omoto, Real-time two-dimensional blood flow imaging using an autocorrelation technique, IEEE Trans. Sonics Ultrason. 32 (1985) 458–464. T. Baldeweck, P. Laugier, A. Herment, G. Berger, Application of autoregressive spectral analysis for ultrasound attenuation estimation: interest in highly

S.W. Heo, H. Kim / Ultrasonics 50 (2010) 592–599

[28]

[29]

[30]

[31] [32] [33]

attenuating medium, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 42 (1995) 99–109. J. Ophir, R.E. McWhirt, N.F. Maklad, P.M. Jaeger, A narrowband pulse-echo technique for in vivo ultrasonic attenuation estimation, IEEE Trans. Biomed. Eng. 32 (1985) 205–212. P. Welch, The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms, IEEE Trans. Audio Electroacoust. 15 (1967) 70–73. Murtaza Ali, Dave Magee, Udayan Dasgupta, Signal Processing Overview of Ultrasound Systems for Medical Imaging, Texas Instruments White Paper, 2008, pp. 1–26. R.F. Wagner, S.W. Smith, J.M. Sandrik, H. Lopez, Statistics of speckle ultrasound B-scans, IEEE Trans. Sonics Ultrason. 30 (1983) 156–163. J.H. Winters, Optimum combining in digital mobile radio with cochannel interference, IEEE Trans. Veh. Technol. 33 (3) (1984) 144–155. T.H. Szabo, Diagnostic Ultrasound Imaging: Inside Out, Elsevier Academic Press, 2004. pp. 213–221.

599

[34] L. Denis, F. Tupin, J. Darbon, M. Sigelle, SAR image regularization with fast approximate discrete minimization, IEEE Trans. Image Process. 18 (7) (2009) 588–1599. [35] H. Kim, T. Varghese, Attenuation estimation using spectral cross-correlation, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 54 (2007) 510–519. [36] Q. Chen, Computer Simulations in Parametric Ultrasonic Imaging (Ph.D. Dissertation), University of Wisconsin-Madison, 2004. [37] Y. Li, J.A. Zagzebski, A frequency domain model for generating B-mode images with array transducers, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 46 (1999) 690–699. [38] T. Varghese, K.D. Donohue, Estimating mean scatterer spacing with the frequency-smoothed spectral autocorrelation function, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 42 (1995) 451–463. [39] Thaddeus Wilson, James Zagzebski, Yadong Li, A test phantom for estimating changes in the effective frequency of an ultrasonic scanner, J. Ultrasound Med. 21 (2002) 937–945. [40] H. Kim, T. Varghese, Hybrid spectral domain method for attenuation slope estimation, Ultrasound Med. Biol. 34 (11) (2008) 1808–1819.