diff --git a/src/Stingray.jl b/src/Stingray.jl
index 8f42a49..49a5a74 100644
--- a/src/Stingray.jl
+++ b/src/Stingray.jl
@@ -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,
@@ -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")
diff --git a/src/fourier.jl b/src/fourier.jl
index 0ede925..3817e96 100644
--- a/src/fourier.jl
+++ b/src/fourier.jl
@@ -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
@@ -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
 
@@ -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}, 
                              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;
@@ -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
@@ -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
diff --git a/test/runtests.jl b/test/runtests.jl
index 12ca46a..f8b8038 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -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
diff --git a/test/test_fourier.jl b/test/test_fourier.jl
deleted file mode 100644
index 38de843..0000000
--- a/test/test_fourier.jl
+++ /dev/null
@@ -1,425 +0,0 @@
-function compare_tables(table1, table2; rtol=0.001, discard = [])
-
-    s_discard = Symbol.(discard)
-    test_result = true
-
-    for key in propertynames(table1)
-        if key in s_discard
-            continue
-        end
-        oe, oc = getproperty(table1,key), getproperty(table1,key)
-        if oe isa Integer || oe isa String
-            if !(oe==oc) 
-                test_result = false
-                break
-            end
-        elseif isnothing(oe)
-            if !(isnothing(oc)) 
-                test_result = false
-                break
-            end
-        else
-            if !(≈(oe,oc,rtol=rtol)) 
-                test_result = false
-                break
-            end
-        end
-    end
-
-
-    data_fields_table1 = filter(field -> !(typeof(field) == Symbol && startswith(string(field), "meta")), names(table1))
-    data_fields_table2 = filter(field -> !(typeof(field) == Symbol && startswith(string(field), "meta")), names(table2))
-
-    
-    
-    for field in names(table1)
-        if field in discard
-            continue
-        end
-        oe, oc = table1[!,field], table2[!,field]
-
-        if !(≈(oe,oc,rtol=rtol)) 
-            test_result = false
-            break
-        end
-    end
-    @test test_result
-end
-
-@testset "positive_fft_bins" begin
-    freq = fftfreq(11)
-    goodbins = positive_fft_bins(11)
-    @test filter(x -> x >0, freq) == freq[goodbins]
-end
-
-@testset "test_coherence" begin
-
-    fid = h5open(joinpath(@__DIR__ ,"data","sample_variable_lc.h5"),"r")
-    data = (fid["dataset1"][1:10000]) .* 1000.0
-    close(fid)
-
-    data1 = rand.(Poisson.(data))
-    data2 = rand.(Poisson.(data))
-    ft1 = fft(data1)
-    ft2 = fft(data2)
-    dt = 0.01
-    N = length(data)
-    mean = Statistics.mean(data)
-    meanrate = mean / dt
-    freq = fftfreq(length(data), 1/dt)
-    good = 0 .< freq .< 0.1
-    keepat!(ft1, good)
-    keepat!(ft2, good)
-    cross = normalize_periodograms(
-        ft1 .* conj(ft2), dt, N, mean_flux = mean, norm="abs", power_type="all")
-    pds1 = normalize_periodograms(
-        ft1 .* conj(ft1), dt, N, mean_flux = mean, norm="abs", power_type="real")
-    pds2 = normalize_periodograms(
-        ft2 .* conj(ft2), dt, N, mean_flux = mean, norm="abs", power_type="real")
-
-    p1noise = poisson_level("abs",meanrate=meanrate)
-    p2noise = poisson_level("abs",meanrate=meanrate)
-
-    @testset "test_intrinsic_coherence" begin
-        coh = estimate_intrinsic_coherence.(
-            cross, pds1, pds2, p1noise, p2noise, N)
-        @test all(x -> isapprox(x, 1; atol=0.001), coh)
-    end
-
-    @testset "test_raw_high_coherence" begin
-        coh = raw_coherence.(cross, pds1, pds2, p1noise, p2noise, N)
-        @test all(x -> isapprox(x, 1; atol=0.001), coh)
-    end
-
-    @testset "test_raw_low_coherence" begin
-        nbins = 2
-        C, P1, P2 = @view(cross[1:nbins]), @view(pds1[1:nbins]), @view(pds2[1:nbins])
-        bsq = bias_term.(P1, P2, p1noise, p2noise, N)
-        # must be lower than bsq!
-        low_coh_cross = @. rand(Normal(bsq^0.5 / 10, bsq^0.5 / 100)) + 0.0im
-        coh = raw_coherence.(low_coh_cross, P1, P2, p1noise, p2noise, N)
-        @test all(iszero,coh)
-    end
-
-    @testset "test_raw_high_bias" begin
-        C = [12986.0 + 8694.0im]
-        P1 = [476156.0]
-        P2 = [482751.0]
-        P1noise = 495955
-        P2noise = 494967
-    
-        coh = raw_coherence.(C, P1, P2, P1noise, P2noise, 499; intrinsic_coherence=1)
-        coh_sngl = raw_coherence(C[1], P1[1], P2[1], P1noise, P2noise, 499; intrinsic_coherence=1)
-        @test coh==real(C .* conj(C))./ (P1 .* P2)
-        @test coh_sngl==real(C .* conj(C))[1] / (P1[1] * P2[1])
-    end
-    
-end
-
-@testset "test_fourier" begin
-    dt = 1
-    len = 100
-    ctrate = 10000
-    N = len ÷ dt
-    dt = len / N
-    times = sort(rand(Uniform(0, len), len * ctrate))
-    gti = [[0 len];;]
-    bins = LinRange(0, len, N + 1)
-    counts = fit(Histogram, times, bins).weights
-    errs = fill!(similar(counts),1) * sqrt(ctrate)
-    bin_times = (@view(bins[1:end-1]) + @view(bins[2:end])) / 2
-    segment_size = 20.0
-    times2 = sort(rand(Uniform(0, len), len * ctrate))
-    counts2 = fit(Histogram, times2, bins).weights
-    errs2 = fill!(similar(counts2),1) * sqrt(ctrate)
-
-    @test get_average_ctrate(times, gti, segment_size) == ctrate
-    @test get_average_ctrate(bin_times, gti, segment_size; counts=counts) == ctrate
-
-    @testset "test_fts_from_segments_cts_and_events_are_equal" begin
-        N = round(Int, segment_size / dt)
-        fts_evts = collect(get_flux_iterable_from_segments(times, gti, segment_size, n_bin=N))
-        fts_cts = collect(get_flux_iterable_from_segments(
-                bin_times, gti, segment_size, fluxes=counts))
-        @test fts_evts == fts_cts
-    end
-
-    @testset "test_error_on_averaged_cross_spectrum_low_nave" for common_ref in [true, false]
-        @test_logs (:warn,r"n_ave is below 30."
-        ) error_on_averaged_cross_spectrum([4 + 1.0im], [2], [4], 29, 2, 2; common_ref=common_ref)
-    end
-
-    @testset "test_avg_pds_bad_input" begin
-        _times = rand(Uniform(0,1000),1)
-        out_ev = avg_pds_from_events(_times, gti, segment_size, dt,silent = true)
-        @test isnothing(out_ev)
-    end
-
-    @testset "test_avg_cs_bad_input" for return_auxil in [true, false]
-        _times1 = rand(Uniform(0,1000),1)
-        _times2 = rand(Uniform(0,1000),1)
-        out_ev = avg_cs_from_events(_times1, _times2, gti,
-                                    segment_size, dt, silent = true, return_auxil=return_auxil)
-        @test isnothing(out_ev) 
-    end
-
-    @testset "test_avg_pds_use_common_mean_similar_stats" for norm in ["frac", "abs", "none", "leahy"]
-        out_comm = avg_pds_from_events(
-            times,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=true,
-            silent=true,
-            fluxes=nothing,
-        ).power
-        out = avg_pds_from_events(
-            times,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=false,
-            silent=true,
-            fluxes=nothing,
-        ).power
-        @test std(out_comm)≈std(out) rtol=0.1
-    end
-
-    @testset "test_avg_cs_use_common_mean_similar_stats" for 
-        norm in ["frac", "abs", "none", "leahy"], 
-        return_auxil in [true, false], fullspec in [true,false]
-        out_comm = avg_cs_from_events(
-            times,
-            times2,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=true,
-            silent=true,
-            fullspec=fullspec,
-            return_auxil=return_auxil
-        ).power
-        out = avg_cs_from_events(
-            times,
-            times2,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=false,
-            silent=true,
-            fullspec=fullspec,
-            return_auxil=return_auxil
-        ).power
-        @test std(out_comm)≈std(out) rtol=0.1
-    end
-
-    @testset "test_avg_pds_cts_and_events_are_equal" for use_common_mean in [true,false], norm in ["frac", "abs", "none", "leahy"]
-        out_ev = avg_pds_from_events(
-            times,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=use_common_mean,
-            silent=true,
-            fluxes=nothing,
-        )
-        out_ct = avg_pds_from_events(
-            bin_times,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=use_common_mean,
-            silent=true,
-            fluxes=counts,
-        )
-        compare_tables(out_ev, out_ct)
-    end
-
-    @testset "test_avg_pds_cts_and_err_and_events_are_equal" for use_common_mean in [true,false], norm in ["frac", "abs", "none", "leahy"]       
-        out_ev = avg_pds_from_events(
-            times,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=use_common_mean,
-            silent=true,
-            fluxes=nothing,
-        )
-        out_ct = avg_pds_from_events(
-            bin_times,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=use_common_mean,
-            silent=true,
-            fluxes=counts,
-            errors=errs,
-        )
-        # The variance is not _supposed_ to be equal, when we specify errors
-        if use_common_mean
-            compare_tables(out_ev, out_ct, rtol=0.01, discard=["variance"])
-        else
-            compare_tables(out_ev, out_ct, rtol=0.1, discard=["variance"])
-        end
-    end
-
-    @testset "test_avg_cs_cts_and_events_are_equal" for use_common_mean in [true,false], norm in ["frac", "abs", "none", "leahy"]
-        out_ev = avg_cs_from_events(
-            times,
-            times2,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=use_common_mean,
-            silent=true,
-        )
-        out_ct = avg_cs_from_events(
-            bin_times,
-            bin_times,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=use_common_mean,
-            silent=true,
-            fluxes1=counts,
-            fluxes2=counts2,
-        )
-        if use_common_mean
-            compare_tables(out_ev, out_ct, rtol=0.01)
-        else
-            compare_tables(out_ev, out_ct, rtol=0.1)
-        end
-    end
-
-    @testset "test_avg_cs_cts_and_err_and_events_are_equal" for use_common_mean in [true,false], norm in ["frac", "abs", "none", "leahy"]
-        out_ev = avg_cs_from_events(
-            times,
-            times2,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=use_common_mean,
-            silent=true,
-        )
-        out_ct = avg_cs_from_events(
-            bin_times,
-            bin_times,
-            gti,
-            segment_size,
-            dt,
-            norm=norm,
-            use_common_mean=use_common_mean,
-            silent=true,
-            fluxes1=counts,
-            fluxes2=counts2,
-            errors1=errs,
-            errors2=errs2,
-        )
-        discard = [m for m in propertynames(out_ev) if m == :variance]
-
-        if use_common_mean
-            compare_tables(out_ev, out_ct, rtol=0.01, discard=discard)
-        else
-            compare_tables(out_ev, out_ct, rtol=0.1, discard=discard)
-        end
-    end
-end
-
-@testset "test_norm" begin
-    mean = var = 100000
-    N = 1000000
-    dt = 0.2
-    meanrate = mean / dt
-    lc = rand(Poisson(mean),N)
-    pds = abs2.(fft(lc))
-    freq = fftfreq(N, dt)
-    good = 2:(N÷2)
-
-    pdsabs = normalize_abs(pds, dt, size(lc,1))
-    pdsfrac = normalize_frac(pds, dt, size(lc,1), mean)
-    pois_abs = poisson_level("abs", meanrate=meanrate)
-    pois_frac = poisson_level("frac", meanrate=meanrate)
-
-    @test Statistics.mean(keepat!(pdsabs,good))≈pois_abs rtol=0.02
-    @test Statistics.mean(keepat!(pdsfrac,good))≈pois_frac rtol=0.02
-
-    mean = var = 100000.0
-    N = 800000
-    dt = 0.2
-    df = 1 / (N * dt)
-    freq = fftfreq(N, dt)
-    good = freq .> 0
-    meanrate = mean / dt
-    lc = rand(Poisson(mean),N)
-    nph = sum(lc)
-    pds = keepat!(abs2.(fft(lc)),good)
-    lc_bksub = lc .- mean
-    pds_bksub = keepat!(abs2.(fft(lc_bksub)),good)
-    lc_renorm = lc / mean
-    pds_renorm = keepat!(abs2.(fft(lc_renorm)),good)
-    lc_renorm_bksub = lc_renorm .- 1
-    pds_renorm_bksub = keepat!(abs2.(fft(lc_renorm_bksub)),good)
-
-    @testset "test_leahy_bksub_var_vs_standard" begin
-        leahyvar = normalize_leahy_from_variance(pds_bksub, Statistics.var(lc_bksub), N)
-        leahy =  2 .* pds ./ sum(lc)
-        ratio = Statistics.mean(leahyvar./leahy)
-        @test ratio≈1 rtol=0.01
-    end
-    @testset "test_abs_bksub" begin
-        ratio = normalize_abs(pds_bksub, dt, N) ./ normalize_abs(pds, dt, N)
-        @test Statistics.mean(ratio)≈1 rtol=0.01
-    end
-    @testset "test_frac_renorm_constant" begin
-        ratio = normalize_frac(pds_renorm, dt, N, 1) ./ normalize_frac(pds, dt, N, mean)
-        @test Statistics.mean(ratio)≈1 rtol=0.01
-    end
-    @testset "test_frac_to_abs_ctratesq" begin
-        ratio = (
-            normalize_frac(pds, dt, N, mean)
-            ./ normalize_abs(pds, dt, N)
-            .* meanrate .^ 2
-        )
-        @test Statistics.mean(ratio)≈1 rtol=0.01
-    end
-    @testset "test_total_variance" begin
-        vdk_total_variance = sum((lc .- mean) .^ 2)
-        ratio = Statistics.mean(pds) / vdk_total_variance
-        @test Statistics.mean(ratio)≈1 rtol=0.01
-    end
-    @testset "test_poisson_level_$(norm)" for norm in ["abs", "frac", "leahy","none"]
-        pdsnorm = normalize_periodograms(pds, dt, N; mean_flux=mean, n_ph=nph, norm=norm)
-        @test Statistics.mean(pdsnorm)≈poisson_level(norm; meanrate=meanrate, n_ph=nph) rtol=0.01
-    end
-
-    @testset "test_poisson_level_real_$(norm)" for norm in ["abs", "frac", "leahy","none"]
-        pdsnorm = normalize_periodograms(pds, dt, N; mean_flux=mean, n_ph=nph, norm=norm, power_type = "real")
-        @test Statistics.mean(pdsnorm)≈poisson_level(norm; meanrate=meanrate, n_ph=nph) rtol=0.01
-    end
-
-    @testset "test_poisson_level_absolute_$(norm)" for norm in ["abs", "frac", "leahy","none"]
-        pdsnorm = normalize_periodograms(pds, dt, N; mean_flux=mean, n_ph=nph, norm=norm, power_type = "abs")
-        @test Statistics.mean(pdsnorm)≈poisson_level(norm; meanrate=meanrate, n_ph=nph) rtol=0.01
-    end
-
-    @testset "test_normalize_with_variance" begin
-        pdsnorm = normalize_periodograms(pds, dt, N; mean_flux=mean, variance = var, norm="leahy")
-        @test Statistics.mean(pdsnorm)≈2 rtol=0.01
-    end
-
-    @testset "test_normalize_none" begin
-        pdsnorm = normalize_periodograms(pds, dt, N; mean_flux=mean, n_ph=nph, norm="none")
-        @test Statistics.mean(pdsnorm)≈Statistics.mean(pds) rtol=0.01
-    end
-end
diff --git a/test/test_fourier/test_coherence.jl b/test/test_fourier/test_coherence.jl
new file mode 100644
index 0000000..456b5aa
--- /dev/null
+++ b/test/test_fourier/test_coherence.jl
@@ -0,0 +1,88 @@
+function test_data()
+    fid = h5open(joinpath(@__DIR__, "..", "data", "sample_variable_lc.h5"), "r")
+    data = (fid["dataset1"][1:10000]) .* 1000.0
+    close(fid)
+    dt = 0.01
+    N = length(data)
+    mean = Statistics.mean(data)
+    meanrate = mean / dt
+    freq = fftfreq(length(data), 1/dt)
+    good = 0 .< freq .< 0.1
+    return (data=data, dt=dt, N=N, mean=mean, meanrate=meanrate, good=good)
+end
+
+# Helper function to generate processed data for tests
+function generate_test_data(data; same_poisson=true)
+    data1 = rand.(Poisson.(data.data))
+    data2 = same_poisson ? data1 : rand.(Poisson.(data.data))
+    
+    ft1 = fft(data1)
+    ft2 = fft(data2)
+    keepat!(ft1, data.good)
+    keepat!(ft2, data.good)
+    cross = normalize_periodograms(
+        ft1 .* conj(ft2), data.dt, data.N, 
+        mean_flux = data.mean, norm="abs", power_type="all")
+    pds1 = normalize_periodograms(
+        ft1 .* conj(ft1), data.dt, data.N, 
+        mean_flux = data.mean, norm="abs", power_type="real")
+    pds2 = normalize_periodograms(
+        ft2 .* conj(ft2), data.dt, data.N, 
+        mean_flux = data.mean, norm="abs", power_type="real")
+    p1noise = poisson_level("abs", meanrate=data.meanrate)
+    p2noise = poisson_level("abs", meanrate=data.meanrate)
+    
+    return (cross=cross, pds1=pds1, pds2=pds2, p1noise=p1noise, p2noise=p2noise)
+end
+data = test_data()
+
+# test_coherence
+let
+    test_data = generate_test_data(data, same_poisson=false)
+end
+
+# test_intrinsic_coherence
+let
+    test_data = generate_test_data(data, same_poisson=true)
+    
+    coh = estimate_intrinsic_coherence.(
+        test_data.cross, test_data.pds1, test_data.pds2, 
+        test_data.p1noise, test_data.p2noise, data.N)
+    @test all(x -> isapprox(x, 1; atol=0.001), coh)
+end
+
+# test_raw_high_coherence  
+let
+    test_data = generate_test_data(data, same_poisson=true)
+    
+    coh = raw_coherence.(test_data.cross, test_data.pds1, test_data.pds2, 
+                        test_data.p1noise, test_data.p2noise, data.N)
+    @test all(x -> isapprox(x, 1; atol=0.001), coh)
+end
+
+# test_raw_low_coherence
+let
+    test_data = generate_test_data(data, same_poisson=true)
+    
+    nbins = 2
+    C, P1, P2 = @view(test_data.cross[1:nbins]), @view(test_data.pds1[1:nbins]), @view(test_data.pds2[1:nbins])
+    bsq = bias_term.(P1, P2, test_data.p1noise, test_data.p2noise, data.N)
+    # must be lower than bsq!
+    low_coh_cross = @. rand(Normal(bsq^0.5 / 10, bsq^0.5 / 100)) + 0.0im
+    coh = raw_coherence.(low_coh_cross, P1, P2, test_data.p1noise, test_data.p2noise, data.N)
+    @test all(iszero, coh)
+end
+
+# test_raw_high_bias (this one was already different, so keeping it separate)
+let
+    C = [12986.0 + 8694.0im]
+    P1 = [476156.0]
+    P2 = [482751.0]
+    P1noise = 495955
+    P2noise = 494967
+
+    coh = raw_coherence.(C, P1, P2, P1noise, P2noise, 499; intrinsic_coherence=1)
+    coh_sngl = raw_coherence(C[1], P1[1], P2[1], P1noise, P2noise, 499; intrinsic_coherence=1)
+    @test coh == real(C .* conj(C)) ./ (P1 .* P2)
+    @test coh_sngl == real(C .* conj(C))[1] / (P1[1] * P2[1])
+end
\ No newline at end of file
diff --git a/test/test_fourier/test_fourier.jl b/test/test_fourier/test_fourier.jl
new file mode 100644
index 0000000..a4f53bc
--- /dev/null
+++ b/test/test_fourier/test_fourier.jl
@@ -0,0 +1,495 @@
+function compare_tables(table1, table2; rtol=0.001, discard=[])
+    s_discard = Symbol.(discard)
+    test_result = true
+
+    for key in propertynames(table1)
+        if key in s_discard
+            continue
+        end
+        oe, oc = getproperty(table1, key), getproperty(table2, key)
+        if oe isa Integer || oe isa String
+            if !(oe == oc) 
+                test_result = false
+                break
+            end
+        elseif isnothing(oe)
+            if !(isnothing(oc)) 
+                test_result = false
+                break
+            end
+        else
+            if !(≈(oe, oc, rtol=rtol)) 
+                test_result = false
+                break
+            end
+        end
+    end
+
+    data_fields_table1 = filter(field -> !(typeof(field) == Symbol && startswith(string(field), "meta")), names(table1))
+    data_fields_table2 = filter(field -> !(typeof(field) == Symbol && startswith(string(field), "meta")), names(table2))
+
+    for field in names(table1)
+        if field in discard
+            continue
+        end
+        oe, oc = table1[!, field], table2[!, field]
+
+        if !(≈(oe, oc, rtol=rtol)) 
+            test_result = false
+            break
+        end
+    end
+    @test test_result
+end
+
+const FOURIER_DT = 1
+const FOURIER_LEN = 100
+const FOURIER_CTRATE = 10000
+const FOURIER_SEGMENT_SIZE = 20.0
+
+function setup_fourier_test_data(; dt=FOURIER_DT, len=FOURIER_LEN, ctrate=FOURIER_CTRATE, two_datasets=false)
+    N = len ÷ dt
+    dt = len / N
+    times = sort(rand(Uniform(0, len), len * ctrate))
+    gti = Float64[0 len]
+    bins = LinRange(0, len, N + 1)
+    counts = fit(Histogram, times, bins).weights
+    errs = fill!(similar(counts), 1) * sqrt(ctrate)
+    bin_times = (@view(bins[1:end-1]) + @view(bins[2:end])) / 2
+    
+    if two_datasets
+        times2 = sort(rand(Uniform(0, len), len * ctrate))
+        counts2 = fit(Histogram, times2, bins).weights
+        errs2 = fill!(similar(counts2), 1) * sqrt(ctrate)
+        return (times=times, times2=times2, gti=gti, bins=bins, counts=counts, counts2=counts2, 
+                errs=errs, errs2=errs2, bin_times=bin_times, dt=dt, N=N)
+    else
+        return (times=times, gti=gti, bins=bins, counts=counts, errs=errs, 
+                bin_times=bin_times, dt=dt, N=N)
+    end
+end
+
+function test_pds_common_mean(norm::String)
+    data = setup_fourier_test_data()
+    
+    out_comm = avg_pds_from_events(
+        data.times, data.gti, FOURIER_SEGMENT_SIZE, data.dt,
+        norm=norm, use_common_mean=true, silent=true, fluxes=nothing
+    ).power
+    out = avg_pds_from_events(
+        data.times, data.gti, FOURIER_SEGMENT_SIZE, data.dt,
+        norm=norm, use_common_mean=false, silent=true, fluxes=nothing
+    ).power
+    @test std(out_comm) ≈ std(out) rtol=0.1
+end
+
+function test_cs_common_mean(norm::String, return_auxil::Bool, fullspec::Bool)
+    data = setup_fourier_test_data(two_datasets=true)
+    
+    out_comm = avg_cs_from_events(
+        data.times, data.times2, data.gti, FOURIER_SEGMENT_SIZE, data.dt,
+        norm=norm, use_common_mean=true, silent=true,
+        fullspec=fullspec, return_auxil=return_auxil
+    ).power
+    out = avg_cs_from_events(
+        data.times, data.times2, data.gti, FOURIER_SEGMENT_SIZE, data.dt,
+        norm=norm, use_common_mean=false, silent=true,
+        fullspec=fullspec, return_auxil=return_auxil
+    ).power
+    @test std(out_comm) ≈ std(out) rtol=0.1
+end
+
+function test_events_vs_counts_pds(use_common_mean::Bool, norm::String, include_errors::Bool=false)
+    data = setup_fourier_test_data()
+    
+    out_ev = avg_pds_from_events(
+        data.times, data.gti, FOURIER_SEGMENT_SIZE, data.dt,
+        norm=norm, use_common_mean=use_common_mean, silent=true, fluxes=nothing
+    )
+    
+    if include_errors
+        out_ct = avg_pds_from_events(
+            data.bin_times, data.gti, FOURIER_SEGMENT_SIZE, data.dt,
+            norm=norm, use_common_mean=use_common_mean, silent=true,
+            fluxes=data.counts, errors=data.errs
+        )
+        rtol = use_common_mean ? 0.01 : 0.1
+        compare_tables(out_ev, out_ct, rtol=rtol, discard=["variance"])
+    else
+        out_ct = avg_pds_from_events(
+            data.bin_times, data.gti, FOURIER_SEGMENT_SIZE, data.dt,
+            norm=norm, use_common_mean=use_common_mean, silent=true, fluxes=data.counts
+        )
+        compare_tables(out_ev, out_ct)
+    end
+end
+
+function test_events_vs_counts_cs(use_common_mean::Bool, norm::String, include_errors::Bool=false)
+    data = setup_fourier_test_data(two_datasets=true)
+    
+    out_ev = avg_cs_from_events(
+        data.times, data.times2, data.gti, FOURIER_SEGMENT_SIZE, data.dt,
+        norm=norm, use_common_mean=use_common_mean, silent=true
+    )
+    
+    if include_errors
+        out_ct = avg_cs_from_events(
+            data.bin_times, data.bin_times, data.gti, FOURIER_SEGMENT_SIZE, data.dt,
+            norm=norm, use_common_mean=use_common_mean, silent=true,
+            fluxes1=data.counts, fluxes2=data.counts2, 
+            errors1=data.errs, errors2=data.errs2
+        )
+        discard = [m for m in propertynames(out_ev) if m == :variance]
+        rtol = use_common_mean ? 0.01 : 0.1
+        compare_tables(out_ev, out_ct, rtol=rtol, discard=discard)
+    else
+        out_ct = avg_cs_from_events(
+            data.bin_times, data.bin_times, data.gti, FOURIER_SEGMENT_SIZE, data.dt,
+            norm=norm, use_common_mean=use_common_mean, silent=true,
+            fluxes1=data.counts, fluxes2=data.counts2
+        )
+        rtol = use_common_mean ? 0.01 : 0.1
+        compare_tables(out_ev, out_ct, rtol=rtol)
+    end
+end
+
+let
+    data = setup_fourier_test_data()
+    @test get_average_ctrate(data.times, data.gti, FOURIER_SEGMENT_SIZE) == FOURIER_CTRATE
+    @test get_average_ctrate(data.bin_times, data.gti, FOURIER_SEGMENT_SIZE; counts=data.counts) == FOURIER_CTRATE
+end
+
+let
+    data = setup_fourier_test_data()
+    N = round(Int, FOURIER_SEGMENT_SIZE / data.dt)
+    fts_evts = collect(get_flux_iterable_from_segments(data.times, data.gti, FOURIER_SEGMENT_SIZE, n_bin=N))
+    fts_cts = collect(get_flux_iterable_from_segments(data.bin_times, data.gti, FOURIER_SEGMENT_SIZE, fluxes=data.counts))
+    @test fts_evts == fts_cts
+end
+
+let
+    common_ref = true
+    @test_logs (:warn, r"n_ave is below 30.") error_on_averaged_cross_spectrum([4 + 1.0im], [2], [4], 29, 2, 2; common_ref=common_ref)
+end
+
+let
+    common_ref = false
+    @test_logs (:warn, r"n_ave is below 30.") error_on_averaged_cross_spectrum([4 + 1.0im], [2], [4], 29, 2, 2; common_ref=common_ref)
+end
+
+let
+    gti = Float64[0 FOURIER_LEN]
+    _times = rand(Uniform(0, 1000), 1)
+    out_ev = avg_pds_from_events(_times, gti, FOURIER_SEGMENT_SIZE, FOURIER_DT, silent=true)
+    @test isnothing(out_ev)
+end
+
+let
+    gti = Float64[0 FOURIER_LEN]
+    return_auxil = true
+    _times1 = rand(Uniform(0, 1000), 1)
+    _times2 = rand(Uniform(0, 1000), 1)
+    out_ev = avg_cs_from_events(_times1, _times2, gti, FOURIER_SEGMENT_SIZE, FOURIER_DT, silent=true, return_auxil=return_auxil)
+    @test isnothing(out_ev) 
+end
+
+let
+    gti = Float64[0 FOURIER_LEN]
+    return_auxil = false
+    _times1 = rand(Uniform(0, 1000), 1)
+    _times2 = rand(Uniform(0, 1000), 1)
+    out_ev = avg_cs_from_events(_times1, _times2, gti, FOURIER_SEGMENT_SIZE, FOURIER_DT, silent=true, return_auxil=return_auxil)
+    @test isnothing(out_ev) 
+end
+
+let; test_pds_common_mean("frac"); end
+let; test_pds_common_mean("abs"); end
+let; test_pds_common_mean("none"); end
+let; test_pds_common_mean("leahy"); end
+
+let; test_cs_common_mean("frac", true, true); end
+let; test_cs_common_mean("frac", true, false); end
+let; test_cs_common_mean("frac", false, true); end
+let; test_cs_common_mean("frac", false, false); end
+
+let; test_events_vs_counts_pds(true, "frac"); end
+let; test_events_vs_counts_pds(false, "frac"); end
+let; test_events_vs_counts_pds(true, "frac", true); end
+
+let; test_events_vs_counts_cs(true, "frac"); end
+let; test_events_vs_counts_cs(true, "frac", true); end
+
+# Test LightCurve with vector dt (should use first element)
+let
+    data = setup_fourier_test_data()
+    eventlist = EventList(data.times, nothing, FITSMetadata(
+        "[test]", 1, nothing, Dict{String,Vector}(), 
+        Dict{String,Any}(), data.gti, nothing
+    ))
+    
+    result_eventlist = avg_pds_from_eventlist(eventlist, FOURIER_SEGMENT_SIZE, data.dt, silent=true)
+    result_times = avg_pds_from_events(data.times, data.gti, FOURIER_SEGMENT_SIZE, data.dt, silent=true)
+    
+    @test result_eventlist.power ≈ result_times.power
+    @test result_eventlist.freq ≈ result_times.freq
+end
+
+# Test EventList PDS function equivalence with different norms
+let
+    data = setup_fourier_test_data()
+    eventlist = EventList(data.times, nothing, FITSMetadata(
+        "[test]", 1, nothing, Dict{String,Vector}(), 
+        Dict{String,Any}(), data.gti, nothing
+    ))
+    
+    for norm in ["frac", "abs", "none", "leahy"]
+        result_eventlist = avg_pds_from_eventlist(eventlist, FOURIER_SEGMENT_SIZE, data.dt, norm=norm, silent=true)
+        result_times = avg_pds_from_events(data.times, data.gti, FOURIER_SEGMENT_SIZE, data.dt, norm=norm, silent=true)
+        
+        @test result_eventlist.power ≈ result_times.power
+        @test result_eventlist.freq ≈ result_times.freq
+    end
+end
+
+# Test EventList PDS function with use_common_mean parameter
+let
+    data = setup_fourier_test_data()
+    eventlist = EventList(data.times, nothing, FITSMetadata(
+        "[test]", 1, nothing, Dict{String,Vector}(), 
+        Dict{String,Any}(), data.gti, nothing
+    ))
+    
+    for use_common_mean in [true, false]
+        result_eventlist = avg_pds_from_eventlist(eventlist, FOURIER_SEGMENT_SIZE, data.dt, use_common_mean=use_common_mean, silent=true)
+        result_times = avg_pds_from_events(data.times, data.gti, FOURIER_SEGMENT_SIZE, data.dt, use_common_mean=use_common_mean, silent=true)
+        
+        @test result_eventlist.power ≈ result_times.power
+        @test result_eventlist.freq ≈ result_times.freq
+    end
+end
+
+# Test EventList cross spectrum function equivalence
+let
+    data = setup_fourier_test_data(two_datasets=true)
+    eventlist1 = EventList(data.times, nothing, FITSMetadata(
+        "[test]", 1, nothing, Dict{String,Vector}(), 
+        Dict{String,Any}(), data.gti, nothing
+    ))
+    eventlist2 = EventList(data.times2, nothing, FITSMetadata(
+        "[test]", 1, nothing, Dict{String,Vector}(), 
+        Dict{String,Any}(), data.gti, nothing
+    ))
+    
+    result_eventlists = avg_cs_from_eventlists(eventlist1, eventlist2, FOURIER_SEGMENT_SIZE, data.dt, silent=true)
+    result_times = avg_cs_from_events(data.times, data.times2, data.gti, FOURIER_SEGMENT_SIZE, data.dt, silent=true)
+    
+    @test result_eventlists.power ≈ result_times.power
+    @test result_eventlists.freq ≈ result_times.freq
+end
+
+# Test EventList cross spectrum function with different parameters
+let
+    data = setup_fourier_test_data(two_datasets=true)
+    eventlist1 = EventList(data.times, nothing, FITSMetadata(
+        "[test]", 1, nothing, Dict{String,Vector}(), 
+        Dict{String,Any}(), data.gti, nothing
+    ))
+    eventlist2 = EventList(data.times2, nothing, FITSMetadata(
+        "[test]", 1, nothing, Dict{String,Vector}(), 
+        Dict{String,Any}(), data.gti, nothing
+    ))
+    
+    for norm in ["frac", "abs", "none", "leahy"]
+        for use_common_mean in [true, false]
+            for fullspec in [true, false]
+                for return_auxil in [true, false]
+                    result_eventlists = avg_cs_from_eventlists(
+                        eventlist1, eventlist2, FOURIER_SEGMENT_SIZE, data.dt,
+                        norm=norm, use_common_mean=use_common_mean,
+                        fullspec=fullspec, return_auxil=return_auxil, silent=true
+                    )
+                    result_times = avg_cs_from_events(
+                        data.times, data.times2, data.gti, FOURIER_SEGMENT_SIZE, data.dt,
+                        norm=norm, use_common_mean=use_common_mean,
+                        fullspec=fullspec, return_auxil=return_auxil, silent=true
+                    )
+                    
+                    @test result_eventlists.power ≈ result_times.power
+                    @test result_eventlists.freq ≈ result_times.freq
+                    
+                    if return_auxil
+                        @test result_eventlists.pds1 ≈ result_times.pds1
+                        @test result_eventlists.pds2 ≈ result_times.pds2
+                    end
+                end
+            end
+        end
+    end
+end
+
+# Test LightCurve PDS function equivalence
+let
+    dt = 1.0
+    N = 100
+    times = collect(0:dt:(N-1)*dt)
+    counts = rand(Poisson(100), N)
+    
+    metadata = LightCurveMetadata(
+        "TEST", "INST", "TEST_OBJ", 58000.0, 
+        (times[1], times[end]), dt,
+        [Dict{String,Any}()], Dict{String,Any}()
+    )
+    
+    lc = LightCurve(times, dt, counts, nothing, nothing, 
+                    EventProperty{Float64}[], metadata, :poisson)
+    flux_iterable = [counts]
+    
+    result_lightcurve = avg_pds_from_lightcurve(lc, silent=true)
+    result_iterable = avg_pds_from_iterable(flux_iterable, dt, silent=true)
+    
+    @test result_lightcurve.power ≈ result_iterable.power
+    @test result_lightcurve.freq ≈ result_iterable.freq
+end
+
+# Test LightCurve PDS function with different norms
+let
+    dt = 1.0
+    N = 100
+    times = collect(0:dt:(N-1)*dt)
+    counts = rand(Poisson(100), N)
+    
+    metadata = LightCurveMetadata(
+        "TEST", "INST", "TEST_OBJ", 58000.0, 
+        (times[1], times[end]), dt,
+        [Dict{String,Any}()], Dict{String,Any}()
+    )
+    
+    lc = LightCurve(times, dt, counts, nothing, nothing, 
+                    EventProperty{Float64}[], metadata, :poisson)
+    flux_iterable = [counts]
+    
+    for norm in ["frac", "abs", "none", "leahy"]
+        for use_common_mean in [true, false]
+            result_lightcurve = avg_pds_from_lightcurve(lc, norm=norm, use_common_mean=use_common_mean, silent=true)
+            result_iterable = avg_pds_from_iterable(flux_iterable, dt, norm=norm, use_common_mean=use_common_mean, silent=true)
+            
+            @test result_lightcurve.power ≈ result_iterable.power
+            @test result_lightcurve.freq ≈ result_iterable.freq
+        end
+    end
+end
+
+# Test LightCurve cross spectrum function equivalence
+let
+    dt = 1.0
+    N = 100
+    times = collect(0:dt:(N-1)*dt)
+    counts1 = rand(Poisson(100), N)
+    counts2 = rand(Poisson(100), N)
+    
+    metadata = LightCurveMetadata(
+        "TEST", "INST", "TEST_OBJ", 58000.0, 
+        (times[1], times[end]), dt,
+        [Dict{String,Any}()], Dict{String,Any}()
+    )
+    
+    lc1 = LightCurve(times, dt, counts1, nothing, nothing, 
+                     EventProperty{Float64}[], metadata, :poisson)
+    lc2 = LightCurve(times, dt, counts2, nothing, nothing, 
+                     EventProperty{Float64}[], metadata, :poisson)
+    
+    flux_iterable1 = [counts1]
+    flux_iterable2 = [counts2]
+    
+    result_lightcurves = avg_cs_from_lightcurves(lc1, lc2, silent=true)
+    result_iterables = avg_cs_from_iterables(flux_iterable1, flux_iterable2, dt, silent=true)
+    
+    @test result_lightcurves.power ≈ result_iterables.power
+    @test result_lightcurves.freq ≈ result_iterables.freq
+end
+
+# Test LightCurve cross spectrum function with different parameters
+let
+    dt = 1.0
+    N = 100
+    times = collect(0:dt:(N-1)*dt)
+    counts1 = rand(Poisson(100), N)
+    counts2 = rand(Poisson(100), N)
+    
+    metadata = LightCurveMetadata(
+        "TEST", "INST", "TEST_OBJ", 58000.0, 
+        (times[1], times[end]), dt,
+        [Dict{String,Any}()], Dict{String,Any}()
+    )
+    
+    lc1 = LightCurve(times, dt, counts1, nothing, nothing, 
+                     EventProperty{Float64}[], metadata, :poisson)
+    lc2 = LightCurve(times, dt, counts2, nothing, nothing, 
+                     EventProperty{Float64}[], metadata, :poisson)
+    
+    flux_iterable1 = [counts1]
+    flux_iterable2 = [counts2]
+    
+    for norm in ["frac", "abs", "none", "leahy"]
+        for use_common_mean in [true, false]
+            for fullspec in [true, false]
+                for return_auxil in [true, false]
+                    result_lightcurves = avg_cs_from_lightcurves(
+                        lc1, lc2, norm=norm, use_common_mean=use_common_mean,
+                        fullspec=fullspec, return_auxil=return_auxil, silent=true
+                    )
+                    result_iterables = avg_cs_from_iterables(
+                        flux_iterable1, flux_iterable2, dt,
+                        norm=norm, use_common_mean=use_common_mean,
+                        fullspec=fullspec, return_auxil=return_auxil, silent=true
+                    )
+                    
+                    @test result_lightcurves.power ≈ result_iterables.power
+                    @test result_lightcurves.freq ≈ result_iterables.freq
+                    
+                    if return_auxil
+                        @test result_lightcurves.pds1 ≈ result_iterables.pds1
+                        @test result_lightcurves.pds2 ≈ result_iterables.pds2
+                    end
+                end
+            end
+        end
+    end
+end
+
+# Test EventList with no GTI (should create default GTI)
+let
+    data = setup_fourier_test_data()
+    eventlist_no_gti = EventList(data.times, nothing, FITSMetadata(
+        "[test]", 1, nothing, Dict{String,Vector}(), 
+        Dict{String,Any}(), nothing, nothing
+    ))
+    
+    expected_gti = Float64[minimum(data.times) maximum(data.times)]
+    
+    result_eventlist = avg_pds_from_eventlist(eventlist_no_gti, FOURIER_SEGMENT_SIZE, data.dt, silent=true)
+    result_times = avg_pds_from_events(data.times, expected_gti, FOURIER_SEGMENT_SIZE, data.dt, silent=true)
+    
+    @test result_eventlist.power ≈ result_times.power
+    @test result_eventlist.freq ≈ result_times.freq
+end
+
+# Test LightCurve with vector dt (should use first element)
+let
+    dt_vec = [1.0, 1.0, 1.0]
+    N = 100
+    times = collect(0:dt_vec[1]:(N-1)*dt_vec[1])
+    counts = rand(Poisson(100), N)
+    metadata = LightCurveMetadata(
+        "TEST", "INST", "TEST_OBJ", 58000.0, 
+        (times[1], times[end]), dt_vec[1],
+        [Dict{String,Any}()], Dict{String,Any}()
+    )
+    lc_vec_dt = LightCurve(times, dt_vec, counts, nothing, nothing, 
+                          EventProperty{Float64}[], metadata, :poisson)
+    flux_iterable = [counts]
+    result_lightcurve = avg_pds_from_lightcurve(lc_vec_dt, silent=true)
+    result_iterable = avg_pds_from_iterable(flux_iterable, dt_vec[1], silent=true)
+    @test result_lightcurve.power ≈ result_iterable.power
+    @test result_lightcurve.freq ≈ result_iterable.freq
+end
\ No newline at end of file
diff --git a/test/test_fourier/test_norm.jl b/test/test_fourier/test_norm.jl
new file mode 100644
index 0000000..7b440ce
--- /dev/null
+++ b/test/test_fourier/test_norm.jl
@@ -0,0 +1,130 @@
+const TEST_MEAN = TEST_VAR = 100000.0
+const TEST_N = 800000
+const TEST_DT = 0.2
+
+function setup_test_data(; N=TEST_N, dt=TEST_DT, mean=TEST_MEAN)
+    freq = fftfreq(N, dt)
+    good = freq .> 0
+    meanrate = mean / dt
+    lc = rand(Poisson(mean), N)
+    nph = sum(lc)
+    pds = keepat!(abs2.(fft(lc)), good)
+    return (
+        freq=freq, good=good, meanrate=meanrate, 
+        lc=lc, nph=nph, pds=pds, N=N, dt=dt, mean=mean
+    )
+end
+
+# Helper function for Poisson level tests
+function test_poisson_level(norm::String, power_type::String="all")
+    data = setup_test_data()
+    pdsnorm = normalize_periodograms(
+        data.pds, data.dt, data.N; 
+        mean_flux=data.mean, n_ph=data.nph, 
+        norm=norm, power_type=power_type
+    )
+    expected = poisson_level(norm; meanrate=data.meanrate, n_ph=data.nph)
+    @test Statistics.mean(pdsnorm) ≈ expected rtol=0.01
+end
+
+# test_positive_fft_bins
+let
+    freq = fftfreq(11)
+    goodbins = positive_fft_bins(11)
+    @test filter(x -> x > 0, freq) == freq[goodbins]
+end
+
+# test_norm_setup_and_basic_tests
+let
+    data = setup_test_data(N=1000000)
+    pdsabs = normalize_abs(data.pds, data.dt, size(data.lc, 1))
+    pdsfrac = normalize_frac(data.pds, data.dt, size(data.lc, 1), data.mean)
+    pois_abs = poisson_level("abs", meanrate=data.meanrate)
+    pois_frac = poisson_level("frac", meanrate=data.meanrate)
+
+    @test Statistics.mean(pdsabs) ≈ pois_abs rtol=0.02
+    @test Statistics.mean(pdsfrac) ≈ pois_frac rtol=0.02
+end
+
+# test_leahy_bksub_var_vs_standard
+let
+    data = setup_test_data()
+    lc_bksub = data.lc .- data.mean
+    pds_bksub = keepat!(abs2.(fft(lc_bksub)), data.good)
+    leahyvar = normalize_leahy_from_variance(pds_bksub, Statistics.var(lc_bksub), data.N)
+    leahy = 2 .* data.pds ./ sum(data.lc)
+    ratio = Statistics.mean(leahyvar ./ leahy)
+    @test ratio ≈ 1 rtol=0.01
+end
+
+# test_abs_bksub
+let
+    data = setup_test_data()
+    lc_bksub = data.lc .- data.mean
+    pds_bksub = keepat!(abs2.(fft(lc_bksub)), data.good)
+    ratio = normalize_abs(pds_bksub, data.dt, data.N) ./ normalize_abs(data.pds, data.dt, data.N)
+    @test Statistics.mean(ratio) ≈ 1 rtol=0.01
+end
+# test_frac_renorm_constant
+let
+    data = setup_test_data()
+    lc_renorm = data.lc / data.mean
+    pds_renorm = keepat!(abs2.(fft(lc_renorm)), data.good)
+    ratio = normalize_frac(pds_renorm, data.dt, data.N, 1) ./ normalize_frac(data.pds, data.dt, data.N, data.mean)
+    @test Statistics.mean(ratio) ≈ 1 rtol=0.01
+end
+
+# test_frac_to_abs_ctratesq
+let
+    data = setup_test_data()
+    
+    ratio = (
+        normalize_frac(data.pds, data.dt, data.N, data.mean)
+        ./ normalize_abs(data.pds, data.dt, data.N)
+        .* data.meanrate .^ 2
+    )
+    @test Statistics.mean(ratio) ≈ 1 rtol=0.01
+end
+
+# test_total_variance
+let
+    data = setup_test_data()
+    
+    vdk_total_variance = sum((data.lc .- data.mean) .^ 2)
+    ratio = Statistics.mean(data.pds) / vdk_total_variance
+    @test Statistics.mean(ratio) ≈ 1 rtol=0.01
+end
+
+# All Poisson level tests - much more concise now!
+@testset "Poisson Level Tests" begin
+    test_poisson_level("abs")
+    test_poisson_level("frac") 
+    test_poisson_level("leahy")
+    test_poisson_level("none")
+    
+    # Real power type tests
+    test_poisson_level("abs", "real")
+    test_poisson_level("frac", "real")
+    test_poisson_level("leahy", "real") 
+    test_poisson_level("none", "real")
+    
+    # Absolute power type tests
+    test_poisson_level("abs", "abs")
+    test_poisson_level("frac", "abs")
+    test_poisson_level("leahy", "abs")
+    test_poisson_level("none", "abs")
+end
+
+# test_normalize_with_variance
+let
+    data = setup_test_data()
+    pdsnorm = normalize_periodograms(data.pds, data.dt, data.N; mean_flux=data.mean, variance=TEST_VAR, norm="leahy")
+    @test Statistics.mean(pdsnorm) ≈ 2 rtol=0.01
+end
+
+# test_normalize_none
+let
+    data = setup_test_data()
+    pdsnorm = normalize_periodograms(data.pds, data.dt, data.N; mean_flux=data.mean, n_ph=data.nph, norm="none")
+    @test Statistics.mean(pdsnorm) ≈ Statistics.mean(data.pds) rtol=0.01
+end
\ No newline at end of file