Optimization Approaches in Sideband Calculations and a Non- Iterative Harmonic Suppression Strategy in 4D Arrays

Abstract

In this chapter, the harmonic calculations and suppression methods in 4D arrays are explored. First, the fundamentals of the 4D arrays including the switching strategies and metaheuristic approach to suppression of sidebands are provided, and further, the harmonic reduction strategies using optimizers are discussed. Finally, a novel noniterative algorithmic way to suppress the harmonic radiations is introduced and exemplified by an illustrative example.

1. Introduction

Channel capacity is one of the most critical constraints in wireless communication technology that both academic and industrial foundations are trying to efficiently use and/or to increase the capacity with the help of various multiplexing, duplexing and coding algorithms. Although software-based algorithms allow channel capacity to be efficiently used, when it is compared with the hardware modifications of antenna elements (e.g. adaptive smart antenna, MIMO, switched arrays, etc.), the gain stays at minor levels; however, it does not mean that either one of them is strictly selected. Instead, hardware- and software-based modifications should be mixed in order to obtain the best results [1, 2].

In the late 1950s, Shanks and Bickmore introduced a new concept of antenna array design which is based on switching of individual array elements periodically [3]. The concept so-called time modulation enables time variable, t, as an antenna design parameter alongside of complex excitations and placement of the elements. While having a new design parameter, due to the nature of periodic excitations, infinite number of sideband radiations occurs at the multiplying frequencies of switching frequency. In a short while after the introduction of the time modulation concept, Kummer et al. analyzed the time-modulated arrays (i.e., TMAs) in terms of sideband radiations (i.e., SR) [4]. However, time modulation concept has not attracted much attention in scientific community, until the publication of a study on SR suppression with adjusting the sidelobe level using differential evolution algorithm by Yang et al. in 2002 [5]. Furthermore, with the emerging technology in computers and mobile communication infrastructures, TMAs gained popularity, and many studies have been conducted using different metaheuristic algorithms to suppress the harmonics [6, 7, 8, 9, 10, 11, 12] and using different excitation schemes so as to have more control parameters on array design such as shifting/splitting pulses and using different pulse amplitudes [6, 13, 14, 15, 16, 17]. Additionally, interference suppression and adaptive beamforming using TMAs have been studied in the literature [18, 19, 20, 21]. Moreover, TMAs having different geometries such as planar and conformal are analyzed for harmonic suppression [22, 23, 24], and closed form finite calculations have been done expressed in order to calculate total radiated power in SRs for different excitation schemes and geometries [25, 26, 27, 28, 29, 30, 31, 32].

Although there are many studies that are trying to suppress the all SRs since harmonic radiations consume main radiation power, the emerging concepts such as direction finding, spatial diversity and adaptive beamforming require harmonic radiations to be used in communication. In the literature, TMAs have been used for direction finding and direction of arrival estimation applications [33, 34, 35, 36, 37], harmonic beam steering, adaptive beamforming and communication over sidebands [38, 39, 40]. Furthermore, an overview about TMAs and additional enhancements for TMA concept have been published by Maneiro-Catoria et al. in 2017 [41, 42].

In this chapter, harmonic radiation suppression methods using metaheuristic algorithms will be examined in detail, total power calculations for SRs will be presented and novel non-iterative algorithmic method for harmonic suppression will be introduced.

2. Theoretical background of switched arrays

Antenna arrays having identical elements are examined independently of directional characteristics of the radiating elements (e.g., dipole, horn, etc.) since total electric field is calculated by multiplication of array factor with electric field of single element at reference point according to pattern multiplication. The time-modulation concept introduces a periodic switching for array elements and array factor, naturally, becomes dependent to excitation function unlike conventional arrays. Although it is shown that signal radiation in switched arrays (i.e., time-modulated arrays, TMAs) may be approximated to an unmodulated array under the condition of switching frequency being much smaller than carrier frequency by Bregains et al. [25], relation between conventional arrays’ and switched arrays’ radiation power density should be derived in order to obtain limitations and constraints of calculations [12, 25, 43].

In order to derive an average power density for TMAs, it is convenient to start the derivation from the Poynting vector, and the instantaneous power density is defined as

$$\mathcal{W}=\mathcal{E}\times \mathcal{H}$$ E1

where $\mathcal{E}$ and $\mathcal{H}$ represent the time-dependent electric and magnetic field vectors, and once the time modulation is considered, electric and magnetic fields may be expressed by

$$\mathcal{E}=\mathit{Re}\left\{\mathit{E}{e}^{i\left({\omega }_0+q{\omega }_p\right)t}\right\}=\frac{1}{2}\left[\mathit{E}{e}^{i\left({\omega }_0+q{\omega }_p\right)t}+{\mathit{E}}^{\mathbf{\ast }}{e}^{i\left({\omega }_0+q{\omega }_p\right)t}\right]$$ E2

and

$$\mathcal{H}=\mathit{Re}\left\{\mathit{H}{e}^{i\left({\omega }_0+q{\omega }_p\right)t}\right\}=\frac{1}{2}\left[\mathit{H}{e}^{i\left({\omega }_0+q{\omega }_p\right)t}+{\mathit{H}}^{\mathbf{\ast }}{e}^{i\left({\omega }_0+q{\omega }_p\right)t}\right].$$ E3

Here, $\mathit{Re}\left\{z\right\}$ and $q\in \left(-\infty ,\infty \right)$ stand for real part of $z\in \mathbb{C}$ and harmonic level, respectively. Furthermore, ${\omega }_0=2\pi f_0$ and ${\omega }_p=2\pi f_p$ represent angular carrier frequency and angular switching frequency, respectively. If Eqs. (2) and (3) are substituted into Eq. (1) and after some basic mathematical operations, Eq. (1) may be rewritten as

$$\mathcal{W}=\frac{1}{2}\mathit{Re}\left\{\mathit{E}\times {\mathit{H}}^{\mathbf{\ast }}\right\}+\frac{1}{2}\mathit{Re}\left\{\mathit{E}\times \mathit{H}{e}^{i2\left({\omega }_0+q{\omega }_p\right)t}\right\}.$$ E4

In order to obtain time average radiation power density, Poynting vector can be integrated over one period (i.e., $T_p$ ), and the time average Poynting vector may be expressed as

$${\mathit{W}}_{⟨\mathit{av}⟩}=\frac{1}{2T_p}{\int }_{T_p}\mathit{Re}\left\{\mathit{E}\times {\mathit{H}}^{\mathbf{\ast }}\right\}\mathit{dt}+\frac{1}{2T_p}{\int }_{T_p}\mathit{Re}\left\{\mathit{E}\times \mathit{H}{e}^{i2\left({\omega }_0+q{\omega }_p\right)t}\right\}\mathit{dt}=\frac{1}{2}{Q}_1+\frac{1}{2}{Q}_2.$$ E5

If $Q_1$ is integrated, it may be written as

$$Q_1=\frac{1}{T_p}\mathit{Re}\left\{\mathit{E}\times {\mathit{H}}^{\mathbf{\ast }}\right\}{\int }_{T_p}\mathit{dt}=\mathit{Re}\left\{\mathit{E}\times {\mathit{H}}^{\mathbf{\ast }}\right\}.$$ E6

Once the $Q_2$ is examined, as the term $\mathit{E}\times \mathit{H}$ is time-independent complex constant (i.e., $\mathit{E}\times \mathit{H}=\mathrm{g},\mathrm{g}\in \mathbb{C}$ ), it becomes

$$Q_2=\frac{1}{T_p}\mathit{Re}\left\{\mathit{E}\times \mathit{H}\right\}{\int }_{T_p}\text{cos}\left(2\left({\omega }_0+q{\omega }_p\right)t\right)\mathit{dt}-\frac{1}{T_p}\mathit{Im}\left\{\mathit{E}\times \mathit{H}\right\}{\int }_{T_p}\text{sin}\left(2\left({\omega }_0+q{\omega }_p\right)t\right)={\mathrm{\Gamma }}_1-{\mathrm{\Gamma }}_2$$ E7

where $\mathit{Im}\left\{z\right\}$ denotes imaginer part of $z\in \mathbb{C}$ . If ${\mathrm{\Gamma }}_1$ and ${\mathrm{\Gamma }}_2$ are solved, Eqs. (6) and (7) may be rewritten as

$${\mathrm{\Gamma }}_1=\mathit{Re}\left\{\mathit{E}\times \mathit{H}\right\}\frac{\text{sin}\left(2\left({\omega }_0+q{\omega }_p\right)t\right)}{2T_p{\omega }_p\left(\frac{{\omega }_0}{{\omega }_p}+q\right)}$$ E8

and

$${\mathrm{\Gamma }}_2=\mathit{Im}\left\{\mathit{E}\times \mathit{H}\right\}\left[\frac{\text{cos}\left(2\left({\omega }_0+q{\omega }_p\right)t\right)-1}{2T_p{\omega }_p\left(\frac{{\omega }_0}{{\omega }_p}+q\right)}\right]$$ E9

where $T_p{\omega }_p=2\pi$ . Looking closely to Eq. (8), ${\mathrm{\Gamma }}_1$ tends to zero as $\kappa ={\omega }_0/{\omega }_p$ ratio tends to infinity. The same situation holds true for also Eq. (9). Since ${\mathrm{\Gamma }}_2$ oscillates in range $\left[1/2\pi \left(\kappa +m\right):0\right]$ as $\kappa$ tends to infinity, ${\mathrm{\Gamma }}_2$ also vanishes. Hence, it can be said that both ${\mathrm{\Gamma }}_1$ and ${\mathrm{\Gamma }}_2$ vanish when ${\omega }_0\gg {\omega }_p$ so that under $\kappa \gg 1$ approximation Eq. (7) reduces to

$$Q_2={\mathrm{\Gamma }}_1-{\mathrm{\Gamma }}_2\approx 0$$ E10

Since, $Q_2$ vanishes, the time average Poynting vector may be written as

$${\mathit{W}}_{⟨\mathit{av}⟩}\approx \frac{1}{2}\mathit{Re}\left\{\mathit{E}\times {\mathit{H}}^{\mathbf{\ast }}\right\}.$$ E11

It must be noted that Eq. (11) is only valid under the condition of $\kappa ={\omega }_0/{\omega }_p\gg 1$ , and it clearly defines that the time average power density of switched arrays may be approximated to conventional arrays’ under defined conditions. In the light of this information, the array factor in time domain of an $N$ element time-modulated array with periodic switching function of $g_n\left(t\right)$ may be written as

$$\mu \left(\theta ,\varphi ,t\right)=e^{i{\omega }_{0}t}{\sum }_n{W}_n{g}_n\left(t\right)e^{\mathit{ik}\left({\widehat{a}}_{\mathit{pn}}·{\widehat{a}}_r\right)}$$ E12

where $W_n$ , ${\widehat{a}}_{\mathit{pn}}$ and ${\widehat{a}}_r$ represent the complex excitations, the position vector and the radial unit vector in Cartesian coordinates, respectively. Since $g_n\left(t\right)$ is periodic in time, it can be decomposed into Fourier series, hence, if $G_{\mathit{nq}}$ is being the Fourier coefficients, the array factor becomes

$$\mu \left(\theta ,\varphi ,t\right)=e^{i{\omega }_{0}t}{\sum }_n{W}_n\sum _{q=-\infty }^{\infty }{G}_{\mathit{nq}}{e}^{\mathit{iq}{\omega }_{p}t}{e}^{\mathit{ik}\left({\widehat{a}}_{\mathit{pn}}·{\widehat{a}}_r\right)}.$$ E13

Hence, in the light of the result obtained in Eq. (11), this array factor may be written in phasor form for a specific harmonic number $q$ as

$${\mu }_q\left(\theta ,\varphi \right)={\sum }_n{W}_n{G}_{\mathit{nq}}{e}^{\mathit{ik}\left({\widehat{a}}_{\mathit{pn}}·{\widehat{a}}_r\right)}$$ E14

As shown in Eq. (14), in switched arrays, the excitation strategies directly affect both the main and sideband radiations. Every periodic function is eligible to be used for switching the array elements; however, once the constraint number that shapes the radiation patterns increases, excitation strategies are needed to be evolved in order to meet the requirements. Time modulation concept firstly introduced by binary switching scheme, which switches the elements between on (‘1’) and off (‘0’) (i.e., variable aperture size (VAS) time scheme), then, the ideas shifting (i.e., pulse shift (PS)) and splitting pulses (i.e., pulse split (PSp) time scheme) are raised in order to get more control on antenna array characteristics. Additionally, another scheme so-called variable pulse amplitude (VPA) has been recently introduced, but the definition of this scheme is left to the further sections. An illustrative representation of these schemes is depicted in Figure 1.

Figure 1.

Common excitation schemes - darker bars represent a higher amplitude level.

3. Harmonic suppression via metaheuristics

3.1. A general overview on metaheuristic optimization

Metaheuristics may be defined as the nature-inspired algorithms which mimics a natural phenomenon in order to achieve a certain goal. To be specific, the main purpose of this kind of approach is to find a close-optimal parameter set over a specifically constructed mathematical space so-called a search space with some limiting constraints for the problem at hand. Most of the metaheuristic algorithms operate on a point-based search mechanism in order to keep the dimension space unchanged. In other words, these kinds of algorithms operate on a search space sampling strategy which forms a parallel optimal point search mechanism over the search space. These search space samples form a potential solution set called, usually, the population. If there is no pre information about search space behavior, usually, these optimizers start their search from a uniformly distributed search space sample set so-called the initial population.

Here, the search space mentioned above is a user-defined mathematical function space, usually called as cost function, fitness function, and so on, which varies depending on the problem construction. Definition of the search space is one of the key issues for these kinds of algorithms since it directly defines the problem. A poorly defined search space may lead inappropriate results even the algorithm works properly. A search space may be constructed in several ways but constructing it on an error function such as mean square or absolute error is a popular approach since the optimal point in an error function is predictable or known. However, using an error-based function is not a necessity, and some other well-defined functions may be used as the function that is to be optimized. If an error-based cost function is formed for a problem at hand, typically the problem becomes a minimization process, and the algorithm tries to find the possible solutions close to a desired one with a minimal error. A typical error-based cost function may be defined as

$$\begin{array}{cc}\underset{\overline{x}}{\text{minimize}}& f\left(\overline{x}\right)\\ \text{subject to}:& \overline{y}\end{array}$$ E15

where $\overline{x}$ stands for the parameter vector to be optimized, $\overline{y}$ is the constraint set and $f\left(\overline{x}\right)$ is the cost function may be defined by $f\left(\overline{x}\right)={\left[\sum _i{\left(g_i^d\left(\overline{x}\right)-g_i^a\left(\overline{x}\right)\right)}^2\right]}^{1/2}$ . Here, in the cost function, g, d and a represent the specification that is to be optimized, the desired level and the actual level, respectively.

After constructing and sampling the search space, that is, defining the cost function and constructing the initial population, all point-based metaheuristics move these initial samples to different points by their predefined operators. Actually, at this point, it should be noted that mainly the algorithms differ from each other by their relocation processes. Each algorithm uses a different method to relocate the initial samples, and usually this relocation process gives the algorithms name. By the relocation of the initial samples, a newly generated sample points are created which stands for new solutions. Comparing these new solutions with old ones and keeping the suitable ones and discarding the others cause an overall improvement in the solution set. In other words, after creating new points by relocation process, comparing them with already known possible solutions and keeping the better ones provide an overall convergence to the optimal point. This comparing and selecting process may be conducted in several ways such as using a roulette wheel or a binomial selection. In any way, the basic idea behind this process is comparing the solutions and keeping the better ones in order to improve the solution quality at hand. The relocation and selection process continues until a satisfied convergence rate is achieved or an iteration limit is reached. By this way a close optimal point may be found after the overall process in an iterative manner.

The basic idea behind the metaheuristics may be defined in main four steps as

1. Step 1: Constructing a search space

2. Step 2: Forming a randomly generated initial sample set

3. Step 3: Moving the samples in search space and create new solutions

4. Step 4: Comparing the new set with old one and keeping the suitable ones

As stated above, Step 3 and Step 4 are repeated in an iterative manner to achieve a satisfactory result. Here it should be noted that, because of their randomized nature, none of the metaheuristics guarantees to find an optimal solution. But it can be said that elitist methods (the methods that always keeping the current best solution in the solution set) with enough randomization of the solution set (this process is usually necessary to avoid local optimal points) guarantees to find a close optimal solution with enough convergence rate as the time tends to infinity.

3.2. Application to the harmonic suppression problem

As all metaheuristic approaches, definition of the problem is the key issue in sideband suppression problems. In TMA design which the sidebands are not to be used on specific purpose, the sideband radiation levels are preferred to be as low as possible in order to concentrate the total power in main radiation, lowering the total signal interference and unwanted signal radiations. However, suppressing sideband radiations is not the only problem while designing a switched antenna array. The most important problem needs to be solved is shaping main radiation beam with required side lobe level which may be a ‘must’ of design procedure since it affects communication quality directly. Then, other parameters which affect the performance such as harmonic-level suppression, harmonic beamforming and spatial filtering (i.e., adaptive beamforming to use the power effectively and null-steering to reduce interference) need to be handled to have the best fit design for application. Thus, if all the constraints are taken into consideration at the same time, different optimization algorithms become practical once today’s computers’ calculation speed, performance and technology are getting better day by day.

As mentioned earlier, designing a switching array via a metaheuristic algorithm starts with the definition of the problem at hand which will directly appears as the cost function that the algorithm uses. After the problem definition, the switching strategy should be decided which directly affects the parameter vector construction. These two steps should be taken into account carefully. First, a possible mistake in problem definition causes a defected mathematical search space, actually the search space will not be mathematically defected, but wrong statements will lead completely different search space that does not related to the actual problem. The algorithm will produce some results but completely unrelated to the problem actually at hand. Second, an improper selection of the excitation strategy may lead the problem unsolvable. As a simple example to this situation, for a problem involving only the suppression of specific harmonic levels, the selection of a VAS or PS switching scheme makes the problem unsolvable, since the harmonic levels are bounded to each other in these schemes. Hence, the independent suppression of different harmonics is mathematically impossible for this problem using the mentioned time schemes. At this point, it is beneficial to decide which conventional array design parameters such as element orientation or excitation phase to use in addition to the ‘time’ parameter. Since many problems may be solved via only correct switching strategy, this selection is not crucial and the traditional parameters may be excluded from the problem. However, adding some extra degrees of freedom will relax the solution of the problem and in some specific cases such as main beam steering, inclusion of the excitation phase to the parameter vector is inevitable. From this point on, the rest of the procedure is the application of the optimizer to the problem. The key points in the TMA design may be summarized as

• Defining the problem correctly,

• Selecting the switching strategy correctly,

• Deciding whether the conventional parameters are necessary.

The first study dealing with harmonic suppression via an optimizer has been appeared by Yang et al. in 2002 [5]. In Yang’s study, the switch-on durations of a VAS schemed linear array has been tried to be optimized. Unlike conventional array design, in designing a switched array, the harmonic radiations should be taken into account and the question of ‘how?’ is an issue that needs to be considered. At early researches, the calculation of the sidebands was conducted as a sidelobe-like algorithmic procedure, and the result of the procedure is used in metaheuristic algorithms as in Yang’s study [6, 7]. In these researches, the pattern of the sidebands has been sampled in a certain precision, and the maximum of these sample set has been taken as the maximum sideband level (SBL) which has been used as the parameter to be minimized. The basic idea behind this approach is to lower the maximum level of the infinite harmonics below a certain level usually a certain communication threshold. However, her the main problem is: there exist infinite harmonics in number and the question is how many of them should be considered? For this issue, usually the first M harmonics are taken into account in practice and M depends on the problem at hand. However, there also exists a common opinion that conducting operation only on the first harmonic is enough to extract the maximum harmonic level since the Fourier coefficients should follow an envelope of a ‘sinc’ function. This idea is not entirely true since there exist situations that the first harmonic is not the highest especially in forced cases (e.g., the situation in [23]). SBL suppression is a generally applicable generic way to gather information about sidebands; however, calculating sidebands from array factor definition with a certain accuracy is a time- and system resource-consuming process even conducting operation in the first harmonic. If first M harmonics are considered, the time consumption is M fold for SBL extraction. Over and above, in azimuthal asymmetric cases such as planar and volumetric situations, the calculation time is enormously increasing exponentially in proportion to the azimuth resolution. Hence, there exists a solution time-accuracy trade-off in this technique, which always needed to be considered. An illustrative block diagram of this technique may be found in [31].

In the cost function of an optimizer, this SBL information may be directly added. As an example, for a linear array and a problem involving only sideband suppression, the cost function of the metaheuristic algorithm may be written as

$$f=H\left(\mathrm{\Lambda }\right){\text{max}}_{\overline{x}}\left[{\mu }_{\overline{m}}\left(\overline{x},\theta \right)\right]$$ E16

where $H\left(\bullet \right)$ represents the standard Heaviside step function. Here, $\mathrm{\Lambda }=\left(\underset{\overline{x}}{\text{max}}\left[{\mu }_{\overline{m}}\left(\overline{x},\theta \right)\right]-{\mathrm{\Xi }}_d\right)$ , $\mu \left(\bullet \right)$ is the array factor, ${\mathrm{\Xi }}_d$ is the desired sideband level and $\overline{m}$ represents the harmonic number vector where $\forall m\in {\mathbb{Z}}^{+}$ , $m\le M$ and $M$ represents the harmonic number taken into account. Additionally, it should be kept in mind that in switched array design, the parameter vector $\overline{x}$ should contain the switching parameters. For this kind cost function, the expected value of the cost of the solution is zero. In other words, the optimizer should stop when the cost value reaches zero and the parameter vector $\overline{x}$ that satisfies $f=0$ regarded as a solution to the problem.

Extracting harmonic information to use in an optimizer from harmonic pattern samples is a generic way, but it is not the only way to gather information about sidebands. Calculating the power in sidebands and forming a bound function may also be used in sideband suppression problems. These two methods will be discussed in separate sections.

4. Total harmonic power radiation calculation

4.1. Harmonic power calculations in literature

In 2008, Bregains et al. published a paper that contains a closed form equation which defines an asymptotic approximation of the total power in infinite harmonic radiations for a VAS schemed linear array [25]. The main idea behind this approach is the total power of the radiated field, which may be expressed as the integral of the absolute squared array pattern over elevation under far field and $\kappa \gg 1$ approximations. Since the infinite series of absolute squared complex Fourier coefficients is convergent which appears in the raw equation, the infinite summation reduces to a specific function involving switch-on durations. After taking the integration over elevation, the result becomes an elegant closed form equation which gives an asymptotic approximation of the total power in infinite harmonics. After this first paper, Poli et al. published a paper that extends this equation to a rectangular grid planar array [26], and Aksoy and Afacan rewrite the equation for arbitrarily distributed planar case [27]. The main idea behind these papers is the integration of the array factor in planar case, which may be expressed as ordinary Bessel functions. After this conversion, the rest of the way to writing the equation is similar to the original Bregains’ paper. One handicap of this equation is that it only holds for the VAS scheme. For a shifted case (i.e., PS scheme), the results of the equation have been considered as inappropriate, since the equation does not contain information about starting instants. In 2012, this situation is added and combined with original equation by Aksoy and Afacan [28]. The idea behind this attempt is writing the infinite summation of the absolute squared complex Fourier coefficients as the Riemann’s zeta function. From this moment on, the equation has been valid for VAS and PS switched linear and planar arrays. But the equation still did not include volumetric cases, and in 2015, Aksoy published a paper to close this gap [29]. The key point in volumetric calculations is solving a definite integral over elevation involving an ordinary Bessel function times a complex exponential function. Writing this multiplication using spherical Bessel and Legendre functions makes this integral easily solvable, and the result reduces to a zeroth order spherical Bessel function which is equal to a ‘sinc’ function with argument of Euclidean distances as in linear and planar cases. Hence, by the publication of [29], the equation has been taken its final and more general form for VAS and PS switched arrays. Here, it must be noted that the studies mentioned so far have been conducted using ideal pulses (i.e., rectangular pulses with no transition regions). In addition to the studies given above, a pulse model having transition region has been studied by Bekele et al. and may be a reference to practical pulse models [32]. The summary of the derivation of the general form of power equation (from this point on, it will be referred as SR equation) is given in the next subsection, and for more information, references given above may be followed.

4.2. Derivation of generalized total harmonic power calculation

In this section, the general form of the harmonic power equation will be briefly derived. Before beginning derivation, some remarks should be noted. First of all, the equation derived here is an asymptotic approximation which the far field and $\kappa \gg 1$ approximations are assumed. Since the wave front of a spherically propagating wave may be assumed as planar in Fraunhofer region, this assumption helps to write the array factor in a simple summation consisting of planar wave front oriented phase shifts. Hence, in the region where the maximum phase error is below a certain tolerable level, the wave may be assumed as planar wave. Second, as stated earlier, the Poynting vector of a switched array may be written as in ordinary cases under $\kappa \gg 1$ approximation which makes the power density calculations as in conventional arrays.

Under these approximations, assume that $N$ isotropic radiators are randomly oriented in a three-dimensional space. In this general case, the time average power at harmonics may be represented as

$${\mathcal{P}}_H=\frac{1}{2}{\int }_0^{2\pi }{\int }_0^{\pi }\sum _{\begin{array}{c}q=-\infty \\ q\ne 0\end{array}}^{\infty }{\left|{\mu }_q\left(\theta ,\varphi \right)\right|}^2\text{sin}\mathit{\text{θd}}\mathit{\text{θdϕ}}.$$ E17

Here, ${\mu }_q\left(\theta ,\varphi \right)$ is being the array factor of isotropic sources and can be written for volumetric arrays as

$${\mu }_q\left(\theta ,\varphi \right)=\sum _{n=0}^{N-1}{W}_n{G}_{\mathit{nq}}{e}^{\mathit{ik}{x}_n\text{sin}\theta \text{cos}\varphi }{e}^{\mathit{ik}{y}_n\text{sin}\theta \text{sin}\varphi }{e}^{\mathit{ik}{z}_n\text{cos}\theta }$$ E18

where $W_n$ and $G_{\mathit{nq}}$ represent the complex excitations and complex Fourier coefficients, respectively. $x_n,y_n$ and $z_n$ stand for the coordinates of the array element in three orthogonal axis in standard Cartesian coordinates, $k=2\pi /{\lambda }_o$ is being the wavenumber and $i=\sqrt{-1}$ . In Eq. (18), the absolute squared array factor may be written as ${\left|{\mu }_q\left(\theta ,\varphi \right)\right|}^2={\mu }_q\left(\theta ,\varphi \right){\mu }_q^{\ast }\left(\theta ,\varphi \right),$ and after simple manipulations, ${\mathcal{P}}_H$ turns out to be.

$$\begin{array}{c}{\mathcal{P}}_H=\frac{1}{2}\sum _n{\left|W_n\right|}^2{\int }_0^{2\pi }{\int }_0^{\pi }\text{sin}\mathit{\text{θdθdϕ}}\sum _{\begin{array}{c}q=-\infty \\ q\ne 0\end{array}}^{\infty }{\left|G_{\mathit{nq}}\right|}^2\hfill \\ +\frac{1}{2}\sum _{\begin{array}{c}m,n=0\\ m\ne n\end{array}}^{N-1}{W}_m{W}_n^{\ast }{\int }_0^{2\pi }{\int }_0^{\pi }{e}^{\mathit{ik}{x}_n\text{sin}\theta \text{cos}\varphi }{e}^{\mathit{ik}{y}_n\text{sin}\theta \text{sin}\varphi }{e}^{\mathit{ik}{z}_n\text{cos}\theta }\text{sin}\mathit{\text{θdθdϕ}}\sum _{\begin{array}{c}q=-\infty \\ q\ne 0\end{array}}^{\infty }{G}_{\mathit{mq}}{G}_{\mathit{nq}}^{\ast }.\end{array}$$ E19

Here, the Fourier coefficients depend on the switching strategy. If an ideal PS time scheme consisting of rectangular pulses is considered, the switching function may be modeled as

