CN109870669B - Two-dimensional multi-snapshot meshless compressed beam forming sound source identification method - Google Patents

Abstract

The invention discloses a two-dimensional multi-snapshot gridless compression beamforming sound source identification method, comprising the following steps: step 1, obtaining the measured sound pressure matrix P ^★ ; step 2, iteratively seeking the sound pressure P: 1, using MATLAB The SDPT3 solver in the CVX toolbox solves:

2. Determine the regularization parameter κ ^l+1 of the l+1 iteration; 3. Determine the weight matrix W ^l+1 of the l+1 iteration; when two consecutive iterations

When the relative variation between is less than or equal to ^10-3 or the maximum number of iterations is completed, the iteration is terminated; step 3, estimating the DOA of the sound source; step 4, estimating the intensity of the sound source. The invention can accurately estimate the DOA with small distances between sound sources, and can quantitatively obtain the intensity of the sound sources, thereby improving resolution, noise removal ability and sound source recognition accuracy.

Description

A two-dimensional multi-snapshot gridless compression beamforming method for sound source identification

technical field

The invention belongs to the technical field of sound field identification.

Background technique

范数。上述假设不成立时，其重构结果不准确，该问题称作基不匹配，在实际应用中经常出现。Compressed beamforming based on planar microphone array measurement is an effective way to realize two-dimensional direction-of-arrival (DOA) estimation and intensity quantification of sound sources. Traditional compressed beamforming assumes that the sound source is distributed on a set of preset discrete grid points, each grid point represents an observation direction, and the sparse constraint is reflected by minimizing the distribution vector of the sound source

norm. When the above assumptions are not established, the reconstruction results will be inaccurate. This problem is called base mismatch, which often occurs in practical applications.

In order to fundamentally solve this problem, two-dimensional single-snapshot and multi-snapshot gridless compression beamforming strategies have been developed successively. Compared with single-snapshot, multi-snapshot method is more robust, but the existing atomic norm-based Atomic NormMinimization (ANM) two-dimensional multi-snapshot gridless compression beamforming sound source identification has the defect that the distance between sound sources is small.

Contents of the invention

Aiming at the problems existing in the prior art, the technical problem to be solved by the present invention is to provide a two-dimensional multi-snapshot gridless compression beamforming sound source identification method, which uses iterative re-weighted atomic norm minimization, which can improve the resolution , Noise removal ability and sound source identification accuracy.

The technical problem to be solved in the present invention is realized by such technical scheme, and it comprises steps:

Step 1. Obtain the measured sound pressure matrix P ^★ ;

为：Measuring Sound Pressure Matrix

for:

P ^★ ＝P+N

为噪声干扰，

为复数集，A为矩形传声器阵列的行数，B为矩形传声器阵列的列数，L为快拍次数；The measured sound pressure matrix P is ^obtained through the measurement of the microphone array,

for noise interference,

is a set of complex numbers, A is the number of rows of the rectangular microphone array, B is the number of columns of the rectangular microphone array, and L is the number of snapshots;

Step 2, reconstructing the sound pressure P generated by the sound source at the array microphone;

Step 1), establish the mathematical model of reconstruction P

是该式计算结果，||·||_F表示Frobenius范数，P为待求声压矩阵，定义T_b(·)为二重Toeplitz算子，对任意给定向量u映射为A×A维的块Toeplitz型Hermitian矩阵：In the formula,

is the calculation result of this formula, ||·|| _F represents the Frobenius norm, P is the sound pressure matrix to be sought, define T _b (·) as a double Toeplitz operator, and map any given vector u into A×A dimension The block Toeplitz-type Hermitian matrix:

In T _b (u), each block T _a (0≤a≤A-1) is a B×B-dimensional Toeplitz matrix:

都是辅助量，κ＞0为规则化参数，

为单位矩阵，“tr(·)”表示求矩阵的迹，上标“H”表示共轭转置，“≥0”表示半正定；

are auxiliary quantities, κ>0 is a regularization parameter,

is the identity matrix, "tr( )" means to seek the trace of the matrix, the superscript "H" means conjugate transpose, and "≥0" means positive semi-definite;

Step 2), solve for P

κ⁰＝1，则W⁰＝I，基于第l次迭代的结果

κ^l，第l+1次迭代的步骤为：initialization

κ ⁰ =1, then W ⁰ =I, based on the result of the lth iteration

κ ^l , the steps of the l+1th iteration are:

1. Use the SDPT3 solver in the CVX toolbox of MATLAB to solve:

2. Determine the regularization parameter κ ^l+1 :

为

的最大特征值；

for

The largest eigenvalue of ;

当

小于等于10^-3或者最大迭代次数被完成时，迭代终止；3. The weight matrix W ^l+1 determined in the l+1 iteration:

when

When less than or equal to 10 ^-3 or the maximum number of iterations is completed, the iteration terminates;

Step 3, estimating the sound source DOA;

Step 4. Estimate the intensity of the sound source.

Technical effect of the present invention is:

The present invention can accurately estimate the DOA with small distance between sound sources, and can quantitatively obtain the intensity of the sound source, overcomes the defects of existing ANM's two-dimensional multi-snapshot gridless compression beamforming, improves the resolution and denoises ability and accuracy of sound source identification.

Description of drawings

The accompanying drawings of the present invention are as follows:

Figure 1 is the measurement layout of the microphone array;

Figure 2 is a comparison chart of simulation reconstruction results when the sound source frequency is 2000Hz, 3000Hz, 4000Hz and 4900Hz;

Fig. 2 (a), (c), (e), (g) are ANM; Fig. 2 (b), (d), (f), (h) are the present invention;

Fig. 3 is the comparison chart of the reconstruction accuracy of the sound pressure P, the estimation accuracy of the sound source DOA and the quantization accuracy of the sound source intensity between the present invention and the existing ANM;

随Δ_min的变化曲线；Figure 3(a)

Variation curve with Δ _min ;

随Δ_min的变化曲线；Figure 3(b)

Variation curve with Δ _min ;

随Δ_min的变化曲线；Figure 3(c)

Variation curve with Δ _min ;

Figure 4 is a test layout;

Figure 5 is a comparison chart of the experimental sound field reconstruction results when the sound source frequency is 2000Hz, 3000Hz, 4000Hz and 4900Hz;

Fig. 5 (a), (c), (e), (g) are ANM; Fig. 5 (b), (d), (f), (h) are the present invention.

Detailed ways

Below in conjunction with accompanying drawing and embodiment the present invention will be further described:

The present invention comprises the following steps:

Step 1. Obtain the measured sound pressure matrix P ^★

为各快拍下i号声源的强度(声源在(0,0)号传声器处产生的声压)组成的行向量，l＝1,2,…,L为快拍索引，s_i,l为第l快拍下i号声源的强度，

为复数集。Two-dimensional gridless compressed beamforming sound source identification is to use a rectangular microphone array to measure the sound signal. The measurement layout of the microphone array is shown in Figure 1, the symbol "●" indicates the microphone, a=0,1,...,A-1, b=0,1,...,B-1 are the microphone indexes of x and y dimensions respectively , Δx, Δy are the microphone spacing in x and y dimensions, respectively, θ _i , φ _i are the elevation angle and azimuth angle of sound source i DOA (0°≤θ≤90°, 0°≤φ≤360°), respectively. remember

It is a row vector composed of the intensity of the sound source i in each snapshot (the sound pressure produced by the sound source at the (0,0) microphone), l=1,2,...,L is the index of the snapshot, _{si, l} is the intensity of the sound source i in the first snapshot,

is a complex set.

可表示为：Assuming that the sound source radiates a plane sound wave with a wavelength of λ, t _1i ≡sinθ _i cosφ _i Δx/λ, t _2i ≡sinθ _i sinφ _i Δy/λ, then each snapshot of the sound source is generated at the microphones (a,b) row vector of sound pressures

Can be expressed as:

为虚数单位。In formula (1), k is the total number of sound sources,

is an imaginary unit.

Build the matrix:

标量

行向量

其中，上标“T”表示转置，符号

表示Kronecker乘积，“||·||₂”表示

范数，

为正实数集，||ψ_i||₂＝1，与式(1)对应，有：Column vector

scalar

row vector

Among them, the superscript "T" means transpose, and the symbol

Indicates the Kronecker product, "||·|| ₂ "indicates

norm,

is a set of positive real numbers, ||ψ _i || ₂ = 1, corresponding to formula (1), we have:

时，测量声压矩阵

可表示为：There is noise interference

When the sound pressure matrix is measured

Can be expressed as:

P ^★ ＝P+N (3)

In the simulation experiment, the added noise is independent and identically distributed Gaussian white noise, and the Signal-to-noise Ratio (SNR) is defined as SNR=20log ₁₀ (||P|| _F /||N|| _F ), given by This can determine ||N|| _F =||P|| _F 10 ^−SNR/20 , where “||·|| _F ” denotes the Frobenius norm.

由传声器阵列测量得到。When testing the test,

Measured by the microphone array.

Step 2. Reconstruct the sound pressure P generated by the sound source at the array microphone

Step 1), establish the mathematical model of reconstruction P

The first step in the post-processing of 2D multi-snapshot ^gridless compression beamforming is to filter out the noise in the measured sound pressure P and reconstruct the sound pressure P produced by the sound source at the array microphone. This step is applied to the sound source Sparse constraints are implemented. Build a measure of source sparsity:

是关于P的函数，

都是辅助量，κ＞0为规则化参数，

为单位矩阵，“tr(·)”表示求矩阵的迹，上标“H”表示共轭转置，“≥0”表示半正定。In formula (4),

is a function of P,

are auxiliary quantities, κ>0 is a regularization parameter,

is the identity matrix, "tr(·)" means to seek the trace of the matrix, the superscript "H" means conjugate transpose, and "≥0" means positive semi-definite.

Define T _b ( ) as a double Toeplitz operator, and any given vector u is mapped to an A×A-dimensional block Toeplitz-type Hermitian matrix:

In formula (5), each block T _a (0≤a≤A-1) is a B×B-dimensional Toeplitz matrix:

The reconstruction problem of P can be written as:

In formula (7), ε is the noise control parameter, usually taken as ||N|| _F .

Simultaneous formula (4), (7) can get:

表示计算结果，P为待求声压矩阵。In formula (8),

Indicates the calculation result, and P is the sound pressure matrix to be sought.

度量声源稀疏性，P的重构问题可写为：The existing ANM two-dimensional multi-snapshot gridless compressed beamforming sound source identification method adopts the atomic norm of P

To measure the sparsity of the sound source, the reconstruction problem of P can be written as:

为凸函数，该式为凸优化问题，可转化为如下半正定规划后利用CVX工具箱中的SDPT3求解器求解：

is a convex function, and this formula is a convex optimization problem, which can be transformed into the following semi-positive definite programming and solved by the SDPT3 solver in the CVX toolbox:

不能高概率准确包含声源DOA信息或对声源数目的预测不稳健可靠而失效。When the separation between sound sources is small, the existing ANM two-dimensional multi-snapshot gridless compression beamforming will be obtained by the above formula

Failure to accurately contain sound source DOA information with a high probability, or the prediction of the number of sound sources is not robust and reliable.

Step 2), use this method invention to solve P

κ^l为第l次迭代确定的最优解和规则化参数，

第l次迭代确定的权矩阵，第l+1次迭代的替代函数式为：remember

κ ^l is the optimal solution and regularization parameters determined in the lth iteration,

The weight matrix determined in the lth iteration, the substitution function of the l+1th iteration is:

是

处的切平面，c^l为与变量无关的常数。忽略c^l，第l+1次迭代的最小化问题可写为：In formula (9),

Yes

The tangent plane at , c ^l is a constant that has nothing to do with variables. Neglecting c ^l , the minimization problem of the l+1th iteration can be written as:

Equation (10) is a convex optimization problem, which can be solved by the SDPT3 solver in the CVX toolbox.

κ⁰＝1，则W⁰＝I，为增强稀疏，令κ在迭代过程中逐渐减小，按下式变化：initialization

κ ⁰ ＝1, then W ⁰ ＝I, in order to enhance the sparseness, let κ gradually decrease in the iterative process, and change according to the following formula:

为

的最大特征值，l为迭代次数。In formula (11),

for

The largest eigenvalue of , l is the number of iterations.

κ⁰＝1，则W⁰＝I，基于第l次迭代的结果

κ^l，第l+1次迭代的步骤为：Therefore, initializing

κ ⁰ =1, then W ⁰ =I, based on the result of the lth iteration

κ ^l , the steps of the l+1th iteration are:

1. Use the SDPT3 solver in the CVX toolbox of MATLAB to solve:

2. Determine the regularization parameter κ ^l+1 :

为

的最大特征值；

for

The largest eigenvalue of ;

当连续两次迭代的

的相对变化量，即

小于等于10^-3或者最大迭代次数被完成时，迭代终止，最大迭代次数选为20。3. The weight matrix W ^l+1 determined in the l+1 iteration:

When two consecutive iterations of

The relative change of

When it is less than or equal to 10 ^-3 or the maximum number of iterations is completed, the iteration terminates, and the maximum number of iterations is selected as 20.

Step 3, estimating the sound source DOA;

The operation of this step is the same as the prior art ANM two-dimensional multi-snapshot gridless compression beamforming sound source identification method. References Z.Yang, L.Xie, P.Stoica.Vandermonde decomposition of multilevelToeplitz matrices with application to multidimensional super-resolution[J].IEEE Transactions on Information Theory,2016,62(6):3685-3701. (“Multiple Vandermonde decomposition of Toeplitz matrix and its application in multidimensional super-resolution[J].”, Yang Zai, Xie Lihua, P.Stoica, "IEEE Transactions on Information Theory", 2016,62(6):3685-3701).

及声源

Get the number of sound sources

and sound source

Step 4, estimating the intensity of the sound source;

为根据估计的声源DOA计算的感知矩阵，

为量化的各快拍下各声源的强度构成的矩阵，根据式(2)得到声源强度的计算式：remember

is the perceptual matrix calculated from the estimated sound source DOA,

To quantify the matrix formed by the intensity of each sound source under each snapshot, the formula for calculating the intensity of the sound source is obtained according to formula (2):

为使用本发明迭代求解所得声压。In the formula, A is the perceptual matrix, and the superscript "+" indicates the pseudo-inverse,

The resulting sound pressure is iteratively solved for using the present invention.

simulation test

In order to verify and establish the accuracy of the present invention and compare its performance improvement, a sound source identification simulation is carried out.

)依次为100dB、97dB、97dB、94dB、94dB、90dB(参考2×10^-5Pa))。The simulation settings are as follows: assume a point sound source with a specific intensity radiating sound waves of a specific frequency at a specific position (assuming six sound sources, the elevation angles are 60°, 60°, 50°, 40°, 30°, 70°, and the azimuth angles are sequentially 180°, 190°, 90°, 90°, 200°, 300°, intensity (RMS:

) are 100dB, 97dB, 97dB, 94dB, 94dB, 90dB in turn (refer to 2×10 ^-5 Pa)).

The specific flow of the simulation is:

1. Calculate the sound pressure signal received by each microphone according to formula (2) and formula (3) (measure the sound signal with a microphone array of A=B=8, Δx=Δy=0.035m, the signal-to-noise ratio SNR is 20dB, and the total number of snapshots set to 10);

2. According to the method invention of step 2, iteratively solve the sound pressure

3. Estimate the sound source DOA according to step 3;

4. Estimate the sound source intensity according to step 4.

