Contrast Improvement in Sub- and Ultraharmonic Ultrasound Contrast Imaging by Combining Several Hammerstein Models

Sub- and ultraharmonic (SUH) ultrasound contrast imaging is an alternative modality to the second harmonic imaging, since, in specific conditions it could produce high quality echographic images. This modality enables the contrast enhancement of echographic images by using SUH present in the contrast agent response but absent from the nonperfused tissue. For a better access to the components generated by the ultrasound contrast agents, nonlinear techniques based on Hammerstein model are preferred. As the major limitation of Hammerstein model is its capacity of modeling harmonic components only, in this work we propose two methods allowing to model SUH. These new methods use several Hammerstein models to identify contrast agent signals having SUH components and to separate these components from harmonic components. The application of the proposed methods for modeling simulated contrast agent signals shows their efficiency in modeling these signals and in separating SUH components. The achieved gain with respect to the standard Hammerstein model was 26.8 dB and 22.8 dB for the two proposed methods, respectively.


Introduction
Introduction of contrast agents in ultrasound medical imaging has strongly improved the image contrast leading to a better medical diagnosis [1][2][3]. By adapting the transmitting ultrasound sequences composed of short wave trains to longer sinusoidal wave trains, it has been possible to enhance the harmonics detection witnessing of the presence of nonlinear explored media [3][4][5]. The most prominent example in echographic imaging is the second harmonic imaging (SHI) [3,6] which consists to send a sinusoidal wave train of frequency 0 and to receive the backscattered signal at twice the transmitted frequency, that is, 2 0 (see Figure 1).
Although the second harmonic imaging possesses undoubted advantages compared to standard echographic imaging, contrast harmonic imaging, however, has image contrast limitations related to the presence of harmonic components of nonlinear nonperfused tissues [7]. This contrast reduction can be overcome by proposing no more contrast harmonic imaging but rather contrast subultraharmonic (SUH) imaging [8,9]. Under certain conditions of incident frequency and pressure levels, this solution has been envisaged as a serious alternative [10,11] since it has been shown that only contrast agent is capable of supplying SUH components sufficient to construct perfused tissue images with a strong contrast. Contrast SUH imaging consists to send a sinusoidal wave train of frequency 0 and to extract from the backscattered signal only SUH frequencies at 0 /2, (3/2) 0 , (5/2) 0 , . . . (see Figure 2).
To extract such SUH components from the whole spectrum, a certain number of approaches called "black box methods" has been proposed such as those based on the multiple input and single output (MISO) Volterra filtering [12][13][14]. These recent methods are capable of accurately modeling the signal backscattered by the contrast agent with adequate values of order and memory of such models. However, these methods are quite complex methods, and they do not give an extraction of harmonic components 0 , 2 0 , 3 0 , . . . and SUH components ( 0 /2)2, (3/2) 0 , (5/2) 0 , . . . since they model all spectral components ( 0 /2)2, 0 , (3/2) 0 , 2 0 , (5/2) 0 , 3 0 , . . ..  In order to reduce the complexity of such methods and to extract SUH components from all spectral components, we propose two new original approaches neither based on the Volterra filtering but based on the Hammerstein filtering. This paper is organized as follows: after recalling standard Hammerstein model, the new methods are presented. To validate our methods, we propose realistic simulations of contrast agents signals. Then a quantitative and qualitative comparison is made between the two proposed methods with respect to standard Hammerstein model. Finally a discussion completed by a conclusion closed the paper.

Methods
Polynomial Hammerstein model is a special type of nonlinear filters in which a static nonlinear system is followed by a dynamic linear system [15]. From this model, the nonlinear system is approximated by a polynomial function, and the linear part is a finite impulse response (FIR) filter. The block diagram of Hammerstein model is shown in Figure 3.
As for the Volterra decomposition, the Hammerstein decomposition is able to model harmonic components, but it was unable to model sub-and ultraharmonic (SUH) components. Before explaining how it was possible to model SUH components, we recall the Hammerstein decomposition.

Hammerstein Decomposition.
The Hammerstein modeled signal̂( ) can be seen as the summation of signals ( ) coming from parallel branches: In each branch, the signal ( ) is the output of a linear filter ℎ ( ) of input ( ): where is the output signal of the non linear system and E is the symbol of the mathematical expectation. The corresponding solution is if (W W) is invertible. Otherwise, regularization techniques can be used. Consequently, identifying a nonlinear system of input ( ) and output ( ) with an Hammerstein model is equivalent to calculate the signal̂( ) using (4) and (5).

Sub-and Ultraharmonics
Modeling. As previously mentioned, Hammerstein model is able to model harmonic components only. This is justified by the steady state theorem reported in [16] which stipulates that the output response of the model to a periodic input of frequency 0 is also a periodic signal of fundamental frequency 0 . Suppose that the transmitted signal ( ) is periodic of frequency 0 , and that the signal backscattered by the contrast agent has the following spectral components ( 0 /2)2, 0 , (3/2) 0 , 2 0 , (5/2) 0 , 3 0 , . . .. Let̂1 1 ( ) be the output signal of the Hammerstein model of input ( ) (see Figure 4). According to the theorem reported in [16], the spectral content of the output signal̂1 1 ( ) is composed of the following harmonic components 0 , 2 0 , 3 0 , . . . if and only if the input frequency is 0 .
Thus, modeling SUH components of frequency 0 /2, (3/2) 0 , (5/2) 0 , . . . using Hammerstein model is possible if and only if these components could be seen as integer multiples of the input frequency to the Hammerstein model. To do this, some modifications must be performed. Two types of modifications are proposed in this work, either at the input or at the output. Based on these two types of modifications, two methods for modeling and separating the sub-and ultraharmonic components using Hammerstein models are described in this section.
Each of the two proposed methods consists of two steps: one step for harmonic modeling and another step for SUH modeling. The first method is based on the modification of the input frequency, while the second one is based on the modification of the output frequency.

Method 1: Modeling by Input Frequency Shifting.
As previously mentioned, this method consists of two steps; each step uses one Hammerstein model as presented in Figure 4.
(1) Harmonic modeling: harmonic modeling is done by identifying the system of input 11 ( ) = ( ) and output ( ) with an Hammerstein model. The obtained signal̂1 1 ( ) has the harmonic components only. Referring to (4), the harmonic signal̂1 1 ( ) could be written aŝ with 11 and 11 defined as in Section 1.
(2) sub-and ultraharmonic modeling: the SUH information is found in the difference signal 21 ( ) between the output signal ( ) and the harmonic signal̂1 1 ( ): The spectral content of 21 ( ) is composed of 0 /2, (3/2) 0 , (5/2) 0 , . . . by referring to the previous theorem; these components could be modeled using Hammerstein model if the input frequency is 0 /2. Consequently, the initial input signal ( ) must be modified in such a way to bring up the subharmonic frequency 0 /2. To do this, the spectrum of ( ) is downshifted by 0 /2 to shift the frequency 0 toward the position of 0 /2. The modified input 12 ( ) is calculated according to the following equation: where R is the real part,̃( ) = H( ( )) is the Hilbert transform of ( ), and is the sampling frequency. The modified input signal 12 ( ) could be written in the following vectorial form: Then, the SUH signal̂1 2 is the output of a Hammerstein model that identifies the system of input 12 ( ) and output 12 ( ). In the same way,̂1 2 ( ) is calculated according to the following equation:̂1 Finally, the total modeled signal̂( ) is the sum of the harmonic signal̂1 1 ( ) and the sub-and ultraharmonic signal For this method, note that the maximal frequency that could be modeled is limited by the order of the Hammerstein model. It is well known that the Hammerstein model of order excited by a signal of frequency 0 can model harmonic components until the th harmonic of

Method 2: Modeling by Output Frequency
Shifting. This method also consists of two steps, on step dedicated for harmonic modeling and another step dedicated for SUH modeling. Each step uses one Hammerstein model as presented in Figure 5.
(1) Harmonic modeling: this step is the same as the first step of the method 1. The signal ( ) is modeled with a Hammerstein model of input 21 ( ) = ( ). The obtained signal is the harmonic signal̂2 1 ( ) calculated according tô (2) Sub-and Ultraharmonic modeling: based on the same idea as reported in method 1, SUH components could be modeled when they are considered as integer multiples of the input frequency. In this method, we propose to keep the input signal ( ) and to change the output signal ( ) by upshifting its spectrum of 0 /2. Now, the SUH components are shifted toward the harmonics position, and then they could be modeled with a Hammerstein model. The modified output signal mod ( ) is calculated according to the following equation: In vector form, Then the signal mod ( ) is modeled with a Hammerstein model of input ( ). The obtained signal̂2 2 mod ( ) is calculated according to the following equation: has SUH components upshifted by 0 /2. To recover the sub-and ultraharmonic signal̂2 2 ( ), the signal̂2 2 mod ( ) is downshifted of 0 /2 according to the following equation: The final signal̂( ) is the sum of the harmonic signal̂1 1 ( ) and the sub-and ultraharmonic signal̂2 2 ( ):

Simulations and Results
To validate the two proposed methods and to quantify their performance in ultrasound medical imaging, realistic simulations are proposed. To achieve these simulations, the free simulation program Bubblesim developed by Hoff [17] was used to calculate the scattered echoes for a specified    [18] and then modified by Hoff [17] is based on the theoretical description of contrast agents as air-filled particles with surface layers of elastic solids. In order to simulate the mean behavior of a contrast agent cloud, we hypothesized that the response of a cloud of contrast agents was times the response of a single contrast agent with the mean properties. The incident burst is a sinusoidal wave composed of 18 cycles of frequency 0 = 4.5 MHz and pressure of 0.6 MPa [13]. (The resonance frequency of the encapsulated contrast agent of 1.5 m is about 2.5 MHz. The emission frequency at 4.5 MHz is nearly the double of the resonance frequency.) Under the previous conditions of frequency and pressure, the oscillation of the contrast agent is nonlinear with sub-and ultraharmonic generation. The sampling frequency is = 1/ = 60 MHz. The parameters of the contrast agent are given in the Table 1.
In this research work, the performances of the different methods are evaluated qualitatively and quantitatively. Figure 6 represents a qualitative comparison in both time and frequency domains between the signal backscattered by the contrast agent ( ) and the signal obtained with the standard Hammerstein model, method 1, and method 2. Method 1 is applied with a Hammerstein model of order = 3 and memory = 30 for the first step and a Hammerstein model of order = 5 and memory = 30 for the second step. Method 2 is applied with a Hammerstein model of order = 3 for the two steps and memory = 30. Figure 6(a) (top) shows that the modeled signal with the standard Hammerstein model does not describe correctly the signal backscattered by the contrast agent. Corresponding spectra in Figure 6(b) (top) show that the signal modeled with the standard Hammerstein model has the harmonic components only ( 0 , 2 0 , and 3 0 ). Figure 6(b) (middle, bottom) shows that the signals modeled with methods 1 and 2 perfectly describe the contrast agent signal. Corresponding spectra on Figure 6(b) (middle, bottom) show that all the frequency components are modeled: harmonics at ( 0 , 2 0 , 3 0 ), subharmonic 0 /2, and ultraharmonics ((3/2) 0 , (5/2) 0 ). Figure 7(a) (top) shows the different signal between ( ) and the harmonic signal (in black) and the sub-and ultraharmonic signal obtained with method 1 (top) and method 2 (bottom). We can see the good agreement between the signals. Spectra on Figure 7(b). (top) confirms that subharmonic 0 /2, first ultra-harmonic (3/2) 0 , and second ultraharmonic (5/2) 0 are well modeled and separated from other harmonic components. These results confirm the efficiency of the two proposed methods in modeling and separating the sub-and ultraharmonics present in contrast agent response.

Quantitative Evaluation.
To quantify the performance of each method and to know which method provides the best performance, a quantitative study is necessary. The relative mean square error EQMR defined as  is evaluated for different noise levels at the system output. The noise level adjusted as the function of SNR (signal to noise ratio) is Gaussian and white. Ten realizations are made to evaluate the fluctuations of RMSE. RMSE for SNR = ∞, 20, 15, 10 dB is reported in Figure 8. These simulations show that the RMSE achieved with the two proposed methods 1 and 2 is always less than the RMSE achieved with the standard Hammerstein model for the different SNR values.
The gap between the standard model and the two methods 1 and 2 decreases when the SNR value increases. A gap ranging from 4 to 26 dB could be obtained depending on the SNR conditions. These results confirm that the standard Hammerstein model is not adapted for sub-and ultraharmonic modeling. A zoom on Figure 8(d) shows that the RMSE varies slightly around a mean value. This result shows that the two methods 1 and 2 are robust toward noise. Note that the curves of variation of RMSE obtained with the two methods have the same trend, indicating that the two methods tend toward the optimal solution.

Discussions and Conclusions
In this research work the problem of modeling sub and ultraharmonics with Hammerstein model is presented. Usually, the standard Hammerstein model is able to model harmonics only, which are integer multiples of the input frequency. Sub and ultra-harmonics could not be modeled.
In this work, we propose for the first time two new methods that use Hammerstein filters that model sub and ultra-harmonics. The two methods are based on the same idea stipulating that modeling SUH with Hammerstein model is possible if the input signal or the output signal is modified In such a way that the SUH components become in the position of integer multiples of the input frequency. Each method uses two Hammerstein models successively. The first one is dedicated to model harmonic components and the second one to model the SUH components.
The first method (method 1) applies a spectral downshifting of 0 /2 on the input signal of the Hammerstein model. Now, SUH are seen as integer multiples of the input frequency, and therefore they can be modeled with Hammerstein model. In this step, the order of Hammerstein model is an important parameter that need to be adjusted to ensure the modeling of all SUH. For the second method (method 2), a spectral upshifting of 0 /2 is applied on the output signal to move the SUH components toward the harmonic positions. Then, SUH components can be modeled with Hammerstein model excited with the input signal of frequency 0 . Finally, a last spectral downshifting is performed to recover the SUH signal.
The two proposed methods are characterized by its simplicity. The originality of these methods is that they allow for the first time both the modeling of contrast agents signals and the separation of SUH components of the contrast agent response.
However, the two methods do not present the same advantages and disadvantages. Method 1, which is based on the modification of the input signal, is less sensitive to the noise compared to method 2, which is based on the modification of the output signal. This is due to the fact that all the noise generated in the different parts of the non linear procedure are added to the output.
On the other hand, although method 1 has a more simple structure, it is slower than method 2. This is due to the fact that the second step of method 1 requires an order higher than method 2, the order of the first step being fixed. And as the computation time is related to the order of the model, higher the order, the slower the method.
The application of the proposed methods for modeling the contrast agents response shows their efficiency in modeling and in separating SUH components from other harmonic components. Gains of 25.8 dB and 22.8 dB in term of the RMSE are achieved with methods 1 and 2, respectively, compared to the standard Hammerstein model. The achieved error (RMSE) gain by the two methods is related to the sub and ultraharmonics energy initially presented in the output signal of the non linear system. The more important the energy of the sub-and ultraharmonics, the more important the gain. Although in this paper, the two methods works well 8 International Journal of Biomedical Imaging for 0 /2, (3/2) 0 , (5/2) 0 , these methods can be extended for other orders. The two proposed methods find theirs applications in the field of sub and ultra-harmonic contrast imaging in order to produce high contrast images. This work opens a new research axis for new modeling techniques of SUH using Hammerstein model or any other non linear models.