CN104021293A - DOA and frequency combined estimation method based on structure least square method - Google Patents

Abstract

The invention provides a DOA and frequency combined estimation method based on a structure least square method. The problem that DOA and frequency combined estimation with a traditional ESPRIT algorithm is disabled for coherent signals is solved. Different from a sampling covariance matrix used in the traditional ESPRIT algorithm, a front-and-back average covariance matrix is adopted in the DOA and frequency combined estimation method to replace the sampling covariance matrix, risks generated when the coherent signals are processed with the ESPRIT algorithm are avoided, and meanwhile the problem that the optimal time domain factor is difficult to obtain is solved. Errors of signal subspace are considered, then a rotation invariant equation is solved through an SLS, the signal subspace higher in accuracy is obtained, and robust estimation of the DOA and the frequency of signals is completed.

Description

A kind of associating angle of arrival-frequency estimating methods based on structure least square method

Technical field

The present invention relates to array signal process technique field, relate in particular to a kind of associating angle of arrival-frequency estimating methods.

Background technology

Array Signal Processing is an important branch in signal process field, and the development of passing through decades has reached its maturity and all had a wide range of applications in a plurality of military affairs such as radar, biologic medical, exploration and astronomy and national economy field.Its principle of work is that a plurality of sensors are formed to sensor array, and utilizes this array that spacing wave is received and processed, and object is to suppress to disturb and noise, extracts the useful information of signal.Different from general signal processing mode, Array Signal Processing is to receive signal by being arranged in the sensor group in space, and utilizes the spatial domain characteristic of signal to come filtering and information extraction.Therefore, Array Signal Processing is also often called as spatial domain signal processing.In addition, Array Signal Processing has that wave beam is flexibly controlled, very strong antijamming capability and the high advantages such as space hyperresolution, thereby has received numerous scholars' concern, and its range of application is constantly increase also.

In recent years, direction of arrival, frequency are combined estimation has important application background and has caused discussion widely in fields such as radar, mobile communication.Between decades, direction of arrival is combined estimation research and development with frequency parameter is rapid in the past, and both accurate estimations can, in the situation that promoting link performance, better guarantee channel information.

Since minimum variance (MV) method of frequency-wave number is come out, orientation-frequency is combined linear prediction method, multidimensional MUSIC method, maximum likelihood method, the ESPRIT method scheduling algorithm of estimation and is come out one after another.In these methods, ML (maximal possibility estimation) is as theoretical optimum, estimates the white Gaussian noise in the situation that it is of equal value with least square method.Although maximal possibility estimation has outstanding statistic property, it need to calculate huge multi-dimensional optimization.Although the contradiction that the algorithm of class ESPRIT has been compromised and existed between estimated accuracy and computation complexity, they are when processing coherent signal, and performance but decays seriously.For overcoming aforesaid shortcoming, when data from the sample survey is done-empty pre-service, and then adopt ESPRIT algorithm to make joint parameter estimation to DOA and frequency, however optimum time domain factor m in practical application ₀be difficult to obtain, because m ₀the linear function about number of sampling points N, i.e. m ₀=(3N+2)=5, along with the increase of number of sampling points, its complexity will increase to O (M ³).

In the situation that signal frequency approaches or is relevant, how DOA and frequency are combined while estimating, improve to greatest extent estimated accuracy, the present invention need to set forth and solve just.

Summary of the invention

The object of the invention is to solve traditional E SPRIT algorithm and combine while estimating DOA and frequency, the problem that coherent signal was lost efficacy.Be different from the sampling covariance matrix using in traditional E SPRIT algorithm, before and after adopting, the present invention replaced to average covariance matrix, evade the risk that ESPRIT algorithm has in processing coherent signal, also overcome the problem that the optimum time domain factor is difficult to obtain simultaneously.Consider the error of signal subspace, then use SLS to solve invariable rotary equation, obtain the signal subspace that precision is higher, complete the Robust Estimation to signal DOA and frequency.

The present invention is achieved through the following technical solutions:

The associating angle of arrival-frequency estimating methods based on structure least square method, comprises the following steps:

1) the data from the sample survey matrix of M * 1 the observation vector x (t) of P the narrow band signal of even linear array that obtains M omnidirectional array element after down-converting to baseband signal

Wherein, P<M, F is sampling rate, m is a sampling sampling subsequence number, the sampling number that N is each sequence;

2) calculate sampling covariance matrix represent X _msampling covariance matrix, before and after then calculating to mean matrix R:

$R=\frac{1}{2}(\hat{R}+\mathrm{\&Pi;}\hat{R}*\mathrm{\&Pi;}),$