$$g_n\left(t\right)=\left\{\begin{array}{c}1,0\le t_n^1<t<t_n^2\le T_p\\ 0,\mathit{\text{otherwise}}\hfill \end{array}\right.$$ E20

where $t_n^1$ and $t_n^2$ represent the starting and finishing instants of the corresponding pulse. For this time scheme, the Fourier coefficients may be written as

$$G_{\mathit{nq}}=\frac{i}{2\mathit{\pi q}}\left[e^{-\mathit{iq}2\pi t_n^2}-e^{-\mathit{iq}2\pi t_n^1}\right]$$ E21

where $q$ stands for the harmonic number and the infinite summation of $G_{\mathit{mq}}{G}_{\mathit{nq}}^{\ast }$ product may be examined for $m=n$ and $m\ne n$ cases.

For $m=n$ , this product may be written as

$$G_{\mathit{nq}}{G}_{\mathit{nq}}^{\ast }=\frac{1}{2{\pi }^2{q}^2}\left[1-\text{cos}\left(2\mathit{\pi q}\left(t_n^2-t_n^1\right)\right)\right].$$ E22

Hence, the infinite summation of this product turns out to be

$$\sum _{\begin{array}{c}q=-\infty \\ q\ne 0\end{array}}^{\infty }{G}_{\mathit{nq}}{G}_{\mathit{nq}}^{\ast }=\sum _{q=1}^{\infty }\frac{1}{{\pi }^2{q}^2}-\sum _{q=1}^{\infty }\frac{1}{{\pi }^2{q}^2}\text{cos}\left(2\mathit{\pi q}{\tau }_n\right)$$ E23

where ${\tau }_n=\left(t_n^2-t_n^1\right)/T_p$ represents the normalized on-time duration of nth element. Since, this series is convergent, using the Riemann’s zeta function where ${\sum }_{x=1}^{\infty }1/x^n=\zeta \left(n\right)$ , the infinite series of coefficient product reduces to

$$\sum _{\begin{array}{c}q=-\infty \\ q\ne 0\end{array}}^{\infty }{G}_{\mathit{nq}}{G}_{\mathit{nq}}^{\ast }=\frac{1}{{\pi }^2}\frac{{\pi }^2}{6}-\frac{1}{{\pi }^2}\left\{\frac{{\pi }^2}{6}-\frac{\pi }{2}2\pi {\tau }_n+\frac{{\left(2\pi {\tau }_n\right)}^2}{4}\right\}={\tau }_n\left(1-{\tau }_n\right).$$ E24

For $m\ne n$ by using the fact that cosine is an even function and sine is odd, this infinite series turns out to be

$$\begin{array}{c}\sum _{\begin{array}{c}q=-\infty \\ q\ne 0\end{array}}^{\infty }{G}_{\mathit{nq}}{G}_{\mathit{mq}}^{\ast }=\frac{1}{2{\left(\mathit{\pi q}\right)}^2}\sum _{q=1}^{\infty }\left\{\text{cos}\left[2\mathit{\pi q}\left(t_n^1-t_m^1\right)\right]-\text{cos}\left[2\mathit{\pi q}\left(t_n^1-t_m^2\right)\right]\right\}\\ +\frac{1}{2{\left(\mathit{\pi q}\right)}^2}\sum _{q=1}^{\infty }\left\{\text{cos}\left[2\mathit{\pi q}\left(t_n^2-t_m^2\right)\right]-\text{cos}\left[2\mathit{\pi q}\left(t_n^2-t_m^1\right)\right]\right\}.\end{array}$$ E25

By using relation for summable series, see [25, Appendix], Eq. (25) reduces to

$$\sum _{\begin{array}{c}q=-\infty \\ q\ne 0\end{array}}^{\infty }{G}_{\mathit{nq}}{G}_{\mathit{mq}}^{\ast }=\frac{1}{2}\left(-\left|t_n^1-t_m^1\right|+\left|t_n^1-t_m^2\right|-\left|t_n^2-t_m^2\right|+\left|t_n^2-t_m^1\right|-2{\tau }_n{\tau }_m\right).$$ E26

For above equation, there exist 53 different situations in terms of equality and inequality conditions of switching instants. However, considering the fact that $t_n^2>t_n^1$ and $t_m^2>t_m^1$ should be satisfied always, these possibilities reduce to six. If these six equality and inequality possibility of starting and finishing instants of the pulses are considered, these six conditions may be written in one form given by

$$\sum _{\begin{array}{c}q=-\infty \\ q\ne 0\end{array}}^{\infty }{G}_{\mathit{nq}}{G}_{\mathit{mq}}^{\ast }=\overline{{\tau }_{\mathit{mn}}}-{\tau }_n{\tau }_m.$$ E27

In this expression, the bar represents the duration of the intersection between n^th and m^th element’s pulses. For more detailed calculations, [28] may be followed. By substituting this result into Eq. (19), it becomes

$${\mathcal{P}}_H=2\pi \sum _n{\left|W_n\right|}^2\left[{\tau }_n\left(1-{\tau }_n\right)\right]+\frac{1}{2}\sum _{\begin{array}{c}m,n=0\\ m\ne n\end{array}}^{N-1}{W}_m{W}_n^{\ast }\left(\overline{{\tau }_{\mathit{mn}}}-{\tau }_m{\tau }_n\right)\mathrm{\Omega }\left(\theta ,\varphi \right)$$ E28

where $\mathrm{\Omega }\left(\theta ,\varphi \right)$ is being

$$\mathrm{\Omega }\left(\theta ,\varphi \right)={\int }_0^{2\pi }{\int }_0^{\pi }{e}^{\mathit{ik}{x}_n\text{sin}\theta \text{cos}\varphi }{e}^{\mathit{ik}{y}_n\text{sin}\theta \text{sin}\varphi }{e}^{\mathit{ik}{z}_n\text{cos}\theta }\text{sin}\mathit{\text{θdθdϕ}}$$ E29

and it should be solved. Solving this integral is not an easy task, and some manipulations should be conducted. It can be started from to write some exponential terms in terms of Bessel functions. To do that, the following fact may be used

$$\text{acos}x+b\text{sin}x\sqrt{a^2+b^2}\text{cos}\left(x-{\text{tan}}^{-1}\left[\frac{b}{a}\right]\right).$$ E30

By using this conversion and integral definition of ordinary Bessel functions with $a=k\left(x_m-x_n\right)\text{sin}\theta$ and $b=k\left(y_m-y_n\right)\text{sin}\theta$ , the $\mathrm{\Omega }\left(\theta ,\varphi \right)$ can be written as

$$\mathrm{\Omega }\left(\theta ,\varphi \right)=2\pi {\int }_0^{\pi }{J}_0\left(\sqrt{a^2+b^2}\right)e^{\mathit{ik}\left(z_m-z_n\right)\text{cos}\theta }\text{sin}\mathit{\theta d\theta }.$$ E31

Here, $J_0\left(\bullet \right)$ represents the ordinary Bessel function of order zero. By letting $c_1=k\sqrt{{\left(x_m-x_n\right)}^2+{\left(y_m-y_n\right)}^2}$ and $c_2=k\sqrt{{\left(z_m-z_n\right)}^2}$ with $\mathrm{\Omega }\left(\theta ,\varphi \right)=2\pi {\mathrm{\Omega }}^{\prime }\left(\theta ,\varphi \right)$ , Eq. (31) becomes

$${\mathrm{\Omega }}^{\prime }\left(\theta ,\varphi \right)={\int }_0^{\pi }{J}_0\left(c_1\text{sin}\theta \right)e^{i{c}_2\text{cos}\theta }\text{sin}\mathit{\theta d\theta }.$$ E32

Here, $c_1$ and $c_2$ are orthogonal to each other and can be represented as $c_1=\mathit{\text{Rcos}}\psi$ and $c_2=\mathit{\text{Rsin}}\psi$ where $R=k{R}^{\prime }$ and $\psi ={\text{cos}}^{-1}{\widehat{a}}_z·\left({\mathbf{r}}_m-{\mathbf{r}}_n\right)$ . Here, $R^{\prime }$ represents the Euclidean distance between corresponding elements (i.e., n^th and m^th elements under process) and ${\mathbf{r}}_m\left(x_m,y_m,z_m\right)$ and ${\mathbf{r}}_n\left(x_n,y_n,z_n\right)$ stand for the position vectors m^th and n^th elements. By using this conversion, Eq. (32) can be written as

$${\mathrm{\Omega }}^{\prime }\left(\theta ,\varphi \right)={\int }_0^{\pi }{J}_0\left(R\text{sin}\theta \text{sin}\psi \right)e^{\mathit{iR}\text{cos}\theta \text{cos}\psi }\text{sin}\mathit{\theta d\theta }.$$ E33

Once the ${\mathrm{\Omega }}^{\prime }\left(\theta ,\varphi \right)$ is written as in Eq. (33), the following relation may be used instead

$$J_0\left(R\text{sin}\theta \text{sin}\psi \right)e^{\mathit{iR}\text{cos}\theta \text{cos}\psi }=2\sum _{n=0}^{\infty }{i}^n\left(n+\frac{1}{2}\right)j_n\left(R\right)P_n\left(\text{cos}\theta \right)P_n\left(\text{cos}\psi \right).$$ E34

Here, $j_n\left(·\right)$ and $P_n\left(·\right)$ represent the n^th order spherical Bessel and Legendre functions, respectively. Multiplying both sides of Eq. (34) by $P_l\left(\text{cos}\theta \right)\text{sin}\theta$ and integrate over $\left[0,\pi \right]$ results in the following relation

$$\begin{array}{c}{\int }_0^{\pi }{J}_0\left(R\text{sin}\theta \text{sin}\psi \right)e^{\mathit{iR}\text{cos}\theta \text{cos}\psi }{P}_l\left(\text{cos}\theta \right)\text{sin}\mathit{\theta d\theta }\hfill \\ =\sum _{n=0}^{\infty }{i}^n\left(2n+1\right)j_n\left(R\right)P_n\left(\text{cos}\psi \right){\int }_0^{\pi }{P}_n\left(\text{cos}\theta \right)P_l\left(\text{cos}\theta \right)\text{sin}\mathit{\theta d\theta }\end{array}$$ E35

where

$${\int }_0^{\pi }{P}_n\left(\text{cos}\theta \right)P_l\left(\text{cos}\theta \right)\text{sin}\mathit{\theta d\theta }=\frac{2}{2n+1}{\delta }_{n,l}.$$ E36

In this expression, ${\delta }_{n,l}$ represents the Kronecker delta function defined by

$${\delta }_{n,l}=\left\{\begin{array}{c}1,n=l\\ 0,n\ne l\end{array}\right..$$ E37

Hence, for $l=0$ , the equation Eq. (35) becomes

$${\int }_0^{\pi }{J}_0\left(R\text{sin}\theta \text{sin}\psi \right)e^{\mathit{iR}\text{cos}\theta \text{cos}\psi }\text{sin}\mathit{\theta d\theta }=2j_0\left(R\right)$$ E38

and $\mathrm{\Omega }\left(\theta ,\varphi \right)$ can be written as

$$\mathrm{\Omega }\left(\theta ,\varphi \right)=2\pi {\mathrm{\Omega }}^{\prime }\left(\theta ,\varphi \right)=4\pi j_0\left(R\right)=4\pi \frac{\text{sin}\left(R\right)}{R}.$$ E39

If this result is substituted to Eq. (19) with Eq. (24) and Eq. (27), total power approximation may be finally written as

$${\mathcal{P}}_H=2\pi \sum _n{\left|W_n\right|}^2\left[{\tau }_n\left(1-{\tau }_n\right)\right]+2\pi \sum _{\begin{array}{c}m,n=0\\ m\ne n\end{array}}^{N-1}{W}_m{W}_n^{\ast }\left(\overline{{\tau }_{\mathit{mn}}}-{\tau }_m{\tau }_n\right)\frac{\text{sin}\left(R\right)}{R}.$$ E40

Once more, in Eq. (40), $R=\sqrt{{\left(x_m-x_n\right)}^2+{\left(y_m-y_n\right)}^2+{\left(z_m-z_n\right)}^2}$ represents the Euclidean distance between corresponding elements, $I$ represents the complex excitations and $\tau$ represents the normalized switch-on durations for an $N$ element volumetric array. For a more detailed proof, [29] may be followed, and for linear and planar case, see [25] and [27], respectively.

4.3. Using total power in suppression problems

The first usage of total power in harmonic suppression problems was conducted by Poli et al. in 2010 [10]. The main idea behind this technique is based on the idea of the total power reduction in harmonics will concentrate the power to main radiation, hence the harmonics become suppressed. Since the result of the equation does not meet the actual power in harmonics, the ratio of the harmonic power to the actual power is more meaningful in practice. Hence, usually the direct result of the Eq. (40) is not preferred to be used in an optimizer. Instead of using the direct result, the ratio of the harmonic power to the actual power is more meaningful which may be written by

$$\acute{{\mathcal{P}}_H}=\frac{{\mathcal{P}}_H}{{\mathcal{P}}_H+{\mathcal{P}}_0}$$ E41

where ${\mathcal{P}}_0$ represents the power in main radiation given by

$${\mathcal{P}}_0=\frac{1}{2}{\int }_0^{2\pi }{\int }_0^{\pi }{\left|{\left.{\mu }_q\left(\theta ,\varphi \right)\right|}_{q=0}\right|}^2\text{sin}\mathit{\text{θd}}\mathit{\text{θdϕ}}$$ E42

Without loss of generality for a problem concerning only harmonic reduction, the usage of Eq. (41) in a cost function may be written as

$$f=H\left(\mathrm{\varpi }\right)\acute{{\mathcal{P}}_H}\left(\overline{I},\overline{\tau },\overline{R}\right)$$ E43

where $\mathrm{\varpi }=\left(\acute{{\mathcal{P}}_H}\left(\overline{I},\overline{\tau },\overline{R}\right)-{\mathrm{\Theta }}_d\right)$ while ${\mathrm{\Theta }}_d$ is being the desired harmonic to main power ratio.

5. Harmonic level bound (HLB)

5.1. HLB concept

As shown in [30] that neither harmonic level reduction nor power reduction ensures the total suppression in terms of both power and communication level. Hence, a combined way may be a more suitable approach if the both level and power reduction is necessary. At this point, the question of ‘is there any simple method exist for the level calculation in order to avoid time consuming processes?’ may be emerged in mind, and it can be said that such studies exist in the literature as an answer to this question.

As mentioned earlier, the SBL method extracts the harmonic level information in an algorithmic way similar to sidelobe calculations. Since there exists no analytical solution to find a maximum point in an unknown sidelobe region, the sidelobe calculations are being conducted such an operation involving sampling the pattern in sidelobe region and finding its maximum. In contrast, the sideband calculations are being conducted in whole visible region, which makes a difference in both SLL and SBL calculations. In other words, since the SBL calculations are being operated over a complete elevation and azimuth space, in some cases, the maximum of a pattern may be extracted analytically. Since the maximum points of all individual harmonic patterns can be calculated, they form a maximum sideband level set, and the maximum point of this set bounds the whole harmonic maxima. Hence, this maximum appears as a bound function covering all individual maximum harmonic levels which are infinite in number. By this way, lowering the total bound ensures that all harmonic levels are below this level.

First attempt to calculate a bound is conducted by Aksoy for linear arrays with a VAS time scheme [23]. In this attempt, the Poynting vector of a TMA has been considered, and the pattern equations are written as the normalized time average power densities. After this first attempt, Aksoy and Afacan published a proof including the planar case using the same idea [24]. By this way, an overall harmonic bound for a VAS switched time modulated linear and planar arrays has been written as

$$\mathrm{\Psi }={\left(\frac{{\sum }_{n=0}^{N-1}{I}_n\text{sin}\left(\pi {\tau }_n\right)}{\pi {\sum }_{n=0}^{N-1}{\xi }_n}\right)}^2$$ E44

where $\mathrm{\Psi }$ stands for the normalized upper sideband bound, $I_n$ and ${\tau }_n$ represent the excitation amplitude and the normalized on-time durations of the n^th element, respectively. Also, ${\xi }_n$ is the dynamic excitation defined by ${\xi }_n=I_n{\tau }_n$ and $N$ represents the total element number in the array. Here it must be noted that, to clear possible misunderstandings, the excitation vector represents only the amplitude of the complex excitations and defined in region $\forall I\in {\mathbb{R}}^{+},0<I\le 1$ . Similarly, the normalized switch-on durations are also defined in region $\forall \tau \in {\mathbb{R}}^{+}$ , where $0<\tau \le 1$ . Since this equation is a closed form expression, the total computation time for harmonic information is tremendously short as compared to the SBL method. Since it is in a shorter form, the computation time is slightly better compared to SR equality. A detailed comparison may be found in [31].

Like all switched array calculations, this technique also uses an asymptotic approach and valid only under far-field and $\kappa \gg 1$ approximations and has its own advantages and disadvantages. The major advantage of this technique is that it ensures the harmonic level suppression among infinite sidebands in a very short time. However, currently, it is only applicable to linear and planar arrays with VAS time scheme, which may be considered as a handicap of this technique.

5.2. HLB applications

If the harmonic levels are wanted to be suppressed via the bound function in a metaheuristic approach, the result produced by the bound function may be used directly in the cost function of the optimizer. The aim of the optimization is to find a switch-on duration sequence to satisfy a desired level. Hence, the parameter vector should contain at least switch-on durations $\overline{\tau }$ . For a uniform array, the excitation $\overline{I}$ may be taken as a vector whose elements are set to unity. Without loss of generality for a linear case with only sideband level, reduction is necessary for such a cost function that may be written as

$$f=H\left(\mathrm{\lambda }\right)\mathrm{\Psi }\left(\overline{I},\overline{\tau }\right)$$ E45

where $H\left(\bullet \right)$ stands for the standard Heaviside step function and $\mathrm{\lambda }=\left(\mathrm{\Psi }\left(\overline{I},\overline{\tau }\right)-{\mathrm{\Xi }}_d\right)$ while ${\mathrm{\Xi }}_d$ is being the predefined desired harmonic level. Here, the $\left\{\overline{\tau },\overline{I}\right\}$ set stands for a solution for $f=0$ .

Moreover, the equality has been derived for the VAS scheme, but it can be used for any switching strategy involving one pulse in one period. On the other hand, the efficiency of the bound function in shifted cases does not appear as good as in the VAS scheme, since the equality does not contain information about starting instants. In practice, the actual harmonic levels are usually less than, or at least equal to, the results produced by the bound function. In VAS scheme, the difference between actual and the bound is usually getting smaller as the iteration counts, but the same statement does not hold true for the PS scheme. Hence, using the bound function in shifted cases causes a random walk on starting instants, and this situation may lead to some impractical results.

6. A non-iterative algorithmic approach to suppress the harmonics using variable pulse amplitude

6.1. Variable pulse amplitude (VPA)

Instead of using “on–off” switching scheme, a novel approach is proposed by Aksoy in 2014 named as “Variable Pulse Amplitude (VPA),” which is based on switching each element between amplifiers [12]. The main purpose of the newly proposed method is preventing the array silencing, which means all the elements are switched to “off” position and no communication at that instant (e.g., see Figure 1). The binary “On – off” switching with modifications such as shifting and splitting allows the array to be silenced even for a very short time instant and the received or transmitted data may be missed at defined instants that all the elements are switched off in every period of switching. It may be seemed as minor possibility, but it results in additional optimization constraint to be checked whether there is an instant that all the elements are turned-off or not.

6.2. Harmonic suppression method

It is stated in previous sections that metaheuristic optimization algorithms do not guarantee to find the proper solution; however, they generally produce acceptable results. In addition, optimization time may be longer than usual according to host computer performance, optimization constraints and iteration count. On the contrary, a noniterative approach offers a fast and reliable solution with respect to metaheuristic methods, and it is preferable if exists.

In this part of section, a novel approach in order to suppress the sideband levels using variable pulse amplitude excitation scheme will be presented with an explanatory example.

6.2.1. Mathematical expressions

The array factor of time-modulated array is given in Eq. (14), and it is easy to understand that complex Fourier coefficients (CFC), $G_{\mathit{nq}}$ , affect the individual radiation patterns since excitation amplitude, excitation phase and placement of the elements are the same for main radiation and sideband radiations. In order to calculate $G_{\mathit{nq}}$ values, VPA excitation strategy may be expressed as

$$g_n\left(t\right)=\left\{\begin{array}{c}{K}_n^1,0<t<t_n\\ K_n^2,t_n<t<T_p\end{array}\right.$$ E46

where $K_n^1\text{and}{K}_n^2$ are the adjacent pulse amplitudes and $t_n$ represents amplitude reversal time. Hence, the $g_n\left(t\right)$ in Eq. (14) should be replaced with Eq. (46) whose CFCs are defined by

$$G_{\mathit{nq}}=\frac{1}{T_p}\underset{-\infty }{\overset{\infty }{\int }}{g}_n\left(t\right)e^{-\mathit{iq}{\omega }_{p}t}\mathit{dt}=\frac{1}{T_p}\underset{0}{\overset{t_n}{\int }}{K}_n^1{e}^{-\mathit{iq}{\omega }_{p}t}\mathit{dt}+\frac{1}{T_p}\underset{t_n}{\overset{T_p}{\int }}{K}_n^2{e}^{-\mathit{iq}{\omega }_{p}t}\mathit{dt}.$$ E47

In this equation, while $q=0$ represents main radiation, $\left|q\right|>0$ defines sideband radiations, and they can be calculated individually. The CFCs for main radiation may be obtained as follows:

$$G_{n0}=\frac{1}{T_p}\underset{0}{\overset{t_n}{\int }}{K}_n^1\mathit{dt}+\frac{1}{T_p}\underset{t_n}{\overset{T_p}{\int }}{K}_n^2\mathit{dt}={\tau }_n{∆}_n+K_n^2$$ E48

where ${\tau }_n=t_n/T_P$ and ${∆}_n=\left(K_n^1-K_n^2\right)$ . Here, ${\tau }_n$ and ${∆}_n$ represent normalized amplitude reversal time with respect to switching period and pulse amplitude difference of adjacent pulses, respectively. For the sideband radiations where $\left|q\right|>0$ , CFC values are defined by

$$G_{\mathit{nq}}=\frac{1}{T_p}\underset{0}{\overset{t_n}{\int }}{K}_n^1{e}^{-\mathit{iq}{\omega }_{p}t}\mathit{dt}+\frac{1}{T_p}\underset{t_n}{\overset{T_p}{\int }}{K}_n^2{e}^{-\mathit{iq}{\omega }_{p}t}\mathit{dt}=\frac{{∆}_n}{\mathit{iq}2\pi }\left[1-e^{-\mathit{iq}{\omega }_p{\tau }_n}\right].$$ E49

To sum up, after some basic manipulations, CFCs for all radiation levels are expressed as

$$G_{\mathit{nq}}=\left\{\begin{array}{c}{∆}_n{\tau }_n+K_n^2,q=0\\ \frac{{∆}_n}{\mathit{q\pi }}\text{sin}\left(\mathit{q\pi }{\tau }_n\right)e^{-\mathit{iq\pi }{\tau }_n},\left|q\right|>0\end{array}\right.$$ E50

Without loss of generality, since discrete windowing functions (e.g., Chebyshev or Taylor $\overline{n}$ distribution) are eligible to be used for controlling the sidelobe level of main radiation, CFC values may be selected equal to distribution coefficients for $q=0$ . In addition, for sideband radiations, $\left({∆}_n/\mathit{q\pi }\right)\text{sin}\left(\mathit{q\pi }{\tau }_n\right)$ expression appears as an amplitude term and may be used for suppressing the harmonic level.

In order to shape the main radiation pattern and to suppress the harmonic level, according to the abovementioned definitions, steps to be followed are given as:

1. Step 1: Use discrete window distribution coefficients, ${\alpha }_n$ , as main radiation CFCs (i.e., ${\alpha }_n={∆}_n{\tau }_n+K_n^2)$ .

2. Step 2: Set a constant ${\tau }_n$ value where ${\tau }_n\in \left(0,1\right)$ .

3. Step 3: Select amplitude term of fundamental radiation, which is equal to ${\alpha }_n/\gamma$ in order to suppress the radiation level where $\gamma$ is the suppression ratio and $\gamma \in {\mathbb{R}}^{+}$ . (i.e., ${\alpha }_n/\gamma =\left({\mathrm{\Delta }}_n/\pi \right)\text{sin}\left(\pi {\tau }_n\right)$ ).

4. Step 4: Calculate ${∆}_n$ values using step 3 (i.e., ${∆}_n=\frac{\pi {\alpha }_n}{\gamma \text{sin}\left(\pi {\tau }_n\right)}$ ).

5. Step 5: Calculate $K_n^2={\alpha }_n-{∆}_n{\tau }_n$ .

6. Step 6: Calculate $K_n^1={∆}_n+K_n^2$ .

7. Step 7: Normalize ${∆}_n,K_n^1\text{and}{K}_n^2$ values individually with respect to $\rho =\text{max}\left({∆}_n,K_n^1,K_n^2\right)$ to obtain final design parameters.

Using this approach, main radiation sidelobe level is directly adjusted according to windowing function parameters, and fundamental harmonic radiation is suppressed with the ratio of $\gamma$ . It is obvious that once the first harmonic radiation is suppressed, the other sideband radiations are self-suppressed for single pulsed schemes. In the next subsection, the method will be exemplified with an explanatory scenario.

6.2.2. A simple example

Let us assume a 10-element -30 dB Chebyshev array whose elements are located along z-axis with interelement spacing of $d=0.5\lambda$ . For simplicity, excitation amplitudes of each element are selected equal to unity (i.e., $I_n=1$ ), and no progressive phase shift is applied to the elements (i.e., ${\beta }_n=0$ ). Normalized pulse amplitude reversal time is selected as ${\tau }_n=0.45$ for each element, but it does not have any effect on maximum harmonic level since normalized pulse difference values (i.e., ${∆}_n$ ) are adjusted according to ${\tau }_n$ values.

Harmonic suppression ratio, $\gamma$ , should be selected using $20{\text{log}}_{10}\left(1/\gamma \right)=\delta$ where $\delta$ represents desired the maximum sideband radiation level in $\mathit{dB}$ . It should be noted again that with the adjustment of $\gamma$ , the maximum level of fundamental harmonic radiation is set to related dB value, thus, further sideband radiation levels are self-suppressed due to the nature of harmonic radiations.

In this simple example, maximum harmonic levels are selected as $-30\mathit{dB}$ and $-40\mathit{dB}$ so that harmonic suppression ratios should be $\gamma =31.6$ and $\gamma =100$ , respectively. Once $-30\mathit{dB}$ maximum level of sideband radiation is considered, main and first harmonic radiation patterns are given together in Figure 2a. Moreover, maximum levels of first 30 harmonic radiations are shown in Figure 2b. According to the calculation steps given in previous subsection, normalized pulse difference and adjacent pulse amplitude values that are used to obtain the results are given in Table 1.

Figure 2.

Results for $\gamma =31.6$ (harmonic level $=-30\mathit{dB}$ ) (a) Radiation patterns of $q=0$ and $q=1$ (b) Maximum levels of first 30 harmonic radiations.

n		1	2	3	4	5	6	7	8	9	10
$\gamma =31.6$	$K_n^1$	0.258	0.430	0.669	0.878	1.000	1.000	0.878	0.669	0.430	0.258
	$K_n^2$	0.233	0.389	0.605	0.794	0.905	0.905	0.794	0.605	0.389	0.233
	${∆}_n$	0.025	0.041	0.064	0.084	0.095	0.095	0.084	0.064	0.041	0.025
$\gamma =100$	$K_n^1$	0.258	0.430	0.669	0.878	1.000	1.000	0.878	0.669	0.430	0.258
	$K_n^2$	0.249	0.417	0.648	0.851	0.969	0.969	0.851	0.648	0.417	0.249
	${∆}_n$	0.008	0.013	0.021	0.027	0.031	0.031	0.027	0.021	0.013	0.008

Table 1.

Normalized pulse difference and adjacent pulse amplitude values for $\gamma$ =31.6 and $\gamma$ =100.

If $-40\mathit{dB}$ of maximum harmonic level is taken into consideration (i.e., $\gamma =100$ ), Figure 3a and Figure 3b display the radiation patterns for $q=0$ together with $q=1$ and maximum levels of first 30 harmonic radiations, respectively. ${∆}_n,K_n^1$ and $K_n^2$ values are also shown in Table 1.

Figure 3.

Results for $\gamma =100$ (harmonic level $=-40\mathit{dB}$ ) (a) Radiation patterns of $q=0$ and $q=1$ (b) Maximum levels of first 30 harmonic radiations.

7. Conclusions

Nowadays, communication technology moving rapidly toward 5G, and frequency spectrum is one of the most important issues in terms of operational and capital expenses for industry. In order not to use redundant frequency bands and to make efficient use of channel capacity, it is preferred to suppress the unwanted harmonic radiations. In this study, a general overview on harmonic suppression in 4D arrays using optimization methods is given. A brief mathematical background for switched arrays and the optimization basics are explained. Furthermore, common excitation strategies and the techniques used in harmonic calculations are summarized. More importantly, a noniterative algorithmic suppression strategy is introduced and exemplified via a simple example of the harmonic suppression of a 10-element linear array with sidelobe control. The results of the introduced noniterative strategy seem quite satisfactory in terms of calculation complexity as compared to suppression via an optimizer.