Powerspectrum and averagedpowerspectrum

[kashish2210]

this pr contains the powerspectrum and averagepowerspectrum functions for light curve[binned edges] and event[unbinned edges]

Here are some examples with output :

events[unbinned]

events = readevents("ni1200120104_0mpu7_cl.evt", load_gti=true, sort=true)
ps = Powerspectrum(events, norm="leahy")
avg_ps = AveragedPowerspectrum(events, 1024.0, norm="leahy", dt=0.1)

output: data [ik redunddant :)]

lightcurve[binned]

lc = create_lightcurve(events, 0.1)
ps_lc = Powerspectrum(lc, norm="leahy")
ps_avg_lc = AveragedPowerspectrum(lc, 1024.0, norm="leahy")
ps_avg_lc = AveragedPowerspectrum(lc, 1024.0, norm="leahy")

output the same type as the above one:)

println(ps.df)
println(avg_ps.df)
println(ps_lc.df)
println(ps_avg_lc.df)

output:

1.3662135391761732e-5
0.0009765625
1.3662135391761732e-5
0.0009765625

println(ps.m)
println(avg_ps.m)
println(ps_lc.m)
println(ps_avg_lc.m)

output:

1
7
1
7

println(ps.n)
println(avg_ps.n)
println(ps_lc.n)
println(ps_avg_lc.n)

output:

36597
5119
365974
5119

println(ps.nphots)
println(avg_ps.nphots)
println(ps_lc.nphots)
println(ps_avg_lc.nphots)

output

21244574
19285950
21244574
19544935

println(avg_ps.mean_rate)
println(ps_avg_lc.mean_rate)

output

2690.5622209821427
2726.692940848214

# Access frequency and power arrays
freqs = ps.freq
power = ps.power
freqs_lc = ps_lc.freq
power_lc = ps_lc.power

# For averaged spectruōm
avg_freqs = avg_ps.freq  
avg_power = avg_ps.power
avg_freqs_lc = ps_avg_lc.freq  
avg_power_lc = ps_avg_lc.power

# noise reduction
println("Single spectrum noise level: ", std(ps.power[end-100:end]))
println("Averaged spectrum noise level: ", std(avg_ps.power[end-100:end]))
println("Single spectrum noise level[lc]: ", std(ps_lc.power[end-100:end]))
println("Averaged spectrum noise level[lc]: ", std(ps_avg_lc.power[end-100:end]))

output

Single spectrum noise level: 222.6690624733685
Averaged spectrum noise level: 14.04983825672826
Single spectrum noise level[lc]: 11.118772184090005
Averaged spectrum noise level[lc]: 4.547165020012018

tagging @fjebaker @matteobachetti @stefanocovino

[codecov-commenter]

Codecov Report

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

Files with missing lines	Patch %	Lines
src/powerspectrum.jl	85.93%	27 Missing ⚠️
src/lightcurve.jl	80.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #57      +/-   ##
==========================================
- Coverage   88.06%   87.78%   -0.29%     
==========================================
  Files           5        6       +1     
  Lines        1081     1277     +196     
==========================================
+ Hits          952     1121     +169     
- Misses        129      156      +27     

☔ 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.

[kashish2210]

i have tested like :

events = readevents("ni1200120104_0mpu7_cl.evt", load_gti=true, sort=true)

# Powerspectrum(events, dt, segment_size; norm="leahy")
ps_direct = Powerspectrum(events, 0.01, 1000, norm="frac")

# Events method with AveragedPowerspectrum - same as before
ps_events = AveragedPowerspectrum(events, 1000, norm="frac", dt=0.01)

# Compare with light curve method
lc = create_lightcurve(events, 0.01)  # Using same dt=1.0 for fair comparison
ps_lc = AveragedPowerspectrum(lc, 1000, norm="frac")

# Check results
println("Direct method nphots: ", ps_direct.nphots)
println("Events method nphots: ", ps_events.nphots)  
println("LC method nphots: ", ps_lc.nphots)

# All three should now have mean_rate since they're all AveragedPowerspectrum
println("\nMean rate comparison:")
println("Direct method mean_rate: ", ps_direct.mean_rate)
println("Events method mean_rate: ", ps_events.mean_rate)
println("LC method mean_rate: ", ps_lc.mean_rate)

# Additional comparisons
println("\nFrequency resolution comparison:")
println("Direct method df: ", ps_direct.df)
println("Events method df: ", ps_events.df)
println("LC method df: ", ps_lc.df)

println("\nNumber of segments:")
println("Direct method m: ", ps_direct.m)
println("Events method m: ", ps_events.m)
println("LC method m: ", ps_lc.m)

println("\nSegment size:")
println("Direct method segment_size: ", ps_direct.segment_size)
println("Events method segment_size: ", ps_events.segment_size)
println("LC method segment_size: ", ps_lc.segment_size)

output :

Found GTI data: 16 intervals
GTI time range: 1.3253976595089495e8 to 1.3261337476368216e8
Direct method nphots: 19163002
Events method nphots: 19163002
LC method nphots: 19497872

Mean rate comparison:
Direct method mean_rate: 2737.5717142857143
Events method mean_rate: 2737.5717142857143
LC method mean_rate: 2785.4102857142857

Frequency resolution comparison:
Direct method df: 0.001
Events method df: 0.001
LC method df: 0.001

Number of segments:
Direct method m: 7
Events method m: 7
LC method m: 7

Segment size:
Direct method segment_size: 1000.0
Events method segment_size: 1000.0
LC method segment_size: 1000.0

The main problem lies in the number of photons. Using the events method and the light curve method should yield the same result, but they exhibit a ~1.3% difference. I am not sure about the numbers, but I have some doubts, like whether it is valid or not?

tagging @matteobachetti

[kashish2210]

#33

[matteobachetti]

Why this double definition of dt?

[kashish2210]

So our constructor signature was like this

 LightCurve(::Vector{T}, ::Union{Vector{T}, T}, ::Vector{Int64}, ::Union{Nothing, Vector{T}}, 
           ::Union{Nothing, Vector{T}}, ::Array{EventProperty{T}, 1}, ::LightCurveMetadata, ::Symbol)

For exposure and count_error, we have the same struct data i.e:Union{Nothing, Vector{T}}
So, I have used this to create a vector where every element is dt for the previous test cases in #53, which we then concluded to remove. I forgot to remove it; by the way, nice catch! :)

[kashish2210]

I have removed that redundant dt, you may review again