CN105785320A - Function type delay summation method for identifying solid sphere array three-dimensional sound source - Google Patents

Abstract

The invention discloses a function type delay summation method for identifying a sold sphere array three-dimensional sound source, wherein a novel FDAS sound source identification method is established. In identifying incoherent sound sources in a three-dimensional space, the function type delay summation method supplies an effective approach for quickly obtaining a clean sound source imaging graph and realizes an important meaning in accurate analysis on a sound source identification result.

Description

Functional delay summation method for solid sphere array three-dimensional sound source identification

Technical Field

The invention relates to the field of three-dimensional space sound source identification

Background

The SHB is a sound source identification method commonly used in a solid sphere array, and theoretically, an output function is constructed as an ideal dirac function based on the orthogonality of spherical harmonics, and only a maximum value is output at a real sound source position, and zeros are output at other non-sound source positions. However, in practical application, due to the spherical harmonic order truncation, the discrete sampling of the microphone and other factors, the method outputs a large amount of high-level side lobes at a non-sound source position besides a main lobe at a real sound source position, and the side lobes can cover a weak sound source and cannot be identified, and can be mistakenly considered as a sound source, so that an analysis result is subjected to serious uncertainty. In order to obtain definite three-dimensional sound source identification results, many scholars have devoted themselves to exploring SHB side lobe attenuation paths, and typical methods are filtering summation, deconvolution and the like. The former is introduced from Hald in 2013, the calculation efficiency is high, but the side lobe attenuation effect is not obvious; the latter was introduced by the inventors in 2015, which can attenuate the side lobe more significantly, but has a serious time-consuming problem in practical application. Therefore, there is a need to explore new side lobe attenuation approaches that are both powerful and fast.

In 2014, aiming at the problem of sidelobe attenuation of planar array two-dimensional sound source identification, the Dougherty provides a novel FB method with dual advantages of strong sidelobe attenuation capability and high calculation efficiency. The method is based on the idea that an exponential function with the base number larger than 0 and smaller than 1 decreases with the increase of an exponent to attenuate a side lobe, the process is simple, only a cross-spectrum matrix of an array microphone receiving sound pressure signals needs to be decomposed by characteristic values, exponential parameters are introduced, the original characteristic values are replaced by the inverse power of the exponential parameters of the original characteristic values to construct a new cross-spectrum matrix, the new cross-spectrum matrix is processed by adopting the existing sound source identification method, and finally the exponential parameter power of the obtained result is calculated. If the FB method is expanded to the three-dimensional sound source identification of the solid sphere array, a clean sound source imaging image can be expected to be obtained quickly, and the FB method has important significance for quickly and accurately identifying the sound source in the three-dimensional space. However, the SHB method does not satisfy the premise that the output value of the point spread function at the sound source position is 1 and the output value at the non-sound source position is less than 1.

Disclosure of Invention

The invention aims to solve the problem that the SHB method does not meet the premise that the output value of a point spread function at a sound source position is 1 and the output value at a non-sound source position is less than 1, and provides a functional delay summation method for solid sphere array three-dimensional sound source identification.

The technical scheme adopted for realizing the purpose of the invention is as follows, namely a functional delay summation method for solid sphere array three-dimensional sound source identification:

1): constructing a measuring system:

referring to fig. 1, an array of solid spheres in three-dimensional space; the solid sphere array comprises a solid sphere with the radius of a and Q microphones distributed on the surface of the solid sphere; the microphones are numbered Q, Q1, 2, …, Q; the sphere center of the solid sphere is used as an origin, the position coordinates of any point in the three-dimensional space where the solid sphere array is located are described by (r, omega), r represents the distance between the described position and the origin,denotes the direction of the depicted position, theta,Respectively its elevation angle and azimuth angle; (a, Ω)_q) Position coordinates of a microphone with the number q;

2) obtaining sound pressure

Each microphone receives a sound pressure signal p (ka, Ω)_q) (ii) a The monopole point sound source in the three-dimensional space radiates sound waves as spherical waves; wave number k is 2 pi f/c, sound frequency f, sound velocity c, (a, omega)_q) Position coordinates of a microphone with the number q;

3) receiving sound pressure signals through the microphones to obtain a cross-spectrum matrix C:

$C=\stackrel{\‾}{{\mathrm{pp}}^H},$

wherein: b represents a set of all sound source position coordinates, (r)₀,Ω₀) And (r)₀',Ω₀') each represent the position coordinates of the sound source,

$p={\mathrm{\Σ}}_{(r_0,{\mathrm{\Ω}}_0)\∈B}s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0),$

s(kr₀,Ω₀) Which is indicative of the intensity of the sound source,