依次为0.025、0.038、0.050、0.062。This method invention compares with the result of ANM method, as shown in Figure 2, the frequency of sound source radiation sound wave: Fig. 2 (a) and Fig. 2 (b) are 2000Hz, and Fig. 2 (c) and Fig. (d) are 3000Hz, Figure 2(e) and Figure 2(f) are 4000Hz; Figure 2(g) and Figure 2(h) are 4900Hz; in each figure, "*" indicates the reconstructed sound source distribution, and "○" indicates the real The distribution of sound sources, and for the convenience of comparison, each figure is scaled in dB with reference to the maximum output value, and the displayed dynamic range is 0~-20dB. At the same time, the standard sound pressure of 2.0×10 ^-5 Pa is the maximum output value for reference. Minimum separation between sound sources at 2000Hz, 3000Hz, 4000Hz, 4900Hz

They are 0.025, 0.038, 0.050, 0.062 in turn.

It can be seen from Figure 2 that the two-dimensional multi-snapshot gridless compression beamforming using ANM (Figure 2(a), Figure 2(c), Figure 2(e), Figure 2(g)) only accurately reconstructs The sound source distribution at 4900Hz is obtained, and the two-dimensional multi-snapshot gridless compression beamforming (Fig. 2(b), Fig. 2(d), Fig. 2(f), Fig. 2(h)) of the present invention is The reconstruction results at the four frequencies are all accurate, so it is proved that the present invention overcomes the defects of ANM and has higher resolution.

平均Frobenius范数误差

和标准

范数误差

来衡量P的重构精度、声源DOA的估计精度和声源强度的量化精度，其中，

为估计的和真实的声源仰角构成的列向量，

为估计的和真实的声源方位角构成的列向量，

为量化的和真实的声源强度的均方根构成的列向量。误差值越小，精度越高，计算结果越准确Define the standard Frobenius norm error

Average Frobenius norm error

and standard

Norm error

To measure the reconstruction accuracy of P, the estimation accuracy of sound source DOA and the quantization accuracy of sound source intensity, where,

is a column vector of estimated and true source elevation angles,

is a column vector of estimated and true source azimuths,

Column vector of quantized and true source intensities rms. The smaller the error value, the higher the precision and the more accurate the calculation result

It can be seen from Fig. 3 that compared with the two-dimensional multi-snapshot gridless compression beamforming using ANM, the present invention has stronger denoising ability and better resolution, and can more accurately reconstruct the sound source in Sound pressure generated at array microphones, more accurately estimating the DOA and quantifying the intensity of sound sources that are less separated from each other.

Test verification

In order to check the correctness of the simulation conclusion and the effectiveness and superiority of using the two-dimensional multi-snapshot gridless compression beamforming of the present invention in practical applications, a verification test is carried out.

公司的64个4958型传声器均匀分布的矩形阵列(A＝B＝8，Δx＝Δy＝0.035m)，声源为稳态白噪声信号激励的3个扬声器，从左至右，笛卡尔坐标依次为(2.24,0,5)m、(-1.24,0,5)m、(-2.24,0,5)m，仰角依次为24.13°、13.93°、24.13°，方位角依次为0°、180°、180°，由于试验在半消声室内进行，地面存在反射，扬声器关于地面还存在3个镜像声源，笛卡尔坐标依次为(2.24,-2.2,5)m、(-1.24,-2.2,5)m、(-2.24,-2.2,5)m，仰角依次为32.13°、26.80°、32.13°，方位角依次为315.52°、240.59°、224.48°。各传声器测得的声压经PULSE 3560D型数据采集系统同时采集并传输到PULSELABSHOP中进行快速傅里叶变换，得声压频谱，采样频率为16384Hz，每个快拍长1s、包含2¹⁴个采样点，10个快拍被采用。进一步，使用ANM(基于矩阵束与配对)和本发明后处理声压频谱来重构声源分布，其余设置均与仿真一致。Figure 4 is the test layout diagram, the measurement array is

