CN102147458B - Method and device for estimating direction of arrival (DOA) of broadband sound source - Google Patents

Abstract

The invention relates to a method and device for estimating the direction of arrival (DOA) of a broadband sound source. The method comprises the following steps: selecting one section of sound data X(n)=[x1(n),...,xP(n)]T from sound data collected from a microphone array, wherein n is a subset of Z*; evenly dividing the selected sound data X(n) into L frames; carrying out N-point quick Fourier transform on each frame of sound data Xl (n); presenting the frequency domain of the obtained sound data as follows: Xl(k)=[xl1(k),...,xlP(k)]T, wherein k=1,..., N/2; then, calculating the modulus square value eta (k) of a coherence function between two microphones in the microphone array; randomly selecting the modulus square values of q coherence functions from the eta (k); and estimating the DOA only by using a frequency band corresponding to the modulus squares of the q coherence functions. According to the method and device, the computation amount of an algorithm is greatly lowered on the premise of ensuring that the precision is hardly lowered.

Description

A method and device for estimating the direction of arrival of a broadband sound source

technical field

The invention relates to the technical field of target direction finding, in particular to a method and device for estimating the direction of arrival of a broadband sound source.

Background technique

Target direction finding has important applications in sonar, radar, navigation and wireless sensor networks and other fields. The sensor array is the main way to find the direction of the target.

In the target direction finding technology, an important technical problem is how to locate the broadband sound source. There are already many algorithms for this. The wideband MUSIC algorithm first transforms the wideband signal into the frequency domain through Fourier transform, then processes each frequency band with the narrowband MUSIC method, and finally combines the results of each frequency band to obtain the direction estimation of the sound source. That is, the Incoherent Sub-space Method (ISM, Incoherent Sub-space Method). This type of method has a small amount of calculation, but poor accuracy, and can only deal with incoherent signals; the fitting algorithm represented by the maximum likelihood (ML, Maximum Likelihood) estimation algorithm, with its good gradualness and flexible applicability, Its low computational complexity has attracted widespread attention. For broadband signals, Professor Yao Kung and others proposed the AML (Approximated Maximum Likelihood) direction of arrival estimation algorithm for broadband sound sources. This type of algorithm has high precision and can deal with coherent signals, but it has a large amount of computation. At the same time, the method of focusing matrix can also be used for coherent signal sources. First of all, it is necessary to use a low-complexity algorithm to estimate the direction of the sound source, construct a focusing matrix through this estimated angle, focus signals of different frequency bands on a center frequency through the focusing matrix, and then treat it as a narrow-band signal , the disadvantage of this method is that the estimated angle of the sound source needs to be known, and the selection of the estimated angle has a great influence on the performance of the algorithm.

Contents of the invention

The purpose of the present invention is to greatly reduce the calculation amount of the algorithm under the condition that the accuracy is basically not reduced.

In order to achieve the above purpose, a DOA estimation method for broadband sound sources is proposed. The specific steps of the method include:

Step 1): From the sound data collected by the microphone array, select a piece of sound data X(n)=[x ₁ (n), ..., x _P (n)] ^T , where P represents the number of microphones in the microphone array number, P≥2; the section of sound data X (n) is evenly divided into L frames, the length of the data X _l (n) of each frame is N, and the data length of X (n) is N×L, l=1 ,...,L; n∈Z ^* ;

Perform N-point fast Fourier transform on each frame of sound data X _l (n), and obtain the frequency domain representation of sound data X _l (k)=[x _l1 (k), ..., x _lP (k)] ^T , where k=1, 2, ..., N/2;

Step 2): X _l (k) that obtains according to described step 1) calculates the modular square η (k) of the coherence function between two microphones in the microphone array according to the following formula (1);

$\mathrm{<span\; class="notranslate">\; \&eta;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\frac{{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

Φ_j，j(k)为第j个传声器的自功率谱，

Among them, k=1, 2,..., N/2; i, j∈[1, 2,..., P]; Φ _{i, j} (k) is the i-th microphone and the j-th microphone cross power spectrum, Φ _{i, i} (k) is the self-power spectrum of the i-th microphone,

Φ _{j, j} (k) is the self-power spectrum of the jth microphone,

Step 3): according to the η (k) that described step 2) obtains, randomly select the square of q coherent function modules;

Step 4): According to the sound data X _l (k) corresponding to the squares of the q coherence function modules selected in the above step 3), the direction of arrival is estimated, where k∈[k ₁ ,...,k _q ], k _q represents the serial number of the frequency band, and [k ₁ ,...,k _q ]∈[1,...,N/2].

The step 1) also includes: saving the continuous L frames of sound data collected so far by the microphone array; The oldest sound data of 1 frame is deleted.

In the step 2), the two microphones are the two microphones with the largest spatial distance in the microphone array.

In the described step 2), the data of consecutive (L-1) frames between the two microphones is used as a data buffer, and the N × (L-1) data obtained according to the N data of the Lth frame currently obtained and the cache are calculated Take the square of the magnitude of the coherence function.

In the step 4), direction of arrival estimation is performed on the L-th frame sound data X _L (k) corresponding to the squares of the q coherence function modules selected in the step 3), wherein, k∈[k ₁ ,...,k _q ], k _q represents the frequency band number, and [k ₁ ,...,k _q ]∈[1,...,N/2].

In the step 3), compare the η(k) values obtained in the step 2), and select the first q η(k) values sequentially from the largest, that is, η(k ₁ )≥η(k ₂ )≥...≥η(k _q )≥η(k _m ), where k _q represents the frequency band number, [k ₁ ,...,k _q ]∈[1,2,...,N/ 2],

In order to achieve the above object, a kind of direction-of-arrival estimation device for broadband sound source is proposed, it is characterized in that, the device comprises: a Fourier transform module, a modulus square module of a coherent function and a direction-of-arrival estimation module;

The Fourier transform module includes: a data selection unit, a data framing unit and a Fourier transform unit; for the data selection unit to select a section of sound data X(n) from the sound data collected by the microphone array )=[x ₁ (n), ..., x _P (n)] ^T , the data framing unit divides the sound data X(n) selected by the data selection unit into L frames evenly, and the The Fourier transform unit performs N-point fast Fourier transform on each frame of sound data X _l (n), and obtains the frequency domain representation of sound data X _l (k)=[x _l1 (k),..., x _lP (k)] ^T ; and the sound data of the frequency domain is output to the modular square module of the coherence function; wherein, P represents the number of microphones in the microphone array, and P≥2; the data x of each frame The length of _l (n) is N, l=1,...,L, the data length of X(n) is N×L; n∈Z ^* ; k=1,2,...,N/2;

The module square module of described coherence function is used for calculating the module square η(k) of the coherence function between two microphones in the microphone array according to the X _l (k) that described Fourier transform module obtains );

$\mathrm{<span\; class="notranslate">\; \&eta;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\frac{{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

Φ_i，i(k)为第i个传声器的自功率谱，

Φ_j，j(k)为第j个传声器的自功率谱，

Among them, k=1, 2,..., N/2; i, j∈[1, 2,..., P]; Φ _{i, j} (k) is the i-th microphone and the j-th microphone cross power spectrum,

Φ _{i, i} (k) is the self-power spectrum of the i-th microphone,

Φ _{j, j} (k) is the self-power spectrum of the jth microphone,

The direction of arrival estimation module is used to randomly select the sound data X _l (k) corresponding to the square of the q coherent function modulus from the modulus square of the coherent function obtained from the modulus square module of the coherent function. Direction of arrival estimation, where k∈[k ₁ ,...,k _q ], k _q represents the frequency band number, and [k ₁ ,...,k _q ]∈[1,...,N/2 ].

In the modular square module of the coherence function, the two microphones are the two microphones with the largest spatial distance in the microphone array.

The device also includes a sound data buffer module; the sound data buffer module is used to save the continuous L frames of sound data collected by the microphone array so far; and after saving the latest 1 frame of sound data received from the microphone array, Deleting the earliest frame of sound data among the saved consecutive L frames of sound data; and outputting the buffered sound data to the Fourier transform module.

The modular square module of the coherence function uses the data of consecutive (L-1) frames between the selected two microphones as a data buffer, and obtains N data of the L-th frame and the N×(L-1) of the buffer according to the current data to find the modulus square of the coherence function.

The direction of arrival estimation module performs direction of arrival estimation on the L-th frame sound data X _L (k) corresponding to the squares of the q coherence function modules obtained by the comparison module, where k∈[k ₁ ,. .., k _q ], k _q represents the serial number of the frequency band, and [k ₁ , ..., k _q ]∈[1, ..., N/2].

；将挑选出的相干函数的模平方输出至所述的波达方向估计模块。The device also includes a comparison module, which is used to select the first q η(k) values sequentially from the largest one, i.e. η(k ₁ ) ≥η(k ₂ )≥...≥η(k _q )≥η(k _m ), where k _q represents the frequency band number, [k ₁ ,...,k _q ]∈[1, 2, .. ..,N/2],

; Output the modulus square of the selected coherence function to the DOA estimation module.

The advantage of the present invention is that the present invention can greatly reduce the calculation amount. At the same time, this simplified algorithm can be applied to any algorithm discussed in the background art section. Taking the AML algorithm as an example, since the AML algorithm has the largest amount of calculation, it can be concluded that the calculation amount is greatly reduced when the orientation effect is equivalent.

Description of drawings

Fig. 1 is a schematic diagram of the spatial structure of a microphone array;

Fig. 2 is a flow chart of a method for estimating the direction of arrival of a broadband sound source proposed by the present invention;

Fig. 3 is a comparison diagram of the results of DOA estimation between the embodiment scheme of the present invention and the full-band scheme under the white noise situation in the simulation test;

Fig. 4 is a comparison diagram of the results of direction of arrival estimation between the embodiment scheme of the present invention and the full-band scheme under the condition of colored noise in the simulation test;

Fig. 5 is a comparison diagram of DOA estimation results between the solution of the embodiment of the present invention and the full-band solution in the field test.

Detailed ways

The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

There are R omni-directional microphones, and the signal received by the pth microphone is:

${\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

根据阵型及目标声源位置的几何关系获得。对于一个远场条件下的R元圆阵，即R个传声器等距排列在一个圆周上，第p个传声器的用采样点数表示的时间延迟为：Among them: n=0,..., L-1, L is the number of sampling points of a frame signal; p=1,..., R, R is the number of microphones in the receiving array; m=1,... ., M, M is the number of sound sources (M<R). is the source signal, is the time delay expressed in number of sampling points. w _p (n) is the noise component in the received signal of the array. Among them, the time delay

It is obtained according to the geometric relationship between the formation and the position of the target sound source. For an R-element circular array under far-field conditions, that is, R microphones are arranged equidistantly on a circle, the time delay of the p-th microphone represented by the number of sampling points is:

${\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}\mathrm{<span\; class="notranslate">\; r</span>}}{\mathrm{<span\; class="notranslate">\; c</span>}}\mathrm{<span\; class="notranslate">\; cos</span>}<span\; class="notranslate">\; (</span>\frac{\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; \&pi;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; R</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

Among them, f _s is the signal sampling frequency, r is the radius of the circular array, c is the wave velocity in the isotropic medium, and θ _m is the far-field azimuth of the mth source.

Perform N-point FFT transformation on the received signal to obtain the frequency domain signal model:

X(k)=D(Θ,k)S ₀ (k)+η(k) (3)

方向矩阵(Steering Matrix)为：D(Θ，k)＝[d⁽¹⁾(k)，...，d^(M)(k)]；方向矢量(Steering Vector)为：

η(k)为噪声频谱。In the formula, k=0,...,N/2; frequency domain array output signal: X(k)=[X ₁ (k),...,X _R (k)] ^T ; M sound sources Spectrum:

The direction matrix (Steering Matrix) is: D(Θ, k)=[d ⁽¹⁾ (k), ..., d ^(M) (k)]; the direction vector (Steering Vector) is:

η(k) is the noise spectrum.

The noise spectrum vector η(k) is the zero-mean complex Gaussian white noise with variance ^Lσ2 . Note that even if w(n) is any independent and identically distributed random variable, according to the central limit theorem, it can be deduced that η(k) obeys the Gaussian distribution.

Definition S(k)＝D(k)S ₀ (k) (4)

Get NR/2-dimensional space-time-frequency vector:

X＝G(θ)+ξ ₁ (5)

where G(θ)=[S(0) ^T , ..., S(N/2) ^T ] ^T , and there are:

${\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; [</span>{\mathrm{<span\; class="notranslate">\; \&xi;\&xi;</span>}}^{\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; ]</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; L</span>}\left[\begin{array}{cccc}{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; R</span>}}& & & \mathrm{<span\; class="notranslate">\; 0</span>}\\ & {\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; R</span>}}& & \\ & & \mathrm{<span\; class="notranslate">\; L</span>}& \\ \mathrm{<span\; class="notranslate">\; 0</span>}& & & {\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; R</span>}}\end{array}\right]<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 6</span>}<span\; class="notranslate">\; )</span>$

where the superscript ^H represents the conjugate transpose of a complex matrix.

Suppose the initial unknown parameter space is:

$\mathrm{<span\; class="notranslate">\; \&Theta;</span>}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; [</span>{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; S</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}^{{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; T</span>}}}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; S</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}^{{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; T</span>}}}<span\; class="notranslate">\; ]</span>}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 7</span>}<span\; class="notranslate">\; )</span>$

源频谱

where the source azimuth

source spectrum

So there is the maximum likelihood function:

$\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&Theta;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; \&pi;</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; NR</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 2</span>}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 2</span>}}}\mathrm{<span\; class="notranslate">\; exp</span>}<span\; class="notranslate">\; \{</span><span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}{<span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; G</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&Theta;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span>}^{\mathrm{<span\; class="notranslate">\; h</span>}}{{\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}}^{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; G</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&Theta;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span><span\; class="notranslate">\; \}</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 8</span>}<span\; class="notranslate">\; )</span>$

Logarithmizing this likelihood function and omitting the constant term gives:

$\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&Theta;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\frac{{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; G</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&Theta;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 9</span>}<span\; class="notranslate">\; )</span>$

The optimization criterion can be obtained:

$\underset{\mathrm{<span\; class="notranslate">\; \&Theta;</span>}}{\mathrm{<span\; class="notranslate">\; max</span>}}\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&Theta;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; \&Theta;</span>}}{\mathrm{<span\; class="notranslate">\; min</span>}}\underset{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\frac{{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; D.</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; S</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>}$

$<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; \&Theta;</span>}}{\mathrm{<span\; class="notranslate">\; min</span>}}\underset{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\frac{{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; D.</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; S</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 10</span>}<span\; class="notranslate">\; )</span>$

This is equivalent to, for all k functions, finding And because σ ² is a constant, f(k) can be abbreviated as:

f(k)＝||X(k)-D(k)S ₀ (k)|| ² (11)

The source signal vector S ₀ (k) that makes f(k) take the minimum value must satisfy So in order to have the smallest residual error for any source location, the source vector estimate is given by:

是方向矩阵D(k)的伪逆矩阵(pseudo-inverse)。定义正交规划那么AML算法的波达方向估计可以通过下式得到。in,

is the pseudo-inverse of the direction matrix D(k). Define an Orthogonal Program Then the direction of arrival estimation of the AML algorithm can be obtained by the following formula.

${\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}<span\; class="notranslate">\; =</span>\underset{{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}}{\mathrm{<span\; class="notranslate">\; max</span>}}\underset{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 13</span>}<span\; class="notranslate">\; )</span>$

Although the sound signal is broadband, it is generally impossible for the signal to be white. It generally has some harmonic frequency structure, that is, the signal is mainly concentrated in certain frequency bands; or because the noise is colored, the signal in certain frequency bands is completely blocked. The noise is annihilated, so when calculating, the information on these frequency bands is redundant and may even affect the final result. Therefore, it is necessary to use a certain method to select certain frequency bands. In the solution of the present invention, the coherence function on the frequency band is used as the criterion for selecting the frequency band.

The coherence function is used to describe the correlation between two sets of data. In the scheme of the present invention, the definition of coherence function is:

${\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; )</span>}{\sqrt{{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; xxx</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; yy</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; )</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 14</span>}<span\; class="notranslate">\; )</span>$

Among them, x and y represent two different sets of data, G _xy (f) is the cross-spectrum of x and y, G _xx (f) and G _yy (f) are the auto-spectrum of x and y respectively. G _xy (f) = E{X(f)Y ^* (f)}, G _xx (f) = E{X(f)X ^* (f)}, G _yy (f) = E{Y(f) Y ^* (f)}.

The magnitude squared coherence function (MSC) is the square of the modulus of the coherence function γ _xy :

${<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\frac{{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; xxx</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; yy</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 15</span>}<span\; class="notranslate">\; )</span>$

From this, it can be seen that the coherence function refers to operations in the frequency domain.

The value of the coherence function is between 0 and 1, namely:

0≤| _γxy | ^2≤1 (16)

When the value is 0, it means that the spectrum on the frequency band x and y is completely uncorrelated, and the value is 1, which means that the spectrum on the frequency band x and y is completely correlated.

由此可见MSC可以用来衡量信号某个频带上的信噪比的大小。当|γ_xy(k)|²较大时，表明第k个频带上有用的信号较多；反之，则说明第k频带上噪声占主要部分。In specific cases there are: SNR

It can be seen that MSC can be used to measure the signal-to-noise ratio in a certain frequency band of the signal. When |γ _xy (k)| ² is large, it indicates that there are more useful signals in the kth frequency band; otherwise, it indicates that the noise in the kth frequency band is the main part.

A schematic diagram of the spatial structure of the microphone array is shown in Figure 1. Use the cross-shaped microphone array shown in Figure 1 to collect wind noise signals, circle 101 represents 12 B&K 4189 standard microphones, and is equipped with a spherical sponge windshield with a diameter of 10cm. 102 denotes brackets for fixing these microphones. The apertures of the three 4-element arrays on the cross array are 0.2m, 0.6m and 1.2m respectively. Using SONY SIR-1000 16-channel synchronous digital recorder, the sampling frequency is 48kHz. During data analysis, low-pass filtering, normalization, and downsampling to 1kHz were used for analysis.

Using the measured data, the spatial correlation of the wind noise of the array is counted. The average correlation coefficient between channels was calculated using the following formula:

$\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; R</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; R</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; R</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; R</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; \&gamma;</span>}}_{\mathrm{<span\; class="notranslate">\; ij</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 17</span>}<span\; class="notranslate">\; )</span>$

