Target acoustic field and transducer state optimization using Diff-PAT

Acoustic holograms are ubiquitous in ultrasonic applications in medicine, biology, physics, and chemistry. Consequently, the accuracy and precision of acoustic holograms are becoming a significant point of interest. In mid-air ultrasonics, the phased array transducer (PAT) is the most popular method to render acoustic holograms; PATs have been designed by multiple authors.^1–3 They have been applied in many applications, such as in airborne ultrasonic tactile displays (AUTDs),⁴ in acoustic levitation,⁵ and for display purposes.^6,7 Acoustic hologram optimizers for PATs have also been investigated. Long et al. first developed an eigensolver and Tikhonov regularization based method (ES) in 2014,⁸ and Marzo and Drinkwater developed an iterative back propagation (IBP) method in 2018.⁵ Plasencia et al. developed the Gerchberg–Saxton (GS-PAT) and the “corrected version of the eigensolver and Tikhonov regularization based method” (CES) in 2020.⁹ The GS-PAT and IBP are modified versions of the Gerchberg–Saxton (GS) algorithm. Moreover, Sakiyama et al. developed the Levenberg–Marquardt (LM) algorithm based method in 2020,¹⁰ and Suzuki et al. demonstrated a Greedy Algorithm (GA) with brute-force search method in 2021.¹¹ Most recently, Fushimi et al. developed a new acoustic hologram optimization method based on the stochastic gradient descent algorithm and automatic differentiation called Diff-PAT.¹² 

While many optimization methods and PATs have been developed, these optimizers can only optimize for phase or amplitude at the target point. In addition, only “phase” and “phase and amplitude” control methods have been explored for the side of transducers. The currently explored and unexplored combination of transducers and target point control methods in PAT systems is summarized in Fig. 1; much of the combinations are yet to be explored. Thus, in this paper, we developed a unifying methodological framework for sound-field computation that allows all prior approaches in PATs to be modeled (e.g., IBP and GS-PAT) while also exploring other approaches that, even if, in principle, applicable,^13–18 had never been adapted to PATs. Achieving acoustic holograms with “amplitude” only control of transducers may contribute to further simplification of PATs. In addition, more advanced access to pressure amplitude and phase at the target point could offer a new opportunity in ultrasonic applications such as acoustic manipulation,^5,19 acoustic streaming,²⁰ and ultrasonic tactile display.^4,8

Here, Diff-PAT was selected as the basis of the generalized acoustic hologram optimizer because Diff-PAT was shown to be more effective than a conventional algorithm based on eigensolver or GS methods. Furthermore, it is simple to implement different modes of wave propagation and loss functions in Diff-PAT. By exploring the unchartered opportunities in the field, this paper also aims to demonstrate the flexibility of the newly developed optimizer. The basic components for the gradient-descent based hologram optimizer are similar to the original Diff-PAT. The acoustic field at the target point (p_a) is calculated using linear superposition of pressure contribution from each transducer,

where T is the total number of transducers, d(x_c, x_t) is the Euclidean distance between x_c and x_t, and P_ref is the reference pressure amplitude. Here, a single sided PAT with 16 × 16 (T = 256) transducers is assumed; $D(θ)=2J1(kr⁡sin⁡θ)kr⁡sin⁡θ$ is the far-field directivity function for the piston source based on the angle θ, which is evaluated as the angle between the transducer normal and target point x_c; r = 5 mm is the transducer radius, and J₁ is the Bessel function of the first kind of order 1; $k=2πfc0$ is the wavenumber given the acoustic frequency, f = 40 kHz, and the speed of sound c₀; P_ref and c₀ are assumed to be 1 Pa at 1 m and 346 ms⁻¹, respectively; p_a is a complex number; the amplitude $Aac=pa$⁠, and the phase $ϕac=Arg(pa)$ at the target point are calculated accordingly. The optimization function is defined as

The contents and size of the optimization variable τ depend on the transducer drive mode (phase, amplitude, and P&A). For phase only control, τ = ϕ_t with size T and A_t = 1 were defined for all transducers. For amplitude only control, τ = A_t with size T and ϕ_t = 0 for all transducers. For P&A control, τ = [Ψ] with $Ψ∈C$ and size T were defined. Since passing on the complex number to the optimizer was considered unfavorable, the input values were separated into real and imaginary numbers; $τpa=[R(Ψ)I(Ψ)]$ with size 2T in the actual implementation. The amplitude and complex angle of Ψ were used to determine A_t and ϕ_t. For both amplitude only and P&A control, the amplitude of the transducer is limited to values between 0 and 1. In order to enforce this, the amplitude of the transducer was normalized as $At=τmax(τ)$⁠.

