COMPRESSED BEAMFORMING APPLIED TO B-MODE ULTRASOUND IMAGING

Diagnostic sonography allows visualization of body tissues, by radiating them with acoustic energy pulses, transmitted from a transducer array. Applying appropriate delays to the transmitting elements, in a process known as “beamforming”, the interfering wave pattern forms a narrow beam, along which most energy propagates. As it propagates in the tissue, echoes are scattered by density and propagation-velocity perturbations [1]. These reflections are detected by the transducer array in a reception phase, where a second beamforming process is applied, aimed at localizing the scattering structures while improving Signal to Noise Ratio (SNR) [2]. The reception beamforming process involves averaging the detected signals, after aligning them with time-varying delays. Ultrasound systems perform this process digitally, requiring that the analog signals first be sampled. Confined to classic Nyquist-Shannon sampling theorem [3], these systems must sample the signals at twice their baseband bandwidth, in order to avoid aliasing.

The goal of our work is to produce ultrasound images from samples taken at rates far below these dictated by the classical framework. A first step in this direction was proposed in [4], where the signals detected in each element were assumed to satisfy the FRI property [5]. Applying the FRI framework, [4, 6] derive schemes, which allow reconstruction of ultrasound signals from samples taken at sub-Nyquist rates. These approaches follow the spirit of Xampling [7], a term which refers to the combination of classical sampling theory with recent developments in Compressed Sensing (CS). However, applying such schemes to data reflected from biological tissues, the resulting signal parameters are often erroneous. This is mainly due to the noisy nature of multiple reflections, also known as “speckles” [2], which are scattered by microscopic perturbations in the tissue.

In order to cope with the noisy nature of the signals detected in individual elements, [8] first introduces a multi-sensor Xampling scheme, which enhances the SNR of the low-rate samples. This is achieved, by integrating traditional beamforming with FRI Xampling. The approach proposed in [8] is farther developed in [9], resulting in a more straightforward sampling mechanism. Both approaches, referred to by the term “Compressed Beamforming”, are derived and simulated for the special case of linear, B-mode scan.

In the current work, we generalize the method introduced in [8, 9] to the case of polar scan. Additionally, we propose reconstructing the sampled signal by applying CS techniques, rather than spectral analysis methods, traditionally used in the FRI framework [5]. Finally, we apply our approach to real cardiac ultrasound data, where we successfully image macroscopic perturbations, while achieving a nearly eight-fold reduction in sampling rate, compared to standard imaging techniques.

This paper is organized as follows: in Section 2, we outline the principles of beamforming in ultrasound imaging. In Section 3 we discuss the FRI model in the context of ultrasound imaging, and its contribution to rate reduction when sampling signals at individual receivers. Our multiple-receiver Xampling scheme, which integrates beamforming into the Xampling process, is presented in Section 4. In Section 5 we introduce experimental results, obtained for cardiac ultrasound data.

We begin by outlining the beamforming process, lying at the heart of ultrasound imaging. The process is carried out both in transmission and reception. During transmission, modulated pulses are transmitted from all active array elements. Applying appropriate delays to the transmitting elements, the propagating energy may be directed along a narrow beam, in any desired direction, the result of a constructive and destructive interference pattern. The acoustic reciprocity theorem allows to perform a similar process during reception. In this case, we regard the perturbations, which reflect the transmitted energy, as point transmitters. Applying appropriate delays to the detected signals, the array may be focused to any coordinate, such that if, indeed, energy was scattered from that point, then the combined signals will undergo constructive interference. Knowing that most energy propagates along a narrow beam, and knowing its propagation velocity, we may associate time instances to radiated points in the medium, thus optimizing the reception focal point at these instances. The resulting process, in which time-varying delays are applied to the detected signals, prior to combining them for SNR enhancement, is referred to as “Dynamic Focusing”.

