Broadband quasi-stationary signal direction-of-arrival estimation method and device based on frequency pair

Abstract

The invention provides a broadband quasi-stationary signal direction of arrival estimation method and device based on frequency pairs. The method comprises the following steps: constructing a second order difference composite array on different frequencies in an interested frequency range according to a far-field broadband quasi-stationary signal received by a physical array; the interested frequency is combined with itself, every two frequencies are mutually combined, a plurality of frequency pairs are obtained, the difference synthesis array of the second-order difference synthesis array of each frequency pair is calculated, and a double-difference synthesis array corresponding to the frequency pair is obtained; the double-difference composite array corresponding to the self-combination frequency is called as a self-double-difference composite array, and the double-difference composite array corresponding to the mutual-combination frequency is called as a mutual-double-difference composite array; and solving the double-difference composite array by using group sparsity to realize direction of arrival estimation. Compared with the existing similar method, the DOA estimation performance of the quasi-stationary signal can be further improved by using the time domain characteristic and the frequency domain characteristic of the quasi-stationary signal at the same time.

Description

Technical Field

[0001] The present invention relates to the technical field of array signal processing, and in particular to a method and device for estimating the direction of arrival of a broadband quasi-stationary signal based on a frequency pair. Background Art

[0002] Underdetermined Direction of Arrival (DOA) estimation is an important problem in array signal processing, where the number of targets to be estimated exceeds the number of physical array elements. For traditional DOA estimation methods, such as the Multiple Signal Classification (MUSIC) algorithm and the Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) algorithm, a uniform linear array (ULA) with M elements can resolve at most M-1 targets. An important idea to solve underdetermined DOA estimation is to generate a virtual array with more virtual elements (far more than the number of physical elements). Based on this idea, high-order cumulants, sparse array structures, coprime frequencies, and quasi-stationary signals can all be used to create virtual arrays.

[0003] Quasi-stationary signal (QSS) is a special type of signal. Its statistical characteristics remain basically unchanged within a time frame, but change between different time frames. Common speech and audio signals in life can be regarded as quasi-stationary signals. In the past decade, scholars have studied the underdetermined DOA estimation of quasi-stationary signals. In the 2010 issue of IEEE Transactions on Signal Processing, Vol. 58, No. 4, pp. 2168 to 2180, Wing-Kin Ma proposed the KR-MUSIC method in "DOAEstimation of Quasi-Stationary SignalsWith Less SensorsThan Sources and Unknown Spatial Noise Covariance:AKhatri–Rao Subspace Approach". This method constructs a virtual array containing 2M-1 virtual elements by calculating the Khatri-Rao product of the local covariance matrix, where M represents the number of elements in the ULA used. Therefore, this method can estimate 2M-2 targets at most, achieving underdetermined DOA estimation. Considering that sparse arrays can further increase the number of virtual array elements, Yuexian Wang published “Sparsity-aware DOA estimation of quasi-stationary signals using nested arrays” on pages 87 to 98 of Volume 144 of Signal Processing in 2018. A two-level nested array (TL-NA) was used, and the KR-MUSIC method was used for DOA estimation. It can estimate up to 2 (M1M2+M2-1) targets (M1 and M2 are the number of array elements in the two sub-arrays of TL-NA, respectively).

[0004] In addition to utilizing time domain features, effectively utilizing information at different frequencies is crucial to improving DOA estimation performance. In “DOA estimation exploiting a uniform linear array with multiple co-prime frequencies” published by Si Qin in Volume 130, pages 37 to 46 of Signal Processing, a broadband DOA estimation algorithm based on coprime frequency pairs was proposed. This algorithm constructs a difference co-array with higher degrees of freedom by treating two array models with coprime frequencies as subarrays of a coprime array. This method simultaneously utilizes the autocorrelation information within the frequency and the co-correlation information between frequencies, and can achieve O(M 2 ), where M represents the number of physical elements of the ULA used. It is worth noting that this method is not designed for quasi-stationary signals, so the time domain characteristics of quasi-stationary signals are not used. Summary of the invention

[0005] In view of this, the present invention provides a broadband quasi-stationary signal DOA estimation method based on frequency pairs. Considering that the existing quasi-stationary signal DOA estimation techniques only utilize the time domain characteristics or frequency domain characteristics of the quasi-stationary signal alone, and the improvement of the DOA estimation performance is limited, the present invention simultaneously utilizes the time domain characteristics and the frequency domain characteristics to further improve the DOA estimation performance of the quasi-stationary signal.