t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),…,t_Q(kr₀,Ω₀)]^T，t_q(kr₀,Ω₀) Representing the sound field transfer function from the sound source to the q microphone,

$t_q({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\≡\underset{n=0}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-n}{\overset{n}{\Σ}}{R}_n({\mathrm{kr}}_0,ka)Y_n^{m*}\left({\mathrm{\Ω}}_0\right)Y_n^m\left({\mathrm{\Ω}}_q\right),$

are respectively omega₀、Ω_qThe spherical harmonics of the directions, n and m are the orders of the spherical harmonics,

the superscripts "-", "T" and "H" respectively denote an averaging operation, a transposition operation and a transposition conjugate operation;

4) determining the spherical harmonic order truncation length corresponding to each frequency:

indicating rounding the value to the nearest integer toward positive infinity;

5) substituting the cross-spectrum matrix C and the spherical harmonic order truncation length N into:

$W_F({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0){\left(\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{\left|\right|v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)\left|{|}_2^2\right|\left|t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\right|{|}_2^2}\right)}^{\mathrm{\ξ}}---\left(I\right)$

wherein: p_AC(kr₀,Ω₀) Is the trace of the cross-spectrum matrix C, (r)_f,Ω_f) As a position coordinate of the focal point, v (kr)_f,Ω_f)＝[v₁(kr_f,Ω_f),v₂(kr_f,Ω_f),…,v_Q(kr_f,Ω_f)]^T，v_q(kr_f,Ω_f) The focus component of the microphone of the q-th order is represented,

6) according to W_F(kr_f,Ω_f) The method can locate the position of the sound source and draw a sound source imaging graph to realize three-dimensional sound field visualization.

The derivation process of the method of the invention is as follows:

step 1: obtaining cross-spectrum matrix of each microphone receiving signal

FIG. 1 shows a spherical array beamforming coordinate system, where the origin is located at the center of the array, and any position in three-dimensional space can be described by (r, Ω), where r represents the distance between the described position and the origin,denotes the direction of the depicted position, theta,Respectively its elevation angle and azimuth angle. The symbol "●" represents microphones embedded in a solid spherical surface, where a is the sphere radius, Q is the total number of microphones, Q is 1,2, …, Q is the microphone index, and (a, Ω)_q) And the spherical array microphone spatially samples sound pressure signals for the position coordinates of the q microphone. The sound source is assumed to be a monopole point sound source, and the radiated sound wave is spherical wave. The acoustic frequency is f, the sound velocity is c, the wave number k is 2 pi f/c, (r)₀,Ω₀) Indicating the position coordinates of the sound source, t_q(kr₀,Ω₀) Representing sound field transfer function from sound source to q-microphoneThe number of the particles is, according to the scattering theory,

$t_q({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\≡\underset{n=0}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-n}{\overset{n}{\Σ}}{R}_n({\mathrm{kr}}_0,ka)Y_n^{m*}\left({\mathrm{\Ω}}_0\right)Y_n^m\left({\mathrm{\Ω}}_q\right)---\left(1\right)$

wherein,are respectively omega₀、Ω_qDirectional spherical harmonics, n, m being the order of the spherical harmonics, superscript "+" denoting the conjugate operation, R_n(kr₀Ka) is the radial function of acoustic signal propagation:

$R_n({\mathrm{kr}}_0,ka)=-4{\mathrm{\πih}}_n^{\left(2\right)}\left({\mathrm{kr}}_0\right)\left(j_n\right(ka)-\frac{j_n^{\′}\left(ka\right)}{h_n^{(2{)}^{\′}}\left(ka\right)}{h}_n^{\left(2\right)}(ka\left)\right)---\left(2\right)$

wherein, the subtractive term is a radial function of spherical wave propagating in a free field, the subtractive term is a spherical wave scattering radial function caused by a solid sphere,is an imaginary unit, j_n(ka) is an n-th order Bessel function of the first kind,is a second class of global Hank function of order n, j'_n(ka)、Are respectively j_n(ka)、The first derivative of (a).

With s (kr)₀,Ω₀) Representing the intensity of the sound source, B representing the set of position coordinates of all sound sources, the sound pressure signal p (ka, Ω) received by the q microphone_q) Can be written as:

$p(ka,{\mathrm{\Ω}}_q)={\mathrm{\Σ}}_{(r_0,{\mathrm{\Ω}}_0)\∈B}s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t_q({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)---\left(3\right)$

let p ═ p (ka, Ω)₁),p(ka,Ω₂),…,p(ka,Ω_Q)]^T、t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),…,t_Q(kr₀,Ω₀)]^TRespectively are microphone sound pressure signal column vector and sound field transfer function column vector, the superscript "T" represents transposition operation, corresponding to formula (3),

