Introduction

Ultrasonic methods are widely used to investigate the internal structures of several kinds of media, particularly in biological subjects. This is mainly due to ultrasound (US) being a nonionizing radiation able to generate real-time images at low cost. Its most attractive application is to quantitatively describe the US properties of tissues with the objective of differentiating normal and pathological ones.

The mean scattering space (MSS) is also considered to be an important tool for tissue characterization in addition to classical parameters that can be estimated from US data of biological tissues, such as scattering cross-section, wave speed, and attenuation coefficient ( ^{Bridal et al., 1998} ). The MSS is related to the natural periodicity of structures that exist in some organs, such as the liver and kidney ( ^{Simon et al., 1997} ).

There are several techniques available in the literature for estimating MSS. For example, ^{Wear et al. (1993)} compared methods using the Burg algorithm and Fast Fourier Transform (FFT). ^{Varghese and Donohue (1993} , ^{1994} , ^{1995} ) proposed a spectral correlation function to calculate the mean scattering space. ^{Simon et al. (1997)} presented an algorithm that uses the magnitude information of the spectrum of radio frequency (RF) signals.

The above mentioned techniques are applied to RF signals, or to their envelopes, which contains a mixture of contributions from periodic and non-periodic scatterers. To separate periodic and non-periodic signals, ^{Pereira and Maciel (2001)} used singular spectrum analysis (SSA) to estimate the MSS. Other studies considering SSA used phantoms ( ^{Pereira et al., 2002} ) and bone *in vitro* ( ^{Pereira et al., 2004} ).

^{Kauati et al. (2012)} evaluated the Burg, Wiener, and MUSIC spectral estimation methods for automatic MSS calculation without requiring operator intervention. The evaluation of these methods was carried out using 10,000 simulated signals generated by a model where the variables of interest can be controlled. Signals from a periodic phantom using nylon wires were also used. The subspace method MUSIC provided the best performance.

Based on RF ultrasound signals, ^{Donohue et al. (2001)} used the generalized spectrum (GS) to analyze scatterer properties and compared them with the conventional texture analysis (CTA) for classifying breast tumors. ^{Huang et al. (2000)} used methods based on the GS and cepstrum to detect and estimate duct wall spacings in signals from ultrasonic phantom experiments and Monte Carlo simulations. ^{Chen et al. (2011)} proposed a method combining the Simon’s method ( ^{Simon et al., 1997} ) and the temporal autocorrelation to avoid spurious spectral peaks related to highly reflecting targets. Their proposed method outperformed Simon’s method alone.

Since then, the publication rate of studies related to periodicity estimation has slowed down. Recently, ^{Pan et al. (2014} , ^{2015} ) presented a method of mode decomposition-based cepstrum spacing measurement of MSS changes in quasi-periodic trabecular bones. In 2013, ^{Rubert and Varghese (2013)} compared multi-taper coherence and single-taper to estimate MSS and concluded that the multi-taper significantly improves over single-taper. Periodicity estimation remains an open problem as can be seen in the recent paper from ^{Zhou et al. (2017)} , that makes a detailed comparison of the 13 most common methods used for MSS estimation, regarding their advantages and limitations. They proposed that a comparative analysis of state-of-the-art MSS estimation methods is, thus, needed in the future, considering their estimation accuracy, time efficiency, and robustness.

This paper introduces a method to estimate the MSS by using SSA in association with the concept of entropy. SSA is a means of optimizing the separation between noise and signal subspaces and entropy is used as a quantitative criterion to establish a boundary between the two. The results of Monte Carlo simulations and real signals from a phantom are presented. These results were compared with the SSA method used by ^{Pereira and Maciel (2001)} and the SIMON method proposed by ^{Simon et al. (1997)} . The SIMON method was chosen because it performed better than several other methods that we tested, such as Burg ( ^{Kauati et al., 2012} ), MUSIC ( ^{Kauati et al., 2012} ), SAC ( ^{Varghese and Donohue, 1993} , ^{1994} , ^{1995} ), ^{Tufts and Kumaresan (1982)} , and Wiener ( ^{Kauati et al., 2012} ).

Methods

The analysis was carried out on simulated US signals backscattered from biological soft tissues and subsequently on phantom signals. In this section, we present the main concepts regarding the tissue model, SSA, and Entropy with the steps followed for their implementation. We also describe the real phantom signals.

Model of the simulated tissue

A unidimensional model ( ^{Pereira and Maciel, 2001} ) was utilized to simulate backscattered RF signals from biological soft tissues. This model is based on the hypothesis that the medium possesses linear properties, so the signal (echo) received r(t) can be characterized as:

where: *p (t)* = transmitted ultrasound pulse;

