Bayesian noise wave calibration for 21-cm global experiments

The noise necessitating calibration emerges during measurement-taking. In an averaged or global experiment, the sky temperature $T_{\mathrm{sky}}(\Omega,\nu,t)$ is a function of the direction $\Omega$, frequency $\nu$ and time $t$. This can be broken down into two primary components: the global 21-cm signal $T_{21}$ and astrophysical foregrounds $T_{\mathrm{f}}$

The antenna measures the sky signal convolved with the normalised antenna directivity $B$. The process of measurement introduces the random noise term $N_{\mathrm{data}}$.

Our 21-cm signature can thus be represented as

Here, the integral is assessed through foreground and beam modelling techniques such as those discussed in Anstey et al. (2020) while modelling of $N_{\mathrm{data}}$ from the statistical properties of $D(\nu,t)$ is accomplished by a calibration algorithm as articulated in this paper and outlined in Fig. 1. Having a fully Bayesian framework when modelling the beam, the sky and the systematics has major advantages for global 21-cm experiments such as REACH (de Lera Acedo et al., 2021), as it provides the greatest flexibility in being able to model all effects and jointly fit for them.

The standard calibration strategy follows the method introduced by Dicke to characterise systematic features in radio frequency instruments (Dicke, 1946) and is widely used in experiments such as EDGES (Monsalve et al., 2017) and LOFAR (Bilous, A. V. et al., 2016) to evaluate the spectral index of the sky’s diffuse radio background (Rogers & Bowman, 2012). This technique involves measurements of two internal reference standards; a load and a noise source, in addition to a series of external calibration sources attached to the receiver input in lieu of the antenna. These include an ambient-temperature ‘cold’ load, a ‘hot’ load heated to $\sim 400$ K, an open-ended cable and a shorted cable. A block diagram showing this arrangement is presented in Fig. 2.

When calibrating the receiver, reflection coefficients are taken of the calibration source connected to the receiver input ($\Gamma_{\mathrm{cal}}$) and of the receiver itself ($\Gamma_{\mathrm{rec}}$) as well as power spectral densities (PSDs) of the input ($P_{\mathrm{cal}}$), the internal reference load ($P_{\mathrm{L}}$) and the internal reference noise source ($P_{\mathrm{NS}}$) (Monsalve et al., 2017). These measurements are used to calculate a preliminary ‘uncalibrated’ antenna temperature $T_{\mathrm{cal}}^{*}$

where $T_{\mathrm{L}}$ and $T_{\mathrm{NS}}$ are assumptions for the noise temperature of the internal reference load and excess noise temperature of the internal noise source above ambient, respectively. This initial calculation is used to calibrate out any time-dependent system gain that emerges from a series of filters, amplifiers and cables, as well as the analogue-to-digital converter within the experimental apparatus (Monsalve et al., 2017). Each PSD measurement can be expressed in terms of specific response contributions as detailed in Bowman et al. (2018)

Here, $g_{\mathrm{sys}}$ is the system gain referenced to the receiver input and $T_{\mathrm{cal}}$ is our calibrated input temperature. $T_{\mathrm{unc}}$, $T_{\mathrm{cos}}$, and $T_{\mathrm{sin}}$ are the ‘noise wave parameters’ introduced by Meys (1978) to calibrate the instrument. $T_{\mathrm{unc}}$ represents the portion of noise reflected by the antenna that is uncorrelated with the output noise of the low noise amplifier (LNA). $T_{\mathrm{cos}}$ and $T_{\mathrm{sin}}$ are the portions of reflected noise correlated with noise from the LNA (Monsalve et al., 2017; Rogers & Bowman, 2012). In the EDGES experiment, these calibration quantities are modelled using seven-term polynomials in frequency.

The PSDs for the internal reference load and noise source can similarly be expressed as in equation 6. However, since the reflection coefficients of the internal references are typically less than 0.005, they are taken to be zero in order to simplify the equations

As shown in Fig. 2, the internal references may be on a separate reference plane than the receiver input, resulting in a system gain $g_{\mathrm{sys}}^{*}$ and a noise offset $T_{0}^{*}$ different from those defined in equation 6. This effect is taken into account by two additional scale and offset parameters, $C_{1}$ and $C_{2}$, introduced by EDGES (Monsalve et al., 2017).

Since $C_{1}$ and $C_{2}$ also correct for first-order assumptions in the noise temperatures of the internal reference load and noise source, we have chosen to absorb these terms into $T_{\mathrm{L}}$ and $T_{\mathrm{NS}}$. This adjustment allows all calibration parameters, $T_{\mathrm{unc}}$, $T_{\mathrm{cos}}$, $T_{\mathrm{sin}}$, and an ‘effective’ $T_{\mathrm{NS}}$ and $T_{\mathrm{L}}$, to be solved for in units of kelvin, facilitating a joint solution of parameters. Expanding equation 5 using equations 6, 7 and 8 yields a linear identity providing a relationship between the uncalibrated input temperature and a final calibrated temperature of any device connected to the receiver input

where all parameters are frequency-dependent. This is not explicitly shown for simplicity of notation. For estimation of the noise wave parameters, $T_{\mathrm{cal}}$, $\Gamma_{\mathrm{cal}}$ and $\Gamma_{\mathrm{rec}}$ are measured along with the PSDs while $g_{\mathrm{sys}}$ and $T_{\mathrm{0}}$ are calibrated out. The cold and hot loads exhibit the main temperature references needed for $T_{\mathrm{L}}$ and $T_{\mathrm{NS}}$. The cables facilitate the derivation of the noise wave parameters describing spectral ripples from the noise properties of the receiver by acting as antennas looking at an isotropic sky with temperatures equal to the cables’ physical temperatures (Rogers & Bowman, 2012).

One possible source of systematics in the calibration methodology used by EDGES comes from measuring the response of the four external calibrators along with the receiver reflection coefficient in a laboratory away from where the instrument is actually deployed (Bowman et al., 2018). This process, especially with regards to how calibration parameters change, can be non-trivial. Furthermore, the fixed polynomial order used by EDGES for all noise wave parameters may underfit or overfit individual parameters and thus ‘fit out’ data useful for determining systematics or potentially even the 21-cm signal itself if a joint fit is performed.

In response to these issues, we have developed a calibration pipeline that improves on the strategies presented in section 2.2. We introduce a novel Bayesian methodology using conjugate priors for a dynamic application of our algorithm to be run with data collection regardless of system complexity. Also included are model selection methods using machine learning techniques for the optimisation of individual noise wave parameters to combat overfitting and underfitting, the results of which converge with that of a least-squares approach when wide priors are adopted. Our pipeline easily incorporates many more calibrators than the standard four shown in Fig. 2 to increase constraints on noise wave parameters while identifying possible correlations between them. A schematic of the improved calibration method is shown in Fig. 3.

In order to simplify our calibration approach, we first define the following terms

which represent initial calibration measurements on $D$ in the frequency domain for the characterisation of $N_{\mathrm{data}}$ from equation 3 via our noise wave parameters. It is expected that calibration-related deviations of $D$ in the time domain are sufficiently curtailed through practical strategies such as temperature control of the receiver environment. Incorporating these into equation 9, with some rearrangement, then gives the equation

at each frequency. Here, there are no squared or higher-order terms, allowing us to take advantage of the linear form by grouping the data and noise wave parameters into separate matrices

In these equations, all of our data; the reflection coefficient measurements and power spectral densities, are grouped in an X vector which forms a matrix where one of the axes is frequency. The calibration parameters as frequency-dependent polynomials of varying degree are collected into a $\boldsymbol{\boldsymbol{\Theta}}$ vector which serves as our model describing $N_{\mathrm{data}}$. Applying these definitions condenses the calibration equation into

where $\textbf{{T}}_{\mathrm{cal}}$ is a vector over frequency and $\sigma$ is a noise vector representing our error. Since EDGES assumes that each power spectral density measurement is frequency independent, we have assumed that $\sigma$ is a multivariate normal distribution. This assumption is implicit in the EDGES analysis in which they use a least-squares minimisation approach for solving model parameters.

For calibration of the receiver, we are concerned with the construction of predictive models of the noise wave parameters, $\boldsymbol{\Theta}$, in the context of some dataset, T. We can use $\boldsymbol{\Theta}$ to calculate the probability of observing the data given a specific set of noise wave parameters:

where, $N$ is the number of measurements. This distribution on the data is the likelihood. For the purposes of calibration, T may be $\textbf{{T}}_{\mathrm{cal}}$ measurements or alternatively, $\textbf{{T}}_{\mathrm{sky}}$ for prediction of a sky signal. Our model must also specify a prior distribution, quantifying our initial assumptions on the values and spread of our noise wave parameters which we specify as a multivariate normal inverse gamma distribution:

which is proportional up to an integration constant. Here, $a$ and $b$, which are greater than zero, along with $\textbf{{V}}_{\boldsymbol{\Theta}}$ and $\boldsymbol{\mu}_{\boldsymbol{\Theta}}$ represent our prior knowledge on the noise wave parameters. $d$ is the length of our vector $\boldsymbol{\Theta}$.

Equation 18 is determined by a set of values for our model $\boldsymbol{\Theta}$. We can marginalise out the dependence on $\boldsymbol{\Theta}$ and our noise term by integrating over the prior distribution by both $\boldsymbol{\Theta}$ and $\sigma^{2}$ at once. Following the steps in Banerjee (2009)

where

and $\Gamma\left(x\right)$ represents the Gamma function, not to be confused with the notation for our reflection coefficients. Equation 20 is the evidence, which gives the probability of observing the data $\textbf{{T}}_{\mathrm{cal}}$ given our model.¹¹1It is in fact better to use the equivalent more numerically stable expression $b^{*}=b+\boldsymbol{q}^{\top}\boldsymbol{q}+\boldsymbol{q}^{\top}\textbf{{X}}\textbf{{V}}_{\boldsymbol{\Theta}}\textbf{{X}}^{\top}\boldsymbol{q}$, where $\boldsymbol{q}=\textbf{{T}}_{\mathrm{cal}}-\textbf{{X}}\boldsymbol{\mu}^{*}$ to avoid cancellation of large terms.

With the prior distribution specified, we use Bayes’ equation to invert the conditioning of the likelihood and find the posterior using the likelihood, prior and evidence:

Similarly from Banerjee (2009), this can be written as

The posterior distribution represents the uncertainty of our parameters after analysis, reflecting the increase in information (Nagel, 2017). We highlight the difference between the ‘likelihood-only’ least-squares approach versus the Bayesian approach with the former being a special case of the latter with very wide priors demonstrable when $\textbf{{V}}_{\boldsymbol{\Theta}}\rightarrow\infty\Rightarrow\textbf{{V}}_{\boldsymbol{\Theta}}^{-1}\rightarrow 0$, and $\boldsymbol{\mu}^{*}$ becomes $\boldsymbol{\Theta}$. The transition from ‘non-starred’ variables to ‘starred’ variables represents our ‘Bayesian update’ of the prior to the posterior noise wave parameters in light of the calibration data $\textbf{{T}}_{\mathrm{cal}}$.

As we can see, the posterior distribution is in the same probability distribution family as equation 19, making our prior a conjugate prior on the likelihood distribution. The use of conjugate priors gives a closed-form solution for the posterior distribution through updates of the prior hyperparameters via the likelihood function (Banerjee, 2009; Orloff & Bloom, 2013). The resulting numerical computation is many orders of magnitude faster than MCMC methods relying on full numerical sampling and permits an in-place calculation in the same environment as the data acquisition. This becomes particularly useful for the speed of the algorithm as frequency dependence is introduced in which the computations would not be manageable without conjugate gradients.

To allow for a smooth frequency dependency, we promote each of our noise wave parameters in section 2.3 to a vector of polynomial coefficients

where $i$ is our noise wave parameter label; $i\in\{\mathrm{unc,\ cos,\ sin,\ NS,\ L}\}$, modelled using $n+1$ polynomial coefficients. Likewise

where $\nu$ is a vector of input frequencies which are raised to powers up to $n$. For a vector of $n$’s attributed to our calibration parameters, under this notation multiplication in equation 17 is element-wise and equation 20 is effectively $p\left(\textbf{{T}}_{\mathrm{cal}}|\textbf{{n}}\right)$. Assuming a uniform prior on n, inverting Bayes’ theorem gives $p\left(\textbf{{n}}|\textbf{{T}}_{\mathrm{cal}}\right)$ for use in model comparison in which the relative probabilities of models can be evaluated in light of the data and priors. Occam’s razor advises whether the extra complexity of a model is needed to describe the data (Trotta, 2008), permitting optimisation of the polynomial orders for individual noise wave parameters as detailed in section 3.3. By taking a random sampling of the resulting posterior, we characterise the noise wave parameters as multivariate distributions depicted in contour plots which exhibit a peak value accompanied by $1\sigma$ and $2\sigma$ variance as well as correlation between parameters inferred from a covariance matrix.

Following characterisation of the receiver, we next apply the $\textbf{{T}}_{\mathrm{cal}}$ from our calibration to a set of raw antenna data $\hat{\textbf{{X}}}$ for prediction of our sky signal, $\textbf{{T}}_{\mathrm{sky}}$, from equation 3. The predictions for the data follow from the posterior predictive distribution

The first probability in the integral is the likelihood for our antenna measurement $\textbf{{T}}_{\mathrm{sky}}$ and the second is our posterior from equation 23. Following the steps in Banerjee (2009), this can be shown to be a multivariate Student’s t-distribution written as:

where $I$ is the $N\times N$ identity matrix and $a^{*}$, $b^{*}$, $\boldsymbol{\mu}^{*}$ and $\textbf{{V}}^{*}$ are defined in equation 21. This new distribution on $\textbf{{T}}_{\mathrm{sky}}$ corresponds to a set of points with error bars and represents the calibrated sky temperature as the output of the receiver.

In this section, we show the first major feature of our Bayesian pipeline; the consideration of correlation between noise wave parameters when deriving their values. This is best demonstrated when noise is introduced in an idealised way as to retain a form matching the Gaussian form of our mathematical model. To do this, empirical models of power spectral densities are calculated from equations 6, 7 and 8 using measurements of $\Gamma_{\mathrm{rec}}$, $\Gamma_{\mathrm{cal}}$ and $T_{\mathrm{cal}}$ for the cold and hot loads, as well as a set of realistic noise wave parameters. Gaussian noise of one unit variation is then added to the $T_{\mathrm{cal}}$ measurements after the calculation to conserve its Gaussian form. This data is submitted to our algorithm and the resulting posterior distributions for coefficients of the polynomial noise wave parameters are compared to the initial values.

Such posterior distributions can be seen in Fig. 5 showing the results of models using only the cold load (grey posterior), only the hot load (red posterior) and using both loads in tandem (blue posterior). For these calculations we chose a set of model noise wave parameters as constants across the frequency band;

In Fig. 5, a strong correlation between the $T_{\mathrm{L}}$ and $T_{\mathrm{NS}}$ is evident as the hot-load posterior is highly skewed as expected from equations 11 and 14. The resulting intersection of posteriors from the individual loads facilitate the derivation of noise wave parameters as the dual-load posterior is found within the region of posterior overlap crossing with the values of the model shown in the inset of Fig. 5. Retrieval of the noise wave parameter values using correlations between them found in the data demonstrate the relevance of this information which is not taken into account in previous calibration techniques.

A nice feature of our pipeline is the ability to include as many calibrators as required to constrain the calibration parameters. For analysis, six more calibrators are introduced in pairs following the order presented in table 1. We include data generated from measurements of multiple resistors terminating a high quality 25 m cable made by GORE^®. Data for these calibrators is once again generated using fixed terms and Gaussian noise of one unit variation added to $T_{\mathrm{cal}}$ as discussed above. Fig. 6 shows the results of models using four, six, and eight calibrators.

As shown, the inclusion of more calibrators increases the constraint on the resulting noise wave parameters. However, we note that after the inclusion of four calibrators, the relative additional constraint decreases with each additional calibrator and thus the use of more than eight calibrators would be unnecessary. The values of noise wave parameters used to generate the data as indicated by the cross hairs in Fig. 6 all fall within $1\sigma$ of our pipeline’s resulting posterior averages for models using all eight calibrators.

The final highlight of our Bayesian pipeline is a the use of machine learning techniques to optimise individual noise wave parameters. This is advantageous as a blanket set of order-seven polynomials applied to all noise wave parameters, such as done in the EDGES experiment, may underfit or overfit individual parameters and misidentify systematics or information about the signal being measured.

The optimisation procedure compares the evidences (equation 20) of different models to determine the vector of noise wave parameter polynomial coefficients n that best describes the data as briefly mentioned at the end of section 2.3. Since the model favoured by the data will have the highest evidence, we use a steepest descent procedure to compare models in ‘n-space’ and determine the direction of the gradient in ‘evidence-space’. After multiple iterations, this brings us to the model with the maximal evidence. Since n consists of five numbers corresponding to the number of polynomial coefficients for each of the five noise wave parameters, models are generated by individually increasing each index of n by 1. We expect the evidence to follow an ‘Occam’s cliff,’ in which the evidence sharply increases preceding the optimal n with a slow fall off following the maximum.

To demonstrate this, data is generated using measurements from all eight calibrators of table 1 and noise wave parameters as second-order polynomials

where $x$ is our normalised frequency. Gaussian noise of one unit variation is applied to the calibrator input temperatures as before. The evidences of various models are plotted in Fig. 7 in which an Occam’s cliff can be seen peaking at polynomial order two. As expected from the plot, the steepest descent algorithm finds that noise wave parameters modelled as second-order polynomials best describe the data.

To demonstrate the robustness of our pipeline, we conducted self-consistency checks using empirically modelled data with a more complicated noise model. This data was generated using reflection coefficients of eight calibrators and the receiver measured in the laboratory. These reflection coefficients were then smoothed using a cubic smoothing spline (Weinert, 2009) in order to maintain their approximate shape over frequency. The same second-order noise wave parameters detailed in section 3.3 are used with the reflection coefficients to generate our model power spectral densities. Following this, we added of order 1% Gaussian noise independently to the smoothed $\Gamma_{\mathrm{rec}}$ and $\Gamma_{\mathrm{cal}}$ as well as $P_{\mathrm{cal}}$ to more accurately represent the instrument noise from measurement equipment such as vector network analysers. No noise was added to the calibrator input temperatures. This results in a model that does not match the Gaussian form of our mathematical model as in the previous sections and thus does not demonstrate the features of our pipeline as explicitly, but is more representative of data set expected from measurements in the field. Data for the receiver and the cold load generated using this noise model are shown in Fig. 8.

Using data generated for all eight calibrators with our realistic noise model, the calibration algorithm selects optimal polynomial orders matching those of the model noise wave parameters whose values fall within within $1\sigma$ of the posterior peak values as shown in Fig. 9. For these higher order tests, we use fgivenx plots which condense noise wave parameter posteriors into samples that can be compared to the model parameter values instead of comparing each individual coefficient (Handley, 2018).

When this calibration model is used to calibrate an ambient-temperature $50\ \Omega$ load, the RMS error between the calibrated temperature and the measured temperature is 8 mK, well within the $1\sigma$ noise level (bottom right panel of Fig. 9). This level of accuracy is comparable to the 26 mK noise floor estimated for the EDGES pipeline in 2016 (Monsalve et al., 2017).

By individually adjusting each component of noise arising in our realistic noise model, we may determine what kind of noise our calibration algorithm is most sensitive to, as well as calculate the maximum amount of noise permissible for a specified level of systematic feature reduction. These topics are intended to be explored in a future work.