[0006] In order to solve the above technical problems, the present invention is implemented as follows.

[0007] A method for estimating the direction of arrival of a broadband quasi-stationary signal based on frequency pairs, comprising:

[0008] Step 1: Based on the far-field broadband quasi-stationary signal received by the physical array, construct a second-order difference synthesis array corresponding to different frequencies in the frequency range of interest;

[0009] Step 2: Combine the frequency of interest with itself, and combine the frequencies with each other to obtain multiple frequency pairs; calculate the difference synthesis array of the second-order difference synthesis array of each frequency pair to obtain the double difference synthesis array model corresponding to the frequency pair; the double difference synthesis array model corresponding to the self-combination frequency is called the self-double difference synthesis array, and the double difference synthesis array model corresponding to the mutual combination frequency is called the mutual-double difference synthesis array;

[0010] Step 3: Combine the self-double difference synthesis array and the mutual-double difference synthesis array to construct an optimization problem, and use group sparse solution to achieve wave direction estimation.

[0011] Preferably, step 2 comprises:

[0012] Assume that K far-field targets are distributed at the incident angle θ k , k = 1, 2, ..., K; each source signal contains I frequencies of interest f i ,i=1,2,...,I;

[0013] In step 1, let the received signal be at frequency f i The autocovariance matrix of the pth frame is R (i,i) [p]; quantize the matrix to obtain the equivalent received signal z corresponding to the second-order difference synthesis matrix (i,i) [p];

[0014]

[0015] Among them, vec(·) is a vectorization operator, is the steering matrix corresponding to the second-order difference synthesis matrix, where a i (θ k ) represents the steering vector of the k-th source signal; and is the equivalent source signal vector; is the Kronecker product, * indicates taking the conjugate, I M represents the M×M identity matrix; represents the noise vector n i The autocovariance matrix of [l];

[0016] Calculate the frequency f i Up (i,i) [p] is the expectation between frames and is converted into a zero-mean variable

[0017]

[0018] in is the corresponding source signal vector, and E{·} represents taking the mean of the input variables;

[0019] Similarly, the frequency f j Where for

[0020] Calculate z (i,i) [p] and z (j,j) The cross-covariance matrix and auto-covariance matrix of [p] are uniformly expressed as R (i,j) , then R (i,j) Vectorize and get a new virtual array model:

[0021] y (i,j) =vec(R (i,j) )=B(i,j) u (i,j) (1)

[0022] Among them, y (i,j) is the array receiving signal vector; B (i,j) =[b(θ1),b(θ2),…,b(θ K )] is the equivalent steering matrix, where u (i,j) is the equivalent source signal vector of a single snapshot;

[0023] When i≠j, R (i,j) is the cross-covariance matrix, and the virtual array model represents the cross-double difference matrix; when i=j, R (i,j) is the autocovariance matrix, and the virtual array model represents the self-double difference composite array.

[0024] Preferably, the step 3 comprises:

[0025] Let K g Indicates the number of angle search grids, represents the angle search grid in sparse recovery; in the compressed sensing framework, formula (1) is converted into a sparse representation:

[0026]

[0027] in As an overcomplete representation of the equivalent steering matrix, the sparse vector Each element in represents the direction of arrival result on the corresponding angle search grid;

[0028] Based on group sparsity, multiple frequency pairs (f i ,f j ) is expressed as:

[0029]

[0030] in, 1≤i<j≤I, Represents the corresponding predefined K g The DOA estimation result of the angle search grid, ∈ is the defined error limit;

[0031] Solve the optimization problem and obtain K g ×1 vector, The angle in the angle search grid corresponding to the mid-peak value is the target wave direction.

[0032] The present invention also provides a broadband quasi-stationary signal direction of arrival estimation device based on frequency pairs, the device comprising:

[0033] A second-order difference synthesis array construction module is used to construct a second-order difference synthesis array corresponding to different frequencies in a frequency range of interest according to a far-field broadband quasi-stationary signal received by a physical array;

[0034] A double difference synthesis array construction module is used to combine the frequency of interest with itself, and to combine two frequencies with each other, to obtain multiple frequency pairs; calculate the difference synthesis array of the second-order difference synthesis array of each frequency pair, to obtain a double difference synthesis array model corresponding to the frequency pair; the double difference synthesis array model corresponding to the self-combination frequency is called a self-double difference synthesis array, and the double difference synthesis array model corresponding to the mutual combination frequency is called a mutual-double difference synthesis array;

[0035] The solution module is used to construct an optimization problem by combining the self-double difference synthesis array and the mutual-double difference synthesis array, and realize the direction of arrival estimation by using group sparse solution.

[0036] Beneficial effects:

[0037] The present invention first constructs a second-order difference synthesis array at different frequencies in the frequency range of interest. Due to the quasi-stationary characteristics of the signal, the signal can be processed in segments. At this time, the second-order difference synthesis array obtained is multi-fast, and the difference synthesis array can be calculated again. Therefore, the present invention further calculates the difference synthesis array of the second-order difference synthesis array in the frequency domain, which is called a double difference synthesis array. The double difference synthesis array has a higher array degree of freedom after two difference synthesis array calculations. The higher degree of freedom can improve the DOA estimation performance of the quasi-stationary signal. Therefore, the degree of freedom of the double difference synthesis array finally obtained by the present invention is greatly improved compared to the second-order difference synthesis array, which makes the DOA estimation result also have a higher estimation accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 It is a principle diagram of a broadband quasi-stationary signal direction of arrival estimation method based on frequency pairs of the present invention;

[0039] Figure 2 It is a flow chart of a method for estimating the direction of arrival of a broadband quasi-stationary signal based on frequency pairs according to the present invention;

[0040] Figure 3 The DOA estimation result of an example provided by the present invention;

[0041] Figure 4 It is a structural diagram of a broadband quasi-stationary signal direction of arrival estimation device based on frequency pairs according to the present invention. DETAILED DESCRIPTION

[0042] The present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0043] The present invention provides an underdetermined DOA estimation scheme for broadband quasi-stationary signals based on frequency pairs. The scheme constructs a second-order difference synthesis array model at different frequencies in a frequency range of interest according to the far-field broadband quasi-stationary signal received by a physical array. Due to the quasi-stationary characteristics of the signal, the signal can be processed in segments. At this time, the obtained second-order difference synthesis array model is multi-fast, and the difference synthesis array can be calculated again. Subsequently, the difference synthesis array of the second-order difference synthesis array is calculated in the frequency domain, and a virtual array named Dual Difference Co-Array (DDCA) is obtained, which is located in the space-spectral domain (Spatio-Spectral Domain) and has greatly improved array freedom. Finally, group sparsity is used to solve the dual difference synthesis array model to achieve wave direction estimation. Since the dual difference synthesis array has a higher array freedom, the DOA estimation performance of the quasi-stationary signal can be improved.

[0044] See also Figure 1 and Figure 2 The present invention provides a method for estimating the direction of arrival of a broadband quasi-stationary signal based on a frequency pair, which comprises the following steps:

[0045] Step 1: Establish an array receiving signal model.

[0046] In this step, the array element spacing is calculated in advance according to the target frequency of interest. A far-field broadband quasi-stationary signal is incident on the array, and the L-point discrete Fourier transform of the received signal is calculated to obtain a frequency domain multi-frequency sub-band model.

[0047] Step 2: Generate a signal model based on the second-order difference synthesis array.

[0048] In this step, the different frequency components of the multi-frequency sub-band model are processed in sections, and the covariance matrix is ​​calculated and vectorized on each section to obtain a multi-fast second-order difference synthesis array model at different frequencies.

[0049] Step 3: Generate a signal model based on a double difference synthesis array.

[0050] In this step, the frequency of interest is combined with itself, and the frequencies are combined with each other to obtain multiple frequency pairs, and the difference synthesis matrix of the second-order difference synthesis matrix of each frequency pair is calculated to obtain a double difference synthesis matrix model corresponding to the frequency pair. The double difference synthesis matrix has a higher degree of freedom than the second-order difference synthesis matrix.

[0051] The double difference synthesis array corresponding to the self-combination frequency is called the self-double difference synthesis array. The self-combination frequency here refers to the combination of a frequency and itself, for example (f i ,f i); The double difference synthesis array corresponding to the mutual combination frequency is called the mutual-double difference synthesis array, where the mutual combination frequency refers to the combination of different frequencies, for example (f i ,f j ),i≠j.

[0052] Step 4: Direction of Arrival estimation based on group sparseness.

[0053] In this step, the optimization problem is constructed by combining the self-double difference synthesis array and the mutual-double difference synthesis array, and the direction of arrival estimation is realized by using the group sparse solution.

[0054] The above process is described in detail below.

[0055] Step 1. Build array receiving signal model

[0056] Consider a random linear array with M physical elements and a unit element spacing of d. Define the set of element positions of the linear array as for:

[0057]

[0058] Where P m =η m d represents the position of the mth array element, where η m is an arbitrary integer.

[0059] Assume that K far-field targets are distributed at the incident angle θ k (k=1,2,...,K), each source signal contains I frequencies of interest f i ,i=1,2,...,I, these frequencies satisfy:

[0060] f i =ζ i f0 (2)

[0061] Where ζ i is an integer, and f0 is the unit frequency, i.e., the equivalent operating frequency, with a wavelength of λ0 = 2d to avoid spatial aliasing. In actual use, the array element spacing d is designed according to the unit frequency. The array element spacing is related to the frequency of interest. For example, if you are interested in f0, you can design the array element spacing of d = c / (2f0).

[0062] Let x i,m [l], l = 0, 1, 2, ... means the received signal of the mth sensor is converted to digital at frequency f i The weight on i,k [l] represents the kth source signal at frequency f i The frequency f iThe source signal vector and the received signal vector on are expressed as follows:

[0063]

[0064] in[·] T is the transpose operator.

[0065] The array output is represented as

[0066]

[0067] Among them, n i [l] is a noise vector, each element of which follows a zero-mean Gaussian distribution. The noise is uncorrelated with the source signal. i =[a i (θ1),a i (θ2),…,a i (θ K )] is the array steering matrix, whose kth column represents the steering vector of the kth source signal, where λ i is the frequency f i wavelength.

[0068] The broadband array model involved in the present invention uses a broadband quasi-stationary signal as the source signal and is required to meet the following assumptions:

[0069] 1. Each source signal s i,k [l] are all generalized second-order quasi-stationary signals, and the frame length L is random. The second-order statistics of the p-th frame (p∈[1,P]) are:

[0070]

[0071] Among them, E{·} means taking the mean of the input variables.

[0072] In addition, x i The local autocovariance matrix of [l] at the pth frame (l∈[(p-1)L,pL-1]) is:

[0073]

[0074] in{·} H represents the conjugate transpose operator, I M represents the M×M identity matrix, Indicates n i The autocovariance matrix of [l].

[0075] In addition D (i,i) [p] is defined as:

[0076]

[0077] where diag{·} returns a diagonal matrix with the input vector as its diagonal elements.

[0078] 2. At each frequency, k=1,2,…,K, is wide-sense stationary and uncorrelated with each other. Definition

[0079]

[0080] where var{·} returns the variance of the input values.

[0081] 3. For different sources (k1≠k2), and i,j∈[1,2,…,I], we have

[0082]

[0083] in(·) * Returns the conjugate of the input value.

[0084] 4. For each source, assume that its variances at different frequencies show partial correlation. They satisfy

[0085]

[0086] Step 2. Generate a signal model based on the second-order difference synthesis array.

[0087] Considering the time domain characteristics of the quasi-stationary signal, we first calculate the autocovariance matrix R in each frame according to formula (6): (i,i) [p]. Then, the R of the pth frame (i,i) [p] Vectorization yields:

[0088]

[0089] The vectorization operator vec(·) converts the input matrix into a vector by column stacking. (i,i) [p] is the equivalent received signal corresponding to the second-order difference synthesis array. is the steering matrix of the equivalent array, where and is the equivalent source signal vector. is the Kronecker product, and

[0090] In addition, the elements of the second-order difference synthesis array are located at:

[0091]

[0092] in

[0093] Step 3. Generate a signal model based on a double difference synthesis array.

[0094] Next, calculate z (i,i) [p] is the expectation between frames and is converted into a zero-mean process

[0095]

[0096] in is the corresponding source signal vector.

[0097] Likewise, the frequency f j Where Can be obtained in the same way:

[0098]

[0099] z (i,i) [p] and z (j,j) [p] represents the position of the array element and The equivalent received signal vector of the virtual array of . According to the assumption, their corresponding source signal vectors and There is partial correlation. Calculate z (i,i) [p] and z (j,j) The cross-covariance matrix and autocovariance matrix of [p] are then vectorized to obtain a virtual array with higher degrees of freedom.

[0100] Specifically, the covariance matrix R is calculated by (i,j) (i,j∈[1,2,…,I]):

[0101]

[0102] R (i,j) Vectorized to get:

[0103] y (i,j) =vec(R (i,j) )=B (i,j) u (i,j) , (16)

[0104] Equation (16) describes a virtual array model (double difference synthesis array model), where y (i,j) is the received signal vector. (i,j) =[b(θ1),b(θ2),…,b(θ K )] is an equivalent steering matrix, where its kth column in addition, is the equivalent source signal vector of a single snapshot.

[0105] Furthermore, the set of virtual array elements associated with the double difference synthesis array is represented as:

[0106]

[0107] When i≠j, the virtual array represents a mutual-double-difference synthesis array, which utilizes the information across frequencies in the frequency pair; when i=j, it represents a self-double-difference synthesis array, which only utilizes the information of a single frequency in the frequency pair.

[0108] Step 4. Combine the self-double difference synthesis array and the mutual-double difference synthesis array to construct the optimization problem, and use group sparsity to estimate the direction of arrival.

[0109] Let K g Indicates the number of angle search grids, represents the angle search grid in sparse recovery. Under the compressed sensing framework, the expression in (16) can be converted into a sparse representation, that is,

[0110]

[0111] in As an overcomplete representation of the equivalent steering matrix, the sparse vector Each element in represents the direction of arrival result on the corresponding angle search grid.

[0112] Based on group sparsity, a single frequency pair (f i ,f j ) can be expressed as:

[0113]

[0114] in and Return a K g ×1 vector, where the i-th element is The i-th row vector of norm. In addition, ∈ is a user-defined error limit.

[0115] When there are multiple frequency pairs of interest, equation (19) can be expressed as

[0116]

[0117] in 1≤i<j≤I, Represents the corresponding predefined K g The DOA estimation results for the search grid.

[0118] Solve the optimization problem and obtain It is a K g ×1 vector, The angle in the angle search grid corresponding to the mid-peak value is the target wave direction.

[0119] The broadband quasi-stationary signal direction of arrival estimation method based on frequency pairs provided by the present invention can greatly improve the degree of freedom of the array. The degree of freedom is analyzed below.

[0120] For a sparse linear array of any shape, assume that the continuous part of its second-order difference composite array is Where φ0 is an arbitrary integer. Then, at frequency f i and f j The corresponding set of second-order difference composite matrices is

[0121]

[0122] Afterwards, with and The subset of the relevant cross-double difference matrix is

[0123]

[0124] in For simplicity, the symbol d is omitted in the following text, and only

[0125] The elements in are symmetric about the origin, that is, if but Assumptions Then ψ2=-ψ1 satisfies

[0126] ψ2=-ψ1=ζ i (-μ1)-ζ j (-μ2) = ζ i μ3-ζ j μ4, (23)

[0127] where μ3=-μ1, μ4=-μ2, and This means

[0128] Use N u express The number of elements in . Taking redundancy into account, The maximum number of elements in is N t =(2φ0+1) 2 . To calculate the number of elements N excluding redundancy u , first define the weight function. The weight function ω(ψ c ) is a specific interval The number of virtual array element pairs:

[0129]

[0130] in |·|Returns the cardinality of the input set, and 1≤ζ i ≤ζ j .

[0131] Define any element The redundancy is:

[0132] ρ(ψ c )=ω(ψ c )-1. (25)

[0133] Then put The redundancy of all elements in is added up, and the result is defined as The number of non-repeating elements in can be obtained by the following formula:

[0134] N u =N t -N o =(2φ0+1) 2 -N o . (26)

[0135] If you get N o , then N can also be calculated u .

[0136] Assume that there are any two elements where μ1,μ2,μ3,μ4∈[-φ0,φ0]. If And μ1>μ3,μ2>μ4 holds, then N o =|(μ1,μ2)|.

[0137] Because i and j are mutually prime, the above requirement is equivalent to

[0138]

[0139] Where k is a positive integer.

[0140] From formula (27), we can get

[0141]

[0142] in is the maximum value of k, and Represents a floor function that returns the largest integer not greater than the input value.

[0143] N o The derivation of is divided into the following two cases.

[0144] First consider 1≤ζ j <2φ0+1, at this time k max ≥1. Since μ1,μ2∈[-φ0,φ0], equation (28) is transformed into

[0145]

[0146] For any k∈[1,k max ],have

[0147]

[0148] The numbers of μ1 and μ2 are given by and OK. There are 2φ0-ζ j +1 integer, There are 2φ0-ζ i +1 integer. So N o =(2φ0-ζ j +1)(2φ0-ζ i +1).