In order to derive an explicit expression for dynamic focusing, consider an array of $M$ receivers, aligned along the $\hat{x}$ axis. We denote by $\delta_{m}$ the distance from the $m$’th receiver to the origin. Assume a pulse of energy, transmitted at $t=0$ from the origin, at direction $\theta$. The setup is depicted in Fig. 1. At time $t$, the pulse, propagating at velocity $c$, crosses $\left(x,z\right)=\left(ct\cos\theta,ct\sin\theta\right)$. A reflection originating from this coordinate will reach the $m$’th element at time:

$$\tau_{m}(t)=t+\sqrt{\left(t\cos\theta\right)^{2}+\left(\gamma_{m}-t\sin\theta\right)^{2}},$$

(1)

with $\gamma_{m}=\delta_{m}/c$. Denoting by $m_{0}$ the index of a reference receiver, satisfying $\delta_{m_{0}}=0$, it is readily seen that $\tau_{m_{0}}\left(t\right)=2t$. Hence, in order to align the reflection detected in the $m$’th receiver with the one detected in the reference receiver, we need to apply a delay to the detected signal, $\varphi_{m}\left(t\right)$, such that the resulting signal, ${\hat{\varphi}}_{m}\left(t;\theta\right)$, satisfies ${\hat{\varphi}}_{m}\left(2t;\theta\right)=\varphi_{m}\left(\tau_{m}\left(t\right)\right)$. This implies a time-varying time delay, applied to $\varphi_{m}\left(t\right)$, resulting in the following distorted signal:

$\displaystyle\begin{split}{\hat{\varphi}}_{m}\left(t;\theta\right)&=\varphi_{m}\left(\frac{1}{2}\left(t+\sqrt{t^{2}-4\gamma_{m}t\sin\theta+4\gamma_{m}^{2}}\right)\right).\end{split}$

(2)

The aligned signals may be combined into the beamformed signal

$$\Phi\left(t;\theta\right)=\frac{1}{M}\sum_{m=1}^{M}{{\hat{\varphi}}_{m}\left(t;\theta\right)},$$

(3)

which exhibits enhanced SNR compared to $\left\{{\hat{\varphi}}_{m}\left(t;\theta\right)\right\}_{m=1}^{M}$.

Although defined over continuous time, ultrasound systems perform the process formulated in (2)-(3) in the digital domain, requiring that the analog signals $\left\{\varphi_{m}\left(t\right)\right\}_{m=1}^{M}$ first be sampled. Confined to classic Nyquist-Shannon sampling theorem, these systems sample the signals at twice their baseband bandwidth, in order to avoid aliasing. The detected signals typically occupy only a portion of their baseband bandwidth. Exploiting this fact, some state of the art systems manage to reduce the amount of samples transmitted from the front-end, by down-sampling the data, after demodulation and low-pass filtering. However, since such operations are carried out digitally, the preliminary sampling-rate remains unchanged.

Our goal is to produce an ultrasound image from samples taken far below the Nyquist rate. A first step in this direction was proposed in [4], where the signals detected in each element were assumed to satisfy the FRI property. More specifically, they were approximated as $L$ replicas of a known-shape pulse, $h\left(t\right)$, caused by scattering of the transmitted pulse from $L$ reflectors, located along the narrow transmission beam. Explicitly, the signal detected in the $m$’th receiver was modeled as

$$\varphi_{m}\left(t\right)=\sum_{l=1}^{L}{a_{l,m}h\left(t-t_{l,m}\right)},$$

(4)

