Joint Estimation for Time Delay and Direction of Arrival in Reconfigurable Intelligent Surface with OFDM

Abstract

Recently, the joint estimation for time delay (TD) and direction of arrival (DOA) has suffered from the high complexity of processing multi-dimensional signal models and the ineffectiveness of correlated/coherent signals. In order to improve this situation, a joint estimation method using orthogonal frequency division multiplexing (OFDM) and a uniform planar array composed of reconfigurable intelligent surface (RIS) is proposed. First, the time-domain coding function of the RIS is combined with the multi-carrier characteristic of the OFDM signal to construct the coded channel frequency response in tensor form. Then, the coded channel frequency response covariance matrix is decomposed by CANDECOMP/PARAFAC (CPD) to separate the signal subspaces of TD and DOA. Finally, we perform a one-dimensional (1D) spectral search for TD values and a two-dimensional (2D) spectral search for DOA values. Compared to previous efforts, this algorithm not only enhances the adaptability of coherent signals, but also greatly decreases the complexity. Simulation results indicate the robustness and effectiveness for the proposed algorithm in independent, coherent, and mixed multipath environments and low signal-to-noise ratio (SNR) conditions.

1. Introduction

Reconfigurable, intelligent surface (RIS) is a planar array consisting of reflective elements equipped with a time-domain encoding that can alter the frequency characteristics of incoming signals [1]. At present, the RIS has been widely researched in the fields of wireless communication [2], beamforming [3], and direction of arrival (DOA) estimation [4]. The RIS has the function of time-domain encoding, which can change the frequency characteristics of received signals [1]. Orthogonal frequency division multiplexing (OFDM) is a multi-carrier digital modulation technology which utilizes multiple parallel subcarriers to realize serial high-speed data communication. OFDM provides data transmission and positioning services for users and is widely used in wireless local area networks (WLAN) [5], 5G mobile communications [6], and time delay (TD) estimation [7]. A reasonable combination RIS functions and OFDM features can improve performance in a variety of wireless applications [8,9,10].

TD and DOA are important components of positioning systems, such as for indoor localization [11] and radar [12]. Like the time-domain narrowband signal model, the frequency-domain DOA estimation algorithm is studied by the multiple signal classification (MUSIC) method [13]. However, the capability of this algorithm is restricted to the size of the array aperture. An Algorithm [14] is used for TD estimation of OFDM signals, but it is not available in multipath environments. As researchers choose super-resolution methods such as the propagator method (PM) [15], MUSIC [16] and estimating signal parameters via the rotational invariance techniques (ESPRIT) [17], it is difficult to improve the accuracy of the estimations due to signal bandwidth constraints.

The joint estimation for TD and DOA by space–time parameter coupling not only improves the estimation performance, but also reduces the number of receiver nodes and improves the efficiency of the location system. Thus, the joint estimation has obvious strengths. Some researchers proposed a method [18] to solve TD under wideband signal conditions and then estimate DOA based on the trigonometric geometry TD inequality, but its performance does not improve significantly. In addition, a new method [19] was proposed to provide the frequency response of extended channels by using a Hadamard product, which achieves high accuracy joint estimation on OFDM systems. Nevertheless, the complexity of this approach is extremely high due to the demand for a full-field search.

The above algorithms are performed in an independent environment, which is necessary for the subspace approaches. For the coherent signal estimation, an approach [20] employing cyclostationary signals and a method [21] for separating coherent signal components have been proposed, both of which lead to more complex systems that are hard to implement in a joint estimation model. In addition, the method in [22] applies smoothing preprocessing to joint estimation to improve the adaptability in a coherent environment. However, smoothing preprocessing results in aperture loss and limited estimation accuracy.

The main efforts and results of this paper are summarized in the following points.

(i)

We present a joint estimation algorithm for TD and DOA. The algorithmic model integrates the RIS array response and the OFDM subcarrier response to build a coded channel frequency response (Coded-Response). The algorithm achieves excellent estimation performance under low signal-to-noise ratio (SNR) conditions.

(ii)

Due to the reasonable combination of the RIS time-domain coding function and the OFDM multi-subcarrier features, we construct frequency asymmetry of the space–time phase. This method reduces the singularity of the signal correlation matrix and thus effectively suppresses the coherence signal.

(iii)

The covariance matrix of Coded-Response is reconstructed by using the structural advantages of tensor. Further use of CANDECOMP/PARAFAC to decompose the covariance matrix results in corresponding TD and DOA signal subspaces, avoiding multi-dimensional spectral peak search and greatly reducing complexity.

(iv)

Compared with the existing RIS- and OFDM-based localization algorithms [23,24], the proposed algorithm can obtain the required parameters for localization based on a single station node in a coherent environment. Simulation results show that the proposed algorithm avoids the aperture loss of current smoothing algorithms [22,25] when processing coherent signals and thus has more advantages in terms of estimation accuracy.

The rest of the paper is summarized below. First, we describe the details of the signal model in Section 2. In Section 3, we describe the joint estimation algorithm under tensor structure and list the steps of the proposed algorithm. We perform computational complexity analysis in Section 4. For Section 5, we analyze the simulation performance. Finally, we summarize our efforts in Section 6. The notations used in this paper are described in

Table 1. Notations.

Transpose	${(•)}^{\mathrm{T}}$
Conjugate	${(•)}^{\ast }$
Hermitian transpose	${(•)}^{\mathrm{H}}$
Khatri-Rao product	⊙
Kronecker product	⊗
Hadamard product	⊕
Tensor outer product	∘
Convolution	∗
Identity matrix	I
Statistical expectation	$\mathrm{E}\{•\}$
Orthogonalization	$\mathrm{orth}\left[•\right]$
Tensor contraction along the $q\mathrm{th}$ dimension	${\langle •,•\rangle }_{\left\{q\right\}}$

.

2. Signal Model

2.1. Time-Domain Model

In order to construct the space–time asymmetry, the antenna is designed as in Figure 1 with $M\times M$ RIS elements. The model is used to estimate TD and DOA without designing signal source locations. Multipath wireless propagation channels are usually modeled as complex low-pass equivalent impulse responses. Following the channel estimation, the channel impact response of the RIS element at location $(x,y)$ in the $s\mathrm{th}$ time interval is denoted by

$$h_{x,y}^{\left(s\right)}\left(t\right)=\sum _{k=1}^K{\alpha }_k^{\left(s\right)}{e}^{j{\beta }_k^{\left(s\right)}}\delta \left(t-{\tau }_k\right)\ast \left[\delta \left(t-{\xi }_{k,x,y}\right)\mathrm{\Gamma }\left(t\right)\right],$$

(1)

where K indicates the number of multipaths; ${\alpha }_k^{\left(s\right)}{e}^{j{\beta }_k^{\left(S\right)}}$ is the complex attenuation of the $k\mathrm{th}$ path; ${\alpha }_k^{\left(s\right)}$ is the amplitude; ${\beta }_k^{\left(s\right)}$ is the phase that is consistent with the uniform distribution of the density function $U(0,2\pi )$.

Suppose the RIS element at the origin of the coordinate axis is the reference RIS element, and the propagation delay of the source reaching the reference RIS element through the $k\mathrm{th}$ path is ${\tau }_k$ in the multipath environment. The relative delay of the RIS element at location $(x,y)$ with regard to the reference position is ${\xi }_{k,x,y}$, which is indicated by

$${\xi }_{k,x,y}=\frac{\lambda sin{\theta }_k\left(xcos{\phi }_k+ysin{\phi }_k\right)}{2c},$$

(2)

in which $\lambda$ is the impinging signal wavelength; $\left({\phi }_k,{\theta }_k\right)$ is the direction of the incoming wave; $\mathrm{\Gamma }\left(t\right)$ is the coding coefficient of RIS, and we set its cycle to T with the ODFM symbol period. Furthermore, we have

$$\mathrm{\Gamma }\left(t\right)=\sum _{w=0}^{W-1}{\mathrm{\Gamma }}_{w}g(t-m\psi ),(0\le t\le T),$$

(3)

where the period T is uniformly divided into W; the coding coefficient of the $w\mathrm{th}$ time segment is a constant ${\mathrm{\Gamma }}_w$; $g\left(t\right)$ is a periodic unit pulse signal with a width of $\psi =T/W$, which can be expressed as

$$g(t)=\left\{\begin{array}{cc}& 1,(0\le t\le T)\\ & 0,(0\le t\le T)\end{array}\right.,$$

(4)

which can be expressed in terms of the Fourier series as

$$g\left(t\right)=\sum _{l=0}^{L-1}{c}_l{e}^{j2\pi \Delta flt},$$

(5)

where $\Delta f=1/T$ is the OFDM subcarrier spacing, and L is the quantity of OFDM subcarriers. Moreover, $c_l$ is the Fourier series coefficient, which can be expressed as

$$c_l=\frac{1}{W}\frac{sin\left(\frac{l\pi }{W}\right)}{\frac{l\pi }{W}}{e}^{-j\frac{l\pi }{W}}=\frac{1}{W}Sa\left(\frac{l\pi }{W}\right)e^{-j\frac{l\pi }{W}}.$$

(6)

We take (5) into (3) to obtain

$$\begin{array}{cc}\hfill \mathrm{\Gamma }\left(t\right)& =\sum _{l=0}^L\left[c_l\sum _{w=0}^{W-1}{\mathrm{\Gamma }}_w{e}^{-j{l}_g\frac{2w\pi }{W}}\right]e^{j2\pi \Delta flt}\hfill \\ \hfill & =\sum _{l=0}^L{\alpha }_{\mathrm{\Gamma }}{e}^{j{\beta }_{\mathrm{\Gamma }}}{e}^{j2\pi \Delta flt}\hfill \end{array},$$

(7)

where ${\alpha }_{\mathrm{\Gamma }}{e}^{j{\beta }_{\mathrm{\Gamma }}}$ is the composite attenuation caused by RIS; ${\alpha }_{\mathrm{\Gamma }}$ is the amplitude, and $e^{j{\beta }_{\mathrm{\Gamma }}}$ is the phase, whose effects can be neglected by properly setting W and ${\mathrm{\Gamma }}_w$ to achieve ${\alpha }_{\mathrm{\Gamma }}{e}^{j{\beta }_{\mathrm{\Gamma }}}\alpha e^{j\beta }\approx \alpha e^{j\beta }$.

2.2. Frequency-Domain Model

Performing the Fourier transform to (1), we obtain the Coded-Response of the $l\mathrm{th}$ subcarrier at the $(x,y)$ RIS element as

$$\begin{array}{cc}\hfill H_{l,x,y}^{\left(s\right)}& =\sum _{k=1}^K{\alpha }_k^{\left(s\right)}{e}^{j{\beta }_k^{\left(s\right)}}{e}^{-j2\pi (f_c+\Delta fl){\tau }_k}\hfill \\ \hfill & (e^{-j2\pi (f_c+\Delta fl){\xi }_{k,x,y}}\ast \mathrm{\Gamma }\left(f\right))+n_{l,x,y}^{\left(s\right)}\hfill \\ \hfill & =\sum _{k=1}^K{\alpha }_k^{\left(s\right)}{e}^{j{\beta }_k^{\left(s\right)}}{e}^{-j2\pi (f_c+\Delta fl){\tau }_k}{e}^{-j2\pi f_c{\xi }_{k,x,y}}+n_{l,x,y}^{\left(s\right)}\hfill \end{array},$$

(8)

where $f_c$ is the carrier frequency, and $n_{l,x,y}^{\left(s\right)}$ is additive white Gaussian noise with power ${\sigma }^2$. According to (8), the Coded-Response of $l\mathrm{th}$ subcarrier is given by

$$\begin{array}{cc}\hfill {\mathit{H}}_l^{\left(s\right)}& =\sum _{k=1}^K\left[\begin{array}{c}{H}_{l,0,0}^{\left(s\right)}\\ H_{l,0,1}^{\left(s\right)}\\ ⋮\\ H_{l,M-1,M-1}^{\left(s\right)}\end{array}\right]\hfill \\ \hfill & =\sum _{k=1}^K{A}_{\tau l}\left({\tau }_k\right)\odot {\mathit{A}}_{\xi }\left({\mathit{\xi }}_k\right){\rho }_k^{\left(s\right)}+n_l^{\left(s\right)}\hfill \end{array},$$

(9)

where

$$A_{{\tau }_l}\left({\tau }_k\right)=e^{-j2\pi (f_c+\Delta fl){\tau }_k},$$

(10)

$${\mathit{A}}_{\xi }\left({\mathit{\xi }}_k\right)={\left[\begin{array}{ccc}{e}^{-j2\pi f_c{\xi }_{k,0,0}}& \cdots & e^{-j2\pi f_c{\xi }_{k,M-1,M-1}}\end{array}\right]}^{\mathrm{T}},$$

(11)

$${\rho }_k^{\left(s\right)}={\alpha }_k^{\left(s\right)}{e}^{j{\beta }_k^{\left(s\right)}},$$

(12)

$${\mathit{n}}_l^{\left(s\right)}={\left[\begin{array}{cccc}{n}_{l,0,0}^{\left(s\right)}& n_{l,0,1}^{\left(s\right)}& \cdots & n_{l,M-1,M-1}^{\left(s\right)}\end{array}\right]}^{\mathrm{T}}.$$

(13)

The subcarriers in the OFDM signal can be similar to the RIS arrays based on space–time equivalence. Based on the signal subspace algorithms, a space–frequency Coded-Response matrix can be constructed with the expression

$$\begin{array}{cc}\hfill & {\mathit{H}}^{\left(s\right)}=\sum _{k=1}^K\left[\begin{array}{c}{A}_{\tau 0}\left({\tau }_k\right)\\ A_{\tau 1}\left({\tau }_k\right)\\ ⋮\\ A_{\tau (L-1)}\left({\tau }_k\right)\end{array}\right]\odot {\mathit{A}}_{\xi }\left({\mathit{\xi }}_k\right){\rho }_k^{\left(s\right)}+{\mathit{n}}^{\left(s\right)}\hfill \\ \hfill & =\sum _{k=1}^K{\mathit{A}}_{\tau }\left({\mathit{\tau }}_k\right)\odot {\mathit{A}}_{\xi }\left({\mathit{\xi }}_k\right){\mathit{\rho }}_k^{\left(s\right)}+{\mathit{n}}^{\left(s\right)}\hfill \end{array},$$

(14)

where

$${\mathit{n}}^{\left(s\right)}={\left(\begin{array}{cccc}{{\mathit{n}}^{\left(s\right)}}_0^{\mathrm{T}}& {{\mathit{n}}^{\left(s\right)}}_1^{\mathrm{T}}& \cdots & {{\mathit{n}}^{\left(s\right)}}_{L-1}^{\mathrm{T}}\end{array}\right)}^{\mathrm{T}}.$$

(15)

In addition, if there are an overall S time intervals, then (14) is represented as

$$\mathit{H}={\mathit{A}}_{\tau }\left(\mathit{\tau }\right)\odot {\mathit{A}}_{\xi }\left(\mathit{\xi }\right)\mathit{\rho }+\mathit{n},$$

(16)

where

$$\mathit{H}=\left[\begin{array}{cccc}{\mathit{H}}^{\left(1\right)}& {\mathit{H}}^{\left(2\right)}& \cdots & {\mathit{H}}^{\left(S\right)}\end{array}\right],$$

(17)

$${\mathit{A}}_{\tau }\left(\mathit{\tau }\right)=\left[\begin{array}{cccc}{\mathit{A}}_{\tau }\left({\tau }_1\right)& {\mathit{A}}_{\tau }\left({\tau }_2\right)& \cdots & {\mathit{A}}_{\tau }\left({\tau }_K\right)\end{array}\right],$$

(18)

$${\mathit{A}}_{\xi }\left(\mathit{\xi }\right)=\left[\begin{array}{cccc}{\mathit{A}}_{\xi }\left({\mathit{\xi }}_1\right)& {\mathit{A}}_{\xi }\left({\mathit{\xi }}_2\right)& \cdots & {\mathit{A}}_{\xi }\left({\mathit{\xi }}_K\right)\end{array}\right],$$

(19)

$$\mathit{\rho }={\left[\begin{array}{cccc}{\mathit{\rho }}_1& {\mathit{\rho }}_1& \cdots & {\mathit{\rho }}_K\end{array}\right]}^{\mathrm{T}},$$

(20)

$${\mathit{\rho }}_k=\left[\begin{array}{cccc}{\rho }_k^{\left(1\right)}& {\rho }_k^{\left(2\right)}& \cdots & {\rho }_k^{\left(S\right)}\end{array}\right],$$

(21)

$$\mathit{n}=\left[\begin{array}{cccc}{\mathit{n}}^{\left(1\right)}& {\mathit{n}}^{\left(2\right)}& \cdots & {\mathit{n}}^{\left(S\right)}\end{array}\right].$$

(22)

The time-domain encoding of RIS transforms the combination of ${\mathit{A}}_{\tau }\left(\mathit{\tau }\right)$ and ${\mathit{A}}_{\xi }\left(\mathit{\xi }\right)$ from Hadmard to Khatri–Rao, providing the possibility of solving TD and DOA separately. The coded signal model converts the time-domain structure into a frequency-domain structure so that it is still suitable for the time-domain method. Additionally, there are two clear benefits to building a Coded-Response for joint TD and DOA estimation. On the one hand, the simultaneous acquisition of TD and DOA parameters enables single-station localization to reduce the localization system overhead. On the other hand, the proposed joint estimation approach leads to an increase in L times aperture and $M^2$ times bandwidth compared to reality, which is better than the individual estimation performance.

3. The Joint TD and DOA Estimation

3.1. Tensor Approach

In this section, we represent the Coded-Response (16) into a tensor. The tensorial form can reveal signal model structure clearly, which enables us to realize signal and noise subspaces more accurately. In addition, the proposed method improves the robustness in a coherent multipath environment with a tensor form of

$$\mathcal{H}={\mathit{A}}_{\tau }\left(\mathit{\tau }\right)\circ {\mathit{A}}_{\xi }\left(\mathit{\xi }\right)\circ {\mathit{\rho }}^{\mathrm{T}}+\mathcal{N},$$

(23)

where $\mathcal{H}\in {\mathbb{C}}^{M^2\times L\times S}$, and $\mathcal{N}$ is additive white Gaussian noise in the same dimension as $\mathcal{H}$ with power of ${\sigma }^2$. In (23), the tensor covariance matrix $\mathcal{R}{}_{\mathcal{H}}\in {\mathbb{C}}^{M^2\times L\times M^2\times L}$ can be obtained as

$$\begin{array}{cc}\hfill & \mathcal{R}{}_{\mathcal{H}}=E\left[<\mathcal{H},{\mathcal{H}}^{\ast }{>}_{\left\{3\right\}}\right]\hfill \\ \hfill & =\sum _{k=1}^K{\beta }_k^2{\mathit{A}}_{\tau }\left({\tau }_k\right)\circ {\mathit{A}}_{\xi }\left({\mathit{\xi }}_k\right)\circ {\mathit{A}}_{\tau }^{\ast }\left({\tau }_k\right)\circ {\mathit{A}}_{\xi }^{\ast }\left({\mathit{\xi }}_k\right)+{\mathcal{N}}_{\mathcal{R}}\hfill \end{array}\prime$$

(24)

where ${\beta }_k^2=\mathrm{E}\left[{\mathit{\rho }}_k^{\mathrm{T}}{\left({\mathit{\rho }}_k^{\mathrm{T}}\right)}^{\ast }\right]$ is the complex decay of the $k\mathrm{th}$ path, and ${\mathcal{N}}_{\mathcal{R}}=\mathrm{E}\left[{〈\mathcal{N},{\mathcal{N}}^{\ast }〉}_{\left\{3\right\}}\right]$ is the noise term. In practice, the tensor covariance matrix $\mathcal{R}{}_{\mathcal{H}}$ can be estimated as

$${\widehat{\mathcal{R}}}_{\mathcal{H}}=\frac{1}{S}{〈\mathcal{H},{\mathcal{H}}^{\ast }〉}_{\left\{3\right\}}.$$

(25)

The CANDECOMP/PARAFAC decomposition (CPD) [26] is a common method for splitting a high-dimensional tensor. We perform CPD on the four-dimensional tensor ${\widehat{\mathcal{R}}}_{\mathcal{H}}$ to obtain a sum ${\mathcal{R}}_{\mathcal{C}\mathcal{P}}\in {\mathbb{C}}^{M^2\times K,L\times K,M^2\times K,L\times K}$ of component rank-one tensors as

$${\mathcal{R}}_{CP}=\left[{\mathcal{T}}_{CP}\right],\left[{\mathcal{A}}_{CP}\right],\left[{\mathcal{T}}_{CP}^{\ast }\right],\left[{\mathcal{A}}_{CP}^{\ast }\right].$$

(26)

Due to the milder uniqueness condition of tensor decomposition and the asymmetry of the constituent factor vectors in (23), there is

$$\begin{array}{cc}\hfill & \mathrm{rank}\left({\mathcal{T}}_{CP}\right)+\mathrm{rank}\left({\mathcal{A}}_{CP}\right)+\mathrm{rank}\left(\mathcal{T}{}_{CP}{}^{\ast }\right)+\mathrm{rank}\left({\mathcal{A}}_{CP}{}^{\ast }\right)\ge 2K+N-1\hfill \\ \hfill & \Rightarrow K+K+K+K\ge 2K+N-1\hfill \\ \hfill & \Rightarrow 2K\ge 3\hfill \end{array}.$$

(27)

N is the amount of matrices obtained by CPD or the dimension of the decomposed matrix, which is four in this algorithm. Therefore, when the number of sources is not unique, the proposed algorithm satisfies the uniqueness of CPD [26]. Further, we extract the signal subspace ${\mathcal{T}}_{CP}$ from ${\mathcal{R}}_{CP}$ for TD estimation and the signal subspace ${\mathcal{A}}_{CP}$ for DOA estimation.

3.2. TD Estimation

Since the phases of TD and DOA are asymmetric in frequency, it allows TD to be estimated separately while suppressing the coherent signal. To demonstrate this process, we will fold the $M^{2}L\times 1$ received vector (14) into an $L\times M^2$ matrix as follows

$${\mathit{H}}_{Tfold}^{\left(s\right)}(p,q)={\mathit{H}}^{\left(s\right)}(M^2(p-1)+q),$$

(28)

where ${\mathit{H}}_{Tfold}^{\left(s\right)}\in {\mathbb{C}}^{L\times M^2}$ is a new observation matrix of $p=1,2,\dots ,L$ and $q=1,2,\dots ,M^2$. Thus, TD and DOA can be decoupled in phase as

$${\mathit{H}}_{Tfold}^{\left(s\right)}=\sum _{k=1}^K{\mathit{A}}_{\tau }\left({\tau }_k\right){\mathit{A}}_{\xi }^{\mathrm{T}}\left({\mathit{\xi }}_k\right){\rho }_k^{\left(s\right)}+{\mathit{n}}_{{Tfold}^{\prime }}^{\left(s\right)}$$

(29)

where ${\mathit{n}}_{Tfold}^{\left(s\right)}$ denotes the noise term and still follows the white Gaussian distribution. Further, ${\mathit{H}}_{Tfold}^{\left(s\right)}$ can be considered as a sample collected by the RIS at $s\mathrm{th}$ time interval. Here, we define

$${\mathit{F}}^{\left(s\right)}={\left[\begin{array}{cccc}{\mathit{A}}_{\xi }^{\mathrm{T}}\left({\mathit{\xi }}_1\right){\rho }_1^{\left(s\right)}& {\mathit{A}}_{\xi }^{\mathrm{T}}\left({\mathit{\xi }}_2\right){\rho }_2^{\left(s\right)}& \cdots & {\mathit{A}}_{\xi }^{\mathrm{T}}\left({\mathit{\xi }}_K\right){\rho }_K^{\left(s\right)}\end{array}\right]}^{\mathrm{T}}.$$

(30)

Then the covariance matrix corresponding to (29) is

$$\begin{array}{cc}\hfill {\mathit{R}}_{Tfold}& =\frac{1}{S}\sum _{s=1}^S{\mathit{H}}_{Tfold}^{\left(s\right)}{\left({\mathit{H}}_{Tfold}^{\left(s\right)}\right)}^{\mathrm{H}}\hfill \\ \hfill & ={\mathit{A}}_{\tau }\left(\mathit{\tau }\right){\mathit{R}}_F\mathit{A}{\left(\mathit{\tau }\right)}^{\mathrm{H}}+{\sigma }^2{\mathit{I}}_L\hfill \end{array},$$

(31)

where

$$\begin{array}{cc}\hfill & {\mathit{R}}_F=\frac{1}{S}\sum _{s=1}^S{\mathit{F}}^{\left(s\right)}{\left({\mathit{F}}^{\left(s\right)}\right)}^{\mathrm{H}}\hfill \\ \hfill & ={\mathit{R}}_{\xi }\oplus {\mathit{R}}_{\rho }\hfill \end{array},$$

(32)

where

$${\mathit{R}}_{\xi }=\left[\begin{array}{ccc}{\mathit{A}}_{\xi }^{\mathrm{T}}\left({\mathit{\xi }}_1\right){\mathit{A}}_{\xi }^{\ast }\left({\mathit{\xi }}_1\right)& \cdots & {\mathit{A}}_{\xi }^{\mathrm{T}}\left({\mathit{\xi }}_1\right){\mathit{A}}_{\xi }^{\ast }\left({\mathit{\xi }}_K\right)\\ ⋮& \ddots & ⋮\\ {\mathit{A}}_{\xi }^{\mathrm{T}}\left({\mathit{\xi }}_K\right){\mathit{A}}_{\xi }^{\ast }\left({\mathit{\xi }}_1\right)& \cdots & {\mathit{A}}_{\xi }^{\mathrm{T}}\left({\mathit{\xi }}_K\right){\mathit{A}}_{\xi }^{\ast }\left({\mathit{\xi }}_K\right)\end{array}\right],$$

(33)

$${\mathit{R}}_{\rho }=\frac{1}{S}\sum _{s=1}^S{\mathit{\rho }}^{\left(s\right)}{\left({\mathit{\rho }}^{\left(s\right)}\right)}^{\mathrm{H}}.$$

(34)

Note that ${\mathit{A}}_{\xi }^{\mathrm{T}}\left({\mathit{\xi }}_k\right){\mathit{A}}_{\xi }^{\ast }\left({\mathit{\xi }}_i\right)$, $k,i=1,2,\cdots ,K$ can only achieve maximum value if $k=i$. Since the relative delays of the signals are different from each other, ${\mathit{R}}_{\xi }$ is a diagonally dominant matrix which reduces the singularity of the signal correlation matrix when the signals are correlated. Combining with (32), it can be seen that in the coherent environment even if $\mathrm{rank}\left({\mathit{R}}_{\rho }\right)<K$ there is still $\mathrm{rank}\left({\mathit{R}}_{\xi }\right)<K$ enabling $\mathrm{rank}\left({\mathit{R}}_F\right)<K$. Therefore, if the proposed method can obtain the signal or noise subspace of ${\mathit{A}}_{\tau }\left(\mathit{\tau }\right)$ alone for TD estimation, the method is adaptive to coherent signals.

Then we obtain the corresponding signal subspace ${\mathcal{T}}_{CP}$ of ${\mathit{A}}_{\tau }\left(\mathit{\tau }\right)$ from (26) and construct the corresponding noise subspace ${\mathit{U}}_{nt}{\mathit{U}}_{nt}^{\mathrm{H}}$ as follows

$${\mathit{U}}_{nt}{\mathit{U}}_{nt}^{\mathrm{H}}=\mathit{I}-\mathrm{orth}\left({\mathcal{T}}_{CP}\right)\mathrm{orth}{\left({\mathcal{T}}_{CP}\right)}^{\mathrm{H}}.$$

(35)

Then the spatial spectrum expression of TD estimation is

$$P\left(\tau \right)=\frac{1}{\mathcal{T}{\left(\tau \right)}^{\mathrm{H}}\left({\mathit{U}}_{nt}{\mathit{U}}_{nt}^{\mathrm{H}}\right)\mathcal{T}\left(\tau \right)}.$$

(36)

We can use the one-dimensional (1D) MUSIC algorithm to solve the value of $\widehat{\tau }$, which has a higher estimation accuracy and significantly reduces the complexity compared to the three-dimensional (3D) MUSIC algorithm. Further, in the proposed tensor model, the separate TD estimates are not affected by the coherent multipath environment.

3.3. DOA Estimation

As the phases of TD and DOA are asymmetric in frequency, it allows DOA to be estimated separately while suppressing the coherent signal. To demonstrate this process, we will fold the $M^{2}L\times 1$ received vector (14) into an $M^2\times L$ matrix as follows

$${\mathit{H}}_{Afold}^{\left(s\right)}(p,q)={\mathit{H}}^{\left(s\right)}(p+M^2(q-1)),$$

(37)

where ${\mathit{H}}_{Afold}^{\left(s\right)}\in {\mathbb{C}}^{M^2\times L}$ is a new observation matrix of $p=1,2,\dots ,M^2$ and $q=1,2,\dots ,L$. Thus, DOA and TD can be decoupled in phase as

$${\mathit{H}}_{Afold}^{\left(s\right)}=\sum _{k=1}^K{\mathit{A}}_{\xi }\left({\mathit{\xi }}_k\right){\mathit{A}}_{\tau }^{\mathrm{T}}\left({\tau }_k\right){\rho }_k^{\left(s\right)}+{\mathit{n}}_{{Afold}^{\prime }}^{\left(s\right)}$$

(38)

where ${\mathit{n}}_{Afold}^{\left(s\right)}\in {\mathbb{C}}^{M^2\times L}$ stands for the complex white Gaussian noise. The remaining proof process is similar to (31)–(34), and the proposed method is adaptive to coherent signals if the signal or noise subspace of ${\mathit{A}}_{\xi }\left(\mathit{\xi }\right)$ can be obtained separately for DOA estimation. Then we obtain the corresponding signal subspace ${\mathcal{A}}_{CP}$ of ${\mathit{A}}_{\xi }\left(\mathit{\xi }\right)$ from (26) and construct the corresponding noise subspace ${\mathit{U}}_{na}{\mathit{U}}_{na}^{\mathrm{H}}$ as follows

$${\mathit{U}}_{na}{\mathit{U}}_{na}^{\mathrm{H}}=\mathit{I}-\mathrm{orth}\left({\mathcal{A}}_{CP}\right)\mathrm{orth}{\left({\mathcal{A}}_{CP}\right)}^{\mathrm{H}}.$$

(39)

Then the spatial spectrum expression of the DOA estimate can be obtained as

$$P\left(\theta ,\phi \right)=\frac{1}{\mathcal{A}{(\theta ,\phi )}^{\mathrm{H}}\left({\mathit{U}}_{na}{\mathit{U}}_{na}^{\mathrm{H}}\right)\mathcal{A}(\theta ,\phi )}.$$

(40)

We can use the two-dimensional (2D) MUSIC algorithm to solve the values of $(\widehat{\theta },\widehat{\phi })$, which has a higher estimation accuracy and greatly reduces the complexity compared to the 3D MUSIC algorithm. Further, in the proposed tensor model, the separate TD estimates are not affected by the coherent multipath environment.

3.4. Algorithm Steps

The major steps of the proposed algorithm are presented in Algorithm 1.

Algorithm 1 Algorithm steps.
step1:	The Coded-Response $\mathcal{H}$ in tensor form is constructed according to (23).
step2:	The tensor covariance matrix $\stackrel{⌢}{\mathcal{R}}{}_{\mathcal{H}}$ is constructed according to (25).
step3:	Perform a CPD of ${\widehat{\mathcal{R}}}_{\mathcal{H}}$, which solves the signal subspaces ${\mathcal{T}}_{CP}$ and ${\mathcal{V}}_{CP}$, and then the corresponding noise subspaces and are obtained by (35) and (39), respectively. ${\mathit{U}}_{nt}{\mathit{U}}_{nt}^{\mathrm{H}}$ and ${\mathit{U}}_{na}{\mathit{U}}_{na}^{\mathrm{H}}$.
step4:	Proceed to conduct a 1D spectral peak search for ${\mathit{U}}_{nt}{\mathit{U}}_{nt}^{\mathrm{H}}$ to solve the $\widehat{\tau }$ by (36), and a 2D spectral peak search for ${\mathit{U}}_{na}{\mathit{U}}_{na}^{\mathrm{H}}$ to solve the $(\widehat{\theta },\widehat{\phi })$ by (40).

4. Algorithms Complexity

In this part, we analyze the complexity of the proposed algorithm under the Coded-Response model and compare it with the 3D spectral peak search method (3D-MUSIC). To analyze the accuracy in the next section, we also compare the complexity of the algorithm (joint-smooth) employing smoothing preprocessing [22], which reduces the dimensions of the joint estimation model.

The related complexity is mainly split into vector-based covariance matrix computation, tensor-based covariance matrix computation, eigenvalue decomposition, and CPD and 1D spectrum peak search, which are $\mathcal{O}\left(S{M}^4{L}^2\right)$, $\mathcal{O}\left(S{M}^2{L}^2\right)$, $\mathcal{O}\left(M^3\right)$, $\mathcal{O}(2^{N}K{M}^4{L}^2+N{K}^3)$, and $\mathcal{O}\left(M\right(M-K\left)G\right)$, separately, where G represents the quantity of spectral points in the 1D search. Hence, the proposed complexity is $\mathcal{O}\left(M^2{L}^2\right)(S+2^{N}K{M}^2)+N{K}^3+M^2(M^2-K)G_{\phi }{G}_{\theta }+L(L-K)G_{\tau }$. The cost of the 3D-MUSIC is $\mathcal{O}((S+M^{2}L)M^4{L}^2+M^{2}L(M^{2}L-K)G_{\phi }{G}_{\theta }{G}_{\tau })$, as well as the smoothing preprocessing method (joint-smooth) is $\mathcal{O}(S{M}^4{L}^2+M^6+L^3+M^2(M^2-K)G_{\phi }{G}_{\theta }+L(L-K)G_{\tau })$, in which $G_{\phi }$, $G_{\theta }$, and $G_{\tau }$ denote the quantity of searches for azimuth, elevation, and propagation delay, separately. To make the comparison clear,

Table 2. Comparison of complexity.

Algorithm	Complexity
Proposed	$\mathcal{O}\left(M^2{L}^2\right)(S+2^{N}K{M}^2)+N{K}^3+L(L-K)G_{\tau }+M^2(M^2-K)G_{\phi }{G}_{\theta }$
3D-Music	$\mathcal{O}((S+M^{2}L)M^4{L}^2+M^{2}L(M^{2}L-K)G_{\phi }{G}_{\theta }{G}_{\tau })$
Joint-Smooth	$\mathcal{O}(S{M}^4{L}^2+M^6+L^3+M^2(M^2-K)G_{\phi }{G}_{\theta }+L(L-K)G_{\tau })$

summarizes the complexity of all methods. Additionally, we have compared the complexity of the algorithms in terms of time intervals $\left(S\right)$, the quantity of RIS elements $\left(M\right)$, and the searching step ($\Delta \tau$, $\Delta \phi =\Delta \theta$, in which $G_{\phi }=360/\Delta \phi$, $G_{\theta }=90/\Delta \theta$) in Figure 2a–d, separately.

From Figure 2, the 3D-MUSIC complexity is substantial because of the large quantity of spectral points, specifically in the small spectral steps. The high complexity of the 3D-MUSIC algorithm brings a very high estimation accuracy, but it cannot be used for coherent signals. In contrast, the proposed and joint-smooth methods are used to reduce the spectral peak search dimension by tensor CPD and full-dimensional smoothing preprocessing, which leads to a significant reduction in complexity. Thus, both algorithms have similar complexity, and both are applicable in coherent environments. Due to the large-scale smoothing preprocessing that leads to more aperture loss, the estimation performance of the joint-smooth method is lower than that of the proposed method.

5. Simulation Results

Here, we present simulation results which focus on the analysis of the proposed algorithm with joint-smooth. In addition, the proposed algorithm is compared with 3D-MUSIC in the independent multipath environment.

5.1. Performance at Low SNR

Three cases are designed in this paper. First, suppose that the coherent signals are two, the corresponding delays are $3.5\mathrm{ns}$ and $23.5\mathrm{ns}$, related azimuth angles are $-{20}^{\circ }$ and ${30}^{\circ }$, and corresponding elevation angles are ${20}^{\circ }$ and ${45}^{\circ }$, respectively. Secondly, we append an independent signal with $23.5\mathrm{ns}$ time delay, ${30}^{\circ }$ azimuth angle, and ${45}^{\circ }$ elevation angle. Thirdly, we assume all three signals are independent. Moreover, $Q=100$ is set and the distribution of TD and DOA under $S=500$, and $\mathrm{SNR}=-5\mathrm{dB}$ is determined, as illustrated in Figure 3. Figure 3 shows that this algorithm has the ability to estimate TD and DOA values in various situations. In addition, the proposed algorithm has good robustness and performance under low SNR.

5.2. Performance versus SNR

This part analyzes the capabilities of the proposed algorithm, Joint-smooth, 3D-MUSIC, TD estimation [27], DOA estimation [28], and Cramer–Rao bound (Appendix A) in the coherent and independent multipath environments. Suppose the multipath signals are three with the same parameters as simulated in Section 5.1. We set $Q=100$, $S=500$, and the spectral steps of $\Delta \tau =0.001\mathrm{ns}$ and $\Delta \phi =\Delta \theta ={0.05}^{\circ }$. Further, the RMSE performance versus SNR with ranges of $-15\mathrm{dB}$ to $20\mathrm{dB}$ in $5\mathrm{dB}$ intervals is shown in Figure 4 and Figure 5.

As shown in Figure 4, the joint estimation algorithms (joint-smooth and proposed) in the coherent multipath environment work significantly better than the performance of the independent estimation algorithms (DOA-Smooth and Smooth-TD) due to the aperture and bandwidth enhancement brought by the joint approach. Combined with the complexity comparison of Section 4, it can be concluded that the proposed has a higher estimation accuracy than the joint-smooth, while the complexity is very close for $\mathcal{O}\left(2.88\times {10}^9\right)$ and $\mathcal{O}\left(3.32\times {10}^9\right)$. The the proposed avoids the aperture loss in processing coherent signals by the time-domain coding function of the RIS.

As illustrated in Figure 5, the RMSE of 3D-MUSIC is closer to CRB than the proposed in the independent multipath environment because the CPD process of the proposed loses part of the extended structure of the subspace. However, the complexity of 3D-MUSIC is $\mathcal{O}\left(3.39\times {10}^{17}\right)$, which is much higher than the proposed and difficult to accept in practical scenarios.

5.3. Performance versus Time Intervals

To highlight the impact of time intervals on RMSE, we set $\mathrm{SNR}=15\mathrm{dB}$, the quantity of time intervals varies at the range of $S=\{20,50,100,200,500,1000,2000,5000\}$, and the rest of the simulation parameters are identical to those simulated in Section 5.2.

As shown in Figure 6, when processing coherent signals, the joint-smooth and the proposed have significantly better performance than DOA-Smooth and Smooth-TD due to the aperture and bandwidth improvement brought by the joint approach. Further, the proposed avoids the aperture loss of the current smoothing method when processing coherent signals by the time-domain coding function of RIS. Therefore, its estimation performance is better than the joint-smooth.

As shown in Figure 7, in the independent multipath environment, 3D-MUSIC has a lower RMSE than the proposed, which is due to the fact that the CPD process of the proposed loses part of the subspace expansion structure. Other conclusions are the same as in Figure 6.

In summary, the RMSE of the proposed algorithm reduces as the time interval increases, although the reduction gradually achieves a plateau. Additionally, more algorithms are adapted to independent environments, and the estimated performance is typically better.

5.4. Performance versus RIS Elements

For highlighting the effect of RIS elements on RMSE, we set $\mathrm{SNR}=15\mathrm{dB}$, the number of RIS elements varies in the range of $M=\{3,4,5,6,7,8,9,10\}$, and the rest of the simulation parameters are the same as those simulated in Section 5.2.

In Figure 8, we simulate the performance of the joint estimation algorithm in a coherent environment. As illustrated in (a) and (b), the performance of the joint estimation DOA gets better as the number of RIS elements increases. However, the trend flattens out, which is because the algorithm performance is mainly limited by the SNR when there are enough RIS elements.

As shown in (c), the performance of the joint estimation TD also gets better as the number of RIS elements increases. However, it is much less affected due to the fact that TD estimation is mainly limited by the signal bandwidth. Even though the increase of RIS elements in the joint estimation model expands the signal bandwidth, its effect is not infinitely increasing.

6. Conclusions

For the joint estimation for TD and DOA of coherent signals on a single station system, we propose an algorithm based on RIS and OFDM techniques in this paper. The proposed algorithmic model integrates the RIS array response and the channel frequency response of OFDM subcarriers in tensor form to extend the bandwidth and aperture and achieve better estimation. The appropriate combination of the RIS time-domain coding function and the tensor structure allows the algorithm to estimate coherent signals with reduced complexity. Compared with the current algorithm, the proposed algorithm has better estimation performance in processing coherent signals without the aperture loss of the smoothing process and can rely on the advantage of joint estimation for single station localization systems. In this paper, we begin by presenting the signal model and the algorithmic procedure, which has detailed derivations and proofs. Then, to reflect the fast estimation of the proposed algorithm, we compare the related algorithms in terms of complexity. Finally, the validity and robustness of the algorithm are verified by simulation experiments. In comparison with current algorithms, the proposed algorithm provides more effective parameter estimation with lower complexity in coherent or mixed multipath environments.