From Photons to Power Spectra: Journey with Stingray.jl

-

From Photons to Power Spectra: Journey with Stingray.jl

The final week of Google Summer of Code is here. As I sit back and look at the codebase, it feels surreal. Three months ago, I was an undergraduate with a deep fascination for black holes, neutron stars, and the extreme physics of the cosmos. I knew I wanted to write code that would help decode the universe, but I never imagined the ride would be this thrilling.

When I started my journey with Stingray.jl under the OpenAstronomy umbrella, the mission was clear: take the powerful spectral-timing capabilities of the Python Stingray library and bring them to the Julia ecosystem. Python is fantastic, but when you are dealing with millions of X-ray photons and computing averaged cross-spectra over thousands of segments, the blazing speed of Julia becomes a game-changer.

This is the story of how we built a high-performance spectral timing engine, the promises we kept, the fun we had, and where we go from here.

What was Promised vs. What was Delivered

At the start of the summer, my proposal outlined an ambitious set of goals: building high-level spectral types, implementing time lags, porting the Lomb-Scargle periodogram, and laying the foundation for a variability-vs-energy framework.

As we got into the weeds of the Fourier transforms and spectral analysis, here is what we successfully shipped into the core library:

High-Level Spectral Types: We implemented the Powerspectrum and Crossspectrum struct types from scratch, giving users a clean, type-safe API rather than dealing with raw DataFrames.

Mission-Specific FITS Support: We added telescope-specific event file readers for missions like NICER, NuSTAR, XMM-Newton, and RXTE.

Frequency Rebinning: Implemented both geometric (logarithmic) and linear frequency rebinning for proper broadband noise studies.

Robust Normalizations: Added support for Leahy, Fractional RMS, and Absolute RMS normalizations, rigorously tested against Poisson noise statistics.

Lomb-Scargle Periodograms: Ported the Lomb-Scargle algorithm to handle unevenly sampled astrophysical data.

Fourier Utilities: Successfully ported complex utility functions like shift_and_add for dynamic spectral analysis.

Codecov & Testing: Integrated Codecov and dramatically improved the test suite to ensure our Julia port produces scientifically valid results matching its Python sibling.

Code In Action: Power Spectra and Time Lags

To give you a taste of what the new Stingray.jl API feels like, let's look at a real-world workflow. Imagine you have raw photon event data from the NICER telescope and you are hunting for a Quasi-Periodic Oscillation (QPO) from an accreting black hole.

Before GSoC, this would have required manual DataFrame manipulation. Now? It’s beautifully concise.

Generating and Rebinning a Power Spectrum

julia

using Stingray

using CairoMakie

# 1. Load up mission data (e.g., NICER)

ev = readevents("nicer_observation.fits")

# 2. Generate an Averaged Power Spectrum

# (128s segments, Leahy normalization sets Poisson noise at power = 2)

ps = Powerspectrum(ev, segment_size=128.0, dt=1/4096, norm="leahy")

# 3. Logarithmic Rebinning for broadband noise studies

# (Bins get 3% wider at each step, perfect for log-log plots!)

ps_rebinned = rebin(ps, factor=1.03, method=:geometric)

With our plotting utilities and Makie integration, we can easily visualize what a QPO looks like on top of a red-noise continuum:


Review

Simulated QPO Plot

A simulated power spectrum showing a 5 Hz Quasi-Periodic Oscillation (QPO) sitting above the Poisson (white) noise level, correctly normalized using the Leahy formulation. The red points show the geometrically rebinned power.

Computing Time Lags via Cross-Spectra

What if we want to measure the light-travel time delay between soft X-ray photons reflecting off the accretion disk and hard X-ray photons from the corona? We use the Crossspectrum:

julia

# Filter events into soft (0.3-1.0 keV) and hard (2.0-10.0 keV) bands

soft_ev = filter_energy(e -> 0.3 <= e <= 1.0, ev)

