Open AccessArticle

Joint Estimation of DOA and Frequency of Multiple Sources with Orthogonal Coprime Arrays

Kai-Chieh Hsu

^†

and

Jean-Fu Kiang

^*,†

Department of Electrical Engineering, National Taiwan University, Taipei 106, Taiwan

Author to whom correspondence should be addressed.

^†

These authors contributed equally to this work.

Sensors 2019, 19(2), 335; https://doi.org/10.3390/s19020335

Submission received: 22 November 2018 / Revised: 3 January 2019 / Accepted: 12 January 2019 / Published: 15 January 2019

(This article belongs to the Section Physical Sensors)

Download

Browse Figures

Versions Notes

Abstract

A two-stage method is proposed to jointly estimate the direction-of-arrival (DOA) and carrier frequency (CF) of multiple sources, by using two orthogonal coprime arrays (CPAs). The DOAs of CF-known sources are estimated first by applying a spatial smoothing MUSIC algorithm. The contribution of these source signals is then removed from the originally received signal by applying an orthogonal complement projector. Next, a joint-ESPRIT algorithm is applied to estimate the DOAs and CFs of the remaining CF-unknown sources. With two orthogonal CPA(5, 6), the RMSE of DOA and CF of applying the proposed method to 30 sources, 13 of which have unknown CF, is less than 1% at SNR

> 5

dB.

Keywords:

array signal processing; direction-of-arrival; carrier frequency; coprime array

1. Introduction

The direction-of-arrival (DOA) of source signals is a prerequisite in many applications like cognitive radar and antenna beamforming. Typical algorithms for estimating the DOAs include subspace-based MUSIC [1], ESPRIT [2] and sparsity-based compressive sensing (CS) approach [3], among others. However, these algorithms are applicable only when the carrier frequency (CF) of signal sources is known [1,2,3,4]. For subspace-based algorithms, MUSIC [1] can generate multiple steering vectors and calculate their projection norm in the null space of receiving signals. On the other hand, ESPRIT [2] can calculate the phase delays if the sensor array possesses rotational invariance, then the DOA can be estimated under a known carrier frequency. As for compressed sensing-based algorithms [4], minimization of

ℓ_{1}

-norm is imposed on the sensing matrix to determine the sparse coding. However, due to Doppler shift [5] or in applications of radar tracking [6] and cognitive radios [7,8], the carrier frequency of signal sources can be ambiguous.

Recently, several algorithms were proposed to jointly estimate the DOA and CF of source signals [7,8,9,10]. In [9], a joint angle-frequency estimation (JAFE) algorithm was proposed, in which both the CF and DOA were iteratively updated from their previous values. In [7,8], a compressed carrier and DOA estimation (CaSCADE) algorithm was proposed. A modulated wideband converter (MWC) was used to estimate the CF at sub-Nyquist rate. Then, a CS approach, a parallel factor (PARAFAC) approach and a joint-ESPRIT (JE) were applied to jointly estimate the CF and DOA. In [10], a two-stage method on an uniform linear array (ULA) was proposed. The frequency components of the received signal were determined by applying Fourier transform and spectrum peak searching. Candidate DOAs at each frequency component were estimated with subspace-based methods like MUSIC. After all the CFs and DOAs were estimated, an exhaustive matching algorithm was applied to match each CF to one DOA, at high computational cost. In practice, the received signal is contributed by a combination of CF-known sources and CF-unknown sources. The performance of estimation algorithms like JE degrade as the number of signal sources increases [7,8].

A Khatri–Rao (KR) subspace was proposed to extend the degrees of freedom (DOFs) of an N-element ULA from

N - 1

2 N - 2

so that DOAs of more signal sources can be estimated [11]. In [12,13], KR subspace-based algorithms were proposed to reduce the uncertainty in gain and phase of sensors. Different array configurations such as nested array [14], coprime array (CPA) [15,16] were proposed to further increase the DOFs. In [4], two types of generalized CPA structure were reviewed and their DOFs were estimated. In [17], the DOF of a CPA was extended by scavenging the fourth-order statistics of the received signals. In [18], a triply primed array (TPA), based on fourth-order statistics, was proposed to increase the DOF beyond CPAs. The effects of mutual interference between sensors in an array were also discussed in [19,20,21].

In applications of beamforming [15,16] or cognitive radios [22], a large number of source signals are usually involved. Previous algorithms for joint estimation of DOA and CF [7,8,9,10] were applicable to a limited number of sources. In this work, a two-stage algorithm is proposed to estimate the DOAs of signal sources with known CF first, followed by the JE algorithm to estimate the DOAs and CFs of the remaining signal sources. Two orthogonal CPAs are used to increase the DOFs. This paper is organized as follows. The signal model is presented in Section 2, the proposed two-stage joint estimation algorithm is presented in Section 3, the simulation results are discussed in Section 4 and, finally, some conclusions are drawn in Section 5.

2. Signal Model

2.1. Sub-Nyquist Sampling Scheme

In this work, a sub-Nyquist sampling scheme [7,8] is used, where the sampling frequency is on the same order as the bandwidth of source signals, which is much lower than the carrier frequency. Figure 1 shows two orthogonal CPAs deployed along x and z axes, respectively, to receive signals from M sources, with the CF and DOA of the m-th source (

s_{m} (t)

) being

f_{m}

and

θ_{m}

, respectively. Without loss of generality, assume that all the signals propagate in the

x z

plane and are uncorrelated to one another. Also assume the first

M_{c}

sources have the same CF of

f_{0}

. The bandwidth of each source signal is limited to

B = 1 / T

, where T is the period of baseband signals. The other signals with unknown CF are assumed to fall in disjoint bands, namely,

min_{m, n > M_{c}, m \neq n} | f_{m} - f_{n} | > B

Figure 2 shows the configuration of a CPA (

N_{1}

N_{2}

), which is composed of two ULAs at spacings of

N_{1} d

and

N_{2} d

, respectively; where

N_{1}

and

N_{2}

are coprime integers and d is the unit spacing. The first ULA has

N_{2}

sensors at spacing of

N_{1} d

and the second ULA has

2 N_{1}

sensors at spacing of

N_{2} d

. Both ULAs share the same first sensor, which is located at the origin, leading to a total of

N = 2 N_{1} + N_{2} - 1

physical sensors in the CPA (

N_{1}

N_{2}

) configuration. The locations of the physical sensors in these two ULAs and the whole CPA are specified by

Γ_{1}

Γ_{2}

and

Γ_{p}

, respectively, with

\begin{matrix} Γ_{1} = \{n_{2} N_{1} d | 0 \leq n_{2} \leq N_{2} - 1\} \\ Γ_{2} = \{n_{1} N_{2} d | 0 \leq n_{1} \leq 2 N_{1} - 1\} \\ Γ_{p} = Γ_{1} \cup Γ_{2} = {χ_{1} d, χ_{2} d, \dots, χ_{N} d} \end{matrix}

where

χ_{1} < χ_{2} < \dots < χ_{N}

The received signal by the n-th sensor along the x axis can be represented as

\begin{matrix} u_{n} (t) = \sum_{m = 1}^{M} s_{m} (t + χ_{n} τ_{x} (θ_{m})) e^{j 2 π f_{m} (t + χ_{n} τ_{x} (θ_{m}))} ≃ \sum_{m = 1}^{M} s_{m} (t) e^{j 2 π f_{m} (t + χ_{n} τ_{x} (θ_{m}))} \end{matrix}

where

τ_{x} (θ) = d cos θ / c

, and

χ_{n} τ_{x} (θ)

is the time advance of the n-th sensor with respect to the first one at the origin. The Fourier transform of

u_{n} (t)

\begin{matrix} U_{n} (f) = \int_{- \infty}^{\infty} d t u_{n} (t) e^{- j 2 π f t} = \sum_{m = 1}^{M} e^{j 2 π f_{m} χ_{n} τ_{x} (θ_{m})} \int_{- \infty}^{\infty} d t s_{m} (t) e^{- j 2 π (f - f_{m}) t} \\ = \sum_{m = 1}^{M} S_{m} (f - f_{m}) e^{j 2 π f_{m} χ_{n} τ_{x} (θ_{m})} \end{matrix}

where

S_{m} (f)

is the spectrum of the m-th baseband signal, with

| f | \leq B / 2

Consider a periodic signal

p (t) = \sum_{v = - \infty}^{\infty} c_{v} e^{j 2 π f_{p} v t}

with period

T_{p} = 1 / f_{p}

, where

c_{v} = \frac{1}{T_{p}} \int_{0}^{T_{p}} p (t) e^{- j 2 π f_{p} v t} d t

. The received signal

u_{n} (t)

is multiplied with

p (t)

to derive

\begin{matrix} {\tilde{x}}_{n} (t) = p (t) u_{n} (t) = \sum_{m = 1}^{M} e^{j 2 π f_{m} χ_{n} τ_{x} (θ_{m})} \sum_{v = - \infty}^{\infty} c_{v} s_{m} (t) e^{j 2 π (f_{m} + v f_{p}) t} \end{matrix}

which is Fourier transformed to have

\begin{matrix} {\tilde{X}}_{n} (f) = \int_{- \infty}^{\infty} d t {\tilde{x}}_{n} (t) e^{- j 2 π f t} = \sum_{m = 1}^{M} e^{j 2 π f_{m} χ_{n} τ_{x} (θ_{m})} \sum_{v = - \infty}^{\infty} c_{v} S_{m} (f - f_{m} - v f_{p}) \end{matrix}

By applying an ideal low-pass filter (LPF) with cutoff frequency

f_{c} / 2

{\tilde{X}}_{n} (f)

, we obtain

