Modelling features in the power spectrum

While this package is designed to model broadband power spectral densities that are power-law-like, we can add narrow features using basis functions.

This feature is currently experimental and may not work properly. Use it at your own risk!

Defining the model and simulating form the model

we defined a power spectral model as sum of two QPOs and a bending power-law.

using Plots
using Pioran
using Random
using Distributions
cont = SingleBendingPowerLaw(0.5,1., 3.6)
q1 =  QPO(7.3,0.1,2)
q2 = QPO(9e-4,2,11)
𝓟 = cont + q1 + q2
f = 10 .^range(-3,1,1000)
p1 = Plots.plot(f,𝓟.(f),xlabel="Frequency (day^-1)",ylabel="Power Spectral Density",legend=true,framestyle = :box,xscale=:log10,yscale=:log10,lw=2,ylim=(1e-5,1e3),label="total")
p1 = Plots.plot!(f,q1.(f),label="QPO1",linestyle=:dash)
Plots.plot!(p1,f,q2.(f),label="QPO2",linestyle=:dot)
Plots.plot!(f,cont.(f),label="cont",linestyle=:dash)
p1

We now simulate time series from this process. To do so we will use the GP method by approximating the bending power-law with basis functions.

t = range(0.,1e3,step=0.1)
yerr = 0.3*ones(length(t));

10001-element Vector{Float64}:
 0.3
 0.3
 0.3
 0.3
 0.3
 0.3
 0.3
 0.3
 0.3
 0.3
 ⋮
 0.3
 0.3
 0.3
 0.3
 0.3
 0.3
 0.3
 0.3
 0.3

f_min,f_max = 1e-3,1e2
𝓡 = approx(𝓟, f_min, f_max, 35, 1., basis_function="SHO")
GP = ScalableGP(0.0, 𝓡)
GPc = GP(t,yerr.^2)

