DSPSR: Digital Signal Processing Software for Pulsar Astronomy

W. van Straten; M. Bailes

Source-linked AI summary

DSPSR: Digital Signal Processing Software for Pulsar Astronomy

W. van Straten, M. Bailes

arXiv:1008.3973v1Aug 24, 2010astro-ph.IM

TL;DR

Radio pulsar observations need greater sensitivity from wider bandwidth without losing information to interstellar dispersion, while practical processing must handle large data streams. The paper presents DSPSR, a modular software suite using coherent-dedispersion and other signal-processing algorithms across CPUs and GPUs. DSPSR is shown to process contemporary signals in real time using four GPU-equipped workstations, while the paper also identifies limits of a deprecated filterbank method and future computational scope.

Problem
Wider bandwidth improves pulsar sensitivity but introduces interstellar dispersion, requiring processing methods that preserve temporal information while handling demanding data rates.
Method
The paper describes DSPSR’s modular radio-astronomy algorithms, including coherent dedispersion, and evaluates their implementation on multi-core CPUs and GPUs.
Results
Four workstations, each equipped with a single GPU, can comfortably process the demonstrated signals in real time.
Takeaways & Limitations
DSPSR is presented as a mature, extensively functional, open-source library for processing radio-astronomical data from diverse observatories and formats.
Takeaways & Limitations
The deprecated synthetic-filterbank method reduces coherent-dedispersion cost but suffers spectral leakage that can produce delayed pulse-profile images when inter-channel delays are large.

Abstract

from arXiv · show

DSPSR is a high-performance, open-source, object-oriented, digital signal processing software library and application suite for use in radio pulsar astronomy. Written primarily in C++, the library implements an extensive range of modular algorithms that can optionally exploit both multiple-core processors and general-purpose graphics processing units. After over a decade of research and development, DSPSR is now stable and in widespread use in the community. This paper presents a detailed description of its functionality, justification of major design decisions, analysis of phase-coherent dispersion removal algorithms, and demonstration of performance on some contemporary microprocessor architectures.

1 Introduction

Pulsar astronomy seeks greater sensitivity through wider bandwidth, but interstellar-medium dispersion increasingly distorts the signal. DSPSR builds on evolving baseband processing, coherent-dedispersion, and GPU technologies to address these demands with modular algorithms.

Wider instrumental bandwidth increases sensitivity but also increases signal distortion, primarily from plasma dispersion in the interstellar medium.
Post-detection dispersion removal is limited by residual intra-channel smearing and reduced temporal resolution, motivating phase-coherent techniques.
Baseband systems evolved from tape-based recorders toward quasi-real-time processing of increasingly wide bandwidths using commercial computing and storage.
DSPSR implements phase-coherent dispersion removal, synthetic filterbanks, polarimetric processing, folding, sub-integrations, and simultaneous multi-pulsar processing.
GPU developments led DSPSR to adopt CUDA, extending its processing capabilities alongside contemporary multi-core systems.

2 Algorithms and Features

DSPSR processes continuous radio-pulsar data streams into integrated statistics through a modular signal-processing chain. Its salient features include invalid-data excision and correction of digitisation distortions.

DSPSR processes continuous streams of radio-pulsar data to produce integrated statistics such as phase-resolved average polarisation.
Its processing features include automatic invalid-data excision and dynamic output-level setting to correct digitisation distortions.
The paper describes these features in detail while outlining a typical pulsar signal-processing path.

2.1 Invalid Data Excision

DSPSR propagates validity information through the signal-processing chain so later operations can ignore corrupted data. This prevents invalid segments from contaminating integrated results.

Any operation can flag data sections as invalid for exclusion by subsequent signal-processing components.The mechanism uses weights maintained in parallel with the astronomical signal data.
Weights typically flag time-domain segments across all frequency channels and polarisations for a given epoch.
Conversion classes monitor signal statistics and flag segments when values such as noise power fall outside acceptable operating ranges.
Integration operations use the weights array to prevent invalid data from corrupting final results.Examples include data recorded before instrumentation has properly initialized and broadband impulsive interference.

2.2 Dynamic Output Level Setting

DSPSR corrects quantisation distortion during conversion from 2-bit samples to floating-point voltages by estimating undigitised power from segment statistics. It uses precomputed lookup values to avoid repeated iterative solutions and can flag unacceptable segments.

Dynamic output-level setting restores the linear relationship between digitised and undigitised power to correct quantisation distortion.The correction addresses noise and signal distortion introduced by non-linear analog-to-digital conversion.
Less than 0.05% of APSR data was discarded during a lightning storm, substantially improving the resulting 2-minute integration.
During 2-bit-to-floating-point conversion, DSPSR divides samples into segments, counts low-voltage states M, and estimates undigitised power σ^2.
DSPSR solves for σ with Newton-Raphson but replaces per-segment iteration with a lookup table because only L+1 values of M are possible.
Segments with M outside the acceptable range are flagged invalid, while an M histogram is archived as a diagnostic of baseband-recording quality.

