Acoustic hologram optimisation using automatic differentiation

Acoustic holograms are the keystone of modern acoustics. They encode three-dimensional acoustic fields in two dimensions, and their quality determines the performance of acoustic systems. Optimisation methods that control only the phase of an acoustic wave are considered inferior to methods that control both the amplitude and phase of the wave. In this paper, we present Diff-PAT, an acoustic hologram optimisation platform with automatic differentiation. We show that in the most fundamental case of optimizing the output amplitude to match the target amplitude; our method with only phase modulation achieves better performance than conventional algorithm with both amplitude and phase modulation. The performance of Diff-PAT was evaluated by randomly generating 1000 sets of up to 32 control points for single-sided arrays and single-axis arrays. This optimisation platform for acoustic hologram can be used in a wide range of applications of PATs without introducing any changes to existing systems that control the PATs. In addition, we applied Diff-PAT to a phase plate and achieved an increase of > 8 dB in the peak noise-to-signal ratio of the acoustic hologram.

Acoustic hologram is becoming an imperative part of a wide range of acoustics applications such as in the fields of medicine 1-3 , biology [4][5][6][7][8] , and engineering [9][10][11] . The reconstruction accuracy of the acoustic field from the hologram plays a significant role in determining the performance of the system. Thus, developing an acoustic hologram optimiser that can generate holograms with accurate field reconstruction is of significant interest of the field. Recent advancement of airborne phased array transducers (PATs) have enabled new applications such as airborne ultrasound tactile display (AUTD) 12,13 and acoustic levitation [14][15][16][17] . Acoustic levitation in particular is used for digital fabrication 18,19 , display applications [20][21][22][23][24][25] and sample holding in medicine 26,27 , physics 28,29 and chemistry 30 . In all of these applications, acoustic holograms must generate a pressure control/focal point at a specified position. The control point is modulated at low frequencies to create haptic sensations that are sensible by human hands in an AUTD 12,13 . Marzo et al. demonstrated that the control points can be converted into levitation points by adding a twin-trap hologram in acoustic levitation 14,15 .
Generating an acoustic hologram for a single control point in space using PATs is trivial; however, it has been a significant challenge to generate a hologram that can create more than one control point in space. Multiple control points are becoming a necessary part of PAT applications, and low-quality acoustic hologram can lead to poor performances during practical use in applications. To address this challenge, Long 32 . Most recently, Sakiyama et al. 33 demonstrated acoustic hologram optimisation with the Levenberg-Marquardt algorithm (LMA); they examined the accuracy of ultrasonic stimulation in real and simulated environment. It should be noted that, while there are a number of hologram optimisers available within the research community, the usage of acoustic holograms to generate multiple focal points is a recent trend in the PAT community, and this field is yet to be developed. The listed acoustic hologram optimisers above are the most cited or considered state of the art in the field. LMA and IBP optimise only the phase of the transducer array, and GS-PAT and Eigensolver can optimise both the amplitude and phase of the acoustic hologram. Placensia et al. 31 demonstrated that hologram optimisation methods with amplitude control (i.e. GS-PAT and Eigensolver) achieve higher-quality holograms than phase-only methods such as IBP, and that the Eigensolver method achieves the best optimised result because it seeks the global solution.
In this paper, we propose Diff-PAT, a phase-only, gradient-descent algorithm with automatic differentiation as a new platform for optimizing acoustic holograms. We demonstrate that in the most fundamental case of matching the output pressure to the target value; Diff-PAT achieve higher accuracy than the conventional www.nature.com/scientificreports/ Eigensolver and GS-PAT with both amplitude and phase control. Automatic differentiation is different from other differentiation operations (such as numerical differentiation and symbolic differentiation). This method calculates a derivative of a function at high accuracy by applying chain rules to each elementary operation of a given function. It is commonly used in machine learning. Our work was inspired by Peng et al. 34 , who demonstrated this automatic differentiation based method for optical holograms. In comparison to acoustic holograms, optical holograms have been studied for a longer period, and a number of hologram optimisation methods have been developed. However, Peng et al. demonstrated that the method using automatic differentiation outperforms conventional methods such as the Gerchberg-Saxton method. The challenge in optimising an optical or acoustic holograms is based on the fact that the loss function converts a complex number to a real number (nonholomorphic objective 35 ). Whilst this function is challenging to differentiate analytically (the derivative of the function is not just a constant 35 ), latest automatic differentiation packages [36][37][38] can differentiate these functions with high precision 39 and ease 34 . Diff-PAT does have multiple modes of implementation, and the gradient of the loss function does not necessarily need to be solved using automatic differentiation. Analytical and numerical differentiation can be just as effective for identifying the necessary information; however, it is cumbersome to analytically calculate the loss function (Wirtinger calculus is an advanced concept in methematics 34,40 ), and the implementation with numerical differentiation is not as efficient as automatic differentiation (details in section on computational time). Considering the development cycle of the ultrasonic application, the implementation and execution time should be kept minimal to allow the user to test out different array configurations and loss functions to fit their purpose. Thus, we argue that automatic differentiation is the best mode of implementation for Diff-PAT considering its potential usage in early to late stages of development. Achieving a high-quality acoustic hologram without amplitude control in PATs represents a significant contribution to the literature, as it allows PAT controllers to retain a simple design. Furthermore, we demonstrated the versatility of Diff-PAT by employing it for the optimisation of a phase plate [41][42][43][44] . With regard to phase plates, the iterative angular spectrum approach (IASA) developed by Melde et al. 41 is the standard method for optimizing acoustic holograms for arbitrary acoustic fields. We achieved enhanced performance by simply applying Diff-PAT to the optimisation process, without altering the fundamental design proposed by Melde et al 41 .
The system overview of Diff-PAT is as shown in Fig. 1. We used JAX (ver. 0.2.5) on Python 3.6.9 to perform automatic differentiation and optimisation 38 . To identify a suitable acoustic hologram ( φ n ) for a given control point x c , target amplitude A c , and transducer position x t , we adopted the Adam optimiser 45 , which enables efficient stochastic optimisation with only first-order gradients. Before the optimisation of hologram φ n , we defined the following optimisation problem and objective function : C is the total number of the control points, and p t (φ n , x c , x t ) is the total acoustic pressure at x c (see "Methods" section for details on p t ). In order to evaluate the effectiveness of Diff-PAT in PATs applications, acoustic hologram was generated for three setups of arrays (see Fig. 2a for  [13][14][15]21,31 . Some discrepancy between numerical and experimental results has been reported; however its effect is limited to specific application requirements such as accurate positioning 24,47,52 and is otherwise is rarely considered. The initial phase estimate was randomly generated for all transducers. The gradient is calculated by reverse-mode automatic differentiation according to the computational graph of the objective function L , and the phase is optimised by updating itself with the Adam optimiser on JAX. The Adam optimiser updates the current phase estimate according to the following 45 : are bias-corrected first and second raw moment estimates.
, where t is the iteration number, and m t and v t are both initially zero vectors. The required hyperparameters for the Adam optimiser, α , β 1 , β 2 , and ε were set as; 0.1, 0.9, 0.999, and 1 × 10 −8 , respectively.

Results and discussions
Convergence rate of Diff-PAT. In the following section, the performance of Diff-PAT is evaluated. The best optimiser for PATs can be defined as the most versatile and accurate optimiser. In order to assess this, the performance of algorithms are evaluated using three array configurations and five levels of control point numbers. Control points are randomly generated using the process described in the "Methods" section, and this framework covers basic expectations of PATs in a wide range of applications and configurations available to date.
We show that the optimisation function converges quickly to the target value (evaluated by the ratio between the target and current acoustic pressure amplitude), as shown in Fig. 2b (see Data Availability for results at each iteration). Increasing the number of control points (N) has a negative effect on the convergence rate of the solution, and the length of the error bar (which shows the standard deviation of R p ) increases. However, the increase in the number of transducers has a negligible effect on the convergence rate  www.nature.com/scientificreports/ (Fig. 2b). We confirm that our algorithm converges sufficiently between 100 and 150 iterations in Fig. 2b, and the number of iterations was set to 150 in all evaluations of Diff-PAT in this study. The convergence rate of the algorithm can be further improved by hyperparameter tuning; however, this is beyond the scope of the current study. Here, the sum of squared error was chosen as the loss function; however, it could be evaluated using a number of different loss functions (such as Huber loss and sum of absolute error). The choice of loss function has an impact on the final results, and the sum of squared error gave the best results in the preliminary analysis. Further improvement in the loss function may be possible; however, such investigation is beyond the scope of the current study.
Accuracy of Diff-PAT and comparison to conventional solvers. The performance of Diff-PAT was compared with that of conventional algorithms in PATs which were made available by the authors of GS-PAT in C ++ 31 . Here, we compare Diff-PAT with the Eigensolver, corrected Eigensolver, and GS-PAT. The results are shown in the box-and-whisker plot in Fig. 3. Each algorithm was tasked to optimise the acoustic hologram for the same dataset, which included 1000 sets of randomised control points and amplitude (GS-PAT's iteration number was set to 100, as in their study 31 ). Raw data for Fig. 3 can be accessed as stated in the Data Availability section. For readers wishing to inspect the optimised acoustic pressure field in detail, the exemplary field from each optimiser is available in the supplementary material. The visual inspection of the acoustic pressure field shows that the target is achieved as specified in Diff-PAT. Furthermore, we performed a simple case study to demonstrate that the identified solutions by Diff-PAT are optimised solutions and that Diff-PAT strictly follows limits imposed by physical laws; details of which can be found in the supplementary material. The optimised  www.nature.com/scientificreports/ phases from Diff-PAT were exported to a performance analyser developed by the author of GS-PAT, and the performance was analysed on it to independently verify our solution.
Interestingly, when the control number is small (N = 2), Eigensolver performs poorly (median values for M = 196, 512 and 1024 are 0.737, 0.728 and 0.721, respectively, in Fig. 3). This large deviation is considered to be caused by the tendency of the regularisation policy to "homogenise transducer's output rather than reconstruction accuracy" 31 . This performance can be improved by using the corrected Eigensolver method, and the median value improves to unity for M = 196, 512 and 1024 when N = 2.
In addition, the performance of Eigensolver and the corrected Eigensolver improves when the control point number increases from N = 2 to N = 32 in all types of arrays. Specifically, the interquartile (IQ) range value for the corrected Eigensolver method when N = 2 and M = 196 is 0.058; however increasing N to 32 improves the IQ range value to 0.015. We attribute this to the variability of the target amplitude in a given control point geometry. As the number of control points increases, the average difference between the minimum and maximum amplitude decreases. This is partially because of how the target amplitude is assigned when control points are randomly generated. The sum of the amplitudes assigned to the control points is always equal to a constant (see the "Methods" section), and the difference between the minimum and maximum amplitude for a given control point geometry decreases from 523, 283, and 70 Pa for N = 2, 8, and 32, respectively (M = 196 with corrected Eigensolver). This statement is supported by further analysis of the dataset of the corrected Eigensolver when N = 2 and M = 196. When the control point geometries lie between the IQ range, the average amplitude difference between the control points is 443 Pa. However, beyond the IQ range, the average amplitude difference between the control points increases to 603 Pa.
As shown in Fig. 3, the performance of GS-PAT decreases as the number of control points N increases. When N = 2 and M = 196, the IQ range is 0.012; however, when N is increased to 32, the IQ range increases to 0.238. Conventional solvers (i.e. Eigensolver, corrected Eigensolver, and GS-PAT) show a trade-off between the numbers of transducers and control points. In contrast, Diff-PAT outperforms all solvers in all cases evaluated in Fig. 3. Diff-PAT is robust against amplitude variability within the control point geometries and can consistently achieve the target amplitude despite the reduction in transducer numbers. In addition, while the Eigensolver type method and GS-PAT use both amplitude and phase control, Diff-PAT outperforms both methods with only phase control. It should be noted that Diff-PAT does have an increased background noise in comparison to the other algorithms (see acoustic pressure distributions from each solver in supplementary material). However, the increased background noise as obtained in this case study does not invalidate the effectiveness of Diff-PAT, as this is the direct results of the clearly set objectives in the loss function (i.e. only acoustic pressure at the target point is considered). Whether the increased background noise become an issue is largely dependent on each application. For example, in ultrasonic tactile displays increased background noise could lead to blunt haptic sensation, whereas in acoustic levitation, high background noise away from the control point may not arise as significant problem. It is difficult to ascertain the effect of the increased background noise, sole from the numerical simulation. Thus, the purpose of the paper is to set to introduce the effectiveness of Diff-PAT as the platform for optimizing acoustic systems, and we hope researchers from each field will develop suitable loss function which meet their purpose in their system.

Comparing computational time of each solver. The computational time of the solver is insignificant
for non-active PATs; however, some applications of PATs require consideration of the computational efficiency. Figure 4 compares the computational time of each solver, and all of the solvers were executed on the same highend desktop computer (4.2 GHz Core i7-7700 K, 64 GB RAM). Eigensolver and corrected Eigensolver ( Fig. 4a  and b)  Diff-PAT does not have a computational efficiency that is as high as that of GS-PAT; however, its computational time is comparable to that of Eigensolver type optimiser with M = 512. When M = 512 and N ≤ 4 , Diff-PAT is faster than the Eigensolver type method and when M = 1024, Diff-PAT is faster than the Eigensolver method in all cases of N. In addition, Diff-PAT scales well with an increasing number of transducers as shown in Fig. 4d. Given that the demand for high transducer number PATs is increasing 18,[49][50][51]53 : Diff-PAT is already more suited to optimise large transducer number PATs than an Eigensolver type method, both in terms of accuracy and computational efficiency. Furthermore, there are abundant applications of PATs where the accurate reconstruction of the acoustic field is prioritised over computational speed 18,[20][21][22]54 . Diff-PAT can also be implemented using numerical differentiation as stated earlier. The accuracy of Diff-PAT with numerical differentiation is comparable to that of automatic differentiation; however, its computational efficiency is lower than that of automatic differentiation, www.nature.com/scientificreports/ as shown in the supplementary material. This is due to the fact that numerical differentiation is an O(n) process, whereas automatic differentiation scales better with the increasing number of variables in the loss function 39 .
Application of Diff-PAT for Phase Plate. In this section, we further explore the versatility of Diff-PAT by employing Diff-PAT for the optimisation of a phase plate [41][42][43][44] which is used underwater. It has a significantly higher number of elements than PATs 43 and is well-suited for testing the capabilities of Diff-PAT. The IASA developed by Melde et al. is one of the most popular optimisation methods for phase plates. Similar to the GS-PAT and IBP methods, IASA also applies the Gerchberg-Saxton algorithm 32 . To create a phase plate version of Diff-PAT, the acoustic pressure field was assumed to be a plane wave 41 , and the angular spectrum approach 55 was used to calculate the propagating acoustic field in the loss function (see "Methods" section for details). Figure 5 compares the holograms optimised using IASA and Diff-PAT (see Supplementary Information for the acoustic holograms of IASA and Diff-PAT). The reconstruction accuracy of Diff-PAT is clearly higher than that of IASA. A simple visual inspection verifies that the acoustic hologram optimised using IASA, shown in Fig. 5a-c, has many artefacts and does not resolve the test image well (Fig. 5c). In contrast, the acoustic hologram optimised using Diff-PAT has a significantly improved image with minimal artefacts, as shown in Fig. 5a-b, and the USAF 1951 resolution test chart is clearly resolved in Fig. 5c. Quantitatively, the peak signal-to-noise ratio (PSNR) values for Diff-PAT are 16.4, 20.7 and 16.4 dB (Fig. 5a-c, respectively), and they are at least 8 dB higher the respective PSNRs of IASA. In practical applications of phase plates, these optimised acoustic holograms are encoded and fabricated into a plate 41 using a 3D printer. As shown in the supplementary material, Diff-PAT's acoustic hologram has a spatially higher frequency component in comparison with IASA. This could potentially cause issues during manufacture of the phase plates. The effect of such manufacturing tolerances was evaluated, as presented in the supplementary material, and it was found that the PSNR drops by 3 dB when the standard deviation of normally distributed manufacturing tolerance drops to σ≈20%. If manufacturing tolerances become an issue, some weighted regularization term or a low pass filter can be added in the loss function to reduce the high frequency component of the acoustic hologram. The detailed view of acoustic pressure distribution demonstrates that Diff-PAT achieves this high accuracy in the region of interest (ROI) by diffracting the acoustic wave outside the ROI as shown in the supplementary material. The improved accuracy of Diff-PAT with significantly www.nature.com/scientificreports/ reduced artefacts in the ROI can enhance performance of a wide range of acoustic systems used in medicine 1-3 , biology 4-8 , and engineering 9-11 applications.
In conclusion, we presented Diff-PAT, which optimises only the phase of an acoustic hologram using automatic differentiation. We demonstrated the effectiveness of Diff-PAT as a new platform for optimising acoustic holograms in PATs and phase plates. In addition, contrary to common belief, the phase-only optimiser can achieve better accuracy than an optimiser with both amplitude and phase control in the simplest test case for acoustic holograms. For PATs with high transducer numbers (M = 1024), Diff-PAT is more suitable than  www.nature.com/scientificreports/ Eigensolver type methods, both in terms of accuracy and computational efficiency. In future, this improved acoustic hologram optimiser will improve the performance of many applications by providing higher quality acoustic holograms. In addition, this study, along with that of Peng et al. 31 , demonstrated the effectiveness of implementing automatic differentiation in the context of hologram optimization. We hope that these results will further facilitate the use of Diff-PAT platform in a wide range of fields and applications in acoustics.

Methods
Total acoustic pressure calculation. The total pressure was evaluated by p t = M m=0 P ref d(x c ,x t ) D(θ)e j(kd(x c ,x t )+φn,m) where M is the total number of transducers, P ref is the reference pressure amplitude; d(x c , x t ) is the Euclidean distance between x c and x t , and D(θ) = 2J 1 (krsin(θ)) krsin(θ) is the far-field directivity function for the piston source based on the angle, θ . Moreover, θ is evaluated as the angle between the transducer normal and x c ; r = 5 mm is the transducer radius, and J 1 is the Bessel function of the first kind of order 1. k = 2πf c 0 is the wavenumber given the acoustic frequency, f = 40kHz and the speed of sound is c 0 . P ref and c 0 are assumed to be 1.98 Pa at 1 m (12 V Pk-Pk 25 ) and 346ms −1 , respectively. In the implementation of Diff-PAT, the acoustic pressure calculation is split into pre-calculated (A) and re-calculated (B) components: Pa, respectively. Amplitude A c for each control point is randomly assigned such that the minimum pressure amplitude is 10 Pa, and the amplitude at each control point sums up to the average acoustic pressure for the single focus of each transducer array. This target amplitude assignment method is technically inaccurate as it assumes the acoustic amplitude to be conserved (i.e. acoustic power/energy should be conserved). However, the assignment of target via acoustic energy is impractical as it requires assignment of both target pressure and velocity field. Thus, this target assignment method is employed in this manuscript, and it is considered appropriate since the optimizer can follow the targets specified.
Application of Diff-PAT to the phase plate. For the optimisation of the phase plate, an underwater transducer with a resonance frequency of 2 MHz, diameter of 35 mm, and input amplitude of 1 Pa was assumed at the surface. The input pressure field was assumed to be a plane wave, and the initial estimate of the phase was set to zero for all elements. The speed of sound in water was assumed to be 1480 ms −1 , and the pixel size was assumed to be 150 µm . For simplicity, the acoustic transmission loss through the plate was considered negligible, and the angular spectrum was solved using methods described by Zeng & McGough 55 (following the implementation of k-Wave 56 ). Both IASA and Diff-PAT were programmed in Python 3.6.9, and the phase plate version of Diff-PAT used Tensorflow 36 (ver. 2.3.0) to differentiate the loss function automatically. We also used the Adam optimiser in TensorFlow with the same hyperparameter setting as in the PAT version. IASA was calculated by following the steps described by Melde et al. 41  where A c x, y represents the target image with a screen resolution (N) of 256 × 256 pixels. The optimised acoustic hologram was exported from Python and the resultant acoustic pressure fields shown in Fig. 5 were calculated using k-Wave (ver. 1.3 angularSpectrumCW) 56 on MATLAB. Finally, the PSNR was calculated using the MAT-LAB Image Processing Toolbox.

Data availability
The data that support the findings of this study are available within this article. The data for transducer arrangements, control point geometries and amplitude, output phase by Diff-PAT, and the results from all solvers are available in Zenodo at https:// doi. org/ 10. 5281/ zenodo. 49063 51.

Code availability
The program code for optimisation and analysis is also made available at the link provided in Data Availability.