AbstractGPs.FiniteGP{ScalableGP{AbstractGPs.GP{AbstractGPs.ConstMean{Float64}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}(
f: ScalableGP{AbstractGPs.GP{AbstractGPs.ConstMean{Float64}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}(AbstractGPs.GP{AbstractGPs.ConstMean{Float64}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}(AbstractGPs.ConstMean{Float64}(0.0), Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}(Celerite[Celerite(0.0022640560425049115, 0.0022640560425049115, 0.0002221441469079183, 0.0002221441469079183), Celerite(0.0018258374381283412, 0.0018258374381283412, 0.00037172765722227083, 0.00037172765722227083), Celerite(0.002622586160855605, 0.002622586160855605, 0.0006220350752758573, 0.0006220350752758573), Celerite(0.0033513472441581253, 0.0033513472441581253, 0.0010408900907851527, 0.0010408900907851527), Celerite(0.004318114150934439, 0.004318114150934439, 0.0017417863142431935, 0.0017417863142431935), Celerite(0.005632943777700911, 0.005632943777700911, 0.0029146396832315435, 0.0029146396832315435), Celerite(0.007203952150881374, 0.007203952150881374, 0.004877248382077916, 0.004877248382077916), Celerite(0.009458218124187118, 0.009458218124187118, 0.008161403935222519, 0.008161403935222519), Celerite(0.012001743513903906, 0.012001743513903906, 0.013656986270911955, 0.013656986270911955), Celerite(0.01591565679409391, 0.01591565679409391, 0.0228530869791819, 0.0228530869791819)  …  Celerite(4.382326422024152e-6, 4.382326422024152e-6, 241.86045634856222, 241.86045634856222), Celerite(1.6423844695175336e-6, 1.6423844695175336e-6, 404.7201876105593, 404.7201876105593), Celerite(3.253764625721934e-7, 3.253764625721934e-7, 677.2435342777359, 677.2435342777359), Celerite(1.0782099277320264e-7, 1.0782099277320264e-7, 1133.2738488507068, 1133.2738488507068), Celerite(2.3420179514461132e-8, 2.3420179514461132e-8, 1896.3778190345968, 1896.3778190345968), Celerite(7.261317849496646e-9, 7.261317849496646e-9, 3173.327290816335, 3173.327290816335), Celerite(1.5266208011970595e-9, 1.5266208011970595e-9, 5310.126491442592, 5310.126491442592), Celerite(6.789431574780224e-10, 6.789431574780224e-10, 8885.765876316733, 8885.765876316733), Celerite(0.8108997419288619, 0.20937341306229162, 0.15707963267948966, 0.6083668013960418), Celerite(0.01099713348643251, 0.0005003868998290508, 0.5711986642890533, 12.553382114129409)], [0.0022640560425049115, 0.0018258374381283412, 0.002622586160855605, 0.0033513472441581253, 0.004318114150934439, 0.005632943777700911, 0.007203952150881374, 0.009458218124187118, 0.012001743513903906, 0.01591565679409391  …  4.382326422024152e-6, 1.6423844695175336e-6, 3.253764625721934e-7, 1.0782099277320264e-7, 2.3420179514461132e-8, 7.261317849496646e-9, 1.5266208011970595e-9, 6.789431574780224e-10, 0.8108997419288619, 0.01099713348643251], [0.0022640560425049115, 0.0018258374381283412, 0.002622586160855605, 0.0033513472441581253, 0.004318114150934439, 0.005632943777700911, 0.007203952150881374, 0.009458218124187118, 0.012001743513903906, 0.01591565679409391  …  4.382326422024152e-6, 1.6423844695175336e-6, 3.253764625721934e-7, 1.0782099277320264e-7, 2.3420179514461132e-8, 7.261317849496646e-9, 1.5266208011970595e-9, 6.789431574780224e-10, 0.20937341306229162, 0.0005003868998290508], [0.0002221441469079183, 0.00037172765722227083, 0.0006220350752758573, 0.0010408900907851527, 0.0017417863142431935, 0.0029146396832315435, 0.004877248382077916, 0.008161403935222519, 0.013656986270911955, 0.0228530869791819  …  241.86045634856222, 404.7201876105593, 677.2435342777359, 1133.2738488507068, 1896.3778190345968, 3173.327290816335, 5310.126491442592, 8885.765876316733, 0.15707963267948966, 0.5711986642890533], [0.0002221441469079183, 0.00037172765722227083, 0.0006220350752758573, 0.0010408900907851527, 0.0017417863142431935, 0.0029146396832315435, 0.004877248382077916, 0.008161403935222519, 0.013656986270911955, 0.0228530869791819  …  241.86045634856222, 404.7201876105593, 677.2435342777359, 1133.2738488507068, 1896.3778190345968, 3173.327290816335, 5310.126491442592, 8885.765876316733, 0.6083668013960418, 12.553382114129409])), Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}(Celerite[Celerite(0.0022640560425049115, 0.0022640560425049115, 0.0002221441469079183, 0.0002221441469079183), Celerite(0.0018258374381283412, 0.0018258374381283412, 0.00037172765722227083, 0.00037172765722227083), Celerite(0.002622586160855605, 0.002622586160855605, 0.0006220350752758573, 0.0006220350752758573), Celerite(0.0033513472441581253, 0.0033513472441581253, 0.0010408900907851527, 0.0010408900907851527), Celerite(0.004318114150934439, 0.004318114150934439, 0.0017417863142431935, 0.0017417863142431935), Celerite(0.005632943777700911, 0.005632943777700911, 0.0029146396832315435, 0.0029146396832315435), Celerite(0.007203952150881374, 0.007203952150881374, 0.004877248382077916, 0.004877248382077916), Celerite(0.009458218124187118, 0.009458218124187118, 0.008161403935222519, 0.008161403935222519), Celerite(0.012001743513903906, 0.012001743513903906, 0.013656986270911955, 0.013656986270911955), Celerite(0.01591565679409391, 0.01591565679409391, 0.0228530869791819, 0.0228530869791819)  …  Celerite(4.382326422024152e-6, 4.382326422024152e-6, 241.86045634856222, 241.86045634856222), Celerite(1.6423844695175336e-6, 1.6423844695175336e-6, 404.7201876105593, 404.7201876105593), Celerite(3.253764625721934e-7, 3.253764625721934e-7, 677.2435342777359, 677.2435342777359), Celerite(1.0782099277320264e-7, 1.0782099277320264e-7, 1133.2738488507068, 1133.2738488507068), Celerite(2.3420179514461132e-8, 2.3420179514461132e-8, 1896.3778190345968, 1896.3778190345968), Celerite(7.261317849496646e-9, 7.261317849496646e-9, 3173.327290816335, 3173.327290816335), Celerite(1.5266208011970595e-9, 1.5266208011970595e-9, 5310.126491442592, 5310.126491442592), Celerite(6.789431574780224e-10, 6.789431574780224e-10, 8885.765876316733, 8885.765876316733), Celerite(0.8108997419288619, 0.20937341306229162, 0.15707963267948966, 0.6083668013960418), Celerite(0.01099713348643251, 0.0005003868998290508, 0.5711986642890533, 12.553382114129409)], [0.0022640560425049115, 0.0018258374381283412, 0.002622586160855605, 0.0033513472441581253, 0.004318114150934439, 0.005632943777700911, 0.007203952150881374, 0.009458218124187118, 0.012001743513903906, 0.01591565679409391  …  4.382326422024152e-6, 1.6423844695175336e-6, 3.253764625721934e-7, 1.0782099277320264e-7, 2.3420179514461132e-8, 7.261317849496646e-9, 1.5266208011970595e-9, 6.789431574780224e-10, 0.8108997419288619, 0.01099713348643251], [0.0022640560425049115, 0.0018258374381283412, 0.002622586160855605, 0.0033513472441581253, 0.004318114150934439, 0.005632943777700911, 0.007203952150881374, 0.009458218124187118, 0.012001743513903906, 0.01591565679409391  …  4.382326422024152e-6, 1.6423844695175336e-6, 3.253764625721934e-7, 1.0782099277320264e-7, 2.3420179514461132e-8, 7.261317849496646e-9, 1.5266208011970595e-9, 6.789431574780224e-10, 0.20937341306229162, 0.0005003868998290508], [0.0002221441469079183, 0.00037172765722227083, 0.0006220350752758573, 0.0010408900907851527, 0.0017417863142431935, 0.0029146396832315435, 0.004877248382077916, 0.008161403935222519, 0.013656986270911955, 0.0228530869791819  …  241.86045634856222, 404.7201876105593, 677.2435342777359, 1133.2738488507068, 1896.3778190345968, 3173.327290816335, 5310.126491442592, 8885.765876316733, 0.15707963267948966, 0.5711986642890533], [0.0002221441469079183, 0.00037172765722227083, 0.0006220350752758573, 0.0010408900907851527, 0.0017417863142431935, 0.0029146396832315435, 0.004877248382077916, 0.008161403935222519, 0.013656986270911955, 0.0228530869791819  …  241.86045634856222, 404.7201876105593, 677.2435342777359, 1133.2738488507068, 1896.3778190345968, 3173.327290816335, 5310.126491442592, 8885.765876316733, 0.6083668013960418, 12.553382114129409]))
x: 0.0:0.1:1000.0
Σy: [0.09 0.0 … 0.0 0.0; 0.0 0.09 … 0.0 0.0; … ; 0.0 0.0 … 0.09 0.0; 0.0 0.0 … 0.0 0.09]
)

We can sample realisations from this process:

rng = MersenneTwister(12)
y = [rand(rng,GPc) for i in 1:100]
Plots.plot(t,y[1:3],xlabel="Time",ylabel="Value",legend=false,framestyle = :box,ms=3)

As it is not very informative to look at the time series, let's compute the periodogram using Tonari.jl.

using Tonari
x_GP = mapreduce(permutedims,vcat,y)
fP,I = periodogram(t,x_GP',apply_end_matching=false)

([0.000999900009999, 0.001999800019998, 0.0029997000299970002, 0.003999600039996, 0.004999500049995, 0.0059994000599940004, 0.006999300069993001, 0.007999200079992, 0.008999100089991, 0.00999900009999  …  4.99050094990501, 4.991500849915008, 4.992500749925007, 4.993500649935006, 4.9945005499450055, 4.995500449955005, 4.996500349965004, 4.997500249975002, 4.998500149985001, 4.999500049995], [33.33070329338245, 21.62041921319835, 17.479122056968748, 14.863564559267934, 16.40865407592501, 13.591805393994216, 13.514361448514356, 10.556124876804024, 12.751340993052017, 11.209468492725884  …  0.022866277287100126, 0.021220394763766033, 0.024548521191515096, 0.018368991425055662, 0.024147200946210267, 0.021012794373673222, 0.022959600604151215, 0.02098017770866714, 0.02111036102466429, 0.023293889178581426])

Δt = t[2]-t[1]
noise_level = 2Δt*mean(yerr.^2)

Plots.plot(fP,I,yscale=:log10,xscale=:log10,xlabel="Frequency (Hz)",ylabel="Power",label="Periodogram",framestyle=:box)
Plots.plot!(f,𝓟.(f),label="Model",linewidth=2)
hline!([noise_level],label="Noise level",linewidth=2,linestyle=:dash,ylim=(noise_level/10,1e3),)

We can even compare the periodogram to the one of simulations using Timmer and König. See the documentation of Tonari.jl for more details:

T, Δt = 504.2, 0.132
simu = Simulation(𝓟, 1e3, 0.1, 20, 20)
rng = MersenneTwister(42)
N = 100
t_TK, x_TK, yerr_TK = sample(rng, simu, N, error_size = 0.25)
fP_TK,I_TK = periodogram(t_TK,x_TK,apply_end_matching=false)

([0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009000000000000001, 0.01  …  4.9910000000000005, 4.992, 4.993, 4.994, 4.995, 4.996, 4.997, 4.998, 4.999, 5.0], [50.60840630598973, 20.281774070297036, 30.834240939618464, 25.807963999220195, 15.697843014051099, 19.179527178870774, 16.5902728983377, 21.37396707249356, 29.145020347851812, 20.11771612020057  …  0.023900335603256564, 0.022882969906295212, 0.024639147281892895, 0.027820229883572162, 0.02987817536149938, 0.030568122632746844, 0.028762870394919517, 0.02454477593557762, 0.03355210991990242, 0.030121193090363232])

At the moment there is a discrepancy in the normalisations but it will be solved eventually...

Δt = t[2]-t[1]
noise_level = 2Δt*mean(yerr.^2)
noise_level_TK = 2Δt*mean(yerr_TK.^2)

Plots.plot(fP,I,yscale=:log10,xscale=:log10,xlabel="Frequency (Hz)",ylabel="Power",framestyle=:box,label="Periodogram GP",alpha=1,lw=.5)
Plots.plot!(fP_TK,I_TK,label="Periodogram TK",alpha=0.2)
Plots.plot!(f,𝓟.(f)/2,label="Model",linewidth=2)
hline!([noise_level],label="Noise level",linewidth=2,linestyle=:dash,ylim=(noise_level/4,2e2),)