A Fast Method to Calculate the Spatial Impulse Response for 1-D Linear Ultrasonic Phased Array Transducers

A method is developed to accurately determine the spatial impulse response at the specifically discretized observation points in the radiated field of 1-D linear ultrasonic phased array transducers with great efficiency. In contrast, the previously adopted solutions only optimize the calculation procedure for a single rectangular transducer and required approximation considerations or nonlinear calculation. In this research, an algorithm that follows an alternative approach to expedite the calculation of the spatial impulse response of a rectangular linear array is presented. The key assumption for this algorithm is that the transducer apertures are identical and linearly distributed on an infinite rigid plane baffled with the same pitch. Two points in the observation field, which have the same position relative to two transducer apertures, share the same spatial impulse response that contributed from corresponding transducer, respectively. The observation field is discretized specifically to meet the relationship of equality. The analytical expressions of the proposed algorithm, based on the specific selection of the observation points, are derived to remove redundant calculations. In order to measure the proposed methodology, the simulation results obtained from the proposed method and the classical summation method are compared. The outcomes demonstrate that the proposed strategy can speed up the calculation procedure since it accelerates the speed-up ratio which relies upon the number of discrete points and the number of the array transducers. This development will be valuable in the development of advanced and faster linear ultrasonic phased array systems.


Introduction
Ultrasonic phased array transducers are attracting attention across the board with respect to applications in both non-destructive evaluation (NDE) and clinical diagnostics. Therefore, acquiring more adequate information and enhancing inspection efficiency is significant. Ultrasound field simulation is an essential tool for probe design and optimization. The accuracy and data processing speeds are the major challenges involved in the probe design process. Various semi-analytical methods have been developed to address this problem, such as Field II [1,2], FOCUS [3,4], DREAM [5,6], Ultrasim [7], k-wave [8], and the angular spectrum technique [9][10][11][12]. Additionally, the spatial impulse response (SIR) is an excellent method to predict the transient ultrasound field radiated by an ultrasonic phased array probe in an isotropic, homogeneous, and non-dissipative medium [13]. The SIR method was first proposed by Oberhettinger [14], and then further developed by Tupholme et al. [15] and the specific characteristic of the array, i.e., the shapes of the transducers in an array are assumed to be the same, meaning that the spatial relationship between the field points and different transducers are also the same, only if the relative position of the field points are identical in the local coordinate system of the transducers. With this consideration, the computational cost can be reduced if the radiated field points are elaborately designed to meet the specific characteristics of array transducers.
The main content is outlined as follows: first, the classical theory of spatial impulse response is summarized; second, the complexity of the spatial impulse response of a rectangular aperture is discussed and the relationship between the aperture shape and size, and the discontinuity possibilities, are analyzed; third, the improved computational method of the spatial impulse response for 1-D linear array transducers is discussed and, later, several computing simulations are carried out to compare the time required for the proposed method and the classical summation method.

Spatial Impulse Response Function
A rectangular transducer, with length 2b and width 2a (2b ≥ 2a), is located on the plane z = 0, surrounded by an infinite rigid baffle, as shown in Figure 1. The radiated space is a homogeneous lossless fluid medium, with density ρ and speed of sound c. The center of the aperture is located at the origin of the coordinate system, with length along the y axis and width along the x axis. The transient pressure p(x,y,z;t) radiated from a planar transducer at point P(x,y,z), and the time instant t, can be written as the convolution of the first-order derivative of the normal velocity of the aperture surface, u(t) with a spatial impulse response function h(x,y,z;t): p(x, y, z; t) = ρ ∂u(t) ∂t * h(x, y, z; t), where the asterisk (*) denotes the convolution operation, and the spatial impulse response function h(x,y,z;t) is derived from the Huygens' principle [17] as: h(x, y, z; t) = S δ(t − R/c) 2πR ds, where S denotes the active radiating surface of the aperture, δ(t) is the Kronecker function, and R is the distance between the infinitesimal ds and the field point P(x,y,z).
Sensors 2016, 16, 1873 3 of 16 transducers are also the same, only if the relative position of the field points are identical in the local coordinate system of the transducers. With this consideration, the computational cost can be reduced if the radiated field points are elaborately designed to meet the specific characteristics of array transducers.
The main content is outlined as follows: first, the classical theory of spatial impulse response is summarized; second, the complexity of the spatial impulse response of a rectangular aperture is discussed and the relationship between the aperture shape and size, and the discontinuity possibilities, are analyzed; third, the improved computational method of the spatial impulse response for 1-D linear array transducers is discussed and, later, several computing simulations are carried out to compare the time required for the proposed method and the classical summation method.

