LSL also provides a high-level interface to LWA1 data through the LSL Developer Primitives (LDP) module. First, download a small portion of the various LWA1 data formats:
# This may take a bit...
import os
import urllib2
from lsl.reader import tbw, tbn, drx
if not os.path.exists('temp.tbw'):
fh1 = urllib2.urlopen('http://fornax.phys.unm.edu/lwa/data/firstlight/firstLightTBW.dat')
fh2 = open('temp.tbw', 'wb')
fh2.write(fh1.read(tbw.FrameSize*1000))
fh1.close()
fh2.close()
print "TBW Size: %.1f kB" % (os.path.getsize('temp.tbw')/1024.)
if not os.path.exists('temp.tbn'):
fh1 = urllib2.urlopen('http://fornax.phys.unm.edu/lwa/data/TBN/055987_000496599_DEFAULT_TBN_00000_30.dat')
fh2 = open('temp.tbn', 'wb')
fh2.write(fh1.read(tbn.FrameSize*1000))
fh1.close()
fh2.close()
print "TBN Size: %.1f kB" % (os.path.getsize('temp.tbn')/1024.)
if not os.path.exists('temp.drx'):
fh1 = urllib2.urlopen('http://lda10g.alliance.unm.edu/data1/055917_000007389_Jupiter.dat')
fh2 = open('temp.drx', 'wb')
fh2.write(fh1.read(drx.FrameSize*300))
fh1.close()
fh2.close()
print "DRX Size: %.1f kB" % (os.path.getsize('temp.drx')/1024.)
if not os.path.exists('temp.drspec'):
fh1 = urllib2.urlopen('http://lda10g.alliance.unm.edu/tutorial/B0329+54/056770_000044687')
fh2 = open('temp.drspec', 'wb')
fh2.write(fh1.read(drx.FrameSize*300))
fh1.close()
fh2.close()
print "DR Spectrometer Size: %.1f kB" % (os.path.getsize('temp.drspec')/1024.)
TBW Size: 1195.3 kB TBN Size: 1023.4 kB DRX Size: 1209.4 kB DR Spectrometer Size: 1209.4 kB
Getting started with LDP is as easy as:
from lsl.reader.ldp import LWA1DataFile
The LWA1DataFile() function takes either a filename or an open file handle and uses the contents of the file to select the correct data format handler. For example:
idfD = LWA1DataFile('temp.drx')
print "'temp.drx' is of type:", type(idfD)
idfS = LWA1DataFile('temp.drspec')
print "'temp.drspec' is of type:", type(idfS)
fh = open('temp.tbw', 'rb')
idfW = LWA1DataFile(fh=fh)
print "'temp.tbw' is of type:", type(idfW)
idfN = LWA1DataFile('temp.tbn')
print "'temp.tbn' is of type:", type(idfN)
'temp.drx' is of type: <class 'lsl.reader.ldp.DRXFile'> 'temp.drspec' is of type: <class 'lsl.reader.ldp.DRSpecFile'> 'temp.tbw' is of type: <class 'lsl.reader.ldp.TBWFile'> 'temp.tbn' is of type: <class 'lsl.reader.ldp.TBNFile'>
As part of the initialization process LDP also takes care of preparing the file for reading, such as aligning the time tags for DRX files, and creating the ring buffer for TBN and DRX files.
LDP provides access to the file metadata through the getInfo() method. Called without an argument, getInfo() returns a dictionary of various metadata parameters that have been determined. Specific metadata values can be returned by specifying the parameter:
# Poll the TBN file for its specifics
print "TBN Metadata:"
for key,value in idfN.getInfo().iteritems():
print " %s:" % key, value
print " "
# Poll the DRX file for the sample rate
print "DRX Sample Rate %.3f MS/s" % (idfD.getInfo('sampleRate')/1e6,)
print " "
# Poll the DR spectrometer file for its specifics
print "DR Spectrometer Metadata:"
for key,value in idfS.getInfo().iteritems():
print " %s:" % key, value
print " "
TBN Metadata: dataBits: 8 FrameSize: 1048 tStart: 1330573596.0 nAntenna: 520 sampleRate: 100000.0 nFrames: 1000 freq1: 74030000.4482 size: 1048000 DRX Sample Rate 19.600 MS/s DR Spectrometer Metadata: tInt: 0.0401240816327 FrameSize: 16460 dataBits: 32 nFrames: 75 beam: 2 nProducts: 2 tStart: 1398276221.03 beampols: 4 nInt: 768 dataProducts: ['XX', 'YY'] freq2: 73999999.7616 freq1: 54400000.9298 LFFT: 1024 sampleRate: 19600000.0 size: 1238400
To access the data through LDP, the most convenient method is to use the read() method of the LDP instance returned by LWA1DataFile(). This method takes in a single argument for the amount of time to read in and returns a multi-element tuple of the timestamp for the first data sample, the duration of data returned, and the data itself as a NumPy array. For example, to read in 10 ms of raw data from a DRX file:
duration, tStart, data = idfD.read(0.010)
print "Duration read: %.3f ms" % (duration*1e3,)
print "Time of first sample: %.6f" % tStart
print "Data specifics:", data.shape, data.dtype
Duration read: 10.031 ms Time of first sample: 1324515670.601953 Data specifics: (4, 196608) complex64
Here the duration variable holds the exact amount of data read in, the tStart variable holds the UNIX timestamp for the first sample, and the data variable is an array of the I/Q values. The timestamp can be converted to a more easily understandable format using the datetime module:
from datetime import datetime
tStartDT = datetime.utcfromtimestamp(tStart)
print tStartDT
2011-12-22 01:01:10.601953
LDP also keeps track of the flow of time inside a data file and will raise a RuntimeError if a jump is found in the file while reading. The behavior can be disabled by setting the "ignoreTimeTagErrors" keyword to True.
The data can also be worked with in a similar fashion to that was presented in the TBN and DRX tutorials. The LDP data handlers always return data that is explicitly ordered. In the case of TBW and TBN data the order is by the digitizer number and for DRX it is by tuning/polarization. For DRX data:
%matplotlib inline
import numpy
from lsl.correlator import fx as fxc
from lsl.misc.mathutil import to_dB
from matplotlib import pyplot as plt
# Compute the spectra
freq, spec = fxc.SpecMaster(data, LFFT=512, SampleRate=idfD.getInfo('sampleRate'), CentralFreq=0.0)
freq1 = freq + idfD.getInfo('freq1')
freq2 = freq + idfD.getInfo('freq2')
# Plot
fig = plt.figure()
ax = fig.gca()
ax.plot(freq1/1e6, to_dB(spec[0,:]), label='B%i,T1,P0' % idfD.getInfo('beam'))
ax.plot(freq1/1e6, to_dB(spec[1,:]), label='B%i,T1,P1' % idfD.getInfo('beam'))
ax.set_xlabel('Frequency [MHz]')
ax.set_ylabel('PSD [arb. dB]')
ax.legend(loc=0)
plt.show()
Here the order along the first axis of data is:
Tuning 1, Pol. X
Tuning 1, Pol. Y
Tuning 2, Pol. X
Tuning 2, Pol. Y
The interface to the DR spectrometer files is slightly different since the data have already been Fourier transformed and averaged:
durationS, tStartS, dataS = idfS.read(0.100)
print "Duration read: %.3f ms" % (durationS*1e3,)
print "Time of first sample: %.6f" % tStartS
print "Data specifics:", dataS.shape, dataS.dtype
Duration read: 80.248 ms Time of first sample: 1398276221.034108 Data specifics: (4, 2, 1024) float32
The data array here is three dimensional with dimensions that index the tuning/polarization pair, the integrations, and the channels. In order to have the plot contain frequencies rather than channel numbers there is an additional step required to determine the frequency bins. This is accomplished through the fftfreq and fftshift functions in the numpy.fft module:
%matplotlib inline
# Get the metadata
b = idfS.getInfo('beam')
dp = idfS.getInfo('dataProducts')
# Compute the frequency bins
freq = numpy.fft.fftfreq(dataS.shape[2], d=1.0/idfS.getInfo('sampleRate'))
freq = numpy.fft.fftshift(freq)
freq1 = freq + idfS.getInfo('freq1')
freq2 = freq + idfS.getInfo('freq2')
# Intgrate
spec = dataS.mean(axis=1)
# Plot
fig = plt.figure()
ax = fig.gca()
ax.plot(freq1/1e6, to_dB(spec[0,:]), label='B%i,T1,%s' % (b, dp[0 % len(dp)]))
ax.plot(freq1/1e6, to_dB(spec[1,:]), label='B%i,T1,%s' % (b, dp[1 % len(dp)]))
ax.plot(freq2/1e6, to_dB(spec[2,:]), label='B%i,T1,%s' % (b, dp[2 % len(dp)]))
ax.plot(freq2/1e6, to_dB(spec[3,:]), label='B%i,T1,%s' % (b, dp[3 % len(dp)]))
ax.set_xlabel('Frequency [MHz]')
ax.set_ylabel('PSD [arb. dB]')
ax.legend(loc=0)
plt.show()
The LDP handlers also provide an easy-to-use mechanism to skip forward within a data file using the offset() method. This method takes an offset in seconds and adjusts the location of the read() method within the file accordingly. For example:
# Offset 10 ms into the file from the current location
skip = idfD.offset(0.010)
print "Skipped: %.3f ms" % (skip*1e3,)
# Read in another 50 ms of data
duration2, tStart2, data2 = idfD.read(0.050)
print "Duration read: %.3f ms" % (duration*1e3,)
print "Time of first sample: %.6f" % tStart
print "Data specifics:", data.shape, data.dtype
print " "
print "Time skip: %.3f ms" % ((tStart2-tStart)*1e3,)
Skipped: 9.822 ms Duration read: 10.031 ms Time of first sample: 1324515670.601953 Data specifics: (4, 196608) complex64 Time skip: 21.734 ms
At any point you can also skip back to the beginning of the file using the reset() method.
# Go back to the beginning
idfD.reset()
# Read in the first 10 ms again
duaration3, tStart3, data3 = idfD.read(0.010)
# Are the start times the same?
print "Time good?", True if tStart3-tStart == 0 else False
# Are the data the same?
print "Data good?", True if (numpy.abs(data-data3)**2).sum() == 0 else False
Time good? True Data good? True
When you are done using an LDP instance, just close it out:
idfD.close()
idfS.close()
idfN.close()
idfW.close()
Once a LDP instance is closed you will not be able to access the data in the file:
try:
idfD.read(0.010)
except Exception, e:
print "ERROR: %s" % str(e)
ERROR: I/O operation on closed file