Wherein, Π represents switching matrix, and on its back-diagonal, element is 1, and all the other are 0;

3) calculate signal subspace U _s, U _sby p the proper vector formation corresponding to P maximum eigenwert;

4) use under least square LS Algorithm for Solving and establish an equation, obtain Φ _θ, Φ _finitial estimate this initial estimate and U _sall for step 6) initialization of iteration:

$\begin{array}{cc}{J}_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_{\mathrm{\&theta;}}=J_{\mathrm{\&theta;}}^{\&DownArrow;}{U}_s& {\mathrm{\&Psi;}}_{\mathrm{\&theta;}}=\mathrm{T\&Theta;}{T}^{-1}\end{array}$

$J_f^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_f=J_f^{\&DownArrow;}{U}_s,{\mathrm{\&Psi;}}_f=\mathrm{T\&Theta;}{T}^{-1},$

Wherein, $\mathrm{\&Theta;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{1}d/\mathrm{\&lambda;}},...,e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{P}d/\mathrm{\&lambda;}}\},$ $\mathrm{\&Phi;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}{f}_1/F},...,e^{j2\mathrm{\&pi;}{f}_P/F}\},$ λ is signal wavelength, array element distance $d=\mathrm{\&lambda;}/2,J_{\mathrm{\&theta;}}^{\&UpArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{I}_{M-1}& 0_1\end{array}\right],J_{\mathrm{\&theta;}}^{\&DownArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{0}_1& I_{M-1}\end{array}\right],$ $J_f^{\&UpArrow;}=\left[\begin{array}{cc}{I}_{m-1}& 0_1\end{array}\right]\&CircleTimes;I_m,$ $J_f^{\&UpArrow;}=\left[\begin{array}{cc}{0}_1& I_{m-1}\end{array}\right]\&CircleTimes;I_M,$ represent that Kronecker is long-pending, 0 ₁the null vector of (M-1) * 1, I _mbe the unit matrix of M * M, T is the nonsingular matrix of a P * P;

5) be defined as follows two error matrixes:

$E_0=J_0^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;}}-J_{\mathrm{\&theta;}}^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s$

$E_f=J_f^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_f-J_f^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s;$

6) for the K time iteration, K>1, solves following formula with gauss-newton method, obtains Δ Ψ _{θ, k}, Δ Ψ _{f, k}with Δ Ψ _{s, k}:

$\underset{{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k},{\mathrm{\&Delta;\&Psi;}}_{f,k}}{\min }\left|\right|{H_k\&CenterDot;\left[\begin{array}{c}\mathrm{vec}\left\{{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{{\mathrm{\&Delta;\&Psi;}}_{f,k}\right\}\\ \mathrm{vec}\left\{{\mathrm{\&Delta;U}}_{s,k}\right\}\end{array}\right]+\left[\begin{array}{c}\mathrm{vec}\left\{E_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{E_{f,k}\right\}\\ \mathrm{\&kappa;}\&CenterDot;\mathrm{vec}\left\{{\mathrm{\&Delta;U}}_k\right\}\end{array}\right]\left|\right|}^2,$

Wherein, vec{} represents vector calculus, with reformed matrix be respectively ${\stackrel{\&OverBar;}{U}}_{s,k}=U_{s,k-1}+{\mathrm{\&Delta;U}}_{s,k-1},{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;},k}={\mathrm{\&Psi;}}_{\mathrm{\&theta;},k-1}+{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k-1}$ With ${\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{f,k}={\mathrm{\&Psi;}}_{f,k-1}+{\mathrm{\&Delta;\&Psi;}}_{f,k-1},$ E _{θ, k}and E _{f, k}the error matrix that represents the k time iteration, be the evaluated error matrix of the signal subspace of the k time iteration,

$H_k=\left[\begin{array}{ccc}{I}_P\&CircleTimes;\left(J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}\right)& 0& {\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}^T\&CircleTimes;J_{\mathrm{\&theta;}}^{\&UpArrow;}-I_P\&CircleTimes;J_{\mathrm{\&theta;}}^{\&DownArrow;}\\ 0& I_P\&CircleTimes;\left(J_f^{\&UpArrow;}{U}_{s,k}\right)& {\mathrm{\&Psi;}}_{f,k}^T\&CircleTimes;J_f^{\&UpArrow;}-I_P\&CircleTimes;J_f^{\&DownArrow;}\\ 0& 0& \mathrm{\&kappa;}{I}_{\mathrm{mMP}}\end{array}\right],$

$E_{\mathrm{\&theta;},k+1}\&ap;E_{\mathrm{\&theta;},k}+J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k}+J_{\mathrm{\&theta;}}^{\&UpArrow;}{\mathrm{\&Delta;U}}_{s,k}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}-J_{\mathrm{\&theta;}}^{\&DownArrow;}{\mathrm{\&Delta;U}}_{s,k},$

