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')