\begin{matrix} X_{n} (f) = \{\begin{matrix} {\tilde{X}}_{n} (f), & | f | \leq f_{c} / 2 \\ 0, & otherwise \end{matrix} \end{matrix}

By choosing

f_{c} > B

, the baseband information in

S_{m} (f)

’s is preserved in

X_{n} (f)

. Explicitly,

X_{n} (f)

can be represented as

\begin{matrix} X_{n} (f) = \sum_{m = 1}^{M} {\tilde{S}}_{m} (f) e^{j 2 π f_{m} χ_{n} τ_{x} (θ_{m})} \end{matrix}

with

{\tilde{S}}_{m} (f) = \sum_{v = - V_{1}}^{V_{1}} c_{v} S_{m} (f - f_{m} - v f_{p})

, where

V_{1} = ⌈ (f_{c} + f_{m}) / (2 f_{p}) ⌉ - 1

, which is slightly different from that in [8].

Next, take the inverse Fourier transform of

X_{n} (f)

to have

\begin{matrix} x_{n} (t) = \sum_{m = 1}^{M} e^{j 2 π f_{m} χ_{n} τ_{x} (θ_{m})} {\tilde{s}}_{m} (t) \end{matrix}

which is sampled at the rate of

f_{s} = 1 / T_{s}

to obtain

\begin{matrix} x_{n} [ℓ] = \sum_{m = 1}^{M} e^{j 2 π f_{m} χ_{n} τ_{x} (θ_{m})} w_{m} [ℓ] \end{matrix}

(1)

where

x_{n} [ℓ] = x_{n} (ℓ T_{s})

and

w_{m} [ℓ] = {\tilde{s}}_{m} (ℓ T_{s})

. Note that we may choose

f_{s} = f_{c}

Equation (1) can be put in a matrix form as

\begin{matrix} \bar{x} [ℓ] = {\bar{\bar{A}}}_{x} (\bar{f}, \bar{θ}) \cdot \bar{w} [ℓ] \end{matrix}

(2)

where

\bar{x} [ℓ] = {[x_{1} [ℓ], x_{2} [ℓ], \dots, x_{N} [ℓ]]}^{t}

\bar{w} [ℓ] = {[w_{1} [ℓ], w_{2} [ℓ], \dots, w_{M} [ℓ]]}^{t}

\bar{f} = {[f_{1}, f_{2}, \dots, f_{M}]}^{t}

\bar{θ} = {[θ_{1}, θ_{2}, \dots, θ_{M}]}^{t}

and

\begin{matrix} {\bar{\bar{A}}}_{x} (\bar{f}, \bar{θ}) = [\begin{matrix} e^{j 2 π f_{1} χ_{1} τ_{x} (θ_{1})} & \dots & e^{j 2 π f_{M} χ_{1} τ_{x} (θ_{M})} \\ e^{j 2 π f_{1} χ_{2} τ_{x} (θ_{1})} & \dots & e^{j 2 π f_{M} χ_{2} τ_{x} (θ_{M})} \\ ⋮ & ⋮ \\ e^{j 2 π f_{1} χ_{N} τ_{x} (θ_{1})} & \dots & e^{j 2 π f_{M} χ_{N} τ_{x} (θ_{M})} \end{matrix}] \end{matrix}

Similarly, the received signals in the array along the z-axis are processed to have

\begin{matrix} \bar{z} [ℓ] = {\bar{\bar{A}}}_{z} (\bar{f}, \bar{θ}) \cdot \bar{w} [ℓ] \end{matrix}

(3)

where

\bar{z} [ℓ] = {[z_{1} [ℓ], z_{2} [ℓ], \dots, z_{N} [ℓ]]}^{t}

z_{n} [ℓ] = \sum_{m = 1}^{M} e^{j 2 π f_{m} χ_{n} τ_{z} (θ_{m})} w_{m} [ℓ]

and

\begin{matrix} {\bar{\bar{A}}}_{z} (\bar{f}, \bar{θ}) = [\begin{matrix} e^{j 2 π f_{1} χ_{1} τ_{z} (θ_{1})} & \dots & e^{j 2 π f_{M} χ_{1} τ_{z} (θ_{M})} \\ e^{j 2 π f_{1} χ_{2} τ_{z} (θ_{1})} & \dots & e^{j 2 π f_{M} χ_{2} τ_{z} (θ_{M})} \\ ⋮ & ⋮ \\ e^{j 2 π f_{1} χ_{N} τ_{z} (θ_{1})} & \dots & e^{j 2 π f_{M} χ_{N} τ_{z} (θ_{M})} \end{matrix}] \end{matrix}

with

τ_{z} (θ) = d sin θ / c

By taking noise into consideration, (2) and (3) are modified as

\begin{matrix} \bar{x} [ℓ] = {\bar{\bar{A}}}_{x} \cdot \bar{w} [ℓ] + {\bar{n}}_{x} [ℓ], \bar{z} [ℓ] = {\bar{\bar{A}}}_{z} \cdot \bar{w} [ℓ] + {\bar{n}}_{z} [ℓ] \end{matrix}

(4)

where

{\bar{\bar{A}}}_{x}

and

{\bar{\bar{A}}}_{z}

are the abbreviations of

{\bar{\bar{A}}}_{x} (\bar{f}, \bar{θ})

and

{\bar{\bar{A}}}_{z} (\bar{f}, \bar{θ})

, respectively;

{\bar{n}}_{x} [ℓ]

and

{\bar{n}}_{z} [ℓ]

are uncorrelated additive white Gaussian noise vectors at the ℓth sampling time instant. The noise vectors are assumed to have zero mean and covariance matrix

σ_{n}^{' 2} {\bar{\bar{I}}}_{N}

, where

{\bar{\bar{I}}}_{N}

is an

N \times N

identity matrix.

2.2. Second-Order Statistics

Assume the source signals are wide-sense quasi-stationary, which are observed over Q non-overlapped time frames, with each time frame lasting for

L T_{s}

[4,11,18]. The average power of the m-th source signal in the q-th time frame is represented as

\begin{matrix} P_{m}^{q} = E \{{|w_{m} [ℓ]|}^{2}\}, 1 \leq q \leq Q, (q - 1) L \leq ℓ \leq q L - 1 \end{matrix}

Assume

P_{m}^{q}

’s from different sources are uncorrelated, with their first and second moments determined as

\begin{matrix} E_{q} \{P_{m}^{q}\} = μ_{m}, E_{q} \{{(P_{m}^{q} - μ_{m})}^{2}\} = σ_{m}^{2} \end{matrix}

where

E_{q}

means the expectation value is taken over the q-th time frame. The covariance matrices of the received signals in the q-th time frame are determined as

\begin{matrix} {\bar{\bar{R}}}_{x}^{q} = E {{\bar{x}}^{q} [ℓ] {\bar{x}}^{q †} [ℓ]} = {\bar{\bar{A}}}_{x} \cdot {\bar{\bar{D}}}^{q} \cdot {\bar{\bar{A}}}_{x}^{†} + σ_{n}^{' 2} {\bar{\bar{I}}}_{N} \\ {\bar{\bar{R}}}_{z}^{q} = E {{\bar{z}}^{q} [ℓ] {\bar{z}}^{q †} [ℓ]} = {\bar{\bar{A}}}_{z} \cdot {\bar{\bar{D}}}^{q} \cdot {\bar{\bar{A}}}_{z}^{†} + σ_{n}^{' 2} {\bar{\bar{I}}}_{N} \end{matrix}

(5)

where

{\bar{x}}^{q} [ℓ]

is the received signal in the q-th time frame and

{\bar{\bar{D}}}^{q} = diag \{P_{1}^{q}, P_{2}^{q}, \dots, P_{M}^{q}\}

In practice,

{\bar{\bar{R}}}_{x}^{q}

and

{\bar{\bar{R}}}_{z}^{q}

are estimated as

\begin{matrix} {\tilde{\bar{\bar{R}}}}_{x}^{q} = \frac{1}{L} \sum_{ℓ = (q - 1) L}^{q L - 1} {\bar{x}}^{q} [ℓ] {\bar{x}}^{q †} [ℓ], {\tilde{\bar{\bar{R}}}}_{z}^{q} = \frac{1}{L} \sum_{ℓ = (q - 1) L}^{q L - 1} {\bar{z}}^{q} [ℓ] {\bar{z}}^{q †} [ℓ] \end{matrix}

(6)

By concatenating all the columns of

{\bar{\bar{R}}}_{x}^{q}

and

{\bar{\bar{R}}}_{z}^{q}

, respectively, into

N^{2} \times 1

vectors, (5) is transformed to

\begin{matrix} {\bar{y}}_{x}^{'} [q] = vec \{{\bar{\bar{R}}}_{x}^{q}\} = {\bar{\bar{B}}}_{x} \cdot {\bar{P}}^{'} [q] + σ_{n}^{' 2} {\bar{ν}}_{N} \\ {\bar{y}}_{z}^{'} [q] = vec \{{\bar{\bar{R}}}_{z}^{q}\} = {\bar{\bar{B}}}_{z} \cdot {\bar{P}}^{'} [q] + σ_{n}^{' 2} {\bar{ν}}_{N} \end{matrix}

where

{\bar{P}}^{'} [q] = {[P_{1}^{q}, P_{2}^{q}, \dots, P_{M}^{q}]}^{t}

{\bar{ν}}_{N} = vec {{\bar{\bar{I}}}_{N}}

{\bar{\bar{B}}}_{x} = {\bar{\bar{A}}}_{x}^{*} ⊙ {\bar{\bar{A}}}_{x}

and

{\bar{\bar{B}}}_{z} = {\bar{\bar{A}}}_{z}^{*} ⊙ {\bar{\bar{A}}}_{z}

are second-order difference manifold matrices, and ⊙ means column-wise Kronecker product. By subtracting the expectation value over time-frame index q from

{\bar{y}}_{x}^{'} [q]

and

{\bar{y}}_{z}^{'} [q]

, respectively, we have

\begin{matrix} {\bar{y}}_{x} [q] = {\bar{y}}_{x}^{'} [q] - E_{q} {{\bar{y}}_{x}^{'} [q]} = {\bar{\bar{B}}}_{x} \cdot \bar{P} [q] \\ {\bar{y}}_{z} [q] = {\bar{y}}_{z}^{'} [q] - E_{q} {{\bar{y}}_{z}^{'} [q]} = {\bar{\bar{B}}}_{z} \cdot \bar{P} [q] \end{matrix}

where

\bar{P} [q]

{[P_{1}^{q} - μ_{1}, P_{2}^{q} - μ_{2}, \dots, P_{M}^{q} - μ_{M}]}^{t}

Next, by taking the average of entries in

{\bar{y}}_{x} [q]

and

{\bar{y}}_{z} [q]

, respectively, that have the same lags (of same or opposite sign), we derive the observation vectors as

\begin{matrix} {\bar{y}}_{x}^{Φ} [q] = {\bar{\bar{B}}}_{x}^{Φ} \cdot \bar{P} [q], {\bar{y}}_{z}^{Φ} [q] = {\bar{\bar{B}}}_{z}^{Φ} \cdot \bar{P} [q] \end{matrix}

(7)

where

Φ

contains the index of entries having lags in

Γ_{s} = {- N_{s} d, \dots, N_{s} d}

, with

N_{s} = N_{1} N_{2} + N_{1} - 1

. In other words,

{\bar{\bar{B}}}_{x}^{Φ}

and

{\bar{\bar{B}}}_{z}^{Φ}

serve as steering matrices of effective ULAs deployed along x and z axes, respectively, with virtual sensor locations specified by

Γ_{s}

. The set

Φ

associated with the x-direction ULA is the same as that with the z-direction ULA because both ULAs have the same physical configuration.

3. Joint Estimation of DOA and CF

In this Section, a two-stage algorithm is proposed, in which the DOA of CF-known sources are estimated first, then the DOA and CF of the remaining sources are jointly estimated.

3.1. DOA Estimation of Signal Sources with Known CF

The second-order statistics of the received signals can be used to derive Q-manifold signals, which are then processed with subspace-based algorithms like ESPRIT or MUSIC to estimate the DOA and CF of the original sources. If the fourth-order statistics of the source signals are available, the fourth-order covariance matrices can be estimated as

\begin{matrix} {\tilde{\bar{\bar{R}}}}_{x x} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{x}^{Φ} [q] {\bar{y}}_{x}^{Φ †} [q] ≃ E {{\bar{y}}_{x}^{Φ} [q] {\bar{y}}_{x}^{Φ †} [q]} = {\bar{\bar{B}}}_{x}^{Φ} \cdot \bar{\bar{D}} \cdot {\bar{\bar{B}}}_{x}^{Φ †} \\ {\tilde{\bar{\bar{R}}}}_{z z} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{z}^{Φ} [q] {\bar{y}}_{z}^{Φ †} [q] ≃ E {{\bar{y}}_{z}^{Φ} [q] {\bar{y}}_{z}^{Φ †} [q]} = {\bar{\bar{B}}}_{z}^{Φ} \cdot \bar{\bar{D}} \cdot {\bar{\bar{B}}}_{z}^{Φ †} \end{matrix}

where

\bar{\bar{D}} = diag {σ_{1}^{2}, σ_{2}^{2}, \dots, σ_{M}^{2}}

. The covariance matrices

{\tilde{\bar{\bar{R}}}}_{x x}

and

{\tilde{\bar{\bar{R}}}}_{z z}

are then vectorized as

\begin{matrix} {\bar{y}}_{x x} = vec \{{\bar{\bar{R}}}_{x x}\} = {\bar{\bar{C}}}_{x} \cdot \bar{P} \\ {\bar{y}}_{z z} = vec \{{\bar{\bar{R}}}_{z z}\} = {\bar{\bar{C}}}_{z} \cdot \bar{P} \end{matrix}

(8)

where

{\bar{\bar{C}}}_{x} = {\bar{\bar{B}}}_{x}^{*} ⊙ {\bar{\bar{B}}}_{x}

{\bar{\bar{C}}}_{z} = {\bar{\bar{B}}}_{z}^{*} ⊙ {\bar{\bar{B}}}_{z}

, and

\bar{P} = {[σ_{1}^{2}, σ_{2}^{2}, \dots, σ_{M}^{2}]}^{t}

. By taking the average of entries with same lags in

{\bar{y}}_{x x}

and

{\bar{y}}_{z z}

, respectively, we obtain

\begin{matrix} {\bar{y}}_{x x}^{Φ^{'}} = {\bar{\bar{C}}}_{x}^{Φ^{'}} \cdot \bar{P}, {\bar{y}}_{z z}^{Φ^{'}} = {\bar{\bar{C}}}_{z}^{Φ^{'}} \cdot \bar{P} \end{matrix}

(9)

where

Φ^{'}

contains the index of entries having lags in

Γ_{f} = {- N_{f} d, \dots, N_{f} d}

, with

N_{f} = 3 N_{1} N_{2} + N_{1} - N_{2} - 1

[17]. The set

Φ^{'}

associated with the ULA in x-direction is the same as that with the ULA in z-direction because both ULAs have the same physical configuration. Note that subspace-based algorithms are applicable to CPA with contiguous lags in both x and z directions. Next, construct spatial smoothing covariance matrices of the fourth-order as [23]

\begin{matrix} {\bar{\bar{R}}}_{x x}^{'} = [\begin{matrix} {\bar{y}}_{x x}^{Φ^{'}} [1] & {\bar{y}}_{x x}^{Φ^{'}} [2] & \dots & {\bar{y}}_{x x}^{Φ^{'}} [N_{f} + 1] \\ {\bar{y}}_{x x}^{Φ^{'}} [2] & {\bar{y}}_{x x}^{Φ^{'}} [3] & \dots & {\bar{y}}_{x x}^{Φ^{'}} [N_{f} + 2] \\ ⋮ & ⋮ & ⋱ & ⋮ \\ {\bar{y}}_{x x}^{Φ^{'}} [N_{f} + 1] & {\bar{y}}_{x x}^{Φ^{'}} [N_{f} + 2] & \dots & {\bar{y}}_{x x}^{Φ^{'}} [2 N_{f} + 1] \end{matrix}] \end{matrix}

\begin{matrix} {\bar{\bar{R}}}_{z z}^{'} = [\begin{matrix} {\bar{y}}_{z z}^{Φ^{'}} [1] & {\bar{y}}_{z z}^{Φ^{'}} [2] & \dots & {\bar{y}}_{z z}^{Φ^{'}} [N_{f} + 1] \\ {\bar{y}}_{z z}^{Φ^{'}} [2] & {\bar{y}}_{z z}^{Φ^{'}} [3] & \dots & {\bar{y}}_{z z}^{Φ^{'}} [N_{f} + 2] \\ ⋮ & ⋮ & ⋱ & ⋮ \\ {\bar{y}}_{z z}^{Φ^{'}} [N_{f} + 1] & {\bar{y}}_{z z}^{Φ^{'}} [N_{f} + 2] & \dots & {\bar{y}}_{z z}^{Φ^{'}} [2 N_{f} + 1] \end{matrix}] \end{matrix}

(10)

Algorithm 1 lists the major steps of the spatial smoothing MUSIC (SS-MUSIC). In steps 5–7, the conventional MUSIC algorithm is applied to

{\bar{\bar{R}}}_{x x}^{'}

and

{\bar{\bar{R}}}_{z z}^{'}

, respectively, to estimate the DOAs to the nearest grid points, and record them in

{\bar{ϕ}}_{1 x}

and

{\bar{ϕ}}_{1 z}

, respectively. Next, for every entry in

{\bar{ϕ}}_{1 z}

, find the best match in

{\bar{ϕ}}_{1 x}

with DOA difference smaller than

0.5 °

, then record the average of these two DOAs in

{\bar{θ}}_{1}^{'}

Algorithm 1 SS-MUSIC based on two orthogonal CPAs

Input:

Q snapshots of manifold vectors, ${\bar{y}}_{x}^{Φ}$ and ${\bar{y}}_{z}^{Φ}$ in (7)
Total number of signal sources, M
Known carrier frequency, $f_{0}$

Output:

Estimated DOA vector, ${\bar{θ}}_{1}^{'}$

Algorithm:

Estimate fourth-order covariance matrices as
${\tilde{\bar{\bar{R}}}}_{x x} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{x}^{Φ} [q] {\bar{y}}_{x}^{Φ †} [q],$ ${\tilde{\bar{\bar{R}}}}_{z z} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{z}^{Φ} [q] {\bar{y}}_{z}^{Φ †} [q]$
Vectorize, average and remove non-contiguous lags to form ${\bar{y}}_{x x}^{Φ^{'}}$ and ${\bar{y}}_{z z}^{Φ^{'}}$ , as in (8) and (9)
Construct spatial smoothing covariance matrices ${\bar{\bar{R}}}_{x x}^{'}$ and ${\bar{\bar{R}}}_{z z}^{'}$ , as in (10)
Apply singular-value decomposition (SVD) on ${\bar{\bar{R}}}_{x x}^{'}$ and ${\bar{\bar{R}}}_{z z}^{'}$ to have
${\bar{\bar{R}}}_{x x}^{'} = {\bar{\bar{U}}}_{x} \cdot {\bar{\bar{Σ}}}_{x} \cdot {\bar{\bar{V}}}_{x}^{†},$ ${\bar{\bar{R}}}_{z z}^{'} = {\bar{\bar{U}}}_{z} \cdot {\bar{\bar{Σ}}}_{z} \cdot {\bar{\bar{V}}}_{z}^{†}$
Construct null spaces ${\bar{\bar{U}}}_{x n}$ and ${\bar{\bar{U}}}_{z n}$ out of ${\bar{\bar{R}}}_{x x}^{'}$ and ${\bar{\bar{R}}}_{z z}^{'}$ , respectively, from the last $N_{f} - M$ columns of ${\bar{\bar{U}}}_{x}$ and ${\bar{\bar{U}}}_{z}$
Construct matrices ${\bar{\bar{B}}}_{x G} = [{\bar{b}}_{x} (ϕ_{1}), \dots, {\bar{b}}_{x} (ϕ_{N_{g}})]$ and ${\bar{\bar{B}}}_{z G} = [{\bar{b}}_{z} (ϕ_{1}), \dots, {\bar{b}}_{z} (ϕ_{N_{g}})]$ , where $ϕ_{g} = - 70^{°} + (g - 1) Δ ϕ_{g}$ , $Δ ϕ_{g} = 140^{°} / (N_{g} - 1)$ and
${\bar{b}}_{x} (ϕ) = [1, e^{j 2 π f_{0} τ_{x} (ϕ)}, \dots, e^{j 2 π f_{0} N_{f} τ_{x} (ϕ)}],$ ${\bar{b}}_{z} (ϕ) = [1, e^{j 2 π f_{0} τ_{z} (ϕ)}, \dots, e^{j 2 π f_{0} N_{f} τ_{z} (ϕ)}]$
Compute the DOA spectra
$P_{SSM, x} (ϕ_{g}) = {∥ {\bar{\bar{U}}}_{x n}^{†} \cdot {\bar{b}}_{x} (ϕ_{g}) ∥}^{- 2}, 1 \leq g \leq N_{g},$ $P_{SSM, z} (ϕ_{g}) = {∥ {\bar{\bar{U}}}_{z n}^{†} \cdot {\bar{b}}_{z} (ϕ_{g}) ∥}^{- 2}, 1 \leq g \leq N_{g}$
Pick M highest peaks from each spectrum and record their associated angles in ${\bar{ϕ}}_{1 x}$ and ${\bar{ϕ}}_{1 z}$ , respectively.
For every entry in ${\bar{ϕ}}_{1 z}$ , find a best match in ${\bar{ϕ}}_{1 x}$ with DOA difference smaller than $0.5 °$ , then record the average of both DOAs in ${\bar{θ}}_{1}^{'}$ .

3.2. Joint-ESPRIT in Projected Subspace

The DOAs of the signal sources with known CF (

f_{0}

) are estimated and stored in

{\bar{θ}}_{1}^{'}

. The steering vectors associated with these sources are then removed by applying the orthogonal complement projectors

\begin{matrix} {\bar{\bar{Π}}}_{x} = {\bar{\bar{I}}}_{2 N_{s} + 1} - {\bar{\bar{B}}}_{x}^{'} \cdot {({\bar{\bar{B}}}_{x}^{' †} \cdot {\bar{\bar{B}}}_{x}^{'})}^{- 1} \cdot {\bar{\bar{B}}}_{x}^{' †} \\ {\bar{\bar{Π}}}_{z} = {\bar{\bar{I}}}_{2 N_{s} + 1} - {\bar{\bar{B}}}_{z}^{'} \cdot {({\bar{\bar{B}}}_{z}^{' †} \cdot {\bar{\bar{B}}}_{z}^{'})}^{- 1} \cdot {\bar{\bar{B}}}_{z}^{' †} \end{matrix}

(11)

where

N_{s}

is the largest lag in

Γ_{s}

of the second-order virtual arrays,

{\bar{\bar{B}}}_{x}^{'}

and

{\bar{\bar{B}}}_{z}^{'}

are the steering matrices of CPAs in x and z directions, respectively, with contiguous lags specified in

Γ_{s}

Next, apply the joint-ESPRIT (JE) algorithm to estimate the DOA and CF of the remaining sources embedded in the residual signals [8]. Construct two subarrays out of each second-order virtual array in x and z direction, respectively, as

\begin{matrix} {\bar{y}}_{x 1}^{″} [q] = {\bar{\bar{Π}}}_{x} \cdot {\bar{y}}_{x 1}^{Φ} [q] = {\bar{\bar{Π}}}_{x} \cdot {\bar{\bar{B}}}_{x 1}^{Φ} \cdot \bar{P} [q] = {\bar{\bar{B}}}_{x 1}^{″} \cdot \bar{P} [q] \\ {\bar{y}}_{x 2}^{″} [q] = {\bar{\bar{Π}}}_{x} \cdot {\bar{y}}_{x 2}^{Φ} [q] = {\bar{\bar{Π}}}_{x} \cdot {\bar{\bar{B}}}_{x 2}^{Φ} \cdot \bar{P} [q] = {\bar{\bar{B}}}_{x 2}^{″} \cdot \bar{P} [q] \\ {\bar{y}}_{z 1}^{″} [q] = {\bar{\bar{Π}}}_{z} \cdot {\bar{y}}_{z 1}^{Φ} [q] = {\bar{\bar{Π}}}_{z} \cdot {\bar{\bar{B}}}_{z 1}^{Φ} \cdot \bar{P} [q] = {\bar{\bar{B}}}_{z 1}^{″} \cdot \bar{P} [q] \\ {\bar{y}}_{z 2}^{″} [q] = {\bar{\bar{Π}}}_{z} \cdot {\bar{y}}_{z 2}^{Φ} [q] = {\bar{\bar{Π}}}_{z} \cdot {\bar{\bar{B}}}_{z 2}^{Φ} \cdot \bar{P} [q] = {\bar{\bar{B}}}_{z 2}^{″} \cdot \bar{P} [q] \end{matrix}

(12)

where

{\bar{y}}_{x 1}^{Φ} [q]

and

{\bar{y}}_{x 2}^{Φ} [q]

are the observation vectors from the first

2 N_{s}

and the last

2 N_{s}

second-order virtual sensors, respectively, in x direction, at the q-th sampling instant;

{\bar{\bar{B}}}_{x 1}^{Φ}

and

{\bar{\bar{B}}}_{x 2}^{Φ}

are the first

2 N_{s}

rows and the last

2 N_{s}

rows, respectively, of

{\bar{\bar{B}}}_{x}^{Φ}

;

{\bar{\bar{B}}}_{x 1}^{''} = {\bar{\bar{Π}}}_{x} \cdot {\bar{\bar{B}}}_{x 1}^{Φ}

and

{\bar{\bar{B}}}_{x 2}^{''} = {\bar{\bar{Π}}}_{x} \cdot {\bar{\bar{B}}}_{x 2}^{Φ}

. The vectors

{\bar{y}}_{z 1}^{Φ} [q]

{\bar{y}}_{z 2}^{Φ} [q]

and the matrices

{\bar{\bar{B}}}_{z 1}^{Φ}

{\bar{\bar{B}}}_{z 2}^{Φ}

{\bar{\bar{B}}}_{z 1}^{″}

{\bar{\bar{B}}}_{z 2}^{″}

are defined in the same way as their counterparts in x-direction. The matrices

{\bar{\bar{B}}}_{x 1}^{″}

{\bar{\bar{B}}}_{x 2}^{″}

{\bar{\bar{B}}}_{z 1}^{″}

and

{\bar{\bar{B}}}_{z 2}^{″}

satisfy rotational invariance property (RIP), namely,

\begin{matrix} {\bar{\bar{B}}}_{x 2}^{″} = {\bar{\bar{B}}}_{x 1}^{″} \cdot {\bar{\bar{Ψ}}}_{x}, {\bar{\bar{B}}}_{z 2}^{″} = {\bar{\bar{B}}}_{z 1}^{″} \cdot {\bar{\bar{Ψ}}}_{z} \end{matrix}

(13)

with

\begin{matrix} {\bar{\bar{Ψ}}}_{x} = diag \{e^{j 2 π f_{M - M_{c}^{'} + 1} τ_{x} (θ_{M - M_{c}^{'} + 1})}, \dots, e^{j 2 π f_{M} τ_{x} (θ_{M})}\} \\ {\bar{\bar{Ψ}}}_{z} = diag \{e^{j 2 π f_{M - M_{c}^{'} + 1} τ_{z} (θ_{M - M_{c}^{'} + 1})}, \dots, e^{j 2 π f_{M} τ_{z} (θ_{M})}\} \end{matrix}

which is prerequisite to applying the JE algorithm.

Define covariance matrices

\begin{matrix} {\bar{\bar{R}}}_{1} = E \{{\bar{y}}_{x 1}^{″} [q] {\bar{y}}_{z 1}^{″ †} [q]\} = {\bar{\bar{B}}}_{x 1}^{″} \cdot \bar{\bar{D}} \cdot {\bar{\bar{B}}}_{z 1}^{″ †} \\ {\bar{\bar{R}}}_{2} = E \{{\bar{y}}_{x 2}^{″} [q] {\bar{y}}_{z 1}^{″ †} [q]\} = {\bar{\bar{B}}}_{x 1}^{″} \cdot {\bar{\bar{Ψ}}}_{x} \cdot \bar{\bar{D}} \cdot {\bar{\bar{B}}}_{z 1}^{″ †} \\ {\bar{\bar{R}}}_{3} = E \{{\bar{y}}_{x 1}^{″} [q] {\bar{y}}_{z 2}^{″ †} [q]\} = {\bar{\bar{B}}}_{x 1}^{″} \cdot \bar{\bar{D}} \cdot ({\bar{\bar{Ψ}}}_{z}^{†} \cdot {\bar{\bar{B}}}_{z 1}^{″ †}) \\ {\bar{\bar{R}}}_{4} = E \{{\bar{y}}_{x 2}^{″} [q] {\bar{y}}_{z 2}^{″ †} [q]\} = ({\bar{\bar{B}}}_{x 1}^{″} \cdot {\bar{\bar{Ψ}}}_{x}) \cdot \bar{\bar{D}} \cdot ({\bar{\bar{Ψ}}}_{z}^{†} \cdot {\bar{\bar{B}}}_{z 1}^{″ †}), \end{matrix}

(14)

which can be implemented as

\begin{matrix} {\bar{\bar{R}}}_{1} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{x 1}^{″} [q] {\bar{y}}_{z 1}^{″ †} [q], {\bar{\bar{R}}}_{2} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{x 2}^{″} [q] {\bar{y}}_{z 1}^{″ †} [q], \\ {\bar{\bar{R}}}_{3} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{x 1}^{″} [q] {\bar{y}}_{z 2}^{″ †} [q], {\bar{\bar{R}}}_{4} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{x 2}^{″} [q] {\bar{y}}_{z 2}^{″ †} [q] \end{matrix}

(15)

The four matrices in (14) are concatenated to form

\begin{matrix} {\bar{\bar{R}}}_{c} = [\begin{matrix} {\bar{\bar{R}}}_{1} \\ {\bar{\bar{R}}}_{2} \\ {\bar{\bar{R}}}_{3} \\ {\bar{\bar{R}}}_{4} \end{matrix}] = {\bar{\bar{B}}}_{c} \cdot \bar{\bar{D}} \cdot {\bar{\bar{B}}}_{z 1}^{″ †}, {\bar{\bar{B}}}_{c} = [\begin{matrix} {\bar{\bar{B}}}_{x 1}^{″} \\ {\bar{\bar{B}}}_{x 1}^{″} \cdot {\bar{\bar{Ψ}}}_{x} \\ {\bar{\bar{B}}}_{x 1}^{″} \cdot {\bar{\bar{Ψ}}}_{z}^{†} \\ {\bar{\bar{B}}}_{x 1}^{″} \cdot {\bar{\bar{Ψ}}}_{x} \cdot {\bar{\bar{Ψ}}}_{z}^{†} \end{matrix}] \end{matrix}

(16)

Then, apply SVD on

{\bar{\bar{R}}}_{c}

to have

\begin{matrix} {\bar{\bar{R}}}_{c} = [{\bar{\bar{U}}}_{s} {\bar{\bar{U}}}_{n}] [\begin{matrix} {\bar{\bar{Σ}}}_{s} & 0 \\ 0 & 0 \end{matrix}] [\begin{matrix} {\bar{\bar{V}}}_{s}^{†} \\ {\bar{\bar{V}}}_{n}^{†} \end{matrix}], \end{matrix}

where

\begin{matrix} {\bar{\bar{U}}}_{s} = {[\begin{matrix} {\bar{\bar{U}}}_{s 1} \\ {\bar{\bar{U}}}_{s 2} \\ {\bar{\bar{U}}}_{s 3} \\ {\bar{\bar{U}}}_{s 4} \end{matrix}]}_{8 N_{s} \times (M - M_{c}^{'})} = {\bar{\bar{B}}}_{c} \cdot \bar{\bar{T}} \end{matrix}

(17)

spans the same subspace as

{\bar{\bar{B}}}_{c}

does, and

M_{c}^{'}

is the number of signals with known CF, estimated in the first stage.

By substituting (16) into (17), we have

\begin{matrix} {\bar{\bar{U}}}_{s 1} = {\bar{\bar{B}}}_{x 1}^{″} \cdot \bar{\bar{T}} \\ {\bar{\bar{U}}}_{s 2} = {\bar{\bar{B}}}_{x 1}^{″} \cdot {\bar{\bar{Ψ}}}_{x} \cdot \bar{\bar{T}} = {\bar{\bar{U}}}_{s 1} \cdot {\bar{\bar{T}}}^{- 1} \cdot {\bar{\bar{Ψ}}}_{x} \cdot \bar{\bar{T}} \\ {\bar{\bar{U}}}_{s 3} = {\bar{\bar{B}}}_{x 1}^{″} \cdot {\bar{\bar{Ψ}}}_{z}^{†} \cdot \bar{\bar{T}} = {\bar{\bar{U}}}_{s 1} \cdot {\bar{\bar{T}}}^{- 1} \cdot {\bar{\bar{Ψ}}}_{z}^{†} \cdot \bar{\bar{T}} \\ {\bar{\bar{U}}}_{s 4} = {\bar{\bar{B}}}_{x 1}^{″} \cdot {\bar{\bar{Ψ}}}_{x} \cdot {\bar{\bar{Ψ}}}_{z}^{†} \cdot \bar{\bar{T}} = {\bar{\bar{U}}}_{s 1} \cdot {\bar{\bar{T}}}^{- 1} \cdot {\bar{\bar{Ψ}}}_{x} \cdot {\bar{\bar{Ψ}}}_{z}^{†} \cdot \bar{\bar{T}} \end{matrix}

Next, take the inverse of

{\bar{\bar{U}}}_{s 1}

and multiply it to

{\bar{\bar{U}}}_{s 2}

{\bar{\bar{U}}}_{s 3}

and

{\bar{\bar{U}}}_{s 4}

, respectively, to derive

\begin{matrix} {\bar{\bar{U}}}_{12} = {\bar{\bar{U}}}_{s 1}^{- 1} \cdot {\bar{\bar{U}}}_{s 2} = {\bar{\bar{T}}}^{- 1} \cdot {\bar{\bar{Ψ}}}_{x} \cdot \bar{\bar{T}} \\ {\bar{\bar{U}}}_{13} = {\bar{\bar{U}}}_{s 1}^{- 1} \cdot {\bar{\bar{U}}}_{s 3} = {\bar{\bar{T}}}^{- 1} \cdot {\bar{\bar{Ψ}}}_{z}^{†} \cdot \bar{\bar{T}} \\ {\bar{\bar{U}}}_{14} = {\bar{\bar{U}}}_{s 1}^{- 1} \cdot {\bar{\bar{U}}}_{s 4} = {\bar{\bar{T}}}^{- 1} \cdot {\bar{\bar{Ψ}}}_{x} \cdot {\bar{\bar{Ψ}}}_{z}^{†} \cdot \bar{\bar{T}} \end{matrix}

(18)

By applying eigenvalue decomposition (EVD) on the average of these three matrices, we obtain

\begin{matrix} \bar{\bar{Ω}} = \frac{1}{3} ({\bar{\bar{U}}}_{12} + {\bar{\bar{U}}}_{13} + {\bar{\bar{U}}}_{14}) = {\bar{\bar{T}}}^{- 1} \cdot \bar{\bar{Λ}} \cdot \bar{\bar{T}} \end{matrix}

where

\bar{\bar{Λ}}

is a diagonal matrix containing all the eigenvalues of

\bar{\bar{Ω}}

. The first two equations in (19) imply that

\begin{matrix} {\bar{\bar{Ψ}}}_{x} = \bar{\bar{T}} \cdot {\bar{\bar{U}}}_{12} \cdot {\bar{\bar{T}}}^{- 1}, {\bar{\bar{Ψ}}}_{z} = {(\bar{\bar{T}} \cdot {\bar{\bar{U}}}_{13} \cdot {\bar{\bar{T}}}^{- 1})}^{†} \end{matrix}

(19)

Finally, the DOAs and CFs are estimated as

\begin{matrix} {\tilde{θ}}_{m} = {tan}^{- 1} \frac{∠ Ψ_{z, m}}{∠ Ψ_{x, m}}, \\ {\tilde{f}}_{m} = \frac{c}{2 π d} \sqrt{{(∠ Ψ_{x, m})}^{2} + {(∠ Ψ_{z, m})}^{2}}, \\ m = M_{c}^{'} + 1, M_{c}^{'} + 2, \dots, M \end{matrix}

(20)

where

Ψ_{x, m}

and

Ψ_{z, m}

are the (

m - M_{c}^{'}

)-th diagonal entries of

{\bar{\bar{Ψ}}}_{x}

and

{\bar{\bar{Ψ}}}_{z}

, respectively. The proposed projected Joint-ESPRIT (PJE) algorithm is listed in Algorithm 2.

Algorithm 2 Projected Joint-ESPRIT (PJE)

Input:

Q snapshots of manifold vectors, ${\bar{y}}_{x}^{Φ}$ and ${\bar{y}}_{z}^{Φ}$ , in (7)
Total number of signal sources, M
Known carrier frequency, $f_{0}$
Estimated DOA vector ( ${\bar{θ}}_{1}^{'}$ ) from Algorithm 1

Output:

Estimated DOA vector, ${\bar{θ}}^{'}$
Estimated CF vector, ${\bar{f}}^{'}$

Algorithm:

Generate orthogonal complement projector ${\bar{\bar{Π}}}_{x}$ and ${\bar{\bar{Π}}}_{z}$ as in (11):
${\bar{\bar{Π}}}_{x} = {\bar{\bar{I}}}_{2 N_{s} + 1} - {\bar{\bar{B}}}_{x}^{'} \cdot {({\bar{\bar{B}}}_{x}^{' †} \cdot {\bar{\bar{B}}}_{x}^{'})}^{- 1} \cdot {\bar{\bar{B}}}_{x}^{' †},$ ${\bar{\bar{Π}}}_{z} = {\bar{\bar{I}}}_{2 N_{s} + 1} - {\bar{\bar{B}}}_{z}^{'} \cdot {({\bar{\bar{B}}}_{z}^{' †} \cdot {\bar{\bar{B}}}_{z}^{'})}^{- 1} \cdot {\bar{\bar{B}}}_{z}^{' †}$
Define ${\bar{y}}_{x 1}^{Φ}$ and ${\bar{y}}_{x 2}^{Φ}$ as the first $2 N_{s}$ elements and last $2 N_{s}$ elements, respectively, of ${\bar{y}}_{x}^{Φ}$ ; ${\bar{y}}_{z 1}^{Φ}$ and ${\bar{y}}_{z 2}^{Φ}$ are the first $2 N_{s}$ elements and last $2 N_{s}$ elements, respectively, of ${\bar{y}}_{z}^{Φ}$
Retrieve residual signals as in (12):
${\bar{y}}_{x 1}^{″} = {\bar{\bar{Π}}}_{x} \cdot {\bar{y}}_{x 1}^{Φ},$ ${\bar{y}}_{x 2}^{″} = {\bar{\bar{Π}}}_{x} \cdot {\bar{y}}_{x 2}^{Φ},$ ${\bar{y}}_{z 1}^{″} = {\bar{\bar{Π}}}_{z} \cdot {\bar{y}}_{z 1}^{Φ},$ ${\bar{y}}_{z 2}^{″} = {\bar{\bar{Π}}}_{z} \cdot {\bar{y}}_{z 2}^{Φ}$
Estimate cross-covariance matrices as in (15):
${\bar{\bar{R}}}_{1} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{x 1}^{″} [q] {\bar{y}}_{z 1}^{″ †} [q],$ ${\bar{\bar{R}}}_{2} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{x 2}^{″} [q] {\bar{y}}_{z 1}^{″ †} [q],$ ${\bar{\bar{R}}}_{3} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{x 1}^{″} [q] {\bar{y}}_{z 2}^{″ †} [q],$
${\bar{\bar{R}}}_{4} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{y}}_{x 2}^{″} [q] {\bar{y}}_{z 2}^{″ †} [q],$
Perform SVD to have $\bar{\bar{R}} = {[{\bar{\bar{R}}}_{1}, {\bar{\bar{R}}}_{2}, {\bar{\bar{R}}}_{3}, {\bar{\bar{R}}}_{4}]}^{t} = \bar{\bar{U}} \cdot \bar{\bar{Σ}} \cdot \bar{\bar{V}}$
Define ${\bar{\bar{U}}}_{s}$ as the first $M - M_{c}^{'}$ columns of $\bar{\bar{U}}$ , where $M_{c}^{'}$ is the size of ${\bar{θ}}_{1}^{'}$
Define ${\bar{\bar{U}}}_{s 1}$ as the first $2 N_{s}$ rows of ${\bar{\bar{U}}}_{s}$ , ${\bar{\bar{U}}}_{s 2}$ as the second $2 N_{s}$ rows of ${\bar{\bar{U}}}_{s}$ . ${\bar{\bar{U}}}_{s 3}$ and ${\bar{\bar{U}}}_{s 4}$ are defined similarly.
Compute (19):
${\bar{\bar{U}}}_{12} = {\bar{\bar{U}}}_{s 1}^{- 1} \cdot {\bar{\bar{U}}}_{s 2},$ ${\bar{\bar{U}}}_{13} = {\bar{\bar{U}}}_{s 1}^{- 1} \cdot {\bar{\bar{U}}}_{s 3},$ ${\bar{\bar{U}}}_{14} = {\bar{\bar{U}}}_{s 1}^{- 1} \cdot {\bar{\bar{U}}}_{s 4}$
Perform eigenvalue decomposition (EVD) to have $\bar{\bar{Ω}} = ({\bar{\bar{U}}}_{12} + {\bar{\bar{U}}}_{13} + {\bar{\bar{U}}}_{14}) / 3 = {\bar{\bar{T}}}^{- 1} \cdot \bar{\bar{Λ}} \cdot \bar{\bar{T}}$
Retrieve matrices ${\bar{\bar{Ψ}}}_{x}$ and ${\bar{\bar{Ψ}}}_{z}$ as in (19):
${\bar{\bar{Ψ}}}_{x} = \bar{\bar{T}} \cdot {\bar{\bar{U}}}_{12} \cdot {\bar{\bar{T}}}^{- 1},$ ${\bar{\bar{Ψ}}}_{z} = {(\bar{\bar{T}} \cdot {\bar{\bar{U}}}_{13} \cdot {\bar{\bar{T}}}^{- 1})}^{†}$
Compute remaining DOAs and CFs as in (21) to have ${\bar{θ}}_{2}^{'}$ and ${\bar{f}}_{2}^{'}$ .
Record ${\bar{θ}}^{'} = {[{\bar{θ}}_{1}^{' t}, {\bar{θ}}_{2}^{' t}]}^{t}$ and ${\bar{f}}^{'} = {[f_{0}, \dots, f_{0}, {\bar{f}}_{2}^{' t}]}^{t}$

4. Simulations and Discussions

4.1. Simulation Setup

Two CPA(5, 6), each composed of 15 physical sensors, are deployed along x and z axes, respectively. The unit spacing is

d = 3

cm, which is

λ / 2

at 5 GHz. A total of

M = 30

sources are considered, including

M_{u}

CF-unknown sources and

M_{c}

CF-known sources. For convenience, we assume the first

M_{c}

sources have known carrier frequency of

f_{0}

. Half of the DOAs are distributed over

[20^{°}, 70^{°}]

and the other half over

[- 70^{°}, - 20^{°}]

, all at uniform spacing.

Each source signal is assumed to display quasi-stationary property [11]. The time-frame length of each source signal is set to

L = 500

. A total of 200 Monte Carlo realizations are simulated in each scenario. The signal-to-noise ratio (SNR) is defined as

\begin{matrix} SNR = \frac{\frac{1}{L Q} \sum_{ℓ = 1}^{L Q} (E \{{∥{\bar{\bar{A}}}_{x} \cdot \bar{w} [ℓ]∥}^{2}\} + E \{{∥{\bar{\bar{A}}}_{z} \cdot \bar{w} [ℓ]∥}^{2}\})}{E \{{∥{\bar{n}}_{x} [ℓ]∥}^{2}\} + E \{{∥{\bar{n}}_{z} [ℓ]∥}^{2}\}} \end{matrix}

(21)

The CF-known sources emit at

f_{0} = 5

GHz. The CFs of the remaining sources are chosen from

[0, f_{N} - B]

, with

f_{N} = 6

GHz and

B = 50

MHz, where the spacing between two adjacent channels is wider than B. To evaluate the performance of the proposed algorithm, define root-mean-square error (RMSE) as

\begin{matrix} {RMSE}_{θ} = \frac{1}{180} \sqrt{\frac{1}{M} \sum_{m = 1}^{M} {| {\tilde{θ}}_{m} - θ_{m} |}^{2}} \\ {RMSE}_{f} = \frac{1}{f_{N}} \sqrt{\frac{1}{M_{u}} \sum_{m = M_{c} + 1}^{M} {| {\tilde{f}}_{m} - f_{m} |}^{2}} \end{matrix}

where

{\tilde{θ}}_{m}

and

θ_{m}

are the estimated DOA and actual DOA, respectively, of source m,

{\tilde{f}}_{m}

and

f_{m}

are the estimated CF and actual CF, respectively, of the m-th CF-unknown source.

4.2. Cramer–Rao Bound

In [24,25], the Cramer–Rao bound (CRB) of DOA estimation was analyzed, assuming that there are fewer sources than antennas. In [26,27], the CRB of 1D DOA estimation in [24,25] was extended to 2D DOA estimation. In this work, the approach in [24,25] is extended to derive the CRB of joint DOA and frequency estimation. By concatenating

\bar{x} [ℓ]

and

\bar{z} [ℓ]

in (4), a signal vector is formed as

\begin{matrix} \bar{g} [ℓ] = [\begin{matrix} \bar{x} [ℓ] \\ \bar{z} [ℓ] \end{matrix}] = [\begin{matrix} {\bar{\bar{A}}}_{x} \\ {\bar{\bar{A}}}_{z} \end{matrix}] \cdot \bar{w} [ℓ] + [\begin{matrix} {\bar{n}}_{x} [ℓ] \\ {\bar{n}}_{z} [ℓ] \end{matrix}] = {\bar{\bar{A}}}_{g} \cdot \bar{w} [ℓ] + {\bar{n}}_{g} [ℓ] \end{matrix}

(22)

from which its covariance matrix is estimated as

\begin{matrix} {\bar{\bar{R}}}_{g} = \frac{1}{L Q} \sum_{ℓ = 1}^{L Q} \bar{g} [ℓ] {\bar{g}}^{†} [ℓ] = \frac{1}{Q} \sum_{q = 1}^{Q} \frac{1}{L} \sum_{ℓ = 1}^{L} \bar{g} [(q - 1) L + ℓ] {\bar{g}}^{†} [(q - 1) L + ℓ] \\ = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{\bar{A}}}_{g} \cdot {\bar{\bar{D}}}^{q} \cdot {\bar{\bar{A}}}_{g}^{†} + σ_{n}^{' 2} {\bar{\bar{I}}}_{2 N} = {\bar{\bar{A}}}_{g} \cdot {\bar{\bar{D}}}_{g} \cdot {\bar{\bar{A}}}_{g}^{†} + σ_{n}^{' 2} {\bar{\bar{I}}}_{2 N} \end{matrix}

(23)

where

{\bar{\bar{D}}}_{g} = \frac{1}{Q} \sum_{q = 1}^{Q} {\bar{\bar{D}}}^{q} =

diag

{μ_{1}, \dots, μ_{M}}

. The parameters are listed in a vector as

\begin{matrix} \bar{η} = {[θ_{1}^{'}, \dots, θ_{M}^{'}, f_{0}^{'}, f_{M_{c} + 1}^{'}, \dots, f_{M}^{'}, μ_{1}, \dots, μ_{M}, σ_{n}^{' 2}]}^{t} \end{matrix}

(24)

where

θ_{m}^{'} = θ_{m} / 180

and

f_{m}^{'} = f_{m} / f_{N}

The CRB matrix is formulated as [24,25]

\begin{matrix} \bar{\bar{B}} = \frac{σ_{n}^{' 2}}{L Q} {\{Re \{({\bar{\bar{F}}}^{†} \cdot {\bar{\bar{Π}}}_{g} \cdot \bar{\bar{F}}) ° {\bar{\bar{D}}}_{g}^{' t}\}\}}^{- 1} \end{matrix}

(25)

where

σ_{n}^{' 2}

is the noise power defined in (4),

L Q

is the total number of snapshots, Re

{}

denotes the real part of its complex argument,

{\bar{\bar{Π}}}_{g} = {\bar{\bar{I}}}_{2 N} - {\bar{\bar{A}}}_{g} \cdot {({\bar{\bar{A}}}_{g}^{†} \cdot {\bar{\bar{A}}}_{g})}^{- 1} \cdot {\bar{\bar{A}}}_{g}^{†}

is the orthogonal complement projection based on

{\bar{\bar{A}}}_{g}

, ∘ is the Hadamard product and

{\bar{\bar{D}}}_{g}^{'} = [\begin{matrix} {\bar{\bar{D}}}_{g} & {\bar{\bar{D}}}_{g} \\ {\bar{\bar{D}}}_{g} & {\bar{\bar{D}}}_{g} \end{matrix}]

\bar{\bar{F}} = [{\bar{\bar{A}}}_{θ}^{'}, {\bar{\bar{A}}}_{f}^{'}]

, with

\begin{matrix} {\bar{\bar{A}}}_{θ}^{'} = [\frac{\partial \bar{α} (f_{0}, θ_{1})}{\partial θ_{1}}, \dots, \frac{\partial \bar{α} (f_{0}, θ_{M_{c}})}{\partial θ_{M_{c}}}, \frac{\partial \bar{α} (f_{M_{c} + 1}, θ_{M_{c} + 1})}{\partial θ_{M_{c} + 1}}, \dots, \frac{\partial \bar{α} (f_{M}, θ_{M})}{\partial θ_{M}}] \\ {\bar{\bar{A}}}_{f}^{'} = [\frac{\partial \bar{α} (f_{0}, θ_{1})}{\partial f_{0}}, \dots, \frac{\partial \bar{α} (f_{0}, θ_{M_{c}})}{\partial f_{0}}, \frac{\partial \bar{α} (f_{M_{c} + 1}, θ_{M_{c} + 1})}{\partial f_{M_{c} + 1}}, \dots, \frac{\partial \bar{α} (f_{M}, θ_{M})}{\partial f_{M}}] \end{matrix}

(26)

where

\bar{α} (f, θ) = {[{\bar{α}}_{x}^{t} (f, θ), {\bar{α}}_{z}^{t} (f, θ)]}^{t}

{\bar{α}}_{x} (f, θ) = {[e^{j 2 π f χ_{1} τ_{x} (θ)}, \dots, e^{j 2 π f χ_{N} τ_{x} (θ)}]}^{t}

{\bar{α}}_{z} (f, θ) = {[e^{j 2 π f χ_{1} τ_{z} (θ)}, \dots, e^{j 2 π f χ_{N} τ_{z} (θ)}]}^{t}

The SNR defined in (21) can be rewritten as

\begin{matrix} SNR = \frac{\sum_{m = 1}^{M_{c}} μ_{m} {∥\bar{α} (f_{0}, θ_{m})∥}_{2}^{2} + \sum_{m = M_{c} + 1}^{M} μ_{m} {∥\bar{α} (f_{m}, θ_{m})∥}_{2}^{2}}{2 σ_{n}^{' 2}} \end{matrix}

(27)

and the CRBs in DOA and CF estimations can be represented as

\begin{matrix} {CRB}_{θ} = \sqrt{\frac{1}{M} \sum_{m = 1}^{M} B_{m m}} \\ {CRB}_{f} = \sqrt{\frac{1}{M_{u}} \sum_{m = M + M_{c} + 1}^{2 M} B_{m m}} \end{matrix}

(28)

4.3. Maximum Detectable Number of Sources

In the first scenario, we compare the maximum number of detectable signal sources when there is only one source with unknown CF (1.15 GHz). Figure 3 shows the success probability of DOA and CF estimation under different numbers of signal sources, where the success probability is defined as the number of Monte Carlo realizations with RMSE

< 5

% in both DOA and CF estimation, divided by the total number of realizations. It is observed that the success probability with conventional JE [7,8] is acceptable when

M \leq 5

, and drops dramatically when

M > 5

. On the other hand, the success probability with the proposed PJE remains about 90% even when

M > 5

Figure 4 shows the RMSE of DOA and CF estimation, respectively, under different numbers of signal sources. The RMSE with conventional JE algorithm is higher than 5% when

M > 2

, while that with the proposed PJE algorithm remains about 1% even at

M = 8

In the second scenario, we try to find the maximum detectable number of CF-unknown signal sources, under a fixed number of sources. The success probability and RMSE with the proposed PJE algorithm are shown in Figure 5 and Figure 6, respectively. The RMSE with the conventional JE algorithm is higher than 7%. The proposed two-stage algorithm works reasonably well when

M_{u} \leq 11

, and the RMSE increases dramatically when

M_{u} \geq 13

, indicating that the maximum number of CF-unknown sources is about 13 in this case.

Next, we try to find the maximum allowable number of signal sources, including 10 CF-unknown sources. Figure 7 shows that the success probability decreases rapidly when

M > 36

and drops to 0 at

M = 40

. Figure 8 shows that the RMSE of DOA and CF estimation increases significantly when M is increased from 30 to 32, and stays at the same level before M reaches 39.

4.4. Detail of Proposed Two-Stage Algorithm

Figure 9 shows the normalized SS-MUSIC spectrum with Algorithm 1. It is observed that the blue peaks are aligned with the green peaks at the actual DOAs of CF-known sources in both cases of

M_{u} = 13

and

M_{u} = 14

. All the CF-known sources and some of the CF-unknown sources are accurately estimated.

Figure 10 shows the estimated DOA and CF, with

M_{u} = 13

and 14, respectively. In both cases, the DOA and CF of CF-known sources are accurately estimated. Figure 10a shows that all the CF-unknown sources are accurately estimated, while Figure 10b shows that some CF-unknown sources are poorly estimated, which is consistent with the observation in Figure 6 that the maximum number of detectable CF-unknown sources is 13.

Figure 11a,b show the estimated DOA and CF by using the proposed PJE algorithm and the conventional JE algorithm [8], respectively, with

M_{u} = 10

and

M = 30

. It is observed that the proposed PJE algorithm can detect all the sources accurately, while the RMSE of either DOA or CF with the conventional JE algorithm exceeds 20%. It is interesting to observe in Figure 11b that the error of CF-known sources is higher than that of CF-unknown sources. Since the received signals are contributed by several CF-known sources directionally nearby the CF-unknown sources, a high resolution of signal phase is required to resolve these sources. The joint diagonalization in the conventional JE algorithm probably degrade such phase resolution. Figure 11c shows the results by using the proposed PJE algorithm with

M_{u} = 13

. All the sources are accurately estimated, with slightly larger error than in the case with

M_{u} = 10

The SS-MUSIC algorithm in [17] has been extended to deal with multiple sources, some with known CF and others with unknown CF. By matching the results from x-CPA and z-CPA, the DOAs of CF-known sources are estimated first. Then, the CF-known signals are removed from the received signals by using (11)–(13) in the proposed PJE algorithm. Finally, the conventional JE algorithm is applied to jointly estimate the DOA and CF of the remaining CF-unknown signals.

4.5. Robustness of Proposed Two-Stage Algorithm

Figure 12 shows the RMSE of DOA and CF estimation with the proposed PJE algorithm, under different SNR. The RMSE is well below 1% at SNR

> 0

dB if

M_{u} = 10

, and is below 1% at SNR

> 5

dB if

M_{u} = 13

In the third scenario, we try to find the minimum angular spacing between two signal sources that can be detected by using the proposed two-stage algorithm. Consider the case of

M = 3

, with the incident angles of

- 60^{°}

5^{°}

and

{(5 + x)}^{°}

, respectively, where x is in the range of

[0.1, 1]

. Figure 13 shows that the proposed PJE algorithm can detect two signals as close as

0 . 3^{°}

apart, with success probability higher than 95%, as compared to

{0.5}^{°}

if conventional JE algorithm is used.

Figure 14 shows that the RMSEs of both DOA and CF estimation by using the proposed PJE algorithm are smaller than their counterparts with conventional JE algorithm when the angular spacing is smaller than

{0.5}^{°}

5. Conclusions

A two-stage method is proposed to estimate the DOA of multiple signal sources by processing the received signals contributed by multiple CF-known sources and multiple CF-unknown sources. In the first stage, the intermediate results derived from the received signals of two orthogonal CPAs are matched to estimate the DOA of CF-known sources. In the second stage, the proposed PJE algorithm is applied to the residual signals, after removing the contribution of CF-known sources, to estimate the DOA and CF of the remaining sources. Simulation results show that by applying the proposed two-stage method to received signals from two orthogonal CPA(5, 6), both DOA and CF of 30 signal sources, 13 of which with unknown CF, can be accurately estimated.

Author Contributions

Conceptualization: K.-C.H. and J.-F.K.; methodology: K.-C.H. and J.-F.K.; Software: K.-C.H.; validation: K.-C.H. and J.-F.K.; formal analysis: K.-C.H. and J.-F.K.; writing-original draft preparation: K.-C.H.; writing-review and editing: J.-F.K.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

Schmidt, R. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propag. 1986, 34, 276–280. [Google Scholar] [CrossRef]
Roy, R.; Kailath, T. ESPRIT-estimation of signal parameters via rotational invariance techniques. IEEE Trans. Acoust. Speech Signal Process. 1989, 37, 984–995. [Google Scholar] [CrossRef]
Kim, J.M.; Lee, O.K.; Ye, J.C. Compressive MUSIC: Revisiting the link between compressive sensing and array signal processing. IEEE Trans. Inf. Theory 2012, 58, 278–301. [Google Scholar] [CrossRef]
Qin, S.; Zhang, Y.D.; Amin, M.G. Generalized coprime array configurations for direction-of-arrival estimation. IEEE Trans. Signal Process. 2015, 63, 1377–1390. [Google Scholar] [CrossRef]
Haardt, M.; Nossek, J.A. 3-D unitary ESPRIT for joint angle and carrier estimation. IEEE Int. Conf. Acoust. Speech Signal Process. 1997, 255–258. [Google Scholar] [CrossRef]
Dong, W.; Zhang, X.; Li, J.; Bai, J. Improved ESPRIT method for joint direction-of-arrival and frequency estimation using multiple-delay output. Int. Antennas J. Propag. 2012, 2012, 309269. [Google Scholar]
Stein, S.; Yair, O.; Cohen, D.; Eldar, Y.C. Joint spectrum sensing and direction of arrival recovery from sub-Nyquist samples. In Proceedings of the 2015 IEEE 16th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC), Stockholm, Sweden, 28 June–1 July 2015; pp. 331–335. [Google Scholar]
Stein, S.; Yair, O.; Cohen, D.; Eldar, Y.C. CaSCADE: Compressed carrier and DOA estimation. IEEE Trans. Signal Process. 2017, 65, 2645–2658. [Google Scholar] [CrossRef]
Lemma, A.N.; Van Der Veen, A.J.; Deprettere, E.F. Analysis of joint angle-frequency estimation using ESPRIT. IEEE J. Sel. Top. Signal Process. 2003, 51, 1264–1283. [Google Scholar] [CrossRef] [Green Version]
Cui, C.; Wu, W.; Wang, W. Carrier frequency and DOA estimation of sub-Nyquist sampling multi-band sensor signals. IEEE Sens. J. 2017, 17, 7470–7478. [Google Scholar] [CrossRef]
Mao, W.K.; Hsieh, T.H.; Chi, C.Y. DOA estimation of quasi-stationary signals with less sensors than sources and unknown spatial noise covariance: A Khatri–Rao subspace approach. IEEE Trans. Signal Process. 2010, 58, 2168–2180. [Google Scholar] [CrossRef]
Wang, B.; Wang, W.; Gu, Y.; Lei, S. Underdetermined DOA estimation of quasi-stationary signals using a partly-calibrated array. Sensors 2017, 17, 702. [Google Scholar] [CrossRef] [PubMed]
Hsu, K.C.; Kiang, J.F. DOA estimation of quasi-stationary signals using a partly-calibrated uniform linear array with fewer sensors than sources. Prog. Electromagn. Res. M 2018, 63, 185–193. [Google Scholar] [CrossRef]
Pal, P.; Vaidyanathan, P.P. Nested arrays: A novel approach to array processing with enhanced degrees of freedom. IEEE Trans. Signal Process. 2010, 58, 4167–4181. [Google Scholar] [CrossRef]
Vaidyanathan, P.P.; Pal, P. Sparse sensing with co-prime samplers and arrays. IEEE Trans. Signal Process. 2011, 59, 573–586. [Google Scholar] [CrossRef]
Pal, P.; Vaidyanathan, P.P. Coprime sampling and the MUSIC algorithm. In Proceedings of the 2011 Digital Signal Processing and Signal Processing Education Meeting (DSP/SPE), Sedona, AZ, USA, 4–7 January 2011; pp. 289–294. [Google Scholar]
Shen, Q.; Liu, W.; Cui, W.; Wu, S. Extension of co-prime arrays based on the fourth-order difference co-array concept. IEEE Signal Process. Lett. 2016, 23, 615–619. [Google Scholar] [CrossRef]
Hsu, K.C.; Kiang, J.F. DOA estimation using triply primed arrays based on fourth-order statistics. Prog. Electromagn. Res. M 2018, 67, 55–64. [Google Scholar] [CrossRef]
Chen, K.H.; Kiang, J.F. Coupling characterization of a linear dipole array to improve direction-of-arrival estimation. IEEE Trans. Antennas Propag. 2015, 63, 5056–5062. [Google Scholar] [CrossRef]
Adve, R.S.; Sarkar, T.K. Compensation for the effects of mutual coupling on direct data domain adaptive algorithms. IEEE Trans. Antennas Propag. 2000, 48, 86–94. [Google Scholar] [CrossRef] [Green Version]
Dandekar, K.R.; Ling, H.; Xu, G. Experimental study of mutual coupling compensation in smart antenna applications. IEEE Trans. Wirel. Commun. 2002, 1, 480–487. [Google Scholar] [CrossRef]
Morales-Jimenez, D.; Raymond Louie, H.Y.; McKay, M.R.; Chen, Y. Analysis and design of multiple-antenna cognitive radios with multiple primary user signals. IEEE Trans. Signal Process. 2015, 63, 4925–4939. [Google Scholar] [CrossRef]
Liu, C.L.; Vaidyanathan, P. Remarks on the spatial smoothing step in coarray MUSIC. IEEE Signal Process. Lett. 2015, 22, 1438–1442. [Google Scholar] [CrossRef]
Stoica, P.; Nehorai, A. MUSIC, maximum likelihood, and Cramer–Rao bound. IEEE Trans. Acoust. Speech Signal Process. 1989, 37, 720–741. [Google Scholar] [CrossRef]
Stoica, P.; Nehorai, A. Performance study of conditional and unconditional direction-of-arrival estimation. IEEE Trans. Acoust. Speech Signal Process. 1990, 38, 1783–1795. [Google Scholar] [CrossRef]
Chen, H.; Hou, P.; Wang, Q.; Huang, L.; Yan, W.-Q. Cumulants-based Toeplitz matrices reconstruction method for 2-D coherent DOA estimation. IEEE Sens. J. 2014, 14, 2824–2832. [Google Scholar] [CrossRef]
Wang, X.; Mao, X.; Wang, Y.; Zhang, N.; Li, B. A novel 2-D coherent DOA estimation method based on dimension reduction sparse reconstruction for orthogonal arrays. Sensors 2016, 16, 1496. [Google Scholar] [CrossRef] [PubMed]

Figure 1. Orthogonal CPAs deployed along x and z axes, respectively, to estimate DOAs of M sources.

Figure 2. Configuration of CPA(

N_{1}

N_{2}

N_{1} = 5

N_{2} = 6

, ■: first ULA at spacing

N_{1} d

, ▲: second ULA at spacing

N_{2} d

, •: physical array of CPA (

N_{1}

N_{2}

), •: second-order virtual array, ×: second-order holes, •: fourth-order virtual array, ×: fourth-order holes.

Figure 2. Configuration of CPA(

N_{1}

N_{2}

N_{1} = 5

N_{2} = 6

, ■: first ULA at spacing

N_{1} d

, ▲: second ULA at spacing

N_{2} d

, •: physical array of CPA (

N_{1}

N_{2}

), •: second-order virtual array, ×: second-order holes, •: fourth-order virtual array, ×: fourth-order holes.

Figure 3. Success probability of DOA and CF estimation under different numbers of signals, SNR = 10 dB,

Q = 300

M_{u} = 1

. ───: proposed PJE algorithm, ───: conventional JE algorithm.

Figure 3. Success probability of DOA and CF estimation under different numbers of signals, SNR = 10 dB,

Q = 300

M_{u} = 1

. ───: proposed PJE algorithm, ───: conventional JE algorithm.

Figure 4. RMSE of DOA and CF estimation under different numbers of signals, SNR

= 10

dB,

Q = 300

M_{u} = 1

. ───: RMSE of DOA with proposed PJE algorithm,

- - -

: RMSE of CF with proposed PJE algorithm, ───: RMSE of DOA with conventional JE algorithm,

- - -

: RMSE of CF with conventional JE algorithm, ───: CRB of DOA estimation,

- - -

: CRB of CF estimation.

Figure 4. RMSE of DOA and CF estimation under different numbers of signals, SNR

= 10

dB,

Q = 300

M_{u} = 1

. ───: RMSE of DOA with proposed PJE algorithm,

- - -

: RMSE of CF with proposed PJE algorithm, ───: RMSE of DOA with conventional JE algorithm,

- - -

: RMSE of CF with conventional JE algorithm, ───: CRB of DOA estimation,

- - -

: CRB of CF estimation.

Figure 5. Success probability of DOA and CF estimation with proposed PJE algorithm, under different numbers of signals with unknown CF, SNR

= 10

dB,

Q = 300

M = 30

Figure 5. Success probability of DOA and CF estimation with proposed PJE algorithm, under different numbers of signals with unknown CF, SNR

= 10

dB,

Q = 300

M = 30

Figure 6. RMSE of DOA and CF estimation with proposed PJE algorithm, under different numbers of signals with unknown CF, SNR

= 10

dB,

Q = 300

M = 30

. ───: RMSE of DOA,

- - -

: RMSE of CF.

Figure 6. RMSE of DOA and CF estimation with proposed PJE algorithm, under different numbers of signals with unknown CF, SNR

= 10

dB,

Q = 300

M = 30

. ───: RMSE of DOA,

- - -

: RMSE of CF.

Figure 7. Success probability of DOA and CF estimation with proposed PJE algorithm, under different numbers of signal sources, SNR

= 10

dB,

Q = 300

M_{u} = 10

Figure 7. Success probability of DOA and CF estimation with proposed PJE algorithm, under different numbers of signal sources, SNR

= 10

dB,

Q = 300

M_{u} = 10

Figure 8. RMSE of DOA and CF estimation with proposed PJE algorithm, under different numbers of signal sources, SNR

= 10

dB,

Q = 300

M_{u} = 10

. ───: RMSE of DOA,

- - -

: RMSE of CF.

Figure 8. RMSE of DOA and CF estimation with proposed PJE algorithm, under different numbers of signal sources, SNR

= 10

dB,

Q = 300

M_{u} = 10

. ───: RMSE of DOA,

- - -

: RMSE of CF.

Figure 9. Normalized SS-MUSIC spectrum with Algorithm 1, (a)

M_{u} = 13

and (b)

M_{u} = 14

, SNR = 10 dB,

Q = 300

M = 30

. ───: actual DOA of CF-known sources, ───: actual DOA of CF-unknown sources,

- - -

: estimated DOA with x-axis CPA,

- - -

: estimated DOA with z-axis CPA, ×: matched DOA.

Figure 9. Normalized SS-MUSIC spectrum with Algorithm 1, (a)

M_{u} = 13

and (b)

M_{u} = 14

, SNR = 10 dB,

Q = 300

M = 30

. ───: actual DOA of CF-known sources, ───: actual DOA of CF-unknown sources,

- - -

: estimated DOA with x-axis CPA,

- - -

: estimated DOA with z-axis CPA, ×: matched DOA.

Figure 10. Joint estimation of DOA and CF with proposed PJE algorithm, (a)

M_{u} = 13

and (b)

M_{u} = 14

, SNR = 10 dB,

Q = 300

M = 30

. •: CF-known sources, •: CF-unknown sources, ×: estimated sources.

Figure 10. Joint estimation of DOA and CF with proposed PJE algorithm, (a)

M_{u} = 13

and (b)

M_{u} = 14

, SNR = 10 dB,

Q = 300

M = 30

. •: CF-known sources, •: CF-unknown sources, ×: estimated sources.

Figure 11. Joint estimation of DOA and CF, SNR = 10 dB,

Q = 300

M = 30

. (a) Proposed PJE algorithm with

M_{u} = 10

, (b) conventional JE algorithm with

M_{u} = 10

, (c) proposed PJE algorithm with

M_{u} = 13

. •: CF-known sources, •: CF-unknown sources, ×: estimated sources.

Figure 11. Joint estimation of DOA and CF, SNR = 10 dB,

Q = 300

M = 30

. (a) Proposed PJE algorithm with

M_{u} = 10

, (b) conventional JE algorithm with

M_{u} = 10

, (c) proposed PJE algorithm with

M_{u} = 13

. •: CF-known sources, •: CF-unknown sources, ×: estimated sources.

Figure 12. RMSE of DOA and CF estimation versus SNR, with proposed PJE algorithm,

Q = 300

M = 30

. ───: RMSE of DOA with

M_{u} = 10

- - -

: RMSE of CF with

M_{u} = 10

- \cdot -

: RMSE of DOA with

M_{u} = 13

, ⋯: RMSE of CF with

M_{u} = 13

Figure 12. RMSE of DOA and CF estimation versus SNR, with proposed PJE algorithm,

Q = 300

M = 30

. ───: RMSE of DOA with

M_{u} = 10

- - -

: RMSE of CF with

M_{u} = 10

- \cdot -

: RMSE of DOA with

M_{u} = 13

, ⋯: RMSE of CF with

M_{u} = 13

Figure 13. Success probability of DOA and CF estimation versus angular spacing, SNR

= 10

dB,

Q = 300

M = 3

M_{u} = 1

. ───: proposed PJE algorithm, ───: conventional JE algorithm.

Figure 13. Success probability of DOA and CF estimation versus angular spacing, SNR

= 10

dB,

Q = 300

M = 3

M_{u} = 1

. ───: proposed PJE algorithm, ───: conventional JE algorithm.

Figure 14. RMSE of DOA and CF estimation versus angular spacing, SNR

= 10

dB,

Q = 300

M = 3

M_{u} = 1

. ───: RMSE of DOA with proposed PJE algorithm,

- - -

: RMSE of CF with proposed PJE algorithm, ───: RMSE of DOA with conventional JE algorithm,

- - -

: RMSE of CF with conventional JE algorithm, ───: CRB of DOA estimation,

- - -

: CRB of CF estimation.

Figure 14. RMSE of DOA and CF estimation versus angular spacing, SNR

= 10

dB,

Q = 300

M = 3

M_{u} = 1

. ───: RMSE of DOA with proposed PJE algorithm,

- - -

: RMSE of CF with proposed PJE algorithm, ───: RMSE of DOA with conventional JE algorithm,

- - -

: RMSE of CF with conventional JE algorithm, ───: CRB of DOA estimation,

- - -

: CRB of CF estimation.

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Hsu, K.-C.; Kiang, J.-F. Joint Estimation of DOA and Frequency of Multiple Sources with Orthogonal Coprime Arrays. Sensors 2019, 19, 335. https://doi.org/10.3390/s19020335

AMA Style

Hsu K-C, Kiang J-F. Joint Estimation of DOA and Frequency of Multiple Sources with Orthogonal Coprime Arrays. Sensors. 2019; 19(2):335. https://doi.org/10.3390/s19020335

Chicago/Turabian Style

Hsu, Kai-Chieh, and Jean-Fu Kiang. 2019. "Joint Estimation of DOA and Frequency of Multiple Sources with Orthogonal Coprime Arrays" Sensors 19, no. 2: 335. https://doi.org/10.3390/s19020335

APA Style

Hsu, K. -C., & Kiang, J. -F. (2019). Joint Estimation of DOA and Frequency of Multiple Sources with Orthogonal Coprime Arrays. Sensors, 19(2), 335. https://doi.org/10.3390/s19020335

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Joint Estimation of DOA and Frequency of Multiple Sources with Orthogonal Coprime Arrays

Abstract

1. Introduction

2. Signal Model

2.1. Sub-Nyquist Sampling Scheme

2.2. Second-Order Statistics

3. Joint Estimation of DOA and CF

3.1. DOA Estimation of Signal Sources with Known CF

3.2. Joint-ESPRIT in Projected Subspace

4. Simulations and Discussions

4.1. Simulation Setup

4.2. Cramer–Rao Bound

4.3. Maximum Detectable Number of Sources

4.4. Detail of Proposed Two-Stage Algorithm

4.5. Robustness of Proposed Two-Stage Algorithm

5. Conclusions

Author Contributions

Funding

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI