A Coherent Harmonic Summing Pulsar Search Code

The gist: CoherentSearch.jl is an open-source Julia code that reconstructs the pulse profile at every trial spin frequency directly from the Fourier transform, offering superior sensitivity and noise resilience compared to standard incoherent harmonic summing methods.

The core problem addressed

Pulsar searches typically use incoherent harmonic summing, which adds powers of many harmonics to improve sensitivity but discards the pulse shape by discarding phases, which throws away some of the signal. Folding the time series at each trial period keeps both the pulse shape and its phase, but this is computationally expensive even with algorithms like Fast Folding Algorithm (FFA). The paper describes CoherentSearch.jl, which reconstructs the pulse profile at every trial spin frequency directly from the Fourier transform.

The Coherent Search Algorithm

The goal of the algorithm is to produce a pulse profile as a function of trial spin frequency f, and then to ask two things of each profile: is there a significant pulsation, and is it pulsar-like. The process involves several steps:

1. Start with a normalized FFT of a long, evenly sampled time series which is stepped finely through Fourier frequencies r = fT.

2. At each trial r, determine the complex Fourier amplitudes of each of H harmonics hr.

3. Inverse FFT them into a real-valued pulse profile.

4. Compute a detection statistic for both significance and whether the pulse is pulsarlike.

5. Reuse subsets of the same H amplitudes to build lower-resolution profiles at integer multiples of r, and test those with the same statistic.

Key algorithmic details

The algorithm utilizes Fourier interpolation to obtain the complex amplitude of each harmonic at each trial frequency, solving both the phase accuracy and amplitude loss problems inherent in simple nearest-bin reading. The interpolation formula is given by Ar ≃ [r]+ Sum m/2 k=[r]−m/2 Ak e −iπ(r−k) sinc [π(r − k)], where Ak are the complex Fourier amplitudes at the nearby integer Fourier frequencies k, sinc(x) = sin(x)/x is the normalized sinc function, [r] is the nearest integer to r, and m is the total number of nearby FFT bins used for the interpolation. The trial grid spacing δr = 1/(2H) ensures that harmonic H advances by half a Fourier bin per step, which is analogous to how FFA folds a base period into b bins and produces trials spaced by 1/b Fourier bins.

Efficiency and Performance

CoherentSearch.jl leverages Julia's capabilities, compiling to native machine code for fast numerical loops while supporting SIMD instructions, threads across cores, and GPUs through CUDA.jl. The search is computationally efficient because the input spectrum is normalized, meaning the noise in every profile is known in closed form. The cost of interpolation per harmonic involves m realtimes-complex multiply-adds and reads only m consecutive bins of the FFT. Harmonic decimation, which selects every kth amplitude to create a profile with nk = 2Hk bins, costs one short inverse FFT and one matchedfilter pass per decimation factor.

Comparison with Existing Codes

The paper compares CoherentSearch.jl against PRESTO’s accelsearch, which uses an incoherent harmonic summing search, and riptide’s rseek, which is the state-of-the-art FFA implementation. In an injection study of over 860,000 simulated pulsars at a matched false-alarm rate, CoherentSearch.jl recovers 76% of the pulsars in white noise, against 64% for riptide’s FFA in its recommended configuration and 42% for PRESTO’s accelsearch. On one CPU core it is 1.4–2.2 times faster than rptide at matched frequency coverage, scaling to 27 times faster on 48 cores.

Noise Resilience and Thresholds

A key advantage of CoherentSearch.jl is that because the input spectrum is normalized, one threshold gives one false-alarm rate across the whole search. The analytic noise scale σ tracks the exact noise level to about 1%, and this removes a source of variance from the statistic that we return to in §3.2.2. Furthermore, CoherentSearch.jl’s detection threshold does not move under red noise, its candidate lists are an order of magnitude shorter, and it is faster than rseek single-threaded.

Conclusion

CoherentSearch.jl recovers 76% of the population in white noise against a false-alarm rate of 0.1 per observation, demonstrating superior sensitivity and resilience to red noise compared to existing codes. The method is not specific to pulsars, as it can be adapted for any periodic signal whose shape is narrow compared with its period.

How it works

The search walks a grid of Fourier frequencies r = fT at a grid of Fourier frequencies rlo, rlo + δr,..., rhi do ▷ Fourier frequency trials. At each trial, the code interpolates the complex Fourier amplitude of each of H harmonics Ah ← interpolate(F, h r, m) ▷ Fourier interpolation. It then stacks those amplitudes into a half-complex spectrum with the DC (i.e. zerofrequency) term held at zero, and inverse transforms them into a real pulse profile of n = 2H phase bins. For each profile, it computes a detection statistic for both significance and whether the pulse is pulsarlike.

