CN107153172B - Cross-spectrum generalized inverse beam forming method based on cross-spectrum optimization - Google Patents
Cross-spectrum generalized inverse beam forming method based on cross-spectrum optimization Download PDFInfo
- Publication number
- CN107153172B CN107153172B CN201710315666.XA CN201710315666A CN107153172B CN 107153172 B CN107153172 B CN 107153172B CN 201710315666 A CN201710315666 A CN 201710315666A CN 107153172 B CN107153172 B CN 107153172B
- Authority
- CN
- China
- Prior art keywords
- matrix
- spectrum
- cross
- vector
- sound
- 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.)
- Expired - Fee Related
Links
- 238000001228 spectrum Methods 0.000 title claims abstract description 84
- 238000000034 method Methods 0.000 title claims abstract description 58
- 238000005457 optimization Methods 0.000 title claims abstract description 27
- 239000011159 matrix material Substances 0.000 claims abstract description 125
- 230000008901 benefit Effects 0.000 claims abstract description 10
- 239000013598 vector Substances 0.000 claims description 55
- 230000008569 process Effects 0.000 claims description 29
- 238000009826 distribution Methods 0.000 claims description 14
- 238000005259 measurement Methods 0.000 claims description 13
- 238000012546 transfer Methods 0.000 claims description 13
- 238000012805 post-processing Methods 0.000 claims description 11
- 238000001914 filtration Methods 0.000 claims description 9
- 230000006870 function Effects 0.000 claims description 9
- 230000004044 response Effects 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000010606 normalization Methods 0.000 claims description 4
- 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 claims description 3
- 238000004458 analytical method Methods 0.000 claims description 3
- 230000007423 decrease Effects 0.000 claims description 3
- 230000005284 excitation Effects 0.000 claims description 3
- 230000005404 monopole Effects 0.000 claims description 3
- 230000005540 biological transmission Effects 0.000 claims 1
- 238000001093 holography Methods 0.000 description 6
- 238000003491 array Methods 0.000 description 4
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 description 2
- 125000004432 carbon atom Chemical group C* 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000001788 irregular Effects 0.000 description 2
- 230000004807 localization Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000004075 alteration Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000002955 isolation Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 230000000452 restraining effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Images
Classifications
-
- 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
- G01S5/20—Position of source determined by a plurality of spaced direction-finders
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Circuit For Audible Band Transducer (AREA)
Abstract
The invention discloses a cross-spectrum generalized inverse beam forming method based on cross-spectrum optimization. And secondly, optimizing a sound pressure cross-spectrum matrix by combining a linear optimization method. And finally, constructing a new cross-spectrum generalized inverse beam forming algorithm based on the optimized cross-spectrum matrix by solving the obtained steering matrix and the optimized sound pressure cross-spectrum matrix and combining the traditional cross-spectrum matrix beam forming idea. The algorithm can fully utilize the anti-noise performance of the cross-spectrum beam forming based on the cross-spectrum optimization while applying the advantage of high resolution formed by the generalized inverse beam forming.
Description
Technical Field
The invention relates to a sound source identification algorithm, in particular to a method for accurately identifying and positioning a sound source on the basis of eliminating background noise by applying a cross-spectrum optimization generalized inverse beam forming algorithm based on a microphone array.
Background
At present, there are two main sound source identification methods based on microphone arrays, namely beam forming and acoustic holography. Both have advantages and disadvantages, wherein the identification performance of the beam-formed medium-long distance sound source is better than that of the acoustic holography, and the low-frequency performance is poorer than that of the acoustic holography.
Beamforming is a sound source identification method that has been rapidly developed in recent years, and it is based on regular or irregular microphone arrays (such as cross-axis and wheel-type arrays) to perform signal acquisition and then perform sound field reconstruction using, for example, a simple delay-and-sum algorithm. Compared with the traditional methods such as a subjective evaluation method, a selective isolation method and the like, the beam forming based on the multi-channel signal processing technology has the advantages of convenience in measurement, visual visualization of a reconstructed sound field and the like. Compared with a sound intensity measurement method, beam forming has the advantages of high measurement efficiency, wide application environment and the like. Compared with the acoustic holography method developed in recent years, the beam forming has higher spatial resolution at medium and high frequencies, can realize the reconstruction of a sound source at medium and long distances, and can adapt to different operating conditions. The noise source identification method based on beam forming has the advantages of flexible dynamic display range, strong interference capability, high medium-high frequency and high resolution, and the like, so the method is highly concerned by researchers and scholars in various countries in recent years. The application range of the system is continuously expanded, and the system is gradually developed to civil fields such as vehicle noise test from military fields such as aviation and sonar. Meanwhile, various beam forming algorithms are derived, such as traditional cross-spectrum beam forming, deconvolution algorithm and generalized inverse beam forming algorithm, and various continuously-developed beam forming algorithms of various colors have characteristics and are mutually connected, so that the application convenience and the practicability of the method are continuously promoted.
In the prior art, a sound source identification method based on a microphone array exists, and a broadband acoustic holography algorithm is disclosed on the basis of a traditional equivalent source algorithm. A series of equivalent sources are constructed on an equivalent source plane close to a sound source plane, and the equivalent sources are used for replacing the actual sound source distribution. And establishing an acoustic transfer equation from an equivalent source to the microphone array by a sound source identification theory. And solving the transfer equation through a traditional gradient descent algorithm so as to iteratively solve the equivalent source strength. From this equivalent source intensity, the sound source distribution between the equivalent source plane and the array plane is then reconstructed based on the sound propagation process. In addition, an iterative filtering process is introduced in the equivalent source solving process, so that the convergence efficiency and accuracy of the iterative process are ensured. Due to the introduction of a gradient descent algorithm and a filtering process, the resolution of the traditional equivalent source method at medium and low frequencies can be improved. Meanwhile, the algorithm can be combined with or switched over with other acoustic holography algorithms, so that the reconstruction of the medium-low frequency sound source is realized. But the algorithm is based on the traditional equivalent source algorithm reconstruction, the spatial resolution of the algorithm is higher in medium and low frequency sound sources and closer measurement distances, but the spatial resolution of sound source identification is lower for high frequency sound sources or longer distance measurement.
Disclosure of Invention
The invention aims to provide a sound source identification algorithm with high sound source reconstruction accuracy and good robustness, which is constructed on the basis of general generalized inverse and cross-spectrum beam forming, has high reconstruction accuracy, is fully combined with the optimization of a sound pressure cross-spectrum matrix, and has strong anti-noise performance.
The technical proposal adopted for realizing the aim of the invention is that the method for forming the cross-spectrum generalized inverse beam based on the cross-spectrum optimization is characterized in that,
building a sound source reconstruction system:
the system comprises a camera, a microphone array, a multi-channel signal collector and a post-processing computer, wherein the camera, the microphone array, the multi-channel signal collector and the post-processing computer are positioned near a sound source. The camera is used for collecting a space image of the whole sound field; the microphone array is used for measuring sound field signals; the multi-channel signal collector is used for measuring a time domain analog signal, converting the time domain analog signal into a digital signal and transmitting the digital signal to a computer for post-processing; the computer mainly has the main functions of storing and post-processing the collected sound field signals and efficiently displaying the distribution condition of the reconstructed sound field;
the beamforming process comprises the steps of:
1) computing a steering matrix based on generalized inverse beamforming:
with N unit sound sources arranged on the sound source plane, the measured sound pressure signals at the microphone ri points in the microphone array are as follows:
in the above formula, the first and second carbon atoms are,
is r in the sound fieldiAt a unit excitation of r0Complex response of (a) q (r)i) Is the assumed monopole sound source intensity. The equation (1) expresses the transfer relationship between the sound source and the microphone array, and can be converted into a vector-matrix form for the subsequent analysis:
p＝Aq (2)
where A is an NxM sound field transfer matrix whose element composition is determined by equation (1). p is the M-dimensional measured sound pressure vector whose elements are the frequency domain response of the corresponding sensor at a single frequency f. q is an N-dimensional sound source intensity column vector whose element components represent intensity distributions at corresponding points in a sound source plane.
To solve the reconstructed sound field distribution, equation (2) can be transformed into the minimization problem:
solving the minimization problem (3) by an iterative algorithm, and finally obtaining an expression of the steering matrix in the step (N):
w(n)＝L(n)B(n)(4)
in the formula (4) B(n)For an N × M matrix, the row vector can be expressed as follows:
B(i,:)(n+1)＝G(:,i)HALF/(G(:,i)HALF G(:,i)) (5)
the main function of equation (5) is to normalize the transfer matrix, thereby ensuring the accuracy of the output. In the above formula, the superscript H represents conjugate transpose, G (: i) is the ith column vector of the transfer matrix G after normalization, wherein the expression of G is as follows:
G＝AL(n)(6)
the ALF in equation (5) is an invertible M × M dimensional matrix, and its expression has the following form:
ALF＝(GGH+βI)-1(7)
the regularization parameter β in the above equation (7) can be obtained by an L curve or a GVC method.
L(n)＝diag(|Sp|/||Sp||∞) (8)
In the above formula (8), diag is to convert the vector into a diagonal matrix form, | | | is to take the absolute value of the vector, | | | the non-woven phosphor∞To find the infinite norm of the vector. Sp is generalized inverse beamforming after normalizationAnd (3) outputting:
Sp＝w(n)p (9)
the following filtering process is also added in the solving process for solving the beam output Sp:
in the formula (10), (| Sp |)maxTaking the largest element in the vector | Sp | where Dk increases with increasing number of iterations (n), such that the filtering threshold gradually decreases with increasing number of iterations:
Dk＝0.2+0.5n (11)
solving to obtain a steering matrix w based on generalized inverse beam forming(n)；
2) Optimizing a sound pressure cross-spectrum matrix:
based on an M-dimensional sound pressure vector p obtained by sensor array measurement, an M multiplied by M-dimensional sound pressure cross spectrum matrix is obtained, and the expression is as follows:
C＝ppH(12)
the optimized sound pressure cross-spectrum matrix form is expressed as follows:
Csuperior food＝C-diag(x) (13)
Wherein, CSuperior foodTo eliminate the matrix of diagonal elements, diag (x) is a pair of diagonal matrices, x is the noise component;
the element x (i,1) where x is defined as [0, C (i, i)]This maintains some of the advantages of the diagonalization element while reducing channel self-noise. To obtain the noise component x, the sound pressure cross-spectral matrix may be optimized by an iterative optimization algorithm. Assume a series of normalized unit vectors uk＝ek(k ═ 1, 2.., m), the following inequality holds:
uk H(C-diag(x))uk≥0 (14)
from the above equation, the following linear inequality can be derived:
wherein D (k, i) ═ u (i, k) <' > ceiling2(where u (i, k) is the vector ukElement No. i), d (k,1) ═ uk HCuk. A constraint Dx ≦ d for the noise component x may be obtained from the above equation, and since the core idea of the optimization of the sound pressure cross-spectral matrix is to minimize the noise component-x, an objective function min F ═ x may be obtained. Summarizing the above process of minimizing the noise component as follows, a linear programming problem can be achieved:
solving the minimized noise component-x and outputting the optimized sound pressure cross-spectrum matrix CSuperior food。
3) Reconstructing a beam output according to the steering matrix and the optimized sound pressure cross-spectrum matrix: the beam output is in the form:
Ssuperior food＝w(n)CSuperior food(w(n))H(18)
Sound pressure cross spectrum matrix C after optimization in the above formulaSuperior foodMay alternatively be in the form of:
Ssuperior food＝w(n)(C-diag(x))(w(n))H(19)
Taking diagonal elements from the above formula (19) and squaring to finally obtain the output of the cross-spectrum generalized inverse beam forming based on the optimized cross-spectrum matrix as follows:
Soutput of＝(diag(SSuperior food))1/2(20)。
The technical effect of the present invention is undoubtedly that the steering matrix is first constructed based on conventional generalized inverse beamforming. And secondly, optimizing a sound pressure cross-spectrum matrix by combining a linear optimization method. And finally, constructing a new cross-spectrum generalized inverse beam forming algorithm based on the optimized cross-spectrum matrix by solving the obtained steering matrix and the optimized sound pressure cross-spectrum matrix and combining the traditional cross-spectrum matrix beam forming idea. The algorithm can fully utilize the anti-noise performance of the cross-spectrum beam forming based on the cross-spectrum optimization while applying the advantage of high resolution formed by the generalized inverse beam forming.
Drawings
FIG. 1 is a schematic view of sound source reconstruction;
fig. 2 is a schematic diagram of a sound pressure cross-spectrum matrix optimization process.
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.
Since the present embodiment relates to a beam forming algorithm based on sound source identification, which is based on a sound source identification system, the sound source identification system will be briefly described first. Wherein figure 1 depicts a schematic diagram of a sound source reconstruction system comprising a camera positioned near the sound source, a microphone array, a multi-channel signal collector, and a post-processing computer.
The camera is used for collecting a space image of the whole sound field; and the beam output calculated by the beam forming algorithm is coincided with the space influence, so that the sound source distribution condition of the whole sound field can be obtained. The camera can be a common digital camera and can be transmitted to a computer through a common USB interface.
The microphone array is used for measuring sound field signals; the structure of the microphone comprises a sound pressure or sound velocity sensor, a microphone bracket and an associated wiring harness. In the microphone array, the distribution of the sensors may be a regular or an irregular one-, two-, three-dimensional spatial distribution, which is selected according to the actual measurement requirements. The choice of the number of sensors and the choice of the microphone spacing is then determined according to the frequency range to be measured.
The multi-channel signal collector is used for measuring a time domain analog signal, converting the time domain analog signal into a digital signal and transmitting the digital signal to a computer for post-processing; it should include a multi-channel preamplifier and a multi-channel analog/digital converter with high accuracy, and even associated filters, to ensure that the digital signal output to the computer is sufficiently accurate.
The computer mainly has the main functions of storing and post-processing the collected sound field signals and efficiently displaying the distribution condition of the reconstructed sound field; . The memory required by the computer is more than 512MB, 2G is recommended, and the CPU dominant frequency is more than 2.0G, so that the high efficiency of calculation can be fully ensured.
The beamforming process comprises the steps of:
supposing that the sound source reconstruction system comprises M sound pressure microphones, time domain signals collected by the multi-channel signal collector are transmitted and stored in a computer, and the collected time domain signals are converted into a frequency domain through fast Fourier transform. In the embodiment, the beam forming algorithm mainly relates to the reconstruction problem of the medium-high frequency sound source, so that a complex signal corresponding to a single frequency f component in an obtained array frequency domain signal is extracted and converted into a measurement sound pressure signal p to be processed by a later algorithm.
Since the algorithm is based on generalized inverse beamforming, cross-spectral matrix optimization and cross-spectral beamforming.
1) Computing a steering matrix based on generalized inverse beamforming:
wherein the steering matrix is a steering matrix after normalization, which basically acts like steering vectors in conventional cross-spectrum beam forming. To solve the normalized steering matrix, N unit sound sources are arranged on the sound source plane, and the measured sound pressure signals (frequency domain single frequency components) at the microphone ri points in the microphone array are as follows:
in the above formula, the first and second carbon atoms are,
is r in the sound fieldiAt a unit excitation of r0Complex response of (a) q (r)i) Is the assumed monopole sound source intensity. The expression (1) is from sound source to microphoneThe transfer relationship between the arrays, necessary for subsequent analysis, can be converted into a vector-matrix form:
p＝Aq (2)
where A is an NxM sound field transfer matrix whose element composition is determined by equation (1). p is the M-dimensional measured sound pressure vector whose elements are the frequency domain response of the corresponding sensor at a single frequency f. q is an N-dimensional sound source intensity column vector whose element components represent intensity distributions at corresponding points in a sound source plane.
To solve the reconstructed sound field distribution, equation (2) can be transformed into the minimization problem:
solving the minimization problem (3) by an iterative algorithm, and finally obtaining an expression of the steering matrix in the step (N):
w(n)＝L(n)B(n)(4)
in the formula (4) B(n)For an N × M matrix, the row vector can be expressed as follows:
B(i,:)(n+1)＝G(:,i)HALF/(G(:,i)HALF G(:,i)) (5)
the main function of equation (5) is to normalize the transfer matrix, thereby ensuring the accuracy of the output. In the above formula, the superscript H represents conjugate transpose, G (: i) is the ith column vector of the transfer matrix G after normalization, wherein the expression of G is as follows:
G＝AL(n)(6)
the ALF in equation (5) is an invertible M × M dimensional matrix, and its expression has the following form:
ALF＝(GGH+βI)-1(7)
the regularization parameter β in the above formula (7) can be obtained by an L curve or a GVC method, and the regularization parameter β recommended to be used in the patent is GGHOf maximum eigenvalue
Multiple (where f is the focused sound source frequency).
Regularization matrix L in equation (6)(n)Is an important parameter, and changes along with the change of the output of the generalized inverse beam forming in the iteration process, thereby ensuring the accuracy of the reconstruction process. It can be generally defined as follows:
L(n)＝diag(|Sp|/||Sp||∞) (8)
in the above formula (8), diag is to convert the vector into a diagonal matrix form, | | | is to take the absolute value of the vector, | | | the non-woven phosphor∞To find the infinite norm of the vector. Sp is the normalized generalized inverse beamforming output:
Sp＝w(n)p (9)
in the conventional generalized inverse beamforming output, interference side lobes will occur in the non-acoustic source direction due to the influence of noise. In order to further suppress the influence of these interference side lobes and thus achieve more accurate sound source localization, the following filtering process is also added in the solving process for solving the beam output Sp:
the core idea of the filtering process is that elements smaller than a threshold value in a beam output vector are continuously forced to be zero in the iteration process by assuming the threshold value changing along with the iteration, and then the influence of noise is inhibited. In the formula (10) (| Sp |)maxTaking the largest element in the vector | Sp | where Dk increases with increasing number of iterations (n), such that the filtering threshold gradually decreases with increasing number of iterations:
Dk＝0.2+0.5n (11)
it should be noted that the final purpose of step 1) is not to solve the beam output, but to solve for a steering matrix w based on generalized inverse beamforming(n)Namely, the steering matrix of the formula (2) is obtained through solving.
2) Optimizing a sound pressure cross-spectrum matrix:
based on an M-dimensional sound pressure vector p obtained by sensor array measurement, an M multiplied by M-dimensional sound pressure cross spectrum matrix is obtained, and the expression is as follows:
C＝ppH(12)
traditional beam forming based on sound pressure cross-spectrum matrix in order to remove the influence of channel self-noise, usually by eliminating the diagonal elements of the cross-spectrum matrix, i.e. making C ═ C- (C)Opposite angle. Since the diagonal elements in the sound pressure cross-spectrum matrix are the power spectrum of the same channel, and the power spectrum of the signal in the same channel may introduce channel self-noise, the anti-noise capability of the beam output can be improved by removing the diagonal elements. However, this method also has major side effects, such as negative power in the output spectrum of the beam, and even large deviation of sound source location. Of course, in practical measurement, the influence of negative situations such as positioning deviation can be reduced while the noise-resistant advantage of the diagonalization removing method is kept. For the situation, the sound pressure cross-spectrum matrix is optimized and improved based on the generalized inverse beam forming output of the cross-spectrum matrix.
If active components and noise components are included in diagonal elements in the sound pressure cross-spectrum matrix, and the noise components can be removed on the basis of storing the active components in the calculation process of the sound pressure cross-spectrum matrix, the sound field reconstruction can be performed by fully using the optimized sound pressure cross-spectrum matrix. Based on the above assumptions, we can express the optimized sound pressure cross-spectrum matrix form as follows:
Csuperior food＝C-diag(x) (13)
The noise component can be removed by the optimization method while the active component is maintained. The above equation is similar to the sound pressure cross-spectrum matrix with diagonal elements removed, where C is the matrix with diagonal elements removed, diag (x) is a diagonal matrix, which represents the noise component and is the portion that needs to be removed in the diagonal elements, which is the element whose diagonal elements correspond to the vector x, for example
The element x (i,1) where x is defined as [0, C (i, i)]In between, C (i, i) refers to the ith diagonal element of the matrix C, which maintains the diagonals of the matrix CSome of the advantages of the elements are combined with the ability to reduce channel self-noise. To obtain the noise component x, we can optimize the sound pressure cross-spectrum matrix by an iterative optimization algorithm. Assume a series of normalized unit vectors
uk＝ek(k ═ 1, 2.., m), the following inequality holds:
uk H(C-diag(x))uk≥0 (14)
from the above equation, the following linear inequality can be derived:
wherein D (k, i) ═ u (i, k) <' > ceiling2(where u (i, k) is the vector ukElement No. i), d (k,1) ═ uk HCuk. A constraint Dx ≦ D for the noise component x is obtained by the above equation, (D is a matrix whose element values are represented by D (k, i) ═ u (i, k) < >2Determining; d is a vector with the element value d (k,1) ═ uk HCukDetermining; ) Since the core idea of the optimization of the sound pressure cross-spectrum matrix is to minimize the noise component-x, an objective function min F ═ x can be obtained. Summarizing the above process of minimizing the noise component as follows, a linear programming problem can be achieved:
of course, the method for solving the linear programming problem (17) is many, and different optimization methods can be selected according to different requirements. Finally, solving the minimized noise component-x by selecting a general linear iterative algorithm, and outputting an optimized sound pressure cross-spectrum matrix CSuperior food。
The optimization process of the sound pressure cross-spectrum matrix is summarized as shown in fig. 2, and mainly comprises four processes, namely initialization, a first cycle, a second cycle and a third cycle. The initialization process that is first initiated is primarily to assign values to some variables. For optimized sound pressure cross-spectrum matrix CSuperior foodOfThe starting value we assume as a sound pressure cross-spectral matrix without any processing and let u be a matrix consisting of a series of normalized unit column vectors. Where b is the magnitude factor of the vector x, xLower partAnd xOn the upper partRespectively representing the lower and upper limit values of the vector x in the iterative solution process, mainly for restricting the convergence range of the vector x and ensuring the accuracy of the solution, wherein t1 and t2 are constraint variables, and mainly for improving the convergence speed of the optimization process, wherein t1 is more than or equal to 0<t2 ≦ 1, most commonly t1 may be made 0 and t2 may be set to 1. By setting the values of t1 and t2, it is possible to effectively switch between optimizing the sound pressure cross-spectral matrix and the method of diagonalization. The iteration step number N1 is a very important parameter, and its value varies with different measurement environments, and theoretically, the larger the value is, the better the value is, but the larger the value is, the lower the calculation efficiency is, so that an acceptable reasonable value must be selected in actual measurement.
Wherein the second and third loops of the optimization process are mainly for calculating the relevant parameters in the constraints, wherein the third loop calculates the matrix D by u. The second loop is mainly to compute the vector d by the matrix u and the pressure cross-spectrum matrix C. The most important of them belongs to the first loop, the beginning part of the first loop is mainly to calculate the loop times of the second loop, and to reassign the dimension and initial value of the matrix D and the vector D according to the dimension of the column of the matrix u. When the second and third loops are completed, a new matrix D and vector D have been calculated. At this time, the second half of the first cycle is optimized and solved by a conventional linear programming method on the basis of restraining the upper and lower limits of the noise component x, so as to calculate a new noise component x. And then converting the noise component vector x into a diagonal matrix, and subtracting the diagonal matrix from the sound pressure cross-spectrum matrix to obtain an optimized sound pressure cross-spectrum matrix, which is equivalent to subtracting the noise component from the diagonal elements of the sound pressure cross-spectrum matrix, thereby improving the anti-noise performance of the beam output result. Of course, an important step is added at the end of the first cycle, namely, the column vector of the matrix u obtained in the previous step and the feature vector of the optimized sound pressure cross-spectrum matrix are merged, and the obtained union is assigned to the new matrix u again, so that the number of columns of the matrix u is gradually increased along with the increase of the iteration times.
After the circulation process is completed, the optimized sound pressure cross-spectrum matrix C can be outputSuperior food. The optimized sound pressure cross-spectrum matrix can be applied to beam forming based on cross-spectrum, so that the output accuracy is ensured, and in the following step 3), the output is performed through the cross-spectrum beam forming based on the generalized inverse, so that sound source positioning identification is performed.
3) Reconstructing a beam output according to the steering matrix and the optimized sound pressure cross-spectrum matrix: according to conventional cross-spectral beam forming theory, the beam output can be expressed in the form:
Sout＝wCwH(17)
wherein w is a steering matrix composed of steering vectors, and the column vectors are the steering vectors. The steering vector mainly plays a role in increasing the beam output from the sound source direction by focusing the beam to the sound source direction, and the beam output in the non-sound source direction is suppressed, so that the sound source direction can be identified and located by observing the focused main lobe direction. Based on the above idea, on the basis of the obtained turning matrix based on generalized inverse beamforming and the optimized sound pressure cross-spectrum matrix, a cross-spectrum beamforming based on generalized inverse beamforming can be reconstructed, and the beam output form is as follows:
Ssuperior food＝w(n)CSuperior food(w(n))H(18)
Sound pressure cross spectrum matrix C after optimization in the above formulaSuperior foodMay alternatively be in the form of:
Ssuperior food＝w(n)(C-diag(x))(w(n))H(19)
Upper formula w(n)That is, the steering matrix obtained by the generalized inverse beamforming in equation (4) can focus the beam output direction by the steering matrix, thereby realizing the sound source identification and localization. S output in the above formulaSuperior foodSimilar in form to the output form of conventional cross-spectral beamforming, but is superior because its steering matrix undergoes generalized inverse beamformingAnd the reconstruction precision of the algorithm is further improved. Meanwhile, due to the introduction of the optimized sound pressure cross-spectrum matrix, the anti-noise performance of the sound source reconstruction method is improved. In order to simplify the output of the above formula, the diagonal elements are taken from the above formula (19) and are squared, and finally the output of the cross-spectrum generalized inverse beam forming based on the optimized cross-spectrum matrix is obtained as follows:
Soutput of＝(diag(SSuperior food))1/2(20)。
Claims (1)
1. A cross-spectrum generalized inverse beam forming method based on cross-spectrum optimization is characterized in that,
building a sound source reconstruction system:
the system comprises a camera, a microphone array, a multi-channel signal collector and a post-processing computer, wherein the camera, the microphone array, the multi-channel signal collector and the post-processing computer are positioned near a sound source; the camera is used for collecting a space image of the whole sound field; the microphone array is used for measuring sound field signals; the multi-channel signal collector is used for measuring time domain analog signals, converting the time domain analog signals into digital signals and transmitting the digital signals to the computer for post-processing; the main functions of the computer are mainly storing and post-processing the collected sound field signals and efficiently displaying the distribution condition of the reconstructed sound field;
the beamforming process comprises the steps of:
1) computing a steering matrix based on generalized inverse beamforming:
the sound source surface is provided with N unit sound sources, and then the microphones r in the microphone arrayiThe measured sound pressure signal at a point is as follows:
in the above formula, i is 1,2, …, M, M is the number of sensors,
|ri-rjl is r in the sound fieldiAt a unit excitation of r0Complex response of (a) q (r)j) Is the assumed monopole sound source intensity;the equation (1) expresses the transfer relationship between the sound source and the microphone array, and can be converted into a vector-matrix form for the subsequent analysis:
p＝Aq (2)
wherein A is an NxM sound field transfer matrix, the element composition of which is determined by formula (1); p is an M-dimensional measurement sound pressure vector, and the element of the vector is the frequency domain response of the corresponding sensor at the single frequency f; q is an N-dimensional sound source intensity column vector, and the element components of the N-dimensional sound source intensity column vector represent the intensity distribution at the corresponding point in the sound source plane;
to solve the reconstructed sound field distribution, equation (2) can be transformed into a minimization problem:
β is GiHonokov regularization parameter, L is N multiplied by N regularization diagonal matrix, the minimization problem (3) is solved through an iterative algorithm, and finally the expression of the steering matrix in the step N can be obtained:
w(n)＝L(n)B(n)(4)
in the formula (4), n is the number of iterations, B(n)For an N × M matrix, the row vector can be expressed as follows:
B(i,:)(n+1)＝G(:,i)HALF/(G(:,i)HALF G(:,i)) (5)
the main function of the formula (5) is to normalize the transmission matrix, so as to ensure the accuracy of output; in the above formula, the superscript H represents conjugate transpose, G (: i) is the ith column vector of the transfer matrix G after normalization, wherein the expression of G is as follows:
G＝AL(n)(6)
the ALF in equation (5) is an invertible M × M dimensional matrix, and its expression has the following form:
ALF＝(GGH+βI)-1(7)
the regularization parameter β in the above formula (7) can be obtained by an L curve or a GVC method;
L(n)＝diag(|Sp|/||Sp||∞) (8)
in the above formula (8), diag is the general directionThe quantity is converted into a diagonal matrix form, | is the absolute value of the vector, | | | non-calculation∞Calculating infinite norm of vector; sp is the normalized generalized inverse beamforming output:
Sp＝w(n)p (9)
the following filtering process is also added in the solving process for solving the beam output Sp:
in the formula (10), (| Sp |)maxTaking the largest element in the vector | Sp | where Dk increases with increasing number of iterations (n), such that the filtering threshold gradually decreases with increasing number of iterations:
Dk＝0.2+0.5n (11)
solving to obtain a steering matrix w based on generalized inverse beam forming(n)；
2) Optimizing a sound pressure cross-spectrum matrix:
based on an M-dimensional sound pressure vector p obtained by sensor array measurement, an M multiplied by M-dimensional sound pressure cross spectrum matrix is obtained, and the expression is as follows:
C＝ppH(12)
the optimized sound pressure cross-spectrum matrix form is expressed as follows:
Csuperior food＝C-diag(x) (13)
Wherein, CSuperior foodTo eliminate the matrix of diagonal elements, diag (x) is a pair of diagonal matrices, x is the noise component;
the element x (i,1) in the vector x is limited to [0, C (i, i)]While maintaining some of the advantages of the diagonalization removal element, the method can reduce channel self-noise; in order to obtain the noise component x, a sound pressure cross-spectrum matrix can be optimized through an iterative optimization algorithm; assume a series of normalized unit vectors uk＝ek(k 1, 2.., m), uk denotes a normalized vector, the initial value of which is determined by ek;
ek represents a unit vector, and the following inequality holds:
uk H(C-diag(x))uk≥0 (14)
from the above equation, the following linear inequality can be derived:
wherein D (k, i) ═ u (i, k) <' > ceiling2Where u (i, k) is a vector ukI element of (1), d (k,1) ═ uk HCuk(ii) a A constraint Dx ≦ d for the noise component x is obtained by equation (15),
the core idea of the optimization of the sound pressure cross-spectrum matrix is to minimize a noise component-x, so that an objective function min F ═ x can be obtained; the above process of minimizing the noise component is summarized as follows, i.e. a linear programming problem:
solving the minimized noise component-x and outputting the optimized sound pressure cross-spectrum matrix CSuperior food；
3) Reconstructing beam output according to the steering matrix and the optimized sound pressure cross-spectrum matrix: the beam output is in the form:
Ssuperior food＝w(n)CSuperior food(w(n))H(18)
Sound pressure cross spectrum matrix C after optimization in the above formulaSuperior foodMay alternatively be in the form of:
Ssuperior food＝w(n)(C-diag(x))(w(n))H(19)
Taking diagonal elements from the above formula (19) and squaring to finally obtain the output of the cross-spectrum generalized inverse beam forming based on the optimized cross-spectrum matrix as follows:
Soutput of＝(diag(SSuperior food))1/2(20)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710315666.XA CN107153172B (en) | 2017-05-08 | 2017-05-08 | Cross-spectrum generalized inverse beam forming method based on cross-spectrum optimization |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710315666.XA CN107153172B (en) | 2017-05-08 | 2017-05-08 | Cross-spectrum generalized inverse beam forming method based on cross-spectrum optimization |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107153172A CN107153172A (en) | 2017-09-12 |
| CN107153172B true CN107153172B (en) | 2020-04-21 |
Family
ID=59793447
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710315666.XA Expired - Fee Related CN107153172B (en) | 2017-05-08 | 2017-05-08 | Cross-spectrum generalized inverse beam forming method based on cross-spectrum optimization |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107153172B (en) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109270492B (en) * | 2018-09-27 | 2022-11-25 | 重庆大学 | Regularization parameter selection method for large holographic distance |
| CN109541548B (en) * | 2018-11-22 | 2021-06-25 | 西安联丰迅声信息科技有限责任公司 | Air sonar positioning method based on matching field |
| CN109343003B (en) * | 2018-11-29 | 2022-11-11 | 重庆大学 | Method for identifying sound source formed by fast iterative shrinking wave beams |
| CN109631756B (en) * | 2018-12-06 | 2020-07-31 | 重庆大学 | Rotary sound source identification method based on mixed time-frequency domain |
| CN109613481A (en) * | 2019-01-10 | 2019-04-12 | 重庆大学 | A kind of Wave beam forming identification of sound source method adapting to wind tunnel test environment |
| CN111562316B (en) * | 2019-12-20 | 2021-03-23 | 襄阳达安汽车检测中心有限公司 | Material sound absorption coefficient measuring method based on double-sided array and generalized inverse algorithm |
| CN111580049B (en) * | 2020-05-20 | 2023-07-14 | 陕西金蝌蚪智能科技有限公司 | Dynamic target sound source tracking and monitoring method and terminal equipment |
| CN111736050B (en) * | 2020-08-28 | 2020-11-17 | 杭州兆华电子有限公司 | Partial discharge fault monitoring and evaluating device and method |
| CN113253197B (en) * | 2021-04-26 | 2023-02-07 | 西北工业大学 | Method for recognizing directivity of noise source of engine and part thereof |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103278799A (en) * | 2013-05-10 | 2013-09-04 | 中国计量学院 | Reverse beamforming method based on Toeplitz improvement of uniform linear array |
| CN103837871A (en) * | 2012-11-23 | 2014-06-04 | 中国科学院声学研究所 | Inverse beamforming method and system |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10078006B2 (en) * | 2013-07-22 | 2018-09-18 | Brüel & Kjær Sound & Vibration Measurement A/S | Wide-band acoustic holography |
-
2017
- 2017-05-08 CN CN201710315666.XA patent/CN107153172B/en not_active Expired - Fee Related
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103837871A (en) * | 2012-11-23 | 2014-06-04 | 中国科学院声学研究所 | Inverse beamforming method and system |
| CN103278799A (en) * | 2013-05-10 | 2013-09-04 | 中国计量学院 | Reverse beamforming method based on Toeplitz improvement of uniform linear array |
Non-Patent Citations (2)
| Title |
|---|
| Iterative regularization method in generalized inverse beamforming;Zhifei Zhang等;《Journal of Sound and Vibration》;20170228;第396卷;108-121 * |
| 基于反问题的正则化波束形成改进算法;张志飞等;《仪器仪表学报》;20150831;第36卷(第8期);1752-1758 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107153172A (en) | 2017-09-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107153172B (en) | Cross-spectrum generalized inverse beam forming method based on cross-spectrum optimization | |
| KR101238362B1 (en) | Method and apparatus for filtering the sound source signal based on sound source distance | |
| KR101415026B1 (en) | Method and apparatus for acquiring the multi-channel sound with a microphone array | |
| CN105556260B (en) | broadband acoustical holography | |
| US20150163577A1 (en) | Low noise differential microphone arrays | |
| CN107564539B (en) | Acoustic echo cancellation method and device facing microphone array | |
| CN109254261B (en) | Coherent signal null deepening method based on uniform circular array EPUMA | |
| CN111024209B (en) | Line spectrum detection method suitable for vector hydrophone | |
| CN109343003B (en) | Method for identifying sound source formed by fast iterative shrinking wave beams | |
| CN113063490B (en) | Sound field separation method based on sound pressure and particle vibration velocity double-sided measurement | |
| JP6763332B2 (en) | Sound collectors, programs and methods | |
| CN110444220B (en) | Multi-mode remote voice perception method and device | |
| CN111812581B (en) | Spherical array sound source direction-of-arrival estimation method based on atomic norms | |
| CN107170441B (en) | Circular ring array optimal frequency invariant response super-directivity beam forming method | |
| Yu et al. | Adaptive imaging of sound source based on total variation prior and a subspace iteration integrated variational bayesian method | |
| TWI429885B (en) | Method for visualizing sound source energy distribution in reverberant environment | |
| CN114001816A (en) | Acoustic imager audio acquisition system based on MPSOC | |
| CN107861114B (en) | Noise suppression method based on underwater acoustic array airspace reversal | |
| Morgenstern et al. | Joint spherical beam forming for directional analysis of reflections in rooms | |
| CN116309921A (en) | Delay summation acoustic imaging parallel acceleration method based on CUDA technology | |
| Samarasinghe et al. | On room impulse response between arbitrary points: An efficient parameterization | |
| CN110632579B (en) | Iterative beam forming method using subarray beam domain characteristics | |
| CN113381793A (en) | Coherent information source estimation-oriented non-grid direction-of-arrival estimation method | |
| Noh et al. | Identification of low-frequency noise sources in high-speed train via resolution improvement | |
| Gur | Modal beamforming for small circular arrays of particle velocity sensors |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200421 |