Sound field reproduction device and program

Abstract

To reproduce a sound field formed by a moving directional sound source when a drive signal for a speaker array is generated by SDM.SOLUTION: A drive signal calculator 10 of a sound field reproducing device 1 calculates a reproduced sound field angle spectrum using a preset reproduction boundary yref, a calculated wave number kx in an x-axis direction, a wave number k, and an angular frequency ω, calculates a desired sound field angle spectrum by using a preset buffer length N, spherical harmonic spectra Glmi, a truncation order Nu and a desired boundary ydes, the calculated wave number kx in the x-axis direction, the wave number k and the angular frequency ω, and virtual sound source directions (α(n), β(n), γ(n)), virtual sound source coordinates (rs(n), θs(n ), φs(n)) and a sound source signal s(n) for each time sample n in a buffering section of the buffer length N, and divides the desired sound field angle spectrum by the reproduced sound field angle spectrum, thereby calculating the angle spectrum of the drive signal.SELECTED DRAWING: Figure 3

Description

[Technical field]

[0001] The present invention relates to a sound field reproducing device and a program for reproducing a sound field formed by a directional, moving sound source using a speaker array. [Background technology]

[0002] Research has been conducted on sound field reproduction technology that uses multiple speakers to create an arbitrary sound field in a certain space. Known technologies for sound field reproduction include Wave Field Synthesis (WFS), Boundary Surface Control (BoSC), Higher Order Ambisonics (HOA), and Spectral Division Method (SDM).

[0003] 7 is a diagram showing an example of speaker array arrangement in wave field synthesis (WFS). WFS is a method based on the Rayleigh integral, and uses a speaker array 100-1 arranged linearly or planarly to reproduce the sound pressure or sound pressure gradient on the boundary plane of a desired sound field, thereby forming a wavefront due to a sound source coming from outside the boundary.

[0004] Fig. 8 shows an example of speaker array placement in boundary sound field control (BoSC). BoSC is a method based on the Kirchhoff-Helmholtz integral equation. BoSC places speaker array 100-2 outside the area where the sound field is to be reproduced so as to surround the listener toward the inside of the area, and reproduces the desired sound field inside the area by controlling the sound pressure and sound pressure gradient on the boundary surface of the area using an inverse system.

[0005] Fig. 9 is a diagram showing an example of speaker array arrangement in Higher Order Ambisonics (HOA). HOA is a method of expressing a desired sound field using spherical harmonic functions, and reproducing the sound field by matching the spherical harmonic coefficients of the reproduced sound field with the spherical harmonic coefficients of the desired sound field using a spherical speaker array 100-3 arranged to surround the listener toward the inside of the sphere.

[0006] Fig. 10 is a diagram showing an example of speaker array arrangement in the spectral division method (SDM). SDM is a method for expressing a desired sound field using an angular spectrum, and reproducing the sound field by matching the angular spectrum of the reproduced sound field with the angular spectrum of the desired sound field using a speaker array 100-4 arranged in a line or a plane. For details of SDM, see, for example, Non-Patent Document 1.

[0007] [Calculation method of drive signal in SDM] A method for calculating a drive signal in a typical SDM will be described below. Fig. 11 is a diagram for explaining a typical SDM.

[0008] y=y in xyz space 0 ,z=z 0 The coordinates r0 = (x 0 , y 0 ,z 0 At the point r = (x, y), the driving signal with angular frequency ω is D(r0,ω). ref ,z 0 The sound pressure P(r,ω) at the point r is expressed by the following equation:

number

[0009] Here, G(r-r0,ω) is the transfer function of direction and distance expressed by the vector (r-r0). When a free sound field is assumed as the reproduced sound field, it is expressed by the three-dimensional free-field Green's function and is given by the following formula.

number

[0010] The sound pressure P(r,ω) can be interpreted as being obtained by convolving the drive signal D(r0,ω) and the transfer function G(r-r0,ω) in space, and y=y ref By performing a Fourier transform along the x-axis direction in x ,y ref ,z 0 ,ω) is obtained.

number

[0011] P^(k x ,y ref ,z 0 ,ω) is the sound pressure P(r,ω) at y=y ref The angular spectrum is obtained by performing a spatial Fourier transform on D^(k x ,y 0 ,z 0 ,ω) converts the driving signal D(r0,ω) into y=y 0 The angular spectrum is obtained by performing a spatial Fourier transform on G^(k x ,y ref -y 0 ,z 0 ,ω) is the transfer function G(r-r0,ω) as y=y ref The angular spectrum is obtained by performing a spatial Fourier transform with k x is the wave number along the x-axis.

[0012] Here, the desired sound field to be reproduced is y=y des The angular spectrum obtained by spatial Fourier transform of the sound pressure distribution at d (k x ,y des ,z 0 ,ω), the angular spectrum D^(k x ,y 0 ,z 0 ,ω) is expressed by the following formula:

number

[0013] Coordinate r s =(x s ,y s ,z 0 ) driven by S(ω) at frequency ω, the desired sound field y=y des Angular spectrum P^ at d (k x ,y des ,z 0 , ω), i.e., the desired sound field angular spectrum, is expressed by the following equation:

number

[0014] Here, H 0 (2) is the Hankel function of the second kind, and K 0 is the 0th order modified Bessel function. Also, 0≦k 2 -k x 2 In the case of 0>k, the sound wave propagates far away. 2 -k x 2 In this case, the sound wave is called an evanescent wave, and the amplitude of the wave decays rapidly and the wave does not propagate far.

[0015] Similarly, the angular spectrum G^(k x ,y ref -y 0 ,z 0 , ω), i.e., the reproduced sound field angular spectrum, is expressed by the following equation:

number

[0016] Therefore, the angular spectrum of the drive signal D^(k x ,y 0 ,z 0 ,ω) is expressed by the following formula.

number

[0017] 12 is a diagram for explaining the reproduction area and the reproduction boundary. When the speaker array 100 arranged in a line is an infinite linear sound source, the infinite linear sound source is driven using the drive signal of the above-mentioned formula (7), and thereby, y≧y ref Therefore, in this case, the sound field can be reproduced in the area of ​​y=y ref The reproduction boundary, y ≧ y ref This area is called the clipping region.

[0018] The angular spectrum D^(k x ,y 0 ,z 0 , ω) is assumed to perform wave field synthesis using an infinite linear sound source, but in reality it is realized using a speaker array 100 consisting of a finite number of speaker units arranged in a line.

[0019] Therefore, the angular spectrum D^(k x ,y 0 ,z 0, ω), it is necessary to discretize the frequency and angular spectrum, truncate it to a finite length, and then perform an Inverse Discrete Fourier Transform (IDFT) to calculate the drive signal in the time-frequency domain for each speaker unit.

[0020] The number of speaker units constituting the speaker array 100 is M (a natural number), the interval between the speaker units is Δx, and the length of the speaker array is L (=(M-1)Δx). x ,y 0 ,z 0 By performing an inverse discrete Fourier transform on x, ω), the driving signal D~(x,ω) in the time-frequency domain for the speaker unit placed at x is expressed by the following equation.

number

[0021] Then, for each speaker unit that constitutes the speaker array 100, a driving signal in the time domain is calculated.

[0022] By using such sound field reproduction technology as SDM, it is possible to reproduce a sound field formed by a point sound source located at coordinates specified by the user. However, while conventional speaker driving signal calculation methods can reproduce a sound field formed by a stationary sound source, they do not take into account a situation in which the sound source has directionality and moves continuously.

[0023] In order to solve this problem, Non-Patent Document 2 proposes a method for reproducing a sound field formed by a moving sound source, based on SDM. [Prior art documents] [Non-patent literature]

[0024] [Non-Patent Document 1] J. Ahrens, and S. Spors, “Sound field reproduction using planar and linear arrays of loudspeakers”, IEEE Trans. Audio, Speech, Lang. Process., vol.18, no.8, pp.2038-2050, Nov.2010. [Non-Patent Document 2] Firtha Gergely and Fiala Peter, “Sound Field Synthesis of Uniformly Moving Virtual Monopoles”, JAES Volume 63 Issue 1 / 2 pp. 46-53, January 2015 Summary of the Invention [Problem to be solved by the invention]

[0025] However, since the situation assumed here is only when the sound source is omnidirectional and moves linearly, the above-mentioned sound field reproduction method cannot be directly applied to a directional sound source.

[0026] Therefore, the present invention has been made to solve the above-mentioned problems, and its object is to provide a sound field reproduction device and program capable of reproducing a sound field formed by a moving directional sound source when generating a drive signal for a speaker array using SDM. [Means for solving the problem]

[0027] In order to solve the above problem, the sound field reproducing device of claim 1 is a sound field reproducing device that generates a drive signal for reproducing a sound field formed by a moving directional sound source using a speaker array consisting of a plurality of speaker units, the speaker array being located at z=z in an xyz space. 0 y=y 0 It is assumed that the y ref The reproduced boundary, f s is the sampling frequency, M is the number of speaker units, Δx is the speaker unit interval, N is the buffer length, F is the number of DFT points, Nu is the truncation order, y des the desired boundary, G l m i, the sampling frequency f s and the l (l=0 to L)-th order spherical harmonic spectrum (i is the frequency bin number, l is the order, and m is the rank) representing the radiation characteristic of the directional sound source for each frequency bin determined by the DFT point number F, and k x Let be the wave number in the x-axis direction, k be the wave number, ω be the angular frequency, c be the sound speed, n be the time sample, (α(n), β(n), γ(n)) be the virtual sound source direction which is the direction of the directional sound source for each time sample n shown in ZYZ Euler angles, and (r s (n),θ s (n),φ s (n)) is defined as the virtual sound source coordinates, which are the coordinates of the directional sound source for each time sample n, (r s (n) is the distance from the origin to the virtual sound source, θ s (n) is the zenith angle, φ s (n) indicates an azimuth angle. With s(n) being a sound source signal that is a signal of the directional sound source for each time sample n, the number M of speaker units, the speaker unit interval Δx, and an index m=0 to M−1 (m is an integer) that are set in advance, the formula: k x = 2πm / (ΔxM), the wave number k in the x-axis direction x and a wave number calculation unit in the x-axis direction that calculates the sampling frequency f s From the number of DFT points F and the frequency index l = 0 to F-1 (l is an integer), the formula: ω = 2πf sa wave number calculation unit that calculates the angular frequency ω from l / F and calculates the wave number k from the equation: k=ω / c; 0 (2) is a second kind of Hankel function with a degree of 0, i is an imaginary unit, and the wave number k in the x-axis direction calculated by the x-axis direction wave number calculation unit is x , the wave number k calculated by the wave number calculation unit, and the preset reproduction boundary y ref and the value y 0 From the formula: a reproduction sound field angle spectrum calculation unit that calculates a reproduction sound field angle spectrum according to TIFF0007674954000009.tif17170, a desired sound field angle spectrum calculation unit that calculates a desired sound field angle spectrum, and a division unit that obtains an angular spectrum of the drive signal by dividing the desired sound field angle spectrum calculated by the desired sound field angle spectrum calculation unit by the reproduced sound field angle spectrum calculated by the reproduced sound field angle spectrum calculation unit, l m i, the truncation order Nu and the desired boundary y des , the wave number k in the x-axis direction calculated by the x-axis direction wave number calculation unit x , the wave number k and the angular frequency ω calculated by the wave number calculation unit, and the virtual sound source direction (α(n), β(n), γ(n)) and the virtual sound source coordinates (r s (n),θ s (n),φ s (n)) and the source signal s(n), the formula: TIFF0007674954000010.tif22170(P^ d (k x ,y des , ω) are the desired sound field angular spectrum, ν and μ are the degree and order of a predetermined spherical harmonic function, respectively, i is the imaginary unit, H μ (2) μ is the second kind of Hankel function, D μμ’ ν is the Wigner-D function, Let TIFF0007674954000011.tif18170 be a given spherical harmonic spectrum, and the sigma Σ operation of order ν be truncated by the truncation order Nu. formula: TIFF0007674954000012.tif39170 (l,m are the degree and order of the given spherical harmonic function, respectively, TIFF0007674954000013.tif17170, the spherical harmonic spectrum G l m i, the spherical harmonic spectrum, D mm’ l is the Wigner D function, S lν mμ ,S^ lν mμ Let R be a given function. s (τ) is the virtual sound source coordinate (r s (n),θ s (n),φ s (n)), s(τ) is the sound source signal s(n), j ν is the νth order spherical Bessel function, h ν (2) Let ν be the two-kind spherical Hankel function, Y ν μ is a μ-th order spherical harmonic function, R is the coordinate of the sound receiving point, and r is the distance from the origin to the sound receiving point. The integral of dτ and the sigma Σ operation of the order l are truncated by the truncation order Nu. formula: TIFF0007674954000014.tif23170(h q (2) q is a two-kind spherical Hankel function, Y q μ-m ￣ is the complex conjugate of the qth μ-mth order spherical harmonic function, g(l,m;ν,-μ;q) is the gaunt coefficient, and b=R s (τ).) formula: TIFF0007674954000015.tif24170(j q is the qth order spherical Bessel function, Y q μ-m￣ is the complex conjugate of the q-th order μ-m spherical harmonic function, g(l,m;ν,-μ;q) is the Gaunt coefficient, and b=R s (τ). The desired sound field angular spectrum is calculated by:

[0028] The program of claim 2 further comprises a program for controlling a computer constituting a sound field reproducing device that generates a drive signal for reproducing a sound field formed by a moving directional sound source using a speaker array consisting of a plurality of speaker units, the program causing the speaker array to move in a direction z=z in the xyz space. 0 y=y 0 It is assumed that it is placed on the y ref The reproduced boundary, f s is the sampling frequency, M is the number of speaker units, Δx is the speaker unit interval, N is the buffer length, F is the number of DFT points, Nu is the truncation order, y des the desired boundary, G l m i, the sampling frequency f s and the l (l=0 to L)-th order spherical harmonic spectrum (i is the frequency bin number, l is the order, and m is the rank) representing the radiation characteristic of the directional sound source for each frequency bin determined by the DFT point number F, and k x Let be the wave number in the x-axis direction, k be the wave number, ω be the angular frequency, c be the sound speed, n be the time sample, (α(n), β(n), γ(n)) be the virtual sound source direction which is the direction of the directional sound source for each time sample n shown in ZYZ Euler angles, and (r s (n),θ s (n),φ s (n)) is defined as the virtual sound source coordinates, which are the coordinates of the directional sound source for each time sample n, (r s (n) is the distance from the origin to the virtual sound source, θ s (n) is the zenith angle, φ s (n) indicates an azimuth angle. With s(n) being a sound source signal that is a signal of the directional sound source for each time sample n, the number M of speaker units, the speaker unit interval Δx, and an index m=0 to M−1 (m is an integer) that are set in advance, the formula: k x = 2πm / (ΔxM), the wave number k in the x-axis direction xa wave number calculation unit in the x-axis direction that calculates the sampling frequency f s From the number of DFT points F and the frequency index l = 0 to F-1 (l is an integer), the formula: ω = 2πf s a wave number calculation unit that calculates the angular frequency ω from l / F and calculates the wave number k from the equation: k=ω / c; 0 (2) is a second kind of Hankel function with a degree of 0, i is an imaginary unit, and the wave number k in the x-axis direction calculated by the x-axis direction wave number calculation unit is x , the wave number k calculated by the wave number calculation unit, and the preset reproduction boundary y ref and the value y 0 From the formula: According to TIFF0007674954000016.tif17170, a reproduced sound field angle spectrum calculation unit calculates a reproduced sound field angle spectrum, a desired sound field angle spectrum calculation unit calculates a desired sound field angle spectrum, and a division unit that divides the desired sound field angle spectrum calculated by the desired sound field angle spectrum calculation unit by the reproduced sound field angle spectrum calculated by the reproduced sound field angle spectrum calculation unit to obtain the angular spectrum of the drive signal, and the desired sound field angle spectrum calculation unit calculates the spherical harmonic spectrum G l m i, the truncation order Nu and the desired boundary y des , the wave number k in the x-axis direction calculated by the x-axis direction wave number calculation unit x , the wave number k and the angular frequency ω calculated by the wave number calculation unit, and the virtual sound source direction (α(n), β(n), γ(n)) and the virtual sound source coordinates (r s (n),θ s (n),φ s (n)) and the source signal s(n), the formula: TIFF0007674954000017.tif22170(P^ d (k x ,y des, ω) are the desired sound field angular spectrum, ν and μ are the degree and order of a predetermined spherical harmonic function, respectively, i is the imaginary unit, H μ (2) μ is the second kind of Hankel function, D μμ’ ν is the Wigner-D function, Let TIFF0007674954000018.tif18170 be a given spherical harmonic spectrum, and the sigma Σ operation of order ν be truncated by the truncation order Nu. formula: TIFF0007674954000019.tif37170 (l,m are the degree and order of a given spherical harmonic function, respectively, TIFF0007674954000020.tif15170, the spherical harmonic spectrum G l m i, the spherical harmonic spectrum, D mm’ l is the Wigner D function, S lν mμ ,S^ lν mμ Let R be a given function. s (τ) is the virtual sound source coordinate (r s (n),θ s (n),φ s (n)), s(τ) is the sound source signal s(n), j ν is the νth order spherical Bessel function, h ν (2) Let ν be the two-kind spherical Hankel function, Y ν μ is a μ-th order spherical harmonic function, R is the coordinate of the sound receiving point, and r is the distance from the origin to the sound receiving point. The integral of dτ and the sigma Σ operation of the order l are truncated by the truncation order Nu. formula: TIFF0007674954000021.tif22170(h q (2) q is a two-kind spherical Hankel function, Y q μ-m￣ is the complex conjugate of the qth μ-mth order spherical harmonic function, g(l,m;ν,-μ;q) is the gaunt coefficient, and b=R s (τ).) formula: TIFF0007674954000022.tif23170(j q is the qth order spherical Bessel function, Y q μ-m ￣ is the complex conjugate of the q-th order μ-m spherical harmonic function, g(l,m;ν,-μ;q) is the Gaunt coefficient, and b=R s (τ). The desired sound field angular spectrum is calculated by: Effect of the Invention

[0029] As described above, according to the present invention, when generating driving signals for a speaker array using SDM, it is possible to reproduce a sound field formed by a moving directional sound source. [Brief description of the drawings]

[0030] [Figure 1] FIG. 2 is a diagram illustrating an example of a configuration of a sound field assumed in an embodiment of the present invention. [Diagram 2] 1 is a block diagram showing an example of the configuration of a sound field reproduction device according to an embodiment of the present invention. [Diagram 3] 4 is a block diagram showing an example of the configuration of a drive signal calculation unit; FIG. [Figure 4] 10 is a flowchart showing an example of processing by a drive signal calculation unit. [Diagram 5] 11 is a block diagram showing an example of the configuration of a reproduced sound field angular spectrum calculation unit. FIG. [Figure 6] 4 is a block diagram showing an example of the configuration of a desired sound field angular spectrum calculation unit. FIG. [Figure 7] FIG. 1 is a diagram showing an example of speaker array arrangement in wave field synthesis (WFS). [Figure 8] FIG. 1 is a diagram showing an example of speaker array arrangement in boundary sound field control (BoSC). [Figure 9]FIG. 1 is a diagram showing an example of speaker array arrangement in Higher Order Ambisonics (HOA). [Figure 10] FIG. 1 is a diagram showing an example of speaker array arrangement in the spectral division method (SDM). [Figure 11] FIG. 1 is a diagram for explaining a general SDM. [Figure 12] FIG. 2 is a diagram illustrating a reproduction area and a reproduction boundary. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0031] The present invention is characterized in that, when generating a driving signal for a speaker array by SDM, a sound field formed by a directional sound source moving along an arbitrary trajectory is formulated.

[0032] This makes it possible to reproduce the sound field formed by a smoothly moving directional sound source, including the Doppler effect.

[0033] [Speaker array drive signal] First, a speaker array driving signal generated by the sound field reproduction device according to the embodiment of the present invention will be described.

[0034] Consider the convolution of a moving directional sound source with respect to the time waveform s(t) of the sound source signal. As the sound source moves, the impulse response between the sound source and the sound receiving point changes over time, so this convolution can be said to be the convolution of the time-varying impulse response and the sound source.

[0035] FIG. 1 is a diagram for explaining an example of the configuration of a sound field assumed in an embodiment of the present invention. The coordinates of a sound source at time t in a spherical coordinate system are R s (t)=(r s (t),θ s (t),φ s (t)), and the impulse response from the sound source to the sound receiving point coordinate R = (r, θ, φ) is g(RR s Let (t),t). s(t) is the distance from the origin to the coordinates of the sound source at time t, and r is the distance from the origin to the sound receiving point at time t. θ s (t),θ is the zenith angle,φ s (t),φ is the azimuth angle.

[0036] The response at the coordinates R(x, y) of the sound receiving point can be expressed as the convolution of the time waveform s(t) of the sound source signal and the time-varying impulse response, as shown in the following equation.

number

[0037] And the impulse response g(RR s (τ), t-τ) is expressed by the following equation by inverse Fourier transform on the time-frequency axis.

number

[0038] Furthermore, the transfer function G(RR s (τ),ω') are R s By expanding the spherical harmonic spectrum with (τ) as the expansion center, the spherical harmonic spectrum G^ from the sound source to the sound receiving point is obtained as follows: l m It can be expressed using (τ,ω'). The spherical harmonic spectrum G^ l m The l and m in (τ, ω') are the degree and the order of the spherical harmonic spectrum.

number

[0039] Here, h l (2) is the l-order spherical Hankel function, where l is the order of the spherical Hankel function. l m(·) is the lth order spherical harmonic function of degree m, where l and m are the degree and order of the spherical harmonic function, respectively. k is the wave number.

[0040] Furthermore, the following equation holds true from the addition theorem of spherical Bessel functions.

number

[0041] where j ν (·) is the ν-th order spherical Bessel function, where ν is the order of the spherical Bessel function. Y ν μ (·) is the μ-th order spherical harmonic function of degree ν, where ν and μ are the degree and order of the spherical harmonic function, respectively. ν (2) (·) is the ν-order spherical Hankel function, where ν is the order of the spherical Hankel function.

[0042] Function S lν mμ (·),S^ lν mμ (·) is the transfer function G(RR s This function is derived by applying the additive theorem of spherical Bessel functions to (τ), ω') and re-expanding it into spherical harmonics, and is expressed by the following formula. The function S lν mμ (·),S^ lν mμ In (·), l and ν are the orders before and after the redeployment, respectively, and m and μ are the orders before and after the redeployment, respectively.

number

number

[0043] where i is the imaginary unit. h q (2) is a q-ordered second-order spherical Hankel function, where q is the order of the spherical Hankel function. Y q μ-m￣(·) is the complex conjugate of the qth order μ-m spherical harmonic function, where q, μ-m are the degree and order of the spherical harmonic function, respectively. g(l,m;ν,-μ;q) is the gaunt coefficient. j q (·) is the qth order spherical Bessel function, where q is the order of the spherical Bessel function.

[0044] For details of the addition theorem of spherical Bessel functions, please refer to the following non-patent literature. [Non-patent literature] PA Martin, Multiple Scattering: Interaction of Time-Harmonic Waves with N obstacles, Cambridge university press, 2006.

[0045] Also, the spherical harmonic spectrum G^ l m Let (τ,ω') be the Wigner-D function D mm’ l Using (α(τ), β(τ), γ(τ)), this can be expressed as follows:

number

[0046] Here, Wigner's D function D mm’ l In (α(τ), β(τ), γ(τ)), α(τ), β(τ), γ(τ) are ZYZ Euler angles, l is the order of the spherical harmonic spectrum, and m, m' are the orders of the spherical harmonic spectrum. The above formula (15) is the spherical harmonic spectrum that shows the radiation characteristics. Represents the rotation of TIFF0007674954000030.tif13170.

[0047] For details about the Wigner-D function, please refer to the following non-patent literature. [Non-patent literature] B. Rafaely, Fundamentals of Spherical Array Processing, Springer, 2015.

[0048] The sound pressure p(R,t) at time t of the sound receiving point coordinates R=(r,θ,φ) can be transformed into the following equation by Fourier transforming it in the time direction using the equations (12) and (15).

number

[0049] From this, the spherical harmonic spectrum is expanded with the origin of the coordinate system of the sound receiving point coordinate R = (r, θ, φ) and the sound pressure p(R, t) at time t as the center. TIFF0007674954000032.tif16170 can be expressed as follows: where v, μ are the order and degree of this spherical harmonic spectrum.

number

[0050] Here, the spherical harmonic spectrum TIFF0007674954000034.tif16170 is included in the above formula (15). Corresponds to TIFF0007674954000035.tif12170. D mm’ l (·) is the Wigner D function. The function S lν mμ (·),S^ lν mμ (·) correspond to the above equations (13) and (14), respectively. s(τ) is the sound source signal at time τ. j ν (·) is the νth order spherical Bessel function, Y ν μ (·) is the μ-th order spherical harmonic function, and h ν (2) (·) is the ν-dependent second-order spherical Hankel function.

[0051] The following non-patent literature proposes a method of converting any spherical harmonic spectrum into an angular spectrum, applying it to SDM, and reproducing a sound field using a line array speaker. [Non-patent literature] Takuma Okamoto, “Angular spectrum decomposition based 2.5D higher order spherical harmonic sound field synthesis with a linear loudspeaker arrays”, Proc WASPAA 2017, pp.180-184, 2017

[0052] According to this non-patent document, the spherical harmonic spectrum is TIFF0007674954000036.tif16170 to angular spectrum P^ d (k x ,y des ,z 0 ,ω), where ν,μ are the degree and order of this spherical harmonic spectrum.

number

[0053] where ν and μ are spherical harmonic functions Y ν μ is the degree and order of (·), i is the imaginary unit, and H μ (2) (·) is the μ-dependent Hankel function of the second kind. D μμ’ ν (·) is the Wigner D function. In the above formula (18), z 0 is omitted, and P^ d (k x ,y des ,ω)=P^ d (k x ,y des ,z 0 ,ω).

[0054] In the above formula (18), the desired boundary y des(>0) is the distance to a line parallel to the x-coordinate of the desired sound field, where the sound pressure distribution is to be reproduced in the reproduced sound field. Therefore, by applying the above formulas (18) and (6) to the above formula (4), it is possible to reproduce a sound field formed by a directional sound source that moves and rotates in an arbitrary trajectory using a line array speaker.

[0055] However, when actually implementing the formula (17) included in the formula (18), the integral for τ needs to be truncated to a finite length, and the data actually inputted is a discrete value, so it needs to be appropriately discretized before calculation. Also, when implementing the formula (17) and the formula (18), the orders l and ν of the spherical harmonic spectrum need to be truncated to a finite order.

[0056] In this way, the sound field reproduction device according to the embodiment of the present invention uses the above equation (18) to calculate the desired sound field angular spectrum P^ d (k x ,y des ,ω)=P^ d (k x ,y des ,z 0 ,ω) and use the formula (6) to obtain the reproduced sound field angular spectrum G^(k x ,y ref -y 0 ,z 0 , ω). Then, the sound field reproduction device calculates the desired sound field angular spectrum P^ d (k x ,y des ,z 0 ,ω) and the reproduced sound field angular spectrum G^(k x ,y ref -y 0 ,z 0 ,ω) to obtain the angular spectrum of the driving signal D^(k x ,y 0 ,z 0 ,ω) is calculated.

[0057] The sound field reproduction device uses the angular spectrum D^(k x ,y 0 ,z 0, ω) is subjected to a discrete inverse Fourier transform to calculate a driving signal D̂(x, ω) in the time-frequency domain for each speaker unit constituting the speaker array 100.

[0058] The sound field reproduction device calculates a speaker drive signal in the time domain for each speaker unit by performing a discrete inverse Fourier transform on the drive signal D~(x,ω) in the time frequency domain, thereby generating a speaker drive signal for reproducing a sound field formed by a moving directional sound source.

[0059] [Sound field reproduction device] Next, a sound field reproduction device according to an embodiment of the present invention will be described. Fig. 2 is a block diagram showing an example of the configuration of a sound field reproduction device according to an embodiment of the present invention. This sound field reproduction device 1 includes a drive signal calculation unit 10, a spatial frequency domain discrete inverse Fourier transform unit 11, and a time frequency domain discrete inverse Fourier transform unit 12.

[0060] The drive signal calculation unit 10 calculates the reproduction boundary y ref , sampling frequency f s , number of speaker units M, speaker unit interval Δx, buffer length N, number of DFT points F, sampling frequency f s The l (l=0~L)th order spherical harmonic spectrum G represents the radiation characteristics of the sound source for each frequency bin determined by the DFT point number F, which is the number of DFT (or FFT) taps. l m i (i is the frequency bin number, l is the order, m is the rank), the truncation order Nu, and the desired boundary y des Enter.

[0061] In addition, the excitation signal calculation unit 10 calculates the virtual sound source direction (α(n), β(n), γ(n)) and virtual sound source coordinates (r s (n),θ s (n),φ s The virtual sound source direction (α(n), β(n), γ(n)) is expressed by Euler angles, and the Wigner D function D mm’l It corresponds to (α(τ),β(τ),γ(τ)) in (α(τ),β(τ),γ(τ)). s (n),θ s (n),φ s (n)) is the coordinate R of the sound source shown in Figure 1. s (t)=(r s (t),θ s (t),φ s (t)).

[0062] The virtual sound source (point sound source) moves and rotates along an arbitrary trajectory, and its direction and coordinates are the virtual sound source direction (α(n), β(n), γ(n)) and virtual sound source coordinates (r s (n),θ s (n),φ s (n)).

[0063] The drive signal calculation unit 10 calculates the drive signal for the SDM based on the above formula (17) and formula (18). Specifically, the drive signal calculation unit 10 calculates the desired sound field angular spectrum P^ in the above formula (18). d (k x ,y des ,z 0 ,ω) and the reproduced sound field angular spectrum G^(k x ,y ref -y 0 ,z 0 , ω). Then, the drive signal calculation unit 10 calculates the desired sound field angular spectrum P^ d (k x ,y des ,z 0 ,ω) and the reproduced sound field angular spectrum G^(k x ,y ref -y 0 ,z 0 ,ω), the wave number k in the x-axis direction x For each speaker unit (number of speaker units M), the angular spectrum D^(k x ,y 0 ,z 0 ,ω) is calculated.

[0064] The drive signal calculation unit 10 calculates the wave number k in the x-axis direction. x Angular spectrum of the drive signal for each D^(k x ,y 0 ,z 0 , ω) to the spatial frequency domain discrete inverse Fourier transform unit 11. The drive signal calculation unit 10 will be described in detail later.

[0065] The spatial frequency domain discrete inverse Fourier transform unit 11 receives the number M of speaker units and the speaker unit interval Δx that are set in advance, and also receives the wave number k in the x-axis direction from the drive signal calculation unit 10. x Angular spectrum of the drive signal for each D^(k x ,y 0 ,z 0 ,ω).

[0066] The spatial frequency domain discrete inverse Fourier transform unit 11 calculates the angular spectrum D^(k x ,y 0 ,z 0 , ω) to obtain the driving signals D~(x,ω) in the time frequency domain for each speaker unit constituting the speaker array 100. Specifically, the spatial frequency domain discrete inverse Fourier transform unit 11 calculates the speaker array length L=(M-1)Δx from the number of speaker units M and the speaker unit interval Δx, and performs the calculation of the above formula (8).

[0067] The spatial frequency domain discrete inverse Fourier transform unit 11 outputs the time frequency domain drive signal D~(x, ω) for each speaker unit to the time frequency domain discrete inverse Fourier transform unit 12.

[0068] The time-frequency domain discrete inverse Fourier transform unit 12 inputs the time-frequency domain drive signals D~(x,ω) for each speaker unit from the spatial frequency domain discrete inverse Fourier transform unit 11. Then, the time-frequency domain discrete inverse Fourier transform unit 12 performs a discrete inverse Fourier transform on the time-frequency domain drive signals D~(x,ω) to obtain time-domain speaker drive signals for each speaker unit constituting the speaker array 100.

[0069] The time-frequency domain discrete inverse Fourier transform unit 12 outputs the time domain speaker drive signals for each speaker unit to the corresponding speaker units constituting the speaker array 100 .

[0070] (Drive signal calculation unit 10) Next, a detailed description will be given of the drive signal calculation section 10 shown in Fig. 2. Fig. 3 is a block diagram showing an example of the configuration of the drive signal calculation section 10, and Fig. 4 is a flowchart showing an example of the processing performed by the drive signal calculation section 10.

[0071] The drive signal calculation unit 10 includes an x-axis direction wave number calculation unit 20, a wave number calculation unit 21, a reproduced sound field angular spectrum calculation unit 22, a desired sound field angular spectrum calculation unit 23, and a division unit 24. These components will be described with reference to the flowchart of FIG.

[0072] The drive signal calculation unit 10 calculates the reproduction boundary y ref , sampling frequency f s , number of speaker units M, speaker unit interval Δx, buffer length N, number of DFT points F, spherical harmonic spectrum G l m i, the truncation order Nu and the desired boundary y des In addition, the virtual sound source direction (α(n), β(n), γ(n)) and virtual sound source coordinates (r s (n),θ s (n),φ s s(n) and a sound source signal s(n) are input (step S401).

[0073] The x-axis direction wave number calculation unit 20 receives the number of speaker units M and the speaker unit interval Δx that are set in advance. Then, the x-axis direction wave number calculation unit 20 calculates the wave number k in the x-axis direction using the number of speaker units M, the speaker unit interval Δx, and the index m. x (=2πm / (ΔxM)) is calculated (step S402). Here, the index m is, as explained in the above formula (8), m=0 to M-1 (m is an integer).

[0074] The x-axis wave number calculation unit 20 calculates the wave number k x to the reproduced sound field angular spectrum calculation unit 22 and the desired sound field angular spectrum calculation unit 23.

[0075] The wave number calculation unit 21 calculates the wave number at a preset sampling frequency f s and the number of DFT points F. Then, the wave number calculation unit 21 inputs the sampling frequency f s , using the DFT point number F and frequency index l, the wave number k (= ω / c) and the angular frequency ω (= 2πf s l / F) is calculated (step S403), where c is the speed of sound. The frequency index l is l=0 to F-1 (l is an integer).

[0076] The wave number calculation unit 21 outputs the wave number k and the angular frequency ω to the reproduced sound field angular spectrum calculation unit 22 and the desired sound field angular spectrum calculation unit 23.

[0077] (Reproduced sound field angle spectrum calculation unit 22) The reproduced sound field angular spectrum calculation unit 22 calculates the reproduced sound field angular spectrum based on the preset reproduction boundary y ref is input, and the wave number k in the x-axis direction is calculated from the x-axis wave number calculation unit 20. x The wave number k and the angular frequency ω are input from the wave number calculation unit 21.

[0078] The reproduced sound field angular spectrum calculation unit 22 calculates the reproduction boundary y ref , wave number k in the x-axis direction x , the wave number k, the angular frequency ω, and a preset value y0 are used to calculate a reproduced sound field angular spectrum corresponding to the denominator of the above equation (4) (step S404).

[0079] As described above, speaker array 100 is arranged on y=y0 of preset z=z0 in the xyz space, and the value y0 indicates the y value in the xyz space where speaker array 100 is arranged. Then, reproduced sound field angular spectrum calculation unit 22 outputs the reproduced sound field angular spectrum to division unit 24.

[0080] 5 is a block diagram showing an example of the configuration of the reproduced sound field angular spectrum calculation unit 22. The reproduced sound field angular spectrum calculation unit 22 includes a calculation unit 40 and a memory 41.

[0081] The calculation unit 40 calculates the reproduction boundary y ref is input, and the wave number k in the x-axis direction (for the number M of speaker units) is calculated from the x-axis wave number calculation unit 20. x The wave number k and the angular frequency ω are input from the wave number calculation unit 21.

[0082] The calculation unit 40 calculates the reproduced sound field angular spectrum by the following formula, and stores the reproduced sound field angular spectrum in the memory 41.

number

[0083] Specifically, the calculation unit 40 calculates the wave number k in the x-axis direction from the square of the wave number k. x Subtract the squared value of y and find the square root of the result. ref to the value y 0 The calculation unit 40 subtracts the square root of the multiplication result from the second kind of Hankel function H 0 (2) The reproduced sound field angular spectrum is calculated by multiplying the result by (-i / 4). i is the imaginary unit. This reproduced sound field angular spectrum is calculated by dividing the reproduced sound field boundary y ref The angular spectrum G^(k x ,y ref -y 0 ,z 0 ,ω).

[0084] In this way, the reproduced sound field angular spectrum is ref Therefore, the reproduced sound field angular spectrum is calculated using only the angular spectrum D^(k x ,y 0 ,z 0, ω) are calculated in advance and stored in the memory 41.

[0085] (Desired sound field angle spectrum calculation unit 23) Returning to FIG. 3 and FIG. 4, the desired sound field angular spectrum calculation unit 23 calculates the desired sound field angular spectrum G l m i, the truncation order Nu and the desired boundary y des is input, and the wave number k in the x-axis direction (for the number M of speaker units) is calculated from the x-axis wave number calculation unit 20. x The wave number k and the angular frequency ω are input from the wave number calculation unit 21.

[0086] In addition, the desired sound field angular spectrum calculation unit 23 calculates the virtual sound source direction (α(n), β(n), γ(n)) and the virtual sound source coordinates (r s (n),θ s (n),φ s s(n) and sound source signal s(n) are input (step S405).

[0087] The desired sound field angular spectrum calculation unit 23 calculates the desired sound field angular spectrum G l m i, truncation order Nu, desired boundary y des , wave number k in the x-axis direction x , wave number k and angular frequency ω, as well as the virtual sound source direction (α(n), β(n), γ(n)) and virtual sound source coordinates (r s (n),θ s (n),φ s The desired sound field angular spectrum corresponding to the numerator of the formula (4) is calculated by performing the calculations of the formulas (17) and (18) using the sound source signal s(n) and the sound source signal s(n) (step S406).

[0088] The desired sound field angular spectrum calculation unit 23 outputs the desired sound field angular spectrum to the division unit 24 .

[0089] 6 is a block diagram showing an example of the configuration of the desired sound field angular spectrum calculation unit 23. The desired sound field angular spectrum calculation unit 23 includes an input unit 30, a buffer 31, and a calculation unit 32.

[0090] The input unit 30 inputs a preset buffer length N, and also inputs the virtual sound source direction (α(n), β(n), γ(n)) and virtual sound source coordinates (r s (n),θ s (n),φ s The input unit 30 receives the virtual sound source direction (α(n), β(n), γ(n)) and virtual sound source coordinates (r s (n),θ s (n),φ s The signal s(n) and the sound source signal s(n) are stored in the buffer 31.

[0091] In this case, the buffer 31 always stores the virtual sound source direction (α(n), β(n), γ(n)) and virtual sound source coordinates (r s (n),θ s (n),φ s The input unit 30 performs a storage process so that the input signal s(n) and the sound source signal s(n) are stored.

[0092] The input unit 30 receives the virtual sound source direction (α(n), β(n), γ(n)) and virtual sound source coordinates (r s (n),θ s (n),φ s Each time a signal s(n) and a sound source signal s(n) are input, the data stored in the buffer section 0 to L-2 in the buffer 31 is shifted to the buffer section 1 to L-1.

[0093] Then, the input unit 30 calculates the virtual sound source direction (α(n), β(n), γ(n)) and virtual sound source coordinates (r s (n),θ s (n),φ s (n)) and sound source signal s(n) are expressed as the virtual sound source direction (α(0), β(0), γ(0)) and the virtual sound source coordinates (rs (0),θ s (0),φ s The signal s(0) and the sound source signal s(0) are stored in the buffer section 0. That is, the buffer 31 is a buffer that is updated every time sample L corresponding to the buffer length N.

[0094] The calculation unit 32 calculates a predetermined spherical harmonic spectrum G l m i, the truncation order Nu and the desired boundary y des is input, and the wave number k in the x-axis direction is calculated from the x-axis wave number calculation unit 20. x The wave number k and the angular frequency ω are input from the wave number calculation unit 21.

[0095] The calculation unit 32 also calculates from the buffer 31 the virtual sound source direction (α(n), β(n), γ(n)) and the virtual sound source coordinates (r s (n),θ s (n),φ s (n) and the source signal s(n).

[0096] The calculation unit 32 calculates the spherical harmonic spectrum G l m i, truncation order Nu, desired boundary y des , wave number k in the x-axis direction x , wave number k and angular frequency ω, as well as the virtual sound source direction (α(n), β(n), γ(n)) and virtual sound source coordinates (r s (n),θ s (n),φ s By performing the calculations of equations (17) and (18) using sound source signal s(n), the desired sound field angular spectrum is obtained for each buffering section of buffer length N, and the desired sound field angular spectrum is output to division unit 24.

[0097] In this case, when performing the calculations of the formulas (17) and (18), the calculation unit 32 truncates the integral of τ in the formula (17) to a finite length at the value indicated by the truncation order Nu. In addition, the calculation unit 32 calculates the virtual sound source direction (α(n), β(n), γ(n)), the virtual sound source coordinates (r s (n),θ s (n),φ s The calculation unit 32 receives discrete values, which are the input signal s(n) and the sound source signal s(n), and performs a discretized calculation. The calculation unit 32 also truncates the calculation of sigma Σ for the orders l and v to a finite order at a value indicated by a preset truncation order Nu.

[0098] Here, the virtual sound source direction (α(n), β(n), γ(n)) is expressed by the Wigner D function D mm’ l It corresponds to (α(τ),β(τ),γ(τ)) in (α(τ),β(τ),γ(τ)). s (n),θ s (n),φ s (n)) is the function S included in the above formula (17). lν mμ R in (Rs(τ)) s The sound source signal s(n) corresponds to s(τ) included in the above equation (17).

[0099] Also, the spherical harmonic spectrum G l m i corresponds to the left side of the above equation (15), and the spherical harmonic spectrum G l m i, the spherical harmonic spectrum in the above formula (17) TIFF0007674954000039.tif15170 is uniquely determined (given).

[0100] The above-mentioned equation (17) includes an integral with respect to τ, but the virtual sound source direction (α(n), β(n), γ(n)) and the virtual sound source coordinate (r s (n),θ s (n),φ sSince the sound source signal s(n) and the sound source signal s(n) are discrete values ​​for each time sample n, the integral is discretized and truncated as described above.

[0101] In other words, the integral for τ included in the above equation (17) is the virtual sound source direction (α(n), β(n), γ(n)), the virtual sound source coordinates (r s (n),θ s (n),φ s The calculation is performed by discretizing the signal s(n) and the discrete values ​​of the sound source signal s(n).

[0102] In addition, the truncation order Nu indicates a finite interval used instead of the interval [-∞,∞] in the integration for τ included in the above equation (17), and also indicates the truncation order of the calculation of sigma Σ for orders l and ν included in the above equations (17) and (18).

[0103] (Calculation of formula (18)) Specifically, the calculation unit 32 obtains the desired sound field angular spectrum of the above formula (18) by the following calculation. That is, the calculation unit 32 obtains the wave number k in the x-axis direction from the square value of the wave number k. x Subtract the squared value of , find the square root of the subtraction result, and add the square root to the desired boundary y des Multiply by μ and use the second kind Hankel function H μ (2) Then, the calculation unit 32 raises the imaginary unit i to the μ-ν power and applies the μ-order Hankel function H μ (2) Then, the sigma Σ operation from μ=-ν to ν is performed on the multiplication result to obtain the 18-1 sigma Σ operation result.

[0104] The calculation unit 32 calculates the Wigner D function D μμ’ ν (0,π / 2,0) and apply the spherical harmonic spectrum of the above formula (17) to it. TIFF0007674954000040.tif16170 is multiplied, and the sigma Σ operation from μ'=-ν to ν is performed on the multiplication result to obtain the 18-2nd sigma Σ operation result.

[0105] The calculation unit 32 multiplies the 18-1st sigma Σ calculation result and the 18-2nd sigma Σ calculation result, and performs sigma Σ calculation on the multiplication result from ν=0 to the order indicated by the truncation order Nu to obtain the 18-3rd sigma Σ calculation result. The calculation unit 32 then multiplies the result of dividing π by the wave number k by the 18-3rd sigma Σ calculation result to obtain the desired sound field angular spectrum indicated in the above formula (18).

[0106] (Calculation of formula (17)) Here, the calculation unit 32 calculates the spherical harmonic spectrum of the above formula (17). TIFF0007674954000041.tif17170 is calculated by the following calculation. <r s In the case of (τ), (R is the coordinate of the receiving point, R S (τ) indicates the coordinates of the sound source. ) The Wigner D function D of (α(τ),β(τ),γ(τ)) where (α(n),β(n),γ(n)) is the virtual sound source direction. mm’ l The virtual sound source coordinates (r s (n),θ s (n),φ s (n)) s (τ)=(r s (τ),θ s (τ),φ s (τ)) function S of the above formula (13) lν mμ (In the above formula (13), b = R s (τ)). The calculation unit 32 calculates the Wigner D function D mm’ l , the function S lν mμ , s(τ), which is the sound source signal s(n), and a complex number having the multiplication result of the angular frequency ω and the time τ as the argument, to obtain a 17-1-1 multiplication result.

[0107] The calculation unit 32 integrates the 17-1-1 multiplication result in the interval of time τ indicated by the truncation order Nu. Then, the calculation unit 32 multiplies the wave number k and the distance r from the origin to the sound receiving point, and calculates the ν-th order spherical Bessel function jν Then, divide the coordinate R of the sound receiving point by the distance r from the origin to the sound receiving point, and calculate the ν-th order μ-th order spherical harmonic function Y ν μ Request.

[0108] The calculation unit 32 calculates the spherical harmonic spectrum G l m The spherical harmonic spectrum corresponding to i TIFF0007674954000042.tif17170, the integration result, the ν-th order spherical Bessel function j ν and the ν-th order μ-order spherical harmonic function Y ν μ The 17th-1-2 multiplication result is obtained by multiplying by .

[0109] The calculation unit 32 performs a sigma Σ calculation from m'=-l to l, a sigma Σ calculation from m=-l to l, and a sigma Σ calculation from l=0 to the order indicated by the truncation order Nu on the 17-1-2 multiplication result, thereby obtaining r <r s We obtain the spherical harmonic spectrum for (τ).

[0110] On the other hand, the calculation unit 32 determines whether r>r s In the case of (τ), the Wigner D function D of (α(τ),β(τ),γ(τ)) with the virtual sound source direction (α(n),β(n),γ(n)) mm’ l The virtual sound source coordinates (r s (n),θ s (n),φ s (n)) s (τ)=(r s (τ),θ s (τ),φ s (τ)) function S^ of the above equation (14) lν mμ (In the above formula (14), b = R s (τ)). The calculation unit 32 calculates the Wigner D function D mm’ l , the function S lν mμ, s(τ), which is the sound source signal s(n), and a complex number having the multiplication result of the angular frequency ω and the time τ as the argument, to obtain a 17-2-1 multiplication result.

[0111] The calculation unit 32 integrates the 17-2-1 multiplication result in the interval of time τ indicated by the truncation order Nu. Then, the calculation unit 32 multiplies the wave number k and the distance r from the origin to the sound receiving point, and calculates the ν-ordered two-dimensional spherical Hankel function h ν (2) Then, divide the coordinate R of the sound receiving point by the distance r from the origin to the sound receiving point, and calculate the ν-th order μ-th order spherical harmonic function Y ν μ Request.

[0112] The calculation unit 32 calculates the spherical harmonic spectrum G l m The spherical harmonic spectrum corresponding to i TIFF0007674954000043.tif16170, the integration result, the ν-dependent second-order spherical Hankel function h ν (2) and the ν-th order μ-order spherical harmonic function Y ν μ The 17-2-2 multiplication result is obtained by multiplying by .

[0113] The calculation unit 32 performs a sigma Σ operation from m'=-l to l, a sigma Σ operation from m=-l to l, and a sigma Σ operation from l=0 to the order indicated by the truncation order Nu on the 17-2-2 multiplication result, thereby obtaining r>r shown in the lower part of the above formula (17). s We obtain the spherical harmonic spectrum for (τ).

[0114] (Calculation of formula (13)) In addition, the calculation unit 32 uses the function S lν mμ (b=R s (τ) is calculated by the following calculation: That is, the calculation unit 32 raises the imaginary unit i to the ν-μ power and multiplies the result by 4π to obtain the 13-1st multiplication result.

[0115] The calculation unit 32 obtains the 13-2 multiplication result by raising the imaginary unit i to the q power, and obtains the 13-3 multiplication result by raising -1 to the m power. Then, the calculation unit 32 calculates the virtual sound source coordinates (r s (n),θ s (n),φ s (n)) s The q-ordered two-kind spherical Hankel function h for the multiplication result obtained by multiplying the absolute value of (τ) q (2) The virtual sound source coordinates (r s (n),θ s (n),φ s (n)) s (τ) is the R s The q-th order μ-m-order spherical harmonic function Y for the division result obtained by dividing by the absolute value of (τ) q μ-m The calculation unit 32 also calculates a Gaunt coefficient g with the degree and order of the spherical harmonic function l, m, v, -μ, and q as variables.

[0116] The calculation unit 32 calculates the 13-2 multiplication result, the 13-3 multiplication result, and the q-order binary spherical Hankel function h q (2) , the q-order μ-m spherical harmonic function Y q μ-m and the gauntlet coefficient g to obtain the 13-4th multiplication result.

[0117] The calculation unit 32 performs a sigma Σ calculation from q=0 to |m|+|μ| on the 13-4th multiplication result to obtain a 13-1st sigma Σ calculation result.

[0118] The calculation unit 32 multiplies the 13-1st multiplication result by the 13-1st sigma Σ calculation result to obtain the function S lν mμ get.

[0119] (Calculation of formula (14)) In addition, the calculation unit 32 uses the function S^ in the formula (14) in the calculation of the spherical harmonic spectrum in the formula (17). lν mμ(b=R s (τ) is calculated by the following calculation: That is, the calculation unit 32 raises the imaginary unit i to the ν-μ power and multiplies the result by 4π to obtain the 14-1st multiplication result.

[0120] The calculation unit 32 obtains the 14-2 multiplication result by raising the imaginary unit i to the q power, and obtains the 14-3 multiplication result by raising -1 to the m power. Then, the calculation unit 32 calculates the virtual sound source coordinates (r s (n),θ s (n),φ s (n)) s The qth order spherical Bessel function j for the multiplication result obtained by multiplying the absolute value of (τ) q The virtual sound source coordinates (r s (n),θ s (n),φ s (n)) s (τ) is the R s The q-th order μ-m-order spherical harmonic function Y for the division result obtained by dividing by the absolute value of (τ) q μ-m The calculation unit 32 also calculates a Gaunt coefficient g with the degree and order of the spherical harmonic function l, m, v, -μ, and q as variables.

[0121] The calculation unit 32 calculates the 14-2 multiplication result, the 14-3 multiplication result, and the q-th order spherical Bessel function j q , the q-order μ-m spherical harmonic function Y q μ-m and the gauntlet coefficient g to obtain the 14-4th multiplication result.

[0122] The calculation unit 32 performs a sigma Σ calculation from q=0 to |m|+|μ| on the 14-4th multiplication result to obtain a 14-1st sigma Σ calculation result.

[0123] The calculation unit 32 multiplies the 14-1st multiplication result by the 14-1st sigma Σ calculation result to obtain the function S^ in the above equation (14). lν mμ get.

[0124] (Division part 24) 3 and 4, the division unit 24 inputs the reproduced sound field angular spectrum from the reproduced sound field angular spectrum calculation unit 22, and also inputs the desired sound field angular spectrum for each buffer length N from the desired sound field angular spectrum calculation unit 23. Then, the division unit 24 divides the desired sound field angular spectrum for each buffer length N by the reproduced sound field angular spectrum, as in the above formula (4), to obtain the angular spectrum D^(k x ,y 0 ,z 0 ,ω). The division unit 24 calculates the angular spectrum D^(k x ,y 0 ,z 0 , ω) to the spatial frequency domain discrete inverse Fourier transform unit 11 (step S407).

[0125] Unless a predetermined end condition is satisfied (step S408: N), the drive signal calculation unit 10 proceeds to step S405. Then, the drive signal calculation unit 10 calculates the angular spectrum D^(k x ,y 0 ,z 0 The process of steps S405 to S407 for obtaining ω is repeated. On the other hand, if a predetermined end condition is satisfied (step S408: Y), the drive signal calculation section 10 ends the process.

[0126] As described above, according to the sound field reproduction device 1 of the embodiment of the present invention, the drive signal calculation unit 10 calculates the wave number k in the x-axis direction based on the preset number M of speaker units and the speaker unit interval Δx. x Calculate the preset sampling frequency f s Based on the number of DFT points F and the like, the wave number k and the angular frequency ω are calculated.

[0127] The drive signal calculation unit 10 calculates the reproduction boundary y ref , the calculated wave number in the x-axis direction k x Using the wave number k and the angular frequency ω, the reproduced sound field angular spectrum G^(k x ,y ref-y 0 ,z 0 ,ω) is calculated.

[0128] The drive signal calculation unit 10 calculates the output of the drive signal by calculating a predetermined buffer length N and a spherical harmonic spectrum G l m i, the truncation order Nu and the desired boundary y des , the calculated wave number in the x-axis direction k x , wave number k and angular frequency ω, as well as the virtual sound source direction (α(n), β(n), γ(n)) and virtual sound source coordinates (r s (n),θ s (n),φ s (n)) and the sound source signal s(n), the desired sound field angular spectrum P^ is calculated using the equations (17) and (18). d (k x ,y des ,z 0 ,ω) is calculated.

[0129] The drive signal calculation unit 10 calculates the desired sound field angular spectrum P^ d (k x ,y des ,z 0 ,ω) and the reproduced sound field angular spectrum G^(k x ,y ref -y 0 ,z 0 , ω) in the above formula (4), the wave number k in the x-axis direction is x Angular spectrum of the drive signal for each D^(k x ,y 0 ,z 0 ,ω) is calculated.

[0130] The spatial frequency domain discrete inverse Fourier transform unit 11 calculates the wave number k in the x-axis direction. x Angular spectrum of the drive signal for each D^(k x ,y 0 ,z 0 By performing a spatial frequency domain discrete inverse Fourier transform on x, ω), a drive signal D~(x, ω) in the time frequency domain for each speaker unit is obtained.

[0131] The time-frequency domain discrete inverse Fourier transform unit 12 performs a time-frequency domain discrete inverse Fourier transform on the time-frequency domain drive signal D~(x,ω) for each speaker unit to obtain a time-domain speaker drive signal for each speaker unit. The time-frequency domain discrete inverse Fourier transform unit 12 then outputs the speaker drive signal to the speaker units that constitute the speaker array 100.

[0132] In this manner, in the embodiment of the present invention, when generating drive signals for the speaker array 100 using SDM, it is possible to reproduce a sound field formed by a moving directional sound source.

[0133] Although the present invention has been described above with reference to the embodiments, the present invention is not limited to the above-described embodiments and can be modified in various ways without departing from the technical concept thereof.

[0134] A normal computer can be used as the hardware configuration of the sound field reproduction device 1 according to the embodiment of the present invention. The sound field reproduction device 1 is configured by a computer equipped with a CPU, a volatile storage medium such as a RAM, a non-volatile storage medium such as a ROM, an interface, etc.

[0135] The functions of the drive signal calculation unit 10, the spatial frequency domain discrete inverse Fourier transform unit 11 and the time frequency domain discrete inverse Fourier transform unit 12 provided in the sound field reproduction device 1 are each realized by causing a CPU to execute a program in which these functions are written.

[0136] These programs are stored in the storage medium and are read and executed by the CPU. These programs can also be distributed by storing them in a storage medium such as a magnetic disk (floppy disk, hard disk, etc.), an optical disk (CD-ROM, DVD, etc.), or a semiconductor memory, and can also be transmitted and received via a network. [Explanation of symbols]

[0137] 1. Sound field reproduction device 10 Drive signal calculation unit 11 Spatial frequency domain discrete inverse Fourier transform section 12 Time-frequency domain discrete inverse Fourier transform section 20 Wave number calculation section in the x-axis direction 21 Wave number calculation section 22 Reproduced sound field angle spectrum calculation unit 23 Desired sound field angle spectrum calculation unit 24 Division section 30 Input section 31 Buffer 32,40 Calculation section 41 Memory 100-1,100-2,100-3,100-4,100 Speaker array 101 Infinite Linear Sound Source R Coordinates of the sound receiving point R s (t) Coordinates of the sound source D^(k x ,y 0 ,z 0 ,ω) Angular spectrum of the drive signal P^ d (k x ,y des ,z 0 ,ω) Desired sound field angular spectrum G^(k x ,y ref -y 0 ,z 0 ,ω) Reproduced sound field angular spectrum y ref Reproduced Boundary f s Sampling Frequency M Number of speaker units Δx Speaker unit spacing N buffer length F DFT score (α(n),β(n),γ(n)) Virtual sound source direction (r s (n),θ s (n),φ s (n)) Virtual sound source coordinates s(n) sound source signal G l mi Spherical harmonic spectrum Nu truncation order y des desired boundary D~(x,ω) driving signal in the time-frequency domain n time samples k x Wave number along the x-axis ω angular frequency k wavenumber

Claims

1. A sound field reproduction device that generates a drive signal for reproducing a sound field formed by a moving directional sound source using a speaker array consisting of a plurality of speaker units, The speaker array is located at z=z in the xyz space. 0 y=y 0 and y ref is the reproduction boundary, f s is the sampling frequency, M is the number of speaker units, Δx is the speaker unit interval, N is the buffer length, F is the number of DFT points, Nu is the truncation order, y des Let G be the desired boundary. l m i is the sampling frequency f s and an l (l = 0 to L)-th order spherical harmonic spectrum (i is a frequency bin number, l is an order, and m is an order) representing the radiation characteristics of the directional sound source for each frequency bin determined by the DFT point number F; k ｘ is the wave number in the x-axis direction, k is the wave number, ω is the angular frequency, and c is the sound speed. n is a time sample, (α(n), β(n), γ(n)) is a virtual sound source direction, which is the direction of the directional sound source for each time sample n, expressed in ZYZ Euler angles, and (r s (n), θ s (n), φ s (n)) is a virtual sound source coordinate, which is the coordinate of the directional sound source for each time sample n, (r s (n) is the distance from the origin to the virtual sound source, θ s (n) is the zenith angle, φ s s(n) is a sound source signal that is a signal of the directional sound source for each time sample n, and s(n) is an azimuth angle. From the preset number of speaker units M, the speaker unit interval Δx, and an index m=0 to M−1 (m is an integer), the formula: k ｘ = 2πm / (ΔxM), the wave number k in the x-axis direction ｘ An x-axis direction wave number calculation unit that calculates The preset sampling frequency f s From the number of DFT points F and the frequency index l = 0 to F-1 (l is an integer), the formula: ω = 2πf s A wave number calculation unit that calculates the angular frequency ω by l / F and calculates the wave number k by the equation: k = ω / c; H 0 (2) is the Hankel function of the second kind with 0, and i is the imaginary unit. The wave number k in the x-axis direction calculated by the x-axis direction wave number calculation unit ｘ , the wave number k calculated by the wave number calculation unit, and the preset reproduction boundary y ref and the value y 0 from, formula: a reproduced sound field angular spectrum calculation unit for calculating a reproduced sound field angular spectrum by a desired sound field angular spectrum calculation unit that calculates a desired sound field angular spectrum; a division unit that obtains an angular spectrum of the drive signal by dividing the desired sound field angular spectrum calculated by the desired sound field angular spectrum calculation unit by the reproduced sound field angular spectrum calculated by the reproduced sound field angular spectrum calculation unit, The desired sound field angular spectrum calculation unit is The predetermined spherical harmonic spectrum G l m i, the truncation order Nu and the desired boundary y des , the wave number k in the x-axis direction calculated by the x-axis direction wave number calculation unit ｘ , the wave number k and the angular frequency ω calculated by the wave number calculation unit, and the virtual sound source direction (α(n), β(n), γ(n)) and the virtual sound source coordinates (r s (n), θ s (n), φ s (n)) and the sound source signal s(n), formula: (P^ d (k x , y des , ω) is the desired sound field angular spectrum, ν and μ are the degree and order of a predetermined spherical harmonic function, i is the imaginary unit, H μ (2) Let μ be the second kind of Hankel function, D μμ’ ν is the Wigner-D function, is a predetermined spherical harmonic spectrum, and the sigma Σ operation of order v is truncated by the truncation order Nu. formula: (l, m are the degree and order of a given spherical harmonic function, respectively, The spherical harmonic spectrum G l m i, the spherical harmonic spectrum, D mm’ l is the Wigner D function, S lν mμ , S^ lν mμ Let R be a given function. s (τ) is the virtual sound source coordinate (r s (n), θ s (n), φ s (n)), s(τ) is the sound source signal s(n), j at time τ ν is the ν-th order spherical Bessel function, h ν (2) Let us use the ν-ordered two-kind spherical Hankel function, Y ν μ is a μ-th order spherical harmonic function, R is the coordinate of the sound receiving point, and r is the distance from the origin to the sound receiving point. The integral of dτ and the sigma Σ operation of the order l are truncated by the truncation order Nu. formula: (h q (2) The two-kind spherical Hankel function, Y q μ-ｍ ̂ is the complex conjugate of the qth order μ-mth order spherical harmonic function, g(l,m;ν,-μ;q) is the gaunt coefficient, and b = R s (Let τ be the square root of the formula: (j q is the qth order spherical Bessel function, Y q μ-ｍ Let ̂ be the complex conjugate of the qth μ-mth order spherical harmonic function, g(l, m; ν, -μ; q) be the Gaunt coefficient, and b = R s (Let τ be the square root of the and calculating the desired sound field angular spectrum by:

2. A computer constituting a sound field reproduction device that generates a drive signal for reproducing a sound field formed by a moving directional sound source using a speaker array consisting of a plurality of speaker units, The speaker array is located at z=z in the xyz space. 0 y=y 0 and y ref is the reproduction boundary, f s is the sampling frequency, M is the number of speaker units, Δx is the speaker unit interval, N is the buffer length, F is the number of DFT points, Nu is the truncation order, y des Let G be the desired boundary. l m i is the sampling frequency f s and an l (l = 0 to L)-th order spherical harmonic spectrum (i is a frequency bin number, l is an order, and m is an order) representing the radiation characteristics of the directional sound source for each frequency bin determined by the DFT point number F; k ｘ is the wave number in the x-axis direction, k is the wave number, ω is the angular frequency, and c is the sound speed. n is a time sample, (α(n), β(n), γ(n)) is a virtual sound source direction, which is the direction of the directional sound source for each time sample n, expressed in ZYZ Euler angles, and (r s (n), θ s (n), φ s (n)) is a virtual sound source coordinate, which is the coordinate of the directional sound source for each time sample n, (r s (n) is the distance from the origin to the virtual sound source, θ s (n) is the zenith angle, φ s s(n) is a sound source signal that is a signal of the directional sound source for each time sample n, and s(n) is an azimuth angle. From the preset number of speaker units M, the speaker unit interval Δx, and an index m=0 to M−1 (m is an integer), the formula: k ｘ = 2πm / (ΔxM), the wave number k in the x-axis direction ｘ An x-axis direction wave number calculation unit that calculates The preset sampling frequency f s From the number of DFT points F and the frequency index l = 0 to F-1 (l is an integer), the formula: ω = 2πf s a wave number calculation unit that calculates the angular frequency ω by ω / F and calculates the wave number k by the equation: k=ω / c; H 0 (2) is the Hankel function of the second kind with 0, and i is the imaginary unit. The wave number k in the x-axis direction calculated by the x-axis direction wave number calculation unit ｘ , the wave number k calculated by the wave number calculation unit, and the preset reproduction boundary y ref and the value y 0 from, formula: a reproduced sound field angular spectrum calculation unit for calculating a reproduced sound field angular spectrum by A desired sound field angular spectrum calculation unit that calculates a desired sound field angular spectrum; and a division unit that divides the desired sound field angular spectrum calculated by the desired sound field angular spectrum calculation unit by the reproduced sound field angular spectrum calculated by the reproduced sound field angular spectrum calculation unit to obtain the angular spectrum of the drive signal; The desired sound field angular spectrum calculation unit is The predetermined spherical harmonic spectrum G l m i, the truncation order Nu and the desired boundary y des , the wave number k in the x-axis direction calculated by the x-axis direction wave number calculation unit ｘ , the wave number k and the angular frequency ω calculated by the wave number calculation unit, and the virtual sound source direction (α(n), β(n), γ(n)) and the virtual sound source coordinates (r s (n), θ s (n), φ s (n)) and the sound source signal s(n), formula: (P^ d (k x , y des , ω) is the desired sound field angular spectrum, ν and μ are the degree and order of a predetermined spherical harmonic function, i is the imaginary unit, H μ (2) Let μ be the second kind of Hankel function, D μμ’ ν is the Wigner-D function, is a predetermined spherical harmonic spectrum, and the sigma Σ operation of order v is truncated by the truncation order Nu. formula: (l, m are the degree and order of a given spherical harmonic function, respectively, The spherical harmonic spectrum G l m i, the spherical harmonic spectrum, D mm’ l is the Wigner D function, S lν mμ , S^ lν mμ Let R be a given function. s (τ) is the virtual sound source coordinate (r s (n), θ s (n), φ s (n)), s(τ) is the sound source signal s(n), j at time τ ν is the ν-th order spherical Bessel function, h ν (2) Let us use the ν-ordered two-kind spherical Hankel function, Y ν μ is a μ-th order spherical harmonic function, R is the coordinate of the sound receiving point, and r is the distance from the origin to the sound receiving point. The integral of dτ and the sigma Σ operation of the order l are truncated by the truncation order Nu. formula: (h q (2) The two-kind spherical Hankel function, Y q μ-ｍ ̂ is the complex conjugate of the qth order μ-mth order spherical harmonic function, g(l,m;ν,-μ;q) is the gaunt coefficient, and b = R s (Let τ be the square root of the formula: (j q is the qth order spherical Bessel function, Y q μ-ｍ Let ̂ be the complex conjugate of the qth μ-mth order spherical harmonic function, g(l, m; ν, -μ; q) be the Gaunt coefficient, and b = R s (Let τ be the square root of the and calculating the desired sound field angular spectrum by the above-mentioned method.

Images