DOA estimation method based on VSSLMS and FFT algorithm

Abstract

The invention discloses a DOA estimation method based on a VSSLMS and an FFT algorithm, and the method comprises the steps: firstly obtaining a weight value containing target direction information based on array receiving data and a VSSLMS beam forming method; carrying out phase rotation search on the obtained weight; performing small-range phase rotation search to obtain an optimal phase offset; and finally, obtaining a final high-precision direction estimation value based on the optimal phase deviation. Spatial spectrum estimation operation with high calculation complexity is avoided, the complexity lower than that of a spatial spectrum estimation method and the DOA estimation performance higher than that of the spatial spectrum estimation method are obtained through an FFT calculation method, and the method can be used for fast estimation of a phase control array on the target orientation.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of signal processing, and particularly relates to a DOA estimation method based on VSSLMS (fast variable step size least mean square algorithm) and FFT algorithm. BACKGROUND

[0002] The DOA estimation has been widely applied in many fields such as wireless communication, radar, sonar, medical imaging and radio astronomy. In the past few decades, the subspace decomposition-based methods, such as multiple signal classification (MUSIC) and estimation of signal parameters via rotational invariance techniques (ESPRIT), have been widely concerned. However, the subspace decomposition method needs to perform eigenvalue decomposition on the sample covariance matrix, and the calculation amount is large. In order to reduce the calculation amount, the document "A Fast Adaptive Beamforming Technique for Efficient Direction-of-Arrival Estimation" uses the fast variable step size least mean square (VSSLMS) beamforming method in the document "Fast and Robust Variable-Step-Size LMS Algorithm for Adaptive Beamforming", and proposes a fast DOA estimation method. The method adaptively updates the step size according to the estimated signal power and the cumulative square instantaneous error, and obtains the spatial spectrum by solving the reciprocal of the antenna array directivity pattern, and then estimates the DOA of the signal by using the peak position of the spatial spectrum. The DOA estimation performance provided by the algorithm is equivalent to that of the MUSIC algorithm, and the sample covariance matrix information and the eigenvalue decomposition calculation are not needed, so that the calculation efficiency is improved. However, the high-efficiency DOA estimation method still needs to perform spatial spectrum estimation. As with all DOA algorithms based on spectrum estimation, the vector multiplication involved in the spatial spectrum estimation needs a large number of sampling points, and the calculation complexity is O(MS), where M represents the number of array elements, S represents the number of sampling points used in the spectrum estimation process, and S>>M is satisfied, which leads to a huge calculation amount requirement, and limits the application in engineering scenarios. SUMMARY

[0003] The application proposes a DOA estimation method based on VSSLMS and FFT algorithm to avoid high calculation amount and improve the performance of DOA estimation.

[0004] The technical scheme for achieving the object of the application is as follows: a DOA estimation method based on VSSLMS and FFT algorithm, comprising:

[0005] According to the array receiving data at time t and the number of snapshots N, the weight vector at time t+1 is iteratively calculated;

[0006] According to the weight w(N) at the last moment obtained in step 1, a new weight vector is constructed;

[0007] According to the FFT point number L, a numerical interval [-1 / L, 1 / L] is constructed, G equidistant values from -1 / L to 1 / L are generated, and a phase rotation vector is constructed;

[0008] According to the constructed weight vector and phase rotation vector, the gth equal-part FFT spectrum of L points is obtained by zero padding;

[0009] The gth equal-part FFT spectrum of L points is searched for a spectrum peak, and the amplitude A corresponding to the signal q in the gth search is obtained qg And the position l qg ;

[0010] According to the amplitude A obtained in the gth search qg , the optimal phase offset of the signal q is obtained;

[0011] According to the optimal phase offset and the position l qg , the accurate DOA estimation of the signal q is obtained.

[0012] Compared with the prior art, the present application has the following advantages.

[0013] 1. The present application is based on the FFT algorithm method to accelerate the estimation weight of DOA, avoids complex spatial spectrum estimation, and reduces the calculation amount;

[0014] 2. The present application realizes more accurate DOA estimation by using phase search calculation, and also has high calculation efficiency.

[0015] The present application will be described in further detail below with reference to the accompanying drawings. BRIEF DESCRIPTION OF DRAWINGS

[0016] Figure 1 is a schematic diagram of the principle of the present application.

[0017] Figure 2 is a flowchart of the implementation of the present application.

[0018] Figure 3 is a comparison diagram of the number of complex multiplications and the number of array elements of the present application and the spatial spectrum estimation method.

[0019] Figure 4 is a comparison diagram of the spatial spectrum diagram and the spatial angle of the present application and the spatial spectrum estimation method.

[0020] Figure 5 is a comparison diagram of the root mean square error and the signal-to-noise ratio of the present application and the spatial spectrum estimation method. DETAILED DESCRIPTION

[0021] As Figure 1 shown in FIG. 2, a DOA estimation method based on VSSLMS and FFT algorithm, first, based on array receiving data and VSSLMS beam forming method, the weight containing target direction information is obtained; then the obtained weight is phase-rotated to search; then the optimal phase offset is obtained by doing a small range of phase-rotation search; finally, based on the optimal phase offset, the final high-precision direction estimation value is obtained. The implementation steps include the following:

[0022] Step 1: according to the array receiving data x(t) at time t and the number of snapshots N, the weight w(t+1) at time t+1 is iteratively calculated.

[0023] 1a) the signal q is incident to the steering vector a(θ q ) of the antenna array from the direction θ q

[0024]

[0025] wherein represents the array factor of the qth signal at the mth element, q=1,2,…,Q-1, Q is the total number of signals, m=0,1,…,M-1, the variable M represents the number of elements, and the variable symbol The variable d represents the distance between adjacent elements, the variable λ represents the working wavelength, and the symbol T represents transposition;

[0026] 1b) according to the steering vector a(θ q ), the steering matrix A is obtained:

[0027] A=[a(θ0),a(θ1),…,a(θ q ),…,a(θ Q-1 )];

[0028] 1c) according to the steering matrix A, the receiving data x(t) of the antenna array is represented as:

[0029]

[0030] wherein s(t)=[s0(t),s1(t),…,s q (t),…,s Q-1 (t)] T represent Q uncorrelated narrowband signals, s q (t) represents the qth signal, each signal is independent and uncorrelated in time, and obeys complex Gaussian distribution represents the power of the qth signal, and n(t) represents the mean value of 0 and the variance of σ 2 ​Gaussian white noise, which is independent and identically distributed and uncorrelated with the signal, t = 1, 2,..., N represents the sampling time, and N represents the number of snapshots;

[0031] 1d) Calculate the array output signal y(t) and signal power

[0032] y(t) = w H (t)x(t),

[0033]

[0034] where the symbol represents the complex field, and the symbol H represents the conjugate transpose;

[0035] 1e) Subtract the output signal y(t) obtained in 1d) from the reference signal d(t) to obtain the instantaneous error e(t):

[0036] e(t) = d(t) - y(t) = x0(t) - w H (t)x(t);

[0037] 1f) Obtain the mean value of the instantaneous error at time t+1 from the instantaneous error e(t) obtained in 1e): and the cumulative sum δ(t+1):

[0038]

[0039] δ(t+1) = δ(t) + |e(t)|,

[0040] where the symbol || represents the modulus value, e(1) = 0;

[0041] 1g) Combine the output signal y(t) obtained in 1d), the signal power and the mean value of the instantaneous error and the cumulative sum δ(t) obtained in 1f) to obtain the step size coefficient μ(t):

[0042]

[0043] 1h) Iteratively calculate the weight w(t+1) from the step size coefficient μ(t) obtained in 1g), the array receiving data x(t), and the instantaneous error e(t):

[0044] w(t+1) = w(t) + μ(t)x(t)e * (t),

[0045] where here Symbol * denotes the conjugate.

[0046] Step 2: Construct a new weight vector

[0047]

[0048] where the weight vector

[0049] Step 3: According to the FFT point number L, divide the range [-1 / L, 1 / L] into G equal parts, and the g-th equal part value is denoted by symbol g , construct the phase rotation vector g :

[0050]

[0051] where g ∈ [-1 / L, 1 / L], g = 1, 2,..., G.

[0052] Step 4: According to the weight vector constructed in step 2 and the phase rotation vector g constructed in step 3, obtain the FFT spectrum P g of the g-th equal part by zero padding.

[0053] 4a) Perform point multiplication between the weight vector constructed in step 2 and the phase rotation vector g constructed in step 3 to obtain the g-th search weight

[0054]

[0055] where the symbol · denotes point multiplication calculation;

[0056] 4b) Zero-pad the search weight obtained in 4a) to obtain a zero-padded weight w g with length L points:

[0057]

[0058] where

[0059] 4c) According to the zero-padded weight w g obtained in 4b), calculate the FFT spectrum:

[0060] P g = shift(FFT(w g )),

[0061] where the symbol shift() represents a shift operation for exchanging the values of the FFT spectrum in the interval and , the symbol represents a rounding operation, and the symbol FFT() represents an FFT calculation.

[0062] Step 5: According to the FFT spectrum P obtained in step 4, the amplitude A and the position l corresponding to the signal q in the gth search are obtained through a spectrum peak search. g qg qg .

[0063] Step 6: According to the amplitude A obtained in step 5, the optimal phase offset of the signal q is obtained. qg

[0064]

[0065] where the symbol arg represents a condition satisfying an amplitude operation, and the symbol max represents a maximum value operation.

[0066] Step 7: According to the optimal phase offset obtained in step 6 and the position l, the accurate DOA estimation of the signal q is obtained. qg

[0067]

[0068] where the symbol arcsin represents an inverse sine calculation.

[0069] The effects of the present application are further illustrated by the following calculation simulation:

[0070] Simulation 1: The simulation of the relationship between the number of complex multiplications and the number of array elements of the present application and the spatial spectrum estimation method

[0071] 1.1) Simulation conditions: Consider a uniform linear array composed of M array elements with a half-wavelength spacing. In particular, we consider Q=3 narrowband signals, the number of snapshots K=150, the number of sampling points used in the spectrum estimation process S=1024, and the number of FFT points L=1024.

[0072] 1.2) Simulation content and results:

[0073] The computational complexity of the present application and the spatial spectrum estimation method is determined by the number of complex multiplications, the number of complex multiplications of the spatial spectrum estimation method is O(3MK+5K+MS), and the number of complex multiplications of the present application is O(3MK+5K+GMlogL).​​​​​

[0074] The computational complexity of the existing spatial spectrum estimation method and the present invention is compared. Under the simulation conditions of 1.1) above, the relationship between the number of complex multiplications and the number of array elements of the present invention and the spatial spectrum estimation method is simulated. The simulation results are as follows: Figure 3 As shown, the horizontal axis represents the number of array elements varying from 32, 64, 128, 256, and 512, and the vertical axis represents the number of complex multiplications.

[0075] from Figure 3 As can be seen, the number of complex multiplications for the present invention is lower than that for the spatial spectrum estimation method for both G = 5 and G = 10, indicating that the computational complexity of the present invention is lower than that of the spatial spectrum estimation method. Furthermore, the computational complexity of the present invention for G = 10 is only slightly higher than that for G = 5, indicating that a small change in G has little impact on the computational complexity.

[0076] Simulation 2: Simulate the relationship between the spatial spectrum and the spatial angle according to the present invention and the spatial spectrum estimation method.

[0077] 2.1) Simulation Conditions: Consider a uniform linear array consisting of M = 20 elements, with a half-wavelength spacing. Assume that three signals of equal power are incident on the array from angles θ1 = -40°, θ2 = -10°, and θ3 = 40°, respectively. The SNR is set to 10 dB, the number of snapshots is 150, and the uniform sampling points and FFT length for spatial spectrum estimation are both 1024.

[0078] 2.2) Simulation content and results:

[0079] The relationship between the spatial spectrum and the spatial angle of the existing spatial spectrum estimation method and the present invention is compared. Under the simulation conditions of 2.1) above, the relationship between the spatial spectrum and the spatial angle of the present invention and the existing spatial spectrum estimation method is simulated. The simulation results are as follows: Figure 4 As shown, the present invention has an optimal phase shift at signal θ2 = -10°, the horizontal axis represents the spatial angle changing from -70° to 70°, and the vertical axis represents the spatial spectrum.

[0080] from Figure 4 As can be seen, the present invention has a significant peak at angle θ2 = -10° because it has the optimal phase offset at θ2 = -10°. Furthermore, the present invention's estimate at angle θ2 = -10° is closer to the true target than the method based on spatial spectrum estimation, indicating that the present invention has higher DOA estimation accuracy.

[0081] Simulation 3: Simulate the relationship between the root mean square error and the signal-to-noise ratio of the present invention and the spatial spectrum estimation method.

[0082] 3.1) Simulation Conditions: Consider a uniform linear array consisting of M = 20 elements, with a half-wavelength spacing. Assume that three signals of equal power are incident on the array from angles θ1 = -40°, θ2 = -10°, and θ3 = 40°, respectively. The SNR is set to 10 dB, the number of snapshots is 150, the number of uniform sampling points for spatial spectrum estimation, and the FFT length are both 1024, and 200 independent replicates are performed.

[0083] 3.2) Simulation content and results:

