WO2024119366A1

WO2024119366A1 - Method for calculating of sound wave diffraction at straight edge based on path tracking

Info

Publication number: WO2024119366A1
Application number: PCT/CN2022/136882
Authority: WO
Inventors: 曹春晓; 任重; 安自立; 周昆
Original assignee: 浙江大学
Priority date: 2022-12-06
Filing date: 2022-12-06
Publication date: 2024-06-13

Abstract

Disclosed is a calculation method for diffraction of a sound wave at a straight edge, based on path tracking. The present method expresses diffraction at a convex straight edge as the integral of a radial function on a geometric straight edge, and the integral is calculated using a path tracking method, thereby achieving efficient calculation of diffraction at a straight edge under a path tracking framework. In addition, the diffraction calculation method of the present invention provides a series of matched optimization technologies, so as to improve the quality of a calculation result, and reduce resource consumption of a calculation result. Experimental results show that the result of the method of the present invention in a simple scene is consistent with the accurate solution for diffraction at a convex straight edge, and in a large complex scene, the method can provide a diffraction simulation result for a user at nearly real-time speeds.

Description

A method for calculating straight-edge diffraction of acoustic waves based on path tracking

Technical Field

The invention relates to the technical field of acoustic simulation, and in particular to a method for calculating straight-edge diffraction of sound waves based on path tracking.

Background technique

Sound propagation simulation technology is widely used in fields such as architectural design and virtual reality. Among the sound propagation simulation algorithms, the method based on geometric acoustics is more common in practical applications due to its high computational efficiency. The geometric acoustics method simplifies the wave equation and assumes that sound waves propagate in a straight line. This simplification requires independent processing of phenomena such as diffraction and interference. The lack of diffraction components means that when the listener enters the "shadow" area blocked by the object, the sound propagation simulation results will show obvious discontinuity, which is inconsistent with the daily acoustic experience and greatly reduces the realism of the simulation results.

It is a complicated matter to put the diffraction phenomenon into the framework of geometric acoustics. In actual engineering, an impulse response representing diffraction is usually superimposed on the impulse response calculated using ordinary geometric acoustics theory. The impulse response representing diffraction here can be calculated using the wave equation. But in more cases, it is hoped that the calculation of diffraction is similar to the ordinary geometric method. A natural idea is to first split the scene into some simpler structures that are easy to calculate diffraction, and then piece together the calculation results to form the final diffraction solution. Generally, triangular meshes are used to represent objects in the scene, and the shadow boundaries caused by object occlusion are related to the edges in the triangular mesh. Therefore, you can try to calculate the diffraction near each edge of the triangular mesh. This is the motivation for geometric acoustics to care about the representation of straight-edge diffraction.

Straight edge diffraction, or wedge diffraction, refers to diffraction near two infinite half-planes connected by a straight edge and forming a certain angle. In the middle of the last century, Keller first proposed a method to describe convex straight edge diffraction in geometric optics, called the geometric diffraction theory. Keller's results are close to the exact solution in most directions, but this model fails near the "shadow boundary", that is, the discontinuity of the geometric optics solution. Later, others made improvements to GTD and solved its discontinuity problem near the shadow boundary. The improved new theory is called the uniform geometric diffraction theory. For details, please refer to MCNAMARA D A, PISTORIUS C W I, MALHERBE J A G. Introduction to the Uniform Geometrical Theory of Diffraction [M]. Artech House, 1990. Another attempt to express straight edge diffraction began with the Biot-Tolstoy solution. It gives an analytical form of the echo when a spherical wave is incident on a convex dihedral angle. When applying the Biot-Tolstoy solution to acoustic research, Medwin decomposed the echo into a combination of a series of radial functions centered on each point on the straight edge based on the intuition of the Huygens principle of the wave equation, and found that this decomposition helps to explain the diffraction problem of a finite-length straight edge. Combining the Biot-Tolstoy solution and the Medwin decomposition, another description model of straight edge diffraction can be obtained, called the Biot-Tolstoy-Medwin (BTM) model. For details, please refer to SVENSSON U P, FRED R I, VANDERKOOY J. An analytic secondary source model of edge diffraction impulse responses [J]. Journal of the Acoustical Society of America, 1999, 106: 2331-2344.

UTD and BTM are the two most common geometric diffraction models in sound propagation simulation. The UTD (Consistent Geometric Theory of Diffraction) model has been extensively studied because it appeared earlier. For plane and spherical incident waves, ideal half-planes or curved half-planes, various boundary conditions, and other situations, relevant research articles can be found. However, the Kirchhoff approximation is used in the construction of UTD, and the straight edge length is assumed to be infinite. Therefore, in practical applications, the results of UTD are bound to be different from the exact solution. In addition, the representation of the UTD solution is not in the time domain like the general impulse response, but in the frequency domain. This makes it difficult to coordinate with other parts of the geometric sound propagation simulator. Common applications often only select the results of the UTD solution at a few frequency points and discard all phase information. This makes the error of UTD applied to the actual sound propagation simulator much greater than the theoretical error. The Biot-Tolstoy solution on which the UTD model relies is an analytical solution to the infinite length straight edge diffraction problem, and the Medwin decomposition is not an approximation. Therefore, the BTM method can achieve a higher degree of accuracy than UTD in theory. In addition, since the Medwin decomposition decomposes straight edges into a series of radial functions, the BTM method can also handle the diffraction problem of segmented straight edges and even curved edges. BTM is inferior to UTD in that there is less relevant research data and its implementation in existing sound propagation simulators is not as efficient as UTD.

Both UTD and BTM are widely used in sound propagation simulation. At the beginning of this century, an article realized single diffraction based on UTD in an acoustic virtual reality system, and later researchers quickly extended it to the case of high-order diffraction. For more complex scenes, there are also methods to quickly find feasible diffraction paths, but the scene needs to be simplified in advance. The original BTM expression is directly used for sound propagation simulation, but its computational efficiency is too low and it is only suitable for offline calculations of small-scale scenes. There are also methods based on BTM to calculate diffraction solutions for common geometric bodies, and use these solutions to approximate objects in the actual environment. Such as volume diffraction and propagation method.

In addition to UTD and BTM, there are some diffraction models that are not popular. Some methods use the Heisenberg uncertainty principle to explain the diffraction of sound waves and use the diffraction angle probability density function to describe the characteristics of diffracted waves propagating along non-straight lines. This method lacks a physical basis and it is difficult to guarantee the correctness of the results; the recently proposed directional line source model is a potential model. It uses a small number of approximations to explain a large number of straight edge diffraction in different situations. See MENOUNOU P, NIKOLAOU P. Analytical model for predicting edge diffraction in the time domain. [J]. The Journal of the Acoustical Society of America, 2017, 142 6: 3580. However, the application of DLSM in actual sound propagation simulators has not yet appeared.

Therefore, we urgently need a solution to calculate high-precision straight-edge diffraction in large-scale scenes within a geometric framework.

Summary of the invention

The purpose of the present invention is to provide a method for simulating acoustic diffraction phenomena based on path tracking in view of the deficiencies of the prior art.

The present invention is achieved through the following technical solutions:

A method for calculating straight-edge diffraction of acoustic waves based on path tracking comprises the following steps:

(1) reading an acoustic environment state related to diffraction, wherein the acoustic environment state includes the positions of the sound source and the listener, a scene geometry description, and interface features of a geometric surface;

(2) randomly generating a large number of diffraction paths in the scene, wherein the paths are composed of a series of nodes and line segments connecting the nodes in the scene, the starting node of the path is the sound source, the ending node is the listener, and the remaining intermediate nodes are located on the convex straight edges of the geometric body in the scene; except for the nodes constituting the paths, the other parts of the paths do not intersect with the surfaces of the geometric bodies in the scene;

(3) Calculate the final contribution of each path. Assume that the path contains n+1 nodes. The nodes are numbered from _x0 to _xn in sequence from the sound source to the listener. The line segment between the nodes is _xi → _xi+1 . The complete path is _x0 → _x1 →…→ _xn . The total contribution f( _x0 → _x1 →…→ _xn ) is calculated as follows:

Wherein, L ₀ (x ₀ →x) ₁ is the sound source characteristic, M ( _xi → _xi+1 ) is the sound medium characteristic at _xi → _xi+1 , and ρ (xi _-1 → _xi → _xi+1 ) is the radial function representation of the straight edge diffraction at _xi ;

(4) Add the contributions of all paths to the impulse response function: For the path x ₀ →x ₁ →…→x _n , its contribution is known to be f(x ₀ →x ₁ →…→x _n ) in step (3), and its generation probability is recorded as p(x ₀ →x ₁ →…→x _n ), and the propagation delay is t _S , then add t S to the original impulse response function.

Where δ(x) is the unit impulse function;

(5) Output the impulse response function generated in step (4) as the diffraction simulation result.

Specifically, the step (2) includes the following sub-steps:

(2.1) Initialize the path so that the only node it contains is the starting point of the path;

(2.2) Assume that the current path state is x ₀ →x ₁ →…→x _n , and randomly emit a ray from x _n ; if x _n is not the starting point of the path, use a sampling function that matches the straight-edge diffraction radial function to determine the ray's exit angle;

(2.3) The first intersection point of the ray and the scene geometric surface is

And it is called a pseudo-intersection point; if the point does not exist, the intersection fails and go to step (2.8);

(2.4) In step (2.3), the pseudo intersection

Randomly select one of the three edges of the triangle where it is located; let the triangle be Δabc, the selected edge be ab, and the vertex opposite to the edge be c, then use a straight line to connect c and

And extend it to intersect ab at a point, denoted as x _n+1 , called the intersection point;

(2.5) If the edge ab in step (2.4) is not a convex straight edge, or the line segment _xn → _xn+1 intersects the scene geometry at a place other than the two endpoints, the intersection fails and go to step (2.8); otherwise, add the intersection point _xn+1 as a node to the end of the path _x0 → _x1 →…→ _xn to form a new path _x0 → _x1 →…→ _xn+1 ;

(2.6) Calculate the generation probability of the path x ₀ →x ₁ →…→x _n+1 : If the generation probability of the path x ₀ →x ₁ →…→x _n is known to be p(x ₀ →x ₁ →…→x _n ), then the generation probability of the path x ₀ →x ₁ →…→x _n+1 can be calculated as

in,

When the incident path direction at x _n is x _n-1 →x _n , the direction of the generated outgoing ray is

If n = 0, this term must be replaced by the direction of the ray source emitting rays.

The probability;

For rays

Intersects with scene geometry

The probability of; p _C (x _n →x _n+1 ) is the probability correction term;

(2.7) If the path length does not meet the user's requirements, go to step (2.2);

(2.8) The current path is taken as the diffraction path in the scene.

Furthermore, in step (2), paths are generated simultaneously from both ends of the sound source and the listener, the path starting from the sound source is called the forward subpath, and the path starting from the listener is called the backward subpath; thereafter, connections are randomly established between the nodes of the forward and backward subpaths to form a complete path from the sound source to the listener; assuming that the probability of generating the part of the forward subpath from the sound source to the connected point is p _F , the probability of generating the part of the backward subpath from the listener to the connected point is p _B , and the probability of establishing a connection is p _L , then the probability of generating the complete path is p _F ·p _B ·p _L .

Furthermore, in step (5), the contribution of the path needs to be multiplied by a factor C _s for suppressing extreme values:

in

is a variable defined by the user, _Co is the product of the _Co corresponding to each randomly generated segment in the path, and the _Co of each path segment is a parameter that measures the possibility of the path containing the path segment generating an extreme value.

Furthermore, in step (3), the contribution calculation needs to be divided by the product of the elements in [F _j,i ] corresponding to the positions of the nodes connected in the forward and reverse paths; [F _j,i ] is an array of M rows and N columns, where N is the maximum number of nodes in the forward and backward paths; the nodes of the path are stored in an array with the same number of rows and columns, where each row corresponds to a separate path, and the elements in the first column of the array [F _j,i ] are set to 1, and the first column of the node array is the starting point of the path; when the node in the i-1th column has been generated, the generation process of the node in the i-th column of the array includes the following sub-steps:

(5.1) Starting from the i-1th node x _i-2 in each row of the node array, a ray is randomly emitted; if i>1, a sampling function matching the straight-edge diffraction radial function is used to determine the ray's exit angle;

(5.2) The first intersection point of the ray and the scene geometric surface is

And it is called a pseudo-intersection point; if the point does not exist, the intersection fails and go to step (5.7);

(5.3) In step (5.2), the pseudo intersection

And extend it to intersect ab at a point, denoted as x _i-1 , called the intersection point;

(5.4) If the edge ab in step (5.3) is not a convex straight edge, or the line segment x _i-2 → _xi-1 intersects the scene geometry at a place other than the two end points, the intersection fails and go to step (5.7); otherwise, add the intersection point x _i-1 as a node to the end of the path x ₀ →x ₁ →…→ _xi-2 in the current row to form a new path x ₀ →x ₁ →…→ _xi-1 ;

(5.5) Calculate the generation probability of the path x ₀ →x ₁ →…→ _i- x: If the generation probability of the path x ₀ →x ₁ →…→ _i- x is known to be

Then the path

The generation probability calculation method is:

in,

When the incident path direction at x _i-2 is x _i-3 →x _i-2 , the direction of the generated outgoing ray is

If i = 1, this item must be replaced by the direction of the ray source.

The probability;

For rays

Intersects with scene geometry

The probability of; p _C ( _xi-2 → _xi-1 ) is the probability correction term;

(5.6) If the path length does not meet the user's requirements, go to step (5.1) and prepare to generate the next column of nodes;

(5.7) For each node x _i-2 in the i-1th column, count the set of rows S ₊ (x _i-2 ) whose successor nodes are successfully generated and the set of rows S _- (x _i-2 ) whose successor nodes are unsuccessful;

(5.8) Copy the nodes in the i-th column of each row in S ₊ ( _xi-2 ) to the i-th column of each row in S _- ( _xi-2 ); if the node in the j-th row is copied k times, let Fj _,i = (k+1)Fj _,i-1 .

The beneficial effects of the present invention are:

Compared with traditional methods, the method in this paper has higher computational efficiency and does not require pre-processing calculations such as pre-solving wave equations and geometric simplification when dealing with large-scale scenes with complex occlusion relationships. The diffraction model used in this method does not have approximate operations such as Kirchhoff approximation and infinite straight edge assumption, so the result is more accurate. For simpler geometric bodies, the result is completely consistent with the theoretical diffraction solution. In addition, the diffraction algorithm based on path tracking is a natural extension of the traditional sound propagation algorithm based on path tracking, so the method of the present invention is also easier to add to the existing sound propagation simulator.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG1 is a schematic diagram of a rectangular coordinate system in the present invention;

FIG2 is an absolute value image of the straight edge diffraction interface term function of the present invention;

FIG3 is a schematic diagram of the sampling probability distribution of the straight edge diffraction interface item of the present invention;

FIG4 is a schematic diagram of the arrangement of a first-order diffraction test scene;

FIG5 is a schematic diagram of calculation results of the present invention and the control method under the first-order diffraction test scene arrangement;

FIG6 is a schematic diagram of the arrangement of a second-order diffraction test scene;

FIG. 7 is a schematic diagram of calculation results of the present invention and the control method under the second-order diffraction test scene arrangement.

Detailed ways

The present invention provides a method for calculating straight-edge diffraction of sound waves based on path tracking, comprising the following steps:

Wherein, L ₀ (x ₀ →x ₁ ) is the sound source characteristic, M ( _xi → _xi+1 ) is the sound medium characteristic at _xi → _xi+1 , ρ (xi _-1 → _xi →xi ₊₁ ) is the radial function representation of the straight edge diffraction at _xi ;

Where δ(x) is the unit impulse function;

The path tracking-based straight-edge diffraction calculation method for acoustic waves is specifically as follows:

Path tracing: The technique of path tracing in sound propagation simulation is the path tracing algorithm used in sound propagation simulation and the application of diffraction representation in the algorithm.

In the sound propagation simulation, the path tracing algorithm solves the following rendering equation:

Where t represents time. The moment t=0 is the moment when the sound source starts to make sound; L(x′→x, t) is the energy response or impulse response excited by the sound source x′ at position x, that is, the objective function to be solved in this equation; L ₀ (x′→x, t) represents the direct component, that is, the response excited directly at x by the sound source without interacting with the scene geometry; V(x″, x) is the visibility function, which is 1 when x″ and x are visible to each other, and 0 otherwise; ρ(x′→x″→x) is the interface term, which represents the outgoing sound intensity generated near a point x″ on the interface in the direction x″→x when the incident sound direction is x′→x″; M _t (x″, x′, t) is the medium term, which represents the influence of the sound medium between x″ and x′, including energy absorption and propagation delay. The asterisk * is the convolution symbol; Ω represents the point set corresponding to the surface of the scene geometry, and A _Ω is the area measure on the set. When calculating diffraction, for homogeneous media, the following medium terms are used:

Where a is the acoustic attenuation coefficient of the medium, ||x″-x′|| is the length of the vector x″-x′, δ(·) is the unit impulse response function, and c is the speed of sound in the medium.

Let the above integral be the operator T, then the rendering equation can be simplified as L=L ₀ +TL; the Neumann series expansion is used to express the equation:

The operator T represents an interaction between a sound wave and the scene geometry. Since T is essentially an integral, its value can be calculated using the Monte Carlo integration method.

In path tracking, the present invention represents ^TnL0 as the mathematical expectation of the path contribution; a path is _composed of a series of nodes and line segments connecting the nodes in the scene, with the starting node being the sound source, the end node being the listener, and the intermediate nodes being located on the surface of the scene geometry; a path containing n intermediate nodes (n+2 nodes) is called an n-order path, corresponding to _the term ^TnL0 in the series.

From the expression of the dielectric term, it is not difficult to know that the total propagation delay of a path is the sum of the time delays of each segment of the path (length divided by the speed of sound). Therefore, it can be calculated separately. Ignore the time term and only calculate the intensity of the path contribution, using the dielectric term without time parameters

Replace the original medium term M _t in the rendering equation. At this time, the contribution expression of an n-order path x ₀ →x ₁ →…→x _n+1 is as follows:

If the above path is randomly generated, then the mathematical expectation of its contribution needs to be divided by its generation probability when calculating it. Considering the time delay, the contribution expectation can be calculated as follows:

Where t _S is the total time delay of the path. When multiple paths are used to calculate the expectation, the expectation value is the average of the results of each path.

There are many ways to generate the random path; either using unidirectional path tracing or using bidirectional path tracing.

Interface term expression of straight edge diffraction: The key to using path tracking to calculate diffraction is to construct the corresponding interface term. Here, the present invention will give the form of the interface term. As shown in Figure 1, the coordinate system near the straight edge and the representation method and related symbols of the incident and outgoing path directions are given. The present invention calls n in the figure the normal, b the binormal, and t the tangent (the t direction coincides with the straight edge direction, both positive and negative). The three constitute a frame at a point on the straight edge. Half of the dihedral angle at the straight edge is denoted as φ. Since this article only considers convex straight edges, the value of φ should be less than π/2. The incident direction vector _vi and the outgoing direction vector v _o are represented by a series of angles:

In this coordinate system, write the interface term form of straight edge diffraction; first define the following auxiliary variables:

ω _lu =θ _i -π

ω _lr = -θ _i -π + 2φ

ω _ru = θ _i + π

_ωrr ＝ _-θi +π-2φ

At this time, if the geometric surface near a point on the straight edge satisfies the Dirichlet boundary condition, the interface term at this point is:

If the geometric surface satisfies the Neumann boundary condition, the interface term is:

The form of the above interface term can be proved by mathematics. Under other conditions, an interface term of the form aρ _D +bρ _N can be considered, where a and b are two constants.

In path tracing, it is necessary to launch a ray from a point on a straight edge so that the probability distribution of its launch direction is close to the absolute value of the interface term function. This probability distribution needs to be specially designed. To describe this distribution, some auxiliary variables need to be defined. Let

When using Dirichlet boundary conditions, let D = D _D ; when using Neumann boundary conditions, let D = D _N . The remaining variables are defined as follows:

Then, the probability density function of the ray emission direction is

It can be expressed as follows:

Where p _θ is the probability density function of θ _o ,

yes

The probability density function of . The expressions of these two functions are as follows:

The function constructed above satisfies the non-negativity and normalization requirements of the probability density function.

Figure 2 shows that under the Dirichlet condition,

θ _i = 0.4636,

FIG3 is a schematic diagram of the sampling probability distribution of the straight edge diffraction interface term of the present invention, showing the sampling probability density function under the corresponding direction. It can be seen that the shapes of the functions in the two figures are very similar.

Diffraction path generation: In the path tracing algorithm, the common path generation method is to generate each node in the path in sequence. Generally, a ray with a random direction is emitted from the end node of the current path. If it intersects with the scene geometry, the intersection point is used as the successor node of the original end node. However, when simulating diffraction, the middle node of the path must be located on the edge of the geometry, and the general ray intersection cannot guarantee that the intersection point meets this requirement. Therefore, it needs to be modified.

Assume that the scene geometry is described by a triangular mesh. After the ray intersection is successful, the intersection point between the ray and the geometry is called a pseudo-intersection point. Check the edges of the triangle where the pseudo-intersection point is located and determine whether these edges are convex edges. Next, randomly select one of these convex edges. Let the triangle face be Δabc ₁ and the selected edge be ab, then the probability of selecting this edge is recorded as

After that, use a straight line to connect the pseudo intersection point and point _c1 and extend it to intersect ab at one point. This is called the (actual) intersection point corresponding to the pseudo intersection point, and it is used as the successor node of the original end node.

Now let’s calculate the probability of node generation under this method. Let the original path be x ₀ →x ₁ →…→x _n , and the pseudo intersection be

The corresponding intersection point is _xn+1 . ab has two adjacent faces, let the other face different from Δabc ₁ be Δabc _2. Then, all pseudo-intersection points corresponding to _xn+1 are located on the line segments _xn+1c1 _and xn ₊ _1c2 . Next, calculate two auxiliary variables S( _xn , _xn+ _1c1 ) and S( _xn , _xn+ _1c2 ). The calculation method of S( _xn , xn _+1c1 ₎ is as follows:

(3.1) Find the part of x _n+1 c ₁ that is visible (not blocked) to x _n . Since the scene is composed of triangles, the visible part must be a set of non-intersecting line segments; calculate all blocked line segments by finding the triangles that intersect with the triangle Δx _n x _n+1 c ₁ in the scene geometry, and then remove these blocked line segments from x _n+1 c _1. The remaining is the set of visible line segments. These visible line segments are denoted as

in

compared to

Closer to the _c1 end.

(3.2) Calculate S(x _n ,x _n+1 c ₁ ). The formula is

The calculation method of S(x _n ,x _n+1 c ₂ ) is the same as above.

Finally, an estimate of the node generation probability is given:

in,

When the incident path direction at x _n is x _n-1 →x _n , the outgoing direction is generated

The specific probability distribution is determined by the algorithm implementer. If n = 0, this item must be replaced by the direction of the ray source emitting rays.

The probability.

is the intersection of the outgoing ray and the scene

The probability is calculated as follows:

in

For rays

and

The angle between the normal lines of the triangle surface. _{p C} (x _n →x _n+1 ) is the probability correction term, and the calculation formula is as follows:

in

is the area of Δabc ₁ ,

Similarly, since x _n+1 corresponds to

More than one, the above probability itself is not directly equivalent to the probability of generating node _xn+1 , but its expectation is indeed equal to the probability of generating xn ₊₁ .

The path generation probability after adding nodes is as follows:

In the bidirectional path tracking, it is necessary to generate paths from both ends of the sound source and the listener at the same time, and then establish connections between the nodes of the two paths to generate a complete path. At this time, the probability of generating a complete path is _pF · _pB · _pL . Among them, _pF is the probability of generating the part of the forward subpath from the sound source to the connected part, _pB is the probability of generating the part of the backward subpath from the listener to the connected part, and _pL is the probability of establishing a connection.

In addition, the present invention can also adopt multiple importance sampling, which is a method for improving sample quality and is usually used in conjunction with bidirectional path tracking. For a path, multiple importance sampling requires giving the generation probability of this path under all generation methods. Suppose a path _x0 → _x1 →…→ _xn is generated through bidirectional path tracking. Its connecting segment may be any segment in the path, so there are a total of n-1 possible generation methods. The method for calculating the path generation probability under these generation methods. After combining the multiple importance sampling technology on the basis of the straight edge diffraction calculation of the sound wave, in order to calculate the probabilities of various strategies in the multiple importance sampling technology, corresponding forward and backward pseudo-intersections must be generated for each node that is not the starting point or end point, where the forward pseudo-intersection is visible to the predecessor node of the node, and the backward pseudo-intersection is visible to the successor node of the node; specifically as follows:

Assume that the path x ₀ →x ₁ →…→x _n is generated by bidirectional path tracking, and its connecting segment is x _i → _xi+1 , then the generation probability of this path is p(x ₀ →x ₁ →…→ _xi )p( _xi+1 ←xi ₊₂ ←…←x _n )p _L ( _xi →xi ₊₁ ). Here, x _i+1 ← _xi+2 ←…←x _n means the same as x _n →x _n-1 →…→ _xi+1 . The scene positions involved in calculating the path generation probability are listed as follows:

Here is similar to

Points and

The same is true for the pseudo-intersections, which are called "backward pseudo-intersections". They correspond to the nodes on the reverse path. From the generation process of the reverse path, we can see that

is visible to xi ₊₁ , but not guaranteed to be visible to _xi-1 , which is similar to

Such a false intersection (called a forward false intersection) is just the opposite.

By changing the value of i in the above diagram, we can see that in order to obtain the generation probability when the connecting segment is another path segment, each intermediate node needs to have corresponding forward and backward pseudo-intersections. Therefore, it is necessary to actively fill in the missing pseudo-intersections in the diagram. Here we will give a method to fill in the forward pseudo-intersections, and the method to fill in the reverse pseudo-intersections is similar.

Assuming that x _i+1 has no corresponding forward pseudo-intersection, first find the edge where x _i+1 is located, which is denoted as ab. The two adjacent faces adjacent to this edge are denoted as Δabc ₁ and Δabc _2. Then find the visible parts of the line segments x _i+1 c ₁ and x _i+1 c ₂ relative to x _i , and represent them as a set of line segments, each of which is represented by

Represents. Define the function:

Next, for each line segment

If it is a subset of x _i+1 c ₁ , then

Map it to the set [0,1]; if it is a subset of x _i+1 c ₂ , then

Map it to the set [0,1]; secondly, randomly select a point in the mapped line segment set according to uniform distribution; finally, map the selected point back to the original line segment set through _flerp and use it as the pseudo intersection point.

After completing the pseudo-intersection points, the path generation probability under all generation methods can be calculated.

In addition, the present invention can also be used in other optimizations:

When the probability of generating a path is too small, or the path contribution is too large relative to the generation probability, extreme values will appear in the results that may affect the user experience. Therefore, these extreme values need to be identified and removed.

Given a path x ₀ →x ₁ →…→x _n+1 , define auxiliary variables as follows:

Where p _V and

The definition of is the same as in Section 3. If i = 0, p _V also needs to be replaced according to the method in Chapter 3. This variable represents the possibility of the path containing this segment to produce an extreme value. The _Co value of the entire path is the product of the _Co values of each segment of the path. If the path is composed of multiple sub-paths (such as bidirectional path tracking), its _Co value is the product of the _Co values of each sub-path. Next, multiply the final contribution of the path by a value _Cs :

here

is a user-defined variable (usually between 100-1000). The multiplier above will lower

The contribution of all paths is reduced, thereby reducing the generation of extreme values.

If the efficiency of the method for generating the successor node of the path is not high enough, in order to improve the utilization efficiency of the nodes in the bidirectional path tracking, the following method can be used: First, assume that there are M paths generated at a time in the simulation algorithm, and the maximum order of the path is N-2. Therefore, these paths have a total of M×N nodes. These nodes, regardless of whether they are actually successfully generated, are arranged in an array of M rows and N columns. For each node in the array, define a variable F _j,i ,j＝1,2, _… M,i＝1,2,…N, called the splitting degree, and let F _j,1 ＝1. The nodes in the array need to be generated row by row. For the nodes in the i-th column of the array, the specific generation process is as follows: try to generate all the nodes in the i-th column by the method in Section 3; for each node x in the i-1th column, count the row set S ₊ (x) where its successor node is successfully generated and the row set S _- (x) where it fails to generate; copy the nodes in the i-th column of each row in S ₊ (x) to the i-th column of each row in S _- (x). If the node in the jth row is replicated k times, let F _j,i = (k+1)F _j,i-1 . Let the expected contribution of the path to this node be L _j,i , and it is recommended that the number of replications of each node be proportional to L _j,i /F _j,i-1 .

When establishing forward and reverse path connections and generating the final path, the expected contribution of the final path should be divided by the product of the two split degrees corresponding to the forward and reverse path connection nodes. In this way, the empty spaces in the array can be fully utilized, storage space can be saved, and the node connection efficiency can be improved.

Implementation Examples

The inventor implemented the algorithm described above on a computer equipped with an Intel i7 2.60GHz quad-core CPU and 12GB memory. The specific technology selected was bidirectional path tracking combined with multiple importance sampling, and all optimization techniques in Section 5 of the implementation process description were used. The inventor compared its results with the simulation results using the UTD and BTM models.

Figures 4 and 6 show the first-order and second-order diffraction scenarios for comparison. The scenario in Figure 4 contains a straight dihedral with a straight side length of 2m, while the scenario in Figure 6 contains two parallel right dihedrals. The path from the sound source must undergo two diffractions before it can reach the listener. Both scenario geometries use Dirichlet boundary conditions. Figures 5 and 7 show the impulse responses calculated using different methods in the scenario. The methods used here include BTM, UTD, the method of this paper, and the results reconstructed by selecting 8 frequency points in UTD (located between 62.5-8000Hz, distributed in an exponential interval). The error of the UTD method in calculating the second-order and above diffraction is very large, which is not shown in Figure 7. It can be seen that the quality of the calculation results of the algorithm in this paper in the case of first to second order diffraction is not significantly different from that of the existing methods, and is significantly better than the common simplified UTD scheme (only calculating the results at 8 frequencies). Since the BTM method is accurate for the solution of the first-order diffraction, this comparison also verifies the correctness of the method in the present invention. For large-scale complex scenarios, this method shows near real-time computing efficiency. When 12,000 paths are used, for a model containing 22,974 triangular faces, the time overhead for calculating the diffraction impulse response of the present invention is about 140 ms. For a model containing 279,133 faces, the time overhead is about 205 ms. It can be seen that the method of the present invention is suitable for fast calculation of diffraction results in complex scenes.

The above-described embodiments express specific implementation methods of the present invention, and the descriptions thereof are relatively specific and detailed, and are intended to help understand the method of the present invention and its core concept, but they cannot be understood as limiting the scope of the present invention. It should be pointed out that, for those of ordinary skill in the art, several variations and improvements can be made without departing from the concept of the present invention, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the present invention patent shall be subject to the attached claims.

Claims

A method for calculating straight-edge diffraction of acoustic waves based on path tracking, characterized in that it comprises the following steps:

(1) reading an acoustic environment state related to diffraction, wherein the acoustic environment state includes the positions of the sound source and the listener, a scene geometry description, and interface features of a geometric surface;

(2) randomly generating a large number of diffraction paths in the scene, wherein the paths are composed of a series of nodes and line segments connecting the nodes in the scene, the starting node of the path is the sound source, the ending node is the listener, and the remaining intermediate nodes are located on the convex straight edges of the geometric body in the scene; except for the nodes constituting the paths, the other parts of the paths do not intersect with the surfaces of the geometric bodies in the scene;

(3) Calculate the final contribution of each path. Assume that the path contains n+1 nodes. The nodes are numbered from x0 to xn in sequence from the sound source to the listener. The line segment between the nodes is xi → xi+1 . The complete path is x0 → x1 →…→ xn . The total contribution f( x0 → x1 →…→ xn ) is calculated as follows:

Wherein, L 0 (x 0 →x) 1 is the sound source characteristic, M ( xi → xi+1 ) is the sound medium characteristic at xi → xi+1 , and ρ (xi -1 → xi → xi+1 ) is the radial function representation of the straight edge diffraction at xi ;

(4) Add the contributions of all paths to the impulse response function: For the path x 0 →x 1 →…→x n , its contribution is known to be f(x 0 →x 1 →…→x n ) in step (3), and its generation probability is recorded as p(x 0 →x 1 →…→x n ), and the propagation delay is t S , then add t S to the original impulse response function.
Where δ(x) is the unit impulse function;

(5) Output the impulse response function generated in step (4) as the diffraction simulation result.
The method for calculating straight-edge diffraction of acoustic waves based on path tracking according to claim 1, characterized in that step (2) comprises the following sub-steps:

(2.1) Initialize the path so that the only node it contains is the starting point of the path;

(2.2) Assume that the current path state is x 0 →x 1 →…→x n , and randomly emit a ray from x n ; if x n is not the starting point of the path, use a sampling function that matches the straight-edge diffraction radial function to determine the ray's exit angle;

(2.3) The first intersection point of the ray and the scene geometric surface is
And it is called a pseudo-intersection point; if the point does not exist, the intersection fails and go to step (2.8);

(2.4) In step (2.3), the pseudo intersection
Randomly select one of the three edges of the triangle where it is located; let the triangle be Δabc, the selected edge be ab, and the vertex opposite to the edge be c, then use a straight line to connect c and
And extend it to intersect ab at a point, denoted as x n+1 , which is called the intersection point;

(2.5) If the edge ab in step (2.4) is not a convex straight edge, or the line segment xn → xn+1 intersects the scene geometry at a place other than the two endpoints, the intersection fails and go to step (2.8); otherwise, add the intersection point xn+1 as a node to the end of the path x0 → x1 →…→ xn to form a new path x0 → x1 →…→ xn+1 ;

(2.6) Calculate the generation probability of the path x 0 →x 1 →…→x n+1 : If the generation probability of the path x 0 →x 1 →…→x n is known to be p(x 0 →x 1 →…→x n ), then the generation probability of the path x 0 →x 1 →…→x n+1 can be calculated as
in,
When the incident path direction at x n is x n-1 →x n , the direction of the generated outgoing ray is
If n = 0, this term must be replaced by the direction of the ray source emitting rays.
The probability;
For rays
Intersects with scene geometry
The probability of; p C (x n →x n+1 ) is the probability correction term;

(2.7) If the path length does not meet the user's requirements, go to step (2.2);

(2.8) The current path is taken as the diffraction path in the scene.
The method for calculating straight-edge diffraction of sound waves based on path tracing according to claim 1 is characterized in that, in the step (2), paths are generated simultaneously from both ends of the sound source and the listener, the path starting from the sound source is called the forward subpath, and the path starting from the listener is called the backward subpath; thereafter, connections are randomly established between the nodes of the forward and backward subpaths to form a complete path from the sound source to the listener; assuming that the probability of generating the part of the forward subpath from the sound source to the connected part is pF , the probability of generating the part of the backward subpath from the listener to the connected part is pB , and the probability of establishing the connection is pL , then the probability of generating the complete path is pF · pB · pL .
The method for calculating straight-edge diffraction of acoustic waves based on path tracking according to claim 1 is characterized in that, in step (5), the contribution of the path needs to be multiplied by a factor C s for suppressing extreme values:

in
is a variable defined by the user, Co is the product of the Co corresponding to each randomly generated segment in the path, and the Co of each path segment is a parameter that measures the possibility of the path containing the path segment generating an extreme value.
According to the path tracking-based straight-edge diffraction calculation method of sound waves described in claim 3, it is characterized in that in the step (3), the calculated contribution needs to be divided by the product of the elements in [F j,i ] corresponding to the positions of the nodes connecting the forward and reverse paths; the [F j,i ] is an array of M rows and N columns, where N is the maximum number of nodes in the forward and backward paths; the nodes of the path are stored in an array with the same number of rows and columns, each row of which corresponds to a separate path, and the elements in the first column of the array [F j,i ] are set to 1, and the first column of the node array is the starting point of the path; when the nodes in the i-1th column have been generated, the generation process of the nodes in the i-th column in the array includes the following sub-steps:

(5.1) Starting from the i-1th node x i-2 in each row of the node array, a ray is randomly emitted; if i>1, a sampling function matching the straight-edge diffraction radial function is used to determine the ray's exit angle;

(5.2) The first intersection point of the ray and the scene geometric surface is
And it is called a pseudo-intersection point; if the point does not exist, the intersection fails and go to step (5.7);

(5.3) In step (5.2), the pseudo intersection
Randomly select one of the three edges of the triangle where it is located; let the triangle be Δabc, the selected edge be ab, and the vertex opposite to the edge be c, then use a straight line to connect c and
And extend it to intersect ab at a point, denoted as x i-1 , called the intersection point;

(5.4) If the edge ab in step (5.3) is not a convex straight edge, or the line segment x i-2 → xi-1 intersects the scene geometry at a place other than the two end points, the intersection fails and go to step (5.7); otherwise, add the intersection point x i-1 as a node to the end of the path x 0 →x 1 →…→ xi-2 in the current row to form a new path x 0 →x 1 →…→ xi-1 ;

(5.5) Calculation path
The generation probability of: If the path is known
The probability of generating
Then the path
The generation probability calculation method is:

in,
When the incident path direction at x i-2 is x i-3 →x i-2 , the direction of the generated outgoing ray is
If i = 1, this item must be replaced by the direction of the ray source.
The probability;
For rays
Intersects with scene geometry
The probability of; p C ( xi-2 → xi-1 ) is the probability correction term;

(5.6) If the path length does not meet the user's requirements, go to step (5.1) and prepare to generate the next column of nodes;

(5.7) For each node x i-2 in the i-1th column, count the set of rows S + (x i-2 ) whose successor nodes are successfully generated and the set of rows S - (x i-2 ) whose successor nodes are unsuccessful;

(5.8) Copy the nodes in the i-th column of each row in S + ( xi-2 ) to the i-th column of each row in S - ( xi-2 ); if the node in the j-th row is copied k times, let Fj ,i = (k+1)Fj ,i-1 .