where $t_{l,m}$ denotes the time, in which the reflection, originating from the $l$’th reflector, arrives at the $m$’th receiver, and $a_{l,m}$ denotes the reflection amplitude. Hence, the detected signal is completely defined by $2L$ parameters, $\left\{t_{l,m},a_{l,m}\right\}_{l=1}^{L}$. As shown in [5], spectral analysis techniques (e.g. annihilating filter [10]) may be applied, in order to extract the signal parameters, from its projection onto a $2L$-dimensional subspace, corresponding to its Fourier series coefficients. The sampling scheme is thus reduced to the problem of extracting $K$ frequency samples of the detected signal, where $K\geq 2L$. Denote the temporal support of the detected signal, $\varphi_{m}\left(t\right)$, by $\left[0,T\right)$. Imaging to a nominal depth of $16\mbox{cm}$, $T$ is approximately $210\mu\mbox{sec}$. The signal’s baseband bandwidth requires a nominal sampling rate of $f_{s}=16\mbox{Mhz}$, resulting in an overall number of $Tf_{s}=3360$ real-valued samples. Still within classical sampling framework, assuming that the signal’s passband bandwidth is only $4\mbox{MHz}$, the data sampled at Nyquist rate may be down-sampled to 1680 real-valued samples. Thus, as long as the number of reflected pulses satisfies $L\ll 1680$, substantial rate reduction may be achieved by an FRI Xampling scheme.

The scheme proposed in [4] enables sub-Nyquist sampling of the signal detected in an individual array element. Sampling signals reflected from phantom targets, it managed to achieve good signal reconstruction. However, when applied to data reflected from biological tissues, the resulting signal parameters are often erroneous. This is mainly the result of multiple reflections, scattered from microscopic perturbations, manifested as noise. In extreme scenarios, where the noise masks the FRI component, the extracted parameters will be meaningless, such that any attempt to cope with errors in the parametric space will turn out useless.

The novelty of our approach regards the assumption, that the FRI scheme applied to an individual element may be equally applied to the beamformed signal defined in (3), which obviously exhibits much better SNR. Indeed, it may be shown that, if the signals $\left\{\varphi_{m}\left(t\right)\right\}_{m=1}^{M}$ obey the FRI model (4), then $\Phi\left(t;\theta\right)$ is approximately of the form:

$$\Phi\left(t;\theta\right)=\sum_{l=1}^{L}{b_{l}h\left(t-t_{l}\right)},$$

(5)

and may thus be regarded as FRI. The parameters $\left\{t_{l},b_{l}\right\}_{l=1}^{L}$ in (5) depend on $\theta$, yet we omit this dependency from our notation, for convenience. An obvious problem is that, as stated in Section 2, the signal $\Phi\left(t;\theta\right)$ does not exist in the analog domain. In the next section, we derive our scheme, which manages to estimate its necessary low-rate samples from low-rate samples of filtered versions of $\left\{\varphi_{m}\left(t\right)\right\}_{m=1}^{M}$.

The FRI Xampling scheme applies spectral analysis techniques, in order to estimate the $2L$ unknown signal parameters, from a subset of its consecutive Fourier coefficients. Denoting the cardinality of this subset by $K$, perfect reconstruction requires that $K\geq 2L$.

Denote by $\kappa$ the set of indices $\left\{k_{j}\right\}_{j=1}^{K}$, corresponding to the selected Fourier coefficients. Having obtained these coefficients, other mechanisms may be applied, in order to reconstruct the signal. In our work, for instance, we use CS methodology. Since each of the signals $\left\{\varphi_{m}\left(t\right)\right\}_{m=1}^{M}$ is supported on $\left[0,T\right)$, then, by (2)-(3), $\Phi\left(t;\theta\right)$ may only be constructed on a support contained in $\left[0,T\right)$. Considering its Fourier series, calculated with respect to this interval, let us explicitly write the $k_{j}$’th Fourier coefficient, $c_{j}$. We will additionally assume, that the time delays, $\left\{t_{l}\right\}_{L=1}^{L}$, are quantized with $\Delta_{s}$ quantization step, such that $t_{l}=q_{l}\Delta_{s}$, $q_{l}\in\mathbb{Z}$:

	$\displaystyle c_{j}$	$\displaystyle=\frac{1}{T}\int_{0}^{T}{\sum_{l=1}^{L}{b_{l}h\left(t-t_{l}\right)}e^{-j\frac{2\pi}{T}k_{j}t}dt}$
		$\displaystyle=\frac{1}{T}\sum_{l=1}^{L}{b_{l}\int_{0}^{T}{h\left(t-q_{l}\Delta_{s}\right)}e^{-j\frac{2\pi}{T}k_{j}t}dt}$
		$\displaystyle=\frac{1}{T}H\left(\frac{2\pi}{T}k_{j}\right)\sum_{l=1}^{L}{b_{l}e^{-j\frac{2\pi}{T}\Delta_{s}k_{j}q_{l}}}.$		(6)

Here $H\left(\omega\right)$ denotes the CTFT of $h\left(t\right)$. Let $N$ be the ratio $\lfloor T/\Delta_{s}\rfloor$. Then (4) may be brought into the following matrix form:

$${\bf{c}}=\frac{1}{T}{\bf{HAx}},$$

(7)

where ${\bf{H}}$ is the $K\times K$ diagonal matrix with $H\left(\frac{2\pi}{T}k_{j}\right)$ as its $j$’th element, and ${\bf{x}}$ is a length $N$ vector, whose $j$’th element equals $b_{l}$ for $j=q_{l}$, and $0$ otherwise. Finally, ${\bf{A}}$ is a $K\times N$ matrix, formed by taking the set $\kappa$ of rows from a $N\times N$ FFT matrix. Assuming that $L\ll N$, this is a CS problem, aimed at estimating the $L$-sparse vector ${\bf{x}}$ from $K$ samples, taken with a measurement matrix $\frac{1}{T}{\bf{HA}}$. A common technique for recovering ${\mathbf{x}}$ is the Orthogonal Matching Pursuit (OMP) algorithm, which we use in our work.

We now return to the problem of obtaining the set $\kappa$ of $\Phi\left(t;\theta\right)$’s Fourier coefficients. Let $c_{j}$ denote the $j$’th Fourier coefficient of $\Phi\left(t;\theta\right)$. Then,

$\displaystyle\begin{split}c_{j}&=\frac{1}{M}\sum_{m=1}^{M}{\frac{1}{T}\int_{0}^{T}{e^{-j\frac{2\pi}{T}k_{j}t}{\hat{\varphi}}_{m}\left(t;\theta\right)dt}}=\frac{1}{M}\sum_{m=1}^{M}{c_{j,m}},\end{split}$

(8)

where, from (2):

$\displaystyle\begin{split}c_{j,m}&=\\ &\frac{1}{T}\int_{0}^{T}{e^{-j\frac{2\pi}{T}k_{j}t}\varphi_{m}\left(\frac{1}{2}\left(t+\sqrt{t^{2}-4\gamma_{m}t\sin\theta+4\gamma_{m}^{2}}\right)\right)dt}.\end{split}$

(9)

Applying a change of variable to the integral,we can write (9) as

$\displaystyle c_{j,m}$

$\displaystyle=\frac{1}{T}{\int_{0}^{T}g_{j,m}(t;\theta)\varphi_{m}\left(t\right)dt},$

(10)

with

$\displaystyle\begin{split}g_{j,m}(t;\theta)=&q_{j,m}(t;\theta)e^{-i\frac{2\pi}{T}k_{j}t},\\ q_{j,m}(t;\theta)=&I_{\left[|\gamma_{m}|,T_{m}\left(\theta\right)\right)}\left(t\right)\left(1+\frac{\gamma_{m}^{2}\cos^{2}\theta}{\left(t-\gamma_{m}\sin\theta\right)^{2}}\right)\times\\ &\exp\left\{i\frac{2\pi}{T}k_{j}\frac{\gamma_{m}-t\sin\theta}{t-\gamma_{m}\sin\theta}\gamma_{m}\right\}.\end{split}$

(11)

