A multi-sound source localization method based on signal time-frequency correlation

Abstract

The present invention proposes a multi-sound source localization method that utilizes the correlation between signal time-frequency points to detect single source points through the time-domain correlation and local frequency stability of speech signals. First, the guiding time-frequency point is detected. Secondly, the direction of arrival information of the guiding point is used to detect single source points within the frame by measuring the directional deviation between the guiding point and the adjacent time-frequency points, and the frequency-domain correlated single sound source point of the guiding point is obtained. Thirdly, the time-domain correlated single sound source point of the guiding point is obtained by combining the inter-frame frequency sub-band correlation coefficient and the directional deviation measure between the guiding point and the adjacent time-frequency points. Then, the direction of arrival estimation values ​​of all detected points are calculated, and the angular interval of each sound source is estimated using kernel density estimation and peak search. Finally, the angular estimation value of each sound source is obtained through statistically weighted fine positioning.

Description

Technical Field

[0001] The present invention belongs to the field of sound source localization in the field of speech signal processing, and in particular to the problem of multi-sound source localization in high reverberation scenes. Background Art

[0002] Estimating the direction of arrival (DOA) of speech signals is crucial in acoustic signal processing. Its goal is to obtain the spatial location information of all sound sources, without any prior knowledge of the sound sources or the recording environment, using only the listening signals recorded by microphones placed in the listening environment. The application of microphone array networks has enabled multi-source DOA estimation under reverberant and noisy conditions. Related technologies, such as those used in urban monitoring, target detection, and human-computer interaction, are therefore being applied to more complex environments. However, multi-source DOA estimation in highly reverberant and noisy environments remains to be improved.

[0003] In practical applications, the acquisition of sound source location information can be negatively impacted by the aliasing of recorded signals caused by multiple sound sources simultaneously, limitations in microphone array structure, reverberation effects, and non-stationary noise in the scene. This can ultimately lead to errors in sound source location estimation and even incorrect sound source counting. To address these issues, scientists and researchers at home and abroad have conducted extensive research, proposing various sound source localization technologies, including those based on time difference of arrival, high-resolution spectrum estimation, beamforming, and sparse component analysis.

[0004] Among them, the method based on sparse component analysis mainly searches for time-frequency points in the time-frequency domain where the direct component of a single sound source is dominant, that is, single sound source points. The above positioning method transforms the multi-source localization problem into a single sound source localization problem by screening single sound source points and using the obtained single sound source components for subsequent positioning, thus realizing multi-source localization under underdetermined conditions. Under the same conditions, the positioning performance of this method is better than other positioning technologies of the same period. However, since the W-DO assumption is difficult to meet under the conditions of high reverberation and multiple sound sources, the detected time-frequency points are relatively few and contain a large proportion of outliers. Moreover, the proportion of outliers in the detection points always increases with the increase of reverberation time and the number of sound sources, which ultimately affects the sound source count and positioning accuracy. Summary of the Invention

[0005] In order to solve the problem of a small number of single sound source points detected by traditional single sound source point detection methods in reverberation and multi-sound source environments, the present invention designs a multi-source localization method based on the correlation of adjacent time-frequency points. The proposed method is applicable to any microphone array that can perform DOA estimation for each time-frequency point. This method utilizes the time domain correlation of speech signals, the local stability of frequency coefficients and the continuous distribution of single source points in the time-frequency domain to predict the adjacent points of single source points detected by traditional methods and determine whether they are single source points. This design is applicable to any microphone array that can perform DOA estimation for each time-frequency point. Taking into account the characteristics of sound field microphones that are light and flexible and can accurately capture sound field information, this design uses sound field microphones to achieve multi-sound source localization.

[0006] The overall design process is briefly described as follows:

[0007] First, the input four-channel sound field microphone signals are framed, windowed, and short-time Fourier transformed to obtain the time-frequency coefficients of each frame signal. Then, phase consistency detection is performed on the four microphone signals to obtain the time-frequency domain guidance points. The original recorded signals are then converted into B-format signals with strong spatial characteristics. The active sound intensity vectors of the guidance points and their adjacent points are calculated, and this is used to perform single-source point detection within and between frames. For single-source point detection within a frame, the specific steps are as follows: the directional deviation coefficient between the guidance point and its adjacent time-frequency points within the frame is calculated, and a unified deviation threshold is designed based on the actual situation. Time-frequency points that meet the threshold are detected as single-source points. Similarly, for single-source point detection between frames, the specific steps are as follows: the directional deviation coefficient between the guidance point and its adjacent time-frequency points and the correlation coefficient between the adjacent frames are calculated, and a unified deviation threshold and correlation measure validity threshold are designed based on the actual situation. Time-frequency points that meet the conditions are detected as single-source points. The DOA estimates of the single source point and the guide point are calculated, and the sound source angle interval is estimated by kernel density estimation and peak search for all the DOA estimates. Finally, the sound source angle estimate is obtained by statistical weighted fine positioning.

[0008] The technical solution of the present invention is to solve the problem of multi-sound source localization under reverberation conditions, which is mainly divided into the following steps:

[0009] Step 1: First, perform a short-time Fourier transform on the four acquired signals from the sound field microphones to obtain a time-frequency domain representation of each signal frame. Next, calculate the angular spacing between the real and imaginary vectors of the four microphone signals, set an angular spacing threshold, and perform time-frequency domain guidance point (abbreviated as "guide point"). Simultaneously, convert the A-format signal to B-format, and use this to calculate the activity intensity vector of the time-frequency point.

