Computational Complexity Reduction of Synthetic-aperture Focus in Ultrasound Imaging Using Frequency-domain Reconstruction

Abstract

A new frequency-domain implementation of a synthetic aperture focusing technique is presented in the paper. The concept is based on synthetic aperture radar (SAR) and sonar that is a developed version of the convolution model in the frequency domain. Compared with conventional line-by-line imaging, synthetic aperture imaging has a better resolution and contrast at the cost of more computational load. To overcome this problem, point-by-point reconstruction methods have been replaced by block-processing algorithms in radar and sonar; however, these techniques are relatively unknown in medical imaging. In this paper, we extended one of these methods called wavenumber to medical ultrasound imaging using a simple model of synthetic aperture focus. The model, derived here for monostatic mode, can be generalized to multistatic as well. The method consists of 4 steps: a 2D fast Fourier transform of the data, frequency shift of the data to baseband, interpolation to convert polar coordinates to rectangular ones, and returning the data to the spatial-domain using a 2D inverse Fourier transform. We have also used chirp pulse excitation followed by matched filtering and spotlighting algorithm to compensate the effect of differences in parameters between radar and medical imaging. Computational complexities of the two methods, wavenumber and delay-and-sum (DAS), have been calculated. Field II simulated point data have been used to evaluate the results in terms of resolution and contrast. Evaluations with simulated data show that for typical phantoms, reconstruction by the wavenumber algorithm is almost 20 times faster than classical DAS while retaining the resolution.

Keywords

fast Fourier transform frequency domain synthetic aperture focus ultrasound imaging wavenumber

Introduction

Synthetic aperture imaging is based on the idea that a sequence of pulses recorded from a single moving real aperture can be treated as the output of a much longer array. The recorded data can be dynamically focused to reconstruct an image. The concept was first introduced by C. A. Wiley in 1951. He analyzed along-track spatial modulation (Doppler) in radar.^1,2 However, L. J. Cutrona and C. W. Sherwin studied the same ideas from different points of view.^3,4 First digital processor’s invention also heralded the introduction of modern synthetic aperture radar (SAR) systems.^5,6 In 1975, Cutrona a radar specialist, pointed out how the various aspects of SAR could be translated to synthetic aperture sonar (SAS) systems.^7,8 Research was simultaneously started in medical and nondestructive testing (NDT) with strong correspondence with SAS.^9,10

Synthetic aperture imaging provides an image with enhanced lateral resolution, axial resolution, and a higher contrast ratio than standard line-by-line imaging. The drawbacks are motion effects and more computational load. There are a lot of motion compensation algorithms to solve the first problem; however, the second one remains unresolved. Usually, the medical synthetic aperture implementations are performed using a delay-and-sum (DAS) processing in the time domain. The straightforward DAS algorithm, though physically understandable, has been shown to be quite time-consuming on general purpose computers due to a large number of operations. Hence, unlike medical applications, most synthetic aperture implementations in the framework of SAR and SAS are based on a developed version of the convolution model in the frequency domain. The basic concept of these techniques is based on the following model.

Image reconstruction is an inverse problem to create an image of the medium reflectivity from measurements of the echo signals recorded by set of apertures. Assuming that superposition applies, most algorithms try to invert the system model. Correlation or time-domain correlation algorithm is the simplest one. The algorithm is based on the correlation between echo signals and records the peak value for each image pixel. In broadband systems like medical ultrasound imaging, it is mathematically identical to DAS beamforming and backprojection. Time-domain algorithms are so compatible with different field parameters, array geometries, and arbitrary platform histories. The drawback is the computation time required to reconstruct a typical scene because of point-by-point processing.

In practice, the computational load makes it impossible to use synthetic aperture for real-time ultrasound imaging. To overcome this problem, point-by-point reconstruction methods have been replaced by block-processing algorithms in radar and sonar for many years. Although these techniques have been widely used in radar and sonar, they are relatively unknown in medical imaging. Conventional medical ultrasound imaging is performed using line-by-line transmission (real aperture) with classical DAS algorithm to reconstruct the image. Although DAS algorithm is so compatible with different parameters and platform geometry, it is time-consuming. However, synthetic aperture has essentially more computational load than real aperture imaging. Hence, the point-by-point DAS method is not an appropriate choice to create a real-time synthetic aperture imaging system. Due to the matrix processing which is simultaneously done over the whole data matrix in frequency-domain algorithms, the block-processing term is used in contrast with the point-by-point process in the DAS method.

Reviewing the literature, it is obvious that recent developments of frequency-domain methods in SAR and SAS based on efficient fast Fourier transform (FFT), for example,^11-17 have had no significant impact on the implementations of synthetic aperture imaging in medical imaging. Although some publications have presented different implementations of frequency-domain synthetic aperture, for example,^18-21 it seems that the recently developed algorithms using the framework for SAR and SAS are much less popular than the classical time-domain DAS schemes.^22-24 One possible explanation might be that the specific radar terminology used by radar specialists in SAR and SAS obscured this technology to specialists in the medical field. Another reason may be the parameter differences among those fields that make it hard to extend the methods or affect the result. The most effective parameters are depth of imaging, carrier frequency, bandwidth and beamwidth. Side effects of these differences include increasing side lobe levels, grating lobes, or contrast reduction.