$p={\mathrm{\Σ}}_{(r_0,{\mathrm{\Ω}}_0)\∈B}s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)---\left(4\right)$

the cross-spectrum matrix C of each microphone receiving the sound pressure signal is written as:

$C=\stackrel{\‾}{{\mathrm{pp}}^H}={\Σ}_{(r_0,{\mathrm{\Ω}}_0)\∈B}{\Σ}_{({r_0}^{\′},{{\mathrm{\Ω}}_0}^{\′})\∈B}\stackrel{\‾}{s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)s^{*}({{\mathrm{kr}}_0}^{\′},{{\mathrm{\Ω}}_0}^{\′})}t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({{\mathrm{kr}}_0}^{\′},{{\mathrm{\Ω}}_0}^{\′})---\left(5\right)$

wherein, the superscripts "-", respectively, represent the average operation and the transposition conjugate operation, (r)₀',Ω₀') are also sound source position coordinates.

Step 2: obtaining a solid sphere delay summation output function

Process 201 defines solid sphere DAS output functions

The output function of the inventive setup method is defined as:

$W({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=\frac{1}{Q}\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)Cv({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{\left|\right|v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)|{|}_2^2}---\left(6\right)$

wherein (r)_f,Ω_f) Is the position coordinate of the focal point,

v(kr_f,Ω_f)＝[v₁(kr_f,Ω_f),v₂(kr_f,Ω_f),…,v_Q(kr_f,Ω_f)]^Tto focus the column vector, | | | | non-conducting phosphor₂Representing a 2 norm. v (kr)_f,Ω_f) In, element v_q(kr_f,Ω_f) Represents the focus component of a microphone with q, defined as:

$v_q({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=\underset{n=0}{\overset{N}{\Σ}}\underset{m=-n}{\overset{n}{\Σ}}{R}_n({\mathrm{kr}}_f,ka)Y_n^{m*}\left({\mathrm{\Ω}}_f\right)Y_n^m\left({\mathrm{\Ω}}_q\right)---\left(7\right)$

wherein,is omega_fSpherical harmonics of direction, R_n(kr_fKa) is a focusing radial function, and the expression is similar to the expression (2), only r in the expression (2) is required₀Is replaced by r_fThat is, N is the spherical harmonic order truncation length and needs to be set manually. Formula (6) can also be written as:

$W({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=\underset{q,q^{\′}=1}{\overset{Q}{\Σ}}\frac{C_{{\mathrm{qq}}^{\′}}{v}_q^{*}({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)v_{q^{\′}}({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{Q\left|\right|v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)|{|}_2^2}---\left(8\right)$

where q' is also the microphone index, C_qqThe 'is a cross spectrum of a sound pressure signal received by a q-type microphone relative to a sound pressure signal received by a q' -type microphone. Obviously, the method of the invention is firstly usedFor cross spectrum C_qq' Compensation of phase delay, reuseAnd weighting the compensation result, and finally adding all cross-spectrum processing results. The process is similar to the data processing process of the delay summation method for identifying the sound source of the two-dimensional plane array and the three-dimensional hollow sphere array, so the method can be called as delay summation for identifying the three-dimensional sound source of the solid sphere array.

Process 202 deformed solid sphere DAS output function

Defining the average sound pressure contribution as the average value of the self-spectrum of the sound pressure signal generated by the sound source at each microphone of the array, and recording the average value as P_AC(kr₀,Ω₀) Equal to the trace average of the cross-spectral matrix of the sound pressure signal received by the microphone. For a single sound source, the cross-spectral matrix C can be written as:

C＝S(kr₀,Ω₀)t(kr₀,Ω₀)t^H(kr₀,Ω₀)(9)

wherein,is defined as the intensity of the sound source with a power meaning. Accordingly, the number of the first and second electrodes,

$P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)=\frac{1}{Q}tr\left(C\right)=\frac{1}{Q}S({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\left|\right|t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)|{|}_2^2---\left(10\right)$

where "tr" represents the trace of the matrix. According to formula (10), with P_AC(kr₀,Ω₀) Denotes S (kr)₀,Ω₀) The progeny of formula (9) can be obtained:

$C=\frac{P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)Q}{\left|\right|t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)|{|}_2^2}t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)---\left(11\right)$

substituting formula (11) for formula (6) to obtain:

$W({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{\left|\right|v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)\left|{|}_2^2\right|\left|t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\right|{|}_2^2}---\left(12\right)$

the PSF is defined as the response of the sound source identification method to a monopole point sound source contributing to a unit mean sound pressure, and the PSF expression of the DAS method is as follows according to equation (12):

$psf\left(\right({\mathrm{kr}}_f,{\mathrm{\Ω}}_f\left)\right|\left({\mathrm{kr}}_0,{\mathrm{\Ω}}_0\right))=\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{\left|\right|v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)\left|{|}_2^2\right|\left|t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\right|{|}_2^2}---\left(13\right)$

for a single sound source, the output of this DAS method is equal to the product of the average sound pressure contribution of the sound source and the corresponding point spread function.

Process 203 determines spherical harmonic order truncation length

According to the cauchy inequality:

$v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)\≤\left|\right|v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)\left|{|}_2^2\right|\left|t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\right|{|}_2^2---\left(14\right)$

the condition that "═" holds is v (kr)_f,Ω_f)＝t(kr₀,Ω₀). Contrast vector v (kr)_f,Ω_f) And t (kr)₀,Ω₀) The expression of (2), i.e. formula (7) and formula (1), is evident: at a specific frequency, v (kr)_f,Ω_f) Whether or not to be equal to t (kr)₀,Ω₀) Not only with (r)_f,Ω_f) Whether or not to be equal to (r)₀,Ω₀) And the method also relates to the value of the spherical harmonic order truncation length N in the reverse focusing. If a reasonable value of N can be determined, then (r)_f,Ω_f)≠(r₀,Ω₀) When, v (kr)_f,Ω_f)≠t(kr₀,Ω₀) Is established while v (kr) is being established₀,Ω₀)＝t(kr₀,Ω₀) If so, according to the formula (13) and the formula (14), the PSF outputs a value smaller than 1 at a non-sound source position, and outputs 1 at the sound source position, at the moment, the established DAS method not only can accurately position the sound source and quantize the average sound pressure contribution of the sound source, but also can meet the precondition that the FB method is established, in order to solve the problem, a PSF simulation is carried out by taking a 36-channel solid sphere array with the radius of 0.0975m shown in FIG. 1 as an example, the specific flow is that ① assumes a point sound source radiating sound waves with specific frequency at a specific position, ② calculates the sound field transfer function from the sound source to each microphone of the array according to the formula (1) and forms a column vector t (kr)₀,Ω₀) It should be noted that in the formula (1), the upper limit of the order N of the spherical harmonic is written to infinity, and considering that infinite terms cannot be calculated by simulation, 50 is used instead of infinity, ③ is used to set the focusing sound source plane, and the focusing components of the microphones are calculated for each focusing point by the formula (7) so as to form a column vector v (kr) when N is different from each other_f,Ω_f) ④ the PSF output value is calculated and plotted according to equation (12) where the focused sound source plane is set to be a spherical surface concentric with the array and having a radius equal to the distance from the sound source to the center of the array, the elevation angle is from 0 DEG to 180 DEG, the azimuth angle is from 0 DEG to 360 DEG, the intervals are all 5 DEG, and 37 × 73 focal points are counted, FIG. 2 shows that when N takes different values under 3000Hz conditions (r is a value of r₀＝1m,θ₀＝90°,) The PSF three-dimensional direction graph of the sound source at the position is characterized in that in each graph, three coordinate axes respectively represent components of the square root of the PSF in the x direction, the y direction and the z direction and are sequentially marked as PSF_x、psf_y、psf_zThe units are Pa, the coordinates of the points indicate the focus direction, the distance of the points to the origin is equal to the square root of the PSF in the focus direction, and the change in color also reflects the change in this square root. Psf notation_maxIs the maximum square root of the PSF, it is apparent that in each figure, it appears_xNegative, i.e. (theta)₀＝90°,) The true direction of the sound source of (c),the output values for the other directions are relatively small. When N is 1, psf is as shown in fig. 2(a)_max0.3111 only, much less than 1; further increase in N value, psf_maxAs shown in fig. 2(b) to (j), the psf corresponding to N-9 in fig. 2(i) increases with the increase in the number of bits_maxEqual to 0.9999, very close to 1, fig. 2(j) N10 corresponds to psf_maxExactly equal to 1; then continue to increase the value of N, psf_maxRemains unchanged and is always equal to 1, as shown in fig. 2(k) - (l). The above phenomena prove that: for the array, N is more than or equal to 10 under the condition of 3000Hz, and v (kr) can be ensured_f,Ω_f)≠t(kr₀,Ω₀) (when (r)_f,Ω_f)≠(r₀,Ω₀) Time) and v (kr)₀,Ω₀)＝t(kr₀,Ω₀) And meanwhile, the method is established and is a reasonable value. Further, the thick solid line in fig. 3 gives the radial function R of the acoustic signal propagation of the acoustic source_n(kr₀Ka) in which the symbol "◆" marks a point where n is 10, it is obvious that R takes any value greater than 10_n(kr₀Ka) is extremely low, almost 0, in combination with equation (1), which means that only spherical harmonic order pairs t (kr) of 10 or less are present₀,Ω₀) Contributes to the value of (c), which is v (kr) when N is more than or equal to 10₀,Ω₀)＝t(kr₀,Ω₀) The root cause of the establishment. Note N₀To make v (kr)₀,Ω₀)＝t(kr₀,Ω₀) Minimum established spherical harmonic order truncation length, accordingly, N > N₀When R is_n(kr₀Ka) is almost 0. Except for n, R_n(kr₀Ka) is also the wave number k or frequency f, the distance r of the sound source to the center of the array₀Array radius a as a function of test N₀Whether these factors are related to the values of (b) is determined, and FIG. 3 also shows R in 6 other cases_n(kr₀Ka) of the amplitude of the curve with n, ① f changing to 1000Hz and other parameters not changing, ② f changing to 5000Hz and other parameters not changing, ③ r₀To 0.5m, with other parameters unchanged ④ r₀2m for no change in other parameters, ⑤ a for 0.06m for no change in other parameters, ⑥ a for 0.2m for no change in other parametersAnd (6) changing. There is one N per curve₀Marked with the symbol "◆", N in the case of ③, ④₀Value of 10, N in the case of ①, ②, ⑤, ⑥₀The values are mutually different and are respectively 5, 14, 7 and 16, which shows that N is₀Is taken from the value of (a) and r₀Independently, it is only related to k and a. Analysis of R_n(kr₀The expression of ka), i.e., formula (2), can be found: r_n(kr₀Ka) is 0 is equivalent to j_n(ka) is equal toObtained N₀Necessarily only dependent on ka, proving that the simulation result is correct. To explore N₀Relationship with ka, changing the frequency of the sound source shown in FIG. 2, and counting N corresponding to different ka values₀The value obtained is shown by the symbol "◆" in FIG. 4, considering that N should be not less than N₀Value of (1) in fitting N₀The relation with ka should make the fitting value slightly larger than the actual value, and based on the principle, the obtained fitting relation is as follows:

wherein, the symbolMeaning rounding the value to the nearest integer toward positive infinity. The results of the fitting are shown as the solid line with "●" in FIG. 4. In practical application, the value of N is not less than N calculated by the formula (15)₀The value can ensure that the PSF outputs a value less than 1 at a non-sound source position and outputs 1 at a sound source position, so that the established DAS method not only can accurately position a quantitative sound source, but also can meet the precondition that an FB method is established.

Step 3, expanding solid sphere function beam forming

The FB method firstly decomposes a cross-spectrum matrix C of sound pressure signals received by an array microphone by using characteristic values to obtain:

C＝UΣU^H(16)

wherein U is [ U ]₁,u₂,…,u_q,…,u_Q]Is a unitary matrix, sigma ═ diag (sigma)₁,σ₂,…,σ_q,…,σ_Q) As a diagonal matrix, σ_q(Q is 1,2, …, Q) is the eigenvalue of matrix C, satisfying σ₁≤σ₂≤…≤σ_Q，u_qIs σ_qThe corresponding feature vector then introduces an exponential parameter ξ, let:

$C^{\frac{1}{\mathrm{\ξ}}}\≡Udiag({\mathrm{\σ}}_1^{\frac{1}{\mathrm{\ξ}}},{\mathrm{\σ}}_2^{\frac{1}{\mathrm{\ξ}}},...,{\mathrm{\σ}}_q^{\frac{1}{\mathrm{\ξ}}},...,{\mathrm{\σ}}_Q^{\frac{1}{\mathrm{\ξ}}})U^H=\underset{q=1}{\overset{Q}{\Σ}}{\mathrm{\σ}}_q^{\frac{1}{\mathrm{\ξ}}}{u}_q{u}_q^H---\left(17\right)$

finally, the existing sound source identification method is adopted for post-processingAnd calculating the ξ th power of the obtained result, expanding the FB method into the three-dimensional sound source identification of the solid sphere array according to the step 2, and constructing an output function of the FDAS method as follows:

$W_F({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=\frac{1}{Q}{\left(\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)C^{\frac{1}{\mathrm{\ξ}}}v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{\left|\right|v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)|{|}_2^2}\right)}^{\mathrm{\ξ}},(\mathrm{\ξ}\≥2)---\left(18\right)$

obviously, when ξ ═ 1, equation (18) is an expression of the DAS method output function in step 2. Assuming that only one monopole point sound source exists in the three-dimensional space, the average sound pressure contribution is

P_AC(kr₀,Ω₀) Theoretically, the following relationship holds:

σ₁＝σ₂＝…＝σ_Q-1＝0(19)

σ_Q＝tr(C)＝QP_AC(kr₀,Ω₀)(20)

$u_Q{u}_Q^H=\frac{t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)}{\left|\right|t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)|{|}_2^2}---\left(21\right)$

the combined type (17) to (21) can obtain:

$\begin{array}{c}{W}_F({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0){\left(\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{\left|\right|v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)\left|{|}_2^2\right|\left|t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\right|{|}_2^2}\right)}^{\mathrm{\ξ}}\\ =P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0){\mathrm{psf}}^{\mathrm{\ξ}}\left(\right({\mathrm{kr}}_f,{\mathrm{\Ω}}_f\left)\right|\left({\mathrm{kr}}_0,{\mathrm{\Ω}}_0\right))\end{array}---\left(22\right)$

the first situation is that the sound source position is at the focus point, the PSF output value at the sound source position is equal to 1 and the output value at the non-sound source position is less than 1 as long as the value of the spherical harmonic order truncation length N is large enough.By the nature of the exponential function, the arbitrary power of 1 remains equal to 1, and exponential functions with base numbers less than 1 decrease with increasing exponential parameters, so that, whatever the value ξ takes, the output value of the FDAS method at the sound source position will always equal the average sound pressure contribution of the sound source, while the output value at non-sound source positions will decrease with increasing ξ₀,Ω₀) Monopole point source with intensity s at position and (r)_f,Ω₀) Intensity at position r_fs/r₀The initial phase differs from the former by k (r)_f-r₀) Based on the fact that the sound pressure distributions generated by the monopole point sound sources on the surface of the spherical array are approximately the same, equation (9) can be written as follows:

$C=S({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\≈{\left(\frac{r_f}{r_0}\right)}^{2}S({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t({\mathrm{kr}}_f,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_0)---\left(23\right)$

formula (10) can be written as:

$P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)=\frac{1}{Q}S({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\left|\right|t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)|{|}_2^2\≈\frac{1}{Q}{\left(\frac{r_f}{r_0}\right)}^{2}S({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)|\left|t({\mathrm{kr}}_f,{\mathrm{\Ω}}_0)\right|{|}_2^2---\left(24\right)$

the joint type (23) and the formula (24) can obtain:

$C\≈\frac{P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)Q}{\left|\right|t({\mathrm{kr}}_f,{\mathrm{\Ω}}_0)|{|}_2^2}t({\mathrm{kr}}_f,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_0)---\left(25\right)$

substituting formula (25) for formula (6) to obtain:

W(kr_f,Ω_f)≈P_AC(kr₀,Ω₀)psf((kr_f,Ω_f)|(kr_f,Ω₀))(26)

it is clear that when the sound source position is not at the focus point, the main lobe peaks output by the FDAS method are all smaller than the average sound pressure contribution of the sound source, and the deviation increases with ξ, mainly because v (kr) is at the main lobe peak position when the sound source is not at the focus point_f,Ω_f)≠t(kr₀,Ω₀)，psf((kr_f,Ω_f)|(kr₀,Ω₀) To compensate for this deviation, a scale integration method is introduced, with the specific steps of ① building a scaling factor:

$\mathrm{\ϵ}=\frac{tr\left(C\right)}{{\mathrm{Q\Σ}}_{(r_f,{\mathrm{\Ω}}_f)\∈J}{W}_F({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}---\left(27\right)$

j is a set formed by a series of focus position coordinates, and the level of output values of the FDAS method at the focus positions is not less than the peak level of the main lobe and the set limit dynamic range D_rThe difference between, namely:

J＝{(r_f,Ω_f)|10lg|W_F(kr_f,Ω_f)|≥10lg|W_Fmax|-D_r}(28)

wherein the symbol "|" represents modulo, W_FmaxIs the main lobe peak. Introduction of D_rIn order to avoid side lobe integration, the strong side lobe attenuation capability of the FDAS method is taken into account, D_rThe larger value may be taken, which is set herein to 50dB ② scaling the output of the FDAS method to the scaled amount W_SF(kr_f,Ω_f)：

W_SF(kr_f,Ω_f)＝W_F(kr_f,Ω_f)(29)

③ integrating the scaling in the main lobe region to quantize the average sound pressure contribution of the sound source, and before integration, the level difference with the peak value is larger than D_rIs set to 0.

It should be noted that in practical applications, besides that the sound source does not fall on the focus point, the FDAS method suffers from a large quantization deviation due to factors such as mismatch of the frequency response of the microphone and unstable environmental conditions, and accordingly, the deviation compensation may also be performed by using a scaling integration method. Therefore, when the FDAS method identifies the sound source in a complex three-dimensional structure, the focusing sound source surface can be set as a simple spherical surface, so that the sound source can be accurately positioned, side lobes can be obviously attenuated, and the spatial resolution capability can be improved.

Drawings

FIG. 1 is a diagram of a ball array beamforming coordinate system of the present invention;

FIG. 2 shows the difference of N (r) at 3000Hz in the present invention₀＝1m,θ₀＝90°,) A PSF three-dimensional orientation map of the sound source at the location;

FIG. 3 shows the invention R_n(kr₀Ka) as a function of n.

FIG. 4 shows a truncation length N of the present invention₀Variation with ka.

Fig. 5 is a schematic view of a measurement arrangement of the present invention.

Fig. 6 is an imaging diagram of an incoherent sound source-recognition sound source of the present invention.

Detailed Description

The present invention is further illustrated by the following examples, but it should not be construed that the scope of the above-described subject matter is limited to the following examples. Various substitutions and alterations can be made without departing from the technical idea of the invention and the scope of the invention is covered by the present invention according to the common technical knowledge and the conventional means in the field.

The functional delay summation method for solid sphere array three-dimensional sound source identification comprises the following steps:

1): constructing a measuring system:

referring to fig. 1, an array of solid spheres in three-dimensional space; the solid sphere array comprises a solid sphere with the radius of a and Q microphones distributed on the surface of the solid sphere; the microphones are numbered Q, Q1, 2, …, Q; in an embodiment, Q-36; the sphere center of the solid sphere is used as an origin, the position coordinates of any point in the three-dimensional space where the solid sphere array is located are described by (r, omega), r represents the distance between the described position and the origin,denotes the direction of the depicted position, theta,Respectively its elevation angle and azimuth angle; (a, Ω)_q) Position coordinates of a microphone with the number q; as a verification, in this example, the frequency of the sound source was 3000Hz, and the position (r) was₀＝1m,θ₀＝90°,)。

2) Obtaining sound pressure

Each microphone receives a sound pressure signal p (ka, Ω)_q) (ii) a Wherein, wave number k is 2 pi f/c, sound wave frequency is f, sound velocity is c, (a, omega)_q) Position coordinates of a microphone with the number q;