Here $I_{\left[a,b\right)}\left(t\right)$ denotes an indicator function having the value 1 for $a\leq t<b$ and 0 otherwise, and $T_{m}\left(\theta\right)$ is the maximal time value which is mapped, under the temporal distortion defined in (2), to the support of the beamformed signal.

The process defined in (10)-(11) can be translated into a multi-channel Xampling scheme, in which each signal, $\varphi_{m}\left(t\right)$, is multiplied by a bank of kernels $\left\{g_{j,m}\left(t;\theta\right)\right\}_{j=1}^{K}$ defined by (11), and integrated over $\left[0,T\right)$. This results in a vector ${\bf{c_{m}}}=\left[\begin{array}[]{cccc}c_{1,m}&c_{2,m}&...&c_{K,m}\end{array}\right]^{T}$. The vectors $\left\{{\bf{c_{m}}}\right\}_{m=1}^{M}$ are then averaged in ${\bf{c}}=\left[\begin{array}[]{cccc}c_{1}&c_{2}&...&c_{K}\end{array}\right]^{T}$, which has the desired improved SNR property, and provides a basis for extracting the $2L$ parameters which define $\Phi\left(t;\theta\right)$. Nevertheless, the complexity of the analog kernels, together with their dependency on $\theta$, makes hardware implementation of such a scheme extremely complicated. We thereby take an additional step, approximating the samples $\left\{c_{j,m}\right\}_{j=1}^{K}$ from low-rate samples of $\varphi_{m}\left(t\right)$, taken in a much more straightforward manner.

First, we substitute $\varphi_{m}\left(t\right)$ of (10) by its Fourier series, calculated with respect to $\left[0,T\right)$. Denoting the $n$’th Fourier coefficient by $\phi_{m}\left[n\right]$, we get:

$\displaystyle\begin{split}c_{j,m}&=\sum_{n}{\phi_{m}\left[n\right]\frac{1}{T}\int_{0}^{T}{q_{j,m}(t;\theta)e^{-i\frac{2\pi}{T}\left(k_{j}-n\right)t}dt}}\\ &=\sum_{n}{\phi_{m}\left[k_{j}-n\right]Q_{j,m;\theta}\left[n\right]},\end{split}$

(12)

where $Q_{j,m;\theta}\left[n\right]$ are the Fourier series coefficients of $q_{j,m}(t;\theta)$, also defined on $\left[0,T\right)$. Replacing the infinite summation of (12) by its finite approximation:

$\displaystyle\begin{split}{\hat{c}}_{j,m}=\sum_{n=N_{1}}^{N_{2}}{\phi_{m}[k_{j}-n]Q_{j,m;\theta}\left[n\right]},\end{split}$

(13)

it can be shown that, for any $\epsilon>0$, and for any selection $\left(j,m;\theta\right)$, there exist $N_{1}\left(\epsilon,k_{j},m;\theta\right)$ and $N_{2}\left(\epsilon,k_{j},m;\theta\right)$ such that $|{\hat{c}}_{j,m}-c_{j,m}|<\epsilon$. Therefore, with the finite sum (13), we may approximate $c_{j,m}$ to any accuracy. Furthermore, setting an upper bound on the energy of $\varphi_{m}\left(t\right)$, $N_{1}$ and $N_{2}$ may be chosen off-line, subject to the decay properties of the sequence $\left\{Q_{j,m;\theta}\left[n\right]\right\}_{n=-\infty}^{\infty}$.

With this result we can, for each $\theta$ and receiver index $m$, find the minimal set of $\varphi_{m}\left(t\right)$’s Fourier series coefficients, from which we can simultaneously approximate the set $\left\{c_{j,m}\right\}_{j=1}^{K}$. Denoting the set of indices corresponding to these coefficients by $\kappa_{m}$, arrange the set $\left\{\phi_{m}\left[k\right]\right\}_{k\in\kappa_{m}}$ in a length $K_{m}$ vector ${\bf{\Phi_{m}}}$. Using (13), we may define a linear transformation ${\bf{A_{m}}}\left(\theta\right)$, such that:

$\displaystyle\begin{split}{\bf{{\hat{c}}_{m}}}={\bf{A_{m}}}\left(\theta\right){\bf{\Phi_{m}}},\end{split}$

(14)

where ${\bf{\hat{c}_{m}}}$ is the length $K$ vector containing $\hat{c}_{j,m}$ as its $j$’th element. Instead of extracting only $K$ samples from each element using the complicated analog kernels, we now extract a larger number of samples, $K_{m}$. As we show in our simulation, good approximation is obtained with just a small sampling overhead. Our final Xampling scheme is illustrated in Fig. 2. Note the simple mechanism used for obtaining the Fourier coefficients in each individual element: a linear transformation, ${\bf{V_{m}}}$, is applied to point-wise samples of the signal, taken at a sub-Nyquist rate, after filtering it with an appropriate kernel, $s_{m}^{*}\left(-t\right)$, such as the Sum of Sincs  [4].

We examine the result of applying our scheme of Fig. 2 to raw RF data, acquired and stored for cardiac images of a healthy consenting volunteer. The acquisition was performed using a breadboard ultrasonic scanner of 64 acquisition channels. The transducer employed was a 64-element phased array probe, with $2.5\mbox{MHz}$ central frequency, operating in second harmonic imaging mode: 3 half cycle pulses are transmitted at $1.7\mbox{MHz}$, resulting in a signal characterized by a rather narrow bandpass bandwidth, centered at $1.7\mbox{MHz}$. The corresponding second harmonic signal, centered at $3.4\mbox{MHz}$, is then acquired. The signal detected in each acquisition channel is amplified and digitized at a sampling-rate of $50\mbox{MHz}$. Data from all acquisition channels were acquired along 120 beams, forming a $60^{\circ}$ sector, where imaging to a depth of $z=16\mbox{cm}$, we have $T=207{\mu\mbox{sec}}$. The results are illustrated in Fig. 3.

The first image (a) was generated using the standard technique, applying beamforming to data first sampled at the Nyquist rate, and then down-sampled, exploiting its limited essential bandwidth. For a single image line, sampling at $50{\mbox{MHz}}$, we acquire $10389$ real-valued samples from each element, which are then down-sampled, to $1662$ real-valued samples, used for beamforming. The resulting image is used as reference, which we aim at reproducing with our Xampling scheme. We begin by applying our scheme without approximation, simulating the modulation with the analog kernels defined in (11). Assuming $L=25$ reflectors, and using two-fold oversampling, $\kappa$ comprises $K=100$ indices. Since each sample is complex, we get an eight-fold reduction in sample-rate. The resulting image (b) well depicts strong perturbations observed in (a). Isolated reflectors at the proximity of the array ($z\approx 6{\mbox{cm}}$) remain well in focus. Applying our approximated scheme, for every $k_{j}\in\kappa$, $1\leq m\leq M$ and $\theta$, we choose $N_{1}$ and $N_{2}$ of (13) such that $\sum_{n=N_{1}}^{N_{2}}{\left|Q_{j,m;\theta}\left(\frac{2\pi}{T}n\right)\right|^{2}}\approx 0.95\sum_{n=-\infty}^{\infty}{\left|Q_{j,m;\theta}\left(\frac{2\pi}{T}n\right)\right|^{2}}$. At the current stage, this choice, and consequently the choice of $\left\{\kappa_{m}\right\}_{m=1}^{M}$ and $\left\{{\bf{A_{m}}}\right\}_{m=1}^{M}$, are performed off-line, based on the imaging setup. This results in an average number of $116$ complex samples per receiving element, required for constructing each image line (with at most $133$ samples required for some combinations of $\theta$ and element indices). We arrive at a seven-fold reduction in sample rate, and the resulting image (c) appears very similar to (b).