One of the above-mentioned block-processing algorithms is wavenumber (also called Stolt mapping, Omega-K, or range migration algorithm).^16,25,26 In the present paper, we extended the wavenumber algorithm to the medical field and tried to compensate the consequent side effects. Hence, we used chirp excitation followed by matched filtering to compensate the effect of different frequency bands between radar and medical field that reduces the contrast and resolution. Furthermore, parameter differences increase the approximation error in the image reconstruction step that leads to much higher side lobes. To reduce this effect, we used spotlighting method, which is a method for aperture upsampling in radar. Finally, computational complexities of reconstruction using DAS and the proposed method have been calculated and compared. Field II simulated point data have been used to evaluate the results in terms of resolution, side lobes, and reconstruction time. Evaluations with simulated data show that for typical phantoms, reconstruction by wavenumber algorithm is almost 20 times faster than classical DAS while retaining the resolution.

The paper is organized as follows: in the next section, we present the models of the imaging system in the time domain and its Fourier transform version. This presentation is followed by the proposed algorithm including chirp pulse excitation and matched filtering, wavenumber algorithm, and side lobes and grating lobes suppression. In the “Proposed Method” section, computational complexities of the reconstruction algorithms, DAS, and wavenumber are calculated. The section “Computational Complexity” presents the results of simulations performed for the DAS and the proposed algorithm and a comparison between their performances. Moreover, we show experimental results of processing times for three phantoms with different sizes and compare them with their theoretical value. In the section “Simulations and Results,” we have a discussion about the proposed method, advantages, and drawbacks.

Theory

System Modeling

Assume a stationary target region composed of a set of point reflectors with reflectivity $σ_{n}$ located at coordinates $(x_{n}, y_{n}), n = 1, 2, \dots$ in the spatial $(x, y)$ domain.

The variable $x$ is used for the range domain and $y$ identifies the cross-range (also called azimuth or along-track). A transducer element located at $(0, u)$ in the spatial domain illuminates the target area with a multifrequency (large-bandwidth) signal $p (t)$ The measured echo signal is

s (t, u) = \sum_{n} σ_{n} p [t - \frac{2 \sqrt{x_{n}^{2} + {(y_{n} - u)}^{2}}}{c}],

where $c$ is the wave propagation speed and $t$ represents the time. The Fourier transform of $s (t, u)$ is

s (ω, u) = P (ω) \sum_{n} σ_{n} \underset{Spherical PM signal}{\underset{︸}{\exp [- j 2 k \sqrt{x_{n}^{2} + {(y_{n} - u)}^{2}}]}},

where $k = ω / c$ is the wavenumber. As one can see in the above model, the signal in the $(ω, u)$ domain is composed of a linear combination of the spherical pulse modulated (PM) signals.

For the generic model, we assume that the measurements are made for $u \in (- \infty, \infty)$ The Fourier transform of the spherical PM signal with respect to slow-time $u$ is

ℱ_{(u)} [\exp [- j 2 k \sqrt{x_{n}^{2} + {(y_{n} - u)}^{2}}]] = \exp (- j \sqrt{4 k^{2} - k_{u}^{2}} x_{n} - j k_{u} y_{n}),

for $k_{u} \in [- 2 k, 2 k]$ where $k_{u}$ is called spatial or Doppler frequency. Using the spherical PM signal property, the Fourier transform of $s (ω, u)$ is

S (ω, k_{u}) = P (ω) \sum_{n} σ_{n} \underset{Linear phase function of (x_{n}, y_{n})}{\underset{︸}{\exp (- j \sqrt{4 k^{2} - k_{u}^{2}} x_{n} - j k_{u} y_{n}) .}}

We can rewrite the equation with new functions as

S (ω, k_{u}) = P (ω) \sum_{n} σ_{n} \exp [- j k_{x} (ω, k_{u}) x_{n} - j k_{y} (ω, k_{u}) y_{n}],

where the two new functions are defined to be

\begin{matrix} k_{x} (ω, k_{u}) = \sqrt{4 k^{2} - k_{u}^{2}}, \\ k_{y} (ω, k_{u}) = k_{u} . \end{matrix}

These two functions are called the spatial frequency mapping or transformation.

We define the ideal target function in the spatial domain via

f_{0} (x, y) = \sum_{n} σ_{n} δ (x - x_{n}, y - y_{n}),

which has the following two-dimensional (2D) spatial Fourier transform:

F_{0} (k_{x}, \underset{Linear combination}{\underset{︸}{k_{y}) = \underset{Linear phase function}{\underset{︸}{\sum_{n} σ_{n} \exp (- j k_{x} x_{n} - j k_{y} x_{n})}}}} .

Note that $F_{0} (k_{x}, k_{y})$ is also composed of a linear combination of linear phase functions. This is due to the Fourier shift property.

Using Equation (8) in Equation (5), we have

S (ω, k_{u}) = P (ω) F_{0} [k_{x} (ω, k_{u}), k_{y} (ω, k_{u})] .

Then the ideal target function is

F_{0} [k_{x} (ω, k_{u}), k_{y} (ω, k_{u})] = \frac{S (ω, k_{u})}{P (ω)} .

As $p (t)$ is a bandlimited signal, the above equation is a theoretical reconstruction. In practical reconstruction, we can estimate the ideal function via matched filtering, that is,

\begin{matrix} F [k_{x} (ω, k_{u}), k_{y} (ω, k_{u})] = P^{*} (ω) S (ω, k_{u}) \\ = {| P (ω) |}^{2} \sum_{n} σ_{n} \exp (- j k_{x} x_{n} - j k_{y} y_{n}), \end{matrix}

for $k_{u} \in [- 2 k, 2 k]$ and $ω \in [ω_{c} - ω_{0}, ω_{c} + ω_{0}]$ .

Proposed Method

We divided ultrasound image reconstruction into three steps: pulse transmission, reconstruction, and error reduction. In the proposed method, we try to select the best signal to achieve the maximum information and as a result a better resolution and contrast. Thus, different signals and their effects in the reconstructed image are introduced in the first part. As mentioned before, computational load is a problem in real-time synthetic aperture ultrasound imaging. To overcome this problem, a block-processing algorithm called wavenumber is developed for image reconstruction in the second part. Finally, a spotlighting method is introduced to reduce the side lobes and grating lobes caused by differences in parameters between radar and medical imaging.

Chirp Pulse Excitation and Matched Filtering

Chirp pulses make up a special class of signals that can possess long duration and wide bandwidth simultaneously. The long duration provides more energy and better contrast. Also, the signal bandwidth determines the resolvability of the targets in the range domain or axial resolution:

Δ_{x} = \frac{π}{K_{0}} = \frac{π c}{2 ω_{0}},

where $[- ω_{0}, ω_{0}]$ is the baseband signal bandwidth in Hertz.

Assume a rectangular pulse with duration $τ_{0}$ as the baseband signal. Then the spectral shape of the signal is a sinc function with first zero in $ω_{0}$ :

ω_{0} = \frac{2 π}{τ_{0}},

which is a half bandwidth. Then the axial resolution is

Δ_{x} = \frac{c τ_{0}}{4} .

It is obvious that we cannot increase the two variables, duration and bandwidth, simultaneously due to their inverse relationship.

Chirp pulse is a linear frequency modulated signal that has been widely used in radar.²⁷ The most important property of the signal is the capability to increase the pulse duration and bandwidth simultaneously. Chirp pulse transmission also used in time-domain ultrasound imaging to improve the signal-to-noise ratio and contrast especially in higher depth.^19,22,23

A chirp signal is defined as

