Nội dung

Hindawi Publishing Corporation EURASIP Journal on Advances in Signal Processing Volume 2011, Article ID 146175, 10 pages doi:10.1155/2011/146175 Research Article Bayesian Spectral Estimation Applied to Echo Signals from Nonlinear Ultrasound Scatterers Yan Yan,1 James R. Hopgood (EURASIP Member),1 and Vassilis Sboros2 1 Institute for Digital Communications, School of Engineering, University of Edinburgh, EH9 3JL Edinburgh, UK of Medical Physics, University of Edinburgh, EH16 4TJ Edinburgh, UK 2 Department Correspondence should be addressed to James R. Hopgood, james.hopgood@ed.ac.uk Received 15 July 2010; Accepted 27 October 2010 Academic Editor: Elias Aboutanios Copyright © 2011 Yan Yan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The understanding and exploitation of acoustic echo signals from nonlinear ultrasound scatterers is an active research area that aims to improve the sensitivity and specificity of diagnostic imaging. Discriminating between acoustic echoes from linear scatterers, such as tissue, and nonlinear scatterers, such as contrast microbubbles, based on their frequency content is also an important topic in ultrasound contrast imaging. In order to achieve these objectives, a fundamental preliminary stage is to extract information about the reflected signals in the frequency domain with high accuracy: this is essentially a feature extraction and estimation problem. In this paper, a parametric Bayesian spectral estimation method is utilised for the analysis of the backscattered echo signals from microbubbles. In contrast to existing nonparametric discrete-Fourier-transform- (DFT-) based spectral estimation techniques used in the ultrasonic literature, this method is able to estimate the number of spectral components as well as their amplitudes and frequencies. The Bayesian spectral analysis technique has improved frequency resolution compared with the DFT for short multiple-component signals at low signal-to-noise ratios. The performance of the method is demonstrated with simulated signals, as well as analysing experimentally measured echo signals from nonlinear microbubble scatterers. 1. Introduction Ultrasound contrast agents (UCAs) used for enhancing ultrasound images were first discovered accidentally by cardiologist Charles Joiner in the 1960’s [1]. In particular, microbubbles (MBs) have been widely used as UCAs in biomedical research since the 1990’s [2]. Whereas imaging techniques such as positron emission tomography (PET), magnetic resonance imaging (MRI), and computed tomography (CT), have disadvantages due to toxicity levels, leakage from the vascular system, and a high cost of repetitive imaging, ultrasound contrast imaging (UCI) provides a nondestructive, noninvasive, and cost-effective paradigm. In general, UCAs greatly improve resolution and sensitivity of ultrasound imaging. Nevertheless, a major limitation remains since imaging results are highly dependent on the scanner operator and transmit (insonification) waveform selection. Though some adaptive waveform technologies have recently been introduced into commercial systems, this imaging modality would benefit greatly from more sophisticated adaptive waveform technologies such as those used in radar imaging which, in turn, are being biologically inspired by the waveform selection used in nature by, for example, bats and dolphins. In order to analyse the performance of ultrasound contrast imaging systems and the ability to detect the presence of MBs in blood and being able to distinguish them from tissue, it is important to analyse the acoustic responses of individual components, for example, a single MB scatterer [3–5]. This paper analyses the characteristic backscattered waveforms received from a single MB as a result of waveform insonification. The analysis in this paper focuses on the spectral content of these waveforms, while future work will focus on the spectral-temporal content. The Fourier transform is ubiquitous in the spectral analysis of ultrasound contrast imaging and, though very useful, has inherent limitations. In particular, as a direct result of the time-frequency uncertainty principle, it is known that traditional nonparametric methods suffer from the trade-off between pulse duration and the width of 2 EURASIP Journal on Advances in Signal Processing the spectral peak. Since Fourier-based theory makes very weak assumptions regarding the structure of a particular signal, the Fourier transform (FT) may provide false spectral content and cannot provide time domain information: both are important in MB detection. In addition to these fundamental limitations, FT-based signal detection often relies on, for example, ad hoc signal-dependent criteria such as thresholds or other heuristics for peak detection. An example of a situation where the FT alone cannot provide all the required information is when determining the number of dominant frequency components present in a signal; this number has been shown to characterise the type of scatterer [2, 6–8]. Thus, distinguishing tissue and MB responses in UCI using spectral and temporal features based on FT analysis is limited. Fourier methods can be improved upon by explicitly incorporating strong prior information about the structure of the signals and systems involved. Bayesian spectral analysis [9–11] leads to frequency estimators in which the achievable spectral resolution is directly dependent on the signal-tonoise ratio (SNR) and can be orders of magnitude better than that of a conventional Fourier power spectrum or periodogram [12, 13]. A key advantage of Bayesian spectral estimators is its ability to produce a measure of uncertainty of the estimated parameter, either through the full probability density function (PDF) or via characteristic features such as an estimate of variance of the frequencies. The pulse-echo returns from discrete MB events can be analysed by explicitly modelling the signal using parametric representations. Although the ultrasound environment and the MB response exhibit nonlinearities, the returned signal is still a finite-duration pulse typically with a welldetermined structure. In this paper, a Bayesian spectraltemporal estimator is proposed in which the echo responses are explicitly modelled as a linear combination of an unknown number of sinusoids with unknown frequencies, amplitudes, and phases. Using Bayesian inference, a posterior PDF, p(k, Φk | y), for the pulse parameters, Φk , can be derived, where y is the observed data and k is the number of sinusoids in the signal pulse. Although this PDF encapsulates all the information needed to characterise the ultrasound echo signal, it is still necessary to find a point estimate for the “optimal” parameter set. This could, for example, be a marginal maximum a posteriori (MMAP) or minimum mean-square error (MMSE) estimate which reduces to either an optimisation or integration problem, respectively:   Φ k k,    Φ k k, MMAP  MMSE   = arg max p k, Φk | y , k,Φk   (1) = E p(k, Φk |y) k, Φk | y , where E p(x| y) [x | y] denotes the conditional expectation with respect to the PDF p(x | y). Unfortunately, both optimisation and integration of the PDF in (1) are extremely difficult numerical problems. However, Markov chain Monte Carlo (MCMC) and other numerical Bayesian methods are simulation-based techniques that provide enormous scope for realistic statistical modelling [14] and allow for parameter estimates to be obtained from analytically intractable PDFs. The number of sinusoids (model order) and parameters which maximise the posterior PDF is estimated using a reversible jump MCMC (rjMCMC) algorithm. Such techniques have gained popularity in a wide range of subject areas over the past few years, from economics to biology, primarily due to the rapid increase in computation power. The improvements provided by the Bayesian approach compared with classical analysis are significant. This algorithm is used to analyse the spectral characteristics of pulse echo returns from nonlinearly scattering MBs [3–5]. This technique is the starting point for automatic data classification, thereby leading to the creation of a framework for the generation of an active dynamically adaptive pulsing scheme for ultrasound imaging. The outline of this paper is as follows. Section 2 gives a background discussion on existing spectral analysis methods used in ultrasound contrast imaging (UCI). Section 3 proposes a physically inspired parametric model for experimentally measured echo signals from ultrasound MBs and describes the detailed procedure of Bayesian spectral estimation, together with the rjMCMC algorithm. A comparison between parametric estimation and nonparametric estimation when applied to simulated signals is presented in Section 4.1, while in Section 4.2 we apply Bayesian spectral estimation to the analysis of ultrasound MB echoes and displays estimation results. A discussion and conclusions are also presented in Section 5. 2. Spectral Analysis in UCI Ultrasound MBs are stable subcapillary sized microspheres gas-filled and encapsulated in a thin shell, usually with diameter below 7 μm, which can go through microcirculations in the human body [2]. Encapsulation is necessary to prevent rapid dissolution of the gas content into the blood [15]. When exposed to ultrasound, a contrast MB starts to oscillate under the pressure of the sound field. This oscillatory behavior results in a high scattering strength of the contrast MB [16]. Compared to other ultrasound scatterers, including tissue, MBs are more compressible and expandable when insonificated with ultrasound. This nonlinear behavior results in the spectra of the echo signals containing a greater number of harmonics than those originally present in the transmit pulse wave [2, 7, 8]. Therefore, MBs are usually referred to as nonlinear scatterers. Moreover, the nonlinear acoustic signature of MBs forms the basis for discrimination between themselves and other scatterers, for example, soft tissue. Current pulse sequence design research attempts to maximise the difference between responses from MBs and tissue and attempts to increase the contrast-to-tissue ratio (CTR) [17]. An important step prior to pulse-waveform design is therefore an exploration of the spectral content of responses from contrast MBs in the frequency domain from a signal processing perspective; this is the focus of this paper. Traditional spectral estimation techniques in ultrasonics are nonparametric, which are typically based on the Fourier transform (FT). Using these techniques, the spectra are usually derived directly from the observed data, EURASIP Journal on Advances in Signal Processing which is either explicitly or implicitly windowed in time. It is often assumed that the unknown data outside the window is zero or that the entire signal is periodic. In many practical applications, this is an unreasonably strong assumption and the window effect will limit the frequency resolution obtained by the estimator [18]. It would be better to make much weaker assumptions on the nature of the unobserved data. Moreover, for signals with a relatively short duration, which are very common in ultrasound echo signals, the peaks in their Fourier spectra are broad and sufficiently closely spaced that their frequency components can become indistinguishable. Therefore, it is not easy to determine the exact positions of peaks if nonparametric methods are used. For example, various optimisation methods, such as thresholding, need to be considered to determine the accuracy to which the peak frequency can be determined. Furthermore, the number of frequency components is often assumed fixed and known in nonparametric methods. However, this is not practical in analysing echo signals from ultrasound contrast MBs because determining the number of harmonics in ultrasound nonlinear echoes is exactly the question that needs to addressed. In contrast to nonparametric spectral estimation methods, parametric estimation relies on specific signal models which are assumed to generate all data samples, including both observed data and unseen or hidden data. In this case, no windowing occurs and the resolution limitations may be overcome. After a signal model is chosen, Bayesian inference [19] is used to estimate the model parameters. If the parameters that need to be estimated are the frequency components only, this parametric estimation can specifically be called Bayesian spectral estimation [9, 12]. Gregory [12] describes in detail how to apply Bayesian inference to spectral estimation of astronomy signals. Bayesian inference is one popular category of statistical inference. It is able to produce various measures of the estimated parameter by constructing the joint PDF of all unobserved quantities on the basis of all those that are known. Hence, for example, not only can the MMAP estimate or MMSE estimate be found, but the variance or other statistic on the parameter estimate can also be derived. The literature on Bayesian inference is vast, and introductions can be found in [13, 19]. If the posterior PDF of the proposed signal model is multimodal and generally has a complicated shape, it can be difficult to find a closed-form expression for an estimator. However, if the posterior PDF can be approximated by MCMC algorithms [19, 20], it is possible to numerically evaluate the estimators. Moreover, when the number of spectral components is unknown, the joint posterior PDF of the model order and the parameters themselves can be approximated using a rjMCMC technique proposed by Green [14]. This approach is based on reversible samplers that can jump between parameters subspaces of differing dimensionality. Andrieu and Doucet [10] first applied a rjMCMC algorithm to well-constructed sinusoidal signal models, and this work has been extended by others [11, 21]. In parametric estimation methods, the most important decision is finding a proposed signal model that fits the observed data. 3 Otherwise, the estimation of model parameters has no meaning at all. In the ultrasound literature, nonparametric methods are the most popular choices in ultrasound. In this paper, a parametric Bayesian spectral estimation for the analysis of ultrasound MB echoes is presented as an alternative. It has seldom been applied in the ultrasound literature yet has been widely used in other applications. 3. Parametric Bayesian Spectral Analysis In Section 3.1, model structure and model order selection for the acquired echo signals are discussed; detailsof Bayesian spectral estimation technique for the proposed model is given in Section 3.2; a comparison between Bayesian spectral analysis and conventional DFT estimation techniques addressed for simulated signals is given in Section 4.1. 3.1. Selection of an Appropriate Model. In statistical signal modeling the process of signal model design and selection usually consists of choosing a model structure, choosing a model order, and then performing model parameter estimation and model evaluation. The first step of model structure selection depends on observing or understanding the nature of the physical system that generates the signal. The model order can be regarded as another model parameter and thus may be included in the parameter estimation step. An analysis of the physical mechanism that generates ultrasound echo signals indicates that a typical transmit signal is composed of one single frequency component, and the acoustic response is predicted by MB theory to have several more harmonics compared with the transmit pulse. Since the measured ultrasound incident or transmit waveforms are composed of a six-cycle sinusoids, the corresponding echoes from the scatterers are also likely to have the same periodic structure. An analysis of the transmit and receive waveforms suggests that a sum-of-sinusoids (SoSs) signal model is appropriate, and therefore the following model is proposed: Y0 : y(t) = n(t), (2a) Yk : y(t) = s(t) + n(t), (2b) where Y0 and Yk represent a noise-only signal and a signal embedded in Gaussian noise, respectively. The model for the embedded sum-of-sinusoids signal, s(t) for t ∈ {0, . . . , N − 1}, where N is the length of the signal, is given by k s(t) =     ac, j cos ω j t + as, j sin ω j t , (3) j =1 where k is the number of frequencies, or model order, and ac, j and as, j are amplitudes for each frequency k. The noise sequence, n(t), is zero-mean Gaussian noise, given by   n(t) ∼ N 0, σk2 , (4) 4 EURASIP Journal on Advances in Signal Processing where σk2 is the variance. The model can also be written in a vector-matrix form: y = s + n = Dk ak + n, (5) where ak = [ac,1 as,1 ac,2 as,2 · · · ac,k as,k ]T represents an amplitude vector for each frequency component. The matrix Dk contains the frequency information in a signal, which can be defined as ⎡ Dk  ⎤ B(ω1 t1 ) · · · B(ωk t1 ) ⎢ ⎥ ⎢ B(ω t ) · · · B(ω t ) ⎥ 1 2 k 2 ⎥ ⎢ ⎢ ⎥ ⎢ ⎥, .. .. .. ⎢ ⎥ ⎢ ⎥ . . . ⎣ ⎦ B(ω1 tN ) · · · B(ωk tN ) (6) where B(·)  [cos(·), sin(·)] and it operates on an element by element basis; for notational clarity, tn = n − 1. For a fixed model order, k, the harmonic model in (3), is defined by the model parameter set Φk  {ωk , ak , σk2 }, where ωk is the vector of unknown frequencies ωk  {ω1 , . . . , ωk }. The complete set of unknown parameters in the case when the model order is unknown is thus given by Ψ = {k, Φk }. 3.2. Bayesian Spectral Estimation Using a rjMCMC Algorithm. Once a signal model is proposed, a parametric estimation approach using Bayesian inference is applied, taking into account all parameters of interest. According to Bayes theorem [13, 19], the joint distribution of the parameters, Ψ, conditional on the observed data sequence y, p(Ψ | y), is termed the joint posterior density and can be calculated as   p Ψ|y =   p y | Ψ p(Ψ)   , p y (7) 3.2.2. Joint Prior Distribution for Model Parameters, p(Ψ). The joint prior distribution, p(Ψ), reflects the degree of belief of the relevant values of the parameters of a model. It is assumed that the prior distributions for the amplitudes, frequencies, and model order are independent of one another, such that the joint prior distribution simply becomes the product of all priors. Most of these priors are chosen as uninformative Jeffrey’s priors or conjugate priors. The selection of uninformative Jeffrey’s priors reduces the dependence of the results on the exact choice of prior; the selection of conjugate priors allows some unknown parameters to be integrated out analytically and thus reduces the computational complexity. Furthermore, they have been demonstrated to have good performance in [10] when Bayesian inference is applied to spectral estimation. Typical priors for scalar multipliers, such as ak , is a Gaussian density; the inverse-Gamma PDF is used for scalar amplitudes such as variances, σk2 , and the Poisson distribution is used for the discrete-model order, k. These priors themselves then depend on hyperparameters which are usually assumed known. The frequencies ωk are assumed to be uniform over the range [0, π). If 0 < k ≤ kmax , where kmax is a fixed maximum number of frequency, the joint prior can be expressed as p(k, Φk )       = p(k | Λ)p(Λ)p(ωk | k)p ak | σk2 , δ 2 p δ 2 p σk2  ∝  Λk exp[−Λ] k!    poisson      k Λ1 −1/2 exp[−2 Λ]   Gamma (hyperparam )  1 ×   exp − 2πσk  1/2 k   aTk Σ−k 1 ak 2σk2 1 π   uniforms   Gaussian where p(Ψ) is termed as the joint prior density and p(y | Ψ) is the likelihood function. The parameter Ψ is then estimated by, for example, finding the MMAP or MMSE as described in (1). However, the resulting posterior density is often multimodal, and finding analytic expressions for the MMAP or MMSE values becomes complicated. Therefore, the corresponding inference will be approximated using MCMC algorithms, which draw samples from target distribution and find statistical averages using the law of large numbers. 3.2.1. Likelihood Function for the Proposed Model, p(y | Ψ). As the noise is assumed to be independent and identically distributed (i. i. d) Gaussian noise with variance σk2 , the likelihood function based on the chosen signal model for measured ultrasound echo signals discussed in Section 3.1 is obtained as     −N/2 y − Dk ak 2    2 × exp − , (8) p y | k, Φk = 2πσk 2 2σk where A2  AT A, Φk  {ωk , ak , σk2 } contains model parameters, and Ψ = {k, Φk }.  ! ⎧  ⎫ ⎨ exp −γ/ 2σk2 exp −βδ /δ 2 ⎬ × . 2 , 2(ν/2+1) ⎩ δ (αδ + 1) ⎭ σk    inverse-Gammas (variance and hyperparam) (9a) where 1 , 2 are the hyper-hyperparameters of Λ, which itself is the hyperparameter for the model order k, and αδ , βδ are the hyperparameters of δ 2 , which itself is a hyperparameter for the scaler multipliers ak . Moreover, for analytic tractability, Σ−k 1  δ −2 DTk Dk . If there are no harmonic components present in the signal, such that k = 0, the joint prior distribution simplifies to p(Φk | k = 0) % & γ ∝ exp − · 2σk 1 2(ν/2+1) σk · exp[−Λ(1 + 2 )] · Λ1 −1/2 . (9b) Note, the following results are adopted in this calculation: aT0 Σ−0 1 a0  0 and |2πσk2 Σ0 |  1 [10, 11]. EURASIP Journal on Advances in Signal Processing 5 3.2.3. Joint Posterior Distribution for Model Parameters, p(Ψ | y). According to Bayes’ theorem, the corresponding posterior distribution for the unknown parameters can be found, up to a normalising constant, by multiplying (8) and (9a), and then marginalising over the linear parameters ak and the variance term σk2    p k, ωk , Λ, δ 2 | y ∝ γ + yT Pk y −(N+ν)/2 · Λk + 1)π]k k![(δ 2   βδ Λ1 −1/2 e−2 Λ 1 × · 2(α +1) · exp − 2 , eΛ δ δ δ (10) where Pk = IN − Dk Mk DTk , (11) M−k 1 = DTk Dk + Σ−k 1 , In the case of k = 0, the corresponding posterior distribution is presented as:   ' yT y + γ p k, ωk , Λ, δ | y, k = 0 ∝ N+2(ν/2+1) . exp − 2δ02 σ0 2 1 ( × Λ1 −1/2 . exp[−Λ(1 + 2 )]. (12) More details about choosing prior distributions for model parameters and the derivation of joint posterior distributions can be found in [6, 10, 11]. Although the resulting joint posterior distribution has been simplified, it is still not in an analytical form for finding either a MMAP or MMSE estimate. Therefore, statistical sampling techniques are used to construct an empirical approximation of the posterior  y), by sampling distributions, p(k | y) and p(ωk | k, the parameters of interest. Given these empirical approximations, corresponding MMAP estimates of the desired parameters, k and ωk , in the model can then easily be found from   k = arg max p k | y , k    y , ω  k = arg max p ωk | k, (13) ωk where the marginal p(k | y) can be found by simply histogramming the model order samples and ignoring the frequency estimates. Similarly, MMSE can be found from   k = E k | y ,   ω  k = E ωk | k, y , (14) where the ensemble expectation is approximated to the sample average. 3.2.4. A rjMCMC Algorithm for Sampling from the Posterior Distribution. Markov chain Monte Carlo (MCMC) techniques are a set of algorithms that are able to sample from any probability distribution. Reviews can be found in [13, 19, 20]. Using these algorithms, instead of sampling from a probability distribution directly, an ergodic Markov chain is constructed. After a number of iterations, the chain will become stable and its equilibrium distribution can be regarded as the desired probability distribution. Having obtained samples from the desired distributions, an empirical PDF can be constructed using simple histograms and a MMAP criterion adopted to obtain the mode of the  y). estimated posterior distribution p(k | y) and p(ωk | k, Since the model order k is unknown, a reversible jump MCMC technique proposed in [14] must be incorporated in order to search the joint (variable) dimensional state-space. The rjMCMC method combines model order selection and parameter estimation and regards the model order as one of the unknown model parameters. This technique allows the proposals to jump between subspaces of different dimensions and visit all relevant model orders. With a certain ratio that ensures the reversibility condition, and hence the invariance of a Markov chain, samples from the proposal distribution are accepted [19]. Due to this calculated acceptance ratio, the most likely model orders are visited most often and the least likely ones are visited with lower probability. In this way, the computational complexity can be reduced. In general, there are three simple and commonly used candidate moves for frequency estimation: birth, death, and update moves. Birth and death moves are widely used complementary moves. In a birth move, the algorithm proposes a candidate of higher dimension whereas in a death move, the algorithm proposes a candidate in the model of lower dimension. If neither move is chosen, an update move will be carried out. Let bk , dk , and uk denote the probabilities of each move for frequency estimation, satisfying the relationship of bk + dk + uk = 1 for all k. In each iteration step, the three moves can be calculated as follow: ⎧   ⎪ ⎪ ⎨c min 1, p(k + 1) , p(k) bk = ⎪ ⎪ ⎩0, ⎧ ⎪ 0, ⎪ ⎨   dk = ⎪ p(k − 1) ⎪ , ⎩c min 1, p(k) 0 ≤ k < kmax k = kmax k=0 (15) 1 ≤ k ≤ kmax where p(k) represents the prior distribution of the model parameter k [10, 11]. The constant c is a tuning factor which determines the ratio of the update move to the jump moves; c = 0.5 is chosen so that the probability of a jump move is between 0.5 and 1 at every iteration. This choice also ensures bk p(k) = dk+1 p(k+1), which could guarantee certain acceptance in the corresponding Metropolis-Hastings (MH) sampler [14]. The procedure of approximation of Bayesian spectral estimation is summarized in Algorithm 1. In this algorithm, the acceptance probabilities of birth move, death move, and update move, denoted as αB , αD , and αU , respectively, can be obtained following the MH sampling 6 EURASIP Journal on Advances in Signal Processing (1) Initialization: set ({k, Φk }(0) ) (2) Iteration: (3) for i = 1 to numIteration do (4) Update each frequency component: (5) — Sample hyperparameters Λ and δ 2 . (6) — Sample u from U(0,1) . (uniform distribution) (7) if u ≤ bk(i) then (8) perform birth move: Propose a new frequency randomly on [0, π) and accept it with a probability of αB (9) else if(u ≤ bk(i) + dk(i) )then (10) perform death move: Remove an existing frequency randomly from ωk and accept it with a probability of αD (11) else (12) perform update move: Update for all k frequencies according to a proposal distribution and accept it with a probability αU (13) end if (14) — Sample nuisance parameters ak and σk2 (15) end for Algorithm 1: A rjMCMC algorithm for frequency estimation. procedure, which is described in detail in [19]. The general form is as follows: *  α  min 1, posterior distribution ratio   + × proposal distribution ratio , (16) where the posterior distribution ratio is the same for three different moves and can be derived from (10). The proposal distribution ratio is different according to different moves. For the update move, the selection of proposal distribution follows that described in [10]. In summary, there are two proposal distributions for estimating the frequencies present in the signal, namely: (1) q1 for a global distribution which occurs with probability λ, (2) q2 for a local distribution which occurs with probability 1 − λ. Together, these proposals form a hybrid mixture kernel in the MH algorithm [19, 20]. When proposing new frequency estimates, the DFT spectrum of the observed signal is used to define the global proposal distribution up to a normalising constant, and a random walk perturbation with a variance of 2 is taken as the local proposal distribution. The factor λ, σRW which determines the ratio of how often q1 or q2 is sampled from, is determined in a heuristic way and set to λ = 0.2 [10]. However, if there are a large number of frequency components in a very short period in the time domain, using the DFT to define the proposal distribution tends to overestimate the number of frequencies. Therefore, rather than choosing DFT spectrum, the multitaper power spectrum is adopted. This reduces the variance without losing the frequency resolution and thus provides a clearer spectrum [22]. With respect to the proposal distribution for birth move and death move, they are related to the probabilities of birth move, bk , and death move, dk , respectively. When a birth move is performed, the probability of proposing a new frequency component on [0, π) is 1/π. When a death move is performed, the probability of choosing one frequency component from the existing frequencies is 1/(k + 1). Therefore, the proposal distributions for a birth move starting from k frequencies can be expressed as bk /π, while a death move starting from k + 1 frequencies is dk+1 /(k + 1), respectively. 4. Results In order to increase the penetration depth in conventional pulsed ultrasound contrast imaging systems, ultrasound excitation pulses usually have short-time duration and large bandwidth in the frequency domain. The performance of nonparametric spectral estimation methods depends on the length of a signal in the time domain since, as stated by the time-frequency uncertainty principle, there is a fundamental trade-off between temporal and spectral localisation. The DFT spectra of the short-duration broadband excitation pulses exhibit low-frequency resolutions and thus DFTbased methods are not adequate for analysis of ultrasound MB echoes where high-frequency resolution is required. In contrast, the frequency resolution using Bayesian spectral estimation in a parametric framework mainly depends on the SNR of a signal, in addition to the signal duration in samples [12]. Section 4.1 compares the Bayesian spectral estimation method of Section 3.2 with the DFT estimate. In Section 4.2 we apply Bayesian spectral estimation to backscattered MB echo signals. 4.1. Comparison to Nonparametric Spectral Estimation Methods. In order to compare traditional nonparametric spectral estimation and parametric Bayesian spectral estimation with respect to frequency resolution, an example signal with two closely spaced frequency components is provided. EURASIP Journal on Advances in Signal Processing 7 35 Magnitude response in log 140 Magnitude response in log Simulated signal with two different frequencies (10 cycles) Amplitude 20 0 −20 0 5 10 15 20 25 30 Simulated signal in the frequency domain (10 cycles) 102 0 1 2 Time Simulated signal with two different frequencies (35 cycles) Amplitude 20 0 −20 0 20 40 60 80 100 120 3 4 Radial frequency 5 6 Simulated signal in the frequency domain (35 cycles) 104 102 0 1 2 Time 3 4 Radial frequency (a) 5 6 (b) Figure 1: Two simulated signals with different lengths in the time and frequency domain. (a) Two simulated signals with different lengths in the time domain. (b) Two simulated signals with different lengths in the frequency domain. Table 2: Parameters of the algorithm. Table 1: Parameter settings for the simulated signal with two closely spaced frequencies. Component k 1 2 Amplitude Ek a2s,k + a2c,k Phase −tan−1 (as,k /ac,k ) Frequency ωk 20 20 0 π/4 0.6π 0.62π This signal is simulated according to (2b) and (3), with the parameter setting displayed in Table 1. With these parameters, Figure 1(a) displays two signals with different lengths, which can also be modelled by sinusoidal cycles. The first one in the upper panel consists of 32 data points, which corresponds to 10 cycles of a sinusoidal signal. It is also a typical example of broadband ultrasound echo signals. The second one in the lower panel has 128 data points, for example, 35 cycles, which is much longer than most pulses in ultrasound echo signals. The SNRs of two simulated signals are both 5 dB. Figure 1(b) also shows their corresponding frequency spectra using DFT estimation. As discussed in any standard text on spectral analysis, for example [22], the frequency resolution using DFT can be approximated by Δω ≈ 2π/(N − 1), in which N represents the number of data points in a signal. Therefore, in order to differentiate two frequency components in this example, the signal length must satisfy N − 1 > 2π/(Δω) = 2π/(0.62π − 0.6π) = 100. This indicates that only if the number of data points in a signal is larger than 101, can the DFT discriminate two frequencies with frequency difference of 0.02π. This result is demonstrated in Figure 1(b). It can be clearly seen that the upper panel displays only one obvious peak in the spectrum of the signal with 10 cycles, whereas the lower panel Parameters ν γ 1 2 αδ βδ λ σRW Values 0 0 0.01 0.01 2 10 0.2 1/(5N) exhibits two obscure peaks in the spectrum of the signal with 35 cycles. The same examples of the simulated signals are now analysed using the proposed Bayesian spectral estimation. The parameter setup for the MCMC algorithm, as presented in (10), is displayed in Table 2. Note that ν and γ, which are the hyperparameters of the noise variance σk2 , are both set to zero in order to ensure that an uninformative Jeffrey’s prior is selected, that is, p(σk2 ) = 1/σk2 . Moreover, Λ and δ 2 are hyperparameters of the number of frequency k; 1 and 2 are hyperparameters of Λ; αδ and βδ are hyperparameters of δ 2 . Therefore, 1 , 2 , αδ , and βδ are hyper-hyperparameters of the signal model. It has been demonstrated in [10] that the estimation procedure is very insensitive to the specification of the hyper-hyperparameters when it is applied to frequency estimation. Their values are assigned according to the discussion in [10]. Other parameters, such as λ and σRW , are chosen in a rather heuristic way: λ = 0.2 and σRW = 1/(5N ), where N is the signal length, are used in [10], and achieve good estimated results. Therefore, they are also adopted in this example. In this way, the number of frequency components in the simulated signal can be estimated by p(k | y). In order to compare the performance of nonparametric DFT estimation and parametric Bayesian spectral estimation when applied to a signal with short-time duration, for example, a sinusoidal signal with fewer than 15 cycles, 8 EURASIP Journal on Advances in Signal Processing Estimates (Freq.) 1st 2nd 3rd 4th 5th 6th 7th 8th Means (MHz) 1.8081 2.1639 2.4560 2.7567 2.9575 3.2091 3.4039 3.5924 SD values (MHz) 0.0951 0.2550 0.2388 0.2198 0.1949 0.1573 0.1629 0.1447 the probabilities of correct detection of the number of frequencies present can be calculated. This is equivalent to evaluating their ability of discriminating two closely spaced frequencies. Figure 2 depicts the probabilities of correct detection of the number of frequencies with an increase of the number of sinusoidal cycles using Bayesian spectral estimation and DFT estimation. In this figure, there are three different SNRs, 5 dB, 10 dB, and 20 dB. It is seen that the DFT estimation is independent of SNR, whereas Bayesian spectral estimation highly depends on SNR. Moreover, for a sinusoidal signal with a length shorter than 15 cycles, the DFT estimation cannot correctly estimate its frequency numbers, no matter what the SNR is. However, Bayesian spectral estimation outperforms DFT estimation when the sinusoidal signal has a length shorter than 15 cycles. Furthermore, for a sinusoidal signal with a SNR higher than 20 dB, only 10 cycles are required to estimate frequency numbers correctly with a 100% percentage. Ultrasound echo signals from MBs usually have fewer sinusoidal cycles and relatively high SNRs. Therefore, with respect to the spectral analysis of ultrasound MB echoes, the parametric Bayesian spectral estimation is more appropriate than the traditional nonparametric DFT estimation. 4.2. Analysis of Backscattered MB Signals. The experimentally measured echo signals from ultrasound MBs used in this paper are all acquired by a modified ultrasound transducer (Sonos5500 Philips Medical Systems, Andover, MA, USA), which is performed in an in vitro environment. This transducer transmits an ultrasound wave and receives its reflection from scatterers. The single MBs used in this paper are Definity (Bristol-Myers Squibb Inc, MA, USA), which were lipid coated and the gas was perfluorocarbon. The mean diameter of the MBs ranges from 1.1 μm to 3.3 μm. They are released from a micropipette with a diameter of 100 μm, which is put at the bottom of a water tank, as described in [4, 5]. The tank is filled with anechoic water, which does not produce echoes. The raw echo signals measured from these MBs are preamplified, collected, and digitally stored and are the subject of spectral analysis in this section. As mentioned in Section 1, ultrasound transmit signals usually have a single fundamental frequency whereas their echo signals from MBs usually have more harmonics. Figure 3(a) displays a typical example of measured ultrasound MB echoes. It has a length of 80 data points, which is Table 4: Spread of MAP estimates of number of frequency for the measured MB response. 0∼6 0 K MAP Probability of correct detection of frequency numbers Table 3: Estimated frequency values for the measured MB response. 7 6 8 80  10 9 14 0 Comparison between Bayesian and DFT spectral estimation against SNR 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 5 10 15 20 25 Number of cycles Bayesian (20 dB) FFT (20 dB) Bayesian (10 dB) 30 35 FFT (10 dB) Bayesian (5 dB) FFT (5 dB) Figure 2: Comparison of probability of correct detection of frequency numbers using Bayesian inference and DFT estimation. relatively short duration in the time domain. Moreover, its corresponding DFT spectrum estimation is shown in Figure 3(b). The peaks in this figure are obscured and need to be explored further to determine whether there is a frequency component present or not. When applying Bayesian spectral estimation to this typical MB echo signal, results are obtained by the analysis of 100 replications of Bayesian inference on the same MB echo signal; each run of the simulation, which involves a finite number of samples in the rjMCMC algorithm, will be subject to statistical variations due to the stochastic sampling strategy. The spread of MMAP values of the number of frequency for the pulse is shown in Table 4. The results indicate that for 80% out of 100 replications of the estimation algorithm, the estimated number of frequency component is 8 for the measured MB response, which perhaps is not immediately obvious from the DFT in Figure 3(b). Moreover, the correspondingly estimated frequency values are presented in Table 3. These estimates are given with both means and standard deviation (SD) values. 5. Discussion and Conclusions This paper introduces Bayesian spectral estimation for the analysis of echo signals from ultrasound contrast MBs. It assumes a parametric model and estimates model parameters EURASIP Journal on Advances in Signal Processing 9 Frequency spectrum of the 4th pulse in the echo signal from MB Multiple pulse echo signal from MB 103 8 Magnitude response in log scale Pressure (Pascal) 6 4 2 0 −2 −4 −6 5.5 5.6 5.7 5.8 5.9 Time (seconds) 6 6.1 ×10−5 102 101 100 10−1 10−2 0 (a) 2 4 6 Frequency (Hz) 8 10 ×106 (b) Figure 3: Display the MB echo with transmit frequency at 1.83 MHz in the time and the frequency domains. (a) Enlarged version of the 4th pulse in the MB echo. (b) Frequency spectrum of the 4th pulse in the MB echo. within a Bayesian framework. As the posterior density function cannot be solved in a closed form and the dimension of the parameters changes, a reversible jump MCMC algorithm is adopted to approximate Bayesian inference. This Bayesian spectral estimation is also compared to traditional nonparametric DFT estimation. With respect to measured ultrasound MB echo signals, which usually have short-time duration and high SNRs, Bayesian spectral estimation is able to differentiate two closely spaced frequency components and thus correctly estimate frequency numbers and their values. Furthermore, these attributes detected using a parametric method can lead to an extension to discriminating echo signals from MBs and soft tissue, which is not within the scope of this paper but is a prosperous field for MB behavior analysis. Despite the advantages of this newly introduced Bayesian spectral estimation in a parametric framework, its limitation lies in that the proposed signal model should fit the measured ultrasound MB responses well. In this paper, rather than being justified in a strict sense, the proposed model is only demonstrated from the perspective of physical understanding of MB characteristics. Therefore, more investigations of model validation are needed in the future work. Moreover, the computational complexity is increased when Bayesian spectral estimation is used. As a result, Bayesian spectral estimation applied to ultrasound contrast MB responses reveals more characteristics in the frequency domain, which may broaden the research field in ultrasound contrast agents. Acknowledgments The authors would like to thank Biological Inspired Acoustic Systems (BIAS) project for providing the financial support. Vassilis Sboros is supported by British Heart Foundation (BHF). The Institute of Digital Communications is a member of a Joint Research Institute with the HeriotWatt University, which is a part of the Edinburgh Research Partnership in Engineering and Mathematics (ERPem), and acknowledges support from the Scottish Funding Council. References [1] R. Gramiak and P. M. Shah, “Echocardiography of the aortic root,” Investigative Radiology, vol. 3, no. 5, pp. 356–366, 1968. [2] B. B. Goldberg, Ed., Ultrasound Contrast Agents, Martin Dunitz, London, UK, 1997. [3] V. Sboros, S. D. Pye, T. Anderson et al., “The measurement of backscatter from individual contrast agent microbubbles,” in Proceedings of IEEE Ultrasonics Symposium, pp. 1945–1947, October 2002. [4] V. Sboros, S. D. Pye, C. A. MacDonald, J. Gomatam, C. M. Moran, and W. N. McDicken, “Absolute measurement of ultrasonic backscatter from single microbubbles,” Ultrasound in Medicine and Biology, vol. 31, no. 8, pp. 1063–1072, 2005. [5] D. H. Thomas, M. B. Butler, A. Dermitzakis, T. Anderson, W. N. McDicken, and V. Sboros, “The acoustic scatter from single bisphere microbubbles,” Ultrasound in Medicine and Biology, vol. 36, no. 11, pp. 1884–1892, 2010. [6] Y. Yan, Statistical signal processing for echo signals from ultrasound linear and nonlinear scatterers, Ph.D. thesis, University of Edinburgh, Edinburgh, Scotland, November 2009. [7] T. Albrecht and J. Hohmann, “Contrast agents in sonography,” in Toshiba VISIONS, 2004. [8] E. Stride and N. Saffari, “Microbubble ultrasound contrast agents: a review,” Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, vol. 217, no. 6, pp. 429–447, 2003. [9] G. Bretthorst, Bayesian Spectrum Analysis and Parameters Estimation, Springer, New York, NY, USA, 1988. [10] C. Andrieu and A. Doucet, “Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC,” IEEE Transactions on Signal Processing, vol. 47, no. 10, pp. 2667–2676, 1999. 10 [11] K. Copsey, N. Gordon, and A. Marrs, “Bayesian analysis of generalized frequency-modulated signals,” IEEE Transactions on Signal Processing, vol. 50, no. 3, pp. 725–735, 2002. [12] P. C. Gregory, “A Bayesian revolution in spectral analysis,” in Bayesian Inference and Maximum Entropy Methods in Science and Engineering, 2001. [13] A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin, Bayesian Data Analysis (Texts in Statistical Science), Chapman & Hall/CRC, Boca Raton, Fla, USA, 2004. [14] P. J. Green, “Reversible jump MCMC computation and Bayesian model determination,” Biometrika, vol. 82, pp. 711– 732, 1995. [15] A. Bouakaz and N. de Jong, “New contrast imaging method using double frequency exposure,” in Proceedings of the IEEE Ultrasonics Symposium, pp. 339–342, August 2004. [16] S. Hilgenfeldt, D. Lohse, and M. Zomack, “Sound scattering and localized heat deposition of pulse-driven microbubbles,” Journal of the Acoustical Society of America, vol. 107, no. 6, pp. 3530–3539, 2000. [17] N. de Jong, R. Cornet, and C. T. Lancée, “Higher harmonics of vibrating gas-filled microspheres. Part one: simulations,” Ultrasonics, vol. 32, no. 6, pp. 447–453, 1994. [18] P. M. Clarkson, Optimal and Adaptive Signal Processing, CRC Press, Boca Raton, Fla, USA, 1993. [19] D. Gamerman, Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference, Chapman & Hall/CRC, Boca Raton, Fla, USA, 1997. [20] C. Andrieu, N. de Freitas, A. Doucet, and M. I. Jordan, “An introduction to MCMC for machine learning,” Machine Learning, vol. 50, no. 1-2, pp. 5–43, 2003. [21] J.-R. Larocque and J. P. Reilly, “Reversible jump MCMC for joint detection and estimation of sources in colored noise,” IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 231– 240, 2002. [22] D. Manolakis, V. Ingle, and S. Kogon, Statistical and Adaptive Signal Processing: Spectral Estimation, Signal Modelling, Adaptive Filtering and Array Processing, McGraw-Hill, New York, NY, USA, 2003. EURASIP Journal on Advances in Signal Processing