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, few small size filterbank files are included in the /tests/data/ 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:~sigpyproc.readers.FilReader class from the :automodule:sigpyproc.readers module.

import numpy as np
from rich.pretty import Pretty
from sigpyproc.readers import FilReader

myFil = FilReader("../../tests/data/tutorial.fil")

myFil

<sigpyproc.readers.FilReader at 0x7f4cf19b87c0>

myFil now contains an instance of the :class:sigpyproc.readers.FilReader class. We can access obervational metadata through the myFil.header attribute:

Pretty(myFil.header)

Header(
    filename='../../tests/data/tutorial.fil',
    data_type='filterbank',
    nchans=64,
    foff=-1.09,
    fch1=1510.0,
    nbits=2,
    tsamp=0.00032,
    tstart=50000.0,
    nsamples=187520,
    nifs=1,
    coord=<SkyCoord (ICRS): (ra, dec) in deg
    (0., 0.)>,
    azimuth=<Angle 0. deg>,
    zenith=<Angle 0. deg>,
    telescope='Fake',
    backend='FAKE',
    source='P: 250.000000000000 ms, DM: 30.000',
    frame='topocentric',
    ibeam=0,
    nbeams=0,
    dm=0,
    period=0,
    accel=0,
    signed=False,
    rawdatafile='',
    hdrlens=[244],
    datalens=[3000320],
    filenames=['../../tests/data/tutorial.fil'],
    nsamples_files=[187520],
    tstart_files=[50000.0]
)

where

type(myFil.header)

sigpyproc.header.Header

All values stored in the myFil.header attribute may be accessed as attributes, i.e.:

myFil.header.nchans

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:

myTim = myFil.dedisperse(30)

/home/docs/checkouts/readthedocs.org/user_builds/sigpyproc3/envs/stable/lib/python3.9/site-packages/rich/live.py:23
1: UserWarning: install "ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

myTim

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

type(myTim)

sigpyproc.timeseries.TimeSeries

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

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

myTim.sum()

TimeSeries(19636992., dtype=float32)

myTim.max()

TimeSeries(121., dtype=float32)

myTim.min()

TimeSeries(88., dtype=float32)

np.median(myTim)

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.

myFS = myTim.rfft()

type(myFS)

sigpyproc.fourierseries.FourierSeries

myFS

FourierSeries([ 1.9636884e+07   +0.j     , -2.9424850e+02 +429.87863j,
                5.4652838e+02 -577.44696j, ...,
               -1.3942198e+03+1670.1677j , -1.2117781e+03-1779.332j  ,
               -2.7670000e+03   +0.j     ], dtype=complex64)

The :autoclass:sigpyproc.fourierseries.FourierSeries is also a subclass of numpy.ndarray, where array elements are numpy.complex64.

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

myFS_red = myFS.remove_rednoise()

myFS_red

FourierSeries([ 1.0000000e+00+0.0000000e+00j,
               -1.9937569e+20+2.9127539e+20j,
                3.7031446e+20-3.9126415e+20j, ...,
               -8.7243772e-01+1.0451130e+00j,
               -7.5827420e-01-1.1134230e+00j,
               -1.7314595e+00+0.0000000e+00j], dtype=complex64)

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

mySpec = myFS_red.form_spec(interpolated=True)

mySpec

PowerSpectrum([1.       ,       inf,       inf, ..., 1.3613995, 1.5284487,
               1.7314595], 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:

mySpec.period2bin(0.25)

mySpec.freq2bin(5.0)

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

folds = mySpec.harmonic_fold(5)

folds

[PowerSpectrum([1.       ,       inf,       inf, ..., 2.7979643, 2.9650135,
                1.7314595], dtype=float32),
 PowerSpectrum([1.       ,       inf,       inf, ..., 2.7979643, 2.9650135,
                1.7314595], dtype=float32),
 PowerSpectrum([1.       ,       inf,       inf, ..., 2.7979643, 2.9650135,
                1.7314595], dtype=float32),
 PowerSpectrum([1.       ,       inf,       inf, ..., 2.7979643, 2.9650135,
                1.7314595], dtype=float32),
 PowerSpectrum([1.       ,       inf,       inf, ..., 2.7979643, 2.9650135,
                1.7314595], dtype=float32)]

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

Folding data#

Both the :autoclass:sigpyproc.timeseries.TimeSeries and the :autoclass: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}$.

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

type(myFold)

sigpyproc.foldedcube.FoldedData

myFold.shape

(32, 16, 128)

The the :autoclass:sigpyproc.foldedcube.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:

myFold.update_dm(dm=100)
myFold.update_period(period=0.2502)

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:

myFil   # then press tab

<sigpyproc.readers.FilReader at 0x7f4cf19b87c0>

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

myFil.downsample?

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

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:

spectrum = FilReader("../../tests/data/tutorial.fil").collapse().rfft().remove_rednoise().form_spec(True)

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.