Frequency Interpolation of LOFAR Embedded Element Patterns Using Spherical Wave Expansion

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


Introduction
For the calibration of phased arrays, a priori knowledge of the embedded element patterns (EEPs) is needed. ese patterns can be modeled by an analytic function and are often assumed to be the same for each antenna element. To improve the calibration, a full-wave electromagnetic simulation can be used to provide more accurate EEPs [1]. However, if the array is large, the amount of data also becomes large. A way to reduce the amount of data is the use of spherical wave expansion (SWE) to describe the EEPs. e use of SWE related to radiation patterns is mainly seen in near-field to far-field transformations [2]. In [3], the use of SWE to describe the embedded element patterns of phased array antennas is described. One of the advantages of the use of SWE is that the correlation between two embedded element patterns can be computed efficiently [4].
In this paper, we will extend the use of spherical wave expansion with an interpolation technique in the frequency domain.
is technique uses the Fast Fourier Transform (FFT) to interpolate the spherical wave coefficients in the frequency domain.
ese interpolated coefficients can be used to calculate the embedded element patterns at nonsimulated frequencies. e interpolation method can be used for phased arrays and Phased Array Feeds (PAFs). Application of PAF's can be found in radio astronomy [5][6][7] and future 5G communication systems [8].
e LOFAR radio telescope [9] will be used as a test case to verify our proposed concept. e structure of this paper is as follows: Section 2 gives an overview of spherical wave expansion. Special attention is paid to the calculation of the coefficients which is done by using the Moore-Penrose pseudoinverse. e interpolation method is also discussed in this section. A description of the LOFAR radio telescope is given in Section 3. is section gives a general overview of the LOFAR system and a detailed description of the low-band antenna. e results are described in Section 4. Conclusions are given in Section 5. E → r, u r � e → u r e − jβr r , (1) where β is the free-space wave number, r the distance between the antenna and the observation point, e → (u r ) the normalized electric component of the far field, and u r a unit vector pointing towards the observation point defined as u r � u x sin θ cos φ + u y sin θ sin φ + u z cos θ, (2) where u x , u y , and u z are the unit vectors in the x-, y-, and z-direction, respectively. θ is the angle between u r and the z-axis measured from the z-axis and φ is the angle between the x-axis and the projection of u r on the xy-plane measured from the x-axis (see Figure 1). e normalized electric far field can be expanded in terms of spherical wave functions: where f → smn (u r ) is the spherical function (s, m, n) and Q smn the corresponding coefficient. It is convenient to use a socalled "compressed index" p defined as [2, p. 15] p � 2(n(n + 1) + m − 1) + s, (4) to replace the indices (s, m, n). e numbering scheme up to n � 2 is shown in Table 1.
Using the compressed index, equation (3) can be written as Each spherical wave function f → p (u r ) is normalized such that 2π φ�0 π θ�0 f → p u r (θ, φ) 2 sin θ dθ dφ � Z F , (6) where Z F is the wave impedance in free space (≈377Ω). e spherical wave functions f → p (u r (θ, φ)) and f → q (u r (θ, φ)) are orthogonal with respect to each other: Note that we used the notation u r (θ, φ) to denote that u r depends on θ and φ according to equation (2). e complex conjugate is denoted by * .

Calculation of the Coefficients Q p .
In practice, the summation over n in (3) is limited to a finite number N max . By setting n � m � N max and s � 2 in (4), one can show that the maximum number of modes N modes should be at least N modes should be chosen such that the contribution of the higher-order modes (>N modes ) can be neglected. As a rule of thumb, N max will be chosen as [2, p. 17] which corresponds with N modes � 2 βr 0 + 10 βr 0 + 12 � 2β 2 r 2 0 + 44βr 0 + 240, where r 0 is the radius of the smallest sphere enclosing the antenna. If the radiation pattern e → (u r ) of an antenna is known in N dir directions and we want to describe it in N modes spherical modes, it can be written as a vector e with length 2N dir : where e θ (u r,n ) and e φ (u r,n ) are, respectively, the θ-and φ-component of the normalized electric far field in the direction u r,n :   International Journal of Antennas and Propagation e → u r,n � e θ u r,n u θ + e φ u r,n u φ , where u θ and u φ are the unit vector in θ-and φ-direction, respectively. e far-field spherical vector functions for the first N modes modes specified in N dir directions can be written in a 2N dir × N modes matrix F: e scalars f θ,p (u r ) and f φ,p (u r ) are the θ-and φ-component of the spherical wave function f → p (u r ), respectively: where we used the compressed index p. e spherical mode coefficients Q p can be grouped in a vector q with length N modes : Equation (5) can now be written as e solution of equation (16) is given by [10] where the operator † represents the Moore-Penrose pseudoinverse, U Nmodes is an N modes × N modes identity matrix, and w is a column vector with length N modes in which elements can have arbitrary values. If e � FF † e, at least one solution exists. For the cases when no solutions exist, equation (17) gives the "closest" solution(s). Mathematically, these are the value or values of q for which the Euclidean norm ‖e − Fq‖ 2 has its absolute minimum. If FF † � U 2Ndi r (with U 2Ndi r a 2N di r × 2N di r identity matrix), equation (16) has at least one solution for arbitrary e. is requirement is only fulfilled if the number of rows of F is equal to or less than its number of columns and F has a full (row) rank. In practical cases, the number of rows of F is larger than its number of columns, so equation (16) has no solution. In practice, the coefficients Q p ≈ 0 for p > N modes , so the norm ‖e − Fq‖ 2 will be very close to zero, and equation (17) can practically be considered as a solution of equation (16). If F † F � U Nmodes , there exists only one solution or value for which the norm ‖e − Fq‖ 2 has its absolute minimum. is is only the case if the number of rows of F is equal to or larger than its number of columns and F has a full (column) rank. If F † F ≠ U Nmodes , there exist infinite solutions or values for which the norm ‖e − Fq‖ 2 has its absolute minimum. ese solutions or values are parametric on the vector w in equation (17). Figure 2 shows the rank of matrix F as a function of the number of columns (coefficients) for a different number of rows (angles). e number of columns represents the number of spherical wave coefficients. e number of rows is directly related to the number of directions for which the far field is given. is relation is shown in Table 2. e number of rows is two times the product of the number of θ-and φ-angles for which the far field is given (note that the far field has a θ-and φ-component for each direction). e angles θ � 0°and θ � 180°are not included because the spherical wave functions are not defined for these angles.
For an increasing number of columns (number of coefficients) for a given number of rows, the matrix F has full column rank until a certain number of columns. e highest number of columns where F has full column rank are denoted by * in Figure 2. e column rank still increases if the number of columns is further increased. From a certain point where the number of columns is higher than the number of rows, F has a full row rank. ese points are denoted by°in Figure 2. For the cases when F has 5040 and 25 776 rows, this number could not be determined because the size of F becomes too large. Table 3 shows for which number of columns the matrix F has full row and column rank. Table 3 shows, for example, that if matrix F has 288 rows (θ-and φ-angles in steps of 20°) it has full column rank if the number of columns is less than or equal to 160. is means that, for this case, equation (16) has one value for which the norm ‖e − Fq‖ 2 has its absolute minimum if the number of coefficients Q is chosen equal to or less than 160. If the number of coefficients is 510 or more, there exists more than one vector q for which the norm has its absolute minimum. Note that there is no number of coefficients for which there exists one exact solution. Figure 3 shows the condition number of the matrix F as a function of the number of columns for different numbers of rows.
e point with the maximum number of columns where F has full column rank and the point with the minimum number of columns where F has full row rank are marked in the same way as in Figure 2.
One can see that the condition number remains below 10 until the point where F reaches full column rank. After this point, the condition number increases to values greater than 10 15 . A large condition number means that calculations involving F become numerically unstable; that is, a small variation of the far field e leads to a large variation of the computed coefficients q. After F has reached full row rank, the condition number reduces.
e result of Figure 3 shows that, for solving equation (16) for a given number of directions where e is specified, the number of columns of F should be chosen such that the matrix has full column rank.

Interpolation.
For frequency points in between the simulated frequencies, one has to interpolate the spherical wave coefficients and determine the embedded element patterns from these interpolated coefficients.

International Journal of Antennas and Propagation
Suppose that the spherical wave coefficients Q n (f) are calculated in steps Δf between f min and f max (f min and f max included) which are both a multiple of the frequency step Δf, then f min and f max can be written as e number of frequencies is N 2 − N 1 + 1. If the bins of the FFT are numbered from 0 to M − 1, the values of Q n (f) are stored in the bins as shown in Table 4. e bins (M/2) to M − 1 represent the negative frequencies. In these bins, the complex conjugate of the coefficients (denoted by * in Table 4) at the corresponding positive frequency is stored. e value of M must be larger than 2N 2 + 2 and should be a power of 2 for an FFT. e spherical wave coefficients can be transformed to the time domain by using the inverse FFT: , the time domain version of the spherical wave coefficient q n (mΔt) is a real function. e length of the time domain window of the spherical wave coefficient is given by T � (1/Δf). e time step is given by Δt � 1/(MΔf). q n (mΔt) will decay to zero because the time domain response of an antenna element also decays to zero. If T is large enough, the values of q n (mΔt) can be neglected for larger values of m. e value of T, for which this happens, determines the largest frequency step Δf that can be used for interpolation; that is, the sampling in the frequency domain is too coarse and aliasing in the time domain will happen if T is large enough. Due to the aliasing in the time domain, q n (mΔt) does not decay to zero for large values of m.
If the value of T is large enough to neglect the effect of aliasing in the time domain, the spherical wave coefficients can be interpolated in the frequency domain by padding q n (mΔf) with zeros for values of m > M − 1 and using the FFT:    International Journal of Antennas and Propagation where M ′ is the new number of bins and Δf ′ is the new frequency step which is equal to Δf ′ � 1/(M ′ Δt). By choosing M ′ > M the frequency step, Δf ′ will be smaller than the original frequency step Δf. q n (mΔt) � 0 for values of m from M to M ′ − 1. So, in summary, the following procedure can be used to determine the values of M and M ′ : (i) M should be larger than 2N 2 + 2 and a power of 2.
is determines a minimum value of M.
(ii) Δf should be chosen small enough such that the time domain spherical wave coefficients q n (mΔt) decay to zero within the time window T � (1/Δf). (iii) If a finer resolution in the time domain is needed, the value of M can be increased further (Δt � 1/(MΔf). (iv) e value of M ′ is determined by the required frequency step after interpolation:

Test Case: The LOFAR System
3.1. General Description of the LOFAR Telescope. LOFAR (LOw Frequency ARray) is a radio telescope for radio astronomical observations in the frequency range of 10 MHz to 240 MHz. It consists of 52 stations distributed over Europe.
ere are three kinds of stations: core stations, remote stations, and international stations. e core stations and remote stations are located in the Netherlands. e international stations are located in France, Germany, Ireland, Latvia, Poland, Sweden, and the United Kingdom. Each station consists of two antenna types: low-band antennas (LBA) and high-band antennas (HBA). e low-band antennas are designed for the frequency range of 10-90 MHz. e high-band antennas are designed for the frequency range of 110-240 MHz. Each low-band antenna consists of two orthogonal inverted-V dipoles. e highband antennas consist of two orthogonal bowtie antennas. A group of 16 high-band antennas is combined in a unit called a tile. e number of LBA and HBA tiles for the different kinds of stations is shown in Table 5. More information about LOFAR and its science cases can be found in [9].
Each LOFAR station is rotated by a particular angle such that the side lobes of each station beam point towards a different point in the sky. is minimizes the effect of the side lobes on the sky image because the effect of the side lobes is averaged [11]. e orientation of the individual antennas with respect to a common coordinate system for all stations is the same. e result is that the configuration of each individual station is unique.

Detailed Description of the Low-Band Antenna.
e electromagnetic simulations are performed using FEKO [12]. e simulation model of a low-band antenna is shown in Figure 1.
Each low-band antenna consists of two inverted-V dipoles above a metal ground plane of 3×3 m 2 . In this paper, we will use the term "element" to denote one of these dipoles. e term "dual-pol element" will be used to denote a group of two elements as shown in Figure 1. Each arm of the two inverted-V dipoles has a length of 1.38 meters and is excited separately. e height of the feed point above-ground plane is 1.7 m. e dipole in the φ � 45°plane is referred to as polarization 1. e dipole in the φ � 135°plane is referred to as polarization 2. e outer conductor of the coaxial cable from the low-noise amplifier (LNA, positioned at the terminals marked by V rec 1,a , V rec 1,b , V rec 2,a , and V rec 2,b in Figure 1) to the metal ground plane is modeled by a vertical wire with a Table 4: Values of Q n (f) stored in bins.

Bin index
Frequency Condition number of F as a function of #columns for a given #rows International Journal of Antennas and Propagation length of 1.7 meters. e coupling between the outer conductor and the metal ground plane is modeled as a perfect connection. In practice, this is a capacitive coupling. e metal plate is modeled as an infinite thin perfect conducting sheet. e wires are modeled as perfect conducting wires with a radius of 0.1 mm. e metal ground plane is located at z � 0 mm. e space at z < 0 mm is modeled as a dielectric with a relative dielectric constant ε r of 3 and a conductivity σ of 0.01 S/m. is dielectric represents the soil. e space above z � 0 mm is modeled as a vacuum.
Each arm is excited by a voltage source. is means that each arm has its own embedded element pattern which is calculated by exciting its corresponding arm with 1 V while the other arms are short-circuited. ese embedded element patterns are transformed to voltages over the single-ended LNAs in receive situation for the case the antenna is excited with a plane wave with an electrical field strength of 1 V/m. is transformation is described in the appendix and follows an approach similar to [13]. Exciting each dipole arm separately instead of each dipole gives the possibility to investigate the common mode behaviour of the antenna.
is is however not discussed in this paper. If the single-ended voltages across each single-ended LNA in receive situation are equal to the voltages V rec 1,a , V rec 1,b , V rec 2,a , and V rec 2,b as shown in Figure 1, the open-circuit differential voltages of polarizations 1 and 2 can be written as where V rec diff,n with n � 1, 2 being the differential voltage for polarization n. Note that the load impedance used to calculate the single-ended voltages should be equal to the single-ended LNA impedance (half the differential mode impedance).

Description of Embedded Element Patterns.
In this section, we will discuss the amount of data that is needed to store all embedded element patterns of a LOFAR station and the possibility of reducing this. e dataset of the embedded element patterns of each dual-pol element consists of four complex numbers per direction per frequency point: a θ-and φ-component for each element (each dual-pol element consists of two orthogonal elements).
For this paper, we simulated the embedded element patterns from θ � 0°to θ � 90°in steps of 1°and from φ � 0°t o φ � 355°in steps of 5°. Table 6 shows the amount of complex numbers that is needed in this case for the different kind of LOFAR stations.
In general, the element patterns are not independent. is means that an element pattern of an arbitrary element can be written as a linear combination of the element patterns of other elements. If we can write the element patterns as a sum of orthogonal element patterns, the amount of data can be reduced because the number of orthogonal element patterns is less than the number of element patterns. To describe the element patterns only the coefficients of the linear combination have to be stored besides the orthogonal element patterns. e element patterns of an array of N elem elements can be written as a 2N dir × N elem matrix E: where e i is the element pattern as defined in equation (11) of element i. Each column of E represents an element pattern. Using a singular value decomposition, the matrix E can be written as where D is a 2N dir × N elem diagonal matrix containing the singular values sorted from large to small, U and W are unitary matrices of size 2N dir × 2N dir and N elem × N elem , respectively, and H denotes the Hermitian operator. e matrix E can be approximated by where D basis is a N basis × N basis diagonal matrix containing the N basis largest singular values. It is assumed that smaller singular values can be neglected. e matrices U basis and W basis are constructed by taking the first N basis columns of U and W and therefore have sizes 2N elem × N basis and N basis × N basis , respectively. e columns of U basis are a set of orthonormal basis functions for the element patterns. Figure 4 shows the number of dimensions N basis for the LOFAR low-band array as a function of frequency for frequencies from 10 to 90 MHz in steps of 5 MHz (17 frequencies) and 1 MHz (81 frequencies). It is assumed that singular values that are less than 0.01 times the largest singular value can be neglected.
Without using the singular value decomposition, one would need 2 × 96 × 17 �3264 embedded element patterns for 17 frequencies and 2 × 96 × 81 �15 552 embedded element patterns for 81 frequencies. If singular value decomposition is used, only 543 (the sum of the number of dominant dimensions N b for each frequency) element patterns are needed as a basis function for 17 frequency points. is means that only 543/3264�16.6% of the amount of data is needed to describe the element patterns for 17 frequency points. For 81 frequency points, these numbers are 2610 dimensions and 2610/15 552�16.8%, respectively.
is is a reduction of the number of embedded element  6 International Journal of Antennas and Propagation patterns that is needed to describe the low-band array. e following section demonstrates that the technique using spherical wave expansion, proposed in this work, reduces the data needed to describe each embedded element pattern.

Verification.
As a test case, the embedded element patterns of one of the core stations (CS302) are described using spherical wave expansion. e layout of the simulated station is shown in Figure 5. Element 1 is the center element. Elements 46 and 47 are the outer elements. e distance between these elements and the other ones is more than 20 meters. To apply spherical wave expansion, one has to determine N max in equation (9). In theory, the radius r 0 should be chosen such that the whole station will fit in a sphere with that radius. at would mean that r 0 should be equal to around 65 m. In practice, the embedded element pattern of an element is only determined by its closest neighbours. is means that r 0 can have a lower value than 65 m. For a given N max , r 0 is equal to It turned out that the coupling between an element for which the embedded element pattern is calculated and the elements further away than r 0 is less than − 35 dB for all simulated frequencies if N max � 31. is value is used for the calculations in this paper.
According to equation (8), N max � 31 corresponds with 2046 spherical wave coefficients. e patterns are calculated for the same θ-and φ-angles as in Section 3.3. e matrix F defined in (13) is calculated for θ � 1°t o θ � 179°in steps of 1°a nd from φ � 0°t o φ � 355°in steps of 5°. e normalized electric field e is set to zero for θ > 90°. e resulting size of F is 25 776 × 2046. Its condition number is 4.8, which shows that the sensitivity of the solution of (16) for deviations in the vector e is small. e rank of F is 2046 which means that the matrix has a full column rank. So, there is one solution for which the norm ‖e − Fq‖ 2 has its minimum. Figures 6 and 7 show the simulated embedded element pattern (V sim ) of element 47 (one of the two "isolated: elements) in the φ � 45°plane at 57 MHz together with the same embedded element pattern reconstructed from the calculated spherical mode coefficients (V recon ) and the difference between the simulated and reconstructed patterns (EES). e difference in amplitude is expressed as an Equivalent Error Signal (EES). e EES expressed in dB is defined as EES � 20 log 10 V recon − V sim .
(26) e frequency of 57 MHz is chosen because this is the resonance frequency of the inverted-V antenna used in the LBA. e embedded element patterns are plotted as the voltage at the input of the LNA for the case of an incoming plane wave with a field strength of 1 V/m.
One can see that the EES is less than − 30 dB and the phase difference is less than 0.5°in phase for most angles. e simulated and reconstructed embedded element patterns for element 1 in the φ � 45°plane at 57 MHz are shown in Figures 8 and 9. is element is the center element of the LBA station, within the densest part of the array. Due to the mutual coupling, the embedded element patterns have an irregular shape in this part of the array, especially around the resonance frequency of 57 MHz. So, the embedded element pattern of the element at 57 MHz can be considered as the worst case for the description of the embedded element patterns by spherical wave expansion.
It is clearly visible that the element pattern has a minimum around θ � 20°. is minimum is caused by the mutual coupling between the antennas and is also validated by a measurement [14]. e EES is less than − 20 dB. e phase difference is 7°around the amplitude dip at for θ � 20°a nd less than 2°for other angles. e embedded element patterns are described by 2046 coefficients (i.e., complex numbers) per simulated frequency.  Comparing this number with the amount of numbers needed to describe the embedded element patterns directly (in our case 90 θ-angles × 72 φ-angles × 2 field components � 12 960) does not make sense because this number can be chosen arbitrarily. A better number to compare with is the minimum amount of complex numbers needed to describe the embedded element patterns to fulfill the Nyquist criterion in the spatial domain. To determine this number, we interpolated the embedded element patterns to get data points on a regular u, v-grid. ese data points were transformed to the spatial frequency domain using a discrete Fourier transform. For the embedded element pattern of the center element at 57 MHz, it turned out that 2449 sample points were needed to have 99.9% of the energy included. Because every sample point contains two complex numbers (one for each field component), this means that 4898 complex numbers are needed to describe each embedded element pattern. ere were 2046 coefficients needed to describe the embedded element pattern. Because spherical wave expansion considers the total field (both field components), the amount of complex numbers to describe each embedded element pattern is the same. is means that by using spherical wave expansion 2046/4898 � 41.8% of the data is needed to describe the embedded element pattern at 57 MHz compared to direct storage of it sampled at the Nyquist rate.

Interpolation.
e LBA is simulated for frequencies between 10 and 90 MHz in discrete steps. In this subsection, we will show the application of the interpolation method described in Section 2.3. Figure 10 shows the time domain response (impulse response) of the electric field at 10 m distance at θ � 1°and φ � 45°if 33 frequency points between 10 and 90 MHz are used. is corresponds with a frequency resolution Δf of 2.5 MHz and a period T of 400 ns in the time domain. e size M of the used inverse FFT is 1024 points which corresponds with a time domain step of Δt � 1/(MΔf) � 1/(1024 · 2.5 · 10 6 ) Hz � 0.39 ns. It can be seen that the time domain response does not decay to zero. is means that the frequency step is too coarse. e time domain response for the case that the frequency step Δf is decreased to 1 MHz (81 frequency points) is shown in Figure 11. e corresponding length of the period T in the time domain is 1000 ns.
One can see that the time domain response decays to almost zero (less than 1% of the maximum amplitude). is means that a frequency resolution of 1 MHz is sufficient to interpolate the spherical wave coefficients. To illustrate the interpolation method, the frequency step is reduced to Δf ′ � 0.5 MHz. is means that the number of bins in the FFT has to be increased to M ′ � (MΔf/Δf ′ ) � (1024· 1 · 10 6 /0.5 · 10 6 ) � 2048. e spherical wave coefficients in the time domain are padded with 1024 zeros. Figures 12 and 13 show the amplitude and phase of the embedded element pattern at 57.5 MHz in the φ � 45°plane. V sim is the simulated pattern, V recon is the embedded element pattern obtained from the calculated spherical wave coefficients, and V interp,fft is the embedded element pattern obtained from the spherical wave coefficients obtained by FFT-interpolation. As a reference, the amplitude and phase of the embedded element patterns obtained by linear interpolation (V interp,linear ) and cubic spline interpolation (V interp,spline ) of the spherical wave coefficients are also included. Figure 14 shows several Equivalent Error Signals (EES). EES β α is the EES based on the difference between V β and V α , where β stands for recon, interp, fft, interp, linear, or interp, spline and α stands for sim or recon. Figure 15 shows the phase differences corresponding to the Equivalent Error Signals of Figure 14.
One can see that the differences between the interpolated and reconstructed embedded element patterns, on one hand, and the simulated embedded element pattern, on the other hand, are comparable. e EEP has a maximum of − 18 dB and the difference in phase has a maximum of 5°except around θ � 20°. e difference between the interpolated and reconstructed embedded element pattern is less (EEP less than − 30 dB, phase difference less than 1°(except around θ � 20°)). is means that the error due to the interpolation is less than the error due to the use of spherical wave expansion and it shows that the error due to the interpolation is not significant. So, using the FFT for interpolation eliminates the need to perform a full-wave electromagnetic simulation of the embedded element patterns with a frequency resolution of less than 1 MHz. From Figures 14 and  15, it is also clear that using the FFT for interpolation of the spherical wave coefficients gives a lower EES and phase difference compared to linear interpolation and cubic spline interpolation.
If the embedded element patterns between 10 and 90 MHz would be simulated with a frequency resolution of 0.5 MHz, a full-wave electromagnetic simulation of 161 frequency points would be needed. By using FFT-interpolation only 81 frequency points have to be simulated, which is a reduction of 80 frequency points. e cost for this reduction is to perform FFTs, which can be performed in less than 1 minute on a modern laptop computer. is is considerably less than a full-wave electromagnetic simulation of a complete LBA station at 80 additional frequency points on a dedicated server (which takes approximately 7 hours per frequency point).

Conclusions
is paper shows that it is possible to use spherical wave expansion to describe the embedded element patterns of the LOFAR low-band array by the coefficients of spherical wave  functions. e example showed that 41.8% of the data is needed to describe the embedded element patterns, compared to direct storage of the embedded element patterns sampled at the Nyquist rate. It turned out that the original embedded element patterns could be reconstructed with an EES of less than − 20 dB and a phase error better than 2°.
Only at the dip in the pattern, the accuracy is less. By using an FFT, it is possible to interpolate the embedded element patterns in the frequency domain. It turned out that a frequency resolution of 1 MHz is sufficient to interpolate the embedded element patterns. e main difference with respect to the simulated pattern is caused by the use of spherical wave expansion. It is also shown that the error due to interpolation by using the FFT is lower compared to the use of linear interpolation or cubic spline interpolation.

Transformation of the Embedded Element Patterns in Transmit Situation to Embedded Element Patterns in Receive Situation for a Given Load Impedance
e embedded element patterns of an array in receive situation depend on the load impedance terminating the array elements. e simulation software used for the calculation of the low-band array calculates the Z-matrix and the embedded element patterns in transmit situation using voltage sources. is means that, for the calculation of each embedded element pattern, one element is excited with a voltage source of 1 V while the other elements are shortcircuited.
ese embedded element patterns and the impedance matrix Z of the antenna will be used to calculate the voltages across the load impedances in receive situation. e advantage of this approach is that these voltages can be calculated in postprocessing, so the embedded element patterns for different load impedances can be calculated without rerunning the electromagnetic simulation. is appendix describes this postprocessing procedure. e embedded element patterns in transmit situation for the case that the array is excited with voltage sources can be written as an N ports × 2 matrix p y (u r ) where N ports is the number of elements: where e y i,θ (u r ) and e y i,φ (u r ) are θ-and φ-component of the normalized electric field in the direction u r as defined in (1) for the case element i is excited while the other ports are shorted. e total field p tot of the array can be written as where the vector V ports is an N ports -elements column vector containing the excitation voltages of each element. e row vector p tot (u r ) contains the normalized electric field in a direction u r : p tot u r � e tot,θ u r e tot,φ u r .
In a similar way, the embedded element patterns for the case the elements are excited with current sources (while the nonexcited elements are open) can be written in a matrix p z (u r ). e total electric field p tot (u r ) can now be written as where A → (u r ) is the (complex) amplitude of an incident plane wave with an electric field E →inc ( r → , u r ) at position r → coming from direction u r : and g → z,i (u r ) is the "voltage receive pattern" given by e relation between voltages and currents at the ports of an array antenna can now be written as V rec � ZI rec + g z u r A u r , (A.9) where g z (u r ) is an N ports × 2 matrix containing the voltage receive patterns of all array elements for a specific direction u r : g z u r � g z 1,θ u r g z 1,φ u r ⋮ ⋮ g z N ports ,θ u r g z N ports ,φ u r ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .
(A.10) e vector A(u r ) is a two-element column vector containing the complex amplitude of the θ-and φ-component of an incoming plane wave from direction u r : where A θ (u r ) and A φ (u r ) are the θ-and φ-components of A → (u r ): A → u r � A θ u r u θ + A φ u r u φ . (A.12) International Journal of Antennas and Propagation e relation between p z and g z is given by [15] g z � 4π jωμ 0 p z . (A. 13) If the array in receive situation is described in terms of incoming and reflected waves a rec and b rec , respectively, and a scattering matrix S, the reflected wave can be written as b rec � Sa rec + g s A, (A.14) where g s is an N ports × 2 matrix containing the "voltage wave receiving patterns" in a similar way to (A.10). e voltage wave receiving pattern is the reflected voltage wave for the case all the ports of the array are terminated with the reference impedance for which a rec and b rec are calculated (i.e., a rec � 0). e relation between the voltages V rec and I rec and the voltage waves a rec and b rec is given by [16] a rec � F V rec + Z 0 I rec , where Z 0 is a diagonal matrix containing the (complex) reference impedances for each antenna port. We will use the notation Z 0 � R 0 + jX 0 where R 0 is a diagonal matrix with the real part of the reference impedances and X 0 is a diagonal matrix with the imaginary part of the reference impedances. In this paper, it is assumed that the real part of the reference impedance is positive. e matrix F in equations (A.15) and (A.16) is equal to Data Availability e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.