hard_ev = filter_energy(e -> 2.0 <= e <= 10.0, ev)

# Compute the averaged cross spectrum

cs = Crossspectrum(soft_ev, hard_ev, segment_size=64.0, dt=1/256)

# Calculate the frequency-dependent time lags (and intrinsic coherence)

lags, lag_err = time_lag(cs)

coh, coh_err = coherence(cs)

This API mimics the flexibility of Python but compiles down to heavily optimized Julia machine code!

Fun Times, Mentorship, and the Open Source Spirit

None of this would have been possible without my mentors, Fergus and Matteo.

Our weekly sync-up meetings were the absolute highlight of my GSoC experience.
Whether we were diving deep into the math behind intrinsic coherence, squashing sneaky bugs related to phase wrapping in complex cross-spectra, or just having a good laugh about timezone differences, those meetings were something I looked forward to every week.

They taught me more than just how to write good Julia code; they taught me how to think like a scientific software engineer.
They pushed me to rigorously test edge cases (like making sure our Leahy normalization perfectly averaged to 2.0 for Poisson noise!) and showed me what it takes to maintain an open-source astrophysics package. Thank you, Fergus and Matteo, for an unforgettable summer!

Beyond GSoC: My Future with Stingray.jl

Google Summer of Code might be ending, but my journey with Stingray is definitely not over. I have fallen in love with this codebase and the OpenAstronomy community.

Moving forward, I plan to stay actively involved as a core contributor. Here are some of the goals I still want to accomplish:

  1. The Variability-vs-Energy Framework: Completing the LagEnergySpectrum and CovarianceSpectrum implementations to allow researchers to easily map out the geometry of accretion disks.

  2. Native Plot Recipes: Finalizing our Makie.jl and Plots.jl recipes so users can generate publication-ready plots with a simple plot(ps) command.

  3. Community Mentorship: I want to help onboard the next wave of students who discover Stingray.jl, just as Fergus and Matteo helped me.

To everyone who followed along with my progress, reviewed my PRs, and supported me, thank you. The universe is vast, but with tools like Julia and Stingray, we are getting a little bit better at understanding it every single day.

Signing off (for now!),




 

Comments

Post a Comment

Popular posts from this blog

Core Spectral Types

-
Image
  Lift off: Weeks 1 & 2 The first month of Google Summer of Code has been nothing short of a whirlwind. Here's a quick breakdown of how the first two weeks went down: Navigating the Codebase: If you've ever stared at a massive codebase and wondered where to begin, you know how my first few days felt. Getting my hands dirty made things click into place rapidly. Finishing Power Colors: My immediate target was getting the Power Colors implementation over the finish line. Classifying spectral states based on power spectral shapes is a huge deal in X-ray binaries. Porting from Python to Julia: Bringing this capability from the established Python ecosystem to our native Julia environment felt like unlocking a new superpower. Rigorous Testing: It required strict testing against established scientific benchmarks, but finally seeing those PR checks turn green was an absolute thrill. Leveling Up: Core Spectral Types: Weeks 3 & 4 With the confidence boost from the Power Colors mer...
Read more

Rebinning Updates!

-
Image
  Updates on Rebinning in Stingray.jl Welcome back! Over the past couple of weeks, we have been working hard on expanding the capabilities of Stingray.jl. One of the key features we've recently added is data rebinning. Rebinning is an essential part of time series analysis, especially when dealing with noisy light curves or power spectra where we want to trade off some resolution for a higher signal-to-noise ratio. In this post, we'll walk through the two new rebinning functions we've introduced: rebin_data for linear rebinning and rebin_data_log for logarithmic rebinning. Linear Rebinning with rebin_data The rebin_data function allows you to linearly rebin your data array to a new, coarser resolution. You can choose to either sum or average the data points falling into the new bins using the method keyword argument. It also automatically propagates uncertainties if you provide them! Here is a quick example of how you can use it to rebin a light curve or a signal:
Read more