Basics

LSL provides readers for the two DP transient buffer outputs: the transient buffer - wideband (TBW) and the transient buffer - narrowband (TBW) through the lsl.reader module. First, download a small portion of a TBW data file:

In [5]:
# This may take a bit...
import os
import urllib2
from lsl.reader import tbw

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.)
TBW Size: 1195.3 kB

To read in the TBW data:

In [5]:
from lsl.reader import tbw, errors

fh = open('temp.tbw', 'rb')
for i in xrange(5):
    try:
        frame = tbw.readFrame(fh)
    except errors.syncError:
        continue
    except errors.eofError:
        break
        
    print "DP Stand Number: %i" % frame.parseID()
    print "-> Frame Count: %i" % frame.header.frameCount
    print "-> Time tag: %.6f s" % frame.getTime()
    print "-> Mean value: %.3f" % frame.data.xy.mean()
fh.close()
DP Stand Number: 22
-> Frame Count: 1
-> Time tag: 1302134255.000000 s
-> Mean value: 12.954
DP Stand Number: 21
-> Frame Count: 2
-> Time tag: 1302134255.000002 s
-> Mean value: 14.061
DP Stand Number: 22
-> Frame Count: 2
-> Time tag: 1302134255.000002 s
-> Mean value: 8.617
DP Stand Number: 102
-> Frame Count: 1
-> Time tag: 1302134255.000000 s
-> Mean value: 7.439
DP Stand Number: 101
-> Frame Count: 2
-> Time tag: 1302134255.000002 s
-> Mean value: 3.271

The try...except blocks help deal with problems in the file as the come up and keep things running smoothly.

One of the most confusing aspects of this interface is that the stand values returned by the "frame.parseID" method are DP stand numbers, which are used internally in DP to keep track of signals, and not LWA1 stand numbers. However, the DP stand numbers can be related to LWA1 stand numbers using the lsl.common.stations.lwa1 instance:

In [7]:
from lsl.common.stations import lwa1
from lsl.reader import tbw, errors

# Get the antenna list
antennas = lwa1.getAntennas()

# Read in the first 5 frames
fh = open('temp.tbw', 'rb')
for i in xrange(5):
    try:
        frame = tbw.readFrame(fh)
    except errors.syncError:
        continue
    except errors.eofError:
        break
        
    # Calculate the digitizer
    dig = 2*(frame.parseID()-1) + 1
    
    print "DP Stand Number: %i" % frame.parseID()
    print "-> LWA1 Stand: %i" % antennas[dig-1].stand.id
    print "-> Frame Count: %i" % frame.header.frameCount
    print "-> Time tag: %.6f s" % frame.getTime()
    print "-> Mean value: %.3f" % frame.data.xy.mean()
fh.close()
DP Stand Number: 22
-> LWA1 Stand: 85
-> Frame Count: 1
-> Time tag: 1302134255.000000 s
-> Mean value: 12.954
DP Stand Number: 21
-> LWA1 Stand: 242
-> Frame Count: 2
-> Time tag: 1302134255.000002 s
-> Mean value: 14.061
DP Stand Number: 22
-> LWA1 Stand: 85
-> Frame Count: 2
-> Time tag: 1302134255.000002 s
-> Mean value: 8.617
DP Stand Number: 102
-> LWA1 Stand: 112
-> Frame Count: 1
-> Time tag: 1302134255.000000 s
-> Mean value: 7.439
DP Stand Number: 101
-> LWA1 Stand: 187
-> Frame Count: 2
-> Time tag: 1302134255.000002 s
-> Mean value: 3.271

The interface for TBN data is similar. To download a small section of data:

In [8]:
# This may take a bit...
import os
import urllib2
from lsl.reader import tbn

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.)
TBN Size: 1023.4 kB
In [29]:
from lsl.reader import tbn, errors

# Read in the first 5 frames
fh = open('temp.tbn', 'rb')
for i in xrange(5):
    try:
        frame = tbn.readFrame(fh)
    except errors.syncError:
        continue
    except errors.eofError:
        break
        
    print "DP Stand Number: %i, pol. %i" % frame.parseID()
    print "-> Time tag: %.6f s" % frame.getTime()
    print "-> Mean value: %.3f + %.3f j" % (frame.data.iq.mean().real, frame.data.iq.mean().imag)
fh.close()
DP Stand Number: 101, pol. 0
-> Time tag: 1330573596.004160 s
-> Mean value: -0.254 + 0.369 j
DP Stand Number: 191, pol. 0
-> Time tag: 1330573596.004160 s
-> Mean value: 0.064 + -0.055 j
DP Stand Number: 151, pol. 0
-> Time tag: 1330573596.004160 s
-> Mean value: -0.371 + 0.217 j
DP Stand Number: 81, pol. 0
-> Time tag: 1330573596.004160 s
-> Mean value: -0.188 + -0.326 j
DP Stand Number: 171, pol. 0
-> Time tag: 1330573596.004160 s
-> Mean value: -0.158 + 0.129 j

Similar to the TBW case, the stand and polarizations returned by "frame.parseID" are from the interal DP structure. To relate this to LWA1 stand/polarizations you use a similar method:

In [2]:
from lsl.common.stations import lwa1
from lsl.reader import tbn, errors

# Download a small bit of TBN data (1000 frames to be exact)
# This may take a bit...
import os
import urllib2
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()
    
# Get the antenna list
antennas = lwa1.getAntennas()

# Read in the first 5 frames
fh = open('temp.tbn', 'rb')
for i in xrange(5):
    try:
        frame = tbn.readFrame(fh)
    except errors.syncError:
        continue
    except errors.eofError:
        break
        
    # Calculate the digitizer
    stand,pol = frame.parseID()
    dig = 2*(stand-1) + pol + 1
    
    print "DP Stand Number: %i, pol. %i" % (stand,pol)
    print "-> LWA1 Stand, pol.: %i, %i" % (antennas[dig-1].stand.id, antennas[dig-1].pol)
    print "-> Time tag: %.6f s" % frame.getTime()
    print "-> Mean value: %.3f + %.3f j" % (frame.data.iq.mean().real, frame.data.iq.mean().imag)
fh.close()
DP Stand Number: 101, pol. 0
-> LWA1 Stand, pol.: 187, 0
-> Time tag: 1330573596.004160 s
-> Mean value: -0.254 + 0.369 j
DP Stand Number: 191, pol. 0
-> LWA1 Stand, pol.: 13, 0
-> Time tag: 1330573596.004160 s
-> Mean value: 0.064 + -0.055 j
DP Stand Number: 151, pol. 0
-> LWA1 Stand, pol.: 90, 0
-> Time tag: 1330573596.004160 s
-> Mean value: -0.371 + 0.217 j
DP Stand Number: 81, pol. 0
-> LWA1 Stand, pol.: 125, 0
-> Time tag: 1330573596.004160 s
-> Mean value: -0.188 + -0.326 j
DP Stand Number: 171, pol. 0
-> LWA1 Stand, pol.: 130, 0
-> Time tag: 1330573596.004160 s
-> Mean value: -0.158 + 0.129 j

Generating and Plotting Spectra

Although computing the mean I/Q value for each frame isn't particularlly interesting this example shows how data can be extracted from within each TBN frame. A more complicated example that generates spectra would look like:

In [37]:
%matplotlib inline
from lsl.misc.mathutil import to_dB
from matplotlib import pyplot as plt

fh = open('temp.tbn', 'rb')
srate = tbn.getSampleRate(fh)   # Data sample rate in Hz
frame = tbn.readFrame(fh)
cFreq = frame.getCentralFreq()  # Data tuning in Hz
frame.data.iq.shape = (1,frame.data.iq.size)  # SpecMaster expects 2-D data

from lsl.correlator import fx as fxc
freq, spec = fxc.SpecMaster(frame.data.iq, LFFT=512, SampleRate=srate, CentralFreq=cFreq) 
print freq.shape, spec.shape

fig = plt.figure()
ax = fig.gca()
ax.plot((freq-freq.mean())/1e3, to_dB(spec[0,:]))
ax.set_xlabel('Relative Freq [kHz]')
ax.set_ylabel('PSD [arb. dB]')
plt.show()
(511,) (1, 511)

The "SpecMaster" function used above also works for TBW data.

Post-Acquisition Beam Formings with TBN

For post-acquisition beam forming, you need need an azimuth (in degrees) and elevation (in degrees) to point the beam towards. For planets, this can be accomplished using the pyephem package that is required by lsl. For example, compute the location of Jupiter at LWA-1 on 12/17/2010 at 21:18 UTC (JD 2,455,548.38787):

In [20]:
import math
import ephem
from lsl.common import stations

lwa1 = stations.lwa1
lwaObserver = lwa1.getObserver(2455548.38787, JD=True)
jove = ephem.Jupiter()
jove.compute(lwaObserver)
print "Jupiter:  az -> %.1f, el -> %.1f" % (jove.az*180/math.pi, jove.alt*180/math.pi)
Jupiter:  az -> 112.4, el -> 24.4

For fixed positions, use:

In [35]:
lwaObserver.date = '2010/12/17 21:18:34'

cyga = ephem.FixedBody()
cyga._ra = '19:59:28.30'
cyga._dec = '+40:44:02'
cyga.compute(lwaObserver)
print "Cygnus A:  az -> %.1f, el -> %.1f" % (cyga.az*180/math.pi, cyga.alt*180/math.pi)
Cygnus A:  az -> 10.0, el -> 83.2

After TBN data have been read in and a pointing position has been found, a beam can be formed through phase-and-sum beamforming:

In [39]:
from lsl.misc import beamformer

antennas = [lwa1.getAntennas()[0],]

beamdata = beamformer.phaseAndSum(antennas, frame.data.iq, sampleRate=1e5, azimuth=10.9, elevation=83.2)
print beamdata.shape
(1, 512)