CN117634144B - Rapid deconvolution sound source positioning algorithm based on beam focus neighborhood imaging evaluation - Google Patents
Rapid deconvolution sound source positioning algorithm based on beam focus neighborhood imaging evaluation Download PDFInfo
- Publication number
- CN117634144B CN117634144B CN202311409206.5A CN202311409206A CN117634144B CN 117634144 B CN117634144 B CN 117634144B CN 202311409206 A CN202311409206 A CN 202311409206A CN 117634144 B CN117634144 B CN 117634144B
- Authority
- CN
- China
- Prior art keywords
- sound source
- grid point
- grid
- damas
- point
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 85
- 238000003384 imaging method Methods 0.000 title claims abstract description 27
- 238000011156 evaluation Methods 0.000 title claims abstract description 12
- 239000011159 matrix material Substances 0.000 claims abstract description 24
- 230000003044 adaptive effect Effects 0.000 claims abstract description 10
- 238000001228 spectrum Methods 0.000 claims abstract description 6
- 238000005259 measurement Methods 0.000 claims abstract description 4
- 238000000034 method Methods 0.000 claims description 18
- 230000004807 localization Effects 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 9
- 238000001514 detection method Methods 0.000 claims description 9
- 230000003595 spectral effect Effects 0.000 claims description 2
- 239000000523 sample Substances 0.000 claims 1
- 241001077419 Damas Species 0.000 abstract description 29
- 230000006835 compression Effects 0.000 abstract description 4
- 238000007906 compression Methods 0.000 abstract description 4
- 238000004088 simulation Methods 0.000 description 23
- 238000010586 diagram Methods 0.000 description 13
- XOFYZVNMUHMLCC-ZPOLXVRWSA-N prednisone Chemical compound O=C1C=C[C@]2(C)[C@H]3C(=O)C[C@](C)([C@@](CC4)(O)C(=O)CO)[C@@H]4[C@@H]3CCC2=C1 XOFYZVNMUHMLCC-ZPOLXVRWSA-N 0.000 description 6
- 230000000694 effects Effects 0.000 description 5
- 238000012795 verification Methods 0.000 description 5
- 230000001133 acceleration Effects 0.000 description 2
- 230000001427 coherent effect Effects 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 230000001747 exhibiting effect Effects 0.000 description 1
- 238000002826 magnetic-activated cell sorting Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
- G01S5/18—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- General Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Operations Research (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Circuit For Audible Band Transducer (AREA)
Abstract
The invention discloses a rapid deconvolution sound source positioning algorithm based on beam focus neighborhood imaging evaluation, which relates to the technical field of sound source positioning and comprises the following steps: s1, establishing a signal measurement model; s2, establishing delay summation beam forming and point spread functions; s3, according to the cross spectrum matrix and the delay summation beam forming, solving a point spread function, and expressing the output of the delay summation beam forming as a linear equation set of all grid points; s4, obtaining effective grid points through an adaptive threshold algorithm; s5, solving a linear equation set of the effective grid points through deconvolution iteration, and finally finding out the sound source distribution of the scanning plane. The invention obtains the threshold value of each grid point by summing the area table for each grid point so as to judge whether the grid point is reserved or not, thereby realizing the compression of the grid point. A new point spread function matrix is calculated for the reserved grid points and is input into DAMAS algorithm to achieve the purpose of accelerating DAMAS algorithm.
Description
Technical Field
The invention relates to the technical field of sound source localization, in particular to a rapid deconvolution sound source localization algorithm based on beam focus neighborhood imaging evaluation.
Background
In recent years, acoustic photography techniques based on microphone arrays are increasingly exhibiting their indispensable roles in many fields of noise control, defect detection, aerospace, and the like. The core method of the acoustic photographing technology is to collect sound field information in the environment by using a microphone array and combine a sound source positioning algorithm so as to accurately determine the position of a sound source. Since acoustic photography generally requires real-time, the computation time of the algorithm is one of the key considerations in designing sound source localization algorithms. Existing acoustic photography techniques mostly use beamforming as a sound source localization algorithm. The delay-and-sum beamforming algorithm (delay and sum beamforming, DAS) is one of the most commonly applied beamforming algorithms in acoustic photography. The DAS has the advantages of simple implementation and high operation speed, but has the disadvantages of higher side lobe, low spatial resolution, narrow dynamic range and the like. Following the advent of the fast fourier transform, the DAS calculated in the time domain can also be implemented in the frequency domain based on cross-spectral matrix (CSM). However, the frequency domain DAS only improves some of the computational efficiency and flexibility, and does not improve the disadvantages of the DAS.
In order to obtain better sound source localization effect and improve the spatial resolution and dynamic range of imaging results, the academic community proposes a plurality of deconvolution sound source localization algorithms with better performance. In 2004, brooks proposed DAMAS algorithm, which removes the influence of channel propagation in beam forming by modeling a point spread function, improves spatial resolution, and can more truly reflect the spatial distribution of noise source intensity. Then, DAMAS-C algorithm applicable to multiple coherent noise sources was also proposed. DAMAS because of the significant effect in improving spatial resolution, this algorithm has become one of the dominant techniques of acoustic imaging technology, and is considered to be a major breakthrough in deconvolution sound source localization. Although the DAMAS algorithm performs well in terms of high spatial resolution and dynamic range and effectively attenuates side lobe effects, it results in a significant computational burden as it involves iteratively solving a system of linear equations. This makes algorithms highly demanding computational resources, and it is difficult to meet the real-time requirements of acoustic imaging.
To reduce the time complexity of DAMAS algorithms, researchers have made many searches. DAMAS2 and DAMAS use a fast fourier transform to transform the convolution between the sound source intensity distribution and the point spread function into the product of the beam domain, based on the assumption that the point spread function has spatial transfer invariant, improving the efficiency of the algorithm to some extent. But the spatial transfer invariance assumption does not hold when the scan plane size is much larger than the array aperture.
In 2008, yardibi introduces sparsity constraint into deconvolution sound source localization algorithm for the first time, solves through convex optimization, and proposes MACS beam forming algorithm applicable to coherent sound source based on sparsity in 2010. To further improve the algorithm efficiency, padois proposes an orthogonal matching pursuit method to solve, but is prone to falling into local optima. Bergh respectively introducing a non-negative matrix factorization method and a hierarchical clustering method, so that the speed of orthogonal matching pursuit calculation is ensured, and meanwhile, the global optimal solution can be obtained.
Another strategy for accelerating DAMAS algorithm, which is called a grid point compression algorithm, is to select a specific grid point according to the beam forming result, and form a linear equation set with smaller dimensions to perform iterative solution to improve efficiency. The DAMAS-CG1 algorithm performs grid point compression based on the concept of wavelet scaling, but the DAMAS-CG1 algorithm has a drawback that is difficult to avoid: when the sound field environment is complex, the probability of spatial aliasing increases greatly. In response to this problem, a DAMAS-CG2 algorithm was proposed, and the DAMAS-CG2 algorithm employed delay-and-sum beamforming with diagonal removal, with active grid points selected based on the non-negativity of the beamforming results. The DAMAS-CG2 algorithm not only avoids the occurrence of aliasing, but also accelerates the DAMAS algorithm. However, DAMAS-CG2 does not accelerate very well for the DAMAS algorithm. To further increase the speed of DAMAS-CG2 algorithm, DAMAS-CG3 algorithm introduces functional beamforming (functional beamforming, FB), which further increases algorithm speed through the advantage of the large dynamic range of FB algorithm. However, the DAMAS-CG3 algorithm may lead to reduced positioning accuracy in the case of array flow pattern errors due to, for example, array element position disturbances. Meanwhile, when the power parameter of DAMAS-CG3 is large, the grid point where the weaker sound source on the scan plane is located may be discarded.
Disclosure of Invention
The invention aims to provide a rapid deconvolution sound source positioning algorithm based on beam focus neighborhood imaging evaluation, which aims to solve the problem that the efficiency of the DAMAS improved algorithm is not enough to be applied to the acoustic photographing technology.
In order to achieve the above purpose, the present invention provides the following technical solutions:
a rapid deconvolution sound source localization algorithm based on beam focus neighborhood imaging evaluation comprises the following steps:
S1, establishing a signal measurement model: a planar microphone array comprising M microphones and a two-dimensional scanning plane comprising N grid points, said scanning plane being parallel to the plane in which the microphone array is located;
Establishing a calculation expression of a frequency domain form of a microphone array receiving signal at a detection frequency omega;
S2, establishing a response expression of delay summation beam forming aiming at an nth grid point on a scanning plane;
S3, according to the cross spectrum matrix and the delay summation beam forming, a point spread function is obtained, and based on the point spread function, the output of the delay summation beam forming is expressed as a linear equation set of all grid points;
S4, obtaining effective grid points through an adaptive threshold algorithm in DAMAS-AdaThr, and deforming the linear equation set in S3 into the linear equation set of the effective grid points;
S5, solving a linear equation set of the effective grid points in the S4 through deconvolution iteration, and finally finding the sound source distribution of the scanning plane.
Preferably, the calculation expression in S1 is:
Wherein,
Where ρ n (ω) is the component of the sound source intensity at the nth grid point at the detection frequency ω; τ nm is the time delay of the signal received by the mth microphone relative to the center of the microphone array when the nth grid point has a sound source; omega is the detection frequency; j is an imaginary unit.
Preferably, the expression in S2 is:
wherein c=xx H;
Where w n is the weight vector for delay-and-sum beamforming; matrix array A cross spectral matrix of signals received for the microphone array; For the conjugate transpose of w n, x H is the conjugate transpose of x.
Preferably, in S3, the point spread function is:
Where a n′ is the steering vector for grid point n', Is the conjugate transpose of a n′;
The linear equation set for all grid points is:
b＝Pq；
Where b is a vector of outputs of beam forming for all grid points, b= [ b 1,b2,...,bn]T, Q is a vector of acoustic power of the acoustic source at each grid point, q= [ q 1,q2,...,qn]T,
Preferably, in S4, the set of linear equations of the effective grid points is:
Wherein the method comprises the steps of N r is the number of active grid points reserved.
Preferably, in S4, the adaptive threshold algorithm includes the following steps:
S61, calculating the average value of all beam forming responses in a window with the size of r multiplied by r and taking each grid point as the center;
S62. if the beamforming response of the current grid point is greater than the threshold epsilon n of the grid point, the grid point is an active grid point, otherwise the grid point is regarded as an inactive grid point; the threshold is the product of this average value and the adjustment coefficient ζ:
Wherein, (i 1,j1) and (i 2,j2) are required to satisfy both of i 2-i1+1＝j2-j1+1＝r,(i1,j1) and (i 2,j2) as grid point coordinates.
Preferably, in S5, the calculation expression of the linear equation set for solving the effective grid point by deconvolution iteration is:
wherein i is a positive integer, AndRespectively are AndIs a combination of the elements of (1), Selection ofAs an initial value for the start of the iteration,Also beingBut n' noteqn.
Compared with the prior art, the invention has the following beneficial effects:
1. The invention obtains the threshold value of each grid point by summing the area table for each grid point so as to judge whether the grid point is reserved or not, thereby realizing the compression of the grid point. A new point spread function matrix is calculated for the reserved grid points and is input into DAMAS algorithm to achieve the purpose of accelerating DAMAS algorithm.
2. The invention uses the adaptive threshold algorithm to select the grid points containing effective information to reduce the dimension of the linear equation system based on the adaptive threshold algorithm, thereby realizing the purpose of acceleration. In terms of the running time, DAMAS-AdaThr provided by the invention can keep the running time within 10% of the DAMAS algorithm under the same scene. And in some complex scenes, the running time and the positioning accuracy are better than DAMAS-CG3 algorithm. Although the positioning accuracy is not as good as DAMAS-CG2 algorithm, the operation time of DAMAS-CG2 is far longer than DAMAS-AdaThr proposed by the invention.
Drawings
Fig. 1 is a schematic representation of the spatial location of a planar square microphone array and a two-dimensional scan of the present invention.
FIG. 2 (a) is a schematic diagram of the imaging results of delay-and-sum beamforming in a single-sound-source scene simulation experiment according to the present invention, and the simulated sound source placement positions are marked with an "X";
Fig. 2 (b) is a schematic diagram of imaging results of DAMAS algorithm in the single sound source scene simulation experiment, and the simulated sound source placement position is marked by "x";
FIG. 2 (c) is a schematic diagram of grid points (small circles in the drawing) reserved based on delay-sum beamforming results by using DAMAS-AdaThr in a single-sound-source scene simulation experiment, and simulated sound source placement positions are marked by using an 'X';
Fig. 2 (d) is a schematic diagram of an imaging result using DAMAS-AdaThr algorithm in the single sound source scene simulation experiment, and the simulated sound source placement position is marked by an "x".
FIG. 3 (a) is a schematic diagram of the imaging results of delay-and-sum beamforming in a three-source (different source intensities) scene simulation experiment of the present invention, where the simulated sound source placement locations are marked with an "x";
Fig. 3 (b) is a schematic diagram of imaging results of DAMAS algorithm in a single three sound source (different sound source intensities) scene simulation experiment, and the simulated sound source placement positions are marked by using an 'x';
FIG. 3 (c) is a schematic diagram of grid points (small circles in the drawing) reserved based on delay and sum beam forming results by DAMAS-AdaThr in a three-sound source (different sound source intensity) scene simulation experiment, and simulated sound source placement positions are marked by an X;
Fig. 3 (d) is a schematic diagram of imaging results using DAMAS-AdaThr algorithm in the three-sound source (different sound source intensities) scene simulation experiment, and the simulated sound source placement positions are marked by an 'x'.
FIG. 4 (a) is a schematic diagram of the imaging results of delay-and-sum beamforming in a "DAS" graphical scene simulation experiment of the present invention, with simulated sound source placement locations marked with "X";
FIG. 4 (b) is a schematic diagram of imaging results of DAMAS algorithm in "DAS" graphic scene simulation experiment of the present invention, where simulated sound source placement positions are marked with "X";
FIG. 4 (c) is a schematic diagram of grid points (small circles in the figure) reserved based on delay-and-sum beamforming results by DAMAS-AdaThr in a "DAS" graphical scene simulation experiment of the invention, and simulated sound source placement positions are marked by "X";
FIG. 4 (d) is a schematic diagram of imaging results using DAMAS-AdaThr algorithm in the "DAS" graphical scene simulation experiment of the present invention, and simulated sound source placement positions are marked with "X".
FIG. 5 is a schematic diagram of the variation relationship between the "DAS" graphic scene and the number of iterations of the present invention.
Detailed Description
The following description of the embodiments of the present invention will be made clearly and fully with reference to the accompanying drawings, in which it is evident that the embodiments described are only some, but not all embodiments of the invention; all other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
The invention provides a rapid deconvolution sound source localization algorithm based on beam focus neighborhood imaging evaluation, which comprises the following steps:
s1, establishing a signal measurement model;
A planar uniform microphone array comprising M microphones and a two-dimensional scan plane are shown in fig. 1. The scanning plane is parallel to the plane of the microphone array, and the distance between the two planes is z 0. A total of N grid points are arranged on the scanning plane, and each horizontal or vertical axis is provided with Each grid point is arranged on the shaft, and the distance between each adjacent grid point on the shaft is deltax. Thus, the size of the scan plane isA spatial reference point is selected as the origin of coordinates to establish a three-dimensional cartesian coordinate system. In this coordinate system, the coordinates of the mth microphone are:
pm＝[pmx,pmy,pmz]T，m＝1，...，M (1)；
Wherein T is a matrix transpose symbol;
the coordinates of the nth grid point on the scanning plane are:
gn＝[gnx,gny,gnz]T，n＝1，...，N (2)；
the center of the microphone array is selected as a spatial reference point, the plane in which the microphone lies is the x-y plane, and the z-axis passes through the centers of the microphone array and the scan plane and is perpendicular to both the plane in which the microphone array lies and the scan plane. The coordinates of each microphone satisfy:
Defining the distance from the nth grid point to the center of the microphone array as d n, and the distance from the nth grid point to the center of the mth microphone array as follows:
Where g nx is the x-axis coordinate of the nth grid point on the scan plane, g ny is the y-axis coordinate of the nth grid point on the scan plane, g nz is the z-axis coordinate of the nth grid point on the scan plane, p mx is the x-axis coordinate of the mth microphone, p my is the y-axis coordinate of the mth microphone, and p mz is the z-axis coordinate of the mth microphone.
Assuming that the nth grid point has a sound source, the delay τ nm of the signal received by the mth microphone with respect to the center of the microphone array is related to d n and d nm. The propagation speed of sound in air is c, and then the time delay τ nm is:
τnm＝(dnm-dn)/c (6)；
calculating a steering vector of the nth grid point:
Wherein j is an imaginary unit, and ω is a detection frequency;
Assuming that the waveform of the signal received by the mth microphone in the microphone array is s m (t), the spectrum of the signal received by the mth microphone can be obtained The frequency domain form of the microphone array receiving signal on the detection frequency is x (omega) = [ x 1(ω),x2(ω),...,xM(ω)]T ] after the vector is combined. The frequency domain version of the microphone array received signal at the probing frequency ω without any noise or interference being considered can be expressed as:
Where ρ n (ω) is the component of the sound source intensity at the nth grid point at the detection frequency ω.
S2, establishing a response expression of delay summation beam forming aiming at an nth grid point on a scanning plane;
the response of delay-and-sum beamforming for the nth grid point on the scan plane is:
Where w n＝an/M is the weight vector, matrix, of delay-and-sum beamforming A cross-spectral matrix (CSM) of signals received for the microphone array, the cross-spectral matrix expressed as:
C＝xxH (10)；
In order to eliminate self-noise for each channel of the microphone array, the diagonal elements of the cross-spectral matrix C may be removed when delay-and-sum beamforming is applied. Correspondingly, the weight vector of delay-sum beamforming applying diagonal removal
S3, according to the cross spectrum matrix and the delay summation beam forming, a point spread function is obtained, and based on the point spread function, the output of the delay summation beam forming is expressed as a linear equation set of all grid points;
The cross-spectral matrix expression can be obtained by taking the formula (8) into the formula (10):
taking equation (11) into equation (9) yields the response of the delay-and-sum beamforming for the nth grid point on the scan plane under ideal conditions:
Where N' =1,..n, Being the conjugate transpose of w n, a n′ is the steering vector for grid point n',Is the conjugate transpose of a n′,The acoustic power for the nth grid point.
From equation (12), it can be seen that the response b n of the beam formation for the nth grid point is not only related to the sound intensity ρ n of the nth grid point, but is obtained by multiplying the sound intensities of all grid points on the scan plane by a coefficient and finally summing. This coefficient is defined as the point spread function (point spread function, PSF). Thus, the point spread function is expressed as:
Calculating the values of the point spread function for all possible pairs of grid points (n, n') may form a point spread matrix The unfolding is as follows:
Based on the point spread function, the output of delay-and-sum beamforming can be expressed in the form of a system of linear equations:
b＝Pq (15)；
Where b is a vector of outputs of beam forming for all grid points, b= [ b 1,b2,...,bn]T, Q is a vector of the acoustic power of the acoustic source at each grid point, q= [ q 1,q2,...,qn]T,
Although delay-and-sum beamforming has the disadvantage that the position and number of sound sources cannot be accurately distinguished due to low spatial resolution, narrow dynamic range, etc., grid points without sound sources or with extremely low sound source intensity can be effectively identified according to the result of delay-and-sum beamforming, and the influence of the grid points on the result of beamforming by a point spread function is small or even negligible. Thus, the acoustic power q at these grid points can be assumed to be zero and considered as redundant grid points, without participating in the iterative solution process of the DAMAS algorithm to equation (15) thereafter. After the effective grid points are obtained through the effective grid point selection algorithm, selecting the values corresponding to the effective grid points from the vector b to form a new vectorSelecting elements corresponding to the effective grid point pairs from the matrix P to form a new matrixThe equation (15) can be modified accordingly:
In the method, in the process of the invention, N r is the number of active grid points reserved.
S4, obtaining effective grid points through an adaptive threshold algorithm in DAMAS-AdaThr, and deforming the linear equation set in S3 into the linear equation set of the effective grid points;
By means of the adaptive thresholding algorithm DAMAS-AdaThr, the average of all the beamformed responses within a window of size r x r centered at each grid point is first calculated. If the beamforming response of the current grid point is greater than the product of the average value and the adjustment coefficient, then the grid point is a valid grid point, otherwise the grid point is considered as an invalid grid point.
Calculating the average value within the sliding window of the beamformed response requires calculating a sum area table (summed area table, SAT) b s of the beamformed response b, with which b s the process of solving the average value can be accelerated. B s (i, j) of each element of the sum total area table b s is expressed as:
bs(i,j)＝b(i,j)+bs(i-1,j)+bs(i,j-1)-bs(i-1,j-1) (17)；
Where (i, j) is a grid point located in the ith row and jth column on the scanning plane.
After calculation of the sum area table b s, the beamformed response sum of any square region consisting of the upper left corner grid point (i 1,j1) and the lower right corner grid point (i 2,j2) can be calculated:
The threshold epsilon n for the nth grid point is recalculated:
Where, (for the adjustment coefficients, (i 1,j1) and (i 2,j2) need to satisfy i 2-i1+1＝j2-j1 +1=r. Only grid points n are reserved in the adaptive thresholding algorithm, which satisfy:
bn＞εn (20)；
S5, solving a linear equation set of the effective grid points in the S4 through deconvolution iteration, and finally finding the sound source distribution of the scanning plane.
In the formula (16), the values are calculated by the formulas (9) and (14)AndSolving by deconvolution sound source localization algorithmWherein the dot diffusion matrixUsually singular matrices, not by solving for matricesIs used to solve this system of linear equations. Furthermore, the matrixThe rank is typically much smaller than its dimension, soMany solutions are possible. Therefore, a Gaussian-Seidel iterative method is applied to solve the system of linear equations. In each round of iteration, a guaranteed vector is requiredAll elements of (2) are positive numbers. From an initial valueTo the ith iteration resultIs expressed as:
wherein i is a positive integer, AndRespectively are AndIs a combination of the elements of (1), General choicesAs an initial value for the start of the iteration.
By setting three simulated acoustic scenes, the invention respectively applies DAMAS, DAMAS-AdaThr, DAMAS-CG2 and DAMAS-CG3 for comparison verification. All simulation experiments were run on a computer loaded with AMD Ryzen 5800H 3.20GHz processor. The microphone array is a uniform circular array with an array aperture d=1m for a total of 65 microphones, with 64 microphones on the circle and one microphone in the center of the array. The microphone array is at a distance of 2m from the scan plane. The size of the scanning plane is 4m×4m, and is divided into 51×51 total 2601 grid points. The frequency of the simulated sound source used was 2000Hz. The propagation speed of sound in air was 340m/s. The signal of each microphone was added with white gaussian noise with a signal to noise ratio of 20dB. Typically, the sliding window should be sized slightly larger than the-3 dB main lobe width. The acoustic power of the simulated sound source is ρ 0. DAMAS-AdaThr use a sliding window size r=7. The DAS used by DAMAS-CG2 algorithm applies a diagonal removal algorithm of the cross-spectrum matrix.
The performance indexes of the evaluation positioning algorithm are respectively the relative calculation time taking the running time of DAMAS as a reference, the standard deviation sigma of grid points and three acoustic power errors, namely the total acoustic power error delta E o, the specific acoustic power error delta E s and the acoustic power inverse error delta E I. The grid-point-wise standard deviation σ of the result q (i) of the ith iteration of deconvolution relative to the true value q e is defined as:
Where q e,n is the true sound source power at the nth bin, and if there is no sound source at the bin, the value is 0.
Defining the sum of the acoustic powers of simulated sound sources placed on a scan plane in each simulated scene asDeconvolution iteration result acoustic power sum for each grid point And summing deconvolution results at grid points each located within a circle centered at the sound source placement position to q sum,specific. The reason for summing the deconvolution results at grid points within the circle is that the real sound source location is not necessarily on the scan grid. In the simulation experiment, the diameter of this circle was 0.1d.
The calculation of the total acoustic power error Δe o, the specific acoustic power error Δe s, and the acoustic power inverse error Δe I is shown by equations (23) to (25).
ΔEo＝|qsum,exact-qsum,deconv| (23)；
ΔEs＝|qsum,exact-qsum,specific| (24)；
ΔEI＝|qsum,deconv-qsum,specific| (25)；
The overall acoustic power error Δe o indicates the ability of the algorithm to quantify all of the sound sources across the scan plane, the specific acoustic power error Δe s indicates the ability of the algorithm to quantify each sound source, and the acoustic power inverse error Δe I indicates the ability of the algorithm to accurately locate the sound source.
Verification example 1
Single sound source: the coordinates of a simulated sound source on the scan plane are (0 m,0 m). The DAMAS-AdaThr adjustment coefficient ζ=1.6 of the present invention. The imaging results of DAS are shown in fig. 2 (a), DAMAS in fig. 2 (b), DAMAS-AdaThr with reserved grid points in fig. 2 (c), and DAMAS-AdaThr in fig. 2 (d). Simulation experiment results of the four algorithms for the single sound source are shown in table 1.
Table 1 single sound source simulation experiment result table
The DAMAS-AdaThr algorithm of the invention can accurately deconvolute the sound source distribution on the scan plane in the simplest acoustic scene. DAMAS is 3.72s and DAMAS-AdaThr is 4.8% of DAMAS, indicating DAMAS-AdaThr are effective in accelerating the deconvolution iterative solution process. And the delta E I of the four algorithms aiming at the simulation scene is 0, which indicates that the sound sources are accurately positioned.
Verification example 2
Three sound sources (different source intensities):
The coordinates of the three simulated sound sources are (0 m,0 m), (-0.28 m,0.28 m) and (0.28 m, -0.28 m), respectively. Wherein the power of the simulated sound source with coordinates (0 m,0 m) is ρ 0, and the acoustic power of the other two simulated sound sources is 0.25 ρ 0. DAMAS-AdaThr has an adjustment coefficient ζ=1.6. The imaging results of DAS are shown in fig. 3 (a), DAMAS in fig. 3 (b), DAMAS-AdaThr reserved grid points in fig. 3 (c), and DAMAS-AdaThr in fig. 3 (d). Simulation experiment results of the four algorithms for scene 2 are shown in table 2.
Table 2 three sound source simulation experiment result table
Simulation experiments prove that sound source distribution on a scanning plane can still be accurately deconvoluted and calculated in scenes of sound sources with different intensities by DAMAS-AdaThr algorithm, and sound sources with weaker intensities can not be ignored. In scenario 2, DAMAS runs 3.63s, DAMAS-AdaThr is the shortest of the four algorithms, 0.040s, 1.1% of DAMAS run. The delta E I of the four algorithms are all O, which indicates that sound sources with different intensities are accurately positioned.
Verification example 3
"DAS" graphics: 24 sound sources are placed on the scanning plane to form a DAS figure. DAMAS-AdaThr has an adjustment coefficient. The imaging results of DAS are shown in fig. 4 (a), DAMAS in fig. 4 (b), DAMAS-AdaThr reserved grid points in fig. 4 (c), and DAMAS-AdaThr in fig. 4 (d). Simulation experiment results of the four algorithms for scene 3 are shown in table 3.
TABLE 3 simulation experiment result table of DAS
DAMAS has an operation time of 3.54s and DAMAS-AdaThr has an operation time of 0.25s, which is the shortest of four algorithms, and is 7.2% of DAMAS. Note that the lattice point preservation idea of DAMAS-AdaThr is to preserve local peak points in the neighborhood, so more lattice points with higher DAS output power are discarded in this scenario, resulting in a lower power level for the deconvolution result than for the real power, and therefore a larger Δe o for the deconvolution result. In terms of positioning accuracy, DAMAS-CG3 and DAMAS-AdaThr with good acceleration effects are large in delta E I, which indicates that certain deviation occurs in positioning. Q sum,specific of DAMAS-CG3 is 19.0897 ρ 0, q sum,specific of DAMAS-AdaThr is 9.0441 ρ 0. The effect of source quantization error on the accuracy of the Δe I estimate positioning can be eliminated by calculating the ratio of Δe I to q sum,specific. The ratio of DAMAS-CG3 is 21.1% and the ratio of DAMAS-AdaThr is 16.1%, so that the positioning accuracy of DAMAS-AdaThr in this scene is slightly better than DAMAS-CG3.
In this verification example, the variation curve of σ with respect to the iteration number is shown in fig. 5, and when σ is not changed any more with the increase of the iteration number, it can be determined that the iteration reaches the convergence state. From the data of fig. 5, it can be found that: DAMAS-AdaThr reached the converging state on iteration 38, while DAMAS-CG3 required 82 iterations to converge. Therefore, in practical application, whether the iteration converges or not can be judged by judging whether sigma of the previous i rounds of iteration is equal or not in the iteration solving process. After judging that the iteration has converged, the iteration process can be terminated in advance, so that the algorithm can output an acoustic imaging result more quickly. The relative run times for simulation experiments using the new iteration termination scheme are shown in table 4. DAMAS-AdaThr has less run time after applying the early convergence method, which is only 0.11s, compared to the run time of DAMAS-CG3, which is 0.70s.
The present invention is not limited to the preferred embodiments, and the present invention is described above in any way, but is not limited to the preferred embodiments, and any person skilled in the art will appreciate that the present invention is not limited to the embodiments described above, while the above disclosure is directed to various equivalent embodiments, which are capable of being modified or varied in several ways, any simple modification, equivalent variation and variation of the above embodiments according to the technical principles of the present invention will still fall within the scope of the present invention.
Claims (4)
1. The rapid deconvolution sound source localization algorithm based on beam focus neighborhood imaging evaluation is characterized by comprising the following steps of:
S1, establishing a signal measurement model: a planar microphone array comprising M microphones and a two-dimensional scanning plane comprising N grid points, said scanning plane being parallel to the plane in which the microphone array is located;
establishing a computational expression of a frequency domain form of the microphone array received signal at the probe frequency omega:
Wherein,
Where ρ n (ω) is the component of the sound source intensity at the nth grid point at the detection frequency ω; τ nm is the time delay of the signal received by the mth microphone relative to the center of the microphone array when the nth grid point has a sound source; omega is the detection frequency; j is an imaginary unit;
S2 establishes a response expression for the nth grid point on the scan plane for delay-sum beamforming:
wherein c=xx H;
Where w n is the weight vector for delay-and-sum beamforming; matrix array A cross spectral matrix of signals received for the microphone array; For the conjugate transpose of w n, x H is the conjugate transpose of x;
S3, according to the cross spectrum matrix and the delay summation beam forming, a point spread function is obtained, and based on the point spread function, the output of the delay summation beam forming is expressed as a linear equation set of all grid points;
S4, obtaining effective grid points through an adaptive threshold algorithm in DAMAS-AdaThr, and deforming the linear equation set in S3 into the linear equation set of the effective grid points; the adaptive threshold algorithm comprises the following steps:
S61, calculating the average value of all beam forming responses in a window with the size of r multiplied by r and taking each grid point as the center; b s (i, j) of each element of the sum total area table b s is expressed as:
bs(i,j)＝b(i,j)+bs(i-1,j)+bs(i,j←1)-bs(i←1，j-1)；
Where (i, j) is a grid point located in the ith row and jth column on the scanning plane;
After calculation of the sum area table b s, the beamformed response sum of any square region consisting of the upper left corner grid point (i 1,j1) and the lower right corner grid point (i 2,j2) can be calculated:
S62. if the beamforming response of the current grid point is greater than the threshold epsilon n of the grid point, the grid point is an active grid point, otherwise the grid point is regarded as an inactive grid point; the threshold is the product of this average value and the adjustment coefficient ζ:
wherein, (i 1,j1) and (i 2,j2) are required to satisfy both of i 2-i1+1＝j2-j1+1＝r,(i1,j1) and (i 2,j2) as grid point coordinates;
S5, solving a linear equation set of the effective grid points in the S4 through deconvolution iteration, and finally finding the sound source distribution of the scanning plane.
2. The rapid deconvolution sound source localization algorithm based on beam focus neighborhood imaging evaluation of claim 1, wherein in S3 the point spread function is:
Where a n′ is the steering vector for grid point n', Is the conjugate transpose of a n′;
The linear equation set for all grid points is:
b＝Pq；
Where b is a vector of outputs of beam forming for all grid points, b= [ b 1,b2,...,bn]T, Q is a vector of acoustic power of the acoustic source at each grid point, q= [ q 1,q2,...,qn]T,
3. The rapid deconvolution sound source localization algorithm based on beam focus neighborhood imaging evaluation of claim 1, wherein in S4, the set of linear equations for the effective grid points is:
Wherein the method comprises the steps of N r is the number of active grid points reserved.
4. The rapid deconvolution sound source localization algorithm based on beam focus neighborhood imaging evaluation of claim 1, wherein in S5, a calculation expression of a system of linear equations for deconvolution iteration solving the effective grid points is:
wherein i is a positive integer, AndRespectively are AndIs a combination of the elements of (1), Selection ofAs an initial value for the start of the iteration,Also beingBut n' noteqn.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311409206.5A CN117634144B (en) | 2023-10-27 | 2023-10-27 | Rapid deconvolution sound source positioning algorithm based on beam focus neighborhood imaging evaluation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311409206.5A CN117634144B (en) | 2023-10-27 | 2023-10-27 | Rapid deconvolution sound source positioning algorithm based on beam focus neighborhood imaging evaluation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117634144A CN117634144A (en) | 2024-03-01 |
| CN117634144B true CN117634144B (en) | 2024-09-13 |
Family
ID=90036700
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311409206.5A Active CN117634144B (en) | 2023-10-27 | 2023-10-27 | Rapid deconvolution sound source positioning algorithm based on beam focus neighborhood imaging evaluation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117634144B (en) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111640157A (en) * | 2020-05-28 | 2020-09-08 | 华中科技大学 | Checkerboard corner detection method based on neural network and application thereof |
| CN116148770A (en) * | 2023-04-21 | 2023-05-23 | 湖南工商大学 | Sound source positioning method, device and system based on array signal processing |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8170234B2 (en) * | 2005-05-10 | 2012-05-01 | The United States of America by the Administrator of the National Aeronautics and Space Adminstration | Deconvolution methods and systems for the mapping of acoustic sources from phased microphone arrays |
| US10871548B2 (en) * | 2015-12-04 | 2020-12-22 | Fazecast, Inc. | Systems and methods for transient acoustic event detection, classification, and localization |
| CN112180329B (en) * | 2020-09-07 | 2023-04-11 | 黑龙江工程学院 | Automobile noise source acoustic imaging method based on array element random uniform distribution spherical array deconvolution beam forming |
| CN113219409B (en) * | 2021-04-15 | 2023-08-18 | 华南理工大学 | Acoustic imaging and multi-sound source positioning method based on focusing grid screening |
| CN115792921A (en) * | 2022-11-01 | 2023-03-14 | 天津大学 | Composite material structure collision source positioning method based on acoustic emission time difference approximation |
-
2023
- 2023-10-27 CN CN202311409206.5A patent/CN117634144B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111640157A (en) * | 2020-05-28 | 2020-09-08 | 华中科技大学 | Checkerboard corner detection method based on neural network and application thereof |
| CN116148770A (en) * | 2023-04-21 | 2023-05-23 | 湖南工商大学 | Sound source positioning method, device and system based on array signal processing |
Non-Patent Citations (1)
| Title |
|---|
| "混响背景下基于波束成形的声成像技术";庞凯元;《舰船科学技术》;20230331;第45卷(第6期);133-139 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117634144A (en) | 2024-03-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US10607358B2 (en) | Ear shape analysis method, ear shape analysis device, and ear shape model generation method | |
| CN106788653B (en) | Adaptive beam forming method based on covariance matrix reconstruction | |
| CN110531313B (en) | Near-field signal source positioning method based on deep neural network regression model | |
| CN112180329B (en) | Automobile noise source acoustic imaging method based on array element random uniform distribution spherical array deconvolution beam forming | |
| Chu et al. | Deconvolution for three-dimensional acoustic source identification based on spherical harmonics beamforming | |
| Yang et al. | Functional delay and sum beamforming for three-dimensional acoustic source identification with solid spherical arrays | |
| Tourbabin et al. | Theoretical framework for the optimization of microphone array configuration for humanoid robot audition | |
| CN107765221B (en) | Deconvolution sound source imaging method suitable for identifying coherent and incoherent sound sources | |
| CN109343003B (en) | Method for identifying sound source formed by fast iterative shrinking wave beams | |
| CN107390216A (en) | High speed super-resolution stationary point scan imaging method based on wave-number domain coherence factor | |
| Padois et al. | Time domain localization technique with sparsity constraint for imaging acoustic sources | |
| CN109884627B (en) | Short-range millimeter wave rapid three-dimensional imaging method of any linear array configuration | |
| Xia et al. | Noise reduction method for acoustic sensor arrays in underwater noise | |
| CN113593596B (en) | Robust self-adaptive beam forming directional pickup method based on subarray division | |
| CN117634144B (en) | Rapid deconvolution sound source positioning algorithm based on beam focus neighborhood imaging evaluation | |
| CN111487594A (en) | Circular array beam forming method based on particle swarm optimization | |
| CN108614235B (en) | Single-snapshot direction finding method for information interaction of multiple pigeon groups | |
| Yan et al. | Low-degree root-MUSIC algorithm for fast DOA estimation based on variable substitution technique | |
| CN110736976A (en) | sonar beam former performance estimation method of arbitrary array | |
| Merino-Martınez et al. | Three–dimensional acoustic imaging using asynchronous microphone array measurements | |
| CN106908754A (en) | L-type acoustic vector-sensor array row ESPRIT decorrelation LMS method for parameter estimation | |
| CN116301195B (en) | Function beam optimization method and device | |
| CN111830465B (en) | Two-dimensional Newton orthogonal matching pursuit compressed beam forming method | |
| CN111551892B (en) | Steady self-adaptive beam forming method and device | |
| Chu et al. | Two-dimensional total variation norm constrained deconvolution beamforming algorithm for acoustic source identification |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |