An Efficient Method for Joint Delay-Doppler Estimation of Moving Targets in Passive Radar

At the beginning of this part, we briefly give a blanket of settings on passive radar systems, and then develop an accurately matched filtering formula. The settings include: (i) the received OFDM signal comes from an un-coordinated but synchronized illuminator; (ii) the passive radar system only performs demodulation on data symbols, but not on forward error correction decoding.

The OFDM is a multicarrier modulation scheme which can be described as follows: suppose that there are $M$ blocks in transmitted data and there are $N$ orthogonal subcarriers in each block. Accordingly, the data symbols in $m$-th block and $n$-th subcarrier are denoted by $s_{m}(n)$. As the demodulation type is usually known, then the data-symbol $s_{m}(n)$ can be estimated and the matched filter can be implemented. Given a rectangular window of length $T$ at the receiver, the frequency spacing is denoted as $\bigtriangleup f:=1/T$, and the duration of each transmission block is denoted as $\bar{T}:=T+T_{\text{cp}}$, where $T_{\text{cp}}$ is the length of a cyclic prefix. Using these notations, the baseband signal over the $M$ blocks can be expressed as

where $x_{m}(t)$ is the baseband signal in the $m$-th block with $m\bar{T}-T_{\text{cp}}\leq t\leq m\bar{T}+T$, and $\xi(t)=1$ if $t\in[-T_{\text{cp}},T]$ and $0$ otherwise. Indexing the return of the $k$-th arrival and its associated Doppler shift and delay, the received signal and the received signal $y(t)$ takes the following form:

where $K$ is the number of paths, $w(t)$ is an additive noise, $A_{k}$ is an attenuations coefficient, $\tau_{k}$ is a delay, and $f_{k}$ is a Doppler frequency.

If taking Fourier transform on $y(t)$ in $m$-th block, the resulting signal in $n$-th subcarrier is formulated as:

Assuming that the largest possible delay is smaller than the cyclic prefix, then the term $\xi(t-\tau_{k}-m^{\prime}\bar{T})$ is $1$ in the integration interval $[m\bar{T},m\bar{T}+T]$ and $0$ otherwise. Usually, the integration time is almost on the order of a second, which means that when the product between $\bar{T}$ and the Doppler frequency is small compared to unity, we can approximate the phase rotation within one OFDM block as constant [3], that is,

For notational simplicity, we denote $\alpha_{k}:=A_{k}T$, $\phi_{k}:=f_{k}\bar{T}\in[0,1)$ and $\psi_{k}:=\Delta f\tau_{k}\in[0,1)$. Besides, we note that $\int_{m\bar{T}}^{m\bar{T}+T}e^{i2\pi(l-n)\Delta ft}dt=0$ if $l\neq n$. Then $r_{m}(n)$ can further be rewritten as

where

in which, $v_{m}(n)$ is assumed to be a complex Gaussian variable with zero mean and variance $\sigma^{2}$. For convenience, we denote the estimated data-symbol in $m$-th block as $\hat{s}_{m}(n)$. Then the demodulation error, denoted by $e_{m}(n)$, takes the following form

It should be noted that the demodulation error rate of a communication system is typically low under some normal operating conditions. Hence, most of $e_{m}(n)$ are zero and the noise caused by demodulation error is always sparse.

Let $\bm{\hat{S}}\in\mathbb{C}^{M\times N}$, $\bm{R}\in\mathbb{C}^{M\times N}$, $\bm{Z}\in\mathbb{C}^{M\times N}$, $\bm{E}\in\mathbb{C}^{M\times N}$ and $\bm{V}\in\mathbb{C}^{M\times N}$ be matrices whose $(m,n)$-th element are denoted as $\hat{s}_{m}(n)$, $r_{m}(n)$, $z_{m}(n)$, $e_{m}(n)$ and $v_{m}(n)$, respectively. Let $\bm{\phi}:=(\phi_{1},\phi_{2},...,\phi_{K})^{\top}$ and $\bm{\psi}:=(\psi_{1},\psi_{2},...,\psi_{K})^{\top}$. Accordingly, the response matrices are defined as $\bm{B}(\bm{\phi}):=(\bm{b}(\phi_{1}),\bm{b}(\phi_{2}),...,\bm{b}(\phi_{K}))\in\mathbb{C}^{M\times K}$ and $\bm{G}(\bm{\psi}):=(\bm{g}(\psi_{1}),\bm{g}(\psi_{2}),...,\bm{g}(\psi_{K}))\in\mathbb{C}^{N\times K}$, in which, $\bm{b}(\phi):=(1,e^{i2\pi\phi},\ldots,e^{i2\pi(M-1)\phi})^{\top}\in\mathbb{C}^{M}$ and $\bm{g}(\psi):=(1,e^{i2\pi\psi},\dots,e^{i2\pi(N-1)\psi})^{\top}\in\mathbb{C}^{N}$. Using these notations, the relation (2.5) can be rewritten in a matrix form, that is,

where $\bm{\alpha}:=(\alpha_{1},\alpha_{2},...,\alpha_{K})^{\top}\in\mathbb{C}^{K}$ and “$\odot$” is a Hadamard product. Utilizing the structure of the signal, the matrices in (2.6) can be vectorized to obtain a more concise and intuitive representation, that is,

where $\bm{\tilde{S}}=\text{diag}(\text{vec}(\bm{\hat{S}}))\in\mathbb{C}^{MN\times MN}$, $\bm{r}=\text{vec}(\bm{R})\in\mathbb{C}^{MN}$, $\bm{z}=\text{vec}(\bm{Z})\in\mathbb{C}^{MN}$, $\bm{e}=\text{vec}(\bm{E})\in\mathbb{C}^{MN}$ and $\bm{v}=\text{vec}(\bm{V})\in\mathbb{C}^{MN}$. In this formula, the term $\left(\bm{G}(\bm{\psi})^{C}\circ\bm{B}(\bm{\phi})\right)\in\mathbb{C}^{MN\times K}$ is a matrix whose $k$-th column has the form $\bm{g}(\psi_{k})^{C}\otimes\bm{b}(\phi_{k})$ in which $\otimes$ denotes Kronecker product.

From the work of Chandrasekaran et al. [8] that, the atomic-norm is used to characterize a sparse combination of sinusoids. Hence, define an atom as $\bm{a}(\phi,\psi)=\bm{g}(\psi)\otimes\bm{b}(\phi)$ where ‘$\otimes$’ denotes Kronecker product. Obviously, we have that $\bm{z}:=\sum_{k=1}^{K}\alpha_{k}\bm{a}(\phi_{k},\psi_{k})$. The atomic-norm can enforce sparsity in the atom set in the form of ${\cal A}:=\{\left.\bm{a}(\phi,\psi)\right|\phi\in[0,1),\psi\in[0,1)\}$ when it has been normalized. The definition of atomic-norm is given as follows:

At the end of this section, we quickly review the definition of proximal point operator defined in real space [24] which is used frequently at the subsequent analysis.

It is known that $\text{Prox}_{f}^{\mu}(y)$ is properly defined for every $y\in\mathbb{R}^{n}$ if $f$ is convex and proper. Specially, if $f(x)$ is an $\ell_{1}$-norm function, it can be expressed explicitly as $\text{Prox}_{f}^{\mu}(y)=\text{sign}(y)\odot\max\left\{|y|-\mu,0\right\}$. If $f(x)$ is an indicator function over a closed set $\mathcal{K}$, denoted by $\delta_{\mathcal{K}}(x)$, i.e., $0$ if $x\in\mathcal{K}$ and $\infty$ otherwise, then the proximal point mapping $\text{Prox}_{\delta_{\mathcal{K}}(\cdot)}^{\mu}(y)$ is equivalent to an orthogonal projection over ${\mathcal{K}}$, i.e., $\Pi_{\mathcal{K}}(y)$.

For convenience, we let $\bm{\Theta}:=\begin{pmatrix}\mathcal{T}(\bm{U})&\bm{z}\\ \bm{z}^{H}&\epsilon\end{pmatrix}\in\mathbb{C}^{(MN+1)\times(MN+1)}$, then model (1.4) is reformulated as

where ${\bm{z}}\in\mathbb{C}^{MN}$, ${\bm{e}}\in\mathbb{C}^{MN}$, $\bm{U}\in\mathbb{C}^{(2M-1)\times(2N-1)}$, $\epsilon\in\mathbb{C}$, the notation ${\mathcal{S}}_{+}^{MN+1}$ represents a symmetric and semi-positive matrix space with dimension $MN+1$. From the work of Tang et al. [29] that, the dual of (3.1) takes the following form

where $\bm{\nu}\in\mathbb{C}^{MN}$ is a dual variable, and $\|\bm{\nu}\|_{\cal A}^{H}=\sup_{\|\bm{z}\|_{\cal A}\leq 1}\langle\bm{\nu},\bm{z}\rangle_{\mathbb{R}}$ is called dual norm of $\|\cdot\|_{{\mathcal{A}}}$.

Let $\widetilde{\mathcal{L}}_{\rho}({\bm{e}},{\bm{z}},\epsilon,\bm{U},\bm{\Theta};\bm{\Gamma})$ be an augmented Lagrangian function associated with (3.1), where $\bm{\Gamma}\in\mathbb{C}^{(MN+1)\times(MN+1)}$ is a multiplier. The directly-extended ADMM employed by Zheng et al. [36] minimizes $\widetilde{\mathcal{L}}_{\rho}({\bm{z}},\epsilon,\bm{U},\bm{e},\bm{\Theta};\bm{\Gamma})$ in a simple Gauss-Seidel manner with respect to each variable by ignoring the relationship between them. More explicitly, with given $(\bm{z}^{k},\epsilon^{k},\bm{U}^{k},\bm{e}^{k},\bm{\Theta}^{k})$, the next point $({\bm{z}}^{k+1},\epsilon^{k+1},\bm{U}^{k+1},\bm{e}^{k+1},\bm{\Theta}^{k+1})$ is generated in the order of $\bm{z}\to\epsilon\to\bm{U}\to\bm{e}\to\bm{\Theta}$. As shown by Zheng et al. [36] that, each step involves solving a convex minimization problem which admits closed-form solutions by taking its favorable structures. Even though better performance has been achieved experimentally, the convergence of ADMM is still lack of certificate. For more details on this argument, one may refer to a counter-example of Chen et al. [9]. Because of this, it is ideal to design a convergent variant method which is at least as efficient as ADMM.

This part is devoted to designing a sGS based ADMM for solving (3.1). For convenience, we let $\bm{g}:={\bm{r}}-{\bm{e}}-\tilde{S}{\bm{z}}$, then (3.1) is rewritten equivalently as the following from

We note that although problem (3.3) has separable structures with regarding to objective function and constraints, it actually contains six blocks, which indicates that the traditional two-block ADMM will no longer be used. To address this issue, we partition these variables into two groups, for example, we can view $(\bm{e},\bm{g})$ as one group and $(\bm{z},\epsilon,\bm{U},\bm{\Theta})$ as another. For clearly catching this partition, we let

where ‘$\bm{0}$’ denotes a zero vector, and ‘$\bm{0}_{MN}$’ denotes a zero matrix with order $MN$. Using both notations, model (3.3) is transformed into the following form with respect to two variables $\bm{Y}$ and $\bm{W}$, that is,

where $\bm{I}$ is an identity matrix with an appropriate size. Obviously, we see that there are two nonsmooth terms in the objective function of (3.5), one is about $\bm{Y}_{1}$ in group $\bm{Y}$ and the other one is about $\bm{\Theta}$ in group $\bm{W}$, and each group possesses a “nonsmooth+smooth” structure.

Let $\rho>0$ be a penalty parameter and $\mathcal{L}_{\rho}({\bm{e}},\bm{g},{\bm{z}},\epsilon,\bm{U},\bm{\Theta};\bm{\beta},\bm{\Gamma})$ be an augmented Lagrangian function associated with (3.3), that is,

where $\bm{\beta}\in\mathbb{C}^{MN\times 1}$ and $\bm{\Gamma}\in\mathbb{C}^{(MN+1)\times(MN+1)}$ are multipliers. The sGS based ADMM minimizes $\mathcal{L}_{\rho}({\bm{e}},\bm{g},{\bm{z}},\epsilon,\bm{U},\bm{\Theta};\bm{\beta},\bm{\Gamma})$ alternatively with respect to two groups $(\bm{e},\bm{g})$ and $(\bm{z},\epsilon,\bm{U},\bm{\Theta})$, and then adopts sGS technique in each group. For example, in the $(\bm{e},\bm{g})$-group, it takes the order of $\bm{g}\to\bm{e}\to\bm{g}$ in the sense of computing $\bm{g}$ twice, while in the $(\bm{z},\epsilon,\bm{U},\bm{\Theta})$-group, it takes the order of $\bm{z}\to\epsilon\to\bm{U}\to\bm{\Theta}\to\bm{U}\to\epsilon\to\bm{z}$, i.e., computing $(\bm{z},\epsilon,\bm{U})$ twice. In light of these analysis, we get that, with given $(\bm{e}^{k},\bm{g}^{k})$, the next $(\bm{e}^{k+1},\bm{g}^{k+1})$ is obtained using the following scheme

Then, using $(\bm{e}^{k+1},\bm{g}^{k+1})$, the next point $(\bm{z}^{k+1},\epsilon^{k+1},\bm{U}^{k+1},\bm{\Theta}^{k+1})$ is achieved via the following scheme

While the latest iterates $(\bm{e}^{k+1},\bm{g}^{k+1},\bm{z}^{k+1},\epsilon^{k+1},\bm{U}^{k+1},\bm{\Theta}^{k+1})$ are obtained, the Lagrangian multipliers $(\bm{\beta},\bm{\Gamma})$ are updated according to the following schemes:

where $\varrho\in(0,(1+\sqrt{5})/2)$ is a step length.

Comparing the scheme (3.7)-(3.8) to ADMM of Zheng et al. [36], we see that the visible difference between them is that there are some additional subproblems. One may naturally think that this approach must lead to more computational burdens, and thus make it more inefficient. Next, we will show that this iterative scheme converges globally and that the total iterations would be obviously reduced, so that, these extra computational costs are almost negligible.

This part is devoted to solving the subproblems in (3.7)-(3.8). For convenience, we partition the variables $\bm{\Theta}\in\mathbb{C}^{(MN+1)\times(MN+1)}$ and $\bm{\Gamma}\in\mathbb{C}^{(MN+1)\times(MN+1)}$ into the following forms

where $\bm{\Theta}_{0}\in\mathbb{C}^{MN\times MN}$ and $\bm{\Gamma}_{0}\in\mathbb{C}^{MN\times MN}$ are submatrices, $\bm{\theta}_{1}\in\mathbb{C}^{MN}$ and $\bm{\gamma}\in\mathbb{C}^{MN}$ are column vectors, $\bar{\Theta}$ and $\bar{\Gamma}$ are scalars.

At the first place, we focus on solving subproblems regarding to the variables $\bm{g}$, $\bm{z}$, $\epsilon$, and $\bm{U}$ because ${\mathcal{L}}_{\rho}(\cdot)$ is smooth if $\bm{e}$ and $\bm{\Theta}$ are fixed. For simplicity, we abbreviate the augmented Lagrangian function $\mathcal{L}_{\rho}({\bm{e}},\bm{g},{\bm{z}},\epsilon,\bm{U},\bm{\Theta};\bm{\beta},\bm{\Gamma})$ as ${\mathcal{L}}_{\rho}$, and omit the superscripts appeared in (3.7) and (3.8). It is trivial to deduce that the partial derivative of ${\mathcal{L}}_{\rho}(\cdot)$ with respect to variables $\bm{g}$, $\bm{z}$, $\epsilon$, and $\bm{U}$ takes the following explicit form

For any $\bm{P}\in\mathbb{C}^{MN\times MN}$, we use the operator $\mathcal{S}_{l,j}(\bm{P})$ to return the $(l,j)$–th $M\times M$ submatrix $\bm{P}_{l,j}$, where the position of matrix $\bm{P}_{l,j}$ is corresponding to the $j$–th block matrix Toep($\bm{u}_{l}$) in $\mathcal{T}(\bm{U})$ (from left to right), and then use the operator ${\rm Tr}_{m}(\cdot)$ to output the sum of the elements locating in $m$-th sub-diagonal of an input matrix. We note that (3.11a)-(3.11c) are actually linear systems, which indicates that the iterates $(\bm{g}^{k+1},\bm{z}^{k+1},\epsilon^{k+1})$ can be got easily. We now pay our attention to the linear system (3.11d) to get the latest $\bm{U}^{k+1}$ in which its value is related to the position of the elements in $\bm{U}$. Let $\bm{I}_{(M,N)}$ be a $(2M-1)\times(2N-1)$ complex matrix whose element is $1$ at $(M,N)$-th position and $0$ otherwise. Using this matrix, the solution to $\bigtriangledown_{\bm{U}}\mathcal{L}_{\rho}=0$ can be expressed as the following unified form: