Wave calculation program, device and method for calculating boundaries for each frequency band

Abstract

To provide a wave motion calculation program which can suppress a deviation from a measured value of a calculation value that is generated depending on a frequency band.SOLUTION: A wave motion calculation program of calculating a wave motion relevant amount being an amount relevant to a wave motion in a space having a prescribed boundary by using the finite differential time region method causes a computer to function as: calculation means for each frequency band which calculates a calculation object amount related to calculation of a wave motion in a frequency band in a space by using the boundary calculating formula set for the frequency band and the discretization formula relevant to the wave motion for each frequency band in the plurality of frequency bands in a frequency range relevant to the wave motion; and wave motion relevant amount calculation means which calculates the wave motion relevant amount in the space from the calculation object amount in the frequency band in the space calculated for each frequency band. Here, the wave motion is the sound wave; the calculation object amount is the amount relevant to the sound pressure; and the wave motion relevant amount can be an impulse response between the set sound source position and the sound reception position.SELECTED DRAWING: Figure 1

Description

[Technical field]

[0001] The present invention relates to a numerical simulation technique for waves propagating through space. [Background technology]

[0002] Numerous techniques have been developed to calculate wave-related quantities, which are quantities related to waves propagating through space, by using numerical simulations as a form of numerical analysis using computers. For example, the impulse response of a sound field through which sound waves propagate is an important indicator of spatial acoustic characteristics, and techniques to calculate this impulse response by numerical simulations are being actively developed.

[0003] As examples of such sound field numerical simulation techniques, Patent Document 1 and Non-Patent Document 1 disclose methods for calculating impulse responses using the ray method and image method. These methods do not take into account the wave properties of sound waves, and therefore the discrepancy between the calculated values ​​and the measured values ​​becomes significant, especially in the low frequency band.

[0004] In response to this, Non-Patent Documents 2 and 3 disclose a technique for calculating impulse responses using the Finite Difference Time Domain (FDTD) method. The FDTD method is a numerical simulation method that has been widely used in electromagnetic field analysis, as disclosed in Patent Document 2, for example. Non-Patent Documents 2 and 3 use the Compact Explicit (CE)-FDTD method, which applies a discretization formula that refers to more directions than the normal FDTD method. This improves the accuracy of the calculation results, making it possible to calculate impulse responses up to the sampling frequency, for example. [Prior art documents] [Patent documents]

[0005] [Patent Document 1] Japanese Patent Application Publication No. 9-292885 [Patent Document 2] JP 2006-139723 A [Non-patent literature]

[0006] [Non-Patent Document 1] Koji Ishida, "Visualization of hall acoustics using imaginary sound sources", Noise Control, Vol. 22, No. 1, pp.31-33, 1998. [Non-Patent Document 2] Takao Tsuchiya, "Sound Field Simulation and Sound Space Rendering", IEICE Fundamentals Review, Vol. 10, No. 3, pp. 206-218, 2017 [Non-Patent Document 3] Konrad Kowalczyk and Maarten van Walstijn, “Room Acoustics Simulation Using 3-D Compact Explicit FDTD Schemes”, IEEE Transactions on Audio, Speech, and Language Processing, Vol.19, Issue.1, pp.34-36, 2011 Summary of the Invention [Problem to be solved by the invention]

[0007] However, in the conventional techniques disclosed in the above-mentioned Non-Patent Documents 2 and 3, although the accuracy is improved by using the CE-FDTD method, there is still a problem that the reproduction accuracy (degree of agreement with actual measured values) decreases depending on the frequency band.

[0008] In fact, the inventors of the present application calculated the reverberation time of the impulse response of a certain soundproof room using the CE-FDTD method and confirmed that the calculated value deviated from the actual measured value in frequency bands lower and higher than the frequency band used to set the boundary conditions, with the deviation being particularly large in the low frequency band.

[0009] Such deviations from actual measurements are often observed in the calculation results of sound-wave-related quantities using the conventional FDTD method, and also in the calculation results of FDTD methods related to waves other than sound waves, such as electromagnetic waves.

[0010] Therefore, an object of the present invention is to provide a wave calculation program, device, and method that can suppress the deviation of calculated values ​​from actual measured values ​​that occurs depending on the frequency band. [Means for solving the problem]

[0011] According to the present invention, there is provided a program for calculating a wave-related quantity, which is a quantity related to a wave in a space having a predetermined boundary, by using a finite-difference time-domain method, comprising the steps of: A frequency band calculation means for calculating a calculation target quantity related to the calculation of the wave in the frequency band in the space by using a boundary calculation formula set for each frequency band among a plurality of frequency bands included in the frequency range related to the wave and a discretization formula related to the wave; a wave-related quantity calculation means for calculating the wave-related quantity in the space from the calculation target quantity in the frequency band in the space calculated for each frequency band; A wave calculation program is provided that causes a computer to function as such.

[0012] In this wave calculation program according to the present invention, it is also preferable that the frequency band calculation means specifically applies, for each frequency band, the value of the boundary quantity set for the frequency band to the boundary calculation formula which includes as a parameter a boundary quantity, which is a quantity related to the absorption or impedance of the boundary, and calculates, using the boundary calculation formula for the frequency band to which the value of the boundary quantity has been applied, and the discretization formula, the calculation target quantity in the frequency band in the space, including the calculation target quantity in the frequency band at a position outside the boundary adjacent to the boundary, and the calculation target quantity at the boundary calculated using the calculation target quantity at a position outside the boundary.

[0013] In addition, as one embodiment of the wave calculation program according to the present invention, it is also preferable that the frequency band calculation means calculates, for each frequency band, the calculation target quantity in the frequency band in the space using the boundary calculation formula and the discretization formula corresponding to the sampling frequency set for that frequency band.

[0014] Furthermore, as another embodiment of the wave calculation program according to the present invention, it is also preferable that the wave-related quantity calculation means calculates, for each frequency band, the wave-related quantity in that frequency band from the calculation target quantity in that frequency band, and calculates the wave-related quantity in the space by adding up the calculated wave-related quantities in each frequency band, or by adding up the wave-related quantities in a predetermined frequency band after assigning a larger or smaller weight to the wave-related quantity in that frequency band.

[0015] As yet another embodiment of the wave calculation program according to the present invention, it is also preferable that the wave-related quantity calculation means calculates, for each frequency band, the wave-related quantity in that frequency band from the calculation target quantity in that frequency band, applies band-pass filtering to the calculated wave-related quantity in each frequency band with a passband set to match the frequency band, and calculates the wave-related quantity in the space by adding up the wave-related quantities in each frequency band that have been subjected to band-pass filtering, or by adding up the wave-related quantities in a predetermined frequency band after applying a larger or smaller weight to the wave-related quantities.

[0016] Furthermore, as yet another embodiment of the wave calculation program according to the present invention, it is also preferable that the wave-related quantity calculation means calculates, for each frequency band, the wave-related quantity in that frequency band from the calculation target quantity in that frequency band, performs Fourier transform processing on the calculated wave-related quantity in each frequency band, determines a Fourier-transformed wave-related quantity in the space from the wave-related quantity in each wavenumber band subjected to the Fourier transform processing, and performs inverse Fourier transform processing on the determined Fourier-transformed wave-related quantity to calculate the wave-related quantity in the space.

[0017] Furthermore, it is also preferable that the finite-difference time-domain method in the wave calculation program according to the present invention is a Compact-Explicit Finite-Difference Time-Domain (CE-FDTD) method, the boundary calculation formula is set for each reference direction among a plurality of reference directions, and the discretization formula is an formula corresponding to the plurality of reference directions.

[0018] Furthermore, it is also preferable that the finite-difference time domain method in the wave calculation program according to the present invention is a method in which the target variable in the governing equation includes a vector quantity, or a method in which the target variable in the governing equation is a scalar quantity.

[0019] In addition, in the wave calculation program according to the present invention, it is also preferable that the wave is a sound wave, the quantity to be calculated is a quantity related to sound pressure, and the wave-related quantity is an impulse response between a set sound source position and a sound receiving position.

[0020] According to the present invention, there is also provided an apparatus for calculating a wave-related quantity, which is a quantity related to a wave in a space having a predetermined boundary, by using a finite-difference time-domain method, comprising: A frequency band calculation means for calculating a calculation target quantity related to the calculation of the wave in the frequency band in the space by using a boundary calculation formula set for each frequency band among a plurality of frequency bands included in the frequency range related to the wave and a discretization formula related to the wave; a wave-related quantity calculation means for calculating the wave-related quantity in the space from the calculation target quantity in the frequency band in the space calculated for each frequency band; A wave calculation device is provided having:

[0021] According to the present invention, there is further provided a method for calculating a wave-related quantity, which is a quantity related to a wave in a space having a predetermined boundary, by using a finite-difference time-domain method, comprising the steps of: A step of calculating a calculation target quantity related to the calculation of the wave in the frequency band in the space, for each frequency band among a plurality of frequency bands included in the frequency range related to the wave, using a boundary calculation formula set for the frequency band and a discretization formula related to the wave; calculating the wave-related quantity in the space from the calculation target quantity in the frequency band in the space calculated for each frequency band; There is provided a computer-implemented wave calculation method comprising the steps of: Effect of the Invention

[0022] According to the wave calculation program, device, and method of the present invention, it is possible to suppress the deviation of the calculated value from the actual measured value that occurs depending on the frequency band. [Brief description of the drawings]

[0023] [Figure 1] FIG. 1 is a functional block diagram showing a functional configuration of a wave calculation device according to an embodiment of the present invention. [Diagram 2] FIG. 1 is a schematic diagram for explaining a reference direction in the CE-FDTD method. [Diagram 3] FIG. 1 is a schematic diagram for explaining a boundary calculation formula in the CE-FDTD method. [Figure 4] FIG. 1 is a schematic diagram showing a staggered grid and a colcade grid in a vector variable FDTD method. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0024] Hereinafter, an embodiment of the present invention will be described in detail with reference to the drawings.

[0025] [Wave calculation device] FIG. 1 is a functional block diagram showing the functional configuration of a wave calculation device according to an embodiment of the present invention.

[0026] The wave calculation device 1 according to one embodiment of the present invention shown in Fig. 1 is a device that calculates a "wave-related quantity", which is a quantity related to waves in a space having a predetermined boundary, by using the finite difference time domain (FDTD) method. Here, the wave calculation device 1 of this embodiment employs a sound wave as a wave, calculates sound pressure as a "quantity to be calculated" related to the wave calculation, and further uses the calculation result of this sound pressure to calculate an impulse response as a "wave-related quantity" in a sound field through which this sound wave propagates.

[0027] In order to perform such a wave calculation process, the wave calculation device 1 specifically: (A) a frequency band (quantity to be calculated) calculation unit 111 that calculates a "quantity to be calculated" (sound pressure in this embodiment) for each frequency band in a space having a predetermined boundary, using a "boundary calculation formula" set for that frequency band and a "discretization formula" for the wave (sound wave) for each frequency band included in a frequency range related to the wave (sound wave in this embodiment); (B) a wave-related quantity calculation unit 112 that calculates a "wave-related quantity" (impulse response in this embodiment) in this space from the "quantity to be calculated" (sound pressure) in the frequency band in this space calculated for each frequency band; It has.

[0028] Conventionally, the calculated values ​​of wave-related quantities derived using the FDTD method tend to deviate from the measured values ​​in frequency bands lower and higher than the frequency band used to set the boundary conditions. In contrast, as can be seen from the above configurations (A) and (B), the wave calculation device 1 uses a "boundary calculation formula" set for each frequency band to calculate the "quantity to be calculated" in each frequency band, and then calculates the "wave-related quantity" in a space having a predetermined boundary, which is the target. As a result, it is possible to suppress the deviation of the calculated value from the measured value caused by the frequency band, and improve the accuracy of the wave calculation process.

[0029] In this embodiment, as described above, the wave is a sound wave, and the calculation target quantity and the wave-related quantity are sound pressure (or a quantity related to sound pressure) and impulse response, respectively, but of course the wave, calculation target quantity and wave-related quantity according to the present invention are not limited to these. For example, it is also possible to calculate acoustic intensity as a wave-related quantity by using sound pressure and particle velocity as the calculation target quantity. Furthermore, it is also possible to adopt an electromagnetic wave as a wave, and calculate an electromagnetic field gradient or electromagnetic field energy (Poynting vector) as a wave-related quantity by using electromagnetic field strength (electric field strength, magnetic field strength, etc.) as the calculation target quantity.

[0030] In any case, the wave calculation process according to the present invention can handle a wide variety of other waves, calculation target quantities, and wave-related quantities, but is particularly suitable for calculating wave-related quantities related to waves propagating in a space whose boundary conditions change significantly depending on the frequency.

[0031] Furthermore, the space in which the wave propagates is not limited to three dimensions, and for example, the FDTD method for waves propagating in a two-dimensional space (plane) can also be used. Furthermore, the FDTD method used in this embodiment is a CE (Compact explicit)-FDTD method that applies a discretization formula with reference to more directions than the normal FDTD method, as described later. However, the wave calculation process according to the present invention can also be performed using the normal (not multi-directional) FDTD method.

[0032] Furthermore, in the wave calculation process according to the present invention, (a) "A method in which the target variable in the governing equation includes a vector quantity" (hereinafter referred to as the vector variable type), and (b) "A method in which the variable of interest in the governing equation is a scalar quantity" (hereinafter referred to as the scalar variable type) It is possible to use either of the above FDTD methods. The CE-FDTD method is an improved version of the scalar variable FDTD method using multi-directional references.

[0033] Hereinafter, the FDTD method including the above-mentioned "boundary calculation formula" and "discretization formula" in (A) will be described with reference to the drawings. Incidentally, the FDTD method including the CE-FDTD method is described in detail in Non-Patent Document 2 and Non-Patent Document: Osamu Yamashita, Takao Tsuchiya, Yukio Iwaya, Makoto Otani, and Yasushi Inoguchi, "Reflective boundary condition with arbitrary boundary shape for compact-explicit finite-difference time-domain method", Japanese Journal of Applied Physics, Volume 54, Number 7S1, 2015.

[0034] <Explanation of the FDTD method> First, a simulation of the sound pressure (distribution) in a three-dimensional space of a certain medium (e.g., air) with a predetermined boundary using the CE-FDTD method (scalar variable type) when the wave is a sound wave and the quantity to be calculated is the sound pressure will be explained.

[0035] In this case, as the governing equation, the following equation (1) c -2 ×∂ 2 p / ∂t 2 =∇ 2 p is used. Here, p is the sound pressure and c is the speed of sound (in that medium). Incidentally, the above equation (1) is the wave equation of the sound pressure derived from the wave equation and the continuity equation of the sound wave (sound field) described later using the state equation p = c 2 ×ρ (ρ is the medium density).

[0036] Next, by discretizing the above equation (1), the discretization formula of the CE-FDTD method for the sound pressure (2) p i,j,k n+1 = d1×(p i+1,j,k n +p i-1,j,kn +p i,j+1,k n +p i,j-1,k n +p i,j,k+1 n +p i,j,k-1 n ) +d2×(p i+1,j+1,k n +p i+1,j-1,k n +p i+1,j,k+1 n +p i+1,j,k-1 n +p i,j+1,k+1 n +p i,j+1,k-1 n +p i,j-1,k+1 n +p i,j-1,k-1 n +p i-1,j+1,k n +p i-1,j-1,k n +p i-1,j,k+1 n +p i-1,j,k-1 n ) +d3×(p i+1,j+1,k+1 n +p i+1,j-1,k+1 n +p i+1,j+1,k-1 n +p i+1,j-1,k-1 n +p i-1,j+1,k+1 n +p i-1,j-1,k+1 n +p i-1,j+1,k-1 n +p i-1,j-1,k-1 n ) +d4×p i,j,k n -p i,j,k n-1 Here, d1=χ 2 ×(1-4a+4b) d2=χ 2 ×(a-2b) d3=χ 2 ×b d4=2×(1-3×χ 2 +6a×χ 2-4b×χ 2 ) is obtained.

[0037] In the above formula (2), p i,j,k n is the sound pressure at time n×Δt at a three-dimensional grid position (i×Δ, j×Δ, k×Δ) in space, where Δ is the grid interval and Δt is the time step. Also, χ is expressed by the following formula: (3) χ=c×Δt / Δ is the CFL (Courant-Friedrichs-Lewy) number that represents the stability condition, defined by the following. Furthermore, a and b are parameters for controlling the accuracy, and are set to appropriate values ​​depending on the method. For example, if a=1 / 4 and b=1 / 16, the cutoff frequency will match the Nyquist frequency, improving accuracy.

[0038] In addition, in the above formula (2), The term d1 corresponds to the reference along the coordinate axis in grid space (hereafter abbreviated as SLF (Standard Leapflog)), The term d2 corresponds to the diagonal reference (within the plane) in the grid space (hereafter abbreviated as CCP (Cubic Close packed)), The term d3 corresponds to the reference of the solid diagonal direction in the grid space (hereafter abbreviated as OCTA (OCTAhedral)).

[0039] Here, the reference directions in the grid space for SLF (Fig. 2(A)), CCP (Fig. 2(B)) and OCTA (Fig. 2(C)) are shown in Fig. 2. As shown above, the discretization equation for the CE-FDTD method (Equation (2) above) corresponds to the multiple reference directions that have been set.

[0040] Incidentally, if a and b are set to zero (a = 0, b = 0), then both d2 and d3 become zero, and as a result, the above equation (2) becomes a discretized equation for the normal scalar variable FDTD method (with one type of reference direction), including the term d1 corresponding to the SLF.

[0041] Next, a method for calculating the sound pressure at the boundary of this three-dimensional space will be described. Figure 3 is a schematic diagram for explaining the boundary calculation formula in the CE-FDTD method.

[0042] As shown in Figure 3, to calculate the sound pressure value at a boundary position (for example, the position of the black mark in Figure 3) in the grid space (the space shown by the solid line in Figure 3) where the sound pressure should be calculated, the sound pressure value (as an estimated value) at the virtual grid position outside the boundary adjacent to this boundary (for example, the position of the white mark in Figure 3) is also required. Therefore, a boundary calculation formula that includes the sound pressure at these virtual grid positions is determined by applying the impedance boundary condition, and the sound pressure at the boundary is calculated using this boundary calculation formula.

[0043] Specifically, the following equation is one of the boundary calculation formulas corresponding to SLF (see white circle in Figure 3): (4) p i+1,j,k n+1 =p i,j,k n +A×(p i,j,k n+1 -p i+1,j,k n ) Here, A=(Z * ×χ-1) / (Z * ×χ+1) Z * =Z / ρc Here, the above equation (4) expresses the sound pressure p i+1,j,k n+1 However, similar boundary calculation formulas can be defined for sound pressure at other SLF-equivalent positions. Incidentally, the boundary calculation formulas for the normal scalar variable FDTD method (with one reference direction) are composed only of those corresponding to such SLFs.

[0044] In the above formula (4), Z * is the normalized impedance at the boundary, and Z is the acoustic impedance at the boundary. The relationship between this acoustic impedance Z and the normal incidence sound absorption coefficient α and the normal incidence sound reflection coefficient r of the boundary is expressed by the following equation: (5) r = ((Z-Z0) / (Z+Z0)) 2 (6) α=1-r Here, Z0 is the acoustic impedance of the medium (e.g., air). Therefore, if you specify the acoustic impedance Z of the boundary, or the sound absorption coefficient α or reflectance r as the boundary condition, the above parameter A is determined, and the boundary calculation formula (4) can be used.

[0045] Furthermore, the following equation is used as one of the boundary calculation formulas for CCP (see the white square in Figure 3). (7) p i+1,j+1,k n+1 =p i,j,k n +(2 0.5 xZ * χ-1) -1 ×(p i+1,j,k n +p i,j+1,k n -p i+1,j,k n+1 -p i,j+1,k n+1 ) +(2 0.5 xZ * χ-1)×(2 0.5 xZ * χ-1) -1 ×(p i,j,k n+1 -p i+1,j+1,k n ) Here, the above equation (7) expresses the sound pressure p i+1,j+1,k n+1 However, similar boundary calculation formulas can be prescribed for sound pressure at other CCP-equivalent positions.

[0046] Furthermore, as one of the boundary calculation formulas corresponding to COTA (see the double white circle in Figure 3), the following formula is used: (8) p i+1,j+1,k+1 n+1 =p i,j,k n +(3 0.5 -Z* χ)×(3 0.5 +3Z * χ) -1 ×(p i,j+1,k+1 n +p i+1,j+1,k n +p i+1,j,k+1 n -p i+1,j,k n+1 -p i,j+1,k n+1 -p i,j,k+1 n+1 ) +(3 0.5 -Z * χ)×(3 0.5 +3Z * χ) -1 ×(p i+1,j,k n +p i,j+1,k n +p i,j,k+1 n -p i+1,j+1,k n+1 -p i+1,j,k+1 n+1 -p i,j+1,k+1 n+1 ) +(3 0.5 xZ * χ-1)×(3 0.5 xZ * χ+1) -1 × ×(p i,j,k n+1 -p i+1,j+1,k+1 n ) Here, the above equation (8) expresses the sound pressure p i+1,j+1,k+1 n+1 However, similar boundary calculation formulas can be prescribed for sound pressure at other COTA-equivalent positions.

[0047] As described above, it can be seen that the boundary calculation formulas for the CE-FDTD method (for example, the above formulas (4), (7), and (8)) are set for each of multiple reference directions, and are uniquely determined by setting the boundary quantity, i.e., the acoustic impedance Z, or the boundary normal incidence sound absorption coefficient α or the boundary reflectance coefficient r, to a predetermined value in the space of a predetermined medium (i.e., a predetermined ρ value or c value). Note that the fact that the boundary calculation formula is uniquely determined by setting the acoustic impedance Z naturally also applies to the normal FDTD method (with one type of reference direction).

[0048] Next, a simulation using the above-mentioned vector variable type FDTD method will be explained. Specifically, the wave motion is assumed to be a sound wave, and the wave-related quantity is assumed to be sound pressure (scalar quantity) and particle velocity (vector quantity). In addition, the space in which the sound wave propagates (sound field) is assumed to be a two-dimensional space with a predetermined boundary in order to simplify the expression. In this case, the governing equation is the following equation: (9) ∂p / ∂t+ρc 2 ×∇u=0 (10) ∂u / ∂t+ρ -1 ×∇ρ=0 Incidentally, the above equations (9) and (10) are the continuity equation and the wave equation of the sound wave (sound field), respectively, and u in the equations is the particle velocity vector (u=(u_x, u_y)).

[0049] Then, by discretizing the governing equations (9) and (10), the discretized equations for the sound pressure p and particle velocity u of the vector variable type FDTD method are obtained: (11) p i,j n+1 =p i,j n -(ρc 2 Δt / Δ)×(u_x i+0.5,j n+0.5 -u_x i-0.5,j n+0.5 +u_y i,j+0.5 n+0.5 -u_y i,j-0.5 n+0.5 ) (12) u_x i+0.5,j n+0.5 =u_xi+0.5,j n-0.5 -(Δt / ρΔ)×(p i+1,j n -p i,j n ) (13) u_y i,j+0.5 n+0.5 =u_y i,j+0.5 n-0.5 -(Δt / ρΔ)×(p i,j+1 n -p i,j n ) Further, the boundary calculation formula is, for example, the following formula: (14) u_x i+0.5,j n+0.5 =p i,j n / Z (15) u_y i,j+0.5 n+0.5 =p i,j n / Z It can be obtained as follows.

[0050] Incidentally, the above-mentioned discretization formula is the result of discretization on the staggered grid shown in Figure 4(A). In this staggered grid, the sound pressure p and the particle velocity u are defined alternately in the spatial direction with a shift of Δ×0.5, and also in the time direction with a shift of Δt×0.5. In the vector variable FDTD method, this staggered grid is usually used to prevent a decrease in accuracy. However, it is of course possible to implement the vector variable FDTD method using the colcade grid shown in Figure 4(B). In addition, the boundary calculation formula of the vector variable FDTD method using this colcade grid also includes the medium density ρ, the sound speed c in the medium, and the acoustic impedance Z of the boundary.

[0051] In this way, it can be understood that the boundary calculation formulas (e.g., the above formulas (14) and (15)) of the vector variable type FDTD method are also uniquely determined by setting the boundary quantity, i.e., the boundary acoustic impedance Z, or the boundary normal incidence sound absorption coefficient α or the boundary normal incidence sound reflectance r, to a predetermined value in the space of a predetermined medium (i.e., a predetermined ρ value or c value).

[0052] Various types of FDTD methods including their boundary conditions have been described in detail above, but any of these types of FDTD methods can be adopted in the wave calculation process according to the present invention. In other words, regardless of which type of FDTD method is adopted, a boundary calculation formula is set for each frequency band using the acoustic impedance Z or the boundary normal incidence sound absorption coefficient α or boundary reflection coefficient r in that frequency band, and the calculation target quantity and the wave-related quantity are calculated using the boundary calculation formula set for each frequency band and the discretization formula.

[0053] [Device configuration, wave calculation program and method] The functional configuration of the wave calculation device 1 according to an embodiment of the present invention will be described in more detail below. In the functional block diagram of Fig. 1, the wave calculation device 1 includes a communication interface unit 101, an input data storage unit 102, a calculation result storage unit 103, an input / output interface (IF) unit 104, and a processor / memory (a calculation processing system with a memory function). Here, the processor / memory stores a wave calculation program according to the present invention, has a computer function, and executes the wave calculation program to perform wave calculation processing.

[0054] For this reason, the wave calculation device 1 may be a device dedicated to wave calculation processing, but it may also be a general-purpose cloud server or non-cloud server equipped with the wave calculation program according to the present invention, and it may also be a personal computer (PC), a notebook or tablet computer, a smartphone, or even a wearable device such as a head-mounted display (HMD).

[0055] In addition, the processor and memory are functional components, (a) a frequency band (quantity to be calculated) calculation unit 111 including a sampling frequency setting unit 111a, a boundary calculation formula setting unit 111b, a discretization formula setting unit 111c, and a boundary / space interior calculation unit 111d; (b) a wave-related quantity calculation unit 112 including a sound source / sound receiving position setting unit 112a and a band-pass filter processing unit 112b; (c) Input / output control unit 121 These functional components can be regarded as functions realized by executing the wave calculation program according to the present invention stored in the processor memory. The flow of processing shown by connecting the functional components of the wave calculation device 1 in the functional block diagram of FIG. 1 with arrows can also be understood as one embodiment of the wave calculation method according to the present invention.

[0056] In the embodiment described below, the sound pressure (quantity to be calculated) of the sound wave propagating in a target room (space, sound field) surrounded by a certain room wall (boundary), and the impulse response (quantity related to wave motion) are calculated using the CE-FDTD method. However, as described above, the wave motion calculation process of the present invention is not limited to these conditions, objects, and methods.

[0057] <Calculation method for each frequency band> Also in the functional block diagram of FIG. 1, in this embodiment, the frequency band calculation unit 111 calculates the sound pressure in a frequency band in a target room for each frequency band among a plurality of frequency bands related to sound waves, using a boundary calculation formula set for that frequency band and a discretization formula based on the three-dimensional CE-FDTD method.

[0058] In this embodiment, the above-mentioned frequency bands are set to the 125 Hz band, 250 Hz band, 500 Hz band, 1000 Hz band, 2000 Hz band, 4000 Hz band, 8000 Hz band, and 16000 Hz band, which are commonly used in the field. Furthermore, in this embodiment, the frequency band calculation unit 111 acquires the value of the (normal incidence) sound absorption coefficient α of the room wall in each frequency band from the external DB 3 via the communication interface unit 101 and the input / output control unit 121 under the instruction of the external PC 2, for example. From these values, the following equation is used. (16) Z = ρc × (r p +1) / (r p -1) r p =(1-α)0.5 The acoustic impedance Z of the room wall in each frequency band is calculated using the above formula and stored in, for example, the input data storage unit 102. p is the (normal incidence) sound pressure reflection coefficient, and the above equation (16) is the formula for calculating the acoustic impedance from the (normal incidence) sound absorption coefficient α.

[0059] Incidentally, when calculating the acoustic impedance Z in each frequency band using the above formula (16), it is preferable to use the value of the acoustic impedance Z0 of the medium (e.g., air) in each frequency band. In addition, as a DB of the sound absorption coefficient α, for example, Non-Patent Document: "Sound Absorption Coefficient Data", [online], [searched on August 30, 2022], Internet<URL: http: / / www.ymec.com / products / ymcad / doc2.htm> Examples include:

[0060] Furthermore, in this embodiment, the frequency band calculation unit 111 acquires data related to the shape and size of the target room and room walls (sound field) and other parameter values ​​(e.g., Δ, Δt, ρ, c, etc.) required to implement the CE-FDTD method, for example, from the input / output IF unit 104 or the external PC 2 (via the communication interface unit 101 or the input / output control unit 121), stores them in the input data storage unit 102, and then applies them to the wave calculation process to be implemented. Note that the time step Δt as a parameter to be applied may be set by the sampling frequency setting unit 111a of the frequency band calculation unit 111 as the reciprocal of the sampling frequency fs to be used.

[0061] 1, in this embodiment, the boundary calculation formula setting unit 111b of the frequency band calculation unit 111 applies the value of the acoustic impedance Z as a boundary amount set for each frequency band to the boundary calculation formula to be used (e.g., the above formulas (4), (7), and (8)) and sets the boundary calculation formula for each frequency band. For example, a boundary calculation formula in which the acoustic impedance Z_125 for the 125 Hz band is substituted for Z in the formula (e.g., the above formulas (4), (7), and (8)) is prepared as the boundary calculation formula for the 125 Hz band.

[0062] Furthermore, the discretization formula setting unit 111c of the frequency band calculation unit 111 applies (substitutes) necessary parameter values ​​to the discretization formula to be used, and sets (prepares) the discretization formula. Furthermore, the boundary / inside space calculation unit 111d calculates, for each of the above frequency bands, the sound pressure (quantity to be calculated) in the frequency band at a position outside the room wall (boundary) adjacent to the room wall (boundary), and the sound pressure (quantity to be calculated) in the frequency band in the target room, including the calculation quantity at the room wall (boundary) calculated using the sound pressure (quantity to be calculated) at the position outside the room wall (boundary).

[0063] <Means for calculating wave-related quantities> Similarly, in the functional block diagram of Fig. 1, the wave-related quantity calculation unit 112 calculates an impulse response (wave-related quantity) in a target room from the sound pressure (quantity to be calculated) in each of the above frequency bands. Specifically, in this embodiment, the wave-related quantity calculation unit 112 (a) for each of the frequency bands, calculating an impulse response in that frequency band between the set sound source position and sound receiving position from the calculated sound pressure in that frequency band; (b) The calculated impulse responses in each frequency band are summed, or a greater or lesser weight is assigned to the impulse responses in a specific frequency band, and then the summed impulse responses are calculated to calculate the impulse response in the target room (between the set sound source position and sound receiving position).

[0064] Here, the sound source position and sound receiving position used in the above (a) can be those acquired by the sound source / sound receiving position setting unit 112a in this embodiment, for example, from the input / output IF unit 104 or the external PC 2 (via the communication interface unit 101 or the input / output control unit 121) and stored in, for example, the input data storage unit 102. Also, in the above (a), it is preferable that, at the timing when the sound pressure (distribution) in one frequency band is calculated in the frequency band calculation unit 111, the wave-related quantity calculation unit 112 calculates an impulse response in this one frequency band using the sound pressure calculation result. Of course, it is also possible to receive the sound pressure calculation results in all frequency bands and set to calculate the impulse response in all frequency bands.

[0065] Furthermore, in the above (b), it is also preferable that the impulse responses in each frequency band to be added are band-pass filtered by the band-pass filter processing unit 112b before being added. Here, this band-pass filtering can be filtering in which a passband is set according to the corresponding frequency band, as will be described in a specific example below. This makes it possible to reduce or remove components related to other frequency bands contained in an impulse response in a certain frequency band, thereby improving the accuracy of the impulse response finally calculated.

[0066] Furthermore, in the above (b), it is also preferable that the calculated impulse responses in each frequency band are weighted more or less for the impulse responses in a given frequency band (for example, after the bandpass filter process described above) and then summed up. This weighting process can also be realized by adjusting the gain for each frequency band in the bandpass filter process described above.

[0067] In any case, such weighting processing can correct errors that tend to occur in the FDTD method, for example. In fact, in the CE-FDTD method, there is a tendency for sound to be estimated slightly stronger around the sampling frequency, but by appropriately performing this weighting processing, it is possible to reduce or eliminate such errors. In addition, it is also possible to intentionally generate impulse responses that strengthen or weaken sound in a specific frequency band. This can contribute to generating sound signals that strengthen sound in a specific low frequency band in order to produce a more natural sound in content production, for example.

[0068] It is also preferable that the impulse response calculation results in the wave-related quantity calculation unit 112 and the sound pressure calculation results in the frequency band calculation unit described above are stored in the calculation result storage unit 103 and are read out from this calculation result storage unit 103 as appropriate.

[0069] Hereinafter, a specific example 1 of the calculation process of the impulse response for each frequency band described above and the calculation process of the impulse response in the target room, which is the target, will be described. Note that this specific example 1 sums up the impulse responses in each frequency band on the time axis.

[0070] (Example 1) In this specific example 1, the above-mentioned 125 Hz band, 250 Hz band, ..., 16000 Hz band are adopted as the multiple frequency bands related to the boundary conditions, and the sampling frequency fs is set to 48 kHz. The sampling frequency is usually set to twice (or more) the frequency to be used. Therefore, the sampling frequency fs (= 48 kHz) in this specific example 1 is a value that allows a margin for adopting frequency bands up to 24 kHz.

[0071] First, data on the shape and size of the target room and room walls, other parameter values ​​required to implement the CE-FDTD method, and the sound source and sound receiving positions in the target room are set. Here, the medium is air, and the room walls (boundaries) are made of rock wool sound absorbing panels. The (normal incidence) sound absorption coefficients α of these rock wool sound absorbing panels are α_125=0.60, α_250=0.65, α_500=0.66, α_1000=0.70, α_2000=0.79, α_4000=0.95, .... From this, the acoustic impedances of the boundaries in each frequency band (Z_125, Z_250, Z_500, ..., Z_16000) are determined using the above equation (16). In general, the sound absorption coefficient α of a boundary, and therefore the acoustic impedance Z of the boundary, are quantities that depend greatly on the frequency band of the sound waves.

[0072] Next, for each of the above frequency bands, the determined acoustic impedance for that frequency band (for example, Z_125 if the frequency band is 125 Hz) is applied (substituted) to the boundary calculation formula (for example, formulas (4), (7), and (8) above), and a sound pressure simulation is performed using the boundary calculation formula for that frequency band and the discretization formula (formula (2) above) using the CE-FDTD method. Specifically, an impulse sound pressure is applied to the sound source position, and the time change in sound pressure (distribution) in the target room is calculated, and the sound pressure value at the sound receiving position at each time point t is taken as the impulse response x(t) to which the boundary conditions for that frequency band are applied. As a result, the impulse response for each frequency band (x_125(t), x_250(t), x_500(t), ..., x_16000(t)) is calculated.

[0073] Next, for each of the above frequency bands, a bandpass filter b(t) corresponding to that frequency band is designed, and the impulse response x(t) in that frequency band is convoluted with this bandpass filter b(t) to calculate the impulse response y(t) after filtering. In this specific example 1, the bandpass filters (b_125(t), b_250(t), b_500(t), . . . , b_16000(t)) are used, with the upper cutoff frequency being 2 times the center frequency. 0.5(=√2) times the center frequency and the lower limit cutoff frequency is 1 / 2 0.5 In this case, for example, a bandpass filter b_125(t) for the 125 Hz band is an octave bandpass filter with a passband of 90 to 177 Hz. When designing such a bandpass filter, it is preferable to make the number of samples in all frequency bands the same in order to suppress deviations on the time axis in the summation process that is performed later.

[0074] In any case, in the above-described filtering process, for example, the impulse response y_125(t) after filtering in the 125 Hz band is expressed by the following equation, with * as the convolution integral operator: (17) y_125(t)=x_125(t)*b_125(t) This allows the impulse response z(t) in the target room to be calculated using the following formula: (18) z(t)=y_125(t)+y_500(t)+···+y_16000(t) That is, by adding up the impulse responses y(t) in each frequency band. As mentioned above, in the above formula (18), it is also possible to add up the impulse responses y(t) in each frequency band after applying a predetermined weighting.

[0075] Similarly, in the functional block diagram of FIG. 1, the wave-related quantity calculation unit 112 is configured as follows, as an embodiment different from the embodiment described above (and specific example 1): (a) For each of the above frequency bands, calculate the impulse response (wave-related quantity) in that frequency band from the sound pressure (quantity to be calculated) in that frequency band, (b) performing a Fourier transform process on the calculated impulse response in each frequency band, and determining a Fourier transformed impulse response in the target room from the impulse response in each wavenumber band that has been subjected to the Fourier transform process; (c) The impulse response in the target room is calculated by performing an inverse Fourier transform on the determined Fourier transformed impulse response. It is also preferable.

[0076] In this way, by calculating the impulse response on the frequency axis using the Fourier transform, it is possible to extract only the frequency components to be applied from the impulse response in each frequency band, and to accurately calculate the impulse response in the target room. In other words, according to this impulse response calculation process using the Fourier transform, it is possible to adjust the frequency components without using a bandpass filter, which actually leaves some amount of filtering target just before and just after the passband, and this also improves the calculation accuracy of the target impulse response.

[0077] Hereinafter, a specific example 2 of such an impulse response calculation process using a Fourier transform will be described. Note that this specific example 2 adds up the impulse responses in each frequency band on the frequency axis.

[0078] (Example 2) In this specific example 2, first, impulse responses (x_125(t), x_250(t), x_500(t), . . . , x_16000(t)) in each frequency band are calculated in the same manner as in the above specific example 1. Next, a Fourier transform is performed on each of these impulse responses to obtain impulse responses after the Fourier transform (X_125(k), X_250(k), X_500(k), . . . , x_16000(k)). Here, k is a frequency bin as a frequency point determined by the size of the Fourier transform and the sampling frequency.

[0079] Next, the Fourier-transformed impulse response Z(k) in the target room is determined from these Fourier-transformed impulse responses. Specifically, (19) Z(k)=X_125(k) (90 <k≦177) =X_250(k) (177 <k≦354) ... =X_16000(k) (11314 <k≦22627) The impulse response Z(k) is determined as shown below. Finally, the impulse response z(t) in the target room is calculated by performing an inverse Fourier transform on this impulse response Z(k).

[0080] <Other embodiments regarding sampling frequency> Next, another preferred embodiment regarding the sampling frequency in the FDTD method will be described. In the functional block diagram of Fig. 1, the sampling frequency setting unit 111a of the frequency band calculation unit 111 calculates the sound pressure (quantity to be calculated) in each of the above frequency bands by using a boundary calculation formula and a discretization formula corresponding to the sampling frequency fs set for each of the frequency bands in this embodiment.

[0081] For example, in each frequency band, the sampling frequency fs is expressed as follows: (20) fs = (upper cutoff frequency of the frequency band) × 2 (Center frequency of the frequency band) x 2 x 2 0.5 In this case, for example, the sampling frequency fs_125 in the 125 Hz band is 354 (≒125××2×2 0.5 ) Hz, which is a sufficiently large value for sampling, and there is no need to increase it to, for example, 48 kHz (set in the above specific examples 1 and 2) to match the 16,000 Hz band.

[0082] In this way, if the sampling frequency (fs_125, fs_250, fs_500, ..., fs_16000) in each frequency band is set as small as possible, the time step Δt, which is the reciprocal of the sampling frequency, is set to a larger value. In the FDTD method, the CFL number χ (=c×Δt / Δ, see formula (3) above), which represents the stability condition, is usually fixed to a predetermined value. Therefore, the grid spacing Δ in each frequency band is also set to a larger value in proportion to the time step Δt, which allows the number of grids to be smaller. As a result, the burden of wave calculations is reduced, and it is also possible to shorten the calculation time compared to the case where a uniform sampling frequency is used.

[0083] Incidentally, the above-mentioned setting of the sampling frequency for each frequency band can be implemented in, for example, the embodiment shown in Specific Example 1, the embodiment shown in Specific Example 2, and even in the case of using an FDTD method other than the CE-FDTD method. Of course, in any case, it is possible to shorten the calculation time.

[0084] 1, in this embodiment, the wave-related quantity calculation unit 112 may store the calculation result of the impulse response in the target room in, for example, the calculation result storage unit 103, and then output the result to the input / output IF unit 104 via the input / output control unit 121. In this case, the input / output IF unit 104 can display the impulse response calculation result, although it may have a display. In addition, the impulse response calculation result may be transmitted to the external PC 2 via the input / output control unit 121 and the communication interface unit 101, and used by the external PC 2.

[0085] As described above in detail, according to the present invention, in the wave-related quantity calculation process using the FDTD method, the deviation of the calculated value from the actual measured value caused by the frequency band can be suppressed, and the accuracy of the wave-related quantity calculation process can be improved. For example, it is possible to calculate the impulse response as the wave-related quantity in a predetermined sound field with higher accuracy.

[0086] Furthermore, as one application example, the present invention can contribute to estimating the impulse response of the sound collection space with greater accuracy, thereby reproducing the sound in this sound collection space more faithfully to the listener.

[0087] Furthermore, for example, it is possible to calculate more accurate impulse responses in a specific concert hall, theater, or speech venue, and create music and audio teaching materials for performances, plays, and speeches that have rich acoustic characteristics as if they were recorded in these concert halls, theaters, and speech venues. In addition, it is possible to use the music and audio teaching materials created in this way, for example, in online classes, to provide opportunities for higher quality education to children and students not only in urban areas but also in rural areas. In other words, according to the present invention, it is also possible to contribute to Goal 4 of the United Nations-led Sustainable Development Goals (SDGs) "Ensure inclusive and equitable quality education and promote lifelong learning opportunities for all."

[0088] Regarding the various embodiments of the present invention described above, various changes, modifications, and omissions within the scope of the technical idea and perspective of the present invention can be easily made by those skilled in the art. The above description is merely an example and is not intended to be restrictive in any way. The present invention is limited only by the scope of the claims and their equivalents. [Explanation of symbols]

[0089] 1 Wave calculation device 101 Communication interface section 102 Input data storage unit 103 Calculation result storage section 104 Input / Output Interface (IF) Section 111 Frequency band (quantity to be calculated) calculation section 111a Sampling frequency setting section 111b Boundary calculation formula setting part 111c Discretization formula setting part 111d Boundary / space calculation part 112 Wave-related quantity calculation unit 112a Sound source / sound receiving position setting section 112b Bandpass filter processing section 121 Input / Output Control Unit 2 External PC 3 External Database (DB)

Claims

1. A program for calculating a wave-related quantity, which is a quantity related to a wave in a space having a predetermined boundary, using the finite-difference time-domain method, A frequency band calculation means for calculating a calculation target quantity related to the calculation of the wave in the frequency band in the space by using a boundary calculation formula set for each frequency band among a plurality of frequency bands included in the frequency range related to the wave and a discretization formula related to the wave; a wave-related quantity calculation means for calculating the wave-related quantity in the space from the calculation target quantity in the frequency band in the space calculated for each frequency band; A wave calculation program characterized by causing a computer to function as described above.

2. The wave calculation program according to claim 1, characterized in that the frequency band calculation means applies, for each frequency band, the value of the boundary quantity set for the frequency band to the boundary calculation formula including a boundary quantity, which is a quantity related to the absorption or impedance of the boundary, as a parameter, and calculates, using the boundary calculation formula for the frequency band to which the value of the boundary quantity has been applied and the discretization formula, the calculation target quantity in the frequency band at a position outside the boundary adjacent to the boundary, and the calculation target quantity at the boundary calculated using the calculation target quantity at the position outside the boundary, in the frequency band in the space.

3. The wave calculation program according to claim 1 or 2, characterized in that the frequency band calculation means calculates the calculation target quantity in the frequency band in the space using the boundary calculation formula and the discretization formula corresponding to a sampling frequency set for each frequency band.

4. 3. The wave calculation program according to claim 1 or 2, wherein the wave-related quantity calculation means calculates, for each frequency band, the wave-related quantity in that frequency band from the calculation target quantity in that frequency band, and calculates the wave-related quantity in the space by adding up the wave-related quantities calculated in each frequency band, or by adding up the wave-related quantities in a predetermined frequency band after assigning a larger or smaller weight to the wave-related quantity in that frequency band.

5. 3. The wave-related quantity calculation means calculates, for each frequency band, the wave-related quantity in the frequency band from the calculation target quantity in the frequency band, applies band-pass filtering to the calculated wave-related quantity in each frequency band, with a passband set to match the frequency band, and calculates the wave-related quantity in the space by adding up the wave-related quantities in each frequency band that have been subjected to band-pass filtering, or by adding up the wave-related quantities in a predetermined frequency band after applying a larger or smaller weight to the wave-related quantities.

6. The wave-related quantity calculation means calculates, for each frequency band, the wave-related quantity in the frequency band from the calculation target quantity in the frequency band, performs Fourier transform processing on the calculated wave-related quantity in each frequency band, determines a Fourier-transformed wave-related quantity in the space from the wave-related quantity in each wave number band subjected to the Fourier transform processing, and performs inverse Fourier transform processing on the determined Fourier-transformed wave-related quantity to calculate the wave-related quantity in the space.

7. The wave calculation program according to claim 1 or 2, characterized in that the finite-difference time-domain method is a Compact-Explicit Finite-Difference Time-Domain (CE-FDTD) method, the boundary calculation formula is set for each reference direction in a plurality of reference directions, and the discretization formula is a formula corresponding to the plurality of reference directions.

8. 3. The wave calculation program according to claim 1 or 2, characterized in that the finite difference time domain method is a method in which the target variable in the governing equation includes a vector quantity, or a method in which the target variable in the governing equation is a scalar quantity.

9. 3. The wave calculation program according to claim 1 or 2, characterized in that the wave is a sound wave, the calculation target quantity is a quantity related to sound pressure, and the wave-related quantity is an impulse response between a set sound source position and a sound receiving position.

10. An apparatus for calculating a wave-related quantity, which is a quantity related to a wave in a space having a predetermined boundary, using a finite-difference time-domain method, comprising: A frequency band calculation means for calculating a calculation target quantity related to the calculation of the wave in the frequency band in the space by using a boundary calculation formula set for each frequency band among a plurality of frequency bands included in the frequency range related to the wave and a discretization formula related to the wave; a wave-related quantity calculation means for calculating the wave-related quantity in the space from the calculation target quantity in the frequency band in the space calculated for each frequency band; A wave calculation device comprising:

11. A method for calculating a wave-related quantity, which is a quantity related to a wave in a space having a predetermined boundary, using a finite-difference time-domain method, comprising the steps of: A step of calculating a calculation target quantity related to the calculation of the wave in the frequency band in the space, for each frequency band among a plurality of frequency bands included in the frequency range related to the wave, using a boundary calculation formula set for the frequency band and a discretization formula related to the wave; calculating the wave-related quantity in the space from the calculation target quantity in the frequency band in the space calculated for each frequency band; A computer-implemented wave calculation method comprising:

Images