Different mode of control at the target point is achieved by customizing the loss function for each purpose. In general, the loss function is defined as

This method is similar to methods proposed by Li et al.,¹⁷ who first demonstrated P&A control at the target plane using a single phase plate. Eq. (3) determines the loss function by evaluating the difference of real and imaginary parts separately. Similar to the optimization variable τ, the target variable Ω depends on the mode at the target point. Here, we assume G number of target points. For phase only control at the target point, Ω = ϕ_g with size G is set, while A_a = A_g = 1. For amplitude only control at the target point, Ω = A_g with size G is set, while ϕ_a = ϕ_g = 0. For P&A control, Ω = [A_g, ϕ_g] with size 2C. The methods proposed above are examples of the many different implementation approaches. In all cases considered in this article, the optimizers are tasked to optimize the variables until the iteration number reached n = 1000. These loss functions as defined above were programmed in Python 3.7.10, and the derivative of the loss function with respect to the optimization variables τ was automatically differentiated using JAX²¹ (ver. 0.2.17). The Adam optimizer²² was used to update the current guess of τ_n: $τn=τn−1−α⋅m̂n(v̂n+ε),$ where $m̂n=mt1−β1n$ and $v̂n=vn1−β2n$ are the bias-corrected first and second raw moment estimates, respectively. $mn=β1mn−1+(1−β1)δLδτ$⁠, and $vn=β2vn−1+(1−β2)δLδτ2$⁠. The hyperparameters α, β₁, β₂, and ɛ are 0.1, 0.9, 0.999, and 1 × 10⁻⁸, respectively. These values are arbitrary; further improvements may be achieved by tuning these parameters. The initial guess of τ was generated randomly. The sets of codes necessary to fully reproduce all results and graphs in this paper are made available as shown in the Data Availability statement.²³ 

The generalized acoustic hologram optimizer can combine any propagation modes and loss functions and allows advanced acoustic field control. Figure 2 shows an example field when P&A control is specified for both transducer and target point sides. This can generate a phase surfer¹⁹ using PATs and effectively transport particles at the surface of the water using the phase gradient.⁵ This offers a new opportunity in particle manipulation or acoustic streaming.

To further assess the effectiveness of the optimizer in each combination, we assigned a set of target positions, amplitudes, and phases. Here, two cases were studied with a varying number of focal points, N = 2 and 4. For each case, 1000 sets of target points were randomly generated such that they lie within region of interest (ROI) (x = [−0.05, 0.05], y = [−0.05, 0.05], and z = [0.05, 0.15] m). The pressure amplitudes were assigned such that the minimum pressure amplitude for each point is 10 Pa and such that the sum of each target pressure adds up to the maximum acoustic pressure amplitude of the PAT (identified by averaging the acoustic pressure amplitude of a single focus at the vertex of ROI; in this case, it is 950 Pa). The target phases were assigned such that they lie between [0, 2π]. The performance of the optimizer was measured using two indexes: target pressure accuracy $(Ra=AgAac)$ and phase error $(Rp=cosϕg−cos⁡ϕac)$⁠. The result of performance assessments is shown in Figs. 3(a) and 3(b). As it can be confusing to have phase control for a target point and transducer; a combination of letters and numbers will be used to specify each combination. For example, [A], [B], and [C] specify “phase” only, “amplitude” only, and “P&A” control of transducers, respectively; [i], [ii], and [iii] specify “phase” only, “amplitude” only, and “P&A” control at target points, respectively. Furthermore, A:i is a combination of phase only transducer control with phase only control at the target point, and C:iii denotes a combination of P&A for both the transducer and target point.

It is evident that [i] can be achieved well in all transducer driving modes and control number points, as shown in Figs. 3(a) and 3(b). The average phase error (R_p) for A:i, B:i, and C:i is 9.8 × 10⁻⁴, 1.61 × 10⁻⁵, and 9.7 × 10⁻⁶, respectively, when N = 2 and 1.26 × 10⁻⁵, 6.73 × 10⁻⁶, and 8.08 × 10⁻⁵, respectively, when N = 4. [B] achieves good phase error in comparison with [A] and [C].

