Basics

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:

In [1]:
# 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:

In [2]:
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:

In [3]:
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:

In [4]:
# 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:

In [5]:
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:

In [6]:
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.

Working with Data

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:

In [7]:
%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:

1) Tuning 1, Pol. X

2) Tuning 1, Pol. Y

3) Tuning 2, Pol. X

4) Tuning 2, Pol. Y

The interface to the DR spectrometer files is slightly different since the data have already been Fourier transformed and averaged:

In [8]:
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:

In [9]:
%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()

Offsetting

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:

In [10]:
# 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.

In [11]:
# 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

Closing

When you are done using an LDP instance, just close it out:

In [12]:
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:

In [13]:
try:
    idfD.read(0.010)
except Exception, e:
    print "ERROR: %s" % str(e)
ERROR: I/O operation on closed file