为平均相关系数。Wherein, k=0,...,N/2; R is the number of channels, and γ _ij (k) is the cross-correlation coefficient between the i-th channel and the j-th channel,

is the average correlation coefficient.

Table 1 Average correlation coefficient between wind noise channels

Table 1 shows that, unlike the spatial correlation of low-frequency wind noise, the inter-channel correlation coefficient of wind noise at medium and high frequencies is small, and the spatial correlation is weak. The array wind noise can basically be considered independent of the channels.

The solution of the present invention is to filter the sound data according to the coherence coefficient after collecting the sound data, to filter out the sound data of the frequency bands carrying more useful signals, and to discard the sound data of other frequency bands, thus reducing the wave arrival rate. The amount of data required for direction estimation, thereby reducing the amount of computation for subsequent DOA estimation.

A flow chart of a method for estimating the direction of arrival of a broadband sound source proposed by the present invention; as shown in Figure 2, it includes the following steps:

Step 201: Select a piece of sound data x(n)=[x ₁ (n), . . . , x _P (n)] ^T from the sound data collected by a microphone array composed of P microphones, where n Z ^* .

Step 202: the selected sound data x (n) is evenly divided into L segments x _l (n), l=1, ..., L, the length of each segment of data x _l (n) is N, x (n ) data data length is N×L.

Step 203: do N-point Fast Fourier Transform (FFT) to each section of sound data x _l (n), and obtain the frequency domain representation of the sound data as X _l (k)=[X _l1 (k), ..., X _lP (k)] ^T , where k=1, 2, . . . , N/2.

其中，i，j∈[1，2，...，P]，Φ_i，j(k)为第i个传声器和第j个传声器的互功率谱，Φ_i，i(k)和Φ_j，j(k)分别为第i个传声器和第j个传声器的自功率谱：Step 204: Calculate the square of the coherence function modulus (MSC) between two microphones in the microphone array

Among them, i, j∈[1, 2, ..., P], Φ _{i, j} (k) is the cross-power spectrum of the i-th microphone and the j-th microphone, Φ _{i, i} (k) and Φ _{j , j} (k) are the self-power spectra of the i-th microphone and the j-th microphone respectively:

${\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; L</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; li</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ;</span>$ ${\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; L</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; lj</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ;</span>$ ${\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; L</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; li</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; lj</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; .</span>$

Preferably, the two microphones may be the two microphones with the largest spatial distance in the microphone array.

Step 205: Select larger q coherent function modulus square values from the coherent function modulus square η(k), that is, η(k ₁ )≥η(k ₂ )≥...≥η(k _q ) ≥η(k _m ), where k _q is the frequency band number, [k ₁ ,...,k _q ]∈[1,2,...,N/2],

Step 206: Use the selected sound data X _l (k) corresponding to the squares of the q coherence function moduli to perform DOA estimation. Among them, k∈[k ₁ ,...,k _q ]; taking the AML algorithm as an example, the calculation formula of the angle can be known from (18),

${\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}<span\; class="notranslate">\; =</span>\underset{{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}}{\mathrm{<span\; class="notranslate">\; max</span>}}\underset{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}{\overset{{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; q</span>}}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; l</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 18</span>}<span\; class="notranslate">\; )</span>$

Since the calculation of the coherence function requires more data, that is, N×L is larger, it may affect the real-time performance of the method. In the calculation of the coherence function, the data of the previous consecutive (L-1) frames are saved as a data cache, and the MSC is calculated according to the currently obtained N data and the cached N×(L-1) data. And each time a new frame of data is obtained, the data of the previous frame is discarded, and the entire data cache is updated. For example, all the data in the system are x _l (n), l=1,..., L, where x _L (n) is a frame of data currently obtained, and x _l (n), l=1, .. ., L-1 is the previous data. MSC is calculated using x _l (n), l=1, . . . , L, while only x _L (n) is used for angle estimation.

The present invention also proposes a device for estimating the direction of arrival of a broadband sound source, including:

The Fourier transform module is used to perform fast Fourier transform on the sound data collected by the microphone array composed of P microphones to obtain the frequency domain representation of the sound data, and output the frequency domain sound data to the modular square module of the coherence function;

The module square of the coherence function is used to calculate the module square of the coherence function of the sound data between two microphones in the microphone array, and the module square of the calculated coherence function is output to the comparison module;

A comparison module is used to select larger q values from the modulus square of the coherence function;

The direction of arrival estimation module is used to estimate the direction of arrival according to the frequency band of the sound data corresponding to the modulus square of the coherence function selected by the comparison module.

Preferably, the module square of the coherence function is used to calculate the module square of the coherence function of the sound data between the two farthest microphones in the microphone array, and output the calculated module square of the coherence function to the comparison module.

Preferably, the Fourier transform module includes:

The data selection unit is used to select a piece of sound data with a length of N×L from the sound data collected by the microphone array composed of P microphones as the sound data to be processed x(n)=[x ₁ (n),... , x _P (n)] ^T , where n∈Z ^* ;

A data framing unit for dividing the sound data x(n) to be processed selected by the data selection unit into L segments, wherein the sound data of any segment is x _l (n), l=1,...,L , and output the segmented sound data; and

The Fourier transform unit is used to do N-point fast Fourier transform to each section of sound data x _l (n) output by the data framing unit, and obtain the frequency domain expression of the sound data as X _l (k)=[X _l1 (k), . . . , X _lP (k)] ^T , where k=1, 2, . . . , N/2.

其中，i，j∈[1，2，...，P]，Φ_i，j(k)为第i个传声器和第j个传声器的互功率谱，Φ_i，i(k)和Φ_j，j(k)分别为第i个传声器和第j个传声器的自功率谱。Preferably, the module square of the coherence function calculates the module square of the coherence function between the i-th microphone and the j-th microphone of the microphone array as:

Among them, i, j∈[1, 2, ..., P], Φ _{i, j} (k) is the cross-power spectrum of the i-th microphone and the j-th microphone, Φ _{i, i} (k) and Φ _{j , j} (k) are the self-power spectra of the i-th microphone and the j-th microphone, respectively.

Preferably, the comparison module selects the larger q modulus square values of coherent functions from the modulus squares η(k) of N/2 coherent functions, that is, η(k ₁ )≥η(k ₂ )≥ ...≥η(k _q ), k _q represents the serial number of the frequency band, k=1, 2, ..., N/2, [k ₁ , ..., k _q ]∈[1, 2, ... ., N/2].

Preferably, the DOA estimating module performs direction-of-arrival estimation according to the frequency band of the sound data corresponding to the modulus square of the coherence function selected by the comparing module.

Preferably, the device further includes:

The sound data cache module is used to save the continuous L frames of sound data collected so far by the microphone array; and after receiving the latest 1 frame of sound data from the microphone array, the earliest continuous L frames of sound data will be saved 1 frame of sound data is deleted.

The Fourier transform module is used to perform row-fast Fourier transform on the sound data stored in the sound data buffer module to obtain the frequency domain representation of the sound data, and output the frequency domain sound data to the modulus square module of the coherence function.

The technical effect of the solution of the present invention will be described below through the results of simulation experiments and field experiments. In the simulation test, the microphone group array is a circular array, the number of array elements is 4, the radius of the circular array is 0.2m, the sound source is white noise, the number is 1, and the incident angle is 70°. The comparison of the results of the AML algorithm using the full frequency band and the AML algorithm with frequency band selection is shown in Figure 3 and Figure 4, respectively, where Figure 3 is the case of white noise, and Figure 4 is the case of colored noise. The colored noise is generated by passing white noise through a 50-order FIR filter.

It can be seen from Figure 3 that when the background is white noise, the full-band AML algorithm with low SNR is better than the AML algorithm with frequency band selection; the deviations of the two algorithms are similar and tend to be zero. It can be seen from Figure 4 that when the background is colored noise, it is not optimal to use the full-band AML algorithm for low SNRs. It can be seen from the figure that different numbers of frequency bands should be selected under different SNRs. Because the background is colored noise, the signal-to-noise ratio will be large in some frequency bands, and the signal-to-noise ratio will be small in some frequency bands, so selecting too many frequency bands will make these frequency bands with low signal-to-noise ratio also participate in the calculation , so that the deviation of the algorithm increases; under high signal-to-noise ratio, the deviation of the two algorithms also tends to 0.

In the actual field experiment, the array is a circular array, the number of array elements is 4, the radius of the circular array is 0.2m, and the sound source is a tank moving in a circle around the array. The experimental results are shown in Figure 5. The icon "*" is the result of the full-band search, and the icon "o" is the result calculated by taking only the 20 frequency bands with the largest coherence coefficients. It can be seen from the figure that the positioning accuracy of the two is basically the same. The same, but after using the coherence coefficient to select the frequency band, the amount of calculation is 1/15 of the original.

Through the description of the above embodiments, those skilled in the art can clearly understand that the present invention can be realized by means of software plus a necessary hardware platform, and of course all can be implemented by hardware, but in many cases the former is better implementation. Based on this understanding, all or part of the contribution made by the technical solution of the present invention to the background technology can be embodied in the form of software products, and the computer software products can be stored in storage media, such as ROM/RAM, magnetic disks, optical disks, etc. , including several instructions for causing a computer device to execute the methods described in various embodiments or some parts of the embodiments of the present invention. Computer equipment such as: personal computer, server, or network equipment, etc.

Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than limit them. Although the present invention has been described in detail with reference to the embodiments, those skilled in the art should understand that modifications or equivalent replacements to the technical solutions of the present invention do not depart from the spirit and scope of the technical solutions of the present invention, and all of them should be included in the scope of the present invention. within the scope of the claims.

Claims (12)

1. Wave arrival direction estimating method for broadband sound source, the method concrete steps comprise:

Step 1): from the voice data that microphone array gathers, choose one section voice data X (n)=[x ₁(n) ..., x _P(n)] ^T, wherein, P represents the number of microphone in the microphone array, P 〉=2; This section voice data X (n) evenly is divided into the L frame, the data X of every frame _l(n) length is N, and the data length of X (n) is N * L, l=1 ..., L; N ∈ Z ^*

To each frame voice data X _l(n) do the conversion of N point quick Fourier, obtain the frequency domain representation X of voice data _l(k)=[x _L1(k) ..., x _LP(k)] ^T, k=1 wherein, 2 ..., N/2;

Step 2): the X that obtains according to described step 1) _l(k) calculate the mould square η (k) of coherence function between two microphones in the microphone array according to following formula (1);

