The LSL interface for DRX (beam formed) data is similar to that of TBN.
First, download a snippet of DRX data:
# This may take a bit...
import os
import urllib2
from lsl.reader import drx
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.)
DRX Size: 1209.4 kB
To read in a chunk of data:
from lsl.reader import drx, errors
# Read in the first 5 frames
fh = open('temp.drx', 'rb')
for i in xrange(5):
try:
frame = drx.readFrame(fh)
except errors.syncError:
continue
except errors.eofError:
break
print "Beam: %i, Tuning: %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()
Beam: 1, Tuning: 1, Pol.: 1 -> Time tag: 1324515670.601744 s -> Mean value: -0.009 + 0.000 j Beam: 1, Tuning: 1, Pol.: 0 -> Time tag: 1324515670.601953 s -> Mean value: -0.018 + 0.010 j Beam: 1, Tuning: 1, Pol.: 1 -> Time tag: 1324515670.601953 s -> Mean value: 0.000 + 0.010 j Beam: 1, Tuning: 1, Pol.: 0 -> Time tag: 1324515670.602162 s -> Mean value: 0.011 + -0.014 j Beam: 1, Tuning: 1, Pol.: 1 -> Time tag: 1324515670.602162 s -> Mean value: -0.020 + -0.006 j
The time tags reported in the above example are in seconds since the Unix epoch. These can be converted to Python datetime instances through:
from datetime import datetime
from lsl.reader import drx, errors
# Read in the first 5 frames
fh = open('temp.drx', 'rb')
for i in xrange(5):
try:
frame = drx.readFrame(fh)
except errors.syncError:
continue
except errors.eofError:
break
print "Beam: %i, Tuning: %i, Pol.: %i" % frame.parseID()
print "-> Time tag: %.6f s (%s)" % (frame.getTime(), datetime.utcfromtimestamp(frame.getTime()))
print "-> Mean value: %.3f + %.3f j" % (frame.data.iq.mean().real, frame.data.iq.mean().imag)
fh.close()
Beam: 1, Tuning: 1, Pol.: 1 -> Time tag: 1324515670.601744 s (2011-12-22 01:01:10.601744) -> Mean value: -0.009 + 0.000 j Beam: 1, Tuning: 1, Pol.: 0 -> Time tag: 1324515670.601953 s (2011-12-22 01:01:10.601953) -> Mean value: -0.018 + 0.010 j Beam: 1, Tuning: 1, Pol.: 1 -> Time tag: 1324515670.601953 s (2011-12-22 01:01:10.601953) -> Mean value: 0.000 + 0.010 j Beam: 1, Tuning: 1, Pol.: 0 -> Time tag: 1324515670.602162 s (2011-12-22 01:01:10.602162) -> Mean value: 0.011 + -0.014 j Beam: 1, Tuning: 1, Pol.: 1 -> Time tag: 1324515670.602162 s (2011-12-22 01:01:10.602162) -> Mean value: -0.020 + -0.006 j
Similar to TBW and TBN data, you can easily create spectra with SpecMaster. The following example reads and processes the first two frames in a DRX file:
%matplotlib inline
import numpy
from lsl.correlator import fx as fxc
from lsl.misc.mathutil import to_dB
from matplotlib import pyplot as plt
fh = open('temp.drx', 'rb')
frame1 = drx.readFrame(fh)
frame2 = drx.readFrame(fh)
srate = frame1.getSampleRate() # Data sample rate in Hz
cFreq = frame1.getCentralFreq() # Data tuning in Hz
# SpecMaster expects 2-D data
data = numpy.zeros((2,frame1.data.iq.size), dtype=frame1.data.iq.dtype)
data[0,:] = frame1.data.iq
data[1,:] = frame2.data.iq
freq, spec = fxc.SpecMaster(data, LFFT=512, SampleRate=srate, CentralFreq=cFreq)
fig = plt.figure()
ax = fig.gca()
ax.plot(freq/1e6, to_dB(spec[0,:]), label='B%i,T%i,P%i' % frame1.parseID())
ax.plot(freq/1e6, to_dB(spec[1,:]), label='B%i,T%i,P%i' % frame2.parseID())
ax.set_xlabel('Frequency [MHz]')
ax.set_ylabel('PSD [arb. dB]')
ax.legend(loc=0)
plt.show()
ERROR: Magic function `matplotlib` not found.
To integrate down you need to read in more data and save it to a NumpyArray. To read in 50 frames from each tuning/polarization use:
data = numpy.zeros((4,50*4096), dtype=numpy.complex64)
count = [0,0,0,0]
for i in xrange(50):
for j in xrange(2):
frame = drx.readFrame(fh)
srate = frame.getSampleRate() # Data sample rate in Hz
cFreq = frame.getCentralFreq() # Data tuning in Hz
beam,tune,pol = frame.parseID()
k = 2*(tune-1) + pol
data[k, count[k]*4096:(count[k]+1)*4096] = frame.data.iq
count[k] += 1
freq, spec = fxc.SpecMaster(data, LFFT=512, SampleRate=srate, CentralFreq=cFreq)
fig = plt.figure()
ax = fig.gca()
ax.plot(freq/1e6, to_dB(spec[0,:]), label='B%i,T%i,P%i' % (beam,tune,0))
ax.plot(freq/1e6, to_dB(spec[1,:]), label='B%i,T%i,P%i' % (beam,tune,1))
ax.set_xlabel('Frequency [MHz]')
ax.set_ylabel('PSD [arb. dB]')
ax.legend(loc=0)
plt.show()
LSL also supports computing all four Stokes parameters from data using the StokesMaster function. Using this function requires a little more setup since the function needs to know which singals are from which beam/tuning/polarization. This is done by first creating "dummy" lsl.common.stations.Antenna instances to store the mapping:
from lsl.common import stations
antennas = []
for i in xrange(2):
if i / 2 == 0:
newAnt = stations.Antenna(1)
else:
newAnt = stations.Antenna(2)
if i % 2 == 0:
newAnt.pol = 0
else:
newAnt.pol = 1
newAnt.digitizer = i+1
antennas.append(newAnt)
for ant in antennas:
print str(ant)
Antenna 1: stand=0, polarization=0; digitizer 1; status is 0 Antenna 1: stand=0, polarization=1; digitizer 2; status is 0
Now the data can be passed into StokesMaster and the results plotted:
freq, spec = fxc.StokesMaster(data, antennas, LFFT=512, SampleRate=srate, CentralFreq=cFreq)
print freq.shape, spec.shape
fig = plt.figure()
ax = fig.gca()
ax.plot(freq/1e6, spec[0,0,:], label='I')
ax.plot(freq/1e6, spec[1,0,:], label='Q')
ax.plot(freq/1e6, spec[2,0,:], label='U')
ax.plot(freq/1e6, spec[3,0,:], label='V')
ax.set_xlabel('Frequency [MHz]')
ax.set_ylabel('PSD [arb. dB]')
ax.legend(loc=0)
plt.show()
(512,) (4, 1, 512)
Note that the shape of the output spectra is 3-D in the case of StokesMaster as opposed to 2-D for SpecMaster. The extra dimension stores the four Stokes parameters in the other of I, Q, U, and V.