===== Read raw data (FFTS format) =====
import numpy as np
import os
def getRawData(filename,metainfo=None):
"""read raw data from FFTS into numpy record arrays
metainfo ... if provided (e.g. if known from a previous file
of the same type, this prevents parsing it again
returns alldata,metainfo
alldata contains the following fields (use alldata['FIELD'] to access)
IEEE, DATFMT, DATALENGTH, BE_IDENT, TIMESTAMP, INTTIME, PHASENUM, BENUMSEC,
BLOCKINGFACTOR, DATABLOCK
DATABLOCK is itself a record, containing:
SECTION, NUMCHANNELS, DATA
to access the raw data use 'alldata['DATABLOCK']['DATA']' which is a 3D array,
of shape (dump,baseband,specchannel), e.g., to access first dump, first baseband
spectrum 'alldata['DATABLOCK']['DATA'][0,0]', last dump, second baseband:
'alldata['DATABLOCK']['DATA'][-1,1]', etc.
use:
import dateutil.parser
alltimestamps=[dateutil.parser.parse(a) for a in alldata['TIMESTAMP']]
allmjd=np.array([cal_mjd(dt.month,
dt.day+(dt.hour+dt.minute/60.+dt.second/3600.+dt.microsecond/3.6e9)/24.,
dt.year) for dt in alltimestamps])
to compute MJD values from the timestamps
"""
totalsize=os.path.getsize(filename)
if metainfo is None:
# getting endian
ieee=np.fromfile(filename , dtype = np.dtype([('IEEE','|S4')]) , count = 1 )
endian='>' if ieee=="IEEE" else '<'
metainfo=np.fromfile(filename , dtype = np.dtype([ ('IEEE','|S4'), ('DATFMT','|S4'),
('DATALENGTH','%si4'%endian), ('BE_IDENT','|S8'), ('TIMESTAMP','|S28'),
('INTTIME','%si4'%endian), ('PHASENUM','%si4'%endian), ('BENUMSEC','%si4'%endian),
('BLOCKINGFACTOR','%si4'%endian), ('SECTION','%si4'%endian),
('NUMCHANNELS','%si4'%endian) ]) , count = 1 )
else:
endian='>' if metainfo['IEEE'][0]=="IEEE" else '<'
#
intorfloat=np.int32 if metainfo['DATFMT'][0]=="I" else np.float32
datfmt=('%si4'%endian) if metainfo['DATFMT'][0]=="I" else ('%sf4'%endian)
datalength=metainfo['DATALENGTH'][0]
numchannels=metainfo['NUMCHANNELS'][0]
datablock=np.dtype([('SECTION','%si4'%endian), ('NUMCHANNELS','%si4'%endian),
('DATA',datfmt,numchannels)])
record=np.dtype([ ('IEEE','|S4'), ('DATFMT','|S4'), ('DATALENGTH','%si4'%endian),
('BE_IDENT','|S8'), ('TIMESTAMP','|S28'), ('INTTIME','%si4'%endian),
('PHASENUM','%si4'%endian), ('BENUMSEC','%si4'%endian), ('BLOCKINGFACTOR','%si4'%endian),
('DATABLOCK',datablock,metainfo['BLOCKINGFACTOR'][0]*metainfo['BENUMSEC'][0]) ])
record_length=record.itemsize
num_records=totalsize/record_length
alldata=np.fromfile(filename , dtype =record, count = num_records )
return alldata,metainfo
#
fname="/be4-daten/AFFTS/2788.0001"
alldata,metainfo=getRawData(fname)
print alldata['DATABLOCK']['SECTION']
print alldata['DATABLOCK']['NUMCHANNELS']
print alldata['DATABLOCK']['DATA'].shape
If you want to convert the timestamps to MJD, you can do the following:
import dateutil.parser
def cal_mjd(mn,dy,yr):
"""calculate MJD from a data
mn ... month (int)
yr ... year (int)
dy ... day (float, i.e. including fractions of a day)
returns MJD
(code snippet adapted from pyephem, which is using libastro)"""
m=int(mn)
y=int(yr)
if y<0: y+=1
if mn<3: mn+=12
if mn<3: y-=1
b=0
if not ( (yr < 1582) | ((yr == 1582) & ((mn < 10) | ((mn == 10) & (dy < 15))))):
a=y//100
b=2 - a + a//4
c=0
if y<0:
c=int((365.25*y) - 0.75) - 694025L
else:
c=int(365.25*y) - 694025L
d=int(30.6001*(m+1))
return b + c + d + dy - 0.5 +15020L -0.5
#
alltimestamps=[dateutil.parser.parse(a) for a in alldata['TIMESTAMP']]
allmjd=np.array([cal_mjd(dt.month,
dt.day+(dt.hour+dt.minute/60.+dt.second/3600.+dt.microsecond/3.6e9)/24.,
dt.year) for dt in alltimestamps])