Method, apparatus, and system for controlling sound generation

Abstract

The present technology provides methods, apparatuses, and systems for high-speed acoustic levitation, for example, for high-speed acoustic holography or other applications. A novel technology is presented that enables high-speed multi-point levitation even in the presence of any sound-scattering surface and demonstrates a process that works well in the presence of any physical object. Embodiments provide a simplified approach for determining the location of traps within a working volume, which may sometimes also be referred to as an acoustic volume or an acoustic chamber. Further, embodiments provide an approach for determining the location of traps by determining the contribution from the scattering surface within the working volume and the contribution from the target object within the working volume.

Description

Technical Field

[0001] The present technology generally relates to methods and systems for applications of acoustic levitation, for example, acoustic holography.

Background Art

[0002] Acoustic levitation is a technology that utilizes the mechanical energy of sound to levitate and manipulate substances. This is described in many papers, for example, "Acoustic levitation in mid-air: Recent advances, challenges and future perspectives" by Andrade et al. published in Applied Physics Letters 116 (2020). This technology has advanced significantly in the past decade with the introduction of two basic technologies, namely, the phased array of transducer (PAT) and acoustic holography. PAT enables dynamic control of a high-density array of sound sources (e.g., a 16×16 ultrasonic transducer), while holography, a wavefront processing technology originally developed in optics, enables PAT to accurately control the sound field in 3D space. Thanks to its ability to levitate almost any type of substance, acoustic holography using PAT has many potential applications in laboratory-on-chip, biology, computational fabrication, and aerial displays. Also, acoustic levitation is emerging as a strong candidate for creating a new mixed reality (MR) interface that can seamlessly fuse the digital and physical worlds, as foreseen in "The Ultimate Display" by Ivan Sutherland published in Proc IFIPS Congr. 65 506-601 (1965).

[0003] Recent advances in high-speed acoustic holography have enabled volumetric displays based on levitation that involve tactile and auditory sensations. Many acoustic levitation techniques assume that the volume in which levitation occurs is either empty or contains nearly flat objects. In other words, current methods usually do not calculate the scattering of sound from physical objects within the volume, even though the presence of such physical objects within the working volume is likely to distort the sound field. It is difficult to perform these calculations in real time.

Prior Art Documents

Non-Patent Documents

[0004]

Non-Patent Document 1

Non-Patent Document 2

Non-Patent Document 3

Non-Patent Document 4

Non-Patent Document 5

Non-Patent Document 6

Non-Patent Document 7

Summary of the Invention

Problems to be Solved by the Invention

[0005] Therefore, the applicant has found a need for an improved technology for acoustic levitation.

Means for Solving the Problems

[0006] In a first approach of the present technique, a method implemented by a computer for controlling the location of a target object within an acoustic volume using an array of transducers, the acoustic volume including scattering objects, is provided. The method includes obtaining a static matrix (H) representing the contribution of each transducer of the array of transducers to each of a plurality of locations on the scattering object, the static matrix not changing when controlling the location of the target object. The method further includes defining a plurality of control points within the acoustic volume, calculating in real time a direct transmission matrix (F) representing the direct contribution from each transducer of the array of transducers to each of the plurality of control points, calculating in real time a scattered transmission matrix (G) representing the scattered contribution from each of the plurality of locations on the scattering object to each of the plurality of control points, determining in real time an extended transmission matrix (E) representing the direct and scattered contributions from each transducer of the array of transducers to each of the plurality of control points, the extended transmission matrix (sometimes also referred to as a complete transmission matrix) being determined using the static matrix, the direct transmission matrix, and the scattered transmission matrix. The method further includes determining control commands for each transducer of the array of transducers for generating an acoustic trap at at least one of the plurality of control points using the extended transmission matrix, the acoustic trap trapping the target object for controlling the location of the target object within the acoustic volume.

[0007] In other words, there is a two-step modeling of the extended transfer matrix where the static matrix is obtained before the direct transfer matrix and the scattering transfer matrix are calculated in real time. Once the extended transfer matrix is obtained, the next step may be to solve for the activation (τ) of the transducer that generates a levitation trap at the target position (i.e., the control point). Then, the transducer can be controlled with the determined control commands, and controlling the transducer uses the sound pressure created by the array of transducers to control the location of the target object. In other words, the location of the target object is controlled using acoustic levitation that binds a single particle to a node of a standing wave using the radiation pressure of sound waves. The transducer may typically be an ultrasonic transducer.

[0008] The extended transfer matrix E may be expressed as E = F + GH, where F is the direct transfer matrix, G is the scattering transfer matrix, and H is the static matrix. There may be L control points, N transducers, and M points on the scattering object within the acoustic volume. Thus, the sizes of these matrices are L×N for E and F, L×M for G, and M×N for H. Typically, taking into account the fact that the inequality L≪N≪M is satisfied in acoustic levitation, the determination of H takes more time than the other matrices. While matrices F and G depend on the positions of the control points, the largest and most computationally costly element of the claimed model, matrix H, does not. Thus, after the geometry of the setup (i.e., the transducers and the scattering object) is determined, H remains constant and need not be recalculated every time the binding position is updated (i.e., the part related to the setup). On the other hand, F and G must be calculated every time for an interactive application (i.e., the part related to the application), but these calculations are very suitable for parallel computing. Thus, once matrix H is pre-calculated once, the extended transfer matrix can be calculated very quickly. This two-part modeling means that it is possible to calculate the extended transfer matrix in real time, i.e., when the target object is being controlled.

[0009] The static matrix may be obtained by calculating the static matrix in the setup phase. In other words, the static matrix may be calculated before the target object is placed within the acoustic volume. For example, the method may include defining a plurality of locations on a scattering object within the acoustic volume, obtaining location information for each of the plurality of locations, obtaining position information for each transducer of an array of transducers, calculating, for each of the plurality of locations, a set of acoustic pressure contributions from each transducer of the array of transducers, and storing each set of acoustic pressure contributions in the static matrix. The position information for each transducer may include the position and normal of each transducer.

[0010] The plurality of locations may be a plurality of mesh elements. In other words, the surface of the scattering object may be covered by a plurality of mesh elements. In this example, the location information may include one or more of the position, area, and normal of each mesh. Each mesh element may have a maximum length that is less than λ, λ / 2, λ / 4, or λ / 6, where λ is the wavelength of the sound generated by each transducer. Each mesh element has a maximum length of λ / 2 because this is the mesh size that best balances speed and accuracy.

[0011] The scattering object may change location and / or shape over time. Although the static matrix is fixed, a plurality of static matrices may be calculated, one for each time step. This calculation may be performed in advance. The method may include using these plurality of static matrices to determine an extended transfer matrix for a plurality of time steps. Since the computational load of the static matrix is performed in the setup phase, i.e., offline, the extended transfer matrix may still be determined in real time.