Spatial Impulse Response Function
A rectangular transducer, with length 2b and width 2a (2b ≥ 2a), is located on the plane z = 0, surrounded by an infinite rigid baffle, as shown in Figure 1. The radiated space is a homogeneous lossless fluid medium, with density ρ and speed of sound c. The center of the aperture is located at the origin of the coordinate system, with length along the y axis and width along the x axis. The transient pressure p(x,y,z;t) radiated from a planar transducer at point P(x,y,z), and the time instant t, can be written as the convolution of the first-order derivative of the normal velocity of the aperture surface, u(t) with a spatial impulse response function h(x,y,z;t): where the asterisk (*) denotes the convolution operation, and the spatial impulse response function h(x,y,z;t) is derived from the Huygens' principle [17] as: where S denotes the active radiating surface of the aperture, δ(t) is the Kronecker function, and R is the distance between the infinitesimal ds and the field point P(x,y,z). From Equation (2), the spatial impulse response at the time instant t only varies with the position of the field point for a specific aperture. The surface integral in Equation (2) can be expressed instead by using polar coordinates with the origin at the field point projection on the aperture plane P′(x,y,0): where Θ1 and Θ2 are the minimum and maximum angle components, and r2 and r1 are the minimum and maximum radius components. From Equation (2), the spatial impulse response at the time instant t only varies with the position of the field point for a specific aperture. The surface integral in Equation (2) can be expressed instead by using polar coordinates with the origin at the field point projection on the aperture plane P (x,y,0): where Θ 1 and Θ 2 are the minimum and maximum angle components, and r 2 and r 1 are the minimum and maximum radius components. Deriving from the relationship R 2 = r 2 + z 2 , the equation 2RdR = 2rdr can be obtained, and then substituting the equation with τ = R/c into Equation (3), the spatial impulse response function can be simplified as the integral of angle at time instant t multiplying with a constant term: where t min and t max determine the time range of integration, and Θ 1 (t) and Θ 2 (t) are determined by the intersection of the transducer and the projection spherical wave centered at the field projection as shown in Figure 2 with radius R = ct at the time instant t. The intersection between the spherical wave and the plane z = 0 is a projection circle centered at P (x,y,0) with radius r (t) = (ct) 2 − z 2 , when t > z/c. The integral Equation (4) can be expressed as superimposition of the central angles of the arcs intersected between the projection circle and the transducer aperture. If more than one arc exists, the spatial impulse response can be expressed as the summation of the central angles of all the intersected arcs inside the transducer aperture: where k is the index of the intersected arcs, Θ k1 (x,y,z;t) and Θ k2 (x,y,z;t) denote the start and the end angle of intersected arc seperately, Θ k (x,y,z;t) denotes the central angle of the k-th arc, and Θ Σ (x,y,z;t) denotes the summation of all of the central angles of the intersected arcs.
where tmin and tmax determine the time range of integration, and Θ1(t) and Θ2(t) are determined by the intersection of the transducer and the projection spherical wave centered at the field projection as shown in Figure 2 with radius R = ct at the time instant t. The intersection between the spherical wave and the plane z = 0 is a projection circle centered at P′(x,y,0) with radius ( ) = ( ) − , when t > z/c. The integral Equation (4) can be expressed as superimposition of the central angles of the arcs intersected between the projection circle and the transducer aperture. If more than one arc exists, the spatial impulse response can be expressed as the summation of the central angles of all the intersected arcs inside the transducer aperture: where k is the index of the intersected arcs, Θk1(x,y,z;t) and Θk2(x,y,z;t) denote the start and the end angle of intersected arc seperately, Θk(x,y,z;t) denotes the central angle of the k-th arc, and ΘΣ(x,y,z;t) denotes the summation of all of the central angles of the intersected arcs.

Spatial Impulse Response Function of a Rectangular Piston
Based on the previous derivation, the spatial impulse response function is expressed as a summation of the central angles of the intersected arcs. The calculation of the central angles become the key point for computing the spatial impulse response. Consider the Cartesian coordinates in the transducer aperture plane, as shown in Figure 2. Due to the geometrical symmetry of rectangular aperture, the projection points can be transformed to the first quadrant x ≥ 0, y ≥ 0:

Spatial Impulse Response Function of a Rectangular Piston
Based on the previous derivation, the spatial impulse response function is expressed as a summation of the central angles of the intersected arcs. The calculation of the central angles become the key point for computing the spatial impulse response. Consider the Cartesian coordinates in the transducer aperture plane, as shown in Figure 2. Due to the geometrical symmetry of rectangular aperture, the projection points can be transformed to the first quadrant x ≥ 0, y ≥ 0: The points in the first quadrant will be considered. The four vertices are named as o i (i = 1-4). The four edges are named as s i (i = 1-4). The distance between the projection point P and the vertices and the edges of the rectangle are defined as r i (i = 1-4) and d i (i = 1-4), separately: Define td i and tr i as the time constants when the projection circle expands encountering edges and the vertices: Define t 0 = z/c as the spherical wave encounter with the aperture plane z = 0.
To reduce the calculation expressions, the angles α i (i = 1-4) are defined as critical angles of the intersected arc, with radius r(t) intersected with the i-th edge: where r(t) denotes the radius that varies with t for the projection point at (x,y,z). The spatial coordinate symbols are omitted here. It should be noted that Equations (7) and (8) hold when t ≥ t 0 , and Equation (10) holds when the arc intersects with corresponding edge, d i ≤ r(t). The complementary angle of α i is defined as: Due to the complexity of the rectangular aperture geometry, expression changes may occur when the intersected arcs step over the vertices or the edges of the rectangle. The expression changes then cause discontinuities in the temporal slope of Θ Σ (x,y,z;t) and finally affect the temporal slope of p(x,y,z;t). San Emeterio and Ullate [18] presented a detailed solution about this problem, and the methods proposed by them will be implemented as the classical method in the following sections. The calculation procedure varies with the position of the projection points.
As shown in Figure 3, the first quadrant of the aperture plane is divided into four regions by the lines x = a and y = b: I (x ≥ a, y ≤ b), II (x ≥ a, y ≥ b), III (x ≤ a, y ≥ b) and IV (x ≤ a, y ≤ b). The points in the first quadrant will be considered. The four vertices are named as oi (i = 1-4). The four edges are named as si (i = 1-4). The distance between the projection point P′ and the vertices and the edges of the rectangle are defined as ri (i = 1-4) and di (i = 1-4), separately: Define tdi and tri as the time constants when the projection circle expands encountering edges and the vertices: Define t0 = z/c as the spherical wave encounter with the aperture plane z = 0.
To reduce the calculation expressions, the angles αi (i = 1-4) are defined as critical angles of the intersected arc, with radius r(t) intersected with the i-th edge: where r(t) denotes the radius that varies with t for the projection point at (x,y,z). The spatial coordinate symbols are omitted here. It should be noted that Equations (7) and (8) hold when t ≥ t0, and Equation (10) holds when the arc intersects with corresponding edge, di ≤ r(t). The complementary angle of αi is defined as: Due to the complexity of the rectangular aperture geometry, expression changes may occur when the intersected arcs step over the vertices or the edges of the rectangle. The expression changes then cause discontinuities in the temporal slope of ΘΣ(x,y,z;t) and finally affect the temporal slope of p(x,y,z;t). San Emeterio and Ullate [18] presented a detailed solution about this problem, and the methods proposed by them will be implemented as the classical method in the following sections. The calculation procedure varies with the position of the projection points.
As shown in Figure 3, the first quadrant of the aperture plane is divided into four regions by the lines x = a and y = b: Figure 3. General geometry of the divided regions and sub-regions of the aperture plane.
The curves l1, l2, and l3 denote r1 = d2, r2 = r4, and r4 = d3 separately, which are expressed as: The curves l 1 , l 2 , and l 3 denote r 1 = d 2 , r 2 = r 4 , and r 4 = d 3 separately, which are expressed as: The dashed lines are illustrated as the extension of Equation (12). The curves l 1 and x = a intersect at the point q 1 (a,b−2a). The curves l 2 and x = a intersect at the point q 2 (a,a 2 /b). The curves l 1 , l 3 , and y = 0 intersect at the point q 3 (b 2 /4a,0). The curves l 2 and y = b intersect at the point q 4 (b 2 /a,b).
The four regions are then divided into several sub-regions by the curves. The projection points located in the same sub-region share the same calculation procedure for the central angle summation Θ Σ (x,y,z;t).
As presented in Figure 3, the number of sub-regions in region I relies on the position of the three intersection points q 1 , q 2 , and q 3 , which has direct relation to the calculation procedure of the central angle of intersected arcs. It can be derived that the position of the three intersection points are related to the ratio of the width and height of the rectangle. The relationships between the number of sub-regions in region I and the ratio of the width and height are presented in Table 1. Table 1. Relationship between the rectangular geometry and the possibility of discontinuities.