*g(t)* = function characterizing the medium (impulse response);

*n(t)* = experimental system noise;

*** = convolution.

The function *g(t)* is modeled as a sum of regularly and randomly spaced particles in the medium:

where: *N* = total number of regularly spaced particles;

*M* = total number of irregularly spaced particles;

*C _{i}* = amplitude of the

*i*

^{th}regularly spaced particle;

τ_{i} = time-delay of the *i*^{th} regularly spaced particle (regarding the position);

*b _{i}* = amplitude of the

*i*

^{th}irregularly spaced particle;

*θ _{i}* = Time-delay of the

*i*

^{th}irregularly spaced particle (regarding the position);

*δ* = impulse function.

The simulation algorithm of the RF signals is described below:

1. Definition of the transmitted signal, medium, and US pulse.

1. − The emitted US is defined by the sampling rate, transducer exciting frequency, and excitation signal band;

1. − The medium is characterized by the velocity of the US in it, mean distance between the regular particles and their percentile variance (jitter) around their exact periodic position, maximum and minimum amplitudes of the regular particles, ratio between the mean amplitudes of the echoes from the irregular and regular particles (Ad), signal-to-noise ratio between the regular signal and white noise in dB, and number of regular and diffuse particles;

1. − The US pulse has a Gaussian-shaped envelope.

2. Definition of the scattering characteristics.

2. − The unidimensional model is equivalent to aligning the scattering of regular and irregular particles;

2. − The position of the irregular scatterers follows a uniform distribution and the position of the regular scatterers follows a gamma distribution;

2. − The amplitudes of the regular and irregular scatterers are random with a uniform distribution limited by the parameter Ad for the irregular particles and by the maximum and minimum amplitudes for the regular particles.

3. Building the RF signal by the convolution of the pulse and medium particles along with the noise ( Equation 1 ).

Figure 1 shows an example of a simulated signal with an MSS of 1.25 mm, Ad of 11.7% and sampling frequency of 25 MHz.

The standard periodicity of the time signal can be estimated from the first peak of the PSD (power spectral density). The analysis uses the envelope signals only. This envelope can be obtained as the modulus of the Hilbert transform of the RF signal with its mean subtracted from it and divided by its variance for normalization purposes.

In this example, periodicity is evident by a simple visual inspection, for didactic purposes. In fact, we simulated 20 combinations of the jitter and Ad parameters, obtaining much more complex signals where the periodicity could not be identified visually. See Table 1 for the details.

Group |
Jitter (%) |
Ad (%) |
Group |
Jitter (%) |
Ad (%) |
Group |
Jitter (%) |
Ad (%) |
---|---|---|---|---|---|---|---|---|

1 | 1 | 11.7 | 9 | 10 | 11.7 | 17 | 20 | 11.7 |

2 | 1 | 41.2 | 10 | 10 | 41.2 | 18 | 20 | 41.2 |

3 | 1 | 58.8 | 11 | 10 | 58.8 | 19 | 20 | 58.8 |

4 | 1 | 76.5 | 12 | 10 | 76.5 | 20 | 20 | 76.5 |

5 | 5 | 11.7 | 13 | 15 | 11.7 | 21 | 30 | 11.7 |

6 | 5 | 41.2 | 14 | 15 | 41.2 | 22 | 30 | 41.2 |

7 | 5 | 58.8 | 15 | 15 | 58.8 | 23 | 30 | 58.8 |

8 | 5 | 76.5 | 16 | 15 | 76.5 | 24 | 30 | 76.5 |

Singular spectrum analysis

SSA is based on the well-known principal component analysis (PCA) which projects an *N* data vector , with respect to an optimized orthonormal basis (*E ^{k}, 1 ≤ K ≤ M*).

*E*are the eigenvectors obtained from diagonalization of the correlation matrix of , which is constructed by moving a window of

^{k}*M*points along

*x*. The

*M*data points inside the window can be described by a linear combination of the vectors

*E*:

^{k}where: 1 ≤ j ≤ M;

M < N;

1 ≤ i+j ≤ N;

1 ≤ i ≤ N-M+1.

The projection coefficients
*E ^{k}
* are the orthonormal functions. The vectors

*E*are the eigenvectors of the cross-correlation matrix

^{k}*C*of the sequence

_{x}*x*.

According to ^{Vautard and Ghil (1989)} , if a periodic component is present in a time series, its periodicity can be found from the power spectrum of the signal reconstructed with the eigenvectors corresponding to the periodic component. Thus, SSA can be used to separate the periodic and non-periodic structures of a medium.

As an example, SSA was applied to the envelope signal in Figure 2 A (cross-correlation matrix with dimensions of 200 × 200). Figures 2 2B show the signal reconstructed using the first four eigenvectors (i.e., first two pairs) and its power spectrum, respectively. A peak was found at 610 kHz, which corresponds to an estimated MSS of 1.26 mm (the original simulated MSS was 1.25 mm). Figures 2 2D show the residual signal obtained by the subtraction of the reconstructed signal from the original one and its power spectrum, respectively.

Figure 3 shows a set of the first 50 eigenvalues of the 200 × 200 correlation matrix arranged in descending order according to their magnitude. The arrows show the eigenvalue pairs responsible for the periodicity of the signal. In fact, the first pair corresponds to the peak in the power spectrum of Figure 3 , and the second pair corresponds to the 2^{nd} peak in the same power spectrum.

The next step is to seek a quantitative criterion that can be used to define how many eigenvalue pairs should be considered to reconstruct the signals to include all the periodic information and minimize the noise.

If the periodic part of the signal is considered as the “organized signal” compared to the aperiodic and noisy part, the concept of entropy–associated with the level of order of the medium–may help with the separation.

Thus, as periodicity is always related to eigenvalue pairs ( ^{Vautard and Ghil, 1989} ), it can be recovered by inspecting the corresponding eigenvectors. The question that remains is how many eigenvectors should be used to reconstruct the periodic part of the signal. This is treated in the next section.

Entropy

The concept of entropy is applied to determine how many eigenvectors should be used in the linear combination for signal reconstruction. Entropy is a measure of the uncertainty of a random variable, i.e., how much information, on average, is necessary to describe a random variable. One of its basic definitions includes a function that is not dependent on values of the random variable *x*, but on its probability distribution, as shown in Equation 4 ( ^{Cover and Thomas, 1991} ).

where: *H(x)* = absolute entropy function;

*x* = discrete random variable;

*p(x) =* probability density function.

Another useful concept is the relative entropy, which is the measure of the distance between two distributions:

where:

*p(x)* = probability distribution;

*q(x) =* probability distribution.

In this work *q(x)* was calculated from the envelope of the RF signal and *p(x)* from the signal reconstructed using the PCs.

To test the concept of relating entropy to periodicity, simple simulated signals composed of 10 sinusoids (0.25, 0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00, 2.25, and 2.50 MHz) with added noise (SNR: 7.0, 3.5, and 2.0 dB) were generated. SSA was applied to the signals and the respective eigenvalue/eigenvector sets were obtained from the 200 × 200 correlation matrix.

The 10 sinusoidal frequencies correspond to 10 eigenvalue/eigenvector pairs. The next step was to reconstruct each signal by sequentially adding the two eigenvectors corresponding to the eigenvalues in descending order of magnitude. Each subsequent reconstruction was performed by adding the two eigenvectors corresponding to the next two largest eigenvalues. The process was repeated until the last eigenvalue pair was used. The relative entropy was calculated for each residual signal (RD) that remained after the reconstructed signal was subtracted from the original signal as each new eigenvector pair was added. A curve of RD as function of the number of eigenpairs was built. The minimum relative entropy value was obtained for the residual signal RD when the reconstructed signal had 10 sinusoids, which was probably due to the original signal being composed of exactly ten sinusoids. The RD curve seemed to be an indication that the position of the minimum value may be used as a numerical criterion to estimate the number of periodic signals present in the analysis. Thus, the next step was to verify if this relative entropy behavior is also exhibited by the simulated RF signals.

Estimating the MSS based on entropy: Monte Carlo Simulation

The main objective of this work was to estimate the MSS of a regular structure from US signals. Here, we propose to use SSA and apply entropy as a criterion to define the number of eigenvector pairs used to reconstruct the periodic part of the signal. The basic steps to estimate the MSS are as follows:

1. Apply SSA to the US scattered signal whose envelope is obtained by the magnitude of the Hilbert transformation without the DC component and divided by its variance;

2. The cutoff frequencies of the 6th-order bandpass filter were chosen to limit the MSS window between 0.05–5.0 mm, which are enough for biological tissues. The filter order was determined empirically. The sixth-order filter has a flat band around the central frequency so the filtered signal has basically no distorting and tests with higher order filters did not improve the results;

3. Obtain the set of eigenvalues and eigenvectors of the envelope of the RF signal;

4. Reconstruct the signal by sequentially adding pairs of eigenvectors according to the order of the magnitude of the corresponding eigenvalues;

5. Obtain the relative entropy estimation for the residual signal as each pair is added to the reconstructed signal;

6. Find the minimum on the relative entropy curve. The number of pairs that correspond to this minimum is used to make the final reconstruction of the signal;

7. Reconstruct the signal with the number of eigenvectors pointed out by the entropy curve, obtain the FFT and estimate the frequency with the highest amplitude;

8. Use this frequency to estimate the MSS.

where: ν is the speed of ultrasound in the propagation medium (1540 m/s for water).

These eight steps were applied to 24,000 simulated US signals with different levels of Ad and jitter ( Table 1 ).

Phantom

One wire phantom was used with four 0.5 mm diameter nylon wires aligned in a vertical plane with a regular spacing of 1.2 mm. The phantom was immersed in a water bath (no diffuse scatterer) and the backscattered RF signals were acquired with a 20 MHz center-frequency transducer. For data acquisition, the transducer was placed parallel to the length of the wires on a B-scan plane at an angle of 10–15° to the vertical plane of the wires to assure that all wires were insonified. The acquired RF signals were amplified and digitized at a 100 MHz sampling rate with an 8-bit oscilloscope. Eight bits was a limitation of the equipment but the results presented small MSS error acceptable for the application. We collected 152 signals from each phantom. The US speed in the water was assumed to be 1,498 m/s. A detailed explanation of the experimental setup was given by ^{Pereira and Maciel (2001)} .

Results

One thousand US signals were simulated for each of the 24 groups ( Table 1 ) and 24,000 MSS estimates were obtained. A “correct” estimate was defined within a 5% window around the real value (1.25 mm). Figure 4 presents the results as a function of the jitter and Ad parameters for: (A) the SSA method and entropy, (B) the SSA method with two pairs of eigenvalues, and (C) SIMON method. SSA with two pairs of eigenvalues was used for the comparison because of its performance according ^{Pereira and Maciel (2001)} .

The wire phantom was used to record 152 signals. The SSA with entropy method identified the most frequent MSS to be 1.096 mm for the wire phantom spaced at 1.2 mm ( Figure 5 A). When two pairs of eigenvalues were used, the MSS characterized as 0.959 mm for 13 signals, while for the other 139 signals the MSS was estimated to be equal to 1.096 mm ( Figure 5 B). The SIMON method estimated the MSS as smaller than 1.096 mm for 30 and equal to 1.096 mm for 122 signals ( Figure 5 C).

Discussion

A method was presented to estimate the MSS of US backscattered signals from biological tissues based on SSA combined with the concept of entropy. The relative entropy was used as a means to establish a numerical criterion to define the number of eigenvalues for the reconstruction of the signal subspace.

The reconstructed signal was used to estimate the MSS of periodic structures for simulated signals with the Monte Carlo approach giving consistent results even for high levels of noise and jitter.

Regarding the real phantom signals, the results of the SSA + Entropy method were compared to the SSA implemented by ^{Pereira and Maciel (2001)} and the method of ^{Simon et al. (1997)} , as these two achieved superior periodicity estimation performance compared to others found in the literature. The smaller MSS value obtained by every method can be explained by the fact that the B-scan plane enclosed an angle of 10–15° with the vertical plane of the wires, so the MSS seen being in fact smaller.

In terms of MSS estimation, the results of the ^{Pereira and Maciel (2001)} and ^{Simon et al. (1997)} were similar to those of our method for the same 24,000 simulated US signals ( Figure 4 ). The SIMON method was chosen because it performed better than several other methods tested using the same data ( ^{Kauati et al., 2012} ). Our method performed slightly better for the wire-phantom signals ( Figure 5 ). The method of ^{Simon et al. (1997)} does not keep phase information, and hence, it is not capable of the characterization of the irregular signals.

Comparing our results with other works in the literature is not straightforward, as the signals used to test the MSS estimation differed. ^{Zhou et al. (2017)} concluded that a comparative analysis of MSS estimators with respect to estimation accuracy, time efficiency and robustness should be carried out.

In our case, one comparable work is the technique from ^{Chen et al. (2011)} as they proposed a method based on the SIMON method ( ^{Simon et al., 1997} ) in combination with temporal signal autocorrelation. They have simulated signals similar to ours and–for a jitter of 1%, 5%, 10% and 15%, and Ad of 10% ( Figure 4 , p. 605 of their paper)–it is possible to see that the coefficient of variation (standard deviation/mean) is higher than 10%. Meanwhile, in our case ( Table 2 ), the coefficient of variation is equal to or less than 6.0%. Therefore, our proposed method outperforms their method.

This work proposed the combination of SSA with entropy to estimate the MSS of a periodic or quasi-periodic medium that resulted in estimation capabilities similar or better compared to two other methods found in the literature. The novelty of our method is the application of entropy as a quantitative criterion for selecting the SSA periodic components, achieving independency of heuristic criteria.