However, when [ii] is assigned, the performance of [B] drops significantly in comparison to that of [A] and [C]. The transducer amplitude only method is considerably more difficult than phase-only control since the acoustic field can only be defined via “the relative position of the focal point and the transducer position” and “the amplitude of each transducer.” The average target pressure accuracies (R_a) for A:ii, B:ii, and C:ii are 1.00, 0.603, and 1.00, respectively, when N = 2 and 1.00, 0.859, and 1.00, respectively, when N = 4. The lack of performance by [B] is attributed to the fact that on average, 45.4% and 44.0% of the transducers in PATs were off (off state is defined when transducer amplitude <1%) when N = 2 and N = 4, respectively (see the supplementary material for an example of transducer amplitude distribution). Thus, the achievable target pressure amplitudes by PATs drop significantly and cannot achieve high pressure amplitudes when assigned. The improved performance of B:ii when N = 4 (in comparison to B:ii in N = 2) is attributed to the pressure assignment method requiring each target point to “share” the maximum pressure amplitude (thus reducing the assigned target amplitude average). Since some transducers are already at 100% capacity (see the supplementary material), uniform scaling of the transducers to improve the performance is not possible.

Similar trends are observed in [iii] when both phase and amplitude are specified at the target point. The average target pressure accuracy (R_a) for A:iii, B:iii, and C:iii is 1.00, 0.501, and 1.00 when N = 2 and 1.00, 0.665, and 1.00, respectively, when N = 4. The average phase error (R_p) is 8.48 × 10⁻⁴, 0.0907, and 9.03 × 10⁻⁴, respectively, when N = 2 and 9.66 × 10⁻⁵, 0.0957, and 1.35 × 10⁻⁴, respectively, when N = 4.

The performance difference between A:ii and C:ii or A:iii and C:iii is not statistically significant (i.e., equivalent performance) in matching the target specifications (see the supplementary material for details on statistical analysis). Hence, changing from [A] to [C] does not significantly change the performance, and the benefit of adding amplitude control in transducer side is not apparent solely from Fig. 3.

On the other hand, a close inspection of the convergence graph demonstrates that in some circumstances, the loss value of [A] is lower than that of [C] (see the supplementary material). Since the solution set of [A] is a subset of [C], [C] should be better than [A]. We believe the increased number of optimization variables and complexity of optimizing both amplitude and phase as optimization variables decrease the performance of the [C] optimizer. However, as demonstrated in Fig. 3, the performance between [A] and [C] is not significant. In addition, the current performance of [C] is sufficiently better than that of traditional optimizers in the field.^9,12 Further exploration to improve [C] optimizer may not be necessary.

Figures 4(a) and 4(b) are examples of acoustic field as optimized by A:ii and C:ii; they show that the output field is significantly different between (a) A:ii and (b) C:ii. The ambient noise surrounding the focal point is qualitatively reduced in Fig. 4(b) despite using the same loss function (i.e., the loss function only evaluates the acoustic amplitude at the target point). To assess this ambient noise more quantitatively, we calculated the sum of the pressure amplitude in the x–y plane (square area with a side length of 0.03 m and focal point at the center) for A:ii and C:ii (excluding the circular region with a radius of 2.5λ around the focal point, see Fig. 4) using 1000 sets of optimized field data from Fig. 3 (N = 2). The results are as shown in Fig. 5, and we quantitatively confirm the reduction in ambient noise in C:ii (this trends are continued to be seen in N = 4, 8, 16, and 32 as shown in the supplementary material). The reason behind this reduction in ambient noise is unclear, and further study is required to understand the difference in results. The reduction in ambient noise around the focal point is suitable for applications such as airborne ultrasonic tactile displays (AUTDs),^4,8 and Diff-PAT with [C] could be better suited for AUTDs than [A].

In conclusion, we demonstrated a generalized acoustic hologram optimizer that can control phase, amplitude, and phase and amplitude at the transducer side and target point. The optimizer was based on Diff-PAT, and various loss functions and optimization variables were developed to achieve the target combination. The performance of the optimizer with each combination was quantitatively assessed; we discussed the potential and limitations. We showed that [A] and [C] perform equally when assessed at the target point. However, around the focal point, [C] suppresses the ambient noise and makes it more suitable for use in AUTDs than [A]. PAT has already been demonstrated to be a versatile and useful tool for ultrasonic applications. A more advanced control of PATs and the insights we obtained in this study will aid further development of applications in the field.

See the supplementary material for the transducer amplitude distribution, statistical analysis, convergence graph, and broader assessment of noise reduction in P&A.

This work was funded by the Strategic Research Platform toward Digital Nature powered by Pixie Dust Technologies, Inc. This work was partially supported by JSPS KAKENHI, Grant No. 21K14103. We would like to thank Editage (www.editage.com) for English language editing.