Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3879804

Formats

Article sections

Authors

Related links

Ultrason Imaging. Author manuscript; available in PMC 2014 January 3.

Published in final edited form as:

PMCID: PMC3879804

NIHMSID: NIHMS532651

Ivan M. Rosado-Mendez, Department of Medical Physics, University of Wisconsin, 1111 Highland Ave., Madison, WI 53705, USA. Email: ude.csiw@zednemodasor

The publisher's final edited version of this article is available at Ultrason Imaging

See other articles in PMC that cite the published article.

Reported here is a phantom-based comparison of methods for determining the power spectral density of ultrasound backscattered signals. Those power spectral density values are then used to estimate parameters describing *α(f)*, the frequency dependence of the acoustic attenuation coefficient. Phantoms were scanned with a clinical system equipped with a research interface to obtain radiofrequency echo data. Attenuation, modeled as a power law *α(f)=α _{0}f^{β}*, was estimated using a reference phantom method. The power spectral density as estimated using the short-time Fourier transform (STFT), Welch's periodogram, and Thomson's multitaper technique, and performance was analyzed when limiting the size of the parameter estimation region. Errors were quantified by the bias and standard deviation of the

Quantitative Ultrasound (QUS) methods attempt to characterize tissue microstructure to aid in disease diagnosis and treatment monitoring.^{1,2} Two QUS parameters with demonstrated potential are the frequency dependent attenuation coefficient, *α(f)* and the backscatter coefficient, *σ _{b}(f)*.

Most methods for estimating *α(f)* and *σ _{b}(f)* rely on a spectral analysis of radiofrequency (RF) echo signals by computing their power spectral density.

Each periodogram is the squared magnitude of the Fourier transform of the gated scanline signal.^{22,23} Other methods for estimating the power spectral density based on the periodogram include Welch's method and Thomson's multitaper method. In Welch's method, the periodogram resulting from one RF echo-signal segment is the average of the periodograms of shorter and overlapping subsections of the original segment.^{28,29} This method was recently suggested as optimum for *σ _{b}(f)* estimations.

In practice, the frequency response of the gating window is convolved with the true RF signal spectrum. The frequency response of the windowing operation is characterized by a central main lobe bracketed by lower-amplitude side lobes. The main-lobe width, which is inversely proportional to the length of the estimation window, defines the frequency resolution, i.e., the amount of correlation between values of the power spectral density at neighboring frequency points. The side-lobe amplitude defines the bias introduced by spectral leakage. Bias can be reduced by multiplying each RF signal segment within the estimation window by a tapering function such as the Hamming or Hann, but this leads to wider main lobes, and affects the accuracy and precision of power spectral density estimates.

The area defined by the estimation window is illustrated in figure 1. The region of interest (ROI) outlined on the image might include a lesion or a section of it and the surrounding tissue. Within the ROI, a “parameter estimation region” defines the total amount of echo data applied to obtain an estimate of the QUS parameter of interest. This parameter estimation region can include one or more power spectral density estimation windows, depending on the estimated QUS parameter. The number of parameter estimation regions within an ROI depends on the task to be performed. In **bulk parameter estimation** the ROI includes usually one parameter estimation region. This can be thought of as a “large parameter estimation region” mode aimed at comprehensively gathering values of tissue properties to be used in tissue classification.^{22} In contrast, in **parametric imaging** small parameter estimation regions are positioned at different locations within the ROI to display spatial variations of the QUS parameter.^{22}

Description of a parameter estimation region (PER) from which one estimate of the attenuation coefficient is obtained. Each power spectral density (PSD) estimate is assigned to the center of the spectral estimation window, indicated by a dot. The PER **...**

Algorithms for local estimation of *α*(*f*) characterize the variation of a spectral property along the propagation depth. At least three methods are quantitatively meaningful because they compensate for diffraction and system settings. These are the spectral difference method (included as part of the reference phantom method),^{13} the spectral log difference method, and the hybrid method.^{19} In general, the reference phantom method requires smaller parameter estimation regions, although still fairly long (some cases suggesting up to 50mm), compared to the other methods.^{31-36} Although this does not limit bulk-parameter estimation, it is a major limitation in terms of the spatial resolution of parametric images. In this sense, it has been suggested that Thomson's multitaper method may outperform the conventional short-time Fourier transform and Welch's power spectral density estimate in terms of the overall diagnostic performance of parametric images in which the size of the parameter estimation region is severely reduced.^{26}

