Analog-Domain Suppression of Strong Interference Using Hybrid Antenna Array

Abstract

The proliferation of wireless applications, the ever-increasing spectrum crowdedness, as well as cell densification makes the issue of interference increasingly severe in many emerging wireless applications. Most interference management/mitigation methods in the literature are problem-specific and require some cooperation/coordination between different radio frequency systems. Aiming to seek a more versatile solution to counteracting strong interference, we resort to the hybrid array of analog subarrays and suppress interference in the analog domain so as to greatly reduce the required quantization bits of the analog-to-digital converters and their power consumption. To this end, we design a real-time algorithm to steer nulls towards the interference directions and maintain flat in non-interference directions, solely using constant-modulus phase shifters. To ensure sufficient null depth for interference suppression, we also develop a two-stage method for accurately estimating interference directions. The proposed solution can be applicable to most (if not all) wireless systems as neither training/reference signal nor cooperation/coordination is required. Extensive simulations show that more than 65 dB of suppression can be achieved for 3 spatially resolvable interference signals yet with random directions.

1. Introduction

Due to the proliferation of wireless systems, the ever-increasing spectrum crowdedness, as well as the densification of cellular networks, radio frequency (RF) interference has reappeared to be a major challenge in many emerging wireless applications [1,2,3,4,5,6,7,8]. This, for instance, happens in wirelessly powered communications (WPC) [2], where the range difference between power and information transfers makes the information signal weaker than the power transfer signal by many orders of magnitude [3]. Another example where the interference is highly detrimental can be found in the convergence of radar sensing and wireless communications that has attracted extensive attention recently [4,5,6]. Given its high transmission power, a radar can cause strong interference to the nearby communication system [7]. A further example showing the destructive impact of the strong interference is in the increasingly popular full-duplex communications; refer to [8] for details.

Numerous interference management/mitigation methods have been proposed. The resource scheduling is applied for relieving the mutual interference in WPC [9,10] and also in the co-existence of radar and communications (CRC) [11]. As pointed out in [2], this method can reduce spectral efficiency and puts strict requirement on time synchronization between different RF systems. Moreover, this method relies on the information exchange between participant RF systems [12]. Another popular method of suppressing RF interference is the null-space projection-based transmission, as applied in CRC [13,14] and full-duplex communications [15]. In particular, this method designs the transmitting beamforming/precoding to project transmitted signals onto the null-space of the interference matrix between the involved RF systems. Despite the effectiveness of the above methods, they require specially designed training signals between the participating RF systems for estimating the aforementioned inference matrix. In the case of massive MIMO systems, the training overhead can be heavy. To this end, it is desirable if a receiver itself has the capability of suppressing strong interference signals without requiring lengthy training from nearby transmitters or RF sources.

Antenna-array-based receivers can suppress interference signals by performing adaptive beamforming [16,17]. The well-known optimum criteria for adaptive beamforming include the minimum variance distortion-less response (also known as Capon), the maximum signal-to-interference-plus-noise ratio, and the linear constraints minimum variance [18]. These adaptive beamforming criteria generally rely on the fully digitized signal samples and hence require a high-dynamic-range front end and ADC to be equipped for each antenna. This requirement further leads to the high power consumption and cost of the digital-array-based receiver [19]. Moreover, given the limited dynamic range, information-bearing signals can be severely corrupted by clipping or quantization noises [2,20]. This makes the conventional adaptive beamforming-based interference suppression not sufficiently effective due to the potential loss/distortion of the direction information in the expected signals.

In contrast to the fully digital array, the hybrid antenna array of analog subarrays, as illustrated in Figure 1, has been recognized as a cost-effective and energy-efficient solution to counteracting strong RF interference [21,22]. This is attributed to its two significant benefits. One is that the analog subarrays can suppress strong interference signals prior to frequency conversion and digitization, hence easing the dynamic range requirements and reducing the power consumption of RF chains [23]. The other benefit of using a hybrid array is that the spatial degrees of freedom (DoF) are preserved given multiple analog subarrays. However, to enjoy these benefits of hybrid arrays, we need an efficient design of the analog beamformer that flexibly steers nulls towards strong interference signals while maintaining a flat amplitude response in the spatial passband, i.e., non-interference directions [2,23,24]. Hereafter, we refer to such a beamformer as the analog interference-nullifying beamformer (AINB).

So far, only a few works have specifically considered the AINB design in hybrid arrays. In [2], the AINB is constructed as a Kronecker product of component vectors, and each vector is designed based on the truncated Fourier/Hardamard matrices. However, this design only focuses on forming nulls, which may lead to severe attenuation in the spatial passband. In [24], the AINB is designed by successively approximating the desired beamformer as a linear combination of implementable analog beamformers. This method requires a reference signal which, however, would be fully corrupted by the strong quantization noises in the presence of interference signals, and the reference signal may not always be available. In [23], the AINB is achieved by optimizing the RF input impedance with the phase shifters used for shifting the null to the interference direction. Since only a single null is steered, this design can be incompetent in the presence of multiple interference signals.

To address the issues in the above works, we design a novel AINB by minimizing the spatial responses towards the interference directions and approximately maximizing the beam flatness in the spatial passband. Such a design can quickly adapt to the dynamic interference environment and form multiple nulls as required. However, due to the constant modulus constraints on the phase shifters in an analog subarray, the AINB design problem is highly non-convex and non-trivial to be solved efficiently in real time. Moreover, to steer nulls towards interference signals, the accurate estimations of the interference directions are required. In this paper, we consider the case where the power difference between interference and information-bearing signals is greater than the dynamic range of a receiver. In such a case, the AGCs in RF chains will reduce the receiving gains to prevent signal clipping, as illustrated in Figure 1. As a result, information-bearing signals will have too large quantization noises, which makes it infeasible to estimate the directions of interference signals. The AINB sequentially designed can substantially suppress interference signals, triggering the AGCs to adjust for the reception of information-bearing signals. The key contributions and novelties of this paper are summarized as follows.

• We develop an efficient solver for designing the phase-only AINB under the framework of majorization–minimization (MM). In particular, we design the objective function in such a way that we are able to simplify it substantially based on the newly unveiled relation between the spatial responses in the interference directions and the spatial passband. Thanks to the simplification of the objective function, we then propose a low-complexity method of constructing the majorization function, where we manage to remove the need of computationally intensive eigenvalue decomposition (EVD) required in the conventional construction. We further derive an iterative solver for the AINB design, where a closed-form solution with low complexity is achieved in each iteration. In addition, the impact of the initial solution to the proposed solver is investigated, based on which a high-quality initialization for the solver is established;
• We develop a two-stage angle of arrival (AoA) estimation method based on the conventional ESPRIT (estimation of signal parameters via invariance techniques). A major innovation of the method is the design of subarray beamforming in the two stages. In particular, an omnidirectionally flat beam is produced at each subarray in the first stage, while in the second stage the beam is created towards each of the AoA estimates obtained previously. To the best of our knowledge, it has not been investigated in the literature to use specially optimized flat beams for improving the AoA estimation performance of ESPRIT in hybrid antenna arrays.
• We provide extensive simulation results to validate the effectiveness of the proposed designs. As for the AINB, we cannot find similar designs in the literature. Thus, we comprehensively evaluate and observe numerous performance metrics, including spatial responses, interference suppression capability, and convergence curves for designing AINBs over tens of thousands of independent trials. As for the AoA estimation method, we employ the state-of-the-art [25] as a benchmark for the reasons to be explained at the beginning of Section 5. Due to the proposed use of deliberately optimized flat beams, the AoA estimation performance is substantially improved over the prior art [25]. Moreover, thanks to the high accuracy of the proposed AoA estimation method, the proposed AINB design can efficiently steer deep nulls towards interference signals.

It should be pointed out that it is not the first effort to apply ESPRIT to the hybrid arrays. The state of the art, such as that in [25], does so by focusing on augmenting the signal subspace at the expense of estimation time. Our study was motivated to achieve a large dimension of the signal subspace in the hybrid array that is not limited to the number of RF chains. This, however, is not a severe issue in our case, since the number of strong interference signals that are required to be suppressed in the analog domain is generally small [23]. We also remark that the hybrid array, which is also called a massive MIMO arrays, has been extensively studied in the context of the millimeter-wave (mmWave) and also terahertz communications for 5G and beyond [26,27,28,29,30]. Given the sparse nature of the mmWave channels, the AoA estimation methods recently developed for hybrid arrays, e.g., [31,32,33], mainly target at the LoS-dominated scenarios. In contrast, our design here assumes several strong interference signals.

The remainder of the paper is organized as follows. Section 2 formulates the AINB design problem, along with the signal model established. Section 3 elaborates the MM-based algorithm for designing AINB, with the objective simplification, algorithm development, and initialization detailed in three subsections. Section 4 illustrates the proposed two-stage AoA estimation method. The simulation results are provided in Section 5, followed by Section 6, which concludes the paper.

2. Problem Formulation

As shown in Figure 1, we consider a uniform linear array of M analog subarrays. Each subarray, having N antennas, performs analog beamforming using phase shifters. The spacing between any adjacent antennas is identical, as denoted by d. Note that the antenna coupling can be very small when $d\ge \lambda /2$ [21], where $\lambda$ denotes the wavelength. Thus, to focus on introducing the proposed interference suppression scheme, we assume $d\ge \lambda /2$ and ignore antenna coupling in this work. In the normal receiving array, the signals, impinging on the antennas, are first combined through subarray beamforming in the analog domain; then, each subarray output is processed by an individual radio frequency chain, including amplifying, down-converting in frequency, and filtering; and, finally, the baseband (or intermediate frequency, IF) signal is sampled and digitized through an analog-to-digital converter (ADC). A central digital signal processor collects the digitized subarray outputs for the sequential processing. In the presence of strong interference signals, AINB needs to be designed and performed across analog subarrays so that low-bit ADC can be used for better power efficiency, as illustrated in Section 1.

Considering the far-field receiving, the subarrays “see” interference signals from the same direction. Thus, the same AINB can be applied for all the subarrays. Letting $\mathbf{w}$ collect the beamforming weights of the AINB, we have:

$$\begin{array}{c}\mathbf{w}={[e^{\mathsf{j}{\varphi }_0},\cdots ,e^{\mathsf{j}{\varphi }_n},\cdots ,e^{\mathsf{j}{\varphi }_{N-1}}]}^{\mathrm{T}},\end{array}$$

(1)

where ${\varphi }_n$ is the value of the phase shifter connected to antenna n. As commonly performed in the beamformer design, we divide the angular region of [−90°, 90°] evenly into the following set:

$$\begin{array}{c}\mathsf{\Theta }=\{-{90}^{\circ },-{90}^{\circ }+{\delta }_{\theta }\cdots ,0^{\circ },\cdots ,{90}^{\circ }\},\end{array}$$

(2)

where ${\delta }_{\theta }$ is the minimum angular interval. Then, we design $\mathbf{w}$ to make the spatial responses at the discrete angles approach the desired spatial responses. Let ${\mathsf{\Theta }}_{\mathrm{j}}$ denote the set of angles in $\mathsf{\Theta }$ that are closest to interference AoAs and ${\mathsf{\Theta }}_{\mathrm{s}}=\mathsf{\Theta }/{\mathsf{\Theta }}_{\mathrm{j}}$ denote the set difference between $\mathsf{\Theta }$ and ${\mathsf{\Theta }}_{\mathrm{j}}$. For notation clarity, we add subscripts “${(·)}_{\mathrm{j}}$” and “${(·)}_{\mathrm{s}}$” for the variables related to interference and useful signals, respectively.

Steering deep nulls towards the directions of interference signals can be fulfilled by the following optimization:

$$\begin{array}{c}f\left(\mathbf{w}\right)\approxeq \mathbf{0}\hspace{1em}\mathrm{s}.\mathrm{t}.f\left(\mathbf{w}\right)={\mathbf{A}}_{\mathrm{j}}^{\mathrm{H}}\mathbf{w};{\mathbf{A}}_{\mathrm{j}}=[{\mathbf{a}}_{\mathrm{j}}\left({\left[{\mathsf{\Theta }}_{\mathrm{j}}\right]}_0\right),\cdots ,{\mathbf{a}}_{\mathrm{j}}\left({\left[{\mathsf{\Theta }}_{\mathrm{j}}\right]}_{L_{\mathrm{j}}-1}\right)],\end{array}$$

(3)

where $f\left(\mathbf{w}\right)$ collects the beamforming gains achieved by $\mathbf{w}$ towards the directions in ${\mathsf{\Theta }}_{\mathrm{j}}$, ${\left[{\mathsf{\Theta }}_{\mathrm{j}}\right]}_{l_{\mathrm{j}}}$ takes the $l_{\mathrm{j}}$-th entry in the set ${\mathsf{\Theta }}_{\mathrm{j}}$, and $L_{\mathrm{j}}$ is the total number of angles in ${\mathsf{\Theta }}_{\mathrm{j}}$. The dimensions of zero/one vector and identity matrix can be readily deduced given the context and hence are not explicitly noted in this paper. The steering vector ${\mathbf{a}}_{\mathrm{j}}\left({\left[{\mathsf{\Theta }}_{\mathrm{j}}\right]}_{l_{\mathrm{j}}}\right)$ can be written as:

$$\begin{array}{c}\hfill {\mathbf{a}}_{\mathrm{j}}\left({\left[{\mathsf{\Theta }}_{\mathrm{j}}\right]}_{l_{\mathrm{j}}}\right)={[e^{-\mathsf{j}n\pi \sin {\left[{\mathsf{\Theta }}_{\mathrm{j}}\right]}_0},\cdots ,e^{-\mathsf{j}n\pi \sin {\left[{\mathsf{\Theta }}_{\mathrm{j}}\right]}_{L_{\mathrm{j}}-1}}]}^{\mathrm{T}},\end{array}$$

(4)

where the antenna spacing d has taken half the wavelength to simplify the phase terms. To achieve a flat beam response in the spatial passband, i.e., non-interference directions, we generally minimize the ${\ell }_{\infty }$-norm of the spatial responses. However, we use the ${\ell }_2$-norm here to make it tractable to jointly achieve the interference suppression and flat beam response in non-interference directions. This will be seen shortly in Section 3.1. In particular, we formulate the following problem:

$$\begin{array}{c}\underset{\mathbf{w}}{\max }\hspace{1em}g\left(\mathbf{w}\right)\hspace{1em}\hspace{1em}\mathrm{s}.\mathrm{t}.g\left(\mathbf{w}\right)=\frac{1}{L_{\mathrm{s}}}{\parallel {\mathbf{A}}_{\mathrm{s}}^{\mathrm{H}}\mathbf{w}\parallel }_2^2;{\mathbf{A}}_{\mathrm{s}}=[{\mathbf{a}}_{\mathrm{s}}\left({\left[{\mathsf{\Theta }}_{\mathrm{s}}\right]}_0\right),\cdots ,{\mathbf{a}}_{\mathrm{s}}\left({\left[{\mathsf{\Theta }}_{\mathrm{s}}\right]}_{L_{\mathrm{s}}-1}\right)],\hfill \end{array}$$

(5)

where $L_{\mathrm{s}}$ denotes the cardinality of ${\mathsf{\Theta }}_{\mathrm{s}}$, and ${\mathbf{a}}_{\mathrm{s}}\left({\left[{\mathsf{\Theta }}_{\mathrm{s}}\right]}_{l_{\mathrm{s}}}\right)$ can be obtained by replacing ${\left[{\mathsf{\Theta }}_{\mathrm{j}}\right]}_{l_{\mathrm{j}}}$ in (4) with ${\left[{\mathsf{\Theta }}_{\mathrm{s}}\right]}_{l_{\mathrm{s}}}$. Combining (3) and (5), the AINB design can be formulated into the following optimization problem:

$$\begin{array}{c}\hfill \underset{\mathbf{w}}{\min }h\left(\mathbf{w}\right),\hspace{1em}\mathrm{s}.\mathrm{t}.h\left(\mathbf{w}\right)=\frac{\frac{1}{L_{\mathrm{j}}}{\parallel f\left(\mathbf{w}\right)\parallel }_2^2}{g\left(\mathbf{w}\right)};\left|\mathbf{w}\right|=\mathbf{1},\end{array}$$

(6)

where the numerator is the normalized power of the spatial responses in ${\mathsf{\Theta }}_{\mathrm{j}}$, and $|·|$ takes the point-wise modulus. Due to the constant modulus constraint, problem (6) is highly non-convex and hence non-trivial to solve efficiently (in real time).

Using ${\ell }_2$-norm in (5) can greatly simplify the non-convex problem (6), as will be illustrated in Section 3.1. The simplification further enables us to develop a real-time solver for problem (6), as will be presented in Section 3.2. Moreover, we unveil in Section 3.3 that initializing the solver properly can help enhance the flatness of the AINB in the non-jamming angular region. The AINB design to be detailed is based on the accurate AoAs of interference signals. Thus, an accurate AoA estimation method will be developed in Section 4.

3. AINB Design

In this section, an efficient solver to problem (6) is proposed under the framework of MM. Specifically, we first simplify the objective function using signal processing techniques/properties, then develop MM-based solver, and finally propose a high-quality initialization for the solver.

3.1. Simplifying Beamformer Design Problem

We notice that there is a relationship between the numerator and denominator of the objective function given in (6). In particular, based on (3)–(5), we have:

$$\begin{array}{c}{\parallel f\left(\mathbf{w}\right)\parallel }_2^2+L_{\mathrm{s}}g\left(\mathbf{w}\right)={\parallel {\mathbf{A}}^{\mathrm{H}}\mathbf{w}\parallel }_2^2={\displaystyle \sum _{l=0,\cdots ,L-1}}{\left|\frac{1}{N}{\displaystyle \sum _{n=0}^{N-1}}{w}_n{e}^{-\mathsf{j}n\pi \sin {\theta }_l}\right|}^2\hfill \\ ={\displaystyle \sum _{l=0,\cdots ,L-1}}{\left|{\displaystyle \sum _{n=0}^{N-1}}{w}_n{e}^{-\mathsf{j}n{u}_l}\right|}^2,u_l=\pi \sin {\theta }_l\hfill \end{array}$$

where $\mathbf{A}=[{\mathbf{A}}_{\mathrm{j}},{\mathbf{A}}_{\mathrm{s}}]$, $w_n$ is the n-th entry of $\mathbf{w}$ given in (1), $L=\left|\mathsf{\Theta }\right|$ denotes the cardinality of $\mathsf{\Theta }$ given in (2), and ${\theta }_l$ is the l-th element in $\mathsf{\Theta }$. In the final result of the above derivation, ${\displaystyle {\sum }_{n=0}^{N-1}}{w}_n{e}^{-\mathsf{j}n{u}_l}$ can be seen as the sampling of the discrete time Fourier transform (DTFT) of $\mathbf{w}$ at $u_l$. When L is sufficiently large, $u_l$ approaches a continuous u in $[-\pi ,\pi ]$. Thus, provided a large L, we have:

$$\begin{array}{c}\left({\parallel f\left(\mathbf{w}\right)\parallel }_2^2+L_{\mathrm{s}}g\left(\mathbf{w}\right)\right)/L\approxeq \frac{1}{2\pi }\underset{{\int }_{-\pi }^{\pi }\mathrm{d}u}{\underbrace{\underset{L\to \infty }{\lim }\frac{2\pi }{L}\sum _{l=0,\cdots ,L-1}}}{\left|\sum _{n=0}^{N-1}{w}_n{e}^{-\mathsf{j}n{u}_l}\right|}^2\hfill \\ =\frac{1}{2\pi }{\int }_{-\pi }^{\pi }|\tilde{w}{\left(u\right)|}^2\mathrm{d}u\stackrel{\left(a\right)}{=}\sum _{n=0}^{N-1}{\left|w_n\right|}^2\stackrel{\left(b\right)}{=}N,\hfill \end{array}$$

(7)

where $\tilde{w}\left(u\right)$ denotes the DTFT of $w_n$, the equation $\stackrel{\left(a\right)}{=}$ is due to the Parseval’s theorem of DTFT [34], and the equation $\stackrel{\left(b\right)}{=}$ is a result of the constant-modulus constraint on $\mathbf{w}$; see (6).

The derivation in (7) further leads to:

$$g\left(\mathbf{w}\right)=\left(NL-{\parallel f\left(\mathbf{w}\right)\parallel }_2^2\right)/L_{\mathrm{s}}.$$

Substituting this into (6), the objective function becomes:

$$h\left(\mathbf{w}\right)=\frac{\frac{L_{\mathrm{s}}}{L_{\mathrm{j}}}{\parallel f\left(\mathbf{w}\right)\parallel }_2^2}{NL-{\parallel f\left(\mathbf{w}\right)\parallel }_2^2},$$

which is a monotonically increasing function of ${\parallel f\left(\mathbf{w}\right)\parallel }_2^2$. To this end, $h\left(\mathbf{w}\right)$ is minimized when ${\parallel f\left(\mathbf{w}\right)\parallel }_2^2$ takes the minimum. The above derivation and analysis are formally summarized into the following lemma.

Lemma 1. 

Provided a large L, the problem (6) is approximately equivalent to the following simplified minimization:

$$\begin{array}{c}\underset{\mathbf{w}}{\min }{\parallel f\left(\mathbf{w}\right)\parallel }_2^2={\mathbf{w}}^{\mathrm{H}}{\mathbf{R}}_{\mathrm{j}}\mathbf{w},\hspace{1em}\mathrm{s}.\mathrm{t}.{\mathbf{R}}_{\mathrm{j}}={\mathbf{A}}_{\mathrm{j}}{\mathbf{A}}_{\mathrm{j}}^{\mathrm{H}};\left|\mathbf{w}\right|=\mathbf{1},\hfill \end{array}$$

(8)

where ${\mathbf{A}}_{\mathrm{j}}$ is given in (3).

3.2. MM-Based Iterative Solver for Problem (8)

Under the framework of the MM technique, problem (8) can be iteratively solved by: majorizing $f\left(\mathbf{w}\right)$ at ${\mathbf{w}}_i$ (the solution at iteration i), leading to the majorization function denoted by $\tilde{f}(\mathbf{w},{\mathbf{w}}_i)$ and minimizing $\tilde{f}(\mathbf{w},{\mathbf{w}}_i)$ subject to the constraint in (8). We see from (8) that the objective function is in a quadratic form. According to [35], a quadratic function ${\mathbf{x}}^{\mathrm{H}}\mathbf{L}\mathbf{x}$ can be majorized at any ${\mathbf{x}}_0$ by the following function:

$$\begin{array}{c}{\mathbf{x}}^{\mathrm{H}}\mathbf{M}\mathbf{x}+2\Re \left\{{\mathbf{x}}_0^{\mathrm{H}}(\mathbf{L}-\mathbf{M})\mathbf{x}\right\}+{\mathbf{x}}_0^{\mathrm{H}}(\mathbf{M}-\mathbf{L}){\mathbf{x}}_0,\end{array}$$

(9)

where $\mathbf{L}$ and $\mathbf{M}⪰\mathbf{L}$ are Hermitian matrices. Therefore, to majorize the objective function of (8), we need to find a Hermitian matrix ${\tilde{\mathbf{R}}}_{\mathrm{j}}$ such that ${\tilde{\mathbf{R}}}_{\mathrm{j}}⪰{\mathbf{R}}_{\mathrm{j}}$. A straightforward selection of ${\tilde{\mathbf{R}}}_{\mathrm{j}}$ is ${\lambda }_{\max }\left\{{\mathbf{R}}_{\mathrm{j}}\right\}\mathbf{I}$, where ${\lambda }_{\max }\left\{{\mathbf{R}}_{\mathrm{j}}\right\}$ denotes the maximum eigenvalue. Such a selection, however, needs the eigenvalue decomposition of the N-dimensional matrix ${\mathbf{R}}_{\mathrm{j}}$. This can be avoided applying the following result.

Lemma 2. 

The maximum eigenvalue of ${\mathbf{R}}_{\mathrm{j}}$ is upper bounded by:

$$\begin{array}{c}{\lambda }_{\max }\left\{{\mathbf{R}}_{\mathrm{j}}\right\}\le \overline{\lambda }=N{L}_{\mathrm{j}},\end{array}$$

(10)

where $L_{\mathrm{j}}$ is the cardinality of Θ given in (3).

Proof. 

The maximum eigenvalue of a Hermitian matrix can be expressed as the following Rayleigh quotient [Theorem 7.16] [36]:

$$\begin{array}{c}{\lambda }_{\max }\left\{{\mathbf{R}}_{\mathrm{j}}\right\}=\underset{{\parallel \mathbf{x}\parallel }_2^2=1}{\max }\frac{{\mathbf{x}}^{\mathrm{H}}{\mathbf{R}}_{\mathrm{j}}\mathbf{x}}{{\mathbf{x}}^{\mathrm{H}}\mathbf{x}}.\hfill \end{array}$$

(11)

Substituting ${\mathbf{R}}_{\mathrm{j}}={\mathbf{A}}_{\mathrm{j}}{\mathbf{A}}_{\mathrm{j}}^{\mathrm{H}}$ into the above definition, we have:

$$\begin{array}{c}{\lambda }_{\max }\left\{{\mathbf{R}}_{\mathrm{j}}\right\}=\underset{{\parallel \mathbf{x}\parallel }_2^2=1}{\max }\sum _{\begin{array}{c}{l}_{\mathrm{j}}=0,\\ \cdots ,L_{\mathrm{j}}-1\end{array}}{\left|{\mathbf{x}}^{\mathrm{H}}{\mathbf{a}}_{\mathrm{j}}\left({\left[{\mathsf{\Theta }}_{\mathrm{j}}\right]}_{l_{\mathrm{j}}}\right)\right|}^2\hfill \\ \le \sum _{\begin{array}{c}{l}_{\mathrm{j}}=0,\\ \cdots ,L_{\mathrm{j}}-1\end{array}}\left({\mathbf{x}}^{\mathrm{H}}\mathbf{x}\right)\times \left({\mathbf{a}}_{\mathrm{j}}^{\mathrm{H}}\left({\left[{\mathsf{\Theta }}_{\mathrm{j}}\right]}_{l_{\mathrm{j}}}\right){\mathbf{a}}_{\mathrm{j}}\left({\left[{\mathsf{\Theta }}_{\mathrm{j}}\right]}_{l_{\mathrm{j}}}\right)\right)=N{L}_{\mathrm{j}},\hfill \end{array}$$

(12)

where “≤” is based on the Cauchy–Schwartz inequality ([36], Sec. 1.4).    □

Substituting $\mathbf{x}=\mathbf{w}$, $\mathbf{M}={\mathbf{R}}_{\mathrm{j}}={\mathbf{A}}_{\mathrm{j}}{\mathbf{A}}_{\mathrm{j}}^{\mathrm{H}}$, ${\mathbf{x}}_0={\mathbf{w}}_i$ and $\mathbf{L}=\overline{\lambda }\mathbf{I}$ into (9), the objective function of (8) is majorized at a known point ${\mathbf{w}}_i\in {\mathbb{C}}^{N\times 1}$ by the following function:

$$\begin{array}{c}\tilde{f}(\mathbf{w},{\mathbf{w}}_i)=-2\Re \left\{{\mathbf{r}}_i^{\mathrm{H}}\mathbf{w}\right\},\mathrm{s}.\mathrm{t}.{\mathbf{r}}_i={(\overline{\lambda }\mathbf{I}-{\mathbf{A}}_{\mathrm{j}}{\mathbf{A}}_{\mathrm{j}}^{\mathrm{H}})}^{\mathrm{H}}{\mathbf{w}}_i\end{array}$$

(13)

where $\overline{\lambda }$ is given in (10), and the terms independent of $\mathbf{w}$ are dropped for brevity. Taking ${\mathbf{w}}_i$ as the solution to problem (8) in the i-th iteration, $\tilde{f}(\mathbf{w},{\mathbf{w}}_i)$ is minimized when the real-taking operator achieves the maximum. Due to the constant modulus constraint on $\mathbf{w}$ given in (8), the maximum is achieved when the phases of $\mathbf{w}$ are aligned with those of ${\mathbf{r}}_i$ in a pointwise manner. To this end, we achieve the following efficient solver for problem (8).

Proposition 1. 

The iterative solver to problem (8) is given by:

$$\begin{array}{c}\hfill {\mathbf{w}}_{i+1}=e^{\mathsf{j}\arg \left\{{\mathbf{r}}_i\right\}},\mathrm{s}.\mathrm{t}.{\mathbf{r}}_i={\left(\frac{L_{\mathrm{j}}}{N}\mathbf{I}-{\mathbf{A}}_{\mathrm{j}}{\mathbf{A}}_{\mathrm{j}}^{\mathrm{H}}\right)}^{\mathrm{H}}{\mathbf{w}}_i,\end{array}$$

(14)

where ${\mathbf{w}}_i$ is the solution in iteration i, $e^{(·)}$, and $\arg \{·\}$ are pointwise operators, and $\overline{\lambda }$ is derived in Lemma 2.

Note that the iterative solver has a low computational complexity in each iteration. In particular, the complexity is in the order of $\mathcal{O}\left(N^2\right)$, which is dominated by computing ${\mathbf{r}}_i$ given in (14). This owes to the simplification of the objective function achieved in Lemma 1 and the construction of the majorization based on the upper bound derived in Lemma 2.

3.3. Initializing ${\mathbf{w}}_0$

Given the non-convex nature of problem (8), the iterative solution achieved in Proposition 1 tends to converge to a local minimum. Namely, the convergence performance of (14) is closely related to the initialization, i.e., ${\mathbf{w}}_0$. Below, we first investigate the impact of ${\mathbf{w}}_0$ on the quality of AINB that is measured by the null depth and beam flatness in the spatial passband, and accordingly design a high-quality ${\mathbf{w}}_0$. The investigation is performed through a set of simulations.

Figure 2 observes the value of the objective function given in (8) under the iterative solution given in (14), where $N=16$, ${\delta }_{\theta }={0.1}^{\circ }$, ${\mathsf{\Theta }}_{\mathrm{j}}=\{{30.5}^{\circ },{60.9}^{\circ },-{50.3}^{\circ }\}$, and $I=800$ denotes the total iteration number. Three independent trials are performed, where the initial beamforming weight vector, i.e., ${\mathbf{w}}_0$, is independently and randomly generated for each trial. We see that all the three trials have monotonically decreasing value of the objective function before convergence. This validates the effectiveness of the solver developed in Proposition 1. From Figure 2, we also see that the convergence speed is affected by the initialization of ${\mathbf{w}}_0$. In particular, the number of iterations required for convergence increases from about 180 to 300 when comparing the second trial with the third one.

Figure 2 also observes the features of the optimized beams in the spatial passband, where $10{\log }_{10}\{{\mathbf{w}}_i^{\mathrm{H}}{\mathbf{R}}_{\mathrm{s}}{\mathbf{w}}_i/L_{\mathrm{s}}\}$ denotes the average power of the beam steered by ${\mathbf{w}}_i$ in the spatial passband, while $20{\log }_{10}\left\{\frac{\max \left\{\right|{\mathbf{w}}_i^{\mathrm{H}}{\mathbf{A}}_{\mathrm{s}}\left|\right\}}{\min \left\{\right|{\mathbf{w}}_i^{\mathrm{H}}{\mathbf{A}}_{\mathrm{s}}\left|\right\}}\right\}$ depicts the flatness of the beam, where ${\mathbf{R}}_{\mathrm{s}}={\mathbf{A}}_{\mathrm{s}}{\mathbf{A}}_{\mathrm{s}}^{\mathrm{H}}$ with ${\mathbf{A}}_{\mathrm{s}}$ given in (5). Note that the larger $20{\log }_{10}\left\{\frac{\max \left\{\right|{\mathbf{w}}_i^{\mathrm{H}}{\mathbf{A}}_{\mathrm{s}}\left|\right\}}{\min \left\{\right|{\mathbf{w}}_i^{\mathrm{H}}{\mathbf{A}}_{\mathrm{s}}\left|\right\}}\right\}$ is, the less flat the beam becomes. We see from Figure 2 that the average powers of the optimized beams in the spatial passband are almost identical. This validates the proposed simplification of the objective function for AINB design, as illustrated in Lemma 1. In contrast to the average power, we see from Figure 2 that the flatness of the optimized beams in the spatial passband can substantially vary given different initializations of ${\mathbf{w}}_0$. In particular, the difference of the beam flatness is larger than 10 dB between trials 2 and 3.

Figure 3 plots the amplitude responses of the beams steered by ${\mathbf{w}}_I$’s obtained from the three independent trials. We see that three deep nulls are produced by the beams obtained in the three trials, and the null depth is below $-300$ dB. This validates the strong interference suppression ability achieved by the beamformer design problem formulated in Lemma 1 and the iterative solver derived in Proposition 1. We also see from the figure (particular the zoomed-in sub-figure in the middle) that the passband variations in the obtained beams can be subatantially different. This is consistent with what has been observed in Figure 2.

Figure 4 observes the impact of the initialization on the convergence performance of the proposed solver (14). In particular, the left y-axis gives the amplitude variance, i.e., the beam flatness, of the beam steered by ${\mathbf{w}}_I$ with respect to that of the beam steered by ${\mathbf{w}}_0$, while the converging value of the objective function is shown on the right y-axis. We see that similar converging values of the objective function keep steady regardless of the initialized ${\mathbf{w}}_0$. We also see that the flatness of the beam steered by the final ${\mathbf{w}}_I$ increases linearly overall with that of the beam steered by the initial ${\mathbf{w}}_0$.

Given the observations made above, we conclude that an initialized ${\mathbf{w}}_0$ leading to a flat amplitude response in the spatial passband is better in the sense that the optimized beam can also have better flatness in the spatial passband in addition to the strong ability of interference suppression. On the other hand, we notice that a flat beam under the constant modulus constraint given in (8) is difficult to realize. The spatial response of the beam steered by $\mathbf{w}$ is given by:

$$\begin{array}{c}\hfill P\left(u_l\right)=\sum _{n=0}^{N-1}{w}_n{e}^{-\mathsf{j}n{u}_l},l=0,1,\cdots ,L-1,\end{array}$$

(15)

where $w_n$ is the n-th entry of $\mathbf{w}$, $u_l=\pi \sin {\theta }_l$, and ${\theta }_l$ is the l-th element in $\mathsf{\Theta }$ given in (2). Regarding $u_l\in [-\pi ,\pi ]$ as the angular frequency, we obtain that $P\left(u_l\right)$ is the (discrete) Fourier transform of $w_n$. According to the uncertainty principle of Fourier transform [34], a signal spreading out in one domain corresponds to a signal localized (such as the Dirac delta function) in the dual domain. In our case, the beam flatness implies that $P\left(u_l\right)$ spreads out across $u_l$, and hence, $w_n$ is expected to be a localized signal sequence. This, however, contradicts with the constant modulus constraint on $\mathbf{w}$, i.e., $|w_n|=1\forall n$. To steer a beam as flat as possible in the spatial passband, we formulate the following problem to reduce the maximum amplitude variation:

$$\begin{array}{c}\hfill {\mathbf{w}}_0^{*}:{\mathrm{argmin}}_{{\mathbf{w}}_0}\frac{\max \left\{\right|{\mathbf{w}}_0^{\mathrm{H}}\mathbf{A}\left|\right\}}{\min \left\{\right|{\mathbf{w}}_0^{\mathrm{H}}\mathbf{A}\left|\right\}},\mathrm{s}.\mathrm{t}.\left|{\mathbf{w}}_0\right|=1,\end{array}$$

(16)

where $\mathbf{A}$ is the union of ${\mathbf{A}}_{\mathrm{j}}$ and ${\mathbf{A}}_{\mathrm{s}}$ given in (3) and (5), respectively. Since the initial ${\mathbf{w}}_0$ is designed offline, we can resort to the bio-inspired optimization tools, e.g., the genetic algorithm (GA) and the currently popular grey wolf optimizer [37], etc., to solve the highly non-convex problem (16).

3.4. Speeding up Convergence

As seen from Figure 2, directly solving (14) can lead to a number of iterations in the order of ${10}^2$. To improve the convergence speed, we employ the widely used squared expectation maximization (SQUAREM) [38] to solve (14) and summarize the overall iterative processing in Algorithm 1. SQUAREM is a gradient descent method in essence. It performs Cauchy and Barzilai–Borwein (BB) methods in succession in each iteration. Both Cauchy and BB methods perform gradient descent using different step sizes; refer to [38] for details. This leads to the update of $\mathbf{w}$, as enclosed in the projector $\mathcal{P}(·)$; see Steps 9 and 12 of Algorithm 1. The BB step size [39] is adopted in SQUAREM, which results in the update of $\alpha$, as given in Step 8. Updating $\alpha$ uses the intermediate outcomes from Steps 6 and 7. Steps 11 and 12 perform backtracking, a widely used trick in adjusting step size [40], to keep the objective non-increasing. Finally, Algorithm 1 terminates when the condition in either Step 3 or 15 is reached.

Algorithm 1 MM-based Analog Beamformer Design.

1: Input: $\mathsf{\Theta }$, ${\mathsf{\Theta }}_{\mathrm{j}}$, N and ${\mathbf{w}}_0^{*}$; 2: Initialize the iteration index $i=0$ and set $ϵ$, such as ${10}^{-6}$, for example; ▹$ϵ$ is used for examining convergence. 3: while$i\le I_{\max }$do                   ▹$I_{\max }$ is the maximum iteration number. 4:     Solving (14) based on ${\mathbf{w}}_i$ gives ${\mathbf{w}}_{i+1}^{\left(1\right)}$; 5:     Solving (14) again based on ${\mathbf{w}}_{i+1}^{\left(1\right)}$ gives ${\mathbf{w}}_{i+1}^{\left(2\right)}$; 6:     ${\mathbf{v}}_1={\mathbf{w}}_{i+1}^{\left(1\right)}-{\mathbf{w}}_i$; 7:     ${\mathbf{v}}_2={\mathbf{w}}_{i+1}^{\left(2\right)}-{\mathbf{w}}_{i+1}^{\left(1\right)}-{\mathbf{v}}_1$; 8:     $\alpha =-\frac{\parallel {\mathbf{v}}_1\parallel }{\parallel {\mathbf{v}}_2\parallel }$; 9:     $\mathbf{w}=\mathcal{P}({\mathbf{w}}_i-2\alpha {\mathbf{v}}_1+{\alpha }^2{\mathbf{v}}_2)$;                    ▹$\mathcal{P}\left(\mathbf{x}\right)=\frac{\mathbf{x}}{\left|\mathbf{x}\right|}$ (point-wise). 10:     while ${\parallel f\left(\mathbf{w}\right)\parallel }_2^2>{\parallel f\left({\mathbf{w}}_i\right)\parallel }_2^2$ do 11:         Update $\alpha =(\alpha -1)/2$; 12:         $\mathbf{w}=\mathcal{P}({\mathbf{w}}_i-2\alpha {\mathbf{v}}_1+{\alpha }^2{\mathbf{v}}_2)$; 13:     end while 14:     Obtain ${\mathbf{w}}_{i+1}=\mathbf{w}$ and $i=i+1$; 15:     Terminate iteration if $\parallel f\left({\mathbf{w}}_{i+1}\right){\parallel }_2^2-{\parallel f\left({\mathbf{w}}_i\right)\parallel }_2^2<ϵ$; 16: end while 17: Return $\mathbf{w}={\mathbf{w}}_i$.

Algorithm 1 is developed based on the MM framework, which converges to a finite value for the problem such as (8) with the quadratic objective bounded by zero. Following the convergence analysis in ([41], Section V), one can readily show that solving (14) iteratively for problem (8) certainly converges to a stationary point. However, due to the use of SQUAREM, the convergence of Algorithm 1 can only be achieved when the algorithm starts in the vicinity of a stationary point already [42]. As will be shown in Section 5, enabled by the proposed initialization method, Algorithm 1 converges within a few tens of iterations with satisfactory AINBs achieved.

4. Estimation of Interference AoAs

As illustrated in Section 3.3, the uncertainty principle determines that the constant-modulus $\mathbf{w}$ can only yield a highly localized AINB. This is beneficial in the sense that a deep null can be achieved, as shown in Figure 3. On the other hand, this limits the null width to be as small as possible. The limitation of null width further requires the AoA estimates of the interference signals to be highly accurate. In this section, a two-stage AoA estimation method is developed for the hybrid antenna arrays illustrated in Figure 1. At the core of the method is the conventional ESPRIT. A key problem of applying ESPRIT to the hybrid antenna array is how to design subarray beamforming. This is solved below.

Assume that there are P strong interference signals to be suppressed. Denote the p-th interference signal as $s_p\left(k\right)$ and its AoA as ${\theta }_p$, where k denotes the index of discrete time, also known as snapshot. We consider a short time period during which ${\theta }_p\forall p$ does not change. The signals impinging on the N antennas of subarray m can be collected by the following vector:

$$\begin{array}{c}\hfill {\mathbf{x}}_m\left(k\right)=\sum _{p=0}^{P-1}{s}_p\left(k\right)e^{\mathsf{j}mN{u}_p}{\mathbf{a}}_N\left(u_p\right),u_p=\pi \sin {\theta }_p,\end{array}$$

(17)

where $e^{\mathsf{j}mN{u}_p}$ is caused by the phase offset between the m-th and the first (i.e., reference) subarray, and ${\mathbf{a}}_N\left({\theta }_p\right)$ is an $N\times 1$ steering vector. Specifically, ${\mathbf{a}}_N\left({\theta }_p\right)$ is written as:

$$\begin{array}{c}\hfill {\mathbf{a}}_N\left(u_p\right)={\left[1,e^{\mathsf{j}{u}_p},\cdots ,e^{\mathsf{j}\pi (N-1)u_p}\right]}^{\mathrm{T}}.\end{array}$$

(18)

Let ${\mathbf{w}}_m\left(k\right)$ denote the analog beamforming weight vector of subarray m at the k-th snapshot. The subarray output is the inner product between ${\mathbf{w}}_m\left(k\right)$ and ${\mathbf{x}}_m$, as given by:

$$\begin{array}{c}\hfill y_m\left(k\right)=\sum _{p=0}^{P-1}{s}_p\left(k\right)e^{\mathsf{j}mN{u}_p}{\mathbf{w}}_m^{\mathrm{H}}\left(k\right){\mathbf{a}}_N\left(u_p\right)+{\xi }_m\left(k\right),\end{array}$$

(19)

where ${\xi }_m\left(k\right)$ is an AWGN. Stacking the M subarray outputs into a vector gives:

$$\begin{array}{cc}\hfill & \tilde{\mathbf{y}}\left(k\right)=\sum _{p=0}^{P-1}{\mathsf{\Lambda }}_p{\mathbf{W}}^{\mathrm{H}}\left(k\right){\mathbf{a}}_N\left(u_p\right)s_p\left(k\right)+\mathbf{z}\left(k\right),\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {\mathsf{\Lambda }}_p=\mathrm{diag}\left\{{\mathbf{a}}_M\left(N{u}_p\right)\right\},\mathbf{W}\left(k\right)=\left[\cdots ,{\mathbf{w}}_m\left(k\right),\cdots \right];\hfill \\ \hfill & {\mathbf{a}}_M\left(N{u}_p\right)={\left[\cdots ,e^{\mathsf{j}mN{u}_p},\cdots \right]}^{\mathrm{T}};\hfill \\ \hfill & \mathbf{z}\left(k\right)={\left[\cdots ,{\xi }_m\left(k\right),\cdots \right]}^{\mathrm{T}},\hfill \end{array}$$

(20)

where $\mathrm{diag}\{·\}$ generates a diagonal matrix. Provided that the same beamformer is applied across the M subarrays over K snapshots, i.e.:

$$\begin{array}{c}\hfill {\mathbf{w}}_m\left(k\right)={\mathbf{w}}_{m^{\prime }}\left(k^{\prime }\right)=\mathbf{w}(\forall m\ne m^{\prime },\forall k\ne k^{\prime }),\end{array}$$

(21)

we can simplify $\tilde{\mathbf{y}}\left(k\right)$ into:

$$\begin{array}{cc}\hfill \mathbf{y}\left(k\right)& =\sum _{p=0}^{P-1}{\mathbf{a}}_M\left(N{u}_p\right)g_p{s}_p\left(k\right)+\mathbf{z}\left(k\right)=\mathbf{A}\mathbf{s}\left(k\right)+\mathbf{z}\left(k\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{A}=\left[\cdots ,{\mathbf{a}}_M\left(N{u}_p\right),\cdots \right];\mathbf{s}\left(k\right)={\left[\cdots ,g_p{s}_p\left(k\right),\cdots \right]}^{\mathrm{T}},\hfill \end{array}$$

(22)

where $g_p={\mathbf{w}}^{\mathrm{H}}{\mathbf{a}}_N\left(u_p\right)$ is the subarray beamforming gain on the p-th interference signal.

From (22), we see that ESPRIT is now applicable provided that (C1) $M>P$ and (C2) the P interference signals are not coherent. For (C1), we remark that only those interference signals that are strong enough to cause signal clipping are of interest, and the number of such interference signals can be limited. If $P>M$, we can trade the time-domain degree of freedom (DoF) for the spatial ones, as designed in [25]. As for the second condition (C2), we remark that techniques such as spatial smoothing [18] can be employed to de-correlate incident signals, as derived in (22). Given the above reasons, we proceed to assume that both conditions are satisfied and focus on the selection of $\mathbf{w}$ for accurate AoA estimation.

From the signal model derived in (22), we see that the subarray beamforming introduces a coefficient $g_p$ to the p-th interference signal. To ensure high SNRs for the AoA estimation, we expect that $g_p$, if not enhancing, does not attenuate the incident signal too much. However, in the initial stage, we do not have any a priori information on the AoAs of interference signals. Thus, an omnidirectionally flat beam is suggested for a first-stage AoA estimation. Such a beam has been designed in Section 3.3, where the analog beamformer is initialized for the proposed AINB design. In particular, the beamformer is represented by ${\mathbf{w}}_0^{*}$ that is optimized via solving (16) (offline).

By taking $\mathbf{w}={\mathbf{w}}_0^{*}$ in (21), a first stage of AoA estimation is performed, as summarized in Algorithm 2, where the well-developed ESPRIT is encapsulated as a function based on ([18], Section 9.3). Note that performing ESPRIT based on subarray outputs obtained in (22) only estimates $N{u}_p$. According to (17), $u_p$ has the range of $[-\pi ,\pi ]$, given ${\theta }_p\in [-{90}^{\circ },{90}^{\circ }]$. Thus, we have $N{u}_p\in [-N\pi ,N\pi ]$. Since the angle-taking in Step 25 of Algorithm 2 only returns the angle between $\pm \pi$. Thus, there is an estimation ambiguity in ${\widehat{Nu}}_p$, which leads to the following N possible estimates of $u_p$:

$$\begin{array}{c}\hfill {\widehat{u}}_p\left(n\right)=\left({\widehat{Nu}}_p+2n\pi \right)/N,n=0,1,\cdots ,N-1.\end{array}$$

(23)

Algorithm 2 An Accurate Two-Stage AoA Estimation Method.

1: procedureStage 1(Initial Estimation) 2:     Take $\mathbf{w}={\mathbf{w}}_0^{*}$ and perform ${\widehat{Nu}}_p=$ ESPRIT$\left({\mathbf{w}}_0^{*},\mathbf{y}\left(k\right),K\right)$; 3:     Extract ${\widehat{u}}_p\left(n\right)$ from ${\widehat{Nu}}_p$, as performed in (23); 4:     The final estimate is ${\widehat{u}}_p^{\left(1\right)}={\widehat{u}}_p\left(n^{*}\right)$ with $n^{*}$ given in (24); 5: end procedure 6: procedureStage 2(AoA Refinement) 7:     for all ${\widehat{u}}_p^{\left(1\right)}(p=0,1,\cdots )$ do 8:         Take $\mathbf{w}={\mathbf{a}}_N\left({\widehat{u}}_p^{\left(1\right)}\right)$; see (18) and (21); 9:         Perform ${\widehat{Nu}}_p^{\left(2\right)}=$ ESPRIT$\left({\mathbf{a}}_N\left({\widehat{u}}_p^{\left(1\right)}\right),\mathbf{y}\left(k\right),K\right)$; 10:         Get ${\widehat{u}}_p^{\left(2\right)}\left(n\right)(n=0,\cdots ,N-1)$ from ${\widehat{Nu}}_p^{\left(2\right)}$; see (23); 11:         The final estimate is: ${\widehat{u}}_p={\min }_{{\widehat{u}}_p^{\left(2\right)}\left(n\right)}|{\widehat{u}}_p^{\left(2\right)}\left(n\right)-{\widehat{u}}_p^{\left(1\right)}|$. 12:     end for 13: end procedure 14: function${\widehat{Nu}}_p=$ ESPRIT$\left(\mathbf{w},\mathbf{y}\left(k\right),K\right)$ 15:     Perform subarray beamforming using $\mathbf{w}$, as performed in (21); 16:     Compute $\mathbf{R}\left(=\frac{1}{K}{\sum }_{k=0}^{K-1}\mathbf{y}\left(k\right){\mathbf{y}}^{\mathrm{H}}\left(k\right)\right)$; see (22) for $\mathbf{y}\left(k\right)$; 17:     Estimate the number of paths based on $\mathbf{R}$, yielding $\widehat{P}$; 18:     The EVD of $\mathbf{R}$ gives $\mathbf{R}{\lambda }_i={\mathbf{v}}_i{\lambda }_i$ with ${\lambda }_i>{\lambda }_{i+1}$; 19:     Construct ${\mathbf{U}}_{\mathrm{s}}=\left[{\mathbf{v}}_0,\cdots ,{\mathbf{v}}_{\widehat{P}-1}\right]$; 20:     Denote ${\mathbf{U}}_1={\left[{\mathbf{U}}_{\mathrm{s}}\right]}_{0:M-2,:}$ and ${\mathbf{U}}_2={\left[{\mathbf{U}}_{\mathrm{s}}\right]}_{2:M-1,:}$; 21:     Compute $\mathbf{C}=\left[\begin{array}{c}{\mathbf{U}}_1^{\mathrm{H}}\\ {\mathbf{U}}_2^{\mathrm{H}}\end{array}\right]\left[{\mathbf{U}}_1,{\mathbf{U}}_2\right]$. 22:     The EVD of $\mathbf{C}$ leads to $\mathbf{C}=\left[\begin{array}{c}{\mathbf{V}}_{11}{\mathbf{V}}_{12}\\ {\mathbf{V}}_{21}{\mathbf{V}}_{22}\end{array}\right]\mathsf{\Lambda }\left[\begin{array}{c}{\mathbf{V}}_{11}^{\mathrm{H}}{\mathbf{V}}_{12}^{\mathrm{H}}\\ {\mathbf{V}}_{21}^{\mathrm{H}}{\mathbf{V}}_{22}^{\mathrm{H}}\end{array}\right]$; 23:     Compute $\mathsf{\Phi }=-{\mathbf{V}}_{12}{\mathbf{V}}_{22}^{-1}$; 24:     Find the eigenvalues of $\mathsf{\Phi }$: ${\lambda }_p(p=0,\cdots ,\widehat{P}-1)$; 25:     Return ${\widehat{Nu}}_p=\arg \left\{{\lambda }_p\right\}$.                ▹$\arg \left\{\right\}$ takes angle. 26: end function

This ambiguity can be suppressed by enumerating the possible estimates and identifying the one leading to the largest receiving power as the final estimate [25]. In particular, the following problem is formulated for removing the ambiguity:

$$\begin{array}{c}\hfill n^{*}:\underset{n}{\max }\frac{{\sum }_{k=0}^{\tilde{K}}{\left|y_m\left(k\right)\right|}^2}{\tilde{K}},\mathrm{s}.\mathrm{t}.{\mathbf{w}}_m\left(n\right)={\mathbf{a}}_N\left({\widehat{u}}_p\left(n\right)\right),\end{array}$$

(24)

where $y_m\left(k\right)$ is obtained by plugging the above ${\mathbf{w}}_m\left(n\right)$ into (19), while ${\mathbf{w}}_m\left(n\right)$ is obtained by taking $u_p={\widehat{u}}_p\left(n\right)$ in (18). For ease of exposition, we use the m-th subarray in (24) to illustrate the idea of removing the estimation ambiguity in (23). Given M subarrays, we can simultaneously test M estimates. Being straightforward, the details are suppressed. Let the final estimate of the first stage be denoted by ${\widehat{u}}_p^{\left(1\right)}$. As performed in Steps 7–12 of Algorithm 2, the subarrays point towards each of ${\widehat{u}}_p^{\left(1\right)}$’s, and ESPRIT is performed again with the a greater beamforming gain exploited for better estimation performance, as compared with the first stage.

Next, we remark that the computational complexity of Algorithm 2 is in the order of $\mathcal{O}\left(K{M}^2+8{(M-1)}^3\right)$, which is generally low due to the small value of M. Note that the first part $\mathcal{O}\left(K{M}^2\right)$ is for computing $\mathbf{R}$ in Step 16. Normally, we do not include this part in evaluating the complexity of ESPRIT. Here, due to the small value of M, $\mathcal{O}\left(K{M}^2\right)$ can be comparable to or even larger than the second part, i.e., $\mathcal{O}\left(8{(M-1)}^3\right)$. This second part is for the EVD in Step 22 and is an upper limit obtained based on the maximum value of $\widehat{P}(=M-1)$. The complexity of the other EVD in Step 18 is in the order of $\mathcal{O}\left(M^3\right)$, which is generally lower than $\mathcal{O}\left(8{(M-1)}^3\right)$ for small M’s. Thus, the later is used in the overall computational complexity.

5. Simulation Results

In this section, the simulation results are provided to validate the proposed designs. Unless otherwise specified, the following parameters are used: $N=16$, $M=4$, $K=30$, $u_p\sim {\mathcal{U}}_{[-\pi ,\pi ]}$, $s_p\sim \mathcal{C}\mathcal{N}(0,{\sigma }_{\mathrm{i}}^2)$, and ${\sigma }_{\mathrm{n}}^2=0.01$ Watt, where $s_p$ denotes the p-th $(p=0,1,\cdots ,P-1)$ interference signal, as given in (17), and ${\sigma }_{\mathrm{n}}^2$ is the power of the AWGN in the subarray output, i.e., ${\xi }_m\left(k\right)$ given in (19). When applicable, $N_{\mathrm{trial}}=2\times {10}^4$ independent trials are performed to obtain an averaged performance. When designing AINB, the whole angular region given in (2) is considered, and the angle step ${\delta }_{\theta }$ takes ${0.1}^{\circ }$. Moreover, to evaluate the proposed AoA estimation method, the state-of-the-art hybrid ESPRIT (H-ESPRIT) [25] is simulated for comparison. Note that H-ESPRIT runs Stage 1 of Algorithm 2 by taking $\mathbf{w}={\mathbf{w}}_{\mathrm{r}}$, where ${\mathbf{w}}_{\mathrm{r}}$ denotes the subarray beamforming vector with randomly generated phases. The random subarray beamforming used by H-ESPRIT can handle the case that the number of subarrays is smaller than that of incident sources, while the proposed algorithm cannot. However, from the open literature, H-ESPRIT can be most related to the proposed Algorithm 2 in the sense that they are both for hybrid antenna arrays and can estimate the whole angular region. Thus, we employ H-ESPRIT as a benchmark.

For reproducibility, we provide in (25) the omnidirectionally flat beams for $N=8,16$ and 24, respectively.

These beams are obtained by solving problem (16) using the GA toolbox in MATLAB [43]. They are used in the following simulations. Figure 5 plots the amplitude response of the beams steered by ${\mathbf{w}}_0^{*}$ given above, in comparison with the beams steered by ${\mathbf{w}}_{\mathrm{r}}$ with randomly generated phases in $[0,2\pi ]$. We see that the optimized ${\mathbf{w}}_0^{*}$ has substantially better beam flatness than ${\mathbf{w}}_{\mathrm{r}}$.

$$\begin{array}{c}{\mathbf{w}}_0^{*}=\exp \left[\mathrm{j}\left(0.2988,1.1812,2.6686,1.3361,0.0484,\right.\right.\left.\left.-3.0455,-0.5214,3.1145\right)\right];\\ {\mathbf{w}}_0^{*}=\exp \left[\mathrm{j}\left(2.9280,-1.4330,-2.7477,0.3015,-0.4253,\right.\right.2.1168,0.9070,1.5538,-1.3083,1.5321,1.8232,\cdots \\ \left.\left.1.3972,-0.0852,0.4551,-2.0871,-2.5000\right)\right];\mathrm{and}\\ {\mathbf{w}}_0^{*}=\exp \left[\mathrm{j}\left(2.4481,-2.6471,2.2376,-1.6242,-2.4315,\right.\right.-2.7266,-0.5203,-2.6060,-1.1905,3.0591,\cdots \\ 0.5462,-1.8140,2.7608,1.9835,1.6974,-1.5860,0.3268,-1.8279,-1.2872,1.1163,\cdots \\ \left.\left.2.4170,-2.5172,-0.5398,-1.0528\right)\right].\end{array}$$

(25)

Figure 6 plots the MSE of AoA estimates as the INR, the ratio between ${\sigma }_{\mathrm{i}}^2$ and ${\sigma }_{\mathrm{n}}^2$, increases. A single source is considered in this simulation to focus on illustrating the benefits of using a flat beam for subarrays and adding a second stage of AoA refinement. We see that the AoA estimates from the first stage of the proposed Algorithm 2, using ${\mathbf{w}}_0^{*}$ given in (25), improve over H-ESPRIT [25]. In particular, to achieve the same estimation accuracy, the proposed method reduces the SNR requirement over 10 dB, compared with H-ESPRIT. This manifests the importance of steering an omnidirectionally flat beam for the accurate AoA estimation. We also see from the figure that the AoA estimates from the second stage are further improved over those from the first stage. This is owed to the greater beamforming gain achieved by the directional beam used in the second stage of Algorithm 2.

Figure 7 plots the CDF of the interference suppression after applying the proposed AINB at a subarray. The AINB is designed by running the proposed Algorithm 1, where the AoA estimates obtained for plotting Figure 6 are used for constructing ${\mathsf{\Theta }}_{\mathrm{j}}$ required by the algorithm. We see that using the proposed AoA estimation method achieves better interference suppression, as compared with using H-ESPRIT [25]. Take ${\sigma }_{\mathrm{i}}^2=-10$ dB for instance. More than 95% of the independent trials suppress interference to lower than $-30$ dB when using the proposed method for AoA estimation, while only less than 80% of the trials can achieve the equivalent interference suppression when using H-ESPRIT.

Next, we consider the case of multiple interference signals. Figure 8 (upper) plots the CDF of the sum of the absolute AoA estimation errors of multiple signals. We see that the proposed Algorithm 2 achieves better estimation performance than H-ESPRIT. For the proposed method, we see that the second-stage AoA estimation substantially outperforms the first stage. In particular, the absolute estimation error after the second stage is almost always smaller than ${0.0531}^{\circ }$. This is because the mutual interference between different incident signals can be smaller in the second stage that use a directional beam, as compared with the first stage that uses omnidirectional beams illustrated in Figure 5.

Figure 8 (lower) illustrates the interference suppression ability achieved under different sets of AoA estimates from the upper sub-figure. We see that a better AoA estimation results in a stronger interference suppression. This is reasonable, since the estimates used for the upper sub-figure are passed to Algorithm 2 for designing the AINBs. It is worth pointing out that jointly running Algorithms 1 and 2, the interference suppression is always greater than $-65$ dB. This not only validates the high performance of the proposed AoA estimation and AINB design but also demonstrates the feasibility of using hybrid antenna arrays to counteract strong interference signals. In addition, from Figure 8, we can observe the impact of angle quantization error (which is resembled by the angle estimation error) on interference suppression.

Figure 9 illustrates the amplitude responses of the AINBs designed using Algorithm 1, where the target and array parameters set for Figure 2 are used here, and ${10}^3$ trials are carried out with independently and randomly generated interference signals and noises. From Figure 9 (top), we see that the absolute AoA estimation errors of all the three sources are consistently smaller than $0.05$ over ${10}^3$ independent trials. Thus, the AoA estimations do not have dependencies over the trials. However, in each trial, the obtained AoA estimates are used for designing the AINB by performing the proposed Algorithm 1. The spatial amplitude responses of all AINBs are plotted in Figure 9 (bottom). We see that, due to the consistently accurate AoA estimates, the AINBs obtained in the ${10}^3$ trials overlap.

Moreover, we see from Figure 9 (bottom) that the nulls deeper than $-200$ dB are achieved in the directions of three interference signals in all the ${10}^3$ independent trials. Moreover, we see that the beam flatness in the spatial passband is much better than that achieved in Figure 3. This validates the high quality of the initial solution (for solver (14)), as designed in Section 3.3. Lastly, given the overlapping beams obtained over ${10}^3$ trials, we conclude that the proposed Algorithm 1 is reliable in terms of the convergence performance.

Figure 10 plots the convergence curves of performing Algorithm 1 when designing AINBs for Figure 7 and Figure 9. From Figure 10a, we see that, for single-source cases, the proposed Algorithm 1 can converge within 10 iterations. Note that the difference among the convergence curves is caused by the randomly changing AoAs over independent trials in Figure 7. In Figure 10b, we see that, when there are 3 interference signals to nullify, Algorithm 1 converges within 15 iterations. We also see that the convergence curves overlap with each other. This is owed to the high accuracy of the AoA estimates illustrated in Figure 9. Combining all the simulation results presented above, the proposed interference suppression scheme, consisting of Algorithms 1 and 2, can quickly and effectively suppress strong interference signals using hybrid antenna arrays.

Figure 11 illustrates the impact of the quantization bit of the phase shifters in analog subarrays on the proposed design, where the same $2\times {10}^4$ independent trials from Figure 8 are rerun here. In the method development, we employ the continuous phase shifters, while this can only be approximated in practice. For a phase shifter with B bits, it can take $2^B$ states, i.e., $e^{\mathsf{j}\frac{2\pi b}{B}},b=0,1,\cdots ,2^B-1$. To count for the phase shifter quantization, we replace each entry of ${\mathbf{w}}^{\star }$, which is the optimal result of Algorithm 1, with the closest value the phase shifter can take. From Figure 11, we see that the number of quantization bits of the phase shifters has a substantial impact on the interference suppression ability. As the number of quantization bits increases, the optimal ${\mathbf{w}}^{\star }$ can be better approximated, hence also leading to stronger interference suppression. From Figure 11, we also see that, even for the small 5-bit quantization, the proposed design can achieve greater than 38 dB interference suppression for all $2\times {10}^4$ independent trials (with random interference directions over trials). Though the interference signals may not be thoroughly suppressed, the power difference between interference and information-bearing signals can be substantially reduced. Consequently, AGCs will readjust, allowing the reception of information-bearing signals without being clipped. This further enables the suppression digital-domain interference, such as the null projection reusing the estimated interference AoAs.

6. Conclusions and Future Works

In this paper, we develop a versatile solution for suppressing strong interference signals using the hybrid array of subarrays. This is achieved by a real-time solver for designing the phase-only AINB that steers deep nulls towards the interference directions and maintains flat in the spatial passband. This is also accomplished by an accurate AoA estimation method using the ESPRIT algorithm. In particular, to achieve the same estimation accuracy, the proposed method reduces the SNR requirement by more than 10 dB, compared with state-of-the-art H-ESPRIT [25]. However, we would point out that H-ESPRIT has the capability of augmenting signal subspace dimension, while the proposed method does not have this capability. Moreover, simulation results show that, employing a uniform linear array of four subarrays each with 16 antennas, the proposed solution can provide the interference suppression of 65 dB or higher in the presence of three resolvable interference signals with randomly distributed directions. In addition, the interference suppression larger than 200 dB is also observed when deliberately setting the three interference directions apart. Since the hybrid array has been extensively studied in numerous use cases of mmWave and THz communications, our design is promising to counteract strong interference in many scenarios.

In this work, we rely on the uncertainty principle of the Fourier transform and minimize the null depth to approximately maintain the beam flatness in the spatial passband. However, we notice from Figure 9 that the flatness in the spatial passband can be uneven. This may affect some tasks, e.g., radar/radio sensing, that have stringent requirements on signal strength. As a future work, we will consider to impose the equi-ripple constraints on the beam response in the spatial passband. This is in light of the Parks–McClellan filter design [34], suggesting that the equi-ripple constraint can help increase the minimum mainlobe level. However, due to the non-convexity of the design problem, such a constraint may affect the efficiency and convergence speed/performance of a solver.