Skip to content

Spectra #58

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Spectra #58

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1e2cb44
integrating eventlist and lightcurve in spectras
kashish2210 Jul 28, 2025
ca89280
updating tests and module
kashish2210 Jul 28, 2025
60df7f0
adding test cases
kashish2210 Aug 6, 2025
2902dee
adding condition in posetie_ftt_bins
kashish2210 Aug 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 18 additions & 15 deletions src/Stingray.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,21 +7,6 @@ using DocStringExtensions
using LinearAlgebra
using Random

include("fourier.jl")
export positive_fft_bins
export poisson_level
export normalize_abs
export normalize_frac
export normalize_leahy_from_variance
export normalize_periodograms
export bias_term
export raw_coherence
export estimate_intrinsic_coherence
export error_on_averaged_cross_spectrum
export get_average_ctrate
export get_flux_iterable_from_segments
export avg_pds_from_events
export avg_cs_from_events

include("events.jl")
export FITSMetadata,
Expand Down Expand Up @@ -59,6 +44,24 @@ export AbstractLightCurve,
extract_metadata,
create_lightcurve,
rebin

include("fourier.jl")
export positive_fft_bins
export poisson_level
export normalize_abs
export normalize_frac
export normalize_leahy_from_variance
export normalize_periodograms
export bias_term
export raw_coherence
export estimate_intrinsic_coherence
export error_on_averaged_cross_spectrum
export get_average_ctrate
export get_flux_iterable_from_segments
export avg_pds_from_events,avg_pds_from_iterable
export avg_cs_from_events,avg_cs_from_iterables,avg_cs_from_iterables_quick
export avg_pds_from_eventlist,avg_cs_from_eventlists,avg_pds_from_lightcurve,avg_cs_from_lightcurves

include("utils.jl")

include("gti.jl")
Expand Down
120 changes: 110 additions & 10 deletions src/fourier.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,12 @@ function positive_fft_bins(n_bin::Integer; include_zero::Bool = false)
if include_zero
minbin = 1
end
return (minbin : (n_bin+1) ÷ 2)
if n_bin % 2 == 0
return minbin:(n_bin ÷ 2)
else
return minbin:((n_bin + 1) ÷ 2)
end
end

function poisson_level(norm::String; meanrate = nothing, n_ph = nothing, backrate::Real = 0.0)
if norm == "abs"
return 2.0 * meanrate
Expand Down Expand Up @@ -333,7 +336,7 @@ function avg_pds_from_iterable(flux_iterable, dt::Real; norm::String="frac",
results[!,"freq"] = freq
results[!,"power"] = cross
results[!,"unnorm_power"] = unnorm_cross
results = (; results , n = n_bin, m = n_ave, dt, norm, df = 1 / (dt * n_bin), nphots = n_ph, setment_size = dt * n_bin, mean = common_mean, variance = common_variance)
results = (; results , n = n_bin, m = n_ave, dt, norm, df = 1 / (dt * n_bin), nphots = n_ph, segment_size = dt * n_bin, mean = common_mean, variance = common_variance)
results = DataFrame(results.results)
return results

Expand Down Expand Up @@ -674,17 +677,115 @@ function avg_cs_from_iterables(


return results
end

function avg_pds_from_eventlist(eventlist::EventList, segment_size::Real, dt::Real;
norm::String="frac", use_common_mean::Bool=true,
silent::Bool=false, fluxes=nothing, errors=nothing)
# Extract times and GTI from EventList
times = eventlist.times
gti = eventlist.meta.gti

# If no GTI available, create one from the full time range
if isnothing(gti)
gti = reshape([minimum(times), maximum(times)], 1, 2)
end

# Call the original function with extracted data
return avg_pds_from_events(times, gti, segment_size, dt;
norm=norm, use_common_mean=use_common_mean,
silent=silent, fluxes=fluxes, errors=errors)
end

function avg_cs_from_eventlists(eventlist1::EventList, eventlist2::EventList,
segment_size::Real, dt::Real; norm::String="frac",
use_common_mean::Bool=true, fullspec::Bool=false,
silent::Bool=false, power_type::String="all",
fluxes1=nothing, fluxes2=nothing, errors1=nothing,
errors2=nothing, return_auxil=false)
# Extract times and GTI from EventLists
times1 = eventlist1.times
times2 = eventlist2.times

# Use GTI from first eventlist, or create from time range if not available
gti = eventlist1.meta.gti
if isnothing(gti)
min_time = min(minimum(times1), minimum(times2))
max_time = max(maximum(times1), maximum(times2))
gti = reshape([min_time, max_time], 1, 2)
end

# Call the original function with extracted data
return avg_cs_from_events(times1, times2, gti, segment_size, dt;
norm=norm, use_common_mean=use_common_mean,
fullspec=fullspec, silent=silent,
power_type=power_type, fluxes1=fluxes1,
fluxes2=fluxes2, errors1=errors1,
errors2=errors2, return_auxil=return_auxil)
end

function avg_pds_from_lightcurve(lc::LightCurve, segment_size::Union{Real,Nothing}=nothing;
norm::String="frac", use_common_mean::Bool=true,
silent::Bool=false)
# Extract relevant data from LightCurve
times = lc.time
dt_lc = isa(lc.dt, Vector) ? lc.dt[1] : lc.dt # Use first dt if vector
counts = lc.counts

# Use segment_size or entire light curve
if isnothing(segment_size)
segment_size = (maximum(times) - minimum(times)) + dt_lc
end

# Create flux iterable from light curve data
flux_iterable = [counts] # Simple case: use entire light curve as one segment

return avg_pds_from_iterable(flux_iterable, dt_lc; norm=norm,
use_common_mean=use_common_mean, silent=silent)
end

function avg_cs_from_lightcurves(lc1::LightCurve, lc2::LightCurve,
segment_size::Union{Real,Nothing}=nothing;
norm::String="frac", use_common_mean::Bool=true,
fullspec::Bool=false, silent::Bool=false,
power_type::String="all", return_auxil=false)
# Extract relevant data from LightCurves
dt1 = isa(lc1.dt, Vector) ? lc1.dt[1] : lc1.dt
dt2 = isa(lc2.dt, Vector) ? lc2.dt[1] : lc2.dt

# Ensure both light curves have the same time resolution
@assert dt1 ≈ dt2 "Light curves must have the same time resolution"

counts1 = lc1.counts
counts2 = lc2.counts

# Ensure both light curves have the same length
min_length = min(length(counts1), length(counts2))
if length(counts1) != length(counts2)
@warn "Light curves have different lengths, truncating to shorter one"
counts1 = counts1[1:min_length]
counts2 = counts2[1:min_length]
end

# Create flux iterables from light curve data
flux_iterable1 = [counts1] # Simple case: use entire light curves as one segment
flux_iterable2 = [counts2]

return avg_cs_from_iterables(flux_iterable1, flux_iterable2, dt1;
norm=norm, use_common_mean=use_common_mean,
fullspec=fullspec, silent=silent,
power_type=power_type, return_auxil=return_auxil)
end

function avg_pds_from_events(times:: AbstractVector{<:Real}, gti::AbstractMatrix{<:Real},
# Original functions (kept for backward compatibility)
function avg_pds_from_events(times::AbstractVector{<:Real}, gti::AbstractMatrix{<:Real},
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this functions

segment_size::Real, dt::Real; norm::String="frac",
use_common_mean::Bool=true, silent::Bool=false,
fluxes=nothing, errors=nothing)
if isnothing(segment_size)
segment_size = max(gti) - min(gti)
segment_size = maximum(gti) - minimum(gti)
end
n_bin = round(Int,segment_size / dt)
n_bin = round(Int, segment_size / dt)
dt = segment_size / n_bin

flux_iterable = get_flux_iterable_from_segments(times, gti, segment_size;
Expand All @@ -703,14 +804,14 @@ function avg_pds_from_events(times:: AbstractVector{<:Real}, gti::AbstractMatrix

end

function avg_cs_from_events(times1:: AbstractVector{<:Real}, times2:: AbstractVector{<:Real},
function avg_cs_from_events(times1::AbstractVector{<:Real}, times2::AbstractVector{<:Real},
gti::AbstractMatrix{<:Real}, segment_size::Real, dt::Real;
norm::String="frac", use_common_mean::Bool=true,
fullspec::Bool=false, silent::Bool=false,
power_type::String="all", fluxes1=nothing, fluxes2=nothing,
errors1=nothing, errors2=nothing, return_auxil=false)
if isnothing(segment_size)
segment_size = max(gti) - min(gti)
segment_size = maximum(gti) - minimum(gti)
end
n_bin = round(Int, segment_size / dt)
# adjust dt
Expand All @@ -723,8 +824,7 @@ function avg_cs_from_events(times1:: AbstractVector{<:Real}, times2:: AbstractVe
times2, gti, segment_size; n_bin, fluxes=fluxes2, errors=errors2
)

is_events = all(isnothing,(fluxes1, fluxes2, errors1,
errors2))
is_events = all(isnothing, (fluxes1, fluxes2, errors1, errors2))

if (is_events
&& silent
Expand Down
6 changes: 5 additions & 1 deletion test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,11 @@ using Logging ,LinearAlgebra
using CFITSIO
using Random

include("test_fourier.jl")
@testset "Fourier" begin
include("test_fourier/test_fourier.jl")
include("test_fourier/test_coherence.jl")
include("test_fourier/test_norm.jl")
end
@testset "GTI" begin
include("test_gti.jl")
end
Expand Down
Loading