The company's 64 4958-type microphones are uniformly distributed rectangular arrays (A=B=8, Δx=Δy=0.035m), and the sound source is 3 loudspeakers excited by steady-state white noise signals. From left to right, the Cartesian coordinates are sequential (2.24,0,5)m, (-1.24,0,5)m, (-2.24,0,5)m, the elevation angles are 24.13°, 13.93°, 24.13°, the azimuth angles are 0°, 180 °, 180°, because the test is carried out in a semi-anechoic chamber, there are reflections on the ground, and there are still three image sound sources of the speaker on the ground, and the Cartesian coordinates are (2.24,-2.2,5)m, (-1.24,-2.2 ,5)m, (-2.24,-2.2,5)m, the elevation angles are 32.13°, 26.80°, 32.13°, and the azimuth angles are 315.52°, 240.59°, 224.48°. The sound pressure measured by each microphone is simultaneously collected by the PULSE 3560D data acquisition system and transmitted to PULSELABSHOP for fast Fourier transformation to obtain the sound pressure spectrum. The sampling frequency is 16384Hz. Each snapshot is 1s long and contains ²¹⁴ samples. point, 10 snapshots are taken. Further, ANM (based on matrix beam and pairing) and the post-processing sound pressure spectrum of the present invention are used to reconstruct the sound source distribution, and the rest of the settings are consistent with the simulation.

Figure 5 shows the results of sound source frequencies of 2000Hz, 3000Hz, 4000Hz and 4900Hz. In each figure, "*" indicates the reconstructed sound source distribution, with reference to the maximum output value for dB scaling, and "○" indicates the real sound source DOA , without intensity information. Under the four frequencies, the minimum separation between sound sources is 0.032, 0.048, 0.065, 0.079 in turn.

It can be seen from Fig. 5 that using ANM, the reconstructed sound source distribution at 2000Hz (Fig. 5(a)) and 3000Hz (Fig. 5(c)) seriously deviates from the real sound source distribution, only at 4000Hz (Fig. 5(e) ) and 4900Hz (Fig. 5(g)) accurately estimated the DOA of the six main sources. Using the present invention, the sound source DOA estimation results at the four frequencies in Figs. 5(b), 5(d), 5(f) and 5(h) are all accurate.

The test conclusion is consistent with the simulation conclusion, which proves that the simulation conclusion is correct, and the two-dimensional multi-snapshot gridless compression beamforming of the present invention is more effective in practical application.

Claims (2)

1. a two-dimensional multi-snapshot gridless compression beamforming sound source identification method, is characterized in that, comprises the following steps: Step 1. Obtain the measured sound pressure matrix P ^★ ; The measured sound pressure matrix P ^★ ∈ C ^{AB × L} is: P ^★ ＝P+N Measured sound pressure matrix P ^★ Obtained by microphone array measurement, N∈C ^AB×L is noise interference, C is complex number set, A is the number of rows of rectangular microphone array, B is the number of columns of rectangular microphone array, L is the number of snapshots ; Step 2, reconstructing the sound pressure P generated by the sound source at the array microphone; Step 1), establish the mathematical model of reconstruction P

In the formula,

is the calculation result of this formula, ||·|| _F represents the Frobenius norm, P is the sound pressure matrix to be sought, define T _b (·) as a double Toeplitz operator, and map any given vector u into A×A dimension The block Toeplitz-type Hermitian matrix:

In T _b (u), each block T _a , 0≤a≤A-1, is a B×B-dimensional Toeplitz matrix:

Both are auxiliary quantities, κ>0 is the regularization parameter, I∈ΡAB ^×AB is the unit matrix, "tr( )" means the trace of the matrix, and the superscript " ^H " means the conjugate transpose,

Indicates positive semi-definite; Step 2), solve for P initialization

κ ⁰ =1, then W ⁰ =I, based on the result of the lth iteration

κ ^l , W ^l , the steps of the l+1th iteration are: 1. Use the SDPT3 solver in the CVX toolbox of MATLAB to solve:

2. Determine the regularization parameter κ ^l+1 :

for

The largest eigenvalue of ; 3. The weight matrix W ^l+1 determined in the l+1 iteration:

,when

When less than or equal to 10 ^-3 or the maximum number of iterations is completed, the iteration terminates; Step 3, estimating the sound source DOA; Step 4. Estimate the intensity of the sound source. 2. A method for identifying sound sources with two-dimensional multi-snapshot gridless compression beamforming according to claim 1, characterized in that the maximum number of iterations is 20.

Images