$E_{f,k+1}\&ap;E_{f,k}+J_f^{\&UpArrow;}{U}_{s,k}{\mathrm{\&Delta;\&Psi;}}_{f,k}+J_f^{\&UpArrow;}{\mathrm{\&Delta;U}}_{s,k}{\mathrm{\&Psi;}}_{f,k}-J_f^{\&DownArrow;}{\mathrm{\&Delta;U}}_{s,k};$

7) judge whether to meet: if met, iteration stops, and obtains Ψ _θ, Ψ _ffinal value; If do not met, K=K+1, then gets back to step 5;

8) according to Φ _θ, Φ _ffinal value, obtain the estimation of angle and frequency, be specially: by Ψ _θ, Ψ _ffinal value try to achieve Θ and Φ; Θ and Φ are made to feature decomposition, obtain the estimation of following DOA and frequency:

${\widehat{\mathrm{\&theta;}}}_i={\sin }^{-1}\left(\frac{\mathrm{\&lambda;}\&CenterDot;\&angle;\left({\mathrm{\&alpha;}}_i\right)}{2\mathrm{\&pi;d}}\right)$

${\hat{f}}_i=\frac{F\&CenterDot;\&angle;\left({\mathrm{\&beta;}}_i\right)}{2\mathrm{\&pi;}},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,P$

In above formula, ∠ represents angle computing, α _iand β _irepresent respectively Ψ _θand Ψ _fi eigenwert.

Associating DOA-frequency estimating methods based on structure least square method of the present invention, has solved traditional E SPRIT algorithm and has combined while estimating DOA and frequency, the problem that coherent signal was lost efficacy.Be different from the sampling covariance matrix using in traditional E SPRIT algorithm, before and after adopting, the present invention replaced to average covariance matrix, evade the risk that ESPRIT algorithm has in processing coherent signal, also overcome the problem that the optimum time domain factor is difficult to obtain simultaneously.Consider the error of signal subspace, then use SLS to solve invariable rotary equation, obtain the signal subspace that precision is higher, complete the Robust Estimation to signal DOA and frequency.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of the associating DOA-frequency estimating methods based on structure least square method of the present invention;

Fig. 2 is associating DOA-frequency estimating methods of the present invention and the RMSE of traditional algorithm and the graph of relation of SNR;

Fig. 3 is associating DOA-frequency estimating methods of the present invention and the RMSE of traditional algorithm and the graph of relation of angle intervals;

Fig. 4 is associating DOA-frequency estimating methods of the present invention and the RMSE of traditional algorithm and the graph of relation of frequency interval.

Embodiment

Below in conjunction with accompanying drawing explanation and embodiment, the present invention is further described.

Consider the even linear array of a M omnidirectional array element, the centre frequency f of signal of interest bandwidth _c.(P<M) individual narrow band signal (d that supposes there is P _p(t) }, centre frequency is f _c+ f _p, p=1 ..., P, direction is { θ ₁..., θ _pfar-field signal, act on even linear array.Down-convert to after baseband signal, M * 1 observation vector is

$\begin{array}{c}x\left(t\right)=\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}a\left({\mathrm{\&theta;}}_i\right)d_i\left(t\right)e^{j2\mathrm{\&pi;}{f}_{i}t/F}+n\left(t\right)\\ =A{\mathrm{\&Phi;}}^{t}d\left(t\right)+n\left(t\right)\end{array}---\left(1\right)$

In formula (1), F is sampling rate, be source signal vector, n (t) is that average is zero, variance is additive white Gaussian noise process, I _mthe unit matrix of M * M, array manifold A=[a (θ ₁) ..., a (θ _p)] in p steering vector be

$a\left({\mathrm{\&theta;}}_p\right)={[1,e^{-j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{p}d/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{-j2\mathrm{\&pi;}(M-1)\sin {\mathrm{\&theta;}}_{p}d/\mathrm{\&lambda;}}]}^T---\left(2\right)$

And

$\mathrm{\&Phi;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}{f}_1/F},...,e^{j2\mathrm{\&pi;}{f}_P/F}\}.---\left(3\right)$

Wherein, λ is signal wavelength, array element distance d=λ/2.

The channel of supposing above-mentioned narrow band signal is piece decline, thinks { d within the of short duration sampling interval _p(k) } remain unchanged,