[0012] The step of determining the control command is the trapping stiffness (∇ 2 U jTo maximize

[0013]

Number

[0014] it may include optimizing the phase of each transducer in the transducer array. The Laplacian (∇ 2 U j ) of the Gor'kov potential at point j may be used as an indicator for optimizing the trapping stiffness. The optimization may include maximizing the trapping stiffness using a cost function, and calculating the cost function may include sampling the sound pressure at several control points around each acoustic trap. Alternatively, the optimization step may include a simplified cost function. In other words, calculating the cost function may include sampling the sound pressure at only two control points per acoustic trap. The two control points may be along the main axis of the transducer array.

[0015] Accordingly, the method may include defining each position of the acoustic trap using a plurality of control points, determining the main axis of the transducer array, sampling the sound pressure values at two locations along the main axis around each position of the acoustic trap, calculating an indicator of the trapping stiffness using these sampled sound pressures, and maximizing the calculated indicator of the trapping stiffness using a cost function. Such a method represents a simplified solver, and it will be understood that the simplified solver may be used independently of or together with the model consisting of the two parts of the transfer matrix described above.

[0016] Accordingly, in another aspect, there is provided a method implemented by a computer for controlling the location of a target object within an acoustic volume using an array of transducers that generate sound, where the acoustic volume includes scattering objects. The method includes defining a plurality of control points within the acoustic volume, determining an extended transfer matrix that represents the direct contribution from each transducer of the array of transducers to each of the plurality of control points and the scattered contribution from each transducer to each of the plurality of control points via the scattering objects, and using the extended transfer matrix to determine control commands for each transducer of the array of transducers for generating an acoustic trap at at least one of the plurality of control points. The determining step is performed by defining each position of the acoustic trap using the plurality of control points, determining the main axis of the array of transducers, sampling the sound pressure values at two locations along the main axis around each position of the acoustic trap, calculating an indicator of the confinement stiffness using these sampled sound pressures, and maximizing the calculated indicator of the confinement stiffness using a cost function. The acoustic trap confines the target object to control the location of the target object within the acoustic volume. The extended transfer matrix may be determined using the two-step process described above.

[0017] The indicator may be the simplified Gor'kov metric U j ', and

[0018]

Number

[0019] may be defined as follows, where V represents the volume of the target object, ω represents the angular frequency of the target object, c and ρ represent the speed of sound and density, the subscripts 0 and p refer to the host medium (i.e., air) and the particle material, respectively, and p jrepresents the sound pressure at the control point from the j-th transducer, and z is the main axis. The constants K1 and K2 are determined by the physical properties of the particles and air, and are constant values determined by the Gor'kov equation, not weights. The cost function may be described as

[0020]

Number

[0021] and may be described as w s is a weight coefficient and may be fixed, for example, at 0.0001. Any suitable optimization algorithm such as Broyden-Fletcher-Goldfarb-Shanno (BFGS) or gradient descent method may be used to minimize this cost. There may be multiple iterations of the optimization algorithm, for example, 100 iterations.

[0022] There may be multiple target objects, and the number of traps may be selected to match the number of target objects. The number of target objects (and thus the generated traps) may vary and may range, for example, from 1 to 16 depending on the application. It will be understood that more may be used in some cases.

[0023] The properties of the target object also depend on the application. For example, the target object may be solid particles, such as expanded polystyrene (EPS) particles. Alternatively, the target object may be liquid particles, such as printing fluid, resin, or water. Thus, the method may be used in a printing process to control the location of the droplets to be printed. In other words, there may be a printing method, and the method includes controlling the first location of the target object in the form of a plurality of printing droplets to change the state of each droplet from liquid to solid, and controlling the second location of each of the plurality of solid printing droplets to deposit each printing droplet at a desired location. The step of controlling the first location and the step of controlling the second location are performed using the method described above. In this way, each printing droplet may be deposited, resulting in an additive assembly, sometimes also referred to as 3D printing. As a real-time example of a practical application of particle manipulation, there may be 50 frames per second (fps) for manipulating particles at a step size of 0.2 mm at 1 cm / s.

[0024] Another application is a display in which a volumetric image in the air is created by utilizing the principle of persistence of vision. This may be achieved, for example, by using a projection screen on which an image is displayed, and the projection screen is supported by several particles, for example, four particles - one at each corner - and the movement of the projection screen may be controlled by controlling the movement of each particle as the target object. Alternatively, the volumetric image may be generated by moving a plurality of particles to show a volumetric image (i.e., a 3D shape). The plurality of particles must be moved fast enough for the persistence of vision to be utilized. For applications using persistence of vision to create an image, the update rate may be as high as 10,000 fps.

[0025] In other words, there may be a method for generating a moving volumetric image, the method including providing a plurality of particles, or a screen supported by a plurality of particles, and controlling the location of each of the plurality of particles as a target object using the above-described method, wherein the moving volumetric image is generated by the movement of the plurality of particles or the movement of the screen.

[0026] In a related approach, there may be provided an apparatus including an array of transducers for generating a sound pressure, an acoustic volume defined by the sound pressure generated by the array of transducers within which the location of a target object is controllable, and a processor for performing the above-described method to control the movement of the target object within the acoustic volume.

[0027] In the related method of the present technology, it includes an array of transducers for generating sound pressure, an acoustic volume defined by the sound pressure generated by the array of transducers, in which the location of the target object is controllable, and a processor. The processor is to obtain a static matrix (H) representing the contribution of each transducer of the array of transducers to each of a plurality of locations on the scattering object, where the static matrix does not change when controlling the location of the target object; to define a plurality of control points within the acoustic volume; to calculate in real time a direct transfer matrix (F) representing the direct contribution from each transducer of the array of transducers to each of the plurality of control points; to calculate in real time a scattering transfer matrix (G) representing the scattering contribution from each of the plurality of locations on the scattering object to each of the plurality of control points; to determine in real time an extended transfer matrix (E) representing the direct and scattered contributions from each transducer of the array of transducers to each of the plurality of control points, where the extended transfer matrix is determined from E = F + GH using the static matrix, the direct transfer matrix, and the scattering transfer matrix; and to determine control commands for each transducer of the array of transducers for generating an acoustic trap at at least one of the plurality of control points using the extended transfer matrix, where the acoustic trap is configured to confine the target object to control the location of the target object within the acoustic volume. An apparatus is provided that is configured to perform the above.

[0028] In a method related to this technology, it includes an array of transducers for generating sound pressure, an acoustic volume defined by the sound pressure generated by the array of transducers and within which the location of a target object is controllable, and a processor. The processor is configured to define a plurality of control points within the acoustic volume, determine an extended transfer matrix representing the direct contribution from each transducer of the array of transducers to each of the plurality of control points and the scattered contribution from each transducer through a scattering object to each of the plurality of control points, use the extended transfer matrix to determine control commands for each transducer of the array of transducers for generating an acoustic trap at at least one of the plurality of control points, define each position of the acoustic trap using the plurality of control points, determine the main axis of the array of transducers, sample the sound pressure values at two locations along the main axis around each position of the acoustic trap, use these sampled sound pressures to estimate an index of confinement stiffness, and determine by maximizing the calculated index of confinement stiffness using a cost function, where the acoustic trap is configured to confine a target object to control the location of the target object within the acoustic volume. An apparatus is provided that is configured to perform the determining.

[0029] In a method related to this technology, a non-transitory data carrier carrying processor control code for implementing any of the methods, processes, and techniques described herein is provided.

[0030] As will be understood by those skilled in the art, the present technology may be embodied as a system, method, or computer program product. Accordingly, the present technology may take the form of all hardware embodiments, all software embodiments, or embodiments combining software aspects and hardware aspects.

[0031] Furthermore, the present technology may take the form of a computer program product embodied in a computer-readable medium embodying computer-readable program code. The computer-readable medium may be a computer-readable signal medium or a computer-readable storage medium. The computer-readable medium may be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination thereof.

[0032] The computer program code for carrying out the operations of the present technology may be written in any combination of one or more programming languages, including object-oriented programming languages and conventional procedural programming languages. The code components may be embodied as procedures, methods, etc., and may include sub-components that take the form of instructions or sequences of instructions at any level of abstraction, from direct machine instructions of a native instruction set to constructs of a high-level compiler-type or interpreter-type language.

[0033] Embodiments of the present technology also provide a non-transitory data carrier carrying code that, when executed on a processor, causes the processor to perform any of the methods described herein.

[0034] The technology further provides processor control code for implementing the above-described method, for example, on a general-purpose computer system or a digital signal processor (DSP). The technology also provides a carrier carrying processor control code for implementing any of the above methods when executed, particularly on a non-transitory data carrier. The code may be provided on a carrier such as a programmed memory such as a disk, a microprocessor, a CD-ROM or a DVD-ROM, a non-volatile memory (e.g., flash) or a read-only memory (firmware), or on a data carrier such as an optical or electrical signal carrier. The code (and / or data) for implementing embodiments of the technology described herein may include source code, object code, or executable code in a normal programming language such as C (interpreted or compiled), assembly code, code for setting up or controlling an ASIC (application-specific integrated circuit) or an FPGA (field-programmable gate array), or code in a hardware description language such as Verilog (RTM) or VHDL (Very high speed integrated circuit Hardware Description Language). As will be understood by those skilled in the art, such code and / or data may be distributed among a plurality of coupled components that communicate with each other. The technology may include a controller including a microprocessor, a working memory, and a program memory coupled to one or more of the components of the system.

[0035] All or part of the logical method according to the embodiments of the present technology may be preferably embodied in a logic device including logic elements for executing the steps of the above-described method, and it will be apparent to those skilled in the art that such logic elements may include components such as programmable logic arrays or logic gates of application-specific integrated circuits. Further, such a logical arrangement may be embodied, for example, in an enabling element for temporarily or permanently establishing a logical structure within such an array or circuit using a virtual hardware descriptor language, and those enabling elements may be stored and transmitted using a fixed medium or a transmissible carrier medium.

[0036] In embodiments, the present technology may be implemented using a plurality of processors or control circuits. The present technology may be executed on the operating system of the device or adapted to be integrated into the operating system of the device.

[0037] In embodiments, the present technology may be realized in the form of a data carrier having functional data, and the functional data includes a functional computer data structure that enables the computer system to execute all steps of the above-described method when loaded into and operated by the computer system or network.

[0038] As merely an example, aspects of the present technology will be described hereinafter with reference to the accompanying drawings.

Brief Description of the Drawings

[0039]

Figure 1a

Figure 1b

Figure 1c

Figure 1d

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8a

Figure 8b

Figure 8c

Figure 8d

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13a

Figure 13b

Figure 14a

Figure 14b

Figure 15a

Figure 15b

Figure 15c

Figure 15d

Figure 15e

Figure 15f

Figure 16a

Figure 16b

Figure 16c

Figure 16d

Figure 16e

Figure 17a

Figure 17b

Figure 17c

[0040] Broadly speaking, embodiments of the present technology provide methods, apparatus, and systems for high-speed acoustic levitation, for example, for high-speed acoustic holography or other applications. A novel technique is presented that enables high-speed multi-point levitation even in the presence of any sound scattering surface and demonstrates a process that works well in the presence of any physical object. As will be described in more detail below, embodiments provide a simplified approach for determining the location of traps within a working volume, sometimes also referred to as an acoustic volume or acoustic chamber. Further, embodiments provide an approach for determining the location of traps by determining the contribution from the scattering surface within the working volume and the contribution from the target object within the working volume.

[0041] Figure 1a shows a system 100 for generating and controlling sound to provide acoustic levitation. System 100 may include a control device 110 which may be any suitable computing device, such as a personal computer or computing device, laptop, or server, or a combination thereof. Figure 1a shows a part of the components of the control device, and it will be understood that they may further be other standard components not shown. Control device 100 includes at least one processor 112 coupled to a memory 114. The at least one processor 112 may include one or more of a microprocessor, a microcontroller, and an integrated circuit. Memory 114 may include volatile memory such as random access memory (RAM) for use as temporary memory, and / or non-volatile memory such as flash, read-only memory (ROM), or electrically erasable programmable ROM (EEPROM) for storing data, programs, or instructions, for example.

[0042] Control device 100 typically includes at least one input / output interface 116 for at least the user to input commands and / or receive information. The at least one input / output interface 116 may take any suitable form, such as a keyboard, mouse, touchpad, or other input device for inputting commands from the user, and / or a display or other output device for providing results and / or data generated between the methods described below.

[0043] Also, the processor 112 is coupled to the transducer array 122 to control the generation of sound from the transducer array 122. The transducer array 122 is placed within the acoustic chamber 120 (which may also be referred to as an acoustic volume). Within the acoustic chamber 120, there is also at least one target object 124 whose movement within the acoustic chamber 120 is controlled by the generation of sound. A scattering object 126 is also placed within the acoustic chamber 120, and the scattering object 126 affects the generation of sound within the acoustic chamber 120.

[0044] As will be described in more detail below, the calculations within the control device may be divided into two stages. In the first stage, a two-step scattering model 118 is used, and a part of this model known as matrix H may be pre-generated and stored in the memory 114 as shown. The next stage is applied using a simplified floating solver 119. Combining both techniques enables updates in excess of 10,000 times per second to create volume images above and below the sound scattering object.

[0045] Figures 1b through 1d show different setups of the acoustic chamber. In Figure 1b, the acoustic chamber 220 is defined between an upper plane and a lower plane that are generally parallel to each other. An array 222 of transducers is mounted on the upper surface to generate sound toward the lower surface of the acoustic chamber as indicated by the direct sound waves. In this example, the array is 16×16. The scattering object 226 (which may also be called a physical object) is placed on the lower plane, and the direct sound is scattered from the scattering object 226 as indicated by the scattered sound waves. In this example, the target object 224 is a holographic display that is a projection screen (i.e., a light fabric) attached to four corner particles. The four particles, and thus the screen, are levitated, and the movement of the holographic display (also called digital content) within the acoustic chamber is controlled by the sound generated by the transducer array. Such an arrangement shows a mixed reality display that creates digital content in the presence of a 3D printed physical object. The high computational speed of the proposed approach enables the digital content to be interactive with respect to the user's input (i.e., the levitated screen moves according to keyboard input).

[0046] Figure 1c shows an arrangement in which an acoustic chamber 320 is again defined between an upper plane and a lower plane that are generally parallel to each other. In this example, there is an array 322 of transducers mounted on both the upper and lower surfaces to direct sound toward each other and into the acoustic chamber 320. In this example, the scattering object 326 is a sphere suspended within the acoustic chamber 320. As shown by the heat map within the acoustic chamber 320, the process described below can create a plurality of levitation traps 328 in the presence of a physical object that scatters sound. P max represents the maximum amplitude of the pressure within the sound field, and thus the traps occur at locations of maximum amplitude.

[0047] Figure 1d shows an arrangement where two planes are arranged to define a V-shape, and the acoustic chamber 420 is defined as the volume between and above the planes. In this example, there are arrays 422 of transducers attached to each of the planes such that the two arrays direct sound towards each other and into the acoustic chamber 420. Similar to Figure 1c, the scattering object 426 is a sphere suspended within the acoustic chamber 420.

[0048] Model and Solver Figures 2 and 7 are flow diagrams of steps executed by the system to achieve fast multi-point levitation with minimal loss of accuracy even within a non-empty working volume. As described above, the method utilizes the two-step scattering model shown in Figure 2 and the simplified levitation model shown in Figure 7.

[0049] Acoustic levitation, including acoustic holography using PAT, relies on a linear model represented as the transfer matrix F. The matrix F describes how the complex activation of N transducers

[0050]

Number

[0051] contributes to the complex acoustic pressure at L points of interest within the sound field

[0052]

Number

[0053] using the linear system ζ = Fτ, where L ≪ N. Each element of this matrix (F l,n) is equal to the pressure at the l-th point of interest generated by the n-th transducer when the activation of the n-th transducer is 1 (i.e., the maximum amplitude with a phase lag of 0 rad), and can be approximated as a piston model when considering only the direct contribution. Using this ordinary linear model, existing methods use different solvers to obtain the activation (τ) of the transducer that generates an ideal sound field (ζ) that, for example, creates a focus for providing a tactile sensation or provides the maximum confinement stiffness (i.e., an acoustic trap) to levitate particles at a desired position.

[0054] The proposed scattering model is based on the boundary element method BEM. First, how the ordinary BEM works for a general scattering problem is explained, and then how it is reformulated for the method shown in Fig. 2 is explained.

[0055] Ordinary BEM for scattering problems: In BEM, the sound pressure at a point x can be expressed as a boundary integral equation (i.e., the Helmholtz-Kirchhoff integral equation) obtained by Green's theorem. In a scattering problem, BEM can be calculated by discretizing the surface of the scattering object into M mesh elements. The size of the elements is small enough so that the pressure (p m ) can be considered constant within each element. Then, under specific impedance boundary conditions parameterized by β m , the complex pressure (p(x)) in the domain of propagation (i.e., the region where the wave propagates) is given by the contribution of the direct incidence (p inc (x)) and the scattered contributions from all mesh elements

[0056]

Equation

[0057] as given below.

[0058] Here, s m represents the surface area, k is the wave number, and β m is the relative surface admittance at the boundary, calculated as the ratio of the acoustic impedance Z0 of the propagation medium to the acoustic impedance Z s of the scattering object (i.e., β m = Z0 / Z s ), and when the surface is acoustically rigid, β m = 0). G(y, x) is the so-called free-field Green's function defined by

[0059]

Number

[0060] in the 3D case.

[0061] Here, d(x, y) is the Euclidean distance between two points x and y. In equation (a), ∂ / ∂n represents the normal derivative on the boundary (i.e., the rate of increase in the direction of the normal n m of the mesh). Let ψ(x, y) represent the angle between the normal of the mesh at y and the vector x - y, and ∇ y represent the gradient with respect to the components of y. The normal derivative of the Green's function at y can be expressed as

[0062]

Number

[0063] as follows.

[0064] On the other hand, the sound pressure (p m ) on each mesh can be derived from the Helmholtz-Kirchhoff integral equation under the same impedance boundary condition when the surface is smooth around x m .

[0065]

Number

[0066] Equation (d) leads to a set of M linear equations for determining M unknown pressure values (p m ) in the mesh element. The equations can be represented in the form of a simple equation system Ap = b, and the elements of matrix A and vector b are given as follows.

[0067]

Number

[0068]

Number

[0069] When the set of pressure values in the mesh element (p = [p1…p M ) T ) is obtained by solving this equation system, the sound pressure p(x) can be calculated at any position in the propagation field (i.e., at any point within the acoustic volume) using Equation (a). Matrix A depends only on the geometry of the boundary, while vector b depends on the incident wave (i.e., the contribution of the direct sound from the transducer). It should be noted that solving this equation requires an enormous amount of time and memory for large M. In other words, BEM can simulate the sound-scattering field and has been used to levitate objects several times larger than the wavelength, but BEM is generally considered not suitable for real-time applications, especially for high requirements of afterimage display applications (i.e., 10,000 fbps), and no dynamic operations using BEM are shown.

[0070] Conventional BEM Adaptation for the Scattering Problem: BEM can model a sound-scattering object by modeling it as a mesh of M boundary elements (i.e., in this example, a mesh with 3,000 to 6,000 elements may be used). In the first step S200 of adapting the standard BEM method to this technology, a transfer matrix E that captures both the direct contribution of the transducer to the target point and the contribution of scattering is obtained. The matrix E is defined as E = F + GH by three matrices. The first matrix F represents the contribution from the transducer to the point of interest (i.e., F is a normal transfer matrix that captures only the contribution of the transducer and may be called the direct contribution matrix). The second matrix G represents the contribution from the scattering object to the point of interest (i.e., G may be called the scattering contribution matrix). The third matrix H represents the contribution from the transducer to the scattering object (i.e., H may be called the static contribution matrix).

[0071] The elements of each of these matrices F l,n , G l,m , H m,n are schematically shown in FIG. 3. l represents a control point, m represents a point on the mesh, and n represents a transducer within the transducer array. There are L control points, N transducers, and M points on the mesh within the acoustic volume. Therefore, the sizes of these matrices are L×N for E and F, L×M for G, and M×N for H. It is noted that the matrix E may be calculated by repeating the BEM calculation for each of the N single transducers (i.e., N = 256). However, each BEM calculation involves solving a large and dense system of linear equations, and thus it is not realistic to repeat this process for all transducers in real time. Usually, considering the fact that the inequality L≪N≪M is satisfied in acoustic levitation, the determination of H takes more time than the other matrices.

[0072] For a static setup, the Applicants have found that H is constant and thus can be pre - calculated after the setup is defined. In other words, to calculate the transfer matrix (E) at a high update rate, the new process reformulates the BEM into two parts, namely the part / phase related to the setup and the part / phase related to the application, and both of these two parts are shown in Figure 2. In the setup phase, the contribution of each transducer to the mesh is pre - calculated, and these pre - calculated values are used to update the transfer matrix in real - time as the trap position moves. Thus, the extended version of the transfer matrix maintains the efficiency of the empty - volume method while providing accuracy equivalent to that of the conventional BEM.

[0073] Each element of the matrix (E l,n ) is equal to the pressure (p n ) generated at the l - th point with the n - th transducer having a complex activation τ l,n = 1. In this case, since the acoustic impedance of all the sound - scattering surfaces used (i.e., plastic, water) is very high when compared to air, for these sound - scattering surfaces, it is assumed that β m = 0 in equations (a) and (f). Then, p l,n can be expressed as, by using the BEM,

[0074]

Number

[0075] as follows.

[0076] Here,

[0077]

Number

[0078] represents the direct contribution of the n - th transducer to the l - th point, and pm,n represents the pressure in the m-th mesh generated by the n-th transducer. And, as shown in FIG. 3, the complete transfer matrix E (which may also be referred to as the complete transfer matrix or the extended transfer matrix) can be represented as

[0079]

Number

[0080] as follows,

[0081]

Number

[0082] and is

[0083] Direct incidence contribution

[0084]

Number

[0085] is

[0086]

Number

[0087] capable of being represented as follows, P l,n represents the scalar directivity of the sound source approximated as a piston model, and Φ l,n represents the complex phase propagation approximated as a spherical sound source.

[0088]

Number

[0089] Here, P refrepresents the reference pressure of the transducer at a distance of 1 m, r represents the radius of the transducer, and θ(x l , x n ) is the angle between the normal of the transducer and the point l, and J1 represents the Bessel function of the first kind.

[0090] One method for pre-computing the third matrix H is shown in FIG. 2. The process includes defining a plurality of locations on the scattering object, for example, by obtaining mesh information (e.g., the position, area, and normal of each mesh) regarding the reflecting surface of the scattering object S202. The process also includes obtaining the location information (e.g., position and normal) of each transducer S204. Although these steps are shown in parallel, it will be understood that they can be sequentially performed in either order.

[0091] Taking into account the mesh information, i.e., considering the geometry of the sound scattering object, the next step S206 is to construct an intermediate matrix A of size M×M of the matrix A using the formula f. The next step S208 is to construct the incident vector b (n) for the nth transducer, where each element of the vector defines the incident pressure at that point,

[0092]

Equation

[0093] is defined by, P m,n represents the scalar directivity of the sound source approximated as a piston model, and Φ m,n represents the complex phase propagation approximated as a spherical sound source. Using the results of the previous two steps, in step S210, p (n) can be obtained, which represents the set of pressure values at each mesh element due to the nth transducer. This can be obtained by solving Ap (n) = b (n) .

[0094] And the results for the transducers are stored in the appropriate elements of the third matrix in step S212, i.e.,

[0095]

Number

[0096] which is. Then, in step S214, there is a determination as to whether there are any further transducers within the array for which values have not yet been calculated. If there are further transducers, steps S208 through S212 are repeated for all N transducers. If there are no further transducers, the third matrix is complete and can be used in real time to obtain the complete transfer matrix.

