Novel time-varying metasurface design for enhanced DOA estimation based on sparse signal recovery

Introduction: Metasurfaces can help innovate traditional estimation algorithms as an excellent alternative to phased arrays. Direction of arrival (DoA) estimate is an important research area in array signal processing, and various progressive direction-finding methods have already been created. Existing time-domain DoA estimation methods have the problems of spectral aliasing in target signals and constraints of temporal modulation.

Methods: Thus, this paper proposes a novel time-varying metasurface design with prominent feature of asynchronous control-based DoA estimation. The proposed metasurface mitigates the spectral aliasing, enhances the signal processing bandwidth and improves the modulation rate. It also utilizes the persistence of unit states over a period of time which is inherent and it does not require amendments to the hardware constraints. The accuracy of DoA estimation is improved by deploying asynchronous modulation which effectively increases the number of virtual channels under the constraints of materials characteristics.

Results and Discussion: Simulation results shows that the proposed metasurface has better performance as compared with existing methods. Also, the results obtained have closed approximation to the theoretical optimum which further validated its effectiveness.

1 Introduction

Direction of arrival (DOA) estimation is widely used to obtain user location and thus improve communication quality [1, 2]. As a key technology in the field of radar and wireless communication, DOA estimation has been widely studied in the fields of array signal processing, sensor networks, remote sensing, etc., and plays an important role in practical applications. For example, in satellite mobile communications, accurate DOA estimation can assist in forming high-precision directional beams, improving communication quality and efficiency [3, 4]. In terms of smart antenna technology, DOA estimation can be used for beamforming to guide signals to achieve high-capacity communication, which is crucial to improving the performance of communication systems [5].

In past studies, the multiple signal classification (MUSIC) algorithm [6–9], the estimating signal parameter via rotational invariance technique (ESPRIT) algorithm [10–12] and related improved subspace algorithms [13, 14] have been used to provide high-precision DOA estimation results. At the same time, research on compressed sensing algorithms represented by orthogonal matching pursuit [15] is emerging. These classic DOA estimation techniques usually rely on multi-channel antenna arrays [16] to ensure performance, and the most direct way to improve resolution is to increase the number of antennas. However, the increase in the number of antennas will bring about problems such as large structure size, complex feeding circuit, increased manufacturing cost and power consumption.

In recent years, metasurfaces, as a new type of material that can actively control the amplitude, phase and polarization of electromagnetic waves, thereby affecting the propagation characteristics of electromagnetic waves [17–19], have attracted widespread attention due to their advantages such as small size and easy processing. In the research based on metasurfaces, Ref. [20] introduced the time dimension and proposed the concept of time varying metasurface (TVM), which provides a new degree of freedom for electromagnetic wave control. The harmonic effect brought by time control can well realize the construction of multi-dimensional receiving space and achieve the effect of virtual multi-channel reception. The TVM receives multi-order harmonics of the signal by generating multi-modal and low-correlation directional patterns, and only needs a single channel to realize DOA estimation [21–23], which fundamentally reduces the hardware complexity.

Reference [24] proposed a DOA estimation method based on TVM. This method realizes DOA estimation by analyzing the spectrum scattering characteristics of TVM. Reference [25] realizes DOA estimation by analyzing the amplitude imbalance characteristics of the received signal at two first-order harmonic frequencies. Reference [26] proposed a DOA estimation method based on time control theory. This method uses the harmonic characteristic matrix, combined with time control and phase control, to reconstruct the manifold vector of the array. Reference [27] modulates the incident signal in the time domain by designing the coding sequence of the TVM unit, and uses harmonics to restore the multi-channel receiving signal matrix, thereby realizing DOA estimation.

Most of the existing TVM-based DOA estimation methods do not consider the impact of the incident signal bandwidth, and directly assume that the incident signal can obtain virtual multi-channels after time modulation, thereby realizing DOA estimation. However, when the incident signal irradiates the metasurface with a frequency of f_c, the periodic time modulation will cause the metasurface to generate a large number of harmonic components, which are emitted in beams in different directions. The frequency interval of the harmonic components is determined by the control frequency of the TVM. When the control frequency is lower than the incident signal bandwidth, the generated harmonics will be aliased, resulting in the inability to perform DOA estimation. Metasurfaces made of different materials have different electromagnetic wave control capabilities. Metasurfaces based on varactor diodes have the ability to continuously adjust amplitude, phase, and polarization, but their adjustment speed is slow, and a state change requires hundreds of nanoseconds. In contrast, metasurfaces based on PIN diodes can only switch between several discrete adjustment parameters, but their state change speed is very fast, and can be achieved within 100 ns [28]. Therefore, even for the metasurface with the highest control rate, the upper limit of the incident signal bandwidth is less than 5 MHz. However, with the development of wireless communication technology, the bandwidth requirements of communication signals will become higher.

In view of the above problems, this paper designs a time-controlled metasurface DOA estimation method based on asynchronous control. The main contributions are outlined as follows:

1) Without changing the hardware constraints of TVM, the property that the unit state will last for a period of time is used to stagger the change start time of different columns of units to obtain multiple different responses within a response duration.

2) This asynchronous control method enables TVM to increase the number of virtual multi-channels equivalently under material constraints and improve the DOA estimation accuracy. The simulation results verify the effectiveness of this method.

3) Compared with the existing synchronous control method, the DOA estimation performance has been greatly improved and can approach the theoretical optimal DOA estimation result.

The remaining of this paper is organized as follows. In Section 2, the metasurface model is designed and discussed. In Section 3, the proposed methodology is discussed. In Section 4, the simulation experimentations are discussed. In Section 5, the conclusion is described.

2 Metasurface model for received signal

Consider a TVM with M columns for signal reception, where the spacing between adjacent column units is d. The single-channel TVM signal reception model is shown in Figure 1. Starting from the one-dimensional metasurface, assuming that each column unit is controlled by the same signal, the entire metasurface can be equivalent to an antenna array with an array unit spacing of d and an array element number of M.

Figure 1

Figure 1. Metasurface model for single channel received signal.

Assume that the far-field transmitter transmits L independent narrowband signals of the same frequency $S\left(t\right)={\displaystyle \sum _{l=1}^L}{s}_l\left(t\right)e^{2\pi f_{c}t}$, where $s_l\left(t\right)$ is the baseband signal from ${\theta }_l$, $j$ is the imaginary unit, and $f_c$ is the carrier frequency. The received signal of the mth column unit can be expressed in Equation 1:

$x_m\left(t\right)=S\left(t\right)e^{j\beta \left(m-1\right)d\sin {\theta }_l}+n_m\left(t\right)(1)$

Among them, $\beta =2\pi /\lambda$, which represents the wave number of frequency $f_c$; $n_m\left(t\right)$ represents zero-mean Gaussian white noise; The variance is ${\sigma }^2$. Therefore, in the single-channel model, the received signal $y\left(t\right)$ can be expressed in Equation 2:

$y\left(t\right)={\displaystyle \sum _{m=1}^M}{\mathrm{\Gamma }}_m\left(t\right)x_m\left(t\right)(2)$

Among them, ${\mathrm{\Gamma }}_m\left(t\right)$ is the periodic time control function in the mth column, which can be expressed as a linear combination of harmonics of various orders which is expressed by Equation 3:

${\mathrm{\Gamma }}_m\left(t\right)={\displaystyle \sum _{k=-\infty }^{+\infty }}{c}_{mk}{e}^{j2k\pi f_{0}t}(3)$

Where, $f_0$ represents the control frequency, $c_{mk}$ is the kth order Fourier coefficient of ${\mathrm{\Gamma }}_m\left(t\right)$, expressed in Equation 4:

$c_{mk}=\frac{1}{T_0}{\int }_0^{T_0}{\mathrm{\Gamma }}_m\left(t\right)e^{-j2k\pi f_{0}t}\mathrm{d}t(4)$

Among them, $T_0=1/f_0$, represents the control period. In summary, the received signal can be expressed by Equation 5:

$y\left(t\right)={\displaystyle \sum _{k=-\infty }^{+\infty }}{\displaystyle \sum _{m=1}^M}{x}_m\left(t\right)c_{mk}{e}^{j2k\pi f_{0}t}+n\left(t\right)(5)$

It can be seen from Equation 5 that due to the time control of the metasurface, the incident signal will produce a nonlinear phenomenon, which is characterized by the generation of a large number of harmonics in addition to the fundamental component. That is, the signal spectrum is distributed on an infinite number of frequency points $f_c+k{f}_0,k\in \left\{-\infty ,+\infty \right\}$. The interval between each frequency component is determined by $f_0$. According to the Nyquist sampling theorem, $f_0$ must be greater than the incident signal bandwidth $B_x$ to ensure that spectrum aliasing does not occur. The phase angle is expressed in Equation 6:

${\phi }_k\left(\theta \right)={\displaystyle \sum _{m=1}^M}{c}_{mk}{e}^{j\beta \left(m-1\right)d\sin \theta }(6)$

The schematic diagram of TVM’s virtual multi-channel principle is shown in Figure 2. In one symbol period, the TVM unit changes N times, which is equivalent to N heterogeneous array patterns observing the signal separately, thus constructing an N-dimensional virtual channel. Under ideal conditions, in order to fully utilize the TVM array’s degrees of freedom (DoF), the number of states N needs to be no less than the number of arrays M.

Figure 2

Figure 2. Metasurface model for virtual multi-channel received signal.

For an M-column TVM, when N ≥ M, that is, when the number of virtual channels is not less than the number of metasurface columns, the theoretically optimal DOA estimation result can be obtained. When 2 ≤ N < M, the number of channels is insufficient, resulting in insufficient estimation freedom, and the DOA estimation accuracy will decrease as the number of channels N decreases, and the maximum number of estimated sources becomes N - 1. When the number of states N is the minimum value 2, the TVM control rate achieves the theoretical maximum value ${\left[f_0\right]}_{\max }=f_v/2$, where $f_v$ is the unit control frequency. When N < 2, the TVM unit cannot change the response within the duration of a symbol [29]. At this time, TVM loses its meaning and the traditional single-channel algorithm needs to be used to achieve DOA estimation [30]. Therefore, DOA estimation based on TVM can be divided into the following three cases according to the relationship between the control rate and the incident signal bandwidth.

(1) When $B_x>f_v/2$, the spectrum will be aliased when TVM is used to receive the signal, and the traditional single-channel algorithm needs to be used to achieve DOA estimation.

(2) When $f_v/M<B_x\le f_v/2$, the TVM can meet the minimum rate requirement of periodic control, but the number of virtual channels is less than the number of TVM columns, and the spatial resources of multiple antenna units cannot be fully utilized, resulting in a decrease in estimation performance.

(3) When $B_x\le f_v/M$, the TVM unit can change $M$ times within the duration of a code element, and the theoretically optimal number of virtual channels can be obtained. The DOA estimation performance can approach that of a multi-channel array.

3 Proposed methodology

According to the analysis in Section 1, when the synchronous control method switches the state of the TVM unit, the starting time of each column unit is the same. When the control rate is determined, the number of states N is also determined. When $f_v/M<B_x\le f_v/2$, the number of states N is less than the number of TVM columns M, and the number of equivalent channels is insufficient, and the best DOA estimation effect cannot be achieved. Therefore, in order to make full use of the airspace resources of multiple antenna units of TVM, this paper proposes an asynchronous control method for the DOA estimation scenario when $f_v/2\le B_x<f_v/M$.

3.1 Asynchronous control method

Considering the nature that each unit response of TVM will last for a period of time, this paper staggers the start time of unit changes in different columns so that the start time of each unit control of TVM is spaced out within a response duration. At this time, the duration of a single state is shortened, and multiple different responses are realized within a response duration, which is equivalent to increasing the number of virtual channels and improving the DOA estimation accuracy. The schematic diagram of the asynchronous control method is shown in Figure 3.

Figure 3

Figure 3. Schematic diagram of the asynchronous control method principle.

In the asynchronous control method, it is assumed that each TVM unit has N different states in each control cycle, the duration of a state is $T_0/N$, and the control start time difference between adjacent columns is $∆t$. To ensure encoding independence, Equations 7, 8 are deployed:

$M∆t=k\frac{T_0}{N},k\in \left\{1,2,\dots ,M/2\right\}(7)$

$∆t=\frac{T_v}{k},k\in \left\{2,\dots ,M\right\}(8)$

By using the asynchronous control method, the number of original states can be increased to $k$ times within one control cycle.

For ease of understanding, the subsequent analysis is based on the PIN diode TVM for theoretical explanation. Each unit controls the reflection coefficient of the unit by applying a control voltage to the PIN diode to achieve 1 bit of control capability. At this time, the time control function of the TVM can be expressed as: A periodic square wave signal, including high level “1” and low level “0”. The schematic diagram of Hadamard matrix asynchronous control is shown in Figure 4. It shows the state of the TVM after asynchronous control when $N=4,M=8,k=4$, and the Hadamard matrix is ​​selected as the original state encoding matrix. In Figure 4, the dotted line represents the level duration, and the origin represents the state. In synchronous control, the state matrix is ​​encoded as “1 0 1 0 1 0 1 0”, “1 0 0 1 10 0 1”, “1 0 1 0 0 1 0 1”, and “1 0 0 1 0 1 1 0”. In asynchronous control, four groups of codes are added: “1 0 1 1 1 1 1 0”, “10 0 0 1 0 0 0”, “1 0 1 1 0 0 0 1”, and “1 0 0 0 0 1 11”. As shown in Figure 4, the state encoding matrix after asynchronous control has four more columns than before control, the duration of the state is shortened, and the number of virtual multi-channels is equivalently increased.

Figure 4

Figure 4. Schematic diagram of Hadamard matrix asynchronous control.

3.2 Sparse representation of received signal

The single-channel TVM receiving signal $y_n\left(t\right)$ in the nth state is expressed in Equation 9:

$y_n\left(t\right)={\displaystyle \sum _{l=1}^L}{\displaystyle \sum _{m=1}^M}{B}_{mn}S\left(t\right)e^{j\beta \left(m-1\right)d\sin {\theta }_l}+n\left(t\right)(9)$

Among them, $B_{mn}={\left(-1\right)}^{b_{mn}}$ represents the response coefficient of the mth unit in the nth state. When the state is low level “0”, $b_{mn}=0$, which means it is in the opposite direction of the signal. When the state is high level “1”, $b_{mn}=1$, which means it is in the same direction as the signal. The received signal is expressed in Equation 10:

$\mathit{Y}\left(n\right)={\mathit{B}}^{\mathrm{T}}\mathit{A}\mathit{S}\left(n\right)+\mathit{N}\left(n\right)(10)$

Among them, $\mathit{Y}\left(n\right)={\left[y_1\left(n\right),\dots ,y_N\left(n\right)\right]}^{\mathrm{T}}$ is the received signal matrix composed of the received signals in each state; $\mathit{B}=\left\{{\left(-1\right)}^{b_{mn}}\right\}\in {\mathbb{C}}^{M\times N}$ is the reflection coefficient matrix; $\mathit{A}={\left[\mathit{\alpha }\left({\theta }_1\right),\dots ,\mathit{\alpha }\left({\theta }_L\right)\right]}^{\mathrm{T}}$ is the array flow pattern matrix; $\mathit{\alpha }\left(\theta \right)={\left[1,\dots ,e^{j\beta d\left(M-1\right)\sin \theta }\right]}^{\mathrm{T}}$ is the far-field signal steering vector; $\mathit{S}\left(n\right)={\left[s_1\left(n\right),\dots ,s_L\left(n\right)\right]}^{\mathrm{T}}$ is the incident signal matrix; $\mathit{N}\left(n\right)={\left[n_1\left(n\right),\dots ,n_N\left(n\right)\right]}^{\mathrm{T}}$ is the noise signal matrix.

The array flow matrix $\mathit{A}$ contains the directional information $\theta$ of all far-field signals. At the same time, for the entire signal space, the actual signal’s incoming direction is sparse. The schematic diagram of the incident signal’s spatial sparsity is shown in Figure 5. Among them, the hollow represents the possible incident direction of the signal, the solid represents the signal incident, and the signal’s incoming direction is ${\theta }_l$.

Figure 5

Figure 5. Schematic diagram of spatial sparsity of incident signal.

Therefore, by evenly dividing the space into $\left\{{\theta }_1,\dots ,{\theta }_K\right\}$ and assuming that there is a possible potential signal in each of the above directions, we can obtain an overcomplete basis matrix $\overline{\mathit{A}}=\left[\mathit{\alpha }\left({\theta }_1\right),\dots ,\mathit{\alpha }\left({\theta }_K\right)\right]$ and a sparse signal vector $\overline{\mathit{S}}={\left[s_1\left(n\right),\dots ,s_K\left(n\right)\right]}^{\mathrm{T}}$. In order to ensure that spatial sparsity holds, the number of potential signals $K$ must be much larger than the number of actual signals $L$, satisfying $K\gg L$. At this time, Equation 10 can be further expressed by Equation 11:

$\mathit{Y}\left(n\right)={\mathit{B}}^{\mathrm{T}}\overline{\mathit{A}}\overline{\mathit{S}}\left(n\right)+\mathit{N}\left(n\right)(11)$

At this time, $\mathit{B}\overline{\mathit{A}}$ is the measurement matrix in the theory of compressed sensing. When the spatial division angle is determined, the measurement matrix $\mathit{B}\overline{\mathit{A}}$ is determined accordingly. Compressed sensing is used in the field of DOA estimation, which is the process of reconstructing the original incident signal $S$ by receiving the signal $Y$ and the determined measurement matrix. Ideally, in the sparse signal vector $\overline{\mathit{S}}\left(n\right)$, only the element value signals of $L$ positions are not zero, and the element values ​​of the remaining $K-L$ positions are zero. At this time, the matrix composed of $L$ non-zero vectors is the original incident signal matrix, and the corresponding angular position is the DOA estimation value of the incident signal. This typical compressed sensing (CS) problem can be solved by solving the $l_1$-norm minimization model [21] which is defined by Equation 12:

$\widehat{\mathit{S}}=\arg {\min ‖\overline{\mathit{S}}\left(n\right)‖}_1$

$\mathrm{s}.\mathrm{t}.‖\mathit{Y}-\mathit{B}\overline{\mathit{A}}\overline{\mathit{S}}‖\le \delta (12)$

Among them, $\delta$ is a parameter related to noise, which indicates the error caused by noise.

3.3 ${\mathit{l}}_{\mathbf{1}}$-SVD algorithm based on asynchronous control

In the field of spatial signal processing, the sparsity of the signal in space must exist and the sparsity is high enough, so the compressed sensing theory can be applied to solve the DOA estimation value. The $l_1$-SVD algorithm is a classic $l_1$-norm minimization problem. It relies on singular value decomposition (SVD) of multi-snapshot data to achieve dimensionality reduction of the data matrix and solve the optimal solution of the target through a second-order convex programming method. First, the received signal matrix $\mathit{Y}$ is decomposed by singular value, which is expressed by Equation 13:

$\mathit{Y}=\mathit{U}\mathit{L}{\mathit{V}}^{\mathrm{T}}(13)$

Among them, $\mathit{U}$ and $\mathit{V}$ are both orthogonal matrices, respectively containing left singular vectors and right singular vectors, and $\mathit{L}$ is a diagonal matrix whose diagonal elements are singular values.

Define ${\mathit{D}}_K={\left[{\mathit{I}}_K\mathbf{0}\right]}^{\mathrm{T}}$, where ${\mathit{I}}_K$ represents the K-order identity matrix, 0 represents the $K\times \left(NN-K\right)$-dimensional zero matrix, and $NN$ represents the number of snapshots. Therefore, the received signal matrix can be reduced to an $M\times K$-dimensional matrix ${\mathit{Y}}_{SV}$ which is defined in Equation 14:

${\mathit{Y}}_{SV}=\mathit{U}\mathit{L}{\mathit{D}}_K=\mathit{Y}\mathit{V}{\mathit{D}}_K(14)$

By performing dimensionality reduction processing on the sparse signal vector Sˉ(n) and the noise vector $N\left(n\right)$ in turn, we can obtain the following Equations 15, 16:

${\overline{\mathit{S}}}_{SV}=\overline{\mathit{S}}\mathit{V}{\mathit{D}}_K(15)$

${\overline{\mathit{N}}}_{SV}=\mathit{N}\mathit{V}{\mathit{D}}_K(16)$

At this time, the received signal matrix after dimension reduction can be expressed by Equation 17:

${\mathit{Y}}_{SV}=\mathit{B}\mathit{A}{\overline{\mathit{S}}}_{SV}+N_{SV}(17)$

From Equation 17, we can see that the singular value decomposition processes the columns of the received signal matrix without changing the position of the non-zero rows. The dimension of the received signal matrix is ​​reduced from $M\times NN$ to $M\times L$, and the number of matrix columns is changed from the number of aliasing to the number of signal sources. In actual situations, $NN\gg L$ can greatly reduce the computational complexity. At the same time, the optimization model is transformed into Equation 18:

$\widehat{\mathit{S}}=\arg {\min ‖{\overline{\mathit{S}}}_{SV}‖}_1$

$\mathrm{s}.\mathrm{t}.{‖{\mathit{Y}}_{SV}-\mathit{B}\overline{\mathit{A}}{\overline{\mathit{S}}}_{SV}‖}_{\mathrm{F}}^2\le {\sigma }^2(18)$

Among them, ${‖·‖}_{\mathrm{F}}^2$ represents the Frobenius norm and ${\sigma }^2$ represents the Gaussian white noise power.

For the convenience of solving, the constrained optimization problem is transformed into an unconstrained optimization problem through the penalty function method expressed by Equation 19:

$\min {‖{\mathit{Y}}_{SV}-\mathit{B}\overline{\mathit{A}}{\overline{\mathit{S}}}_{SV}‖}_{\mathrm{F}}^2+\lambda {‖S^{l_2}‖}_1(19)$

Among them, $\lambda$ is the regularization parameter, which is related to the size of the noise. The selection of $\lambda$ is extremely important [22]. When the parameter is too small, a large number of pseudo peaks will appear in the spatial spectrum function. When the parameter is too large, the peak corresponding to the actual incoming wave direction will disappear.

4 Simulation results

This section evaluates the feasibility of estimating DOA by the asynchronous control method proposed in this paper through numerical simulation. Then, by comparing with the existing methods, the effectiveness of the asynchronous control method proposed in this paper is proved. All simulation experiments were carried out 1,000 times of Monte Carlo simulation. In the simulation experiment, this paper sets the incident signal bandwidth to 5 MHz and the control frequency of the TVM unit to 20 MHz. In the synchronous control method, up to four different states can be achieved in one code element period. In the asynchronous control method, in order to achieve the estimation performance of eight columns of TVM, let $∆t=1/2f_v$, and the simulation related parameters are shown in Table 1. The root-mean square error (RMSE) expressed is defined by Equation 20:

Table 1

Table 1. Simulation parameters.

$\text{RMSE}=\sqrt{\frac{1}{LJ}{\displaystyle \sum _{j=1}^J}{\displaystyle \sum _{l=1}^L}{\left[{\widehat{\theta }}_l\left(j\right)-{\theta }_l\right]}^2}(20)$

Where $J$ is the number of Monte Carlo experiments, $L$ is the number of signal sources, and ${\widehat{\theta }}_l\left(j\right)$ is the $j$ th Monte Carlo estimate of the lth signal.

4.1 Algorithm effectiveness

The spatial spectra with incident directions of −25° and 20° and the spatial spectra with incident directions of -1° and 1° are shown in Figures 6, 7, respectively. It compares the spatial spectra using the MUSIC algorithm and the asynchronous control method proposed in this paper. As can be seen from Figures 6, 7, the proposed method not only has a sharper spatial spectrum peak, but also has lower side lobes generated by interference signals, and can also effectively distinguish signals with close intervals in the incident direction. The MUSIC algorithm can distinguish signals that are far apart, but cannot distinguish adjacent signals. Therefore, the proposed asynchronous control method has higher resolution.

Figure 6

Figure 6. Spatial spectrum comparison at incident directions of −25° and 20°.

Figure 7

Figure 7. Spatial spectrum comparison at incident directions of -1° and 1°.

The spatial spectra corresponding to seven and three incident signals are shown in Figures 8, 9 respectively. Under the simulation parameters in Table 1, the number of virtual channels N = 4. Figure 8 shows the spatial spectrum function corresponding to the incident angle uniformly distributed between −45° and 45° when the number of incident signals L = 7. It can be seen from Figure 8 that both the asynchronous control method and the multi-channel array can obtain correct DOA estimation results, but the spatial spectrum sidelobes of the asynchronous control method are lower and the anti-interference ability is stronger, which is due to the gain brought by the heterogeneous directional pattern reception. In particular, when the synchronous control method is used, the spatial spectrum cannot be obtained because the number of virtual multi-channels is less than the number of incident signals.

Figure 8

Figure 8. 7-way incident signal corresponding to the spatial spectrum.

Figure 9

Figure 9. 3-way incident signal corresponding to the spatial spectrum.

Figure 9 shows the corresponding spatial spectrum function when the incident direction is (−30°, −10°, 25°). It can be seen from Figure 9 that although the synchronous control method can obtain an effective spatial spectrum function, the estimation error is large and DOA estimation cannot be achieved. The spatial sidelobe power of the asynchronous control method and the synchronous control method is lower than that of the multi-channel array, which proves that the DOA estimation method based on TVM has better anti-interference performance.

4.2 DOA estimation accuracy

The asynchronous control method proposed in this paper is compared with the synchronous control method and the classic MUSIC and l₁-SVD algorithms in the multi-channel array. The change of DOA estimation RMSE with SNR under different methods is shown in Figure 10, where the number of aliases is 1,000 and the incoming wave direction is (−20°, 30°). As shown in Figure 10, the DOA estimation accuracy of the asynchronous control method proposed in this paper and the l₁-SVD algorithm implemented by the multi-channel array is basically the same. At the same time, the RMSE estimated by the synchronous control method is higher and the DOA estimation accuracy is poor, which proves that when the control rate is insufficient, the DOA estimation accuracy of the synchronous control method is low, but the asynchronous control method can effectively improve the DOA estimation accuracy and obtain the theoretically optimal DOA estimation result, which proves the effectiveness of the asynchronous control method proposed in this paper.

Figure 10

Figure 10. The variation of RMSE of DOA estimation with SNR under different methods.

The variation of DOA estimation RMSE with the number of sources under different methods is shown in Figure 11, where the number of aliases is 1,000 and the SNR is 0 dB. As shown in Figure 11, the synchronous control method can only obtain effective DOA estimation results when there are less than three sources, and the estimated RMSE is high and the DOA estimation accuracy is poor. When the number of sources exceeds 4, no effective DOA estimation results can be obtained. This is because the number of virtual multi-channels obtained by the synchronous control method is insufficient when the control rate is insufficient, and the full information of the array cannot be fully utilized. In addition, the DOA estimation accuracy of the asynchronous control method and the l₁-SVD algorithm implemented by the multi-channel array is basically the same, which proves the effectiveness of the asynchronous control method proposed in this paper.

Figure 11

Figure 11. The variation of RMSE of DOA estimation with the number of sources under different methods.

5 Conclusion

This paper proposes a novel time-varying metasurface based DOA estimation method using the idea of asynchronous control. When the TVM unit control rate is limited by materials, the start time of the change of different column units is staggered to obtain multiple different responses within a state duration. In view of the problem of insufficient control rate, the TVM unit state duration is fully utilized to increase the number of virtual multi-channels. Then, a new measurement matrix is ​​obtained according to the asynchronous control codeword, which is used to evenly divide the signal space and obtain a sparse representation of the received signal. Finally, the l₁-SVD algorithm based on asynchronous control is used to realize DOA estimation. The simulation results verify the effectiveness of the proposed method. Compared with the existing synchronous control method, the DOA estimation performance has been greatly improved and can approach the theoretical optimal DOA estimation result. The future work will be to consider THz frequency and other operating conditions.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

Author contributions

KA: Methodology, Writing – review and editing, Software, Supervision, Writing – original draft, Investigation, Visualization, Validation, Conceptualization, Data curation. NS: Writing – review and editing, Validation, Conceptualization, Investigation, Writing – original draft, Software, Formal Analysis, Methodology, Resources, Visualization, Project administration. JH: Supervision, Writing – review and editing, Investigation, Writing – original draft, Conceptualization, Data curation, Visualization, Resources, Validation, Project administration. JE-G: Writing – review and editing, Formal Analysis, Methodology, Validation, Supervision, Writing – original draft, Software, Resources, Conceptualization. IH: Project administration, Visualization, Funding acquisition, Conceptualization, Resources, Validation, Formal Analysis, Software, Writing – review and editing, Methodology, Supervision, Data curation, Writing – original draft, Investigation.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported by Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2025R66), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia.

Conflict of interest

Author IH was employed by Deep Tech As. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Wang T, Han J, Ma X, Liu H, Li L. Frequency-diverse MIMO metasurface antenna for computational imaging with aperture rotation technique. Front Mater (2023) 9:1–13. doi:10.3389/fmats.2022.1112339

CrossRef Full Text | Google Scholar

2. Venneri F, Costanzo S, Borgia A. Fractal metasurfaces and antennas: an overview of advanced applications in wireless communications. Appl Sci (2024) 14(7):1–17. doi:10.3390/app14072843

CrossRef Full Text | Google Scholar

3. Awan W, Hussain N, Park S, Kim N, “Intelligent metasurface based antenna with pattern and beam reconfigurability for internet of things applications,” Alexandria Eng J (2024) 92, pp. 50–62. doi:10.1016/j.aej.2024.02.034

CrossRef Full Text | Google Scholar

4. Alibakhshikenari M, Babaeian F, Virdee B, Aissa S, Azpilicueta L, See C, et al. A comprehensive survey on various decoupling mechanisms with focus on metamaterial and metasurface principles applicable to SAR and MIMO antenna systems. IEEE Access (2020) 8:192965–3004. doi:10.1109/access.2020.3032826

CrossRef Full Text | Google Scholar

5. Alibakhshikenari M, Virdee B, Azpilicueta L, Moghadasi M, Akinsolu M, See C, et al. A comprehensive survey of Metamaterial transmission-line based antennas: design, challenges, and applications. IEEE Access (2020) 8:144778–808. doi:10.1109/access.2020.3013698

CrossRef Full Text | Google Scholar

6. Althuwayb A, Alibakhshikenari M, Virdee B, Rashid N, Kaaniche K, Atitallah A, et al. Metasurface-inspired flexible wearable MIMO antenna array for wireless body area network applications and biomedical telemetry devices. IEEE Access (2022) 11:1039–56. doi:10.1109/access.2022.3233388

CrossRef Full Text | Google Scholar

7. Alibakhshikenari M, Ali E, Soruri M, Dalarsson M, Moghadasi M, Virdee B, et al. A comprehensive survey on antennas on-chip based on metamaterial, metasurface, and substrate integrated waveguide principles for millimeter-waves and terahertz integrated circuits and systems. IEEE Access (2022) 10:3668–92. doi:10.1109/access.2021.3140156

CrossRef Full Text | Google Scholar

8. Kenari MA. Printed planar patch antennas based on metamaterial. Int J Electronics Lett (2013) 2(1):37–42. doi:10.1080/21681724.2013.874042

CrossRef Full Text | Google Scholar

9. Alibakhsikenari M, Movahhedi M, Naderian H. A new miniature ultra wide band planar microstrip antenna based on the metamaterial transmission line. In: IEEE asia-pacific conference on applied electromagnetics (APACE) (2012). p. 293–7.

Google Scholar

10. Hou L, Jin L, Huang K, Xiao S, Lou Y, Chen Y. Beamspace spatial music DOA estimation method using dynamic metasurface antenna. Entropy (2025) 27(4):1–15. doi:10.3390/e27040335

CrossRef Full Text | Google Scholar

11. Yildirim H, Storrer L, Doncker P, Louveaux J, Horlin F. A multi-antenna super-resolution passive Wi-Fi radar algorithm: combined model order selection and parameter estimatoin. IET Radar, Sonar and Navigation (2022) 16(8):1376–87. doi:10.1049/rsn2.12267

CrossRef Full Text | Google Scholar

12. Jung Y, Jeon H, Lee S, Jung Y. Scalable ESPRIT processor for direction-of-arrival estimation of frequency modulation continuous wave radar. Electronics (2021) 10(6):1–15. doi:10.3390/electronics10060695

CrossRef Full Text | Google Scholar

13. Zhu Y, Zhang W, Yi H, Xu H. Enhanced root-music algorithm based on matrix reconstruction for frequency estimation. Sensors (2023) 23(8):1829–15. doi:10.3390/s23041829

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Chen G, Su X, He L, Guan D, Liu Z. Coherent signal DOA estimation method based on space-time-coding. Remote Sensing (2025) 17(2):1–13. doi:10.3390/rs17020218

CrossRef Full Text | Google Scholar

15. Yao L, Zhang R, Hu C, Wu W. Off-grid DOA estimation for metasurface antenna systems using sparse Bayesian learning. AEU – Int J Electronics Commun (2025) 190:155615–2. doi:10.1016/j.aeue.2024.155615

CrossRef Full Text | Google Scholar

16. Zhou C, Gu Y, Fan X, Shi Z, Mao G, Zhang YD. Direction-of-arrival estimation for coprime array via virtual array interpolation. IEEE Trans Signal Process (2018) 66(22):5956–71. doi:10.1109/tsp.2018.2872012

CrossRef Full Text | Google Scholar

17. Hu D, He S, Li S, Zhu W. A dynamic beam switching metasurface based on angular mode-hopping effect. Front Phys (2024) 12:1–13. doi:10.3389/fphy.2024.1392115

CrossRef Full Text | Google Scholar

18. Pertsch T, Xiao S, Majumdar A, Li G. Optical metasurfaces: fundamentals and applications. Photon Res (2023) 11(5):1–15. doi:10.1364/prj.487440

CrossRef Full Text | Google Scholar

19. Hou J, Zhang X, Guo Y, Zhang R, Guo M. Design of electromagnetic metasurface using two dimensional crystal nets. Scientific Rep (2023) 13(7248):7248–14. doi:10.1038/s41598-023-32660-y

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Yang W, Qin J, Long J, Yan W, Yang Y, Li C, et al. A self-biased non-reciprocal magnetic metasurface for bidirectional phase modulation. Nat Electronics (2023) 6:225–34. doi:10.1038/s41928-023-00936-w

CrossRef Full Text | Google Scholar

21. Hoang T, Fusco V, Abbasi M, Yurduseven O. Single-pixel polarimetric direction of arrival estimation using programmable coding metasurface aperture. Scientific Rep (2021) 11:23830–14. doi:10.1038/s41598-021-03228-5

PubMed Abstract | CrossRef Full Text | Google Scholar

22. He C, Cao A, Chen J, Liang X, Zhu W, Geng J, et al. Direction finding by time-modulated linear array. IEEE Trans Antennas Propagation (2018) 66(7):3642–52. doi:10.1109/tap.2018.2835164

CrossRef Full Text | Google Scholar

23. Chen X, Zhang L, Liu S, Cui TJ. Artificial neural network for direction-of-arrival estimation and secure wireless communications via space-time-coding digital metasurfaces. Adv Opt Mater (2022) 10(23):1–16. doi:10.1002/adom.202201900

CrossRef Full Text | Google Scholar

24. Dai J, Tang W, Wang M, Chen MZ, Cheng Q, Jin S, et al. Simultaneous in situ direction finding and field manipulation based on space-time-coding digital metasurface. IEEE Trans Antennas Propagation (2022) 70(6):4774–83. doi:10.1109/tap.2022.3145445

CrossRef Full Text | Google Scholar

25. Fang X, Li M, Han J, Ramaccia D, Toscano A, Bilotti F, et al. Accurate direction-of-arrival estimation method based on space-time modulated metasurface. IEEE Trans Antennas Propagation¸ (2022) 70(11):10951–64. doi:10.1109/tap.2022.3184556

CrossRef Full Text | Google Scholar

26. Fu H, Dai F, Hong L. Two-dimensional off-grid DOA estimation with metasurface aperture based on MMV sparse Bayesian learning. IEEE Trans Instrumentation Meas (2023) 72:1–18. doi:10.1109/tim.2023.3318716

CrossRef Full Text | Google Scholar

27. Huang M, Zheng B, Cai T, Li X, Liu J, Qian C, et al. Machine-learning-enabled metasurface direction of arrival estimation. Nanophotonics (2022) 11(9):1–18. doi:10.1515/nanoph-2021-0663

CrossRef Full Text | Google Scholar

28. Carlon J, Castellanos M, Heath R. Hierarchical codebook design with dynamic metasurface antennas for energy-efficient arrays. IEEE Trans Wireless Commun (2024) 23(10):14790–804. doi:10.1109/TWC.2024.3419107

CrossRef Full Text | Google Scholar

29. Qian C, Kaminer I, Chen H. A guidance to intelligent metamaterials and metamaterials intelligence. Nat Commun (2025) 16:1154–26. doi:10.1038/s41467-025-56122-3

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zuo M, Xie S, Zhang X, Yang M. DOA estimation based on weighted l₁-norm sparse signal representation for low SNR scenarios. Sensors (2021) 21(13):1–15. doi:10.3390/s21134614

CrossRef Full Text | Google Scholar