Listening to HERA Data

D. Jacobs
16 Dec 2021

Lets make some interesting sounds with HERA data? I've tried this before with PAPER data. There are many ways to do it, they all sound different.

I can't find the script I used on the PAPER data, so I will have to re-remember. This notebook reduces data into a reasonable selection and time ordering suitable for audio.

%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
# import hera_pspec as hp
from pyuvdata import UVData
import os,sys
import glob
# from hera_pspec.data import DATA_PATH
# import copy
# import time
#import uvtools
import wavio

#make plots safe for pasting into dark background
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('retina')

#Lets first try the redcal vsibility solutions. These are nice because the files have no redundant baselines 
IDR3_omnivis_glob = '/lustre/aoc/projects/hera/H1C_IDR3/IDR3_2/2458041/*sum.omni_vis.uvh5'
#this amounts to about 7GB of data

#lets add more data
mychan = 200
uvd = UVData()
files = glob.glob(IDR3_omnivis_glob)
print("reading {n} files".format(n=(len(files))))
uvd.read(files,freq_chans=mychan)

reading 73 files

#how much time and how many baselines do we have
print("Nbls = ",uvd.Nbls,"Ntimes  = ",uvd.Ntimes)
print(uvd.data_array.shape)
#How many data points do we have?
uvd.data_array[:,:,:,0].size
#at a sample rate of 44.1k, how long of an audio recording will we get
rate = 44100
print(uvd.data_array[:,:,:,0].size/rate,'s')

Nbls =  113 Ntimes  =  4354
(492002, 1, 1, 2)
11.156507936507937 s

# select the only spw and one pol (probably X)
D = np.ma.masked_array(uvd.data_array[:,0,:,0],mask=uvd.flag_array[:,0,0,0]) 

# order by frequency, time, baseline
D = D.reshape((uvd.Ntimes,uvd.Nbls))

#go from 3D to 1D. 'F' order goes fastest along the first axis. 
# This should result in each baseline being played sequentially in time order
D_f = D.ravel(order='F')

rate = 44100
wavio.write("HERA_H1C_omni_12h_ch200_32bit_22k_tblf_F.wav",D_f.real,rate/2,sampwidth=4)

plt.figure()
print(uvd.uvw_array[0])
plt.plot(D[:,0])
plt.plot(D_f[:uvd.Ntimes]+0.05)

[ 1.46078427e+01  5.57885268e-02 -9.77542666e-03]

Casting complex values to real discards the imaginary part
Casting complex values to real discards the imaginary part

[<matplotlib.lines.Line2D at 0x7f3b8c398d00>]

#time vs baseline
plt.figure()
plt.imshow(D.real,aspect='auto')

<matplotlib.image.AxesImage at 0x7f3b8ab49550>

Conclusion

Method

The following algorithm makes for interesting sound

1. Pick a channel, use select on read
2. apply the flags
3. reshape to TIME, BLS
4. Ravel in F order. This is the fortran "row fastest" ordering.

This results in each time series for each baseline being played in sequence.

Description

The sounds you hear are the cross correlations between HERA antennas. The signal is 99.999% foregrounds like our galaxy and billions of other galaxies. In fact, its mostly one or two very bright ones, the biggest being Fornax A. The tonality you hear derives from the rate at the sources appear to move. Pairs of antennas experience sky sources as tones with a pitch proportional to dish separation. The sound clip plays pairs of antennas roughly in distance sorted order, close together antennas turn into longer antennas.

The sound itself is a strange wobbly moaning, punctuated by scratchy pops like an old record. I hear a recording from the Palagic Abyss as a Deep Old One activates its arcane machinery.