[0010] Step 2: Single source point detection within the frame. The directional offset of the activity intensity vector between the guide point and its adjacent time-frequency points is defined as the vector angular deviation. The vector angular deviation between two adjacent time-frequency points within the frame is mapped to angle information and normalized to obtain the normalized direction of arrival deviation coefficient between adjacent frequency points. The normalized direction of arrival deviation coefficients of all guide points and their adjacent time-frequency points within the frame are calculated, and the time-frequency points that meet the direction of arrival deviation threshold are determined to be single source points, which are called frequency-domain related single sound source points of the guide points.

[0011] Step 3: Calculate the normalized direction of arrival deviation coefficient between adjacent frames. Calculate the vector angle deviation between the guidance point and the same frequency point in its adjacent frames, map it to angle information, and normalize it to obtain the normalized direction of arrival deviation coefficient between adjacent frames at the same frequency.

[0012] Step 4: Single source point detection in adjacent frames. Calculate the frequency subband correlation coefficient of adjacent frames and combine it with the normalized direction of arrival deviation coefficient of adjacent frames for joint judgment. The time-frequency points that meet the judgment criteria are determined to be single source points, which are called time-domain correlated single sound source points of the guidance points.

[0013] Step 5: Calculate the DOA estimation values ​​of all detected single sound source points, and calculate the angle estimation interval of each sound source through kernel density estimation and peak detection.

[0014] In step 6, the DOA estimation value of each sound source is finally obtained through a fine positioning method based on statistical weighting.

[0015] 1. The implementation method of step 1 is to first perform short-time Fourier transform on the sound field microphone collected signal to obtain the time-frequency domain representation S of the collected signals of the four microphones in the sound field microphone array: left front upper, right front lower, rear left lower, and rear right upper FLU (n,k),S FRD (n,k),S BLD (n,k),S BRU (n,k), where frame index n = 1, 2, 3, ..., N, where N is the number of frames; frequency point index k = 1, 2, 3, ..., L, where L is the frame length, i.e. the number of samples or frequency domain points per frame. The A format signal vector is: X A (n,k)=[S FLU (n,k),S FRD (n,k),S BLD (n,k),S BRU (n,k)].

[0016] Calculate the angular distance between the real and imaginary vectors of the A-format signal at the time-frequency point (n, k) as follows:

[0017]

[0018] Where Re{·} and Im{·} are real and imaginary part operations respectively, T is the transpose operator, |·| is the absolute value operation, and ||·||2 represents the two-norm operation. The angular spacing of the time-frequency point (n, k) vector is The value range of is [0,1]. When , the time-frequency point (n, k) is a single sound source point; The larger the value, the higher the proportion of single sound source components in the time-frequency point. According to this criterion, the time-frequency point detected as a single sound source is called a time-frequency domain guidance point, referred to as a guidance point. The guidance point is used to determine whether its adjacent time-frequency point is a single sound source time-frequency point. The set of detected guidance points is as follows;

[0019]

[0020] Where ε1 is the phase consistency detection threshold, and the value range of this threshold is [0.85, 0.98]. g All conditions are met The set of time-frequency points (n,k).

[0021] Then, the B-format signal vector X is obtained by the following formula: B (n,k)=[S W (n,k),S X (n,k),S Y (n,k),S Z (n,k)]:

[0022]

[0023] The B format signal contains an omnidirectional channel signal S W (n, k) and three channel signals pointing to the positive direction of the Cartesian coordinate system {S X (n,k),S Y (n,k),S Z (n,k)}. Using X B (n,k) Calculate the activity intensity vector of the time-frequency point (n,k):

[0024]

[0025] Among them, I B (n,k) is the calculated activity intensity vector of the time-frequency point (n,k), * is the conjugate operator, Re{} is the real part operation, I X (n,k),I Y (n,k),I Z (n, k) is the component of activity intensity corresponding to the positive directions of the X, Y, and Z axes in the Cartesian coordinate system.

[0026] 2. In step 2, the activity intensity vector IB (n,k) contains the direction of arrival information of the sound source signal at the time-frequency point (n,k). Therefore, the vector angle deviation between two adjacent time-frequency points (n,k) and (n,k±1) in the nth frame is defined as:

[0027]

[0028] Among them, I B (n,k) is the calculated activity intensity vector of the time-frequency point (n,k), <·> is the dot product operation, and ||·||2 represents the two-norm operation. The vector angle deviation between the two time-frequency points is converted into angle information and normalized to obtain the normalized direction of arrival deviation coefficient C between adjacent frequency points. F (n,k,k±1), as follows:

[0029]

[0030] Where π is 180 degrees in radians.

[0031] C F (n,k,k±1) reflects the difference in the direction of arrival of the sound source corresponding to the time-frequency point (n,k) and the time-frequency point (n,k±1), and its value range is between 0 and 1. F The smaller (n,k,k±1) is, the greater the direction deviation is; on the contrary, C F The larger (n,k,k±1) is, the smaller the direction deviation is. F When (n,k,k±1)=1, it means that the estimated direction of arrival value of the sound source at the time-frequency point (n,k±1) is the same as the DOA estimated value at the time-frequency point (n,k).

[0032] If the time-frequency point (n, k) is the guidance point obtained in step 1 (i.e., (n, k)∈Q g ), and the time-frequency point (n,k+1) or (n,k-1) is not a guidance point (i.e., ), then we can use formula (6) to determine whether they are single-source time-frequency points. The new single-source time-frequency points detected by this step are called frequency-domain related single sound source points of the guidance points, and their set Q fu and Q fd for:

[0033]

[0034]

[0035] Among them, ε2 is the frequency domain correlation detection threshold, and the value range of this threshold is [0.85, 0.98]. fu All conditions are met The set of time-frequency points (n, k); Q fdAll conditions are met The set of time-frequency points (n,k).

