Comparative study of periodicity estimation methods using ultrasonic signals

Introduction: Various signal-processing techniques have been proposed to extract quantitative information about internal structures of tissues from the original radio frequency (RF) signals instead of an ultrasound image. The quantifiable parameter called the mean scatterer spacing (MSS) can be useful to detect changes in the quasi-periodic microstructure of tissues such as the liver or the spleen, using ultrasonic signals. Methods: We evaluate and compare the performance of three classic methods of spectral estimation to calculate the MSS without operator intervention: Tufts-Kumaresan, SAC (Spectral Autocorrelation) and MUSIC (MUltiple SIgnal Classification). Initially the evaluations were performed with 10,000 signals simulated from a model in which the variables of interest are controlled, and then, real signals from sponge phantoms were used. Results: For the simulated signals, the performance of all three methods decreased with increasing Ad or jitter levels. For the sponges, none of the methods accurately estimated the pore size. Conclusion: For the simulated signals, Tufts-Kumaresan had the lowest performance, whereas SAC and MUSIC had similar results. For sponges, only Tufts-Kumaresan was able to detect the increase in the size of the pores of the sponge, although most often, it estimated sizes larger than expected.


Introduction
Ultrasound is a well-established method of medical imaging, which uses non-ionizing radiation and has a low cost of implementation.However, the ultrasound image is usually interpreted visually and qualitatively.Conventional ultrasound images (grayscale) are formed from echoes that return from the irradiated structures.These images represent specular reflections from interfaces between different acoustic impedance structures and backscattered echoes from the interaction between the ultrasound and microstructures.Specular reflections indicate the contours of organs, vessels, cysts, and macroscopic structures in general, while scattering generates the gray granulation (speckle) that fills the interiors of organs and tissues.These images provide a visualization of anatomical structures in motion in real time, but the more subtle changes in tissue microstructure, such as the early stages of cirrhosis or fibrosis or fetal maturation are difficult to identify visually (Shung and Thieme, 1993).
MSS represents the spatial organization and is therefore of particular interest for detecting changes in tissue that have quasi-periodic microstructures, such as the liver or spleen tissue (Fellingham and Sommer, 1984).To this end, various techniques of power spectrum and cepstrum (power spectrum of the logarithm of the power spectrum) several techniques have been proposed for estimating MSS such as the ones described by: Fellingham and Sommer (1984), Narayanan et al. (1997) and Wear et al. (1993).The studies of Varghese and Donohue (1993, 1994, 1995) were based on spectral autocorrelation, and Simon et al. (1997) proposed an algorithm to estimate the MSS based on the spectral redundancy generated by an RF signal quadratic transformation.All these techniques are applied to RF signals or their envelope, which contains a mixture of contributions from periodic and non-periodic scatterers.
Using human calcanei, Pereira et al. ( 2004) estimated MSS using Singular Spectral Analysis and concluded it may be useful for providing information associated with tissue microarchitecture.In 2006  Machado et al. (2006) studied in vitro the periodicity of human liver tissue and found that MSS alone could not characterize tissues, concluding that additional parameters are needed to improve its accuracy.Also using human tissue, Huang et al. (2008) wrote an article estimating the mean spacing within trabecular bone using a simplified inverse filter-tracking algorithm and concluded it has the potential to be a method for estimating average spacing in low signal-to-noise environments with large spacing variations.
In 2011, Chen et al. (2011) proposed a time-frequency approach to estimate the MSS.They used simulated and real backscattered echoes of liver tissues in vivo and concluded that their method is potentially superior to Simon's (1997), avoiding false peaks in the spectrum.More recently, Pan et al. (2015) proposed a Golay code-based cepstrum estimation and the results show that the method is robust against noise.
Comparing the techniques, Kauati et al. ( 2012) evaluated the Burg, Wiener, and MUSIC (MUltiple SIgnal Classification) spectral estimation methods to calculate the MSS (without operator intervention).Initially the evaluation was performed using 10,000 simulated signals with the aim of studying each method's behavior using a model in which the variables of interest can be controlled.Then, the methods were applied to real signals of nylon phantoms immersed in water.The Burg method could not estimate the spacing of phantom signals, presenting results similar to the other methods only for simulated signals.The Wiener method for simulated signals was second best in terms of percentage of success.When considering the phantom signals, the MUSIC method had the best performance of all three methods.
In this work, we will compare the MUSIC method, which had the best performance in Kauati et al. (2012), with two other methods also used to estimate periodicity: Tufts-Kumaresan and SAC (Spectral Autocorrelation).

Methods
In this section, the equations and algorithms of the MUSIC, Tufts-Kumaresan and SAC methods will be presented.The study was initially performed using simulated ultrasonic (US) signals based on a linear model of US interaction with biological tissue, in order to have more control of the variables involved.Then, we used real signals from sponge phantoms immersed in water.Examples of both kinds of signals (simulated and real) will also be presented.The approaches used to evaluate the methods, the phantoms used to generate the real signals, and the signal treatments proposed in this work will be described at the end of this section.

