Epoch Folding Basics

The folding algorithm makes use of the fact that the pulse from a pulsar has a period which is extremely regular.  By sampling the pulsar signal at a stable, known rate it is possible to 'synchronously integrate' the pulsar signal in software.  This can be done in a limited fashion in real-time, but is more conveniently done on saved files which are analysed after the observation session - i.e., "off-air". 

We will assume that the RF signal (say @ 400 MHz) has been down-converted to a lower frequency (e.g. - 10 MHz) at, say, 1 MHz bandwidth.  This signal is fed to total power detector (true RMS) and the output is sampled and recorded.

Imaginary One Second Pulsar

Imagine we are observing a pulsar whose pulse period is exactly 1 second (1000 ms) - a fairly typical value.  As most pulsars have a duty-cycle of about 5%, our imaginary typical pulsar has a pulse that lasts for 50 ms.

Lowpass Filtering of Total Power Detector Output

For a 1 second pulsar with 50 ms pulse width, the output of the power detector can be lowpass-filtered with a corner frequency about 100 Hz without affecting the pulse.  Note that pulsars with pulse periods < 1 second will need a proportionally higher cut-off frequency for the lowpass filter to ensure the pulse shape is preserved.

Sampling Frequency

As the output signal has been lowpass-filtered with a corner frequency of 100 Hz, a sampling frequency at the Nyquist frequency or higher is required (e.g. - > 200 sps).  For convenience (governed by the sampling hardware) the output can be sampled at a higher rate - but this produces larger files for a given observation time.  A higher corner frequency can be used for the lowpass filter if desired - but the sampling frequency needs to be raised to the Nyquist frequency or higher (related to the lowpass cutoff frequency - not the pulsar pulse frequency) to avoid aliasing, which degrades S/N.

The Concept of a 'Record'

In the above figure of our simulated 1 second, 5% duty-cycle pulsar, we can see that that part of the sampled data which covers the time between consecutive pulses can be viewed as a single 'record' of a pulse.  In the above figure the data contains 6 pulses, and so we have 6 independent 'records' of the pulse.

Referring back to the general radiometer equation ("Radio Astronomy" - Kraus)...

where

ΔSmin = minimum detectable flux density (watts per square metre per hertz)

k = Boltzmann's constant (1.38064852 × 10-23 Joules • K-1)

Ks = factor for a particular observatory system, hopefully near to 1

Tsys = System noise temperature (Ko)

Ae = Effective aperture of the antenna (m2)

Δf = pre-detection bandwidth (Hz)

t = post-detection integration time (seconds)

n = number of records averaged

How to Fold Consecutive Records into One Record

Although the above figure shows the pulses neatly placed at the beginning of each record, in practice this is not possible for small systems - simply because the pulses are hidden in the noise in the original recorded data and so their position in time is unknown.  Fortunately this does not affect the folding process in terms of sensitivity, but merely means the single integrated record will show the pulse somewhere between the start and end of the record.  The folding software can re-position the pulse (by rotating the data in the record) and the standard procedure is to place it in the middle of the displayed integrated record.

Basic Folding Algorithm

The output record display is divided into 'time bins'.  Typically the number of time bins, 'N', is 256.  To make the maths neater I will use N = 250. This means that for our 1 second (1000 ms) imaginary pulsar, each bin is 4 ms wide.  For this example the sampling frequency is, say, 1000 sps (one sample every 1 ms).  This means each time bin contains 4 samples each out of the 1000 samples covering one record (which is 1 second long).  Note that the 50 ms pulse will be spread over about 12 time bins (e.g. - 50 ms/4 ms).

The process goes as follows...

Let's assume we have collected 1 hour (3600 seconds) of data of our 1 second pulsar.   This means there should be 3600 'records' of the pulse buried in the data.   We take the first 4 samples of the data and add them to the first time bin.  We then take the next 4 samples of the data and add them to the next time bin.  We keep on adding 4 samples to consecutive time bins until we have reached the last time bin (#250).   As we have run out of time bins we go back to the first time bin and add the next 4 samples to the samples already added to that time bin, and then repeat the process until we have reached the end of the data.  In this way the data is 'folded' into time bins exactly matching the period of the pulse.  The noise averages-down and, hopefully, our pulsar pulse emerges.

Fast Folding Algorithm

The folding operation is CPU-intensive and a 'Fast Folding Algorithm' has been developed over half-a-century ago (Staelin - 1969) to reduce computation times.  The 'Fast Folding Algorithm' (FFA) is to Epoch Folding what the 'Fast Fourier Transform' is to Fourier Transforms - that is, both simply speed up the calculations.

The FFA does NOT make the radio telescope more sensitive - that sensitivity is still governed by the pulsar radiometer equation - you simply get the results more quickly.

How Many Time Bins ?

The number of time bins used has a significant effect on the results.  All the sample values from the data file get distributed across the time bins.  So, for a file of 106 samples, and 100 time bins, each bin would have 104 samples in it.  However, for 200 time bins, the number of samples integrated in each time bin is halved (0.5 x 104).  Therefore the S/N in each time bin is reduced by a factor of √2 = 1.414 or about 1.5 dB.

So why not set the number of time bins to some low value ?  Well - there is a lower limit to the time bins because if the width (in time) of the time bins is greater than the pulsar pulse width, not only is the pulse profile detail lost, but also the S/N starts to be degraded.  To get the best S/N the optimum situation is to just fit the pulsar pulse into one time bin.  In this condition you get all the energy from the pulsar pulse into one time bin and at the same time exclude, from that time bin, those parts of the period where there is just noise.

If the strongest pulsars @ 400 MHz have their P0 to W50 ratio values calculated (C1), an estimate of the optimum number of time bins for best S/N can be made...

Note that the number of time bins indicated by C1 (rounded to an integer) is for best S/N.  In terms of getting detail of the pulse shape these numbers are too low - more plotting points (time bins) are needed to detail the rise and fall of the pulse as well as any pre- or post-pulse detail.

Usually the number of time bins is chosen to be around 250 as a compromise between best S/N and a reasonable amount of detail of the pulse profile.

Note also - if it is decided to go for best S/N and choose the number of time bins such that a time bin width exactly matches the W50 of the pulsar pulse, then it would be useful to be able to adjust the folding phase by ±0.5 time bins to ensure the pulse peak sits right in the middle of a time bin to extract the last bit of S/N.

Some real data (Vela data supplied by Guillermo Gancio & Peter East) was tested to find the optimum number of time bins.  Time bin numbers from 10 bins to 256 bins were tested in increments of 1.   At each 'time bin number' test, the starting phase was varied over a time equivalent to the width of one time bin for that test, and the best S/N found was recorded.  This phase search was done to ensure, for each test, the pulse peak is central to a time bin - rather than perhaps straddling two time bins which would reduce the calculated S/N.

The results were plotted producing the following graph...

...which agrees roughly with estimate from the table above for Vela (B0833-45).  While the estimate from the table for best number of time bins is 43, the graph shows a peak S/N at about 52.  Nonetheless, the shape of the curve follows the principle described above where the S/N improves with a smaller number of time bins - up to the point where the width of the time bin grows larger than roughly the width of the pulsar pulse (i.e., number of time bins < 43).

Note: this plot agrees in general shape with a prior independent test done by Peter East on the same data using different software.

Note also the decrease in S/N when going from 100 bins to 200 bins as predicted above.  The calculation above gave a drop of 1.5 dB, the graph shows a drop of about 1.3 dB - a reasonable agreement.

Some 'Flies in the Ointment' (Ecclesiastes 10:1)

The above example of the folding algorithm uses neat numbers and a few assumptions - which need to be addressed in the folding software.

Accuracy and Stability of the Sampling Hardware

It should be obvious that if the sampling frequency is not known accurately and the period of the pulsar is not known accurately then samples are not going to be put in the right bins over the course of the 1 hour data analysis and the pulse will not rise out of the noise.  This can be overcome by making the folding software able to 'search' for the best match by adjusting the software sampling time value and/or the software assumed pulsar pulse period.  For a wide pulse like our imaginary 1 second pulsar (50 ms), the search steps for a 1 hour record can be of the order of 10 ppm.  However for some pulsars the pulse width is much narrower (e.g., 5 ms) and so the step must be correspondingly smaller (1 ppm).  This can be time-consuming if done manually, and automation is useful - let the computer run while you have a cup of coffee.

The accuracy side of things can be catered for by some sort of search process as described, but stability is a bigger problem.  If at the start of the 1 hour data run the hardware is sampling at 1000.00 sps and after 15 minutes (after heating up for example) it is sampling at 1000.01 sps, then after 1 hour of data the folding mechanism will be out by about 9 time bins.  As the pulse only occupies about 12 time bins, this misalignment begins to degrade the sensitivity.  Once again for pulsars with narrow pulses the degradation is much more severe.

Note that the adjustment detailed above for compensating for inaccuracies will not correct the stability (drift) problem - because any adjustment to line up with the sampling rate at one point of the data will cause a misalignment somewhere else.

So - bottom line - there are enough variables to chase without chasing drift.  A stable (and hopefully accurate) clock source for the sampling hardware is critical.  Some use TCXOs - the author uses one of the those secondhand Rubidium sources on eBay.

Spindown and Doppler Shift on the Pulsar Signal

After ensuring you know the sampling frequency accurately and also ensuring it is stable (at least to 0.1 ppm is advised), there remains the information about the current period of the target pulsar's pulse.  It is simply not sufficient to look up a database (say the ATNF Pulsar database) and use that value.

Firstly, the value of the period of the pulse quoted there may have been measured a decade or more ago.  All pulsars 'spin down', i.e., the pulse period lengthens over time and so invariably the the pulsar is spinning slower than the value quoted in the database.  For example, the value of the spin frequency (1/period of the pulse) of the Vela given in the ATNF database is 11.194650 Hz (dated MJD 51559.32 - about 15 years ago), while the current rotational frequency is 11.187242 Hz - a difference of about 600 ppm.  If you are using 0.1 ppm steps in your search (Vela has a pulse width of 2.1 ms and so needs small steps) to match the pulse period, then it is going to take a long time to find the pulse !!!

Secondly - even if you have used the coefficients in the database to calculate the spindown frequency, that result is the intrinsic value - that is, the spin frequency that would be observed from a fixed point in space.  Any observatory on Earth is not fixed in space, but moving through space with a velocity which is the sum of the Earth's velocity around its orbit of the Sun and the velocity imparted to the observatory by the daily rotation of the Earth on its axis.  This in turn imparts a Doppler shift on the intrinsic spin frequency - the exact magnitude and sign of which is determined by the current state of those velocity vectors as well as the direction in space of the pulsar and the observatory location.

The graphic above shows the combination of intrinsic spindown and the annual doppler shift of Vela due to Earth's orbit around the Sun as viewed from 34oS (HawkRAO).  The peak to peak annual doppler swing is about 100 ppm.  Note that the rate of change of the observed topocentric annual doppler swing for Vela is zero around June and November.  Note also that the topocentric frequency actually increases between June and November (after the diurnal doppler shift is averaged out over a day).

The diurnal doppler shift on Vela due to the daily rotation of the Earth on its axis is about 1.8 ppm peak-to-peak as observed from 34oS (HawkRAO).  The figure above shows an example where the calculation was done on a rising part of the annual doppler shift curve (October).  Note how the topocentric frequency is falling during the time from 6 hours before transit until 6 hours after transit.

For narrow pulse widths and long observation times, a S/N degradation occurs if the frequency is assumed to be constant over the duration of the observation.  This amounts to about 0.6 dB for Vela (pulse width 2.1 ms) and an observation time of 120 minutes.  Narrower pulse widths and longer observation times are affected to a greater degree.  Naturally, the converse is true.

The typical solution to this is to run TEMPO which gives the current spin frequency at any time and from any location on Earth.  While the intrinsic spindown and annual doppler shift (a function of the position of the pulsar w.r.t. the plane of the Earth's orbit) for a particular pulsar is the same for all Earth-bound observers, the magnitude of the diurnal doppler shift depends on the combination of the declination of the pulsar and the latitude of the observatory and is a function of the product Cos(DEC)*Cos(LAT).  The values of pulsar RA and observatory longitude combine to determine the phase of the diurnal doppler shift.

The result of all the calculations is called the 'Topocentric Frequency'.

Most of the successful amateurs listed on the home page use some version of TEMPO (linux or Windows) to make predictions for their target pulsars.

The use of TEMPO for predictions of topocentric frequency is more than accurate for amateur purposes, except where the pulsar exhibits significant glitch behaviour.  A discussion of this topic is found here "Accounting for Vela's Glitch Behaviour".

The above command line is called from a "Command Window" in a Windows machine.  The RTL SDR raw I/Q data file needs to reside in the same directory as the executable.

The output text file ("test.txt") can be imported into a spreadsheet and the first two columns can be plotted to get the pulse profile as shown in the top plot on the right (period = 89.39 ms).

The sample rate for the data file is 2M sps and 256 time bins has been selected.   The expected pulse period is 89.39 ms is input.  This pulse period information was recorded along with the file documentation at the time of the observation after calculations were made to obtain the current 'Topocentric Period'.

A deviation from the optimum input period of 89.39 ms is tolerated well for this high S/N data file.  A deviation of 0.01 ms (i.e., 89.38 ms - about 112 ppm) from the optimum pulse period input on the command line still produces a pulse as shown in the bottom plot on the right, but at a reduced S/N because of the spreading of the pulse across adjacent time bins.  For data where the S/N is lower (small dishes) an input period closer to the optimum will be needed to dig out a pulse profile.

The RTL SDR raw I/Q data file (20150420_185141_1413405_157_100.bin - about 190 MB size) can be downloaded from here...

https://drive.google.com/open?id=0B3WMi2fu54AMQXE5Q3p6NjFpOVU

The raw data file was obtained from a 50 second observation of the Vela pulsar (J0835-4510) by the 26 metre diameter dish at IAR by Guillermo Gancio.

The source code for 'rapulsar2.exe' can be found in the *.ZIP file referenced above (currently called "NewSW.zip").

The "RTL SDR Large File Software Tools" document resides on Peter's website, a site which contains a lot of useful radio astronomy information.

Thanks to Peter and Guillermo for making this information available. 

I hope that others can use this information to develop their own software.