- Email: [email protected]
A method of improving overall resolution in ultrasonic array imaging using spatio-temporal deconvolution
Ultrasonics 42 (2004) 961–968 www.elsevier.com/locate/ultras
A method of improving overall resolution in ultrasonic array imaging using spatio-temporal deconvolution Fredrik Lingvall Department of Engineering Sciences, Signals and Systems Group, Uppsala University, Box 528, 751 20 Uppsala, Sweden
Abstract In this paper a beamforming method for ultrasonic array imaging is presented that performs both spatial and temporal deconvolution based on a minimum mean square error (MMSE) criteria. The presented MMSE receive mode beamformer performs a regularized inversion of the propagation operator for the ultrasonic array system at hand. The MMSE beamformer accounts for the transmit and receive processes, defined in terms of finite array element sizes, transmit focusing laws and electrical transducer characteristics. The MMSE beamformer is compared to the traditional delay-and-sum (DAS) beamformer with respect to both resolution and signal-to-noise ratio. The two algorithms are compared using both simulated and measured data. The simulated data was obtained using ultrasonic field simulations and the measured data was acquired using a linear phased array imaging wire targets in water. The results show that the MMSE beamformer has superior temporal and lateral resolution compared to DAS. It is also shown that the MMSE beamformer can be expressed as a filter bank, which enables parallel processing at high frame rates. 2003 Elsevier B.V. All rights reserved. Keywords: Ultrasonic array imaging; Spatio-temporal deconvolution; Beamforming
1. Introduction In conventional ultrasonic array imaging focusing in receptions is performed by means of delay-and-sum (DAS) operations on the signals received by the array elements. To successfully apply DAS receive mode beamforming the array elements must be small compared to the wavelength, otherwise, the diffraction effects will degrade the processed images [1]. There is, however, a trade-off between element size and signal-tonoise ratio (SNR) and the elements must therefore be large enough to generate a sufficient amount of acoustic energy to obtain a sufficiently high SNR. Another important issue that also limits the element size in DAS beamforming is that, if undersampled arrays with element spacing larger than half the wavelength are used, the images will be degraded by the grating lobes. Avoiding grating lobes by limiting the element spacing to less than half the wavelength also limits the element size and in consequence also the transmitted acoustic power. In order to improve the SNR it is, therefore,
E-mail address: fl@signal.uu.se (F. Lingvall). 0041-624X/$ - see front matter 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.ultras.2003.12.016
desirable to increase the size of the array elements. Even if advanced methods such as dynamic receive focusing and parallel beamforming [2,3], also based on DAS processing, are used to improve image quality and frame-rate the fundamental limitations of DAS approach can not be circumvented. Another deficiency of DAS beamforming is that the electrical impulse response of the transducer is not taken into consideration and the DAS processing will, therefore, not improve temporal resolution. Recently, methods that are based on inverse filtering of the propagation operator have been introduced to overcome the limitations of the DAS approach [4,5]. Since the propagation operator describes both the transmission and reception processes, such an approach could in theory provide full information of the scattering objects in the insonified region of interest. However, the inverse operator is generally ill-posed and regularized solutions must be sought, for instance sub-space methods [4,5]. A serious disadvantage of sub-space regularization methods is that they do not incorporate any prior knowledge, such as, the statistical distribution of the measurement noise and the fact that the target strength is limited.
962
F. Lingvall / Ultrasonics 42 (2004) 961–968
In this paper a recently introduced inverse filtering approach for synthetic aperture imaging [6], is generalized and applied to ultrasonic array imaging. This new approach, that uses a priori statistical information about the scattering targets and the measurement noise, uses a model based method and a minimum mean squared error (MMSE) criteria to find a solution that minimizes the reconstruction error for the array system at hand. The linear model takes both the acoustic transmission and reception processes into account, using the general concept of spatial impulse responses [7], as well as the electrical impulse response of the transducer. The application of MMSE method is, therefore, not limited to any particular constraints regarding the element size, element spacing, or transmit focusing method. The MMSE method can, for example, be applied to undersampled arrays and arrays with finite sized elements. Moreover, since the linear model includes the electrical impulse response the temporal resolution will also be improved. This paper is organized as follows: The linear discrete model and the MMSE filter is described in Section 2. The DAS and the MMSE methods are then compared, using both simulated and measured data, in Section 3. Finally, the conclusions are given in Section 4.
2. Theory 2.1. A discrete linear model for ultrasonic array imaging Two assumptions are made to model the array system: first all array elements are regarded as a rigid pistons with infinite baffles, and second an ultrasonic response from a target is modeled as a sum of contributions from single point targets. Linear systems theory is then used to model the array transmission and reception processes. The pressure pðr; tÞ at an observation point r is modeled as a convolution of the transmit spatial (velocity potential) impulse response hSIR t ðr; tÞ, the transmit aperture’s electrical impulse response hte ðtÞ, and the input signal ui ðtÞ, t pðr; tÞ ¼ hSIR t ðr; tÞ he ðtÞ ui ðtÞ
ð1Þ
where denotes temporal convolution [8]. Note that for simplicity it is assumed that the same input signal is applied to all transmit elements in Eq. (1) and that all transmit elements have the same (forward) electrical impulse response. Focusing is modeled by applying a different time delays to the transmit SIRs for each individual transmit element and then combining them into one transmit SIR hSIR t ðr; tÞ. The signal from a single receiving element, received for a point target at r, is modeled as a convolution of the pressure waveform pðr; tÞ with the receive SIR hSIR r ðr; tÞ
and the (backward) electrical impulse response hre ðtÞ of the receiving element r uo ðtÞ ¼ se hSIR r ðr; tÞ he ðtÞ pðr; tÞ;
ð2Þ
where se is the elementary surface of the point target [8]. Now, define the electrical impulse response, he ðtÞ, as D
he ðtÞ ¼ hte ðtÞ hre ðtÞ ui ðtÞ
ð3Þ
and let T be the set of points where the point targets are located. Then, using Eqs. (1)–(3) the output signal uo ðtÞ can be expressed as transmit
receive
X zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflffl{ t SIR r uo ðtÞ ¼ se hSIR t ðr; tÞ he ðtÞ ui ðtÞ hr ðr; tÞ he ðtÞ r2T
¼ se
X
SIR hSIR t ðr; tÞ hr ðr; tÞ he ðtÞ:
ð4Þ
r2T
In a discrete model of the imaging system the output signals from each receiving element, the SIRs, and the electrical impulse response are represented by vectors. The region-of-interest (ROI) is defined by a set of image points (zm , xn ) where m ¼ 0; 1; . . . ; M 1 and n ¼ 0; 1; . . . ; N 1. Here a rectangular ROI is used which is represented by a M N image matrix O. The output waveform received from a point target at point ðm; nÞ by the lth receive element T ðm;nÞ yl ¼ uðlÞ ; ð5Þ uðlÞ uðlÞ o ð0Þ o ð1Þ o ðK 1Þ will be given by a discrete convolution of the discrete double-path (transmit and receive) SIR corresponding to point ðm; nÞ and the discrete electrical impulse response scaled with ðOÞm;n , where ðOÞm;n denotes element ðm; nÞ in O and T denotes the transpose. If the system impulse response is defined as the combined transmit and receive response for point ðm; nÞ, D
SIR hðr; tÞ¼hSIR t ðr; tÞ hr ðr; tÞ he ðtÞ
then the sampled received signal can be expressed as ðm;nÞ
yl
ðm;nÞ
¼ hl
ðOÞm;n ;
ðm;nÞ yl
ð6Þ from element l ð7Þ
ðm;nÞ hl
where the vector is the discrete system impulse response corresponding to element l and point ðm; nÞ. The lth output for all targets in the ROI can now be found by super-position yl ¼
M 1 X N 1 X
ðm;nÞ
hl
ðOÞm;n :
ð8Þ
m¼0 n¼0
Eq. (8) can be expressed using matrix formalism by rearranging O in lexicographic order (i.e. vectorizing O) according to yl ¼ hlð0;0Þ hlð1;0Þ hlðM1;0Þ hlð0;1Þ hlðM1;N 1Þ o ¼ Pl o: ð9Þ
F. Lingvall / Ultrasonics 42 (2004) 961–968
By appending the received signals from all L elements T into a column vector y ¼ ½ yT0 yT1 yTL1 the array imaging system can finally be expressed as 2 32 3 2 3 o0 e0 P0 6 P1 76 o1 7 6 e1 7 6 76 7 6 7 ð10Þ y ¼ 6 .. 76 .. 7 þ 6 .. 7 ¼ Po þ e; 4 . 54 . 5 4 . 5 PL1
oN 1
eL1
where e is the measurement noise vector. The KL MN so called propagation matrix P in Eq. (10) describes both the transmission and the reception process for the array system. Note that Eq. (10) is a generalization of the model used in [6], since in the model in [6] the transmit and receive SIRs are the same because a single scanned transducer is normally used in synthetic aperture imaging. The model (10) is, however, valid for any transmit and receive aperture (provided that the corresponding receive and transmit SIRs can be computed).
963
linear array geometries). In a sampled system, where for simplicity apodization is not used, one can illustrate DAS processing by applying a binary mask over the measured vector y, such that, it has ones on positions corresponding the respectively time-delays and zeros otherwise. Then the data points corresponding to the positions with ones are summed. This mask is then computed (matched) for a single target located at the focus point. Thus, a discrete implementation of a DAS receive beamformer, focused at point ðm; nÞ, can be expressed as b DAS Þ ¼ dT y; ðO m;n m;n
ð11Þ
where d is a vector of length KL containing only ones and zeros. A parallel DAS beamformer can be obtained by extending Eq. (11) to handle focusing delays for all points in the ROI according to ^oDAS ¼ ½ d0;0 ¼ DT y:
d1;0
dM1;0
d0;1
T
dM1;N 1 y ð12Þ
2.2. Receive focusing 2.2.1. Delay-and-sum focusing The most common approach to focusing an ultrasonic array is the DAS method that performs focusing by compensating for the different time-of-flight from each transducer element to a focusing point (time shifting) and then performs a summation. The transducer elements are then implicitly treated as a point elements with SIRs in the form of Dirac pulses. This procedure is equivalent to a two-dimensional spatiotemporal filter that has all delta functions in the temporal direction and usually a hyperbolic shape (for
where D is a KL MN matrix and ^oDAS is the vectorized b DAS . The operation of the parallel DAS matrix O beamformer is illustrated in Fig. 1(a). 2.2.2. Minimum mean squared error focusing The MMSE approach, previously used for synthetic aperture imaging [6], consists in designing a linear filter which minimizes the overall reconstruction error. That is, the filter matrix K is the best linear filter in the sense that when the data y is filtered with K the resulting error in the estimate ^o ¼ Ky should be minimal measured with the mean squared error (MSE) criteria
Fig. 1. Illustration of the operation of the parallel delay-and-sum receive beamformer and the MMSE beamformer. (a) The parallel delay-and-sum Beamformer. (b) The MMSE beamformer.
964
F. Lingvall / Ultrasonics 42 (2004) 961–968 16 Element Linear Array
2
JMSE ¼ Efko Kyk g
T
¼ trfCoo g 2trfKPCoo g þ tr KðPCoo PT þ Cee ÞK ; ð13Þ where tr{} is the trace operator and Efg is the expectation operator. In Eq. (13) it is assumed that o and e are mutually independent zero-mean Gaussian processes with covariance matrices Coo and Cee , respectively. Thus, the prior information of the variation and correlation of the target strengths is given by the covariance matrix Coo and the statistical properties of the measurement noise is given by the noise covariance matrix Cee . The filter matrix K that minimizes (13) is the MMSE estimator for the linear model (10) [9]
x φ Steering Angle
Wire target ROI
z
Fig. 2. Illustration of the measurement setup.
KMMSE ¼ arg min JMSE ¼ Coo PT ðPCoo PT þ Cee Þ
1
k
1
T 1 1 ¼ ðCoo þ PT C1 ee PÞ P Cee :
ð14Þ
The MMSE filter (14) is the optimal filter, for the model (10), for all points in the ROI. That is, each row in KMMSE holds the optimal filter for compensating both the transmit and receive acoustic propagation effects as well as the electrical impulse response for the corresponding point in O. The MMSE filter KMMSE is therefore the optimal receive beamformer, in the mean squared sense, regardless how the array is focused in transmit, the size and shape of the elements, and what input signal that is used. Note that since the MMSE filter computes an estimate of o for all points simultaneously it is possible to beamform a full image from one emission. This operation can be viewed as a filter bank of optimal filters for all points which is illustrated in Fig. 1(b). Note also that there is a conceptual difference between the DAS beamforming approach and the MMSE approach. The DAS approach is aimed at maximizing the energy at the focusing point while the MMSE approach minimizes the overall reconstruction error. Thus, a DAS processed image may have a higher SNR at the focusing point compared to MMSE. However, an image processed by the MMSE filter will always have better (or equal) temporal and lateral resolution than a DAS processed image [6]. In particular, the MMSE method can compensate for finite sized transducer elements as well as for undersampled arrays.
array used in the experiments has 64 elements that are cylindrically focused at z ¼ 190 mm, where each element is 0.9 · 33 mm. Here, 16 of the array’s 64 elements were used and the beam was focused at 50 mm depth and steered )20 in transmit (see illustration in Fig. 2). The (transmit) electrical impulse response was measured using a Precision Acoustics PVDF hydrophone. The double-path electrical impulse response, shown in Fig. 3, was found by assuming that the input signal ui ðtÞ was the delta function (e.g. high voltage transient) and that the transmit and receive electrical impulse responses were identical; hence the double-path electrical impulse response can be expressed as ðtÞ he ðtÞ ¼ hðtÞ e ðtÞ he ðtÞ:
ð15Þ
The center frequency of the electrical impulse response was approximately 3 MHz and the distance between the array elements was 1 mm. The room temperature wavelength, k, in water at 3 MHZ is about 0.5 mm, hence the width of the array elements and the element spacing is approximately 2k, which means that the array is undersampled. The SIRs used in the model (10) were computed using the numerical DREAM method which outputs SIRs directly in a discrete time format [11]. The filter matrix for the MMSE beamformer was computed using the a priori assumption that all points in o were uncorrelated with variance ko ðCoo ¼ ko IÞ, and that the measurement noise e was white with variance ke ðCee ¼ ke IÞ, where I is the identity matrix. Using these assumptions Eq. (14) reduces to
3. Experiments
KMMSE ¼ PT ðPPT þ lIÞ1 ¼ ðPT P þ lIÞ1 PT ;
The MMSE and DAS beamformers have been compared using both simulated data, generated using the model (10), and measured data acquired using an ALLIN phased array system [10]. All measurements were performed in water using a 0.3 mm wire target. The
where l ¼ ke =ko . The l parameter in Eq. (16) can be seen as a measure of the ‘‘confidence’’ that the MMSE filter has to the measurements. If the measurement noise is very large then l will also be large and the factor lI 1 will dominate the inverse ðPPT þ lIÞ in Eq. (16). The
ð16Þ
F. Lingvall / Ultrasonics 42 (2004) 961–968
965
0 1
−10 Normalized Amplitude [dB]
Normalized Amplitude
0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6
−20 −30 −40 −50 −60 −70 −80
−0.8
−90
−1 0
(a)
0.5
1
1.5
2
2.5
3
3.5
−100 0
(b)
t [µs]
2
4
6
8
10
12
f [MHz]
Fig. 3. Double-path electrical impulse response of the ALLIN array measured with a hydrophone. (a) Double-path electrical impulse response. (b) Double-path electrical frequency response.
MMSE filter will then reduce to (a matched filter) PT =l. That is, the gain of the MMSE filter decreases when ke (and l) increases and the output will tend toward the a priori mean value of o which in this case is zero. In other words, the measurement vector y contains no information of about o. The other extreme case is when the measurement noise is zero, then Eq. (16) reduces to KMMSE ¼ PT ðPPT Þ
1
ð17Þ
and no prior information is used, that is, only the data is used for estimating o. The performance of the filter (17) is, however, most likely poor due to the ill-posed nature of the propagation operator P [4]. That is, when the matrix PPT has singular values close to zero even a small error in the measurement data will result in a large output ^ o. In fact, if KL < MN then PPT will have a number of zero singular values and the inverse will not even exist. I practice, modeling errors will also be present since the linear model does not perfectly models the true imaging system, and the l parameter can then be used as a tuning parameter to stabilize (16) and compensate for both modeling errors and the measurement noise. 3.1. Simulated data The simulated data was obtained by using the linear model (10) were the noise vector e was white with Gaussian amplitude distribution. The SNR is not uniquely defined since the received signal amplitude depends on the position of the target. Here the SNR is defined as the total energy received by all elements divided by the total noise energy yT y : ð18Þ eT e The signal energy yT y is computed using (10) with e ¼ 0 and is given by D
g¼
yT y ¼ oT PT Po:
ð19Þ
Note that, even for a constant noise energy eT e, g is a function of target position since y is a function of o. 3.2. Beamforming results for a single line As mentioned above the transmit beam is focused at 50 mm and steered )20. This will give a focus point at (z ¼ 50 mm, x 18 mm) but also a grating lobe at an angle of approximately 10 due to the fact that the element spacing is larger than k=2. In the first experiment the max amplitude of the beamformed image at the line corresponding to x ¼ 18 mm in the the ROI was plotted as a function of the target position. Both beamformers processed the same data and the results are shown in shown in Fig. 4. The receive DAS beamformer used a fixed receive focus at (z ¼ 50 mm, x ¼ 18 mm) and the MMSE beamformer used the optimal filters that corresponded to the points in the image line defined by x ¼ 18 mm. The noise level used in the simulation was kept a constant level for all target positions where the SNR, defined by Eq. (18), was approximately 15.8 dB when the target was at the focal point. To show the influence of the l parameter in the MMSE beamformer (16) the max amplitude was plotted, as shown in Fig. 4(a), as a function of target position on simulated data for five values of l, l ¼ 1010 , 1011 , 1012 , 1015 , and 1018 . For low values of l the MMSE filter relies strongly on the measured data resulting in a strong noise amplification. This is seen in Fig. 4(a) for the l ¼ 1010 where the noise floor only is about 7 dB lower than the target signal at focus. Increasing l lowers the noise sensitivity but it also results in a poorer suppression of the grating lobe which can be seen for l ¼ 1015 and 1018 . A good trade-off between the noise sensitivity and grating lobe suppression is in this case l ¼ 1012 where the noise and grating lobe levels are roughly 25 dB lower than the target signal. Note that the MMSE beamformer suppressed the
966
F. Lingvall / Ultrasonics 42 (2004) 961–968
(a)
(b)
Fig. 4. Normalized maximum amplitude as a function of horizontal (point) target position for a 16 element linear array focused at 50 mm and steered at an angle of )20. The target was located at depth z ¼ 50 mm and the SNR when the target is at x ¼ 18 mm was 15.8 dB. (a) Simulated data (SNR ¼ 15.8 [dB]). (b) Measured data (0.3 mm wire target).
grating lobe better for all the examined l values. It is also noticeable that the MMSE beamformed data has higher lateral resolution than the DAS beamformed data. The same behavior can be seen for the measured results shown in Fig. 4(b) where was l ¼ 1010 . The MMSE beamformer has both higher resolution and better grating lobe suppression even though model errors degraded the MMSE reconstructions somewhat. 3.3. Parallel beamforming results In the second experiment the parallel DAS beamformer (12) and the MMSE beamformer (16) was used
to beamform data from a single emission. Again the array was focused and steered at (z ¼ 50 mm, x 18 mm) in transmission. The performance for two target positions was examined, the first with the target in the main lobe, shown in Fig. 5, and the second with the target located in the grating lobe, which is shown in Fig. 6. The SNR in the simulated data was 15.8 dB for both target positions and l ¼ 1012 was used for the MMSE beamformer at both target positions. Fig. 5(a) shows the simulated data, received by the 16 elements, when the target is in the main lobe and Fig. 5(b) and (c) show the processed images from simulated data with the
Fig. 5. Beamformed data for a point target in the main lobe using a single emission of a 16 element linear array focused at 50 mm and steered at an angle of )20. (a) Simulated data (SNR ¼ 15.8 [dB]). (b) Parallel DAS on simulated data. (c) MMSE on simulated data (l ¼ 1012 ). (d) Measured data. (e) Parallel DAS on measured data. (f) MMSE on measured data (l ¼ 1012 ).
F. Lingvall / Ultrasonics 42 (2004) 961–968
967
Fig. 6. Beamformed data for a point target in the grating lobe using a single emission of a 16 element linear array focused at 50 mm and steered at an angle of )20. (a) Simulated data (SNR ¼ 15.8 [dB]). (b) Parallel DAS on simulated data. (c) MMSE on simulated data (l ¼ 1012 ). (d) Measured data. (e) Parallel DAS on measured data. (f) MMSE on measured data (l ¼ 1012 ).
parallel DAS and the MMSE beamformer, respectively. Fig. 5(d) and (e) show the corresponding images for the measured data. Similarly, Fig. 6(a) presents the simulated data for the target in the grating lobe and Fig. 6(b) and (c) the processed images with the parallel DAS and the MMSE for the simulated data. Finally, Fig. 6(d) and (e) shows the corresponding images obtained for the measured data. As can be seen in Fig. 5(b), (c), (e), and (f) both beamformers perform rather well when the target is in the main lobe. However, the results when the target is in the grating lobe, shown in Fig. 6(b), (c), (e), and (f), the parallel DAS beamformer produces images of poor quality with strong artifacts from the grating lobe. The MMSE beamformer, on the other hand, performs well, with a high temporal and lateral resolution, regardless if the target is in the main or the grating lobe.
4. Conclusions A MMSE based ultrasonic array receive beamformer has been presented and compared to the conventional DAS beamformer in this paper. The two beamformers have been compared using both simulated and measured data obtained from a linear phased array system. It was shown that the MMSE algorithm has superior performance in terms of higher resolution and contrast compared to the DAS beamformer. The DAS beamformer was in particular unable to suppress grating lobes in the transmit beam which resulted in strong artifacts in the
beamformed images. The MMSE beamformer on the another hand, accounting for both the acoustic transmission and reception effects, produced images with high temporal and lateral resolution characterized by low levels of artifacts despite that the transmit beam had strong acoustic energy in the grating lobe. It was also shown that the MMSE beamformer can be viewed as a bank of optimal filters for every image point and it is thus possible to perform the MMSE filtering in parallel to obtain a high frame rate. The extra cost of the MMSE method compared to the DAS approach is that all data points are used in processing for every image point whereas the DAS method only uses those data points that corresponds to the respective time-delays. Acknowledgements This work was sponsored by the Swedish Nuclear Fuel and Waste Management Co. (SKB). References [1] G.S. Kino, Acoustic Waves: Devices, Imaging and Analog Signal Processing, volume 6 of Prentice-Hall Signal Processing Series. Prentice-Hall, 1987. [2] D.P. Shattuck, M.D. Weinshenker, S.W. Smith, O.T. von Ramm, Explososcan: a parallel processing technique for high speed ultrasound imaging with linear phased arrays, J. Acoust. Soc. Am. 75 (4) (1984) 1273–1282. [3] O.T. von Ram, S.W. Smith, H.G. Pavy Jr., High-speed ultrasound volumetric imaging system. II. Parallel processing and image
968
[4]
[5]
[6]
[7]
F. Lingvall / Ultrasonics 42 (2004) 961–968 display, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 38 (2) (1991) 109–115. M. Tanter, J.-F. Aubry, J. Gerber, J.-L. Thomas, M. Fink, Optimal focusing by spatio-temporal inverse filter. I. Basic principles, J. Acoust. Soc. Am. (2001) 37–47. J. Shen, E.S. Ebbini, Filter-based coded-excitation system for high-speed ultrasonic imaging, IEEE Trans. Med. Imaging 17 (6) (1998) 923–934. F. Lingvall, T. Olofsson, T. Stepinski, Synthetic aperture imaging using sources with finite aperture: deconvolution of the spatial impulse response, J. Acoust. Soc. Am. 114 (1) (2003) 225–234. P.R. Stepanishen, Transient radiation from pistons in an infinite planar baffle, J. Acoust. Soc. Am. 49 (1971) 1629–1638.
[8] A. Lhemery, Impulse-response method to predict echo-responses from targets of complex geometry. Part I: Theory, J. Acoust. Soc. Am. (1991) 2799–2807. [9] S.M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory, vol. I, Prentice Hall, 1993. [10] P. Wu, T. Stepinski, Spatial impulse response method for predicting pulse-echo fields from a linear array with cylindrically concave surface, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 46 (5) (1999) 1283–1290. [11] B. Piwakowski, K. Sbai, A new approach to calculate the field radiated from arbitrarily structured transducer arrays, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 46 (2) (1999) 422– 440.
A method of improving overall resolution in ultrasonic array imaging using spatio-temporal deconvolution Fredrik Lingvall Department of Engineering Sciences, Signals and Systems Group, Uppsala University, Box 528, 751 20 Uppsala, Sweden
Abstract In this paper a beamforming method for ultrasonic array imaging is presented that performs both spatial and temporal deconvolution based on a minimum mean square error (MMSE) criteria. The presented MMSE receive mode beamformer performs a regularized inversion of the propagation operator for the ultrasonic array system at hand. The MMSE beamformer accounts for the transmit and receive processes, defined in terms of finite array element sizes, transmit focusing laws and electrical transducer characteristics. The MMSE beamformer is compared to the traditional delay-and-sum (DAS) beamformer with respect to both resolution and signal-to-noise ratio. The two algorithms are compared using both simulated and measured data. The simulated data was obtained using ultrasonic field simulations and the measured data was acquired using a linear phased array imaging wire targets in water. The results show that the MMSE beamformer has superior temporal and lateral resolution compared to DAS. It is also shown that the MMSE beamformer can be expressed as a filter bank, which enables parallel processing at high frame rates. 2003 Elsevier B.V. All rights reserved. Keywords: Ultrasonic array imaging; Spatio-temporal deconvolution; Beamforming
1. Introduction In conventional ultrasonic array imaging focusing in receptions is performed by means of delay-and-sum (DAS) operations on the signals received by the array elements. To successfully apply DAS receive mode beamforming the array elements must be small compared to the wavelength, otherwise, the diffraction effects will degrade the processed images [1]. There is, however, a trade-off between element size and signal-tonoise ratio (SNR) and the elements must therefore be large enough to generate a sufficient amount of acoustic energy to obtain a sufficiently high SNR. Another important issue that also limits the element size in DAS beamforming is that, if undersampled arrays with element spacing larger than half the wavelength are used, the images will be degraded by the grating lobes. Avoiding grating lobes by limiting the element spacing to less than half the wavelength also limits the element size and in consequence also the transmitted acoustic power. In order to improve the SNR it is, therefore,
E-mail address: fl@signal.uu.se (F. Lingvall). 0041-624X/$ - see front matter 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.ultras.2003.12.016
desirable to increase the size of the array elements. Even if advanced methods such as dynamic receive focusing and parallel beamforming [2,3], also based on DAS processing, are used to improve image quality and frame-rate the fundamental limitations of DAS approach can not be circumvented. Another deficiency of DAS beamforming is that the electrical impulse response of the transducer is not taken into consideration and the DAS processing will, therefore, not improve temporal resolution. Recently, methods that are based on inverse filtering of the propagation operator have been introduced to overcome the limitations of the DAS approach [4,5]. Since the propagation operator describes both the transmission and reception processes, such an approach could in theory provide full information of the scattering objects in the insonified region of interest. However, the inverse operator is generally ill-posed and regularized solutions must be sought, for instance sub-space methods [4,5]. A serious disadvantage of sub-space regularization methods is that they do not incorporate any prior knowledge, such as, the statistical distribution of the measurement noise and the fact that the target strength is limited.
962
F. Lingvall / Ultrasonics 42 (2004) 961–968
In this paper a recently introduced inverse filtering approach for synthetic aperture imaging [6], is generalized and applied to ultrasonic array imaging. This new approach, that uses a priori statistical information about the scattering targets and the measurement noise, uses a model based method and a minimum mean squared error (MMSE) criteria to find a solution that minimizes the reconstruction error for the array system at hand. The linear model takes both the acoustic transmission and reception processes into account, using the general concept of spatial impulse responses [7], as well as the electrical impulse response of the transducer. The application of MMSE method is, therefore, not limited to any particular constraints regarding the element size, element spacing, or transmit focusing method. The MMSE method can, for example, be applied to undersampled arrays and arrays with finite sized elements. Moreover, since the linear model includes the electrical impulse response the temporal resolution will also be improved. This paper is organized as follows: The linear discrete model and the MMSE filter is described in Section 2. The DAS and the MMSE methods are then compared, using both simulated and measured data, in Section 3. Finally, the conclusions are given in Section 4.
2. Theory 2.1. A discrete linear model for ultrasonic array imaging Two assumptions are made to model the array system: first all array elements are regarded as a rigid pistons with infinite baffles, and second an ultrasonic response from a target is modeled as a sum of contributions from single point targets. Linear systems theory is then used to model the array transmission and reception processes. The pressure pðr; tÞ at an observation point r is modeled as a convolution of the transmit spatial (velocity potential) impulse response hSIR t ðr; tÞ, the transmit aperture’s electrical impulse response hte ðtÞ, and the input signal ui ðtÞ, t pðr; tÞ ¼ hSIR t ðr; tÞ he ðtÞ ui ðtÞ
ð1Þ
where denotes temporal convolution [8]. Note that for simplicity it is assumed that the same input signal is applied to all transmit elements in Eq. (1) and that all transmit elements have the same (forward) electrical impulse response. Focusing is modeled by applying a different time delays to the transmit SIRs for each individual transmit element and then combining them into one transmit SIR hSIR t ðr; tÞ. The signal from a single receiving element, received for a point target at r, is modeled as a convolution of the pressure waveform pðr; tÞ with the receive SIR hSIR r ðr; tÞ
and the (backward) electrical impulse response hre ðtÞ of the receiving element r uo ðtÞ ¼ se hSIR r ðr; tÞ he ðtÞ pðr; tÞ;
ð2Þ
where se is the elementary surface of the point target [8]. Now, define the electrical impulse response, he ðtÞ, as D
he ðtÞ ¼ hte ðtÞ hre ðtÞ ui ðtÞ
ð3Þ
and let T be the set of points where the point targets are located. Then, using Eqs. (1)–(3) the output signal uo ðtÞ can be expressed as transmit
receive
X zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ zfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflffl{ t SIR r uo ðtÞ ¼ se hSIR t ðr; tÞ he ðtÞ ui ðtÞ hr ðr; tÞ he ðtÞ r2T
¼ se
X
SIR hSIR t ðr; tÞ hr ðr; tÞ he ðtÞ:
ð4Þ
r2T
In a discrete model of the imaging system the output signals from each receiving element, the SIRs, and the electrical impulse response are represented by vectors. The region-of-interest (ROI) is defined by a set of image points (zm , xn ) where m ¼ 0; 1; . . . ; M 1 and n ¼ 0; 1; . . . ; N 1. Here a rectangular ROI is used which is represented by a M N image matrix O. The output waveform received from a point target at point ðm; nÞ by the lth receive element T ðm;nÞ yl ¼ uðlÞ ; ð5Þ uðlÞ uðlÞ o ð0Þ o ð1Þ o ðK 1Þ will be given by a discrete convolution of the discrete double-path (transmit and receive) SIR corresponding to point ðm; nÞ and the discrete electrical impulse response scaled with ðOÞm;n , where ðOÞm;n denotes element ðm; nÞ in O and T denotes the transpose. If the system impulse response is defined as the combined transmit and receive response for point ðm; nÞ, D
SIR hðr; tÞ¼hSIR t ðr; tÞ hr ðr; tÞ he ðtÞ
then the sampled received signal can be expressed as ðm;nÞ
yl
ðm;nÞ
¼ hl
ðOÞm;n ;
ðm;nÞ yl
ð6Þ from element l ð7Þ
ðm;nÞ hl
where the vector is the discrete system impulse response corresponding to element l and point ðm; nÞ. The lth output for all targets in the ROI can now be found by super-position yl ¼
M 1 X N 1 X
ðm;nÞ
hl
ðOÞm;n :
ð8Þ
m¼0 n¼0
Eq. (8) can be expressed using matrix formalism by rearranging O in lexicographic order (i.e. vectorizing O) according to yl ¼ hlð0;0Þ hlð1;0Þ hlðM1;0Þ hlð0;1Þ hlðM1;N 1Þ o ¼ Pl o: ð9Þ
F. Lingvall / Ultrasonics 42 (2004) 961–968
By appending the received signals from all L elements T into a column vector y ¼ ½ yT0 yT1 yTL1 the array imaging system can finally be expressed as 2 32 3 2 3 o0 e0 P0 6 P1 76 o1 7 6 e1 7 6 76 7 6 7 ð10Þ y ¼ 6 .. 76 .. 7 þ 6 .. 7 ¼ Po þ e; 4 . 54 . 5 4 . 5 PL1
oN 1
eL1
where e is the measurement noise vector. The KL MN so called propagation matrix P in Eq. (10) describes both the transmission and the reception process for the array system. Note that Eq. (10) is a generalization of the model used in [6], since in the model in [6] the transmit and receive SIRs are the same because a single scanned transducer is normally used in synthetic aperture imaging. The model (10) is, however, valid for any transmit and receive aperture (provided that the corresponding receive and transmit SIRs can be computed).
963
linear array geometries). In a sampled system, where for simplicity apodization is not used, one can illustrate DAS processing by applying a binary mask over the measured vector y, such that, it has ones on positions corresponding the respectively time-delays and zeros otherwise. Then the data points corresponding to the positions with ones are summed. This mask is then computed (matched) for a single target located at the focus point. Thus, a discrete implementation of a DAS receive beamformer, focused at point ðm; nÞ, can be expressed as b DAS Þ ¼ dT y; ðO m;n m;n
ð11Þ
where d is a vector of length KL containing only ones and zeros. A parallel DAS beamformer can be obtained by extending Eq. (11) to handle focusing delays for all points in the ROI according to ^oDAS ¼ ½ d0;0 ¼ DT y:
d1;0
dM1;0
d0;1
T
dM1;N 1 y ð12Þ
2.2. Receive focusing 2.2.1. Delay-and-sum focusing The most common approach to focusing an ultrasonic array is the DAS method that performs focusing by compensating for the different time-of-flight from each transducer element to a focusing point (time shifting) and then performs a summation. The transducer elements are then implicitly treated as a point elements with SIRs in the form of Dirac pulses. This procedure is equivalent to a two-dimensional spatiotemporal filter that has all delta functions in the temporal direction and usually a hyperbolic shape (for
where D is a KL MN matrix and ^oDAS is the vectorized b DAS . The operation of the parallel DAS matrix O beamformer is illustrated in Fig. 1(a). 2.2.2. Minimum mean squared error focusing The MMSE approach, previously used for synthetic aperture imaging [6], consists in designing a linear filter which minimizes the overall reconstruction error. That is, the filter matrix K is the best linear filter in the sense that when the data y is filtered with K the resulting error in the estimate ^o ¼ Ky should be minimal measured with the mean squared error (MSE) criteria
Fig. 1. Illustration of the operation of the parallel delay-and-sum receive beamformer and the MMSE beamformer. (a) The parallel delay-and-sum Beamformer. (b) The MMSE beamformer.
964
F. Lingvall / Ultrasonics 42 (2004) 961–968 16 Element Linear Array
2
JMSE ¼ Efko Kyk g
T
¼ trfCoo g 2trfKPCoo g þ tr KðPCoo PT þ Cee ÞK ; ð13Þ where tr{} is the trace operator and Efg is the expectation operator. In Eq. (13) it is assumed that o and e are mutually independent zero-mean Gaussian processes with covariance matrices Coo and Cee , respectively. Thus, the prior information of the variation and correlation of the target strengths is given by the covariance matrix Coo and the statistical properties of the measurement noise is given by the noise covariance matrix Cee . The filter matrix K that minimizes (13) is the MMSE estimator for the linear model (10) [9]
x φ Steering Angle
Wire target ROI
z
Fig. 2. Illustration of the measurement setup.
KMMSE ¼ arg min JMSE ¼ Coo PT ðPCoo PT þ Cee Þ
1
k
1
T 1 1 ¼ ðCoo þ PT C1 ee PÞ P Cee :
ð14Þ
The MMSE filter (14) is the optimal filter, for the model (10), for all points in the ROI. That is, each row in KMMSE holds the optimal filter for compensating both the transmit and receive acoustic propagation effects as well as the electrical impulse response for the corresponding point in O. The MMSE filter KMMSE is therefore the optimal receive beamformer, in the mean squared sense, regardless how the array is focused in transmit, the size and shape of the elements, and what input signal that is used. Note that since the MMSE filter computes an estimate of o for all points simultaneously it is possible to beamform a full image from one emission. This operation can be viewed as a filter bank of optimal filters for all points which is illustrated in Fig. 1(b). Note also that there is a conceptual difference between the DAS beamforming approach and the MMSE approach. The DAS approach is aimed at maximizing the energy at the focusing point while the MMSE approach minimizes the overall reconstruction error. Thus, a DAS processed image may have a higher SNR at the focusing point compared to MMSE. However, an image processed by the MMSE filter will always have better (or equal) temporal and lateral resolution than a DAS processed image [6]. In particular, the MMSE method can compensate for finite sized transducer elements as well as for undersampled arrays.
array used in the experiments has 64 elements that are cylindrically focused at z ¼ 190 mm, where each element is 0.9 · 33 mm. Here, 16 of the array’s 64 elements were used and the beam was focused at 50 mm depth and steered )20 in transmit (see illustration in Fig. 2). The (transmit) electrical impulse response was measured using a Precision Acoustics PVDF hydrophone. The double-path electrical impulse response, shown in Fig. 3, was found by assuming that the input signal ui ðtÞ was the delta function (e.g. high voltage transient) and that the transmit and receive electrical impulse responses were identical; hence the double-path electrical impulse response can be expressed as ðtÞ he ðtÞ ¼ hðtÞ e ðtÞ he ðtÞ:
ð15Þ
The center frequency of the electrical impulse response was approximately 3 MHz and the distance between the array elements was 1 mm. The room temperature wavelength, k, in water at 3 MHZ is about 0.5 mm, hence the width of the array elements and the element spacing is approximately 2k, which means that the array is undersampled. The SIRs used in the model (10) were computed using the numerical DREAM method which outputs SIRs directly in a discrete time format [11]. The filter matrix for the MMSE beamformer was computed using the a priori assumption that all points in o were uncorrelated with variance ko ðCoo ¼ ko IÞ, and that the measurement noise e was white with variance ke ðCee ¼ ke IÞ, where I is the identity matrix. Using these assumptions Eq. (14) reduces to
3. Experiments
KMMSE ¼ PT ðPPT þ lIÞ1 ¼ ðPT P þ lIÞ1 PT ;
The MMSE and DAS beamformers have been compared using both simulated data, generated using the model (10), and measured data acquired using an ALLIN phased array system [10]. All measurements were performed in water using a 0.3 mm wire target. The
where l ¼ ke =ko . The l parameter in Eq. (16) can be seen as a measure of the ‘‘confidence’’ that the MMSE filter has to the measurements. If the measurement noise is very large then l will also be large and the factor lI 1 will dominate the inverse ðPPT þ lIÞ in Eq. (16). The
ð16Þ
F. Lingvall / Ultrasonics 42 (2004) 961–968
965
0 1
−10 Normalized Amplitude [dB]
Normalized Amplitude
0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6
−20 −30 −40 −50 −60 −70 −80
−0.8
−90
−1 0
(a)
0.5
1
1.5
2
2.5
3
3.5
−100 0
(b)
t [µs]
2
4
6
8
10
12
f [MHz]
Fig. 3. Double-path electrical impulse response of the ALLIN array measured with a hydrophone. (a) Double-path electrical impulse response. (b) Double-path electrical frequency response.
MMSE filter will then reduce to (a matched filter) PT =l. That is, the gain of the MMSE filter decreases when ke (and l) increases and the output will tend toward the a priori mean value of o which in this case is zero. In other words, the measurement vector y contains no information of about o. The other extreme case is when the measurement noise is zero, then Eq. (16) reduces to KMMSE ¼ PT ðPPT Þ
1
ð17Þ
and no prior information is used, that is, only the data is used for estimating o. The performance of the filter (17) is, however, most likely poor due to the ill-posed nature of the propagation operator P [4]. That is, when the matrix PPT has singular values close to zero even a small error in the measurement data will result in a large output ^ o. In fact, if KL < MN then PPT will have a number of zero singular values and the inverse will not even exist. I practice, modeling errors will also be present since the linear model does not perfectly models the true imaging system, and the l parameter can then be used as a tuning parameter to stabilize (16) and compensate for both modeling errors and the measurement noise. 3.1. Simulated data The simulated data was obtained by using the linear model (10) were the noise vector e was white with Gaussian amplitude distribution. The SNR is not uniquely defined since the received signal amplitude depends on the position of the target. Here the SNR is defined as the total energy received by all elements divided by the total noise energy yT y : ð18Þ eT e The signal energy yT y is computed using (10) with e ¼ 0 and is given by D
g¼
yT y ¼ oT PT Po:
ð19Þ
Note that, even for a constant noise energy eT e, g is a function of target position since y is a function of o. 3.2. Beamforming results for a single line As mentioned above the transmit beam is focused at 50 mm and steered )20. This will give a focus point at (z ¼ 50 mm, x 18 mm) but also a grating lobe at an angle of approximately 10 due to the fact that the element spacing is larger than k=2. In the first experiment the max amplitude of the beamformed image at the line corresponding to x ¼ 18 mm in the the ROI was plotted as a function of the target position. Both beamformers processed the same data and the results are shown in shown in Fig. 4. The receive DAS beamformer used a fixed receive focus at (z ¼ 50 mm, x ¼ 18 mm) and the MMSE beamformer used the optimal filters that corresponded to the points in the image line defined by x ¼ 18 mm. The noise level used in the simulation was kept a constant level for all target positions where the SNR, defined by Eq. (18), was approximately 15.8 dB when the target was at the focal point. To show the influence of the l parameter in the MMSE beamformer (16) the max amplitude was plotted, as shown in Fig. 4(a), as a function of target position on simulated data for five values of l, l ¼ 1010 , 1011 , 1012 , 1015 , and 1018 . For low values of l the MMSE filter relies strongly on the measured data resulting in a strong noise amplification. This is seen in Fig. 4(a) for the l ¼ 1010 where the noise floor only is about 7 dB lower than the target signal at focus. Increasing l lowers the noise sensitivity but it also results in a poorer suppression of the grating lobe which can be seen for l ¼ 1015 and 1018 . A good trade-off between the noise sensitivity and grating lobe suppression is in this case l ¼ 1012 where the noise and grating lobe levels are roughly 25 dB lower than the target signal. Note that the MMSE beamformer suppressed the
966
F. Lingvall / Ultrasonics 42 (2004) 961–968
(a)
(b)
Fig. 4. Normalized maximum amplitude as a function of horizontal (point) target position for a 16 element linear array focused at 50 mm and steered at an angle of )20. The target was located at depth z ¼ 50 mm and the SNR when the target is at x ¼ 18 mm was 15.8 dB. (a) Simulated data (SNR ¼ 15.8 [dB]). (b) Measured data (0.3 mm wire target).
grating lobe better for all the examined l values. It is also noticeable that the MMSE beamformed data has higher lateral resolution than the DAS beamformed data. The same behavior can be seen for the measured results shown in Fig. 4(b) where was l ¼ 1010 . The MMSE beamformer has both higher resolution and better grating lobe suppression even though model errors degraded the MMSE reconstructions somewhat. 3.3. Parallel beamforming results In the second experiment the parallel DAS beamformer (12) and the MMSE beamformer (16) was used
to beamform data from a single emission. Again the array was focused and steered at (z ¼ 50 mm, x 18 mm) in transmission. The performance for two target positions was examined, the first with the target in the main lobe, shown in Fig. 5, and the second with the target located in the grating lobe, which is shown in Fig. 6. The SNR in the simulated data was 15.8 dB for both target positions and l ¼ 1012 was used for the MMSE beamformer at both target positions. Fig. 5(a) shows the simulated data, received by the 16 elements, when the target is in the main lobe and Fig. 5(b) and (c) show the processed images from simulated data with the
Fig. 5. Beamformed data for a point target in the main lobe using a single emission of a 16 element linear array focused at 50 mm and steered at an angle of )20. (a) Simulated data (SNR ¼ 15.8 [dB]). (b) Parallel DAS on simulated data. (c) MMSE on simulated data (l ¼ 1012 ). (d) Measured data. (e) Parallel DAS on measured data. (f) MMSE on measured data (l ¼ 1012 ).
F. Lingvall / Ultrasonics 42 (2004) 961–968
967
Fig. 6. Beamformed data for a point target in the grating lobe using a single emission of a 16 element linear array focused at 50 mm and steered at an angle of )20. (a) Simulated data (SNR ¼ 15.8 [dB]). (b) Parallel DAS on simulated data. (c) MMSE on simulated data (l ¼ 1012 ). (d) Measured data. (e) Parallel DAS on measured data. (f) MMSE on measured data (l ¼ 1012 ).
parallel DAS and the MMSE beamformer, respectively. Fig. 5(d) and (e) show the corresponding images for the measured data. Similarly, Fig. 6(a) presents the simulated data for the target in the grating lobe and Fig. 6(b) and (c) the processed images with the parallel DAS and the MMSE for the simulated data. Finally, Fig. 6(d) and (e) shows the corresponding images obtained for the measured data. As can be seen in Fig. 5(b), (c), (e), and (f) both beamformers perform rather well when the target is in the main lobe. However, the results when the target is in the grating lobe, shown in Fig. 6(b), (c), (e), and (f), the parallel DAS beamformer produces images of poor quality with strong artifacts from the grating lobe. The MMSE beamformer, on the other hand, performs well, with a high temporal and lateral resolution, regardless if the target is in the main or the grating lobe.
4. Conclusions A MMSE based ultrasonic array receive beamformer has been presented and compared to the conventional DAS beamformer in this paper. The two beamformers have been compared using both simulated and measured data obtained from a linear phased array system. It was shown that the MMSE algorithm has superior performance in terms of higher resolution and contrast compared to the DAS beamformer. The DAS beamformer was in particular unable to suppress grating lobes in the transmit beam which resulted in strong artifacts in the
beamformed images. The MMSE beamformer on the another hand, accounting for both the acoustic transmission and reception effects, produced images with high temporal and lateral resolution characterized by low levels of artifacts despite that the transmit beam had strong acoustic energy in the grating lobe. It was also shown that the MMSE beamformer can be viewed as a bank of optimal filters for every image point and it is thus possible to perform the MMSE filtering in parallel to obtain a high frame rate. The extra cost of the MMSE method compared to the DAS approach is that all data points are used in processing for every image point whereas the DAS method only uses those data points that corresponds to the respective time-delays. Acknowledgements This work was sponsored by the Swedish Nuclear Fuel and Waste Management Co. (SKB). References [1] G.S. Kino, Acoustic Waves: Devices, Imaging and Analog Signal Processing, volume 6 of Prentice-Hall Signal Processing Series. Prentice-Hall, 1987. [2] D.P. Shattuck, M.D. Weinshenker, S.W. Smith, O.T. von Ramm, Explososcan: a parallel processing technique for high speed ultrasound imaging with linear phased arrays, J. Acoust. Soc. Am. 75 (4) (1984) 1273–1282. [3] O.T. von Ram, S.W. Smith, H.G. Pavy Jr., High-speed ultrasound volumetric imaging system. II. Parallel processing and image
968
[4]
[5]
[6]
[7]
F. Lingvall / Ultrasonics 42 (2004) 961–968 display, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 38 (2) (1991) 109–115. M. Tanter, J.-F. Aubry, J. Gerber, J.-L. Thomas, M. Fink, Optimal focusing by spatio-temporal inverse filter. I. Basic principles, J. Acoust. Soc. Am. (2001) 37–47. J. Shen, E.S. Ebbini, Filter-based coded-excitation system for high-speed ultrasonic imaging, IEEE Trans. Med. Imaging 17 (6) (1998) 923–934. F. Lingvall, T. Olofsson, T. Stepinski, Synthetic aperture imaging using sources with finite aperture: deconvolution of the spatial impulse response, J. Acoust. Soc. Am. 114 (1) (2003) 225–234. P.R. Stepanishen, Transient radiation from pistons in an infinite planar baffle, J. Acoust. Soc. Am. 49 (1971) 1629–1638.
[8] A. Lhemery, Impulse-response method to predict echo-responses from targets of complex geometry. Part I: Theory, J. Acoust. Soc. Am. (1991) 2799–2807. [9] S.M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory, vol. I, Prentice Hall, 1993. [10] P. Wu, T. Stepinski, Spatial impulse response method for predicting pulse-echo fields from a linear array with cylindrically concave surface, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 46 (5) (1999) 1283–1290. [11] B. Piwakowski, K. Sbai, A new approach to calculate the field radiated from arbitrarily structured transducer arrays, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 46 (2) (1999) 422– 440.