[0084] The relationship between the root mean square error RMSE and SNR of the existing spatial spectrum estimation method and the present invention is compared. Under the simulation conditions of 3.1) above, the relationship between the root mean square error RMSE and the signal-to-noise ratio of the present invention and the existing spatial spectrum estimation method is simulated. The simulation results are as follows: Figure 5 As shown, the horizontal axis represents the change of signal-to-noise ratio from -15dB to 25dB, and the vertical axis represents the root mean square error RMSE.

[0085] from Figure 5 It can be seen that the DOA estimation performance improves with the increase of SNR. The two methods with different values ​​of G in the present invention are superior to the method based on spatial spectrum estimation. In addition, when the SNR is greater than 5dB, a larger value of G will have better DOA estimation performance.

Claims

1. A DOA estimation method based on VSSLMS and FFT algorithm, characterized in that: include: Iteratively calculate the weight vector at time t+1 based on the array received data and the number of snapshots N at time t; Construct a new weight vector based on the weight w(N) at the last moment obtained in step 1; According to the number of FFT points L, construct the numerical interval [-1 / L, 1 / L], generate G equally spaced numerical values ​​from -1 / L to 1 / L, and construct the phase rotation vector; According to the constructed weight vector and phase rotation vector, the FFT spectrum of the g-th equal part of point L is obtained by zero padding; Perform a peak search on the FFT spectrum of the g-th equal division at point L to obtain the amplitude A corresponding to the signal q under the g-th search. qg and location l qg ; According to the amplitude A corresponding to the signal q obtained in the g-th search qg , get the optimal phase shift of signal q; According to the obtained optimal phase shift and location l qg , and obtain the accurate DOA estimate of signal q.

2. The DOA estimation method based on VSSLMS and FFT algorithm according to claim 1, characterized in that: Based on the array received data and the number of snapshots N at time t, the specific formula for iteratively calculating the weight vector at time t+1 is: w(t+1)=w(t)+μ(t)x(t)e * (t) Where μ(t) represents the step coefficient, e(t) represents the instantaneous error between the reference signal and the array output signal at time t, w(t) is the weight vector at time t, w(t+1) is the weight vector at time t+1, and the symbol () * represents the conjugate, and x(t) is the received data of the antenna array at time t.

3. The DOA estimation method based on VSSLMS and FFT algorithm according to claim 2, characterized in that: The instantaneous error between the reference signal d(t) and the array output signal y(t) at time t is specifically: e(t)=d(t)-y(t)=x0(t)-w H (t) x (t) In the formula, the symbol () H represents the conjugate transpose, and x0(t) represents the received data of the 0th reference array element.

4. The DOA estimation method based on VSSLMS and FFT algorithm according to claim 2, characterized in that: The step size coefficient is specifically: Where, is the mean of the instantaneous error at time t+1, is the signal power at time t, and δ(t+1) is the instantaneous cumulative error at time t+1.

5. The DOA estimation method based on VSSLMS and FFT algorithm according to claim 1, characterized in that: The constructed new weight vector for: Where w(N) represents the weight obtained at the last moment t=N.

6. The DOA estimation method based on VSSLMS and FFT algorithm according to claim 1, characterized in that: The constructed phase rotation vector is specifically: where η g ∈[-1 / L,1 / L], L is the number of FFT points, g=1,2,…,G, M is the number of array elements.

7. The DOA estimation method based on VSSLMS and FFT algorithm according to claim 1, characterized in that: According to the constructed weight vector and the phase rotation vector Φ(η g ), and the FFT spectrum P of the g-th equal part of point L is obtained by filling zeros g The specific process is: 4a) The constructed weight vector and the phase rotation vector Φ(η g ) point multiplication to obtain the g-th search weight The symbol · represents the dot product calculation; 4b) The obtained search weight Pad zero to get the zero-padding weight w of length L g : 4c) The zero-filling weight w obtained according to 4b) g , calculate the FFT spectrum: P g =shift(FFT(w g )), The symbol shift() represents the shift operation, which is used to exchange the FFT spectrum in the interval and The value between the symbol Indicates a rounding operation, and the symbol FFT() indicates FFT calculation.

8. The DOA estimation method based on VSSLMS and FFT algorithm according to claim 1, characterized in that: The optimal phase shift of signal q Specifically: Among them, arg represents the amplitude operation that meets the conditions, and max represents the maximum value operation.

9. The DOA estimation method based on VSSLMS and FFT algorithm according to claim 1, characterized in that: Accurate DOA estimation of signal q Specifically: Wherein, arcsin represents inverse sine calculation, λ represents the operating wavelength, M represents the number of array elements, and d represents the distance between two adjacent array elements.

Images