$\mathrm{\&eta;}\left(k\right)={\left|{\mathrm{\&gamma;}}_{i,j\left(k\right)}\right|}^2=\frac{{\left|{\mathrm{\&Phi;}}_{i,j}\left(k\right)\right|}^2}{{\mathrm{\&Phi;}}_{i,i}\left(k\right){\mathrm{\&Phi;}}_{j,j}\left(k\right)}---\left(1\right)$

Wherein, k=1,2 ..., N/2; I, j ∈ [1,2 ..., P]; Φ _{I, j}(k) be the cross-power spectrum of i microphone and j microphone,

Φ _{I, i}(k) be the auto-power spectrum of i microphone, ${\mathrm{\&Phi;}}_{i,i}\left(k\right)=\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}{|x_{\mathrm{li}}\left(k\right)|}^2;$ Φ _{J, j}(k) be the auto-power spectrum of j microphone, ${\mathrm{\&Phi;}}_{j,j}\left(k\right)=\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}{\left|x_{\mathrm{lj}}\right|}^2;$

Step 3): according to described step 2) η (k) that obtains from the beginning of maximum select successively q coherence function mould square;

Step 4): the voice data X of square correspondence of q the coherence function mould of selecting according to described step 3) _l(k) carry out direction of arrival and estimate, wherein, k ∈ [k ₁..., k _q], k _qExpression band number, and [k ₁..., k _q] ∈ [1 ..., N/2].

2. the Wave arrival direction estimating method for broadband sound source according to claim 1 is characterized in that, described step 1) also comprises: preserve the continuous L frame voice data that described microphone array gathers till current; And preserve from behind the up-to-date 1 frame voice data of described microphone array, with the earliest 1 frame voice data deletion in the continuous L frame voice data of preserving.

3. the Wave arrival direction estimating method for broadband sound source according to claim 1 is characterized in that, described step 2) in, two microphones are two microphones of space length maximum in the microphone array.

4. the Wave arrival direction estimating method for broadband sound source according to claim 2, it is characterized in that, described step 2) in, between two microphones continuously the data of (L-1) frame ask for the mould square of coherence function as data buffer storage according to N data of the current acquisition L frame individual data of the N that deposits * (L-1) of easing up.

5. the Wave arrival direction estimating method for broadband sound source according to claim 4 is characterized in that, in the described step 4), and the L frame voice data X of square correspondence of q the coherence function mould of selecting according to described step 3) _L(k) carry out direction of arrival and estimate, wherein, k ∈ [k ₁..., k _q], k _qExpression band number, and [k ₁..., k _q] ∈ [1 ..., N/2].

6. the Wave arrival direction estimating method for broadband sound source according to claim 1, it is characterized in that in the described step 3), more described step 2) in η (k) size that obtains, select successively a front q η (k) value, i.e. η (k from the beginning of maximum ₁) 〉=η (k ₂) 〉=... 〉=η (k _q) 〉=η (k _m), wherein, k _qThe expression band number, [k ₁..., k _q∈ [1,2 ...., N/2],

7. the direction of arrival estimation unit for broadband sound source is characterized in that this device comprises: the mould of Fourier transform module, coherence function square module and direction of arrival estimation module;

Described Fourier transform module comprises: data selection unit, data framing unit and Fourier transform unit; Described data selection unit is chosen one section voice data X (n)=[x from the voice data that microphone array gathers ₁(n) ..., x _P(n)] ^T, the voice data X (n) that described data selection unit is selected in described data framing unit evenly is divided into the L frame, and described Fourier transform unit is to each frame voice data X _l(n) do the conversion of N point quick Fourier, obtain the frequency domain representation X of voice data _l(k)=[x _L1(k) ..., x _LP(k)] ^TAnd described frequency domain voice data is exported to the mould square module of described coherence function; Wherein, P represents the number of microphone in the microphone array, P 〉=2; The data X of every frame _l(n) length is N, l=1 ..., L, the data length of X (n) is N * L; N ∈ Z ^*K=1,2 ..., N/2;

The mould of described coherence function square module is used for the X that obtains according to described Fourier transform module _l(k) calculate the mould square η (k) of coherence function between two microphones in the microphone array according to formula (2);

$\mathrm{\&eta;}\left(k\right)={\left|{\mathrm{\&gamma;}}_{i,j\left(k\right)}\right|}^2=\frac{{\left|{\mathrm{\&Phi;}}_{i,j}\left(k\right)\right|}^2}{{\mathrm{\&Phi;}}_{i,i}\left(k\right){\mathrm{\&Phi;}}_{j,j}\left(k\right)}---\left(2\right)$

Wherein, k=1,2 ..., N/2; I, j ∈ [1,2 ..., P]; Φ _{I, j}(k) be the cross-power spectrum of i microphone and j microphone,

Φ _{I, i}(k) be the auto-power spectrum of i microphone,

${\mathrm{\&Phi;}}_{i,i}\left(k\right)=\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}|{{X_{\mathrm{lj}}\left(k\right)|}_{;}}^2$ Φ _{J, j}(k) be the auto-power spectrum of j microphone, ${\mathrm{\&Phi;}}_{j,j}\left(k\right)=\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}{\left|X_{\mathrm{lj}}\left(k\right)\right|}^2;$

Described direction of arrival estimation module, the mould square that is used for the coherence function that obtains according to the mould square module from described coherence function is selected the voice data X of square correspondence of q coherence function mould successively from the beginning of maximum _l(k) carry out direction of arrival and estimate, wherein, k ∈ [k ₁..., k _q], k _qExpression band number, and [k ₁..., k _q] ∈ [1 ..., N/2].

8. the direction of arrival estimation unit for broadband sound source according to claim 7 is characterized in that, in the mould of the described coherence function square module, two microphones are two microphones of space length maximum in the microphone array.

9. the direction of arrival estimation unit for broadband sound source according to claim 7 is characterized in that described device also comprises the voice data cache module;

Described voice data cache module be used for to be preserved the continuous L frame voice data that microphone array gathers till current; And preserve and to receive from behind the up-to-date 1 frame voice data of microphone array, with the earliest 1 frame voice data deletion in the continuous L frame voice data of preserving; And export the voice data of buffer memory to described Fourier transform module.

10. the direction of arrival estimation unit for broadband sound source according to claim 9, it is characterized in that, between two microphones that the mould of described coherence function square module will be chosen continuously the data of (L-1) frame ask for the mould square of coherence function as data buffer storage according to N data of the current acquisition L frame individual data of the N that deposits * (L-1) of easing up.

11. the direction of arrival estimation unit for broadband sound source according to claim 10 is characterized in that, the L frame voice data X of square correspondence of q the coherence function mould that described direction of arrival estimation module obtains according to comparison module _L(k) carry out direction of arrival and estimate, wherein, k ∈ [k ₁..., k _q].

12. the direction of arrival estimation unit for broadband sound source according to claim 7, it is characterized in that, this device also comprises a comparison module, be used for η (k) size that the mould square module according to more described coherence function obtains, select successively a front q η (k) value, i.e. η (k from the beginning of maximum ₁) 〉=η (k ₂) 〉=... 〉=η (k _q) 〉=η (k _m), wherein, k _qThe expression band number, [k ₁..., k _q] ∈ [1,2 ...., N/2],

Square export the mould of the coherence function picked out to described direction of arrival estimation module.

Images