[0149] Next, when j =2φ0+1,k max =0, which contradicts k>0, meaning that equation (27) has no solution. Then

[0150] N o =0=(2φ0-ζ j +1)(2φ0-ζ i +1). (31)

[0151] In short, when 1≤ζ j ≤2φ0+1,N o =(2φ0-ζ j +1)(2φ0-ζ i +1), and

[0152]

[0153] When j =2φ0+1, we can get N u The maximum value N u max .

[0154] When j >2φ0+1,k max =0, which contradicts k>0. Therefore, equation (27) has no solution, and N o =0, which means N u =Nt =(2φ0+1) 2 .

[0155] Next, we will take a TL-NA with M (odd number) array elements. The two sub-arrays have (M-1) / 2 and (M+1) / 2 array elements respectively. We can deduce

[0156]

[0157] Assumption i and j (1≤ζ i ≤ζ j ≤2φ0+1) are relatively prime. Substituting φ0 into equation (32) yields:

[0158]

[0159] Therefore, the degrees of freedom of the mutual-double difference matrix can reach This exceeds the degrees of freedom provided by the above methods.

[0160] Finally, because the self-double-differenced array and the cross-double-differenced array can be combined by group sparsity, the virtual element positions are located at because The number of elements in is greater than N u More can bring higher degrees of freedom, so N u It can be regarded as the lower limit of the degrees of freedom of this method. The actual degrees of freedom will not be lower than this value, which is better than the existing methods.

[0161] In short, by using the quasi-stationary characteristics of the signal, a multi-fast second-order difference synthesis array model can be generated at different frequencies. Based on this model, the difference synthesis array can be calculated again between frequencies. Generally speaking, calculating the difference synthesis array will generate more degrees of freedom. As proved above, the degree of freedom of the final virtual array is greatly improved compared to the second-order difference synthesis array.

[0162] An example is given below to introduce the scheme and effect of the present invention.

[0163] (1) System model design

[0164] In this example, a 4-element microphone array is used as the receiving array. The array structure adopts a nested array, and the array element positions are located at {1,2,3,6}d, where d is the unit array element spacing. The signal is stored on the computer using a data acquisition card (PXUA216MB-DL2-M). Let the system sampling rate f s =44.1KHz, the number of discrete Fourier transform time points N FFT =512.

[0165] There are K = 9 far-field targets distributed at the incident angle θ k =-60+15k (k = 1, 2, ..., 9). Set the target source signal s(t) to a broadband quasi-stationary signal: the source signal in each frame obeys a zero-mean Gaussian distribution (with a variance of σ s 2 ), σ of different frames s It follows a uniform distribution within [0,1], and the number of sampling points in each frame follows [18N FFT ,42N FFT ] are evenly distributed within the source signal, and the total length of the source signal is 29s.

[0166] The signal received by the nth array element can be expressed as:

[0167]

[0168] The discrete form of the received signal column vector is shown as follows:

[0169] x[l]=[x1[l],x2[l],x3[l],x6[l]] T . (36)

[0170] The received signal of each array element is divided into non-overlapping groups of length I = 512, and the N of each group of signals is calculated. FFT =512-point discrete Fourier transform, the number of frequency domain data sets L is calculated based on the signal length f =2580. The lth frequency band of the ith frequency band can be obtained f The frequency domain observation signal vector of the group:

[0171] X[i,l f ]=[X1[i,l f ],X2[i,l f ],X3[i,l f ],X6[i,l f ]] T , (37)

[0172] In the formula,

[0173]

[0174] Among them, l f =0,1,…,2579 and l=0,1,…,511.

[0175] The frequency domain model of the array output signal can be expressed as:

[0176]

[0177] In the formula, A(i,θ)=[a(i,θ1),...,a(i,θ K )] indicates the frequency is f i The steering matrix of the i-th frequency band corresponds to the incident angle θ k The steering vector a(i,θ k )for:

[0178]

[0179] The frequency spacing between adjacent frequency sub-bands is f Δ =f s / 512. Select unit frequency f0 = 33f Δ , let the array element spacing d be equal to 1 / 2 of the wavelength of the frequency, then the steering vector at frequency f0 can be expressed as:

[0180]

[0181] Where φ k =2πsin(θ k ).

[0182] Frequency i f0 and ζ j f0, The steering vectors are:

[0183]

