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:
# 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:
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:
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:
# 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
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:
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
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:
%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.
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):
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:
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:
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)