Relationship of Width and Height
Sub-Regions The order of time instant precedence and the central angle summation Θ Σ (t) for the projection points located in the sub-regions of region I are presented in Table 2. Typically, 1-D array transducers employ rectangular transducers with a < √ 2 − 1 b, which means there are five sub-regions in the region I. Table 2. Central angle summation of the intersected arc for projection points located in sub-regions I 1 , I 2 , I 3 , I 4 , and I 5 .
Region II is always divided into two sub-regions by the line l 2 , which are denoted as II 1 and II 2 . The central angle summation calculation expressions for the projection points located in the sub-regions of region II are shown in Table 3. Table 3. Central angle summation of the intersected arc for projection points located in sub-regions II 1 and II 2 .
Region III maintains itself without sub-regions. There is only one possibility of the discontinuity order of the time instant for the projection points located in region III. The computation procedure for projection points in region III is: Since the possibilities of the discontinuity order of the time instant are complicated when the projection points are located inside the rectangle, region IV, the calculation procedure is far more complicated when compared to the procedure for points outside the rectangle [18]. The alternative solution for this condition is implemented in this paper. The projection point P is used to divide the rectangle aperture into four sub-rectangles, as shown in Figure 4. The shape of the rectangle, the position of the projection point P and the circle shown in Figure 4 are selected so as to illustrate the condition that the projection circle intersects with all of the four edges. Region III maintains itself without sub-regions. There is only one possibility of the discontinuity order of the time instant for the projection points located in region III. The computation procedure for projection points in region III is: , Since the possibilities of the discontinuity order of the time instant are complicated when the projection points are located inside the rectangle, region IV, the calculation procedure is far more complicated when compared to the procedure for points outside the rectangle [18]. The alternative solution for this condition is implemented in this paper. The projection point P′ is used to divide the rectangle aperture into four sub-rectangles, as shown in Figure 4. The shape of the rectangle, the position of the projection point P′ and the circle shown in Figure 4 are selected so as to illustrate the condition that the projection circle intersects with all of the four edges. The final central angle summation can be expressed as the sum of central angles created by intersection of the projection circle and each sub-rectangle: where the symbols x, y, z, and t are omitted for simplicity.  The final central angle summation can be expressed as the sum of central angles created by intersection of the projection circle and each sub-rectangle: where the symbols x, y, z, and t are omitted for simplicity.
where α i is defined from Equation (10) when the arc intersects the i-th edge. In addition, α i equals 0 and α i equals π/2 when r(t) ≤ d i .
The expressions introduced are direct, although, they are complicated due to the complexity of the rectangular geometry. However, they can be easily implemented.
Moreover, since the intersected arcs only exist when the wave sphere intersects with the rectangular aperture, the spatial impulse response function is continuous nonzero from the first time instant, t min , the projection circle intersects the rectangle or is located within the rectangle to the last time instant, t max , the radius exceeds r 3 . The range of for each region and sub-region is shown in Table 4. The values beyond this range are zeros and summation procedure of these zero values can be avoided.

Spatial Impulse Response for Linear Array Transducers
Previously, the classical methodology for computing the transient pressure at field point (x,y,z) radiated from a linear array transducer with N planar rectangular transducers was a simple summation of the pressure from each single element [25]: where i is the element index number, from 0 to N − 1, h i (x,y,z;t) is the spatial impulse response at field point (x,y,z) radiated from the i-th element, u i (t) is the aperture surface vibration velocity of the i-th element, and a i is the weighting level of the i-th element. T i is the delayed time of the i-th element. Generally, the aperture surface vibration velocities are the same, so Equation (16) can be expressed as: where h Σ (x,y,z;t) is the total spatial impulse response at field point P: The general geometry of linear phased array transducers and the discretized field points are shown in Figure 5. where hΣ(x,y,z;t) is the total spatial impulse response at field point P: The general geometry of linear phased array transducers and the discretized field points are shown in Figure 5.  Figure 5. Geometry of the discretized field points and auxiliary points for linear phased array transducers. The field points are discretized with N x , N y , and N z along with the x, y, and z axes. The x component of each point meet the relationship: where k is the x component index, and d is the pitch between two neighboring transducers. Although the y and z components do not have a mandatory rule, they are typically evenly distributed. The spatial impulse response at field points (x k+j ,y,z) and (x k ,y,z) will meet the relationship: where j is an arbitrary integer. This can be easily understood under the assumption that the geometry shape and the surface vibration velocity of the transducer elements are identical, because if two different field points have the same relative distance from two different transducers, they share the same spatial impulse response function contributed from the corresponding transducers. This means that repeated computations exist if the traditional summation method is utilized. The calculation process can be optimized for the field points which meet the distribution requirement introduced previously. The spatial impulse response h i (x k ,y,z;t) at point (x k ,y,z) contributed from transducer 1 to N−1 can be derived from the expression: which means the final result can be derived only from computing the result contributed from transducer 0. The redundant calculation procedures can be avoided. From Equation (21), if k < N − 1, the condition that k − i < 0 may occur, which means that the right part of Equation (21) may not exist since the x indices of the field points are not less than zero. In order to fix this problem, N − 1 auxiliary points are created in the x direction for every y and z: Thus, the total spatial impulse response function at field point (x k ,y,z) can be expressed as: Finally, the pressure at field point (x k ,y,z) can be obtained from Equation (17) by only considering the field contributed from transducer 0.
Generally, a field point gap in the x direction set equal to d lowers the obtained precision. However, if higher accuracy needs to be obtained, the proposed method can be easily modified. For more accuracy, the field point gap in the x direction can be set as d/M, where M > 1, the field points can be divided into M groups: where x m,0 is the start field point on direction of x for group m. Consequently, N − 1 auxiliary points are required to be created on direction of x for each group of field points: The spatial impulse response at field points (x m,k+j ,y,z) and (x m,k ,y,z) will meet the relationship: Correspondingly, Equation (21) can be expressed as: h i (x m,k , y, z; t) = h 0 (x m,k−i , y, z; t), for i = 1, 2, · · ·, N − 1.
The spatial impulse response function at field point (x m,i , y, z) can be expressed as: The computation procedure of Equation (28) is illustrated in Figure 6 for m = 0, . . . , M − 1 loops.
The computation procedure of Equation (28) is illustrated in Figure 6 for m = 0,…, M − 1 loops. The block shown in Figure 6 denotes one step of the calculation of spatial impulse response function for a field point at the time constant t contributed from the one single transducer. Since the functions are related to the x variable, the symbols y, z, and t are omitted and the symbol xm,k is simplified as xk for simplicity of illustration. Ti denotes the time shift index for summation. It can be seen from Figure 6 that the calculation of the spatial impulse response involves the same transducer, which is called the basic transducer. Moreover, the basic transducer can be any one from the array transducers because the relationship shown in Equation (20) implies that the identical characteristic holds only if the field points are distributed as required. The block shown in Figure 6 denotes one step of the calculation of spatial impulse response function for a field point at the time constant t contributed from the one single transducer. Since the functions are related to the x variable, the symbols y, z, and t are omitted and the symbol x m,k is simplified as x k for simplicity of illustration. T i denotes the time shift index for summation. It can be seen from Figure 6 that the calculation of the spatial impulse response involves the same transducer, which is called the basic transducer. Moreover, the basic transducer can be any one from the array transducers because the relationship shown in Equation (20) implies that the identical characteristic holds only if the field points are distributed as required.

Optimization of the Computation Procedure
Furthermore, from Figure 6, the computation procedure contains noteworthy redundancy. For example, the calculation of h 0 (x 2−N ), . . . , h 0 (x −1 ) and h 0 (x 0 ) are used to calculate h Σ (x 1 ) again. Only h 0 (x 1 ) needs to be calculated if the results obtained during the last step can be stored in a lookup memory. Consequently, the computational cost can be reduced by reusing the results which have been calculated for the last field point. The procedure of the optimized method is illustrated in Figure 7. The block shown in Figure 6 denotes one step of the calculation of spatial impulse response function for a field point at the time constant t contributed from the one single transducer. Since the functions are related to the x variable, the symbols y, z, and t are omitted and the symbol xm,k is simplified as xk for simplicity of illustration. Ti denotes the time shift index for summation. It can be seen from Figure 6 that the calculation of the spatial impulse response involves the same transducer, which is called the basic transducer. Moreover, the basic transducer can be any one from the array transducers because the relationship shown in Equation (20) implies that the identical characteristic holds only if the field points are distributed as required.

Optimization of the Computation Procedure
Furthermore, from Figure 6, the computation procedure contains noteworthy redundancy. For example, the calculation of h0(x2−N), …, h0(x−1) and h0(x0) are used to calculate hΣ(x1) again. Only h0(x1) needs to be calculated if the results obtained during the last step can be stored in a lookup memory. Consequently, the computational cost can be reduced by reusing the results which have been calculated for the last field point. The procedure of the optimized method is illustrated in Figure 7.
h 0 (x 0 )  The classical method for calculating the spatial impulse response for linear phased array transducers implements a simple summation of the radiated field from each transducer and suffers N x × N y × N z × N times the number of calculations of the spatial impulse response functions and N x × N y × N z × (N − 1) times the number of calculations of summations. The proposed method reduces the calculation to only (N x + N − 1) × N y × N z times the number of calculations of the spatial impulse response functions and N x × N y × N z × (N − 1) times the number of calculations of summations. The estimated speed-up ratio of computation times between the classical method and the proposed method can be expressed as: From Equation (29), the speed-up ratio is dependent on two factors: the number of discretized field points in the direction of the x axis and the transducer number. The speed-up ratio also increases along with both of the two factors.
Since the anti-trigonometric calculation process costs much more time than the summation process, this can greatly reduce the impulse response calculation process without any nonlinear approximation. Examples will be given to demonstrate how much the efficiency the proposed method can be improved. It is expected that the efficiency can be improved even more when the transducer number increases.