2.3 Coherent Dedispersion

Coherent dedispersion removes interstellar dispersion by deconvolving the observed signal in the frequency domain, while overlap-save convolution handles cyclic wrap-around and FFT length is selected for real-time efficiency.

The ISM imposes frequency-dependent delays, so dspsr obtains a dispersion response and deconvolves the observed radio signal with its inverse.The dispersion measure is the line-of-sight integral of free-electron density.
Frequency-domain deconvolution multiplies the signal by the inverse response, producing cyclical convolution whose wrap-around samples must be discarded.The discarded samples arise because the discrete Fourier transform assumes periodicity.
Because the impulse response is asymmetric, the numbers of discarded samples differ on the two sides of the convolution result.The asymmetry follows from the quadratic frequency dependence of dispersion delay.
N must exceed the discarded-sample count nd, and its value is chosen by minimizing the estimated real-time floating-point operations.For power-of-two N, the FFT operation count asymptotically approaches 5N log2 N; actual benchmarks can improve the choice.
Implementations can use FFTW, Intel Math Kernel Library, or Intel IPP, with aligned arrays enabling vector instruction sets.The software also flags invalid input segments so subsequent convolution outputs can be excluded from integrations.
The overlap-save method processes overlapping N-point segments, transforms each segment, removes polluted outputs, and copies the remaining samples into the output stream.The overlap equals the number of non-zero points in the dispersion impulse response.

2.4 Synthetic Filterbank

Synthetic filterbanks trade some time resolution for frequency resolution to reduce coherent-dedispersion cost, with the convolving method addressing spectral leakage that affects the deprecated approach.

2.4 Synthetic Filterbank: Dividing the band into sub-bands reduces sampling rate, smearing time, and the transform length required for subsequent coherent dedispersion.The resulting efficiency gain trades some temporal resolution for frequency resolution.
2.4.1 Deprecated Filterbank: The deprecated filterbank separately dedisperses Nc non-overlapping channels using channel-specific response functions, but spectral leakage mixes power between neighboring channels.When inter-channel delays are large, leakage produces delayed pulse-profile images.
2.4.2 Convolving Filterbank: The convolving filterbank uses longer K-point transforms, divides their spectra into Nc sub-bands, and applies a unique response function followed by an inverse FFT to each sub-band.Here K = Nc × N′, and overlap-save discards polluted samples before producing channel outputs.
2.4.2 Convolving Filterbank: Coherent dedispersion can run simultaneously with filterbank formation, with N′ chosen large enough for the maximum sub-band smearing time.Each sub-band is multiplied by its own dedispersion frequency response in the Fourier domain.
2.4 Synthetic Filterbank: Both deprecated and convolving filterbank methods require the same number of floating-point operations per second for phase-coherent dispersion removal.The number of channels should be as large as possible without sacrificing required time resolution, subject to efficiency limits.
2.4.2 Convolving Filterbank: The convolving filterbank eliminates the spectral-leakage artifacts produced by the deprecated method, although upstream filterbank leakage can remain.In CASPSR, neighboring-channel delays of approximately 0.11–0.27 µs are below the folded profile's temporal resolution.

2.5 Detection

DSPSR computes pulse-phase-resolved polarisation from dual-receptor voltage signals using the coherency matrix and its relationship to the four Stokes parameters.

DSPSR computes electromagnetic polarisation when voltage signals from both receptors of a dual-polarisation receiver are available.
The coherency matrix ρ represents second-order statistics of the transverse electric-field vector and is computed independently in each frequency channel for synthetic filterbank data.
A geometric relationship maps the complex 2 × 2 coherency matrix to the real four-dimensional space of the Stokes parameters.
The mapping uses the identity matrix, Pauli matrices, Einstein summation, and the matrix trace operator.

2.6 Folding

DSPSR folds continuous pulsar data into phase-resolved profiles using a timing model, while supporting parallel execution, flexible sub-integrations, single pulses, and multiple pulsars.

2.6 Folding: Folding integrates samples into n equally spaced phase bins according to apparent topocentric pulse phase from a polynomial pulsar timing model.The phase model may be generated with tempo or tempo2.
2.6 Folding: The start and end observation times determine the mean time assigned to the rising edge of phase bin zero, assuming sufficiently accurate phase calculations.
2.6 Folding: DSPSR creates profiles for arbitrary time segments, enabling regular sub-integrations that can later be recombined after timing updates or corrupted-segment deletion.
2.6 Folding: Phase-aligned segment boundaries enable single-pulse profiles; dispersion-delay removal and phase-reference adjustment keep each pulse aligned across channels and integrations.
2.6 Folding: For short-period pulsars, DSPSR supports reducing or conditionally retaining single-pulse outputs when full frequency-resolution storage is impractical.
2.6 Folding: User-supplied psrchive scripts support post-processing such as polarisation calibration, frequency integration, and threshold-based giant-pulse selection.
2.6 Folding: Multiple pulsars can be folded simultaneously by branching the signal path at the final folding stage, with a shared mean or median dispersion measure for coherent removal.
2.6 Folding: Multiple threads fold separate data blocks into independent integrations before combining results, while invalid flagged data can be omitted from integration.

2.7 Thread-safe Input Buffering

DSPSR gives each sample-losing operation independent overlap management through thread-safe buffering, preserving contiguous inputs while reducing pipeline interdependency and reprocessing.

2.7 Thread-safe Input Buffering: Signal-processing operations may lose samples, creating incompatible block-size and overlap constraints when multiple such operations share a chain.
2.7 Thread-safe Input Buffering: Each operation independently manages its overlap requirements with a thread-safe input-buffering scheme.
2.7 Thread-safe Input Buffering: The design increases modularity and efficiency by reducing interdependency between components and avoiding reprocessing of data lost later in the signal path.
2.7 Thread-safe Input Buffering: Operations buffer the appropriate tail of each input block and prefix it to the next block, with bookkeeping preserving block contiguity.
2.7 Thread-safe Input Buffering: With multiple processing threads, input blocks are distributed across separate signal-processing paths.

3 Performance

DSPSR performance depends on FFT implementation, cache and memory behavior, processor threading, and workload configuration. Benchmarks show how these factors guide CPU/GPU processing choices and estimate hardware requirements for real-time pulsar signal processing.

3.1 Central Processing Units: FFT execution dominates coherent dedispersion and filterbank synthesis, making FFT-library performance central to DSPSR throughput.Benchmarks compare FFTW and IPP across transform lengths and thread counts.
3.1 Central Processing Units: Cache capacity creates distinct performance transitions: FFT speed drops after L1 and L2 data limits, with eight threads exposing the shared-L2 constraint at NFFT = 217 = 128 k.The E5345 uses private L1 caches but shared L2 caches, so thread scaling depends on transform size.
3.1 Central Processing Units: Measured FFT timings let DSPSR choose the optimal transform length and FFT library instead of relying on the O(N log2 N) estimate, improving shorter-transform performance.The selection minimizes the processing equations using measured N-point FFT times.
3.1 Central Processing Units: At low dispersion measure, CPU throughput scales roughly linearly with thread count, whereas at high dispersion measure eight threads can underperform four because of L2-cache competition.The high-DM regime uses larger FFTs, making shared-cache effects more consequential.
3.1 Central Processing Units: At the highest dispersion measure, processing on the E5520 is as much as 4 times faster than on the E5345, indicating memory bandwidth limits the latter when FFTs exceed L2 cache.At the lowest dispersion measure, performance is not significantly improved, consistent with processor-speed limitation in that regime.
3.2 Graphics Processing Units: GPU filterbank performance improves when Nc × NFFT is large, while a 512 MHz band processed to DM 1000 in real time requires only 4 GPU-equipped workstations.Larger combined workloads improve GPU utilization; the benchmarked pulsar had a spin period of approximately 5.76 ms.

4 Conclusion

dspsr supports a diverse range of pulsar and radio-astronomy experiments, handles common voltage- and power-level data formats, and provides real-time processing for contemporary observing configurations. Its future development includes more portable GPU programming and potentially new statistical measurements.

dspsr has enabled experiments spanning giant micropulses, sub-millisecond pulsar surveys, precision timing, polarimetry, gravitational-wave searches, and interstellar-plasma studies.
dspsr can comfortably process the demonstrated signals in real time using four workstations, each equipped with a single GPU.
dspsr processes voltage- or power-level signals and understands data recorded at most major radio observatories.It supports formats including VDIF, PSRFITS, and sigproc.
The library’s algorithms can also be applied to single-beam radio observations beyond pulsar analysis, including high-resolution spectra for extragalactic neutral-hydrogen studies.
Future developments may adopt OpenCL for more portable GPU programming and make computationally intensive statistics feasible to observe regularly.The paper gives cyclostationary statistics and cyclic power spectra as examples connected to scintillation and multipath-scattering studies.

A Floating Point Operations

The floating-point cost of synthetic filterbank formation and phase-coherent dispersion removal depends on channel count and dispersion-response size, while practical speed depends on FFT implementation and hardware architecture. The convolving filterbank is recommended despite cache-dependent performance trade-offs because the deprecated method introduces spectral-leakage artifacts.

The operation count depends on the number of synthetic filterbank channels, Nc, and the dispersion-response size N′ in the lowest-frequency output channel.This applies to both real-valued and complex-valued input data.
Both filterbank methods require a total of 5NcN′(log2 Nc + 2 log2 N′) floating-point operations.
FFT implementation details and computing-device architecture determine how quickly the required floating-point operations are performed.FFT libraries typically become more efficient with larger transforms until the problem exceeds L1-cache capacity.
When Nc is below the L1-cache-limited transform length, the convolving filterbank may outperform the deprecated method.
When Nc × N′ exceeds L2-cache capacity, the deprecated method may perform better, but the convolving filterbank is recommended because the deprecated method introduces spectral-leakage artifacts.