cross spectrum and ploting recpies

[kashish2210]

this pr consists of cross-spectrum, timing analysis, and plotting recipes:

So the methods included here as :

event1 = readevents("ni1050360108_0mpu7_cl.evt", load_gti=true, sort=true)
ev_lowE  = filter_energy(x -> x >= 20 && x < 500, event1)
ev_highE = filter_energy(x -> x >= 500 && x <= 1500, event1)
println("Low energy events: ", length(ev_lowE))
println("High energy events: ", length(ev_highE))

Low energy events: 28841910
High energy events: 2941370

lc1 = create_lightcurve(ev_lowE,0.01)
lc2 = create_lightcurve(ev_highE,0.01)
cs = CrossSpectrum(ev_lowE, ev_highE, 64.0, 0.01; norm="frac")
acs = AveragedCrossSpectrum(ev_lowE, ev_highE, 64.0, 0.01)

CrossSpectrum{Float64} (Single)
  Frequencies: 3199
  Normalization: frac
CrossSpectrum{Float64} (Averaged)
  Frequencies: 3199
  Segments averaged: 53
  Segment size: 64.0
  Mean rates: 8205.026827830188, 828.4949882075472
  Normalization: frac

i am using these functions to determine whether it's an averaged or a single cross-spectrum

is_averaged(cs::CrossSpectrum) = cs.m > 1
is_single(cs::CrossSpectrum) = cs.m == 1

now plotting:

simple plot direct cs and acs

plot(cs)
plot(acs)

or

plot(acs, plot_type=:amplitude)
plot(cs, plot_type=:amplitude)

u can add a frequency range to:

plot(cs, plot_type=:amplitude, freq_range=(0.1, 10.0), color=:green)

Cross power spectrum

plot(acs, plot_type=:power, color=:purple)
plot(cs, plot_type=:power, color=:purple)

Phase lag spectrum

plot(cs, plot_type=:phase_lag, color=:red)
plot(acs, plot_type=:phase_lag, color=:red)

Time lag spectrum

plot(cs, plot_type=:time_lag, color=:orange)
plot(acs, plot_type=:time_lag, color=:orange)

Coherence function

plot(cs, plot_type=:coherence, color=:green)
plot(acs, plot_type=:coherence, color=:green)

Real and imaginary parts

plot(cs, plot_type=:real_imaginary)
plot(acs, plot_type=:real_imaginary)

PDS comparison

plot(cs, plot_type=:pds_comparison)
plot(acs, plot_type=:pds_comparison)

You can have a comprehensive 6-panel analysis of all plots via:

plot(cs, Val(:analysis))
plot(acs, Val(:analysis))

Coherence with confidence levels

plot(cs, Val(:coherence_confidence))
plot(acs, Val(:coherence_confidence))

Time lag with error bars and reference values

plot(cs, Val(:lag_frequency_errors), 
     freq_range=(1.0, 5.0),
     reference_lag=Dict(:frequency => 3.0, :expected_lag => 0.15))
plot(acs, Val(:lag_frequency_errors), 
     freq_range=(1.0, 5.0),
     reference_lag=Dict(:frequency => 3.0, :expected_lag => 0.15))

Phase lag with error bars and reference values

plot(cs, Val(:phase_frequency_errors),
     freq_range=(2.0, 5.0), 
     reference_phase=Dict(:frequency => 3.0, :expected_phase => π/3))
plot(acs, Val(:phase_frequency_errors),
     freq_range=(2.0, 5.0), 
     reference_phase=Dict(:frequency => 3.0, :expected_phase => π/3))

Time lag with rebinning

plot(cs, Val(:lag_frequency_errors), rebin_factor=4, freq_range=(1.0, 10.0))
plot(acs, Val(:lag_frequency_errors), rebin_factor=4, freq_range=(1.0, 10.0))

Compare two light curves

plot(lc1, lc2)

#using just plots, not recipes 
dt = 0.01
test_times = 0:dt:100  # 100 seconds
test_signal1 = 300 * sin.(2π * test_times / 0.5) .+ 1000
test_signal2 = 200 * sin.(2π * test_times / 0.5 .+ π/4) .+ 900

# Simple plot test
plot(test_times, test_signal1, label="Signal 1", linewidth=2)
plot!(test_times, test_signal2, label="Signal 2", linewidth=2)

Normalization comparison (requires an array of differently normalized CrossSpectra)

cs_leahy = CrossSpectrum(lc1, lc2, norm="leahy")
cs_frac = CrossSpectrum(lc1, lc2, norm="frac") 
cs_abs = CrossSpectrum(lc1, lc2, norm="abs")

plot([cs_leahy, cs_frac, cs_abs], Val(:normalization_comparison))

plot([cs_leahy, cs_frac, cs_abs], Val(:normalization_comparison),
     norm_labels=["Leahy Norm", "Fractional RMS", "Absolute RMS"])

acs_leahy = AveragedCrossSpectrum(lc1, lc2,62, norm="leahy")
acs_frac = AveragedCrossSpectrum(lc1, lc2,62, norm="frac") 
acs_abs = AveragedCrossSpectrum(lc1, lc2,62, norm="abs")

plot([acs_leahy, acs_frac, acs_abs], Val(:normalization_comparison))

plot([acs_leahy, acs_frac, acs_abs], Val(:normalization_comparison),
     norm_labels=["Leahy Norm", "Fractional RMS", "Absolute RMS"])

Rebinning comparison (original vs rebinned)

cs_rebinned = rebin(cs, 0.25)  # Your rebinning function
plot(cs, cs_rebinned, Val(:rebinning_comparison), rebin_type="Linear")
acs_rebinned = rebin(acs, 0.25)  # Your rebinning function
plot(acs, acs_rebinned, Val(:rebinning_comparison), rebin_type="Linear")

Log rebinning comparison

cs_log_rebinned = rebin_log(cs, f=0.02)
plot(cs, cs_log_rebinned, Val(:rebinning_comparison), rebin_type="Logarithmic")
acs_log_rebinned = rebin_log(acs, f=0.02)   
plot(acs, acs_log_rebinned, Val(:rebinning_comparison), rebin_type="Logarithmic")

geometric rebinning comparison

cs_geom_rebinned = geometric_rebin(cs, 1.3)
plot(cs, cs_geom_rebinned, Val(:rebinning_comparison), rebin_type="Geometric")
acs_geom_rebinned = geometric_rebin(acs, 1.3)
plot(acs, acs_geom_rebinned, Val(:rebinning_comparison), rebin_type="Geometric")

integer rebinning comparison

cs_int_rebinned = rebin(cs, 4)  # Rebin by factor of 4
plot(cs, cs_int_rebinned, Val(:rebinning_comparison), rebin_type="Integer")
acs_int_rebinned = rebin(acs, 4)  # Rebin by factor of 4
plot(acs, acs_int_rebinned, Val(:rebinning_comparison), rebin_type="Integer")

noise analysis

plot(cs, Val(:noise_analysis))
plot(acs, Val(:noise_analysis))

signal to noise

plot(cs, Val(:frequency_snr), rebin_factor=5)
plot(acs, Val(:frequency_snr), rebin_factor=5)

noise compare

plot(cs, Val(:noise_comparison))
plot(acs, Val(:noise_comparison))

white noise

plot(cs, Val(:white_noise), high_freq_fraction=0.3)
plot(acs, Val(:white_noise), high_freq_fraction=0.3)

output: same as amplitude :)

phase lag errors

plot(cs, Val(:phase_lag_errors))
plot(acs, Val(:phase_lag_errors))

noise timeline (I created for study purposes)

plot(cs, Val(:noise_timeline))
plot(acs, Val(:noise_timeline))

noise properties

plot(acs, Val(:noise_properties))
plot(cs, Val(:noise_properties))

significant_detections

plot(acs, Val(:significant_detections), threshold=3.0)
plot(cs, Val(:significant_detections), threshold=3.0)

tagging @matteobachetti , @fjebaker , @stefanocovino

I know I have made a huge pr (apology for this), but this pr covers full cross-spectrum analysis and docs, and tests, that's why much bigger

These functions are functions I really liked since they are automatic in nature : )

[codecov-commenter]

Codecov Report

❌ Patch coverage is 58.55513% with 654 lines in your changes missing coverage. Please review.
✅ Project coverage is 70.53%. Comparing base (2a7829d) to head (9507769).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/plotting/plots_recipes_crossspectrum.jl	48.13%	473 Missing ⚠️
src/crossspectrum.jl	74.95%	144 Missing ⚠️
src/fourier.jl	61.90%	24 Missing ⚠️
src/utils.jl	45.83%	13 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##             main      #65       +/-   ##
===========================================
- Coverage   88.06%   70.53%   -17.54%     
===========================================
  Files           5        7        +2     
  Lines        1081     2654     +1573     
===========================================
+ Hits          952     1872      +920     
- Misses        129      782      +653     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.