This work presents a phantom-based analysis of constraints in the size of the parameter estimation region when estimating the attenuation coefficient for two different tasks: either bulk-parameter estimation or parametric image formation. The analysis is based on a power-law model for the frequency dependence of attenuation.^{37-44} We present a comparison of the performance of two spectral estimation methods commonly applied in QUS when estimating the attenuation coefficient (the short-time Fourier transform and Welch's periodogram) with Thomson's multitaper method. Results quantify the accuracy and precision of the estimates as the amount of data incorporated into the estimation is reduced.

The reference phantom method^{13} was used to estimate the attenuation coefficient of a sample from RF echo signals obtained from a clinical scanner. The reference phantom method accounts for system-dependent parameters (transducer and pulser-receiver transfer functions) using power spectral density estimates from a reference material (with speed of sound equal to that of the sample), generated at the same depths as estimates from the sample. The sample consisted of a homogeneous tissue-mimicking phantom composed of water-based gel containing graphite powder (50mg/cm^{3} of agar) to control the attenuation and 3000E glass beads (Potters Industries, Inc., Valley Forge, PA; 5 to 40μm diameter, 4mg/cm^{3} of agar) to provide scattering.^{45} The reference phantom was a similar homogeneous material, consisting of an emulsion of 70% safflower oil in gelatin and also containing 3000E glass-bead scatterers (4mg/cm^{3} of emulsion). Both the sample and the reference phantoms were in acrylic boxes that had 25μm thick Saran wrap® (Dow Chemical, Midland, MI) scanning windows.

During phantom fabrication, 2.5cm-thick test samples were prepared to measure the sound speed and attenuation of the phantom materials. Laboratory estimates of *α(f)* and the speed of sound of these materials were performed using a narrowband substitution technique^{46} at frequencies from 2.25 to 10MHz. The attenuation coefficient was modeled as a power of frequency *f*, using:

$$\alpha \left(f\right)={\alpha}_{0}{f}^{\beta}.$$

(1)

Table I shows the resultant *α _{0}* and

The phantoms were scanned with an 18L6 linear array transducer on a Siemens Acuson S2000 clinical imaging system (Siemens Medical Solutions USA, Inc, Malvern, PA). The transducer was driven by single cycle 9MHz pulses. The system was equipped with the Axius Direct Ultrasound Research Interface^{48} that allows acquisition of beamformed RF echo signals after time-gain compensation. The transmit focus was placed as shallow as possible (2.5mm) to estimate attenuation in regions distal to the focus, in accordance with previous results.^{47,49} RF echoes for 20 parallel image planes from the sample and 25 parallel planes from the reference phantom were acquired by sweeping the transducer elevationally between successive planes in steps equal to half the height of the transducer's elements. Each plane scanned provided a 2D RF echo-signal frame consisting of 456 acoustic scanlines spaced at 0.128mm. The RF data were digitized at 40MHz for a total of 2076 samples spanning 39mm axially.

The size of the scattering volume determines the correlation among neighboring RF echo samples. This volume is affected by the pulse bandwidth and the diffraction field. Therefore, the number of uncorrelated axial samples and acoustic scanlines within a spectral estimation window or a parameter estimation region depends on the experimental conditions. To generalize our results, we parameterized the size of the spectral estimation window and the area spanned by the parameter estimation region in terms of the number of uncorrelated axial (time) and lateral (acoustic scanlines) samples of RF echo signals. For this purpose we computed the separation between the two closest axial samples and closest lateral samples considered to be uncorrelated.^{20,26,50-53} For segments including *m* axial samples of RF echo signals from *N* adjacent acoustic scanlines, the lateral correlation coefficient is defined as:^{26}

$${\rho}_{L}\left(\Delta \right)=\frac{1}{N}\sum _{j=1}^{N}\frac{\sum _{i=1}^{m}\left[({X}_{j,i}-\overline{{X}_{j}})({X}_{j+\Delta ,i}-\overline{{X}_{j+\Delta}})\right]}{\sqrt{\sum _{i=1}^{m}\left[{({X}_{j,i}-\overline{{X}_{j}})}^{2}{({X}_{j+\Delta ,i}-\overline{{X}_{j+\Delta}})}^{2}\right]}}$$

(2)

where *X _{j,i}* is the

The parameter estimation region required to obtain one *α(f)* estimate includes a set of axially-aligned spectral estimation windows, as illustrated in figure 1. The size of the parameter estimation region is defined axially by the number of spectral estimation windows, and laterally by their widths.

It is important to notice the trade-offs in the axial size of a spectral estimation window and the parameter estimation region. Each spectral estimation window should be long enough so that the power spectral density estimate is not dominated by windowing and tapering-function effects while also short enough that nonstationarity effects (diffraction, spatial variations in scattering or attenuation) can be considered negligible. Following this criterion, the maximum spectral estimation window length that was explored in this work was smaller than 10mm. The axial span of the parameter estimation region should be long enough to be sensitive to attenuation-dependent variations in the power spectral density amplitude over depth.

The size of the spectral estimation window was varied both axially and laterally, and the effects of these variations on the accuracy and precision of estimates of the power-law parameters of *α(f)* were studied. The window length was reduced from 20 pulse-echo correlation lengths down to 2.5 pulse-echo correlation lengths leading to window lengths from 7.5mm down to 0.9mm. The number of uncorrelated acoustic lines included in the spectral estimation window, and therefore in the parameter estimation region, was varied from 22 down to 4, resulting in widths from 10.9mm down to 2.0mm. Fourier transforms were computed after zero-padding the echo signals segments to 1024 points.

For each spectral estimation window size, a periodogram was computed by each of three methods: The Short-time Fourier transform (STFT), Welch's method, and Thomson's multitaper method. The STFT estimate of the periodogram is computed as the squared magnitude of the Fourier transform of the RF echo segment multiplied by a tapering function. Two tapering functions were used: the rectangular and Hann functions.^{28,29} In Welch's method, the RF echo segment is first divided into N_{SS} overlapped subsegments. The power spectral density estimate from Welch's method is the average of the periodograms of each of the subsegments after tapering.^{54} The subsegment length and the overlap ratio were 50% of the original segment, and 50% of the subsegment length, values commonly used in QUS.^{26,29,55-59} The Hann tapering function was used for subsegment tapering. The multitaper method requires defining a value for the time-half bandwidth product *NW* (described below). Then, a set of tapered signal segments is obtained by multiplying an RF echo-signal segment by each one of the first N_{DPSS}=2*NW*-1 Discrete Prolate Spheroidal Sequences (DPSS). These tapering functions contain the greatest amount of energy within a frequency band 2*W*.^{29,60,61} The squared magnitude of the Fourier transform of each tapered signal segment is referred to as an eigenspectrum, and the final power spectral density estimate is given by the weighted average of the set of N_{DPSS} eigenspectra.^{29,60} The first three DPSS for 256 time samples and *NW*=3 are displayed in figure 2.

First three Discrete Prolate Spheroidal Sequences for 256 time samples and a time –half bandwidth product *NW*=3.^{61}

The rectangular and Hann are common tapering functions used in STFT-based *α(f)* estimations and offer good frequency resolution.^{28} Other tapering functions with lower side-lobe levels (Blackman-Harris, Chebyshev, and Kaiser) were not found to improve the accuracy and precision of the estimates over the Hann function. Assuming an underlying white Gaussian noise process, the 50%-overlap, Hann-tapered Welch's method is expected to reduce the variance of the power spectral density estimate by a factor k_{Weclh}= 18N_{SS}^{2}/(19N_{SS}-1).^{29} Thus, k_{Welch}= 2.89 for our case. Under the same assumptions, the multitaper method would reduce the spectral variance by a factor k_{Multitaper}~N_{DPSS}=5,^{29,60} although this factor varies with the power spectral density value. When applied to RF echo signals from tissues or phantom materials, noise properties of the RF echo signals differing from white noise, such as colored noise from the pulse-tissue spectrum and coherent noise from the interference of scattered waves, affect the variance reduction performance of these methods.^{29}

The multitaper method allows the user to conveniently choose the time-bandwidth product (*NW*) for the windowing function, where *N* is the axial length of the spectral estimation window and *W* is half the main-lobe bandwidth of the frequency response of the multitaper method, *W*_{Multitaper}, normalized to the sampling frequency *f*_{s} (*W*=*W*_{Multitaper}/*f*_{s}). It thus expresses the main-lobe width of the multitaper method as a multiple *k* of the main-lobe width of the STFT method with a rectangular taper (*W*_{Rect}=1/*N*Δ*t*), where Δ*t* =1/*f*_{s} is the sampling period:

$$NW=N\frac{{W}_{\text{Multitaper}}}{{f}_{s}}=N\Delta t{W}_{\text{Multitaper}}=N\Delta tk{W}_{\mathrm{Rect}}=N\Delta tk\frac{1}{N\Delta t}k$$

(3)

As previously mentioned, the variance of the power spectral density estimate can be reduced by using a larger number of tapers, and this number is related to the selected time- half bandwidth product (*N _{DPSS}*=2

Attenuation estimates were made using power spectral density estimates from the sample and from the reference phantom computed using various window sizes. Three non-overlapped parameter estimation regions were defined side by side in each RF data frame. Each one of these regions included 22 uncorrelated scanlines. Estimates of the power spectral density were computed within each parameter estimation region at 34 axial locations, separated by one pulse-echo correlation length. This axial span accommodated at least four uncorrelated power spectral density estimates (see Discussion). The length and width of the spectral estimation windows centered at these locations, as well as the spectral estimation methods, were varied (see above), and attenuation estimates were obtained for each combination of these parameters. Reference spectra were obtained at the same axial locations as in the sample. To reduce spectral coherence noise in the reference spectra, a total of 4225 estimates from 25 parallel frames were averaged.

The sample-to-reference power spectral density ratio at each spectral estimation window depth was obtained. Assuming negligible differences in attenuation losses over the spectral estimation window length, this is given by:^{13}

$$R(f;z)\frac{PS{D}_{\text{sample}}(f;z)}{PS{D}_{\mathrm{ref}}(f;z)}=\frac{{\sigma}_{b,\text{sample}}(f;z)}{{\sigma}_{b,\mathrm{ref}}(f;z)}A(f;z),$$

(4)

where

$$A(f;z)=\mathrm{exp}\left[-4{\int}_{0}^{z}({\alpha}_{\text{sample}}(f;z)-{\alpha}_{\mathrm{ref}}(f;z))dz\right],$$

*σ _{b}(f,z)* is the backscatter coefficient,

Assuming that the backscatter and attenuation coefficients do not vary over the parameter estimation region the logarithm *L(f ;z)* of the power spectral density ratio can be expressed as:^{13}

$$\begin{array}{c}\hfill L(f;z)\mathrm{ln}\left[R(f;z)\right],\hfill & \hfill L(f;z)=\mathrm{ln}\left(\frac{{\sigma}_{b,\text{sample}}\left(f\right)}{{\sigma}_{b,\mathrm{ref}}\left(f\right)}\right)-4[{\alpha}_{\text{sample}}\left(f\right)-{\alpha}_{\mathrm{ref}}\left(f\right)]z.\hfill \end{array}$$

(5)

A least squares linear fit is applied to *L(f ;z)* vs. *z.* The slope of this fit is proportional to the difference between the attenuation coefficients of the sample and the reference phantom at frequency *f*. The latter is known so the attenuation coefficient for the sample at *f* is derived. Three estimates of *α*_{sample}*(f)* were available per frame, one for each parameter estimation region. Their average was taken, and a power-law least-squares fit (vs. frequency) was applied to produce one estimate of *α*_{0} and *β* per frame [see Eq. (1)].

In addition to the effects of reducing the spectral estimation window size, we also analyzed the accuracy and precision of *α*_{0} and *β* estimates when shortening the distance over which the linear fit of *L(f ;z)* vs. depth was applied. In figure 1 this distance is referred to as the “α-estimation length”. This was achieved by gradually reducing the number of power spectral density estimates within the parameter estimation region from 34 down to 2.

Two characteristics of the power spectral density estimate that affect the accuracy and precision of QUS parameters are its frequency resolution, and variance. These were quantified in terms of Δ*f,* the spectral correlation bandwidth, and *SNR*_{PSD}, the signal to noise ratio of the power spectral density estimate. Δ*f* is defined by Lizzi *et al.*^{22} as the frequency band over which the spectral correlation coefficient varies by no more than 70%. The power spectral density correlation coefficient *C _{PSD}* is defined as:

$${C}_{PSD}(f,f+\Delta f)=\frac{\Delta PSD\left(f\right)\Delta PSD(f+\Delta f)\sqrt{{\Delta PSD\left(f\right)2}^{}}}{}$$

(6)

where

$$\Delta PSD\left(f\right)=PSD\left(f\right)-PSD\left(f\right).$$

The value of Δ*f* was assessed using the power spectral density estimates at the deepest locations (34^{th} row of the 34×3 array) because they present the most severe bandwidth constraints caused by attenuation. The values of *C*_{PSD}*(f, f+*Δ*f)* were averaged across all frequency bins within the analysis bandwidth previously defined (4.8-8.0MHz). The average and standard deviation of Δ*f* were computed from each of the three deepest *PSD*_{sample}*(f)* estimates in each of the 20 sample frames.

The value of *SNR*_{PSD} at frequency *f* and depth *z* was computed as

$$SN{R}_{\mathrm{PSD}}(f,z)=\frac{PS{D}_{sample,r}(f,z)\sqrt{{[PS{D}_{sample,r}(f,z)-PS{D}_{sample,r}(f,z)]2}^{}}}{}$$

(7)

where *r* refers to each estimate of *PSD(f,z)* from the 20 sample frames. Brackets indicate averaging across the 20 estimates. For simplicity, we report the values of *SNR _{PSD}(f,z)* at the frequency of the maximum power spectral density value, i.e., 6.2MHz. The average and standard deviation of

The accuracy and the precision of *α*_{0} and *β* estimates as independent QUS parameters were quantified when constraining the spectral estimation window size and the *α*-estimation length. For this purpose, the normalized bias and standard deviation of *α*_{0} and *β* estimates were computed. The normalized bias (Bias* _{Q}*) was defined as the expected value of the difference between estimated quantity

$${\mathrm{Bias}}_{Q}=\frac{1}{{N}_{R}}\sum _{r=1}^{{N}_{R}}\left[\frac{{Q}_{r}-{Q}_{\mathrm{exp}}}{{Q}_{\mathrm{exp}}}\right]$$

(8)

where the subindex *r* refers to each of the *N _{R}*=20 sample frames. The normalized standard deviation (STD

$${\mathrm{STD}}_{Q}=\sqrt{\frac{1}{{N}_{R}}\sum _{r=1}^{{N}_{R}}{[\frac{{Q}_{r}-{Q}_{r}{Q}_{\mathrm{exp}}]}{2}}^{}.}$$

(9)

We also analyzed an overall “fit error”, *FE*, defined as the expected value of the absolute difference between the estimated power-law fit *α*_{Fit,r}*(f)* and the laboratory-estimated power law fit α_{exp}*(f)*, normalized with respect to the latter at each frequency point, and averaged across all frequency points included in the analysis bandwidth. It was estimated by computing for each frame *r:*

$$F{E}_{r}=\frac{1}{{N}_{F}}\sum _{i=1}^{{N}_{F}}[\frac{{\alpha}_{Fit,r}\left({f}_{i}\right)-{\alpha}_{\mathrm{exp}}\left({f}_{i}\right){\alpha}_{\mathrm{exp}}\left({f}_{i}\right)]}{}$$

(10)

and then obtaining the average and standard deviation across the 20 sample frames.

Error bars in plots of Bias* _{Q}*, STD

RESULTS

Figure 3(a) shows the frequency resolution measured by the spectral correlation bandwidth Δ*f* as a function of the power spectral density estimation window length in pulse-echo (axial) correlation lengths. Each curve corresponds to a different spectral estimation method. Shortening the window increases Δ*f*, as expected. Furthermore, Δ*f* for the short-time Fourier transform (STFT) with a Hann taper, Welch's periodogram, and Thomson's multitaper method was larger than for the rectangular-taper STFT (or no taper) by a factor of 1.5, 2.2, and 2.8, respectively, because of the wider frequency-response main lobe of the latter spectral estimation methods with respect to the former. Thus, the STFT methods lead to better frequency resolution than the Welch and the multitaper methods.

The *SNR*_{PSD} is shown in figure 3(b) as a function of the number of uncorrelated scanline signals within the spectral estimation window. A guide curve corresponding to the square root of the number of uncorrelated scanlines is also provided. *SNR*_{PSD}=1 when no spectral averaging is performed and increases with N^{1/2} when averaging N estimates from uncorrelated RF echo scanlines.^{63,64} The agreement with the N^{1/2} curve corroborates our estimation of the lateral correlation length. *SNR*_{PSD} is also expected to increase by factors k_{Welch} and k_{Multitaper} (defined in the Methods section) when using the Welch and multitaper methods. From figure 3(b), k_{Welch}=2.7 and k_{Multitaper}=4.5. These values are slightly smaller than those predicted for a white Gaussian noise process because of the particular noise properties of the RF echo signals and the use of Thomson's version of the multitaper method.^{29}

Figure 4 shows power spectral density estimates of the sample from all spectral estimation methods using either (a, left) a severely length-constrained but wide spectral-analysis window (2.5 pulse-echo correlation lengths, 22 uncorrelated scanlines), or (b, right) a severely width-constrained spectral estimation window with its length in the middle of the analyzed range (10 pulse-echo correlation lengths, 4 uncorrelated scanlines). In the former case, the Welch and multitaper methods have wider power spectral density estimates than the STFT-based methods. In the latter case, the smooth appearance of the multitaper estimate is due to the ability of this method to reduce coherent noise, as quantified in figure 3(b).

Figure 5 shows *α*_{sample}*(f)* estimates plotted versus frequency for an *α*-estimation length of 33 pulse-echo correlation lengths (12.4mm). Results are divided into two rows (i and ii) and three columns (a-c). The rectangular-taper STFT was used for spectral estimation in the top (i) row, while the multitaper was used in the bottom (ii) row. Columns (a), (b), and (c) correspond to different window sizes: 20 pulse-echo correlation lengths and 22 uncorrelated scanlines (largest window), 10 pulse-echo correlation lengths and 22 uncorrelated scanlines, and 10 pulse-echo correlation lengths and 2 uncorrelated scanlines, respectively. In each subfigure, the **dark** gray continuous curve is the estimated *α*_{sample}*(f)* averaged across 20 frames, the light gray continuous curves are the average ± one standard deviation, and the black curve is the laboratory estimate. The *α*_{sample}*(f)* estimates are centered about the laboratory estimates, indicating their accuracy. Furthermore, the standard deviation of the estimated *α*_{sample}*(f)* does not depend on frequency, which agrees with theoretical derivations.^{27,63,64}

Laboratory and “reference phantom method” estimates of *α*_{sample}*(f)* [dB/cm] vs. frequency for (a) the largest window size (20 pulse-echo correlation lengths (PECLs) and 22 uncorrelated acoustic scanlines (SL)), (b) 10 PECLs and 22 **...**

The reduction of the spectral coherence noise by the multitaper method leads to a less fluctuating frequency dependence of the estimated *α*_{sample}*(f)* (figure 5, bottom row), compared to the use of the STFT method (figure 5, top row). The reduction of the spectral estimation window length also reduces the spectral coherence noise, although a slight increase in variance expressed as an increase in the overall separation between the lightly shaded curves (average ± one standard deviation) is observed. This can be attributed to the loss of backscatter signal information as the window length is reduced.^{29}

Figures 6 and and77 show the normalized bias as defined in Eq. (8) [*α*_{0} in figure 6(a), *β* in figure 7(a)] and the normalized standard deviation as defined in Eq. (9) [*α*_{0} in figure 6(b), *β*] in figure 7(b)] as functions of the *α*-estimation length expressed in pulse-echo correlation lengths. Subfigures (a) and (b) have four components (i-iv), one for each spectral estimation method. Different curves in each component correspond to a different spectral estimation window length. In all cases, the number of uncorrelated scanlines used for spectral estimation was 22. To facilitate image interpretation, error bars for window lengths only of 20 and 7.5 pulse-echo correlation lengths at two different *α*-estimation length values are included. Figure 8 shows the overall fit error defined in Eq. (10), and subfigures (a)-(d) correspond to different spectral estimators.

(a) Normalized bias and (b) normalized standard deviation of the estimates of *α*_{0} vs. the *α*-estimation length. Each curve corresponds to a different spectral estimation window length, but a fixed width of 22 uncorrelated acoustic scanlines. **...**

(a) Normalized bias and (b) normalized standard deviation of the estimates of *β* vs. the *α*-estimation length. Each curve corresponds to a different spectral estimation window length, but a fixed width of 22 uncorrelated acoustic scanlines. **...**

Normalized fit error of the power-law fit applied to *α*_{sample}(*f*) with respect to the laboratory-estimated power-law model vs the *α*-estimation length. Each curve corresponds to a different spectral estimation window length, but a fixed width **...**

In figures 6--8,8, combinations of spectral estimation window lengths and *α*-estimation lengths shorter than 5 pulse-echo correlation lengths led to errors in *α _{0}* and

In some cases, an apparent decrease in fit error when using *α*-estimation blocks shorter than ~5 pulse-echo correlation lengths combined with spectral estimation window lengths shorter than 7.5 pulse-echo correlation lenghts was observed. This might have been caused by the inability to measure a signal decrease over such short data lengths. Under these circumstances, the estimate of the slope in Eq. (5) is effectively zero, and therefore the value assigned to *α*_{sample}*(f)* is that of *α*_{ref}*(f)*, which happened to be similar to the expected sample attenuation in this experiment.

Figure 9 depicts (a) the normalized bias and (b) the normalized standard deviation in estimates of (i) *α _{0}* and (ii)

(a) Normalized bias and (b) normalized standard deviation of the estimates of (i) *α*_{0} and (ii) *β* vs. the spectral estimation window width expressed in the number of uncorrelated acoustic scanlines, with a fixed window length of 10 pulse-echo **...**

The problem of constraints in the size of the parameter-estimation region on *α(f)* estimates was addressed by Parker,^{27} Yao *et al.,*^{63} and Huisman *et al.*^{64} These authors show that the variance of the slope of a linear fit of *α* versus frequency is proportional to the ratio of the square root of the frequency resolution divided by the signal-to-noise ratio of the power spectral density estimate (*SNR*_{PSD}). Considering this and the results presented in figures 3(a) and 3(b), the multitaper method is expected to offer the best precision for attenuation coefficient estimations from spectral data. However, theoretical expressions for the variance of descriptors of the attenuation coefficient^{27,63,64} cannot be directly applied to our experiment because they were derived assuming the estimates of the power spectral density are normally distributed random variables. As discussed by Lizzi *et al.*,^{22} this assumption does not apply to severely constrained parameter estimation regions like some of those considered here.

Averaging periodograms from adjacent scanlines was used for computing the sample and reference power spectral density.^{22} When a small number of independent realizations of the periodogram are available, the Welch and multitaper methods notably reduced the contribution of spectral coherence noise. The increase in *SNR*_{PSD} for the Welch and multitaper methods, with respect to the rectangular-taper STFT method (k_{Welch} and k_{Multitaper}), was close to theoretical predictions for white Gaussian noise. Thus the multitaper method outperformed the other spectral estimation methods analyzed in terms of *SNR*_{PSD}.

Average estimates of *α*_{sample}*(f)* shown in figure 5 were in agreement with laboratory estimates done on test samples of the phantom material. However, it is important to note that an estimate of the parameters of the power-law model applied to *α*_{sample}*(f)* is obtained from a single estimate of *α*_{sample}*(f)* and not from the inter-frame averaged *α*_{sample}*(f)* presented in figure 5. Assuming a Gaussian distribution,^{63} a single estimate of *α*_{sample}*(f)* is expected to lie between the two “average ± one standard deviation” curves in each subfigure of figure 5 with a probability of 68%. Although this assumption does not hold when the number of uncorrelated scanlines used in spectral estimation is severely constrained, it gives an idea of the expected variability of *α*_{0} and *β* estimates.

Figures 6--1010 show that the shortest reliable spectral estimation window length was about 8 pulse-echo correlation lengths. The pulse-echo correlation lengths estimated from intensity (squared RF) data can be considered an estimate of the pulse length.^{65,66} The intensity pulse-echo correlation length was 1.47 RF pulse-echo correlation lengths in our experiment. Therefore, 8 RF pulse-echo correlation lengths ~ 11 pulse lengths. This is within the range of 7-15 pulse lengths that is frequently presented as the window length offering the best tradeoff between spatial resolution, accuracy and precision when estimating attenuation and effective scatterer size.^{20,24-26,31,36,64-69} Tapering functions such as the Hann or Hamming, commonly used in QUS,^{14,15,24,25,36,62,67,70} were found to increase the bias and standard deviation of the attenuation estimates. This indicates a less important role of side-lobe reduction in the accuracy and precision of the estimated parameters than the frequency resolution and SNR of the power spectral density estimate.

The bias and standard deviation of both *α _{0}* and

When constraining the lateral dimension of the spectral estimation window (figures 9 and and10),10), our results show that the standard deviation of the *α*_{0} and *β* estimates are significantly larger than the bias. In the results shown, the window length was already reduced to 10 pulse-echo correlation lengths and the *α*-estimation length was 33 pulse-echo correlation lengths. According to the discussion in the previous paragraph, these parameter estimation region length values are within the ranges where the *α*_{0} and *β* estimates are not expected to vary significantly. The average standard deviations of *α*_{0} and *β* in these ranges, shown in figures 6(b) and 7(b), were about 37% and 16%. To keep these standard deviation levels, figures 9(b.i) and 9(b.ii) show that the window width can be reduced to about 15 uncorrelated scanlines. This is within the range of 10-20 independent scanlines suggested elsewhere.^{22,25,26,32} For our conditions, this parameter estimation region width spans 7.4mm. This is smaller than the parameter estimation region length of ~9mm stated in the previous paragraph. These results agree with previous findings that the total area spanned by the parameter estimation region would be limited by errors caused by its axial length rather than those caused by its lateral width.^{26,27,63}

If tissue homogeneity allows the selection of a parameter estimation region that can axially accommodate at least two uncorrelated power spectral density estimates and at least 15 uncorrelated scanlines, all spectral estimation methods shown here are expected to perform similarly. This parameter estimation region size (15 pulse-echo correlation lengths × 15 uncorrelated scanlines) fits the task of bulk parameter estimation where large homogeneous parameter estimation regions are available to obtain accurate and precise estimates of *α*_{0} and *β*.

The diagnostic performance of a parametric image is related to the contrast-to-noise ratio (CNR) exhibited by the structure of interest. The CNR is affected by the loss of accuracy and precision of the estimates when the parameter estimation region is gradually reduced below ~1cm × 1cm to create parametric images with sufficient spatial resolution to depict spatial fluctuations of the tissue property. In the case of *α*_{0}, the minimum *α*-estimation and spectral estimation window lengths for which the bias is less than 10% are ~10 and ~8 pulse-echo correlation lengths, respectively, spanning a total length of 18 pulse-echo correlation lengths or ~6.8mm (considering a parameter estimation region width of 22 uncorrelated scanlines or 10.9mm). This parameter estimation region length would allow, for example, to distinguish fatty tissue from invasive ductal carcinoma in the breast, whose *α _{0}* values differ from healthy parenchyma by about 72% and 53%, respectively.

The advantage of the multitaper method has been attributed to its ability to reduce spectral coherence noise without severely compromising the accuracy of the power spectral density estimates. The influence of spectral coherence noise may be more important when the analysis bandwidth is further reduced, such as when dealing with low SNR values caused by strongly attenuating or deep lesions. In these scenarios, the multitaper method would offer additional advantages compared to other estimation methods. As an example, in figure 11 the overall fit error is plotted versus the spectral estimation window width for the rectangular-taper STFT (gray curves) and multitaper (black curves) methods. As observed, the multitaper method maintains the same error level under the bandwidth reductions tested, while the error from rectangular-taper STFT significantly increased when the bandwidth was reduced. Another method to reduce spectral coherence noise is spectral smoothing, achieved by convolving the power spectral density estimate with a smoothing function, such as that applied by Kim *et al.*^{19} for estimating attenuation. However, spectral smoothing introduces bias that would affect the smooth spectral components that the multitaper method was designed to resolve by means of the *NW* selection criterion.^{29} Alternatives to spectral estimation based on the Fourier transform include autoregressive algorithms. A comprehensive comparison between the multitaper method and the autoregressive methods is beyond the scope of this manuscript and will be a topic of future research.

This work has investigated the performance of various methods to compute the power spectral density in the context of estimating parameters describing the frequency dependence of the acoustic attenuation coefficient. For the task of obtaining bulk estimates of these parameters, our results support the use of any of the power spectral density estimations tested, with expected error levels of the power-law fit to *α(f)* on the order of 4%. When creating parametric images, the selection of the size of the parameter estimation region must consider its strong influence on the bias and standard deviation of the *α*_{0} and *β* estimates. However, in general, the multitaper method reduced the standard deviation of the parameter estimates with respect to the other methods (21% reduction for a parameter estimation region offering a 10% bias in *α*_{0}). These results, which are specific to the materials investigated here, encourage use of Thomson's multitaper method for estimating acoustic parameters under spatially constrained conditions.

1. Insana MF, Oelze ML. Advanced ultrasonic imaging techniques for breast cancer research. In: Suri JS, Rangayyan RM, Laxminarayan S, editors. Emerging technologies in breast imaging and mammography. American Scientific Publishers; Valencia, CA: 2006.

2. Thijssen JM. Ultrasonic speckle formation, analysis and processing applied to tissue characterization. Pattern Recognition Letters. 2003;24:659–675.

3. Zagzebski JA. Essentials of Ultrasound Physics. Mosby; St. Louis, MO: 1996.

4. Chen JF, Zagzebski JA, Madsen EL. Test of backscatter coefficient measurement using broadband pulses. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 1993:40, 603–607. [PubMed]

5. Campbell JA, Waag RC. Measurements of calf liver ultrasonic differential and total scattering cross sections. J. Acoust. Soc. Am. 1984;75:603–611.

6. Lizzi FL, King DL, Rorke MC, Hui J, Ostromogilsky M, Yaremko MM, Feleppa EJ, Wai P. Comparison of theoretical scattering results and ultrasonic data from clinical liver examinations. Ultrasound Med. Biol. 1988;1988:377–385. [PubMed]

7. Insana MF, Hall TJ, Fischback JL. Identifying acoustic scattering sources in normal renal parenchyma from the anisotropy in acoustic properties. Ultrasound Med. Biol. 1991;17:613–626. [PubMed]

8. Insana MF, Wood JG, Hall TJ. Identifying acoustic scattering sources in normal renal parenchyma in vitro by varying arterial and ureteral pressures. Ultrasound Med. Biol. 1992;18:587–599. [PubMed]

9. Oelze ML, Zachary JF, O'Brien WD. Parametric imaging of rat mammary tumors in vivo for the purpose of tissue characterization. J. Ultrasound Med. 2002;21:1201–1210. [PubMed]

10. Wilson T, Chen Q, Zagzebski JA, Varghese T, Van Middlesworth L. Initial clinical experience imaging scatterer size and strain in thyroid nodules. J. Ultrasound Med. 2006;25:1021–1029. [PubMed]

11. Insana MF, Hall TJ. Characterizing the microstructure of random media using ultrasound. Phys. Med. Biol. 1990;35:1373–1386. [PubMed]

12. Insana MF, Wagner RF, Brown DG, Hall TJ. Describing small-scale structure in random media using pulse-echo ultrasound. J. Acoust. Soc. Am. 1990;87:179–192. [PMC free article] [PubMed]

13. Yao LX, Zagzebski JA, Madsen EL. Backscatter coefficient measurements using a reference phantom to extract depth-dependent instrumentation factors. Ultrasonic Imaging. 1990;12:58–79. [PubMed]

14. Oelze ML, O'Brien WD. Method for improved scatterer size estimation and application to parametric imaging using ultrasound. J. Acoust. Soc. Am. 2002;112:3053–3063. [PubMed]

15. Gerig A, Zagzebski JA, Varghese T. Statistics of ultrasonic scatterer size estimation with a reference phantom. J. Acoust. Soc. Am. 2003;113:3430–3437. [PubMed]

16. Bigelow TA, O'Brien WD. Scatterer size estimation in pulse-echo ultrasound using focused sources: Theoretical approximations and simulation analysis. J. Acoust. Soc. Am. 2004;116:578–593. [PubMed]

17. Liu W, Zagzebski JA, Hall TJ, Madsen EL, Varghese T, Kliewer MA, Panda S, Lowery C, Barnes S. Acoustic backscatter and effective scatterer size estimates using a 2D CMUT transducer. Phys. Med. Biol. 2008;53:4169–4183. [PMC free article] [PubMed]

18. Bigelow TA, O'Brien WD. Impact of local attenuation approximations when estimating correlation length from backscattered ultrasound echoes. J. Acoust. Soc. Am. 2006;120:546–552. [PubMed]

19. Kim H, Varghese T. Hybrid spectral domain method for attenuation slope estimation. Ultrasound Med. Biol. 2008;34:1808–1819. [PubMed]

20. Labyed Y, Bigelow TA. Estimating the total ultrasound attenuation along the propagation path by using a reference phantom. J. Acoust. Soc. Am. 2010;128:3232–3238. [PubMed]

21. Heo SW, Kim H. A novel power spectral density calculation method using phase-compensation and weighted averaging for the estimation of ultrasound attenuation. Ultrasonics. 2010;50:592–599. [PubMed]

22. Lizzi FL, Alam K, Mikaelian S, Lee P, Feleppa EJ. On the statistics of ultrasonic spectral parameters. Ultrasound Med. Biol. 2006;32:1671–1685. [PubMed]

23. Wear KA, Wagner RF, Garra BS. High resolution ultrasonic backscatter coefficient estimation based on autoregressive spectral estimation using Burg's algorithm. IEEE Trans. Med. Imag. 1994;13:500–507. [PubMed]

24. Oelze ML, O'Brien WD. Improved scatterer property estimates from ultrasound backscatter for small gate lengths using a gate-edge correction factor. J. Acoust. Soc. Am. 2004;116:3212–3223. [PubMed]

25. Oelze ML, O'Brien WD. Defining optimal axial and lateral resolution for estimating scatterer properties from volumes using ultrasound backscatter. J. Acoust. Soc. Am. 2004;115:3226–3234. [PubMed]

26. Liu W, Zagzebski JA. Trade-offs in data acquisition and processing parameters for backscatter and scatterer size estimations. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2010;57:340–341. [PMC free article] [PubMed]

27. Parker KJ. Attenuation measurement uncertainities caused by speckle statistics. J. Acoust. Soc. Am. 1986;80:727–734. [PubMed]

28. Harris FJ. On the use of windows for harmonic analysis with the Discrete Fourier Transform. Proc. IEEE. 1978;66:51–83.

29. Percival DB, Walden AT. Spectral Analysis for Physics Applications: Multitaper and Conventional Univariate Techniques. Cambridge University Press; New York, NY: 1993.

30. Bronez TP. On the performance advantage of multitaper spectral analysis. IEEE Trans. Signal Process. 1992;40:2941–2946.

31. Labyed Y, Bigelow TA. A theoretical comparison of attenuation measurement techniques from backscattered ultrasound echoes. J. Acoust. Soc. Am. 2011;129:2316–2324. [PubMed]

32. Bigelow TA, Labyed Y, McFarlin BL, Sen-Gupta E, O'Brien WD. Comparison of algorithms for estimating ultrasound attenuation when predicting cervical remodeling in a rat model.. Biomedical Imaging: From Nano to Macro, 2011 IEEE International Symposium 2011.pp. 883–886.

33. Labyed Y, Bigelow TA, McFarlin BL. Estimate of the attenuation coefficient using a clinical array transducer for the detection of cervical ripening in human pregnancy. Ultrasonics. 2011;51:34–39. [PMC free article] [PubMed]

34. Bigelow TA. Estimation of ultrasound tissue attenuation along the propagation path by applying multiple filters to the backscattered echoes. Ultrasonics Symposium. Proc. IEEE. 2009;2009:1954–1957.

35. Bigelow TA. Estimating the total ultrasound attenuation along the propagation path by applying multiple filters to backscattered echoes from a single spherically focused source. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2010;57:900–907. [PMC free article] [PubMed]

36. Zagzebski JA, Liu ZF, Yao LX. Quantitative ultrasound imaging: In vivo results in normal liver. Ultrasonic Imaging. 1993;15:335–351. [PubMed]

37. Goss SA, Frizzel LA, Dunn F. Ultrasonic absorption and attenuation in mammalian tissues. Ultrasound Med. Biol. 1979;5:181–186. [PubMed]

38. Narayana PA, Ophir J. On the frequency dependence of attenuation in normal and fatty liver. IEEE Trans. Sonics. Ultrason. 1983;SU-30:379–383.

39. Narayana PA, Ophir J. The measurement of attenuation in nonlinearly attenuation media by the zero crossing method. Ultrasound Med. Biol. 1984;10:715–718. [PubMed]

40. Ophir J, McWhirt RE, Jaeger PM. A narrowband pulse-echo technique for in vivo ultrasonic attenuation estimation. IEEE Trans. Biomed. Eng. 1985;32:205–212. [PubMed]

41. Lin T, Ophir J, Potter G. Correlation of ultrasonic attenuation with pathologic fat and fibrosis in liver disease. Ultrasound Med. Biol. 1988;14:729–734. [PubMed]

42. D'Astous FT, Foster FS. Frequency dependence of ultrasound attenuation and backscatter in breast tissue. Ultrasound Med. Biol. 1986;12:795–808. [PubMed]

43. Cloostermans MJ, Verhoef WA, Thijssen JM. Generalized description and tracking estimation of the frequency dependent attenuation of ultrasound in biological tissues. Ultrasonic Imaging. 1985;7:133–141. [PubMed]

44. Daft CMW, Briggs GAD, O'Brien WD. Frequency dependence of tissue attenuation measured by acoustic microscopy. Ultrasonics Symposium. 1988. Proc. IEEE. 1988;2:971–974.

45. Anderson JJ, Herd MT, King MR, Haak A, Hafez ZT, Song J, Oelze ML, Madsen EL, Zagzebski JA, O'Brien WD. Interlaboratory comparison of backscatter coefficient estimates for tissue-mimicking phantoms. Ultrasonic Imaging. 2010;32:48–64. [PMC free article] [PubMed]

46. Madsen EL, Hobson AM, Shi H, Varghese T, Frank GR. Stability of heterogeneous elastography phantoms made from oil dispersions in aqueous gels. Ultrasound Med. Biol. 2006;32:261–270. [PMC free article] [PubMed]

47. Nam K, Rosado-Mendez IM, Rubert NC, Madsen EL, Zagzebski JA, Hall TJ. Effects of sound speed mismatch on ultrasound attenuation measurements using a reference phantom. Ultrasonic Imaging. 2011;33:251–263. [PMC free article] [PubMed]

48. Brunke SS, Insana MF, Dahl JJ, Hansen C, Ashfaq M, Emert H. An ultrasound research interface for a clinical system. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2007;54:198–210. [PubMed]

49. Rosado-Mendez IM, Madsen EL, Hall TJ, Zagzebski JA. Susceptibility of backscatter coefficient estimates to nonlinear wave propagation. Annual Convention of the American Institute of Ultrasound in Medicine; Phoenix, AZ: 2012.

50. O'Donnel M, Silverstein SD. Optimum displacement for compound image generation in medical ultrasound. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 1998;35:470–476. [PubMed]

51. Trahey GE, Smith SW, von Ramm OT. Speckle pattern correlation with lateral aperture translation: Experimental results and implications for spatial compounding. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 1986;33:257–264. [PubMed]

52. Hall TJ, Insana MF, Harrison LA, Cox GG. Ultrasonic measurement of glomerular diameters in normal adult humans. Ultrasound Med. Biol. 1996;22:987–997. [PubMed]

53. Herd MT, Hall TJ, Jiang J, Zagzebski JA. Improving the statistics of quantitative ultrasound techniques with deformation compounding: an experimental study. Ultrasound Med. Biol. 2011;37:2066–2074. [PMC free article] [PubMed]

54. Welch PD. 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. 1967;15:70–73.

55. Kim H, Varghese T. Attenuation estimation using spectral cross-correlation. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2007;54:510–519. [PubMed]

56. Donohue KD, Forsberg F, Piccoli CW, Goldberg BB. Analysis and classification of tissue with scatterer structure templates. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 1999;46:300–310. [PubMed]

57. Varghese T, Donohue KD. Estimation of mean scatterer spacing with spectral correlation. J. Acoust. Soc. Am. 1994;96:3504–3515. [PubMed]

58. Goshal G, Oelze ML. Improved scatterer property estimates form ultrasound backscatter using gate-edge correction in a pseudo-Welch technique. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2010;57:2828–2832. [PMC free article] [PubMed]

59. Stoica P, Moses RL. Introduction to Spectral Analysis. Prentice Hall; New Jersey, NJ: 1997.

60. Thomson DJ. Spectral density estimation and harmonic analysis. Proc. IEEE. 1982;70:1055–1082.

61. Gruenbacher DM, Hummels DR. A simple algorithm for generating discrete prolate spheroidal sequences. IEEE Trans. Signal Process. 1994;42:3276–3278.

62. Oelze ML, O'Brien WD. Frequency-dependent attenuation-compensation functions for ultrasonic signals backscatterered from random media. J. Acoust. Soc. Am. 2002;11:2308–2319. [PubMed]

63. Yao LX, Zagzebski JA, Madsen EL. Statistical uncertainty in ultrasonic backscatter and attenuation coefficients determined with a reference phantom. Ultrasound Med. Biol. 1991;17:187–194. [PubMed]

64. Huisman JH, Thijssen TM. Precision and accuracy of acoustospectrographic parameters. Ultrasound Med. Biol. 1996;22:855–871. [PubMed]

65. Wagner RF, Insana MF, Brown DG. Statistical properties of radio-frequency and envelope-detected signals with applications to medical ultrasound. J. Opt. Soc. Am. 1987;4:910–922. [PubMed]

66. Wagner RF, Insana MF, Smith SW. Fundamental Correlation Lengths of Coherent Speckle in Medical Ultrasonic Images. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2002;35:34–44. [PubMed]

67. Bigelow TA, Oelze ML, O'Brien WD. Estimation of total attenuation and scatterer size from backscattered ultrasound waveforms. J. Acoust. Soc. Am. 2005;117:1431–1429. [PubMed]

68. Bigelow TA. Ultrasound attenuation estimation using backscattered echoes from multiple sources. J. Acoust. Soc. Am. 2008;124:1367–1373. [PubMed]

69. Bigelow TA. Improved algorithm for estimation of attenuation along propagation path using backscattered echoes from multiple sources. Ultrasonics. 2010;50:496–501. [PMC free article] [PubMed]

70. Wear KA. Fundamental precision limitations for measurements of frequency dependence of backscatter: Application in tissue-mimicking phantoms and trabecular bone. J. Acoust. Soc. Am. 2001;110:3275–3282. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |