Singular spectrum analysis of time series data from low frequency radiometers, with an application to SITARA data

In this section, we describe the notation employed and define certain matrices that are useful for the subsequent analysis. Many of these definitions can be found in Davis (1979) as well as Olson et al. (2014). We use bold capital Roman letters such as $\bm{X}$ to denote matrices, and small Roman letters such as $n$ are used for indexing. We use zero-based numbering such that the indices start at $0$. Small Roman letters with an arrow such as $\vec{u}$ denote vectors. Integers are denoted by capital Roman letters such as $L$. We now proceed to define some basic matrices and associated linear algebra.

Consider a time series $x[n],n=0,1,....N-1$ of length $N$ that is strictly periodic with a period given by $L$.

We also assume that the length $N$ of the time series obeys the relation $N=zL-1$, where $z$ is an integer. Forming a trajectory matrix out of this sequence gives the following matrix.

As $K$ is an integer multiple of $L$, the trajectory matrix can be interpreted as a block matrix, and partitioned into $M=\frac{K}{L}$ square submatrices as given in Eq.10, where $\bm{X_{ac}}$ is an anti-circulant matrix of dimensions $L\times L$.

Before proceeding further, some observations can be made.

1. 1.

   We find that $rank(\bm{X})\leq L$, as there can only be a maximum of $L$ linearly independent columns in $\bm{X}$. Consequently there can be a maximum of $L$ non-zero singular values for $\bm{X}$.

2. 2.

   The anti-circulant matrix form given in Eq.10 is real symmetric. Therefore, it has an eigendecomposition given as $\bm{X_{ac}}=\bm{Q\Lambda Q^{T}}$, where $\bm{Q}$ is an $L\times L$ orthonormal matrix of eigenvectors and $\bm{\Lambda}$ is an $L\times L$ matrix of the corresponding eigenvalues.

This enables us to write the trajectory matrix as a block matrix $\bm{X}=\begin{bmatrix}\bm{Q\Lambda Q^{T}}~{}~{}\bm{Q\Lambda Q^{T}}..\end{bmatrix}$. Without going through the pedagogical details, we state that $\bm{X}$ can be decomposed into three matrices as given below.

where, in the last step we recall the SVD of $\bm{X}$ to facilitate a direct comparison.

We now attempt to find the nature of this decomposition - specifically what one can expect from the singular vectors if the time series is periodic in nature. From Eq.12, it can be noted that the singular values are a scaled version of the eigenvalues of anti-circulant matrix $\bm{X_{ac}}$. Inspecting the block matrix representation in Eq.12, it can also be seen that the matrix that forms to the left and right singular vectors is the matrix of eigenvectors $\bm{Q}$. Therefore, all that we need to know is the nature of $\bm{Q}$. A calculation of this is provided in Appendix.A, where it is shown that $\bm{Q}$ consists of sinusoidal eigenvectors if the data are periodic. Therefore, the left singular vectors - which are exactly the eigenvectors contained in the matrix $\bm{Q}$ - are all sinusoidal. Besides, the vectors in $\bm{Q}$ are all periodic in $L$ and thus when arranged as blocks in the right singular vector matrix, they form continuous sinusoids. Thus, the equivalence of decompositions in Eq.12 shows that when SSA is applied to a strictly periodic sequence, the singular vectors obtained are sinusoidal.

In Appendix.A, it is also shown that the eigendecomposition of anti-circulant matrices such as $\bm{X_{ac}}$ has a deep relation to discrete Fourier transforms (DFT). Indeed, when applied to an ideal periodic series, SSA results are directly comparable to DFT of one cycle of periodic data; this is also demonstrated with simulations in Sec.4.1.1. However, if the data have aperiodic structures, DFT of the entire span of data is not a suitable choice as aperiodic structures are not well represented in a sine/cosine Fourier basis. In such a scenario, SSA is better suited to reveal aperiodic structures in the data. Therefore, regardless of the nature of data, SSA can be applied to them and if the data are periodic, SSA reduces to a form of Fourier analysis.

Even though the decomposition given may be treated as an SVD of the trajectory matrix, some caution has to be exercised. Since the singular vector matrices have to be unitary, the right singular matrix gets divided by a scaling factor of $\sqrt{\frac{K}{L}}$ while the singular values get multiplied the same factor. Also, the left singular matrix in SVD is unitary, while the eigenvector form given in Appendix.A has unit amplitude. The corresponding scaling applies to the singular value, however it is inconsequential when the eigendecomposition and SVD are computed with numerical packages. Moreover, eigenvalues can be positive and negative, and as given in Appendix.A, they occur in positive and negative pairs in this context. However, singular values are always non-negative and therefore the signs of the eigenvalues get applied to the singular vectors when equated with the SVD given in Eq.12. Also, we assume that both the singular values and eigenvalues (and corresponding vectors) have been ordered in the same fashion - typically in descending magnitude.

The calculations and simulations so far assume a case where the radiometer system has been fully calibrated to yield data calibrated to a reference plane outside the Earth’s ionosphere. Observational radiometric data have multiple nonidealties, the dominant ones are additive components such as receiver noise temperatures and mutliplicative factors that are often called gains. Both of them can be time-varying and therefore are required to be known to correct the data. We now incorporate them into our model.

Assuming a simple system model for SITARA data (see Sec. 5.2 of Thekkeppattu et al. (2022)), the trajectory matrix of the measured data can be expressed as a Hadamard product, given in Eq.14.

where $\bm{X_{s}}$ is the sky temperature trajectory matrix, $\bm{R}$ is the receiver noise temperature trajectory matrix and $\bm{G}$ is the trajectory matrix of time-varying system gains. Since the low frequency radio sky is bright, the dominant time-varying component in the data is due to sky temperature multiplied with time-varying system gains. Therefore, we can apply a simplifying assumption that the receiver noise temperature is constant with time and absorb it into the “ideal” sky matrix, and the sum can be written as $\bm{X}$. Thus we write

While Eq.14 and 15 may be written without resorting to a matrix formulation, expressing them with trajectory matrices enables SSA.

Based on the summation of $\bm{X}$ from Eq.5 and Eq.15, the following relations can be derived.

where we have used the distributive property of Hadamard products. Thus we have,

Though Eq.17 appears trivial, it reveals a useful aspect of the decomposition. Applying diagonal averaging to each $\bm{X^{\prime}_{i}}$ leads to gain multiplied reconstructed series, as the Hadamard product is an element-by-element multiplication. However, we already know the “true” reconstructed series to be purely sinusoidal from Sec.4. The multiplication of a pure sinusoid waveform by a signal of much lower frequency results in amplitude modulation of the sinusoid, a process that is called amplitude modulation (AM) in communication systems (Haykin & Moher, 2008). Therefore, we may treat the gain multiplied reconstructed series as a form of AM signals with the sinusoids acting analogous to “carriers”; and using AM demodulation techniques to extract the temporal variations in gains.

Before proceeding further, it is instructive to verify the above calculations with simulations. For this, we begin with the basic time series from Sec.4. To this series, we add a constant of $100~{}K$ as the receiver noise temperature and simulated radiometric noise (integration time of 15 minutes and bandwidth of 61 kHz) for the overall system temperature. The gain profile is modelled as a random-walk, which is generated as a cumulative sum of Gaussian white noise. The resulting profile is subsequently smoothed with box-car averaging to suppress variations of time cadence less than one sidereal day. The gains are also normalised to avoid negative excursions. The time series is then multiplied with the simulated gains and the resulting series and the associated trajectory matrix are shown in Fig.6. It may be noted that the dimensions of the matrices and time series have all been kept the same.

The trajectory matrix is then decomposed with SVD, and the eigentriples grouped in the same manner as in Sec.4 and diagonally averaged. It is easy to notice that the resulting reconstructed series shown in Fig.7 are indeed amplitude modulated sinusoids.

Several algorithms exist to demodulate such AM signals to obtain their envelopes. Here we use an algorithm that is mathematically simple to interpret. The first step in this algorithm is converting a real valued signal into an analytic signal, which can be accomplished with Hilbert transforms. Taking the magnitude of this analytic signal yields the modulation envelope. We apply this procedure to the series given in Fig.7 . The resulting gain profiles are plotted in Fig.8, in which the curves are also normalised according to the method described in Sec.4.3.

Fig.8 validates the interpretation of Eq.17 that the reconstructed series have identical gain templates that can be recovered with amplitude demodulation. Indeed, a major outcome of this exercise is that the gain template for the DC component is identical to those of the periodic components. Therefore the time-variations in the DC term may be corrected for by using the gain templates from periodic components. This leads us to a novel calibration strategy for global 21–cm radiometers that can correct long-term drift and related systematics.

The recovered gains require normalisation to enable a comparison with the injected gains and make them useful for calibration. We begin with Eq.17 and incorporate Eq.5 as shown in Eq.18, where we also include the eigentriple grouping. The indices have also been modified for convenience.

where we used the property that singular values are identical for the eigentriples in a grouped pair. The vectors $\vec{u}$ and $\vec{v}$ are from unitary matrices and therefore their individual inner products equal to unity. However, for subsequent analysis, these vectors need to be normalised to have unity amplitude when diagonally averaged after grouping. The normalisation can be achieved with the following operations, where the sinusoidal terms in the brackets have been normalised to have unity amplitude.

as the lengths of the vectors $\vec{u}$ and $\vec{v}$ are $L$ and $K$ respectively. The grouped and diagonally averaged reconstructed series is given as

where $g[n]$ is the gain series and $c_{i}[n]$ is the “carrier” sinusoid. The normalisation in Eq.19 ensures unity amplitude for $c_{i}[n]$. Upon demodulation of a series $\widetilde{x_{i}}[n]$, the carrier is removed and a modulation envelope $g_{i,t}[n]$ is obtained as

From these envelopes, gain templates $g^{\prime}_{i}[n]$ can be recovered as

For Eq.22 to be useful for calibration, the only auxiliary information required is the singular value ${\sigma_{i}}$. This is the singular value of the component with the same periodicity of that we expect from the true sky, and can be obtained from the SSA of the simulated ideal sky time series from Sec.4. Using those singular values, we obtain the recovered and appropriately normalised gains that are plotted in Fig.8.

We are now in a position to apply those gains to the $0^{th}$ reconstructed series to achieve calibration and compare the calibration with the expected levels. The results are given in Fig.9.

We find that the mean-sky component has been calibrated solely using the periodic component(s) of the sky with the aid of SSA. It has to be noted that the prescription provided here to calibrate the mean-sky relies only on the singular values of the simulated periodic components. Even without using the singular values, and only using the data, it is possible to achieve relative gain calibration, while applying gains with appropriate singular values to the mean-sky component establishes its proper brightness temperature scale.

Let us now try to obtain an intuitive appreciation of this calibration. The trend component obtained from SSA is essentially a box-car averaging of time series data. If the window used for averaging encompasses exactly one sidereal day, the periodic variations caused by the Galaxy transiting through the beam are averaged out. If the system were perfect the result would be a constant temperature, which is the sum of mean-sky and receiver noise temperatures. Thus the trend that we observe in the non-ideal system may be modelled by a constant value multiplied with the time-varying gains.

The gain variations that we have studied so far have smooth evolution over time with no periodicity. However, for real radio telescope systems that do not have temperature regulation, gain variations with local temperature are expected - details given in Sec.5. Indeed, the rising and setting of the Sun can heat and cool the components in the analog signal chain and induce diurnal variations in the gain that may be correlated with the sky-drift. Yet another potential origin for diurnal patterns in gains is the ionosphere. Therefore it is imperative to study such a scenario, which we perform next.

To the gains simulated in the previous section, we add a sinusoidal diurnal gain component. For a time series that spans a month, the diurnal gain variations are expected to show a strong correlation with Galaxy transit. For longer time series, the difference between sidereal time and civil time reduces this correlation. The modified gains and the resulting time series are shown in Fig.10.

We perform SSA in the exact same manner as in the previous case, and recover normalised gains, as shown in Fig.11.

It is interesting to note that the correlation between gains and sky patterns reduces the separability of the two, and the recovered gains from the reconstructed series differ from one another. We attempt a calibration of the mean component with these gains and the results are given in Fig.12.

The differences in the recovered gains in Fig.11 as well as the calibrated temperatures in Fig.12 demonstrate a fundamental limitation of the SSA technique in calibrating data when the gain variations are correlated with sky patterns. Lacking auxiliary information related to the cause of such gain variations, for example temperature measurements, \colorblack the calibration technique fails and yields erroneous results as demonstrated. On the other hand, a disagreement between the recovered gains can be used as a practical indicator of gain variations that have frequencies coinciding with the periodicity induced by sky drift, thus pointing to a necessity for subsidiary information. In other words, differences in recovered gains between the groupings indicate presence of correlations between the gains and sky drift. It may be also noted that no a priori information is required to infer such a correlation. Having demonstrated the various cases of SSA, we now apply these techniques to data from SITARA observations.

We use a time series of SITARA spectral data from the concatenated 30 days of June 2021 data. The individual data files in miriad format are converted to hdf5 format and concatenated using custom tools written in python. From the concatenated dataset, auto-correlation data for a frequency of 111 MHz are extracted, as it is a frequency with relatively low RFI occupancy. The major source of RFI observed during this period was Solar radio bursts that are transient in nature. This raw data time series as a function of Julian Date (JD) is shown in Fig.1.

The raw SITARA data spanning a month have about 749,000 points, which makes calculations difficult. Therefore, the data are binned into the same number of bins (2879) as done in the simulations, giving one sample per 15 minutes or 96 samples per sidereal day. Once again, in performing binning it is important to use sidereal time - and not Julian time - to ensure an equal number of samples for each cycle of sky-drift. The data are then embedded into a trajectory matrix and SSA is applied in the same manner as before. The eigentriples are grouped and reconstructed series formed.

Fig.13 shows first 4 of the reconstructed series. The trend plot also has a low time cadence record of the temperature recorded within the Murchison shire, which is available from NCEI-NOAA ²²2https://www.ncei.noaa.gov/access/search/index. The temperature scale has been inverted, as the trend shows an anti-correlation with temperature. Further, the radio data and temperature data fall into two different families with different sampling, and therefore the temperature data have been interpolated to the same Julian day bins as the radio data. The Pearson correlation coefficient computed between the trend and the temperature data is -0.83, which denotes a strong anti-correlation between the trend and temperature.

Apart from component ageing effects, there are two major reasons for such temperature-induced variations in radiometric data.

1. 1.

   The noise temperatures of active devices used in amplifiers increase with physical temperature³³3It may be noted that this is one of the reasons for cryogenic cooling of radio telescope front-ends to achieve a low overall system temperature.. The same holds true for passive components such as attenuators. If this is the cause of drift, the pattern so obtained is expected to be correlated with physical temperature.

2. 2.

   The gain of amplifiers reduces as temperature is increased. In this case, the trend pattern and physical temperature would be anti-correlated.

It is therefore evident from Fig.13 that the major contributor to the trend is temperature-induced gain variations. The receiver noise would inevitably vary as a function of physical temperature, however when the overall system temperature is sky-dominated the impact of this would be secondary to gain variations. This informs our choice of calibration model given in Sec.4.3, where we assume a constant receiver noise temperature and a time-varying gain.

Subsequently, the gain patterns are recovered from the reconstructed series with the demodulation technique. These recovered gains for the first two periodic components are shown in Fig.14. As can be seen, the normalised gains are different between the series, thus pointing to potential diurnal gain variations. Therefore, we refrain from applying these gains to the mean component to establish its brightness temperature scale.

Since this is the first work exploring SSA for radiometric data analysis and calibration, we have kept the analysis simple. Specifically, low frequency radio data can have RFI and we ignored the flagging that has been performed. Flagging considerations lead to gaps in a time series dataset and may necessitate data in-painting or interpolation. Such explorations will be taken up in future work.

We have developed the mathematical background and demonstrated an application of SSA for auto-correlations at a single frequency, while SITARA has nearly two thousand usable frequency channels in auto-correlations and cross-correlations each. To use all the frequency channels and to consider complex cross-correlations, the mathematical framework developed in this paper has to be expanded. For a proper treatment of multiple frequency channels, multivariate SSA techniques will be explored in future work.

While we demonstrated via simulations a calibration technique to remove gain variations that evolve smoothly with time, the application of it to SITARA data is hampered by diurnal gain variations. This shows the necessity to maintain temperature regulation or inclusion of thermometers at points in the signal chain to track such gain variations - this will be incorporated in future revisions of the SITARA system. Additionally, the data in frequency channels below 50 MHz, where receiver temperature dominates over sky temperature, show diurnal periodicity arising from such gain variations which may be exploited as a template to improve calibration.

As we demonstrated, a major limitation with the SSA technique for calibration is the confusion between diurnal gain variations and sky drift, as both have approximately the same periodicity. If the periodicities can be made to differ, the calibration can be improved significantly. Specifically, if the sky-drift can be made faster than the gain variations, substantially better calibration can be expected. While it is difficult for ground-based experiments to introduce such a separation, a radiometer payload on a spin-stabilized satellite can have beams that rapidly scan the sky, thereby increasing the sky-drift rate. For example, a broadband dipole antenna could be placed on a satellite spinning such that the nulls sweep the Galactic plane at a rate much faster than any gain variations, including the flicker noise of the electronic systems. Spinning a satellite at a rate faster than than the $1/f$ knee frequency of the radiometric system, to reduce the deleterious impact of gain fluctuations on the images, has been employed in CMB missions such as Planck (Bersanelli et al., 2010). As there are a few projects proposed, planned or launched targeting the global 21–cm signal from a satellite platform such as DARE (Burns et al., 2012), Longjiang/Chang’e-4 (Jia et al., 2018), DSL (Chen et al., 2021), PRATUSH etc., we opine that SSA would be an ideal tool for analysis and calibration of time series data from such radiometers in space.