# Frequency Interpolation of LOFAR Embedded Element Patterns Using Spherical Wave Expansion

Authors: M. J. Arts; D. S. Prinsloo; M. J. Bentum; A. B. Smolders
Journal: International Journal of Antennas and Propagation
 (2021)

---

## Abstract

This paper describes the use of spherical wave expansion (SWE) to model the embedded element patterns of the LOFAR low-band array. The goal is to reduce the amount of data needed to store the embedded element patterns. The coefficients are calculated using the Moore–Penrose pseudoinverse. The Fast Fourier Transform (FFT) is used to interpolate the coefficients in the frequency domain. It turned out that the embedded element patterns can be described by only 41.8% of the data needed to describe them directly if sampled at the Nyquist rate. The presented results show that a frequency resolution of 1 MHz is needed for proper interpolation of the spherical wave coefficients over the 80 MHz operating frequency band of the LOFAR low-band array. It is also shown that the error due to interpolation using the FFT is less than the error due to linear interpolation or cubic spline interpolation.


---

## Body

## 1. Introduction

For the calibration of phased arrays, a priori knowledge of the embedded element patterns (EEPs) is needed. These 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.The 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. This technique uses the Fast Fourier Transform (FFT) to interpolate the spherical wave coefficients in the frequency domain. These interpolated coefficients can be used to calculate the embedded element patterns at nonsimulated frequencies. The interpolation method can be used for phased arrays and Phased Array Feeds (PAFs). Application of PAF’s can be found in radio astronomy [5–7] and future 5G communication systems [8]. The LOFAR radio telescope [9] will be used as a test case to verify our proposed concept.The structure of this paper is as follows: Section2 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. The interpolation method is also discussed in this section. A description of the LOFAR radio telescope is given in Section 3. This section gives a general overview of the LOFAR system and a detailed description of the low-band antenna. The results are described in Section 4. Conclusions are given in Section 5.


## 2. Far-Field Spherical Wave Expansion

### 2.1. Spherical Wave Functions

The electric field component of the far field of an antenna can be written as(1)E⟶r,u^r=e⟶u^re−jβrr,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(2)u^r=u^xsinθcosφ+u^ysinθsinφ+u^zcosθ,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). The normalized electric far field can be expanded in terms of spherical wave functions:(3)e⟶u^r=∑n=1∞∑m=−nn∑s=12Qsmnf⟶smnu^r,where f⟶smnu^r is the spherical function s,m,n and Qsmn the corresponding coefficient. It is convenient to use a so-called “compressed index” p defined as [2, p. 15](4)p=2nn+1+m−1+s,to replace the indices s,m,n. The numbering scheme up to n=2 is shown in Table 1.Figure 1
Simulation model of a low-band antenna.Table 1
Numbering scheme of the compressed indexp.
p12345678910111213141516s1212121212121212n1111112222222222m−1−10011−2−2−1−1001122Using the compressed index, equation (3) can be written as(5)e⟶u^r=∑p=1∞Qpf⟶pu^r.Each spherical wave functionf⟶pu^r is normalized such that(6)∫φ=02π∫θ=0πf⟶pu^rθ,φ2sinθdθdφ=ZF,where ZF is the wave impedance in free space (≈377Ω). The spherical wave functions f⟶pu^rθ,φ and f⟶qu^rθ,φ are orthogonal with respect to each other:(7)∫φ=02π∫θ=0πf⟶pu^rθ,φ⋅f⟶q∗u^rθ,φsinθdθdφ=0,forp≠q.Note that we used the notationu^rθ,φ to denote that u^r depends on θ and φ according to equation (2). The complex conjugate is denoted by ∗.


### 2.2. Calculation of the CoefficientsQp

In practice, the summation overn in (3) is limited to a finite number Nmax. By setting n=m=Nmax and s=2 in (4), one can show that the maximum number of modes Nmodes should be at least(8)Nmodes=2NmaxNmax+1+Nmax−1+2=2NmaxNmax+2,Nmodes should be chosen such that the contribution of the higher-order modes (>Nmodes) can be neglected. As a rule of thumb, Nmax will be chosen as [2, p. 17](9)Nmax=βr0+10,which corresponds with(10)Nmodes=2βr0+10βr0+12=2β2r02+44βr0+240,where r0 is the radius of the smallest sphere enclosing the antenna.If the radiation patterne⟶u^r of an antenna is known in Ndir directions and we want to describe it in Nmodes spherical modes, it can be written as a vector e with length 2Ndir:(11)e=eθu^r,1⋮eθu^r,Ndireφu^r,1⋮eφu^r,Ndir,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:(12)e⟶u^r,n=eθu^r,nu^θ+eφu^r,nu^φ,where u^θ and u^φ are the unit vector in θ-and φ-direction, respectively. The far-field spherical vector functions for the first Nmodes modes specified in Ndir directions can be written in a 2Ndir×Nmodes matrix F:(13)F=fθ,1u^r,1…fθ,Nmodesu^r,1⋮⋱⋮fθ,1u^r,Ndir…fθ,Nmodesu^r,Ndirfφ,1u^r,1…fφ,Nmodesu^r,1⋮⋱⋮fφ,1u^r,Ndir…fφ,Nmodesu^r,Ndir.The scalarsfθ,pu^r and fφ,pu^r are the θ- and φ-component of the spherical wave function f⟶pu^r, respectively:(14)f⟶pu^r=fθ,pu^ru^θ+fφ,pu^ru^φ,where we used the compressed index p. The spherical mode coefficients Qp can be grouped in a vector q with length Nmodes:(15)q=Q1⋮QNmodes.Equation (5) can now be written as(16)e=Fq.The solution of equation (16) is given by [10](17)q=F†e+UNmodes−F†Fw,where the operator † represents the Moore–Penrose pseudoinverse, UNmodes is an Nmodes×Nmodes identity matrix, and w is a column vector with length Nmodes 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−Fq2 has its absolute minimum. If FF†=U2Ndir (with U2Ndir a 2Ndir×2Ndir identity matrix), equation (16) has at least one solution for arbitrary e. This 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 Qp≈0 for p>Nmodes, so the norm e−Fq2 will be very close to zero, and equation (17) can practically be considered as a solution of equation (16). If F†F=UNmodes, there exists only one solution or value for which the norm e−Fq2 has its absolute minimum. This 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≠UNmodes, there exist infinite solutions or values for which the norm e−Fq2 has its absolute minimum. These solutions or values are parametric on the vector w in equation (17).Figure2 shows the rank of matrix F as a function of the number of columns (coefficients) for a different number of rows (angles).Figure 2
The rank of the matrixF as a function of the number of columns.The number of columns represents the number of spherical wave coefficients. The number of rows is directly related to the number of directions for which the far field is given. This relation is shown in Table2. The 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). The angles θ=0° and θ=180° are not included because the spherical wave functions are not defined for these angles.Table 2
Relation between the number of rows and directions for which the far field is given.
Number of rowsNumber ofθ-anglesθ-rangeNumber ofφ-anglesφ-range288820° to 160° in steps of 20°180° to 340° in steps of 20°12241710° to 170° in steps of 10°360° to 350° in steps of 10°5040355° to 175° in steps of 5°720° to 355° in steps of 5°25 7761791° to 179° in steps of 1°720° to 355° in steps of 5°For an increasing number of columns (number of coefficients) for a given number of rows, the matrixF has full column rank until a certain number of columns. The highest number of columns where F has full column rank are denoted by ∗ in Figure 2. The 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. These 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.Table3 shows for which number of columns the matrix F has full row and column rank.Table 3
Values for which the matrixF has full row and column rank.
Number of rowsFull column rank up toFull row rank from288160 columns510 columns1224646 columns2310 columns50402590 columns—25 7762736 columns—Table3 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. This means that, for this case, equation (16) has one value for which the norm e−Fq2 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.Figure3 shows the condition number of the matrix F as a function of the number of columns for different numbers of rows. The 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.Figure 3
The condition number of the matrixF as a function of the number of columns. The horizontal black line denotes a condition number of 10.One can see that the condition number remains below 10 until the point whereF reaches full column rank. After this point, the condition number increases to values greater than 1015. 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.The result of Figure3 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.


### 2.3. 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.Suppose that the spherical wave coefficientsQnf are calculated in steps Δf between fmin and fmax (fmin and fmax included) which are both a multiple of the frequency step Δf, then fmin and fmax can be written as(18)fmin=N1Δf,fmax=N2Δf.The number of frequencies isN2−N1+1. If the bins of the FFT are numbered from 0 to M−1, the values of Qnf are stored in the bins as shown in Table 4.Table 4
Values ofQnf stored in bins.
Bin indexFrequencyStored value0001Δf0⋮⋮⋮N1−1N1−1Δf0N1N1ΔfQnN1Δf⋮⋮⋮N2N2ΔfQnN2ΔfN2+1N2+1Δf0⋮⋮⋮M/2−1M/2−1Δf0M/2−M/2Δf0⋮⋮⋮M−N2−1−N2+1Δf0M−N2−N2ΔfQn∗N2Δf⋮⋮⋮M−N1−N1ΔfQn∗N1ΔfM−N1+1−N1−1Δ0⋮⋮⋮M−1−Δf0The binsM/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. The value of M must be larger than 2N2+2 and should be a power of 2 for an FFT. The spherical wave coefficients can be transformed to the time domain by using the inverse FFT:(19)qnmΔt=1M∑k=0M−1QnkΔfej2πkm/M.BecauseQnM−kΔf=Qn−kΔf=Qn∗kΔf, the time domain version of the spherical wave coefficient qnmΔt is a real function.The length of the time domain window of the spherical wave coefficient is given byT=1/Δf. The time step is given by Δt=1/MΔf. qnmΔ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 qnmΔt can be neglected for larger values of m. The 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, qnmΔt does not decay to zero for large values of m.If the value ofT 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 qnmΔf with zeros for values of m>M−1 and using the FFT:(20)QnkΔf′=∑m=0M′−1qnmΔte−j2πkm/M′,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. qnmΔt=0 for values of m from M to M′−1.So, in summary, the following procedure can be used to determine the values ofM and M′:(i)
M should be larger than 2N2+2 and a power of 2. This determines a minimum value of M.(ii)
Δf should be chosen small enough such that the time domain spherical wave coefficients qnmΔt decay to zero within the time window T=1/Δf.(iii)
If a finer resolution in the time domain is needed, the value ofM can be increased further (Δt=1/MΔf.(iv)
The value ofM′ is determined by the required frequency step after interpolation: M′=1/Δf′Δt=MΔf/Δf′.


## 2.1. Spherical Wave Functions

The electric field component of the far field of an antenna can be written as(1)E⟶r,u^r=e⟶u^re−jβrr,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(2)u^r=u^xsinθcosφ+u^ysinθsinφ+u^zcosθ,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). The normalized electric far field can be expanded in terms of spherical wave functions:(3)e⟶u^r=∑n=1∞∑m=−nn∑s=12Qsmnf⟶smnu^r,where f⟶smnu^r is the spherical function s,m,n and Qsmn the corresponding coefficient. It is convenient to use a so-called “compressed index” p defined as [2, p. 15](4)p=2nn+1+m−1+s,to replace the indices s,m,n. The numbering scheme up to n=2 is shown in Table 1.Figure 1
Simulation model of a low-band antenna.Table 1
Numbering scheme of the compressed indexp.
p12345678910111213141516s1212121212121212n1111112222222222m−1−10011−2−2−1−1001122Using the compressed index, equation (3) can be written as(5)e⟶u^r=∑p=1∞Qpf⟶pu^r.Each spherical wave functionf⟶pu^r is normalized such that(6)∫φ=02π∫θ=0πf⟶pu^rθ,φ2sinθdθdφ=ZF,where ZF is the wave impedance in free space (≈377Ω). The spherical wave functions f⟶pu^rθ,φ and f⟶qu^rθ,φ are orthogonal with respect to each other:(7)∫φ=02π∫θ=0πf⟶pu^rθ,φ⋅f⟶q∗u^rθ,φsinθdθdφ=0,forp≠q.Note that we used the notationu^rθ,φ to denote that u^r depends on θ and φ according to equation (2). The complex conjugate is denoted by ∗.


## 2.2. Calculation of the CoefficientsQp

In practice, the summation overn in (3) is limited to a finite number Nmax. By setting n=m=Nmax and s=2 in (4), one can show that the maximum number of modes Nmodes should be at least(8)Nmodes=2NmaxNmax+1+Nmax−1+2=2NmaxNmax+2,Nmodes should be chosen such that the contribution of the higher-order modes (>Nmodes) can be neglected. As a rule of thumb, Nmax will be chosen as [2, p. 17](9)Nmax=βr0+10,which corresponds with(10)Nmodes=2βr0+10βr0+12=2β2r02+44βr0+240,where r0 is the radius of the smallest sphere enclosing the antenna.If the radiation patterne⟶u^r of an antenna is known in Ndir directions and we want to describe it in Nmodes spherical modes, it can be written as a vector e with length 2Ndir:(11)e=eθu^r,1⋮eθu^r,Ndireφu^r,1⋮eφu^r,Ndir,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:(12)e⟶u^r,n=eθu^r,nu^θ+eφu^r,nu^φ,where u^θ and u^φ are the unit vector in θ-and φ-direction, respectively. The far-field spherical vector functions for the first Nmodes modes specified in Ndir directions can be written in a 2Ndir×Nmodes matrix F:(13)F=fθ,1u^r,1…fθ,Nmodesu^r,1⋮⋱⋮fθ,1u^r,Ndir…fθ,Nmodesu^r,Ndirfφ,1u^r,1…fφ,Nmodesu^r,1⋮⋱⋮fφ,1u^r,Ndir…fφ,Nmodesu^r,Ndir.The scalarsfθ,pu^r and fφ,pu^r are the θ- and φ-component of the spherical wave function f⟶pu^r, respectively:(14)f⟶pu^r=fθ,pu^ru^θ+fφ,pu^ru^φ,where we used the compressed index p. The spherical mode coefficients Qp can be grouped in a vector q with length Nmodes:(15)q=Q1⋮QNmodes.Equation (5) can now be written as(16)e=Fq.The solution of equation (16) is given by [10](17)q=F†e+UNmodes−F†Fw,where the operator † represents the Moore–Penrose pseudoinverse, UNmodes is an Nmodes×Nmodes identity matrix, and w is a column vector with length Nmodes 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−Fq2 has its absolute minimum. If FF†=U2Ndir (with U2Ndir a 2Ndir×2Ndir identity matrix), equation (16) has at least one solution for arbitrary e. This 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 Qp≈0 for p>Nmodes, so the norm e−Fq2 will be very close to zero, and equation (17) can practically be considered as a solution of equation (16). If F†F=UNmodes, there exists only one solution or value for which the norm e−Fq2 has its absolute minimum. This 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≠UNmodes, there exist infinite solutions or values for which the norm e−Fq2 has its absolute minimum. These solutions or values are parametric on the vector w in equation (17).Figure2 shows the rank of matrix F as a function of the number of columns (coefficients) for a different number of rows (angles).Figure 2
The rank of the matrixF as a function of the number of columns.The number of columns represents the number of spherical wave coefficients. The number of rows is directly related to the number of directions for which the far field is given. This relation is shown in Table2. The 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). The angles θ=0° and θ=180° are not included because the spherical wave functions are not defined for these angles.Table 2
Relation between the number of rows and directions for which the far field is given.
Number of rowsNumber ofθ-anglesθ-rangeNumber ofφ-anglesφ-range288820° to 160° in steps of 20°180° to 340° in steps of 20°12241710° to 170° in steps of 10°360° to 350° in steps of 10°5040355° to 175° in steps of 5°720° to 355° in steps of 5°25 7761791° to 179° in steps of 1°720° to 355° in steps of 5°For an increasing number of columns (number of coefficients) for a given number of rows, the matrixF has full column rank until a certain number of columns. The highest number of columns where F has full column rank are denoted by ∗ in Figure 2. The 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. These 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.Table3 shows for which number of columns the matrix F has full row and column rank.Table 3
Values for which the matrixF has full row and column rank.
Number of rowsFull column rank up toFull row rank from288160 columns510 columns1224646 columns2310 columns50402590 columns—25 7762736 columns—Table3 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. This means that, for this case, equation (16) has one value for which the norm e−Fq2 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.Figure3 shows the condition number of the matrix F as a function of the number of columns for different numbers of rows. The 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.Figure 3
The condition number of the matrixF as a function of the number of columns. The horizontal black line denotes a condition number of 10.One can see that the condition number remains below 10 until the point whereF reaches full column rank. After this point, the condition number increases to values greater than 1015. 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.The result of Figure3 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.


## 2.3. 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.Suppose that the spherical wave coefficientsQnf are calculated in steps Δf between fmin and fmax (fmin and fmax included) which are both a multiple of the frequency step Δf, then fmin and fmax can be written as(18)fmin=N1Δf,fmax=N2Δf.The number of frequencies isN2−N1+1. If the bins of the FFT are numbered from 0 to M−1, the values of Qnf are stored in the bins as shown in Table 4.Table 4
Values ofQnf stored in bins.
Bin indexFrequencyStored value0001Δf0⋮⋮⋮N1−1N1−1Δf0N1N1ΔfQnN1Δf⋮⋮⋮N2N2ΔfQnN2ΔfN2+1N2+1Δf0⋮⋮⋮M/2−1M/2−1Δf0M/2−M/2Δf0⋮⋮⋮M−N2−1−N2+1Δf0M−N2−N2ΔfQn∗N2Δf⋮⋮⋮M−N1−N1ΔfQn∗N1ΔfM−N1+1−N1−1Δ0⋮⋮⋮M−1−Δf0The binsM/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. The value of M must be larger than 2N2+2 and should be a power of 2 for an FFT. The spherical wave coefficients can be transformed to the time domain by using the inverse FFT:(19)qnmΔt=1M∑k=0M−1QnkΔfej2πkm/M.BecauseQnM−kΔf=Qn−kΔf=Qn∗kΔf, the time domain version of the spherical wave coefficient qnmΔt is a real function.The length of the time domain window of the spherical wave coefficient is given byT=1/Δf. The time step is given by Δt=1/MΔf. qnmΔ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 qnmΔt can be neglected for larger values of m. The 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, qnmΔt does not decay to zero for large values of m.If the value ofT 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 qnmΔf with zeros for values of m>M−1 and using the FFT:(20)QnkΔf′=∑m=0M′−1qnmΔte−j2πkm/M′,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. qnmΔt=0 for values of m from M to M′−1.So, in summary, the following procedure can be used to determine the values ofM and M′:(i)
M should be larger than 2N2+2 and a power of 2. This determines a minimum value of M.(ii)
Δf should be chosen small enough such that the time domain spherical wave coefficients qnmΔt decay to zero within the time window T=1/Δf.(iii)
If a finer resolution in the time domain is needed, the value ofM can be increased further (Δt=1/MΔf.(iv)
The value ofM′ is determined by the required frequency step after interpolation: M′=1/Δf′Δt=MΔf/Δf′.


## 3. 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. There are three kinds of stations: core stations, remote stations, and international stations. The core stations and remote stations are located in the Netherlands. The 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). The low-band antennas are designed for the frequency range of 10–90 MHz. The high-band antennas are designed for the frequency range of 110–240 MHz. Each low-band antenna consists of two orthogonal inverted-V dipoles. The high-band antennas consist of two orthogonal bowtie antennas. A group of 16 high-band antennas is combined in a unit called a tile. The number of LBA and HBA tiles for the different kinds of stations is shown in Table5. More information about LOFAR and its science cases can be found in [9].Table 5
Number of LBA and HBA tiles for the different kinds of stations.
Kind of stationNumber of LBANumber of HBA tilesCore9648 (2 clusters of 24 tiles)Remote9648 (1 cluster of 48 tiles)International9696 (1 cluster of 96 tiles)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. This minimizes the effect of the side lobes on the sky image because the effect of the side lobes is averaged [11]. The orientation of the individual antennas with respect to a common coordinate system for all stations is the same. The result is that the configuration of each individual station is unique.


### 3.2. Detailed Description of the Low-Band Antenna

The electromagnetic simulations are performed using FEKO [12]. The 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 m2. In this paper, we will use the term “element” to denote one of these dipoles. The 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. The height of the feed point above-ground plane is 1.7 m. The dipole in the φ=45° plane is referred to as polarization 1. The dipole in the φ=135° plane is referred to as polarization 2. The outer conductor of the coaxial cable from the low-noise amplifier (LNA, positioned at the terminals marked by V1,arec, V1,brec, V2,arec, and V2,brec in Figure 1) to the metal ground plane is modeled by a vertical wire with a length of 1.7 meters. The coupling between the outer conductor and the metal ground plane is modeled as a perfect connection. In practice, this is a capacitive coupling. The metal plate is modeled as an infinite thin perfect conducting sheet. The wires are modeled as perfect conducting wires with a radius of 0.1 mm. The metal ground plane is located at z=0 mm. The 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. This dielectric represents the soil. The space above z=0 mm is modeled as a vacuum.Each arm is excited by a voltage source. This 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. These 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. This 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. This is however not discussed in this paper.If the single-ended voltages across each single-ended LNA in receive situation are equal to the voltagesV1,arec, V1,brec, V2,arec, and V2,brec as shown in Figure 1, the open-circuit differential voltages of polarizations 1 and 2 can be written as(21)Vdiff,1rec=V1,arec−V1,brec,Vdiff,2rec=V2,arec−V2,brec,where Vdiff,nrec 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).


### 3.3. 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. The 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° to φ=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.Table 6
Number of complex numbers needed to describe element patterns of one station.
Kind of stationAntenna typeNumber of dual-pol elementsNumber of complex numbers per frequency pointCore and remoteLBA962 515 968Core and remoteHBA76820 127 744InternationalLBA962 515 968InternationalHBA1 53640 255 488In general, the element patterns are not independent. This 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. The element patterns of an array ofNelem elements can be written as a 2Ndir×Nelem matrix E:(22)E=e1…eNelem,where ei 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 matrixE can be written as(23)E=UDWH,where D is a 2Ndir×Nelem diagonal matrix containing the singular values sorted from large to small, U and W are unitary matrices of size 2Ndir×2Ndir and Nelem×Nelem, respectively, and H denotes the Hermitian operator. The matrix E can be approximated by(24)E≈UbasisDbasisWbasisH,where Dbasis is a Nbasis×Nbasis diagonal matrix containing the Nbasis largest singular values. It is assumed that smaller singular values can be neglected. The matrices Ubasis and Wbasis are constructed by taking the first Nbasis columns of U and W and therefore have sizes 2Nelem×Nbasis and Nbasis×Nbasis, respectively. The columns of Ubasis are a set of orthonormal basis functions for the element patterns.Figure4 shows the number of dimensions Nbasis 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.Figure 4
Number of dominant dimensions of the embedded element patterns of the LOFAR low-band array as a function of frequency.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 Nb for each frequency) element patterns are needed as a basis function for 17 frequency points. This 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. This is a reduction of the number of embedded element patterns that is needed to describe the low-band array. The following section demonstrates that the technique using spherical wave expansion, proposed in this work, reduces the data needed to describe each embedded element pattern.


## 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. There are three kinds of stations: core stations, remote stations, and international stations. The core stations and remote stations are located in the Netherlands. The 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). The low-band antennas are designed for the frequency range of 10–90 MHz. The high-band antennas are designed for the frequency range of 110–240 MHz. Each low-band antenna consists of two orthogonal inverted-V dipoles. The high-band antennas consist of two orthogonal bowtie antennas. A group of 16 high-band antennas is combined in a unit called a tile. The number of LBA and HBA tiles for the different kinds of stations is shown in Table5. More information about LOFAR and its science cases can be found in [9].Table 5
Number of LBA and HBA tiles for the different kinds of stations.
Kind of stationNumber of LBANumber of HBA tilesCore9648 (2 clusters of 24 tiles)Remote9648 (1 cluster of 48 tiles)International9696 (1 cluster of 96 tiles)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. This minimizes the effect of the side lobes on the sky image because the effect of the side lobes is averaged [11]. The orientation of the individual antennas with respect to a common coordinate system for all stations is the same. The result is that the configuration of each individual station is unique.


## 3.2. Detailed Description of the Low-Band Antenna

The electromagnetic simulations are performed using FEKO [12]. The 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 m2. In this paper, we will use the term “element” to denote one of these dipoles. The 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. The height of the feed point above-ground plane is 1.7 m. The dipole in the φ=45° plane is referred to as polarization 1. The dipole in the φ=135° plane is referred to as polarization 2. The outer conductor of the coaxial cable from the low-noise amplifier (LNA, positioned at the terminals marked by V1,arec, V1,brec, V2,arec, and V2,brec in Figure 1) to the metal ground plane is modeled by a vertical wire with a length of 1.7 meters. The coupling between the outer conductor and the metal ground plane is modeled as a perfect connection. In practice, this is a capacitive coupling. The metal plate is modeled as an infinite thin perfect conducting sheet. The wires are modeled as perfect conducting wires with a radius of 0.1 mm. The metal ground plane is located at z=0 mm. The 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. This dielectric represents the soil. The space above z=0 mm is modeled as a vacuum.Each arm is excited by a voltage source. This 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. These 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. This 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. This is however not discussed in this paper.If the single-ended voltages across each single-ended LNA in receive situation are equal to the voltagesV1,arec, V1,brec, V2,arec, and V2,brec as shown in Figure 1, the open-circuit differential voltages of polarizations 1 and 2 can be written as(21)Vdiff,1rec=V1,arec−V1,brec,Vdiff,2rec=V2,arec−V2,brec,where Vdiff,nrec 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).


## 3.3. 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. The 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° to φ=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.Table 6
Number of complex numbers needed to describe element patterns of one station.
Kind of stationAntenna typeNumber of dual-pol elementsNumber of complex numbers per frequency pointCore and remoteLBA962 515 968Core and remoteHBA76820 127 744InternationalLBA962 515 968InternationalHBA1 53640 255 488In general, the element patterns are not independent. This 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. The element patterns of an array ofNelem elements can be written as a 2Ndir×Nelem matrix E:(22)E=e1…eNelem,where ei 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 matrixE can be written as(23)E=UDWH,where D is a 2Ndir×Nelem diagonal matrix containing the singular values sorted from large to small, U and W are unitary matrices of size 2Ndir×2Ndir and Nelem×Nelem, respectively, and H denotes the Hermitian operator. The matrix E can be approximated by(24)E≈UbasisDbasisWbasisH,where Dbasis is a Nbasis×Nbasis diagonal matrix containing the Nbasis largest singular values. It is assumed that smaller singular values can be neglected. The matrices Ubasis and Wbasis are constructed by taking the first Nbasis columns of U and W and therefore have sizes 2Nelem×Nbasis and Nbasis×Nbasis, respectively. The columns of Ubasis are a set of orthonormal basis functions for the element patterns.Figure4 shows the number of dimensions Nbasis 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.Figure 4
Number of dominant dimensions of the embedded element patterns of the LOFAR low-band array as a function of frequency.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 Nb for each frequency) element patterns are needed as a basis function for 17 frequency points. This 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. This is a reduction of the number of embedded element patterns that is needed to describe the low-band array. The following section demonstrates that the technique using spherical wave expansion, proposed in this work, reduces the data needed to describe each embedded element pattern.


## 4. Results

### 4.1. Verification

As a test case, the embedded element patterns of one of the core stations (CS302) are described using spherical wave expansion. The layout of the simulated station is shown in Figure5. Element 1 is the center element. Elements 46 and 47 are the outer elements. The distance between these elements and the other ones is more than 20 meters.Figure 5
Layout of the low-band array of the simulated LOFAR station (CS302).To apply spherical wave expansion, one has to determineNmax in equation (9). In theory, the radius r0 should be chosen such that the whole station will fit in a sphere with that radius. That would mean that r0 should be equal to around 65 m. In practice, the embedded element pattern of an element is only determined by its closest neighbours. This means that r0 can have a lower value than 65 m. For a given Nmax, r0 is equal to(25)r0=Nmax−10β.It turned out that the coupling between an element for which the embedded element pattern is calculated and the elements further away thanr0 is less than −35 dB for all simulated frequencies if Nmax=31. This value is used for the calculations in this paper.According to equation (8), Nmax=31 corresponds with 2046 spherical wave coefficients. The patterns are calculated for the same θ- and φ-angles as in Section 3.3. The matrix F defined in (13) is calculated for θ=1° to θ=179° in steps of 1° and from φ=0° to φ=355° in steps of 5°. The normalized electric field e is set to zero for θ>90°. The 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. The 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−Fq2 has its minimum.Figures6 and 7 show the simulated embedded element pattern (Vsim) 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 (Vrecon) and the difference between the simulated and reconstructed patterns (EES). The difference in amplitude is expressed as an Equivalent Error Signal (EES). The EES expressed in dB is defined as(26)EES=20log10Vrecon−Vsim.Figure 6
Simulated and reconstructed EEP of element 47 (amplitude).Figure 7
Simulated and reconstructed EEP of element 47 (phase).The frequency of 57 MHz is chosen because this is the resonance frequency of the inverted-V antenna used in the LBA. The 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 than0.5° in phase for most angles.The simulated and reconstructed embedded element patterns for element 1 in theφ=45° plane at 57 MHz are shown in Figures 8 and 9. This 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.Figure 8
Simulated and reconstructed EEP of element 1 (amplitude).Figure 9
Simulated and reconstructed EEP of element 1 (phase).It is clearly visible that the element pattern has a minimum aroundθ=20°. This minimum is caused by the mutual coupling between the antennas and is also validated by a measurement [14]. The EES is less than −20 dB. The phase difference is 7° around the amplitude dip at for θ=20° and less than 2° for other angles.The 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. These 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. There 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. This 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.


### 4.2. Interpolation

The 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 Section2.3.Figure10 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. This corresponds with a frequency resolution Δf of 2.5 MHz and a period T of 400 ns in the time domain. The 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⋅106Hz=0.39 ns. It can be seen that the time domain response does not decay to zero. This means that the frequency step is too coarse.Figure 10
Impulse response of the E-field using 33 frequency points.The time domain response for the case that the frequency stepΔf is decreased to 1 MHz (81 frequency points) is shown in Figure 11. The corresponding length of the period T in the time domain is 1000 ns.Figure 11
Impulse response of the E-field using 81 frequency points.One can see that the time domain response decays to almost zero (less than 1% of the maximum amplitude). This 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. This means that the number of bins in the FFT has to be increased to M′=MΔf/Δf′=1024⋅1⋅106/0.5⋅106=2048. The spherical wave coefficients in the time domain are padded with 1024 zeros.Figures12 and 13 show the amplitude and phase of the embedded element pattern at 57.5 MHz in the φ=45° plane. Vsim is the simulated pattern, Vrecon is the embedded element pattern obtained from the calculated spherical wave coefficients, and Vinterp,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 (Vinterp,linear) and cubic spline interpolation (Vinterp,spline) of the spherical wave coefficients are also included.Figure 12
Simulated, reconstructed, and interpolated EEP of element 1 at 57.5 MHz (amplitude).Figure 13
Simulated, reconstructed, and interpolated EEP of element 1 at 57.5 MHz (phase).Figure14 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 14
EES of the EEP of element 1 at 57.5 MHz (amplitude).Figure15 shows the phase differences corresponding to the Equivalent Error Signals of Figure 14.Figure 15
Phase differences of EEP of element 1 at 57.5 MHz (phase).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. The EEP has a maximum of −18 dB and the difference in phase has a maximum of5° except around θ=20°. The difference between the interpolated and reconstructed embedded element pattern is less (EEP less than −30 dB, phase difference less than 1° (except around θ=20°)). This 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. The cost for this reduction is to perform FFTs, which can be performed in less than 1 minute on a modern laptop computer. This 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).


## 4.1. Verification

As a test case, the embedded element patterns of one of the core stations (CS302) are described using spherical wave expansion. The layout of the simulated station is shown in Figure5. Element 1 is the center element. Elements 46 and 47 are the outer elements. The distance between these elements and the other ones is more than 20 meters.Figure 5
Layout of the low-band array of the simulated LOFAR station (CS302).To apply spherical wave expansion, one has to determineNmax in equation (9). In theory, the radius r0 should be chosen such that the whole station will fit in a sphere with that radius. That would mean that r0 should be equal to around 65 m. In practice, the embedded element pattern of an element is only determined by its closest neighbours. This means that r0 can have a lower value than 65 m. For a given Nmax, r0 is equal to(25)r0=Nmax−10β.It turned out that the coupling between an element for which the embedded element pattern is calculated and the elements further away thanr0 is less than −35 dB for all simulated frequencies if Nmax=31. This value is used for the calculations in this paper.According to equation (8), Nmax=31 corresponds with 2046 spherical wave coefficients. The patterns are calculated for the same θ- and φ-angles as in Section 3.3. The matrix F defined in (13) is calculated for θ=1° to θ=179° in steps of 1° and from φ=0° to φ=355° in steps of 5°. The normalized electric field e is set to zero for θ>90°. The 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. The 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−Fq2 has its minimum.Figures6 and 7 show the simulated embedded element pattern (Vsim) 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 (Vrecon) and the difference between the simulated and reconstructed patterns (EES). The difference in amplitude is expressed as an Equivalent Error Signal (EES). The EES expressed in dB is defined as(26)EES=20log10Vrecon−Vsim.Figure 6
Simulated and reconstructed EEP of element 47 (amplitude).Figure 7
Simulated and reconstructed EEP of element 47 (phase).The frequency of 57 MHz is chosen because this is the resonance frequency of the inverted-V antenna used in the LBA. The 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 than0.5° in phase for most angles.The simulated and reconstructed embedded element patterns for element 1 in theφ=45° plane at 57 MHz are shown in Figures 8 and 9. This 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.Figure 8
Simulated and reconstructed EEP of element 1 (amplitude).Figure 9
Simulated and reconstructed EEP of element 1 (phase).It is clearly visible that the element pattern has a minimum aroundθ=20°. This minimum is caused by the mutual coupling between the antennas and is also validated by a measurement [14]. The EES is less than −20 dB. The phase difference is 7° around the amplitude dip at for θ=20° and less than 2° for other angles.The 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. These 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. There 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. This 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.


## 4.2. Interpolation

The 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 Section2.3.Figure10 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. This corresponds with a frequency resolution Δf of 2.5 MHz and a period T of 400 ns in the time domain. The 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⋅106Hz=0.39 ns. It can be seen that the time domain response does not decay to zero. This means that the frequency step is too coarse.Figure 10
Impulse response of the E-field using 33 frequency points.The time domain response for the case that the frequency stepΔf is decreased to 1 MHz (81 frequency points) is shown in Figure 11. The corresponding length of the period T in the time domain is 1000 ns.Figure 11
Impulse response of the E-field using 81 frequency points.One can see that the time domain response decays to almost zero (less than 1% of the maximum amplitude). This 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. This means that the number of bins in the FFT has to be increased to M′=MΔf/Δf′=1024⋅1⋅106/0.5⋅106=2048. The spherical wave coefficients in the time domain are padded with 1024 zeros.Figures12 and 13 show the amplitude and phase of the embedded element pattern at 57.5 MHz in the φ=45° plane. Vsim is the simulated pattern, Vrecon is the embedded element pattern obtained from the calculated spherical wave coefficients, and Vinterp,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 (Vinterp,linear) and cubic spline interpolation (Vinterp,spline) of the spherical wave coefficients are also included.Figure 12
Simulated, reconstructed, and interpolated EEP of element 1 at 57.5 MHz (amplitude).Figure 13
Simulated, reconstructed, and interpolated EEP of element 1 at 57.5 MHz (phase).Figure14 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 14
EES of the EEP of element 1 at 57.5 MHz (amplitude).Figure15 shows the phase differences corresponding to the Equivalent Error Signals of Figure 14.Figure 15
Phase differences of EEP of element 1 at 57.5 MHz (phase).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. The EEP has a maximum of −18 dB and the difference in phase has a maximum of5° except around θ=20°. The difference between the interpolated and reconstructed embedded element pattern is less (EEP less than −30 dB, phase difference less than 1° (except around θ=20°)). This 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. The cost for this reduction is to perform FFTs, which can be performed in less than 1 minute on a modern laptop computer. This 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).


## 5. Conclusions

This 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. The 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 than2°. 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. The 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.



---

*Source: 5598380-2021-06-16.xml*