Data format
Baseband
EPN
PSRFITS
The “PSRFITS” format supports mean pulse profile (‘fold-mode’) and streamed (‘search-mode’) multi-channel full-polarisation data. A PSRFITS format file consists of a primary header-data unit (HDU) followed by a series of binary extension HDUs, in which are stored specific information about an observation.
Usually, there are three types of pulsar data products are available:
Fold-mode observations
Fold-mode observations are recorded at the telescope for a particular known pulsar, where the data are stacked or ‘folded’ at the rotation period of the pulsar, to form a pulse profile averaged over the length of the observation.
Files of this type have the extension ‘.rf’. All fold-mode files are also processed by averaging over frequency channels, polarisation and time - files of this type have the extension ‘.FTp’.
Search-mode observations
Search-mode observations are essentially comprised of a multi-channel full-polarisation data stream for the length of the observation.
Files of this type have the extension ‘.sf’.
Calibration files
Before and after an observation, a signal from a linear noise diode is injected into the feed. This allows the associated pulsar observation to be polarimetrically calibrated.
Files of this type have the extension ‘.cf’.
Psrchive
Checking that it’s working:
import psrchive as pr
Fundamental data classes:
Profile is a single pulse profile – data as a function of pulse phase only.
Integration is a set of pulse profiles recorded simultaneously – usually
profiles as a function of frequency channel and/or polarization.
Archive is a set of Integration as a function of time. Represents a
single data file.
 Fig 1, Profile
|
 Fig 2, Integration
|
 Fig 3, Archive
|
Accessing data with these classes
Archive:
An Archive object represents a single data file, and contains a collection of Integration objects, also known as subintegrations.
arch = pr.Archive_load('file.fits')
arch.get_source()
# '1744-1134'
The get_data() method is an extra Python function that returns the raw data from the file as a NumPy array:
data = arch.get_data()
data.shape
# (10, 4, 2048, 256)
The dimensions of the data array are subintegration, polarization, channel, and profile bin. Important note: The get_data() function returns a copy of the Archive data. So changes made to the returned array will not be reflected in the Archive, and vice-versa.
Integrations within an Archive can be accessed in Python by indexing the Archive object use archive.get_Integration(isub) or archive[isub] methods. So the following two lines are equivalent:
subint = arch.get_Integration(0)
#等价
subint = arch[0]
Integration:
An Integration object represents a collection of pulse profiles recorded simultaneously, typically from a number of frequency channels and polarization states (Stokes parameters):
subint.get_duration()
# 60.230205439999999
subint.get_nchan()
# 2048
The existing baseline_stats() and cal_levels() methods have been modified for Python to return tuples of arrays, for example:
(b_mean, b_var) = subint.baseline_stats()
b_mean.shape
# (4, 2048)
An Integration object contains a number of pulse profiles that can be accessed via integration.get_Profile(ipol,ichan):
prof = subint.get_Profile(0,0)
Profile:
A Profile object represents a single pulse profile. The profile.get_amps() method returns the profile data as a NumPy array:
pdata = prof.get_amps()
pdata.shape
# (256,)
In contrast to the arch.get_data() example above, the pdata array here points into the actual C++ data structure. So changes made to the array values will be reflected in the underlying Profile object.
The profile data values can also be read and/or altered by indexing the Profile object directly, for example:
pdata[16]
# 91.366058
prof[16]
#91.366058349609375
prof[16] = 500.0
pdata[16]
# 500.0
Shortcut to get all data:
Use archive.get_data() to return entire (Nsub, Npol, Nchan, Nbin) data cube as a NumPy array.
Examples:
import psrchive as pr
import matplotlib.pyplot as plt
arch = psrchive.Archive_load('file.fits')
data = arch[0].get_Profile(0,200).get_amps()
print(data.shape)
# (1024,)
plt.plot(data)
arch = psrchive.Archive_load(‘file.fits’)
data = arch.get_data()
print(data.shape)
(64,4.25,1024)
plt.plot(data[0,0,200:],’g’)
<div align="center"><img src="https://github.com/AstronDog/blog-pic/raw/master/20200110psrpy/profile_example2.png" alt="" width="550"></div>
#### Important: view and copy
One data access subtlety / gotcha:
**profile.get_amps()** returns a **view** of the original data
**archive.get_data()** returns a **copy** of the original data
This means that if you want to modify the data in a file you need to change the values in the results of profile.get_amps()
Modified data files can be saved to disk using
**archive.unload(“new_filename”)**
#### Plot profile
Archive has a large number of methods (functions) for performing common data processing steps.
Common examples: dedisperse(), remove_baseline(), fscrunch(), tscrunch(), pscrunch(), convert_state(), …
archive.execute(“[psrsh code…]”) will run any psrsh command on the archive.
import psrchive as pr
import matplotlib.pyplot as plt
arch = psrchive.Archive_load(‘file.fits’)
arch.dedisperse()
arch.fscrunch()
arch.pscrunch()
arch.remove_baseline()
data = a[0].get_Profiel(0,0).get_amps()
plt.plot(data)
<div align="center"><img src="https://github.com/AstronDog/blog-pic/raw/master/20200110psrpy/preprocess.png" alt="" width="550"></div>
Looping over profiles, extracting data:
import psrchive as pr
arch = psrchive.Archive_load(‘file.fits’)
for isub in range(a.get_nsubint())
i = a[isub]
for ichan in range(a.get_nchan()):
print isub,ichan,i.get_Profiel(0,ichan).get_amps()[0]
#### MJD
The MJD class is also now available in Python to deal with high-precision date/time values. Return values for date/time functions such as **Integration.get_epoch()** were previously converted to doubles for convenience. To reproduce the old behavior and get a double-precision MJD value, use the **MJD.in_days()** method, for example:
mjd_dbl = arch.get_Integration(0).get_start_time().in_days()
#### Generating TOAs with Arrival Time
The ArrivalTime class can be used to generate pulse times of arrival (TOAs) directly in Python, as an alternative to using the command-line utility pat. Here is an example illustrating the basic usage:
import psrchive
arrtim = psrchive.ArrivalTime()
arrtim.set_shift_estimator(‘PGS’) # Set algorithm (see ‘pat -A’ help)
arrtim.set_format(‘Tempo2’) # set TOA format
arrtim.set_format_flags(‘IPTA’) # set some TOA flags
arrtim.set_attributes([‘chan’,’subint’]) # More TOA flags
Load template profile
std = psrchive.Archive_load(‘J1713+0747.Rcvr1_2.GUPPI.9y.x.sum.sm’)
std.pscrunch()
arrtim.set_standard(std)
Load observation profiles
obs = psrchive.Archive_load(‘guppi_55616_J1713+0747_0009.12y.x.ff’)
obs.pscrunch()
arrtim.set_observation(obs)
Result is a tuple of TOA strings:
toas = arrtim.get_toas()
#### Simple Example
In this example, we recreate the helpful "pav -GTdp" plot using python, and demonstrate some common psrchive functionality:
#! /usr/bin/env python
import pylab
import psrchive
arch = psrchive.Archive_load(‘1744_0001_0001.fits’)
arch.bscrunch_to_nbin(256)
arch.dedisperse()
arch.fscrunch_to_nchan(512)
arch.remove_baseline()
arch.convert_state(‘Stokes’)
data = arch.get_data()
freq_lo = arch.get_centre_frequency() - arch.get_bandwidth()/2.0
freq_hi = arch.get_centre_frequency() + arch.get_bandwidth()/2.0
pylab.imshow(data[:,0,:,:].mean(0),extent=(0,1,freq_lo,freq_hi))
pylab.xlabel(‘Pulse phase’)
pylab.ylabel(‘Frequency (MHz)’)
pylab.savefig(‘python_plot_example.png’)
This results in the following plot:
<div align="center"><img src="https://github.com/AstronDog/blog-pic/raw/master/20200110psrpy/frequency.png" alt="" width="550"></div>
#### 3 Timing
Most use cases only need Archive, Integration, and Profile classes.
But some of the PSRCHIVE algorithm classes are also included in the Python interface.
For example, ProfileShiftFit for doing template-matching:
No very comprehesive list of these unfortunately. Browse the C++ class docs, let me know if you want something added.
import psrchive as pr
import matplotlib.pyplot as plt
a = psrchive.Archive_load(‘file.fits’).total()
s = psrchive.Archive_load(‘file.std’)
aprof = a[0].get_Profile(0,0)
sprof = s[0].get_Profile(0,0)
psf = pr.ProfileShiftFit()
psf.set_standard(sprof)
psr.set_Profile(aprof)
psr.get_shifts()
subplot(211); plot(aprof.get_a,ps())
subplot(212); plot(sprof.get_a,ps())
<div align="center"><img src="https://github.com/AstronDog/blog-pic/raw/master/20200110psrpy/advance.png" alt="" width="550"></div>
### Archive Functions
#### Public Member Functions
##### Archive
void **bscrunch** (unsigned nscrunch)
Integrate pulse profiles in phase. More...
void **pscrunch** ()
Integrate profiles in polarization. More...
void **fscrunch** (unsigned nscrunch=0)
Integrate profiles in frequency. More...
void **tscrunch** (unsigned nscrunch=0)
Integrate profiles in time. More...
void **remove_chan** (unsigned first, unsigned last)
Delete the specified inclusive channel range from the Archive.
#### Signal::Basis
void **get_basis** () const
Convenience interface to Receiver::get_basis.
void **dedisperse** ()
Rotate the Profiles to remove dispersion delays b/w chans. More...
void defaraday ()
Correct the Faraday rotation of Q into U. More...
#### PhaseWeight
void **baseline** () const
Return a new PhaseWeight instance with the baseline phase bins masked.
void **remove_baseline** ()
Remove the baseline from all profiles.
void **uniform_weight** (float new_weight=1.0)
Set the weight of each profile to the given number.
void **bscrunch_to_nbin** (unsigned new_nbin)
Call bscrunch with the appropriate value. More...
void **fscrunch_to_nchan** (unsigned new_nchan)
Call fscrunch with the appropriate value. More...
void **tscrunch_to_nsub** (unsigned new_nsub)
Call tscrunch with the appropriate value. More...
MJD **start_time** () const
Return the MJD at the start of the first sub-integration.
MJD **end_time** () const
Return the MJD at the end of the last sub-integration.
double **integration_length** () const
Returns the total time integrated into all Integrations (in seconds)
### Common Attributes
These pure virtual methods provide access to the common attributes stored by all derived classes.
virtual std::string **get_telescope** () const =0
Get the name of the telescope used.
virtual void **set_telescope** (const std::string &code)=0
Set the name of the telescope used.
virtual Signal::State **get_state** () const =0
Get the state of the profile data.
virtual void **set_state** (Signal::State state)=0
Set the state of the profile data.
virtual Signal::Scale **get_scale** () const =0
Get the scale in which flux density is measured.
virtual void **set_scale** (Signal::Scale scale)=0
Set the scale in which flux density is measured.
virtual Signal::Source **get_type** () const =0
Get the observation type (psr, cal)
virtual void **set_type** (Signal::Source type)=0
Set the observation type (psr, cal)
virtual std::string **get_source** () const =0
Get the source name.
virtual void **set_source** (const std::string &source)=0
Set the source name.
virtual sky_coord **get_coordinates** () const =0
Get the coordinates of the source.
virtual void **set_coordinates** (const sky_coord &coordinates)=0
Set the coordinates of the source.
virtual double **get_centre_frequency** () const =0
Get the centre frequency of the observation.
virtual void **set_centre_frequency** (double cf)=0
Set the centre frequency of the observation.
virtual double **get_bandwidth** () const =0
Get the overall bandwidth of the observation.
virtual void **set_bandwidth** (double bw)=0
Set the overall bandwidth of the observation.
virtual double **get_dispersion_measure** () const =0
Get the dispersion measure (in ${\rm pc\, cm}^{-3}$)
virtual void **set_dispersion_measure** (double dm)=0
Set the dispersion measure (in ${\rm pc\, cm}^{-3}$)
virtual bool **get_dedispersed** () const =0
Inter-channel dispersion delay has been removed.
virtual void **set_dedispersed** (bool done=true)=0
Set the value to be returned by get_dedispersed.
virtual double **get_rotation_measure** () const =0
Get the rotation measure (in ${\rm rad\, m}^{-2}$)
virtual void **set_rotation_measure** (double rm)=0
Set the rotation measure (in ${\rm rad\, m}^{-2}$)
virtual bool **get_faraday_corrected** () const =0
Data has been corrected for ISM faraday rotation.
virtual void **set_faraday_corrected** (bool done=true)=0
Set the value to be returned by get_ism_rm_corrected.
virtual bool **get_poln_calibrated** () const =0
Data has been calibrated for polarimetric response of instrument.
virtual void **set_poln_calibrated** (bool done=true)=0
Set the value to be returned by get_poln_calibrated.
#### Dimension Attributes
These pure virtual methods are used by the Archive class to set the dimension attributes stored by the derived classes.
virtual unsigned **get_nbin** () const =0
Get the number of pulsar phase bins used.
virtual unsigned **get_nchan** () const =0
Get the number of frequency channels used.
virtual unsigned **get_npol** () const =0
Get the number of polarizations.
virtual void **set_nbin** (unsigned numbins)=0
Set the number of pulsar phase bins.
virtual void **set_nchan** (unsigned numchan)=0
Set the number of frequency channels.
virtual void **set_npol** (unsigned numpol)=0
Set the number of polarization measurements.
virtual void **resize** (unsigned nsubint, unsigned npol=0, unsigned nchan=0, unsigned nbin=0)
Resize the Integration vector with new_Integration instances. More...
virtual void **erase** (unsigned isubint)
Remove the specified sub-integration.
## Astropy
## Pypulse
## pdat
## Pint
## References
[ANTF](https://www.atnf.csiro.au/observers/data/ppdu_guide.html)
