Coherent Signal DOA Estimation Method Based on Space–Time–Coding Metasurface

Abstract

A novel method for the direction of arrival (DOA) estimation of coherent signals under a space–time–coding metasurface (STCM) is proposed in this paper. Noticeably, the STCM can replace multi-channel arrays with a single channel, which can be utilized to modulate incident electromagnetic waves and generate harmonics. However, coherent signals are overlapping in the frequency spectrum and cannot achieve DOA estimation through subspace methods. Therefore, the proposed method transforms the angle information in the time domain into amplitude and phase information at harmonics in the frequency domain by modulating incident coherent signals using the STCM and performing a fast Fourier transform (FFT) on these signals. Based on the harmonics in the frequency spectrum of the coherent signals, appropriate harmonics are selected. Finally, the ℓ₁ norm singular value decomposition (ℓ₁-SVD) algorithm is utilized for achieving high-precision DOA estimation. Simulation experiments are conducted to show the performance of the proposed method under the condition of different incident angles, harmonic numbers, signal-to-noise ratios (SNRs), etc. Compared to the traditional algorithms, the performance of the proposed algorithm can achieve more accurate DOA estimation under a low SNR.

1. Introduction

Over the past few decades, remarkable progress has been made in radar and wireless sensing technologies, which can be adopted in fields such as radar detection, remote sensing, navigation, and indoor positioning [1,2,3,4]. Direction of arrival (DOA) estimation is a pivotal technology in the domain of perception and localization techniques, which is essential for mapping wireless electromagnetic environments [5,6,7]. Traditionally, the study of DOA estimation problems has primarily involved array signal processing techniques. A robust and mature algorithmic framework has been established in this area, which can be broadly categorized into several classes. Beamforming methods have been extensively studied [8], including the conventional beamforming (CBF) method and the minimum variance distortionless response (MVDR) algorithm. Subspace-based methods have also been developed, such as the multiple signal classification (MUSIC) algorithm [9] and the estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm [10]. Fitting methods have been developed, such as the maximum likelihood (ML) method [11]. Compressed sensing (CS) algorithms have been explored as well, which are exemplified by the orthogonal matching pursuit (OMP) algorithm [12]. Currently, high-resolution methods based on subspace algorithms are the most commonly used. However, these methods often require the use of multiple sensors and receiving channels to capture signals, leading to complex computations and expensive hardware requirements for achieving high-precision DOA estimation.

In recent years, metasurfaces have emerged as a novel class of electromagnetic materials distinguished by their cost savings, low power requirements, and simple deployment [13,14,15,16]. These materials have been utilized to manipulate the scattering field through space–time–coding metasurface (STCM) and subsequent reception via a single channel, achieving performance comparable to that of multi-channel arrays while significantly reducing hardware complexity [17,18]. Consequently, DOA estimation methods based on metasurfaces have attracted considerable attention. In [19], a set of random radiation patterns is generated using a dynamic metasurface, and the OMP algorithm is employed to recover DOA information. In [20], an off-grid DOA estimation approach utilizing the atomic norm minimization of a coded metasurface is proposed. In [21], the orthogonality of the sensing matrix is enhanced by utilizing random dual beams produced by a digital programmable metasurface, thereby increasing the precision of DOA estimation. In [22], a method based on sparse Bayesian learning for two-dimensional off-grid DOA estimation using a metasurface is introduced. Beyond leveraging the spatial modulation capabilities of metasurfaces for DOA estimation, these materials have also demonstrated excellent modulation capabilities in the time and frequency domains [23,24,25]. Numerous studies have exploited these capabilities to address DOA estimation challenges. An STCM can modulate the amplitude and phase information of electromagnetic waves according to the encoding configuration strategy. This modulation results in the production of a series of harmonics that contain the incident angle information. Subsequently, the spectral harmonic properties of the STCM can be analyzed to achieve DOA estimation [26,27]. Several studies have shown that DOA estimation can be effectively achieved by analyzing the harmonic properties generated by the STCM [28,29]. In [30,31], artificial neural network algorithms are introduced to metasurface DOA estimation problems, establishing a mapping relationship between harmonic amplitudes and source signal directions to solve for DOA. However, in the condition of coherent signals, the frequency spectrum is overlapping and cannot achieve DOA estimation through the above-mentioned methods.

In practical applications, antenna arrays may receive signals that become mutually coherent due to factors such as multipath reflections or intentional interference [32]. Current metasurface-based methods for the DOA estimation of coherent signals exhibit certain limitations. References [33,34] employ spatial smoothing techniques to address the direction finding of coherent sources, which suffer from limitations in resolution and accuracy that need to be enhanced. This paper introduces a novel approach that utilizes the STCM for modulating coherent signals, which is followed by the application of the ℓ₁ norm singular value decomposition (ℓ₁-SVD) algorithm for DOA estimation. In detail, the proposed method transforms the angle information in the time domain into amplitude and phase information at harmonics in the frequency domain by modulating incident coherent signals using the STCM and performing a fast Fourier transform (FFT) on these signals. Based on the harmonics in the frequency spectrum of the coherent signals, appropriate harmonics are selected. Finally, the ℓ₁-SVD algorithm processes the harmonics in the frequency spectrum, constructs an observation vector, and reconstructs the sparse vector via ℓ₁ minimization. The final step involves solving a second-order cone program to determine the DOA of coherent signals. This method offers several advantages: the use of the STCM enables single-channel reception, which reduces system costs; the modulation capability of the STCM captures rich spectral information, enhancing spectral efficiency; and compared to existing algorithms, it improves the accuracy and resolution of DOA estimation for coherent signals.

The content of this paper is organized as follows: Section 2 offers a comprehensive explanation of the construction of a coherent signal reception model based on the STCM, elucidating the theoretical framework and mathematical formulation underlying the proposed approach. Section 3 delineates the application of the algorithm to the pre-processed signals obtained from the previous step, explicating the methodology and theoretical underpinnings for resolving the DOA estimation of coherent signals. Section 4 presents the simulation experiments conducted to validate the proposed method, encompassing a comprehensive description and in-depth analysis of the experimental setup, parameters, and results. Section 5 concludes the paper.

2. Coherent Signal Reception Model Based on STCM

As shown in Figure 1, a coherent signal reception model for the STCM is established. We consider the scenario where multiple signals of the same frequency arrive from various angles.

When coherent signals are illuminated by the STCM at incident angles ${\theta }_k$ and are subsequently reflected, they are captured by a horn receiving antenna. For one-dimensional angle estimation, each column of the STCM-reflective elements can be regarded as a subarray, collectively forming the entire STCM array. Consequently, the STCM model can be equivalently represented as a metasurface linear array model. Here, D represents the distance between consecutive elements in the linear array. h represents the distance from the center of the STCM to the horn receiving antenna. n indicates the index of the linear array element where $n\in \{1,2,\cdots ,N\}$. N is the total number of elements in the linear array. The index of the coherent signal is denoted by k, where $k\in \{1,2,\cdots ,K\}$. K represents the total number of coherent signals. The path differences of each subarray are represented as $d_{n,k}$. The path differences from each subarray to the horn receiving antenna are represented as $g_n$.

$$d_{n,k}=-(n-1)Dsin{\theta }_k$$

(1)

$$g_n=\sqrt{{\left((n-1)D\right)}^2+h^2}$$

(2)

The model of the received signals can be mathematically formulated as

$$s\left(t\right)=\sum _{n=1}^N\sum _{k=1}^K{A}_0{U}_n\left(t\right)e^{j[{\omega }_0(d_{n,k}+g_n)+2\pi F_{c}t+{\varphi }_0]}$$

(3)

where $F_c$ represents the frequency of the received signals. $A_0$ and ${\varphi }_0$ denote the initial amplitude and phase of the received signals, respectively. ${\omega }_0$ is the wave number. t is the time variable of the received signals. The STCM sequence for the metasurface linear array model is denoted as $U_n\left(t\right)$, which is expressed as

$$U_n\left(t\right)=\left\{\begin{array}{cc}1,\hfill & m{T}_p+{\tau }_{n,ON}\le t<m{T}_p+{\tau }_{n,OFF}\hfill \\ -1,\hfill & m{T}_p\le t<m{T}_p+{\tau }_{n,ON}\hfill \\ & orm{T}_p+{\tau }_{n,OFF}\le t<(m+1)T_p\hfill \end{array}\right.m\in M$$

(4)

where $T_p$ represents the modulation period, and ${\tau }_{n,ON}$ and ${\tau }_{n,OFF}$ denote the moments when the coding is set to “1” at the beginning and end of each modulation period, respectively. “1” indicates that the phase delay in the subarrays is set to $0^{\circ }$, while “−1” signifies that the phase delay is set to ${180}^{\circ }$. m indicates the number of modulation periods. M represents the total number of modulation periods. Figure 2 illustrates the STCM sequence diagram for $U_n\left(t\right)$.

The Fourier series expansion of the space–time–coding modulation sequence for the metasurface linear array model is given by

$$U_n\left(t\right)=\sum _{q=-\infty }^{\infty }{a}_{n,q}{e}^{j2\pi q{F}_{p}t}$$

(5)

where $F_p$ represents the modulation frequency, which is the reciprocal of the modulation period $T_p$. $a_{n,q}$ denotes the qth-order harmonic coefficient generated by the nth linear array element:

$$\begin{array}{cc}\hfill a_{n,q}& =\left\{\begin{array}{c}{\displaystyle \frac{{\tau }_{n,OFF}-{\tau }_{n,ON}}{T_p}},q=0\hfill \\ \\ {\displaystyle \frac{sin\left[\pi q{F}_p({\tau }_{n,OFF}-{\tau }_{n,ON})\right]}{\pi q}}{e}^{-j\pi q{F}_p({\tau }_{n,OFF}+{\tau }_{n,ON})},q\ne 0\hfill \end{array}\right.\hfill \end{array}$$

(6)

By substituting Equation (6) into Equation (3), we derive the coherent signal reception model for the STCM-based DOA estimation system. The resulting model can be expressed as follows:

$$s\left(t\right)=\sum _{q=-\infty }^{\infty }\sum _{n=1}^N\sum _{k=1}^K{A}_0{a}_{n,q}{e}^{j[{\omega }_0(d_{n,k}+g_n)+2\pi (F_c+q{F}_p)t+{\varphi }_0]}$$

(7)

3. Proposed Algorithm

After modulation by the STCM, the scattered field of coherent signals is received by a horn antenna and then connected to a spectrum analyzer for signal acquisition and sampling. The frequency-domain signal obtained on the spectrum analyzer exhibits the following characteristics: a fundamental spectral line with frequency $F_c$ and harmonic spectral components located at integer multiples of the modulation frequency $F_p$ on both sides of the fundamental spectral line. The amplitudes of the fundamental and harmonic spectral components are represented by the fundamental component coefficient ${\gamma }_0$ and harmonic component coefficients ${\gamma }_q$, respectively. When $q=0$, it denotes the fundamental component coefficient. Combining these with the mathematical relationships in the coherent signal reception model (7), ${\gamma }_q$ is expressed as follows:

$$\sum _{n=1}^N\sum _{k=1}^K{A}_0{a}_{n,q}{e}^{j[{\varphi }_0+{\omega }_0(d_{n,k}+g_n)]}={\gamma }_q$$

(8)

Utilizing the mathematical relationships in the coherent signal reception model and the space–time–coding matrix, we can recover the observation vector containing the angular information of coherent signals. By taking Q harmonics on either side of the fundamental wave, where Q is the number of harmonics on each side, there are a total of $2Q+1$ spectral amplitude values, which are indexed as $q=-Q,-Q+1,\dots ,0,\dots ,Q-1,Q$. The mathematical relations (8) and (9) in the coherent signal reception model can be expressed in matrix form as follows:

$$\left[\begin{array}{cccc}{a}_{1,-Q}& a_{2,-Q}& \cdots & a_{N,-Q}\\ a_{1,-Q+1}& a_{2,-Q+1}& \cdots & a_{N,-Q+1}\\ ⋮& ⋮& \ddots & ⋮\\ a_{1,0}& a_{2,0}& \cdots & a_{N,0}\\ ⋮& ⋮& \ddots & ⋮\\ a_{1,Q-1}& a_{2,Q-1}& \cdots & a_{N,Q-1}\\ a_{1,Q}& a_{2,Q}& \cdots & a_{N,Q}\end{array}\right]\left[\begin{array}{c}{\displaystyle \sum _{k=1}^K}{e}^{j{\omega }_0(d_{1,k}+g_1)}\\ {\displaystyle \sum _{k=1}^K}{e}^{j{\omega }_0(d_{2,k}+g_2)}\\ ⋮\\ {\displaystyle \sum _{k=1}^K}{e}^{j{\omega }_0(d_{N,k}+g_N)}\end{array}\right]={\displaystyle \frac{e^{-j{\varphi }_0}}{A_0}}\left[\begin{array}{c}{\gamma }_{-Q}\\ {\gamma }_{-Q+1}\\ ⋮\\ {\gamma }_0\\ ⋮\\ {\gamma }_{Q-1}\\ {\gamma }_Q\end{array}\right]$$

(9)

From Equation (10), it can be observed that the matrix on the left side of the equation is composed of Fourier coefficients representing the fundamental and harmonic components. We refer to this matrix as the space–time–coding matrix. The values of each element in this space–time–coding matrix are decided by the space–time–coding modulation sequence of the metasurface linear array model. The observation vector on the left side of the equation contains the angular information of coherent signals, while the vector on the right side of the equation is composed of the amplitude values of the fundamental and harmonic components of the spectrum obtained by the spectrum analyzer. The space–time–coding matrix $\mathit{A}$ is given by

$$\mathit{A}=\left[\begin{array}{cccc}{a}_{1,-Q}& a_{2,-Q}& \cdots & a_{N,-Q}\\ a_{1,-Q+1}& a_{2,-Q+1}& \cdots & a_{N,-Q+1}\\ ⋮& ⋮& \ddots & ⋮\\ a_{1,0}& a_{2,0}& \cdots & a_{N,0}\\ ⋮& ⋮& \ddots & ⋮\\ a_{1,Q-1}& a_{2,Q-1}& \cdots & a_{N,Q-1}\\ a_{1,Q}& a_{2,Q}& \cdots & a_{N,Q}\end{array}\right]$$

(10)

Let

$$\psi ={\left[K,\sum _{k=1}^K{e}^{j{\omega }_0{d}_{2,k}},\dots ,\sum _{k=1}^K{e}^{j{\omega }_0{d}_{N,k}}\right]}^{\mathrm{T}}$$

(11)

$$\delta ={\left[\begin{array}{c}{\displaystyle \sum _{k=1}^K}{e}^{j{\omega }_0{g}_1},{\displaystyle \sum _{k=1}^K}{e}^{j{\omega }_0{g}_2},\cdots ,{\displaystyle \sum _{k=1}^K}{e}^{j{\omega }_0{g}_N}\end{array}\right]}^{\mathrm{T}}$$

(12)

$$\gamma ={\left[\begin{array}{c}{\gamma }_{-Q},{\gamma }_{-Q+1},\cdots ,{\gamma }_0,\cdots ,{\gamma }_{Q-1},{\gamma }_Q\end{array}\right]}^{\mathrm{T}}$$

(13)

In the following matrix operations, the element-wise product is defined as ⊙ and the element-wise division is defined as ⊘. These operations apply to corresponding elements of the matrices. Therefore, Equation (10) can be simplified as

$$\mathit{A}\psi \odot \delta ={\displaystyle \frac{e^{-j{\varphi }_0}}{A_0}}\gamma$$

(14)

The space–time coding matrix $\mathit{A}$ is an $(2Q+1)\times N$ matrix, and its generalized inverse matrix is ${\mathit{A}}^{-1}$. The vector of the path differences from each subarray to the horn receiving antenna is $\delta$. Thus, the observation vector $\psi$, which contains the angular information of coherent signals, is given by

$$\psi ={\displaystyle \frac{e^{-j{\varphi }_0}}{A_0}}{\mathit{A}}^{-1}\gamma \oslash \delta$$

(15)

Building upon the coherent signal reception model developed in the previous section, we have derived the observation vector that requires a solution. This vector encapsulates the angular information of multiple coherent signals with the angles being sparsely distributed in the spatial domain. A compressed sensing algorithm, specifically the ℓ₁-SVD algorithm, can be employed for the DOA estimation of coherent signals.

From a sparse perspective, considering the observation vector, we can write

$$\psi =\mathit{Zx}+\mathit{b}$$

(16)

The observation vector $\psi$ is of size $N\times 1$, where N is the total number of linear array elements. K is the total number of coherent signals, satisfying $N>>K$, for instance, $N=10K$. The overcomplete dictionary, denoted as $\mathit{Z}$, is composed of a set of grids in the spatial domain, with dimensions $N\times N_{\theta }$, where $N_{\theta }$ represents the total number of spatial domain grids. The predefined angle range is given by

$$S_{\theta }=\left\{-{90}^{\circ }:0.1^{\circ }:{90}^{\circ }\right\}$$

(17)

Consequently, $N_{\theta }=1801$. The sparse signal vector $\mathit{x}$ of size $N_{\theta }\times 1$ corresponds to the sparse signal. $\mathit{b}$ represents the noise term of size $N\times 1$. Given that the angles to be estimated are sparse in the spatial domain relative to the divided angles, most elements of the sparse signal vector $\mathit{x}$ are zero with non-zero elements at the positions of the sparse signals. By identifying the index positions of the non-zero elements, we can estimate the DOA of the signals in $S_{\theta }$.

Due to the influence of noise, additional non-zero terms are inevitably introduced. To address this, we minimize the noise parameter and employ ℓ₁ norm minimization in the context of compressed sensing algorithms. Thus, the sparse signal vector $\mathit{x}$ can be reconstructed using the following Lagrangian form of the objective equation:

$$min{∥\psi -\mathit{Zx}∥}_2^2+\lambda {∥\mathit{x}∥}_1$$

(18)

Here, $\mathit{x}$ is the reconstructed sparse vector, $\lambda$ is the regularization parameter used to balance the sparsity between $\mathit{x}$ and the residual $\psi -\mathit{Zx}$, and ${∥*∥}_2$ denotes the ℓ₂ norm, while ${∥*∥}_1$ denotes the ℓ₁ norm.

The covariance matrix $\mathit{Y}$ of the observation vector $\psi$ is calculated as

$$\mathit{Y}=\psi {\psi }^{\mathrm{H}}$$

(19)

Dimensionality reduction is performed on $\mathit{Y}$ using Singular Value Decomposition (SVD):

Here, $\mathit{U}$ is of size $N\times N$, $\mathit{V}$ is of size $T\times T$, both are unitary matrices, and $\mathit{D}$ is of size $N\times T$, which is a singular value matrix. The selection matrix ${\mathit{D}}_{\mathit{K}}$ of size $N\times K$ is constructed from a K-dimensional identity matrix and a $(T-K)\times K$ zero matrix. N and T denote the dimensions of the matrix after singular value decomposition, and K represents the assumed number of coherent signals, which is typically less than N. This results in the dimensionally reduced observation matrix ${\mathit{Y}}_{\mathit{SV}}$ of size $N\times K$:

$${\mathit{Y}}_{\mathit{SV}}={\mathit{YVD}}_{\mathit{K}}={\mathit{UDD}}_{\mathit{K}}$$

(20)

The dimensionally reduced covariance matrices of the sparse vector $\mathit{x}$ and the noise, denoted as ${\mathit{X}}_{\mathit{SV}}$ and ${\mathit{b}}_{\mathit{SV}}$, respectively, are defined as

$${\mathit{X}}_{\mathit{SV}}={\mathit{xx}}^{\mathrm{H}}{\mathit{VD}}_{\mathit{K}}$$

(21)

$${\mathit{B}}_{\mathit{SV}}={\mathit{bb}}^{\mathrm{H}}{\mathit{VD}}_{\mathit{K}}$$

(22)

Based on the sparse angle notation in Equation (14), the observation matrix ${\mathit{Y}}_{\mathit{SV}}$ can be expressed as

$${\mathit{Y}}_{\mathit{SV}}={\mathit{ZX}}_{\mathit{SV}}+{\mathit{B}}_{\mathit{SV}}$$

(23)

The aforementioned dimensionality reduction process transforms the optimization problem in Equation (16) as follows:

$$min{∥{\mathit{Y}}_{\mathit{SV}}-\mathit{Z}{\mathit{X}}_{\mathit{SV}}∥}_2^2+\lambda {∥{\mathit{X}}_{\mathit{SV}}∥}_{2,1}$$

(24)

Here, ${\mathit{X}}_{\mathit{SV}}$ now represents a matrix relative to $\mathit{x}$, and ${∥*∥}_{2,1}$ denotes the ℓ₁ norm of the ℓ₂ norm of the column vector in the matrix, serving to replace each element in $\mathit{x}$ [35]. This transformation enhances the efficiency of the problem-solving process. The objective equation constitutes a second-order cone programming problem, which is a type of convex optimization problem. It can be solved using MATLAB’s convex optimization toolbox, cvx, specifically for second-order cone programming. By minimizing the ℓ₂ norm of ${\mathit{Y}}_{\mathit{SV}}-\mathit{Z}{\mathit{X}}_{\mathit{SV}}$, we find ${\mathit{X}}_{\mathit{SV}}$ that satisfies the conditions. Subsequently, by calculating the ℓ₂ norm of each row vector of ${\mathit{X}}_{\mathit{SV}}$, we obtain an N-dimensional column vector. The optimization process of a convex function can be achieved by solving the following equation. The problem includes the following three constraints:

$$\begin{array}{c}{\displaystyle \underset{p,q,\mathit{r},{\mathit{X}}_{\mathit{SV}}}{min}}p+\lambda q\\ s.t\parallel \mathrm{vec}\{{\mathit{Y}}_{\mathit{SV}}-{\mathit{ZX}}_{\mathit{SV}}\}{\parallel }_2^2\le p\\ {\mathbf{1}}^T\mathit{r}\le q\\ \parallel {\mathit{X}}_{\mathit{SV}}\left(i,:\right){\parallel }_2\le {\mathit{r}}_i\end{array}$$

(25)

At this point, we have obtained the optimized ${\mathit{X}}_{\mathit{SV}}$ matrix. Each row of the ${\mathit{X}}_{\mathit{SV}}$ matrix represents a potential DOA result. Specifically, the first column of the ${\mathit{X}}_{\mathit{SV}}$ matrix contains the primary signal directional information. Therefore, by plotting the spatial spectrum of only the first column of the ${\mathit{X}}_{\mathit{SV}}$ matrix, we can obtain the DOA estimation results for coherent signals. The angles within the predefined angle range $S_{\theta }$ corresponding to the indices of the K non-zero values in this vector are identified as the directions of arrival of the coherent signals. We summarize the steps of the algorithm proposed in Algorithm 1.

Algorithm 1 STCM-Based ℓ1-SVD DOA Estimation

Require:  Time domain sampling signal s Ensure:  Estimated DoAs ${\widehat{\theta }}_k,k=1,2,\dots ,K$.   1:   Calculate FFT of received signal s to obtain $s_{FFT}$;   2:   Construct coding matrix $\mathit{A}$ using space–time–coding sequence;   3:   Extract frequency components $\gamma$ at $F_c+q{F}_p(q=-Q,\cdots ,Q)$ and path differences $\delta$;   4:   Calculate equivalent array manifold $\psi$ using $\mathit{A}$, $\gamma$ and $\delta$;   5:   Form signal covariance matrix $\mathit{Y}$ = $\psi {\psi }^{\mathrm{H}}$;   6:   Perform SVD on $\mathit{Y}$ to obtain $\mathit{U}$, $\mathit{D}$, ${\mathit{D}}_{\mathit{K}}$, $\mathit{V}$, ${\mathit{Y}}_{\mathit{SV}}$ and ${\mathit{X}}_{\mathit{SV}}$;   7:   Construct overcomplete dictionary $\mathit{Z}$;   8:   Obtain the optimization problem $min{∥{\mathit{Y}}_{\mathit{SV}}-\mathit{Z}{\mathit{X}}_{\mathit{SV}}∥}_2^2+\lambda {∥{\mathit{X}}_{\mathit{SV}}∥}_{2,1}$;   9:   Solve ℓ1-SVD optimization problem in (25);   10: Calculate spatial spectrum from the first column of the ${\mathit{X}}_{\mathit{SV}}$;   11: Find K highest peaks in spatial spectrum to obtain ${\widehat{\theta }}_k$.

4. Simulation Results

This section displays the simulation results of the proposed DOA estimation method for coherent signals based on the STCM. The simulation experiments are run to validate the performance of the algorithm in estimating the DOA of coherent signals. Firstly, the performance of the proposed method in DOA estimation for coherent signals at different angles is verified through simulation. Then, the performance of the proposed method is evaluated by studying the impact of different parameters on the accuracy of DOA estimation, such as the number of elements, snapshots, and harmonics. Finally, the performance advantages of the proposed method are verified by comparing it with existing methods.

4.1. Performance Verification

Firstly, simulation experiments are conducted according to the parameter settings in

Table 1. Parameter settings.

Parameters	Values
Equivalent Number of Elements N	8
Modulation Frequency $F_p$	50 MHz
Number of Snapshots $N_s$	400
One-Sided Harmonic Number Q	4
Signal-to-Noise Ratio $SNR$	20 dB
Incident Angles ${\theta }_k$	$0^{\circ }$ and $5^{\circ }$

, simulating two coherent signal sources illuminating the metasurface from the far field. The size of the metasurface is set to $8\times 8$ elements with an element spacing which is half a wavelength ${\lambda }_0/2$. The receiving horn antenna is placed over the center of the metasurface at a height of 0.175 m. To simplify the calculations, a one-dimensional spatial DOA estimation experiment is initially conducted, treating the two-dimensional metasurface as a one-dimensional linear array. The modulation frequency $F_p$ is set to 50 MHz, resulting in a modulation period $T_p$ of 0.02 $\mathsf{\mu }$s. The number of sampling snapshots $N_s$ is set to 400, the SNR is set to 20 dB, and the angles of the two incident coherent signals are $0^{\circ }$ and $5^{\circ }$.

The STCM is modulated column-wise, as detailed in Equation (4) and shown in Figure 2. The spectrum of the modulated coherent signals is obtained by the FFT, as shown in Figure 3. The fundamental spectral line can be clearly visible at the carrier frequency, with harmonic spectral lines located on both sides of the fundamental spectral line, spaced at integer multiples of the modulation frequency $F_p$. This indicates that the modulation of the STCM has been successfully implemented, which is essential for the subsequent processing steps.

With four single-sided harmonics considered, we have a total of nine amplitude values. Using the algorithm mentioned in Section 3, the relationship between the spectral amplitude values and the coherent incident angles is evaluated.

To assess the feasibility of the proposed coherent signal DOA estimation method, we evaluated its performance across various angles. Firstly, we set several sets of coherent signal angle values with different intervals. The angles are set to be (−48°, −27°), (−36°, −19°), (−8°, 10°), (5°, 18°), and (29°, 41°) with unequal intervals between each pair.

The estimation results are shown in Figure 4, and the absolute errors of these results are provided in

Table 2. Absolute errors of the estimation results.

Preset Angles [deg]	Estimation Results [deg]	Errors [deg]
(−48, −27)	(−49.2, −25.7)	(1.2, 1.3)
(−36, −19)	(−37.7, −17.2)	(1.7, 1.8)
(−8, 10)	(−8.3, 10.2)	(0.3, 0.2)
(5, 18)	(2.5, 20.4)	(2.5, 2.4)
(29, 41)	(28.9, 41.1)	(0.1, 0.1)

. It is apparent from the data that the absolute angular estimation error for both coherent signals is within 2.4°. This confirms that the proposed method performs well in DOA estimation.

To further investigate the DOA estimation performance of the proposed method for coherent signals at various angles, a series of experiments are conducted. The coherent signals with different angles are tested, while the angle interval is kept fixed. The performance is evaluated using the root mean square error (RMSE), where ${\widehat{\theta }}_k$ represents the DOA estimation result for each experiment, K denotes the number of incident signals, and the number of Monte Carlo trials is 100.

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{100K}}\sum _{i=1}^{100}\sum _{k=1}^K{\left({\widehat{\theta }}_k-{\theta }_k\right)}^2}$$

(26)

The results are presented in Figure 5, where the angle interval is fixed at 5°, and one of the angles varies from −75° to 75°. It is found that the DOA estimation error increases with the absolute value of the incident angle, which is due to the characteristics of the metasurface as an array. At larger incident angles, the modulation performance deviation of the STCM is exacerbated, whereas areas closer to the metasurface’s normal line exhibit lesser effects. Therefore, under conditions of smaller incident angles, the error of the DOA estimation is smaller.

4.2. Performance Analysis

After verifying the preliminary feasibility of the proposed coherent signal DOA estimation method, the DOA estimation performance characteristics of the method will be analyzed in the following sections. The main aspects of analysis include the number of elements, the number of snapshots, and the number of harmonics.

4.2.1. Performance Analysis for the Number of Elements

In the first part, Figure 6 shows the impact of different numbers of elements on DOA estimation performance. In the experiment, the number of harmonics is fixed at 4, the number of snapshots is fixed at 400, the SNR is set from −10 dB to 20 dB, and other parameters remain unchanged. The performance of DOA estimation is measured by RMSE. An increase in the number of elements from 8 to 20 results in a notable decrease in RMSE. This suggests that enhancing the number of elements significantly boosts the DOA estimation accuracy. Specifically, when the number of elements increases from 8 to 10, there is a significant decrease in RMSE. As the number of elements further increases to 16 and 20, the RMSE continues to decrease, but the magnitude of the decrease gradually slows down. This trend underscores that while augmenting elements can markedly enhance DOA estimation performance, the marginal gains from additional elements will incrementally diminish. However, adding more elements will also increase the complexity of the system.

4.2.2. Performance Analysis for the Number of Snapshots

In the second part, Figure 7 shows the impact of different snapshot numbers on DOA estimation performance. In the experiment, the number of harmonics is fixed at 4, the number of elements is fixed at 8, and other parameters remain unchanged. As the number of snapshots increased from 400 to 2000, the RMSE shows a significant downward trend. This indicates that increasing the number of snapshots can significantly improve the accuracy of DOA estimation. Specifically, when the number of snapshots increases from 400 to 1600, there is a significant decrease in RMSE. As the number of snapshots further reaches 1600 and 2000, the RMSE continues to decrease, but the magnitude of the decrease gradually reduces. This indicates that augmenting the snapshot count can markedly improve DOA estimation performance, yet the incremental benefits of additional snapshots will eventually experience diminishing returns. However, increasing the number of snapshots also lengthens the sampling time of the system, which increases the computational cost.

4.2.3. Performance Analysis for the Number of Harmonics

In the third part, the impact of the number of harmonics on the DOA estimation performance of the proposed method is explored. The results are presented in Figure 8. We varied the number of single-sided harmonics from one to eight, keeping all other parameters constant. As the number of harmonics increases from one to three, the RMSE shows a significant decreasing trend. It is observed that when the number of harmonics reached four, the RMSE no longer decreases significantly, and accurate DOA estimation results are obtained. In addition, when the SNR is low, the influence of noise causes the harmonic error of the 5th order and higher to slightly increase the RMSE of DOA estimation. This indicates that even with a moderate number of harmonics, the proposed method maintains high accuracy. The simulation results support the effectiveness of the proposed direction estimation method based on the selection of harmonics.

4.3. Performance Comparison

Firstly, to compare the DOA estimation resolution of the algorithm proposed in this paper with the existing spatial smoothing MUSIC (SS-MUSIC) algorithm, a comparative experiment was conducted at two angles, 0° and 5°, under the same angular interval.

The results, as shown in the Figure 9 and Figure 10, indicate that the method proposed in this paper can resolve the DOA of two coherent signals with a 5° interval, while the spatial smoothing MUSIC algorithm struggles to do so. Therefore, this experiment demonstrates that the method proposed in this paper has an advantage in angular resolution.

Proceeding to the subsequent part, the performance of the proposed method for the DOA estimation of coherent signals is compared with OMP and SS-MUSIC algorithms. In the experiment, the number of harmonics is fixed at four, the number of elements is fixed at eight, the number of snapshots is fixed at 400, and all other parameters are kept constant. The results are presented in Figure 11. The SNR is adjusted from −10 dB to 20 dB in increments of 5 dB with all other parameters remaining constant. The RMSE is employed to assess the DOA estimation performance. As the SNR ranged from −10 dB to 20 dB, the RMSE of the proposed method showed a decreasing trend. In comparison, the RMSE of the other algorithms did not decrease as significantly. This indicates that even at lower SNR levels, the proposed method maintains high accuracy, outperforming the other algorithms. The simulation results confirm the effectiveness and superiority of the proposed direction estimation method for coherent signals based on the STCM.

Similarly, we will continue to compare the performance of different algorithms in the DOA estimation of coherent signals. The variable in focus is now the number of snapshots, which is increased from 400 to 2000, while all other parameters remain unchanged. The results are presented in Figure 12.

As the number of snapshots increases, the RMSE of different methods shows a slight downward trend. In comparison, the proposed method also maintains high accuracy and outperforms other algorithms with varying numbers of snapshots. The simulation results further validate the effectiveness and superiority of the algorithm presented in this paper for coherent signal direction estimation based on the STCM.

Finally, we compared the running times of the three algorithms, as shown in

Table 3. The comparison of the processing time for these three methods.

Algorithm	Processing Time [s]
SS-MUSIC	0.027269
OMP	0.026063
ℓ1-SVD	1.412667

. This paper’s proposed algorithm provides high precision and resolution for coherent signal direction estimation using the STCM, but it necessitates a longer optimization duration.

5. Conclusions

In conclusion, this paper presents a thorough analysis of the DOA estimation performance of the proposed method for coherent signals under various conditions. The results confirm that the proposed method exhibits good performance across different angular spacings, numbers of harmonics, and SNR levels, highlighting its robustness and accuracy. This paper introduces a novel and efficient method for estimating the DOA of coherent signals by leveraging the modulation capabilities of the STCM. The proposed method modulates incident electromagnetic waves and processes the received signals to achieve high-precision DOA estimation. Simulation results validate the effectiveness of the proposed method under various conditions, including different incident angles and SNR levels, and show its superiority over traditional algorithms. Future work will focus on extending the proposed method to more complex scenarios, exploring its potential in real-world communication systems, enhancing its performance, and integrating it with other advanced signal processing techniques.