Ultraleap Patent | Ultrasound acoustic field manipulation techniques
Patent: Ultrasound acoustic field manipulation techniques
Drawings: Click to check drawins
Publication Number: 20220252550
Publication Date: 20220811
Applicants: Ultraleap
Abstract
Two methods for rendering haptic surfaces through spatial modulation of ultrasound are described: Singular Value Decomposition (SVD) with Tikhonov Regularization and Mini-Batch, Stochastic Gradient Descent with Momentum (MSGDM). Further, adjustment in the placement of transducers, arrays can be generated which allow for variable focus points using a limited set of driving signals. This is accomplished by generating a set of mutually-exclusive transducer placement layouts which when driven together produce a set of sufficiently orthogonal fields at various points of interest. Further, using basic transformations on sampling locations a field can be translated or rotated in 3-dimensional space. This enables the creation of spatiotemporal modulated haptics with reduced computation. Further, mid-air haptics using ultrasound can be generated via amplitude modulation or via spatial modulation of control points defined at spatial locations in an acoustic field above a phased array device.
Claims
1. An apparatus comprising: Singular Value Decomposition (SVD) with Tikhonov Regularization and Mini-Batch, Stochastic Gradient Descent with Momentum (MSGDM).
2. An method comprising: adjusting in the placement of transducers, generating arrays which allow for variable focus points using a limited set of driving signals.
Description
PRIOR APPLICATIONS
[0001] This application claims the benefit of:
[0002] (1) U.S. Provisional Patent Application No. 63/141,897, filed Jan. 26, 2021;
[0003] (2) U.S. Provisional Patent Application No. 63/156,829, filed on Mar. 4, 2021;
[0004] (3) U.S. Provisional Patent Application No. 63/167,855, filed on Mar. 30, 2021; and
[0005] (4) U.S. Provisional Patent Application No. 63/266,972, filed on Jan. 20, 2022.
[0006] All of these applications are incorporated by reference in their entirety.
FIELD OF THE DISCLOSURE
[0007] The present disclosure relates generally to improved manipulation techniques in acoustic transducer structures used in mid-air haptic systems.
BACKGROUND
[0008] Ultrasonic phased arrays can be used to create arbitrary acoustic fields. These can be used for haptic feedback, parametric audio, acoustic trapping, etc. The rendering of haptic surfaces is desirable for various applications areas such as; virtual reality, augmented reality and gesture control. The current method for rendering haptic surfaces involves the spatiotemporal modulation of single or multiple focal points. The problem we are trying to solve is the development of alternative methods for rendering haptic surfaces.
[0009] Gavrilov et al. pioneered early research into the possibility of using focused ultrasound as a non-invasive method for stimulating nerve structures in humans. Gavrilov, Leonid R., et al. "A study of reception with the use of focused ultrasound. I. Effects on the skin and deep receptor structures in man." Brain research 135.2 (1977): 265-277; Gavrilov, Leonid R., et al. "A study of reception with the use of focused ultrasound. II. Effects on the animal receptor structures." Brain research 135.2 (1977): 279-285; Gavrilov, L. R. "Use of focused ultrasound for stimulation of nerve structures." Ultrasonics 22.3 (1984): 132-138.
[0010] Gavrilov et al. went on to describe the mechanism by which focused ultrasound could be used to elicit tactile sensations. The non-linear, acoustic radiation force of focused ultrasound induces a shear wave in the skin. This shear wave produces an associated displacement that triggers certain mechanoreceptors in the skin. Gavrilov, L., and E. Tsirulnikov. "Mechanisms of stimulation effects of focused ultrasound on neural structures: Role of nonlinear effects." Nonlinear Acoust. at the Beginning of the 21st Cent (2002): 445-448.
[0011] Hoshi et al. developed a prototype two dimensional, ultrasonic phased-array capable of translating a single high intensity, focal point in a volume above the array. This system enabled users to interact with virtual objects through haptic feedback. Algorithms for producing multiple high intensity focal points using ultrasonic phased arrays had already been developed in the medical ultrasound community. Hoshi, Takayuki, et al. "Noncontact tactile display based on radiation pressure of airborne ultrasound." IEEE Transactions on Haptics 3.3 (2010): 155-165.
[0012] In particular, the conjugate field method of Ibinni et al. and the pseudo-inverse method of Ebbini et al. Ibbini, Mohammed S., and Charles A. Cain. "A field conjugation method for direct synthesis of hyperthermia phases-array heating patterns." IEEE transactions on ultrasonics, ferroelectrics, and frequency control 36.1 (1989): 3-9. Ebbini, Emad S., and Charles A. Cain. "Multiple-focus ultrasound phased-array pattern synthesis: optimal driving-signal distributions for hyperthermia." IEEE transactions on ultrasonics, ferroelectrics, and frequency control 36.5 (1989): 540-548.
[0013] Drawing on research from the field of computer-generated holography. Hertzberg et al. demonstrated that the iterative, weighted Gerchberg-Saxton (GSW) algorithm could be adapted to efficiently produce multiple high intensity focal points in 2D. Hertzberg, Yoni, et al. "Towards multifocal ultrasonic neural stimulation: pattern generation algorithms." Journal of neural engineering 7.5 (2010): 056002.
[0014] In an attempt to develop a more stable and more efficient algorithm, Long et al. extended the pseudo-inverse method with a power iteration step. The power iteration step aims to find the set of phases that should be assigned to the focal points in order to maximize the constructive interference and minimize the destructive interference. Long, Benjamin. et al. "Rendering volumetric haptic shapes in mid-air using ultrasound." ACM Transactions on Graphics (TOG) 33.6 (2014): 181.
[0015] This produced higher intensity focal points with a lower computational cost than Ebbini et al.'s weighted generalized inverse algorithm.
[0016] Kappus et al. presented a more efficient implementation of this method by exploiting a property of the pseudo-inverse. Kappus, Brian, and Ben Long. "Spatiotemporal modulation for mid-air haptic feedback from an ultrasonic phased array." The Journal of the Acoustical Society of America 143.3 (2018): 1836-1836.
[0017] Marzo et al. provided the iterative backpropagation algorithm which is essentially an iterative extension of Ibinni et al.'s conjugate field method. Marzo Perez, Asier, and Bruce W. Drinkwater. "Holographic acoustic tweezers." Proceedings of The National Academy of Sciences, 116 (1), 84-89(2019). More recently, Inoue et al. formulated the problem as a semi-definite programming problem and solved it using a block coordinate descent method. Inoue, Seki, Yasutoshi Makino, and Hiroyuki Shinoda. "Active touch perception produced by airborne ultrasonic haptic hologram." 2015 IEEE World Haptics Conference (WHC). IEEE, 2015.
[0018] Further, a phased array of acoustic transducers requires individually addressable elements in order to have flexibility in focus locations. This requires many separate drive circuits which themselves need to be supplied with variable phase and amplitude. This adds cost and sophistication. In many applications, this level of flexibility is not necessary. For instance, one application of phased acoustic arrays is mid-air haptics from focused ultrasound. Adding haptics to a fixed interface may only require 3-4 selectable focus locations and having an array capable of thousands of different focus locations will leave that functionality wasted. This invention shows how to design a phased array system with minimal drive signals to achieve only the necessary number of focal points using specific, derived, arrangements of transducers.
[0019] There have been 2 prior different arrangements which provide a focus using a specific layout.
[0020] 1. A bowl arrangement provides a focus at the geometric center of the bowl with all transducers driven in parallel.
[0021] 2. A ring arrangement where each ring is spaced at one or one-half wavelength from a focus above the center of the rings.
[0022] Further, in order to generate mid-air haptics from a phased array of ultrasonic transducers, a focus must be formed and intersect the hand or body part for which haptics is desired. Generating this field requires a set of driving (phase and amplitude) values to be produced for each transducer in the phased array. In addition, this solution must be updated dynamically as the target moves relative to the array. This is a computationally expensive task.
[0023] There are many different `solvers` which exist to approach this problem. In most cases, they use some kind of transducer model which, along with linear algebra, is used to create a set of activation coefficients. This requires a costly level of computation. This invention instead uses basic transformations on a pre-computed solution stored in memory to deliver a similar level of flexibility using less computation.
[0024] Further, a continuous distribution of sound energy, which will be referred to as an "acoustic field", can be used for a range of applications including haptic feedback in mid-air, sound-from-ultrasound systems and producing encoded wave fields for use by tracking and imaging systems.
[0025] By defining one of more control points in space, the acoustic field can be controlled. Each point can be assigned a value equating to a desired amplitude at the control point. A physical set of transducers can then be controlled to create an acoustic field exhibiting the desired amplitude at the control points.
[0026] By using an ultrasonic carrier frequency, the points of control can be made small enough to be useful and with no consequent audible sound effect. However, a static acoustic field at a power level suitable for a wide range of use cases imparts too small a force to be detectable by human skin, so consequently a fluctuation in time must be added that sensitizes the skin to the acoustic field by exploiting frequencies and skin movement speeds that create more noticeable haptic sensations. Creating extra fluctuations in time using traditional methods involves complicated systems to control the system through time.
[0027] Typically, the amplitude or position must be modified in time to generate a frequency that excites the mechanoreceptors in human skin to produce a strong haptic response. There are a number of ways to provoke such a suitable frequency component, although mostly this is achieved by directly modulating either the position or amplitude at an appropriate frequency.
[0028] This can be achieved indirectly however by using beat frequencies--frequencies generated as the sum or difference between two time-harmonic acoustic fields at different frequencies, such that given an acoustic field at a single frequency .omega..sub.1 and a second acoustic field at a frequency .omega..sub.2, when overlaid, an amplitude modulation effect is created which may be a frequency suitable for creating a haptic effect. This may be written as:
s.sub.1(x,t)=(.phi..sub.1(x)exp(i(.omega..sub.1t+.PHI..sub.1))),
s.sub.2(x,t)=(.phi..sub.2(x)exp(i(.omega..sub.2t+.PHI..sub.2))),
where .phi. is a complex valued function that varies spatially defining the wavefunction, which includes phase and amplitude, with the frequency .omega. and phase offset .PHI..
[0029] By then defining the combined system as the superposition, that is:
s .function. ( x , t ) = ( k = 1 2 .phi. k ( x ) .times. exp .function. ( i .function. ( .omega. k .times. t + .PHI. k ) ) ) = k = 1 2 "\[LeftBracketingBar]" .phi. k ( x ) "\[RightBracketingBar]" .times. cos .function. ( arg .function. ( .phi. k ( x ) ) + .omega. k .times. t + .PHI. k ) , ##EQU00001##
where the equal-amplitude component can be written as a product, which is equivalent to an amplitude modulation of the average by the difference:
= 2 .times. "\[LeftBracketingBar]" .phi. 1 ( x ) "\[RightBracketingBar]" .times. cos .function. ( .rho. 1 + .rho. 2 2 ) .times. cos .times. ( .rho. 1 - .rho. 2 2 ) + ( "\[LeftBracketingBar]" .phi. 2 ( x ) "\[RightBracketingBar]" - "\[LeftBracketingBar]" .phi. 1 ( x ) "\[RightBracketingBar]" ) .times. cos .function. ( .rho. 2 ) ##EQU00002## = 2 .times. "\[LeftBracketingBar]" .phi. 2 ( x ) "\[RightBracketingBar]" .times. cos .function. ( .rho. 2 + .rho. 1 2 ) .times. cos .function. ( .rho. 2 - .rho. 1 2 ) + ( "\[LeftBracketingBar]" .phi. 1 ( x ) "\[RightBracketingBar]" - "\[LeftBracketingBar]" .phi. 2 ( x ) "\[RightBracketingBar]" ) .times. cos .function. ( .rho. 1 ) , ##EQU00002.2##
where the sinusoidal angle is defined by:
.rho..sub.k=arg(.phi..sub.k(x))+.omega..sub.kt+.PHI..sub.k.
[0030] When using this to create an amplitude modulation for acoustic mid-air haptics for example, one could use a 39950 Hz and 40050 Hz tone with equal amplitude to create a 50 Hz amplitude modulation of a carrier frequency at 40000 Hz. This is especially useful for transducers that exhibit strong resonant properties or are otherwise designed for use at a specific frequency; by driving them at slightly offset frequencies to generate the difference tone with high efficiency, this can create specific haptic effects with electronics that only needs to generate simple tones.
[0031] This has previously been exploited to generate single or multiple points of amplitude modulated haptic feedback.
SUMMARY
[0032] Two methods for rendering haptic surfaces through spatial modulation of ultrasound are shown: Singular Value Decomposition (SVD) with Tikhonov Regularization and Mini-Batch, Stochastic Gradient Descent with Momentum (MSGDM). These methods are variations of extant algorithms from the medical ultrasound and human-computer interaction (HCI) literature.
[0033] Further, traditional phased arrays of acoustic transducers require independent drive circuitry for each element to adjust the output field in a useful way. This invention introduces a method whereby allowing for adjustment in the placement of transducers, arrays can be generated which allow for variable focus points using a limited set of driving signals. This is accomplished by generating a set of mutually-exclusive transducer placement layouts which when driven together produce a set of sufficiently orthogonal fields at various points of interest. When driven at specific relative phases the resulting fields sum to produce a specific high pressure acoustic focus at one point of interest while minimizing all others. Focus points can be activated one at a time or in sequence using amplitude modulation or multiple frequencies.
[0034] Further, by sampling a set of activation coefficients on a pre-computed focusing solution, one can reproduce that field. Then, using basic transformations on the sampling locations, this field can be translated or rotated in 3-dimensional space. This enables the creation of spatiotemporal modulated haptics with reduced computation.
[0035] Further, mid-air haptics using ultrasound can be generated via amplitude modulation or via spatial modulation of control points defined at spatial locations in an acoustic field above a phased array device. The amplitude modulation effect can be generated using a previously disclosed method to generate a difference tone at one or more points, however, there has been no prior disclosure of a method to achieve spatial modulation of points in any equivalent fashion. In this disclosure, methods to use the same effect in conjunction with phased arrays at multiple frequencies to generate complex patterns that mimic spatially modulated motion using the phase of the difference tone will be described.
BRIEF DESCRIPTION OF THE DRAWINGS
[0036] The accompanying figures, where like reference numerals refer to identical or functionally similar elements throughout the separate views, together with the detailed description below, are incorporated in and form part of the specification, serve to further illustrate embodiments of concepts that include the claimed invention and explain various principles and advantages of those embodiments.
[0037] FIGS. 1, 2, and 3 are plots depicting the activation coefficients required to generate the Ultrahaptics with the Chebyshev directivity model.
[0038] FIG. 4 is a plot where amplitude utilization is defined as the sum of the amplitudes of all the transducers divided by the total number of transducers.
[0039] FIGS. 5 and 6 are plots demonstrating the measured acoustic field obtained by propagating the activation coefficients with the Chebyshev directivity model.
[0040] FIG. 7 depicts a solution field using reciprocity with array simulation.
[0041] FIG. 8 depicts a solution field using a pseudoinverse with array simulation.
[0042] FIG. 9 depicts a real projection of a set of simulated pressure fields.
[0043] FIG. 10 depicts an imaginary projection of a set of simulated pressure fields.
[0044] FIG. 11 depicts a structured acoustic field generation using a solution field to guide transducer placement.
[0045] FIG. 12 shows a decision flow chart for populating transducers for multiple solutions fields simultaneously.
[0046] FIG. 13 depicts an example execution of the algorithm in FIG. 12.
[0047] FIG. 14 depicts the final populated transducer array built using solution fields from FIG. 13.
[0048] FIG. 15 shows plots of the normalized real projection of the simulated acoustic field produced by each group from FIG. 14.
[0049] FIG. 16 shows plots of normalized pressure field simulation produced by transducer arrangement and grouping from FIG. 14.
[0050] FIG. 17 shows plots of example populated transducer arrays using pre-rotation of source solution fields.
[0051] FIG. 18 depicts an example layout with a center-excluded region.
[0052] FIG. 19 shows a time-domain simulation of a 6-point, 117 transducer at 40 kHz embodiment.
[0053] FIGS. 20, 21, 22, and 23 show a time-domain simulation of a 6-point, 117 transducer at 40 kHz embodiment.
[0054] FIG. 24 is an example 256-element solution sampling.
[0055] FIGS. 25, 26, 27, 28, and 29 show example 256-element solution samplings with offset.
[0056] FIG. 30 shows an example 256-element solution sampling with rotation offset.
[0057] FIG. 31 shows an example 256-element solution sampling with a z-offset transformation.
[0058] FIG. 32 shows an example 256-element solution sampling with a translation in x followed by a z-offset transformation.
[0059] FIG. 33 shows an example 256-element solution sampling, illustrating the sum of two solutions.
[0060] FIGS. 34, 35, 36, 37, and 38 show slices of the acoustic field generated from one of the carrier frequencies, where shading is linearly proportional to amplitude and phase angle.
[0061] FIG. 39 shows a grouped-transducer layout design based upon solution fields generated to form a heart-shaped pattern in a plane at z=10 cm.
[0062] FIG. 40 shows the populated transducer array produced by data from FIG. 39.
[0063] FIG. 41 shows a simulation of the normalized pressure amplitude at z=10 cm generated by data from FIG. 40.
[0064] FIG. 42 shows a normalized pressure produced by an acoustic field simulation of the array shown in FIG. 40.
[0065] Skilled artisans will appreciate that elements in the figures are illustrated for simplicity and clarity and have not necessarily been drawn to scale. For example, the dimensions of some of the elements in the figures may be exaggerated relative to other elements to help to improve understanding of embodiments of the present invention.
[0066] The apparatus and method components have been represented where appropriate by conventional symbols in the drawings, showing only those specific details that are pertinent to understanding the embodiments of the present invention so as not to obscure the disclosure with details that will be readily apparent to those of ordinary skill in the art having the benefit of the description herein.
DETAILED DESCRIPTION
1) Algorithms for Rendering Haptic Surfaces
I. Defining the Problem
[0067] We can approximate the far-field acoustic pressure field emitted by a single ultrasonic transducer as:
.PSI..sub.q(r.sub.kq,X.sub.q)=H(r.sub.kQ)X.sub.q
where [0068] .PSI..sub.q(r.sub.kq,X.sub.q) is the complex-valued, far-field acoustic pressure field emitted by a single ultrasonic transducer. [0069] r.sub.kq is the vector oriented from the transducer q to the sampling point k. [0070] X.sub.q=A.sub.qe.sup.i.phi..sup.q is the complex-valued activation coefficient. [0071] H(r.sub.kq) is the directivity function which models the transducer.
[0072] We then assume that the total field emitted by a phased array of an ultrasonic transducer, excited by independent signals, can be modelled as a linear combination of the fields emitted by each individual transducer yielding
.PSI. .function. ( r k * , x ) = q = 1 n .PSI. q ( r kq , X q ) , ##EQU00003##
where: [0073] r.sub.k*=r.sub.k1, r.sub.k2, . . . , r.sub.kn is a list of vectors from each transducer q to a given point k. [0074] r.sub.**=r.sub.11, r.sub.12, . . . , r.sub.1n, . . . , r.sub.m1, r.sub.m2, . . . , r.sub.mn is a list of all the vectors from each transducer q to each point k. [0075] x=X.sub.1, X.sub.2, . . . , X.sub.n is a vector containing the activation coefficients for each transducer.
[0076] Let .PSI..sup.D be the discretised, desired field defined on a set of points r.sub.1, r.sub.2, . . . , r.sub.m. We can now formalize the problem as an optimization problem of the form:
min x [ k = 1 m .PSI. .function. ( r k * , x ) - .PSI. D ( r k ) ] ##EQU00004##
such that
0.ltoreq.|X.sub.q|.ltoreq.1,.A-inverted.q.di-elect cons.(1,2, . . . ,n)
II. Singular Value Decomposition (SVD) with Tikhonov Regularization
[0077] A. Formulating the Linear System
[0078] Building on the notation above, we can formulate the minimization problem as finding a solution to a complex-valued linear system Ax=b where A is the complex-valued, forward propagation operator, x is the complex-valued vector of activation coefficients and b is the complex-valued, desired field collapsed into a vector. In particular we have,
A .di-elect cons. m .times. n , x .di-elect cons. n .times. 1 , b .di-elect cons. m .times. 1 ##EQU00005## A = [ H .function. ( r 11 ) H .function. ( r 1 .times. n ) H .function. ( r m .times. 1 ) H .function. ( r mn ) ] ##EQU00005.2## x = [ X 1 X n ] ##EQU00005.3## b = [ .PSI. D ( r 1 ) .PSI. D ( r m ) ] ##EQU00005.4##
[0079] B. Column Weighing
[0080] In order to improve the numerical accuracy of any solution to the linear system, we can normalize the operator A using the 2-norm of each column. We can construct a weighting operator, W.sub.c.di-elect cons..sup.n.times.n, defined as:
W c = [ 1 A * 1 2 0 0 1 A * n 2 ] ##EQU00006##
[0081] We can then define a new linear system to solve:
AW.sub.cx=b
[0082] C. Row Weighing
[0083] Standard techniques used to find the least-squares solution to the linear system assume homoscedasticity in the error terms. For this linear system, this means that the error between the generated field and the desired field would be assumed to have the same variance for each point in the desired field:
E[ .sub.k.sup.2|x]=.sigma..sup.2
where,
.sub.k.sup.2=|.PSI.(r.sub.k*,x)-.PSI..sup.D(r.sub.k)|
[0084] When this assumption does not hold, one can improve the results of the solution by weighting the operator A with a weighting operator, W.sub.r.di-elect cons..sup.m.times.m, defined as:
W r = [ W 11 0 0 W m .times. m ] ##EQU00007##
[0085] The diagonal elements, W.sub.kk, reflect the influence that each point in the desired field has in generating the activation coefficients. We can then define a new linear system:
W.sub.rAx=b
[0086] This weighting operator allows a user more control over the field they would like to produce. For instance, they may care more or less about the value of the field in the area outside the desired shape and therefore choose to weight these field points accordingly.
[0087] D. Tikhonov Regularization
[0088] For what follows, let us assume that the user has decided to use neither row weighting nor column weighting. It is worth noting that this does not change the derivation significantly. The complex-valued forward propagation operator can also be interpreted as a mapping from one space to another A:.sup.n.fwdarw..sup.n. Given this mapping, regularisation can be defined as a parametric family of approximate inverse operators R.sub..alpha.:.sup.n.fwdarw..sup.n.
[0089] Tikhonov's regularization seeks to find a solution to the least-squares problem with a small 2-norm. The problem can be expressed as a minimization problem of the form:
min(.parallel.Ax-b.parallel..sub.2+.alpha..parallel.x.parallel..sub.2)
[0090] The solution can be written as:
x.sub..alpha.=R.sub..alpha.b
=(A.sup.HA+.alpha.I).sup.-1A.sup.Hb
[0091] E. SVD Solution
[0092] In general, the singular value decomposition or SVD can be used to solve a linear system. The SVD of a complex-valued operator A is given by A=U.SIGMA.V.sup.H. The solution of a linear system Ax=b can then be computed using the Moore-Penrose pseudoinverse, A.sup.+, as
x=A.sup.+b
where the Moore-Penrose pseudoinverse can be expressed in terms of the SVD
A.sup.+=V.SIGMA..sup.+U.sup.H.
[0093] We can then exploit an identity of the pseudoinverse to reduce the computational cost of the algorithm. Given the relationship:
A.sup.+=(A.sup.HA).sup.+A.sup.H
we can express the solution of the linear system Ax=b as
x=(A.sup.HA).sup.+A.sup.Hb
=V.sub.p.SIGMA..sub.p.sup.+U.sub.p.sup.HA.sup.Hb
where:
(A.sup.HA).sup.+=V.sub.p.SIGMA..sub.p.sup.+U.sub.p.sup.H
[0094] Given that A is an (m.times.n) matrix and A.sup.H is an (n.times.m) matrix the product (A.sup.HA).sup.+ is an (n.times.n) matrix. Since the computation of the SVD is the bottleneck in this algorithm, performing the SVD on an (n.times.n) provides a significant speed up over computing the SVD on an (m.times.n) matrix. This also makes the algorithm more scalable to larger problems, as the computation of the SVD now scales with O(n.sup.3(k+k')) as opposed to O(knm.sup.2+k'n.sup.3).
[0095] F. SVD with Tikhonov Regularization
[0096] The solution to the Tikhonov regularization problem posed in II.D can be expressed using the SVD:
x .alpha. = R a .times. b = ( A H .times. A + .alpha. .times. I ) - 1 .times. A H .times. b = V .function. ( .SIGMA. + .times. .SIGMA. + .alpha. .times. I ) - 1 .times. .SIGMA. + .times. U H .times. b ##EQU00008##
This solution can be rewritten using the identity described in II.E:
x .alpha. = V .function. ( .SIGMA. + .times. .SIGMA. + .alpha. .times. I ) - 1 .times. .SIGMA. + .times. U h .times. b = V p ( .SIGMA. p + .times. .SIGMA. p + .alpha. .times. I ) - 1 .times. .SIGMA. p + .times. U p H .times. A H .times. b ##EQU00009##
[0097] The parameter .alpha. reduces the contribution of the singular values of A to the solution vector x.sub..alpha. as a function of the size of the singular value. The reciprocal of the i-th singular value,
1 .sigma. i , ##EQU00010##
is scaled by the function:
.sigma. i 2 .sigma. i 2 + .alpha. ##EQU00011##
[0098] This means that the reciprocal of smaller singular values is attenuated more than the reciprocal of larger singular values. An (m.times.n) matrix A can be thought of as a linear transformation operator which maps a vector in .sup.n to a vector in .sup.m. The SVD shows us that the linear transformation operator A can be decomposed into the weighted sum of k linear transformations, where k is the rank of A.
[0099] The singular values determine the influence of each of these linear transformations. Large singular values indicate that the contribution of the corresponding linear transformation is large and vice versa. In addition, the magnitude of the singular values determines the magnitude of the norm of the solution vector x.sub..alpha.. Scaling the reciprocal of the singular values with the function defined above, reduces the magnitude of the norm of the solution vector by attenuating the contribution of the less influential linear transformations more than the more influential linear transformations. The idea is that we can reduce the power output of the transducers without significantly affecting the quality of the solution.
[0100] We have implemented this algorithm in Python using numpy. It runs quite fast but could be accelerated using the randomized SGD algorithm and porting the code to OpenCL so that it can run on a GPU.
[0101] G. Optimal Value for the Tikhonov Regularization Parameter .alpha.
[0102] The choice of the Tikhonov regularization parameter .alpha. affects the quality and magnitude of the solution vector obtained x.sub..alpha.. In particular, there is a trade-off between the accuracy and magnitude of the solution. For large values of .alpha. the solution x.sub..alpha. is inaccurate but has a small magnitude and vice versa. There are various methods that can be used to determine the optimal value of .alpha. for the given problem such as; the discrepancy principle, the L-curve and the GCV. These methods are generally computationally expensive as they involve solving a minimization problem.
[0103] To obviate this computational cost, we suggest determining reasonable values of .alpha. with a data-based approach. An example of such an approach could involve determining the best .alpha. for a range of different surfaces and then finding the average .alpha..
III. Mini-Batch, Stochastic Gradient Descent with Momentum (MSGDM)
[0104] We can use the mini-batch, stochastic gradient descent (SGD) with momentum algorithm to solve a minimization problem similar to that posed in Section I. We begin by defining an error function at a given point k as:
E(r.sub.k*,x)=(.PSI.(r.sub.k*,x)-.PSI..sup.D(r.sub.k))(.PSI.(r.sub.k*,x)- -.PSI..sup.D(r.sub.k)).
[0105] Given this error function, we can define a real-valued objective function that we would like to minimize:
F .function. ( r B * , x ) = 1 b .times. k .di-elect cons. B E .function. ( r k * , x ) , ##EQU00012##
where [0106] B is a set of b uniformly selected points from the range {1, 2, . . . , m} where b
All transformations within the plane (translation or rotation) must be done before applying the z-transformation and the post-translation or rotation then fed into the z-offset transformation.
[0266] FIG. 32 illustrates an example of a translation followed by a z-offset transformation 3200. The z-offset transformation requires more computation than translation and could be beyond the capability of some hardware implementations to execute quickly. One operation that is possible after a z-offset transformation which is computationally inexpensive is to offset the solution from one transducer to an adjacent transducer. This shifts the resulting field by the pitch of the transducers in the offset-z plane. This could be done with a shift of one pitch unit or many or even including rotation or mirroring. Transducers near the edge or in some way do not have a solution to shift to can either be disabled or a new solution be generated for them via a new lookup-location calculation or direct computation.
[0267] Shown in FIG. 32 is an example 256-element solution sampling with a translation in x followed by a z-offset transformation 3200. Pre-generated solutions are given in graph 3210 having scale intensity 3215 A (phase in radians) and graph B 3220 having scale intensity 3225 (normalized amplitude). These are sampled at 1 mm pitch and generated to produce a focus at [x,y,z]=[0,0, 10 cm] relative to the center of the array at 40 kHz. The 256-element simulated array has 1.03 cm pitch and the sampled values are given by `x`s in Graphs A 3210 and B 3220. In this case the sampling is first translated in x by -4 cm and then adjusted by x.sub.new.sup.2+y.sub.new.sup.2=z.sub.new.sup.2-z.sub.old.sup.2+x.sup.2+y- .sup.2, with z.sub.new given by 15 cm. Graph C 3230 having intensity scale 3235 in radians shows the resulting sampling with phase given by shade and amplitude given by the radius of the inner, filled, circle. In this case the amplitude is rescaled so that at least one transducer is at full drive and the rest are proportional to that value as described by the solution in Graph B 3220. Graph D 3240 having intensity scale 3245 in normalized pressure shows an x-z cross-section of the resulting simulated field pressure normalized to maximum=1.0. As expected, the focus appears at [4 cm,0,15 cm].
[0268] While transducer locations and stored solutions have largely been discussed in the context of rectilinear coordinates, this is by no means a requirement. Any coordinate system can be used as long as it maps to a 2-dimensional space. In this way, certain solutions can exploit symmetries to reduce computation, memory usage, or both. A solution to a single focus point, for instance, is radially symmetric, and can be stored as a function only of radius from the origin. After transforming the location of each transducer (through any of the above presented transformations), the radius from the origin then determines the driving phase and amplitude. While this reduces solution memory storage, it might not be the fastest implementation in hardware. In another embodiment, the solution is kept in rectilinear coordinates but only one quadrant is stored. Since the solution is symmetric about both the x and y axis, it is relatively easy to map a negative-x or negative-y coordinate to the positive equivalent, thereby reducing memory usage.
[0269] Certain transducer layouts can be used to simplify the transformations. Rectilinear arrangements make translation in x or y simple as the offset for one transducer applies to all, but changes in z take more computation.
[0270] Arrangements with many transducers at equal radius from the origin can simplify changes in z when x and y are zero, as one radius transformation will apply to multiple transducers.
[0271] The sampled solution amplitude represents only one possible amplitude in the solution. If, for example, none of the transducers are being driven at maximum, the solution can be scaled up to the point where at least one transducer is driven at full-scale. This scaling factor can be found by searching the array of sampled values and finding the max drive from the solution. To preserve the shape of the field, the remaining transducers can be scaled by the same factor.
[0272] Another application of scaling is to modulate the field. After searching the solution array for the maximum value, this determines then the range of the possible amplitudes. If less total amplitude is desired, for instance if amplitude modulation of the field is desired, all activation coefficients can be scaled together by a time-varying function. This can happen on a fixed sample arrangement or dynamically as the sampled locations are adjusted.
[0273] In another embodiment, the solution amplitude is ignored and only the phase sampled in the solution while leaving the amplitude at a fixed value (full-drive, for instance). For a single focus point solution, this will preserve any focus location, but may change the field away from the foci depending on the transformations applied to the sampled locations. It does, however, greatly simplify amplitude modulation as all transducers will receive the same amplitude.
[0274] Two (or more) sets of sampled activation coefficients can be summed (using complex phasor notation) to produce all output fields simultaneously. For this to function, first, the summed coefficients together need to be under the maximum drive for every transducer. If they exceed the max, the solution must be scaled (by the maximum value, or more) or clipped (all values over the maximum drive, are set to maximum while preserving phase). Second, it is not a guarantee that the field produced by the first solution will not interfere and change the other solutions in some meaningful way. Using the example of a single focus solution, sampled to produce two spatially distinct foci, depending on the spatial relationship of each focus to each other and the physical layout of the array, the fringing field from one focus can reduce the second or vice versa. In other arrangements of foci location, one could constructively interfere and increase the pressure of one or both foci beyond what was intended. In one embodiment, this is ignored-interference can be a rare occurrence in many arrangements.
[0275] FIG. 33 illustrates an example simulation of two solutions being produced simultaneously without any phase rotation 3300. In another embodiment, the array output could be simulated and all possible combinations of points which cause interference could be predicted and avoided. In another embodiment, in cases of interference, the phase of one or more solutions are rotated to avoid or mitigate the interference. The locations which cause interference and what rotation values to use can be stored and reference during operation. In another embodiment, a machine learning model is trained using acoustic simulations or real data to predict locations of possible interference and return phases and amplitudes needed to remedy the interference for each solution.
[0276] FIG. 33 shows an example 256-element solution sampling, illustrating the sum of two solutions 3300. Pre-generated solutions are given in graph A 3310 having intensity scale 3315 (phase in radians) and graph B 3320 having intensity scale 3325 (normalized amplitude). These are sampled at 1 mm pitch and generated to produce a focus at [x,y,z]=[0,0, 10 cm] relative to the center of the array at 40 kHz. The 256-element simulated array has 1.03 cm pitch and one of the two sampled solutions are given by `x`s in graphs A 3310 and B 3320. In this case the sampling is first translated in x by +2 cm and then adjusted by x.sub.new.sup.2+y.sub.new.sup.2=z.sub.new.sup.2-z.sub.old.sup.2+x.sup.2+y- .sup.2, with z.sub.new given by 13 cm. The other solution (not shown) is rectilinear sampling translated by -4 cm without a z-offset transformation (same as FIGS. 25 and 26). These two sampled solutions are summed (in complex phasor notation) with normalization coefficients 0.65 and 0.35 respectively. Graph C 3330 having intensity scale 3335 in radians shows the resulting solution with phase given by shade and amplitude given by the radius of the inner, filled, circle. In this case the amplitude is rescaled so that at least one transducer is at full drive and the rest are proportional to that value. Graph D 3340 having intensity scale 3345 in normalized pressure shows an x-z cross-section of the resulting simulated field pressure normalized to maximum=1.0. As expected, two foci appear at [-2 cm,0,13 cm] and [4 cm,0,10 cm].
[0277] Other methods also exist to manage possible interference between summed fields. The key to mitigating the interference is to predict the output of each field at key points of interest. With all fields estimated at all points of interest, the results can be summed to predict final field. Points of interest can include focus points, points of known objects nearby the array, a reflective surface, or similar. Any one field can be driven by a different phase than is stored in the solution by rotating the phase of all activation coefficients by the same amount. This affects every point in the produced field by the same amount and the sum at each point of interest will be similarly affected. This divides the problem into two separate problems: 1. How to estimate the field at arbitrary points of interest and 2. How to choose phases of individual field solutions so that interference is minimized.
[0278] Estimating the field has a number of solutions. In one embodiment, a mathematical model of each transducer's acoustic field is used to estimate the contribution of each transducing element to each point of interest. The sum of all these contributions multiplied by their individual activation coefficients (all in complex phasor notation) is then the estimate of the field at the point of interest. In another embodiment of the invention, a complete simulation of the field when it is focused at each location in the interactive volume is stored in memory and then referenced as at each point of interest for each focus location. Interpolation may need to be used if the stored solution is at low resolution and higher resolution is needed. In another embodiment, a mathematical model of the field around a focal point is used rather than a model from individual transducer. This model includes the steering direction and distance from the array and can be validated against measured data. Yet another embodiment uses a machine-learning model which performs the same calculation and is trained with simulation and/or experimental data.
[0279] After the field is estimated at each point of interest a phase-oracle algorithm can be used to rotate phases and scale amplitudes of each solution to optimize the output. Power iteration on the phase of each solution's activation and resulting phase of each point of interest is one method. Other methods exist in the literature.
[0280] Transitioning from one set of sampled coefficients to another requires care to avoid sharp changes in the acoustic field which can lead to audible artifacts. Note that this transition includes both transformations within one solution as well as transitions to a completely different stored solution which can include a different transformation. In one embodiment, the system limits the magnitude of the change of each coefficient in the complex phasor domain. If the desired change is larger than this magnitude, the system steps the activation coefficient by the maximum amount allowed until the desired value is reached. If the desired change is modified before reaching the previously desired value, the system starts stepping in the new direction from its last state. In another embodiment, this step is dynamically sized so that the interpolated value contains a minimum bandwidth of modulation. An example of this is a sinusoidal profile with the start and finish of the transition approaching a zero time-derivative. In another embodiment, a proportional-integral-derivative (PID)-type controller framework is used to manage the transitions for each coefficient. In another embodiment, the real and imaginary part of the activation coefficient are fed to a low-pass frequency filter of arbitrary architecture. The output of this filter is used as the driving coefficient. In another embodiment, the stored solution itself is used as a moderator of transition speed by limiting the next-sampled solution value to be within a certain maximum distance within the solution. If the desired solution is more than one distance-step away from the current, then it chooses the value most towards the desired value within the limits of the distance-limit. For example, if the solution is in cartesian coordinates, one could limit the change to be an adjacent cell within the 2-dimensional solution and the direction to step would proceed along the direction in the difference vector with the highest magnitude. This direction is updated with each change in desired solution. In yet another embodiment, the amplitude of the drive is reduced to zero (via a linear or variable-step ramp) and then ramped back up at the new phases, thus avoiding a sharp transition. This method is the most extreme as the ramp would take time and during this period the fields would be lower in magnitude. This technique could be used in concert with others and only utilized when some criteria is entered such as the average degree of change required.
[0281] To change to a new field, in most cases every transducer requires a new driving condition. When all transducers are changed simultaneously (such as in the same acoustic cycle), the changing fields act in unison to produce audible artifacts. To mitigate this effect, the system can limit the number of transducers which receive changes per acoustic cycle. In one embodiment, only one transducer is changed per acoustic cycle. This transducer could be chosen to maximize the change in the field (such as magnitude change), at random, or some other metric. In another embodiment, the system would limit the number of transducers allowed to change to some value n where n is less than the total number of transducers in the system. The selection of transducers could be made at random or with a quality metric. An example quality metric would be magnitude of the complex driving difference sum. Another could be the n largest magnitude changes. The change in sampling could be limited so that every transducer has a chance to transition before choosing a new set of samplings, but this is not required. If requested changes build up faster than the coefficients are allowed to change, however, then the field will no longer be manifesting as desired, and could be difficult to predict, but will tend preferentially to uniform noise.
[0282] The solution density used is a balance between memory and field accuracy. The more dense the solution field stored, the finer the system will be able to manipulate the field. This includes all of the transforms given above. In particular, when a solution is modified with a z-offset transformation, this necessarily groups the sampled values closer together. As long as the sampled separation is much larger than the pitch of the solution, the field produced will be close to ideal. If the sampling starts to give the same value to multiple transducers, then the system will tend towards uniform beaming behavior. In one embodiment, a higher-resolution can be simulated using interpolation of adjacent solution values. In another embodiment, losses compression is used to reduce the memory usage of a stored solution which is then decompressed in real-time or into memory before use.
[0283] Using a set of basis functions with an independent time variable and solving a least squares system to achieve lossy compression of the transducer drive conditions may be considered as another method to achieve a further trade-off, in this case memory and decoding complexity. A set of basis functions may be sampled at a number of points in time to form an A matrix, the intended drive conditions for a given transducer form a b vector and the compressed set of coefficients x that may be reconstructed into a set of functions whose summations approximate the original drive conditions when the linear system Ax=b is solved off-line. This linear system may be written as:
A .times. x = b = [ f 1 ( t 1 ) f m .times. ( t 1 ) f 1 .times. ( t N ) f m .times. ( t N ) ] [ x 1 x m ] = [ .alpha. .function. ( t 1 ) .alpha. .times. ( t N ) ] ##EQU00029##
where there are m basis functions f.sub.1(t), . . . , f.sub.m(t), N samples of the transducer drive conditions where at least part of the transducer drive condition under consideration is .alpha.(t).
[0284] As harmonic fields that repeat in time (but more slowly than the carrier frequency) are often the most useful, these can be represented as basis functions with periodic boundary conditions and can have a simple trigonometric or wavelet decomposition that allow for example series of sines, sawtooth waves or square waves with different frequencies to reconstruct a set of transducer drive conditions that may be then applied to the array by evaluating the summation over time. While this kind of implementation would generally require multiplications due to the presence of different coefficients, other methods may be used to take advantage of specific choices of basis function such as CORDIC for trigonometric functions, additive or subtractive counters for sawtooth waves and bit-wise manipulation for square wave generation among others.
[0285] Normalizing the sampled solution so that at least one transducer (or all transducers) have maximum drive generates the maximum pressure in the field possible from the array. If, however, a lower pressure is desired, the system can normalize the array activation coefficients with a lower value. In some instances, such as controlling the pressure produced at a focal point solution, a specific pressure is desired. In those cases, the pressure at one or many points of interest from the full-power solution must be estimated, this then can be used as a scaling factor to proportionally scale the solution's amplitudes down to the desired value. Various methods of field estimation have been discussed above and all could be used in this context. Other solutions exist, however, if points of interest follow reduced degrees of freedom.
[0286] One example of reduced degrees-of-freedom is if the system only needs to control the pressure at a single focal point solution. In one embodiment, the resulting focus pressure for every translation and z-offset transformation is simulated or measured in advance and stored on the device. Sampling of this solution provides a scaling value to use to scale the sampled solution to the desired pressure. Resolution of the focus solution pressure need not be higher than lambda/10 where lambda is the wavelength of ultrasound used. If memory is an issue, a lower resolution can be used and sampling could be rounded or interpolated. In another embodiment, the only pressure stored in memory is the default solution pressure, without translation or transformation. When the solution is resampled, the offset from the default distance from the origin of the array is calculated and the default pressure is scaled by the reciprocal of the new distance. For example, if the original solution has a pressure P.sub.0 at a distance of r.sub.0, a reasonable approximation of the pressure P at a new distance r is given by
P = P 0 .times. r 0 r . ##EQU00030##
This approach works for single focus solutions or more complicated field solutions. In another embodiment, an alternate equation is used as a pressure estimate which better matches simulation or experiment such as a polynomial or more sophisticated fit to real or simulated data. In another embodiment, a machine-learning model is trained with real or simulated data to predict pressure for a given offset from the original solution and then is used to scale the solution.
[0287] For multipoint solution summation, the pressure estimate for each solution is used to scale each solution before summation to assure appropriate levels. For instance, if one solution provides a pressure of 1 (in scaled units) and the other pressure provides a pressure of 0.5 and the desire is to have each solution at the same pressure, the solution would be to weight the second solution with double the amplitude of the first one. In that way, the resulting output would be similar.
[0288] In summary, this invention provides a method of sampling a stored solutions within an acoustic phased array which manipulates the output field in a controlled way. This reduces computation requirements and increases flexibility of simple acoustic phased array systems.
[0289] II. Additional Disclosure
[0290] 1. An acoustic phased array system comprising,
A plurality of transducers with known positions and orientations
Memory
[0291] A set of driving phases stored in memory A mapping of the stored driving phases to each transducer A second, distinct mapping of the stored driving phases A electronic diving circuit The electronic driving circuit uses the first mapping to power the transducers at the specified phases The electronic driving unit then uses the second mapping to power the transducers at the specified phases 2) A system as in 1, where amplitude is stored alongside phase in memory and is included in the mapping. Subclaims involving various transforms of the mapping Subclaims involving pressure estimation and scaling Subclaims involving multipoint possibilities and how to sum them. Subclaims involving interpolation of the solution Subclaims involving storing multiple solutions in memory. Subclaims involving transition techniques.
4. Spatio-Temporal Modulation Using Multiple Carrier Frequencies
[0292] I. Introduction
[0293] Spatial modulation may be generated by generating paths through space for control points of relatively high acoustic pressure to travel along. Haptics may be created in this way along both open and closed paths. By creating a diffraction pattern from a phased array at a single carrier frequency with the phase changing progressively along the curve, when the same pattern is produced by phased arrays with multiple carrier frequencies, the phase movement along the path can be made to generate the emergent effect of a moving high acoustic pressure region.
[0294] Referring back to the beat frequencies equation:
= 2 .times. "\[LeftBracketingBar]" .phi. 1 ( x ) "\[RightBracketingBar]" .times. cos .function. ( .rho. 1 + .rho. 2 2 ) .times. cos .function. ( .rho. 1 - .rho. 2 2 ) + ( "\[LeftBracketingBar]" .phi. 2 ( x ) "\[RightBracketingBar]" - "\[LeftBracketingBar]" .phi. 1 ( x ) "\[RightBracketingBar]" ) .times. cos .function. ( .rho. 2 ) ##EQU00031##
where 1 and 2 are interchangeable due to the symmetry of cosine around the origin, and:
.rho..sub.k=arg(.phi..sub.k(X))+.omega..sub.kt+.PHI..sub.k.
[0295] Given that in each case, a phased array may be assumed to be creating the spatially defined complex-valued functions .phi..sub.k(x), as these are controlled in phase, .PHI..sub.k may be absorbed into this function and so set to zero. The phase angle of the low frequency component may then be written as:
.rho.'=(arg(.phi..sub.1(x))-arg((.phi..sub.2(x)))+(.omega..sub.1-.omega.- .sub.2)t,
where (.omega..sub.1-.omega..sub.2)t can be described as the progressive part of the phase angle in the low frequency, with arg(.phi..sub.1(x))-arg(.phi..sub.2(x)) describing the phase offset for each spatial location at that low frequency.
[0296] Both carrier frequencies are expected in most embodiments to be high frequencies with a small difference, so components that are only high frequency, such as the component proportional to the difference between the absolute amplitudes (|.phi..sub.2(x)|-|.phi..sub.1(x)|)cos(.rho..sub.2) can be neglected. However, the objective here is to maximize the effect of the low frequency component cos 1/2(.rho..sub.1-.rho..sub.2). The most efficient approach is therefore two-fold, first maximising each |.phi..sub.k(x)| while minimising the difference between them (as difference would result in wasted high frequency signal that would not be modulated), and second ensuring spatially neighbouring locations differ in phase .rho.' to create movement over each low frequency period.
[0297] One way to describe the target field given a desired speed of movement v(x) is by writing the phase angle in the form of an Eikonal equation:
|.gradient..rho.'|=|.gradient.(arg(.phi..sub.1(x))-arg((.phi..sub.2(x)))- |=v(x),
with the implication that this may be achieved in multiple dimensions yielding moving lines or surfaces, where initial conditions can denote times where the haptics is present and may be solved for using a differential equation formulation coupled with an algorithm such as fast marching.
[0298] Equally, for simple scenarios constructing .gradient..rho.' as the direction and magnitude of the apparent low frequency movement is straightforward and does not require a solver algorithm to yield a set of phases that can be used to generate a high acoustic pressure point with apparent motion at the low frequency modulation. This can be achieved by heuristics, such as for example taking a closed curve and a desired number of points and for a constant speed of movement, determining the length of the curve and generating a set of control points that generate the phases required at each point along the curve by dividing the segment of the curve travelled by the total curve length and multiplying through by the number of 2.pi. rotations desired. In general, painting a contiguous set of points with similar amplitude magnitude and moving phase generates a point of high acoustic pressure with apparent movement.
[0299] II. Solving for Acoustic Sources
[0300] Once the sets of points have been set up with amplitude and phase, then to create these points a `phase plate` or some region of space which is intended to act as a source must be populated with virtual acoustic sources which are solved for which will then create the field specified. This has also been referred to as a `solution field`. One method to achieve this is by using an iteratively re-weighted quadratic maximization.
[0301] Taking the set of unit amplitude and zero phase sources, these can be evaluated at any location (x,y,z) yielding a column vector of linear complex-valued acoustic outputs (such as pressure or medium particle velocity) from acoustic sources as .PSI.(x,y,z):
.PSI. = [ .PSI. 1 ( x , y , z ) .PSI. n .times. ( x , y , z ) ] , ##EQU00032##
[0302] Taking a spatially defined complex-valued weighting function w(x,y,z) (initialised to one to begin with at any target point) the basis vector of source contributions for any target point may be reweighted, yielding:
f.sub.q(x.sub.j,y.sub.j,z.sub.j)=v.sub.qw(x.sub.j,y.sub.j,z.sub.j).PSI..- sub.q(x.sub.j,y.sub.j,z.sub.j),
where the overline denotes the operation of complex-conjugate and v.sub.q is a real weighting on the acoustic source function that may be used to control the proportion of each source used, in the case that sources are over- or under-utilised. Then a column vector may be defined as for m control points:
f = [ v 1 j = 1 m f 1 ( x j , y j , z j ) v q j = 1 m f q ( x j , y j , z j ) v n j = 1 m f n .times. ( x j , y j , z j ) ] , ##EQU00033##
where the summation over m may be replaced by a continuous integral if required.
[0303] Then using linearity, this may be multiplied component-wise by complex activation values `driving` each source (x.sub.q) leading to the vector:
m = [ x 1 v 1 j = 1 m f 1 .times. ( x j , y j , z j ) x q v q j = 1 m f q .times. ( x j , y j , z j ) x n v n j = 1 m f n .times. ( x j , y j , z j ) ] , ##EQU00034##
where n is the number of acoustic sources. Driving an array of acoustic transducers at a fixed frequency with such an input x vector to maximise an acoustic quantity encoded in a matrix M may be expressed as a linear algebra problem, so that the intention is to: [0304] maximize x.sup.HMx [0305] subject to x.sup.Hx=1, [0306] and x.di-elect cons..sup.n.
[0307] This may be solved by taking the dominant eigenvector of the matrix M, as the statement of the problem is also equivalently the definition of an eigenvector, here x where:
x = [ x 1 x q x n ] , ##EQU00035##
[0308] This may be solved by taking the dominant eigenvector of the matrix M, as the statement of the problem is also the definition of an eigenvector, here x. Building the matrix M can be written now as:
m.sup.Hm=x.sup.HMx=x.sup.Hff.sup.Hx,
so,
M=ff.sup.H.
[0309] Therefore, the iterative method may be obtained by on each iteration determining the dominant eigenvector, weighting the output vectors to yield on average the correct amplitude level and then reweighting the amplitude between each control point using the weight update equation:
w t + 1 ( x j , y j , z j ) := w t ( x j , y j , z j ) .PSI. r ' ( x j , y j , z j ) .PSI. t ' ( x j , y j , z j ) , ##EQU00036##
where .PSI.'.sub.r(x.sub.j,y.sub.j,z.sub.j) is the desired total complex-valued linear acoustic quantity and .PSI.'.sub.y(x.sub.j,y.sub.j,z.sub.j) is the sum generated by the solution at iteration t as:
.PSI. t ' ( x j , y j , z j ) = q = 1 n x q , t v q , t f q ( x j , y j , z j ) , ##EQU00037##
so although this iteratively re-weighted maximization method cannot easily produce zero control points, it is well suited to generating fields for this multiple frequency method for generating high acoustic output points with apparent movement at the beat frequency.
[0310] If the two frequencies are close enough, then the solution may be reused for both acoustic source systems, although separate solutions may always be obtained for each carrier frequency.
[0311] As shown in FIG. 34, constructing movement for closed curves is subject to some restrictions due to the winding number of the phase only taking integer values.
[0312] Specifically, FIG. 34 shows slices of the acoustic field 3400 generated from one of the carrier frequencies, where shading is linearly proportional to amplitude and phase angle. This creates a haptic point running around a heart-shaped closed curve. Presented here are four different phase winding counts to create different perceived movement speeds on the skin using the same difference frequency. The top left 3410 completes 2.pi., top right 3420 completes 4.pi., bottom left 3430 completes 8.pi. and bottom right 3440 completes 12.pi., where each corresponds to a haptic frequency of the point completing the trajectory of 1, 1/2, 1/4 and 1/6 of the difference frequency respectively. Note the linear shading in phase angle may distort the appreciation of the amplitude of the result.
[0313] In FIG. 35, it is shown that fewer restrictions are present when dealing with open curves and that asymmetry, multiple curve segments and different speeds are possible within the same field.
[0314] Specifically, FIG. 35 shows a slice of the acoustic field 3500 generated from one of the carrier frequencies, where shading is linearly proportional to amplitude and phase angle. Here it is shown that the movements need not be restricted to continuous curves (top left 3510), need not be symmetrical (3520 top right), need not be a round multiple of 2.pi. (bottom left 3530, 16.pi./6) and need not have each segment of discontinuous curve have the same point speed (bottom right 3540).
[0315] In FIG. 36, it is shown that different directions of travel are possible by manipulating the direction of .gradient..rho.'.
[0316] Specifically, FIG. 36 shows slices of the acoustic field 3600 generated from one of the carrier frequencies, where shading is linearly proportional to amplitude and phase angle. Here it shown that there is also freedom to choose the apparent movement direction of the generated points. While the left sensation 3610 is generated radially, the right sensation 3620 has two orthogonal lines that cross together through the center point.
[0317] FIG. 37 shows a variety of different field partitioning schemes that yield the same scenario of a small heart shape in the center of the field with four high acoustic pressure regions, as shown in many of the later figures. FIG. 38 shows that the method is not limited to simple paths, as wavefronts and other features may be constructed that are not simple control points.
[0318] Specifically, FIG. 37 shows slices of the acoustic field 3700 generated from each of the carrier frequencies, where shading is linearly proportional to amplitude and phase angle. The left column 3710 3720 3730 is the first frequency and the right column 3740 3750 3760 is the second frequency, subtracting each phase configuration across the top 3710 3740, center 3720 3750 and bottom 3730 3760 rows generates the same four point system moving around a small heart shape.
[0319] FIG. 38 shows an acoustic field 300 with slices 3810 3820 generated from one of the carrier frequencies, where shading is linearly proportional to amplitude and phase angle. These show that the approach is not limited to control points but can create higher dimensional phase features.
[0320] One way to generate phased arrays using the `phase plates` or solution fields described above is to populate transducers on locations of similar amplitude and phase. This has been described elsewhere, specifically in U.S. Ser. No. 63/156,829, filed on Mar. 4, 2021, which is incorporated by reference in its entirety.
[0321] Starting with a complex-valued solution, first a projection is made to a particular phase. Next, transducers are placed starting with the largest values and proceeding down to lower values in order. When a transducer is placed, an exclusion region based on the size of the transducer and its minimum spacing to other transducers is formed. Solution field values in this exclusion region are not considered for the next placement. When multiple solution fields are considered simultaneously, the exclusion region formed is applied to all fields. Fields are prioritized for population based on the sum of the solution field values for each placed transducer for each respective field. By populating based on a given projection, and driving all the transducers in phase, the desired acoustic field based on the solution field is produced.
[0322] Error! Reference source not found. FIG. 39 shows the construction of a transducer layout populated using a solution field generated to form a heart-shaped high-pressure region at z=10 cm. As shown in FIG. 33, the field which results from this solution produces a shape of high pressure with one or more high-pressure regions moving around the shape at the carrier frequency. Various methods exist to generate solutions which change the rate and direction of these high-pressure regions along the same curve. For instance, adding more phase winding slows the movement of said regions, but also adds more. The example shown in Error! Reference source not found.ure 39 changes the movement by mirroring the solution field along the y-axis. Due to the symmetry of the field, this reverses the helicity of the high-pressure region's movement.
[0323] Specifically, FIG. 39 shows a grouped-transducer layout design 3900 based upon solution fields generated to form a heart-shaped pattern in a plane at z=10 cm. This corresponds to the upper-left simulation 3410 in FIG. 34. The solution fields shown are real projections with light values positive and dark values negative. Transducers for each of group 1 3910 and group 2 3920 are shown as circles. Some circles are truncated due to the extent of the solution field in this particular visualization. The fields are mirrored about they-axis to produce fields of opposite helicity which yields discrete high-pressure points when summed.
[0324] FIG. 40 shows the populated transducer array 4000 produced by solution fields group 1 4012 and group 2 4014 from FIG. 39.
[0325] FIG. 41 shows a simulation of the normalized pressure amplitude at z=10 cm generated by each of group 1 4110 and group 2 4120 as shown in FIG. 40 at 40 kHz. A heart-shape is clearly visible in each plot. Each group is driven independently and in phase to form each field. The amplitude is shown rather than a real projection for clarity.
[0326] Producing these two fields simultaneously results in 4 high-pressure points which lie on the heart-shape as in FIG. 42, upper-left. Adjusting the relative phase of each group, the 4 points process around the heart-shape. This can yield a spatiotemporal haptic by adjusting the relative phase in a steady manner, repeating at a frequency for which human skin is most sensitive such as 100 Hz. One way this can be realized by driving one group at 50 Hz above the carrier frequency and one group 50 Hz below, as although due to the two on the denominator creating a modulation of 50 Hz, the modulation of the envelope by both positive and negative portions of the low-frequency sinusoid creates high acoustic pressure regions, which then have an apparent motion at 100 Hz, since as shown in FIG. 42 they are not distinct. The relatively small difference relative to an ultrasonic carrier frequency such as 40 kHz will preserve the heart-shaped solution while providing a continually, smoothly changing relative phase.
[0327] Specifically, FIG. 42 shows a normalized pressure field 4200 produced by an acoustic field simulation of the array shown in FIG. 40 at 40 kHz at z=10 cm with relative phase between groups 1 and 2 above each plot 4210 4220 4230 4240. The 4 points shown in the upper-left FIG. 4210 smoothly move around the heart-shape clockwise as the phase progresses.
5) Conclusion
[0328] In the foregoing specification, specific embodiments have been described. However, one of ordinary skill in the art appreciates that various modifications and changes can be made without departing from the scope of the invention as set forth in the claims below. Accordingly, the specification and figures are to be regarded in an illustrative rather than a restrictive sense, and all such modifications are intended to be included within the scope of present teachings.
[0329] Moreover, in this document, relational terms such as first and second, top and bottom, and the like may be used solely to distinguish one entity or action from another entity or action without necessarily requiring or implying any actual such relationship or order between such entities or actions. The terms "comprises." "comprising," "has", "having," "includes", "including," "contains", "containing" or any other variation thereof, are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises, has, includes, contains a list of elements does not include only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. An element proceeded by "comprises . . . a", "has . . . a", "includes . . . a", "contains . . . a" does not, without more constraints, preclude the existence of additional identical elements in the process, method, article, or apparatus that comprises, has, includes, contains the element. The terms "a" and "an" are defined as one or more unless explicitly stated otherwise herein. The terms "substantially", "essentially", "approximately", "about" or any other version thereof, are defined as being close to as understood by one of ordinary skill in the art. The term "coupled" as used herein is defined as connected, although not necessarily directly and not necessarily mechanically. A device or structure that is "configured" in a certain way is configured in at least that way but may also be configured in ways that are not listed.
[0330] The Abstract of the Disclosure is provided to allow the reader to quickly ascertain the nature of the technical disclosure. It is submitted with the understanding that it will not be used to interpret or limit the scope or meaning of the claims. In addition, in the foregoing Detailed Description, various features are grouped together in various embodiments for the purpose of streamlining the disclosure. This method of disclosure is not to be interpreted as reflecting an intention that the claimed embodiments require more features than are expressly recited in each claim. Rather, as the following claims reflect, inventive subject matter lies in less than all features of a single disclosed embodiment. Thus, the following claims are hereby incorporated into the Detailed Description, with each claim standing on its own as a separately claimed subject matter.