The search metric

For a profile Pj with n bins, the S/N statistic is given by S/N = max w,p Sw(p) − δ Stot σ p w(1 − δ), where Stot is the sum of the whole profile and σ the noise per bin. The denominator uses the exact variance of the boxcar sum instead of assuming independence between bins, which is more accurate for band-limited noise.

Parallelism

The frequency trials are independent, so the search parallelizes in three ways: SIMD instructions within a core, threads across the cores of a machine, and GPUs through the CUDA.jl package. Groups of consecutive r are dealt to threads round-robin, with each having a private workspace.

Conclusion

CoherentSearch.jl takes the FFT search’s input, one long transform, and does the FFA’s job in the Fourier domain by interpolating complex amplitudes of 60 harmonics at each trial frequency. It then inverts them into a pulse profile and matched filters that profile against boxcars, reusing the same harmonics to fold at two to six times each trial frequency almost for free. Because the input spectrum is normalized, the noise in every profile is known in closed form, and one threshold gives one false-alarm rate across all fold depths and all frequencies. On the same injections it recovers 76% of the population, against 64% for riptide’s FFA in its recommended configuration at about three times the cost. The one place its default configuration loses is to a much deeper FFA fold at duty cycles below 0.5%, and a second pass at twice the depth over the slowest frequencies makes it the best performing search for a fifth of the FFA’s cost.

1. No longer discard pulse shape information in initial searches by reconstructing it directly from harmonics: Folding the time series at each trial period keeps both [powers and phases], but is expensive even with the Fast Folding Algorithm (FFA). This allows for a more complete characterization of signals that are not purely sinusoidal.

2. Enable closed-form noise statistics for every profile, eliminating per-profile scatter: Because the input spectrum is normalized, the noise in every profile is known in closed form, so one threshold gives one false-alarm rate across the whole search. This provides a stable and reproducible detection threshold across all trial frequencies and fold depths.

3. Improve sensitivity by using Fourier interpolation for phase accuracy: Fourier interpolation solves both problems [phase inaccuracy and amplitude loss] effectively, ensuring that the reconstructed pulse profile accurately reflects the true signal at every trial frequency, unlike methods that lose amplitude when a signal sits between bins.

4. Implement efficient harmonic decimation for high-frequency searches: Selecting every kth amplitude gives the harmonic stack of a signal at kr with Hk = ⌊H/k⌋ harmonics, and one inverse FFT of length Hk turns it into a profile with nk = 2Hk bins, enabling fast access to spin frequencies beyond the base grid.

5. Optimize computational efficiency via SIMD and GPU utilization: The code leverages Julia's parallelism, allowing it to run on a GPU through the CUDA.jl package and utilizes SIMD instructions, achieving speedups up to nine times faster than a 20-core workstation. This enables rapid screening of large datasets across modern hardware.

6. Use analytic noise scaling for robust detection statistics: The search uses an analytic noise scale rather than profile estimation: We compute σ rather than estimate it from each profile, which turns out to be both faster and better, removing a source of variance that otherwise inflates the reported signal-to-noise ratio.

7. Develop a dynamic thresholding mechanism for real-world noise: Instead of using a single threshold, the system can match the rate within bins of knee frequency or use a dynamic threshold like the one its authors describe (V. Morello et al. 2020) to manage false alarms effectively under varying red noise conditions.

8. Create parallel search strategies for large surveys: The code includes a helper script that partitions DMs across multiple jobs, allowing the start-up costs being paid only once, enabling efficient searching of large, unknown DM ranges using single-threaded processes.

Pulsars have narrow pulses, so their signals spread over many harmonics in a Fourier transform. Standard pulsar searches add the powers of those harmonics and discard their phases, which throws away the pulse shape and some of the signal. Folding the time series at each trial period keeps both, but is expensive even with the Fast Folding Algorithm (FFA). We describe CoherentSearch.jl, an open-source Julia code that reconstructs the pulse profile at every trial spin frequency directly from the Fourier transform. It interpolates the complex amplitudes of typically 60 harmonics, inverse transforms them into a profile, and tests that profile against boxcar templates of many widths. Because the input spectrum is normalized, the noise in every profile is known in closed form, so one threshold gives one false-alarm rate across the whole search. In an injection study of over 860,000 simulated pulsars at a matched false-alarm rate, the new search recovers 76% of the pulsars in white noise, against 64% for riptide's FFA in its recommended configuration and 42% for PRESTO's accelsearch, and its threshold does not move under red noise. On one CPU core it is 1.4-2.2 times faster than riptide at matched frequency coverage, it scales to 27 times faster on 48 cores, and modern GPUs run it up to nine times faster than a 20-core workstation.