$d\left(t\right)\&ap;d(t+\frac{1}{F})\&ap;...\&ap;d(t+\frac{m-1}{F}).---\left(4\right)$

This means, for same sampling period, the individual sampling of the m of beginning (m < < F) is approximately uniform.Collect m sampling subsequence, wherein each sequence comprises N sampling.Finally, data matrix structure is as follows

By in (1) and (4) substitution formula (5),

X≈A _mD _m+N _m (6)

Wherein

$\begin{array}{c}{A}_m={\left[\begin{array}{cccc}{\left(A\right)}^T& {\left(\mathrm{A\&Phi;}\right)}^T& \&CenterDot;\&CenterDot;\&CenterDot;& {\left({\mathrm{A\&Phi;}}^{m-1}\right)}^T\end{array}\right]}^T\\ D_m=\left[\begin{array}{cccc}d\left(0\right)& \mathrm{\&Phi;d}\left(\frac{1}{F}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&Phi;}}^{N-1}d\left(\frac{N-1}{F}\right)\end{array}\right]\\ N_m=\left[\begin{array}{cccc}n\left(0\right)& n\left(\frac{1}{F}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& n\left(\frac{N-1}{F}\right)\end{array}\right]\end{array}---\left(7\right)$

Make

Like this, formula (6) can be written as

In formula (9), represent that Khatri-Rao is long-pending.

Order represent X _msampling covariance matrix.When source signal is while being relevant, to be difficult to estimate.Therefore, the present invention replaces to mean matrix before and after adopting

$R=\frac{1}{2}(\hat{R}+\mathrm{\&Pi;}\hat{R}*\mathrm{\&Pi;})---\left(10\right)$

In formula, Π represents switching matrix, and on its back-diagonal, element is 1, and all the other are 0.

Notice that P the proper vector corresponding to P maximum eigenwert formed signal subspace U _s, namely, span{U _s}=span{A _m.

Be defined as follows two selection matrixs:

$J_{\mathrm{\&theta;}}^{\&UpArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{I}_{M-1}& 0_1\end{array}\right]---\left(11a\right)$

$J_{\mathrm{\&theta;}}^{\&DownArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{0}_1& I_{M-1}\end{array}\right]---\left(11b\right)$

And

$J_f^{\&UpArrow;}=\left[\begin{array}{cc}{I}_{m-1}& 0_1\end{array}\right]\&CircleTimes;I_M---\left(12a\right)$

$J_f^{\&DownArrow;}=\left[\begin{array}{cc}{0}_1& I_{m-1}\end{array}\right]\&CircleTimes;I_M---\left(12b\right)$

In formula represent that Kronecker is long-pending, 0 ₁it is the null vector of (M-1) * 1.Then, the invariable rotary equation for DOA and frequency can be expressed as

$J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_{\mathrm{\&theta;}}=J_{\mathrm{\&theta;}}^{\&DownArrow;}{U}_s$

$J_f^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_f=J_f^{\&DownArrow;}{U}_s---\left(13\right)$

Wherein

Ψ _θ=TΘT ^-1

Ψ _f＝TΦT ^-1 (14)

T is the nonsingular matrix of a P * P, simultaneously

$\mathrm{\&Theta;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{1}d/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{P}d/\mathrm{\&lambda;}}\}.---\left(15\right)$

By solution formula (14) and Θ and Φ are made to feature decomposition, obtain the estimation of following DOA and frequency

${\widehat{\mathrm{\&theta;}}}_i={\sin }^{-1}\left(\frac{\mathrm{\&lambda;}\&CenterDot;\&angle;\left({\mathrm{\&alpha;}}_i\right)}{2\mathrm{\&pi;d}}\right)---\left(16\right)$

${\hat{f}}_i=\frac{F\&CenterDot;\&angle;\left({\mathrm{\&beta;}}_i\right)}{2\mathrm{\&pi;}},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,P---\left(17\right)$

In formula, ∠ represents angle computing, α _iand β _irepresent respectively Ψ _θand Ψ _fi eigenwert.

Because formula (13) is the equation of a highly structural and overdetermination, along with the wherein increase of overlay elements, the performance that LS (Least square, least square) separates will decline, and this is because hypothesis exists with in be not have error, only need to minimize with in error.Yet actual conditions are that each all exists error, i.e. U in formula (13) _s, Ψ _θand Ψ _fin all there is error.Their accurate estimation can be expressed as with in order to minimize Δ U _s, Δ Ψ _θwith Δ Ψ _f, need to use an iteration minimization process.

Be defined as follows two error matrixes:

$E_{\mathrm{\&theta;}}=J_{\mathrm{\&theta;}}^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;}}-J_{\mathrm{\&theta;}}^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s---\left(18a\right)$

$E_f=J_f^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_f-J_f^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s---\left(18b\right)$

The k time iteration, allow ${\stackrel{\&OverBar;}{U}}_{s,k}=U_{s,k-1}+\mathrm{\&Delta;}{U}_{s,k-1},$ ${\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;},k}={\mathrm{\&Psi;}}_{\mathrm{\&theta;},k-1}+\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k-1}$ With represent respectively with reformed matrix.Allow E _{θ, k}and E _{f, k}the error matrix that represents the k time iteration, therefore in (k+1) inferior iteration, ignores quadratic term Δ U _{s, k}Δ Ψ _{θ, k}with Δ U _{s, k}Δ Ψ _{f, k}, can obtain

$E_{\mathrm{\&theta;},k+1}\&ap;E_{\mathrm{\&theta;},k}+J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}+J_{\mathrm{\&theta;}}^{\&UpArrow;}\mathrm{\&Delta;}{U}_{s,k}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}-J_{\mathrm{\&theta;}}^{\&DownArrow;}\mathrm{\&Delta;}{U}_{s,k}---\left(19a\right)$

$E_{f,k+1}\&ap;E_{f,k}+J_f^{\&UpArrow;}{U}_{s,k}\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}+J_f^{\&UpArrow;}\mathrm{\&Delta;}{U}_{s,k}{\mathrm{\&Psi;}}_{f,k}-J_f^{\&DownArrow;}\mathrm{\&Delta;}{U}_{s,k}---\left(19b\right)$

From above formula, can obtain

$\begin{array}{c}\mathrm{vec}\left\{E_{\mathrm{\&theta;},k+1}\right\}\&ap;\mathrm{vec}\left\{E_{\mathrm{\&theta;},k}\right\}+[I_P\&CircleTimes;\left(J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}\right)]\&times;\mathrm{vec}\left\{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}\right\}\\ +[{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}^T\&CircleTimes;J_{\mathrm{\&theta;}}^{\&UpArrow;}-I_P\&CircleTimes;J_{\mathrm{\&theta;}}^{\&DownArrow;}]\&times;\mathrm{vec}\left\{{\mathrm{\&Delta;\&Psi;}}_{s,k}\right\}\\ \mathrm{vec}\left\{E_{f,k+1}\right\}\&ap;\mathrm{vec}\left\{E_{f,k}\right\}+[I_P\&CircleTimes;\left(J_f^{\&UpArrow;}{U}_{s,k}\right)]\&times;\mathrm{vec}\left\{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}\right\}\\ +[{\mathrm{\&Psi;}}_{f,k}^T\&CircleTimes;J_f^{\&UpArrow;}-I_P\&CircleTimes;J_f^{\&DownArrow;}]\&times;\mathrm{vec}\left\{\mathrm{\&Delta;}{U}_{s,k}\right\}\end{array}---\left(20\right)$

In formula, vec{} represents vector calculus.Meanwhile, definition it is the evaluated error matrix of the signal subspace of the k time iteration.(20) formula is arranged and to obtain following SLS problem for matrix form

$\underset{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}}{\min }\left|\right|{H_k\&CenterDot;\left[\begin{array}{c}\mathrm{vec}\left\{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{{\mathrm{\&Delta;\&Psi;}}_{f,k}\right\}\\ \mathrm{vec}\left\{\mathrm{\&Delta;}{U}_{s,k}\right\}\end{array}\right]+\left[\begin{array}{c}\mathrm{vec}\left\{E_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{E_{f,k}\right\}\\ \mathrm{\&kappa;}\&CenterDot;\mathrm{vec}\left\{\mathrm{\&Delta;}{U}_k\right\}\end{array}\right]\left|\right|}^2---\left(21\right)$

Wherein k>1 is a user-defined parameter, and its effect is to keep Δ U _selement compare E _θand E _fin large.

And

$H_k=\left[\begin{array}{ccc}{I}_p\&CircleTimes;\left(J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}\right)& 0& {\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}^T\&CircleTimes;J_{\mathrm{\&theta;}}^{\&UpArrow;}-I_P\&CircleTimes;J_{\mathrm{\&theta;}}^{\&DownArrow;}\\ 0& I_P\&CircleTimes;\left(J_f^{\&UpArrow;}{U}_{s,k}\right)& {\mathrm{\&Psi;}}_{f,k}^T\&CircleTimes;J_f^{\&UpArrow;}-I_P\&CircleTimes;J_f^{\&DownArrow;}\\ 0& 0& \mathrm{\&kappa;}{I}_{\mathrm{mMP}}\end{array}\right]---\left(22\right)$

Use least square solution and U _srespectively as Ψ _θ, Ψ _fand U _sinitial estimate, when K=1, Δ Ψ _{θ, k}, Δ Ψ _{f, k}with Δ U _{s, k}it can be an initial value arbitrarily.The condition that iteration stops is

$\min \{{\left|\right|{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k}\left|\right|}_F^2,{\left|\right|{\mathrm{\&Delta;\&Psi;}}_{f,k}\left|\right|}_F^2,{\left|\right|{\mathrm{\&Delta;u}}_{s,k}\left|\right|}_F^2\}\&le;\mathrm{\&epsiv;}---\left(23\right)$

In formula, ε >0 is a predefined little constant.When algorithm is at the k time iteration convergence, so with represent respectively Ψ _θand Ψ _ffinal estimated value.The estimated value of DOA and frequency finally can be calculated by formula (16) and (17).

Accompanying drawing 1 is the algorithm flow chart of method provided by the invention.Associating DOA-frequency estimating methods based on structure least square method is specially:

Below, data verification method of the present invention by experiment.Consider the even linear array of an array number M=7, array pitch is d=λ/2.Noise is that average is white Gauss noise zero, that variance is 1, and uncorrelated mutually with signal.Experiment, under RMSE (Root Mean Square Error, root-mean-square error) meaning, is compared ESPIRT algorithm and JAFE algorithm.In JAFE algorithm, L ≈ 0.4M.In follow-up experiment, set L=3.In addition sampling number m=3, fast umber of beats N=64, convergency value ε=10, ^-7, the number of simultaneously supposing signal is known, all simulation results obtain by 1000 Monte Carlo Experiments.

Test the relation of 1 RMSE and signal to noise ratio (S/N ratio)

In this emulation, shown the variation that the performance of RMSE presents with signal to noise ratio (S/N ratio).Steady noise power situation under, from-6dB, to 30dB, carry out variable signal power.The identical signal function of three power is to array, and incident angle is θ ₁=10 °, θ ₂=19 ° and θ ₃=30 °.The centre frequency of signal is respectively f ₁=2MHz, f ₂=2.06MHz and f ₃=2.2MHz.For JAFE algorithm, time domain and spatial domain smoothing factor are all close to 3.From accompanying drawing 2, can observe, in illustrated SNR variation range, the performance of method provided by the invention (proposed in accompanying drawing), is that DOA estimates or Frequency Estimation is all better than ESPRIT and JAFE algorithm.

Test the relation of 2 RMSE and angle intervals

In this experiment, managed to compare the variation tendency that the evaluated error of DOA and frequency presents with angle intervals.In this example, only consider two signals, signal to noise ratio (S/N ratio) is fixed as SNR=0dB, and the centre frequency of two signals is respectively f ₁=2MHz and f ₂=2.15MHz, the DOA of first signal is set as θ ₁=0 °, the DOA of second signal is θ ₂=0 °+Δ θ, wherein Δ θ ∈ [0 °, 16 °].From accompanying drawing 3, can find out, when angle intervals is less than 6 °, method provided by the invention is all better than other two kinds of algorithms in this performance aspect two of DOA and Frequency Estimation.

Test the relation of 3 RMSE and frequency interval

What in this emulation, present is that RMSE is with the variation of frequency interval.The DOA that sets two signals is θ ₁=0 ° and θ ₂=6 °.The frequency of first signal is f ₁=2MHz, the frequency of another one signal is f ₂=(2+ Δ f), Δ f ∈ [0,200] kHz.When frequency interval is very little, as Δ f<100kHz, two signals are that time domain is relevant.Because front and back are to average and SLS algorithm, be easy to find out that the performance of method provided by the invention is better than other two kinds of methods from accompanying drawing 4.When Δ f enough large, the slightly inferior ESPRIT method of JAFE, and together with method provided by the invention fits to ESPRIT.

Method provided by the invention is a kind of associating DOA-frequency estimation algorithm based on SLS-ESPRIT, can not fine processing coherent signal different from traditional ESPRIT algorithm, the processing of the algorithm coherent signal of mentioning also has suitable precision, on estimated performance, is significantly improved.

1. the processing of pair coherent signal.Concrete implementation step is as follows:

First, when sampling frequency is F, extraction number is the signal of m (m < < F), now thinks that the signal in the sampling interval is constant, and when snap is N, the array output data matrix obtaining is secondly, set for X _msampling covariance matrix, because signal is correlated with, with front and back to mean matrix replace the not good sampling covariance matrix of estimated performance reduce the susceptibility of algorithm to coherent signal, the required calculated amount of this step is O (m ²m ²n), add the follow-up calculated amount O (m that the EVD computing of R is produced ³m ³), the calculated amount of this process is O (m ²m ²n+m ³m ³).

2. consider the error existing in subspace, by search iteration, solved signal subspace.Concrete implementation step is as follows:

First, when using ESPRIT algorithm to combine estimation DOA and frequency, also considered the error delta U existing in subspace _s, build the SLS problem as shown in formula (21).Then, utilize SLS to carry out search iteration to constant equation, the calculated amount of every single-step iteration is O (10m ³m ³p ³), with minimum value all meet and be less than after predefine value ε, iteration will stop.Finally, the estimated value of signal subspace is the estimated value of DOA and frequency also respectively by ${\mathrm{\&Psi;}}_{\mathrm{\&theta;}}^{\mathrm{SLS}}={\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}+\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}$ With ${\mathrm{\&Psi;}}_f^{\mathrm{SLS}}={\mathrm{\&Psi;}}_{f,k}+\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}$ Provide, the calculated amount that this process produces is O (10Km ³m ³p ³).

By simulation analysis, method estimated performance provided by the invention, higher than traditional ESPRIT algorithm and JAFE algorithm, has not only completed the estimation of combining to DOA and frequency, also obtains comparatively accurate signal subspace simultaneously.

Above content is in conjunction with concrete preferred implementation further description made for the present invention, can not assert that specific embodiment of the invention is confined to these explanations.For general technical staff of the technical field of the invention, without departing from the inventive concept of the premise, can also make some simple deduction or replace, all should be considered as belonging to protection scope of the present invention.

Claims (2)

1. the associating angle of arrival-frequency estimating methods based on structure least square method, is characterized in that:

Said method comprising the steps of:

1) the data from the sample survey matrix of M * 1 the observation vector x (t) of P the narrow band signal of even linear array that obtains M omnidirectional array element after down-converting to baseband signal

Wherein, P<M, F is sampling rate, m is a sampling sampling subsequence number, the sampling number that N is each sequence;

2) calculate sampling covariance matrix represent X _msampling covariance matrix, before and after then calculating to mean matrix R:

$R=\frac{1}{2}(\hat{R}+\mathrm{\&Pi;}\widehat{R*}\mathrm{\&Pi;}),$

Wherein, Π represents switching matrix, and on its back-diagonal, element is 1, and all the other are 0;

3) calculate signal subspace U _s, U _sby p the proper vector formation corresponding to P maximum eigenwert;

4) use under least square LS Algorithm for Solving and establish an equation, obtain Φ _θ, Φ _finitial estimate this initial estimate and U _sall for step 6) initialization of iteration:

$J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_{\mathrm{\&theta;}}=J_{\mathrm{\&theta;}}^{\&DownArrow;}{U}_s,{\mathrm{\&Psi;}}_{\mathrm{\&theta;}}=\mathrm{T\&Theta;}{T}^{-1}$

$J_f^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_f=J_f^{\&DownArrow;}{U}_s,{\mathrm{\&Psi;}}_f=\mathrm{T\&Theta;}{T}^{-1},$

Wherein, $\mathrm{\&Theta;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{1}d/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{P}d/\mathrm{\&lambda;}}\},$ $\mathrm{\&Phi;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}{f}_1/F},\&CenterDot;\&CenterDot;\&CenterDot;,e^{j2\mathrm{\&pi;}{f}_P/F}\},$ λ is signal wavelength, array element distance $d=\mathrm{\&lambda;}/2,J_{\mathrm{\&theta;}}^{\&UpArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{I}_{M-1}& 0_1\end{array}\right],J_{\mathrm{\&theta;}}^{\&DownArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{0}_1& I_{M-1}\end{array}\right],$ $J_f^{\&UpArrow;}=\left[\begin{array}{cc}{I}_{m-1}& 0_1\end{array}\right]\&CircleTimes;I_M,J_f^{\&DownArrow;}=\left[\begin{array}{cc}{0}_1& I_{m-1}\end{array}\right]\&CircleTimes;I_M,$ represent that Kronecker is long-pending, 0 ₁the null vector of (M-1) * 1, I _mbe the unit matrix of M * M, T is the nonsingular matrix of a P * P;

5) be defined as follows two error matrixes:

$E_{\mathrm{\&theta;}}=J_{\mathrm{\&theta;}}^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;}}-J_{\mathrm{\&theta;}}^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s$

$E_f=J_f^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_f-J_f^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s;$

