Initial commit

.gitignore

+*.dat
+*.fil
+*.pyc

README.md

+This contains some scripts to add metadata to the camras pulsar format, i.e. the format that comes out of `pulsar_record`.

filterbank.py

+"""
+A module for reading filterbank files.
+
+Patrick Lazarus, June 26, 2012
+(Minor modification from file originally from June 6th, 2009)
+"""
+
+import sys
+import warnings
+import os
+import os.path
+import numpy as np
+import sigproc
+
+
+DEBUG = False
+
+def create_filterbank_file(outfn, header, spectra=None, nbits=8, \
+                           verbose=False, mode='append'):
+    """Write filterbank header and spectra to file.
+
+        Input:
+            outfn: The outfile filterbank file's name.
+            header: A dictionary of header paramters and values.
+            spectra: Spectra to write to file. (Default: don't write
+                any spectra - i.e. write out header only)
+            nbits: The number of bits per sample of the filterbank file.
+                This value always overrides the value in the header dictionary.
+                (Default: 8 - i.e. each sample is an 8-bit integer)
+            verbose: If True, be verbose (Default: be quiet)
+            mode: Mode for writing (can be 'append' or 'write')
+
+        Output:
+            fbfile: The resulting FilterbankFile object opened
+                in read-write mode.
+    """
+    dtype = get_dtype(nbits) # Get dtype. This will check to ensure
+                             # 'nbits' is valid.
+    header['nbits'] = nbits
+    outfile = open(outfn, 'wb')
+    outfile.write(sigproc.addto_hdr("HEADER_START", None))
+    for paramname in list(header.keys()):
+        if paramname not in sigproc.header_params:
+            # Only add recognized parameters
+            continue
+        if verbose:
+            print("Writing header param (%s)" % paramname)
+        value = header[paramname]
+        outfile.write(sigproc.addto_hdr(paramname, value))
+    outfile.write(sigproc.addto_hdr("HEADER_END", None))
+    if spectra is not None:
+        spectra.flatten().astype(dtype).tofile(outfile)
+    outfile.close()
+    return FilterbankFile(outfn, mode=mode)
+
+
+def is_float(nbits):
+    """For a given number of bits per sample return
+        true if it corresponds to floating-point samples
+        in filterbank files.
+
+        Input:
+            nbits: Number of bits per sample, as recorded in the filterbank
+                file's header.
+
+        Output:
+            isfloat: True, if 'nbits' indicates the data in the file
+                are encoded as floats.
+    """
+    check_nbits(nbits)
+    if nbits == 32:
+        return True
+    else:
+        return False
+
+
+def check_nbits(nbits):
+    """Given a number of bits per sample check to make
+        sure 'filterbank.py' can cope with it.
+
+        An exception is raise if 'filterbank.py' cannot cope.
+
+        Input:
+            nbits: Number of bits per sample, as recorded in the filterbank
+                file's header.
+
+        Output:
+            None
+    """
+    if nbits not in [32, 16, 8]:
+        raise ValueError("'filterbank.py' only supports " \
+                                    "files with 8- or 16-bit " \
+                                    "integers, or 32-bit floats " \
+                                    "(nbits provided: %g)!" % nbits)
+
+
+def get_dtype(nbits):
+    """For a given number of bits per sample return
+        a numpy-recognized dtype.
+
+        Input:
+            nbits: Number of bits per sample, as recorded in the filterbank
+                file's header.
+
+        Output:
+            dtype: A numpy-recognized dtype string.
+    """
+    check_nbits(nbits)
+    if is_float(nbits):
+        dtype = 'float%d' % nbits
+    else:
+        dtype = 'uint%d' % nbits
+    return dtype
+
+
+def read_header(filename, verbose=False):
+    """Read the header of a filterbank file, and return
+        a dictionary of header paramters and the header's
+        size in bytes.
+
+        Inputs:
+            filename: Name of the filterbank file.
+            verbose: If True, be verbose. (Default: be quiet)
+
+        Outputs:
+            header: A dictionary of header paramters.
+            header_size: The size of the header in bytes.
+    """
+    header = {}
+    filfile = open(filename, 'rb')
+    filfile.seek(0)
+    paramname = ""
+    while (paramname != 'HEADER_END'):
+        if verbose:
+            print("File location: %d" % filfile.tell())
+        paramname, val = sigproc.read_hdr_val(filfile, stdout=verbose)
+        if verbose:
+            print("Read param %s (value: %s)" % (paramname, val))
+        if paramname not in ["HEADER_START", "HEADER_END"]:
+            header[paramname] = val
+    header_size = filfile.tell()
+    filfile.close()
+    return header, header_size
+
+
+class FilterbankFile(object):
+    def __init__(self, filfn, mode='readonly'):
+        self.filename = filfn
+        self.filfile = None
+        if not os.path.isfile(filfn):
+            raise ValueError("ERROR: File does not exist!\n\t(%s)" % filfn)
+        self.header, self.header_size = read_header(self.filename)
+        self.frequencies = self.fch1 + self.foff*np.arange(self.nchans)
+        self.is_hifreq_first = (self.foff < 0)
+        self.bytes_per_spectrum = self.nchans*self.nbits / 8
+        data_size = os.path.getsize(self.filename)-self.header_size
+        self.nspec = data_size/self.bytes_per_spectrum
+       
+        # Check if this file is a folded-filterbank file
+        if 'npuls' in self.header and 'period' in self.header and \
+                'nbins' in self.header and 'tsamp' not in self.header:
+            # Foleded file
+            self.isfold = True
+            self.dt = self.period/self.nbins
+        else:
+            self.isfold = False
+            self.dt = self.tsamp
+
+        # Get info about dtype
+        self.dtype = get_dtype(self.nbits)
+        if is_float(self.nbits):
+            tinfo = np.finfo(self.dtype)
+        else:
+            tinfo = np.iinfo(self.dtype)
+        self.dtype_min = tinfo.min
+        self.dtype_max = tinfo.max
+
+        if mode.lower() in ('read', 'readonly'):
+            self.filfile = open(self.filename, 'rb')
+        elif mode.lower() in ('write', 'readwrite'):
+            self.filfile = open(self.filename, 'r+b')
+        elif mode.lower() == 'append':
+            self.filfile = open(self.filename, 'a+b')
+        else:
+            raise ValueError("Unrecognized mode (%s)!" % mode)
+
+    def close(self):
+        if self.filfile is not None:
+            self.filfile.close()
+
+    def get_timeslice(self, start, stop):
+        startspec = int(np.round(start/self.tsamp))
+        stopspec = int(np.round(stop/self.tsamp))
+        return self.get_spectra(startbins, stopbins)
+
+    def get_spectra(self, start, stop):
+        stop = min(stop, self.nspec)
+        pos = self.header_size+start*self.bytes_per_spectrum
+        # Compute number of elements to read
+        nspec = int(stop) - int(start)
+        num_to_read = nspec*self.nchans
+        num_to_read = max(0, num_to_read)
+        self.filfile.seek(pos, os.SEEK_SET)
+        spectra = np.fromfile(self.filfile, dtype=self.dtype, 
+                              count=num_to_read)
+        spectra.shape = nspec, self.nchans 
+        return spectra
+
+    def append_spectra(self, spectra):
+        """Append spectra to the file if is not read-only.
+            
+            Input:
+                spectra: The spectra to append. The new spectra
+                    must have the correct number of channels (ie
+                    dimension of axis=1.
+
+            Outputs:
+                None
+        """
+        if self.filfile.mode.lower() in ('r', 'rb'):
+            raise ValueError("FilterbankFile object for '%s' is read-only." % \
+                        self.filename)
+        nspec, nchans = spectra.shape
+        if nchans != self.nchans:
+            raise ValueError("Cannot append spectra. Incorrect shape. " \
+                        "Number of channels in file: %d; Number of " \
+                        "channels in spectra to append: %d" % \
+                        (self.nchans, nchans))
+        data = spectra.flatten()
+        np.clip(data, self.dtype_min, self.dtype_max, out=data)
+        # Move to end of file
+        self.filfile.seek(0, os.SEEK_END)
+        self.filfile.write(data.astype(self.dtype))
+        self.nspec += nspec
+        #self.filfile.flush()
+        #os.fsync(self.filfile)
+
+    def write_spectra(self, spectra, ispec):
+        """Write spectra to the file if is writable.
+            
+            Input:
+                spectra: The spectra to append. The new spectra
+                    must have the correct number of channels (ie
+                    dimension of axis=1.
+                ispec: The index of the spectrum of where to start writing.
+
+            Outputs:
+                None
+        """
+        if 'r+' not in self.filfile.mode.lower():
+            raise ValueError("FilterbankFile object for '%s' is not writable." % \
+                        self.filename)
+        nspec, nchans = spectra.shape
+        if nchans != self.nchans:
+            raise ValueError("Cannot write spectra. Incorrect shape. " \
+                        "Number of channels in file: %d; Number of " \
+                        "channels in spectra to write: %d" % \
+                        (self.nchans, nchans))
+        if ispec > self.nspec:
+            raise ValueError("Cannot write past end of file! " \
+                             "Present number of spectra: %d; " \
+                             "Requested index of write: %d" % \
+                             (self.nspec, ispec))
+        data = spectra.flatten()
+        np.clip(data, self.dtype_min, self.dtype_max, out=data)
+        # Move to requested position
+        pos = self.header_size + ispec*self.bytes_per_spectrum
+        self.filfile.seek(pos, os.SEEK_SET)
+        self.filfile.write(data.astype(self.dtype))
+        if nspec+ispec > self.nspec:
+            self.nspec = nspec+ispec
+
+    def __getattr__(self, name):
+        if name in self.header:
+            if DEBUG:
+                print("Fetching header param (%s)" % name)
+            val = self.header[name]
+        else:
+            raise ValueError("No FilterbankFile attribute called '%s'" % name)
+        return val
+
+    def print_header(self):
+        """Print header parameters and values.
+        """
+        for param in sorted(self.header.keys()):
+            if param in ("HEADER_START", "HEADER_END"):
+                continue
+            print("%s: %s" % (param, self.header[param]))
+
+
+def main():
+    fil = FilterbankFile(sys.argv[1])
+    fil.print_header()
+
+
+if __name__ == '__main__':
+    main()

psr_constants.py

+## Automatically adapted for numpy Apr 14, 2006 by convertcode.py
+
+ARCSECTORAD = float('4.8481368110953599358991410235794797595635330237270e-6')
+RADTOARCSEC = float('206264.80624709635515647335733077861319665970087963')
+SECTORAD    = float('7.2722052166430399038487115353692196393452995355905e-5')
+RADTOSEC    = float('13750.987083139757010431557155385240879777313391975')
+RADTODEG    = float('57.295779513082320876798154814105170332405472466564')
+DEGTORAD    = float('1.7453292519943295769236907684886127134428718885417e-2')
+RADTOHRS    = float('3.8197186342054880584532103209403446888270314977710')
+HRSTORAD    = float('2.6179938779914943653855361527329190701643078328126e-1')
+PI          = float('3.1415926535897932384626433832795028841971693993751')
+TWOPI       = float('6.2831853071795864769252867665590057683943387987502')
+PIBYTWO     = float('1.5707963267948966192313216916397514420985846996876')
+SECPERDAY   = float('86400.0')
+SECPERJULYR = float('31557600.0')
+KMPERPC     = float('3.0856776e13')
+KMPERKPC    = float('3.0856776e16')
+Tsun        = float('4.925490947e-6') # sec
+Msun        = float('1.9891e30')      # kg
+Mjup        = float('1.8987e27')      # kg
+Rsun        = float('6.9551e8')       # m
+Rearth      = float('6.378e6')        # m
+SOL         = float('299792458.0')    # m/s
+MSUN        = float('1.989e+30')      # kg
+G           = float('6.673e-11')      # m^3/s^2/kg 
+C           = SOL

sigproc.py

+#!/usr/bin/env python
+import os
+import struct
+import sys
+import math
+import warnings
+from psr_constants import ARCSECTORAD
+
+telescope_ids = {"Fake": 0, "Arecibo": 1, "ARECIBO 305m": 1, 
+                 "Ooty": 2, "Nancay": 3,
+                 "Parkes": 4, "Jodrell": 5, "GBT": 6, "GMRT": 7,
+                 "Effelsberg": 8, "ATA": 9, "UTR-2": 10, "LOFAR": 11}
+ids_to_telescope = dict(list(zip(list(telescope_ids.values()), list(telescope_ids.keys()))))
+
+machine_ids = {"FAKE": 0, "PSPM": 1, "Wapp": 2, "WAPP": 2, "AOFTM": 3,
+               "BCPM1": 4, "BPP": 4, "OOTY": 5, "SCAMP": 6, 
+               "GBT Pulsar Spigot": 7, 
+               "SPIGOT": 7, "BG/P": 11, "PDEV": 12}
+ids_to_machine = dict(list(zip(list(machine_ids.values()), list(machine_ids.keys()))))
+
+header_params = {
+    "HEADER_START": 'flag',
+    "telescope_id": 'i',
+    "machine_id": 'i',
+    "data_type": 'i', 
+    "rawdatafile": 'str',
+    "source_name": 'str', 
+    "barycentric": 'i', 
+    "pulsarcentric": 'i', 
+    "az_start": 'd',  
+    "za_start": 'd',  
+    "src_raj": 'd',  
+    "src_dej": 'd',  
+    "tstart": 'd',  
+    "tsamp": 'd',  
+    "nbits": 'i', 
+    "nsamples": 'i', 
+    "nbeams": "i",
+    "ibeam": "i",
+    "fch1": 'd',  
+    "foff": 'd',
+    "FREQUENCY_START": 'flag',
+    "fchannel": 'd',  
+    "FREQUENCY_END": 'flag',
+    "nchans": 'i', 
+    "nifs": 'i', 
+    "refdm": 'd',  
+    "period": 'd',  
+    "npuls": 'q',
+    "nbins": 'i', 
+    "HEADER_END": 'flag'}
+
+def dec2radians(src_dej):
+    """
+    dec2radians(src_dej):
+       Convert the SIGPROC-style DDMMSS.SSSS declination to radians
+    """
+    sign = 1.0
+    if (src_dej < 0): sign = -1.0;
+    xx = math.fabs(src_dej)
+    dd = int(math.floor(xx / 10000.0))
+    mm = int(math.floor((xx - dd * 10000.0) / 100.0))
+    ss = xx - dd * 10000.0 - mm * 100.0
+    return sign * ARCSECTORAD * (60.0 * (60.0 * dd + mm) + ss)
+
+def ra2radians(src_raj):
+    """
+    ra2radians(src_raj):
+       Convert the SIGPROC-style HHMMSS.SSSS right ascension to radians
+    """
+    return 15.0 * dec2radians(src_raj)
+
+def read_doubleval(filfile, stdout=False):
+    dblval = struct.unpack('d', filfile.read(8))[0]
+    if stdout:
+        print("  double value = '%20.15f'"%dblval)
+    return dblval
+
+def read_intval(filfile, stdout=False):
+    intval = struct.unpack('i', filfile.read(4))[0]
+    if stdout:
+        print("  int value = '%d'"%intval)
+    return intval
+
+def read_longintval(filfile, stdout=False):
+    longintval = struct.unpack('q', filfile.read(8))[0]
+    if stdout:
+        print("  long int value = '%d'"%longintval)
+    return longintval
+
+def read_string(filfile, stdout=False):
+    strlen = struct.unpack('i', filfile.read(4))[0]
+    strval = filfile.read(strlen)
+    if stdout:
+        print("  string = '%s'"%strval)
+    return strval
+
+def read_paramname(filfile, stdout=False):
+    paramname = read_string(filfile, stdout=False)
+    if stdout:
+        print("Read '%s'"%paramname)
+    return paramname
+
+def read_hdr_val(filfile, stdout=False):
+    paramname = read_paramname(filfile, stdout)
+    if header_params[paramname] == 'd':
+        return paramname, read_doubleval(filfile, stdout)
+    elif header_params[paramname] == 'i':
+        return paramname, read_intval(filfile, stdout)
+    elif header_params[paramname] == 'q':
+        return paramname, read_longintval(filfile, stdout)
+    elif header_params[paramname] == 'str':
+        return paramname, read_string(filfile, stdout)
+    elif header_params[paramname] == 'flag':
+        return paramname, None
+    else:
+        print("Warning:  key '%s' is unknown!" % paramname)
+        return None, None
+
+def prep_string(string):
+    return struct.pack('i', len(string))+string
+
+def prep_double(name, value):
+    return prep_string(name)+struct.pack('d', float(value))
+
+def prep_int(name, value):
+    return prep_string(name)+struct.pack('i', int(value))
+
+def addto_hdr(paramname, value):
+    if header_params[paramname] == 'd':
+        return prep_double(paramname, value)
+    elif header_params[paramname] == 'i':
+        return prep_int(paramname, value)
+    elif header_params[paramname] == 'str':
+        return prep_string(paramname) + prep_string(value)
+    elif header_params[paramname] == 'flag':
+        return prep_string(paramname)
+    else:
+        warnings.warning("key '%s' is unknown!" % paramname)
+    return hdr
+
+def read_header(infile):
+    """
+    read_header(infile):
+       Read a SIGPROC-style header and return the keys/values in a dictionary,
+          as well as the length of the header: (hdrdict, hdrlen)
+    """
+    hdrdict = {}
+    if type(infile) == type("abc"):
+        infile = open(infile)
+    param = ""
+    while (param != "HEADER_END"):
+        param, val = read_hdr_val(infile, stdout=False)
+        hdrdict[param] = val
+    hdrlen = infile.tell()
+    infile.close()
+    return hdrdict, hdrlen
+
+def samples_per_file(infile, hdrdict, hdrlen):
+    """
+    samples_per_file(infile, hdrdict, hdrlen):
+       Given an input SIGPROC-style filterbank file and a header
+           dictionary and length (as returned by read_header()),
+           return the number of (time-domain) samples in the file.
+    """
+    numbytes = os.stat(infile)[6] - hdrlen
+    bytes_per_sample = hdrdict['nchans'] * (hdrdict['nbits']/8)
+    if numbytes % bytes_per_sample:
+        print("Warning!:  File does not appear to be of the correct length!")
+    numsamples = numbytes / bytes_per_sample
+    return numsamples
+
+if __name__ == "__main__":
+    if len(sys.argv)==1:
+        print("\nusage:  mod_filterbank_hdr.py infile.fil [outfile.fil]\n")
+        sys.exit()
+    filhdr = {}
+    newhdr = ""
+    infile = open(sys.argv[1], 'rb')
+
+    # Loop over the values in the .fil file
+    while 1:
+        param, val = read_hdr_val(infile, stdout=False)
+        filhdr[param] = val
+
+        # Add lines here to correct stuff
+        #if param=="nchans":  val = 768
+
+        # Append to the new hdr string
+        newhdr += addto_hdr(param, val)
+
+        # Break out of the loop if the header is over
+        if param=="HEADER_END":  break
+            
+    if len(sys.argv) > 2:
+        print("Writing new header to '%s'"%sys.argv[2])
+        outfile = open(sys.argv[2], 'wb')
+        outfile.write(newhdr)
+        outfile.close()
+    else:
+        print(filhdr)

tofil.py

+#!/usr/bin/env python2
+from __future__ import print_function
+
+import numpy as np
+import filterbank
+import sigproc
+from astropy.time import Time
+
+import os.path
+import sys
+import re
+
+import argparse
+
+def create_filterbank(infile, outfile=None, date=None, source=None):
+    """
+    Create a filterbank file from the CAMRAS-format, which is generated
+    by pulsar_record.
+
+    infile:   name of the file generated by pulsar_record
+    outfile: name of the filterbank file to be created; if not given, a name
+              will be created from the infile
+    date:     the date to be used; default is the file creation date of infile
+              if given, the date must be in a format the astropy.time understands
+    source:   name of the source to be used; default is guessed from filename
+    """
+
+    if not outfile:
+        outfile = re.split('[+-\.]', infile)[0] + ".fil"
+
+    if not source:
+        match = re.match("[BJ][0-9]*[+-][0-9]*", infile)
+        if match:
+            source = "PSR " + match.group(0)
+        else:
+            raise ValueError("Could not guess source name from filename")
+
+    if not date:
+        date = Time(os.path.getctime(infile), format='unix')
+    else:
+        date = Time(date, format='isot')
+
+    print("Creating filterbank file {} from {}".format(outfile, infile))
+    print("Source name: {}".format(source))
+    print("Observation time: {}".format(date))
+
+    date = date.mjd
+ 
+    # Generate filterbank header
+    fil_header = {}
+    fil_header["telescope_id"]  = sigproc.telescope_ids["Effelsberg"]
+    fil_header["machine_id"]    = sigproc.machine_ids["FAKE"]
+    fil_header["data_type"]     = 1
+    fil_header["rawdatafile"]   = infile
+    fil_header["source_name"]   = source
+    fil_header["barycentric"]   = 0
+    fil_header["pulsarcentric"] = 0
+    fil_header["az_start"]      = 0.0
+    fil_header["za_start"]      = 0.0
+    fil_header["src_raj"]       = 0.0
+    fil_header["src_dej"]       = 0.0
+    fil_header["tstart"]        = date
+    fil_header["tsamp"]         = 0.00046811428571414528 # 1.0/35e6*256 chans*64 decimation
+    fil_header["nbits"]         = 32
+    fil_header["fch1"]          = 441.4
+    fil_header["foff"]          = -35.0/256.0
+    fil_header["nchans"]        = 256
+    fil_header["nifs"]          = 1
+
+    # Write header
+    out = filterbank.create_filterbank_file(outfile, fil_header, nbits=32)
+
+    # Read file
+    data = np.flipud(np.fromfile(infile, dtype='>u4'))
+    print(data.shape)
+
+    data = np.flipud(data.reshape(-1, 256))
+
+    # Check that no packets were dropped
+    if np.sum(data[:-1,-1]-data[1:,-1]-1) > 0.1:
+        print("Packets were dropped:", np.count_nonzero(data[:-1,0]-data[1:,0]-1))
+
+    # Write data
+    out.append_spectra(data)
+
+    # Close file