[0036] 3. In step 3, based on the short-term continuous characteristics of the single-source time-frequency point in the time domain, the guidance point obtained in step 1 is used to detect whether the adjacent time-frequency point in the next frame is a single-source point. If the time-frequency point (n-1, k) is a guidance point (i.e., (n-1, k)∈Q g ), and the time-frequency point (n, k) is not detected as a guide point and a frequency-domain related single sound source point of the guide point (i.e., ), according to the time domain correlation of the speech signal, the time-frequency point (n, k) is most likely also a single sound source time-frequency point.

[0037] Similarly, the present invention uses the direction of arrival deviation between adjacent time-frequency points as the single sound source point judgment criterion. The vector angle deviation of the activity intensity vector between two time-frequency points (n-1, k) and (n, k) at the kth frequency of adjacent frames is defined as:

[0038]

[0039] Among them, C DT (n-1,n,k) is the vector angle deviation of the activity intensity vector between two time-frequency points, I B (n, k) is the calculated activity intensity vector of the time-frequency point (n, k), <·> is the dot product operation, and ||·||2 represents the two-norm operation. The vector angle deviation between the two time-frequency points is converted into angle information and normalized to obtain the normalized direction of arrival deviation coefficient C between the same frequency points in adjacent frames. T (n-1,n,k), as follows:

[0040]

[0041] Where π is 180 degrees in radians.

[0042] 4. In step 4, it is worth noting that the speech signal has obvious time-varying characteristics, and the speech signal may change suddenly between adjacent frames. Therefore, it is not possible to detect a single source point between frames only by the direction of arrival deviation coefficient. To this end, the present invention uses the frequency sub-band correlation measure between adjacent frames to perform sudden change discrimination. Define the frequency sub-band (Ω k =[kL F ,…,k,…,k+L F ]) is:

[0043] R F (n,k)=[|S W (n,kL F )|,…,|S W (n,k)|,…,|SW (n,k+L F )|]#(11)

[0044] Among them, S W (n, k) is the omnidirectional channel signal in B format, |·| is the absolute value operation, and the frequency subband contains 2 F +1 frequency point. The present invention defines the inter-frame frequency sub-band correlation coefficient between the time-frequency point (n-1, k) and the time-frequency point (n, k) as:

[0045]

[0046] Among them, |·| is the absolute value operation, ||·||2 represents the two-norm operation, and T is the transposition operator. T The range of (n,k) values ​​is [0,1], r T (n,k) reflects the difference between the n-1th frame signal and the nth frame signal in the sub-band Ω k If r T A larger (n, k) means that adjacent frames have a stronger correlation in this sub-band, that is, the probability of a sudden change at this frequency point between two frames is lower.

[0047] If the time-frequency point (n-1, k) is the guidance point obtained in step 1 (i.e., (n-1, k)∈Q g ), and the time-frequency point (n, k) is neither a guide point nor a frequency-domain related single sound source point of the guide point (i.e., Then, we can use formula (12) to determine whether (n, k) is a single source time-frequency point. The new single source time-frequency point detected by this step is called the time domain related single sound source point of the guidance point, and its set Q t for:

[0048]

[0049] Among them, ε3 is the time domain correlation detection threshold, and its value range is [0.85, 0.98]; ε4 is the frequency sub-band correlation coefficient threshold, and its value range is [0.6, 0.7]. t All satisfying condition C T (n-1,n,k)>ε3, The set of time-frequency points (n,k).

[0050] In summary, the time-frequency domain guidance points, the frequency-domain related single sound source points of the guidance points, and the time-domain related single sound source points of the guidance points are combined to obtain the final detected single sound source time-frequency points, whose set is defined as:

[0051] Q=Q g ∪Q fu ∪Q fd∪Q t #(14)

[0052] Q is the set of all time-frequency points (n, k) detected.

[0053] 5. In step 5, the DOA estimation is performed for the single sound source time-frequency point (i.e., (n, k) ∈ Q) detected in steps 2, 3, and 4. The purpose of DOA estimation is to obtain the direction of arrival information of the incident signal relative to the microphone array at the single sound source time-frequency point. The DOA estimate θ(n, k) of the time-frequency point (n, k) can be calculated using the activity intensity vector of the B-format signal, and its calculation formula is as follows:

[0054]

[0055] Among them, θ(n,k) is the DOA estimation value of the time-frequency point (n,k), I X (n,k),I Y (n, k) is the component of activity intensity corresponding to the positive direction of the X and Y axes in the Cartesian coordinate system.

[0056] The DOA estimation values ​​of all single sound source points are obtained using formula (15). The kernel density function is used to calculate the probability density function of the DOA estimation value in the angle range [0,359]. For the DOA estimation value set {θ(n,k)|(n,k)∈Q} obtained for a single sound source point, its probability density function is It is obtained by the following formula:

[0057]

[0058] Among them, K(·) is the kernel function, The smoothing parameter is selected as 20, and M is the number of elements in the set Q. h = 0, 1, ..., 359 represents the angle value, with a total of 360 angles.

[0059] The obtained probability density function Also known as the kernel density envelope estimate, it searches for the local maxima of the probability density function through peak search. The angle values ​​corresponding to these maxima are the estimated values ​​of the direction of arrival of each sound source, recorded as I E is the estimated number of sound sources.

[0060] In order to reduce the influence of the probability density peak between adjacent sound sources on the positioning accuracy, the present invention will As the rough estimate of the angle of the sound source, the angle estimation interval of the i-th sound source is defined as

[0061]

[0062] Among them, εθ is an offset constant whose value range is [2,5]. and The left and right side compensation coefficients are obtained by the following formula:

[0063]

[0064]

[0065] u and v are calculated using the following two formulas:

[0066]

[0067]

[0068] Among them I E is the estimated number of sound sources.

[0069] 6. Step 6: Using the DOA estimation value θ(n,k) of the detected single sound source time-frequency point (n,k), the histogram statistical vector of the DOA distribution is obtained, which is recorded as: L = [l0,l1,…,l 359 ]. Each element in L is obtained by the following formula:

[0070] l h =card(Q h )#(twenty two)

[0071] Here, h = 0, 1, ..., 359 represents the angle value, with a total of 360 angle values. The card(·) function is used to count the number of elements in a set. h is the set of time-frequency points of a single sound source with DOA estimated as h, defined as:

[0072] Q h ={(n,k)|θ(n,k)=h,(n,k)∈Q}#(23)

[0073] Among them, Q h is the set of all time-frequency points that satisfy the condition θ(n,k) = h, (n,k)∈Q. h = 0, 1, …, 359 represents the angle value, with a total of 360 angle values.

[0074] A refined search is performed by extracting DOA statistical information within the angle estimation interval of the sound source. The present invention designs a weighted sound source direction of arrival refined estimation method. Accordingly, the weight vector of the i-th sound source is defined as:

[0075]

[0076] in, The definition is as follows:

[0077]

[0078] Among them, θ i is the angle estimation interval of the i-th sound source, h = 0, 1, ..., 359, a total of 360 angle values. Then the estimated value of the direction of arrival of each sound source is obtained by the following formula:

[0079]

[0080] in is the final estimated value of the i-th sound source angle, with sound source index i=1,2,…,I E ,Θ i =[0,1,…,359] T is the angle matrix and 0 is the zero matrix of dimension 360×1.

[0081] Beneficial effects

[0082] This method uses the time-domain correlation and local frequency stability of speech signals to detect single sources. Compared to traditional single-source detection methods, it can obtain more single source points with higher accuracy. This improves the precision of multi-source localization, making the positioning performance more stable under reverberant conditions and applicable to complex scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] Figure 1 This is the overall framework diagram of this design method. Specific implementation methods

