A Better Pulsar Search Code

CoherentSearch.jl // Oxford, 5 October 2026

A Better Pulsar
Search Code

my first real research project with Claude

Scott Ransom

02 / Pulsar signals

Pulsar signals via an FFT

Pulses are narrow: duty cycles of a few %, some well under 1%
A narrow pulse is built from many harmonics, each with a fixed phase
Repeat it every P seconds and an FFT of the time series puts those harmonics at f, 2f, 3f, …

A periodic pulse train in a long time series

P = 1/f

PSR J2019+2425 and its first 10 harmonics: intrinsic (a, b) and observed (c, d). Nice, Fruchter & Taylor 1995, ApJ 449, 156

03 / FFT searches

Traditional FFT searches

Convolution / modulation theorem: harmonics land at multiples of the spin frequency
Convert the complex harmonics to powers, then sum 1, 2, 4, 8, 16 of them
“Incoherent” harmonic summing: phases are lost, and so is the profile shape
Well-understood χ² statistics
Very fast, but other issues (e.g. scalloping: up to 36% of power lost between integer bins)

Courtesy of Ewan Barr

04 / Folding

The Fast Folding Algorithm

Known since the early days (Staelin 1969)
Keeps phase, and so pulse shape: more sensitive to pulsar-like signals

goodRe-uses partial sumsMuch faster than brute-force folding: N log N, like an FFT

badTricky and slowHard to implement and to apply to real data, and much slower than incoherent FFT searches

state of the artriptideMorello et al. 2020, MNRAS 497, 4654. An FFA search that is fast and easy to use

Song et al. 2023, Chinese Journal of Aeronautics 36, 309

05 / Goals

What would a new search do?

01SensitiveAs sensitive as an FFA, or more. So it has to keep the pulse shape.

02FastFFT-based, and built for modern compute (SIMD, many cores, GPUs). So we need the harmonic phases: coherent harmonic summing.

03Red-noise resilientAt least as good as the best current techniques.

04Few false positivesA short, trustworthy candidate list.

~20 years in my headCoherent harmonic summing is not a new idea.PRESTO’s accelsearch has a “coherent power”. Nice et al. 1995 used phases to check candidates.I recently added fourier_fold.py to PRESTO for the same reason.Plus I wanted to learn how to use LLMs!

06 / Parameter space

Which pulsars do we want?

FFT searches struggle with slow pulsars and narrow duty cycles, especially in red noise
MeerKAT TPA: ~13% of 1176 pulsars are below 1% FWHM duty, ~2% below 0.5%
So 60 harmonics (a 120-bin profile) covers most of the population
But we still want to be very sensitive to the fastest MSPs (5–30% duty, up to 716 Hz)

Posselt et al. 2021, MNRAS 508, 4249 (supplementary W50 table)

07 / The idea

What is coherent harmonic summing?

01Grab the harmonicsInterpolate the exact complex amplitudes of H harmonics, phases and all, from one long FFT.

02Build a profileStack them in a vector and take an inverse real FFT: a pulse profile of 2H bins.

03Test for pulsar-nessMatch-filter the profile with a bank of boxcars at every width and phase.

It is folding — done in the Fourier domain, at every trial frequency.

08 / Step 01

Fourier interpolation

A_r complex amplitude at any frequency r = fTA_k FFT amplitudes at integer bins k[r] nearest integer to rm bins used (16): loses ~1.3% of the power

e.g. Ransom, Eikenberry & Middleditch 2002, AJ 124, 1788

Coefficients factor as A(Δr)/(Δr−j): a real weight per bin
Trial grid is rational, so only 2H offsets ever occur: tabulate the weights once
32 trials at a time in SIMD lanes, bins read from cache. Accurate to ~10⁻⁷

Interpolated power around a 1.158 Hz pulsar

09 / Step 02+

Harmonic decimation

Every k^th harmonic of r is the harmonic stack of a signal at k r
Speed: one short inverse FFT, no new interpolation. k = 1…6 gives 120…20 bins, and 0.1–125 Hz fundamentals reach 750 Hz
Sensitivity: it is the harmonic-sum ladder. A broad pulse is best seen in a shallow fold; one depth alone is ~8% worse
Higher harmonics nearly for free!