[0184] For simplicity, the above two steering vectors are defined as: and The corresponding signal is defined as and

[0185] (2) Double difference synthesis array model for a single frequency pair

[0186] The frame length of the fixed signal is L = 30. First, calculate the p-th frame (l∈[(p-1)L,pL-1]) The local autocovariance matrix of :

[0187]

[0188] Then, R of the pth frame (i,i) [p] Vectorization yields:

[0189]

[0190] The equivalent steering matrix is the steering matrix of the equivalent array, where and is the equivalent source signal vector.

[0191] Next, calculate z (i,i) [p] is the expectation between frames and is converted into a zero-mean process

[0192]

[0193] in is the corresponding source signal vector.

[0194] Likewise, the frequency f j Where Can be obtained in the same way:

[0195]

[0196] Next, the covariance matrix R is calculated as follows (i,j) :

[0197]

[0198] R (i,j) Vectorization gives the frequency ζ i f0 and ζ j Double difference synthesis array model under f0:

[0199] y (i,j) =vec(R (i,j) )=B (i,j) u (i,j) , (48)

[0200] Among them B (i,j) =[b (i,j) (θ1),b (i,j) (θ2),…,b (i,j) (θ K )] is the equivalent steering matrix, where

[0201] (3) Direction of Arrival Estimation Based on Group Sparsity for a Single Frequency Pair

[0202] If we use the frequency pair 6f0 and 7f0, that is, ζ i =6,ζ j =7. In the compressed sensing framework, first, a total K is generated in the full incident angle space from -90° to 90° with a step of 0.1°. g = 1801 search grids θ g,0 ,...,θ g,1800 . Construct an overcomplete representation of each equivalent steering matrix in the virtual array model:

[0203]

[0204] by Represents the 1801×1 dimensional column signal vector to be estimated, then the underdetermined DOA estimation based on group sparsity for a single frequency pair can be expressed as:

[0205]

[0206] in and Returns a 1801×1-dimensional column vector where the i-th element is The i-th row vector of Norm. The angle in the angle search grid corresponding to the peak is the target wave direction. In addition, ∈ is the user-defined error limit.

[0207] (4) Direction of Arrival Estimation Based on Joint Sparseness for Multiple Frequency Pairs

[0208] If the two frequency pairs 6f0, 7f0 and 2f0, 7f0 are used at the same time in the example, the optimization problem is transformed into:

[0209]

[0210] in represents the direction of arrival estimation result corresponding to the predefined 1801 search grids. The error limit ∈ is determined by the trial-and-error method. The azimuth estimation result is as follows Figure 3 , where the dotted line represents the real target's orientation, and the solid line represents the orientation spectrum estimated by this algorithm. It can be seen that this algorithm accurately estimates the orientation of all targets.

[0211] Based on the above method, the present invention also provides a broadband quasi-stationary signal direction of arrival estimation device based on frequency pairs, such as Figure 4 As shown, the device comprises:

[0212] The second-order difference synthesis array construction module is used to construct the second-order difference synthesis array corresponding to different frequencies in the frequency range of interest according to the far-field broadband quasi-stationary signal received by the physical array.

[0213] The double difference synthesis array construction module is used to perform self-combination and mutual combination of the frequencies of interest in pairs, obtain multiple frequency pairs, calculate the difference synthesis array of the second-order difference synthesis array of each frequency pair, and obtain the double difference synthesis array corresponding to the frequency pair; the double difference synthesis array corresponding to the self-combination frequency is called the self-double difference synthesis array, and the double difference synthesis array corresponding to the mutual combination frequency is called the mutual-double difference synthesis array. The way in which the double difference synthesis array construction module constructs the double difference synthesis array is the same as step 3 above, which will not be repeated here.

[0214] The solution module is used to solve the double difference matrix by group sparseness to realize the direction of arrival estimation. The solution of the solution module is the same as step 4 above, which will not be repeated here.

[0215] In summary, the present invention provides an underdetermined DOA estimation method and device for broadband quasi-stationary signals, which greatly improves the degree of freedom in DOA estimation by designing an array structure and selecting a suitable frequency pair, and is superior to similar methods. Therefore, the number of estimable targets of this method is better than other methods, and it has higher estimation accuracy.

[0216] The above specific embodiments only describe the design principle of the present invention. The shapes and names of the components in the description may be different and are not limited. Therefore, those skilled in the art in the field of the present invention may modify or replace the technical solutions recorded in the above embodiments; and these modifications and replacements do not deviate from the creative purpose and technical solutions of the present invention and should all fall within the protection scope of the present invention.

Claims

1. A method for estimating the direction of arrival of a broadband quasi-stationary signal based on frequency pairs, characterized in that: include: Step 1: Based on the far-field broadband quasi-stationary signal received by the physical array, construct a second-order difference synthesis array corresponding to different frequencies in the frequency range of interest; Step 2: Combine the frequency of interest with itself, and combine the frequencies with each other to obtain multiple frequency pairs; calculate the difference synthesis array of the second-order difference synthesis array of each frequency pair to obtain the double difference synthesis array model corresponding to the frequency pair; the double difference synthesis array corresponding to the self-combination frequency is called the self-double difference synthesis array, and the double difference synthesis array corresponding to the mutual combination frequency is called the mutual-double difference synthesis array; Step 3: Combine the self-double difference synthesis array and the mutual-double difference synthesis array to construct an optimization problem, and use group sparse solution to achieve wave direction estimation.

2. The method according to claim 1, characterized in that The step 2 comprises: Assume that K far-field targets are distributed at the incident angle θ k , k = 1, 2, ..., K; each source signal contains I frequencies of interest f i ,i=1,2,...,I; In step 1, let the received signal be at frequency f i The autocovariance matrix of the pth frame is R (i,i) [p]; quantize the matrix to obtain the equivalent received signal z corresponding to the second-order difference synthesis matrix (i,i) [p]; Among them, vec(·) is a vectorization operator, is the steering matrix corresponding to the second-order difference synthesis matrix, where a i (θ k ) represents the steering vector of the k-th source signal; and is the equivalent source signal vector; is the Kronecker product, * indicates taking the conjugate, I M represents the M×M identity matrix; represents the noise vector n i The autocovariance matrix of [l]; Calculate the frequency f i Up (i,i) [p] is the expectation between frames and is converted into a zero-mean variable in is the corresponding source signal vector, and E{·} represents taking the mean of the input variables; Similarly, the frequency f j Where for Calculate z (i,i) [p] and z (j,j) The cross-covariance matrix and auto-covariance matrix of [p] are uniformly expressed as R (i,j) , then R (i,j) Vectorize and get a new virtual array model: y (i,j) =vec(R (i,j) )=B (i,j) u (i,j) (1) Among them, y (i,j) is the array receiving signal vector; B (i,j) =[b(θ1),b(θ2),…,b(θ K )] is the equivalent steering matrix, where u (i,j) is the equivalent source signal vector of a single snapshot; When i≠j, R (i,j) is the cross-covariance matrix, and the corresponding virtual array represents the cross-double difference matrix; when i=j, R (i,j) is the autocovariance matrix, and the corresponding virtual array represents the self-double difference matrix.

3. The method according to claim 2, characterized in that The step 3 comprises: Let K g Indicates the number of angle search grids, represents the angle search grid in sparse recovery; in the compressed sensing framework, formula (1) is converted into a sparse representation: in As an overcomplete representation of the equivalent steering matrix, the sparse vector Each element in represents the direction of arrival result on the corresponding angle search grid; Based on group sparsity, multiple frequency pairs (f i ,f j ) is expressed as: in, Represents the corresponding predefined K g The DOA estimation result of the angle search grid, ∈ is the defined error limit; Solve the optimization problem and obtain K g ×1 vector, The angle in the angle search grid corresponding to the mid-peak value is the target wave direction.

4. A broadband quasi-stationary signal direction of arrival estimation device based on frequency pairs, characterized in that: include: A second-order difference synthesis array construction module is used to construct a second-order difference synthesis array corresponding to different frequencies in a frequency range of interest according to a far-field broadband quasi-stationary signal received by a physical array; A double difference synthesis array construction module is used to combine the frequency of interest with itself, and to combine two frequencies with each other to obtain multiple frequency pairs; Calculate the difference synthesis matrix of the second-order difference synthesis matrix of each frequency pair to obtain a double difference synthesis matrix model corresponding to the frequency pair; the double difference synthesis matrix corresponding to the self-combination frequency is called a self-double difference synthesis matrix, and the double difference synthesis matrix corresponding to the mutual-combination frequency is called a mutual-double difference synthesis matrix; The solution module is used to construct an optimization problem by combining the self-double difference synthesis array and the mutual-double difference synthesis array, and realize the direction of arrival estimation by using group sparse solution.

Images