Comparison between the Proposed Method and the Classical Approach
In order to validate the speed-up efficiency of this new method, a series of simulations are carried out. One-dimensional linear arrays with 16, 32, 64, 128, and 256 elements are used for comparison. The arrays share the same pitch of 0.15 mm. The dimension of each rectangular transducer is the same, 0.14 mm in width and 14 mm in length. The arrays are located on the 2-D plane defined by (x, 0, z), centered at (0, 0) [mm] and linearly distributed along the x axis.
The field points are located on the 2-D plane defined by (x, 0, z). Four rectangular fields with different widths in the x direction and the same length, 60 mm, in the z direction. The points are discretized with a gap, dx = dz = 0.15 mm, in both the x and z directions. All of the fields start from z 0 = 0 mm. The discretized numbers, N x , on direction of x are, 65, 129, 257, and 513, which form the widths of 9.6 mm, 19.2 mm, 38.55 mm, and 76.8 mm of the rectangular fields. The center of the width is at x = 0 mm.
The time sampling frequency is 100 MHz. The time range for calculation is from 0 to 80 µs, which involves 8001 time instants. The simulations were carried out on a Windows © 7 64-bit PC with i7 4-cores 3.4 GHz CPU and 32 GB DDR3 1600 MHz memory. The proposed method and the classical method were both programed and compiled using the C++ language. The programs were established without using multithreading optimization. Both programs implemented the previously mentioned optimization algorithm to avoid redundant addition procedure of zero values.
The time cost for calculating the spatial impulse response at the 2-D plane field points radiated from the different transducers, with 16, 32, 64, 128, and 256 elements, are presented in Figure 8. The solid lines and dashed lines, with different marks, represent the results carried out using the classical method and the proposed method, respectively, for different widths of field points. Both of the axes are logarithmic coordinates for clear illustration. It can be perceived from Figure 8 that the time cost increases along with the number of the array transducers, and also along with the number of the discretized points in the direction of x. Additionally, compared with the results carried out using the classical method, the computation time for the proposed method has reduced. For instance, the time reduced from about 110 s to about 6 s for 256 array transducers and 513 points in the x direction.  The speed-up ratio is calculated through dividing the computation time cost using the classical method by that using the proposed method for each simulation pair of the same transducer number and discretized points in the direction of x. The results are given in Figure 9. The speed-up ratio is both dependent on the number of the array transducers and the number of discretized point in the direction of x. Although, the speed-up ratio does not increase constantly along the two factors, it can be distinguished that, with more linear distributed transducers or a higher number of discretized points in the x direction, the speed-up ratio increases. The increment of the speed-up ratio reduces as the transducer number increases. This is due to the fact that, in the proposed method, auxiliary points are created with the same number as the transducer number in the x direction, which can be seen from Figure 5. This can also increase the cost calculation time if the number of auxiliary points are comparable to, or greater than, the number of the discretized points in the x direction. Although the speed-up efficiency will rebate under this condition, the proposed method could speed up the calculation procedure with significant ratios when the transducer number and the number of the discretized points in the x direction increase. For instance, the speed-up ratio reaches about 18.18 when the number of the array transducers is 256 and the number of discretized points in the x direction is 513. The speed-up ratio is calculated through dividing the computation time cost using the classical method by that using the proposed method for each simulation pair of the same transducer number and discretized points in the direction of x. The results are given in Figure 9. The speed-up ratio is both dependent on the number of the array transducers and the number of discretized point in the direction of x. Although, the speed-up ratio does not increase constantly along the two factors, it can be distinguished that, with more linear distributed transducers or a higher number of discretized points in the x direction, the speed-up ratio increases. The increment of the speed-up ratio reduces as the transducer number increases. This is due to the fact that, in the proposed method, auxiliary points are created with the same number as the transducer number in the x direction, which can be seen from Figure 5. This can also increase the cost calculation time if the number of auxiliary points are comparable to, or greater than, the number of the discretized points in the x direction. Although the speed-up efficiency will rebate under this condition, the proposed method could speed up the calculation procedure with significant ratios when the transducer number and the number of the discretized points in the x direction increase. For instance, the speed-up ratio reaches about 18.18 when the number of the array transducers is 256 and the number of discretized points in the x direction is 513.
direction of x. Although, the speed-up ratio does not increase constantly along the two factors, it can be distinguished that, with more linear distributed transducers or a higher number of discretized points in the x direction, the speed-up ratio increases. The increment of the speed-up ratio reduces as the transducer number increases. This is due to the fact that, in the proposed method, auxiliary points are created with the same number as the transducer number in the x direction, which can be seen from Figure 5. This can also increase the cost calculation time if the number of auxiliary points are comparable to, or greater than, the number of the discretized points in the x direction. Although the speed-up efficiency will rebate under this condition, the proposed method could speed up the calculation procedure with significant ratios when the transducer number and the number of the discretized points in the x direction increase. For instance, the speed-up ratio reaches about 18.18 when the number of the array transducers is 256 and the number of discretized points in the x direction is 513. Figure 9. The speed-up ratios for the time cost using the proposed method in comparison with that using the classical method.