6) for the K time iteration, K>1, solves following formula with gauss-newton method, obtains Δ Ψ _{θ, k}, Δ Ψ _{f, k}

With Δ U _{s, k}:

$\underset{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k},\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}}{\min }{\left|\right|H_k\&CenterDot;\left[\begin{array}{c}\mathrm{vec}\left\{\mathrm{\&Delta;}{\mathrm{\&psi;}}_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}\right\}\\ \mathrm{vec}\left\{\mathrm{\&Delta;}{U}_{s,k}\right\}\end{array}\right]+\left[\begin{array}{c}\mathrm{vec}\left\{E_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{E_{f,k}\right\}\\ k\&CenterDot;\mathrm{vec}\left\{{\mathrm{\&Delta;U}}_k\right\}\end{array}\right]\left|\right|}^2,$

Wherein, vec{} represents vector calculus, with reformed matrix be respectively ${\stackrel{\&OverBar;}{U}}_{s,k}=U_{s,k-1}+{\mathrm{\&Delta;U}}_{s,k-1},{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;},k}={\mathrm{\&Psi;}}_{\mathrm{\&theta;},k-1}+{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k-1}$ With ${\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{f,k}={\mathrm{\&Psi;}}_{f,k-1}+{\mathrm{\&Delta;\&Psi;}}_{f,k-1},E_{\mathrm{\&theta;},k}$ And E _{f, k}the error matrix that represents the k time iteration, be the evaluated error matrix of the signal subspace of the k time iteration,

$H_k=\left[\begin{array}{ccc}{I}_P\&CircleTimes;\left(J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}\right)& 0& {\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}^T\&CircleTimes;J_{\mathrm{\&theta;}}^{\&UpArrow;}-I_P\&CircleTimes;J_{\mathrm{\&theta;}}^{\&DownArrow;}\\ 0& I_P\&CircleTimes;\left(J_f^{\&UpArrow;}{U}_{s,k}\right)& {\mathrm{\&Psi;}}_{f,k}^T\&CircleTimes;J_f^{\&UpArrow;}-I_P\&CircleTimes;J_f^{\&DownArrow;}\\ 0& 0& {\mathrm{kI}}_{\mathrm{mMP}}\end{array}\right]$

$E_{\mathrm{\&theta;},k}\&ap;E_{\mathrm{\&theta;},k-1}+J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k-1}{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k-1}+J_{\mathrm{\&theta;}}^{\&UpArrow;}{\mathrm{\&Delta;U}}_{s,k-1}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k-1}-J_{\mathrm{\&theta;}}^{\&DownArrow;}{\mathrm{\&Delta;U}}_{s,k-1}$

$E_{f,k}\&ap;E_{f,k-1}+J_f^{\&UpArrow;}{U}_{s,k-1}{\mathrm{\&Delta;\&Psi;}}_{f,k-1}+J_f^{\&UpArrow;}{\mathrm{\&Delta;U}}_{s,k-1}{\mathrm{\&Psi;}}_{\mathrm{f\&theta;},k-1}-J_f^{\&DownArrow;}{\mathrm{\&Delta;U}}_{s,k-1};$

7) judge whether to meet: $\min \{{\left|\right|{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k}\left|\right|}_F^2,{\left|\right|{\mathrm{\&Delta;\&Psi;}}_{f,k}\left|\right|}_F^2,{\left|\right|{\mathrm{\&Delta;U}}_{s,k}\left|\right|}_F^2\}\&le;\&Element;$ If met, iteration stops, and obtains Ψ _θ, Ψ _ffinal value; If do not met, K=K+1, then gets back to step 5;

8) according to Φ _θ, Φ _ffinal value, obtain the estimation of angle and frequency, be specially: by Ψ _θ, Ψ _ffinal value try to achieve Θ and Φ; Θ and Φ are made to feature decomposition, obtain the estimation of following DOA and frequency:

${\widehat{\mathrm{\&theta;}}}_i={\sin }^{-1}\left(\frac{\mathrm{\&lambda;}\&CenterDot;\&angle;\left({\mathrm{\&alpha;}}_i\right)}{2\mathrm{\&pi;d}}\right)$

${\hat{f}}_i=\frac{F\&CenterDot;\&angle;\left({\mathrm{\&beta;}}_i\right)}{2\mathrm{\&pi;}},i=1,...,P$

In above formula, ∠ represents angle computing, α _iand β _irepresent respectively Ψ _θand Ψ _fi eigenwert.

2. method of estimation according to claim 1, is characterized in that: described Δ Ψ _{θ, k}, Δ Ψ _{f, k}with Δ U _{s, k}iteration initial value be appoint to a value.