Tufts-Kumaresan
Based on modified covariance, Tufts and Kumaresan (1982) determined a prediction filter of N samples in the forward and backward directions, and the prediction equations can be written in matrix form as: where, y(n): input signal; g: filter impulse response; L: predictor filter order; N: number of samples.
In other words, where, A: matrix formed by the input signal; h: vector formed by the input signal.Thus, g is given by where, A # : pseudoinverse of A # .The correlation matrix R can be calculated by: Where, * is the transpose of the complex conjugate.
The correlation vector is defined by and the prediction filter g can be calculated as The Tufts-Kumaresan filter is defined as The filter order is crucial to the frequency estimation performance.According to Tufts and Kumaresan (1982) the best order is defined by the window size minus the number of periodicities in the signal divided by two.
The Tufts-Kumaresan method, with values of example parameters, is presented below.The size of the window was defined according to Tufts and Kumaresan (1982) and experimentally confirmed using a few windows of different sizes and by checking which one had the best result.
• number of periodicities = 1 • order of the filter = window size -number of periodicities • overlap = half of the window size 2. Thus, for each window, we calculate the matrix A and h, according to Equation 1, for the RF signal envelope.3. The vector g will be calculated for each window as: 4. At the end, we calculate the average of g with all windows 5. From the filter g tk frequency response, we can find the frequency within the specified window that presents the largest amplitude value.This will be the requested frequency (Figure 1A).

SAC
The SAC (Spectral Autocorrelation) method performs an analysis on the spectrum calculated by: The relative frequency of the mean spacing of regular scatterers is determined by the off-diagonal dominant peak, which is equivalent to the PSD.The analysis outside the main diagonal was used, because according to Varghese and Donohue (1993), the expected value of diffuse component is zero in this region.
The steps for the implementation of SAC method can be summarized as follows: 1. Calculate the FFT of the envelope of the RF signal and its conjugate.2. Calculate the product of both of them (Equation 8). 3. Normalize through the main diagonal.4. Remove the main diagonal.5.The frequency related to the position of maximum amplitude in the spectrum is the one related to MSS (Figure 1B).

MUSIC
The MUSIC (MUltiple SIgnal Classification) method (Marple, 1987) is based on obtaining the matrix eigenvalues and eigenvectors formed by covariances according to equation: where, c (i) = matrix elements; Y = original signal.
From a generic estimator of frequencies: being, a k = eigenvector; e H = eigenvector transposed; α k = weight.
The estimator of frequencies of MUSIC is based on the noise subspace eigenvectors with uniform weight, so the frequency regarding the mean spacing of scatterers is defined by the point of maximum amplitude in the MUSIC spectrum.

( ) ( ) (
) where ( ) e f is a vector formed by: 2) where, T is half the period limited by the analysis window for M eigenvectors forming the signal space.
The steps for the MUSIC method implementation can be summarized as follows: 1. Draw the MUSIC spectrum according to Equation 11, considering the RF signal envelope.2. The frequency related to the position of maximum amplitude in the spectrum is the one related to MSS. (Figure 1C).

Simulated tissue
The RF signals simulated in this study followed the model used by Maciel (2000) of backscattered echoes, considering a one-dimensional medium and linear behavior.Thus, the echo signal received as a function of time, r(t), can be written as: where, ( ) p t = transmitted ultrasound pulse; ( ) g t = function of characterization of medium (Impulse Response); ( ) n t = experimental system noise; * = convolution.
The function g(t) can be modeled as the sum of echoes scattered by each mediumparticle where a group of particles has a regular spatial distribution (periodic) and the other group has an irregular distribution (aperiodic): where: N = total number of regular particles; M = total number of dispersed particles; i c = regular signal amplitude; i τ = regular part delay (regarding the position); i b = dispersed signal; i θ = dispersed part delay (regarding the position); δ = impulse function.
This model is compatible with biological tissues that have regular structures mixed with diffuse ones, such as liver tissue, and is similar to the ones found in the literature.

Evaluation
First, the methods were evaluated using 10,000 signals simulated with different levels of Ad (the ratio between the echoes' mean amplitudes from the diffuse and regular particles) and jitter (percentile variance of mean space between the regularly spaced particles).The signals were divided into 20 groups of 500, according to the characteristics in Table 1.The sampling frequency was 25 MHz, the transducer excitation frequency was 20 MHz, and the excitation signal band was 1.5 MHz.
The signals were simulated with MSS = 1.25 mm (compatible with the biological tissues), whose frequency related to periodicity (FreqMSS) is calculated according to the following equation: where, FreqMSS = frequency of maximum magnitude; c = ultrasound velocity; MSS = mean scatterer spacing.

Sponge phantom
After simulating signals we used a sponge phantom signals.The experiment was conducted by immersing the sponges in water.We used a transducer with a central frequency and bandwidth of 20 MHz to acquire the backscattered RF signals.The sponges were placed with the surface parallel to the scanning plan (Figure 3).The amplified RF signal was scanned with an 8-bit resolution oscilloscope (details on Table 2) (Pereira et al., 2002).
Figure 4 shows a section of each sponge signal with different spacing.Each of the 176 signals for each sponge was reduced to 512 points, a size sufficient to identify frequencies of the expected order of magnitude.

Signal treatment
To evaluate the three methods, all the signals were filtered by a Butterworth filter of order 6, with cutoff frequencies related to spacings of 0.02 mm (38.5 MHz) to 3.00 mm (256.67 kHz) for sponge signals, and 0.10 mm (7.70 MHz) to 5.00 mm (154.00 kHz) for the simulated signals.Cutoff frequencies were calculated for each signal part depending on the estimated US speed.For these values, the speed used was 1540 m/s (average speed for biological tissues).As the RF signal envelope has the same properties as the RF signal average spacing (Maciel, 2000), the analysis was performed only on the envelope so that it is more robust in estimating the periodic component.The signal envelope was calculated as the analytic signal obtained using an RF signal Hilbert transform, subtracting the average (to remove the DC level), and dividing by the variance for the energy normalization.

Results
For the three methods studied, an initial analysis using 10,000 simulated signals with various levels of Ad and jitter combined was performed.By using such a signal model, it was possible to control the signal characteristics independently.Sponges, which mimic the almost quasi-periodic spacing of tissues, were used to obtain signals closer to real biological tissues.
For the simulated signals, Figure 5 presents the results for the percentage of correct estimates for different levels of Ad.For each method, different curves are presented with a tolerance range of 10% around the value of the correct frequency in the simulated spacing -in this case, 616.00 kHz, in a medium with a US speed of 1.540 m/s.It was observed that for all methods the performance decreased with the increase of Ad or jitter levels, as expected.
The summarized results for the four phantoms for the three methods are presented in Table 3.The mean velocity for each sponge was obtained as the estimated mean velocity.Their values were 1.494 m/s sponge 0.1-0.2mm 1.486 m/s sponge 0.2-0.3mm, 1.499 m/s sponge 0.3-0.5 mm and 1.494 m/s for the sponge of 0.5 to 1.0 mm (Pereira et al., 2002).

Discussion
For the simulated signals the Tufts-Kumaresan method had the lowest performance, whereas SAC and MUSIC obtained similar results.It can be noted from Figure 5 that each method maintains a reasonable precision percentage even when the average amplitude of the aperiodic signal is about 40% of the average amplitude of the periodic signal (Ad = 41.2%) and jitter reached 10%.From this point, performance rapidly deteriorates indicating that the frequency is no longer evident.It is important to note that for the simulated signals, the dominant feature is the use of frequencies with greater or lesser amount of associated noise.In real cases, there are large isolated reflectors and/or highly attenuating structures that generate more complex signals.
For sponges, only the Tufts-Kumaresan method was able to identify an increase in the size of pores (frequency decreases as the size of pores increases), although it estimated higher frequencies than expected.
No outlier was excluded, although they clearly exist.If there was such an exclusion, the results would be closer to the expected values.This technique was used by some authors, such as Pereira et al. (2002).
It should be noted that the sponges do not have a preferred frequency, but a distribution of unknown frequencies.It is also observed that the sponges are composed of structures with various shapes and thicknesses, which cause a variation in the reflected energy that make up the backscattering signal under study.
The results with simulated signals show that SAC and MUSIC are potential methods to estimate the MSS.

Sponges Sampling rate
Pores between 0.1 and 0.2 mm 100 MHz Pores between 0.2 and 0.3 mm 100 MHz Pores between 0.3 and 0.5 mm 250 MHz Pores between 0.5 and 1.0 mm 250 MHz  Regarding sponge phantoms, SAC and MUSIC were able to estimate MSS values inside the periodicity range for one sponge (2 nd sponge), while Tufts-Kumaresan estimated values inside the periodicity range only once (1 st sponge).However, Tufts-Kumaresan was the only method able to estimate progressively higher MSS values; hence, it was able to sense the increase in the pore sizes of the sponges.A more detailed study is needed to investigate the sensibility of each method and their contour conditions.
transform of the envelope RF signal * ( ) Y f = Conjugate of the Fourier transform of the envelope RF signal

Figure 1 .
Figure 1.Example of application of methods, with simulated signal of center frequency of 3.5 MHz, sampling frequency of 25 MHz, and signal window of 512 points.(A) Tufts-Kumaresan (the estimated FreqMSS was 610.35 kHz); (B) SAC (the estimated FreqMSS was 610.35 kHz) and (C) MUSIC (the estimated FreqMSS was 611.55 kHz).

Figure 3 .
Figure 3. RF signal acquisition set-up using a commercial transducer (model M316; 3.2-mm diameter, Panametrics, Waltham, MA) with central frequency and bandwidth of 20 MHz.

Table 1 .
Characteristics of simulated signals.

Table 2 .
The sponge phantom with respective sampling rate.

Table 3 .
Results for four different sponge phantoms for the three methods for a US speed in the medium of 1.494 m/s.