[0084] This embodiment is used to detect the direction of arrival of multiple sound sources under 400ms reverberation. The sound source is located in a 6.0m×4.0m×3.0m silent room environment. The sound field microphone is 1.5m above the ground, the sound source and the sound field microphone are located on the same horizontal plane, the sound source and the microphone are 1.7m apart, the angle interval between adjacent sound sources is 80°, and the number of sound sources is set to 4. The phase consistency detection threshold ε1 is selected as 0.984, the frequency domain correlation detection threshold ε2 and the time domain correlation detection threshold ε3 are selected as 0.944, the frequency sub-band correlation coefficient threshold ε4 is selected as 0.65, and the offset constant ε θ Select 4. The signal processing software is Matlab 2014a.

[0085] During implementation, the present invention embeds the algorithm into the software to realize the automatic operation of each process. The present invention is further explained below with specific implementation steps and accompanying drawings: The specific workflow is as follows:

[0086] Step 1: Perform time-frequency transformation on the microphone collected signal to obtain the time-frequency domain guidance point (guidance point), and calculate the activity intensity vector of the guidance point and its adjacent time-frequency points.

[0087] First, the short-time Fourier transform is performed on the sound field microphone acquisition signal to obtain the time-frequency domain representation S of the acquisition signal of the four microphones in the sound field microphone array: the upper left front, the lower right front, the lower left rear, and the upper right rear. FLU (n,k),S FRD (n,k),S BLD (n,k),S BRU (n, k), where frame index n = 1, 2, 3, ..., N, where N is the number of frames; frequency point index k = 1, 2, 3, ..., L, where L is the frame length (i.e., the number of samples or frequency domain points per frame). The A-format signal vector is: X A (n,k)=[S FLU (n,k),S FRD (n,k),S BLD (n,k),S BRU (n,k)]. Among them, S FLU (n,k),S FRD (n,k), S BLD (n,k),S BRU (n, k) is the time-frequency domain representation of the signals collected by the four microphones in the sound field microphone array: the upper left front, the lower right front, the lower left rear, and the upper right rear.

[0088] Calculate the angular distance between the real and imaginary vectors of the A-format signal at the time-frequency point (n, k) as follows:

[0089]

[0090] Among them, Re{·} and Im{·} are the real part and imaginary part operations respectively, T is the transposition operator, |·| is the absolute value operation, and ||·||2 represents the two-norm operation. is the angular spacing of the vectors, The value range is [0,1], which reflects the phase consistency of each channel signal in the microphone array. When , the phases of the signals in each channel are completely consistent, and the time-frequency point (n, k) is a single sound source point; The larger the value, the higher the phase similarity of the signals in each channel, and the higher the proportion of single sound source components in this time-frequency point. Based on this criterion, the time-frequency point detected as a single sound source is called a guidance time-frequency point, or simply a guidance point. The guidance point is used to determine whether the time-frequency point adjacent to it is a single sound source time-frequency point. The set of detected guidance points is as follows;

[0091]

[0092] Among them, ε1 is the phase consistency detection threshold, which is selected as 0.984 in this implementation. g All conditions are met The set of time-frequency points (n,k).

[0093] The B-format signal vector X is obtained by the following formula: B (n,k)=[S W (n,k),S X (n,k),S Y (n,k),S Z (n,k)]:

[0094]

[0095] The B format signal contains an omnidirectional channel signal S w (,k) and three channel signals pointing to the positive direction of the Cartesian coordinate system {S x (n,k),S y (n,k),S y (n,k)}. Using X B (n,k) calculates the activity intensity vector I at the time-frequency point (n,k) B (n,k):

[0096]

[0097] Among them, I B (n,k) is the calculated activity intensity vector of the time-frequency point (n,k), * is the conjugate operator, Re{} is the real part operation, I X (n,k),I Y (n,k),I Z (n, k) is the activity intensity component corresponding to the positive direction of the X, Y, and Z axes in the Cartesian coordinate system. Finally, the activity intensity vector of the guidance point and its adjacent time-frequency points is calculated and obtained by formula (4).

[0098] Step 2: Calculate the normalized direction of arrival deviation coefficients of all guide points and their adjacent time-frequency points in the frame, and use this to perform frequency domain correlation single sound source point detection of the guide points.

[0099] Activity intensity vector I B (n,k) contains the direction of arrival information of the sound source signal at the time-frequency point (n,k). The vector angle deviation between two adjacent time-frequency points (n,k) and (n,k±1) in the nth frame is calculated by the following formula:

[0100]

[0101] Among them, I B (n,k) is the calculated activity intensity vector of the time-frequency point (n,k), <·> is the dot product operation, and ||·||2 represents the two-norm operation. The vector angle deviation between the two time-frequency points is converted into angle information and normalized to obtain the normalized direction of arrival deviation coefficient C between adjacent frequency points. F(n,k,k±1), as follows:

[0102]

[0103] Where π is 180 degrees in radians.

[0104] C F (n,k,k±1) reflects the difference in the direction of arrival of the sound source corresponding to the time-frequency point (n,k) and the time-frequency point (n,k±1), and its value range is between 0 and 1. F The smaller (n,k,k±1), the greater the direction deviation. On the contrary, C F The larger (n,k,k±1) is, the smaller the direction deviation is. F When (n,k,k±1)=1, the estimated direction of arrival (DOA) of the sound source at the time-frequency point (n,k±1) is the same as the estimated DOA at the time-frequency point (n,k).

[0105] If the time-frequency point (n, k) is the guidance point obtained in step 1 (i.e., (n, k)∈Q g ), and the time-frequency point (n,k+1) or (n,k-1) is not a guidance point (i.e., ), then we can use formula (6) to determine whether they are single-source time-frequency points. The new single-source time-frequency points detected by this step are called frequency-domain related single sound source points of the guidance points, and their set Q fu and Q fd for:

[0106]

[0107]

[0108] Among them, ε2 is the frequency domain correlation detection threshold, and in this implementation, the threshold is selected as 0.944. fu All conditions are met The set of time-frequency points (n, k); Q fd All conditions are met The set of time-frequency points (n,k).

[0109] Step 3: Calculate the normalized direction of arrival deviation coefficient between frames.

[0110] According to the short-term continuous characteristics of the single-source time-frequency point in the time domain, the guidance point obtained in step 1 is used to detect whether the adjacent time-frequency point in the next frame is a single-source point. If the time-frequency point (n-1, k) is a guidance point (i.e., (n-1, k)∈Q g ), and the time-frequency point (n, k) is not detected as a guide point and a frequency-domain related single sound source point of the guide point (i.e., ), according to the time domain correlation of the speech signal, the time-frequency point (n, k) is most likely also a single sound source time-frequency point.

[0111] Similarly, the present invention uses the direction of arrival deviation between adjacent time-frequency points as the single sound source point determination criterion. The vector angle deviation of the activity intensity vector between two time-frequency points (n-1, k) and (n, k) at the kth frequency in adjacent frames is calculated using the following formula:

[0112]

[0113] Among them, C DT (n-1,n,k) is the vector angle deviation of the activity intensity vector between two time-frequency points, I B (n,k) is the calculated activity intensity vector of the time-frequency point (n,k), <·> is the dot product operation, and ||·||2 represents the two-norm operation. The vector angle deviation between the two time-frequency points is converted into angle information and normalized to obtain the normalized direction of arrival deviation coefficient C between the same frequency points in adjacent frames. T (n-1,n,k), as follows:

[0114]

[0115] Where π is expressed in radians of 180 degrees. The inter-frame normalized direction of arrival deviation coefficient C is calculated by the above formula T (n-1,n,k).

[0116] Step 4: Calculate the inter-frame correlation coefficient and combine it with the inter-frame normalized direction of arrival deviation coefficient to perform time-domain correlation single sound source point detection of the guidance point.

[0117] The frequency sub-band (Ω k =[kL F ,…,k,…,k+L F ])’s power coefficient vector:

[0118] R F (n,k)=[|S W (n,kL F )|,…,|S W (n,k)|,…,|S W (n,k+L F )|]#(11)

[0119] Among them, S W (n, k) is the omnidirectional channel signal in B format, |·| is the absolute value operation, and the frequency subband contains 2L F +1 frequency point. The present invention calculates the inter-frame correlation coefficient between the time-frequency point (n-1, k) and the time-frequency point (n, k) by the following formula:

[0120]

[0121] Among them, |·| is the absolute value operation, ||·||2 represents the two-norm operation, and T is the transposition operator. T The range of (n,k) values ​​is [0,1], r T (n,k) reflects the difference between the n-1th frame signal and the nth frame signal in the sub-band Ω k If r T A larger (n, k) means that adjacent frames have a stronger correlation in this sub-band, that is, the probability of a sudden change at this frequency point between two frames is lower.

[0122] If the time-frequency point (n-1, k) (i.e. (n, k)∈Q g ) is the guidance point obtained in step 1, and the time-frequency point (n, k) is neither a guidance point nor a frequency-domain related single sound source point of the guidance point, that is, Then, we can use formula (12) to determine whether (n, k) is a single source time-frequency point. The new single source time-frequency point detected by this step is called the time domain related single sound source point of the guidance point, and its set Q t for:

[0123]

[0124] Among them, ε3 is the time domain correlation detection threshold, which is 0.944 in this implementation; ε4 is the frequency sub-band correlation coefficient threshold, which is 0.65 in this implementation. t All satisfying condition C T (n-1,n,k)>ε3, The set of time-frequency points (n,k).

[0125] In summary, the time-frequency domain guidance points, the frequency-domain related single sound source points of the guidance points, and the time-domain related single sound source points of the guidance points are combined to obtain the final detected single sound source time-frequency points, whose set is defined as:

[0126] Q=Q g ∪Q fu ∪Q fd ∪Q t #(14)

[0127] Among them, Q is the set of all time-frequency points (n, k) detected.

[0128] Step 5: Calculate the DOA estimation values ​​of all detected single sound source points, and calculate the angle estimation interval of each sound source through kernel density estimation and peak detection.

[0129] Perform DOA estimation on the single sound source time-frequency points (i.e., (n, k)∈Q) detected in steps 2, 3, and 4. The purpose of DOA estimation is to obtain the direction of arrival information of the incident signal relative to the microphone array at the measured single sound source time-frequency point. The DOA estimate θ(n, k) at the time-frequency point (n, k) can be calculated using the activity intensity vector of the B-format signal, and its calculation formula is as follows:

[0130]

[0131] Among them, θ(n,k) is the DOA estimation value of the time-frequency point (n,k), I X (n,k),I Y (n, k) is the component of activity intensity corresponding to the positive direction of the X and Y axes in the Cartesian coordinate system.

[0132] The DOA estimation values ​​of all single sound source points are obtained using formula (15). The kernel density function is used to calculate the probability density function of the DOA estimation value in the angle range [0,359]. For the DOA estimation value set {θ(n,k)|(n,k)∈Q} obtained for a single sound source point, its probability density function is It is obtained by the following formula:

[0133]

[0134] Among them, K(·) is the kernel function, The smoothing parameter is selected as 20, and M is the number of elements in the set Q. h = 0, 1, ..., 359 represents the angle value, with a total of 360 angles.

[0135] The obtained probability density function Also known as the kernel density envelope estimate, it searches for the local maxima of the probability density function through peak search. The angle values ​​corresponding to these maxima are the estimated values ​​of the direction of arrival of each sound source, recorded as I E is the estimated number of sound sources.

[0136] In order to reduce the influence of the probability density peak between adjacent sound sources on the positioning accuracy, the present invention will As the rough estimate of the angle of the sound source, the angle estimation interval of the i-th sound source is calculated by the following formula:

[0137]

[0138] Among them, ε θ is an offset constant, which is set to 4 in this implementation. and The left and right side compensation coefficients are obtained by the following formula:

[0139]

[0140]

[0141] u and v are calculated using the following two formulas:

[0142]

[0143]

[0144] Among them I E is the estimated number of sound sources.

[0145] Step 6: Calculate the estimated direction of arrival of each sound source through a statistically weighted fine positioning method.

[0146] Using the DOA estimation value θ(n,k) of the detected single sound source time-frequency point (n,k), the histogram statistical vector of the DOA distribution is obtained, which is recorded as: L = [l0,l1,…,l 359 ]. Where l h Obtained by the following formula:

[0147] l h =card(Q h )#(twenty two)

[0148] Where h = 0, 1, ..., 359 represents the angle value, a total of 360 angles. The card(·) function is used to count the number of elements in a set. h is the set of time-frequency points of a single sound source with DOA estimated as h, which is obtained by the following formula:

[0149] Q h ={(n,k)|θ(n,k)=h,(n,k)∈Q}#(23)

[0150] Among them, Q h is the set of all time-frequency points that satisfy the condition θ(n,k) = h, (n,k)∈Q. h = 0, 1, …, 359 represents the angle value, with a total of 360 angle values.

[0151] A refined search is performed by extracting DOA statistical information within the angle estimation interval of the sound source. The present invention designs a weighted method for refining the direction of arrival of the sound source. Accordingly, the weight vector of the i-th sound source is obtained by the following formula:

[0152]

[0153] in Calculated by the following formula:

[0154]

[0155] Among them, θi is the angle estimation interval of the i-th sound source, h = 0, 1, ..., 359, a total of 360 angle values. Then the estimated value of the direction of arrival of each sound source is obtained by the following formula:

[0156]

[0157] in, is the final estimated value of the i-th sound source angle, with sound source index i=1,2,…,I E ,Θ i =[0,1,…,359] T is the angle matrix and 0 is the zero matrix of dimension 360×1.

[0158] The specific embodiments described herein are merely illustrative of the spirit of the present invention. Persons skilled in the art may make various modifications, additions, or substitutions to the described specific embodiments without departing from the spirit of the present invention or exceeding the scope of the appended claims.

Claims

1. A method for localizing multiple sound sources using signal time-frequency correlation, characterized in that The following steps are involved: Step 1: Perform time-frequency transformation on the signal collected by the microphone to obtain the guiding time-frequency point, and calculate the activity intensity vector of the guiding point and its adjacent time-frequency points; Step 2: Calculate the normalized direction of arrival deviation coefficients of all guide points and their adjacent time-frequency points in the frame, and use this to perform frequency domain correlation single sound source point detection of the guide points; Step 3, calculate the inter-frame normalized direction of arrival deviation coefficient; Step 4: Calculate the inter-frame correlation coefficient and combine it with the inter-frame normalized direction of arrival deviation coefficient to perform time domain correlation single sound source point detection of the guidance point; Step 5: Calculate the DOA estimation values ​​of all detected single sound source points, and calculate the angle estimation interval of each sound source through kernel density estimation and peak detection; Step 6: Calculate the estimated direction of arrival of each sound source through a statistically weighted fine positioning method; Step 1: Get the guidance point and calculate the activity intensity vector: Acquire the guidance point; the A format signal vector is: X A (n,k)=[S FLU (n,k),S FRD (n,k),S BLD (n,k),S BRU (n,k)]; where S FLU (n,k),S FRD (n,k),S BLD (n,k),S BRU (n, k) is the time-frequency domain representation of the signals collected by the four microphones in the sound field microphone array: the upper left front, lower right front, lower left rear, and upper right rear. The following is the vector angular spacing of the microphone signals collected by each channel at the time-frequency point (n, k): Among them, Re{·} and Im{·} are the real part and imaginary part operations respectively, T is the transpose operator, |·| is the absolute value operation, and ||·||2 represents the two-norm operation; is the angular spacing of the vectors, The value range of is [0,1], which reflects the phase consistency of each channel signal in the microphone array; When , the phases of the signals in each channel are completely consistent, and the time-frequency point (n, k) is a single sound source point; The larger the value, the higher the phase similarity of the signals in each channel, and the higher the proportion of single sound source components in the time-frequency point. According to this criterion, the time-frequency point detected as a single sound source is called a guidance time-frequency point, or simply a guidance point. The guidance point is used to determine whether the adjacent time-frequency point is a single sound source time-frequency point. The set of detected guidance points is as follows: Among them, ε1 is the phase consistency detection threshold, and the value range of this threshold is [0.85, 0.98]. In this formula, Q g All conditions are met The set of time-frequency points (n,k); The B-format signal vector X is obtained by the following formula: B (n,k)=[S W (n,k),S X (n,k),S Y (n,k),S Z (n,k)]: The B format signal contains an omnidirectional channel signal S w (n, k) and three channel signals pointing to the positive direction of the Cartesian coordinate system {S X (n,k),S Y (n,k),S Z (n,k)}; using X B (n,k) Calculate the activity intensity vector of the time-frequency point (n,k): Among them, I B (n,k) is the calculated activity intensity vector of the time-frequency point (n,k), * is the conjugate operator, Re{} is the real part operation, I X (n,k),I Y (n,k),I Z (n, k) are the components of activity intensity corresponding to the positive directions of the X, Y, and Z axes in the Cartesian coordinate system; Calculate the normalized direction of arrival deviation coefficients of all guide points and their adjacent time-frequency points in the frame, and use them to perform frequency domain correlation single sound source point detection of the guide points; The vector angle deviation between two adjacent time-frequency points (n, k) and (n, f±1) in the nth frame is calculated by the following formula: Among them, I B (n, k) is the calculated activity intensity vector of the time-frequency point (n, k), <·> is the dot product operation, and ||·||2 represents the two-norm operation; the vector angle deviation between the two time-frequency points is converted into angle information and normalized to obtain the normalized direction of arrival deviation coefficient C between adjacent frequency points F (n,k,k±1), as follows: Where π is 180 degrees in radians; C F (n,k,k±1) reflects the difference in the direction of arrival of the sound source corresponding to the time-frequency point (n,k) and the time-frequency point (n,k±1), and its value range is between 0 and 1; C F The smaller (n,k,k±1) is, the greater the direction deviation is; on the contrary, C F The larger (n, k, k±1) is, the smaller the direction deviation is; when C F When (n,k,k±1)=1, the estimated direction of arrival (DOA) of the sound source at the time-frequency point (n,k±1) is the same as the DOA estimated at the time-frequency point (n,k); If the time-frequency point (n, k) is the guidance point obtained in step 1, that is, (n, k)∈Q g , and the time-frequency point (n,k+1) or (n,k-1) is not a guide point, that is The normalized direction of arrival deviation coefficient between adjacent frequency points is used to determine whether they are single-source time-frequency points; the new single-source time-frequency points detected by this step are called frequency-domain related single sound source points of the guidance points, and their set Q fu and Q fd for: Among them, ε2 is the frequency domain correlation detection threshold, and the value range of this threshold is [0.85, 0.98]. In this formula, Q fu All satisfying condition C F (n,k,k+1)>ε2,(n,k)∈Q g , The set of time-frequency points (n, k); Q fd All satisfying condition C F (n,k,k-1)>ε2,(n,k)∈Q g , The set of time-frequency points (n,k).

2. The method for localizing multiple sound sources using signal time-frequency correlation as claimed in claim 1, wherein: Calculate the normalized direction of arrival deviation coefficient between frames; Using the guidance point obtained in step 1, detect whether the adjacent time-frequency point in the next frame is a single source point; if the time-frequency point (n-1, k) is a guidance point, that is, (n-1, k)∈Q g , while the time-frequency point (n, k) is not detected as the guidance point and the frequency domain related single sound source point of the guidance point, According to the time domain correlation of the speech signal, the time-frequency point (n, k) is most likely also a single sound source time-frequency point; The vector angle deviation of the activity intensity vector between two time-frequency points (n-1, k) and (n, k) at the kth frequency in adjacent frames is calculated by the following formula: Among them, C DT (n-1,n,k) is the vector angle deviation of the activity intensity vector between two time-frequency points, I B (n, k) is the calculated activity intensity vector of the time-frequency point (n, k), <·> is the dot product operation, and ||·||2 represents the two-norm operation; the vector angle deviation between the two time-frequency points is converted into angle information and normalized to obtain the normalized direction of arrival deviation coefficient C between the same frequency points in adjacent frames T (n-1,n,k), as follows: Where π is expressed in radians of 180 degrees; the inter-frame normalized direction of arrival deviation coefficient C is calculated by the above formula T (n-1,n,k); Calculate the estimated direction of arrival of each sound source through a fine positioning method based on statistical weighting; Using the DOA estimation value θ(n,k) of the detected single sound source time-frequency point (n,k), the histogram statistical vector of the DOA distribution is obtained, which is recorded as: L = [l0,l1,…,l 359 ]; where l h Obtained by the following formula: l h =card(Q h ) Where h = 0, 1, ..., 359 represents the angle value, a total of 360 angles; the function card(·) is used to count the number of elements in the set; Q h is the set of time-frequency points of a single sound source with DOA estimated as h, which can be obtained by the following formula: Q h {(n,k)|θ(n,k)=h,(n,k)∈Q} Among them, Q h is the set of all time-frequency points that satisfy the condition θ(n,k)=h,(n,k)∈Q; h=0,1,…,359 represents the angle value, a total of 360 angle values; A fine search is performed by extracting the DOA statistical information within the angle estimation interval of the sound source; the weight vector of the i-th sound source is calculated by the following formula: in, Obtained by the following formula: Among them, θ i is the angle estimation interval of the i-th sound source, h = 0, 1, ..., 359, a total of 360 angle values; then the estimated value of the direction of arrival of each sound source is obtained by the following formula: in, is the final estimated value of the direction of arrival of the i-th sound source, where the sound source index i = 1, 2, ..., I E ,Θ i =[0,1,…,359] T is the angle matrix and 0 is the zero matrix of dimension 360×1.

3. The method for localizing multiple sound sources using signal time-frequency correlation as claimed in claim 1, wherein: Calculate the inter-frame correlation coefficient and combine it with the inter-frame normalized direction of arrival deviation coefficient to perform time-domain correlation single sound source point detection of the guidance point; First, calculate the frequency sub-band (Ω k =[kL F ,…,k,…,k+L F ])’s power coefficient vector: R F (n,k)=[|S W (n,k-L F )|,…,|S W (n,k)|,…,|S W (n,k+L F )|] Among them, S W (n, k) is the omnidirectional channel signal in B format, |·| is the absolute value operation, and the frequency subband contains 2L F +1 frequency point; then the inter-frame correlation coefficient between the time-frequency point (n-1, k) and the time-frequency point (n, k) is obtained by the following formula: Among them, |·| is the absolute value operation, ||·||2 represents the two-norm operation, T is the transposition operator; r T The range of (n,k) values ​​is [0,1], r T (n,k) reflects the difference between the n-1th frame signal and the nth frame signal in the sub-band Ω k The correlation in T A larger (n,k) means that adjacent frames have a stronger correlation in this subband, that is, the probability of a sudden change at this frequency point between two frames is low; If the time-frequency point (n-1, k) is the guidance point obtained in step 1, that is, (n-1, k)∈Q g , and the time-frequency point (n, k) is neither a guiding point nor a frequency-domain related single sound source point of the guiding point, that is, The inter-frame normalized direction of arrival deviation coefficient and inter-frame correlation coefficient are used to determine whether (n, k) is a single-source time-frequency point; the new single-source time-frequency point detected by this step is called the time-domain correlated single sound source point of the guidance point, and its set Q t for: Q t ={(n,k)|C T (n-1,n,k)>ε3,r T (n,k)>ε4,(n-1,k)∈Q g , Among them, ε3 is the time domain correlation detection threshold, and its value range is [0.85, 0.98]; ε4 is the frequency sub-band correlation coefficient threshold, and its value range is [0.6, 0.7]; in this formula, Q t All satisfying condition C T (n-1,n,k)>ε3,r T (n,k)>ε4,(n-1,k)∈Q g , The set of time-frequency points (n,k); In summary, the guiding time-frequency points, the frequency-domain-correlated single sound source points of the guiding points, and the time-domain-correlated single sound source points of the guiding points are combined to obtain the final detected single sound source time-frequency points, the set of which is as follows: Q=Q g ∪Q fu ∪Q fd ∪Q t Q is the set of all time-frequency points (n, k) detected.

4. The method for localizing multiple sound sources using signal time-frequency correlation as claimed in claim 1, wherein: Calculate the DOA estimation values ​​of all detected single sound source points, and calculate the angle estimation interval of each sound source through kernel density estimation and peak detection: First, the DOA estimate of the guidance point and its adjacent points is calculated using the following formula: Among them, θ(n,k) is the DOA estimation value of the time-frequency point (n,k), I X (n,k),I Y (n, k) is the component of activity intensity corresponding to the positive directions of the X and Y axes in the Cartesian coordinate system; The above formula is used to obtain the DOA estimation values ​​of all single sound source points; the kernel density function is used to calculate the probability density function of the DOA estimation value in the angle range [0,359]; for the DOA estimation value set {θ(n,k)|(n,k)∈Q} of the single sound source point, its probability density function It is obtained by the following formula: Among them, K(·) is the kernel function, 20 is selected as the smoothing parameter, M is the number of elements in the set Q; h = 0, 1, ..., 359 represents the angle value, a total of 360 angles; The obtained probability density function Also known as the kernel density envelope estimate, it searches for the local maximum of the probability density function through peak search; the angle values ​​corresponding to these maximum values ​​are the estimated values ​​of the direction of arrival of each sound source, recorded as I E is the estimated number of sound sources; In order to reduce the impact of the probability density peak between adjacent sound sources on the positioning accuracy, As a rough estimate of the angle of the sound source, the angle estimation interval of the i-th sound source is calculated by the following formula: Among them, ε θ is an offset constant whose value range is [2,5]; and The left and right side compensation coefficients are obtained by the following formula: u and v are calculated using the following two formulas: Among them I E is the estimated number of sound sources.

Images