CC BY
65Electronic Journal «Technical Acoustics» http://www .ejta.org
2008, 10
M. Malekzadeh Shafaroudi*, H. Behnam**
* Electrical Engineering Dept., Iran University of Science and Technology, Tehran, Iran; e-mail: [email protected]
Professor of Electrical Engineering Dept., Iran University of Science and Technology, Tehran, Iran; e-mail: [email protected]
Ultrasound field simulation for array transducers based on impulse response calculation
Received 14.09.2007, published 26.04.2008
Array ultrasound transducers have very complicated structures and their field measurements consume a lot of time and cost; therefore using simulation programs is very helpful. In this paper a method for calculation of ultrasound field from continuous wave and pulse excitation for array transducers is presented. The method has been implemented as a simulator software. By defining array-transducer parameters such as elements dimensions, their distances and environment features, this software is capable of calculation of ultrasound pressure field and plotting the results in various methods. A special method (Ring Function) is used that computes field of arcs with the same distance on the transducer from the observation point for calculating of impulse response, and then considering the excitation function, with a simple convolution, pressure field can be calculated. Therefore it prevents repetitive and time-consuming calculation because after impulse response calculation the results can be used for any kind of excitation waveform. Also with some examples, the results of this software were compared with some available software such as Field II and Ultrasim.
Key words: ultrasound, ultrasound transducer, array transducer, ultrasound field, simulation.
INTRODUCTION
Today sound waves are utilized in many medical and non-medical applications like treatment and imaging. Various types of probes are used in ultrasound imaging. Some of the most important and the most common of them are 1D (one-dimensional) and 2D array transducers that are very important in two and three dimensional imaging. 2D array transducers are vastly used in 3D and real-time imaging such as echocardiography, urology, sonography and so on.
The role of the acoustic transducer in high quality image acquisition from medical ultrasound scanners is very important. In fact it determines the quality of the image. There are considerable efforts in designing transducers and determining the characteristics of the emitted field.
The goal of the 3D imaging is observing the real structure of the body organs. 2D-Array transducers are useful in acquiring real-time 3D images; in echocardiography, in monitoring
fetus and in urology. An array transducer removes the physical movement of 1D array in order to acquiring 3D images. It is able to focus the acoustic beam in the elevation direction in comparison with 1D array that uses lenses to do this. Calculating and simulating the acoustic fields are simpler and more economical than experimental tasks in order to determine the field pattern emitted from acoustic transducers.
There are different methods [1] for simulating the acoustic fields that most of them are based on the Rayleigh integral. (Field II by Jensen et al, [2-7] & Ultrasim by Holm et. al., [811]). Far-field approximation is mostly applied in simulation of acoustic fields. So there is not enough accuracy in near-field. In 2004 a new software was presented (Science and Research Branch of Iran Islamic Azad university [12]) capable of calculating the pressure field from array transducers stimulated by continuous wave.
In sonography sometimes imaging the tissues very near to the transducer, specially in intravascular sonography, is interesting so it will be very helpful to develop a calculating method that has a good accuracy both in near-field and in far-field.
In this paper a new software will be introduced, which is more efficient and has no approximation for near-field, that calculates the ultrasound field from array transducers. Here the Rayleigh integral is calculated by considering the points on the transducer surface that have the same distance from the observation point, so it reduces the double integral to a one-dimensional integral expression. This technique has been described by Schoch [13] and Ohtsuki [14] for impulse response and continuous wave excitation, respectively.
By the software introduced here that is named SimATra (Simulator of Array Transducers), both continuous and pulse wave fields can be calculated and simulated. The Ring Function method is used in SimAtra and it is possible to alter and modify the distance between the elements of array in each direction and after impulse response calculation the results can be used for deriving the results for any kind of excitation waveform.
1. BASIC THEORY
1.1. Rayleigh Integral
Figure 1 shows an aperture placed in an infinite, rigid baffle on which the velocity normal to the plane is zero, except at the aperture. The field point is denoted by r and the aperture by r2, c is the speed of sound, and t is time. The special impulse response generated by the aperture is then found by the Rayleigh integral [15]
The integral is a statement of Huygens' Principle that the field is found by integrating the contributions from all the infinitesimally small area elements that make up the aperture. This integral formulation assumes linearity and propagation in a homogeneous medium without attenuation. Further, the radiating aperture is assumed flat, so no re-radiation from scattering and reflection takes place.
(1)
Figure 1. Position of aperture placed in an infinite, rigid baffle, field point, and coordinate system. The field point is denoted by r1 and the aperture by r2
Using the spatial impulse response the pressure is written as
p(fl, t) = p ^dr * hfo t), (2)
dt
where vn (t) is the velocity normal to the plane.
Equation (1) is called the spatial impulse response and characterizes the three-dimensional extent of the field for a particular transducer geometry. Note that this is a function of the relative position between the aperture and the field.
The continuous wave field can be found from the Fourier transform of equation (2) (section 3.2). Thus, the calculation of the spatial impulse response makes it possible to find all ultrasound fields of interest.
1.2. Calculation of spatial impulse response
A short review of the calculation of spatial impulse responses is given in this section to facilitate the development of the new calculation procedure.
As mentioned before equation (1) is essentially a statement of Huygens's principle that the field is found by summing the radiated spherical waves from all parts of the aperture. This can also be reformulated, due to acoustic reciprocity, as finding the part of the spherical wave emanating from the field point that intersects the aperture. The task is, thus, to project the field point onto the plane coinciding with the aperture, and then find the intersection of the projected spherical wave (the circle) with the active aperture as shown in Fig. 2.
Figure 2. Intersection of spherical waves from the field point by the aperture, when the field point is projected onto the plane of the aperture
Rewriting the integral into polar coordinates gives
^2 d15(t - R) " " c rdrdG,
h(r,t) =H
(3)
ex di
2nR
where r is the radius of the projected circle and R is the distance from the field point to the aperture given by R2 = r2 + zp2. Here zp is the field point height above the x-y plane of the Aperture. The projected distancesd1, d2 are determined by the aperture and are the distance closest to and furthest away from the aperture, and ei, e2 are the corresponding angles for a given time (see Fig. 3).
Figure 3.
Definition of distances and angles in the aperture plan for evaluating the Rayleigh integral
Introducing the substitution 2RdR = 2rdr and finally using the substitution t' = R/c gives
e2 12
h(r1, t) = \S(t -1')dt'de.
(4)
1 n
For a given time instance the contribution along the arc is constant and the integral gives h(r„ t) = c. (5)
The angles 61 and 62 are determined by the intersection of the aperture and the projected spherical wave, and the spatial impulse response is, thus, solely determined by these intersections, when no apodization of the aperture is used. The response can therefore be evaluated by keeping track of the intersections as a function of time [16].
2. CALCULATION FOR ARRAY TRANSDUCERS
An array transducer, generally two dimensional, is consisting of a number of rectangular elements. Each element of the array is a rectangle that restricted with four lines and the distance between the elements can be the same or different (in x and y directions). In all cases the coordination of the rectangular elements can be determined.
Figure (4) illustrates the first element and the position of observation point and its image on (x-y) plane.
Figure 4. One element transducer, the circle is the observation point and the cross is its image
on transducer plane
By assuming a and b, the length and width of the rectangle respectively on x and y direction, rectangle's coordinates would be determined:
{y = 0;0 < x < a} {y = b;0 < x < a},
(6)
{x = 0;0 < y < b} {x = a;0 < y < b}.
The coordinates of the observation point and its image are
O Уo, zo), O':^ Уo, ^0). (7)
In this condition, it is convenient to find intersection points of the r radius circles and the transducer on x-y plane (Fig. 5).
Y
S
0
__l_I_I_I_1_I_I_I_
0 5 10 15 20 25 30 35
a(rinn) K
Figure 5. One r radius circle and its intersection with transducer
The algorithm that is chosen to find the intersection points is dependent on the position of the image point O in any regions of figure (6).
Figure 6. Regions around one element of transducer; the image of observation point would be
in one of these regions
The minimum and maximum of r (radius of image circle) are calculated with different algorithms if the image point would be in any of 16 regions. These r magnitudes are the upper and lower limits of circles' radiuses that have at least one intersection point with aperture and it would be possible to calculate the moment of receiving the first and last part of the wave.
As shown in figure (5) after determination of intersection points to find effective arcs, the arcs inside the aperture's limitation, (ticked on figure) using of different methods is possible. Here, in a particular case, the method of calculation of arcs' angles is described. By assuming a new coordination system centered at O, the arc is being calculated like this:
0 = tan
f
y A - yo
\
(
V XA XO ' J
02 = tan
yB - yo
\
V XB xo '
arc = 0 - 02.
(8)
In the case of multi-element transducer there are two different methods to calculate the response of all elements. In the first method (one observation point and «-element transducer), the impulse response of each element of transducer at the observation point is calculated. It is sufficient to repeat the above algorithm to find intersection points and the effective arcs between them and finally the impulse response. This method is complicated and time consuming. In the second method that is used in SimATra, the case "one observation point and «-element transducer" is changed to the case "one element and n observation point". The other n-1 observation points are the result of displacement of the main point related to the transducer elements. The number of imaginary points is the same as the other elements of transducer. The total impulse response can be calculated by summing the response of all points in any instance.
3. FIELD CALCULATION AND SIMULATION
3.1. Pulse Field
As it is shown in equation (2), to calculate the pulsed field for each point, medium density should be multiplied into the time domain convolution of derivative of excitation function and impulse response to determine the pressure for the specific point.
3.2. The Continuous Wave Field
The major characteristics of linear systems have been utilized in this software to calculate the continuous wave field, as the transducer is assumed as a linear system with a determined impulse response. Therefore, it would be easy to determine the response of this system to a sinusoidal input in a given frequency by calculating of its transfer function in the frequency domain. This is shown in the figure 7:
Figure 7. The scheme of calculation mono-frequency response from impulse response
In figure 7, (0 is sinusoidal input frequency of the system. This equation is describing that
it would be convenient to calculate the continuous wave response using the amplitude and phase of transfer function at the frequency of excitation. Utilizing this method would increase the performance of the software as it is based on the impulse response calculation which is already obtained. In addition, it would be easy to calculate the continuous wave pressure field for any given sinusoidal input or even a combined input of several sinusoidal vectors with variant amplitudes.
4. EXAMPLES
4.1. Amplitude and Waveform Comparison
In order to validate performance and calculations of the software more precisely, in terms of amplitude and impulse response, some of the models already detailed in other papers have been used. In the first example (figure 8), the impulse response is calculated for a single element rectangular transducer of 4 x 5 mm in different observation points in front of the transducer surface and distanced by 10 mm from it. As figure 8-a illustrates this points are in 1 mm distance from each other in the x direction. Sound velocity is assumed to be 1540 m/s. Figure 8-b illustrates the impulse response relevant to these points plotted in [17] and 8-c is the same response calculated and drawn by SimATra.
(b)
Figure 8. (a) one element transducer with dimension 4x 5 mm and observation points are located at 10 mm from the surface, (b) impulse response at observation points plotted in [17], (c) SimATra response that is exactly the same as (b). There is no difference
4.2. Focusing and Steering
By giving proper delays to the stimulating pulse applied on array transducer elements (linear or 2D), it is possible to direct the consequent of the waves distributed in the space in front of the transducer to a specified direction and focus it on a particular point. The delays pertaining to different elements are separately calculated according to the position of the canonical point.
For example, having the array transducer of figure 9-a in mind, if focusing of the waves resulted from the elements on (8.5, 4, 45) and (18.5, 4, 65) mm is desired, the required delays has been calculated by SimATra after giving the coordinates of the canonical points. As expected by applying these delays, the results are illustrated in figure 9-b and 9-c.
(b)
(a)
(c)
Non-naJiMd Pmnn
Figure 9. (a) A transducer with nine 1 x 8 mm elements and two focal points. (b and c) Total wave at the focal point before (upper figure) and after (lower figure) applying the delays; as expected after applying them, the waves from different elements arrive at the
focal point simultaneously
4.3. Apodization
In SimATra, it is also possible to apply Apodization on the linear transducer. For example when stair function is chosen as apodization function, for the 9-element transducer of figure 9-a, weighing coefficients are illustrated on figure 10-a. Figures 10-b is the pressure magnitude for a continuous wave excitation on the plane normal to the transducer surface without apodization and figure 10-c illustrates the result after applying apodization. As expected with apodization the side lobes are smaller and the main lobe is bigger.
Stair Apodizalion Faction
(a)
(b)
(c)
Figure 10. (a) Stair weighing coefficients for the 9-element transducer of figure 9-a, (b) pressure on the plane normal to the transducer surface without apodization, (c) the result after applying apodization
4.4. Ultrasim and SimATra Comparison in Near-Field
Figure 11 shows the results of simulating the near field (without focusing) for the 2D array defined in table 1 with a continuous-wave stimulation. The observation plane is assumed to be normal to the center of the transducer. For comparison, the results of the same simulation using Ultrasim are also shown. Both figures are similar and the little difference between the qualities is because of the lower number of observation points in SimATra.
Table 1. The parameters of 2D array transducer
Number of elements 8 x 8
Dimension of elements 0.4 x 0.4 mm
Distance between elements 0 mm
Sound speed in medium 1540 m/s
Medium density 1000 Kg / m3
Frequency 3.5 MHz
Figure 11. The results of simulating for the 2D array defined in table 1. (a) Using SimATra, (b) using Ultrasim
4.5. Field II and SimATra Comparison
In field II it is possible to increase and decrease the number of divisions of transducer elements. In calculation near field if the number of such divisions is low, it could not be done with good approximation, but in far-field calculation there is a good approximation. In SimATra there is no need to such settings and in all observation point, field is calculated with no approximation. Above mentioned is shown in Figure 12. The figures on the right and those on the left respectively show the pressure in observation points distanced by 3 mm (near-field) and 50 mm (far-field). As shown in figures the more the number of division the more accuracy, while the figures of SimATra have good accuracy and precision in any distance.
Figure 12. Results of Field II in near-field, 3mm (right hand) and far-field, 50mm (left hand).
Rows 1-4 belong to Field II and their divisions are 2 x 2, 5 x 5, 15 x 15, 25 x 25, 50 x 50 respectively. Increasing accuracy in the right column is clear. (i) & (j) are related to SimATra
CONCLUSIONS
A method for calculation of ultrasound field from continuous wave and pulse excitation for array transducers is presented. The method has been implemented as a simulator software. By defining array transducer parameters such as elements dimensions, their distances and environment features, this software is capable of calculation of ultrasound pressure field and plotting the results in various methods. A special method is used that computes field of arcs with the same distance on the transducer from the observation point for calculating of impulse response, then considering the excitation function, with a simple convolution, pressure field can be calculated. In the case of multi-element transducer that is used in SimATra, the case "one observation point and «-element transducer" is changed to the case "one element and n observation points". The other n-1 observation points are the result of displacement of the main point related to the transducer elements. The number of imaginary points is the same as the other elements of transducer. The total impulse response can be calculated by summing the response of all points in any instance. In SimATra the transducer is assumed as a linear system with a determined impulse response. Therefore, it would be easy to determine the response of this system to a sinusoidal input in a given frequency by calculating of its transfer function in the frequency domain.
In spite of almost all other softwares that use many approximations in their calculations to reduce the computational time and complication, the method used in SimATra (linear system) applies no approximation. According to technological development of computers' hardware and by using the method described in this paper, the application of approximations seems no necessary, specially in near-field calculations.
In order to validate performance and calculations of the software more precisely, with some examples the accuracy and preciseness of SimATra were illustrated in the last section but comparison of the computational time needs more evaluations and examples. SimATra can be used for optimizing array parameters in designing and testing stages and also for education and research purposes.
REFERENCES
1. IEEE Ultrasonics, Ferroelectrics and Frequency Control web site for software package and utilities for ultrasound research: http://www.ieeeuffc.org.
2. J. A. Jensen. Users' guide for the Field II program, August 2001.
3. http://www.ifi.uio.no/~ultrasim.
4. J. A. Jensen, N. B Svendsen. Calculation of Pressure Fields from Arbitrarily Shaped Apodized, and Excited Ultrasound Transducers. IEEE Trans. Ultrason., Ferroelec., Freq., contr., 1992, 39, 262-267.
5. J. A. Jensen, S. I Nikolov. Fast Simulation of Ultrasound Images. IEEE International Ultrasonic Symposium, Puerto Rico, 2000.
6. J. A. Jensen. A New Approach to Calculating Spatial Impulse Responses. IEEE International Ultrasonic Symposium, Toronto, Canada, 1997.
7. J. A. Jensen. Simulating Arbitrary-Geometry Ultrasound Transducers Using Triangles. IEEE International Ultrasonic Symposium, Antonio, Texas, 1996.
8. J. A. Jensen. A program for simulating Ultrasound Systems", Published in Medical & Biological Engineering & Computing, Volume 34, Supplement 1, Part1, 1996, pp.351353.
9. Severre Holm, Frode Teigen, Lars Odegaard, Vebjorn Berre, Jan Ove Erstad, Kapila Epasinghe. Ultrasim User's Manual, ver. 2.1. Program of Simulation of Ultrasonic Fields. April, 1998.
10. Severe Holm. Simulation of Acoustic Fields from Medical Ultrasound Transducers of Arbitrary Shape. Nordic Symposium in Physical Acoustics, Norway, January, 1995.
11. L. Odegard, S. Holm, F. Teingen, T. Kleveland. Acoustic Field Simulation for Arbitrarily Shaped Transducers in a Stratified Medium. University of Oslo.
12. Hamid Behnam, Elham Mokhtar Tajvidi. Simulation of acoustic fields from 2D-array transducers with continuous wave excitation in homogeneous medium. Electronic Journal "Technical Acoustics", http://www.ejta.org, 2007, 14.
13. A. Austeng, S. Holm. Sparse 2D Arrays for Real-Time 3D Ultrasound. University of Oslo.
14. S. Ohtsuki. Ring function method for calculating near field of sound source. Bull. Of the Tokyo Inst. of Tech., 1974, N°123, 23-27.
15. A. D. Pierce. An Introduction to Physical Principles and Applications. Acoustical Society of America, New York, 1989.
16. J. A. Jensen. Linear description of ultrasound imaging systems. Notes for the International Summer School on Advanced Ultrasound Imaging. Technical University of Denmark. Release 1.01, June 29, 2001.
17. J. A. Jensen. A new Calculation Procedure for Spatial Impulse Responses in Ultrasound. Journal of the Acoustical Society of America, vol. 105, October, 1999, pp. 3266-3274.