[0097] In this application, in step S210, to solve the linear system, the MATLAB function gmres that uses the Generalized Minimum Residual (GMRES) algorithm is used. Another way to represent steps S206 through S210 is like AH = B, where B = [b (1) …b (N) . Instead of using GMRES, the matrix A (e.g., LU decomposition) could be decomposed to calculate H more quickly. It will be understood that other methods may also be used to calculate H.

[0098] In contrast to calculating the third transfer matrix H, calculating the first and second transfer matrices F and G requires the position of the point of interest in addition to the setup information. For interactive applications, these points of interest are usually not known in advance and, therefore, these matrices F and G need to be created in real time depending on the application logic and / or user input. The next step S216 is to obtain the location of the point of interest. For example, the point of interest may be the point where the trap should be placed. Once the point of interest is known, a complete transfer matrix can be obtained in step S218 by using the stored matrix H and by calculating the other matrices F and G in real time. The calculation of these matrices has a direct representation given by equation (i) and is, therefore, well-suited for parallel computation using a graphics processing unit (GPU).

[0099] As already mentioned, in the present application, it is assumed that β m = 0 for all sound scattering surfaces. When solving the matrix H, the term ikβ m' G(x m' , x m ) is left and equation (i) is adjusted to have that term when calculating the matrix G, so that an extension to other values of β m is also possible. This extension can be adopted without increasing the computational effort significantly.

[0100] Determining the phase of the transducer for acoustic 3D manipulation Once it is known how to model the extended transfer matrix (E = F + GH), the next step is to determine the activation (τ) of the transducer that generates the levitation trap at the target position in the presence of the sound-scattering object, and how to control the transducer to achieve the desired trap, as shown in step S220. Then, the transducer can be controlled with the determined control commands. Acoustic levitation utilizes the radiation pressure of sound waves (typically ultrasonic waves) to confine a single particle to a node of a standing wave. Thus, the trap is air-pressure dependent. Only phase optimization is assumed (i.e., the amplitude of the transducer is always maximum), and thus the goal of this optimization is to find the optimal phase of the transducer that maximizes the trapping stiffness (∇ 2 U j ) at every trap position

[0101]

Number

[0102] .

[0103] Three different levitation solvers, namely, BASELINE, HEURISTIC, and SIMPLIFIED, are considered. The BASELINE solver uses stiffness as a physically accurate and widely accepted metric for the trapping quality but is the slowest. The HEURISTIC solver is the fastest but not accurate enough. The SIMPLIFIED solver uses a simplified Gor'kov potential as a new metric instead of stiffness and represents the most balanced optimal solution that enables accurate and fast acoustic manipulation, as explained below.

[0104] BASELINE levitation solver: One straightforward approach in this optimization problem is, as proposed

[0105]

Number

[0106] The cost function determined as follows

[0107]

Number

[0108] is to directly maximize the confinement stiffness (∇ 2 U j ).

[0109] The confinement stiffness is a common metric for evaluating (and optimizing) the quality of acoustic traps, and is calculated as the Laplacian of the Gor'kov potential (∇ 2 U j ) at point j. Here, the bar (￣) represents the average value over all J traps, and w s is the weighting factor. Thus, this conventional method creates levitation traps by maximizing the confinement stiffness at the desired location using an optimization algorithm such as the gradient descent method. The second term of this cost function is added to equalize the quality (i.e., stiffness) of all J traps by minimizing the standard deviation in the same way as described in "Acoustic levitation with optimized reflective metamaterials" by Polychronopoulos et al. published in Sci. Rep 10 4254 (2020). The BASELINE solver uses, for example, the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm described in "On the limited memory BFGS method for large scale optimization" by Liu et al. published in Math. Program 45 503-528 (1989) to minimize this cost function of Equation 1

[0110]

Number

[0111] Minimize it.

[0112] However, calculating the confinement stiffness (∇ 2 U j ) is computationally expensive because it is necessary to sample the pressure values at many points (e.g., 55 points as described in relation to FIG. 4, i.e., L = 55J) around each trap. The reason why calculating the confinement stiffness (∇ 2 U j ) requires so many points is that the second spatial derivative of U j requires up to the third derivative of the pressure value at the trap position, as follows.

[0113]

Equation

[0114] Here, a represents x, y, or z, and the dot operator (·) is defined as p f ·p g = Re[p f Re[p g + Im[p f Im[p g . To numerically obtain the derivative of Equation (b), this index requires sampling the pressure values at many points. Since this index needs to serve as a baseline, a second-order centred difference approximation was used to calculate these derivatives for accuracy.

[0115] FIG. 4 shows how the pressure values at the points around the trap in the ab plane were sampled, where ab ∈ {xy, yz, zx}. Here, p 10 in the xy plane overlaps with that in the other two planes, and p9 and p 11 in the xy, yz, and zx planes overlap with p6 and p 14 in the yz, zx, and xy planes, respectively.Note that they are identical to each other. This means that in this application, a total of 55 points were used per trap (i.e., 21 for each of the xy, yz, and zx planes, excluding 2 + 2×3 = 8 overlapping points).

[0116] HEURISTIC Floating Solver: To simplify this optimization problem, a heuristic approach proposed for the setup of up-down floating was adapted. This approach uses two points of interest per trap, at λ / 4 above and λ / 4 below the position where the trap needs to be placed, such that there is a π - radian offset in the target phase (i.e., L = 2J). Simply backpropagate those points with the conjugate transpose of the transfer matrix (E * ) and then, without any iteration (i.e., K = 0), by constraining the amplitude of the transducer to its maximum value, the phase

[0117]

Number

[0118] can be calculated. This HEURISTIC approach is the simplest and will work well for single - point operations, but there is a high likelihood of destructive interference between traps for multi - point operations.

[0119] This HEURISTIC approach will still work even if the solver uses slightly different positions for two control points at the trap position (x j , y j , z j ) and its slightly upper position (x j , y j , z j + h). The modified version of this HEURISTIC Floating Solver is also used to obtain initial guesses for the BASELINE Solver and the SIMPLIFIED Solver.

[0120] SIMPLIFIED Levitation Solver: The SIMPLIFIED solver developed by the applicant uses a gradient descent method to minimize a simplified metric (U j ') at every trap position, enabling us to create multiple stable traps at high computational speed. The metric (U j ') is based on the Gor'kov potential (U j ) at point j, and the Gor'kov potential (U j ) is

[0121]

Number

[0122] such as the acoustic radiation force applied to a small particle (i.e., much smaller than the acoustic wavelength) at point j

[0123]

Number

[0124] and can be used to calculate. Here, U j can be determined as follows by the complex sound pressure (p j ) at the trap position, its spatial derivative, and constant values (K1 and K2).

[0125]

Number

[0126] In this application, the process of solving to create J traps is accelerated by using the simplified Gor'kov potential (U j ') as a new metric (i.e., the cost function in the gradient descent method).

[0127]

Number

[0128] In the formula, V represents the volume of the floating particles, ω represents the angular frequency, c and ρ represent the sound speed and density, and the subscripts 0 and p refer to the host medium (i.e., air) and the particle material, respectively.

[0129] The important point here is that, in order to numerically calculate the differentiation along the main axis, i.e., the z-axis in this example, the pressure values of only two points on the periphery of each trap at the trap position (x j , y j , z j ) and its slightly upper position (x j , y j , z j + h) (for example, in this application, h = λ / 32 was used) are sampled (i.e., L = 2J, see Figure 5), whereby U j ' can be calculated. For further simplification, instead of the central difference approximation, a first-order forward difference approximation was used instead. Adding the differentiation along the x-axis and y-axis (i.e., using the original Gor'kov potential shown in Equation 3) requires sampling the pressure values of four points on the periphery of each trap (i.e., L = 4J). The simplified indicator enables an update rate that is approximately twice as fast as using the original Gor'kov potential, but the (slower) solution using the original Gor'kov potential requires minimal changes.

[0130] The advantage of this new indicator is that it can be calculated by sampling the pressure values of only two points per trap (i.e., the total number of points of interest is L = 2J). This simplified indicator is suitable for experimental setups such as those shown in Figures 1b and 1c, where the transducer faces downward (i.e., in the -z direction) and the sound-scattering object is placed below it, or the transducer faces both downward and upward (i.e., in the -z direction and the z direction) and the sound-scattering object is between the transducer arrays. In both of these arrangements, acoustic traps such as standing waves are created along the z-axis.

[0131] The simplification of Equation (4) is that the derivative of pressure along the z-axis is dominant over the derivatives along the other axes, i.e., since the z-axis is the principal axis in this example, the potential U j is approximated well. Also, the Gor'kov potential along the z-axis behaves locally as a sine-wave pattern. Therefore, the second derivative of such a sine-wave pattern (i.e., the binding rigidity) should also be sinusoidal with the opposite sign, supporting the assumption that the negative relationship between U j and its Laplacian (∇ 2 U j ) still holds.

[0132] This is demonstrated in Figure 6, which shows the relationship between this new metric and the binding rigidity in this setup, and a very good correlation (i.e., R 2 = 0.940), as well as the experimental evaluation of this assumption. In this evaluation, a sound-scattering object with l max = λ / 2 was placed at the origin (x, y, z) = (0, 0, 0), and the PAT was placed 12 cm above the object. For each of the four objects, we created a single trap in 2,000 random configurations (i.e., a total of 8,000 samples). Here, the x and y coordinates of the trap positions ranged from -5 cm to 5 cm, and z was set from 2 to 9 cm. Trap positions that were too close to the object (i.e., at a distance of less than 2λ) were excluded. The BASELINE solver was used to create the traps, and U j ' and ∇ 2 U j were calculated and plotted together (see Figure 6). The data obtained are

[0133]

Number

[0134] (b1 = -7.23×10 -7 and b2 = -1.69×10 -8) can be fitted with a straight line as shown, and the square of the correlation R 2 = 0.940. This correlation indicates that minimizing U j ' results in maximizing the binding rigidity (∇ 2 U j ).

[0135] The simplified indicator (Equation (4)) may not be directly used in setups where this assumption is not valid, but it should be noted that it can be easily adjusted to other setups such as the V - shape in Figure 1d. Here, it is shown how this technique can be adjusted to three other PAT setups, namely, the up - down type, the V - shape, and the single - sided type without a reflector. Here, it is assumed that the same 16×16 PAT is used, but note that the up - down type and the V - shape use two PATs. First, since sound waves propagating in the +z and -z directions from the upper and lower arrays can create a trap like a standing wave perpendicular to each other, the same simplification (i.e., Equation (4)) can be used for the up - down type setup.

[0136] In the V - shape setup of Figure 1d where there is an angle (φ = 90°) between the PATs, the propagation directions of the two PATs are (sinφ / 2, 0, cosφ / 2) and (-sinφ / 2, 0, cosφ / 2) respectively. Thus, thanks to the waves propagating in opposite directions along the x - axis, the lower indicator enables the creation of a strong floating trap. In other words, in this arrangement, the major axis is the x - axis.

[0137]

Number

[0138] Here, it should be noted that the constants K1 and K2 are determined by the physical properties of the particles and air (see Equation (4)). Therefore, K1 and K2 are constant values determined by the Gor'kov equation, rather than weights. The single-sided setup without a reflector is the most difficult among the three because there are no sound waves propagating in the opposite direction. However, by using the following indicators, a vortex trap very similar to what has already been shown can still be created.

[0139]

Number

[0140] The two-step scattering model works well in any floating setup. Therefore, by combining the two-step scattering model with a floating solver that uses appropriate indicators, floating traps can be created in top-bottom, V-shaped, and single-sided setups, even in the presence of sound scattering objects such as the shown sphere that may have a radius of 3 cm.

[0141] Instead of directly using the bound rigidity (∇ 2 U j ), this SIMPLIFIED solver uses the proposed simplified Gor'kov potential (U j ') as the target cost function. The differential of U

[0142]

Number

[0143] with respect to the phase of each transducer is calculated as j ', and can be calculated as

[0144]

Number

[0145] like this.

[0146] Here, Re[ ] and Im[ ] represent the real and imaginary parts, and p j,n represents the complex pressure value at the j-th trap position created by the single transducer n.

[0147] ∇ 2 U j and U j ' have a negative correlation, so the cost function

[0148]

Number

[0149] can be obtained as follows to maximize the binding rigidity.

[0150]

Number

[0151] The weight coefficient w s was fixed at 0.0001. The gradient of this cost function

[0152]

Number

[0153] can be calculated as follows.

[0154]

Number

[0155] Also in this case, calculating this gradient requires sampling only the pressure values at two points per trap, enabling fast calculation.

[0156] Any optimization algorithm such as BFGS for this cost function

[0157]

Number

[0158] Although it can be used to minimize, since the gradient descent method is suitable for parallel calculation, the gradient descent method was used. Further, for simplification, the step size of the gradient descent algorithm is

[0159]

Number

[0160] set to, which can be determined without using any line searching algorithm. For all evaluations in this application, the number of iterations (K) was set to K = 100, but it will be understood that this is adjustable.

[0161] Figure 7 summarizes the steps in implementing the simplified solver. In the first step S700, the desired location of the floating trap is determined. In the next step S702, the main axis of the transducer array is determined. In the above example, the main axis is the z-axis, but it will be understood that this is not always the case. It will also be understood that these steps can be performed in any order or simultaneously.

[0162] Then, in step S704, the pressure values at two points around each trap are sampled. These sampled values are used in step S706 to calculate, for example, the simplified Gor'kov potential using equation (d). And the pressure (i.e., the binding rigidity) at each floating trap may be maximized in step S708 using the simplified Gor'kov potential calculated in the cost function as shown, for example, in equation (f). Any suitable optimization algorithm can be used to minimize the cost function.

[0163] Evaluation of the Above-mentioned Technology In this section, how the above-mentioned technology was evaluated is explained. In the evaluation, four 3D models of scattering objects, flat, smooth, brick, and rabbit were used. These are shown in FIGS. 8a to 8d. The maximum length (l max ) of the mesh elements was always less than λ, λ / 2, λ / 4, or λ / 6, and a polygon mesh processing library was used to uniformly remesh the 3D model. The following table shows the number of mesh elements of each size.

[0164]

Table 1

[0165] The program detects edges with a dihedral angle greater than a specific degree as features of the object and retains those features while remeshing. l max = λ / 2 model is the mesh size with the best balance between speed and accuracy as shown below, so in most of the evaluations, l max = λ / 2 model was used.

[0166] Mesh size dependence of trap quality: FIG. 9 summarizes the computational performance of the above-mentioned technology. When the number of transducers (N) and the number of solver iterations (K) were fixed (i.e., N = 256, K = 100), how the number of traps (J) and the number of mesh elements (M) affected the calculation speed was evaluated. The results show a linear relationship between them as expected, and a high update rate exceeding 10,000 fps (i.e., a calculation time of less than 0.1 ms) can be achieved in some scenarios (e.g., M = ~8,000 and J = 4). For example, the rabbit and flat reflector used in the application of the four traps in FIG. 6 (i.e., 12×12 cm 2) The 3D model consists of a total of 4,134 elements and achieves over 15,000 fps. The plot also shows that even in the slowest scenario among the plots (i.e., J = 16 and M = 32,000), it is still possible to obtain over 700 fps, which is sufficient to manipulate particles in real time. The part related to the setup cannot be calculated in real time as described above, but this part can be pre-calculated once the setup is defined.

[0167] As shown in Fig. 9, the number of mesh elements (M) is an important parameter that greatly affects the calculation speed of the newly described technology. The total number of mesh elements depends on the mesh resolution of the 3D model (i.e., the size of the elements), and the mesh resolution also affects the accuracy of the BEM. In a normal scattering problem using BEM, usually six boundary elements per wavelength are required for an accurate scattering simulation. However, the purpose of this study is not to accurately simulate the sound field but to obtain the phase of the transducer that provides sufficient binding stiffness. Therefore, such a high degree of freedom per wavelength may not be necessary in the proposed scattering model.

[0168] To find the best-sized balance of mesh elements, the mesh size dependence of the trap quality (i.e., stiffness) was evaluated using 3D models with different maximum lengths of mesh elements (l max = {λ, λ / 2, λ / 4, λ / 6}). In this evaluation, a single trap was created using the BASELINE solver at the same trap positions as described in the Metric Validity test explained below, and the binding stiffness (∇ 2 U j ) was simulated using the finest mesh (i.e., λ / 6). Fig. 10 summarizes the average stiffness, and l max= The use of λ is insufficient for the two-step scattering model and cannot provide sufficient rigidity compared to the maximum element size smaller than the wavelength (e.g., especially with respect to smoothness and bricks). Considering the balance between speed and accuracy, for the rest of the evaluation, l max = It was decided to use λ / 2 to obtain the phase of the transducer.

[0169] Figure 10 shows the mesh size dependence of the trap quality in the proposed method. The maximum length of the mesh elements is l max = When λ, the model cannot provide sufficient rigidity at the trap position.

[0170] Computational performance: The computational performance of the proposed technique was evaluated using a consumer-grade laptop PC (2.60 GHz Intel Core i7-9750H CPU) equipped with a single GPU (NVIDIA GeForce RTX 2080), C++, and OpenCL. The positions of the traps and mesh elements were randomly generated for testing since the computation time does not depend on them. 100 tests were conducted for each combination of the number of traps (J = {1, 2, 4, 8, 16}) and the number of mesh elements (M = {1,000, 2,000, 4,000, 8,000, 16,000, 32,000}), and the average computation time was reported. In the implementation described here, the maximum number of frames (i.e., the activation of the transducer) that the GPU can compute simultaneously depends on the number of the GPU's workgroup size (i.e., in this case N w = 1,024), and the number of points required to compute each frame (i.e., L = 2J in the proposed simplified solver), and N wIt is determined like / 2J. Since this is directly related to the available update rate by selecting an index with a small L, it shows the importance of selecting an index with a small L. For example, using the proposed simplified index (L = 2J) enables the solver to calculate approximately twice as fast as using the original Gor'kov potential (L = 4J).

[0171] As described above, Figure 9 summarizes the overall computational performance of the proposed technique (i.e., the combination of the proposed two-step model and the proposed simplified solver, excluding pre-computation) for a given number of transducers (N = 256) and a given number of solver iterations (K = 100). Further, to show the breakdown of the computation time, how fast the proposed scattering model can be computed alone was tested (see Figure 11). In these plots, the solid line represents the computation time of the model only, and the dashed line represents the overall computation time (i.e., the same plot as Figure 9). These plots show that the solving process becomes more dominant when the number of traps (J) is larger. This is more prominent when the number of iterations (K) is larger (see Figure 12). The number of transducers (N) and the number of traps (N) are determined by the hardware and the application respectively, and thus cannot be changed. To reduce the overall computation time while maintaining sufficient accuracy, the number of mesh elements (M) and the number of iterations (K) become the key to balancing speed and accuracy, and these will be explored next.

[0172] In these performance evaluations, since the main focus is on how fast the application-related part can achieve the goal, the setup-related part (i.e., pre-computation) is excluded. Different from the application-related part, the computation time of the setup-related part depends not only on N, L, and M, but also on the geometry of the object. That is, even when two objects have the same number of mesh elements (M), the computation time of these objects may be different (for example, a flat reflector is easily solved). For reference, l maxThe pre - calculation of the 3D model with λ / 2 takes about 9s for flat, 12s for smooth, 21s for brick, and 17s for rabbit using a naive CPU implementation.

[0173] Convergence and initialization: Next, it is shown how well the SIMPLIFIED floating - point solver works for multi - point floating (i.e., the number of traps J={1, 2, 4, 8, 16}) in the presence of the four scattering objects used in the previous evaluation. 1,000 random combinations of trap positions were used for each condition. To avoid cases where the traps are too close to each other, the minimum distance between traps was set to 2λ. Figure 13a shows the average stiffness and their standard deviations for different numbers of traps (J) with K = {10, 20, 40, 80, 100, 200, 400, 800}, and shows the increase in stiffness with iteration when the phase of the transducer is randomly initialized. Even with the maximum number of traps (i.e., J = 16), after several iterations, the positive stiffness required to bind the particles can be achieved.

[0174] Figure 13b shows the results when the phase obtained using the modified HEURISTIC solver is used instead of a random initial phase. The plot shows that the use of such a HEURISTIC initial guess reduces the number of iterations (K) required by the SIMPLIFIED solver. Here, note that even though the HEURISTIC solver already provides a relatively high average stiffness without iteration (i.e., K = 0), iterations are still required to reduce the standard deviation. This is because in multi - point acoustic levitation, weak traps may not be able to hold the particles in the air, and the goal is to generate equally strong traps (see further discussion in the next section). The advantage of using the modified version of the HEURISTIC solver is that it starts from exactly the same points as the SIMPLIFIED solver (i.e., the trap positions (x j , y j , z j ) and its slightly above position (xj , y j , z j The pressure value of + h)) is used, and thus the same transfer matrix can be used for both this initial step and the iterative step without any additional modeling process being required. Along these results, this HEURISTIC initialization and K = 100 were used for all applications and for other evaluations.

[0175] In other words, Figure 13a shows that the trap quality of the SIMPLIFIED solver improves as a function of the number of iterations (K). Figure 13b shows that using the modified HEURISTIC solver as an initial guess reduces the number of iterations required to be converged. Error bars represent the standard deviation.

[0176] Comparison between solvers: As described above, three solvers, namely, BASELINE, HEURISTIC, and SIMPLIFIED were considered. Figures 14a and 14b evaluate the three solvers in terms of the binding stiffness. As in the previous evaluation, 1,000 random combinations of trap positions were used for each condition (i.e., four scattering objects with different numbers of traps J = {1, 2, 4, 8, 16}). The number of transducers (N = 256) and the number of iterations (K = 100) were fixed. Bars represent the standard deviation. The symbol "<" indicates that there is a significant difference between homogeneous groups represented by the symbol "{}".

[0177] The BASELINE solver is a physically accurate and widely accepted metric of binding quality (i.e., binding stiffness ∇ 2 U j) is used, which is a conventional method with low speed. The HEURISTIC solver is an extension of the HAE framework described in detail by Marzo et al. in "Holographic acoustic tweezers" published in Proc. Natl. Acad. Sci. 201813047 (2018). By introducing a π-radian offset to the target phase, it enables trap creation by creating two foci around each trap. The HAE framework simplifies the calculation of levitation traps by encoding them as a combination of a holographic acoustic lens that creates the foci and a determined levitation signature. This framework supports a wide range of symmetric transducer arrangements (e.g., single-sided, top-bottom, V-shaped) and has been extended to multi-point levitation. This method is fast and will work well for single-point operations, but there is a high possibility of destructive interference between traps in multi-point operations. Figures 14a and 14b show that only the SIMPLIFIED solver provides both high computational speed and trap quality. Here, it should be noted that all three solvers use the two-step scattering model described in the process of Figure 2.

[0178] Figure 14a shows the average binding stiffness and its standard deviation obtained by different solvers. The average values indicate that the BASELINE is slightly better than the HEURISTIC overall, and the performance of the SIMPLIFIED tends to be intermediate between these two. This relationship was also statistically confirmed using statistical software (IBM SPSS Statistics 25) as shown in Figure 14a. The plot also shows that the SIMPLIFIED provides the smallest standard deviation among the solvers. Providing a small standard deviation is important in multi-point acoustic levitation to avoid weak traps and achieve stable particle manipulation.

[0179] To emphasize this point, the same evaluation was performed, focusing on the weakest trap among the J traps (see Figure 14b). The plot makes the difference between HEURISTIC and the other two more evident, indicating that HEURISTIC is likely unable to create traps when the number of traps is large (i.e., negative stiffness at J = 16). This is why HEURISTIC is not sufficient even if it provides the fastest computational performance. Figure 14b also shows that SIMPLIFIED is slightly better than BASELINE in terms of minimum stiffness, indicating that SIMPLIFIED is more suitable for uniformly providing sufficient stiffness to all traps in multi-point floating.

[0180] In other words, as shown in Figures 14a and 14b, the SIMPLIFIED solver avoids destructive interference between multiple traps compared to the HEURISTIC solver while achieving a similar quality (i.e., binding stiffness) to the BASELINE solver that directly maximizes ∇ 2 U j . Furthermore, with appropriate initialization, the SIMPLIFIED solver can converge within 100 iterations. Therefore, the proposed simplified solver represents the best-balanced solution for achieving accurate and fast acoustic operations.

[0181] Sound field distortion and correction: To show how the sound field is distorted by sound-scattering objects and how they are corrected by the proposed two-step BEM model, attempts were made to create four traps without using and using the above process. Figures 15a to 15f show the evaluation results. In Figures 15a to 15c, the arrangement is as shown in Figure 8a (i.e., smooth - without sound-scattering objects), and in Figures 15d to 15f, the arrangement is as shown in Figure 8c (i.e., bricks as sound-scattering objects).

[0182] Figure 15a plots the acoustic field simulation when creating four traps in a smooth configuration using standard methods of imaging technology, as described, for example, in "Automatic contactless injection, transportation, merging and ejection of droplets with a multifocal point acoustic levitator" by Andrade et al. published in Rev. Sci. Instrum 89 (2018). Such methods can be calculated by assuming that sound waves scattered from a flat reflector are emitted from virtual sound sources located at the mirror image positions of the actual sound sources (i.e., transducers). Such simulations do not take into account the scattering of sound from objects (i.e., assume the presence of only a flat reflector), and thus the generated acoustic field can be distorted due to ignoring the presence of objects. Figure 15b plots the acoustic field simulation when creating four traps in a smooth configuration using the above-described technique. Then, the trap positions are rotated horizontally, and the binding rigidity at the points of the four traps corresponding to the rotation angle is plotted as shown in Figure 15c (∇ 2 U j ). The shaded regions represent the minimum and maximum binding rigidities in each case.

[0183] Similarly, Figure 15d plots the acoustic field simulation when creating four traps in a further "brick" configuration using standard methods of imaging technology. Figure 15e plots the acoustic field simulation when creating four traps in a smooth configuration using the above-described technique. Figure 15f plots the binding rigidity corresponding to the rotation angle.

[0184] Figures 15a through 15f show that the sound field is greatly distorted by smooth and brick-shaped objects (e.g., the average confinement stiffness decreases by an average of 77% and 75%, respectively). The brick-shaped object is more difficult because it has a non-smooth surface. The minimum confinement stiffness in the brick using known techniques becomes even negative (especially in Figure 15f), suggesting that at least one of the four traps cannot levitate particles (e.g., the trap in the lower right of the uncorrected image in Figure 15d). On the other hand, the new two-step scattering model can correct such distortion and improve the confinement stiffness by considering the scattering of sound from the object.

[0185] Application ability of new technology The combination of the two-step scattering model and the simplified levitation solver enables real-time manipulation of materials in 3D space in the presence of sound-scattering objects. Examples of such uses and applications are shown in Figures 16a through 16d. All applications used the same levitation setup. The applications were created using a single PAT of a 16×16 transducer designed as an extension of the Ultraino platform modified for higher communication speeds. The array used Murata's MA40S4S transducer (40 kHz, 10.5 mm diameter (~1.2λ), delivering ~8.1 Pa at a distance of 1 m when driven at 20 Vpp). A Waveshare CoreEP4CE10 field-programmable gate array (FPGA) board was used to receive phase and amplitude updates from the CPU, using an 8 Mbyte / sec USB FT245 asynchronous FIFO interface, enabling more than 10,000 phase and amplitude updates per second. The PAT and the underlying flat acrylic reflector were aligned to overlap each other with an adjustable spacing (e.g., fixed at 12 cm in this study). The square portion of the flat reflector (12×12 cm 2) can be replaced by any scattering surface such as a 3D-printed object, a set of bricks, and a glass container filled with water. We used a LulzBot Mini 3D printer with eSUN's PLA+ filament to 3D-print the objects.

[0186] In Fig. 16a, a plurality (e.g., 10) of expanded polystyrene (EPS) particles are levitated on a 3D-printed smooth surface. The simulated sound field may be plotted on the xy plane at λ / 4 above the trap positions, showing 10 high-pressure points. This application demonstrates that acoustic 3D manipulation is possible even with non-flat reflectors.

[0187] Fig. 16b shows that particles can be levitated even under an obstacle that scatters sound. Usually, such an obstacle blocks most of the direct sound contribution from the transducer. Thus, Fig. 16b demonstrates the operating ability in scenarios that were previously impossible. The simulated sound field may be plotted on the xz plane to show the trap positions.

[0188] Unlike other levitation techniques such as electromagnetics, the acoustic approach can levitate almost all types of substances, including solids and liquids. Fig. 16c shows the manipulation of water droplets in the presence of a 3D-printed cactus. Acoustic manipulation of droplets is particularly difficult because the acoustic velocity of air particles at the trap position needs to be carefully adjusted to keep it within a range determined by the droplet radius and surface tension to avoid atomization of the droplets. The high computational speed of the proposed technique enables real-time estimation of the acoustic velocity and dynamically adjusts the amplitude of the transducer to keep the acoustic velocity constant along the operation path (see Fig. 16e). In addition, by modulating the amplitude of all transducers at a specific frequency, it is possible to induce oscillatory vibrations in the levitated droplets, which is useful for mixing multiple materials non-contact without causing any cross-contamination.

[0189] Furthermore, the two - scattering model works well even if the scattering surface is liquid, as shown in Fig. 16d. It is possible to manipulate a mixture of water droplets while taking into account the liquid surface of a container filled with water. The liquid surface is approximated as acoustically rigid (i.e., β m = 0), and continues to show correct droplet manipulation. Such material independence brings versatility to the proposed technology, which can be applied in fields such as computational manufacturing, lab - on - a - chip, and biomedical imaging. As detailed above, the use of other β m values is also possible.

[0190] In the acoustic manipulation of droplets, the ratio of the acoustic force to the surface force of a levitated droplet is described by the acoustic Bond number

[0191]

Equation

[0192] , where σ is the surface tension of the liquid, R s is the droplet radius, and v rms is the root - mean - square of the acoustic velocity of air particles (see, for example, "Acoustophoretic contactless transport and handling of matter in air" by Foresti et al. published in Proc Natl. Acad. Sci 110 (2013)). To avoid atomization (i.e., droplet breakup) of the levitated droplets, this acoustic Bond number needs to be between 2.5 and 3.6, as determined experimentally. Therefore, it is important to keep the acoustic velocity constant along the manipulation path. The high computational speed of the proposed technology enables real - time estimation of the acoustic velocity and adjustment of the transducer amplitude to keep the acoustic velocity constant along the manipulation path (see Fig. 16e).

[0193] Creation of POV images using a high update rate An important aspect of the proposed method is its computational speed (see Figures 9 and 12). The high update rate of the transducer's phased array (PAT), ideally exceeding 10,000 fps, enables the manipulation of high-speed expanded polystyrene (EPS) particles (i.e., a maximum speed of 8.75 m / s was reported in an up-down setup). This allows for the creation of volumetric images in air (i.e., acoustic holography), even in the presence of sound-scattering objects, by using the persistence of vision (POV) effect achieved by scanning the particles in 0.1 s. Furthermore, due to the high update rate of the new technology, the created POV images can be interactive with user input (e.g., keyboard, hand gestures) without any noticeable delay. For example, a LeapMotion sensor may be used to detect the position of the user's fingertips. Thus, we have the ability to create composite reality applications that utilize acoustic holography.

[0194] This technology enables, for example, the creation of a butterfly flying around a 3D-printed rabbit (M = 4,134) controlled by hand gestures by using a single particle colored by a full-color LED. Other examples of volumetric shapes, such as two particles above a plastic brick (M = 5,010) or a single particle below a sound-scattering obstacle (M = 3,792), are also possible. These are the first demonstrations of the creation of digital volumetric images with physical objects as a new MR human-computer interface that blurs the boundary between the digital and physical worlds.

[0195] However, since the particles must scan all geometric shapes within the POV time (i.e., 0.1 seconds), the volumetric geometric shapes that can be created by such point scanning techniques are limited to simple shapes. Therefore, furthermore, a display based on surface scanning of free space for creating more complex volumetric shapes with many voxels (volume elements) was first demonstrated. In this technique shown in Figure 1b, a lightweight fabric was levitated in the same floating setup used for the point scanning technique, using a high-speed projector (i.e., 1,440 fps) and a mirror.

[0196] The lightweight fabric screen may be prepared in any suitable manner. For example, we first laser cut a lightweight and acoustically transparent fabric (Super Organza) into 3×3 cm 2 squares. Four EPS particles were adhered to the fabric and served as anchors to enable 6-degree-of-freedom manipulation of the fabric. For projection mapping onto this levitated fabric, we used a projector (Texas Instruments, DLP LightCrafter Evaluation Module) with a native resolution of 608×684 pixels. We pre-acquired the internal parameters of this projector by using the functions of OpenCV together with a checkerboard and a webcam, and then obtained the external parameters (i.e., the position and orientation with respect to the coordinates of the levitation device) by using the manually collected combinations of the trap positions within the coordinates of the levitation device and the pixel positions within the coordinates of the projector. Then, we used such parameters in our OpenGL camera (projection matrix and view matrix) to enable real-time projection mapping.

[0197] In the application of volume display based on butterfly-like point scanning, we used a high-brightness full-color LED (OptoSupply, OSTCWBTHC1S) to illuminate the levitated EPS particles. The LED was directly controlled by the FPGA, which also controlled the transducer, resulting in synchronization between the color of the illumination and the movement of the levitated particles. All scan paths were generated to be scanned by the particles in the POV time (i.e., 0.1 s). Thus, we were able to create a volume POV image. Here, in the point-scanning-based method, the maximum number of voxels (N v ) is determined as N v = J·f l / f POV by the update rate (f l ) of the levitation device, the number of traps (J), and the POV rate (f POV = 10 Hz) (for example, when f l = 10,000 and J = 4, N v = 4,000). Note that since the paths created by these voxels need to be scanned by single or multiple points, there are further constraints on the voxel arrangement. That is, the voxels need to be continuous, and the movement of the particles along the voxel path needs to be within the range of the system's capabilities (i.e., maximum speed and acceleration). These constraints make it difficult to create complex volume shapes with point-scanning-based methods.

[0198] Regarding the volume display based on the surface scan as shown in Fig. 1b, we reused the same fabric, projector, and calibration method used in the air screen application. However, in this application, we used the projector in the high-speed binary mode of 1,440 fps. We placed a mirror in the system to cover the angle when the projector cannot project directly onto the fabric (when the fabric and the projection direction are parallel). In other words, we used the mirror as a second projector. We created 144 cross-sectional binary images of a 3D model (i.e., a rabbit) at every 1.25 degrees, mapped those images onto a rotating screen, and encoded them into 24-bit images as in (46). Then, the system floated the fabric and rotated it at 5 revolutions per second while updating the images encoded at 60 Hz. Our OpenGL-based software can adjust the timing of projecting the cross-sectional images to match the timing of the rotation of the fabric. Also, the software receives the VSYNC signal and automatically adjusts the rotation speed of the fabric corresponding to the update of the projector.

[0199] In the surface-based approach, the maximum number of voxels (N v ) of the created image is determined by the update rate (f p ) of the projector and the number of pixels (N p ) of the projected 2D image as N v = N p ·f p / f POV as follows. Therefore, ideally, N v= 608×684×1,440 / 10 ≒ 60,000,000, which is approximately 15,000 times larger than the point - based method. Assuming full utilization of pixels with a stationary projector like our current system is not realistic, but by using a projection engine with a rotating mirror, it is possible to increase the pixel utilization rate to almost 100%. Furthermore, since the voxel arrangement is fixed and does not depend on the content, if the floating device can rotate the fabric at 5 revolutions per second, the displayed content does not need to consider the capabilities (i.e., speed and acceleration) of the floating device.

[0200] The proposed technology can rotate the fabric at 5 revolutions per second in the presence of a sound - scattering object while synchronously projecting cross - sectional images of a 3D model onto the rotating fabric, showing a complete volume image in the air due to the POV effect. The reason the mirror is used is to project the image even when the projection direction and the fabric are parallel. Two photos taken from different viewpoints show a digital 3D image of a rabbit projected onto the rotating fabric. The digital 3D image was created on top of a physical rabbit (M = 4,134) 3D - printed using the exact same 3D model of the digital rabbit. Therefore, the proposed technology can project complex volume shapes in the air, which can be viewed from any direction.

[0201] Sound field simulation Similar to the normal BEM, the main purpose for which we developed our model in this study was to solve for the transducer activation (τ) for creating multiple traps at high speed, but our model can be used for the broader purpose of simulating the sound field. The sound field simulated by the normal BEM and the sound field generated by our model are equivalent to the sound field simulated by the normal BEM.

[0202] Conventional BEM requires solving the linear equation Ap = b every time to simulate the sound field with different τ even with the same setup. In contrast, in our model, once we calculate the transfer matrix (E), that transfer matrix (E) can be used to simulate the sound field with different τ as long as the same setup is not used. Here, it should be noted that E can be calculated very quickly once we obtain the data from the pre-calculation (i.e., the matrix H). Our model is particularly useful for simulating and evaluating the sound field multiple times with different τ but the same setup.

[0203] Lifting of particles A potential limitation of the proposed technique is the manipulation of particles near the scattering surface. When we try to create a trap near the surface, the strong sound reflection from the surface tends to create a sound field like a standing wave on the surface, resulting in creating traps at specific discrete heights (z) from the surface (i.e., z = λ / 4, 3λ / 4). Therefore, it is difficult to manipulate the particles from z = λ / 4 to z = 3λ / 4 or vice versa. Figure 17a shows an experimental setup to illustrate this limitation. We tried to create a single trap at a specific height (z) from a flat surface with our solver.

[0204] Figure 17b plots how far the simulated trap positions (i.e., the positions where the Gor'kov potential is minimum) deviated from the target trap positions for both the BASELINE solver and the SIMPLIFIED solver. The plot shows a very large position error in the region of λ / 2 < z < 3λ / 4, indicating that we failed to create a trap in this region. In this study, the wavelength was λ = 8.65 mm. The difficulty of this operation near the scattering surface was also confirmed experimentally. To realize an acoustic holographic system with this feature, further research efforts are required in both the algorithm area and the hardware area (e.g., transducer placement).

[0205] Figure 17c shows a practical way to avoid this problem by using a prop that scatters sound in the form of a wedge-shaped ramp in this example. Our two-step scattering model allows us to manipulate particles along the ramp by creating traps on the λ / 4 of the ramp's surface. When the particles are high enough from the surface (e.g., z≧3λ / 4 = 6.49mm), we can move the particles away from the ramp and manipulate them in 3D without constraints. We have experimentally confirmed that this technique works well for lifting particles from the ground.

[0206] Movement or change of sound scattering object In our scattering model, the mesh model does not change over time. This assumption allows us to pre-compute the scattering model (i.e., matrix H). In other words, it is difficult to handle dynamic scattering objects with the techniques proposed above. If we know in advance the nature of the dynamic evolution of the sound scattering object, different H matrices can be pre-computed, and the other two matrices F and G can be computed in real time. If the sound scattering object changes in such a way that it cannot be predicted in advance, we need to repeatedly solve the linear equation Ap (n) = b (n) for all N transducers to compute H in real time, where A is an M×M matrix. Note here that, as shown in equation (f), matrix A depends only on the geometry of the scattering object and not on the position of the transducer or trap.

[0207] One common scenario is when the shape of the scattering object is constant, but the position of the object or the arrangement of the transducers changes. In such a scenario, we can assume that the object is relatively static by assuming instead that the position of the transducer changes. Thus, matrix A remains constant even while the actual position of the object is moving. Therefore, if we decompose this matrix (e.g., using LU decomposition), we can reuse the decomposed matrix and easily solve the linear equation to obtain different Hs at high speed during the movement of the object.

[0208] Summary Prior to this proposed new approach, 3D manipulation of materials using acoustic holography had only been achieved within an empty volume. This limitation has hitherto forced the technology to be used in limited scenarios (i.e., where there are no scattering objects in the vicinity). Here, this limitation is overcome by reformulating and simplifying the model and solver for acoustic holography. The proposed approach expands the possibilities of acoustic levitation, enables 3D printing for non-contact manufacturing, and enables the mixing of physical and digital artifacts for new MR applications. This application assumes that only sound scattering objects (e.g., plastic, water) with a high acoustic impedance compared to air are present within a single propagation medium (i.e., air). However, BEM can also be used to calculate the scattering of sound from a sound soft boundary even through multiple media. The same two-step approach can potentially be applied to such more complex scenarios, increasing the calculation speed and opening the way to real-time utilization beyond the demonstrated environment.

[0209] Also, when the current limitation of using only static scattering objects (i.e., a single pre-computed matrix H) is relaxed, this range of potential scenarios expands, but so do the issues that need to be considered. That is, by eliminating the need for an empty volume, this method already enables ultrasonic-based solutions to be applied to many more realistic environments, such as inside household appliances or inside a car dashboard.

[0210] The obvious step for supporting dynamic (i.e., moving / changing) objects is to pre-compute different H matrices, one for each state of the object. This requires prior knowledge of the nature of the dynamic evolution of the object, but even this simple step is sufficient to enable many new applications such as 3D printing and non-contact assembly, since in all these cases the geometric evolution is known in advance.

[0211] Moving on to a fully interactive scenario opens up new challenges and possibilities. Regarding objects with fixed shapes that interactively change their position and orientation, the LU decomposition technique discussed above may make it possible for the matrix H to be computed in real time (e.g., >60fps). The most difficult scenarios are when objects change their shape, position, and orientation in unpredictable ways (e.g., an MR application where the user's hand interacts within the working volume). A new approach is needed here to compute H in real time, and one potential solution is to utilize the local nature of the changes. That is, if the position and / or geometric shape of the object do not change drastically during the update, the solution for the previous geometric shape can be used to reduce the computational cost and as a good initial estimate for the next geometric shape. It is worth noting that the computational speed of this setup part does not need to achieve 10,000fps and may be sufficient at a more normal speed (e.g., >30fps).

[0212] Also, the new two-step scattering model can be adapted without modification to various PAT arrangements (see, e.g., single-sided type, V-shaped type, up-and-down type, Figures 1b, 1c, 1d). This provides great flexibility in the design of new applications using the proposed technique. However, the simplified metrics should not be used by default as in Equation (4), but rather should be adjusted according to the geometric relationship between the relevant PAT and the trap positions. The option to modify the simplified metrics according to the experimental setup and the content to be displayed provides additional adjustment parameters for an acoustic holographic system to achieve both optimal speed and accuracy. This suggests that dynamically selecting the simplification of the metrics most suitable for the setup and content used can extract the highest accuracy and speed from the device.

[0213] Point scan-based techniques have been employed and explored to realize a free-space volume display by using several levitation techniques such as acoustic traps, optical tweezer traps, and electromagnetic traps. In this application, for the first time, a surface scan-based technique is introduced into these levitation techniques, realizing a free-space volume display that can represent more voxels with minimal constraints on voxel placement compared to point-based ones. Compared with volume displays that use mechanically rotatable screens or emitters, the advantage of the proposed technique is that the rotating screen itself can be operated within a space directly accessible to the user, highlighting the MR aspect of the acoustic holographic technology proposed in this application.

[0214] Those skilled in the art have described above what is considered to be the best mode and other modes of practicing the technology as appropriate, but will understand that the technology should not be limited to the specific configurations and methods disclosed in this description of the preferred embodiments. Those skilled in the art will recognize that the technology has broad applications and that embodiments may be subject to extensive modifications without departing from any inventive concept defined in the appended claims.

Explanation of Reference Numerals

[0215] 100 System 110 Control Device 112 Processor 114 Memory 116 Input / Output Interface 118 Two-Step Scattering Model 119 Levitation Solver 120 Acoustic Chamber 122 Transducer Array 124 Target Object 126 Scattering Object 220 Acoustic Chamber 222 Transducer Array 224 Target Object 226 Scattering Object 320 Acoustic Chamber 322 Transducer Array 326 Scattering Object 420 Acoustic Chamber 422 Transducer Array 426 Scattering Object

Claims

1. A computer-implemented method for controlling the location of a target object within an acoustic volume using an array of sound-generating transducers, wherein the acoustic volume includes a scattering object, and the method A step of obtaining a static matrix (H) representing the contribution of each transducer in the array of transducers to each of a plurality of locations on the scattering object, wherein the static matrix does not change when controlling the locations on the target object. The steps include defining a plurality of control points within the aforementioned acoustic volume, The steps include: calculating in real time a direct transfer matrix (F) representing the direct contribution from each transducer in the array of transducers to each of the plurality of control points; The steps include: calculating in real time a scattering transfer matrix (G) representing the contribution of scattering from the plurality of locations on the scattering object to each of the plurality of control points; A step of determining in real time an extended transfer matrix (E) representing the direct and scattered contributions from each transducer of the array of transducers to each of the plurality of control points, wherein the extended transfer matrix is ​​determined from E = F + GH using the static matrix, the direct transfer matrix, and the scattered transfer matrix, A method comprising the steps of using the extended transfer matrix to determine a control command for each transducer in the array of transducers to generate an acoustic trap at at least one of the plurality of control points, wherein the acoustic trap is configured to restrain the target object in order to control the location of the target object within the acoustic volume.

2. The aforementioned static matrix, in the setup phase, Defining multiple locations on the scattering object within the aforementioned acoustic volume, To obtain location information for each of the aforementioned multiple locations, To obtain positional information for each transducer in the array of transducers, For each of the aforementioned multiple locations, the set of sound pressure contributions from each transducer in the transducer array is calculated using the location information and the position information. The method according to claim 1, wherein the calculation is performed by storing each set of sound pressure contributions in the static matrix.

3. The method according to claim 1, wherein the plurality of locations are a plurality of mesh elements.

4. The method according to claim 3, wherein each mesh element has a maximum length of λ / 2, where λ is the wavelength of the sound generated by each transducer.

5. The scattering object changes over time, and the method A step of obtaining a plurality of static matrices, each static matrix representing the contribution of each transducer in the array of transducers to each of a plurality of locations on the scattering object at a particular time step, The method according to claim 1, further comprising the step of determining the extended transfer matrix using the static matrix of the particular time step for each time step.

6. The method according to claim 1, wherein the step of determining a control command includes optimizing the phase of each transducer in the array of transducers to maximize the binding stiffness at each location of the acoustic trap.

7. Optimizing the phase of each transducer The position of each acoustic trap is defined using the aforementioned multiple control points, Determining the principal axis of the transducer array, Sampling sound pressure values ​​at two locations along the main axis around each position of the acoustic trap, Using these sampled sound pressures, we can calculate an index of the binding stiffness, The method according to claim 6, comprising maximizing the calculated index of the constraint stiffness using a cost function.

8. A computer-implemented method for controlling the location of a target object within an acoustic volume using an array of sound-generating transducers, wherein the acoustic volume includes a scattering object, and the method A step of defining a plurality of control points within the acoustic volume, A step of determining an extended transfer matrix representing the direct contribution from each transducer in the array of transducers to each of the plurality of control points, and the scattering contribution from each transducer to each of the plurality of control points via the scattering object. Using the extended transfer matrix, control commands are issued to each transducer in the transducer array to generate an acoustic trap at at least one of the plurality of control points. The position of each acoustic trap is defined using the aforementioned multiple control points. Determining the principal axis of the array of transducers, Sampling sound pressure values ​​at two locations along the main axis around each position of the acoustic trap, These sampled sound pressures are used to estimate an index of binding stiffness, and A method comprising the step of determining by maximizing an estimated index of the binding stiffness using a cost function, wherein the acoustic trap is configured to bind the target object to control the location of the target object within the acoustic volume.

9. The aforementioned constraint stiffness index is the simplified Gor'kov index U j 'and, [Math 1] It is defined as follows, where V represents the volume of the target object, ω represents the angular frequency of the target object, c and ρ represent the speed of sound and density, and the subscripts 0 and p refer to the host medium (i.e., air) and particulate material, respectively, p j The method according to claim 7 or 8, wherein z represents the sound pressure at the control point from the j-th transducer, and z is the principal axis.

10. The aforementioned cost function is, [Math 2] Defined as, w s is the weight coefficient, bar (￣) represents the average value of all traps, J is the number of traps, and U j ' is a simplified Gor'kov index, [Math 3] The method according to claim 9, wherein the phase of each transducer in the array is the phase of the array.

11. The method according to claim 1 or claim 8, wherein the number of traps generated is in the range of 1 to 16.

12. A step of controlling the first location of a plurality of printing droplets in order to change the state of each of the plurality of printing droplets from liquid to solid, The process includes the step of controlling a second location for each of a plurality of solid printed objects in order to deposit each of a plurality of solid printed droplets at a desired location, A printing method comprising the steps of controlling the first location and controlling the second location using the method described in claim 1 or claim 8.

13. A method for generating volumetric motion images, The steps of providing multiple particles or a projection screen supported by multiple particles, A method comprising the steps of controlling the location of each of the plurality of particles as a target object using the method of claim 1 or claim 8, wherein a volumetric motion image is generated by the movement of the plurality of particles or the movement of a screen.

14. A non-transient data carrier that carries code causing the processor to perform the method according to claim 1 or claim 8 when executed on the processor.

15. An array of transducers configured to generate sound pressure, An acoustic volume defined by the sound pressure generated by the array of transducers, in which the location of a target object can be controlled, A processor is included, and the processor is Obtaining a static matrix (H) representing the contribution of each transducer in the array of transducers to each of a plurality of locations on a scattering object, wherein the static matrix does not change when controlling the locations on the target object. Define multiple control points within the aforementioned acoustic volume, Calculating in real time a direct transfer matrix (F) representing the direct contribution from each transducer in the transducer array to each of the multiple control points, Calculating in real time a scattering transfer matrix (G) representing the contribution of scattering from the plurality of locations on the scattering object to each of the plurality of control points, Determining in real time an extended transfer matrix (E) representing the direct and scattered contributions from each transducer of the array of transducers to each of the plurality of control points, wherein the extended transfer matrix is ​​determined from E = F + GH using the static matrix, the direct transfer matrix, and the scattered transfer matrix, and An apparatus configured to determine, using the extended transfer matrix, a control command to each transducer in the array of transducers for generating an acoustic trap at at least one of the plurality of control points, wherein the acoustic trap is configured to bind the target object to control the location of the target object within the acoustic volume.

16. An array of transducers configured to generate sound pressure, An acoustic volume defined by the sound pressure generated by the array of transducers, in which the location of a target object can be controlled, A processor is included, and the processor is Define multiple control points within the aforementioned acoustic volume, Determining an extended transfer matrix that represents the direct contribution from each transducer in the transducer array to each of the plurality of control points, and the scattering contribution from each transducer to each of the plurality of control points via a scattering object. Using the extended transfer matrix, control commands are issued to each transducer in the transducer array to generate an acoustic trap at at least one of the plurality of control points. The position of each acoustic trap is defined using the aforementioned multiple control points. Determining the principal axis of the array of transducers, Sampling sound pressure values ​​at two locations along the main axis around each position of the acoustic trap, These sampled sound pressures are used to estimate an index of binding stiffness, and An apparatus configured to make a determination, which is determined by maximizing an estimated index of the binding stiffness using a cost function, wherein the acoustic trap is configured to bind the target object to control the location of the target object within the acoustic volume.

Images