3) the sound pressure signals are received by the respective microphones to obtain a cross-spectrum matrix C (which can be directly obtained by using pulse software of B & K corporation):

$C=\stackrel{\‾}{{\mathrm{pp}}^H},$

wherein: b represents a set of all sound source position coordinates, (r)₀,Ω₀) And (r)₀',Ω₀') each represent the position coordinates of the sound source,

$p={\mathrm{\Σ}}_{(r_0,{\mathrm{\Ω}}_0)\∈B}s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0),$

s(kr₀,Ω₀) Which is indicative of the intensity of the sound source,

t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),…,t_Q(kr₀,Ω₀)]^T，t_q(kr₀,Ω₀) Representing the sound field transfer function from the sound source to the q microphone,

$t_q({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\≡\underset{n=0}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-n}{\overset{n}{\Σ}}{R}_n({\mathrm{kr}}_0,ka)Y_n^{m*}\left({\mathrm{\Ω}}_0\right)Y_n^m\left({\mathrm{\Ω}}_q\right),$

are respectively omega₀、Ω_qSpherical harmonics of direction, n, m being orders of the spherical harmonics

The superscripts "-", "H" represent the average operation and the transposed conjugate operation, respectively;

4) determining the spherical harmonic order truncation length corresponding to each frequency:

5) substituting the cross-spectrum matrix C and the spherical harmonic order truncation length N into:

$W_F({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0){\left(\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{\left|\right|v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)\left|{|}_2^2\right|\left|t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\right|{|}_2^2}\right)}^{\mathrm{\ξ}}---\left(1\right)$

wherein: p_AC(kr₀,Ω₀) Is the trace of the cross-spectrum matrix C, (r)_f,Ω_f) As a position coordinate of the focal point, v (kr)_f,Ω_f)＝[v₁(kr_f,Ω_f),v₂(kr_f,Ω_f),…,v_Q(kr_f,Ω_f)]^T，v_q(kr_f,Ω_f) The focus component of the microphone of the q-th order is represented,ξ is an index parameter, the recommended value is 16, and the superscript "T" represents transposition operation;

6) according to W_F(kr_f,Ω_f) And drawing a sound source imaging graph, realizing sound field visualization and accurately positioning the sound source position.

Referring to FIG. 2, when N is different values (1-50) (r) under 3000Hz₀＝1m,θ₀＝90°,) A PSF three-dimensional orientation map of the sound source at the location;

in the novel FDAS sound source identification method established in this embodiment, when the value of the spherical harmonic order truncation length N of the PSF is large enough, the PSF converges, and outputs 1 at the sound source position and outputs a value smaller than 1 at the non-sound source position, which satisfies the establishment premise of function beam forming.

Referring to fig. 3, the curve of the amplitude of the radial function of the simulated acoustic signal propagation along with the order of the spherical harmonics shows that: the minimum value of N for which the PSF converges depends only on the wave number and the spherical array radius. Referring to fig. 4, the statistical value and the fitting value of the spherical harmonic truncation length are plotted as a function of the wave number and the product of the array radius, and the fitting relation is shown in formula (15).

Referring to fig. 5, a schematic view of the measurement arrangement of the present embodiment is shown.

Referring to fig. 6, it can be seen that the value principle of the method N of the present invention is correct, the sound source position can be accurately and clearly located, and the sound source average sound pressure contribution can be quantized, for coherent sound sources, the result is similar, and the overall performance is not inferior to that of the existing SHB method. The embodiment constructs the FDAS method, and successfully expands the function beam forming method for identifying the planar array two-dimensional sound source into the identification of the solid sphere array three-dimensional sound source.

Claims (1)

1. The functional delay summation method for solid sphere array three-dimensional sound source identification is characterized by comprising the following steps:

1) constructing a measuring system:

an array of solid spheres in three-dimensional space; the solid sphere array comprises a solid sphere with the radius of a and Q microphones distributed on the surface of the solid sphere; the microphones are numbered Q, Q1, 2, …, Q; the sphere center of the solid sphere is used as an origin, the position coordinates of any point in the three-dimensional space where the solid sphere array is located are described by (r, omega), and r represents the position between the described position and the originThe distance between the first and second electrodes,denotes the direction of the depicted position, theta,Respectively its elevation angle and azimuth angle; (a, Ω)_q) Position coordinates of a microphone with the number q;

2) obtaining sound pressure

Each microphone receives a sound pressure signal p (ka, Ω)_q) (ii) a Wherein, wave number k is 2 pi f/c, sound wave frequency is f, sound velocity is c, (a, omega)_q) Position coordinates of a microphone with the number q;

3) receiving sound pressure signals through the microphones to obtain a cross-spectrum matrix C:

$C=\stackrel{\‾}{{\mathrm{pp}}^H},$

wherein: b represents a set of all sound source position coordinates, (r)₀,Ω₀) And (r)₀',Ω₀') each represent the position coordinates of the sound source,

$p={\mathrm{\Σ}}_{(r_0,{\mathrm{\Ω}}_0)\∈B}s({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0),$

s(kr₀,Ω₀) Which is indicative of the intensity of the sound source,

t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),…,t_Q(kr₀,Ω₀)]^T，t_q(kr₀,Ω₀) Representing the sound field transfer function from the sound source to the q microphone,

$t_q({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\≡\underset{n=0}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-n}{\overset{n}{\Σ}}{R}_n({\mathrm{kr}}_0,ka)Y_n^{m*}\left({\mathrm{\Ω}}_0\right)Y_n^m\left({\mathrm{\Ω}}_q\right),$

are respectively omega₀、Ω_qThe spherical harmonics of the directions, n and m are the orders of the spherical harmonics,

the superscripts "-", "H" represent the average operation and the transposed conjugate operation, respectively;

4) determining the spherical harmonic order truncation length corresponding to each frequency:

indicating rounding the value to the nearest integer toward positive infinity;

5) substituting the cross-spectrum matrix C and the spherical harmonic order truncation length N into:

$W_F({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)=P_{AC}({\mathrm{kr}}_0,{\mathrm{\Ω}}_0){\left(\frac{v^H({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)t^H({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)}{\left|\right|v({\mathrm{kr}}_f,{\mathrm{\Ω}}_f)\left|{|}_2^2\right|\left|t({\mathrm{kr}}_0,{\mathrm{\Ω}}_0)\right|{|}_2^2}\right)}^{\mathrm{\ξ}}---\left(I\right)$

wherein: p_AC(kr₀,Ω₀) Is the trace of the cross-spectrum matrix C, (r)_f,Ω_f) As a position coordinate of the focal point, v (kr)_f,Ω_f)＝[v₁(kr_f,Ω_f),v₂(kr_f,Ω_f),…,v_Q(kr_f,Ω_f)]^T，v_q(kr_f,Ω_f) The focus component of the microphone of the q-th order is represented,ξ is an exponential parameter, the superscript "T" indicates the transpose operation;

6) according to W_F(kr_f,Ω_f) The method can locate the position of the sound source and draw a sound source imaging graph to realize three-dimensional sound field visualization.