This tutorial covers some of the basic functionality of the sigpyproc package. For a guide on how to extend the package, see the

Getting started

For test puproses a small 2-bit filterbank file is included in the /examples/ directory of the sigpyproc package.

Loading data into sigpyproc

Lets start by loading our filterbank file into sigpyproc. To do this, we require the :class:FilReader class from the :automodule:sigpyproc.Readers module.

[1]:

import numpy as np
from sigpyproc.Readers import FilReader

[2]:

myFil = FilReader("tutorial.fil")

[3]:

myFil

[3]:

<sigpyproc.Readers.FilReader at 0x7f29d0537ee0>

myFil now contains an instance of the :class:sigpyproc.Readers.FilReader class. We can access obervational meta-data through the myFil.header attribute:

[4]:

myFil.header

[4]:

{'telescope_id': 0,
 'nbits': 2,
 'fch1': 1510.0,
 'data_type': 1,
 'nchans': 64,
 'tsamp': 0.00032,
 'foff': -1.09,
 'tstart': 50000.0,
 'source_name': 'P: 250.000000000000 ms, DM: 30.000',
 'nifs': 1,
 'machine_id': 0,
 'hdrlen': 244,
 'filelen': 3000564,
 'nbytes': 3000320,
 'nsamples': 187520,
 'filename': 'tutorial.fil',
 'basename': 'tutorial',
 'extension': '.fil',
 'bandwidth': 69.76,
 'ftop': 1510.545,
 'fbottom': 1440.785,
 'fcenter': 1475.665,
 'tobs': 60.006400000000006,
 'src_raj': 0,
 'src_dej': 0,
 'ra': '00:00:00.0000',
 'dec': '00:00:00.0000',
 'ra_rad': 0.0,
 'dec_rad': 0.0,
 'ra_deg': 0.0,
 'dec_deg': 0.0,
 'obs_date': '09/10/1995',
 'obs_time': '00:00:00.00000',
 'dtype': '<u1'}

where

[5]:

type(myFil.header)

[5]:

sigpyproc.Header.Header

All values stored in the myFil.header attribute may be accessed both as dictionary items and/or as attributes, i.e.:

[6]:

myFil.header["nchans"]

[6]:

64

[7]:

myFil.header.nchans

[7]:

64

Now that we know how to load a file into sigpyproc, let’s look at at doing something with the loaded data.

Dedispersing the data

One of the most used techniques in pulsar processing is dedispersion, wherein we add or remove frequency dependent time delays to the data.

To dedisperse our myFil instance, we simply call the dedisperse method:

[8]:

myTim = myFil.dedisperse(30)

dedisperse : 100%|██████████| 19/19 [00:00<00:00, 339.45it/s]


Filterbank reading plan:
------------------------
Called on file:       tutorial.fil
Called by:            dedisperse
Number of samps:      187520
Number of reads:      18
Nsamps per read:      10000
Nsamps of final read: 7826
Nsamps to skip back:  17




[9]:

myTim

[9]:

TimeSeries([108., 100., 102., ..., 105., 111., 107.], dtype=float32)

[10]:

type(myTim)

[10]:

sigpyproc.TimeSeries.TimeSeries

Here we have dedispersed to a DM of 30 pc cm\(^{-3}\) with the result being an instance of the :automodule:sigpyproc.TimeSeries.TimeSeries class, which we have called myTim.

The :automodule:sigpyproc.TimeSeries.TimeSeries class in a subclass of numpy.ndarray, and is capable of using all standard numpy functions. For example:

[11]:

myTim.sum()

[11]:

TimeSeries(19636992., dtype=float32)

[12]:

myTim.max()

[12]:

TimeSeries(121., dtype=float32)

[13]:

myTim.min()

[13]:

TimeSeries(88., dtype=float32)

[14]:

np.median(myTim)

[14]:

TimeSeries(105., dtype=float32)

The use of numpy.ndarray subclasses is important in allowing sigpyproc to easily interface with many 3rd party python libraries.

Performing a Fourier transform

To perform a discrete fourier transform of the data contained in the myTim instance we may invoke the myTim.rFFT method.

[15]:

myFS = myTim.rFFT()

[16]:

type(myFS)

[16]:

sigpyproc.FourierSeries.FourierSeries

[17]:

myFS

[17]:

FourierSeries([ 1.9636880e+07,  0.0000000e+00, -2.9421130e+02, ...,
               -1.7792952e+03, -2.7670000e+03,  0.0000000e+00],
              dtype=float32)

The :automodule:sigpyproc.FourierSeries.FourierSeries is also a subclass of numpy.ndarray, where array elements are [real, imaginary, real, imaginary, real, imaginary, real,…].

Using the rednoise method of myFS, we can de-redden the Fourier series:

[18]:

myFS_red = myFS.rednoise()

[19]:

myFS_red

[19]:

FourierSeries([ 1.        ,  0.        , -0.24654439, ..., -1.641717  ,
               -2.5530508 ,  0.        ], dtype=float32)

with the dereddened fourier series, we can now form the power spectrum of the observation:

[20]:

mySpec = myFS_red.formSpec(interpolated=True)

[21]:

mySpec

[21]:

PowerSpectrum([1.        , 0.91751015, 0.78774583, ..., 2.0073905 ,
               2.2536902 , 2.5530508 ], dtype=float32)

Here we have set the interpolated flag to True, causing the formSpec function to perform nearest bin interpolation.

mySpec contains several convenience methods to help with navigating the power spectrum. For instance:

[22]:

mySpec.period2bin(0.25)

[22]:

240

[23]:

mySpec.freq2bin(5.0)

[23]:

300

We can also perofrm Lyne-Ashworth harmonic folding to an arbitrary number of harmonics:

[24]:

folds = mySpec.harmonicFold(5)

[25]:

folds

[25]:

[PowerSpectrum([1.       , 1.8350203, 1.705256 , ..., 3.460951 , 3.7072508,
                2.5530508], dtype=float32),
 PowerSpectrum([1.       , 3.7525306, 3.622766 , ..., 3.460951 , 3.7072508,
                2.5530508], dtype=float32),
 PowerSpectrum([1.       , 7.5875506, 7.457786 , ..., 3.460951 , 3.7072508,
                2.5530508], dtype=float32),
 PowerSpectrum([ 1.       , 15.257591 , 15.127827 , ...,  3.460951 ,
                 3.7072508,  2.5530508], dtype=float32),
 PowerSpectrum([ 1.       , 30.59768  , 30.467915 , ...,  3.460951 ,
                 3.7072508,  2.5530508], dtype=float32)]

Where the variable folds is a python list containing each of the requested harmonic folds.

Folding data

Both the :automodule:sigpyproc.TimeSeries.TimeSeries and the :automodule:sigpyproc.FourierSeries.FourierSeries have methods to phase fold their data. Using our earlier myFil instance, we will fold our filterbank file with a period of 250 ms and a DM of pc cm\(^{-3}\) and acceleration of 0 ms\(^{-2}\).

[26]:

myFold = myFil.fold(0.25,30.,accel=0,nbins=128,nints=32,nbands=16)

fold : 100%|██████████| 19/19 [00:00<00:00, 288.27it/s]


Filterbank reading plan:
------------------------
Called on file:       tutorial.fil
Called by:            fold
Number of samps:      187520
Number of reads:      18
Nsamps per read:      10000
Nsamps of final read: 7826
Nsamps to skip back:  17




[27]:

type(myFold)

[27]:

sigpyproc.FoldedData.FoldedData

[28]:

myFold.shape

[28]:

(32, 16, 128)

The the :automodule:sigpyproc.FoldedData.FoldedData has several functions to enable simple slicing and summing of the folded data cube. These include:

  • getSubband: select all data in a single frequency band

  • getSubint: select all data in a single subintegration

  • getFreqPhase: sum the data in the time axis

  • getTimePhase: sum the data in the frequency axis *getProfile: get the pulse profile of the fold

We can also tweak the DM and period of the fold using the updateParams method:

[29]:

myFold.updateParams(period=0.2502,dm=100)

Tips and tricks

There are several tips and tricks to help speed up sigpyproc and also make it more user friendly. For people who are familiar with Python and IPython these will be old news, but for newbies these may be of use.

Tab completion: One of the many nice things about IPython is that it allows for tab completion:

[30]:

myFil   # then press tab

[30]:

<sigpyproc.Readers.FilReader at 0x7f29d0537ee0>

Docstrings: by using question marks or double question marks we can access both information about a function and its raw source:

[31]:

myFil.downsample?

Signature:
myFil.downsample(
    tfactor=1,
    ffactor=1,
    gulp=512,
    filename=None,
    back_compatible=True,
    **kwargs,
)
Docstring:
Downsample data in time and/or frequency and write to file.

:param tfactor: factor by which to downsample in time
:type tfactor: int

:param ffactor: factor by which to downsample in frequency
:type ffactor: int

:param gulp: number of samples in each read
:type gulp: int

:param filename: name of file to write to (defaults to ``basename_tfactor_ffactor.fil``)
:type filename: str

:param back_compatible: sigproc compatibility flag (legacy code)
:type back_compatible: bool

:return: output file name
:rtype: :func:`str`
File:      ~/anaconda3/lib/python3.8/site-packages/sigpyproc/Filterbank.py
Type:      method

Note that all docstrings are written in reStructuredText. This is to facilitate automatic documentation creation with the Sphinx package.

Threading: The use of OpenMP means that several of C++ library calls in sigpyproc can be sped up:

[32]:

myFil.setNthreads(2) # enable OpenMP to use 2 threads

It should be noted that threading will not always speed up a process, and that the performance does not scale linearly with number of threads used.

Chaining: The ability to chain together methods, combined with history recall in IPython means that it is simple to condense a sigpyproc request into a single line:

[33]:

spectrum = FilReader("tutorial.fil").collapse().rFFT().rednoise().formSpec(True)

collapse : 100%|██████████| 367/367 [00:00<00:00, 8174.90it/s]


Filterbank reading plan:
------------------------
Called on file:       tutorial.fil
Called by:            collapse
Number of samps:      187520
Number of reads:      366
Nsamps per read:      512
Nsamps of final read: 128
Nsamps to skip back:  0

Here we create a FilReader instance, which is then collapsed in frequency, FFT’d, cleaned of rednoise and interpolated to form a power spectrum. In the intrests of readability, this is not always a good idea, however for testing code quickly, it is invaluable.