Visibility Analyser

In the SKA, visibilities (i.e. correlated data) are produced by correlators. This stream of data is carried over using SPEAD software protocol over UDP. In order to analyse visibilities, the final and integrated system would use software from SDP. In the meantime the visibility analyser offer an alternative for visibilities analysis more tailored for the analysis of Low CBF output and thus focusing more on the immediate needs of Low CBF.

Initialisation

The visibility analyser class only takes a single json configuration string, it should contain the following for 2 stations and 2 channels:

{
    "coarse_channels": [124, 125],
    "stations": [345, 350]
}

Using this json string, let’s now create the object

visibilities = VisibilityAnalyser(json.dumps({ "coarse_channels": [channel for channel in range(330, 330+96)],
                                                "stations": [1, 2, 3, 4, 5, 6]  }))

Once the object initialised, the visibility analyser can now extract visibilities from the pcap file. This is done as follows:

visibilities.extract_spead_data("/path/to/pcap/visibilities.pcap")

At this point, the VisibilityAnalyser object would contain the following variables accessible publicly:

coarse_channels: list of coarse channels to analyse
stations: list of stations to analyse
nb_station: the number of stations
visibilities: a data structure containing all visibility baselines (complex numbers)
time_sequence: the set containing all timestamps from the data

Use cases

As mentioned above, Visibility Analyser is tailored for various test cases in Low CBF. As such, we will go through those use-cases one by one and will link them to the Perentie Test Cases when possible.

Correlation Matrix

In order to check that the correlator, we are going to populate the correlation matrix. In the analyser, this is done by constructing two data structures:

corr_abs_matrixes_time, containing correlation matrix of sum of norm of visibilities in various coarse channels

corr_abs_matrixes_time = {
        time1: {
            channel1: [2 * self.nb_station, 2 * self.nb_station] correlation matrix of sum of norm of visibilities in channel1
            channel2: [2 * self.nb_station, 2 * self.nb_station] correlation matrix of sum of norm of visibilities in channel2
        }
        time2: {
            channel1: [2 * self.nb_station, 2 * self.nb_station] correlation matrix of sum of norm of visibilities in channel1
            channel2: [2 * self.nb_station, 2 * self.nb_station] correlation matrix of sum of norm of visibilities in channel2
        }
    }

corr_abs_matrixes_time_fine_channels, containing correlation matrix of norm of visibilities in various fine channels

corr_abs_matrixes_time_fine_channels = {
        time1: {
            channel1: [2 * self.nb_station, 2 * self.nb_station] correlation matrix of norm of visibilities in channel1
            channel2: [2 * self.nb_station, 2 * self.nb_station] correlation matrix of norm of visibilities in channel2
        }
        time2: {
            channel1: [2 * self.nb_station, 2 * self.nb_station] correlation matrix of norm of visibilities in channel1
            channel2: [2 * self.nb_station, 2 * self.nb_station] correlation matrix of norm of visibilities in channel2
        }
    }

Those 2 data structure are populated using the following method:

visibilities.populate_correlation_matrix()

Once the two correlation matrix data structures are completed, we can generate heatmap figures using the following method:

visibilities.save_heatmap_to_disk(mode, frequency_to_check=None, case="")

Where

mode represents 4 different modes so far:
- average_all: heat map of the sum of correlation across time for all coarse channels
- average_frequencies: produces 1 figure per integration time across all frequencies
- average_times: produces 1 figure per channel across all integration times
- average_times_fine: produces 1 figure per visibility channel for frequency_to_check coarse channel
frequency_to_check is the coarse channel to further analyse when using the average_times_fine mode
case is the name of case study for the various figures

Angle Analysis

Next we analyse the phase of correlated signal between pairs of stations. To do so the VisibilityAnalyser offers the following method

visibilities.angle_analysis(station1: int, station2: int, polarisation: str, nb_time_period: int, coarse_channels, case="")

where:

station1: the first reference station to generate angle against
station2: the second reference station to generate angle against
polarisation: the polarisation of interest, 4 are possible XX, XY, YX, YY
nb_time_period: the number of integration we want to analyse
coarse_channels: the list of coarse channels to plot
case: case study for the various figures

As a result of this function, the VisibilityAnalyser would produce a figure plotting the frequency vs phase plot where the fine frequencies are on the x-axis and the phase of the correlated signal is on the y-axis.

PTC 14: Visibility Biases

A typical example of the VisibilityAnalyser purpose is the first Perentie Test Case (PTC) that was tackled during PI#19, namely PTC#14 Visibilities Biaises. This test aims at demonstrating Low CBF compliance to our Requirement 199 that states that when processing white noise uncorrelated signals, the average of any visibility shall be proportional to the integration time raised to power of -0.5, with and accuracy of 10% for integration up to 1000 hours.

While this test requires multiple steps for the data production as explained in more details in documentation, the analysis part relies on the use of 2 methods to calculate the average power in the uncorrelated signal from the correlator (i.e. only cross-correlation between stations)

The first, and maybe less important, method is

visibilities.rms_total()

This function calculate the entire dataset calculate its RMS - square, sum, average and take square root. As a result two RMS numbers come out, one for the auto-correlations and one for the cross-correlations. Along the way this function will also print some useful information such as the ratio of cross vs auto RMS calculation.

The second, and maybe main, method for PTC#14 is

visibilities.rms_slice_total(slice_size: int, index_rms: int, case="")

Where

slice_size: the number of cross correlation baseline to group together. Rule of thumb slice = nb_baseline/sum_1 vector (power of 2)
index_rms: if not using the full default 64k (2^16) vector, needs to indicate the maximum power we are looking into
case: case study for the various figures

In a nutshell, this function, for the entire dataset, calculates two giant vectors sum_1 and rms1. From these two vectors we can generate various figures including the main figure of PTC 14 Those figure are the following

{case}ptc_14_rms.png: plot of the RMS Standard Deviation vs number of emulated hours
{case}ptc_14_sum.png: plot of the Sum of baselines Standard Deviation vs number of emulated hours
{case}ptc_sum_real.png: plot of mean values of real component of the cross correlations, with its standard deviation vs number of emulated hours
{case}ptc_sum_dev_real.png: plot of the standard deviation of real component of the cross correlations vs number of emulated hours
{case}ptc_sum_imag.png: plot of mean values of imaginary component of the cross correlations, with its standard deviation vs number of emulated hours
{case}ptc_sum_dev_imag.png: plot of the standard deviation of imaginary component of the cross correlations vs number of emulated hours
{case}rms_ptc_14.png: plot of the RMS for a given integration period vs number of emulated hours

The main figure being the last one.

Phase And Amplitude closure

Calculate the phase and amplitude closure following documentation from the doc

We then need a list of at least 4 stations to perform as follows:

Closure (triple) phase: arg( V12*V23*conjg(V13) )
Closure (quad) amplitude: abs( (V12*V34)/(V14*conjg(V34)) )
Triple amplitude: abs( V12*V23*conjg(V13) )**0.3333
Quad phase: arg( (V12*V34)/(V14*conjg(V34)) )

The closure phase, quad phase and closure amplitude should be independent of antenna-based errors, and for a point source should have values of zero phase or unit amplitude. The triple amplitude is independent of antenna-based phase errors (but not amplitude errors), and for a point source is a measure of the flux density.

We plan to implement this functionality in the following method

visibilities.phase_and_amplitude_closure( stations_closure, coarse_channels)

Where * stations is the the list of station to goes through in circle * coarse_channels are the channels of interest

At the moment, we are only reporting closure phase and amplitude. Those two are returned from the function once called.

For example, when we configured an experiment with a CNIC generating traffic equivalent to 6 stations with stations 1, 2, 3, and 4 seeing the same noise source with no delay, we were able to produce the following result by calling:

import numpy as np
import matplotlib.pyplot as plt

(closure_angle, closure_amp) = visibilities.phase_and_amplitude_closure([0,1,2,3],range(171, 171+24))
plt.clf()
plt.plot(closure_angle, ".-")
plt.ylabel("phase closure")
plt.xlabel("visibitility channel")
plt.show()

plt.clf()
plt.plot(closure_amp, ".-")
plt.ylabel("amplitude closure")
plt.xlabel("visibitility channel")
plt.show()

Autocorrelation analysis

Generate the 3D auto amplitude-freq-time analysis figure for a station. This is done through the call of the following function

visibilities.auto_station_analysis(station_nb,
    nb_time_period,
    polarisation,
    coarse_channels)

Where:

station_nb is the station number.
polarisation is the polarisation of interest, 4 are possible XX, XY, YX, YY
nb_time_period is the number of integration we want to analyse
coarse_channels is the list of coarse channels to plot

For example, when we configured an experiment with a CNIC generating traffic equivalent to 6 stations with stations 1, 2, 3, and 4 seeing the same noise source with no delay, we were able to produce the following result by calling:

visibilities.auto_station_analysis(1, 20, "XX", [channel for channel in range(171, 171+24)])

The call above aims at generated the 3D mesh plot for station 1, XX polarisation, for 24 coarse channels equivalent to 24x144 = 3456 visibility channels, and for 20 integration period. As a result we obtain

Amplitude analysis

Another set of analysis we would like to do when looking at the visibility data is to have a look at the amplitude for various correlation across time and frequency. To do so, we have implemented two similar functions that would generate the amplitude vs freq and vs time analysis figures.

visibilities.amplitude_freq_analysis(
    polarisation,
    nb_time_period,
    coarse_channels,
    sps_correction,
)
visibilities.amplitude_time_analysis(
    polarisation,
    nb_time_period,
    coarse_channels,
    sps_correction,
)

Where:

polarisation is the polarisation of interest, 4 are possible XX, XY, YX, YY
nb_time_period is the number of integration we want to analyse
coarse_channels is the list of coarse channels to plot
sps_correction is the correction to be applied on the data due to ripples from SPS

For example, when we configured an experiment with a CNIC generating traffic equivalent to 6 stations with stations 1, 2, 3, and 4 seeing the same noise source with no delay, we were able to produce the following result by calling:

visibilities.amplitude_freq_analysis("XX", 20, range(171, 171+24), false)

The call above aims at generated the 3D mesh plot for XX polarisation, for 24 coarse channels equivalent to 24x144 = 3456 visibility channels, and for 20 integration period. As a result we obtain

visibilities.amplitude_time_analysis("XX", 20, range(171, 171+24), false)

The call above aims at generated the 3D mesh plot for XX polarisation, for 24 coarse channels equivalent to 24x144 = 3456 visibility channels, and for 20 integration period. As a result we obtain