METHOD OF DIRECTION ESTIMATION FOR NONCIRCULAR SIGNALS VIA DECOUPLED OPTIMIZATION

Information

  • Patent Application
  • 20240230819
  • Publication Number
    20240230819
  • Date Filed
    December 26, 2022
    2 years ago
  • Date Published
    July 11, 2024
    5 months ago
  • Inventors
    • CHOI; YANG HO
Abstract
The disclosure presents a method that estimates the arrival angles of noncircular signals incident on a sensor array based on the ML principle. The complex multidimensional problem according to the ML is transformed into a series of simple problems through the separation of one from the superimposed signals, which leads to the optimization of the function of the two-variable θ and ϕ associated with, respectively, the direction and the initial phase of the separated signal. The optimum value of ϕ is theoretically solved. Then the arrival angle of the separated signal is efficiently estimated through a simple 1-D search with respect to θ. Such an optimization process, updating signal separations, is iterated until the direction estimates converge. The proposed method can be applied also in the presence of Doppler effect.
Description
BACKGROUND
Description of Related Art

The arrival angles of the signals incident onto a sensor array can be estimated by relying on the maximum likelihood (ML) principle, which can provide better performance than on the multiple signal classification (MUSIC) principle that the signal subspace is orthogonal to the noise subspace. However, the former requires solving a nonlinear multidimensional problem. Suboptimal solutions can be obtained by transforming the multidimensional problem into iterative one-dimensional (1-D) problems through signal separations. The separation of one signal from the superimposed signals is made by removing the contributions of the others (A. J. Weiss et al, J. Yin et al). Then, the parameters associated with the separated signal are optimized. Such an optimization is repeated from signal to signal. The basic idea of separation has also been applied to the estimation of time delays based on the least square in the frequency domain (R. Wu et al, J. Li et al).


If the incident signals are noncircular so that the averages of their squared complex envelopes, as in binary phase shift keying (BPSK) and amplitude shift keying (ASK), are nonzero, the noncircularity can be exploited to improve estimation performance (Y.-H. Choi; P. Charge et al.). Various methods for direction finding have been suggested that use noncircularity, but without the consideration of the Doppler effect. If there is relative motion between receivers and transmitters Doppler shifts take place, which can cause severe performance degradation with no consideration of the effect. Recently, based on signal separation, a direction estimation method has been presented that deals with the Doppler effect when noncircular signals impinge on a moving array (B. Yang et al.). It eigen-decomposes a matrix of size N at each search point to find the optimal one, where Nis the number of snapshots used for the estimation. The eigen-decomposition would be computationally very intensive, particularly when N is large. This existing method is termed decoupled estimation method 2 (DEM2).


SUMMARY

This disclosure presents a method that estimates the arrival angles of noncircular signals incident on a sensor array based on the maximum likelihood (ML) principle. The ML estimation of the directions of noncircular signals is formulated as a complex multidimensional problem, which is transformed into a series of simple problems through signal separations. The separation of one from the superimposed signals impinging on the sensor array is made by removing the contributions of the others, which are deduced from the parameters estimated through the previous steps. The formulation for the separated signal results in the optimization of a function of two variables θ and ϕ associated with its direction and initial phase, respectively. Given θ, the optimum of ϕ is theoretically solved, which enables us to efficiently estimate the arrival angles. Then the optimization with respect to θ is carried out through a simple one-dimensional (1-D) search. Such an optimization process, updating signal separations, is iterated until the direction estimates converge. In other words, the iteration is terminated when the differences between the present and the previous estimates are very small. The proposed method, which is referred to as DEM1, can also be applied in the presence of Doppler effect.





BRIEF DESCRIPTION OF THE DRAWINGS

The above and other aspects, and advantages of certain embodiments of the disclosure will be more apparent from the following description taken in conjunction with the accompanying drawings, in which:



FIG. 1 shows the computational procedure of the proposed method, DEM1.



FIG. 2 is a graph comparing the computational complexities of DEM1 and DEM2 with respect to the number of snapshots.



FIG. 3 is a graph comparing the cost functions gk(i)(θ) and hk(i)(θ), k=1, 2, of DEM1 and DEM2, respectively, at the 2nd iteration.



FIG. 4 compares the root mean square error (RMSE) of DEM1 with those of other methods as functions of signal-to-noise ratio (SNR).



FIG. 5 compares the RMSE of DEM1 with those of other methods as functions of the Doppler constant.



FIG. 6 is a block diagram of direction estimator according to an embodiment of the disclosure.





DETAILED DESCRIPTION

Before specifically describing the disclosure, a method for demonstrating the present specification and drawings will be described.


First of all, the terms used in the present specification and the claims are general terms identified in consideration of the functions of the various embodiments of the disclosure. However, these terms may vary depending on intention, legal or technical interpretation, emergence of new technologies, and the like of those skilled in the related art. Also, there may be some terms arbitrarily identified by an applicant. Unless there is a specific definition of a term, the term may be construed based on the overall contents and technological common sense of those skilled in the related art.


Further, like reference numerals indicate like components that perform substantially the same functions throughout the specification. For convenience of descriptions and understanding, the same reference numerals or symbols are used and described in different exemplary embodiments. In other words, although elements having the same reference numerals are all illustrated in a plurality of drawings, the plurality of drawings do not mean one exemplary embodiment.


In the disclosure, relational terms such as first and second, and the like, may be used to distinguish one entity from another entity, without necessarily implying any actual relationship or order between such entities. In embodiments of the disclosure, relational terms such as first and second, and the like, may be used to distinguish one entity from another entity, without necessarily implying any actual relationship or order between such entities.


The terms used herein are solely intended to explain a specific exemplary embodiment, and not to limit the scope of the disclosure. It is to be understood that the singular forms include plural referents unless the context clearly dictates otherwise. The terms “include”, “comprise”, “is configured to,” etc., of the description are used to indicate that there are features, numbers, steps, operations, elements, parts or combination thereof, and they should not exclude the possibilities of combination or addition of one or more features, numbers, steps, operations, elements, parts or a combination thereof.


The term such as “module,” “unit,” “part”, and so on is used to refer to an element that performs at least one function or operation, and such element may be implemented as hardware or software, or a combination of hardware and software. Further, except for when each of a plurality of “modules”, “units”, “parts”, and the like needs to be realized in an individual hardware, the components may be integrated in at least one module or chip and be realized in at least one processor.


Also, when any part is connected to another part, this includes a direct connection and an indirect connection through another medium. Further, when a certain portion includes a certain element, unless specified to the contrary, this means that another element may be additionally included, rather than precluding another element.


Hereinafter, the disclosure is described in detail. A linear array may consists of M sensors, and K narrowband signals impinge on the array from θ={θi, . . . , θK} where θk is the arrival angle of the kth signal. The array response vector is denoted by a(θ) for a direction θ. The received signal vector can be expressed as










x

(
t
)

=



A

(
θ
)



s

(
t
)


+

n

(
t
)






(
1
)







where A(θ)=[a(θ1), . . . , a(θK)], s(t) is a complex envelope vector of the received signals, and n(t) is the noise vector. Noise is assumed to be a circularly symmetric complex Gaussian random process with zero mean and variance σ2 and to be uncorrelated from element to element so that










E
[


n

(
t
)




n
H

(
t
)


]

=


σ
2


I





(
2
)







where E, H, and I designate expectation, complex conjugate transpose, and an identity matrix respectively.


The incoming signals are fully noncircular and their initial phases are ϕ={ϕ1, . . . , ϕK} where ϕk is the initial phase of the kth signal. Then the complex envelope vector s(t) can be written as










s

(
t
)

=


B

(
ϕ
)



γ

(
t
)






(
3
)







where










B

(
ϕ
)

=

diag
[


e

j


ϕ
1



,


,

e

j


ϕ
K




]





(
4
)










γ

(
t
)

=


[



γ
1

(
t
)

,


,


γ
K

(
t
)


]

T





with T standing for the transpose. If sk(t) is a BPSK signal γk(+) has a value of +A or −A where A is an amplitude. In (3), γ(t) is a real vector. Substituting (3) into (1) yields










x

(
t
)

=



A

(
θ
)



B

(
ϕ
)



γ

(
t
)


+

n

(
t
)






(
6
)







For the sake of simplicity, the Doppler effect is not considered in this step and is described later.


N snapshots x(1), . . . , x(N) are available. The arrival angles can be estimated based on the deterministic ML criterion, in which the cost function can be written as










f
ML

=




n
=
1

N






x

(
n
)

-


A

(
θ
)



B

(
ϕ
)



γ

(
n
)





2






(
7
)







where ∥⋅∥ denotes the Euclidean norm. The minimization of ƒML is a nonlinear multidimensional problem. Applying signal separation, we can transform it into iterative 1-D problems. The parameters related to the kth signal are θk, ϕk, and γn(n), n=1, . . . , N. At the kth step in the ith iteration, the estimates θk(i), ϕk(i), and γk(i)(n) for them are attained. To this end, the other signal components, which are deduced from the parameters estimated through the previous steps, are removed from the snapshot data. The resulting decoupled vector can be expressed as











z
k

(
i
)


(
n
)

=


x

(
n
)

-


A

(

θ
k

(
i
)


)



B

(

ϕ
k

(
i
)


)




γ
k

(
i
)


(
n
)







(
8
)







where











γ
k

(
i
)


(
n
)

=


[



y
1

(
i
)


(
n
)

,


,


y

k
-
1


(
i
)


(
n
)

,


y

k
+
1


(
i
)


(
n
)

,


,


y
K

(
i
)


(
n
)


]

T





(
9
)













θ
k

(
i
)


=

{


θ
1

(
i
)


,


,

θ

k
-
1


(
i
)


,

θ

k
+
1


(

i
-
1

)


,


,


θ
K

(

i
-
1

)




{
}








(
10
)







and ϕk(i) is similarly defined. Then the following cost function is minimized to find the estimates:










f
k

(
i
)


=




n
=
1

N







z
k

(
i
)


(
n
)

-


e

j

ϕ




a

(
θ
)



γ

(
n
)





2






(
11
)







Minimizing ƒk(i) with respect to a real number γ(n) leads to










γ

(
n
)

=


Re

(


e


-
j


ϕ





a
H

(
θ
)




z
k

(
i
)


(
n
)


)





a

(
θ
)



2






(
12
)







Re(⋅) and Im(⋅) represent the real and imaginary parts of a complex number, respectively. When replacing γ(n) in (11) by (12), the minimization of ƒk(i) becomes equivalent to the maximization of











g
k

(
i
)


(

θ
,
ϕ

)

=




n
=
1

N




[

Re

(


e


-
j


ϕ





a
H

(
θ
)




z
k

(
i
)


(
n
)


)

]

2





a

(
θ
)



2







(
13
)







It is straightforward to see that gk(i)(θ,ϕ) is written as











g
k

(
i
)


(

θ
,
ϕ

)

=




d
1



e

j

2

ϕ



+


d
1
*



e


-
j


2

ϕ



+

2


d
0




4





a

(
θ
)



2







(
14
)







where










d
0

=




n
=
1

N





"\[LeftBracketingBar]"




a
T

(
θ
)




z
k


(
i
)

*


(
n
)




"\[RightBracketingBar]"


2






(

15

a

)













d
1

=




n
=
1

N




a
T

(
θ
)





z
k


(
i
)

*


(
n
)

2







(

15

b

)







with * denoting the complex conjugate. For notational convenience, the dependency of d0 and d1 on θ were omitted. Note that the real number d0 is independent of ϕ. Let d1=|d1|e. Given θ, clearly the maximum of gk(i)(θ,ϕ) with respect to ϕ becomes











g
k

(
i
)


(
θ
)

=





"\[LeftBracketingBar]"


d
1



"\[RightBracketingBar]"


+

d
0



2





a

(
θ
)



2







(
16
)







when









ϕ
=

-

ψ
2






(
17
)







The estimate θk(i) is obtained as










θ
k

(
i
)


=

arg

max
θ



g
k

(
i
)


(
θ
)






(
18
)







Once θk(i) is discovered, the estimate ϕk(i) is given by (17) with θ=θk(i) and then γk(i)(n) by (12). If, θk(i), ϕk(i)(n) and γk(i)(n) are obtained, zk+1(i)(n) can be calculated as












z

k
+
1


(
i
)


(
n
)

=



z
k

(
i
)


(
n
)

+


y
k

(
i
)


(
n
)

-


y

k
+
1


(

i
-
1

)


(
n
)



,

n
=
1

,


,
N




(
19
)







where











y
p

(
q
)


(
n
)

=


e

j


ϕ
p

(
q
)






a

(

θ
p

(
q
)


)




γ
p

(
q
)


(
n
)






(
20
)







The next step proceeds to find θk+1(i). The iteration is terminated if













"\[LeftBracketingBar]"



θ
k

(
i
)


-

θ
k

(

i
-
1

)





"\[RightBracketingBar]"



ε

,

k
=
1

,


,
K




(
21
)







where ε is a small constant. The computational procedure of the proposed method, DEM1, is shown in FIG. 1 where θ(0)={θ1(0), . . . , θK(0)}, s0(0)(n)=[s0(0)(n), . . . , sK(0)(n)]T, and d0 and d1 are explicitly expressed as functions of θ. With θ(0) obtained, s0(0)(n) can be computed as s0(0)(n)=A(θ(0))(AH(0))A(θ(0)))−1AH(0))x(n). When i=1, zk(i)(n) is calculated as zk(1)(n)=x(n)−A(θk(1))sk(1)(n) where sk(1)(n)=[s1(1)(n), . . . , sk−1(1)(n), sk+1(0)(n), . . . , sK(0)(n)]T.


In the following, for comparison, we briefly introduce the conventional method, DEM2, by B. Yang et al. The direction estimates in DEM2 are also obtained via the minimization of ƒk(i), which can be represented as










f
k

(
i
)


=





z
k

(
i
)


-

μ


C

(
θ
)


u




2





(
22
)







where zk(i)=[zk(i)T(1), . . . , zk(i)T(N)]T, C(θ) is the N×N block diagonal matrix given as C(θ)=blkdiag[a(θ), . . . , a(θ)], μ is a complex number, and u is a real vector with unit norm. Comparing (11) and (22), we see










μ

u

=


e

j

ϕ



γ





(
23
)







where γ=[γ(1), . . . , γ(N)]T is a real vector. It is obvious that for arbitrary μ, ϕ, u, and γ, the set of the complex vectors that has the form of μu is the same as that of eγ. The complex number μ that minimizes (22) for a given θ is









μ
=



u
H



h

(
θ
)







C

(
θ
)


u



2






(
24
)







where










h

(
θ
)

=



C
H

(
θ
)



z
k

(
i
)







(
25
)







The minimization of (22) after the substitution of (24) is reduced to the maximization of











h
k

(
i
)


(

θ
,
u

)

=




"\[LeftBracketingBar]"



u
H



h

(
θ
)




"\[RightBracketingBar]"


2





(
26
)







subject to ∥u∥=1, where ∥a(θ)∥ is assumed to be a constant. The complex vector h(θ) is written as










h

(
θ
)

=



h
r

(
θ
)

+

j



h
i

(
θ
)







(
27
)







where hr(θ) and hi(θ) are real vectors. The maximum of hk(i)(θ, u) with respect to u is equal to the maximal eigenvalue of H′(θ) given as











H


(
θ
)

=


H

(
θ
)




H
T

(
θ
)






(
28
)







where










H

(
θ
)

=

[



h
r

(
θ
)

,



h
i

(
θ
)


]





(
29
)







The estimate of θk(i) is given by










θ
k

(
i
)


=

arg

max
θ



h
k

(
i
)


(
θ
)






(
30
)







where











h
k

(
i
)


(
θ
)

=


λ
max

(


H


(
θ
)

)





(
31
)







The maximization of (30) is solved through the eigen-decomposition of N×N matrices H′(θ) at every search point θ. Once θk(i) is obtained, μk(i) is given by (24) with θ=θk(i) and zk+1(i) can be calculated as










z

k
+
1


(
i
)


=


z
k

(
i
)


+

y
k

(
i
)


-

y

k
+
1


(

i
-
1

)







(
32
)







where yp(q)p(q)C(θp(1))up(q)


The estimates of both (18) and (30) are found by minimizing the cost function (11). As mentioned above, the sets of complex vectors μu and of complex vectors eγ are the same. Hence, (18) and (30) are identical. When γ(n) is given by (12), gk(i)(θ,ϕ) is related to ƒk(i) by











g
k

(
i
)


(

θ
,
ϕ

)

=





z
k

(
i
)




2

-

f
k

(
i
)







(
33
)







In addition, hk(i)(θ, u), when μ is given by (24), is expressed as











h
k

(
i
)


(

θ
,
u

)

=


(





z
k

(
i
)




2

-

f
k

(
i
)



)






a

(
θ
)



2






(
34
)







As a result, under the assumption of the constant











h
k

(
i
)


(
θ
)

=





a

(
θ
)



2




g
k

(
i
)


(
θ
)






(
35
)







It is seen from (35) that when ∥a(θ)∥ is independent of θ the estimates are the same.












TABLE I







DEM1
DEM2










text missing or illegible when filed : O(4NM + 4N)


text missing or illegible when filed : O(8NM + 2N)




dtext missing or illegible when filed dtext missing or illegible when filed : O(4NM + 3N)

text missing or illegible when filed : O(4NM + Ntext missing or illegible when filed  + N)




gtext missing or illegible when filed (θ): O(0)
htext missing or illegible when filed (θ): O(Ntext missing or illegible when filed )







Overall load for the computation of θtext missing or illegible when filed










O(Ntext missing or illegible when filed (4NM +
O(Ntext missing or illegible when filed (4NM + Ntext missing or illegible when filed  +



3N) + 4NM + 4N)
Ntext missing or illegible when filed  + Ntext missing or illegible when filed  + 8NM + 2N)








text missing or illegible when filed indicates data missing or illegible when filed







The computational loads for θk(i) of DEM1 and DEM2 are compared in terms of the number of real multiplications in Table I. Ng designates the number of search points and O(0) indicates that the complexity is neither dependent on N nor M. The decoupled vector zk(i), which does not depend on θ, is calculated, with the replacement of k by k−1, as (19) in DEM1 and as (32) in DEM2. The computation of CH(θ)zk(i) requires multiplication of O(4NM). At each grid point θ, the quantities d0, d1, and gk(i)(θ) in DEM1 or H′(θ) and hk(i)(θ) in DEM2 are evaluated, the complexities of which are O(4NM) and O(N3+N2+4NM), respectively. In FIG. 2 the complexities at each search point are compared when M=5. The computational load of DEM2 is huge for large N as it is essentially proportional to N3. For example, when N=200, its complexity is larger than O(8×106) whereas that of DEM1 is less than O(5×103). As can be seen from FIG. 2 and Table I, the computational cost of the proposed method is much less than that of the existing method.


Relative motion between receivers and transmitters brings about the Doppler effect. If the array is moving along its baseline, the Doppler shift is given by ƒDCv sin θ/c where ƒC and θ are the center frequency and the arrival angle of an incoming signal, respectively, and v and c are the velocities of the array and the electromagnetic wave, respectively. The phase shift by the Doppler effect is φD(θ,t)=φ0t sin θ where φ0=2πƒCv/c. We have tacitly considered that the unit of time is a sampling period Ts. Otherwise φ0=2πƒCTsv/c. The Doppler effect can be included in the array response vector instead of the complex envelope so that a(θ) is replaced by










a

(

θ
,

t

)

=


e

j



ϕ
D

(

θ
,
t

)





a

(
θ
)






(
36
)







The equations presented above can be applied in the presence of the Doppler effect via such replacement. For example, the nth diagonal entry of the block diagonal matrix C(θ) is a(θ, n), rather than a(θ).


We investigate the estimation performance of the proposed method of the disclosure in comparison with those of other existing ones. To this end, a uniform linear array of five sensors with half a wavelength interelement spacing is employed, on which two BPSK signals are incident from θ1=15° and θ2=25° relative to broadside. They have the same signal-to-noise ratio (SNR). A total of 100 independent runs have been performed to find the average root mean square error (RMSE) for the arrival angles. When N=200 the performances of DEM1, MUSIC, and noncircular MUSIC (NC-MUSIC) are presented together with the CRB (Cramér-Rao bound). The MSE of θk is obtained as MSE(θk)=Σl=1Lk,l−θk)2/L where L=100 and θk,l is the estimate for θl at the lth trial. The average RMSE is defined as the square root of Σk=1KMSE(θk)/K. Accordingly the corresponding CRB is obtained by taking the square root of the average of the bounds. As explained above, DEM1 has the same performance as DEM2. The direction estimates of MUSIC are used for the initial values of the DEM. The termination threshold is set at ε=0.01°. In FIGS. 3 and 4, φ0=0.005.


In FIG. 3, where σ2 is set at 0.001, gk(i)(θ) and hk(i)(θ), k=1, 2 at the 2nd iteration in one simulation run are illustrated. As the sensors are isotropic, ∥a(θ)∥2=5 regardless of θ. These results show that hk(2)(θ) is identical to 5gk(2)(θ), which confirms (35) and ensures that DEM1 and DEM2 have the same performance.


In FIG. 4, the RMSE performances are displayed as functions of SNR. One can see from (37) that when all signals have the same power, the CRB is inversely proportional to SNR, which is shown in FIG. 4. DEM and NC-MUSIC that exploit the noncircularity of received signals outperform MUSIC. However, NC-MUSIC suffers from the Doppler effect so that its performance improvement due to an increase in SNR is limited. In contrast, the RMSE of DEM, which is close to the CRB, decreases with an increase in SNR without such limitation.


When SNR is 10 dB FIG. 5 exhibits the RMSEs against the Doppler constant φ0, which varies from 0 to 10−2. Although the RMSE of NC-MUSIC when φ0 is near zero is very close to the CRB, its performance is severely degraded as φ0 increases. However, the performance of DEM, which takes the Doppler effect into consideration, is not deteriorated even though it is increased. The signal subspace corresponds to the column space of A(θ) regardless of whether or not φ0 is zero, and the performance of MUSIC is also not degraded.



FIG. 6 is a block diagram of direction estimator according to an embodiment of the disclosure. The process of FIG. 1 described above may be performed in the direction estimator, and the direction estimator may be implemented in software and/or hardware and may be included in one device or in separate devices.

Claims
  • 1. A method, when K noncircular signals are incident on a linear array from θ={θ1, . . . , θK} where θk is the arrival angle of the kth signal, comprising: sampling the received signals to obtain snapshot data x(n), n=1 . . . , N;obtaining a decoupled vector for one signal which is found through signal separation from the snapshot data;finding a function of two variables for the separated signal by using the decoupled vector, where the two variables are θ and ϕ related to its arrival angle and initial phase, respectively;obtaining a cost function of θ from the two-variable function by theoretically solving the optimum value of ϕ; andestimating the angle of arrival of the separated signal through the optimization of the cost function of θ.
  • 2. The method of claim 1, after obtaining all K direction estimates, comprising: iterating the process of optimization with new decoupled vectors that are found using the previous estimates; andterminating the iteration if the condition that |θk(i)−θk(i-1)|≤ε, k=1, . . . , K, is satisfied where θk(j) is the estimate of θk at the jth iteration and ε is a small constant.
  • 3. The method of claim 2, wherein the received signal vector at time t corresponds to the following equation,
  • 4. The method of claim 3, at the kth step in the ith iteration to estimate θk, wherein the decoupled vector is calculated as
  • 5. The method of claim 4, comprising: finding the two-variable function as
  • 6. The method of claim 5, at each θ, comprising: obtaining the optimum value of ϕ(θ) as
  • 7. The method of claim 6, comprising: obtaining the estimates of θk and ϕk as
  • 8. The method of claim 7, comprising: calculating γk(i)(n) as
  • 9. The method of claim 2, wherein in the presence of Doppler effect by motion of the linear array, a(θ) and A(θ) are replaced by a(θ, t)=ejφD(θ,t)a(θ) and A(θ, t)=[a(θ1, t), . . . , a(θK, t)], where, when the center frequency is ƒC and the velocity of the linear array is v in its baseline direction, φD(θ,t)=φ0t sin θ with φ0=2πƒCv/c and c denoting the speed of light.
  • 10. The method of claim 9, wherein the received signal vector at time t corresponds to the following equation,
  • 11. The method of claim 10, wherein the decoupled vector is calculated as
  • 12. The method of claim 11, comprising: finding the two-variable function, using the equality of ∥a(θ, n)∥=∥a(θ)∥, as
  • 13. The method of claim 12, at each θ, comprising: obtaining the optimum value of ϕ(θ) as
  • 14. The method of claim 13, comprising: obtaining the estimates of θk and ϕk as
  • 15. The method of claim 14, comprising: calculating the γk(i)(n) as