A pulsar at 3× the trial fundamental: the k = 3 fold finds it

10 / The search

Making this into a search

01Step over frequenciesFourier frequency trials r, spaced 1/(2H): the FFA’s own rule.

02Examine profilesInterpolate H harmonics once, then fold at r, 2r, … 6r.

03Compute S/NBest boxcar over width and phase, against one threshold.

The paper’s pseudocode. bin/toy_coherent_search.jl implements it line by line

11 / Implementation

CoherentSearch.jl: why Julia?

Compiles to native code: C/Fortran speed, readable source, works well with Python
SIMD, threads, and GPUs (CUDA.jl) all built in
One-command install. No compiler, no system libraries
And I wanted to learn it a bit!

validationPinned to PythonEverything is referenced to coherent_search, the Python/NumPy version I wrote: same amplitudes, profiles and statistic to ~10⁻¹⁶.

Who is doing the coding?All of the Julia was written by ClaudeThe majority by some version of Opus 5, under my direction, against that Python reference.

github.com/scottransom/CoherentSearch.jl

12 / Nutshell

So how does it do?

76%Sensitiveof injected pulsars found, vs 64% for riptide (recommended) and 42% for accelsearch

1.4–2.2×Fastfaster than riptide on one core. 27× on 48 cores; GPUs up to 9× a 20-core workstation

6.75Red-noise resilientour matched threshold, 6.70–6.75 out to a 50 Hz knee. rseek C (whitened): 7.55–7.65. riptide unwhitened: up to 250

12Few false positivesmedian candidates per observation (riptide config. A: ~160). Gaussian noise tail

All at a matched false-alarm rate of 0.1 per observation.

13 / Monte Carlo

The simulations

~870,000injected pulsarsin ~161,000 simulated observations: 2²⁴ samples at 60 µs (1006 s), six pulsars each

2/3 : 1/3slow : MSPslow: 20 ms–5 s with TPA widths and shapes. MSPs: 1.5–10 ms, 2–30% duty. S/N 5.5–11.5

0.1false alarms / obsEvery code gets the threshold that gives this. Identical noise for all codes; red paired with white

rseek AMatched work20–120 bins everywhere: our coverage and fold depths. 1.7× our time

rseek CRecommended120–140 bins above 7.2 ms. Suggested by Vincent Morello. 3.2× our time

rseek BDeep foldsUp to 960–1040 bins for slow pulsars, like riptide’s pipeline. 7.3× our time

alsoaccelsearch & prepfoldaccelsearch with 16 harmonics. prepfold at the known period: the best case

14 / Results

White noise

$white-noise detection fractions$

CoherentSearch.jl76.3%

rseek B71.4%

rseek C64.5%

rseek A50.2%

accelsearch41.6%

Only loss: rseek B’s 1000-bin folds below 0.5% duty
Deep pass: 120 harmonics, decimation to 12, below 5 Hz. +3.1 s; beats rseek B at every duty

Threshold matched to 0.1 false alarms per observation. prepfold, told the period, finds 99.4%

15 / Statistics

A calculable noise statistic

Normalized FFT ⇒ noise per profile bin is known:
σ = √(2n_low + ½n_nyq)
Exact boxcar variance for band-limited profiles: every depth and frequency on one scale
No per-profile σ estimate: no scatter (stdev of recovered S/N 1.0 vs 1.6 rseek, 1.9 prepfold)
Gaussian tail ⇒ threshold 6.75 vs 7.20 (accelsearch), 7.55 (rseek A, C)
Against rseek C, our whole lead is this lower threshold

False alarms in white noise, 22,368 simulated observations

16 / Red noise

What is red noise?

Slow baseline wander: gain drifts, atmosphere. Power spectrum 1 + (f/f_knee)^−α
Real Arecibo, GBT, & Parkes data: α ≈ 2 (random walk), knees of 1.6–31.5 Hz
It hurts every search, mostly for slow pulsars (Lazarus et al. 2015; van Heerden et al. 2017; Singh et al. 2022)

Exactly the red noise we injected: knee 1 Hz, α = 2, five draws

17 / Results

Red noise

Knee: where the red power falls to the white level. Higher knee = worse noise
Drawn log-uniform from 0.1 to 50 Hz, α = 2.0 ± 0.3
Whitened with PRESTO’s rednoise in the Fourier domain
We had to change how riptide is called: give it the whitened time series, running median off
Then both lose sensitivity in the same proportion (keep 63%): ours 76.1→48.1%, rseek C 61.3→38.3%

$red-noise detection fractions$

One threshold per observation, matched in each knee bin

18 / Sidenote

riptide’s default de-reddening

riptide unwhitened detection and threshold

riptide subtracts a running median (4-s window): a high-pass at ~0.25 Hz
Above a ~0.3 Hz knee its slow trials fill with false positives
Threshold for 0.1 false alarms climbs 7.6 → 250; one threshold per observation finds nothing
This is a detrending issue, not the FFA

riptide as its authors run it, on the unwhitened time series

19 / Speed

Speed on CPUs

One core, idle Xeon w5-3433: 10.8 s per 1006 s observation. rseek A 1.7×, C 3.2×, B 7.3× longer; accelsearch 14× faster

26.8× on 48 cores (2× EPYC 7413), serial fraction 0.017. On a real Parkes observation, 1.4–2.2× faster than rseek on three CPUs

20 / GPUs

The GPU version

card	wall (s)	ns / trial	vs 20 CPU cores
L40	4.30	10.2	9.3×
A100 (80 GB)	4.52	10.7	8.9×
RTX 4500 Ada	7.30	17.2	5.5×
RTX A4000	12.07	28.5	3.3×
GTX 1080	29.9	70.7	1.3×
RTX A400	61.4	145	0.65×

7.3-hr GBT observation, 0.1–133 Hz: 423M trials. CPU: 2× Xeon Silver 4114, 40.1 s. Same 285 candidates on every card

--gpu, NVIDIA only, via CUDA.jl
A weak dependency: CPU-only users download nothing
Why so well? Every trial is independent, identical work: one thread per trial, no recursion
A block of threads shares one window of FFT bins; batched cuFFT for the profiles
Big CPU boxes still compete: 48 EPYC cores beat their own A4000

21 / Summary

CoherentSearch.jl: it works really well!

sensitive76% vs 64%of a realistic population, vs riptide recommended. 42% for accelsearch.

fast1.4–2.2× riptideon one core; 27× on 48 cores; ~150–200× on an L40 or A100

red-noise resilientFlat threshold6.70–6.75 to a 50 Hz knee, vs 7.55–7.65 for whitened rseek C.

few false positives12 per obs.median candidates, vs ~160 for riptide A. Gaussian tail.

This has been a success!We need to write a paper!

https://github.com/scottransom/CoherentSearch.jl

22 / Process

Writing a paper (and talk!) with Claude

Letting Claude write the first draft has been contentious among colleagues.
For this project it has been an excellent tool, and I will keep using it this way.
We iterated on the analysis many times. I wrote a detailed outline, described the tone and style I wanted, and gave examples.
Claude (Fable 5.1) wrote the first draft. I edited heavily and fed the corrections back. That version went to colleagues.
Same process for this talk, in HTML since LLMs are better at that.
Is writing really sacrosanct? Or is tool use OK?

paper title and abstract (click for the draft)

23 / Lessons

Lessons learned

A long-standing idea became the best^* available pulsar search code very quickly, in a language I had never used.
Some of the most fun I’ve had as a scientist in years.
How we do cutting-edge research is changing extremely rapidly, and it isn’t all slop.
Deep domain knowledge (currently) really helps: course corrections kept us out of dead ends and quagmires.
Things I loved doing, like coding, I will probably never do again. Claude et al. are way better at it. (NOTE: I did not really learn Julia)
Probably true for paper and talk writing too (or will be soon).
How do we train our students for this? ANS: I have no idea.
I’m both in awe and a bit frightened.

* based on subjective metrics