p (t) = {\begin{matrix} \exp [j (ω_{c} t + α t^{2})] & t \in [0, τ_{0}] \\ 0 & otherwise \end{matrix},

$α$ is called the chirp rate. The instantaneous frequency of the chirp pulse is the derivative of its phase function with respect to time:

ω_{i} (t) = \frac{d}{d t} [ω_{c} t + α t^{2}] = ω_{c} + 2 α t .

Note that with $α > 0$ , the instantaneous frequency is an increasing function of time and is bounded by

ω_{c} < ω_{i} (t) < ω_{c} + 2 α τ_{0} .

The spectral bandwidth of the signal in the frequency domain is $[\frac{2 π}{τ_{0}}, \frac{2 π}{τ_{0}} + 2 α τ_{0}]$ . Thus, its baseband bandwidth is

2 ω_{0} = 2 α τ_{0} .

Clearly, the baseband bandwidth of a chirp pulse increases with its duration $τ_{0}$ .

The axial resolution is then

Δ_{x} = \frac{π c}{2 ω_{0}} = \frac{π c}{2 α τ_{0}} .

Note that $ω_{0}$ has a direct relationship with $α$ . Hence, we can use it to adjust bandwidth and axial resolution. A wide-bandwidth signal has the low and high frequencies simultaneously that makes it possible to have a good lateral resolution in different depths.

If $p (t)$ is not a sharp signal (similar to an impulse function) and has relatively large duration (e.g., chirp signal), then the hyperbolic signatures of the targets may overlap with each other. One way to improve the visualization is matched filtering, that is,

s_{M} (t, u) = s (t, u) \times p^{*} (- t) = \sum_{n} σ_{n} {psf}_{t} (t - \frac{2 R_{n}}{c}),

where the point spread function, ${psf}_{t}$ is

{psf}_{t} (t) = ℱ_{ω}^{- 1} [{| P (ω) |}^{2}],

where $ω$ presents the frequency domain.

As mentioned above, we can use the long-duration chirp signal to have a high transmitting energy (better signal-to-noise ratio and contrast) and wide bandwidth (better resolution) simultaneously. It is not the case for rectangular signal (i.e., a chirp signal with zero chirp rate, $α = 0$ ).

Wavenumber Algorithm

Among the frequency-domain SAR algorithms, wavenumber is the most popular one. The algorithm was first introduced in geophysics named seismic migration, whereas it is called Omega-K or range migration algorithm in SAR; however, we will refer to all these collectively as the wavenumber algorithm. Some publications in ultrasound imaging also have presented different implementations of the frequency-domain synthetic aperture.^18-21

The algorithm starts with either the raw or matched-filtered echo data and immediately performs a 2D FFT to calculate a 2D spectrum of the ultrasonic data. After compensating the phase delay to shift the data to the baseband, the spectrum is interpolated to convert the polar coordinate system to the rectangular one. The transformed spectrum is then returned to the time domain by a 2D inverse Fourier to achieve a final image. Given that the transfer function of the data collection process is as shown in Equation (9), the matched filter wavenumber image reconstruction algorithm can be summarized by

F (k_{x}, k_{y}) = S {P^{*} (ω) S (ω, k_{u})},

where the coordinate reformatting (the Stolt mapping) is given by Equation (6). Due to the nonlinear nature of the 2D Stolt mapping, the resultant database of $F (k_{x}, k_{y})$ is unevenly spaced. Moreover, we require the knowledge of $F (k_{x}, k_{y})$ on a (uniform) rectangular grid to retrieve a sampled version of $f (x, y)$ via 2D FFTs. Thus, we need an interpolation step to achieve the measurements in new, evenly spaced coordinates.

Assume that the illuminated target area in the depth domain is defined via the region $x \in [X_{c} - X_{0}, X_{c} + X_{0}]$ , where $X_{c}$ is its mean range and $2 X_{0}$ is the size of imaging area in the depth domain. Thus, the signal $F (k_{x}, k_{y})$ is a bandpass signal in the $k_{x}$ domain, as its inverse Fourier transform is centered at $x = X_{c}$ The inverse Fourier transform of a lowpass signal in the $k_{x}$ domain has its spectrum centered at $x = 0$ For the interpolation step, we need to transform this signal to the lowpass one. For this purpose, we perform the following baseband conversion:

F_{b} (k_{x}, k_{y}) = F (k_{x}, k_{y}) \exp (j k_{x} X_{c}) .

Note that the origin in the spatial $(x, y)$ domain is at the center of the baseband target function $f_{b} (x, y)$ The reconstruction equation for the lowpass target function becomes

F_{b} (k_{x}, k_{y}) = S {P^{*} (ω) \exp (j \sqrt{4 k^{2} - k_{u}^{2}} X_{c}) S (ω, k_{u})} .

If we assume $P (ω) \exp (- j \sqrt{4 k^{2} - k_{u}^{2}} X_{c})$ as a response of a point reflector in the middle of imaging area, then the above baseband version of image reconstruction has some similarities to matched filtering in the cross-range domain.

A block diagram of the method is shown in Figure 1.

Side Lobes and Grating Lobes Suppression

As it was shown in Equation (11), spatial frequency is bounded by $k_{u} \in [- 2 k, 2 k]$ for synthetic aperture imaging. Sampling constraint in the spatial domain based on the Nyquist theorem is

Δ_{u} \leq \frac{π}{2 k} = \frac{λ}{4} .

In the case of a larger sampling distance, grating lobes and side lobes smear into a main lobe target. It is obvious that synthetic aperture imaging has a twice wider beamwidth in comparison with the classical real aperture imaging, where $k_{u} \in [- k, k]$ . Thus, sampling constraint is much restricted for synthetic aperture imaging than real aperture. Also, the transmission of a wide-bandwidth signal causes frequency-dependent grating lobe targets.

As mentioned before, we selected a wide-bandwidth chirp signal in transmission that has a wide frequency spectrum. In Equation (26), carrier frequency is used to determine the sampling distance. Note that for higher frequencies, this constraint is not enough and $λ$ must be replaced by $λ_{m i n}$ in Equation (26). In practice, spatial sampling distance is determined by transducer pitch that is limited by manufacturing technology of making trancducer.

Windowing is a common way to control the side lobes that usually accompany with loss in resolution. Adaptive beamformer is a nonlinear weighting method introduced for time-domain image reconstruction that makes it possible to reduce the side lobes level without loss in resolution.

There are very few researches that have investigated a method to control the side lobes and grating lobes in wide-band frequency-domain applications. Some linear and nonlinear window functions are introduced in Hawkins²⁸ and Vu.²⁹ Note that we want a block-processing method; however, some of investigated filters are pixel based that is in contrast with our goal to reduce the computational complexity.³⁰ Moreover, applying 2D weighting functions to the spectrum of SAR images usually accompanies with loss in resolution.

In this paper, we used spotlighting to control the undesired lobes. This is a frequency-domain upsampling method in radar,¹⁶ where, spatial sampling constraint is not usually fulfilled. Then, spotlighting is an algorithm to produce alias-free data from the aliased one. We used the method to control the side lobes and grating lobes.

For the small extended target on boresight, the linear spatial frequency dependence of $u$ is approximately known and may be removed to produce a spatially compressed signal $s_{c} (ω, u)$ given by

s_{c} (ω, u) = s (ω, u) \exp (j 2 k \sqrt{X_{c}^{2} + {(Y_{c} - u)}^{2}} - j 2 k R_{c}),

where $(X_{c}, Y_{c})$ is the center of the imaging area and $R_{c}$ is its distance to the center of transducer. Assume a wide-band excitation signal with $f_{m a x} = 2 f_{c}$ , where $f_{c}$ and $f_{m a x}$ are carrier frequency and maximum frequency in the bandwidth, respectively. A standard medical transducer’s element pitch is usually $λ_{c} / 2$ that is also assumed as spatial sampling frequency. According to above assumptions, we should extend the frequency band four times wider to eliminate aliasing (based on Nyquist theorem). A new temporal frequency domain can be defined as

S_{c d} (ω, k_{u}) = {\begin{matrix} 0 & - \frac{8 π}{λ_{c}} \leq k_{u} \leq - \frac{2 π}{λ_{c}} \\ S_{c} (ω, k_{u}) & f o r - \frac{2 π}{λ_{c}} \leq k_{u} \leq \frac{2 π}{λ_{c}} \\ 0 & \frac{2 π}{λ_{c}} \leq k_{u} \leq \frac{8 π}{λ_{c}} \end{matrix},

where $λ_{m i n} / 4 = λ_{c} / 8$ is the sampling rate that should be applied.

Accordingly, we can calculate the unaliased signal $s_{d} (ω, u)$ from the compressed one $S_{c d} (ω, k_{u})$ :

s_{d} (ω, u) = s_{c d} (ω, u) \exp (- j 2 k \sqrt{X_{c}^{2} + {(Y_{c} - u)}^{2}} + j 2 k R_{c}) .

Note that $s_{d} (ω, u)$ is a spatially upsampled version of $s (ω, u)$ without the effects of the original undersampling. Now, we can use $s_{d} (ω, u)$ instead of the original aliased signal $s (ω, u)$ as if it had been adequately sampled, that is, if the sampling rate was indeed $λ_{c} / 8$ .

Computational Complexity

In this section, computational complexities of the classical DAS and our proposed method are calculated and compared. Note that we did not consider matched filtering and spotlighting algorithms in the following discussion for any of the methods. As we use a monostatic configuration, the delays should be calculated once for each pixel for each echo received signal. The delay from the pixel to each transmitter/receiver is different and also it is particularly different for individual pixels.

Assuming an N elements transducer with M samples in each line of data, there are $M \times N$ samples in the data matrix. For simplicity, we assume that the number of pixels in the final image is equal to the data samples. To be realistic, an optimized sampling frequency can be assumed for two methods.

Each pixel in the DAS method is formed using one delay and a one-dimensional (1D) interpolation. Number of calculations for the whole image is then,

M N^{2} [Distance + 1 D interpolation] .

If we use Euler distance in a 2D space, that is,

Distance = \sqrt{{(x_{1} - x_{2})}^{2} + {(y_{1} - y_{2})}^{2}},

then we need three additions/subtractions, two multiplications, and one square root to calculate the distance. The 1D interpolation is also given by

\begin{matrix} \hat{y} (n + η) = (1 - η) y (n) + η y (n + 1) \\ = y (n) + η [y (n + 1) - y (n)], \end{matrix}

which includes two additions/subtractions and one multiplication. Assuming $n$ digit numbers, it needs $n$ calculations for additions/subtractions and $n^{2}$ calculations for multiplications. (Note that there are several algorithms for multiplication with different complexities. The simplest one is used in this paper.) In the following discussion, square root is put in the multiplication group due to the same calculation complexity. Moreover, the shortened form “ $m u l t$ ” for multiplication, “ $s u b$ ” for addition/subtraction, and “ $c o m p l e x_m u l t$ ” for complex numbers multiplication are used. Assuming 8-digit numbers, the calculation for DAS is as follows:

M N^{2} [4 m u l t + 5 s u b] = M N^{2} [4 \times n^{2} + 5 \times n] = M N^{2} [4 \times 8^{2} + 5 \times 8] = M N^{2} \times 296 .

However, our proposed method that is a block-processing algorithm needs a FFT transform for whole data matrix, phase correction, and an interpolation for each pixel and an inverse fast Fourier transform (IFFT). The calculations are as follows:

FFT + M N [phase correction + 1 D complex interpolation] + IFFT .

FFT and IFFT need $M N \log_{2} (M N) \times n$ computations. The $n$ is 8 for FFT and 16 for IFFT because of the complex number in the wavenumber domain. Phase correction is given by

\hat{S} (ω, k_{u}) = \exp (j \sqrt{4 k^{2} - k_{u}^{2}} X_{c}) S (ω, k_{u}),

which needs one addition/subtraction, four multiplications, a square root, and an exponential. The number of computations is $n^{2} \log_{2} n$ for exponential function where $n$ is 8 due to having real numbers in this step. However, it needs a complex multiplication for final phase correction. Thus, the phase correction computation number is given by

\begin{matrix} M N [5 m u l t + 1 s u b + n^{2} \log_{2} n + c o m p l e x_m u l t] = M N [5 n^{2} + n + n^{2} \log_{2} n + {(2 n)}^{2}] \\ = M N [9 n^{2} + n + n^{2} \log_{2} n] . \end{matrix}

Complex interpolation also needs three additions/subtractions and four multiplications. Then the overall computations are

\begin{array}{l} M N \log_{2} (M N) \times 8 + M N [13 n^{2} + 4 n + n^{2} \log_{2} n] + M N \log_{2} (M N) \times 16 \\ = M N [13 \times 8^{2} + 4 \times 8 + 8^{2} \times \log_{2} 8 + 24 \log_{2} (M N)] \\ = M N [1056 + 24 \log_{2} (M N)] . \end{array}

The relative computation number of DAS to the proposed method for 8-digit real numbers (16-digit complex numbers) is

296 \times \frac{N}{[1056 + 24 \log_{2} (M N)]} .

Note that the calculations are just for wavenumber and DAS algorithms, and any other pre- or post-processing algorithm’s computational load must be calculated, separately.

Simulations and Results

The simulations were performed for the synthetic aperture created from the measurements made by a linear transducer. A 96-element transducer having $λ_{c} / 2$ pitch was assumed to be excited by a 4 MHz carrier frequency and the echoes were recorded in the 100 MHz sampling frequency ( $λ_{c}$ is the wavelength in the carrier frequency). As a monostatic (or synthetic aperture focus [SAF]) configuration is used, one element is active in each transmit and receive. The transducer was excited using two cycles of sinusoid and in the case of chirp excitation, the bandwidth was assumed to be twice of carrier frequency. We used matched filtering as mentioned in section “Chirp Pulse Excitation and Matched Filtering” for either of sinusoidal or chirp pulses. The data were simulated using ultrasound simulation software Field II^31,32 and the DAS beamformation has done using beamformation toolbox III (BFT3).³³

As mentioned in the “Chirp Pulse Excitation and Matched Filtering” section, chirp pulse excitation provides a possibility to have a long-duration pulse as well as having a wide bandwidth, which corresponds to more transmitted energy (better contrast) and better resolution, respectively. Chirp rate variable $α$ should be set to 0 for sinusoidal and more for chirp signal. Two simulated data with chirp rates 0 and 4 are shown in Figure 2(a) and (b). As shown, exciting a pulse with a higher chirp rate results in an image with a better resolution and lower side lobes level. Note that a much higher chirp rate may result in the undesired effects.

In section “Side Lobes and Grating Lobes Suppression” we introduced a method called spotlighting to suppress the side lobes and grating lobes. Due to the wide bandwidth of the chirp pulse, the spatial sampling frequency selected from the carrier frequency is not efficient for higher frequencies and causes the grating lobes to be formed. The shape and location of the grating lobes are dependent on the target’s distances from the center of imaging region and signal bandwidth.

Figure 3(a) and (b) shows the structure of grating lobes in a typical phantom before and after the spotlighting algorithm. For better presentation of the method’s efficiency, a 400-element transducer was used to enlarge the field of view. It makes the grating lobes come into the region of interest (ROI).

As it can be seen, the algorithm is able to suppress the grating lobes. Although, reducing the ROI width results in the lower grating lobes in the field of view, spotlighting is still useful to reduce the side lobes level.

Results of reconstruction algorithms, DAS and wavenumber, are compared in Figure 4(a) and (b) on a point target phantom. The excitation pulse was two cycles of sinusoid and a Hanning apodization window is used for both of them. The axial and lateral beam patterns are also compared in Figure 4(c) and (d). For quantitative evaluation, the values of full width at half maximum are compared in Table 1. As it can be seen, values are the same in the case of the lateral pattern and slightly better for the axial resolution of the wavenumber algorithm in comparison with DAS.

Table 1.

Full Width at Half Maximum Values of Beam Patterns in Figure 4(c) and (d).

	Full Width at Half Maximum (mm)
Beam Patterns	Wavenumber	DAS
Axial beam pattern	0.3162	0.3570
Lateral beam pattern	0.8278	0.8278

DAS = delay-and-sum.

The algorithms are also simulated on two different phantoms to study their effect on off-axis boresight targets in Figure 5(a) and (b) and different depth in Figure 6(a) and (b). The lateral beam patterns are also compared in Figure 5(c) and Figure 6(c) to (e). The results of first phantom (Figure 5) show a lower side lobe level in boresight targets for the wavenumber algorithm than DAS. Compared with Figure 4 (where the side lobe levels of DAS are lower than the wavenumber), here the point targets are on the off-axis (lateral distance of +5 mm, −5 mm) resulting in the lower side lobe levels for the wavenumber. The explanation could be using different implementation and also phase information in the wavenumber algorithm that makes it more efficient in this condition. This is also true for the targets in the higher depth (Figure 6).

To study more realistic conditions, we also simulated a cyst phantom of size 10 (mm) × 10 (mm) × 10 (mm) centering at a 30-mm depth with a 3-mm radius. The results are shown in Figure 7(a) and (b). The values for signal-to-noise and contrast ratios are also shown in Table 2. As it can be seen, the results are better for the wavenumber algorithm than DAS.

Table 2.

Contrast and Signal-to-Noise Ratios for Figure 7(a) and (b).

	Wavenumber	DAS
Contrast ratio	0.9611	0.9494
Signal-to-noise	0.9800	0.9737

DAS = delay-and-sum.

As it was already mentioned, the motivation of selecting a wavenumber algorithm for image reconstruction was to reduce the computational load and reconstruction time. We tend to achieve a real-time, high-resolution synthetic aperture imaging by extending the method to multistatic modality in future works. The reconstruction times of DAS and the proposed method are compared for three different phantoms in Table 3. These phantoms numbered 1 to 3 have 80-, 40-, and 10-mm depths, respectively. The times were measured on a PC with Core i7 2.8 GHz CPU and 3 GB RAM.

Table 3.

Reconstruction Times Compared for DAS and the Proposed Method.

	Processing Time (s)
Phantom No. $(M \times N)$	DAS	Proposed Method	Time Ratio (DAS/Proposed Method)
No. 1 $(10, 391 \times 96)$	10.89	0.60	18.15
No. 2 $(5196 \times 96)$	5.41	0.28	19.32
No. 3 $(1300 \times 96)$	1.39	0.71	19.57

DAS = delay-and-sum.

The times were just measured for the reconstruction step, and any other pre- and post-processing algorithms such as matched filtering, windowing, or spotlighting time were not considered. As it is shown, the processing time is about 20 times better for the proposed method in comparison with the classical DAS method. Based on the theory developed in Equation (37), these time ratios (DAS/proposed method) should be 18.52, 18.81, and 19.43, respectively, that are so close to the experimental values shown in Table 3. As it was expected, the ratio decreases as the imaging depth increases. Note that the interpolation used in the above calculations is a linear interpolation.

Conclusion

A frequency-domain reconstruction algorithm called wavenumber was proposed in the paper. As previously mentioned, synthetic aperture imaging in its classical forms has a huge computational load that makes it impossible to use in real-time imaging. To overcome this problem, block-processing algorithms are used in radar and sonar instead of point-by-point DAS method. Although these techniques are popular in radar and sonar, they are almost unknown in the medical field.

Wavenumber algorithm was extended to the medical ultrasound imaging in this paper. At first, a 2D FFT was used to calculate a 2D spectrum of the ultrasonic data. After compensating the phase delay to shift the data to baseband, the spectrum was interpolated to convert the polar coordinate system to the rectangular one. The transformed spectrum was then returned to the time domain with a 2D inverse Fourier. Furthermore, we used the chirp pulse excitation followed by matched filtering to increase the frequency band and pulse duration simultaneously (corresponding to contrast and resolution enhancement). Matched filtering was also used to eliminate the distortion effect of long pulse excitation. Results showed that using chirp pulse excitation with matched filtering enhances the resolution and contrast. Moreover, the wavenumber algorithm reduces the computation time by about 20 times when compared with the classical DAS method retaining the resolution and contrast. It is also shown to perform better in some conditions like boreside and higher depth targets. In the case of a more complicated cyst phantom, the wavenumber algorithm showed a better performance in terms of signal-to-noise ratio and contrast. The explanation could be much more phase information that is used in the wavenumber algorithm in comparison with DAS that just uses an amplitude of delayed signals. Finally, we introduced a method to suppress grating lobes in the case of using wide-bandwidth chirp pulse excitation called spotlighting that showed an acceptable efficiency in simulated results. The computational load of the proposed method was compared with DAS and the values were so close to the theoretical ones.

The method could be extended to multistatic modality in the future work to achieve a real-time, high-resolution configuration for ultrasound imaging. Moreover, we can pursue the research on the nonlinear frequency-domain windowing to reduce the side lobes level. Other frequency-domain block-processing algorithms such as chirp scaling with less computational load can be a topic for future research, the problem of which is more approximations in equations that make it more difficult to extend to medical imaging; however, the approximations are acceptable for radar imaging.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by Tarbiat Modares University, Tehran, Iran.

References

Wiley

. Synthetic aperture radars. IEEE Trans Aerosp Electron Syst. 1985;21:440-3.

Wiley

. Pulsed Doppler radar methods and apparatus. U.S. Patent; 3196436; 1965.

Sherwin

Ruina

Rawcliffe

. Some early developments in synthetic aperture radar systems. IRE Trans Mil Electron. 1962;6:111-5.

Cutrona

Vivian

Lieth

Hall

. A high resolution radar combat-surveillance system. IRE Trans Mil Electron. 1961;5:127-31.

Curlander

McDonough

. Synthetic Aperture Radar: Systems and Signal Processing. New York: Wiley; 1991.

Kirk

. A discussion of digital processing in synthetic aperture radar. IEEE Trans Aerosp Electron Syst. 1975;AES-11:338-48.

Cutrona

. Comparison of sonar system performance achievable using synthetic-aperture techniques with the performance achievable with more conventional means. J Acoust Soc Am. 1975;58:336-48.

Cutrona

. Additional characteristics of synthetic-aperture sonar systems and a further comparison with nonsynthetic-aperture sonar systems. J Acoust Soc Am. 1977;61:1213-7.

Burckhardt

Grandchamp

Hoffman

. An experimental 2MHz synthetic aperture sonar system intended for medical use. IEEE Trans Sonics Ultrason. 1974;SU-21:1-6.

10.

Thomson

. Transverse and longitudinal resolution of the synthetic aperture focusing technique. Ultrasonics. 1984;22:9-15.

11.

Bamler

. A comparison of range-Doppler and wave number domain SAR focusing algorithms. IEEE Trans Geosci Remote Sensing. 1992;30:706-13.

12.

Carrara

Goodman

Majewski

. Spotlight Synthetic Aperture Radar: Signal Processing Algorithms. Boston, MA: Artech House; 1995.

13.

Cumming

Wong

Raney

. A SAR processing algorithm with no interpolation. Int Geosci Remote Sensing Symp. 1992;1:376-9.

14.

Raney

Runge

Bamler

Cumming

Wong

. Precision SAR processing using chirp scaling. IEEE Trans Geosci Remote Sensing. 1994;32:786-99.

15.

Runge

Bamler

. A novel high precision SAR focusing algorithm based on chirp scaling. Int Geosci Remote Sensing Symp. 1992;1:372-5.

16.

Soumekh

. Synthetic Aperture Radar Signal Processing with MATLAB Algorithms. New York: Wiley; 1999.

17.

Zaugg

Long

. Generalized frequency-domain SAR processing. IEEE Trans Geosci Remote Sensing. 2009;47:3761-73.

18.

Busse

. Three-dimensional imaging using a frequency domain synthetic aperture focusing technique. IEEE Trans Ultrason Ferroelec Freq Contr. 1992;39:174-9.

19.

Gran

Jensen

. Frequency division transmission imaging and synthetic aperture reconstruction. IEEE Trans Ultrason Ferroelec Freq Contr. 2006;53:900-11.

20.

J-Y

Cheng

Wang

. High frame rate imaging system for limited diffraction array beam imaging with square-wave aperture weightings. IEEE Trans Ultrason Ferroelec Freq Contr. 2006;53:1796-812.

21.

Garcia

Tarnec

Muth

Montagnon

Porée

Cloutier

. Stolt’s f-k migration for plane wave ultrasound imaging. IEEE Trans Ultrason Ferroelec Freq Contr. 2013;60:1853-67.

22.

Nikolov

. Synthetic aperture tissue and flow ultrasound imaging. PhD dissertation, Technical University of Denmark, Kongens Lyngby, 2001.

23.

Gammelmark

. Improving the image quality of synthetic transmit aperture ultrasound images. PhD dissertation, Technical University of Denmark, Kongens Lyngby, 2004.

24.

Kortbek

. Synthetic aperture sequential beamforming and other beamforming techniques in ultrasound imaging. PhD dissertation, Technical University of Denmark, Kongens Lyngby, 2007.

25.

Dendal

Marchand

. Ω-K techniques advantages and weaker aspects. Proc IEEE Int Geosci Remote Sens Symp. 1992;1:366-8.

26.

Bamler

. A comparison of range-Doppler and wave number domain SAR focusing algorithms. IEEE Trans Geosci Remote Sensing. 1992;30:706-13.

27.

Soumekh

. Fourier Array Imaging. Englewood, CA: Prentice Hall; 1994.

28.

Hawkins

. Synthetic aperture imaging algorithms: with application to wide bandwidth sonar. PhD dissertation, University of Canterbury, Christchurch, New Zealand, 1996.

29.

. Ultrawideband-ultrawidebeam synthetic aperture radar signal processing and applications. PhD dissertation, School of Engineering Blekinge Institute of Technology, Karlskrona, Sweden, 2011.

30.

Moreira

. Suppressing the azimuth ambiguities in synthetic aperture radar images. IEEE Trans Geosci Remote Sensing. 1993;31:885-95.

31.

Jensen

Svendsen

. Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers. IEEE Trans Ultrason Ferroelec Freq Contr. 1992;39:262-7.

32.

Jensen

, Field: a program for simulating ultrasound systems. Med Biol Eng Comp, vol. 10th Nordic-Baltic Conference on Biomedical Imaging, vol. 4, Suppl. 1, Part 1, pp. 351-3., 1996.

33.

Hansen

Jensen

. An object-oriented multi-threaded soft-ware beamformation toolbox. In: SPIE Medical Imaging, Lake Buena Vista (Orlando), FL, 12-17 February 2011.