Comparison between the Proposed Method and a Popular Simulation Toolbox
An example, examined in MATLAB © R2016a (MathWorks, Inc., Natick, MA, USA) by using mex C++, is presented to be compared to the popular ultrasound simulation toolbox, FOCUS [3,4], which Figure 9. The speed-up ratios for the time cost using the proposed method in comparison with that using the classical method.

Comparison between the Proposed Method and a Popular Simulation Toolbox
An example, examined in MATLAB © R2016a (MathWorks, Inc., Natick, MA, USA) by using mex C++, is presented to be compared to the popular ultrasound simulation toolbox, FOCUS [3,4], which is based on fast near-field. A set of flat linear arrays with 32, 64, and 128 rectangular elements distributed along the x-axis is implemented. The width of the transducer element is 1 mm and the height of the transducer element is 5 mm. The spacing gap along the x-axis is 0.4 mm. Thus, the spacing between each neighboring transducer element is 1.4 mm. The ultrasound waves are propagated in a lossless water medium with a speed of sound of 1500 m/s. The center frequency is 1 MHz, and the sample frequency is 10 MHz. The width of the radiated field along the x-axis is twice that of the array width. The radiated field is discretized as a regular grid with the same spacing between the array element, dx = 1.4 mm. The focal point is set at (0, 0, 6λ). The radiated field is calculated by using the proposed method, based on classical spatial impulse response for each flat rectangular array element, and compared to the results calculated by using FOCUS. Figure 10 shows an example in which the results calculated by using the proposed method is compared with a popular simulation toolbox, FOCUS, which is based on far near field method. The simulation condition is described previously. The accuracy of the spatial impulse response depends on the sampling time. In this example, only 10 times that of the center frequency is used for sampling. It can be seen that FOCUS has a good result in accuracy when using a lower sampling frequency, which is the main contribution of FOCUS. In order to obtain a more accurate result, the classical method of spatial impulse response requires a higher sampling frequency. However, this does not affect the acceleration ability since the proposed method is frequency independent as long as the method used to calculate each single transducer only depends on the relative spatial position between the field points and the center of array transducer.
The time cost to calculate the result for different array element numbers, respectively using the proposed method and FOCUS, is presented in Table 5. Since the time cost depends both on the method, intrinsically, and also on the programming methodology, the comparison is just qualitative so as to show that the proposed method can effectively reduce the calculation time. It can be seen that the quantity of repeating calculations seriously increases if the hidden redundant calculation steps are not taken into account, especially when array aperture expands. However, the time cost by using the proposed method increases marginally. compared with a popular simulation toolbox, FOCUS, which is based on far near field method. The simulation condition is described previously. The accuracy of the spatial impulse response depends on the sampling time. In this example, only 10 times that of the center frequency is used for sampling. It can be seen that FOCUS has a good result in accuracy when using a lower sampling frequency, which is the main contribution of FOCUS. In order to obtain a more accurate result, the classical method of spatial impulse response requires a higher sampling frequency. However, this does not affect the acceleration ability since the proposed method is frequency independent as long as the method used to calculate each single transducer only depends on the relative spatial position between the field points and the center of array transducer. The time cost to calculate the result for different array element numbers, respectively using the proposed method and FOCUS, is presented in Table 5. Since the time cost depends both on the method, intrinsically, and also on the programming methodology, the comparison is just qualitative so as to show that the proposed method can effectively reduce the calculation time. It can be seen that the quantity of repeating calculations seriously increases if the hidden redundant calculation steps are not taken into account, especially when array aperture expands. However, the time cost by using the proposed method increases marginally.   Since the proposed speed-up method is based on solving the repeating calculation expressions which only depends on the relative position between the space field point and the array element, any method conforming to this principle can be used as the method to calculate the wave field radiated from each transducer.
The repeating of array transducers becomes crucial when the array aperture expands. The previous calculation methods can be divided into two approaches, as illustrated in Figure 11. First, like the spatial impulse response, the impulse response, which is dependent on the relative spatial position, was calculated and then the transient radiated field is calculated from the convolution between the superimposed impulse response and the first derivation of the impulse signal as expressed in Equation (17). Second, like FOCUS, which calculates the harmonic wave field for a series of frequencies, and then calculates the transient field by using the inverse Fourier transform at each field point. Both methods involving the spatial impulse response or harmonic waves are dependent on the spatial relative position between the field points and each array transducer. As described previously, the repeating calculation is hidden if the field grid is discretized casually, which can be solved by using a specific regular grid method as presented in this paper. The proposed method can be beneficial for any other analytical or semi-analytical calculation method in which the main steps depend on the relative spatial position between field points and the array elements, no matter whether or not it breaks the transducer elements into subsections or has any approximation for the single element.
Intrinsically, the proposed method does not change the calculation steps for each single transducer, although it improves the calculation steps, so the accuracy depends just on the calculation method of a single transducer element and the sampling frequency.
On the other hand, a parallel computing step is also hidden in the first step and the last step, which is not referred to in this paper, however this can also be a beneficial tool to improve the calculation speed by designing an effective parallel program by using multi-core computer.
Moreover, other than the advantages of the proposed method described previously, it, in itself, does have some disadvantages. First, the proposed method cannot deal with imaging problems [26], since the imaging functions are non-linear related to spatial components. Second, the proposed method cannot deal with irregular grids that do not conform to the distribution rule of the array elements, since regular grids are the key to solving the repeating calculation problem. Although it has some limitations, it is accurate enough if the grid spacing is equal or less than the element spacing, since the majority of array elements are distributed with comparative spacing to the wavelength. method conforming to this principle can be used as the method to calculate the wave field radiated from each transducer.
The repeating of array transducers becomes crucial when the array aperture expands. The previous calculation methods can be divided into two approaches, as illustrated in Figure 11. First, like the spatial impulse response, the impulse response, which is dependent on the relative spatial position, was calculated and then the transient radiated field is calculated from the convolution between the superimposed impulse response and the first derivation of the impulse signal as expressed in Equation (17). Second, like FOCUS, which calculates the harmonic wave field for a series of frequencies, and then calculates the transient field by using the inverse Fourier transform at each field point. Both methods involving the spatial impulse response or harmonic waves are dependent on the spatial relative position between the field points and each array transducer. As described previously, the repeating calculation is hidden if the field grid is discretized casually, which can be solved by using a specific regular grid method as presented in this paper. The proposed method can be beneficial for any other analytical or semi-analytical calculation method in which the main steps depend on the relative spatial position between field points and the array elements, no matter whether or not it breaks the transducer elements into subsections or has any approximation for the single element. Intrinsically, the proposed method does not change the calculation steps for each single transducer, although it improves the calculation steps, so the accuracy depends just on the calculation method of a single transducer element and the sampling frequency.
On the other hand, a parallel computing step is also hidden in the first step and the last step, which is not referred to in this paper, however this can also be a beneficial tool to improve the calculation speed by designing an effective parallel program by using multi-core computer.
Moreover, other than the advantages of the proposed method described previously, it, in itself, does have some disadvantages. First, the proposed method cannot deal with imaging problems [26], since the imaging functions are non-linear related to spatial components. Second, the proposed method cannot deal with irregular grids that do not conform to the distribution rule of the array elements, since regular grids are the key to solving the repeating calculation problem. Although it has some limitations, it is accurate enough if the grid spacing is equal or less than the element spacing, since the majority of array elements are distributed with comparative spacing to the wavelength.

Conclusions
A method is proposed to compute the spatial impulse response for 1-D linear phased array transducers. This method is based on the classical spatial impulse response calculation method for rectangular transducers, without aperture discretization and superposition, for rectangular apertures. The rectangular aperture surface is divided into several sub-regions. The relationship between the possibilities of discontinuities and the geometry of the aperture is proposed. The spatial impulse response functions are derived for each sub-region. Under the assumption that the transducer elements have the same rectangular geometry, an analytical expression is derived for the spatial impulse responses from two different field points sharing the same relative distance from two

Conclusions
A method is proposed to compute the spatial impulse response for 1-D linear phased array transducers. This method is based on the classical spatial impulse response calculation method for rectangular transducers, without aperture discretization and superposition, for rectangular apertures. The rectangular aperture surface is divided into several sub-regions. The relationship between the possibilities of discontinuities and the geometry of the aperture is proposed. The spatial impulse response functions are derived for each sub-region. Under the assumption that the transducer elements have the same rectangular geometry, an analytical expression is derived for the spatial impulse responses from two different field points sharing the same relative distance from two different transducers. Simulation with different numbers of array transducers for different dimensions of field points have been carried out. The time cost for the calculation of the spatial impulse response is reduced by implementing the proposed method. The speed-up ratio varies with the number of array transducers, and also the number of discretized points on direction of x. Since this method is based on analytical derivation without approximation or nonlinear transformation, the results are accurate. The effort to solve the repeating calculation problem regarding the spatial-dependent expressions is important because unnecessary calculations can seriously reduce the calculation speed.
Moreover, it may be beneficial if combined with other speed-up methods, for instance, the approximation method proposed in [13] or the nonlinear transformation in the z-axis direction proposed in [19]. This method may have potential to improve the design speed of ultrasound linear phased array systems.